GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PythiaDecayer.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 <vector>
12 #include <cassert>
13 
14 #include <TClonesArray.h>
15 #include <TLorentzVector.h>
16 #include <TDecayChannel.h>
17 #include <RVersion.h>
18 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6)
19 #include <TMCParticle.h>
20 #else
21 #include <TMCParticle6.h>
22 #endif
23 
36 
37 using std::vector;
38 
39 using namespace genie;
40 
41 // actual PYTHIA calls:
42 extern "C" void py1ent_(int *, int *, double *, double *, double *);
43 extern "C" void pydecy_(int *);
44 //____________________________________________________________________________
46 Decayer("genie::PythiaDecayer")
47 {
48  this->Initialize();
49 }
50 //____________________________________________________________________________
52 Decayer("genie::PythiaDecayer", config)
53 {
54  this->Initialize();
55 }
56 //____________________________________________________________________________
58 {
59 
60 }
61 //____________________________________________________________________________
63 {
64  LOG("ResonanceDecay", pINFO)
65  << "Running PYTHIA6 particle decayer "
66  << ((fRunBefHadroTransp) ? "*before*" : "*after*") << " FSI";
67 
68  // Loop over particles, find unstable ones and decay them
69  TObjArrayIter piter(event);
70  GHepParticle * p = 0;
71  int ipos = -1;
72 
73  while( (p = (GHepParticle *) piter.Next()) ) {
74  ipos++;
75 
76  LOG("Pythia6Decay", pDEBUG) << "Checking: " << p->Name();
77 
78  int pdg_code = p->Pdg();
79  GHepStatus_t status_code = p->Status();
80 
81  if(!this->IsHandled (pdg_code)) continue;
82  if(!this->ToBeDecayed(pdg_code, status_code)) continue;
83 
84  LOG("Pythia6Decay", pNOTICE)
85  << "Decaying unstable particle: " << p->Name();
86 
87  bool decayed = this->Decay(ipos, event);
88  assert(decayed); // handle this more graciously and throw an exception
89  }
90 
91  LOG("Pythia6Decay", pNOTICE)
92  << "Done finding & decaying unstable particles";
93 }
94 //____________________________________________________________________________
95 bool PythiaDecayer::Decay(int decay_particle_id, GHepRecord * event) const
96 {
97  fWeight = 1.; // reset previous decay weight
98 
99  // Get particle to be decayed
100  GHepParticle * decay_particle = event->Particle(decay_particle_id);
101  if(!decay_particle) return 0;
102 
103  // Get the particle 4-momentum, 4-position and PDG code
104  TLorentzVector decay_particle_p4 = *(decay_particle->P4());
105  int decay_particle_pdg_code = decay_particle->Pdg();
106 
107  // Convert to PYTHIA6 particle code and check whether decay is inhibited
108  int kc = fPythia->Pycomp(decay_particle_pdg_code);
109  int mdcy = fPythia->GetMDCY(kc, 1);
110  if(mdcy == 0) {
111  LOG("Pythia6Decay", pNOTICE)
112  << (PDGLibrary::Instance())->Find(decay_particle_pdg_code)->GetName()
113  << " decays are inhibited!";
114  return false;
115  }
116 
117  // Get sub of BRs and compute weight if decay channels were inhibited
118  double sumbr = this->SumOfBranchingRatios(kc);
119  if(sumbr <= 0) {
120  LOG("Pythia6Decay", pWARN)
121  << "The sum of enabled "
122  << (PDGLibrary::Instance())->Find(decay_particle_pdg_code)->GetName()
123  << " decay channel branching ratios is non-positive!";
124  return false;
125  }
126  fWeight = 1./sumbr; // update weight to account for inhibited channels
127 
128  // Run PYTHIA6 decay
129  int ip = 0;
130  double E = decay_particle_p4.Energy();
131  double theta = decay_particle_p4.Theta();
132  double phi = decay_particle_p4.Phi();
133  fPythia->SetMSTJ(22,1);
134  py1ent_(&ip, &decay_particle_pdg_code, &E, &theta, &phi);
135 
136  // Get decay products
137  fPythia->GetPrimaries();
138  TClonesArray * impl = (TClonesArray *) fPythia->ImportParticles("All");
139  if(!impl) {
140  LOG("Pythia6Decay", pWARN)
141  << "TPythia6::ImportParticles() returns a null list!";
142  return false;
143  }
144 
145  // Copy the PYTHIA6 container to the GENIE event record
146 
147  // Check whether the interaction is off a nuclear target or free nucleon
148  // Depending on whether this module is run before or after the hadron
149  // transport module it would affect the daughters status code
150  GHepParticle * target_nucleus = event->TargetNucleus();
151  bool in_nucleus = (target_nucleus!=0);
152 
153  // the values of space coordinates from pythia are in mm.
154  // our conventions want it in fm
155  constexpr double space_scale = units::mm / units::fm ;
156 
157  // the values of time coordinate from pythia is in mm/c.
158  // our conventions want it in ys
159  constexpr double time_scale = 1e21 * units::m / units::s ;
160 
161  TMCParticle * p = 0;
162  TIter particle_iter(impl);
163  while( (p = (TMCParticle *) particle_iter.Next()) ) {
164  // Convert from TMCParticle to GHepParticle
166  p->GetKF(), // pdg
167  GHepStatus_t(p->GetKS()), // status
168  p->GetParent(), // first parent
169  0, // second parent
170  p->GetFirstChild(), // first daughter
171  p->GetLastChild(), // second daughter
172  p->GetPx(), // px
173  p->GetPy(), // py
174  p->GetPz(), // pz
175  p->GetEnergy(), // e
176  p->GetVx() * space_scale , // x
177  p->GetVy() * space_scale , // y
178  p->GetVz() * space_scale , // z
179  p->GetTime() * time_scale // t
180  );
181 
182  if(mcp.Status()==kIStNucleonTarget) continue; // mother particle, already in GHEP
183 
184  int daughter_pdg_code = mcp.Pdg();
185  SLOG("Pythia6Decay", pINFO)
186  << "Adding daughter particle wit PDG code = "
187  << daughter_pdg_code << ", m = " << mcp.Mass()
188  << " GeV, E = " << mcp.Energy() << " GeV)";
189 
190  bool is_hadron = pdg::IsHadron(daughter_pdg_code);
191  bool hadron_in_nuc = (in_nucleus && is_hadron && fRunBefHadroTransp);
192 
193  GHepStatus_t daughter_status_code = (hadron_in_nuc) ?
195 
196  TLorentzVector daughter_p4(
197  mcp.Px(),mcp.Py(),mcp.Pz(),mcp.Energy());
198 
199  event->AddParticle(
200  daughter_pdg_code, daughter_status_code,
201  decay_particle_id,-1,-1,-1,
202  daughter_p4, * mcp.X4() );
203  }
204 
205  // Update the event weight for each weighted particle decay
206  double weight = event->Weight() * fWeight;
207  event->SetWeight(weight);
208 
209  // Mark input particle as a 'decayed state' & add its daughter links
210  decay_particle->SetStatus(kIStDecayedState);
211 
212  return true;
213 }
214 //____________________________________________________________________________
216 {
217  fPythia = TPythia6::Instance();
218  fWeight = 1.;
219 
220  // sync GENIE/PYTHIA6 seeds
222 }
223 //____________________________________________________________________________
224 bool PythiaDecayer::IsHandled(int pdg_code) const
225 {
226 // does not handle requests to decay baryon resonances
227 
228  bool is_handled = (!utils::res::IsBaryonResonance(pdg_code));
229 
230  LOG("Pythia6Decay", pDEBUG)
231  << "Can decay particle with PDG code = " << pdg_code
232  << "? " << ((is_handled)? "Yes" : "No");
233 
234  return is_handled;
235 }
236 //____________________________________________________________________________
237 void PythiaDecayer::InhibitDecay(int pdg_code, TDecayChannel * dc) const
238 {
239  if(! this->IsHandled(pdg_code)) return;
240 
241  int kc = fPythia->Pycomp(pdg_code);
242 
243  if(!dc) {
244  LOG("Pythia6Decay", pINFO)
245  << "Switching OFF ALL decay channels for particle = " << pdg_code;
246  fPythia->SetMDCY(kc, 1,0);
247  return;
248  }
249 
250  LOG("Pythia6Decay", pINFO)
251  << "Switching OFF decay channel = " << dc->Number()
252  << " for particle = " << pdg_code;
253 
254  int ichannel = this->FindPythiaDecayChannel(kc, dc);
255  if(ichannel != -1) {
256  fPythia->SetMDME(ichannel,1,0); // switch-off
257  }
258 }
259 //____________________________________________________________________________
260 void PythiaDecayer::UnInhibitDecay(int pdg_code, TDecayChannel * dc) const
261 {
262  if(! this->IsHandled(pdg_code)) return;
263 
264  int kc = fPythia->Pycomp(pdg_code);
265 
266  if(!dc) {
267  LOG("Pythia6Decay", pINFO)
268  << "Switching ON all PYTHIA decay channels for particle = " << pdg_code;
269 
270  fPythia->SetMDCY(kc, 1,1);
271 
272  int first_channel = fPythia->GetMDCY(kc,2);
273  int last_channel = fPythia->GetMDCY(kc,2) + fPythia->GetMDCY(kc,3) - 1;
274 
275  for(int ichannel = first_channel;
276  ichannel < last_channel; ichannel++) {
277  fPythia->SetMDME(ichannel,1,1); // switch-on
278  }
279  return;
280  }//!dc
281 
282  LOG("Pythia6Decay", pINFO)
283  << "Switching OFF decay channel = " << dc->Number()
284  << " for particle = " << pdg_code;
285 
286  int ichannel = this->FindPythiaDecayChannel(kc, dc);
287  if(ichannel != -1) {
288  fPythia->SetMDME(ichannel,1,1); // switch-on
289  }
290 }
291 //____________________________________________________________________________
293 {
294 // Sum of branching ratios for enabled channels
295 //
296  double sumbr=0.;
297 
298  int first_channel = fPythia->GetMDCY(kc,2);
299  int last_channel = fPythia->GetMDCY(kc,2) + fPythia->GetMDCY(kc,3) - 1;
300 
301  bool has_inhibited_channels=false;
302 
303  // loop over pythia decay channels
304  for(int ichannel = first_channel;
305  ichannel < last_channel; ichannel++) {
306 
307  bool enabled = (fPythia->GetMDME(ichannel,1) == 1);
308  if (!enabled) {
309  has_inhibited_channels = true;
310  } else {
311  sumbr += fPythia->GetBRAT(ichannel);
312  }
313  }
314 
315  if(!has_inhibited_channels) return 1.;
316  LOG("Pythia6Decay", pINFO) << "Sum{BR} = " << sumbr;
317 
318  return sumbr;
319 }
320 //____________________________________________________________________________
321 int PythiaDecayer::FindPythiaDecayChannel(int kc, TDecayChannel* dc) const
322 {
323  if(!dc) return -1;
324 
325  int first_channel = fPythia->GetMDCY(kc,2);
326  int last_channel = fPythia->GetMDCY(kc,2) + fPythia->GetMDCY(kc,3) - 1;
327 
328  bool found_match = false;
329 
330  // loop over pythia decay channels
331  for(int ichannel = first_channel;
332  ichannel < last_channel; ichannel++) {
333 
334  // does the current pythia channel matches the input TDecayChannel?
335  LOG("Pythia6Decay", pINFO)
336  << "\nComparing PYTHIA's channel = " << ichannel
337  << " with TDecayChannel = " << dc->Number();
338 
339  found_match = this->MatchDecayChannels(ichannel,dc);
340  if(found_match) {
341  LOG("Pythia6Decay", pNOTICE)
342  << " ** TDecayChannel id = " << dc->Number()
343  << " corresponds to PYTHIA6 channel id = " << ichannel;
344  return ichannel;
345  }//match?
346  }//loop pythia decay ch.
347 
348  LOG("Pythia6Decay", pWARN)
349  << " ** No PYTHIA6 decay channel match found for "
350  << "TDecayChannel id = " << dc->Number();
351 
352  return -1;
353 }
354 //____________________________________________________________________________
355 bool PythiaDecayer::MatchDecayChannels(int ichannel, TDecayChannel* dc) const
356 {
357  // num. of daughters in the input TDecayChannel & the input PYTHIA ichannel
358  int nd = dc->NDaughters();
359 
360  int py_nd = 0;
361  for (int i = 1; i <= 5; i++) {
362  if(fPythia->GetKFDP(ichannel,i) != 0) py_nd++;
363  }
364 
365  LOG("Pythia6Decay", pDEBUG)
366  << "NDaughters: PYTHIA = " << py_nd << ", ROOT's TDecayChannel = " << nd;
367 
368  if(nd != py_nd) return false;
369 
370  //
371  // if the two channels have the same num. of daughters, then compare them
372  //
373 
374  // store decay daughters for the input TDecayChannel
375  vector<int> dc_daughter(nd);
376  int k=0;
377  for( ; k < nd; k++) {
378  dc_daughter[k] = dc->DaughterPdgCode(k);
379  }
380  // store decay daughters for the input PYTHIA's ichannel
381  vector<int> py_daughter(nd);
382  k=0;
383  for(int i = 1; i <= 5; i++) {
384  if(fPythia->GetKFDP(ichannel,i) == 0) continue;
385  py_daughter[k] = fPythia->GetKFDP(ichannel,i);
386  k++;
387  }
388 
389  // sort both daughter lists
390  sort( dc_daughter.begin(), dc_daughter.end() );
391  sort( py_daughter.begin(), py_daughter.end() );
392 
393  // compare
394  for(int i = 0; i < nd; i++) {
395  if(dc_daughter[i] != py_daughter[i]) return false;
396  }
397  return true;
398 }
399 //____________________________________________________________________________
virtual bool ToBeDecayed(int pdgc, GHepStatus_t ist) const
Definition: Decayer.cxx:51
bool IsHandled(int pdgc) const
void UnInhibitDecay(int pdgc, TDecayChannel *ch=0) const
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
static constexpr double s
Definition: Units.h:95
double Mass(void) const
Mass that corresponds to the PDG code.
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
GHepStatus_t Status(void) const
Definition: GHepParticle.h:64
double Energy(void) const
Get energy.
Definition: GHepParticle.h:92
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
int Pdg(void) const
Definition: GHepParticle.h:63
string Name(void) const
Name that corresponds to the PDG code.
bool IsHadron(int pdgc)
Definition: PDGUtils.cxx:392
bool MatchDecayChannels(int ic, TDecayChannel *ch) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void Initialize(void) const
bool fRunBefHadroTransp
is invoked before or after FSI?
Definition: Decayer.h:57
bool Decay(int ip, GHepRecord *event) const
TPythia6 * fPythia
PYTHIA6 wrapper class.
Definition: PythiaDecayer.h:51
#define pINFO
Definition: Messenger.h:62
#define pWARN
Definition: Messenger.h:60
void pydecy_(int *)
double SumOfBranchingRatios(int kc) const
void py1ent_(int *, int *, double *, double *, double *)
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:126
Base class for decayer classes. Implements common configuration, allowing users to toggle on/off flag...
Definition: Decayer.h:34
void InhibitDecay(int pdgc, TDecayChannel *ch=0) const
int FindPythiaDecayChannel(int kc, TDecayChannel *ch) const
static constexpr double mm
Definition: Units.h:65
void ProcessEventRecord(GHepRecord *event) const
bool IsBaryonResonance(int pdgc)
is input a baryon resonance?
static constexpr double fm
Definition: Units.h:75
#define pNOTICE
Definition: Messenger.h:61
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:84
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
static constexpr double m
Definition: Units.h:71
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
enum genie::EGHepStatus GHepStatus_t
#define pDEBUG
Definition: Messenger.h:63
double Py(void) const
Get Py.
Definition: GHepParticle.h:89