src/hadronization/ColoredHadronization.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 "ColoredHadronization.h"
0017 #include "JetScapeXML.h"
0018 #include "JetScapeLogger.h"
0019 #include "tinyxml2.h"
0020
0021 using namespace Jetscape;
0022 using namespace Pythia8;
0023
0024 // Register the module with the base class
0025 RegisterJetScapeModule<ColoredHadronization>
0026     ColoredHadronization::reg("ColoredHadronization");
0027
0028 Pythia8::Pythia ColoredHadronization::pythia("IntentionallyEmpty", false);
0029
0030 ColoredHadronization::ColoredHadronization() {
0031   SetId("MyHadroTest");
0032   VERBOSE(8);
0033 }
0034
0035 ColoredHadronization::~ColoredHadronization() { VERBOSE(8); }
0036
0037 void ColoredHadronization::Init() {
0038
0039   std::string s = GetXMLElementText({"JetHadronization", "name"});
0040   JSDEBUG << s << " to be initializied ...";
0041
0042   double p_read_xml =
0043       GetXMLElementDouble({"JetHadronization", "eCMforHadronization"});
0044   p_fake = p_read_xml;
0045
0046   /*std::string weak_decays =
0047       GetXMLElementText({"JetHadronization", "weak_decays"});*/
0048
0049   VERBOSE(2) << "Start Hadronizing using the PYTHIA module...";
0050
0051   // Show initialization at DEBUG or high verbose level
0052   pythia.readString("Init:showProcesses = off");
0053   pythia.readString("Init:showChangedSettings = off");
0054   pythia.readString("Init:showMultipartonInteractions = off");
0055   pythia.readString("Init:showChangedParticleData = off");
0056   if (JetScapeLogger::Instance()->GetDebug() ||
0057       JetScapeLogger::Instance()->GetVerboseLevel() > 2) {
0058     pythia.readString("Init:showProcesses = on");
0059     pythia.readString("Init:showChangedSettings = on");
0060     pythia.readString("Init:showMultipartonInteractions = on");
0061     pythia.readString("Init:showChangedParticleData = on");
0062   }
0063
0064   // No event record printout.
0065   pythia.readString("Next:numberShowInfo = 0");
0066   pythia.readString("Next:numberShowProcess = 0");
0067   pythia.readString("Next:numberShowEvent = 0");
0068   if (JetScapeLogger::Instance()->GetDebug() ||
0069       JetScapeLogger::Instance()->GetVerboseLevel() > 2) {
0070     pythia.readString("Next:numberShowInfo = 1");
0071     pythia.readString("Next:numberShowProcess = 1");
0072     pythia.readString("Next:numberShowEvent = 1");
0073   }
0074
0075   pythia.readString("ProcessLevel:all = off");
0076   pythia.readString("PartonLevel:FSR=off");
0077
0078   // General settings for hadron decays
0079   std::string pythia_decays = GetXMLElementText({"JetHadronization", "pythia_decays"});
0080   double tau0Max = 10.0;
0081   double tau0Max_xml = GetXMLElementDouble({"JetHadronization", "tau0Max"});
0082     if(tau0Max_xml >= 0){tau0Max = tau0Max_xml;}
0083   else{JSWARN << "tau0Max should be larger than 0. Set it to 10.";}
0084   if(pythia_decays == "on"){
0085     JSINFO << "Pythia decays are turned on for tau0Max < " << tau0Max;
0086     pythia.readString("HadronLevel:Decay = on");
0087     pythia.readString("ParticleDecays:limitTau0 = on");
0088     pythia.readString("ParticleDecays:tau0Max = " + std::to_string(tau0Max));
0089   } else {
0090     JSINFO << "Pythia decays are turned off";
0091     pythia.readString("HadronLevel:Decay = off");
0092   }
0093
0094   // Settings for decays (old flag, will be depracted at some point)
0095   // This overwrites the previous settings if the user xml file contains the flag
0096   std::string weak_decays =
0097     GetXMLElementText({"JetHadronization", "weak_decays"});
0098   if (weak_decays == "off") {
0099     JSINFO << "Hadron decays are turned off.";
0100     JSWARN << "This parameter will be depracted at some point. Use 'pythia_decays' instead.\nOverwriting 'pythia_decays'.";
0101     pythia.readString("HadronLevel:Decay = off");
0102   } else if(weak_decays == "on") {
0103     JSINFO << "Hadron decays inside a range of 10 mm/c are turned on.";
0104     JSWARN << "This parameter will be depracted at some point. Use 'pythia_decays' and 'tau0Max' for more control on decays.\nOverwriting 'pythia_decays' and fix 'tau0Max' to 10.";
0105     pythia.readString("HadronLevel:Decay = on");
0106     pythia.readString("ParticleDecays:limitTau0 = on");
0107     pythia.readString("ParticleDecays:tau0Max = 10.0");
0108   }
0109
0110   std::stringstream lines;
0111   lines << GetXMLElementText({"JetHadronization", "LinesToRead"}, false);
0112   while (std::getline(lines, s, '\n')) {
0113     if (s.find_first_not_of(" \t\v\f\r") == s.npos)
0114       continue; // skip empty lines
0115     JSINFO << "Also reading in: " << s;
0116     pythia.readString(s);
0117   }
0118
0119   pythia.init();
0120 }
0121
0122 void ColoredHadronization::WriteTask(weak_ptr<JetScapeWriter> w) {
0123   VERBOSE(8);
0124   auto f = w.lock();
0125   if (!f)
0126     return;
0127   f->WriteComment("Hadronization Module : " + GetId());
0128   f->WriteComment("Hadronization to be implemented accordingly ...");
0129 }
0130
0131 void ColoredHadronization::DoHadronization(
0132     vector<vector<shared_ptr<Parton>>> &shower,
0133     vector<shared_ptr<Hadron>> &hOut, vector<shared_ptr<Parton>> &pOut) {
0134
0135   Event &event = pythia.event;
0136   event.reset();
0137   double pz = p_fake;
0138
0139   JSDEBUG << "&&&&&&&&&&&&&&&&&&& the number of showers are: " << shower.size();
0140   for (unsigned int ishower = 0; ishower < shower.size(); ++ishower) {
0141     JSDEBUG << "&&&&&&&&&&&&&&&&&&& there are " << shower.at(ishower).size()
0142             << " partons in the shower number " << ishower;
0143     for (unsigned int ipart = 0; ipart < shower.at(ishower).size(); ++ipart) {
0144       double onshellE = pow(pow(shower.at(ishower).at(ipart)->px(), 2) +
0145                                 pow(shower.at(ishower).at(ipart)->py(), 2) +
0146                                 pow(shower.at(ishower).at(ipart)->pz(), 2),
0147                             0.5);
0148
0149       if (shower.at(ishower).at(ipart)->pid() == 22) {
0150
0151         VERBOSE(1) << BOLDYELLOW
0152                    << " photon found in colored hadronization with ";
0153         VERBOSE(1) << BOLDYELLOW
0154                    << "px = " << shower.at(ishower).at(ipart)->px();
0155         //cin >> blurb;
0156       }
0157       event.append(shower.at(ishower).at(ipart)->pid(), 23,
0158                    shower.at(ishower).at(ipart)->color(),
0159                    shower.at(ishower).at(ipart)->anti_color(),
0160                    shower.at(ishower).at(ipart)->px(),
0161                    shower.at(ishower).at(ipart)->py(),
0162                    shower.at(ishower).at(ipart)->pz(), onshellE);
0163     }
0164
0165     //first, find unpaired color and anticolor tags.
0166     std::vector<int> cols;
0167     std::vector<int> acols;
0168     for (unsigned int ipart = 0; ipart < shower.at(ishower).size(); ++ipart) {
0169       if (shower.at(ishower).at(ipart)->pid() == 22) {
0170         continue;
0171       }
0172       if (shower.at(ishower).at(ipart)->color() != 0) {
0173         cols.push_back(shower.at(ishower).at(ipart)->color());
0174       }
0175       if (shower.at(ishower).at(ipart)->anti_color() != 0) {
0176         acols.push_back(shower.at(ishower).at(ipart)->anti_color());
0177       }
0178     }
0179     //the outcomes are: 1-unpaired color tag, 2-unpaired anticolor tag, 3-both an unpaired color & anticolor tag, 4-no unpaired tags
0180     //1-add an antiquark, 2-add a quark, 3-add a gluon, 4-add nothing (possibly photon only event)
0181     int icol = 0;
0182     while (icol < cols.size()) {
0183       bool foundpair = false;
0184       for (int iacol = 0; iacol < acols.size(); ++iacol) {
0185         if (cols[icol] == acols[iacol]) {
0186           cols.erase(cols.begin() + icol);
0187           acols.erase(acols.begin() + iacol);
0188           foundpair = true;
0189           continue;
0190         }
0191       }
0192       if (!foundpair) {
0193         ++icol;
0194       }
0195     }
0196
0197     int pid = 0;
0198     int color = 0;
0199     int anti_color = 0;
0200     if ((cols.size() > 0) && (acols.size() > 0)) {
0201       pid = 21;
0202       color = cols[0];
0203       anti_color = acols[0];
0204     } else if ((cols.size() > 0) && (acols.size() == 0)) {
0205       pid = -1;
0206       color = cols[0];
0207       anti_color = 0;
0208     } else if ((cols.size() == 0) && (acols.size() > 0)) {
0209       pid = 1;
0210       color = 0;
0211       anti_color = acols[0];
0212     }
0213
0214     if (pid != 0) {
0215       pz = -1 * pz;
0216       event.append(pid, 23, anti_color, color, 0.2, 0.2, pz,
0217                    sqrt(pz * pz + 0.08));
0218     }
0219
0220     VERBOSE(2) << "There are " << hOut.size() << " Hadrons and " << pOut.size()
0221                << " partons after Hadronization";
0222   }
0223
0224   //there still may be color tag duplicates - will SegFault if color_reconnections is ever invoked.
0225   //this should be fixed *here*, before pythia.next() below, if that's ever a concern.
0226   //scan over list of color & anticolor tags - if we find a duplicate, set it to a new value >0, <101 (will fail for more than 100 duplicates)
0227   std::vector<std::vector<int>> col_instances;
0228   for (unsigned int i = 0; i < event.size(); ++i) {
0229     bool newcol = true;
0230     bool newacol = true;
0231     for (int icols = 0; icols < col_instances.size(); ++icols) {
0232       if (event[i].col() == col_instances[icols][0]) {
0233         ++col_instances[icols][1];
0234         newcol = false;
0235       }
0236       if (event[i].acol() == col_instances[icols][0]) {
0237         ++col_instances[icols][1];
0238         newacol = false;
0239       }
0240     }
0241     if (newcol && (event[i].col() != 0)) {
0242       std::vector<int> tmpcol;
0243       tmpcol.push_back(event[i].col());
0244       tmpcol.push_back(1);
0245       col_instances.push_back(tmpcol);
0246     }
0247     if (newacol && (event[i].acol() != 0)) {
0248       std::vector<int> tmpcol;
0249       tmpcol.push_back(event[i].acol());
0250       tmpcol.push_back(1);
0251       col_instances.push_back(tmpcol);
0252     }
0253   }
0254   col_instances.erase(
0255       std::remove_if(col_instances.begin(), col_instances.end(),
0256                      [](const std::vector<int> &val) { return val[1] <= 2; }),
0257       col_instances.end());
0258   int updcol = 1;
0259   while (col_instances.size() > 0) {
0260     int nupd = 2;
0261     for ( int i = event.size() - 1; i >= 0; --i) {
0262       if (col_instances[0][0] == event[i].col()) {
0263         event[i].col(updcol);
0264         --nupd;
0265       }
0266       if (col_instances[0][0] == event[i].acol()) {
0267         event[i].acol(updcol);
0268         --nupd;
0269       }
0270       if (nupd == 0) {
0271         break;
0272       }
0273     }
0274     ++updcol;
0275     col_instances[0][1] -= 2;
0276     col_instances.erase(
0277         std::remove_if(col_instances.begin(), col_instances.end(),
0278                        [](const std::vector<int> &val) { return val[1] <= 2; }),
0279         col_instances.end());
0280   }
0281
0282   pythia.next();
0283   // event.list();
0284
0285   unsigned int ip = hOut.size();
0286   for (unsigned int i = 0; i < event.size(); ++i) {
0287     if (!event[i].isFinal())
0288       continue;
0289     //if ( !event[i].isHadron() )  continue;
0290     if (fabs(event[i].eta()) > 20)
0291       continue; //To prevent "nan" from propagating, very rare though
0292
0293     double x[4] = {0, 0, 0, 0};
0294     hOut.push_back(make_shared<Hadron>(ip, event[i].id(), event[i].status(),
0295                                        event[i].pT(), event[i].eta(),
0296                                        event[i].phi(), event[i].e(), x));
0297     ++ip;
0298   }
0299
0300   shower.clear();
0301 }