gEvPick.cxx File Reference

#include <cassert>
#include <string>
#include <sstream>
#include <TSystem.h>
#include <TFile.h>
#include <TTree.h>
#include <TChain.h>
#include <TChainElement.h>
#include "Framework/Conventions/GBuild.h"
#include "Framework/EventGen/EventRecord.h"
#include "Framework/GHEP/GHepStatus.h"
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/GHEP/GHepUtils.h"
#include "Framework/Ntuple/NtpMCFormat.h"
#include "Framework/Ntuple/NtpMCTreeHeader.h"
#include "Framework/Ntuple/NtpMCEventRecord.h"
#include "Framework/Ntuple/NtpWriter.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/Utils/AppInit.h"
#include "Framework/Utils/CmdLnArgParser.h"
#include "Framework/Utils/SystemUtils.h"
#include "Framework/Utils/RunOpt.h"

Go to the source code of this file.

Typedefs
typedef enum EGPickType	GPickType_t

Enumerations
enum	EGPickType {   kPtUndefined = 0, kPtAll, kPtTopoNumuCC1pip, kPtTopoNumuCC1pi0,   kPtTopoNumuCC1pim, kPtTopoNumuNC1pip, kPtTopoNumuNC1pi0, kPtTopoNumuNC1pim,   kPtTopoNumuCChyperon, kPtTopoNumubarCChyperon, kPtTopoCChyperon, kPtReacModeCCQE,   kPtReacModeNumuCCQE, kPtReacModeCCMEC, kPtReacModeNumuCCMEC, kPtReacModeCCQEMEC,   kPtReacModeNumuCCQEMEC, kPtReacModeNotCCQEMEC, kPtReacModeCCNotQEMEC, kPtReacModeNotNumuCCQEMEC,   kPtReacModeNumuNotCCQEMEC }

Functions
void	GetCommandLineArgs (int argc, char **argv)
void	RunCherryPicker (void)
bool	AcceptEvent (const EventRecord &event)
void	PrintSyntax (void)
string	DefaultOutputFile (void)
int	main (int argc, char **argv)

Variables
string	gOptInpFileNames
	input file name More...
string	gOptOutFileName
	output file name More...
string	gPickedTypeStr
	output file name More...
GPickType_t	gPickedType
	output file format id More...

Typedef Documentation

typedef enum EGPickType GPickType_t

Enumeration Type Documentation

enum EGPickType

Enumerator
kPtUndefined
kPtAll
kPtTopoNumuCC1pip
kPtTopoNumuCC1pi0
kPtTopoNumuCC1pim
kPtTopoNumuNC1pip
kPtTopoNumuNC1pi0
kPtTopoNumuNC1pim
kPtTopoNumuCChyperon
kPtTopoNumubarCChyperon
kPtTopoCChyperon
kPtReacModeCCQE
kPtReacModeNumuCCQE
kPtReacModeCCMEC
kPtReacModeNumuCCMEC
kPtReacModeCCQEMEC
kPtReacModeNumuCCQEMEC
kPtReacModeNotCCQEMEC
kPtReacModeCCNotQEMEC
kPtReacModeNotNumuCCQEMEC
kPtReacModeNumuNotCCQEMEC

Definition at line 183 of file gEvPick.cxx.

                        {
  kPtUndefined = 0,  
  kPtAll,
  //
  // Types based on final state topology, regardless of reaction mode
  //
  kPtTopoNumuCC1pip,
  kPtTopoNumuCC1pi0,
  kPtTopoNumuCC1pim,
  kPtTopoNumuNC1pip,
  kPtTopoNumuNC1pi0,
  kPtTopoNumuNC1pim,
  kPtTopoNumuCChyperon,
  kPtTopoNumubarCChyperon,
  kPtTopoCChyperon,
  //
  // Types based on true reaction mode, regardless of final state topology
  //
  kPtReacModeCCQE,
  kPtReacModeNumuCCQE,
  kPtReacModeCCMEC,
  kPtReacModeNumuCCMEC,
  kPtReacModeCCQEMEC,
  kPtReacModeNumuCCQEMEC,
  kPtReacModeNotCCQEMEC,
  kPtReacModeCCNotQEMEC,
  kPtReacModeNotNumuCCQEMEC,
  kPtReacModeNumuNotCCQEMEC

} GPickType_t;

Function Documentation

bool AcceptEvent

(

const EventRecord &

event

)

Definition at line 331 of file gEvPick.cxx.

