diff --git a/cmake/MINERvASetup.cmake b/cmake/MINERvASetup.cmake
new file mode 100644
index 0000000..0ccfcd9
--- /dev/null
+++ b/cmake/MINERvASetup.cmake
@@ -0,0 +1,19 @@
+# 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 .
+################################################################################
+LIST(APPEND EXTRA_CXX_FLAGS -D__MINERVA_RW_ENABLED__)
diff --git a/cmake/ReweightEnginesSetup.cmake b/cmake/ReweightEnginesSetup.cmake
index 44c80ff..87f47d4 100644
--- a/cmake/ReweightEnginesSetup.cmake
+++ b/cmake/ReweightEnginesSetup.cmake
@@ -1,82 +1,87 @@
# 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 .
################################################################################
################################## NEUT ######################################
if(USE_NEUT)
include(${CMAKE_SOURCE_DIR}/cmake/NEUTSetup.cmake)
cmessage(STATUS "Using NEUT Reweight engine.")
set(USE_NEUT TRUE CACHE BOOL "Whether to enable NEUT (reweight) support. Requires external libraries. " FORCE)
endif()
################################# NuWro ######################################
if(USE_NuWro)
include(${CMAKE_SOURCE_DIR}/cmake/NuWroSetup.cmake)
cmessage(STATUS "Using NuWro Reweight engine.")
set(USE_NuWro TRUE CACHE BOOL "Whether to enable NuWro support. " FORCE)
endif()
################################## GENIE #####################################
if(USE_GENIE)
include(${CMAKE_SOURCE_DIR}/cmake/GENIESetup.cmake)
cmessage(STATUS "Using GENIE Reweight engine.")
set(USE_GENIE TRUE CACHE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. " FORCE)
endif()
################################## NIWG ######################################
if(USE_NIWG)
include(${CMAKE_SOURCE_DIR}/cmake/NIWGSetup.cmake)
cmessage(STATUS "Using NIWG Reweight engine.")
set(USE_NIWG TRUE CACHE BOOL "Whether to enable (T2K) NIWG ReWeight support. Requires external libraries. " FORCE)
endif()
################################## T2K ######################################
if(USE_T2K)
include(${CMAKE_SOURCE_DIR}/cmake/T2KSetup.cmake)
cmessage(STATUS "Using T2K Reweight engine.")
set(USE_T2K TRUE CACHE BOOL "Whether to enable T2KReWeight support. Requires external libraries. " FORCE)
endif()
-
+################################## MINERvA ######################################
+if(USE_MINERvA_RW)
+ include(${CMAKE_SOURCE_DIR}/cmake/MINERvASetup.cmake)
+ cmessage(STATUS "Using MINERvA Reweight engine.")
+ set(USE_MINERvA_RW TRUE CACHE BOOL "Whether to enable MINERvA ReWeight support. " FORCE)
+endif()
################################################################################
################################ Prob3++ ####################################
include(${CMAKE_SOURCE_DIR}/cmake/Prob3++Setup.cmake)
################################################################################
cmessage(STATUS "Reweight engine include directories: ${RWENGINE_INCLUDE_DIRECTORIES}")
if(NEED_ROOTEVEGEN)
cmessage(STATUS "Require ROOT eve generation libraries")
LIST(REVERSE ROOT_LIBS)
LIST(APPEND ROOT_LIBS
Gui
Ged
Geom
TreePlayer
EG
Eve)
LIST(REVERSE ROOT_LIBS)
endif()
if(NEED_ROOTPYTHIA6)
cmessage(STATUS "Require ROOT Pythia6 libraries")
LIST(APPEND ROOT_LIBS
EGPythia6
Pythia6)
endif()
LIST(APPEND EXTRA_LIBS ${ROOT_LIBS})
diff --git a/cmake/cacheVariables.cmake b/cmake/cacheVariables.cmake
index 5e88e9c..1c75a96 100644
--- a/cmake/cacheVariables.cmake
+++ b/cmake/cacheVariables.cmake
@@ -1,206 +1,208 @@
# 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 .
################################################################################
function(CheckAndSetDefaultEnv VARNAME DEFAULT CACHETYPE DOCSTRING ENVNAME)
#cmessage(DEBUG "Trying to assign variable ${VARNAME} into the cache.")
if(NOT DEFINED ${VARNAME})
if(DEFINED ENV{${ENVNAME}} AND NOT $ENV{${ENVNAME}} STREQUAL "")
set(${VARNAME} $ENV{${ENVNAME}} CACHE ${CACHETYPE} ${DOCSTRING})
cmessage(DEBUG " Read ${VARNAME} from ENVVAR ${ENVNAME} as $ENV{${ENVNAME}}.")
else()
set(${VARNAME} ${DEFAULT} CACHE ${CACHETYPE} ${DOCSTRING})
endif()
else()
set(${VARNAME} ${${VARNAME}} CACHE ${CACHETYPE} ${DOCSTRING})
unset(${VARNAME})
endif()
cmessage(CACHE "--Set cache variable: \"${VARNAME}\" to \"${${VARNAME}}\", in cache ${CACHETYPE}.")
endfunction()
function(CheckAndSetDefaultCache VARNAME DEFAULT CACHETYPE DOCSTRING)
# cmessage(DEBUG "Trying to assign variable ${VARNAME} into the cache.")
if(NOT DEFINED ${VARNAME})
set(${VARNAME} ${DEFAULT} CACHE ${CACHETYPE} ${DOCSTRING})
else()
set(${VARNAME} ${${VARNAME}} CACHE ${CACHETYPE} ${DOCSTRING})
unset(${VARNAME})
endif()
cmessage(CACHE "--Set cache variable: \"${VARNAME}\" to \"${${VARNAME}}\", in cache ${CACHETYPE}.")
endfunction()
function(CheckAndSetDefault VARNAME DEFAULT)
# cmessage(DEBUG "Trying to assign variable ${VARNAME}.")
if(NOT DEFINED ${VARNAME})
set(${VARNAME} ${DEFAULT} PARENT_SCOPE)
set(${VARNAME} ${DEFAULT})
endif()
cmessage(CACHE "--Set variable: \"${VARNAME}\" to \"${${VARNAME}}\".")
endfunction()
CheckAndSetDefaultCache(VERBOSE TRUE BOOL "Whether to configure loudly.")
set (CMAKE_SKIP_BUILD_RPATH TRUE)
#Changes default install path to be a subdirectory of the build dir.
#Can set build dir at configure time with -DCMAKE_INSTALL_PREFIX=/install/path
if(CMAKE_INSTALL_PREFIX STREQUAL "" OR CMAKE_INSTALL_PREFIX STREQUAL
"/usr/local")
set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}")
elseif(NOT DEFINED CMAKE_INSTALL_PREFIX)
set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}")
endif()
if(CMAKE_BUILD_TYPE STREQUAL "")
set(CMAKE_BUILD_TYPE DEBUG)
elseif(NOT DEFINED CMAKE_BUILD_TYPE)
set(CMAKE_BUILD_TYPE DEBUG)
endif()
CheckAndSetDefaultCache(USE_MINIMIZER TRUE INTERNAL "Whether we are using the ROOT minimization libraries. ")
CheckAndSetDefaultCache(USE_HEPMC FALSE BOOL "Whether to enable HepMC input support. ")
CheckAndSetDefaultEnv(HEPMC "" PATH "Path to HepMC source tree root directory. Overrides environment variable \$HEPMC <>" HEPMC)
CheckAndSetDefaultCache(HEPMC_MOMUNIT "GEV" STRING "HepMC momentum units [MEV|GEV]. ")
CheckAndSetDefaultCache(HEPMC_LENUNIT "CM" STRING "HepMC momentum units [MM|CM]. ")
CheckAndSetDefaultCache(HEPMC_USED_EP FALSE INTERNAL "Whether we built HepMC or not. ")
CheckAndSetDefaultCache(USE_NEUT FALSE BOOL "Whether to enable NEUT (reweight) support. Requires external libraries. ")
CheckAndSetDefaultEnv(NEUT_ROOT "" PATH "Path to NEUT source tree root directory. Overrides environment variable \$NEUT_ROOT <>" NEUT_ROOT)
CheckAndSetDefaultEnv(CERN "" PATH "Path to CERNLIB source tree root directory that NEUT was built against. Overrides environment variable \$CERN <>" CERN)
CheckAndSetDefaultEnv(CERN_LEVEL "" STRING "CERNLIB Library version. Overrides environment variable \$CERN_LEVEL <>" CERN_LEVEL)
CheckAndSetDefaultCache(USE_NuWro FALSE BOOL "Whether to enable NuWro support. ")
CheckAndSetDefaultEnv(NUWRO "" PATH "Path to NuWro source tree root directory. Overrides environment variable \$NUWRO <>" NUWRO)
CheckAndSetDefaultEnv(NUWRO_INC "" PATH "Path to NuWro installed includes directory, needs to contain \"params_all.h\". Overrides environment variable \$NUWRO_INC <>" NUWRO_INC)
CheckAndSetDefaultCache(NUWRO_INPUT_FILE "" FILEPATH "Path to an input NuWro event vector, which can be used to build NuWro i/o libraries. <>")
CheckAndSetDefaultCache(NUWRO_BUILT_FROM_FILE FALSE INTERNAL "Whether the NuWro libraries were built by NUISANCE. ")
CheckAndSetDefaultCache(USE_NuWro_RW FALSE BOOL "Whether to try and build support for NuWro reweighting. ")
CheckAndSetDefaultCache(USE_NuWro_SRW_Event FALSE BOOL "Whether to use cut down NuWro reweight event format. Requires NuWro reweight. ")
CheckAndSetDefaultCache(USE_GENIE FALSE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. ")
CheckAndSetDefaultEnv(GENIE "" PATH "Path to GENIE source tree root directory. Overrides environment variable \$GENIE <>" GENIE)
CheckAndSetDefaultEnv(LHAPDF_LIB "" PATH "Path to pre-built LHAPDF libraries. Overrides environment variable \$LHAPDF_LIB. <>" LHAPDF_LIB)
CheckAndSetDefaultEnv(LHAPDF_INC "" PATH "Path to installed LHAPDF headers. Overrides environment variable \$LHAPDF_INC. <>" LHAPDF_INC)
CheckAndSetDefaultEnv(LHAPATH "" PATH "Path to LHA PDF inputs. Overrides environment variable \$LHAPATH. <>" LHAPATH)
CheckAndSetDefaultEnv(LIBXML2_LIB "" PATH "Path to pre-built LIBXML2 libraries. Overrides environment variable \$LIBXML2_LIB. <>" LIBXML2_LIB)
CheckAndSetDefaultEnv(LIBXML2_INC "" PATH "Path to installed LIBXML2 headers. Overrides environment variable \$LIBXML2_INC. <>" LIBXML2_INC)
CheckAndSetDefaultEnv(LOG4CPP_LIB "" PATH "Path to pre-built LOG4CPP libraries. Overrides environment variable \$LOG4CPP_LIB. <>" LOG4CPP_LIB)
CheckAndSetDefaultEnv(LOG4CPP_INC "" PATH "Path to installed LOG4CPP headers. Overrides environment variable \$LOG4CPP_INC. <>" LOG4CPP_INC)
CheckAndSetDefaultCache(BUILD_GEVGEN FALSE BOOL "Whether to build nuisance_gevgen app.")
CheckAndSetDefaultCache(USE_T2K FALSE BOOL "Whether to enable T2KReWeight support. Requires external libraries. ")
CheckAndSetDefaultEnv(T2KREWEIGHT "" PATH "Path to installed T2KREWEIGHTReWeight. Overrides environment variable \$T2KREWEIGHT. <>" T2KREWEIGHT)
CheckAndSetDefaultCache(USE_NIWG FALSE BOOL "Whether to enable (T2K) NIWG ReWeight support. Requires external libraries. ")
CheckAndSetDefaultEnv(NIWG_ROOT "" PATH "Path to installed NIWGReWeight. Overrides environment variable \$NIWG. <>" NIWG)
+CheckAndSetDefaultCache(USE_MINERvA_RW FALSE BOOL "Whether to enable MINERvA ReWeight support. ")
+
CheckAndSetDefaultEnv(PYTHIA6 "" PATH "Path to directory containing libPythia6.so. Overrides environment variable \$PYTHIA6 <>" PYTHIA6)
CheckAndSetDefaultEnv(PYTHIA8 "" PATH "Path to directory containing libPythia8.so. Overrides environment variable \$PYTHIA8 <>" PYTHIA8)
CheckAndSetDefaultCache(USE_PYTHIA8 FALSE BOOL "Whether to enable PYTHIA8 event support. ")
CheckAndSetDefaultCache(USE_GiBUU TRUE BOOL "Whether to enable GiBUU event support. ")
CheckAndSetDefaultCache(BUILD_GiBUU FALSE BOOL "Whether to build supporting GiBUU event tools along with a patched version of GiBUU. ")
CheckAndSetDefaultCache(USE_NUANCE TRUE BOOL "Whether to enable NUANCE event support. ")
CheckAndSetDefaultCache(USE_PROB3PP FALSE BOOL "Whether to download and compile in Prob3++ support. ")
CheckAndSetDefaultCache(NO_EXTERNAL_UPDATE TRUE BOOL "Whether to perform the update target for external dependencies. ")
CheckAndSetDefaultCache(USE_GPERFTOOLS FALSE BOOL "Whether to compile in google performance tools. ")
CheckAndSetDefault(NEED_PYTHIA6 FALSE)
CheckAndSetDefault(NEED_PYTHIA8 FALSE)
CheckAndSetDefault(NEED_ROOTEVEGEN FALSE)
CheckAndSetDefault(NEED_ROOTPYTHIA6 FALSE)
CheckAndSetDefaultCache(USE_OMP FALSE BOOL "Whether to enable multicore features (there currently are none...). ")
CheckAndSetDefault(NO_EXPERIMENTS FALSE)
cmessage(STATUS "NO_EXPERIMENTS: ${NO_EXPERIMENTS}")
CheckAndSetDefaultCache(NO_ANL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build ANL samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_ArgoNeuT ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build ArgoNeuT samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_BEBC ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build BEBC samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_BNL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build BNL samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_FNAL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build FNAL samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_GGM ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build GGM samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_K2K ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build K2K samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_MINERvA ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build MINERvA samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_MiniBooNE ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build MiniBooNE samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_T2K ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build T2K samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_SciBooNE ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build SciBooNE samples. <-DNO_EXPERIMENTS=FALSE>")
function(SAYVARS)
LIST(APPEND VARS
USE_HEPMC
HEPMC
HEPMC_MOMUNIT
HEPMC_LENUNIT
HEPMC_USED_EP
USE_NEUT
NEUT_ROOT
CERN
CERN_LEVEL
USE_NuWro
NUWRO
NUWRO_INC
NUWRO_INPUT_FILE
NUWRO_BUILT_FROM_FILE
USE_GENIE
GENIE
LHAPDF_LIB
LHAPDF_INC
LIBXML2_LIB
LIBXML2_INC
LOG4CPP_LIB
GENIE_LOG4CPP_INC
BUILD_GEVGEN
USE_T2K
USE_NIWG
USE_GiBUU
BUILD_GiBUU
USE_NUANCE
NO_EXTERNAL_UPDATE
USE_GPERFTOOLS
NO_ANL
NO_ArgoNeuT
NO_BEBC
NO_BNL
NO_FNAL
NO_GGM
NO_K2K
NO_MINERvA
NO_MiniBooNE
NO_T2K
NO_SciBooNE)
foreach(v ${VARS})
if(DEFINED ${v})
cmessage(DEBUG "VARIABLE: \"${v}\" = \"${${v}}\"")
endif()
endforeach(v)
endfunction()
diff --git a/src/Reweight/CMakeLists.txt b/src/Reweight/CMakeLists.txt
index 1f9270c..685133c 100644
--- a/src/Reweight/CMakeLists.txt
+++ b/src/Reweight/CMakeLists.txt
@@ -1,79 +1,81 @@
# 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 .
################################################################################
set(IMPLFILES
GlobalDialList.cxx
FitWeight.cxx
WeightEngineBase.cxx
NEUTWeightEngine.cxx
NuWroWeightEngine.cxx
GENIEWeightEngine.cxx
WeightUtils.cxx
SampleNormEngine.cxx
LikelihoodWeightEngine.cxx
SplineWeightEngine.cxx
NUISANCESyst.cxx
T2KWeightEngine.cxx
NUISANCEWeightEngine.cxx
NUISANCEWeightCalcs.cxx
NIWGWeightEngine.cxx
OscWeightEngine.cxx
+MINERvAWeightCalcs.cxx
)
set(HEADERFILES
GlobalDialList.h
FitWeight.h
WeightEngineBase.h
NEUTWeightEngine.h
NuWroWeightEngine.h
GENIEWeightEngine.h
WeightUtils.h
SampleNormEngine.h
LikelihoodWeightEngine.h
SplineWeightEngine.h
NUISANCESyst.h
T2KWeightEngine.h
NUISANCEWeightEngine.h
NUISANCEWeightCalcs.h
NIWGWeightEngine.h
OscWeightEngine.h
+MINERvAWeightCalcs.h
)
set(LIBNAME Reweight)
if(CMAKE_BUILD_TYPE MATCHES DEBUG)
add_library(${LIBNAME} STATIC ${IMPLFILES})
else(CMAKE_BUILD_TYPE MATCHES RELEASE)
add_library(${LIBNAME} SHARED ${IMPLFILES})
endif()
include_directories(${MINIMUM_INCLUDE_DIRECTORIES})
set_target_properties(${LIBNAME} PROPERTIES VERSION
"${ExtFit_VERSION_MAJOR}.${ExtFit_VERSION_MINOR}.${ExtFit_VERSION_REVISION}")
#set_target_properties(${LIBNAME} PROPERTIES LINK_FLAGS ${ROOT_LD_FLAGS})
if(DEFINED PROJECTWIDE_EXTRA_DEPENDENCIES)
add_dependencies(${LIBNAME} ${PROJECTWIDE_EXTRA_DEPENDENCIES})
endif()
install(TARGETS ${LIBNAME} DESTINATION lib)
#Can uncomment this to install the headers... but is it really neccessary?
install(FILES ${HEADERFILES} DESTINATION include)
set(MODULETargets ${MODULETargets} ${LIBNAME} PARENT_SCOPE)
diff --git a/src/Reweight/MINERvAWeightCalcs.cxx b/src/Reweight/MINERvAWeightCalcs.cxx
index 5fefde9..ce70e83 100644
--- a/src/Reweight/MINERvAWeightCalcs.cxx
+++ b/src/Reweight/MINERvAWeightCalcs.cxx
@@ -1,345 +1,403 @@
#ifdef __MINERVA_RW_ENABLED__
#ifdef __GENIE_ENABLED__
#include "MINERvAWeightCalcs.h"
#include "BaseFitEvt.h"
namespace nuisance {
namespace reweight {
+//*******************************************************
MINERvAReWeight_MEC::MINERvAReWeight_MEC() {
+ //*******************************************************
fTwk_NormCCMEC = 0.0;
fDef_NormCCMEC = 1.0;
fCur_NormCCMEC = fDef_NormCCMEC;
}
-virtual ~MINERvAReWeight_MEC::MINERvAReWeight_MEC(){};
+MINERvAReWeight_MEC::~MINERvAReWeight_MEC(){};
double MINERvAReWeight_MEC::CalcWeight(BaseFitEvt* evt) {
// Check GENIE
if (evt->fType != kGENIE) return 1.0;
// Extract the GENIE Record
GHepRecord* ghep = static_cast(evt->genie_event->event);
const Interaction* interaction = ghep->Summary();
const InitialState& init_state = interaction->InitState();
const ProcessInfo& proc_info = interaction->ProcInfo();
const Target& tgt = init_state.Tgt();
// If the event is not MEC this Calc doesn't handle it
if (!proc_info.IsMEC()) return 1.0;
// WEIGHT CALCULATIONS -------------
double w = 1.0;
// CCMEC Dial
- if (!proc_info.IsCC()) w *= fCur_NormCCMEC;
+ if (!proc_info.IsWeakCC()) w *= fCur_NormCCMEC;
// Return Combined Weight
return w;
}
void MINERvAReWeight_MEC::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void MINERvAReWeight_MEC::SetDialValue(int rwenum, double val) {
// Check Handled
int curenum = rwenum % 1000;
if (!IsHandled(curenum)) return;
// Set Values
- if (curenum == kMINERvARW_NormCCMEC) fTwk_NormCCMEC = val;
+ if (curenum == kMINERvARW_NormCCMEC) {
+ fTwk_NormCCMEC = val;
+ fCur_NormCCMEC = fDef_NormCCMEC + fTwk_NormCCMEC;
+ }
// Define Tweaked
fTweaked = ((fTwk_NormCCMEC != 0.0));
}
bool MINERvAReWeight_MEC::IsHandled(int rwenum) {
int curenum = rwenum % 1000;
switch (curenum) {
case kMINERvARW_NormCCMEC:
return true;
default:
return false;
}
}
+
+//*******************************************************
+MINERvAReWeight_RES::MINERvAReWeight_RES() {
+ //*******************************************************
+ fTwk_NormCCRES = 0.0;
+ fDef_NormCCRES = 1.0;
+ fCur_NormCCRES = fDef_NormCCRES;
+}
+
+MINERvAReWeight_RES::~MINERvAReWeight_RES(){};
+
+double MINERvAReWeight_RES::CalcWeight(BaseFitEvt* evt) {
+
+ // std::cout << "Caculating RES" << std::endl;
+ // Check GENIE
+ if (evt->fType != kGENIE) return 1.0;
+
+ // Extract the GENIE Record
+ GHepRecord* ghep = static_cast(evt->genie_event->event);
+ const Interaction* interaction = ghep->Summary();
+ const InitialState& init_state = interaction->InitState();
+ const ProcessInfo& proc_info = interaction->ProcInfo();
+ const Target& tgt = init_state.Tgt();
+
+ // If the event is not RES this Calc doesn't handle it
+ if (!proc_info.IsResonant()) return 1.0;
+
+ // WEIGHT CALCULATIONS -------------
+ double w = 1.0;
+
+ // CCRES Dial
+ if (proc_info.IsWeakCC()) w *= fCur_NormCCRES;
+
+ // Return Combined Weight
+ return w;
+}
+
+void MINERvAReWeight_RES::SetDialValue(std::string name, double val) {
+ SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
+}
+
+void MINERvAReWeight_RES::SetDialValue(int rwenum, double val) {
+ // Check Handled
+ int curenum = rwenum % 1000;
+ if (!IsHandled(curenum)) return;
+
+ // Set Values
+ if (curenum == kMINERvARW_NormCCRES) {
+ fTwk_NormCCRES = val;
+ fCur_NormCCRES = fDef_NormCCRES + fTwk_NormCCRES;
+ }
+
+ // Define Tweaked
+ fTweaked = ((fTwk_NormCCRES != 0.0));
+}
+
+bool MINERvAReWeight_RES::IsHandled(int rwenum) {
+ int curenum = rwenum % 1000;
+
+ switch (curenum) {
+ case kMINERvARW_NormCCRES:
+ return true;
+ default:
+ return false;
+ }
+}
+
+//*******************************************************
+RikRPA::RikRPA() {
+ //*******************************************************
+
+ // - Syst : kMINERvA_RikRPA_ApplyRPA
+ // - Type : Binary
+ // - Limits : 0.0 (false) -> 1.0 (true)
+ // - Default : 0.0
+ fApplyDial_RPACorrection = false;
+
+ // - Syst : kMINERvA_RikRPA_LowQ2
+ // - Type : Absolute
+ // - Limits : 1.0 -> 1.0
+ // - Default : 0.0
+ // - Frac Error : 100%
+ fDefDial_RPALowQ2 = 0.0;
+ fCurDial_RPALowQ2 = fDefDial_RPALowQ2;
+ fErrDial_RPALowQ2 = 0.0;
+
+ // - Syst : kMINERvA_RikRPA_HighQ2
+ // - Type : Absolute
+ // - Limits : 1.0 -> 1.0
+ // - Default : 0.0
+ // - Frac Error : 100%
+ fDefDial_RPAHighQ2 = 0.0;
+ fCurDial_RPAHighQ2 = fDefDial_RPAHighQ2;
+ fErrDial_RPAHighQ2 = 1.0;
+
+ fEventWeights = new double[5];
+
+ for (size_t i = 0; i < kMaxCalculators; i++) {
+ // fRPACalculators[i] = NULL;
+ }
+
+ fTweaked = false;
+}
+
+RikRPA::~RikRPA() {
+ delete fEventWeights;
+
+ for (size_t i = 0; i < kMaxCalculators; i++) {
+ // if (fRPACalculators[i]) delete fRPACalculators[i];
+ // fRPACalculators[i] = NULL;
+ }
+}
+
+double RikRPA::CalcWeight(BaseFitEvt* evt) {
+ // LOG(FIT) << "Calculating RikRPA" << std::endl;
+ // Return 1.0 if not tweaked
+ if (!fTweaked) return 1.0;
+
+ double w = 1.0;
+
+ // Extract the GENIE Record
+ GHepRecord* ghep = static_cast(evt->genie_event->event);
+ const Interaction* interaction = ghep->Summary();
+ const InitialState& init_state = interaction->InitState();
+ const ProcessInfo& proc_info = interaction->ProcInfo();
+ // const Kinematics & kine = interaction->Kine();
+ // const XclsTag & xcls = interaction->ExclTag();
+ const Target& tgt = init_state.Tgt();
+
+ // If not QE return 1.0
+ // LOG(FIT) << "RikRPA : Event QE = " << proc_info.IsQuasiElastic() <<
+ // std::endl;
+ if (!tgt.IsNucleus()) return 1.0;
+ if (!proc_info.IsQuasiElastic()) return 1.0;
+
+ // Extract Beam and Target PDG
+ GHepParticle* neutrino = ghep->Probe();
+ int bpdg = neutrino->Pdg();
+
+ GHepParticle* target = ghep->Particle(1);
+ assert(target);
+ int tpdg = target->Pdg();
+
+ // Find the enum we need
+ int calcenum = GetRPACalcEnum(bpdg, tpdg);
+ if (calcenum == -1) return 1.0;
+
+ // Check we have the RPA Calc setup for this enum
+ // if not, set it up at that point
+ // if (!fRPACalculators[calcenum]) SetupRPACalculator(calcenum);
+ // weightRPA* rpacalc = fRPACalculators[calcenum];
+ // if (!rpacalc) {
+ // THROW("Failed to grab the RPA Calculator : " << calcenum);
+ // }
+
+ // Extract Q0-Q3
+ GHepParticle* fsl = ghep->FinalStatePrimaryLepton();
+ const TLorentzVector& k1 = *(neutrino->P4());
+ const TLorentzVector& k2 = *(fsl->P4());
+ double q0 = fabs((k1 - k2).E());
+ double q3 = fabs((k1 - k2).Vect().Mag());
+
+ // Now use q0-q3 and RPA Calculator to fill fWeights
+ // LOG(FIT) << "Getting Weights = " << q0 << " " << q3 << std::endl;
+ // rpacalc->getWeight(q0, q3, fEventWeights);
+
+ // Apply Interpolation (for the time being simple linear)
+
+ // Syst Application : kMINERvA_RikRPA_ApplyRPA
+ if (fApplyDial_RPACorrection) {
+ w *= fEventWeights[0]; // CV
+ }
+
+ LOG(FIT) << " fCurDial_RPALowQ2 = " << fCurDial_RPALowQ2
+ << " fCurDial_RPAHighQ2 = " << fCurDial_RPAHighQ2 << " Weights "
+ << fEventWeights[0] << " " << fEventWeights[1] << " "
+ << fEventWeights[2] << " " << fEventWeights[3] << " "
+ << fEventWeights[4] << std::endl;
+
+ // Syst Application : kMINERvA_RikRPA_LowQ2
+ if (fabs(fCurDial_RPALowQ2) > 0.0) {
+ double interpw = fEventWeights[0];
+
+ if (fCurDial_RPALowQ2 > 0.0) {
+ interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[1]) *
+ fCurDial_RPALowQ2; // WLow+ } else if
+ } else if (fCurDial_RPALowQ2 < 0.0) {
+ interpw = fEventWeights[0] - (fEventWeights[2] - fEventWeights[0]) *
+ fCurDial_RPALowQ2; // WLow-
+ }
+ w *= interpw / fEventWeights[0]; // Div by CV again
+ }
+
+ // Syst Application : kMINERvA_RikRPA_HighQ2
+ if (fabs(fCurDial_RPAHighQ2) > 0.0) {
+ double interpw = fEventWeights[0];
+ if (fCurDial_RPAHighQ2 > 0.0) {
+ interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) *
+ fCurDial_RPAHighQ2; // WHigh+ } else
+ if (fCurDial_RPAHighQ2 < 0.0) {
+ interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) *
+ fCurDial_RPAHighQ2; // WHigh-
+ }
+ w *= interpw / fEventWeights[0]; // Div by CV again
+ }
+ }
+
+ // LOG(FIT) << "RPA Weight = " << w << std::endl;
+ return w;
} // namespace reweight
-} // namespace nuisance
-#endif
-#endif
+void RikRPA::SetDialValue(std::string name, double val) {
+ SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
+}
+
+void RikRPA::SetDialValue(int rwenum, double val) {
+ int curenum = rwenum % 1000;
+ // Check Handled
+ if (!IsHandled(curenum)) return;
+ if (curenum == kMINERvARW_RikRPA_ApplyRPA)
+ fApplyDial_RPACorrection = (val > 0.5);
+ if (curenum == kMINERvARW_RikRPA_LowQ2) fCurDial_RPALowQ2 = val;
+ if (curenum == kMINERvARW_RikRPA_HighQ2) fCurDial_RPAHighQ2 = val;
+
+ // Assign flag to say stuff has changed
+ fTweaked = (fApplyDial_RPACorrection ||
+ fabs(fCurDial_RPAHighQ2 - fDefDial_RPAHighQ2) > 0.0 ||
+ fabs(fCurDial_RPALowQ2 - fDefDial_RPALowQ2) > 0.0);
+}
+bool RikRPA::IsHandled(int rwenum) {
+ int curenum = rwenum % 1000;
+ switch (curenum) {
+ case kMINERvARW_RikRPA_ApplyRPA:
+ return true;
+ case kMINERvARW_RikRPA_LowQ2:
+ return true;
+ case kMINERvARW_RikRPA_HighQ2:
+ return true;
+ default:
+ return false;
+ }
+}
+void RikRPA::SetupRPACalculator(int calcenum) {
+ std::string rwdir = FitPar::GetDataBase() + "/reweight/MINERvA/RikRPA/";
+ std::string fidir = "";
+ switch (calcenum) {
+ case kNuMuC12:
+ fidir = "outNievesRPAratio-nu12C-20GeV-20170202.root";
+ break;
+
+ case kNuMuO16:
+ fidir = "outNievesRPAratio-nu16O-20GeV-20170202.root";
+ break;
+
+ case kNuMuAr40:
+ fidir = "outNievesRPAratio-nu40Ar-20GeV-20170202.root";
+ break;
+
+ case kNuMuCa40:
+ fidir = "outNievesRPAratio-nu40Ca-20GeV-20170202.root";
+ break;
+
+ case kNuMuFe56:
+ fidir = "outNievesRPAratio-nu56Fe-20GeV-20170202.root";
+ break;
+
+ case kNuMuBarC12:
+ fidir = "outNievesRPAratio-anu12C-20GeV-20170202.root";
+ break;
+
+ case kNuMuBarO16:
+ fidir = "outNievesRPAratio-anu16O-20GeV-20170202.root";
+ break;
+
+ case kNuMuBarAr40:
+ fidir = "outNievesRPAratio-anu40Ar-20GeV-20170202.root";
+ break;
+
+ case kNuMuBarCa40:
+ fidir = "outNievesRPAratio-anu40Ca-20GeV-20170202.root";
+ break;
+
+ case kNuMuBarFe56:
+ fidir = "outNievesRPAratio-anu56Fe-20GeV-20170202.root";
+ break;
+ }
+ LOG(FIT) << "Loading RPA CALC : " << fidir << std::endl;
+ TDirectory* olddir = gDirectory;
-// RikRPA::RikRPA() {
+ // fRPACalculators[calcenum] = new weightRPA(rwdir + "/" + fidir);
+ olddir->cd();
+ return;
+}
-// // - Syst : kMINERvA_RikRPA_ApplyRPA
-// // - Type : Binary
-// // - Limits : 0.0 (false) -> 1.0 (true)
-// // - Default : 0.0
-// fApplyDial_RPACorrection = false;
+int RikRPA::GetRPACalcEnum(int bpdg, int tpdg) {
+ if (bpdg == 14 && tpdg == 1000060120)
+ return kNuMuC12;
+ else if (bpdg == 14 && tpdg == 1000080160)
+ return kNuMuO16;
+ else if (bpdg == 14 && tpdg == 1000180400)
+ return kNuMuAr40;
+ else if (bpdg == 14 && tpdg == 1000200400)
+ return kNuMuCa40;
+ else if (bpdg == 14 && tpdg == 1000280560)
+ return kNuMuFe56;
+ else if (bpdg == -14 && tpdg == 1000060120)
+ return kNuMuBarC12;
+ else if (bpdg == -14 && tpdg == 1000080160)
+ return kNuMuBarO16;
+ else if (bpdg == -14 && tpdg == 1000180400)
+ return kNuMuBarAr40;
+ else if (bpdg == -14 && tpdg == 1000200400)
+ return kNuMuBarCa40;
+ else if (bpdg == -14 && tpdg == 1000280560)
+ return kNuMuBarFe56;
+ else {
+ ERROR(WRN, "Unknown beam and target combination for RPA Calcs! "
+ << bpdg << " " << tpdg);
+ }
-// // - Syst : kMINERvA_RikRPA_LowQ2
-// // - Type : Absolute
-// // - Limits : 1.0 -> 1.0
-// // - Default : 0.0
-// // - Frac Error : 100%
-// fDefDial_RPALowQ2 = 0.0;
-// fCurDial_RPALowQ2 = fDefDial_RPALowQ2;
-// fErrDial_RPALowQ2 = 0.0;
+ return -1;
+}
-// // - Syst : kMINERvA_RikRPA_HighQ2
-// // - Type : Absolute
-// // - Limits : 1.0 -> 1.0
-// // - Default : 0.0
-// // - Frac Error : 100%
-// fDefDial_RPAHighQ2 = 0.0;
-// fCurDial_RPAHighQ2 = fDefDial_RPAHighQ2;
-// fErrDial_RPAHighQ2 = 1.0;
-
-// fEventWeights = new double[5];
-
-// for (size_t i = 0; i < kMaxCalculators; i++) {
-// fRPACalculators[i] = NULL;
-// }
-
-// fTweaked = false;
-// }
-
-// RikRPA::~RikRPA() {
-
-// delete fEventWeights;
-
-// for (size_t i = 0; i < kMaxCalculators; i++) {
-// if (fRPACalculators[i]) delete fRPACalculators[i];
-// fRPACalculators[i] = NULL;
-// }
-
-// }
-
-// double RikRPA::CalcWeight(BaseFitEvt* evt) {
-// // LOG(FIT) << "Calculating RikRPA" << std::endl;
-// // Return 1.0 if not tweaked
-// if (!fTweaked) return 1.0;
-
-// double w = 1.0;
-
-// // Extract the GENIE Record
-// GHepRecord* ghep = static_cast(evt->genie_event->event);
-// const Interaction* interaction = ghep->Summary();
-// const InitialState& init_state = interaction->InitState();
-// const ProcessInfo& proc_info = interaction->ProcInfo();
-// // const Kinematics & kine = interaction->Kine();
-// // const XclsTag & xcls = interaction->ExclTag();
-// const Target& tgt = init_state.Tgt();
-
-// // If not QE return 1.0
-// // LOG(FIT) << "RikRPA : Event QE = " << proc_info.IsQuasiElastic() <<
-// std::endl;
-// if (!tgt.IsNucleus()) return 1.0;
-// if (!proc_info.IsQuasiElastic()) return 1.0;
-
-// // Extract Beam and Target PDG
-// GHepParticle* neutrino = ghep->Probe();
-// int bpdg = neutrino->Pdg();
-
-// GHepParticle* target = ghep->Particle(1);
-// assert(target);
-// int tpdg = target->Pdg();
-
-// // Find the enum we need
-// int calcenum = GetRPACalcEnum(bpdg, tpdg);
-// if (calcenum == -1) return 1.0;
-
-// // Check we have the RPA Calc setup for this enum
-// // if not, set it up at that point
-// if (!fRPACalculators[calcenum]) SetupRPACalculator(calcenum);
-// weightRPA* rpacalc = fRPACalculators[calcenum];
-// if (!rpacalc) {
-// THROW("Failed to grab the RPA Calculator : " << calcenum);
-// }
-
-// // Extract Q0-Q3
-// GHepParticle* fsl = ghep->FinalStatePrimaryLepton();
-// const TLorentzVector& k1 = *(neutrino->P4());
-// const TLorentzVector& k2 = *(fsl->P4());
-// double q0 = fabs((k1 - k2).E());
-// double q3 = fabs((k1 - k2).Vect().Mag());
-
-// // Now use q0-q3 and RPA Calculator to fill fWeights
-// // LOG(FIT) << "Getting Weights = " << q0 << " " << q3 << std::endl;
-// rpacalc->getWeight(q0, q3, fEventWeights);
-
-// // Apply Interpolation (for the time being simple linear)
-
-// // Syst Application : kMINERvA_RikRPA_ApplyRPA
-// if (fApplyDial_RPACorrection) {
-// w *= fEventWeights[0]; // CV
-// }
-
-// LOG(FIT) << " fCurDial_RPALowQ2 = " << fCurDial_RPALowQ2
-// << " fCurDial_RPAHighQ2 = " << fCurDial_RPAHighQ2 << " Weights "
-// << fEventWeights[0] << " " << fEventWeights[1] << " "
-// << fEventWeights[2] << " " << fEventWeights[3] << " "
-// << fEventWeights[4] << std::endl;
-
-// // Syst Application : kMINERvA_RikRPA_LowQ2
-// if (fabs(fCurDial_RPALowQ2) > 0.0) {
-// double interpw = fEventWeights[0];
-
-// if (fCurDial_RPALowQ2 > 0.0) {
-// interpw = fEventWeights[0] -
-// (fEventWeights[0] - fEventWeights[1]) *
-// fCurDial_RPALowQ2; // WLow+ } else if
-// (fCurDial_RPALowQ2 < 0.0) {
-// interpw = fEventWeights[0] - (fEventWeights[2] - fEventWeights[0]) *
-// fCurDial_RPALowQ2; // WLow-
-// }
-// w *= interpw / fEventWeights[0]; // Div by CV again
-// }
-
-// // Syst Application : kMINERvA_RikRPA_HighQ2
-// if (fabs(fCurDial_RPAHighQ2) > 0.0) {
-// double interpw = fEventWeights[0];
-// if (fCurDial_RPAHighQ2 > 0.0) {
-// interpw = fEventWeights[0] -
-// (fEventWeights[0] - fEventWeights[3]) *
-// fCurDial_RPAHighQ2; // WHigh+ } else
-// if
-// (fCurDial_RPAHighQ2 < 0.0) {
-// interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0])
-// *
-// fCurDial_RPAHighQ2; // WHigh-
-// }
-// w *= interpw / fEventWeights[0]; // Div by CV again
-// }
-
-// // LOG(FIT) << "RPA Weight = " << w << std::endl;
-
-// return w;
-// }
-
-// void RikRPA::SetDialValue(std::string name, double val) {
-// SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
-// }
-
-// void RikRPA::SetDialValue(int rwenum, double val) {
-// int curenum = rwenum % 1000;
-
-// // Check Handled
-// if (!IsHandled(curenum)) return;
-// if (curenum == kMINERvA_RikRPA_ApplyRPA)
-// fApplyDial_RPACorrection = (val > 0.5);
-// if (curenum == kMINERvA_RikRPA_LowQ2) fCurDial_RPALowQ2 = val;
-// if (curenum == kMINERvA_RikRPA_HighQ2) fCurDial_RPAHighQ2 = val;
-
-// // Assign flag to say stuff has changed
-// fTweaked = (fApplyDial_RPACorrection ||
-// fabs(fCurDial_RPAHighQ2 - fDefDial_RPAHighQ2) > 0.0 ||
-// fabs(fCurDial_RPALowQ2 - fDefDial_RPALowQ2) > 0.0);
-// }
-
-// bool RikRPA::IsHandled(int rwenum) {
-// int curenum = rwenum % 1000;
-// switch (curenum) {
-// case kMINERvA_RikRPA_ApplyRPA:
-// return true;
-// case kMINERvA_RikRPA_LowQ2:
-// return true;
-// case kMINERvA_RikRPA_HighQ2:
-// return true;
-// default:
-// return false;
-// }
-// }
-
-// void RikRPA::SetupRPACalculator(int calcenum) {
-// std::string rwdir = FitPar::GetDataBase() +
-// "/reweight/MINERvA/RikRPA/"; std::string fidir = ""; switch (calcenum)
-// {
-// case kNuMuC12:
-// fidir = "outNievesRPAratio-nu12C-20GeV-20170202.root";
-// break;
-
-// case kNuMuO16:
-// fidir = "outNievesRPAratio-nu16O-20GeV-20170202.root";
-// break;
-
-// case kNuMuAr40:
-// fidir = "outNievesRPAratio-nu40Ar-20GeV-20170202.root";
-// break;
-
-// case kNuMuCa40:
-// fidir = "outNievesRPAratio-nu40Ca-20GeV-20170202.root";
-// break;
-
-// case kNuMuFe56:
-// fidir = "outNievesRPAratio-nu56Fe-20GeV-20170202.root";
-// break;
-
-// case kNuMuBarC12:
-// fidir = "outNievesRPAratio-anu12C-20GeV-20170202.root";
-// break;
-
-// case kNuMuBarO16:
-// fidir = "outNievesRPAratio-anu16O-20GeV-20170202.root";
-// break;
-
-// case kNuMuBarAr40:
-// fidir = "outNievesRPAratio-anu40Ar-20GeV-20170202.root";
-// break;
-
-// case kNuMuBarCa40:
-// fidir = "outNievesRPAratio-anu40Ca-20GeV-20170202.root";
-// break;
-
-// case kNuMuBarFe56:
-// fidir = "outNievesRPAratio-anu56Fe-20GeV-20170202.root";
-// break;
-// }
-
-// LOG(FIT) << "Loading RPA CALC : " << fidir << std::endl;
-// TDirectory* olddir = gDirectory;
-
-// fRPACalculators[calcenum] = new weightRPA(rwdir + "/" + fidir);
-// olddir->cd();
-// return;
-// }
-
-// int RikRPA::GetRPACalcEnum(int bpdg, int tpdg) {
-// if (bpdg == 14 && tpdg == 1000060120)
-// return kNuMuC12;
-// else if (bpdg == 14 && tpdg == 1000080160)
-// return kNuMuO16;
-// else if (bpdg == 14 && tpdg == 1000180400)
-// return kNuMuAr40;
-// else if (bpdg == 14 && tpdg == 1000200400)
-// return kNuMuCa40;
-// else if (bpdg == 14 && tpdg == 1000280560)
-// return kNuMuFe56;
-// else if (bpdg == -14 && tpdg == 1000060120)
-// return kNuMuBarC12;
-// else if (bpdg == -14 && tpdg == 1000080160)
-// return kNuMuBarO16;
-// else if (bpdg == -14 && tpdg == 1000180400)
-// return kNuMuBarAr40;
-// else if (bpdg == -14 && tpdg == 1000200400)
-// return kNuMuBarCa40;
-// else if (bpdg == -14 && tpdg == 1000280560)
-// return kNuMuBarFe56;
-// else {
-// ERROR(WRN, "Unknown beam and target combination for RPA Calcs! "
-// << bpdg << " " << tpdg);
-// }
-
-// return -1;
-// }
-// } // namespace reweight
-// } // namespace nuisance
-
-// #endif // __GENIE_ENABLED__
-// #endif //__MINERVA_RW_ENABLED__
+} // namespace reweight
+} // namespace nuisance
+#endif
+#endif
diff --git a/src/Reweight/MINERvAWeightCalcs.h b/src/Reweight/MINERvAWeightCalcs.h
index ec059dc..8bcba23 100644
--- a/src/Reweight/MINERvAWeightCalcs.h
+++ b/src/Reweight/MINERvAWeightCalcs.h
@@ -1,96 +1,117 @@
#ifndef MINERVA_WEIGHT_CALCS
#define MINERVA_WEIGHT_CALCS
#include
-// #ifdef __MINERVA_RW_ENABLED__
-// #ifdef __GENIE_ENABLED__
+#ifdef __MINERVA_RW_ENABLED__
+#ifdef __GENIE_ENABLED__
#include "Conventions/Units.h"
#include "EVGCore/EventRecord.h"
#include "FitEvent.h"
#include "FitParameters.h"
#include "GHEP/GHepParticle.h"
#include "GHEP/GHepRecord.h"
#include "GHEP/GHepUtils.h"
#include "GeneralUtils.h"
#include "NUISANCESyst.h"
#include "NUISANCEWeightCalcs.h"
#include "Ntuple/NtpMCEventRecord.h"
#include "PDG/PDGUtils.h"
#include "WeightUtils.h"
// #include "weightRPA.h"
using namespace genie;
class BaseFitEvt;
namespace nuisance {
namespace reweight {
// MEC Dials
class MINERvAReWeight_MEC : public NUISANCEWeightCalc {
public:
MINERvAReWeight_MEC();
virtual ~MINERvAReWeight_MEC();
double CalcWeight(BaseFitEvt* evt);
void SetDialValue(std::string name, double val);
void SetDialValue(int rwenum, double val);
bool IsHandled(int rwenum);
double fTwk_NormCCMEC;
double fCur_NormCCMEC;
double fDef_NormCCMEC;
-}
-
-// /// RPA Weight Calculator that applies RPA systematics
-// /// to GENIE events. GENIE EVENTS ONLY!
-// class RikRPA : public NUISANCEWeightCalc {
-// public:
-// RikRPA();
-// ~RikRPA();
-
-// double CalcWeight(BaseFitEvt* evt);
-// void SetDialValue(std::string name, double val);
-// void SetDialValue(int rwenum, double val);
-// bool IsHandled(int rwenum);
-
-// void SetupRPACalculator(int calcenum);
-// int GetRPACalcEnum(int bpdg, int tpdg);
-
-// bool fApplyDial_RPACorrection;
-
-// double fTwkDial_RPALowQ2;
-// double fDefDial_RPALowQ2;
-// double fCurDial_RPALowQ2;
-// double fErrDial_RPALowQ2;
-
-// double fTwkDial_RPAHighQ2;
-// double fDefDial_RPAHighQ2;
-// double fCurDial_RPAHighQ2;
-// double fErrDial_RPAHighQ2;
-
-// double* fEventWeights;
-// bool fTweaked;
-
-// const static int kMaxCalculators = 10;
-// enum rpacalcenums {
-// kNuMuC12,
-// kNuMuO16,
-// kNuMuAr40,
-// kNuMuCa40,
-// kNuMuFe56,
-// kNuMuBarC12,
-// kNuMuBarO16,
-// kNuMuBarAr40,
-// kNuMuBarCa40,
-// kNuMuBarFe56
-// };
-// weightRPA* fRPACalculators[kMaxCalculators];
-// };
+ bool fTweaked;
+
+};
+
+// RES Dials
+class MINERvAReWeight_RES : public NUISANCEWeightCalc {
+ public:
+ MINERvAReWeight_RES();
+ virtual ~MINERvAReWeight_RES();
+
+ double CalcWeight(BaseFitEvt* evt);
+ void SetDialValue(std::string name, double val);
+ void SetDialValue(int rwenum, double val);
+ bool IsHandled(int rwenum);
+
+ double fTwk_NormCCRES;
+ double fCur_NormCCRES;
+ double fDef_NormCCRES;
+ bool fTweaked;
+
+ };
+
+/// RPA Weight Calculator that applies RPA systematics
+/// to GENIE events. GENIE EVENTS ONLY!
+class RikRPA : public NUISANCEWeightCalc {
+public:
+ RikRPA();
+ ~RikRPA();
+
+ double CalcWeight(BaseFitEvt* evt);
+ void SetDialValue(std::string name, double val);
+ void SetDialValue(int rwenum, double val);
+ bool IsHandled(int rwenum);
+
+ void SetupRPACalculator(int calcenum);
+ int GetRPACalcEnum(int bpdg, int tpdg);
+
+ bool fApplyDial_RPACorrection;
+
+ double fTwkDial_RPALowQ2;
+ double fDefDial_RPALowQ2;
+ double fCurDial_RPALowQ2;
+ double fErrDial_RPALowQ2;
+
+ double fTwkDial_RPAHighQ2;
+ double fDefDial_RPAHighQ2;
+ double fCurDial_RPAHighQ2;
+ double fErrDial_RPAHighQ2;
+
+ double* fEventWeights;
+ bool fTweaked;
+
+ const static int kMaxCalculators = 10;
+ enum rpacalcenums {
+ kNuMuC12,
+ kNuMuO16,
+ kNuMuAr40,
+ kNuMuCa40,
+ kNuMuFe56,
+ kNuMuBarC12,
+ kNuMuBarO16,
+ kNuMuBarAr40,
+ kNuMuBarCa40,
+ kNuMuBarFe56
+ };
+ // weightRPA* fRPACalculators[kMaxCalculators];
+};
+
}; // namespace reweight
}; // namespace nuisance
-// #endif // __GENIE_ENABLED__
-// #endif //__MINERVA_RW_ENABLED__
+#endif // __GENIE_ENABLED__
+#endif //__MINERVA_RW_ENABLED__
#endif
diff --git a/src/Reweight/NUISANCESyst.cxx b/src/Reweight/NUISANCESyst.cxx
index 44366e1..40aef5e 100644
--- a/src/Reweight/NUISANCESyst.cxx
+++ b/src/Reweight/NUISANCESyst.cxx
@@ -1,52 +1,60 @@
#include "NUISANCESyst.h"
int Reweight::ConvertNUISANCEDial(std::string type) {
for (int i = kUnkownNUISANCEDial + 1; i < kNUISANCEDial_LAST; i++) {
if (!type.compare(ConvNUISANCEDial(i).c_str())) {
return i;
}
}
return kUnkownNUISANCEDial;
};
std::string Reweight::ConvNUISANCEDial(int type) {
switch (type) {
case kGaussianCorr_CCQE_norm: { return "GaussianCorr_CCQE_norm"; }
case kGaussianCorr_CCQE_tilt: { return "GaussianCorr_CCQE_tilt"; }
case kGaussianCorr_CCQE_Pq0: { return "GaussianCorr_CCQE_Pq0"; }
case kGaussianCorr_CCQE_Wq0: { return "GaussianCorr_CCQE_Wq0"; }
case kGaussianCorr_CCQE_Pq3: { return "GaussianCorr_CCQE_Pq3"; }
case kGaussianCorr_CCQE_Wq3: { return "GaussianCorr_CCQE_Wq3"; }
case kGaussianCorr_2p2h_norm: { return "GaussianCorr_2p2h_norm"; }
case kGaussianCorr_2p2h_tilt: { return "GaussianCorr_2p2h_tilt"; }
case kGaussianCorr_2p2h_Pq0: { return "GaussianCorr_2p2h_Pq0"; }
case kGaussianCorr_2p2h_Wq0: { return "GaussianCorr_2p2h_Wq0"; }
case kGaussianCorr_2p2h_Pq3: { return "GaussianCorr_2p2h_Pq3"; }
case kGaussianCorr_2p2h_Wq3: { return "GaussianCorr_2p2h_Wq3"; }
case kGaussianCorr_2p2h_PPandNN_norm: { return "GaussianCorr_2p2h_PPandNN_norm"; }
case kGaussianCorr_2p2h_PPandNN_tilt: { return "GaussianCorr_2p2h_PPandNN_tilt"; }
case kGaussianCorr_2p2h_PPandNN_Pq0: { return "GaussianCorr_2p2h_PPandNN_Pq0"; }
case kGaussianCorr_2p2h_PPandNN_Wq0: { return "GaussianCorr_2p2h_PPandNN_Wq0"; }
case kGaussianCorr_2p2h_PPandNN_Pq3: { return "GaussianCorr_2p2h_PPandNN_Pq3"; }
case kGaussianCorr_2p2h_PPandNN_Wq3: { return "GaussianCorr_2p2h_PPandNN_Wq3"; }
case kGaussianCorr_2p2h_NP_norm: { return "GaussianCorr_2p2h_NP_norm"; }
case kGaussianCorr_2p2h_NP_tilt: { return "GaussianCorr_2p2h_NP_tilt"; }
case kGaussianCorr_2p2h_NP_Pq0: { return "GaussianCorr_2p2h_NP_Pq0"; }
case kGaussianCorr_2p2h_NP_Wq0: { return "GaussianCorr_2p2h_NP_Wq0"; }
case kGaussianCorr_2p2h_NP_Pq3: { return "GaussianCorr_2p2h_NP_Pq3"; }
case kGaussianCorr_2p2h_NP_Wq3: { return "GaussianCorr_2p2h_NP_Wq3"; }
case kGaussianCorr_CC1pi_norm: { return "GaussianCorr_CC1pi_norm"; }
case kGaussianCorr_CC1pi_tilt: { return "GaussianCorr_CC1pi_tilt"; }
case kGaussianCorr_CC1pi_Pq0: { return "GaussianCorr_CC1pi_Pq0"; }
case kGaussianCorr_CC1pi_Wq0: { return "GaussianCorr_CC1pi_Wq0"; }
case kGaussianCorr_CC1pi_Pq3: { return "GaussianCorr_CC1pi_Pq3"; }
case kGaussianCorr_CC1pi_Wq3: { return "GaussianCorr_CC1pi_Wq3"; }
case kGaussianCorr_AllowSuppression: { return "GaussianCorr_AllowSuppression"; }
+
case kModeNorm_NormRES: { return "NormRES"; }
+ case kMINERvARW_NormCCMEC: { return "MINERvARW_NormCCMEC"; }
+ case kMINERvARW_NormCCRES: { return "MINERvARW_NormCCRES"; }
+
+ case kMINERvARW_RikRPA_ApplyRPA: { return "MINERvARW_RikRPA_ApplyRPA"; }
+ case kMINERvARW_RikRPA_LowQ2: { return "MINERvARW_RikRPA_LowQ2"; }
+ case kMINERvARW_RikRPA_HighQ2: { return "MINERvARW_RikRPA_HighQ2"; }
+
default: return "unknown_nuisance_dial";
}
};
diff --git a/src/Reweight/NUISANCESyst.h b/src/Reweight/NUISANCESyst.h
index 296c956..a3cfb74 100644
--- a/src/Reweight/NUISANCESyst.h
+++ b/src/Reweight/NUISANCESyst.h
@@ -1,50 +1,57 @@
#ifndef NUISANCESyst_H
#define NUISANCESyst_H
#include "GeneralUtils.h"
namespace Reweight {
enum NUISANCESyst {
kUnkownNUISANCEDial = 0,
kGaussianCorr_CCQE_norm,
kGaussianCorr_CCQE_tilt,
kGaussianCorr_CCQE_Pq0,
kGaussianCorr_CCQE_Wq0,
kGaussianCorr_CCQE_Pq3,
kGaussianCorr_CCQE_Wq3,
kGaussianCorr_2p2h_norm,
kGaussianCorr_2p2h_tilt,
kGaussianCorr_2p2h_Pq0,
kGaussianCorr_2p2h_Wq0,
kGaussianCorr_2p2h_Pq3,
kGaussianCorr_2p2h_Wq3,
kGaussianCorr_2p2h_PPandNN_norm,
kGaussianCorr_2p2h_PPandNN_tilt,
kGaussianCorr_2p2h_PPandNN_Pq0,
kGaussianCorr_2p2h_PPandNN_Wq0,
kGaussianCorr_2p2h_PPandNN_Pq3,
kGaussianCorr_2p2h_PPandNN_Wq3,
kGaussianCorr_2p2h_NP_norm,
kGaussianCorr_2p2h_NP_tilt,
kGaussianCorr_2p2h_NP_Pq0,
kGaussianCorr_2p2h_NP_Wq0,
kGaussianCorr_2p2h_NP_Pq3,
kGaussianCorr_2p2h_NP_Wq3,
kGaussianCorr_CC1pi_norm,
kGaussianCorr_CC1pi_tilt,
kGaussianCorr_CC1pi_Pq0,
kGaussianCorr_CC1pi_Wq0,
kGaussianCorr_CC1pi_Pq3,
kGaussianCorr_CC1pi_Wq3,
kGaussianCorr_AllowSuppression,
kModeNorm_NormRES,
+ kMINERvARW_NormCCMEC,
+ kMINERvARW_NormCCRES,
+
+ kMINERvARW_RikRPA_ApplyRPA,
+ kMINERvARW_RikRPA_LowQ2,
+ kMINERvARW_RikRPA_HighQ2,
+
kNUISANCEDial_LAST
};
int ConvertNUISANCEDial(std::string type);
std::string ConvNUISANCEDial(int type);
};
#endif
diff --git a/src/Reweight/NUISANCEWeightEngine.cxx b/src/Reweight/NUISANCEWeightEngine.cxx
index e7344ca..116385c 100644
--- a/src/Reweight/NUISANCEWeightEngine.cxx
+++ b/src/Reweight/NUISANCEWeightEngine.cxx
@@ -1,129 +1,110 @@
#include "NUISANCEWeightEngine.h"
#include "NUISANCEWeightCalcs.h"
-NUISANCEWeightEngine::NUISANCEWeightEngine(std::string name) {
+#ifdef __MINERVA_RW_ENABLED__
+#ifdef __GENIE_ENABLED__
+#include "MINERvAWeightCalcs.h"
+#endif
+#endif
+
+NUISANCEWeightEngine::NUISANCEWeightEngine(std::string name) {
// Setup the NUISANCE Reweight engine
fCalcName = name;
LOG(FIT) << "Setting up NUISANCE Custom RW : " << fCalcName << std::endl;
// Load in all Weight Calculations
- fWeightCalculators.push_back( new GaussianModeCorr() );
- fWeightCalculators.push_back( new ModeNormCalc() );
+ fWeightCalculators.push_back(new GaussianModeCorr());
+ fWeightCalculators.push_back(new ModeNormCalc());
+
+#ifdef __MINERVA_RW_ENABLED__
+#ifdef __GENIE_ENABLED__
+ fWeightCalculators.push_back(new nuisance::reweight::MINERvAReWeight_MEC());
+ fWeightCalculators.push_back(new nuisance::reweight::MINERvAReWeight_RES());
+ fWeightCalculators.push_back(new nuisance::reweight::RikRPA());
+#endif
+#endif
+
// Set Abs Twk Config
fIsAbsTwk = true;
-
};
-
void NUISANCEWeightEngine::IncludeDial(std::string name, double startval) {
-
// Get First enum
int nuisenum = Reweight::ConvDial(name, kCUSTOM);
// Setup Maps
- fEnumIndex[nuisenum];// = std::vector(0);
- fNameIndex[name]; // = std::vector(0);
+ 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
- int singleenum = Reweight::ConvDial(singlename, kCUSTOM);
+ int singleenum = Reweight::ConvDial(singlename, kCUSTOM);
// Fill Maps
int index = fValues.size();
fValues.push_back(0.0);
fNUISANCEEnums.push_back(singleenum);
// Initialize dial
std::cout << "Registering " << singlename << " from " << name << std::endl;
// Setup index
fEnumIndex[nuisenum].push_back(index);
fNameIndex[name].push_back(index);
}
// Set Value if given
if (startval != -999.9) {
SetDialValue(nuisenum, startval);
}
};
-
void NUISANCEWeightEngine::SetDialValue(int nuisenum, double val) {
std::vector indices = fEnumIndex[nuisenum];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
}
}
void NUISANCEWeightEngine::SetDialValue(std::string name, double val) {
std::vector indices = fNameIndex[name];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
}
}
void NUISANCEWeightEngine::Reconfigure(bool silent) {
-
for (size_t i = 0; i < fNUISANCEEnums.size(); i++) {
-
- for (std::vector::iterator calciter = fWeightCalculators.begin();
+ for (std::vector::iterator calciter =
+ fWeightCalculators.begin();
calciter != fWeightCalculators.end(); calciter++) {
-
- NUISANCEWeightCalc* nuiscalc = static_cast(*calciter);
+ NUISANCEWeightCalc* nuiscalc =
+ static_cast(*calciter);
if (nuiscalc->IsHandled(fNUISANCEEnums[i])) {
nuiscalc->SetDialValue(fNUISANCEEnums[i], fValues[i]);
}
}
}
}
-
-
-double NUISANCEWeightEngine::CalcWeight(BaseFitEvt * evt) {
+double NUISANCEWeightEngine::CalcWeight(BaseFitEvt* evt) {
double rw_weight = 1.0;
// Cast as usable class
- for (std::vector::iterator iter = fWeightCalculators.begin();
+ for (std::vector::iterator iter =
+ fWeightCalculators.begin();
iter != fWeightCalculators.end(); iter++) {
NUISANCEWeightCalc* nuiscalc = static_cast(*iter);
rw_weight *= nuiscalc->CalcWeight(evt);
}
// Return rw_weight
return rw_weight;
}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-