diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx
index 60ef26f..7c9c12d 100644
--- a/src/InputHandler/GENIEInputHandler.cxx
+++ b/src/InputHandler/GENIEInputHandler.cxx
@@ -1,636 +1,635 @@
// Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#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;
}
}
bool GENIEInputHandler::IsPrimary(GHepParticle *p) {
// If initial state nucleon, or nucleon target, it is definintely primary!
if (p->Status() == genie::kIStInitialState ||
p->Status() == genie::kIStNucleonTarget) return true;
// Reject intermediate resonant state or pre-DIS state
if (p->Status() == genie::kIStDISPreFragmHadronicState || // DIS before fragmentation
p->Status() == genie::kIStPreDecayResonantState || // pre decay resonance state
p->Status() == genie::kIStIntermediateState || // intermediate state
p->Status() == genie::kIStDecayedState || // Decayed state
p->Status() == genie::kIStUndefined) return false; // undefined state
// Check if the mother is the neutrino or IS nucleon
if (p->FirstMother() < 2) return true;
// We've now filtered off intermediate states and obvious initial states
// Loop over this particle's mothers, grandmothers, and so on cleaning out intermediate particles
// Set the starting particle to be our particle
GHepParticle *mother = fGenieGHep->Particle(p->FirstMother());
while (mother->FirstMother() > 1) {
// It could be that the mother's status is actually a decayed state linked back to the vertex
if (mother->Status() == genie::kIStDecayedState || // A decayed state
mother->Status() == genie::kIStDISPreFragmHadronicState || // A DIS state before fragmentation
mother->Status() == genie::kIStPreDecayResonantState ) { // A pre-decay resonant state
mother = fGenieGHep->Particle(mother->FirstMother());
} else { // If not, move out of the loop
break;
}
}
// Then do a simple check of mother is associated with the primary
int MotherID = mother->FirstMother();
if (MotherID > 2) return false;
// Finally, this should mean that our partcile is marked for transport through the nucleus
// Could also be interactions of free proton
if (p->Status() == genie::kIStHadronInTheNucleus || // Then require the particle to be paseed to FSI
(p->Status() == genie::kIStStableFinalState && // Can also have interaction on free proton
fGenieGHep->Summary()->InitState().TgtPtr()->A() == 1 &&
fGenieGHep->Summary()->InitState().TgtPtr()->Z() == 1) ) {
return true;
}
return false;
}
GENIEInputHandler::GENIEInputHandler(std::string const &handle,
std::string const &rawinputs) {
NUIS_LOG(SAM, "Creating GENIEInputHandler : " << handle);
// Plz no shouting
StopTalking();
genie::Messenger::Instance()->SetPriorityLevel("GHepUtils", pFATAL);
StartTalking();
// Shout all you want
// Run a joint input handling
fName = handle;
// Setup the TChain
fGENIETree = new TChain("gtree");
fSaveExtra = FitPar::Config().GetParB("SaveExtraGenie");
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
// 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]);
}
// 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");
}
// Get N Events
TTree *genietree = (TTree *)inp_file->Get("gtree");
if (!genietree) {
NUIS_ERR(FTL, "gtree not located in GENIE file: " << inputs[inp_it]);
NUIS_ABORT(
"Check your inputs, they may need to be completely regenerated!");
}
int nevents = genietree->GetEntries();
if (nevents <= 0) {
NUIS_ABORT("Trying to a TTree with "
<< nevents << " to TChain from : " << inputs[inp_it]);
}
// 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());
}
// Registor all our file inputs
SetupJointInputs();
// Assign to tree
fEventType = kGENIE;
fGenieNtpl = NULL;
fGENIETree->SetBranchAddress("gmcrec", &fGenieNtpl);
// Libraries should be seen but not heard...
StopTalking();
fGENIETree->GetEntry(0);
StartTalking();
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetGenieEvent(fGenieNtpl);
if (fSaveExtra) {
fGenieInfo = new GENIEGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fGenieInfo);
}
fNUISANCEEvent->HardReset();
};
GENIEInputHandler::~GENIEInputHandler() {
// if (fGenieGHep) delete fGenieGHep;
// if (fGenieNtpl) delete fGenieNtpl;
// if (fGENIETree) delete fGENIETree;
// if (fGenieInfo) delete fGenieInfo;
}
void GENIEInputHandler::CreateCache() {
if (fCacheSize > 0) {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 1);
fGENIETree->SetCacheSize(fCacheSize);
}
}
void GENIEInputHandler::RemoveCache() {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 0);
fGENIETree->SetCacheSize(0);
}
FitEvent *GENIEInputHandler::GetNuisanceEvent(const UInt_t ent,
const bool lightweight) {
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 (kRemoveNuclearParticles) {
if (p->Pdg() > 1000000 &&
pdg::IsNeutronOrProton(
fGenieGHep->Summary()->InitState().TgtPtr()->HitNucPdg())) { // Check that it isn't coherent! Want to keep initial state if so
if (state == kInitialState)
state = kNuclearInitial;
else if (state == kFinalState)
state = kNuclearRemnant;
}
}
return state;
}
#endif
#ifdef __GENIE_ENABLED__
int GENIEInputHandler::ConvertGENIEReactionCode(GHepRecord *gheprec) {
// I randomly picked 53 here because NEUT doesn't have an appropriate mode...
if (gheprec->Summary()->ProcInfo().IsNuElectronElastic()){
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg())) return 53;
else return -53;
}
// And the same story for 54
- if (gheprec->Summary()->ProcInfo().IsIMDAnnihilation()){
+ if (gheprec->Summary()->ProcInfo().IsIMDAnnihilation() || gheprec->Summary()->ProcInfo().IsInverseMuDecay()){
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg())) return 54;
else return -54;
}
-
// Electron Scattering
if (gheprec->Summary()->ProcInfo().IsEM()) {
if (gheprec->Summary()->InitState().ProbePdg() == 11) {
if (gheprec->Summary()->ProcInfo().IsQuasiElastic())
return 1;
else if (gheprec->Summary()->ProcInfo().IsMEC())
return 2;
else if (gheprec->Summary()->ProcInfo().IsResonant())
return 13;
else if (gheprec->Summary()->ProcInfo().IsDeepInelastic())
return 26;
else {
NUIS_ERR(WRN,
"Unknown GENIE Electron Scattering Mode!"
<< std::endl
<< "ScatteringTypeId = "
<< gheprec->Summary()->ProcInfo().ScatteringTypeId() << " "
<< "InteractionTypeId = "
<< gheprec->Summary()->ProcInfo().InteractionTypeId()
<< std::endl
<< genie::ScatteringType::AsString(
gheprec->Summary()->ProcInfo().ScatteringTypeId())
<< " "
<< genie::InteractionType::AsString(
gheprec->Summary()->ProcInfo().InteractionTypeId())
<< " " << gheprec->Summary()->ProcInfo().IsMEC());
return 0;
}
}
// Weak CC
} else if (gheprec->Summary()->ProcInfo().IsWeakCC()) {
// CC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 2;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -2;
#ifndef GENIE_PRE_R3
} else if (gheprec->Summary()->ProcInfo().IsDiffractive()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 15;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -15;
#endif
// CC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
// Weak NC
} else if (gheprec->Summary()->ProcInfo().IsWeakNC()) {
// NC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 32;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -32;
#ifndef GENIE_PRE_R3
} else if (gheprec->Summary()->ProcInfo().IsDiffractive()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 35;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -35;
#endif
// NC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
}
return 0;
}
void GENIEInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Check for GENIE Event
if (!fGenieNtpl)
return;
if (!fGenieNtpl->event)
return;
// Cast Event Record
fGenieGHep = static_cast(fGenieNtpl->event);
if (!fGenieGHep)
return;
// Convert GENIE Reaction Code
fNUISANCEEvent->Mode = ConvertGENIEReactionCode(fGenieGHep);
if (!fNUISANCEEvent->Mode) {
NUIS_ERR(WRN, "Failed to determine mode for GENIE event: ");
std::cout << *fGenieGHep << std::endl;
}
// Set Event Info
fNUISANCEEvent->fEventNo = 0.0;
fNUISANCEEvent->fTotCrs = fGenieGHep->XSec();
// Have a bool storing if interaction happened on free or bound nucleon
bool IsFree = false;
// Set the TargetPDG
if (fGenieGHep->TargetNucleus() != NULL) {
fNUISANCEEvent->fTargetPDG = fGenieGHep->TargetNucleus()->Pdg();
IsFree = false;
// Sometimes GENIE scatters off free nucleons, electrons, photons
// In which TargetNucleus is NULL and we need to find the initial state
// particle
} else {
// Check the particle is an initial state particle
// Follows GHepRecord::TargetNucleusPosition but doesn't do check on
// pdg::IsIon
GHepParticle *p = fGenieGHep->Particle(1);
// Check that particle 1 actually exists
if (!p) {
NUIS_ABORT("Can't find particle 1 for GHepRecord");
}
// If not an ion but is an initial state particle
if (!pdg::IsIon(p->Pdg()) && p->Status() == kIStInitialState) {
IsFree = true;
fNUISANCEEvent->fTargetPDG = p->Pdg();
// Catch if something strange happens:
// Here particle 1 is not an initial state particle OR
// particle 1 is an ion OR
// both
} else {
if (pdg::IsIon(p->Pdg())) {
NUIS_ABORT(
"Particle 1 in GHepRecord stack is an ion but isn't an initial "
"state particle");
} else {
NUIS_ABORT(
"Particle 1 in GHepRecord stack is not an ion but is an initial "
"state particle");
}
}
}
// Set the A and Z and H from the target PDG
// Depends on if we scattered off a free or bound nucleon
if (!IsFree) {
fNUISANCEEvent->fTargetA =
TargetUtils::GetTargetAFromPDG(fNUISANCEEvent->fTargetPDG);
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();
// 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] = IsPrimary(p);
// Add to N particle count
fNUISANCEEvent->fNParticles++;
// Extra Check incase GENIE fails.
if ((UInt_t)fNUISANCEEvent->fNParticles == kmax) {
NUIS_ERR(WRN, "Number of GENIE Particles exceeds maximum!");
NUIS_ERR(WRN, "Extend kMax, or run without including FSI particles!");
break;
}
}
// Fill Extra Stack
if (fSaveExtra)
fGenieInfo->FillGeneratorInfo(fGenieNtpl);
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
FitParticle *ISAnyLepton = fNUISANCEEvent->GetHMISAnyLeptons();
if (ISAnyLepton) {
fNUISANCEEvent->probe_E = ISAnyLepton->E();
fNUISANCEEvent->probe_pdg = ISAnyLepton->PDG();
}
return;
}
void GENIEInputHandler::Print() {}
#endif
diff --git a/src/InputHandler/GiBUUNativeInputHandler.cxx b/src/InputHandler/GiBUUNativeInputHandler.cxx
index 4a20db1..f5c7335 100644
--- a/src/InputHandler/GiBUUNativeInputHandler.cxx
+++ b/src/InputHandler/GiBUUNativeInputHandler.cxx
@@ -1,494 +1,511 @@
#ifdef __GiBUU_ENABLED__
#include "GiBUUNativeInputHandler.h"
#include "InputUtils.h"
#include "PhysConst.h"
bool GiBUUEventReader::SetBranchAddresses(TChain* tn){
tn->SetBranchAddress("weight", &weight);
tn->SetBranchAddress("evType", &mode);
tn->SetBranchAddress("process_ID", &process_ID);
tn->SetBranchAddress("flavor_ID", &flavor_ID);
tn->SetBranchAddress("numEnsembles", &num_ensembles);
tn->SetBranchAddress("numRuns", &num_runs);
tn->SetBranchAddress("nucleus_A", &nucleus_A);
tn->SetBranchAddress("nucleus_Z", &nucleus_Z);
tn->SetBranchAddress("barcode", &pdg, &b_pdg);
tn->SetBranchAddress("Px", &Px, &b_Px);
tn->SetBranchAddress("Py", &Py, &b_Py);
tn->SetBranchAddress("Pz", &Pz, &b_Pz);
tn->SetBranchAddress("E", &E, &b_E);
tn->SetBranchAddress("x", &x, &b_x);
tn->SetBranchAddress("y", &y, &b_y);
tn->SetBranchAddress("z", &z, &b_z);
tn->SetBranchAddress("lepIn_Px", &lepIn_Px);
tn->SetBranchAddress("lepIn_Py", &lepIn_Py);
tn->SetBranchAddress("lepIn_Pz", &lepIn_Pz);
tn->SetBranchAddress("lepIn_E", &lepIn_E);
tn->SetBranchAddress("lepOut_Px", &lepOut_Px);
tn->SetBranchAddress("lepOut_Py", &lepOut_Py);
tn->SetBranchAddress("lepOut_Pz", &lepOut_Pz);
tn->SetBranchAddress("lepOut_E", &lepOut_E);
tn->SetBranchAddress("nuc_Px", &nuc_Px);
tn->SetBranchAddress("nuc_Py", &nuc_Py);
tn->SetBranchAddress("nuc_Pz", &nuc_Pz);
tn->SetBranchAddress("nuc_E", &nuc_E);
tn->SetBranchAddress("nuc_charge", &nuc_chrg);
return 0;
}
int ConvertModeGiBUUtoNEUT(int gibuu_mode, int process_ID, int struck_nucleon_pdg, int first_part_pdg){
bool IsCC = (abs(process_ID) == 2) ? true : false;
bool IsNu = (process_ID > 0) ? true : false;
switch(gibuu_mode){
// QE/elastic
case 1:
if (IsCC) {
return (IsNu ? 1 : -1);
} else {
return (IsNu ? 1 : -1) * ((struck_nucleon_pdg == 2212) ? 51 : 52);
}
// Different resonances
case 2:
case 3:
case 4:
case 5:
case 6:
case 7:
case 8:
case 9:
case 10:
case 11:
case 12:
case 13:
case 14:
case 15:
case 16:
case 17:
case 18:
case 19:
case 20:
case 21:
case 22:
case 23:
case 24:
case 25:
case 26:
case 27:
case 28:
case 29:
case 30:
case 31: {
if (IsCC){
if (IsNu){
if (struck_nucleon_pdg == 2212) return 11;
if (first_part_pdg == 111) return 12;
return 13;
} else {
if (struck_nucleon_pdg == 2112) return -11;
if (first_part_pdg == 111) return -12;
return -13;
}
} else {
if (struck_nucleon_pdg == 2212) {
if (first_part_pdg == 111) return 32 * (IsNu ? 1 : -1);
return 34 * (IsNu ? 1 : -1);
} else {
if (first_part_pdg == 111) return 31 * (IsNu ? 1 : -1);
return 33 * (IsNu ? 1 : -1);
}
}
}
case 32:
case 33: { // 1Pi Bkg
return (IsNu ? 1 : -1) * (10 + 20 * (!IsCC));
}
case 34: { // DIS
return (IsNu ? 1 : -1) * (26 + 20 * (!IsCC));
}
case 35:
case 36: { // MEC/2p-2h
return (IsNu ? 1 : -1) * (2 + 40 * (!IsCC));
}
case 37: { // MultiPi
return (IsNu ? 1 : -1) * (21 + 20 * (!IsCC));
}
default:
NUIS_ERR(WRN, "Unable to map GiBUU code " << gibuu_mode << " to a NEUT mode");
}
return 0;
}
int GetGiBUUNuPDG(int flavor_ID, int process_ID){
int pdg = 0;
switch(flavor_ID){
case 1:
pdg = 12;
break;
case 2:
pdg = 14;
break;
case 3:
pdg = 16;
break;
}
if (process_ID < 0) pdg*=-1;
return pdg;
}
int GetGiBUULepPDG(int flavor_ID, int process_ID){
int nuPDG = GetGiBUUNuPDG(flavor_ID, process_ID);
// Check for EM
if (abs(process_ID) == 1){
NUIS_ABORT("GiBUU file includes electron scattering events, not currently supported!");
} else if (abs(process_ID) == 3){
return nuPDG;
} else if (abs(process_ID) == 2){
return (nuPDG > 0) ? nuPDG - 1 : nuPDG + 1;
}
NUIS_ABORT("Unknown GiBUU process_ID " << process_ID);
}
// Check whether the particle is on-shell or not
int CheckGiBUUParticleStatus(double E, int pdg, double dist){
double onshell_mass = PhysConst::GetMass(pdg);
// Hard code some other values used in GiBUU...
if (pdg == 2212 || pdg == 2112) onshell_mass = 0.938;
if (abs(pdg) == 211 || pdg == 111) onshell_mass = 0.138;
// If the particle is still within the nucleus, it's not final state
// Not that this depends on the nucleus, and too few time steps could cause an issue
// 6 fm is Ulrich's guess for Argon (and will be fine for smaller nuclei)
if (dist < 6) {
return kFSIState;
}
-
+
// Check for unknown particles and default to on shell
if (onshell_mass == -1) return kFinalState;
-
+
// Ulrich's second criteria
if (E < onshell_mass) return kFSIState;
return kFinalState;
}
GiBUUNativeInputHandler::~GiBUUNativeInputHandler() {
jointtype .clear();
jointrequested.clear();
}
GiBUUNativeInputHandler::GiBUUNativeInputHandler(std::string const &handle,
std::string const &rawinputs) {
NUIS_LOG(SAM, "Creating GiBUUNativeInputHandler : " << handle);
// Run a joint input handling
fName = handle;
fEventType = kGiBUU;
fGiBUUTree = new TChain("RootTuple");
// 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!");
}
// 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, "evtrt").c_str());
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 "
"use PrepareGiBUU!");
}
// Get N Events
TChain *tempChain = new TChain("RootTuple");
tempChain->Add(inputs[inp_it].c_str());
if (!tempChain){
NUIS_ERR(FTL, "RootTuple tree missing from GiBUU file "
<< inputs[inp_it]);
NUIS_ABORT(
"Check your inputs, they may need to be completely regenerated!");
}
int nevents = tempChain->GetEntries();
if (nevents <= 0) {
NUIS_ABORT("Trying to a TTree with "
<< nevents << " to TChain from : " << inputs[inp_it]);
}
// Get specific information about the config from the file
GiBUUEventReader *tempReader = new GiBUUEventReader();
tempReader->SetBranchAddresses(tempChain);
tempChain ->GetEntry(0);
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], tempReader->process_ID, tempReader->flavor_ID,
tempReader->nucleus_A, nevents,
tempReader->nucleus_A*tempReader->num_runs*tempReader->num_ensembles,
fluxhist, eventhist);
// Add To TChain
fGiBUUTree->AddFile(inputs[inp_it].c_str());
delete tempChain;
delete tempReader;
}
// Registor all our file inputs
SetupJointInputs();
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fGiReader = new GiBUUEventReader();
fGiReader->SetBranchAddresses(fGiBUUTree);
fNUISANCEEvent->HardReset();
};
FitEvent *GiBUUNativeInputHandler::GetNuisanceEvent(const UInt_t ent,
const bool lightweight) {
UInt_t entry = ent + fSkip;
// Check out of bounds
if (entry >= (UInt_t)fNEvents)
return NULL;
fGiBUUTree->GetEntry(entry);
fNUISANCEEvent->ResetEvent();
fNUISANCEEvent->fEventNo = entry;
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
else {
fNUISANCEEvent->probe_E = fGiReader->lepIn_E*1E3;
fNUISANCEEvent->probe_pdg = GetGiBUUNuPDG(fGiReader->flavor_ID, fGiReader->process_ID);
}
#endif
fNUISANCEEvent->InputWeight *= GetInputWeight(entry);
return fNUISANCEEvent;
}
void GiBUUNativeInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
// fNUISANCEEvent->ResetEvent();
FitEvent *evt = fNUISANCEEvent;
evt->Mode = ConvertModeGiBUUtoNEUT(fGiReader->mode, fGiReader->process_ID,
(fGiReader->nuc_chrg) ? 2212 : 2112,
(*fGiReader->pdg)[0]);
// evt->fEventNo = 0.0;
evt->fTotCrs = 0;
evt->fTargetA = fGiReader->nucleus_A;
evt->fTargetZ = fGiReader->nucleus_Z;
evt->fTargetH = 0;
evt->fBound = 0.0;
// Extra GiBUU Input Weight
evt->InputWeight = fGiReader->weight;
// Check number of particles in the stack. Note the additional 3 particles!
uint npart = fGiReader->pdg->size();
int kmax = evt->kMaxParticles;
if ((UInt_t)npart+3 > (UInt_t)kmax) {
NUIS_ERR(WRN, "GiBUU has too many particles (" << npart << ") Expanding Stack.");
fNUISANCEEvent->ExpandParticleStack(npart+2);
}
// Stuff the leptons in
evt->fParticleState[0] = kInitialState;
evt->fParticleMom[0][0] = fGiReader->lepIn_Px * 1E3;
evt->fParticleMom[0][1] = fGiReader->lepIn_Py * 1E3;
evt->fParticleMom[0][2] = fGiReader->lepIn_Pz * 1E3;
evt->fParticleMom[0][3] = fGiReader->lepIn_E * 1E3;
evt->fParticlePDG[0] = GetGiBUUNuPDG(fGiReader->flavor_ID, fGiReader->process_ID);
evt->fParticleState[1] = kFinalState;
evt->fParticleMom[1][0] = fGiReader->lepOut_Px * 1E3;
evt->fParticleMom[1][1] = fGiReader->lepOut_Py * 1E3;
evt->fParticleMom[1][2] = fGiReader->lepOut_Pz * 1E3;
evt->fParticleMom[1][3] = fGiReader->lepOut_E * 1E3;
evt->fParticlePDG[1] = GetGiBUULepPDG(fGiReader->flavor_ID, fGiReader->process_ID);
// Also the struck nucleon (note only one!)
evt->fParticleState[2] = kInitialState;
evt->fParticleMom[2][0] = fGiReader->nuc_Px * 1E3;
evt->fParticleMom[2][1] = fGiReader->nuc_Py * 1E3;
evt->fParticleMom[2][2] = fGiReader->nuc_Pz * 1E3;
evt->fParticleMom[2][3] = fGiReader->nuc_E * 1E3;
evt->fParticlePDG[2] = (fGiReader->nuc_chrg) ? 2212 : 2112;
// Create Stack
evt->fNParticles = 3;
for (uint i = 0; i < npart; i++) {
int curpart = evt->fNParticles;
double dist = sqrt((*fGiReader->x)[i]*(*fGiReader->x)[i] + (*fGiReader->y)[i]*(*fGiReader->y)[i] + (*fGiReader->z)[i]*(*fGiReader->z)[i]);
// Set State
evt->fParticleState[curpart] = CheckGiBUUParticleStatus((*fGiReader->E)[i], (*fGiReader->pdg)[i], dist);
-
// Mom
evt->fParticleMom[curpart][0] = (*fGiReader->Px)[i] * 1E3;
evt->fParticleMom[curpart][1] = (*fGiReader->Py)[i] * 1E3;
evt->fParticleMom[curpart][2] = (*fGiReader->Pz)[i] * 1E3;
evt->fParticleMom[curpart][3] = (*fGiReader->E)[i] * 1E3;
// PDG
evt->fParticlePDG[curpart] = (*fGiReader->pdg)[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 GiBUUNativeInputHandler::Print() {}
void GiBUUNativeInputHandler::RegisterJointInput(std::string input, int process_ID,
int flavor_ID, int nnucleons, int n,
int nrequested, 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());
// This has to correspond to the real number of events
jointindexlow .push_back(fNEvents);
jointindexhigh.push_back(fNEvents + n);
fNEvents += n;
// This should be nnucleons*num_ensembles*num_runs (the number of requested events)
jointrequested .push_back(nrequested);
// To keep track of how many of each type there are
jointtype .push_back(nnucleons*1000+process_ID*10 + flavor_ID);
}
void GiBUUNativeInputHandler::SetupJointInputs() {
// Always apply the scaling
jointinput = true;
if (jointeventinputs.size() > 1) {
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!");
}
- // Loop over the samples
- for (size_t i = 0; i < jointeventinputs.size(); i++) {
+ // Need to know what the total number of nucleons is for correct normalization...
+ int total_unique_nucl = 0;
+ std::vector unique_nucl;
+ std::vector nsame_vect;
- // Assume all of the fluxes are the same... is this safe?
+ for (size_t i = 0; i < jointeventinputs.size(); i++) {
+
+ // Assume all of the fluxes are the same
if (!fFluxHist) fFluxHist = (TH1D*)jointfluxinputs[i]->Clone();
- // For this sample, we need to get the total event rate and number of requested events
+ // Get the list of unique nuclei nuclei
+ int this_nucl = jointtype[i]/1000;
+ if (std::find(unique_nucl.begin(), unique_nucl.end(), this_nucl) == unique_nucl.end()){
+ unique_nucl.push_back(this_nucl);
+ total_unique_nucl += this_nucl;
+ }
+
+ // Get the total event rate and number of requested events
TH1D *this_evt_total = NULL;
- int requested_total = 0;
int nsame = 0;
+ // How many files were set up in the same way?
for (size_t j = 0; j < jointeventinputs.size(); j++) {
- // Only consider the same type
if (jointtype[j] != jointtype[i]) continue;
-
nsame++;
-
- // Get the totals
- requested_total += jointrequested[j];
-
+
if (!this_evt_total) this_evt_total = (TH1D*)jointeventinputs[j]->Clone();
else this_evt_total->Add(jointeventinputs[j]);
}
+ nsame_vect.push_back(nsame);
- this_evt_total->Scale(1./double(nsame));
- double scale = double(jointrequested[i]) / this_evt_total->Integral("width")/ double(requested_total);
- scale *= jointfluxinputs.at(i)->Integral("width");
- jointindexscale.push_back(scale);
-
- // Keep track of the total event rate
+ // Need to average events with the same run config, and add others
+ this_evt_total->Scale(jointtype[i]/1000/double(nsame));
if (!fEventHist) fEventHist = (TH1D*)this_evt_total->Clone();
else fEventHist ->Add(this_evt_total);
delete this_evt_total;
}
+ // Get the correct / nucleon event rate
+ fEventHist->Scale(1/double(total_unique_nucl));
+
+ // Loop over the samples to get the event scaling correct
+ for (size_t i = 0; i < jointeventinputs.size(); i++) {
+
+ // This simply reverses the usual scaling to get to a flux averaged XSEC used everywhere else in NUISANCE
+ double scale = fFluxHist->Integral("width")*fNEvents/fEventHist->Integral("width");
+
+ // Need to correctly weight by the number of nucleons to arrive at the /nucleon XSEC for a compound target
+ scale *= jointtype[i]/1000/double(total_unique_nucl)/nsame_vect[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;
}
}
#endif