src/initialstate/PGun.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 //Parton Gun Test
0016
0017 #include "PGun.h"
0018 #include "JetScapeParticles.h"
0019 #include "Pythia8/Pythia.h"
0020
0021 using namespace Jetscape;
0022
0023 // Register the module with the base class
0024 RegisterJetScapeModule<PGun> PGun::reg("PGun");
0025
0026 Pythia8::Pythia PGun::InternalHelperPythia("IntentionallyEmpty", false);
0027
0028 PGun::PGun() : HardProcess() {
0029   fixed_pT = 0;
0030   parID = 21;
0031   SetId("PGun");
0032   VERBOSE(8);
0033 }
0034
0035 PGun::~PGun() { VERBOSE(8); }
0036
0037 void PGun::InitTask() {
0038   JSDEBUG << "Initialize PGun Brick (Test) ...";
0039   VERBOSE(8);
0040
0041   std::string s = GetXMLElementText({"Hard", "PGun", "name"});
0042   JSDEBUG << s << " to be initialized ...";
0043
0044   fixed_pT = GetXMLElementDouble({"Hard", "PGun", "pT"});
0045   JSDEBUG << s << " with fixed pT = " << fixed_pT;
0046   JSINFO << "Parton Gun with fixed pT = " << fixed_pT;
0047
0048   parID = GetXMLElementDouble({"Hard", "PGun", "parID"});
0049   JSINFO << "Parton Gun with parID = " << parID;
0050 }
0051
0052 void PGun::Exec() {
0053   VERBOSE(2) << "Run Hard Process : " << GetId() << " ...";
0054
0055   double p[4], xLoc[4];
0056
0057   double pT, rapidity, phi, pseudorapidity;
0058   double eta_cut = 1.0;
0059   double tempRand;
0060   const double maxN = 1.0 * RAND_MAX;
0061   const double PI = 3.1415926;
0062
0063   double ppx, ppy, ppz, pp0, mass;
0064
0065   // for (int i=0;i<1;i++)
0066   //    {
0067   //      tempRand = rand()/maxN;
0068   //    if(tempRand < 0.25) parID = 21;
0069   //    else if(tempRand < 0.50) parID = 1;
0070   //    else if(tempRand < 0.75) parID = 2;
0071   //    else parID = 3;
0072   //    if (parID != 21) {
0073   //     tempRand = rand()/maxN;
0074   //     if(tempRand < 0.50) parID = -parID;
0075   //      }
0076   //     mass = 0.0;
0077   mass = InternalHelperPythia.particleData.m0(parID);
0078   //JSINFO << BOLDYELLOW << " Mass = " << mass ;
0079   pT = fixed_pT; //max_pT*(rand()/maxN);
0080
0081   phi = 2.0 * PI * (rand() / maxN);
0082   rapidity = 0; //2.0*eta_cut*(rand()/maxN)-eta_cut;
0083   phi = 0.0;
0084
0085   p[1] = pT * cos(phi);
0086   p[2] = pT * sin(phi);
0087   p[3] = sqrt(pT * pT + mass * mass) * sinh(rapidity);
0088   p[0] = sqrt(pT * pT + mass * mass) * cosh(rapidity);
0089
0090   double p_abs = std::sqrt(pT*pT + p[3]*p[3]);
0091   if(std::abs(p_abs - p[3]) > rounding_error){
0092     pseudorapidity = 0.5 * std::log((p_abs + p[3]) / (p_abs - p[3]));
0093   } else {
0094     JSWARN << "Particle in PGun has infinite pseudorapidity.";
0095   }
0096
0097   // Roll for a starting point
0098   // See: https://stackoverflow.com/questions/15039688/random-generator-from-vector-with-probability-distribution-in-c
0099   for (int i = 0; i <= 3; i++) {
0100     xLoc[i] = 0.0;
0101   };
0102
0103   if (!ini) {
0104     VERBOSE(1)
0105         << "No initial state module, setting the starting location to 0. ";
0106   } else {
0107     double x, y;
0108     ini->SampleABinaryCollisionPoint(x, y);
0109     xLoc[1] = x;
0110     xLoc[2] = y;
0111   }
0112
0113   xLoc[1] = 0.0;
0114   xLoc[2] = 0.0;
0115
0116   auto ptn = make_shared<Parton>(0, parID, 0, pT, pseudorapidity, phi, p[0], xLoc);
0117   ptn->set_color((parID > 0) ? 100 : 0);
0118   ptn->set_anti_color(((parID > 0) || (parID == 21)) ? 0 : 101);
0119   ptn->set_max_color(102);
0120   AddParton(ptn);
0121
0122   VERBOSE(8) << GetNHardPartons();
0123 }