PythiaDecayer.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 
6  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
7  University of Liverpool & STFC Rutherford Appleton Lab
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  TLorentzVector decay_particle_x4 = *(decay_particle->X4());
106  int decay_particle_pdg_code = decay_particle->Pdg();
107 
108  // Convert to PYTHIA6 particle code and check whether decay is inhibited
109  int kc = fPythia->Pycomp(decay_particle_pdg_code);
110  int mdcy = fPythia->GetMDCY(kc, 1);
111  if(mdcy == 0) {
112  LOG("Pythia6Decay", pNOTICE)
113  << (PDGLibrary::Instance())->Find(decay_particle_pdg_code)->GetName()
114  << " decays are inhibited!";
115  return false;
116  }
117 
118  // Get sub of BRs and compute weight if decay channels were inhibited
119  double sumbr = this->SumOfBranchingRatios(kc);
120  if(sumbr <= 0) {
121  LOG("Pythia6Decay", pWARN)
122  << "The sum of enabled "
123  << (PDGLibrary::Instance())->Find(decay_particle_pdg_code)->GetName()
124  << " decay channel branching ratios is non-positive!";
125  return false;
126  }
127  fWeight = 1./sumbr; // update weight to account for inhibited channels
128 
129  // Run PYTHIA6 decay
130  int ip = 0;
131  double E = decay_particle_p4.Energy();
132  double theta = decay_particle_p4.Theta();
133  double phi = decay_particle_p4.Phi();
134  fPythia->SetMSTJ(22,1);
135  py1ent_(&ip, &decay_particle_pdg_code, &E, &theta, &phi);
136 
137  // Get decay products
138  fPythia->GetPrimaries();
139  TClonesArray * impl = (TClonesArray *) fPythia->ImportParticles("All");
140  if(!impl) {
141  LOG("Pythia6Decay", pWARN)
142  << "TPythia6::ImportParticles() returns a null list!";
143  return false;
144  }
145 
146  // Copy the PYTHIA6 container to the GENIE event record
147 
148  // Check whether the interaction is off a nuclear target or free nucleon
149  // Depending on whether this module is run before or after the hadron
150  // transport module it would affect the daughters status code
151  GHepParticle * target_nucleus = event->TargetNucleus();
152  bool in_nucleus = (target_nucleus!=0);
153 
154  TMCParticle * p = 0;
155  TIter particle_iter(impl);
156  while( (p = (TMCParticle *) particle_iter.Next()) ) {
157  // Convert from TMCParticle to GHepParticle
159  p->GetKF(), // pdg
160  GHepStatus_t(p->GetKS()), // status
161  p->GetParent(), // first parent
162  0, // second parent
163  p->GetFirstChild(), // first daughter
164  p->GetLastChild(), // second daughter
165  p->GetPx(), // px
166  p->GetPy(), // py
167  p->GetPz(), // pz
168  p->GetEnergy(), // e
169  p->GetVx(), // x
170  p->GetVy(), // y
171  p->GetVz(), // z
172  p->GetTime() // t
173  );
174 
175  if(mcp.Status()==kIStNucleonTarget) continue; // mother particle, already in GHEP
176 
177  int daughter_pdg_code = mcp.Pdg();
178  SLOG("Pythia6Decay", pINFO)
179  << "Adding daughter particle wit PDG code = "
180  << daughter_pdg_code << ", m = " << mcp.Mass()
181  << " GeV, E = " << mcp.Energy() << " GeV)";
182 
183  bool is_hadron = pdg::IsHadron(daughter_pdg_code);
184  bool hadron_in_nuc = (in_nucleus && is_hadron && fRunBefHadroTransp);
185 
186  GHepStatus_t daughter_status_code = (hadron_in_nuc) ?
188 
189  TLorentzVector daughter_p4(
190  mcp.Px(),mcp.Py(),mcp.Pz(),mcp.Energy());
191  event->AddParticle(
192  daughter_pdg_code, daughter_status_code,
193  decay_particle_id,-1,-1,-1,
194  daughter_p4, decay_particle_x4);
195  }
196 
197  // Update the event weight for each weighted particle decay
198  double weight = event->Weight() * fWeight;
199  event->SetWeight(weight);
200 
201  // Mark input particle as a 'decayed state' & add its daughter links
202  decay_particle->SetStatus(kIStDecayedState);
203 
204  return true;
205 }
206 //____________________________________________________________________________
208 {
209  fPythia = TPythia6::Instance();
210  fWeight = 1.;
211 
212  // sync GENIE/PYTHIA6 seeds
214 }
215 //____________________________________________________________________________
216 bool PythiaDecayer::IsHandled(int pdg_code) const
217 {
218 // does not handle requests to decay baryon resonances
219 
220  bool is_handled = (!utils::res::IsBaryonResonance(pdg_code));
221 
222  LOG("Pythia6Decay", pDEBUG)
223  << "Can decay particle with PDG code = " << pdg_code
224  << "? " << ((is_handled)? "Yes" : "No");
225 
226  return is_handled;
227 }
228 //____________________________________________________________________________
229 void PythiaDecayer::InhibitDecay(int pdg_code, TDecayChannel * dc) const
230 {
231  if(! this->IsHandled(pdg_code)) return;
232 
233  int kc = fPythia->Pycomp(pdg_code);
234 
235  if(!dc) {
236  LOG("Pythia6Decay", pINFO)
237  << "Switching OFF ALL decay channels for particle = " << pdg_code;
238  fPythia->SetMDCY(kc, 1,0);
239  return;
240  }
241 
242  LOG("Pythia6Decay", pINFO)
243  << "Switching OFF decay channel = " << dc->Number()
244  << " for particle = " << pdg_code;
245 
246  int ichannel = this->FindPythiaDecayChannel(kc, dc);
247  if(ichannel != -1) {
248  fPythia->SetMDME(ichannel,1,0); // switch-off
249  }
250 }
251 //____________________________________________________________________________
252 void PythiaDecayer::UnInhibitDecay(int pdg_code, TDecayChannel * dc) const
253 {
254  if(! this->IsHandled(pdg_code)) return;
255 
256  int kc = fPythia->Pycomp(pdg_code);
257 
258  if(!dc) {
259  LOG("Pythia6Decay", pINFO)
260  << "Switching ON all PYTHIA decay channels for particle = " << pdg_code;
261 
262  fPythia->SetMDCY(kc, 1,1);
263 
264  int first_channel = fPythia->GetMDCY(kc,2);
265  int last_channel = fPythia->GetMDCY(kc,2) + fPythia->GetMDCY(kc,3) - 1;
266 
267  for(int ichannel = first_channel;
268  ichannel < last_channel; ichannel++) {
269  fPythia->SetMDME(ichannel,1,1); // switch-on
270  }
271  return;
272  }//!dc
273 
274  LOG("Pythia6Decay", pINFO)
275  << "Switching OFF decay channel = " << dc->Number()
276  << " for particle = " << pdg_code;
277 
278  int ichannel = this->FindPythiaDecayChannel(kc, dc);
279  if(ichannel != -1) {
280  fPythia->SetMDME(ichannel,1,1); // switch-on
281  }
282 }
283 //____________________________________________________________________________
285 {
286 // Sum of branching ratios for enabled channels
287 //
288  double sumbr=0.;
289 
290  int first_channel = fPythia->GetMDCY(kc,2);
291  int last_channel = fPythia->GetMDCY(kc,2) + fPythia->GetMDCY(kc,3) - 1;
292 
293  bool has_inhibited_channels=false;
294 
295  // loop over pythia decay channels
296  for(int ichannel = first_channel;
297  ichannel < last_channel; ichannel++) {
298 
299  bool enabled = (fPythia->GetMDME(ichannel,1) == 1);
300  if (!enabled) {
301  has_inhibited_channels = true;
302  } else {
303  sumbr += fPythia->GetBRAT(ichannel);
304  }
305  }
306 
307  if(!has_inhibited_channels) return 1.;
308  LOG("Pythia6Decay", pINFO) << "Sum{BR} = " << sumbr;
309 
310  return sumbr;
311 }
312 //____________________________________________________________________________
313 int PythiaDecayer::FindPythiaDecayChannel(int kc, TDecayChannel* dc) const
314 {
315  if(!dc) return -1;
316 
317  int first_channel = fPythia->GetMDCY(kc,2);
318  int last_channel = fPythia->GetMDCY(kc,2) + fPythia->GetMDCY(kc,3) - 1;
319 
320  bool found_match = false;
321 
322  // loop over pythia decay channels
323  for(int ichannel = first_channel;
324  ichannel < last_channel; ichannel++) {
325 
326  // does the current pythia channel matches the input TDecayChannel?
327  LOG("Pythia6Decay", pINFO)
328  << "\nComparing PYTHIA's channel = " << ichannel
329  << " with TDecayChannel = " << dc->Number();
330 
331  found_match = this->MatchDecayChannels(ichannel,dc);
332  if(found_match) {
333  LOG("Pythia6Decay", pNOTICE)
334  << " ** TDecayChannel id = " << dc->Number()
335  << " corresponds to PYTHIA6 channel id = " << ichannel;
336  return ichannel;
337  }//match?
338  }//loop pythia decay ch.
339 
340  LOG("Pythia6Decay", pWARN)
341  << " ** No PYTHIA6 decay channel match found for "
342  << "TDecayChannel id = " << dc->Number();
343 
344  return -1;
345 }
346 //____________________________________________________________________________
347 bool PythiaDecayer::MatchDecayChannels(int ichannel, TDecayChannel* dc) const
348 {
349  // num. of daughters in the input TDecayChannel & the input PYTHIA ichannel
350  int nd = dc->NDaughters();
351 
352  int py_nd = 0;
353  for (int i = 1; i <= 5; i++) {
354  if(fPythia->GetKFDP(ichannel,i) != 0) py_nd++;
355  }
356 
357  LOG("Pythia6Decay", pDEBUG)
358  << "NDaughters: PYTHIA = " << py_nd << ", ROOT's TDecayChannel = " << nd;
359 
360  if(nd != py_nd) return false;
361 
362  //
363  // if the two channels have the same num. of daughters, then compare them
364  //
365 
366  // store decay daughters for the input TDecayChannel
367  vector<int> dc_daughter(nd);
368  int k=0;
369  for( ; k < nd; k++) {
370  dc_daughter[k] = dc->DaughterPdgCode(k);
371  }
372  // store decay daughters for the input PYTHIA's ichannel
373  vector<int> py_daughter(nd);
374  k=0;
375  for(int i = 1; i <= 5; i++) {
376  if(fPythia->GetKFDP(ichannel,i) == 0) continue;
377  py_daughter[k] = fPythia->GetKFDP(ichannel,i);
378  k++;
379  }
380 
381  // sort both daughter lists
382  sort( dc_daughter.begin(), dc_daughter.end() );
383  sort( py_daughter.begin(), py_daughter.end() );
384 
385  // compare
386  for(int i = 0; i < nd; i++) {
387  if(dc_daughter[i] != py_daughter[i]) return false;
388  }
389  return true;
390 }
391 //____________________________________________________________________________
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
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
const Var weight
enum genie::EGHepStatus GHepStatus_t
const char * p
Definition: xmltok.h:285
Definition: config.py:1
TString ip
Definition: loadincs.C:5
double Mass(void) const
Mass that corresponds to the PDG code.
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:91
GHepStatus_t Status(void) const
Definition: GHepParticle.h:65
double Energy(void) const
Get energy.
Definition: GHepParticle.h:93
double Px(void) const
Get Px.
Definition: GHepParticle.h:89
int Pdg(void) const
Definition: GHepParticle.h:64
string Name(void) const
Name that corresponds to the PDG code.
bool IsHadron(int pdgc)
Definition: PDGUtils.cxx:355
const char * Find(const char *fname)
Definition: Icons.cxx:12
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:97
Float_t E
Definition: plot.C:20
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:63
#define pWARN
Definition: Messenger.h:61
int ichannel
Definition: plotTbData.C:13
void pydecy_(int *)
double SumOfBranchingRatios(int kc) const
void py1ent_(int *, int *, double *, double *, double *)
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
void SetStatus(GHepStatus_t s)
Definition: GHepParticle.h:127
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
assert(nhit_max >=nhit_nbins)
void ProcessEventRecord(GHepRecord *event) const
bool IsBaryonResonance(int pdgc)
is input a baryon resonance?
#define pNOTICE
Definition: Messenger.h:62
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
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
#define pDEBUG
Definition: Messenger.h:64
double Py(void) const
Get Py.
Definition: GHepParticle.h:90