Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HybridHadronization.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HybridHadronization.h
1 /*******************************************************************************
2  * Copyright (c) The JETSCAPE Collaboration, 2019
3  *
4  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
5  *
6  * For the list of contributors see AUTHORS.
7  *
8  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
9  *
10  * or via email to bugs.jetscape@gmail.com
11  *
12  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
13  * See COPYING for details.
14  ******************************************************************************/
15 
16 #ifndef HYBRIDHADRONIZATION_H
17 #define HYBRIDHADRONIZATION_H
18 
19 #include "HadronizationModule.h"
20 #include "JetScapeLogger.h"
21 #include "Pythia8/Pythia.h"
22 
23 #include <cmath>
24 #include <random>
25 #include <vector>
26 
27 using namespace Jetscape;
28 
29 class HybridHadronization : public HadronizationModule<HybridHadronization>
30 {
31  public:
32 
33  HybridHadronization();
34  virtual ~HybridHadronization();
35 
36  void Init();
37  void DoHadronization(vector<vector<shared_ptr<Parton>>>& shower, vector<shared_ptr<Hadron>>& hOut, vector<shared_ptr<Parton>>& pOut);
38  void WriteTask(weak_ptr<JetScapeWriter> w);
39 
40  private:
41  // Allows the registration of the module so that it is available to be used by the Jetscape framework.
42  static RegisterJetScapeModule<HybridHadronization> reg;
43 
44  double SigM2_calc(double R2chg, double qm1, double qm2, double qq1, double qq2);
45  double SigBR2_calc(double R2chg, double qm1, double qm2, double qm3, double qq1, double qq2, double qq3);
46  double SigBL2_calc(double SigBR2, double qm1, double qm2, double qm3);
47 
48  //double sigma_pi, sigma_k, sigma_nuc, maxE_level, gmax, xmq, xms, hbarc, dist2cut, sh_recofactor, th_recofactor, SigRB, SigLB;
49  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;
50  int number_p_fake;
51  double SigNucR2,SigNucL2,SigOmgR2,SigOmgL2,SigXiR2,SigXiL2,SigSigR2,SigSigL2,
52  SigOcccR2,SigOcccL2,SigOccR2,SigOccL2,SigXiccR2,SigXiccL2,SigOcR2,SigOcL2,SigXicR2,SigXicL2,SigSigcR2,SigSigcL2,
53  SigObbbR2,SigObbbL2,SigObbcR2,SigObbcL2,SigObbR2,SigObbL2,SigXibbR2,SigXibbL2,
54  SigObccR2,SigObccL2,SigObcR2,SigObcL2,SigXibcR2,SigXibcL2,SigObR2,SigObL2,SigXibR2,SigXibL2,SigSigbR2,SigSigbL2;
55  double SigPi2, SigPhi2, SigK2, SigJpi2, SigDs2, SigD2, SigUps2, SigBc2, SigB2;
56  bool inbrick, inhydro; int nreusehydro; double brickL, brickT;
57  const double pi = 3.1415926535897932384626433832795;
58  int attempts_max;
59  unsigned int rand_seed;
60  int reco_hadrons_pythia;
61  int additional_pythia_particles;
62  bool goldstonereco;
63  bool torder_reco;
64  bool afterburner_frag_hadrons = false;
65  std::string pythia_decays;
66 
67  //variables for recombination color structure
68  vector<vector<vector<int>>> Tempjunctions; // vector of all tempjunctions
69  vector<vector<int>> JunctionInfo; // vector of one junction's color tag and particle info
70  vector<int> IdColInfo1; vector<int> IdColInfo2; vector<int> IdColInfo3; vector<int> IdColInfo4;
71 
72  std::mt19937_64 eng; //RNG - Mersenne Twist - 64 bit
73  double ran();
74 
75  //4-vector boost
76  static FourVector HHboost(FourVector B, FourVector vec_in){
77  double xlam[4][4], beta2;
78  beta2 = B.x()*B.x() + B.y()*B.y() + B.z()*B.z();
79  B.Set(B.x(),B.y(),B.z(),1./(sqrt(1. - beta2)));
80 
81  double beta2inv;
82  if(beta2 > 0.){beta2inv = 1./beta2;}
83  else{beta2inv = 0.;}
84 
85  xlam[0][0] = B.t();
86  xlam[0][1] = -B.t()*B.x();
87  xlam[0][2] = -B.t()* B.y();
88  xlam[0][3] = -B.t()*B.z();
89  xlam[1][0] = xlam[0][1];
90  xlam[1][1] = 1.+(B.t() - 1.)*(B.x()*B.x())*beta2inv;
91  xlam[1][2] = (B.t()-1.)*B.x()*B.y()*beta2inv;
92  xlam[1][3] = (B.t()-1.)*B.x()*B.z()*beta2inv;
93  xlam[2][0] = xlam[0][2];
94  xlam[2][1] = xlam[1][2];
95  xlam[2][2] = 1.+(B.t()-1.)*(B.y()*B.y())*beta2inv;
96  xlam[2][3] = (B.t()-1.)*B.y()*B.z()*beta2inv;
97  xlam[3][0] = xlam[0][3];
98  xlam[3][1] = xlam[1][3];
99  xlam[3][2] = xlam[2][3];
100  xlam[3][3] = 1.+(B.t()-1.)*(B.z()*B.z())*beta2inv;
101 
102  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];
103  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];
104  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];
105  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];
106  FourVector vec_out(x_out,y_out,z_out,t_out);
107 
108  return vec_out;
109  }
110 
111  //3-vec (4-vec w/3 components) diff^2 function
112  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());}
113 
114  //parton class
115  //class parton{
116  class HHparton : public Parton{
117  protected:
118  //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
119  bool is_shower_, is_thermal_, is_used_, is_decayedglu_, is_remnant_, used_reco_, used_str_, is_strendpt_, used_junction_, is_fakep_;
120  //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')
121  //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
122  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_;
123  //parton mass
124  double alt_mass_;
125 
126  public:
127  //default constructor
128  HHparton() : Parton::Parton(1,1,1,0.,0.,0.,0.) {
129  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;
130  alt_id_ = 0; orig_ = 0; par_ = -1; string_id_ = 0; pos_str_ = 0; endpt_id_ = 0; sibling_ = 0; alt_mass_ = 0.;
131  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;
132  set_color(0); set_anti_color(0); set_stat(0);
133  }
134 
135  //getter functions
136  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();}
137  double px() {return p_in().x();} double py() {return p_in().y();} double pz() {return p_in().z();} double e() {return p_in().t();}
138  bool is_shower() {return is_shower_;} bool is_thermal() {return is_thermal_;} bool is_used() {return is_used_;} bool is_decayedglu() {return is_decayedglu_;}
139  bool is_remnant() {return is_remnant_;} bool used_reco() {return used_reco_;} bool used_str() {return used_str_;} bool is_strendpt() {return is_strendpt_;}
140  bool used_junction() {return used_junction_;} bool is_fakeparton() {return is_fakep_;}
141  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_;}
142  int endpt_id() {return endpt_id_;} int sibling() {return sibling_;} int PY_par1() {return PY_par1_;} int PY_par2() {return PY_par2_;}
143  int PY_dau1() {return PY_dau1_;} int PY_dau2() {return PY_dau2_;} int PY_stat() {return PY_stat_;} int PY_origid() {return PY_origid_;}
144  int PY_tag1() {return PY_tag1_;} int PY_tag2() {return PY_tag2_;} int PY_tag3() {return PY_tag3_;}
145  double mass() {return alt_mass_;}
146  FourVector pos() {return x_in();}
147  FourVector P() {return p_in();}
148  int status() {return pstat();} int col() {return color();} int acol() {return anti_color();}
149 
150  //setter functions
151  void x(double val) {x_in_.Set(val,x_in().y(),x_in().z(),x_in().t());}
152  void y(double val) {x_in_.Set(x_in().x(),val,x_in().z(),x_in().t());}
153  void z(double val) {x_in_.Set(x_in().x(),x_in().y(),val,x_in().t());}
154  void x_t(double val){x_in_.Set(x_in().x(),x_in().y(),x_in().z(),val);}
155  void pos(FourVector val) {x_in_.Set(val.x(),val.y(),val.z(),val.t());}
156  void px(double val) {reset_momentum(val,py(),pz(),e());}
157  void py(double val) {reset_momentum(px(),val,pz(),e());}
158  void pz(double val) {reset_momentum(px(),py(),val,e());}
159  void e(double val) {reset_momentum(px(),py(),pz(),val);}
160  void P(FourVector val) {reset_momentum(val);}
161 
162  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;}
163  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;}
164  void used_junction(bool val) {used_junction_ = val;} void is_fakeparton(bool val) {is_fakep_ = val;}
165  void id(int val) {alt_id_ = val;} void orig(int val) {orig_ = val;} void par(int val) {par_ = val;}
166  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;}
167  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;}
168  void PY_stat(int val) {PY_stat_ = val;} void PY_origid(int val) {PY_origid_ = val;}
169  void PY_tag1(int val) {PY_tag1_ = val;} void PY_tag2(int val) {PY_tag2_ = val;} void PY_tag3(int val) {PY_tag3_ = val;}
170  void mass(double val) {alt_mass_ = val;}
171  void status(int val) {set_stat(val);} void col(int val) {set_color(val);} void acol(int val) {set_anti_color(val);}
172 
173  //boost functions
174  FourVector boost_P(FourVector B){return HHboost(B,p_in());}
175  FourVector boost_P(double vx, double vy, double vz){FourVector B(vx,vy,vz,0.); return HHboost(B,p_in());}
176  FourVector boost_pos(FourVector B){return HHboost(B,x_in());}
177  FourVector boost_pos(double vx, double vy, double vz){FourVector B(vx,vy,vz,0.); return HHboost(B,x_in());}
178 
179  double pDif2(HHparton comp){return dif2(p_in(),comp.p_in());}
180  double posDif2(HHparton comp){return dif2(x_in(),comp.x_in());}
181 
182  };
183 
184  //hadron class
185  class HHhadron : public Hadron{
186  protected:
187  bool is_excited_, is_shsh_, is_shth_, is_thth_, is_recohad_, is_strhad_, is_final_;
188  //id: particle id, ?orig: particle origin(sh-sh, sh-th, th-th, recombined, stringfragmented...),
189  //par_str: parent string number, parh: parent hadron (decayed), status: status of particle (final, used, decayed, etc...)
190  int orig_, parstr_, parh_;
191  //hadron mass
192  double alt_mass_;
193  //vector of partonic parents
194  //std::vector<int> parents;
195 
196  public:
197  //default constructor
198  HHhadron() : Hadron::Hadron(1,1,1,0.,0.,0.,0.) {
199  is_excited_ = false; is_shsh_ = false; is_shth_ = false; is_thth_ = false; is_recohad_ = false; is_strhad_ = false; is_final_ = false;
200  orig_ = 0; parstr_ = 0; parh_ = -1; alt_mass_ = 0.; set_stat(0);
201  }
202 
203  //vector of partonic parents
204  std::vector<int> parents;
205 
206  //getter/setter for parents
207  int par(int i){if((i>=0) && (i<parents.size())){return parents[i];}else{return 999999;}}
208  void add_par(int i){parents.push_back(i);}
209 
210  //vector of colors (from partons that formed it)
211  std::vector<int> cols;
212 
213  //getter/setter for colors
214  int col(int i){if((i>=0) && (i<cols.size())){return cols[i];}else{return -1;}}
215  void add_col(int i){cols.push_back(i);}
216 
217  //getter functions
218  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();}
219  double px() {return p_in().x();} double py() {return p_in().y();} double pz() {return p_in().z();} double e() {return p_in().t();}
220  bool is_excited() {return is_excited_;} bool is_shsh() {return is_shsh_;} bool is_shth() {return is_shth_;} bool is_thth() {return is_thth_;}
221  bool is_recohad() {return is_recohad_;} bool is_strhad() {return is_strhad_;} bool is_final() {return is_final_;}
222  int orig() {return orig_;} int parstr() {return parstr_;} int parh() {return parh_;}
223  double mass() {return alt_mass_;}
224  FourVector pos() {return x_in();}
225  FourVector P() {return p_in();}
226  int id() {return pid();} int status() {return pstat();}
227 
228  //setter functions
229  void x(double val) {x_in_.Set(val,x_in().y(),x_in().z(),x_in().t());}
230  void y(double val) {x_in_.Set(x_in().x(),val,x_in().z(),x_in().t());}
231  void z(double val) {x_in_.Set(x_in().x(),x_in().y(),val,x_in().t());}
232  void x_t(double val){x_in_.Set(x_in().x(),x_in().y(),x_in().z(),val);}
233  void pos(FourVector val) {x_in_.Set(val.x(),val.y(),val.z(),val.t());}
234  void px(double val) {reset_momentum(val,py(),pz(),e());}
235  void py(double val) {reset_momentum(px(),val,pz(),e());}
236  void pz(double val) {reset_momentum(px(),py(),val,e());}
237  void e(double val) {reset_momentum(px(),py(),pz(),val);}
238  void P(FourVector val) {reset_momentum(val);}
239 
240  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;}
241  void is_recohad(bool val) {is_recohad_ = val;} void is_strhad(bool val) {is_strhad_ = val;} void is_final(bool val) {is_final_ = val;}
242  void orig(int val) {orig_ = val;} void parstr(int val) {parstr_ = val;} void parh(int val) {parh_ = val;}
243  void mass(double val) {alt_mass_ = val;}
244  void id(int val) {set_id(val);} void status(int val) {set_stat(val);}
245 
246  //Lorentz boost functions
247  FourVector boost_P(FourVector B){return HHboost(B,p_in());}
248  FourVector boost_P(double vx, double vy, double vz){FourVector B(vx,vy,vz,0.); return HHboost(B,p_in());}
249  FourVector boost_pos(FourVector B){return HHboost(B,x_in());}
250  FourVector boost_pos(double vx, double vy, double vz){FourVector B(vx,vy,vz,0.); return HHboost(B,x_in());}
251  };
252 
253  //class holding a collection of partons
254  class parton_collection{
255  public:
256  //the actual collection of partons
257  std::vector<HHparton> partons;
258  //default constructor
259  parton_collection() {partons.clear();}
260  //constructor given a single initial parton
261  parton_collection(HHparton par) {partons.push_back(par);}
262  //add a parton to the collection
263  void add(HHparton par) {partons.push_back(par);}
264  //add a collection of partons to the collection
265  void add(parton_collection pars) {for (int i=0; i<pars.num(); ++i){partons.push_back(pars[i]);}}
266  //get the number of partons in the collection
267  int num() {return partons.size();}
268  //empty the collection
269  void clear() {partons.clear();}
270  //insert a parton at position i
271  void insert(int i, HHparton newparton) {partons.insert(partons.begin()+i, newparton);}
272  //remove the parton at position i
273  void remove(int i) {partons.erase(partons.begin()+i);}
274  //swap partons at positions i, j
275  void swap(int i, int j) {if((i >= 0) && (i < partons.size()) && (j >= 0) && (j < partons.size())) {std::swap(partons[i],partons[j]);}}
276  //random access iterator pointing to first element in partons
277  std::vector<HHparton>::iterator begin() {return partons.begin();}
278  //random access iterator pointing to last element in partons
279  std::vector<HHparton>::iterator end() {return partons.end();}
280 
281  //overloading a few operators; ++, --, []
282  //++ adds an empty particle, -- removes last particle, [i] allows access to i'th particle directly
283 /* parton_collection& operator++(){HHparton par; partons.push_back(par);}
284  parton_collection operator++(int){
285  parton_collection temp(*this);
286  ++(*this);
287  return temp;
288  }
289  parton_collection& operator--(){partons.pop_back();}
290  parton_collection operator--(int){
291  parton_collection temp(*this);
292  --(*this);
293  return temp;
294  }*/
295  HHparton& operator[](int i) {return partons[i];}
296  const HHparton& operator[](int i) const {return partons[i];}
297  };
298 
299  //class holding a collection of hadrons
300  class hadron_collection{
301  public:
302  //the actual collection of hadrons
303  std::vector<HHhadron> hadrons;
304  //default constructor
305  hadron_collection() {hadrons.clear();}
306  //constructor given an initial hadron
307  hadron_collection(HHhadron had) {hadrons.push_back(had);}
308  //add a hadron to the collection
309  void add(HHhadron had) {hadrons.push_back(had);}
310  //add a collection of partons to the collection
311  void add(hadron_collection hads) {for (int i=0; i<hads.num(); ++i){hadrons.push_back(hads[i]);}}
312  //get the number of hadrons in the collection
313  int num() {return hadrons.size();}
314  //empty the collection
315  void clear() {hadrons.clear();}
316 
317  //overloading a few operators; ++, --, []
318  //++ adds an empty particle, -- removes last particle, [i] allows access to i'th particle directly
319 /* hadron_collection& operator++(){HHhadron had; hadrons.push_back(had);}
320  hadron_collection operator++(int){
321  hadron_collection temp(*this);
322  ++(*this);
323  return temp;
324  }
325  hadron_collection& operator--(){hadrons.pop_back();}
326  hadron_collection operator--(int){
327  hadron_collection temp(*this);
328  --(*this);
329  return temp;
330  }*/
331  HHhadron& operator[](int i) {return hadrons[i];}
332 // const HHhadron& operator[](int i) {return hadrons[i];}
333  const HHhadron& operator[](int i) const {return hadrons[i];}
334  };
335 
336  //used classes
337  parton_collection HH_shower, HH_thermal;
338  parton_collection HH_recomb_extrapartons;
339  parton_collection HH_showerptns, HH_remnants, HH_pyremn;
340  hadron_collection HH_hadrons, HH_pythia_hadrons;
341 
342  //function to form strings out of original shower
343  void stringform();
344 
345  //recombination module
346  void recomb();
347 
348  //functions to set hadron id based on quark content, mass, and if it's in an excited state
349  void set_baryon_id(parton_collection& qrks,HHhadron& had);
350  void set_meson_id(parton_collection& qrks,HHhadron& had, int l, int k);
351 
352  //gluon to q-qbar splitting function - for recombination use
353  void gluon_decay(HHparton& glu, parton_collection& qrks);
354 
355  //finding a sibling thermal parton for the "ithm'th" thermal parton
356  int findthermalsibling(int ithm, parton_collection& therm);
357  int findcloserepl(HHparton ptn, int iptn, bool lbt, bool thm, parton_collection& sh_lbt, parton_collection& therm);
358  void findcloserepl_glu(HHparton ptn, int iptn, bool lbt, bool thm, parton_collection& sh_lbt, parton_collection& therm, int sel_out[]);
359 
360  //function to prepare strings for input into Pythia8
361  void stringprep(parton_collection& SP_remnants, parton_collection& SP_prepremn, bool cutstr);
362 
363  //function to hand partons/strings and hadron resonances (and other color neutral objects) to Pythia8
364  bool invoke_py();
365 
366  // function to set the spacetime information for the hadrons coming from pythia
367  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);
368 
369  // function to 'shake' the event a bit to bring the hadrons to their mass shell
370  void bring_hadrons_to_mass_shell(hadron_collection& HH_hadrons);
371 
372  void set_initial_parton_masses(parton_collection& HH_showerptns);
373 
374  void convert_color_tags_to_int_type(vector<vector<shared_ptr<Parton>>>& shower);
375 
376  // scale the energies/momenta of the negative hadrons to have energy conservation
377  void scale_kinematics_negative_hadrons(hadron_collection& HH_hadrons, double shower_energy, double positive_hadrons_energy);
378 
379  protected:
380  static Pythia8::Pythia pythia;
381 
382 };
383 
384 
385 #endif // HYBRIDHADRONIZATION_H