Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-09 09:09:45 

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4QSS3
0027 //
0028 // G4QSS3 simulator
0029 
0030 // Authors: Lucio Santi, Rodrigo Castro (Univ. Buenos Aires), 2018-2021
0031 // --------------------------------------------------------------------
0032 #ifndef G4QSS3_HH
0033 #define G4QSS3_HH
0034 
0035 #include "G4Types.hh"
0036 #include "G4qss_misc.hh"
0037 
0038 #include <cmath>
0039 
0040 /**
0041  * @brief G4QSS3 defines the QSS3 simulator engine used in QSS field stepper.
0042  */
0043 
0044 class G4QSS3
0045 {
0046   public:
0047 
0048     G4QSS3(QSS_simulator);
0049 
0050     inline QSS_simulator getSimulator() const { return this->simulator; }
0051 
0052     inline G4int order() const { return 3; }
0053 
0054     inline void full_definition(G4double coeff)
0055     {
0056       G4double* const x = simulator->q;
0057       G4double* const dx = simulator->x;
0058       G4double* const alg = simulator->alg;
0059 
0060       dx[1] = x[12];
0061       dx[2] = 0;
0062       dx[3] = 0;
0063 
0064       dx[5] = x[16];
0065       dx[6] = 0;
0066       dx[7] = 0;
0067 
0068       dx[9] = x[20];
0069       dx[10] = 0;
0070       dx[11] = 0;
0071 
0072       dx[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
0073       dx[14] = 0;
0074       dx[15] = 0;
0075 
0076       dx[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
0077       dx[18] = 0;
0078       dx[19] = 0;
0079 
0080       dx[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
0081       dx[22] = 0;
0082       dx[23] = 0;
0083     }
0084 
0085     inline void dependencies(G4int i, G4double coeff)
0086     {
0087       G4double* const x = simulator->q;
0088       G4double* const der = simulator->x;
0089       G4double* const alg = simulator->alg;
0090 
0091       switch (i)
0092       {
0093         case 0:
0094           der[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
0095           der[14] = ((alg[2] * x[17] - x[21] * alg[1]) * coeff) / 2;
0096           der[15] = (coeff * (alg[2] * x[18] - x[22] * alg[1])) / 3;
0097 
0098           der[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
0099           der[18] = ((alg[0] * x[21] - alg[2] * x[13]) * coeff) / 2;
0100           der[19] = (coeff * (alg[0] * x[22] - alg[2] * x[14])) / 3;
0101 
0102           der[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
0103           der[22] = (coeff * (x[13] * alg[1] - alg[0] * x[17])) / 2;
0104           der[23] = (coeff * (alg[1] * x[14] - x[18] * alg[0])) / 3;
0105           return;
0106         case 1:
0107           der[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
0108           der[14] = ((alg[2] * x[17] - x[21] * alg[1]) * coeff) / 2;
0109           der[15] = (coeff * (alg[2] * x[18] - x[22] * alg[1])) / 3;
0110 
0111           der[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
0112           der[18] = ((alg[0] * x[21] - alg[2] * x[13]) * coeff) / 2;
0113           der[19] = (coeff * (alg[0] * x[22] - alg[2] * x[14])) / 3;
0114 
0115           der[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
0116           der[22] = (coeff * (x[13] * alg[1] - alg[0] * x[17])) / 2;
0117           der[23] = (coeff * (alg[1] * x[14] - x[18] * alg[0])) / 3;
0118           return;
0119         case 2:
0120           der[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
0121           der[14] = ((alg[2] * x[17] - x[21] * alg[1]) * coeff) / 2;
0122           der[15] = (coeff * (alg[2] * x[18] - x[22] * alg[1])) / 3;
0123 
0124           der[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
0125           der[18] = ((alg[0] * x[21] - alg[2] * x[13]) * coeff) / 2;
0126           der[19] = (coeff * (alg[0] * x[22] - alg[2] * x[14])) / 3;
0127 
0128           der[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
0129           der[22] = (coeff * (x[13] * alg[1] - alg[0] * x[17])) / 2;
0130           der[23] = (coeff * (alg[1] * x[14] - x[18] * alg[0])) / 3;
0131           return;
0132         case 3:
0133           der[1] = x[12];
0134           der[2] = x[13] / 2;
0135           der[3] = x[14] / 3;
0136 
0137           der[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
0138           der[18] = ((alg[0] * x[21] - alg[2] * x[13]) * coeff) / 2;
0139           der[19] = (coeff * (alg[0] * x[22] - alg[2] * x[14])) / 3;
0140 
0141           der[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
0142           der[22] = (coeff * (x[13] * alg[1] - alg[0] * x[17])) / 2;
0143           der[23] = (coeff * (alg[1] * x[14] - x[18] * alg[0])) / 3;
0144           return;
0145         case 4:
0146           der[5] = x[16];
0147           der[6] = x[17] / 2;
0148           der[7] = x[18] / 3;
0149 
0150           der[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
0151           der[14] = ((alg[2] * x[17] - x[21] * alg[1]) * coeff) / 2;
0152           der[15] = (coeff * (alg[2] * x[18] - x[22] * alg[1])) / 3;
0153 
0154           der[21] = coeff * (alg[1] * x[12] - alg[0] * x[16]);
0155           der[22] = (coeff * (x[13] * alg[1] - alg[0] * x[17])) / 2;
0156           der[23] = (coeff * (alg[1] * x[14] - x[18] * alg[0])) / 3;
0157           return;
0158         case 5:
0159           der[9] = x[20];
0160           der[10] = x[21] / 2;
0161           der[11] = x[22] / 3;
0162 
0163           der[13] = coeff * (alg[2] * x[16] - alg[1] * x[20]);
0164           der[14] = ((alg[2] * x[17] - x[21] * alg[1]) * coeff) / 2;
0165           der[15] = (coeff * (alg[2] * x[18] - x[22] * alg[1])) / 3;
0166 
0167           der[17] = coeff * (alg[0] * x[20] - alg[2] * x[12]);
0168           der[18] = ((alg[0] * x[21] - alg[2] * x[13]) * coeff) / 2;
0169           der[19] = (coeff * (alg[0] * x[22] - alg[2] * x[14])) / 3;
0170           return;
0171       }
0172     }
0173 
0174     void recompute_next_times(G4int* inf, G4double t);  // __attribute__((hot));
0175 
0176     inline void recompute_all_state_times(G4double t)
0177     {
0178       G4double mpr;
0179       G4double* const x = simulator->x;
0180       G4double* const lqu = simulator->lqu;
0181       G4double* const time = simulator->nextStateTime;
0182 
0183       for (G4int var = 0, icf0 = 0; var < 6; var++, icf0 += 4)
0184       {
0185         const G4int icf1 = icf0 + 1;
0186 
0187         if (x[icf1] == 0)
0188         {
0189           time[var] = Qss_misc::INF;
0190         }
0191         else
0192         {
0193           mpr = lqu[var] / x[icf1];
0194           if (mpr < 0) { mpr *= -1; }
0195           time[var] = t + mpr;
0196         }
0197       }
0198     }
0199 
0200     inline void next_time(G4int i, G4double t)
0201     {
0202       const G4int cf3 = 4 * i + 3;
0203       G4double* const x = simulator->x;
0204       G4double* const lqu = simulator->lqu;
0205       G4double* const time = simulator->nextStateTime;
0206 
0207       if (likely(x[cf3])) {
0208         time[i] = t + std::cbrt(lqu[i] / std::fabs(x[cf3]));
0209       } else {
0210         time[i] = Qss_misc::INF;
0211       }
0212     }
0213 
0214     inline void update_quantized_state(G4int i)
0215     {
0216       const G4int cf0 = i * 4, cf1 = cf0 + 1, cf2 = cf1 + 1;
0217       G4double* const q = simulator->q;
0218       G4double* const x = simulator->x;
0219 
0220       q[cf0] = x[cf0];
0221       q[cf1] = x[cf1];
0222       q[cf2] = x[cf2];
0223     }
0224 
0225     inline void reset_state(G4int i, G4double value)
0226     {
0227       G4double* const x = simulator->x;
0228       G4double* const q = simulator->q;
0229       G4double* const tq = simulator->tq;
0230       G4double* const tx = simulator->tx;
0231       const G4int idx = 4 * i;
0232 
0233       x[idx] = value;
0234 
0235       simulator->lqu[i] = simulator->dQRel[i] * std::fabs(value);
0236       if (simulator->lqu[i] < simulator->dQMin[i])
0237       {
0238         simulator->lqu[i] = simulator->dQMin[i];
0239       }
0240       q[idx] = value;
0241       q[idx + 1] = q[idx + 2] = tq[i] = tx[i] = 0;
0242     }
0243 
0244     inline G4double evaluate_x_poly(G4int i, G4double dt, G4double* p)
0245     {
0246       return ((p[i + 3] * dt + p[i + 2]) * dt + p[i + 1]) * dt + p[i];
0247     }
0248 
0249     inline void advance_time_q(G4int i, G4double dt)  //  __attribute__((hot))
0250     {
0251       G4double* const p = simulator->q;
0252       const G4int i0 = i, i1 = i0 + 1, i2 = i1 + 1;
0253       p[i0] = (p[i2] * dt + p[i1]) * dt + p[i0];
0254       p[i1] = p[i1] + 2 * dt * p[i2];
0255     }
0256 
0257     inline void advance_time_x(G4int i, G4double dt)  // __attribute__((hot))
0258     {
0259       G4double* const p = simulator->x;
0260       const G4int i0 = i, i1 = i0 + 1, i2 = i1 + 1, i3 = i2 + 1;
0261       p[i0] = ((p[i3] * dt + p[i2]) * dt + p[i1]) * dt + p[i0];
0262       p[i1] = (3 * p[i3] * dt + 2 * p[i2]) * dt + p[i1];
0263       p[i2] = p[i2] + 3 * dt * p[i3];
0264     }
0265 
0266     G4double min_pos_root(G4double* coeff, G4int order);
0267 
0268     inline G4double min_pos_root_2(G4double* coeff)
0269     {
0270       G4double mpr = Qss_misc::INF;
0271 
0272       if (coeff[2] == 0 || (1000 * std::fabs(coeff[2])) < std::fabs(coeff[1]))
0273       {
0274         if (coeff[1] == 0) {
0275           mpr = Qss_misc::INF;
0276         } else {
0277           mpr = -coeff[0] / coeff[1];
0278         }
0279 
0280         if (mpr < 0) { mpr = Qss_misc::INF; }
0281       }
0282       else
0283       {
0284         G4double disc;
0285         disc = coeff[1] * coeff[1] - 4 * coeff[2] * coeff[0];
0286         if (disc < 0)   // no real roots
0287         {
0288           mpr = Qss_misc::INF;
0289         }
0290         else
0291         {
0292           G4double sd, r1;
0293           G4double cf2_d2 = 2 * coeff[2];
0294 
0295           sd = std::sqrt(disc);
0296           r1 = (-coeff[1] + sd) / cf2_d2;
0297           if (r1 > 0) {
0298             mpr = r1;
0299           } else {
0300             mpr = Qss_misc::INF; 
0301           }
0302           r1 = (-coeff[1] - sd) / cf2_d2;
0303           if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
0304         }
0305       }
0306 
0307       return mpr;
0308     }  // __attribute__((hot))
0309 
0310     inline G4double min_pos_root_3(G4double* coeff)
0311     {
0312       G4double mpr = Qss_misc::INF;
0313       static const G4double sqrt3 = std::sqrt(3);
0314 
0315       if ((coeff[3] == 0) || (1000 * std::fabs(coeff[3]) < std::fabs(coeff[2])))
0316       {
0317         mpr = min_pos_root_2(coeff);
0318       }
0319       else if (coeff[0] == 0)
0320       {
0321         if (coeff[1] == 0)
0322         {
0323           mpr = -coeff[2] / coeff[3];
0324         }
0325         else
0326         {
0327           coeff[0] = coeff[1];
0328           coeff[1] = coeff[2];
0329           coeff[2] = coeff[3];
0330           mpr = min_pos_root_2(coeff);
0331         }
0332       }
0333       else
0334       {
0335         G4double q, r, disc, q3;
0336         G4double val = coeff[2] / 3 / coeff[3];
0337         G4double cf32 = coeff[3] * coeff[3];
0338         G4double cf22 = coeff[2] * coeff[2];
0339         G4double denq = 9 * cf32;
0340         G4double denr = 6 * coeff[3] * denq;
0341         G4double rcomm = 9 * coeff[3] * coeff[2] * coeff[1] - 2 * cf22 * coeff[2];
0342 
0343         q = (3 * coeff[3] * coeff[1] - cf22) / denq;
0344         q3 = q * q * q;
0345 
0346         r = (rcomm - 27 * cf32 * coeff[0]) / denr;
0347         disc = q3 + r * r;
0348         mpr = Qss_misc::INF;
0349 
0350         if (disc >= 0)
0351         {
0352           G4double sd, sx, t, r1, rsd;
0353           sd = std::sqrt(disc);
0354           rsd = r + sd;
0355           if (rsd > 0) {
0356             sx = std::cbrt(rsd);
0357           } else {
0358             sx = -std::cbrt(std::fabs(rsd));
0359           }
0360 
0361           rsd = r - sd;
0362           if (rsd > 0) {
0363             t = std::cbrt(rsd);
0364           } else {
0365             t = -std::cbrt(std::fabs(rsd));
0366           }
0367 
0368           r1 = sx + t - val;
0369 
0370           if (r1 > 0) { mpr = r1; }
0371         }
0372         else
0373         {
0374           // three real roots
0375           G4double rho, th, rho13, costh3, sinth3, spt, smti32, r1;
0376           rho = std::sqrt(-q3);
0377           th = std::acos(r / rho);
0378           rho13 = std::cbrt(rho);
0379           costh3 = std::cos(th / 3);
0380           sinth3 = std::sin(th / 3);
0381           spt = rho13 * 2 * costh3;
0382           smti32 = -rho13 * sinth3 * sqrt3;
0383           r1 = spt - val;
0384           if (r1 > 0) { mpr = r1; }
0385           r1 = -spt / 2 - val + smti32;
0386           if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
0387           r1 = r1 - 2 * smti32;
0388           if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
0389         }
0390       }
0391 
0392       return mpr;
0393     }  // __attribute__((hot))
0394 
0395     inline G4double min_pos_root_2_alt(G4double* coeff, G4double cf0Alt)
0396     {
0397       G4double mpr = Qss_misc::INF;
0398       G4double mpr2;
0399 
0400       if (coeff[2] == 0 || (1000 * std::fabs(coeff[2])) < std::fabs(coeff[1]))
0401       {
0402         if (coeff[1] == 0)
0403         {
0404           mpr = Qss_misc::INF;
0405         }
0406         else
0407         {
0408           mpr = -coeff[0] / coeff[1];
0409           mpr2 = -cf0Alt / coeff[1];
0410           if (mpr < 0 || (mpr2 > 0 && mpr2 < mpr)) { mpr = mpr2; }
0411         }
0412 
0413         if (mpr < 0) { mpr = Qss_misc::INF; }
0414       }
0415       else
0416       {
0417         G4double cf1_2 = coeff[1] * coeff[1];
0418         G4double cf2_4 = 4 * coeff[2];
0419         G4double disc1 = cf1_2 - cf2_4 * coeff[0];
0420         G4double disc2 = cf1_2 - cf2_4 * cf0Alt;
0421         G4double cf2_d2 = 2 * coeff[2];
0422 
0423         if (unlikely(disc1 < 0 && disc2 < 0))
0424         {
0425           mpr = Qss_misc::INF;
0426         }
0427         else if (disc2 < 0)
0428         {
0429           G4double sd, r1;
0430           sd = std::sqrt(disc1);
0431           r1 = (-coeff[1] + sd) / cf2_d2;
0432           if (r1 > 0) {
0433             mpr = r1;
0434           } else {
0435             mpr = Qss_misc::INF;
0436           }
0437           r1 = (-coeff[1] - sd) / cf2_d2;
0438           if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
0439         }
0440         else if (disc1 < 0)
0441         {
0442           G4double sd, r1;
0443           sd = std::sqrt(disc2);
0444           r1 = (-coeff[1] + sd) / cf2_d2;
0445           if (r1 > 0) {
0446             mpr = r1;
0447           } else {
0448             mpr = Qss_misc::INF;
0449           }
0450           r1 = (-coeff[1] - sd) / cf2_d2;
0451           if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
0452         }
0453         else
0454         {
0455           G4double sd1, r1, sd2, r2;
0456           sd1 = std::sqrt(disc1);
0457           sd2 = std::sqrt(disc2);
0458           r1 = (-coeff[1] + sd1) / cf2_d2;
0459           r2 = (-coeff[1] + sd2) / cf2_d2;
0460 
0461           if (r1 > 0) {
0462             mpr = r1;
0463           } else {
0464             mpr = Qss_misc::INF;
0465           }
0466           r1 = (-coeff[1] - sd1) / cf2_d2;
0467           if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
0468 
0469           if (r2 > 0 && r2 < mpr) { mpr = r2; }
0470           r2 = (-coeff[1] - sd2) / cf2_d2;
0471           if ((r2 > 0) && (r2 < mpr)) { mpr = r2; }
0472         }
0473       }
0474 
0475       return mpr;
0476     }  // __attribute__((hot))
0477 
0478     inline G4double min_pos_root_3_alt(G4double* coeff, G4double cf0Alt)
0479     {
0480       G4double mpr = Qss_misc::INF;
0481       static const G4double sqrt3 = std::sqrt(3);
0482 
0483       if ((coeff[3] == 0) || (1000 * std::fabs(coeff[3]) < std::fabs(coeff[2])))
0484       {
0485         mpr = min_pos_root_2_alt(coeff, cf0Alt);
0486       }
0487       else if (coeff[0] == 0)
0488       {
0489         G4double mpr2;
0490         coeff[0] = cf0Alt;
0491         mpr = min_pos_root_3(coeff);
0492 
0493         if (coeff[1] == 0)
0494         {
0495           mpr2 = -coeff[2] / coeff[3];
0496         }
0497         else
0498         {
0499           coeff[0] = coeff[1];
0500           coeff[1] = coeff[2];
0501           coeff[2] = coeff[3];
0502           mpr2 = min_pos_root_2(coeff);
0503         }
0504 
0505         if (mpr2 > 0 && mpr2 < mpr) { mpr = mpr2; }
0506       }
0507       else if (cf0Alt == 0)
0508       {
0509         G4double mpr2;
0510         mpr = min_pos_root_3(coeff);
0511 
0512         if (coeff[1] == 0)
0513         {
0514           mpr2 = -coeff[2] / coeff[3];
0515         }
0516         else
0517         {
0518           coeff[0] = coeff[1];
0519           coeff[1] = coeff[2];
0520           coeff[2] = coeff[3];
0521           mpr2 = min_pos_root_2(coeff);
0522         }
0523 
0524         if (mpr2 > 0 && mpr2 < mpr) { mpr = mpr2; }
0525       }
0526       else
0527       {
0528         G4double q, r, rAlt, disc, discAlt, q3;
0529         G4double val = coeff[2] / 3 / coeff[3];
0530         G4double cf32 = coeff[3] * coeff[3];
0531         G4double cf22 = coeff[2] * coeff[2];
0532         G4double denq = 9 * cf32;
0533         G4double denr = 6 * coeff[3] * denq;
0534         G4double rcomm = 9 * coeff[3] * coeff[2] * coeff[1] - 2 * cf22 * coeff[2];
0535 
0536         q = (3 * coeff[3] * coeff[1] - cf22) / denq;
0537         q3 = q * q * q;
0538 
0539         r = (rcomm - 27 * cf32 * coeff[0]) / denr;
0540         rAlt = (rcomm - 27 * cf32 * cf0Alt) / denr;
0541 
0542         disc = q3 + r * r;
0543         discAlt = q3 + rAlt * rAlt;
0544         mpr = Qss_misc::INF;
0545 
0546         if (disc >= 0)
0547         {
0548           G4double sd, sx, t, r1, rsd;
0549           sd = std::sqrt(disc);
0550           rsd = r + sd;
0551           if (rsd > 0) {
0552             sx = std::cbrt(rsd);
0553           } else {
0554             sx = -std::cbrt(std::fabs(rsd));
0555           }
0556 
0557           rsd = r - sd;
0558           if (rsd > 0) {
0559             t = std::cbrt(rsd);
0560           } else {
0561             t = -std::cbrt(std::fabs(rsd));
0562           }
0563 
0564           r1 = sx + t - val;
0565 
0566           if (r1 > 0) { mpr = r1; }
0567 
0568           if (discAlt >= 0)
0569           {
0570             G4double sdAlt, sAlt, tAlt, r1Alt, rsdAlt;
0571             sdAlt = std::sqrt(discAlt);
0572             rsdAlt = rAlt + sdAlt;
0573             if (rsdAlt > 0) {
0574               sAlt = std::cbrt(rsdAlt);
0575             } else {
0576               sAlt = -std::cbrt(std::fabs(rsdAlt));
0577             }
0578 
0579             rsdAlt = rAlt - sdAlt;
0580             if (rsdAlt > 0) {
0581               tAlt = std::cbrt(rsdAlt);
0582             } else {
0583               tAlt = -std::cbrt(std::fabs(rsdAlt));
0584             }
0585 
0586             r1Alt = sAlt + tAlt - val;
0587 
0588             if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
0589           }
0590           else
0591           {
0592             G4double rho, th, rho13, costh3, sinth3, spt, smti32, r1Alt;
0593 
0594             rho = std::sqrt(-q3);
0595             th = std::acos(rAlt / rho);
0596             rho13 = std::cbrt(rho);
0597             costh3 = std::cos(th / 3);
0598             sinth3 = std::sin(th / 3);
0599             spt = rho13 * 2 * costh3;
0600             smti32 = -rho13 * sinth3 * sqrt3;
0601             r1Alt = spt - val;
0602             if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
0603             r1Alt = -spt / 2 - val + smti32;
0604             if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
0605             r1Alt = r1Alt - 2 * smti32;
0606             if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
0607           }
0608         }
0609         else
0610         {
0611           G4double rho, th, rho13, costh3, sinth3, spt, smti32, r1;
0612 
0613           rho = std::sqrt(-q3);
0614           th = std::acos(r / rho);
0615           rho13 = std::cbrt(rho);
0616           costh3 = std::cos(th / 3);
0617           sinth3 = std::sin(th / 3);
0618           spt = rho13 * 2 * costh3;
0619           smti32 = -rho13 * sinth3 * sqrt3;
0620           r1 = spt - val;
0621           if (r1 > 0) { mpr = r1; }
0622           r1 = -spt / 2 - val + smti32;
0623           if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
0624           r1 = r1 - 2 * smti32;
0625           if ((r1 > 0) && (r1 < mpr)) { mpr = r1; }
0626 
0627           if (discAlt >= 0)
0628           {
0629             G4double sdAlt, sAlt, tAlt, r1Alt, rsdAlt;
0630             sdAlt = std::sqrt(discAlt);
0631             rsdAlt = rAlt + sdAlt;
0632             if (rsdAlt > 0) {
0633               sAlt = std::cbrt(rsdAlt);
0634             } else {
0635               sAlt = -std::cbrt(std::fabs(rsdAlt));
0636             }
0637 
0638             rsdAlt = rAlt - sdAlt;
0639             if (rsdAlt > 0) {
0640               tAlt = std::cbrt(rsdAlt);
0641             } else {
0642               tAlt = -std::cbrt(std::fabs(rsdAlt));
0643             }
0644 
0645             r1Alt = sAlt + tAlt - val;
0646 
0647             if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
0648           }
0649           else
0650           {
0651             G4double thAlt, costh3Alt, sinth3Alt, sptAlt, smti32Alt, r1Alt;
0652             thAlt = std::acos(rAlt / rho);
0653             costh3Alt = std::cos(thAlt / 3);
0654             sinth3Alt = std::sin(thAlt / 3);
0655             sptAlt = rho13 * 2 * costh3Alt;
0656             smti32Alt = -rho13 * sinth3Alt * sqrt3;
0657             r1Alt = sptAlt - val;
0658             if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
0659             r1Alt = -sptAlt / 2 - val + smti32Alt;
0660             if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
0661             r1Alt = r1Alt - 2 * smti32Alt;
0662             if (r1Alt > 0 && r1Alt < mpr) { mpr = r1Alt; }
0663           }
0664         }
0665       }
0666 
0667       return mpr;
0668     }
0669 
0670   private:
0671 
0672     QSS_simulator simulator;
0673 };
0674 
0675 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!