Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LBT.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file LBT.h
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 #ifndef LBT_H
17 #define LBT_H
18 
19 #include "JetEnergyLossModule.h"
20 #include <iostream>
21 #include <string>
22 #include <sstream>
23 #include <fstream>
24 
25 using namespace Jetscape;
26 
27 //class LBTUserInfo: public Parton::PseudoJet::UserInfoBase {
28 class LBTUserInfo : public fjcore::PseudoJet::UserInfoBase {
29 public:
30  LBTUserInfo(double ttt) : _lrf_T_tot(ttt){};
31  double lrf_T_tot() const { return _lrf_T_tot; }
32  double _lrf_T_tot;
33 };
34 
35 //variables for unit test
36 //double de1[41]={0.0};
37 //double de2[41]={0.0};
38 //int ctEvt=0;
39 
40 class LBT : public JetEnergyLossModule<
41  LBT> //, public std::enable_shared_from_this<Matter>
42 {
43 public:
44  LBT();
45  virtual ~LBT();
46 
47  void Init();
48  //void Exec();
49  //void DoEnergyLoss(double deltaT, double Q2, const vector<Parton>& pIn, vector<Parton>& pOut);
50  void DoEnergyLoss(double deltaT, double time, double Q2, vector<Parton> &pIn,
51  vector<Parton> &pOut);
52  void WriteTask(weak_ptr<JetScapeWriter> w);
53 
54 private:
56  //
57  // Define variables and functions for LBT
58  //
60 
61  //...constant
62  const double pi = 3.1415926;
63  const double epsilon = 1e-6;
64  const double CA = 3.0;
65  const double CF = 4.0 / 3.0;
66  const double sctr = 0.1973; //fm to GeV^-1
67 
68  //...fixed parameters (cannot change using input parameter file)
69  const int lightOut = 1; // 1 -> write out light parton information
70  const int heavyOut = 0; // 1 -> write out heavy quark information
71  const int outFormat =
72  2; // different output formats: 1 - my preferred, 2 - for JETSCAPE patch code
73  const double cutOut =
74  epsilon; // neglect light partons in output files for E < cutOut;
75 
76  //...paramters for the bulk, can be changed from input parameter file
77  int vacORmed = 1; // 0-vacuum; 1-medium
78  int bulkFlag = 0; // 0-static medium; 1-OSU hydro; 2-CCNU hydro
79  double temp0 = 0.25; // medium temperature
80  double hydro_Tc = 0.165;
81  //...derived quantities
82  double temp00;
83 
84  //...input parameters for LBT
85  double KPamp = 0.0;
86  double KPsig = 5.0;
87  double KTamp = 0.0;
88  double KTsig = 0.05 * hydro_Tc;
89  double preKT = 0.3 / 0.3;
90  double scaleAK = 2.0;
91  //...derived quantities
92  double KPfactor, KTfactor, Kfactor, runKT;
93  double fixedLog, runLog, scaleMu2, runAlphas;
94 
95  //...initialization parameters for jet partons
96  int initHardFlag =
97  2; //1 initialize by the code itself; 2 initialize by reading in particle list
98  int fixMomentum = 0;
99  int fixPosition = 1;
100  int run_alphas = 1;
101  int flagJetX =
102  0; // 0: do nothing; 1: keep momentum but reset jet position within LBT
103  int Kjet = 21; //initial flavor of the jet parton
104  double ipTmin = 0.0;
105  double ipTmax = 800.0;
106  double eta_cut = 0.5;
107  double ener = 50.0; //initial energy of the jet parton/heavy quark
108  double amss = 0.0; //initial mass of the jet parton/heavy quark
109  int nj = 1000; //number of initial partons per event
110  int ncall = 1; // number of events
111  //...derived quantities
112  int np; //number of partons at this time step ti
113 
114  //...clock information
115  int nprint = 100; // print out infomation per # of events
116  int tauswitch =
117  0; // 0 for t-z coordinate, 1 (not included in this version) for tau-eta
118  double tau0 = 0.0; // initial time
119  double dtau = 0.1; // time interval
120  double tauend = 4.0; // time when program ends
121  //...derived quantities
122  double time0;
123  double
124  dt; //dtime when tauswitch is turned off (0) / dtau when tauswitch is turned on (1)
125  double timend; //end time or tau RENAME
126 
127  //...switches and cuts
128  int Kprimary =
129  0; // 0 keep all partons, 1 keep leading jet parton only (switch off other partons)
130  int KINT0 = 1; // 0 no radiation, 1 elastic + inelastic
131  int Kqhat0 = 2; //Debye screening mass switch RENAME
132  int Kalphas = 1; //alphas switch
133  double Ecut = 0.0; //energy cut of the recoiled partons
134  double fixAlphas = 0.3;
135  //...derived quantities
136  int KINT; //radiation switch
137  double alphas;
138  double qhat0; //Debye mass RENAME
139  double qhat00;
140 
141  std::uniform_real_distribution<double> ZeroOneDistribution;
142 
143  // flag to make sure initialize only once
144  static bool flag_init;
145 
146  // scattering rate
147  static double Rg
148  [60]
149  [20]; //total gluon scattering rate as functions of initial energy and temperature
150  static double Rg1[60][20]; //gg-gg CT1
151  static double Rg2[60][20]; //gg-qqbar CT2
152  static double Rg3[60][20]; //gq-qg CT3
153  static double Rq
154  [60]
155  [20]; //total gluon scattering rate as functions of initial energy and temperature
156  static double Rq3[60][20]; //qg-qg CT13
157  static double Rq4[60][20]; //qiqj-qiqj CT4
158  static double Rq5[60][20]; //qiqi-qiqi CT5
159  static double Rq6[60][20]; //qiqibar-qjqjbar CT6
160  static double Rq7[60][20]; //qiqibar-qiqibar CT7
161  static double Rq8[60][20]; //qqbar-gg CT8
162  static double qhatLQ[60][20];
163  static double qhatG[60][20];
164 
165  static double RHQ[60][20]; //total scattering rate for heavy quark
166  static double RHQ11[60][20]; //Qq->Qq
167  static double RHQ12[60][20]; //Qg->Qg
168  static double qhatHQ[60][20]; //qhat of heavy quark
169  double qhat_over_T3; //qhat/T^3 for heavy quark as fnc of (T,p)
170  double D2piT;
171 
172  // for heavy quark radiation table
173  static const int HQener_gn = 500;
174  static const int t_gn_1 = 100;
175  static const int t_gn_2 = 125;
176  static const int t_gn = t_gn_1 + t_gn_2;
177  static const int temp_gn = 100;
178 
179  static double dNg_over_dt_c[t_gn + 2][temp_gn + 1][HQener_gn + 1];
180  static double dNg_over_dt_q[t_gn + 2][temp_gn + 1][HQener_gn + 1];
181  static double dNg_over_dt_g[t_gn + 2][temp_gn + 1][HQener_gn + 1];
182  static double max_dNgfnc_c[t_gn + 2][temp_gn + 1][HQener_gn + 1];
183  static double max_dNgfnc_q[t_gn + 2][temp_gn + 1][HQener_gn + 1];
184  static double max_dNgfnc_g[t_gn + 2][temp_gn + 1][HQener_gn + 1];
185 
186  const double HQener_max = 1000.0;
187  const double t_max_1 = 20.0;
188  const double t_max_2 = 270.0;
189  const double t_max = t_max_2;
190  const double temp_max = 0.65;
191  const double temp_min = 0.15;
192  double delta_tg_1 = t_max_1 / t_gn_1;
193  double delta_tg_2 = (t_max_2 - t_max_1) / t_gn_2;
194  double delta_temp = (temp_max - temp_min) / temp_gn;
195  double delta_HQener = HQener_max / HQener_gn;
196 
197  // for MC initialization of jet partons
198  static const int maxMC = 2000000;
199  static double initMCX[maxMC], initMCY[maxMC];
200 
201  //...radiation block
202  int icl22;
203  int icl23; //the numerical switch in colljet23
204  int iclrad; //the numerical switch in radiation
205  int isp; //the splitting function switch
206 
207  //...global variable qhat
208  int counth100 = 0;
209 
210  double qhat; //transport parameter
211 
212  double dng0[101][101] = {{0.0}}; //table of dn/dkperp2/dx
213 
214  double Vtemp[4] = {0.0};
215 
216  //...time system
217 
218  static const int dimParList = 50; // originally 500000
219 
220  double tirad[dimParList] = {0.0};
221  double tiscatter[dimParList] = {0.0};
222  double tiform[dimParList] = {0.0}; //pythia undone
223  double Tint_lrf[dimParList] = {0.0}; //for heavy quark
224  double eGluon = 0.0;
225  double nGluon = 0.0;
226  double dEel = 0.0;
227 
228  //....radiated gluon
229  double radng[dimParList] = {0.0};
230  double xwm[3] = {0.0};
231  double wkt2m;
232 
233  double vf[4] = {0.0}; //flow velocity in tau-eta coordinate
234  double vfcar[4] = {0.0}; //flow velocity in t-z coordinate
235 
236  double vp[4] = {0.0}; //position of particle
237  double vc0[4] = {0.0}; //flow velocity
238 
239  //...dimensions in subroutine colljet and twcoll
240  double vc[4] = {0.0};
241  double pc0[4] = {0.0};
242  double pc2[4] = {0.0};
243  double pc3[4] = {0.0};
244  double pc4[4] = {0.0};
245  double p0[4] = {0.0};
246  double p2[4] = {0.0};
247  double p3[4] = {0.0};
248  double p4[4] = {0.0};
249 
250  double pc00[4] = {0.0};
251  double pc30[4] = {0.0};
252 
253  double pc01[4] = {0.0};
254  double pb[4] = {0.0};
255 
256  double Pj0[4] = {0.0};
257 
258  double V[4][dimParList] = {{0.0}}; //parton position
259  double P[7][dimParList] = {{0.0}}; //parton 4-momentum
260  double V0[4][dimParList] = {{0.0}}; //negative parton position
261  double P0[7][dimParList] = {{0.0}}; //negative parton 4-momentum
262 
263  double Prad[4][dimParList] = {{0.0}};
264 
265  double WT[dimParList] = {0};
266  double WT0[dimParList] = {0};
267 
268  int NR[dimParList] = {0}; //scattering rank
269  int KATT1[dimParList] = {0}; //parton flavor
270  int KATT10[dimParList] = {0}; //negative parton flavor
271 
272  double PGm[4] = {0.0};
273  double tjp[dimParList] = {0.0};
274  double Vfrozen[4][dimParList] = {{0.0}}; //parton final 4 coordinate
275  double Vfrozen0[4][dimParList] = {{0.0}}; //negative parton final 4 coordinate
276  double Ecmcut = 2.0; //energy cut for free streaming
277  double Tfrozen[dimParList] = {0.0};
278  double Tfrozen0[dimParList] = {0.0};
279  double vcfrozen[4][dimParList] = {{0.0}};
280  double vcfrozen0[4][dimParList] = {{0.0}};
281 
282  double VV[4][dimParList] = {{0.0}};
283  double VV0[4][dimParList] = {{0.0}};
284  double PP[4][dimParList] = {{0.0}};
285  double PP0[4][dimParList] = {{0.0}};
286  int CAT[dimParList] = {0};
287  int CAT0[dimParList] = {0};
288 
289  int ncut;
290  int ncut0;
291 
292  int n_coll22 = 0;
293  int n_coll23 = 0;
294  int ng_coll23 = 0;
295  int ng_nrad = 0;
296  int n_radiation = 0;
297  int ng_radiation = 0;
298  int n_gluon = 0;
299  int n_sp1 = 0;
300  int n_sp2 = 0;
301 
302  // Variables for HQ 2->2
303  static const int N_p1 = 500;
304  static const int N_T = 60;
305  static const int N_e2 = 75;
306  static double distFncB[N_T][N_p1][N_e2], distFncF[N_T][N_p1][N_e2],
307  distMaxB[N_T][N_p1][N_e2], distMaxF[N_T][N_p1][N_e2];
308  static double distFncBM[N_T][N_p1], distFncFM[N_T][N_p1];
309  double min_p1 = 0.0;
310  double max_p1 = 1000.0;
311  double bin_p1 = (max_p1 - min_p1) / N_p1;
312  double min_T = 0.1;
313  double max_T = 0.7;
314  double bin_T = (max_T - min_T) / N_T;
315  double min_e2 = 0.0;
316  double max_e2 = 15.0;
317  double bin_e2 = (max_e2 - min_e2) / N_e2;
318 
319  //int loopN=10000;
320  int loopN = 1000;
321 
322  int numInitXY = 0;
323  int flagScatter, flag_update, flag_update0;
324  double Q00, Q0;
325 
326  double tStart;// = 0.6;
327 
328  //...functions
329 
330  void LBT0(int &n, double &ti);
331 
332  void trans(double v[4], double p[4]);
333  void transback(double v[4], double p[4]);
334 
335  float ran0(long *idum);
336 
337  double alphas0(int &Kalphas, double temp0);
338  double DebyeMass2(int &Kqhat0, double alphas, double temp0);
339 
340  void lam(int KATT0, double &RTE, double E, double T, double &T1, double &T2,
341  double &E1, double &E2, int &iT1, int &iT2, int &iE1, int &iE2);
342 
343  void flavor(int &CT, int &KATT0, int &KATT2, int &KATT3, double RTE, double E,
344  double T, double &T1, double &T2, double &E1, double &E2,
345  int &iT1, int &iT2, int &iE1, int &iE2);
346  void linear(int KATT, double E, double T, double &T1, double &T2, double &E1,
347  double &E2, int &iT1, int &iT2, int &iE1, int &iE2, double &RTEg,
348  double &RTEg1, double &RTEg2, double &RTEg3, double &RTEq,
349  double &RTEq3, double &RTEq4, double &RTEq5, double &RTEq6,
350  double &RTEq7, double &RTEq8, double &RTEHQ, double &RTEHQ11,
351  double &RTEHQ12, double &qhatTP);
352  void twflavor(int &CT, int &KATT0, int &KATT2, double E, double T);
353  void twlinear(int KATT, double E, double T, double &RTEg1, double &RTEg2,
354  double &RTEq6, double &RTEq7, double &RTEq8);
355  void colljet22(int CT, double temp, double qhat0ud, double v0[4],
356  double p0[4], double p2[4], double p3[4], double p4[4],
357  double &qt);
358  void twcoll(int CT, double qhat0ud, double v0[4], double p0[4], double p2[4]);
359 
360  void titau(double ti, double vf[4], double vp[4], double p0[4], double &Vx,
361  double &Vy, double &Veta, double &Xtau);
362 
363  void rotate(double px, double py, double pz, double p[4], int icc);
364 
365  int KPoisson(double alambda);
366 
367  void radiationHQ(int parID, double qhat0ud, double v0[4], double P2[4],
368  double P3[4], double P4[4], double Pj0[4], int &ic,
369  double Tdiff, double HQenergy, double max_Ng,
370  double temp_med, double xLow, double xInt);
371  void collHQ22(int CT, double temp, double qhat, double v0[4], double p0[4],
372  double p2[4], double p3[4], double p4[4], double &qt);
373  double Mqc2qc(double s, double t, double M);
374  double Mgc2gc(double s, double t, double M);
375 
376  void collHQ23(int parID, double temp_med, double qhat0ud, double v0[4],
377  double p0[4], double p2[4], double p3[4], double p4[4],
378  double qt, int &ic, double Tdiff, double HQenergy,
379  double max_Ng, double xLow, double xInt);
380  double dNg_over_dxdydt(int parID, double x0g, double y0g, double HQenergy,
381  double HQmass, double temp_med, double Tdiff);
382  double tau_f(double x0g, double y0g, double HQenergy, double HQmass);
383  double splittingP(int parID, double z0g);
384  double lambdas(double kTFnc);
385  double nflavor(double kTFnc);
386  double alphasHQ(double kTFnc, double tempFnc);
387  double nHQgluon(int parID, double dtLRF, double &time_gluon, double &temp_med,
388  double &HQenergy, double &max_Ng);
389 
390  void read_xyMC(int &numXY);
391  void jetInitialize(int numXY);
392  void setJetX(int numXY);
393  void read_tables();
394  void jetClean();
395  void setParameter(string fileName);
396  int checkParameter(int nArg);
397 
398  // extern "C" {
399  // void read_ccnu_(char *dataFN_in, int len1);
400  // void hydroinfoccnu_(double *Ct, double *Cx, double *Cy, double *Cz, double *Ctemp, double *Cvx, double *Cvy, double *Cvz, int *Cflag);
401  //
402  // void sethydrofilesez_(int *dataID_in, char *dataFN_in, int *ctlID_in, char *ctlFN_in, int *bufferSize, int len1, int len2);
403  // void readhydroinfoshanshan_(double *t, double *x, double *y, double *z, double *e, double *s, double *temp, double *vx, double *vy, double *vz, int *flag);
404  // }
405 
406  // Allows the registration of the module so that it is available to be used by the Jetscape framework.
407  static RegisterJetScapeModule<LBT> reg;
408 };
409 
410 #endif // LBT_H