src/framework/JetEnergyLoss.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 <iostream>
0017 #include <thread>
0018 //#include <mutex>
0019 //#include <condition_variable>
0020 //#include <future>
0021
0022 #include "JetEnergyLoss.h"
0023 #include "JetScapeLogger.h"
0024 #include "JetScapeXML.h"
0025 #include <string>
0026 #include "tinyxml2.h"
0027 #include "JetScapeSignalManager.h"
0028 #include "JetScapeWriterStream.h"
0029 #include "HardProcess.h"
0030 #include "JetScapeModuleMutex.h"
0031 #include "LiquefierBase.h"
0032 #include "MakeUniqueHelper.h"
0033 #include "FluidDynamics.h"
0034 #include <GTL/dfs.h>
0035
0036 #ifdef USE_HEPMC
0037 #include "JetScapeWriterHepMC.h"
0038 #endif
0039
0040 #define BOLDCYAN "\033[1m\033[36m" /* Bold Cyan */
0041
0042 using namespace std;
0043
0044 namespace Jetscape {
0045
0046 JetEnergyLoss::JetEnergyLoss() {
0047   qhat = -99.99;
0048   SetId("JetEnergyLoss");
0049   jetSignalConnected = false;
0050   edensitySignalConnected = false;
0051   GetHydroCellSignalConnected = false;
0052   GetHydroTau0SignalConnected = false;
0053   SentInPartonsConnected = false;
0054
0055   deltaT = 0;
0056   maxT = 0;
0057
0058   inP = nullptr;
0059   pShower = nullptr;
0060
0061   VERBOSE(8);
0062 }
0063
0064 JetEnergyLoss::~JetEnergyLoss() {
0065   VERBOSE(8);
0066
0067   //pShower->clear();
0068   disconnect_all();
0069 }
0070
0071 JetEnergyLoss::JetEnergyLoss(const JetEnergyLoss &j) {
0072   qhat = j.GetQhat();
0073   SetActive(j.GetActive());
0074   SetId(j.GetId());
0075   SetJetSignalConnected(false);
0076   SetEdensitySignalConnected(false);
0077   SetGetHydroCellSignalConnected(false);
0078   SetGetHydroTau0SignalConnected(false);
0079   SetSentInPartonsConnected(false);
0080
0081   deltaT = j.deltaT;
0082   maxT = j.maxT;
0083
0084   inP = nullptr;
0085   pShower = nullptr;
0086
0087   VERBOSE(8) << "To be copied : # Subtasks = " << j.GetTaskList().size();
0088   for (auto it : j.GetTaskList()) {
0089     // Working via CRTP JetEnergyLossModule Clone function !
0090     auto st = dynamic_pointer_cast<JetEnergyLoss>(it)
0091                   ->Clone(); //shared ptr with clone !!????
0092     Add(st);
0093   }
0094 }
0095
0096 void JetEnergyLoss::Clear() {
0097   VERBOSESHOWER(8);
0098   if (pShower)
0099     pShower->clear();
0100
0101   this->final_Partons.clear();
0102
0103   inP = nullptr;
0104 }
0105
0106 void JetEnergyLoss::Init() {
0107   JetScapeModuleBase::Init();
0108
0109   JSINFO << "Initialize JetEnergyLoss ...";
0110
0111   deltaT = GetXMLElementDouble({"Eloss", "deltaT"});
0112
0113   maxT = GetXMLElementDouble({"Eloss", "maxT"});
0114   JSINFO << "Eloss shower with deltaT = " << deltaT << " and maxT = " << maxT;
0115
0116   std::string mutexOnString = GetXMLElementText({"Eloss", "mutex"}, false);
0117   if (!mutexOnString.compare("ON"))
0118   //Check mutual exclusion of Eloss Modules
0119   {
0120     if (GetNumberOfTasks() > 1) {
0121       for (auto elossModule : GetTaskList()) {
0122         shared_ptr<JetScapeModuleMutex> mutex_ptr = elossModule->GetMutex();
0123         if (mutex_ptr) {
0124           if (!(mutex_ptr->CheckMutex(GetTaskList()))) {
0125             JSWARN << "Mutually exclusive Energy-Loss modules attached together!";
0126             throw std::runtime_error("Fix it by attaching one of them.");
0127           }
0128         }
0129       }
0130     }
0131   }
0132
0133   if (GetNumberOfTasks() < 1) {
0134     JSWARN << " : No valid Energy Loss modules found ...";
0135     exit(-1);
0136   }
0137
0138   inP = nullptr;
0139   pShower = nullptr;
0140
0141   JSINFO << "Found " << GetNumberOfTasks()
0142          << " Eloss Tasks/Modules Initialize them ... ";
0143
0144   JetScapeTask::InitTasks();
0145 }
0146
0147 void JetEnergyLoss::DoShower() {
0148   double tStart = 0;
0149   double currentTime = 0;
0150
0151   VERBOSESHOWER(8) << "Hard Parton from Initial Hard Process ...";
0152   VERBOSEPARTON(6, *GetShowerInitiatingParton());
0153
0154   // consider pointers for speed up ...
0155   vector<Parton> pIn;
0156   // DEBUG this guy isn't linked to anything - put in test particle for now
0157   pIn.push_back(*GetShowerInitiatingParton());
0158
0159   vector<node> vStartVec;
0160   // Add here the Hard Shower emitting parton ...
0161   vStart = pShower->new_vertex(make_shared<Vertex>());
0162   vEnd = pShower->new_vertex(make_shared<Vertex>());
0163   // Add original parton later, after it had a chance to acquire virtuality
0164   // pShower->new_parton(vStart,vEnd,make_shared<Parton>(*GetShowerInitiatingParton()));
0165
0166   // start then the recursive shower ...
0167   vStartVec.push_back(vEnd);
0168
0169   // cerr << " ---------------------------------------------- " << endl;
0170   // cerr << "Start with " << *GetShowerInitiatingParton()
0171   //      << "  -> " << GetShowerInitiatingParton()->t() << endl;
0172   bool foundchangedorig = false;
0173   int droplet_stat = -11;
0174   int miss_stat = -13;
0175   int neg_stat = -17;
0176   if (!weak_ptr_is_uninitialized(liquefier_ptr)) {
0177     droplet_stat = liquefier_ptr.lock()->get_drop_stat();
0178     miss_stat = liquefier_ptr.lock()->get_miss_stat();
0179     neg_stat = liquefier_ptr.lock()->get_neg_stat();
0180   }
0181   do {
0182     vector<Parton> pOut;
0183     vector<Parton> pInTemp;
0184
0185     vector<node> vStartVecOut;
0186     vector<node> vStartVecTemp;
0187
0188     VERBOSESHOWER(7) << "Current time = " << currentTime << " with #Input "
0189                      << pIn.size();
0190     currentTime += deltaT;
0191
0192     for (int i = 0; i < pIn.size(); i++) {
0193       vector<Parton> pInTempModule;
0194       vector<Parton> pOutTemp;
0195       // JSINFO << pIn.at(i).edgeid();
0196       pInTempModule.push_back(pIn[i]);
0197       SentInPartons(deltaT, currentTime, pIn[i].pt(), pInTempModule, pOutTemp);
0198
0199       // apply liquefier
0200       if (!weak_ptr_is_uninitialized(liquefier_ptr)) {
0201         liquefier_ptr.lock()->add_hydro_sources(pInTempModule, pOutTemp);
0202       }
0203
0204       // stuffs related to vertex
0205       if (!foundchangedorig) {
0206         // cerr << " End with "<< pInTempModule.at(0) << "  -> "
0207         //      << pInTempModule.at(0).t() << endl;
0208         // cerr << " ---------------------------------------------- "
0209         //      << endl;
0210         pShower->new_parton(vStart, vEnd,
0211                             make_shared<Parton>(pInTempModule.at(0)));
0212         foundchangedorig = true;
0213       }
0214
0215       vStart = vStartVec[i];
0216       if (pOutTemp.size() == 0) {
0217         // no need to generate a vStart for photons and liquefied
0218         // partons
0219         if (pInTempModule[0].pstat() != droplet_stat &&
0220             pInTempModule[0].pstat() != miss_stat &&
0221             pInTempModule[0].pstat() != neg_stat &&
0222             !pInTempModule[0].isPhoton(pInTempModule[0].pid())) {
0223           vStartVecTemp.push_back(vStart);
0224         }
0225       } else if (pOutTemp.size() == 1) {
0226         // no need to generate a vStart for photons and liquefied
0227         // partons
0228         if (pOutTemp[0].pstat() != droplet_stat &&
0229             pOutTemp[0].pstat() != miss_stat &&
0230             pOutTemp[0].pstat() != neg_stat &&
0231             !pOutTemp[0].isPhoton(pOutTemp[0].pid())) {
0232           vStartVecTemp.push_back(vStart);
0233         }
0234       } else {
0235         for (int k = 0; k < pOutTemp.size(); k++) {
0236           int edgeid = 0;
0237           if (pOutTemp[k].pstat() == neg_stat) {
0238             node vNewRootNode = pShower->new_vertex(
0239                 make_shared<Vertex>(0, 0, 0, currentTime - deltaT));
0240             edgeid = pShower->new_parton(vNewRootNode, vStart,
0241                                          make_shared<Parton>(pOutTemp[k]));
0242           } else {
0243             vEnd =
0244                 pShower->new_vertex(make_shared<Vertex>(0, 0, 0, currentTime));
0245             edgeid = pShower->new_parton(vStart, vEnd,
0246                                          make_shared<Parton>(pOutTemp[k]));
0247           }
0248           pOutTemp[k].set_shower(pShower);
0249           pOutTemp[k].set_edgeid(edgeid);
0250
0251           // no need to generate a vStart for photons and liquefied
0252           // partons
0253           if (pOutTemp[k].pstat() != droplet_stat &&
0254               pOutTemp[k].pstat() != miss_stat &&
0255               pOutTemp[k].pstat() != neg_stat &&
0256               !pOutTemp[k].isPhoton(pOutTemp[k].pid())) {
0257             vStartVecOut.push_back(vEnd);
0258           }
0259
0260           // --------------------------------------------
0261           // Add new roots from ElossModules ...
0262           // (maybe add for clarity a new vector in the signal!???)
0263           // Otherwise keep track of input size (so far always 1
0264           // and check if size > 1 and create additional root nodes to that vertex ...
0265           // Simple Test here below:
0266           // DEBUG:
0267           //cout<<"In JetEnergyloss : "<<pInTempModule.size()<<end;
0268           if (pInTempModule.size() > 1) {
0269             VERBOSE(7) << pInTempModule.size() - 1
0270                        << " new root node(s) to be added ...";
0271             //cout << pInTempModule.size()-1
0272             //     << " new root node(s) to be added ..." << endl;
0273
0274             for (int l = 1; l < pInTempModule.size(); l++) {
0275               node vNewRootNode = pShower->new_vertex(
0276                   make_shared<Vertex>(0, 0, 0, currentTime - deltaT));
0277               pShower->new_parton(vNewRootNode, vEnd,
0278                                   make_shared<Parton>(pInTempModule[l]));
0279             }
0280           }
0281         }
0282       }
0283
0284       // update parton shower
0285       if (pOutTemp.size() == 0) {
0286         // this is the free-streaming case for MATTER
0287         // do not push back droplets
0288         if (!weak_ptr_is_uninitialized(liquefier_ptr)) {
0289           if (pInTempModule[0].pstat() == droplet_stat)
0290             continue;
0291           if (pInTempModule[0].pstat() == miss_stat)
0292             continue;
0293           if (pInTempModule[0].pstat() == neg_stat)
0294             continue;
0295         }
0296         // do not push back photons
0297         if (pInTempModule[0].isPhoton(pInTempModule[0].pid()))
0298           continue;
0299         pInTemp.push_back(pInTempModule[0]);
0300       } else if (pOutTemp.size() == 1) {
0301         // this is the free-streaming case for MARTINI or LBT
0302         // do not push back droplets
0303         if (!weak_ptr_is_uninitialized(liquefier_ptr)) {
0304           if (pOutTemp[0].pstat() == droplet_stat)
0305             continue;
0306           if (pOutTemp[0].pstat() == miss_stat)
0307             continue;
0308           if (pOutTemp[0].pstat() == neg_stat)
0309             continue;
0310         }
0311         // do not push back photons
0312         if (pOutTemp[0].isPhoton(pOutTemp[0].pid()))
0313           continue;
0314         pInTemp.push_back(pOutTemp[0]);
0315       } else {
0316         for (int k = 0; k < pOutTemp.size(); k++) {
0317           // do not push back droplets
0318           if (!weak_ptr_is_uninitialized(liquefier_ptr)) {
0319             if (pOutTemp[k].pstat() == droplet_stat)
0320               continue;
0321             if (pOutTemp[k].pstat() == miss_stat)
0322               continue;
0323             if (pOutTemp[k].pstat() == neg_stat)
0324               continue;
0325           }
0326           // do not push back missing (from AdSCFT)
0327           if (pOutTemp[k].pstat() == miss_stat)
0328             continue;
0329           // do not push back photons
0330           if (pOutTemp[k].isPhoton(pOutTemp[k].pid()))
0331             continue;
0332
0333           pOut.push_back(pOutTemp[k]);
0334         }
0335       }
0336     }
0337
0338     // one time step is finished, now update parton shower to pIn
0339     pIn.clear();
0340     pIn.insert(pIn.end(), pInTemp.begin(), pInTemp.end());
0341     pIn.insert(pIn.end(), pOut.begin(), pOut.end());
0342
0343     // update vertex vector
0344     vStartVec.clear();
0345     vStartVec.insert(vStartVec.end(), vStartVecTemp.begin(),
0346                      vStartVecTemp.end());
0347     vStartVec.insert(vStartVec.end(), vStartVecOut.begin(), vStartVecOut.end());
0348   } while (currentTime < maxT); // other criteria (how to include; TBD)
0349
0350   pIn.clear();
0351   vStartVec.clear();
0352 }
0353
0354 void JetEnergyLoss::Exec() {
0355   VERBOSE(1) << "Run JetEnergyLoss ...";
0356   VERBOSE(1) << "Found " << GetNumberOfTasks()
0357              << " Eloss Tasks/Modules Execute them ... ";
0358   //DEBUGTHREAD<<"Task Id = "<<this_thread::get_id()<<" | Run JetEnergyLoss ...";
0359   //DEBUGTHREAD<<"Task Id = "<<this_thread::get_id()<<" | Found "<<GetNumberOfTasks()<<" Eloss Tasks/Modules Execute them ... ";
0360
0361   if (GetShowerInitiatingParton()) {
0362     pShower = make_shared<PartonShower>();
0363
0364     /*
0365        //Check Memory ...
0366        VERBOSE(8)<<"Use PartonShowerGenerator to do Parton shower stored in PartonShower Graph class";
0367        JSDEBUG<<"Use PartonShowerGenerator to do Parton shower stored in PartonShower Graph class";
0368
0369        PartonShowerGenerator PSG;
0370        PSG.DoShower(*shared_from_this()); //needed otherwise all signal slots have to be recreated for shower module ....
0371        // Overall not the nicest logic though .... Just to make changing and expanding the shower code in the future ...
0372        // (basically, just now to remove the code out of the jet energy loss class ...) TBD
0373        // also not really nice, since now the energy loss part in the parton shower and not really visible in this class ...
0374        // Keep both codes so far ...
0375        */
0376
0377     // Shower handled in this class ...
0378     DoShower();
0379
0380     pShower->PrintNodes();
0381     pShower->PrintEdges();
0382
0383     weak_ptr<HardProcess> hproc =
0384         JetScapeSignalManager::Instance()->GetHardProcessPointer();
0385
0386     for (unsigned int ipart = 0; ipart < pShower->GetNumberOfPartons();
0387          ipart++) {
0388       //   Uncomment to dump the whole parton shower into the parton container
0389       // auto hp = hproc.lock();
0390       // if ( hp ) hp->AddParton(pShower->GetPartonAt(ipart));
0391     }
0392
0393     shared_ptr<PartonPrinter> pPrinter =
0394         JetScapeSignalManager::Instance()->GetPartonPrinterPointer().lock();
0395     if (pPrinter) {
0396       pPrinter->GetFinalPartons(pShower);
0397     }
0398
0399     shared_ptr<JetEnergyLoss> pEloss =
0400         JetScapeSignalManager::Instance()->GetEnergyLossPointer().lock();
0401     if (pEloss) {
0402       pEloss->GetFinalPartonsForEachShower(pShower);
0403     }
0404   } else {
0405     JSWARN << "NO Initial Hard Parton for Parton shower received ...";
0406   }
0407
0408   //DEBUGTHREAD<<"Task Id = "<<this_thread::get_id()<<" Finished!";
0409   //JetScapeTask::ExecuteTasks(); // prevent Further modules to be execute, everything done by JetEnergyLoss ... (also set the no active flag ...!?)
0410 }
0411
0412 void JetEnergyLoss::WriteTask(weak_ptr<JetScapeWriter> w) {
0413   VERBOSE(8);
0414   VERBOSE(4) << "In JetEnergyLoss::WriteTask";
0415   auto f = w.lock();
0416   if (!f)
0417     return;
0418
0419   f->WriteComment("Energy loss Shower Initating Parton: " + GetId());
0420   f->Write(inP);
0421
0422   VERBOSE(4) << " writing partons... found " << pShower->GetNumberOfPartons();
0423   f->Write(pShower);
0424 }
0425
0426 void JetEnergyLoss::PrintShowerInitiatingParton() {
0427   //JSDEBUG<<inP->pid();
0428 }
0429
0430 void JetEnergyLoss::GetFinalPartonsForEachShower(
0431     shared_ptr<PartonShower> shower) {
0432
0433   this->final_Partons.push_back(shower.get()->GetFinalPartons());
0434 }
0435
0436 } // end namespace Jetscape