Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Martini.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Martini.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 MARTINI_H
17 #define MARTINI_H
18 
19 #include <fstream>
20 #include <math.h>
21 #include "JetEnergyLossModule.h"
22 #include "JetScapeConstants.h"
23 
24 using namespace Jetscape;
25 
26 class MARTINIUserInfo : public fjcore::PseudoJet::UserInfoBase {
27 private:
28  bool _coherent{false};
29  int _sibling{-1};
30  FourVector _pAtSplit{FourVector(0., 0., 0., 0.)};
31 
32 public:
33  MARTINIUserInfo() {}
34  MARTINIUserInfo(bool coherent, int sibling, FourVector pAtSplit)
35  : _coherent(coherent), _sibling(sibling), _pAtSplit(pAtSplit){};
36 
37  bool coherent() const { return _coherent; }
38  int coherent_with() const { return _sibling; }
39  FourVector p_at_split() const { return _pAtSplit; }
40 };
41 
42 //Basic.h//
43 struct RateRadiative {
44  double qqg;
45  double ggg;
46  double gqq;
47  double qqgamma;
48 };
49 
50 struct RateElastic {
51  double qq;
52  double gq;
53  double qg;
54  double gg;
55 };
56 
57 struct RateConversion {
58  double qg;
59  double gq;
60  double qgamma;
61 };
62 
63 class Martini : public JetEnergyLossModule<
64  Martini> //, public std::enable_shared_from_this<Martini>
65 {
66 private:
67  // AMY rates are calculated in p/T > AMYpCut
68  static constexpr double AMYpCut = 4.01;
69  // minimum energy of parton that can lose energy
70  const double eLossCut = 1.0;
71 
72  double Q0; // Separation scale between Matter and Martini
73  double alpha_s0, alpha_s;
74  double alpha_em;
75  double g;
76  double pcut; // below this scale, no further Eloss
77  double hydro_Tc; // critical temperature
78  double tStart; // initilization time of hydro
79  int run_alphas; // running alpha_s or not
80  int recoil_on; // turn on recoil
81 
82  //Import.h//
83  static const int NP = 230;
84  static const int NK = 381;
85 
86  static const int Nalphas = 11;
87  static const int Nomega = 120;
88  static const int Nq = 60;
89 
90  static constexpr double omegaStep = 0.2;
91  static constexpr double qStep = 0.2;
92  static constexpr double alphaMin = 0.15;
93  static constexpr double alphaStep = 0.03;
94 
95  typedef struct {
96  double ddf;
97  double dda;
98  double dcf;
99  double dca;
100  int include_gluons;
101  int Nc;
102  int Nf;
103  int BetheHeitler;
104  int BDMPS;
105  double alpha;
106  double delta_x;
107  double dx_max;
108  double dp;
109  double p_min;
110  double p_max;
111  long n_p;
112  long n_pmin;
113  double k_min;
114  double k_max;
115  long n_k;
116  long n_kmin;
117  } Gamma_info;
118 
119  // Structure to store information about splitting functions
120  typedef struct {
121  double qqg[NP][NK];
122  double gqq[NP][NK];
123  double ggg[NP][NK];
124  double qqgamma[NP][NK];
125 
126  double tau_qqg[NP][NK];
127  double tau_gqq[NP][NK];
128  double tau_ggg[NP][NK];
129  double tau_qqgamma[NP][NK];
130  } dGammas;
131 
132  Gamma_info dat;
133  dGammas Gam;
134 
135  vector<double> *dGamma_qq;
136  vector<double> *dGamma_qg;
137  vector<double> *dGamma_qq_q;
138  vector<double> *dGamma_qg_q;
139 
140  static int pLabelNew;
141 
142  // Allows the registration of the module so that it is available to be used by the Jetscape framework.
143  static RegisterJetScapeModule<Martini> reg;
144 
145  double RunningAlphaS(double muSquare);
146 
147 public:
148  Martini();
149  virtual ~Martini();
150 
151  //main//
152  void Init();
153  void DoEnergyLoss(double deltaT, double Time, double Q2, vector<Parton> &pIn,
154  vector<Parton> &pOut);
155  int DetermineProcess(double p, double T, double deltaTRest, int id);
156  void WriteTask(weak_ptr<JetScapeWriter> w){};
157 
158  //Radiative//
159  RateRadiative getRateRadTotal(double p, double T);
160  RateRadiative getRateRadPos(double u, double T);
161  RateRadiative getRateRadNeg(double u, double T);
162 
163  double getNewMomentumRad(double p, double T, int process);
164  double area(double y, double u, int posNegSwitch, int process);
165  double function(double u, double y, int process);
166 
167  bool isCoherent(Parton &pIn, int sibling, double T);
168 
169  //Elastic//
170  RateElastic getRateElasTotal(double p, double T);
171  RateElastic getRateElasPos(double u, double T);
172  RateElastic getRateElasNeg(double u, double T);
173 
174  RateConversion getRateConv(double p, double T);
175 
176  double getEnergyTransfer(double p, double T, int process);
177  double getMomentumTransfer(double p, double omega, double T, int process);
178 
179  double areaOmega(double u, int posNegSwitch, int process);
180  double areaQ(double u, double omega, int process);
181  double functionOmega(double u, double y, int process) {
182  return getElasticRateOmega(u, y, process);
183  }
184  double functionQ(double u, double omega, double q, int process) {
185  return getElasticRateQ(u, omega, q, process);
186  }
187  FourVector getNewMomentumElas(FourVector pVec, double omega, double q);
188  FourVector getThermalVec(FourVector qVec, double T, int id);
189  double getThermal(double k_min, double T, int kind);
190 
191  //Rate table//
192  void readRadiativeRate(Gamma_info *dat, dGammas *Gam);
193  void readElasticRateOmega();
194  void readElasticRateQ();
195 
196  double getRate_qqg(double p, double k);
197  double getRate_gqq(double p, double k);
198  double getRate_ggg(double p, double k);
199  double getRate_qqgamma(double p, double k);
200  double use_table(double p, double k, double dGamma[NP][NK], int which_kind);
201 
202  double getElasticRateOmega(double u, double omega, int process);
203  double getElasticRateQ(double u, double omega, double q, int process);
204  double use_elastic_table_omega(double omega, int which_kind);
205  double use_elastic_table_q(double u, double omega, int which_kind);
206 
207  static void IncrementpLable() { pLabelNew++; }
208 
209 protected:
210  uniform_real_distribution<double> ZeroOneDistribution;
211 
212  string PathToTables;
213 };
214 
215 double LambertW(double z);
216 
217 #endif // MARTINI_H