{
  if ( gPickedType == kPtAll       ) return true;
  if ( gPickedType == kPtUndefined ) return false;

  const Interaction * interaction = event.Summary();

  int  nupdg     = event.Probe()->Pdg();
  bool isnumu    = (nupdg == kPdgNuMu);
  bool isnumubar = (nupdg == kPdgAntiNuMu);
  bool iscc      = interaction->ProcInfo().IsWeakCC();
  bool isnc      = interaction->ProcInfo().IsWeakNC();
  bool isqe      = interaction->ProcInfo().IsQuasiElastic();
  bool ismec     = interaction->ProcInfo().IsMEC();
  bool isstr     = interaction->ExclTag().IsStrangeEvent();
  bool ischm     = interaction->ExclTag().IsCharmEvent();

  int NfP        = 0; // number of protons         in final state
  int NfPbar     = 0; // number of anti-protons    in final state
  int NfN        = 0; // number of neutrons        in final state
  int NfNbar     = 0; // number of anti-neutrons   in final state
  int NfPip      = 0; // number of \pi^+'s         in final state
  int NfPim      = 0; // number of \pi^-'s         in final state
  int NfPi0      = 0; // number of \pi^0's         in final state
  int NfKp       = 0; // number of \K^+'s          in final state
  int NfKm       = 0; // number of \K^-'s          in final state
  int NfK0       = 0; // number of \K^0's          in final state
  int NfK0bar    = 0; // number of \bar{\K^0}'s    in final state
  int NfSigmap   = 0; // number of \Sigma^+'s      in final state
  int NfSigma0   = 0; // number of \Sigma^0's      in final state
  int NfSigmam   = 0; // number of \Sigma^-'s      in final state
  int NfLambda0  = 0; // number of \Lambda^0's     in final state
  int NfXi0      = 0; // number of \Xi^0's         in final state
  int NfXim      = 0; // number of \Xi^-'s         in final state
  int NfOmegam   = 0; // number of \Omega^-'s      in final state
  int NfOther    = 0; // number of other particles in final state

  TObjArrayIter piter(&event);
  GHepParticle * p = 0;
  int ip=-1;
  while( (p = (GHepParticle *) piter.Next())) {
    ip++;
    int pdgc = p->Pdg();
    int ist  = p->Status();
    // only final state particles
    if(ist!=kIStStableFinalState) continue;
    // don't count final state lepton as part of the hadronic system
    if(event.Particle(ip)->FirstMother()==0) continue;
    // skip pseudo-particles
    if(pdg::IsPseudoParticle(pdgc)) continue;
    // count ...
    if      (pdgc == kPdgProton     ) NfP++;
    else if (pdgc == kPdgAntiProton ) NfPbar++;
    else if (pdgc == kPdgNeutron    ) NfN++;
    else if (pdgc == kPdgAntiNeutron) NfNbar++;
    else if (pdgc == kPdgPiP        ) NfPip++;
    else if (pdgc == kPdgPiM        ) NfPim++;
    else if (pdgc == kPdgPi0        ) NfPi0++;
    else if (pdgc == kPdgKP         ) NfKp++;
    else if (pdgc == kPdgKM         ) NfKm++;
    else if (pdgc == kPdgK0         ) NfK0++;
    else if (pdgc == kPdgAntiK0     ) NfK0bar++;
    else if (pdgc == kPdgSigmaP     ) NfSigmap++;
    else if (pdgc == kPdgSigma0     ) NfSigma0++;
    else if (pdgc == kPdgSigmaM     ) NfSigmam++;
    else if (pdgc == kPdgLambda     ) NfLambda0++;
    else if (pdgc == kPdgXi0        ) NfXi0++;
    else if (pdgc == kPdgXiM        ) NfXim++;
    else if (pdgc == kPdgOmegaM     ) NfOmegam++;
    else                              NfOther++;
  }

  bool is1pipX  = (NfPip==1 && NfPi0==0 && NfPim==0);
  bool is1pi0X  = (NfPip==0 && NfPi0==1 && NfPim==0);
  bool is1pimX  = (NfPip==0 && NfPi0==0 && NfPim==1);
  bool has_hype = (NfSigmap+NfSigma0+NfSigmam+NfLambda0+NfXi0+NfXim+NfOmegam > 0);

  if ( gPickedType == kPtTopoNumuCC1pip ) {
    if(isnumu && iscc && is1pipX) return true;
  }
  else
  if ( gPickedType == kPtTopoNumuCC1pi0 ) {
    if(isnumu && iscc && is1pi0X) return true;
  }
  else
  if ( gPickedType == kPtTopoNumuCC1pim ) {
    if(isnumu && iscc && is1pimX) return true;
  }
  else
  if ( gPickedType == kPtTopoNumuNC1pip ) {
    if(isnumu && isnc && is1pipX) return true;
  }
  else
  if ( gPickedType == kPtTopoNumuNC1pi0 ) {
    if(isnumu && isnc && is1pi0X) return true;
  }
  else
  if ( gPickedType == kPtTopoNumuNC1pim ) {
    if(isnumu && isnc && is1pimX) return true;
  }
  else
  if ( gPickedType == kPtTopoNumuCChyperon ) {
    if(isnumu && iscc && has_hype) return true;
  }
  else
  if ( gPickedType == kPtTopoNumubarCChyperon ) {
    if(isnumubar && iscc && has_hype) return true;
  }
  else
  if ( gPickedType == kPtTopoCChyperon ) {
    if(iscc && has_hype) return true;
  }
  else 
  if ( gPickedType == kPtReacModeCCQE ) {
    if(isstr || ischm) return false;
    if(iscc && isqe) return true;
  }
  else 
  if ( gPickedType == kPtReacModeNumuCCQE ) {
    if(isstr || ischm) return false;
    if(isnumu && iscc && isqe) return true;
  }
  else 
  if ( gPickedType == kPtReacModeCCMEC ) {
    if(isstr || ischm) return false;
    if(iscc && ismec) return true;
  }
  else 
  if ( gPickedType == kPtReacModeNumuCCMEC ) {
    if(isstr || ischm) return false;
    if(isnumu && iscc && ismec) return true;
  }
  else 
  if ( gPickedType == kPtReacModeCCQEMEC ) {
    if(isstr || ischm) return false;
    if(iscc && (isqe || ismec)) return true;
  }
  else 
  if ( gPickedType == kPtReacModeNumuCCQEMEC ) {
    if(isstr || ischm) return false;
    if(isnumu && iscc && (isqe || ismec)) return true;
  }
  else 
  if ( gPickedType == kPtReacModeNotCCQEMEC ) {
    if(isstr || ischm) return false;
    if(!(iscc && (isqe || ismec))) return true;
  }
  else 
  if ( gPickedType == kPtReacModeNotNumuCCQEMEC ) {
    if(isstr || ischm) return false;
    if(!(isnumu && iscc && (isqe || ismec))) return true;
  }
  else 
  if ( gPickedType == kPtReacModeCCNotQEMEC ) {
    if(isstr || ischm) return false;
    if(iscc && !(isqe || ismec)) return true;
  }
  else 
  if ( gPickedType == kPtReacModeNumuNotCCQEMEC ) {
    if(isstr || ischm) return false;
    if(isnumu && !(iscc && (isqe || ismec))) return true;
  }

  return false;
}

