GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Pythia8Hadro2019.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 
9  Shivesh Mandalia <s.p.mandalia@qmul.ac.uk>
10  Queen Mary University of London
11 */
12 //____________________________________________________________________________
13 
14 #include <RVersion.h>
15 #include <TClonesArray.h>
16 // Avoid the inclusion of dlfcn.h by Pythia.h that CINT is not able to process
17 #ifdef __CINT__
18 #define _DLFCN_H_
19 #define _DLFCN_H
20 #endif
21 
35 
36 #ifdef __GENIE_PYTHIA8_ENABLED__
37 #include "Pythia8/Pythia.h"
38 #endif
39 
40 using namespace genie;
41 using namespace genie::constants;
42 
43 //____________________________________________________________________________
45 PythiaBaseHadro2019("genie::Pythia8Hadro2019")
46 {
47  this->Initialize();
48 }
49 //____________________________________________________________________________
51 PythiaBaseHadro2019("genie::Pythia8Hadro2019", config)
52 {
53  this->Initialize();
54 }
55 //____________________________________________________________________________
57 {
58 #ifdef __GENIE_PYTHIA8_ENABLED__
59  delete fPythia;
60 #endif
61 }
62 //____________________________________________________________________________
64 #ifdef __GENIE_PYTHIA8_ENABLED__
65  event // avoid unused variable warning if PYTHIA6 is not enabled
66 #endif
67 ) const
68 {
69 #ifdef __GENIE_PYTHIA8_ENABLED__
71 #else
72  LOG("Pythia8Had", pFATAL)
73  << "Calling GENIE/PYTHIA8 hadronization modules without enabling PYTHIA8";
74  gAbortingInErr = true;
75  std::exit(1);
76 #endif
77 }
78 //____________________________________________________________________________
80 #ifdef __GENIE_PYTHIA8_ENABLED__
81  event // avoid unused variable warning if PYTHIA6 is not enabled
82 #endif
83 ) const
84 {
85 #ifdef __GENIE_PYTHIA8_ENABLED__
86  LOG("Pythia8Had", pNOTICE) << "Running PYTHIA8 hadronizer";
87 
88  const Interaction * interaction = event->Summary();
89  const Kinematics & kinematics = interaction->Kine();
90  double W = kinematics.W();
91 
92  LOG("Pythia8Had", pNOTICE)
93  << "Fragmentation: "
94  << "q = " << fLeadingQuark << ", qq = " << fRemnantDiquark
95  << ", W = " << W << " GeV";
96 
97  // Hadronize
98 
99  LOG("Pythia8Had", pDEBUG) << "Reseting PYTHIA8 event";
100  fPythia->event.reset();
101 
102  // Get quark/diquark masses
103  double mA = fPythia->particleData.m0(fLeadingQuark);
104  double mB = fPythia->particleData.m0(fRemnantDiquark);
105 
106  LOG("Pythia8Had", pINFO)
107  << "Leading quark mass = " << mA
108  << " GeV, remnant diqurak mass = " << mB << ", GeV";
109 
110  // Calculate quark/diquark energy/momentum
111  double pzAcm = 0.5 * Pythia8::sqrtpos(
112  (W + mA + mB) * (W - mA - mB) * (W - mA + mB) * (W + mA - mB) ) / W;
113  double pzBcm = -pzAcm;
114  double eA = sqrt(mA*mA + pzAcm*pzAcm);
115  double eB = sqrt(mB*mB + pzBcm*pzBcm);
116 
117  LOG("Pythia8Had", pINFO)
118  << "Quark: (pz = " << pzAcm << ", E = " << eA << ") GeV, "
119  << "Diquark: (pz = " << pzBcm << ", E = " << eB << ") GeV";
120 
121  // Pythia8 status code for outgoing particles of the hardest subprocesses is 23
122  // anti/colour tags for these 2 particles must complement each other
123  LOG("Pythia8Had", pDEBUG) << "Appending quark/diquark into the PYTHIA8 event";
124  fPythia->event.append(fLeadingQuark, 23, 101, 0, 0., 0., pzAcm, eA, mA);
125  fPythia->event.append(fRemnantDiquark, 23, 0, 101, 0., 0., pzBcm, eB, mB);
126  fPythia->event.list();
127 
128  LOG("Pythia8Had", pDEBUG) << "Generating next PYTHIA8 event";
129  fPythia->next();
130 
131  // List the event information
132  fPythia->event.list();
133  fPythia->stat();
134 
135  // Get LUJETS record
136  LOG("Pythia8Had", pDEBUG) << "Copying PYTHIA8 event record into GENIE's";
137  Pythia8::Event &fEvent = fPythia->event;
138  int np = fEvent.size();
139  assert(np>0);
140 
141  TClonesArray * particle_list = new TClonesArray("genie::GHepParticle", np);
142  particle_list->SetOwner(true);
143 
144  // Hadronic 4vec
145  TLorentzVector p4Had = kinematics.HadSystP4();
146 
147  // Vector defining rotation from LAB to LAB' (z:= \vec{phad})
148  TVector3 unitvq = p4Had.Vect().Unit();
149 
150  // Boost velocity LAB' -> HCM
151  TVector3 beta(0,0,p4Had.P()/p4Had.Energy());
152 
153  // Check target and decide appropriate status code for f/s particles
154  bool is_nucleus = interaction->InitState().Tgt().IsNucleus();
156 
157  // Get the index of the mother of the hadronic system
158  int mom = event->FinalStateHadronicSystemPosition();
159  assert(mom!=-1);
160 
161  // Get the neutrino vertex position (all hadrons positions set to this point)
162  GHepParticle * neutrino = event->Probe();
163  const TLorentzVector & vtx = *(neutrino->X4());
164 
165  // Loop over PYTHIA8 event particles and copy relevant entries
166  for (int i = 0; i < np; i++) {
167 
168  if (fEvent[i].id() == 90) continue; // ignore (system) pseudoparticle
169 
170  // Get PYTHIA8 particle ID and status codes
171  int particle_pdg_code = fEvent[i].id();
172  int pythia_particle_status = fEvent[i].status();
173 
174  // Sanity check
175  if(pythia_particle_status > 0) {
176  if( pdg::IsQuark (particle_pdg_code) ||
177  pdg::IsDiQuark(particle_pdg_code) )
178  {
179  LOG("Pythia8Had", pERROR)
180  << "Hadronization failed! Bare quarks appear in final state!";
181  return false;
182  }
183  }
184 
185  // Copy the initial q, qq (status = -23) and all undecayed particles
186  // Ignore particles we asked PYTHIA to decay (eg omega, Delta) but
187  // record their decay products
188  bool copy = (pythia_particle_status==-23) || (pythia_particle_status > 0);
189  if(copy) {
190  // The fragmentation products are generated in the hadronic CM frame
191  // where the z>0 axis is the \vec{phad} direction. For each particle
192  // returned by the hadronizer:
193  // - boost it back to LAB' frame {z:=\vec{phad}} / doesn't affect pT
194  // - rotate its 3-momentum from LAB' to LAB
195  TLorentzVector p4o(
196  fEvent[i].px(), fEvent[i].py(), fEvent[i].pz(), fEvent[i].e());
197  p4o.Boost(beta);
198  TVector3 p3 = p4o.Vect();
199  p3.RotateUz(unitvq);
200  TLorentzVector p4(p3,p4o.Energy());
201 
202  // Set the proper GENIE status according to a number of things:
203  // interaction on a nucleus or nucleon, particle type
204  GHepStatus_t ist = (pythia_particle_status > 0) ?
206  // Handle gammas, and leptons that might come from internal pythia decays
207  // mark them as final state particles
208  bool is_gamma = (particle_pdg_code == kPdgGamma);
209  bool is_nu = pdg::IsNeutralLepton(particle_pdg_code);
210  bool is_lchg = pdg::IsChargedLepton(particle_pdg_code);
211  bool not_hadr = is_gamma || is_nu || is_lchg;
212  if(not_hadr) { ist = kIStStableFinalState; }
213 
214  // Set mother/daugher indices
215  // int mother1 = mom+ fEvent[i].mother1();
216  // int mother2 = (pythia_particle_status > 0) ? mom + fEvent[i].mother2() : -1;
217  int mother1 = mom; // fEvent[i].mother1();
218  int mother2 = -1; //(pythia_particle_status > 0) ? mom + fEvent[i].mother2() : -1;
219  if(pythia_particle_status > 0) {
220  mother1 = mom+1;
221  mother2 = mom+2;
222  }
223  int daughter1 = -1;//(fEvent[i].daughter1() <= 0 ) ? -1 : mom + fEvent[i].daughter1();
224  int daughter2 = -1;//(fEvent[i].daughter1() <= 0 ) ? -1 : mom + fEvent[i].daughter2();
225 
226  // Create GHepParticle
227  GHepParticle particle = GHepParticle(
228  particle_pdg_code, // pdg
229  ist, // status
230  mother1, // first parent
231  mother2, // second parent
232  daughter1, // first daughter
233  daughter2, // second daughter
234  p4.Px(), // px
235  p4.Py(), // py
236  p4.Pz(), // pz
237  p4.Energy(), // e
238  vtx.X(), // x
239  vtx.Y(), // y
240  vtx.Z(), // z
241  vtx.T() // t
242  );
243 
244  LOG("Pythia8Had", pDEBUG)
245  << "Adding final state particle pdgc = " << particle.Pdg()
246  << " with status = " << particle.Status();
247 
248  // Insert the particle in the list
249  event->AddParticle(particle);
250  }// copy?
251  }// loop over particles
252 
253  return true;
254 
255 #else
256  return false;
257 #endif
258 }
259 //____________________________________________________________________________
261 {
262 #ifdef __GENIE_PYTHIA8_ENABLED__
263  fOriDecayFlag_pi0 = fPythia->particleData.canDecay(kPdgPi0);
264  fOriDecayFlag_K0 = fPythia->particleData.canDecay(kPdgK0);
265  fOriDecayFlag_K0b = fPythia->particleData.canDecay(kPdgAntiK0);
266  fOriDecayFlag_L0 = fPythia->particleData.canDecay(kPdgLambda);
267  fOriDecayFlag_L0b = fPythia->particleData.canDecay(kPdgAntiLambda);
268  fOriDecayFlag_Dm = fPythia->particleData.canDecay(kPdgP33m1232_DeltaM);
269  fOriDecayFlag_D0 = fPythia->particleData.canDecay(kPdgP33m1232_Delta0);
270  fOriDecayFlag_Dp = fPythia->particleData.canDecay(kPdgP33m1232_DeltaP);
271  fOriDecayFlag_Dpp = fPythia->particleData.canDecay(kPdgP33m1232_DeltaPP);
272 
273 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
274  LOG("Pythia8Had", pDEBUG)
275  << "Original PYTHIA6 decay flags:"
276  << "\n pi0 = " << fOriDecayFlag_pi0
277  << "\n K0 = " << fOriDecayFlag_K0
278  << "\n \bar{K0} = " << fOriDecayFlag_K0b
279  << "\n Lambda = " << fOriDecayFlag_L0
280  << "\n \bar{Lambda0} = " << fOriDecayFlag_L0b
281  << "\n D- = " << fOriDecayFlag_Dm
282  << "\n D0 = " << fOriDecayFlag_D0
283  << "\n D+ = " << fOriDecayFlag_Dp
284  << "\n D++ = " << fOriDecayFlag_Dpp;
285 #endif
286 
287 #endif
288 }
289 //____________________________________________________________________________
291 {
292 #ifdef __GENIE_PYTHIA8_ENABLED__
293  fPythia->particleData.mayDecay(kPdgPi0, fReqDecayFlag_pi0 );
294  fPythia->particleData.mayDecay(kPdgK0, fReqDecayFlag_K0 );
295  fPythia->particleData.mayDecay(kPdgAntiK0, fReqDecayFlag_K0b );
296  fPythia->particleData.mayDecay(kPdgLambda, fReqDecayFlag_L0 );
297  fPythia->particleData.mayDecay(kPdgAntiLambda, fReqDecayFlag_L0b );
298  fPythia->particleData.mayDecay(kPdgP33m1232_DeltaM, fReqDecayFlag_Dm );
299  fPythia->particleData.mayDecay(kPdgP33m1232_Delta0, fReqDecayFlag_D0 );
300  fPythia->particleData.mayDecay(kPdgP33m1232_DeltaP, fReqDecayFlag_Dp );
301  fPythia->particleData.mayDecay(kPdgP33m1232_DeltaPP, fReqDecayFlag_Dpp );
302 #endif
303 }
304 //____________________________________________________________________________
306 {
307 #ifdef __GENIE_PYTHIA8_ENABLED__
308  fPythia->particleData.mayDecay(kPdgPi0, fOriDecayFlag_pi0 );
309  fPythia->particleData.mayDecay(kPdgK0, fOriDecayFlag_K0 );
310  fPythia->particleData.mayDecay(kPdgAntiK0, fOriDecayFlag_K0b );
311  fPythia->particleData.mayDecay(kPdgLambda, fOriDecayFlag_L0 );
312  fPythia->particleData.mayDecay(kPdgAntiLambda, fOriDecayFlag_L0b );
313  fPythia->particleData.mayDecay(kPdgP33m1232_DeltaM, fOriDecayFlag_Dm );
314  fPythia->particleData.mayDecay(kPdgP33m1232_Delta0, fOriDecayFlag_D0 );
315  fPythia->particleData.mayDecay(kPdgP33m1232_DeltaP, fOriDecayFlag_Dp );
316  fPythia->particleData.mayDecay(kPdgP33m1232_DeltaPP, fOriDecayFlag_Dpp );
317 #endif
318 }
319 //____________________________________________________________________________
321 {
322  Algorithm::Configure(config);
323  this->LoadConfig();
324 }
325 //____________________________________________________________________________
326 void Pythia8Hadro2019::Configure(string config)
327 {
328  Algorithm::Configure(config);
329  this->LoadConfig();
330 }
331 //____________________________________________________________________________
333 {
335 
336 #ifdef __GENIE_PYTHIA8_ENABLED__
337  fPythia->settings.parm("StringFlav:probStoUD", fSSBarSuppression);
338  fPythia->settings.parm("Diffraction:primKTwidth", fGaussianPt2);
339  fPythia->settings.parm("StringPT:enhancedFraction", fNonGaussianPt2Tail);
340  fPythia->settings.parm("StringFragmentation:stopMass", fRemainingECutoff);
341  fPythia->settings.parm("StringFlav:probQQtoQ", fDiQuarkSuppression);
342  fPythia->settings.parm("StringFlav:mesonUDvector", fLightVMesonSuppression);
343  fPythia->settings.parm("StringFlav:mesonSvector", fSVMesonSuppression);
344  fPythia->settings.parm("StringZ:aLund", fLunda);
345  fPythia->settings.parm("StringZ:bLund", fLundb);
346  fPythia->settings.parm("StringZ:aExtraDiquark", fLundaDiq);
347 
348  fPythia->init(); // needed again to read the above?
349 
350 #endif
351 
352  LOG("Pythia8Had", pDEBUG) << this->GetConfig();
353 }
354 //____________________________________________________________________________
356 {
357 #ifdef __GENIE_PYTHIA8_ENABLED__
358  fPythia = new Pythia8::Pythia();
359 
360  fPythia->readString("ProcessLevel:all = off");
361  fPythia->readString("Print:quiet = on");
362 
363  // sync GENIE and PYTHIA8 seeds
364  RandomGen * rnd = RandomGen::Instance();
365  long int seed = rnd->GetSeed();
366  fPythia->readString("Random:setSeed = on");
367  fPythia->settings.mode("Random:seed", seed);
368  LOG("Pythia8Had", pINFO)
369  << "PYTHIA8 seed = " << fPythia->settings.mode("Random:seed");
370 
371  fPythia->init();
372 
373 #endif
374 }
375 //____________________________________________________________________________
double fSSBarSuppression
ssbar suppression
const int kPdgP33m1232_DeltaPP
Definition: PDGCodes.h:107
double W(bool selected=false) const
Definition: Kinematics.cxx:157
void ProcessEventRecord(GHepRecord *event) const
#define pERROR
Definition: Messenger.h:59
const int kPdgLambda
Definition: PDGCodes.h:85
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
#define pFATAL
Definition: Messenger.h:56
bool IsNucleus(void) const
Definition: Target.cxx:272
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
virtual void ProcessEventRecord(GHepRecord *event) const
double fGaussianPt2
gaussian pt2 distribution width
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:101
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:66
double fDiQuarkSuppression
di-quark suppression parameter
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
GHepStatus_t Status(void) const
Definition: GHepParticle.h:64
void Configure(const Registry &config)
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:246
const int kPdgP33m1232_DeltaP
Definition: PDGCodes.h:106
double fSVMesonSuppression
strange vector meson suppression
const int kPdgP33m1232_DeltaM
Definition: PDGCodes.h:104
const int kPdgK0
Definition: PDGCodes.h:174
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
int Pdg(void) const
Definition: GHepParticle.h:63
bool Hadronize(GHepRecord *event) const
Summary information for an interaction.
Definition: Interaction.h:56
const double e
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
double fLightVMesonSuppression
light vector meson suppression
double fLundaDiq
adjustment of Lund a for di-quark
const int kPdgGamma
Definition: PDGCodes.h:189
long int GetSeed(void) const
Definition: RandomGen.h:82
bool IsDiQuark(int pdgc)
Definition: PDGUtils.cxx:231
const Kinematics & Kine(void) const
Definition: Interaction.h:71
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
double fNonGaussianPt2Tail
non gaussian pt2 tail parameterization
double fLunda
Lund a parameter.
const int kPdgPi0
Definition: PDGCodes.h:160
Base class for the Pythia (6 and 8) hadronization modules in GENIE. In particular, the base class provides common checks and basic assignments of quark/diquark codes for a no frills interface to Pythia hadronization routines.
#define pINFO
Definition: Messenger.h:62
const int kPdgAntiK0
Definition: PDGCodes.h:175
const int kPdgP33m1232_Delta0
Definition: PDGCodes.h:105
void CopyOriginalDecayFlags(void) const
bool IsNeutralLepton(int pdgc)
Definition: PDGUtils.cxx:95
void SetDesiredDecayFlags(void) const
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
bool IsQuark(int pdgc)
Definition: PDGUtils.cxx:250
const InitialState & InitState(void) const
Definition: Interaction.h:69
double fLundb
Lund b parameter.
#define pNOTICE
Definition: Messenger.h:61
const Target & Tgt(void) const
Definition: InitialState.h:66
const int kPdgAntiLambda
Definition: PDGCodes.h:86
bool gAbortingInErr
Definition: Messenger.cxx:34
void RestoreOriginalDecayFlags(void) const
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
double fRemainingECutoff
remaining E cutoff stopping fragmentation
enum genie::EGHepStatus GHepStatus_t
#define pDEBUG
Definition: Messenger.h:63