src/jet/Martini.cc

0001 /*******************************************************************************

0002  * Copyright (c) The JETSCAPE Collaboration, 2018

0003  *

0004  * Modular, task-based framework for simulating all aspects of heavy-ion collisions

0005  *

0006  * For the list of contributors see AUTHORS.

0007  *

0008  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues

0009  *

0010  * or via email to bugs.jetscape@gmail.com

0011  *

0012  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).

0013  * See COPYING for details.

0014  ******************************************************************************/
0015
0016 #include "Martini.h"
0017 #include "JetScapeLogger.h"
0018 #include "JetScapeXML.h"
0019 #include <string>
0020
0021 #include "tinyxml2.h"
0022 #include <iostream>
0023
0024 #include "FluidDynamics.h"
0025 #include "MartiniMutex.h"
0026 #define hbarc 0.197327053
0027
0028 #define MAGENTA "\033[35m"
0029
0030 using namespace Jetscape;
0031
0032 // Register the module with the base class

0033 RegisterJetScapeModule<Martini> Martini::reg("Martini");
0034
0035 using std::ofstream;
0036 using std::ifstream;
0037 using std::ostream;
0038 using std::ios;
0039
0040 int Martini::pLabelNew = 0;
0041
0042 Martini::Martini() {
0043   SetId("Martini");
0044   VERBOSE(8);
0045
0046   //vectors for elastic rates:

0047   dGamma_qq = new vector<double>;
0048   dGamma_qg = new vector<double>;
0049   dGamma_qq_q = new vector<double>;
0050   dGamma_qg_q = new vector<double>;
0051
0052   // create and set Martini Mutex

0053   auto martini_mutex = make_shared<MartiniMutex>();
0054   SetMutex(martini_mutex);
0055 }
0056
0057 Martini::~Martini() { VERBOSE(8); }
0058
0059 void Martini::Init() {
0060   JSINFO << "Initialize Martini ...";
0061
0062   double deltaT = 0.0;
0063   double Martini_deltaT_Max = 0.03 + rounding_error;
0064
0065   deltaT = GetXMLElementDouble({"Eloss", "deltaT"});
0066
0067   if (deltaT > Martini_deltaT_Max) {
0068     JSWARN << "Timestep for Martini ( deltaT = " << deltaT
0069            << " ) is too large. "
0070            << "Please choose a detaT smaller than or equal to 0.03 in the XML "
0071               "file.";
0072     throw std::runtime_error(
0073         "Martini not properly initialized in XML file ...");
0074   }
0075
0076   string s = GetXMLElementText({"Eloss", "Martini", "name"});
0077   JSDEBUG << s << " to be initialized ...";
0078
0079   tStart = GetXMLElementDouble({"Eloss", "tStart"});
0080   Q0 = GetXMLElementDouble({"Eloss", "Martini", "Q0"});
0081   alpha_s0 = GetXMLElementDouble({"Eloss", "Martini", "alpha_s"});
0082   pcut = GetXMLElementDouble({"Eloss", "Martini", "pcut"});
0083   hydro_Tc = GetXMLElementDouble({"Eloss", "Martini", "hydro_Tc"});
0084   recoil_on = GetXMLElementInt({"Eloss", "Martini", "recoil_on"});
0085   run_alphas = GetXMLElementInt({"Eloss", "Martini", "run_alphas"});
0086
0087   alpha_em = 1. / 137.;
0088
0089   // Path to additional data

0090   PathToTables = GetXMLElementText({"Eloss", "Martini", "path"});
0091
0092   // Initialize random number distribution

0093   ZeroOneDistribution = uniform_real_distribution<double>{0.0, 1.0};
0094
0095   readRadiativeRate(&dat, &Gam);
0096   readElasticRateOmega();
0097   readElasticRateQ();
0098 }
0099
0100 void Martini::DoEnergyLoss(double deltaT, double Time, double Q2,
0101                            vector<Parton> &pIn, vector<Parton> &pOut) {
0102   VERBOSESHOWER(5) << MAGENTA << "SentInPartons Signal received : " << deltaT
0103                    << " " << Q2 << " " << pIn.size();
0104
0105   // particle info

0106   int Id, newId;
0107   int pStat;               // status of parton:

0108                            // 0: normal parton, 1: recoil parton,

0109                            // -1: sampled thermal parton (negative)

0110   int pLabel;              // particle number

0111   double p0, px, py, pz; // momentum of initial parton (pIn)

0112   double pAbs;
0113   double velx, vely, velz;
0114   double pRest, pxRest;    // momentum in the rest frame of fluid cell (pIn)

0115   double pyRest, pzRest;
0116   double k, kRest;    // momentum of radiated parton (pOut)

0117   double pNew, pxNew; // momentum of final parton (pOut)

0118   double pyNew, pzNew;
0119   double pNewRest;            // momentum in the rest frame of fluid cell (pOut)

0120   double omega, q;            // transferred energy/momentum for scattering

0121   double pThermal, pxThermal; // momentum of thermal parton (pOut)

0122   double pyThermal, pzThermal;
0123   double pRecoil, pxRecoil; // momentum of recoil parton (pOut)

0124   double pyRecoil, pzRecoil;
0125   double pRecoilRest;
0126   double xx, yy, zz, tt;    // position of initial parton (pIn)

0127   FourVector pVec, pVecNew; // 4-vectors of momenta before & after process

0128   FourVector kVec;          // 4-vector of momentum of radiated parton

0129   FourVector pVecRest;      // 4-vectors in the rest frame of fluid cell

0130   FourVector pVecNewRest;
0131   FourVector qVec; // 4-vector of momentum transfer (Note: space-like)

0132   FourVector
0133       pVecThermalRest;       // 4-vectors of thermal parton in the rest frame of

0134   FourVector pVecThermal;    // fluid cell and lab frame

0135   FourVector pVecRecoilRest; // 4-vectors of recoil parton in the rest frame of

0136   FourVector pVecRecoil;     // fluid cell and lab frame

0137   FourVector xVec;           // 4-vector of position (for next time step!)

0138   double velocity_jet[4];    // jet velocity for MATTER

0139   double eta;                // pseudo-rapidity

0140   double SpatialRapidity;    // space-time rapidity

0141   bool coherent;             // whether particle is coherent with its

0142                              // mother or daughter.

0143       // while coherent, additional radiation is prohibited.

0144   int sibling;         // counter parton in radiation process

0145   FourVector pAtSplit; // mother's momentum when splitting happens

0146   bool userInfo;       // whether user information is stored

0147   bool active;         // wheter this parton is active in MARTINI

0148
0149   // flow info

0150   double vx, vy, vz;   // 3 components of flow velocity

0151   double T;            // Temperature of fluid cell

0152   double beta, gamma;  // flow velocity & gamma factor

0153   double cosPhi;       // angle between flow and particle

0154   double cosPhiRest;   // angle between flow and particle in rest frame

0155   double boostBack;    // factor for boosting back to lab frame

0156   double cosPhiRestEl; // angle between flow and scat. particle in rest frame

0157   double boostBackEl;
0158
0159   for (int i = 0; i < pIn.size(); i++) {
0160
0161     Id = pIn[i].pid();
0162     pStat = pIn[i].pstat();
0163     // do nothing for negative particles

0164     if (pStat < 0)
0165       continue;
0166
0167     px = pIn[i].px();
0168     py = pIn[i].py();
0169     pz = pIn[i].pz();
0170     p0 = pIn[i].e();
0171
0172     pAbs = sqrt(px * px + py * py + pz * pz);
0173
0174     // In MARTINI, particles are all massless and on-shell

0175     pVec = FourVector(px, py, pz, pAbs);
0176
0177     velx = px / p0;
0178     vely = py / p0;
0179     velz = pz / p0;
0180     double velocityMod = std::sqrt(std::pow(velx, 2) + std::pow(vely, 2) +
0181                                    std::pow(velz, 2));
0182     tt = Time;
0183     xx = pIn[i].x_in().x() + (Time - pIn[i].x_in().t()) * velx / velocityMod;
0184     yy = pIn[i].x_in().y() + (Time - pIn[i].x_in().t()) * vely / velocityMod;
0185     zz = pIn[i].x_in().z() + (Time - pIn[i].x_in().t()) * velz / velocityMod;
0186
0187     eta = pIn[i].eta();
0188     SpatialRapidity = 0.5 * std::log((tt + zz) / (tt - zz));
0189     double boostedTStart = tStart * cosh(SpatialRapidity);
0190
0191     // Extract fluid properties

0192     std::unique_ptr<FluidCellInfo> check_fluid_info_ptr;
0193     GetHydroCellSignal(Time, xx, yy, zz, check_fluid_info_ptr);
0194     VERBOSE(8) << MAGENTA << "Temperature from Brick (Signal) = "
0195                << check_fluid_info_ptr->temperature;
0196
0197     vx = check_fluid_info_ptr->vx;
0198     vy = check_fluid_info_ptr->vy;
0199     vz = check_fluid_info_ptr->vz;
0200     T = check_fluid_info_ptr->temperature;
0201
0202     beta = sqrt(vx * vx + vy * vy + vz * vz);
0203
0204     // Only accept low t particles

0205     if (pIn[i].t() > Q0 * Q0 + rounding_error || Time <= boostedTStart ||
0206         T < hydro_Tc)
0207       continue;
0208     TakeResponsibilityFor(
0209         pIn[i]); // Generate error if another module already has responsibility.

0210
0211     pLabel = pIn[i].plabel();
0212     if (pLabel == 0) {
0213       IncrementpLable();
0214       pIn[i].set_label(pLabelNew);
0215       pLabel = pLabelNew;
0216     }
0217
0218     userInfo = pIn[i].has_user_info<MARTINIUserInfo>();
0219     // set user information to every parton

0220     if (!userInfo)
0221       pIn[i].set_user_info(new MARTINIUserInfo());
0222
0223     coherent = pIn[i].user_info<MARTINIUserInfo>().coherent();
0224     sibling = pIn[i].user_info<MARTINIUserInfo>().coherent_with();
0225     pAtSplit = pIn[i].user_info<MARTINIUserInfo>().p_at_split();
0226
0227     // Set momentum in fluid cell's frame

0228     if (beta < 1e-10) {
0229       // 1: for static medium

0230       gamma = 1.;
0231       cosPhi = 1.;
0232       pRest = pAbs;
0233       pVecRest = pVec;
0234
0235       cosPhiRest = 1.;
0236       boostBack = 1.;
0237     } else {
0238       // 2: for evolving medium

0239       gamma = 1. / sqrt(1. - beta * beta);
0240       cosPhi = (px * vx + py * vy + pz * vz) / (pAbs * beta);
0241
0242       // boost particle to the local rest frame of fluid cell

0243       pRest = pAbs * gamma * (1. - beta * cosPhi);
0244
0245       pxRest = -vx * gamma * pAbs +
0246                (1. + (gamma - 1.) * vx * vx / (beta * beta)) * px +
0247                (gamma - 1.) * vx * vy / (beta * beta) * py +
0248                (gamma - 1.) * vx * vz / (beta * beta) * pz;
0249       pyRest = -vy * gamma * pAbs +
0250                (1. + (gamma - 1.) * vy * vy / (beta * beta)) * py +
0251                (gamma - 1.) * vx * vy / (beta * beta) * px +
0252                (gamma - 1.) * vy * vz / (beta * beta) * pz;
0253       pzRest = -vz * gamma * pAbs +
0254                (1. + (gamma - 1.) * vz * vz / (beta * beta)) * pz +
0255                (gamma - 1.) * vx * vz / (beta * beta) * px +
0256                (gamma - 1.) * vy * vz / (beta * beta) * py;
0257
0258       pVecRest = FourVector(pxRest, pyRest, pzRest, pRest);
0259
0260       cosPhiRest = (pxRest * vx + pyRest * vy + pzRest * vz) / (pRest * beta);
0261       boostBack = gamma * (1. + beta * cosPhiRest);
0262     }
0263
0264     // free-streaming for too soft partons

0265     if (pRest < eLossCut) continue;
0266
0267     xVec = FourVector(xx + px / pAbs * deltaT, yy + py / pAbs * deltaT,
0268                       zz + pz / pAbs * deltaT, Time + deltaT);
0269
0270     velocity_jet[0] = 1.0;
0271     velocity_jet[1] = pIn[i].jet_v().x();
0272     velocity_jet[2] = pIn[i].jet_v().y();
0273     velocity_jet[3] = pIn[i].jet_v().z();
0274
0275     // set alpha_s and g

0276     alpha_s = alpha_s0;
0277     if (run_alphas) {
0278       double mu = sqrt(2.*pAbs*T);
0279       alpha_s = RunningAlphaS(mu);
0280     }
0281     g = sqrt(4. * M_PI * alpha_s);
0282
0283     double deltaTRest = deltaT / gamma;
0284     // determine the energy-loss process

0285     int process = DetermineProcess(pRest, T, deltaTRest, Id);
0286
0287     //cout << "process = " << process << endl;

0288
0289     VERBOSE(8) << MAGENTA << "Time = " << Time << " Id = " << Id
0290                << " process = " << process << " T = " << setprecision(3) << T
0291                << " pAbs = " << pAbs << " " << px << " " << py << " " << pz
0292                << " | pRest = " << pRest << "/" << pcut
0293                << " | position = " << xx << " " << yy << " " << zz
0294                << " | stat = " << pStat << " " << pLabel << " " << coherent
0295                << " | " << pAtSplit.x();
0296
0297     // update the status of parton if this parton is still coherent with its sibling

0298     bool coherent_temp = coherent;
0299     // de-activate formation time of in-medium radiation at this moment..

0300     //if ( coherent ) coherent = isCoherent(pIn[i], sibling, T);

0301     //if ( coherent_temp != coherent )

0302     //  JSDEBUG << "Coherent status changed: " << coherent_temp << " -> " << coherent;

0303
0304     // Do nothing for this parton at this timestep

0305     // If this parton is in coherent state with its sibling,

0306     // further radiation is prohibited until it becomes de-coherent

0307     if (process == 0 || (coherent && (process < 5))) {
0308       pOut.push_back(Parton(pLabel, Id, pStat, pVec, xVec));
0309       pOut[pOut.size() - 1].set_form_time(0.);
0310       pOut[pOut.size() - 1].set_jet_v(velocity_jet); // use initial jet velocity

0311
0312       if (coherent) {
0313         pOut[pOut.size() - 1].set_user_info(
0314             new MARTINIUserInfo(coherent, sibling, pAtSplit));
0315       } else {
0316         pOut[pOut.size() - 1].set_user_info(new MARTINIUserInfo());
0317       }
0318
0319       return;
0320     }
0321
0322     if (std::abs(Id) == 1 || std::abs(Id) == 2 || std::abs(Id) == 3) {
0323
0324       // quark radiating gluon (q->qg)

0325       if (process == 1) {
0326         if (pRest / T < AMYpCut)
0327           return;
0328
0329         // sample radiated parton's momentum

0330         kRest = getNewMomentumRad(pRest, T, process);
0331         if (kRest > pRest)
0332           return;
0333
0334         // final state parton's momentum

0335         pNewRest = pRest - kRest;
0336
0337         pNew = pNewRest * boostBack;
0338         pVecNew.Set((px / pAbs) * pNew, (py / pAbs) * pNew, (pz / pAbs) * pNew,
0339                     pNew);
0340         pOut.push_back(Parton(pLabel, Id, pStat, pVecNew, xVec));
0341         pOut[pOut.size() - 1].set_form_time(0.);
0342         pOut[pOut.size() - 1].set_jet_v(
0343             velocity_jet); // use initial jet velocity

0344         pOut[pOut.size() - 1].set_user_info(new MARTINIUserInfo());
0345
0346         if (kRest > pcut) {
0347           IncrementpLable();
0348           k = kRest * boostBack;
0349           kVec.Set((px / pAbs) * k, (py / pAbs) * k, (pz / pAbs) * k, k);
0350           pOut.push_back(Parton(pLabelNew, 21, pStat, kVec, xVec));
0351           pOut[pOut.size() - 1].set_form_time(0.);
0352           pOut[pOut.size() - 1].set_jet_v(
0353               velocity_jet); // use initial jet velocity

0354         }
0355
0356         //// If original and radiated partons are above pcut, set them coherent

0357         //if (pOut.size() == 2) {

0358         //  bool coherentNew = true;

0359         //  FourVector pAtSplitNew = pOut[0].p_in();

0360         //  pOut[0].set_user_info(new MARTINIUserInfo(coherentNew, pLabelNew, pAtSplitNew));

0361         //  pOut[1].set_user_info(new MARTINIUserInfo(coherentNew, pLabel, pAtSplitNew));

0362         //}

0363
0364         return;
0365       } else if (process == 2) {
0366         // quark radiating photon (q->qgamma)

0367         if (pRest / T < AMYpCut)
0368           return;
0369
0370         // sample radiated parton's momentum

0371         kRest = getNewMomentumRad(pRest, T, process);
0372         if (kRest > pRest)
0373           return;
0374
0375         // final state parton's momentum

0376         pNewRest = pRest - kRest;
0377
0378         pNew = pNewRest * boostBack;
0379         pVecNew.Set((px / pAbs) * pNew, (py / pAbs) * pNew, (pz / pAbs) * pNew,
0380                     pNew);
0381         pOut.push_back(Parton(pLabel, Id, pStat, pVecNew, xVec));
0382         pOut[pOut.size() - 1].set_form_time(0.);
0383         pOut[pOut.size() - 1].set_jet_v(
0384             velocity_jet); // use initial jet velocity

0385         pOut[pOut.size() - 1].set_user_info(new MARTINIUserInfo());
0386
0387         // photon doesn't have energy threshold; No absorption into medium

0388         // However, we only keep positive energy photons

0389         if (kRest > 0.) {
0390           IncrementpLable();
0391           k = kRest * boostBack;
0392           kVec.Set((px / pAbs) * k, (py / pAbs) * k, (pz / pAbs) * k, k);
0393           pOut.push_back(Parton(pLabelNew, 21, pStat, kVec, xVec));
0394           pOut[pOut.size() - 1].set_form_time(0.);
0395           pOut[pOut.size() - 1].set_jet_v(
0396               velocity_jet); // use initial jet velocity

0397         }
0398
0399         return;
0400       } else if (process == 5 || process == 6) {
0401         // quark scattering with either quark (qq->qq) or gluon (qg->qg)

0402         omega = getEnergyTransfer(pRest, T, process);
0403         q = getMomentumTransfer(pRest, omega, T, process);
0404
0405         // momentum transfer is always space-like

0406         if (q < fabs(omega))
0407           return;
0408
0409         pVecNewRest = getNewMomentumElas(pVecRest, omega, q);
0410
0411         pNewRest = pVecNewRest.t();
0412
0413         // Boost scattered particle to lab frame

0414         if (beta < 1e-10) {
0415           // 1: for static medium

0416           pVecNew = pVecNewRest;
0417         } else {
0418           // 2: for evolving medium

0419           pxNew =
0420               vx * gamma * pVecNewRest.t() +
0421               (1. + (gamma - 1.) * vx * vx / (beta * beta)) * pVecNewRest.x() +
0422               (gamma - 1.) * vx * vy / (beta * beta) * pVecNewRest.y() +
0423               (gamma - 1.) * vx * vz / (beta * beta) * pVecNewRest.z();
0424           pyNew =
0425               vy * gamma * pVecNewRest.t() +
0426               (1. + (gamma - 1.) * vy * vy / (beta * beta)) * pVecNewRest.y() +
0427               (gamma - 1.) * vx * vy / (beta * beta) * pVecNewRest.x() +
0428               (gamma - 1.) * vy * vz / (beta * beta) * pVecNewRest.z();
0429           pzNew =
0430               vz * gamma * pVecNewRest.t() +
0431               (1. + (gamma - 1.) * vz * vz / (beta * beta)) * pVecNewRest.z() +
0432               (gamma - 1.) * vx * vz / (beta * beta) * pVecNewRest.x() +
0433               (gamma - 1.) * vy * vz / (beta * beta) * pVecNewRest.y();
0434
0435           pNew = sqrt(pxNew * pxNew + pyNew * pyNew + pzNew * pzNew);
0436           pVecNew.Set(pxNew, pyNew, pzNew, pNew);
0437         }
0438
0439         pOut.push_back(Parton(pLabel, Id, pStat, pVecNew, xVec));
0440         pOut[pOut.size() - 1].set_form_time(0.);
0441         pOut[pOut.size() - 1].set_jet_v(
0442             velocity_jet); // use initial jet velocity

0443
0444         if (coherent) {
0445           pOut[pOut.size() - 1].set_user_info(
0446               new MARTINIUserInfo(coherent, sibling, pAtSplit));
0447         } else {
0448           pOut[pOut.size() - 1].set_user_info(new MARTINIUserInfo());
0449         }
0450
0451         if (recoil_on) {
0452           // momentum transfer between elastic scattering

0453           qVec = pVecRest;
0454           qVec -= pVecNewRest;
0455
0456           pVecThermalRest = getThermalVec(qVec, T, Id);
0457           pVecRecoilRest = qVec;
0458           pVecRecoilRest += pVecThermalRest;
0459           double pRecoilRest = pVecRecoilRest.t();
0460
0461           if (pRecoilRest > pcut) {
0462
0463             // Boost recoil particle to lab frame

0464             if (beta < 1e-10) {
0465               // 1: for static medium

0466               pVecThermal = pVecThermalRest;
0467               pVecRecoil = pVecRecoilRest;
0468             } else {
0469               // 2: for evolving medium

0470               pxThermal =
0471                   vx * gamma * pVecThermalRest.t() +
0472                   (1. + (gamma - 1.) * vx * vx / (beta * beta)) *
0473                       pVecThermalRest.x() +
0474                   (gamma - 1.) * vx * vy / (beta * beta) * pVecThermalRest.y() +
0475                   (gamma - 1.) * vx * vz / (beta * beta) * pVecThermalRest.z();
0476               pyThermal =
0477                   vy * gamma * pVecThermalRest.t() +
0478                   (1. + (gamma - 1.) * vy * vy / (beta * beta)) *
0479                       pVecThermalRest.y() +
0480                   (gamma - 1.) * vx * vy / (beta * beta) * pVecThermalRest.x() +
0481                   (gamma - 1.) * vy * vz / (beta * beta) * pVecThermalRest.z();
0482               pzThermal =
0483                   vz * gamma * pVecThermalRest.t() +
0484                   (1. + (gamma - 1.) * vz * vz / (beta * beta)) *
0485                       pVecThermalRest.z() +
0486                   (gamma - 1.) * vx * vz / (beta * beta) * pVecThermalRest.x() +
0487                   (gamma - 1.) * vy * vz / (beta * beta) * pVecThermalRest.y();
0488
0489               pThermal = sqrt(pxThermal * pxThermal + pyThermal * pyThermal +
0490                               pzThermal * pzThermal);
0491               pVecThermal.Set(pxThermal, pyThermal, pzThermal, pThermal);
0492
0493               pxRecoil =
0494                   vx * gamma * pVecRecoilRest.t() +
0495                   (1. + (gamma - 1.) * vx * vx / (beta * beta)) *
0496                       pVecRecoilRest.x() +
0497                   (gamma - 1.) * vx * vy / (beta * beta) * pVecRecoilRest.y() +
0498                   (gamma - 1.) * vx * vz / (beta * beta) * pVecRecoilRest.z();
0499               pyRecoil =
0500                   vy * gamma * pVecRecoilRest.t() +
0501                   (1. + (gamma - 1.) * vy * vy / (beta * beta)) *
0502                       pVecRecoilRest.y() +
0503                   (gamma - 1.) * vx * vy / (beta * beta) * pVecRecoilRest.x() +
0504                   (gamma - 1.) * vy * vz / (beta * beta) * pVecRecoilRest.z();
0505               pzRecoil =
0506                   vz * gamma * pVecRecoilRest.t() +
0507                   (1. + (gamma - 1.) * vz * vz / (beta * beta)) *
0508                       pVecRecoilRest.z() +
0509                   (gamma - 1.) * vx * vz / (beta * beta) * pVecRecoilRest.x() +
0510                   (gamma - 1.) * vy * vz / (beta * beta) * pVecRecoilRest.y();
0511
0512               pRecoil = sqrt(pxRecoil * pxRecoil + pyRecoil * pyRecoil +
0513                              pzRecoil * pzRecoil);
0514               pVecRecoil.Set(pxRecoil, pyRecoil, pzRecoil, pRecoil);
0515             }
0516
0517             // determine id of recoil parton

0518             if (process == 5) {
0519               // choose the Id of new qqbar pair. Note that we only deal with nf = 3

0520               double r = ZeroOneDistribution(*GetMt19937Generator());
0521               if (r < 1. / 3.)
0522                 newId = 1;
0523               else if (r < 2. / 3.)
0524                 newId = 2;
0525               else
0526                 newId = 3;
0527             } else {
0528               newId = 21;
0529             }
0530
0531             if (pVecThermal.t() > 1e-5) {
0532               IncrementpLable();
0533               pOut.push_back(Parton(pLabelNew, newId, -1, pVecThermal, xVec));
0534               pOut[pOut.size() - 1].set_form_time(0.);
0535               pOut[pOut.size() - 1].set_jet_v(
0536                   velocity_jet); // use initial jet velocity

0537               //cout << "process == 5&6::newLabel:" << pLabelNew << endl;

0538             }
0539
0540             IncrementpLable();
0541             pOut.push_back(Parton(pLabelNew, newId, 1, pVecRecoil, xVec));
0542             pOut[pOut.size() - 1].set_form_time(0.);
0543             pOut[pOut.size() - 1].set_jet_v(
0544                 velocity_jet); // use initial jet velocity

0545             //cout << "process == 5&6::newLabel:" << pLabelNew << endl;

0546           }
0547         }
0548         return;
0549       } else if (process == 9) {
0550         // quark converting to gluon

0551         pOut.push_back(Parton(pLabel, 21, pStat, pVec, xVec));
0552         pOut[pOut.size() - 1].set_form_time(0.);
0553         pOut[pOut.size() - 1].set_jet_v(
0554             velocity_jet); // use initial jet velocity

0555
0556         if (coherent) {
0557           pOut[pOut.size() - 1].set_user_info(
0558               new MARTINIUserInfo(coherent, sibling, pAtSplit));
0559         } else {
0560           pOut[pOut.size() - 1].set_user_info(new MARTINIUserInfo());
0561         }
0562
0563         return;
0564       } else if (process == 10) {
0565         // quark converting to photon

0566         pOut.push_back(Parton(pLabel, 22, pStat, pVec, xVec));
0567         pOut[pOut.size() - 1].set_form_time(0.);
0568         pOut[pOut.size() - 1].set_jet_v(
0569             velocity_jet); // use initial jet velocity

0570
0571         return;
0572       }
0573     } else if (Id == 21) {
0574
0575       // gluon radiating gluon (g->gg)

0576       if (process == 3) {
0577         if (pRest / T < AMYpCut)
0578           return;
0579
0580         // sample radiated parton's momentum

0581         kRest = getNewMomentumRad(pRest, T, process);
0582         if (kRest > pRest)
0583           return;
0584
0585         // final state parton's momentum

0586         pNewRest = pRest - kRest;
0587
0588         pNew = pNewRest * boostBack;
0589         pVecNew.Set((px / pAbs) * pNew, (py / pAbs) * pNew, (pz / pAbs) * pNew,
0590                     pNew);
0591         pOut.push_back(Parton(pLabel, Id, pStat, pVecNew, xVec));
0592         pOut[pOut.size() - 1].set_form_time(0.);
0593         pOut[pOut.size() - 1].set_jet_v(
0594             velocity_jet); // use initial jet velocity

0595         pOut[pOut.size() - 1].set_user_info(new MARTINIUserInfo());
0596
0597         if (kRest > pcut) {
0598           IncrementpLable();
0599           k = kRest * boostBack;
0600           kVec.Set((px / pAbs) * k, (py / pAbs) * k, (pz / pAbs) * k, k);
0601           pOut.push_back(Parton(pLabelNew, Id, pStat, kVec, xVec));
0602           pOut[pOut.size() - 1].set_form_time(0.);
0603           pOut[pOut.size() - 1].set_jet_v(
0604               velocity_jet); // use initial jet velocity

0605         }
0606
0607         //// If original and radiated partons are above pcut, set them coherent

0608         //if (pOut.size() == 2) {

0609         //  bool coherentNew = true;

0610         //  FourVector pAtSplitNew = pOut[0].p_in();

0611         //  pOut[0].set_user_info(new MARTINIUserInfo(coherentNew, pLabelNew, pAtSplitNew));

0612         //  pOut[1].set_user_info(new MARTINIUserInfo(coherentNew, pLabel, pAtSplitNew));

0613         //}

0614
0615         return;
0616       } else if (process == 4) {
0617         // gluon split into quark-antiquark pair (g->qqbar)

0618         if (pRest / T < AMYpCut)
0619           return;
0620
0621         // sample radiated parton's momentum

0622         kRest = getNewMomentumRad(pRest, T, process);
0623         if (kRest > pRest)
0624           return;
0625
0626         // final state parton's momentum

0627         pNewRest = pRest - kRest;
0628
0629         // choose the Id of new qqbar pair. Note that we only deal with nf = 3

0630         double r = ZeroOneDistribution(*GetMt19937Generator());
0631         if (r < 1. / 3.)
0632           newId = 1;
0633         else if (r < 2. / 3.)
0634           newId = 2;
0635         else
0636           newId = 3;
0637
0638         // *momentum of quark is usually larger than that of anti-quark

0639         pNew = pNewRest * boostBack;
0640         pVecNew.Set((px / pAbs) * pNew, (py / pAbs) * pNew, (pz / pAbs) * pNew,
0641                     pNew);
0642         pOut.push_back(Parton(pLabel, newId, pStat, pVecNew, xVec));
0643         pOut[pOut.size() - 1].set_form_time(0.);
0644         pOut[pOut.size() - 1].set_jet_v(
0645             velocity_jet); // use initial jet velocity

0646         pOut[pOut.size() - 1].set_user_info(new MARTINIUserInfo());
0647
0648         if (kRest > pcut) {
0649           IncrementpLable();
0650           k = kRest * boostBack;
0651           kVec.Set((px / pAbs) * k, (py / pAbs) * k, (pz / pAbs) * k, k);
0652           pOut.push_back(Parton(pLabelNew, -newId, pStat, kVec, xVec));
0653           pOut[pOut.size() - 1].set_form_time(0.);
0654           pOut[pOut.size() - 1].set_jet_v(
0655               velocity_jet); // use initial jet velocity

0656         }
0657
0658         //// If original and radiated partons are above pcut, set them coherent

0659         //if (pOut.size() == 2) {

0660         //  bool coherentNew = true;

0661         //  FourVector pAtSplitNew = pOut[0].p_in();

0662         //  pOut[0].set_user_info(new MARTINIUserInfo(coherentNew, pLabelNew, pAtSplitNew));

0663         //  pOut[1].set_user_info(new MARTINIUserInfo(coherentNew, pLabel, pAtSplitNew));

0664         //}

0665
0666         return;
0667       } else if (process == 7 || process == 8) {
0668         // gluon scattering with either quark (gq->gq) or gluon (gg->gg)

0669         omega = getEnergyTransfer(pRest, T, process);
0670         q = getMomentumTransfer(pRest, omega, T, process);
0671
0672         // momentum transfer is always space-like

0673         if (q < fabs(omega))
0674           return;
0675
0676         pVecNewRest = getNewMomentumElas(pVecRest, omega, q);
0677
0678         pNewRest = pVecNewRest.t();
0679
0680         // Boost scattered particle to lab frame

0681         if (beta < 1e-10) {
0682           // 1: for static medium

0683           pVecNew = pVecNewRest;
0684         } else {
0685           // 2: for evolving medium

0686           pxNew =
0687               vx * gamma * pVecNewRest.t() +
0688               (1. + (gamma - 1.) * vx * vx / (beta * beta)) * pVecNewRest.x() +
0689               (gamma - 1.) * vx * vy / (beta * beta) * pVecNewRest.y() +
0690               (gamma - 1.) * vx * vz / (beta * beta) * pVecNewRest.z();
0691           pyNew =
0692               vy * gamma * pVecNewRest.t() +
0693               (1. + (gamma - 1.) * vy * vy / (beta * beta)) * pVecNewRest.y() +
0694               (gamma - 1.) * vx * vy / (beta * beta) * pVecNewRest.x() +
0695               (gamma - 1.) * vy * vz / (beta * beta) * pVecNewRest.z();
0696           pzNew =
0697               vz * gamma * pVecNewRest.t() +
0698               (1. + (gamma - 1.) * vz * vz / (beta * beta)) * pVecNewRest.z() +
0699               (gamma - 1.) * vx * vz / (beta * beta) * pVecNewRest.x() +
0700               (gamma - 1.) * vy * vz / (beta * beta) * pVecNewRest.y();
0701
0702           pNew = sqrt(pxNew * pxNew + pyNew * pyNew + pzNew * pzNew);
0703           pVecNew.Set(pxNew, pyNew, pzNew, pNew);
0704         }
0705
0706         pOut.push_back(Parton(pLabel, Id, pStat, pVecNew, xVec));
0707         pOut[pOut.size() - 1].set_form_time(0.);
0708         pOut[pOut.size() - 1].set_jet_v(
0709             velocity_jet); // use initial jet velocity

0710
0711         if (coherent) {
0712           pOut[pOut.size() - 1].set_user_info(
0713               new MARTINIUserInfo(coherent, sibling, pAtSplit));
0714         } else {
0715           pOut[pOut.size() - 1].set_user_info(new MARTINIUserInfo());
0716         }
0717
0718         if (recoil_on) {
0719           // momentum transfer between elastic scattering

0720           qVec = pVecRest;
0721           qVec -= pVecNewRest;
0722
0723           pVecThermalRest = getThermalVec(qVec, T, Id);
0724           pVecRecoilRest = qVec;
0725           pVecRecoilRest += pVecThermalRest;
0726           double pRecoilRest = pVecRecoilRest.t();
0727
0728           if (pRecoilRest > pcut) {
0729
0730             // Boost recoil particle to lab frame

0731             if (beta < 1e-10) {
0732               // 1: for static medium

0733               pVecThermal = pVecThermalRest;
0734               pVecRecoil = pVecRecoilRest;
0735             } else {
0736               // 2: for evolving medium

0737               pxThermal =
0738                   vx * gamma * pVecThermalRest.t() +
0739                   (1. + (gamma - 1.) * vx * vx / (beta * beta)) *
0740                       pVecThermalRest.x() +
0741                   (gamma - 1.) * vx * vy / (beta * beta) * pVecThermalRest.y() +
0742                   (gamma - 1.) * vx * vz / (beta * beta) * pVecThermalRest.z();
0743               pyThermal =
0744                   vy * gamma * pVecThermalRest.t() +
0745                   (1. + (gamma - 1.) * vy * vy / (beta * beta)) *
0746                       pVecThermalRest.y() +
0747                   (gamma - 1.) * vx * vy / (beta * beta) * pVecThermalRest.x() +
0748                   (gamma - 1.) * vy * vz / (beta * beta) * pVecThermalRest.z();
0749               pzThermal =
0750                   vz * gamma * pVecThermalRest.t() +
0751                   (1. + (gamma - 1.) * vz * vz / (beta * beta)) *
0752                       pVecThermalRest.z() +
0753                   (gamma - 1.) * vx * vz / (beta * beta) * pVecThermalRest.x() +
0754                   (gamma - 1.) * vy * vz / (beta * beta) * pVecThermalRest.y();
0755
0756               pThermal = sqrt(pxThermal * pxThermal + pyThermal * pyThermal +
0757                               pzThermal * pzThermal);
0758               pVecThermal.Set(pxThermal, pyThermal, pzThermal, pThermal);
0759
0760               pxRecoil =
0761                   vx * gamma * pVecRecoilRest.t() +
0762                   (1. + (gamma - 1.) * vx * vx / (beta * beta)) *
0763                       pVecRecoilRest.x() +
0764                   (gamma - 1.) * vx * vy / (beta * beta) * pVecRecoilRest.y() +
0765                   (gamma - 1.) * vx * vz / (beta * beta) * pVecRecoilRest.z();
0766               pyRecoil =
0767                   vy * gamma * pVecRecoilRest.t() +
0768                   (1. + (gamma - 1.) * vy * vy / (beta * beta)) *
0769                       pVecRecoilRest.y() +
0770                   (gamma - 1.) * vx * vy / (beta * beta) * pVecRecoilRest.x() +
0771                   (gamma - 1.) * vy * vz / (beta * beta) * pVecRecoilRest.z();
0772               pzRecoil =
0773                   vz * gamma * pVecRecoilRest.t() +
0774                   (1. + (gamma - 1.) * vz * vz / (beta * beta)) *
0775                       pVecRecoilRest.z() +
0776                   (gamma - 1.) * vx * vz / (beta * beta) * pVecRecoilRest.x() +
0777                   (gamma - 1.) * vy * vz / (beta * beta) * pVecRecoilRest.y();
0778
0779               pRecoil = sqrt(pxRecoil * pxRecoil + pyRecoil * pyRecoil +
0780                              pzRecoil * pzRecoil);
0781               pVecRecoil.Set(pxRecoil, pyRecoil, pzRecoil, pRecoil);
0782             }
0783
0784             // determine id of recoil parton

0785             if (process == 7) {
0786               // choose the Id of new qqbar pair. Note that we only deal with nf = 3

0787               double r = ZeroOneDistribution(*GetMt19937Generator());
0788               if (r < 1. / 3.)
0789                 newId = 1;
0790               else if (r < 2. / 3.)
0791                 newId = 2;
0792               else
0793                 newId = 3;
0794             } else {
0795               newId = 21;
0796             }
0797
0798             if (pVecThermal.t() > 1e-5) {
0799               IncrementpLable();
0800               pOut.push_back(Parton(pLabelNew, newId, -1, pVecThermal, xVec));
0801               pOut[pOut.size() - 1].set_form_time(0.);
0802               pOut[pOut.size() - 1].set_jet_v(
0803                   velocity_jet); // use initial jet velocity

0804               //cout << "process == 7&8::newLabel:" << pLabelNew << endl;

0805             }
0806
0807             IncrementpLable();
0808             pOut.push_back(Parton(pLabelNew, newId, 1, pVecRecoil, xVec));
0809             pOut[pOut.size() - 1].set_form_time(0.);
0810             pOut[pOut.size() - 1].set_jet_v(
0811                 velocity_jet); // use initial jet velocity

0812             //cout << "process == 7&8::newLabel:" << pLabelNew << endl;

0813           }
0814         }
0815
0816         return;
0817       } else if (process == 11) {
0818         // gluon converting to quark

0819         // choose the Id of new qqbar pair. Note that we only deal with nf = 3

0820         double r = ZeroOneDistribution(*GetMt19937Generator());
0821         if (r < 1. / 3.)
0822           newId = 1;
0823         else if (r < 2. / 3.)
0824           newId = 2;
0825         else
0826           newId = 3;
0827
0828         double antiquark = ZeroOneDistribution(*GetMt19937Generator());
0829         if (antiquark < 0.5)
0830           newId *= -1;
0831
0832         pOut.push_back(Parton(pLabel, newId, pStat, pVec, xVec));
0833         pOut[pOut.size() - 1].set_form_time(0.);
0834         pOut[pOut.size() - 1].set_jet_v(
0835             velocity_jet); // use initial jet velocity

0836
0837         if (coherent) {
0838           pOut[pOut.size() - 1].set_user_info(
0839               new MARTINIUserInfo(coherent, sibling, pAtSplit));
0840         } else {
0841           pOut[pOut.size() - 1].set_user_info(new MARTINIUserInfo());
0842         }
0843
0844         return;
0845       }
0846     } // Id==21

0847   }   // particle loop

0848 }
0849
0850 double Martini::RunningAlphaS(double mu) {
0851
0852   if (mu < 1.0) return alpha_s0;
0853
0854   double beta0 = 9./(4.*M_PI);
0855   double LambdaQCD = exp( -1./(2.*beta0*alpha_s0) );
0856   double running_alphas = 1./(beta0 * log( pow(mu/LambdaQCD, 2.) ));
0857   return running_alphas;
0858 }
0859
0860 int Martini::DetermineProcess(double pRest, double T, double deltaTRest,
0861                               int Id) {
0862
0863   double dT = deltaTRest / hbarc; // time step in [GeV^(-1)]

0864
0865   // get the rates for each process

0866   // total Probability = dT*Rate

0867   RateRadiative rateRad;
0868   rateRad = getRateRadTotal(pRest, T);
0869   RateElastic rateElas;
0870   rateElas = getRateElasTotal(pRest, T);
0871   RateConversion rateConv;
0872   rateConv = getRateConv(pRest, T);
0873
0874   // evolution for quark (u, d, s)

0875   if (std::abs(Id) == 1 || std::abs(Id) == 2 || std::abs(Id) == 3) {
0876     // multiplying by (ef/e)^2

0877     if (abs(Id) == 1) {
0878       rateRad.qqgamma *= 4. / 9.;
0879       rateConv.qgamma *= 4. / 9.;
0880     } else {
0881       rateRad.qqgamma *= 1. / 9.;
0882       rateConv.qgamma *= 1. / 9.;
0883     }
0884
0885     double totalQuarkProb = 0.;
0886
0887     /* block the photon process at this moment */
0888     //if (pRest/T > AMYpCut) totalQuarkProb += (rateRad.qqg + rateRad.qqgamma)*dT;

0889     //totalQuarkProb += (rateElas.qq + rateElas.qg + rateConv.qg + rateConv.qgamma)*dT;

0890
0891     if (pRest / T > AMYpCut)
0892       totalQuarkProb += rateRad.qqg * dT;
0893     totalQuarkProb += (rateElas.qq + rateElas.qg + rateConv.qg) * dT;
0894
0895     // warn if total probability exceeds 1

0896     if (totalQuarkProb > 1.) {
0897       JSWARN << " : Total Probability for quark processes exceeds 1 ("
0898              << totalQuarkProb << "). "
0899              << " : Most likely this means you should choose a smaller deltaT "
0900                 "in the xml"
0901              << " (e.g. 0.01). "
0902              << "pRest=" << pRest << " T=" << T << " pRest/T=" << pRest / T;
0903       //throw std::runtime_error ("Martini probability problem.");

0904     }
0905
0906     double accumProb = 0.;
0907     double nextProb = 0.;
0908     double Prob = 0.;
0909
0910     if (ZeroOneDistribution(*GetMt19937Generator()) < totalQuarkProb) {
0911       /* label for process

0912          [1-4  : Radiation ] 1: q->qg , 2 : q->qgamma, 3 : g->gg , 4: g->qqbar

0913          [5-8  : Elastic   ] 5: qq->qq, 6 : qg->qg   , 7 : gq->gq, 8: gg->gg

0914          [9-11 : Conversion] 9: q->g  , 10: q->gamma , 11: g->q                */
0915       double randProb = ZeroOneDistribution(*GetMt19937Generator());
0916
0917       // AMY radiation only happens if energy scale is above certain threshold.

0918       // but elastic/conversion processes doesn't have threshold.

0919       if (pRest / T > AMYpCut) {
0920         Prob = rateRad.qqg * dT / totalQuarkProb;
0921         if (accumProb <= randProb && randProb < (accumProb + Prob))
0922           return 1;
0923
0924         accumProb += Prob;
0925         Prob = rateRad.qqgamma * dT / totalQuarkProb;
0926         if (accumProb <= randProb && randProb < (accumProb + Prob))
0927           return 2;
0928       }
0929
0930       accumProb += Prob;
0931       Prob = rateElas.qq * dT / totalQuarkProb;
0932       if (accumProb <= randProb && randProb < (accumProb + Prob))
0933         return 5;
0934
0935       accumProb += Prob;
0936       Prob = rateElas.qg * dT / totalQuarkProb;
0937       if (accumProb <= randProb && randProb < (accumProb + Prob))
0938         return 6;
0939
0940       accumProb += Prob;
0941       Prob = rateConv.qg * dT / totalQuarkProb;
0942       if (accumProb <= randProb && randProb < (accumProb + Prob))
0943         return 9;
0944
0945       accumProb += Prob;
0946       Prob = rateConv.qgamma * dT / totalQuarkProb;
0947       if (accumProb <= randProb && randProb < (accumProb + Prob))
0948         return 10;
0949     } else {
0950       // nothing happens to quark

0951       return 0;
0952     }
0953   } else if (Id == 21) {
0954     // evolution for gluon

0955     double totalGluonProb = 0.;
0956
0957     if (pRest / T > AMYpCut)
0958       totalGluonProb += (rateRad.ggg + rateRad.gqq) * dT;
0959     totalGluonProb += (rateElas.gq + rateElas.gg + rateConv.gq) * dT;
0960
0961     // warn if total probability exceeds 1

0962     if (totalGluonProb > 1.) {
0963       JSWARN << " : Total Probability for gluon processes exceeds 1 ("
0964              << totalGluonProb << "). "
0965              << " : Most likely this means you should choose a smaller deltaT "
0966                 "in the xml"
0967              << " (e.g. 0.01). "
0968              << "pRest=" << pRest << " T=" << T << " pRest/T=" << pRest / T;
0969       //throw std::runtime_error ("Martini probability problem.");

0970     }
0971
0972     double accumProb = 0.;
0973     double nextProb = 0.;
0974     double Prob = 0.;
0975
0976     if (ZeroOneDistribution(*GetMt19937Generator()) < totalGluonProb) {
0977       /* label for process

0978          [1-4  : Radiation ] 1: q->qg, 2 : q->qgamma, 3 : g->gg, 4: g->qq

0979          [5-8  : Elastic   ] 5: q->q , 6 : q->g     , 7 : g->q , 8: g->g

0980          [9-11 : Conversion] 9: q->g , 10: q->gamma , 11: g->q            */
0981       double randProb = ZeroOneDistribution(*GetMt19937Generator());
0982
0983       // AMY radiation only happens if energy scale is above certain threshold.

0984       // but elastic/conversion processes doesn't have threshold.

0985       if (pRest / T > AMYpCut) {
0986         Prob = rateRad.ggg * dT / totalGluonProb;
0987         if (accumProb <= randProb && randProb < (accumProb + Prob))
0988           return 3;
0989
0990         accumProb += Prob;
0991         Prob = rateRad.gqq * dT / totalGluonProb;
0992         if (accumProb <= randProb && randProb < (accumProb + Prob))
0993           return 4;
0994       }
0995
0996       accumProb += Prob;
0997       Prob = rateElas.gq * dT / totalGluonProb;
0998       if (accumProb <= randProb && randProb < (accumProb + Prob))
0999         return 7;
1000
1001       accumProb += Prob;
1002       Prob = rateElas.gg * dT / totalGluonProb;
1003       if (accumProb <= randProb && randProb < (accumProb + Prob))
1004         return 8;
1005
1006       accumProb += Prob;
1007       Prob = rateConv.gq * dT / totalGluonProb;
1008       if (accumProb <= randProb && randProb < (accumProb + Prob))
1009         return 11;
1010     } else {
1011       // nothing happens to gluon

1012       return 0;
1013     }
1014   }
1015
1016   // if parton is other than u,d,s,g, do nothing

1017   return 0;
1018 }
1019
1020 RateRadiative Martini::getRateRadTotal(double pRest, double T) {
1021   RateRadiative rate;
1022
1023   double u = pRest / T; // making arguments in log to be dimensionless

1024
1025   if (u > AMYpCut) {
1026     rate.qqg = (0.8616 - 3.2913 / (u * u) + 2.1102 / u - 0.9485 / sqrt(u)) *
1027                pow(g, 4.) * T;
1028     rate.ggg = (1.9463 + 61.7856 / (u * u * u) - 30.7877 / (u * u) +
1029                 8.0409 / u - 2.6249 / sqrt(u)) *
1030                pow(g, 4.) * T;
1031     rate.gqq = (2.5830 / (u * u * u) - 1.7010 / (u * u) + 1.4977 / u -
1032                 1.1961 / pow(u, 0.8) + 0.1807 / sqrt(u)) *
1033                pow(g, 4.) * T * nf;
1034     rate.qqgamma =
1035         (0.0053056 + 2.3279 / pow(u, 3.) - 0.6676 / u + 0.3223 / sqrt(u)) *
1036         pow(g, 4.) * alpha_em / alpha_s * T;
1037
1038     double runningFactor =
1039         log(g * T * pow(10., 0.25) / .175) / log(g * T * pow(u, 0.25) / .175);
1040     if (runningFactor < 1.) {
1041       rate.qqg *= runningFactor;
1042       rate.gqq *= runningFactor;
1043       rate.ggg *= runningFactor;
1044       rate.qqgamma *= runningFactor;
1045     }
1046   } else {
1047     rate.qqg = 0.;
1048     rate.ggg = 0.;
1049     rate.gqq = 0.;
1050     rate.qqgamma = 0.;
1051   }
1052
1053   return rate;
1054 }
1055
1056 RateRadiative Martini::getRateRadPos(double u, double T) {
1057   RateRadiative rate;
1058
1059   rate.qqg = (0.5322 - 3.1037 / (u * u) + 2.0139 / u - 0.9417 / sqrt(u)) *
1060              pow(g, 4.) * T;
1061   rate.ggg = (1.1923 - 11.5250 / (u * u * u) + 3.3010 / u - 1.9049 / sqrt(u)) *
1062              pow(g, 4.) * T;
1063   rate.gqq =
1064       (0.0004656 - 0.04621 / (u * u) + 0.0999 / u - 0.08171 / pow(u, 0.8) +
1065        0.008090 / pow(u, 0.2) - 1.2525 * pow(10., -8.) * u) *
1066       pow(g, 4.) * T * nf;
1067   rate.qqgamma = 0.;
1068
1069   return rate;
1070 }
1071
1072 RateRadiative Martini::getRateRadNeg(double u, double T) {
1073   RateRadiative rate;
1074
1075   rate.qqg = (0.3292 - 0.6759 / (u * u) + 0.4871 / pow(u, 1.5) - 0.05393 / u +
1076               0.007878 / sqrt(u)) *
1077              pow(g, 4.) * T;
1078   rate.ggg = (0.7409 + 1.8608 / (u * u * u) - 0.1353 / (u * u) + 0.1401 / u) *
1079              pow(g, 4.) * T;
1080   rate.gqq = (0.03215 / (u * u * u) + 0.01419 / (u * u) + 0.004338 / u -
1081               0.00001246 / sqrt(u)) *
1082              pow(g, 4.) * T * nf;
1083   rate.qqgamma = 0.;
1084
1085   return rate;
1086 }
1087
1088 double Martini::getNewMomentumRad(double pRest, double T, int process) {
1089   double u = pRest / T; // making arguments in log to be dimensionless

1090
1091   double kNew = 0.; // momentum of radiated gluon (dimentionless)

1092   double y;         // kNew candidate

1093   double x;         // random number, uniform on [0,1]

1094   double randA; // uniform random number on [0, Area under the envelop function]

1095   double fy;    // total area under the envelop function

1096   double fyAct; // actual rate

1097
1098   RateRadiative Pos, Neg;
1099   Pos = getRateRadPos(u, T);
1100   Neg = getRateRadNeg(u, T);
1101
1102   // this switch will hold the decision whether k is positive or negative:

1103   // 0 : negative, 1 : positive

1104   int posNegSwitch = 1;
1105
1106   /* process == 1 : quark radiating gluon

1107      process == 2 : quark radiating photon

1108      process == 3 : gluon radiating gluon

1109      process == 4 : gluon split into quark-antiquark pair */
1110
1111   if (process == 1) {
1112     // decide whether k shall be positive or negative

1113     // if x (uniform on [0,1]) < area(k<0)/area(all k) then k < 0

1114     if (ZeroOneDistribution(*GetMt19937Generator()) <
1115         Neg.qqg / (Neg.qqg + Pos.qqg))
1116       posNegSwitch = 0;
1117
1118     if (posNegSwitch == 1) // if k > 0

1119     {
1120       do {
1121         randA = ZeroOneDistribution(*GetMt19937Generator()) *
1122                 area(u + 12., u, posNegSwitch, 1);
1123         y = 2.5 / (LambertW(2.59235 * pow(10., 23.) * exp(-100. * randA)));
1124
1125         fy = 0.025 / (y * y) + 0.01 / y;
1126         fyAct = function(u, y, process);
1127
1128         x = ZeroOneDistribution(*GetMt19937Generator());
1129
1130       } while (x > fyAct / fy);
1131       // reject if x is larger than the ratio fyAct/fy

1132       kNew = y;
1133     } else // if k < 0

1134     {
1135       do {
1136         randA = ZeroOneDistribution(*GetMt19937Generator()) *
1137                 area(-0.05, u, posNegSwitch, 1);
1138         y = -12. / (1. + 480. * randA);
1139
1140         fy = 0.025 / (y * y);
1141         fyAct = function(u, y, process);
1142
1143         x = ZeroOneDistribution(*GetMt19937Generator());
1144
1145       } while (x > fyAct / fy);
1146       // reject if x is larger than the ratio fyAct/fy

1147       kNew = y;
1148     }
1149   } else if (process == 2) {
1150     do {
1151       randA = ZeroOneDistribution(*GetMt19937Generator()) *
1152               area(1.15 * u, u, posNegSwitch, 2);
1153       y = 83895.3 * pow(pow(u, 0.5) * randA, 10. / 3.);
1154
1155       fy = (0.01 / (pow(y, 0.7))) / pow(u, 0.5);
1156       fyAct = function(u, y, process);
1157
1158       x = ZeroOneDistribution(*GetMt19937Generator());
1159
1160     } while (x > fyAct / fy);
1161     // reject if x is larger than the ratio fyAct/fy

1162     kNew = y;
1163   } else if (process == 3) {
1164     // decide whether k shall be positive or negative

1165     // if x (uniform on [0,1]) < area(k<0)/area(all k) then k < 0

1166     if (ZeroOneDistribution(*GetMt19937Generator()) <
1167         Neg.ggg / (Neg.ggg + Pos.ggg))
1168       posNegSwitch = 0;
1169
1170     if (posNegSwitch == 1) // if k > 0

1171     {
1172       do {
1173         randA = ZeroOneDistribution(*GetMt19937Generator()) *
1174                 area(u / 2., u, posNegSwitch, 3);
1175         y = 5. / (LambertW(2.68812 * pow(10., 45.) * exp(-50. * randA)));
1176
1177         fy = 0.1 / (y * y) + 0.02 / y;
1178         fyAct = function(u, y, process);
1179
1180         x = ZeroOneDistribution(*GetMt19937Generator());
1181
1182       } while (x > fyAct / fy);
1183       // reject if x is larger than the ratio fyAct/fy

1184       kNew = y;
1185     } else // if k < 0

1186     {
1187       do {
1188         randA = ZeroOneDistribution(*GetMt19937Generator()) *
1189                 area(-0.05, u, posNegSwitch, 3);
1190         y = -12. / (1. + 120. * randA);
1191
1192         fy = 0.1 / (y * y);
1193         fyAct = function(u, y, process);
1194
1195         x = ZeroOneDistribution(*GetMt19937Generator());
1196
1197       } while (x > fyAct / fy);
1198       // reject if x is larger than the ratio fyAct/fy

1199       kNew = y;
1200     }
1201   } else if (process == 4) {
1202     // decide whether k shall be positive or negative

1203     // if x (uniform on [0,1]) < area(k<0)/area(all k) then k < 0

1204     if (ZeroOneDistribution(*GetMt19937Generator()) <
1205         Neg.gqq / (Neg.gqq + Pos.gqq))
1206       posNegSwitch = 0;
1207
1208     if (posNegSwitch == 1) // if k > 0

1209     {
1210       do {
1211         randA = ZeroOneDistribution(*GetMt19937Generator()) *
1212                 area(u / 2., u, posNegSwitch, 4);
1213         y = 0.83333 * (0.06 * function(u, 0.05, process) + randA) /
1214             function(u, 0.05, process);
1215
1216         fy = 1.2 * function(u, 0.05, process);
1217         fyAct = function(u, y, process);
1218
1219         x = ZeroOneDistribution(*GetMt19937Generator());
1220
1221       } while (x > fyAct / fy);
1222       // reject if x is larger than the ratio fyAct/fy

1223       kNew = y;
1224     } else // if k < 0

1225     {
1226       do {
1227         randA = ZeroOneDistribution(*GetMt19937Generator()) *
1228                 area(-0.05, u, posNegSwitch, 4);
1229         y = (2.5 - u * log(7.81082 * pow(10., -6.) * exp(14.5 / u) +
1230                            (-115.883 + 113.566 * u) * randA)) /
1231             (1. - 0.98 * u);
1232
1233         fy = 0.98 * exp((1. - 1. / u) * (-2.5 + y)) / u;
1234         fyAct = function(u, y, process);
1235
1236         x = ZeroOneDistribution(*GetMt19937Generator());
1237
1238       } while (x > fyAct / fy);
1239       // reject if x is larger than the ratio fyAct/fy

1240       kNew = y;
1241     }
1242   } else {
1243     JSWARN << "Invalid process number (" << process << ")";
1244   }
1245
1246   return kNew * T; // kNew*T is in [GeV]

1247 }
1248
1249 // calculates the area under the envelope function when using the rejection method