string DefaultOutputFile

(

void

)

Definition at line 577 of file gEvPick.cxx.

{
  string tp = "";

  if      (gPickedType == kPtAll                    ) { tp = "all";                }

  else if (gPickedType == kPtTopoNumuCC1pip         ) { tp = "numu_cc_1pip";       }
  else if (gPickedType == kPtTopoNumuCC1pi0         ) { tp = "numu_cc_1pi0";       }
  else if (gPickedType == kPtTopoNumuCC1pim         ) { tp = "numu_cc_1pim";       }
  else if (gPickedType == kPtTopoNumuNC1pip         ) { tp = "numu_nc_1pip";       }
  else if (gPickedType == kPtTopoNumuNC1pi0         ) { tp = "numu_nc_1pi0";       }
  else if (gPickedType == kPtTopoNumuNC1pim         ) { tp = "numu_nc_1pim";       }
  else if (gPickedType == kPtTopoNumuCChyperon      ) { tp = "numu_cc_hyperon";    }
  else if (gPickedType == kPtTopoNumubarCChyperon   ) { tp = "numubar_cc_hyperon"; }
  else if (gPickedType == kPtTopoCChyperon          ) { tp = "cc_hyperon";         } 

  else if (gPickedType == kPtReacModeCCQE           ) { tp = "cc_qe";              }
  else if (gPickedType == kPtReacModeNumuCCQE       ) { tp = "numu_cc_qe";         }
  else if (gPickedType == kPtReacModeCCMEC          ) { tp = "cc_mec";             }
  else if (gPickedType == kPtReacModeNumuCCMEC      ) { tp = "numu_cc_mec";        }
  else if (gPickedType == kPtReacModeCCQEMEC        ) { tp = "cc_qe_mec";          }
  else if (gPickedType == kPtReacModeNumuCCQEMEC    ) { tp = "numu_cc_qe_mec";     }
  else if (gPickedType == kPtReacModeNotCCQEMEC     ) { tp = "not_cc_qe_mec";      }
  else if (gPickedType == kPtReacModeCCNotQEMEC     ) { tp = "cc_not_qe_mec";      }
  else if (gPickedType == kPtReacModeNotNumuCCQEMEC ) { tp = "not_numu_cc_qe_mec"; }
  else if (gPickedType == kPtReacModeNumuNotCCQEMEC ) { tp = "numu_not_cc_qe_mec"; }

  ostringstream fnm;
  fnm << "gntp." << tp << ".ghep.root";

  return fnm.str();
}

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

