src/jet/AdSCFT.cc

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
0017 #include "AdSCFT.h"
0018 #include "JetScapeLogger.h"
0019 #include "JetScapeXML.h"
0020 #include <string>
0021
0022 #include "tinyxml2.h"
0023 #include <iostream>
0024 #include <fstream>
0025
0026 #include "FluidDynamics.h"
0027 #include "AdSCFTMutex.h"
0028 #define MAGENTA "\033[35m"
0029
0030 // Register the module with the base class
0031 RegisterJetScapeModule<AdSCFT> AdSCFT::reg("AdSCFT");
0032
0033 AdSCFT::AdSCFT() {
0034   SetId("AdSCFT");
0035   kappa = -99;
0036   T0 = -99;
0037   Q0 = -99;
0038
0039   VERBOSE(8);
0040
0041   // create and set Martini Mutex
0042   auto adscft_mutex = make_shared<AdSCFTMutex>();
0043   SetMutex(adscft_mutex);
0044 }
0045
0046 AdSCFT::~AdSCFT() { VERBOSE(8); }
0047
0048 void AdSCFT::Init() {
0049   JSINFO << "Initialize AdSCFT ...";
0050
0051   std::string s = GetXMLElementText({"Eloss", "AdSCFT", "name"});
0052   JSDEBUG << s << " to be initialized ...";
0053
0054   //Kappa
0055   kappa = GetXMLElementDouble({"Eloss", "AdSCFT", "kappa"});
0056   JSINFO << "AdSCFT kappa = " << kappa;
0057
0058   //T0 [GeV]
0059   T0 = GetXMLElementDouble({"Eloss", "AdSCFT", "T0"});
0060   JSINFO << "AdSCFT T0 = " << T0;
0061
0062   //Q0 [GeV]
0063   Q0 = GetXMLElementDouble({"Eloss", "AdSCFT", "Q0"});
0064   JSINFO << "AdSCFT Q0 = " << Q0;
0065
0066   //Vac or Med
0067   in_vac = GetXMLElementDouble({"Eloss", "AdSCFT", "in_vac"});
0068   ;
0069 }
0070
0071 void AdSCFT::WriteTask(weak_ptr<JetScapeWriter> w) {
0072   VERBOSE(8);
0073   auto f = w.lock();
0074   if (!f)
0075     return;
0076   f->WriteComment("ElossModule Parton List: " + GetId());
0077   f->WriteComment("Energy loss to be implemented accordingly ...");
0078 }
0079
0080 void AdSCFT::DoEnergyLoss(double deltaT, double time, double Q2,
0081                           vector<Parton> &pIn, vector<Parton> &pOut) {
0082
0083   VERBOSESHOWER(8) << MAGENTA << "SentInPartons Signal received : " << deltaT
0084                    << " " << Q2 << " " << &pIn;
0085
0086   for (int i = 0; i < pIn.size(); i++) {
0087     //Skip photons
0088     if (pIn[i].pid() == photonid) {
0089       pOut.push_back(pIn[i]);
0090       return;
0091     }
0092
0093     JSDEBUG << " in AdS/CFT";
0094     JSDEBUG << " Parton Q2= " << pIn[i].t();
0095     JSDEBUG << " Parton Id= " << pIn[i].pid()
0096             << " and mass= " << pIn[i].restmass();
0097
0098     //Parton 4-momentum
0099     double p[4];
0100     p[0] = pIn[i].px();
0101     p[1] = pIn[i].py();
0102     p[2] = pIn[i].pz();
0103     p[3] = pIn[i].e();
0104     double pmod = std::sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]);
0105
0106     //Parton velocity
0107     vector<double> w;
0108     for (unsigned int j = 0; j < 4; j++)
0109       w.push_back(p[j] / p[3]);
0110     double w2 = std::pow(w[0], 2.) + std::pow(w[1], 2.) + std::pow(w[2], 2.);
0111     for (unsigned int j = 0; j < 3; j++)
0112       w[j] /= std::sqrt(w2);
0113
0114     //Parton 4-position
0115     double initR0 = pIn[i].x_in().t(); //Time when the parton was last modified
0116     double x[4];
0117     //x[0] = pIn[i].x_in().x() + (time - initR0) * w[0] / std::sqrt(w2);
0118     //x[1] = pIn[i].x_in().y() + (time - initR0) * w[1] / std::sqrt(w2);
0119     //x[2] = pIn[i].x_in().z() + (time - initR0) * w[2] / std::sqrt(w2);
0120     x[0] = pIn[i].x_in().x() + (time - initR0) * w[0];
0121     x[1] = pIn[i].x_in().y() + (time - initR0) * w[1];
0122     x[2] = pIn[i].x_in().z() + (time - initR0) * w[2];
0123     x[3] = pIn[i].x_in().t() + (time - initR0) * w[3];
0124
0125     //Extract fluid properties
0126     std::unique_ptr<FluidCellInfo> check_fluid_info_ptr;
0127     double tau = std::sqrt(x[3] * x[3] - x[2] * x[2]);
0128     double temp, vx, vy, vz;
0129     //Only get a temp!=0 if in_vac=0
0130     //if (x[3]>=tStart && !in_vac) //Can use time if there is preequilibrium eloss
0131     if (tau >= tStart &&
0132         !in_vac) //Should use tau, not t, in absence of preequilibrium eloss
0133     {
0134       GetHydroCellSignal(x[3], x[0], x[1], x[2], check_fluid_info_ptr);
0135       if (!GetJetSignalConnected()) {
0136         JSWARN << "Couldn't find a hydro module attached!";
0137         throw std::runtime_error(
0138             "Please attach a hydro module or set in_vac to 1 in the XML file");
0139       }
0140
0141       VERBOSE(8) << MAGENTA << "Temperature from Brick (Signal) = "
0142                  << check_fluid_info_ptr->temperature;
0143
0144       temp = check_fluid_info_ptr->temperature;
0145       vx = check_fluid_info_ptr->vx;
0146       vy = check_fluid_info_ptr->vy;
0147       vz = check_fluid_info_ptr->vz;
0148     } else
0149       temp = 0., vx = 0., vy = 0., vz = 0.;
0150     JSDEBUG << " system time= " << time << " parton time= " << x[3]
0151             << " tau= " << tau << " temp= " << temp << " vz= " << w[2];
0152     JSDEBUG << " energy= " << pIn[i].e();
0153
0154     //Do eloss in AdS/CFT if:
0155     // *Virtuality below Qcut ( QS[GeV^2] , NOT taken from XML yet )
0156     // *Fluid temperature above Tcut ( T0 from XML )
0157     // *Parton is not completely quenched ( Ecut = 0.00001 )
0158     double QS = Q0 * Q0;
0159     if (pIn[i].t() <= QS + rounding_error && temp >= T0 && pmod > 0.00001 &&
0160         pIn[i].pstat() >= 0) {
0161       //cout << " ADS Q= " << pIn[i].t() << " Q0= " << Q0 << " temp= " << temp << " T0= " << T0 << endl;
0162       //cout << " ADS tau= " << tau << " x= " << x[0] << " y= " << x[1] << " z= " << x[2] << " t= " << x[3] << endl;
0163       TakeResponsibilityFor(
0164           pIn[i]); // Generate error if another module already has responsibility.
0165
0166       VERBOSE(8) << " ************ \n \n";
0167       VERBOSE(8) << " DOING ADSCFT \n \n";
0168       VERBOSE(8) << " ************ \n \n";
0169
0170       //Fluid quantities
0171       vector<double> v;
0172       v.push_back(vx), v.push_back(vy), v.push_back(vz), v.push_back(1.);
0173       double v2 = std::pow(v[0], 2.) + std::pow(v[1], 2.) + std::pow(v[2], 2.);
0174       double lore = 1. / sqrt(1. - v2); //Gamma Factor
0175
0176       //Color Factor (wrt quark)
0177       double CF;
0178       if (pIn[i].pid() == 21)
0179         CF = std::pow(9. / 4., 1. / 3.); //Gluon
0180       else
0181         CF = 1.; //Quark
0182
0183       //Energy (three-momentum) of parton as it entered this module for the first time
0184       double ei = pmod;
0185       double l_dist = 0., f_dist = 0.;
0186       if (pIn[i].has_user_info<AdSCFTUserInfo>()) {
0187         ei = pIn[i].user_info<AdSCFTUserInfo>().part_ei();
0188         l_dist = pIn[i].user_info<AdSCFTUserInfo>().l_dist();
0189         f_dist = pIn[i].user_info<AdSCFTUserInfo>().f_dist();
0190       } else {
0191         pIn[i].set_user_info(new AdSCFTUserInfo(ei, f_dist, l_dist));
0192       }
0193
0194       //JSDEBUG << " ei= " << ei;
0195       //JSDEBUG << " px= " << p[0] << " py= " << p[1] << " pz= " << p[2] << " en= " << p[3];
0196       //JSDEBUG << " x= " << x[0] << " y= " << x[1] << " z= " << x[2] << " t= " << x[3];
0197
0198       double restmass = pIn[i].restmass();
0199       if (abs(pIn[i].pid()) <= 3)
0200         restmass = 0.;
0201       double virttwo = p[3] * p[3] - p[0] * p[0] - p[1] * p[1] - p[2] * p[2] -
0202                        restmass * restmass;
0203       double virt = std::sqrt(virttwo);
0204       //JSDEBUG << " virt= " << virt;
0205
0206       //Needed for boosts (v.w)
0207       double vscalw = v[0] * w[0] + v[1] * w[1] + v[2] * w[2];
0208
0209       //Distance travelled in LAB frame
0210       l_dist += deltaT;
0211
0212       //Distance travelled in FRF - accumulating steps from previous, different fluid cells
0213       double insqrt = w2 + lore * lore * (v2 - 2. * vscalw + vscalw * vscalw);
0214       if (insqrt <= 0.)
0215         insqrt = 0.;
0216       f_dist += deltaT * std::sqrt(insqrt);
0217
0218       //JSDEBUG << " l_dist= " << l_dist << " f_dist= " << f_dist;
0219       //Initial energy of the parton in the FRF
0220       double Efs = ei * lore * (1. - vscalw);
0221
0222       double newEn = pmod;
0223       if (temp >= 0.)
0224         newEn = pmod - AdSCFT::Drag(f_dist, deltaT, Efs, temp, CF);
0225       if (newEn < 0.)
0226         newEn = pmod / 10000000.;
0227       double lambda = newEn / pmod;
0228       //JSDEBUG << " lambda= " << lambda;
0229       //JSDEBUG << " Elost= " << pmod-newEn;
0230
0231       //Update 4-momentum (don't modify mass)
0232       for (unsigned a = 0; a < 3; a++)
0233         p[a] *= lambda;
0234       virt *= lambda;
0235       p[3] = std::sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2] + virt * virt +
0236                        restmass * restmass);
0237       virttwo = p[3] * p[3] - p[0] * p[0] - p[1] * p[1] - p[2] * p[2] -
0238                 restmass * restmass;
0239       virt = std::sqrt(virttwo);
0240
0241       //DON'T Update 4-position here, already done at beginning
0242       //for (unsigned a=0; a<4; a++) x[a]+=w[a]*deltaT;
0243
0244       double fx[4];
0245       fx[0] = x[3], fx[1] = x[0], fx[2] = x[1], fx[3] = x[2];
0246       //Feed into parton list
0247       int pLabel = pIn[i].plabel();
0248       int Id = pIn[i].pid();
0249       int pStat = pIn[i].pstat();
0250       //cout << " pStat= " << pStat << endl;
0251       FourVector pVec(p[0], p[1], p[2], p[3]);
0252       FourVector xVec;
0253       pOut.push_back(Parton(pLabel, Id, pStat, pVec, xVec));
0254       pOut[pOut.size() - 1].set_x(fx);
0255       pOut.back().set_user_info(new AdSCFTUserInfo(ei, f_dist, l_dist));
0256
0257       //Copy variables needed in case parton returns to MATTER in future steps
0258       double velocity_jet[4];
0259       velocity_jet[0] = 1.0;
0260       velocity_jet[1] = pIn[i].jet_v().x();
0261       velocity_jet[2] = pIn[i].jet_v().y();
0262       velocity_jet[3] = pIn[i].jet_v().z();
0263       pOut[pOut.size() - 1].set_jet_v(velocity_jet);
0264       pOut[pOut.size() - 1].set_mean_form_time();
0265       double ft = pOut[pOut.size() - 1].mean_form_time();
0266       pOut[pOut.size() - 1].set_form_time(ft);
0267
0268       //Add missing momentum
0269       FourVector pVecM(pIn[i].px() - p[0], pIn[i].py() - p[1],
0270                        pIn[i].pz() - p[2], pIn[i].e() - p[3]);
0271       pOut.push_back(Parton(0, 21, -13, pVecM, xVec));
0272       pOut[pOut.size() - 1].set_x(fx);
0273
0274     } //End if do-eloss
0275
0276   } //End pIn loop
0277 }
0278
0279 double AdSCFT::Drag(double f_dist, double deltaT, double Efs, double temp,
0280                     double CF) {
0281   if (kappa != 0.) {
0282     double tstop = 0.2 * std::pow(Efs, 1. / 3.) /
0283                    (2. * std::pow(temp, 4. / 3.) * kappa) /
0284                    CF; //Stopping distance in fm
0285     double beta =
0286         tstop / f_dist; //Ratio of stopping distance to distance travelled
0287     if (beta > 1.)      //If did not get to stopping distance
0288     {
0289       double intpiece = Efs * deltaT * 4. / (3.141592) *
0290                         (1. / (beta * tstop * std::sqrt(beta * beta - 1.)));
0291       return intpiece;
0292     } else //If reached stopping distance, subtract all energy
0293     {
0294       return 100000.;
0295     }
0296   } else
0297     return 0.;
0298 }
0299
0300 void AdSCFT::Clear() {}