diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx
index 36b63d2..60ef26f 100644
--- a/src/InputHandler/GENIEInputHandler.cxx
+++ b/src/InputHandler/GENIEInputHandler.cxx
@@ -1,635 +1,636 @@
// Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifdef __GENIE_ENABLED__
#include "GENIEInputHandler.h"
#ifdef GENIE_PRE_R3
#include "Messenger/Messenger.h"
#else
#include "Framework/Messenger/Messenger.h"
#endif
#include "InputUtils.h"
GENIEGeneratorInfo::~GENIEGeneratorInfo() { DeallocateParticleStack(); }
void GENIEGeneratorInfo::AddBranchesToTree(TTree *tn) {
tn->Branch("GenieParticlePDGs", &fGenieParticlePDGs, "GenieParticlePDGs/I");
}
void GENIEGeneratorInfo::SetBranchesFromTree(TTree *tn) {
tn->SetBranchAddress("GenieParticlePDGs", &fGenieParticlePDGs);
}
void GENIEGeneratorInfo::AllocateParticleStack(int stacksize) {
fGenieParticlePDGs = new int[stacksize];
}
void GENIEGeneratorInfo::DeallocateParticleStack() {
delete fGenieParticlePDGs;
}
void GENIEGeneratorInfo::FillGeneratorInfo(NtpMCEventRecord *ntpl) {
Reset();
// Check for GENIE Event
if (!ntpl)
return;
if (!ntpl->event)
return;
// Cast Event Record
GHepRecord *ghep = static_cast<GHepRecord *>(ntpl->event);
if (!ghep)
return;
// Fill Particle Stack
GHepParticle *p = 0;
TObjArrayIter iter(ghep);
// Loop over all particles
int i = 0;
while ((p = (dynamic_cast<genie::GHepParticle *>((iter).Next())))) {
if (!p)
continue;
// Get PDG
fGenieParticlePDGs[i] = p->Pdg();
i++;
}
}
void GENIEGeneratorInfo::Reset() {
for (int i = 0; i < kMaxParticles; i++) {
fGenieParticlePDGs[i] = 0;
}
}
+bool GENIEInputHandler::IsPrimary(GHepParticle *p) {
+
+ // If initial state nucleon, or nucleon target, it is definintely primary!
+ if (p->Status() == genie::kIStInitialState ||
+ p->Status() == genie::kIStNucleonTarget) return true;
+
+ // Reject intermediate resonant state or pre-DIS state
+ if (p->Status() == genie::kIStDISPreFragmHadronicState || // DIS before fragmentation
+ p->Status() == genie::kIStPreDecayResonantState || // pre decay resonance state
+ p->Status() == genie::kIStIntermediateState || // intermediate state
+ p->Status() == genie::kIStDecayedState || // Decayed state
+ p->Status() == genie::kIStUndefined) return false; // undefined state
+
+ // Check if the mother is the neutrino or IS nucleon
+ if (p->FirstMother() < 2) return true;
+
+ // We've now filtered off intermediate states and obvious initial states
+
+ // Loop over this particle's mothers, grandmothers, and so on cleaning out intermediate particles
+ // Set the starting particle to be our particle
+ GHepParticle *mother = fGenieGHep->Particle(p->FirstMother());
+ while (mother->FirstMother() > 1) {
+ // It could be that the mother's status is actually a decayed state linked back to the vertex
+ if (mother->Status() == genie::kIStDecayedState || // A decayed state
+ mother->Status() == genie::kIStDISPreFragmHadronicState || // A DIS state before fragmentation
+ mother->Status() == genie::kIStPreDecayResonantState ) { // A pre-decay resonant state
+ mother = fGenieGHep->Particle(mother->FirstMother());
+ } else { // If not, move out of the loop
+ break;
+ }
+ }
+
+ // Then do a simple check of mother is associated with the primary
+ int MotherID = mother->FirstMother();
+ if (MotherID > 2) return false;
+
+ // Finally, this should mean that our partcile is marked for transport through the nucleus
+ // Could also be interactions of free proton
+ if (p->Status() == genie::kIStHadronInTheNucleus || // Then require the particle to be paseed to FSI
+ (p->Status() == genie::kIStStableFinalState && // Can also have interaction on free proton
+ fGenieGHep->Summary()->InitState().TgtPtr()->A() == 1 &&
+ fGenieGHep->Summary()->InitState().TgtPtr()->Z() == 1) ) {
+ return true;
+ }
+
+ return false;
+}
+
GENIEInputHandler::GENIEInputHandler(std::string const &handle,
- std::string const &rawinputs) {
+ std::string const &rawinputs) {
NUIS_LOG(SAM, "Creating GENIEInputHandler : " << handle);
// Plz no shouting
StopTalking();
genie::Messenger::Instance()->SetPriorityLevel("GHepUtils", pFATAL);
StartTalking();
// Shout all you want
// Run a joint input handling
fName = handle;
// Setup the TChain
fGENIETree = new TChain("gtree");
fSaveExtra = FitPar::Config().GetParB("SaveExtraGenie");
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
- // Are we running with NOvA weights
- fNOvAWeights = FitPar::Config().GetParB("NOvA_Weights");
- MAQEw = 1.0;
- NonResw = 1.0;
- RPAQEw = 1.0;
- RPARESw = 1.0;
- MECw = 1.0;
- DISw = 1.0;
- NOVAw = 1.0;
-
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector<std::string> inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
TFile *inp_file = new TFile(
InputUtils::ExpandInputDirectories(inputs[inp_it]).c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
NUIS_ABORT(
"GENIE File IsZombie() at : '"
<< inputs[inp_it] << "'" << std::endl
<< "Check that your file paths are correct and the file exists!"
<< std::endl
<< "$ ls -lh " << inputs[inp_it]);
}
// Get Flux/Event hist
TH1D *fluxhist = (TH1D *)inp_file->Get("nuisance_flux");
TH1D *eventhist = (TH1D *)inp_file->Get("nuisance_events");
if (!fluxhist or !eventhist) {
NUIS_ERR(FTL, "Input File Contents: " << inputs[inp_it]);
inp_file->ls();
NUIS_ABORT("GENIE FILE doesn't contain flux/xsec info."
- << std::endl
- << "Try running the app PrepareGENIE first on :"
- << inputs[inp_it] << std::endl
- << "$ PrepareGENIE -h");
+ << std::endl
+ << "Try running the app PrepareGENIE first on :"
+ << inputs[inp_it] << std::endl
+ << "$ PrepareGENIE -h");
}
// Get N Events
TTree *genietree = (TTree *)inp_file->Get("gtree");
if (!genietree) {
NUIS_ERR(FTL, "gtree not located in GENIE file: " << inputs[inp_it]);
NUIS_ABORT(
"Check your inputs, they may need to be completely regenerated!");
}
int nevents = genietree->GetEntries();
if (nevents <= 0) {
NUIS_ABORT("Trying to a TTree with "
- << nevents << " to TChain from : " << inputs[inp_it]);
- }
-
- // Check for precomputed weights
- TTree *weighttree = (TTree *)inp_file->Get("nova_wgts");
- if (fNOvAWeights) {
- if (!weighttree) {
- NUIS_ABORT("Did not find nova_wgts tree in file "
- << inputs[inp_it] << " but you specified it" << std::endl);
- } else {
- NUIS_LOG(FIT, "Found nova_wgts tree in file " << inputs[inp_it]);
- }
+ << nevents << " to TChain from : " << inputs[inp_it]);
}
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add To TChain
fGENIETree->AddFile(inputs[inp_it].c_str());
- if (weighttree != NULL)
- fGENIETree->AddFriend(weighttree);
}
// Registor all our file inputs
SetupJointInputs();
// Assign to tree
fEventType = kGENIE;
fGenieNtpl = NULL;
fGENIETree->SetBranchAddress("gmcrec", &fGenieNtpl);
- // Set up the custom weights
- if (fNOvAWeights) {
- fGENIETree->SetBranchAddress("MAQEwgt", &MAQEw);
- fGENIETree->SetBranchAddress("nonResNormWgt", &NonResw);
- fGENIETree->SetBranchAddress("RPAQEWgt", &RPAQEw);
- fGENIETree->SetBranchAddress("RPARESWgt", &RPARESw);
- fGENIETree->SetBranchAddress("MECWgt", &MECw);
- fGENIETree->SetBranchAddress("DISWgt", &DISw);
- fGENIETree->SetBranchAddress("nova2018CVWgt", &NOVAw);
- }
-
// Libraries should be seen but not heard...
StopTalking();
fGENIETree->GetEntry(0);
StartTalking();
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetGenieEvent(fGenieNtpl);
if (fSaveExtra) {
fGenieInfo = new GENIEGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fGenieInfo);
}
fNUISANCEEvent->HardReset();
};
GENIEInputHandler::~GENIEInputHandler() {
// if (fGenieGHep) delete fGenieGHep;
// if (fGenieNtpl) delete fGenieNtpl;
// if (fGENIETree) delete fGENIETree;
// if (fGenieInfo) delete fGenieInfo;
}
void GENIEInputHandler::CreateCache() {
if (fCacheSize > 0) {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 1);
fGENIETree->SetCacheSize(fCacheSize);
}
}
void GENIEInputHandler::RemoveCache() {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 0);
fGENIETree->SetCacheSize(0);
}
FitEvent *GENIEInputHandler::GetNuisanceEvent(const UInt_t ent,
- const bool lightweight) {
+ const bool lightweight) {
UInt_t entry = ent + fSkip;
if (entry >= (UInt_t)fNEvents)
return NULL;
// Clear the previous event (See Note 1 in ROOT TClonesArray documentation)
if (fGenieNtpl) {
fGenieNtpl->Clear();
}
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fGENIETree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
else {
// Check for GENIE Event
if (!fGenieNtpl)
return NULL;
if (!fGenieNtpl->event)
return NULL;
// Cast Event Record
fGenieGHep = static_cast<GHepRecord *>(fGenieNtpl->event);
if (!fGenieGHep)
return NULL;
TObjArrayIter iter(fGenieGHep);
genie::GHepParticle *p;
while ((p = (dynamic_cast<genie::GHepParticle *>((iter).Next())))) {
if (!p) {
continue;
}
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
if (state != genie::kIStInitialState) {
continue;
}
fNUISANCEEvent->probe_E = p->E() * 1.E3;
fNUISANCEEvent->probe_pdg = p->Pdg();
break;
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
return fNUISANCEEvent;
}
int GENIEInputHandler::GetGENIEParticleStatus(genie::GHepParticle *p,
- int mode) {
+ int mode) {
/*
kIStUndefined = -1,
kIStInitialState = 0, / generator-level initial state /
kIStStableFinalState = 1, / generator-level final state:
particles to be tracked by detector-level MC /
kIStIntermediateState = 2,
kIStDecayedState = 3,
kIStCorrelatedNucleon = 10,
kIStNucleonTarget = 11,
kIStDISPreFragmHadronicState = 12,
kIStPreDecayResonantState = 13,
kIStHadronInTheNucleus = 14, / hadrons inside the nucleus: marked
for hadron transport modules to act on /
kIStFinalStateNuclearRemnant = 15, / low energy nuclear fragments
entering the record collectively as a 'hadronic blob' pseudo-particle /
kIStNucleonClusterTarget = 16, // for composite nucleons before
phase space decay
*/
int state = kUndefinedState;
switch (p->Status()) {
- case genie::kIStNucleonTarget:
- case genie::kIStInitialState:
- case genie::kIStCorrelatedNucleon:
- case genie::kIStNucleonClusterTarget:
- state = kInitialState;
- break;
-
- case genie::kIStStableFinalState:
- state = kFinalState;
- break;
-
- case genie::kIStHadronInTheNucleus:
- if (abs(mode) == 2)
+ case genie::kIStNucleonTarget:
+ case genie::kIStInitialState:
+ case genie::kIStCorrelatedNucleon:
+ case genie::kIStNucleonClusterTarget:
state = kInitialState;
- else
+ break;
+
+ case genie::kIStStableFinalState:
+ state = kFinalState;
+ break;
+
+ case genie::kIStHadronInTheNucleus:
+ if (abs(mode) == 2)
+ state = kInitialState;
+ else
+ state = kFSIState;
+ break;
+
+ case genie::kIStPreDecayResonantState:
+ case genie::kIStDISPreFragmHadronicState:
+ case genie::kIStIntermediateState:
state = kFSIState;
- break;
-
- case genie::kIStPreDecayResonantState:
- case genie::kIStDISPreFragmHadronicState:
- case genie::kIStIntermediateState:
- state = kFSIState;
- break;
-
- case genie::kIStFinalStateNuclearRemnant:
- case genie::kIStUndefined:
- case genie::kIStDecayedState:
- default:
- break;
+ break;
+
+ case genie::kIStFinalStateNuclearRemnant:
+ case genie::kIStUndefined:
+ case genie::kIStDecayedState:
+ default:
+ break;
}
// Flag to remove nuclear part in genie
- if (p->Pdg() > 1000000) {
- if (state == kInitialState)
- state = kNuclearInitial;
- else if (state == kFinalState)
- state = kNuclearRemnant;
+ if (kRemoveNuclearParticles) {
+ if (p->Pdg() > 1000000 &&
+ pdg::IsNeutronOrProton(
+ fGenieGHep->Summary()->InitState().TgtPtr()->HitNucPdg())) { // Check that it isn't coherent! Want to keep initial state if so
+ if (state == kInitialState)
+ state = kNuclearInitial;
+ else if (state == kFinalState)
+ state = kNuclearRemnant;
+ }
}
return state;
}
#endif
#ifdef __GENIE_ENABLED__
int GENIEInputHandler::ConvertGENIEReactionCode(GHepRecord *gheprec) {
// I randomly picked 53 here because NEUT doesn't have an appropriate mode...
if (gheprec->Summary()->ProcInfo().IsNuElectronElastic()){
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg())) return 53;
else return -53;
}
// And the same story for 54
if (gheprec->Summary()->ProcInfo().IsIMDAnnihilation()){
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg())) return 54;
else return -54;
}
// Electron Scattering
if (gheprec->Summary()->ProcInfo().IsEM()) {
if (gheprec->Summary()->InitState().ProbePdg() == 11) {
if (gheprec->Summary()->ProcInfo().IsQuasiElastic())
return 1;
else if (gheprec->Summary()->ProcInfo().IsMEC())
return 2;
else if (gheprec->Summary()->ProcInfo().IsResonant())
return 13;
else if (gheprec->Summary()->ProcInfo().IsDeepInelastic())
return 26;
else {
NUIS_ERR(WRN,
- "Unknown GENIE Electron Scattering Mode!"
- << std::endl
- << "ScatteringTypeId = "
- << gheprec->Summary()->ProcInfo().ScatteringTypeId() << " "
- << "InteractionTypeId = "
- << gheprec->Summary()->ProcInfo().InteractionTypeId()
- << std::endl
- << genie::ScatteringType::AsString(
- gheprec->Summary()->ProcInfo().ScatteringTypeId())
- << " "
- << genie::InteractionType::AsString(
- gheprec->Summary()->ProcInfo().InteractionTypeId())
- << " " << gheprec->Summary()->ProcInfo().IsMEC());
+ "Unknown GENIE Electron Scattering Mode!"
+ << std::endl
+ << "ScatteringTypeId = "
+ << gheprec->Summary()->ProcInfo().ScatteringTypeId() << " "
+ << "InteractionTypeId = "
+ << gheprec->Summary()->ProcInfo().InteractionTypeId()
+ << std::endl
+ << genie::ScatteringType::AsString(
+ gheprec->Summary()->ProcInfo().ScatteringTypeId())
+ << " "
+ << genie::InteractionType::AsString(
+ gheprec->Summary()->ProcInfo().InteractionTypeId())
+ << " " << gheprec->Summary()->ProcInfo().IsMEC());
return 0;
}
}
// Weak CC
} else if (gheprec->Summary()->ProcInfo().IsWeakCC()) {
// CC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 2;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -2;
#ifndef GENIE_PRE_R3
} else if (gheprec->Summary()->ProcInfo().IsDiffractive()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 15;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -15;
#endif
// CC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
// Weak NC
} else if (gheprec->Summary()->ProcInfo().IsWeakNC()) {
// NC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 32;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -32;
#ifndef GENIE_PRE_R3
} else if (gheprec->Summary()->ProcInfo().IsDiffractive()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 35;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -35;
#endif
// NC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
}
return 0;
}
void GENIEInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Check for GENIE Event
if (!fGenieNtpl)
return;
if (!fGenieNtpl->event)
return;
// Cast Event Record
fGenieGHep = static_cast<GHepRecord *>(fGenieNtpl->event);
if (!fGenieGHep)
return;
// Convert GENIE Reaction Code
fNUISANCEEvent->Mode = ConvertGENIEReactionCode(fGenieGHep);
if (!fNUISANCEEvent->Mode) {
NUIS_ERR(WRN, "Failed to determine mode for GENIE event: ");
std::cout << *fGenieGHep << std::endl;
}
// Set Event Info
fNUISANCEEvent->fEventNo = 0.0;
fNUISANCEEvent->fTotCrs = fGenieGHep->XSec();
+
// Have a bool storing if interaction happened on free or bound nucleon
bool IsFree = false;
+
// Set the TargetPDG
if (fGenieGHep->TargetNucleus() != NULL) {
fNUISANCEEvent->fTargetPDG = fGenieGHep->TargetNucleus()->Pdg();
IsFree = false;
// Sometimes GENIE scatters off free nucleons, electrons, photons
// In which TargetNucleus is NULL and we need to find the initial state
// particle
} else {
// Check the particle is an initial state particle
// Follows GHepRecord::TargetNucleusPosition but doesn't do check on
// pdg::IsIon
GHepParticle *p = fGenieGHep->Particle(1);
// Check that particle 1 actually exists
if (!p) {
NUIS_ABORT("Can't find particle 1 for GHepRecord");
}
// If not an ion but is an initial state particle
if (!pdg::IsIon(p->Pdg()) && p->Status() == kIStInitialState) {
IsFree = true;
fNUISANCEEvent->fTargetPDG = p->Pdg();
// Catch if something strange happens:
// Here particle 1 is not an initial state particle OR
// particle 1 is an ion OR
// both
} else {
if (pdg::IsIon(p->Pdg())) {
NUIS_ABORT(
"Particle 1 in GHepRecord stack is an ion but isn't an initial "
"state particle");
} else {
NUIS_ABORT(
"Particle 1 in GHepRecord stack is not an ion but is an initial "
"state particle");
}
}
}
// Set the A and Z and H from the target PDG
// Depends on if we scattered off a free or bound nucleon
if (!IsFree) {
fNUISANCEEvent->fTargetA =
- TargetUtils::GetTargetAFromPDG(fNUISANCEEvent->fTargetPDG);
+ TargetUtils::GetTargetAFromPDG(fNUISANCEEvent->fTargetPDG);
fNUISANCEEvent->fTargetZ =
- TargetUtils::GetTargetZFromPDG(fNUISANCEEvent->fTargetPDG);
+ TargetUtils::GetTargetZFromPDG(fNUISANCEEvent->fTargetPDG);
fNUISANCEEvent->fTargetH = 0;
} else {
// If free proton scattering
if (fNUISANCEEvent->fTargetPDG == 2212) {
fNUISANCEEvent->fTargetA = 1;
fNUISANCEEvent->fTargetZ = 1;
fNUISANCEEvent->fTargetH = 1;
// If free neutron scattering
} else if (fNUISANCEEvent->fTargetPDG == 2112) {
fNUISANCEEvent->fTargetA = 0;
fNUISANCEEvent->fTargetZ = 1;
fNUISANCEEvent->fTargetH = 0;
// If neither
} else {
fNUISANCEEvent->fTargetA = 0;
fNUISANCEEvent->fTargetZ = 0;
fNUISANCEEvent->fTargetH = 0;
}
}
fNUISANCEEvent->fBound = !IsFree;
fNUISANCEEvent->InputWeight =
- 1.0; //(1E+38 / genie::units::cm2) * fGenieGHep->XSec();
-
- // And the custom weights
- if (fNOvAWeights) {
- fNUISANCEEvent->CustomWeight = NOVAw;
- fNUISANCEEvent->CustomWeightArray[0] = MAQEw;
- fNUISANCEEvent->CustomWeightArray[1] = NonResw;
- fNUISANCEEvent->CustomWeightArray[2] = RPAQEw;
- fNUISANCEEvent->CustomWeightArray[3] = RPARESw;
- fNUISANCEEvent->CustomWeightArray[4] = MECw;
- fNUISANCEEvent->CustomWeightArray[5] = NOVAw;
- } else {
- fNUISANCEEvent->CustomWeight = 1.0;
- fNUISANCEEvent->CustomWeightArray[0] = 1.0;
- fNUISANCEEvent->CustomWeightArray[1] = 1.0;
- fNUISANCEEvent->CustomWeightArray[2] = 1.0;
- fNUISANCEEvent->CustomWeightArray[3] = 1.0;
- fNUISANCEEvent->CustomWeightArray[4] = 1.0;
- fNUISANCEEvent->CustomWeightArray[5] = 1.0;
- }
+ 1.0; //(1E+38 / genie::units::cm2) * fGenieGHep->XSec();
// Get N Particle Stack
unsigned int npart = fGenieGHep->GetEntries();
unsigned int kmax = fNUISANCEEvent->kMaxParticles;
if (npart > kmax) {
NUIS_ERR(WRN, "GENIE has too many particles, expanding stack.");
fNUISANCEEvent->ExpandParticleStack(npart);
}
// Fill Particle Stack
GHepParticle *p = 0;
TObjArrayIter iter(fGenieGHep);
fNUISANCEEvent->fNParticles = 0;
// Loop over all particles
while ((p = (dynamic_cast<genie::GHepParticle *>((iter).Next())))) {
if (!p)
continue;
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
// Remove Undefined
if (kRemoveUndefParticles && state == kUndefinedState)
continue;
// Remove FSI
if (kRemoveFSIParticles && state == kFSIState)
continue;
if (kRemoveNuclearParticles &&
(state == kNuclearInitial || state == kNuclearRemnant))
continue;
// Fill Vectors
int curpart = fNUISANCEEvent->fNParticles;
fNUISANCEEvent->fParticleState[curpart] = state;
// Mom
fNUISANCEEvent->fParticleMom[curpart][0] = p->Px() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][1] = p->Py() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][2] = p->Pz() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][3] = p->E() * 1.E3;
// PDG
fNUISANCEEvent->fParticlePDG[curpart] = p->Pdg();
// Set if the particle was on the fundamental vertex
- fNUISANCEEvent->fPrimaryVertex[curpart] = (p->FirstMother() < 2);
+ fNUISANCEEvent->fPrimaryVertex[curpart] = IsPrimary(p);
// Add to N particle count
fNUISANCEEvent->fNParticles++;
// Extra Check incase GENIE fails.
if ((UInt_t)fNUISANCEEvent->fNParticles == kmax) {
NUIS_ERR(WRN, "Number of GENIE Particles exceeds maximum!");
NUIS_ERR(WRN, "Extend kMax, or run without including FSI particles!");
break;
}
}
// Fill Extra Stack
if (fSaveExtra)
fGenieInfo->FillGeneratorInfo(fGenieNtpl);
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
FitParticle *ISAnyLepton = fNUISANCEEvent->GetHMISAnyLeptons();
if (ISAnyLepton) {
fNUISANCEEvent->probe_E = ISAnyLepton->E();
fNUISANCEEvent->probe_pdg = ISAnyLepton->PDG();
}
return;
}
void GENIEInputHandler::Print() {}
#endif
diff --git a/src/InputHandler/GENIEInputHandler.h b/src/InputHandler/GENIEInputHandler.h
index a09cda0..44d5d1a 100644
--- a/src/InputHandler/GENIEInputHandler.h
+++ b/src/InputHandler/GENIEInputHandler.h
@@ -1,130 +1,121 @@
// Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef GENIEINPUTHANDLER_H
#define GENIEINPUTHANDLER_H
/*!
* \addtogroup InputHandler
* @{
*/
#ifdef __GENIE_ENABLED__
#include "InputHandler.h"
#include "InputUtils.h"
#include "PlotUtils.h"
#ifdef GENIE_PRE_R3
#include "GHEP/GHepParticle.h"
#include "PDG/PDGUtils.h"
#include "GHEP/GHepUtils.h"
#include "Conventions/Units.h"
#include "EVGCore/EventRecord.h"
#include "GHEP/GHepRecord.h"
#include "Ntuple/NtpMCEventRecord.h"
#else
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Framework/GHEP/GHepUtils.h"
#include "Framework/Conventions/Units.h"
#include "Framework/EventGen/EventRecord.h"
#include "Framework/GHEP/GHepRecord.h"
#include "Framework/Ntuple/NtpMCEventRecord.h"
#endif
using namespace genie;
/// GENIE Generator Container to save extra particle status codes.
class GENIEGeneratorInfo : public GeneratorInfoBase {
-public:
- GENIEGeneratorInfo() {};
- virtual ~GENIEGeneratorInfo();
+ public:
+ GENIEGeneratorInfo() {};
+ virtual ~GENIEGeneratorInfo();
- /// Assigns information to branches
- void AddBranchesToTree(TTree* tn);
+ /// Assigns information to branches
+ void AddBranchesToTree(TTree* tn);
- /// Setup reading information from branches
- void SetBranchesFromTree(TTree* tn);
+ /// Setup reading information from branches
+ void SetBranchesFromTree(TTree* tn);
- /// Allocate any dynamic arrays for a new particle stack size
- void AllocateParticleStack(int stacksize);
+ /// Allocate any dynamic arrays for a new particle stack size
+ void AllocateParticleStack(int stacksize);
- /// Clear any dynamic arrays
- void DeallocateParticleStack();
+ /// Clear any dynamic arrays
+ void DeallocateParticleStack();
- /// Read extra genie information from the event
- void FillGeneratorInfo(NtpMCEventRecord* ntpl);
+ /// Read extra genie information from the event
+ void FillGeneratorInfo(NtpMCEventRecord* ntpl);
- /// Reset extra information to default/empty values
- void Reset();
+ /// Reset extra information to default/empty values
+ void Reset();
- int kMaxParticles; ///< Number of particles in stack
- int* fGenieParticlePDGs; ///< GENIE Particle PDGs (example)
+ int kMaxParticles; ///< Number of particles in stack
+ int* fGenieParticlePDGs; ///< GENIE Particle PDGs (example)
};
/// Main GENIE InputHandler
class GENIEInputHandler : public InputHandlerBase {
-public:
+ public:
- /// Standard constructor given a name and input files
- GENIEInputHandler(std::string const& handle, std::string const& rawinputs);
- virtual ~GENIEInputHandler();
+ /// Standard constructor given a name and input files
+ GENIEInputHandler(std::string const& handle, std::string const& rawinputs);
+ virtual ~GENIEInputHandler();
- /// Create a TTree Cache to speed up file read
- void CreateCache();
+ /// Create a TTree Cache to speed up file read
+ void CreateCache();
- /// Remove TTree Cache to save memory
- void RemoveCache();
-
- /// Returns a NUISANCE format event from the GENIE TTree. If !lightweight
- /// then CalcNUISANCEKinematics() is called to convert the GENIE event into
- /// a standard NUISANCE format.
- FitEvent* GetNuisanceEvent(const UInt_t entry, const bool lightweight = false);
+ /// Remove TTree Cache to save memory
+ void RemoveCache();
- /// Converts GENIE event into standard NUISANCE FitEvent by looping over all
- /// particles in the event and adding them to stack in fNUISANCEEvent.
- void CalcNUISANCEKinematics();
+ /// Returns a NUISANCE format event from the GENIE TTree. If !lightweight
+ /// then CalcNUISANCEKinematics() is called to convert the GENIE event into
+ /// a standard NUISANCE format.
+ FitEvent* GetNuisanceEvent(const UInt_t entry, const bool lightweight = false);
- /// Placeholder for GENIE related event printing.
- void Print();
+ /// Converts GENIE event into standard NUISANCE FitEvent by looping over all
+ /// particles in the event and adding them to stack in fNUISANCEEvent.
+ void CalcNUISANCEKinematics();
- /// Converts GENIE particle status codes into NUISANCE status codes.
- int GetGENIEParticleStatus(genie::GHepParticle* part, int mode = 0);
+ /// Placeholder for GENIE related event printing.
+ void Print();
- /// Converts GENIE event reaction codes into NUISANCE reaction codes.
- int ConvertGENIEReactionCode(GHepRecord* gheprec);
+ /// Converts GENIE particle status codes into NUISANCE status codes.
+ int GetGENIEParticleStatus(genie::GHepParticle* part, int mode = 0);
- GHepRecord* fGenieGHep; ///< Pointer to actual event record
- NtpMCEventRecord* fGenieNtpl; ///< Ntpl Wrapper Class
+ /// Converts GENIE event reaction codes into NUISANCE reaction codes.
+ int ConvertGENIEReactionCode(GHepRecord* gheprec);
- TChain* fGENIETree; ///< Main GENIE Event TTree
- bool fSaveExtra; ///< Flag to save Extra GENIE info into Nuisance Event
- GENIEGeneratorInfo* fGenieInfo; ///< Extra GENIE Generator Info Writer
+ GHepRecord* fGenieGHep; ///< Pointer to actual event record
+ NtpMCEventRecord* fGenieNtpl; ///< Ntpl Wrapper Class
- bool fNOvAWeights; ///< Flag to save nova weights or not
+ TChain* fGENIETree; ///< Main GENIE Event TTree
+ bool fSaveExtra; ///< Flag to save Extra GENIE info into Nuisance Event
+ GENIEGeneratorInfo* fGenieInfo; ///< Extra GENIE Generator Info Writer
- // Extra weights from Jeremy for NOvA weights
- double MAQEw;
- double NonResw;
- double RPAQEw;
- double RPARESw;
- double MECw;
- double DISw;
- double NOVAw;
+ bool IsPrimary(GHepParticle *p);
};
/*! @} */
#endif
#endif