gtestEventLoop.cxx

Go to the documentation of this file.

//____________________________________________________________________________
/*!
 
\program gtestEventLoop
 
\brief   Example event loop. Use that as a template for your analysis code.
 
\author  Costas Andreopoulos <c.andreopoulos \at cern.ch>
 University of Liverpool
 
\created May 4, 2004
 
\cpright Copyright (c) 2003-2025, The GENIE Collaboration
         For the full text of the license visit http://copyright.genie-mc.org
         
*/
//____________________________________________________________________________
 
#include <string>
 
#include <TSystem.h>
#include <TFile.h>
#include <TTree.h>
#include <TIterator.h>
 
#include "Framework/EventGen/EventRecord.h"
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/Ntuple/NtpMCFormat.h"
#include "Framework/Ntuple/NtpMCTreeHeader.h"
#include "Framework/Ntuple/NtpMCEventRecord.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/Utils/CmdLnArgParser.h"
 
using std::string;
using namespace genie;
 
void GetCommandLineArgs (int argc, char ** argv);
 
int    gOptNEvt;
string gOptInpFilename;
 
//___________________________________________________________________
int main(int argc, char ** argv)
{
  GetCommandLineArgs (argc, argv);
 
  //-- open the ROOT file and get the TTree & its header
  TTree *           tree = 0;
  NtpMCTreeHeader * thdr = 0;
 
  TFile file(gOptInpFilename.c_str(),"READ");
 
  tree = dynamic_cast <TTree *>           ( file.Get("gtree")  );
  thdr = dynamic_cast <NtpMCTreeHeader *> ( file.Get("header") );
 
  if(!tree) return 1;
 
  NtpMCEventRecord * mcrec = 0;
  tree->SetBranchAddress("gmcrec", &mcrec);
 
  int nev = (gOptNEvt > 0) ?
        TMath::Min(gOptNEvt, (int)tree->GetEntries()) :
        (int) tree->GetEntries();
 
  //
  // Loop over all events
  //
  for(int i = 0; i < nev; i++) {
 
    // get next tree entry
    tree->GetEntry(i);
 
    // get the GENIE event
    EventRecord &  event = *(mcrec->event);
 
    LOG("myAnalysis", pNOTICE) << event;
 
    //
    // Put your event analysis code here 
    //
    // ... ... ... ... ...
    // ... ... ... ... ...
    //
    // 
 
    
 
    //
    // Loop over all particles in this event
    //
 
    GHepParticle * p = 0;
    TIter event_iter(&event);
 
    while((p=dynamic_cast<GHepParticle *>(event_iter.Next())))
    {
       //
       // Put your event analysis code here 
       //
       // ... ... ... ... ...
       // ... ... ... ... ...
       //
       //
 
       // EXAMPLE: Print out the energy of all final state pions.
       if (p->Status() == kIStStableFinalState ) {
          if (p->Pdg() == kPdgPi0 ||
              p->Pdg() == kPdgPiP ||
              p->Pdg() == kPdgPiM) 
           {
            LOG("myAnalysis", pNOTICE)  
               << "Got a : " << p->Name() << " with E = " << p->E() << " GeV"; 
          }
       }
 
    }// end loop over particles 
    
    // clear current mc event record
    mcrec->Clear();
 
  }//end loop over events
 
  // close input GHEP event file
  file.Close();
 
  LOG("myAnalysis", pNOTICE)  << "Done!";
 
  return 0;
}
//___________________________________________________________________
void GetCommandLineArgs(int argc, char ** argv)
{
  LOG("myAnalysis", pINFO) << "Parsing commad line arguments";
 
  CmdLnArgParser parser(argc,argv);
 
  // get GENIE event sample
  if( parser.OptionExists('f') ) {
    LOG("myAnalysis", pINFO) 
       << "Reading event sample filename";
    gOptInpFilename = parser.ArgAsString('f');
  } else {
    LOG("myAnalysis", pFATAL) 
        << "Unspecified input filename - Exiting";
    exit(1);
  }
 
  // number of events to analyse
  if( parser.OptionExists('n') ) {
    LOG("myAnalysis", pINFO) 
      << "Reading number of events to analyze";
    gOptNEvt = parser.ArgAsInt('n');
  } else {
    LOG("myAnalysis", pINFO)
      << "Unspecified number of events to analyze - Use all";
    gOptNEvt = -1;
  }
}
//_________________________________________________________________________________