src/hadronization/HybridHadronization.h

0001 /*******************************************************************************
0002  * Copyright (c) The JETSCAPE Collaboration, 2019
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 #ifndef HYBRIDHADRONIZATION_H
0017 #define HYBRIDHADRONIZATION_H
0018
0019 #include "HadronizationModule.h"
0020 #include "JetScapeLogger.h"
0021 #include "Pythia8/Pythia.h"
0022
0023 #include <cmath>
0024 #include <random>
0025 #include <vector>
0026
0027 using namespace Jetscape;
0028
0029 class HybridHadronization : public HadronizationModule<HybridHadronization>
0030 {
0031  public:
0032
0033   HybridHadronization();
0034   virtual ~HybridHadronization();
0035
0036   void Init();
0037   void DoHadronization(vector<vector<shared_ptr<Parton>>>& shower, vector<shared_ptr<Hadron>>& hOut, vector<shared_ptr<Parton>>& pOut);
0038   void WriteTask(weak_ptr<JetScapeWriter> w);
0039
0040  private:
0041   // Allows the registration of the module so that it is available to be used by the Jetscape framework.
0042   static RegisterJetScapeModule<HybridHadronization> reg;
0043
0044   double SigM2_calc(double R2chg, double qm1, double qm2, double qq1, double qq2);
0045   double SigBR2_calc(double R2chg, double qm1, double qm2, double qm3, double qq1, double qq2, double qq3);
0046   double SigBL2_calc(double SigBR2, double qm1, double qm2, double qm3);
0047
0048   //double sigma_pi, sigma_k, sigma_nuc, maxE_level, gmax, xmq, xms, hbarc, dist2cut, sh_recofactor, th_recofactor, SigRB, SigLB;
0049   double maxM_level, maxB_level, gmax, xmq, xms, xmc, xmb, hbarc, dist2cut, sh_recofactor, th_recofactor, p_fake, had_prop, part_prop, delta_t, hydro_Tc, eta_max_boost_inv;
0050   int number_p_fake;
0051   double SigNucR2,SigNucL2,SigOmgR2,SigOmgL2,SigXiR2,SigXiL2,SigSigR2,SigSigL2,
0052     SigOcccR2,SigOcccL2,SigOccR2,SigOccL2,SigXiccR2,SigXiccL2,SigOcR2,SigOcL2,SigXicR2,SigXicL2,SigSigcR2,SigSigcL2,
0053     SigObbbR2,SigObbbL2,SigObbcR2,SigObbcL2,SigObbR2,SigObbL2,SigXibbR2,SigXibbL2,
0054     SigObccR2,SigObccL2,SigObcR2,SigObcL2,SigXibcR2,SigXibcL2,SigObR2,SigObL2,SigXibR2,SigXibL2,SigSigbR2,SigSigbL2;
0055   double SigPi2, SigPhi2, SigK2, SigJpi2, SigDs2, SigD2, SigUps2, SigBc2, SigB2;
0056   bool inbrick, inhydro; int nreusehydro; double brickL, brickT;
0057   const double pi = 3.1415926535897932384626433832795;
0058   int attempts_max;
0059   unsigned int rand_seed;
0060   int reco_hadrons_pythia;
0061   int additional_pythia_particles;
0062   bool goldstonereco;
0063   bool torder_reco;
0064   bool afterburner_frag_hadrons = false;
0065   std::string pythia_decays;
0066
0067   //variables for recombination color structure
0068   vector<vector<vector<int>>> Tempjunctions; // vector of all tempjunctions
0069   vector<vector<int>> JunctionInfo; // vector of one junction's color tag and particle info
0070   vector<int> IdColInfo1; vector<int> IdColInfo2; vector<int> IdColInfo3; vector<int> IdColInfo4;
0071
0072   std::mt19937_64 eng; //RNG - Mersenne Twist - 64 bit
0073   double ran();
0074
0075   //4-vector boost
0076   static FourVector HHboost(FourVector B, FourVector vec_in){
0077     double xlam[4][4], beta2;
0078         beta2 = B.x()*B.x() + B.y()*B.y() + B.z()*B.z();
0079         B.Set(B.x(),B.y(),B.z(),1./(sqrt(1. - beta2)));
0080
0081     double beta2inv;
0082     if(beta2 > 0.){beta2inv = 1./beta2;}
0083     else{beta2inv = 0.;}
0084
0085     xlam[0][0] = B.t();
0086     xlam[0][1] = -B.t()*B.x();
0087     xlam[0][2] = -B.t()* B.y();
0088     xlam[0][3] = -B.t()*B.z();
0089     xlam[1][0] = xlam[0][1];
0090     xlam[1][1] = 1.+(B.t() - 1.)*(B.x()*B.x())*beta2inv;
0091     xlam[1][2] = (B.t()-1.)*B.x()*B.y()*beta2inv;
0092     xlam[1][3] = (B.t()-1.)*B.x()*B.z()*beta2inv;
0093     xlam[2][0] = xlam[0][2];
0094     xlam[2][1] = xlam[1][2];
0095     xlam[2][2] = 1.+(B.t()-1.)*(B.y()*B.y())*beta2inv;
0096     xlam[2][3] = (B.t()-1.)*B.y()*B.z()*beta2inv;
0097     xlam[3][0] = xlam[0][3];
0098     xlam[3][1] = xlam[1][3];
0099     xlam[3][2] = xlam[2][3];
0100     xlam[3][3] = 1.+(B.t()-1.)*(B.z()*B.z())*beta2inv;
0101
0102     double t_out = vec_in.t()*xlam[0][0] + vec_in.x()*xlam[0][1] + vec_in.y()*xlam[0][2] + vec_in.z()*xlam[0][3];
0103     double x_out = vec_in.t()*xlam[1][0] + vec_in.x()*xlam[1][1] + vec_in.y()*xlam[1][2] + vec_in.z()*xlam[1][3];
0104     double y_out = vec_in.t()*xlam[2][0] + vec_in.x()*xlam[2][1] + vec_in.y()*xlam[2][2] + vec_in.z()*xlam[2][3];
0105     double z_out = vec_in.t()*xlam[3][0] + vec_in.x()*xlam[3][1] + vec_in.y()*xlam[3][2] + vec_in.z()*xlam[3][3];
0106     FourVector vec_out(x_out,y_out,z_out,t_out);
0107
0108     return vec_out;
0109   }
0110
0111   //3-vec (4-vec w/3 components) diff^2 function
0112   static double dif2(FourVector vec1, FourVector vec2){return (vec2.x()-vec1.x())*(vec2.x()-vec1.x()) + (vec2.y()-vec1.y())*(vec2.y()-vec1.y()) + (vec2.z()-vec1.z())*(vec2.z()-vec1.z());}
0113
0114   //parton class
0115   //class parton{
0116   class HHparton : public Parton{
0117   protected:
0118     //is shower/thermal originating parton; has been used; is a decayed gluon; is a remnant parton; used for reco; used for stringfrag; is endpoint of string
0119     bool is_shower_, is_thermal_, is_used_, is_decayedglu_, is_remnant_, used_reco_, used_str_, is_strendpt_, used_junction_, is_fakep_;
0120     //id: particle id, ?orig: particle origin(shower, thermal...)?, par: parent parton (perm. set for gluon decays), status: status of particle (is being considered in loop 'i')
0121     //string_id: string identifier, pos_str: position of parton in string, endpt_id: denotes which endpoint this parton is (of the string), if it is one
0122     int alt_id_, orig_, par_, string_id_, pos_str_, endpt_id_, sibling_, PY_par1_, PY_par2_, PY_dau1_, PY_dau2_, PY_stat_, PY_origid_, PY_tag1_, PY_tag2_, PY_tag3_;
0123     //parton mass
0124     double alt_mass_;
0125
0126   public:
0127     //default constructor
0128     HHparton() : Parton::Parton(1,1,1,0.,0.,0.,0.) {
0129         is_shower_ = false; is_thermal_ = false; is_used_ = false; is_decayedglu_ = false; is_remnant_ = false; used_reco_ = false; used_str_ = false; is_strendpt_ = false; used_junction_ = false; is_fakep_ = false;
0130         alt_id_ = 0; orig_ = 0; par_ = -1; string_id_ = 0; pos_str_ = 0; endpt_id_ = 0; sibling_ = 0; alt_mass_ = 0.;
0131         PY_origid_ = 0; PY_par1_ = -1; PY_par2_ = -1; PY_dau1_ = -1; PY_dau2_ = -1; PY_stat_ = 23; /*PY_stat_ = 11;*/ PY_tag1_ = 0; PY_tag2_ = 0; PY_tag3_ = 0;
0132         set_color(0); set_anti_color(0); set_stat(0);
0133     }
0134
0135     //getter functions
0136     double  x() {return x_in().x();} double  y() {return x_in().y();} double  z() {return x_in().z();} double x_t() {return x_in().t();}
0137     double px() {return p_in().x();} double py() {return p_in().y();} double pz() {return p_in().z();} double   e() {return p_in().t();}
0138     bool is_shower() {return is_shower_;} bool is_thermal() {return is_thermal_;} bool is_used() {return is_used_;} bool is_decayedglu() {return is_decayedglu_;}
0139     bool is_remnant() {return is_remnant_;} bool used_reco() {return used_reco_;} bool used_str() {return used_str_;} bool is_strendpt() {return is_strendpt_;}
0140     bool used_junction() {return used_junction_;} bool is_fakeparton() {return is_fakep_;}
0141     int id() {return alt_id_;} int orig() {return orig_;} int par() {return par_;} int string_id() const {return string_id_;} int pos_str() const {return pos_str_;}
0142     int endpt_id() {return endpt_id_;} int sibling() {return sibling_;} int PY_par1() {return PY_par1_;} int PY_par2() {return PY_par2_;}
0143     int PY_dau1() {return PY_dau1_;} int PY_dau2() {return PY_dau2_;} int PY_stat() {return PY_stat_;} int PY_origid() {return PY_origid_;}
0144     int PY_tag1() {return PY_tag1_;} int PY_tag2() {return PY_tag2_;} int PY_tag3() {return PY_tag3_;}
0145     double mass() {return alt_mass_;}
0146     FourVector pos() {return x_in();}
0147     FourVector P()   {return p_in();}
0148     int status() {return pstat();} int col() {return color();} int acol() {return anti_color();}
0149
0150     //setter functions
0151     void  x(double val) {x_in_.Set(val,x_in().y(),x_in().z(),x_in().t());}
0152     void  y(double val) {x_in_.Set(x_in().x(),val,x_in().z(),x_in().t());}
0153     void  z(double val) {x_in_.Set(x_in().x(),x_in().y(),val,x_in().t());}
0154     void x_t(double val){x_in_.Set(x_in().x(),x_in().y(),x_in().z(),val);}
0155     void pos(FourVector val) {x_in_.Set(val.x(),val.y(),val.z(),val.t());}
0156     void px(double val) {reset_momentum(val,py(),pz(),e());}
0157     void py(double val) {reset_momentum(px(),val,pz(),e());}
0158     void pz(double val) {reset_momentum(px(),py(),val,e());}
0159     void  e(double val) {reset_momentum(px(),py(),pz(),val);}
0160     void  P(FourVector val) {reset_momentum(val);}
0161
0162     void is_shower(bool val) {is_shower_ = val;} void is_thermal(bool val) {is_thermal_ = val;} void is_used(bool val) {is_used_ = val;} void is_decayedglu(bool val) {is_decayedglu_ = val;}
0163     void is_remnant(bool val) {is_remnant_ = val;} void used_reco(bool val) {used_reco_ = val;} void used_str(bool val) {used_str_ = val;} void is_strendpt(bool val) {is_strendpt_ = val;}
0164     void used_junction(bool val) {used_junction_ = val;} void is_fakeparton(bool val) {is_fakep_ = val;}
0165     void id(int val) {alt_id_ = val;} void orig(int val) {orig_ = val;} void par(int val) {par_ = val;}
0166     void string_id(int val) {string_id_ = val;} void pos_str(int val) {pos_str_ = val;} void endpt_id(int val) {endpt_id_ = val;} void sibling(int val) {sibling_ = val;}
0167     void PY_par1(int val) {PY_par1_ = val;} void PY_par2(int val) {PY_par2_ = val;} void PY_dau1(int val) {PY_dau1_ = val;} void PY_dau2(int val) {PY_dau2_ = val;}
0168     void PY_stat(int val) {PY_stat_ = val;} void PY_origid(int val) {PY_origid_ = val;}
0169     void PY_tag1(int val) {PY_tag1_ = val;} void PY_tag2(int val) {PY_tag2_ = val;} void PY_tag3(int val) {PY_tag3_ = val;}
0170     void mass(double val) {alt_mass_ = val;}
0171     void status(int val) {set_stat(val);} void col(int val) {set_color(val);} void acol(int val) {set_anti_color(val);}
0172
0173     //boost functions
0174     FourVector boost_P(FourVector B){return HHboost(B,p_in());}
0175     FourVector boost_P(double vx, double vy, double vz){FourVector B(vx,vy,vz,0.); return HHboost(B,p_in());}
0176     FourVector boost_pos(FourVector B){return HHboost(B,x_in());}
0177     FourVector boost_pos(double vx, double vy, double vz){FourVector B(vx,vy,vz,0.); return HHboost(B,x_in());}
0178
0179     double pDif2(HHparton comp){return dif2(p_in(),comp.p_in());}
0180     double posDif2(HHparton comp){return dif2(x_in(),comp.x_in());}
0181
0182   };
0183
0184     //hadron class
0185   class HHhadron : public Hadron{
0186   protected:
0187     bool is_excited_, is_shsh_, is_shth_, is_thth_, is_recohad_, is_strhad_, is_final_;
0188     //id: particle id, ?orig: particle origin(sh-sh, sh-th, th-th, recombined, stringfragmented...),
0189     //par_str: parent string number, parh: parent hadron (decayed), status: status of particle (final, used, decayed, etc...)
0190     int orig_, parstr_, parh_;
0191     //hadron mass
0192     double alt_mass_;
0193     //vector of partonic parents
0194     //std::vector<int> parents;
0195
0196   public:
0197     //default constructor
0198     HHhadron() : Hadron::Hadron(1,1,1,0.,0.,0.,0.) {
0199         is_excited_ = false; is_shsh_ = false; is_shth_ = false; is_thth_ = false; is_recohad_ = false; is_strhad_ = false; is_final_ = false;
0200         orig_ = 0; parstr_ = 0; parh_ = -1; alt_mass_ = 0.; set_stat(0);
0201     }
0202
0203     //vector of partonic parents
0204     std::vector<int> parents;
0205
0206     //getter/setter for parents
0207     int par(int i){if((i>=0) && (i<parents.size())){return parents[i];}else{return 999999;}}
0208     void add_par(int i){parents.push_back(i);}
0209
0210     //vector of colors (from partons that formed it)
0211     std::vector<int> cols;
0212
0213     //getter/setter for colors
0214     int col(int i){if((i>=0) && (i<cols.size())){return cols[i];}else{return -1;}}
0215     void add_col(int i){cols.push_back(i);}
0216
0217     //getter functions
0218     double  x() {return x_in().x();} double  y() {return x_in().y();} double  z() {return x_in().z();} double x_t() {return x_in().t();}
0219     double px() {return p_in().x();} double py() {return p_in().y();} double pz() {return p_in().z();} double   e() {return p_in().t();}
0220     bool is_excited() {return is_excited_;} bool is_shsh() {return is_shsh_;} bool is_shth() {return is_shth_;} bool is_thth() {return is_thth_;}
0221     bool is_recohad() {return is_recohad_;} bool is_strhad() {return is_strhad_;} bool is_final() {return is_final_;}
0222     int orig() {return orig_;} int parstr() {return parstr_;} int parh() {return parh_;}
0223     double mass() {return alt_mass_;}
0224     FourVector pos() {return x_in();}
0225     FourVector P()   {return p_in();}
0226     int id() {return pid();} int status() {return pstat();}
0227
0228     //setter functions
0229     void  x(double val) {x_in_.Set(val,x_in().y(),x_in().z(),x_in().t());}
0230     void  y(double val) {x_in_.Set(x_in().x(),val,x_in().z(),x_in().t());}
0231     void  z(double val) {x_in_.Set(x_in().x(),x_in().y(),val,x_in().t());}
0232     void x_t(double val){x_in_.Set(x_in().x(),x_in().y(),x_in().z(),val);}
0233     void pos(FourVector val) {x_in_.Set(val.x(),val.y(),val.z(),val.t());}
0234     void px(double val) {reset_momentum(val,py(),pz(),e());}
0235     void py(double val) {reset_momentum(px(),val,pz(),e());}
0236     void pz(double val) {reset_momentum(px(),py(),val,e());}
0237     void  e(double val) {reset_momentum(px(),py(),pz(),val);}
0238     void  P(FourVector val) {reset_momentum(val);}
0239
0240     void is_excited(bool val) {is_excited_ = val;} void is_shsh(bool val) {is_shsh_ = val;} void is_shth(bool val) {is_shth_ = val;} void is_thth(bool val) {is_thth_ = val;}
0241     void is_recohad(bool val) {is_recohad_ = val;} void is_strhad(bool val) {is_strhad_ = val;} void is_final(bool val) {is_final_ = val;}
0242     void orig(int val) {orig_ = val;} void parstr(int val) {parstr_ = val;} void parh(int val) {parh_ = val;}
0243     void mass(double val) {alt_mass_ = val;}
0244     void id(int val) {set_id(val);} void status(int val) {set_stat(val);}
0245
0246     //Lorentz boost functions
0247     FourVector boost_P(FourVector B){return HHboost(B,p_in());}
0248     FourVector boost_P(double vx, double vy, double vz){FourVector B(vx,vy,vz,0.); return HHboost(B,p_in());}
0249     FourVector boost_pos(FourVector B){return HHboost(B,x_in());}
0250     FourVector boost_pos(double vx, double vy, double vz){FourVector B(vx,vy,vz,0.); return HHboost(B,x_in());}
0251   };
0252
0253     //class holding a collection of partons
0254   class parton_collection{
0255   public:
0256     //the actual collection of partons
0257     std::vector<HHparton> partons;
0258     //default constructor
0259     parton_collection() {partons.clear();}
0260     //constructor given a single initial parton
0261     parton_collection(HHparton par) {partons.push_back(par);}
0262     //add a parton to the collection
0263     void add(HHparton par) {partons.push_back(par);}
0264     //add a collection of partons to the collection
0265     void add(parton_collection pars) {for (int i=0; i<pars.num(); ++i){partons.push_back(pars[i]);}}
0266     //get the number of partons in the collection
0267     int num() {return partons.size();}
0268     //empty the collection
0269     void clear() {partons.clear();}
0270     //insert a parton at position i
0271     void insert(int i, HHparton newparton) {partons.insert(partons.begin()+i, newparton);}
0272     //remove the parton at position i
0273     void remove(int i) {partons.erase(partons.begin()+i);}
0274     //swap partons at positions i, j
0275     void swap(int i, int j) {if((i >= 0) && (i < partons.size()) && (j >= 0) && (j < partons.size())) {std::swap(partons[i],partons[j]);}}
0276     //random access iterator pointing to first element in partons
0277     std::vector<HHparton>::iterator begin() {return partons.begin();}
0278     //random access iterator pointing to last element in partons
0279     std::vector<HHparton>::iterator end() {return partons.end();}
0280
0281     //overloading a few operators; ++, --, []
0282     //++ adds an empty particle, -- removes last particle, [i] allows access to i'th particle directly
0283 /*  parton_collection& operator++(){HHparton par; partons.push_back(par);}
0284     parton_collection operator++(int){
0285         parton_collection temp(*this);
0286         ++(*this);
0287         return temp;
0288     }
0289     parton_collection& operator--(){partons.pop_back();}
0290     parton_collection operator--(int){
0291         parton_collection temp(*this);
0292         --(*this);
0293         return temp;
0294     }*/
0295     HHparton& operator[](int i) {return partons[i];}
0296     const HHparton& operator[](int i) const {return partons[i];}
0297   };
0298
0299     //class holding a collection of hadrons
0300   class hadron_collection{
0301   public:
0302     //the actual collection of hadrons
0303     std::vector<HHhadron> hadrons;
0304     //default constructor
0305     hadron_collection() {hadrons.clear();}
0306     //constructor given an initial hadron
0307     hadron_collection(HHhadron had) {hadrons.push_back(had);}
0308     //add a hadron to the collection
0309     void add(HHhadron had) {hadrons.push_back(had);}
0310     //add a collection of partons to the collection
0311     void add(hadron_collection hads) {for (int i=0; i<hads.num(); ++i){hadrons.push_back(hads[i]);}}
0312     //get the number of hadrons in the collection
0313     int num() {return hadrons.size();}
0314     //empty the collection
0315     void clear() {hadrons.clear();}
0316
0317     //overloading a few operators; ++, --, []
0318     //++ adds an empty particle, -- removes last particle, [i] allows access to i'th particle directly
0319 /*  hadron_collection& operator++(){HHhadron had; hadrons.push_back(had);}
0320     hadron_collection operator++(int){
0321         hadron_collection temp(*this);
0322         ++(*this);
0323         return temp;
0324     }
0325     hadron_collection& operator--(){hadrons.pop_back();}
0326     hadron_collection operator--(int){
0327         hadron_collection temp(*this);
0328         --(*this);
0329         return temp;
0330     }*/
0331     HHhadron& operator[](int i) {return hadrons[i];}
0332 //  const HHhadron& operator[](int i) {return hadrons[i];}
0333     const HHhadron& operator[](int i) const {return hadrons[i];}
0334   };
0335
0336   //used classes
0337   parton_collection HH_shower, HH_thermal;
0338   parton_collection HH_recomb_extrapartons;
0339   parton_collection HH_showerptns, HH_remnants, HH_pyremn;
0340   hadron_collection HH_hadrons, HH_pythia_hadrons;
0341
0342   //function to form strings out of original shower
0343   void stringform();
0344
0345   //recombination module
0346   void recomb();
0347
0348   //functions to set hadron id based on quark content, mass, and if it's in an excited state
0349   void set_baryon_id(parton_collection& qrks,HHhadron& had);
0350   void set_meson_id(parton_collection& qrks,HHhadron& had, int l, int k);
0351
0352   //gluon to q-qbar splitting function - for recombination use
0353   void gluon_decay(HHparton& glu, parton_collection& qrks);
0354
0355   //finding a sibling thermal parton for the "ithm'th" thermal parton
0356   int findthermalsibling(int ithm, parton_collection& therm);
0357   int findcloserepl(HHparton ptn, int iptn, bool lbt, bool thm, parton_collection& sh_lbt, parton_collection& therm);
0358   void findcloserepl_glu(HHparton ptn, int iptn, bool lbt, bool thm, parton_collection& sh_lbt, parton_collection& therm, int sel_out[]);
0359
0360   //function to prepare strings for input into Pythia8
0361   void stringprep(parton_collection& SP_remnants, parton_collection& SP_prepremn, bool cutstr);
0362
0363   //function to hand partons/strings and hadron resonances (and other color neutral objects) to Pythia8
0364   bool invoke_py();
0365
0366   // function to set the spacetime information for the hadrons coming from pythia
0367   void set_spacetime_for_pythia_hadrons(Pythia8::Event &event, int &size_input, std::vector<int> &eve_to_had, int pythia_attempt, bool find_positions, bool is_recohadron, bool recohadron_shsh);
0368
0369   // function to 'shake' the event a bit to bring the hadrons to their mass shell
0370   void bring_hadrons_to_mass_shell(hadron_collection& HH_hadrons);
0371
0372   void set_initial_parton_masses(parton_collection& HH_showerptns);
0373
0374   void convert_color_tags_to_int_type(vector<vector<shared_ptr<Parton>>>& shower);
0375
0376   // scale the energies/momenta of the negative hadrons to have energy conservation
0377   void scale_kinematics_negative_hadrons(hadron_collection& HH_hadrons, double shower_energy, double positive_hadrons_energy);
0378
0379   protected:
0380     static Pythia8::Pythia pythia;
0381
0382 };
0383
0384
0385 #endif // HYBRIDHADRONIZATION_H