Definition at line 497 of file gEvPick.cxx.

{
  // Common run options. Set defaults and read.
  RunOpt::Instance()->ReadFromCommandLine(argc,argv);

  // Parse run options for this app

  CmdLnArgParser parser(argc,argv);

  // get input ROOT file (containing a GENIE GHEP event tree)
  if( parser.OptionExists('i') ) {
    gOptInpFileNames = parser.ArgAsString('i');
  } else {
    LOG("gevpick", pFATAL)
       << "Unspecified input filename - Exiting";
    PrintSyntax();
    gAbortingInErr = true;
    exit(1);
  }

  // requested event type
  string evtype = "";
  if( parser.OptionExists('t') ) {
    evtype = parser.ArgAsString('t');
    if      ( evtype == "all"                ) { gPickedType = kPtAll;                    }

    else if ( evtype == "numu_cc_1pip"       ) { gPickedType = kPtTopoNumuCC1pip;         }
    else if ( evtype == "numu_cc_1pi0"       ) { gPickedType = kPtTopoNumuCC1pi0;         }
    else if ( evtype == "numu_cc_1pim"       ) { gPickedType = kPtTopoNumuCC1pim;         }
    else if ( evtype == "numu_nc_1pip"       ) { gPickedType = kPtTopoNumuNC1pip;         }
    else if ( evtype == "numu_nc_1pi0"       ) { gPickedType = kPtTopoNumuNC1pi0;         }
    else if ( evtype == "numu_nc_1pim"       ) { gPickedType = kPtTopoNumuNC1pim;         }
    else if ( evtype == "numu_cc_hyperon"    ) { gPickedType = kPtTopoNumuCChyperon;      }
    else if ( evtype == "numubar_cc_hyperon" ) { gPickedType = kPtTopoNumubarCChyperon;   }
    else if ( evtype == "cc_hyperon"         ) { gPickedType = kPtTopoCChyperon;          }

    else if ( evtype == "cc_qe"              ) { gPickedType = kPtReacModeCCQE;           }
    else if ( evtype == "numu_cc_qe"         ) { gPickedType = kPtReacModeNumuCCQE;       }
    else if ( evtype == "cc_mec"             ) { gPickedType = kPtReacModeCCMEC;          }
    else if ( evtype == "numu_cc_mec"        ) { gPickedType = kPtReacModeNumuCCMEC;      }
    else if ( evtype == "cc_qe_mec"          ) { gPickedType = kPtReacModeCCQEMEC;        }
    else if ( evtype == "numu_cc_qe_mec"     ) { gPickedType = kPtReacModeNumuCCQEMEC;    }
    else if ( evtype == "not_cc_qe_mec"      ) { gPickedType = kPtReacModeNotCCQEMEC;     }
    else if ( evtype == "cc_not_qe_mec"      ) { gPickedType = kPtReacModeCCNotQEMEC;     }
    else if ( evtype == "not_numu_cc_qe_mec" ) { gPickedType = kPtReacModeNotNumuCCQEMEC; }
    else if ( evtype == "numu_not_cc_qe_mec" ) { gPickedType = kPtReacModeNumuNotCCQEMEC; }

    else                                       { gPickedType = kPtUndefined;              }

    if(gPickedType == kPtUndefined) {
      LOG("gevpick", pFATAL) << "Unknown event type (" << evtype << ")";
      gAbortingInErr = true;
      exit(1);
    }
    gPickedTypeStr = evtype;

  } else {
    LOG("gevpick", pFATAL) << "Unspecified event type";
    gAbortingInErr = true;
    exit(1);
  }

  // get output file name 
  if( parser.OptionExists('o') ) {
    gOptOutFileName = parser.ArgAsString('o');
  } else {
    LOG("gevpick", pINFO)
       << "Unspecified output filename - Using default";
    gOptOutFileName = DefaultOutputFile();
  }

  // Summarize
  LOG("gevpick", pNOTICE) 
    << "\n\n gevpick job info: "
    << "\n - input file(s)          : " << gOptInpFileNames
    << "\n - output file            : " << gOptOutFileName
    << "\n - cherry-picked topology : " << evtype
    << "\n";
}

