Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-03 08:19:58 

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 "ColorlessHadronization.h"
0017 #include "JetScapeXML.h"
0018 #include "JetScapeLogger.h"
0019 #include "tinyxml2.h"
0020 #include "JetScapeConstants.h"
0021 #include <sstream>
0022 #include <iostream>
0023 #include <fstream>
0024 #include <sstream>
0025 #include <random>
0026 
0027 using namespace Jetscape;
0028 using namespace Pythia8;
0029 
0030 // Hadrons output file
0031 //ofstream hadfile;
0032 
0033 // Register the module with the base class
0034 RegisterJetScapeModule<ColorlessHadronization>
0035     ColorlessHadronization::reg("ColorlessHadronization");
0036 
0037 // Initialize static helper here
0038 Pythia8::Pythia ColorlessHadronization::pythia("IntentionallyEmpty", false);
0039 
0040 ColorlessHadronization::ColorlessHadronization() {
0041   SetId("ColorlessHadronization");
0042   VERBOSE(8);
0043 }
0044 
0045 ColorlessHadronization::~ColorlessHadronization() { VERBOSE(8); }
0046 
0047 void ColorlessHadronization::Init() {
0048   // Open output file
0049   //hadfile.open("CH_myhad.dat");
0050 
0051   std::string s = GetXMLElementText({"JetHadronization", "name"});
0052   JSDEBUG << s << " to be initializied ...";
0053 
0054   // Read sqrts to know remnants energies
0055   double p_read_xml =
0056       GetXMLElementDouble({"JetHadronization", "eCMforHadronization"});
0057   p_fake = p_read_xml;
0058 
0059   take_recoil = GetXMLElementInt({"JetHadronization", "take_recoil"});
0060 
0061   JSDEBUG << "Initialize ColorlessHadronization";
0062   VERBOSE(8);
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 
0069   // Standard settings
0070   pythia.readString("ProcessLevel:all = off");
0071   pythia.readString("PartonLevel:FSR=off");
0072 
0073   // General settings for hadron decays
0074   std::string pythia_decays = GetXMLElementText({"JetHadronization", "pythia_decays"});
0075   double tau0Max = 10.0;
0076   double tau0Max_xml = GetXMLElementDouble({"JetHadronization", "tau0Max"});
0077     if(tau0Max_xml >= 0){tau0Max = tau0Max_xml;}
0078   else{JSWARN << "tau0Max should be larger than 0. Set it to 10.";}
0079   if(pythia_decays == "on"){
0080     JSINFO << "Pythia decays are turned on for tau0Max < " << tau0Max;
0081     pythia.readString("HadronLevel:Decay = on");
0082     pythia.readString("ParticleDecays:limitTau0 = on");
0083     pythia.readString("ParticleDecays:tau0Max = " + std::to_string(tau0Max));
0084   } else {
0085     JSINFO << "Pythia decays are turned off";
0086     pythia.readString("HadronLevel:Decay = off");
0087   }
0088 
0089   // Settings for decays (old flag, will be depracted at some point)
0090   // This overwrites the previous settings if the user xml file contains the flag
0091   std::string weak_decays =
0092     GetXMLElementText({"JetHadronization", "weak_decays"});
0093   if (weak_decays == "off") {
0094     JSINFO << "Hadron decays are turned off.";
0095     JSWARN << "This parameter will be depracted at some point. Use 'pythia_decays' instead.\nOverwriting 'pythia_decays'.";
0096     pythia.readString("HadronLevel:Decay = off");
0097   } else if(weak_decays == "on") {
0098     JSINFO << "Hadron decays inside a range of 10 mm/c are turned on.";
0099     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.";
0100     pythia.readString("HadronLevel:Decay = on");
0101     pythia.readString("ParticleDecays:limitTau0 = on");
0102     pythia.readString("ParticleDecays:tau0Max = 10.0");
0103   }
0104 
0105   std::stringstream lines;
0106   lines << GetXMLElementText({"JetHadronization", "LinesToRead"}, false);
0107   while (std::getline(lines, s, '\n')) {
0108     if (s.find_first_not_of(" \t\v\f\r") == s.npos)
0109       continue; // skip empty lines
0110     JSINFO << "Also reading in: " << s;
0111     pythia.readString(s);
0112   }
0113 
0114   // Initialize random number distribution
0115   ZeroOneDistribution = std::uniform_real_distribution<double> { 0.0, 1.0 };
0116   // And initialize
0117   pythia.init();
0118 }
0119 
0120 void ColorlessHadronization::WriteTask(weak_ptr<JetScapeWriter> w) {
0121   VERBOSE(8);
0122   auto f = w.lock();
0123   if (!f)
0124     return;
0125   f->WriteComment("Hadronization Module : " + GetId());
0126 }
0127 
0128 void ColorlessHadronization::DoHadronization(
0129     vector<vector<shared_ptr<Parton>>> &shower,
0130     vector<shared_ptr<Hadron>> &hOut, vector<shared_ptr<Parton>> &pOut) {
0131   VERBOSE(1) << "Start Hadronizing using PYTHIA Lund string model (does NOT "
0132                 "use color flow, needs to be tested)...";
0133   Event &event = pythia.event;
0134   ParticleData &pdt = pythia.particleData;
0135 
0136   // Hadronize positive (status = 0) and negative (status = -1) partons in a different space
0137   for (int want_pos = 1; want_pos >= 0; --want_pos) {
0138     event.reset();
0139 
0140     if (!take_recoil && want_pos == 0)
0141       continue; // SC: don't need negative if don't take recoil
0142 
0143     // Set remnants momentum
0144     double random_direction = ZeroOneDistribution(*GetMt19937Generator());
0145     if (random_direction>0.5) {
0146       random_direction = 1.0;
0147     }else{
0148       random_direction = -1.0;
0149     };
0150 
0151     double rempx = 0.2;
0152     double rempy = 0.2;
0153     double rempz = random_direction*p_fake;
0154     double reme = std::sqrt(std::pow(rempx, 2.) + std::pow(rempy, 2.) +
0155                             std::pow(rempz, 2.));
0156 
0157     // Hadronize all showers together
0158     vector<shared_ptr<Parton>> pIn;
0159     vector<vector<Parton>> shower_in;
0160     for (unsigned int ishower = 0; ishower < shower.size(); ++ishower) {
0161       vector<Parton> p_sh;
0162       for (unsigned int ipart = 0; ipart < shower.at(ishower).size(); ++ipart) {
0163         p_sh.push_back(*(shower[ishower][ipart]));
0164       }
0165       shower_in.push_back(p_sh);
0166     }
0167     for (unsigned int ishower = 0; ishower < shower_in.size(); ++ishower) {
0168       for (unsigned int ipart = 0; ipart < shower_in[ishower].size(); ++ipart) {
0169         //if (shower_in.at(ishower).at(ipart)->pstat()==0 && want_pos==1) pIn.push_back(shower_in.at(ishower).at(ipart));  // Positive
0170         if (want_pos == 1) { // Positive
0171           if (take_recoil && shower_in[ishower][ipart].pstat() == 1) {
0172             pIn.push_back(make_shared<Parton>(shower_in[ishower][ipart]));
0173           }
0174           if (shower_in[ishower][ipart].pstat() == 0) {
0175             pIn.push_back(make_shared<Parton>(shower_in[ishower][ipart]));
0176           }
0177           if (shower_in[ishower][ipart].pstat() == 22) {
0178             pIn.push_back(make_shared<Parton>(shower_in[ishower][ipart])); //Allow photons with status code 22 to pass through Colorless hadronization.
0179           }
0180         }
0181         if (take_recoil && shower_in[ishower][ipart].pstat() == -1 &&
0182             want_pos == 0) {
0183           pIn.push_back(make_shared<Parton>(shower_in[ishower][ipart]));
0184         } // Negative
0185       }
0186       JSDEBUG << "Shower#" << ishower + 1
0187               << ". Number of partons to hadronize so far: " << pIn.size();
0188     }
0189     if (want_pos == 1)
0190       VERBOSE(1) << "# Positive Partons to hadronize: " << pIn.size();
0191     else
0192       VERBOSE(1) << "# Negative Partons to hadronize: " << pIn.size();
0193 
0194     // Check whether event is empty (specially important for negative partons case)
0195     if (pIn.size() == 0)
0196       continue;
0197 
0198     int col[pIn.size() + 2], acol[pIn.size() + 2], isdone[pIn.size() + 2];
0199     memset(col, 0, (pIn.size() + 2) * sizeof(int)),
0200         memset(acol, 0, (pIn.size() + 2) * sizeof(int)),
0201         memset(isdone, 0, (pIn.size() + 2) * sizeof(int));
0202 
0203     // Find number of quarks
0204     int nquarks = 0;
0205     int isquark[pIn.size() + 2];
0206     memset(isquark, 0, (pIn.size() + 2) * sizeof(int));
0207     for (unsigned int ipart = 0; ipart < pIn.size(); ++ipart) {
0208       if (abs(pIn[ipart]->pid()) <= 6) {
0209         isquark[nquarks] = ipart;
0210         nquarks += 1;
0211       }
0212     }
0213     JSDEBUG << "#Quarks = " << nquarks;
0214 
0215     // Find number of strings
0216     int nstrings = max(int(double(nquarks) / 2. + 0.6), 1);
0217     JSDEBUG << "#Strings = " << nstrings;
0218 
0219     // If there are no quarks, need to attach two of them
0220     int istring = 0;
0221     int one_end[nstrings], two_end[nstrings];
0222     if (nquarks == 0) { // Only attach remnants if event is not empty
0223       // First quark
0224       FourVector p1(rempx, rempy, rempz, reme);
0225       FourVector x1;
0226       pIn.push_back(std::make_shared<Parton>(0, 1, 0, p1, x1));
0227       isquark[nquarks] = pIn.size() - 1;
0228       nquarks += 1;
0229       isdone[pIn.size() - 1] = 1;
0230       one_end[0] = pIn.size() - 1;
0231       VERBOSE(1) << "Attached quark remnant flying down +Pz beam";
0232       // Second quark
0233       FourVector p2(rempx, rempy, -rempz, reme);
0234       FourVector x2;
0235       pIn.push_back(std::make_shared<Parton>(0, 1, 0, p2, x2));
0236       isquark[nquarks] = pIn.size() - 1;
0237       nquarks += 1;
0238       isdone[pIn.size() - 1] = 1;
0239       two_end[istring] = pIn.size() - 1;
0240       VERBOSE(1) << "Attached quark remnant flying down -Pz beam";
0241     }
0242 
0243     // Assign ends of strings (order matters in this algo)
0244     for (unsigned int iquark = 0; iquark < nquarks; iquark++) {
0245       if (isdone[isquark[iquark]] == 0) {
0246         isdone[isquark[iquark]] = 1;
0247         one_end[istring] = isquark[iquark];
0248         double min_delR = 10000.;
0249         int partner = -2;
0250         for (unsigned int jquark = 0; jquark < nquarks; jquark++) {
0251           if (iquark == jquark)
0252             continue;
0253           int d_jquark = isquark[jquark];
0254           if (isdone[d_jquark] == 0) {
0255             fjcore::PseudoJet pf(pIn[d_jquark]->px(), pIn[d_jquark]->py(),
0256                                  pIn[d_jquark]->pz(), pIn[d_jquark]->e());
0257             double delR = pIn[isquark[iquark]]->delta_R(pf);
0258             if (delR < min_delR)
0259               min_delR = delR, partner = jquark;
0260           }
0261         }
0262         if (partner != -2) {
0263           isdone[isquark[partner]] = 1;
0264           two_end[istring] = isquark[partner];
0265           istring += 1;
0266         } else {
0267           FourVector p(rempx, rempy, rempz, reme);
0268           FourVector x;
0269           pIn.push_back(std::make_shared<Parton>(0, 1, 0, p, x));
0270           isquark[nquarks] = pIn.size() - 1;
0271           nquarks += 1;
0272           isdone[pIn.size() - 1] = 1;
0273           two_end[istring] = pIn.size() - 1;
0274           VERBOSE(1) << "Attached quark remnant flying down +Pz beam";
0275         }
0276       }
0277     }
0278 
0279     // Assign gluons to a certain string
0280     int my_string[pIn.size()];
0281     memset(my_string, 0, pIn.size() * sizeof(int));
0282     for (unsigned int ipart = 0; ipart < pIn.size(); ++ipart) {
0283       if (pIn[ipart]->pid() == 21) {
0284         double min_delR = 100000.;
0285         for (int ns = 0; ns < nstrings; ns++) {
0286           int fq = one_end[ns];
0287           int sq = two_end[ns];
0288           fjcore::PseudoJet pfq(pIn[fq]->px(), pIn[fq]->py(), pIn[fq]->pz(),
0289                                 pIn[fq]->e());
0290           double f_delR = pIn[ipart]->delta_R(pfq);
0291           fjcore::PseudoJet psq(pIn[sq]->px(), pIn[sq]->py(), pIn[sq]->pz(),
0292                                 pIn[sq]->e());
0293           double s_delR = pIn[ipart]->delta_R(psq);
0294           double delR = (f_delR + s_delR) / 2.;
0295           if (delR < min_delR)
0296             my_string[ipart] = ns, min_delR = delR;
0297         }
0298       }
0299     }
0300 
0301     // Build up chain using gluons assigned to each string, in a closest pair order
0302     int lab_col = 102;
0303     for (int ns = 0; ns < nstrings; ns++) {
0304       int tquark = one_end[ns];
0305       if (pIn[tquark]->pid() > 0)
0306         col[tquark] = lab_col;
0307       else
0308         acol[tquark] = lab_col;
0309       lab_col += 1;
0310       int link = tquark;
0311       int changes = 1;
0312       do {
0313         changes = 0;
0314         double min_delR = 100000.;
0315         int next_link = 0;
0316         for (unsigned int ipart = 0; ipart < pIn.size(); ++ipart) {
0317           if (pIn[ipart]->pid() == 21 && isdone[ipart] == 0 &&
0318               my_string[ipart] == ns) {
0319             changes = 1;
0320             fjcore::PseudoJet pf(pIn[ipart]->px(), pIn[ipart]->py(),
0321                                  pIn[ipart]->pz(), pIn[ipart]->e());
0322             double delR = pIn[link]->delta_R(pf);
0323             if (delR < min_delR)
0324               min_delR = delR, next_link = ipart;
0325           }
0326         }
0327         if (changes == 1) {
0328           isdone[next_link] = 1;
0329           if (col[link] == lab_col - 1)
0330             col[next_link] = lab_col, acol[next_link] = lab_col - 1;
0331           else
0332             col[next_link] = lab_col - 1, acol[next_link] = lab_col;
0333           lab_col += 1;
0334           JSDEBUG << " Linked parton= " << next_link;
0335           link = next_link;
0336         }
0337       } while (changes == 1);
0338       // Attach second end
0339       if (col[link] == lab_col - 1)
0340         col[two_end[ns]] = 0, acol[two_end[ns]] = lab_col - 1;
0341       else
0342         col[two_end[ns]] = lab_col - 1, acol[two_end[ns]] = 0;
0343     }
0344     // Changing identity of quarks to be consistent with color charge
0345     for (int iq = 0; iq < nquarks; ++iq) {
0346       if (col[isquark[iq]] != 0) {
0347         if (pIn[isquark[iq]]->pid() < 0)
0348           pIn[isquark[iq]]->set_id(-pIn[isquark[iq]]->pid());
0349       } else {
0350         if (pIn[isquark[iq]]->pid() > 0)
0351           pIn[isquark[iq]]->set_id(-pIn[isquark[iq]]->pid());
0352       }
0353     }
0354 
0355     // Introduce partons into PYTHIA
0356     /*
0357     for (unsigned int ipart=0; ipart <  pIn.size(); ++ipart)
0358     {
0359       JSDEBUG << "Parton #" << ipart << " is a " << pIn[ipart]->pid() << "with energy = " << pIn[ipart]->e() << " with phi= " << pIn[ipart]->phi() << " and has col= " << col[ipart] << " and acol= " << acol[ipart];
0360     }
0361     */
0362 
0363     /***************************************************************************************************************/
0364     //
0365     // Making collinear partons not collinear
0366     //
0367     // Could have dangerous effects, yet to be tested....
0368     //
0369 
0370     for (unsigned int ipart = 0; ipart < pIn.size(); ++ipart) {
0371       double px = pIn[ipart]->px();
0372       double py = pIn[ipart]->py();
0373       double pz = pIn[ipart]->pz();
0374       double ee = pIn[ipart]->e();
0375       for (unsigned int j = ipart + 1; j < pIn.size(); j++) {
0376         double p2x = pIn[j]->px();
0377         double p2y = pIn[j]->py();
0378         double p2z = pIn[j]->pz();
0379         double e2e = pIn[j]->e();
0380 
0381         double diff =
0382             sqrt(pow(px - p2x, 2) + pow(py - p2y, 2) + pow(pz - p2z, 2));
0383         double f = 4.0;
0384         if (diff < f * Lambda_QCD) {
0385           if ((pz >= 0) && (p2z >= 0)) {
0386             if (pz >= p2z) {
0387               pIn[ipart]->reset_momentum(px, py, pz + f * Lambda_QCD,
0388                                          sqrt(ee * ee +
0389                                               2 * pz * f * Lambda_QCD +
0390                                               f * Lambda_QCD * f * Lambda_QCD));
0391             } else
0392               pIn[j]->reset_momentum(p2x, p2y, p2z + f * Lambda_QCD,
0393                                      sqrt(e2e * e2e + 2 * p2z * f * Lambda_QCD +
0394                                           f * Lambda_QCD * f * Lambda_QCD));
0395           } else if ((pz >= 0) && (p2z < 0)) {
0396             if (abs(pz) >= abs(p2z)) {
0397               pIn[ipart]->reset_momentum(px, py, pz + f * Lambda_QCD,
0398                                          sqrt(ee * ee +
0399                                               2 * pz * f * Lambda_QCD +
0400                                               f * Lambda_QCD * f * Lambda_QCD));
0401             } else
0402               pIn[j]->reset_momentum(p2x, p2y, p2z - f * Lambda_QCD,
0403                                      sqrt(e2e * e2e - 2 * p2z * f * Lambda_QCD +
0404                                           f * Lambda_QCD * f * Lambda_QCD));
0405           } else if ((pz < 0) && (p2z >= 0)) {
0406             if (abs(pz) >= abs(p2z)) {
0407               pIn[ipart]->reset_momentum(px, py, pz - f * Lambda_QCD,
0408                                          sqrt(ee * ee -
0409                                               2 * pz * f * Lambda_QCD +
0410                                               f * Lambda_QCD * f * Lambda_QCD));
0411             } else
0412               pIn[j]->reset_momentum(p2x, p2y, p2z + f * Lambda_QCD,
0413                                      sqrt(e2e * e2e + 2 * p2z * f * Lambda_QCD +
0414                                           f * Lambda_QCD * f * Lambda_QCD));
0415           } else {
0416             if (abs(pz) >= abs(p2z)) {
0417               pIn[ipart]->reset_momentum(px, py, pz - f * Lambda_QCD,
0418                                          sqrt(ee * ee -
0419                                               2 * pz * f * Lambda_QCD +
0420                                               f * Lambda_QCD * f * Lambda_QCD));
0421             } else
0422               pIn[j]->reset_momentum(p2x, p2y, p2z - f * Lambda_QCD,
0423                                      sqrt(e2e * e2e - 2 * p2z * f * Lambda_QCD +
0424                                           f * Lambda_QCD * f * Lambda_QCD));
0425           }
0426         }
0427       }
0428     }
0429 
0430     /**************************************************************************************************************************/
0431 
0432     for (unsigned int ipart = 0; ipart < pIn.size(); ++ipart) {
0433       int ide = pIn[ipart]->pid();
0434       double px = pIn[ipart]->px();
0435       double py = pIn[ipart]->py();
0436       double pz = pIn[ipart]->pz();
0437       double ee = pIn[ipart]->e();
0438       double mm = pdt.m0(int(ide));
0439       ee = std::sqrt(px * px + py * py + pz * pz + mm * mm);
0440       if (col[ipart] == 0 && acol[ipart] == 0 && (ide == 21 || abs(ide) <= 6)) {
0441         JSINFO << "Stopping because of colorless parton trying to be "
0442                   "introduced in PYTHIA string";
0443         exit(0);
0444       }
0445       event.append(int(ide), 23, col[ipart], acol[ipart], px, py, pz, ee, mm);
0446     }
0447 
0448     pythia.next();
0449     for (unsigned int ipart = 0; ipart < event.size(); ++ipart) {
0450       if (event[ipart].isFinal()) {
0451         int ide = pythia.event[ipart].id();
0452         FourVector p(pythia.event[ipart].px(), pythia.event[ipart].py(),
0453                      pythia.event[ipart].pz(), pythia.event[ipart].e());
0454         FourVector x;
0455         if (want_pos == 1)
0456           hOut.push_back(
0457               std::make_shared<Hadron>(Hadron(0, ide, 0, p, x))); // Positive
0458         else
0459           hOut.push_back(
0460               std::make_shared<Hadron>(Hadron(0, ide, -1, p, x))); // Negative
0461         //JSINFO << "Produced Hadron has id = " << pythia.event[ipart].id();
0462         // Print on output file
0463         //hadfile << pythia.event[ipart].px() << " " << pythia.event[ipart].py() << " " << pythia.event[ipart].pz() << " " << pythia.event[ipart].e() << " " << pythia.event[ipart].id() << " " << pythia.event[ipart].charge() << endl;
0464       }
0465     }
0466     VERBOSE(1) << "#Showers hadronized together: " << shower.size()
0467                << ". There are " << hOut.size() << " hadrons and "
0468                << pOut.size() << " partons after PYTHIA Hadronization";
0469     //hadfile << "NEXT" << endl;
0470 
0471   } // End of positive or negative loop
0472 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!