src/jet/Martini.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 MARTINI_H
0017 #define MARTINI_H
0018
0019 #include <fstream>
0020 #include <math.h>
0021 #include "JetEnergyLossModule.h"
0022 #include "JetScapeConstants.h"
0023
0024 using namespace Jetscape;
0025
0026 class MARTINIUserInfo : public fjcore::PseudoJet::UserInfoBase {
0027 private:
0028   bool _coherent{false};
0029   int _sibling{-1};
0030   FourVector _pAtSplit{FourVector(0., 0., 0., 0.)};
0031
0032 public:
0033   MARTINIUserInfo() {}
0034   MARTINIUserInfo(bool coherent, int sibling, FourVector pAtSplit)
0035       : _coherent(coherent), _sibling(sibling), _pAtSplit(pAtSplit){};
0036
0037   bool coherent() const { return _coherent; }
0038   int coherent_with() const { return _sibling; }
0039   FourVector p_at_split() const { return _pAtSplit; }
0040 };
0041
0042 //Basic.h//
0043 struct RateRadiative {
0044   double qqg;
0045   double ggg;
0046   double gqq;
0047   double qqgamma;
0048 };
0049
0050 struct RateElastic {
0051   double qq;
0052   double gq;
0053   double qg;
0054   double gg;
0055 };
0056
0057 struct RateConversion {
0058   double qg;
0059   double gq;
0060   double qgamma;
0061 };
0062
0063 class Martini : public JetEnergyLossModule<
0064                     Martini> //, public std::enable_shared_from_this<Martini>
0065 {
0066 private:
0067   // AMY rates are calculated in p/T > AMYpCut
0068   static constexpr double AMYpCut = 4.01;
0069   // minimum energy of parton that can lose energy
0070   const double eLossCut = 1.0;
0071
0072   double Q0; // Separation scale between Matter and Martini
0073   double alpha_s0, alpha_s;
0074   double alpha_em;
0075   double g;
0076   double pcut;         // below this scale, no further Eloss
0077   double hydro_Tc;     // critical temperature
0078   double tStart; // initilization time of hydro
0079   int run_alphas;      // running alpha_s or not
0080   int recoil_on;       // turn on recoil
0081
0082   //Import.h//
0083   static const int NP = 230;
0084   static const int NK = 381;
0085
0086   static const int Nalphas = 11;
0087   static const int Nomega = 120;
0088   static const int Nq = 60;
0089
0090   static constexpr double omegaStep = 0.2;
0091   static constexpr double qStep = 0.2;
0092   static constexpr double alphaMin = 0.15;
0093   static constexpr double alphaStep = 0.03;
0094
0095   typedef struct {
0096     double ddf;
0097     double dda;
0098     double dcf;
0099     double dca;
0100     int include_gluons;
0101     int Nc;
0102     int Nf;
0103     int BetheHeitler;
0104     int BDMPS;
0105     double alpha;
0106     double delta_x;
0107     double dx_max;
0108     double dp;
0109     double p_min;
0110     double p_max;
0111     long n_p;
0112     long n_pmin;
0113     double k_min;
0114     double k_max;
0115     long n_k;
0116     long n_kmin;
0117   } Gamma_info;
0118
0119   // Structure to store information about splitting functions
0120   typedef struct {
0121     double qqg[NP][NK];
0122     double gqq[NP][NK];
0123     double ggg[NP][NK];
0124     double qqgamma[NP][NK];
0125
0126     double tau_qqg[NP][NK];
0127     double tau_gqq[NP][NK];
0128     double tau_ggg[NP][NK];
0129     double tau_qqgamma[NP][NK];
0130   } dGammas;
0131
0132   Gamma_info dat;
0133   dGammas Gam;
0134
0135   vector<double> *dGamma_qq;
0136   vector<double> *dGamma_qg;
0137   vector<double> *dGamma_qq_q;
0138   vector<double> *dGamma_qg_q;
0139
0140   static int pLabelNew;
0141
0142   // Allows the registration of the module so that it is available to be used by the Jetscape framework.
0143   static RegisterJetScapeModule<Martini> reg;
0144
0145   double RunningAlphaS(double muSquare);
0146
0147 public:
0148   Martini();
0149   virtual ~Martini();
0150
0151   //main//
0152   void Init();
0153   void DoEnergyLoss(double deltaT, double Time, double Q2, vector<Parton> &pIn,
0154                     vector<Parton> &pOut);
0155   int DetermineProcess(double p, double T, double deltaTRest, int id);
0156   void WriteTask(weak_ptr<JetScapeWriter> w){};
0157
0158   //Radiative//
0159   RateRadiative getRateRadTotal(double p, double T);
0160   RateRadiative getRateRadPos(double u, double T);
0161   RateRadiative getRateRadNeg(double u, double T);
0162
0163   double getNewMomentumRad(double p, double T, int process);
0164   double area(double y, double u, int posNegSwitch, int process);
0165   double function(double u, double y, int process);
0166
0167   bool isCoherent(Parton &pIn, int sibling, double T);
0168
0169   //Elastic//
0170   RateElastic getRateElasTotal(double p, double T);
0171   RateElastic getRateElasPos(double u, double T);
0172   RateElastic getRateElasNeg(double u, double T);
0173
0174   RateConversion getRateConv(double p, double T);
0175
0176   double getEnergyTransfer(double p, double T, int process);
0177   double getMomentumTransfer(double p, double omega, double T, int process);
0178
0179   double areaOmega(double u, int posNegSwitch, int process);
0180   double areaQ(double u, double omega, int process);
0181   double functionOmega(double u, double y, int process) {
0182     return getElasticRateOmega(u, y, process);
0183   }
0184   double functionQ(double u, double omega, double q, int process) {
0185     return getElasticRateQ(u, omega, q, process);
0186   }
0187   FourVector getNewMomentumElas(FourVector pVec, double omega, double q);
0188   FourVector getThermalVec(FourVector qVec, double T, int id);
0189   double getThermal(double k_min, double T, int kind);
0190
0191   //Rate table//
0192   void readRadiativeRate(Gamma_info *dat, dGammas *Gam);
0193   void readElasticRateOmega();
0194   void readElasticRateQ();
0195
0196   double getRate_qqg(double p, double k);
0197   double getRate_gqq(double p, double k);
0198   double getRate_ggg(double p, double k);
0199   double getRate_qqgamma(double p, double k);
0200   double use_table(double p, double k, double dGamma[NP][NK], int which_kind);
0201
0202   double getElasticRateOmega(double u, double omega, int process);
0203   double getElasticRateQ(double u, double omega, double q, int process);
0204   double use_elastic_table_omega(double omega, int which_kind);
0205   double use_elastic_table_q(double u, double omega, int which_kind);
0206
0207   static void IncrementpLable() { pLabelNew++; }
0208
0209 protected:
0210   uniform_real_distribution<double> ZeroOneDistribution;
0211
0212   string PathToTables;
0213 };
0214
0215 double LambertW(double z);
0216
0217 #endif // MARTINI_H