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