Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Matter.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Matter.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 "Matter.h"
17 #include "JetScapeLogger.h"
18 #include "JetScapeParticles.h"
19 #include "Pythia8/Pythia.h"
20 
21 #include <string>
22 
23 #include <iostream>
24 
25 #include "FluidDynamics.h"
26 #include <GTL/dfs.h>
27 
28 #define MAGENTA "\033[35m"
29 
30 using namespace Jetscape;
31 using namespace std;
32 
33 const double QS = 0.9;
34 
35 // Register the module with the base class
36 RegisterJetScapeModule<Matter> Matter::reg("Matter");
37 
38 static Pythia8::Pythia PythiaFunction("IntentionallyEmpty", false);
39 
40 bool Matter::flag_init = 0;
41 
42 double Matter::RHQ[60][20] = {{0.0}}; //total scattering rate for heavy quark
43 double Matter::RHQ11[60][20] = {{0.0}}; //Qq->Qq
44 double Matter::RHQ12[60][20] = {{0.0}}; //Qg->Qg
45 double Matter::qhatHQ[60][20] = {{0.0}}; //qhat of heavy quark
46 
47 double Matter::distFncB[N_T][N_p1][N_e2] = {{{0.0}}};
48 double Matter::distFncF[N_T][N_p1][N_e2] = {{{0.0}}};
49 double Matter::distMaxB[N_T][N_p1][N_e2] = {{{0.0}}};
50 double Matter::distMaxF[N_T][N_p1][N_e2] = {{{0.0}}};
51 double Matter::distFncBM[N_T][N_p1] = {{0.0}};
52 double Matter::distFncFM[N_T][N_p1] = {{0.0}};
53 
54 Matter::Matter() {
55  SetId("Matter");
56  VERBOSE(8);
57  qhat = 0.0;
58  ehat = 0.0;
59  e2hat = 0.0;
60  length = 0.0;
61  MaxColor = 0;
62  matter_on = true;
63  in_vac = false;
64  brick_med = true;
65  recoil_on = false;
66  broadening_on = false;
67  hydro_Tc = 0.;
68  qhat0 = 0.;
69  alphas = 0.;
70  tscale=1;
71  QhatParametrizationType=-1;
72  qhatA=0.;
73  qhatB=0.;
74  qhatC=0.;
75  qhatD=0.;
76  brick_length = 0.;
77  vir_factor = 0.;
78  initR0 = 0.;
79  initRx = 0.;
80  initRy = 0.;
81  initRz = 0.;
82  initVx = 0.;
83  initVy = 0.;
84  initVz = 0.;
85  initRdotV = 0.;
86  initVdotV = 0.;
87  initEner = 0.;
88  Q00 = 0.;
89  Q0 = 0.;
90  T0 = 0.;
91  iEvent = 0;
92  NUM1 = 0;
93 }
94 
95 Matter::~Matter() { VERBOSE(8); }
96 
97 void Matter::Init() {
98  JSINFO << "Initialize Matter ...";
99 
100  in_vac = false;
101  brick_med = true;
102  recoil_on = false;
103 
104  qhat = 0.0;
105  Q00 = 1.0; // virtuality separation scale
106  qhat0 = 2.0; // GeV^2/fm for gluon at s = 96 fm^-3
107  alphas = 0.3; // only useful when qhat0 is a negative number
108  tscale=1;
109  QhatParametrizationType=-1;
110  qhatA=1;
111  qhatB=1;
112  qhatC=1;
113  qhatD=1;
114  hydro_Tc = 0.16;
115  brick_length = 4.0;
116  vir_factor = 1.0;
117 
118  double m_qhat = GetXMLElementDouble({"Eloss", "Matter", "qhat0"});
119  SetQhat(m_qhat);
120  //qhat = GetQhat()/fmToGeVinv ;
121  //qhat0 = GetQhat()/fmToGeVinv ;
122  qhat0 = GetQhat();
123 
124  int flagInt = -100;
125  double inputDouble = -99.99;
126 
127  matter_on = GetXMLElementInt({"Eloss", "Matter", "matter_on"});
128  in_vac = GetXMLElementInt({"Eloss", "Matter", "in_vac"});
129  recoil_on = GetXMLElementInt({"Eloss", "Matter", "recoil_on"});
130  broadening_on = GetXMLElementInt({"Eloss", "Matter", "broadening_on"});
131  brick_med = GetXMLElementInt({"Eloss", "Matter", "brick_med"});
132  Q00 = GetXMLElementDouble({"Eloss", "Matter", "Q0"});
133  T0 = GetXMLElementDouble({"Eloss", "Matter", "T0"});
134  alphas = GetXMLElementDouble({"Eloss", "Matter", "alphas"});
135  qhatA = GetXMLElementDouble({"Eloss", "Matter", "qhatA"});
136  qhatB = GetXMLElementDouble({"Eloss", "Matter", "qhatB"});
137  qhatC = GetXMLElementDouble({"Eloss", "Matter", "qhatC"});
138  qhatD = GetXMLElementDouble({"Eloss", "Matter", "qhatD"});
139  tStart = GetXMLElementDouble({"Eloss", "tStart"});
140  //run_alphas = GetXMLElementInt({"Eloss", "Matter", "run_alphas"});
141  QhatParametrizationType=GetXMLElementInt({"Eloss", "Matter", "QhatParametrizationType"});
142  hydro_Tc = GetXMLElementDouble({"Eloss", "Matter", "hydro_Tc"});
143  brick_length = GetXMLElementDouble({"Eloss", "Matter", "brick_length"});
144  vir_factor = GetXMLElementDouble({"Eloss", "Matter", "vir_factor"});
145 
146  if (vir_factor < 0.0) {
147  cout << "Reminder: negative vir_factor is set, initial energy will be used "
148  "as initial t_max"
149  << endl;
150  }
151 
152  MaxColor = 101; // MK:recomb
153 
154  JSINFO << MAGENTA << "MATTER input parameter";
155  JSINFO << MAGENTA << "matter shower on: " << matter_on;
156  JSINFO << MAGENTA << "in_vac: " << in_vac << " brick_med: " << brick_med
157  << " recoil_on: " << recoil_on<<", tStart ="<<tStart;
158  JSINFO << MAGENTA << "Q0: " << Q00 << " vir_factor: " << vir_factor
159  << " qhat0: " << qhat0 << " alphas: " << alphas << ", QhatParametrizationType="<<QhatParametrizationType
160  << " qhatA: " << qhatA << " qhatB: " <<qhatB << " qhatC: " << qhatC << " qhatD: " <<qhatD
161  << " hydro_Tc: " << hydro_Tc << " brick_length: " << brick_length;
162  if(QhatParametrizationType ==1){ JSINFO <<"Running alphas will be used as 4*pi/(9.0*log(2*enerLoc*tempLoc/LambdaQCDS))"; }
163 
164  if (recoil_on && !flag_init) {
165  JSINFO << MAGENTA
166  << "Reminder: download LBT tables first and cmake .. if recoil is "
167  "switched on in MATTER.";
168  read_tables(); // initialize various tables
169  flag_init = true;
170  }
171 
172  // Initialize random number distribution
173  ZeroOneDistribution = uniform_real_distribution<double>{0.0, 1.0};
174 
175  //...initialize the random number generator
176  srand((unsigned)time(NULL));
177  NUM1 = -1 * rand();
178  // NUM1=-33;
179  iEvent = 0;
180 }
181 
182 void Matter::WriteTask(weak_ptr<JetScapeWriter> w) {
183  VERBOSE(8);
184  auto f = w.lock();
185  if (!f)
186  return;
187  f->WriteComment("ElossModule Parton List: " + GetId());
188  f->WriteComment("Energy loss to be implemented accordingly ...");
189 }
190 
191 void Matter::Dump_pIn_info(int i, vector<Parton> &pIn) {
192  JSWARN << "i=" << i << " MATTER -- status: " << pIn[i].pstat()
193  << " color: " << pIn[i].color() << " " << pIn[i].anti_color();
194  JSWARN << "pid = " << pIn[i].pid() << " E = " << pIn[i].e()
195  << " px = " << pIn[i].p(1) << " py = " << pIn[i].p(2)
196  << " pz = " << pIn[i].p(3) << " virtuality = " << pIn[i].t()
197  << " form_time in fm = " << pIn[i].form_time()
198  << " split time = " << pIn[i].form_time() + pIn[i].x_in().t();
199 }
200 
201 void Matter::DoEnergyLoss(double deltaT, double time, double Q2,
202  vector<Parton> &pIn, vector<Parton> &pOut) {
203 
204  if (std::isnan(pIn[0].e()) || std::isnan(pIn[0].px()) ||
205  std::isnan(pIn[0].py()) || std::isnan(pIn[0].pz()) ||
206  std::isnan(pIn[0].t()) || std::isnan(pIn[0].form_time())) {
207  JSINFO << BOLDYELLOW << "Parton on entry busted on time step " << time;
208  Dump_pIn_info(0, pIn);
209  }
210 
211  double z = 0.5;
212  double blurb, zeta, tQ2;
213  int iSplit, pid_a, pid_b;
214  unsigned int max_color, min_color, min_anti_color;
215  double velocity[4], xStart[4], velocity_jet[4];
216  bool photon_brem = false;
217 
218  VERBOSE(8) << MAGENTA << "SentInPartons Signal received : " << deltaT << " "
219  << Q2 << " " << pIn.size();
220 
221  VERBOSE(8) << BOLDYELLOW
222  << " ********************************************************** ";
223 
224  double rNum;
225 
226  double delT = deltaT;
227  double Time = time * fmToGeVinv;
228  double deltaTime = delT * fmToGeVinv;
229  double ehat = 0;
230  double ehat_over_T2 = 10.0;
231 
232  std::unique_ptr<FluidCellInfo> check_fluid_info_ptr;
233 
234  VERBOSE(8) << MAGENTA << " the time in fm is " << time
235  << " The time in GeV-1 is " << Time;
236  VERBOSE(8) << MAGENTA << "pid = " << pIn[0].pid() << " E = " << pIn[0].e()
237  << " px = " << pIn[0].p(1) << " py = " << pIn[0].p(2)
238  << " pz = " << pIn[0].p(3) << " virtuality = " << pIn[0].t()
239  << " form_time in fm = " << pIn[0].form_time()
240  << " split time = " << pIn[0].form_time() + pIn[0].x_in().t();
241  VERBOSE(8) << " color = " << pIn[0].color()
242  << " anti-color = " << pIn[0].anti_color();
243 
244  unsigned int ShowerMaxColor = pIn[0].max_color();
245  unsigned int CurrentMaxColor;
246 
247  if (pIn[0].max_color() < MaxColor) {
248  pIn[0].set_max_color(MaxColor);
249 
250  } else {
251  MaxColor = pIn[0].max_color();
252  }
253 
254  //VERBOSE(8) << MAGENTA << " Max color = " << MaxColor;
255  //JSDEBUG << " For MATTER, the qhat in GeV^-3 = " << qhat ;
256 
257  double qhatbrick;
258  if (brick_med)
259  qhatbrick = qhat0 / 3.0;
260  qhat = qhat0;
261 
262  VERBOSE(8) << " qhat0 = " << qhat0 << " qhat = " << qhat;
263 
264  for (int i = 0; i < pIn.size(); i++) {
265 
266  // Reject photons
267  if (pIn[i].pid() == photonid) {
268  if(pIn[i].pstat() != 22) {
269  pIn[i].set_stat(22);
270  VERBOSE(1) << BOLDYELLOW
271  << " A photon was RECEIVED with px = " << pIn[i].px()
272  << " from framework and sent back ";
273 
274  pOut.push_back(pIn[i]);
275  }
276  return;
277  }
278 
279  // Reject photons
280 
281  if (std::abs(pIn[i].pstat()) == 1) {
282 
283  // JSINFO << BOLDYELLOW << " A recoil was RECEIVED with px = " << pIn[i].px() << " py = " << pIn[i].py() << " pz = " << pIn[i].pz() << " E = " << pIn[i].e() << " from framework and sent back " ;
284  // JSINFO << BOLDYELLOW << "t=" << " * parton formation spacetime point= "<< pIn[i].x_in().t() << " " << pIn[i].x_in().x() << " " << pIn[i].x_in().y() << " " << pIn[i].x_in().z();
285  // Dump_pIn_info(i,pIn);
286 
287  // pOut.push_back(pIn[i]);
288 
289  return;
290  }
291 
292  VERBOSE(2) << BOLDYELLOW
293  << " * parton formation spacetime point= " << pIn[i].x_in().t()
294  << " " << pIn[i].x_in().x() << " " << pIn[i].x_in().y() << " "
295  << pIn[i].x_in().z();
296 
297  //JSINFO << MAGENTA << " particle rest mass = " << pIn[i].restmass();
298 
299  int jet_stat = pIn[i].pstat(); // daughter of recoil will always be recoil
300 
301  JSDEBUG << MAGENTA << "MATTER -- status: " << pIn[i].pstat()
302  << " energy: " << pIn[i].e() << " color: " << pIn[i].color()
303  << " " << pIn[i].anti_color() << " clock: " << time;
304 
305  velocity[0] = 1.0;
306  // Define 3 velocity of the parton
307  for (int j = 1; j <= 3; j++) {
308  velocity[j] = pIn[i].p(j) / pIn[i].e();
309  }
310  // velocityMod will be 1 for most partons except maybe heavy quarks, we say "maybe", as for very
311  // energetic heavy quarks, the velocity may be very close to 1.
312  double velocityMod =
313  std::sqrt(std::pow(velocity[1], 2) + std::pow(velocity[2], 2) +
314  std::pow(velocity[3], 2));
315 
316  if (velocityMod > 1.0 + rounding_error) {
317  JSINFO << BOLDRED
318  << " tachyonic propagation detected for parton passed from hard "
319  "scattering, velocity mod = "
320  << velocityMod;
321  JSWARN << "velocityMod=" << std::setprecision(20) << velocityMod;
322  Dump_pIn_info(i, pIn);
323  //assert(velocityMod < 1.0 + rounding_error);
324  }
325  VERBOSE(2) << BOLDYELLOW << " velocityMod = " << velocityMod;
326 
327  if (pIn[i].form_time() < 0.0)
328  pIn[i].set_jet_v(velocity); // jet velocity is set only once
329  // Notice the assumption that partons passed from hard scattering are on shell.
330  // If they had virtuality then this would not be correct.
331 
332  // Define a vector in the direction of the jet originating parton.
333  // there is some amount of redundancy here, pIn[i].jet_v() is basically the jet velocity
334  // we are defining a local copy in velocity_jet
335  velocity_jet[0] = 1.0;
336  velocity_jet[1] = pIn[i].jet_v().x();
337  velocity_jet[2] = pIn[i].jet_v().y();
338  velocity_jet[3] = pIn[i].jet_v().z();
339 
340  // Modulus of the vector pIn[i].jet_v
341  double mod_jet_v =
342  std::sqrt(pow(pIn[i].jet_v().x(), 2) + pow(pIn[i].jet_v().y(), 2) +
343  pow(pIn[i].jet_v().z(), 2));
344 
345  for (int j = 0; j <= 3; j++) {
346  xStart[j] = pIn[i].x_in().comp(j);
347  }
348 
349  // SC: read in hydro
350  initR0 = xStart[0];
351  initRx = xStart[1];
352  initRy = xStart[2];
353  initRz = xStart[3];
354  initVx = velocity[1] / velocityMod;
355  initVy = velocity[2] / velocityMod;
356  initVz = velocity[3] / velocityMod;
357 
358  if (std::abs(pIn[i].pid()) == 4 || std::abs(pIn[i].pid()) == 5) {
359  double OnShellEnergy = std::sqrt(
360  pIn[i].px() * pIn[i].px() + pIn[i].py() * pIn[i].py() +
361  pIn[i].pz() * pIn[i].pz() + pIn[i].restmass() * pIn[i].restmass());
362 
363  initVx = pIn[i].px() / OnShellEnergy;
364  initVy = pIn[i].py() / OnShellEnergy;
365  initVz = pIn[i].pz() / OnShellEnergy;
366 
367  velocityMod =
368  std::sqrt(initVx * initVx + initVy * initVy + initVz * initVz);
369  }
370 
371  initRdotV = (initRx * pIn[i].jet_v().x() + initRy * pIn[i].jet_v().y() +
372  initRz * pIn[i].jet_v().z()) /
373  mod_jet_v;
374  initVdotV = (initVx * pIn[i].jet_v().x() + initVy * pIn[i].jet_v().y() +
375  initVz * pIn[i].jet_v().z()) /
376  mod_jet_v;
377  // Note: jet_v()/mod_jet_v is a unit 3 vector in the direction of the jet originating parton.
378 
379  double now_R0 = time;
380  double now_Rx = initRx + (time - initR0) * initVx;
381  double now_Ry = initRy + (time - initR0) * initVy;
382  double now_Rz = initRz + (time - initR0) * initVz;
383  double now_temp;
384 
385  double SpatialRapidity =
386  0.5 * std::log((now_R0 + now_Rz) / (now_R0 - now_Rz));
387 
388  //JSINFO << MAGENTA << " Particle Rapidity = " << SpatialRapidity ;
389 
390  if (std::isnan(velocityMod) || std::isnan(velocity[1]) ||
391  std::isnan(velocity[2]) || std::isnan(velocity[3]) ||
392  std::isinf(velocityMod) || std::isinf(velocity[1]) ||
393  std::isinf(velocity[2]) || std::isinf(velocity[3])) {
394  JSINFO << BOLDYELLOW << "Zeroth instance";
395  JSINFO << BOLDYELLOW << "time, initR0, initRx, initRy, initRz=" << time
396  << ", " << initR0 << ", " << initRx << ", " << initRy << ", "
397  << initRz;
398  JSINFO << BOLDYELLOW << "Vx, Vy, Vz =" << velocity[1] << ", "
399  << velocity[2] << ", " << velocity[3];
400  JSINFO << BOLDYELLOW << "initVx, initVy, initVz =" << initVx << ", "
401  << initVy << ", " << initVz;
402  Dump_pIn_info(i, pIn);
403  }
404 
405  initEner = pIn[i].e(); // initial Energy of parton
406  if (!in_vac) {
407  if (GetJetSignalConnected())
408  length = fillQhatTab(SpatialRapidity);
409  else {
410  JSWARN << "Couldn't find a hydro module attached!";
411  throw std::runtime_error(
412  "Please attach a hydro module or set in_vac to 1 in the XML file");
413  }
414  }
415  if (brick_med)
416  length = brick_length *
417  fmToGeVinv;
418  //if(brick_med) length = 5.0*fmToGeVinv; /// length in GeV-1 will have to changed for hydro
419 
420  // SC
421  zeta = ((xStart[0] + initRdotV) / std::sqrt(2)) * fmToGeVinv;
422 
423  VERBOSE(8) << BOLDYELLOW << " zeta = " << zeta;
424 
425  // if(now_R0^2-now_Ri^2<0) print out pIn info and exit
426 
427  if (std::isinf(now_R0) || std::isnan(now_R0) || std::isinf(now_Rz) ||
428  std::isnan(now_Rz) || std::abs(now_Rz) > now_R0) {
429  JSINFO << BOLDYELLOW << "First instance";
430  JSINFO << BOLDYELLOW << "now_R for vector is:" << now_R0 << ", " << now_Rx
431  << ", " << now_Ry << ", " << now_Rz;
432  JSINFO << BOLDYELLOW << "time, initR0, initRx, initRy, initRz=" << time
433  << ", " << initR0 << ", " << initRx << ", " << initRy << ", "
434  << initRz;
435  JSINFO << BOLDYELLOW << "initVx, initVy, initVz =" << initVx << ", "
436  << initVy << ", " << initVz;
437  JSINFO << BOLDYELLOW << "velocityMod=" << velocityMod;
438  JSINFO << BOLDYELLOW << "initVMod="
439  << std::sqrt(initVx * initVx + initVy * initVy + initVz * initVz);
440  Dump_pIn_info(i, pIn);
441  //exit(0);
442  }
443  double boostedTStart = tStart * cosh(SpatialRapidity);
444  if (!in_vac && now_R0 >= boostedTStart) {
445  if (now_R0 * now_R0 < now_Rz * now_Rz)
446  cout << "Warning 1: " << now_R0 << " " << now_Rz << endl;
447  GetHydroCellSignal(now_R0, now_Rx, now_Ry, now_Rz, check_fluid_info_ptr);
448  //VERBOSE(8)<<MAGENTA<<"Temperature from medium = "<<check_fluid_info_ptr->temperature;
449  now_temp = check_fluid_info_ptr->temperature;
450  //JSINFO << BOLDYELLOW << "MATTER time = " << now_R0 << " x = " << now_Rx << " y = " << now_Ry << " z = " << now_Rz << " temp = " << now_temp;
451  //JSINFO << BOLDYELLOW << "MATTER initVx, initVy, initVz =" << initVx << ", " << initVy << ", " << initVz;
452  //JSINFO << BOLDYELLOW << "MATTER velocityMod=" << velocityMod;
453  } else {
454  now_temp = 0.0;
455  }
456 
457  int pid = pIn[i].pid();
458 
459  if (pIn[i].form_time() <
460  0.0)
461  {
462 
463  VERBOSE(8) << BOLDYELLOW << " pid = " << pIn[i].pid()
464  << " E = " << pIn[i].e();
465 
466  if ((pIn[i].t() < 0.0) &&
467  ((pIn[i].form_time() < -0.1 - rounding_error) ||
468  (pIn[i].form_time() > -0.1 + rounding_error))) {
469  JSWARN << " parton with a negative virtuality was sent to MATTER and "
470  "will now have its virtuality reset!, press 1 and return to "
471  "proceed... ";
472  // cin >> blurb; //remove the input to prevent an error caused by heavy quark from pythia (by Chathuranga)
473  }
474 
475  iSplit = 0; // (anti)quark splitting into (anti)quark + gluon
476  if (pIn[i].pid() == gid) {
477  JSDEBUG << " parton is a gluon ";
478  iSplit = 1; // gluon splitting to two gluons
479  } else {
480  JSDEBUG << " parton is a quark ";
481  }
482 
483  //tQ2 = generate_vac_t(pIn[i].pid(), pIn[i].nu(), QS/2.0, pIn[i].e()*pIn[i].e() ,zeta , iSplit);
484 
485  // SC:
486  double pT2 = pIn[i].p(1) * pIn[i].p(1) + pIn[i].p(2) * pIn[i].p(2);
487  double max_vir;
488  if (vir_factor < 0.0)
489  max_vir =
490  std::abs(vir_factor) * (pIn[i].e() * pIn[i].e() - pIn[i].restmass() * pIn[i].restmass());
491  else
492  max_vir = pT2 * vir_factor;
493 
494  if (max_vir <= QS * QS) {
495  tQ2 = 0.0;
496  } else {
497  JSDEBUG << BOLDYELLOW << " at x,y,z,t = " << pIn[i].x_in().x() << " "
498  << pIn[i].x_in().y() << " " << pIn[i].x_in().z() << " "
499  << pIn[i].x_in().t();
500  if (abs(pIn[i].pid()) == 4 || abs(pIn[i].pid()) == 5) {
501 
502  double min_vir =
503  (QS * QS / 2.0) *
504  (1.0 + std::sqrt(1.0 + 4.0 * pIn[i].restmass() *
505  pIn[i].restmass() / QS / QS));
506  if (max_vir > min_vir) {
507  tQ2 =
508  generate_vac_t_w_M(pIn[i].pid(), pIn[i].restmass(), pIn[i].nu(),
509  QS * QS / 2.0, max_vir, zeta, iSplit);
510  } else {
511  tQ2 = QS * QS;
512  }
513  // std::ofstream tdist;
514  // tdist.open("tdist_heavy.dat", std::ios::app);
515  // tdist << tQ2 << endl;
516  // tdist.close();
517 
518  VERBOSE(8) << BOLDYELLOW << " virtuality calculated as = " << tQ2;
519  } else if (pIn[i].pid() == gid) {
520  tQ2 = generate_vac_t_w_M(pIn[i].pid(), pIn[i].restmass(), pIn[i].nu(),
521  QS * QS / 2.0, max_vir, zeta, iSplit);
522  } else {
523  tQ2 = generate_vac_t(pIn[i].pid(), pIn[i].nu(), QS * QS / 2.0,
524  max_vir, zeta, iSplit);
525  }
526  }
527 
528  // SC: if matter_on = false, set zero virtuality and MATTER will not do parton shower
529  if (matter_on) {
530  pIn[i].set_t(tQ2); // Also resets momentum!
531  VERBOSE(8) << BOLDYELLOW << " virtuality set to " << tQ2
532  << " max virtuality allowed = " << max_vir;
533  VERBOSE(8) << BOLDYELLOW << " ID = " << pIn[i].pid()
534  << " Color = " << pIn[i].color()
535  << " Anti-Color = " << pIn[i].anti_color();
536  VERBOSE(8) << BOLDYELLOW << " E = " << pIn[i].e()
537  << " px = " << pIn[i].px() << " py = " << pIn[i].py()
538  << " pz = " << pIn[i].pz();
539 
540  } else
541  pIn[i].set_t(rounding_error);
542  //else pIn[i].set_t(0.0);
543 
544  pIn[i].set_mean_form_time();
545  double ft = generate_L(pIn[i].mean_form_time());
546  pIn[i].set_form_time(ft);
547 
548 
549  pIn[i].set_min_color(pIn[i].color());
550  pIn[i].set_min_anti_color(pIn[i].anti_color());
551  MaxColor = pIn[i].max_color();
552 
553 
554  // VERBOSE OUTPUT ON INITIAL STATUS OF PARTICLE:
555  VERBOSE(8);
556  VERBOSE(8) << " *********************************************************"
557  "******************** ";
558  VERBOSE(8) << BOLDYELLOW << " ID = " << pIn[i].pid()
559  << " Color = " << pIn[i].color()
560  << " Anti-Color = " << pIn[i].anti_color();
561  VERBOSE(8) << BOLDYELLOW << " E = " << pIn[i].e()
562  << " px = " << pIn[i].px() << " py = " << pIn[i].py()
563  << " pz = " << pIn[i].pz();
564  VERBOSE(8) << BOLDYELLOW << " * New generated virtuality = " << tQ2
565  << " Mean formation time = " << pIn[i].mean_form_time();
566  VERBOSE(8) << BOLDYELLOW << " * set new formation time to "
567  << pIn[i].form_time();
568  VERBOSE(8) << BOLDYELLOW << " * Maximum allowed virtuality = "
569  << pIn[i].e() * pIn[i].e() -
570  pIn[i].restmass() * pIn[i].restmass()
571  << " Minimum Virtuality = " << QS * QS;
572  VERBOSE(8) << " * Qhat = " << qhat << " Length in fm = " << length / 5.0;
573  VERBOSE(8) << " * Jet velocity = " << pIn[i].jet_v().comp(0) << " "
574  << pIn[i].jet_v().comp(1) << " " << pIn[i].jet_v().comp(2)
575  << " " << pIn[i].jet_v().comp(3);
576  VERBOSE(8) << " * reset location of parton formation = "
577  << pIn[i].x_in().t() << " " << pIn[i].x_in().x() << " "
578  << pIn[i].x_in().y() << " " << pIn[i].x_in().z();
579  VERBOSE(8) << " *********************************************************"
580  "******************** ";
581  VERBOSE(8);
582  // end VERBOSE OUTPUT:
583  }
584 
585  // SC: Q0 can be changed based on different setups
586  if (in_vac) { // for vaccuum
587  qhat = 0.0;
588  if (Q00 < 0.0)
589  Q0 = 1.0; // set Q0 = 1 if Q00 < 0
590  else
591  Q0 = Q00;
592  } else { // for medium
593  double tempEner = initEner;
594  qhat = fncQhat(zeta);
595  ehat = 0.0;
596 
597  if (now_temp > 0.0)
598  ehat = 0.0 * qhat / 4.0 / now_temp;
599  VERBOSE(8) << BOLDYELLOW << "at Origin of parton, qhat = " << qhat
600  << " ehat = " << ehat;
601 
602  // set the minimum Q0 for MATTER using Q^2 = qhat * tau ELSE use the positive value in XML
603  if (Q00 < 0.0) { // use dynamical Q0 if Q00 < 0
604  if (pid == gid)
605  Q0 = sqrt(sqrt(2.0 * tempEner * qhat * sqrt(2.0)));
606  else
607  Q0 = sqrt(sqrt(2.0 * tempEner * qhat * sqrt(2.0) / Ca * Cf));
608  if (Q0 < 1.0)
609  Q0 = 1.0;
610  if (zeta > length)
611  Q0 = 1.0;
612  } else {
613  Q0 = Q00;
614  }
615  }
616 
617  // I dont care what you say, we are not doing pQCD below 1 GeV
618  if (Q0 < 1.0)
619  Q0 = 1.0;
620 
621  //if (pIn[i].t() > QS + rounding_error)
622  if (pIn[i].t() > Q0 * Q0 + rounding_error ||
623  ((!in_vac) && now_temp <= T0 &&
624  pIn[i].t() > QS * QS + rounding_error)) {
625 
626  TakeResponsibilityFor(
627  pIn[i]); // Generate error if another module already has responsibility.
628  double decayTime = pIn[i].mean_form_time();
629 
630  //JSDEBUG << " deltaT = " << deltaT;
631  VERBOSE(8) << " parton origin time = " << pIn[i].x_in().t()
632  << " parton formation time = " << pIn[i].form_time();
633  VERBOSE(8) << " parton id " << pIn[i].pid()
634  << " parton virtuality = " << pIn[i].t();
635  //JSDEBUG << " parton momentum " << pIn[i].e() << " " << pIn[i].px() << " " << pIn[i].py() << " " << pIn[i].pz();
636 
637  double splitTime = pIn[i].form_time() + pIn[i].x_in().t();
638 
639  VERBOSE(8) << " splitTime = " << splitTime;
640  VERBOSE(8) << " qhat before splitime loop = " << qhat;
641 
642  if (splitTime <
643  time) // it is time to split and calculate the effect of scattering
644  {
645 
646  VERBOSE(8) << "SPLIT in MATTER";
647 
648  // SC: add elastic scattering that generates recoiled and back-reaction partons
649  double el_dt = 0.1;
650  double el_p0[5];
651  double el_CR;
652  double el_rand;
653  double HQ_mass;
654 
655  if (pid == gid)
656  el_CR = Ca;
657  else
658  el_CR = Cf;
659 
660  for (int j = 1; j <= 3; j++)
661  el_p0[j] = pIn[i].p(j);
662  el_p0[0] = pIn[i].e();
663  el_p0[4] = pIn[i].t();
664  HQ_mass = pIn[i].restmass();
665 
666  for (double el_time = initR0; el_time < time + rounding_error;
667  el_time = el_time + el_dt) {
668 
669  if (in_vac)
670  continue;
671  if (!recoil_on)
672  continue;
673 
674  double boostedTStart = tStart * std::cosh(SpatialRapidity);
675  if (el_time < boostedTStart)
676  continue;
677 
678  if (abs(pid) == 5)
679  continue; // recoil not ready for b quark yet
680 
681  double el_rx = initRx + (el_time - initR0) * initVx;
682  double el_ry = initRy + (el_time - initR0) * initVy;
683  double el_rz = initRz + (el_time - initR0) * initVz;
684 
685  double tempLoc, sdLoc, vxLoc, vyLoc, vzLoc, qhatLoc, enerLoc;
686  double betaLoc, gammaLoc, flowFactor;
687  int hydro_ctl;
688  double dt_lrf;
689 
690  double pc0[4] = {0.0};
691  double vc0[4] = {0.0};
692  double soln_alphas, prob_el;
693  double muD2;
694  double recordE0;
695 
696  // Convert hard parton momentum to onshell
697  pc0[1] = el_p0[1];
698  pc0[2] = el_p0[2];
699  pc0[3] = el_p0[3];
700  if (abs(pid) == 4 || abs(pid) == 5)
701  pc0[0] = sqrt(pc0[1] * pc0[1] + pc0[2] * pc0[2] + pc0[3] * pc0[3] +
702  HQ_mass * HQ_mass);
703  else
704  pc0[0] = sqrt(pc0[1] * pc0[1] + pc0[2] * pc0[2] + pc0[3] * pc0[3]);
705 
706  recordE0 = pc0[0];
707 
708  if (std::isinf(el_time) || std::isnan(el_time) || std::isinf(el_rz) ||
709  std::isnan(el_rz) || std::abs(el_rz) > el_time) {
710  JSWARN << "Second instance";
711  JSWARN << "el_vector for vector is:" << el_time << ", " << el_rx
712  << ", " << el_ry << ", " << el_rz;
713  JSWARN << "initR0, initRx, initRy, init Rz=" << initR0 << ", "
714  << initRx << ", " << initRy << ", " << initRz;
715  JSWARN << "initVx, initVy, initVz =" << initVx << ", " << initVy
716  << ", " << initVz;
717  JSWARN << "velocityMod=" << std::setprecision(20) << velocityMod;
718  JSWARN << "initVMod=" << std::setprecision(20)
719  << std::sqrt(initVx * initVx + initVy * initVy +
720  initVz * initVz);
721  Dump_pIn_info(i, pIn);
722  //exit(0);
723  }
724  GetHydroCellSignal(el_time, el_rx, el_ry, el_rz,
725  check_fluid_info_ptr);
726  VERBOSE(8) << MAGENTA << "Temperature from medium = "
727  << check_fluid_info_ptr->temperature;
728 
729  tempLoc = check_fluid_info_ptr->temperature;
730  sdLoc = check_fluid_info_ptr->entropy_density;
731  vxLoc = check_fluid_info_ptr->vx;
732  vyLoc = check_fluid_info_ptr->vy;
733  vzLoc = check_fluid_info_ptr->vz;
734 
735  vc0[1] = vxLoc;
736  vc0[2] = vyLoc;
737  vc0[3] = vzLoc;
738 
739  hydro_ctl = 0;
740 
741  if (hydro_ctl == 0 && tempLoc >= hydro_Tc) {
742 
743  trans(vc0, pc0);
744  enerLoc = pc0[0];
745  transback(vc0, pc0);
746 
747  betaLoc = sqrt(vxLoc * vxLoc + vyLoc * vyLoc + vzLoc * vzLoc);
748  gammaLoc = 1.0 / sqrt(1.0 - betaLoc * betaLoc);
749  flowFactor =
750  gammaLoc *
751  (1.0 - (initVx * vxLoc + initVy * vyLoc + initVz * vzLoc));
752 
753  //if(run_alphas==1){ alphas= 4*pi/(9.0*log(2*enerLoc*tempLoc/0.04));}
754 
755  /* if (qhat0 < 0.0) { // calculate qhat with alphas
756  muD2 = 6.0 * pi * alphas * tempLoc * tempLoc;
757  if (enerLoc > 2.0 * pi * tempLoc)
758  qhatLoc = Ca * 50.4864 / pi * pow(alphas, 2) * pow(tempLoc, 3) *
759  log(5.7 * enerLoc * tempLoc / 4.0 / muD2);
760  else
761  qhatLoc = Ca * 50.4864 / pi * pow(alphas, 2) * pow(tempLoc, 3) *
762  log(5.7 * 2.0 * pi * tempLoc * tempLoc / 4.0 / muD2);
763  if (qhatLoc < 0.0)
764  qhatLoc = 0.0;
765  } else { // use input qhat
766  if (brick_med)
767  qhatLoc = qhat0 * 0.1973;
768  else
769  qhatLoc = qhat0 / 96.0 * sdLoc * 0.1973;
770  }*/
771 
772  //GeneralQhatFunction(int QhatParametrizationType, double Temperature, double EntropyDensity, double FixAlphas, double Qhat0, double E, double muSquare)
773  double muSquare= pIn[i].t(); //Virtuality of the parent; Revist this when q-hat is virtuality dependent
774  qhatLoc= GeneralQhatFunction(QhatParametrizationType, tempLoc, sdLoc, alphas, qhat0, enerLoc, muSquare);
775 
776  } else { // outside the QGP medium
777  continue;
778  }
779 
780  // Calculating the delta t in the fluid rest frame
781  if (el_time + el_dt <= time)
782  dt_lrf = el_dt * flowFactor;
783  else
784  dt_lrf = (time - el_time) * flowFactor;
785 
786  //if(run_alphas==1){ alphas= 4*pi/(9.0*log(2*enerLoc*tempLoc/0.04));}
787 
788  // solve alphas
789  if (qhat0 < 0.0 || QhatParametrizationType==0 || QhatParametrizationType==1 || QhatParametrizationType==5 ||
790  QhatParametrizationType==6 || QhatParametrizationType==7 )
791  soln_alphas = alphas;
792  else
793  soln_alphas = solve_alphas(qhatLoc, enerLoc, tempLoc);
794 
795  // Calculate the proability of elastic scattering in time delta t in fluid rest frame
796  muD2 = 6.0 * pi * soln_alphas * tempLoc * tempLoc;
797  prob_el = 42.0 * zeta3 * el_CR * tempLoc / 6.0 / pi /
798  pi * dt_lrf / 0.1973;
799 
800  prob_el=prob_el*ModifiedProbability(QhatParametrizationType, tempLoc, sdLoc, enerLoc, pIn[i].t());
801 
802  el_rand = ZeroOneDistribution(*GetMt19937Generator());
803 
804  //cout << " qhat: " << qhatLoc << " alphas: " << soln_alphas << " ener: " << enerLoc << " prob_el: " << prob_el << " " << el_rand << endl;
805 
806  if (el_rand < prob_el) { // elastic scattering happens
807 
808  //cout << "elastic scattering happens" << endl;
809  int CT = -1;
810  int pid0 = -999;
811  int pid2 = -999;
812  int pid3 = -999;
813  double pc2[4] = {0.0}; // final recoil thermal parton
814  double pc3[4] = {0.0}; // initial thermal parton
815  double pc4[4] = {0.0}; // not used
816  double qt = 0.0; // not used
817  unsigned int el_max_color, el_color0, el_anti_color0, el_color2,
818  el_anti_color2, el_color3, el_anti_color3;
819 
820  el_max_color = pIn[i].max_color();
821  el_color0 = pIn[i].color();
822  el_anti_color0 = pIn[i].anti_color();
823 
824  pid0 = pid;
825 
826  // deterimine channel
827  //flavor(CT,pid0,pid2,pid3);
828  flavor(CT, pid0, pid2, pid3, el_max_color, el_color0,
829  el_anti_color0, el_color2, el_anti_color2, el_color3,
830  el_anti_color3);
831 
832  //cout << "color: " << el_color0 << " " << el_anti_color0 << " " << el_color2 << " " << el_anti_color2 << " " << el_color3 << " " << el_anti_color3 << " max: " << el_max_color << endl;
833 
834  // do scattering
835  if (CT == 11 || CT == 12) { // for heavy quark scattering
836  collHQ22(CT, tempLoc, muD2, vc0, pc0, pc2, pc3, pc4, qt);
837  } else { // for light parton scattering
838  colljet22(CT, tempLoc, muD2, vc0, pc0, pc2, pc3, pc4, qt);
839  }
840 
841  if (pc0[0] < pc2[0] && abs(pid0) != 4 &&
842  pid0 ==
843  pid2) { //disable switch for heavy quark, only allow switch for identical particles
844  double p0temp[4] = {0.0};
845  for (int k = 0; k <= 3; k++) {
846  p0temp[k] = pc2[k];
847  pc2[k] = pc0[k];
848  pc0[k] = p0temp[k];
849  }
850  }
851 
852  //Amit: do not add recoil parton with E=0,Px=0,Py=0,Pz=0 into pOut vector
853  if (pc2[0] < rounding_error || pc3[0] < rounding_error)
854  continue;
855 
856  // push out recoied and back-reaction (negative) partons
857  // need to add color information later!!!
858  double el_vertex[4];
859  el_vertex[0] = el_time;
860  el_vertex[1] = el_rx;
861  el_vertex[2] = el_ry;
862  el_vertex[3] = el_rz;
863 
864  int iout;
865  double ft;
866 
867  pc2[0] = sqrt(pc2[1] * pc2[1] + pc2[2] * pc2[2] + pc2[3] * pc2[3] +
868  rounding_error);
869 
870  if (std::isnan(pc2[1]) || std::isnan(pc2[2]) ||
871  std::isnan(pc2[3]) || std::isinf(pc2[1]) ||
872  std::isinf(pc2[2]) || std::isinf(pc2[3])) {
873  JSWARN << "recoil in MATTER instance 1: pc[0]=" << pc2[0]
874  << ", pc2[1]=" << pc2[1] << ", pc2[2]=" << pc2[2]
875  << ", pc2[3]=" << pc2[3];
876  }
877 
878  pOut.push_back(Parton(0, pid2, 1, pc2, el_vertex)); // recoiled
879  iout = pOut.size() - 1;
880  pOut[iout].set_jet_v(velocity_jet); // use initial jet velocity
881  pOut[iout].set_mean_form_time();
882  ft = 10000.0;
883  pOut[iout].set_form_time(ft);
891  pOut[iout].set_color(0);
892  pOut[iout].set_anti_color(0);
893  pOut[iout].set_max_color(pIn[i].max_color());
894  pOut[iout].set_min_color(pIn[i].min_color());
895  pOut[iout].set_min_anti_color(pIn[i].min_anti_color());
896 
897  pOut.push_back(
898  Parton(0, pid3, -1, pc3, el_vertex)); // back reaction
899  iout = pOut.size() - 1;
900  pOut[iout].set_jet_v(velocity_jet);
901  pOut[iout].set_mean_form_time();
902  ft = 10000.0;
903  pOut[iout].set_form_time(ft);
909  pOut[iout].set_color(0);
910  pOut[iout].set_anti_color(0);
911  pOut[iout].set_max_color(pIn[i].max_color());
912  pOut[iout].set_min_color(pIn[i].min_color());
913  pOut[iout].set_min_anti_color(pIn[i].min_anti_color());
914 
919 
920  } // end if scattering
921 
922  // E1'+E2 = E3'+E4 (all massless in scattering)
923  // Now I want E1+E2 = E3+E4 where 1 and 3 have mass
924  // -> E1-E1' = E3-E3' -> E3 = E1-E1'+E3'
925  // m3^2 = E3^2-E3'^2
926  for (int j = 1; j <= 3; j++)
927  el_p0[j] = pc0[j];
928 
929  el_p0[0] = el_p0[0] - recordE0 + pc0[0];
930  el_p0[4] = el_p0[0] * el_p0[0] - pc0[0] * pc0[0];
931  if (el_p0[4] < 0)
932  cout << "complain negative virt" << endl;
933 
934  } // end time loop for elastic scattering
935 
936  // do split
937  double t_used = pIn[i].t();
938  //if (t_used<QS) t_used = QS; // SC: not necessary
939  double tau_form = 2.0 * pIn[i].nu() / t_used;
940  double z_low = QS * QS / t_used / 2.0;
941  double z_hi = 1.0 - z_low;
942 
943  VERBOSE(8) << " zeta = " << zeta;
944 
945  if (pid == gid) { // gluon
946  double val1 =
947  P_z_gg_int(z_low, z_hi, zeta, t_used, tau_form, pIn[i].nu());
948  double val2 =
949  nf * P_z_qq_int(z_low, z_hi, zeta, t_used, tau_form, pIn[i].nu());
950  double M = PythiaFunction.particleData.m0(cid);
951  z_low = (QS * QS + 2.0 * M * M) / t_used / 2.0;
952  z_hi = 1.0 - z_low;
953  double val3 = 0.0;
954  if (z_hi > z_low)
955  val3 = P_z_qq_int_w_M_vac_only(M, z_low, z_hi, zeta, t_used,
956  tau_form, pIn[i].nu());
957  if (t_used < (QS * QS + 2.0 * M * M))
958  val3 = 0.0;
959 
960  M = PythiaFunction.particleData.m0(bid);
961  z_low = (QS * QS + 2.0 * M * M) / t_used / 2.0;
962  z_hi = 1.0 - z_low;
963  double val4 = 0.0;
964  if (z_hi > z_low)
965  val4 = P_z_qq_int_w_M_vac_only(M, z_low, z_hi, zeta, t_used,
966  tau_form, pIn[i].nu());
967  if (t_used < (QS * QS + 2.0 * M * M))
968  val4 = 0.0;
969 
970  VERBOSE(8) << BOLDYELLOW << " val1 = " << val1 << " val2 = " << val2
971  << " val3 = " << val3 << " val4 = " << val4;
972 
973  if (val1 < 0.0 || val2 < 0.0 || val3 < 0.0 || val4 < 0.0) {
974  cerr << " minus log of sudakov negative val1 , val2 , val3, val4 = "
975  << val1 << " " << val2 << " " << val3 << " " << val4
976  << endl;
977  throw std::runtime_error("minus log of sudakov negative");
978  // cin >> blurb ;
979  }
980 
981  double ratio1 = val1 / (val1 + val2 + val3 + val4);
982  double ratio2 = (val1 + val2) / (val1 + val2 + val3 + val4);
983  double ratio3 = (val1 + val2 + val3) / (val1 + val2 + val3 + val4);
984  double r = ZeroOneDistribution(*GetMt19937Generator());
985  if (r >= ratio1 && r < ratio2) { // qqbar
986 
987  double r2 = ZeroOneDistribution(*GetMt19937Generator());
988 
989  // assign flavors
990  if (r2 > 0.6666) {
991  pid_a = uid;
992  pid_b = -1 * uid;
993  } else if (r2 > 0.3333) {
994  pid_a = did;
995  pid_b = -1 * did;
996 
997  } else {
998  pid_a = sid;
999  pid_b = -1 * sid;
1000  }
1001  iSplit = 2;
1002  } else if (r >= ratio2 && r < ratio3) {
1003  pid_a = cid;
1004  pid_b = -cid;
1005  iSplit = 4;
1006 
1007  VERBOSE(1) << BOLDYELLOW << " split to c c-bar";
1008 
1009  } else if (r >= ratio3) {
1010  pid_a = bid;
1011  pid_b = -bid;
1012  iSplit = 5;
1013 
1014  VERBOSE(1) << BOLDYELLOW << " Split to b b-bar ";
1015  } else { // gg
1016  pid_a = gid;
1017  pid_b = gid;
1018  iSplit = 1;
1019  }
1020  } else if (std::abs(pIn[i].pid()) < 4) { // we had a light quark
1021 
1022  double relative_charge = 0.0;
1023  if ((std::abs(pIn[i].pid()) > 0) && (std::abs(pIn[i].pid()) < 4))
1024  relative_charge = 1.0 / 9.0;
1025  if (std::abs(pIn[i].pid()) == 2)
1026  relative_charge = relative_charge * 4.0;
1027  if (std::abs(pIn[i].pid()) == 4)
1028  relative_charge = 4.0 / 9.0;
1029  if (std::abs(pIn[i].pid()) == 5)
1030  relative_charge = 1.0 / 9.0;
1031 
1032  double ProbGluon =
1033  1.0 - sudakov_Pqg(QS * QS / 2, pIn[i].t(), zeta, pIn[i].nu());
1034  double ProbPhoton =
1035  1.0 -
1036  std::pow(sudakov_Pqp(QS * QS / 2, pIn[i].t(), zeta, pIn[i].nu()),
1037  relative_charge);
1038 
1039  double val = ProbGluon / (ProbGluon + ProbPhoton);
1040 
1041  VERBOSE(8) << MAGENTA
1042  << " probability of gluon radiation from quark = " << val;
1043 
1044  double r2 = ZeroOneDistribution(*GetMt19937Generator());
1045 
1046  if (r2 <= val) { // light quark decay to quark and gluon
1047  pid_a = pid;
1048  pid_b = gid;
1049  iSplit = 0; // iSplit for quark radiating a gluon
1050  } else { // light quark decay to quark and photon
1051  pid_a = pid;
1052  pid_b = photonid;
1053  iSplit = 3; // iSplit for quark radiating a photon
1054  photon_brem = true;
1055  }
1056 
1057  } else {
1058  // we had a heavy quark
1059  pid_a = pid;
1060  pid_b = gid;
1061  iSplit = 0;
1062  }
1063 
1064  // daughter virtualities
1065  double tQd1 = QS * QS;
1066  double tQd2 = QS * QS;
1067 
1068  double new_parent_p0 = el_p0[0];
1069  double new_parent_px = el_p0[1];
1070  double new_parent_py = el_p0[2];
1071  double new_parent_pz = el_p0[3];
1072  double new_parent_t = el_p0[4];
1073  double new_parent_pl = (new_parent_px * pIn[i].jet_v().x() +
1074  new_parent_py * pIn[i].jet_v().y() +
1075  new_parent_pz * pIn[i].jet_v().z()) /
1076  mod_jet_v;
1077  if (new_parent_pl < 0.0) {
1078  VERBOSE(8) << BOLDYELLOW
1079  << " parton traversing opposite to jet direction ... Just "
1080  "letting you know ! ";
1081  // cin >> blurb ;
1082  }
1083  //JSINFO << BOLDYELLOW << " old virtuality = " << pIn[i].t() << " new virtuality = " << new_parent_t ;
1084  //double new_parent_pl = sqrt(pow(new_parent_px,2)+pow(new_parent_py,2)+pow(new_parent_pz,2));
1085  //double new_parent_vx = new_parent_px/new_parent_pl;
1086  //double new_parent_vy = new_parent_py/new_parent_pl;
1087  //double new_parent_vz = new_parent_pz/new_parent_pl;
1088  double new_parent_nu = (new_parent_p0 + new_parent_pl) / sqrt(2.0);
1089 
1090  //set color of daughters here
1091  unsigned int d1_col, d1_acol, d2_col, d2_acol, color, anti_color;
1092  //std::uniform_int_distribution<short> uni(101,103);
1093  //color = pIn[i].color();
1094 
1095  if (iSplit != 3) // not photon radiation, generate new colors
1096  {
1097  max_color = pIn[i].max_color();
1098  //if (pIn[i].anti_color()>maxcolor) color = pIn[i].anti_color();
1099  JSDEBUG << " old max color = " << max_color;
1100  max_color = ++MaxColor;
1101  color = max_color;
1102  anti_color = max_color;
1103  pIn[i].set_max_color(max_color);
1104  JSDEBUG << " new color = " << color;
1105  }
1106 
1107  if (iSplit == 1)
1108  {
1109  d1_col = pIn[i].color();
1110  d2_col = color;
1111  d1_acol = anti_color;
1112  d2_acol = pIn[i].anti_color();
1113  } else if (
1114  iSplit ==
1115  0)
1116  {
1117  if (pIn[i].pid() > 0) // parent is a quark
1118  {
1119  d1_col = color;
1120  d1_acol = 0;
1121  d2_col = pIn[i].color();
1122  d2_acol = anti_color;
1123  } else {
1124  d1_col = 0; // parent is an anti-quark
1125  d1_acol = anti_color;
1126  d2_col = color;
1127  d2_acol = pIn[i].anti_color();
1128  }
1129  } else if (
1130  iSplit == 2 || iSplit == 4 ||
1131  iSplit ==
1132  5) // gluon splits into quark anti-quark, c anti-c, b anti-b
1133  {
1134  d1_col = pIn[i].color();
1135  d1_acol = 0;
1136  d2_acol = pIn[i].anti_color();
1137  d2_col = 0;
1138  } else if (
1139  iSplit ==
1140  3) // radiating a photon has col = acol = 0, all color remains in quark(anti-quark)
1141  {
1142  d1_col = pIn[i].color();
1143  d1_acol = pIn[i].anti_color();
1144  d2_col = 0;
1145  d2_acol = 0;
1146  } else {
1147  throw std::runtime_error("error in iSplit");
1148  }
1149 
1150  double l_perp2 = -1.0; // SC: initialization
1151  int ifcounter = 0;
1152 
1153  while ((l_perp2 <= Lambda_QCD * Lambda_QCD) && (ifcounter < 100)) {
1154 
1155  // if(abs(pid) == 4 || abs(pid) == 5)
1156  // {
1157  z = generate_vac_z_w_M(pid, pIn[i].restmass(), QS * QS / 2.0,
1158  pIn[i].t(), zeta, pIn[i].nu(), iSplit);
1159  // }
1160  // else if (pid == 21 && ( iSplit == 4 || iSplit == 5 ) )
1161  // {
1162  // z = generate_vac_z(pid, QS/2.0, pIn[i].t(), zeta, pIn[i].nu(), iSplit);
1163  // }
1164  VERBOSE(8) << MAGENTA << " generated z = " << z;
1165 
1166  int iSplit_a = 0;
1167  if (pid_a == gid)
1168  iSplit_a = 1;
1169 
1170  // use pIn information to sample z above, but use new_parent to calculate daughter partons below
1171 
1172  if (std::abs(pid_a) == 4 || std::abs(pid_a) == 5) {
1173  double M = PythiaFunction.particleData.m0(pid_a);
1174  if (QS * QS * (1.0 + std::sqrt(1.0 + 4.0 * M * M / QS / QS)) / 2.0 <
1175  z * z * new_parent_t) {
1176  tQd1 = generate_vac_t_w_M(
1177  pid_a, pIn[i].restmass(), z * new_parent_nu, QS * QS / 2.0,
1178  z * z * new_parent_t,
1179  zeta + std::sqrt(2) * pIn[i].form_time() * fmToGeVinv,
1180  iSplit_a);
1181 
1182  } else {
1183  tQd1 = z * z * new_parent_t;
1184  }
1185  } else if (pid_a == 21) {
1186  double M = 0.0;
1187  if ((QS * QS) < z * z * new_parent_t) {
1188  tQd1 = generate_vac_t_w_M(
1189  pid_a, M, z * new_parent_nu, QS * QS / 2.0,
1190  z * z * new_parent_t,
1191  zeta + std::sqrt(2) * pIn[i].form_time() * fmToGeVinv,
1192  iSplit_a);
1193  } else {
1194  tQd1 = z * z * new_parent_t;
1195  }
1196  VERBOSE(8) << BOLDYELLOW << " daughter 1 virt generated = " << tQd1;
1197 
1198  } else // light quark anti-quark
1199  {
1200  if (z * z * new_parent_t > QS * QS) {
1201  tQd1 = generate_vac_t(
1202  pid_a, z * new_parent_nu, QS * QS / 2.0, z * z * new_parent_t,
1203  zeta + std::sqrt(2) * pIn[i].form_time() * fmToGeVinv,
1204  iSplit_a);
1205  } else { // SC
1206  tQd1 = z * z * new_parent_t;
1207  }
1208  }
1209 
1210  int iSplit_b = 0;
1211  if (pid_b == gid)
1212  iSplit_b = 1;
1213 
1214  if (std::abs(pid_b) == 4 || std::abs(pid_b) == 5) {
1215  double M = PythiaFunction.particleData.m0(pid_b);
1216 
1217  if (QS * QS * (1.0 + std::sqrt(1.0 + 4.0 * M * M / QS / QS)) / 2.0 <
1218  (1.0 - z) * (1.0 - z) * new_parent_t) {
1219  tQd2 = generate_vac_t_w_M(
1220  pid_b, pIn[i].restmass(), (1.0 - z) * new_parent_nu,
1221  QS * QS / 2.0, (1.0 - z) * (1.0 - z) * new_parent_t,
1222  zeta + std::sqrt(2) * pIn[i].form_time() * fmToGeVinv,
1223  iSplit_b);
1224 
1225  } else {
1226  tQd2 = (1.0 - z) * (1.0 - z) * new_parent_t;
1227  }
1228 
1229  } else if (pid_b == 21) {
1230  double M = 0.0;
1231  if (QS * QS < (1.0 - z) * (1.0 - z) * new_parent_t) {
1232  tQd2 = generate_vac_t_w_M(
1233  pid_b, M, (1.0 - z) * new_parent_nu, QS * QS / 2.0,
1234  (1.0 - z) * (1.0 - z) * new_parent_t,
1235  zeta + std::sqrt(2) * pIn[i].form_time() * fmToGeVinv,
1236  iSplit_b);
1237 
1238  } else {
1239  tQd2 = (1.0 - z) * (1.0 - z) * new_parent_t;
1240  }
1241  VERBOSE(8) << BOLDYELLOW << " daughter virt generated = " << tQd2;
1242 
1243  } else {
1244  if (((1.0 - z) * (1.0 - z) * new_parent_t > QS * QS) &&
1245  (iSplit != 3)) {
1246  tQd2 = generate_vac_t(
1247  pid_b, (1.0 - z) * new_parent_nu, QS * QS / 2.0,
1248  (1.0 - z) * (1.0 - z) * new_parent_t,
1249  zeta + std::sqrt(2) * pIn[i].form_time() * fmToGeVinv,
1250  iSplit_b);
1251  } else { // SC
1252  tQd2 = (1.0 - z) * (1.0 - z) * new_parent_t;
1253  }
1254 
1255  if (iSplit == 3) {
1256  tQd2 = rounding_error; // forcing the photon to have no virtuality
1257  }
1258  }
1259 
1260  l_perp2 = new_parent_t * z * (1.0 - z) - tQd2 * z -
1261  tQd1 * (1.0 - z);
1262 
1263  if (abs(pid) == 4 || abs(pid) == 5) {
1264  l_perp2 = new_parent_t * z * (1.0 - z) - tQd2 * z -
1265  tQd1 * (1.0 - z) -
1266  std::pow((1.0 - z) * pIn[i].restmass(),
1267  2);
1268  } else if ((pid == 21) &&
1269  (iSplit > 3)) // gluon decay into heavy quark anti-quark
1270  {
1271 
1272  double M = PythiaFunction.particleData.m0(pid_a);
1273 
1274  l_perp2 = new_parent_t * z * (1.0 - z) - tQd2 * z -
1275  tQd1 * (1.0 - z) - M * M;
1276  } else {
1277  l_perp2 = new_parent_t * z * (1.0 - z) - tQd2 * z -
1278  tQd1 * (1.0 - z);
1279  }
1280 
1281  ifcounter++;
1282  }
1283 
1284  // std::ofstream zdist;
1285  // zdist.open("zdist_heavy.dat", std::ios::app);
1286  // std::ofstream tdist;
1287  // tdist.open("tdist_heavy.dat", std::ios::app);
1288  // if (std::abs(pid_a) == 4 || std::abs(pid_a) == 5)
1289  // {
1290  // tdist << tQd1<< endl;
1291  // zdist << z << endl;
1292  // }
1293  // else if (std::abs(pid_b) == 4 || std::abs(pid_b) == 5)
1294  // {
1295  // tdist << tQd2<< endl;
1296  // zdist << z << endl;
1297  // }
1298  // tdist.close();
1299  // zdist.close();
1300 
1301  if (l_perp2 <= Lambda_QCD * Lambda_QCD)
1302  l_perp2 = Lambda_QCD * Lambda_QCD;
1303  double l_perp = std::sqrt(
1304  l_perp2);
1305  VERBOSE(8) << BOLDYELLOW << " after ifcounter = " << ifcounter
1306  << " l_perp2 = " << l_perp2;
1307  VERBOSE(8) << BOLDYELLOW << " z = " << z << " tQd1 = " << tQd1
1308  << " tQd2 = " << tQd2;
1309  VERBOSE(8) << BOLDYELLOW << " pid_a = " << pid_a
1310  << " pid_b = " << pid_b;
1311 
1312  // axis of split
1313  double angle = generate_angle();
1314 
1315  // KK: changed to x,y,z
1316  //double parent_perp = std::sqrt( pow(pIn[i].px(),2) + pow(pIn[i].py(),2) + pow(pIn[i].pz(),2) - pow(pIn[i].pl(),2) );
1317  //double mod_jet_v = std::sqrt( pow(pIn[i].jet_v().x(),2) + pow(pIn[i].jet_v().y(),2) + pow(pIn[i].jet_v().z(),2) ) ;
1318  double c_t =
1319  pIn[i].jet_v().z() /
1320  mod_jet_v; // c_t is cos(theta) for the jet_velocity unit vector
1321  double s_t = std::sqrt(1.0 - c_t * c_t); // s_t is sin(theta)
1322 
1323  double s_p = pIn[i].jet_v().y() / std::sqrt(pow(pIn[i].jet_v().x(), 2) +
1324  pow(pIn[i].jet_v().y(), 2));
1325  double c_p = pIn[i].jet_v().x() / std::sqrt(pow(pIn[i].jet_v().x(), 2) +
1326  pow(pIn[i].jet_v().y(), 2));
1327 
1328  VERBOSE(8) << BOLDYELLOW
1329  << " Jet direction w.r.t. beam: theta = " << std::acos(c_t)
1330  << " phi = " << std::acos(c_p);
1331 
1332  //double c_t = new_parent_vz;
1333  //double s_t = std::sqrt( 1.0 - c_t*c_t) ;
1334  //
1335  //double s_p = new_parent_vy/std::sqrt( pow(new_parent_vx,2) + pow(new_parent_vy,2) ) ;
1336  //double c_p = new_parent_vx/std::sqrt( pow(new_parent_vx,2) + pow(new_parent_vy,2) ) ;
1337 
1338  // First daughter
1339  double k_perp1[4];
1340  k_perp1[0] = 0.0;
1341  k_perp1[1] = z * (new_parent_px - new_parent_pl * s_t * c_p) +
1342  l_perp * std::cos(angle) * c_t * c_p -
1343  l_perp * std::sin(angle) * s_p;
1344  k_perp1[2] = z * (new_parent_py - new_parent_pl * s_t * s_p) +
1345  l_perp * std::cos(angle) * c_t * s_p +
1346  l_perp * std::sin(angle) * c_p;
1347  k_perp1[3] = z * (new_parent_pz - new_parent_pl * c_t) -
1348  l_perp * std::cos(angle) * s_t;
1349  double k_perp1_2 =
1350  pow(k_perp1[1], 2) + pow(k_perp1[2], 2) + pow(k_perp1[3], 2);
1351 
1352  double M = 0.0;
1353 
1354  if ((std::abs(pid_a) == 4) || (std::abs(pid_a) == 5))
1355  M = PythiaFunction.particleData.m0(pid_a);
1356 
1357  double energy = (z * new_parent_nu + (tQd1 + k_perp1_2 + M * M) /
1358  (2.0 * z * new_parent_nu)) /
1359  std::sqrt(2.0);
1360  double plong = (z * new_parent_nu - (tQd1 + k_perp1_2 + M * M) /
1361  (2.0 * z * new_parent_nu)) /
1362  std::sqrt(2.0);
1363  if (energy < 0.0) {
1364  JSWARN << " Energy negative after rotation, press 1 and return to "
1365  "continue ";
1366  cin >> blurb;
1367  }
1368 
1369  double newp[4];
1370  newp[0] = energy;
1371  newp[1] = plong * s_t * c_p + k_perp1[1];
1372  newp[2] = plong * s_t * s_p + k_perp1[2];
1373  newp[3] = plong * c_t + k_perp1[3];
1374 
1375  VERBOSE(8) << MAGENTA << " D1 px = " << newp[1] << " py = " << newp[2]
1376  << " pz = " << newp[3] << " E = " << newp[0];
1377  double newx[4];
1378  //newx[0] = time + deltaT;
1379  newx[0] = time;
1380  for (int j = 1; j <= 3; j++) {
1381  //newx[j] = pIn[i].x_in().comp(j) + (time + deltaT - pIn[i].x_in().comp(0))*velocity[j]/velocityMod;
1382  newx[j] = pIn[i].x_in().comp(j) +
1383  (time - pIn[i].x_in().comp(0)) * velocity[j];
1384  }
1385 
1386  VERBOSE(8) << BOLDRED << " PiD - a = " << pid_a;
1387 
1388  pOut.push_back(Parton(0, pid_a, jet_stat, newp, newx));
1389  int iout = pOut.size() - 1;
1390 
1391  if (std::isnan(newp[1]) || std::isnan(newp[2]) || std::isnan(newp[3])) {
1392  JSINFO << MAGENTA << plong << " " << s_t << " " << c_p << " "
1393  << k_perp1[1];
1394 
1395  JSINFO << MAGENTA << newp[0] << " " << newp[1] << " " << newp[2]
1396  << " " << newp[3];
1397  cin >> blurb;
1398  }
1399  pOut[iout].set_jet_v(velocity_jet); // use initial jet velocity
1400  pOut[iout].set_mean_form_time();
1401  double ft = generate_L(pOut[iout].mean_form_time());
1402  pOut[iout].set_form_time(ft);
1403  pOut[iout].set_color(d1_col);
1404  pOut[iout].set_anti_color(d1_acol);
1405  pOut[iout].set_max_color(max_color);
1406  pOut[iout].set_min_color(pIn[i].min_color());
1407  pOut[iout].set_min_anti_color(pIn[i].min_anti_color());
1408 
1409  VERBOSE(8) << BOLDRED << " virtuality of D 1 = " << pOut[iout].t();
1410  VERBOSE(8) << BOLDRED << " mass of parton = " << pOut[iout].restmass();
1411 
1412  // Second daughter
1413  double k_perp2[4];
1414  k_perp2[0] = 0.0;
1415  k_perp2[1] = (1.0 - z) * (new_parent_px - new_parent_pl * s_t * c_p) -
1416  l_perp * std::cos(angle) * c_t * c_p +
1417  l_perp * std::sin(angle) * s_p;
1418  k_perp2[2] = (1.0 - z) * (new_parent_py - new_parent_pl * s_t * s_p) -
1419  l_perp * std::cos(angle) * c_t * s_p -
1420  l_perp * std::sin(angle) * c_p;
1421  k_perp2[3] = (1.0 - z) * (new_parent_pz - new_parent_pl * c_t) +
1422  l_perp * std::cos(angle) * s_t;
1423  double k_perp2_2 =
1424  pow(k_perp2[1], 2) + pow(k_perp2[2], 2) + pow(k_perp2[3], 2);
1425 
1426  M = 0.0;
1427  if ((std::abs(pid_b) == 4) || (std::abs(pid_b) == 5))
1428  M = PythiaFunction.particleData.m0(pid_b);
1429  ;
1430 
1431  energy =
1432  ((1.0 - z) * new_parent_nu +
1433  (tQd2 + k_perp2_2 + M * M) / (2.0 * (1.0 - z) * new_parent_nu)) /
1434  std::sqrt(2.0);
1435  plong =
1436  ((1.0 - z) * new_parent_nu -
1437  (tQd2 + k_perp2_2 + M * M) / (2.0 * (1.0 - z) * new_parent_nu)) /
1438  std::sqrt(2.0);
1439 
1440  //parent_perp = std::sqrt( pow(pIn[i].p(1),2) + pow(pIn[i].p(2),2) + pow(pIn[i].p(3),2) - pow(pIn[i].pl(),2) );
1441  //mod_jet_v = std::sqrt( pow(pIn[i].jet_v().x(),2) + pow(pIn[i].jet_v().y(),2) + pow(pIn[i].jet_v().z(),2) ) ;
1442 
1443  if (energy < 0.0) {
1444  JSWARN << " Energy of 2nd daughter negative after rotation, press 1 "
1445  "and return to continue ";
1446  cin >> blurb;
1447  }
1448 
1449  newp[0] = energy;
1450  newp[1] = plong * s_t * c_p + k_perp2[1];
1451  newp[2] = plong * s_t * s_p + k_perp2[2];
1452  newp[3] = plong * c_t + k_perp2[3];
1453 
1454  if (std::isnan(newp[1]) || std::isnan(newp[2]) || std::isnan(newp[3])) {
1455  JSINFO << MAGENTA << "THIS IS THE SECOND DAUGHTER";
1456  JSINFO << MAGENTA << plong << " " << s_t << " " << c_p << " "
1457  << k_perp1[1];
1458  JSINFO << MAGENTA << newp[0] << " " << newp[1] << " " << newp[2]
1459  << " " << newp[3];
1460  cin >> blurb;
1461  }
1462 
1463  VERBOSE(8) << MAGENTA << " D1 px = " << newp[1] << " py = " << newp[2]
1464  << " pz = " << newp[3] << " E = " << newp[0];
1465 
1466  //newx[0] = time + deltaT;
1467  newx[0] = time;
1468  for (int j = 1; j <= 3; j++) {
1469  //newx[j] = pIn[i].x_in().comp(j) + (time + deltaT - pIn[i].x_in().comp(0))*velocity[j]/velocityMod;
1470  newx[j] = pIn[i].x_in().comp(j) +
1471  (time - pIn[i].x_in().comp(0)) * velocity[j];
1472  }
1473 
1474  if (iSplit != 3) // not a photon
1475  {
1476 
1477  VERBOSE(8) << BOLDRED << " PiD - b = " << pid_b;
1478  pOut.push_back(Parton(0, pid_b, jet_stat, newp, newx));
1479  iout = pOut.size() - 1;
1480  pOut[iout].set_jet_v(velocity_jet); // use initial jet velocity
1481  pOut[iout].set_mean_form_time();
1482  ft = generate_L(pOut[iout].mean_form_time());
1483  pOut[iout].set_form_time(ft);
1484  pOut[iout].set_color(d2_col);
1485  pOut[iout].set_anti_color(d2_acol);
1486  pOut[iout].set_max_color(max_color);
1487  pOut[iout].set_min_color(pIn[i].min_color());
1488  pOut[iout].set_min_anti_color(pIn[i].min_anti_color());
1489 
1490  } else // is a photon
1491  {
1492  VERBOSE(8) << BOLDRED << " is a photon PiD - b = " << pid_b;
1493  pOut.push_back(Photon(0, pid_b, jet_stat, newp, newx));
1494  iout = pOut.size() - 1;
1495  pOut[iout].set_jet_v(velocity_jet); // use initial jet velocity
1496  pOut[iout].set_mean_form_time();
1497  ft = generate_L(pOut[iout].mean_form_time());
1498  pOut[iout].set_form_time(1.0 / rounding_error);
1499  pOut[iout].set_color(0);
1500  pOut[iout].set_anti_color(0);
1501  pOut[iout].set_max_color(max_color);
1502  pOut[iout].set_min_color(pIn[i].min_color());
1503  pOut[iout].set_min_anti_color(pIn[i].min_anti_color());
1504 
1505  //JSINFO << BOLDYELLOW << " A photon was made in MATTER with px = " << pOut[iout].px() << " and sent to the framework " ;
1506  }
1507 
1508  VERBOSE(8) << BOLDRED << " virtuality of D 2 = " << pOut[iout].t();
1509 
1510  } else { // not time to split yet broadening it
1511 
1512  if (broadening_on) {
1513 
1514  double now_zeta =
1515  ((time + initRdotV + (time - initR0)) / std::sqrt(2)) *
1516  fmToGeVinv;
1517  qhat = fncQhat(now_zeta);
1518  if (now_temp > 0.1) {
1519  ehat = 0.0 * qhat / 4.0 / now_temp;
1520  } else {
1521  ehat = 0.0 * qhat / 4.0 / 0.3;
1522  }
1523 
1524  VERBOSE(8) << BOLDRED << " between splits broadening qhat = " << qhat
1525  << " ehat = " << ehat << " and delT = " << delT;
1526  VERBOSE(8) << BOLDBLUE << " zeta at formation = " << zeta
1527  << " zeta now = " << now_zeta;
1528 
1529  if ((!recoil_on) && (qhat > 0.0)) {
1530  double kt = 0;
1531  if (pIn[i].pid() == 21) // particle is a gluon
1532  {
1533  kt = generate_kt(qhat * 1.414 / 0.197, delT);
1534  } else if ((pIn[i].pid() < 6) &&
1535  (pIn[i].pid() > -6)) // particle is a quark
1536  {
1537  kt = generate_kt(qhat * 1.414 / 0.197 * Cf / Ca,
1538  delT); // scale down q-hat by Cf/Ca
1539  } else // a photon, or something else that does not have color
1540  {
1541  kt = 0;
1542  }
1543 
1544  JSDEBUG << " kt generated = " << kt
1545  << " for qhat = " << qhat * 1.414 / 0.197
1546  << " and delT = " << delT;
1547 
1548  double ktx, kty, ktz;
1549  ktx = kty = ktz = 0.0;
1550  double vx = initVx;
1551  double vy = initVy;
1552  double vz = initVz;
1553 
1554  bool solved = false;
1555  int trials = 0;
1556  while ((!solved) && (trials < 1000)) {
1557 
1558  if ((abs(vy) > approx) || (abs(vz) > approx)) {
1559  ktx =
1560  kt * (1 - 2 * ZeroOneDistribution(*GetMt19937Generator()));
1561 
1562  JSDEBUG << " vx = " << vx << " vy = " << vy << " vz = " << vz;
1563 
1564  double rad = sqrt(
1565  4 * ktx * ktx * vx * vx * vy * vy -
1566  4 * (vy * vy + vz * vz) *
1567  (ktx * ktx * (vx * vx + vz * vz) - kt * kt * vz * vz));
1568 
1569  double sol1 =
1570  (-2 * ktx * vx * vy + rad) / 2 / (vy * vy + vz * vz);
1571  double sol2 =
1572  (-2 * ktx * vx * vy - rad) / 2 / (vy * vy + vz * vz);
1573 
1574  kty = sol1;
1575 
1576  if ((ktx * ktx + sol1 * sol1) > kt * kt)
1577  kty = sol2;
1578 
1579  if ((ktx * ktx + kty * kty) < kt * kt) {
1580  ktz = sqrt(kt * kt - ktx * ktx - kty * kty);
1581 
1582  double sign = ZeroOneDistribution(*GetMt19937Generator());
1583 
1584  if (sign > 0.5)
1585  ktz = -1 * ktz;
1586 
1587  if (vz != 0)
1588  ktz = (-1 * ktx * vx - kty * vy) / vz;
1589 
1590  solved = true;
1591  }
1592 
1593  } else {
1594  ktx = 0;
1595  kty =
1596  kt * (1 - 2 * ZeroOneDistribution(*GetMt19937Generator()));
1597  double sign = ZeroOneDistribution(*GetMt19937Generator());
1598  if (sign > 0.5)
1599  kty = -1 * kty;
1600  ktz = sqrt(kt * kt - kty * kty);
1601  sign = ZeroOneDistribution(*GetMt19937Generator());
1602  if (sign > 0.5)
1603  ktz = -1 * ktz;
1604  }
1605  trials++;
1606  }
1607 
1608  VERBOSE(8) << MAGENTA << " ktx = " << ktx << " kty = " << kty
1609  << " ktz = " << ktz;
1610 
1611  double px = pIn[i].px();
1612  double py = pIn[i].py();
1613  double pz = pIn[i].pz();
1614  double energy = pIn[i].e();
1615 
1616  double p = sqrt(px * px + py * py + pz * pz);
1617 
1618  VERBOSE(8) << BOLDBLUE << " p before b & d, E = " << energy
1619  << " pz = " << pz << " px = " << px << " py = " << py;
1620 
1621  px += ktx;
1622  py += kty;
1623  pz += ktz;
1624 
1625  double np = sqrt(px * px + py * py + pz * pz);
1626  JSDEBUG << " p = " << p << " np = " << np;
1627 
1628  double correction = 0.0;
1629 
1630  if (np * np > p * p)
1631  correction = sqrt(np * np - p * p);
1632 
1633  px -= correction * vx;
1634  py -= correction * vy;
1635  pz -= correction * vz;
1636 
1637  // double nnp = sqrt(px*px + py*py + pz*pz);
1638 
1639  /* if (abs(nnp-p)>abs(np-p))
1640  {
1641 
1642  VERBOSE(8) << MAGENTA << " negative condition invoked ! " ;
1643 
1644  px += 2*(np-p)*vx;
1645  py += 2*(np-p)*vy;
1646  pz += 2*(np-p)*vz;
1647  }
1648 */
1649  np = sqrt(px * px + py * py + pz * pz);
1650  JSDEBUG << "FINAL p = " << p << " np = " << np;
1651 
1652  pOut.push_back(pIn[i]);
1653  int iout = pOut.size() - 1;
1654 
1655  pOut[iout].reset_p(px, py, pz);
1656 
1657  double drag = ehat * delT;
1658 
1659  VERBOSE(8) << MAGENTA << " drag = " << drag
1660  << " temperature = " << now_temp;
1661 
1662  if ((np > drag) && (energy > drag) &&
1663  (energy > sqrt(px * px + py * py + pz * pz))) {
1664  px -= drag * vx;
1665  py -= drag * vy;
1666  pz -= drag * vz;
1667  energy -= drag;
1668  pOut[iout].reset_momentum(px, py, pz, energy);
1669  }
1670 
1671  VERBOSE(8) << BOLDYELLOW << " p after b & d, E = " << energy
1672  << " pz = " << pz << " px = " << px << " py = " << py;
1673  }
1674 
1675  } // end if(broadening_on)
1676  pOut.push_back(pIn[i]);
1677  }
1678  } else { // virtuality too low lets broaden it
1679 
1680  if (broadening_on) {
1681 
1682  double now_zeta =
1683  ((time + initRdotV + (time - initR0)) / std::sqrt(2)) * fmToGeVinv;
1684  //double now_zeta = ( ( time + initRdotV )/std::sqrt(2) )*fmToGeVinv;
1685  qhat = fncQhat(now_zeta);
1686  if (now_temp > 0.1) {
1687  ehat = 0.0 * qhat / 4.0 / now_temp;
1688  } else {
1689  ehat = 0.0 * qhat / 4.0 / 0.3;
1690  }
1691 
1692  VERBOSE(8) << BOLDRED << " after splits broadening qhat = " << qhat
1693  << " ehat = " << ehat << " and delT = " << delT;
1694  VERBOSE(8) << BOLDBLUE << " zeta at formation = " << zeta
1695  << " zeta now = " << now_zeta;
1696 
1697  //JSINFO << " broadening qhat = " << qhat << " and delT = " << delT ;
1698 
1699  if ((!recoil_on) && (qhat > 0.0)) {
1700  double kt = 0.0;
1701 
1702  if (pIn[i].pid() == 21) {
1703  kt = generate_kt(qhat * 1.414 / 0.197, delT);
1704  } else {
1705  kt = generate_kt(qhat * 1.414 / 0.197 * Cf / Ca, delT);
1706  }
1707 
1708  JSDEBUG << " kt generated = " << kt
1709  << " for qhat = " << qhat * 1.414 / 0.197
1710  << " and delT = " << delT;
1711 
1712  double ktx, kty, ktz;
1713  ktx = kty = ktz = 0;
1714  double vx = initVx;
1715  double vy = initVy;
1716  double vz = initVz;
1717 
1718  bool solved = false;
1719 
1720  int trials = 0;
1721  while ((!solved) && (trials < 1000)) {
1722  if ((abs(vy) > approx) || (abs(vz) > approx)) {
1723 
1724  ktx = kt * (1 - 2 * ZeroOneDistribution(*GetMt19937Generator()));
1725 
1726  JSDEBUG << " vx = " << vx << " vy = " << vy << " vz = " << vz;
1727 
1728  double rad = sqrt(
1729  4 * ktx * ktx * vx * vx * vy * vy -
1730  4 * (vy * vy + vz * vz) *
1731  (ktx * ktx * (vx * vx + vz * vz) - kt * kt * vz * vz));
1732 
1733  double sol1 =
1734  (-2 * ktx * vx * vy + rad) / 2 / (vy * vy + vz * vz);
1735  double sol2 =
1736  (-2 * ktx * vx * vy - rad) / 2 / (vy * vy + vz * vz);
1737 
1738  kty = sol1;
1739 
1740  if ((ktx * ktx + sol1 * sol1) > kt * kt)
1741  kty = sol2;
1742 
1743  if ((ktx * ktx + kty * kty) < kt * kt) {
1744  ktz = sqrt(kt * kt - ktx * ktx - kty * kty);
1745 
1746  double sign = ZeroOneDistribution(*GetMt19937Generator());
1747  if (sign > 0.5)
1748  ktz = -1 * ktz;
1749 
1750  if (vz != 0)
1751  ktz = (-1 * ktx * vx - kty * vy) / vz;
1752 
1753  solved = true;
1754  }
1755  } else {
1756  ktx = 0;
1757  kty = kt * (1 - 2 * ZeroOneDistribution(*GetMt19937Generator()));
1758  double sign = ZeroOneDistribution(*GetMt19937Generator());
1759  if (sign > 0.5)
1760  kty = -1 * kty;
1761  ktz = sqrt(kt * kt - kty * kty);
1762  sign = ZeroOneDistribution(*GetMt19937Generator());
1763  if (sign > 0.5)
1764  ktz = -1 * ktz;
1765  }
1766  trials++;
1767  }
1768 
1769  VERBOSE(8) << MAGENTA << " ktx = " << ktx << " kty = " << kty
1770  << " ktz = " << ktz;
1771 
1772  double px = pIn[i].px();
1773  double py = pIn[i].py();
1774  double pz = pIn[i].pz();
1775  double energy = pIn[i].e();
1776 
1777  double p = sqrt(px * px + py * py + pz * pz);
1778 
1779  VERBOSE(8) << BOLDBLUE << " p before b & d, E = " << energy
1780  << " pz = " << pz << " px = " << px << " py = " << py;
1781 
1782  px += ktx;
1783  py += kty;
1784  pz += ktz;
1785 
1786  double np = sqrt(px * px + py * py + pz * pz);
1787  JSDEBUG << " p = " << p << " np = " << np;
1788  double correction = 0.0;
1789  if (np * np > p * p)
1790  correction = sqrt(np * np - p * p);
1791 
1792  px -= correction * vx;
1793  py -= correction * vy;
1794  pz -= correction * vz;
1795 
1796  // double nnp = sqrt(px*px + py*py + pz*pz);
1797 
1798  /* if (abs(nnp-p)>abs(np-p))
1799  {
1800  px += 2*(np-p)*vx;
1801  py += 2*(np-p)*vy;
1802  pz += 2*(np-p)*vz;
1803  }
1804 */
1805  np = sqrt(px * px + py * py + pz * pz);
1806  JSDEBUG << "FINAL p = " << p << " nnnp = " << np;
1807 
1808  pOut.push_back(pIn[i]);
1809  int iout = pOut.size() - 1;
1810 
1811  pOut[iout].reset_p(px, py, pz);
1812  np = sqrt(px * px + py * py + pz * pz);
1813  double drag = ehat * delT;
1814 
1815  VERBOSE(8) << MAGENTA << " drag = " << drag
1816  << " temperature = " << now_temp;
1817 
1818  if ((np > drag) && (energy > drag) &&
1819  (energy > sqrt(px * px + py * py + pz * pz))) {
1820  px -= drag * vx;
1821  py -= drag * vy;
1822  pz -= drag * vz;
1823  energy -= drag;
1824  pOut[iout].reset_momentum(px, py, pz, energy);
1825  }
1826 
1827  VERBOSE(8) << BOLDYELLOW << " p after b & d, E = " << energy
1828  << " pz = " << pz << " px = " << px << " py = " << py;
1829  }
1830 
1831  //pOut.push_back(pIn[i]);
1832  }
1833  }
1834 
1835  } // particle loop
1836 }
1837 
1838 double Matter::generate_kt(double local_qhat, double dzeta) {
1839  double width = local_qhat * dzeta;
1840 
1841  double x, r;
1842 
1843  r = ZeroOneDistribution(*GetMt19937Generator());
1844  x = -log(1.0 - r);
1845  if (x < 0.0)
1846  throw std::runtime_error(" k_t^2 < 0.0 ");
1847  double kt = sqrt(x * local_qhat * dzeta);
1848 
1849  return (kt);
1850 }
1851 
1852 double Matter::generate_angle() {
1853  double ang, r, blurb;
1854 
1855  // r = double(random())/ (maxN );
1856  r = ZeroOneDistribution(*GetMt19937Generator());
1857  // r = mtrand1();
1858 
1859  // cout << " r = " << r << endl;
1860  // cin >> blurb ;
1861 
1862  ang = r * 2.0 * pi;
1863 
1864  return (ang);
1865 }
1866 
1867 double Matter::generate_vac_t(int p_id, double nu, double t0, double t,
1868  double loc_a, int is) {
1869  double r, z, ratio, diff, scale, t_low, t_hi, t_mid, numer, denom, test;
1870 
1871  // r = double(random())/ (maxN );
1872  r = ZeroOneDistribution(*GetMt19937Generator());
1873  // r = mtrand1();
1874 
1875  if ((r >= 1.0) || (r <= 0.0)) {
1876  throw std::runtime_error(
1877  "error in random number in t *GetMt19937Generator()");
1878  }
1879 
1880  ratio = 1.0;
1881 
1882  diff = (ratio - r) / r;
1883 
1884  t_low = 2.0 * t0;
1885 
1886  t_hi = t;
1887 
1888  tscale = t; //for virtuality dependent q-hat
1889  // cout << " in gen_vac_t : t_low , t_hi = " << t_low << " " << t_hi << endl;
1890  // cin >> test ;
1891 
1892  if (p_id == gid) {
1893  numer = sudakov_Pgg(t0, t, loc_a, nu) *
1894  std::pow(sudakov_Pqq(t0, t, loc_a, nu), nf);
1895 
1896  if ((is != 1) &&
1897  (is !=
1898  2)) // there is almost no use of is = 2, the `is' in this function is redundant, just for consistency
1899  {
1900  throw std::runtime_error(" error in isp ");
1901  }
1902  } else {
1903  if ((is != 0) &&
1904  (is !=
1905  3)) // there is almost no use of is = 3, the `is' in this function is redundant, just for consistency
1906  {
1907  throw std::runtime_error("error in isp in quark split");
1908  }
1909  double relative_charge = 0.0;
1910  if ((std::abs(p_id) > 0) && (std::abs(p_id) < 4))
1911  relative_charge = 1.0 / 9.0;
1912  if (std::abs(p_id) == 2)
1913  relative_charge = relative_charge * 4.0;
1914 
1915  numer = sudakov_Pqg(t0, t, loc_a, nu) *
1916  std::pow(sudakov_Pqp(t0, t, loc_a, nu), relative_charge);
1917  }
1918 
1919  t_mid = t_low;
1920 
1921  if (numer > r) {
1922  // cout << " numer > r, i.e. ; " << numer << " > " << r << endl ;
1923  return (t_mid);
1924  }
1925 
1926  //t_mid = (t_low+t_hi)/2.0 ;
1927 
1928  scale = t0;
1929 
1930  // cout << " s_approx, s_error = " << s_approx << " " << s_error << endl;
1931 
1932  do {
1933 
1934  t_mid = (t_low + t_hi) / 2.0;
1935 
1936  if (p_id == gid) {
1937  denom = sudakov_Pgg(t0, t_mid, loc_a, nu) *
1938  std::pow(sudakov_Pqq(t0, t_mid, loc_a, nu), nf);
1939  if ((is != 1) && (is != 2)) {
1940  throw std::runtime_error(" error in isp numerator");
1941  }
1942 
1943  } else {
1944  if ((is != 0) && (is != 3)) {
1945  throw std::runtime_error(" error in isp in quark split numerator ");
1946  }
1947  double relative_charge = 0.0;
1948  if ((std::abs(p_id) > 0) && (std::abs(p_id) < 4))
1949  relative_charge = 1.0 / 9.0;
1950  if (std::abs(p_id) == 2)
1951  relative_charge = relative_charge * 4.0;
1952 
1953  denom = sudakov_Pqg(t0, t_mid, loc_a, nu) *
1954  std::pow(sudakov_Pqp(t0, t_mid, loc_a, nu), relative_charge);
1955  }
1956 
1957  ratio = numer / denom;
1958 
1959  diff = (ratio - r) / r;
1960 
1961  // cout << "num, den, r = " << numer << " "<< denom << " " << r << " " << endl;
1962  // cout << "diff, t_mid = " << diff << " " << t_mid << endl;
1963  // cout << " t_low, t_hi = " << t_low << " " << t_hi << endl;
1964  // cin >> test ;
1965 
1966  if (diff < 0.0) {
1967  t_low = t_mid;
1968  //t_mid = (t_low + t_hi)/2.0;
1969  } else {
1970  t_hi = t_mid;
1971  //t_mid = (t_low + t_hi)/2.0;
1972  }
1973 
1974  } while ((abs(diff) > s_approx) && (abs(t_hi - t_low) / t_hi > s_error));
1975 
1976  return (t_mid);
1977 }
1978 
1979 double Matter::generate_vac_t_w_M(int p_id, double M, double nu, double t0,
1980  double t, double loc_a, int is) {
1981  double r, z, ratio, diff, scale, t_low_M0, t_low_MM, t_low_00, t_hi_M0,
1982  t_hi_MM, t_hi_00, t_mid_M0, t_mid_MM, t_mid_00, numer, denom, test;
1983 
1984  double M_charm = PythiaFunction.particleData.m0(cid);
1985  double M_bottom = PythiaFunction.particleData.m0(bid);
1986 
1987  // r = double(random())/ (maxN );
1988  r = ZeroOneDistribution(*GetMt19937Generator());
1989  // r = mtrand1();
1990 
1991  if ((r >= 1.0) || (r <= 0.0)) {
1992  throw std::runtime_error(
1993  "error in random number in t *GetMt19937Generator()");
1994  }
1995 
1996  ratio = 1.0;
1997 
1998  diff = (ratio - r) / r;
1999 
2000  // if (t0<M*M) t0 = M*M;
2001 
2002  t_low_M0 = t0 * (1.0 + std::sqrt(1.0 + 2.0 * M * M / t0));
2003  t_low_MM = 2.0 * (M * M + t0);
2004  t_low_00 = 2.0 * t0;
2005 
2006  t_hi_M0 = t;
2007  t_hi_MM = t;
2008  t_hi_00 = t;
2009 
2010  tscale = t; //for virtuality dependent q-hat
2011 
2012  VERBOSE(1) << MAGENTA << " in gen_vac_t_w_M : t_low , t_hi = " << t_low_M0 << " " << t ;
2013  //cin >> test ;
2014 
2015  if (p_id == gid) {
2016  if (t < 2.0 * (M_charm * M_charm + t0)) {
2017  numer = sudakov_Pgg(t0, t, loc_a, nu) *
2018  std::pow(sudakov_Pqq(t0, t, loc_a, nu), 3.0);
2019  } else if (t < 2.0 * (M_bottom * M_bottom + t0)) {
2020  numer = sudakov_Pgg(t0, t, loc_a, nu) *
2021  std::pow(sudakov_Pqq(t0, t, loc_a, nu), 3.0) *
2022  sudakov_Pqq_w_M_vac_only(M_charm, t0, t, loc_a, nu);
2023  } else {
2024  numer = sudakov_Pgg(t0, t, loc_a, nu) *
2025  std::pow(sudakov_Pqq(t0, t, loc_a, nu), 3.0) *
2026  sudakov_Pqq_w_M_vac_only(M_charm, t0, t, loc_a, nu) *
2027  sudakov_Pqq_w_M_vac_only(M_bottom, t0, t, loc_a, nu);
2028  }
2029  // JSINFO << BOLDYELLOW << " numer = " << numer ;
2030 
2031  // double blurb;
2032  // cin >> blurb;
2033 
2034  if ((is != 1) && (is != 2) && (is != 4) && (is != 5)) {
2035  throw std::runtime_error(" error in isp ");
2036  }
2037  } else {
2038  if (is != 0) {
2039  throw std::runtime_error("error in isp in quark split");
2040  }
2041  if (((int)std::abs((double)p_id)) == 4 ||
2042  ((int)std::abs((double)p_id)) == 5) {
2043  numer = sudakov_Pqg_w_M(M, t0, t, loc_a, nu);
2044 
2045  // std::ofstream Sud_dist;
2046  // Sud_dist.open("Sud_dist.dat", std::ios::app);
2047  // Sud_dist << t << " " << numer << " " << sudakov_Pqg(t0,t,loc_a,nu) << " " << r << " t_low_MO = " << t_low_M0 << endl;
2048  // Sud_dist.close();
2049  } else {
2050  numer = sudakov_Pqg(t0, t, loc_a, nu);
2051  }
2052  }
2053 
2054  t_mid_M0 = t_low_M0;
2055  t_mid_MM = t_low_MM;
2056  t_mid_00 = t_low_00;
2057 
2058  if (numer > r) {
2059  // cout << " numer > r, i.e. ; " << numer << " > " << r << endl ;
2060  if (std::fabs(p_id) == cid || std::fabs(p_id) == bid)
2061  return (t_mid_M0);
2062  else
2063  return (t_mid_00);
2064  }
2065 
2066  //t_mid = (t_low+t_hi)/2.0 ;
2067 
2068  scale = t0;
2069 
2070  // cout << " s_approx, s_error = " << s_approx << " " << s_error << endl;
2071 
2072  bool exit_condition = true;
2073 
2074  do {
2075  t_mid_M0 = (t_low_M0 + t_hi_M0) / 2.0;
2076  t_mid_MM = (t_low_MM + t_hi_MM) / 2.0;
2077  t_mid_00 = (t_low_00 + t_hi_00) / 2.0;
2078 
2079  if (p_id == gid) {
2080  if (t_mid_00 < 2.0 * (M_charm * M_charm + t0)) {
2081  denom = sudakov_Pgg(t0, t_mid_00, loc_a, nu) *
2082  std::pow(sudakov_Pqq(t0, t_mid_00, loc_a, nu), 3.0);
2083  } else if (t_mid_00 < 2.0 * (M_bottom * M_bottom + t0)) {
2084  denom = sudakov_Pgg(t0, t_mid_00, loc_a, nu) *
2085  std::pow(sudakov_Pqq(t0, t_mid_00, loc_a, nu), 3.0) *
2086  sudakov_Pqq_w_M_vac_only(M_charm, t0, t_mid_00, loc_a, nu);
2087  } else {
2088  denom = sudakov_Pgg(t0, t_mid_00, loc_a, nu) *
2089  std::pow(sudakov_Pqq(t0, t_mid_00, loc_a, nu), 3.0) *
2090  sudakov_Pqq_w_M_vac_only(M_charm, t0, t_mid_00, loc_a, nu) *
2091  sudakov_Pqq_w_M_vac_only(M_bottom, t0, t_mid_00, loc_a, nu);
2092  }
2093  if ((is != 1) && (is != 2) && (is != 4) && (is != 5)) {
2094  throw std::runtime_error(" error in isp numerator");
2095  }
2096  } else {
2097  if (is != 0) {
2098  throw std::runtime_error(" error in isp in quark split numerator ");
2099  }
2100  if (((int)std::abs((double)p_id)) == 4 ||
2101  ((int)std::abs((double)p_id)) == 5) {
2102  VERBOSE(8) << BOLDYELLOW << " In generate_vac_t_w_M, t0 = " << t0
2103  << " t_mid_M0 = " << t_mid_M0;
2104  denom = sudakov_Pqg_w_M(M, t0, t_mid_M0, loc_a, nu);
2105  } else {
2106  denom = sudakov_Pqg(t0, t_mid_00, loc_a, nu);
2107  }
2108  }
2109 
2110  ratio = numer / denom;
2111 
2112  diff = (ratio - r) / r;
2113 
2114  if (diff < 0.0) {
2115  t_low_M0 = t_mid_M0;
2116  t_low_MM = t_mid_MM;
2117  t_low_00 = t_mid_00;
2118 
2119  //t_mid = (t_low + t_hi)/2.0;
2120  } else {
2121  t_hi_M0 = t_mid_M0;
2122  t_hi_MM = t_mid_MM;
2123  t_hi_00 = t_mid_00;
2124 
2125  //t_mid = (t_low + t_hi)/2.0;
2126  }
2127 
2128  if (std::abs(p_id) == cid || std::abs(p_id) == bid)
2129  exit_condition = (std::abs(diff) < s_approx) ||
2130  (std::abs(t_hi_M0 - t_low_M0) / t_hi_M0 < s_error);
2131  if (p_id == gid)
2132  exit_condition = (std::abs(diff) < s_approx) ||
2133  (std::abs(t_hi_00 - t_low_00) / t_hi_00 < s_error);
2134  // need to think about the second statement in the gluon exit condition.
2135 
2136  } while (!exit_condition);
2137 
2138  if (std::fabs(p_id) == cid || std::fabs(p_id) == bid)
2139  return (t_mid_M0);
2140  else
2141  return (t_mid_00);
2142 }
2143 
2144 /*
2145 
2146 
2147 
2148  New function
2149 
2150 
2151 
2152 */
2153 
2154 double Matter::generate_vac_z(int p_id, double t0, double t, double loc_b,
2155  double nu, int is) {
2156  double r, z, ratio, diff, e, numer1, numer2, numer, denom, z_low, z_hi, z_mid,
2157  test;
2158 
2159  r = ZeroOneDistribution(*GetMt19937Generator());
2160 
2161  if ((r > 1) || (r < 0)) {
2162  throw std::runtime_error(
2163  " error in random number in z *GetMt19937Generator()");
2164  }
2165 
2166  ratio = 1.0;
2167 
2168  diff = (ratio - r) / r;
2169 
2170  e = t0 / t;
2171 
2172  if (e > 0.5) {
2173  throw std::runtime_error(" error in epsilon");
2174  }
2175 
2176  z_low = e;
2177 
2178  z_hi = double(1.0) - e;
2179 
2180  if (p_id == gid) {
2181  if (is == 1) {
2182  denom = P_z_gg_int(z_low, z_hi, loc_b, t, 2.0 * nu / t, nu);
2183  } else {
2184  denom = P_z_qq_int(z_low, z_hi, loc_b, t, 2.0 * nu / t, nu);
2185  }
2186 
2187  } else if ((p_id != gid) && (is == 0)) {
2188  denom = P_z_qg_int(z_low, z_hi, loc_b, t, 2.0 * nu / t, nu);
2189  } else if ((p_id != gid) && (is == 3)) {
2190  denom = P_z_qp_int(z_low, z_hi, loc_b, t, 2.0 * nu / t, nu);
2191  } else {
2192  throw std::runtime_error(
2193  " I do not understand your combination of particle id and split id in "
2194  "denominator of generate_z ");
2195  }
2196 
2197  //z_mid = (z_low + z_hi)/2.0 ;
2198 
2199  int itcounter = 0;
2200  // cout << " generate_vac_z called with p_id = " << p_id << " t0 = " << t0 << " t = " << t << " loc_b=" << loc_b<< " nu = " << nu << " is = " << is << endl;
2201 
2202  do { // Getting stuck in here for some reason
2203  if (itcounter++ > 10000) {
2204  cout << " in here"
2205  << " abs(diff) = " << abs(diff) << " approx = " << approx
2206  << " r = " << r << " zmid = " << z_mid << " denom = " << denom
2207  << " numer = " << numer << " e = " << e << " " << numer / denom
2208  << endl;
2209  throw std::runtime_error("Stuck in endless loop");
2210  }
2211 
2212  z_mid = (z_low + z_hi) / 2.0;
2213 
2214  if (p_id == gid) {
2215  if (is == 1) {
2216  numer = P_z_gg_int(e, z_mid, loc_b, t, 2.0 * nu / t, nu);
2217  } else {
2218  numer = P_z_qq_int(e, z_mid, loc_b, t, 2.0 * nu / t, nu);
2219  }
2220  } else if ((p_id != gid) && (is == 0)) {
2221  numer = P_z_qg_int(e, z_mid, loc_b, t, 2.0 * nu / t, nu);
2222  } else if ((p_id != gid) && (is == 3)) {
2223  numer = P_z_qp_int(e, z_mid, loc_b, t, 2.0 * nu / t, nu);
2224  } else {
2225  throw std::runtime_error(
2226  " I do not understand your combination of particle id and split id "
2227  "in numerator of generate_z ");
2228  }
2229  ratio = numer / denom;
2230  diff = (ratio - r) / r;
2231 
2232  if (diff > 0.0) {
2233  z_hi = z_mid;
2234  //z_mid = (z_low + z_hi)/2.0;
2235  } else {
2236  z_low = z_mid;
2237  //z_mid = (z_low + z_hi)/2.0 ;
2238  }
2239 
2240  } while ((abs(diff) > approx) && (abs(z_hi - z_low) / z_hi > s_error));
2241 
2242  return (z_mid);
2243 }
2244 
2245 double Matter::generate_vac_z_w_M(int p_id, double M, double t0, double t,
2246  double loc_b, double nu, int is) {
2247  double r, z, ratio, diff, e, numer1, numer2, numer, denom, z_low_00, z_low_M0,
2248  z_low_MM, z_hi_00, z_hi_M0, z_hi_MM, z_mid_00, z_mid_M0, z_mid_MM, test,
2249  z_min_00, z_min_M0, z_min_MM;
2250 
2251  r = ZeroOneDistribution(*GetMt19937Generator());
2252 
2253  // if (t0<M*M) t0 = M*M;
2254 
2255  if ((r > 1) || (r < 0)) {
2256  throw std::runtime_error(
2257  " error in random number in z *GetMt19937Generator()");
2258  }
2259 
2260  ratio = 1.0;
2261 
2262  diff = (ratio - r) / r;
2263 
2264  e = t0 / t;
2265 
2266  double M2 = M * M;
2267 
2268  if (e > 0.5) {
2269  throw std::runtime_error(" error in epsilon");
2270  }
2271 
2272  z_low_M0 = e + M2 / (t + M2);
2273  z_low_00 = e;
2274  z_low_MM = e + M2 / t;
2275 
2276  z_hi_00 = double(1.0) - e;
2277  z_hi_M0 = double(1.0) - e;
2278  z_hi_MM = double(1.0) - e - M2 / t;
2279 
2280  z_min_00 = z_low_00;
2281  z_min_M0 = z_low_M0;
2282  z_min_MM = z_low_MM;
2283 
2284  // JSINFO << BOLDYELLOW << " r = " << r << " 1st z_low_00 = " << z_low_00 << "1st z_hi_00 = " << z_hi_00 << " M = " << M << " is = " << is << " pid = " << p_id;
2285 
2286  if (p_id == gid) {
2287  if (is == 1) // for a gluon to 2 gluons
2288  {
2289  denom = P_z_gg_int(z_min_00, z_hi_00, loc_b, t, 2.0 * nu / t, nu);
2290  //JSINFO << MAGENTA << " denom = " << denom ;
2291  } else {
2292  if (is > 3) //THIS IS FOR GLUON-> HEAVY Q-QBAR
2293  {
2294  denom = P_z_qq_int_w_M_vac_only(M, z_min_MM, z_hi_MM, loc_b, t,
2295  2.0 * nu / t, nu);
2296  } else if (is == 2) // For a gluon to light quark anti-quark
2297  {
2298  denom = P_z_qq_int(z_min_00, z_hi_00, loc_b, t, 2.0 * nu / t, nu);
2299  }
2300  }
2301  } else {
2302  if ((std::abs(p_id) == 4) || (std::abs(p_id) == 5)) {
2303  denom = P_z_qg_int_w_M(M, z_min_M0, z_hi_M0, loc_b, t, 2.0 * nu / t, nu);
2304  } else {
2305  denom = P_z_qg_int(z_min_00, z_hi_00, loc_b, t, 2.0 * nu / t, nu);
2306  }
2307  }
2308 
2309  //z_mid = (z_low + z_hi)/2.0 ;
2310 
2311  int itcounter = 0;
2312  //JSINFO << BOLDYELLOW << " generate_vac_z called with p_id = " << p_id << " t0 = " << t0 << " t = " << t << " loc_b=" << loc_b<< " nu = " << nu << " is = " << is ;
2313 
2314  do { // Getting stuck in here for some reason
2315 
2316  //JSINFO << BOLDYELLOW << " r = " << r << " z_low_00 in loop = " << z_low_00 << " z_hi_00 in loop = " << z_hi_00 << " M = " << M << " is = " << is << " pid = " << p_id;
2317 
2318  if (itcounter++ > 10000) {
2319  cout << " in here"
2320  << " abs(diff) = " << abs(diff) << " approx = " << approx
2321  << " r = " << r << " zmid = " << z_mid_M0 << " denom = " << denom
2322  << " numer = " << numer << " e = " << e << " " << numer / denom
2323  << endl;
2324  throw std::runtime_error("Stuck in endless loop");
2325  }
2326 
2327  z_mid_00 = (z_low_00 + z_hi_00) / 2.0;
2328  z_mid_M0 = (z_low_M0 + z_hi_M0) / 2.0;
2329  z_mid_MM = (z_low_MM + z_hi_MM) / 2.0;
2330 
2331  if (p_id == gid) {
2332  if (is == 1) {
2333  numer = P_z_gg_int(z_min_00, z_mid_00, loc_b, t, 2.0 * nu / t, nu);
2334  //JSINFO << MAGENTA << " numer = " << numer ;
2335  } else {
2336  if (is >
2337  3) // 3 is quark radiating a photon, 4,5 is gluon radiating a c cbar, b bbar
2338  {
2339  numer = P_z_qq_int_w_M_vac_only(M, z_min_MM, z_mid_MM, loc_b, t,
2340  2.0 * nu / t, nu);
2341  } else if (is == 2) // gluon to light quark anti-quark
2342  {
2343  numer = P_z_qq_int(z_min_00, z_mid_00, loc_b, t, 2.0 * nu / t, nu);
2344  }
2345  }
2346  } else {
2347 
2348  if ((std::abs(p_id) == 4) || (std::abs(p_id) == 5)) {
2349  numer =
2350  P_z_qg_int_w_M(M, z_min_M0, z_mid_M0, loc_b, t, 2.0 * nu / t, nu);
2351  } else {
2352  numer = P_z_qg_int(z_min_00, z_mid_00, loc_b, t, 2.0 * nu / t, nu);
2353  }
2354  }
2355  ratio = numer / denom;
2356  diff = (ratio - r) / r;
2357  //JSINFO<< BOLDYELLOW << "num = " << numer << " denom = "<< denom << " ratio = " << numer/denom << " r = " << r ;
2358  // JSINFO << BOLDYELLOW << " diff = " << diff << " z_mid = " << z_mid_00 ;
2359  // cin >> test ;
2360 
2361  if ((diff == 0.0) && ((p_id == gid) || (std::abs(p_id) < 4)))
2362  return (z_mid_00);
2363  if ((diff == 0.0) && ((std::abs(p_id) == 4) || (std::abs(p_id) == 5)))
2364  return (z_mid_M0);
2365 
2366  if (diff > 0.0) {
2367  z_hi_M0 = z_mid_M0;
2368  z_hi_00 = z_mid_00;
2369  z_hi_MM = z_mid_MM;
2370  //z_mid = (z_low + z_hi)/2.0;
2371  } else {
2372  z_low_M0 = z_mid_M0;
2373  z_low_00 = z_mid_00;
2374  z_low_MM = z_mid_MM;
2375  //z_mid = (z_low + z_hi)/2.0 ;
2376  }
2377 
2378  } while (
2379  ((std::abs(p_id) == cid || std::abs(p_id) == bid) &&
2380  (abs(diff) > approx) && (abs(z_hi_M0 - z_low_M0) / z_hi_M0 > s_error)) ||
2381  ((std::abs(p_id) == uid || std::abs(p_id) == did ||
2382  std::abs(p_id) == sid) &&
2383  (abs(diff) > approx) && (abs(z_hi_00 - z_low_00) / z_hi_00 > s_error)) ||
2384  ((std::abs(p_id) == gid && is >= 1 && is < 3) && (abs(diff) > approx) &&
2385  (abs(z_hi_00 - z_low_00) / z_hi_00 > s_error)) ||
2386  ((std::abs(p_id) == gid && is > 3) && (abs(diff) > approx) &&
2387  (abs(z_hi_MM - z_low_MM) / z_hi_MM > s_error)));
2388 
2389  // JSINFO << BOLDRED << " z_mid_00 = " << z_mid_00 ;
2390 
2391  if (p_id == gid && (is == 1 || is == 2))
2392  return (z_mid_00);
2393  else if (p_id == gid && (is == 4 || is == 5))
2394  return (z_mid_MM);
2395  else if (std::abs(p_id) == uid || std::abs(p_id) == did ||
2396  std::abs(p_id) == sid)
2397  return (z_mid_00);
2398  else
2399  return (z_mid_M0);
2400 }
2401 
2402 double Matter::generate_L(double form_time) {
2403  double r, x_low, x_high, x, diff, span, val, arg, norm;
2404 
2405  // r = double(random())/ (maxN );
2406  r = ZeroOneDistribution(*GetMt19937Generator());
2407  // r = mtrand1();
2408 
2409  if ((r > 1) || (r < 0)) {
2410  throw std::runtime_error(
2411  " error in random number in z *GetMt19937Generator()");
2412  }
2413 
2414  x_low = 0;
2415 
2416  x_high = 8.0 * form_time;
2417  // the value of x_high is slightly arbitrary, the erf function is more or less zero at this distance.
2418  // picking 10*form_time will not lead to any different results
2419 
2420  x = (x_low + x_high) / 2.0;
2421 
2422  span = (x_high - x_low) / x_high;
2423 
2424  arg = x / form_time / std::sqrt(pi);
2425 
2426  val = std::erf(arg);
2427 
2428  diff = std::abs(val - r);
2429 
2430  while ((diff > approx) && (span > error)) {
2431  if ((val - r) > 0.0) {
2432  x_high = x;
2433  } else {
2434  x_low = x;
2435  }
2436 
2437  x = (x_low + x_high) / 2.0;
2438 
2439  arg = x / form_time / std::sqrt(pi);
2440 
2441  val = std::erf(arg);
2442 
2443  diff = std::abs(val - r);
2444 
2445  span = (x_high - x_low) / x_high;
2446  }
2447 
2448  // cout << " random number for dist = " << r << " distance generated = " << x << endl;
2449 
2450  return (x);
2451 }
2452 
2453 double Matter::sudakov_Pgg(double g0, double g1, double loc_c, double E) {
2454  double sud, g;
2455  int blurb;
2456 
2457  sud = 1.0;
2458 
2459  if (g1 < 2.0 * g0) {
2460  cerr << " warning: the lower limit of the sudakov > 1/2 upper limit, "
2461  "returning 1 "
2462  << endl;
2463  cerr << " in sudakov_P glue glue, g0, g1 = " << g0 << " " << g1 << endl;
2464  throw std::runtime_error(" warning: the lower limit of the sudakov > 1/2 "
2465  "upper limit, returning 1");
2466 
2467  return (sud);
2468  }
2469  g = 2.0 * g0;
2470 
2471  if (g1 > g) {
2472 
2473  sud = exp(-1.0 * (Ca / 2.0 / pi) * sud_val_GG(g0, g, g1, loc_c, E));
2474  }
2475  return (sud);
2476 }
2477 
2478 double Matter::sud_val_GG(double h0, double h1, double h2, double loc_d,
2479  double E1) {
2480  double val, h, intg, hL, hR, diff, intg_L, intg_R, t_form, span;
2481 
2482  val = 0.0;
2483 
2484  h = (h1 + h2) / 2.0;
2485 
2486  span = (h2 - h1) / h2;
2487 
2488  t_form = 2.0 * E1 / h;
2489 
2490  val = alpha_s(h) * sud_z_GG(h0, h, loc_d, t_form, E1);
2491 
2492  intg = val * (h2 - h1);
2493 
2494  hL = (h1 + h) / 2.0;
2495 
2496  t_form = 2.0 * E1 / hL;
2497 
2498  intg_L = alpha_s(hL) * sud_z_GG(h0, hL, loc_d, t_form, E1) * (h - h1);
2499 
2500  hR = (h + h2) / 2.0;
2501 
2502  t_form = 2.0 * E1 / hR;
2503 
2504  intg_R = alpha_s(hR) * sud_z_GG(h0, hR, loc_d, t_form, E1) * (h2 - h);
2505 
2506  diff = std::abs((intg_L + intg_R - intg) / intg);
2507 
2508  // cout << " iline, gap, diff = " << i_line << " " << h2 << " " << h1 << " " << diff << endl ;
2509  // cout << " intg, Left , right = " << intg << " " << intg_L << " " << intg_R << endl;
2510 
2511  if ((diff > approx) || (span > error)) {
2512  intg = sud_val_GG(h0, h1, h, loc_d, E1) + sud_val_GG(h0, h, h2, loc_d, E1);
2513  }
2514 
2515  // cout << " returning with intg = " << intg << endl;
2516 
2517  return (intg);
2518 }
2519 
2520 double Matter::sud_z_GG(double cg, double cg1, double loc_e, double l_fac,
2521  double E2) {
2522 
2523  double t2, t3, t7, t11, t12, t15, t21, t25, q2, q3, q8, q12, qL, tau, res,
2524  z_min, limit_factor, lz, uz, mz, m_fac;
2525  double t_q1, t_q3, t_q4, t_q6, t_q8, t_q9, t_q12, q_q1, q_q4, q_q6, q_q9,
2526  q_q11;
2527 
2528  z_min = std::sqrt(2) * E_minimum / E2;
2529 
2530  if (cg1 < 2.0 * cg) {
2531 
2532  // cout << " returning with cg, cg1 = " << cg << " " << cg1 << " " << E_minimum << " " << E2 << endl ;
2533  return (0.0);
2534  };
2535 
2536  t2 = std::pow(cg1, 2);
2537  t3 = t2 * cg1;
2538  t7 = std::log(cg);
2539  t11 = std::abs(cg - cg1);
2540  t12 = std::log(t11);
2541  t15 = std::pow(cg, 2);
2542  t21 = t2 * t2;
2543  t25 = -(5.0 * t3 - 12.0 * cg * t2 + 6.0 * t7 * t3 - 6.0 * t12 * t3 -
2544  4.0 * t15 * cg + 6.0 * t15 * cg1) /
2545  t21 / 3.0;
2546 
2547  res = t25;
2548 
2549  limit_factor = 2.0 * std::sqrt(2.0) * cg1 / E2 / 0.1;
2550 
2551  if (limit_factor < 0.0) {
2552  cerr << " error in z limit factor for medium calculation in sud-z-gg = "
2553  << limit_factor << endl;
2554  throw std::runtime_error(
2555  "error in z limit factor for medium calculation in sud-z-gg");
2556  }
2557 
2558  q2 = 1.0 / cg1;
2559  q3 = cg * q2;
2560  q8 = 1.0 - q3;
2561  q12 = (2.0 - 4.0 * q3 + 2.0 / cg * cg1 - 2.0 / q8) * q2;
2562 
2563  if (q12 < 0.0) {
2564  cerr << "ERROR: medium contribution negative in sud_z_GG: q12 = " << q12
2565  << endl;
2566  cerr << "cg, cg1 = " << cg << " " << cg1 << endl;
2567  cerr << " t25 = " << t25 << endl;
2568  throw std::runtime_error("ERROR: medium contribution negative in sud_z_GG");
2569  }
2570 
2571  tau = l_fac;
2572 
2573  if ((length - loc_e) < tau)
2574  tau = (length - loc_e);
2575 
2576  if (loc_e > length)
2577  tau = 0.0;
2578 
2579  m_fac = 1.0;
2580 
2581  // SC
2582  //qL = m_fac*qhat*0.6*tau*profile(loc_e+tau) ;
2583  if (tau < rounding_error) {
2584  qL = 0.0;
2585  } else {
2586  qhat = fncAvrQhat(loc_e, tau);
2587  qL = qhat * tau;
2588  }
2589 
2590  res = t25 + 2.0 * qL * q12 / cg1;
2591  // }
2592  // else{
2593  // cout << " z trap for medium enabled in sud-val-z " << endl ;
2594  // }
2595  // }
2596 
2597  return (res);
2598 }
2599 
2600 double Matter::P_z_gg_int(double cg, double cg1, double loc_e, double cg3,
2601  double l_fac, double E2) {
2602 
2603  double t3, t4, t5, t10, t11, t12, t15, t9, qL, tau, res, limit_factor, lz, uz,
2604  m_fac;
2605 
2606  //if ((cg< cg1/(2.0*E2*E2/cg1+1.0) )) cg = cg1/( 2.0*E2*E2/cg1 + 1.0 );
2607 
2608  t3 = std::log((1.0 - cg1));
2609  t4 = std::log(cg1);
2610  t5 = std::pow(cg1, 2);
2611  t10 = std::log((1.0 - cg));
2612  t11 = std::log(cg);
2613  t12 = std::pow(cg, 2);
2614  t15 = -(2.0 * cg1) - t3 + t4 - (2.0 / 3.0) * t5 * cg1 + t5 + (2.0 * cg) +
2615  t10 - t11 + (2.0 / 3.0) * t12 * cg - t12;
2616 
2617  res = t15;
2618 
2619  limit_factor = 2.0 * std::sqrt(2.0) * cg3 / E2 / 0.1;
2620 
2621  if (limit_factor < 0.0) {
2622  cerr << " error in z limit factor for medium calculation = " << limit_factor
2623  << endl;
2624  throw std::runtime_error(" error in z limit factor for medium calculation");
2625  }
2626 
2627  tau = l_fac;
2628 
2629  if ((length - loc_e) < tau)
2630  tau = (length - loc_e);
2631 
2632  // if ((length - loc_e) < tau) tau = 0;
2633 
2634  if (loc_e > length)
2635  tau = 0.0;
2636 
2637  m_fac = 1.0;
2638 
2639  // if ((qhat*tau<1.0)&&(in_vac==false)) m_fac = 1.0/qhat/tau;
2640 
2641  // SC
2642  //qL = m_fac*qhat*0.6*tau*profile(loc_e + tau) ;
2643  if (tau < rounding_error) {
2644  qL = 0.0;
2645  } else {
2646  qhat = fncAvrQhat(loc_e, tau);
2647  qL = qhat * tau;
2648  }
2649 
2650  t9 = 2.0 * cg1 - 1.0 / (-1.0 + cg1) - 1.0 / cg1 - 2.0 * cg +
2651  1.0 / (-1.0 + cg) + 1.0 / cg;
2652 
2653  if (t9 < 0.0) {
2654  cerr << "ERROR: medium contribution negative in P_z_gg_int : t9 = " << t9
2655  << endl;
2656 
2657  cerr << " cg, cg1 = " << cg << " " << cg1 << endl;
2658  throw std::runtime_error(
2659  "ERROR: medium contribution negative in P_z_gg_int");
2660  }
2661 
2662  res = t15 + 2.0 * t9 * qL / cg3;
2663 
2664  return (res);
2665 }
2666 
2667 double Matter::sudakov_Pqq(double q0, double q1, double loc_c, double E) {
2668  double sud, q;
2669 
2670  sud = 1.0;
2671 
2672  if (q1 < 2.0 * q0)
2673  // if (g1<g0)
2674  {
2675  JSWARN << " warning: the lower limit of the sudakov > 1/2 upper limit, "
2676  "returning 1 ";
2677  JSWARN << " in sudakov_Pquark quark, q0, q1 = " << q0 << " " << q1;
2678  return (sud);
2679  }
2680  q = 2.0 * q0;
2681 
2682  // g = g0 ;
2683 
2684  sud = exp(-1.0 * (Tf / 2.0 / pi) * sud_val_QQ(q0, q, q1, loc_c, E));
2685 
2686  return (sud);
2687 }
2688 
2689 double Matter::sudakov_Pqq_w_M_vac_only(double M, double q0, double q1,
2690  double loc_c, double E) {
2691  double sud, q;
2692 
2693  sud = 1.0;
2694 
2695  if (q1 < 2.0 * (q0 + M * M)) {
2696  JSWARN << " warning: the upper limit of the sudakov q1<2.0*(q0+M*M), "
2697  "returning 1 ";
2698  JSWARN << " in sudakov_Pquark quark, q0, q1 = " << q0 << " " << q1;
2699  return (sud);
2700  } else {
2701  q = 2.0 * (q0 + M * M);
2702  sud = exp(-1.0 * (Tf / 2.0 / pi) *
2703  sud_val_QQ_w_M_vac_only(M, q0, q, q1, loc_c, E));
2704  return (sud);
2705  }
2706 }
2707 
2708 double Matter::sud_val_QQ(double h0, double h1, double h2, double loc_d,
2709  double E1) {
2710  double val, h, intg, hL, hR, diff, intg_L, intg_R, t_form, span;
2711 
2712  h = (h1 + h2) / 2.0;
2713 
2714  span = (h2 - h1) / h2;
2715 
2716  t_form = 2.0 * E1 / h;
2717 
2718  val = alpha_s(h) * sud_z_QQ(h0, h, loc_d, t_form, E1);
2719 
2720  intg = val * (h2 - h1);
2721 
2722  hL = (h1 + h) / 2.0;
2723 
2724  t_form = 2.0 * E1 / hL;
2725 
2726  intg_L = alpha_s(hL) * sud_z_QQ(h0, hL, loc_d, t_form, E1) * (h - h1);
2727 
2728  hR = (h + h2) / 2.0;
2729 
2730  t_form = 2.0 * E1 / hR;
2731 
2732  intg_R = alpha_s(hR) * sud_z_QQ(h0, hR, loc_d, t_form, E1) * (h2 - h);
2733 
2734  diff = std::abs((intg_L + intg_R - intg) / intg);
2735 
2736  //JSINFO << BOLDYELLOW << " h2, h1, diff = " << h2 << " , " << h1 << " , " << diff ;
2737  //JSINFO << BOLDYELLOW << " intg, Left , right = " << intg << " " << intg_L << " " << intg_R ;
2738 
2739  if ((diff > approx) || (span > error)) {
2740  intg = sud_val_QQ(h0, h1, h, loc_d, E1) + sud_val_QQ(h0, h, h2, loc_d, E1);
2741  }
2742 
2743  // cout << " returning with intg = " << intg << endl;
2744 
2745  return (intg);
2746 }
2747 
2748 double Matter::sud_val_QQ_w_M_vac_only(double M, double h0, double h1,
2749  double h2, double loc_d, double E1) {
2750  double val, h, intg, hL, hR, diff, intg_L, intg_R, t_form, span;
2751 
2752  h = (h1 + h2) / 2.0;
2753 
2754  span = (h2 - h1) / h2;
2755 
2756  t_form = 2.0 * E1 / h;
2757 
2758  val = alpha_s(h) * sud_z_QQ_w_M_vac_only(M, h0, h, loc_d, t_form, E1);
2759 
2760  intg = val * (h2 - h1);
2761 
2762  hL = (h1 + h) / 2.0;
2763 
2764  t_form = 2.0 * E1 / hL;
2765 
2766  intg_L = alpha_s(hL) * sud_z_QQ_w_M_vac_only(M, h0, hL, loc_d, t_form, E1) *
2767  (h - h1);
2768 
2769  hR = (h + h2) / 2.0;
2770 
2771  t_form = 2.0 * E1 / hR;
2772 
2773  intg_R = alpha_s(hR) * sud_z_QQ_w_M_vac_only(M, h0, hR, loc_d, t_form, E1) *
2774  (h2 - h);
2775 
2776  diff = std::abs((intg_L + intg_R - intg) / intg);
2777 
2778  // cout << " iline, gap, diff = " << i_line << " " << h2 << " " << h1 << " " << diff << endl ;
2779  // cout << " intg, Left , right = " << intg << " " << intg_L << " " << intg_R << endl;
2780 
2781  if ((diff > approx) || (span > error)) {
2782  intg = sud_val_QQ_w_M_vac_only(M, h0, h1, h, loc_d, E1) +
2783  sud_val_QQ_w_M_vac_only(M, h0, h, h2, loc_d, E1);
2784  }
2785 
2786  // cout << " returning with intg = " << intg << endl;
2787 
2788  return (intg);
2789 }
2790 
2791 double Matter::sud_z_QQ(double cg, double cg1, double loc_e, double l_fac,
2792  double E2) {
2793 
2794  double t2, t4, t5, t7, t9, t14, q2, q3, q5, q6, q8, q15, qL, tau, res, z_min;
2795 
2796  z_min = std::sqrt(2) * E_minimum / E2;
2797 
2798  // if (cg<cg1*z_min) cg = cg1*z_min;
2799 
2800  // if ((cg< cg1/(2.0*E2*E2/cg1+1.0) )) cg = cg1/( 2.0*E2*E2/cg1 + 1.0 );
2801 
2802  if (cg1 < 2.0 * cg) {
2803 
2804  // cout << " returning with cg, cg1 = " << cg << " " << cg1 << " " << E_minimum << " " << E2 << endl ;
2805  return (0.0);
2806  };
2807 
2808  t2 = 1.0 / cg1;
2809  t4 = 1.0 - cg * t2;
2810  t5 = t4 * t4;
2811  t7 = std::pow(cg, 2.0);
2812  t9 = std::pow(cg1, 2.0);
2813  t14 = ((t5 * t4) - t7 * cg / t9 / cg1) * t2 / 3.0;
2814 
2815  // return(t25);
2816 
2817  q2 = 1.0 / cg1;
2818  q3 = (cg * q2);
2819  q5 = 1.0 - q3;
2820  q6 = std::log(std::abs(q5));
2821  q8 = std::log(q3);
2822  /* q10 = std::log(q3);
2823  q12 = std::log(q5);*/
2824  q15 = (-1.0 + (2.0 * q3) + q6 - q8) * q2;
2825 
2826  if (q15 < 0.0) {
2827  cerr << "ERROR: medium contribution negative in sud_z_QQ: q15 = " << q15
2828  << endl;
2829  cerr << "cg, cg1 = " << cg << " " << cg1 << endl;
2830  cerr << " t14 = " << t14 << endl;
2831  throw std::runtime_error("ERROR: medium contribution negative in sud_z_QQ");
2832  }
2833 
2834  tau = l_fac;
2835 
2836  if ((length - loc_e) < tau)
2837  tau = (length - loc_e);
2838 
2839  if (loc_e > length)
2840  tau = 0.0;
2841 
2842  // SC
2843  //qL = qhat*0.6*tau*profile(loc_e + tau) ;
2844  if (tau < rounding_error) {
2845  qL = 0.0;
2846  } else {
2847  qhat = fncAvrQhat(loc_e, tau);
2848  qL = qhat * tau;
2849  }
2850 
2851  res = t14 + 2.0 * qL * q15 / cg1;
2852 
2853  return (res);
2854 }
2855 
2856 double Matter::sud_z_QQ_w_M_vac_only(double M, double cg, double cg1,
2857  double loc_e, double l_fac, double E2) {
2858 
2859  double q2, q3, q5, q6, q8, q15, qL, tau, res, z_min;
2860 
2861  z_min = std::sqrt(2) * E_minimum / E2;
2862 
2863  // if (cg<cg1*z_min) cg = cg1*z_min;
2864 
2865  // if ((cg< cg1/(2.0*E2*E2/cg1+1.0) )) cg = cg1/( 2.0*E2*E2/cg1 + 1.0 );
2866 
2867  if (cg1 < 2.0 * (cg + M * M)) {
2868 
2869  // cout << " returning with cg, cg1 = " << cg << " " << cg1 << " " << E_minimum << " " << E2 << endl ;
2870  return (0.0);
2871  };
2872 
2873  double t1 = M * M;
2874  double t2 = t1 + cg;
2875  double t3 = 1.0 / cg1;
2876  double t5 = -t2 * t3 + 1.0;
2877  double t6 = t5 * t5;
2878  double t8 = t2 * t2;
2879  double t10 = pow(cg1, 2.0);
2880  double t15 = 2.0 / 3.0 * (t6 * t5 - t8 * t2 / t10 / cg1) * t3;
2881 
2882  q2 = 1.0 / cg1;
2883  q3 = (cg * q2);
2884  q5 = 1.0 - q3;
2885  q6 = std::log(std::abs(q5));
2886  q8 = std::log(q3);
2887  /* q10 = std::log(q3);
2888  q12 = std::log(q5);*/
2889  q15 = (-1.0 + (2.0 * q3) + q6 - q8) * q2;
2890 
2891  if (q15 < 0.0) {
2892  cerr << "ERROR: medium contribution negative in sud_z_QQ: q15 = " << q15
2893  << endl;
2894  cerr << "cg, cg1 = " << cg << " " << cg1 << endl;
2895  cerr << " t15 = " << t15 << endl;
2896  throw std::runtime_error("ERROR: medium contribution negative in sud_z_QQ");
2897  }
2898 
2899  tau = l_fac;
2900 
2901  if ((length - loc_e) < tau)
2902  tau = (length - loc_e);
2903 
2904  if (loc_e > length)
2905  tau = 0.0;
2906 
2907  // SC
2908  //qL = qhat*0.6*tau*profile(loc_e + tau) ;
2909  if (tau < rounding_error) {
2910  qL = 0.0;
2911  } else {
2912  qhat = fncAvrQhat(loc_e, tau);
2913  qL = qhat * tau;
2914  }
2915 
2916  res = t15 + 2.0 * qL * q15 / cg1;
2917 
2918  return (res);
2919 }
2920 
2921 double Matter::P_z_qq_int(double cg, double cg1, double loc_e, double cg3,
2922  double l_fac, double E2) {
2923  double t_q1, t_q3, t_q4, t_q6, t_q8, t_q9, t_q12, q_q1, q_q4, q_q6, q_q9,
2924  q_q11, qL, tau, res;
2925 
2926  if ((cg < cg1 / (2.0 * E2 * E2 / cg1 + 1.0)))
2927  cg = cg1 / (2.0 * E2 * E2 / cg1 + 1.0);
2928 
2929  t_q1 = std::pow(cg1, 2);
2930  t_q3 = 1.0 - cg1;
2931  t_q4 = t_q3 * t_q3;
2932  t_q6 = std::pow(cg, 2);
2933  t_q8 = 1.0 - cg;
2934  t_q9 = t_q8 * t_q8;
2935  t_q12 = t_q1 * cg1 / 6.0 - t_q4 * t_q3 / 6.0 - t_q6 * cg / 6.0 +
2936  t_q9 * t_q8 / 6.0;
2937 
2938  tau = l_fac;
2939 
2940  if ((length - loc_e) < tau)
2941  tau = (length - loc_e);
2942 
2943  if (loc_e > length)
2944  tau = 0.0;
2945 
2946  // SC
2947  //qL = qhat*0.6*tau*profile(loc_e + tau) ;
2948  if (tau < rounding_error) {
2949  qL = 0.0;
2950  } else {
2951  qhat = fncAvrQhat(loc_e, tau);
2952  qL = qhat * tau;
2953  }
2954 
2955  q_q1 = std::log(cg1);
2956  q_q4 = std::log(1.0 - cg1);
2957  q_q6 = std::log(cg);
2958  q_q9 = std::log(1.0 - cg);
2959  q_q11 = -cg1 + q_q1 / 2.0 - q_q4 / 2.0 + cg - q_q6 / 2.0 + q_q9 / 2.0;
2960 
2961  if (q_q11 < 0.0) {
2962  cerr << "ERROR: medium contribution negative in P_z_gg_int : q_q11 = "
2963  << q_q11 << endl;
2964  throw std::runtime_error(
2965  "ERROR: medium contribution negative in P_z_gg_int");
2966  }
2967 
2968  res = t_q12 * Tf / Ca + 2.0 * qL * q_q11 / cg3 * (Tf * Cf / Ca / Ca);
2969 
2970  return (res);
2971 }
2972 //
2973 //
2974 
2975 // Sudakov for a quark to radiate a quark + photon
2976 double Matter::sudakov_Pqp(double g0, double g1, double loc_c, double E) {
2977  double sud, g;
2978  int blurb;
2979 
2980  sud = 1.0;
2981 
2982  if (g1 < 2.0 * g0) {
2983  JSWARN << " warning: the lower limit of the sudakov > 1/2 upper limit, "
2984  "returning 1 ";
2985  JSWARN << " in sudakov_Pquark Photon, g0, g1 = " << g0 << " " << g1;
2986  return (sud);
2987  }
2988  g = 2.0 * g0;
2989 
2990  double logsud = sud_val_QP(g0, g, g1, loc_c, E);
2991 
2992  sud = exp((-1.0 / 2.0 / pi) * logsud);
2993 
2994  return (sud);
2995 }
2996 
2997 double Matter::sud_val_QP(double h0, double h1, double h2, double loc_d,
2998  double E1) {
2999  double val, h, intg, hL, hR, diff, intg_L, intg_R, t_form, span;
3000  int blurb;
3001 
3002  double alphaEM = 1.0 / 137.0;
3003 
3004  val = 0.0;
3005 
3006  h = (h1 + h2) / 2.0;
3007 
3008  span = (h2 - h1) / h2;
3009 
3010  t_form = 2.0 * E1 / h;
3011 
3012  val = alphaEM * sud_z_QP(h0, h, loc_d, t_form, E1);
3013 
3014  intg = val * (h2 - h1);
3015 
3016  hL = (h1 + h) / 2.0;
3017 
3018  t_form = 2.0 * E1 / hL;
3019 
3020  intg_L = alphaEM * sud_z_QP(h0, hL, loc_d, t_form, E1) * (h - h1);
3021 
3022  hR = (h + h2) / 2.0;
3023 
3024  t_form = 2.0 * E1 / hR;
3025 
3026  intg_R = alphaEM * sud_z_QP(h0, hR, loc_d, t_form, E1) * (h2 - h);
3027 
3028  diff = std::abs((intg_L + intg_R - intg) / intg);
3029 
3030  if ((diff > approx) || (span > error)) {
3031  intg = sud_val_QP(h0, h1, h, loc_d, E1) + sud_val_QP(h0, h, h2, loc_d, E1);
3032  }
3033 
3034  return (intg);
3035 }
3036 
3037 double Matter::P_z_qq_int_w_M_vac_only(double M, double cg, double cg1,
3038  double loc_e, double cg3, double l_fac,
3039  double E2) {
3040  double t_q1, t_q3, t_q4, t_q6, t_q8, t_q9, t_q12, q_q1, q_q4, q_q6, q_q9,
3041  q_q11, qL, tau, res;
3042 
3043  if ((cg < cg1 / (2.0 * E2 * E2 / cg1 + 1.0)))
3044  cg = cg1 / (2.0 * E2 * E2 / cg1 + 1.0);
3045 
3046  t_q1 = std::pow(cg1, 2);
3047  t_q3 = 1.0 - cg1;
3048  t_q4 = t_q3 * t_q3;
3049  t_q6 = std::pow(cg, 2);
3050  t_q8 = 1.0 - cg;
3051  t_q9 = t_q8 * t_q8;
3052  t_q12 = t_q1 * cg1 / 6.0 - t_q4 * t_q3 / 6.0 - t_q6 * cg / 6.0 +
3053  t_q9 * t_q8 / 6.0;
3054 
3055  tau = l_fac;
3056 
3057  if ((length - loc_e) < tau)
3058  tau = (length - loc_e);
3059 
3060  if (loc_e > length)
3061  tau = 0.0;
3062 
3063  // SC
3064  //qL = qhat*0.6*tau*profile(loc_e + tau) ;
3065  if (tau < rounding_error) {
3066  qL = 0.0;
3067  } else {
3068  qhat = fncAvrQhat(loc_e, tau);
3069  qL = qhat * tau;
3070  }
3071 
3072  q_q1 = std::log(cg1);
3073  q_q4 = std::log(1.0 - cg1);
3074  q_q6 = std::log(cg);
3075  q_q9 = std::log(1.0 - cg);
3076  q_q11 = -cg1 + q_q1 / 2.0 - q_q4 / 2.0 + cg - q_q6 / 2.0 + q_q9 / 2.0;
3077 
3078  if (q_q11 < 0.0) {
3079  cerr << "ERROR: medium contribution negative in P_z_gg_int_w_M : q_q11 = "
3080  << q_q11 << endl;
3081  cout << " z_low = " << cg << " z_hi = " << cg1 << endl;
3082  throw std::runtime_error(
3083  "ERROR: medium contribution negative in P_z_gg_int");
3084  }
3085 
3086  res = t_q12 * Tf / Ca + 2.0 * qL * q_q11 / cg3 * (Tf * Cf / Ca / Ca);
3087 
3088  return (res);
3089 }
3090 
3091 double Matter::sud_z_QP(double cg, double cg1, double loc_e, double l_fac,
3092  double E2) {
3093 
3094  double t2, t6, t10, t11, t17, q2, q3, q4, q5, q6, q10, q14, qL, tau, res,
3095  z_min;
3096  int blurb;
3097 
3098  z_min = std::sqrt(2) * E_minimum / E2;
3099 
3100  if (cg1 < 2.0 * cg) {
3101  return (0.0);
3102  };
3103 
3104  t2 = std::pow(cg1, 2);
3105  t6 = std::log(cg);
3106  t10 = std::abs(cg - cg1);
3107  t11 = std::log(t10);
3108  t17 = -1.0 / t2 * (3.0 * cg1 - 6.0 * cg + 4.0 * t6 * cg1 - 4.0 * t11 * cg1) /
3109  2.0;
3110 
3111  // return(t17);
3112 
3113  q14 = 0.0;
3114 
3115  tau = l_fac;
3116 
3117  if ((length - loc_e) < tau)
3118  tau = (length - loc_e);
3119 
3120  if (loc_e > length)
3121  tau = 0.0;
3122 
3123  // SC
3124  //qL = qhat*0.6*tau*profile(loc_e + tau) ;
3125  if (tau < rounding_error) {
3126  qL = 0.0;
3127  } else {
3128  qhat = fncAvrQhat(loc_e, tau) * Cf /
3129  Ca; //for photon production, only the quark scatters
3130  qL = qhat * tau;
3131  }
3132 
3133  //JSINFO << BOLDRED << " qhat L = " << qL << " location = " << loc_e << " tau = " << tau << " length = " << length;
3134 
3135  res = t17 + 2.0 * qL * q14 / cg1;
3136 
3137  // cout << " t0 , t , res = " << cg << " " << cg1 << " " << res << endl ;
3138 
3139  if (q14 < 0.0) {
3140  cerr << "ERROR: medium contribution negative in sud_z_QG : q14 = " << q14
3141  << endl;
3142  throw std::runtime_error("ERROR: medium contribution negative in sud_z_QG");
3143  }
3144 
3145  return (res);
3146 }
3147 double Matter::P_z_qp_int(double cg, double cg1, double loc_e, double cg3,
3148  double l_fac, double E2) {
3149 
3150  double t2, t5, t7, t10, t12, q2, q6, q10, tau, qL, res;
3151 
3152  if ((cg < cg1 / (2.0 * E2 * E2 / cg1 + 1.0)))
3153  cg = cg1 / (2.0 * E2 * E2 / cg1 + 1.0);
3154 
3155  t2 = std::pow(cg1, 2);
3156  t5 = std::log(1.0 - cg1);
3157  t7 = std::pow(cg, 2);
3158  t10 = std::log(1.0 - cg);
3159  t12 = -cg1 - t2 / 2.0 - 2.0 * t5 + cg + t7 / 2.0 + 2.0 * t10;
3160 
3161  // return(t12);
3162 
3163  q10 = 0.0;
3164  tau = l_fac;
3165 
3166  if ((length - loc_e) < tau)
3167  tau = (length - loc_e);
3168 
3169  if (loc_e > length)
3170  tau = 0.0;
3171 
3172  // SC
3173  //qL = qhat*0.6*tau*profile(loc_e + tau) ;
3174  if (tau < rounding_error) {
3175  qL = 0.0;
3176  } else {
3177  qhat = fncAvrQhat(loc_e, tau);
3178  qL = qhat * tau;
3179  }
3180 
3181  res = t12 + 2.0 * qL * q10 / cg3;
3182 
3183  return (res);
3184 }
3185 
3186 //
3187 //
3188 // Sudakov for a quark to radiate a quark + gluon.
3189 double Matter::sudakov_Pqg(double g0, double g1, double loc_c, double E) {
3190  double sud, g;
3191  int blurb;
3192 
3193  sud = 1.0;
3194 
3195  if (g1 < 2.0 * g0) {
3196  JSWARN << " warning: the lower limit of the sudakov > 1/2 upper limit, "
3197  "returning 1 ";
3198  JSWARN << " in sudakov_Pquark gluon, g0, g1 = " << g0 << " " << g1;
3199  return (sud);
3200  }
3201  g = 2.0 * g0;
3202 
3203  sud = exp(-1.0 * (Cf / 2.0 / pi) * sud_val_QG(g0, g, g1, loc_c, E));
3204 
3205  return (sud);
3206 }
3207 
3208 double Matter::sudakov_Pqg_w_M(double M, double g0, double g1, double loc_c,
3209  double E) {
3210  double sud, g;
3211  int blurb;
3212 
3213  sud = 1.0;
3214 
3215  if (g1 < g0 * (1.0 + std::sqrt(1.0 + 2.0 * M * M / g0))) {
3216  JSWARN << " warning: Not enough separation between upper and lower limits "
3217  "of Sudakov to have resolvable radiation ";
3218  JSWARN << " in sudakov_Pquark gluon, g0*( 1.0 + std::sqrt( 1.0 + "
3219  "2.0*M*M/g0 ) ) = "
3220  << g0 * (1.0 + std::sqrt(1.0 + 2.0 * M * M / g0)) << " g1 = " << g1;
3221  JSWARN << " M = " << M;
3222 
3223  return (sud);
3224  }
3225  g = g0 * (1.0 + std::sqrt(1.0 + 2.0 * M * M / g0));
3226 
3227  sud = exp(-1.0 * (Cf / 2.0 / pi) * sud_val_QG_w_M(M, g0, g, g1, loc_c, E));
3228 
3229  return (sud);
3230 }
3231 
3232 double Matter::sud_val_QG(double h0, double h1, double h2, double loc_d,
3233  double E1) {
3234  double val, h, intg, hL, hR, diff, intg_L, intg_R, t_form, span;
3235  int blurb;
3236 
3237  val = 0.0;
3238 
3239  h = (h1 + h2) / 2.0;
3240 
3241  span = (h2 - h1) / h2;
3242 
3243  t_form = 2.0 * E1 / h;
3244 
3245  val = alpha_s(h) * sud_z_QG(h0, h, loc_d, t_form, E1);
3246 
3247  intg = val * (h2 - h1);
3248 
3249  hL = (h1 + h) / 2.0;
3250 
3251  t_form = 2.0 * E1 / hL;
3252 
3253  intg_L = alpha_s(hL) * sud_z_QG(h0, hL, loc_d, t_form, E1) * (h - h1);
3254 
3255  hR = (h + h2) / 2.0;
3256 
3257  t_form = 2.0 * E1 / hR;
3258 
3259  intg_R = alpha_s(hR) * sud_z_QG(h0, hR, loc_d, t_form, E1) * (h2 - h);
3260 
3261  diff = std::abs((intg_L + intg_R - intg) / intg);
3262 
3263  if ((diff > approx) || (span > error)) {
3264  intg = sud_val_QG(h0, h1, h, loc_d, E1) + sud_val_QG(h0, h, h2, loc_d, E1);
3265  }
3266 
3267  return (intg);
3268 }
3269 
3270 double Matter::sud_val_QG_w_M(double M, double h0, double h1, double h2,
3271  double loc_d, double E1) {
3272  double val, h, intg, hL, hR, diff, intg_L, intg_R, t_form, span;
3273  int blurb;
3274 
3275  val = 0.0;
3276 
3277  h = (h1 + h2) / 2.0;
3278 
3279  span = (h2 - h1) / h2;
3280 
3281  t_form = 2.0 * E1 / h;
3282 
3283  val = alpha_s(h) * sud_z_QG_w_M(M, h0, h, loc_d, t_form, E1);
3284 
3285  intg = val * (h2 - h1);
3286 
3287  hL = (h1 + h) / 2.0;
3288 
3289  t_form = 2.0 * E1 / hL;
3290 
3291  intg_L = alpha_s(hL) * sud_z_QG_w_M(M, h0, hL, loc_d, t_form, E1) * (h - h1);
3292 
3293  hR = (h + h2) / 2.0;
3294 
3295  t_form = 2.0 * E1 / hR;
3296 
3297  intg_R = alpha_s(hR) * sud_z_QG_w_M(M, h0, hR, loc_d, t_form, E1) * (h2 - h);
3298 
3299  diff = std::abs((intg_L + intg_R - intg) / intg);
3300 
3301  // cout << " iline, gap, diff = " << i_line << " " << h2 << " " << h1 << " " << diff << endl ;
3302  // cout << " intg, Left , right = " << intg << " " << intg_L << " " << intg_R << endl;
3303 
3304  if ((diff > approx) || (span > error)) {
3305  intg = sud_val_QG_w_M(M, h0, h1, h, loc_d, E1) +
3306  sud_val_QG_w_M(M, h0, h, h2, loc_d, E1);
3307  }
3308 
3309  // cout << " returning with intg = " << intg << endl;
3310 
3311  return (intg);
3312 }
3313 
3314 double Matter::sud_z_QG(double cg, double cg1, double loc_e, double l_fac,
3315  double E2) {
3316 
3317  double t2, t6, t10, t11, t17, q2, q3, q4, q5, q6, q10, q14, qL, tau, res,
3318  z_min;
3319  int blurb;
3320 
3321  z_min = std::sqrt(2) * E_minimum / E2;
3322 
3323  if (cg1 < 2.0 * cg) {
3324  return (0.0);
3325  };
3326 
3327  t2 = std::pow(cg1, 2);
3328  t6 = std::log(cg);
3329  t10 = std::abs(cg - cg1);
3330  t11 = std::log(t10);
3331  t17 = -1.0 / t2 * (3.0 * cg1 - 6.0 * cg + 4.0 * t6 * cg1 - 4.0 * t11 * cg1) /
3332  2.0;
3333 
3334  // return(t17);
3335 
3336  q2 = 1.0 / cg1;
3337  q3 = cg * q2;
3338  q4 = q3 - 1.0;
3339  q5 = std::abs(q4);
3340  q6 = std::log(q5);
3341  q10 = std::log(q3);
3342  q14 = (q6 + 2.0 / cg * cg1 - q10 + 2.0 / q4) * q2;
3343 
3344  tau = l_fac;
3345 
3346  if ((length - loc_e) < tau)
3347  tau = (length - loc_e);
3348 
3349  if (loc_e > length)
3350  tau = 0.0;
3351 
3352  // SC
3353  //qL = qhat*0.6*tau*profile(loc_e + tau) ;
3354  if (tau < rounding_error) {
3355  qL = 0.0;
3356  } else {
3357  qhat = fncAvrQhat(loc_e, tau);
3358  if (qhat * sqrt(2) > 0.6) {
3359  // JSINFO << BOLDYELLOW << " length = " << length << " loc = " << loc_e << " tau = " << tau ;
3360  //JSINFO << BOLDYELLOW << " parton formed at x = " << initRx << " y = " << initRy << " z = " << initRz << " t = " << initR0 ;
3361  // JSINFO << BOLDYELLOW << " mean qhat for sudakov in GeV^2/fm = " << qhat*5*sqrt(2) ;
3362  }
3363  qL = qhat * tau;
3364  }
3365 
3366  //JSINFO << BOLDRED << " qhat L = " << qL << " location = " << loc_e << " tau = " << tau << " length = " << length;
3367 
3368  res = t17 + 2.0 * qL * q14 / cg1;
3369 
3370  // cout << " t0 , t , res = " << cg << " " << cg1 << " " << res << endl ;
3371 
3372  if (q14 < 0.0) {
3373  cerr << "ERROR: medium contribution negative in sud_z_QG : q14 = " << q14
3374  << endl;
3375  throw std::runtime_error("ERROR: medium contribution negative in sud_z_QG");
3376  }
3377 
3378  return (res);
3379 }
3380 
3381 double Matter::sud_z_QG_w_M(double M, double cg, double cg1, double loc_e,
3382  double l_fac, double E2) //(t0,t,loc)
3383 {
3384 
3385  double qL, tau, res, z_min;
3386  int blurb;
3387 
3388  z_min = std::sqrt(2) * E_minimum / E2;
3389 
3390  if (cg1 < 2.0 * cg + M * M / (1.0 + M * M / cg1)) {
3391 
3392  JSINFO << MAGENTA << " returning with cg, cg1 = " << cg << " " << cg1
3393  << " " << E_minimum << " " << E2;
3394  return (M * M);
3395  };
3396 
3397  double t1 = 1.0 / cg1;
3398  double t2 = t1 * cg;
3399  double t4 = std::pow(1.0 - t2, 2.0);
3400  double t7 = std::log(t2);
3401  double t9 = M * M;
3402  double t10 = t1 * t9;
3403  double t13 = 1.0 / (t10 + 1.0) * t10;
3404  double t15 = std::pow(t2 + t13, 2.0);
3405  double t18 = std::log(1.0 - t2 - t13);
3406  double t21 = t1 * (-t4 / 2.0 - 1.0 + 2.0 * t2 - 2.0 * t7 + t15 / 2.0 + t13 +
3407  2.0 * t18);
3408 
3409  double q1 = M * M;
3410  double q2 = 1.0 / cg1;
3411  double q3 = q2 * q1;
3412  double q5 = 4.0 * q3 + 1.0;
3413  double q6 = q2 * cg;
3414  double q7 = std::log(q6);
3415  double q9 = q1 * q1;
3416  double q10 = std::pow(cg1, 2.0);
3417  double q12 = 1.0 / q10 * q9;
3418  double q14 = q12 - 2.0 * q3 + 1.0 / 2.0;
3419  double q15 = 1.0 - q6;
3420  double q16 = std::log(q15);
3421  double q18 = q3 + 1.0;
3422  double q28 = q15 * q15;
3423  double q33 = 1.0 / q18 * q3;
3424  double q34 = 1.0 - q6 - q33;
3425  double q35 = std::log(q34);
3426  double q37 = q6 + q33;
3427  double q38 = std::log(q37);
3428  double q48 = q37 * q37;
3429  double q52 = q7 * q5 + q16 * q14 + q15 * q18 / 2.0 + 2.0 / cg * cg1 +
3430  3.0 / 2.0 * q2 / q15 * q1 - 1.0 / q28 * q12 / 2.0 - q35 * q5 -
3431  q38 * q14 - q37 * q18 / 2.0 - 2.0 / q34 -
3432  3.0 / 2.0 * q2 / q37 * q1 + 1.0 / q48 * q12 / 2.0;
3433  double q53 = q2 * q52;
3434 
3435  tau = l_fac;
3436 
3437  if ((length - loc_e) < tau)
3438  tau = (length - loc_e);
3439 
3440  if (loc_e > length)
3441  tau = 0.0;
3442 
3443  // SC
3444  //qL = qhat*0.6*tau*profile(loc_e + tau) ;
3445  if (tau < rounding_error) {
3446  qL = 0.0;
3447  } else {
3448  qhat = fncAvrQhat(loc_e, tau);
3449  // if (qhat*sqrt(2)>0.6)
3450  // {
3451  // JSINFO << MAGENTA << " Big q-hat warning ! ";
3452  // JSINFO << BOLDYELLOW << " length = " << length << " loc = " << loc_e << " tau = " << tau ;
3453  // JSINFO << BOLDYELLOW << " parton formed at x = " << initRx << " y = " << initRy << " z = " << initRz << " t = " << initR0 ;
3454  // JSINFO << BOLDYELLOW << " mean qhat for sudakov in GeV^2/fm = " << qhat*5*sqrt(2) ;
3455  // }
3456  qL = qhat * 2.0 * tau;
3457  }
3458 
3459  double e1 = M * M;
3460  double e2 = 1.0 / cg1;
3461  double e3 = e2 * e1;
3462  double e4 = e2 * cg;
3463  double e5 = std::log(e4);
3464  double e8 = std::log(1.0 - e4);
3465  double e13 = 1.0 / (e3 + 1.0) * e3;
3466  double e15 = std::log(1.0 - e4 - e13);
3467  double e18 = std::log(e4 + e13);
3468  double e22 = e2 * (-(2.0 * e5 - e8 + 1.0 - e4) * e3 +
3469  (2.0 * e15 - e18 + e4 + e13) * e3);
3470 
3471  double eL;
3472 
3473  if (tau < rounding_error) {
3474  eL = 0.0;
3475  } else {
3476  ehat = 0.0; //fncAvrEhat(loc_e,tau);
3477  if (ehat * sqrt(2) > 0.6) {
3478  // JSINFO << BOLDYELLOW << " length = " << length << " loc = " << loc_e << " tau = " << tau ;
3479  //JSINFO << BOLDYELLOW << " parton formed at x = " << initRx << " y = " << initRy << " z = " << initRz << " t = " << initR0 ;
3480  // JSINFO << BOLDYELLOW << " mean qhat for sudakov in GeV^2/fm = " << qhat*5*sqrt(2) ;
3481  }
3482  eL = ehat * 4.0;
3483  }
3484 
3485  double f1 = M * M;
3486  double f2 = 1.0 / cg1;
3487  double f3 = f2 * f1;
3488  double f4 = f2 * cg;
3489  double f5 = std::log(f4);
3490  double f8 = f1 * f1;
3491  double f9 = std::pow(cg1, 2.0);
3492  double f11 = 1.0 / f9 * f8;
3493  double f14 = 13.0 / 4.0 * f11 - 15.0 / 4.0 * f3 + 1.0 / 2.0;
3494  double f15 = 1.0 - f4;
3495  double f16 = std::log(f15);
3496  double f24 = f15 * f15;
3497  double f32 = 1.0 / (f3 + 1.0) * f3;
3498  double f33 = 1.0 - f4 - f32;
3499  double f34 = std::log(f33);
3500  double f37 = f4 + f32;
3501  double f38 = std::log(f37);
3502  double f45 = f37 * f37;
3503  double f52 = f2 * ((15.0 / 2.0 * f5 * f3 + f16 * f14 + 1.0 / cg * cg1 +
3504  15.0 / 4.0 * f2 / f15 * f1 - 13.0 / 8.0 / f24 * f11) *
3505  f3 -
3506  (15.0 / 2.0 * f34 * f3 + f38 * f14 + 1.0 / f33 +
3507  15.0 / 4.0 * f2 / f37 * f1 - 13.0 / 8.0 / f45 * f11) *
3508  f3);
3509  double e2L;
3510 
3511  if (tau < rounding_error) {
3512  e2L = 0.0;
3513  } else {
3514  e2hat = qhat / 2.0; //fncAvrE2hat(loc_e,tau);
3515  if (e2hat * sqrt(2) > 0.6) {
3516  // JSINFO << BOLDYELLOW << " length = " << length << " loc = " << loc_e << " tau = " << tau ;
3517  //JSINFO << BOLDYELLOW << " parton formed at x = " << initRx << " y = " << initRy << " z = " << initRz << " t = " << initR0 ;
3518  // JSINFO << BOLDYELLOW << " mean qhat for sudakov in GeV^2/fm = " << qhat*5*sqrt(2) ;
3519  }
3520  e2L = e2hat * 8.0 / (tau * cg1);
3521  }
3522 
3523  // JSINFO << BOLDRED << " qhat L = " << qL << " location = " << loc_e << " tau = " << tau << " length = " << length;
3524 
3525  res = t21 + qL * q53 / cg1;
3526 
3527  //+ eL*e22/cg1 +e2L*f52/cg1;
3528  // Uncomment only if you have an eL larger than 2 times e2L for charm, and derive expression for bottom.
3529  // MC simulation is not valid for all choices of e-hat and e2-hat.
3530 
3531  if (res < 0.0) {
3532  cerr << "ERROR: medium contribution negative in sud_z_QG : res = " << res
3533  << endl;
3534 
3535  throw std::runtime_error("ERROR: medium contribution negative in sud_z_QG");
3536  }
3537 
3538  return (res);
3539 }
3540 
3541 double Matter::P_z_qg_int(double cg, double cg1, double loc_e, double cg3,
3542  double l_fac, double E2) {
3543 
3544  double t2, t5, t7, t10, t12, q2, q6, q10, tau, qL, res;
3545 
3546  if ((cg < cg1 / (2.0 * E2 * E2 / cg1 + 1.0)))
3547  cg = cg1 / (2.0 * E2 * E2 / cg1 + 1.0);
3548 
3549  t2 = std::pow(cg1, 2);
3550  t5 = std::log(1.0 - cg1);
3551  t7 = std::pow(cg, 2);
3552  t10 = std::log(1.0 - cg);
3553  t12 = -cg1 - t2 / 2.0 - 2.0 * t5 + cg + t7 / 2.0 + 2.0 * t10;
3554 
3555  // return(t12);
3556 
3557  q2 = std::log(cg1);
3558  q6 = std::log(cg);
3559  q10 = q2 - 2.0 / (cg1 - 1.0) - q6 + 2.0 / (cg - 1.0);
3560 
3561  tau = l_fac;
3562 
3563  if ((length - loc_e) < tau)
3564  tau = (length - loc_e);
3565 
3566  if (loc_e > length)
3567  tau = 0.0;
3568 
3569  // SC
3570  //qL = qhat*0.6*tau*profile(loc_e + tau) ;
3571  if (tau < rounding_error) {
3572  qL = 0.0;
3573  } else {
3574  qhat = fncAvrQhat(loc_e, tau);
3575  qL = qhat * tau;
3576  }
3577 
3578  res = t12 + 2.0 * qL * q10 / cg3;
3579 
3580  return (res);
3581 }
3582 
3583 double Matter::P_z_qg_int_w_M(double M, double cg, double cg1, double loc_e,
3584  double cg3, double l_fac, double E2) {
3585 
3586  double t2, t5, t7, t10, t12, tau, qL, res;
3587 
3588  if (std::abs(cg - cg1) < rounding_error)
3589  return (cg);
3590  //if ((cg< cg1/(2.0*E2*E2/cg1+1.0) )) cg = cg1/( 2.0*E2*E2/cg1 + 1.0 );
3591 
3592  t2 = std::pow(cg1, 2);
3593  t5 = std::log(1.0 - cg1);
3594  t7 = std::pow(cg, 2);
3595  t10 = std::log(1.0 - cg);
3596  t12 = -cg1 - t2 / 2.0 - 2.0 * t5 + cg + t7 / 2.0 + 2.0 * t10;
3597 
3598  // return(t12);
3599 
3600  double q1 = M * M;
3601  double q2 = 1.0 / cg3;
3602  double q3 = q2 * q1;
3603  double q5 = 4.0 * q3 + 1.0;
3604  double q6 = 1.0 - cg1;
3605  double q7 = std::log(q6);
3606  double q9 = q1 * q1;
3607  double q10 = cg3 * cg3;
3608  double q12 = 1.0 / q10 * q9;
3609  double q14 = q12 - 2.0 * q3 + 1.0 / 2.0;
3610  double q15 = std::log(cg1);
3611  double q17 = q3 + 1.0;
3612  double q26 = std::pow(cg1, 2.0);
3613  double q30 = 1.0 - cg;
3614  double q31 = std::log(q30);
3615  double q33 = std::log(cg);
3616  double q43 = std::pow(cg, 2.0);
3617  double q47 = q7 * q5 + q15 * q14 + cg1 * q17 / 2.0 + 2.0 / q6 +
3618  3.0 / 2.0 * q2 / cg1 * q1 - 1.0 / q26 * q12 / 2.0 - q31 * q5 -
3619  q33 * q14 - cg * q17 / 2.0 - 2.0 / q30 -
3620  3.0 / 2.0 * q2 / cg * q1 + 1.0 / q43 * q12 / 2.0;
3621 
3622  tau = l_fac;
3623 
3624  if ((length - loc_e) < tau)
3625  tau = (length - loc_e);
3626 
3627  if (loc_e > length)
3628  tau = 0.0;
3629 
3630  // SC
3631  //qL = qhat*0.6*tau*profile(loc_e + tau) ;
3632  if (tau < rounding_error) {
3633  qL = 0.0;
3634  } else {
3635  qhat = fncAvrQhat(loc_e, tau);
3636  qL = qhat * 2.0 * tau;
3637  }
3638 
3639  double e1 = M * M;
3640  double e3 = 1.0 / cg3 * e1;
3641  double e4 = 1.0 - cg1;
3642  double e5 = std::log(e4);
3643  double e10 = 1.0 - cg;
3644  double e11 = std::log(e10);
3645  double e17 = 2.0 * (e5 + 1.0 / e4 + cg1 / 2.0) * e3 -
3646  2.0 * (e11 + 1.0 / e10 + cg / 2.0) * e3;
3647 
3648  double eL;
3649 
3650  if (tau < rounding_error) {
3651  eL = 0.0;
3652  } else {
3653  ehat = 0.0; //fncAvrEhat(loc_e,tau);
3654  eL = ehat * 4.0;
3655  }
3656 
3657  double f1 = M * M;
3658  double f2 = 1.0 / cg3;
3659  double f3 = f2 * f1;
3660  double f4 = 13.0 * f3;
3661  double f6 = f1 * f1;
3662  double f7 = f6 * (f4 + 15.0);
3663  double f8 = 1.0 - cg1;
3664  double f9 = std::log(f8);
3665  double f10 = cg3 * cg3;
3666  double f11 = 1.0 / f10;
3667  double f15 = std::log(cg1);
3668  double f23 = f1 * (39.0 / 4.0 * f11 * f6 + 15.0 / 2.0 * f3 + 1.0);
3669  double f28 = f6 * (f4 + 15.0 / 2.0);
3670  double f29 = f8 * f8;
3671  double f37 = 1.0 / f10 / cg3 * f6 * f1;
3672  double f42 = 1.0 - cg;
3673  double f43 = std::log(f42);
3674  double f47 = std::log(cg);
3675  double f54 = f42 * f42;
3676  double f63 = f11 * f9 * f7 / 4.0 + f2 * f1 * f15 / 2.0 + 1.0 / f8 * f2 * f23 -
3677  1.0 / f29 * f11 * f28 / 2.0 + 13.0 / 6.0 / f29 / f8 * f37 -
3678  f11 * f43 * f7 / 4.0 - f2 * f1 * f47 / 2.0 -
3679  1.0 / f42 * f2 * f23 + 1.0 / f54 * f11 * f28 / 2.0 -
3680  13.0 / 6.0 / f54 / f42 * f37;
3681 
3682  double e2L;
3683 
3684  if (tau < rounding_error) {
3685  e2L = 0.0;
3686  } else {
3687  e2hat = qhat / 2.0; //fncAvrE2hat(loc_e,tau);
3688  e2L = e2hat * 8.0 / (tau * cg3);
3689  }
3690 
3691  res = t12 + qL * q47 / cg3;
3692  //+ eL*e17/cg3 + e2L*f63/cg3;
3693  // Uncomment only if you have an eL larger than 2 times e2L for charm, and derive expression for bottom.
3694  // MC simulation is not valid for all choices of e-hat and e2-hat.
3695 
3696  return (res);
3697 }
3698 
3699 double Matter::alpha_s(double q2) {
3700  double a, L2, q24, c_nf;
3701 
3702  L2 = std::pow(Lambda_QCD, 2);
3703 
3704  q24 = q2 / 4.0;
3705 
3706  c_nf = nf;
3707 
3708  if (q24 > 4.0) {
3709  c_nf = 4;
3710  }
3711 
3712  if (q24 > 64.0) {
3713  c_nf = 5;
3714  }
3715 
3716  if (q24 > L2) {
3717  a = 12.0 * pi / (11.0 * Nc - 2.0 * c_nf) / std::log(q24 / L2);
3718  } else {
3719  JSWARN << " alpha too large ";
3720  a = 0.6;
3721  }
3722 
3723  return (a);
3724 }
3725 
3726 double Matter::profile(double zeta) {
3727  double prof;
3728 
3729  /* Modify the next set of lines to get a particular profile in brick test or further modify the profile for hydro to introduce coherence effects */
3730 
3731  prof = 1.0;
3732 
3733  return (prof);
3734 }
3735 
3737 
3738 double Matter::fillQhatTab(double y) {
3739 
3740  double xLoc, yLoc, zLoc, tLoc;
3741  double vxLoc, vyLoc, vzLoc, gammaLoc, betaLoc;
3742  double edLoc, sdLoc;
3743  double tempLoc;
3744  double flowFactor, qhatLoc;
3745  int hydro_ctl;
3746  double lastLength = initR0;
3747 
3748  double tStep = 0.1;
3749 
3750  std::unique_ptr<FluidCellInfo> check_fluid_info_ptr;
3751 
3752  for (int i = 0; i < dimQhatTab; i++) {
3753  tLoc = tStep * i;
3754 
3755  //if(tLoc<initR0-tStep) { // potential problem of making t^2<z^2
3756 
3757  double boostedTStart = tStart * std::cosh(y);
3758  if (tLoc < initR0 || tLoc < boostedTStart) {
3759  qhatTab1D[i] = 0.0;
3760  continue;
3761  }
3762 
3763  xLoc = initRx + (tLoc - initR0) * initVx;
3764  yLoc = initRy + (tLoc - initR0) * initVy;
3765  zLoc = initRz + (tLoc - initR0) * initVz;
3766 
3767  // if(bulkFlag == 1) { // read OSU hydro
3768  // readhydroinfoshanshan_(&tLoc,&xLoc,&yLoc,&zLoc,&edLoc,&sdLoc,&tempLoc,&vxLoc,&vyLoc,&vzLoc,&hydro_ctl);
3769  // } else if(bulkFlag == 2) { // read CCNU hydro
3770  // hydroinfoccnu_(&tLoc, &xLoc, &yLoc, &zLoc, &tempLoc, &vxLoc, &vyLoc, &vzLoc, &hydro_ctl);
3771  // } else if(bulkFlag == 0) { // static medium
3772  // vxLoc = 0.0;
3773  // vyLoc = 0.0;
3774  // vzLoc = 0.0;
3775  // hydro_ctl = 0;
3776  // tempLoc = T;
3777  // }
3778 
3779  if (std::isinf(tLoc) || std::isnan(tLoc) || std::isinf(zLoc) ||
3780  std::isnan(zLoc) || std::abs(zLoc) > tLoc) {
3781  JSWARN << "Third instance";
3782  JSWARN << "Loc for vector is:" << tLoc << ", " << xLoc << ", " << yLoc
3783  << ", " << zLoc;
3784  JSWARN << "initR0, initRx, initRy, initRz="
3785  << ", " << initR0 << ", " << initRx << ", " << initRy << ", "
3786  << initRz;
3787  JSWARN << "initVx, initVy, initVz =" << initVx << ", " << initVy << ", "
3788  << initVz;
3789  JSWARN << "initVMod=" << std::setprecision(20)
3790  << std::sqrt(initVx * initVx + initVy * initVy + initVz * initVz);
3791  JSWARN << "Can't dump pIn_info as we are in fillQhatTab. But it should "
3792  "be dumped right before this."; //Dump_pIn_info(i, pIn);
3793  //exit(0);
3794  }
3795 
3796  GetHydroCellSignal(tLoc, xLoc, yLoc, zLoc, check_fluid_info_ptr);
3797  VERBOSE(8) << MAGENTA << "Temperature from medium = "
3798  << check_fluid_info_ptr->temperature;
3799 
3800  tempLoc = check_fluid_info_ptr->temperature;
3801  sdLoc = check_fluid_info_ptr->entropy_density;
3802  vxLoc = check_fluid_info_ptr->vx;
3803  vyLoc = check_fluid_info_ptr->vy;
3804  vzLoc = check_fluid_info_ptr->vz;
3805 
3806  hydro_ctl = 0;
3807 
3808  if (hydro_ctl == 0 && tempLoc >= hydro_Tc) {
3809  lastLength = tLoc;
3810  betaLoc = sqrt(vxLoc * vxLoc + vyLoc * vyLoc + vzLoc * vzLoc);
3811  gammaLoc = 1.0 / sqrt(1.0 - betaLoc * betaLoc);
3812  flowFactor =
3813  gammaLoc * (1.0 - (initVx * vxLoc + initVy * vyLoc + initVz * vzLoc));
3814 
3815  //if(run_alphas==1){ alphas= 4*pi/(9.0*log(2*initEner*tempLoc/0.04));}
3816 
3817  /* if (qhat0 < 0.0) {
3818  // calculate qhat with alphas
3819  double muD2 = 6.0 * pi * alphas * tempLoc * tempLoc;
3820  // if(initEner > pi*tempLoc) qhatLoc = Ca*alphas*muD2*tempLoc*log(6.0*initEner*tempLoc/muD2);
3821  // else qhatLoc = Ca*alphas*muD2*tempLoc*log(6.0*pi*tempLoc*tempLoc/muD2);
3822  // fitted formula from https://arxiv.org/pdf/1503.03313.pdf
3823  if (initEner > pi * tempLoc)
3824  qhatLoc = Ca * 50.4864 / pi * pow(alphas, 2) * pow(tempLoc, 3) *
3825  log(5.7 * initEner * tempLoc / 4.0 / muD2);
3826  else
3827  qhatLoc = Ca * 50.4864 / pi * pow(alphas, 2) * pow(tempLoc, 3) *
3828  log(5.7 * pi * tempLoc * tempLoc / 4.0 / muD2);
3829  qhatLoc = qhatLoc * flowFactor;
3830  if (qhatLoc < 0.0)
3831  qhatLoc = 0.0;
3832  } else { // use input qhat
3833  if (brick_med) {
3834  qhatLoc = qhat0 * 0.1973 * flowFactor;
3835  } else {
3836  qhatLoc = qhat0 / 96.0 * sdLoc * 0.1973 *
3837  flowFactor; // qhat0 at s = 96fm^-3
3838  }
3839  }*/
3840 
3841  // GeneralQhatFunction(int QhatParametrizationType, double Temperature, double EntropyDensity, double FixAlphas, double Qhat0, double E, double muSquare);
3842  double muSquare=-1;//For virtuality dependent cases, we explicitly modify q-hat inside Sudakov, due to which we set here scale=-1; Alternatively one could extend the dimension of the q-hat table
3843 
3844  qhatLoc= GeneralQhatFunction(QhatParametrizationType, tempLoc, sdLoc, alphas, qhat0, initEner, muSquare);
3845  qhatLoc = qhatLoc * flowFactor;
3846 
3847  //JSINFO << "check qhat -- ener, T, qhat: " << initEner << " , " << tempLoc << " , " << qhatLoc;
3848  } else { // outside the QGP medium
3849  qhatLoc = 0.0;
3850  }
3851 
3852  qhatTab1D[i] =
3853  qhatLoc / sqrt(2.0); // store qhat value in light cone coordinate
3854  }
3855 
3856  for (int i = 0; i < dimQhatTab; i++) { // dim of loc
3857 
3858  double totValue = 0.0;
3859 
3860  for (int j = 0; i + j < dimQhatTab; j++) { // dim of tau_f
3861 
3862  totValue = totValue + qhatTab1D[i + j];
3863  qhatTab2D[i][j] = totValue / (j + 1);
3864  }
3865  }
3866 
3867  //return(lastLength*sqrt(2.0)*5.0); // light cone + GeV unit
3868  return ((2.0 * lastLength + initRdotV - initR0) / sqrt(2.0) *
3869  5.0); // light cone + GeV unit
3870 }
3871 
3873 // E is the energy and muSquare is the virtuality of the parton
3874 double Matter::GeneralQhatFunction(int QhatParametrization, double Temperature, double EntropyDensity, double FixAlphas, double Qhat0, double E, double muSquare)
3875 {
3876  int ActiveFlavor=3; qhat=0.0;
3877  double DebyeMassSquare = FixAlphas*4*pi*pow(Temperature,2.0)*(6.0 + ActiveFlavor)/6.0;
3878  double ScaleNet=2*E*Temperature;
3879  if(ScaleNet < 1.0){ ScaleNet=1.0; }
3880  switch(QhatParametrization)
3881  {
3882  //HTL formula with all alpha_s as constant and controlled by XML
3883  case 0:
3884  qhat = (Ca*50.4864/pi)*pow(FixAlphas,2)*pow(Temperature,3)*log(ScaleNet/DebyeMassSquare);
3885  break;
3886 
3887  //alpha_s at scale muS=2ET and second alpha_s at muS=DebyeMassSquare is fit parameter
3888  case 1:
3889  qhat = (Ca*50.4864/pi)*RunningAlphaS(ScaleNet)*FixAlphas*pow(Temperature,3)*log(ScaleNet/DebyeMassSquare);
3890  break;
3891 
3892  //Constant q-hat
3893  case 2:
3894  qhat = Qhat0*0.1973;
3895  break;
3896 
3897  //Scale with T^3
3898  case 3:
3899  qhat = Qhat0*pow(Temperature/0.3,3)*0.1973; // w.r.t T=0.3 GeV
3900  break;
3901 
3902  //Scale with entropy density
3903  case 4:
3904  qhat = Qhat0*(EntropyDensity/96.0)*0.1973; // w.r.t S0=96 fm^-3
3905  break;
3906 
3907  //HTL q-hat multiplied by Virtuality dependent function to mimic PDF-Scale dependent q-hat
3908  //Function is 1/(1+A*pow(log(Q^2),2)+B*pow(log(Q^2),4))
3909  case 5:
3910  qhat = (Ca*50.4864/pi)*RunningAlphaS(ScaleNet)*FixAlphas*pow(Temperature,3)*log(ScaleNet/DebyeMassSquare);
3911  qhat = qhat*VirtualityQhatFunction(5, E, muSquare);
3912  break;
3913 
3914  //HTL q-hat multiplied by Virtuality dependent function to mimic PDF-Scale dependent q-hat
3915  //Function is int^{1}_{xB} e^{-ax} / (1+A*pow(log(Q^2),1)+B*pow(log(Q^2),2))
3916  case 6:
3917  qhat = (Ca*50.4864/pi)*RunningAlphaS(ScaleNet)*FixAlphas*pow(Temperature,3)*log(ScaleNet/DebyeMassSquare);
3918  qhat = qhat*VirtualityQhatFunction(6, E, muSquare);
3919  break;
3920 
3921  //HTL q-hat multiplied by Virtuality dependent function to mimic PDF-Scale dependent q-hat
3922  //Function is int^{1}_{xB} x^{a}(1-x)^{b} / (1+A*pow(log(Q^2),1)+B*pow(log(Q^2),2))
3923  case 7:
3924  qhat = (Ca*50.4864/pi)*RunningAlphaS(ScaleNet)*FixAlphas*pow(Temperature,3)*log(ScaleNet/DebyeMassSquare);
3925  qhat = qhat*VirtualityQhatFunction(7, E, muSquare);
3926  break;
3927 
3928  default:
3929  JSINFO<<"q-hat Parametrization "<<QhatParametrization<<" is not used, qhat will be set to zero";
3930  }
3931  return qhat;
3932 }
3933 
3935 double Matter::RunningAlphaS(double muSquare)
3936 {
3937  int ActiveFlavor=3;
3938  double Square_Lambda_QCD_HTL = exp( -12.0*pi/( (33 - 2*ActiveFlavor)*alphas) );
3939  double ans = 12.0*pi/( (33.0- 2.0*ActiveFlavor)*log(muSquare/Square_Lambda_QCD_HTL) );
3940  if(muSquare < 1.0) {ans=alphas; }
3941 
3942  VERBOSE(8)<<"Fixed-alphaS="<<alphas<<", Lambda_QCD_HTL="<<sqrt(Square_Lambda_QCD_HTL)<<", mu2="<<muSquare<<", Running alpha_s"<<ans;
3943  return ans;
3944 }
3947 double Matter::VirtualityQhatFunction(int QhatParametrization, double enerLoc, double muSquare)
3948 {
3949  double ans=0;double xB=0, xB0=0, IntegralNorm=0;
3950 
3951  if( muSquare <= Q00*Q00) {ans=1;}
3952  else
3953  {
3954  switch(QhatParametrization)
3955  {
3956  case 0:
3957  ans = 1.0;
3958  break;
3959 
3960  case 1:
3961  ans = 1.0;
3962  break;
3963 
3964  case 2:
3965  ans = 1.0;
3966  break;
3967 
3968  case 3:
3969  ans = 1.0;
3970  break;
3971 
3972  case 4:
3973  ans = 1.0;
3974  break;
3975 
3976  case 5:
3977  ans = 1.0 + qhatA*log(Q00*Q00)*log(Q00*Q00) + qhatB*pow(log(Q00*Q00),4);
3978  ans = ans/( 1.0 + qhatA*log(muSquare)*log(muSquare) + qhatB*pow(log(muSquare),4) );
3979  break;
3980 
3981  case 6:
3982  xB = muSquare/(2.0*enerLoc);
3983  xB0 = Q00*Q00/(2.0*enerLoc); if(xB<=xB0){ans=1.0; break; }
3984  if ( qhatC > 0.0 && xB < 0.99)
3985  {
3986  ans = ( exp(qhatC*(1.0-xB)) - 1.0 )/(1.0 + qhatA*log(muSquare/0.04) + qhatB*log(muSquare/0.04)*log(muSquare/0.04) );
3987  ans = ans*(1.0 + qhatA*log(Q00*Q00/0.04) + qhatB*log(Q00*Q00/0.04)*log(Q00*Q00/0.04) )/( exp(qhatC*(1.0-xB0)) - 1.0 );
3988  //JSINFO<<"K xB="<<xB<<", and (E,muSquare)=("<<enerLoc<<","<<muSquare<<"), and Virtuality dep Qhat="<<ans;
3989  }
3990  else if( qhatC == 0.0 && xB < 0.99)
3991  {
3992  ans = (1.0-xB)/(1.0 + qhatA*log(muSquare/0.04) + qhatB*log(muSquare/0.04)*log(muSquare/0.04) );
3993  ans = ans*(1.0 + qhatA*log(Q00*Q00/0.04) + qhatB*log(Q00*Q00/0.04)*log(Q00*Q00/0.04) )/(1-xB0);
3994  }
3995  else {ans=0.0;}
3996  //JSINFO<<"L xB="<<xB<<", and (E,muSquare)=("<<enerLoc<<","<<muSquare<<"), and Virtuality dep Qhat="<<ans;
3997  break;
3998 
3999  case 7:
4000  xB = muSquare/(2.0*enerLoc);
4001  xB0 = Q00*Q00/(2.0*enerLoc); if(xB<=xB0){ans=1.0; break; }
4002  ans = IntegralPDF(xB, qhatC, qhatD)/(1.0 + qhatA*log(muSquare/0.04) + qhatB*log(muSquare/0.04)*log(muSquare/0.04) );
4003  IntegralNorm = IntegralPDF(xB0, qhatC, qhatD);
4004  if( IntegralNorm > 0.0 )
4005  {
4006  ans=ans*(1.0 + qhatA*log(muSquare/0.04) + qhatB*log(muSquare/0.04)*log(muSquare/0.04) )/IntegralNorm;
4007  }
4008  else {ans=0;}
4009  //JSINFO<<"L xB="<<xB<<", and (E,muSquare)=("<<enerLoc<<","<<muSquare<<"), and Virtuality dep Qhat="<<ans;
4010  break;
4011 
4012  default:
4013  JSINFO<<"q-hat Parametrization "<<QhatParametrization<<" is not used, VirtualityQhatFunction is set to zero or one";
4014  }
4015  }
4016 
4017  //JSINFO<<"Qhat Type="<<QhatParametrization<<", and (E,muSquare)=("<<enerLoc<<","<<muSquare<<"), and Virtuality dep Qhat="<<ans;
4018  return ans;
4019 }
4021 double Matter::ModifiedProbability(int QhatParametrization, double tempLoc, double sdLoc, double enerLoc, double muSquare)
4022 {
4023  double ModifiedAlphas=0;double qhatLoc=0;
4024  double ScaleNet=2*enerLoc*tempLoc;
4025  if(ScaleNet <1.0) { ScaleNet=1.0; }
4026 
4027  switch(QhatParametrization)
4028  {
4029  //For HTL q-hat formula with all alpha_s as constant and controlled by XML
4030  case 0:
4031  ModifiedAlphas = alphas;
4032  break;
4033 
4034  //For HTL q-hat where one alpha_s at scale muS=2ET and second alpha_s at muS=DebyeMassSquare is fit parameter
4035  case 1:
4036  ModifiedAlphas = RunningAlphaS(ScaleNet);
4037  break;
4038 
4039  //Constant q-hat case, alphas is computed using HTL formula with all alpha_s fixed
4040  case 2:
4041  qhatLoc = qhat0*0.1973;
4042  ModifiedAlphas = solve_alphas(qhatLoc, enerLoc, tempLoc);
4043  break;
4044 
4045  //For q-hat goes as T^3, alphas is computed using HTL formula with all alpha_s fixed
4046  case 3:
4047  qhatLoc = qhat0*pow(tempLoc/0.3,3)*0.1973; // w.r.t T=0.3 GeV
4048  ModifiedAlphas = solve_alphas(qhatLoc, enerLoc, tempLoc);
4049  break;
4050 
4051  //For q-hat goes as entropy density
4052  case 4:
4053  qhatLoc = qhat0*(sdLoc/96.0)*0.1973; // w.r.t S0=96 fm^-3
4054  ModifiedAlphas = solve_alphas(qhatLoc, enerLoc, tempLoc);
4055  break;
4056 
4057  //For HTL q-hat multiplied by Virtuality dependent function to mimic PDF-Scale dependent q-hat
4058  //Function is 1 / (1+A*pow(log(Q^2),2)+B*pow(log(Q^2),4))
4059  case 5:
4060  ModifiedAlphas = RunningAlphaS(ScaleNet)*VirtualityQhatFunction(5, enerLoc, muSquare) ;
4061  break;
4062  //HTL q-hat multiplied by Virtuality dependent function to mimic PDF-Scale dependent q-hat
4063  //Function is int^{1}_{xB} e^{-ax} / (1+A*pow(log(Q^2),1)+B*pow(log(Q^2),2))
4064  case 6:
4065  ModifiedAlphas = RunningAlphaS(ScaleNet)*VirtualityQhatFunction(6, enerLoc, muSquare) ;
4066  break;
4067 
4068  //HTL q-hat multiplied by Virtuality dependent function to mimic PDF-Scale dependent q-hat
4069  //Function is int^{1}_{xB} x^{a}(1-x)^{b} / (1+A*pow(log(Q^2),1)+B*pow(log(Q^2),2))
4070  case 7:
4071  ModifiedAlphas = RunningAlphaS(ScaleNet)*VirtualityQhatFunction(7, enerLoc, muSquare) ;
4072  break;
4073 
4074  default:
4075  JSINFO<<"q-hat Parametrization "<<QhatParametrization<<" is not used,Elastic scattering alphas will be set to zero";
4076  }
4077  //JSINFO<<"q-hat Parametrization "<<QhatParametrization<<" modified alphas="<<ModifiedAlphas;
4078  return ModifiedAlphas;
4079 }
4080 
4081 //x integration of QGP-PDF from xB to 1
4082 double Matter::IntegralPDF(double xB, double a, double b)
4083 {
4084  double Xmin=0.01, Xmax=0.99, X=0, dX=0, ans=0; int N=100, ix=0;
4085  dX = (1.0-0.0)/N;
4086  if(xB > Xmax) {ans=0;}
4087  else
4088  {
4089  ix = xB/dX;
4090  for(int i=ix; i<N; i++)
4091  {
4092  X= (i+0.5)*dX;
4093  if (X<Xmin){X=Xmin;}
4094  ans = ans + pow(X,a)*pow(1-X,b)*dX;
4095  }
4096  }
4097  //JSINFO<<"(xB,a,b)=("<<xB<<","<<a<<","<<b<<"), \t Area="<<ans;
4098  return ans;
4099 }
4100 
4101 
4102 double Matter::fncQhat(double zeta) {
4103  if (in_vac)
4104  return (0.0);
4105 
4106  double tStep = 0.1;
4107  //int indexZeta = (int)(zeta/sqrt(2.0)/5.0/tStep+0.5); // zeta was in 1/GeV and light cone coordinate
4108  int indexZeta =
4109  (int)((sqrt(2.0) * zeta / 5.0 - initRdotV + initR0) / 2.0 / tStep +
4110  0.5); // zeta was in 1/GeV and light cone coordinate
4111 
4112  //if(indexZeta >= dimQhatTab) indexZeta = dimQhatTab-1;
4113  if (indexZeta >= dimQhatTab)
4114  return (0);
4115 
4116  double avrQhat = qhatTab1D[indexZeta]*VirtualityQhatFunction(QhatParametrizationType, initEner, tscale);
4117  return (avrQhat);
4118 }
4119 
4121 
4122 double Matter::fncAvrQhat(double zeta, double tau) {
4123 
4124  if (in_vac)
4125  return (0.0);
4126 
4127  double tStep = 0.1;
4128  //int indexZeta = (int)(zeta/sqrt(2.0)/5.0/tStep+0.5); // zeta was in 1/GeV and light cone coordinate
4129  int indexZeta =
4130  (int)((sqrt(2.0) * zeta / 5.0 - initRdotV + initR0) / 2.0 / tStep +
4131  0.5); // zeta was in 1/GeV and light cone coordinate
4132  int indexTau = (int)(tau / sqrt(2.0) / 5.0 / tStep +
4133  0.5); // tau was in 1/GeV and light cone coordinate
4134 
4135  // if(indexZeta >= dimQhatTab) indexZeta = dimQhatTab-1;
4136  if (indexZeta >= dimQhatTab)
4137  return (0);
4138  if (indexTau >= dimQhatTab)
4139  indexTau = dimQhatTab - 1;
4140 
4141  double avrQhat = qhatTab2D[indexZeta][indexTau]*VirtualityQhatFunction(QhatParametrizationType, initEner, tscale);
4142  return (avrQhat);
4143 }
4144 
4146 
4147 void Matter::flavor(int &CT, int &KATT0, int &KATT2, int &KATT3,
4148  unsigned int &max_color, unsigned int &color0,
4149  unsigned int &anti_color0, unsigned int &color2,
4150  unsigned int &anti_color2, unsigned int &color3,
4151  unsigned int &anti_color3) {
4152 
4153  int vb[7] = {0};
4154  int b = 0;
4155  int KATT00 = KATT0;
4156  unsigned int backup_color0 = color0;
4157  unsigned int backup_anti_color0 = anti_color0;
4158 
4159  vb[1] = 1;
4160  vb[2] = 2;
4161  vb[3] = 3;
4162  vb[4] = -1;
4163  vb[5] = -2;
4164  vb[6] = -3;
4165 
4166  if (KATT00 == 21) { //.....for gluon
4167  double R1 = 16.0; // gg->gg DOF_g
4168  double R2 = 0.0; // gg->qqbar don't consider this channel in MATTER
4169  double R3 = 6.0 * 6 * 4 / 9; // gq->gq or gqbar->gqbar flavor*DOF_q*factor
4170  double R0 = R1 + R3;
4171 
4172  double a = ran0(&NUM1);
4173 
4174  if (a <= R1 / R0) { // gg->gg
4175  CT = 1;
4176  KATT3 = 21;
4177  KATT2 = 21;
4178  //KATT0=KATT0;
4179  max_color++;
4180  color0 = backup_color0;
4181  anti_color0 = max_color;
4182  max_color++;
4183  color2 = anti_color0;
4184  anti_color2 = max_color;
4185  color3 = backup_anti_color0;
4186  anti_color3 = max_color;
4187  } else { // gq->gq
4188  CT = 3;
4189  b = floor(ran0(&NUM1) * 6 + 1);
4190  if (b == 7)
4191  b = 6;
4192  KATT3 = vb[b];
4193  KATT2 = KATT3;
4194  //KATT0=KATT0;
4195  if (KATT3 > 0) { // gq->gq
4196  max_color++;
4197  color0 = backup_color0;
4198  anti_color0 = max_color;
4199  color2 = max_color;
4200  anti_color2 = 0;
4201  color3 = backup_anti_color0;
4202  anti_color3 = 0;
4203  } else { // gqbar->gqbar
4204  max_color++;
4205  color0 = max_color;
4206  anti_color0 = backup_anti_color0;
4207  color2 = 0;
4208  anti_color2 = max_color;
4209  color3 = 0;
4210  anti_color3 = backup_color0;
4211  }
4212  }
4213  } else if (abs(KATT00) == 4) { // for charm quarks
4214  double R1 = 6.0 * 6 * 4 / 9; // Qq->Qq
4215  double R2 = 16.0; // Qg->Qg DOF_ag
4216  double R00 = R1 + R2;
4217 
4218  double a = ran0(&NUM1);
4219 
4220  if (a <= R2 / R00) { // Qg->Qg
4221  CT = 12;
4222  KATT3 = 21;
4223  KATT2 = 21;
4224  if (KATT00 > 0) { // Qg->Qg
4225  max_color++;
4226  color0 = max_color;
4227  anti_color0 = 0;
4228  max_color++;
4229  color2 = max_color;
4230  anti_color2 = color0;
4231  color3 = max_color;
4232  anti_color3 = backup_color0;
4233  } else { // Qbarg->Qbarg
4234  max_color++;
4235  color0 = 0;
4236  anti_color0 = max_color;
4237  max_color++;
4238  color2 = anti_color0;
4239  anti_color2 = max_color;
4240  color3 = backup_anti_color0;
4241  anti_color3 = max_color;
4242  }
4243  } else { // Qq->Qq
4244  CT = 11;
4245  b = floor(ran0(&NUM1) * 6 + 1);
4246  if (b == 7)
4247  b = 6;
4248  KATT3 = vb[b];
4249  KATT2 = KATT3;
4250  if (KATT00 > 0 && KATT2 > 0) { // qq->qq
4251  max_color++;
4252  color0 = max_color;
4253  anti_color0 = 0;
4254  color2 = backup_color0;
4255  anti_color2 = 0;
4256  color3 = max_color;
4257  anti_color3 = 0;
4258  } else if (KATT00 > 0 && KATT2 < 0) { //qqbar->qqbar
4259  max_color++;
4260  color0 = max_color;
4261  anti_color0 = 0;
4262  color2 = 0;
4263  anti_color2 = max_color;
4264  color3 = 0;
4265  anti_color3 = backup_color0;
4266  } else if (KATT00 < 0 && KATT2 > 0) { //qbarq->qbarq
4267  max_color++;
4268  color0 = 0;
4269  anti_color0 = max_color;
4270  color2 = max_color;
4271  anti_color2 = 0;
4272  color3 = backup_anti_color0;
4273  anti_color3 = 0;
4274  } else { //qbarqbar->qbarqbar
4275  max_color++;
4276  color0 = 0;
4277  anti_color0 = max_color;
4278  color2 = 0;
4279  anti_color2 = backup_anti_color0;
4280  color3 = 0;
4281  anti_color3 = max_color;
4282  }
4283  }
4284 
4285  } else { //.....for quark and antiquark (light)
4286  double R3 = 16.0; // qg->qg DOF_g
4287  double R4 = 4.0 * 6 * 4 / 9; // qq'->qq' scatter with other species
4288  double R5 = 1.0 * 6 * 4 / 9; // qq->qq scatter with itself
4289  double R6 =
4290  0.0; // qqbar->q'qbar' to other final state species, don't consider in MATTER
4291  double R7 = 1.0 * 6 * 4 / 9; // qqbar->qqbar scatter with its anti-particle
4292  double R8 = 0.0; // qqbar->gg don't consider in MATTER
4293  double R00 = R3 + R4 + R5 + R7;
4294 
4295  double a = ran0(&NUM1);
4296  if (a <= R3 / R00) { // qg->qg
4297  CT = 13;
4298  KATT3 = 21;
4299  KATT2 = 21;
4300  //KATT0=KATT0;
4301  if (KATT00 > 0) { // qg->qg
4302  max_color++;
4303  color0 = max_color;
4304  anti_color0 = 0;
4305  max_color++;
4306  color2 = max_color;
4307  anti_color2 = color0;
4308  color3 = max_color;
4309  anti_color3 = backup_color0;
4310  } else { // qbarg->qbarg
4311  max_color++;
4312  color0 = 0;
4313  anti_color0 = max_color;
4314  max_color++;
4315  color2 = anti_color0;
4316  anti_color2 = max_color;
4317  color3 = backup_anti_color0;
4318  anti_color3 = max_color;
4319  }
4320  } else if (a <= (R3 + R4) / R00) { // qq'->qq'
4321  CT = 4;
4322  do {
4323  b = floor(ran0(&NUM1) * 6 + 1);
4324  if (b == 7)
4325  b = 6;
4326  KATT3 = vb[b];
4327  } while (KATT3 == KATT0 || KATT3 == -KATT0);
4328  KATT2 = KATT3;
4329  //KATT0=KATT0;
4330  if (KATT00 > 0 && KATT2 > 0) { // qq->qq
4331  max_color++;
4332  color0 = max_color;
4333  anti_color0 = 0;
4334  color2 = backup_color0;
4335  anti_color2 = 0;
4336  color3 = max_color;
4337  anti_color3 = 0;
4338  } else if (KATT00 > 0 && KATT2 < 0) { //qqbar->qqbar
4339  max_color++;
4340  color0 = max_color;
4341  anti_color0 = 0;
4342  color2 = 0;
4343  anti_color2 = max_color;
4344  color3 = 0;
4345  anti_color3 = backup_color0;
4346  } else if (KATT00 < 0 && KATT2 > 0) { //qbarq->qbarq
4347  max_color++;
4348  color0 = 0;
4349  anti_color0 = max_color;
4350  color2 = max_color;
4351  anti_color2 = 0;
4352  color3 = backup_anti_color0;
4353  anti_color3 = 0;
4354  } else { //qbarqbar->qbarqbar
4355  max_color++;
4356  color0 = 0;
4357  anti_color0 = max_color;
4358  color2 = 0;
4359  anti_color2 = backup_anti_color0;
4360  color3 = 0;
4361  anti_color3 = max_color;
4362  }
4363  } else if (a <= (R3 + R4 + R5) / R00) { // scatter with itself
4364  CT = 5;
4365  KATT3 = KATT0;
4366  KATT2 = KATT0;
4367  //KATT0=KATT0;
4368  if (KATT00 > 0) { // qq->qq
4369  max_color++;
4370  color0 = max_color;
4371  anti_color0 = 0;
4372  color2 = backup_color0;
4373  anti_color2 = 0;
4374  color3 = max_color;
4375  anti_color3 = 0;
4376  } else { //qbarqbar->qbarqbar
4377  max_color++;
4378  color0 = 0;
4379  anti_color0 = max_color;
4380  color2 = 0;
4381  anti_color2 = backup_anti_color0;
4382  color3 = 0;
4383  anti_color3 = max_color;
4384  }
4385  } else { // scatter with its anti-particle
4386  CT = 7;
4387  KATT3 = -KATT0;
4388  KATT2 = KATT3;
4389  //KATT0=KATT0;
4390  if (KATT00 > 0) { //qqbar->qqbar
4391  max_color++;
4392  color0 = max_color;
4393  anti_color0 = 0;
4394  color2 = 0;
4395  anti_color2 = max_color;
4396  color3 = 0;
4397  anti_color3 = backup_color0;
4398  } else { //qbarq->qbarq
4399  max_color++;
4400  color0 = 0;
4401  anti_color0 = max_color;
4402  color2 = max_color;
4403  anti_color2 = 0;
4404  color3 = backup_anti_color0;
4405  anti_color3 = 0;
4406  }
4407  }
4408  }
4409 }
4410 
4411 void Matter::colljet22(int CT, double temp, double qhat0ud, double v0[4],
4412  double p0[4], double p2[4], double p3[4], double p4[4],
4413  double &qt) {
4414  //
4415  // p0 initial jet momentum, output to final momentum
4416  // p2 final thermal momentum,p3 initial termal energy
4417  //
4418  // amss=sqrt(abs(p0(4)**2-p0(1)**2-p0(2)**2-p0(3)**2))
4419  //
4420  //************************************************************
4421  p4[1] = p0[1];
4422  p4[2] = p0[2];
4423  p4[3] = p0[3];
4424  p4[0] = p0[0];
4425  //************************************************************
4426 
4427  // transform to local comoving frame of the fluid
4428  // cout << endl;
4429  // cout << "flow "<< v0[1] << " " << v0[2] << " " << v0[3] << " "<<" Elab " << p0[0] << endl;
4430 
4431  trans(v0, p0);
4432  // cout << p0[0] << " " << sqrt(qhat0ud) << endl;
4433 
4434  // cout << sqrt(pow(p0[1],2)+pow(p0[2],2)+pow(p0[3],2)) << " " << p0[1] << " " << p0[2] << " " << p0[3] << endl;
4435 
4436  //************************************************************
4437  trans(v0, p4);
4438  //************************************************************
4439 
4440  // sample the medium parton thermal momentum in the comoving frame
4441 
4442  double xw;
4443  double razim;
4444  double rcos;
4445  double rsin;
4446 
4447  double ss;
4448  double tmin;
4449  double tmid;
4450  double tmax;
4451 
4452  double rant;
4453  double tt;
4454 
4455  double uu;
4456  double ff = 0.0;
4457  double rank;
4458 
4459  double mmax;
4460  double msq = 0.0;
4461 
4462  double f1;
4463  double f2;
4464 
4465  double p0ex[4] = {0.0};
4466  double vc[4] = {0.0};
4467 
4468  int ct1_loop, ct2_loop, flag1, flag2;
4469 
4470  flag1 = 0;
4471  flag2 = 0;
4472 
4473  // Initial 4-momentum of jet
4474  //
4475  //************************************************************
4476  p4[1] = p0[1];
4477  p4[2] = p0[2];
4478  p4[3] = p0[3];
4479  p4[0] = p0[0];
4480  //************************************************************
4481 
4482  int ic = 0;
4483 
4484  ct1_loop = 0;
4485  do {
4486  ct1_loop++;
4487  if(flag2 == 1 || ct1_loop > 1e6){
4488  flag1 = 1;
4489  break;
4490  }
4491  ct2_loop = 0;
4492  do {
4493  ct2_loop++;
4494  if(ct2_loop > 1e6){
4495  flag2 = 1;
4496  break;
4497  }
4498  xw = 15.0 * ran0(&NUM1);
4499  razim = 2.0 * pi * ran0(&NUM1);
4500  rcos = 1.0 - 2.0 * ran0(&NUM1);
4501  rsin = sqrt(1.0 - rcos * rcos);
4502  //
4503  p2[0] = xw * temp;
4504  p2[3] = p2[0] * rcos;
4505  p2[1] = p2[0] * rsin * cos(razim);
4506  p2[2] = p2[0] * rsin * sin(razim);
4507 
4508  //
4509  // cms energy
4510  //
4511  ss =
4512  2.0 * (p0[0] * p2[0] - p0[1] * p2[1] - p0[2] * p2[2] - p0[3] * p2[3]);
4513 
4514  // if(ss.lt.2.d0*qhat0ud) goto 14
4515 
4516  tmin = qhat0ud;
4517  tmid = ss / 2.0;
4518  tmax = ss - qhat0ud;
4519 
4520  // use (s^2+u^2)/(t+qhat0ud)^2 as scattering cross section in the
4521  //
4522  rant = ran0(&NUM1);
4523  tt = rant * ss;
4524 
4525  // ic+=1;
4526  // cout << p0[0] << " " << p2[0] << endl;
4527  // cout << tt << " " << ss << "" << qhat0ud <<endl;
4528  // cout << ic << endl;
4529 
4530  } while ((tt < qhat0ud) || (tt > (ss - qhat0ud)));
4531 
4532  f1 = pow(xw, 3) / (exp(xw) - 1) / 1.4215;
4533  f2 = pow(xw, 3) / (exp(xw) + 1) / 1.2845;
4534 
4535  uu = ss - tt;
4536 
4537  if (CT == 1) {
4538  ff = f1;
4539  mmax =
4540  4.0 / pow(ss, 2) *
4541  (3.0 - tmin * (ss - tmin) / pow(ss, 2) +
4542  (ss - tmin) * ss / pow(tmin, 2) + tmin * ss / pow((ss - tmin), 2));
4543  msq = pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
4544  (3.0 - tt * uu / pow(ss, 2) + uu * ss / pow(tt, 2) +
4545  tt * ss / pow(uu, 2)) /
4546  mmax;
4547  }
4548 
4549  if (CT == 2) {
4550  ff = f1;
4551  mmax = 4.0 / pow(ss, 2) *
4552  (4.0 / 9.0 * (pow(tmin, 2) + pow((ss - tmin), 2)) / tmin /
4553  (ss - tmin) -
4554  (pow(tmin, 2) + pow((ss - tmin), 2)) / pow(ss, 2));
4555  msq = pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
4556  (4.0 / 9.0 * (pow(tt, 2) + pow(uu, 2)) / tt / uu -
4557  (pow(tt, 2) + pow(uu, 2)) / pow(ss, 2)) /
4558  (mmax + 4.0);
4559  }
4560 
4561  if (CT == 3) {
4562  ff = f2;
4563  if (((pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
4564  4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / ss / (ss - tmin)) >
4565  ((pow(ss, 2) + pow((ss - tmax), 2) / pow(tmax, 2) +
4566  4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmax), 2)) / ss /
4567  (ss - tmax)))) {
4568  mmax =
4569  4.0 / pow(ss, 2) *
4570  ((pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
4571  4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / ss / (ss - tmin));
4572  } else {
4573  mmax =
4574  4.0 / pow(ss, 2) *
4575  ((pow(ss, 2) + pow((ss - tmax), 2)) / pow(tmax, 2) +
4576  4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmax), 2)) / ss / (ss - tmax));
4577  }
4578  //
4579  msq = pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
4580  ((pow(ss, 2) + pow(uu, 2)) / pow(tt, 2) +
4581  4.0 / 9.0 * (pow(ss, 2) + pow(uu, 2)) / ss / uu) /
4582  mmax;
4583  }
4584 
4585  if (CT == 13) {
4586  ff = f1;
4587 
4588  if (((pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
4589  4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / ss / (ss - tmin)) >
4590  ((pow(ss, 2) + pow((ss - tmax), 2) / pow(tmax, 2) +
4591  4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmax), 2)) / ss /
4592  (ss - tmax)))) {
4593  mmax =
4594  4.0 / pow(ss, 2) *
4595  ((pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
4596  4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / ss / (ss - tmin));
4597  } else {
4598  mmax =
4599  4.0 / pow(ss, 2) *
4600  ((pow(ss, 2) + pow((ss - tmax), 2)) / pow(tmax, 2) +
4601  4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmax), 2)) / ss / (ss - tmax));
4602  }
4603  //
4604  msq = pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
4605  ((pow(ss, 2) + pow(uu, 2)) / pow(tt, 2) +
4606  4.0 / 9.0 * (pow(ss, 2) + pow(uu, 2)) / ss / uu) /
4607  mmax;
4608  }
4609 
4610  if (CT == 4) {
4611  ff = f2;
4612  mmax = 4.0 / pow(ss, 2) *
4613  (4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2));
4614  msq = pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
4615  (4.0 / 9.0 * (pow(ss, 2) + pow(uu, 2)) / pow(tt, 2)) / mmax;
4616  }
4617 
4618  if (CT == 5) {
4619  ff = f2;
4620  mmax = 4.0 / pow(ss, 2) *
4621  (4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
4622  (pow(ss, 2) + pow(tmin, 2)) / pow((ss - tmin), 2) -
4623  2.0 / 3.0 * pow(ss, 2) / tmin / (ss - tmin));
4624  msq = pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
4625  (4.0 / 9.0 *
4626  ((pow(ss, 2) + pow(uu, 2)) / pow(tt, 2) +
4627  (pow(ss, 2) + pow(tt, 2)) / pow(uu, 2) -
4628  2.0 / 3.0 * pow(ss, 2) / tt / uu)) /
4629  mmax;
4630  }
4631 
4632  if (CT == 6) {
4633  ff = f2;
4634  mmax = 4.0 / pow(ss, 2) *
4635  (4.0 / 9.0 * (pow(tmin, 2) + pow((ss - tmin), 2)) / pow(ss, 2));
4636  msq = pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
4637  (4.0 / 9.0 * (pow(tt, 2) + pow(uu, 2)) / pow(ss, 2)) / (mmax + 0.5);
4638  }
4639 
4640  if (CT == 7) {
4641  ff = f2;
4642  mmax = 4.0 / pow(ss, 2) *
4643  (4.0 / 9.0 * (pow(ss, 2) + pow((ss - tmin), 2)) / pow(tmin, 2) +
4644  (pow(tmin, 2) + pow((ss - tmin), 2)) / pow(ss, 2) +
4645  2.0 / 3.0 * pow((ss - tmin), 2) / ss / tmin);
4646  msq = (pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
4647  (4.0 / 9.0 *
4648  (((pow(ss, 2) + pow(uu, 2)) / pow(tt, 2)) +
4649  (pow(tt, 2) + pow(uu, 2)) / pow(ss, 2) +
4650  2.0 / 3.0 * pow(uu, 2) / ss / tt))) /
4651  mmax;
4652  }
4653 
4654  if (CT == 8) {
4655  ff = f2;
4656  mmax = 4.0 / pow(ss, 2) *
4657  (4.0 / 9.0 * (pow(tmin, 2) + pow((ss - tmin), 2)) / tmin /
4658  (ss - tmin) -
4659  (pow(tmin, 2) + pow((ss - tmin), 2)) / pow(ss, 2));
4660  msq = pow((1.0 / p0[0] / p2[0] / 2.0), 2) *
4661  (4.0 / 9.0 * (pow(tt, 2) + pow(uu, 2)) / tt / uu -
4662  (pow(tt, 2) + pow(uu, 2)) / pow(ss, 2)) /
4663  (mmax + 4.0);
4664  }
4665 
4666  rank = ran0(&NUM1);
4667  } while (rank > (msq * ff));
4668 
4669  if(flag1 == 1 || flag2 == 1){ // scatterings cannot be properly sampled
4670  transback(v0, p0);
4671  transback(v0, p4);
4672  qt = 0;
4673  p2[0] = 0;
4674  p2[1] = 0;
4675  p2[2] = 0;
4676  p2[3] = 0;
4677  p3[0] = 0;
4678  p3[1] = 0;
4679  p3[2] = 0;
4680  p3[3] = 0;
4681  return;
4682  }
4683 
4684  //
4685  p3[1] = p2[1];
4686  p3[2] = p2[2];
4687  p3[3] = p2[3];
4688  p3[0] = p2[0];
4689 
4690  // velocity of the center-of-mass
4691  //
4692  vc[1] = (p0[1] + p2[1]) / (p0[0] + p2[0]);
4693  vc[2] = (p0[2] + p2[2]) / (p0[0] + p2[0]);
4694  vc[3] = (p0[3] + p2[3]) / (p0[0] + p2[0]);
4695  //
4696  // transform into the cms frame
4697  //
4698  trans(vc, p0);
4699  trans(vc, p2);
4700  //
4701  // cm momentum
4702  //
4703  double pcm = p2[0];
4704  //
4705  // sample transverse momentum transfer with respect to jet momentum
4706  // in cm frame
4707  //
4708  double ranp = 2.0 * pi * ran0(&NUM1);
4709  //
4710  // transverse momentum transfer
4711  //
4712  qt = sqrt(pow(pcm, 2) - pow((tt / 2.0 / pcm - pcm), 2));
4713  double qx = qt * cos(ranp);
4714  double qy = qt * sin(ranp);
4715 
4716  //
4717  // longitudinal momentum transfer
4718  //
4719  double qpar = tt / 2.0 / pcm;
4720  //
4721  // qt is perpendicular to pcm, need to rotate back to the cm frame
4722  //
4723  double upt = sqrt(p2[1] * p2[1] + p2[2] * p2[2]) / p2[0];
4724  double upx = p2[1] / p2[0];
4725  double upy = p2[2] / p2[0];
4726  double upz = p2[3] / p2[0];
4727  //
4728  // momentum after collision in cm frame
4729  //
4730  p2[1] = p2[1] - qpar * upx;
4731  p2[2] = p2[2] - qpar * upy;
4732  if (upt != 0.0) {
4733  p2[1] = p2[1] + (upz * upx * qy + upy * qx) / upt;
4734  p2[2] = p2[2] + (upz * upy * qy - upx * qx) / upt;
4735  }
4736  p2[3] = p2[3] - qpar * upz - upt * qy;
4737 
4738  p0[1] = -p2[1];
4739  p0[2] = -p2[2];
4740  p0[3] = -p2[3];
4741  //
4742  // transform from cm back to the comoving frame
4743  //
4744  transback(vc, p2);
4745  transback(vc, p0);
4746 
4747  //************************************************************
4748  //
4749  // calculate qt in the rest frame of medium
4750  //
4751  // if(p0[4]>p2[4])
4752  // {
4753  rotate(p4[1], p4[2], p4[3], p0, 1);
4754  qt = sqrt(pow(p0[1], 2) + pow(p0[2], 2));
4755  rotate(p4[1], p4[2], p4[3], p0, -1);
4756  // }
4757  // else
4758  // {
4759  // rotate(p4[1],p4[2],p4[3],p2,1);
4760  // qt=sqrt(pow(p2[1],2)+pow(p2[2],2));
4761  // rotate(p4[1],p4[2],p4[3],p2,-1);
4762  // }
4763  //************************************************************
4764 
4765  //
4766  // transform from comoving frame to the lab frame
4767  //
4768  transback(v0, p2);
4769  transback(v0, p0);
4770  transback(v0, p3);
4771 
4772  //************************************************************
4773  transback(v0, p4);
4774  //************************************************************
4775 }
4776 
4777 //.........................................................................
4778 void Matter::collHQ22(int CT, double temp, double qhat0ud, double v0[4],
4779  double p0[4], double p2[4], double p3[4], double p4[4],
4780  double &qt) {
4781  //
4782  // HQ 2->2 scatterings
4783  // p0 initial HQ momentum, output to final momentum
4784  // p2 final thermal momentum, p3 initial thermal energy
4785  //
4786  // amss=sqrt(abs(p0(4)**2-p0(1)**2-p0(2)**2-p0(3)**2))
4787  //
4788  //************************************************************
4789 
4790  // transform to local comoving frame of the fluid
4791  trans(v0, p0);
4792 
4793  //************************************************************
4794 
4795  // sample the medium parton thermal momentum in the comoving frame
4796 
4797  double xw;
4798  double razim;
4799  double rcos;
4800  double rsin;
4801 
4802  double ss;
4803 
4804  double rant;
4805  double tt;
4806 
4807  double uu;
4808  double ff = 0.0;
4809  double rank;
4810 
4811  double msq = 0.0;
4812 
4813  double e2, theta2, theta4, phi24; // the four independent variables
4814  double e1, e4, p1, cosTheta24, downFactor,
4815  sigFactor; // other useful variables
4816  double HQmass, fBmax, fFmax, fB, fF, maxValue;
4817  int index_p1, index_T, index_e2;
4818  int ct1_loop, ct2_loop, flag1, flag2;
4819 
4820  flag1 = 0;
4821  flag2 = 0;
4822 
4823  // continue this function for HQ scattering
4824 
4825  HQmass = p0[0] * p0[0] - p0[1] * p0[1] - p0[2] * p0[2] - p0[3] * p0[3];
4826  if (HQmass > 1e-12) {
4827  HQmass = sqrt(HQmass);
4828  } else {
4829  HQmass = 0.0;
4830  }
4831 
4832  // Initial 4-momentum of HQ
4833  //
4834  //************************************************************
4835  p4[1] = p0[1];
4836  p4[2] = p0[2];
4837  p4[3] = p0[3];
4838  p4[0] = p0[0];
4839  //************************************************************
4840 
4841  p1 = sqrt(p0[1] * p0[1] + p0[2] * p0[2] + p0[3] * p0[3]);
4842  index_p1 = (int)((p1 - min_p1) / bin_p1);
4843  index_T = (int)((temp - min_T) / bin_T);
4844  if (index_p1 >= N_p1) {
4845  index_p1 = N_p1 - 1;
4846  cout << "warning: p1 is over p_max: " << p1 << endl;
4847  }
4848  if (index_T >= N_T) {
4849  index_T = N_T - 1;
4850  cout << "warning: T is over T_max: " << temp << endl;
4851  }
4852  if (index_T < 0) {
4853  index_T = 0;
4854  cout << "warning: T is below T_min: " << temp << endl;
4855  }
4856 
4857  fBmax = distFncBM[index_T][index_p1];
4858  fFmax = distFncFM[index_T][index_p1]; // maximum of f(xw) at given p1 and T
4859 
4860  maxValue = 10.0; // need actual value later
4861 
4862  ct1_loop = 0;
4863  do { // sample p2 (light parton) using distribution integrated over 3 angles
4864  ct1_loop++;
4865  if (ct1_loop > 1e6) {
4866  // cout << "cannot sample light parton for HQ scattering ..." << endl;
4867  flag1 = 1;
4868  break;
4869  }
4870  xw = max_e2 * ran0(&NUM1);
4871  index_e2 = (int)((xw - min_e2) / bin_e2);
4872  if (index_e2 >= N_e2)
4873  index_e2 = N_e2 - 1;
4874  if (CT == 11) { // qc->qc
4875  ff = distFncF[index_T][index_p1][index_e2] / fFmax;
4876  maxValue = distMaxF[index_T][index_p1][index_e2];
4877  } else if (CT == 12) { // gc->gc
4878  ff = distFncB[index_T][index_p1][index_e2] / fBmax;
4879  maxValue = distMaxB[index_T][index_p1][index_e2];
4880  } else {
4881  cout << "Wrong HQ channel ID" << endl;
4882  exit(EXIT_FAILURE);
4883  }
4884  } while (ran0(&NUM1) > ff);
4885 
4886  e2 = xw * temp;
4887  e1 = p0[0];
4888 
4889  // now e2 is fixed, need to sample the remaining 3 variables
4890  ct2_loop = 0;
4891  do {
4892  ct2_loop++;
4893  if (ct2_loop > 1e6) {
4894  cout << "cannot sample final states for HQ scattering ..." << endl;
4895  flag2 = 1;
4896  break;
4897  }
4898 
4899  theta2 = pi * ran0(&NUM1);
4900  theta4 = pi * ran0(&NUM1);
4901  phi24 = 2.0 * pi * ran0(&NUM1);
4902 
4903  cosTheta24 =
4904  sin(theta2) * sin(theta4) * cos(phi24) + cos(theta2) * cos(theta4);
4905  downFactor = e1 - p1 * cos(theta4) + e2 - e2 * cosTheta24;
4906  e4 = (e1 * e2 - p1 * e2 * cos(theta2)) / downFactor;
4907  sigFactor = sin(theta2) * sin(theta4) * e2 * e4 / downFactor;
4908 
4909  // calculate s,t,u, different definition from light quark -- tt, uu are negative
4910  ss = 2.0 * e1 * e2 + HQmass * HQmass - 2.0 * p1 * e2 * cos(theta2);
4911  tt = -2.0 * e2 * e4 * (1.0 - cosTheta24);
4912  uu = 2.0 * HQmass * HQmass - ss - tt;
4913 
4914  // re-sample if the kinematic cuts are not satisfied
4915  if (ss <= 2.0 * qhat0ud || tt >= -qhat0ud || uu >= -qhat0ud) {
4916  rank = ran0(&NUM1);
4917  sigFactor = 0.0;
4918  msq = 0.0;
4919  continue;
4920  }
4921 
4922  if (CT == 11) { // qc->qc
4923  ff =
4924  (1.0 / (exp(e2 / temp) + 1.0)) * (1.0 - 1.0 / (exp(e4 / temp) + 1.0));
4925  sigFactor = sigFactor * ff;
4926  msq = Mqc2qc(ss, tt, HQmass) / maxValue;
4927  }
4928 
4929  if (CT == 12) { // gc->gc
4930  ff =
4931  (1.0 / (exp(e2 / temp) - 1.0)) * (1.0 + 1.0 / (exp(e4 / temp) - 1.0));
4932  sigFactor = sigFactor * ff;
4933  msq = Mgc2gc(ss, tt, HQmass) / maxValue;
4934  }
4935 
4936  rank = ran0(&NUM1);
4937 
4938  } while (rank > (msq * sigFactor));
4939 
4940  if (flag1 == 0 && flag2 == 0) {
4941 
4942  // pass p2 value to p3 for initial thermal parton
4943  p3[1] = e2 * sin(theta2);
4944  p3[2] = 0.0;
4945  p3[3] = e2 * cos(theta2);
4946  p3[0] = e2;
4947 
4948  // calculate momenta of outgoing particles
4949  // here p2 is for p4 (light parton) in my note
4950 
4951  p2[1] = e4 * sin(theta4) * cos(phi24);
4952  p2[2] = e4 * sin(theta4) * sin(phi24);
4953  p2[3] = e4 * cos(theta4);
4954  p2[0] = e4;
4955 
4956  // rotate randomly in xy plane (jet is in z), because p3 is assigned in xz plane with bias
4957  double th_rotate = 2.0 * pi * ran0(&NUM1);
4958  double p3x_rotate = p3[1] * cos(th_rotate) - p3[2] * sin(th_rotate);
4959  double p3y_rotate = p3[1] * sin(th_rotate) + p3[2] * cos(th_rotate);
4960  double p2x_rotate = p2[1] * cos(th_rotate) - p2[2] * sin(th_rotate);
4961  double p2y_rotate = p2[1] * sin(th_rotate) + p2[2] * cos(th_rotate);
4962  p3[1] = p3x_rotate;
4963  p3[2] = p3y_rotate;
4964  p2[1] = p2x_rotate;
4965  p2[2] = p2y_rotate;
4966 
4967  // Because we treated p0 (p1 in my note for heavy quark) as the z-direction, proper rotations are necessary here
4968  rotate(p4[1], p4[2], p4[3], p2, -1);
4969  rotate(p4[1], p4[2], p4[3], p3, -1);
4970 
4971  p0[1] = p4[1] + p3[1] - p2[1];
4972  p0[2] = p4[2] + p3[2] - p2[2];
4973  p0[3] = p4[3] + p3[3] - p2[3];
4974  p0[0] =
4975  sqrt(p0[1] * p0[1] + p0[2] * p0[2] + p0[3] * p0[3] + HQmass * HQmass);
4976 
4977  // Debug
4978  if (fabs(p0[0] + p2[0] - p3[0] - p4[0]) > 0.00001) {
4979  cout << "Violation of energy conservation in HQ 2->2 scattering: "
4980  << fabs(p0[0] + p2[0] - p3[0] - p4[0]) << endl;
4981  }
4982 
4983  // calculate qt in the rest frame of medium
4984  rotate(p4[1], p4[2], p4[3], p0, 1);
4985  qt = sqrt(pow(p0[1], 2) + pow(p0[2], 2));
4986  rotate(p4[1], p4[2], p4[3], p0, -1);
4987 
4988  // transform from comoving frame to the lab frame
4989  transback(v0, p2);
4990  transback(v0, p0);
4991  transback(v0, p3);
4992  transback(v0, p4);
4993 
4994  } else { // no scattering
4995  transback(v0, p0);
4996  transback(v0, p4);
4997  qt = 0;
4998  p2[0] = 0;
4999  p2[1] = 0;
5000  p2[2] = 0;
5001  p2[3] = 0;
5002  p3[0] = 0;
5003  p3[1] = 0;
5004  p3[2] = 0;
5005  p3[3] = 0;
5006  }
5007 }
5008 
5009 double Matter::Mqc2qc(double s, double t, double M) {
5010 
5011  double m2m = M * M;
5012  double u = 2.0 * m2m - s - t;
5013  double MM;
5014 
5015  MM = 64.0 / 9.0 * (pow((m2m - u), 2) + pow((s - m2m), 2) + 2.0 * m2m * t) /
5016  t / t;
5017 
5018  return (MM);
5019 }
5020 
5021 double Matter::Mgc2gc(double s, double t, double M) {
5022 
5023  double m2m = M * M;
5024  double u = 2.0 * m2m - s - t;
5025  double MM;
5026 
5027  MM = 32.0 * (s - m2m) * (m2m - u) / t / t;
5028  MM = MM + 64.0 / 9.0 * ((s - m2m) * (m2m - u) + 2.0 * m2m * (s + m2m)) /
5029  pow((s - m2m), 2);
5030  MM = MM + 64.0 / 9.0 * ((s - m2m) * (m2m - u) + 2.0 * m2m * (u + m2m)) /
5031  pow((u - m2m), 2);
5032  MM = MM + 16.0 / 9.0 * m2m * (4.0 * m2m - t) / ((s - m2m) * (m2m - u));
5033  MM = MM + 16.0 * ((s - m2m) * (m2m - u) + m2m * (s - u)) / (t * (s - m2m));
5034  MM = MM + 16.0 * ((s - m2m) * (m2m - u) - m2m * (s - u)) / (t * (u - m2m));
5035 
5036  return (MM);
5037 }
5038 
5039 void Matter::trans(double v[4], double p[4]) {
5040  double vv = sqrt(v[1] * v[1] + v[2] * v[2] + v[3] * v[3]);
5041  double ga = 1.0 / sqrt(1.0 - vv * vv);
5042  double ppar = p[1] * v[1] + p[2] * v[2] + p[3] * v[3];
5043  double gavv = (ppar * ga / (1.0 + ga) - p[0]) * ga;
5044  p[0] = ga * (p[0] - ppar);
5045  p[1] = p[1] + v[1] * gavv;
5046  p[2] = p[2] + v[2] * gavv;
5047  p[3] = p[3] + v[3] * gavv;
5048 }
5049 
5050 void Matter::transback(double v[4], double p[4]) {
5051  double vv = sqrt(v[1] * v[1] + v[2] * v[2] + v[3] * v[3]);
5052  double ga = 1.0 / sqrt(1.0 - vv * vv);
5053  double ppar = p[1] * v[1] + p[2] * v[2] + p[3] * v[3];
5054  double gavv = (-ppar * ga / (1.0 + ga) - p[0]) * ga;
5055  p[0] = ga * (p[0] + ppar);
5056  p[1] = p[1] - v[1] * gavv;
5057  p[2] = p[2] - v[2] * gavv;
5058  p[3] = p[3] - v[3] * gavv;
5059 }
5060 
5061 void Matter::rotate(double px, double py, double pz, double pr[4], int icc) {
5062  // input: (px,py,pz), (wx,wy,wz), argument (i)
5063  // output: new (wx,wy,wz)
5064  // if i=1, turn (wx,wy,wz) in the direction (px,py,pz)=>(0,0,E)
5065  // if i=-1, turn (wx,wy,wz) in the direction (0,0,E)=>(px,py,pz)
5066 
5067  double wx, wy, wz, E, pt, w, cosa, sina, cosb, sinb;
5068  double wx1, wy1, wz1;
5069 
5070  wx = pr[1];
5071  wy = pr[2];
5072  wz = pr[3];
5073 
5074  E = sqrt(px * px + py * py + pz * pz);
5075  pt = sqrt(px * px + py * py);
5076 
5077  w = sqrt(wx * wx + wy * wy + wz * wz);
5078 
5079  if (pt == 0) {
5080  cosa = 1;
5081  sina = 0;
5082  } else {
5083  cosa = px / pt;
5084  sina = py / pt;
5085  }
5086 
5087  cosb = pz / E;
5088  sinb = pt / E;
5089 
5090  if (icc == 1) {
5091  wx1 = wx * cosb * cosa + wy * cosb * sina - wz * sinb;
5092  wy1 = -wx * sina + wy * cosa;
5093  wz1 = wx * sinb * cosa + wy * sinb * sina + wz * cosb;
5094  }
5095 
5096  else {
5097  wx1 = wx * cosa * cosb - wy * sina + wz * cosa * sinb;
5098  wy1 = wx * sina * cosb + wy * cosa + wz * sina * sinb;
5099  wz1 = -wx * sinb + wz * cosb;
5100  }
5101 
5102  wx = wx1;
5103  wy = wy1;
5104  wz = wz1;
5105 
5106  pr[1] = wx;
5107  pr[2] = wy;
5108  pr[3] = wz;
5109 
5110  // pr[0]=sqrt(pr[1]*pr[1]+pr[2]*pr[2]+pr[3]*pr[3]);
5111 }
5112 
5113 float Matter::ran0(long *idum)
5114 
5115 {
5116  const int IM1 = 2147483563;
5117  const int IM2 = 2147483399;
5118  const double AM = (1.0 / IM1);
5119  const int IMM1 = (IM1 - 1);
5120  const int IA1 = 40014;
5121  const int IA2 = 40692;
5122  const int IQ1 = 53668;
5123  const int IQ2 = 52774;
5124  const int IR1 = 12211;
5125  const int IR2 = 3791;
5126  const int NTAB = 32;
5127  const int NDIV = (1 + IMM1 / NTAB);
5128  const double EPS = 1.2e-7;
5129  const double RNMX = (1.0 - EPS);
5130 
5131  int j;
5132  long k;
5133  static long idum2 = 123456789;
5134  static long iy = 0;
5135  static long iv[NTAB];
5136  float temp;
5137 
5138  if (*idum <= 0) {
5139  if (-(*idum) < 1)
5140  *idum = 1;
5141  else
5142  *idum = -(*idum);
5143  for (j = NTAB + 7; j >= 0; j--) {
5144  k = (*idum) / IQ1;
5145  *idum = IA1 * (*idum - k * IQ1) - k * IR1;
5146  if (*idum < 0)
5147  *idum += IM1;
5148  if (j < NTAB)
5149  iv[j] = *idum;
5150  }
5151  iy = iv[0];
5152  }
5153  k = (*idum) / IQ1;
5154  *idum = IA1 * (*idum - k * IQ1) - k * IR1;
5155  if (*idum < 0)
5156  *idum += IM1;
5157  k = idum2 / IQ2;
5158  idum2 = IA2 * (idum2 - k * IQ2) - k * IR2;
5159  if (idum2 < 0)
5160  idum2 += IM2;
5161  j = iy / NDIV;
5162  iy = iv[j] - idum2;
5163  iv[j] = *idum;
5164  if (iy < 1)
5165  iy += IMM1;
5166  if ((temp = AM * iy) > RNMX)
5167  return RNMX;
5168  else
5169  return temp;
5170 }
5171 
5173 
5174 double Matter::solve_alphas(double var_qhat, double var_ener, double var_temp) {
5175 
5176  double preFactor = 42.0 * Ca * zeta3 / pi;
5177 
5178  // reference: qhatLoc = Ca*50.4864/pi*pow(alphas,2)*pow(tempLoc,3)*log(5.7*2.0*pi*tempLoc*tempLoc/4.0/muD2);
5179  double max_qhat =
5180  preFactor * pow(0.5, 2) * pow(var_temp, 3) *
5181  log(5.7 * max(var_ener, 2.0 * pi * var_temp) / 24 / pi / 0.5 / var_temp);
5182 
5183  if (max_qhat < var_qhat) {
5184  JSINFO << "qhat exceeds HTL calculation, use alpha_s = 0.5";
5185  return (0.5);
5186  }
5187 
5188  double solution = sqrt(
5189  var_qhat / preFactor / pow(var_temp, 3) /
5190  log(5.7 * max(var_ener, 2.0 * pi * var_temp) / 24 / pi / 0.2 / var_temp));
5191  double fnc_value, fnc_derivative;
5192  fnc_value = fnc0_alphas(solution, var_qhat, var_ener, var_temp);
5193  fnc_derivative =
5194  fnc0_derivative_alphas(solution, var_qhat, var_ener, var_temp);
5195 
5196  //cout << "initial guess: " << solution << " " << fnc_value << endl;
5197 
5198  while (fabs(fnc_value / var_qhat) > 0.001) {
5199 
5200  solution = solution - fnc_value / fnc_derivative;
5201  fnc_value = fnc0_alphas(solution, var_qhat, var_ener, var_temp);
5202  fnc_derivative =
5203  fnc0_derivative_alphas(solution, var_qhat, var_ener, var_temp);
5204  }
5205 
5206  if (solution < 0.0 || solution > 0.5) {
5207  JSINFO << "unreasonable alpha_s: " << solution << " use alpha_s = 0.5";
5208  solution = 0.5;
5209  }
5210 
5211  return (solution);
5212 }
5213 
5214 double Matter::fnc0_alphas(double var_alphas, double var_qhat, double var_ener,
5215  double var_temp) {
5216 
5217  double preFactor = 42.0 * Ca * zeta3 / pi;
5218  return (preFactor * var_alphas * var_alphas * pow(var_temp, 3) *
5219  log(5.7 * max(var_ener, 2.0 * pi * var_temp) / 24 / pi /
5220  var_alphas / var_temp) -
5221  var_qhat);
5222 }
5223 
5224 double Matter::fnc0_derivative_alphas(double var_alphas, double var_qhat,
5225  double var_ener, double var_temp) {
5226 
5227  double preFactor = 42.0 * Ca * zeta3 / pi;
5228  return (preFactor * pow(var_temp, 3) *
5229  (2.0 * var_alphas *
5230  log(5.7 * max(var_ener, 2.0 * pi * var_temp) / 24 / pi /
5231  var_alphas / var_temp) -
5232  var_alphas));
5233 }
5234 
5235 void Matter::read_tables() { // intialize various tables for LBT
5236 
5237  //...read scattering rate
5238  int it, ie;
5239  int n = 450;
5240  //ifstream f1("LBT-tables/ratedata");
5241  //if(!f1.is_open())
5242  // {
5243  // cout<<"Erro openning date file1!\n";
5244  // }
5245  //else
5246  // {
5247  // for(int i=1;i<=n;i++)
5248  // {
5249  // f1>>it>>ie;
5250  // f1>>qhatG[it][ie]>>Rg[it][ie]>>Rg1[it][ie]>>Rg2[it][ie]>>Rg3[it][ie]>>qhatLQ[it][ie]>>Rq[it][ie]>>Rq3[it][ie]>>Rq4[it][ie]>>Rq5[it][ie]>>Rq6[it][ie]>>Rq7[it][ie]>>Rq8[it][ie];
5251  // }
5252  // }
5253  //f1.close();
5254 
5255  // duplicate for heavy quark
5256  ifstream f11("LBT-tables/ratedata-HQ");
5257  if (!f11.is_open()) {
5258  cout << "Erro openning HQ data file!\n";
5259  } else {
5260  for (int i = 1; i <= n; i++) {
5261  f11 >> it >> ie;
5262  f11 >> RHQ[it][ie] >> RHQ11[it][ie] >> RHQ12[it][ie] >> qhatHQ[it][ie];
5263  }
5264  }
5265  f11.close();
5266 
5267  // preparation for HQ 2->2
5268  ifstream fileB("LBT-tables/distB.dat");
5269  if (!fileB.is_open()) {
5270  cout << "Erro openning data file distB.dat!" << endl;
5271  } else {
5272  for (int i = 0; i < N_T; i++) {
5273  for (int j = 0; j < N_p1; j++) {
5274  double dummy_T, dummy_p1;
5275  fileB >> dummy_T >> dummy_p1;
5276  if (fabs(min_T + (0.5 + i) * bin_T - dummy_T) > 1.0e-5 ||
5277  fabs(min_p1 + (0.5 + j) * bin_p1 - dummy_p1) > 1.0e-5) {
5278  cout << "Erro in reading data file distB.dat!" << endl;
5279  exit(EXIT_FAILURE);
5280  }
5281  fileB >> distFncBM[i][j];
5282  for (int k = 0; k < N_e2; k++)
5283  fileB >> distFncB[i][j][k];
5284  for (int k = 0; k < N_e2; k++)
5285  fileB >> distMaxB[i][j][k];
5286  }
5287  }
5288  }
5289  fileB.close();
5290 
5291  ifstream fileF("LBT-tables/distF.dat");
5292  if (!fileF.is_open()) {
5293  cout << "Erro openning data file distF.dat!" << endl;
5294  } else {
5295  for (int i = 0; i < N_T; i++) {
5296  for (int j = 0; j < N_p1; j++) {
5297  double dummy_T, dummy_p1;
5298  fileF >> dummy_T >> dummy_p1;
5299  if (fabs(min_T + (0.5 + i) * bin_T - dummy_T) > 1.0e-5 ||
5300  fabs(min_p1 + (0.5 + j) * bin_p1 - dummy_p1) > 1.0e-5) {
5301  cout << "Erro in reading data file distF.dat!" << endl;
5302  exit(EXIT_FAILURE);
5303  }
5304  fileF >> distFncFM[i][j];
5305  for (int k = 0; k < N_e2; k++)
5306  fileF >> distFncF[i][j][k];
5307  for (int k = 0; k < N_e2; k++)
5308  fileF >> distMaxF[i][j][k];
5309  }
5310  }
5311  }
5312  fileF.close();
5313 }