GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gtestDecay.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gtestDecay
5 
6 \brief test program used for testing / debugging the GENIE particle decayers
7 
8 \author Costas Andreopoulos <c.andreopoulos \at cern.ch>
9  University of Liverpool
10 
11 \created June 20, 2004
12 
13 \cpright Copyright (c) 2003-2024, The GENIE Collaboration
14  For the full text of the license visit http://copyright.genie-mc.org
15 
16 */
17 //____________________________________________________________________________
18 
19 #include <ostream>
20 #include <iomanip>
21 
22 #include <RVersion.h>
23 #include <TClonesArray.h>
24 #include <TParticlePDG.h>
25 #include <TIterator.h>
26 
31 #include "Physics/Decay/Decayer.h"
39 
40 using namespace genie;
41 
42 using std::ostream;
43 using std::endl;
44 using std::setw;
45 using std::setprecision;
46 using std::setfill;
47 using std::ios;
48 
49 ostream & operator<< (ostream & stream, const TClonesArray * particle_list);
50 ostream & operator<< (ostream & stream, const GHepParticle * particle);
51 
52 void TestPythiaTauDecays(void);
53 void Decay(const Decayer * decayer, int pdgc, double E, int ndecays);
54 
55 //__________________________________________________________________________
56 int main(int /*argc*/, char ** /*argv*/)
57 {
59 
60  return 0;
61 }
62 //__________________________________________________________________________
64 {
65  // Get the pythia decayer
66  LOG("test",pINFO)
67  << "Asking the AlgFactory for a genie::PythiaDecayer\\Default instance";
69  const Decayer * pdecayer =
70  dynamic_cast<const Decayer *> (
71  algf->GetAlgorithm("genie::PythiaDecayer","Default"));
72 
73  // Decayer config print-out
74  LOG("test",pINFO) << "Algorithm name = " << pdecayer->Id().Name();
75  LOG("test",pINFO) << "Parameter set = " << pdecayer->Id().Config();
76 
77  const Registry & conf_registry = pdecayer->GetConfig();
78  LOG("test", pINFO) << conf_registry;
79 
80  double E = 10;
81  int ndec = 3;
82 
83  // inhibit pi0 decay
84  pdecayer->InhibitDecay(kPdgPi0);
85 
86  // decay tau-
87  Decay(pdecayer, kPdgTau, E, ndec);
88 
89  // now inhibit all but the tau- --> nu_mu_bar + mu- + nu_tau decay channel
90  // (see $GENIE/src/contrib/misc/print_decay_channels.C)
91  // and perform some more decays
92  LOG("test",pINFO)
93  << "\n\n"
94  << " **** Inhibiting all but the `tau- --> nu_mu_bar + mu- + nu_tau' decay channel"
95  << "\n\n";
96 
97  PDGLibrary * pdglib = PDGLibrary::Instance();
98 
99  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(0));
100  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(2));
101  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(3));
102  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(4));
103  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(5));
104  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(6));
105  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(7));
106  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(8));
107  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(9));
108  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(10));
109  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(11));
110  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(12));
111  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(13));
112  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(14));
113  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(15));
114  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(16));
115  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(17));
116  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(18));
117  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(19));
118  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(20));
119  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(21));
120  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(22));
121  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(23));
122  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(24));
123  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(25));
124  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(26));
125  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(27));
126  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(28));
127  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(29));
128  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(30));
129  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(31));
130  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(32));
131  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(33));
132  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(34));
133  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(35));
134  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(36));
135  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(37));
136  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(38));
137  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(39));
138  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(40));
139  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(41));
140  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(42));
141  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(43));
142  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(44));
143  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(45));
144  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(46));
145  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(47));
146  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(48));
147  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(49));
148  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(50));
149  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(51));
150  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(52));
151  pdecayer->InhibitDecay(kPdgTau, pdglib->Find(kPdgTau)->DecayChannel(53));
152 
153  // a few more decays
154  Decay(pdecayer, kPdgTau, E, ndec);
155 
156  // restore all decay channels
157  LOG("test",pINFO)
158  << "\n\n"
159  << " **** Restoring all tau- decay channels"
160  << "\n\n";
161 
162  pdecayer->UnInhibitDecay(kPdgTau);
163 
164  // a few more decays
165  Decay(pdecayer, kPdgTau, E, ndec);
166 
167  // now inhibit all decay channels and try to decay!
168  LOG("test",pINFO)
169  << "\n\n"
170  << " **** Inhibit all tau- decay channels"
171  << "\n\n";
172 
173  pdecayer->InhibitDecay(kPdgTau);
174 
175  // a few more decays
176  Decay(pdecayer, kPdgTau, E, ndec);
177 }
178 //__________________________________________________________________________
179 void Decay(const Decayer * decayer, int pdgc, double E, int ndecays)
180 {
181  DecayerInputs_t dinp;
182 
183  TLorentzVector p4;
184  p4.SetE(E);
185  p4.SetTheta(0.);
186  p4.SetPhi(0.);
187 
188  dinp.PdgCode = pdgc;
189  dinp.P4 = &p4;
190 
191  PDGLibrary * pdglib = PDGLibrary::Instance();
192  TParticlePDG * pp = pdglib->Find(pdgc);
193  if(!pp) return;
194 
195  LOG("test",pINFO)
196  << "Decaying a " << pp->GetName() << " with E = " << p4.Energy();
197 
198  // Perform the decay a few times & print-out the decay products
199  for(int idec = 0; idec < ndecays; idec++) {
200 
201  // Decay
202  TClonesArray * particle_list = decayer->Decay(dinp);
203 
204  if(!particle_list) {
205  LOG("test",pWARN)
206  << "\n ** Decay nu.: " << idec << " ==> NULL particle list";
207  continue;
208  }
209 
210  // Print-out
211  LOG("test",pINFO)
212  << "\n ** Decay nu.: " << idec
213  << " (weight = " << decayer->Weight() << ") : \n "
214  << particle_list;
215 
216  // Clean-up
217  particle_list->Delete();
218  delete particle_list;
219  }//ndecays
220 
221 }
222 //__________________________________________________________________________
223 ostream & operator<< (ostream & stream, const TClonesArray * particle_list)
224 {
225  GHepParticle * p = 0;
226  TObjArrayIter particle_iter(particle_list);
227 
228 
229  stream
230  << setfill(' ') << setw(10) << "name "
231  << setfill(' ') << setw(10) << "PDG"
232  << setfill(' ') << setw(10) << "Status"
233  << setfill(' ') << setw(15) << "E (GeV)"
234  << setfill(' ') << setw(15) << "Px (GeV/c)"
235  << setfill(' ') << setw(15) << "Py (GeV/c)"
236  << setfill(' ') << setw(15) << "Pz (GeV/c)"
237  << setfill(' ') << setw(15) << "t (mm/c)"
238  << setfill(' ') << setw(15) << "x (mm)"
239  << setfill(' ') << setw(15) << "y (mm)"
240  << setfill(' ') << setw(15) << "z (mm)"
241  << endl;
242 
243  while( (p = (GHepParticle *) particle_iter.Next()) ) stream << p;
244 
245  stream << setfill('-') << setw(100) << "|";
246 
247  return stream;
248 }
249 //__________________________________________________________________________
250 ostream & operator<< (ostream & stream, const GHepParticle * p)
251 {
252  stream
253  << std::scientific << setprecision(6);
254 
255  stream
256  << setfill(' ') << setw(10) << p->Name()
257  << setfill(' ') << setw(10) << p->Pdg()
258  << setfill(' ') << setw(10) << p->Status()
259  << setfill(' ') << setw(15) << p->Energy()
260  << setfill(' ') << setw(15) << p->Px()
261  << setfill(' ') << setw(15) << p->Py()
262  << setfill(' ') << setw(15) << p->Pz()
263  << setfill(' ') << setw(15) << p->Vt() /(units::mm)
264  << setfill(' ') << setw(15) << p->Vx() /(units::mm)
265  << setfill(' ') << setw(15) << p->Vy() /(units::mm)
266  << setfill(' ') << setw(15) << p->Vz() /(units::mm)
267  << endl;
268 
269  return stream;
270 }
271 //__________________________________________________________________________
272 
virtual void UnInhibitDecay(int pdgc, TDecayChannel *dc=0) const =0
void TestPythiaTauDecays(void)
Definition: gtestDecay.cxx:63
void Decay(const Decayer *decayer, int pdgc, double E, int ndecays)
Definition: gtestDecay.cxx:179
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
int main(int argc, char **argv)
Definition: gAtmoEvGen.cxx:327
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:246
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
double Vt(void) const
Get production time.
Definition: GHepParticle.h:97
int Pdg(void) const
Definition: GHepParticle.h:63
string Name(void) const
Name that corresponds to the PDG code.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const int kPdgTau
Definition: PDGCodes.h:39
string Name(void) const
Definition: AlgId.h:44
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:75
const int kPdgPi0
Definition: PDGCodes.h:160
#define pINFO
Definition: Messenger.h:62
#define pWARN
Definition: Messenger.h:60
virtual void InhibitDecay(int pdgc, TDecayChannel *dc=0) const =0
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
static AlgFactory * Instance()
Definition: AlgFactory.cxx:64
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
double Vz(void) const
Get production z.
Definition: GHepParticle.h:96
ostream & operator<<(ostream &stream, const AlgConfigPool &config_pool)
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
Base class for decayer classes. Implements common configuration, allowing users to toggle on/off flag...
Definition: Decayer.h:34
static constexpr double mm
Definition: Units.h:65
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
double Vy(void) const
Get production y.
Definition: GHepParticle.h:95
The GENIE Algorithm Factory.
Definition: AlgFactory.h:39
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 Vx(void) const
Get production x.
Definition: GHepParticle.h:94
string Config(void) const
Definition: AlgId.h:45
double Py(void) const
Get Py.
Definition: GHepParticle.h:89