int main	(	int	argc,
		char **	argv
	)

Definition at line 221 of file gEvPick.cxx.

{
  GetCommandLineArgs(argc, argv);

  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
  GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());

  RunCherryPicker();

  return 0;
}

void PrintSyntax

(

void

)

Definition at line 610 of file gEvPick.cxx.

{
  string basedir  = string( gSystem->Getenv("GENIE") );
  string thisfile = basedir + string("/src/Apps/gEvPick.cxx");
  string cmd      = "less " + thisfile;

  gSystem->Exec(cmd.c_str());
}

void RunCherryPicker

(

void

)

Definition at line 233 of file gEvPick.cxx.

{
  // Create an NtpWriter for writing out a tree with the cherry-picked events
  // Add 2 additional branches to the output event tree to save the original filename
  // and the event number in the original file (so that all info can be traced back 
  // to its source).

  NtpWriter ntpw(kNFGHEP, 0);
  ntpw.CustomizeFilename(gOptOutFileName);
  ntpw.Initialize();
  TObjString* brOrigFilename = new TObjString;
  Long64_t    brOrigEvtNum;
  ntpw.EventTree()->Branch("orig_filename", "TObjString", &brOrigFilename, 5000,0);
  ntpw.EventTree()->Branch("orig_evtnum", &brOrigEvtNum, "brOrigEvtNum/L");
  Long64_t iev_glob = 0;

  // Load input trees. More than one trees can be loaded here if a wildcard was
  // specified with -f (eg -f /data/myfiles/genie/*.ghep.root)

  TChain gchain;
  gchain.Add(gOptInpFileNames.c_str());

  TObjArray * file_array = gchain.GetListOfFiles();
  int nfiles = file_array->GetEntries();
  LOG("gevpick", pNOTICE) 
      << "Processing " << nfiles
      << (nfiles==1 ? " file " : " files ");

  //
  // Loop over input event files
  //

  TIter next_file(file_array);
  TChainElement *chEl=0;

  unsigned int total_events  = 0;
  unsigned int picked_events = 0;

  while (( chEl=(TChainElement*)next_file() )) {

     TFile fin(chEl->GetTitle(),"read");
     TTree * ghep_tree = 
        dynamic_cast <TTree *> ( fin.Get("gtree")  );

     if(!ghep_tree) {
        LOG("gevpick", pWARN) 
           << "No GHEP tree found in " << chEl->GetTitle();
        LOG("gevpick", pWARN) 
           << "Skipping to next file...";
        continue;
     }

     NtpMCEventRecord * mcrec = 0;
     ghep_tree->SetBranchAddress("gmcrec", &mcrec);
     if (!mcrec) {
       LOG("gevpick", pERROR) << "Null MC record";
       return;
     }
     Long64_t nmax = ghep_tree->GetEntries();
     LOG("gevpick", pNOTICE) 
        << "* Analyzing: " << nmax 
        << " events from GHEP tree in file: " << chEl->GetTitle();

     NtpMCTreeHeader * thdr = 
        dynamic_cast <NtpMCTreeHeader *> ( fin.Get("header") );
     LOG("gevpick", pNOTICE) 
          << "Input tree header: " << *thdr;

     //
     // Loop over events in current file
     //

     for(Long64_t iev = 0; iev < nmax; iev++) {
       total_events++;
       ghep_tree->GetEntry(iev);
       NtpMCRecHeader rec_header = mcrec->hdr;
       EventRecord &  event      = *(mcrec->event);
       LOG("gevpick", pDEBUG) << rec_header;
       LOG("gevpick", pDEBUG) << event;
       if(AcceptEvent(event)) {
         picked_events++;
          brOrigFilename->SetString(chEl->GetTitle());
          brOrigEvtNum = iev;
          ntpw.AddEventRecord( iev_glob, &event );
          iev_glob++;
       }
       mcrec->Clear();

    } // event loop (current file)
  }// file loop

  // save the cherry-picked MC events
  ntpw.Save();
  
  LOG("gevpick", pNOTICE) << "Picked " << picked_events << " / " << total_events << " events of type " << gPickedTypeStr;
  LOG("gevpick", pNOTICE) << "Done!";
}

Variable Documentation

string gOptInpFileNames

input file name

Definition at line 215 of file gEvPick.cxx.

string gOptOutFileName

output file name

Definition at line 216 of file gEvPick.cxx.

GPickType_t gPickedType

output file format id

Definition at line 218 of file gEvPick.cxx.

string gPickedTypeStr

output file name

Definition at line 217 of file gEvPick.cxx.