GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gtestNucleonDecay.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gtestEventLoop
5 
6 \brief Example event loop. Use that as a template for your analysis code.
7 
8 \author Costas Andreopoulos <c.andreopoulos \at cern.ch>
9  University of Liverpool
10 
11 \created May 4, 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 <string>
20 
21 #include <TSystem.h>
22 #include <TFile.h>
23 #include <TTree.h>
24 #include <TIterator.h>
25 #include <TH1.h>
26 #include <TText.h>
27 
38 
39 
40 using std::string;
41 using namespace genie;
42 
43 void GetCommandLineArgs (int argc, char ** argv);
44 
47 
48 //___________________________________________________________________
49 int main(int argc, char ** argv)
50 {
51  GetCommandLineArgs (argc, argv);
52 
53  //-- open the ROOT file and get the TTree & its header
54  TTree * tree = 0;
55  NtpMCTreeHeader * thdr = 0;
56 
57  TFile file(gOptInpFilename.c_str(),"READ");
58 
59  tree = dynamic_cast <TTree *> ( file.Get("gtree") );
60  thdr = dynamic_cast <NtpMCTreeHeader *> ( file.Get("header") );
61 
62  if(!tree) return 1;
63 
64  // Output histogram file
65  std::string filename = std::string(file.GetName());
66  std::string histofilename = filename.substr(0,filename.size()-5) +
67  ".histo.root";
68  TFile *histofile = new TFile(histofilename.c_str(),"RECREATE");
69 
70  // Decayed nucleon histograms
71  TH1D* dnPdgHisto = new TH1D("DecayedNucleonPdg", "DecayedNucleonPdg", 5000, -2500, 2500);
72 
73  TH1D* dnPHisto = new TH1D("DecayedNucleonMomentum", "DecayedNucleonMomentum [GeV/c]", 100, 0., 0.5); // [GeV/c]
74 
75  TH1D* dnRemovalEnergyHisto = new TH1D("DecayedNucleonRemovalEnergy", "DecayedNucleonRemovalEnergy [GeV]", 100, 0., 0.05); // [GeV/c]
76 
77  // Histograms for decayed nucleon daughters
78  TH1D* dndaughtersNHisto = new TH1D("DecayedNucleonNDaughters", "DecayedNucleonNDaughters", 6, 0, 6);
79 
80  TH1D* dndaughter0PdgHisto = new TH1D("DecayedNucleonDaughter0Pdg", "DecayedNucleonDaughter0Pdg", 1000, -500, 500);
81  TH1D* dndaughter1PdgHisto = new TH1D("DecayedNucleonDaughter1Pdg", "DecayedNucleonDaughter1Pdg", 1000, -500, 500);
82  TH1D* dndaughter2PdgHisto = new TH1D("DecayedNucleonDaughter2Pdg", "DecayedNucleonDaughter2Pdg", 1000, -500, 500);
83 
84  TH1D* dndaughter0PHisto = new TH1D("DecayedNucleonDaughter0Momentum", "DecayedNucleonDaughter0Momentum [GeV/c]", 100, 0., 1.); // [GeV/c]
85  TH1D* dndaughter1PHisto = new TH1D("DecayedNucleonDaughter1Momentum", "DecayedNucleonDaughter1Momentum [GeV/c]", 100, 0., 1.); // [GeV/c]
86  TH1D* dndaughter2PHisto = new TH1D("DecayedNucleonDaughter2Momentum", "DecayedNucleonDaughter2Momentum [GeV/c]", 100, 0., 1.); // [GeV/c]
87 
88  // Histograms for stable final state particles
89  TH1D* finalparticlesNHisto = new TH1D("NFinalParticles", "NFinalParticles", 20, 0, 20);
90  TH1D* finalparticlesPdgHisto = new TH1D("FinalParticlesPdg", "FinalParticlesPdg", 5000, -2500, 2500);
91  TH1D* finalparticlesPHisto = new TH1D("FinalParticlesMomentum", "FinalParticlesMomentum [GeV/c]", 100, 0., 1.); // [GeV/c]
92 
93 
94  // address for input file
95  NtpMCEventRecord * mcrec = 0;
96  tree->SetBranchAddress("gmcrec", &mcrec);
97 
98  int nev = (gOptNEvt > 0) ?
99  TMath::Min(gOptNEvt, (int)tree->GetEntries()) :
100  (int) tree->GetEntries();
101 
102  //
103  // Loop over all events
104  //
105  // boolean to set, only for the first event, a few TTexts to be written into the histogram filename
106  bool first = true;
107  TText decayName = TText();
108  decayName.SetName("DecayName");
109  TText targetName = TText();
110  targetName.SetName("TargetName");
111 
112  // decayed nucleon index. In the event record, the index 0 element is typically the target, while the index 1 element is the decayed nucleon within the target
113  int dnIndex = 1;
114 
115  int ndaughters;
116  int nfinalparticles;
117 
118  for(int i = 0; i < nev; i++) {
119 
120  // get next tree entry
121  tree->GetEntry(i);
122 
123  // get the GENIE event
124  EventRecord & event = *(mcrec->event);
125 
126  LOG("testNucleonDecay", pNOTICE) << event;
127 
128  if (first) {
129  first = !first;
130 
131  NucleonDecayMode_t ndm = (NucleonDecayMode_t)event.Summary()->ExclTagPtr()->DecayMode();
132  int npdg = event.Summary()->InitStatePtr()->TgtPtr()->HitNucPdg();
133  string decayNameString = genie::utils::nucleon_decay::AsString(ndm,npdg);
134  decayName.SetTitle(decayNameString.c_str());
135  string targetNameString = event.Summary()->InitStatePtr()->TgtPtr()->AsString();
136  targetName.SetTitle(targetNameString.c_str());
137  }
138 
139  GHepParticle * p = 0;
140  TIter event_iter(&event);
141 
142  // decayed nucleon histograms
143  GHepParticle* dnpart = event.Particle(dnIndex);
144 
145  if (dnpart->Status() != 3) {
146  LOG("testNucleonDecay", pFATAL) << "Unexpected status code for candidate decayed nucleon: " << dnpart->Status() << ", exiting.";
147  exit(1);
148  }
149 
150  if (dnpart->FirstMother() != 0) {
151  LOG("testNucleonDecay", pFATAL) << "Unexpected first mother index for candidate decayed nucleon: " << dnpart->FirstMother() << ", exiting.";
152  exit(1);
153  }
154 
155  dnPdgHisto->Fill(dnpart->Pdg());
156  dnPHisto->Fill(dnpart->P4()->Vect().Mag());
157  dnRemovalEnergyHisto->Fill(dnpart->RemovalEnergy());
158 
159  // histograms for decayed nucleon daughters
160  ndaughters = 0;
161 
162  while((p=dynamic_cast<GHepParticle *>(event_iter.Next())))
163  {
164  if (p->FirstMother() == dnIndex) {
165  if (ndaughters == 0) {
166  dndaughter0PdgHisto->Fill(p->Pdg());
167  dndaughter0PHisto->Fill(p->P4()->Vect().Mag());
168  } else if (ndaughters == 1) {
169  dndaughter1PdgHisto->Fill(p->Pdg());
170  dndaughter1PHisto->Fill(p->P4()->Vect().Mag());
171  } else if (ndaughters == 2) {
172  dndaughter2PdgHisto->Fill(p->Pdg());
173  dndaughter2PHisto->Fill(p->P4()->Vect().Mag());
174  }
175  ndaughters++;
176  }
177  }
178  dndaughtersNHisto->Fill(ndaughters);
179 
180  // histograms for stable final state particles
181  nfinalparticles = 0;
182 
183  event_iter.Reset();
184  while((p=dynamic_cast<GHepParticle *>(event_iter.Next())))
185  {
186  if (p->Status() == kIStStableFinalState ) {
187  finalparticlesPdgHisto->Fill(p->Pdg());
188  finalparticlesPHisto->Fill(p->P4()->Vect().Mag());
189  nfinalparticles++;
190  }
191  }
192  finalparticlesNHisto->Fill(nfinalparticles);
193 
194 
195  // clear current mc event record
196  mcrec->Clear();
197 
198  }//end loop over events
199 
200  // close input GHEP event file
201  file.Close();
202 
203  // write TTexts explicitly. No need to do that for histograms, which are written by default
204  decayName.Write();
205  targetName.Write();
206 
207  histofile->Write();
208  histofile->Close();
209 
210  LOG("testNucleonDecay", pNOTICE) << "Done!";
211 
212  return 0;
213 }
214 //___________________________________________________________________
215 void GetCommandLineArgs(int argc, char ** argv)
216 {
217  LOG("testNucleonDecay", pINFO) << "Parsing commad line arguments";
218 
219  CmdLnArgParser parser(argc,argv);
220 
221  // get GENIE event sample
222  if( parser.OptionExists('f') ) {
223  LOG("testNucleonDecay", pINFO)
224  << "Reading event sample filename";
225  gOptInpFilename = parser.ArgAsString('f');
226  } else {
227  LOG("testNucleonDecay", pFATAL)
228  << "Unspecified input filename - Exiting";
229  exit(1);
230  }
231 
232  // number of events to analyse
233  if( parser.OptionExists('n') ) {
234  LOG("testNucleonDecay", pINFO)
235  << "Reading number of events to analyze";
236  gOptNEvt = parser.ArgAsInt('n');
237  } else {
238  LOG("testNucleonDecay", pINFO)
239  << "Unspecified number of events to analyze - Use all";
240  gOptNEvt = -1;
241  }
242 }
243 //_________________________________________________________________________________
int gOptNEvt
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
MINOS-style ntuple record. Each such ntuple record holds a generated EventRecord object. Ntuples of this type are intended for feeding GENIE events into other applications (for example the GEANT4 based MC generation framework of an experiment) if no direct interface exists.
#define pFATAL
Definition: Messenger.h:56
double RemovalEnergy(void) const
Get removal energy.
Definition: GHepParticle.h:100
GHepStatus_t Status(void) const
Definition: GHepParticle.h:64
int main(int argc, char **argv)
Definition: gAtmoEvGen.cxx:327
int Pdg(void) const
Definition: GHepParticle.h:63
int FirstMother(void) const
Definition: GHepParticle.h:66
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
MINOS-style Ntuple Class to hold an output MC Tree Header.
#define pINFO
Definition: Messenger.h:62
string gOptInpFilename
Definition: gEvDump.cxx:76
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:37
enum genie::ENucleonDecayMode NucleonDecayMode_t
string AsString(NucleonDecayMode_t ndm, int npdg=0)
Command line argument parser.
#define pNOTICE
Definition: Messenger.h:61
void GetCommandLineArgs(int argc, char **argv)
Definition: gAtmoEvGen.cxx:563
void Clear(Option_t *opt="")
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
EventRecord * event
event