Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-03 08:20:05 

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 "LBT.h"
0017 #include "JetScapeLogger.h"
0018 #include "JetScapeXML.h"
0019 
0020 #include "tinyxml2.h"
0021 #include <string>
0022 #include <sstream>
0023 #include <fstream>
0024 #include <iostream>
0025 #include <iomanip>
0026 
0027 #include "FluidDynamics.h"
0028 #include "LBTMutex.h"
0029 #define MAGENTA "\033[35m"
0030 
0031 using namespace Jetscape;
0032 using namespace std;
0033 
0034 // Register the module with the base class
0035 RegisterJetScapeModule<LBT> LBT::reg("Lbt");
0036 
0037 // initialize static members
0038 bool LBT::flag_init = 0;
0039 double LBT::Rg[60][20] = {
0040     {0.0}}; //total gluon scattering rate as functions of initial energy and temperature
0041 double LBT::Rg1[60][20] = {{0.0}}; //gg-gg              CT1
0042 double LBT::Rg2[60][20] = {{0.0}}; //gg-qqbar           CT2
0043 double LBT::Rg3[60][20] = {{0.0}}; //gq-qg              CT3
0044 double LBT::Rq[60][20] = {
0045     {0.0}}; //total gluon scattering rate as functions of initial energy and temperature
0046 double LBT::Rq3[60][20] = {{0.0}}; //qg-qg              CT13
0047 double LBT::Rq4[60][20] = {{0.0}}; //qiqj-qiqj          CT4
0048 double LBT::Rq5[60][20] = {{0.0}}; //qiqi-qiqi          CT5
0049 double LBT::Rq6[60][20] = {{0.0}}; //qiqibar-qjqjbar    CT6
0050 double LBT::Rq7[60][20] = {{0.0}}; //qiqibar-qiqibar    CT7
0051 double LBT::Rq8[60][20] = {{0.0}}; //qqbar-gg           CT8
0052 double LBT::qhatLQ[60][20] = {{0.0}};
0053 double LBT::qhatG[60][20] = {{0.0}};
0054 
0055 double LBT::RHQ[60][20] = {{0.0}};    //total scattering rate for heavy quark
0056 double LBT::RHQ11[60][20] = {{0.0}};  //Qq->Qq
0057 double LBT::RHQ12[60][20] = {{0.0}};  //Qg->Qg
0058 double LBT::qhatHQ[60][20] = {{0.0}}; //qhat of heavy quark
0059 
0060 double LBT::dNg_over_dt_c[t_gn + 2][temp_gn + 1][HQener_gn + 1] = {{{0.0}}};
0061 double LBT::dNg_over_dt_q[t_gn + 2][temp_gn + 1][HQener_gn + 1] = {{{0.0}}};
0062 double LBT::dNg_over_dt_g[t_gn + 2][temp_gn + 1][HQener_gn + 1] = {{{0.0}}};
0063 double LBT::max_dNgfnc_c[t_gn + 2][temp_gn + 1][HQener_gn + 1] = {{{0.0}}};
0064 double LBT::max_dNgfnc_q[t_gn + 2][temp_gn + 1][HQener_gn + 1] = {{{0.0}}};
0065 double LBT::max_dNgfnc_g[t_gn + 2][temp_gn + 1][HQener_gn + 1] = {{{0.0}}};
0066 
0067 double LBT::initMCX[maxMC] = {0.0};
0068 double LBT::initMCY[maxMC] = {0.0};
0069 double LBT::distFncB[N_T][N_p1][N_e2] = {{{0.0}}};
0070 double LBT::distFncF[N_T][N_p1][N_e2] = {{{0.0}}};
0071 double LBT::distMaxB[N_T][N_p1][N_e2] = {{{0.0}}};
0072 double LBT::distMaxF[N_T][N_p1][N_e2] = {{{0.0}}};
0073 double LBT::distFncBM[N_T][N_p1] = {{0.0}};
0074 double LBT::distFncFM[N_T][N_p1] = {{0.0}};
0075 
0076 LBT::LBT() {
0077   SetId("LBT");
0078   VERBOSE(8);
0079 
0080   // create and set Martini Mutex
0081   auto lbt_mutex = make_shared<LBTMutex>();
0082   SetMutex(lbt_mutex);
0083 }
0084 
0085 LBT::~LBT() { VERBOSE(8); }
0086 
0087 void LBT::Init() {
0088   JSINFO << "Initialize LBT ...";
0089 
0090   //...Below is added by Shanshan
0091   //...read parameters from LBT.input first, but can be changed in JETSCAPE xml file if any conflict exists
0092 
0093   setParameter(
0094       "LBT-tables/LBT.input"); // re-set parameters from input parameter file
0095 
0096   //    if(checkParameter(argc)==0) { // check whether the input parameters are all correct
0097   //        cout << "Parameter check passed" << endl;
0098   //    } else {
0099   //        cout << "Parameter check failed" << endl;
0100   //        exit(EXIT_FAILURE);
0101   //    }
0102 
0103   std::string s = GetXMLElementText({"Eloss", "Lbt", "name"});
0104   JSDEBUG << s << " to be initialized ...";
0105 
0106   //double m_qhat=-99.99;
0107   //lbt->FirstChildElement("qhat")->QueryDoubleText(&m_qhat);
0108   //SetQhat(m_qhat);
0109   //
0110   //JSDEBUG  << s << " with qhat = "<<GetQhat();
0111 
0112   int in_vac = GetXMLElementInt({"Eloss", "Lbt", "in_vac"});
0113   if (in_vac == 1) {
0114     vacORmed = 0;
0115     fixPosition = 1;
0116   } else {
0117     vacORmed = 1;
0118   }
0119 
0120   Kprimary = GetXMLElementInt({"Eloss", "Lbt", "only_leading"});
0121   run_alphas = GetXMLElementInt({"Eloss", "Lbt", "run_alphas"});
0122   Q00 = GetXMLElementDouble({"Eloss", "Lbt", "Q0"});
0123   fixAlphas = GetXMLElementDouble({"Eloss", "Lbt", "alphas"});
0124   hydro_Tc = GetXMLElementDouble({"Eloss", "Lbt", "hydro_Tc"});
0125   tStart = GetXMLElementDouble({"Eloss", "tStart"});
0126   JSINFO << MAGENTA << "LBT parameters -- in_med: " << vacORmed
0127          << " Q0: " << Q00 << "  only_leading: " << Kprimary
0128          << "  alpha_s: " << fixAlphas << "  hydro_Tc: " << hydro_Tc<<", tStart="<<tStart;
0129 
0130   if (!flag_init) {
0131     read_tables(); // initialize various tables
0132     flag_init = true;
0133   }
0134 
0135   //...define derived quantities
0136   temp00 = temp0;
0137   dt = dtau;
0138   timend = tauend;
0139   time0 = tau0;
0140   //...alphas
0141   alphas = alphas0(Kalphas, temp0);
0142   //...Debye Mass square
0143   qhat0 = DebyeMass2(Kqhat0, alphas, temp0);
0144 
0145   ZeroOneDistribution = uniform_real_distribution<double>{0.0, 1.0};
0146 }
0147 
0148 void LBT::WriteTask(weak_ptr<JetScapeWriter> w) {
0149   VERBOSE(8);
0150   auto f = w.lock();
0151   if (!f)
0152     return;
0153   f->WriteComment("ElossModule Parton List: " + GetId());
0154   f->WriteComment("Energy loss to be implemented accordingly ...");
0155 }
0156 
0157 void LBT::DoEnergyLoss(double deltaT, double time, double Q2,
0158                        vector<Parton> &pIn, vector<Parton> &pOut) {
0159 
0160   double z = 0.5;
0161 
0162   //  if (Q2>5)
0163   //    {
0164   VERBOSESHOWER(8) << MAGENTA << "SentInPartons Signal received : " << deltaT
0165                    << " " << Q2 << " " << &pIn;
0166 
0167   // hydro information should be moved into the particle list loop
0168 
0169   ////FluidCellInfo* check_fluid_info_ptr = new FluidCellInfo;
0170   //std::unique_ptr<FluidCellInfo> check_fluid_info_ptr;
0171   //GetHydroCellSignal(1, 1.0, 1.0, 0.0, check_fluid_info_ptr);
0172   //VERBOSE(8)<< MAGENTA<<"Temperature from Brick (Signal) = "<<check_fluid_info_ptr->temperature;
0173 
0174   double rNum;
0175   int par_status;
0176 
0177   //DEBUG:
0178   //cout<<" ---> "<<pIn.size()<<endl;
0179   for (int i = 0; i < pIn.size(); i++) {
0180 
0181     // Reject photons
0182 
0183     if (pIn[i].pid() == photonid) {
0184       if(pIn[i].pstat() != 22) {
0185           pIn[i].set_stat(22); //Add status code 22 for photons that pass through LBT if it is already not assigned by one of the other JEL modules.
0186               pOut.push_back(pIn[i]);
0187       }
0188       return;
0189     }
0190 
0191     // pass particle infomation to LBT array (only pass one particle each time)
0192 
0193     jetClean();
0194 
0195     int j = 1;
0196 
0197     //// unit test for LBT, MC vs. semi-analytical
0198     //          KATT1[j] = Kjet;
0199     //          P[1][j] = ener;
0200     //          P[2][j] = 0.0;
0201     //          P[3][j] = 0.0;
0202 
0203     //          cout << "check read in: " << pIn[i].pid() << "  " << pIn[i].p_in().x() << "  " << pIn[i].p_in().y() << "  " << pIn[i].p_in().z() << "  " << pIn[i].p_in().t() << "  " << pIn[i].pt() << endl;
0204 
0205     par_status = pIn[i].pstat();
0206 
0207     if (par_status != -1) { // a positive (active) particle
0208 
0209       KATT1[j] = pIn[i].pid();
0210       P[1][j] = pIn[i].p(1);
0211       P[2][j] = pIn[i].p(2);
0212       P[3][j] = pIn[i].p(3);
0213       P[4][j] = amss;
0214       if (std::abs(pIn[i].pid()) == 4 || std::abs(pIn[i].pid()) == 5)
0215         P[4][j] = pIn[i].restmass();
0216       P[0][j] = sqrt(P[1][j] * P[1][j] + P[2][j] * P[2][j] + P[3][j] * P[3][j] +
0217                      P[4][j] * P[4][j]);
0218       P[5][j] = sqrt(P[1][j] * P[1][j] + P[2][j] * P[2][j]);
0219       P[6][j] = pIn[i].e() * pIn[i].e() - P[0][j] * P[0][j]; // virtuality^2
0220       if (P[6][j] < 0.0)
0221         P[6][j] = 0.0;
0222       //if(std::isnan(P[6][j]) || std::isinf(P[6][j]))
0223       //{JSWARN << "first instance:j=" << j << ", P[0][j]=" << P[0][j] << ", P[1][j]=" << P[1][j] << ", P[2][j]=" << P[2][j] <<", P[3][j]=" << P[3][j] <<", P[4][j]=" << P[4][j] << ", P[6][j]=" << P[6][j];}
0224 
0225       WT[j] = 1.0;
0226 
0227       Vfrozen[1][j] = pIn[i].x_in().x();
0228       Vfrozen[2][j] = pIn[i].x_in().y();
0229       Vfrozen[3][j] = pIn[i].x_in().z();
0230       Vfrozen[0][j] = pIn[i].x_in().t();
0231       V[1][j] = Vfrozen[1][j];
0232       V[2][j] = Vfrozen[2][j];
0233       V[3][j] = Vfrozen[3][j];
0234       V[0][j] = -log(1.0 - ZeroOneDistribution(*GetMt19937Generator()));
0235 
0236       for (int k = 0; k <= 3; k++)
0237         Prad[k][j] = P[k][j];
0238 
0239       Tfrozen[j] = 0.0;
0240       vcfrozen[1][j] = 0.0;
0241       vcfrozen[2][j] = 0.0;
0242       vcfrozen[3][j] = 0.0;
0243 
0244       Tint_lrf[j] =
0245           0.0; // time since last splitting, reset below if there is information from JETSCAPE
0246       if (pIn[i].has_user_info<LBTUserInfo>()) {
0247         Tint_lrf[j] = pIn[i].user_info<LBTUserInfo>().lrf_T_tot();
0248       } else {
0249         pIn[i].set_user_info(new LBTUserInfo(0.0));
0250       }
0251 
0252       if (par_status == 1)
0253         CAT[j] = 2;
0254       else
0255         CAT[j] = 0;
0256 
0257     } else { // negative particle, or particle no longer active, will streamly freely in LBT
0258 
0259       KATT10[j] = pIn[i].pid();
0260       P0[1][j] = pIn[i].p(1);
0261       P0[2][j] = pIn[i].p(2);
0262       P0[3][j] = pIn[i].p(3);
0263       P0[4][j] = amss;
0264       P0[0][j] = sqrt(P0[1][j] * P0[1][j] + P0[2][j] * P0[2][j] +
0265                       P0[3][j] * P0[3][j] + P0[4][j] * P0[4][j]);
0266       P0[5][j] = sqrt(P0[1][j] * P0[1][j] + P0[2][j] * P0[2][j]);
0267       P0[6][j] = pIn[i].e() * pIn[i].e() - P0[0][j] * P0[0][j]; // virtuality^2
0268       if (P0[6][j] < 0.0)
0269         P0[6][j] = 0.0;
0270       WT0[j] = 1.0;
0271 
0272       Vfrozen0[1][j] = pIn[i].x_in().x();
0273       Vfrozen0[2][j] = pIn[i].x_in().y();
0274       Vfrozen0[3][j] = pIn[i].x_in().z();
0275       Vfrozen0[0][j] = pIn[i].x_in().t();
0276       V0[1][j] = Vfrozen0[1][j];
0277       V0[2][j] = Vfrozen0[2][j];
0278       V0[3][j] = Vfrozen0[3][j];
0279       V0[0][j] = -log(1.0 - ZeroOneDistribution(*GetMt19937Generator()));
0280 
0281       // for(int k=0;k<=3;k++) Prad[k][j] = P[k][j];
0282 
0283       Tfrozen0[j] = 0.0;
0284       vcfrozen0[1][j] = 0.0;
0285       vcfrozen0[2][j] = 0.0;
0286       vcfrozen0[3][j] = 0.0;
0287 
0288       Tint_lrf[j] =
0289           0.0; // time since last splitting, reset below if there is information from JETSCAPE
0290       if (pIn[i].has_user_info<LBTUserInfo>()) {
0291         Tint_lrf[j] = pIn[i].user_info<LBTUserInfo>().lrf_T_tot();
0292       } else {
0293         pIn[i].set_user_info(new LBTUserInfo(0.0));
0294       }
0295 
0296     } // end par_status check
0297 
0298     ////if(par_status != -1) {
0299     ////    cout << "LBT -- status: " << par_status << " current time: " << Vfrozen[0][j] << "  accumulated time: " << Tint_lrf[j] << "  energy: " << P[0][1] << "  clock: " << time <<endl;
0300     ////} else {
0301     ////    cout << "LBT -- status: " << par_status << " current time: " << Vfrozen0[0][j] << "  accumulated time: " << Tint_lrf[j] << "  energy: " << P0[0][1] << "  clock: " << time << endl;
0302     ////}
0303 
0304     nj = 1;
0305     np = 1;
0306 
0307     dtau = deltaT;
0308     dt = dtau;
0309 
0310     flag_update = 0;
0311     flag_update0 = 0;
0312 
0313     // do LBT evolution for this particle
0314     int nnn = 1; // dummy variable, not important
0315     //          double systemTime=Vfrozen[0][j]+dt; // how to pass the time of the framework?
0316     double systemTime = time;
0317 
0318     //// for unit test
0319     //          if(systemTime-dtau<0.00001) {
0320     ////              cout << "reset counter ........" << endl;
0321     //              dEel=0.0;
0322     //              nGluon=0.0;
0323     //              eGluon=0.0;
0324     //              ctEvt++;
0325     //          }
0326 
0327     //if(alphas>epsilon && vacORmed!=0) {
0328     if (vacORmed != 0) { // for JETSCAPE, won't change parton if it's in vacuum
0329       flagScatter = 0;
0330       LBT0(
0331           nnn,
0332           systemTime); // most important function of LBT: update particle list for one time step
0333     }
0334 
0335     //// LBT unit test
0336     //          dEel=dEel+(ener-P[0][1]);
0337     ////          dEel=1.0;
0338     //          for(int ik=1; ik<=40; ++ik) {
0339     //              if(fabs(systemTime-ik*0.1) < 0.000001) {
0340     //                  de1[ik] = de1[ik] + dEel/ncall;
0341     //                  de2[ik] = de2[ik] + eGluon/ncall;
0342     //              }
0343     //          }
0344     //          // write out energy loss infomation for LBT unit test
0345     //          if(ctEvt==ncall && fabs(systemTime-tauend)<0.000001) {
0346     //              ofstream dat_de("de.dat");
0347     //              for(int ik=1; ik<=40; ++ik)
0348     //                   dat_de << ik*0.1 << "    " << de1[ik] << "    " << de2[ik] << endl;
0349     //              dat_de.close();
0350     //          }
0351 
0352     double velocity_jet[4];
0353     velocity_jet[0] = 1.0;
0354     velocity_jet[1] = pIn[i].jet_v().x();
0355     velocity_jet[2] = pIn[i].jet_v().y();
0356     velocity_jet[3] = pIn[i].jet_v().z();
0357 
0358     if (par_status !=
0359         -1) { // for particle list initiating with positive particles
0360 
0361       if (flag_update == 0)
0362         continue; // no update on this particle from LBT0 function
0363 
0364       TakeResponsibilityFor(
0365           pIn[i]); // Generate error if another module already has responsibility.
0366 
0367       ////cout << "Evolve in LBT -- flag: " << flagScatter << endl;
0368 
0369       // pass particle list back to JETSCAPE
0370       // how to handle negative particle in the framework?
0371       // only pass positive at this moment
0372 
0373       // may only need to push back particle if flagScatter=1 later
0374       // positive particle stat = 0
0375       for (int j = 1; j <= np; j++) {
0376         // definition Parton (int label, int id, int stat, double p[4], double x[4]);
0377         if (P[0][j] < cutOut)
0378           continue;
0379         int out_stat = 0; // jet parton
0380         if (CAT[j] == 2)
0381           out_stat = 1; // recoil parton
0382         double tempP[4], tempX[4];
0383         tempP[0] = sqrt(P[0][j] * P[0][j] + P[6][j]);
0384         tempP[1] = P[1][j];
0385         tempP[2] = P[2][j];
0386         tempP[3] = P[3][j];
0387         tempX[0] = Vfrozen[0][j];
0388         tempX[1] = Vfrozen[1][j];
0389         tempX[2] = Vfrozen[2][j];
0390         tempX[3] = Vfrozen[3][j];
0391         //        pOut.push_back(Parton(0,21,0,newPt,pIn[i].eta(),pIn[i].phi(),newPt));
0392 
0393         if (std::isnan(tempP[1])) {
0394           JSWARN << "second instance: j=" << j << ", P[0][j]=" << P[0][j]
0395                  << ", P[1][j]=" << P[1][j] << ", P[2][j]=" << P[2][j]
0396                  << ", P[3][j]=" << P[3][j] << ", P[4][j]=" << P[4][j]
0397                  << ", P[6][j]=" << P[6][j];
0398         }
0399 
0400         pOut.push_back(Parton(0, KATT1[j], out_stat, tempP, tempX));
0401         // remember to put Tint_lrf infomation back to JETSCAPE
0402         pOut.back().set_user_info(new LBTUserInfo(Tint_lrf[j]));
0403 
0404         int iout = pOut.size() - 1;
0405         pOut[iout].set_jet_v(velocity_jet); // use initial jet velocity
0406         pOut[iout].set_mean_form_time();
0407         double ft = pOut[iout].mean_form_time();
0408         pOut[iout].set_form_time(ft);
0409       }
0410 
0411       // negative particle stat = -1
0412       for (int j = 2; j <= np; j++) {
0413         if (P0[0][j] < cutOut)
0414           continue;
0415         double tempP[4], tempX[4];
0416         tempP[0] = sqrt(P0[0][j] * P0[0][j] + P0[6][j]);
0417         tempP[1] = P0[1][j];
0418         tempP[2] = P0[2][j];
0419         tempP[3] = P0[3][j];
0420         tempX[0] = Vfrozen0[0][j];
0421         tempX[1] = Vfrozen0[1][j];
0422         tempX[2] = Vfrozen0[2][j];
0423         tempX[3] = Vfrozen0[3][j];
0424         //        pOut.push_back(Parton(0,21,0,newPt,pIn[i].eta(),pIn[i].phi(),newPt));
0425         pOut.push_back(Parton(0, KATT10[j], -1, tempP, tempX));
0426         pOut.back().set_user_info(new LBTUserInfo(0.0));
0427 
0428         int iout = pOut.size() - 1;
0429         pOut[iout].set_jet_v(velocity_jet); // use initial jet velocity
0430         pOut[iout].set_mean_form_time();
0431         double ft = pOut[iout].mean_form_time();
0432         pOut[iout].set_form_time(ft);
0433       }
0434 
0435     } else { // free-streaming daughter particles from negative particles if they are inside a medium
0436 
0437       if (flag_update0 == 0)
0438         continue; // no update on this particle from LBT0 function
0439 
0440       TakeResponsibilityFor(
0441           pIn[i]); // Generate error if another module already has responsibility.
0442 
0443       if (np != 1)
0444         JSDEBUG << "Wrong number of negative partons!";
0445       ////cout << "Evolve in LBT -- flag: " << flagScatter << endl;
0446 
0447       for (int j = 1; j <= np; j++) {
0448         if (P0[0][j] < cutOut)
0449           continue;
0450         double tempP[4], tempX[4];
0451         tempP[0] = P0[0][j];
0452         tempP[1] = P0[1][j];
0453         tempP[2] = P0[2][j];
0454         tempP[3] = P0[3][j];
0455         tempX[0] = Vfrozen0[0][j];
0456         tempX[1] = Vfrozen0[1][j];
0457         tempX[2] = Vfrozen0[2][j];
0458         tempX[3] = Vfrozen0[3][j];
0459         //        pOut.push_back(Parton(0,21,0,newPt,pIn[i].eta(),pIn[i].phi(),newPt));
0460         pOut.push_back(Parton(0, KATT10[j], -1, tempP, tempX));
0461         pOut.back().set_user_info(new LBTUserInfo(0.0));
0462 
0463         int iout = pOut.size() - 1;
0464         pOut[iout].set_jet_v(velocity_jet); // use initial jet velocity
0465         pOut[iout].set_mean_form_time();
0466         double ft = pOut[iout].mean_form_time();
0467         pOut[iout].set_form_time(ft);
0468       }
0469 
0470     } // end par_status check
0471   }
0472 }
0473 
0474 ///////////////////////////////////////////////////////////////////////////////////////////////////
0475 //
0476 // This is the key function of LBT
0477 // Update elastic and inelastic evolution of the particle list for one time step
0478 //
0479 ///////////////////////////////////////////////////////////////////////////////////////////////////
0480 
0481 void LBT::LBT0(int &n, double &ti) {
0482 
0483   int CT;     //collision type 1 2 3 13 4 5 6 7 8
0484   int KATTC0; //flavor code 1/d 2/u 3/s -1/dbar -2/ubar -3/sbar 21/gluon
0485   int KATT2;
0486   int KATT3;
0487   int KATT30;
0488 
0489   double RTE;  //scattering rate (energy temperature)
0490   double E;    //parton energy
0491   double PLen; //parton momentum
0492   double T;    //local temperature
0493 
0494   double T1; //index for difference method
0495   double T2;
0496   double E1;
0497   double E2;
0498   int iT1;
0499   int iT2;
0500   int iE1;
0501   int iE2;
0502 
0503   int nrad;
0504   int idlead;
0505   int idlead1;
0506   int idlead2;
0507   double Xtau; //....................main & titau
0508   double Vx;
0509   double Vy;
0510   double Veta;
0511 
0512   double tcar =
0513       0.0; //....................t x y z in Cartesian coordinate this is for the radiation process
0514   double xcar = 0.0;
0515   double ycar = 0.0;
0516   double zcar = 0.0;
0517 
0518   double tcar0 =
0519       0.0; //....................t x y z in Cartesian coordinate this is for the radiation process
0520   double xcar0 = 0.0;
0521   double ycar0 = 0.0;
0522   double zcar0 = 0.0;
0523 
0524   double rans;
0525 
0526   int KATTx;
0527   int KATTy;
0528 
0529   double Ejp;
0530   double Elab;
0531   double Qinit;
0532 
0533   double qt;
0534 
0535   double px0;
0536   double py0;
0537   double pz0;
0538 
0539   double
0540       Ncoll22; //...................average number of elastic scattering for particle i in the current step
0541   int Nscatter; //...................number of elastic scattering for particle i in the current step
0542 
0543   int nnpp = np; //...................number of particles in the current np loop
0544   int np0 =
0545       np; //...................number of particles in the current step (np0)
0546   //...................number of particles in the beginning of current step (np)
0547 
0548   int free = 0;
0549   int free0 = 0;
0550   double fraction = 0.0;
0551   double fraction0 = 0.0;
0552   double vc0b[4] = {0.0}; //flow velocity
0553   double pMag, vMag, flowFactor;
0554 
0555   double kt2;
0556 
0557   double vp0[4] = {0.0};
0558   double p0temp[4] = {0.0};
0559 
0560   double p0temp1 = 0.0;
0561   double p0temp2 = 0.0;
0562 
0563   double pcx[4] = {0.0};
0564   double pcy[4] = {0.0};
0565 
0566   double pcx1[4] = {0.0};
0567   double pcy1[4] = {0.0};
0568 
0569   // for heavy quark
0570   double dt_lrf, maxFncHQ, Tdiff, lim_low, lim_high, lim_int;
0571   double probCol, probRad, probTot;
0572 
0573   double lrf_rStart[4] = {0.0};
0574   double SpatialRapidity = 0.0;
0575 
0576   probCol = 0.0;
0577   probRad = 0.0;
0578   probTot = 0.0;
0579   flowFactor = 0.0;
0580 
0581   //...........................................................................................................NCUT
0582   ncut = 0;
0583   ncut0 = 0;
0584   //...........................................................................................................NCUTEND
0585 
0586   //...collision of the jet parton (i<=nj), recoiled partons, radiated gluons with the background thermal partons
0587   //...thermal partons
0588   //...np number of partons in current step
0589   //...nj number of jet partons
0590 
0591   for (int i = 1; i <= np; i++) {
0592 
0593     //      // if the production time of a parton is ahead of the current time, do nothing
0594     //      if(Vfrozen[0][i]>=ti) continue;
0595 
0596     icl22 = 0;
0597     qt = 0;
0598     idlead = i;
0599 
0600     //......propagation of each positive particle and its negative counterpart in the current time step (for all particles)
0601 
0602     if (P[0][i] < epsilon) { //anti-overflow NOT NECESSARY
0603       V[1][i] = 0.0;
0604       V[2][i] = 0.0;
0605       V[3][i] = 0.0;
0606       CAT[i] = 1;
0607     }
0608 
0609     if (P0[0][i] < epsilon) { //anti-overflow NOT NECESSARY
0610       V0[1][i] = 0.0;
0611       V0[2][i] = 0.0;
0612       V0[3][i] = 0.0;
0613       CAT0[i] = 1;
0614     }
0615 
0616     if (CAT0[i] != 1) {
0617       //...........................................t-z coordinates
0618       if (tauswitch == 0) {
0619 
0620         //V0[1][i]=V0[1][i]+dt*P0[1][i]/P0[0][i];
0621         //V0[2][i]=V0[2][i]+dt*P0[2][i]/P0[0][i];
0622         //V0[3][i]=V0[3][i]+dt*P0[3][i]/P0[0][i];
0623 
0624         V0[1][i] = V0[1][i] + (ti - Vfrozen0[0][i]) * P0[1][i] / P0[0][i];
0625         V0[2][i] = V0[2][i] + (ti - Vfrozen0[0][i]) * P0[2][i] / P0[0][i];
0626         V0[3][i] = V0[3][i] + (ti - Vfrozen0[0][i]) * P0[3][i] / P0[0][i];
0627 
0628         tcar0 = ti;
0629         zcar0 = V0[3][i];
0630         xcar0 = V0[1][i];
0631         ycar0 = V0[2][i];
0632 
0633         double ed00, sd00, VX00, VY00, VZ00;
0634         int hydro_ctl0;
0635 
0636         free0 = 0;
0637 
0638         //              if(bulkFlag==1) { // read OSU hydro
0639         //                  readhydroinfoshanshan_(&tcar0,&xcar0,&ycar0,&zcar0,&ed00,&sd00,&temp00,&VX00,&VY00,&VZ00,&hydro_ctl0);
0640         //              } else if(bulkFlag==2) { // read CCNU hydro
0641         //                  hydroinfoccnu_(&tcar0, &xcar0, &ycar0, &zcar0, &temp00, &VX00, &VY00, &VZ00, &hydro_ctl0);
0642         //              } else if(bulkFlag==0) { // static medium
0643 
0644         //Extract fluid properties
0645         std::unique_ptr<FluidCellInfo> check_fluid_info_ptr;
0646 
0647         SpatialRapidity = 0.5 * std::log((tcar0 + zcar0) / (tcar0 - zcar0));
0648 
0649         GetHydroCellSignal(tcar0, xcar0, ycar0, zcar0, check_fluid_info_ptr);
0650         //VERBOSE(7)<< MAGENTA<<"Temperature from Brick (Signal) = "<<check_fluid_info_ptr->temperature;
0651 
0652         temp00 = check_fluid_info_ptr->temperature;
0653         sd00 = check_fluid_info_ptr->entropy_density;
0654         VX00 = check_fluid_info_ptr->vx;
0655         VY00 = check_fluid_info_ptr->vy;
0656         VZ00 = check_fluid_info_ptr->vz;
0657 
0658         if (tcar0 < tStart * cosh(SpatialRapidity)) {
0659           temp00 = 0.0;
0660           sd00 = 0.0;
0661           VX00 = 0.0;
0662           VY00 = 0.0;
0663           VZ00 = 0.0;
0664         }
0665 
0666         //JSDEBUG << "check temperature: " << temp00;
0667 
0668         //VX00=0.0;
0669         //VY00=0.0;
0670         //VZ00=0.0;
0671 
0672         hydro_ctl0 = 0;
0673         //              }
0674 
0675         if (hydro_ctl0 == 0 && temp00 >= hydro_Tc) {
0676 
0677           qhat00 = DebyeMass2(Kqhat0, alphas, temp00);
0678           fraction0 = 1.0;
0679 
0680           Vfrozen0[0][i] = ti;
0681           Vfrozen0[1][i] = V0[1][i];
0682           Vfrozen0[2][i] = V0[2][i];
0683           Vfrozen0[3][i] = V0[3][i];
0684           Tfrozen0[i] = temp00;
0685           vcfrozen0[1][i] = VX00;
0686           vcfrozen0[2][i] = VY00;
0687           vcfrozen0[3][i] = VZ00;
0688           flag_update0 =
0689               1; // do free-streaming for negative particle if in medium, no check on virtuality (assume 0)
0690 
0691         } else {
0692           free0 = 1;
0693           fraction0 = 0.0;
0694         }
0695 
0696       } //if(tauswitch==0)
0697 
0698     } //if(CAT0[i] != 1)
0699 
0700     if (CAT[i] != 1) {
0701       //...........................................t-z coordinates
0702       if (tauswitch == 0) {
0703 
0704         //V[1][i]=V[1][i]+dt*P[1][i]/P[0][i];
0705         //V[2][i]=V[2][i]+dt*P[2][i]/P[0][i];
0706         //V[3][i]=V[3][i]+dt*P[3][i]/P[0][i];
0707 
0708         V[1][i] = V[1][i] + (ti - Vfrozen[0][i]) * P[1][i] / P[0][i];
0709         V[2][i] = V[2][i] + (ti - Vfrozen[0][i]) * P[2][i] / P[0][i];
0710         V[3][i] = V[3][i] + (ti - Vfrozen[0][i]) * P[3][i] / P[0][i];
0711 
0712         tcar = ti;
0713         zcar = V[3][i];
0714         xcar = V[1][i];
0715         ycar = V[2][i];
0716 
0717         for (int lrf_i = 0; lrf_i <= 3; lrf_i++)
0718           lrf_rStart[lrf_i] = Vfrozen[lrf_i][i];
0719 
0720         //...........................................Hydro part
0721 
0722         double ed, sd, VX, VY, VZ;
0723         int hydro_ctl;
0724 
0725         free = 0;
0726 
0727         if (free == 0) {
0728 
0729           //                  if(bulkFlag==1) { // read OSU hydro
0730           //                      readhydroinfoshanshan_(&tcar,&xcar,&ycar,&zcar,&ed,&sd,&temp0,&VX,&VY,&VZ,&hydro_ctl);
0731           //                  } else if(bulkFlag==2) { // read CCNU hydro
0732           //                      hydroinfoccnu_(&tcar, &xcar, &ycar, &zcar, &temp0, &VX, &VY, &VZ, &hydro_ctl);
0733           //                  } else if(bulkFlag==0) { // static medium
0734           //VX=0.0;
0735           //VY=0.0;
0736           //VZ=0.0;
0737 
0738           std::unique_ptr<FluidCellInfo> check_fluid_info_ptr;
0739 
0740           SpatialRapidity = 0.5 * std::log((tcar + zcar) / (tcar - zcar));
0741 
0742           GetHydroCellSignal(tcar, xcar, ycar, zcar, check_fluid_info_ptr);
0743           //VERBOSE(7)<< MAGENTA<<"Temperature from Brick (Signal) = "<<check_fluid_info_ptr->temperature;
0744 
0745           if (!GetJetSignalConnected()) {
0746             JSWARN << "Couldn't find a hydro module attached!";
0747             throw std::runtime_error("Please attach a hydro module or set "
0748                                      "in_vac to 1 in the XML file");
0749           }
0750 
0751           temp0 = check_fluid_info_ptr->temperature;
0752           sd = check_fluid_info_ptr->entropy_density;
0753           VX = check_fluid_info_ptr->vx;
0754           VY = check_fluid_info_ptr->vy;
0755           VZ = check_fluid_info_ptr->vz;
0756 
0757           if (tcar < tStart * cosh(SpatialRapidity)) {
0758             temp0 = 0.0;
0759             sd = 0.0;
0760             VX = 0.0;
0761             VY = 0.0;
0762             VZ = 0.0;
0763           }
0764 
0765           // JSINFO << BOLDYELLOW << "LBT time = " << tcar << " x = " << xcar << " y = " << ycar << " z = " << zcar << " temp = " << temp0;
0766           // JSINFO << BOLDYELLOW << "LBT Vel_x, Vel_y, Vel_z =" << P[1][i]/P[0][i] << ", " << P[2][i]/P[0][i]<< ", " << P[3][i]/P[0][i];
0767 
0768           //JSDEBUG << "check temperature: " << temp0;
0769           hydro_ctl = 0;
0770           //                  }
0771 
0772           if (hydro_ctl == 0 && temp0 >= hydro_Tc) {
0773 
0774             vc0[1] = VX;
0775             vc0[2] = VY;
0776             vc0[3] = VZ;
0777 
0778             //                      vc0[1]=0.0;
0779             //                      vc0[2]=0.0;
0780             //                      vc0[3]=0.0;
0781 
0782             vc0b[1] = VX;
0783             vc0b[2] = VY;
0784             vc0b[3] = VZ;
0785 
0786             //...alphas!
0787             alphas = alphas0(Kalphas, temp0);
0788             //...Debye Mass square
0789             qhat0 = DebyeMass2(Kqhat0, alphas, temp0);
0790 
0791             fraction = 1.0;
0792             Vfrozen[0][i] = ti;
0793             Vfrozen[1][i] = V[1][i];
0794             Vfrozen[2][i] = V[2][i];
0795             Vfrozen[3][i] = V[3][i];
0796             Tfrozen[i] = temp0;
0797             vcfrozen[1][i] = VX;
0798             vcfrozen[2][i] = VY;
0799             vcfrozen[3][i] = VZ;
0800 
0801           } else {
0802             free = 1;
0803             fraction = 0.0;
0804           }
0805 
0806           pc0[1] = P[1][i];
0807           pc0[2] = P[2][i];
0808           pc0[3] = P[3][i];
0809           pc0[0] = P[0][i];
0810 
0811           vMag =
0812               sqrt(vc0b[1] * vc0b[1] + vc0b[2] * vc0b[2] + vc0b[3] * vc0b[3]);
0813           flowFactor =
0814               (1.0 - (pc0[1] * vc0b[1] + pc0[2] * vc0b[2] + pc0[3] * vc0b[3]) /
0815                          pc0[0]) /
0816               sqrt(1.0 - vMag * vMag);
0817 
0818         } //if(free==0)
0819       }   //if(tauswitch==0)
0820 
0821       //...........................................tau-eta coordinates
0822       //if(tauswitch==1) { // not ready for JETSCAPE
0823       //} // if(tauswitch==1)
0824 
0825       //......propagation of each positive particle and its negative counterpart in the current time step   end
0826 
0827       //......CAT[i]=1 free streaming particle
0828 
0829       //......free streaming when energy  < Ecut
0830       //......free streaming when energy0 < Ecut0......in the next step
0831 
0832       if (free == 0) {
0833 
0834         double qhatTP, RTE1, RTE2;
0835 
0836         // modification: interplotate rate in l.r.f.
0837         pc0[1] = P[1][i];
0838         pc0[2] = P[2][i];
0839         pc0[3] = P[3][i];
0840         pc0[0] = P[0][i];
0841         trans(vc0, pc0);
0842         E = pc0
0843             [0]; //  p4-the initial 4-momentum of the jet parton in the local rest frame
0844         PLen = sqrt(pc0[1] * pc0[1] + pc0[2] * pc0[2] + pc0[3] * pc0[3]);
0845         transback(vc0, pc0);
0846 
0847         T = temp0;
0848         KATTC0 = KATT1[i];
0849 
0850         //              if(E<T) preKT=4.0*pi/9.0/log(scaleAK*T*T/0.04)/alphas/fixKT;
0851         //              else preKT=4.0*pi/9.0/log(scaleAK*E*T/0.04)/alphas/fixKT;
0852 
0853         if(run_alphas==1) {
0854             //runKT=4.0*pi/9.0/log(2.0*E*T/0.04)/0.3;
0855             fixedLog = log(5.7*E/4.0/6.0/pi/0.3/T);
0856             scaleMu2 = 2.0*E*T;
0857             if(scaleMu2 < 1.0) {
0858                 scaleMu2 = 1.0;
0859                 runAlphas = alphas;
0860             } else {
0861                 double lambdaQCD2 = exp(-4.0*pi/9.0/alphas);
0862                 runAlphas = 4.0*pi/9.0/log(scaleMu2/lambdaQCD2);
0863             }
0864             runKT = runAlphas/0.3;
0865             runLog = log(scaleMu2/6.0/pi/T/T/alphas)/fixedLog;
0866         }    
0867 
0868         lam(KATTC0, RTE, PLen, T, T1, T2, E1, E2, iT1, iT2, iE1,
0869             iE2); //modified: use P instead
0870 
0871         preKT = alphas / 0.3;
0872 
0873         // calculate p,T-dependence K factor
0874         KPfactor = 1.0 + KPamp * exp(-PLen * PLen / 2.0 / KPsig / KPsig);
0875         KTfactor = 1.0 + KTamp * exp(-pow((temp0 - hydro_Tc), 2) / 2.0 / KTsig /
0876                                      KTsig);
0877 
0878         if(run_alphas==1) {
0879            Kfactor = KPfactor * KTfactor * KTfactor * runKT * preKT * runLog; // K factor for qhat
0880         } else {
0881            Kfactor = KPfactor * KTfactor * KTfactor * preKT * preKT; // K factor for qhat
0882         }
0883        
0884 
0885         // get qhat from table
0886         if (KATTC0 == 21) {
0887           RTE1 = (qhatG[iT2][iE1] - qhatG[iT1][iE1]) * (T - T1) / (T2 - T1) +
0888                  qhatG[iT1][iE1];
0889           RTE2 = (qhatG[iT2][iE2] - qhatG[iT1][iE2]) * (T - T1) / (T2 - T1) +
0890                  qhatG[iT1][iE2];
0891         } else if (KATTC0 == 4 || KATTC0 == -4 || KATTC0 == 5 || KATTC0 == -5) {
0892           RTE1 = (qhatHQ[iT2][iE1] - qhatHQ[iT1][iE1]) * (T - T1) / (T2 - T1) +
0893                  qhatHQ[iT1][iE1];
0894           RTE2 = (qhatHQ[iT2][iE2] - qhatHQ[iT1][iE2]) * (T - T1) / (T2 - T1) +
0895                  qhatHQ[iT1][iE2];
0896         } else {
0897           RTE1 = (qhatLQ[iT2][iE1] - qhatLQ[iT1][iE1]) * (T - T1) / (T2 - T1) +
0898                  qhatLQ[iT1][iE1];
0899           RTE2 = (qhatLQ[iT2][iE2] - qhatLQ[iT1][iE2]) * (T - T1) / (T2 - T1) +
0900                  qhatLQ[iT1][iE2];
0901         }
0902 
0903         qhatTP = (RTE2 - RTE1) * (PLen - E1) / (E2 - E1) + RTE1;
0904 
0905         qhatTP = qhatTP * Kfactor;
0906 
0907         ////reset by hand for unit test
0908         ////              RTE=0.09747;
0909         ////              qhatTP=6.469;
0910         ////              RTE=0.2200;
0911         //              qhatTP=14.80;
0912         //              cout << "RTE: " << RTE << endl;
0913 
0914         qhat_over_T3 = qhatTP; // what is read in is qhat/T^3 of quark/gluon
0915 
0916         // D2piT is diffusion coefficient "just for quark"
0917         if (KATTC0 == 21)
0918           D2piT = 8.0 * pi / (qhat_over_T3 / CA * CF);
0919         else
0920           D2piT = 8.0 * pi / qhat_over_T3;
0921 
0922         qhat =
0923             qhatTP *
0924             pow(T,
0925                 3); // for light quark and gluon, need additional table to make it correct
0926 
0927         if (Q00 < 0.0) {
0928           Q0 = sqrt(sqrt(2.0 * E * qhat));
0929         } else {
0930           Q0 = Q00;
0931         }
0932 
0933         if (P[6][i] > Q0 * Q0 + rounding_error) {
0934           continue; // virtuality outside the LBT range, go to the next particle
0935         }
0936 
0937         flag_update =
0938             1; // satisfy T>Tc and Q<Q0, LBT will process this positice particle
0939 
0940         if (E < sqrt(qhat0))
0941           continue; // just do free-streaming
0942         if (i > nj && E < Ecmcut)
0943           continue;
0944 
0945         // update interaction time and average gluon number for heavy quark radiation
0946         // (as long as it's in medium, even though no collision)
0947         dt_lrf =
0948             dt *
0949             flowFactor; // must be put here, flowFactor will be changed below
0950 
0951         // increae time accumulation from the production point of the parton
0952         for (double tLoc = lrf_rStart[0]; tLoc < ti + rounding_error;
0953              tLoc = tLoc + dt) {
0954 
0955           double xLoc =
0956               lrf_rStart[1] + (tLoc - lrf_rStart[0]) * P[1][i] / P[0][i];
0957           double yLoc =
0958               lrf_rStart[2] + (tLoc - lrf_rStart[0]) * P[2][i] / P[0][i];
0959           double zLoc =
0960               lrf_rStart[3] + (tLoc - lrf_rStart[0]) * P[3][i] / P[0][i];
0961           double dtLoc, tempLoc, vxLoc, vyLoc, vzLoc;
0962 
0963           if (tLoc + dt < ti)
0964             dtLoc = dt;
0965           else
0966             dtLoc = ti - tLoc;
0967 
0968           std::unique_ptr<FluidCellInfo> check_fluid_info_ptr;
0969 
0970           SpatialRapidity = 0.5 * std::log((tLoc + zLoc) / (tLoc - zLoc));
0971 
0972           GetHydroCellSignal(tLoc, xLoc, yLoc, zLoc, check_fluid_info_ptr);
0973           tempLoc = check_fluid_info_ptr->temperature;
0974           vxLoc = check_fluid_info_ptr->vx;
0975           vyLoc = check_fluid_info_ptr->vy;
0976           vzLoc = check_fluid_info_ptr->vz;
0977 
0978           if (tLoc < tStart * cosh(SpatialRapidity)) {
0979             tempLoc = 0.0;
0980             vxLoc = 0.0;
0981             vyLoc = 0.0;
0982             vzLoc = 0.0;
0983           }
0984 
0985           if (tempLoc >= hydro_Tc) {
0986             vMag = sqrt(vxLoc * vxLoc + vyLoc * vyLoc + vzLoc * vzLoc);
0987             flowFactor =
0988                 (1.0 -
0989                  (pc0[1] * vxLoc + pc0[2] * vyLoc + pc0[3] * vzLoc) / pc0[0]) /
0990                 sqrt(1.0 - vMag * vMag);
0991             Tint_lrf[i] = Tint_lrf[i] + dtLoc * flowFactor;
0992           }
0993         }
0994 
0995         Tdiff = Tint_lrf[i];
0996 
0997         // get radiation probablity by reading tables -- same for heavy and light partons
0998         if (KATTC0 == 21) {
0999             if (run_alphas==1) {
1000                 radng[i] += nHQgluon(KATT1[i], dt_lrf, Tdiff, temp0, E, maxFncHQ) / 2.0 * KTfactor * runKT;
1001             } else {
1002                 radng[i] += nHQgluon(KATT1[i], dt_lrf, Tdiff, temp0, E, maxFncHQ) / 2.0 * KTfactor * preKT;
1003             }
1004         } else {
1005             if (run_alphas==1) {
1006                 radng[i] += nHQgluon(KATT1[i], dt_lrf, Tdiff, temp0, E, maxFncHQ) * KTfactor * runKT;
1007             } else {
1008                 radng[i] += nHQgluon(KATT1[i], dt_lrf, Tdiff, temp0, E, maxFncHQ) * KTfactor * preKT;
1009             }
1010         }
1011         lim_low = sqrt(6.0 * pi * alphas) * temp0 / E;
1012         if (abs(KATT1[i]) == 4 || abs(KATT1[i]) == 5)
1013           lim_high = 1.0;
1014         else
1015           lim_high = 1.0 - lim_low;
1016         lim_int = lim_high - lim_low;
1017         if (lim_int > 0.0)
1018           probRad = 1.0 - exp(-radng[i]);
1019         else
1020           probRad = 0.0;
1021 
1022         //...........................................t-z coordinates
1023         if (tauswitch == 0) {
1024 
1025           pc0[1] = P[1][i];
1026           pc0[2] = P[2][i];
1027           pc0[3] = P[3][i];
1028           pc0[0] = P[0][i];
1029 
1030           //                V[0][i]=V[0][i]-dt*RTE/0.1970*sqrt(pow(P[1][i],2)+pow(P[2][i],2)+pow(P[3][i],2))/P[0][i];
1031           //          Ncoll22=dt*RTE/0.197*((pc0[0]*1-pc0[1]*vc0[1]-pc0[2]*vc0[2]-pc0[3]*vc0[3])/E);
1032 
1033           //??????@@@
1034           //////////////////
1035           V[0][i] = V[0][i] - fraction * dt * RTE / 0.1970 *
1036                                   sqrt(pow(P[1][i], 2) + pow(P[2][i], 2) +
1037                                        pow(P[3][i], 2)) /
1038                                   P[0][i];
1039           probCol = fraction * dt_lrf * RTE / 0.1970;
1040         }
1041         ////...........................................tau-eta coordinates
1042         //if(tauswitch==1) { // not ready for JETSCAPE
1043 
1044         //  V[0][i]=V[0][i]-dt*RTE*Xtau/0.1970*sqrt(pow(P[1][i],2)+pow(P[2][i],2)+pow(P[3][i],2))/P[0][i];
1045         //  probCol=dt_lrf*RTE*Xtau/0.1970;  // ?? fraction
1046         //  //        Ncoll22=dt*RTE/0.197*Xtau*((pc0[0]*1-pc0[1]*vc0[1]-pc0[2]*vc0[2]-pc0[3]*vc0[3])/E);
1047         //}
1048 
1049         ////////////////////////////////////////
1050 
1051         if (KINT0 == 0)
1052           probRad = 0.0; // switch off radiation
1053         if (run_alphas==1) {
1054             probCol = probCol * KPfactor * KTfactor * runKT;
1055         } else {
1056             probCol = probCol * KPfactor * KTfactor * preKT;
1057         }
1058         probCol = (1.0 - exp(-probCol)) *
1059                   (1.0 - probRad); // probability of pure elastic scattering
1060         if (KINT0 == 2)
1061           probCol = 0.0;
1062         probTot = probCol + probRad;
1063 
1064         if (ZeroOneDistribution(*GetMt19937Generator()) <
1065             probTot) { // !Yes, collision! Either elastic or inelastic.
1066 
1067           flagScatter = 1;
1068 
1069           //                  Nscatter=KPoisson(Ncoll22);
1070           Nscatter = 1;
1071 
1072           for (int nsca = 1; nsca <= Nscatter; nsca++) {
1073 
1074             np0 = np0 + 1;
1075             pc0[1] = P[1][i];
1076             pc0[2] = P[2][i];
1077             pc0[3] = P[3][i];
1078             pc0[0] = P[0][i];
1079 
1080             flavor(CT, KATTC0, KATT2, KATT3, RTE, PLen, T, T1, T2, E1, E2, iT1,
1081                    iT2, iE1, iE2);
1082 
1083             //...........collision between a jet parton or a recoiled parton and a thermal parton from the medium
1084 
1085             if (CT == 11 || CT == 12) { // for heavy quark scattering
1086               collHQ22(CT, temp0, qhat0, vc0, pc0, pc2, pc3, pc4, qt);
1087             } else { // for light parton scattering
1088               colljet22(CT, temp0, qhat0, vc0, pc0, pc2, pc3, pc4, qt);
1089             }
1090 
1091             icl22 = 1;
1092             n_coll22 += 1;
1093 
1094             KATT1[i] = KATTC0;
1095             KATT1[np0] = KATT2;
1096             KATT10[np0] = KATT3;
1097 
1098             //        if(pc0[0]<pc2[0] && abs(KATTC0)!=4) { // for the purpose of unit test
1099             if (pc0[0] < pc2[0] && abs(KATTC0) != 4 && abs(KATTC0) != 5 &&
1100                 KATTC0 ==
1101                     KATT2) { //disable switch for heavy quark, only allow switch for identical particles
1102               for (int k = 0; k <= 3; k++) {
1103                 p0temp[k] = pc2[k];
1104                 pc2[k] = pc0[k];
1105                 pc0[k] = p0temp[k];
1106               }
1107               KATT1[i] = KATT2;
1108               KATT1[np0] = KATTC0;
1109             }
1110 
1111             for (int j = 0; j <= 3; j++) {
1112 
1113               P[j][i] = pc0[j];
1114               P[j][np0] = pc2[j];
1115               V[j][np0] = V[j][i];
1116 
1117               P0[j][np0] = pc3[j];
1118               V0[j][np0] = V[j][i];
1119 
1120               Vfrozen[j][np0] = Vfrozen[j][i];
1121               Vfrozen0[j][np0] = Vfrozen[j][i];
1122               Tfrozen[np0] = temp0;
1123               Tfrozen0[np0] = temp0;
1124 
1125               if (j != 0) {
1126                 vcfrozen[j][np0] = vc0b[j];
1127                 vcfrozen0[j][np0] = vc0b[j];
1128               }
1129             }
1130 
1131             // pass initial pT and weight information from jet parton to recoil and negative parton
1132             P[5][np0] = P[5][i];
1133             P0[5][np0] = P[5][i];
1134             WT[np0] = WT[i];
1135             WT0[np0] = WT[i];
1136             P[4][np0] = 0.0; // recoiled parton is always massless
1137             P[6][np0] = 0.0; // recoiled parton has 0 virtuality
1138             P0[4][np0] = 0.0;
1139             P0[6][np0] = 0.0;
1140 
1141             CAT[np0] = 2;
1142 
1143           } //for(int nsca=1;nsca<=Nscatter;nsca++)
1144 
1145           //...inelastic scattering begins
1146 
1147           if (qt < epsilon) {
1148             KINT = 0;
1149           } else
1150             KINT = 1;
1151 
1152           if (KINT0 != 0 && KINT != 0) { // Switch on the radiation processes
1153             for (int j = 0; j <= 3; j++) {
1154               p0temp[j] = P[j][i];
1155             }
1156 
1157             px0 = pc4[1];
1158             py0 = pc4[2];
1159             pz0 = pc4[3];
1160 
1161             Elab = pc4
1162                 [0]; // p4-the initial 4-momentum of the jet parton in the lab frame
1163             trans(vc0, pc4);
1164             Ejp = pc4
1165                 [0]; //  p4-the initial 4-momentum of the jet parton in the local rest frame
1166             transback(vc0, pc4);
1167 
1168             if (Ejp > 2 * sqrt(qhat0)) { // !radiation or not?
1169 
1170               Vtemp[1] = xcar;
1171               Vtemp[2] = ycar;
1172               Vtemp[3] = zcar;
1173 
1174               if (KATT1[i] != 21) {
1175                 isp = 1;
1176                 n_sp1 += 1;
1177               } else {
1178                 isp = 2;
1179                 n_sp2 += 1;
1180               }
1181 
1182               if (ZeroOneDistribution(*GetMt19937Generator()) <
1183                   probRad /
1184                       probTot) { // radiation -- either heavy or light parton:
1185 
1186                 for (int j = 0; j <= 3; j++) { // prepare for the collHQ23
1187                   pc01[j] = pc4[j];            // 2->3 process
1188                   pb[j] = pc4[j];
1189                 }
1190 
1191                 collHQ23(KATT1[i], temp0, qhat0, vc0, pc01, pc2, pc3, pc4, qt,
1192                          icl23, Tdiff, Ejp, maxFncHQ, lim_low, lim_int);
1193                 n_coll23 += 1;
1194 
1195                 if (icl23 != 1) {
1196 
1197                   int ctGluon = 1;
1198 
1199                   ng_coll23 += 1;
1200                   nrad = KPoisson(radng[i]);
1201                   int nrad0 = nrad;
1202 
1203                   //                          nrad=1; // only allow single gluon emission right now
1204 
1205                   ng_nrad += nrad;
1206                   np0 = np0 + 1;
1207                   KATT1[np0] = 21;
1208 
1209                   for (int j = 0; j <= 3; j++) {
1210                     P[j][i] = pc01[j];      // heavy quark after colljet23
1211                     P[j][np0 - 1] = pc2[j]; // replace the final state of 22
1212                     V[j][np0 - 1] = V[j][i];
1213                     P0[j][np0 - 1] = pc3[j];
1214                     V0[j][np0 - 1] = V[j][i];
1215                     P[j][np0] = pc4[j]; //radiated gluon from colljet23
1216                     V[j][np0] = V[j][i];
1217                     P0[j][np0] = 0.0;
1218                     V0[j][np0] = 0.0;
1219 
1220                     Vfrozen[j][np0] = Vfrozen[j][i];
1221                     Vfrozen0[j][np0] = 0.0;
1222                     if (j != 0) {
1223                       vcfrozen[j][np0] = vc0b[j];
1224                       vcfrozen0[j][np0] = 0.0;
1225                     }
1226                   }
1227 
1228                   Tfrozen[np0] = temp0;
1229                   Tfrozen0[np0] = 0.0;
1230 
1231                   // pass initial pT and weight information
1232                   P[5][np0] = P[5][i];
1233                   P0[5][np0] = 0.0;
1234                   WT[np0] = WT[i];
1235                   WT0[np0] = 0.0;  // doesn't exist
1236                   P[4][np0] = 0.0; // radiated gluons are massless
1237                   P0[4][np0] = 0.0;
1238                   // assign virtuality for daughter partons
1239                   Qinit = P[6][i];
1240                   P[6][i] = pow(P[0][i] / Elab, 2) * Qinit;
1241                   P[6][np0] = pow(P[0][np0] / Elab, 2) * Qinit;
1242                   P0[6][np0] = 0.0;
1243 
1244                   if (CAT[i] == 2)
1245                     CAT[np0] = 2; // daughters of recoil are recoils
1246 
1247                   // comment out for unit test
1248                   Tint_lrf[i] =
1249                       0.0; //reset radiation infomation for heavy quark
1250 
1251                   flagScatter = 2;
1252 
1253                   eGluon = eGluon + pc4[0];
1254                   nGluon = nGluon + 1.0;
1255 
1256                   V[0][np0] = -log(1.0 - ZeroOneDistribution(*GetMt19937Generator()));
1257 
1258                   // add multiple radiation for heavy quark
1259                   while (nrad > 1) {
1260 
1261                     for (int j = 0; j <= 3; j++)
1262                       pc2[j] = pc4[j];
1263 
1264                     radiationHQ(KATT1[i], qhat0, vc0, pc2, pc01, pc4, pb,
1265                                 iclrad, Tdiff, Ejp, maxFncHQ, temp0, lim_low,
1266                                 lim_int);
1267                     n_radiation += 1;
1268 
1269                     if (iclrad != 1) {
1270                       np0 = np0 + 1;
1271                       ng_radiation += 1;
1272                       KATT1[np0] = 21;
1273 
1274                       ctGluon++;
1275 
1276                       for (int j = 0; j <= 3; j++) {
1277                         P[j][i] = pc01[j]; // heavy quark after one more gluon
1278                         P[j][np0 - 1] = pc2[j]; // update the previous gluon
1279                         P[j][np0] = pc4[j];     // record the new gluon
1280                         V[j][np0] = V[j][i];
1281                         P0[j][np0] = 0.0;
1282                         V0[j][np0] = 0.0;
1283 
1284                         Vfrozen[j][np0] = Vfrozen[j][i];
1285                         Vfrozen0[j][np0] = 0.0;
1286                         if (j != 0) {
1287                           vcfrozen[j][np0] = vc0b[j];
1288                           vcfrozen0[j][np0] = 0.0;
1289                         }
1290                       }
1291 
1292                       Tfrozen[np0] = temp0;
1293                       Tfrozen0[np0] = 0.0;
1294 
1295                       // pass initial pT and weight information
1296                       P[5][np0] = P[5][i];
1297                       P0[5][np0] = 0.0;
1298                       WT[np0] = WT[i];
1299                       WT0[np0] = 0.0;  // doesn't exist
1300                       P[4][np0] = 0.0; // radiated gluons are massless
1301                       P0[4][np0] = 0.0;
1302                       // assign virtuality for daughter partons
1303                       P[6][i] = pow(P[0][i] / Elab, 2) * Qinit;
1304                       P[6][np0 - 1] = pow(P[0][np0 - 1] / Elab, 2) * Qinit;
1305                       P[6][np0] = pow(P[0][np0] / Elab, 2) * Qinit;
1306                       P0[6][np0] = 0.0;
1307                       if (CAT[i] == 2)
1308                         CAT[np0] = 2; // daughters of recoil are recoils
1309 
1310                       eGluon = eGluon + pc4[0];
1311                       nGluon = nGluon + 1.0;
1312 
1313                       V[0][np0] = -log(1.0 - ZeroOneDistribution(*GetMt19937Generator()));
1314 
1315                     } else { //end multiple radiation
1316                       break;
1317                     } //if(icl!=1)radiation
1318 
1319                     nrad = nrad - 1;
1320                   } //multiple radiation for heavy quark ends
1321 
1322                   //                          cout<<"radiate! <Ng>: "<<nrad0<<"  real Ng H: "<< ctGluon << endl;
1323                 } // icl23 == 1, 1st gluon radiation from heavy quark
1324 
1325               } // if(ZeroOneDistribution(*GetMt19937Generator())<probRad/probTot)
1326 
1327             } //if(Ejp>2*sqrt(qhat0))
1328 
1329           } //if(KINT!=0)
1330 
1331           //...........collision end
1332           //...........determine the leading partons and set radng[i]
1333 
1334           //....tiscatter information
1335           if (abs(KATT1[i]) == 4 || abs(KATT1[i]) == 5)
1336             radng[i] = 0.0; // do it below
1337           tiscatter[i] = tcar;
1338           V[0][i] = -log(1.0 - ZeroOneDistribution(*GetMt19937Generator()));
1339 
1340           for (unsigned ip = nnpp + 1; ip <= np0; ++ip) {
1341             tiscatter[ip] = tcar;
1342             V[0][ip] = -log(1.0 - ZeroOneDistribution(*GetMt19937Generator()));
1343             V0[0][ip] = -log(1.0 - ZeroOneDistribution(*GetMt19937Generator()));
1344             tirad[ip] = tcar;
1345             Tint_lrf[ip] = 0.0;
1346             radng[ip] = 0.0;
1347           }
1348 
1349           // SC: only do possible switch for g->ggg... here
1350           if (KATT1[i] == 21 && np0 - nnpp >= 2) {
1351             idlead = nnpp + 2;
1352             p0temp1 = P[0][idlead];
1353             for (unsigned ip = nnpp + 2; ip <= np0 - 1; ++ip) {
1354               if (p0temp1 < P[0][ip + 1]) {
1355                 p0temp1 = P[0][ip + 1];
1356                 idlead = ip + 1;
1357               }
1358             }
1359             if (P[0][idlead] > P[0][i]) {
1360               KATTx = KATT1[i];
1361               KATTy = KATT1[idlead];
1362               KATT1[i] = KATTy;
1363               KATT1[idlead] = KATTx;
1364               for (int j = 0; j <= 3; j++) {
1365                 pcx[j] = P[j][i];
1366                 pcy[j] = P[j][idlead];
1367                 P[j][i] = pcy[j];
1368                 P[j][idlead] = pcx[j];
1369               }
1370             }
1371           }
1372 
1373           //...........sorting block i np~np0(positive) np+1(negative)
1374 
1375           //........................................................................................................
1376 
1377           //CAT!!!
1378           /*              
1379        */
1380           //for(int m=nnpp+1;m<=np0;m++) { // only put recoil parton CAT as 2, radiated gluons do not count
1381           //  if(CAT[i]==2) {
1382           //    CAT[m]=2;
1383           //  }
1384           //}
1385 
1386           if (P[0][nnpp + 1] < Ecut) {
1387             //  CAT[nnpp+1]=1;
1388             //  CAT0[nnpp+1]=1;
1389             for (int j = 0; j <= 3; j++) {
1390               P[j][nnpp + 1] = 0.0;
1391               P0[j][nnpp + 1] = 0.0;
1392             }
1393           }
1394 
1395           if (P[0][i] < Ecut) {
1396             //  CAT[i]=1;
1397             for (int j = 0; j <= 3; j++)
1398               P[j][i] = 0.0;
1399           }
1400 
1401           for (int m = nnpp + 2; m <= np0; m++) {
1402             if (P[0][m] < Ecut) {
1403               //    CAT[m]=1;
1404               for (int j = 0; j <= 3; j++)
1405                 P[j][m] = 0.0;
1406             }
1407           }
1408 
1409         } //...!Yes, collision!
1410 
1411         // even if there is no scattering, one should reset last scatter time now in the new framework
1412         tiscatter[i] = tcar;
1413         radng[i] = 0.0;
1414 
1415       } //....for if(free==0)
1416     }   //..........if CAT[i]=0 end
1417 
1418     nnpp = np0;
1419 
1420   } //......for np end
1421 
1422   //........time step end, np: the number of particles at this point
1423   np = np0;
1424 
1425   if (np >= dimParList) {
1426     cout << "np exceeds the grid size ... terminate " << endl;
1427     exit(EXIT_FAILURE);
1428   }
1429 
1430   if (Kprimary == 1) {
1431     np = nj;
1432   }
1433 }
1434 
1435 ///////////////////////////////////////////////////////////////////////////////////////////////////
1436 //
1437 // Define functions for LBT
1438 //
1439 ///////////////////////////////////////////////////////////////////////////////////////////////////
1440 
1441 void LBT::read_tables() { // intialize various tables for LBT
1442 
1443   //     if(bulkFlag==1) { // read in OSU hydro profiles
1444   //         int dataID_in=1;
1445   //         char dataFN_in[]="../hydroProfile/JetData.dat";
1446   //         int ctlID_in=2;
1447   //         char ctlFN_in[]="../hydroProfile/JetCtl.dat";
1448   //         int bufferSize=1000;
1449   //         int len1=strlen(dataFN_in);
1450   //         int len2=strlen(ctlFN_in);
1451   //         sethydrofilesez_(&dataID_in, dataFN_in, &ctlID_in, ctlFN_in, &bufferSize, len1, len2);
1452   //     } else if(bulkFlag==2) { // read in CCNU hydro profiles
1453   //         char dataFN_in[]="../hydroProfile/bulk.dat";
1454   //         int len1=strlen(dataFN_in);
1455   //         read_ccnu_(dataFN_in,len1);
1456   //     }
1457 
1458   if (fixPosition != 1)
1459     read_xyMC(numInitXY);
1460 
1461   //...read scattering rate
1462   int it, ie;
1463   int n = 450;
1464   ifstream f1("LBT-tables/ratedata");
1465   if (!f1.is_open()) {
1466     cout << "Erro openning date file1!\n";
1467   } else {
1468     for (int i = 1; i <= n; i++) {
1469       f1 >> it >> ie;
1470       f1 >> qhatG[it][ie] >> Rg[it][ie] >> Rg1[it][ie] >> Rg2[it][ie] >>
1471           Rg3[it][ie] >> qhatLQ[it][ie] >> Rq[it][ie] >> Rq3[it][ie] >>
1472           Rq4[it][ie] >> Rq5[it][ie] >> Rq6[it][ie] >> Rq7[it][ie] >>
1473           Rq8[it][ie];
1474     }
1475   }
1476   f1.close();
1477 
1478   // duplicate for heavy quark
1479   ifstream f11("LBT-tables/ratedata-HQ");
1480   if (!f11.is_open()) {
1481     cout << "Erro openning HQ data file!\n";
1482   } else {
1483     for (int i = 1; i <= n; i++) {
1484       f11 >> it >> ie;
1485       f11 >> RHQ[it][ie] >> RHQ11[it][ie] >> RHQ12[it][ie] >> qhatHQ[it][ie];
1486     }
1487   }
1488   f11.close();
1489 
1490   // read radiation table for heavy quark
1491   if (KINT0 != 0) {
1492     ifstream f12("LBT-tables/dNg_over_dt_cD6.dat");
1493     ifstream f13("LBT-tables/dNg_over_dt_qD6.dat");
1494     ifstream f14("LBT-tables/dNg_over_dt_gD6.dat");
1495     if (!f12.is_open() || !f13.is_open() || !f14.is_open()) {
1496       cout << "Erro openning HQ radiation table file!\n";
1497     } else {
1498       for (int k = 1; k <= t_gn; k++) {
1499         char dummyChar[100];
1500         long double dummyD;
1501         f12 >> dummyChar >> dummyChar >> dummyChar >> dummyChar;
1502         f13 >> dummyChar >> dummyChar >> dummyChar >> dummyChar;
1503         f14 >> dummyChar >> dummyChar >> dummyChar >> dummyChar;
1504         for (int i = 1; i <= temp_gn; i++) {
1505           dNg_over_dt_c[k + 1][i][0] = 0.0;
1506           dNg_over_dt_q[k + 1][i][0] = 0.0;
1507           dNg_over_dt_g[k + 1][i][0] = 0.0;
1508           max_dNgfnc_c[k + 1][i][0] = 0.0;
1509           max_dNgfnc_q[k + 1][i][0] = 0.0;
1510           max_dNgfnc_g[k + 1][i][0] = 0.0;
1511           for (int j = 1; j <= HQener_gn; j++) {
1512             f12 >> dNg_over_dt_c[k + 1][i][j] >> max_dNgfnc_c[k + 1][i][j];
1513             f13 >> dNg_over_dt_q[k + 1][i][j] >> max_dNgfnc_q[k + 1][i][j];
1514             f14 >> dNg_over_dt_g[k + 1][i][j] >> max_dNgfnc_g[k + 1][i][j];
1515           }
1516         }
1517       }
1518     }
1519     //     cout << dNg_over_dt_c[t_gn+1][temp_gn][HQener_gn] << "    " << max_dNgfnc_c[t_gn+1][temp_gn][HQener_gn] << endl;
1520     //     cout << dNg_over_dt_q[t_gn+1][temp_gn][HQener_gn] << "    " << max_dNgfnc_q[t_gn+1][temp_gn][HQener_gn] << endl;
1521     //     cout << dNg_over_dt_g[t_gn+1][temp_gn][HQener_gn] << "    " << max_dNgfnc_g[t_gn+1][temp_gn][HQener_gn] << endl;
1522     f12.close();
1523     f13.close();
1524     f14.close();
1525 
1526     for (int i = 1; i <= temp_gn; i++) {
1527       for (int j = 1; j <= HQener_gn; j++) {
1528         dNg_over_dt_c[1][i][j] = 0.0;
1529         dNg_over_dt_q[1][i][j] = 0.0;
1530         dNg_over_dt_g[1][i][j] = 0.0;
1531         max_dNgfnc_c[1][i][j] = 0.0;
1532         max_dNgfnc_q[1][i][j] = 0.0;
1533         max_dNgfnc_g[1][i][j] = 0.0;
1534       }
1535     }
1536   }
1537 
1538   // preparation for HQ 2->2
1539   ifstream fileB("LBT-tables/distB.dat");
1540   if (!fileB.is_open()) {
1541     cout << "Erro openning data file distB.dat!" << endl;
1542   } else {
1543     for (int i = 0; i < N_T; i++) {
1544       for (int j = 0; j < N_p1; j++) {
1545         double dummy_T, dummy_p1;
1546         fileB >> dummy_T >> dummy_p1;
1547         if (fabs(min_T + (0.5 + i) * bin_T - dummy_T) > 1.0e-5 ||
1548             fabs(min_p1 + (0.5 + j) * bin_p1 - dummy_p1) > 1.0e-5) {
1549           cout << "Erro in reading data file distB.dat!" << endl;
1550           exit(EXIT_FAILURE);
1551         }
1552         fileB >> distFncBM[i][j];
1553         for (int k = 0; k < N_e2; k++)
1554           fileB >> distFncB[i][j][k];
1555         for (int k = 0; k < N_e2; k++)
1556           fileB >> distMaxB[i][j][k];
1557       }
1558     }
1559   }
1560   fileB.close();
1561 
1562   ifstream fileF("LBT-tables/distF.dat");
1563   if (!fileF.is_open()) {
1564     cout << "Erro openning data file distF.dat!" << endl;
1565   } else {
1566     for (int i = 0; i < N_T; i++) {
1567       for (int j = 0; j < N_p1; j++) {
1568         double dummy_T, dummy_p1;
1569         fileF >> dummy_T >> dummy_p1;
1570         if (fabs(min_T + (0.5 + i) * bin_T - dummy_T) > 1.0e-5 ||
1571             fabs(min_p1 + (0.5 + j) * bin_p1 - dummy_p1) > 1.0e-5) {
1572           cout << "Erro in reading data file distF.dat!" << endl;
1573           exit(EXIT_FAILURE);
1574         }
1575         fileF >> distFncFM[i][j];
1576         for (int k = 0; k < N_e2; k++)
1577           fileF >> distFncF[i][j][k];
1578         for (int k = 0; k < N_e2; k++)
1579           fileF >> distMaxF[i][j][k];
1580       }
1581     }
1582   }
1583   fileF.close();
1584 
1585   cout << "Initialization completed for LBT." << endl;
1586 }
1587 
1588 //..............................................................subroutine
1589 //..............................................................
1590 
1591 float LBT::ran0(long *idum)
1592 
1593 {
1594   const int IM1 = 2147483563;
1595   const int IM2 = 2147483399;
1596   const double AM = (1.0 / IM1);
1597   const int IMM1 = (IM1 - 1);
1598   const int IA1 = 40014;
1599   const int IA2 = 40692;
1600   const int IQ1 = 53668;
1601   const int IQ2 = 52774;
1602   const int IR1 = 12211;
1603   const int IR2 = 3791;
1604   const int NTAB = 32;
1605   const int NDIV = (1 + IMM1 / NTAB);
1606   const double EPS = 1.2e-7;
1607   const double RNMX = (1.0 - EPS);
1608 
1609   int j;
1610   long k;
1611   static long idum2 = 123456789;
1612   static long iy = 0;
1613   static long iv[NTAB];
1614   float temp;
1615 
1616   if (*idum <= 0) {
1617     if (-(*idum) < 1)
1618       *idum = 1;
1619     else
1620       *idum = -(*idum);
1621     for (j = NTAB + 7; j >= 0; j--) {
1622       k = (*idum) / IQ1;
1623       *idum = IA1 * (*idum - k * IQ1) - k * IR1;
1624       if (*idum < 0)
1625         *idum += IM1;
1626       if (j < NTAB)
1627         iv[j] = *idum;
1628     }
1629     iy = iv[0];
1630   }
1631   k = (*idum) / IQ1;
1632   *idum = IA1 * (*idum - k * IQ1) - k * IR1;
1633   if (*idum < 0)
1634     *idum += IM1;
1635   k = idum2 / IQ2;
1636   idum2 = IA2 * (idum2 - k * IQ2) - k * IR2;
1637   if (idum2 < 0)
1638     idum2 += IM2;
1639   j = iy / NDIV;
1640   iy = iv[j] - idum2;
1641   iv[j] = *idum;
1642   if (iy < 1)
1643     iy += IMM1;
1644   if ((temp = AM * iy) > RNMX)
1645     return RNMX;
1646   else
1647     return temp;
1648 }
1649 
1650 //.........................................................................
1651 double LBT::alphas0(int &Kalphas, double temp0) {
1652 
1653   double X;
1654   if (Kalphas == 1) {
1655     X = fixAlphas;
1656   } else {
1657     cout << "un-recoganized Kalphas" << endl;
1658     exit(EXIT_FAILURE);
1659   }
1660   return X;
1661 }
1662 //.........................................................................
1663 double LBT::DebyeMass2(int &Kqhat0, double alphas, double temp0) {
1664 
1665   switch (Kqhat0) {
1666   case 1:
1667     return 4.0 * pi * alphas * pow(temp0, 2);
1668     break;
1669   case 2:
1670     return (3.0 / 2.0) * 4.0 * pi * alphas * pow(temp0, 2);
1671     break;
1672   case 3:
1673     return 1.0;
1674   default:
1675     JSWARN << "Kqhat0 = " << Kqhat0;
1676     throw std::runtime_error("Unexpected value for Kqhat0");
1677     return -1;
1678   }
1679 
1680   // double Y;
1681   // if(Kqhat0==1)
1682   //   {
1683   //     Y=4.0*pi*alphas*pow(temp0,2);
1684   //   }
1685   // if(Kqhat0==2)
1686   //   {
1687   //     Y=(3.0/2.0)*4.0*pi*alphas*pow(temp0,2);
1688   //   }
1689   // if(Kqhat0==3)
1690   //   {
1691   //     Y=1.0;
1692   //   }
1693   // return Y;
1694 }
1695 //.........................................................................
1696 
1697 //.........................................................................
1698 void LBT::titau(double ti, double vf[4], double vp[4], double p0[4], double &Vx,
1699                 double &Vy, double &Veta, double &Xtau) {
1700 
1701   //..............................................................test part
1702   //          cout<<"ti"<<" "<<ti<<" "<<"vf"<<" "<<vf[1]<<" "<<vf[2]<<" "<<vf[3]<<endl;
1703   //          cout<<"vp[4]"<<" "<<vp[1]<<" "<<vp[2]<<" "<<vp[3]<<" "<<vp[0]<<endl;
1704   //          cout<<"p0[4]"<<" "<<p0[1]<<" "<<p0[2]<<" "<<p0[3]<<" "<<p0[0]<<endl;
1705   //..............................................................test part
1706 
1707   //....notice the form of vf
1708   double gamma = 1.0 / sqrt(1 - (vf[1] * vf[1] + vf[2] * vf[2]));
1709   double mt = sqrt(p0[1] * p0[1] + p0[2] * p0[2]);
1710   double Yp = 1.0 / 2.0 * log((p0[0] + p0[3]) / (p0[0] - p0[3]));
1711   double etas = vp[3];
1712   double etaf = atanh(vf[3]) + etas;
1713   double pper = sqrt(p0[1] * p0[1] + p0[2] * p0[2]);
1714   double vper = sqrt(vf[1] * vf[1] + vf[2] * vf[2]);
1715   double pvper = p0[1] * vf[1] + p0[2] * vf[2];
1716 
1717   Vx = p0[1] / pper / cosh(Yp - etas);
1718   Vy = p0[2] / pper / cosh(Yp - etas);
1719   Veta = (1.0 / ti) * tanh(Yp - etas);
1720 
1721   Xtau = (gamma * mt * cosh(Yp - etaf) - pvper * vper) / (mt * cosh(Yp - etas));
1722 
1723   //..............................................................test part
1724   //          cout<<"gamma"<<" "<<gamma<<" "<<"mt"<<" "<<mt<<" "<<"Yp"<<" "<<Yp<<endl;
1725   //          cout<<"etas"<<" "<<etas<<" "<<"etaf"<<" "<<etaf<<" "<<"pper"<<" "<<pper<<endl;
1726   //          cout<<"vper"<<" "<<vper<<" "<<"pvper"<<" "<<pvper<<endl;
1727   //          cout<<"Vx"<<" "<<Vx<<" "<<"Vy"<<" "<<Vy<<" "<<"Veta"<<" "<<Veta<<endl;
1728   //          cout<<"Xtau"<<" "<<Xtau<<endl;
1729   //..............................................................test part
1730 }
1731 //.........................................................................
1732 
1733 //.........................................................................
1734 
1735 void LBT::lam(int KATT0, double &RTE, double E, double T, double &T1,
1736               double &T2, double &E1, double &E2, int &iT1, int &iT2, int &iE1,
1737               int &iE2) {
1738 
1739   double dtemp = 0.02;
1740   iT1 = (int)((T - 0.1) / dtemp);
1741   iT2 = iT1 + 1;
1742   iE1 = (int)(log(E) + 2);
1743   if (iE1 < 1)
1744     iE1 = 1;
1745   iE2 = iE1 + 1;
1746 
1747   T1 = 0.12 + (iT1 - 1) * 0.02;
1748   T2 = T1 + dtemp;
1749   E1 = exp(iE1 - 2.0);
1750   E2 = exp(iE2 - 2.0);
1751 
1752   if (KATT0 == 21) {
1753     double RTE1 =
1754         (Rg[iT2][iE1] - Rg[iT1][iE1]) * (T - T1) / (T2 - T1) + Rg[iT1][iE1];
1755     double RTE2 =
1756         (Rg[iT2][iE2] - Rg[iT1][iE2]) * (T - T1) / (T2 - T1) + Rg[iT1][iE2];
1757     RTE = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
1758   } else if (KATT0 == 4 || KATT0 == -4 || KATT0 == 5 ||
1759              KATT0 == -5) { // add heavy quark channel
1760     double RTE1 =
1761         (RHQ[iT2][iE1] - RHQ[iT1][iE1]) * (T - T1) / (T2 - T1) + RHQ[iT1][iE1];
1762     double RTE2 =
1763         (RHQ[iT2][iE2] - RHQ[iT1][iE2]) * (T - T1) / (T2 - T1) + RHQ[iT1][iE2];
1764     RTE = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
1765   } else {
1766     double RTE1 =
1767         (Rq[iT2][iE1] - Rq[iT1][iE1]) * (T - T1) / (T2 - T1) + Rq[iT1][iE1];
1768     double RTE2 =
1769         (Rq[iT2][iE2] - Rq[iT1][iE2]) * (T - T1) / (T2 - T1) + Rq[iT1][iE2];
1770     RTE = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
1771     //          cout<<"RTE2,RTE1,E,E1,E2,RTE: "<<RTE2<<"  "<<RTE1<<"  "<<E<<"  "<<E1<<"  "<<E2<<"  "<<RTE<<endl;
1772   }
1773 }
1774 
1775 //.........................................................................
1776 
1777 void LBT::flavor(int &CT, int &KATT0, int &KATT2, int &KATT3, double RTE,
1778                  double E, double T, double &T1, double &T2, double &E1,
1779                  double &E2, int &iT1, int &iT2, int &iE1, int &iE2) {
1780 
1781   double RTEg;
1782   double RTEg1;
1783   double RTEg2;
1784   double RTEg3;
1785   double RTEq;
1786   double RTEq3;
1787   double RTEq4;
1788   double RTEq5;
1789   double RTEq6;
1790   double RTEq7;
1791   double RTEq8;
1792   double RTEHQ;
1793   double RTEHQ11;
1794   double RTEHQ12;
1795   double qhatTP;
1796 
1797   int vb[7] = {0};
1798   int b = 0;
1799   int KATT00 = KATT0;
1800 
1801   vb[1] = 1;
1802   vb[2] = 2;
1803   vb[3] = 3;
1804   vb[4] = -1;
1805   vb[5] = -2;
1806   vb[6] = -3;
1807 
1808   //.....
1809   linear(KATT0, E, T, T1, T2, E1, E2, iT1, iT2, iE1, iE2, RTEg, RTEg1, RTEg2,
1810          RTEg3, RTEq, RTEq3, RTEq4, RTEq5, RTEq6, RTEq7, RTEq8, RTEHQ, RTEHQ11,
1811          RTEHQ12, qhatTP);
1812 
1813   if (KATT00 == 21) { //.....for gluon
1814     double R0 = RTE;
1815     double R1 = RTEg1;
1816     double R2 = RTEg2;
1817     double R3 = RTEg3;
1818 
1819     double a = ZeroOneDistribution(*GetMt19937Generator());
1820 
1821     if (a <= R1 / R0) {
1822       CT = 1;
1823       KATT3 = 21;
1824       KATT2 = 21;
1825       KATT0 = 21;
1826     }
1827 
1828     if (a > R1 / R0 && a <= (R1 + R2) / R0) {
1829       CT = 2;
1830       b = floor(ZeroOneDistribution(*GetMt19937Generator()) * 6 + 1);
1831       if (b == 7) {
1832         b = 6;
1833       }
1834       KATT3 = 21;
1835       KATT2 = vb[b];
1836       KATT0 = -KATT2;
1837     }
1838 
1839     if (a > (R1 + R2) / R0 && a <= 1.0) {
1840       CT = 3;
1841       b = floor(ZeroOneDistribution(*GetMt19937Generator()) * 6 + 1);
1842       if (b == 7) {
1843         b = 6;
1844       }
1845       KATT3 = vb[b];
1846       KATT2 = KATT3;
1847       KATT0 = 21;
1848     }
1849   } else if (KATT00 == 4 || KATT00 == -4 || KATT00 == 5 ||
1850              KATT00 == -5) { // for heavy quark
1851     double R0 = RTE;
1852     double R1 = RTEHQ11;
1853     double R2 = RTEHQ12;
1854 
1855     double a = ZeroOneDistribution(*GetMt19937Generator());
1856 
1857     //          qhat_over_T3=qhatTP;  // what is read in is qhat/T^3 of quark
1858     //          D2piT=8.0*pi/qhat_over_T3;
1859 
1860     if (a <= R1 / R0) { //Qq->Qq
1861       CT = 11;
1862       b = floor(ZeroOneDistribution(*GetMt19937Generator()) * 6 + 1);
1863       if (b == 7) {
1864         b = 6;
1865       }
1866       KATT3 = vb[b];
1867       KATT2 = KATT3;
1868     } else { //Qg->Qg
1869       CT = 12;
1870       KATT3 = 21;
1871       KATT2 = KATT3;
1872     }
1873 
1874   } else { //.....for quark and antiquark (light)
1875     double R00 = RTE;
1876     double R3 = RTEq3;
1877     double R4 = RTEq4;
1878     double R5 = RTEq5;
1879     double R6 = RTEq6;
1880     double R7 = RTEq7;
1881     double R8 = RTEq8;
1882 
1883     double a = ZeroOneDistribution(*GetMt19937Generator());
1884     if (a <= R3 / R00) {
1885       CT = 13;
1886       KATT3 = 21;
1887       KATT2 = 21;
1888       //          KATT0=KATT0;
1889     }
1890 
1891     if (a > R3 / R00 && a <= (R3 + R4) / R00) {
1892       CT = 4;
1893     f1:
1894       b = floor(ZeroOneDistribution(*GetMt19937Generator()) * 6 + 1);
1895       if (b == 7) {
1896         b = 6;
1897       }
1898       KATT3 = vb[b];
1899       if (KATT3 == KATT0) {
1900         goto f1;
1901       }
1902       KATT2 = KATT3;
1903       //          KATT0=KATT0
1904     }
1905 
1906     if (a > (R3 + R4) / R00 && a <= (R3 + R4 + R5) / R00) {
1907       CT = 5;
1908       KATT3 = KATT0;
1909       KATT2 = KATT0;
1910     }
1911     //.....the only difference between quark and antiquark
1912 
1913     if (a > (R3 + R4 + R5) / R00 && a <= (R3 + R4 + R5 + R6) / R00) {
1914       CT = 6;
1915       KATT3 = -KATT0;
1916     f2:
1917       b = floor(ZeroOneDistribution(*GetMt19937Generator()) * 3 + 1);
1918       if (b == 4) {
1919         b = 3;
1920       }
1921       KATT2 = -KATT0 / abs(KATT0) * vb[b];
1922       if (abs(KATT2) == abs(KATT3)) {
1923         goto f2;
1924       }
1925       KATT0 = -KATT2;
1926     }
1927 
1928     if (a > (R3 + R4 + R5 + R6) / R00 && a <= (R3 + R4 + R5 + R6 + R7) / R00) {
1929       CT = 7;
1930       KATT3 = -KATT0;
1931       KATT2 = KATT3;
1932       //          KATT0=KATT0
1933     }
1934 
1935     //  if(a>(R00-R8)/R00 && a<=1.0)
1936     if (a > (R3 + R4 + R5 + R6 + R7) / R00 && a <= 1.0) {
1937       CT = 8;
1938       KATT3 = -KATT0;
1939       KATT2 = 21;
1940       KATT0 = 21;
1941     }
1942   }
1943 }
1944 
1945 //.........................................................................
1946 
1947 void LBT::linear(int KATT, double E, double T, double &T1, double &T2,
1948                  double &E1, double &E2, int &iT1, int &iT2, int &iE1, int &iE2,
1949                  double &RTEg, double &RTEg1, double &RTEg2, double &RTEg3,
1950                  double &RTEq, double &RTEq3, double &RTEq4, double &RTEq5,
1951                  double &RTEq6, double &RTEq7, double &RTEq8, double &RTEHQ,
1952                  double &RTEHQ11, double &RTEHQ12, double &qhatTP) {
1953   if (KATT == 21) {
1954     double RTE1 =
1955         (Rg1[iT2][iE1] - Rg1[iT1][iE1]) * (T - T1) / (T2 - T1) + Rg1[iT1][iE1];
1956     double RTE2 =
1957         (Rg1[iT2][iE2] - Rg1[iT1][iE2]) * (T - T1) / (T2 - T1) + Rg1[iT1][iE2];
1958     RTEg1 = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
1959 
1960     RTE1 =
1961         (Rg2[iT2][iE1] - Rg2[iT1][iE1]) * (T - T1) / (T2 - T1) + Rg2[iT1][iE1];
1962     RTE2 =
1963         (Rg2[iT2][iE2] - Rg2[iT1][iE2]) * (T - T1) / (T2 - T1) + Rg2[iT1][iE2];
1964     RTEg2 = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
1965 
1966     RTE1 =
1967         (Rg3[iT2][iE1] - Rg3[iT1][iE1]) * (T - T1) / (T2 - T1) + Rg3[iT1][iE1];
1968     RTE2 =
1969         (Rg3[iT2][iE2] - Rg3[iT1][iE2]) * (T - T1) / (T2 - T1) + Rg3[iT1][iE2];
1970     RTEg3 = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
1971 
1972     //    RTE1=(qhatG[iT2][iE1]-qhatG[iT1][iE1])*(T-T1)/(T2-T1)+qhatG[iT1][iE1];
1973     //    RTE2=(qhatG[iT2][iE2]-qhatG[iT1][iE2])*(T-T1)/(T2-T1)+qhatG[iT1][iE2];
1974     //    qhatTP=(RTE2-RTE1)*(E-E1)/(E2-E1)+RTE1;
1975 
1976   } else if (KATT == 4 || KATT == -4 || KATT == 5 ||
1977              KATT == -5) { // add heavy quark channel
1978     double RTE1 = (RHQ11[iT2][iE1] - RHQ11[iT1][iE1]) * (T - T1) / (T2 - T1) +
1979                   RHQ11[iT1][iE1];
1980     double RTE2 = (RHQ11[iT2][iE2] - RHQ11[iT1][iE2]) * (T - T1) / (T2 - T1) +
1981                   RHQ11[iT1][iE2];
1982     RTEHQ11 = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
1983 
1984     RTE1 = (RHQ12[iT2][iE1] - RHQ12[iT1][iE1]) * (T - T1) / (T2 - T1) +
1985            RHQ12[iT1][iE1];
1986     RTE2 = (RHQ12[iT2][iE2] - RHQ12[iT1][iE2]) * (T - T1) / (T2 - T1) +
1987            RHQ12[iT1][iE2];
1988     RTEHQ12 = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
1989 
1990     //    RTE1=(qhatHQ[iT2][iE1]-qhatHQ[iT1][iE1])*(T-T1)/(T2-T1)+qhatHQ[iT1][iE1];
1991     //    RTE2=(qhatHQ[iT2][iE2]-qhatHQ[iT1][iE2])*(T-T1)/(T2-T1)+qhatHQ[iT1][iE2];
1992     //    qhatTP=(RTE2-RTE1)*(E-E1)/(E2-E1)+RTE1;
1993 
1994   } else {
1995     double RTE1 =
1996         (Rq3[iT2][iE1] - Rq3[iT1][iE1]) * (T - T1) / (T2 - T1) + Rq3[iT1][iE1];
1997     double RTE2 =
1998         (Rq3[iT2][iE2] - Rq3[iT1][iE2]) * (T - T1) / (T2 - T1) + Rq3[iT1][iE2];
1999     RTEq3 = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
2000 
2001     RTE1 =
2002         (Rq4[iT2][iE1] - Rq4[iT1][iE1]) * (T - T1) / (T2 - T1) + Rq4[iT1][iE1];
2003     RTE2 =
2004         (Rq4[iT2][iE2] - Rq4[iT1][iE2]) * (T - T1) / (T2 - T1) + Rq4[iT1][iE2];
2005     RTEq4 = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
2006 
2007     RTE1 =
2008         (Rq5[iT2][iE1] - Rq5[iT1][iE1]) * (T - T1) / (T2 - T1) + Rq5[iT1][iE1];
2009     RTE2 =
2010         (Rq5[iT2][iE2] - Rq5[iT1][iE2]) * (T - T1) / (T2 - T1) + Rq5[iT1][iE2];
2011     RTEq5 = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
2012 
2013     RTE1 =
2014         (Rq6[iT2][iE1] - Rq6[iT1][iE1]) * (T - T1) / (T2 - T1) + Rq6[iT1][iE1];
2015     RTE2 =
2016         (Rq6[iT2][iE2] - Rq6[iT1][iE2]) * (T - T1) / (T2 - T1) + Rq6[iT1][iE2];
2017     RTEq6 = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
2018 
2019     RTE1 =
2020         (Rq7[iT2][iE1] - Rq7[iT1][iE1]) * (T - T1) / (T2 - T1) + Rq7[iT1][iE1];
2021     RTE2 =
2022         (Rq7[iT2][iE2] - Rq7[iT1][iE2]) * (T - T1) / (T2 - T1) + Rq7[iT1][iE2];
2023     RTEq7 = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
2024 
2025     RTE1 =
2026         (Rq8[iT2][iE1] - Rq8[iT1][iE1]) * (T - T1) / (T2 - T1) + Rq8[iT1][iE1];
2027     RTE2 =
2028         (Rq8[iT2][iE2] - Rq8[iT1][iE2]) * (T - T1) / (T2 - T1) + Rq8[iT1][iE2];
2029     RTEq8 = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
2030 
2031     //    RTE1=(qhatLQ[iT2][iE1]-qhatLQ[iT1][iE1])*(T-T1)/(T2-T1)+qhatLQ[iT1][iE1];
2032     //    RTE2=(qhatLQ[iT2][iE2]-qhatLQ[iT1][iE2])*(T-T1)/(T2-T1)+qhatLQ[iT1][iE2];
2033     //    qhatTP=(RTE2-RTE1)*(E-E1)/(E2-E1)+RTE1;
2034   }
2035 }
2036 
2037 //.........................................................................
2038 
2039 void LBT::twflavor(int &CT, int &KATT0, int &KATT2, double E, double T) {
2040 
2041   double RTEg;
2042   double RTEg1;
2043   double RTEg2;
2044   double RTEg3;
2045   double RTEq;
2046   double RTEq3;
2047   double RTEq4;
2048   double RTEq5;
2049   double RTEq6;
2050   double RTEq7;
2051   double RTEq8;
2052 
2053   int vb[7] = {0};
2054   int b = 0;
2055   int KATT00 = KATT0;
2056   int KATT20 = KATT2;
2057 
2058   vb[1] = 1;
2059   vb[2] = 2;
2060   vb[3] = 3;
2061   vb[4] = -1;
2062   vb[5] = -2;
2063   vb[6] = -3;
2064 
2065   twlinear(KATT0, E, T, RTEg1, RTEg2, RTEq6, RTEq7, RTEq8);
2066 
2067   //.....for gluon
2068   if (KATT00 == 21) {
2069     //  R0  =RTE
2070     double R1 = RTEg1;
2071     double R2 = RTEg2;
2072     //  R3  =RTEg3
2073 
2074     if (KATT20 == 21) {
2075       double a = ZeroOneDistribution(*GetMt19937Generator());
2076       if (a <= R1 / (R1 + R2)) {
2077         CT = 1;
2078         //          KATT3=KATT2
2079         //          KATT2=21
2080         //          KATT0=21
2081       }
2082 
2083       if (a > R1 / (R1 + R2)) {
2084         CT = 2;
2085         //          KATT3=KATT2
2086         b = floor(ZeroOneDistribution(*GetMt19937Generator()) * 6 + 1);
2087         if (b == 7) {
2088           b = 6;
2089         }
2090         KATT2 = vb[b];
2091         KATT0 = -KATT2;
2092       }
2093     }
2094 
2095     if (KATT20 != 21) {
2096       CT = 3;
2097       //          KATT3=KATT2
2098       //          KATT2=KATT2
2099       //          KATT0=21
2100     }
2101   }
2102 
2103   //.....for quark and antiquark
2104   if (KATT00 != 21) {
2105 
2106     //  R00 =RTE
2107     //  R3  =RTEq3
2108     //  R4  =RTEq4
2109     //  R5  =RTEq5
2110     double R6 = RTEq6;
2111     double R7 = RTEq7;
2112     double R8 = RTEq8;
2113     double R00 = R6 + R7 + R8;
2114 
2115     if (KATT20 == 21) {
2116       CT = 13;
2117       //          KATT3=KATT2
2118       //          KATT2=21
2119       //          KATT0=KATT0
2120     }
2121 
2122     if (KATT20 != 21) {
2123 
2124       if (abs(KATT20) != abs(KATT00)) {
2125         CT = 4;
2126         //        KATT3=KATT2
2127         //        KATT2=KATT3
2128         //        KATT0=KATT0
2129       }
2130 
2131       if (KATT20 == KATT00) {
2132         CT = 5;
2133         //        KATT3=KATT2
2134         //        KATT0=KATT0
2135       }
2136 
2137       if (KATT20 == -KATT00) {
2138         double a = ZeroOneDistribution(*GetMt19937Generator());
2139         if (a <= (R6) / R00) {
2140           CT = 6;
2141           //             KATT3=KATT2
2142         tf2:
2143           b = floor(ZeroOneDistribution(*GetMt19937Generator()) * 3 + 1);
2144           if (b == 4) {
2145             b = 3;
2146           }
2147           KATT2 = -KATT0 / abs(KATT0) * vb[b];
2148           if (abs(KATT2) == abs(KATT0)) {
2149             goto tf2;
2150           }
2151           KATT0 = -KATT2;
2152         }
2153 
2154         if (a > (R6) / R00 && a <= (R6 + R7) / R00) {
2155           CT = 7;
2156           //             KATT3=KATT2
2157           //             KATT2=KATT3
2158           //             KATT0=KATT0
2159         }
2160 
2161         if (a > (R6 + R7) / R00 && a <= 1.0) {
2162           CT = 8;
2163           //             KATT3=KATT2
2164           KATT2 = 21;
2165           KATT0 = 21;
2166         }
2167       }
2168     }
2169   }
2170 }
2171 
2172 //.........................................................................
2173 
2174 void LBT::twlinear(int KATT, double E, double T, double &RTEg1, double &RTEg2,
2175                    double &RTEq6, double &RTEq7, double &RTEq8) {
2176 
2177   //.....
2178   double dtemp = 0.02;
2179   int iT1 = floor((T - 0.1) / dtemp);
2180   int iT2 = iT1 + 1;
2181   int iE1 = floor(log(E) + 2);
2182   int iE2 = iE1 + 1;
2183   //
2184   double T1 = 0.12 + (iT1 - 1) * 0.02;
2185   double T2 = T1 + dtemp;
2186   double E1 = exp(iE1 - 2.0);
2187   double E2 = exp(iE2 - 2.0);
2188   //
2189   if (KATT == 21) {
2190     double RTE1 =
2191         (Rg1[iT2][iE1] - Rg1[iT1][iE1]) * (T - T1) / (T2 - T1) + Rg1[iT1][iE1];
2192     double RTE2 =
2193         (Rg1[iT2][iE2] - Rg1[iT1][iE2]) * (T - T1) / (T2 - T1) + Rg1[iT1][iE2];
2194     RTEg1 = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
2195 
2196     RTE1 =
2197         (Rg2[iT2][iE1] - Rg2[iT1][iE1]) * (T - T1) / (T2 - T1) + Rg2[iT1][iE1];
2198     RTE2 =
2199         (Rg2[iT2][iE2] - Rg2[iT1][iE2]) * (T - T1) / (T2 - T1) + Rg2[iT1][iE2];
2200     RTEg2 = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
2201 
2202     //     RTE1=(Rg3[iT2][iE1]-Rg3[iT1][iE1])*(T-T1)/(T2-T1)+Rg3[iT1][iE1];
2203     //     RTE2=(Rg3[iT2][iE2]-Rg3[iT1][iE2])*(T-T1)/(T2-T1)+Rg3[iT1][iE2];
2204     //     RTEg3=(RTE2-RTE1)*(E-E1)/(E2-E1)+RTE1;
2205   }
2206 
2207   if (KATT != 21) {
2208     //     RTE1=(Rq3[iT2][iE1]-Rq3[iT1][iE1])*(T-T1)/(T2-T1)+Rq3[iT1][iE1];
2209     //     RTE2=(Rq3[iT2][iE2]-Rq3[iT1][iE2])*(T-T1)/(T2-T1)+Rq3[iT1][iE2];
2210     //     RTEq3=(RTE2-RTE1)*(E-E1)/(E2-E1)+RTE1;
2211 
2212     //     RTE1=(Rq4[iT2][iE1]-Rq4[iT1][iE1])*(T-T1)/(T2-T1)+Rq4[iT1][iE1];
2213     //     RTE2=(Rq4[iT2][iE2]-Rq4[iT1][iE2])*(T-T1)/(T2-T1)+Rq4[iT1][iE2];
2214     //     RTEq4=(RTE2-RTE1)*(E-E1)/(E2-E1)+RTE1
2215 
2216     //     RTE1=(Rq5[iT2][iE1]-Rq5[iT1][iE1])*(T-T1)/(T2-T1)+Rq5[iT1][iE1];
2217     //     RTE2=(Rq5[iT2][iE2]-Rq5[iT1][iE2])*(T-T1)/(T2-T1)+Rq5[iT1][iE2];
2218     //     RTEq5=(RTE2-RTE1)*(E-E1)/(E2-E1)+RTE1;
2219 
2220     double RTE1 =
2221         (Rq6[iT2][iE1] - Rq6[iT1][iE1]) * (T - T1) / (T2 - T1) + Rq6[iT1][iE1];
2222     double RTE2 =
2223         (Rq6[iT2][iE2] - Rq6[iT1][iE2]) * (T - T1) / (T2 - T1) + Rq6[iT1][iE2];
2224     RTEq6 = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
2225 
2226     RTE1 =
2227         (Rq7[iT2][iE1] - Rq7[iT1][iE1]) * (T - T1) / (T2 - T1) + Rq7[iT1][iE1];
2228     RTE2 =
2229         (Rq7[iT2][iE2] - Rq7[iT1][iE2]) * (T - T1) / (T2 - T1) + Rq7[iT1][iE2];
2230     RTEq7 = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
2231 
2232     RTE1 =
2233         (Rq8[iT2][iE1] - Rq8[iT1][iE1]) * (T - T1) / (T2 - T1) + Rq8[iT1][iE1];
2234     RTE2 =
2235         (Rq8[iT2][iE2] - Rq8[iT1][iE2]) * (T - T1) / (T2 - T1) + Rq8[iT1][iE2];
2236     RTEq8 = (RTE2 - RTE1) * (E - E1) / (E2 - E1) + RTE1;
2237   }
2238 }
2239 
2240 //.........................................................................
2241 
2242 void LBT::trans(double v[4], double p[4]) {
2243   double vv = sqrt(v[1] * v[1] + v[2] * v[2] + v[3] * v[3]);
2244   double ga = 1.0 / sqrt(1.0 - vv * vv);
2245   double ppar = p[1] * v[1] + p[2] * v[2] + p[3] * v[3];
2246   double gavv = (ppar * ga / (1.0 + ga) - p[0]) * ga;
2247   p[0] = ga * (p[0] - ppar);
2248   p[1] = p[1] + v[1] * gavv;
2249   p[2] = p[2] + v[2] * gavv;
2250   p[3] = p[3] + v[3] * gavv;
2251 }
2252 
2253 void LBT::transback(double v[4], double p[4]) {
2254   double vv = sqrt(v[1] * v[1] + v[2] * v[2] + v[3] * v[3]);
2255   double ga = 1.0 / sqrt(1.0 - vv * vv);
2256   double ppar = p[1] * v[1] + p[2] * v[2] + p[3] * v[3];
2257   double gavv = (-ppar * ga / (1.0 + ga) - p[0]) * ga;
2258   p[0] = ga * (p[0] + ppar);
2259   p[1] = p[1] - v[1] * gavv;
2260   p[2] = p[2] - v[2] * gavv;
2261   p[3] = p[3] - v[3] * gavv;
2262 }
2263 
2264 //.........................................................................
2265 void LBT::colljet22(int CT, double temp, double qhat0ud, double v0[4],
2266                     double p0[4], double p2[4], double p3[4], double p4[4],
2267                     double &qt) {
2268   //
2269   //    p0 initial jet momentum, output to final momentum
2270   //    p2 final thermal momentum,p3 initial termal energy
2271   //
2272   //    amss=sqrt(abs(p0(4)**2-p0(1)**2-p0(2)**2-p0(3)**2))
2273   //
2274   //************************************************************
2275   p4[1] = p0[1];
2276   p4[2] = p0[2];
2277   p4[3] = p0[3];
2278   p4[0] = p0[0];
2279   //************************************************************
2280 
2281   //    transform to local comoving frame of the fluid
2282   //  cout << endl;
2283   //  cout << "flow  "<< v0[1] << " " << v0[2] << " " << v0[3] << " "<<" Elab " << p0[0] << endl;
2284 
2285   trans(v0, p0);
2286   //  cout << p0[0] << " " << sqrt(qhat0ud) << endl;
2287 
2288   //  cout << sqrt(pow(p0[1],2)+pow(p0[2],2)+pow(p0[3],2)) << " " << p0[1] << " " << p0[2] << " " << p0[3] << endl;
2289 
2290   //************************************************************
2291   trans(v0, p4);
2292   //************************************************************
2293 
2294   //    sample the medium parton thermal momentum in the comoving frame
2295 
2296   double xw;
2297   double razim;
2298   double rcos;
2299   double rsin;
2300 
2301   double ss;
2302   double tmin;
2303   double tmid;
2304   double tmax;
2305 
2306   double rant;
2307   double tt;
2308 
2309   double uu;
2310   double ff = 0.0;
2311   double rank;
2312 
2313   double mmax;
2314   double msq = 0.0;
2315 
2316   double f1;
2317   double f2;
2318 
2319   double p0ex[4] = {0.0};
2320 
2321   int ct1_loop, ct2_loop, flag1, flag2;
2322 
2323   flag1 = 0;
2324   flag2 = 0;
2325 
2326 
2327 
2328   //    Initial 4-momentum of jet
2329   //
2330   //************************************************************
2331   p4[1] = p0[1];
2332   p4[2] = p0[2];
2333   p4[3] = p0[3];
2334   p4[0] = p0[0];
2335   //************************************************************
2336 
2337   int ic = 0;
2338 
2339   ct1_loop = 0;
2340   do {
2341     ct1_loop++;
2342     if(flag2 == 1 || ct1_loop > 1e6){
2343        flag1 = 1;
2344        break;
2345     }
2346     ct2_loop = 0;
2347     do {
2348       ct2_loop++;
2349       if(ct2_loop > 1e6){
2350          flag2 = 1;
2351          break;
2352       }
2353       xw = 15.0 * ZeroOneDistribution(*GetMt19937Generator());
2354       razim = 2.0 * pi * ZeroOneDistribution(*GetMt19937Generator());
2355       rcos = 1.0 - 2.0 * ZeroOneDistribution(*GetMt19937Generator());
2356       rsin = sqrt(1.0 - rcos * rcos);
2357       //
2358       p2[0] = xw * temp;
2359       p2[3] = p2[0] * rcos;
2360       p2[1] = p2[0] * rsin * cos(razim);
2361       p2[2] = p2[0] * rsin * sin(razim);
2362 
2363       //
2364       //    cms energy
2365       //
2366       ss =
2367           2.0 * (p0[0] * p2[0] - p0[1] * p2[1] - p0[2] * p2[2] - p0[3] * p2[3]);
2368 
2369       //    if(ss.lt.2.d0*qhat0ud) goto 14
2370 
2371       tmin = qhat0ud;
2372       tmid = ss / 2.0;
2373       tmax = ss - qhat0ud;
2374 
2375       //    use (s^2+u^2)/(t+qhat0ud)^2 as scattering cross section in the
2376       //
2377       rant = ZeroOneDistribution(*GetMt19937Generator());
2378       tt = rant * ss;
2379 
2380       //        ic+=1;
2381       //        cout << p0[0] << "  " << p2[0] <<  endl;
2382       //        cout << tt << "  " << ss <<  "" << qhat0ud <<endl;
2383       //        cout << ic << endl;
2384 
2385     } while ((tt < qhat0ud) || (tt > (ss - qhat0ud)));
2386 
2387     f1 = pow(xw, 3) / (exp(xw) - 1) / 1.4215;
2388     f2 = pow(xw, 3) / (exp(xw) + 1) / 1.2845;
2389  
2390     uu = ss - tt;
2391 
2392     if (CT == 1) {
2393       ff = f1;
2394       mmax =
2395           4.0 / pow(ss, 2) *
2396           (3.0 - tmin * (ss - tmin) / pow(ss, 2) +
2397            (ss - tmin) * ss / pow(tmin, 2) + tmin * ss / pow((ss - tmin), 2));
2398       msq = pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
2399             (3.0 - tt * uu / pow(ss, 2) + uu * ss / pow(tt, 2) +
2400              tt * ss / pow(uu, 2)) /
2401             mmax;
2402     }
2403 
2404     if (CT == 2) {
2405       ff = f1;
2406       mmax = 4.0 / pow(ss, 2) *
2407              (4.0 / 9.0 * (pow(tmin, 2) + pow((ss - tmin), 2)) / tmin /
2408                   (ss - tmin) -
2409               (pow(tmin, 2) + pow((ss - tmin), 2)) / pow(ss, 2));
2410       msq = pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
2411             (4.0 / 9.0 * (pow(tt, 2) + pow(uu, 2)) / tt / uu -
2412              (pow(tt, 2) + pow(uu, 2)) / pow(ss, 2)) /
2413             (mmax + 4.0);
2414     }
2415 
2416     if (CT == 3) {
2417       ff = f2;
2418       if (((pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
2419            4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / ss / (ss - tmin)) >
2420           ((pow(ss, 2) + pow((ss - tmax), 2) / pow(tmax, 2) +
2421             4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmax), 2)) / ss /
2422                 (ss - tmax)))) {
2423         mmax =
2424             4.0 / pow(ss, 2) *
2425             ((pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
2426              4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / ss / (ss - tmin));
2427       } else {
2428         mmax =
2429             4.0 / pow(ss, 2) *
2430             ((pow(ss, 2) + pow((ss - tmax), 2)) / pow(tmax, 2) +
2431              4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmax), 2)) / ss / (ss - tmax));
2432       }
2433       //
2434       msq = pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
2435             ((pow(ss, 2) + pow(uu, 2)) / pow(tt, 2) +
2436              4.0 / 9.0 * (pow(ss, 2) + pow(uu, 2)) / ss / uu) /
2437             mmax;
2438     }
2439 
2440     if (CT == 13) {
2441       ff = f1;
2442 
2443       if (((pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
2444            4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / ss / (ss - tmin)) >
2445           ((pow(ss, 2) + pow((ss - tmax), 2) / pow(tmax, 2) +
2446             4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmax), 2)) / ss /
2447                 (ss - tmax)))) {
2448         mmax =
2449             4.0 / pow(ss, 2) *
2450             ((pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
2451              4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / ss / (ss - tmin));
2452       } else {
2453         mmax =
2454             4.0 / pow(ss, 2) *
2455             ((pow(ss, 2) + pow((ss - tmax), 2)) / pow(tmax, 2) +
2456              4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmax), 2)) / ss / (ss - tmax));
2457       }
2458       //
2459       msq = pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
2460             ((pow(ss, 2) + pow(uu, 2)) / pow(tt, 2) +
2461              4.0 / 9.0 * (pow(ss, 2) + pow(uu, 2)) / ss / uu) /
2462             mmax;
2463     }
2464 
2465     if (CT == 4) {
2466       ff = f2;
2467       mmax = 4.0 / pow(ss, 2) *
2468              (4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2));
2469       msq = pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
2470             (4.0 / 9.0 * (pow(ss, 2) + pow(uu, 2)) / pow(tt, 2)) / mmax;
2471     }
2472 
2473     if (CT == 5) {
2474       ff = f2;
2475       mmax = 4.0 / pow(ss, 2) *
2476              (4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
2477               (pow(ss, 2) + pow(tmin, 2)) / pow((ss - tmin), 2) -
2478               2.0 / 3.0 * pow(ss, 2) / tmin / (ss - tmin));
2479       msq = pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
2480             (4.0 / 9.0 *
2481              ((pow(ss, 2) + pow(uu, 2)) / pow(tt, 2) +
2482               (pow(ss, 2) + pow(tt, 2)) / pow(uu, 2) -
2483               2.0 / 3.0 * pow(ss, 2) / tt / uu)) /
2484             mmax;
2485     }
2486 
2487     if (CT == 6) {
2488       ff = f2;
2489       mmax = 4.0 / pow(ss, 2) *
2490              (4.0 / 9.0 * (pow(tmin, 2) + pow((ss - tmin), 2)) / pow(ss, 2));
2491       msq = pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
2492             (4.0 / 9.0 * (pow(tt, 2) + pow(uu, 2)) / pow(ss, 2)) / (mmax + 0.5);
2493     }
2494 
2495     if (CT == 7) {
2496       ff = f2;
2497       mmax = 4.0 / pow(ss, 2) *
2498              (4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
2499               (pow(tmin, 2) + pow((ss - tmin), 2)) / pow(ss, 2) +
2500               2.0 / 3.0 * pow((ss - tmin), 2) / ss / tmin);
2501       msq = (pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
2502              (4.0 / 9.0 *
2503               (((pow(ss, 2) + pow(uu, 2)) / pow(tt, 2)) +
2504                (pow(tt, 2) + pow(uu, 2)) / pow(ss, 2) +
2505                2.0 / 3.0 * pow(uu, 2) / ss / tt))) /
2506             mmax;
2507     }
2508 
2509     if (CT == 8) {
2510       ff = f2;
2511       mmax = 4.0 / pow(ss, 2) *
2512              (4.0 / 9.0 * (pow(tmin, 2) + pow((ss - tmin), 2)) / tmin /
2513                   (ss - tmin) -
2514               (pow(tmin, 2) + pow((ss - tmin), 2)) / pow(ss, 2));
2515       msq = pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
2516             (4.0 / 9.0 * (pow(tt, 2) + pow(uu, 2)) / tt / uu -
2517              (pow(tt, 2) + pow(uu, 2)) / pow(ss, 2)) /
2518             (mmax + 4.0);
2519     }
2520 
2521     rank = ZeroOneDistribution(*GetMt19937Generator());
2522   } while (rank > (msq * ff));
2523 
2524   if(flag1 == 1 || flag2 == 1){ // scatterings cannot be properly sampled
2525     transback(v0, p0);
2526     transback(v0, p4);
2527     qt = 0;
2528     p2[0] = 0;
2529     p2[1] = 0;
2530     p2[2] = 0;
2531     p2[3] = 0;
2532     p3[0] = 0;
2533     p3[1] = 0;
2534     p3[2] = 0;
2535     p3[3] = 0;
2536     return;
2537   }
2538 
2539   //
2540   p3[1] = p2[1];
2541   p3[2] = p2[2];
2542   p3[3] = p2[3];
2543   p3[0] = p2[0];
2544 
2545   //    velocity of the center-of-mass
2546   //
2547   vc[1] = (p0[1] + p2[1]) / (p0[0] + p2[0]);
2548   vc[2] = (p0[2] + p2[2]) / (p0[0] + p2[0]);
2549   vc[3] = (p0[3] + p2[3]) / (p0[0] + p2[0]);
2550   //
2551   //    transform into the cms frame
2552   //
2553   trans(vc, p0);
2554   trans(vc, p2);
2555   //
2556   //    cm momentum
2557   //
2558   double pcm = p2[0];
2559   //
2560   //    sample transverse momentum transfer with respect to jet momentum
2561   //    in cm frame
2562   //
2563   double ranp = 2.0 * pi * ZeroOneDistribution(*GetMt19937Generator());
2564   //
2565   //    transverse momentum transfer
2566   //
2567   qt = sqrt(pow(pcm, 2) - pow((tt / 2.0 / pcm - pcm), 2));
2568   double qx = qt * cos(ranp);
2569   double qy = qt * sin(ranp);
2570 
2571   //
2572   //    longitudinal momentum transfer
2573   //
2574   double qpar = tt / 2.0 / pcm;
2575   //
2576   //    qt is perpendicular to pcm, need to rotate back to the cm frame
2577   //
2578   double upt = sqrt(p2[1] * p2[1] + p2[2] * p2[2]) / p2[0];
2579   double upx = p2[1] / p2[0];
2580   double upy = p2[2] / p2[0];
2581   double upz = p2[3] / p2[0];
2582   //
2583   //    momentum after collision in cm frame
2584   //
2585   p2[1] = p2[1] - qpar * upx;
2586   p2[2] = p2[2] - qpar * upy;
2587   if (upt != 0.0) {
2588     p2[1] = p2[1] + (upz * upx * qy + upy * qx) / upt;
2589     p2[2] = p2[2] + (upz * upy * qy - upx * qx) / upt;
2590   }
2591   p2[3] = p2[3] - qpar * upz - upt * qy;
2592 
2593   p0[1] = -p2[1];
2594   p0[2] = -p2[2];
2595   p0[3] = -p2[3];
2596   //
2597   //    transform from cm back to the comoving frame
2598   //
2599   transback(vc, p2);
2600   transback(vc, p0);
2601 
2602   //************************************************************
2603   //
2604   //     calculate qt in the rest frame of medium
2605   //
2606   //  if(p0[4]>p2[4])
2607   //    {
2608   rotate(p4[1], p4[2], p4[3], p0, 1);
2609   qt = sqrt(pow(p0[1], 2) + pow(p0[2], 2));
2610   rotate(p4[1], p4[2], p4[3], p0, -1);
2611   //    }
2612   //  else
2613   //    {
2614   //      rotate(p4[1],p4[2],p4[3],p2,1);
2615   //      qt=sqrt(pow(p2[1],2)+pow(p2[2],2));
2616   //      rotate(p4[1],p4[2],p4[3],p2,-1);
2617   //    }
2618   //************************************************************
2619 
2620   //
2621   //    transform from comoving frame to the lab frame
2622   //
2623   transback(v0, p2);
2624   transback(v0, p0);
2625   transback(v0, p3);
2626 
2627   //************************************************************
2628   transback(v0, p4);
2629   //************************************************************
2630 }
2631 
2632 //.........................................................................
2633 void LBT::collHQ22(int CT, double temp, double qhat0ud, double v0[4],
2634                    double p0[4], double p2[4], double p3[4], double p4[4],
2635                    double &qt) {
2636   //
2637   //    HQ 2->2 scatterings
2638   //    p0 initial HQ momentum, output to final momentum
2639   //    p2 final thermal momentum, p3 initial thermal energy
2640   //
2641   //    amss=sqrt(abs(p0(4)**2-p0(1)**2-p0(2)**2-p0(3)**2))
2642   //
2643   //************************************************************
2644 
2645   // transform to local comoving frame of the fluid
2646   trans(v0, p0);
2647 
2648   //************************************************************
2649 
2650   //    sample the medium parton thermal momentum in the comoving frame
2651 
2652   double xw;
2653   double razim;
2654   double rcos;
2655   double rsin;
2656 
2657   double ss;
2658 
2659   double rant;
2660   double tt;
2661 
2662   double uu;
2663   double ff = 0.0;
2664   double rank;
2665 
2666   double msq = 0.0;
2667 
2668   double e2, theta2, theta4, phi24; // the four independent variables
2669   double e1, e4, p1, cosTheta24, downFactor,
2670       sigFactor; // other useful variables
2671   double HQmass, fBmax, fFmax, fB, fF, maxValue;
2672   int index_p1, index_T, index_e2;
2673   int ct1_loop, ct2_loop, flag1, flag2;
2674 
2675   flag1 = 0;
2676   flag2 = 0;
2677 
2678   // continue this function for HQ scattering
2679 
2680   HQmass = p0[0] * p0[0] - p0[1] * p0[1] - p0[2] * p0[2] - p0[3] * p0[3];
2681   if (HQmass > 1e-12) {
2682     HQmass = sqrt(HQmass);
2683   } else {
2684     HQmass = 0.0;
2685   }
2686 
2687   //    Initial 4-momentum of HQ
2688   //
2689   //************************************************************
2690   p4[1] = p0[1];
2691   p4[2] = p0[2];
2692   p4[3] = p0[3];
2693   p4[0] = p0[0];
2694   //************************************************************
2695 
2696   p1 = sqrt(p0[1] * p0[1] + p0[2] * p0[2] + p0[3] * p0[3]);
2697   index_p1 = (int)((p1 - min_p1) / bin_p1);
2698   index_T = (int)((temp - min_T) / bin_T);
2699   if (index_p1 >= N_p1) {
2700     index_p1 = N_p1 - 1;
2701     cout << "warning: p1 is over p_max: " << p1 << endl;
2702   }
2703   if (index_T >= N_T) {
2704     index_T = N_T - 1;
2705     cout << "warning: T is over T_max: " << temp << endl;
2706   }
2707   if (index_T < 0) {
2708     index_T = 0;
2709     cout << "warning: T is below T_min: " << temp << endl;
2710   }
2711 
2712   fBmax = distFncBM[index_T][index_p1];
2713   fFmax = distFncFM[index_T][index_p1]; // maximum of f(xw) at given p1 and T
2714 
2715   maxValue = 10.0; // need actual value later
2716 
2717   ct1_loop = 0;
2718   do { // sample p2 (light parton) using distribution integrated over 3 angles
2719     ct1_loop++;
2720     if (ct1_loop > 1e6) {
2721       //            cout << "cannot sample light parton for HQ scattering ..." << endl;
2722       flag1 = 1;
2723       break;
2724     }
2725     xw = max_e2 * ZeroOneDistribution(*GetMt19937Generator());
2726     index_e2 = (int)((xw - min_e2) / bin_e2);
2727     if (index_e2 >= N_e2)
2728       index_e2 = N_e2 - 1;
2729     if (CT == 11) { // qc->qc
2730       ff = distFncF[index_T][index_p1][index_e2] / fFmax;
2731       maxValue = distMaxF[index_T][index_p1][index_e2];
2732     } else if (CT == 12) { // gc->gc
2733       ff = distFncB[index_T][index_p1][index_e2] / fBmax;
2734       maxValue = distMaxB[index_T][index_p1][index_e2];
2735     } else {
2736       cout << "Wrong HQ channel ID" << endl;
2737       exit(EXIT_FAILURE);
2738     }
2739   } while (ZeroOneDistribution(*GetMt19937Generator()) > ff);
2740 
2741   e2 = xw * temp;
2742   e1 = p0[0];
2743 
2744   // now e2 is fixed, need to sample the remaining 3 variables
2745   ct2_loop = 0;
2746   do {
2747     ct2_loop++;
2748     if (ct2_loop > 1e6) {
2749       cout << "cannot sample final states for HQ scattering ..." << endl;
2750       flag2 = 1;
2751       break;
2752     }
2753 
2754     theta2 = pi * ZeroOneDistribution(*GetMt19937Generator());
2755     theta4 = pi * ZeroOneDistribution(*GetMt19937Generator());
2756     phi24 = 2.0 * pi * ZeroOneDistribution(*GetMt19937Generator());
2757 
2758     cosTheta24 =
2759         sin(theta2) * sin(theta4) * cos(phi24) + cos(theta2) * cos(theta4);
2760     downFactor = e1 - p1 * cos(theta4) + e2 - e2 * cosTheta24;
2761     e4 = (e1 * e2 - p1 * e2 * cos(theta2)) / downFactor;
2762     sigFactor = sin(theta2) * sin(theta4) * e2 * e4 / downFactor;
2763 
2764     // calculate s,t,u, different definition from light quark -- tt, uu are negative
2765     ss = 2.0 * e1 * e2 + HQmass * HQmass - 2.0 * p1 * e2 * cos(theta2);
2766     tt = -2.0 * e2 * e4 * (1.0 - cosTheta24);
2767     uu = 2.0 * HQmass * HQmass - ss - tt;
2768 
2769     // re-sample if the kinematic cuts are not satisfied
2770     if (ss <= 2.0 * qhat0ud || tt >= -qhat0ud || uu >= -qhat0ud) {
2771       rank = ZeroOneDistribution(*GetMt19937Generator());
2772       sigFactor = 0.0;
2773       msq = 0.0;
2774       continue;
2775     }
2776 
2777     if (CT == 11) { // qc->qc
2778       ff =
2779           (1.0 / (exp(e2 / temp) + 1.0)) * (1.0 - 1.0 / (exp(e4 / temp) + 1.0));
2780       sigFactor = sigFactor * ff;
2781       msq = Mqc2qc(ss, tt, HQmass) / maxValue;
2782     }
2783 
2784     if (CT == 12) { // gc->gc
2785       ff =
2786           (1.0 / (exp(e2 / temp) - 1.0)) * (1.0 + 1.0 / (exp(e4 / temp) - 1.0));
2787       sigFactor = sigFactor * ff;
2788       msq = Mgc2gc(ss, tt, HQmass) / maxValue;
2789     }
2790 
2791     rank = ZeroOneDistribution(*GetMt19937Generator());
2792 
2793   } while (rank > (msq * sigFactor));
2794 
2795   if (flag1 == 0 && flag2 == 0) {
2796 
2797     // pass p2 value to p3 for initial thermal parton
2798     p3[1] = e2 * sin(theta2);
2799     p3[2] = 0.0;
2800     p3[3] = e2 * cos(theta2);
2801     p3[0] = e2;
2802 
2803     // calculate momenta of outgoing particles
2804     // here p2 is for p4 (light parton) in my note
2805 
2806     p2[1] = e4 * sin(theta4) * cos(phi24);
2807     p2[2] = e4 * sin(theta4) * sin(phi24);
2808     p2[3] = e4 * cos(theta4);
2809     p2[0] = e4;
2810 
2811     // rotate randomly in xy plane (jet is in z), because p3 is assigned in xz plane with bias
2812     double th_rotate = 2.0 * pi * ZeroOneDistribution(*GetMt19937Generator());
2813     double p3x_rotate = p3[1] * cos(th_rotate) - p3[2] * sin(th_rotate);
2814     double p3y_rotate = p3[1] * sin(th_rotate) + p3[2] * cos(th_rotate);
2815     double p2x_rotate = p2[1] * cos(th_rotate) - p2[2] * sin(th_rotate);
2816     double p2y_rotate = p2[1] * sin(th_rotate) + p2[2] * cos(th_rotate);
2817     p3[1] = p3x_rotate;
2818     p3[2] = p3y_rotate;
2819     p2[1] = p2x_rotate;
2820     p2[2] = p2y_rotate;
2821 
2822     // Because we treated p0 (p1 in my note for heavy quark) as the z-direction, proper rotations are necessary here
2823     rotate(p4[1], p4[2], p4[3], p2, -1);
2824     rotate(p4[1], p4[2], p4[3], p3, -1);
2825 
2826     p0[1] = p4[1] + p3[1] - p2[1];
2827     p0[2] = p4[2] + p3[2] - p2[2];
2828     p0[3] = p4[3] + p3[3] - p2[3];
2829     p0[0] =
2830         sqrt(p0[1] * p0[1] + p0[2] * p0[2] + p0[3] * p0[3] + HQmass * HQmass);
2831 
2832     // Debug
2833     if (fabs(p0[0] + p2[0] - p3[0] - p4[0]) > 0.00001) {
2834       cout << "Violation of energy conservation in HQ 2->2 scattering:  "
2835            << fabs(p0[0] + p2[0] - p3[0] - p4[0]) << endl;
2836     }
2837 
2838     // calculate qt in the rest frame of medium
2839     rotate(p4[1], p4[2], p4[3], p0, 1);
2840     qt = sqrt(pow(p0[1], 2) + pow(p0[2], 2));
2841     rotate(p4[1], p4[2], p4[3], p0, -1);
2842 
2843     // transform from comoving frame to the lab frame
2844     transback(v0, p2);
2845     transback(v0, p0);
2846     transback(v0, p3);
2847     transback(v0, p4);
2848 
2849   } else { // no scattering
2850     transback(v0, p0);
2851     transback(v0, p4);
2852     qt = 0;
2853     p2[0] = 0;
2854     p2[1] = 0;
2855     p2[2] = 0;
2856     p2[3] = 0;
2857     p3[0] = 0;
2858     p3[1] = 0;
2859     p3[2] = 0;
2860     p3[3] = 0;
2861   }
2862 }
2863 
2864 void LBT::twcoll(int CT, double qhat0ud, double v0[4], double p0[4],
2865                  double p2[4]) {
2866   //
2867   //     p0 initial jet momentum, output to final momentum
2868   //     p2 final thermal momentum,p3 initial thermal energy
2869   //
2870   //     amss=sqrt(abs(p0(4)**2-p0(1)**2-p0(2)**2-p0(3)**2))
2871   //
2872   //     transform to local comoving frame of the fluid
2873 
2874   trans(v0, p0);
2875   trans(v0, p2);
2876 
2877   //    velocity of the center-of-mass
2878 
2879   vc[1] = (p0[1] + p2[1]) / (p0[0] + p2[0]);
2880   vc[2] = (p0[2] + p2[2]) / (p0[0] + p2[0]);
2881   vc[3] = (p0[3] + p2[3]) / (p0[0] + p2[0]);
2882   //
2883   //    transform into the cms frame
2884   //
2885 
2886   trans(vc, p0);
2887   trans(vc, p2);
2888 
2889   //
2890   //     cm momentum
2891   //
2892   double pcm = p2[0];
2893   //
2894   //     sample transverse momentum transfer with respect to jet momentum
2895   //     in cm frame
2896   //
2897   //
2898   //     Gaussian distribution
2899   //
2900   //     qt=sqrt(-dt*qhat0ud*log(1-rant+rant*exp(-scm/(4.d0*dt*qhat0ud))))
2901   //
2902   //     static potential distribution
2903   //
2904   double ss = 4.0 * pow(pcm, 2);
2905   //
2906   double tmin = qhat0ud;
2907   double tmid = ss / 2.0;
2908   double tmax = ss - qhat0ud;
2909 
2910   double rant;
2911   double tt;
2912   double uu;
2913   double mmax;
2914   double msq = 0.0;
2915   double rank;
2916 
2917   /////////////////////////////////////////////
2918   double ranp;
2919   double qt;
2920   double qx;
2921   double qy;
2922   double qpar;
2923   double upt;
2924   double upx;
2925   double upy;
2926   double upz;
2927   /////////////////////////////////////////////
2928 
2929   //
2930   //       CT is a variable notated different collision types.
2931   //
2932 
2933   do {
2934     rant = ZeroOneDistribution(*GetMt19937Generator());
2935     tt = rant * ss;
2936 
2937     if ((tt < qhat0ud) || (tt > (ss - qhat0ud)))
2938       break;
2939     //      if((tt<qhat0ud) || (tt>(ss-qhat0ud)))
2940     //      {
2941     //      goto t59;
2942     //      }
2943     uu = ss - tt;
2944 
2945     //     gg to gg
2946     if (CT == 1) {
2947       //
2948       mmax = 3.0 - tmin * (ss - tmin) / pow(ss, 2) +
2949              (ss - tmin) * ss / pow(tmin, 2) + tmin * ss / pow((ss - tmin), 2);
2950       msq = (3.0 - tt * uu / pow(ss, 2) + uu * ss / pow(tt, 2) +
2951              tt * ss / pow(uu, 2)) /
2952             mmax;
2953       //
2954     }
2955 
2956     //     gg to qqbar
2957     if (CT == 2) {
2958       //
2959       mmax = 4.0 / 9.0 * (pow(tmin, 2) + pow((ss - tmin), 2)) / tmin /
2960                  (ss - tmin) -
2961              (pow(tmin, 2) + pow((ss - tmin), 2)) / pow(ss, 2);
2962       msq = (4.0 / 9.0 * (pow(tt, 2) + pow(uu, 2)) / tt / uu -
2963              (pow(tt, 2) + pow(uu, 2)) / pow(ss, 2)) /
2964             (mmax + 4.0);
2965       //
2966     }
2967 
2968     //     gq to gq, gqbar to gqbar
2969     if (CT == 3) {
2970       //
2971       if (((pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
2972            4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / ss / (ss - tmin)) >
2973           ((pow(ss, 2) + pow((ss - tmax), 2) / pow(tmax, 2) +
2974             4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmax), 2)) / ss /
2975                 (ss - tmax)))) {
2976         mmax =
2977             (pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
2978             4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / ss / (ss - tmin);
2979       } else {
2980         mmax =
2981             4.0 / pow(ss, 2) *
2982             ((pow(ss, 2) + pow((ss - tmax), 2)) / pow(tmax, 2) +
2983              4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmax), 2)) / ss / (ss - tmax));
2984       }
2985       //
2986       msq = ((pow(ss, 2) + pow(uu, 2)) / pow(tt, 2) +
2987              4.0 / 9.0 * (pow(ss, 2) + pow(uu, 2)) / ss / uu) /
2988             mmax;
2989       //
2990     }
2991 
2992     //     qg to qg, qbarg to qbarg
2993     if (CT == 13) {
2994       //
2995       if (((pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
2996            4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / ss / (ss - tmin)) >
2997           ((pow(ss, 2) + pow((ss - tmax), 2) / pow(tmax, 2) +
2998             4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmax), 2)) / ss /
2999                 (ss - tmax)))) {
3000         mmax =
3001             (pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
3002             4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / ss / (ss - tmin);
3003       } else {
3004         mmax =
3005             4.0 / pow(ss, 2) *
3006             ((pow(ss, 2) + pow((ss - tmax), 2)) / pow(tmax, 2) +
3007              4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmax), 2)) / ss / (ss - tmax));
3008       }
3009       //
3010       msq = ((pow(ss, 2) + pow(uu, 2)) / pow(tt, 2) +
3011              4.0 / 9.0 * (pow(ss, 2) + pow(uu, 2)) / ss / uu) /
3012             mmax;
3013       //
3014     }
3015 
3016     //     qiqj to qiqj, qiqjbar to qiqjbar, qibarqj to qibarqj, qibarqjbar to qibarqjbar
3017     //     for i not equal j
3018     if (CT == 4) {
3019       //
3020       mmax = 4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2);
3021       msq = (4.0 / 9.0 * (pow(ss, 2) + pow(uu, 2)) / pow(tt, 2)) / mmax;
3022       //
3023     }
3024 
3025     //     qiqi to qiqi, qibarqibar to qibarqibar
3026     if (CT == 5) {
3027       //
3028       mmax = 4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
3029              (pow(ss, 2) + pow(tmin, 2)) / pow((ss - tmin), 2) -
3030              2.0 / 3.0 * pow(ss, 2) / tmin / (ss - tmin);
3031       msq = (4.0 / 9.0 *
3032              ((pow(ss, 2) + pow(uu, 2)) / pow(tt, 2) +
3033               (pow(ss, 2) + pow(tt, 2)) / pow(uu, 2) -
3034               2.0 / 3.0 * pow(ss, 2) / tt / uu)) /
3035             mmax;
3036       //
3037     }
3038 
3039     //     qiqibar to qjqjbar for i not equal j
3040     if (CT == 6) {
3041       //
3042       mmax = 4.0 / 9.0 * (pow(tmin, 2) + pow((ss - tmin), 2)) / pow(ss, 2);
3043       msq = (4.0 / 9.0 * (pow(tt, 2) + pow(uu, 2)) / pow(ss, 2)) / (mmax + 0.5);
3044       //
3045     }
3046 
3047     //     qiqibar to qiqibar
3048     if (CT == 7) {
3049       //
3050       mmax = 4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
3051              (pow(tmin, 2) + pow((ss - tmin), 2)) / pow(ss, 2) +
3052              2.0 / 3.0 * pow((ss - tmin), 2) / ss / tmin;
3053       msq = (4.0 / 9.0 *
3054              (((pow(ss, 2) + pow(uu, 2)) / pow(tt, 2)) +
3055               (pow(tt, 2) + pow(uu, 2)) / pow(ss, 2) +
3056               2.0 / 3.0 * pow(uu, 2) / ss / tt)) /
3057             mmax;
3058       //
3059     }
3060 
3061     //     qqbar to gg
3062     if (CT == 8) {
3063       //
3064       mmax = 4.0 / 9.0 * (pow(tmin, 2) + pow((ss - tmin), 2)) / tmin /
3065                  (ss - tmin) -
3066              (pow(tmin, 2) + pow((ss - tmin), 2)) / pow(ss, 2);
3067       msq = (4.0 / 9.0 * (pow(tt, 2) + pow(uu, 2)) / tt / uu -
3068              (pow(tt, 2) + pow(uu, 2)) / pow(ss, 2)) /
3069             (mmax + 4.0);
3070       //
3071     }
3072 
3073     rank = ZeroOneDistribution(*GetMt19937Generator());
3074 
3075   } while (rank > msq);
3076 
3077   ///////////////////////////////////////////////////////////////////////
3078   //    transverse momentum transfer
3079 
3080   if ((tt > qhat0ud) && (tt < (ss - qhat0ud))) {
3081 
3082     ranp = 2.0 * pi * ZeroOneDistribution(*GetMt19937Generator());
3083     //
3084     //
3085     //
3086     qt = sqrt(pow(pcm, 2) - pow((tt / 2.0 / pcm - pcm), 2));
3087     qx = qt * cos(ranp);
3088     qy = qt * sin(ranp);
3089 
3090     //
3091     //    longitudinal momentum transfer
3092     //
3093     qpar = tt / 2.0 / pcm;
3094     //
3095     //    qt is perpendicular to pcm, need to rotate back to the cm frame
3096     //
3097     upt = sqrt(p2[1] * p2[1] + p2[2] * p2[2]) / p2[0];
3098     upx = p2[1] / p2[0];
3099     upy = p2[2] / p2[0];
3100     upz = p2[3] / p2[0];
3101     //
3102     //    momentum after collision in cm frame
3103     //
3104 
3105     p2[1] = p2[1] - qpar * upx;
3106     p2[2] = p2[2] - qpar * upy;
3107     if (upt != 0.0) {
3108       p2[1] = p2[1] + (upz * upx * qy + upy * qx) / upt;
3109       p2[2] = p2[2] + (upz * upy * qy - upx * qx) / upt;
3110     }
3111     p2[3] = p2[3] - qpar * upz - upt * qy;
3112 
3113     p0[1] = -p2[1];
3114     p0[2] = -p2[2];
3115     p0[3] = -p2[3];
3116   }
3117   //
3118   //    transform from cm back to the comoving frame
3119   //
3120   transback(vc, p2);
3121   transback(vc, p0);
3122   //
3123   //    transform from comoving frame to the lab frame
3124   //
3125 
3126   transback(v0, p2);
3127   transback(v0, p0);
3128   transback(v0, p3);
3129 }
3130 
3131 void LBT::collHQ23(int parID, double temp_med, double qhat0ud, double v0[4],
3132                    double p0[4], double p2[4], double p3[4], double p4[4],
3133                    double qt, int &ic, double Tdiff, double HQenergy,
3134                    double max_Ng, double xLow, double xInt) {
3135 
3136   //    p0 initial jet momentum, output to final momentum
3137   //    p3 initial thermal momentum
3138   //    p2 initial thermal momentum, output to final thermal momentum
3139   //    p4 radiated gluon momentum
3140   //    qt transverse momentum transfer in the rest frame of medium
3141   //    q0,ql energy and longitudinal momentum transfer
3142   //    i=0: 2->3 finished; i=1: can not find a gluon when nloop<=30
3143   //
3144   //    amss=sqrt(abs(p0(4)**2-p0(1)**2-p0(2)**2-p0(3)**2))
3145 
3146   double randomX, randomY;
3147   int count_sample, nloopOut, flagOut;
3148   double theta_gluon, kperp_gluon;
3149   double kpGluon[4];
3150   double zDirection[4];
3151   double HQmass =
3152       sqrt(p0[0] * p0[0] - p0[1] * p0[1] - p0[2] * p0[2] - p0[3] * p0[3]);
3153 
3154   double p00[4];
3155   p00[0] = p0[0];
3156   p00[1] = p0[1];
3157   p00[2] = p0[2];
3158   p00[3] = p0[3];
3159 
3160   if (abs(parID) != 4 && abs(parID) != 5) {
3161     HQmass = 0.0;
3162     p0[0] =
3163         sqrt(p0[1] * p0[1] + p0[2] * p0[2] + p0[3] * p0[3] + HQmass * HQmass);
3164   }
3165 
3166   ic = 0;
3167   nloopOut = 0;
3168   flagOut = 0;
3169 
3170   //    initial thermal parton momentum in lab frame
3171   p2[1] = p3[1];
3172   p2[2] = p3[2];
3173   p2[3] = p3[3];
3174   p2[0] = p3[0];
3175 
3176   //    transform to local comoving frame of the fluid
3177   trans(v0, p0);
3178   trans(v0, p2);
3179 
3180   if (p0[0] < 2.0 * sqrt(qhat0ud)) {
3181     ic = 1;
3182     return;
3183   }
3184 
3185   // define z-direction
3186   zDirection[1] = p0[1];
3187   zDirection[2] = p0[2];
3188   zDirection[3] = p0[3];
3189   zDirection[0] = p0[0];
3190 
3191   rotate(zDirection[1], zDirection[2], zDirection[3], p2,
3192          1); // rotate p2 into p1 system
3193 
3194   // comment do { for unit test
3195   do {
3196 
3197     do {
3198       randomX = xLow + xInt * ZeroOneDistribution(*GetMt19937Generator());
3199       randomY = ZeroOneDistribution(*GetMt19937Generator());
3200     } while (tau_f(randomX, randomY, HQenergy, HQmass) < 1.0 / pi / temp_med);
3201 
3202     count_sample = 0;
3203     while (max_Ng * ZeroOneDistribution(*GetMt19937Generator()) > dNg_over_dxdydt(parID, randomX, randomY,
3204                                                   HQenergy, HQmass, temp_med,
3205                                                   Tdiff)) {
3206       count_sample = count_sample + 1;
3207       if (count_sample > 1e+5) {
3208         //cout << "give up loop at point 1 ..." << endl;
3209         kpGluon[1] = 0.0;
3210         kpGluon[2] = 0.0;
3211         kpGluon[3] = 0.0;
3212         kpGluon[0] = 0.0;
3213         ic = 1;
3214         //          break;
3215         return;
3216       }
3217 
3218       do {
3219         randomX = xLow + xInt * ZeroOneDistribution(*GetMt19937Generator());
3220         randomY = ZeroOneDistribution(*GetMt19937Generator());
3221       } while (tau_f(randomX, randomY, HQenergy, HQmass) < 1.0 / pi / temp_med);
3222     }
3223 
3224     if (parID == 21 && randomX > 0.5)
3225       randomX = 1.0 - randomX;
3226     theta_gluon = 2.0 * pi * ZeroOneDistribution(*GetMt19937Generator());
3227     kperp_gluon = randomX * randomY * HQenergy;
3228     kpGluon[1] = kperp_gluon * cos(theta_gluon);
3229     kpGluon[2] = kperp_gluon * sin(theta_gluon);
3230     kpGluon[3] = randomX * HQenergy * sqrt(1.0 - randomY * randomY);
3231     kpGluon[0] = sqrt(kpGluon[1] * kpGluon[1] + kpGluon[2] * kpGluon[2] +
3232                       kpGluon[3] * kpGluon[3]);
3233 
3234     // comment for unit test
3235     if (kpGluon[0] > (HQenergy - HQmass)) {
3236       kpGluon[1] = 0.0;
3237       kpGluon[2] = 0.0;
3238       kpGluon[3] = 0.0;
3239       kpGluon[0] = 0.0;
3240       nloopOut++;
3241       continue;
3242     }
3243 
3244     // update mometum
3245 
3246     p4[1] = kpGluon[1];
3247     p4[2] = kpGluon[2];
3248     p4[3] = kpGluon[3];
3249     p4[0] = kpGluon[0];
3250 
3251     //comment for unit test begin
3252     // solve energy-momentum conservation
3253     double sE1 = p0[0];
3254     double sp1x = 0.0;
3255     double sp1y = 0.0;
3256     double sp1z = sqrt(p0[1] * p0[1] + p0[2] * p0[2] + p0[3] * p0[3]);
3257     double sE2 = p2[0];
3258     double sp2x = p2[1];
3259     double sp2y = p2[2];
3260     double sp2z = p2[3];
3261     double sk0 = p4[0];
3262     double skx = p4[1];
3263     double sky = p4[2];
3264     double skz = p4[3];
3265     double sqx, sqy, sqzA, sq0A, sqzB, sq0B, sqz, sq0, sqtheta;
3266     double sAA, sBB, sCC, aaa, bbb, ccc, abc, abc2;
3267     int nloop1 = 0;
3268     int nloop2 = 0;
3269     int flagDone = 0;
3270     int yesA, yesB;
3271 
3272     do {
3273       sqtheta = 2.0 * pi * ZeroOneDistribution(*GetMt19937Generator());
3274       sqx = qt * cos(sqtheta);
3275       sqy = qt * sin(sqtheta);
3276       sAA = (sE1 + sE2 - sk0) / (sp1z + sp2z - skz);
3277       sBB =
3278           (pow(sE2, 2) - pow((sE1 - sk0), 2) + pow((sp1x - sqx - skx), 2) +
3279            pow((sp1y - sqy - sky), 2) + pow((sp1z - skz), 2) + pow(HQmass, 2) -
3280            pow((sp2x + sqx), 2) - pow((sp2y + sqy), 2) - pow(sp2z, 2)) /
3281           2.0 / (sp1z + sp2z - skz);
3282       aaa = sAA * sAA - 1.0;
3283       bbb = 2.0 * (sAA * sp2z + sAA * sBB - sE2);
3284       ccc = pow((sp2x + sqx), 2) + pow((sp2y + sqy), 2) + sp2z * sp2z +
3285             2.0 * sp2z * sBB + sBB * sBB - sE2 * sE2;
3286       abc2 = bbb * bbb - 4.0 * aaa * ccc;
3287 
3288       if (abc2 < 0.0) {
3289         nloop1++;
3290         continue;
3291       } else {
3292         nloop1 = 0;
3293       }
3294 
3295       abc = sqrt(abc2);
3296       sq0A = (-bbb + abc) / 2.0 / aaa;
3297       sq0B = (-bbb - abc) / 2.0 / aaa;
3298       sqzA = sAA * sq0A + sBB;
3299       sqzB = sAA * sq0B + sBB;
3300 
3301       // require space-like and E_final > M;
3302       if (sq0A * sq0A - sqx * sqx - sqy * sqy - sqzA * sqzA < 0 &&
3303           sE1 - sq0A - sk0 > HQmass && sE2 + sq0A > 0) {
3304         yesA = 1;
3305       } else {
3306         yesA = 0;
3307       }
3308       if (sq0B * sq0B - sqx * sqx - sqy * sqy - sqzB * sqzB < 0 &&
3309           sE1 - sq0B - sk0 > HQmass && sE2 + sq0B > 0) {
3310         yesB = 1;
3311       } else {
3312         yesB = 0;
3313       }
3314 
3315       // select appropriate solution
3316       if (yesA == 0 && yesB == 0) {
3317         //          cout << "Solutions fail ..." << endl;
3318         nloop2++;
3319         continue;
3320       } else if (yesA == 1 && yesB == 1) {
3321         //          cout << "Both solutions work!" << "  " << sE1 << "  " << sp1x << "  " << sp1y << "  " << sp1z << "  " << sE2 << "  " << sp2x << "  " << sp2y << "  " << sp2z << "  " << sk0 << "  " << skx << "  " << sky << "  " << skz << "  " << sqx << "  " << sqy << "  " << sq0A << "  " << sqzA << "  " << sq0B << "  " << sqzB << endl;
3322         if (abs(sq0A) < abs(sq0B)) {
3323           sq0 = sq0A;
3324           sqz = sqzA;
3325         } else {
3326           sq0 = sq0B;
3327           sqz = sqzB;
3328         }
3329       } else if (yesA == 1) {
3330         //          cout << "pass A ..." << "  " << sE1 << "  " << sp1x << "  " << sp1y << "  " << sp1z << "  " << sE2 << "  " << sp2x << "  " << sp2y << "  " << sp2z << "  " << sk0 << "  " << skx << "  " << sky << "  " << skz << "  " << sqx << "  " << sqy << "  " << sq0A << "  " << sqzA << "  " << sq0B << "  " << sqzB << endl;
3331         sq0 = sq0A;
3332         sqz = sqzA;
3333       } else {
3334         //          cout << "pass B ..." << "  " << sE1 << "  " << sp1x << "  " << sp1y << "  " << sp1z << "  " << sE2 << "  " << sp2x << "  " << sp2y << "  " << sp2z << "  " << sk0 << "  " << skx << "  " << sky << "  " << skz << "  " << sqx << "  " << sqy << "  " << sq0A << "  " << sqzA << "  " << sq0B << "  " << sqzB << endl;
3335         sq0 = sq0B;
3336         sqz = sqzB;
3337       }
3338 
3339       flagDone = 1;
3340     } while (flagDone == 0 && nloop1 < loopN && nloop2 < loopN);
3341 
3342     if (flagDone == 0) {
3343       //      ic=1;  // no appropriate solution
3344       //      cout << "solution fails ..." << endl;
3345       nloopOut++;
3346       continue;
3347     } else {
3348       p0[0] = sE1 - sq0 - sk0;
3349       p0[1] = sp1x - sqx - skx;
3350       p0[2] = sp1y - sqy - sky;
3351       p0[3] = sp1z - sqz - skz;
3352 
3353       p2[0] = sE2 + sq0;
3354       p2[1] = sp2x + sqx;
3355       p2[2] = sp2y + sqy;
3356       p2[3] = sp2z + sqz;
3357 
3358       //           cout<<"sE1,sE2,p0[0],p2[0],p4[0]: "<<sE1<<"  "<<sE2<<"  "<<p0[0]<<"  "<<p2[0]<<"  "<<p4[0]<<endl;
3359       // comment for unit test end
3360 
3361       rotate(zDirection[1], zDirection[2], zDirection[3], p0,
3362              -1); // rotate p0 into global system
3363       rotate(zDirection[1], zDirection[2], zDirection[3], p2,
3364              -1); // rotate p2 into global system
3365       rotate(zDirection[1], zDirection[2], zDirection[3], p4,
3366              -1); // rotate p4 into global system
3367 
3368       transback(v0, p0);
3369       transback(v0, p2);
3370       transback(v0, p4);
3371 
3372       if (p00[0] + p3[0] - p0[0] - p2[0] - p4[0] > 0.01) {
3373         JSINFO << MAGENTA << "Violation of energy-momentum conservation: "
3374                << p00[0] + p3[0] - p0[0] - p2[0] - p4[0] << "  "
3375                << p00[1] + p3[1] - p0[1] - p2[1] - p4[1] << "  "
3376                << p00[2] + p3[2] - p0[2] - p2[2] - p4[2] << "  "
3377                << p00[3] + p3[3] - p0[3] - p2[3] - p4[3];
3378       }
3379 
3380       // comment for unit test begin
3381       // debug: check on-shell condition
3382       if (abs(p0[0] * p0[0] - p0[1] * p0[1] - p0[2] * p0[2] - p0[3] * p0[3] -
3383               HQmass * HQmass) > 0.01 ||
3384           abs(p2[0] * p2[0] - p2[1] * p2[1] - p2[2] * p2[2] - p2[3] * p2[3]) >
3385               0.01) {
3386         cout << "Wrong solution -- not on shell"
3387              << "  " << sE1 << "  " << sp1x << "  " << sp1y << "  " << sp1z
3388              << "  " << sE2 << "  " << sp2x << "  " << sp2y << "  " << sp2z
3389              << "  " << sk0 << "  " << skx << "  " << sky << "  " << skz << "  "
3390              << sq0 << "  " << sqx << "  " << sqy << "  " << sqz << endl;
3391         cout << abs(p0[0] * p0[0] - p0[1] * p0[1] - p0[2] * p0[2] -
3392                     p0[3] * p0[3] - HQmass * HQmass)
3393              << "  "
3394              << abs(p2[0] * p2[0] - p2[1] * p2[1] - p2[2] * p2[2] -
3395                     p2[3] * p2[3])
3396              << "  " << HQmass << endl;
3397       }
3398 
3399       //  cout<<"in colljet23 -- init thermal, final jet, final thermal, gluon: "<<p3[0]<<"  "<<p0[0]<<"  "<<p2[0]<<"  "<<p4[0]<<endl;
3400       flagOut = 1;
3401     }
3402   } while (flagOut == 0 && nloopOut < loopN);
3403 
3404   if (flagOut == 0)
3405     ic = 1;
3406 
3407   //comment for unit test end
3408 }
3409 
3410 void LBT::radiationHQ(int parID, double qhat0ud, double v0[4], double P2[4],
3411                       double P3[4], double P4[4], double Pj0[4], int &ic,
3412                       double Tdiff, double HQenergy, double max_Ng,
3413                       double temp_med, double xLow, double xInt) {
3414 
3415   //    work in the rest frame of medium
3416   //    return the 4-momentum of final states in 1->3 radiation
3417   //    input: P2(4)-momentum of radiated gluon from 2->3
3418   //           P3(4)-momentum of daughter parton from 2->3
3419   //           Pj0(4)-inital momentum of jet before 2->3
3420   //           v0-local velocity of medium
3421   //    output:P2(4)-momentum of 1st radiated gluon
3422   //           P3(4)-momentum of daughter parton
3423   //           P4(4)-momentum of 2nd radiated gluon
3424   //           i=1: no radiation; 0:successful radiation
3425 
3426   double randomX, randomY;
3427   int count_sample, nloopOut, flagOut;
3428   double theta_gluon, kperp_gluon;
3429   double kpGluon[4];
3430   double HQmass =
3431       sqrt(P3[0] * P3[0] - P3[1] * P3[1] - P3[2] * P3[2] - P3[3] * P3[3]);
3432 
3433   if (abs(parID) != 4 && abs(parID) != 5) {
3434     HQmass = 0.0;
3435     P3[0] = sqrt(P3[1] * P3[1] + P3[2] * P3[2] + P3[3] * P3[3]);
3436   }
3437 
3438   double P50[4] = {0.0};
3439   double P2i[4] = {0.0};
3440   double P3i[4] = {0.0};
3441 
3442   ic = 0;
3443 
3444   for (int i = 0; i <= 3; i++) {
3445     P2i[i] = P2[i];
3446     P3i[i] = P3[i];
3447   }
3448 
3449   // transform to local comoving frame of the fluid
3450   trans(v0, P2);
3451   trans(v0, P3);
3452   trans(v0, Pj0);
3453 
3454   xInt = (P3[0] - HQmass) / HQenergy - xLow;
3455   if (xInt <= 0.0) {
3456     ic = 1;
3457     return;
3458   }
3459 
3460   double px0 = P2[1] + P3[1];
3461   double py0 = P2[2] + P3[2];
3462   double pz0 = P2[3] + P3[3];
3463 
3464   // rotate to the frame in which jet moves along z-axis
3465   rotate(px0, py0, pz0, P3, 1);
3466   rotate(px0, py0, pz0, P2, 1);
3467   rotate(px0, py0, pz0, Pj0, 1);
3468 
3469   for (int i = 0; i <= 3; i++) {
3470     P50[i] = P2[i] + P3[i];
3471   }
3472 
3473   nloopOut = 0;
3474   flagOut = 0;
3475   // sample the second gluon
3476 
3477   // comment do { for unit test
3478   do {
3479     do {
3480       randomX = xLow + xInt * ZeroOneDistribution(*GetMt19937Generator());
3481       randomY = ZeroOneDistribution(*GetMt19937Generator());
3482     } while (tau_f(randomX, randomY, HQenergy, HQmass) < 1.0 / pi / temp_med);
3483 
3484     count_sample = 0;
3485     while (max_Ng * ZeroOneDistribution(*GetMt19937Generator()) > dNg_over_dxdydt(parID, randomX, randomY,
3486                                                   HQenergy, HQmass, temp_med,
3487                                                   Tdiff)) {
3488       count_sample = count_sample + 1;
3489       if (count_sample > 1e+5) {
3490         //cout << "give up loop at point 1 ..." << endl;
3491         kpGluon[1] = 0.0;
3492         kpGluon[2] = 0.0;
3493         kpGluon[3] = 0.0;
3494         kpGluon[0] = 0.0;
3495         ic = 1;
3496         //          break;
3497         return;
3498       }
3499 
3500       do {
3501         randomX = xLow + xInt * ZeroOneDistribution(*GetMt19937Generator());
3502         randomY = ZeroOneDistribution(*GetMt19937Generator());
3503       } while (tau_f(randomX, randomY, HQenergy, HQmass) < 1.0 / pi / temp_med);
3504     }
3505 
3506     if (parID == 21 && randomX > 0.5)
3507       randomX = 1.0 - randomX;
3508     theta_gluon = 2.0 * pi * ZeroOneDistribution(*GetMt19937Generator());
3509     kperp_gluon = randomX * randomY * HQenergy;
3510     kpGluon[1] = kperp_gluon * cos(theta_gluon);
3511     kpGluon[2] = kperp_gluon * sin(theta_gluon);
3512     kpGluon[3] = randomX * HQenergy * sqrt(1.0 - randomY * randomY);
3513     kpGluon[0] = sqrt(kpGluon[1] * kpGluon[1] + kpGluon[2] * kpGluon[2] +
3514                       kpGluon[3] * kpGluon[3]);
3515 
3516     rotate(Pj0[1], Pj0[2], Pj0[3], kpGluon, -1);
3517 
3518     // comment for unit test
3519     if (kpGluon[0] >
3520         (P3[0] -
3521          HQmass)) { // which should be impossible due to reset of xInt above
3522       kpGluon[1] = 0.0;
3523       kpGluon[2] = 0.0;
3524       kpGluon[3] = 0.0;
3525       kpGluon[0] = 0.0;
3526       nloopOut++;
3527       continue;
3528     }
3529 
3530     // update mometum
3531 
3532     P4[1] = kpGluon[1];
3533     P4[2] = kpGluon[2];
3534     P4[3] = kpGluon[3];
3535     P4[0] = kpGluon[0];
3536 
3537     // comment for unit test begin
3538 
3539     // solve energy-momentum conservation
3540     // my notation in the note
3541     // p0 is re-constructed off-shell parton, correponding to P5
3542     // p1 is the final heavy quark from 2->3, corresponding to P3 (input). P3 (output) is the final heavy quark after 1->3.
3543     // k1 is the first gluon, corresponding to P2
3544     // k2 is the second gluon, corresponding to P4
3545     // assume k10 and p10 unchanged and modify their other components while p0 and k2 are fixed
3546     double sp0z = sqrt(P50[1] * P50[1] + P50[2] * P50[2] + P50[3] * P50[3]);
3547     double sk10 = P2[0];
3548     double sk1zOld = P2[3];
3549     double sk1z, sk1p, sk1x, sk1y, sktheta; // unknown
3550     double sp10 = P3[0];
3551     double sk20 = P4[0];
3552     double sk2x = P4[1];
3553     double sk2y = P4[2];
3554     double sk2z = P4[3];
3555     double sk2p = sqrt(sk2x * sk2x + sk2y * sk2y);
3556 
3557     double stheta12, cos12Min2;
3558     double sAA, aaa, bbb, ccc, abc, abc2;
3559     double sk1z1, sk1z2, sk1p1, sk1p2;
3560     int nloop1 = 0;
3561     int nloop2 = 0;
3562     int flagDone = 0;
3563     int yesA, yesB;
3564 
3565     sAA = sk10 * sk10 + sk2p * sk2p + sp0z * sp0z + sk2z * sk2z -
3566           2.0 * sp0z * sk2z + HQmass * HQmass - (sp10 - sk20) * (sp10 - sk20);
3567     cos12Min2 =
3568         (sAA * sAA / 4.0 / sk10 / sk10 - (sp0z - sk2z) * (sp0z - sk2z)) / sk2p /
3569         sk2p;
3570     //  cout << "cos^2 min: " << cos12Min2 << endl;
3571 
3572     if (cos12Min2 > 1.0) {
3573       nloopOut++;
3574       continue;
3575       //      cout << "cos^2 min is outside the kinematic range: " << cos12Min2 << endl;
3576       //      return;
3577     }
3578 
3579     do {
3580       stheta12 = 2.0 * pi * ZeroOneDistribution(*GetMt19937Generator()); // theta between k1 and k2
3581       aaa = 4.0 * ((sp0z - sk2z) * (sp0z - sk2z) +
3582                    sk2p * sk2p * cos(stheta12) * cos(stheta12));
3583       bbb = -4.0 * sAA * (sp0z - sk2z);
3584       ccc = sAA * sAA -
3585             4.0 * sk10 * sk10 * sk2p * sk2p * cos(stheta12) * cos(stheta12);
3586 
3587       abc2 = bbb * bbb - 4.0 * aaa * ccc;
3588 
3589       if (abc2 < 0.0) {
3590         nloop1++;
3591         continue;
3592       } else {
3593         nloop1 = 0;
3594       }
3595 
3596       abc = sqrt(abc2);
3597       sk1z1 = (-bbb + abc) / 2.0 / aaa;
3598       sk1z2 = (-bbb - abc) / 2.0 / aaa;
3599       sk1p1 = sqrt(sk10 * sk10 - sk1z1 * sk1z1);
3600       sk1p2 = sqrt(sk10 * sk10 - sk1z2 * sk1z2);
3601 
3602       // Since we have squared both sides of the original equation during solving, the solutions may not satisfy the original equation. Double check is needed.
3603 
3604       // require time-like of p1 and k10>k1z;
3605       if (2.0 * sk1p1 * sk2p * cos(stheta12) - 2.0 * (sp0z - sk2z) * sk1z1 +
3606                   sAA <
3607               0.000001 &&
3608           sk10 > abs(sk1z1) &&
3609           sp10 * sp10 - (sp0z - sk1z1) * (sp0z - sk1z1) -
3610                   (sk10 * sk10 - sk1z1 * sk1z1) >
3611               HQmass * HQmass) {
3612         yesA = 1;
3613       } else {
3614         yesA = 0;
3615       }
3616       if (2.0 * sk1p2 * sk2p * cos(stheta12) - 2.0 * (sp0z - sk2z) * sk1z2 +
3617                   sAA <
3618               0.000001 &&
3619           sk10 > abs(sk1z2) &&
3620           sp10 * sp10 - (sp0z - sk1z2) * (sp0z - sk1z2) -
3621                   (sk10 * sk10 - sk1z2 * sk1z2) >
3622               HQmass * HQmass) {
3623         yesB = 1;
3624       } else {
3625         yesB = 0;
3626       }
3627 
3628       // select appropriate solution
3629       if (yesA == 0 && yesB == 0) {
3630         //          cout << "Solutions fail ..." << endl;
3631         nloop2++;
3632         continue;
3633       } else if (yesA == 1 && yesB == 1) {
3634         //          cout << "Both solutions work!" << "  " << sE1 << "  " << sp1x << "  " << sp1y << "  " << sp1z << "  " << sE2 << "  " << sp2x << "  " << sp2y << "  " << sp2z << "  " << sk0 << "  " << skx << "  " << sky << "  " << skz << "  " << sqx << "  " << sqy << "  " << sq0A << "  " << sqzA << "  " << sq0B << "  " << sqzB << endl;
3635         if (abs(sk1z1 - sk1zOld) < abs(sk1z2 - sk1zOld)) {
3636           sk1z = sk1z1;
3637         } else {
3638           sk1z = sk1z2;
3639         }
3640       } else if (yesA == 1) {
3641         //          cout << "pass A ..." << "  " << sE1 << "  " << sp1x << "  " << sp1y << "  " << sp1z << "  " << sE2 << "  " << sp2x << "  " << sp2y << "  " << sp2z << "  " << sk0 << "  " << skx << "  " << sky << "  " << skz << "  " << sqx << "  " << sqy << "  " << sq0A << "  " << sqzA << "  " << sq0B << "  " << sqzB << endl;
3642         sk1z = sk1z1;
3643       } else {
3644         //          cout << "pass B ..." << "  " << sE1 << "  " << sp1x << "  " << sp1y << "  " << sp1z << "  " << sE2 << "  " << sp2x << "  " << sp2y << "  " << sp2z << "  " << sk0 << "  " << skx << "  " << sky << "  " << skz << "  " << sqx << "  " << sqy << "  " << sq0A << "  " << sqzA << "  " << sq0B << "  " << sqzB << endl;
3645         sk1z = sk1z2;
3646       }
3647 
3648       flagDone = 1;
3649 
3650       sk1p = sqrt(sk10 * sk10 - sk1z * sk1z);
3651       //           cout << "check solution: " << pow(sp10-sk20,2)-pow(sk1p,2)-pow(sk2p,2)-2.0*sk1p*sk2p*cos(stheta12)-pow(sp0z-sk1z-sk2z,2)-pow(HQmass,2) <<"  "<< aaa*sk1z*sk1z+bbb*sk1z+ccc <<"  "<<2.0*sk1p*sk2p*cos(stheta12)-2.0*(sp0z-sk2z)*sk1z+sAA<<"  "<<2.0*sk1p*sk2p*cos(stheta12)+2.0*(sp0z-sk2z)*sk1z-sAA << endl;
3652 
3653     } while (flagDone == 0 && nloop1 < loopN && nloop2 < loopN);
3654 
3655     if (flagDone == 0) {
3656       //       ic=1;  // no appropriate solution
3657       //       cout << "solution fails ...  " << nloop1 <<"  "<<nloop2<< endl;
3658       nloopOut++;
3659       continue;
3660     } else {
3661       sktheta = atan2(sk2y, sk2x);  // theta for k2
3662       sktheta = sktheta + stheta12; // theta for k1
3663       sk1p = sqrt(sk10 * sk10 - sk1z * sk1z);
3664       sk1x = sk1p * cos(sktheta);
3665       sk1y = sk1p * sin(sktheta);
3666 
3667       P2[0] = sk10;
3668       P2[1] = sk1x;
3669       P2[2] = sk1y;
3670       P2[3] = sk1z;
3671 
3672       P3[0] = sp10 - sk20;
3673       P3[1] = -sk1x - sk2x;
3674       P3[2] = -sk1y - sk2y;
3675       P3[3] = sp0z - sk1z - sk2z;
3676       // comment for unit test end
3677 
3678       // rotate back to the local coordinate
3679       rotate(px0, py0, pz0, P2, -1);
3680       rotate(px0, py0, pz0, P3, -1);
3681       rotate(px0, py0, pz0, P4, -1);
3682       rotate(px0, py0, pz0, Pj0, -1);
3683 
3684       // boost back to the global frame
3685       transback(v0, P2);
3686       transback(v0, P3);
3687       transback(v0, P4);
3688       transback(v0, Pj0);
3689 
3690       // comment for unit test begin
3691       // debug: check on-shell condition of P2, P3 and P4
3692       if (abs(P2[0] * P2[0] - P2[1] * P2[1] - P2[2] * P2[2] - P2[3] * P2[3]) >
3693               0.01 ||
3694           abs(P3[0] * P3[0] - P3[1] * P3[1] - P3[2] * P3[2] - P3[3] * P3[3] -
3695               HQmass * HQmass) > 0.01 ||
3696           abs(P4[0] * P4[0] - P4[1] * P4[1] - P4[2] * P4[2] - P4[3] * P4[3]) >
3697               0.01) {
3698         cout << "Wrong solution -- not on shell"
3699              << "  " << sk10 << "  " << sk1x << "  " << sk1y << "  " << sk1z
3700              << "  " << sk20 << "  " << sk2x << "  " << sk2y << "  " << sk2z
3701              << "  " << stheta12 << "  " << sp10 << "  " << sp10 - sk20 << "  "
3702              << -sk1x - sk2x << "  " << -sk1y - sk2y << "  " << sp0z << "  "
3703              << sp0z - sk1z - sk2z << "  " << HQmass << "  "
3704              << pow(sp10 - sk20, 2) - pow(sk1x + sk2x, 2) -
3705                     pow(sk1y + sk2y, 2) - pow(sp0z - sk1z - sk2z, 2) -
3706                     pow(HQmass, 2)
3707              << endl;
3708         cout << abs(P2[0] * P2[0] - P2[1] * P2[1] - P2[2] * P2[2] -
3709                     P2[3] * P2[3])
3710              << "  "
3711              << abs(P3[0] * P3[0] - P3[1] * P3[1] - P3[2] * P3[2] -
3712                     P3[3] * P3[3] - HQmass * HQmass)
3713              << "  "
3714              << abs(P4[0] * P4[0] - P4[1] * P4[1] - P4[2] * P4[2] -
3715                     P4[3] * P4[3])
3716              << endl;
3717       }
3718 
3719       //       cout<<"in radiation HQ -- final gluon1, gluon2 & HQ: "<<P2[0]<<"  "<<P4[0]<<"  "<<P3[0]<<endl;
3720 
3721       flagOut = 1;
3722 
3723       // SC: check energy-momentum conservation
3724       for (int i = 0; i <= 3; i++) {
3725         if (ic == 0 && abs(P2i[i] + P3i[i] - P2[i] - P3[i] - P4[i]) > 0.01) {
3726           cout << "Warning: Violation of E.M. conservation!  " << i << " "
3727                << abs(P2i[i] + P3i[i] - P2[i] - P3[i] - P4[i]) << endl;
3728         }
3729       }
3730 
3731       // SC: check on-shell condition
3732       double shell2 =
3733           abs(P2[0] * P2[0] - P2[1] * P2[1] - P2[2] * P2[2] - P2[3] * P2[3]);
3734       double shell3 = abs(P3[0] * P3[0] - P3[1] * P3[1] - P3[2] * P3[2] -
3735                           P3[3] * P3[3] - HQmass * HQmass);
3736       double shell4 =
3737           abs(P4[0] * P4[0] - P4[1] * P4[1] - P4[2] * P4[2] - P4[3] * P4[3]);
3738       if (ic == 0 && (shell2 > 0.01 || shell3 > 0.01 || shell4 > 0.01)) {
3739         cout << "Warning: Violation of on-shell: " << shell2 << "  " << shell3
3740              << "  " << shell4 << endl;
3741       }
3742     }
3743 
3744   } while (flagOut == 0 && nloopOut < loopN);
3745 
3746   if (flagOut == 0)
3747     ic = 1;
3748   // comment for unit test end
3749 }
3750 
3751 void LBT::rotate(double px, double py, double pz, double pr[4], int icc) {
3752   //     input:  (px,py,pz), (wx,wy,wz), argument (i)
3753   //     output: new (wx,wy,wz)
3754   //     if i=1, turn (wx,wy,wz) in the direction (px,py,pz)=>(0,0,E)
3755   //     if i=-1, turn (wx,wy,wz) in the direction (0,0,E)=>(px,py,pz)
3756 
3757   double wx, wy, wz, E, pt, w, cosa, sina, cosb, sinb;
3758   double wx1, wy1, wz1;
3759 
3760   wx = pr[1];
3761   wy = pr[2];
3762   wz = pr[3];
3763 
3764   E = sqrt(px * px + py * py + pz * pz);
3765   pt = sqrt(px * px + py * py);
3766 
3767   w = sqrt(wx * wx + wy * wy + wz * wz);
3768 
3769   if (pt == 0) {
3770     cosa = 1;
3771     sina = 0;
3772   } else {
3773     cosa = px / pt;
3774     sina = py / pt;
3775   }
3776 
3777   cosb = pz / E;
3778   sinb = pt / E;
3779 
3780   if (icc == 1) {
3781     wx1 = wx * cosb * cosa + wy * cosb * sina - wz * sinb;
3782     wy1 = -wx * sina + wy * cosa;
3783     wz1 = wx * sinb * cosa + wy * sinb * sina + wz * cosb;
3784   }
3785 
3786   else {
3787     wx1 = wx * cosa * cosb - wy * sina + wz * cosa * sinb;
3788     wy1 = wx * sina * cosb + wy * cosa + wz * sina * sinb;
3789     wz1 = -wx * sinb + wz * cosb;
3790   }
3791 
3792   wx = wx1;
3793   wy = wy1;
3794   wz = wz1;
3795 
3796   pr[1] = wx;
3797   pr[2] = wy;
3798   pr[3] = wz;
3799 
3800   //  pr[0]=sqrt(pr[1]*pr[1]+pr[2]*pr[2]+pr[3]*pr[3]);
3801 }
3802 
3803 int LBT::KPoisson(double alambda) {
3804   //....Generate numbers according to Poisson distribution
3805   //    P(lambda)=lambda**k exp(-lambda)/k!
3806   //    input: average number of radiated gluons lambda
3807   //    output: number of radiated gluons in this event
3808 
3809   double target, p;
3810 
3811   double KKPoisson = 0;
3812   target = exp(-alambda);
3813   p = ZeroOneDistribution(*GetMt19937Generator());
3814 
3815   while (p > target) {
3816     p = p * ZeroOneDistribution(*GetMt19937Generator());
3817     KKPoisson = KKPoisson + 1;
3818   }
3819   return KKPoisson;
3820 }
3821 
3822 double LBT::nHQgluon(int parID, double dtLRF, double &time_gluon,
3823                      double &temp_med, double &HQenergy, double &max_Ng) {
3824   // gluon radiation probability for heavy quark
3825 
3826   int time_num, temp_num, HQenergy_num;
3827   double delta_Ng;
3828   double rate_EGrid_low, rate_EGrid_high, max_EGrid_low, max_EGrid_high;
3829   double rate_T1E1, rate_T1E2, rate_T2E1, rate_T2E2, max_T1E1, max_T1E2,
3830       max_T2E1, max_T2E2;
3831 
3832   if (time_gluon > t_max) {
3833     cout << "accumulated time exceeds t_max" << endl;
3834     cout << time_gluon << "    " << temp_med << "    " << HQenergy << endl;
3835     time_gluon = t_max;
3836   }
3837 
3838   //if(temp_med>temp_max) {
3839   //  cout << "temperature exceeds temp_max -- extrapolation is used" << endl;
3840   //  cout << time_gluon << "    " << temp_med << "    " << HQenergy << endl;
3841   //  //     temp_med=temp_max;
3842   //}
3843 
3844   //if(HQenergy>HQener_max) {
3845   //  cout << "HQenergy exceeds HQener_max -- extrapolation is used" << endl;
3846   //  cout << time_gluon << "    " << temp_med << "    " << HQenergy << endl;
3847   //  //     HQenergy=HQener_max;
3848   //}
3849 
3850   if (temp_med < temp_min) {
3851     cout << "temperature drops below temp_min" << endl;
3852     cout << time_gluon << "    " << temp_med << "    " << HQenergy << endl;
3853     temp_med = temp_min;
3854   }
3855 
3856   if(time_gluon < t_max_1) {
3857       time_num = (int)(time_gluon / delta_tg_1 + 0.5) + 1;
3858   } else {
3859       time_num = (int)((time_gluon - t_max_1) / delta_tg_2 + 0.5) + t_gn_1 + 1;
3860   }
3861   //  temp_num=(int)((temp_med-temp_min)/delta_temp+0.5);
3862   //  HQenergy_num=(int)(HQenergy/delta_HQener+0.5); // use linear interpolation instead of finding nearest point for E and T dimensions
3863   temp_num = (int)((temp_med - temp_min) / delta_temp);
3864   HQenergy_num = (int)(HQenergy / delta_HQener); // normal interpolation
3865   if (HQenergy_num >= HQener_gn)
3866     HQenergy_num = HQener_gn - 1; // automatically become extrapolation
3867   if (temp_num >= temp_gn)
3868     temp_num = temp_gn - 1;
3869 
3870   if (parID == 21) {
3871     rate_T1E1 = dNg_over_dt_g[time_num][temp_num][HQenergy_num];
3872     rate_T1E2 = dNg_over_dt_g[time_num][temp_num][HQenergy_num + 1];
3873     rate_T2E1 = dNg_over_dt_g[time_num][temp_num + 1][HQenergy_num];
3874     rate_T2E2 = dNg_over_dt_g[time_num][temp_num + 1][HQenergy_num + 1];
3875     max_T1E1 = max_dNgfnc_g[time_num][temp_num][HQenergy_num];
3876     max_T1E2 = max_dNgfnc_g[time_num][temp_num][HQenergy_num + 1];
3877     max_T2E1 = max_dNgfnc_g[time_num][temp_num + 1][HQenergy_num];
3878     max_T2E2 = max_dNgfnc_g[time_num][temp_num + 1][HQenergy_num + 1];
3879   } else if (abs(parID) == 4 || abs(parID) == 5) {
3880     rate_T1E1 = dNg_over_dt_c[time_num][temp_num][HQenergy_num];
3881     rate_T1E2 = dNg_over_dt_c[time_num][temp_num][HQenergy_num + 1];
3882     rate_T2E1 = dNg_over_dt_c[time_num][temp_num + 1][HQenergy_num];
3883     rate_T2E2 = dNg_over_dt_c[time_num][temp_num + 1][HQenergy_num + 1];
3884     max_T1E1 = max_dNgfnc_c[time_num][temp_num][HQenergy_num];
3885     max_T1E2 = max_dNgfnc_c[time_num][temp_num][HQenergy_num + 1];
3886     max_T2E1 = max_dNgfnc_c[time_num][temp_num + 1][HQenergy_num];
3887     max_T2E2 = max_dNgfnc_c[time_num][temp_num + 1][HQenergy_num + 1];
3888   } else {
3889     rate_T1E1 = dNg_over_dt_q[time_num][temp_num][HQenergy_num];
3890     rate_T1E2 = dNg_over_dt_q[time_num][temp_num][HQenergy_num + 1];
3891     rate_T2E1 = dNg_over_dt_q[time_num][temp_num + 1][HQenergy_num];
3892     rate_T2E2 = dNg_over_dt_q[time_num][temp_num + 1][HQenergy_num + 1];
3893     max_T1E1 = max_dNgfnc_q[time_num][temp_num][HQenergy_num];
3894     max_T1E2 = max_dNgfnc_q[time_num][temp_num][HQenergy_num + 1];
3895     max_T2E1 = max_dNgfnc_q[time_num][temp_num + 1][HQenergy_num];
3896     max_T2E2 = max_dNgfnc_q[time_num][temp_num + 1][HQenergy_num + 1];
3897   }
3898 
3899   rate_EGrid_low = rate_T1E1 + (temp_med - temp_min - temp_num * delta_temp) /
3900                                    delta_temp * (rate_T2E1 - rate_T1E1);
3901   rate_EGrid_high = rate_T1E2 + (temp_med - temp_min - temp_num * delta_temp) /
3902                                     delta_temp * (rate_T2E2 - rate_T1E2);
3903   max_EGrid_low = max_T1E1 + (temp_med - temp_min - temp_num * delta_temp) /
3904                                  delta_temp * (max_T2E1 - max_T1E1);
3905   max_EGrid_high = max_T1E2 + (temp_med - temp_min - temp_num * delta_temp) /
3906                                   delta_temp * (max_T2E2 - max_T1E2);
3907 
3908   delta_Ng = rate_EGrid_low + (HQenergy - HQenergy_num * delta_HQener) /
3909                                   delta_HQener *
3910                                   (rate_EGrid_high - rate_EGrid_low);
3911   max_Ng = max_EGrid_low + (HQenergy - HQenergy_num * delta_HQener) /
3912                                delta_HQener * (max_EGrid_high - max_EGrid_low);
3913 
3914   delta_Ng = delta_Ng * 6.0 / D2piT * dtLRF;
3915   max_Ng = max_Ng * 6.0 / D2piT;
3916 
3917   //  if(delta_Ng>1) {
3918   //     cout << "Warning: Ng greater than 1   " << time_gluon << "  " << delta_Ng << endl;
3919   //  }
3920 
3921   return delta_Ng;
3922 
3923   //  write(6,*) temp_num,HQenergy_num,delta_Ng
3924 }
3925 
3926 //
3927 double LBT::Mqc2qc(double s, double t, double M) {
3928 
3929   double m2m = M * M;
3930   double u = 2.0 * m2m - s - t;
3931   double MM;
3932 
3933   MM = 64.0 / 9.0 * (pow((m2m - u), 2) + pow((s - m2m), 2) + 2.0 * m2m * t) /
3934        t / t;
3935 
3936   return (MM);
3937 }
3938 
3939 double LBT::Mgc2gc(double s, double t, double M) {
3940 
3941   double m2m = M * M;
3942   double u = 2.0 * m2m - s - t;
3943   double MM;
3944 
3945   MM = 32.0 * (s - m2m) * (m2m - u) / t / t;
3946   MM = MM + 64.0 / 9.0 * ((s - m2m) * (m2m - u) + 2.0 * m2m * (s + m2m)) /
3947                 pow((s - m2m), 2);
3948   MM = MM + 64.0 / 9.0 * ((s - m2m) * (m2m - u) + 2.0 * m2m * (u + m2m)) /
3949                 pow((u - m2m), 2);
3950   MM = MM + 16.0 / 9.0 * m2m * (4.0 * m2m - t) / ((s - m2m) * (m2m - u));
3951   MM = MM + 16.0 * ((s - m2m) * (m2m - u) + m2m * (s - u)) / (t * (s - m2m));
3952   MM = MM + 16.0 * ((s - m2m) * (m2m - u) - m2m * (s - u)) / (t * (u - m2m));
3953 
3954   return (MM);
3955 }
3956 
3957 ////////////////////////////////////////////////////////////////////////////////////////
3958 // functions for gluon radiation from heavy quark
3959 
3960 double LBT::alphasHQ(double kTFnc, double tempFnc) {
3961 
3962   //      double kTEff,error_para,resultFnc;
3963   //
3964   //      error_para=1.0;
3965   //
3966   //      if(kTFnc<pi*tempFnc*error_para) {
3967   //         kTEff=pi*tempFnc*error_para;
3968   //      } else {
3969   //         kTEff=kTFnc;
3970   //      }
3971   //
3972   //      resultFnc=4.0*pi/(11.0-2.0*nflavor(kTEff)/3.0)/2.0/log(kTEff/lambdas(kTEff));
3973   //
3974   //      return(resultFnc);
3975   return (alphas);
3976 }
3977 
3978 ////////////////////////////////////////////////////////////////////////////////////////
3979 
3980 double LBT::nflavor(double kTFnc) {
3981 
3982   double resultFnc;
3983   double cMass = 1.27;
3984   double bMass = 4.19;
3985 
3986   if (kTFnc < cMass) {
3987     resultFnc = 3.0;
3988   } else if (kTFnc < bMass) {
3989     resultFnc = 4.0;
3990   } else {
3991     resultFnc = 5.0;
3992   }
3993 
3994   return (resultFnc);
3995 }
3996 
3997 ////////////////////////////////////////////////////////////////////////////////////////
3998 
3999 double LBT::lambdas(double kTFnc) {
4000 
4001   double resultFnc;
4002   double cMass = 1.27;
4003   double bMass = 4.19;
4004 
4005   if (kTFnc < cMass) {
4006     resultFnc = 0.2;
4007   } else if (kTFnc < bMass) {
4008     resultFnc = 0.172508;
4009   } else {
4010     resultFnc = 0.130719;
4011   }
4012 
4013   return (resultFnc);
4014 }
4015 
4016 ////////////////////////////////////////////////////////////////////////////////////////
4017 
4018 double LBT::splittingP(int parID, double z0g) {
4019 
4020   double resultFnc;
4021 
4022   //      resultFnc = (2.0-2.0*z0g+z0g*z0g)*CF/z0g;
4023 
4024   if (parID == 21)
4025     resultFnc = 2.0 * pow(1.0 - z0g + pow(z0g, 2), 3) / z0g / (1.0 - z0g);
4026   else
4027     resultFnc = (1.0 - z0g) * (2.0 - 2.0 * z0g + z0g * z0g) / z0g;
4028 
4029   return (resultFnc);
4030 }
4031 
4032 ////////////////////////////////////////////////////////////////////////////////////////
4033 
4034 double LBT::tau_f(double x0g, double y0g, double HQenergy, double HQmass) {
4035 
4036   double resultFnc;
4037 
4038   resultFnc = 2.0 * HQenergy * x0g * (1.0 - x0g) /
4039               (pow((x0g * y0g * HQenergy), 2) + x0g * x0g * HQmass * HQmass);
4040 
4041   return (resultFnc);
4042 }
4043 
4044 ////////////////////////////////////////////////////////////////////////////////////////
4045 
4046 double LBT::dNg_over_dxdydt(int parID, double x0g, double y0g, double HQenergy,
4047                             double HQmass, double temp_med, double Tdiff) {
4048 
4049   double resultFnc, tauFnc, qhatFnc;
4050 
4051   qhatFnc = qhat_over_T3 *
4052             pow(temp_med,
4053                 3); // no longer have CF factor, taken out in splittingP too.
4054   tauFnc = tau_f(x0g, y0g, HQenergy, HQmass);
4055 
4056   resultFnc = 4.0 / pi * CA * alphasHQ(x0g * y0g * HQenergy, temp_med) *
4057               splittingP(parID, x0g) * qhatFnc * pow(y0g, 5) *
4058               pow(sin(Tdiff / 2.0 / tauFnc / sctr), 2) *
4059               pow((HQenergy * HQenergy /
4060                    (y0g * y0g * HQenergy * HQenergy + HQmass * HQmass)),
4061                   4) /
4062               x0g / x0g / HQenergy / HQenergy / sctr;
4063 
4064   return (resultFnc);
4065 }
4066 
4067 /////////////////////////////////////////////////////////////////////////////////////////
4068 
4069 void LBT::read_xyMC(int &numXY) {
4070 
4071   numXY = 0;
4072 
4073   ifstream fxyMC("../hydroProfile/mc_glauber.dat");
4074   if (!fxyMC.is_open()) {
4075     cout << "Erro openning date fxyMC!\n";
4076     exit(EXIT_FAILURE);
4077   }
4078 
4079   while (true) {
4080     if (fxyMC.eof()) {
4081       numXY--;
4082       break;
4083     }
4084     if (numXY >= maxMC)
4085       break;
4086     fxyMC >> initMCX[numXY] >> initMCY[numXY];
4087     numXY++;
4088   }
4089 
4090   cout << "Number of (x,y) points from MC-Glauber: " << numXY << endl;
4091 
4092   fxyMC.close();
4093 }
4094 
4095 void LBT::jetClean() {
4096 
4097   for (int i = 0; i < dimParList;
4098        i++) { // clear arrays before initialization for each event
4099     P[1][i] = 0.0;
4100     P[2][i] = 0.0;
4101     P[3][i] = 0.0;
4102     P[0][i] = 0.0;
4103     P[4][i] = 0.0;
4104     P[5][i] = 0.0;
4105     P[6][i] = 0.0;
4106     //
4107     P0[1][i] = 0.0;
4108     P0[2][i] = 0.0;
4109     P0[3][i] = 0.0;
4110     P0[0][i] = 0.0;
4111     P0[4][i] = 0.0;
4112     P0[5][i] = 0.0;
4113     P0[6][i] = 0.0;
4114     //
4115     PP[1][i] = 0.0;
4116     PP[2][i] = 0.0;
4117     PP[3][i] = 0.0;
4118     PP[0][i] = 0.0;
4119     //
4120     PP0[1][i] = 0.0;
4121     PP0[2][i] = 0.0;
4122     PP0[3][i] = 0.0;
4123     PP0[0][i] = 0.0;
4124     //
4125     V[1][i] = 0.0;
4126     V[2][i] = 0.0;
4127     V[3][i] = 0.0;
4128     //
4129     V0[1][i] = 0.0;
4130     V0[2][i] = 0.0;
4131     V0[3][i] = 0.0;
4132     //
4133     VV[1][i] = 0.0;
4134     VV[2][i] = 0.0;
4135     VV[3][i] = 0.0;
4136     //
4137     VV0[1][i] = 0.0;
4138     VV0[2][i] = 0.0;
4139     VV0[3][i] = 0.0;
4140     //
4141     CAT[i] = 0;
4142     CAT0[i] = 0;
4143     //
4144     radng[i] = 0.0;
4145     tirad[i] = 0.0;
4146     tiscatter[i] = 0.0;
4147     //
4148     Vfrozen[0][i] = 0.0;
4149     Vfrozen[1][i] = 0.0;
4150     Vfrozen[2][i] = 0.0;
4151     Vfrozen[3][i] = 0.0;
4152     Vfrozen0[0][i] = 0.0;
4153     Vfrozen0[1][i] = 0.0;
4154     Vfrozen0[2][i] = 0.0;
4155     Vfrozen0[3][i] = 0.0;
4156     Tfrozen[i] = 0.0;
4157     Tfrozen0[i] = 0.0;
4158     vcfrozen[0][i] = 0.0;
4159     vcfrozen[1][i] = 0.0;
4160     vcfrozen[2][i] = 0.0;
4161     vcfrozen[3][i] = 0.0;
4162     vcfrozen0[0][i] = 0.0;
4163     vcfrozen0[1][i] = 0.0;
4164     vcfrozen0[2][i] = 0.0;
4165     vcfrozen0[3][i] = 0.0;
4166 
4167     Tint_lrf[i] = 0.0; //for heavy quark
4168   }
4169 }
4170 
4171 void LBT::jetInitialize(int numXY) {
4172 
4173   double ipx, ipy, ipz, ip0, iwt;
4174   double pT_len, ipT, phi, rapidity;
4175 
4176   //...single jet parton input, only useful when fixMomentum=1
4177   double px0 = sqrt(ener * ener - amss * amss);
4178   double py0 = 0.0; // initial momemtum
4179   double pz0 = 0.0;
4180   double en0 = ener;
4181 
4182   for (int i = 1; i <= nj; i = i + 2) {
4183 
4184     if (fixPosition == 1) { // initialize position
4185       V[1][i] = 0.0;
4186       V[2][i] = 0.0;
4187       V[3][i] = 0.0;
4188     } else {
4189       int index_xy = (int)(ZeroOneDistribution(*GetMt19937Generator()) * numXY);
4190       if (index_xy >= numXY)
4191         index_xy = numXY - 1;
4192       V[1][i] = initMCX[index_xy];
4193       V[2][i] = initMCY[index_xy];
4194       V[3][i] = 0.0;
4195     }
4196     V[0][i] = -log(1.0 - ZeroOneDistribution(*GetMt19937Generator()));
4197 
4198     if (fixMomentum == 1) { // initialize momentum
4199       P[1][i] = px0;
4200       P[2][i] = py0;
4201       P[3][i] = pz0;
4202       P[0][i] = en0;
4203       P[4][i] = amss;
4204       P[5][i] = sqrt(px0 * px0 + py0 * py0);
4205       WT[i] = 1.0;
4206     } else {
4207       pT_len = ipTmax - ipTmin;
4208       ipT = ipTmin + ZeroOneDistribution(*GetMt19937Generator()) * pT_len;
4209       phi = ZeroOneDistribution(*GetMt19937Generator()) * 2.0 * pi;
4210       ipx = ipT * cos(phi);
4211       ipy = ipT * sin(phi);
4212       rapidity = 2.0 * eta_cut * ZeroOneDistribution(*GetMt19937Generator()) - eta_cut;
4213       ipz = sqrt(ipT * ipT + amss * amss) * sinh(rapidity);
4214       ip0 = sqrt(ipT * ipT + ipz * ipz + amss * amss);
4215       P[1][i] = ipx;
4216       P[2][i] = ipy;
4217       P[3][i] = ipz;
4218       P[0][i] = ip0;
4219       P[4][i] = amss;
4220       P[5][i] = ipT;
4221       WT[i] = pT_len;
4222     }
4223 
4224     KATT1[i] = Kjet;
4225 
4226     // let jet partons stream freely for tau_0 in x-y plane;
4227     V[1][i] = V[1][i] + P[1][i] / P[0][i] * tau0;
4228     V[2][i] = V[2][i] + P[2][i] / P[0][i] * tau0;
4229 
4230     for (int j = 0; j <= 3; j++)
4231       Prad[j][i] = P[j][i];
4232 
4233     Vfrozen[0][i] = tau0;
4234     Vfrozen[1][i] = V[1][i];
4235     Vfrozen[2][i] = V[2][i];
4236     Vfrozen[3][i] = V[3][i];
4237     Tfrozen[i] = 0.0;
4238     vcfrozen[1][i] = 0.0;
4239     vcfrozen[2][i] = 0.0;
4240     vcfrozen[3][i] = 0.0;
4241 
4242     // now its anti-particle (if consider pair production)
4243     if (KATT1[i] == 21)
4244       KATT1[i + 1] = 21;
4245     else
4246       KATT1[i + 1] = -KATT1[i];
4247     P[1][i + 1] = -P[1][i];
4248     P[2][i + 1] = -P[2][i];
4249     P[3][i + 1] = -P[3][i];
4250     P[0][i + 1] = P[0][i];
4251     P[4][i + 1] = P[4][i];
4252     P[5][i + 1] = P[5][i];
4253     WT[i + 1] = WT[i];
4254     V[1][i + 1] = V[1][i];
4255     V[2][i + 1] = V[2][i];
4256     V[3][i + 1] = V[3][i];
4257     V[0][i + 1] = V[0][i];
4258     for (int j = 0; j <= 3; j++)
4259       Prad[j][i + 1] = P[j][i];
4260     Vfrozen[0][i + 1] = Vfrozen[0][i];
4261     Vfrozen[1][i + 1] = Vfrozen[1][i];
4262     Vfrozen[2][i + 1] = Vfrozen[2][i];
4263     Vfrozen[3][i + 1] = Vfrozen[3][i];
4264     Tfrozen[i + 1] = Tfrozen[i];
4265     vcfrozen[1][i + 1] = vcfrozen[1][i];
4266     vcfrozen[2][i + 1] = vcfrozen[2][i];
4267     vcfrozen[3][i + 1] = vcfrozen[3][i];
4268   }
4269 }
4270 
4271 void LBT::setJetX(int numXY) {
4272 
4273   double setX, setY, setZ;
4274 
4275   if (fixPosition == 1) {
4276     setX = 0.0;
4277     setY = 0.0;
4278     setZ = 0.0;
4279   } else {
4280     int index_xy = (int)(ZeroOneDistribution(*GetMt19937Generator()) * numXY);
4281     if (index_xy >= numXY)
4282       index_xy = numXY - 1;
4283     setX = initMCX[index_xy];
4284     setY = initMCY[index_xy];
4285     setZ = 0.0;
4286   }
4287 
4288   for (int i = 1; i <= nj; i++) {
4289 
4290     V[1][i] = setX;
4291     V[2][i] = setY;
4292     V[3][i] = setZ;
4293     V[0][i] = -log(1.0 - ZeroOneDistribution(*GetMt19937Generator()));
4294 
4295     V[1][i] = V[1][i] + P[1][i] / P[0][i] * tau0;
4296     V[2][i] = V[2][i] + P[2][i] / P[0][i] * tau0;
4297 
4298     Vfrozen[0][i] = tau0;
4299     Vfrozen[1][i] = V[1][i];
4300     Vfrozen[2][i] = V[2][i];
4301     Vfrozen[3][i] = V[3][i];
4302     Tfrozen[i] = 0.0;
4303     vcfrozen[1][i] = 0.0;
4304     vcfrozen[2][i] = 0.0;
4305     vcfrozen[3][i] = 0.0;
4306   }
4307 }
4308 
4309 ////////////////////////////////////////////////////////////////////////////////////////
4310 
4311 void LBT::setParameter(string fileName) {
4312 
4313   const char *cstring = fileName.c_str();
4314 
4315   ifstream input(cstring);
4316 
4317   if (!input) {
4318     cout << "Parameter file is missing!" << endl;
4319     exit(EXIT_FAILURE);
4320   }
4321 
4322   string line;
4323 
4324   while (getline(input, line)) {
4325     char str[1024];
4326     strncpy(str, line.c_str(), sizeof(str)-1);
4327     //          cout << str << endl;
4328 
4329     char *divideChar[5];
4330     divideChar[0] = strtok(str, " ");
4331 
4332     int nChar = 0;
4333     while (divideChar[nChar] != NULL) {
4334       if (nChar == 3)
4335         break;
4336       //              cout << divideChar[nChar] << endl;
4337       nChar++;
4338       divideChar[nChar] = strtok(NULL, " ");
4339     }
4340 
4341     if (nChar == 3) {
4342       //              cout << divideChar[0] << endl;
4343       //              cout << divideChar[1] << endl;
4344       //              cout << divideChar[2] << endl;
4345       string firstStr(divideChar[0]);
4346       string secondStr(divideChar[1]);
4347       string thirdStr(divideChar[2]);
4348       //              cout << firstStr << "     " << secondStr << endl;
4349       if (secondStr == "=") {
4350         if (firstStr == "flag:fixMomentum")
4351           fixMomentum = atoi(divideChar[2]);
4352         else if (firstStr == "flag:fixPosition")
4353           fixPosition = atoi(divideChar[2]);
4354         else if (firstStr == "flag:initHardFlag")
4355           initHardFlag = atoi(divideChar[2]);
4356         else if (firstStr == "flag:flagJetX")
4357           flagJetX = atoi(divideChar[2]);
4358         else if (firstStr == "flag:vacORmed")
4359           vacORmed = atoi(divideChar[2]);
4360         else if (firstStr == "flag:bulkType")
4361           bulkFlag = atoi(divideChar[2]);
4362         else if (firstStr == "flag:Kprimary")
4363           Kprimary = atoi(divideChar[2]);
4364         else if (firstStr == "flag:KINT0")
4365           KINT0 = atoi(divideChar[2]);
4366         else if (firstStr == "para:nEvent")
4367           ncall = atoi(divideChar[2]);
4368         else if (firstStr == "para:nParton")
4369           nj = atoi(divideChar[2]);
4370         else if (firstStr == "para:parID")
4371           Kjet = atoi(divideChar[2]);
4372         else if (firstStr == "para:tau0")
4373           tau0 = atof(divideChar[2]);
4374         else if (firstStr == "para:tauend")
4375           tauend = atof(divideChar[2]);
4376         else if (firstStr == "para:dtau")
4377           dtau = atof(divideChar[2]);
4378         else if (firstStr == "para:temperature")
4379           temp0 = atof(divideChar[2]);
4380         else if (firstStr == "para:alpha_s")
4381           fixAlphas = atof(divideChar[2]);
4382         else if (firstStr == "para:hydro_Tc")
4383           hydro_Tc = atof(divideChar[2]);
4384         else if (firstStr == "para:pT_min")
4385           ipTmin = atof(divideChar[2]);
4386         else if (firstStr == "para:pT_max")
4387           ipTmax = atof(divideChar[2]);
4388         else if (firstStr == "para:eta_cut")
4389           eta_cut = atof(divideChar[2]);
4390         else if (firstStr == "para:Ecut")
4391           Ecut = atof(divideChar[2]);
4392         else if (firstStr == "para:ener")
4393           ener = atof(divideChar[2]);
4394         else if (firstStr == "para:mass")
4395           amss = atof(divideChar[2]);
4396         else if (firstStr == "para:KPamp")
4397           KPamp = atof(divideChar[2]);
4398         else if (firstStr == "para:KPsig")
4399           KPsig = atof(divideChar[2]);
4400         else if (firstStr == "para:KTamp")
4401           KTamp = atof(divideChar[2]);
4402         else if (firstStr == "para:KTsig")
4403           KTsig = atof(divideChar[2]);
4404       }
4405     }
4406   }
4407 
4408   // calculated derived quantities
4409   if (vacORmed == 0)
4410     fixPosition = 1; // position is not important in vacuumm
4411 
4412   KTsig = KTsig * hydro_Tc;
4413   preKT = fixAlphas / 0.3;
4414 
4415   //// write out parameters
4416   //cout << endl;
4417   //cout << "############################################################" << endl;
4418   //cout << "Parameters for this LBT run" << endl;
4419   ////       cout << "flag:initHardFlag: " << initHardFlag << endl;
4420   ////       cout << "flag:flagJetX: " << flagJetX << endl;
4421   ////       cout << "flag:fixMomentum: " << fixMomentum << endl;
4422   ////       cout << "flag:fixPosition: " << fixPosition << endl;
4423   ////       cout << "flag:vacORmed: " << vacORmed << endl;
4424   //cout << "flag:bulkType: " << bulkFlag << endl;
4425   //cout << "flag:Kprimary: " << Kprimary << endl;
4426   //cout << "flag:KINT0: " << KINT0 << endl;
4427   //cout << "para:nEvent: " << ncall << endl;
4428   ////       cout << "para:nParton: " << nj << endl;
4429   ////       cout << "para:parID: " << Kjet << endl;
4430   ////       cout << "para:tau0: " << tau0 << endl;
4431   ////       cout << "para:tauend: " << tauend << endl;
4432   //cout << "para:dtau: " << dtau << endl;
4433   //cout << "para:temperature: " << temp0 << endl;
4434   //cout << "para:alpha_s: " << fixAlphas << endl;
4435   //cout << "para:hydro_Tc: " << hydro_Tc << endl;
4436   ////       cout << "para:pT_min: " << ipTmin << endl;
4437   ////       cout << "para:pT_max: " << ipTmax << endl;
4438   ////       cout << "para:eta_cut: " << eta_cut << endl;
4439   ////       cout << "para:ener: " << ener << endl;
4440   ////       cout << "para:mass: " << amss << endl;
4441   //cout << "para:KPamp: " << KPamp << endl;
4442   //cout << "para:KPsig: " << KPsig << endl;
4443   //cout << "para:KTamp: " << KTamp << endl;
4444   //cout << "para:KTsig: " << KTsig << endl;
4445   //cout << "############################################################" << endl;
4446   //cout << endl;
4447 
4448   //       exit(EXIT_SUCCESS);
4449 }
4450 
4451 ////////////////////////////////////////////////////////////////////////////////////////
4452 
4453 int LBT::checkParameter(int nArg) {
4454 
4455   int ctErr = 0;
4456 
4457   // better to make sure lightOut,heavyOut,fixPosition,fixMomentum,etc are either 0 or 1
4458   // not necessary at this moment
4459 
4460   if (initHardFlag == 1) { // initialize within LBT, no input particle list
4461     if (lightOut == 1 && heavyOut == 1) { // need both heavy and light output
4462       if (nArg != 5) {
4463         ctErr++;
4464         cout << "Wrong # of arguments for command line input" << endl;
4465         cout << "./LBT parameter_file HQ_output light_positive_output "
4466                 "light_negative_output"
4467              << endl;
4468       }
4469     } else if (heavyOut == 1) { // only heavy output
4470       if (nArg != 3) {
4471         ctErr++;
4472         cout << "Wrong # of arguments for command line input" << endl;
4473         cout << "./LBT parameter_file  HQ_output" << endl;
4474       }
4475     } else if (lightOut == 1) { // only light output
4476       if (nArg != 4) {
4477         ctErr++;
4478         cout << "Wrong # of arguments for command line input" << endl;
4479         cout << "./LBT parameter_file light_positive_output "
4480                 "light_negative_output"
4481              << endl;
4482       }
4483     } else { // no output
4484       if (nArg != 2) {
4485         ctErr++;
4486         cout << "Wrong # of arguments for command line input" << endl;
4487         cout << "./LBT parameter_file" << endl;
4488         cout << "No output specified, both heavyOut and lightOut are set to 0"
4489              << endl;
4490       }
4491     }
4492   } else if (initHardFlag == 2) {         // need input particle list
4493     if (lightOut == 1 && heavyOut == 1) { // need both heavy and light output
4494       if (nArg != 6) {
4495         ctErr++;
4496         cout << "Wrong # of arguments for command line input" << endl;
4497         cout << "./LBT parameter_file input_parton_list HQ_output "
4498                 "light_positive_output light_negative_output"
4499              << endl;
4500       }
4501     } else if (heavyOut == 1) { // only heavy output
4502       if (nArg != 4) {
4503         ctErr++;
4504         cout << "Wrong # of arguments for command line input" << endl;
4505         cout << "./LBT parameter_file input_parton_list HQ_output" << endl;
4506       }
4507     } else if (lightOut == 1) { // only light output
4508       if (nArg != 5) {
4509         ctErr++;
4510         cout << "Wrong # of arguments for command line input" << endl;
4511         cout << "./LBT parameter_file input_parton_list light_positive_output "
4512                 "light_negative_output"
4513              << endl;
4514       }
4515     } else { // no output
4516       if (nArg != 3) {
4517         ctErr++;
4518         cout << "Wrong # of arguments for command line input" << endl;
4519         cout << "./LBT parameter_file input_parton_list" << endl;
4520         cout << "No output specified, both heavyOut and lightOut are set to 0"
4521              << endl;
4522       }
4523     }
4524   } else {
4525     cout << "Wrong input for initHardFlag!" << endl;
4526     ctErr++;
4527   }
4528 
4529   return (ctErr);
4530 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!