src/jet/ElossValidation.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 #include "ElossValidation.h"
0017 #include "JetScapeLogger.h"
0018 #include "JetScapeXML.h"
0019 #include <string>
0020
0021 #include "tinyxml2.h"
0022 #include <iostream>
0023
0024 #include "FluidDynamics.h"
0025
0026 #define MAGENTA "\033[35m"
0027
0028 using namespace Jetscape;
0029 using namespace std;
0030
0031 const double QS = 1.0;
0032
0033 ElossValidate::ElossValidate() {
0034   SetId("ElossValidate");
0035   VERBOSE(8);
0036 }
0037
0038 ElossValidate::~ElossValidate() { VERBOSE(8); }
0039
0040 void ElossValidate::Init() {
0041   JSINFO << "Initialize ElossValidate ...";
0042
0043   std::string s = GetXMLElementText({"Eloss", "ElossValidate", "name"});
0044   JSINFO << s << " to be initializied ...";
0045 }
0046
0047 void ElossValidate::WriteTask(weak_ptr<JetScapeWriter> w) {
0048   VERBOSE(8);
0049   auto f = w.lock();
0050   if (!f)
0051     return;
0052   f->WriteComment("ElossModule Parton List: " + GetId());
0053 }
0054
0055 void ElossValidate::DoEnergyLoss(double deltaT, double time, double Q2,
0056                                  vector<Parton> &pIn, vector<Parton> &pOut) {
0057
0058   VERBOSESHOWER(8) << MAGENTA << "SentInPartons Signal received : " << deltaT
0059                    << " " << Q2 << " " << &pIn;
0060
0061   // Check hydro communication
0062   std::unique_ptr<FluidCellInfo> check_fluid_info_ptr;
0063   GetHydroCellSignal(1, 1.0, 1.0, 0.0, check_fluid_info_ptr);
0064
0065   double delT = deltaT;
0066   double Time = time * fmToGeVinv;
0067   double deltaTime = delT * fmToGeVinv;
0068
0069   JSINFO << " the time in fm is " << time << " The time in GeV-1 is " << Time;
0070   JSINFO << " pid = " << pIn[0].pid() << " E = " << pIn[0].e()
0071          << " px = " << pIn[0].p(1) << " py = " << pIn[0].p(2)
0072          << "  pz = " << pIn[0].p(3) << " virtuality = " << pIn[0].t()
0073          << " form_time in fm = " << pIn[0].form_time() / fmToGeVinv;
0074   JSINFO << " color = " << pIn[0].color()
0075          << " anti-color = " << pIn[0].anti_color();
0076
0077   for (int i = 0; i < pIn.size(); i++) {
0078     TakeResponsibilityFor(
0079         pIn[i]); // Generate error if another module already has responsibility.
0080     JSINFO << " Parton Q2= " << pIn[i].t();
0081     JSINFO << " Parton Id= " << pIn[i].pid()
0082            << " and mass= " << pIn[i].restmass();
0083
0084     double velocity[4];
0085     velocity[0] = 1.0;
0086     for (int j = 1; j <= 3; j++) {
0087       velocity[j] = pIn[i].p(j) / pIn[i].e();
0088     }
0089
0090     if (pIn[i].form_time() <
0091         0.0) { /// A parton without a virtuality or formation time, must set...
0092       pIn[i].set_jet_v(velocity);
0093       pIn[i].set_t(QS * 2.);
0094       pIn[i].set_mean_form_time();
0095       JSINFO << " UPDATED Parton Q2= " << pIn[i].t();
0096     }
0097
0098     if (pIn[i].t() >= QS) {
0099       JSINFO << " ************ ";
0100       JSINFO << " DOING ELOSS  ";
0101       JSINFO << " ************ ";
0102
0103       // Split once
0104       // ----------
0105
0106       // daughter virtualities
0107       double tQd1 = QS - 0.00001;
0108       double tQd2 = QS - 0.00001;
0109
0110       // Always just radiate a gluon
0111       int d1_pid = pIn[i].pid();
0112       int d2_pid = gid;
0113
0114       double newp[4];
0115       double newx[4];
0116
0117       double z1 = 0.6;
0118       double z2 = 1.0 - z1;
0119
0120       newx[0] = Time + deltaTime;
0121       for (int j = 1; j <= 3; j++) {
0122         newx[j] =
0123             pIn[i].x_in().comp(j) + (Time + deltaTime - pIn[i].x_in().comp(0));
0124       }
0125
0126       // First daughter
0127       newp[0] = z1 * pIn[i].e();
0128       newp[1] = z1 * pIn[i].px() * (1.1);
0129       newp[2] = z1 * pIn[i].py() * (0.9);
0130       newp[3] = z1 * pIn[i].pz() * (1.1);
0131       pOut.push_back(Parton(0, d1_pid, 0, newp, newx));
0132       pOut.back().set_jet_v(pIn[i].jet_v());
0133       pOut.back().set_t(tQd1);
0134       pOut.back().set_mean_form_time();
0135       pOut.back().set_form_time(10);
0136
0137       // Second daughter
0138       newp[0] = z2 * pIn[i].e();
0139       newp[1] = z2 * pIn[i].px() * (0.9);
0140       newp[2] = z2 * pIn[i].py() * (1.1);
0141       newp[3] = z2 * pIn[i].pz() * (0.8);
0142       pOut.push_back(Parton(0, d2_pid, 0, newp, newx));
0143       pOut.back().set_jet_v(pIn[i].jet_v());
0144       pOut.back().set_t(tQd2);
0145       pOut.back().set_mean_form_time();
0146       pOut.back().set_form_time(10);
0147     }
0148   }
0149 }