GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Decayer.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 <algorithm>
12 #include <sstream>
13 
14 #include <TParticlePDG.h>
15 
21 #include "Physics/Decay/Decayer.h"
22 
23 using std::count;
24 using std::ostringstream;
25 
26 using namespace genie;
27 //___________________________________________________________________________
30 {
31 
32 }
33 //___________________________________________________________________________
34 Decayer::Decayer(string name) :
36 {
37 
38 }
39 //___________________________________________________________________________
40 Decayer::Decayer(string name, string config) :
41 EventRecordVisitorI(name, config)
42 {
43 
44 }
45 //___________________________________________________________________________
47 {
48 
49 }
50 //___________________________________________________________________________
51 bool Decayer::ToBeDecayed(int pdg_code, GHepStatus_t status_code) const
52 {
53  // Check whether it is "unstable" (definition can vary)
54 
55  bool is_unstable = this->IsUnstable(pdg_code);
56 
57  LOG("Decay", pDEBUG)
58  << "Particle is unstable? "
59  << ((is_unstable) ? "Yes" : "No");
60 
61  if(!is_unstable) return false;
62 
63  // Check whether the given unstable particle
64  // has the appropriate status code to be decayed
65 
66  bool to_be_decayed = false;
67 
68  if(fRunBefHadroTransp) {
69  to_be_decayed =
70  (status_code == kIStHadronInTheNucleus ||
71  status_code == kIStPreDecayResonantState ||
72  status_code == kIStStableFinalState);
73  }
74  else {
75  to_be_decayed =
76  (status_code == kIStStableFinalState);
77  }
78 
79  LOG("Decay", pDEBUG)
80  << "Particle to be decayed "
81  << "[" << ((fRunBefHadroTransp) ? "Before" : "After") << " FSI]? "
82  << ((to_be_decayed) ? "Yes" : "No");
83 
84  return to_be_decayed;
85 }
86 //___________________________________________________________________________
87 bool Decayer::IsUnstable(int pdg_code) const
88 {
89  // ROOT's TParticlepdg::Lifetime() does not work properly
90  // do something else instead (temporarily)
91  //
92  // TParticlePDG * ppdg = PDGLibrary::Instance()->Find(pdg_code);
93  //if( ppdg->Lifetime() < fMaxLifetime ) { /* ... */ };
94  //
95 
96  // <temp/>
97  if( fRunBefHadroTransp ) {
98  //
99  // Run *before* the hadron transport MC
100  // At this point we decay only baryon resonances
101  //
102  bool decay = utils::res::IsBaryonResonance(pdg_code);
103  return decay;
104  }
105  else {
106  //
107  // Run *after* the hadron transport MC
108  // At this point we decay only particles in the fParticlesToDecay
109  // PDGCodeList (filled in from config inputs)
110  //
111  bool decay = fParticlesToDecay.ExistsInPDGCodeList(pdg_code);
112  return decay;
113  }
114  // </temp>
115 
116  return false;
117 }
118 //___________________________________________________________________________
119 void Decayer::Configure(const Registry & config)
120 {
121  Algorithm::Configure(config);
122  this->LoadConfig();
123 
124  fAllowReconfig = false;
125 }
126 //___________________________________________________________________________
127 void Decayer::Configure(string config)
128 {
129  Algorithm::Configure(config);
130  this->LoadConfig();
131 
132  fAllowReconfig = false;
133 }
134 //___________________________________________________________________________
136 {
137  // Get the specified maximum lifetime tmax (decay with lifetime < tmax)
138  //
139  //fMaxLifetime = fConfig->GetDoubleDef("MaxLifetime", 1e-9);
140 
141  // Check whether to generate weighted or unweighted particle decays
142  fGenerateWeighted = false ;
143  //this->GetParam("GenerateWeighted", fGenerateWeighted, false);
144 
145  // Check whether the module is being run before or after the hadron
146  // transport (intranuclear rescattering) module.
147  //
148  // If it is run before the hadron transport (and after the hadronization)
149  // step it should decay only "unstable" particles (marked as hadrons in
150  // the nucleus) which would typically decay within the time required to
151  // exit the nucleus - so, the algorithm wouldn't decay particles that
152  // have to be rescattered first. In case that the generated event is off
153  // a free nucleon target, thi instance of the algorithm should do nothing.
154  //
155  // If it is run after the hadon transport, then it should decay all the
156  // 'unstable' particles marked as 'present in the final state' and which
157  // should be decay before the event is passed to the detector particle
158  // transport MC.
159  //
160 
161  this->GetParam("RunBeforeHadronTransport", fRunBefHadroTransp) ;
162 
163  // Allow user to specify a list of particles to be decayed
164  //
165  RgKeyList klist = GetConfig().FindKeys("DecayParticleWithCode=");
166  RgKeyList::const_iterator kiter = klist.begin();
167  for( ; kiter != klist.end(); ++kiter) {
168  RgKey key = *kiter;
169  bool decay = GetConfig().GetBool(key);
170  vector<string> kv = utils::str::Split(key,"=");
171  assert(kv.size()==2);
172  int pdgc = atoi(kv[1].c_str());
173  TParticlePDG * p = PDGLibrary::Instance()->Find(pdgc);
174  if ( ! p ) {
175  LOG("Decay",pFATAL) << "No PDGLibrary entry for pdgc=" << pdgc
176  << " (" << kv[1].c_str()
177  << "), check CommonDecay.xml";
178  continue;
179  }
180  if(decay) {
181  LOG("Decay", pDEBUG)
182  << "Configured to decay " << p->GetName();
184  this->UnInhibitDecay(pdgc);
185  }
186  else {
187  LOG("Decay", pDEBUG)
188  << "Configured to inhibit decays for " << p->GetName();
190  this->InhibitDecay(pdgc);
191  }// decay?
192  }// key iterator
193 
194  // Allow user to inhibit certain decay channels
195  //
196  klist = GetConfig().FindKeys("InhibitDecay/");
197  kiter = klist.begin();
198  for( ; kiter != klist.end(); ++kiter) {
199  RgKey key = *kiter;
200  if(GetConfig().GetBool(key)) {
201  string filtkey = utils::str::FilterString("InhibitDecay/", key);
202  vector<string> kv = utils::str::Split(filtkey,",");
203  assert(kv.size()==2);
204  int pdgc = atoi(utils::str::FilterString("Particle=",kv[0]).c_str());
205  int dc = atoi(utils::str::FilterString("Channel=", kv[1]).c_str());
206  TParticlePDG * p = PDGLibrary::Instance()->Find(pdgc);
207  if(!p) continue;
208  LOG("Decay", pINFO)
209  << "Configured to inhibit " << p->GetName()
210  << "'s decay channel " << dc;
211  this->InhibitDecay(pdgc, p->DecayChannel(dc));
212  }//val[key]=true?
213  }//key iterator
214 
215 
216  sort(fParticlesToDecay.begin(), fParticlesToDecay.end());
217  sort(fParticlesNotToDecay.begin(), fParticlesNotToDecay.end());
218 
219  // Print-out for only one of the two instances of this module
220  if(!fRunBefHadroTransp) {
221  LOG("Decay", pNOTICE)
222  << "\nConfigured to decay: " << fParticlesToDecay
223  << "\nConfigured to inhibit decays of: " << fParticlesNotToDecay
224  << "\n";
225  }
226 }
227 //___________________________________________________________________________
virtual bool ToBeDecayed(int pdgc, GHepStatus_t ist) const
Definition: Decayer.cxx:51
virtual bool IsUnstable(int pdgc) const
Definition: Decayer.cxx:87
virtual ~Decayer()
Definition: Decayer.cxx:46
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
void Configure(const Registry &config)
Definition: Decayer.cxx:119
#define pFATAL
Definition: Messenger.h:56
virtual void UnInhibitDecay(int pdgc, TDecayChannel *dc=0) const =0
bool ExistsInPDGCodeList(int pdg_code) const
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:246
RgKeyList FindKeys(RgKey key_part) const
create list with all keys containing &#39;key_part&#39;
Definition: Registry.cxx:840
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool fRunBefHadroTransp
is invoked before or after FSI?
Definition: Decayer.h:57
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
#define pINFO
Definition: Messenger.h:62
bool fGenerateWeighted
generate weighted or unweighted decays?
Definition: Decayer.h:56
PDGCodeList fParticlesNotToDecay
list of particles for which decay is inhibited
Definition: Decayer.h:59
vector< RgKey > RgKeyList
Definition: Registry.h:50
PDGCodeList fParticlesToDecay
list of particles to be decayed
Definition: Decayer.h:58
virtual void InhibitDecay(int pdgc, TDecayChannel *dc=0) const =0
virtual void LoadConfig(void)
Definition: Decayer.cxx:135
string FilterString(string filt, string input)
Definition: StringUtils.cxx:79
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:36
string RgKey
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
RgBool GetBool(RgKey key) const
Definition: Registry.cxx:460
bool IsBaryonResonance(int pdgc)
is input a baryon resonance?
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
#define pNOTICE
Definition: Messenger.h:61
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
enum genie::EGHepStatus GHepStatus_t
#define pDEBUG
Definition: Messenger.h:63