src/jet/LBT.h

0001 /*******************************************************************************
0002  * Copyright (c) The JETSCAPE Collaboration, 2018
0003  *
0004  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
0005  *
0006  * For the list of contributors see AUTHORS.
0007  *
0008  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
0009  *
0010  * or via email to bugs.jetscape@gmail.com
0011  *
0012  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
0013  * See COPYING for details.
0014  ******************************************************************************/
0015
0016 #ifndef LBT_H
0017 #define LBT_H
0018
0019 #include "JetEnergyLossModule.h"
0020 #include <iostream>
0021 #include <string>
0022 #include <sstream>
0023 #include <fstream>
0024
0025 using namespace Jetscape;
0026
0027 //class LBTUserInfo: public Parton::PseudoJet::UserInfoBase {
0028 class LBTUserInfo : public fjcore::PseudoJet::UserInfoBase {
0029 public:
0030   LBTUserInfo(double ttt) : _lrf_T_tot(ttt){};
0031   double lrf_T_tot() const { return _lrf_T_tot; }
0032   double _lrf_T_tot;
0033 };
0034
0035 //variables for unit test
0036 //double de1[41]={0.0};
0037 //double de2[41]={0.0};
0038 //int ctEvt=0;
0039
0040 class LBT : public JetEnergyLossModule<
0041                 LBT> //, public std::enable_shared_from_this<Matter>
0042 {
0043 public:
0044   LBT();
0045   virtual ~LBT();
0046
0047   void Init();
0048   //void Exec();
0049   //void DoEnergyLoss(double deltaT, double Q2, const vector<Parton>& pIn, vector<Parton>& pOut);
0050   void DoEnergyLoss(double deltaT, double time, double Q2, vector<Parton> &pIn,
0051                     vector<Parton> &pOut);
0052   void WriteTask(weak_ptr<JetScapeWriter> w);
0053
0054 private:
0055   ///////////////////////////////////////////////////////////////////////////////////////////////////
0056   //
0057   // Define variables and functions for LBT
0058   //
0059   ///////////////////////////////////////////////////////////////////////////////////////////////////
0060
0061   //...constant
0062   const double pi = 3.1415926;
0063   const double epsilon = 1e-6;
0064   const double CA = 3.0;
0065   const double CF = 4.0 / 3.0;
0066   const double sctr = 0.1973; //fm to GeV^-1
0067
0068   //...fixed parameters (cannot change using input parameter file)
0069   const int lightOut = 1; // 1 -> write out light parton information
0070   const int heavyOut = 0; // 1 -> write out heavy quark information
0071   const int outFormat =
0072       2; // different output formats: 1 - my preferred, 2 - for JETSCAPE patch code
0073   const double cutOut =
0074       epsilon; // neglect light partons in output files for E < cutOut;
0075
0076   //...paramters for the bulk, can be changed from input parameter file
0077   int vacORmed = 1;    // 0-vacuum; 1-medium
0078   int bulkFlag = 0;    // 0-static medium; 1-OSU hydro; 2-CCNU hydro
0079   double temp0 = 0.25; // medium temperature
0080   double hydro_Tc = 0.165;
0081   //...derived quantities
0082   double temp00;
0083
0084   //...input parameters for LBT
0085   double KPamp = 0.0;
0086   double KPsig = 5.0;
0087   double KTamp = 0.0;
0088   double KTsig = 0.05 * hydro_Tc;
0089   double preKT = 0.3 / 0.3;
0090   double scaleAK = 2.0;
0091   //...derived quantities
0092   double KPfactor, KTfactor, Kfactor, runKT;
0093   double fixedLog, runLog, scaleMu2, runAlphas;
0094
0095   //...initialization parameters for jet partons
0096   int initHardFlag =
0097       2; //1 initialize by the code itself; 2 initialize by reading in particle list
0098   int fixMomentum = 0;
0099   int fixPosition = 1;
0100   int run_alphas = 1;
0101   int flagJetX =
0102       0; // 0: do nothing; 1: keep momentum but reset jet position within LBT
0103   int Kjet = 21; //initial flavor of the jet parton
0104   double ipTmin = 0.0;
0105   double ipTmax = 800.0;
0106   double eta_cut = 0.5;
0107   double ener = 50.0; //initial energy of the jet parton/heavy quark
0108   double amss = 0.0;  //initial mass of the jet parton/heavy quark
0109   int nj = 1000;      //number of initial partons per event
0110   int ncall = 1;      // number of events
0111   //...derived quantities
0112   int np; //number of partons at this time step ti
0113
0114   //...clock information
0115   int nprint = 100; // print out infomation per # of events
0116   int tauswitch =
0117       0; // 0 for t-z coordinate, 1 (not included in this version) for tau-eta
0118   double tau0 = 0.0;   // initial time
0119   double dtau = 0.1;   // time interval
0120   double tauend = 4.0; // time when program ends
0121   //...derived quantities
0122   double time0;
0123   double
0124       dt; //dtime when tauswitch is turned off (0) / dtau when tauswitch is turned on (1)
0125   double timend; //end time or tau RENAME
0126
0127   //...switches and cuts
0128   int Kprimary =
0129       0; // 0 keep all partons, 1 keep leading jet parton only (switch off other partons)
0130   int KINT0 = 1;     // 0 no radiation, 1 elastic + inelastic
0131   int Kqhat0 = 2;    //Debye screening mass switch RENAME
0132   int Kalphas = 1;   //alphas switch
0133   double Ecut = 0.0; //energy cut of the recoiled partons
0134   double fixAlphas = 0.3;
0135   //...derived quantities
0136   int KINT; //radiation switch
0137   double alphas;
0138   double qhat0; //Debye mass RENAME
0139   double qhat00;
0140
0141   std::uniform_real_distribution<double> ZeroOneDistribution;
0142
0143   // flag to make sure initialize only once
0144   static bool flag_init;
0145
0146   //    scattering rate
0147   static double Rg
0148       [60]
0149       [20]; //total gluon scattering rate as functions of initial energy and temperature
0150   static double Rg1[60][20]; //gg-gg              CT1
0151   static double Rg2[60][20]; //gg-qqbar           CT2
0152   static double Rg3[60][20]; //gq-qg              CT3
0153   static double Rq
0154       [60]
0155       [20]; //total gluon scattering rate as functions of initial energy and temperature
0156   static double Rq3[60][20]; //qg-qg              CT13
0157   static double Rq4[60][20]; //qiqj-qiqj          CT4
0158   static double Rq5[60][20]; //qiqi-qiqi          CT5
0159   static double Rq6[60][20]; //qiqibar-qjqjbar    CT6
0160   static double Rq7[60][20]; //qiqibar-qiqibar    CT7
0161   static double Rq8[60][20]; //qqbar-gg           CT8
0162   static double qhatLQ[60][20];
0163   static double qhatG[60][20];
0164
0165   static double RHQ[60][20];    //total scattering rate for heavy quark
0166   static double RHQ11[60][20];  //Qq->Qq
0167   static double RHQ12[60][20];  //Qg->Qg
0168   static double qhatHQ[60][20]; //qhat of heavy quark
0169   double qhat_over_T3;          //qhat/T^3 for heavy quark as fnc of (T,p)
0170   double D2piT;
0171
0172   // for heavy quark radiation table
0173   static const int HQener_gn = 500;
0174   static const int t_gn_1 = 100;
0175   static const int t_gn_2 = 125;
0176   static const int t_gn = t_gn_1 + t_gn_2;
0177   static const int temp_gn = 100;
0178
0179   static double dNg_over_dt_c[t_gn + 2][temp_gn + 1][HQener_gn + 1];
0180   static double dNg_over_dt_q[t_gn + 2][temp_gn + 1][HQener_gn + 1];
0181   static double dNg_over_dt_g[t_gn + 2][temp_gn + 1][HQener_gn + 1];
0182   static double max_dNgfnc_c[t_gn + 2][temp_gn + 1][HQener_gn + 1];
0183   static double max_dNgfnc_q[t_gn + 2][temp_gn + 1][HQener_gn + 1];
0184   static double max_dNgfnc_g[t_gn + 2][temp_gn + 1][HQener_gn + 1];
0185
0186   const double HQener_max = 1000.0;
0187   const double t_max_1 = 20.0;
0188   const double t_max_2 = 270.0;
0189   const double t_max = t_max_2;
0190   const double temp_max = 0.65;
0191   const double temp_min = 0.15;
0192   double delta_tg_1 = t_max_1 / t_gn_1;
0193   double delta_tg_2 = (t_max_2 - t_max_1) / t_gn_2;
0194   double delta_temp = (temp_max - temp_min) / temp_gn;
0195   double delta_HQener = HQener_max / HQener_gn;
0196
0197   // for MC initialization of jet partons
0198   static const int maxMC = 2000000;
0199   static double initMCX[maxMC], initMCY[maxMC];
0200
0201   //...radiation block
0202   int icl22;
0203   int icl23;  //the numerical switch in colljet23
0204   int iclrad; //the numerical switch in radiation
0205   int isp;    //the splitting function switch
0206
0207   //...global variable qhat
0208   int counth100 = 0;
0209
0210   double qhat; //transport parameter
0211
0212   double dng0[101][101] = {{0.0}}; //table of dn/dkperp2/dx
0213
0214   double Vtemp[4] = {0.0};
0215
0216   //...time system
0217
0218   static const int dimParList = 50; // originally 500000
0219
0220   double tirad[dimParList] = {0.0};
0221   double tiscatter[dimParList] = {0.0};
0222   double tiform[dimParList] = {0.0};   //pythia undone
0223   double Tint_lrf[dimParList] = {0.0}; //for heavy quark
0224   double eGluon = 0.0;
0225   double nGluon = 0.0;
0226   double dEel = 0.0;
0227
0228   //....radiated gluon
0229   double radng[dimParList] = {0.0};
0230   double xwm[3] = {0.0};
0231   double wkt2m;
0232
0233   double vf[4] = {0.0};    //flow velocity in tau-eta coordinate
0234   double vfcar[4] = {0.0}; //flow velocity in t-z coordinate
0235
0236   double vp[4] = {0.0};  //position of particle
0237   double vc0[4] = {0.0}; //flow velocity
0238
0239   //...dimensions in subroutine colljet and twcoll
0240   double vc[4] = {0.0};
0241   double pc0[4] = {0.0};
0242   double pc2[4] = {0.0};
0243   double pc3[4] = {0.0};
0244   double pc4[4] = {0.0};
0245   double p0[4] = {0.0};
0246   double p2[4] = {0.0};
0247   double p3[4] = {0.0};
0248   double p4[4] = {0.0};
0249
0250   double pc00[4] = {0.0};
0251   double pc30[4] = {0.0};
0252
0253   double pc01[4] = {0.0};
0254   double pb[4] = {0.0};
0255
0256   double Pj0[4] = {0.0};
0257
0258   double V[4][dimParList] = {{0.0}};  //parton position
0259   double P[7][dimParList] = {{0.0}};  //parton 4-momentum
0260   double V0[4][dimParList] = {{0.0}}; //negative parton position
0261   double P0[7][dimParList] = {{0.0}}; //negative parton 4-momentum
0262
0263   double Prad[4][dimParList] = {{0.0}};
0264
0265   double WT[dimParList] = {0};
0266   double WT0[dimParList] = {0};
0267
0268   int NR[dimParList] = {0};     //scattering rank
0269   int KATT1[dimParList] = {0};  //parton flavor
0270   int KATT10[dimParList] = {0}; //negative parton flavor
0271
0272   double PGm[4] = {0.0};
0273   double tjp[dimParList] = {0.0};
0274   double Vfrozen[4][dimParList] = {{0.0}};  //parton final 4 coordinate
0275   double Vfrozen0[4][dimParList] = {{0.0}}; //negative parton final 4 coordinate
0276   double Ecmcut = 2.0;                      //energy cut for free streaming
0277   double Tfrozen[dimParList] = {0.0};
0278   double Tfrozen0[dimParList] = {0.0};
0279   double vcfrozen[4][dimParList] = {{0.0}};
0280   double vcfrozen0[4][dimParList] = {{0.0}};
0281
0282   double VV[4][dimParList] = {{0.0}};
0283   double VV0[4][dimParList] = {{0.0}};
0284   double PP[4][dimParList] = {{0.0}};
0285   double PP0[4][dimParList] = {{0.0}};
0286   int CAT[dimParList] = {0};
0287   int CAT0[dimParList] = {0};
0288
0289   int ncut;
0290   int ncut0;
0291
0292   int n_coll22 = 0;
0293   int n_coll23 = 0;
0294   int ng_coll23 = 0;
0295   int ng_nrad = 0;
0296   int n_radiation = 0;
0297   int ng_radiation = 0;
0298   int n_gluon = 0;
0299   int n_sp1 = 0;
0300   int n_sp2 = 0;
0301
0302   // Variables for HQ 2->2
0303   static const int N_p1 = 500;
0304   static const int N_T = 60;
0305   static const int N_e2 = 75;
0306   static double distFncB[N_T][N_p1][N_e2], distFncF[N_T][N_p1][N_e2],
0307       distMaxB[N_T][N_p1][N_e2], distMaxF[N_T][N_p1][N_e2];
0308   static double distFncBM[N_T][N_p1], distFncFM[N_T][N_p1];
0309   double min_p1 = 0.0;
0310   double max_p1 = 1000.0;
0311   double bin_p1 = (max_p1 - min_p1) / N_p1;
0312   double min_T = 0.1;
0313   double max_T = 0.7;
0314   double bin_T = (max_T - min_T) / N_T;
0315   double min_e2 = 0.0;
0316   double max_e2 = 15.0;
0317   double bin_e2 = (max_e2 - min_e2) / N_e2;
0318
0319   //int loopN=10000;
0320   int loopN = 1000;
0321
0322   int numInitXY = 0;
0323   int flagScatter, flag_update, flag_update0;
0324   double Q00, Q0;
0325
0326   double tStart;// = 0.6;
0327
0328   //...functions
0329
0330   void LBT0(int &n, double &ti);
0331
0332   void trans(double v[4], double p[4]);
0333   void transback(double v[4], double p[4]);
0334
0335   float ran0(long *idum);
0336
0337   double alphas0(int &Kalphas, double temp0);
0338   double DebyeMass2(int &Kqhat0, double alphas, double temp0);
0339
0340   void lam(int KATT0, double &RTE, double E, double T, double &T1, double &T2,
0341            double &E1, double &E2, int &iT1, int &iT2, int &iE1, int &iE2);
0342
0343   void flavor(int &CT, int &KATT0, int &KATT2, int &KATT3, double RTE, double E,
0344               double T, double &T1, double &T2, double &E1, double &E2,
0345               int &iT1, int &iT2, int &iE1, int &iE2);
0346   void linear(int KATT, double E, double T, double &T1, double &T2, double &E1,
0347               double &E2, int &iT1, int &iT2, int &iE1, int &iE2, double &RTEg,
0348               double &RTEg1, double &RTEg2, double &RTEg3, double &RTEq,
0349               double &RTEq3, double &RTEq4, double &RTEq5, double &RTEq6,
0350               double &RTEq7, double &RTEq8, double &RTEHQ, double &RTEHQ11,
0351               double &RTEHQ12, double &qhatTP);
0352   void twflavor(int &CT, int &KATT0, int &KATT2, double E, double T);
0353   void twlinear(int KATT, double E, double T, double &RTEg1, double &RTEg2,
0354                 double &RTEq6, double &RTEq7, double &RTEq8);
0355   void colljet22(int CT, double temp, double qhat0ud, double v0[4],
0356                  double p0[4], double p2[4], double p3[4], double p4[4],
0357                  double &qt);
0358   void twcoll(int CT, double qhat0ud, double v0[4], double p0[4], double p2[4]);
0359
0360   void titau(double ti, double vf[4], double vp[4], double p0[4], double &Vx,
0361              double &Vy, double &Veta, double &Xtau);
0362
0363   void rotate(double px, double py, double pz, double p[4], int icc);
0364
0365   int KPoisson(double alambda);
0366
0367   void radiationHQ(int parID, double qhat0ud, double v0[4], double P2[4],
0368                    double P3[4], double P4[4], double Pj0[4], int &ic,
0369                    double Tdiff, double HQenergy, double max_Ng,
0370                    double temp_med, double xLow, double xInt);
0371   void collHQ22(int CT, double temp, double qhat, double v0[4], double p0[4],
0372                 double p2[4], double p3[4], double p4[4], double &qt);
0373   double Mqc2qc(double s, double t, double M);
0374   double Mgc2gc(double s, double t, double M);
0375
0376   void collHQ23(int parID, double temp_med, double qhat0ud, double v0[4],
0377                 double p0[4], double p2[4], double p3[4], double p4[4],
0378                 double qt, int &ic, double Tdiff, double HQenergy,
0379                 double max_Ng, double xLow, double xInt);
0380   double dNg_over_dxdydt(int parID, double x0g, double y0g, double HQenergy,
0381                          double HQmass, double temp_med, double Tdiff);
0382   double tau_f(double x0g, double y0g, double HQenergy, double HQmass);
0383   double splittingP(int parID, double z0g);
0384   double lambdas(double kTFnc);
0385   double nflavor(double kTFnc);
0386   double alphasHQ(double kTFnc, double tempFnc);
0387   double nHQgluon(int parID, double dtLRF, double &time_gluon, double &temp_med,
0388                   double &HQenergy, double &max_Ng);
0389
0390   void read_xyMC(int &numXY);
0391   void jetInitialize(int numXY);
0392   void setJetX(int numXY);
0393   void read_tables();
0394   void jetClean();
0395   void setParameter(string fileName);
0396   int checkParameter(int nArg);
0397
0398   //  extern "C" {
0399   //      void read_ccnu_(char *dataFN_in, int len1);
0400   //      void hydroinfoccnu_(double *Ct, double *Cx, double *Cy, double *Cz, double *Ctemp, double *Cvx, double *Cvy, double *Cvz, int *Cflag);
0401   //
0402   //      void sethydrofilesez_(int *dataID_in, char *dataFN_in, int *ctlID_in, char *ctlFN_in, int *bufferSize, int len1, int len2);
0403   //      void readhydroinfoshanshan_(double *t, double *x, double *y, double *z, double *e, double *s, double *temp, double *vx, double *vy, double *vz, int *flag);
0404   //  }
0405
0406   // Allows the registration of the module so that it is available to be used by the Jetscape framework.
0407   static RegisterJetScapeModule<LBT> reg;
0408 };
0409
0410 #endif // LBT_H