Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LBT.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file LBT.cc
1 /*******************************************************************************
2  * Copyright (c) The JETSCAPE Collaboration, 2018
3  *
4  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
5  *
6  * For the list of contributors see AUTHORS.
7  *
8  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
9  *
10  * or via email to bugs.jetscape@gmail.com
11  *
12  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
13  * See COPYING for details.
14  ******************************************************************************/
15 
16 #include "LBT.h"
17 #include "JetScapeLogger.h"
18 #include "JetScapeXML.h"
19 
20 #include "tinyxml2.h"
21 #include <string>
22 #include <sstream>
23 #include <fstream>
24 #include <iostream>
25 #include <iomanip>
26 
27 #include "FluidDynamics.h"
28 #include "LBTMutex.h"
29 #define MAGENTA "\033[35m"
30 
31 using namespace Jetscape;
32 using namespace std;
33 
34 // Register the module with the base class
35 RegisterJetScapeModule<LBT> LBT::reg("Lbt");
36 
37 // initialize static members
38 bool LBT::flag_init = 0;
39 double LBT::Rg[60][20] = {
40  {0.0}}; //total gluon scattering rate as functions of initial energy and temperature
41 double LBT::Rg1[60][20] = {{0.0}}; //gg-gg CT1
42 double LBT::Rg2[60][20] = {{0.0}}; //gg-qqbar CT2
43 double LBT::Rg3[60][20] = {{0.0}}; //gq-qg CT3
44 double LBT::Rq[60][20] = {
45  {0.0}}; //total gluon scattering rate as functions of initial energy and temperature
46 double LBT::Rq3[60][20] = {{0.0}}; //qg-qg CT13
47 double LBT::Rq4[60][20] = {{0.0}}; //qiqj-qiqj CT4
48 double LBT::Rq5[60][20] = {{0.0}}; //qiqi-qiqi CT5
49 double LBT::Rq6[60][20] = {{0.0}}; //qiqibar-qjqjbar CT6
50 double LBT::Rq7[60][20] = {{0.0}}; //qiqibar-qiqibar CT7
51 double LBT::Rq8[60][20] = {{0.0}}; //qqbar-gg CT8
52 double LBT::qhatLQ[60][20] = {{0.0}};
53 double LBT::qhatG[60][20] = {{0.0}};
54 
55 double LBT::RHQ[60][20] = {{0.0}}; //total scattering rate for heavy quark
56 double LBT::RHQ11[60][20] = {{0.0}}; //Qq->Qq
57 double LBT::RHQ12[60][20] = {{0.0}}; //Qg->Qg
58 double LBT::qhatHQ[60][20] = {{0.0}}; //qhat of heavy quark
59 
60 double LBT::dNg_over_dt_c[t_gn + 2][temp_gn + 1][HQener_gn + 1] = {{{0.0}}};
61 double LBT::dNg_over_dt_q[t_gn + 2][temp_gn + 1][HQener_gn + 1] = {{{0.0}}};
62 double LBT::dNg_over_dt_g[t_gn + 2][temp_gn + 1][HQener_gn + 1] = {{{0.0}}};
63 double LBT::max_dNgfnc_c[t_gn + 2][temp_gn + 1][HQener_gn + 1] = {{{0.0}}};
64 double LBT::max_dNgfnc_q[t_gn + 2][temp_gn + 1][HQener_gn + 1] = {{{0.0}}};
65 double LBT::max_dNgfnc_g[t_gn + 2][temp_gn + 1][HQener_gn + 1] = {{{0.0}}};
66 
67 double LBT::initMCX[maxMC] = {0.0};
68 double LBT::initMCY[maxMC] = {0.0};
69 double LBT::distFncB[N_T][N_p1][N_e2] = {{{0.0}}};
70 double LBT::distFncF[N_T][N_p1][N_e2] = {{{0.0}}};
71 double LBT::distMaxB[N_T][N_p1][N_e2] = {{{0.0}}};
72 double LBT::distMaxF[N_T][N_p1][N_e2] = {{{0.0}}};
73 double LBT::distFncBM[N_T][N_p1] = {{0.0}};
74 double LBT::distFncFM[N_T][N_p1] = {{0.0}};
75 
76 LBT::LBT() {
77  SetId("LBT");
78  VERBOSE(8);
79 
80  // create and set Martini Mutex
81  auto lbt_mutex = make_shared<LBTMutex>();
82  SetMutex(lbt_mutex);
83 }
84 
85 LBT::~LBT() { VERBOSE(8); }
86 
87 void LBT::Init() {
88  JSINFO << "Initialize LBT ...";
89 
90  //...Below is added by Shanshan
91  //...read parameters from LBT.input first, but can be changed in JETSCAPE xml file if any conflict exists
92 
93  setParameter(
94  "LBT-tables/LBT.input"); // re-set parameters from input parameter file
95 
96  // if(checkParameter(argc)==0) { // check whether the input parameters are all correct
97  // cout << "Parameter check passed" << endl;
98  // } else {
99  // cout << "Parameter check failed" << endl;
100  // exit(EXIT_FAILURE);
101  // }
102 
103  std::string s = GetXMLElementText({"Eloss", "Lbt", "name"});
104  JSDEBUG << s << " to be initialized ...";
105 
106  //double m_qhat=-99.99;
107  //lbt->FirstChildElement("qhat")->QueryDoubleText(&m_qhat);
108  //SetQhat(m_qhat);
109  //
110  //JSDEBUG << s << " with qhat = "<<GetQhat();
111 
112  int in_vac = GetXMLElementInt({"Eloss", "Lbt", "in_vac"});
113  if (in_vac == 1) {
114  vacORmed = 0;
115  fixPosition = 1;
116  } else {
117  vacORmed = 1;
118  }
119 
120  Kprimary = GetXMLElementInt({"Eloss", "Lbt", "only_leading"});
121  run_alphas = GetXMLElementInt({"Eloss", "Lbt", "run_alphas"});
122  Q00 = GetXMLElementDouble({"Eloss", "Lbt", "Q0"});
123  fixAlphas = GetXMLElementDouble({"Eloss", "Lbt", "alphas"});
124  hydro_Tc = GetXMLElementDouble({"Eloss", "Lbt", "hydro_Tc"});
125  tStart = GetXMLElementDouble({"Eloss", "tStart"});
126  JSINFO << MAGENTA << "LBT parameters -- in_med: " << vacORmed
127  << " Q0: " << Q00 << " only_leading: " << Kprimary
128  << " alpha_s: " << fixAlphas << " hydro_Tc: " << hydro_Tc<<", tStart="<<tStart;
129 
130  if (!flag_init) {
131  read_tables(); // initialize various tables
132  flag_init = true;
133  }
134 
135  //...define derived quantities
136  temp00 = temp0;
137  dt = dtau;
138  timend = tauend;
139  time0 = tau0;
140  //...alphas
141  alphas = alphas0(Kalphas, temp0);
142  //...Debye Mass square
143  qhat0 = DebyeMass2(Kqhat0, alphas, temp0);
144 
145  ZeroOneDistribution = uniform_real_distribution<double>{0.0, 1.0};
146 }
147 
148 void LBT::WriteTask(weak_ptr<JetScapeWriter> w) {
149  VERBOSE(8);
150  auto f = w.lock();
151  if (!f)
152  return;
153  f->WriteComment("ElossModule Parton List: " + GetId());
154  f->WriteComment("Energy loss to be implemented accordingly ...");
155 }
156 
157 void LBT::DoEnergyLoss(double deltaT, double time, double Q2,
158  vector<Parton> &pIn, vector<Parton> &pOut) {
159 
160  double z = 0.5;
161 
162  // if (Q2>5)
163  // {
164  VERBOSESHOWER(8) << MAGENTA << "SentInPartons Signal received : " << deltaT
165  << " " << Q2 << " " << &pIn;
166 
167  // hydro information should be moved into the particle list loop
168 
170  //std::unique_ptr<FluidCellInfo> check_fluid_info_ptr;
171  //GetHydroCellSignal(1, 1.0, 1.0, 0.0, check_fluid_info_ptr);
172  //VERBOSE(8)<< MAGENTA<<"Temperature from Brick (Signal) = "<<check_fluid_info_ptr->temperature;
173 
174  double rNum;
175  int par_status;
176 
177  //DEBUG:
178  //cout<<" ---> "<<pIn.size()<<endl;
179  for (int i = 0; i < pIn.size(); i++) {
180 
181  // Reject photons
182 
183  if (pIn[i].pid() == photonid) {
184  if(pIn[i].pstat() != 22) {
185  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.
186  pOut.push_back(pIn[i]);
187  }
188  return;
189  }
190 
191  // pass particle infomation to LBT array (only pass one particle each time)
192 
193  jetClean();
194 
195  int j = 1;
196 
198  // KATT1[j] = Kjet;
199  // P[1][j] = ener;
200  // P[2][j] = 0.0;
201  // P[3][j] = 0.0;
202 
203  // 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;
204 
205  par_status = pIn[i].pstat();
206 
207  if (par_status != -1) { // a positive (active) particle
208 
209  KATT1[j] = pIn[i].pid();
210  P[1][j] = pIn[i].p(1);
211  P[2][j] = pIn[i].p(2);
212  P[3][j] = pIn[i].p(3);
213  P[4][j] = amss;
214  if (std::abs(pIn[i].pid()) == 4 || std::abs(pIn[i].pid()) == 5)
215  P[4][j] = pIn[i].restmass();
216  P[0][j] = sqrt(P[1][j] * P[1][j] + P[2][j] * P[2][j] + P[3][j] * P[3][j] +
217  P[4][j] * P[4][j]);
218  P[5][j] = sqrt(P[1][j] * P[1][j] + P[2][j] * P[2][j]);
219  P[6][j] = pIn[i].e() * pIn[i].e() - P[0][j] * P[0][j]; // virtuality^2
220  if (P[6][j] < 0.0)
221  P[6][j] = 0.0;
222  //if(std::isnan(P[6][j]) || std::isinf(P[6][j]))
223  //{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];}
224 
225  WT[j] = 1.0;
226 
227  Vfrozen[1][j] = pIn[i].x_in().x();
228  Vfrozen[2][j] = pIn[i].x_in().y();
229  Vfrozen[3][j] = pIn[i].x_in().z();
230  Vfrozen[0][j] = pIn[i].x_in().t();
231  V[1][j] = Vfrozen[1][j];
232  V[2][j] = Vfrozen[2][j];
233  V[3][j] = Vfrozen[3][j];
234  V[0][j] = -log(1.0 - ZeroOneDistribution(*GetMt19937Generator()));
235 
236  for (int k = 0; k <= 3; k++)
237  Prad[k][j] = P[k][j];
238 
239  Tfrozen[j] = 0.0;
240  vcfrozen[1][j] = 0.0;
241  vcfrozen[2][j] = 0.0;
242  vcfrozen[3][j] = 0.0;
243 
244  Tint_lrf[j] =
245  0.0; // time since last splitting, reset below if there is information from JETSCAPE
246  if (pIn[i].has_user_info<LBTUserInfo>()) {
247  Tint_lrf[j] = pIn[i].user_info<LBTUserInfo>().lrf_T_tot();
248  } else {
249  pIn[i].set_user_info(new LBTUserInfo(0.0));
250  }
251 
252  if (par_status == 1)
253  CAT[j] = 2;
254  else
255  CAT[j] = 0;
256 
257  } else { // negative particle, or particle no longer active, will streamly freely in LBT
258 
259  KATT10[j] = pIn[i].pid();
260  P0[1][j] = pIn[i].p(1);
261  P0[2][j] = pIn[i].p(2);
262  P0[3][j] = pIn[i].p(3);
263  P0[4][j] = amss;
264  P0[0][j] = sqrt(P0[1][j] * P0[1][j] + P0[2][j] * P0[2][j] +
265  P0[3][j] * P0[3][j] + P0[4][j] * P0[4][j]);
266  P0[5][j] = sqrt(P0[1][j] * P0[1][j] + P0[2][j] * P0[2][j]);
267  P0[6][j] = pIn[i].e() * pIn[i].e() - P0[0][j] * P0[0][j]; // virtuality^2
268  if (P0[6][j] < 0.0)
269  P0[6][j] = 0.0;
270  WT0[j] = 1.0;
271 
272  Vfrozen0[1][j] = pIn[i].x_in().x();
273  Vfrozen0[2][j] = pIn[i].x_in().y();
274  Vfrozen0[3][j] = pIn[i].x_in().z();
275  Vfrozen0[0][j] = pIn[i].x_in().t();
276  V0[1][j] = Vfrozen0[1][j];
277  V0[2][j] = Vfrozen0[2][j];
278  V0[3][j] = Vfrozen0[3][j];
279  V0[0][j] = -log(1.0 - ZeroOneDistribution(*GetMt19937Generator()));
280 
281  // for(int k=0;k<=3;k++) Prad[k][j] = P[k][j];
282 
283  Tfrozen0[j] = 0.0;
284  vcfrozen0[1][j] = 0.0;
285  vcfrozen0[2][j] = 0.0;
286  vcfrozen0[3][j] = 0.0;
287 
288  Tint_lrf[j] =
289  0.0; // time since last splitting, reset below if there is information from JETSCAPE
290  if (pIn[i].has_user_info<LBTUserInfo>()) {
291  Tint_lrf[j] = pIn[i].user_info<LBTUserInfo>().lrf_T_tot();
292  } else {
293  pIn[i].set_user_info(new LBTUserInfo(0.0));
294  }
295 
296  } // end par_status check
297 
303 
304  nj = 1;
305  np = 1;
306 
307  dtau = deltaT;
308  dt = dtau;
309 
310  flag_update = 0;
311  flag_update0 = 0;
312 
313  // do LBT evolution for this particle
314  int nnn = 1; // dummy variable, not important
315  // double systemTime=Vfrozen[0][j]+dt; // how to pass the time of the framework?
316  double systemTime = time;
317 
319  // if(systemTime-dtau<0.00001) {
321  // dEel=0.0;
322  // nGluon=0.0;
323  // eGluon=0.0;
324  // ctEvt++;
325  // }
326 
327  //if(alphas>epsilon && vacORmed!=0) {
328  if (vacORmed != 0) { // for JETSCAPE, won't change parton if it's in vacuum
329  flagScatter = 0;
330  LBT0(
331  nnn,
332  systemTime); // most important function of LBT: update particle list for one time step
333  }
334 
336  // dEel=dEel+(ener-P[0][1]);
338  // for(int ik=1; ik<=40; ++ik) {
339  // if(fabs(systemTime-ik*0.1) < 0.000001) {
340  // de1[ik] = de1[ik] + dEel/ncall;
341  // de2[ik] = de2[ik] + eGluon/ncall;
342  // }
343  // }
344  // // write out energy loss infomation for LBT unit test
345  // if(ctEvt==ncall && fabs(systemTime-tauend)<0.000001) {
346  // ofstream dat_de("de.dat");
347  // for(int ik=1; ik<=40; ++ik)
348  // dat_de << ik*0.1 << " " << de1[ik] << " " << de2[ik] << endl;
349  // dat_de.close();
350  // }
351 
352  double velocity_jet[4];
353  velocity_jet[0] = 1.0;
354  velocity_jet[1] = pIn[i].jet_v().x();
355  velocity_jet[2] = pIn[i].jet_v().y();
356  velocity_jet[3] = pIn[i].jet_v().z();
357 
358  if (par_status !=
359  -1) { // for particle list initiating with positive particles
360 
361  if (flag_update == 0)
362  continue; // no update on this particle from LBT0 function
363 
364  TakeResponsibilityFor(
365  pIn[i]); // Generate error if another module already has responsibility.
366 
368 
369  // pass particle list back to JETSCAPE
370  // how to handle negative particle in the framework?
371  // only pass positive at this moment
372 
373  // may only need to push back particle if flagScatter=1 later
374  // positive particle stat = 0
375  for (int j = 1; j <= np; j++) {
376  // definition Parton (int label, int id, int stat, double p[4], double x[4]);
377  if (P[0][j] < cutOut)
378  continue;
379  int out_stat = 0; // jet parton
380  if (CAT[j] == 2)
381  out_stat = 1; // recoil parton
382  double tempP[4], tempX[4];
383  tempP[0] = sqrt(P[0][j] * P[0][j] + P[6][j]);
384  tempP[1] = P[1][j];
385  tempP[2] = P[2][j];
386  tempP[3] = P[3][j];
387  tempX[0] = Vfrozen[0][j];
388  tempX[1] = Vfrozen[1][j];
389  tempX[2] = Vfrozen[2][j];
390  tempX[3] = Vfrozen[3][j];
391  // pOut.push_back(Parton(0,21,0,newPt,pIn[i].eta(),pIn[i].phi(),newPt));
392 
393  if (std::isnan(tempP[1])) {
394  JSWARN << "second instance: j=" << j << ", P[0][j]=" << P[0][j]
395  << ", P[1][j]=" << P[1][j] << ", P[2][j]=" << P[2][j]
396  << ", P[3][j]=" << P[3][j] << ", P[4][j]=" << P[4][j]
397  << ", P[6][j]=" << P[6][j];
398  }
399 
400  pOut.push_back(Parton(0, KATT1[j], out_stat, tempP, tempX));
401  // remember to put Tint_lrf infomation back to JETSCAPE
402  pOut.back().set_user_info(new LBTUserInfo(Tint_lrf[j]));
403 
404  int iout = pOut.size() - 1;
405  pOut[iout].set_jet_v(velocity_jet); // use initial jet velocity
406  pOut[iout].set_mean_form_time();
407  double ft = pOut[iout].mean_form_time();
408  pOut[iout].set_form_time(ft);
409  }
410 
411  // negative particle stat = -1
412  for (int j = 2; j <= np; j++) {
413  if (P0[0][j] < cutOut)
414  continue;
415  double tempP[4], tempX[4];
416  tempP[0] = sqrt(P0[0][j] * P0[0][j] + P0[6][j]);
417  tempP[1] = P0[1][j];
418  tempP[2] = P0[2][j];
419  tempP[3] = P0[3][j];
420  tempX[0] = Vfrozen0[0][j];
421  tempX[1] = Vfrozen0[1][j];
422  tempX[2] = Vfrozen0[2][j];
423  tempX[3] = Vfrozen0[3][j];
424  // pOut.push_back(Parton(0,21,0,newPt,pIn[i].eta(),pIn[i].phi(),newPt));
425  pOut.push_back(Parton(0, KATT10[j], -1, tempP, tempX));
426  pOut.back().set_user_info(new LBTUserInfo(0.0));
427 
428  int iout = pOut.size() - 1;
429  pOut[iout].set_jet_v(velocity_jet); // use initial jet velocity
430  pOut[iout].set_mean_form_time();
431  double ft = pOut[iout].mean_form_time();
432  pOut[iout].set_form_time(ft);
433  }
434 
435  } else { // free-streaming daughter particles from negative particles if they are inside a medium
436 
437  if (flag_update0 == 0)
438  continue; // no update on this particle from LBT0 function
439 
440  TakeResponsibilityFor(
441  pIn[i]); // Generate error if another module already has responsibility.
442 
443  if (np != 1)
444  JSDEBUG << "Wrong number of negative partons!";
446 
447  for (int j = 1; j <= np; j++) {
448  if (P0[0][j] < cutOut)
449  continue;
450  double tempP[4], tempX[4];
451  tempP[0] = P0[0][j];
452  tempP[1] = P0[1][j];
453  tempP[2] = P0[2][j];
454  tempP[3] = P0[3][j];
455  tempX[0] = Vfrozen0[0][j];
456  tempX[1] = Vfrozen0[1][j];
457  tempX[2] = Vfrozen0[2][j];
458  tempX[3] = Vfrozen0[3][j];
459  // pOut.push_back(Parton(0,21,0,newPt,pIn[i].eta(),pIn[i].phi(),newPt));
460  pOut.push_back(Parton(0, KATT10[j], -1, tempP, tempX));
461  pOut.back().set_user_info(new LBTUserInfo(0.0));
462 
463  int iout = pOut.size() - 1;
464  pOut[iout].set_jet_v(velocity_jet); // use initial jet velocity
465  pOut[iout].set_mean_form_time();
466  double ft = pOut[iout].mean_form_time();
467  pOut[iout].set_form_time(ft);
468  }
469 
470  } // end par_status check
471  }
472 }
473 
475 //
476 // This is the key function of LBT
477 // Update elastic and inelastic evolution of the particle list for one time step
478 //
480 
481 void LBT::LBT0(int &n, double &ti) {
482 
483  int CT; //collision type 1 2 3 13 4 5 6 7 8
484  int KATTC0; //flavor code 1/d 2/u 3/s -1/dbar -2/ubar -3/sbar 21/gluon
485  int KATT2;
486  int KATT3;
487  int KATT30;
488 
489  double RTE; //scattering rate (energy temperature)
490  double E; //parton energy
491  double PLen; //parton momentum
492  double T; //local temperature
493 
494  double T1; //index for difference method
495  double T2;
496  double E1;
497  double E2;
498  int iT1;
499  int iT2;
500  int iE1;
501  int iE2;
502 
503  int nrad;
504  int idlead;
505  int idlead1;
506  int idlead2;
507  double Xtau; //....................main & titau
508  double Vx;
509  double Vy;
510  double Veta;
511 
512  double tcar =
513  0.0; //....................t x y z in Cartesian coordinate this is for the radiation process
514  double xcar = 0.0;
515  double ycar = 0.0;
516  double zcar = 0.0;
517 
518  double tcar0 =
519  0.0; //....................t x y z in Cartesian coordinate this is for the radiation process
520  double xcar0 = 0.0;
521  double ycar0 = 0.0;
522  double zcar0 = 0.0;
523 
524  double rans;
525 
526  int KATTx;
527  int KATTy;
528 
529  double Ejp;
530  double Elab;
531  double Qinit;
532 
533  double qt;
534 
535  double px0;
536  double py0;
537  double pz0;
538 
539  double
540  Ncoll22; //...................average number of elastic scattering for particle i in the current step
541  int Nscatter; //...................number of elastic scattering for particle i in the current step
542 
543  int nnpp = np; //...................number of particles in the current np loop
544  int np0 =
545  np; //...................number of particles in the current step (np0)
546  //...................number of particles in the beginning of current step (np)
547 
548  int free = 0;
549  int free0 = 0;
550  double fraction = 0.0;
551  double fraction0 = 0.0;
552  double vc0b[4] = {0.0}; //flow velocity
553  double pMag, vMag, flowFactor;
554 
555  double kt2;
556 
557  double vp0[4] = {0.0};
558  double p0temp[4] = {0.0};
559 
560  double p0temp1 = 0.0;
561  double p0temp2 = 0.0;
562 
563  double pcx[4] = {0.0};
564  double pcy[4] = {0.0};
565 
566  double pcx1[4] = {0.0};
567  double pcy1[4] = {0.0};
568 
569  // for heavy quark
570  double dt_lrf, maxFncHQ, Tdiff, lim_low, lim_high, lim_int;
571  double probCol, probRad, probTot;
572 
573  double lrf_rStart[4] = {0.0};
574  double SpatialRapidity = 0.0;
575 
576  probCol = 0.0;
577  probRad = 0.0;
578  probTot = 0.0;
579  flowFactor = 0.0;
580 
581  //...........................................................................................................NCUT
582  ncut = 0;
583  ncut0 = 0;
584  //...........................................................................................................NCUTEND
585 
586  //...collision of the jet parton (i<=nj), recoiled partons, radiated gluons with the background thermal partons
587  //...thermal partons
588  //...np number of partons in current step
589  //...nj number of jet partons
590 
591  for (int i = 1; i <= np; i++) {
592 
593  // // if the production time of a parton is ahead of the current time, do nothing
594  // if(Vfrozen[0][i]>=ti) continue;
595 
596  icl22 = 0;
597  qt = 0;
598  idlead = i;
599 
600  //......propagation of each positive particle and its negative counterpart in the current time step (for all particles)
601 
602  if (P[0][i] < epsilon) { //anti-overflow NOT NECESSARY
603  V[1][i] = 0.0;
604  V[2][i] = 0.0;
605  V[3][i] = 0.0;
606  CAT[i] = 1;
607  }
608 
609  if (P0[0][i] < epsilon) { //anti-overflow NOT NECESSARY
610  V0[1][i] = 0.0;
611  V0[2][i] = 0.0;
612  V0[3][i] = 0.0;
613  CAT0[i] = 1;
614  }
615 
616  if (CAT0[i] != 1) {
617  //...........................................t-z coordinates
618  if (tauswitch == 0) {
619 
620  //V0[1][i]=V0[1][i]+dt*P0[1][i]/P0[0][i];
621  //V0[2][i]=V0[2][i]+dt*P0[2][i]/P0[0][i];
622  //V0[3][i]=V0[3][i]+dt*P0[3][i]/P0[0][i];
623 
624  V0[1][i] = V0[1][i] + (ti - Vfrozen0[0][i]) * P0[1][i] / P0[0][i];
625  V0[2][i] = V0[2][i] + (ti - Vfrozen0[0][i]) * P0[2][i] / P0[0][i];
626  V0[3][i] = V0[3][i] + (ti - Vfrozen0[0][i]) * P0[3][i] / P0[0][i];
627 
628  tcar0 = ti;
629  zcar0 = V0[3][i];
630  xcar0 = V0[1][i];
631  ycar0 = V0[2][i];
632 
633  double ed00, sd00, VX00, VY00, VZ00;
634  int hydro_ctl0;
635 
636  free0 = 0;
637 
638  // if(bulkFlag==1) { // read OSU hydro
639  // readhydroinfoshanshan_(&tcar0,&xcar0,&ycar0,&zcar0,&ed00,&sd00,&temp00,&VX00,&VY00,&VZ00,&hydro_ctl0);
640  // } else if(bulkFlag==2) { // read CCNU hydro
641  // hydroinfoccnu_(&tcar0, &xcar0, &ycar0, &zcar0, &temp00, &VX00, &VY00, &VZ00, &hydro_ctl0);
642  // } else if(bulkFlag==0) { // static medium
643 
644  //Extract fluid properties
645  std::unique_ptr<FluidCellInfo> check_fluid_info_ptr;
646 
647  SpatialRapidity = 0.5 * std::log((tcar0 + zcar0) / (tcar0 - zcar0));
648 
649  GetHydroCellSignal(tcar0, xcar0, ycar0, zcar0, check_fluid_info_ptr);
650  //VERBOSE(7)<< MAGENTA<<"Temperature from Brick (Signal) = "<<check_fluid_info_ptr->temperature;
651 
652  temp00 = check_fluid_info_ptr->temperature;
653  sd00 = check_fluid_info_ptr->entropy_density;
654  VX00 = check_fluid_info_ptr->vx;
655  VY00 = check_fluid_info_ptr->vy;
656  VZ00 = check_fluid_info_ptr->vz;
657 
658  if (tcar0 < tStart * cosh(SpatialRapidity)) {
659  temp00 = 0.0;
660  sd00 = 0.0;
661  VX00 = 0.0;
662  VY00 = 0.0;
663  VZ00 = 0.0;
664  }
665 
666  //JSDEBUG << "check temperature: " << temp00;
667 
668  //VX00=0.0;
669  //VY00=0.0;
670  //VZ00=0.0;
671 
672  hydro_ctl0 = 0;
673  // }
674 
675  if (hydro_ctl0 == 0 && temp00 >= hydro_Tc) {
676 
677  qhat00 = DebyeMass2(Kqhat0, alphas, temp00);
678  fraction0 = 1.0;
679 
680  Vfrozen0[0][i] = ti;
681  Vfrozen0[1][i] = V0[1][i];
682  Vfrozen0[2][i] = V0[2][i];
683  Vfrozen0[3][i] = V0[3][i];
684  Tfrozen0[i] = temp00;
685  vcfrozen0[1][i] = VX00;
686  vcfrozen0[2][i] = VY00;
687  vcfrozen0[3][i] = VZ00;
688  flag_update0 =
689  1; // do free-streaming for negative particle if in medium, no check on virtuality (assume 0)
690 
691  } else {
692  free0 = 1;
693  fraction0 = 0.0;
694  }
695 
696  } //if(tauswitch==0)
697 
698  } //if(CAT0[i] != 1)
699 
700  if (CAT[i] != 1) {
701  //...........................................t-z coordinates
702  if (tauswitch == 0) {
703 
704  //V[1][i]=V[1][i]+dt*P[1][i]/P[0][i];
705  //V[2][i]=V[2][i]+dt*P[2][i]/P[0][i];
706  //V[3][i]=V[3][i]+dt*P[3][i]/P[0][i];
707 
708  V[1][i] = V[1][i] + (ti - Vfrozen[0][i]) * P[1][i] / P[0][i];
709  V[2][i] = V[2][i] + (ti - Vfrozen[0][i]) * P[2][i] / P[0][i];
710  V[3][i] = V[3][i] + (ti - Vfrozen[0][i]) * P[3][i] / P[0][i];
711 
712  tcar = ti;
713  zcar = V[3][i];
714  xcar = V[1][i];
715  ycar = V[2][i];
716 
717  for (int lrf_i = 0; lrf_i <= 3; lrf_i++)
718  lrf_rStart[lrf_i] = Vfrozen[lrf_i][i];
719 
720  //...........................................Hydro part
721 
722  double ed, sd, VX, VY, VZ;
723  int hydro_ctl;
724 
725  free = 0;
726 
727  if (free == 0) {
728 
729  // if(bulkFlag==1) { // read OSU hydro
730  // readhydroinfoshanshan_(&tcar,&xcar,&ycar,&zcar,&ed,&sd,&temp0,&VX,&VY,&VZ,&hydro_ctl);
731  // } else if(bulkFlag==2) { // read CCNU hydro
732  // hydroinfoccnu_(&tcar, &xcar, &ycar, &zcar, &temp0, &VX, &VY, &VZ, &hydro_ctl);
733  // } else if(bulkFlag==0) { // static medium
734  //VX=0.0;
735  //VY=0.0;
736  //VZ=0.0;
737 
738  std::unique_ptr<FluidCellInfo> check_fluid_info_ptr;
739 
740  SpatialRapidity = 0.5 * std::log((tcar + zcar) / (tcar - zcar));
741 
742  GetHydroCellSignal(tcar, xcar, ycar, zcar, check_fluid_info_ptr);
743  //VERBOSE(7)<< MAGENTA<<"Temperature from Brick (Signal) = "<<check_fluid_info_ptr->temperature;
744 
745  if (!GetJetSignalConnected()) {
746  JSWARN << "Couldn't find a hydro module attached!";
747  throw std::runtime_error("Please attach a hydro module or set "
748  "in_vac to 1 in the XML file");
749  }
750 
751  temp0 = check_fluid_info_ptr->temperature;
752  sd = check_fluid_info_ptr->entropy_density;
753  VX = check_fluid_info_ptr->vx;
754  VY = check_fluid_info_ptr->vy;
755  VZ = check_fluid_info_ptr->vz;
756 
757  if (tcar < tStart * cosh(SpatialRapidity)) {
758  temp0 = 0.0;
759  sd = 0.0;
760  VX = 0.0;
761  VY = 0.0;
762  VZ = 0.0;
763  }
764 
765  // JSINFO << BOLDYELLOW << "LBT time = " << tcar << " x = " << xcar << " y = " << ycar << " z = " << zcar << " temp = " << temp0;
766  // 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];
767 
768  //JSDEBUG << "check temperature: " << temp0;
769  hydro_ctl = 0;
770  // }
771 
772  if (hydro_ctl == 0 && temp0 >= hydro_Tc) {
773 
774  vc0[1] = VX;
775  vc0[2] = VY;
776  vc0[3] = VZ;
777 
778  // vc0[1]=0.0;
779  // vc0[2]=0.0;
780  // vc0[3]=0.0;
781 
782  vc0b[1] = VX;
783  vc0b[2] = VY;
784  vc0b[3] = VZ;
785 
786  //...alphas!
787  alphas = alphas0(Kalphas, temp0);
788  //...Debye Mass square
789  qhat0 = DebyeMass2(Kqhat0, alphas, temp0);
790 
791  fraction = 1.0;
792  Vfrozen[0][i] = ti;
793  Vfrozen[1][i] = V[1][i];
794  Vfrozen[2][i] = V[2][i];
795  Vfrozen[3][i] = V[3][i];
796  Tfrozen[i] = temp0;
797  vcfrozen[1][i] = VX;
798  vcfrozen[2][i] = VY;
799  vcfrozen[3][i] = VZ;
800 
801  } else {
802  free = 1;
803  fraction = 0.0;
804  }
805 
806  pc0[1] = P[1][i];
807  pc0[2] = P[2][i];
808  pc0[3] = P[3][i];
809  pc0[0] = P[0][i];
810 
811  vMag =
812  sqrt(vc0b[1] * vc0b[1] + vc0b[2] * vc0b[2] + vc0b[3] * vc0b[3]);
813  flowFactor =
814  (1.0 - (pc0[1] * vc0b[1] + pc0[2] * vc0b[2] + pc0[3] * vc0b[3]) /
815  pc0[0]) /
816  sqrt(1.0 - vMag * vMag);
817 
818  } //if(free==0)
819  } //if(tauswitch==0)
820 
821  //...........................................tau-eta coordinates
822  //if(tauswitch==1) { // not ready for JETSCAPE
823  //} // if(tauswitch==1)
824 
825  //......propagation of each positive particle and its negative counterpart in the current time step end
826 
827  //......CAT[i]=1 free streaming particle
828 
829  //......free streaming when energy < Ecut
830  //......free streaming when energy0 < Ecut0......in the next step
831 
832  if (free == 0) {
833 
834  double qhatTP, RTE1, RTE2;
835 
836  // modification: interplotate rate in l.r.f.
837  pc0[1] = P[1][i];
838  pc0[2] = P[2][i];
839  pc0[3] = P[3][i];
840  pc0[0] = P[0][i];
841  trans(vc0, pc0);
842  E = pc0
843  [0]; // p4-the initial 4-momentum of the jet parton in the local rest frame
844  PLen = sqrt(pc0[1] * pc0[1] + pc0[2] * pc0[2] + pc0[3] * pc0[3]);
845  transback(vc0, pc0);
846 
847  T = temp0;
848  KATTC0 = KATT1[i];
849 
850  // if(E<T) preKT=4.0*pi/9.0/log(scaleAK*T*T/0.04)/alphas/fixKT;
851  // else preKT=4.0*pi/9.0/log(scaleAK*E*T/0.04)/alphas/fixKT;
852 
853  if(run_alphas==1) {
854  //runKT=4.0*pi/9.0/log(2.0*E*T/0.04)/0.3;
855  fixedLog = log(5.7*E/4.0/6.0/pi/0.3/T);
856  scaleMu2 = 2.0*E*T;
857  if(scaleMu2 < 1.0) {
858  scaleMu2 = 1.0;
859  runAlphas = alphas;
860  } else {
861  double lambdaQCD2 = exp(-4.0*pi/9.0/alphas);
862  runAlphas = 4.0*pi/9.0/log(scaleMu2/lambdaQCD2);
863  }
864  runKT = runAlphas/0.3;
865  runLog = log(scaleMu2/6.0/pi/T/T/alphas)/fixedLog;
866  }
867 
868  lam(KATTC0, RTE, PLen, T, T1, T2, E1, E2, iT1, iT2, iE1,
869  iE2); //modified: use P instead
870 
871  preKT = alphas / 0.3;
872 
873  // calculate p,T-dependence K factor
874  KPfactor = 1.0 + KPamp * exp(-PLen * PLen / 2.0 / KPsig / KPsig);
875  KTfactor = 1.0 + KTamp * exp(-pow((temp0 - hydro_Tc), 2) / 2.0 / KTsig /
876  KTsig);
877 
878  if(run_alphas==1) {
879  Kfactor = KPfactor * KTfactor * KTfactor * runKT * preKT * runLog; // K factor for qhat
880  } else {
881  Kfactor = KPfactor * KTfactor * KTfactor * preKT * preKT; // K factor for qhat
882  }
883 
884 
885  // get qhat from table
886  if (KATTC0 == 21) {
887  RTE1 = (qhatG[iT2][iE1] - qhatG[iT1][iE1]) * (T - T1) / (T2 - T1) +
888  qhatG[iT1][iE1];
889  RTE2 = (qhatG[iT2][iE2] - qhatG[iT1][iE2]) * (T - T1) / (T2 - T1) +
890  qhatG[iT1][iE2];
891  } else if (KATTC0 == 4 || KATTC0 == -4 || KATTC0 == 5 || KATTC0 == -5) {
892  RTE1 = (qhatHQ[iT2][iE1] - qhatHQ[iT1][iE1]) * (T - T1) / (T2 - T1) +
893  qhatHQ[iT1][iE1];
894  RTE2 = (qhatHQ[iT2][iE2] - qhatHQ[iT1][iE2]) * (T - T1) / (T2 - T1) +
895  qhatHQ[iT1][iE2];
896  } else {
897  RTE1 = (qhatLQ[iT2][iE1] - qhatLQ[iT1][iE1]) * (T - T1) / (T2 - T1) +
898  qhatLQ[iT1][iE1];
899  RTE2 = (qhatLQ[iT2][iE2] - qhatLQ[iT1][iE2]) * (T - T1) / (T2 - T1) +
900  qhatLQ[iT1][iE2];
901  }
902 
903  qhatTP = (RTE2 - RTE1) * (PLen - E1) / (E2 - E1) + RTE1;
904 
905  qhatTP = qhatTP * Kfactor;
906 
911  // qhatTP=14.80;
912  // cout << "RTE: " << RTE << endl;
913 
914  qhat_over_T3 = qhatTP; // what is read in is qhat/T^3 of quark/gluon
915 
916  // D2piT is diffusion coefficient "just for quark"
917  if (KATTC0 == 21)
918  D2piT = 8.0 * pi / (qhat_over_T3 / CA * CF);
919  else
920  D2piT = 8.0 * pi / qhat_over_T3;
921 
922  qhat =
923  qhatTP *
924  pow(T,
925  3); // for light quark and gluon, need additional table to make it correct
926 
927  if (Q00 < 0.0) {
928  Q0 = sqrt(sqrt(2.0 * E * qhat));
929  } else {
930  Q0 = Q00;
931  }
932 
933  if (P[6][i] > Q0 * Q0 + rounding_error) {
934  continue; // virtuality outside the LBT range, go to the next particle
935  }
936 
937  flag_update =
938  1; // satisfy T>Tc and Q<Q0, LBT will process this positice particle
939 
940  if (E < sqrt(qhat0))
941  continue; // just do free-streaming
942  if (i > nj && E < Ecmcut)
943  continue;
944 
945  // update interaction time and average gluon number for heavy quark radiation
946  // (as long as it's in medium, even though no collision)
947  dt_lrf =
948  dt *
949  flowFactor; // must be put here, flowFactor will be changed below
950 
951  // increae time accumulation from the production point of the parton
952  for (double tLoc = lrf_rStart[0]; tLoc < ti + rounding_error;
953  tLoc = tLoc + dt) {
954 
955  double xLoc =
956  lrf_rStart[1] + (tLoc - lrf_rStart[0]) * P[1][i] / P[0][i];
957  double yLoc =
958  lrf_rStart[2] + (tLoc - lrf_rStart[0]) * P[2][i] / P[0][i];
959  double zLoc =
960  lrf_rStart[3] + (tLoc - lrf_rStart[0]) * P[3][i] / P[0][i];
961  double dtLoc, tempLoc, vxLoc, vyLoc, vzLoc;
962 
963  if (tLoc + dt < ti)
964  dtLoc = dt;
965  else
966  dtLoc = ti - tLoc;
967 
968  std::unique_ptr<FluidCellInfo> check_fluid_info_ptr;
969 
970  SpatialRapidity = 0.5 * std::log((tLoc + zLoc) / (tLoc - zLoc));
971 
972  GetHydroCellSignal(tLoc, xLoc, yLoc, zLoc, check_fluid_info_ptr);
973  tempLoc = check_fluid_info_ptr->temperature;
974  vxLoc = check_fluid_info_ptr->vx;
975  vyLoc = check_fluid_info_ptr->vy;
976  vzLoc = check_fluid_info_ptr->vz;
977 
978  if (tLoc < tStart * cosh(SpatialRapidity)) {
979  tempLoc = 0.0;
980  vxLoc = 0.0;
981  vyLoc = 0.0;
982  vzLoc = 0.0;
983  }
984 
985  if (tempLoc >= hydro_Tc) {
986  vMag = sqrt(vxLoc * vxLoc + vyLoc * vyLoc + vzLoc * vzLoc);
987  flowFactor =
988  (1.0 -
989  (pc0[1] * vxLoc + pc0[2] * vyLoc + pc0[3] * vzLoc) / pc0[0]) /
990  sqrt(1.0 - vMag * vMag);
991  Tint_lrf[i] = Tint_lrf[i] + dtLoc * flowFactor;
992  }
993  }
994 
995  Tdiff = Tint_lrf[i];
996 
997  // get radiation probablity by reading tables -- same for heavy and light partons
998  if (KATTC0 == 21) {
999  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  //??????@@@
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  }
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 
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 
1436 //
1437 // Define functions for LBT
1438 //
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 
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;
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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
4416  //cout << endl;
4417  //cout << "############################################################" << endl;
4418  //cout << "Parameters for this LBT run" << 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;
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;
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 
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 }