gEvPick.cxx
Go to the documentation of this file.
1 //_____________________________________________________________________________________________
2 /*!
3 
4 \program gevpick
5 
6 \brief Reads a list of GENIE event files (GHEP format), `cherry-picks' events with a given
7  topology and writes them out in a separate file. The output event tree contains two
8  additional branches to aid book-keeping by maintaining a link to the source location
9  of each cherry-picked event. For each such event we store a) the name of the original
10  file and b) its original event number.
11 
12  This is the _only_recommended_ way to obtain event files that contain specific final
13  states (by cherry-picking events from files generated running GENIE in a comprehensive
14  mode). We don't recommend you attempt switching off generator-level reaction modes.
15  No detector measures generator-level reaction modes like CCQE or NCRES.
16  Detectors measure final states / topologies like {1mu-,0pi}, {1mu-,1pi+},
17  {0mu-, 1pi0}, {1 track, 1 shower}, {1 mu-like ring} etc depending on granularity,
18  thresholds and PID capabilities.
19  No final state / topology is a proxy for any particular reaction mode (and vice versa).
20  Intranuclear re-scattering in particular causes significant migration between states
21  (see Table 8.1 in the Physics and User manual).
22  Examples:
23  - {1mu-,0pi} is mostly numuCCQE but this particular final state can also come about
24  by numu resonance production followed by pion absorption.
25  - numuCCQE yields mostly {1mu-,0pi} final states but occasionaly can yield {1mu-,1pi}
26  if the recoil nucleon re-interacts.
27  - NC1pi0 final states can be caused by all
28  a) NC elastic followed by nucleon rescattering,
29  b) NC resonance neutrino-production,
30  c) NC non-resonance background,
31  d) low-W NC DIS
32  e) NC coherent scattering.
33  Each such NC1pi0 source contributes differently to the pion momentum distribution.
34 
35  Synopsis:
36  gevpick -i list_of_input_files -t topology
37  [-o output_file]
38  [--message-thresholds xmfile]
39  [--event-record-print-level level]
40 
41  Options:
42 
43  [] denotes an optional argument
44 
45  -i
46  Specify input file(s).
47  Wildcards accepted, eg `-i "/data/genie/t2k/gntp.*.ghep.root"'
48  -t
49  Specify event topology to cherry-pick.
50  The input topology can be any of
51  - all
52  all (basically merges all files into one)
53  - numu_cc_1pip
54  numu CC with 1 \pi^{+} (and no other pion) in final state
55  - numu_cc_1pi0
56  numu CC with 1 \pi^{0} (and no other pion) in final state
57  - numu_cc_1pim
58  numu CC with 1 \pi^{-} (and no other pion) in final state
59  - numu_nc_1pip
60  numu NC with 1 \pi^{+} (and no other pion) in final state
61  - numu_nc_1pi0
62  numu NC with 1 \pi^{0} (and no other pion) in final state
63  - numu_nc_1pim
64  numu NC with 1 \pi^{-} (and no other pion) in final state
65  - numu_cc_hyperon
66  numu CC with at least one hyperon
67  (\Sigma^{+,0,-}, \Lambda^{0}, \Xi^{0,-}, \Omega^{-}) in final state
68  - numubar_cc_hyperon
69  \bar{numu} CC with at least one hyperon
70  (\Sigma^{+,0,-}, \Lambda^{0}, \Xi^{0,-}, \Omega^{-}) in final state
71  - cc_hyperon
72  any (anti)neutrino CC with at least one hyperon
73  (\Sigma^{+,0,-}, \Lambda^{0}, \Xi^{0,-}, \Omega^{-}) in final state
74  - <can add more / please send request to costas.andreopoulos \at stfc.ac.uk>
75  -o
76  Specify output filename.
77  (optional, default: gntp.<topology>.ghep.root)
78  --message-thresholds
79  Allows users to customize the message stream thresholds.
80  The thresholds are specified using an XML file.
81  See $GENIE/config/Messenger.xml for the XML schema.
82  --event-record-print-level
83  Allows users to set the level of information shown when the event
84  record is printed in the screen. See GHepRecord::Print().
85 
86  Examples:
87 
88  (1) % gevpick -i "*.ghep.root" -t numu_nc_1pi0
89 
90  Will read all events in all *.ghep.root files and will cherry-pick
91  numu NC 1pi0 events. All cherry-picked events will be saved in the
92  output file gntp.numu_nc_1pi0.ghep.root (default name).
93 
94 \author Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
95  University of Liverpool & STFC Rutherford Appleton Lab
96 
97 \created August 09, 2010
98 
99 \cpright Copyright (c) 2003-2019, The GENIE Collaboration
100  For the full text of the license visit http://copyright.genie-mc.org
101  or see $GENIE/LICENSE
102 */
103 //_____________________________________________________________________________________________
104 
105 #include <cassert>
106 #include <string>
107 #include <sstream>
108 
109 #include <TSystem.h>
110 #include <TFile.h>
111 #include <TTree.h>
112 #include <TChain.h>
113 #include <TChainElement.h>
114 
128 #include "Framework/Utils/AppInit.h"
131 #include "Framework/Utils/RunOpt.h"
132 
133 using std::string;
134 using std::ostringstream;
135 
136 using namespace genie;
137 
138 // func prototypes
139 void GetCommandLineArgs (int argc, char ** argv);
140 void RunCherryPicker (void);
141 bool AcceptEvent (const EventRecord & event);
142 void PrintSyntax (void);
143 string DefaultOutputFile (void);
144 
145 // cherry-picked topologies
146 typedef enum EGPickTopo {
158 
159 } GPickTopo_t;
160 
161 // input options (from command line arguments):
162 string gOptInpFileNames; ///< input file name
163 string gOptOutFileName; ///< output file name
164 GPickTopo_t gPickedTopology; ///< output file format id
165 
166 //____________________________________________________________________________________
167 int main(int argc, char ** argv)
168 {
169  GetCommandLineArgs(argc, argv);
170 
171  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
172  GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
173 
174  RunCherryPicker();
175 
176  return 0;
177 }
178 //____________________________________________________________________________________
179 void RunCherryPicker(void)
180 {
181  // Create an NtpWriter for writing out a tree with the cherry-picked events
182  // Add 2 additional branches to the output event tree to save the original filename
183  // and the event number in the original file (so that all info can be traced back
184  // to its source).
185 
186  NtpWriter ntpw(kNFGHEP, 0);
188  ntpw.Initialize();
189  TObjString* brOrigFilename = new TObjString;
190  Long64_t brOrigEvtNum;
191  ntpw.EventTree()->Branch("orig_filename", "TObjString", &brOrigFilename, 5000,0);
192  ntpw.EventTree()->Branch("orig_evtnum", &brOrigEvtNum, "brOrigEvtNum/L");
193  Long64_t iev_glob = 0;
194 
195  // Load input trees. More than one trees can be loaded here if a wildcard was
196  // specified with -f (eg -f /data/myfiles/genie/*.ghep.root)
197 
198  TChain gchain;
199  gchain.Add(gOptInpFileNames.c_str());
200 
201  TObjArray * file_array = gchain.GetListOfFiles();
202  int nfiles = file_array->GetEntries();
203  LOG("gevpick", pFATAL)
204  << "Processing " << nfiles
205  << (nfiles==1 ? " file " : " files ");
206 
207  //
208  // Loop over input event files
209  //
210 
211  TIter next_file(file_array);
212  TChainElement *chEl=0;
213 
214  while (( chEl=(TChainElement*)next_file() )) {
215 
216  TFile fin(chEl->GetTitle(),"read");
217  TTree * ghep_tree =
218  dynamic_cast <TTree *> ( fin.Get("gtree") );
219 
220  if(!ghep_tree) {
221  LOG("gevpick", pWARN)
222  << "No GHEP tree found in " << chEl->GetTitle();
223  LOG("gevpick", pWARN)
224  << "Skipping to next file...";
225  continue;
226  }
227 
228  NtpMCEventRecord * mcrec = 0;
229  ghep_tree->SetBranchAddress("gmcrec", &mcrec);
230  if (!mcrec) {
231  LOG("gevpick", pERROR) << "Null MC record";
232  return;
233  }
234  Long64_t nmax = ghep_tree->GetEntries();
235  LOG("gevpick", pNOTICE)
236  << "* Analyzing: " << nmax
237  << " events from GHEP tree in file: " << chEl->GetTitle();
238 
239  NtpMCTreeHeader * thdr =
240  dynamic_cast <NtpMCTreeHeader *> ( fin.Get("header") );
241  LOG("gevpick", pNOTICE)
242  << "Input tree header: " << *thdr;
243 
244  //
245  // Loop over events in current file
246  //
247 
248  for(Long64_t iev = 0; iev < nmax; iev++) {
249  ghep_tree->GetEntry(iev);
250  NtpMCRecHeader rec_header = mcrec->hdr;
251  EventRecord & event = *(mcrec->event);
252  LOG("gevpick", pDEBUG) << rec_header;
253  LOG("gevpick", pDEBUG) << event;
254  if(AcceptEvent(event)) {
255  brOrigFilename->SetString(chEl->GetTitle());
256  brOrigEvtNum = iev;
257  EventRecord * event_copy = new EventRecord(event);
258  ntpw.AddEventRecord(iev_glob,event_copy);
259  iev_glob++;
260  }
261  mcrec->Clear();
262 
263  } // event loop (current file)
264  }// file loop
265 
266  // save the cherry-picked MC events
267  ntpw.Save();
268 
269  LOG("gevpick", pFATAL) << "Done!";
270 }
271 //____________________________________________________________________________________
272 bool AcceptEvent(const EventRecord & event)
273 {
274  if ( gPickedTopology == kPtAll ) return true;
275  if ( gPickedTopology == kPtUndefined ) return false;
276 
277  const Interaction * interaction = event.Summary();
278 
279  int nupdg = event.Probe()->Pdg();
280  bool isnumu = (nupdg == kPdgNuMu);
281  bool isnumubar = (nupdg == kPdgAntiNuMu);
282  bool iscc = interaction->ProcInfo().IsWeakCC();
283  bool isnc = interaction->ProcInfo().IsWeakNC();
284 
285  int NfP = 0; // number of protons in final state
286  int NfPbar = 0; // number of anti-protons in final state
287  int NfN = 0; // number of neutrons in final state
288  int NfNbar = 0; // number of anti-neutrons in final state
289  int NfPip = 0; // number of \pi^+'s in final state
290  int NfPim = 0; // number of \pi^-'s in final state
291  int NfPi0 = 0; // number of \pi^0's in final state
292  int NfKp = 0; // number of \K^+'s in final state
293  int NfKm = 0; // number of \K^-'s in final state
294  int NfK0 = 0; // number of \K^0's in final state
295  int NfK0bar = 0; // number of \bar{\K^0}'s in final state
296  int NfSigmap = 0; // number of \Sigma^+'s in final state
297  int NfSigma0 = 0; // number of \Sigma^0's in final state
298  int NfSigmam = 0; // number of \Sigma^-'s in final state
299  int NfLambda0 = 0; // number of \Lambda^0's in final state
300  int NfXi0 = 0; // number of \Xi^0's in final state
301  int NfXim = 0; // number of \Xi^-'s in final state
302  int NfOmegam = 0; // number of \Omega^-'s in final state
303  int NfOther = 0; // number of other particles in final state
304 
305  TObjArrayIter piter(&event);
306  GHepParticle * p = 0;
307  int ip=-1;
308  while( (p = (GHepParticle *) piter.Next())) {
309  ip++;
310  int pdgc = p->Pdg();
311  int ist = p->Status();
312  // only final state particles
313  if(ist!=kIStStableFinalState) continue;
314  // don't count final state lepton as part of the hadronic system
315  if(event.Particle(ip)->FirstMother()==0) continue;
316  // skip pseudo-particles
317  if(pdg::IsPseudoParticle(pdgc)) continue;
318  // count ...
319  if (pdgc == kPdgProton ) NfP++;
320  else if (pdgc == kPdgAntiProton ) NfPbar++;
321  else if (pdgc == kPdgNeutron ) NfN++;
322  else if (pdgc == kPdgAntiNeutron) NfNbar++;
323  else if (pdgc == kPdgPiP ) NfPip++;
324  else if (pdgc == kPdgPiM ) NfPim++;
325  else if (pdgc == kPdgPi0 ) NfPi0++;
326  else if (pdgc == kPdgKP ) NfKp++;
327  else if (pdgc == kPdgKM ) NfKm++;
328  else if (pdgc == kPdgK0 ) NfK0++;
329  else if (pdgc == kPdgAntiK0 ) NfK0bar++;
330  else if (pdgc == kPdgSigmaP ) NfSigmap++;
331  else if (pdgc == kPdgSigma0 ) NfSigma0++;
332  else if (pdgc == kPdgSigmaM ) NfSigmam++;
333  else if (pdgc == kPdgLambda ) NfLambda0++;
334  else if (pdgc == kPdgXi0 ) NfXi0++;
335  else if (pdgc == kPdgXiM ) NfXim++;
336  else if (pdgc == kPdgOmegaM ) NfOmegam++;
337  else NfOther++;
338  }
339 
340  bool is1pipX = (NfPip==1 && NfPi0==0 && NfPim==0);
341  bool is1pi0X = (NfPip==0 && NfPi0==1 && NfPim==0);
342  bool is1pimX = (NfPip==0 && NfPi0==0 && NfPim==1);
343  bool has_hype = (NfSigmap+NfSigma0+NfSigmam+NfLambda0+NfXi0+NfXim+NfOmegam > 0);
344 
345  if ( gPickedTopology == kPtNumuCC1pip ) {
346  if(isnumu && iscc && is1pipX) return true;
347  }
348  else
349  if ( gPickedTopology == kPtNumuCC1pi0 ) {
350  if(isnumu && iscc && is1pi0X) return true;
351  }
352  else
353  if ( gPickedTopology == kPtNumuCC1pim ) {
354  if(isnumu && iscc && is1pimX) return true;
355  }
356  else
357  if ( gPickedTopology == kPtNumuNC1pip ) {
358  if(isnumu && isnc && is1pipX) return true;
359  }
360  else
361  if ( gPickedTopology == kPtNumuNC1pi0 ) {
362  if(isnumu && isnc && is1pi0X) return true;
363  }
364  else
365  if ( gPickedTopology == kPtNumuNC1pim ) {
366  if(isnumu && isnc && is1pimX) return true;
367  }
368  else
370  if(isnumu && iscc && has_hype) return true;
371  }
372  else
374  if(isnumubar && iscc && has_hype) return true;
375  }
376  else
377  if ( gPickedTopology == kPtCChyperon ) {
378  if(iscc && has_hype) return true;
379  }
380 
381  return false;
382 }
383 //____________________________________________________________________________________
384 void GetCommandLineArgs(int argc, char ** argv)
385 {
386  // Common run options. Set defaults and read.
388 
389  // Parse run options for this app
390 
391  CmdLnArgParser parser(argc,argv);
392 
393  // get input ROOT file (containing a GENIE GHEP event tree)
394  if( parser.OptionExists('i') ) {
395  gOptInpFileNames = parser.ArgAsString('i');
396  } else {
397  LOG("gevpick", pFATAL)
398  << "Unspecified input filename - Exiting";
399  PrintSyntax();
400  gAbortingInErr = true;
401  exit(1);
402  }
403 
404  // requested topology
405  string topo = "";
406  if( parser.OptionExists('t') ) {
407  topo = parser.ArgAsString('t');
408  if ( topo == "all" ) { gPickedTopology = kPtAll; }
409  else if ( topo == "numu_cc_1pip" ) { gPickedTopology = kPtNumuCC1pip; }
410  else if ( topo == "numu_cc_1pi0" ) { gPickedTopology = kPtNumuCC1pi0; }
411  else if ( topo == "numu_cc_1pim" ) { gPickedTopology = kPtNumuCC1pim; }
412  else if ( topo == "numu_nc_1pip" ) { gPickedTopology = kPtNumuNC1pip; }
413  else if ( topo == "numu_nc_1pi0" ) { gPickedTopology = kPtNumuNC1pi0; }
414  else if ( topo == "numu_nc_1pim" ) { gPickedTopology = kPtNumuNC1pim; }
415  else if ( topo == "numu_cc_hyperon" ) { gPickedTopology = kPtNumuCChyperon; }
416  else if ( topo == "numubar_cc_hyperon" ) { gPickedTopology = kPtNumubarCChyperon; }
417  else if ( topo == "cc_hyperon" ) { gPickedTopology = kPtCChyperon; }
418  else { gPickedTopology = kPtUndefined; }
419 
421  LOG("gevpick", pFATAL) << "Unknown topology (" << topo << ")";
422  gAbortingInErr = true;
423  exit(1);
424  }
425 
426  } else {
427  LOG("gevpick", pFATAL) << "Unspecified event topology";
428  gAbortingInErr = true;
429  exit(1);
430  }
431 
432  // get output file name
433  if( parser.OptionExists('o') ) {
434  gOptOutFileName = parser.ArgAsString('o');
435  } else {
436  LOG("gevpick", pINFO)
437  << "Unspecified output filename - Using default";
439  }
440 
441  // Summarize
442  LOG("gevpick", pNOTICE)
443  << "\n\n gevpick job info: "
444  << "\n - input file(s) : " << gOptInpFileNames
445  << "\n - output file : " << gOptOutFileName
446  << "\n - cherry-picked topology : " << topo
447  << "\n";
448 }
449 //____________________________________________________________________________________
450 string DefaultOutputFile(void)
451 {
452  string tp = "";
453 
454  if (gPickedTopology == kPtAll ) { tp = "all"; }
455  else if (gPickedTopology == kPtNumuCC1pip ) { tp = "numu_cc_1pip"; }
456  else if (gPickedTopology == kPtNumuCC1pi0 ) { tp = "numu_cc_1pi0"; }
457  else if (gPickedTopology == kPtNumuCC1pim ) { tp = "numu_cc_1pim"; }
458  else if (gPickedTopology == kPtNumuNC1pip ) { tp = "numu_nc_1pip"; }
459  else if (gPickedTopology == kPtNumuNC1pi0 ) { tp = "numu_nc_1pi0"; }
460  else if (gPickedTopology == kPtNumuNC1pim ) { tp = "numu_nc_1pim"; }
461  else if (gPickedTopology == kPtNumuCChyperon ) { tp = "numu_cc_hyperon"; }
462  else if (gPickedTopology == kPtNumubarCChyperon ) { tp = "numubar_cc_hyperon"; }
463  else if (gPickedTopology == kPtCChyperon ) { tp = "cc_hyperon"; }
464 
465  ostringstream fnm;
466  fnm << "gntp." << tp << ".ghep.root";
467 
468  return fnm.str();
469 }
470 //____________________________________________________________________________________
471 void PrintSyntax(void)
472 {
473  string basedir = string( gSystem->Getenv("GENIE") );
474  string thisfile = basedir + string("/src/Apps/gEvPick.cxx");
475  string cmd = "less " + thisfile;
476 
477  gSystem->Exec(cmd.c_str());
478 }
479 //____________________________________________________________________________________
int iev
Definition: runWimpSim.h:118
static void SetPrintLevel(int print_level)
Definition: GHepRecord.cxx:992
TString fin
Definition: Style.C:24
virtual GHepParticle * Particle(int position) const
Definition: GHepRecord.cxx:148
void GetCommandLineArgs(int argc, char **argv)
Definition: gEvPick.cxx:384
const int kPdgXi0
Definition: PDGCodes.h:77
bool IsWeakCC(void) const
NtpMCRecHeader hdr
record header
Definition: NtpMCRecordI.h:39
string ArgAsString(char opt)
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
const int kPdgLambda
Definition: PDGCodes.h:69
void CustomizeFilename(string filename)
Definition: NtpWriter.cxx:118
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.
void ReadFromCommandLine(int argc, char **argv)
Definition: RunOpt.cxx:107
const char * p
Definition: xmltok.h:285
#define pFATAL
Definition: Messenger.h:57
const int kPdgNuMu
Definition: PDGCodes.h:30
const int kPdgSigma0
Definition: PDGCodes.h:72
TString ip
Definition: loadincs.C:5
enum EGPickTopo GPickTopo_t
GHepStatus_t Status(void) const
Definition: GHepParticle.h:65
const int kPdgK0
Definition: PDGCodes.h:151
int Pdg(void) const
Definition: GHepParticle.h:64
int FirstMother(void) const
Definition: GHepParticle.h:67
EGPickTopo
Definition: gEvPick.cxx:146
Summary information for an interaction.
Definition: Interaction.h:56
int iscc
GPickTopo_t gPickedTopology
output file format id
Definition: gEvPick.cxx:164
string cmd
Definition: run_hadd.py:52
bool IsWeakNC(void) const
TTree * EventTree(void)
Definition: NtpWriter.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
string gOptOutFileName
output file name
Definition: gEvPick.cxx:163
const int kPdgKM
Definition: PDGCodes.h:150
const int kPdgKP
Definition: PDGCodes.h:149
MINOS-style Ntuple Class to hold an output MC Tree Header.
const int kPdgPiP
Definition: PDGCodes.h:135
const int kPdgPi0
Definition: PDGCodes.h:137
void Save(void)
get the even tree
Definition: NtpWriter.cxx:214
#define pINFO
Definition: Messenger.h:63
int main(int argc, char **argv)
Definition: gEvPick.cxx:167
const int kPdgAntiK0
Definition: PDGCodes.h:152
const int kPdgOmegaM
Definition: PDGCodes.h:81
void AddEventRecord(int ievent, const EventRecord *ev_rec)
save the event tree
Definition: NtpWriter.cxx:69
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
#define pWARN
Definition: Messenger.h:61
void PrintSyntax(void)
Definition: gEvPick.cxx:471
const int kPdgSigmaM
Definition: PDGCodes.h:73
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:38
bool AcceptEvent(const EventRecord &event)
Definition: gEvPick.cxx:272
void Initialize(void)
add event
Definition: NtpWriter.cxx:95
static RunOpt * Instance(void)
Definition: RunOpt.cxx:62
const int kPdgXiM
Definition: PDGCodes.h:78
A utility class to facilitate creating the GENIE MC Ntuple from the output GENIE GHEP event records...
Definition: NtpWriter.h:40
void RunCherryPicker(void)
Definition: gEvPick.cxx:179
const int kPdgAntiNeutron
Definition: PDGCodes.h:68
bool IsPseudoParticle(int pdgc)
Definition: PDGUtils.cxx:25
exit(0)
const int kPdgAntiProton
Definition: PDGCodes.h:66
const int kPdgPiM
Definition: PDGCodes.h:136
string gOptInpFileNames
input file name
Definition: gEvPick.cxx:162
const int kPdgSigmaP
Definition: PDGCodes.h:71
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
void MesgThresholds(string inpfile)
Definition: AppInit.cxx:100
const int kPdgProton
Definition: PDGCodes.h:65
MINOS-style Ntuple Class to hold an MC Event Record Header.
Command line argument parser.
#define pNOTICE
Definition: Messenger.h:62
void Clear(Option_t *opt="")
string DefaultOutputFile(void)
Definition: gEvPick.cxx:450
bool gAbortingInErr
Definition: Messenger.cxx:56
const int kPdgNeutron
Definition: PDGCodes.h:67
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
bool OptionExists(char opt)
was option set?
EventRecord * event
event
#define pDEBUG
Definition: Messenger.h:64
enum BeamMode string