include/Geant4/G4QSS3.hh

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