diff --git a/cmake/GENIESetup.cmake b/cmake/GENIESetup.cmake
index df26b47..4139f82 100644
--- a/cmake/GENIESetup.cmake
+++ b/cmake/GENIESetup.cmake
@@ -1,246 +1,249 @@
# 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 .
################################################################################
# TODO
# check system for libxml2
# check whether we need the includes
# check if we can use a subset of the GENIE libraries
include(${CMAKE_SOURCE_DIR}/cmake/parseConfigApp.cmake)
################################################################################
# Check Dependencies
################################################################################
################################# GENIE ######################################
if(GENIE STREQUAL "")
cmessage(FATAL_ERROR "Variable GENIE is not defined. "
"The location of a pre-built GENIE install must be defined either as"
" $ cmake -DGENIE=/path/to/GENIE or as an environment variable"
" $ export GENIE=/path/to/GENIE")
endif()
execute_process(COMMAND genie-config --version
OUTPUT_VARIABLE GENIE_VER OUTPUT_STRIP_TRAILING_WHITESPACE)
cmessage(STATUS "genie_ver: ${GENIE_VER}")
if(GENIE_VER VERSION_GREATER 3.0.0)
set(GENIE_POST_R3 1)
string(REPLACE "." "" GENIE_VERSION ${GENIE_VER})
cmessage(STATUS "set genie_post_r3")
endif()
if(NOT GENIE_POST_R3)
LIST(APPEND EXTRA_CXX_FLAGS -DGENIE_PRE_R3)
cmessage(STATUS "setting genie_pre_r3 ${EXTRA_CXX_FLAGS}")
if(GENIE_EMPMEC_REWEIGHT)
cmessage(STATUS "Enable EMPMEC dials")
LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_EMP_MECRW_ENABLED)
endif()
else()
cmessage(STATUS "Enable EMPMEC dials")
LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_EMP_MECRW_ENABLED)
+ if(DEFINED USE_GENIE_XSECMEC AND USE_GENIE_XSECMEC)
+ LIST(APPEND EXTRA_CXX_FLAGS -DUSE_GENIE_XSECMEC)
+ endif()
endif()
execute_process (COMMAND genie-config
--topsrcdir OUTPUT_VARIABLE GENIE_INCLUDES_DIR OUTPUT_STRIP_TRAILING_WHITESPACE)
#Allows for external override in the case where genie-config lies.
if(NOT DEFINED GENIE_LIB_DIR OR GENIE_LIB_DIR STREQUAL "")
#This looks like it should call libdir, but it strips the argument with -L from the response of --libs
GETLIBDIRS(genie-config --libs GENIE_LIB_DIR)
endif()
GETLIBS(genie-config --libs GENIE_LIBS)
cmessage(STATUS "GENIE version : ${GENIE_VERSION}")
cmessage(STATUS "GENIE libdir : ${GENIE_LIB_DIR}")
cmessage(STATUS "GENIE libs : ${GENIE_LIBS}")
string(REGEX MATCH "ReinSeghal" WASMATCHED ${GENIE_LIBS})
if(WASMATCHED AND GENIE_VERSION STREQUAL "210")
set(GENIE_SEHGAL ${GENIE_LIBS})
STRING(REPLACE "ReinSeghal" "ReinSehgal" GENIE_LIBS ${GENIE_SEHGAL})
cmessage(DEBUG "Fixed inconsistency in library naming: ${GENIE_LIBS}")
endif()
if(NOT USE_REWEIGHT)
SET(USING_GENIE_RW FALSE)
elseif(NOT GENIE_POST_R3)
LIST(FIND GENIE_LIBS GReWeight FOUND_GENIE_RW)
if(FOUND_GENIE_RW EQUAL -1)
cmessage(DEBUG "Did NOT find ReWeight library. Here are libs: ${GENIE_LIBS}")
SET(USING_GENIE_RW FALSE)
else()
SET(USING_GENIE_RW TRUE)
endif()
elseif(DEFINED GENIE_REWEIGHT AND NOT GENIE_REWEIGHT STREQUAL "")
LIST(FIND GENIE_LIBS GRwFwk FOUND_GENIE_RW)
if(FOUND_GENIE_RW EQUAL -1)
LIST(APPEND GENIE_LIBS GRwClc GRwFwk GRwIO)
cmessage(DEBUG "Force added ReWeight library. Here are libs: ${GENIE_LIBS}")
SET(USING_GENIE_RW TRUE)
else()
SET(USING_GENIE_RW FALSE)
endif()
endif()
if(USING_GENIE_RW)
cmessage(STATUS "Using GENIE ReWeight library.")
else()
cmessage(STATUS "Building without GENIE ReWeight support.")
endif()
LIST(APPEND GENIE_LIBS -Wl,--end-group )
LIST(REVERSE GENIE_LIBS)
LIST(APPEND GENIE_LIBS -Wl,--start-group -Wl,--no-as-needed )
LIST(REVERSE GENIE_LIBS)
cmessage(DEBUG "GENIE_LIBS: ${GENIE_LIBS}")
################################ LHAPDF ######################################
if(LHAPDF_LIB STREQUAL "")
cmessage(FATAL_ERROR "Variable LHAPDF_LIB is not defined. The location of a pre-built lhapdf install must be defined either as $ cmake -DLHAPDF_LIB=/path/to/LHAPDF_libraries or as an environment variable $ export LHAPDF_LIB=/path/to/LHAPDF_libraries")
endif()
if(LHAPDF_INC STREQUAL "")
cmessage(FATAL_ERROR "Variable LHAPDF_INC is not defined. The location of a pre-built lhapdf install must be defined either as $ cmake -DLHAPDF_INC=/path/to/LHAPDF_includes or as an environment variable $ export LHAPDF_INC=/path/to/LHAPDF_includes")
endif()
if(LHAPATH STREQUAL "")
cmessage(FATAL_ERROR "Variable LHAPATH is not defined. The location of a the LHAPATH directory must be defined either as $ cmake -DLHAPATH=/path/to/LHAPATH or as an environment variable $ export LHAPATH=/path/to/LHAPATH")
endif()
################################ LIBXML ######################################
if(LIBXML2_LIB STREQUAL "")
GETLIBDIR(xml2-config --libs LIBXML2_LIB IGNORE_EMPTY_RESPONSE)
if(LIBXML2_LIB STREQUAL "")
message(WARNING "Variable LIBXML2_LIB is not defined, as xml2-config was found and didn't report a library include path, it is likely that libxml2.so can be found in the standard system location, lets hope so. Alternativly, a location can be forced by configering with -DLIBXML2_LIB=/path/to/LIBXML2_libraries or as an environment variable LIBXML2_LIB.")
endif()
endif()
if(LIBXML2_INC STREQUAL "")
GETINCDIR(xml2-config --cflags LIBXML2_INC IGNORE_EMPTY_RESPONSE)
if(LIBXML2_INC STREQUAL "")
message(WARNING "Variable LIBXML2_INC is not defined, as xml2-config was found and didn't report an include path, it is likely that libxml2.so can be found in the standard system location, lets hope so. Alternativly, a location can be forced by configering with -DLIBXML2_INC=/path/to/LIBXML2_includes or as an environment variable LIBXML2_INC.")
endif()
endif()
############################### log4cpp ######################################
if(LOG4CPP_LIB STREQUAL "")
GETLIBDIR(log4cpp-config --libs LOG4CPP_LIB IGNORE_EMPTY_RESPONSE)
if(LOG4CPP_LIB STREQUAL "")
message(WARNING "Variable LOG4CPP_LIB is not defined, as xml2-config was found and didn't report a library include path, it is likely that liblog4cpp.so can be found in the standard system location, lets hope so. Alternativly, a location can be forced by configering with -DLOG4CPP_LIB=/path/to/LOG4CPP_libraries or as an environment variable LOG4CPP_LIB.")
endif()
endif()
if(LOG4CPP_INC STREQUAL "")
GETINCDIR(log4cpp-config --cflags LOG4CPP_INC IGNORE_EMPTY_RESPONSE)
if(LOG4CPP_INC STREQUAL "")
message(WARNING "Variable LOG4CPP_LIB is not defined, as xml2-config was found and didn't report an include path, it is likely that log4cpp headers can be found in the standard system location, lets hope so. Alternativly, a location can be forced by configering with -DLOG4CPP_INC=/path/to/LOG4CPP_includes or as an environment variable LOG4CPP_INC.")
endif()
endif()
################################################################################
# Set the compiler defines
LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_ENABLED__ -D__GENIE_VERSION__=${GENIE_VERSION})
LIST(APPEND EXTRA_LIBS ${GENIE_LIBS})
############################### GSL ######################################
if(GENIE_POST_R3)
if(GSL_LIB STREQUAL "")
GETLIBDIR(gsl-config --libs GSL_LIB)
if(GSL_LIB STREQUAL "")
message(FATAL_ERROR "Variable GSL_LIB is not defined and could not be found with gsl-config. The location of a pre-built gsl install must be defined either as $ cmake -DGSL_LIB=/path/to/GSL_libraries or as an environment variable $ export GSL_LIB=/path/to/GSL_libraries")
endif()
endif()
if(GSL_INC STREQUAL "")
GETINCDIR(gsl-config --cflags GSL_INC)
if(GSL_INC STREQUAL "")
message(FATAL_ERROR "Variable GSL_INC is not defined and could not be found with gsl-config. The location of a pre-built gsl install must be defined either as $ cmake -DGSL_INC=/path/to/GSL_includes or as an environment variable $ export GSL_INC=/path/to/GSL_includes")
endif()
endif()
GETLIBS(gsl-config --libs GSL_LIB_LIST)
if(USING_GENIE_RW AND GENIE_REWEIGHT STREQUAL "")
message(FATAL_ERROR "Variable GENIE_REWEIGHT is not defined. When using GENIE v3+, we require the reweight product to be built and accessible via the environment variable GENIE_REWEIGHT")
endif()
endif()
################################################################################
LIST(APPEND EXTRA_LIBS LHAPDF xml2 log4cpp)
LIST(APPEND EXTRA_LINK_DIRS
${GENIE_LIB_DIR}
${LHAPDF_LIB}
${LOG4CPP_LIB})
# Append only if we have found GENIE ReWeight
if(NOT GENIE_POST_R3)
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES
${GENIE_INCLUDES_DIR}
${GENIE_INCLUDES_DIR}/GHEP
${GENIE_INCLUDES_DIR}/Ntuple)
if(USING_GENIE_RW)
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${GENIE_INCLUDES_DIR}/ReWeight)
endif()
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES
${GENIE_INCLUDES_DIR}/Apps
${GENIE_INCLUDES_DIR}/FluxDrivers
${GENIE_INCLUDES_DIR}/EVGDrivers
${LHAPDF_INC}
${LIBXML2_INC}
${LOG4CPP_INC})
else()
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES
${GENIE_INCLUDES_DIR})
if(USING_GENIE_RW)
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${GENIE_REWEIGHT}/src)
endif()
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${GSL_INC}
${LHAPDF_INC}
${LIBXML2_INC}
${LOG4CPP_INC})
if(USING_GENIE_RW)
LIST(APPEND EXTRA_LINK_DIRS
${GENIE_REWEIGHT}/lib)
endif()
LIST(APPEND EXTRA_LINK_DIRS
${GSL_LIB}
)
LIST(APPEND EXTRA_LIBS ${GSL_LIB_LIST})
endif()
if(USE_PYTHIA8)
set(NEED_PYTHIA8 TRUE)
set(NEED_ROOTPYTHIA8 TRUE)
else()
set(NEED_PYTHIA6 TRUE)
set(NEED_ROOTPYTHIA6 TRUE)
endif()
set(NEED_ROOTEVEGEN FALSE)
SET(USE_GENIE TRUE CACHE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. " FORCE)
diff --git a/src/InputHandler/NEUTInputHandler.cxx b/src/InputHandler/NEUTInputHandler.cxx
index e086bb0..6c43e00 100644
--- a/src/InputHandler/NEUTInputHandler.cxx
+++ b/src/InputHandler/NEUTInputHandler.cxx
@@ -1,516 +1,524 @@
#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]);
+ 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!");
+ "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!");
+ 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]);
+ << 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,
const bool lightweight) {
// 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);
+ << " 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);
+ << " 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.");
+ 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/Reweight/GENIEWeightEngine.cxx b/src/Reweight/GENIEWeightEngine.cxx
index 827a00f..dcc1aa6 100644
--- a/src/Reweight/GENIEWeightEngine.cxx
+++ b/src/Reweight/GENIEWeightEngine.cxx
@@ -1,478 +1,492 @@
#include "GENIEWeightEngine.h"
#ifdef __GENIE_ENABLED__
#ifdef GENIE_PRE_R3
#include "Algorithm/AlgConfigPool.h"
#include "EVGCore/EventRecord.h"
#include "GHEP/GHepRecord.h"
#include "Ntuple/NtpMCEventRecord.h"
#ifndef __NO_REWEIGHT__
#include "ReWeight/GReWeightAGKY.h"
#include "ReWeight/GReWeightDISNuclMod.h"
#include "ReWeight/GReWeightFGM.h"
#include "ReWeight/GReWeightFZone.h"
#include "ReWeight/GReWeightINuke.h"
#include "ReWeight/GReWeightNonResonanceBkg.h"
#include "ReWeight/GReWeightNuXSecCCQE.h"
#include "ReWeight/GReWeightNuXSecCCQEvec.h"
#include "ReWeight/GReWeightNuXSecCCRES.h"
#include "ReWeight/GReWeightNuXSecCOH.h"
#include "ReWeight/GReWeightNuXSecDIS.h"
#include "ReWeight/GReWeightNuXSecNC.h"
#include "ReWeight/GReWeightNuXSecNCEL.h"
#include "ReWeight/GReWeightNuXSecNCRES.h"
#include "ReWeight/GReWeightResonanceDecay.h"
#include "ReWeight/GSystUncertainty.h"
#ifdef __GENIE_EMP_MECRW_ENABLED
#include "ReWeight/GReWeightXSecEmpiricalMEC.h"
#endif
#endif
#if __GENIE_VERSION__ >= 212
#include "ReWeight/GReWeightNuXSecCCQEaxial.h"
#endif
#else // GENIE v3
#include "Framework/Algorithm/AlgConfigPool.h"
#include "Framework/EventGen/EventRecord.h"
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/GHEP/GHepRecord.h"
#include "Framework/Ntuple/NtpMCEventRecord.h"
#include "Framework/Utils/AppInit.h"
#include "Framework/Utils/RunOpt.h"
using namespace genie;
#ifndef __NO_REWEIGHT__
#include "RwCalculators/GReWeightAGKY.h"
#include "RwCalculators/GReWeightDISNuclMod.h"
#include "RwCalculators/GReWeightFGM.h"
#include "RwCalculators/GReWeightFZone.h"
#include "RwCalculators/GReWeightINuke.h"
#include "RwCalculators/GReWeightNonResonanceBkg.h"
#include "RwCalculators/GReWeightNuXSecCCQE.h"
#include "RwCalculators/GReWeightNuXSecCCQEaxial.h"
#include "RwCalculators/GReWeightNuXSecCCQEvec.h"
#include "RwCalculators/GReWeightNuXSecCCRES.h"
#include "RwCalculators/GReWeightNuXSecCOH.h"
#include "RwCalculators/GReWeightNuXSecDIS.h"
#include "RwCalculators/GReWeightNuXSecNC.h"
#include "RwCalculators/GReWeightNuXSecNCEL.h"
#include "RwCalculators/GReWeightNuXSecNCRES.h"
+#ifdef USE_GENIE_XSECMEC
+#include "RwCalculators/GReWeightXSecMEC.h"
+#endif
#include "RwCalculators/GReWeightResonanceDecay.h"
#include "RwCalculators/GReWeightXSecEmpiricalMEC.h"
#include "RwFramework/GSystUncertainty.h"
using namespace genie::rew;
#endif
#endif
#endif
#include "FitLogger.h"
GENIEWeightEngine::GENIEWeightEngine(std::string name) {
#ifdef __GENIE_ENABLED__
#ifndef __NO_REWEIGHT__
// Setup the NEUT Reweight engien
fCalcName = name;
NUIS_LOG(DEB, "Setting up GENIE RW : " << fCalcName);
// Create RW Engine suppressing cout
StopTalking();
fGenieRW = new genie::rew::GReWeight();
// Get List of Vetos (Just for debugging)
std::string rw_engine_list =
FitPar::Config().GetParS("FitWeight_fGenieRW_veto");
bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos;
bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos;
bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos;
bool xsec_nnres = rw_engine_list.find("xsec_nonresbkg") == std::string::npos;
bool xsec_nudis = rw_engine_list.find("nuclear_dis") == std::string::npos;
bool xsec_resdec =
rw_engine_list.find("hadro_res_decay") == std::string::npos;
bool xsec_fzone = rw_engine_list.find("hadro_intranuke") == std::string::npos;
bool xsec_intra = rw_engine_list.find("hadro_fzone") == std::string::npos;
bool xsec_agky = rw_engine_list.find("hadro_agky") == std::string::npos;
bool xsec_qevec = rw_engine_list.find("xsec_ccqe_vec") == std::string::npos;
bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos;
bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos;
bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos;
bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos;
bool xsec_nucqe = rw_engine_list.find("nuclear_qe") == std::string::npos;
bool xsec_qeaxial =
rw_engine_list.find("xsec_ccqe_axial") == std::string::npos;
#ifdef __GENIE_EMP_MECRW_ENABLED
bool xsec_empMEC = rw_engine_list.find("xsec_empMEC") == std::string::npos;
#endif
+#ifdef USE_GENIE_XSECMEC
+ bool xsec_MEC = rw_engine_list.find("xsec_MEC") == std::string::npos;
+#endif
#ifndef GENIE_PRE_R3
genie::RunOpt *grunopt = genie::RunOpt::Instance();
grunopt->EnableBareXSecPreCalc(true);
grunopt->SetEventGeneratorList(Config::GetParS("GENIEEventGeneratorList"));
- if(!Config::HasPar("GENIETune")){
- NUIS_ABORT("GENIE tune was not specified, this is required when reweighting GENIE V3+ events. Add a config parameter like: to your nuisance card.");
+ if (!Config::HasPar("GENIETune")) {
+ NUIS_ABORT(
+ "GENIE tune was not specified, this is required when reweighting GENIE "
+ "V3+ events. Add a config parameter like: to your nuisance card.");
}
grunopt->SetTuneName(Config::GetParS("GENIETune"));
grunopt->BuildTune();
std::string genv =
std::string(getenv("GENIE")) + "/config/Messenger_laconic.xml";
genie::utils::app_init::MesgThresholds(genv);
#endif
// Now actually add the RW Calcs
if (xsec_ncel)
fGenieRW->AdoptWghtCalc("xsec_ncel", new genie::rew::GReWeightNuXSecNCEL);
if (xsec_ccqe) {
fGenieRW->AdoptWghtCalc("xsec_ccqe", new genie::rew::GReWeightNuXSecCCQE);
// (dynamic_cast (fGenieRW->WghtCalc("xsec_ccqe")))
// ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCCQE") );
}
#ifdef __GENIE_EMP_MECRW_ENABLED
if (xsec_empMEC) {
fGenieRW->AdoptWghtCalc("xsec_empMEC",
new genie::rew::GReWeightXSecEmpiricalMEC);
}
#endif
+#ifdef USE_GENIE_XSECMEC
+ if (xsec_MEC) {
+ fGenieRW->AdoptWghtCalc("xsec_MEC", new genie::rew::GReWeightXSecMEC);
+ }
+#endif
if (xsec_coh) {
fGenieRW->AdoptWghtCalc("xsec_coh", new genie::rew::GReWeightNuXSecCOH());
// (dynamic_cast (fGenieRW->WghtCalc("xsec_coh")))
// ->SetXSecModel( FitPar::Config().GetParS("GENIEXSecModelCOH") );
}
if (xsec_nnres)
fGenieRW->AdoptWghtCalc("xsec_nonresbkg",
new genie::rew::GReWeightNonResonanceBkg);
if (xsec_nudis)
fGenieRW->AdoptWghtCalc("nuclear_dis", new genie::rew::GReWeightDISNuclMod);
if (xsec_resdec)
fGenieRW->AdoptWghtCalc("hadro_res_decay",
new genie::rew::GReWeightResonanceDecay);
if (xsec_fzone)
fGenieRW->AdoptWghtCalc("hadro_fzone", new genie::rew::GReWeightFZone);
if (xsec_intra)
fGenieRW->AdoptWghtCalc("hadro_intranuke", new genie::rew::GReWeightINuke);
if (xsec_agky)
fGenieRW->AdoptWghtCalc("hadro_agky", new genie::rew::GReWeightAGKY);
if (xsec_qevec)
fGenieRW->AdoptWghtCalc("xsec_ccqe_vec",
new genie::rew::GReWeightNuXSecCCQEvec);
#if __GENIE_VERSION__ >= 212
if (xsec_qeaxial)
fGenieRW->AdoptWghtCalc("xsec_ccqe_axial",
new genie::rew::GReWeightNuXSecCCQEaxial);
#endif
if (xsec_dis)
fGenieRW->AdoptWghtCalc("xsec_dis", new genie::rew::GReWeightNuXSecDIS);
if (xsec_nc)
fGenieRW->AdoptWghtCalc("xsec_nc", new genie::rew::GReWeightNuXSecNC);
if (xsec_ccres) {
#if __GENIE_VERSION__ < 213
fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES);
#else
fGenieRW->AdoptWghtCalc(
"xsec_ccres",
new genie::rew::GReWeightNuXSecCCRES(
FitPar::Config().GetParS("GENIEXSecModelCCRES"), "Default"));
#endif
}
if (xsec_ncres)
fGenieRW->AdoptWghtCalc("xsec_ncres", new genie::rew::GReWeightNuXSecNCRES);
if (xsec_nucqe)
fGenieRW->AdoptWghtCalc("nuclear_qe", new genie::rew::GReWeightFGM);
// Set the CCQE reweighting style
GReWeightNuXSecCCQE *rwccqe =
dynamic_cast(fGenieRW->WghtCalc("xsec_ccqe"));
// For MaCCQE reweighting
std::string ccqetype = FitPar::Config().GetParS("GENIEWeightEngine_CCQEMode");
if (ccqetype == "kModeMa") {
NUIS_LOG(DEB, "Setting GENIE ReWeight CCQE to kModeMa");
rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa);
} else if (ccqetype == "kModeNormAndMaShape") {
NUIS_LOG(DEB, "Setting GENIE ReWeight CCQE to kModeNormAndMaShape");
rwccqe->SetMode(GReWeightNuXSecCCQE::kModeNormAndMaShape);
// For z-expansion reweighting
} else if (ccqetype == "kModeZExp") {
NUIS_LOG(DEB, "Setting GENIE ReWeight CCQE to kModeZExp");
rwccqe->SetMode(GReWeightNuXSecCCQE::kModeZExp);
} else {
NUIS_ERR(FTL, "Did not find specified GENIE ReWeight CCQE mode");
NUIS_ABORT("You provided: " << ccqetype << " in parameters/config.xml");
}
#if (__GENIE_VERSION__ >= 212) and \
(__GENIE_VERSION__ <= \
300) // This doesn't currently work as is for GENIE v3, but the reweighting
// in v3 supposedly does similar checks anyway.
// Check the UserPhysicsOptions too!
AlgConfigPool *Pool = genie::AlgConfigPool::Instance();
Registry *full = Pool->GlobalParameterList();
std::string name_ax = full->GetAlg("AxialFormFactorModel").name;
std::string config_ax = full->GetAlg("AxialFormFactorModel").config;
if (name_ax == "genie::DipoleAxialFormFactorModel" &&
ccqetype == "kModeZExp") {
NUIS_ERR(
FTL,
"Trying to run Z Expansion reweighting with Llewelyn-Smith model.");
NUIS_ERR(FTL, "Please check your "
<< std::getenv("GENIE")
<< "/config/UserPhysicsOptions.xml to match generated");
NUIS_ERR(FTL, "You're telling me " << name_ax << "/" << config_ax);
NUIS_ABORT("Also check your "
<< std::getenv("NUISANCE")
<< "/parameters/config.xml GENIEWeightEngine_CCQEMode: "
<< ccqetype);
}
if (name_ax == "genie::ZExpAxialFormFactorModel" && ccqetype != "kModeZExp") {
NUIS_ERR(
FTL,
"Trying to run Llewelyn-Smith reweighting with Z Expansion model.");
NUIS_ERR(FTL, "Please change your "
<< std::getenv("GENIE")
<< "/config/UserPhysicsOptions.xml to match generated");
NUIS_ERR(FTL, "You're telling me " << name_ax << "/" << config_ax);
NUIS_ABORT("Also check your "
<< std::getenv("NUISANCE")
<< "/parameters/config.xml GENIEWeightEngine_CCQEMode: "
<< ccqetype);
}
std::string name_qelcc =
full->GetAlg("XSecModel@genie::EventGenerator/QEL-CC").name;
std::string config_qelcc =
full->GetAlg("XSecModel@genie::EventGenerator/QEL-CC").config;
if (config_qelcc == "Default" && ccqetype == "kModeZExp") {
NUIS_ERR(
FTL,
"Trying to run Z Expansion reweighting with Llewelyn-Smith model.");
NUIS_ERR(FTL, "Please change your "
<< std::getenv("GENIE")
<< "/config/UserPhysicsOptions.xml to match generated");
NUIS_ERR(FTL, "You're telling me " << name_qelcc << "/" << config_qelcc);
NUIS_ABORT("Also check your "
<< std::getenv("NUISANCE")
<< "/parameters/config.xml GENIEWeightEngine_CCQEMode: "
<< ccqetype);
}
if (config_qelcc == "ZExp" && ccqetype != "kModeZExp") {
NUIS_ERR(
FTL,
"Trying to run Llewelyn-Smith reweighting with Z Expansion model.");
NUIS_ERR(FTL, "Please change your "
<< std::getenv("GENIE")
<< "/config/UserPhysicsOptions.xml to match generated");
NUIS_ERR(FTL, "You're telling me " << name_qelcc << "/" << config_qelcc);
NUIS_ABORT("Also check your "
<< std::getenv("NUISANCE")
<< "/parameters/config.xml GENIEWeightEngine_CCQEMode: "
<< ccqetype);
}
#endif
if (xsec_ccres) {
// Default to include shape and normalization changes for CCRES (can be
// changed downstream if desired)
GReWeightNuXSecCCRES *rwccres =
dynamic_cast(fGenieRW->WghtCalc("xsec_ccres"));
std::string marestype =
FitPar::Config().GetParS("GENIEWeightEngine_CCRESMode");
if (!marestype.compare("kModeNormAndMaMvShape")) {
rwccres->SetMode(GReWeightNuXSecCCRES::kModeNormAndMaMvShape);
} else if (!marestype.compare("kModeMaMv")) {
rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv);
} else {
NUIS_ABORT("Unkown MARES Mode in GENIE Weight Engine : " << marestype);
}
}
if (xsec_ncres) {
// Default to include shape and normalization changes for NCRES (can be
// changed downstream if desired)
GReWeightNuXSecNCRES *rwncres =
dynamic_cast(fGenieRW->WghtCalc("xsec_ncres"));
rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv);
}
if (xsec_dis) {
// Default to include shape and normalization changes for DIS (can be
// changed downstream if desired)
GReWeightNuXSecDIS *rwdis =
dynamic_cast(fGenieRW->WghtCalc("xsec_dis"));
rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12u);
// Set Abs Twk Config
fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
}
// allow cout again
StartTalking();
#else
NUIS_ERR(FTL,
"GENIE ReWeight is __NOT ENABLED__ in GENIE and you're trying to "
"run NUISANCE with it enabled");
NUIS_ERR(FTL, "Check your genie-config --libs for reweighting");
NUIS_ERR(FTL, "If not present you need to recompile GENIE");
NUIS_ABORT("If present you need to contact NUISANCE authors");
#endif
#endif
};
void GENIEWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef __GENIE_ENABLED__
#ifndef __NO_REWEIGHT__
// Get First enum
int nuisenum = Reweight::ConvDial(name, kGENIE);
// Check ZExp sillyness in GENIE
// If ZExpansion parameters are used we need to set a different mode in GENIE
// ReWeight... GENIE doesn't have a setter either...
#if (__GENIE_VERSION__ >= 212) and (__GENIE_VERSION__ <= 300)
std::string ccqetype = FitPar::Config().GetParS("GENIEWeightEngine_CCQEMode");
if (ccqetype != "kModeZExp" &&
(name == "ZExpA1CCQE" || name == "ZExpA2CCQE" || name == "ZExpA3CCQE" ||
name == "ZExpA4CCQE")) {
NUIS_ERR(FTL, "Found a Z-expansion parameter in GENIE although the GENIE "
"ReWeighting engine is set to use Llewelyn-Smith and MaQE!");
NUIS_ABORT("Change your GENIE UserPhysicsOptions.xml in "
<< std::getenv("GENIE")
<< "/config/UserPhysicsOptions.xml to match requirements");
}
if ((ccqetype != "kModeMa" && ccqetype != "kModeMaNormAndMaShape") &&
(name == "MaCCQE")) {
NUIS_ERR(FTL, "Found MaCCQE parameter in GENIE although the GENIE "
"ReWeighting engine is set to not use this!");
NUIS_ABORT("Change your GENIE UserPhysicsOptions.xml in "
<< std::getenv("GENIE")
<< "/config/UserPhysicsOptions.xml to match requirements");
}
#endif
// Setup Maps
fEnumIndex[nuisenum]; // = std::vector(0);
fNameIndex[name]; // = std::vector(0);
// Split by commas
std::vector allnames = GeneralUtils::ParseToStr(name, ",");
for (uint i = 0; i < allnames.size(); i++) {
std::string singlename = allnames[i];
// Get RW
genie::rew::GSyst_t rwsyst = GSyst::FromString(singlename);
// Fill Maps
int index = fValues.size();
fValues.push_back(0.0);
fGENIESysts.push_back(rwsyst);
// Initialize dial
NUIS_LOG(DEB, "Registering " << singlename << " from " << name);
fGenieRW->Systematics().Init(fGENIESysts[index]);
// If Absolute
if (fIsAbsTwk) {
GSystUncertainty::Instance()->SetUncertainty(rwsyst, 1.0, 1.0);
}
// Setup index
fEnumIndex[nuisenum].push_back(index);
fNameIndex[name].push_back(index);
}
// Set Value if given
if (startval != -999.9) {
SetDialValue(nuisenum, startval);
}
#endif
#endif
}
void GENIEWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef __GENIE_ENABLED__
#ifndef __NO_REWEIGHT__
std::vector indices = fEnumIndex[nuisenum];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val);
}
#endif
#endif
}
void GENIEWeightEngine::SetDialValue(std::string name, double val) {
#ifdef __GENIE_ENABLED__
#ifndef __NO_REWEIGHT__
std::vector indices = fNameIndex[name];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val);
}
#endif
#endif
}
void GENIEWeightEngine::Reconfigure(bool silent) {
#ifdef __GENIE_ENABLED__
#ifndef __NO_REWEIGHT__
// Hush now...
if (silent)
StopTalking();
// Reconf
fGenieRW->Reconfigure();
fGenieRW->Print();
// Shout again
if (silent)
StartTalking();
#endif
#endif
}
double GENIEWeightEngine::CalcWeight(BaseFitEvt *evt) {
double rw_weight = 1.0;
#ifdef __GENIE_ENABLED__
#ifndef __NO_REWEIGHT__
// Make nom weight
if (!evt) {
NUIS_ABORT("evt not found : " << evt);
}
// Skip Non GENIE
if (evt->fType != kGENIE)
return 1.0;
if (!(evt->genie_event)) {
NUIS_ABORT("evt->genie_event not found!" << evt->genie_event);
}
if (!(evt->genie_event->event)) {
NUIS_ABORT("evt->genie_event->event GHepRecord not found!"
<< (evt->genie_event->event));
}
if (!fGenieRW) {
NUIS_ABORT("GENIE RW Not Found!" << fGenieRW);
}
rw_weight = fGenieRW->CalcWeight(*(evt->genie_event->event));
// std::cout << "Returning GENIE Weight for electron scattering = " <<
// rw_weight << std::endl;
// if (rw_weight != 1.0 )std::cout << "mode=" << evt->Mode << " rw_weight = "
// << rw_weight << std::endl;
#endif
#endif
// Return rw_weight
return rw_weight;
}