diff --git a/cmake/NEUTSetup.cmake b/cmake/NEUTSetup.cmake
index ac6fc40..5c8c2ce 100644
--- a/cmake/NEUTSetup.cmake
+++ b/cmake/NEUTSetup.cmake
@@ -1,252 +1,261 @@
# Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
include(cmake/parseConfigApp.cmake)
find_program(NEUTCONFIG NAMES neut-config)
LIST(APPEND EXTRA_CXX_FLAGS -DNEED_FILL_NEUT_COMMONS)
SET(HAVENEUTCONFIG FALSE)
# We are dealing with shiny NEUT
if(NOT "${NEUTCONFIG}" STREQUAL "NEUTCONFIG-NOTFOUND")
SET(HAVENEUTCONFIG TRUE)
cmessage(STATUS "Found neut-config, using it to determine configuration.")
else()
cmessage(STATUS "Failed to find neut-config, assuming older NEUT build.")
endif()
if(HAVENEUTCONFIG)
execute_process (COMMAND neut-config
--version OUTPUT_VARIABLE NEUT_VER
OUTPUT_STRIP_TRAILING_WHITESPACE)
execute_process (COMMAND neut-config
--incdir OUTPUT_VARIABLE NEUT_INCLUDE_DIRS
OUTPUT_STRIP_TRAILING_WHITESPACE)
execute_process (COMMAND neut-config
--libdir OUTPUT_VARIABLE NEUT_LINK_DIRS
OUTPUT_STRIP_TRAILING_WHITESPACE)
GETLIBDIRS(neut-config --cernflags CERN_LIB_DIR)
LIST(APPEND NEUT_LINK_DIRS ${CERN_LIB_DIR})
GETLIBS(neut-config --cernflags CERN_LIBS)
- if(USE_REWEIGHT)
+ if(NOT DEFINED USE_NEUT_REWEIGHT)
+ if(USE_REWEIGHT)
+ SET(USE_NEUT_REWEIGHT ON)
+ else()
+ SET(USE_NEUT_REWEIGHT OFF)
+ endif()
+ endif()
+
+ if(USE_REWEIGHT AND USE_NEUT_REWEIGHT)
execute_process (COMMAND neut-config
--rwlibflags OUTPUT_VARIABLE NEUT_RWLIBS
OUTPUT_STRIP_TRAILING_WHITESPACE)
GETLIBS(neut-config --rwlibflags NEUT_RWLIBS)
LIST(APPEND NEUT_LIBS ${NEUT_RWLIBS})
+ LIST(APPEND EXTRA_CXX_FLAGS ${NEUT_INCLUDE_DIRS} -D__USE_NEUT_REWEIGHT__)
else()
GETLIBS(neut-config --libflags NEUT_GENLIBS)
GETLIBS(neut-config --iolibflags NEUT_LIBS)
LIST(APPEND NEUT_LIBS ${NEUT_IOLIBS})
LIST(APPEND NEUT_LIBS ${NEUT_GENLIBS})
endif()
LIST(APPEND NEUT_LIBS ${CERN_LIBS};gfortran)
LIST(APPEND EXTRA_LIBS ${NEUT_LIBS})
string(REPLACE "." "" NEUT_VERSION ${NEUT_VER})
PrefixList(NEUT_INCLUDE_DIRS "-I" ${NEUT_INCLUDE_DIRS})
LIST(APPEND EXTRA_CXX_FLAGS ${NEUT_INCLUDE_DIRS} -D__NEUT_ENABLED__ -D__NEUT_VERSION__=${NEUT_VERSION})
LIST(APPEND EXTRA_LINK_DIRS ${NEUT_LINK_DIRS})
include(CheckCXXCompilerFlag)
CHECK_CXX_COMPILER_FLAG(-Wl,--allow-multiple-definition COMPILER_SUPPORTS_ALLOW_MULTIPLE_DEFINITION)
IF(COMPILER_SUPPORTS_ALLOW_MULTIPLE_DEFINITION)
LIST(APPEND EXTRA_LINK_FLAGS -Wl,--allow-multiple-definition)
ENDIF()
CHECK_CXX_COMPILER_FLAG(-no-pie COMPILER_SUPPORTS_NO_PIE)
CHECK_CXX_COMPILER_FLAG(-fno-pie COMPILER_SUPPORTS_FNO_PIE)
CHECK_CXX_COMPILER_FLAG(-fno-PIE COMPILER_SUPPORTS_FNO_PIE_CAP)
if(COMPILER_SUPPORTS_NO_PIE)
set(PIE_FLAGS "-no-pie")
elseif(COMPILER_SUPPORTS_FNO_PIE)
set(PIE_FLAGS "-fno-pie")
elseif(COMPILER_SUPPOERTS_FNO_PIE_CAP)
set(PIE_FLAGS "-fno-PIE")
else()
message(STATUS "The compiler ${CMAKE_CXX_COMPILER} has no C++11 support. Please use a different C++ compiler.")
endif()
LIST(APPEND EXTRA_EXE_FLAGS
${PIE_FLAGS})
cmessage(STATUS "NEUT")
cmessage(STATUS " Version : ${NEUT_VER}")
cmessage(STATUS " Flags : ${NEUT_CXX_FLAGS}")
cmessage(STATUS " Includes : ${NEUT_INCLUDE_DIRS}")
cmessage(STATUS " Link Dirs : ${NEUT_LINK_DIRS}")
cmessage(STATUS " Libs : ${NEUT_LIBS}")
cmessage(STATUS " Exe Flags : ${EXTRA_EXE_FLAGS}")
else() # Everything better be set up already
if(NEUT_ROOT STREQUAL "")
cmessage(FATAL_ERROR "Variable NEUT_ROOT is not defined. Please export environment variable NEUT_ROOT or configure with -DNEUT_ROOT=/path/to/NEUT. This must be set to point to a prebuilt NEUT instance.")
endif()
if(CERN STREQUAL "")
cmessage(FATAL_ERROR "Variable CERN is not defined. Please export environment variable CERN or configure with -DCERN=/path/to/CERNLIB. This must be set to point to a prebuilt CERNLIB instance.")
endif()
if(CERN_LEVEL STREQUAL "")
cmessage(FATAL_ERROR "Variable CERN_LEVEL is not defined. Please export environment variable CERN_LEVEL or configure with -DCERN_LEVEL=XXXX (likely to be 2005).")
endif()
if(${NEUT_VERSION} VERSION_LESS 5.4.0)
set(NEUT_LIB_DIR ${NEUT_ROOT}/lib/Linux_pc)
else()
set(NEUT_LIB_DIR ${NEUT_ROOT}/lib)
endif()
set(NEUT_CLASS ${NEUT_ROOT}/src/neutclass)
LIST(APPEND EXTRA_CXX_FLAGS -D__NEUT_ENABLED__ -DNEUT_VERSION=${NEUT_VERSION})
LIST(APPEND EXTRA_CXX_FLAGS
-I${NEUT_ROOT}/include
-I${NEUT_ROOT}/src/neutclass)
LIST(APPEND EXTRA_LINK_DIRS
${NEUT_LIB_DIR}
${CERN}/${CERN_LEVEL}/lib)
if(USE_REWEIGHT)
LIST(APPEND EXTRA_CXX_FLAGS
-I${NEUT_ROOT}/src/reweight)
LIST(APPEND EXTRA_LINK_DIRS
${NEUT_ROOT}/src/reweight)
endif()
if(${NEUT_VERSION} VERSION_EQUAL 5.4.2)
LIST(APPEND EXTRA_LIBS
-Wl,--as-needed)
if(USE_REWEIGHT)
LIST(APPEND EXTRA_LIBS
NReWeight)
endif()
LIST(APPEND EXTRA_LIBS
-Wl,--start-group
neutcore_5.4.2
nuccorspl_5.4.2 #typo in NEUT, may hopefully disappear
nuceff_5.4.2
partnuck_5.4.2
skmcsvc_5.4.2
tauola_5.4.2
HT2p2h_5.4.0
N1p1h_5.4.0
-Wl,--end-group
jetset74
pdflib804
mathlib
packlib
pawlib)
LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_COMMON_QEAV)
elseif(${NEUT_VERSION} VERSION_EQUAL 5.4.0)
LIST(APPEND EXTRA_LIBS
-Wl,--as-needed)
if(USE_REWEIGHT)
LIST(APPEND EXTRA_LIBS
NReWeight)
endif()
LIST(APPEND EXTRA_LIBS
-Wl,--start-group
neutcore_5.4.0
nuccorspl_5.4.0 #typo in NEUT, may hopefully disappear
nuceff_5.4.0
partnuck_5.4.0
skmcsvc_5.4.0
tauola_5.4.0
HT2p2h_5.4.0
N1p1h_5.4.0
specfunc_5.4.0
radcorr_5.4.0
gfortran
-Wl,--end-group
jetset74
pdflib804
mathlib
packlib
pawlib)
else()
LIST(APPEND EXTRA_LIBS
-Wl,--as-needed)
if(USE_REWEIGHT)
LIST(APPEND EXTRA_LIBS
NReWeight)
endif()
LIST(APPEND EXTRA_LIBS
-Wl,--start-group
neutcore
nuccorrspl
nuceff
partnuck
skmcsvc
tauola
-Wl,--end-group
jetset74
pdflib804
mathlib
packlib
pawlib)
endif()
set(NEUT_ROOT_LIBS)
LIST(APPEND NEUT_ROOT_LIBS
${NEUT_CLASS}/neutctrl.so
${NEUT_CLASS}/neutfsivert.so)
# Check for new versions of NEUT with NUCLEON FSI
if(EXISTS "${NEUT_CLASS}/neutnucfsistep.so")
set(NEUT_NUCFSI 1)
LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_NUCFSI_ENABLED)
LIST(APPEND NEUT_ROOT_LIBS
${NEUT_CLASS}/neutnucfsistep.so
${NEUT_CLASS}/neutnucfsivert.so
)
endif()
if(${NEUT_VERSION} VERSION_LESS 5.4.0)
LIST(APPEND NEUT_ROOT_LIBS
${NEUT_CLASS}/neutrootTreeSingleton.so)
endif()
LIST(APPEND NEUT_ROOT_LIBS
${NEUT_CLASS}/neutvtx.so
${NEUT_CLASS}/neutfsipart.so
${NEUT_CLASS}/neutpart.so
${NEUT_CLASS}/neutvect.so
)
foreach(OBJ ${NEUT_ROOT_LIBS})
LIST(APPEND EXTRA_SHAREDOBJS ${OBJ})
endforeach()
endif()
diff --git a/cmake/T2KSetup.cmake b/cmake/T2KSetup.cmake
index a4d3028..19c1582 100644
--- a/cmake/T2KSetup.cmake
+++ b/cmake/T2KSetup.cmake
@@ -1,160 +1,160 @@
# Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
include(cmake/parseConfigApp.cmake)
find_program(T2KRWCONFIG NAMES t2kreweight-config)
SET(HAVET2KRWCONFIG FALSE)
# We are dealing with shiny NEUT
if(NOT "${T2KRWCONFIG}" STREQUAL "T2KRWCONFIG-NOTFOUND")
SET(HAVET2KRWCONFIG TRUE)
- cmessage(STATUS "Found neut-config, using it to determine configuration.")
+ cmessage(STATUS "Found t2kreweight-config, using it to determine configuration.")
else()
- cmessage(STATUS "Failed to find neut-config, assuming older NEUT build.")
+ cmessage(STATUS "Failed to find t2kreweight-config, assuming older NEUT build.")
endif()
if(HAVET2KRWCONFIG)
if(NOT DEFINED CMAKE_CXX_STANDARD OR "${CMAKE_CXX_STANDARD} " STREQUAL " ")
SET(CMAKE_CXX_STANDARD 11)
endif()
LIST(APPEND EXTRA_CXX_FLAGS -DT2KRW_OA2021_INTERFACE -D__T2KREW_ENABLED__)
GETLIBDIRS(t2kreweight-config --linkflags T2KRW_LINK_DIRS)
GETINCDIRS(t2kreweight-config --cflags T2KRW_INC_DIRS)
GETLIBS(t2kreweight-config --linkflags T2KRW_LIBS)
cmessage(STATUS "T2KReWeight:")
cmessage(STATUS " LINK DIRS: ${T2KRW_LINK_DIRS}")
cmessage(STATUS " INC DIRS: ${T2KRW_INC_DIRS}")
cmessage(STATUS " LIBS: ${T2KRW_LIBS}")
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${T2KRW_INC_DIRS})
LIST(APPEND EXTRA_LINK_DIRS ${T2KRW_LINK_DIRS})
LIST(APPEND EXTRA_LIBS ${T2KRW_LIBS})
execute_process (COMMAND t2kreweight-config
--has-feature NIWG RESULT_VARIABLE T2KRW_HAS_NIWG)
if("${T2KRW_HAS_NIWG} " STREQUAL "0 ")
GETLIBDIRS(t2kreweight-config --niwgflags NIWG_LINK_DIRS)
GETINCDIRS(t2kreweight-config --niwgflags NIWG_INC_DIRS)
GETLIBS(t2kreweight-config --niwgflags NIWG_LIBS)
cmessage(STATUS "NIWG:")
cmessage(STATUS " LINK DIRS: ${NIWG_LINK_DIRS}")
cmessage(STATUS " INC DIRS: ${NIWG_INC_DIRS}")
cmessage(STATUS " LIBS: ${NIWG_LIBS}")
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${NIWG_INC_DIRS})
LIST(APPEND EXTRA_LINK_DIRS ${NIWG_LINK_DIRS})
LIST(APPEND EXTRA_LIBS ${NIWG_LIBS})
endif()
execute_process (COMMAND t2kreweight-config
--has-feature NEUT RESULT_VARIABLE T2KRW_HAS_NEUT)
if("${T2KRW_HAS_NEUT} " STREQUAL "0 ")
LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_EVENT_ENABLED)
set(USE_NEUT_EVENT TRUE)
set(USE_NEUT_EVENT TRUE CACHE BOOL "Whether to enable NEUT (event i/o) support. Requires external libraries. " FORCE)
GETLIBDIRS(t2kreweight-config --neutflags NEUT_LINK_DIRS)
GETINCDIRS(t2kreweight-config --neutflags NEUT_INC_DIRS)
GETLIBS(t2kreweight-config --neutflags NEUT_LIBS)
cmessage(STATUS "NEUT:")
cmessage(STATUS " LINK DIRS: ${NEUT_LINK_DIRS}")
cmessage(STATUS " INC DIRS: ${NEUT_INC_DIRS}")
cmessage(STATUS " LIBS: ${NEUT_LIBS}")
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${NEUT_INC_DIRS})
LIST(APPEND EXTRA_LINK_DIRS ${NEUT_LINK_DIRS})
LIST(APPEND EXTRA_LIBS ${NEUT_LIBS})
endif()
execute_process (COMMAND t2kreweight-config
--has-feature GEANT RESULT_VARIABLE T2KRW_HAS_GEANT)
if("${T2KRW_HAS_GEANT} " STREQUAL "0 ")
GETLIBDIRS(t2kreweight-config --geantflags GEANT_LINK_DIRS)
GETINCDIRS(t2kreweight-config --geantflags GEANT_INC_DIRS)
GETLIBS(t2kreweight-config --geantflags GEANT_LIBS)
cmessage(STATUS "GEANT:")
cmessage(STATUS " LINK DIRS: ${GEANT_LINK_DIRS}")
cmessage(STATUS " INC DIRS: ${GEANT_INC_DIRS}")
cmessage(STATUS " LIBS: ${GEANT_LIBS}")
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${GEANT_INC_DIRS})
LIST(APPEND EXTRA_LINK_DIRS ${GEANT_LINK_DIRS})
LIST(APPEND EXTRA_LIBS ${GEANT_LIBS})
endif()
execute_process (COMMAND t2kreweight-config
--has-feature ND280 RESULT_VARIABLE T2KRW_HAS_ND280)
if("${T2KRW_HAS_ND280} " STREQUAL "0 ")
GETLIBDIRS(t2kreweight-config --nd280flags ND280_LINK_DIRS)
GETINCDIRS(t2kreweight-config --nd280flags ND280_INC_DIRS)
GETLIBS(t2kreweight-config --nd280flags ND280_LIBS)
cmessage(STATUS "ND280:")
cmessage(STATUS " LINK DIRS: ${ND280_LINK_DIRS}")
cmessage(STATUS " INC DIRS: ${ND280_INC_DIRS}")
cmessage(STATUS " LIBS: ${ND280_LIBS}")
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${ND280_INC_DIRS})
LIST(APPEND EXTRA_LINK_DIRS ${ND280_LINK_DIRS})
LIST(APPEND EXTRA_LIBS ${ND280_LIBS})
endif()
else()
if(T2KREWEIGHT STREQUAL "")
cmessage(FATAL_ERROR "Variable T2KREWEIGHT is not defined. Either configure with -DT2KREWEIGHT or \"\$ export T2KREWEIGHT=/path/to/T2KReWeight\". This must be set to point to a prebuilt T2KReWeight instance.")
endif()
LIST(APPEND EXTRA_CXX_FLAGS -D__T2KREW_ENABLED__ )
# First check the OAANALYSIS libs (need to grab some headers for T2KReWeight linking if compiled with oaAnalysisReader)
IF(NOT $ENV{OAANALYSISREADERROOT} STREQUAL "")
cmessage(STATUS "Found OAANALYSISREADERROOT $ENV{OAANALYSISREADERROOT}, appending...")
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES $ENV{OAANALYSISREADERROOT}/$ENV{OAANALYSISREADERCONFIG})
LIST(APPEND EXTRA_LINK_DIRS $ENV{OAANALYSISREADERROOT}/$ENV{OAANALYSISREADERCONFIG})
LIST(APPEND EXTRA_LIBS oaAnalysisReader)
# Don't have to append this; should be appeneded in ${T2KReWeight/src/T2KBuild.h}
#LIST(APPEND EXTRA_CXX_FLAGS -D__T2KRW_OAANALYSIS_ENABLED__)
endif()
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${T2KREWEIGHT}/src/)
LIST(APPEND EXTRA_LINK_DIRS ${T2KREWEIGHT}/lib)
LIST(APPEND EXTRA_LIBS T2KReWeight)
endif()
\ No newline at end of file
diff --git a/src/Reweight/FitWeight.cxx b/src/Reweight/FitWeight.cxx
index be3b57d..dfbc55d 100644
--- a/src/Reweight/FitWeight.cxx
+++ b/src/Reweight/FitWeight.cxx
@@ -1,305 +1,314 @@
#include "FitWeight.h"
#include "GENIEWeightEngine.h"
#include "LikelihoodWeightEngine.h"
#include "ModeNormEngine.h"
#include "NEUTWeightEngine.h"
#include "NIWGWeightEngine.h"
#include "NUISANCEWeightEngine.h"
#include "NuWroWeightEngine.h"
#include "OscWeightEngine.h"
#include "SampleNormEngine.h"
#include "SplineWeightEngine.h"
#include "T2KWeightEngine.h"
#ifdef __NOVA_ENABLED__
#include "NOvARwgtEngine.h"
#endif
#ifdef __NUSYST_ENABLED__
#include "nusystematicsWeightEngine.h"
#endif
void FitWeight::AddRWEngine(int type) {
NUIS_LOG(FIT, "Adding reweight engine " << type);
switch (type) {
case kNEUT:
fAllRW[type] = new NEUTWeightEngine("neutrw");
break;
case kNUWRO:
fAllRW[type] = new NuWroWeightEngine("nuwrorw");
break;
case kGENIE:
fAllRW[type] = new GENIEWeightEngine("genierw");
break;
case kNORM:
fAllRW[type] = new SampleNormEngine("normrw");
break;
case kLIKEWEIGHT:
fAllRW[type] = new LikelihoodWeightEngine("likerw");
break;
case kT2K:
fAllRW[type] = new T2KWeightEngine("t2krw");
break;
case kCUSTOM:
fAllRW[type] = new NUISANCEWeightEngine("nuisrw");
break;
case kSPLINEPARAMETER:
fAllRW[type] = new SplineWeightEngine("splinerw");
break;
case kNIWG:
fAllRW[type] = new NIWGWeightEngine("niwgrw");
break;
case kOSCILLATION:
fAllRW[type] = new OscWeightEngine();
break;
case kMODENORM:
fAllRW[type] = new ModeNormEngine();
break;
#ifdef __NOVA_ENABLED__
case kNOvARWGT:
fAllRW[type] = new NOvARwgtEngine();
break;
#endif
#ifdef __NUSYST_ENABLED__
case kNuSystematics:
fAllRW[type] = new nusystematicsWeightEngine();
break;
#endif
default:
NUIS_ABORT("CANNOT ADD RW Engine for unknown dial type: " << type);
break;
}
}
WeightEngineBase *FitWeight::GetRWEngine(int type) {
if (HasRWEngine(type)) {
return fAllRW[type];
}
NUIS_ABORT("CANNOT get RW Engine for dial type: " << type);
}
bool FitWeight::HasRWEngine(int type) {
switch (type) {
case kNEUT:
case kNUWRO:
case kGENIE:
case kNORM:
case kLIKEWEIGHT:
case kT2K:
case kCUSTOM:
case kSPLINEPARAMETER:
case kNIWG:
case kOSCILLATION:
#ifdef __NOVA_ENABLED__
case kNOvARWGT:
#endif
#ifdef __NUSYST_ENABLED__
case kNuSystematics:
#endif
{
return fAllRW.count(type);
}
default: { NUIS_ABORT("CANNOT get RW Engine for dial type: " << type); }
}
}
void FitWeight::IncludeDial(std::string name, std::string type, double val) {
// Should register the dial here.
int typeenum = Reweight::ConvDialType(type);
IncludeDial(name, typeenum, val);
}
void FitWeight::IncludeDial(std::string name, int dialtype, double val) {
// Get the dial enum
int nuisenum = Reweight::ConvDial(name, dialtype);
if (nuisenum == -1) {
NUIS_ERR(FTL, "Could not include dial " << name);
NUIS_ERR(FTL, "With dialtype: " << dialtype);
NUIS_ERR(FTL, "With value: " << val);
NUIS_ABORT("With nuisenum: " << nuisenum);
}
// Setup RW Engine Pointer
if (fAllRW.find(dialtype) == fAllRW.end()) {
AddRWEngine(dialtype);
}
WeightEngineBase *rw = fAllRW[dialtype];
// Include the dial
rw->IncludeDial(name, val);
// Set Dial Value
if (val != -9999.9) {
rw->SetDialValue(name, val);
}
// Sort Maps
fAllEnums[name] = nuisenum;
fAllValues[nuisenum] = val;
// Sort Lists
fNameList.push_back(name);
fEnumList.push_back(nuisenum);
fValueList.push_back(val);
}
void FitWeight::Reconfigure(bool silent) {
// Reconfigure all added RW engines
for (std::map::iterator iter = fAllRW.begin();
iter != fAllRW.end(); iter++) {
(*iter).second->Reconfigure(silent);
}
}
void FitWeight::SetDialValue(std::string name, double val) {
// Add extra check, if name not found look for one with name in it.
int nuisenum = fAllEnums[name];
SetDialValue(nuisenum, val);
}
// Allow for name aswell using GlobalList to determine sample name.
void FitWeight::SetDialValue(int nuisenum, double val) {
// Conv dial type
int dialtype = Reweight::GetDialType(nuisenum);
if (fAllRW.find(dialtype) == fAllRW.end()) {
- NUIS_ERR(FTL, "Can't find RW engine for parameter " << fNameList[dialtype]);
+
+ std::string name = "";
+ for(size_t i = 0; i < fEnumList.size(); ++i){
+ if(fEnumList[i] == nuisenum){
+ name = fNameList[i];
+ break;
+ }
+ }
+
+ NUIS_ERR(FTL, "Can't find RW engine for parameter " << name);
NUIS_ERR(FTL, "With dialtype " << dialtype << ", "
<< Reweight::RemoveDialType(nuisenum));
NUIS_ABORT("Are you sure you enabled the right engines?");
}
// Get RW Engine for this dial
fAllRW[dialtype]->SetDialValue(nuisenum, val);
fAllValues[nuisenum] = val;
// Update ValueList
for (size_t i = 0; i < fEnumList.size(); i++) {
if (fEnumList[i] == nuisenum) {
fValueList[i] = val;
}
}
}
void FitWeight::SetAllDials(const double *x, int n) {
for (size_t i = 0; i < (UInt_t)n; i++) {
int rwenum = fEnumList[i];
SetDialValue(rwenum, x[i]);
}
Reconfigure();
}
double FitWeight::GetDialValue(std::string name) {
// Add extra check, if name not found look for one with name in it.
int nuisenum = fAllEnums[name];
return GetDialValue(nuisenum);
}
double FitWeight::GetDialValue(int nuisenum) { return fAllValues[nuisenum]; }
int FitWeight::GetDialPos(std::string name) {
int rwenum = fAllEnums[name];
return GetDialPos(rwenum);
}
int FitWeight::GetDialPos(int nuisenum) {
for (size_t i = 0; i < fEnumList.size(); i++) {
if (fEnumList[i] == nuisenum) {
return i;
}
}
NUIS_ABORT("No Dial Found! (enum = " << nuisenum << ") ");
}
bool FitWeight::DialIncluded(std::string name) {
return (fAllEnums.find(name) != fAllEnums.end());
}
bool FitWeight::DialIncluded(int rwenum) {
return (fAllValues.find(rwenum) != fAllValues.end());
}
double FitWeight::CalcWeight(BaseFitEvt *evt) {
double rwweight = 1.0;
for (std::map::iterator iter = fAllRW.begin();
iter != fAllRW.end(); iter++) {
double w = (*iter).second->CalcWeight(evt);
rwweight *= w;
}
return rwweight;
}
void FitWeight::UpdateWeightEngine(const double *x) {
size_t count = 0;
for (std::vector::iterator iter = fEnumList.begin();
iter != fEnumList.end(); iter++) {
SetDialValue((*iter), x[count]);
count++;
}
}
void FitWeight::GetAllDials(double *x, int n) {
for (int i = 0; i < n; i++) {
x[i] = GetDialValue(fEnumList[i]);
}
}
// bool FitWeight::NeedsEventReWeight(const double* x) {
// bool haschange = false;
// size_t count = 0;
// // Compare old to new and decide if RW needed.
// for (std::vector::iterator iter = fEnumList.begin();
// iter != fEnumList.end(); iter++) {
// int nuisenum = (*iter);
-// int type = (nuisenum / 1000) - (nuisenum % 1000);
+// int type = (nuisenum / NUIS_DIAL_OFFSET) - (nuisenum % NUIS_DIAL_OFFSET);
// // Compare old to new
// double oldval = GetDialValue(nuisenum);
// double newval = x[count];
// if (oldval != newval) {
// if (fAllRW[type]->NeedsEventReWeight()) {
// haschange = true;
// }
// }
// count++;
// }
// return haschange;
// }
double FitWeight::GetSampleNorm(std::string name) {
if (name.empty())
return 1.0;
// Find norm dial
if (fAllEnums.find(name + "_norm") != fAllEnums.end()) {
if (fAllValues.find(fAllEnums[name + "_norm"]) != fAllValues.end()) {
return fAllValues[fAllEnums[name + "_norm"]];
} else {
return 1.0;
}
} else {
return 1.0;
}
}
void FitWeight::Print() {
NUIS_LOG(REC, "Fit Weight State: ");
for (size_t i = 0; i < fNameList.size(); i++) {
NUIS_LOG(REC,
"|-> Par " << i << ". " << fNameList[i] << " " << fValueList[i]);
}
}
diff --git a/src/Reweight/MINERvAWeightCalcs.cxx b/src/Reweight/MINERvAWeightCalcs.cxx
index 0f555b4..68f7aa7 100644
--- a/src/Reweight/MINERvAWeightCalcs.cxx
+++ b/src/Reweight/MINERvAWeightCalcs.cxx
@@ -1,871 +1,871 @@
#ifdef __MINERVA_RW_ENABLED__
#ifdef __GENIE_ENABLED__
#include "MINERvAWeightCalcs.h"
#include "BaseFitEvt.h"
namespace nuisance {
namespace reweight {
//*******************************************************
MINERvAReWeight_QE::MINERvAReWeight_QE() {
//*******************************************************
fTwk_NormCCQE = 0.0;
fDef_NormCCQE = 1.0;
fCur_NormCCQE = fDef_NormCCQE;
}
//*******************************************************
MINERvAReWeight_QE::~MINERvAReWeight_QE(){};
//*******************************************************
//*******************************************************
double MINERvAReWeight_QE::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 QE this Calc doesn't handle it
if (!proc_info.IsQuasiElastic()) return 1.0;
// WEIGHT CALCULATIONS -------------
double w = 1.0;
// CCQE Dial
if (!proc_info.IsWeakCC()) w *= fCur_NormCCQE;
// Return Combined Weight
return w;
}
//*******************************************************
void MINERvAReWeight_QE::SetDialValue(std::string name, double val) {
//*******************************************************
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//*******************************************************
void MINERvAReWeight_QE::SetDialValue(int rwenum, double val) {
//*******************************************************
// Check Handled
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
if (!IsHandled(curenum)) return;
// Set Values
if (curenum == Reweight::kMINERvARW_NormCCQE) {
fTwk_NormCCQE = val;
fCur_NormCCQE = fDef_NormCCQE + fTwk_NormCCQE;
}
// Define Tweaked
fTweaked = ((fTwk_NormCCQE != 0.0));
}
//*******************************************************
bool MINERvAReWeight_QE::IsHandled(int rwenum) {
//*******************************************************
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::kMINERvARW_NormCCQE:
return true;
default:
return false;
}
}
//*******************************************************
MINERvAReWeight_MEC::MINERvAReWeight_MEC() {
//*******************************************************
fTwk_NormCCMEC = 0.0;
fDef_NormCCMEC = 1.0;
fCur_NormCCMEC = fDef_NormCCMEC;
}
//*******************************************************
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.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;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
if (!IsHandled(curenum))
return;
// Set Values
if (curenum == Reweight::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;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::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;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
if (!IsHandled(curenum)) return;
// Set Values
if (curenum == Reweight::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;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::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;
// - Syst : kMINERvA_RikRESRPA_ApplyRPA
// - Type : Binary
// - Limits : 0.0 (false) -> 1.0 (true)
// - Default : 0.0
fApplyDial_RESRPACorrection = false;
// - Syst : kMINERvA_RikRESRPA_LowQ2
// - Type : Absolute
// - Limits : 1.0 -> 1.0
// - Default : 0.0
// - Frac Error : 100%
fDefDial_RESRPALowQ2 = 0.0;
fCurDial_RESRPALowQ2 = fDefDial_RESRPALowQ2;
fErrDial_RESRPALowQ2 = 0.0;
// - Syst : kMINERvA_RikRESRPA_HighQ2
// - Type : Absolute
// - Limits : 1.0 -> 1.0
// - Default : 0.0
// - Frac Error : 100%
fDefDial_RESRPAHighQ2 = 0.0;
fCurDial_RESRPAHighQ2 = fDefDial_RESRPAHighQ2;
fErrDial_RESRPAHighQ2 = 1.0;
// Setup Calculators
fEventWeights = new double[5];
for (int 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() && !proc_info.IsResonant())
return 1.0;
// Extract Beam and Target PDG
GHepParticle *neutrino = ghep->Probe();
int bpdg = neutrino->Pdg();
GHepParticle *target = ghep->TargetNucleus();
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) {
NUIS_ABORT("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());
double Q2 = fabs((k1 - k2).Mag2());
// Quasielastic
if (proc_info.IsQuasiElastic()) {
// Now use q0-q3 and RPA Calculator to fill fWeights
rpacalc->getWeight(q0, q3, fEventWeights);
if (fApplyDial_RPACorrection) {
w *= fEventWeights[0]; // CV
}
// Syst Application : kMINERvA_RikRPA_LowQ2
if (fabs(fCurDial_RPALowQ2) > 0.0) {
double interpw = fEventWeights[0];
if (fCurDial_RPALowQ2 > 0.0 && Q2 < 2.0) {
interpw =
fEventWeights[0] - (fEventWeights[0] - fEventWeights[1]) *
fCurDial_RPALowQ2; // WLow+
} else if (fCurDial_RPALowQ2 < 0.0 && Q2 < 2.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
}
}
// Resonant Events
if (proc_info.IsResonant()) {
// Now use Q2 and RESRPA Calculator to fill fWeights
double CV = rpacalc->getWeight(Q2);
if (fApplyDial_RESRPACorrection) {
w *= CV; // fEventWeights[0]; // CVa
}
/*
// Syst Application : kMINERvA_RikRESRPA_LowQ2
if (fabs(fCurDial_RESRPAHighQ2) > 0.0) {
double interpw = fEventWeights[0];
if (fCurDial_RESRPAHighQ2 > 0.0) {
interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) *
fCurDial_RESRPAHighQ2; // WHigh+
} else if (fCurDial_RESRPAHighQ2 < 0.0) {
interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) *
fCurDial_RESRPAHighQ2; // WHigh-
}
w *= interpw / fEventWeights[0]; // Div by CV again
}
// Syst Application : kMINERvA_RikRESRPA_HighQ2
if (fabs(fCurDial_RESRPAHighQ2) > 0.0) {
double interpw = fEventWeights[0];
if (fCurDial_RESRPAHighQ2 > 0.0) {
interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) *
fCurDial_RESRPAHighQ2; // WHigh+
} else if (fCurDial_RESRPAHighQ2 < 0.0) {
interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) *
fCurDial_RESRPAHighQ2; // WHigh-
}
w *= interpw / fEventWeights[0]; // Div by CV again
}
*/
}
// LOG(FIT) << "RPA Weight = " << w << std::endl;
return w;
} // namespace reweight
//*******************************************************
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;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
// Check Handled
if (!IsHandled(curenum))
return;
if (curenum == Reweight::kMINERvARW_RikRPA_ApplyRPA)
fApplyDial_RPACorrection = (val > 0.5);
if (curenum == Reweight::kMINERvARW_RikRPA_LowQ2)
fCurDial_RPALowQ2 = val;
if (curenum == Reweight::kMINERvARW_RikRPA_HighQ2)
fCurDial_RPAHighQ2 = val;
if (curenum == Reweight::kMINERvARW_RikRESRPA_ApplyRPA)
fApplyDial_RESRPACorrection = (val > 0.5);
if (curenum == Reweight::kMINERvARW_RikRESRPA_LowQ2)
fCurDial_RESRPALowQ2 = val;
if (curenum == Reweight::kMINERvARW_RikRESRPA_HighQ2)
fCurDial_RESRPAHighQ2 = 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 ||
fApplyDial_RESRPACorrection ||
fabs(fCurDial_RESRPAHighQ2 - fDefDial_RESRPAHighQ2) > 0.0 ||
fabs(fCurDial_RESRPALowQ2 - fDefDial_RESRPALowQ2) > 0.0);
}
//*******************************************************
bool RikRPA::IsHandled(int rwenum) {
//*******************************************************
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::kMINERvARW_RikRESRPA_ApplyRPA:
return true;
case Reweight::kMINERvARW_RikRESRPA_LowQ2:
return true;
case Reweight::kMINERvARW_RikRESRPA_HighQ2:
return true;
case Reweight::kMINERvARW_RikRPA_ApplyRPA:
return true;
case Reweight::kMINERvARW_RikRPA_LowQ2:
return true;
case Reweight::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;
}
NUIS_LOG(FIT, "Loading RPA CALC : " << fidir);
TDirectory *olddir = gDirectory;
NUIS_LOG(FIT, "***********************************************");
NUIS_LOG(FIT, "Loading a new weightRPA calculator");
NUIS_LOG(FIT, "Authors: Rik Gran, Heidi Schellman");
NUIS_LOG(FIT, "Citation: arXiv:1705.02932 [hep-ex]");
NUIS_LOG(FIT, "***********************************************");
// Test the file exists
std::ifstream infile((rwdir + fidir).c_str());
if (!infile.good()) {
NUIS_ERR(FTL, "*** ERROR ***");
NUIS_ERR(FTL, "RikRPA file " << rwdir + fidir << " does not exist!");
NUIS_ERR(FTL,
"These can be found at https://nuisance.hepforge.org/files/RikRPA/");
NUIS_ERR(FTL,
"Please run: wget -r -nH --cut-dirs=2 -np -e robots=off -R "
"\"index.html*\" https://nuisance.hepforge.org/files/RikRPA/ -P "
<< rwdir);
NUIS_ERR(FTL, "And try again");
NUIS_ABORT("*************");
}
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 {
// NUIS_ERR(WRN, "Unknown beam and target combination for RPA Calcs! "
//<< bpdg << " " << tpdg);
}
return -1;
}
//*****************************************************************************
COHBrandon::COHBrandon() {
fApply_COHNorm = false;
fDef_COHNorm = 0.5;
fCur_COHNorm = fDef_COHNorm;
fTwk_COHNorm = 0.0;
fDef_COHCut = 0.450;
fCur_COHCut = fDef_COHCut;
fTwk_COHNorm = 0.0;
}
COHBrandon::~COHBrandon() {};
double COHBrandon::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 QE this Calc doesn't handle it
if (!proc_info.IsCoherent()) return 1.0;
// WEIGHT CALCULATIONS -------------
double w = 1.0;
double pionE = -999.9;
TObjArrayIter piter(ghep);
GHepParticle * p = 0;
//int ip = -1;
while ( (p = (GHepParticle *) piter.Next()) ) {
int pdgc = p->Pdg();
int ist = p->Status();
// if (ghep->Particle(ip)->FirstMother()==0) continue;
if (pdg::IsPseudoParticle(p->Pdg())) continue;
if (ist == kIStStableFinalState) {
if (pdgc == 211 || pdgc == -211 || pdgc == 111) {
pionE = p->Energy();
break;
}
}
}
// std::cout << "Coherent Pion Energy : " << pionE << std::endl;
// Apply Weight
if (fApply_COHNorm && pionE > 0.0 && pionE < fCur_COHCut) {
w *= fCur_COHNorm;
}
// Return Combined Weight
return w;
}
void COHBrandon::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void COHBrandon::SetDialValue(int rwenum, double val) {
// Check Handled
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
if (!IsHandled(curenum)) return;
// Set Values
if (curenum == Reweight::kMINERvARW_NormCOH) {
fTwk_COHNorm = val;
fCur_COHNorm = fDef_COHNorm * (1.0 + 0.1 * fTwk_COHNorm);
}
if (curenum == Reweight::kMINERvARW_CutCOH) {
fTwk_COHCut = val;
fCur_COHCut = fDef_COHCut * (1.0 + 0.1 * fTwk_COHCut);
}
if (curenum == Reweight::kMINERvARW_ApplyCOH) {
fApply_COHNorm = (val > 0.5);
}
// Define Tweaked
fTweaked = (fApply_COHNorm);
}
bool COHBrandon::IsHandled(int rwenum) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::kMINERvARW_NormCOH:
case Reweight::kMINERvARW_CutCOH:
case Reweight::kMINERvARW_ApplyCOH:
return true;
default:
return false;
}
}
//*****************************************************************************
WEnhancement::WEnhancement() {
fApply_Enhancement = false;
fDef_WNorm = 1.0;
fCur_WNorm = fDef_WNorm;
fTwk_WNorm = 0.0;
fDef_WMean = 0.95;
fCur_WMean = fDef_WMean;
fTwk_WMean = 0.0;
fDef_WSigma = 0.225;
fCur_WSigma = fDef_WSigma;
fTwk_WSigma = 0.0;
}
WEnhancement::~WEnhancement() {};
double WEnhancement::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 Kinematics & kine = interaction->Kine();
//const InitialState& init_state = interaction->InitState();
const ProcessInfo& proc_info = interaction->ProcInfo();
//const Target& tgt = init_state.Tgt();
// If the event is not QE this Calc doesn't handle it
if (!proc_info.IsWeakCC()) return 1.0;
if (!proc_info.IsResonant()) return 1.0;
if (!fApply_Enhancement) return 1.0;
// WEIGHT CALCULATIONS -------------
double w = 1.0;
bool isCC1pi0AtVertex = false;
// Get W
GHepParticle* neutrino = ghep->Probe();
GHepParticle* fsl = ghep->FinalStatePrimaryLepton();
const TLorentzVector& k1 = *(neutrino->P4());
const TLorentzVector& k2 = *(fsl->P4());
double E_mu = k2.E();
double p_mu = k2.Vect().Mag();
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = k1.Vect().Angle(k2.Vect());
// The factor of 1000 is necessary for downstream functions
const double m_p = PhysConst::mass_proton;
double E_nu = k1.E();
// double hadMass = kine.W (true);
double hadMass = sqrt(m_p * m_p + m_mu * m_mu - 2 * m_p * E_mu + \
2 * E_nu * (m_p - E_mu + p_mu * cos(th_nu_mu)));
// Determine if event is CC1pi0 at vertex
TObjArrayIter piter(ghep);
GHepParticle * p = 0;
int pi0 = 0;
int piother = 0;
while ( (p = (GHepParticle *) piter.Next()) ) {
int pdgc = p->Pdg();
int ist = p->Status();
if (pdg::IsPseudoParticle(p->Pdg())) continue;
if (ist == genie::kIStStableFinalState) {
if (pdgc == 111) {
pi0++;
} else if (pdgc == 211 || pdgc == -211) {
piother++;
}
}
}
// std::cout << "Npi0 = " << pi0 << std::endl;
isCC1pi0AtVertex = (pi0 == 1 && piother == 0);
// Apply Weight
if (fApply_Enhancement && isCC1pi0AtVertex) {
double enhancement = 1.0 + fCur_WNorm * (exp( -0.5 * pow((fCur_WMean - hadMass) / (fCur_WSigma), 2) ));
w *= enhancement;
}
// Return Combined Weight
if (isnan(w) || w < 0.0 || w > 400.0) {
w = 1.0;
}
return w;
}
void WEnhancement::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void WEnhancement::SetDialValue(int rwenum, double val) {
// Check Handled
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
if (!IsHandled(curenum)) return;
// Set Values
if (curenum == Reweight::kMINERvARW_ApplyWTune) {
fApply_Enhancement = (val > 0.5);
}
if (curenum == Reweight::kMINERvARW_NormWTune) {
fTwk_WNorm = val;
fCur_WNorm = fDef_WNorm * (1.0 + 0.1 * fTwk_WNorm);
}
if (curenum == Reweight::kMINERvARW_MeanWTune) {
fTwk_WMean = val;
fCur_WMean = fDef_WMean * (1.0 + 0.1 * fTwk_WMean);
}
if (curenum == Reweight::kMINERvARW_SigmaWTune) {
fTwk_WSigma = val;
fCur_WSigma = fDef_WSigma * (1.0 + 0.1 * fTwk_WSigma);
}
// Define Tweaked
fTweaked = (fApply_Enhancement);
}
bool WEnhancement::IsHandled(int rwenum) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::kMINERvARW_ApplyWTune:
case Reweight::kMINERvARW_NormWTune:
case Reweight::kMINERvARW_MeanWTune:
case Reweight::kMINERvARW_SigmaWTune:
return true;
default:
return false;
}
}
} // namespace reweight
} // namespace nuisance
#endif
#endif
diff --git a/src/Reweight/NEUTWeightEngine.cxx b/src/Reweight/NEUTWeightEngine.cxx
index 66aed33..c4542b8 100644
--- a/src/Reweight/NEUTWeightEngine.cxx
+++ b/src/Reweight/NEUTWeightEngine.cxx
@@ -1,208 +1,208 @@
#include "NEUTWeightEngine.h"
NEUTWeightEngine::NEUTWeightEngine(std::string name) {
-#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
+#if defined(__NEUT_ENABLED__) and defined(__USE_NEUT_REWEIGHT__)
#if defined(__NEUT_VERSION__) && (__NEUT_VERSION__ >= 541)
std::string neut_card = FitPar::Config().GetParS("NEUT_CARD");
if (!neut_card.size()) {
NUIS_ABORT(
"[ERROR]: When using NEUTReWeight must set NEUT_CARD config option.");
}
// No need to vomit the contents of the card file all over my screen
StopTalking();
neut::CommonBlockIFace::Initialize(neut_card);
StartTalking();
#endif
// Setup the NEUT Reweight engien
fCalcName = name;
NUIS_LOG(FIT, "Setting up NEUT RW : " << fCalcName);
// Create RW Engine suppressing cout
StopTalking();
fNeutRW = new neut::rew::NReWeight();
TDirectory *olddir = gDirectory;
// get list of vetoed calc engines (just for debug really)
std::string rw_engine_list =
FitPar::Config().GetParS("FitWeight_fNeutRW_veto");
bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos;
bool xsec_res = rw_engine_list.find("xsec_res") == std::string::npos;
// Activate each calc engine
if (xsec_ccqe)
fNeutRW->AdoptWghtCalc("xsec_ccqe", new neut::rew::NReWeightNuXSecCCQE);
if (xsec_res)
fNeutRW->AdoptWghtCalc("xsec_res", new neut::rew::NReWeightNuXSecRES);
// Dials removed in NEUT 5.4.1
#if __NEUT_VERSION__ < 541
bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos;
bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos;
bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos;
bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos;
bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos;
bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos;
bool nucl_casc = rw_engine_list.find("nucl_casc") == std::string::npos;
bool nucl_piless = rw_engine_list.find("nucl_piless") == std::string::npos;
if (nucl_casc)
fNeutRW->AdoptWghtCalc("nucl_casc", new neut::rew::NReWeightCasc);
if (xsec_coh)
fNeutRW->AdoptWghtCalc("xsec_coh", new neut::rew::NReWeightNuXSecCOH);
if (xsec_nc)
fNeutRW->AdoptWghtCalc("xsec_nc", new neut::rew::NReWeightNuXSecNC);
if (nucl_piless)
fNeutRW->AdoptWghtCalc("nucl_piless", new neut::rew::NReWeightNuclPiless);
if (xsec_ncres)
fNeutRW->AdoptWghtCalc("xsec_ncres", new neut::rew::NReWeightNuXSecNCRES);
if (xsec_ccres)
fNeutRW->AdoptWghtCalc("xsec_ccres", new neut::rew::NReWeightNuXSecCCRES);
if (xsec_ncel)
fNeutRW->AdoptWghtCalc("xsec_ncel", new neut::rew::NReWeightNuXSecNCEL);
if (xsec_dis)
fNeutRW->AdoptWghtCalc("xsec_dis", new neut::rew::NReWeightNuXSecDIS);
#endif
fNeutRW->Reconfigure();
olddir->cd();
// Set Abs Twk Config
fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
// allow cout again
StartTalking();
#else
NUIS_ABORT("NEUT RW NOT ENABLED!");
#endif
};
void NEUTWeightEngine::IncludeDial(std::string name, double startval) {
-#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
+#if defined(__NEUT_ENABLED__) and defined(__USE_NEUT_REWEIGHT__)
// Get First enum
int nuisenum = Reweight::ConvDial(name, kNEUT);
// 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 Syst
#if __NEUT_VERSION__ < 541
neut::rew::NSyst_t gensyst = NSyst::FromString(singlename);
#else
neut::rew::NSyst_t gensyst = neut::rew::NSyst::FromString(singlename);
#endif
// Fill Maps
int index = fValues.size();
fValues.push_back(0.0);
fNEUTSysts.push_back(gensyst);
// Initialize dial
NUIS_LOG(FIT, "Registering " << singlename << " dial.");
fNeutRW->Systematics().Init(fNEUTSysts[index]);
// If Absolute
if (fIsAbsTwk) {
#if __NEUT_VERSION__ < 541
NSystUncertainty::Instance()->SetUncertainty(fNEUTSysts[index], 1.0, 1.0);
#else
neut::rew::NSystUncertainty::Instance()->SetUncertainty(fNEUTSysts[index], 1.0, 1.0);
#endif
}
// Setup index
fEnumIndex[nuisenum].push_back(index);
fNameIndex[name].push_back(index);
}
// Set Value if given
if (startval != _UNDEF_DIAL_VALUE_) {
SetDialValue(nuisenum, startval);
}
#endif
}
void NEUTWeightEngine::SetDialValue(int nuisenum, double val) {
-#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
+#if defined(__NEUT_ENABLED__) and defined(__USE_NEUT_REWEIGHT__)
std::vector indices = fEnumIndex[nuisenum];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
NUIS_LOG(FIT, "Setting Dial Value for "
<< nuisenum << " " << i << " " << indices[i] << " "
<< fValues[indices[i]]
<< " Enum: " << fNEUTSysts[indices[i]]);
fNeutRW->Systematics().Set(fNEUTSysts[indices[i]], val);
}
#endif
}
void NEUTWeightEngine::SetDialValue(std::string name, double val) {
-#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
+#if defined(__NEUT_ENABLED__) and defined(__USE_NEUT_REWEIGHT__)
std::vector indices = fNameIndex[name];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
NUIS_LOG(FIT, "Setting Dial Value for "
<< name << " = " << i << " " << indices[i] << " "
<< fValues[indices[i]]
<< " Enum: " << fNEUTSysts[indices[i]]);
fNeutRW->Systematics().Set(fNEUTSysts[indices[i]], val);
}
#endif
}
void NEUTWeightEngine::Reconfigure(bool silent) {
-#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
+#if defined(__NEUT_ENABLED__) and defined(__USE_NEUT_REWEIGHT__)
// Hush now...
if (silent)
StopTalking();
// Reconf
fNeutRW->Reconfigure();
// if (LOG_LEVEL(DEB)){
fNeutRW->Print();
// }
// Shout again
if (silent)
StartTalking();
#endif
}
double NEUTWeightEngine::CalcWeight(BaseFitEvt *evt) {
double rw_weight = 1.0;
-#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
+#if defined(__NEUT_ENABLED__) and defined(__USE_NEUT_REWEIGHT__)
// Skip Non NEUT
if (evt->fType != kNEUT)
return 1.0;
// Hush now
StopTalking();
// Fill NEUT Common blocks
NEUTUtils::FillNeutCommons(evt->fNeutVect);
// Call Weight calculation
rw_weight = fNeutRW->CalcWeight();
// Speak Now
StartTalking();
#endif
if (!std::isnormal(rw_weight)){
NUIS_ERR(WRN, "NEUT returned weight: " << rw_weight);
rw_weight = 0;
}
// Return rw_weight
return rw_weight;
}
diff --git a/src/Reweight/NEUTWeightEngine.h b/src/Reweight/NEUTWeightEngine.h
index 17ff6ca..009cd9f 100644
--- a/src/Reweight/NEUTWeightEngine.h
+++ b/src/Reweight/NEUTWeightEngine.h
@@ -1,61 +1,57 @@
#ifndef WEIGHT_ENGINE_NEUT_H
#define WEIGHT_ENGINE_NEUT_H
#include "FitLogger.h"
-#ifdef __NEUT_ENABLED__
-#ifndef __NO_REWEIGHT__
+#if defined(__NEUT_ENABLED__) and defined(__USE_NEUT_REWEIGHT__)
#include "NEUTInputHandler.h"
#include "NReWeight.h"
#include "NReWeightNuXSecCCQE.h"
#include "NReWeightNuXSecRES.h"
// Dials removed in NEUT 5.4.1
#if __NEUT_VERSION__ < 541
#include "NReWeightCasc.h"
#include "NReWeightNuclPiless.h"
#include "NReWeightNuXSecNCRES.h"
#include "NReWeightNuXSecCCRES.h"
#include "NReWeightNuXSecNC.h"
#include "NReWeightNuXSecCOH.h"
#include "NReWeightNuXSecNCEL.h"
#include "NReWeightNuXSecDIS.h"
#endif
#if __NEUT_VERSION__ >= 541
#include "CommonBlockIFace.h"
#endif
#include "NSyst.h"
#include "NSystUncertainty.h"
#include "neutpart.h"
#include "neutvect.h"
#endif
-#endif
#include "FitWeight.h"
#include "GeneratorUtils.h"
#include "WeightEngineBase.h"
class NEUTWeightEngine : public WeightEngineBase {
public:
NEUTWeightEngine(std::string name);
~NEUTWeightEngine(){};
void IncludeDial(std::string name, double startval);
void SetDialValue(std::string name, double val);
void SetDialValue(int nuisenum, double val);
void Reconfigure(bool silent = false);
double CalcWeight(BaseFitEvt *evt);
inline bool NeedsEventReWeight() { return true; };
-#ifdef __NEUT_ENABLED__
-#ifndef __NO_REWEIGHT__
+#if defined(__NEUT_ENABLED__) and defined(__USE_NEUT_REWEIGHT__)
std::vector fNEUTSysts;
neut::rew::NReWeight *fNeutRW;
#endif
-#endif
};
#endif
diff --git a/src/Reweight/NOvARwgtEngine.cxx b/src/Reweight/NOvARwgtEngine.cxx
index 7a60320..e184834 100644
--- a/src/Reweight/NOvARwgtEngine.cxx
+++ b/src/Reweight/NOvARwgtEngine.cxx
@@ -1,317 +1,317 @@
#include "NOvARwgtEngine.h"
#include "NOvARwgt/interfaces/GenieInterface.h"
#include "NOvARwgt/rwgt/tunes/Tunes2017.h"
#include "NOvARwgt/rwgt/tunes/Tunes2018.h"
#include "NOvARwgt/rwgt/tunes/Tunes2020.h"
#include "NOvARwgt/rwgt/tunes/TunesSA.h"
#include "NOvARwgt/rwgt/EventRecord.h"
#include "NOvARwgt/rwgt/ISystKnob.h"
#include "NOvARwgt/rwgt/Tune.h"
// GENIEv3
#include "Framework/Utils/AppInit.h"
#include "Framework/Utils/RunOpt.h"
static size_t const kCVTune2017 = 100;
static size_t const kCVTune2018 = 200;
static size_t const kCVTune2020 = 300;
static size_t const kCVTuneSA = 400;
static size_t const kNoSuchKnob = std::numeric_limits::max();
novarwgt::Tune const *TuneFactory(size_t e) {
switch (e) {
case kCVTune2017: {
return &novarwgt::kCVTune2017;
}
case kCVTune2018: {
return &novarwgt::kCVTune2018;
}
case kCVTune2020: {
return &novarwgt::kCVTune2020;
}
case kCVTuneSA: {
return &novarwgt::kCVTuneSA;
}
default: {
return NULL;
}
}
}
void NOvARwgtEngine::InitializeKnobs() {
fTunes.push_back(TuneFactory(kCVTune2017));
fTuneEnums[kCVTune2017] = 0;
fTunes.push_back(TuneFactory(kCVTune2018));
fTuneEnums[kCVTune2018] = 1;
fTunes.push_back(TuneFactory(kCVTune2020));
fTuneEnums[kCVTune2020] = 2;
fTunes.push_back(TuneFactory(kCVTuneSA));
fTuneEnums[kCVTuneSA] = 3;
fTuneInUse = {false, false, false, false};
fTuneValues = {0, 0, 0, 0};
size_t ctr = 0;
NUIS_LOG(FIT, "NOvARwgt kCVTune2017 sub-knobs:");
size_t tune_ctr = 0;
for (auto &k : novarwgt::kCVTune2017.SystKnobs()) {
NUIS_LOG(FIT, "\t" << k.first);
fKnobs.push_back(k.second);
fKnobTunes.push_back(&novarwgt::kCVTune2017);
fKnobTuneidx.push_back(0);
fKnobEnums[kCVTune2017 + 1 + tune_ctr++] = ctr++;
fKnobInUse.push_back(false);
fKnobValues.push_back(0);
}
tune_ctr = 0;
NUIS_LOG(FIT, "NOvARwgt kCVTune2018 sub-knobs:");
for (auto &k : novarwgt::kCVTune2018.SystKnobs()) {
NUIS_LOG(FIT, "\t" << k.first);
fKnobs.push_back(k.second);
fKnobTunes.push_back(&novarwgt::kCVTune2018);
fKnobTuneidx.push_back(1);
fKnobEnums[kCVTune2018 + 1 + tune_ctr++] = ctr++;
fKnobInUse.push_back(false);
fKnobValues.push_back(0);
}
tune_ctr = 0;
NUIS_LOG(FIT, "NOvARwgt kCVTune2020 sub-knobs:");
for (auto &k : novarwgt::kCVTune2020.SystKnobs()) {
NUIS_LOG(FIT, "\t" << k.first);
fKnobs.push_back(k.second);
fKnobTunes.push_back(&novarwgt::kCVTune2020);
fKnobTuneidx.push_back(2);
fKnobEnums[kCVTune2020 + 1 + tune_ctr++] = ctr++;
fKnobInUse.push_back(false);
fKnobValues.push_back(0);
}
tune_ctr = 0;
NUIS_LOG(FIT, "NOvARwgt kCVTuneSA sub-knobs:");
for (auto &k : novarwgt::kCVTuneSA.SystKnobs()) {
NUIS_LOG(FIT, "\t" << k.first);
fKnobs.push_back(k.second);
fKnobTunes.push_back(&novarwgt::kCVTuneSA);
fKnobTuneidx.push_back(3);
fKnobEnums[kCVTuneSA + 1 + tune_ctr++] = ctr++;
fKnobInUse.push_back(false);
fKnobValues.push_back(0);
}
}
void NOvARwgtEngine::InitializeGENIE() {
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.");
}
grunopt->SetTuneName(Config::GetParS("GENIETune"));
grunopt->BuildTune();
std::string genv =
std::string(getenv("GENIE")) + "/config/Messenger_laconic.xml";
genie::utils::app_init::MesgThresholds(genv);
}
size_t NOvARwgtEngine::GetWeightGeneratorIndex(std::string const &strname) {
size_t upos = strname.find_first_of("_");
if (strname.find("CVTune2017") == 0) {
if (upos == std::string::npos) {
return kCVTune2017;
}
std::string knobname = strname.substr(upos + 1);
if (novarwgt::kCVTune2017.SystKnobs().count(knobname)) {
auto loc = novarwgt::kCVTune2017.SystKnobs().find(knobname);
return kCVTune2017 + 1 +
std::distance(novarwgt::kCVTune2017.SystKnobs().begin(), loc);
}
} else if (strname.find("CVTune2018") == 0) {
if (upos == std::string::npos) {
return kCVTune2018;
}
std::string knobname = strname.substr(upos + 1);
if (novarwgt::kCVTune2018.SystKnobs().count(knobname)) {
auto loc = novarwgt::kCVTune2018.SystKnobs().find(knobname);
return kCVTune2018 + 1 +
std::distance(novarwgt::kCVTune2018.SystKnobs().begin(), loc);
}
} else if (strname.find("CVTune2020") == 0) {
if (upos == std::string::npos) {
return kCVTune2020;
}
std::string knobname = strname.substr(upos + 1);
if (novarwgt::kCVTune2020.SystKnobs().count(knobname)) {
auto loc = novarwgt::kCVTune2020.SystKnobs().find(knobname);
return kCVTune2020 + 1 +
std::distance(novarwgt::kCVTune2020.SystKnobs().begin(), loc);
}
} else if (strname.find("CVTuneSA") == 0) {
if (upos == std::string::npos) {
return kCVTuneSA;
}
std::string knobname = strname.substr(upos + 1);
if (novarwgt::kCVTuneSA.SystKnobs().count(knobname)) {
auto loc = novarwgt::kCVTuneSA.SystKnobs().find(knobname);
return kCVTuneSA + 1 +
std::distance(novarwgt::kCVTuneSA.SystKnobs().begin(), loc);
}
}
return kNoSuchKnob;
}
void NOvARwgtEngine::IncludeDial(std::string name, double startval) {
size_t we_indx = GetWeightGeneratorIndex(name);
if (we_indx == kNoSuchKnob) {
NUIS_ABORT("[ERROR]: Invalid NOvARwgt Engine name: " << name);
}
- bool IsTune = !(we_indx % 100);
+ bool IsTune = !(we_indx % (NUIS_DIAL_OFFSET /10));
if (IsTune) {
auto tune_idx = fTuneEnums[we_indx];
fTuneValues[tune_idx] = startval;
fTuneInUse[tune_idx] = true;
} else {
auto knob_idx = fKnobEnums[we_indx];
fKnobValues[knob_idx] = startval;
fKnobInUse[knob_idx] = true;
}
};
void NOvARwgtEngine::SetDialValue(int nuisenum, double val) {
- size_t we_indx = (nuisenum % 1000);
- bool IsTune = !(we_indx % 100);
+ size_t we_indx = (nuisenum % NUIS_DIAL_OFFSET);
+ bool IsTune = !(we_indx % (NUIS_DIAL_OFFSET /10));
if (IsTune) {
auto tune_idx = fTuneEnums[we_indx];
if (!fTuneInUse[tune_idx]) {
NUIS_ABORT("[ERROR]: SetDialValue for NOvARwgt dial: "
<< we_indx << " but that tune hasn't been included yet.");
}
fTuneValues[tune_idx] = val;
} else {
auto knob_idx = fKnobEnums[we_indx];
if (!fKnobInUse[knob_idx]) {
NUIS_ABORT("[ERROR]: SetDialValue for NOvARwgt dial: "
<< we_indx << " but that Knob hasn't been included yet.");
}
fKnobValues[knob_idx] = val;
}
}
void NOvARwgtEngine::SetDialValue(std::string name, double val) {
SetDialValue(GetWeightGeneratorIndex(name), val);
}
bool NOvARwgtEngine::IsDialIncluded(std::string name) {
return IsDialIncluded(GetWeightGeneratorIndex(name));
}
bool NOvARwgtEngine::IsDialIncluded(int nuisenum) {
- size_t we_indx = (nuisenum % 1000);
- bool IsTune = !(we_indx % 100);
+ size_t we_indx = (nuisenum % NUIS_DIAL_OFFSET);
+ bool IsTune = !(we_indx % (NUIS_DIAL_OFFSET /10));
if (IsTune) {
auto tune_idx = fTuneEnums[we_indx];
return fTuneInUse[tune_idx];
} else {
auto knob_idx = fKnobEnums[we_indx];
return fKnobInUse[knob_idx];
}
}
double NOvARwgtEngine::GetDialValue(std::string name) {
return GetDialValue(GetWeightGeneratorIndex(name));
}
double NOvARwgtEngine::GetDialValue(int nuisenum) {
- size_t we_indx = (nuisenum % 1000);
+ size_t we_indx = (nuisenum % NUIS_DIAL_OFFSET);
if (we_indx == kNoSuchKnob) {
NUIS_ABORT("[ERROR]: Invalid NOvARwgt Engine enum: " << nuisenum);
}
- bool IsTune = !(we_indx % 100);
+ bool IsTune = !(we_indx % (NUIS_DIAL_OFFSET /10));
if (IsTune) {
auto tune_idx = fTuneEnums[we_indx];
if (!fTuneInUse[tune_idx]) {
NUIS_ABORT("[ERROR]: GetDialValue for NOvARwgt dial: "
<< we_indx << " but that tune hasn't been included yet.");
}
return fTuneValues[tune_idx];
} else {
auto knob_idx = fKnobEnums[we_indx];
if (!fKnobInUse[knob_idx]) {
NUIS_ABORT("[ERROR]: GetDialValue for NOvARwgt dial: "
<< we_indx << " but that Knob hasn't been included yet.");
}
return fKnobValues[knob_idx];
}
}
double NOvARwgtEngine::CalcWeight(BaseFitEvt *evt) {
double rw_weight = 1.0;
// 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));
}
novarwgt::EventRecord rcd =
novarwgt::ConvertGenieEvent(evt->genie_event->event);
for (size_t k_it = 0; k_it < fKnobs.size(); ++k_it) {
if (!fKnobInUse[k_it]) {
continue;
}
double wght = fKnobTunes[k_it]->EventSystKnobWeight(
fKnobs[k_it]->Name(), fKnobValues[k_it], rcd, {},
fTuneInUse[fKnobTuneidx[k_it]] && fTuneValues[fKnobTuneidx[k_it]]);
// // have to divide out the CV weight for this, ugly hack because the last
// // parameter doesn't do what I want
// if (fTuneInUse[fKnobTuneidx[k_it]] && fTuneValues[fKnobTuneidx[k_it]]) {
// wght /= fKnobTunes[k_it]->EventSystKnobWeight(fKnobs[k_it]->Name(), 0,
// rcd, {}, false);
// }
rw_weight *= wght;
}
for (size_t k_it = 0; k_it < fTunes.size(); ++k_it) {
if (!fTuneInUse[k_it]) {
continue;
}
if (!fTuneValues[k_it]) {
continue;
}
double wght = fTunes[k_it]->EventWeight(rcd);
rw_weight *= wght;
}
return rw_weight;
}
diff --git a/src/Reweight/NUISANCEWeightCalcs.cxx b/src/Reweight/NUISANCEWeightCalcs.cxx
index 426f582..7f9d8d2 100644
--- a/src/Reweight/NUISANCEWeightCalcs.cxx
+++ b/src/Reweight/NUISANCEWeightCalcs.cxx
@@ -1,792 +1,794 @@
#include "NUISANCEWeightCalcs.h"
#include "FitEvent.h"
#include "GeneralUtils.h"
#include "NUISANCESyst.h"
#include "WeightUtils.h"
+#include "WeightEngineBase.h"
+
using namespace Reweight;
ModeNormCalc::ModeNormCalc() { fNormRES = 1.0; }
double ModeNormCalc::CalcWeight(BaseFitEvt *evt) {
int mode = abs(evt->Mode);
double w = 1.0;
if (mode == 11 or mode == 12 or mode == 13) {
w *= fNormRES;
}
return w;
}
void ModeNormCalc::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void ModeNormCalc::SetDialValue(int rwenum, double val) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
// Check Handled
if (!IsHandled(curenum))
return;
if (curenum == kModeNorm_NormRES)
fNormRES = val;
}
bool ModeNormCalc::IsHandled(int rwenum) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case kModeNorm_NormRES:
return true;
default:
return false;
}
}
//*****************************************************************************
MINOSRPA::MINOSRPA() {
fApply_MINOSRPA = false;
fDef_MINOSRPA_A = 1.010;
fTwk_MINOSRPA_A = 0.000;
fDef_MINOSRPA_B = 0.156;
fTwk_MINOSRPA_B = 0.000;
}
//
double MINOSRPA::CalcWeight(BaseFitEvt* evt) {
if (!fTweaked || !fApply_MINOSRPA) return 1.0;
double w = 1.0;
// If GENIE is enabled, use old code
#ifdef __GENIE_ENABLED__
// 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 not on nucleus, not resonant, or NC
if (!tgt.IsNucleus() || !proc_info.IsResonant() || proc_info.IsWeakNC()) 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();
// 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());
double Q2 = fabs((k1 - k2).Mag2());
w *= GetRPAWeight(Q2);
#else
// Get the Q2 from NUISANCE if not GENIE
FitEvent *fevt = static_cast(evt);
// Check the event is resonant
if (!fevt->IsResonant() || fevt->IsNC()) return 1.0;
int targeta = fevt->GetTargetA();
int targetz = fevt->GetTargetZ();
// Apply only to nuclear targets, ignore free protons
if (targeta == 1 || targetz == 1) return 1.0;
// Q2 in GeV2
double Q2 = fevt->GetQ2();
w *= GetRPAWeight(Q2);
#endif
return w;
}
// Do the actual weight calculation
double MINOSRPA::GetRPAWeight(double Q2) {
if (Q2 > 0.7) return 1.0;
double w = fCur_MINOSRPA_A / (1.0 + TMath::Exp(1.0 - TMath::Sqrt(Q2) / fCur_MINOSRPA_B));
return w;
}
bool MINOSRPA::IsHandled(int rwenum) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::kMINERvARW_MINOSRPA_Apply:
case Reweight::kMINERvARW_MINOSRPA_A:
case Reweight::kMINERvARW_MINOSRPA_B:
return true;
default:
return false;
}
}
//
void MINOSRPA::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//
void MINOSRPA::SetDialValue(int rwenum, double val) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
// Check Handled
if (!IsHandled(curenum)) return;
if (curenum == Reweight::kMINERvARW_MINOSRPA_Apply) fApply_MINOSRPA = (val > 0.5);
if (curenum == Reweight::kMINERvARW_MINOSRPA_A) fTwk_MINOSRPA_A = val;
if (curenum == Reweight::kMINERvARW_MINOSRPA_B) fTwk_MINOSRPA_B = val;
// Check for changes
fTweaked = (fApply_MINOSRPA ||
fabs(fTwk_MINOSRPA_A) > 0.0 ||
fabs(fTwk_MINOSRPA_B) > 0.0);
// Update Values
fCur_MINOSRPA_A = fDef_MINOSRPA_A * (1.0 + 0.1 * fTwk_MINOSRPA_A);
fCur_MINOSRPA_B = fDef_MINOSRPA_B * (1.0 + 0.1 * fTwk_MINOSRPA_B);
}
//
//
//*****************************************************************************
LagrangeRPA::LagrangeRPA() {
fApplyRPA = false;
/*
fI1_Def = 4.18962e-01;
fI1 = fI1_Def;
fI2_Def = 7.39927e-01;
fI2 = fI2_Def;
*/
// Table VIII https://arxiv.org/pdf/1903.01558.pdf
fR1_Def = 0.37;
fR1 = fR1_Def;
fR2_Def = 0.60;
fR2 = fR2_Def;
}
//
double LagrangeRPA::CalcWeight(BaseFitEvt* evt) {
if (!fTweaked || !fApplyRPA) return 1.0;
double w = 1.0;
// If GENIE is enabled, use old code
#ifdef __GENIE_ENABLED__
// 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 not on nucleus, not resonant, or NC
if (!tgt.IsNucleus() || !proc_info.IsResonant() || proc_info.IsWeakNC()) 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();
// 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());
double Q2 = fabs((k1 - k2).Mag2());
w *= GetRPAWeight(Q2);
#else
// Get the Q2 from NUISANCE if not GENIE
FitEvent *fevt = static_cast(evt);
// Check the event is resonant
if (!fevt->IsResonant() || fevt->IsNC()) return 1.0;
int targeta = fevt->GetTargetA();
int targetz = fevt->GetTargetZ();
// Apply only to nuclear targets, ignore free protons
if (targeta == 1 || targetz == 1) return 1.0;
// Q2 in GeV2
double Q2 = fevt->GetQ2();
w *= GetRPAWeight(Q2);
#endif
return w;
}
//
//
double LagrangeRPA::GetRPAWeight(double Q2) {
//std::cout << "Getting RPA Weight : " << Q2 << std::endl;
if (Q2 > 0.7) return 1.0;
// Keep original Lagrange RPA for documentation
/*
double x1 = 0.00;
double x2 = 0.30;
double x3 = 0.70;
double y1 = 0.00;
double y2 = fI2;
double y3 = 1.00;
double xv = Q2;
// Automatically 0 because y1 choice
double W1 = y1 * (xv-x2)*(xv-x3)/((x1-x2)*(x1-x3));
double W2 = y2 * (xv-x1)*(xv-x3)/((x2-x1)*(x2-x3));
double W3 = y3 * (xv-x1)*(xv-x2)/((x3-x1)*(x3-x2));
double P = W1 + W2 + W3;
double A1 = (1.0 - sqrt(1.0 - fI1));
double R = P * (1.0 - A1) + A1;
return 1.0 - (1.0-R)*(1.0-R);
*/
// Equation 7 https://arxiv.org/pdf/1903.01558.pdf
const double x1 = 0.00;
const double x2 = 0.35;
const double x3 = 0.70;
// Equation 6 https://arxiv.org/pdf/1903.01558.pdf
double RQ2 = fR2 *( (Q2-x1)*(Q2-x3)/((x2-x1)*(x2-x3)) )
+ (Q2-x1)*(Q2-x2)/((x3-x1)*(x3-x2));
double weight = 1-(1-fR1)*(1-RQ2)*(1-RQ2);
// Check range of R1 and R2
// Commented out because this is implementation dependent: user may want strange double peaks
/*
if (fR1 > 1) return 1;
if (fR2 > 1 || fR2 < 0.5) return 1;
*/
return weight;
}
//
bool LagrangeRPA::IsHandled(int rwenum) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::kMINERvARW_LagrangeRPA_Apply:
case Reweight::kMINERvARW_LagrangeRPA_R1:
case Reweight::kMINERvARW_LagrangeRPA_R2:
return true;
default:
return false;
}
}
//
void LagrangeRPA::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//
void LagrangeRPA::SetDialValue(int rwenum, double val) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
// Check Handled
if (!IsHandled(curenum)) return;
if (curenum == Reweight::kMINERvARW_LagrangeRPA_Apply) fApplyRPA = (val > 0.5);
if (curenum == Reweight::kMINERvARW_LagrangeRPA_R1) fR1 = val;
if (curenum == Reweight::kMINERvARW_LagrangeRPA_R2) fR2 = val;
// Check for changes
fTweaked = (fApplyRPA);
}
//
BeRPACalc::BeRPACalc()
: fBeRPA_A(0.59), fBeRPA_B(1.05), fBeRPA_D(1.13), fBeRPA_E(0.88),
fBeRPA_U(1.2), nParams(0) {
// A = 0.59 +/- 20%
// B = 1.05 +/- 20%
// D = 1.13 +/- 15%
// E = 0.88 +/- 40%
// U = 1.2
}
double BeRPACalc::CalcWeight(BaseFitEvt *evt) {
int mode = abs(evt->Mode);
double w = 1.0;
if (nParams == 0) {
return w;
}
// Get Q2
// Get final state lepton
if (mode == 1) {
FitEvent *fevt = static_cast(evt);
double Q2 = fevt->GetQ2();
// Only CCQE events
w *= calcRPA(Q2, fBeRPA_A, fBeRPA_B, fBeRPA_D, fBeRPA_E, fBeRPA_U);
}
return w;
}
void BeRPACalc::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void BeRPACalc::SetDialValue(int rwenum, double val) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
// Check Handled
if (!IsHandled(curenum))
return;
// Need 4 or 5 reconfigures
if (curenum == kBeRPA_A)
fBeRPA_A = val;
else if (curenum == kBeRPA_B)
fBeRPA_B = val;
else if (curenum == kBeRPA_D)
fBeRPA_D = val;
else if (curenum == kBeRPA_E)
fBeRPA_E = val;
else if (curenum == kBeRPA_U)
fBeRPA_U = val;
nParams++;
}
bool BeRPACalc::IsHandled(int rwenum) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case kBeRPA_A:
case kBeRPA_B:
case kBeRPA_D:
case kBeRPA_E:
case kBeRPA_U:
return true;
default:
return false;
}
}
SBLOscWeightCalc::SBLOscWeightCalc() {
fDistance = 0.0;
fMassSplitting = 0.0;
fSin2Theta = 0.0;
}
double SBLOscWeightCalc::CalcWeight(BaseFitEvt *evt) {
FitEvent *fevt = static_cast(evt);
FitParticle *pnu = fevt->PartInfo(0);
double E = pnu->fP.E() / 1.E3;
// Extract energy
return GetSBLOscWeight(E);
}
void SBLOscWeightCalc::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void SBLOscWeightCalc::SetDialValue(int rwenum, double val) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
if (!IsHandled(curenum))
return;
if (curenum == kSBLOsc_Distance)
fDistance = val;
if (curenum == kSBLOsc_MassSplitting)
fMassSplitting = val;
if (curenum == kSBLOsc_Sin2Theta)
fSin2Theta = val;
}
bool SBLOscWeightCalc::IsHandled(int rwenum) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case kSBLOsc_Distance:
return true;
case kSBLOsc_MassSplitting:
return true;
case kSBLOsc_Sin2Theta:
return true;
default:
return false;
}
}
double SBLOscWeightCalc::GetSBLOscWeight(double E) {
if (E <= 0.0)
return 1.0 - 0.5 * fSin2Theta;
return 1.0 - fSin2Theta * pow(sin(1.267 * fMassSplitting * fDistance / E), 2);
}
GaussianModeCorr::GaussianModeCorr() {
// Apply the tilt-shift Gauss by Patrick
// Alternatively set in config
fMethod = true;
// Init
fApply_CCQE = false;
fGausVal_CCQE[kPosNorm] = 0.0;
fGausVal_CCQE[kPosTilt] = 0.0;
fGausVal_CCQE[kPosPq0] = 1.0;
fGausVal_CCQE[kPosWq0] = 1.0;
fGausVal_CCQE[kPosPq3] = 1.0;
fGausVal_CCQE[kPosWq3] = 1.0;
fApply_2p2h = false;
fGausVal_2p2h[kPosNorm] = 0.0;
fGausVal_2p2h[kPosTilt] = 0.0;
fGausVal_2p2h[kPosPq0] = 1.0;
fGausVal_2p2h[kPosWq0] = 1.0;
fGausVal_2p2h[kPosPq3] = 1.0;
fGausVal_2p2h[kPosWq3] = 1.0;
fApply_2p2h_PPandNN = false;
fGausVal_2p2h_PPandNN[kPosNorm] = 0.0;
fGausVal_2p2h_PPandNN[kPosTilt] = 0.0;
fGausVal_2p2h_PPandNN[kPosPq0] = 1.0;
fGausVal_2p2h_PPandNN[kPosWq0] = 1.0;
fGausVal_2p2h_PPandNN[kPosPq3] = 1.0;
fGausVal_2p2h_PPandNN[kPosWq3] = 1.0;
fApply_2p2h_NP = false;
fGausVal_2p2h_NP[kPosNorm] = 0.0;
fGausVal_2p2h_NP[kPosTilt] = 0.0;
fGausVal_2p2h_NP[kPosPq0] = 1.0;
fGausVal_2p2h_NP[kPosWq0] = 1.0;
fGausVal_2p2h_NP[kPosPq3] = 1.0;
fGausVal_2p2h_NP[kPosWq3] = 1.0;
fApply_CC1pi = false;
fGausVal_CC1pi[kPosNorm] = 0.0;
fGausVal_CC1pi[kPosTilt] = 0.0;
fGausVal_CC1pi[kPosPq0] = 1.0;
fGausVal_CC1pi[kPosWq0] = 1.0;
fGausVal_CC1pi[kPosPq3] = 1.0;
fGausVal_CC1pi[kPosWq3] = 1.0;
fAllowSuppression = false;
fDebugStatements = FitPar::Config().GetParB("GaussianModeCorr_DEBUG");
// fDebugStatements = true;
}
double GaussianModeCorr::CalcWeight(BaseFitEvt *evt) {
FitEvent *fevt = static_cast(evt);
double rw_weight = 1.0;
// Get Neutrino
if (!fevt->Npart()) {
NUIS_ABORT("NO particles found in stack!");
}
FitParticle *pnu = fevt->GetNeutrinoIn();
if (!pnu) {
NUIS_ABORT("NO Starting particle found in stack!");
}
FitParticle *plep = fevt->GetLeptonOut();
if (!plep) return 1.0;
TLorentzVector q = pnu->fP - plep->fP;
// Extra q0,q3
double q0 = fabs(q.E()) / 1.E3;
double q3 = fabs(q.Vect().Mag()) / 1.E3;
int initialstate = -1; // Undef
if (abs(fevt->Mode) == 2) {
int npr = 0;
int nne = 0;
for (UInt_t j = 0; j < fevt->Npart(); j++) {
if ((fevt->PartInfo(j))->fIsAlive)
continue;
if (fevt->PartInfo(j)->fPID == 2212)
npr++;
else if (fevt->PartInfo(j)->fPID == 2112)
nne++;
}
if (fevt->Mode == 2 && npr == 1 && nne == 1) {
initialstate = 2;
} else if (fevt->Mode == 2 &&
((npr == 0 && nne == 2) || (npr == 2 && nne == 0))) {
initialstate = 1;
}
}
// Apply weighting
if (fApply_CCQE && abs(fevt->Mode) == 1) {
if (fDebugStatements)
std::cout << "Getting CCQE Weight" << std::endl;
double g = GetGausWeight(q0, q3, fGausVal_CCQE);
if (g < 1.0)
g = 1.0;
rw_weight *= g;
}
if (fApply_2p2h && abs(fevt->Mode) == 2) {
if (fDebugStatements) std::cout << "Getting 2p2h Weight" << std::endl;
if (fDebugStatements) std::cout << "Got q0 q3 = " << q0 << " " << q3 << " mode = " << fevt->Mode << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h);
if (fDebugStatements) std::cout << "Returning Weight " << rw_weight << std::endl;
}
if (fApply_2p2h_PPandNN && abs(fevt->Mode) == 2 && initialstate == 1) {
if (fDebugStatements) std::cout << "Getting 2p2h PPandNN Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_PPandNN);
}
if (fApply_2p2h_NP && abs(fevt->Mode) == 2 && initialstate == 2) {
if (fDebugStatements) std::cout << "Getting 2p2h NP Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_NP);
}
if (fApply_CC1pi && abs(fevt->Mode) >= 11 && abs(fevt->Mode) <= 13) {
if (fDebugStatements)
std::cout << "Getting CC1pi Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_CC1pi);
}
return rw_weight;
}
void GaussianModeCorr::SetMethod(bool method) {
fMethod = method;
if (fMethod == true) {
NUIS_LOG(FIT,
" Using tilt-shift Gaussian parameters for Gaussian enhancement...");
} else {
NUIS_LOG(FIT, " Using Normal Gaussian parameters for Gaussian enhancement...");
}
};
double GaussianModeCorr::GetGausWeight(double q0, double q3, double vals[]) {
// The weight
double w = 1.0;
// Use tilt-shift method by Patrick
if (fMethod) {
if (fDebugStatements) {
std::cout << "Using Patrick gaussian" << std::endl;
}
// // CCQE Without Suppression
// double Norm = 4.82788679036;
// double Tilt = 2.3501416116;
// double Pq0 = 0.363964889702;
// double Wq0 = 0.133976806938;
// double Pq3 = 0.431769740224;
// double Wq3 = 0.207666663434;
// // Also add for CCQE at the end
// return (w > 1.0) ? w : 1.0;
// // 2p2h with suppression
// double Norm = 15.967;
// double Tilt = -0.455655;
// double Pq0 = 0.214598;
// double Wq0 = 0.0291061;
// double Pq3 = 0.480194;
// double Wq3 = 0.134588;
double Norm = vals[kPosNorm];
double Tilt = vals[kPosTilt];
double Pq0 = vals[kPosPq0];
double Wq0 = vals[kPosWq0];
double Pq3 = vals[kPosPq3];
double Wq3 = vals[kPosWq3];
double a = cos(Tilt) * cos(Tilt) / (2 * Wq0 * Wq0);
a += sin(Tilt) * sin(Tilt) / (2 * Wq3 * Wq3);
double b = -sin(2 * Tilt) / (4 * Wq0 * Wq0);
b += sin(2 * Tilt) / (4 * Wq3 * Wq3);
double c = sin(Tilt) * sin(Tilt) / (2 * Wq0 * Wq0);
c += cos(Tilt) * cos(Tilt) / (2 * Wq3 * Wq3);
w = Norm;
w *= exp(-a * (q0 - Pq0) * (q0 - Pq0));
w *= exp(+2.0 * b * (q0 - Pq0) * (q3 - Pq3));
w *= exp(-c * (q3 - Pq3) * (q3 - Pq3));
if (fDebugStatements) {
std::cout << "Applied Tilt " << Tilt << " " << cos(Tilt) << " "
<< sin(Tilt) << std::endl;
std::cout << "abc = " << a << " " << b << " " << c << std::endl;
std::cout << "Returning " << Norm << " " << Pq0 << " " << Wq0 << " "
<< Pq3 << " " << Wq3 << " " << w << std::endl;
}
if (w != w || std::isnan(w) || w < 0.0) {
w = 0.0;
}
if (w < 1.0 and !fAllowSuppression) {
w = 1.0;
}
return w;
// Use the MINERvA Gaussian method
} else {
/*
* From MINERvA and Daniel Ruterbories:
* Old notes here: *
* http://cdcvs.fnal.gov/cgi-bin/public-cvs/cvsweb-public.cgi/AnalysisFramework/Ana/CCQENu/ana_common/data/?cvsroot=mnvsoft
* These parameters are slightly altered
*
* FRESH:
* 10.5798
* 0.254032
* 0.50834
* 0.0571035
* 0.129051
* 0.875287
*/
if (fDebugStatements) {
std::cout << "Using MINERvA Gaussian" << std::endl;
}
double norm = vals[kPosNorm];
double meanq0 = vals[kPosTilt];
double meanq3 = vals[kPosPq0];
double sigmaq0 = vals[kPosWq0];
double sigmaq3 = vals[kPosPq3];
double corr = vals[kPosWq3];
double z = (q0 - meanq0) * (q0 - meanq0) / (sigmaq0 * sigmaq0) +
(q3 - meanq3) * (q3 - meanq3) / (sigmaq3 * sigmaq3) -
2 * corr * (q0 - meanq0) * (q3 - meanq3) / (sigmaq0 * sigmaq3);
double ret = 1.0;
if ( fabs(1 - corr*corr) < 1.E-5 ) {
return 1.0;
}
if ( (-0.5 * z / (1 - corr*corr)) > 200 or (-0.5 * z / (1 - corr*corr)) < -200 ) {
return 1.0;
} else {
ret = norm * exp( -0.5 * z / (1 - corr*corr) );
}
if (ret != ret or ret < 0.0 or isnan(ret)) {
return 1.0;
}
if (fAllowSuppression) return ret;
return ret + 1.0;
// Need to add 1 to the results
}
return w;
}
void GaussianModeCorr::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void GaussianModeCorr::SetDialValue(int rwenum, double val) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
// Check Handled
if (!IsHandled(curenum)) return;
// CCQE Setting
for (int i = kGaussianCorr_CCQE_norm; i <= kGaussianCorr_CCQE_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_CCQE_norm;
fGausVal_CCQE[index] = val;
fApply_CCQE = true;
}
}
// 2p2h Setting
for (int i = kGaussianCorr_2p2h_norm; i <= kGaussianCorr_2p2h_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_2p2h_norm;
fGausVal_2p2h[index] = val;
fApply_2p2h = true;
}
}
// 2p2h_PPandNN Setting
for (int i = kGaussianCorr_2p2h_PPandNN_norm; i <= kGaussianCorr_2p2h_PPandNN_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_2p2h_PPandNN_norm;
fGausVal_2p2h_PPandNN[index] = val;
fApply_2p2h_PPandNN = true;
}
}
// 2p2h_NP Setting
for (int i = kGaussianCorr_2p2h_NP_norm; i <= kGaussianCorr_2p2h_NP_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_2p2h_NP_norm;
fGausVal_2p2h_NP[index] = val;
fApply_2p2h_NP = true;
}
}
// CC1pi Setting
for (int i = kGaussianCorr_CC1pi_norm; i <= kGaussianCorr_CC1pi_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_CC1pi_norm;
fGausVal_CC1pi[index] = val;
fApply_CC1pi = true;
}
}
if (curenum == kGaussianCorr_AllowSuppression) {
fAllowSuppression = (val > 0.5);
}
}
bool GaussianModeCorr::IsHandled(int rwenum) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case kGaussianCorr_CCQE_norm:
case kGaussianCorr_CCQE_tilt:
case kGaussianCorr_CCQE_Pq0:
case kGaussianCorr_CCQE_Wq0:
case kGaussianCorr_CCQE_Pq3:
case kGaussianCorr_CCQE_Wq3:
case kGaussianCorr_2p2h_norm:
case kGaussianCorr_2p2h_tilt:
case kGaussianCorr_2p2h_Pq0:
case kGaussianCorr_2p2h_Wq0:
case kGaussianCorr_2p2h_Pq3:
case kGaussianCorr_2p2h_Wq3:
case kGaussianCorr_2p2h_PPandNN_norm:
case kGaussianCorr_2p2h_PPandNN_tilt:
case kGaussianCorr_2p2h_PPandNN_Pq0:
case kGaussianCorr_2p2h_PPandNN_Wq0:
case kGaussianCorr_2p2h_PPandNN_Pq3:
case kGaussianCorr_2p2h_PPandNN_Wq3:
case kGaussianCorr_2p2h_NP_norm:
case kGaussianCorr_2p2h_NP_tilt:
case kGaussianCorr_2p2h_NP_Pq0:
case kGaussianCorr_2p2h_NP_Wq0:
case kGaussianCorr_2p2h_NP_Pq3:
case kGaussianCorr_2p2h_NP_Wq3:
case kGaussianCorr_CC1pi_norm:
case kGaussianCorr_CC1pi_tilt:
case kGaussianCorr_CC1pi_Pq0:
case kGaussianCorr_CC1pi_Wq0:
case kGaussianCorr_CC1pi_Pq3:
case kGaussianCorr_CC1pi_Wq3:
case kGaussianCorr_AllowSuppression:
return true;
default:
return false;
}
}
diff --git a/src/Reweight/OscWeightEngine.cxx b/src/Reweight/OscWeightEngine.cxx
index 1224e02..43d9c91 100644
--- a/src/Reweight/OscWeightEngine.cxx
+++ b/src/Reweight/OscWeightEngine.cxx
@@ -1,331 +1,331 @@
// Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
//#define DEBUG_OSC_WE
#include "OscWeightEngine.h"
#include
enum nuTypes {
kNuebarType = -1,
kNumubarType = -2,
kNutaubarType = -3,
kNueType = 1,
kNumuType = 2,
kNutauType = 3,
};
nuTypes GetNuType(int pdg) {
switch (pdg) {
case 16:
return kNutauType;
case 14:
return kNumuType;
case 12:
return kNueType;
case -16:
return kNutaubarType;
case -14:
return kNumubarType;
case -12:
return kNuebarType;
default: { NUIS_ABORT("Attempting to convert \"neutrino pdg\": " << pdg); }
}
}
OscWeightEngine::OscWeightEngine()
:
#ifdef __PROB3PP_ENABLED__
bp(),
#endif
theta12(0.825),
theta13(0.10),
theta23(1.0),
dm12(7.9e-5),
dm23(2.5e-3),
dcp(0.0),
LengthParam(0xdeadbeef),
TargetNuType(0),
ForceFromNuPDG(0) {
Config();
}
void OscWeightEngine::Config() {
std::vector OscParam = Config::QueryKeys("OscParam");
if (OscParam.size() < 1) {
NUIS_ERR(WRN,
"Oscillation parameters specified but no OscParam element "
"configuring the experimental characteristics found.\nExpect at "
"least . Pausing for "
"10...");
sleep(10);
return;
}
if (OscParam[0].Has("baseline_km")) {
LengthParamIsZenith = false;
LengthParam = OscParam[0].GetD("baseline_km");
constant_density = OscParam[0].Has("matter_density")
? OscParam[0].GetD("matter_density")
: 0xdeadbeef;
} else if (OscParam[0].Has("detection_zenith_deg")) {
LengthParamIsZenith = true;
static const double deg2rad = asin(1) / 90.0;
LengthParam = cos(OscParam[0].GetD("detection_zenith_deg") * deg2rad);
} else {
NUIS_ERR(WRN,
"It appeared that you wanted to set up an oscillation weight "
"branch, but it was not correctly configured. You need to specify "
"either: detection_zenith_deg or baseline_km attributes on the "
"OscParam element, and if baseline_km is specified, you can "
"optionally also set matter_density for oscillations through a "
"constant matter density. Pausing for 10...");
sleep(10);
return;
}
dm23 = OscParam[0].Has("dm23") ? OscParam[0].GetD("dm23") : dm23;
theta23 = OscParam[0].Has("sinsq_theta23") ? OscParam[0].GetD("sinsq_theta23")
: theta23;
theta13 = OscParam[0].Has("sinsq_theta13") ? OscParam[0].GetD("sinsq_theta13")
: theta13;
dm12 = OscParam[0].Has("dm12") ? OscParam[0].GetD("dm12") : dm12;
theta12 = OscParam[0].Has("sinsq_theta12") ? OscParam[0].GetD("sinsq_theta12")
: theta12;
dcp = OscParam[0].Has("dcp") ? OscParam[0].GetD("dcp") : dcp;
TargetNuType = OscParam[0].Has("TargetNuPDG")
? GetNuType(OscParam[0].GetI("TargetNuPDG"))
: 0;
ForceFromNuPDG = OscParam[0].Has("ForceFromNuPDG")
? GetNuType(OscParam[0].GetI("ForceFromNuPDG"))
: 0;
NUIS_LOG(FIT, "Configured oscillation weighter:");
if (LengthParamIsZenith) {
NUIS_LOG(FIT,
"Earth density profile with detection cos(zenith) = " << LengthParam);
} else {
if (constant_density != 0xdeadbeef) {
NUIS_LOG(FIT,
"Constant density with experimental baseline = " << LengthParam);
} else {
NUIS_LOG(FIT,
"Vacuum oscillations with experimental baseline = " << LengthParam);
}
}
params[0] = dm23;
params[1] = theta23;
params[2] = theta13;
params[3] = dm12;
params[4] = theta12;
params[5] = dcp;
NUIS_LOG(FIT, "\tdm23 : " << params[0]);
NUIS_LOG(FIT, "\tsinsq_theta23: " << params[1]);
NUIS_LOG(FIT, "\tsinsq_theta13: " << params[2]);
NUIS_LOG(FIT, "\tdm12 : " << params[3]);
NUIS_LOG(FIT, "\tsinsq_theta12: " << params[4]);
NUIS_LOG(FIT, "\tdcp : " << params[5]);
if (TargetNuType) {
NUIS_LOG(FIT, "\tTargetNuType: " << TargetNuType);
}
if (ForceFromNuPDG) {
NUIS_LOG(FIT, "\tForceFromNuPDG: " << ForceFromNuPDG);
}
#ifdef __PROB3PP_ENABLED__
bp.SetMNS(params[theta12_idx], params[theta13_idx], params[theta23_idx],
params[dm12_idx], params[dm23_idx], params[dcp_idx], 1, true, 2);
bp.DefinePath(LengthParam, 0);
if (LengthParamIsZenith) {
NUIS_LOG(FIT, "\tBaseline : " << (bp.GetBaseline() / 100.0) << " km.");
}
#endif
}
void OscWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef DEBUG_OSC_WE
std::cout << "IncludeDial: " << name << " at " << startval << std::endl;
#endif
int dial = SystEnumFromString(name);
if (!dial) {
NUIS_ABORT("OscWeightEngine passed dial: " << name
<< " that it does not understand.");
}
params[dial - 1] = startval;
}
void OscWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef DEBUG_OSC_WE
- std::cout << "SetDial: " << (nuisenum % 1000) << " at " << val << std::endl;
+ std::cout << "SetDial: " << (nuisenum % NUIS_DIAL_OFFSET) << " at " << val << std::endl;
#endif
- fHasChanged = (params[(nuisenum % 1000) - 1] - val) >
+ fHasChanged = (params[(nuisenum % NUIS_DIAL_OFFSET) - 1] - val) >
std::numeric_limits::epsilon();
- params[(nuisenum % 1000) - 1] = val;
+ params[(nuisenum % NUIS_DIAL_OFFSET) - 1] = val;
}
void OscWeightEngine::SetDialValue(std::string name, double val) {
#ifdef DEBUG_OSC_WE
std::cout << "SetDial: " << name << " at " << val << std::endl;
#endif
int dial = SystEnumFromString(name);
if (!dial) {
NUIS_ABORT("OscWeightEngine passed dial: " << name
<< " that it does not understand.");
}
fHasChanged =
(params[dial - 1] - val) > std::numeric_limits::epsilon();
params[dial - 1] = val;
}
bool OscWeightEngine::IsDialIncluded(std::string name) {
return SystEnumFromString(name);
}
bool OscWeightEngine::IsDialIncluded(int nuisenum) {
- return ((nuisenum % 1000) > 0) && ((nuisenum % 1000) < 6);
+ return ((nuisenum % NUIS_DIAL_OFFSET) > 0) && ((nuisenum % NUIS_DIAL_OFFSET) < 6);
}
double OscWeightEngine::GetDialValue(std::string name) {
int dial = SystEnumFromString(name);
if (!dial) {
NUIS_ABORT("OscWeightEngine passed dial: " << name
<< " that it does not understand.");
}
return params[dial - 1];
}
double OscWeightEngine::GetDialValue(int nuisenum) {
- if (!(nuisenum % 1000) || (nuisenum % 1000) > 6) {
+ if (!(nuisenum % NUIS_DIAL_OFFSET) || (nuisenum % NUIS_DIAL_OFFSET) > 6) {
NUIS_ABORT("OscWeightEngine passed dial enum: "
- << (nuisenum % 1000)
+ << (nuisenum % NUIS_DIAL_OFFSET)
<< " that it does not understand, expected [1,6].");
}
- return params[(nuisenum % 1000) - 1];
+ return params[(nuisenum % NUIS_DIAL_OFFSET) - 1];
}
void OscWeightEngine::Reconfigure(bool silent) { fHasChanged = false; };
bool OscWeightEngine::NeedsEventReWeight() {
if (fHasChanged) {
return true;
}
return false;
}
double OscWeightEngine::CalcWeight(BaseFitEvt* evt) {
static bool Warned = false;
if (evt->probe_E == 0xdeadbeef) {
if (!Warned) {
NUIS_ERR(WRN,
"Oscillation weights asked for but using 'litemode' or "
"unsupported generator input. Pasuing for 10...");
sleep(10);
Warned = true;
}
return 1;
}
return CalcWeight(evt->probe_E * 1E-3, evt->probe_pdg);
}
double OscWeightEngine::CalcWeight(double ENu, int PDGNu, int TargetPDGNu) {
if (LengthParam == 0xdeadbeef) { // not configured.
return 1;
}
#ifdef __PROB3PP_ENABLED__
int NuType = (ForceFromNuPDG != 0) ? ForceFromNuPDG : GetNuType(PDGNu);
bp.SetMNS(params[theta12_idx], params[theta13_idx], params[theta23_idx],
params[dm12_idx], params[dm23_idx], params[dcp_idx], ENu, true,
NuType);
int pmt = 0;
double prob_weight = 1;
TargetPDGNu = (TargetPDGNu == -1) ? (TargetNuType ? TargetNuType : NuType)
: GetNuType(TargetPDGNu);
if (LengthParamIsZenith) { // Use earth density
bp.DefinePath(LengthParam, 0);
bp.propagate(NuType);
pmt = 0;
prob_weight = bp.GetProb(NuType, TargetPDGNu);
} else {
if (constant_density != 0xdeadbeef) {
bp.propagateLinear(NuType, LengthParam, constant_density);
pmt = 1;
prob_weight = bp.GetProb(NuType, TargetPDGNu);
} else {
pmt = 2;
prob_weight =
bp.GetVacuumProb(NuType, TargetPDGNu, ENu * 1E-3, LengthParam);
}
}
#ifdef DEBUG_OSC_WE
if (prob_weight != prob_weight) {
NUIS_ABORT("Calculated bad prob weight: " << prob_weight << "(Osc Type: " << pmt
<< " -- " << NuType << " -> "
<< TargetPDGNu << ")");
}
if (prob_weight > 1) {
NUIS_ABORT("Calculated bad prob weight: " << prob_weight << "(Osc Type: " << pmt
<< " -- " << NuType << " -> "
<< TargetPDGNu << ")");
}
std::cout << NuType << " -> " << TargetPDGNu << ": " << ENu << " = "
<< prob_weight << "%%." << std::endl;
#endif
return prob_weight;
#else
return 1;
#endif
}
int OscWeightEngine::SystEnumFromString(std::string const& name) {
if (name == "dm23") {
return 1;
} else if (name == "sinsq_theta23") {
return 2;
} else if (name == "sinsq_theta13") {
return 3;
} else if (name == "dm12") {
return 4;
} else if (name == "sinsq_theta12") {
return 5;
} else if (name == "dcp") {
return 6;
} else {
return 0;
}
}
void OscWeightEngine::Print() {
std::cout << "OscWeightEngine: " << std::endl;
std::cout << "\t theta12: " << params[theta12_idx] << std::endl;
std::cout << "\t theta13: " << params[theta13_idx] << std::endl;
std::cout << "\t theta23: " << params[theta23_idx] << std::endl;
std::cout << "\t dm12: " << params[dm12_idx] << std::endl;
std::cout << "\t dm23: " << params[dm23_idx] << std::endl;
std::cout << "\t dcp: " << params[dcp_idx] << std::endl;
}
diff --git a/src/Reweight/WeightEngineBase.h b/src/Reweight/WeightEngineBase.h
index 5a6e6cf..df326bd 100644
--- a/src/Reweight/WeightEngineBase.h
+++ b/src/Reweight/WeightEngineBase.h
@@ -1,56 +1,57 @@
#ifndef WEIGHTENGINE_BASE_H
#define WEIGHTENGINE_BASE_H
#include "BaseFitEvt.h"
#include "FitLogger.h"
#include "FitUtils.h"
#include
#include
#include
#include
#include
#include
#include