1250 // integrals had been solved analytically before

1251 double Martini::area(double y, double u, int posNegSwitch, int process) {
1252   if (process == 1) {
1253     if (posNegSwitch == 1)
1254       return (0.5299 - 0.025 / y + 0.01 * log(y));
1255     else
1256       return (-0.002083 - 0.025 / y);
1257   } else if (process == 2) {
1258     return ((0.03333 * pow(y, 0.3)) / pow(u, 0.5));
1259   } else if (process == 3) {
1260     if (posNegSwitch == 1)
1261       return (2.05991 - 0.1 / y + 0.02 * log(y));
1262     else
1263       return (-0.008333 - 0.1 / y);
1264   } else if (process == 4) {
1265     if (posNegSwitch == 1)
1266       return (1.2 * function(u, 0.05, process) * (y - 0.05));
1267     else
1268       return ((6.8778 * pow(10., -8.) * exp(14.5 / u) -
1269                0.008805 * exp((2.5 - y + 0.98 * u * y) / u)) /
1270               (1.0204 - u));
1271   }
1272
1273   return 0.;
1274 }
1275
1276 double Martini::function(double u, double y, int process) {
1277   if (process == 1)
1278     return getRate_qqg(u, y);
1279   else if (process == 2)
1280     return getRate_qqgamma(u, y);
1281   else if (process == 3)
1282     return getRate_ggg(u, y);
1283   else if (process == 4)
1284     return getRate_gqq(u, y);
1285
1286   return 0.;
1287 }
1288
1289 bool Martini::isCoherent(Parton &pIn, int sibling, double T) {
1290   bool coherentNow = false;
1291   const weak_ptr<PartonShower> pShower = pIn.shower();
1292   auto shower = pShower.lock();
1293
1294   PartonShower::edge_iterator eIt, eEnd;
1295   for (eIt = shower->edges_begin(), eEnd = shower->edges_end(); eIt != eEnd;
1296        ++eIt) {
1297     if (eIt->target().outdeg() < 1) {
1298
1299       auto p = shower->GetParton(*eIt);
1300       if (p->plabel() == sibling) {
1301
1302         // synchronize the status of coherence with its sibling

1303         if (p->has_user_info<MARTINIUserInfo>()) {
1304           bool coherentSibling = p->user_info<MARTINIUserInfo>().coherent();
1305           if (!coherentSibling)
1306             return coherentSibling;
1307         } else {
1308           JSWARN << "MARTINIUserInfo not set!!"
1309                  << " plabel:" << p->plabel()
1310                  << " Make this parton in de-coherent state";
1311           return false;
1312         }
1313
1314         // check the sibling is out of energy loss modules

1315         bool activeCoherent = (fabs(p->x_in().t() - pIn.x_in().t()) < 0.015);
1316         if (!activeCoherent)
1317           return activeCoherent;
1318
1319         JSDEBUG << "  [Sibling]" << p->plabel() << " " << p->pid() << " "
1320                 << p->x_in().t() << " " << p->x_in().x();
1321
1322         // momentum at the spliting

1323         FourVector pInitVec = pIn.user_info<MARTINIUserInfo>().p_at_split();
1324         double pxInit = pInitVec.x();
1325         double pyInit = pInitVec.y();
1326         double pzInit = pInitVec.z();
1327         double pInit = pInitVec.t();
1328
1329         if (pInit < 1e-2)
1330           return false;
1331
1332         // spatial distance between two partons

1333         double xDelta = pIn.x_in().x() - p->x_in().x();
1334         double yDelta = pIn.x_in().x() - p->x_in().x();
1335         double zDelta = pIn.x_in().x() - p->x_in().x();
1336
1337         // cross product of pInit and rDelta

1338         double xCrox = yDelta * pzInit - zDelta * pyInit;
1339         double yCrox = zDelta * pxInit - xDelta * pzInit;
1340         double zCrox = xDelta * pyInit - yDelta * pxInit;
1341
1342         // transverse distance with respect to original parton

1343         double rPerp =
1344             sqrt(xCrox * xCrox + yCrox * yCrox + zCrox * zCrox) / pInit;
1345
1346         // momentum difference between two partons

1347         double pxDelta = pIn.px() - p->px();
1348         double pyDelta = pIn.py() - p->py();
1349         double pzDelta = pIn.pz() - p->pz();
1350
1351         // cross product of pInit and pDelta

1352         double pxCrox = pyDelta * pzInit - pzDelta * pyInit;
1353         double pyCrox = pzDelta * pxInit - pxDelta * pzInit;
1354         double pzCrox = pxDelta * pyInit - pyDelta * pxInit;
1355
1356         // transverse momentum with respect to original parton

1357         double pPerp =
1358             sqrt(pxCrox * pxCrox + pyCrox * pyCrox + pzCrox * pzCrox) / pInit;
1359
1360         double Crperp = 0.25 * pow(p->e() / T, 0.11);
1361         if (rPerp * pPerp < Crperp * hbarc)
1362           return true;
1363       }
1364     }
1365   }
1366
1367   return coherentNow;
1368 }
1369 RateElastic Martini::getRateElasTotal(double pRest, double T) {
1370   // compute the total transition rate in GeV, integrated over k, omega

1371   // and q and angles for fixed E=p and temperature T

1372   // using parametrization of numerically computed integral

1373   // then interpolate to get right value for used alpha_s

1374   // IMPORTANT: all computed values below are for a minimal omega of 0.05*T

1375   // so this is the cutoff to use in the calculation also!

1376   // also Nf=3 was used ... scales out though

1377
1378   RateElastic rate;
1379
1380   double u = pRest / T; // making arguments in log to be dimensionless

1381
1382   double alpha0 = 0.15;
1383   double deltaAlpha = 0.03;
1384   double iAlpha = floor((alpha_s - alpha0) / deltaAlpha + 0.001);
1385   double alphaFrac = (alpha_s - alpha0) / deltaAlpha - iAlpha;
1386   double rateLower, rateUpper;
1387   double c[2][6], d[2][6];
1388
1389   for (int i = 0; i < 2; i++)
1390     for (int j = 0; j < 6; j++) {
1391       c[i][j] = 0.;
1392       d[i][j] = 0.;
1393     }
1394
1395   if (alpha_s >= 0.15 && alpha_s < 0.18) {
1396     c[0][0] = 0.18172488396136807;
1397     c[1][0] = 0.224596478395945;
1398     c[0][1] = 0.6004740049060965;
1399     c[1][1] = 1.0874259848101948;
1400     c[0][2] = 0.36559627257898347;
1401     c[1][2] = 0.6436398538984057;
1402     c[0][3] = 0.10607576568373664;
1403     c[1][3] = 0.11585154613692052;
1404     c[0][4] = 0.004322466954618182;
1405     c[1][4] = -0.001719701730785056;
1406     c[0][5] = 0.04731599462749122;
1407     c[1][5] = 0.06734745496415469;
1408   } else if (alpha_s >= 0.18 && alpha_s < 0.21) {
1409     c[0][0] = 0.224596478395945;
1410     c[1][0] = 0.2686436092048326;
1411     c[0][1] = 1.0874259848101948;
1412     c[1][1] = 1.7286136256785387;
1413     c[0][2] = 0.6436398538984057;
1414     c[1][2] = 0.9826325498183079;
1415     c[0][3] = 0.11585154613692052;
1416     c[1][3] = 0.13136670133029682;
1417     c[0][4] = -0.001719701730785056;
1418     c[1][4] = -0.004876376882437649;
1419     c[0][5] = 0.06734745496415469;
1420     c[1][5] = 0.09140316977554151;
1421   } else if (alpha_s >= 0.21 && alpha_s < 0.24) {
1422     c[0][0] = 0.2686436092048326;
1423     c[1][0] = 0.3137234778163784;
1424     c[0][1] = 1.7286136256785387;
1425     c[1][1] = 2.445764079999846;
1426     c[0][2] = 0.9826325498183079;
1427     c[1][2] = 1.3083241146035964;
1428     c[0][3] = 0.13136670133029682;
1429     c[1][3] = 0.18341717903923757;
1430     c[0][4] = -0.004876376882437649;
1431     c[1][4] = 0.006098371807040589;
1432     c[0][5] = 0.09140316977554151;
1433     c[1][5] = 0.12054238276023879;
1434   } else if (alpha_s >= 0.24 && alpha_s < 0.27) {
1435     c[0][0] = 0.3137234778163784;
1436     c[1][0] = 0.3597255453974444;
1437     c[0][1] = 2.445764079999846;
1438     c[1][1] = 3.140669321831845;
1439     c[0][2] = 1.3083241146035964;
1440     c[1][2] = 1.535549334026633;
1441     c[0][3] = 0.18341717903923757;
1442     c[1][3] = 0.30505450230754705;
1443     c[0][4] = 0.006098371807040589;
1444     c[1][4] = 0.04285103618362223;
1445     c[0][5] = 0.12054238276023879;
1446     c[1][5] = 0.1558288379712527;
1447   } else if (alpha_s >= 0.27 && alpha_s < 0.3) {
1448     c[0][0] = 0.3597255453974444;
1449     c[1][0] = 0.40656130602563223;
1450     c[0][1] = 3.140669321831845;
1451     c[1][1] = 3.713430971987352;
1452     c[0][2] = 1.535549334026633;
1453     c[1][2] = 1.5818298058630476;
1454     c[0][3] = 0.30505450230754705;
1455     c[1][3] = 0.5269042544852683;
1456     c[0][4] = 0.04285103618362223;
1457     c[1][4] = 0.11594975218839362;
1458     c[0][5] = 0.1558288379712527;
1459     c[1][5] = 0.1982063104156748;
1460   } else if (alpha_s >= 0.3 && alpha_s < 0.33) {
1461     c[0][0] = 0.40656130602563223;
1462     c[1][0] = 0.45415805200862863;
1463     c[0][1] = 3.713430971987352;
1464     c[1][1] = 4.0758813206143785;
1465     c[0][2] = 1.5818298058630476;
1466     c[1][2] = 1.3775134184861555;
1467     c[0][3] = 0.5269042544852683;
1468     c[1][3] = 0.873527536823307;
1469     c[0][4] = 0.11594975218839362;
1470     c[1][4] = 0.23371456949506658;
1471     c[0][5] = 0.1982063104156748;
1472     c[1][5] = 0.24840524848507203;
1473   } else if (alpha_s >= 0.33 && alpha_s < 0.36) {
1474     c[0][0] = 0.45415805200862863;
1475     c[1][0] = 0.5024541413891354;
1476     c[0][1] = 4.0758813206143785;
1477     c[1][1] = 4.159425815179756;
1478     c[0][2] = 1.3775134184861555;
1479     c[1][2] = 0.8719749565879445;
1480     c[0][3] = 0.873527536823307;
1481     c[1][3] = 1.3606690530660879;
1482     c[0][4] = 0.23371456949506658;
1483     c[1][4] = 0.4010658149846402;
1484     c[0][5] = 0.24840524848507203;
1485     c[1][5] = 0.3067901992139913;
1486   } else if (alpha_s >= 0.36 && alpha_s < 0.39) {
1487     c[0][0] = 0.5024541413891354;
1488     c[1][0] = 0.5513999693402064;
1489     c[0][1] = 4.159425815179756;
1490     c[1][1] = 3.893153859527746;
1491     c[0][2] = 0.8719749565879445;
1492     c[1][2] = 0.009578762778659829;
1493     c[0][3] = 1.3606690530660879;
1494     c[1][3] = 2.0095157488463244;
1495     c[0][4] = 0.4010658149846402;
1496     c[1][4] = 0.6260756501912864;
1497     c[0][5] = 0.3067901992139913;
1498     c[1][5] = 0.37424991045026396;
1499   } else if (alpha_s >= 0.39 && alpha_s <= 0.42) {
1500     c[0][0] = 0.5513999693402064;
1501     c[1][0] = 0.600941593540798;
1502     c[0][1] = 3.893153859527746;
1503     c[1][1] = 3.293344337592684;
1504     c[0][2] = 0.009578762778659829;
1505     c[1][2] = -1.1764805445298645;
1506     c[0][3] = 2.0095157488463244;
1507     c[1][3] = 2.792180001243466;
1508     c[0][4] = 0.6260756501912864;
1509     c[1][4] = 0.8949534049225013;
1510     c[0][5] = 0.37424991045026396;
1511     c[1][5] = 0.44878529934031575;
1512   }
1513
1514   rateLower = T * (c[0][0] + c[0][1] / pow(u, 4.) - c[0][2] / pow(u, 3.) -
1515                    c[0][3] / pow(u, 2.) + c[0][4] / pow(u, 1.5) - c[0][5] / u);
1516   rateUpper = T * (c[1][0] + c[1][1] / pow(u, 4.) - c[1][2] / pow(u, 3.) -
1517                    c[1][3] / pow(u, 2.) + c[1][4] / pow(u, 1.5) - c[1][5] / u);
1518
1519   rate.qq = (1. - alphaFrac) * rateLower + alphaFrac * rateUpper;
1520   rate.qq *= nf / 3.; // adjust number of flavors

1521
1522   rate.gq = rate.qq * 9. / 4.;
1523
1524   if (alpha_s >= 0.15 && alpha_s < 0.18) {
1525     d[0][0] = 0.9364689080337059;
1526     d[1][0] = 1.1485486950080581;
1527     d[0][1] = 2.626076478553979;
1528     d[1][1] = 4.993647646894147;
1529     d[0][2] = 2.1171556605834274;
1530     d[1][2] = 3.7295251994302876;
1531     d[0][3] = 0.13123339226210134;
1532     d[1][3] = -0.0017620287506503757;
1533     d[0][4] = 0.02875811664147147;
1534     d[1][4] = 0.010598257485913224;
1535     d[0][5] = 0.27736469898722244;
1536     d[1][5] = 0.3949856219367327;
1537   } else if (alpha_s >= 0.18 && alpha_s < 0.21) {
1538     d[0][0] = 1.1485486950080581;
1539     d[1][0] = 1.3645568637616001;
1540     d[0][1] = 4.993647646894147;
1541     d[1][1] = 8.174225869366722;
1542     d[0][2] = 3.7295251994302876;
1543     d[1][2] = 5.732101892684938;
1544     d[0][3] = -0.001762028750650376;
1545     d[1][3] = -0.1416811579957863;
1546     d[0][4] = 0.010598257485913224;
1547     d[1][4] = 0.011703596451947428;
1548     d[0][5] = 0.3949856219367327;
1549     d[1][5] = 0.5354757997870718;
1550   } else if (alpha_s >= 0.21 && alpha_s < 0.24) {
1551     d[0][0] = 1.3645568637616001;
1552     d[1][0] = 1.5839378568555678;
1553     d[0][1] = 8.174225869366722;
1554     d[1][1] = 11.785897000063443;
1555     d[0][2] = 5.732101892684938;
1556     d[1][2] = 7.758388282689373;
1557     d[0][3] = -0.1416811579957863;
1558     d[1][3] = -0.13163385415183002;
1559     d[0][4] = 0.011703596451947428;
1560     d[1][4] = 0.09016386041913003;
1561     d[0][5] = 0.5354757997870718;
1562     d[1][5] = 0.7042577279136836;
1563   } else if (alpha_s >= 0.24 && alpha_s < 0.27) {
1564     d[0][0] = 1.5839378568555678;
1565     d[1][0] = 1.8062676019060235;
1566     d[0][1] = 11.785897000063443;
1567     d[1][1] = 15.344112642069764;
1568     d[0][2] = 7.758388282689373;
1569     d[1][2] = 9.384190917330093;
1570     d[0][3] = -0.13163385415183002;
1571     d[1][3] = 0.19709400976261568;
1572     d[0][4] = 0.09016386041913003;
1573     d[1][4] = 0.30577623140224813;
1574     d[0][5] = 0.7042577279136836;
1575     d[1][5] = 0.9066501895009754;
1576   } else if (alpha_s >= 0.27 && alpha_s < 0.3) {
1577     d[0][0] = 1.8062676019060235;
1578     d[1][0] = 2.0312125903238236;
1579     d[0][1] = 15.344112642069764;
1580     d[1][1] = 18.36844006721506;
1581     d[0][2] = 9.384190917330093;
1582     d[1][2] = 10.209988454804193;
1583     d[0][3] = 0.19709400976261568;
1584     d[1][3] = 0.9957025988944573;
1585     d[0][4] = 0.30577623140224813;
1586     d[1][4] = 0.7109302867706849;
1587     d[0][5] = 0.9066501895009754;
1588     d[1][5] = 1.1472148515742653;
1589   } else if (alpha_s >= 0.3 && alpha_s < 0.33) {
1590     d[0][0] = 2.0312125903238236;
1591     d[1][0] = 2.258502734110078;
1592     d[0][1] = 18.36844006721506;
1593     d[1][1] = 20.43444928479894;
1594     d[0][2] = 10.209988454804193;
1595     d[1][2] = 9.896928897847518;
1596     d[0][3] = 0.9957025988944573;
1597     d[1][3] = 2.3867073785159003;
1598     d[0][4] = 0.7109302867706849;
1599     d[1][4] = 1.3473328178504662;
1600     d[0][5] = 1.1472148515742653;
1601     d[1][5] = 1.429497460496924;
1602   } else if (alpha_s >= 0.33 && alpha_s < 0.36) {
1603     d[0][0] = 2.258502734110078;
1604     d[1][0] = 2.4879110920956653;
1605     d[0][1] = 20.43444928479894;
1606     d[1][1] = 21.220550462966102;
1607     d[0][2] = 9.896928897847518;
1608     d[1][2] = 8.20639681844989;
1609     d[0][3] = 2.3867073785159003;
1610     d[1][3] = 4.445222616370339;
1611     d[0][4] = 1.3473328178504662;
1612     d[1][4] = 2.2381176005506016;
1613     d[0][5] = 1.429497460496924;
1614     d[1][5] = 1.7550164762706189;
1615   } else if (alpha_s >= 0.36 && alpha_s < 0.39) {
1616     d[0][0] = 2.4879110920956653;
1617     d[1][0] = 2.7192501243929903;
1618     d[0][1] = 21.220550462966102;
1619     d[1][1] = 20.470583876561985;
1620     d[0][2] = 8.20639681844989;
1621     d[1][2] = 4.954737209403953;
1622     d[0][3] = 4.445222616370339;
1623     d[1][3] = 7.227667929705693;
1624     d[0][4] = 2.2381176005506016;
1625     d[1][4] = 3.401378906197122;
1626     d[0][5] = 1.7550164762706189;
1627     d[1][5] = 2.1251383942923474;
1628   } else if (alpha_s >= 0.39 && alpha_s <= 0.42) {
1629     d[0][0] = 2.7192501243929903;
1630     d[1][0] = 2.9523522354248817;
1631     d[0][1] = 20.470583876561985;
1632     d[1][1] = 18.027772799078463;
1633     d[0][2] = 4.954737209403953;
1634     d[1][2] = 0.050298242947981846;
1635     d[0][3] = 7.227667929705693;
1636     d[1][3] = 10.747352232336384;
1637     d[0][4] = 3.401378906197122;
1638     d[1][4] = 4.8378133911595285;
1639     d[0][5] = 2.1251383942923474;
1640     d[1][5] = 2.5391647730624003;
1641   }
1642
1643   rateLower = T * (d[0][0] + d[0][1] / pow(u, 4.) - d[0][2] / pow(u, 3.) -
1644                    d[0][3] / pow(u, 2.) + d[0][4] / pow(u, 1.5) - d[0][5] / u);
1645   rateUpper = T * (d[1][0] + d[1][1] / pow(u, 4.) - d[1][2] / pow(u, 3.) -
1646                    d[1][3] / pow(u, 2.) + d[1][4] / pow(u, 1.5) - d[1][5] / u);
1647
1648   rate.qg = (1. - alphaFrac) * rateLower + alphaFrac * rateUpper;
1649   rate.qg /= 3.; // historic reasons

1650
1651   rate.gg = rate.qg * 9. / 4.;
1652
1653   return rate;
1654 }
1655
1656 RateElastic Martini::getRateElasPos(double u, double T) {
1657   RateElastic rate;
1658
1659   double alpha0 = 0.15;
1660   double deltaAlpha = 0.03;
1661   double iAlpha = floor((alpha_s - alpha0) / deltaAlpha + 0.001);
1662   double alphaFrac = (alpha_s - alpha0) / deltaAlpha - iAlpha;
1663   double rateLower, rateUpper;
1664   double c[2][6], d[2][6];
1665
1666   for (int i = 0; i < 2; i++)
1667     for (int j = 0; j < 6; j++) {
1668       c[i][j] = 0.;
1669       d[i][j] = 0.;
1670     }
1671
1672   if (alpha_s >= 0.15 && alpha_s < 0.18) {
1673     c[0][0] = 0.12199410313320332;
1674     c[1][0] = 0.15243607717720586;
1675     c[0][1] = 0.23732051765097376;
1676     c[1][1] = 0.5403120875137825;
1677     c[0][2] = -0.03285419708803458;
1678     c[1][2] = 0.06440920730334501;
1679     c[0][3] = 0.2255419254079952;
1680     c[1][3] = 0.2881594349535524;
1681     c[0][4] = 0.03991522899907729;
1682     c[1][4] = 0.04948438583750772;
1683     c[0][5] = 0.05022641428394594;
1684     c[1][5] = 0.07152523367501308;
1685   } else if (alpha_s >= 0.18 && alpha_s < 0.21) {
1686     c[0][0] = 0.15243607717720586;
1687     c[1][0] = 0.15243607717720586;
1688     c[0][1] = 0.5403120875137825;
1689     c[1][1] = 0.5403120875137825;
1690     c[0][2] = 0.06440920730334501;
1691     c[1][2] = 0.06440920730334501;
1692     c[0][3] = 0.2881594349535524;
1693     c[1][3] = 0.2881594349535524;
1694     c[0][4] = 0.04948438583750772;
1695     c[1][4] = 0.04948438583750772;
1696     c[0][5] = 0.07152523367501308;
1697     c[1][5] = 0.07152523367501308;
1698   } else if (alpha_s >= 0.21 && alpha_s < 0.24) {
1699     c[0][0] = 0.15243607717720586;
1700     c[1][0] = 0.21661000995329158;
1701     c[0][1] = 0.5403120875137825;
1702     c[1][1] = 1.4087570376612657;
1703     c[0][2] = 0.06440920730334501;
1704     c[1][2] = 0.2713885880193171;
1705     c[0][3] = 0.2881594349535524;
1706     c[1][3] = 0.48681971936565244;
1707     c[0][4] = 0.04948438583750772;
1708     c[1][4] = 0.09567346780679847;
1709     c[0][5] = 0.07152523367501308;
1710     c[1][5] = 0.12780677622585393;
1711   } else if (alpha_s >= 0.24 && alpha_s < 0.27) {
1712     c[0][0] = 0.21661000995329158;
1713     c[1][0] = 0.2501007467879627;
1714     c[0][1] = 1.4087570376612657;
1715     c[1][1] = 1.8034683081244214;
1716     c[0][2] = 0.2713885880193171;
1717     c[1][2] = 0.228092470920281;
1718     c[0][3] = 0.48681971936565244;
1719     c[1][3] = 0.6841577896561725;
1720     c[0][4] = 0.09567346780679847;
1721     c[1][4] = 0.15430793601338547;
1722     c[0][5] = 0.12780677622585393;
1723     c[1][5] = 0.1648297331159989;
1724   } else if (alpha_s >= 0.27 && alpha_s < 0.3) {
1725     c[0][0] = 0.2501007467879627;
1726     c[1][0] = 0.28440720063047276;
1727     c[0][1] = 1.8034683081244214;
1728     c[1][1] = 2.0448244620634055;
1729     c[0][2] = 0.228092470920281;
1730     c[1][2] = -0.018574547528236382;
1731     c[0][3] = 0.6841577896561725;
1732     c[1][3] = 0.9863974758613413;
1733     c[0][4] = 0.15430793601338547;
1734     c[1][4] = 0.2503738253300167;
1735     c[0][5] = 0.1648297331159989;
1736     c[1][5] = 0.2090067594645225;
1737   } else if (alpha_s >= 0.3 && alpha_s < 0.33) {
1738     c[0][0] = 0.28440720063047276;
1739     c[1][0] = 0.31945943548344036;
1740     c[0][1] = 2.0448244620634055;
1741     c[1][1] = 2.0482495934952256;
1742     c[0][2] = -0.018574547528236382;
1743     c[1][2] = -0.5350999123662686;
1744     c[0][3] = 0.9863974758613413;
1745     c[1][3] = 1.4169725257394696;
1746     c[0][4] = 0.2503738253300167;
1747     c[1][4] = 0.3918202096574105;
1748     c[0][5] = 0.2090067594645225;
1749     c[1][5] = 0.26103455441873036;
1750   } else if (alpha_s >= 0.33 && alpha_s < 0.36) {
1751     c[0][0] = 0.31945943548344036;
1752     c[1][0] = 0.35519799231686516;
1753     c[0][1] = 2.0482495934952256;
1754     c[1][1] = 1.7485135425544152;
1755     c[0][2] = -0.5350999123662686;
1756     c[1][2] = -1.3692232011881413;
1757     c[0][3] = 1.4169725257394696;
1758     c[1][3] = 1.9906086576701993;
1759     c[0][4] = 0.3918202096574105;
1760     c[1][4] = 0.5832315715098879;
1761     c[0][5] = 0.26103455441873036;
1762     c[1][5] = 0.32124694953933486;
1763   } else if (alpha_s >= 0.36 && alpha_s < 0.39) {
1764     c[0][0] = 0.35519799231686516;
1765     c[1][0] = 0.39157507493019383;
1766     c[0][1] = 1.7485135425544152;
1767     c[1][1] = 1.0778995684787331;
1768     c[0][2] = -1.3692232011881413;
1769     c[1][2] = -2.5738838613236457;
1770     c[0][3] = 1.9906086576701993;
1771     c[1][3] = 2.727543221296746;
1772     c[0][4] = 0.5832315715098879;
1773     c[1][4] = 0.8323699786704292;
1774     c[0][5] = 0.32124694953933486;
1775     c[1][5] = 0.3905055907877247;
1776   } else if (alpha_s >= 0.39 && alpha_s <= 0.42) {
1777     c[0][0] = 0.39157507493019383;
1778     c[1][0] = 0.4285382777192131;
1779     c[0][1] = 1.0778995684787331;
1780     c[1][1] = 0.05505396151716547;
1781     c[0][2] = -2.5738838613236457;
1782     c[1][2] = -4.113979132685303;
1783     c[0][3] = 2.727543221296746;
1784     c[1][3] = 3.5992808060371506;
1785     c[0][4] = 0.8323699786704292;
1786     c[1][4] = 1.1252568207814462;
1787     c[0][5] = 0.3905055907877247;
1788     c[1][5] = 0.4667953957378259;
1789   }
1790
1791   rateLower = T * (c[0][0] + c[0][1] / pow(u, 4.) - c[0][2] / pow(u, 3.) -
1792                    c[0][3] / pow(u, 2.) + c[0][4] / pow(u, 1.5) - c[0][5] / u);
1793   rateUpper = T * (c[1][0] + c[1][1] / pow(u, 4.) - c[1][2] / pow(u, 3.) -
1794                    c[1][3] / pow(u, 2.) + c[1][4] / pow(u, 1.5) - c[1][5] / u);
1795
1796   rate.qq = (1. - alphaFrac) * rateLower + alphaFrac * rateUpper;
1797   rate.qq *= nf / 3.; // adjust number of flavors

1798
1799   rate.gq = rate.qq * 9. / 4.;
1800
1801   if (alpha_s >= 0.15 && alpha_s < 0.18) {
1802     d[0][0] = 0.6197775378922895;
1803     d[1][0] = 0.7680959463632293;
1804     d[0][1] = 1.5268694134079064;
1805     d[1][1] = 3.282164035377037;
1806     d[0][2] = 0.6939337312845367;
1807     d[1][2] = 1.6359849897319092;
1808     d[0][3] = 0.5967602676773388;
1809     d[1][3] = 0.6770046238563808;
1810     d[0][4] = 0.17320784052297564;
1811     d[1][4] = 0.22074166337990309;
1812     d[0][5] = 0.28964614117694565;
1813     d[1][5] = 0.4128184793199476;
1814   } else if (alpha_s >= 0.18 && alpha_s < 0.21) {
1815     d[0][0] = 0.7680959463632293;
1816     d[1][0] = 0.9206225398305536;
1817     d[0][1] = 3.282164035377037;
1818     d[1][1] = 5.690562370150853;
1819     d[0][2] = 1.6359849897319092;
1820     d[1][2] = 2.8341906487774318;
1821     d[0][3] = 0.6770046238563808;
1822     d[1][3] = 0.7900156706763937;
1823     d[0][4] = 0.22074166337990309;
1824     d[1][4] = 0.2995126102416747;
1825     d[0][5] = 0.4128184793199476;
1826     d[1][5] = 0.5598645426609049;
1827   } else if (alpha_s >= 0.21 && alpha_s < 0.24) {
1828     d[0][0] = 0.9206225398305536;
1829     d[1][0] = 1.0767954081327265;
1830     d[0][1] = 5.690562370150853;
1831     d[1][1] = 8.378841394880034;
1832     d[0][2] = 2.8341906487774318;
1833     d[1][2] = 3.9338968631891396;
1834     d[0][3] = 0.7900156706763937;
1835     d[1][3] = 1.0874771229885156;
1836     d[0][4] = 0.2995126102416747;
1837     d[1][4] = 0.46570985770548107;
1838     d[0][5] = 0.5598645426609049;
1839     d[1][5] = 0.7360069767362173;
1840   } else if (alpha_s >= 0.24 && alpha_s < 0.27) {
1841     d[0][0] = 1.0767954081327265;
1842     d[1][0] = 1.2361819653856791;
1843     d[0][1] = 8.378841394880034;
1844     d[1][1] = 10.877148035367144;
1845     d[0][2] = 3.9338968631891396;
1846     d[1][2] = 4.526191560392149;
1847     d[0][3] = 1.0874771229885156;
1848     d[1][3] = 1.731930015138816;
1849     d[0][4] = 0.46570985770548107;
1850     d[1][4] = 0.7769917594310469;
1851     d[0][5] = 0.7360069767362173;
1852     d[1][5] = 0.9463662091275489;
1853   } else if (alpha_s >= 0.27 && alpha_s < 0.3) {
1854     d[0][0] = 1.2361819653856791;
1855     d[1][0] = 1.3984393292278847;
1856     d[0][1] = 10.877148035367144;
1857     d[1][1] = 12.72181515837248;
1858     d[0][2] = 4.526191560392149;
1859     d[1][2] = 4.227297031355039;
1860     d[0][3] = 1.731930015138816;
1861     d[1][3] = 2.868526983329731;
1862     d[0][4] = 0.7769917594310469;
1863     d[1][4] = 1.2836917844304823;
1864     d[0][5] = 0.9463662091275489;
1865     d[1][5] = 1.1953148369630755;
1866   } else if (alpha_s >= 0.3 && alpha_s < 0.33) {
1867     d[0][0] = 1.3984393292278847;
1868     d[1][0] = 1.5632880021613935;
1869     d[0][1] = 12.72181515837248;
1870     d[1][1] = 13.502896915302873;
1871     d[0][2] = 4.227297031355039;
1872     d[1][2] = 2.7113406243010467;
1873     d[0][3] = 2.868526983329731;
1874     d[1][3] = 4.615035662049938;
1875     d[0][4] = 1.2836917844304823;
1876     d[1][4] = 2.0259357821768784;
1877     d[0][5] = 1.1953148369630755;
1878     d[1][5] = 1.486253368704046;
1879   } else if (alpha_s >= 0.33 && alpha_s < 0.36) {
1880     d[0][0] = 1.5632880021613935;
1881     d[1][0] = 1.730492163581557;
1882     d[0][1] = 13.502896915302873;
1883     d[1][1] = 12.913294655478987;
1884     d[0][2] = 2.7113406243010467;
1885     d[1][2] = -0.2477159937428581;
1886     d[0][3] = 4.615035662049938;
1887     d[1][3] = 7.042004003229154;
1888     d[0][4] = 2.0259357821768784;
1889     d[1][4] = 3.0253452576771465;
1890     d[0][5] = 1.486253368704046;
1891     d[1][5] = 1.8205651561017433;
1892   } else if (alpha_s >= 0.36 && alpha_s < 0.39) {
1893     d[0][0] = 1.730492163581557;
1894     d[1][0] = 1.8998560359992867;
1895     d[0][1] = 12.913294655478987;
1896     d[1][1] = 10.708892844334745;
1897     d[0][2] = -0.2477159937428581;
1898     d[1][2] = -4.823210983922782;
1899     d[0][3] = 7.042004003229154;
1900     d[1][3] = 10.202109059054063;
1901     d[0][4] = 3.0253452576771465;
1902     d[1][4] = 4.298747764427364;
1903     d[0][5] = 1.8205651561017433;
1904     d[1][5] = 2.199497022778097;
1905   } else if (alpha_s >= 0.39 && alpha_s <= 0.42) {
1906     d[0][0] = 1.8998560359992867;
1907     d[1][0] = 2.071204284004704;
1908     d[0][1] = 10.708892844334745;
1909     d[1][1] = 6.741738604119316;
1910     d[0][2] = -4.823210983922782;
1911     d[1][2] = -11.099716230158746;
1912     d[0][3] = 10.202109059054063;
1913     d[1][3] = 14.106488110189458;
1914     d[0][4] = 4.298747764427364;
1915     d[1][4] = 5.846203546614067;
1916     d[0][5] = 2.199497022778097;
1917     d[1][5] = 2.62230136903594;
1918   }
1919
1920   rateLower = T * (d[0][0] + d[0][1] / pow(u, 4.) - d[0][2] / pow(u, 3.) -
1921                    d[0][3] / pow(u, 2.) + d[0][4] / pow(u, 1.5) - d[0][5] / u);
1922   rateUpper = T * (d[1][0] + d[1][1] / pow(u, 4.) - d[1][2] / pow(u, 3.) -
1923                    d[1][3] / pow(u, 2.) + d[1][4] / pow(u, 1.5) - d[1][5] / u);
1924
1925   rate.qg = (1. - alphaFrac) * rateLower + alphaFrac * rateUpper;
1926   rate.qg /= 3.; // historic reasons

1927
1928   rate.gg = rate.qg * 9. / 4.;
1929
1930   return rate;
1931 }
1932
1933 RateElastic Martini::getRateElasNeg(double u, double T) {
1934   RateElastic rate;
1935
1936   double alpha0 = 0.15;
1937   double deltaAlpha = 0.03;
1938   double iAlpha = floor((alpha_s - alpha0) / deltaAlpha + 0.001);
1939   double alphaFrac = (alpha_s - alpha0) / deltaAlpha - iAlpha;
1940   double rateLower, rateUpper;
1941   double c[2][6], d[2][6];
1942
1943   for (int i = 0; i < 2; i++)
1944     for (int j = 0; j < 6; j++) {
1945       c[i][j] = 0.;
1946       d[i][j] = 0.;
1947     }
1948
1949   if (alpha_s >= 0.15 && alpha_s < 0.18) {
1950     c[0][0] = 0.059730780828164666;
1951     c[1][0] = 0.07216040121873951;
1952     c[0][1] = 0.3631534872548789;
1953     c[1][1] = 0.5471138972952214;
1954     c[0][2] = 0.39845046966687;
1955     c[1][2] = 0.5792306465939813;
1956     c[0][3] = 0.11946615972422633;
1957     c[1][3] = 0.1723078888161528;
1958     c[0][4] = 0.03559276204445307;
1959     c[1][4] = 0.05120408756812135;
1960     c[0][5] = 0.00291041965645416;
1961     c[1][5] = 0.0041777787108426695;
1962   } else if (alpha_s >= 0.18 && alpha_s < 0.21) {
1963     c[0][0] = 0.07216040121873951;
1964     c[1][0] = 0.0846236909779996;
1965     c[0][1] = 0.5471138972952214;
1966     c[1][1] = 0.7725791286875564;
1967     c[0][2] = 0.5792306465939813;
1968     c[1][2] = 0.7931123494736929;
1969     c[0][3] = 0.1723078888161528;
1970     c[1][3] = 0.23406373724706608;
1971     c[0][4] = 0.05120408756812135;
1972     c[1][4] = 0.06935459958589639;
1973     c[0][5] = 0.0041777787108426695;
1974     c[1][5] = 0.005644055718614478;
1975   } else if (alpha_s >= 0.21 && alpha_s < 0.24) {
1976     c[0][0] = 0.0846236909779996;
1977     c[1][0] = 0.09711346786308672;
1978     c[0][1] = 0.7725791286875564;
1979     c[1][1] = 1.0370070423372528;
1980     c[0][2] = 0.7931123494736929;
1981     c[1][2] = 1.036935526583188;
1982     c[0][3] = 0.23406373724706608;
1983     c[1][3] = 0.3034025403259155;
1984     c[0][4] = 0.06935459958589639;
1985     c[1][4] = 0.08957509599955729;
1986     c[0][5] = 0.005644055718614478;
1987     c[1][5] = 0.007264393465593115;
1988   } else if (alpha_s >= 0.24 && alpha_s < 0.27) {
1989     c[0][0] = 0.09711346786308672;
1990     c[1][0] = 0.10962479860948156;
1991     c[0][1] = 1.0370070423372528;
1992     c[1][1] = 1.3372010137066646;
1993     c[0][2] = 1.036935526583188;
1994     c[1][2] = 1.307456863105879;
1995     c[0][3] = 0.3034025403259155;
1996     c[1][3] = 0.37910328734850873;
1997     c[0][4] = 0.08957509599955729;
1998     c[1][4] = 0.111456899829735;
1999     c[0][5] = 0.007264393465593115;
2000     c[1][5] = 0.009000895144744121;
2001   } else if (alpha_s >= 0.27 && alpha_s < 0.3) {
2002     c[0][0] = 0.10962479860948156;
2003     c[1][0] = 0.1221541053951596;
2004     c[0][1] = 1.3372010137066646;
2005     c[1][1] = 1.6686065099273535;
2006     c[0][2] = 1.307456863105879;
2007     c[1][2] = 1.600404353394210;
2008     c[0][3] = 0.37910328734850873;
2009     c[1][3] = 0.4594932213772782;
2010     c[0][4] = 0.111456899829735;
2011     c[1][4] = 0.13442407314203592;
2012     c[0][5] = 0.009000895144744121;
2013     c[1][5] = 0.010800449048880756;
2014   } else if (alpha_s >= 0.3 && alpha_s < 0.33) {
2015     c[0][0] = 0.1221541053951596;
2016     c[1][0] = 0.13469861652518803;
2017     c[0][1] = 1.6686065099273535;
2018     c[1][1] = 2.0276317271182074;
2019     c[0][2] = 1.600404353394210;
2020     c[1][2] = 1.912613330851788;
2021     c[0][3] = 0.4594932213772782;
2022     c[1][3] = 0.5434449889160747;
2023     c[0][4] = 0.13442407314203592;
2024     c[1][4] = 0.15810564016236883;
2025     c[0][5] = 0.010800449048880756;
2026     c[1][5] = 0.012629305933671075;
2027   } else if (alpha_s >= 0.33 && alpha_s < 0.36) {
2028     c[0][0] = 0.13469861652518803;
2029     c[1][0] = 0.14725614907227047;
2030     c[0][1] = 2.0276317271182074;
2031     c[1][1] = 2.4109122726272654;
2032     c[0][2] = 1.912613330851788;
2033     c[1][2] = 2.241198157777867;
2034     c[0][3] = 0.5434449889160747;
2035     c[1][3] = 0.6299396046048817;
2036     c[0][4] = 0.15810564016236883;
2037     c[1][4] = 0.18216575652552597;
2038     c[0][5] = 0.012629305933671075;
2039     c[1][5] = 0.014456750325370632;
2040   } else if (alpha_s >= 0.36 && alpha_s < 0.39) {
2041     c[0][0] = 0.14725614907227047;
2042     c[1][0] = 0.15982489441001274;
2043     c[0][1] = 2.4109122726272654;
2044     c[1][1] = 2.815254291049982;
2045     c[0][2] = 2.241198157777867;
2046     c[1][2] = 2.583462624103292;
2047     c[0][3] = 0.6299396046048817;
2048     c[1][3] = 0.7180274724508857;
2049     c[0][4] = 0.18216575652552597;
2050     c[1][4] = 0.20629432847931367;
2051     c[0][5] = 0.014456750325370632;
2052     c[1][5] = 0.01625568033747704;
2053   } else if (alpha_s >= 0.39 && alpha_s <= 0.42) {
2054     c[0][0] = 0.15982489441001274;
2055     c[1][0] = 0.17240331582158486;
2056     c[0][1] = 2.815254291049982;
2057     c[1][1] = 3.238290376079149;
2058     c[0][2] = 2.583462624103292;
2059     c[1][2] = 2.9374985881586273;
2060     c[0][3] = 0.7180274724508857;
2061     c[1][3] = 0.8071008047950518;
2062     c[0][4] = 0.20629432847931367;
2063     c[1][4] = 0.23030341585944009;
2064     c[0][5] = 0.01625568033747704;
2065     c[1][5] = 0.018010096397556033;
2066   }
2067
2068   rateLower = T * (c[0][0] + c[0][1] / pow(u, 4.) - c[0][2] / pow(u, 3.) -
2069                    c[0][3] / pow(u, 2.) + c[0][4] / pow(u, 1.5) - c[0][5] / u);
2070   rateUpper = T * (c[1][0] + c[1][1] / pow(u, 4.) - c[1][2] / pow(u, 3.) -
2071                    c[1][3] / pow(u, 2.) + c[1][4] / pow(u, 1.5) - c[1][5] / u);
2072
2073   rate.qq = (1. - alphaFrac) * rateLower + alphaFrac * rateUpper;
2074   rate.qq *= nf / 3.; // adjust number of flavors

2075
2076   rate.gq = rate.qq * 9. / 4.;
2077
2078   if (alpha_s >= 0.15 && alpha_s < 0.18) {
2079     d[0][0] = 0.3166913701414167;
2080     d[1][0] = 0.3804527486448292;
2081     d[0][1] = 1.0992070651449564;
2082     d[1][1] = 1.7114836115114735;
2083     d[0][2] = 1.4232219292986843;
2084     d[1][2] = 2.093540209692791;
2085     d[0][3] = 0.4655268754156213;
2086     d[1][3] = 0.6787666526042345;
2087     d[0][4] = 0.1444497238817506;
2088     d[1][4] = 0.21014340589291994;
2089     d[0][5] = 0.012281442189758532;
2090     d[1][5] = 0.017832857383112792;
2091   } else if (alpha_s >= 0.18 && alpha_s < 0.21) {
2092     d[0][0] = 0.3804527486448292;
2093     d[1][0] = 0.44393432393104637;
2094     d[0][1] = 1.7114836115114735;
2095     d[1][1] = 2.483663499207573;
2096     d[0][2] = 2.093540209692791;
2097     d[1][2] = 2.8979112438999044;
2098     d[0][3] = 0.6787666526042345;
2099     d[1][3] = 0.9316968286688833;
2100     d[0][4] = 0.21014340589291994;
2101     d[1][4] = 0.28780901378857465;
2102     d[0][5] = 0.017832857383112792;
2103     d[1][5] = 0.02438874287373154;
2104   } else if (alpha_s >= 0.21 && alpha_s < 0.24) {
2105     d[0][0] = 0.44393432393104637;
2106     d[1][0] = 0.5071424487228405;
2107     d[0][1] = 2.483663499207573;
2108     d[1][1] = 3.4070556051784515;
2109     d[0][2] = 2.8979112438999044;
2110     d[1][2] = 3.824491419496227;
2111     d[0][3] = 0.9316968286688833;
2112     d[1][3] = 1.2191109771387096;
2113     d[0][4] = 0.28780901378857465;
2114     d[1][4] = 0.3755459972857442;
2115     d[0][5] = 0.02438874287373154;
2116     d[1][5] = 0.03174924882247299;
2117   } else if (alpha_s >= 0.24 && alpha_s < 0.27) {
2118     d[0][0] = 0.5071424487228405;
2119     d[1][0] = 0.5700856365203443;
2120     d[0][1] = 3.4070556051784515;
2121     d[1][1] = 4.466964606692036;
2122     d[0][2] = 3.824491419496227;
2123     d[1][2] = 4.857999356928031;
2124     d[0][3] = 1.2191109771387096;
2125     d[1][3] = 1.5348360053714125;
2126     d[0][4] = 0.3755459972857442;
2127     d[1][4] = 0.471215528026891;
2128     d[0][5] = 0.03174924882247299;
2129     d[1][5] = 0.03971601962636114;
2130   } else if (alpha_s >= 0.27 && alpha_s < 0.3) {
2131     d[0][0] = 0.5700856365203443;
2132     d[1][0] = 0.6327732610959403;
2133     d[0][1] = 4.466964606692036;
2134     d[1][1] = 5.646624908846933;
2135     d[0][2] = 4.857999356928031;
2136     d[1][2] = 5.982691423451806;
2137     d[0][3] = 1.5348360053714125;
2138     d[1][3] = 1.8728243844356356;
2139     d[0][4] = 0.471215528026891;
2140     d[1][4] = 0.572761497659723;
2141     d[0][5] = 0.03971601962636114;
2142     d[1][5] = 0.04809998538877525;
2143   } else if (alpha_s >= 0.3 && alpha_s < 0.33) {
2144     d[0][0] = 0.6327732610959403;
2145     d[1][0] = 0.6952147319486842;
2146     d[0][1] = 5.646624908846933;
2147     d[1][1] = 6.931552369487635;
2148     d[0][2] = 5.982691423451806;
2149     d[1][2] = 7.185588273540373;
2150     d[0][3] = 1.8728243844356356;
2151     d[1][3] = 2.228328283532209;
2152     d[0][4] = 0.572761497659723;
2153     d[1][4] = 0.6786029643259804;
2154     d[0][5] = 0.04809998538877525;
2155     d[1][5] = 0.056755908207122875;
2156   } else if (alpha_s >= 0.33 && alpha_s < 0.36) {
2157     d[0][0] = 0.6952147319486842;
2158     d[1][0] = 0.7574189285141091;
2159     d[0][1] = 6.931552369487635;
2160     d[1][1] = 8.307255807497631;
2161     d[0][2] = 7.185588273540373;
2162     d[1][2] = 8.454112812202247;
2163     d[0][3] = 2.228328283532209;
2164     d[1][3] = 2.596781386863294;
2165     d[0][4] = 0.6786029643259804;
2166     d[1][4] = 0.7872276571283385;
2167     d[0][5] = 0.056755908207122875;
2168     d[1][5] = 0.06554867983133447;
2169   } else if (alpha_s >= 0.36 && alpha_s < 0.39) {
2170     d[0][0] = 0.7574189285141091;
2171     d[1][0] = 0.8193940883937045;
2172     d[0][1] = 8.307255807497631;
2173     d[1][1] = 9.761691032241623;
2174     d[0][2] = 8.454112812202247;
2175     d[1][2] = 9.777948193339808;
2176     d[0][3] = 2.596781386863294;
2177     d[1][3] = 2.9744411293541457;
2178     d[0][4] = 0.7872276571283385;
2179     d[1][4] = 0.8973688582323887;
2180     d[0][5] = 0.06554867983133447;
2181     d[1][5] = 0.07435862848596686;
2182   } else if (alpha_s >= 0.39 && alpha_s <= 0.42) {
2183     d[0][0] = 0.8193940883937045;
2184     d[1][0] = 0.8811479514201789;
2185     d[0][1] = 9.761691032241623;
2186     d[1][1] = 11.286034194965852;
2187     d[0][2] = 9.777948193339808;
2188     d[1][2] = 11.15001447311135;
2189     d[0][3] = 2.9744411293541457;
2190     d[1][3] = 3.3591358778545803;
2191     d[0][4] = 0.8973688582323887;
2192     d[1][4] = 1.0083901554550654;
2193     d[0][5] = 0.07435862848596686;
2194     d[1][5] = 0.08313659597360733;
2195   }
2196
2197   rateLower = T * (d[0][0] + d[0][1] / pow(u, 4.) - d[0][2] / pow(u, 3.) -
2198                    d[0][3] / pow(u, 2.) + d[0][4] / pow(u, 1.5) - d[0][5] / u);
2199   rateUpper = T * (d[1][0] + d[1][1] / pow(u, 4.) - d[1][2] / pow(u, 3.) -
2200                    d[1][3] / pow(u, 2.) + d[1][4] / pow(u, 1.5) - d[1][5] / u);
2201
2202   rate.qg = (1. - alphaFrac) * rateLower + alphaFrac * rateUpper;
2203   rate.qg /= 3.; // historic reasons

2204
2205   rate.gg = rate.qg * 9. / 4.;
2206
2207   return rate;
2208 }
2209
2210 RateConversion Martini::getRateConv(double pRest, double T) {
2211   RateConversion rate;
2212
2213   rate.qg = 4. / 3. * 2. * M_PI * alpha_s * alpha_s * T * T / (3. * pRest) *
2214             (0.5 * log(pRest * T / ((1. / 6.) * pow(g * T, 2.))) - 0.36149);
2215   rate.gq = nf * 3. / 8. * 4. / 3. * 2. * M_PI * alpha_s * alpha_s * T * T /
2216             (3. * pRest) *
2217             (0.5 * log(pRest * T / ((1. / 6.) * pow(g * T, 2.))) - 0.36149);
2218   rate.qgamma = 2. * M_PI * alpha_s * alpha_em * T * T / (3. * pRest) *
2219                 (0.5 * log(pRest * T / ((1. / 6.) * pow(g * T, 2.))) - 0.36149);
2220
2221   return rate;
2222 }
2223
2224 double Martini::getEnergyTransfer(double pRest, double T, int process) {
2225   double u = pRest / T; // making arguments in log to be dimensionless

2226
2227   double omega = 0.; // momentum of radiated gluon (dimentionless)

2228   double y;          // omega candidate

2229   double x;          // random number, uniform on [0,1]

2230   double randA; // uniform random number on [0, Area under the envelop function]

2231   double fy;    // total area under the envelop function

2232   double fyAct; // actual rate

2233
2234   RateElastic Pos, Neg;
2235   Pos = getRateElasPos(u, T);
2236   Neg = getRateElasNeg(u, T);
2237
2238   // this switch will hold the decision whether k is positive or negative:

2239   // 0 : negative, 1 : positive

2240   int posNegSwitch = 1;
2241
2242   /* process == 5 : qq -> qq (first q is hard, second q is soft)

2243      process == 6 : qg -> qg

2244      process == 7 : gq -> gq

2245      process == 8 : gg -> gg  */
2246
2247   if (process == 5 || process == 7) // for qq or gq

2248   {
2249     // decide whether k shall be positive or negative

2250     // if x (uniform on [0,1]) < area(k<0)/area(all k) then k < 0

2251     if (process == 5) {
2252       if (ZeroOneDistribution(*GetMt19937Generator()) <
2253           Neg.qq / (Neg.qq + Pos.qq))
2254         posNegSwitch = 0;
2255     } else {
2256       if (ZeroOneDistribution(*GetMt19937Generator()) <
2257           Neg.gq / (Neg.gq + Pos.gq))
2258         posNegSwitch = 0;
2259     }
2260
2261     if (posNegSwitch == 1) // if omega > 0

2262     {
2263       do {
2264         randA = ZeroOneDistribution(*GetMt19937Generator()) *
2265                 areaOmega(u, posNegSwitch, process);
2266         y = exp((-1.41428 * pow(10., 9.) * alpha_s -
2267                  8.08158 * pow(10., 8.) * alpha_s * alpha_s +
2268                  2.02327 * pow(10., 9.) * randA) /
2269                 (alpha_s *
2270                  (4.72097 * pow(10., 8.) + 2.6977 * pow(10., 8.) * alpha_s)));
2271
2272         fy = alpha_s / 0.15 * (0.035 + alpha_s * 0.02) / sqrt(y * y);
2273         fyAct = functionOmega(u, y, process);
2274
2275         x = ZeroOneDistribution(*GetMt19937Generator());
2276
2277       } while (x > fyAct / fy);
2278       // reject if x is larger than the ratio fyAct/fy

2279       omega = y;
2280     } else // if omega < 0

2281     {
2282       do {
2283         randA = ZeroOneDistribution(*GetMt19937Generator()) *
2284                 areaOmega(-0.05, posNegSwitch, process);
2285         y = -12. * exp(-30. * randA / (alpha_s * (7. + 4. * alpha_s)));
2286
2287         fy = alpha_s / 0.15 * (0.035 + alpha_s * 0.02) / sqrt(y * y);
2288         fyAct = functionOmega(u, y, process);
2289
2290         x = ZeroOneDistribution(*GetMt19937Generator());
2291
2292       } while (x > fyAct / fy);
2293       // reject if x is larger than the ratio fyAct/fy

2294       omega = y;
2295     }
2296   } else if (process == 6 || process == 8) // for qg or gg

2297   {
2298     // decide whether k shall be positive or negative

2299     // if x (uniform on [0,1]) < area(k<0)/area(all k) then k < 0

2300     if (process == 6) {
2301       if (ZeroOneDistribution(*GetMt19937Generator()) <
2302           Neg.qg / (Neg.qg + Pos.qg))
2303         posNegSwitch = 0;
2304     } else {
2305       if (ZeroOneDistribution(*GetMt19937Generator()) <
2306           Neg.gg / (Neg.gg + Pos.gg))
2307         posNegSwitch = 0;
2308     }
2309
2310     if (posNegSwitch == 1) // if omega > 0

2311     {
2312       do {
2313         randA = ZeroOneDistribution(*GetMt19937Generator()) *
2314                 areaOmega(u, posNegSwitch, process);
2315         y = exp((-2.32591 * pow(10., 17.) * alpha_s -
2316                  1.32909 * pow(10., 17.) * alpha_s * alpha_s +
2317                  2.2183 * pow(10., 17.) * randA) /
2318                 (alpha_s * (7.76406 * pow(10., 16.) +
2319                             4.43661 * pow(10., 16.) * alpha_s)));
2320
2321         fy = 1.5 * alpha_s / 0.15 * (0.035 + alpha_s * 0.02) / sqrt(y * y);
2322         fyAct = functionOmega(u, y, process);
2323
2324         x = ZeroOneDistribution(*GetMt19937Generator());
2325
2326       } while (x > fyAct / fy);
2327       // reject if x is larger than the ratio fyAct/fy

2328       omega = y;
2329     } else // if omega < 0

2330     {
2331       do {
2332         randA = ZeroOneDistribution(*GetMt19937Generator()) *
2333                 areaOmega(-0.05, posNegSwitch, process);
2334         y = -12. * exp(-2.81475 * pow(10., 15.) * randA /
2335                        (alpha_s * (9.85162 * pow(10., 14.) +
2336                                    5.6295 * pow(10., 14.) * alpha_s)));
2337
2338         fy = 1.5 * alpha_s / 0.15 * (0.035 + alpha_s * 0.02) / sqrt(y * y);
2339         fyAct = functionOmega(u, y, process);
2340
2341         x = ZeroOneDistribution(*GetMt19937Generator());
2342
2343       } while (x > fyAct / fy);
2344       // reject if x is larger than the ratio fyAct/fy

2345       omega = y;
2346     }
2347   } else {
2348     JSWARN << "Invalid process number (" << process << ")";
2349
2350     omega = 0.;
2351   }
2352
2353   return omega * T; // omega*T is in [GeV]

2354 }
2355
2356 double Martini::getMomentumTransfer(double pRest, double omega, double T,
2357                                     int process) {
2358   double u = pRest / T; // making arguments in log to be dimensionless

2359   omega /= T;
2360
2361   double q = 0.; // momentum of radiated gluon (dimentionless)

2362   double y;      // q candidate

2363   double x;      // random number, uniform on [0,1]

2364   double randA; // uniform random number on [0, Area under the envelop function]

2365   double fy;    // total area under the envelop function

2366   double fyAct; // actual rate

2367
2368   double A, B;
2369
2370   /* process == 5 : qq -> qq (first q is hard, second q is soft)

2371      process == 6 : qg -> qg

2372      process == 7 : gq -> gq

2373      process == 8 : gg -> gg  */
2374
2375   // small omega using the rejection method

2376   if (omega < 10. && omega > -3.) {
2377     if (process == 5 || process == 7) // for qq or gq

2378     {
2379       A = (0.7 + alpha_s) * 0.0014 *
2380           (1000. + 40. / sqrt(omega * omega) + 10.5 * pow(omega, 4.)) * alpha_s;
2381       B = 2. * sqrt(omega * omega) + 0.01;
2382     } else if (process == 6 || process == 8) // for qg or gg

2383     {
2384       A = (0.7 + alpha_s) * 0.0022 *
2385           (1000. + 40. / sqrt(omega * omega) + 10.5 * pow(omega, 4.)) * alpha_s;
2386       B = 2. * sqrt(omega * omega) + 0.002;
2387     } else {
2388       JSWARN << "Invalid process number (" << process << ")";
2389
2390       A = 0.;
2391       B = 0.;
2392     }
2393
2394     do {
2395       randA = ZeroOneDistribution(*GetMt19937Generator()) *
2396               areaQ(u, omega, process);
2397       y = pow(B, 0.25) *
2398           sqrt(tan((2. * sqrt(B) * randA + A * atan(omega * omega / sqrt(B))) /
2399                    A));
2400
2401       fy = A * y / (pow(y, 4.) + B);
2402       fyAct = functionQ(u, omega, y, process);
2403
2404       x = ZeroOneDistribution(*GetMt19937Generator());
2405
2406     } while (x > fyAct / fy);
2407     // reject if x is larger than the ratio fyAct/fy

2408     q = y;
2409   }
2410   // large omega using the Metropolis method

2411   else if (omega < 40.) {
2412     double g = 0, g_new = 0;
2413     double ratio;
2414     double y_new;
2415
2416     // the ranges in which the variables u and phi need to be sampled

2417     const double y_min = sqrt(omega * omega);
2418     const double y_max = u;
2419
2420     int count = 0;
2421     // randomly select initial values of q=y, such that

2422     do {
2423       y = y_min + ZeroOneDistribution(*GetMt19937Generator()) * (y_max - y_min);
2424       g = functionQ(u, omega, y, process);
2425
2426       // if no y having non-zero g is found, give up here.

2427       if (count > 100)
2428         return 0.;
2429       count++;
2430     } while (g == 0.);
2431
2432     // number of steps in the Markov chain

2433     const int n_steps = 500;
2434
2435     for (int i = 0; i < n_steps; i++) {
2436       do {
2437         y_new = y_min +
2438                 ZeroOneDistribution(*GetMt19937Generator()) * (y_max - y_min);
2439       } while (y_new < y_min || y_new > y_max);
2440       // check that the new value is in range

2441
2442       // calculate the function at the proposed point

2443       g_new = functionQ(u, omega, y_new, process);
2444       // ratio of g(y_new)/g(y)

2445       ratio = g_new / g;
2446
2447       // accept if probability g(y_new)/g(y) is larger than randon number

2448       if (ZeroOneDistribution(*GetMt19937Generator()) < ratio) {
2449         y = y_new;
2450         g = g_new;
2451       }
2452     }
2453     q = y;
2454   }
2455   // set collinear for too large omega

2456   else {
2457     q = omega;
2458   }
2459
2460   return q * T; // q*T is in [GeV]

2461 }
2462
2463 // calculates the area under the envelope function when using the rejection method

