diff --git a/cmake/setup.sh.in b/cmake/setup.sh.in
index 123aaef..44d8d3f 100644
--- a/cmake/setup.sh.in
+++ b/cmake/setup.sh.in
@@ -1,155 +1,159 @@
# 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 <http://www.gnu.org/licenses/>.
################################################################################
#!/bin/sh
### Adapted from https://unix.stackexchange.com/questions/4965/keep-duplicates-out-of-path-on-source
function add_to_PATH () {
for d; do
d=$(cd -- "$d" && { pwd -P || pwd; }) 2>/dev/null # canonicalize symbolic links
if [ -z "$d" ]; then continue; fi # skip nonexistent directory
if [ "$d" == "/usr/bin" ] || [ "$d" == "/usr/bin64" ] || [ "$d" == "/usr/local/bin" ] || [ "$d" == "/usr/local/bin64" ]; then
case ":$PATH:" in
*":$d:"*) :;;
*) export PATH=$PATH:$d;;
esac
else
case ":$PATH:" in
*":$d:"*) :;;
*) export PATH=$d:$PATH;;
esac
fi
done
}
function add_to_LD_LIBRARY_PATH () {
for d; do
d=$(cd -- "$d" && { pwd -P || pwd; }) 2>/dev/null # canonicalize symbolic links
if [ -z "$d" ]; then continue; fi # skip nonexistent directory
if [ "$d" == "/usr/lib" ] || [ "$d" == "/usr/lib64" ] || [ "$d" == "/usr/local/lib" ] || [ "$d" == "/usr/local/lib64" ]; then
case ":$LD_LIBRARY_PATH:" in
*":$d:"*) :;;
*) export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:$d;;
esac
else
case ":$LD_LIBRARY_PATH:" in
*":$d:"*) :;;
*) export LD_LIBRARY_PATH=$d:$LD_LIBRARY_PATH;;
esac
fi
done
}
if [ "@EXTRA_SETUP_SCRIPT@" ]; then
if [ ! -e @EXTRA_SETUP_SCRIPT@ ]; then
echo "[WARN]: Extra setup script \"@EXTRA_SETUP_SCRIPT@\" requested, but could not be found. Skipping..."
else
echo "[INFO]: Sourcing extra setup from \"@EXTRA_SETUP_SCRIPT@\"."
. @EXTRA_SETUP_SCRIPT@
fi
fi
add_to_PATH "@CMAKE_INSTALL_PREFIX@/bin"
add_to_LD_LIBRARY_PATH "@CMAKE_INSTALL_PREFIX@/lib"
if [ ! "${ROOTSYS}" ]; then
echo "[INFO]: Sourcing ROOT from: @CMAKE_ROOTSYS@"
source "@CMAKE_ROOTSYS@/bin/thisroot.sh"
fi
if [ "@USE_T2K@" != "FALSE" ]; then
echo "[INFO]: Adding T2K paths to the environment."
export T2KREWEIGHT=@T2KREWEIGHT@
add_to_LD_LIBRARY_PATH "@T2KREWEIGHT@/lib"
fi
if [ "@USE_NIWG@" != "FALSE" ]; then
echo "[INFO]: Adding NIWG paths to the environment."
export NIWG=@NIWG_ROOT@
export NIWGREWEIGHT_INPUTS=@NIWG_ROOT@/inputs
add_to_LD_LIBRARY_PATH "@NIWG_ROOT@"
fi
if [ "@USE_NEUT@" != "FALSE" ]; then
echo "[INFO]: Adding NEUT library paths to the environment."
export NEUT_ROOT=@NEUT_ROOT@
export CERN=@CERN@
export CERN_LEVEL=@CERN_LEVEL@
add_to_LD_LIBRARY_PATH "${NEUT_LIB_DIR}" "${NEUT_ROOT}/src/reweight"
fi
if [ "@USE_NuWro@" != "FALSE" ]; then
if [ "@NUWRO_BUILT_FROM_FILE@" == "FALSE" ]; then
echo "[INFO]: Adding NuWro library paths to the environment."
export NUWRO="@NUWRO@"
add_to_PATH "@NUWRO@/bin"
add_to_LD_LIBRARY_PATH "@NUWRO@/build/@CMAKE_SYSTEM_NAME@/lib"
if [ "@NUWRO_INC@" ]; then
export NUWRO_INC=@NUWRO_INC@
fi
else
echo "[INFO]: NuWro support included from input event file."
fi
fi
if [ "@NEED_PYTHIA6@" != "FALSE" ]; then
echo "[INFO]: Adding PYTHIA6 library paths to the environment."
export PYTHIA6="@PYTHIA6@"
add_to_LD_LIBRARY_PATH "@PYTHIA6@"
fi
if [ "@USE_GENIE@" != "FALSE" ]; then
echo "[INFO]: Adding GENIE paths to the environment."
export GENIE="@GENIE@"
export LHAPDF_LIB="@LHAPDF_LIB@"
export LHAPDF_INC="@LHAPDF_INC@"
export LIBXML2_LIB="@LIBXML2_LIB@"
export LIBXML2_INC="@LIBXML2_INC@"
export LOG4CPP_LIB="@LOG4CPP_LIB@"
export LOG4CPP_INC="@LOG4CPP_INC@"
if [ "@LHAPATH@" ]; then
export LHAPATH="@LHAPATH@"
fi
add_to_PATH "@GENIE@/bin"
add_to_LD_LIBRARY_PATH "@GENIE@/lib" "@LHAPDF_LIB@" "@LIBXML2_LIB@" "@LOG4CPP_LIB@"
fi
if [ "@BUILD_GiBUU@" != "FALSE" ]; then
echo "[INFO]: Sourcing GiBUU tools."
source @CMAKE_BINARY_DIR@/GiBUUTools/src/GiBUUTools-build/Linux/setup.sh
fi
+if [ "@NUSYST_ROOT@" != "" ]; then
+ source @NUSYST_ROOT@/setup.sh
+fi
+
export NUISANCE="@CMAKE_SOURCE_DIR@"
diff --git a/src/InputHandler/FitEvent.cxx b/src/InputHandler/FitEvent.cxx
index 7516fee..9c33652 100644
--- a/src/InputHandler/FitEvent.cxx
+++ b/src/InputHandler/FitEvent.cxx
@@ -1,466 +1,466 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is pddrt of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "FitEvent.h"
#include "TObjArray.h"
#include <iostream>
FitEvent::FitEvent() {
fGenInfo = NULL;
kRemoveFSIParticles = true;
kRemoveUndefParticles = true;
- AllocateParticleStack(25);
+ AllocateParticleStack(50);
};
void FitEvent::AddGeneratorInfo(GeneratorInfoBase *gen) {
fGenInfo = gen;
gen->AllocateParticleStack(kMaxParticles);
}
void FitEvent::AllocateParticleStack(int stacksize) {
kMaxParticles = stacksize;
if (!stacksize) {
return;
}
LOG(DEB) << "Allocating particle stack of size: " << stacksize << std::endl;
fParticleList = new FitParticle *[kMaxParticles];
fParticleMom = new double *[kMaxParticles];
fParticleState = new UInt_t[kMaxParticles];
fParticlePDG = new int[kMaxParticles];
fOrigParticleMom = new double *[kMaxParticles];
fOrigParticleState = new UInt_t[kMaxParticles];
fOrigParticlePDG = new int[kMaxParticles];
for (size_t i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
fParticleMom[i] = new double[4];
fOrigParticleMom[i] = new double[4];
}
if (fGenInfo)
fGenInfo->AllocateParticleStack(kMaxParticles);
}
void FitEvent::ExpandParticleStack(int stacksize) {
DeallocateParticleStack();
AllocateParticleStack(stacksize);
}
void FitEvent::DeallocateParticleStack() {
if (!kMaxParticles) {
return;
}
for (size_t i = 0; i < kMaxParticles; i++) {
if (fParticleList[i])
delete fParticleList[i];
delete[] fParticleMom[i];
delete[] fOrigParticleMom[i];
}
delete[] fParticleMom;
delete[] fOrigParticleMom;
delete[] fParticleList;
delete[] fParticleState;
delete[] fParticlePDG;
delete[] fOrigParticleState;
delete[] fOrigParticlePDG;
if (fGenInfo)
fGenInfo->DeallocateParticleStack();
kMaxParticles = 0;
fNParticles = 0;
}
void FitEvent::ClearFitParticles() {
for (size_t i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
}
}
void FitEvent::FreeFitParticles() {
for (size_t i = 0; i < kMaxParticles; i++) {
FitParticle *fp = fParticleList[i];
if (fp)
delete fp;
fParticleList[i] = NULL;
}
}
void FitEvent::ResetParticleList() {
for (unsigned int i = 0; i < kMaxParticles; i++) {
FitParticle *fp = fParticleList[i];
if (fp)
delete fp;
fParticleList[i] = NULL;
}
}
void FitEvent::HardReset() {
for (unsigned int i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
}
}
void FitEvent::ResetEvent() {
Mode = 9999;
fEventNo = -1;
fTotCrs = -1.0;
fTargetA = -1;
fTargetZ = -1;
fTargetH = -1;
fBound = false;
fNParticles = 0;
if (fGenInfo)
fGenInfo->Reset();
for (unsigned int i = 0; i < kMaxParticles; i++) {
if (fParticleList[i])
delete fParticleList[i];
fParticleList[i] = NULL;
continue;
fParticlePDG[i] = 0;
fParticleState[i] = kUndefinedState;
fParticleMom[i][0] = 0.0;
fParticleMom[i][1] = 0.0;
fParticleMom[i][2] = 0.0;
fParticleMom[i][3] = 0.0;
fOrigParticlePDG[i] = 0;
fOrigParticleState[i] = kUndefinedState;
fOrigParticleMom[i][0] = 0.0;
fOrigParticleMom[i][1] = 0.0;
fOrigParticleMom[i][2] = 0.0;
fOrigParticleMom[i][3] = 0.0;
}
}
void FitEvent::OrderStack() {
// Copy current stack
int npart = fNParticles;
for (int i = 0; i < npart; i++) {
fOrigParticlePDG[i] = fParticlePDG[i];
fOrigParticleState[i] = fParticleState[i];
fOrigParticleMom[i][0] = fParticleMom[i][0];
fOrigParticleMom[i][1] = fParticleMom[i][1];
fOrigParticleMom[i][2] = fParticleMom[i][2];
fOrigParticleMom[i][3] = fParticleMom[i][3];
}
// Now run loops for each particle
fNParticles = 0;
int stateorder[6] = {kInitialState, kFinalState, kFSIState,
kNuclearInitial, kNuclearRemnant, kUndefinedState};
for (int s = 0; s < 6; s++) {
for (int i = 0; i < npart; i++) {
if ((UInt_t)fOrigParticleState[i] != (UInt_t)stateorder[s])
continue;
fParticlePDG[fNParticles] = fOrigParticlePDG[i];
fParticleState[fNParticles] = fOrigParticleState[i];
fParticleMom[fNParticles][0] = fOrigParticleMom[i][0];
fParticleMom[fNParticles][1] = fOrigParticleMom[i][1];
fParticleMom[fNParticles][2] = fOrigParticleMom[i][2];
fParticleMom[fNParticles][3] = fOrigParticleMom[i][3];
fNParticles++;
}
}
if (LOG_LEVEL(DEB)) {
LOG(DEB) << "Ordered stack" << std::endl;
for (int i = 0; i < fNParticles; i++) {
LOG(DEB) << "Particle " << i << ". " << fParticlePDG[i] << " "
<< fParticleMom[i][0] << " " << fParticleMom[i][1] << " "
<< fParticleMom[i][2] << " " << fParticleMom[i][3] << " "
<< fParticleState[i] << std::endl;
}
}
if (fNParticles != npart) {
ERR(FTL) << "Dropped some particles when ordering the stack!" << std::endl;
}
return;
}
void FitEvent::Print() {
if (LOG_LEVEL(FIT)) {
LOG(FIT) << "FITEvent print" << std::endl;
LOG(FIT) << "Mode: " << Mode << ", Weight: " << InputWeight << std::endl;
LOG(FIT) << "Particles: " << fNParticles << std::endl;
LOG(FIT) << " -> Particle Stack " << std::endl;
for (int i = 0; i < fNParticles; i++) {
LOG(FIT) << " -> -> " << i << ". " << fParticlePDG[i] << " "
<< fParticleState[i] << " "
<< " Mom(" << fParticleMom[i][0] << ", " << fParticleMom[i][1]
<< ", " << fParticleMom[i][2] << ", " << fParticleMom[i][3]
<< ")." << std::endl;
}
}
return;
}
/* Read/Write own event class */
void FitEvent::SetBranchAddress(TChain *tn) {
tn->SetBranchAddress("Mode", &Mode);
tn->SetBranchAddress("EventNo", &fEventNo);
tn->SetBranchAddress("TotCrs", &fTotCrs);
tn->SetBranchAddress("TargetA", &fTargetA);
tn->SetBranchAddress("TargetH", &fTargetH);
tn->SetBranchAddress("Bound", &fBound);
tn->SetBranchAddress("RWWeight", &SavedRWWeight);
tn->SetBranchAddress("InputWeight", &InputWeight);
}
void FitEvent::AddBranchesToTree(TTree *tn) {
tn->Branch("Mode", &Mode, "Mode/I");
tn->Branch("EventNo", &fEventNo, "EventNo/i");
tn->Branch("TotCrs", &fTotCrs, "TotCrs/D");
tn->Branch("TargetA", &fTargetA, "TargetA/I");
tn->Branch("TargetH", &fTargetH, "TargetH/I");
tn->Branch("Bound", &fBound, "Bound/O");
tn->Branch("RWWeight", &RWWeight, "RWWeight/D");
tn->Branch("InputWeight", &InputWeight, "InputWeight/D");
tn->Branch("NParticles", &fNParticles, "NParticles/I");
tn->Branch("ParticleState", fOrigParticleState,
"ParticleState[NParticles]/i");
tn->Branch("ParticlePDG", fOrigParticlePDG, "ParticlePDG[NParticles]/I");
tn->Branch("ParticleMom", fOrigParticleMom, "ParticleMom[NParticles][4]/D");
}
// ------- EVENT ACCESS FUNCTION --------- //
TLorentzVector FitEvent::GetParticleP4(int index) const {
if (index == -1 or index >= fNParticles)
return TLorentzVector();
return TLorentzVector(fParticleMom[index][0], fParticleMom[index][1],
fParticleMom[index][2], fParticleMom[index][3]);
}
TVector3 FitEvent::GetParticleP3(int index) const {
if (index == -1 or index >= fNParticles)
return TVector3();
return TVector3(fParticleMom[index][0], fParticleMom[index][1],
fParticleMom[index][2]);
}
double FitEvent::GetParticleMom(int index) const {
if (index == -1 or index >= fNParticles)
return 0.0;
return sqrt(fParticleMom[index][0] * fParticleMom[index][0] +
fParticleMom[index][1] * fParticleMom[index][1] +
fParticleMom[index][2] * fParticleMom[index][2]);
}
double FitEvent::GetParticleMom2(int index) const {
if (index == -1 or index >= fNParticles)
return 0.0;
return fabs((fParticleMom[index][0] * fParticleMom[index][0] +
fParticleMom[index][1] * fParticleMom[index][1] +
fParticleMom[index][2] * fParticleMom[index][2]));
}
double FitEvent::GetParticleE(int index) const {
if (index == -1 or index >= fNParticles)
return 0.0;
return fParticleMom[index][3];
}
int FitEvent::GetParticleState(int index) const {
if (index == -1 or index >= fNParticles)
return kUndefinedState;
return (fParticleState[index]);
}
int FitEvent::GetParticlePDG(int index) const {
if (index == -1 or index >= fNParticles)
return 0;
return (fParticlePDG[index]);
}
FitParticle *FitEvent::GetParticle(int const i) {
// Check Valid Index
if (i == -1) {
return NULL;
}
// Check Valid
if (i > fNParticles) {
ERR(FTL) << "Requesting particle beyond stack!" << std::endl
<< "i = " << i << " N = " << fNParticles << std::endl
<< "Mode = " << Mode << std::endl;
throw;
}
if (!fParticleList[i]) {
/*
std::cout << "Creating particle with values i " << i << " ";
std::cout << fParticleMom[i][0] << " " << fParticleMom[i][1] << " " <<
fParticleMom[i][2] << " " << fParticleMom[i][3] << " ";
std::cout << fParticlePDG[i] << " " << fParticleState[i] << std::endl;
*/
fParticleList[i] = new FitParticle(fParticleMom[i][0], fParticleMom[i][1],
fParticleMom[i][2], fParticleMom[i][3],
fParticlePDG[i], fParticleState[i]);
} else {
/*
std::cout << "Filling particle with values i " << i << " ";
std::cout << fParticleMom[i][0] << " " << fParticleMom[i][1] << " " <<
fParticleMom[i][2] << " " << fParticleMom[i][3] << " ";
std::cout << fParticlePDG[i] << " "<< fParticleState[i] <<std::endl;
*/
fParticleList[i]->SetValues(fParticleMom[i][0], fParticleMom[i][1],
fParticleMom[i][2], fParticleMom[i][3],
fParticlePDG[i], fParticleState[i]);
}
return fParticleList[i];
}
bool FitEvent::HasParticle(int const pdg, int const state) const {
bool found = false;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 && fParticleState[i] != (uint)state)
continue;
if (fParticlePDG[i] == pdg)
found = true;
}
return found;
}
int FitEvent::NumParticle(int const pdg, int const state) const {
int nfound = 0;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 and fParticleState[i] != (uint)state)
continue;
if (pdg == 0 or fParticlePDG[i] == pdg)
nfound += 1;
}
return nfound;
}
std::vector<int> FitEvent::GetAllParticleIndices(int const pdg,
int const state) const {
std::vector<int> indexlist;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 and fParticleState[i] != (uint)state)
continue;
if (pdg == 0 or fParticlePDG[i] == pdg) {
indexlist.push_back(i);
}
}
return indexlist;
}
std::vector<FitParticle *> FitEvent::GetAllParticle(int const pdg,
int const state) {
std::vector<int> indexlist = GetAllParticleIndices(pdg, state);
std::vector<FitParticle *> plist;
for (std::vector<int>::iterator iter = indexlist.begin();
iter != indexlist.end(); iter++) {
plist.push_back(GetParticle((*iter)));
}
return plist;
}
int FitEvent::GetHMParticleIndex(int const pdg, int const state) const {
double maxmom2 = -9999999.9;
int maxind = -1;
for (int i = 0; i < fNParticles; i++) {
if (state != -1 and fParticleState[i] != (uint)state)
continue;
if (pdg == 0 or fParticlePDG[i] == pdg) {
double newmom2 = GetParticleMom2(i);
if (newmom2 > maxmom2) {
maxind = i;
maxmom2 = newmom2;
}
}
}
return maxind;
}
int FitEvent::GetBeamNeutrinoIndex(void) const {
for (int i = 0; i < fNParticles; i++) {
if (fParticleState[i] != kInitialState)
continue;
int pdg = abs(fParticlePDG[i]);
if (pdg == 12 or pdg == 14 or pdg == 16) {
return i;
}
}
return 0;
}
int FitEvent::GetBeamElectronIndex(void) const {
return GetHMISParticleIndex(11);
}
int FitEvent::GetBeamPionIndex(void) const {
return GetHMISParticleIndex(PhysConst::pdg_pions);
}
int FitEvent::NumFSMesons() {
int nMesons = 0;
for (int i = 0; i < fNParticles; i++) {
if (fParticleState[i] != kFinalState)
continue;
if (abs(fParticlePDG[i]) >= 111 && abs(fParticlePDG[i]) <= 557)
nMesons += 1;
}
return nMesons;
}
int FitEvent::NumFSLeptons(void) const {
int nLeptons = 0;
for (int i = 0; i < fNParticles; i++) {
if (fParticleState[i] != kFinalState)
continue;
if (abs(fParticlePDG[i]) == 11 || abs(fParticlePDG[i]) == 13 ||
abs(fParticlePDG[i]) == 15)
nLeptons += 1;
}
return nLeptons;
}
diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx
index 278295d..aa9a02a 100644
--- a/src/InputHandler/GENIEInputHandler.cxx
+++ b/src/InputHandler/GENIEInputHandler.cxx
@@ -1,540 +1,543 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifdef __GENIE_ENABLED__
#include "GENIEInputHandler.h"
#include "InputUtils.h"
#ifdef __DUNERWT_ENABLED__
#include "systematicstools/utility/ParameterAndProviderConfigurationUtility.hh"
#include "fhiclcpp/make_ParameterSet.h"
#endif
GENIEGeneratorInfo::~GENIEGeneratorInfo() { DeallocateParticleStack(); }
void GENIEGeneratorInfo::AddBranchesToTree(TTree *tn) {
tn->Branch("GenieParticlePDGs", &fGenieParticlePDGs, "GenieParticlePDGs/I");
}
void GENIEGeneratorInfo::SetBranchesFromTree(TTree *tn) {
tn->SetBranchAddress("GenieParticlePDGs", &fGenieParticlePDGs);
}
void GENIEGeneratorInfo::AllocateParticleStack(int stacksize) {
fGenieParticlePDGs = new int[stacksize];
}
void GENIEGeneratorInfo::DeallocateParticleStack() {
delete fGenieParticlePDGs;
}
void GENIEGeneratorInfo::FillGeneratorInfo(NtpMCEventRecord *ntpl) {
Reset();
// Check for GENIE Event
if (!ntpl)
return;
if (!ntpl->event)
return;
// Cast Event Record
GHepRecord *ghep = static_cast<GHepRecord *>(ntpl->event);
if (!ghep)
return;
// Fill Particle Stack
GHepParticle *p = 0;
TObjArrayIter iter(ghep);
// Loop over all particles
int i = 0;
while ((p = (dynamic_cast<genie::GHepParticle *>((iter).Next())))) {
if (!p)
continue;
// Get PDG
fGenieParticlePDGs[i] = p->Pdg();
i++;
}
}
void GENIEGeneratorInfo::Reset() {
for (int i = 0; i < kMaxParticles; i++) {
fGenieParticlePDGs[i] = 0;
}
}
GENIEInputHandler::GENIEInputHandler(std::string const &handle,
std::string const &rawinputs) {
LOG(SAM) << "Creating GENIEInputHandler : " << handle << std::endl;
// 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<std::string> inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
TFile *inp_file = new TFile(
InputUtils::ExpandInputDirectories(inputs[inp_it]).c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
THROW("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) {
ERROR(FTL, "Input File Contents: " << inputs[inp_it]);
inp_file->ls();
THROW("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) {
ERROR(FTL, "gtree not located in GENIE file: " << inputs[inp_it]);
THROW("Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevents = genietree->GetEntries();
if (nevents <= 0) {
THROW("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();
#ifndef __DUNERWT_ENABLED__
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetGenieEvent(fGenieNtpl);
if (fSaveExtra) {
fGenieInfo = new GENIEGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fGenieInfo);
}
fNUISANCEEvent->HardReset();
#else
std::vector<nuiskey> HandlerOpts = Config::QueryKeys("GENIEInputHandler");
fUseCache = HandlerOpts.size() && HandlerOpts.front().Has("UseCache") &&
HandlerOpts.front().GetB("UseCache");
DUNERwtCachedResponseReader = nullptr;
HaveCachedResponseReader = false;
if (fUseCache && (inputs.size() == 1)) {
std::vector<nuiskey> DuneRwtCacheParams =
Config::QueryKeys("DUNERwtResponseCache");
for (nuiskey &key : DuneRwtCacheParams) {
if (key.Has("Input") && (key.GetS("Input") == inputs.front()) &&
key.Has("CacheFile") && key.Has("ParameterFHiCL")) {
fhicl::ParameterSet ps =
fhicl::make_ParameterSet(key.GetS("ParameterFHiCL"));
fhicl::ParameterSet syst_providers = ps.get<fhicl::ParameterSet>(
"generated_systematic_provider_configuration");
systtools::param_header_map_t configuredParameterHeaders =
systtools::BuildParameterHeaders(syst_providers);
DUNERwtCachedResponseReader =
std::make_unique<systtools::PrecalculatedResponseReader<5>>(
key.GetS("CacheFile"), "resp_tree",
configuredParameterHeaders.size());
HaveCachedResponseReader = true;
break;
}
}
} else {
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetGenieEvent(fGenieNtpl);
if (fSaveExtra) {
fGenieInfo = new GENIEGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fGenieInfo);
}
fNUISANCEEvent->HardReset();
}
#endif
};
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,
const bool lightweight) {
if (entry >= (UInt_t)fNEvents)
return NULL;
// Reduce memory pressure from the cache by clearing out the last entry each
// time.
if (entry && rwEvs[entry - 1].NParticles()) {
rwEvs[entry - 1].DeallocateParticleStack();
}
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fGENIETree->GetEntry(entry);
#ifdef __DUNERWT_ENABLED__
if (entry >= rwEvs.size()) {
rwEvs.push_back(FitEvent());
if (HaveCachedResponseReader) {
rwEvs.back().DUNERwtPolyResponses =
DUNERwtCachedResponseReader->GetEventResponse(entry);
rwEvs.back().HasDUNERwtPolyResponses = true;
}
}
rwEvs[entry].SetGenieEvent(fGenieNtpl);
fNUISANCEEvent = &rwEvs[entry];
#endif
// 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 = fGenieNtpl->event;
if (!fGenieGHep)
return NULL;
TObjArrayIter iter(fGenieGHep);
genie::GHepParticle *p;
while ((p = (dynamic_cast<genie::GHepParticle *>((iter).Next())))) {
if (!p) {
continue;
}
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
if (state != genie::kIStInitialState) {
continue;
}
fNUISANCEEvent->probe_E = p->E() * 1.E3;
fNUISANCEEvent->probe_pdg = p->Pdg();
break;
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
return fNUISANCEEvent;
}
int GENIEInputHandler::GetGENIEParticleStatus(genie::GHepParticle *p,
int mode) {
/*
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 {
ERROR(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;
// 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 = fGenieNtpl->event;
if (!fGenieGHep)
return;
// Convert GENIE Reaction Code
fNUISANCEEvent->Mode = ConvertGENIEReactionCode(fGenieGHep);
// Set Event Info
fNUISANCEEvent->fEventNo = 0.0;
fNUISANCEEvent->fTotCrs = fGenieGHep->XSec();
fNUISANCEEvent->fTargetA = 0.0;
fNUISANCEEvent->fTargetZ = 0.0;
fNUISANCEEvent->fTargetH = 0;
fNUISANCEEvent->fBound = 0.0;
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) {
fNUISANCEEvent->ExpandParticleStack(npart);
}
// Fill Particle Stack
GHepParticle *p = 0;
TObjArrayIter iter(fGenieGHep);
fNUISANCEEvent->fNParticles = 0;
// Loop over all particles
while ((p = (dynamic_cast<genie::GHepParticle *>((iter).Next())))) {
if (!p)
continue;
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
// Remove Undefined
if (kRemoveUndefParticles && state == kUndefinedState)
continue;
// Remove FSI
if (kRemoveFSIParticles && state == kFSIState)
continue;
if (kRemoveNuclearParticles &&
(state == kNuclearInitial || state == kNuclearRemnant))
continue;
// Fill Vectors
int curpart = fNUISANCEEvent->fNParticles;
fNUISANCEEvent->fParticleState[curpart] = state;
// Mom
fNUISANCEEvent->fParticleMom[curpart][0] = p->Px() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][1] = p->Py() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][2] = p->Pz() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][3] = p->E() * 1.E3;
// PDG
fNUISANCEEvent->fParticlePDG[curpart] = p->Pdg();
// Add to N particle count
fNUISANCEEvent->fNParticles++;
// Extra Check incase GENIE fails.
if ((UInt_t)fNUISANCEEvent->fNParticles == kmax) {
- ERR(WRN) << "Number of GENIE Particles exceeds maximum!" << std::endl;
+ ERR(WRN) << "Number of GENIE Particles exceeds maximum (Max: " << kmax
+ << ", GHEP: " << fGenieGHep->GetEntries()
+ << ", Added: " << fNUISANCEEvent->fNParticles << ")!"
+ << std::endl;
ERR(WRN) << "Extend kMax, or run without including FSI particles!"
<< std::endl;
break;
}
}
// Fill Extra Stack
if (fSaveExtra)
fGenieInfo->FillGeneratorInfo(fGenieNtpl);
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
FitParticle *ISNeutralLepton =
fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos);
if (ISNeutralLepton) {
fNUISANCEEvent->probe_E = ISNeutralLepton->E();
fNUISANCEEvent->probe_pdg = ISNeutralLepton->PDG();
}
return;
}
void GENIEInputHandler::Print() {}
#endif
diff --git a/src/Reweight/DUNERwtWeightEngine.cxx b/src/Reweight/DUNERwtWeightEngine.cxx
index 3477010..00f5d3d 100644
--- a/src/Reweight/DUNERwtWeightEngine.cxx
+++ b/src/Reweight/DUNERwtWeightEngine.cxx
@@ -1,167 +1,164 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "DUNERwtWeightEngine.h"
#include <limits>
#include <string>
DUNERwtWeightEngine::DUNERwtWeightEngine() { Config(); }
void DUNERwtWeightEngine::Config() {
std::vector<nuiskey> DuneRwtParam = Config::QueryKeys("DUNERwt");
if (DuneRwtParam.size() < 1) {
ERROR(WRN, "Oscillation parameters specified but no OscParam element "
"configuring the experimental characteristics found.\nExpect at "
"least <OscParam baseline_km=\"XXX\" />. Pausing for "
"10...");
sleep(10);
return;
}
std::string fhicl_name = DuneRwtParam.front().GetS("ConfigFHiCL");
DUNErwt.LoadConfiguration(fhicl_name);
}
int DUNERwtWeightEngine::ConvDial(std::string name) {
if (!DUNErwt.HaveHeader(name)) {
THROW("DUNERwtWeightEngine passed dial: "
<< name << " that it does not understand.");
}
return DUNErwt.GetHeaderId(name);
}
void DUNERwtWeightEngine::IncludeDial(std::string name, double startval) {
EnabledParams.push_back({systtools::paramId_t(ConvDial(name)), startval});
}
void DUNERwtWeightEngine::SetDialValue(int nuisenum, double val) {
systtools::paramId_t DuneRwtEnum = (nuisenum % 1000);
systtools::ParamValue &pval =
GetParamElementFromContainer(EnabledParams, DuneRwtEnum);
fHasChanged = (pval.val - val) > std::numeric_limits<double>::epsilon();
pval.val = val;
}
void DUNERwtWeightEngine::SetDialValue(std::string name, double val) {
if (!IsDialIncluded(name)) {
THROW("DUNERwtWeightEngine passed dial: " << name
<< " that is not enabled.");
}
systtools::ParamValue &pval =
GetParamElementFromContainer(EnabledParams, ConvDial(name));
fHasChanged = (pval.val - val) > std::numeric_limits<double>::epsilon();
pval.val = val;
}
bool DUNERwtWeightEngine::IsDialIncluded(std::string name) {
return IsDialIncluded(ConvDial(name));
}
bool DUNERwtWeightEngine::IsDialIncluded(int nuisenum) {
systtools::paramId_t DuneRwtEnum = (nuisenum % 1000);
return systtools::ContainterHasParam(EnabledParams, DuneRwtEnum);
}
double DUNERwtWeightEngine::GetDialValue(std::string name) {
if (!IsDialIncluded(name)) {
THROW("DUNERwtWeightEngine passed dial: " << name
<< " that is not enabled.");
}
systtools::ParamValue &pval =
GetParamElementFromContainer(EnabledParams, ConvDial(name));
return pval.val;
}
double DUNERwtWeightEngine::GetDialValue(int nuisenum) {
if (!IsDialIncluded(nuisenum)) {
THROW("DUNERwtWeightEngine passed dial: " << nuisenum
<< " that is not enabled.");
}
systtools::paramId_t DuneRwtEnum = (nuisenum % 1000);
systtools::ParamValue &pval =
GetParamElementFromContainer(EnabledParams, DuneRwtEnum);
return pval.val;
}
void DUNERwtWeightEngine::Reconfigure(bool silent) { fHasChanged = false; };
bool DUNERwtWeightEngine::NeedsEventReWeight() {
if (fHasChanged) {
return true;
}
return false;
}
double DUNERwtWeightEngine::CalcWeight(BaseFitEvt *evt) {
double weight = 1;
if (evt->HasDUNERwtPolyResponses) {
for (size_t i = 0; i < EnabledParams.size(); ++i) {
if (DUNErwt.IsSplineParam(EnabledParams[i].pid)) {
if (!ContainterHasParam(evt->DUNERwtPolyResponses,
EnabledParams[i].pid)) {
continue;
}
weight *= GetParamElementFromContainer(evt->DUNERwtPolyResponses,
EnabledParams[i].pid)
.resp.eval(EnabledParams[i].val);
- evt->DUNERwtResponses =
- DUNErwt.GetEventResponses(*evt->genie_event->event);
-
} else {
if (!evt->HasDUNERwtResponses) {
evt->DUNERwtResponses =
DUNErwt.GetEventResponses(*evt->genie_event->event);
evt->HasDUNERwtResponses = true;
}
weight *= DUNErwt.GetDiscreteResponse(EnabledParams[i].pid,
int(EnabledParams[i].val),
evt->DUNERwtResponses);
}
}
} else {
if (!evt->HasDUNERwtResponses) {
evt->DUNERwtResponses =
DUNErwt.GetEventResponses(*evt->genie_event->event);
evt->HasDUNERwtResponses = true;
}
for (size_t i = 0; i < EnabledParams.size(); ++i) {
if (DUNErwt.IsSplineParam(EnabledParams[i].pid)) {
weight *= DUNErwt.GetParameterResponse(
EnabledParams[i].pid, EnabledParams[i].val, evt->DUNERwtResponses);
} else {
weight *= DUNErwt.GetDiscreteResponse(EnabledParams[i].pid,
int(EnabledParams[i].val),
evt->DUNERwtResponses);
}
}
}
return weight;
}
void DUNERwtWeightEngine::Print() {
std::cout << "DUNERwtWeightEngine: " << std::endl;
}