GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gtestEventLoop.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 
34 
35 using std::string;
36 using namespace genie;
37 
38 void GetCommandLineArgs (int argc, char ** argv);
39 
42 
43 //___________________________________________________________________
44 int main(int argc, char ** argv)
45 {
46  GetCommandLineArgs (argc, argv);
47 
48  //-- open the ROOT file and get the TTree & its header
49  TTree * tree = 0;
50  NtpMCTreeHeader * thdr = 0;
51 
52  TFile file(gOptInpFilename.c_str(),"READ");
53 
54  tree = dynamic_cast <TTree *> ( file.Get("gtree") );
55  thdr = dynamic_cast <NtpMCTreeHeader *> ( file.Get("header") );
56 
57  if(!tree) return 1;
58 
59  NtpMCEventRecord * mcrec = 0;
60  tree->SetBranchAddress("gmcrec", &mcrec);
61 
62  int nev = (gOptNEvt > 0) ?
63  TMath::Min(gOptNEvt, (int)tree->GetEntries()) :
64  (int) tree->GetEntries();
65 
66  //
67  // Loop over all events
68  //
69  for(int i = 0; i < nev; i++) {
70 
71  // get next tree entry
72  tree->GetEntry(i);
73 
74  // get the GENIE event
75  EventRecord & event = *(mcrec->event);
76 
77  LOG("myAnalysis", pNOTICE) << event;
78 
79  //
80  // Put your event analysis code here
81  //
82  // ... ... ... ... ...
83  // ... ... ... ... ...
84  //
85  //
86 
87 
88 
89  //
90  // Loop over all particles in this event
91  //
92 
93  GHepParticle * p = 0;
94  TIter event_iter(&event);
95 
96  while((p=dynamic_cast<GHepParticle *>(event_iter.Next())))
97  {
98  //
99  // Put your event analysis code here
100  //
101  // ... ... ... ... ...
102  // ... ... ... ... ...
103  //
104  //
105 
106  // EXAMPLE: Print out the energy of all final state pions.
107  if (p->Status() == kIStStableFinalState ) {
108  if (p->Pdg() == kPdgPi0 ||
109  p->Pdg() == kPdgPiP ||
110  p->Pdg() == kPdgPiM)
111  {
112  LOG("myAnalysis", pNOTICE)
113  << "Got a : " << p->Name() << " with E = " << p->E() << " GeV";
114  }
115  }
116 
117  }// end loop over particles
118 
119  // clear current mc event record
120  mcrec->Clear();
121 
122  }//end loop over events
123 
124  // close input GHEP event file
125  file.Close();
126 
127  LOG("myAnalysis", pNOTICE) << "Done!";
128 
129  return 0;
130 }
131 //___________________________________________________________________
132 void GetCommandLineArgs(int argc, char ** argv)
133 {
134  LOG("myAnalysis", pINFO) << "Parsing commad line arguments";
135 
136  CmdLnArgParser parser(argc,argv);
137 
138  // get GENIE event sample
139  if( parser.OptionExists('f') ) {
140  LOG("myAnalysis", pINFO)
141  << "Reading event sample filename";
142  gOptInpFilename = parser.ArgAsString('f');
143  } else {
144  LOG("myAnalysis", pFATAL)
145  << "Unspecified input filename - Exiting";
146  exit(1);
147  }
148 
149  // number of events to analyse
150  if( parser.OptionExists('n') ) {
151  LOG("myAnalysis", pINFO)
152  << "Reading number of events to analyze";
153  gOptNEvt = parser.ArgAsInt('n');
154  } else {
155  LOG("myAnalysis", pINFO)
156  << "Unspecified number of events to analyze - Use all";
157  gOptNEvt = -1;
158  }
159 }
160 //_________________________________________________________________________________
double E(void) const
Get energy.
Definition: GHepParticle.h:91
int gOptNEvt
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
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
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
MINOS-style Ntuple Class to hold an output MC Tree Header.
const int kPdgPiP
Definition: PDGCodes.h:158
const int kPdgPi0
Definition: PDGCodes.h:160
#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
const int kPdgPiM
Definition: PDGCodes.h:159
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