src/initialstate/PythiaGun.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 // Create a pythia collision at a specified point and return the two inital hard partons
0017
0018 #include "PythiaGun.h"
0019 #include <sstream>
0020 #include <iostream>
0021 #include <fstream>
0022 #define MAGENTA "\033[35m"
0023
0024 using namespace std;
0025
0026 // Register the module with the base class
0027 RegisterJetScapeModule<PythiaGun> PythiaGun::reg("PythiaGun");
0028
0029 PythiaGun::~PythiaGun() { VERBOSE(8); }
0030
0031 void PythiaGun::InitTask() {
0032
0033   JSDEBUG << "Initialize PythiaGun";
0034   VERBOSE(8);
0035
0036   // Show initialization at INFO level
0037   readString("Init:showProcesses = off");
0038   readString("Init:showChangedSettings = off");
0039   readString("Init:showMultipartonInteractions = off");
0040   readString("Init:showChangedParticleData = off");
0041   if (JetScapeLogger::Instance()->GetInfo()) {
0042     readString("Init:showProcesses = on");
0043     readString("Init:showChangedSettings = on");
0044     readString("Init:showMultipartonInteractions = on");
0045     readString("Init:showChangedParticleData = on");
0046   }
0047
0048   // No event record printout.
0049   readString("Next:numberShowInfo = 0");
0050   readString("Next:numberShowProcess = 0");
0051   readString("Next:numberShowEvent = 0");
0052
0053   // For parsing text
0054   stringstream numbf(stringstream::app | stringstream::in | stringstream::out);
0055   numbf.setf(ios::fixed, ios::floatfield);
0056   numbf.setf(ios::showpoint);
0057   numbf.precision(1);
0058   stringstream numbi(stringstream::app | stringstream::in | stringstream::out);
0059
0060   std::string s = GetXMLElementText({"Hard", "PythiaGun", "name"});
0061   SetId(s);
0062   // cout << s << endl;
0063
0064   //other Pythia settings
0065   readString("HadronLevel:Decay = off");
0066   readString("HadronLevel:all = off");
0067   readString("PartonLevel:ISR = on");
0068   readString("PartonLevel:MPI = on");
0069   //readString("PartonLevel:FSR = on");
0070   readString("PromptPhoton:all=on");
0071   readString("WeakSingleBoson:all=off");
0072   readString("WeakDoubleBoson:all=off");
0073
0074   pTHatMin = GetXMLElementDouble({"Hard", "PythiaGun", "pTHatMin"});
0075   pTHatMax = GetXMLElementDouble({"Hard", "PythiaGun", "pTHatMax"});
0076
0077   if(pTHatMin < 0.01){ //assuming low bin where softQCD should be used
0078     //running softQCD - inelastic nondiffrative (min-bias)
0079     readString("HardQCD:all = off");
0080     readString("SoftQCD:nonDiffractive = on");
0081     softQCD = true;
0082   }
0083   else{ //running normal hardQCD
0084     readString("HardQCD:all = on"); // will repeat this line in the xml for demonstration
0085     //  readString("HardQCD:gg2ccbar = on"); // switch on heavy quark channel
0086     //readString("HardQCD:qqbar2ccbar = on");
0087     numbf.str("PhaseSpace:pTHatMin = "); numbf << pTHatMin; readString(numbf.str());
0088     numbf.str("PhaseSpace:pTHatMax = "); numbf << pTHatMax; readString(numbf.str());
0089       softQCD = false;
0090   }
0091
0092   // SC: read flag for FSR
0093   FSR_on = GetXMLElementInt({"Hard", "PythiaGun", "FSR_on"});
0094   if (FSR_on)
0095     readString("PartonLevel:FSR = on");
0096   else
0097     readString("PartonLevel:FSR = off");
0098
0099   JSINFO << MAGENTA << "Pythia Gun with FSR_on: " << FSR_on;
0100   JSINFO << MAGENTA << "Pythia Gun with " << pTHatMin << " < pTHat < " << pTHatMax;
0101
0102   // random seed
0103   // xml limits us to unsigned int :-/ -- but so does 32 bits Mersenne Twist
0104   tinyxml2::XMLElement *RandomXmlDescription = GetXMLElement({"Random"});
0105   readString("Random:setSeed = on");
0106   numbi.str("Random:seed = ");
0107   unsigned int seed = 0;
0108   if (RandomXmlDescription) {
0109     tinyxml2::XMLElement *xmle =
0110         RandomXmlDescription->FirstChildElement("seed");
0111     if (!xmle)
0112       throw std::runtime_error("Cannot parse xml");
0113     xmle->QueryUnsignedText(&seed);
0114   } else {
0115     JSWARN << "No <Random> element found in xml, seeding to 0";
0116   }
0117   VERBOSE(7) << "Seeding pythia to " << seed;
0118   numbi << seed;
0119   readString(numbi.str());
0120
0121   // Species
0122   readString("Beams:idA = 2212");
0123   readString("Beams:idB = 2212");
0124
0125   // Energy
0126   eCM = GetXMLElementDouble({"Hard", "PythiaGun", "eCM"});
0127   numbf.str("Beams:eCM = ");
0128   numbf << eCM;
0129   readString(numbf.str());
0130
0131   //Reading vir_factor from xml for MATTER
0132   vir_factor = GetXMLElementDouble({"Eloss", "Matter", "vir_factor"});
0133   softMomentumCutoff = GetXMLElementDouble({"Hard", "PythiaGun", "softMomentumCutoff"});
0134
0135   std::stringstream lines;
0136   lines << GetXMLElementText({"Hard", "PythiaGun", "LinesToRead"}, false);
0137   int i = 0;
0138   while (std::getline(lines, s, '\n')) {
0139     if (s.find_first_not_of(" \t\v\f\r") == s.npos)
0140       continue; // skip empty lines
0141     VERBOSE(7) << "Also reading in: " << s;
0142     readString(s);
0143   }
0144
0145   // And initialize
0146   if (!init()) { // Pythia>8.1
0147     throw std::runtime_error("Pythia init() failed.");
0148   }
0149
0150     std::ofstream sigma_printer;
0151     sigma_printer.open(printer, std::ios::trunc);
0152
0153
0154 }
0155
0156 void PythiaGun::Exec() {
0157   VERBOSE(1) << "Run Hard Process : " << GetId() << " ...";
0158   VERBOSE(8) << "Current Event #" << GetCurrentEvent();
0159
0160   bool flag62 = false;
0161   vector<Pythia8::Particle> p62;
0162
0163   // sort by pt
0164   struct greater_than_pt {
0165     inline bool operator()(const Pythia8::Particle &p1,
0166                            const Pythia8::Particle &p2) {
0167       return (p1.pT() > p2.pT());
0168     }
0169   };
0170
0171   do {
0172     next();
0173     p62.clear();
0174       if (!printer.empty()){
0175             std::ofstream sigma_printer;
0176             sigma_printer.open(printer, std::ios::out | std::ios::app);
0177
0178             sigma_printer << "sigma = " << GetSigmaGen() << " Err =  " << GetSigmaErr() << endl ;
0179             //sigma_printer.close();
0180
0181 //      JSINFO << BOLDYELLOW << " sigma = " << GetSigmaGen() << " sigma err = " << GetSigmaErr() << " printer = " << printer << " is " << sigma_printer.is_open() ;
0182     };
0183
0184     // pTarr[0]=0.0; pTarr[1]=0.0;
0185     // pindexarr[0]=0; pindexarr[1]=0;
0186
0187     for (int parid = 0; parid < event.size(); parid++) {
0188       if (parid < 3)
0189         continue; // 0, 1, 2: total event and beams
0190       Pythia8::Particle &particle = event[parid];
0191
0192       //replacing diquarks with antiquarks (and anti-dq's with quarks)
0193       //the id is set to the heaviest quark in the diquark (except down quark)
0194       //this technically violates baryon number conservation over the entire event
0195       //also can violate electric charge conservation
0196       if( (std::abs(particle.id()) > 1100) && (std::abs(particle.id()) < 6000) && ((std::abs(particle.id())/10)%10 == 0) ){
0197         if(particle.id() > 0){particle.id( -1*particle.id()/1000 );}
0198         else{particle.id( particle.id()/1000 );}
0199       }
0200
0201       if (!FSR_on) {
0202         // only accept particles after MPI
0203         if (particle.status() != 62)
0204           continue;
0205         // only accept gluons and quarks
0206         // Also accept Gammas to put into the hadron's list
0207         if (fabs(particle.id()) > 5 &&
0208             (particle.id() != 21 && particle.id() != 22))
0209           continue;
0210
0211         // reject rare cases of very soft particles that don't have enough e to get
0212         // reasonable virtuality
0213         if (vir_factor > 0. && (particle.pT() < softMomentumCutoff)) {
0214           // this cutoff was 1.0/sqrt(vir_factor) in versions < 3.6
0215           continue;
0216         } else if(vir_factor < 0. && (particle.pAbs() < softMomentumCutoff)) {
0217           continue;
0218         } else if(vir_factor < rounding_error) {
0219           JSWARN << "vir_factor should not be zero.";
0220           exit(1);
0221         }
0222
0223         //if(particle.id()==22) cout<<"########this is a photon!######" <<endl;
0224         // accept
0225       } else { // FSR_on true: use Pythia vacuum shower instead of MATTER
0226         if (!particle.isFinal())
0227           continue;
0228         // only accept gluons and quarks
0229         // Also accept Gammas to put into the hadron's list
0230         if (fabs(particle.id()) > 5 &&
0231             (particle.id() != 21 && particle.id() != 22))
0232           continue;
0233       }
0234       p62.push_back(particle);
0235     }
0236
0237     // if you want at least 2
0238     if (p62.size() < 2)
0239       continue;
0240     //if ( p62.size() < 1 ) continue;
0241
0242     // Now have all candidates, sort them
0243     // sort by pt
0244     std::sort(p62.begin(), p62.end(), greater_than_pt());
0245     // // check...
0246     // for (auto& p : p62 ) cout << p.pT() << endl;
0247
0248     //skipping event if softQCD is on & pThat exceeds max (where next bin is HardQCD with this as pThatmin)
0249     if(softQCD && (info.pTHat() >= pTHatMax)){continue;}
0250
0251     flag62 = true;
0252
0253   } while (!flag62);
0254
0255   double p[4], xLoc[4];
0256
0257   // This location should come from an initial state
0258   for (int i = 0; i <= 3; i++) {
0259     xLoc[i] = 0.0;
0260   };
0261
0262   // // Roll for a starting point
0263   // // See: https://stackoverflow.com/questions/15039688/random-generator-from-vector-with-probability-distribution-in-c
0264   // std::random_device device;
0265   // std::mt19937 engine(device()); // Seed the random number engine
0266
0267   if (!ini) {
0268     VERBOSE(1) << "No initial state module, setting the starting location to "
0269                   "0. Make sure to add e.g. trento before PythiaGun.";
0270   } else {
0271     double x, y;
0272     ini->SampleABinaryCollisionPoint(x, y);
0273     xLoc[1] = x;
0274     xLoc[2] = y;
0275   }
0276
0277   // Loop through particles
0278
0279   // Only top two
0280   //for(int np = 0; np<2; ++np){
0281
0282   // Accept them all
0283
0284   int hCounter = 0;
0285   for (int np = 0; np < p62.size(); ++np) {
0286     Pythia8::Particle &particle = p62.at(np);
0287
0288     VERBOSE(7) << "Adding particle with pid = " << particle.id()
0289                << " at x=" << xLoc[1] << ", y=" << xLoc[2] << ", z=" << xLoc[3];
0290
0291     VERBOSE(7) << "Adding particle with pid = " << particle.id()
0292                << ", pT = " << particle.pT() << ", eta = " << particle.eta()
0293                << ", phi = " << particle.phi() << ", e = " << particle.e();
0294
0295     VERBOSE(7) << " at x=" << xLoc[1] << ", y=" << xLoc[2] << ", z=" << xLoc[3];
0296
0297     auto ptn = make_shared<Parton>(0, particle.id(), 0, particle.pT(), particle.eta(),particle.phi(), particle.e(), xLoc);
0298     ptn->set_color(particle.col());
0299     ptn->set_anti_color(particle.acol());
0300     ptn->set_max_color(1000 * (np + 1));
0301     AddParton(ptn);
0302   }
0303
0304   VERBOSE(8) << GetNHardPartons();
0305 }