2464 // integrals had been solved analytically before

2465 double Martini::areaOmega(double u, int posNegSwitch, int process) {
2466   if (process == 5 || process == 7) {
2467     if (posNegSwitch == 1)
2468       return (0.0333333 * alpha_s * (7. + 4. * alpha_s) * (2.99573 + log(u)));
2469     else
2470       return (-0.133333 * alpha_s * (1.75 + alpha_s) * log(-0.0833333 * u));
2471   } else if (process == 6 || process == 8) {
2472     if (posNegSwitch == 1)
2473       return (0.05 * alpha_s * (7. + 4. * alpha_s) * (2.99573 + log(u)));
2474     else
2475       return (-0.2 * alpha_s * (1.75 + alpha_s) * log(-0.0833333 * u));
2476   } else {
2477     JSWARN << "Invalid process number (" << process << ")";
2478   }
2479
2480   return 0.;
2481 }
2482
2483 double Martini::areaQ(double u, double omega, int process) {
2484   double A, B;
2485   double areaQ;
2486
2487   if (process == 5 || process == 7) {
2488     A = (0.7 + alpha_s - 0.3) * 0.0014 * alpha_s *
2489         (1000. + 40. / sqrt(omega * omega) + 10.5 * pow(omega, 4.));
2490     B = 2. * sqrt(omega * omega) + 0.01;
2491   } else if (process == 6 || process == 8) {
2492     A = (0.7 + alpha_s - 0.3) * 0.0022 * alpha_s *
2493         (1000. + 40. / sqrt(omega * omega) + 10.5 * pow(omega, 4.));
2494     B = 2. * sqrt(omega * omega) + 0.002;
2495   } else {
2496     JSWARN << "Invalid process number (" << process << ")";
2497
2498     A = 0.;
2499     B = 0.;
2500   }
2501   areaQ = (0.5 * A * (atan(u * u / sqrt(B)) - atan(omega * omega / sqrt(B)))) /
2502           sqrt(B);
2503
2504   return areaQ;
2505 }
2506
2507 FourVector Martini::getNewMomentumElas(FourVector pVec, double omega,
2508                                        double q) {
2509   FourVector pVecNew, pVecNewTemp;
2510   FourVector etVec, qtVec, qlVec;
2511   FourVector r;
2512   double qt, ql;
2513   double cosTheta_pq;
2514   double pAbs = pVec.t();
2515   double phi;
2516   double M[3][3]; //rotation matrix

2517   double u;
2518   double xx, yy, zz, tt;
2519
2520   if (omega == q) {
2521     xx = pVec.x() * (pAbs - omega) / pAbs;
2522     yy = pVec.y() * (pAbs - omega) / pAbs;
2523     zz = pVec.z() * (pAbs - omega) / pAbs;
2524     tt = pVec.t() * (pAbs - omega) / pAbs;
2525
2526     pVecNew.Set(xx, yy, zz, tt);
2527     return pVecNew;
2528   }
2529
2530   cosTheta_pq = (-omega * omega + 2. * pAbs * omega + q * q) / (2. * pAbs * q);
2531   qt = q * sqrt(1. - cosTheta_pq * cosTheta_pq); // transverse momentum transfer

2532   ql = q * cosTheta_pq;
2533
2534   if (pVec.y() * pVec.y() > pVec.x() * pVec.x()) {
2535     xx = 0.;
2536     yy = -pVec.z();
2537     zz = pVec.y();
2538   } else {
2539     xx = pVec.z();
2540     yy = 0.;
2541     zz = -pVec.x();
2542   }
2543
2544   tt = sqrt(xx * xx + yy * yy + zz * zz);
2545   etVec.Set(xx / tt, yy / tt, zz / tt, 1.); // normalized to 1

2546
2547   // the transverse transferred momentum vector

2548   qtVec.Set(etVec.x() * qt, etVec.y() * qt, etVec.z() * qt, etVec.t() * qt);
2549   // the longuitudinal transferred momentum vector

2550   qlVec.Set(pVec.x() / pAbs * ql, pVec.y() / pAbs * ql, pVec.z() / pAbs * ql,
2551             pVec.t() / pAbs * ql);
2552
2553   pVecNewTemp = pVec;
2554   pVecNewTemp -= qtVec; // change transverse momentum

2555   pVecNewTemp -= qlVec; // change longitudinal momentum

2556
2557   phi = 2. * M_PI * ZeroOneDistribution(*GetMt19937Generator());
2558   r.Set(pVec.x() / pVec.t(), pVec.y() / pVec.t(), pVec.z() / pVec.t(), 1.);
2559   u = 1. - cos(phi);
2560
2561   // define the rotation matrix for rotations around pvecRest

2562   M[0][0] = r.x() * r.x() * u + cos(phi);
2563   M[1][0] = r.x() * r.y() * u - r.z() * sin(phi);
2564   M[2][0] = r.x() * r.z() * u + r.y() * sin(phi);
2565
2566   M[0][1] = r.y() * r.x() * u + r.z() * sin(phi);
2567   M[1][1] = r.y() * r.y() * u + cos(phi);
2568   M[2][1] = r.y() * r.z() * u - r.x() * sin(phi);
2569
2570   M[0][2] = r.z() * r.x() * u - r.y() * sin(phi);
2571   M[1][2] = r.z() * r.y() * u + r.x() * sin(phi);
2572   M[2][2] = r.z() * r.z() * u + cos(phi);
2573
2574   xx = M[0][0] * pVecNewTemp.x() + M[0][1] * pVecNewTemp.y() +
2575        M[0][2] * pVecNewTemp.z();
2576   yy = M[1][0] * pVecNewTemp.x() + M[1][1] * pVecNewTemp.y() +
2577        M[1][2] * pVecNewTemp.z();
2578   zz = M[2][0] * pVecNewTemp.x() + M[2][1] * pVecNewTemp.y() +
2579        M[2][2] * pVecNewTemp.z();
2580   tt = sqrt(xx * xx + yy * yy + zz * zz);
2581
2582   pVecNew.Set(xx, yy, zz, tt);
2583   return pVecNew;
2584 }
2585
2586 FourVector Martini::getThermalVec(FourVector qVec, double T, int id) {
2587   int kind = 1; // 1: fermion, 2: boson

2588   if (fabs(id) < 4)
2589     kind = 1;
2590   else if (id == 21)
2591     kind = -1;
2592
2593   FourVector kVec;           // new thermal parton's momentum

2594   double k, k_min;           // (minimum) mangitude of thermal momentum

2595   double cosTh, sinTh;       // angle between vector q and k

2596   double u1[3], u2[3];       // unit vector in rotation axis

2597   double M1[3][3], M2[3][3]; // rotation matrix

2598   double xx, yy, zz, tt;     // components of 4-vector

2599
2600   double q =
2601       sqrt(qVec.x() * qVec.x() + qVec.y() * qVec.y() + qVec.z() * qVec.z());
2602   double omega = qVec.t();
2603
2604   // if momentum transfer is on-shell or time-like

2605   // return null FourVector

2606   if (q - fabs(omega) < 1e-5) {
2607     return FourVector(0., 0., 0., 0.);
2608   }
2609
2610   // minimum momenum of thermal parton k that makes recoil parton on-shell

2611   k_min = (q - omega) / 2.;
2612
2613   // sampled magnitude of thermal momentum

2614   k = getThermal(k_min, T, kind);
2615
2616   cosTh = (2. * k * omega - q * q + omega * omega) / (2. * k * q);
2617   sinTh = sqrt(1. - cosTh * cosTh);
2618
2619   // choose a unit vector perpendicular to qVec with which qVec is rotated by theta

2620   // to get the direction of thermal parton kVec

2621   double norm = sqrt(qVec.x() * qVec.x() + qVec.y() * qVec.y());
2622   if (norm > 1e-10) {
2623     u1[0] = qVec.y() / norm;
2624     u1[1] = -qVec.x() / norm;
2625     u1[2] = 0.;
2626   }
2627   // if momentum transfer is z-direction u1 is set to x-direction

2628   else {
2629     u1[0] = 1.;
2630     u1[1] = 0.;
2631     u1[2] = 0.;
2632   }
2633
2634   // define the rotation matrix for theta rotations

2635   M1[0][0] = u1[0] * u1[0] * (1. - cosTh) + cosTh;
2636   M1[1][0] = u1[0] * u1[1] * (1. - cosTh) - u1[2] * sinTh;
2637   M1[2][0] = u1[0] * u1[2] * (1. - cosTh) + u1[1] * sinTh;
2638
2639   M1[0][1] = u1[1] * u1[0] * (1. - cosTh) + u1[2] * sinTh;
2640   M1[1][1] = u1[1] * u1[1] * (1. - cosTh) + cosTh;
2641   M1[2][1] = u1[1] * u1[2] * (1. - cosTh) - u1[0] * sinTh;
2642
2643   M1[0][2] = u1[2] * u1[0] * (1. - cosTh) - u1[1] * sinTh;
2644   M1[1][2] = u1[2] * u1[1] * (1. - cosTh) + u1[0] * sinTh;
2645   M1[2][2] = u1[2] * u1[2] * (1. - cosTh) + cosTh;
2646
2647   // get thermal parton's momentum by rotating theta from qVec

2648   xx = M1[0][0] * qVec.x() + M1[0][1] * qVec.y() + M1[0][2] * qVec.z();
2649   yy = M1[1][0] * qVec.x() + M1[1][1] * qVec.y() + M1[1][2] * qVec.z();
2650   zz = M1[2][0] * qVec.x() + M1[2][1] * qVec.y() + M1[2][2] * qVec.z();
2651
2652   // momentum is normalized to k

2653   tt = sqrt(xx * xx + yy * yy + zz * zz);
2654   xx /= tt / k;
2655   yy /= tt / k;
2656   zz /= tt / k;
2657
2658   kVec.Set(xx, yy, zz, k);
2659
2660   // rotate kVec in random azimuthal angle with respect to the momentum transfer qVec

2661   double phi = 2. * M_PI * ZeroOneDistribution(*GetMt19937Generator());
2662
2663   xx = qVec.x();
2664   yy = qVec.y();
2665   zz = qVec.z();
2666   tt = sqrt(xx * xx + yy * yy + zz * zz);
2667
2668   u2[0] = xx / tt;
2669   u2[1] = yy / tt;
2670   u2[2] = zz / tt;
2671
2672   M2[0][0] = u2[0] * u2[0] * (1. - cos(phi)) + cos(phi);
2673   M2[1][0] = u2[0] * u2[1] * (1. - cos(phi)) - u2[2] * sin(phi);
2674   M2[2][0] = u2[0] * u2[2] * (1. - cos(phi)) + u2[1] * sin(phi);
2675
2676   M2[0][1] = u2[1] * u2[0] * (1. - cos(phi)) + u2[2] * sin(phi);
2677   M2[1][1] = u2[1] * u2[1] * (1. - cos(phi)) + cos(phi);
2678   M2[2][1] = u2[1] * u2[2] * (1. - cos(phi)) - u2[0] * sin(phi);
2679
2680   M2[0][2] = u2[2] * u2[0] * (1. - cos(phi)) - u2[1] * sin(phi);
2681   M2[1][2] = u2[2] * u2[1] * (1. - cos(phi)) + u2[0] * sin(phi);
2682   M2[2][2] = u2[2] * u2[2] * (1. - cos(phi)) + cos(phi);
2683
2684   xx = M2[0][0] * kVec.x() + M2[0][1] * kVec.y() + M2[0][2] * kVec.z();
2685   yy = M2[1][0] * kVec.x() + M2[1][1] * kVec.y() + M2[1][2] * kVec.z();
2686   zz = M2[2][0] * kVec.x() + M2[2][1] * kVec.y() + M2[2][2] * kVec.z();
2687   tt = sqrt(xx * xx + yy * yy + zz * zz);
2688
2689   kVec.Set(xx, yy, zz, tt);
2690
2691   return kVec;
2692 }
2693
2694 double Martini::getThermal(double k_min, double T, int kind) {
2695   // this algorithm uses 5.5/(1+px2) as envelope function and then uses the rejection method

2696   // tan (Pi/2*x) is the inverse of the integral of the envelope function

2697   // this can be improved

2698   // kind: -1 = Boson

2699   //        1 = Fermion

2700   double p;
2701   double gm = 5.5;
2702   int repeat = 1;
2703   double gx = 0.;
2704   double px;
2705   double px2;
2706   double ex;
2707
2708   do {
2709     px = k_min / T +
2710          tan(M_PI * ZeroOneDistribution(*GetMt19937Generator()) / 2.0);
2711     px2 = px * px;
2712     ex = sqrt(px2);
2713     if (kind == 0)
2714       gx = px2 * (1.0 + px2) * exp(-ex);
2715     else if (kind == -1)
2716       gx = px2 * (1.0 + px2) * 1 / (exp(ex) - 1);
2717     else if (kind == +1)
2718       gx = px2 * (1.0 + px2) * 1 / (exp(ex) + 1);
2719     if (ZeroOneDistribution(*GetMt19937Generator()) < gx / gm) {
2720       repeat = 0;
2721       p = ex;
2722     }
2723   } while (repeat);
2724
2725   return p * T; // p*T is in [GeV]

2726 }
2727
2728 // Reads in the binary stored file of dGamma values

