src/initialstate/epemGun.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 e+e- -> qqbar processes with Pythia and return the two inital partons with set virtualities
0017
0018 #include "epemGun.h"
0019 #include <sstream>
0020
0021 #define MAGENTA "\033[35m"
0022
0023 using namespace std;
0024
0025 // Register the module with the base class
0026 RegisterJetScapeModule<epemGun> epemGun::reg("epemGun");
0027
0028 epemGun::~epemGun() { VERBOSE(8); }
0029
0030 void epemGun::InitTask() {
0031
0032   JSDEBUG << "Initialize epemGun";
0033   VERBOSE(8);
0034
0035   // Show initialization at INFO level
0036   readString("Init:showProcesses = off");
0037   readString("Init:showChangedSettings = off");
0038   readString("Init:showMultipartonInteractions = off");
0039   readString("Init:showChangedParticleData = off");
0040   if (JetScapeLogger::Instance()->GetInfo()) {
0041     readString("Init:showProcesses = on");
0042     readString("Init:showChangedSettings = on");
0043     readString("Init:showMultipartonInteractions = on");
0044     readString("Init:showChangedParticleData = on");
0045   }
0046
0047   // No event record printout.
0048   readString("Next:numberShowInfo = 0");
0049   readString("Next:numberShowProcess = 0");
0050   readString("Next:numberShowEvent = 0");
0051
0052   // Standard settings
0053   readString("HadronLevel:all = off");
0054   //readString("PartonLevel:ISR = on");
0055   //readString("PartonLevel:MPI = on");
0056   readString("PartonLevel:FSR = off");
0057   //readString("PromptPhoton:all=on");
0058   readString("PDF:lepton = off");
0059   readString("WeakSingleBoson:ffbar2gmz=on"); //Scattering f fbar → gamma^*/Z^0, with full interference between the gamma^* and Z^0
0060   readString("WeakDoubleBoson:all=on"); //Common switch for the group of pair production of gamma^*/Z^0 and W^+-
0061   readString("WeakBosonExchange:all=on"); //Common switch for the group of gamma^*/Z^0 or W^+- exchange between two fermions
0062
0063   //Stuff I added. Ask if I'm allowed to just do this.
0064   readString("23:onMode = off");
0065   readString("23:onIfAny = 1 2 3 4 5");
0066
0067   // For parsing text
0068   stringstream numbf(stringstream::app | stringstream::in | stringstream::out);
0069   numbf.setf(ios::fixed, ios::floatfield);
0070   numbf.setf(ios::showpoint);
0071   numbf.precision(1);
0072   stringstream numbi(stringstream::app | stringstream::in | stringstream::out);
0073
0074   std::string s = GetXMLElementText({"Hard", "epemGun", "name"});
0075   SetId(s);
0076   // cout << s << endl;
0077
0078   // SC: read flag for FSR
0079   //FSR_on = GetXMLElementInt({"Hard", "epemGun", "FSR_on"});
0080   //if (FSR_on)
0081   //  readString("PartonLevel:FSR = on");
0082   //else
0083   //  readString("PartonLevel:FSR = off");
0084
0085   //pTHatMin = GetXMLElementDouble({"Hard", "epemGun", "pTHatMin"});
0086   //pTHatMax = GetXMLElementDouble({"Hard", "epemGun", "pTHatMax"});
0087
0088   //JSINFO << MAGENTA << "epem Gun with FSR_on: " << FSR_on;
0089   //JSINFO << MAGENTA << "epem Gun with " << pTHatMin << " < pTHat < "
0090   //       << pTHatMax;
0091
0092   //numbf.str("PhaseSpace:pTHatMin = ");
0093   //numbf << pTHatMin;
0094   //readString(numbf.str());
0095   //numbf.str("PhaseSpace:pTHatMax = ");
0096   //numbf << pTHatMax;
0097   //readString(numbf.str());
0098
0099   // random seed
0100   // xml limits us to unsigned int :-/ -- but so does 32 bits Mersenne Twist
0101   tinyxml2::XMLElement *RandomXmlDescription = GetXMLElement({"Random"});
0102   readString("Random:setSeed = on");
0103   numbi.str("Random:seed = ");
0104   unsigned int seed = 0;
0105   if (RandomXmlDescription) {
0106     tinyxml2::XMLElement *xmle =
0107         RandomXmlDescription->FirstChildElement("seed");
0108     if (!xmle)
0109       throw std::runtime_error("Cannot parse xml");
0110     xmle->QueryUnsignedText(&seed);
0111   } else {
0112     JSWARN << "No <Random> element found in xml, seeding to 0";
0113   }
0114   VERBOSE(7) << "Seeding pythia to " << seed;
0115   numbi << seed;
0116   readString(numbi.str());
0117
0118   // Species
0119   readString("Beams:idA = 11");
0120   readString("Beams:idB = -11");
0121
0122   // Energy
0123   eCM = GetXMLElementDouble({"Hard", "epemGun", "eCM"});
0124   numbf.str("Beams:eCM = ");
0125   numbf << eCM;
0126   readString(numbf.str());
0127
0128   std::stringstream lines;
0129   lines << GetXMLElementText({"Hard", "epemGun", "LinesToRead"}, false);
0130   int i = 0;
0131   while (std::getline(lines, s, '\n')) {
0132     if (s.find_first_not_of(" \t\v\f\r") == s.npos)
0133       continue; // skip empty lines
0134     VERBOSE(7) << "Also reading in: " << s;
0135     readString(s);
0136   }
0137
0138   // And initialize
0139   if (!init()) { // Pythia>8.1
0140     throw std::runtime_error("Pythia init() failed.");
0141   }
0142
0143   // Initialize random number distribution
0144   ZeroOneDistribution = uniform_real_distribution<double>{0.0, 1.0};
0145
0146 }
0147
0148 void epemGun::Exec() {
0149   VERBOSE(1) << "Run Hard Process : " << GetId() << " ...";
0150   VERBOSE(8) << "Current Event #" << GetCurrentEvent();
0151   //Reading vir_factor from xml for MATTER
0152   double vir_factor = GetXMLElementDouble({"Eloss", "Matter", "vir_factor"});
0153
0154   bool flag62 = false;
0155   vector<Pythia8::Particle> p62;
0156
0157   // sort by pt
0158   struct greater_than_pt {
0159     inline bool operator()(const Pythia8::Particle &p1,
0160                            const Pythia8::Particle &p2) {
0161       return (p1.pT() > p2.pT());
0162     }
0163   };
0164
0165   do {
0166     next();
0167     //event.list();
0168     p62.clear();
0169
0170     for (int parid = 0; parid < event.size(); parid++) {
0171       if (parid < 3)
0172         continue; // 0, 1, 2: total event and beams
0173       Pythia8::Particle &particle = event[parid];
0174
0175       //replacing diquarks with antiquarks (and anti-dq's with quarks)
0176       //the id is set to the heaviest quark in the diquark (except down quark)
0177       //this technically violates baryon number conservation over the entire event
0178       //also can violate electric charge conservation
0179       if( (std::abs(particle.id()) > 1100) && (std::abs(particle.id()) < 6000) && ((std::abs(particle.id())/10)%10 == 0) ){
0180         //if(particle.id() > 0){particle.id( -1*particle.id()/1000 );}
0181         //else{particle.id( particle.id()/1000 );}
0182         particle.id( -1*particle.id()/1000 );       //Changed from previous two lines.
0183       }
0184
0185       //if (!FSR_on) {
0186         // only accept particles after MPI
0187         //if (particle.status() != 62) {continue;}
0188           if ( particle.status()<0 ) {continue;}
0189         // only accept gluons and quarks
0190         // Also accept Gammas to put into the hadron's list
0191         if (fabs(particle.id()) > 6 && (particle.id() != 21 && particle.id() != 22)) {
0192           continue;
0193     }
0194
0195         //if(particle.id()==22) cout<<"########this is a photon!######" <<endl;
0196         // accept
0197       /*} else { // FSR_on true: use Pythia vacuum shower instead of MATTER
0198         if (!particle.isFinal())
0199           continue;
0200         // only accept gluons and quarks
0201         // Also accept Gammas to put into the hadron's list
0202         if (fabs(particle.id()) > 5 &&
0203             (particle.id() != 21 && particle.id() != 22))
0204           continue;
0205       }*/
0206       p62.push_back(particle);
0207     }
0208
0209     // if you want at least 2
0210     if (p62.size() < 2)
0211       continue;
0212     //if ( p62.size() < 1 ) continue;
0213
0214     // Now have all candidates, sort them
0215     // sort by pt
0216     std::sort(p62.begin(), p62.end(), greater_than_pt());
0217     // // check...
0218     // for (auto& p : p62 ) cout << p.pT() << endl;
0219
0220     flag62 = true;
0221
0222   } while (!flag62);
0223
0224   double p[4], xLoc[4];
0225
0226   // This location should come from an initial state
0227   for (int i = 0; i <= 3; i++) {
0228     xLoc[i] = 0.0;
0229   };
0230
0231   // Now determine virtualities and rebalance
0232   // Need back-to-back kinematics with a q-qbar pair
0233   // This seems always to be the case for the epem gun but check nevertheless
0234   if(p62.size() == 2 && std::abs(p62[0].e() - 0.5*eCM) < 0.001 && std::abs(p62[1].e() - 0.5*eCM) < 0.001 && std::abs(p62[0].id()) < 6 && std::abs(p62[1].id()) < 6){
0235
0236     // Virtualities of the two partons
0237     double q1 = 0.;
0238     double q2 = 0.;
0239     const double QS = 0.9;
0240
0241     //Find initial virtuality one parton at a time
0242     for(int pass=0; pass<2; ++pass){
0243
0244       double mass = p62[pass].m0();
0245       // this part is for the case that light quarks are considered massless
0246       /*if(std::abs(p62[pass].id()) < 4){
0247         mass = 0.;
0248       }*/
0249       double max_vir = (0.25*eCM*eCM - mass*mass) * std::abs(vir_factor);
0250       double min_vir = (0.5 * QS * QS ) * (1.0 + std::sqrt(1.0 + 4.0 * mass*mass / (QS*QS)));
0251
0252       double tQ2 = 0.;
0253
0254       if (max_vir <= QS * QS){
0255         tQ2 = 0.0;
0256       }else if(max_vir < min_vir){
0257         tQ2 = QS * QS;
0258       }else{
0259
0260         double numer = 1.0;
0261         double random = ZeroOneDistribution(*GetMt19937Generator());
0262         double ratio = 1.0;
0263         double diff = ratio;
0264         if(random > rounding_error){
0265           diff = (ratio - random) / random;
0266         }
0267
0268         if(max_vir >= (QS*QS / 2.) * (1.0 + std::sqrt(1.0 + 2.0 * mass * mass / (QS*QS / 2.)))){
0269           double g = (QS*QS / 2.) * (1.0 + std::sqrt(1.0 + 2.0 * mass * mass / (QS*QS / 2.)));
0270           numer = exp(-1.0 * (Cf / 2.0 / pi) * sud_val_QG_w_M(mass,(QS*QS / 2.), g, max_vir, 0.5*eCM));
0271         }
0272
0273         if (numer > random){tQ2 = min_vir;}
0274         else{
0275
0276           double t_hi = max_vir;
0277           double t_low = min_vir;
0278           double t_mid = t_low;
0279
0280           double denom = 1.0;
0281
0282           do{
0283               t_mid = 0.5*(t_low + t_hi);
0284
0285               if (t_mid < (QS*QS / 2.) * (1.0 + std::sqrt(1.0 + 2.0 * mass * mass / (QS*QS / 2.)))){
0286                 denom = 1.0;
0287               }else{
0288                 double g = (QS*QS / 2.) * (1.0 + std::sqrt(1.0 + 2.0 * mass * mass / (QS*QS / 2.)));
0289                 denom = exp(-1.0 * (Cf / 2.0 / pi) * sud_val_QG_w_M(mass,(QS*QS / 2.), g, t_mid, 0.5*eCM));
0290               }
0291
0292               ratio = numer / denom;
0293
0294               diff = (ratio - random) / random;
0295
0296               if (diff < 0.0){t_low = t_mid;}
0297               else{t_hi = t_mid;}
0298
0299           }while((abs(diff) > s_approx) && (abs(t_hi - t_low) / t_hi > s_error));
0300
0301           tQ2 = t_mid;
0302
0303         }
0304       }
0305
0306       if(pass==0){q1=sqrt(tQ2);}
0307       else if(pass==1){q2=sqrt(tQ2);}
0308     }
0309
0310     double modm_sq1 = q1*q1 + p62[0].m0()*p62[0].m0();
0311     double modm_sq2 = q2*q2 + p62[1].m0()*p62[1].m0();
0312
0313     if(eCM > sqrt(modm_sq1)+sqrt(modm_sq2)){
0314       // Check viability condition; should always be satisfied in the massless case if virfactor < 1
0315       double pnew = 0.5*sqrt((eCM*eCM-modm_sq1-modm_sq2)*(eCM*eCM-modm_sq1-modm_sq2)-4.*modm_sq1*modm_sq2)/eCM;
0316
0317       auto magnitude = [](const Pythia8::Particle &p) {
0318         return std::sqrt(p.px() * p.px() + p.py() * p.py() + p.pz() * p.pz());
0319       };
0320
0321       double scale1 = pnew/magnitude(p62[0]);
0322       double scale2 = pnew/magnitude(p62[1]);
0323       double e1new = sqrt(pnew*pnew + modm_sq1);
0324       double e2new = sqrt(pnew*pnew + modm_sq2);
0325       p62[0].e(e1new);
0326       p62[0].px(p62[0].px()*scale1);
0327       p62[0].py(p62[0].py()*scale1);
0328       p62[0].pz(p62[0].pz()*scale1);
0329       p62[1].e(e2new);
0330       p62[1].px(p62[1].px()*scale2);
0331       p62[1].py(p62[1].py()*scale2);
0332       p62[1].pz(p62[1].pz()*scale2);
0333     }
0334   }
0335
0336   //give partons to the framework
0337   for (int np = 0; np < p62.size(); ++np) {
0338     Pythia8::Particle &particle = p62.at(np);
0339
0340     VERBOSE(7) << "Adding particle with pid = " << particle.id()
0341                << " at x=" << xLoc[1] << ", y=" << xLoc[2] << ", z=" << xLoc[3];
0342
0343     VERBOSE(7) << "Adding particle with pid = " << particle.id()
0344                << ", pT = " << particle.pT() << ", eta = " << particle.eta()
0345                << ", phi = " << particle.phi() << ", e = " << particle.e();
0346
0347     VERBOSE(7) << " at x=" << xLoc[1] << ", y=" << xLoc[2] << ", z=" << xLoc[3];
0348
0349     auto ptn = make_shared<Parton>(0, particle.id(), 0, particle.pT(), particle.eta(), particle.phi(), particle.e(), xLoc);
0350     ptn->set_color(particle.col());
0351     ptn->set_anti_color(particle.acol());
0352     ptn->set_max_color(1000 * (np + 1));
0353
0354     if(p62.size() == 2 && std::abs(particle.id()) < 6){
0355       double mean_form_time = (2.*ptn->e()) / (ptn->e()*ptn->e()
0356                             - ptn->px()*ptn->px() - ptn->py()*ptn->py()
0357                             - ptn->pz()*ptn->pz() - ptn->restmass()*ptn->restmass()
0358                             + rounding_error) / fmToGeVinv;
0359       ptn->set_form_time(mean_form_time);
0360       ptn->set_mean_form_time();
0361
0362       double velocity[4];
0363       velocity[0] = 1.0;
0364       for (int j = 1; j <= 3; j++) {
0365         velocity[j] = ptn->p(j) / ptn->e();
0366       }
0367       ptn->set_jet_v(velocity);
0368     }
0369     AddParton(ptn);
0370   }
0371
0372   //event.list();
0373
0374   VERBOSE(8) << GetNHardPartons();
0375 }
0376
0377 double epemGun::sud_val_QG_w_M(double M, double h0, double h1, double h2, double E1) {
0378   double val, h, intg, hL, hR, diff, intg_L, intg_R, span;
0379
0380   val = 0.0;
0381
0382   h = (h1 + h2) / 2.0;
0383
0384   span = (h2 - h1) / h2;
0385
0386   val = alpha_s(h) * sud_z_QG_w_M(M, h0, h, E1);
0387
0388   intg = val * (h2 - h1);
0389
0390   hL = (h1 + h) / 2.0;
0391
0392   intg_L = alpha_s(hL) * sud_z_QG_w_M(M, h0, hL, E1) * (h - h1);
0393
0394   hR = (h + h2) / 2.0;
0395
0396   intg_R = alpha_s(hR) * sud_z_QG_w_M(M, h0, hR, E1) * (h2 - h);
0397
0398   diff = std::abs((intg_L + intg_R - intg) / intg);
0399
0400   if ((diff > approx) || (span > error)) {
0401     intg = sud_val_QG_w_M(M, h0, h1, h, E1) +
0402            sud_val_QG_w_M(M, h0, h, h2, E1);
0403   }
0404
0405   return intg;
0406 }
0407
0408 double epemGun::sud_z_QG_w_M(double M, double cg, double cg1, double E2){
0409   // this part is for the case that light quarks are considered massless
0410   /*if(M < 1.){
0411     if (cg1 < 2.0 * cg) {
0412       return 0.0;
0413     };
0414
0415     double t2 = std::pow(cg1, 2);
0416     double t6 = std::log(cg);
0417     double t10 = std::abs(cg - cg1);
0418     double t11 = std::log(t10);
0419     double t17 = -1.0 / t2 * (3.0 * cg1 - 6.0 * cg + 4.0 * t6 * cg1 - 4.0 * t11 * cg1) /
0420         2.0;
0421     if (t17 < 0.0) {
0422       cerr << "ERROR: t17 negative in sud_z_QG_w_M = " << t17 << endl;
0423       throw std::runtime_error("ERROR: medium contribution negative in sud_z_QG_w_M");
0424     }
0425     return t17;
0426   }*/
0427
0428   if (cg1 < 2.0 * cg + M * M / (1.0 + M * M / cg1)) {
0429     JSINFO << MAGENTA << " returning with cg, cg1 = " << cg << "   " << cg1
0430            << "    " << E_minimum << "  " << E2;
0431     return M * M;
0432   };
0433
0434   double t1 = 1.0 / cg1;
0435   double t2 = t1 * cg;
0436   double t4 = std::pow(1.0 - t2, 2.0);
0437   double t7 = std::log(t2);
0438   double t9 = M * M;
0439   double t10 = t1 * t9;
0440   double t13 = 1.0 / (t10 + 1.0) * t10;
0441   double t15 = std::pow(t2 + t13, 2.0);
0442   double t18 = std::log(1.0 - t2 - t13);
0443   double t21 = t1 * (-t4 / 2.0 - 1.0 + 2.0 * t2 - 2.0 * t7 + t15 / 2.0 + t13 +
0444                      2.0 * t18);
0445
0446   if (t21 < 0.0) {
0447     cerr << "ERROR: t21 negative in sud_z_QG_w_M = " << t21 << endl;
0448     throw std::runtime_error("ERROR: medium contribution negative in sud_z_QG_w_M");
0449   }
0450
0451   return t21;
0452 }
0453
0454 double epemGun::alpha_s(double q2) {
0455   double a, L2, q24, c_nf;
0456
0457   L2 = std::pow(Lambda_QCD, 2);
0458
0459   q24 = q2 / 4.0;
0460
0461   c_nf = nf;
0462
0463   if (q24 > 4.0) {
0464     c_nf = 4;
0465   }
0466
0467   if (q24 > 64.0) {
0468     c_nf = 5;
0469   }
0470
0471   if (q24 > L2) {
0472     a = 12.0 * pi / (11.0 * Nc - 2.0 * c_nf) / std::log(q24 / L2);
0473   } else {
0474     JSWARN << " alpha too large ";
0475     a = 0.6;
0476   }
0477
0478   return a;
0479   }