diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx
index 5a38a97..24dad08 100644
--- a/src/InputHandler/GENIEInputHandler.cxx
+++ b/src/InputHandler/GENIEInputHandler.cxx
@@ -1,594 +1,599 @@
// Copyright 2016 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 .
*******************************************************************************/
#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(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((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;
}
}
GENIEInputHandler::GENIEInputHandler(std::string const &handle,
std::string const &rawinputs) {
NUIS_LOG(SAM, "Creating GENIEInputHandler : " << handle);
genie::Messenger::Instance()->SetPriorityLevel("GHepUtils", pFATAL);
// 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 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]);
+ 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!");
+ 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]);
+ << 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);
+ << inputs[inp_it] << " but you specified it" << std::endl);
} else {
NUIS_LOG(FIT, "Found nova_wgts tree in file " << 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 entry,
+FitEvent *GENIEInputHandler::GetNuisanceEvent(const UInt_t ent,
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(fGenieNtpl->event);
if (!fGenieGHep)
return NULL;
TObjArrayIter iter(fGenieGHep);
genie::GHepParticle *p;
while ((p = (dynamic_cast((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) {
/*
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)
state = kInitialState;
else
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;
}
// Flag to remove nuclear part in genie
if (p->Pdg() > 1000000) {
if (state == kInitialState)
state = kNuclearInitial;
else if (state == kFinalState)
state = kNuclearRemnant;
}
return state;
}
#endif
#ifdef __GENIE_ENABLED__
int GENIEInputHandler::ConvertGENIEReactionCode(GHepRecord *gheprec) {
// 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;
// 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;
// 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(fGenieNtpl->event);
if (!fGenieGHep)
return;
// Convert GENIE Reaction Code
fNUISANCEEvent->Mode = ConvertGENIEReactionCode(fGenieGHep);
// 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");
+ 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");
+ 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);
fNUISANCEEvent->fTargetZ =
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;
}
// 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((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);
// 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/GIBUUInputHandler.cxx b/src/InputHandler/GIBUUInputHandler.cxx
index ce90bb8..f83c563 100644
--- a/src/InputHandler/GIBUUInputHandler.cxx
+++ b/src/InputHandler/GIBUUInputHandler.cxx
@@ -1,300 +1,303 @@
#ifdef __GiBUU_ENABLED__
#include "GIBUUInputHandler.h"
#include "InputUtils.h"
GIBUUGeneratorInfo::~GIBUUGeneratorInfo() { DeallocateParticleStack(); }
void GIBUUGeneratorInfo::AddBranchesToTree(TTree *tn) {
// tn->Branch("NEUTParticleN", fNEUTParticleN, "NEUTParticleN/I");
// tn->Branch("NEUTParticleStatusCode", fNEUTParticleStatusCode,
// "NEUTParticleStatusCode[NEUTParticleN]/I");
// tn->Branch("NEUTParticleAliveCode", fNEUTParticleAliveCode,
// "NEUTParticleAliveCode[NEUTParticleN]/I");
}
void GIBUUGeneratorInfo::SetBranchesFromTree(TTree *tn) {
// tn->SetBranchAddress("NEUTParticleN", &fNEUTParticleN );
// tn->SetBranchAddress("NEUTParticleStatusCode", &fNEUTParticleStatusCode );
// tn->SetBranchAddress("NEUTParticleAliveCode", &fNEUTParticleAliveCode );
}
void GIBUUGeneratorInfo::AllocateParticleStack(int stacksize) {
// fNEUTParticleN = 0;
// fNEUTParticleStatusCode = new int[stacksize];
// fNEUTParticleStatusCode = new int[stacksize];
}
void GIBUUGeneratorInfo::DeallocateParticleStack() {
// delete fNEUTParticleStatusCode;
// delete fNEUTParticleAliveCode;
}
void GIBUUGeneratorInfo::FillGeneratorInfo(GiBUUStdHepReader *nevent) {
Reset();
// for (int i = 0; i < nevent->Npart(); i++) {
// fNEUTParticleStatusCode[i] = nevent->PartInfo(i)->fStatus;
// fNEUTParticleAliveCode[i] = nevent->PartInfo(i)->fIsAlive;
// fNEUTParticleN++;
// }
}
void GIBUUGeneratorInfo::Reset() {
// for (int i = 0; i < fNEUTParticleN; i++) {
// fNEUTParticleStatusCode[i] = -1;
// fNEUTParticleAliveCode[i] = 9;
// }
// fNEUTParticleN = 0;
}
GIBUUInputHandler::GIBUUInputHandler(std::string const &handle,
std::string const &rawinputs) {
NUIS_LOG(SAM, "Creating GiBUUInputHandler : " << handle);
// Run a joint input handling
fName = handle;
fEventType = kGiBUU;
fGIBUUTree = new TChain("giRooTracker");
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
NUIS_LOG(SAM, "Opening event file " << inputs[inp_it]);
TFile *inp_file = new TFile(inputs[inp_it].c_str(), "READ");
if ((!inp_file) || (!inp_file->IsOpen())) {
- NUIS_ABORT("GiBUU file !IsOpen() at : '"
- << inputs[inp_it] << "'" << std::endl
- << "Check that your file paths are correct and the file exists!");
+ NUIS_ABORT(
+ "GiBUU file !IsOpen() at : '"
+ << inputs[inp_it] << "'" << std::endl
+ << "Check that your file paths are correct and the file exists!");
}
int NFluxes = bool(dynamic_cast(inp_file->Get("numu_flux"))) +
bool(dynamic_cast(inp_file->Get("numub_flux"))) +
bool(dynamic_cast(inp_file->Get("nue_flux"))) +
bool(dynamic_cast(inp_file->Get("nueb_flux"))) +
bool(dynamic_cast(inp_file->Get("e_flux")));
if (NFluxes != 1) {
NUIS_ABORT("Found " << NFluxes << " input fluxes in " << inputs[inp_it]
- << ". The NUISANCE GiBUU interface expects to be "
- "passed multiple species vectors as separate "
- "input files like: "
- "\"GiBUU:(MINERVA_FHC_numu_evts.root,MINERVA_FHC_"
- "numubar_evts.root,[...])\"");
+ << ". The NUISANCE GiBUU interface expects to be "
+ "passed multiple species vectors as separate "
+ "input files like: "
+ "\"GiBUU:(MINERVA_FHC_numu_evts.root,MINERVA_FHC_"
+ "numubar_evts.root,[...])\"");
}
// Get Flux/Event hist
TH1D *fluxhist = dynamic_cast(inp_file->Get("flux"));
TH1D *eventhist = dynamic_cast(inp_file->Get("evt"));
if (!fluxhist || !eventhist) {
NUIS_ERR(FTL, "Input File Contents: " << inputs[inp_it]);
inp_file->ls();
NUIS_ABORT("GiBUU FILE doesn't contain flux/xsec info. You may have to "
- "regenerate your MC!");
+ "regenerate your MC!");
}
// Get N Events
TTree *giRooTracker = dynamic_cast(inp_file->Get("giRooTracker"));
if (!giRooTracker) {
- NUIS_ERR(FTL,
- "giRooTracker Tree not located in NEUT file: " << inputs[inp_it]);
- NUIS_ABORT("Check your inputs, they may need to be completely regenerated!");
+ NUIS_ERR(FTL, "giRooTracker Tree not located in NEUT file: "
+ << inputs[inp_it]);
+ NUIS_ABORT(
+ "Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevents = giRooTracker->GetEntries();
if (nevents <= 0) {
NUIS_ABORT("Trying to a TTree with "
- << nevents << " to TChain from : " << 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
fGIBUUTree->AddFile(inputs[inp_it].c_str());
}
// Registor all our file inputs
SetupJointInputs();
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fGiReader = new GiBUUStdHepReader();
fGiReader->SetBranchAddresses(fGIBUUTree);
fNUISANCEEvent->HardReset();
};
-FitEvent *GIBUUInputHandler::GetNuisanceEvent(const UInt_t entry,
+FitEvent *GIBUUInputHandler::GetNuisanceEvent(const UInt_t ent,
const bool lightweight) {
+ UInt_t entry = ent + fSkip;
// Check out of bounds
if (entry >= (UInt_t)fNEvents)
return NULL;
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fGIBUUTree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
else {
for (int i = 0; i < fGiReader->StdHepN; i++) {
int state = GetGIBUUParticleStatus(fGiReader->StdHepStatus[i],
fGiReader->StdHepPdg[i]);
if (state != kInitialState) {
continue;
}
if (std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4,
fGiReader->StdHepPdg[i])) {
fNUISANCEEvent->probe_E = fGiReader->StdHepP4[i][3] * 1.E3;
fNUISANCEEvent->probe_pdg = fGiReader->StdHepPdg[i];
break;
}
}
}
#endif
fNUISANCEEvent->InputWeight *= GetInputWeight(entry);
fNUISANCEEvent->GiRead = fGiReader;
return fNUISANCEEvent;
}
int GetGIBUUParticleStatus(int status, int pdg) {
int state = kUndefinedState;
switch (status) {
case 0: // Incoming
case 11: // Struck nucleon
state = kInitialState;
break;
case 1: // Good Final State
state = kFinalState;
break;
default: // Other
break;
}
// Set Nuclear States Flag
if (pdg > 1000000) {
if (state == kInitialState)
state = kNuclearInitial;
else if (state == kFinalState)
state = kNuclearRemnant;
else
state = kUndefinedState;
}
return state;
}
void GIBUUInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
FitEvent *evt = fNUISANCEEvent;
evt->Mode = fGiReader->GiBUU2NeutCode;
evt->fEventNo = 0.0;
evt->fTotCrs = 0;
evt->fTargetA = 0.0; // Change to get these from nuclear remnant.
evt->fTargetZ = 0.0;
evt->fTargetH = 0;
evt->fBound = 0.0;
// Extra GiBUU Input Weight
evt->InputWeight = fGiReader->EvtWght;
// Check Stack N
int npart = fGiReader->StdHepN;
int kmax = evt->kMaxParticles;
if ((UInt_t)npart > (UInt_t)kmax) {
NUIS_ERR(WRN, "GiBUU has too many particles. Expanding Stack.");
fNUISANCEEvent->ExpandParticleStack(npart);
}
// Create Stack
evt->fNParticles = 0;
for (int i = 0; i < npart; i++) {
// State
int state = GetGIBUUParticleStatus(fGiReader->StdHepStatus[i],
fGiReader->StdHepPdg[i]);
int curpart = evt->fNParticles;
// Set State
evt->fParticleState[evt->fNParticles] = state;
// Mom
evt->fParticleMom[curpart][0] = fGiReader->StdHepP4[i][0] * 1.E3;
evt->fParticleMom[curpart][1] = fGiReader->StdHepP4[i][1] * 1.E3;
evt->fParticleMom[curpart][2] = fGiReader->StdHepP4[i][2] * 1.E3;
evt->fParticleMom[curpart][3] = fGiReader->StdHepP4[i][3] * 1.E3;
// PDG
evt->fParticlePDG[curpart] = fGiReader->StdHepPdg[i];
// Add to total particles
evt->fNParticles++;
}
// 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 GIBUUInputHandler::Print() {}
void GIBUUInputHandler::SetupJointInputs() {
if (jointeventinputs.size() <= 1) {
jointinput = false;
} else if (jointeventinputs.size() > 1) {
jointinput = true;
jointindexswitch = 0;
}
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
if (fMaxEvents != -1 and jointeventinputs.size() > 1) {
NUIS_ABORT("Can only handle joint inputs when config MAXEVENTS = -1!");
}
for (size_t i = 0; i < jointeventinputs.size(); i++) {
double scale = double(fNEvents) / fEventHist->Integral("width");
scale *= jointfluxinputs.at(i)->Integral("width");
jointindexscale.push_back(scale);
}
fEventHist->SetNameTitle((fName + "_EVT").c_str(), (fName + "_EVT").c_str());
fFluxHist->SetNameTitle((fName + "_FLUX").c_str(), (fName + "_FLUX").c_str());
// Setup Max Events
if (fMaxEvents > 1 && fMaxEvents < fNEvents) {
if (LOG_LEVEL(SAM)) {
std::cout << "\t\t|-> Read Max Entries : " << fMaxEvents << std::endl;
}
fNEvents = fMaxEvents;
}
// Print out Status
if (LOG_LEVEL(SAM)) {
std::cout << "\t\t|-> Total Entries : " << fNEvents << std::endl
<< "\t\t|-> Event Integral : "
<< fEventHist->Integral("width") * 1.E-38 << " events/nucleon"
<< std::endl
<< "\t\t|-> Flux Integral : " << fFluxHist->Integral("width")
<< " /cm2" << std::endl
<< "\t\t|-> Event/Flux : "
<< fEventHist->Integral("width") * 1.E-38 /
fFluxHist->Integral("width")
<< " cm2/nucleon" << std::endl;
}
}
#endif
diff --git a/src/InputHandler/InputHandler.cxx b/src/InputHandler/InputHandler.cxx
index 7ef8800..304b445 100644
--- a/src/InputHandler/InputHandler.cxx
+++ b/src/InputHandler/InputHandler.cxx
@@ -1,299 +1,311 @@
// Copyright 2016 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 .
-*******************************************************************************/
+ * 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 .
+ *******************************************************************************/
#include "InputHandler.h"
#include "InputUtils.h"
InputHandlerBase::InputHandlerBase() {
fName = "";
fFluxHist = NULL;
fEventHist = NULL;
fNEvents = 0;
fNUISANCEEvent = NULL;
fBaseEvent = NULL;
kRemoveUndefParticles = FitPar::Config().GetParB("RemoveUndefParticles");
kRemoveFSIParticles = FitPar::Config().GetParB("RemoveFSIParticles");
kRemoveNuclearParticles = FitPar::Config().GetParB("RemoveNuclearParticles");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
fTTreePerformance = NULL;
+ fSkip = 0;
+ if (FitPar::Config().HasConfig("NSKIPEVENTS")) {
+ fSkip = FitPar::Config().GetParI("NSKIPEVENTS");
+ }
};
InputHandlerBase::~InputHandlerBase() {
- if (fFluxHist) delete fFluxHist;
- if (fEventHist) delete fEventHist;
+ if (fFluxHist)
+ delete fFluxHist;
+ if (fEventHist)
+ delete fEventHist;
// if (fXSecHist) delete fXSecHist;
// if (fNUISANCEEvent) delete fNUISANCEEvent;
jointfluxinputs.clear();
jointeventinputs.clear();
jointindexlow.clear();
jointindexhigh.clear();
jointindexallowed.clear();
jointindexscale.clear();
// if (fTTreePerformance) {
// fTTreePerformance->SaveAs(("ttreeperfstats_" + fName +
// ".root").c_str());
// }
}
void InputHandlerBase::Print(){};
-TH1D* InputHandlerBase::GetXSecHistogram(void) {
- fXSecHist = (TH1D*)fFluxHist->Clone();
+TH1D *InputHandlerBase::GetXSecHistogram(void) {
+ fXSecHist = (TH1D *)fFluxHist->Clone();
fXSecHist->Divide(fEventHist);
return fXSecHist;
};
double InputHandlerBase::PredictedEventRate(double low, double high,
std::string intOpt) {
Int_t minBin = fEventHist->GetXaxis()->FindFixBin(low);
Int_t maxBin = fEventHist->GetXaxis()->FindFixBin(high);
if ((fEventHist->IsBinOverflow(minBin) && (low != -9999.9))) {
minBin = 1;
}
if ((fEventHist->IsBinOverflow(maxBin) && (high != -9999.9))) {
maxBin = fEventHist->GetXaxis()->GetNbins() + 1;
}
// If we are within a single bin
if (minBin == maxBin) {
// Get the contained fraction of the single bin's width
return ((high - low) / fEventHist->GetXaxis()->GetBinWidth(minBin)) *
fEventHist->Integral(minBin, minBin, intOpt.c_str());
}
double lowBinUpEdge = fEventHist->GetXaxis()->GetBinUpEdge(minBin);
double highBinLowEdge = fEventHist->GetXaxis()->GetBinLowEdge(maxBin);
double lowBinfracIntegral =
((lowBinUpEdge - low) / fEventHist->GetXaxis()->GetBinWidth(minBin)) *
fEventHist->Integral(minBin, minBin, intOpt.c_str());
double highBinfracIntegral =
((high - highBinLowEdge) / fEventHist->GetXaxis()->GetBinWidth(maxBin)) *
fEventHist->Integral(maxBin, maxBin, intOpt.c_str());
// If they are neighbouring bins
if ((minBin + 1) == maxBin) {
std::cout << "Get lowfrac + highfrac" << std::endl;
// Get the contained fraction of the two bin's width
return lowBinfracIntegral + highBinfracIntegral;
}
double ContainedIntegral =
fEventHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str());
// If there are filled bins between them
return lowBinfracIntegral + highBinfracIntegral + ContainedIntegral;
};
double InputHandlerBase::TotalIntegratedFlux(double low, double high,
std::string intOpt) {
Int_t minBin = fFluxHist->GetXaxis()->FindFixBin(low);
Int_t maxBin = fFluxHist->GetXaxis()->FindFixBin(high);
if ((fFluxHist->IsBinOverflow(minBin) && (low != -9999.9))) {
minBin = 1;
}
if ((fFluxHist->IsBinOverflow(maxBin) && (high != -9999.9))) {
maxBin = fFluxHist->GetXaxis()->GetNbins();
- high = fFluxHist->GetXaxis()->GetBinLowEdge(maxBin+1);
+ high = fFluxHist->GetXaxis()->GetBinLowEdge(maxBin + 1);
}
// If we are within a single bin
if (minBin == maxBin) {
// Get the contained fraction of the single bin's width
return ((high - low) / fFluxHist->GetXaxis()->GetBinWidth(minBin)) *
fFluxHist->Integral(minBin, minBin, intOpt.c_str());
}
double lowBinUpEdge = fFluxHist->GetXaxis()->GetBinUpEdge(minBin);
double highBinLowEdge = fFluxHist->GetXaxis()->GetBinLowEdge(maxBin);
double lowBinfracIntegral =
((lowBinUpEdge - low) / fFluxHist->GetXaxis()->GetBinWidth(minBin)) *
fFluxHist->Integral(minBin, minBin, intOpt.c_str());
double highBinfracIntegral =
((high - highBinLowEdge) / fFluxHist->GetXaxis()->GetBinWidth(maxBin)) *
fFluxHist->Integral(maxBin, maxBin, intOpt.c_str());
// If they are neighbouring bins
if ((minBin + 1) == maxBin) {
std::cout << "Get lowfrac + highfrac" << std::endl;
// Get the contained fraction of the two bin's width
return lowBinfracIntegral + highBinfracIntegral;
}
double ContainedIntegral =
fFluxHist->Integral(minBin + 1, maxBin - 1, intOpt.c_str());
// If there are filled bins between them
return lowBinfracIntegral + highBinfracIntegral + ContainedIntegral;
}
-std::vector InputHandlerBase::GetFluxList(void) {
- return std::vector(1, fFluxHist);
+std::vector InputHandlerBase::GetFluxList(void) {
+ return std::vector(1, fFluxHist);
};
-std::vector InputHandlerBase::GetEventList(void) {
- return std::vector(1, fEventHist);
+std::vector InputHandlerBase::GetEventList(void) {
+ return std::vector(1, fEventHist);
};
-std::vector InputHandlerBase::GetXSecList(void) {
- return std::vector(1, GetXSecHistogram());
+std::vector InputHandlerBase::GetXSecList(void) {
+ return std::vector(1, GetXSecHistogram());
};
-FitEvent* InputHandlerBase::FirstNuisanceEvent() {
+FitEvent *InputHandlerBase::FirstNuisanceEvent() {
fCurrentIndex = 0;
return GetNuisanceEvent(fCurrentIndex);
};
-FitEvent* InputHandlerBase::NextNuisanceEvent() {
+FitEvent *InputHandlerBase::NextNuisanceEvent() {
fCurrentIndex++;
if ((fMaxEvents != -1) && (fCurrentIndex > fMaxEvents)) {
return NULL;
}
return GetNuisanceEvent(fCurrentIndex);
};
-BaseFitEvt* InputHandlerBase::FirstBaseEvent() {
+BaseFitEvt *InputHandlerBase::FirstBaseEvent() {
fCurrentIndex = 0;
return GetBaseEvent(fCurrentIndex);
};
-BaseFitEvt* InputHandlerBase::NextBaseEvent() {
+BaseFitEvt *InputHandlerBase::NextBaseEvent() {
fCurrentIndex++;
if (jointinput and fMaxEvents != -1) {
while (fCurrentIndex < jointindexlow[jointindexswitch] ||
fCurrentIndex >= jointindexhigh[jointindexswitch]) {
jointindexswitch++;
// Loop Around
if (jointindexswitch == jointindexlow.size()) {
jointindexswitch = 0;
}
}
if (fCurrentIndex >
jointindexlow[jointindexswitch] + jointindexallowed[jointindexswitch]) {
fCurrentIndex = jointindexlow[jointindexswitch];
}
}
return GetBaseEvent(fCurrentIndex);
};
-void InputHandlerBase::RegisterJointInput(std::string input, int n, TH1D* f,
- TH1D* e) {
+void InputHandlerBase::RegisterJointInput(std::string input, int n, TH1D *f,
+ TH1D *e) {
if (jointfluxinputs.size() == 0) {
jointindexswitch = 0;
fNEvents = 0;
}
// Push into individual input vectors
- jointfluxinputs.push_back((TH1D*)f->Clone());
- jointeventinputs.push_back((TH1D*)e->Clone());
+ jointfluxinputs.push_back((TH1D *)f->Clone());
+ jointeventinputs.push_back((TH1D *)e->Clone());
jointindexlow.push_back(fNEvents);
jointindexhigh.push_back(fNEvents + n);
fNEvents += n;
// Add to the total flux/event hist
- if (!fFluxHist) fFluxHist = (TH1D*)f->Clone();
- else fFluxHist->Add(f);
-
- if (!fEventHist) fEventHist = (TH1D*)e->Clone();
- else fEventHist->Add(e);
+ if (!fFluxHist)
+ fFluxHist = (TH1D *)f->Clone();
+ else
+ fFluxHist->Add(f);
+
+ if (!fEventHist)
+ fEventHist = (TH1D *)e->Clone();
+ else
+ fEventHist->Add(e);
}
void InputHandlerBase::SetupJointInputs() {
if (jointeventinputs.size() <= 1) {
jointinput = false;
} else if (jointeventinputs.size() > 1) {
jointinput = true;
jointindexswitch = 0;
}
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
if (fMaxEvents != -1 and jointeventinputs.size() > 1) {
NUIS_ABORT("Can only handle joint inputs when config MAXEVENTS = -1!");
}
if (jointeventinputs.size() > 1) {
- NUIS_ERR(WRN,
- "GiBUU sample contains multiple inputs. This will only work for "
- "samples that expect multi-species inputs. If this sample does, you "
- "can ignore this warning.");
+ NUIS_ERR(
+ WRN,
+ "GiBUU sample contains multiple inputs. This will only work for "
+ "samples that expect multi-species inputs. If this sample does, you "
+ "can ignore this warning.");
}
for (size_t i = 0; i < jointeventinputs.size(); i++) {
double scale = double(fNEvents) / fEventHist->Integral("width");
scale *= jointeventinputs.at(i)->Integral("width");
scale /= double(jointindexhigh[i] - jointindexlow[i]);
jointindexscale.push_back(scale);
}
fEventHist->SetNameTitle((fName + "_EVT").c_str(), (fName + "_EVT").c_str());
fFluxHist->SetNameTitle((fName + "_FLUX").c_str(), (fName + "_FLUX").c_str());
// Setup Max Events
if (fMaxEvents > 1 && fMaxEvents < fNEvents) {
if (LOG_LEVEL(SAM)) {
std::cout << "\t\t|-> Read Max Entries : " << fMaxEvents << std::endl;
}
fNEvents = fMaxEvents;
}
// Print out Status
if (LOG_LEVEL(SAM)) {
std::cout << "\t\t|-> Total Entries : " << fNEvents << std::endl
<< "\t\t|-> Event Integral : "
<< fEventHist->Integral("width") * 1.E-38 << " events/nucleon"
<< std::endl
<< "\t\t|-> Flux Integral : " << fFluxHist->Integral("width")
<< " /cm2" << std::endl
<< "\t\t|-> Event/Flux : "
<< fEventHist->Integral("width") * 1.E-38 /
fFluxHist->Integral("width")
<< " cm2/nucleon" << std::endl;
}
}
-BaseFitEvt* InputHandlerBase::GetBaseEvent(const UInt_t entry) {
+BaseFitEvt *InputHandlerBase::GetBaseEvent(const UInt_t entry) {
// Do some light processing: don't calculate the kinematics
- return static_cast(GetNuisanceEvent(entry, true));
+ return static_cast(GetNuisanceEvent(entry, true));
}
double InputHandlerBase::GetInputWeight(int entry) {
- if (!jointinput) return 1.0;
+ if (!jointinput)
+ return 1.0;
// Find Switch Scale
while (entry < jointindexlow[jointindexswitch] ||
entry >= jointindexhigh[jointindexswitch]) {
jointindexswitch++;
// Loop Around
if (jointindexswitch >= jointindexlow.size()) {
jointindexswitch = 0;
}
}
return jointindexscale[jointindexswitch];
};
diff --git a/src/InputHandler/InputHandler.h b/src/InputHandler/InputHandler.h
index 1110091..c696840 100644
--- a/src/InputHandler/InputHandler.h
+++ b/src/InputHandler/InputHandler.h
@@ -1,145 +1,135 @@
// Copyright 2016 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 .
-*******************************************************************************/
+ * 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 .
+ *******************************************************************************/
#ifndef INPUTHANDLER2_H
#define INPUTHANDLER2_H
/*!
* \addtogroup InputHandler
* @{
*/
-#include "TH1D.h"
-#include "FitEvent.h"
#include "BaseFitEvt.h"
+#include "FitEvent.h"
+#include "TH1D.h"
#include "TTreePerfStats.h"
/// Base InputHandler class defining how events are requested and setup.
class InputHandlerBase {
public:
-
/// Base constructor resets everything to default
InputHandlerBase();
/// Removes flux/event rate histograms
virtual ~InputHandlerBase();
/// Return NUISANCE FitEvent Class from given event entry.
- /// Must be overriden by GeneratorInputHandler. Lightweight allows a faster option
- /// to be given where only RW information is needed.
- virtual FitEvent* GetNuisanceEvent(const UInt_t entry, const bool lightweight=false) = 0;
+ /// Must be overriden by GeneratorInputHandler. Lightweight allows a faster
+ /// option to be given where only RW information is needed.
+ virtual FitEvent *GetNuisanceEvent(const UInt_t entry,
+ const bool lightweight = false) = 0;
/// Calls GetNuisanceEvent(entry, TRUE);
- virtual BaseFitEvt* GetBaseEvent(const UInt_t entry);
-
+ virtual BaseFitEvt *GetBaseEvent(const UInt_t entry);
/// Print current event information
virtual void Print();
-
/// Return handler ID
- inline std::string GetName (void) {return fName; };
+ inline std::string GetName(void) { return fName; };
/// Return Handler Event Type Index
- inline int GetType (void) {return fEventType;};
-
+ inline int GetType(void) { return fEventType; };
/// Get Total Number of Events being Handled
- inline virtual int GetNEvents (void) {return fNEvents; };
+ inline virtual int GetNEvents(void) { return std::max(fNEvents - fSkip, 0); };
/// Get the Total Flux Histogram these events were generated with
- inline virtual TH1D* GetFluxHistogram (void) {return fFluxHist; };
+ inline virtual TH1D *GetFluxHistogram(void) { return fFluxHist; };
/// Get the Total Event Histogram these events were generated with
- inline virtual TH1D* GetEventHistogram (void) {return fEventHist;};
+ inline virtual TH1D *GetEventHistogram(void) { return fEventHist; };
/// Get the Total Cross-section Histogram (EventHist/FluxHist)
- virtual TH1D* GetXSecHistogram(void);
-
+ virtual TH1D *GetXSecHistogram(void);
/// Return all Flux Histograms for all InputFiles.
- virtual std::vector GetFluxList(void);
+ virtual std::vector GetFluxList(void);
/// Return all Event Histograms for all InputFiles
- virtual std::vector GetEventList(void);
+ virtual std::vector GetEventList(void);
/// Return all Xsec Histograms for all InputFiles
- virtual std::vector GetXSecList(void);
-
+ virtual std::vector GetXSecList(void);
/// Placeholder to create a cache to speed up reads in GeneratorInputHandler
inline virtual void CreateCache(){};
/// Placeholder to remove optional cache to free up memory
inline virtual void RemoveCache(){};
-
/// Return starting NUISANCE event pointer (entry=0)
- FitEvent* FirstNuisanceEvent();
+ FitEvent *FirstNuisanceEvent();
/// Iterate to next NUISANCE event. Returns NULL when entry > fNEvents.
- FitEvent* NextNuisanceEvent();
+ FitEvent *NextNuisanceEvent();
/// Returns starting Base Event Pointer (entry=0)
- BaseFitEvt* FirstBaseEvent();
+ BaseFitEvt *FirstBaseEvent();
/// Iterate to next NUISANCE Base Event. Returns NULL when entry > fNEvents.
- BaseFitEvt* NextBaseEvent();
-
+ BaseFitEvt *NextBaseEvent();
/// Register an input file and update event/flux information
- virtual void RegisterJointInput(std::string input, int n, TH1D* f, TH1D* e);
+ virtual void RegisterJointInput(std::string input, int n, TH1D *f, TH1D *e);
/// Finalise setup of Input event/flux information and calculate
/// joint input weights if joint input is provided.
virtual void SetupJointInputs();
/// Calculate a weight for the event given the joint input information.
/// Used to scale the relative proportion of multiple inputs correctly
/// with respect to one another.
virtual double GetInputWeight(int entry);
-
/// Returns the total predicted event rate for this input given the
/// low and high energy ranges. intOpt specifies the option the ROOT
/// TH1D integral should use. e.g. "" or "width"
double PredictedEventRate(double low = -9999.9, double high = -9999.9,
- std::string intOpt = "");
+ std::string intOpt = "");
/// Returns the total generated flux for this input given the
/// low and high energy ranges. intOpt specifies the option the ROOT
/// TH1D integral should use. e.g. "" or "width"
double TotalIntegratedFlux(double low = -9999.9, double high = -9999.9,
std::string intOpt = "");
-
/// Actual data members.
- std::vector jointfluxinputs;
- std::vector jointeventinputs;
+ std::vector jointfluxinputs;
+ std::vector jointeventinputs;
std::vector jointindexlow;
std::vector jointindexhigh;
std::vector jointindexallowed;
size_t jointindexswitch;
bool jointinput;
std::vector jointindexscale;
std::string fName;
- TH1D* fFluxHist;
- TH1D* fEventHist;
- TH1D* fXSecHist;
+ TH1D *fFluxHist;
+ TH1D *fEventHist;
+ TH1D *fXSecHist;
int fNEvents;
int fMaxEvents;
- FitEvent* fNUISANCEEvent;
- BaseFitEvt* fBaseEvent;
+ FitEvent *fNUISANCEEvent;
+ BaseFitEvt *fBaseEvent;
int fEventType;
int fCurrentIndex;
int fCacheSize;
bool kRemoveUndefParticles;
bool kRemoveFSIParticles;
bool kRemoveNuclearParticles;
- TTreePerfStats* fTTreePerformance;
-
-
+ TTreePerfStats *fTTreePerformance;
+ int fSkip;
};
/*! @} */
#endif
diff --git a/src/InputHandler/NEUTInputHandler.cxx b/src/InputHandler/NEUTInputHandler.cxx
index 6c43e00..56b41af 100644
--- a/src/InputHandler/NEUTInputHandler.cxx
+++ b/src/InputHandler/NEUTInputHandler.cxx
@@ -1,524 +1,525 @@
#ifdef __NEUT_ENABLED__
#include "NEUTInputHandler.h"
#include "InputUtils.h"
NEUTGeneratorInfo::~NEUTGeneratorInfo() { DeallocateParticleStack(); }
void NEUTGeneratorInfo::AddBranchesToTree(TTree *tn) {
tn->Branch("NEUTParticleN", fNEUTParticleN, "NEUTParticleN/I");
tn->Branch("NEUTParticleStatusCode", fNEUTParticleStatusCode,
"NEUTParticleStatusCode[NEUTParticleN]/I");
tn->Branch("NEUTParticleAliveCode", fNEUTParticleAliveCode,
"NEUTParticleAliveCode[NEUTParticleN]/I");
}
void NEUTGeneratorInfo::SetBranchesFromTree(TTree *tn) {
tn->SetBranchAddress("NEUTParticleN", &fNEUTParticleN);
tn->SetBranchAddress("NEUTParticleStatusCode", &fNEUTParticleStatusCode);
tn->SetBranchAddress("NEUTParticleAliveCode", &fNEUTParticleAliveCode);
}
void NEUTGeneratorInfo::AllocateParticleStack(int stacksize) {
fNEUTParticleN = 0;
fNEUTParticleStatusCode = new int[stacksize];
fNEUTParticleStatusCode = new int[stacksize];
}
void NEUTGeneratorInfo::DeallocateParticleStack() {
delete fNEUTParticleStatusCode;
delete fNEUTParticleAliveCode;
}
void NEUTGeneratorInfo::FillGeneratorInfo(NeutVect *nevent) {
Reset();
for (int i = 0; i < nevent->Npart(); i++) {
fNEUTParticleStatusCode[i] = nevent->PartInfo(i)->fStatus;
fNEUTParticleAliveCode[i] = nevent->PartInfo(i)->fIsAlive;
fNEUTParticleN++;
}
}
void NEUTGeneratorInfo::Reset() {
for (int i = 0; i < fNEUTParticleN; i++) {
fNEUTParticleStatusCode[i] = -1;
fNEUTParticleAliveCode[i] = 9;
}
fNEUTParticleN = 0;
}
NEUTInputHandler::NEUTInputHandler(std::string const &handle,
std::string const &rawinputs) {
NUIS_LOG(SAM, "Creating NEUTInputHandler : " << handle);
// Run a joint input handling
fName = handle;
// Setup the TChain
fNEUTTree = new TChain("neuttree");
fSaveExtra = FitPar::Config().GetParB("SaveExtraNEUT");
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector 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(inputs[inp_it].c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
NUIS_ABORT(
"NEUT 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(
(PlotUtils::GetObjectWithName(inp_file, "flux")).c_str());
TH1D *eventhist = (TH1D *)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "evt")).c_str());
if (!fluxhist or !eventhist) {
NUIS_ERR(FTL, "Input File Contents: " << inputs[inp_it]);
inp_file->ls();
NUIS_ABORT("NEUT FILE doesn't contain flux/xsec info. You may have to "
"regenerate your MC!");
}
// Get N Events
TTree *neuttree = (TTree *)inp_file->Get("neuttree");
if (!neuttree) {
NUIS_ERR(FTL, "neuttree not located in NEUT file: " << inputs[inp_it]);
NUIS_ABORT(
"Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevents = neuttree->GetEntries();
if (nevents <= 0) {
NUIS_ABORT("Trying to a TTree with "
<< 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
fNEUTTree->AddFile(inputs[inp_it].c_str());
}
// Registor all our file inputs
SetupJointInputs();
// Assign to tree
fEventType = kNEUT;
fNeutVect = NULL;
fNEUTTree->SetBranchStatus("*", false);
fNEUTTree->SetBranchStatus("vectorbranch", true);
fNEUTTree->SetBranchAddress("vectorbranch", &fNeutVect);
fNEUTTree->GetBranch("vectorbranch")->SetAutoDelete(true);
fNEUTTree->SetBranchAddress("vectorbranch", &fNeutVect);
fNEUTTree->GetEntry(0);
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetNeutVect(fNeutVect);
if (fSaveExtra) {
fNeutInfo = new NEUTGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fNeutInfo);
}
fNUISANCEEvent->HardReset();
};
NEUTInputHandler::~NEUTInputHandler(){
// if (fNEUTTree) delete fNEUTTree;
// if (fNeutVect) delete fNeutVect;
// if (fNeutInfo) delete fNeutInfo;
};
void NEUTInputHandler::CreateCache() {
if (fCacheSize > 0) {
// fNEUTTree->SetCacheEntryRange(0, fNEvents);
fNEUTTree->AddBranchToCache("vectorbranch", 1);
fNEUTTree->SetCacheSize(fCacheSize);
}
}
void NEUTInputHandler::RemoveCache() {
// fNEUTTree->SetCacheEntryRange(0, fNEvents);
fNEUTTree->AddBranchToCache("vectorbranch", 0);
fNEUTTree->SetCacheSize(0);
}
-FitEvent *NEUTInputHandler::GetNuisanceEvent(const UInt_t entry,
+FitEvent *NEUTInputHandler::GetNuisanceEvent(const UInt_t ent,
const bool lightweight) {
+ UInt_t entry = ent + fSkip;
// Catch too large entries
if (entry >= (UInt_t)fNEvents)
return NULL;
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fNEUTTree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
else {
UInt_t npart = fNeutVect->Npart();
for (size_t i = 0; i < npart; i++) {
NeutPart *part = fNUISANCEEvent->fNeutVect->PartInfo(i);
if ((part->fIsAlive == false) && (part->fStatus == -1) &&
std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4,
part->fPID)) {
fNUISANCEEvent->probe_E = part->fP.T();
fNUISANCEEvent->probe_pdg = part->fPID;
break;
} else {
continue;
}
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
// Return event pointer
return fNUISANCEEvent;
}
// From NEUT neutclass/neutpart.h
// Bool_t fIsAlive; // Particle should be tracked or not
// ( in the detector simulator )
//
// Int_t fStatus; // Status flag of this particle
// -2: Non existing particle
// -1: Initial state particle
// 0: Normal
// 1: Decayed to the other particle
// 2: Escaped from the detector
// 3: Absorped
// 4: Charge exchanged
// 5: Pauli blocked
// 6: N/A
// 7: Produced child particles
// 8: Inelastically scattered
//
int NEUTInputHandler::GetNeutParticleStatus(NeutPart *part) {
// State
int state = kUndefinedState;
// Remove Pauli blocked events, probably just single pion events
if (part->fStatus == 5) {
state = kFSIState;
// fStatus == -1 means initial state
} else if (part->fIsAlive == false && part->fStatus == -1) {
state = kInitialState;
// NEUT has a bit of a strange convention for fIsAlive and fStatus
// combinations
// for NC and neutrino particle isAlive true/false and status 2 means
// final state particle
// for other particles in NC status 2 means it's an FSI particle
// for CC it means it was an FSI particle
} else if (part->fStatus == 2) {
// NC case is a little strange... The outgoing neutrino might be alive or
// not alive. Remaining particles with status 2 are FSI particles that
// reinteracted
if (abs(fNeutVect->Mode) > 30 &&
(abs(part->fPID) == 16 || abs(part->fPID) == 14 ||
abs(part->fPID) == 12)) {
state = kFinalState;
// The usual CC case
} else if (part->fIsAlive == true) {
state = kFSIState;
}
} else if (part->fIsAlive == true && part->fStatus == 2 &&
(abs(part->fPID) == 16 || abs(part->fPID) == 14 ||
abs(part->fPID) == 12)) {
state = kFinalState;
} else if (part->fIsAlive == true && part->fStatus == 0) {
state = kFinalState;
} else if (!part->fIsAlive &&
(part->fStatus == 1 || part->fStatus == 3 || part->fStatus == 4 ||
part->fStatus == 7 || part->fStatus == 8)) {
state = kFSIState;
// There's one hyper weird case where fStatus = -3. This apparently
// corresponds to a nucleon being ejected via pion FSI when there is "data
// available"
} else if (!part->fIsAlive && (part->fStatus == -3)) {
state = kUndefinedState;
// NC neutrino outgoing
} else if (!part->fIsAlive && part->fStatus == 0 &&
(abs(part->fPID) == 16 || abs(part->fPID) == 14 ||
abs(part->fPID) == 12)) {
state = kFinalState;
// Warn if we still find alive particles without classifying them
} else if (part->fIsAlive == true) {
NUIS_ABORT("Undefined NEUT state "
<< " Alive: " << part->fIsAlive << " Status: " << part->fStatus
<< " PDG: " << part->fPID);
// Warn if we find dead particles that we haven't classified
} else {
NUIS_ABORT("Undefined NEUT state "
<< " Alive: " << part->fIsAlive << " Status: " << part->fStatus
<< " PDG: " << part->fPID);
}
return state;
}
void NEUTInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Fill Globals
fNUISANCEEvent->Mode = fNeutVect->Mode;
fNUISANCEEvent->fEventNo = fNeutVect->EventNo;
fNUISANCEEvent->fTargetA = fNeutVect->TargetA;
fNUISANCEEvent->fTargetZ = fNeutVect->TargetZ;
fNUISANCEEvent->fTargetH = fNeutVect->TargetH;
fNUISANCEEvent->fBound = bool(fNeutVect->Ibound);
if (fNUISANCEEvent->fBound) {
fNUISANCEEvent->fTargetPDG = TargetUtils::GetTargetPDGFromZA(
fNUISANCEEvent->fTargetZ, fNUISANCEEvent->fTargetA);
} else {
fNUISANCEEvent->fTargetPDG = 1000010010;
}
// Check Particle Stack
UInt_t npart = fNeutVect->Npart();
UInt_t kmax = fNUISANCEEvent->kMaxParticles;
if (npart > kmax) {
NUIS_ERR(WRN, "NEUT has too many particles. Expanding stack.");
fNUISANCEEvent->ExpandParticleStack(npart);
}
int nprimary = fNeutVect->Nprimary();
// Fill Particle Stack
for (size_t i = 0; i < npart; i++) {
// Get Current Count
int curpart = fNUISANCEEvent->fNParticles;
// Get NEUT Particle
NeutPart *part = fNeutVect->PartInfo(i);
// State
int state = GetNeutParticleStatus(part);
// Remove Undefined
if (kRemoveUndefParticles && state == kUndefinedState)
continue;
// Remove FSI
if (kRemoveFSIParticles && state == kFSIState)
continue;
// Remove Nuclear
if (kRemoveNuclearParticles &&
(state == kNuclearInitial || state == kNuclearRemnant))
continue;
// State
fNUISANCEEvent->fParticleState[curpart] = state;
// Is the paricle associated with the primary vertex?
bool primary = false;
// NEUT events are just popped onto the stack as primary, then continues to
// be non-primary
if (i < nprimary)
primary = true;
fNUISANCEEvent->fPrimaryVertex[curpart] = primary;
// Mom
fNUISANCEEvent->fParticleMom[curpart][0] = part->fP.X();
fNUISANCEEvent->fParticleMom[curpart][1] = part->fP.Y();
fNUISANCEEvent->fParticleMom[curpart][2] = part->fP.Z();
fNUISANCEEvent->fParticleMom[curpart][3] = part->fP.T();
// PDG
fNUISANCEEvent->fParticlePDG[curpart] = part->fPID;
// Add up particle count
fNUISANCEEvent->fNParticles++;
}
// Save Extra Generator Info
if (fSaveExtra) {
fNeutInfo->FillGeneratorInfo(fNeutVect);
}
// 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 NEUTUtils::FillNeutCommons(NeutVect *nvect) {
// WARNING: This has only been implemented for a neuttree and not GENIE
// This should be kept in sync with T2KNIWGUtils::GetNIWGEvent(TTree)
// NEUT version info. Can't get it to compile properly with this yet
// neutversion_.corev = nvect->COREVer;
// neutversion_.nucev = nvect->NUCEVer;
// neutversion_.nuccv = nvect->NUCCVer;
// Documentation: See nework.h
nework_.modene = nvect->Mode;
nework_.numne = nvect->Npart();
#ifdef NEUT_COMMON_QEAV
nemdls_.mdlqeaf = nvect->QEAVForm;
#else
nemdls_.mdlqeaf = nvect->QEVForm;
#endif
nemdls_.mdlqe = nvect->QEModel;
nemdls_.mdlspi = nvect->SPIModel;
nemdls_.mdldis = nvect->DISModel;
nemdls_.mdlcoh = nvect->COHModel;
neutcoh_.necohepi = nvect->COHModel;
nemdls_.xmaqe = nvect->QEMA;
nemdls_.xmvqe = nvect->QEMV;
nemdls_.kapp = nvect->KAPPA;
// nemdls_.sccfv = SCCFVdef;
// nemdls_.sccfa = SCCFAdef;
// nemdls_.fpqe = FPQEdef;
nemdls_.xmaspi = nvect->SPIMA;
nemdls_.xmvspi = nvect->SPIMV;
nemdls_.xmares = nvect->RESMA;
nemdls_.xmvres = nvect->RESMV;
neut1pi_.xmanffres = nvect->SPIMA;
neut1pi_.xmvnffres = nvect->SPIMV;
neut1pi_.xmarsres = nvect->RESMA;
neut1pi_.xmvrsres = nvect->RESMV;
neut1pi_.neiff = nvect->SPIForm;
neut1pi_.nenrtype = nvect->SPINRType;
neut1pi_.rneca5i = nvect->SPICA5I;
neut1pi_.rnebgscl = nvect->SPIBGScale;
nemdls_.xmacoh = nvect->COHMA;
nemdls_.rad0nu = nvect->COHR0;
// nemdls_.fa1coh = nvect->COHA1err;
// nemdls_.fb1coh = nvect->COHb1err;
// neutdis_.nepdf = NEPDFdef;
// neutdis_.nebodek = NEBODEKdef;
neutcard_.nefrmflg = nvect->FrmFlg;
neutcard_.nepauflg = nvect->PauFlg;
neutcard_.nenefo16 = nvect->NefO16;
neutcard_.nemodflg = nvect->ModFlg;
// neutcard_.nenefmodl = 1;
// neutcard_.nenefmodh = 1;
// neutcard_.nenefkinh = 1;
// neutpiabs_.neabspiemit = 1;
nenupr_.iformlen = nvect->FormLen;
neutpiless_.ipilessdcy = nvect->IPilessDcy;
neutpiless_.rpilessdcy = nvect->RPilessDcy;
neutpiless_.ipilessdcy = nvect->IPilessDcy;
neutpiless_.rpilessdcy = nvect->RPilessDcy;
neffpr_.fefqe = nvect->NuceffFactorPIQE;
neffpr_.fefqeh = nvect->NuceffFactorPIQEH;
neffpr_.fefinel = nvect->NuceffFactorPIInel;
neffpr_.fefabs = nvect->NuceffFactorPIAbs;
neffpr_.fefcx = nvect->NuceffFactorPICX;
neffpr_.fefcxh = nvect->NuceffFactorPICXH;
neffpr_.fefcoh = nvect->NuceffFactorPICoh;
neffpr_.fefqehf = nvect->NuceffFactorPIQEHKin;
neffpr_.fefcxhf = nvect->NuceffFactorPICXKin;
neffpr_.fefcohf = nvect->NuceffFactorPIQELKin;
for (int i = 0; i < nework_.numne; i++) {
nework_.ipne[i] = nvect->PartInfo(i)->fPID;
nework_.pne[i][0] =
(float)nvect->PartInfo(i)->fP.X() / 1000; // VC(NE)WORK in M(G)eV
nework_.pne[i][1] =
(float)nvect->PartInfo(i)->fP.Y() / 1000; // VC(NE)WORK in M(G)eV
nework_.pne[i][2] =
(float)nvect->PartInfo(i)->fP.Z() / 1000; // VC(NE)WORK in M(G)eV
}
// fsihist.h
// neutroot fills a dummy object for events with no FSI to prevent memory leak
// when
// reading the TTree, so check for it here
if ((int)nvect->NfsiVert() ==
1) { // An event with FSI must have at least two vertices
// if (nvect->NfsiPart()!=1 || nvect->Fsiprob!=-1)
// ERR(WRN) << "T2KNeutUtils::fill_neut_commons(TTree) NfsiPart!=1 or
// Fsiprob!=-1 when NfsiVert==1" << std::endl;
fsihist_.nvert = 0;
fsihist_.nvcvert = 0;
fsihist_.fsiprob = 1;
} else { // Real FSI event
fsihist_.nvert = (int)nvect->NfsiVert();
for (int ivert = 0; ivert < fsihist_.nvert; ivert++) {
fsihist_.iflgvert[ivert] = nvect->FsiVertInfo(ivert)->fVertID;
fsihist_.posvert[ivert][0] = (float)nvect->FsiVertInfo(ivert)->fPos.X();
fsihist_.posvert[ivert][1] = (float)nvect->FsiVertInfo(ivert)->fPos.Y();
fsihist_.posvert[ivert][2] = (float)nvect->FsiVertInfo(ivert)->fPos.Z();
}
fsihist_.nvcvert = nvect->NfsiPart();
for (int ip = 0; ip < fsihist_.nvcvert; ip++) {
fsihist_.abspvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomLab;
fsihist_.abstpvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomNuc;
fsihist_.ipvert[ip] = nvect->FsiPartInfo(ip)->fPID;
fsihist_.iverti[ip] = nvect->FsiPartInfo(ip)->fVertStart;
fsihist_.ivertf[ip] = nvect->FsiPartInfo(ip)->fVertEnd;
fsihist_.dirvert[ip][0] = (float)nvect->FsiPartInfo(ip)->fDir.X();
fsihist_.dirvert[ip][1] = (float)nvect->FsiPartInfo(ip)->fDir.Y();
fsihist_.dirvert[ip][2] = (float)nvect->FsiPartInfo(ip)->fDir.Z();
}
fsihist_.fsiprob = nvect->Fsiprob;
}
neutcrscom_.crsx = nvect->Crsx;
neutcrscom_.crsy = nvect->Crsy;
neutcrscom_.crsz = nvect->Crsz;
neutcrscom_.crsphi = nvect->Crsphi;
neutcrscom_.crsq2 = nvect->Crsq2;
neuttarget_.numbndn = nvect->TargetA - nvect->TargetZ;
neuttarget_.numbndp = nvect->TargetZ;
neuttarget_.numfrep = nvect->TargetH;
neuttarget_.numatom = nvect->TargetA;
posinnuc_.ibound = nvect->Ibound;
// put empty nucleon FSI history (since it is not saved in the NeutVect
// format)
// Comment out as NEUT does not have the necessary proton FSI information yet
// nucleonfsihist_.nfnvert = 0;
// nucleonfsihist_.nfnstep = 0;
}
#endif
diff --git a/src/InputHandler/NUANCEInputHandler.cxx b/src/InputHandler/NUANCEInputHandler.cxx
index b2a8668..7790b7a 100644
--- a/src/InputHandler/NUANCEInputHandler.cxx
+++ b/src/InputHandler/NUANCEInputHandler.cxx
@@ -1,937 +1,938 @@
#ifdef __NUANCE_ENABLED__
#include "NUANCEInputHandler.h"
#include "InputUtils.h"
NUANCEGeneratorInfo::~NUANCEGeneratorInfo() { DeallocateParticleStack(); }
void NUANCEGeneratorInfo::AddBranchesToTree(TTree *tn) {
// tn->Branch("NEUTParticleN", fNEUTParticleN, "NEUTParticleN/I");
// tn->Branch("NEUTParticleStatusCode", fNEUTParticleStatusCode,
// "NEUTParticleStatusCode[NEUTParticleN]/I");
// tn->Branch("NEUTParticleAliveCode", fNEUTParticleAliveCode,
// "NEUTParticleAliveCode[NEUTParticleN]/I");
}
void NUANCEGeneratorInfo::SetBranchesFromTree(TTree *tn) {
// tn->SetBranchAddress("NEUTParticleN", &fNEUTParticleN );
// tn->SetBranchAddress("NEUTParticleStatusCode", &fNEUTParticleStatusCode );
// tn->SetBranchAddress("NEUTParticleAliveCode", &fNEUTParticleAliveCode );
}
void NUANCEGeneratorInfo::AllocateParticleStack(int stacksize) {
// fNEUTParticleN = 0;
// fNEUTParticleStatusCode = new int[stacksize];
// fNEUTParticleStatusCode = new int[stacksize];
}
void NUANCEGeneratorInfo::DeallocateParticleStack() {
// delete fNEUTParticleStatusCode;
// delete fNEUTParticleAliveCode;
}
void NUANCEGeneratorInfo::FillGeneratorInfo(NuanceEvent *nevent) {
Reset();
// for (int i = 0; i < nevent->Npart(); i++) {
// fNEUTParticleStatusCode[i] = nevent->PartInfo(i)->fStatus;
// fNEUTParticleAliveCode[i] = nevent->PartInfo(i)->fIsAlive;
// fNEUTParticleN++;
// }
}
void NUANCEGeneratorInfo::Reset() {
// for (int i = 0; i < fNEUTParticleN; i++) {
// fNEUTParticleStatusCode[i] = -1;
// fNEUTParticleAliveCode[i] = 9;
// }
// fNEUTParticleN = 0;
}
NUANCEInputHandler::NUANCEInputHandler(std::string const &handle,
std::string const &rawinputs) {
NUIS_LOG(SAM, "Creating NUANCEInputHandler : " << handle);
// Run a joint input handling
fName = handle;
fSaveExtra = FitPar::Config().GetParB("SaveExtraNUANCE");
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
// Parse Inputs
std::vector inputs = InputUtils::ParseInputFileList(rawinputs);
if (inputs.size() > 1) {
NUIS_ABORT("NUANCE is not currently setup to handle joint inputs sorry!"
- << std::endl
- << "If you know how to correctly normalise the events for this"
- << " please let us know!");
+ << std::endl
+ << "If you know how to correctly normalise the events for this"
+ << " please let us know!");
}
// Read in NUANCE Tree
fNUANCETree = new TChain("h3");
fNUANCETree->AddFile(rawinputs.c_str());
// Get entries and fNuwroEvent
int nevents = fNUANCETree->GetEntries();
double EnuMin = 0.0;
double EnuMax = 1000.0;
TH1D *fluxhist = new TH1D((fName + "_FLUX").c_str(),
(fName + "_FLUX").c_str(), 100, EnuMin, EnuMax);
for (int i = 0; i < fluxhist->GetNbinsX(); i++) {
fluxhist->SetBinContent(i + 1, 1.0);
}
fluxhist->Scale(1.0 / fluxhist->Integral());
TH1D *eventhist = new TH1D((fName + "_EVT").c_str(), (fName + "_EVT").c_str(),
100, EnuMin, EnuMax);
for (int i = 0; i < fluxhist->GetNbinsX(); i++) {
eventhist->SetBinContent(i + 1, 1.0);
}
eventhist->Scale(1.0 / eventhist->Integral());
RegisterJointInput(rawinputs, nevents, fluxhist, eventhist);
SetupJointInputs();
// Setup Reader
fNuanceEvent = new NuanceEvent();
fNuanceEvent->SetBranchAddresses(fNUANCETree);
fNUANCETree->GetEntry(0);
// Setup Event in FitEvent
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetNuanceEvent(fNuanceEvent);
// Setup extra if needed
if (fSaveExtra) {
NUIS_ABORT("NO SAVEExtra Implemented for NUANCE YET!");
// fNuanceInfo = new NUANCEGeneratorInfo();
// fNUISANCEEvent->AddGeneratorInfo(fNuanceInfo);
}
};
NUANCEInputHandler::~NUANCEInputHandler() {
if (fNuanceEvent)
delete fNuanceEvent;
if (fNUANCETree)
delete fNUANCETree;
// if (fNuanceInfo) delete fNuanceInfo;
}
void NUANCEInputHandler::CreateCache() {
if (fCacheSize > 0) {
fNUANCETree->SetCacheEntryRange(0, fNEvents);
fNUANCETree->AddBranchToCache("h3", 1);
fNUANCETree->SetCacheSize(fCacheSize);
}
}
void NUANCEInputHandler::RemoveCache() {
fNUANCETree->SetCacheEntryRange(0, fNEvents);
fNUANCETree->AddBranchToCache("h3", 0);
fNUANCETree->SetCacheSize(0);
}
-FitEvent *NUANCEInputHandler::GetNuisanceEvent(const UInt_t entry,
+FitEvent *NUANCEInputHandler::GetNuisanceEvent(const UInt_t ent,
const bool lightweight) {
+ UInt_t entry = ent + fSkip;
// Check out of bounds
if (entry >= (UInt_t)fNEvents)
return NULL;
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fNUANCETree->GetEntry(entry);
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
return fNUISANCEEvent;
}
void NUANCEInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Get shortened pointer
FitEvent *evt = fNUISANCEEvent;
// Fill Global
evt->Mode = ConvertNuanceMode(fNuanceEvent);
evt->fEventNo = 0.0;
evt->fTotCrs = 1.0;
evt->fTargetA = 0.0;
evt->fTargetZ = 0.0;
evt->fTargetH = 0;
evt->fBound = 0.0;
// Fill particle Stack
evt->fNParticles = 0;
// Check Particle Stack
UInt_t npart = 2 + fNuanceEvent->n_leptons + fNuanceEvent->n_hadrons;
UInt_t kmax = evt->kMaxParticles;
if (npart > kmax) {
NUIS_ERR(FTL, "NUANCE has too many particles");
NUIS_ERR(FTL, "npart=" << npart << " kMax=" << kmax);
throw;
}
// Fill Neutrino
evt->fParticleState[0] = kInitialState;
evt->fParticleMom[0][0] = fNuanceEvent->p_neutrino[0];
evt->fParticleMom[0][1] = fNuanceEvent->p_neutrino[1];
evt->fParticleMom[0][2] = fNuanceEvent->p_neutrino[2];
evt->fParticleMom[0][3] = fNuanceEvent->p_neutrino[3];
evt->fParticlePDG[0] = fNuanceEvent->neutrino;
// Fill Target Nucleon
evt->fParticleState[1] = kInitialState;
evt->fParticleMom[1][0] = fNuanceEvent->p_targ[0];
evt->fParticleMom[1][1] = fNuanceEvent->p_targ[1];
evt->fParticleMom[1][2] = fNuanceEvent->p_targ[2];
evt->fParticleMom[1][3] = fNuanceEvent->p_targ[3];
evt->fParticlePDG[1] = fNuanceEvent->target;
evt->fNParticles = 2;
// Fill Outgoing Leptons
for (int i = 0; i < fNuanceEvent->n_leptons; i++) {
evt->fParticleState[evt->fNParticles] = kFinalState;
evt->fParticleMom[evt->fNParticles][0] = fNuanceEvent->p_lepton[i][0];
evt->fParticleMom[evt->fNParticles][1] = fNuanceEvent->p_lepton[i][1];
evt->fParticleMom[evt->fNParticles][2] = fNuanceEvent->p_lepton[i][2];
evt->fParticleMom[evt->fNParticles][3] = fNuanceEvent->p_lepton[i][3];
evt->fParticlePDG[evt->fNParticles] = fNuanceEvent->lepton[i];
evt->fNParticles++;
}
// Fill Outgoing Hadrons
for (int i = 0; i < fNuanceEvent->n_hadrons; i++) {
evt->fParticleState[evt->fNParticles] = kFinalState;
evt->fParticleMom[evt->fNParticles][0] = fNuanceEvent->p_hadron[i][0];
evt->fParticleMom[evt->fNParticles][1] = fNuanceEvent->p_hadron[i][1];
evt->fParticleMom[evt->fNParticles][2] = fNuanceEvent->p_hadron[i][2];
evt->fParticleMom[evt->fNParticles][3] = fNuanceEvent->p_hadron[i][3];
evt->fParticlePDG[evt->fNParticles] = fNuanceEvent->hadron[i];
evt->fNParticles++;
}
// Save Extra info
if (fSaveExtra) {
// fNuanceInfo->FillGeneratorInfo(fNuanceEvent);
}
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
return;
}
void NUANCEInputHandler::Print() {}
int NUANCEInputHandler::ConvertNuanceMode(NuanceEvent *evt) {
int ch = evt->channel;
int sg = 1;
if (evt->neutrino < 0)
sg = -1;
switch (ch) {
// 1 NUANCE CCQE -> NEUT CCQE 1
case 1:
return sg * 1;
// 2 NUANCE NCEL -> NEUT NCEL 51,52 -> Set from whether target is p or n
case 2:
if (evt->target == 2212)
return sg * 51;
else
return sg * 52;
// 3 NUANCE CCPIP -> NEUT CCPIP 11
case 3:
return sg * 11;
// 4 NUANCE CCPI0 -> NEUT CCPI0 = 12
case 4:
return sg * 12;
// 5 NUANCE CCPIPn -> NEUT CCPIPn 13
case 5:
return sg * 13;
// 6 NUANCE NCpPI0 -> NEUT NCpPI0 32
case 6:
return sg * 32;
// 7 NUANCE NCpPI+ -> NEUT NCpPI+ 34
case 7:
return sg * 34;
// 8 NUANCE NCnPI0 -> NEUT NCnPI0 31
case 8:
return sg * 31;
// 9 NUANCE NCnPIM -> NEUT NCnPIM 33
case 9:
return sg * 33;
// 10 NUANCE CCPIP -> NEUT CCPIP -11
case 10:
return sg * 11;
// 11 NUANCE CCPI0 -> NEUT CCPI0 -12
case 11:
return sg * 12;
// 12 NUANCE CCPIPn -> NEUT CCPIPn 13
case 12:
return sg * 13;
// 13 NUANCE NCpPI0 -> NEUT NCnPI0 -32
case 13:
return sg * 32;
// 14 NUANCE NCpPI+ -> NEUT NCpPI+ -34
case 14:
return sg * 34;
// 15 NUANCE NCnPI0 -> NEUT NCnPI0 -31
case 15:
return sg * 31;
// 16 NUANCE NCnPIM -> NEUT NCnPIM -33
case 16:
return sg * 33;
// 17 NUANCE -> NEUT 21 CC MULTIPI
case 17:
return sg * 21;
// 18 NUANCE -> NEUT 21 CC MULTIPI
case 18:
return sg * 21;
// 19 NUANCE -> NEUT 21 CC MULTIPI
case 19:
return sg * 21;
// 20 NUANCE -> NEUT 21 CC MULTIPI
case 20:
return sg * 21;
// 21 NUANCE -> NEUT 21 CC MULTIPI
case 21:
return sg * 21;
// 22 NUANCE -> NEUT 41 NC MULTIPI
case 22:
return sg * 41;
// 23 NUANCE -> NEUT 41 NC MULTIPI
case 23:
return sg * 41;
// 24 NUANCE -> NEUT 41 NC MULTIPI
case 24:
return sg * 41;
// 25 NUANCE -> NEUT 41 NC MULTIPI
case 25:
return sg * 41;
// 26 NUANCE -> NEUT 41 NC MULTIPI
case 26:
return sg * 41;
// 27 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV)
case 27:
return sg * 41;
// 28 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV)
case 28:
return sg * 21;
// 29 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV)
case 29:
return sg * 21;
// 30 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV)
case 30:
return sg * 21;
// 31 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV)
case 31:
return sg * 21;
// 32 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV)
case 32:
return sg * 21;
// 33 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV)
case 33:
return sg * 41;
// 34 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV)
case 34:
return sg * 41;
// 35 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV)
case 35:
return sg * 41;
// 36 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV)
case 36:
return sg * 41;
// 37 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV)
case 37:
return sg * 41;
// 38 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV)
case 38:
return sg * 41;
// 39 NUANCE -> NEUT 22
case 39:
return sg * 22;
// 40 NUANCE -> NEUT 22
case 40:
return sg * 22;
// 41 NUANCE -> NEUT 22
case 41:
return sg * 22;
// 42 NUANCE -> NEUT 43
case 42:
return sg * 43;
// 43 NUANCE -> NEUT 43
case 43:
return sg * 43;
// 44 NUANCE -> NUET 42
case 44:
return sg * 42;
// 45 NUANCE -> NEUT -42
case 45:
return sg * 42;
// 46 NUANCE -> NEUT -22
case 46:
return sg * 22;
// 47 NUANCE -> NEUT -22
case 47:
return sg * 22;
// 48 NUANCE -> NEUT -22
case 48:
return sg * 22;
// 49 NUANCE -> NEUT -43
case 49:
return sg * 43;
// 50 NUANCE -> NEUT -43
case 50:
return sg * 43;
// 51 NUANCE -> NEUT -42
case 51:
return sg * 42;
// 52 NUANCE -> NEUT -42
case 52:
return sg * 42;
// 53 NUANCE -> NEUT 23 CC 1K
case 53:
return sg * 23;
// 54 NUANCE -> NEUT 23 CC 1K
case 54:
return sg * 23;
// 55 NUANCE -> NEUT 23 CC 1K
case 55:
return sg * 23;
// 56 NUANCE -> NEUT 45 NC 1K
case 56:
return sg * 45;
// 57 NUANCE -> NEUT 44 NC 1K
case 57:
return sg * 44;
// 58 NUANCE -> NEUT 44 NC 1K
case 58:
return sg * 44;
// 59 NUANCE -> NEUT 44 NC 1K
case 59:
return sg * 44;
// 60 NUANCE -> NEUT -23 CC 1K
case 60:
return sg * 23;
// 61 NUANCE -> NEUT -23 CC 1K
case 61:
return sg * 23;
// 62 NUANCE -> NEUT -23 CC 1K
case 62:
return sg * 23;
// 63 NUANCE -> NEUT -23 CC 1K
case 63:
return sg * 23;
// 64 NUANCE -> NEUT -44 NC 1K
case 64:
return sg * 44;
// 65 NUANCE -> NEUT -44 NC 1K
case 65:
return sg * 44;
// 66 NUANCE -> NEUT -45 NC 1K
case 66:
return sg * 45;
// 67 NUANCE -> NEUT 22 CC1eta
case 67:
return sg * 22;
// 68 NUANCE -> NEUT 43 NC p eta
case 68:
return sg * 43;
// 69 NUANCE -> NEUT 43 NC n eta
case 69:
return sg * 43;
// 70 NUANCE -> NEUT -22 CC1eta
case 70:
return sg * 22;
// 71 NUANCE -> NEUT -43 NC p eta
case 71:
return sg * 43;
// 72 NUANCE -> NEUT 42 NC n eta
case 72:
return sg * 42;
// 73 NUANCE -> NEUT 21 CC Multi Pi
case 73:
return sg * 21;
// 74 NUANCE -> NEUT 41 NC Multi Pi
case 74:
return sg * 41;
// 75 NUANCE -> NEUT 41 NC Multi Pi
case 75:
return sg * 41;
// 76 NUANCE -> NEUT -21 CC Multi Pi
case 76:
return sg * 21;
// 77 NUANCE -> NEUT -41 NC Multi Pi
case 77:
return sg * 41;
// 78 NUANCE -> NEUT -41 NC Multi Pi
case 78:
return sg * 41;
// 79 NUANCE -> NEUT 21 CC Multi Pi
case 79:
return sg * 21;
// 80 NUANCE -> NEUT 21 CC Multi Pi
case 80:
return sg * 21;
// 81 NUANCE -> NEUT 41 NC Multi Pi
case 81:
return sg * 41;
// 82 NUANCE -> NEUT 41 NC Multi Pi
case 82:
return sg * 41;
// 83 NUANCE -> NEUT 41 NC Multi Pi
case 83:
return sg * 41;
// 84 NUANCE -> NEUT 41 NC Multi Pi
case 84:
return sg * 84;
// 85 NUANCE -> NEUT -21 CC Multi Pi
case 85:
return sg * 21;
// 86 NUANCE -> NEUT -21 CC Multi Pi
case 86:
return sg * 21;
// 87 NUANCE -> NEUT -41 CC Multi Pi
case 87:
return sg * 41;
// 88 NUANCE -> NEUT -41
case 88:
return sg * 41;
// 89 NUANCE -> NEUT -41
case 89:
return sg * 41;
// 90 NUANCE -> NEUT -41
case 90:
return sg * 41;
// 91 NUANCE -> NEUT 26 CC DIS
case 91:
return sg * 26;
// 92 NUANCE -> NEUT 46 NC DIS
case 92:
return sg * 46;
// 93 NUANCE -> NEUT 17 1#gamma from #Delta
case 93:
return sg * 17;
// 94 NUANCE -> NEUT 39 1#gamma from #Delta
case 94:
return sg * 39;
// 95 -> UNKOWN NEUT MODE
case 95:
return sg * 0;
// 96 NUANCE -> NEUT 36 NC COH
case 96:
return sg * 36;
// 97 NUANCE -> NEUT 16
case 97:
return sg * 16;
// 98 -> UNKNOWN NEUT MODE
case 98:
return sg * 0;
// 99 -> UNKNOWN NEUT MODE
case 99:
return sg * 0;
default:
NUIS_ABORT("Unknown Nuance Channel ID = " << ch);
return 0;
}
return 0;
}
/*
// Notes copied from NuanceChannels.pdf
1 NUANCE CCQE -> NEUT CCQE 1
CC, numu n --> mu- p
Cabibbo-allowed quasi-elastic scattering from nucleons
2 NUANCE NCEL -> NEUT NCEL 51,52 -> Set from whether target is p or n
NC, numu N --> num N, (N=n,p)
(quasi)-elastic scattering from nucleons
3 NUANCE CCPIP -> NEUT CCPIP 11
CC, numu p --> mu- p pi+
resonant single pion production
4 NUANCE CCPI0 -> NEUT CCPI0 = 12
CC, numu n --> mu- p pi0
resonant single pion production
5 NUANCE CCPIPn -> NEUT CCPIPn 13
CC, numu n --> mu- n pi+
resonant single pion production
6 NUANCE NCpPI0 -> NEUT NCpPI0 32
NC, numu p --> numu p pi0
resonant single pion production
7 NUANCE NCpPI+ -> NEUT NCpPI+ 34
NC, numu p --> numu n pi+
resonant single pion production
8 NUANCE NCnPI0 -> NEUT NCnPI0 31
NC, numu n --> numu n pi0
resonant single pion production
9 NUANCE NCnPIM -> NEUT NCnPIM 33
NC, numu n --> numu p pi-
resonant single pion production
10 NUANCE CCPIP -> NEUT CCPIP -11
CC, numubar p --> mu- p pi+
resonant single pion production
11 NUANCE CCPI0 -> NEUT CCPI0 -12
CC, numubar n --> mu- p pi0
resonant single pion production
12 NUANCE CCPIPn -> NEUT CCPIPn -13
CC, numubar n --> mu- n pi+
resonant single pion production
13 NUANCE NCpPI0 -> NEUT NCnPI0 -32
NC, numubar p --> numubar p pi0
resonant single pion production
14 NUANCE NCpPI+ -> NEUT NCpPI+ -34
NC, numubar p --> numubar n pi+
resonant single pion production
15 NUANCE NCnPI0 -> NEUT NCnPI0 -31
NC, numubar n --> numubar n pi0
resonant single pion production
16 NUANCE NCnPIM -> NEUT NCnPIM -33
NC, numubar n --> numubar p pi-
resonant single pion production
17 NUANCE -> NEUT 21 CC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow l^{-} N'
multi-#pi CC, numu p --> mu- Delta+ pi+ resonant processes involving more than a
single pion 18 NUANCE -> NEUT 21 CC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow
l^{-} N' multi-#pi CC, numu p --> mu- Delta++ pi0 resonant processes involving
more than a single pion 19 NUANCE -> NEUT 21 CC (1.3 < W < 2 GeV): #nu_{l} N
#rightarrow l^{-} N' multi-#pi CC, numu n --> mu- Delta+ pi0 resonant processes
involving more than a single pion 20 NUANCE -> NEUT 21 CC (1.3 < W < 2 GeV):
#nu_{l} N #rightarrow l^{-} N' multi-#pi CC, numu n --> mu- Delta0 pi+ resonant
processes involving more than a single pion 21 NUANCE -> NEUT 21 CC (1.3 < W < 2
GeV): #nu_{l} N #rightarrow l^{-} N' multi-#pi CC, numu n --> mu- Delta++ pi-
resonant processes involving more than a single pion
22 NUANCE -> NEUT 41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N
multi-#pi" NC, numu p+ --> numu Delta+ pi0 resonant processes involving more
than a single pion 23 NUANCE -> NEUT 41 "NC (1.3 < W < 2 GeV): #nu_{l} N
#rightarrow #nu_{l} N multi-#pi" NC,numu p --> numu Delta0 pi+ resonant
processes involving more than a single pion 24 NUANCE -> NEUT 41 "NC (1.3 < W <
2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi" NC, numu p --> numu Delta++
pi- resonant processes involving more than a single pion 25 NUANCE -> NEUT 41
"NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi" NC, numu n -->
numu Delta+ pi- resonant processes involving more than a single pion 26 NUANCE
-> NEUT 41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi" NC,
numu n --> numu Delta0 pi0 resonant processes involving more than a single pion
27 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N
multi-#pi" NC, numubar n --> numubar Delta- pi+ resonant processes involving
more than a single pion 28 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV): #nu_{l} N
#rightarrow l^{-} N' multi-#pi CC, numubar p --> mu- Delta+ pi+ resonant
processes involving more than a single pion 29 UANCE -> NEUT -21 CC (1.3 < W < 2
GeV): #nu_{l} N #rightarrow l^{-} N' multi-#pi CC, numubar p --> mu- Delta++ pi0
resonant processes involving more than a single pion
30 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow l^{-} N'
multi-#pi CC, numubar n --> mu- Delta+ pi0 resonant processes involving more
than a single pion 31 NUANCE -> NEUT -21 CC (1.3 < W < 2 GeV): #nu_{l} N
#rightarrow l^{-} N' multi-#pi CC, numubar n --> mu- Delta0 pi+ resonant
processes involving more than a single pion 32 NUANCE -> NEUT -21 CC (1.3 < W <
2 GeV): #nu_{l} N #rightarrow l^{-} N' multi-#pi CC, numubar n --> mu- Delta++
pi- resonant processes involving more than a single pion 33 NUANCE -> NEUT -41
"NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi" NC, numubar p+
--> numubar Delta+ pi0 resonant processes involving more than a single pion 34
NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N
multi-#pi" NC,numubar p --> numubar Delta0 pi+ resonant processes involving more
than a single pion 35 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV): #nu_{l} N
#rightarrow #nu_{l} N multi-#pi" NC, numubar p --> numubar Delta++ pi- resonant
processes involving more than a single pion 36 NUANCE -> NEUT -41 "NC (1.3 < W <
2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi" NC, numubar n --> numubar
Delta+ pi- resonant processes involving more than a single pion 37 NUANCE ->
NEUT -41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi" NC,
numubar n --> numubar Delta0 pi0 resonant processes involving more than a single
pion 38 NUANCE -> NEUT -41 "NC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow #nu_{l}
N multi-#pi" NC, numubar n --> numubar Delta- pi+ resonant processes involving
more than a single pion
// RHO Production lumped in with eta production
22 CCeta
43 NCeta on p
42 NCeta on n
39 NUANCE -> NEUT 22
CC, numu p --> mu- p rho+(770)
resonant processes involving more than a single pion
40 NUANCE -> NEUT 22
CC, numu n --> mu- p rho0(770)
resonant processes involving more than a single pion
41 NUANCE -> NEUT 22
CC, numu n --> mu- n rho+(770)
resonant processes involving more than a single pion
42 NUANCE -> NEUT 43
NC, numu p --> numu p rho0(770)
resonant processes involving more than a single pion
43 NUANCE -> NEUT 43
NC, numu p --> numu n rho+(770)
resonant processes involving more than a single pion
44 NUANCE -> NUET 42
NC, numu n --> numu n rho0(770)
resonant processes involving more than a single pion
45 NUANCE -> NEUT -42
NC, numubar n --> numubar p rho-(770)
resonant processes involving more than a single pion
46 NUANCE -> NEUT -22
CC, numubar p --> mu- p rho+(770)
resonant processes involving more than a single pion
47 NUANCE -> NEUT -22
CC, numubar n --> mu- p rho0(770)
resonant processes involving more than a single pion
48 NUANCE -> NEUT -22
CC, numubar n --> mu- n rho+(770)
resonant processes involving more than a single pion
49 NUANCE -> NEUT -43
NC, numubar p --> numubar p rho0(770)
resonant processes involving more than a single pion
50 NUANCE -> NEUT -43
NC, numubar p --> numubar n rho+(770)
resonant processes involving more than a single pion
51 NUANCE -> NEUT -42
NC, numubar n --> numubar n rho0(770)
resonant processes involving more than a single pion
52 NUANCE -> NEUT -42
NC, numubar n --> numubar p rho-(770)
resonant processes involving more than a single pion
53 NUANCE -> NEUT 23 CC 1K: #nu_{l} n #rightarrow l^{-} #Lambda K^{+}
CC, numu p --> mu- Sigma+ K+
resonant processes involving more than a single pion
54 NUANCE -> NEUT 23 CC 1K: #nu_{l} n #rightarrow l^{-} #Lambda K^{+}
CC, numu n --> mu- Sigma0 K+
resonant processes involving more than a single pion
55 NUANCE -> NEUT 23 CC 1K: #nu_{l} n #rightarrow l^{-} #Lambda K^{+}
CC, numu n --> mu- Sigma+ K0
resonant processes involving more than a single pion
56 NUANCE -> NEUT 45 NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{+}
NC, numu p --> numu Sigma0 K+
resonant processes involving more than a single pion
57 NUANCE -> NEUT 44 NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{0}
NC, numu p --> numu Sigma+ K0
resonant processes involving more than a single pion
58 NUANCE -> NEUT 44 NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{0}
NC, numu n --> numu Sigma0 K0
resonant processes involving more than a single pion
59 NUANCE -> NEUT 45 NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{+}
NC, numu n --> numu Sigma- K+
resonant processes involving more than a single pion
60 NUANCE -> NEUT -23 CC 1K: #nu_{l} n #rightarrow l^{-} #Lambda K^{+}
CC, numubar p --> mu- Sigma+ K+
resonant processes involving more than a single pion
61 NUANCE -> NEUT -23 CC 1K: #nu_{l} n #rightarrow l^{-} #Lambda K^{+}
CC, numubar n --> mu- Sigma0 K+
resonant processes involving more than a single pion
62 NUANCE -> NEUT -23 CC 1K: #nu_{l} n #rightarrow l^{-} #Lambda K^{+}
CC, numubar n --> mu- Sigma+ K0
resonant processes involving more than a single pion
63 NUANCE -> NEUT -45 NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{+}
NC, numubar p --> numubar Sigma0 K+
resonant processes involving more than a single pion
64 NUANCE -> NEUT -44 NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{0}
NC, numubar p --> numubar Sigma+ K0
resonant processes involving more than a single pion
65 NUANCE -> NEUT -44 NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{0}
NC, numubar n --> numubar Sigma0 K0
resonant processes involving more than a single pion
66 NUANCE -> NEUT -45 NC 1K: #nu_{l} n #rightarrow #nu_{l} #Lambda K^{+}
NC, numubar n --> numubar Sigma- K+
resonant processes involving more than a single pion
67 NUANCE -> NEUT 22
ModeStack[22]->SetTitle("CC1#eta^{0} on n");
CC, numu n --> mu- p eta
resonant processes involving more than a single pion
68 NUANCE -> NEUT 43
NC, numu p --> numu p eta
resonant processes involving more than a single pion
69 NUANCE -> NEUT 42
NC, numu n --> numu n eta
resonant processes involving more than a single pion
70 NUANCE -> NEUT -22
ModeStack[22]->SetTitle("CC1#eta^{0} on n");
CC, numubar n --> mu- p eta
resonant processes involving more than a single pion
71 NUANCE -> NEUT -43
ModeStack[43]->SetTitle("NC1#eta^{0} on p");
NC, numubar p --> numubar p eta
resonant processes involving more than a single pion
72 NUANCE -> NEUT -42
ModeStack[42]->SetTitle("NC1#eta^{0} on n");
NC, numubar n --> numubar n eta
resonant processes involving more than a single pion
73 NUANCE -> NEUT 21
ModeStack[21]->SetTitle("Multi #pi (1.3 < W < 2.0)");
CC, numu n --> mu- K+ Lambda
resonant processes involving more than a single pion
74 NUANCE -> NEUT 41
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC, numu p --> numu K+ Lambda
resonant processes involving more than a single pion
75 NUANCE -> NEUT 41
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC, numu n --> numu K0 Lambda
resonant processes involving more than a single pion
76 NUANCE -> NEUT -21
ModeStack[21]->SetTitle("Multi #pi (1.3 < W < 2.0)");
CC, numubar n --> mu- K+ Lambda
resonant processes involving more than a single pion
77 NUANCE -> NEUT -41
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC, numubar p --> numubar K+ Lambda
resonant processes involving more than a single pion
78 NUANCE -> NEUT -41
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC, numubar n --> numubar K0 Lambda
resonant processes involving more than a single pion
CC Multipi ModeStack[21]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC Multipi ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
79 NUANCE -> NEUT 21
CC, numu n --> mu- p pi+ pi-
two pion production
80 NUANCE -> NEUT 21
CC, numu n --> mu- p pi0 pi0
two pion production
81 NUANCE -> NEUT 41
NC, numu p --> numu p pi+ pi-
two pion production
82 NUANCE -> NEUT 41
NC, numu p --> numu p pi0 pi0
two pion production
83 NUANCE -> NEUT 41
NC, numu n --> numu n pi+ pi-
two pion production
84 NUANCE -> NEUT 41
NC, numu n --> numu n pi0 pi0
two pion production
85 NUANCE -> NEUT -21
CC, numubar n --> mu- p pi+ pi-
two pion production
86 NUANCE -> NEUT -21
CC, numubar n --> mu- p pi0 pi0
two pion production
87 NUANCE -> NEUT -41
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC, numubar p --> numubar p pi+ pi-
two pion production
88 NUANCE -> NEUT -41
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC, numubar p --> numubar p pi0 pi0
two pion production
89 NUANCE -> NEUT -41
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC, numubar n --> numubar n pi+ pi-
two pion production
90 NUANCE -> NEUT -41
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
NC, numubar n --> numubar n pi0 pi0
two pion production
91 NUANCE -> NEUT 26
ModeStack[26]->SetTitle("DIS (W > 2.0)");
CC, numu N --> mu- X (where N=n,p)
deep inelastic scattering (nu or nubar)
92 NUANCE -> NEUT 46
ModeStack[46]->SetTitle("DIS (W > 2.0)");
NC, numu N --> numu X (where N=n,p)
deep inelastic scattering (nu or nubar)
93 NUANCE -> NEUT 17 1#gamma from #Delta: #nu_{l} n #rightarrow l^{-} p #gamma
CC, numu n --> mu- p gamma
Delta radiative decay, Delta --> N gamma (only in NUANCE versions v3 and and
higher) 94 NUANCE -> NEUT 39 1#gamma from #Delta: #nu_{l} p #rightarrow #nu_{l}
p #gamma neutModeID[15] = 38; neutModeName[15] = "ncngam"; neutModeTitle[15] =
"1#gamma from #Delta: #nu_{l} n #rightarrow #nu_{l} n #gamma"; neutModeID[16] =
39; neutModeName[16] = "ncpgam"; neutModeTitle[16] = "1#gamma from #Delta:
#nu_{l} p #rightarrow #nu_{l} p #gamma"; NC, numu N --> numu N gamma Delta
radiative decay, Delta --> N gamma (only in NUANCE versions v3 and and higher)
95 -> UNKOWN NEUT MODE
CC, numubar p --> mu+ Lambda, numubar n -- > mu+ Sigma-, numubar p --> mu+
Sigma0 Cabibbo-suppressed QE hyperon production from nucleons
96 NUANCE -> NEUT 36
neutModeID[14] = 36; neutModeName[14] = "nccoh"; neutModeTitle[14] = "NC
coherent-#pi: #nu_{l} ^{16}O #rightarrow #nu_{l} ^{16}O #pi^{0}"; NC, numu A -->
numu pi0 A coherent or diffractive pi0 production 97 NUANCE -> NEUT 16
neutModeID[4] = 16; neutModeName[4] = "cccoh"; neutModeTitle[4] = "CC
coherent-#pi: #nu_{l} ^{16}O #rightarrow l^{-} ^{16}O #pi^{+}"; CC, numu A -->
mu- pi+ A (or numubar A --> coherent or diffractive pi0 production
98 -> UNKNOWN NEUT MODE
NC, numu e- --> numu e- (or numubar e- -->
neutrino + electron elastic scattering
99 -> UNKNOWN NEUT MODE
CC, numu e- --> mu- nue
neutrino + electron inverse muon decay
NEUT Modes:
// CC Modes
neutModeID[0] = 1; neutModeName[0] = "ccqe"; neutModeTitle[0] = "CCQE:
#nu_{l} n #rightarrow l^{-} p"; neutModeID[1] = 11; neutModeName[1] =
"ccppip"; neutModeTitle[1] = "CC 1#pi: #nu_{l} p #rightarrow l^{-} p #pi^{+}";
neutModeID[2] = 12; neutModeName[2] = "ccppi0"; neutModeTitle[2] = "CC 1#pi:
#nu_{l} n #rightarrow l^{-} p #pi^{0}"; neutModeID[3] = 13; neutModeName[3] =
"ccnpip"; neutModeTitle[3] = "CC 1#pi: #nu_{l} n #rightarrow l^{-} n #pi^{+}";
neutModeID[4] = 16; neutModeName[4] = "cccoh"; neutModeTitle[4] = "CC
coherent-#pi: #nu_{l} ^{16}O #rightarrow l^{-} ^{16}O #pi^{+}"; neutModeID[5] =
17; neutModeName[5] = "ccgam"; neutModeTitle[5] = "1#gamma from #Delta:
#nu_{l} n #rightarrow l^{-} p #gamma"; neutModeID[6] = 21; neutModeName[6] =
"ccmpi"; neutModeTitle[6] = "CC (1.3 < W < 2 GeV): #nu_{l} N #rightarrow l^{-}
N' multi-#pi"; neutModeID[7] = 22; neutModeName[7] = "cceta"; neutModeTitle[7]
= "CC 1#eta: #nu_{l} n #rightarrow l^{-} p #eta"; neutModeID[8] = 23;
neutModeName[8] = "cck"; neutModeTitle[8] = "CC 1K: #nu_{l} n #rightarrow
l^{-} #Lambda K^{+}"; neutModeID[9] = 26; neutModeName[9] = "ccdis";
neutModeTitle[9] = "CC DIS (2 GeV < W): #nu_{l} N #rightarrow l^{-} N' mesons";
neutModeID[10] = 31; neutModeName[10] = "ncnpi0"; neutModeTitle[10] = "NC 1#pi:
#nu_{l} n #rightarrow #nu_{l} n #pi^{0}"; neutModeID[11] = 32; neutModeName[11]
= "ncppi0"; neutModeTitle[11] = "NC 1#pi: #nu_{l} p #rightarrow #nu_{l} p
#pi^{0}"; neutModeID[12] = 33; neutModeName[12] = "ncppim"; neutModeTitle[12] =
"NC 1#pi: #nu_{l} n #rightarrow #nu_{l} p #pi^{-}"; neutModeID[13] = 34;
neutModeName[13] = "ncnpip"; neutModeTitle[13] = "NC 1#pi: #nu_{l} p #rightarrow
#nu_{l} n #pi^{+}";
neutModeID[14] = 36; neutModeName[14] = "nccoh"; neutModeTitle[14] = "NC
coherent-#pi: #nu_{l} ^{16}O #rightarrow #nu_{l} ^{16}O #pi^{0}"; neutModeID[15]
= 38; neutModeName[15] = "ncngam"; neutModeTitle[15] = "1#gamma from #Delta:
#nu_{l} n #rightarrow #nu_{l} n #gamma"; neutModeID[16] = 39; neutModeName[16]
= "ncpgam"; neutModeTitle[16] = "1#gamma from #Delta: #nu_{l} p #rightarrow
#nu_{l} p #gamma";
neutModeID[17] = 41; neutModeName[17] = "ncmpi"; neutModeTitle[17] = "NC (1.3
< W < 2 GeV): #nu_{l} N #rightarrow #nu_{l} N multi-#pi";
neutModeID[18] = 42; neutModeName[18] = "ncneta"; neutModeTitle[18] = "NC
1#eta: #nu_{l} n #rightarrow #nu_{l} n #eta"; neutModeID[19] = 43;
neutModeName[19] = "ncpeta"; neutModeTitle[19] = "NC 1#eta: #nu_{l} p
#rightarrow #nu_{l} p #eta";
neutModeID[20] = 44; neutModeName[20] = "nck0"; neutModeTitle[20] = "NC 1K:
#nu_{l} n #rightarrow #nu_{l} #Lambda K^{0}"; neutModeID[21] = 45;
neutModeName[21] = "nckp"; neutModeTitle[21] = "NC 1K: #nu_{l} n #rightarrow
#nu_{l} #Lambda K^{+}";
neutModeID[22] = 46; neutModeName[22] = "ncdis"; neutModeTitle[22] = "NC DIS
(2 GeV < W): #nu_{l} N #rightarrow #nu_{l} N' mesons";
neutModeID[23] = 51; neutModeName[23] = "ncqep"; neutModeTitle[23] = "NC
elastic: #nu_{l} p #rightarrow #nu_{l} p"; neutModeID[24] = 52; neutModeName[24]
= "ncqen"; neutModeTitle[24] = "NC elastic: #nu_{l} n #rightarrow #nu_{l} n";
*/
#endif
diff --git a/src/InputHandler/NuWroInputHandler.cxx b/src/InputHandler/NuWroInputHandler.cxx
index b6c2096..dee7c3b 100644
--- a/src/InputHandler/NuWroInputHandler.cxx
+++ b/src/InputHandler/NuWroInputHandler.cxx
@@ -1,518 +1,519 @@
#ifdef __NUWRO_ENABLED__
#include "NuWroInputHandler.h"
#include "InputUtils.h"
NuWroGeneratorInfo::~NuWroGeneratorInfo() { delete fNuWroParticlePDGs; }
void NuWroGeneratorInfo::AddBranchesToTree(TTree *tn) {
tn->Branch("NuWroParticlePDGs", &fNuWroParticlePDGs, "NuWroParticlePDGs/I");
}
void NuWroGeneratorInfo::SetBranchesFromTree(TTree *tn) {
tn->SetBranchAddress("NuWroParticlePDGs", &fNuWroParticlePDGs);
}
void NuWroGeneratorInfo::AllocateParticleStack(int stacksize) {
fNuWroParticlePDGs = new int[stacksize];
}
void NuWroGeneratorInfo::DeallocateParticleStack() {
delete fNuWroParticlePDGs;
}
void NuWroGeneratorInfo::FillGeneratorInfo(event *e) { Reset(); }
void NuWroGeneratorInfo::Reset() {
for (int i = 0; i < kMaxParticles; i++) {
fNuWroParticlePDGs[i] = 0;
}
}
int event1_nof(event *e, int pdg) {
int c = 0;
for (size_t i = 0; i < e->out.size(); i++)
if (e->out[i].pdg == pdg)
c++;
return c;
}
NuWroInputHandler::NuWroInputHandler(std::string const &handle,
std::string const &rawinputs) {
NUIS_LOG(SAM, "Creating NuWroInputHandler : " << handle);
// Run a joint input handling
fName = handle;
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
fSaveExtra = false; // FitPar::Config().GetParB("NuWroSaveExtra");
// Setup the TChain
fNuWroTree = new TChain("treeout");
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector 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(inputs[inp_it].c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
NUIS_ABORT("nuwro File IsZombie() at " << inputs[inp_it]);
}
// Get Flux/Event hist
TH1D *fluxhist = (TH1D *)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "FluxHist")).c_str());
TH1D *eventhist = (TH1D *)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "EvtHist")).c_str());
if (!fluxhist or !eventhist) {
NUIS_ERR(FTL, "nuwro FILE doesn't contain flux/xsec info");
if (FitPar::Config().GetParB("regennuwro")) {
NUIS_ERR(FTL,
- "Regen NuWro has not been added yet. Email the developers!");
+ "Regen NuWro has not been added yet. Email the developers!");
// ProcessNuWroInputFlux(inputs[inp_it]);
throw;
} else {
NUIS_ABORT("If you would like NUISANCE to generate these for you "
- << "please set parameter regennuwro=1 and re-run.");
+ << "please set parameter regennuwro=1 and re-run.");
}
}
// Get N Events
TTree *nuwrotree = (TTree *)inp_file->Get("treeout");
if (!nuwrotree) {
NUIS_ABORT("treeout not located in nuwro file! " << inputs[inp_it]);
}
int nevents = nuwrotree->GetEntries();
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add to TChain
fNuWroTree->Add(inputs[inp_it].c_str());
}
// Registor all our file inputs
SetupJointInputs();
// Setup Events
fNuWroEvent = NULL;
fNuWroTree->SetBranchAddress("e", &fNuWroEvent);
fNuWroTree->GetEntry(0);
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->fType = kNUWRO;
fNUISANCEEvent->fNuwroEvent = fNuWroEvent;
fNUISANCEEvent->HardReset();
if (fSaveExtra) {
fNuWroInfo = new NuWroGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fNuWroInfo);
}
};
NuWroInputHandler::~NuWroInputHandler() {
if (fNuWroTree)
delete fNuWroTree;
}
void NuWroInputHandler::CreateCache() {
// fNuWroTree->SetCacheEntryRange(0, fNEvents);
// fNuWroTree->AddBranchToCache("*", 1);
// fNuWroTree->SetCacheSize(fCacheSize);
}
void NuWroInputHandler::RemoveCache() {
// fNuWroTree->SetCacheEntryRange(0, fNEvents);
// fNuWroTree->AddBranchToCache("*", 0);
// fNuWroTree->SetCacheSize(0);
}
void NuWroInputHandler::ProcessNuWroInputFlux(const std::string file) {}
-FitEvent *NuWroInputHandler::GetNuisanceEvent(const UInt_t entry,
+FitEvent *NuWroInputHandler::GetNuisanceEvent(const UInt_t ent,
const bool lightweight) {
+ UInt_t entry = ent + fSkip;
// Catch too large entries
if (entry >= (UInt_t)fNEvents)
return NULL;
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fNuWroTree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
for (size_t i = 0; i < fNUISANCEEvent->fNuwroEvent->in.size(); i++) {
if (std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4,
fNUISANCEEvent->fNuwroEvent->in[i].pdg)) {
fNUISANCEEvent->probe_E = fNUISANCEEvent->fNuwroEvent->in[i].t;
fNUISANCEEvent->probe_pdg = fNUISANCEEvent->fNuwroEvent->in[i].pdg;
break;
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
#ifdef __USE_NUWRO_SRW_EVENTS__
if (!rwEvs.size()) {
fNuwroParams = fNuWroEvent->par;
}
if (entry >= rwEvs.size()) {
rwEvs.push_back(BaseFitEvt());
rwEvs.back().fType = kNUWRO;
rwEvs.back().Mode = fNUISANCEEvent->Mode;
rwEvs.back().fNuwroSRWEvent = SRW::SRWEvent(*fNuWroEvent);
rwEvs.back().fNuwroEvent = NULL;
rwEvs.back().fNuwroParams = &fNuwroParams;
rwEvs.back().probe_E = rwEvs.back().fNuwroSRWEvent.NeutrinoEnergy;
rwEvs.back().probe_pdg = rwEvs.back().fNuwroSRWEvent.NeutrinoPDG;
}
fNUISANCEEvent->fNuwroSRWEvent = SRW::SRWEvent(*fNuWroEvent);
fNUISANCEEvent->fNuwroParams = &fNuwroParams;
fNUISANCEEvent->probe_E = fNUISANCEEvent->fNuwroSRWEvent.NeutrinoEnergy;
fNUISANCEEvent->probe_pdg = fNUISANCEEvent->fNuwroSRWEvent.NeutrinoPDG;
#endif
return fNUISANCEEvent;
}
int NuWroInputHandler::ConvertNuwroMode(event *e) {
Int_t proton_pdg, neutron_pdg, pion_pdg, pion_plus_pdg, pion_minus_pdg,
lambda_pdg, eta_pdg, kaon_pdg, kaon_plus_pdg;
proton_pdg = 2212;
eta_pdg = 221;
neutron_pdg = 2112;
pion_pdg = 111;
pion_plus_pdg = 211;
pion_minus_pdg = -211;
// O_16_pdg = 100069; // oznacznie z Neuta
lambda_pdg = 3122;
kaon_pdg = 311;
kaon_plus_pdg = 321;
if (e->flag.qel) // kwiazielastyczne oddziaływanie
{
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc)
return -1;
else {
if (event1_nof(e, proton_pdg))
return -51;
else if (event1_nof(e, neutron_pdg))
return -52; // sprawdzam dodatkowo ?
}
} else // oddziaływanie z neutrinem
{
if (e->flag.cc)
return 1;
else {
if (event1_nof(e, proton_pdg))
return 51;
else if (event1_nof(e, neutron_pdg))
return 52;
}
}
}
if (e->flag.mec) {
if (e->flag.anty)
return -2;
else
return 2;
}
if (e->flag.res) // rezonansowa produkcja: pojedynczy pion, pojed.eta, kaon,
// multipiony
{
Int_t liczba_pionow, liczba_kaonow;
liczba_pionow = event1_nof(e, pion_pdg) + event1_nof(e, pion_plus_pdg) +
event1_nof(e, pion_minus_pdg);
liczba_kaonow = event1_nof(e, kaon_pdg) + event1_nof(e, kaon_pdg);
if (liczba_pionow > 1 || liczba_pionow == 0) // multipiony
{
if (e->flag.anty) {
if (e->flag.cc)
return -21;
else
return -41;
} else {
if (e->flag.cc)
return 21;
else
return 41;
}
}
if (liczba_pionow == 1) {
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc) {
if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_minus_pdg))
return -11;
if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_pdg))
return -12;
if (event1_nof(e, proton_pdg) && event1_nof(e, pion_minus_pdg))
return -13;
} else {
if (event1_nof(e, proton_pdg)) {
if (event1_nof(e, pion_minus_pdg))
return -33;
else if (event1_nof(e, pion_pdg))
return -32;
} else if (event1_nof(e, neutron_pdg)) {
if (event1_nof(e, pion_plus_pdg))
return -34;
else if (event1_nof(e, pion_pdg))
return -31;
}
}
} else // oddziaływanie z neutrinem
{
if (e->flag.cc) {
if (event1_nof(e, proton_pdg) && event1_nof(e, pion_plus_pdg))
return 11;
if (event1_nof(e, proton_pdg) && event1_nof(e, pion_pdg))
return 12;
if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_plus_pdg))
return 13;
} else {
if (event1_nof(e, proton_pdg)) {
if (event1_nof(e, pion_minus_pdg))
return 33;
else if (event1_nof(e, pion_pdg))
return 32;
} else if (event1_nof(e, neutron_pdg)) {
if (event1_nof(e, pion_plus_pdg))
return 34;
else if (event1_nof(e, pion_pdg))
return 31;
}
}
}
}
if (event1_nof(e, eta_pdg)) // produkcja rezonansowa ety
{
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc)
return -22;
else {
if (event1_nof(e, neutron_pdg))
return -42;
else if (event1_nof(e, proton_pdg))
return -43; // sprawdzam dodatkowo ?
}
} else // oddziaływanie z neutrinem
{
if (e->flag.cc)
return 22;
else {
if (event1_nof(e, neutron_pdg))
return 42;
else if (event1_nof(e, proton_pdg))
return 43;
}
}
}
if (event1_nof(e, lambda_pdg) == 1 &&
liczba_kaonow == 1) // produkcja rezonansowa kaonu
{
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc && event1_nof(e, kaon_pdg))
return -23;
else {
if (event1_nof(e, kaon_pdg))
return -44;
else if (event1_nof(e, kaon_plus_pdg))
return -45;
}
} else // oddziaływanie z neutrinem
{
if (e->flag.cc && event1_nof(e, kaon_plus_pdg))
return 23;
else {
if (event1_nof(e, kaon_pdg))
return 44;
else if (event1_nof(e, kaon_plus_pdg))
return 45;
}
}
}
}
if (e->flag.coh) // koherentne oddziaływanie tylko na O(16)
{
Int_t _target;
_target = e->par.nucleus_p + e->par.nucleus_n; // liczba masowa O(16)
if (_target == 16) {
if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem
{
if (e->flag.cc && event1_nof(e, pion_minus_pdg))
return -16;
else if (event1_nof(e, pion_pdg))
return -36;
} else // oddziaływanie z neutrinem
{
if (e->flag.cc && event1_nof(e, pion_plus_pdg))
return 16;
else if (event1_nof(e, pion_pdg))
return 36;
}
}
}
// gleboko nieelastyczne rozpraszanie
if (e->flag.dis) {
if (e->flag.anty) {
if (e->flag.cc)
return -26;
else
return -46;
} else {
if (e->flag.cc)
return 26;
else
return 46;
}
}
return 9999;
}
void NuWroInputHandler::CalcNUISANCEKinematics() {
// std::cout << "NuWro Event Address " << fNuWroEvent << std::endl;
// Reset all variables
fNUISANCEEvent->ResetEvent();
FitEvent *evt = fNUISANCEEvent;
// Sort Event Info
evt->Mode = ConvertNuwroMode(fNuWroEvent);
if (abs(evt->Mode) > 60) {
evt->Mode = 0;
}
evt->fEventNo = 0.0;
evt->fTotCrs = 0.0;
evt->fTargetA = fNuWroEvent->par.nucleus_p + fNuWroEvent->par.nucleus_n;
evt->fTargetZ = fNuWroEvent->par.nucleus_p;
evt->fTargetPDG =
TargetUtils::GetTargetPDGFromZA(evt->fTargetZ, evt->fTargetA);
evt->fTargetH = 0;
evt->fBound = (evt->fTargetA != 1);
// Check Particle Stack
UInt_t npart_in = fNuWroEvent->in.size();
UInt_t npart_out = fNuWroEvent->out.size();
UInt_t npart_post = fNuWroEvent->post.size();
UInt_t npart = npart_in + npart_out + npart_post;
UInt_t kmax = evt->kMaxParticles;
if (npart > kmax) {
NUIS_ERR(WRN, "NUWRO has too many particles. Expanding stack.");
fNUISANCEEvent->ExpandParticleStack(npart);
}
evt->fNParticles = 0;
std::vector::iterator p_iter;
// Get the Initial State
for (p_iter = fNuWroEvent->in.begin(); p_iter != fNuWroEvent->in.end();
p_iter++) {
AddNuWroParticle(fNUISANCEEvent, (*p_iter), kInitialState, true);
}
// Try to find the FSI state particles
// Loop over the primary vertex particles
// If they match the post-FSI they haven't undergone FSI.
// If they don't match post-FSI they have undergone FSI.
for (p_iter = fNuWroEvent->out.begin(); p_iter != fNuWroEvent->out.end();
p_iter++) {
// Get the particle
particle p = (*p_iter);
// Check against all the post particles, match them
std::vector::iterator p2_iter;
bool match = false;
for (p2_iter = fNuWroEvent->post.begin();
p2_iter != fNuWroEvent->post.end(); p2_iter++) {
particle p2 = (*p2_iter);
// Check energy and pdg
// A very small cascade which changes the energy by 1E-5 MeV should be
// matched
match = (fabs(p2.E() - p.E()) < 1E-5 && p2.pdg == p.pdg);
// If we match p to p2 break the loop
if (match)
break;
}
// If we've looped through the whole particle stack of post-FSI and haven't
// found a match it's a primary particle that has been FSIed
if (!match)
AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFSIState, true);
}
// Loop over the final state particles
for (p_iter = fNuWroEvent->post.begin(); p_iter != fNuWroEvent->post.end();
p_iter++) {
particle p = (*p_iter);
// To find if it's primary or not we have to loop through the primary ones
// and match, just like above
bool match = false;
std::vector::iterator p2_iter;
for (p2_iter = fNuWroEvent->out.begin(); p2_iter != fNuWroEvent->out.end();
p2_iter++) {
particle p2 = (*p2_iter);
match = (fabs(p2.E() - p.E()) < 1E-5 && p2.pdg == p.pdg);
if (match)
break;
}
AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFinalState, match);
}
// Fill Generator Info
if (fSaveExtra)
fNuWroInfo->FillGeneratorInfo(fNuWroEvent);
// 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 NuWroInputHandler::AddNuWroParticle(FitEvent *evt, particle &p, int state,
bool primary = false) {
// Add Mom
evt->fParticleMom[evt->fNParticles][0] = static_cast(p).x;
evt->fParticleMom[evt->fNParticles][1] = static_cast(p).y;
evt->fParticleMom[evt->fNParticles][2] = static_cast(p).z;
evt->fParticleMom[evt->fNParticles][3] = static_cast(p).t;
// For NuWro a particle that we've given a FSI state is a pre-FSI particle
// An initial state particle is also a primary vertex praticle
evt->fPrimaryVertex[evt->fNParticles] = primary;
// Status/PDG
evt->fParticleState[evt->fNParticles] = state;
evt->fParticlePDG[evt->fNParticles] = p.pdg;
// Add to particle count
evt->fNParticles++;
}
void NuWroInputHandler::Print() {}
#endif
diff --git a/src/InputHandler/SplineInputHandler.cxx b/src/InputHandler/SplineInputHandler.cxx
index 43b3675..cc110e8 100644
--- a/src/InputHandler/SplineInputHandler.cxx
+++ b/src/InputHandler/SplineInputHandler.cxx
@@ -1,132 +1,133 @@
#include "SplineInputHandler.h"
SplineInputHandler::SplineInputHandler(std::string const &handle,
std::string const &rawinputs) {
NUIS_LOG(SAM, "Creating SplineInputHandler : " << handle);
// Run a joint input handling
fName = handle;
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
// Setup the TChain
fFitEventTree = new TChain("nuisance_events");
// Open File for histogram access
TFile *inp_file = new TFile(rawinputs.c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
NUIS_ABORT("FitEvent File IsZombie() at " << rawinputs);
}
// Get Flux/Event hist
TH1D *fluxhist = (TH1D *)inp_file->Get("nuisance_fluxhist");
TH1D *eventhist = (TH1D *)inp_file->Get("nuisance_eventhist");
if (!fluxhist or !eventhist) {
NUIS_ABORT("FitEvent FILE doesn't contain flux/xsec info");
}
// Get N Events
TTree *eventtree = (TTree *)inp_file->Get("nuisance_events");
if (!eventtree) {
NUIS_ABORT("nuisance_events not located in FitSpline file! " << rawinputs);
}
int nevents = eventtree->GetEntries();
// Register input to form flux/event rate hists
RegisterJointInput(rawinputs, nevents, fluxhist, eventhist);
SetupJointInputs();
// Add to TChain
fFitEventTree->Add(rawinputs.c_str());
// Setup NEvents and the FitEvent
fNEvents = fFitEventTree->GetEntries();
fEventType = kSPLINEPARAMETER;
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetBranchAddress(fFitEventTree);
// Setup Spline Reader
NUIS_LOG(SAM, "Loading Spline Reader.");
fSplRead = new SplineReader();
fSplRead->Read((TTree *)inp_file->Get("spline_reader"));
fNUISANCEEvent->fSplineRead = this->fSplRead;
// Setup Matching Spline TTree
fSplTree = (TTree *)inp_file->Get("spline_tree");
fSplTree->SetBranchAddress("SplineCoeff", fSplineCoeff);
fNUISANCEEvent->fSplineCoeff = this->fSplineCoeff;
// Load into memory
for (int j = 0; j < fNEvents; j++) {
std::vector tempval;
fFitEventTree->GetEntry(j);
fStartingWeights.push_back(GetInputWeight(j));
}
};
SplineInputHandler::~SplineInputHandler() {
if (fFitEventTree)
delete fFitEventTree;
if (fSplTree)
delete fSplTree;
if (fSplRead)
delete fSplRead;
fStartingWeights.clear();
}
void SplineInputHandler::CreateCache() {
if (fCacheSize > 0) {
fFitEventTree->SetCacheEntryRange(0, fNEvents);
fFitEventTree->AddBranchToCache("*", 1);
fFitEventTree->SetCacheSize(fCacheSize);
fSplTree->SetCacheEntryRange(0, fNEvents);
fSplTree->AddBranchToCache("*", 1);
fSplTree->SetCacheSize(fCacheSize);
}
}
void SplineInputHandler::RemoveCache() {
fFitEventTree->SetCacheEntryRange(0, fNEvents);
fFitEventTree->AddBranchToCache("*", 0);
fFitEventTree->SetCacheSize(0);
fSplTree->SetCacheEntryRange(0, fNEvents);
fSplTree->AddBranchToCache("*", 0);
fSplTree->SetCacheSize(0);
}
-FitEvent *SplineInputHandler::GetNuisanceEvent(const UInt_t entry,
+FitEvent *SplineInputHandler::GetNuisanceEvent(const UInt_t ent,
const bool lightweight) {
+ UInt_t entry = ent + fSkip;
// Make sure events setup
if (entry >= (UInt_t)fNEvents)
return NULL;
// Reset all variables before tree read
fNUISANCEEvent->ResetEvent();
// Read NUISANCE Tree
if (!lightweight)
fFitEventTree->GetEntry(entry);
// Get Spline Coefficients
fSplTree->GetEntry(entry);
fNUISANCEEvent->fSplineCoeff = fSplineCoeff;
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = fStartingWeights[entry];
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
return fNUISANCEEvent;
}
double SplineInputHandler::GetInputWeight(int entry) {
double w = InputHandlerBase::GetInputWeight(entry);
return w * fNUISANCEEvent->SavedRWWeight;
}
void SplineInputHandler::Print() {}