2729 void Martini::readRadiativeRate(Gamma_info *dat, dGammas *Gam) {
2730   FILE *rfile;
2731   string filename;
2732   filename = PathToTables + "/radgamma";
2733
2734   JSINFO << "Reading rates of inelastic collisions from file ";
2735   JSINFO << filename.c_str() << " ... ";
2736   size_t bytes_read;
2737
2738   rfile = fopen(filename.c_str(), "rb");
2739   if (rfile == NULL) {
2740     JSWARN << "[readRadiativeRate]: ERROR: Unable to open file " << filename;
2741     throw std::runtime_error(
2742         "[readRadiativeRate]: ERROR: Unable to open Radiative rate file");
2743   }
2744
2745   bytes_read = fread((char *)(&dat->ddf), sizeof(double), 1, rfile);
2746   bytes_read = fread((char *)(&dat->dda), sizeof(double), 1, rfile);
2747   bytes_read = fread((char *)(&dat->dcf), sizeof(double), 1, rfile);
2748   bytes_read = fread((char *)(&dat->dca), sizeof(double), 1, rfile);
2749   bytes_read = fread((char *)(&dat->Nc), sizeof(int), 1, rfile);
2750   bytes_read = fread((char *)(&dat->Nf), sizeof(int), 1, rfile);
2751   bytes_read = fread((char *)(&dat->BetheHeitler), sizeof(int), 1, rfile);
2752   bytes_read = fread((char *)(&dat->BDMPS), sizeof(int), 1, rfile);
2753   bytes_read = fread((char *)(&dat->include_gluons), sizeof(int), 1, rfile);
2754   bytes_read = fread((char *)Gam->qqg, sizeof(double), NP * NK, rfile);
2755   bytes_read = fread((char *)Gam->tau_qqg, sizeof(double), NP * NK, rfile);
2756   bytes_read = fread((char *)Gam->gqq, sizeof(double), NP * NK, rfile);
2757   bytes_read = fread((char *)Gam->tau_gqq, sizeof(double), NP * NK, rfile);
2758   bytes_read = fread((char *)Gam->ggg, sizeof(double), NP * NK, rfile);
2759   bytes_read = fread((char *)Gam->tau_ggg, sizeof(double), NP * NK, rfile);
2760   bytes_read = fread((char *)Gam->qqgamma, sizeof(double), NP * NK, rfile);
2761   bytes_read = fread((char *)Gam->tau_qqgamma, sizeof(double), NP * NK, rfile);
2762   fclose(rfile);
2763
2764   dat->Nf = nf;
2765   dat->dp = 0.05;
2766   dat->p_max = 20;
2767   dat->p_min =
2768       0; // exp(LogEmin) set to zero because my array starts at zero!...

2769   dat->n_p = static_cast<int>(
2770       1.001 + dat->p_max / dat->dp); // np = int(0.4 + 121 / 0.5) = 241

2771   dat->p_max = dat->dp * dat->n_p;   // p_max = 0.5*241 = 120.5

2772   dat->n_pmin = static_cast<int>(
2773       1.001 + dat->p_min / dat->dp);  // np_min = int(0.4 + 3.3 / 0.5) = 7

2774   dat->n_p -= dat->n_pmin - 1;        // n_p = 241 - (7 - 1) = 235

2775   dat->p_min = dat->dp * dat->n_pmin; // p_min = 0.5 * 7 = 3.5

2776   dat->n_kmin = 1 + 2 * (static_cast<int>(2. / dat->dp));
2777   dat->k_min = -dat->dp * dat->n_kmin;
2778   dat->n_k = static_cast<int>((8 + dat->p_max) / (2 * dat->dp));
2779   dat->k_max = 2 * dat->dp * (dat->n_k - 1) + dat->k_min;
2780 }
2781
2782 void Martini::readElasticRateOmega() {
2783   ifstream fin;
2784   string filename[2];
2785
2786   double as, omega;
2787   double dGamma;
2788
2789   // open files with data to read in:

2790   filename[0] = PathToTables + "/logEnDtrqq";
2791   filename[1] = PathToTables + "/logEnDtrqg";
2792
2793   JSINFO << "Reading rates of elastic collisions from files";
2794   JSINFO << filename[0];
2795   JSINFO << filename[1] << " ...";
2796
2797   fin.open(filename[0].c_str(), ios::in);
2798   if (!fin) {
2799     JSWARN << "[readElasticRateOmega]: ERROR: Unable to open file "
2800            << filename[0];
2801     throw std::runtime_error(
2802         "[readElasticRateQ]: ERROR: Unable to open ElasticRateOmega file");
2803   }
2804
2805   int ik = 0;
2806   while (!fin.eof()) {
2807     fin >> as;
2808     fin >> omega;
2809     fin >> dGamma;
2810     dGamma_qq->push_back(dGamma);
2811
2812     ik++;
2813   }
2814   fin.close();
2815
2816   fin.open(filename[1].c_str(), ios::in);
2817   if (!fin) {
2818     JSWARN << "[readElasticRateOmega]: ERROR: Unable to open file "
2819            << filename[1];
2820     throw std::runtime_error(
2821         "[readElasticRateQ]: ERROR: Unable to open ElasticRateOmega file");
2822   }
2823
2824   ik = 0;
2825   while (!fin.eof()) {
2826     fin >> as;
2827     fin >> omega;
2828     fin >> dGamma;
2829     dGamma_qg->push_back(dGamma);
2830
2831     ik++;
2832   }
2833   fin.close();
2834 }
2835
2836 void Martini::readElasticRateQ() {
2837   ifstream fin;
2838   string filename[2];
2839
2840   double as, omega, q;
2841   double dGamma;
2842
2843   // open files with data to read in:

2844   filename[0] = PathToTables + "/logEnDqtrqq";
2845   filename[1] = PathToTables + "/logEnDqtrqg";
2846
2847   JSINFO << "Reading rates of elastic collisions from files";
2848   JSINFO << filename[0];
2849   JSINFO << filename[1] << " ...";
2850
2851   fin.open(filename[0].c_str(), ios::in);
2852   if (!fin) {
2853     JSWARN << "[readElasticRateQ]: ERROR: Unable to open file " << filename[0];
2854     throw std::runtime_error(
2855         "[readElasticRateQ]: ERROR: Unable to open ElasticRateQ file");
2856   }
2857
2858   int ik = 0;
2859   while (!fin.eof()) {
2860     fin >> as;
2861     fin >> omega;
2862     fin >> q;
2863     fin >> dGamma;
2864     dGamma_qq_q->push_back(dGamma);
2865
2866     ik++;
2867   }
2868   fin.close();
2869
2870   fin.open(filename[1].c_str(), ios::in);
2871   if (!fin) {
2872     JSWARN << "[readElasticRateQ]: ERROR: Unable to open file " << filename[1];
2873     throw std::runtime_error(
2874         "[readElasticRateQ]: ERROR: Unable to open ElasticRateQ file");
2875   }
2876
2877   ik = 0;
2878   while (!fin.eof()) {
2879     fin >> as;
2880     fin >> omega;
2881     fin >> q;
2882     fin >> dGamma;
2883     dGamma_qg_q->push_back(dGamma);
2884
2885     ik++;
2886   }
2887   fin.close();
2888 }
2889
2890 double Martini::getRate_qqg(double p, double k) {
2891   return use_table(p, k, Gam.qqg, 0);
2892 }
2893
2894 double Martini::getRate_gqq(double p, double k) {
2895   if (k < p / 2.)
2896     return use_table(p, k, Gam.gqq, 1);
2897   else
2898     return 0.;
2899 }
2900
2901 double Martini::getRate_ggg(double p, double k) {
2902   if (k < p / 2.)
2903     return use_table(p, k, Gam.ggg, 2);
2904   else
2905     return 0.;
2906 }
2907
2908 double Martini::getRate_qqgamma(double p, double k) {
2909   return use_table(p, k, Gam.qqgamma, 3);
2910 }
2911
2912 double Martini::use_table(double p, double k, double dGamma[NP][NK],
2913                           int which_kind)
2914 /* Uses the lookup table and simple interpolation to get the value

2915    of dGamma/dkdx at some value of p,k.

2916    This works by inverting the relations between (p,k) and (n_p,n_k)

2917    used in building the table, to find out what continuous values

2918    of n_p, n_k should be considered; then linearly interpolates.     */
2919 {
2920   double a, b, result; // fraction of way from corner of box

2921   int n_p, n_k;        // location of corner of box

2922
2923   // out of range

2924   if ((p < 4.01) || (p > 46000.) || (k < -12.) || (k > p + 12.))
2925     return 0.;
2926
2927   if ((which_kind % 3) && (k > p / 2))
2928     k = p - k; // Take advantage of symmetry in these cases

2929
2930   a = 24.7743737154026 * log(p * 0.2493765586034912718l);
2931   n_p = (int)a;
2932   a -= n_p;
2933   if (k < 2.) {
2934     if (k < -1) {
2935       if (k < -2)
2936         b = 60. + 5. * k;
2937       else
2938         b = 70. + 10. * k;
2939     } else {
2940       if (k < 1.)
2941         b = 80. + 20. * k;
2942       else
2943         b = 90. + 10. * k;
2944     }
2945   } else if (k < p - 2.) { /* This is that tricky middle ground. */
2946     b = 190. - 10. * log(1.000670700260932956l /
2947                              (0.0003353501304664781l + (k - 2.) / (p - 4.)) -
2948                          1.);
2949   } else {
2950     if (k < p + 1.) {
2951       if (k < p - 1.)
2952         b = 290. + 10. * (k - p);
2953       else
2954         b = 300. + 20. * (k - p);
2955     } else {
2956       if (k < p + 2.)
2957         b = 310. + 10. * (k - p);
2958       else
2959         b = 320. + 5. * (k - p);
2960     }
2961   }
2962
2963   n_k = (int)b;
2964   b -= n_k;
2965   result = (1. - a) * ((1. - b) * dGamma[n_p][n_k] + b * dGamma[n_p][n_k + 1]) +
2966            a * ((1. - b) * dGamma[n_p + 1][n_k] + b * dGamma[n_p + 1][n_k + 1]);
2967
2968   if (std::abs(k) > 0.001) // Avoid division by 0, should never get asked for

2969   {
2970     switch (which_kind) {
2971     case 0:
2972       result /= k;
2973       if (k < 20.)
2974         result /= 1. - exp(-k);
2975       if (k > p - 20.)
2976         result /= 1. + exp(k - p);
2977       break;
2978     case 1:
2979       result /= p;
2980       if (k < 20.)
2981         result /= 1 + exp(-k);
2982       if (k > p - 20.)
2983         result /= 1. + exp(k - p);
2984       break;
2985     case 2:
2986       result /= k * (p - k) / p;
2987       if (k < 20.)
2988         result /= 1. - exp(-k);
2989       if (k > p - 20.)
2990         result /= 1. - exp(k - p);
2991       break;
2992     case 3:
2993       result /= k;
2994       if (k < 0)
2995         result = 0.;
2996       if (k > p - 20.)
2997         result /= 1. + exp(k - p);
2998       break;
2999     }
3000   }
3001
3002   return result;
3003 }
3004
3005 double Martini::getElasticRateOmega(double u, double omega, int process) {
3006   if ((omega > 0. && omega < u) || (omega < 0. && omega > -u))
3007     return use_elastic_table_omega(omega, process);
3008
3009   return 0.;
3010 }
3011
3012 double Martini::getElasticRateQ(double u, double omega, double q, int process) {
3013   if (q > std::abs(omega) &&
3014       ((omega > 0 && omega < u) || (omega < 0 && omega > -u)))
3015     return use_elastic_table_q(omega, q, process);
3016
3017   return 0.;
3018 }
3019
3020 double Martini::use_elastic_table_omega(double omega, int which_kind)
3021 /* Uses the lookup table and simple interpolation to get the value

3022    of dGamma/domega at some value of omega and alpha_s. */
3023 {
3024   double result;
3025   double alphaFrac, omegaFrac;
3026   int iOmega, iAlphas;
3027   int position, positionAlphaUp, positionOmegaUp, positionAlphaUpOmegaUp;
3028   double rate, rateAlphaUp, rateOmegaUp, rateAlphaUpOmegaUp;
3029   double rateOmegaAv, rateAlphaUpOmegaAv;
3030
3031   if (omega > 0.)
3032     iOmega = Nomega / 2 + floor((log(omega) + 5) / omegaStep);
3033   else
3034     iOmega = Nomega / 2 - ceil((log(-omega) + 5) / omegaStep) - 1;
3035   iAlphas = floor((alpha_s - 0.15) / alphaStep);
3036
3037   position = iOmega + Nomega * (iAlphas);
3038   positionAlphaUp = iOmega + Nomega * (iAlphas + 1);
3039   positionOmegaUp = iOmega + 1 + Nomega * (iAlphas);
3040   positionAlphaUpOmegaUp = iOmega + 1 + Nomega * (iAlphas + 1);
3041
3042   alphaFrac = (alpha_s - (floor((alpha_s - alphaMin) / alphaStep) * alphaStep +
3043                           alphaMin)) /
3044               alphaStep;
3045   if (omega > 0.) {
3046     if (exp(ceil((log(omega) + 5) / omegaStep) * omegaStep - 5) !=
3047         exp(floor((log(omega) + 5) / omegaStep) * omegaStep - 5))
3048       omegaFrac =
3049           (omega - (exp(floor((log(omega) + 5) / omegaStep) * omegaStep - 5))) /
3050           ((exp(ceil((log(omega) + 5) / omegaStep) * omegaStep - 5)) -
3051            exp(floor((log(omega) + 5) / omegaStep) * omegaStep - 5));
3052     else
3053       omegaFrac = 0.;
3054   } else {
3055     if (exp(ceil((log(-omega) + 5) / omegaStep) * omegaStep - 5) !=
3056         exp(floor((log(-omega) + 5) / omegaStep) * omegaStep - 5))
3057       omegaFrac =
3058           (-omega -
3059            (exp(floor((log(-omega) + 5) / omegaStep) * omegaStep - 5))) /
3060           ((exp(ceil((log(-omega) + 5) / omegaStep) * omegaStep - 5)) -
3061            exp(floor((log(-omega) + 5) / omegaStep) * omegaStep - 5));
3062     else
3063       omegaFrac = 0.;
3064   }
3065
3066   if (which_kind == 5 || which_kind == 7) {
3067     if (position > 0 && iAlphas < Nalphas && iOmega < Nomega)
3068       rate = dGamma_qq->at(position);
3069     else
3070       rate = 0.;
3071
3072     if (iAlphas + 1 < Nalphas)
3073       rateAlphaUp = dGamma_qq->at(positionAlphaUp);
3074     else
3075       rateAlphaUp = rate;
3076
3077     if (iOmega + 1 < Nomega)
3078       rateOmegaUp = dGamma_qq->at(positionOmegaUp);
3079     else
3080       rateOmegaUp = rate;
3081
3082     if (iAlphas < Nalphas && iOmega < Nomega)
3083       rateAlphaUpOmegaUp = dGamma_qq->at(positionAlphaUpOmegaUp);
3084     else
3085       rateAlphaUpOmegaUp = rate;
3086   } else {
3087     if (position > 0 && iAlphas < Nalphas && iOmega < Nomega)
3088       rate = dGamma_qg->at(position);
3089     else
3090       rate = 0.;
3091
3092     if (iAlphas + 1 < Nalphas)
3093       rateAlphaUp = dGamma_qg->at(positionAlphaUp);
3094     else
3095       rateAlphaUp = rate;
3096
3097     if (iOmega + 1 < Nomega)
3098       rateOmegaUp = dGamma_qg->at(positionOmegaUp);
3099     else
3100       rateOmegaUp = rate;
3101
3102     if (iAlphas < Nalphas && iOmega < Nomega)
3103       rateAlphaUpOmegaUp = dGamma_qg->at(positionAlphaUpOmegaUp);
3104     else
3105       rateAlphaUpOmegaUp = rate;
3106   }
3107
3108   if (omega > 0.) {
3109     rateOmegaAv = (1. - omegaFrac) * rate + omegaFrac * rateOmegaUp;
3110     rateAlphaUpOmegaAv =
3111         (1. - omegaFrac) * rateAlphaUp + omegaFrac * rateAlphaUpOmegaUp;
3112   } else {
3113     rateOmegaAv = (omegaFrac)*rate + (1. - omegaFrac) * rateOmegaUp;
3114     rateAlphaUpOmegaAv =
3115         (omegaFrac)*rateAlphaUp + (1. - omegaFrac) * rateAlphaUpOmegaUp;
3116   }
3117   result = (1. - alphaFrac) * rateOmegaAv + alphaFrac * rateAlphaUpOmegaAv;
3118
3119   // leave out the *9./4. for processes 6 and 8 to use the same envelope later

3120   return result;
3121 }
3122
3123 double Martini::use_elastic_table_q(double omega, double q, int which_kind)
3124 /* Uses the lookup table and simple interpolation to get the value

3125    of dGamma/domegadq at some value of omega, q, and alpha_s. */
3126 {
3127   double result;
3128   double alphaFrac, omegaFrac, qFrac;
3129   int iOmega, iAlphas, iQ;
3130   int position, positionAlphaUp, positionOmegaUp, positionAlphaUpOmegaUp;
3131   int positionQUp, positionAlphaUpQUp, positionOmegaUpQUp,
3132       positionAlphaUpOmegaUpQUp;
3133   int position2QUp;
3134   double rate, rateAlphaUp, rateOmegaUp, rateAlphaUpOmegaUp;
3135   double rateQUp, rateAlphaUpQUp, rateOmegaUpQUp, rateAlphaUpOmegaUpQUp;
3136   double rate2QUp, rateAlphaUp2QUp, rateOmegaUp2QUp, rateAlphaUpOmegaUp2QUp;
3137   double rateOmegaAv, rateAlphaUpOmegaAv, rateQUpOmegaAv, rateAlphaUpQUpOmegaAv;
3138   double rate2QUpOmegaAv, rateAlphaUp2QUpOmegaAv;
3139   double rateQAv, rateAlphaUpQAv;
3140   double slope, slopeAlphaUp;
3141
3142   rate2QUp = 0.;
3143   rateAlphaUp2QUp = 0.;
3144   rateOmegaUp2QUp = 0.;
3145   rateAlphaUpOmegaUp2QUp = 0.;
3146   rateOmegaAv = 0.;
3147   rateAlphaUpOmegaAv = 0.;
3148   rateQUpOmegaAv = 0.;
3149   rateAlphaUpQUpOmegaAv = 0.;
3150   rate2QUpOmegaAv = 0.;
3151   rateAlphaUp2QUpOmegaAv = 0.;
3152
3153   if (omega > 0.)
3154     iOmega = Nomega / 2 + floor((log(omega) + 5) / omegaStep);
3155   else
3156     iOmega = Nomega / 2 - ceil((log(-omega) + 5) / omegaStep) - 1;
3157   iQ = floor((log(q) + 5) / qStep + 0.0001);
3158   iAlphas = floor((alpha_s - 0.15) / alphaStep + 0.0001);
3159
3160   position = iQ + Nq * (iOmega + Nomega * (iAlphas));
3161   positionAlphaUp = iQ + Nq * (iOmega + Nomega * (iAlphas + 1));
3162   positionOmegaUp = iQ + Nq * (iOmega + 1 + Nomega * (iAlphas));
3163   positionQUp = iQ + 1 + Nq * (iOmega + Nomega * (iAlphas));
3164   position2QUp = iQ + 2 + Nq * (iOmega + Nomega * (iAlphas));
3165   positionAlphaUpOmegaUp = iQ + Nq * (iOmega + 1 + Nomega * (iAlphas + 1));
3166   positionAlphaUpQUp = iQ + 1 + Nq * (iOmega + Nomega * (iAlphas + 1));
3167   positionOmegaUpQUp = iQ + 1 + Nq * (iOmega + 1 + Nomega * (iAlphas));
3168   positionAlphaUpOmegaUpQUp =
3169       iQ + 1 + Nq * (iOmega + 1 + Nomega * (iAlphas + 1));
3170
3171   alphaFrac = (alpha_s - (floor((alpha_s - alphaMin) / alphaStep) * alphaStep +
3172                           alphaMin)) /
3173               alphaStep;
3174   if (omega > 0.) {
3175     if (exp(ceil((log(omega) + 5) / omegaStep) * omegaStep - 5) !=
3176         exp(floor((log(omega) + 5) / omegaStep) * omegaStep - 5))
3177       omegaFrac =
3178           (omega - (exp(floor((log(omega) + 5) / omegaStep) * omegaStep - 5))) /
3179           ((exp(ceil((log(omega) + 5) / omegaStep) * omegaStep - 5)) -
3180            exp(floor((log(omega) + 5) / omegaStep) * omegaStep - 5));
3181     else
3182       omegaFrac = 0.;
3183   } else {
3184     if (exp(ceil((log(-omega) + 5) / omegaStep) * omegaStep - 5) !=
3185         exp(floor((log(-omega) + 5) / omegaStep) * omegaStep - 5))
3186       omegaFrac =
3187           (-omega -
3188            (exp(floor((log(-omega) + 5) / omegaStep) * omegaStep - 5))) /
3189           ((exp(ceil((log(-omega) + 5) / omegaStep) * omegaStep - 5)) -
3190            exp(floor((log(-omega) + 5) / omegaStep) * omegaStep - 5));
3191     else
3192       omegaFrac = 0.;
3193   }
3194
3195   // interpolate the logs linearly for large omegas

3196   // since there the spectrum is dropping exp in q

3197   if (omega > 20.) {
3198     qFrac = (log(q) - (floor((log(q) + 5.) / qStep) * qStep - 5.)) / qStep;
3199   }
3200   // linear interpolation in q for small omegas

3201   else {
3202     if (exp(ceil((log(q) + 5.) / qStep) * qStep - 5.) !=
3203         exp(floor((log(q) + 5.) / qStep) * qStep - 5.))
3204       qFrac = (q - (exp(floor((log(q) + 5.) / qStep) * qStep - 5.))) /
3205               ((exp(ceil((log(q) + 5.) / qStep) * qStep - 5.)) -
3206                exp(floor((log(q) + 5.) / qStep) * qStep - 5.));
3207     else
3208       qFrac = 0.;
3209   }
3210
3211   if (which_kind == 5 || which_kind == 7) {
3212     if (position >= 0 && iAlphas < Nalphas && iOmega < Nomega && iQ < Nq)
3213       rate = dGamma_qq_q->at(position);
3214     else
3215       rate = 0.;
3216
3217     if (iAlphas + 1 < Nalphas)
3218       rateAlphaUp = dGamma_qq_q->at(positionAlphaUp);
3219     else
3220       rateAlphaUp = rate;
3221
3222     if (iOmega + 1 < Nomega)
3223       rateOmegaUp = dGamma_qq_q->at(positionOmegaUp);
3224     else
3225       rateOmegaUp = rate;
3226
3227     if (iQ + 1 < Nq)
3228       rateQUp = dGamma_qq_q->at(positionQUp);
3229     else
3230       rateQUp = rate;
3231
3232     if (iAlphas < Nalphas && iOmega < Nomega)
3233       rateAlphaUpOmegaUp = dGamma_qq_q->at(positionAlphaUpOmegaUp);
3234     else
3235       rateAlphaUpOmegaUp = rate;
3236
3237     if (iAlphas < Nalphas && iQ < Nq)
3238       rateAlphaUpQUp = dGamma_qq_q->at(positionAlphaUpQUp);
3239     else
3240       rateAlphaUpQUp = rate;
3241
3242     if (iOmega + 1 < Nomega && iQ + 1 < Nq)
3243       rateOmegaUpQUp = dGamma_qq_q->at(positionOmegaUpQUp);
3244     else
3245       rateOmegaUpQUp = rate;
3246
3247     if (iAlphas < Nalphas && iOmega < Nomega && iQ < Nq)
3248       rateAlphaUpOmegaUpQUp = dGamma_qq_q->at(positionAlphaUpOmegaUpQUp);
3249     else
3250       rateAlphaUpOmegaUpQUp = rate;
3251
3252     // used for extrapolation when the data points are too far apart

3253     if (omega > 20.) {
3254       if (iQ + 2 < Nq)
3255         rate2QUp = dGamma_qq_q->at(position2QUp);
3256       else
3257         rate2QUp = rateQUp;
3258
3259       if (iAlphas < Nalphas && iQ + 2 < Nq)
3260         rateAlphaUp2QUp = dGamma_qq_q->at(positionAlphaUpQUp + 1);
3261       else
3262         rateAlphaUp2QUp = rateAlphaUpQUp;
3263
3264       if (iOmega < Nomega && iQ + 2 < Nq)
3265         rateOmegaUp2QUp = dGamma_qq_q->at(positionOmegaUpQUp + 1);
3266       else
3267         rateOmegaUp2QUp = rateOmegaUpQUp;
3268
3269       if (iAlphas < Nalphas && iOmega < Nomega && iQ + 2 < Nq)
3270         rateAlphaUpOmegaUp2QUp = dGamma_qq_q->at(positionAlphaUpOmegaUpQUp + 1);
3271       else
3272         rateAlphaUpOmegaUp2QUp = rateAlphaUpOmegaUpQUp;
3273     }
3274   } else {
3275     if (position > 0 && iAlphas < Nalphas && iOmega < Nomega && iQ < Nq)
3276       rate = dGamma_qg_q->at(position);
3277     else
3278       rate = 0.;
3279
3280     if (iAlphas + 1 < Nalphas)
3281       rateAlphaUp = dGamma_qg_q->at(positionAlphaUp);
3282     else
3283       rateAlphaUp = rate;
3284
3285     if (iOmega + 1 < Nomega)
3286       rateOmegaUp = dGamma_qg_q->at(positionOmegaUp);
3287     else
3288       rateOmegaUp = rate;
3289
3290     if (iQ + 1 < Nq)
3291       rateQUp = dGamma_qg_q->at(positionQUp);
3292     else
3293       rateQUp = rate;
3294
3295     if (iAlphas < Nalphas && iOmega < Nomega)
3296       rateAlphaUpOmegaUp = dGamma_qg_q->at(positionAlphaUpOmegaUp);
3297     else
3298       rateAlphaUpOmegaUp = rate;
3299
3300     if (iAlphas < Nalphas && iQ < Nq)
3301       rateAlphaUpQUp = dGamma_qg_q->at(positionAlphaUpQUp);
3302     else
3303       rateAlphaUpQUp = rate;
3304
3305     if (iOmega + 1 < Nomega && iQ + 1 < Nq)
3306       rateOmegaUpQUp = dGamma_qg_q->at(positionOmegaUpQUp);
3307     else
3308       rateOmegaUpQUp = rate;
3309
3310     if (iAlphas < Nalphas && iOmega < Nomega && iQ < Nq)
3311       rateAlphaUpOmegaUpQUp = dGamma_qg_q->at(positionAlphaUpOmegaUpQUp);
3312     else
3313       rateAlphaUpOmegaUpQUp = rate;
3314
3315     // used for extrapolation when the data points are too far apart

3316     if (omega > 20.) {
3317       if (iQ + 2 < Nq)
3318         rate2QUp = dGamma_qg_q->at(position2QUp);
3319       else
3320         rate2QUp = rateQUp;
3321
3322       if (iAlphas < Nalphas && iQ + 2 < Nq)
3323         rateAlphaUp2QUp = dGamma_qg_q->at(positionAlphaUpQUp + 1);
3324       else
3325         rateAlphaUp2QUp = rateAlphaUpQUp;
3326
3327       if (iOmega < Nomega && iQ + 2 < Nq)
3328         rateOmegaUp2QUp = dGamma_qg_q->at(positionOmegaUpQUp + 1);
3329       else
3330         rateOmegaUp2QUp = rateOmegaUpQUp;
3331
3332       if (iAlphas < Nalphas && iOmega < Nomega && iQ + 2 < Nq)
3333         rateAlphaUpOmegaUp2QUp = dGamma_qg_q->at(positionAlphaUpOmegaUpQUp + 1);
3334       else
3335         rateAlphaUpOmegaUp2QUp = rateAlphaUpOmegaUpQUp;
3336     }
3337   }
3338
3339   if (omega > 0. && omega <= 20.) {
3340     rateOmegaAv = (1. - omegaFrac) * rate + omegaFrac * rateOmegaUp;
3341     rateAlphaUpOmegaAv =
3342         (1. - omegaFrac) * rateAlphaUp + omegaFrac * rateAlphaUpOmegaUp;
3343     rateQUpOmegaAv = (1. - omegaFrac) * rateQUp + omegaFrac * rateOmegaUpQUp;
3344     rateAlphaUpQUpOmegaAv =
3345         (1. - omegaFrac) * rateAlphaUpQUp + omegaFrac * rateAlphaUpOmegaUpQUp;
3346   } else if (omega > 20.) {
3347     if (rate != 0. && rateOmegaUp != 0.)
3348       rateOmegaAv =
3349           exp((1. - omegaFrac) * log(rate) + omegaFrac * log(rateOmegaUp));
3350     else if (rate == 0.)
3351       rateOmegaAv = rateOmegaUp;
3352     else if (rateOmegaUp == 0.)
3353       rateOmegaAv = rate;
3354     else
3355       rateOmegaAv = 0.;
3356
3357     if (rateAlphaUpOmegaUp != 0. && rateAlphaUp != 0.)
3358       rateAlphaUpOmegaAv = exp((1. - omegaFrac) * log(rateAlphaUp) +
3359                                omegaFrac * log(rateAlphaUpOmegaUp));
3360     else if (rateAlphaUp == 0.)
3361       rateAlphaUpOmegaAv = rateAlphaUpOmegaUp;
3362     else if (rateAlphaUpOmegaUp == 0.)
3363       rateAlphaUpOmegaAv = rateAlphaUp;
3364     else
3365       rateAlphaUpOmegaAv = 0.;
3366
3367     if (rateOmegaUpQUp != 0. && rateQUp != 0.)
3368       rateQUpOmegaAv = exp((1. - omegaFrac) * log(rateQUp) +
3369                            omegaFrac * log(rateOmegaUpQUp));
3370     else if (rateOmegaUpQUp == 0.)
3371       rateQUpOmegaAv = rateQUp;
3372     else if (rateQUp == 0.)
3373       rateQUpOmegaAv = rateOmegaUpQUp;
3374     else
3375       rateQUpOmegaAv = 0.;
3376
3377     if (rateAlphaUpOmegaUpQUp != 0. && rateAlphaUpQUp != 0.)
3378       rateAlphaUpQUpOmegaAv = exp((1. - omegaFrac) * log(rateAlphaUpQUp) +
3379                                   omegaFrac * log(rateAlphaUpOmegaUpQUp));
3380     else if (rateAlphaUpQUp == 0.)
3381       rateAlphaUpQUpOmegaAv = rateAlphaUpOmegaUpQUp;
3382     else if (rateAlphaUpOmegaUpQUp == 0.)
3383       rateAlphaUpQUpOmegaAv = rateAlphaUpQUp;
3384     else
3385       rateAlphaUpQUpOmegaAv = 0.;
3386
3387     rate2QUpOmegaAv = exp((1. - omegaFrac) * log(rate2QUp) +
3388                           omegaFrac * log(rateOmegaUp2QUp));
3389     rateAlphaUp2QUpOmegaAv = exp((1. - omegaFrac) * log(rateAlphaUp2QUp) +
3390                                  omegaFrac * log(rateAlphaUpOmegaUp2QUp));
3391   } else if (omega < 0.) {
3392     rateOmegaAv = (omegaFrac)*rate + (1. - omegaFrac) * rateOmegaUp;
3393     rateQUpOmegaAv = (omegaFrac)*rateQUp + (1. - omegaFrac) * rateOmegaUpQUp;
3394     rateAlphaUpOmegaAv =
3395         (omegaFrac)*rateAlphaUp + (1. - omegaFrac) * rateAlphaUpOmegaUp;
3396     rateAlphaUpQUpOmegaAv =
3397         (omegaFrac)*rateAlphaUpQUp + (1. - omegaFrac) * rateAlphaUpOmegaUpQUp;
3398   }
3399
3400   // interpolate logs for large omega

3401   if (omega > 20.) {
3402     if (rateOmegaAv > 0.) {
3403       rateQAv =
3404           exp((1. - qFrac) * log(rateOmegaAv) + qFrac * log(rateQUpOmegaAv));
3405     } else if (rateOmegaAv < 0.) // use extrapolation

3406     {
3407       slope = (log(rate2QUpOmegaAv) - log(rateQUpOmegaAv)) / qStep;
3408       rateQAv = exp(log(rateQUpOmegaAv) - slope * ((1. - qFrac) * qStep));
3409     } else {
3410       rateQAv = 0.;
3411     }
3412
3413     if (rateAlphaUpOmegaAv > 0.) {
3414       rateAlphaUpQAv = exp((1. - qFrac) * log(rateAlphaUpOmegaAv) +
3415                            qFrac * log(rateAlphaUpQUpOmegaAv));
3416     } else if (rateAlphaUpOmegaAv < 0.) // use extrapolation

3417     {
3418       slopeAlphaUp =
3419           (log(rateAlphaUp2QUpOmegaAv) - log(rateAlphaUpQUpOmegaAv)) / qStep;
3420       rateAlphaUpQAv = exp(log(rateAlphaUpQUpOmegaAv) -
3421                            slopeAlphaUp * ((1. - qFrac) * qStep));
3422     } else {
3423       rateAlphaUpQAv = 0.;
3424     }
3425   }
3426   // interpolate linearly for small omega

3427   else {
3428     rateQAv = (1. - qFrac) * rateOmegaAv + qFrac * rateQUpOmegaAv;
3429     rateAlphaUpQAv =
3430         (1. - qFrac) * rateAlphaUpOmegaAv + qFrac * rateAlphaUpQUpOmegaAv;
3431   }
3432
3433   result = (1. - alphaFrac) * rateQAv + alphaFrac * rateAlphaUpQAv;
3434
3435   // the absolute normalization doesn't matter when sampling the shape

3436   // it only matters in "totalRate" etc.

3437   return result;
3438 }
3439
3440 double LambertW(double z) {
3441   double w_new, w_old, ratio, e_old, tol;
3442   int n;
3443
3444   tol = 1.0e-14;
3445
3446   if (z <= -exp(-1.0)) {
3447     JSWARN << "LambertW is not defined for z = " << z;
3448     JSWARN << "z needs to be bigger than " << -exp(-1.0);
3449     throw std::runtime_error("LambertW small z problem");
3450   }
3451
3452   if (z > 3.0) {
3453     w_old = log(z) - log(log(z));
3454   } else {
3455     w_old = 1.0;
3456   }
3457
3458   w_new = 0.0;
3459   ratio = 1.0;
3460   n = 0;
3461   while (std::abs(ratio) > tol) {
3462     e_old = exp(w_old);
3463     ratio = w_old * e_old - z;
3464     ratio /= (e_old * (w_old + 1.0) -
3465               (w_old + 2.0) * (w_old * e_old - z) / (2.0 * w_old + 2.0));
3466     w_new = w_old - ratio;
3467     w_old = w_new;
3468     n++;
3469     if (n > 99) {
3470       JSWARN << "LambertW is not converging after 100 iterations.";
3471       JSWARN << "LambertW: z = " << z;
3472       JSWARN << "LambertW: w_old = " << w_old;
3473       JSWARN << "LambertW: w_new = " << w_new;
3474       JSWARN << "LambertW: ratio = " << ratio;
3475       throw std::runtime_error("LambertW not conversing");
3476     }
3477   }
3478
3479   return w_new;
3480 } // LambertW by Sangyong Jeon