diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0154ce4..e242e49 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,232 +1,232 @@
# 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 .
################################################################################
cmake_minimum_required (VERSION 2.8 FATAL_ERROR)
#Use the compilers found in the path
find_program(CMAKE_C_COMPILER NAMES $ENV{CC} gcc PATHS ENV PATH NO_DEFAULT_PATH)
find_program(CMAKE_CXX_COMPILER NAMES $ENV{CXX} g++ PATHS ENV PATH NO_DEFAULT_PATH)
project(NUISANCE)
include(ExternalProject)
enable_language(Fortran)
set (NUISANCE_VERSION_MAJOR 2)
set (NUISANCE_VERSION_MINOR 7)
set (NUISANCE_VERSION_REVISION 0)
set (NUISANCE_VERSION_STRING "v${NUISANCE_VERSION_MAJOR}r${NUISANCE_VERSION_MINOR}")
if(${NUISANCE_VERSION_REVISION} STRGREATER "0")
set (NUISANCE_VERSION_STRING "${NUISANCE_VERSION_STRING}p${NUISANCE_VERSION_REVISION}")
endif()
#Set this to TRUE to enable build debugging messages
set(BUILD_DEBUG_MSGS TRUE)
include(${CMAKE_SOURCE_DIR}/cmake/cmessage.cmake)
include(${CMAKE_SOURCE_DIR}/cmake/cacheVariables.cmake)
cmessage(STATUS "CMAKE_INSTALL_PREFIX: \"${CMAKE_INSTALL_PREFIX}\"")
cmessage(STATUS "CMAKE_BUILD_TYPE: \"${CMAKE_BUILD_TYPE}\"")
################################################################################
# Check Dependencies
################################################################################
################################## ROOT ######################################
include(${CMAKE_SOURCE_DIR}/cmake/ROOTSetup.cmake)
################################# HEPMC ######################################
include(${CMAKE_SOURCE_DIR}/cmake/HepMC.cmake)
############################ Reweight Engines ################################
include(${CMAKE_SOURCE_DIR}/cmake/ReweightEnginesSetup.cmake)
############################ Other Generators ################################
include(${CMAKE_SOURCE_DIR}/cmake/GiBUUSetup.cmake)
if(USE_NUANCE)
LIST(APPEND EXTRA_CXX_FLAGS -D__NUANCE_ENABLED__)
endif()
################################# Pythia6/8 ####################################
include(${CMAKE_SOURCE_DIR}/cmake/pythia6Setup.cmake)
include(${CMAKE_SOURCE_DIR}/cmake/pythia8Setup.cmake)
################################# gperftools ###################################
include(${CMAKE_SOURCE_DIR}/cmake/gperfSetup.cmake)
if(NOT NOTEST)
enable_testing()
endif()
SET(GENERATOR_SUPPORT)
foreach(gen NEUT;NuWro;GENIE;GiBUU;NUANCE)
if(USE_${gen})
SET(GENERATOR_SUPPORT "${GENERATOR_SUPPORT}${gen} ")
endif()
endforeach(gen)
cmessage(STATUS "Generator Input Support: ${GENERATOR_SUPPORT}")
set(MINCODE
Routines
FCN)
set(CORE
MCStudies
Genie
FitBase
Config
Logger
InputHandler
Splines
Reweight
Utils
Statistical
#Devel
Smearceptance
)
LIST(APPEND ALLEXPERIMENTS
ANL
ArgoNeuT
BEBC
BNL
Electron
FNAL
GGM
K2K
MINERvA
MiniBooNE
SciBooNE
T2K)
foreach(exp ${ALLEXPERIMENTS})
if(NOT NO_${exp})
LIST(APPEND EXPERIMENTS_TO_BUILD ${exp})
else()
LIST(REVERSE EXTRA_CXX_FLAGS)
LIST(APPEND EXTRA_CXX_FLAGS -D__NO_${exp}__)
LIST(REVERSE EXTRA_CXX_FLAGS)
endif()
endforeach()
################################## COMPILER ####################################
include(${CMAKE_SOURCE_DIR}/cmake/c++CompilerSetup.cmake)
################################### doxygen ###################################
include(${CMAKE_SOURCE_DIR}/cmake/docsSetup.cmake)
################################################################################
set(MINIMUM_INCLUDE_DIRECTORIES)
LIST(APPEND MINIMUM_INCLUDE_DIRECTORIES
${RWENGINE_INCLUDE_DIRECTORIES}
${CMAKE_SOURCE_DIR}/src/FitBase
${CMAKE_SOURCE_DIR}/src/Reweight
${CMAKE_SOURCE_DIR}/src/InputHandler
${CMAKE_SOURCE_DIR}/src/Config
${CMAKE_SOURCE_DIR}/src/Logger
${CMAKE_SOURCE_DIR}/src/Statistical
${CMAKE_SOURCE_DIR}/src/Splines
${CMAKE_SOURCE_DIR}/src/Utils
${CMAKE_SOURCE_DIR}/src/Genie)
cmessage(DEBUG "Base include directories: ${MINIMUM_INCLUDE_DIRECTORIES}")
set(EXP_INCLUDE_DIRECTORIES)
foreach(edir ${EXPERIMENTS_TO_BUILD})
LIST(APPEND EXP_INCLUDE_DIRECTORIES ${CMAKE_SOURCE_DIR}/src/${edir})
endforeach()
cmessage(DEBUG "Included experiments: ${EXP_INCLUDE_DIRECTORIES}")
foreach(mdir ${MINCODE})
cmessage (DEBUG "Configuring directory: src/${mdir}")
add_subdirectory(src/${mdir})
endforeach()
foreach(edir ${EXPERIMENTS_TO_BUILD})
cmessage (DEBUG "Configuring directory: src/${edir}")
add_subdirectory(src/${edir})
endforeach()
foreach(cdir ${CORE})
cmessage (DEBUG "Configuring directory: src/${cdir}")
add_subdirectory(src/${cdir})
endforeach()
cmessage(DEBUG "Module targets: ${MODULETargets}")
+set(MODULETargets "-Wl,--start-group;${MODULETargets};-Wl,--start-group")
+
add_subdirectory(app)
add_subdirectory(src/Tests)
configure_file(cmake/setup.sh.in
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.sh" @ONLY)
install(FILES
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.sh" DESTINATION
${CMAKE_INSTALL_PREFIX})
configure_file(cmake/MakeBinaryBlob.in
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/MakeBinaryBlob" @ONLY)
install(PROGRAMS
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/MakeBinaryBlob" DESTINATION
bin)
if(USE_DYNSAMPLES)
SET(ALL_INCLUDES ${MINIMUM_INCLUDE_DIRECTORIES})
LIST(APPEND ALL_INCLUDES ${CMAKE_SOURCE_DIR}/src/Smearceptance)
LIST(APPEND ALL_INCLUDES ${EXP_INCLUDE_DIRECTORIES})
string(REPLACE ";" " -I" ALL_INCLUDES_STR "${ALL_INCLUDES}")
cmessage(DEBUG ${CMAKE_DEPENDLIB_FLAGS})
string(REPLACE "-levent " "" CMAKE_DEPENDLIB_FLAGS_NEW ${CMAKE_DEPENDLIB_FLAGS})
set(CMAKE_DEPENDLIB_FLAGS ${CMAKE_DEPENDLIB_FLAGS_NEW})
cmessage(DEBUG ${CMAKE_DEPENDLIB_FLAGS})
- string(REPLACE ";" " -l" ALL_MODULETARGETS_STR "${MODULETargets}")
-
configure_file(cmake/BuildDynamicSample.in
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/BuildDynamicSample" @ONLY)
install(PROGRAMS
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/BuildDynamicSample" DESTINATION
bin)
configure_file(cmake/BuildDynamicSmearcepter.in
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/BuildDynamicSmearcepter" @ONLY)
install(PROGRAMS
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/BuildDynamicSmearcepter" DESTINATION
bin)
endif()
install(PROGRAMS
"${PROJECT_SOURCE_DIR}/scripts/nuiscardgen" DESTINATION
bin)
install(PROGRAMS
"${PROJECT_SOURCE_DIR}/scripts/nuissamples" DESTINATION
bin)
diff --git a/cmake/DUNERwtSetup.cmake b/cmake/DUNERwtSetup.cmake
new file mode 100644
index 0000000..1f72c03
--- /dev/null
+++ b/cmake/DUNERwtSetup.cmake
@@ -0,0 +1,37 @@
+# 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 .
+################################################################################
+
+if(NUSYST_ROOT STREQUAL "")
+ cmessage(FATAL_ERROR "Variable NUSYST_ROOT is not defined. Either configure with -DNUSYST_ROOT or \"\$ export NUSYST_ROOT=/path/to/nusystematics\". This must be set to point to a prebuilt NuSystematics instance.")
+endif()
+
+if(SYSTTOOLS_ROOT STREQUAL "")
+ cmessage(FATAL_ERROR "Variable SYSTTOOLS_ROOT is not defined. Either configure with -DSYSTTOOLS_ROOT or \"\$ export SYSTTOOLS_ROOT=/path/to/systematicstools\". This must be set to point to a prebuilt ART Systematics Tools instance.")
+endif()
+
+LIST(APPEND EXTRA_CXX_FLAGS -D__DUNERWT_ENABLED__ -DNO_ART -std=c++1y -Wno-deprecated-declarations -Wno-deprecated)
+
+LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${NUSYST_ROOT}/include ${SYSTTOOLS_ROOT}/include)
+
+LIST(APPEND EXTRA_LINK_DIRS ${NUSYST_ROOT}/lib ${SYSTTOOLS_ROOT}/lib)
+LIST(APPEND EXTRA_LIBS nusystematics_systproviders
+ systematicstools_interface
+ systematicstools_interpreters
+ systematicstools_systproviders
+ systematicstools_utility)
diff --git a/cmake/GENIESetup.cmake b/cmake/GENIESetup.cmake
index 11c64b3..3126bfc 100644
--- a/cmake/GENIESetup.cmake
+++ b/cmake/GENIESetup.cmake
@@ -1,168 +1,169 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
# TODO
# check system for libxml2
# check whether we need the includes
# check if we can use a subset of the GENIE libraries
################################################################################
# Check Dependencies
################################################################################
################################# GENIE ######################################
if(GENIE STREQUAL "")
cmessage(FATAL_ERROR "Variable GENIE is not defined. "
"The location of a pre-built GENIE install must be defined either as"
" $ cmake -DGENIE=/path/to/GENIE or as and environment vairable"
" $ export GENIE=/path/to/GENIE")
endif()
if (BUILD_GEVGEN)
cmessage(STATUS "Building custom gevgen")
LIST(APPEND EXTRA_CXX_FLAGS -D__GEVGEN_ENABLED__)
endif()
# Extract GENIE VERSION
if (GENIE_VERSION STREQUAL "AUTO")
execute_process (COMMAND ${CMAKE_SOURCE_DIR}/cmake/getgenieversion.sh ${GENIE}
OUTPUT_VARIABLE GENIE_VERSION OUTPUT_STRIP_TRAILING_WHITESPACE)
endif()
execute_process (COMMAND genie-config
--libs OUTPUT_VARIABLE GENIE_LD_FLAGS_STR OUTPUT_STRIP_TRAILING_WHITESPACE)
execute_process (COMMAND genie-config
--topsrcdir OUTPUT_VARIABLE GENIE_INCLUDES_DIR OUTPUT_STRIP_TRAILING_WHITESPACE)
string(REGEX MATCH "-L\([^ ]+\) \(.*\)$" PARSE_GENIE_LIBS_MATCH ${GENIE_LD_FLAGS_STR})
cmessage(DEBUG "genie-config --libs: ${GENIE_LD_FLAGS_STR}")
if(NOT PARSE_GENIE_LIBS_MATCH)
cmessage(FATAL_ERROR "Expected to be able to parse the result of genie-config --libs to a lib directory and a list of libraries to include, but got: \"${GENIE_LD_FLAGS_STR}\"")
endif()
set(GENIE_LIB_DIR ${CMAKE_MATCH_1})
set(GENIE_LIBS_RAW ${CMAKE_MATCH_2})
string(REPLACE "-l" "" GENIE_LIBS_STRIPED "${GENIE_LIBS_RAW}")
cmessage(STATUS "GENIE version : ${GENIE_VERSION}")
cmessage(STATUS "GENIE libdir : ${GENIE_LIB_DIR}")
cmessage(STATUS "GENIE libs : ${GENIE_LIBS_STRIPED}")
string(REGEX MATCH "ReinSeghal" WASMATCHED ${GENIE_LIBS_STRIPED})
if(WASMATCHED AND GENIE_VERSION STREQUAL "210")
set(GENIE_SEHGAL ${GENIE_LIBS_STRIPED})
STRING(REPLACE "ReinSeghal" "ReinSehgal" GENIE_LIBS_STRIPED ${GENIE_SEHGAL})
cmessage(DEBUG "Fixed inconsistency in library naming: ${GENIE_LIBS_STRIPED}")
endif()
string(REGEX MATCH "ReWeight" WASMATCHED ${GENIE_LIBS_STRIPED})
if(NOT WASMATCHED)
set(GENIE_LIBS_STRIPED "GReWeight ${GENIE_LIBS_STRIPED}")
cmessage(DEBUG "Force added ReWeight library: ${GENIE_LIBS_STRIPED}")
endif()
string(REPLACE " " ";" GENIE_LIBS_LIST "-Wl,--start-group ${GENIE_LIBS_STRIPED} -Wl,--end-group")
cmessage(DEBUG "genie-config --libs -- MATCH1: ${CMAKE_MATCH_1}")
cmessage(DEBUG "genie-config --libs -- MATCH2: ${CMAKE_MATCH_2}")
cmessage(DEBUG "genie-config --libs -- libs stripped: ${GENIE_LIBS_STRIPED}")
cmessage(DEBUG "genie-config --libs -- libs list: ${GENIE_LIBS_LIST}")
################################ LHAPDF ######################################
if(LHAPDF_LIB STREQUAL "")
cmessage(FATAL_ERROR "Variable LHAPDF_LIB is not defined. The location of a pre-built lhapdf install must be defined either as $ cmake -DLHAPDF_LIB=/path/to/LHAPDF_libraries or as and environment vairable $ export LHAPDF_LIB=/path/to/LHAPDF_libraries")
endif()
if(LHAPDF_INC STREQUAL "")
cmessage(FATAL_ERROR "Variable LHAPDF_INC is not defined. The location of a pre-built lhapdf install must be defined either as $ cmake -DLHAPDF_INC=/path/to/LHAPDF_includes or as and environment vairable $ export LHAPDF_INC=/path/to/LHAPDF_includes")
endif()
if(LHAPATH STREQUAL "")
cmessage(FATAL_ERROR "Variable LHAPATH is not defined. The location of a the LHAPATH directory must be defined either as $ cmake -DLHAPATH=/path/to/LHAPATH or as and environment variable $ export LHAPATH=/path/to/LHAPATH")
endif()
################################ LIBXML ######################################
if(LIBXML2_LIB STREQUAL "")
cmessage(FATAL_ERROR "Variable LIBXML2_LIB is not defined. The location of a pre-built libxml2 install must be defined either as $ cmake -DLIBXML2_LIB=/path/to/LIBXML2_libraries or as and environment vairable $ export LIBXML2_LIB=/path/to/LIBXML2_libraries")
endif()
if(LIBXML2_INC STREQUAL "")
cmessage(FATAL_ERROR "Variable LIBXML2_INC is not defined. The location of a pre-built libxml2 install must be defined either as $ cmake -DLIBXML2_INC=/path/to/LIBXML2_includes or as and environment vairable $ export LIBXML2_INC=/path/to/LIBXML2_includes")
endif()
############################### log4cpp ######################################
if(LOG4CPP_LIB STREQUAL "")
find_program(LOG4CPPCFG log4cpp-config)
if(NOT LOG4CPPCFG STREQUAL "LOG4CPPCFG-NOTFOUND")
execute_process (COMMAND ${LOG4CPPCFG}
--pkglibdir OUTPUT_VARIABLE LOG4CPP_LIB OUTPUT_STRIP_TRAILING_WHITESPACE)
else()
message(FATAL_ERROR "Variable LOG4CPP_LIB is not defined. The location of a pre-built log4cpp install must be defined either as $ cmake -DLOG4CPP_LIB=/path/to/LOG4CPP_libraries or as and environment vairable $ export LOG4CPP_LIB=/path/to/LOG4CPP_libraries")
endif()
endif()
if(LOG4CPP_INC STREQUAL "")
find_program(LOG4CPPCFG log4cpp-config)
if(NOT LOG4CPPCFG STREQUAL "LOG4CPPCFG-NOTFOUND")
execute_process (COMMAND ${LOG4CPPCFG}
--pkgincludedir OUTPUT_VARIABLE LOG4CPP_INC OUTPUT_STRIP_TRAILING_WHITESPACE)
else()
message(FATAL_ERROR "Variable LOG4CPP_INC is not defined. The location of a pre-built log4cpp install must be defined either as $ cmake -DGENIE_LOG4CPP_INC=/path/to/LOG4CPP_includes or as and environment vairable $ export LOG4CPP_INC=/path/to/LOG4CPP_includes")
endif()
endif()
################################################################################
LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_ENABLED__ -D__GENIE_VERSION__=${GENIE_VERSION})
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES
${GENIE_INCLUDES_DIR}
${GENIE_INCLUDES_DIR}/GHEP
${GENIE_INCLUDES_DIR}/Ntuple
${GENIE_INCLUDES_DIR}/ReWeight
${GENIE_INCLUDES_DIR}/Apps
${GENIE_INCLUDES_DIR}/FluxDrivers
${GENIE_INCLUDES_DIR}/EVGDrivers
${LHAPDF_INC}
${LIBXML2_INC}
${LOG4CPP_INC})
SAYVARS()
LIST(APPEND EXTRA_LINK_DIRS
${GENIE_LIB_DIR}
${LHAPDF_LIB}
${LIBXML2_LIB}
${LOG4CPP_LIB})
-LIST(APPEND EXTRA_LIBS -Wl,--start-group)
+#-Wl,--no-as-needed used to ensuring dynamic linking of GENIE algorithm libraries so that TClass::GetClass has access to the dictionaries at runtime
+LIST(APPEND EXTRA_LIBS -Wl,--no-as-needed -Wl,--start-group)
LIST(APPEND EXTRA_LIBS ${GENIE_LIBS_LIST})
-LIST(APPEND EXTRA_LIBS -Wl,--end-group)
+LIST(APPEND EXTRA_LIBS -Wl,--end-group -Wl,-as-needed)
LIST(APPEND EXTRA_LIBS LHAPDF xml2 log4cpp)
if(USE_PYTHIA8)
set(NEED_PYTHIA8 TRUE)
set(NEED_ROOTPYTHIA8 TRUE)
else()
set(NEED_PYTHIA6 TRUE)
set(NEED_ROOTPYTHIA6 TRUE)
endif()
set(NEED_ROOTEVEGEN TRUE)
SET(USE_GENIE TRUE CACHE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. " FORCE)
diff --git a/cmake/ReweightEnginesSetup.cmake b/cmake/ReweightEnginesSetup.cmake
index 16b9e85..231ce12 100644
--- a/cmake/ReweightEnginesSetup.cmake
+++ b/cmake/ReweightEnginesSetup.cmake
@@ -1,87 +1,93 @@
# 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 .
################################################################################
################################## 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()
################################## 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()
+################################## NIWG ######################################
+if(USE_DUNERWT)
+ include(${CMAKE_SOURCE_DIR}/cmake/DUNERwtSetup.cmake)
+ cmessage(STATUS "Using DUNE Reweight engine.")
+ set(USE_DUNERWT TRUE CACHE BOOL "Whether to enable DUNE ReWeight 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()
################################## 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()
################################################################################
################################ 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/c++CompilerSetup.cmake b/cmake/c++CompilerSetup.cmake
index bd9adcb..518dc24 100644
--- a/cmake/c++CompilerSetup.cmake
+++ b/cmake/c++CompilerSetup.cmake
@@ -1,130 +1,113 @@
# 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 .
################################################################################
if(USE_OMP)
LIST(APPEND EXTRA_CXX_FLAGS -fopenmp)
endif()
if(USE_DYNSAMPLES)
LIST(APPEND EXTRA_LIBS dl)
LIST(APPEND EXTRA_CXX_FLAGS -D__USE_DYNSAMPLES__)
endif()
set(CXX_WARNINGS -Wall )
cmessage(DEBUG "EXTRA_CXX_FLAGS: ${EXTRA_CXX_FLAGS}")
string(REPLACE ";" " " STR_EXTRA_CXX_FLAGS "${EXTRA_CXX_FLAGS}")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${STR_EXTRA_CXX_FLAGS} ${CXX_WARNINGS}")
set(CMAKE_Fortran_FLAGS_RELEASE "-fPIC")
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0")
if(USE_DYNSAMPLES)
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fPIC")
set(CMAKE_Fortran_FLAGS_DEBUG "-fPIC")
endif()
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fPIC -O3")
if(CMAKE_BUILD_TYPE MATCHES DEBUG)
set(CURRENT_CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS_DEBUG})
elseif(CMAKE_BUILD_TYPE MATCHES RELEASE)
set(CURRENT_CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS_RELEASE})
else()
cmessage(FATAL_ERROR "[ERROR]: Unknown CMAKE_BUILD_TYPE (\"${CMAKE_BUILD_TYPE}\"): Should be \"DEBUG\" or \"RELEASE\".")
endif()
SET(STR_EXTRA_LINK_DIRS)
if(NOT EXTRA_LINK_DIRS STREQUAL "")
string(REPLACE ";" " -L" STR_EXTRA_LINK_DIRS "-L${EXTRA_LINK_DIRS}")
endif()
SET(STR_EXTRA_LIBS)
if(NOT EXTRA_LIBS STREQUAL "")
SET(STR_EXTRA_LIBS_NO_SCRUB_LINKOPTS)
string(REPLACE ";" " -l" STR_EXTRA_LIBS_NO_SCRUB_LINKOPTS "-l${EXTRA_LIBS}")
string(REPLACE "-l-" "-" STR_EXTRA_LIBS ${STR_EXTRA_LIBS_NO_SCRUB_LINKOPTS})
endif()
SET(STR_EXTRA_SHAREDOBJS)
if(NOT EXTRA_SHAREDOBJS STREQUAL "")
string(REPLACE ";" " " STR_EXTRA_SHAREDOBJS "${EXTRA_SHAREDOBJS}")
endif()
SET(STR_EXTRA_LINK_FLAGS)
if(NOT EXTRA_LINK_FLAGS STREQUAL "")
string(REPLACE ";" " " STR_EXTRA_LINK_FLAGS "${EXTRA_LINK_FLAGS}")
endif()
cmessage(DEBUG "EXTRA_LINK_DIRS: ${STR_EXTRA_LINK_DIRS}")
cmessage(DEBUG "EXTRA_LIBS: ${STR_EXTRA_LIBS}")
cmessage(DEBUG "EXTRA_SHAREDOBJS: ${STR_EXTRA_SHAREDOBJS}")
cmessage(DEBUG "EXTRA_LINK_FLAGS: ${STR_EXTRA_LINK_FLAGS}")
if(NOT STR_EXTRA_LINK_DIRS STREQUAL "" AND NOT STR_EXTRA_LIBS STREQUAL "")
SET(CMAKE_DEPENDLIB_FLAGS "${STR_EXTRA_LINK_DIRS} ${STR_EXTRA_LIBS}")
endif()
-if(USE_NEUT)
- foreach(OBJ ${NEUT_ROOT_LIBS})
- if(NOT CMAKE_DEPENDLIB_FLAGS STREQUAL "")
- SET(CMAKE_DEPENDLIB_FLAGS "${CMAKE_DEPENDLIB_FLAGS} ${OBJ}")
- else()
- SET(CMAKE_DEPENDLIB_FLAGS "${OBJ}")
- endif()
- endforeach()
- foreach(OBJ ${NEUT_ROOT_LIBS})
- if(NOT CMAKE_DEPENDLIB_FLAGS STREQUAL "")
- SET(CMAKE_DEPENDLIB_FLAGS "${CMAKE_DEPENDLIB_FLAGS} ${OBJ}")
- else()
- SET(CMAKE_DEPENDLIB_FLAGS "${OBJ}")
- endif()
- endforeach()
-endif()
-
if(NOT EXTRA_SHAREDOBJS STREQUAL "")
if(NOT STR_EXTRA_LINK_FLAGS STREQUAL "")
SET(STR_EXTRA_LINK_FLAGS "${STR_EXTRA_SHAREDOBJS} ${STR_EXTRA_LINK_FLAGS}")
else()
SET(STR_EXTRA_LINK_FLAGS "${STR_EXTRA_SHAREDOBJS}")
endif()
endif()
if(NOT EXTRA_LINK_FLAGS STREQUAL "")
if(NOT CMAKE_LINK_FLAGS STREQUAL "")
SET(CMAKE_LINK_FLAGS "${CMAKE_LINK_FLAGS} ${STR_EXTRA_LINK_FLAGS}")
else()
SET(CMAKE_LINK_FLAGS "${STR_EXTRA_LINK_FLAGS}")
endif()
endif()
if(USE_OMP)
cmessage(FATAL_ERROR "No OMP features currently enabled so this is a FATAL_ERROR to let you know that you don't gain anything with this declaration.")
endif()
if (VERBOSE)
cmessage (STATUS "C++ Compiler : ${CXX_COMPILER_NAME}")
cmessage (STATUS " flags : ${CMAKE_CXX_FLAGS}")
cmessage (STATUS " Release flags : ${CMAKE_CXX_FLAGS_RELEASE}")
cmessage (STATUS " Debug flags : ${CMAKE_CXX_FLAGS_DEBUG}")
cmessage (STATUS " Link Flags : ${CMAKE_LINK_FLAGS}")
cmessage (STATUS " Lib Flags : ${CMAKE_DEPENDLIB_FLAGS}")
endif()
diff --git a/cmake/cacheVariables.cmake b/cmake/cacheVariables.cmake
index 3d52cee..6ab9bb0 100644
--- a/cmake/cacheVariables.cmake
+++ b/cmake/cacheVariables.cmake
@@ -1,216 +1,220 @@
# 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(EXTRA_SETUP_SCRIPT "" PATH "The path to an extra script to inject into the NUISANCE setup script. <>")
CheckAndSetDefaultCache(USE_MINIMIZER TRUE INTERNAL "Whether we are using the ROOT minimization libraries. ")
CheckAndSetDefaultCache(USE_ROOT6 FALSE INTERNAL "Whether we are using the ROOT 6. ")
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. ")
CheckAndSetDefaultCache(IS_NEUT_54 FALSE BOOL "Whether to enabled NEUT is version 5.4 or greater. ")
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. ")
CheckAndSetDefaultCache(GENIE_VERSION "AUTO" STRING "GENIE Version ")
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_DUNERWT FALSE BOOL "Whether to enable DUNE ReWeight support. Requires external libraries. ")
+CheckAndSetDefaultEnv(NUSYST_ROOT "" PATH "Path to installed NuSystematics. Overrides environment variable \$NUSYST_ROOT. <>" DUNERWT_ROOT)
+CheckAndSetDefaultEnv(SYSTTOOLS_ROOT "" PATH "Path to installed ART Systematics Tools. Overrides environment variable \$SYSTTOOLS_ROOT. <>" SYSTTOOLS_ROOT)
+
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 FALSE BOOL "Whether to perform the update target for external dependencies. Note this may produce errors for CMake < 3.8 where a bug was fixed for the feature that this option invokes. ")
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...). ")
CheckAndSetDefaultCache(USE_DYNSAMPLES TRUE BOOL "Whether to enable the dynamic sample loader. ")
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/InputHandler/BaseFitEvt.cxx b/src/InputHandler/BaseFitEvt.cxx
index 4c0cf07..8101b5b 100644
--- a/src/InputHandler/BaseFitEvt.cxx
+++ b/src/InputHandler/BaseFitEvt.cxx
@@ -1,243 +1,242 @@
// 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 .
-*******************************************************************************/
+ * 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 "BaseFitEvt.h"
#include "FitParticle.h"
#include "GIBUUInputHandler.h"
BaseFitEvt::BaseFitEvt() {
Mode = 0;
probe_E = 0xdeadbeef;
probe_pdg = 0;
Weight = 1.0;
InputWeight = 1.0;
RWWeight = 1.0;
CustomWeight = 1.0;
SavedRWWeight = 1.0;
fSplineCoeff = NULL;
fSplineRead = NULL;
fGenInfo = NULL;
fType = 9999;
#ifdef __NEUT_ENABLED__
fNeutVect = NULL;
#endif
#ifdef __NUWRO_ENABLED__
#ifndef __USE_NUWRO_SRW_EVENTS__
fNuwroEvent = NULL;
#endif
#endif
#ifdef __GENIE_ENABLED__
genie_event = NULL;
- genie_record = NULL;
+ #ifdef __DUNERWT_ENABLED__
+ HasDUNERwtPolyResponses = false;
+ HasDUNERwtResponses = false;
+ #endif
#endif
#ifdef __NUANCE_ENABLED__
nuance_event = NULL;
#endif
#ifdef __GiBUU_ENABLED__
GiRead = NULL;
#endif
};
-BaseFitEvt::~BaseFitEvt() {
-#ifdef __NEUT_ENABLED__
- if (fNeutVect) delete fNeutVect;
-#endif
-
-#ifdef __NUWRO_ENABLED__
-#ifndef __USE_NUWRO_SRW_EVENTS__
- if (fNuwroEvent) delete fNuwroEvent;
-#endif
-#endif
-
-#ifdef __GENIE_ENABLED__
- if (genie_event) delete genie_event;
-#endif
+BaseFitEvt::~BaseFitEvt(){};
-#ifdef __NUANCE_ENABLED__
- if (nuance_event) delete nuance_event;
-#endif
-
-#ifdef __GiBUU_ENABLED__
- if (GiRead) delete GiRead;
-#endif
-};
-
-BaseFitEvt::BaseFitEvt(const BaseFitEvt* obj) {
+BaseFitEvt::BaseFitEvt(const BaseFitEvt *obj) {
Mode = obj->Mode;
probe_E = obj->probe_E;
probe_pdg = obj->probe_pdg;
Weight = obj->Weight;
InputWeight = obj->InputWeight;
RWWeight = obj->RWWeight;
CustomWeight = obj->CustomWeight;
SavedRWWeight = obj->SavedRWWeight;
fSplineCoeff = obj->fSplineCoeff;
fSplineRead = obj->fSplineRead;
fGenInfo = obj->fGenInfo;
fType = obj->fType;
#ifdef __NEUT_ENABLED__
fNeutVect = obj->fNeutVect;
#endif
#ifdef __NUWRO_ENABLED__
fNuwroEvent = obj->fNuwroEvent;
#endif
#ifdef __GENIE_ENABLED__
genie_event = obj->genie_event;
+#ifdef __DUNERWT_ENABLED__
+ HasDUNERwtPolyResponses = obj->HasDUNERwtPolyResponses;
+ DUNERwtPolyResponses = obj->DUNERwtPolyResponses;
+ HasDUNERwtResponses = obj->HasDUNERwtResponses;
+ DUNERwtResponses = obj->DUNERwtResponses;
+#endif
#endif
#ifdef __NUANCE_ENABLED__
nuance_event = obj->nuance_event;
#endif
#ifdef __GiBUU_ENABLED__
GiRead = obj->GiRead;
#endif
};
-BaseFitEvt::BaseFitEvt(BaseFitEvt const& other) {
+BaseFitEvt::BaseFitEvt(BaseFitEvt const &other) {
Mode = other.Mode;
probe_E = other.probe_E;
probe_pdg = other.probe_pdg;
Weight = other.Weight;
InputWeight = other.InputWeight;
RWWeight = other.RWWeight;
CustomWeight = other.CustomWeight;
SavedRWWeight = other.SavedRWWeight;
fSplineCoeff = other.fSplineCoeff;
fSplineRead = other.fSplineRead;
fGenInfo = other.fGenInfo;
fType = other.fType;
#ifdef __NEUT_ENABLED__
fNeutVect = other.fNeutVect;
#endif
#ifdef __NUWRO_ENABLED__
fNuwroEvent = other.fNuwroEvent;
#ifdef __USE_NUWRO_SRW_EVENTS__
- fNuwroSRWEvent = other.fNuwroSRWEvent; ///< Pointer to Nuwro event
+ fNuwroSRWEvent = other.fNuwroSRWEvent; ///< Pointer to Nuwro event
fNuwroParams = other.fNuwroParams;
#endif
#endif
#ifdef __GENIE_ENABLED__
genie_event = other.genie_event;
+#ifdef __DUNERWT_ENABLED__
+ HasDUNERwtPolyResponses = other.HasDUNERwtPolyResponses;
+ DUNERwtPolyResponses = other.DUNERwtPolyResponses;
+ HasDUNERwtResponses = other.HasDUNERwtResponses;
+ DUNERwtResponses = other.DUNERwtResponses;
+#endif
#endif
#ifdef __NUANCE_ENABLED__
nuance_event = other.nuance_event;
#endif
#ifdef __GiBUU_ENABLED__
GiRead = other.GiRead;
#endif
};
-BaseFitEvt BaseFitEvt::operator=(BaseFitEvt const& other) {
+BaseFitEvt BaseFitEvt::operator=(BaseFitEvt const &other) {
Mode = other.Mode;
probe_E = other.probe_E;
probe_pdg = other.probe_pdg;
Weight = other.Weight;
InputWeight = other.InputWeight;
RWWeight = other.RWWeight;
CustomWeight = other.CustomWeight;
SavedRWWeight = other.SavedRWWeight;
fSplineCoeff = other.fSplineCoeff;
fSplineRead = other.fSplineRead;
fGenInfo = other.fGenInfo;
fType = other.fType;
#ifdef __NEUT_ENABLED__
fNeutVect = other.fNeutVect;
#endif
#ifdef __NUWRO_ENABLED__
fNuwroEvent = other.fNuwroEvent;
#ifdef __USE_NUWRO_SRW_EVENTS__
- fNuwroSRWEvent = other.fNuwroSRWEvent; ///< Pointer to Nuwro event
+ fNuwroSRWEvent = other.fNuwroSRWEvent; ///< Pointer to Nuwro event
fNuwroParams = other.fNuwroParams;
#endif
#endif
#ifdef __GENIE_ENABLED__
genie_event = other.genie_event;
+#ifdef __DUNERWT_ENABLED__
+ HasDUNERwtPolyResponses = other.HasDUNERwtPolyResponses;
+ DUNERwtPolyResponses = other.DUNERwtPolyResponses;
+ HasDUNERwtResponses = other.HasDUNERwtResponses;
+ DUNERwtResponses = other.DUNERwtResponses;
+#endif
#endif
#ifdef __NUANCE_ENABLED__
nuance_event = other.nuance_event;
#endif
#ifdef __GiBUU_ENABLED__
GiRead = other.GiRead;
#endif
return *this;
}
void BaseFitEvt::ResetWeight() { InputWeight = 1.0; }
double BaseFitEvt::GetWeight() {
return InputWeight * RWWeight * CustomWeight;
};
#ifdef __NEUT_ENABLED__
-void BaseFitEvt::SetNeutVect(NeutVect* v) {
+void BaseFitEvt::SetNeutVect(NeutVect *v) {
fType = kNEUT;
fNeutVect = v;
}
#endif
#ifdef __GENIE_ENABLED__
-void BaseFitEvt::SetGenieEvent(NtpMCEventRecord* ntpl) {
+void BaseFitEvt::SetGenieEvent(NtpMCEventRecord *ntpl) {
fType = kGENIE;
genie_event = ntpl;
}
#endif
#ifdef __NUANCE_ENABLED__
-void BaseFitEvt::SetNuanceEvent(NuanceEvent* e) {
+void BaseFitEvt::SetNuanceEvent(NuanceEvent *e) {
fType = kNUANCE;
nuance_event = e;
}
#endif
void BaseFitEvt::SetInputFitEvent() { fType = kINPUTFITEVENT; }
void BaseFitEvt::SetInputFitSpline() { fType = kNEWSPLINE; }
void BaseFitEvt::SetInputHepMC() { fType = kHEPMC; }
diff --git a/src/InputHandler/BaseFitEvt.h b/src/InputHandler/BaseFitEvt.h
index 448c130..8f16111 100644
--- a/src/InputHandler/BaseFitEvt.h
+++ b/src/InputHandler/BaseFitEvt.h
@@ -1,139 +1,145 @@
// 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 .
-*******************************************************************************/
+ * 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 .
+ *******************************************************************************/
#ifndef FITEVENTBASE_H_SEEN
#define FITEVENTBASE_H_SEEN
/*!
* \addtogroup InputHandler
* @{
*/
#ifdef __NEUT_ENABLED__
#include "neutpart.h"
#include "neutvect.h"
#endif
#ifdef __NUWRO_ENABLED__
#ifdef __USE_NUWRO_SRW_EVENTS__
#include "NuwroReWeightSimpleEvent.h"
#else
#include "event1.h"
#endif
#endif
#ifdef __GENIE_ENABLED__
#include "EVGCore/EventRecord.h"
-#include "GHEP/GHepRecord.h"
#include "Ntuple/NtpMCEventRecord.h"
using namespace genie;
+#ifdef __DUNERWT_ENABLED__
+#include "systematicstools/interpreters/PrecalculatedResponseReader.hh"
+#endif
#endif
#ifdef __NUANCE_ENABLED__
#include "NuanceEvent.h"
#endif
-#include "StdHepEvt.h"
-#include "SplineReader.h"
-#include "InputTypes.h"
#include "GeneratorInfoBase.h"
+#include "InputTypes.h"
+#include "SplineReader.h"
+#include "StdHepEvt.h"
/// Base Event Class used to store just the generator event pointers
class BaseFitEvt {
- public:
-
+public:
/// Base Constructor
BaseFitEvt(void);
~BaseFitEvt();
/// Copy base fit event pointers
- BaseFitEvt(const BaseFitEvt* obj);
+ BaseFitEvt(const BaseFitEvt *obj);
BaseFitEvt(BaseFitEvt const &);
BaseFitEvt operator=(BaseFitEvt const &);
/// Reset weight to 1.0
void ResetWeight();
/// Return combined weight for this event
double GetWeight();
/// Manually set event type
- inline void SetType(int type){fType = type;};
+ inline void SetType(int type) { fType = type; };
// Global Event Variables/Weights
- int Mode; ///< True interaction mode
- double probe_E; ///< True probe energy
+ int Mode; ///< True interaction mode
+ double probe_E; ///< True probe energy
double probe_pdg;
// Weighting Info
double Weight; ///< Total Weight For Event
double InputWeight; ///< Input Starting Weight (used for GiBUU)
double RWWeight; ///< Saved RW from FitWeight
double CustomWeight; ///< Extra custom weight that samples can set
double SavedRWWeight; ///< Saved RW value for FitEvents
// Spline Info Coefficients and Readers
- float* fSplineCoeff; ///< ND Array of Spline Coefficients
- SplineReader* fSplineRead; ///< Spline Interpretter
+ float *fSplineCoeff; ///< ND Array of Spline Coefficients
+ SplineReader *fSplineRead; ///< Spline Interpretter
// Generator Info
- GeneratorInfoBase* fGenInfo; ///< Generator Variable Box
- UInt_t fType; ///< Generator Event Type
+ GeneratorInfoBase *fGenInfo; ///< Generator Variable Box
+ UInt_t fType; ///< Generator Event Type
#ifdef __NEUT_ENABLED__
/// Setup Event Reading from NEUT Event
- void SetNeutVect(NeutVect* v);
- NeutVect* fNeutVect; ///< Pointer to Neut Vector
+ void SetNeutVect(NeutVect *v);
+ NeutVect *fNeutVect; ///< Pointer to Neut Vector
#endif
#ifdef __NUWRO_ENABLED__
#ifdef __USE_NUWRO_SRW_EVENTS__
- SRW::SRWEvent fNuwroSRWEvent; ///< Pointer to Nuwro event
- params * fNuwroParams;
+ SRW::SRWEvent fNuwroSRWEvent; ///< Pointer to Nuwro event
+ params *fNuwroParams;
#endif
- event* fNuwroEvent; ///< Pointer to Nuwro event
+ event *fNuwroEvent; ///< Pointer to Nuwro event
#endif
#ifdef __GENIE_ENABLED__
/// Setup Event Reading from GENIE Event
- void SetGenieEvent(NtpMCEventRecord* ntpl);
- NtpMCEventRecord* genie_event; ///< Pointer to NTuple Genie Event
- GHepRecord* genie_record; ///< Pointer to actually accessible Genie Record
+ void SetGenieEvent(NtpMCEventRecord *ntpl);
+ NtpMCEventRecord const *genie_event; ///< Pointer to NTuple Genie Event
+#ifdef __DUNERWT_ENABLED__
+ bool HasDUNERwtPolyResponses;
+ std::vector::ParamPolyResponses>
+ DUNERwtPolyResponses;
+ bool HasDUNERwtResponses;
+ systtools::event_unit_response_t DUNERwtResponses;
+#endif
#endif
#ifdef __NUANCE_ENABLED__
/// Setup Event Reading from NUANCE Event
- void SetNuanceEvent(NuanceEvent* e);
- NuanceEvent* nuance_event; ///< Pointer to Nuance reader
+ void SetNuanceEvent(NuanceEvent *e);
+ NuanceEvent *nuance_event; ///< Pointer to Nuance reader
#endif
#ifdef __GiBUU_ENABLED__
- GiBUUStdHepReader* GiRead; ///< Pointer to GiBUU reader
+ GiBUUStdHepReader *GiRead; ///< Pointer to GiBUU reader
#endif
/// Setup Event Type to FitEvent
void SetInputFitEvent();
/// Setup Event Type to FitSpline
void SetInputFitSpline();
/// Setup Event Type to HepMC
void SetInputHepMC();
-
};
/*! @} */
#endif
diff --git a/src/InputHandler/FitEvent.cxx b/src/InputHandler/FitEvent.cxx
index 50ef1bc..7516fee 100644
--- a/src/InputHandler/FitEvent.cxx
+++ b/src/InputHandler/FitEvent.cxx
@@ -1,429 +1,466 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
-* This file is pddrt of NUISANCE.
-*
-* NUISANCE is free software: you can redistribute it and/or modify
-* it under the terms of the GNU General Public License as published by
-* the Free Software Foundation, either version 3 of the License, or
-* (at your option) any later version.
-*
-* NUISANCE is distributed in the hope that it will be useful,
-* but WITHOUT ANY WARRANTY; without even the implied warranty of
-* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-* GNU General Public License for more details.
-*
-* You should have received a copy of the GNU General Public License
-* along with NUISANCE. If not, see .
-*******************************************************************************/
+ * This file is pddrt of NUISANCE.
+ *
+ * NUISANCE is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * NUISANCE is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with NUISANCE. If not, see .
+ *******************************************************************************/
#include "FitEvent.h"
-#include
#include "TObjArray.h"
+#include
FitEvent::FitEvent() {
fGenInfo = NULL;
kRemoveFSIParticles = true;
kRemoveUndefParticles = true;
- AllocateParticleStack(400);
+ AllocateParticleStack(25);
};
-void FitEvent::AddGeneratorInfo(GeneratorInfoBase* gen) {
+void FitEvent::AddGeneratorInfo(GeneratorInfoBase *gen) {
fGenInfo = gen;
gen->AllocateParticleStack(kMaxParticles);
}
void FitEvent::AllocateParticleStack(int stacksize) {
- LOG(DEB) << "Allocating particle stack of size: " << stacksize << std::endl;
+
kMaxParticles = stacksize;
- fParticleList = new FitParticle*[kMaxParticles];
+ if (!stacksize) {
+ return;
+ }
+
+ LOG(DEB) << "Allocating particle stack of size: " << stacksize << std::endl;
- fParticleMom = new double*[kMaxParticles];
+ fParticleList = new FitParticle *[kMaxParticles];
+
+ fParticleMom = new double *[kMaxParticles];
fParticleState = new UInt_t[kMaxParticles];
fParticlePDG = new int[kMaxParticles];
- fOrigParticleMom = new double*[kMaxParticles];
+ fOrigParticleMom = new double *[kMaxParticles];
fOrigParticleState = new UInt_t[kMaxParticles];
fOrigParticlePDG = new int[kMaxParticles];
for (size_t i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
fParticleMom[i] = new double[4];
fOrigParticleMom[i] = new double[4];
}
- if (fGenInfo) fGenInfo->AllocateParticleStack(kMaxParticles);
+ if (fGenInfo)
+ fGenInfo->AllocateParticleStack(kMaxParticles);
}
void FitEvent::ExpandParticleStack(int stacksize) {
DeallocateParticleStack();
AllocateParticleStack(stacksize);
}
void FitEvent::DeallocateParticleStack() {
+
+ if (!kMaxParticles) {
+ return;
+ }
+
for (size_t i = 0; i < kMaxParticles; i++) {
- if (fParticleList[i]) delete fParticleList[i];
- delete fParticleMom[i];
- delete fOrigParticleMom[i];
+ if (fParticleList[i])
+ delete fParticleList[i];
+ delete[] fParticleMom[i];
+ delete[] fOrigParticleMom[i];
}
- delete fParticleMom;
- delete fOrigParticleMom;
+ delete[] fParticleMom;
+ delete[] fOrigParticleMom;
- delete fParticleList;
+ delete[] fParticleList;
- delete fParticleState;
- delete fParticlePDG;
+ delete[] fParticleState;
+ delete[] fParticlePDG;
- delete fOrigParticleState;
- delete fOrigParticlePDG;
+ delete[] fOrigParticleState;
+ delete[] fOrigParticlePDG;
- if (fGenInfo) fGenInfo->DeallocateParticleStack();
+ if (fGenInfo)
+ fGenInfo->DeallocateParticleStack();
kMaxParticles = 0;
+
+ fNParticles = 0;
}
void FitEvent::ClearFitParticles() {
for (size_t i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
}
}
void FitEvent::FreeFitParticles() {
for (size_t i = 0; i < kMaxParticles; i++) {
- FitParticle* fp = fParticleList[i];
- if (fp) delete fp;
+ FitParticle *fp = fParticleList[i];
+ if (fp)
+ delete fp;
fParticleList[i] = NULL;
}
}
void FitEvent::ResetParticleList() {
for (unsigned int i = 0; i < kMaxParticles; i++) {
- FitParticle* fp = fParticleList[i];
- if (fp) delete fp;
+ FitParticle *fp = fParticleList[i];
+ if (fp)
+ delete fp;
fParticleList[i] = NULL;
}
}
void FitEvent::HardReset() {
for (unsigned int i = 0; i < kMaxParticles; i++) {
fParticleList[i] = NULL;
}
}
void FitEvent::ResetEvent() {
Mode = 9999;
fEventNo = -1;
fTotCrs = -1.0;
fTargetA = -1;
fTargetZ = -1;
fTargetH = -1;
fBound = false;
fNParticles = 0;
- if (fGenInfo) fGenInfo->Reset();
+ if (fGenInfo)
+ fGenInfo->Reset();
for (unsigned int i = 0; i < kMaxParticles; i++) {
- if (fParticleList[i]) delete fParticleList[i];
+ if (fParticleList[i])
+ delete fParticleList[i];
fParticleList[i] = NULL;
continue;
fParticlePDG[i] = 0;
fParticleState[i] = kUndefinedState;
fParticleMom[i][0] = 0.0;
fParticleMom[i][1] = 0.0;
fParticleMom[i][2] = 0.0;
fParticleMom[i][3] = 0.0;
fOrigParticlePDG[i] = 0;
fOrigParticleState[i] = kUndefinedState;
fOrigParticleMom[i][0] = 0.0;
fOrigParticleMom[i][1] = 0.0;
fOrigParticleMom[i][2] = 0.0;
fOrigParticleMom[i][3] = 0.0;
}
}
void FitEvent::OrderStack() {
// Copy current stack
int npart = fNParticles;
for (int i = 0; i < npart; i++) {
fOrigParticlePDG[i] = fParticlePDG[i];
fOrigParticleState[i] = fParticleState[i];
fOrigParticleMom[i][0] = fParticleMom[i][0];
fOrigParticleMom[i][1] = fParticleMom[i][1];
fOrigParticleMom[i][2] = fParticleMom[i][2];
fOrigParticleMom[i][3] = fParticleMom[i][3];
}
// Now run loops for each particle
fNParticles = 0;
int stateorder[6] = {kInitialState, kFinalState, kFSIState,
kNuclearInitial, kNuclearRemnant, kUndefinedState};
for (int s = 0; s < 6; s++) {
for (int i = 0; i < npart; i++) {
- if ((UInt_t)fOrigParticleState[i] != (UInt_t)stateorder[s]) continue;
+ if ((UInt_t)fOrigParticleState[i] != (UInt_t)stateorder[s])
+ continue;
fParticlePDG[fNParticles] = fOrigParticlePDG[i];
fParticleState[fNParticles] = fOrigParticleState[i];
fParticleMom[fNParticles][0] = fOrigParticleMom[i][0];
fParticleMom[fNParticles][1] = fOrigParticleMom[i][1];
fParticleMom[fNParticles][2] = fOrigParticleMom[i][2];
fParticleMom[fNParticles][3] = fOrigParticleMom[i][3];
fNParticles++;
}
}
if (LOG_LEVEL(DEB)) {
LOG(DEB) << "Ordered stack" << std::endl;
for (int i = 0; i < fNParticles; i++) {
LOG(DEB) << "Particle " << i << ". " << fParticlePDG[i] << " "
<< fParticleMom[i][0] << " " << fParticleMom[i][1] << " "
<< fParticleMom[i][2] << " " << fParticleMom[i][3] << " "
<< fParticleState[i] << std::endl;
}
}
if (fNParticles != npart) {
ERR(FTL) << "Dropped some particles when ordering the stack!" << std::endl;
}
return;
}
void FitEvent::Print() {
if (LOG_LEVEL(FIT)) {
LOG(FIT) << "FITEvent print" << std::endl;
LOG(FIT) << "Mode: " << Mode << ", Weight: " << InputWeight << std::endl;
LOG(FIT) << "Particles: " << fNParticles << std::endl;
LOG(FIT) << " -> Particle Stack " << std::endl;
for (int i = 0; i < fNParticles; i++) {
LOG(FIT) << " -> -> " << i << ". " << fParticlePDG[i] << " "
<< fParticleState[i] << " "
<< " Mom(" << fParticleMom[i][0] << ", " << fParticleMom[i][1]
<< ", " << fParticleMom[i][2] << ", " << fParticleMom[i][3]
<< ")." << std::endl;
}
}
return;
}
/* Read/Write own event class */
-void FitEvent::SetBranchAddress(TChain* tn) {
+void FitEvent::SetBranchAddress(TChain *tn) {
tn->SetBranchAddress("Mode", &Mode);
tn->SetBranchAddress("EventNo", &fEventNo);
tn->SetBranchAddress("TotCrs", &fTotCrs);
tn->SetBranchAddress("TargetA", &fTargetA);
tn->SetBranchAddress("TargetH", &fTargetH);
tn->SetBranchAddress("Bound", &fBound);
tn->SetBranchAddress("RWWeight", &SavedRWWeight);
tn->SetBranchAddress("InputWeight", &InputWeight);
}
-void FitEvent::AddBranchesToTree(TTree* tn) {
+void FitEvent::AddBranchesToTree(TTree *tn) {
tn->Branch("Mode", &Mode, "Mode/I");
tn->Branch("EventNo", &fEventNo, "EventNo/i");
tn->Branch("TotCrs", &fTotCrs, "TotCrs/D");
tn->Branch("TargetA", &fTargetA, "TargetA/I");
tn->Branch("TargetH", &fTargetH, "TargetH/I");
tn->Branch("Bound", &fBound, "Bound/O");
tn->Branch("RWWeight", &RWWeight, "RWWeight/D");
tn->Branch("InputWeight", &InputWeight, "InputWeight/D");
tn->Branch("NParticles", &fNParticles, "NParticles/I");
tn->Branch("ParticleState", fOrigParticleState,
"ParticleState[NParticles]/i");
tn->Branch("ParticlePDG", fOrigParticlePDG, "ParticlePDG[NParticles]/I");
tn->Branch("ParticleMom", fOrigParticleMom, "ParticleMom[NParticles][4]/D");
}
// ------- EVENT ACCESS FUNCTION --------- //
TLorentzVector FitEvent::GetParticleP4(int index) const {
- if (index == -1 or index >= fNParticles) return TLorentzVector();
+ if (index == -1 or index >= fNParticles)
+ return TLorentzVector();
return TLorentzVector(fParticleMom[index][0], fParticleMom[index][1],
fParticleMom[index][2], fParticleMom[index][3]);
}
TVector3 FitEvent::GetParticleP3(int index) const {
- if (index == -1 or index >= fNParticles) return TVector3();
+ if (index == -1 or index >= fNParticles)
+ return TVector3();
return TVector3(fParticleMom[index][0], fParticleMom[index][1],
fParticleMom[index][2]);
}
double FitEvent::GetParticleMom(int index) const {
- if (index == -1 or index >= fNParticles) return 0.0;
+ if (index == -1 or index >= fNParticles)
+ return 0.0;
return sqrt(fParticleMom[index][0] * fParticleMom[index][0] +
fParticleMom[index][1] * fParticleMom[index][1] +
fParticleMom[index][2] * fParticleMom[index][2]);
}
double FitEvent::GetParticleMom2(int index) const {
- if (index == -1 or index >= fNParticles) return 0.0;
+ if (index == -1 or index >= fNParticles)
+ return 0.0;
return fabs((fParticleMom[index][0] * fParticleMom[index][0] +
fParticleMom[index][1] * fParticleMom[index][1] +
fParticleMom[index][2] * fParticleMom[index][2]));
}
double FitEvent::GetParticleE(int index) const {
- if (index == -1 or index >= fNParticles) return 0.0;
+ if (index == -1 or index >= fNParticles)
+ return 0.0;
return fParticleMom[index][3];
}
int FitEvent::GetParticleState(int index) const {
- if (index == -1 or index >= fNParticles) return kUndefinedState;
+ if (index == -1 or index >= fNParticles)
+ return kUndefinedState;
return (fParticleState[index]);
}
int FitEvent::GetParticlePDG(int index) const {
- if (index == -1 or index >= fNParticles) return 0;
+ if (index == -1 or index >= fNParticles)
+ return 0;
return (fParticlePDG[index]);
}
-FitParticle* FitEvent::GetParticle(int const i) {
+FitParticle *FitEvent::GetParticle(int const i) {
// Check Valid Index
if (i == -1) {
return NULL;
}
// Check Valid
if (i > fNParticles) {
ERR(FTL) << "Requesting particle beyond stack!" << std::endl
<< "i = " << i << " N = " << fNParticles << std::endl
<< "Mode = " << Mode << std::endl;
throw;
}
if (!fParticleList[i]) {
/*
std::cout << "Creating particle with values i " << i << " ";
std::cout << fParticleMom[i][0] << " " << fParticleMom[i][1] << " " <<
fParticleMom[i][2] << " " << fParticleMom[i][3] << " ";
std::cout << fParticlePDG[i] << " " << fParticleState[i] << std::endl;
*/
fParticleList[i] = new FitParticle(fParticleMom[i][0], fParticleMom[i][1],
fParticleMom[i][2], fParticleMom[i][3],
fParticlePDG[i], fParticleState[i]);
} else {
/*
std::cout << "Filling particle with values i " << i << " ";
std::cout << fParticleMom[i][0] << " " << fParticleMom[i][1] << " " <<
fParticleMom[i][2] << " " << fParticleMom[i][3] << " ";
std::cout << fParticlePDG[i] << " "<< fParticleState[i] <SetValues(fParticleMom[i][0], fParticleMom[i][1],
fParticleMom[i][2], fParticleMom[i][3],
fParticlePDG[i], fParticleState[i]);
}
return fParticleList[i];
}
bool FitEvent::HasParticle(int const pdg, int const state) const {
bool found = false;
for (int i = 0; i < fNParticles; i++) {
- if (state != -1 && fParticleState[i] != (uint)state) continue;
- if (fParticlePDG[i] == pdg) found = true;
+ if (state != -1 && fParticleState[i] != (uint)state)
+ continue;
+ if (fParticlePDG[i] == pdg)
+ found = true;
}
return found;
}
int FitEvent::NumParticle(int const pdg, int const state) const {
int nfound = 0;
for (int i = 0; i < fNParticles; i++) {
- if (state != -1 and fParticleState[i] != (uint)state) continue;
- if (pdg == 0 or fParticlePDG[i] == pdg) nfound += 1;
+ if (state != -1 and fParticleState[i] != (uint)state)
+ continue;
+ if (pdg == 0 or fParticlePDG[i] == pdg)
+ nfound += 1;
}
return nfound;
}
std::vector FitEvent::GetAllParticleIndices(int const pdg,
int const state) const {
std::vector indexlist;
for (int i = 0; i < fNParticles; i++) {
- if (state != -1 and fParticleState[i] != (uint)state) continue;
+ if (state != -1 and fParticleState[i] != (uint)state)
+ continue;
if (pdg == 0 or fParticlePDG[i] == pdg) {
indexlist.push_back(i);
}
}
return indexlist;
}
-std::vector FitEvent::GetAllParticle(int const pdg,
- int const state) {
+std::vector FitEvent::GetAllParticle(int const pdg,
+ int const state) {
std::vector indexlist = GetAllParticleIndices(pdg, state);
- std::vector plist;
+ std::vector plist;
for (std::vector::iterator iter = indexlist.begin();
iter != indexlist.end(); iter++) {
plist.push_back(GetParticle((*iter)));
}
return plist;
}
int FitEvent::GetHMParticleIndex(int const pdg, int const state) const {
double maxmom2 = -9999999.9;
int maxind = -1;
for (int i = 0; i < fNParticles; i++) {
- if (state != -1 and fParticleState[i] != (uint)state) continue;
+ if (state != -1 and fParticleState[i] != (uint)state)
+ continue;
if (pdg == 0 or fParticlePDG[i] == pdg) {
double newmom2 = GetParticleMom2(i);
if (newmom2 > maxmom2) {
maxind = i;
maxmom2 = newmom2;
}
}
}
return maxind;
}
int FitEvent::GetBeamNeutrinoIndex(void) const {
for (int i = 0; i < fNParticles; i++) {
- if (fParticleState[i] != kInitialState) continue;
+ if (fParticleState[i] != kInitialState)
+ continue;
int pdg = abs(fParticlePDG[i]);
if (pdg == 12 or pdg == 14 or pdg == 16) {
return i;
}
}
return 0;
}
int FitEvent::GetBeamElectronIndex(void) const {
return GetHMISParticleIndex(11);
}
int FitEvent::GetBeamPionIndex(void) const {
return GetHMISParticleIndex(PhysConst::pdg_pions);
}
int FitEvent::NumFSMesons() {
int nMesons = 0;
for (int i = 0; i < fNParticles; i++) {
- if (fParticleState[i] != kFinalState) continue;
+ if (fParticleState[i] != kFinalState)
+ continue;
if (abs(fParticlePDG[i]) >= 111 && abs(fParticlePDG[i]) <= 557)
nMesons += 1;
}
return nMesons;
}
int FitEvent::NumFSLeptons(void) const {
int nLeptons = 0;
for (int i = 0; i < fNParticles; i++) {
- if (fParticleState[i] != kFinalState) continue;
+ if (fParticleState[i] != kFinalState)
+ continue;
if (abs(fParticlePDG[i]) == 11 || abs(fParticlePDG[i]) == 13 ||
abs(fParticlePDG[i]) == 15)
nLeptons += 1;
}
return nLeptons;
}
diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx
index c1fb236..278295d 100644
--- a/src/InputHandler/GENIEInputHandler.cxx
+++ b/src/InputHandler/GENIEInputHandler.cxx
@@ -1,454 +1,540 @@
// 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 .
-*******************************************************************************/
+ * This file is part of NUISANCE.
+ *
+ * NUISANCE is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * NUISANCE is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with NUISANCE. If not, see .
+ *******************************************************************************/
#ifdef __GENIE_ENABLED__
#include "GENIEInputHandler.h"
#include "InputUtils.h"
+#ifdef __DUNERWT_ENABLED__
+
+#include "systematicstools/utility/ParameterAndProviderConfigurationUtility.hh"
+
+#include "fhiclcpp/make_ParameterSet.h"
+
+#endif
+
GENIEGeneratorInfo::~GENIEGeneratorInfo() { DeallocateParticleStack(); }
-void GENIEGeneratorInfo::AddBranchesToTree(TTree* tn) {
+void GENIEGeneratorInfo::AddBranchesToTree(TTree *tn) {
tn->Branch("GenieParticlePDGs", &fGenieParticlePDGs, "GenieParticlePDGs/I");
}
-void GENIEGeneratorInfo::SetBranchesFromTree(TTree* tn) {
+void GENIEGeneratorInfo::SetBranchesFromTree(TTree *tn) {
tn->SetBranchAddress("GenieParticlePDGs", &fGenieParticlePDGs);
}
void GENIEGeneratorInfo::AllocateParticleStack(int stacksize) {
fGenieParticlePDGs = new int[stacksize];
}
void GENIEGeneratorInfo::DeallocateParticleStack() {
delete fGenieParticlePDGs;
}
-void GENIEGeneratorInfo::FillGeneratorInfo(NtpMCEventRecord* ntpl) {
+void GENIEGeneratorInfo::FillGeneratorInfo(NtpMCEventRecord *ntpl) {
Reset();
// Check for GENIE Event
- if (!ntpl) return;
- if (!ntpl->event) return;
+ if (!ntpl)
+ return;
+ if (!ntpl->event)
+ return;
// Cast Event Record
- GHepRecord* ghep = static_cast(ntpl->event);
- if (!ghep) return;
+ GHepRecord *ghep = static_cast(ntpl->event);
+ if (!ghep)
+ return;
// Fill Particle Stack
- GHepParticle* p = 0;
+ GHepParticle *p = 0;
TObjArrayIter iter(ghep);
// Loop over all particles
int i = 0;
- while ((p = (dynamic_cast((iter).Next())))) {
- if (!p) continue;
+ while ((p = (dynamic_cast((iter).Next())))) {
+ if (!p)
+ continue;
// Get PDG
fGenieParticlePDGs[i] = p->Pdg();
i++;
}
}
void GENIEGeneratorInfo::Reset() {
for (int i = 0; i < kMaxParticles; i++) {
fGenieParticlePDGs[i] = 0;
}
}
-GENIEInputHandler::GENIEInputHandler(std::string const& handle,
- std::string const& rawinputs) {
+GENIEInputHandler::GENIEInputHandler(std::string const &handle,
+ std::string const &rawinputs) {
LOG(SAM) << "Creating GENIEInputHandler : " << handle << std::endl;
// Run a joint input handling
fName = handle;
// Setup the TChain
fGENIETree = new TChain("gtree");
fSaveExtra = FitPar::Config().GetParB("SaveExtraGenie");
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
- TFile* inp_file = new TFile(
+ TFile *inp_file = new TFile(
InputUtils::ExpandInputDirectories(inputs[inp_it]).c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
THROW("GENIE File IsZombie() at : '"
<< inputs[inp_it] << "'" << std::endl
<< "Check that your file paths are correct and the file exists!"
<< std::endl
<< "$ ls -lh " << inputs[inp_it]);
}
// Get Flux/Event hist
- TH1D* fluxhist = (TH1D*)inp_file->Get("nuisance_flux");
- TH1D* eventhist = (TH1D*)inp_file->Get("nuisance_events");
+ TH1D *fluxhist = (TH1D *)inp_file->Get("nuisance_flux");
+ TH1D *eventhist = (TH1D *)inp_file->Get("nuisance_events");
if (!fluxhist or !eventhist) {
ERROR(FTL, "Input File Contents: " << inputs[inp_it]);
inp_file->ls();
THROW("GENIE FILE doesn't contain flux/xsec info."
<< std::endl
<< "Try running the app PrepareGENIE first on :" << inputs[inp_it]
<< std::endl
<< "$ PrepareGENIE -h");
}
// Get N Events
- TTree* genietree = (TTree*)inp_file->Get("gtree");
+ TTree *genietree = (TTree *)inp_file->Get("gtree");
if (!genietree) {
ERROR(FTL, "gtree not located in GENIE file: " << inputs[inp_it]);
THROW("Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevents = genietree->GetEntries();
if (nevents <= 0) {
THROW("Trying to a TTree with "
<< nevents << " to TChain from : " << inputs[inp_it]);
}
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add To TChain
fGENIETree->AddFile(inputs[inp_it].c_str());
}
// Registor all our file inputs
SetupJointInputs();
// Assign to tree
fEventType = kGENIE;
fGenieNtpl = NULL;
fGENIETree->SetBranchAddress("gmcrec", &fGenieNtpl);
// Libraries should be seen but not heard...
StopTalking();
fGENIETree->GetEntry(0);
StartTalking();
+#ifndef __DUNERWT_ENABLED__
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetGenieEvent(fGenieNtpl);
if (fSaveExtra) {
fGenieInfo = new GENIEGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fGenieInfo);
}
fNUISANCEEvent->HardReset();
+#else
+ std::vector HandlerOpts = Config::QueryKeys("GENIEInputHandler");
+ fUseCache = HandlerOpts.size() && HandlerOpts.front().Has("UseCache") &&
+ HandlerOpts.front().GetB("UseCache");
+
+ DUNERwtCachedResponseReader = nullptr;
+ HaveCachedResponseReader = false;
+ if (fUseCache && (inputs.size() == 1)) {
+ std::vector DuneRwtCacheParams =
+ Config::QueryKeys("DUNERwtResponseCache");
+ for (nuiskey &key : DuneRwtCacheParams) {
+ if (key.Has("Input") && (key.GetS("Input") == inputs.front()) &&
+ key.Has("CacheFile") && key.Has("ParameterFHiCL")) {
+
+ fhicl::ParameterSet ps =
+ fhicl::make_ParameterSet(key.GetS("ParameterFHiCL"));
+
+ fhicl::ParameterSet syst_providers = ps.get(
+ "generated_systematic_provider_configuration");
+
+ systtools::param_header_map_t configuredParameterHeaders =
+ systtools::BuildParameterHeaders(syst_providers);
+
+ DUNERwtCachedResponseReader =
+ std::make_unique>(
+ key.GetS("CacheFile"), "resp_tree",
+ configuredParameterHeaders.size());
+ HaveCachedResponseReader = true;
+ break;
+ }
+ }
+ } else {
+ fNUISANCEEvent = new FitEvent();
+ fNUISANCEEvent->SetGenieEvent(fGenieNtpl);
+
+ if (fSaveExtra) {
+ fGenieInfo = new GENIEGeneratorInfo();
+ fNUISANCEEvent->AddGeneratorInfo(fGenieInfo);
+ }
+
+ fNUISANCEEvent->HardReset();
+ }
+#endif
};
GENIEInputHandler::~GENIEInputHandler() {
// if (fGenieGHep) delete fGenieGHep;
// if (fGenieNtpl) delete fGenieNtpl;
// if (fGENIETree) delete fGENIETree;
// if (fGenieInfo) delete fGenieInfo;
}
void GENIEInputHandler::CreateCache() {
if (fCacheSize > 0) {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 1);
fGENIETree->SetCacheSize(fCacheSize);
}
}
void GENIEInputHandler::RemoveCache() {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 0);
fGENIETree->SetCacheSize(0);
}
-FitEvent* GENIEInputHandler::GetNuisanceEvent(const UInt_t entry,
+FitEvent *GENIEInputHandler::GetNuisanceEvent(const UInt_t entry,
const bool lightweight) {
- if (entry >= (UInt_t)fNEvents) return NULL;
+ if (entry >= (UInt_t)fNEvents)
+ return NULL;
+
+ // Reduce memory pressure from the cache by clearing out the last entry each
+ // time.
+ if (entry && rwEvs[entry - 1].NParticles()) {
+ rwEvs[entry - 1].DeallocateParticleStack();
+ }
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fGENIETree->GetEntry(entry);
+#ifdef __DUNERWT_ENABLED__
+ if (entry >= rwEvs.size()) {
+ rwEvs.push_back(FitEvent());
+ if (HaveCachedResponseReader) {
+ rwEvs.back().DUNERwtPolyResponses =
+ DUNERwtCachedResponseReader->GetEventResponse(entry);
+ rwEvs.back().HasDUNERwtPolyResponses = true;
+ }
+ }
+ rwEvs[entry].SetGenieEvent(fGenieNtpl);
+
+ fNUISANCEEvent = &rwEvs[entry];
+#endif
+
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
else {
// Check for GENIE Event
- if (!fGenieNtpl) return NULL;
- if (!fGenieNtpl->event) return NULL;
+ if (!fGenieNtpl)
+ return NULL;
+ if (!fGenieNtpl->event)
+ return NULL;
// Cast Event Record
- fGenieGHep = static_cast(fGenieNtpl->event);
- if (!fGenieGHep) return NULL;
+ fGenieGHep = fGenieNtpl->event;
+ if (!fGenieGHep)
+ return NULL;
TObjArrayIter iter(fGenieGHep);
- genie::GHepParticle* p;
- while ((p = (dynamic_cast((iter).Next())))) {
+ genie::GHepParticle *p;
+ while ((p = (dynamic_cast((iter).Next())))) {
if (!p) {
continue;
}
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
if (state != genie::kIStInitialState) {
continue;
}
fNUISANCEEvent->probe_E = p->E() * 1.E3;
fNUISANCEEvent->probe_pdg = p->Pdg();
break;
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
return fNUISANCEEvent;
}
-int GENIEInputHandler::GetGENIEParticleStatus(genie::GHepParticle* p,
+int GENIEInputHandler::GetGENIEParticleStatus(genie::GHepParticle *p,
int mode) {
/*
kIStUndefined = -1,
kIStInitialState = 0, / generator-level initial state /
kIStStableFinalState = 1, / generator-level final state:
particles to be tracked by detector-level MC /
kIStIntermediateState = 2,
kIStDecayedState = 3,
kIStCorrelatedNucleon = 10,
kIStNucleonTarget = 11,
kIStDISPreFragmHadronicState = 12,
kIStPreDecayResonantState = 13,
- kIStHadronInTheNucleus = 14, / hadrons inside the nucleus: marked
- for hadron transport modules to act on /
+ kIStHadronInTheNucleus = 14, / hadrons inside the nucleus:
+ marked for hadron transport modules to act on /
kIStFinalStateNuclearRemnant = 15, / low energy nuclear fragments
entering the record collectively as a 'hadronic blob' pseudo-particle /
kIStNucleonClusterTarget = 16, // for composite nucleons before
phase space decay
*/
int state = kUndefinedState;
switch (p->Status()) {
- case genie::kIStNucleonTarget:
- case genie::kIStInitialState:
- case genie::kIStCorrelatedNucleon:
- case genie::kIStNucleonClusterTarget:
+ case genie::kIStNucleonTarget:
+ case genie::kIStInitialState:
+ case genie::kIStCorrelatedNucleon:
+ case genie::kIStNucleonClusterTarget:
+ state = kInitialState;
+ break;
+
+ case genie::kIStStableFinalState:
+ state = kFinalState;
+ break;
+
+ case genie::kIStHadronInTheNucleus:
+ if (abs(mode) == 2)
state = kInitialState;
- break;
-
- case genie::kIStStableFinalState:
- state = kFinalState;
- break;
-
- case genie::kIStHadronInTheNucleus:
- if (abs(mode) == 2)
- state = kInitialState;
- else
- state = kFSIState;
- break;
-
- case genie::kIStPreDecayResonantState:
- case genie::kIStDISPreFragmHadronicState:
- case genie::kIStIntermediateState:
+ else
state = kFSIState;
- break;
-
- case genie::kIStFinalStateNuclearRemnant:
- case genie::kIStUndefined:
- case genie::kIStDecayedState:
- default:
- break;
+ break;
+
+ case genie::kIStPreDecayResonantState:
+ case genie::kIStDISPreFragmHadronicState:
+ case genie::kIStIntermediateState:
+ state = kFSIState;
+ break;
+
+ case genie::kIStFinalStateNuclearRemnant:
+ case genie::kIStUndefined:
+ case genie::kIStDecayedState:
+ default:
+ break;
}
// Flag to remove nuclear part in genie
if (p->Pdg() > 1000000) {
if (state == kInitialState)
state = kNuclearInitial;
else if (state == kFinalState)
state = kNuclearRemnant;
}
return state;
}
#endif
#ifdef __GENIE_ENABLED__
-int GENIEInputHandler::ConvertGENIEReactionCode(GHepRecord* gheprec) {
+int GENIEInputHandler::ConvertGENIEReactionCode(GHepRecord *gheprec) {
// Electron Scattering
if (gheprec->Summary()->ProcInfo().IsEM()) {
if (gheprec->Summary()->InitState().ProbePdg() == 11) {
if (gheprec->Summary()->ProcInfo().IsQuasiElastic())
return 1;
else if (gheprec->Summary()->ProcInfo().IsMEC())
return 2;
else if (gheprec->Summary()->ProcInfo().IsResonant())
return 13;
else if (gheprec->Summary()->ProcInfo().IsDeepInelastic())
return 26;
else {
ERROR(WRN,
"Unknown GENIE Electron Scattering Mode!"
<< std::endl
<< "ScatteringTypeId = "
<< gheprec->Summary()->ProcInfo().ScatteringTypeId() << " "
<< "InteractionTypeId = "
<< gheprec->Summary()->ProcInfo().InteractionTypeId()
<< std::endl
<< genie::ScatteringType::AsString(
gheprec->Summary()->ProcInfo().ScatteringTypeId())
<< " "
<< genie::InteractionType::AsString(
gheprec->Summary()->ProcInfo().InteractionTypeId())
<< " " << gheprec->Summary()->ProcInfo().IsMEC());
return 0;
}
}
// Weak CC
} else if (gheprec->Summary()->ProcInfo().IsWeakCC()) {
// CC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 2;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -2;
// CC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
// Weak NC
} else if (gheprec->Summary()->ProcInfo().IsWeakNC()) {
// NC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 32;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -32;
// NC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
}
return 0;
}
void GENIEInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Check for GENIE Event
- if (!fGenieNtpl) return;
- if (!fGenieNtpl->event) return;
+ if (!fGenieNtpl)
+ return;
+ if (!fGenieNtpl->event)
+ return;
// Cast Event Record
- fGenieGHep = static_cast(fGenieNtpl->event);
- if (!fGenieGHep) return;
+ fGenieGHep = fGenieNtpl->event;
+ if (!fGenieGHep)
+ return;
// Convert GENIE Reaction Code
fNUISANCEEvent->Mode = ConvertGENIEReactionCode(fGenieGHep);
// Set Event Info
fNUISANCEEvent->fEventNo = 0.0;
fNUISANCEEvent->fTotCrs = fGenieGHep->XSec();
fNUISANCEEvent->fTargetA = 0.0;
fNUISANCEEvent->fTargetZ = 0.0;
fNUISANCEEvent->fTargetH = 0;
fNUISANCEEvent->fBound = 0.0;
fNUISANCEEvent->InputWeight =
- 1.0; //(1E+38 / genie::units::cm2) * fGenieGHep->XSec();
+ 1.0; //(1E+38 / genie::units::cm2) * fGenieGHep->XSec();
// Get N Particle Stack
unsigned int npart = fGenieGHep->GetEntries();
unsigned int kmax = fNUISANCEEvent->kMaxParticles;
if (npart > kmax) {
- ERR(WRN) << "GENIE has too many particles, expanding stack." << std::endl;
fNUISANCEEvent->ExpandParticleStack(npart);
}
// Fill Particle Stack
- GHepParticle* p = 0;
+ GHepParticle *p = 0;
TObjArrayIter iter(fGenieGHep);
fNUISANCEEvent->fNParticles = 0;
// Loop over all particles
- while ((p = (dynamic_cast((iter).Next())))) {
- if (!p) continue;
+ while ((p = (dynamic_cast((iter).Next())))) {
+ if (!p)
+ continue;
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
// Remove Undefined
- if (kRemoveUndefParticles && state == kUndefinedState) continue;
+ if (kRemoveUndefParticles && state == kUndefinedState)
+ continue;
// Remove FSI
- if (kRemoveFSIParticles && state == kFSIState) continue;
+ if (kRemoveFSIParticles && state == kFSIState)
+ continue;
if (kRemoveNuclearParticles &&
(state == kNuclearInitial || state == kNuclearRemnant))
continue;
// Fill Vectors
int curpart = fNUISANCEEvent->fNParticles;
fNUISANCEEvent->fParticleState[curpart] = state;
// Mom
fNUISANCEEvent->fParticleMom[curpart][0] = p->Px() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][1] = p->Py() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][2] = p->Pz() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][3] = p->E() * 1.E3;
// PDG
fNUISANCEEvent->fParticlePDG[curpart] = p->Pdg();
// Add to N particle count
fNUISANCEEvent->fNParticles++;
// Extra Check incase GENIE fails.
if ((UInt_t)fNUISANCEEvent->fNParticles == kmax) {
ERR(WRN) << "Number of GENIE Particles exceeds maximum!" << std::endl;
ERR(WRN) << "Extend kMax, or run without including FSI particles!"
<< std::endl;
break;
}
}
// Fill Extra Stack
- if (fSaveExtra) fGenieInfo->FillGeneratorInfo(fGenieNtpl);
+ if (fSaveExtra)
+ fGenieInfo->FillGeneratorInfo(fGenieNtpl);
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
- FitParticle* ISNeutralLepton =
+ FitParticle *ISNeutralLepton =
fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos);
if (ISNeutralLepton) {
fNUISANCEEvent->probe_E = ISNeutralLepton->E();
fNUISANCEEvent->probe_pdg = ISNeutralLepton->PDG();
}
return;
}
void GENIEInputHandler::Print() {}
#endif
diff --git a/src/InputHandler/GENIEInputHandler.h b/src/InputHandler/GENIEInputHandler.h
index 1cc6b69..e21cd74 100644
--- a/src/InputHandler/GENIEInputHandler.h
+++ b/src/InputHandler/GENIEInputHandler.h
@@ -1,108 +1,134 @@
// 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 .
-*******************************************************************************/
+ * 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 .
+ *******************************************************************************/
#ifndef GENIEINPUTHANDLER_H
#define GENIEINPUTHANDLER_H
/*!
* \addtogroup InputHandler
* @{
*/
#ifdef __GENIE_ENABLED__
#include "InputHandler.h"
#include "InputUtils.h"
#include "PlotUtils.h"
-#include "GHEP/GHepParticle.h"
-#include "PDG/PDGUtils.h"
-#include "GHEP/GHepUtils.h"
#include "Conventions/Units.h"
#include "EVGCore/EventRecord.h"
-#include "GHEP/GHepRecord.h"
+#include "GHEP/GHepParticle.h"
+#include "GHEP/GHepUtils.h"
#include "Ntuple/NtpMCEventRecord.h"
+#include "PDG/PDGUtils.h"
using namespace genie;
/// GENIE Generator Container to save extra particle status codes.
class GENIEGeneratorInfo : public GeneratorInfoBase {
public:
- GENIEGeneratorInfo() {};
- virtual ~GENIEGeneratorInfo();
+ GENIEGeneratorInfo(){};
+ virtual ~GENIEGeneratorInfo();
- /// Assigns information to branches
- void AddBranchesToTree(TTree* tn);
+ /// Assigns information to branches
+ void AddBranchesToTree(TTree *tn);
- /// Setup reading information from branches
- void SetBranchesFromTree(TTree* tn);
+ /// Setup reading information from branches
+ void SetBranchesFromTree(TTree *tn);
- /// Allocate any dynamic arrays for a new particle stack size
- void AllocateParticleStack(int stacksize);
+ /// Allocate any dynamic arrays for a new particle stack size
+ void AllocateParticleStack(int stacksize);
- /// Clear any dynamic arrays
- void DeallocateParticleStack();
+ /// Clear any dynamic arrays
+ void DeallocateParticleStack();
- /// Read extra genie information from the event
- void FillGeneratorInfo(NtpMCEventRecord* ntpl);
+ /// Read extra genie information from the event
+ void FillGeneratorInfo(NtpMCEventRecord *ntpl);
- /// Reset extra information to default/empty values
- void Reset();
+ /// Reset extra information to default/empty values
+ void Reset();
- int kMaxParticles; ///< Number of particles in stack
- int* fGenieParticlePDGs; ///< GENIE Particle PDGs (example)
+ int kMaxParticles; ///< Number of particles in stack
+ int *fGenieParticlePDGs; ///< GENIE Particle PDGs (example)
};
/// Main GENIE InputHandler
class GENIEInputHandler : public InputHandlerBase {
public:
+#ifdef __DUNERWT_ENABLED__
+ std::vector rwEvs;
+#endif
- /// Standard constructor given a name and input files
- GENIEInputHandler(std::string const& handle, std::string const& rawinputs);
- virtual ~GENIEInputHandler();
-
- /// Create a TTree Cache to speed up file read
- void CreateCache();
+ /// Standard constructor given a name and input files
+ GENIEInputHandler(std::string const &handle, std::string const &rawinputs);
+ virtual ~GENIEInputHandler();
+
+ /// Create a TTree Cache to speed up file read
+ void CreateCache();
+
+ /// Remove TTree Cache to save memory
+ void RemoveCache();
+
+ /// Returns a NUISANCE format event from the GENIE TTree. If !lightweight
+ /// then CalcNUISANCEKinematics() is called to convert the GENIE event into
+ /// a standard NUISANCE format.
+ FitEvent *GetNuisanceEvent(const UInt_t entry,
+ const bool lightweight = false);
+
+#ifdef __DUNERWT_ENABLED__
+ // Returns filled BaseFitEvent for a given entry;
+ BaseFitEvt *GetBaseEvent(const UInt_t entry) {
+ if (rwEvs.size() <= entry) {
+ THROW("Tried to get cached BaseFitEv[" << entry << "], but only have "
+ << rwEvs.size()
+ << " in the cache.");
+ }
+ // Here the GENIE record is no longer valid.
+ if (rwEvs[entry].genie_event) {
+ rwEvs[entry].genie_event = 0;
+ }
+ return &rwEvs[entry];
+ }
+#endif
- /// Remove TTree Cache to save memory
- void RemoveCache();
-
- /// Returns a NUISANCE format event from the GENIE TTree. If !lightweight
- /// then CalcNUISANCEKinematics() is called to convert the GENIE event into
- /// a standard NUISANCE format.
- FitEvent* GetNuisanceEvent(const UInt_t entry, const bool lightweight = false);
+ /// Converts GENIE event into standard NUISANCE FitEvent by looping over all
+ /// particles in the event and adding them to stack in fNUISANCEEvent.
+ void CalcNUISANCEKinematics();
- /// Converts GENIE event into standard NUISANCE FitEvent by looping over all
- /// particles in the event and adding them to stack in fNUISANCEEvent.
- void CalcNUISANCEKinematics();
+ /// Placeholder for GENIE related event printing.
+ void Print();
- /// Placeholder for GENIE related event printing.
- void Print();
+ /// Converts GENIE particle status codes into NUISANCE status codes.
+ int GetGENIEParticleStatus(genie::GHepParticle *part, int mode = 0);
- /// Converts GENIE particle status codes into NUISANCE status codes.
- int GetGENIEParticleStatus(genie::GHepParticle* part, int mode = 0);
+ /// Converts GENIE event reaction codes into NUISANCE reaction codes.
+ int ConvertGENIEReactionCode(GHepRecord *gheprec);
- /// Converts GENIE event reaction codes into NUISANCE reaction codes.
- int ConvertGENIEReactionCode(GHepRecord* gheprec);
+ EventRecord *fGenieGHep; ///< Pointer to actual event record
+ NtpMCEventRecord *fGenieNtpl; ///< Ntpl Wrapper Class
- GHepRecord* fGenieGHep; ///< Pointer to actual event record
- NtpMCEventRecord* fGenieNtpl; ///< Ntpl Wrapper Class
+ TChain *fGENIETree; ///< Main GENIE Event TTree
+ bool fSaveExtra; ///< Flag to save Extra GENIE info into Nuisance Event
+ GENIEGeneratorInfo *fGenieInfo; ///< Extra GENIE Generator Info Writer
- TChain* fGENIETree; ///< Main GENIE Event TTree
- bool fSaveExtra; ///< Flag to save Extra GENIE info into Nuisance Event
- GENIEGeneratorInfo* fGenieInfo; ///< Extra GENIE Generator Info Writer
+#ifdef __DUNERWT_ENABLED__
+ std::unique_ptr>
+ DUNERwtCachedResponseReader;
+ bool HaveCachedResponseReader;
+ bool fUseCache;
+#endif
};
/*! @} */
#endif
-#endif
\ No newline at end of file
+#endif
diff --git a/src/InputHandler/InputHandler.h b/src/InputHandler/InputHandler.h
index 1110091..fa22421 100644
--- a/src/InputHandler/InputHandler.h
+++ b/src/InputHandler/InputHandler.h
@@ -1,145 +1,144 @@
// 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 .
*******************************************************************************/
#ifndef INPUTHANDLER2_H
#define INPUTHANDLER2_H
/*!
* \addtogroup InputHandler
* @{
*/
#include "TH1D.h"
#include "FitEvent.h"
#include "BaseFitEvt.h"
#include "TTreePerfStats.h"
/// Base InputHandler class defining how events are requested and setup.
class InputHandlerBase {
public:
/// Base constructor resets everything to default
InputHandlerBase();
/// Removes flux/event rate histograms
virtual ~InputHandlerBase();
/// Return NUISANCE FitEvent Class from given event entry.
/// Must be overriden by GeneratorInputHandler. Lightweight allows a faster option
/// to be given where only RW information is needed.
virtual FitEvent* GetNuisanceEvent(const UInt_t entry, const bool lightweight=false) = 0;
/// Calls GetNuisanceEvent(entry, TRUE);
virtual BaseFitEvt* GetBaseEvent(const UInt_t entry);
/// Print current event information
virtual void Print();
-
/// Return handler ID
inline std::string GetName (void) {return fName; };
/// Return Handler Event Type Index
inline int GetType (void) {return fEventType;};
/// Get Total Number of Events being Handled
inline virtual int GetNEvents (void) {return fNEvents; };
/// Get the Total Flux Histogram these events were generated with
inline virtual TH1D* GetFluxHistogram (void) {return fFluxHist; };
/// Get the Total Event Histogram these events were generated with
inline virtual TH1D* GetEventHistogram (void) {return fEventHist;};
/// Get the Total Cross-section Histogram (EventHist/FluxHist)
virtual TH1D* GetXSecHistogram(void);
/// Return all Flux Histograms for all InputFiles.
virtual std::vector GetFluxList(void);
/// Return all Event Histograms for all InputFiles
virtual std::vector GetEventList(void);
/// Return all Xsec Histograms for all InputFiles
virtual std::vector GetXSecList(void);
/// Placeholder to create a cache to speed up reads in GeneratorInputHandler
inline virtual void CreateCache(){};
/// Placeholder to remove optional cache to free up memory
inline virtual void RemoveCache(){};
/// Return starting NUISANCE event pointer (entry=0)
FitEvent* FirstNuisanceEvent();
/// Iterate to next NUISANCE event. Returns NULL when entry > fNEvents.
FitEvent* NextNuisanceEvent();
/// Returns starting Base Event Pointer (entry=0)
BaseFitEvt* FirstBaseEvent();
/// Iterate to next NUISANCE Base Event. Returns NULL when entry > fNEvents.
BaseFitEvt* NextBaseEvent();
/// Register an input file and update event/flux information
virtual void RegisterJointInput(std::string input, int n, TH1D* f, TH1D* e);
/// Finalise setup of Input event/flux information and calculate
/// joint input weights if joint input is provided.
virtual void SetupJointInputs();
/// Calculate a weight for the event given the joint input information.
/// Used to scale the relative proportion of multiple inputs correctly
/// with respect to one another.
virtual double GetInputWeight(int entry);
/// Returns the total predicted event rate for this input given the
/// low and high energy ranges. intOpt specifies the option the ROOT
/// TH1D integral should use. e.g. "" or "width"
double PredictedEventRate(double low = -9999.9, double high = -9999.9,
std::string intOpt = "");
/// Returns the total generated flux for this input given the
/// low and high energy ranges. intOpt specifies the option the ROOT
/// TH1D integral should use. e.g. "" or "width"
double TotalIntegratedFlux(double low = -9999.9, double high = -9999.9,
std::string intOpt = "");
/// Actual data members.
std::vector jointfluxinputs;
std::vector jointeventinputs;
std::vector jointindexlow;
std::vector jointindexhigh;
std::vector jointindexallowed;
size_t jointindexswitch;
bool jointinput;
std::vector jointindexscale;
std::string fName;
TH1D* fFluxHist;
TH1D* fEventHist;
TH1D* fXSecHist;
int fNEvents;
int fMaxEvents;
FitEvent* fNUISANCEEvent;
BaseFitEvt* fBaseEvent;
int fEventType;
int fCurrentIndex;
int fCacheSize;
bool kRemoveUndefParticles;
bool kRemoveFSIParticles;
bool kRemoveNuclearParticles;
TTreePerfStats* fTTreePerformance;
};
/*! @} */
#endif
diff --git a/src/Reweight/CMakeLists.txt b/src/Reweight/CMakeLists.txt
index 55b9fef..b5e82f5 100644
--- a/src/Reweight/CMakeLists.txt
+++ b/src/Reweight/CMakeLists.txt
@@ -1,83 +1,84 @@
# 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
-weightRPA.h
+DUNERwtWeightEngine.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
weightRPA.h
+DUNERwtWeightEngine.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/DUNERwtWeightEngine.cxx b/src/Reweight/DUNERwtWeightEngine.cxx
new file mode 100644
index 0000000..3477010
--- /dev/null
+++ b/src/Reweight/DUNERwtWeightEngine.cxx
@@ -0,0 +1,167 @@
+// 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 .
+ *******************************************************************************/
+
+#include "DUNERwtWeightEngine.h"
+
+#include
+#include
+
+DUNERwtWeightEngine::DUNERwtWeightEngine() { Config(); }
+
+void DUNERwtWeightEngine::Config() {
+ std::vector DuneRwtParam = Config::QueryKeys("DUNERwt");
+
+ if (DuneRwtParam.size() < 1) {
+ ERROR(WRN, "Oscillation parameters specified but no OscParam element "
+ "configuring the experimental characteristics found.\nExpect at "
+ "least . Pausing for "
+ "10...");
+ sleep(10);
+ return;
+ }
+
+ std::string fhicl_name = DuneRwtParam.front().GetS("ConfigFHiCL");
+
+ DUNErwt.LoadConfiguration(fhicl_name);
+}
+
+int DUNERwtWeightEngine::ConvDial(std::string name) {
+ if (!DUNErwt.HaveHeader(name)) {
+ THROW("DUNERwtWeightEngine passed dial: "
+ << name << " that it does not understand.");
+ }
+ return DUNErwt.GetHeaderId(name);
+}
+
+void DUNERwtWeightEngine::IncludeDial(std::string name, double startval) {
+ EnabledParams.push_back({systtools::paramId_t(ConvDial(name)), startval});
+}
+
+void DUNERwtWeightEngine::SetDialValue(int nuisenum, double val) {
+
+ systtools::paramId_t DuneRwtEnum = (nuisenum % 1000);
+ systtools::ParamValue &pval =
+ GetParamElementFromContainer(EnabledParams, DuneRwtEnum);
+ fHasChanged = (pval.val - val) > std::numeric_limits::epsilon();
+ pval.val = val;
+}
+void DUNERwtWeightEngine::SetDialValue(std::string name, double val) {
+ if (!IsDialIncluded(name)) {
+ THROW("DUNERwtWeightEngine passed dial: " << name
+ << " that is not enabled.");
+ }
+
+ systtools::ParamValue &pval =
+ GetParamElementFromContainer(EnabledParams, ConvDial(name));
+ fHasChanged = (pval.val - val) > std::numeric_limits::epsilon();
+ pval.val = val;
+}
+
+bool DUNERwtWeightEngine::IsDialIncluded(std::string name) {
+ return IsDialIncluded(ConvDial(name));
+}
+bool DUNERwtWeightEngine::IsDialIncluded(int nuisenum) {
+ systtools::paramId_t DuneRwtEnum = (nuisenum % 1000);
+ return systtools::ContainterHasParam(EnabledParams, DuneRwtEnum);
+}
+
+double DUNERwtWeightEngine::GetDialValue(std::string name) {
+ if (!IsDialIncluded(name)) {
+ THROW("DUNERwtWeightEngine passed dial: " << name
+ << " that is not enabled.");
+ }
+ systtools::ParamValue &pval =
+ GetParamElementFromContainer(EnabledParams, ConvDial(name));
+ return pval.val;
+}
+double DUNERwtWeightEngine::GetDialValue(int nuisenum) {
+ if (!IsDialIncluded(nuisenum)) {
+ THROW("DUNERwtWeightEngine passed dial: " << nuisenum
+ << " that is not enabled.");
+ }
+ systtools::paramId_t DuneRwtEnum = (nuisenum % 1000);
+ systtools::ParamValue &pval =
+ GetParamElementFromContainer(EnabledParams, DuneRwtEnum);
+ return pval.val;
+}
+
+void DUNERwtWeightEngine::Reconfigure(bool silent) { fHasChanged = false; };
+
+bool DUNERwtWeightEngine::NeedsEventReWeight() {
+ if (fHasChanged) {
+ return true;
+ }
+ return false;
+}
+
+double DUNERwtWeightEngine::CalcWeight(BaseFitEvt *evt) {
+
+ double weight = 1;
+
+ if (evt->HasDUNERwtPolyResponses) {
+ for (size_t i = 0; i < EnabledParams.size(); ++i) {
+ if (DUNErwt.IsSplineParam(EnabledParams[i].pid)) {
+ if (!ContainterHasParam(evt->DUNERwtPolyResponses,
+ EnabledParams[i].pid)) {
+ continue;
+ }
+ weight *= GetParamElementFromContainer(evt->DUNERwtPolyResponses,
+ EnabledParams[i].pid)
+ .resp.eval(EnabledParams[i].val);
+
+ evt->DUNERwtResponses =
+ DUNErwt.GetEventResponses(*evt->genie_event->event);
+
+ } else {
+ if (!evt->HasDUNERwtResponses) {
+ evt->DUNERwtResponses =
+ DUNErwt.GetEventResponses(*evt->genie_event->event);
+ evt->HasDUNERwtResponses = true;
+ }
+ weight *= DUNErwt.GetDiscreteResponse(EnabledParams[i].pid,
+ int(EnabledParams[i].val),
+ evt->DUNERwtResponses);
+ }
+ }
+ } else {
+ if (!evt->HasDUNERwtResponses) {
+ evt->DUNERwtResponses =
+ DUNErwt.GetEventResponses(*evt->genie_event->event);
+ evt->HasDUNERwtResponses = true;
+ }
+
+ for (size_t i = 0; i < EnabledParams.size(); ++i) {
+ if (DUNErwt.IsSplineParam(EnabledParams[i].pid)) {
+ weight *= DUNErwt.GetParameterResponse(
+ EnabledParams[i].pid, EnabledParams[i].val, evt->DUNERwtResponses);
+
+ } else {
+ weight *= DUNErwt.GetDiscreteResponse(EnabledParams[i].pid,
+ int(EnabledParams[i].val),
+ evt->DUNERwtResponses);
+ }
+ }
+ }
+
+ return weight;
+}
+
+void DUNERwtWeightEngine::Print() {
+ std::cout << "DUNERwtWeightEngine: " << std::endl;
+}
diff --git a/src/Reweight/DUNERwtWeightEngine.h b/src/Reweight/DUNERwtWeightEngine.h
new file mode 100644
index 0000000..908f1d3
--- /dev/null
+++ b/src/Reweight/DUNERwtWeightEngine.h
@@ -0,0 +1,64 @@
+#ifndef DUNERWTWEIGHTENGINE_SEEN
+#define DUNERWTWEIGHTENGINE_SEEN
+// 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 .
+*******************************************************************************/
+
+#include "WeightEngineBase.h"
+
+#include "systematicstools/interface/types.hh"
+
+#include "nusystematics/artless/response_helper.hh"
+
+#include
+
+class DUNERwtWeightEngine : public WeightEngineBase {
+
+ public:
+ DUNERwtWeightEngine();
+
+ nusyst::response_helper DUNErwt;
+
+ systtools::param_value_list_t EnabledParams;
+
+ void Config();
+
+ int ConvDial(std::string name);
+
+ // Functions requiring Override
+ void IncludeDial(std::string name, double startval);
+
+ void SetDialValue(int nuisenum, double val);
+ void SetDialValue(std::string name, double val);
+
+ bool IsDialIncluded(std::string name);
+ bool IsDialIncluded(int nuisenum);
+
+ double GetDialValue(std::string name);
+ double GetDialValue(int nuisenum);
+
+ void Reconfigure(bool silent);
+
+ bool NeedsEventReWeight();
+
+ double CalcWeight(BaseFitEvt* evt);
+
+ void Print();
+};
+
+#endif
diff --git a/src/Reweight/FitWeight.cxx b/src/Reweight/FitWeight.cxx
index 3559a3c..436f333 100644
--- a/src/Reweight/FitWeight.cxx
+++ b/src/Reweight/FitWeight.cxx
@@ -1,274 +1,293 @@
#include "FitWeight.h"
+#include "DUNERwtWeightEngine.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"
void FitWeight::AddRWEngine(int 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;
- default:
- THROW("CANNOT ADD RW Engine for unknown dial type: " << type);
- break;
+ 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;
+ case kDUNERwt:
+ fAllRW[type] = new DUNERwtWeightEngine();
+ break;
+ default:
+ THROW("CANNOT ADD RW Engine for unknown dial type: " << type);
+ break;
}
}
-WeightEngineBase* FitWeight::GetRWEngine(int type) {
+WeightEngineBase *FitWeight::GetRWEngine(int type) {
if (HasRWEngine(type)) {
return fAllRW[type];
}
THROW("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: {
- return fAllRW.count(type);
- }
- default: { THROW("CANNOT get RW Engine for dial type: " << type); }
+ case kNEUT:
+ case kNUWRO:
+ case kGENIE:
+ case kNORM:
+ case kLIKEWEIGHT:
+ case kT2K:
+ case kCUSTOM:
+ case kSPLINEPARAMETER:
+ case kNIWG:
+ case kOSCILLATION:
+ case kDUNERwt: {
+ return fAllRW.count(type);
+ }
+ default: { THROW("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) {
+
+ // Setup RW Engine Pointer
+ if (fAllRW.find(dialtype) == fAllRW.end()) {
+ AddRWEngine(dialtype);
+ }
+
// Get the dial enum
- int nuisenum = Reweight::ConvDial(name, dialtype);
+ int nuisenum;
+ if (dialtype != kDUNERwt) {
+ nuisenum = Reweight::ConvDial(name, dialtype);
+ } else {
+#ifdef __DUNERWT_ENABLED__
+ DUNERwtWeightEngine *drw =
+ static_cast(fAllRW[kDUNERwt]);
+ nuisenum = drw->ConvDial(name) + dialtype*1000;
+#else
+ THROW("[ERROR]: Not built with DUNEReWeight support but ")
+#endif
+ }
if (nuisenum == -1) {
THROW("NUISENUM == " << nuisenum << " for " << name);
}
- // Setup RW Engine Pointer
- if (fAllRW.find(dialtype) == fAllRW.end()) {
- AddRWEngine(dialtype);
- }
- WeightEngineBase* rw = fAllRW[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();
+ 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()) {
THROW("Cannot find RW Engine for dialtype = "
<< dialtype << ", " << Reweight::RemoveDialType(nuisenum));
}
// 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) {
+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;
}
}
ERR(FTL) << "No Dial Found! " << std::endl;
throw;
return -1;
}
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 FitWeight::CalcWeight(BaseFitEvt *evt) {
double rwweight = 1.0;
- for (std::map::iterator iter = fAllRW.begin();
+ 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) {
+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) {
+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);
// // 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;
+ 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() {
LOG(REC) << "Fit Weight State: " << std::endl;
for (size_t i = 0; i < fNameList.size(); i++) {
LOG(REC) << " -> Par " << i << ". " << fNameList[i] << " " << fValueList[i]
<< std::endl;
}
}
diff --git a/src/Reweight/WeightUtils.cxx b/src/Reweight/WeightUtils.cxx
index 1a1bd77..905ea93 100644
--- a/src/Reweight/WeightUtils.cxx
+++ b/src/Reweight/WeightUtils.cxx
@@ -1,614 +1,631 @@
#include "WeightUtils.h"
#include "FitLogger.h"
#ifdef __T2KREW_ENABLED__
#include "T2KGenieReWeight.h"
#include "T2KNIWGReWeight.h"
#include "T2KNIWGUtils.h"
#include "T2KNeutReWeight.h"
#include "T2KNeutUtils.h"
#include "T2KReWeight.h"
using namespace t2krew;
#endif
#ifdef __NIWG_ENABLED__
#include "NIWGReWeight.h"
#include "NIWGReWeight1piAngle.h"
#include "NIWGReWeight2010a.h"
#include "NIWGReWeight2012a.h"
#include "NIWGReWeight2014a.h"
#include "NIWGReWeightDeltaMass.h"
#include "NIWGReWeightEffectiveRPA.h"
#include "NIWGReWeightHadronMultSwitch.h"
#include "NIWGReWeightMEC.h"
#include "NIWGReWeightPiMult.h"
#include "NIWGReWeightProtonFSIbug.h"
#include "NIWGReWeightRPA.h"
#include "NIWGReWeightSpectralFunc.h"
#include "NIWGReWeightSplineEnu.h"
#include "NIWGSyst.h"
#include "NIWGSystUncertainty.h"
#endif
#ifdef __NEUT_ENABLED__
#include "NReWeight.h"
#include "NReWeightCasc.h"
#include "NReWeightNuXSecCCQE.h"
#include "NReWeightNuXSecCCRES.h"
#include "NReWeightNuXSecCOH.h"
#include "NReWeightNuXSecDIS.h"
#include "NReWeightNuXSecNC.h"
#include "NReWeightNuXSecNCEL.h"
#include "NReWeightNuXSecNCRES.h"
#include "NReWeightNuXSecRES.h"
#include "NReWeightNuclPiless.h"
#include "NSyst.h"
#include "NSystUncertainty.h"
#include "neutpart.h"
#include "neutvect.h"
#endif
#ifdef __NUWRO_ENABLED__
#include "event1.h"
#endif
#ifdef __NUWRO_REWEIGHT_ENABLED__
#include "NuwroReWeight.h"
#include "NuwroReWeight_FlagNorm.h"
#include "NuwroReWeight_QEL.h"
#include "NuwroReWeight_SPP.h"
#include "NuwroSyst.h"
#include "NuwroSystUncertainty.h"
#endif
#ifdef __GENIE_ENABLED__
#include "EVGCore/EventRecord.h"
-#include "EVGCore/EventRecord.h"
#include "GHEP/GHepRecord.h"
#include "GSyst.h"
#include "GSystUncertainty.h"
#include "Ntuple/NtpMCEventRecord.h"
#include "ReWeight/GReWeight.h"
#include "ReWeight/GReWeightAGKY.h"
#include "ReWeight/GReWeightDISNuclMod.h"
#include "ReWeight/GReWeightFGM.h"
#include "ReWeight/GReWeightFZone.h"
#include "ReWeight/GReWeightINuke.h"
#include "ReWeight/GReWeightNonResonanceBkg.h"
#include "ReWeight/GReWeightNuXSecCCQE.h"
#include "ReWeight/GReWeightNuXSecCCQEvec.h"
#include "ReWeight/GReWeightNuXSecCCRES.h"
#include "ReWeight/GReWeightNuXSecCOH.h"
#include "ReWeight/GReWeightNuXSecDIS.h"
#include "ReWeight/GReWeightNuXSecNC.h"
#include "ReWeight/GReWeightNuXSecNCEL.h"
#include "ReWeight/GReWeightNuXSecNCRES.h"
#include "ReWeight/GReWeightResonanceDecay.h"
using namespace genie;
using namespace genie::rew;
#endif
#include "GlobalDialList.h"
#include "ModeNormEngine.h"
#include "NUISANCESyst.h"
#include "OscWeightEngine.h"
//********************************************************************
TF1 FitBase::GetRWConvFunction(std::string const &type,
std::string const &name) {
//********************************************************************
std::string dialfunc = "x";
std::string parType = type;
double low = -10000.0;
double high = 10000.0;
- if (parType.find("parameter") == std::string::npos) parType += "_parameter";
+ if (parType.find("parameter") == std::string::npos)
+ parType += "_parameter";
std::string line;
std::ifstream card(
(GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card")
.c_str(),
std::ifstream::in);
while (std::getline(card >> std::ws, line, '\n')) {
std::vector inputlist = GeneralUtils::ParseToStr(line, " ");
// Check the line length
- if (inputlist.size() < 4) continue;
+ if (inputlist.size() < 4)
+ continue;
// Check whether this is a comment
- if (inputlist[0].c_str()[0] == '#') continue;
+ if (inputlist[0].c_str()[0] == '#')
+ continue;
// Check whether this is the correct parameter type
- if (inputlist[0].compare(parType) != 0) continue;
+ if (inputlist[0].compare(parType) != 0)
+ continue;
// Check the parameter name
- if (inputlist[1].compare(name) != 0) continue;
+ if (inputlist[1].compare(name) != 0)
+ continue;
// inputlist[2] should be the units... ignore for now
dialfunc = inputlist[3];
// High and low are optional, check whether they exist
- if (inputlist.size() > 4) low = GeneralUtils::StrToDbl(inputlist[4]);
- if (inputlist.size() > 5) high = GeneralUtils::StrToDbl(inputlist[5]);
+ if (inputlist.size() > 4)
+ low = GeneralUtils::StrToDbl(inputlist[4]);
+ if (inputlist.size() > 5)
+ high = GeneralUtils::StrToDbl(inputlist[5]);
}
TF1 convfunc = TF1((name + "_convfunc").c_str(), dialfunc.c_str(), low, high);
return convfunc;
}
//********************************************************************
std::string FitBase::GetRWUnits(std::string const &type,
std::string const &name) {
//********************************************************************
std::string unit = "sig.";
std::string parType = type;
if (parType.find("parameter") == std::string::npos) {
parType += "_parameter";
}
std::string line;
std::ifstream card(
(GeneralUtils::GetTopLevelDir() + "/parameters/dial_conversion.card")
.c_str(),
std::ifstream::in);
while (std::getline(card >> std::ws, line, '\n')) {
std::vector inputlist = GeneralUtils::ParseToStr(line, " ");
// Check the line length
- if (inputlist.size() < 3) continue;
+ if (inputlist.size() < 3)
+ continue;
// Check whether this is a comment
- if (inputlist[0].c_str()[0] == '#') continue;
+ if (inputlist[0].c_str()[0] == '#')
+ continue;
// Check whether this is the correct parameter type
- if (inputlist[0].compare(parType) != 0) continue;
+ if (inputlist[0].compare(parType) != 0)
+ continue;
// Check the parameter name
- if (inputlist[1].compare(name) != 0) continue;
+ if (inputlist[1].compare(name) != 0)
+ continue;
unit = inputlist[2];
break;
}
return unit;
}
//********************************************************************
double FitBase::RWAbsToSigma(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.GetX(val);
- if (fabs(conv_val) < 1E-10) conv_val = 0.0;
+ if (fabs(conv_val) < 1E-10)
+ conv_val = 0.0;
QLOG(FIT, "AbsToSigma(" << name << ") = " << val << " -> " << conv_val);
return conv_val;
}
//********************************************************************
double FitBase::RWSigmaToAbs(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.Eval(val);
return conv_val;
}
//********************************************************************
double FitBase::RWFracToSigma(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.GetX((val * f1.Eval(0.0)));
- if (fabs(conv_val) < 1E-10) conv_val = 0.0;
+ if (fabs(conv_val) < 1E-10)
+ conv_val = 0.0;
return conv_val;
}
//********************************************************************
double FitBase::RWSigmaToFrac(std::string const &type, std::string const &name,
double val) {
//********************************************************************
TF1 f1 = GetRWConvFunction(type, name);
double conv_val = f1.Eval(val) / f1.Eval(0.0);
return conv_val;
}
int FitBase::ConvDialType(std::string const &type) {
if (!type.compare("neut_parameter"))
return kNEUT;
else if (!type.compare("niwg_parameter"))
return kNIWG;
else if (!type.compare("nuwro_parameter"))
return kNUWRO;
else if (!type.compare("t2k_parameter"))
return kT2K;
else if (!type.compare("genie_parameter"))
return kGENIE;
else if (!type.compare("custom_parameter"))
return kCUSTOM;
else if (!type.compare("norm_parameter"))
return kNORM;
else if (!type.compare("likeweight_parameter"))
return kLIKEWEIGHT;
else if (!type.compare("spline_parameter"))
return kSPLINEPARAMETER;
else if (!type.compare("osc_parameter"))
return kOSCILLATION;
else if (!type.compare("modenorm_parameter"))
return kMODENORM;
+ else if (!type.compare("dunerwt_parameter"))
+ return kDUNERwt;
else
return kUNKNOWN;
}
std::string FitBase::ConvDialType(int type) {
switch (type) {
- case kNEUT: {
- return "neut_parameter";
- }
- case kNIWG: {
- return "niwg_parameter";
- }
- case kNUWRO: {
- return "nuwro_parameter";
- }
- case kT2K: {
- return "t2k_parameter";
- }
- case kGENIE: {
- return "genie_parameter";
- }
- case kNORM: {
- return "norm_parameter";
- }
- case kCUSTOM: {
- return "custom_parameter";
- }
- case kLIKEWEIGHT: {
- return "likeweight_parameter";
- }
- case kSPLINEPARAMETER: {
- return "spline_parameter";
- }
- case kOSCILLATION: {
- return "osc_parameter";
- }
- case kMODENORM: {
- return "modenorm_parameter";
- }
- default:
- return "unknown_parameter";
+ case kNEUT: {
+ return "neut_parameter";
+ }
+ case kNIWG: {
+ return "niwg_parameter";
+ }
+ case kNUWRO: {
+ return "nuwro_parameter";
+ }
+ case kT2K: {
+ return "t2k_parameter";
+ }
+ case kGENIE: {
+ return "genie_parameter";
+ }
+ case kNORM: {
+ return "norm_parameter";
+ }
+ case kCUSTOM: {
+ return "custom_parameter";
+ }
+ case kLIKEWEIGHT: {
+ return "likeweight_parameter";
+ }
+ case kSPLINEPARAMETER: {
+ return "spline_parameter";
+ }
+ case kOSCILLATION: {
+ return "osc_parameter";
+ }
+ case kMODENORM: {
+ return "modenorm_parameter";
+ }
+ case kDUNERwt: {
+ return "dunerwt_parameter";
+ }
+ default:
+ return "unknown_parameter";
}
}
int FitBase::GetDialEnum(std::string const &type, std::string const &name) {
return FitBase::GetDialEnum(FitBase::ConvDialType(type), name);
}
int FitBase::GetDialEnum(int type, std::string const &name) {
int offset = type * 1000;
- int this_enum = Reweight::kNoDialFound; // Not Found
+ int this_enum = Reweight::kNoDialFound; // Not Found
QLOG(DEB, "Getting dial enum " << type << " " << name);
// Select Types
switch (type) {
- // NEUT DIAL TYPE
- case kNEUT: {
+ // NEUT DIAL TYPE
+ case kNEUT: {
#ifdef __NEUT_ENABLED__
- int neut_enum = (int)neut::rew::NSyst::FromString(name);
- if (neut_enum != 0) {
- this_enum = neut_enum + offset;
- }
+ int neut_enum = (int)neut::rew::NSyst::FromString(name);
+ if (neut_enum != 0) {
+ this_enum = neut_enum + offset;
+ }
#else
- this_enum = Reweight::kNoTypeFound; // Not enabled
+ this_enum = Reweight::kNoTypeFound; // Not enabled
#endif
- break;
- }
+ break;
+ }
- // NIWG DIAL TYPE
- case kNIWG: {
+ // NIWG DIAL TYPE
+ case kNIWG: {
#ifdef __NIWG_ENABLED__
- int niwg_enum = (int)niwg::rew::NIWGSyst::FromString(name);
- if (niwg_enum != 0) {
- this_enum = niwg_enum + offset;
- }
+ int niwg_enum = (int)niwg::rew::NIWGSyst::FromString(name);
+ if (niwg_enum != 0) {
+ this_enum = niwg_enum + offset;
+ }
#else
- this_enum = Reweight::kNoTypeFound;
+ this_enum = Reweight::kNoTypeFound;
#endif
- break;
- }
+ break;
+ }
- // NUWRO DIAL TYPE
- case kNUWRO: {
+ // NUWRO DIAL TYPE
+ case kNUWRO: {
#ifdef __NUWRO_REWEIGHT_ENABLED__
- int nuwro_enum = (int)nuwro::rew::NuwroSyst::FromString(name);
- if (nuwro_enum > 0) {
- this_enum = nuwro_enum + offset;
- }
+ int nuwro_enum = (int)nuwro::rew::NuwroSyst::FromString(name);
+ if (nuwro_enum > 0) {
+ this_enum = nuwro_enum + offset;
+ }
#else
- this_enum = Reweight::kNoTypeFound;
+ this_enum = Reweight::kNoTypeFound;
#endif
- }
+ }
- // GENIE DIAL TYPE
- case kGENIE: {
+ // GENIE DIAL TYPE
+ case kGENIE: {
#ifdef __GENIE_ENABLED__
- int genie_enum = (int)genie::rew::GSyst::FromString(name);
- if (genie_enum > 0) {
- this_enum = genie_enum + offset;
- }
+ int genie_enum = (int)genie::rew::GSyst::FromString(name);
+ if (genie_enum > 0) {
+ this_enum = genie_enum + offset;
+ }
#else
- this_enum = Reweight::kNoTypeFound;
+ this_enum = Reweight::kNoTypeFound;
#endif
- break;
- }
+ break;
+ }
- case kCUSTOM: {
- int custom_enum = 0; // PLACEHOLDER
- this_enum = custom_enum + offset;
- break;
- }
+ case kCUSTOM: {
+ int custom_enum = 0; // PLACEHOLDER
+ this_enum = custom_enum + offset;
+ break;
+ }
- // T2K DIAL TYPE
- case kT2K: {
+ // T2K DIAL TYPE
+ case kT2K: {
#ifdef __T2KREW_ENABLED__
- int t2k_enum = (int)t2krew::T2KSyst::FromString(name);
- if (t2k_enum > 0) {
- this_enum = t2k_enum + offset;
- }
+ int t2k_enum = (int)t2krew::T2KSyst::FromString(name);
+ if (t2k_enum > 0) {
+ this_enum = t2k_enum + offset;
+ }
#else
- this_enum = Reweight::kNoTypeFound;
+ this_enum = Reweight::kNoTypeFound;
#endif
- break;
- }
+ break;
+ }
- case kNORM: {
- if (gNormEnums.find(name) == gNormEnums.end()) {
- gNormEnums[name] = gNormEnums.size() + 1 + offset;
- }
- this_enum = gNormEnums[name];
- break;
+ case kNORM: {
+ if (gNormEnums.find(name) == gNormEnums.end()) {
+ gNormEnums[name] = gNormEnums.size() + 1 + offset;
}
+ this_enum = gNormEnums[name];
+ break;
+ }
- case kLIKEWEIGHT: {
- if (gLikeWeightEnums.find(name) == gLikeWeightEnums.end()) {
- gLikeWeightEnums[name] = gLikeWeightEnums.size() + 1 + offset;
- }
- this_enum = gLikeWeightEnums[name];
- break;
+ case kLIKEWEIGHT: {
+ if (gLikeWeightEnums.find(name) == gLikeWeightEnums.end()) {
+ gLikeWeightEnums[name] = gLikeWeightEnums.size() + 1 + offset;
}
+ this_enum = gLikeWeightEnums[name];
+ break;
+ }
- case kSPLINEPARAMETER: {
- if (gSplineParameterEnums.find(name) == gSplineParameterEnums.end()) {
- gSplineParameterEnums[name] = gSplineParameterEnums.size() + 1 + offset;
- }
- this_enum = gSplineParameterEnums[name];
- break;
+ case kSPLINEPARAMETER: {
+ if (gSplineParameterEnums.find(name) == gSplineParameterEnums.end()) {
+ gSplineParameterEnums[name] = gSplineParameterEnums.size() + 1 + offset;
}
- case kOSCILLATION: {
+ this_enum = gSplineParameterEnums[name];
+ break;
+ }
+ case kOSCILLATION: {
#ifdef __PROB3PP_ENABLED__
- int oscEnum = OscWeightEngine::SystEnumFromString(name);
- if (oscEnum != 0) {
- this_enum = oscEnum + offset;
- }
+ int oscEnum = OscWeightEngine::SystEnumFromString(name);
+ if (oscEnum != 0) {
+ this_enum = oscEnum + offset;
+ }
#else
- this_enum = Reweight::kNoTypeFound; // Not enabled
+ this_enum = Reweight::kNoTypeFound; // Not enabled
#endif
+ }
+ case kMODENORM: {
+ size_t us_pos = name.find_first_of('_');
+ std::string numstr = name.substr(us_pos + 1);
+ int mode_num = std::atoi(numstr.c_str());
+ LOG(FTL) << "Getting mode num " << mode_num << std::endl;
+ if (!mode_num) {
+ ERR(FTL) << "Attempting to parse dial name: \"" << name
+ << "\" as a mode norm dial but failed." << std::endl;
+ throw;
}
- case kMODENORM: {
- size_t us_pos = name.find_first_of('_');
- std::string numstr = name.substr(us_pos + 1);
- int mode_num = std::atoi(numstr.c_str());
- LOG(FTL) << "Getting mode num " << mode_num << std::endl;
- if (!mode_num) {
- ERR(FTL) << "Attempting to parse dial name: \"" << name
- << "\" as a mode norm dial but failed." << std::endl;
- throw;
- }
- this_enum = 60 + mode_num + offset;
- break;
- }
+ this_enum = 60 + mode_num + offset;
+ break;
+ }
}
// If Not Enabled
if (this_enum == Reweight::kNoTypeFound) {
ERR(FTL) << "RW Engine not supported for " << FitBase::ConvDialType(type)
<< std::endl;
ERR(FTL) << "Check dial " << name << std::endl;
}
// If Not Found
if (this_enum == Reweight::kNoDialFound) {
ERR(FTL) << "Dial " << name << " not found." << std::endl;
}
return this_enum;
}
int Reweight::ConvDialType(std::string const &type) {
return FitBase::ConvDialType(type);
}
std::string Reweight::ConvDialType(int type) {
return FitBase::ConvDialType(type);
}
int Reweight::GetDialType(int type) {
int t = (type / 1000);
- return t > kMODENORM ? Reweight::kNoDialFound : t;
+ return t > kDUNERwt ? Reweight::kNoDialFound : t;
}
int Reweight::RemoveDialType(int type) { return (type % 1000); }
int Reweight::NEUTEnumFromName(std::string const &name) {
#ifdef __NEUT_ENABLED__
int neutenum = (int)neut::rew::NSyst::FromString(name);
return (neutenum > 0) ? neutenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::NIWGEnumFromName(std::string const &name) {
#ifdef __NIWG_ENABLED__
int niwgenum = (int)niwg::rew::NIWGSyst::FromString(name);
return (niwgenum != 0) ? niwgenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::NUWROEnumFromName(std::string const &name) {
#ifdef __NUWRO_REWEIGHT_ENABLED__
int nuwroenum = (int)nuwro::rew::NuwroSyst::FromString(name);
return (nuwroenum > 0) ? nuwroenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::GENIEEnumFromName(std::string const &name) {
#ifdef __GENIE_ENABLED__
int genieenum = (int)genie::rew::GSyst::FromString(name);
return (genieenum > 0) ? genieenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::T2KEnumFromName(std::string const &name) {
#ifdef __T2KREW_ENABLED__
int t2kenum = (int)t2krew::T2KSyst::FromString(name);
return (t2kenum > 0) ? t2kenum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::OscillationEnumFromName(std::string const &name) {
#ifdef __PROB3PP_ENABLED__
int oscEnum = OscWeightEngine::SystEnumFromString(name);
return (oscEnum > 0) ? oscEnum : Reweight::kNoDialFound;
#else
return Reweight::kGeneratorNotBuilt;
#endif
}
int Reweight::NUISANCEEnumFromName(std::string const &name, int type) {
int nuisenum = Reweight::DialList().EnumFromNameAndType(name, type);
return nuisenum;
}
int Reweight::CustomEnumFromName(std::string const &name) {
int custenum = Reweight::ConvertNUISANCEDial(name);
return custenum;
}
int Reweight::ConvDial(std::string const &name, std::string const &type,
bool exceptions) {
return Reweight::ConvDial(name, Reweight::ConvDialType(type), exceptions);
}
int Reweight::ConvDial(std::string const &fullname, int type, bool exceptions) {
std::string name =
- GeneralUtils::ParseToStr(fullname, ",")[0]; // Only use first dial given
+ GeneralUtils::ParseToStr(fullname, ",")[0]; // Only use first dial given
// Produce offset seperating each type.
int offset = type * 1000;
int genenum = Reweight::kNoDialFound;
switch (type) {
- case kNEUT:
- genenum = NEUTEnumFromName(name);
- break;
-
- case kNIWG:
- genenum = NIWGEnumFromName(name);
- break;
-
- case kNUWRO:
- genenum = NUWROEnumFromName(name);
- break;
-
- case kGENIE:
- genenum = GENIEEnumFromName(name);
- break;
-
- case kT2K:
- genenum = T2KEnumFromName(name);
- break;
-
- case kCUSTOM:
- genenum = CustomEnumFromName(name);
- break;
-
- case kNORM:
- case kLIKEWEIGHT:
- case kSPLINEPARAMETER:
- case kNEWSPLINE:
- genenum = NUISANCEEnumFromName(name, type);
- break;
-
- case kOSCILLATION:
- genenum = OscillationEnumFromName(name);
- break;
-
- case kMODENORM:
- genenum = ModeNormEngine::SystEnumFromString(name);
- break;
-
- default:
- genenum = Reweight::kNoTypeFound;
- break;
+ case kNEUT:
+ genenum = NEUTEnumFromName(name);
+ break;
+
+ case kNIWG:
+ genenum = NIWGEnumFromName(name);
+ break;
+
+ case kNUWRO:
+ genenum = NUWROEnumFromName(name);
+ break;
+
+ case kGENIE:
+ genenum = GENIEEnumFromName(name);
+ break;
+
+ case kT2K:
+ genenum = T2KEnumFromName(name);
+ break;
+
+ case kCUSTOM:
+ genenum = CustomEnumFromName(name);
+ break;
+
+ case kNORM:
+ case kLIKEWEIGHT:
+ case kSPLINEPARAMETER:
+ case kNEWSPLINE:
+ genenum = NUISANCEEnumFromName(name, type);
+ break;
+
+ case kOSCILLATION:
+ genenum = OscillationEnumFromName(name);
+ break;
+
+ case kMODENORM:
+ genenum = ModeNormEngine::SystEnumFromString(name);
+ break;
+
+ default:
+ genenum = Reweight::kNoTypeFound;
+ break;
}
// Throw if required.
if (exceptions) {
// If Not Enabled
if (genenum == Reweight::kGeneratorNotBuilt) {
ERR(FTL) << "RW Engine not supported for " << FitBase::ConvDialType(type)
<< std::endl;
ERR(FTL) << "Check dial " << name << std::endl;
throw;
}
// If no type enabled
if (genenum == Reweight::kNoTypeFound) {
ERR(FTL) << "Type mismatch inside ConvDialEnum" << std::endl;
throw;
}
// If Not Found
if (genenum == Reweight::kNoDialFound) {
ERR(FTL) << "Dial " << name << " not found." << std::endl;
throw;
}
}
// Add offset if no issue
int nuisenum = genenum;
if ((genenum != Reweight::kGeneratorNotBuilt) &&
(genenum != Reweight::kNoTypeFound) &&
(genenum != Reweight::kNoDialFound)) {
nuisenum += offset;
}
// Now register dial
Reweight::DialList().RegisterDialEnum(name, type, nuisenum);
return nuisenum;
}
std::string Reweight::ConvDial(int nuisenum) {
for (size_t i = 0; i < Reweight::DialList().fAllDialEnums.size(); i++) {
if (Reweight::DialList().fAllDialEnums[i] == nuisenum) {
return Reweight::DialList().fAllDialNames[i];
}
}
LOG(FIT) << "Cannot find dial with enum = " << nuisenum << std::endl;
return "";
}
diff --git a/src/Reweight/WeightUtils.h b/src/Reweight/WeightUtils.h
index c7b044f..30e22ea 100644
--- a/src/Reweight/WeightUtils.h
+++ b/src/Reweight/WeightUtils.h
@@ -1,63 +1,64 @@
#ifndef WEIGHTUTILS_H
#define WEIGHTUTILS_H
#include "FitEvent.h"
#include "TF1.h"
enum extra_reweight_types {
kOSCILLATION = kLast_generator_event_type,
- kMODENORM
+ kMODENORM,
+ kDUNERwt
};
namespace FitBase {
TF1 GetRWConvFunction(std::string const &type, std::string const &name);
std::string GetRWUnits(std::string const &type, std::string const &name);
double RWSigmaToFrac(std::string const &type, std::string const &name,
double val);
double RWSigmaToAbs(std::string const &type, std::string const &name,
double val);
double RWAbsToSigma(std::string const &type, std::string const &name,
double val);
double RWFracToSigma(std::string const &type, std::string const &name,
double val);
int ConvDialType(std::string const &type);
std::string ConvDialType(int type);
int GetDialEnum(std::string const &type, std::string const &name);
int GetDialEnum(int type, std::string const &name);
static std::map gNormEnums;
static std::map gLikeWeightEnums;
static std::map gSplineParameterEnums;
}
namespace Reweight {
int ConvDial(std::string const &name, std::string const &type,
bool exceptions = false);
int ConvDial(std::string const &name, int type, bool exceptions = false);
std::string ConvDial(int nuisenum);
int ConvDialType(std::string const &type);
std::string ConvDialType(int type);
int GetDialType(int type);
int RemoveDialType(int type);
int NEUTEnumFromName(std::string const &name);
int NIWGEnumFromName(std::string const &name);
int NUWROEnumFromName(std::string const &name);
int T2KEnumFromName(std::string const &name);
int GENIEEnumFromName(std::string const &name);
int CustomEnumFromName(std::string const &name);
int NUISANCEEnumFromName(std::string const &name, int type);
int OscillationEnumFromName(std::string const &name);
static const int kNoDialFound = -1;
static const int kNoTypeFound = -2;
static const int kGeneratorNotBuilt = -3;
}
#endif