GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Pythia6Hadro2019.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 //____________________________________________________________________________
10 
11 #include <RVersion.h>
12 #include <TClonesArray.h>
13 
27 
28 #ifdef __GENIE_PYTHIA6_ENABLED__
29 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6)
30 #include <TMCParticle.h>
31 #else
32 #include <TMCParticle6.h>
33 #endif
34 #endif // __GENIE_PYTHIA6_ENABLED__
35 
36 using namespace genie;
37 using namespace genie::constants;
38 
39 // the actual PYTHIA call
40 extern "C" void py2ent_(int *, int *, int *, double *);
41 
42 //____________________________________________________________________________
44 PythiaBaseHadro2019("genie::Pythia6Hadro2019")
45 {
46  this->Initialize();
47 }
48 //____________________________________________________________________________
50 PythiaBaseHadro2019("genie::Pythia6Hadro2019", config)
51 {
52  this->Initialize();
53 }
54 //____________________________________________________________________________
56 {
57 
58 }
59 //____________________________________________________________________________
62  event // avoid unused variable warning if PYTHIA6 is not enabled
63  #endif
64 ) const
65 {
66 #ifdef __GENIE_PYTHIA6_ENABLED__
68 #else
69  LOG("Pythia6Had", pFATAL)
70  << "Calling GENIE/PYTHIA6 hadronization modules without enabling PYTHIA6";
71  gAbortingInErr = true;
72  std::exit(1);
73 #endif
74 }
75 //____________________________________________________________________________
78  event // avoid unused variable warning if PYTHIA6 is not enabled
79 #endif
80 ) const
81 {
82 #ifdef __GENIE_PYTHIA6_ENABLED__
83 
84  LOG("Pythia6Had", pNOTICE) << "Running PYTHIA6 hadronizer";
85 
86  const Interaction * interaction = event->Summary();
87  const Kinematics & kinematics = interaction->Kine();
88  double W = kinematics.W();
89 
90  LOG("Pythia6Had", pNOTICE)
91  << "Fragmentation: "
92  << "q = " << fLeadingQuark << ", qq = " << fRemnantDiquark
93  << ", W = " << W;
94 
95  // Hadronize
96  int ip = 0;
97  py2ent_(&ip, &fLeadingQuark, &fRemnantDiquark, &W); // hadronizer
98 
99  // Get LUJETS record
100  fPythia->GetPrimaries();
101  TClonesArray * pythia_particles =
102  (TClonesArray *) fPythia->ImportParticles("All");
103 
104  // Copy PYTHIA container to a new TClonesArray so as to transfer ownership
105  // of the container and of its elements to the calling method
106 
107  int np = pythia_particles->GetEntries();
108  assert(np>0);
109  TClonesArray * particle_list = new TClonesArray("genie::GHepParticle", np);
110  particle_list->SetOwner(true);
111 
112  // Hadronic 4vec
113  TLorentzVector p4Had = kinematics.HadSystP4();
114 
115  // Vector defining rotation from LAB to LAB' (z:= \vec{phad})
116  TVector3 unitvq = p4Had.Vect().Unit();
117 
118  // Boost velocity LAB' -> HCM
119  TVector3 beta(0,0,p4Had.P()/p4Had.Energy());
120 
121  // Check target and decide appropriate status code for f/s particles
122  bool is_nucleus = interaction->InitState().Tgt().IsNucleus();
124 
125  // Get the index of the mother of the hadronic system
126  int mom = event->FinalStateHadronicSystemPosition();
127  assert(mom!=-1);
128 
129  // Get the neutrino vertex position (all hadrons positions set to this point)
130  GHepParticle * neutrino = event->Probe();
131  const TLorentzVector & vtx = *(neutrino->X4());
132 
133  // Loop over PYTHIA8 event particles and copy relevant entries
134  unsigned int i = 0;
135  TMCParticle * p = 0;
136  TIter particle_iter(pythia_particles);
137  while( (p = (TMCParticle *) particle_iter.Next()) ) {
138 
139  int particle_pdg_code = p->GetKF();
140  int pythia_particle_status = p->GetKS();
141 
142  // Sanity check
143  if(pythia_particle_status == 1) {
144  if( pdg::IsQuark (particle_pdg_code) ||
145  pdg::IsDiQuark(particle_pdg_code) )
146  {
147  LOG("Pythia6Had", pERROR)
148  << "Hadronization failed! Bare quarks appear in final state!";
149  return false;
150  }
151  }
152 
153  // The fragmentation products are generated in the hadronic CM frame
154  // where the z>0 axis is the \vec{phad} direction. For each particle
155  // returned by the hadronizer:
156  // - boost it back to LAB' frame {z:=\vec{phad}} / doesn't affect pT
157  // - rotate its 3-momentum from LAB' to LAB
158  TLorentzVector p4o(p->GetPx(), p->GetPy(), p->GetPz(), p->GetEnergy());
159  p4o.Boost(beta);
160  TVector3 p3 = p4o.Vect();
161  p3.RotateUz(unitvq);
162  TLorentzVector p4(p3,p4o.Energy());
163 
164  // Set the proper GENIE status according to a number of things:
165  // interaction on a nucleus or nucleon, particle type
166  GHepStatus_t ist = (pythia_particle_status == 1) ?
168  // Handle gammas, and leptons that might come from internal pythia decays
169  // mark them as final state particles
170  bool is_gamma = (particle_pdg_code == kPdgGamma);
171  bool is_nu = pdg::IsNeutralLepton(particle_pdg_code);
172  bool is_lchg = pdg::IsChargedLepton(particle_pdg_code);
173  bool not_hadr = is_gamma || is_nu || is_lchg;
174  if(not_hadr) { ist = kIStStableFinalState; }
175 
176  // Set mother/daugher indices
177  int mother1 = mom + p->GetParent();
178  int mother2 = -1;
179  int daughter1 = (p->GetFirstChild() <= 0 ) ? -1 : mom + p->GetFirstChild();
180  int daughter2 = (p->GetLastChild() <= 0 ) ? -1 : mom + p->GetLastChild();
181 
182  // Create GHepParticle
183  GHepParticle particle = GHepParticle(
184  particle_pdg_code, // pdg
185  ist, // status
186  mother1, // first parent
187  mother2, // second parent
188  daughter1, // first daughter
189  daughter2, // second daughter
190  p4.Px(), // px
191  p4.Py(), // py
192  p4.Pz(), // pz
193  p4.Energy(), // e
194  vtx.X(), // x
195  vtx.Y(), // y
196  vtx.Z(), // z
197  vtx.T() // t
198  );
199 
200  LOG("Pythia6Had", pDEBUG)
201  << "Adding final state particle pdgc = " << particle.Pdg()
202  << " with status = " << particle.Status();
203 
204  // Insert the particle in the list
205  event->AddParticle(particle);
206  }
207  return true;
208 
209 #else
210  return false;
211 #endif // __GENIE_PYTHIA6_ENABLED__
212 }
213 //____________________________________________________________________________
215 {
216 #ifdef __GENIE_PYTHIA6_ENABLED__
217  fOriDecayFlag_pi0 = (fPythia->GetMDCY(fPythia->Pycomp(kPdgPi0), 1) == 1);
218  fOriDecayFlag_K0 = (fPythia->GetMDCY(fPythia->Pycomp(kPdgK0), 1) == 1);
219  fOriDecayFlag_K0b = (fPythia->GetMDCY(fPythia->Pycomp(kPdgAntiK0), 1) == 1);
220  fOriDecayFlag_L0 = (fPythia->GetMDCY(fPythia->Pycomp(kPdgLambda), 1) == 1);
221  fOriDecayFlag_L0b = (fPythia->GetMDCY(fPythia->Pycomp(kPdgAntiLambda), 1) == 1);
222  fOriDecayFlag_Dm = (fPythia->GetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaM), 1) == 1);
223  fOriDecayFlag_D0 = (fPythia->GetMDCY(fPythia->Pycomp(kPdgP33m1232_Delta0), 1) == 1);
224  fOriDecayFlag_Dp = (fPythia->GetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaP), 1) == 1);
225  fOriDecayFlag_Dpp = (fPythia->GetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaPP), 1) == 1);
226 
227 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
228  LOG("Pythia6Had", pDEBUG)
229  << "Original PYTHIA6 decay flags:"
230  << "\n pi0 = " << fOriDecayFlag_pi0
231  << "\n K0 = " << fOriDecayFlag_K0
232  << "\n \bar{K0} = " << fOriDecayFlag_K0b
233  << "\n Lambda = " << fOriDecayFlag_L0
234  << "\n \bar{Lambda0} = " << fOriDecayFlag_L0b
235  << "\n D- = " << fOriDecayFlag_Dm
236  << "\n D0 = " << fOriDecayFlag_D0
237  << "\n D+ = " << fOriDecayFlag_Dp
238  << "\n D++ = " << fOriDecayFlag_Dpp;
239 #endif
240 
241 #endif // __GENIE_PYTHIA6_ENABLED__
242 }
243 //____________________________________________________________________________
245 {
246 #ifdef __GENIE_PYTHIA6_ENABLED__
247 
248  fPythia->SetMDCY(fPythia->Pycomp(kPdgPi0),
249  1, (fReqDecayFlag_pi0) ? 1 : 0);
250  fPythia->SetMDCY(fPythia->Pycomp(kPdgK0),
251  1, (fReqDecayFlag_K0) ? 1 : 0);
252  fPythia->SetMDCY(fPythia->Pycomp(kPdgAntiK0),
253  1, (fReqDecayFlag_K0b) ? 1 : 0);
254  fPythia->SetMDCY(fPythia->Pycomp(kPdgLambda),
255  1, (fReqDecayFlag_L0) ? 1 : 0);
256  fPythia->SetMDCY(fPythia->Pycomp(kPdgAntiLambda),
257  1, (fReqDecayFlag_L0b) ? 1 : 0);
258  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaM),
259  1, (fReqDecayFlag_Dm) ? 1 : 0);
260  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_Delta0),
261  1, (fReqDecayFlag_D0) ? 1 : 0);
262  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaP),
263  1, (fReqDecayFlag_Dp) ? 1 : 0);
264  fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaPP),
265  1, (fReqDecayFlag_Dpp) ? 1 : 0);
266 
267 #endif // __GENIE_PYTHIA6_ENABLED__
268 }
269 //____________________________________________________________________________
271 {
272 #ifdef __GENIE_PYTHIA6_ENABLED__
273 
274 fPythia->SetMDCY(fPythia->Pycomp(kPdgPi0),
275  1, (fOriDecayFlag_pi0) ? 1 : 0);
276 fPythia->SetMDCY(fPythia->Pycomp(kPdgK0),
277  1, (fOriDecayFlag_K0) ? 1 : 0);
278 fPythia->SetMDCY(fPythia->Pycomp(kPdgAntiK0),
279  1, (fOriDecayFlag_K0b) ? 1 : 0);
280 fPythia->SetMDCY(fPythia->Pycomp(kPdgLambda),
281  1, (fOriDecayFlag_L0) ? 1 : 0);
282 fPythia->SetMDCY(fPythia->Pycomp(kPdgAntiLambda),
283  1, (fOriDecayFlag_L0b) ? 1 : 0);
284 fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaM),
285  1, (fOriDecayFlag_Dm) ? 1 : 0);
286 fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_Delta0),
287  1, (fOriDecayFlag_D0) ? 1 : 0);
288 fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaP),
289  1, (fOriDecayFlag_Dp) ? 1 : 0);
290 fPythia->SetMDCY(fPythia->Pycomp(kPdgP33m1232_DeltaPP),
291  1, (fOriDecayFlag_Dpp) ? 1 : 0);
292 
293 #endif // __GENIE_PYTHIA6_ENABLED__
294 }
295 //____________________________________________________________________________
297 {
298  Algorithm::Configure(config);
299  this->LoadConfig();
300 }
301 //____________________________________________________________________________
302 void Pythia6Hadro2019::Configure(string config)
303 {
304  Algorithm::Configure(config);
305  this->LoadConfig();
306 }
307 //____________________________________________________________________________
309 {
311 
312 #ifdef __GENIE_PYTHIA6_ENABLED__
313  fPythia->SetPARJ(2, fSSBarSuppression );
314  fPythia->SetPARJ(21, fGaussianPt2 );
315  fPythia->SetPARJ(23, fNonGaussianPt2Tail );
316  fPythia->SetPARJ(33, fRemainingECutoff );
317  fPythia->SetPARJ(1, fDiQuarkSuppression );
318  fPythia->SetPARJ(11, fLightVMesonSuppression );
319  fPythia->SetPARJ(12, fSVMesonSuppression );
320  fPythia->SetPARJ(41, fLunda );
321  fPythia->SetPARJ(42, fLundb );
322  fPythia->SetPARJ(45, fLundaDiq );
323 #endif
324 
325  LOG("Pythia6Had", pDEBUG) << this->GetConfig() ;
326 }
327 //____________________________________________________________________________
329 {
331 #ifdef __GENIE_PYTHIA6_ENABLED__
332  fPythia = TPythia6::Instance();
333  // sync GENIE/PYTHIA6 seed number
334  // PYTHIA6 is a singleton, so do this from RandomGen for all
335  // GENIE algorithms that use PYTHIA6
337 #endif
338 }
339 //____________________________________________________________________________
double fSSBarSuppression
ssbar suppression
const int kPdgP33m1232_DeltaPP
Definition: PDGCodes.h:107
void CopyOriginalDecayFlags(void) const
double W(bool selected=false) const
Definition: Kinematics.cxx:157
void ProcessEventRecord(GHepRecord *event) const
#define pERROR
Definition: Messenger.h:59
bool Hadronize(GHepRecord *event) const
const int kPdgLambda
Definition: PDGCodes.h:85
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
#define __GENIE_PYTHIA6_ENABLED__
#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
void RestoreOriginalDecayFlags(void) const
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:101
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:66
double fDiQuarkSuppression
di-quark suppression parameter
TPythia6 * fPythia
PYTHIA6 wrapper class.
GHepStatus_t Status(void) const
Definition: GHepParticle.h:64
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
Summary information for an interaction.
Definition: Interaction.h:56
#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
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
void py2ent_(int *, int *, int *, double *)
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.
void Configure(const Registry &config)
const int kPdgAntiK0
Definition: PDGCodes.h:175
const int kPdgP33m1232_Delta0
Definition: PDGCodes.h:105
bool IsNeutralLepton(int pdgc)
Definition: PDGUtils.cxx:95
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
void SetDesiredDecayFlags(void) const
const int kPdgAntiLambda
Definition: PDGCodes.h:86
bool gAbortingInErr
Definition: Messenger.cxx:34
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