diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt index 0349363..0c24400 100644 --- a/app/CMakeLists.txt +++ b/app/CMakeLists.txt @@ -1,103 +1,103 @@ # Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret ################################################################################ # This file is part of NUISANCE. # # NUISANCE is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # NUISANCE is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with NUISANCE. If not, see . ################################################################################ set(TARGETS_TO_BUILD) if(USE_MINIMIZER) add_executable(nuismin nuismin.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuismin) add_executable(nuissplines nuissplines.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissplines) endif() include_directories(${RWENGINE_INCLUDE_DIRECTORIES}) include_directories(${CMAKE_SOURCE_DIR}/src/Routines) include_directories(${CMAKE_SOURCE_DIR}/src/InputHandler) include_directories(${CMAKE_SOURCE_DIR}/src/Genie) include_directories(${CMAKE_SOURCE_DIR}/src/FitBase) include_directories(${CMAKE_SOURCE_DIR}/src/Statistical) include_directories(${CMAKE_SOURCE_DIR}/src/Utils) include_directories(${CMAKE_SOURCE_DIR}/src/Config) include_directories(${CMAKE_SOURCE_DIR}/src/Logger) include_directories(${CMAKE_SOURCE_DIR}/src/Splines) include_directories(${CMAKE_SOURCE_DIR}/src/Reweight) include_directories(${CMAKE_SOURCE_DIR}/src/FCN) include_directories(${CMAKE_SOURCE_DIR}/src/MCStudies) include_directories(${CMAKE_SOURCE_DIR}/src/Smearceptance) include_directories(${EXP_INCLUDE_DIRECTORIES}) if (USE_NuWro AND NOT NUWRO_BUILT_FROM_FILE) add_executable(nuwro_nuisance nuwro_NUISANCE.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuwro_nuisance) endif() if (USE_GiBUU) add_executable(DumpGiBUUEvents DumpGiBUUEvents.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};DumpGiBUUEvents) endif() add_executable(nuis_flat_tree_combiner nuis_flat_tree_combiner.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuis_flat_tree_combiner) add_executable(nuiscomp nuiscomp.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuiscomp) add_executable(nuisflat nuisflat.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuisflat) add_executable(nuissmear nuissmear.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissmear) add_executable(nuissyst nuissyst.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissyst) add_executable(nuisbayes nuisbayes.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuisbayes) if(USE_GENIE) add_executable(PrepareGENIE PrepareGENIE.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareGENIE) endif() -if(USE_NEUT) +if(USE_NEUT OR USE_NEUT_EVENT) add_executable(PrepareNEUT PrepareNEUT.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareNEUT) endif() # PREPARE NUWRO # Commented out for the time being until it is finished.. if(USE_NuWro) add_executable(PrepareNuwro PrepareNuwroEvents.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareNuwro) endif() add_executable(nuisbac nuisbac.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuisbac) foreach(targ ${TARGETS_TO_BUILD}) if(NOT "${CMAKE_LINK_FLAGS} ${NUIS_EXE_FLAGS}" STREQUAL " ") set_target_properties(${targ} PROPERTIES LINK_FLAGS "${CMAKE_LINK_FLAGS} ${NUIS_EXE_FLAGS}") endif() target_link_libraries(${targ} ${MODULETargets}) target_link_libraries(${targ} ${CMAKE_DEPENDLIB_FLAGS}) endforeach() install(TARGETS ${TARGETS_TO_BUILD} DESTINATION bin) \ No newline at end of file diff --git a/cmake/ROOTSetup.cmake b/cmake/ROOTSetup.cmake index 00520a2..910931a 100644 --- a/cmake/ROOTSetup.cmake +++ b/cmake/ROOTSetup.cmake @@ -1,172 +1,187 @@ # Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret ################################################################################ # This file is part of NUISANCE. # # NUISANCE is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # NUISANCE is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with NUISANCE. If not, see . ################################################################################ if ( NOT DEFINED ENV{ROOTSYS} ) cmessage (FATAL_ERROR "$ROOTSYS is not defined, please set up ROOT first.") else() cmessage(STATUS "Using ROOT installed at $ENV{ROOTSYS}") set(CMAKE_ROOTSYS $ENV{ROOTSYS}) endif() # Get cflags from ROOT execute_process (COMMAND root-config --cflags OUTPUT_VARIABLE ROOT_CXX_FLAGS_RAW OUTPUT_STRIP_TRAILING_WHITESPACE) string(REPLACE " " ";" ROOT_CXX_FLAGS "${ROOT_CXX_FLAGS_RAW}") # Get libdir from ROOT execute_process (COMMAND root-config --libdir OUTPUT_VARIABLE ROOT_LIBDIR OUTPUT_STRIP_TRAILING_WHITESPACE) # Get version from ROOT execute_process (COMMAND root-config --version OUTPUT_VARIABLE ROOT_VERSION OUTPUT_STRIP_TRAILING_WHITESPACE) # Get features from ROOT execute_process (COMMAND root-config --features OUTPUT_VARIABLE ROOT_FEATURES OUTPUT_STRIP_TRAILING_WHITESPACE) LIST(APPEND EXTRA_LINK_DIRS ${ROOT_LIBDIR}) LIST(APPEND ROOT_LIBS Core Cint RIO XMLIO Net Hist Graf Graf3d Gpad Tree Rint Postscript Matrix Physics MathMore MathCore Thread EG Geom GenVector) cmessage(STATUS "Checking ROOT version: ${ROOT_VERSION}") string(REGEX MATCH "^6.*" ROOTVERSIXMATCH ${ROOT_VERSION}) if(ROOTVERSIXMATCH) cmessage(STATUS "Using ROOT6, We are essentially flying blind here.") LIST(REMOVE_ITEM ROOT_LIBS Cint) LIST(APPEND EXTRA_CXX_FLAGS -DROOT6_USE_FIT_FITTER_INTERFACE -DROOT6) set(USE_ROOT6 True) else() string(REGEX MATCH "5.34/([0-9]+)" ROOTVERSMATCH ${ROOT_VERSION}) if(NOT ROOTVERSMATCH OR ${CMAKE_MATCH_1} LESS "19") cmessage(FATAL_ERROR "ROOT Version: ${ROOT_VERSION} has out of date minimizer interface, but minimizer functionality requested. Please configure with -DUSE_MINIMIZER=FALSE or update to 5.34/19 or greater to enable minimization features.") endif() endif() if(USE_MINIMIZER) if("${ROOT_FEATURES}" MATCHES "minuit2") cmessage(STATUS "ROOT built with MINUIT2 support") LIST(APPEND EXTRA_CXX_FLAGS -D__MINUIT2_ENABLED__) else() cmessage(FATAL_ERROR "ROOT built without MINUIT2 support but minimizer functionality requested. Either configure with -DUSE_MINIMIZER=FALSE or use a version of ROOT with MINUIT2 support.") endif() endif() if("${ROOT_FEATURES}" MATCHES "opengl") cmessage(STATUS "ROOT built with OpenGL support") LIST(APPEND ROOT_LIBS RGL) endif() if(DEFINED NEED_ROOTPYTHIA6 AND NEED_ROOTPYTHIA6) LIST(APPEND ROOT_LIBS EGPythia6 Pythia6) endif() +#Check what ROOT thinks the standard is, set that project-wide +# and then remove it from ROOT_CXX_FLAGS +list (FIND ROOT_CXX_FLAGS "-std=c++11" CPP11_INDEX) +list (FIND ROOT_CXX_FLAGS "-std=c++14" CPP14_INDEX) +list (FIND ROOT_CXX_FLAGS "-std=c++17" CPP17_INDEX) +if (${CPP11_INDEX} GREATER -1) + SET(CMAKE_CXX_STANDARD 11) +elseif (${CPP14_INDEX} GREATER -1) + SET(CMAKE_CXX_STANDARD 14) +elseif (${CPP17_INDEX} GREATER -1) + SET(CMAKE_CXX_STANDARD 17) +elseif (${CPP11_INDEX} GREATER -1) + SET(CMAKE_CXX_STANDARD 11) +endif() +list(REMOVE_ITEM ROOT_CXX_FLAGS "-std=c++11") +list(REMOVE_ITEM ROOT_CXX_FLAGS "-std=c++14") +list(REMOVE_ITEM ROOT_CXX_FLAGS "-std=c++17") + cmessage ( STATUS "[ROOT]: root-config --version: ${ROOT_VERSION} ") cmessage ( STATUS "[ROOT]: root-config --cflags : ${ROOT_CXX_FLAGS} ") cmessage ( STATUS "[ROOT]: root-config --libs : ${ROOT_LD_FLAGS} ") LIST(APPEND EXTRA_CXX_FLAGS ${ROOT_CXX_FLAGS}) #Helper functions for building dictionaries function(GenROOTDictionary OutputDictName Header LinkDef) get_directory_property(incdirs INCLUDE_DIRECTORIES) string(REPLACE ";" ";-I" LISTDIRINCLUDES "-I${incdirs}") string(REPLACE " " ";" LISTCPPFLAGS "${EXTRA_CXX_FLAGS}") - #ROOT5 CINT cannot handle it. - list(REMOVE_ITEM LISTCPPFLAGS "-std=c++11") - message(STATUS "LISTCPPFLAGS: ${LISTCPPFLAGS}") message(STATUS "LISTINCLUDES: ${LISTDIRINCLUDES}") #Learn how to generate the Dict.cxx and Dict.hxx add_custom_command( OUTPUT "${OutputDictName}.cxx" "${OutputDictName}.h" COMMAND rootcint ARGS -f ${OutputDictName}.cxx -c -p ${LISTDIRINCLUDES} ${LISTCPPFLAGS} ${Header} ${LinkDef} DEPENDS ${Header};${LinkDef}) endfunction() function(BuildROOTProject ProjectName InputFile CommaSeparatedClassesToDump LIBLINKMODE) string(REPLACE "," ";" HeadersToDump ${CommaSeparatedClassesToDump}) set(OUTPUTFILES ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectSource.cxx ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}LinkDef.h ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectHeaders.h ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectInstances.h) cmessage(STATUS "As part of ROOT project: ${ProjectName}") foreach (header ${HeadersToDump}) LIST(APPEND OUTPUTFILES "${CMAKE_BINARY_DIR}/${ProjectName}/${header}.h") cmessage(STATUS "Will generate: ${CMAKE_BINARY_DIR}/${ProjectName}/${header}.h") endforeach() add_custom_command( OUTPUT ${OUTPUTFILES} COMMAND ${CMAKE_BINARY_DIR}/src/Utils/DumpROOTClassesFromVector ARGS ${InputFile} ${CMAKE_BINARY_DIR}/${ProjectName} ${CommaSeparatedClassesToDump} VERBATIM DEPENDS DumpROOTClassesFromVector) add_custom_target(${ProjectName}_sources DEPENDS ${OUTPUTFILES}) GenROOTDictionary( ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectDict ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectHeaders.h ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}LinkDef.h ) add_custom_target(${ProjectName}ProjectDict DEPENDS ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectDict.cxx ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectDict.h ) # add_dependencies(${ProjectName}ProjectDict ${ProjectName}_sources) #ProjectSource.cxx includes ProjectDict.cxx, so no need to add to compilation. set(ROAA_SOURCEFILES ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectSource.cxx) add_library(${ProjectName} ${LIBLINKMODE} ${ROAA_SOURCEFILES}) add_dependencies(${ProjectName} ${ProjectName}ProjectDict) endfunction() diff --git a/cmake/T2KSetup.cmake b/cmake/T2KSetup.cmake index c67ef53..a4d3028 100644 --- a/cmake/T2KSetup.cmake +++ b/cmake/T2KSetup.cmake @@ -1,40 +1,160 @@ # Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret ################################################################################ # This file is part of NUISANCE. # # NUISANCE is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # NUISANCE is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with NUISANCE. If not, see . ################################################################################ -if(T2KREWEIGHT STREQUAL "") - cmessage(FATAL_ERROR "Variable T2KREWEIGHT is not defined. Either configure with -DT2KREWEIGHT or \"\$ export T2KREWEIGHT=/path/to/T2KReWeight\". This must be set to point to a prebuilt T2KReWeight instance.") +include(cmake/parseConfigApp.cmake) + +find_program(T2KRWCONFIG NAMES t2kreweight-config) +SET(HAVET2KRWCONFIG FALSE) +# We are dealing with shiny NEUT +if(NOT "${T2KRWCONFIG}" STREQUAL "T2KRWCONFIG-NOTFOUND") + SET(HAVET2KRWCONFIG TRUE) + cmessage(STATUS "Found neut-config, using it to determine configuration.") +else() + cmessage(STATUS "Failed to find neut-config, assuming older NEUT build.") endif() -LIST(APPEND EXTRA_CXX_FLAGS -D__T2KREW_ENABLED__ ) +if(HAVET2KRWCONFIG) + if(NOT DEFINED CMAKE_CXX_STANDARD OR "${CMAKE_CXX_STANDARD} " STREQUAL " ") + SET(CMAKE_CXX_STANDARD 11) + endif() -# First check the OAANALYSIS libs (need to grab some headers for T2KReWeight linking if compiled with oaAnalysisReader) -IF(NOT $ENV{OAANALYSISREADERROOT} STREQUAL "") - cmessage(STATUS "Found OAANALYSISREADERROOT $ENV{OAANALYSISREADERROOT}, appending...") - LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES $ENV{OAANALYSISREADERROOT}/$ENV{OAANALYSISREADERCONFIG}) - LIST(APPEND EXTRA_LINK_DIRS $ENV{OAANALYSISREADERROOT}/$ENV{OAANALYSISREADERCONFIG}) - LIST(APPEND EXTRA_LIBS oaAnalysisReader) - # Don't have to append this; should be appeneded in ${T2KReWeight/src/T2KBuild.h} - #LIST(APPEND EXTRA_CXX_FLAGS -D__T2KRW_OAANALYSIS_ENABLED__) -endif() + LIST(APPEND EXTRA_CXX_FLAGS -DT2KRW_OA2021_INTERFACE -D__T2KREW_ENABLED__) + + GETLIBDIRS(t2kreweight-config --linkflags T2KRW_LINK_DIRS) + GETINCDIRS(t2kreweight-config --cflags T2KRW_INC_DIRS) + GETLIBS(t2kreweight-config --linkflags T2KRW_LIBS) + + cmessage(STATUS "T2KReWeight:") + cmessage(STATUS " LINK DIRS: ${T2KRW_LINK_DIRS}") + cmessage(STATUS " INC DIRS: ${T2KRW_INC_DIRS}") + cmessage(STATUS " LIBS: ${T2KRW_LIBS}") + + LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${T2KRW_INC_DIRS}) + LIST(APPEND EXTRA_LINK_DIRS ${T2KRW_LINK_DIRS}) + LIST(APPEND EXTRA_LIBS ${T2KRW_LIBS}) + + + execute_process (COMMAND t2kreweight-config + --has-feature NIWG RESULT_VARIABLE T2KRW_HAS_NIWG) + + if("${T2KRW_HAS_NIWG} " STREQUAL "0 ") + + GETLIBDIRS(t2kreweight-config --niwgflags NIWG_LINK_DIRS) + GETINCDIRS(t2kreweight-config --niwgflags NIWG_INC_DIRS) + GETLIBS(t2kreweight-config --niwgflags NIWG_LIBS) + + cmessage(STATUS "NIWG:") + cmessage(STATUS " LINK DIRS: ${NIWG_LINK_DIRS}") + cmessage(STATUS " INC DIRS: ${NIWG_INC_DIRS}") + cmessage(STATUS " LIBS: ${NIWG_LIBS}") + + LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${NIWG_INC_DIRS}) + LIST(APPEND EXTRA_LINK_DIRS ${NIWG_LINK_DIRS}) + LIST(APPEND EXTRA_LIBS ${NIWG_LIBS}) + + endif() + + execute_process (COMMAND t2kreweight-config + --has-feature NEUT RESULT_VARIABLE T2KRW_HAS_NEUT) + + if("${T2KRW_HAS_NEUT} " STREQUAL "0 ") + + LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_EVENT_ENABLED) + set(USE_NEUT_EVENT TRUE) + set(USE_NEUT_EVENT TRUE CACHE BOOL "Whether to enable NEUT (event i/o) support. Requires external libraries. " FORCE) + + + GETLIBDIRS(t2kreweight-config --neutflags NEUT_LINK_DIRS) + GETINCDIRS(t2kreweight-config --neutflags NEUT_INC_DIRS) + GETLIBS(t2kreweight-config --neutflags NEUT_LIBS) + + cmessage(STATUS "NEUT:") + cmessage(STATUS " LINK DIRS: ${NEUT_LINK_DIRS}") + cmessage(STATUS " INC DIRS: ${NEUT_INC_DIRS}") + cmessage(STATUS " LIBS: ${NEUT_LIBS}") + + LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${NEUT_INC_DIRS}) + LIST(APPEND EXTRA_LINK_DIRS ${NEUT_LINK_DIRS}) + LIST(APPEND EXTRA_LIBS ${NEUT_LIBS}) + + endif() + + execute_process (COMMAND t2kreweight-config + --has-feature GEANT RESULT_VARIABLE T2KRW_HAS_GEANT) + + if("${T2KRW_HAS_GEANT} " STREQUAL "0 ") + + GETLIBDIRS(t2kreweight-config --geantflags GEANT_LINK_DIRS) + GETINCDIRS(t2kreweight-config --geantflags GEANT_INC_DIRS) + GETLIBS(t2kreweight-config --geantflags GEANT_LIBS) + + cmessage(STATUS "GEANT:") + cmessage(STATUS " LINK DIRS: ${GEANT_LINK_DIRS}") + cmessage(STATUS " INC DIRS: ${GEANT_INC_DIRS}") + cmessage(STATUS " LIBS: ${GEANT_LIBS}") + + LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${GEANT_INC_DIRS}) + LIST(APPEND EXTRA_LINK_DIRS ${GEANT_LINK_DIRS}) + LIST(APPEND EXTRA_LIBS ${GEANT_LIBS}) + + endif() + + execute_process (COMMAND t2kreweight-config + --has-feature ND280 RESULT_VARIABLE T2KRW_HAS_ND280) + + if("${T2KRW_HAS_ND280} " STREQUAL "0 ") + + GETLIBDIRS(t2kreweight-config --nd280flags ND280_LINK_DIRS) + GETINCDIRS(t2kreweight-config --nd280flags ND280_INC_DIRS) + GETLIBS(t2kreweight-config --nd280flags ND280_LIBS) + + cmessage(STATUS "ND280:") + cmessage(STATUS " LINK DIRS: ${ND280_LINK_DIRS}") + cmessage(STATUS " INC DIRS: ${ND280_INC_DIRS}") + cmessage(STATUS " LIBS: ${ND280_LIBS}") + + LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${ND280_INC_DIRS}) + LIST(APPEND EXTRA_LINK_DIRS ${ND280_LINK_DIRS}) + LIST(APPEND EXTRA_LIBS ${ND280_LIBS}) + + endif() + +else() + + if(T2KREWEIGHT STREQUAL "") + cmessage(FATAL_ERROR "Variable T2KREWEIGHT is not defined. Either configure with -DT2KREWEIGHT or \"\$ export T2KREWEIGHT=/path/to/T2KReWeight\". This must be set to point to a prebuilt T2KReWeight instance.") + endif() + + LIST(APPEND EXTRA_CXX_FLAGS -D__T2KREW_ENABLED__ ) -LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${T2KREWEIGHT}/src/) -LIST(APPEND EXTRA_LINK_DIRS ${T2KREWEIGHT}/lib) -LIST(APPEND EXTRA_LIBS T2KReWeight) + # First check the OAANALYSIS libs (need to grab some headers for T2KReWeight linking if compiled with oaAnalysisReader) + IF(NOT $ENV{OAANALYSISREADERROOT} STREQUAL "") + cmessage(STATUS "Found OAANALYSISREADERROOT $ENV{OAANALYSISREADERROOT}, appending...") + LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES $ENV{OAANALYSISREADERROOT}/$ENV{OAANALYSISREADERCONFIG}) + LIST(APPEND EXTRA_LINK_DIRS $ENV{OAANALYSISREADERROOT}/$ENV{OAANALYSISREADERCONFIG}) + LIST(APPEND EXTRA_LIBS oaAnalysisReader) + # Don't have to append this; should be appeneded in ${T2KReWeight/src/T2KBuild.h} + #LIST(APPEND EXTRA_CXX_FLAGS -D__T2KRW_OAANALYSIS_ENABLED__) + endif() + LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${T2KREWEIGHT}/src/) + LIST(APPEND EXTRA_LINK_DIRS ${T2KREWEIGHT}/lib) + LIST(APPEND EXTRA_LIBS T2KReWeight) +endif() \ No newline at end of file diff --git a/cmake/cacheVariables.cmake b/cmake/cacheVariables.cmake index c180ab6..5c81858 100644 --- a/cmake/cacheVariables.cmake +++ b/cmake/cacheVariables.cmake @@ -1,247 +1,248 @@ # Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret ################################################################################ # This file is part of NUISANCE. # # NUISANCE is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # NUISANCE is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with NUISANCE. If not, see . ################################################################################ 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 BOOL "Whether we are using the ROOT minimization libraries. ") CheckAndSetDefaultCache(USE_REWEIGHT TRUE BOOL "Whether we are expect to be able to build the reweighting libraries of enabled generators. ") CheckAndSetDefaultCache(USE_ROOT6 FALSE INTERNAL "Whether we are using the ROOT 6. ") CheckAndSetDefaultCache(USE_HEPMCNUEVT FALSE BOOL "Whether to enable HepMC3 input support. ") CheckAndSetDefaultCache(USE_NUSYST FALSE BOOL "Whether to enable DUNE Reweight support. ") CheckAndSetDefaultEnv(NUSYST_ROOT "" PATH "Path to nusystematics install directory <>" NUSYST_ROOT) CheckAndSetDefaultEnv(SYSTTOOLS_ROOT "" PATH "Path to systematicstools install directory <>" SYSTTOOLS_ROOT) 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(USE_NEUT_EVENT FALSE BOOL "Whether to enable NEUT (event i/o) support. Requires external libraries. ") CheckAndSetDefaultEnv(NEUT_VERSION FALSE STRING "NEUT version string, e.g. 5.4.0. <5.4.0>" NEUT_VERSION) 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 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(GENIE_REWEIGHT "" PATH "Path to GENIE ReWeight directory. Only relevant for GENIE v3+. Overrides environment variable \$GENIE_REWEIGHT <>" GENIE_REWEIGHT) CheckAndSetDefaultCache(GENIE_EMPMEC_REWEIGHT FALSE BOOL "Whether to use GENIE EMP MEC reweight (requires custom GENIE) ") CheckAndSetDefaultCache(USE_GENIE_XSECMEC FALSE BOOL "Whether to use GENIE MEC reweight (requires custom 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) CheckAndSetDefaultEnv(GSL_LIB "" PATH "Path to pre-built gsl libraries. Overrides environment variable \$GSL_LIB. <>" GSL_LIB) CheckAndSetDefaultEnv(GSL_INC "" PATH "Path to installed gsl headers. Overrides environment variable \$GSL_INC. <>" GSL_INC) CheckAndSetDefaultCache(USE_T2K FALSE BOOL "Whether to enable T2KReWeight support. Requires external libraries. ") CheckAndSetDefaultEnv(T2KREWEIGHT "" PATH "Path to installed T2KReWeight. Overrides environment variable \$T2KREWEIGHT. <>" T2KREWEIGHT) CheckAndSetDefaultCache(USE_NIWG FALSE BOOL "Whether to enable (T2K) NIWG ReWeight support. Requires external libraries. ") CheckAndSetDefaultEnv(NIWG_ROOT "" PATH "Path to installed NIWGReWeight. Overrides environment variable \$NIWG. <>" NIWG) CheckAndSetDefaultCache(USE_MINERvA_RW FALSE BOOL "Whether to enable MINERvA ReWeight support. ") CheckAndSetDefaultCache(USE_NOvARwgt FALSE BOOL "Whether to enable NOvA ReWeight support. ") CheckAndSetDefaultEnv(NOVARWGT "" PATH "Path to directory containing libPythia6.so. Overrides environment variable \$NOVARWGT <>" NOVARWGT) 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_VERSION NEUT_ROOT CERN CERN_LEVEL USE_NuWro NUWRO NUWRO_INC NUWRO_INPUT_FILE NUWRO_BUILT_FROM_FILE USE_GENIE GENIE_VERSION GENIE GENIE_REWEIGHT LHAPDF_LIB LHAPDF_INC LHAPATH LIBXML2_LIB LIBXML2_INC LOG4CPP_LIB LOG4CPP_INC GSL_LIB GSL_INC PYTHIA6 PYTHIA8 USE_PYTHIA8 USE_T2K T2KREWEIGHT USE_NIWG NIWG_ROOT USE_MINERvA_RW USE_NOvARwgt NOVARWGT USE_GiBUU BUILD_GiBUU USE_NUANCE USE_PROB3PP NO_EXTERNAL_UPDATE USE_GPERFTOOLS NO_EXPERIMENTS 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 9af163c..17d6155 100644 --- a/src/InputHandler/BaseFitEvt.cxx +++ b/src/InputHandler/BaseFitEvt.cxx @@ -1,263 +1,267 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* -* This file is part of NUISANCE. -* -* NUISANCE is free software: you can redistribute it and/or modify -* it under the terms of the GNU General Public License as published by -* the Free Software Foundation, either version 3 of the License, or -* (at your option) any later version. -* -* NUISANCE is distributed in the hope that it will be useful, -* but WITHOUT ANY WARRANTY; without even the implied warranty of -* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -* GNU General Public License for more details. -* -* You should have received a copy of the GNU General Public License -* along with NUISANCE. If not, see . -*******************************************************************************/ + * 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; for (int i = 0; i < 6; ++i) { CustomWeightArray[i] = 1.0; } fSplineCoeff = NULL; fSplineRead = NULL; fGenInfo = NULL; fType = 9999; -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) || defined(NEUT_EVENT_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; #endif #ifdef __NUANCE_ENABLED__ nuance_event = NULL; #endif #ifdef __GiBUU_ENABLED__ GiRead = NULL; #endif }; BaseFitEvt::~BaseFitEvt() { -#ifdef __NEUT_ENABLED__ - if (fNeutVect) delete fNeutVect; +#if defined(__NEUT_ENABLED__) || defined(NEUT_EVENT_ENABLED) + if (fNeutVect) + delete fNeutVect; #endif #ifdef __NUWRO_ENABLED__ #ifndef __USE_NUWRO_SRW_EVENTS__ - if (fNuwroEvent) delete fNuwroEvent; + if (fNuwroEvent) + delete fNuwroEvent; #endif #endif #ifdef __GENIE_ENABLED__ - if (genie_event) delete genie_event; + if (genie_event) + delete genie_event; #endif #ifdef __NUANCE_ENABLED__ - if (nuance_event) delete nuance_event; + if (nuance_event) + delete nuance_event; #endif #ifdef __GiBUU_ENABLED__ - if (GiRead) delete GiRead; + 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; SavedRWWeight = obj->SavedRWWeight; CustomWeight = obj->CustomWeight; for (int i = 0; i < 6; ++i) { CustomWeightArray[i] = obj->CustomWeightArray[i]; } fSplineCoeff = obj->fSplineCoeff; fSplineRead = obj->fSplineRead; fGenInfo = obj->fGenInfo; fType = obj->fType; -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) || defined(NEUT_EVENT_ENABLED) fNeutVect = obj->fNeutVect; #endif #ifdef __NUWRO_ENABLED__ fNuwroEvent = obj->fNuwroEvent; #endif #ifdef __GENIE_ENABLED__ genie_event = obj->genie_event; #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; for (int i = 0; i < 6; ++i) { CustomWeightArray[i] = other.CustomWeightArray[i]; } - fSplineCoeff = other.fSplineCoeff; fSplineRead = other.fSplineRead; fGenInfo = other.fGenInfo; fType = other.fType; -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) || defined(NEUT_EVENT_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; #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; for (int i = 0; i < 6; ++i) { CustomWeightArray[i] = other.CustomWeightArray[i]; } fSplineCoeff = other.fSplineCoeff; fSplineRead = other.fSplineRead; fGenInfo = other.fGenInfo; fType = other.fType; -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) || defined(NEUT_EVENT_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; #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; +void BaseFitEvt::ResetWeight() { + InputWeight = 1.0; #ifdef __GENIE_ENABLED__ for (int i = 0; i < 6; ++i) { CustomWeightArray[i] = 1.0; } #endif } double BaseFitEvt::GetWeight() { return InputWeight * RWWeight * CustomWeight; }; -#ifdef __NEUT_ENABLED__ -void BaseFitEvt::SetNeutVect(NeutVect* v) { +#if defined(__NEUT_ENABLED__) || defined(NEUT_EVENT_ENABLED) +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 8e655fb..6f2fd5d 100644 --- a/src/InputHandler/BaseFitEvt.h +++ b/src/InputHandler/BaseFitEvt.h @@ -1,147 +1,147 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef FITEVENTBASE_H_SEEN #define FITEVENTBASE_H_SEEN /*! * \addtogroup InputHandler * @{ */ -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) || defined(NEUT_EVENT_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__ #ifdef GENIE_PRE_R3 #include "EVGCore/EventRecord.h" #include "GHEP/GHepRecord.h" #include "Ntuple/NtpMCEventRecord.h" #else #include "Framework/EventGen/EventRecord.h" #include "Framework/GHEP/GHepRecord.h" #include "Framework/Ntuple/NtpMCEventRecord.h" #endif using namespace genie; #endif #ifdef __NUANCE_ENABLED__ #include "NuanceEvent.h" #endif #include "StdHepEvt.h" #include "SplineReader.h" #include "InputTypes.h" #include "GeneratorInfoBase.h" /// Base Event Class used to store just the generator event pointers class BaseFitEvt { public: /// Base Constructor BaseFitEvt(void); ~BaseFitEvt(); /// Copy base fit event pointers 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;}; // Global Event Variables/Weights 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 double CustomWeightArray[6]; ///< For custom tuning using arrays, e.g. NOvA MINERvA WS // Spline Info Coefficients and Readers float* fSplineCoeff; ///< ND Array of Spline Coefficients SplineReader* fSplineRead; ///< Spline Interpretter // Generator Info GeneratorInfoBase* fGenInfo; ///< Generator Variable Box UInt_t fType; ///< Generator Event Type -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) || defined(NEUT_EVENT_ENABLED) /// Setup Event Reading from NEUT Event 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; #endif 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 #endif #ifdef __NUANCE_ENABLED__ /// Setup Event Reading from NUANCE Event void SetNuanceEvent(NuanceEvent* e); NuanceEvent* nuance_event; ///< Pointer to Nuance reader #endif #ifdef __GiBUU_ENABLED__ 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/InputFactory.cxx b/src/InputHandler/InputFactory.cxx index f55676b..eedb2f0 100644 --- a/src/InputHandler/InputFactory.cxx +++ b/src/InputHandler/InputFactory.cxx @@ -1,122 +1,122 @@ // Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #include "FitEventInputHandler.h" #include "GENIEInputHandler.h" #include "GIBUUInputHandler.h" #include "HistogramInputHandler.h" #include "NEUTInputHandler.h" #include "NUANCEInputHandler.h" #include "NuWroInputHandler.h" #include "SigmaQ0HistogramInputHandler.h" #include "SplineInputHandler.h" #include "InputFactory.h" #include "TFile.h" namespace InputUtils { InputHandlerBase *CreateInputHandler(std::string const &handle, InputUtils::InputType inpType, std::string const &inputs) { InputHandlerBase *input = NULL; std::string newinputs = InputUtils::ExpandInputDirectories(inputs); switch (inpType) { case (kNEUT_Input): -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) || defined(NEUT_EVENT_ENABLED) input = new NEUTInputHandler(handle, newinputs); #else NUIS_ERR(FTL, "Tried to create NEUTInputHandler : " << handle << " " << inpType << " " << inputs); NUIS_ABORT("NEUT is not enabled!"); #endif break; case (kGENIE_Input): #ifdef __GENIE_ENABLED__ input = new GENIEInputHandler(handle, newinputs); #else NUIS_ERR(FTL, "Tried to create GENIEInputHandler : " << handle << " " << inpType << " " << inputs); NUIS_ABORT("GENIE is not enabled!"); #endif break; case (kNUWRO_Input): #ifdef __NUWRO_ENABLED__ input = new NuWroInputHandler(handle, newinputs); #else NUIS_ERR(FTL, "Tried to create NuWroInputHandler : " << handle << " " << inpType << " " << inputs); NUIS_ABORT("NuWro is not enabled!"); #endif break; case (kGiBUU_Input): #ifdef __GiBUU_ENABLED__ input = new GIBUUInputHandler(handle, newinputs); #else NUIS_ERR(FTL, "Tried to create GiBUUInputHandler : " << handle << " " << inpType << " " << inputs); NUIS_ABORT("GiBUU is not enabled!"); #endif break; case (kNUANCE_Input): #ifdef __NUANCE_ENABLED__ input = new NUANCEInputHandler(handle, newinputs); #else NUIS_ERR(FTL, "Tried to create NUANCEInputHandler : " << handle << " " << inpType << " " << inputs); NUIS_ABORT("NUANCE is not enabled!"); #endif break; case (kFEVENT_Input): input = new FitEventInputHandler(handle, newinputs); break; case (kEVSPLN_Input): input = new SplineInputHandler(handle, newinputs); break; case (kSIGMAQ0HIST_Input): input = new SigmaQ0HistogramInputHandler(handle, newinputs); break; case (kHISTO_Input): input = new HistoInputHandler(handle, newinputs); break; default: break; } /// Input failed if (!input) { NUIS_ABORT("Input handler creation failed!" << std::endl << "Generator Type " << inpType << " not enabled!"); } return input; }; } // namespace InputUtils diff --git a/src/InputHandler/NEUTInputHandler.cxx b/src/InputHandler/NEUTInputHandler.cxx index 0735945..58cb9b6 100644 --- a/src/InputHandler/NEUTInputHandler.cxx +++ b/src/InputHandler/NEUTInputHandler.cxx @@ -1,552 +1,556 @@ -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) || defined(NEUT_EVENT_ENABLED) #include "NEUTInputHandler.h" #include "InputUtils.h" #include "PlotUtils.h" #include "TTreePerfStats.h" #include "fsihistC.h" #include "necardC.h" #include "nefillverC.h" #include "neutcrsC.h" #include "neutfsipart.h" #include "neutfsivert.h" #include "neutmodelC.h" #include "neutparamsC.h" #include "neutpart.h" #include "neutrootTreeSingleton.h" #include "neutvect.h" #include "neworkC.h" #include "vcworkC.h" #include "posinnucC.h" #ifdef __NEUT_NUCFSI_ENABLED__ #include "neutnucfsistep.h" #include "neutnucfsivert.h" #include "nucleonfsihistC.h" #endif NEUTGeneratorInfo::~NEUTGeneratorInfo() { DeallocateParticleStack(); } void NEUTGeneratorInfo::AddBranchesToTree(TTree *tn) { tn->Branch("NEUTParticleN", fNEUTParticleN, "NEUTParticleN/I"); tn->Branch("NEUTParticleStatusCode", fNEUTParticleStatusCode, "NEUTParticleStatusCode[NEUTParticleN]/I"); tn->Branch("NEUTParticleAliveCode", fNEUTParticleAliveCode, "NEUTParticleAliveCode[NEUTParticleN]/I"); } void NEUTGeneratorInfo::SetBranchesFromTree(TTree *tn) { tn->SetBranchAddress("NEUTParticleN", &fNEUTParticleN); tn->SetBranchAddress("NEUTParticleStatusCode", &fNEUTParticleStatusCode); tn->SetBranchAddress("NEUTParticleAliveCode", &fNEUTParticleAliveCode); } void NEUTGeneratorInfo::AllocateParticleStack(int stacksize) { fNEUTParticleN = 0; fNEUTParticleStatusCode = new int[stacksize]; fNEUTParticleStatusCode = new int[stacksize]; } void NEUTGeneratorInfo::DeallocateParticleStack() { delete fNEUTParticleStatusCode; delete fNEUTParticleAliveCode; } void NEUTGeneratorInfo::FillGeneratorInfo(NeutVect *nevent) { Reset(); for (int i = 0; i < nevent->Npart(); i++) { fNEUTParticleStatusCode[i] = nevent->PartInfo(i)->fStatus; fNEUTParticleAliveCode[i] = nevent->PartInfo(i)->fIsAlive; fNEUTParticleN++; } } void NEUTGeneratorInfo::Reset() { for (int i = 0; i < fNEUTParticleN; i++) { fNEUTParticleStatusCode[i] = -1; fNEUTParticleAliveCode[i] = 9; } fNEUTParticleN = 0; } NEUTInputHandler::NEUTInputHandler(std::string const &handle, std::string const &rawinputs) { NUIS_LOG(SAM, "Creating NEUTInputHandler : " << handle); // Run a joint input handling fName = handle; // Setup the TChain fNEUTTree = new TChain("neuttree"); fSaveExtra = FitPar::Config().GetParB("SaveExtraNEUT"); fCacheSize = FitPar::Config().GetParI("CacheSize"); fMaxEvents = FitPar::Config().GetParI("MAXEVENTS"); // Loop over all inputs and grab flux, eventhist, and nevents std::vector inputs = InputUtils::ParseInputFileList(rawinputs); for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) { // Open File for histogram access TFile *inp_file = new TFile(inputs[inp_it].c_str(), "READ"); if (!inp_file or inp_file->IsZombie()) { NUIS_ABORT( "NEUT File IsZombie() at : '" << inputs[inp_it] << "'" << std::endl << "Check that your file paths are correct and the file exists!" << std::endl << "$ ls -lh " << inputs[inp_it]); } // Get Flux/Event hist TH1D *fluxhist = (TH1D *)inp_file->Get( (PlotUtils::GetObjectWithName(inp_file, "flux")).c_str()); TH1D *eventhist = (TH1D *)inp_file->Get( (PlotUtils::GetObjectWithName(inp_file, "evt")).c_str()); if (!fluxhist or !eventhist) { NUIS_ERR(FTL, "Input File Contents: " << inputs[inp_it]); inp_file->ls(); NUIS_ABORT("NEUT FILE doesn't contain flux/xsec info. You may have to " "regenerate your MC!"); } // Get N Events TTree *neuttree = (TTree *)inp_file->Get("neuttree"); if (!neuttree) { NUIS_ERR(FTL, "neuttree not located in NEUT file: " << inputs[inp_it]); NUIS_ABORT( "Check your inputs, they may need to be completely regenerated!"); throw; } int nevents = neuttree->GetEntries(); if (nevents <= 0) { NUIS_ABORT("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 fNEUTTree->AddFile(inputs[inp_it].c_str()); } // Registor all our file inputs SetupJointInputs(); // Assign to tree fEventType = kNEUT; fNeutVect = NULL; fNEUTTree->SetBranchAddress("vectorbranch", &fNeutVect); #if defined(ROOT6) && defined(__NEUT_VERSION__) && (__NEUT_VERSION__ >= 541) fNEUTTree->SetAutoDelete(true); #endif fNEUTTree->GetEntry(0); // Create Fit Event fNUISANCEEvent = new FitEvent(); fNUISANCEEvent->SetNeutVect(fNeutVect); if (fSaveExtra) { fNeutInfo = new NEUTGeneratorInfo(); fNUISANCEEvent->AddGeneratorInfo(fNeutInfo); } fNUISANCEEvent->HardReset(); }; NEUTInputHandler::~NEUTInputHandler(){ // if (fNEUTTree) delete fNEUTTree; // if (fNeutVect) delete fNeutVect; // if (fNeutInfo) delete fNeutInfo; }; void NEUTInputHandler::CreateCache() { if (fCacheSize > 0) { // fNEUTTree->SetCacheEntryRange(0, fNEvents); fNEUTTree->AddBranchToCache("vectorbranch", 1); fNEUTTree->SetCacheSize(fCacheSize); } } void NEUTInputHandler::RemoveCache() { // fNEUTTree->SetCacheEntryRange(0, fNEvents); fNEUTTree->AddBranchToCache("vectorbranch", 0); fNEUTTree->SetCacheSize(0); } FitEvent *NEUTInputHandler::GetNuisanceEvent(const UInt_t ent, const bool lightweight) { UInt_t entry = ent + fSkip; // Catch too large entries if (entry >= (UInt_t)fNEvents) return NULL; // Read Entry from TTree to fill NEUT Vect in BaseFitEvt; fNEUTTree->GetEntry(entry); // Run NUISANCE Vector Filler if (!lightweight) { CalcNUISANCEKinematics(); } #ifdef __PROB3PP_ENABLED__ else { UInt_t npart = fNeutVect->Npart(); for (size_t i = 0; i < npart; i++) { NeutPart *part = fNUISANCEEvent->fNeutVect->PartInfo(i); if ((part->fIsAlive == false) && (part->fStatus == -1) && std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4, part->fPID)) { fNUISANCEEvent->probe_E = part->fP.T(); fNUISANCEEvent->probe_pdg = part->fPID; break; } else { continue; } } } #endif // Setup Input scaling for joint inputs fNUISANCEEvent->InputWeight = GetInputWeight(entry); // Return event pointer return fNUISANCEEvent; } // From NEUT neutclass/neutpart.h // Bool_t fIsAlive; // Particle should be tracked or not // ( in the detector simulator ) // // Int_t fStatus; // Status flag of this particle // -2: Non existing particle // -1: Initial state particle // 0: Normal // 1: Decayed to the other particle // 2: Escaped from the detector // 3: Absorped // 4: Charge exchanged // 5: Pauli blocked // 6: N/A // 7: Produced child particles // 8: Inelastically scattered // int NEUTInputHandler::GetNeutParticleStatus(NeutPart *part) { // State int state = kUndefinedState; // Remove Pauli blocked events, probably just single pion events if (part->fStatus == 5) { state = kFSIState; // fStatus == -1 means initial state } else if (part->fIsAlive == false && part->fStatus == -1) { state = kInitialState; // NEUT has a bit of a strange convention for fIsAlive and fStatus // combinations // for NC and neutrino particle isAlive true/false and status 2 means // final state particle // for other particles in NC status 2 means it's an FSI particle // for CC it means it was an FSI particle } else if (part->fStatus == 2) { // NC case is a little strange... The outgoing neutrino might be alive or // not alive. Remaining particles with status 2 are FSI particles that // reinteracted if (abs(fNeutVect->Mode) > 30 && (abs(part->fPID) == 16 || abs(part->fPID) == 14 || abs(part->fPID) == 12)) { state = kFinalState; // The usual CC case } else if (part->fIsAlive == true) { state = kFSIState; } } else if (part->fIsAlive == true && part->fStatus == 2 && (abs(part->fPID) == 16 || abs(part->fPID) == 14 || abs(part->fPID) == 12)) { state = kFinalState; } else if (part->fIsAlive == true && part->fStatus == 0) { state = kFinalState; } else if (!part->fIsAlive && (part->fStatus == 1 || part->fStatus == 3 || part->fStatus == 4 || part->fStatus == 7 || part->fStatus == 8)) { state = kFSIState; // There's one hyper weird case where fStatus = -3. This apparently // corresponds to a nucleon being ejected via pion FSI when there is "data // available" } else if (!part->fIsAlive && (part->fStatus == -3)) { state = kUndefinedState; // NC neutrino outgoing } else if (!part->fIsAlive && part->fStatus == 0 && (abs(part->fPID) == 16 || abs(part->fPID) == 14 || abs(part->fPID) == 12)) { state = kFinalState; // Warn if we still find alive particles without classifying them } else if (part->fIsAlive == true) { NUIS_ABORT("Undefined NEUT state " << " Alive: " << part->fIsAlive << " Status: " << part->fStatus << " PDG: " << part->fPID << " Mode: " << fNeutVect->Mode); } else if (abs(fNeutVect->Mode) == 35){ NUIS_ERR(WRN, "Marking nonsensical CC difractive event as undefined " << " Alive: " << part->fIsAlive << " Status: " << part->fStatus << " PDG: " << part->fPID << " Mode: " << fNeutVect->Mode); state = kUndefinedState; // Warn if we find dead particles that we haven't classified } else { NUIS_ABORT("Undefined NEUT state " << " Alive: " << part->fIsAlive << " Status: " << part->fStatus << " PDG: " << part->fPID << " Mode: " << fNeutVect->Mode); } return state; } void NEUTInputHandler::CalcNUISANCEKinematics() { // Reset all variables fNUISANCEEvent->ResetEvent(); // Fill Globals fNUISANCEEvent->Mode = fNeutVect->Mode; fNUISANCEEvent->fEventNo = fNeutVect->EventNo; fNUISANCEEvent->fTargetA = fNeutVect->TargetA; fNUISANCEEvent->fTargetZ = fNeutVect->TargetZ; fNUISANCEEvent->fTargetH = fNeutVect->TargetH; fNUISANCEEvent->fBound = bool(fNeutVect->Ibound); if (fNUISANCEEvent->fBound) { fNUISANCEEvent->fTargetPDG = TargetUtils::GetTargetPDGFromZA( fNUISANCEEvent->fTargetZ, fNUISANCEEvent->fTargetA); } else { fNUISANCEEvent->fTargetPDG = 1000010010; } // Check Particle Stack UInt_t npart = fNeutVect->Npart(); UInt_t kmax = fNUISANCEEvent->kMaxParticles; if (npart > kmax) { NUIS_ERR(WRN, "NEUT has too many particles. Expanding stack."); fNUISANCEEvent->ExpandParticleStack(npart); } UInt_t nprimary = fNeutVect->Nprimary(); // Fill Particle Stack for (size_t i = 0; i < npart; i++) { // Get Current Count int curpart = fNUISANCEEvent->fNParticles; // Get NEUT Particle NeutPart *part = fNeutVect->PartInfo(i); // State int state = GetNeutParticleStatus(part); // Remove Undefined if (kRemoveUndefParticles && state == kUndefinedState) continue; // Remove FSI if (kRemoveFSIParticles && state == kFSIState) continue; // Remove Nuclear if (kRemoveNuclearParticles && (state == kNuclearInitial || state == kNuclearRemnant)) continue; // State fNUISANCEEvent->fParticleState[curpart] = state; // Is the paricle associated with the primary vertex? bool primary = false; // NEUT events are just popped onto the stack as primary, then continues to // be non-primary if (i < nprimary) primary = true; fNUISANCEEvent->fPrimaryVertex[curpart] = primary; // Mom fNUISANCEEvent->fParticleMom[curpart][0] = part->fP.X(); fNUISANCEEvent->fParticleMom[curpart][1] = part->fP.Y(); fNUISANCEEvent->fParticleMom[curpart][2] = part->fP.Z(); fNUISANCEEvent->fParticleMom[curpart][3] = part->fP.T(); // PDG fNUISANCEEvent->fParticlePDG[curpart] = part->fPID; // Add up particle count fNUISANCEEvent->fNParticles++; } // Save Extra Generator Info if (fSaveExtra) { fNeutInfo->FillGeneratorInfo(fNeutVect); } // Run Initial, FSI, Final, Other ordering. fNUISANCEEvent->OrderStack(); FitParticle *ISAnyLepton = fNUISANCEEvent->GetHMISAnyLeptons(); if (ISAnyLepton) { fNUISANCEEvent->probe_E = ISAnyLepton->E(); fNUISANCEEvent->probe_pdg = ISAnyLepton->PDG(); } return; } +#endif + +#ifdef NEED_FILL_NEUT_COMMONS + void NEUTUtils::FillNeutCommons(NeutVect *nvect) { // WARNING: This has only been implemented for a neuttree and not GENIE // This should be kept in sync with T2KNIWGUtils::GetNIWGEvent(TTree) // NEUT version info. Can't get it to compile properly with this yet // neutversion_.corev = nvect->COREVer; // neutversion_.nucev = nvect->NUCEVer; // neutversion_.nuccv = nvect->NUCCVer; // Documentation: See nework.h nework_.modene = nvect->Mode; nework_.numne = nvect->Npart(); #ifdef NEUT_COMMON_QEAV nemdls_.mdlqeaf = nvect->QEAVForm; #else nemdls_.mdlqeaf = nvect->QEVForm; #endif nemdls_.mdlqe = nvect->QEModel; nemdls_.mdlspi = nvect->SPIModel; nemdls_.mdldis = nvect->DISModel; nemdls_.mdlcoh = nvect->COHModel; neutcoh_.necohepi = nvect->COHModel; nemdls_.xmaqe = nvect->QEMA; nemdls_.xmvqe = nvect->QEMV; nemdls_.kapp = nvect->KAPPA; // nemdls_.sccfv = SCCFVdef; // nemdls_.sccfa = SCCFAdef; // nemdls_.fpqe = FPQEdef; nemdls_.xmaspi = nvect->SPIMA; nemdls_.xmvspi = nvect->SPIMV; nemdls_.xmares = nvect->RESMA; nemdls_.xmvres = nvect->RESMV; neut1pi_.xmanffres = nvect->SPIMA; neut1pi_.xmvnffres = nvect->SPIMV; neut1pi_.xmarsres = nvect->RESMA; neut1pi_.xmvrsres = nvect->RESMV; neut1pi_.neiff = nvect->SPIForm; neut1pi_.nenrtype = nvect->SPINRType; neut1pi_.rneca5i = nvect->SPICA5I; neut1pi_.rnebgscl = nvect->SPIBGScale; nemdls_.xmacoh = nvect->COHMA; nemdls_.rad0nu = nvect->COHR0; // nemdls_.fa1coh = nvect->COHA1err; // nemdls_.fb1coh = nvect->COHb1err; // neutdis_.nepdf = NEPDFdef; // neutdis_.nebodek = NEBODEKdef; neutcard_.nefrmflg = nvect->FrmFlg; neutcard_.nepauflg = nvect->PauFlg; neutcard_.nenefo16 = nvect->NefO16; neutcard_.nemodflg = nvect->ModFlg; // neutcard_.nenefmodl = 1; // neutcard_.nenefmodh = 1; // neutcard_.nenefkinh = 1; // neutpiabs_.neabspiemit = 1; nenupr_.iformlen = nvect->FormLen; neutpiless_.ipilessdcy = nvect->IPilessDcy; neutpiless_.rpilessdcy = nvect->RPilessDcy; neutpiless_.ipilessdcy = nvect->IPilessDcy; neutpiless_.rpilessdcy = nvect->RPilessDcy; neffpr_.fefqe = nvect->NuceffFactorPIQE; neffpr_.fefqeh = nvect->NuceffFactorPIQEH; neffpr_.fefinel = nvect->NuceffFactorPIInel; neffpr_.fefabs = nvect->NuceffFactorPIAbs; neffpr_.fefcx = nvect->NuceffFactorPICX; neffpr_.fefcxh = nvect->NuceffFactorPICXH; neffpr_.fefcoh = nvect->NuceffFactorPICoh; neffpr_.fefqehf = nvect->NuceffFactorPIQEHKin; neffpr_.fefcxhf = nvect->NuceffFactorPICXKin; neffpr_.fefcohf = nvect->NuceffFactorPIQELKin; for (int i = 0; i < nework_.numne; i++) { nework_.ipne[i] = nvect->PartInfo(i)->fPID; nework_.pne[i][0] = (float)nvect->PartInfo(i)->fP.X() / 1000; // VC(NE)WORK in M(G)eV nework_.pne[i][1] = (float)nvect->PartInfo(i)->fP.Y() / 1000; // VC(NE)WORK in M(G)eV nework_.pne[i][2] = (float)nvect->PartInfo(i)->fP.Z() / 1000; // VC(NE)WORK in M(G)eV } // fsihist.h // neutroot fills a dummy object for events with no FSI to prevent memory leak // when // reading the TTree, so check for it here if ((int)nvect->NfsiVert() == 1) { // An event with FSI must have at least two vertices // if (nvect->NfsiPart()!=1 || nvect->Fsiprob!=-1) // ERR(WRN) << "T2KNeutUtils::fill_neut_commons(TTree) NfsiPart!=1 or // Fsiprob!=-1 when NfsiVert==1" << std::endl; fsihist_.nvert = 0; fsihist_.nvcvert = 0; fsihist_.fsiprob = 1; } else { // Real FSI event fsihist_.nvert = (int)nvect->NfsiVert(); for (int ivert = 0; ivert < fsihist_.nvert; ivert++) { fsihist_.iflgvert[ivert] = nvect->FsiVertInfo(ivert)->fVertID; fsihist_.posvert[ivert][0] = (float)nvect->FsiVertInfo(ivert)->fPos.X(); fsihist_.posvert[ivert][1] = (float)nvect->FsiVertInfo(ivert)->fPos.Y(); fsihist_.posvert[ivert][2] = (float)nvect->FsiVertInfo(ivert)->fPos.Z(); } fsihist_.nvcvert = nvect->NfsiPart(); for (int ip = 0; ip < fsihist_.nvcvert; ip++) { fsihist_.abspvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomLab; fsihist_.abstpvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomNuc; fsihist_.ipvert[ip] = nvect->FsiPartInfo(ip)->fPID; fsihist_.iverti[ip] = nvect->FsiPartInfo(ip)->fVertStart; fsihist_.ivertf[ip] = nvect->FsiPartInfo(ip)->fVertEnd; fsihist_.dirvert[ip][0] = (float)nvect->FsiPartInfo(ip)->fDir.X(); fsihist_.dirvert[ip][1] = (float)nvect->FsiPartInfo(ip)->fDir.Y(); fsihist_.dirvert[ip][2] = (float)nvect->FsiPartInfo(ip)->fDir.Z(); } fsihist_.fsiprob = nvect->Fsiprob; } neutcrscom_.crsx = nvect->Crsx; neutcrscom_.crsy = nvect->Crsy; neutcrscom_.crsz = nvect->Crsz; neutcrscom_.crsphi = nvect->Crsphi; neutcrscom_.crsq2 = nvect->Crsq2; neuttarget_.numbndn = nvect->TargetA - nvect->TargetZ; neuttarget_.numbndp = nvect->TargetZ; neuttarget_.numfrep = nvect->TargetH; neuttarget_.numatom = nvect->TargetA; posinnuc_.ibound = nvect->Ibound; // put empty nucleon FSI history (since it is not saved in the NeutVect // format) // Comment out as NEUT does not have the necessary proton FSI information yet // nucleonfsihist_.nfnvert = 0; // nucleonfsihist_.nfnstep = 0; } #endif diff --git a/src/InputHandler/NEUTInputHandler.h b/src/InputHandler/NEUTInputHandler.h index 9fcfa6e..3842813 100644 --- a/src/InputHandler/NEUTInputHandler.h +++ b/src/InputHandler/NEUTInputHandler.h @@ -1,77 +1,77 @@ #ifndef NEUTINPUTHANDLER_H #define NEUTINPUTHANDLER_H -#ifdef __NEUT_ENABLED__ +#if defined(__NEUT_ENABLED__) || defined(NEUT_EVENT_ENABLED) #include "InputHandler.h" #include "TargetUtils.h" #include "neutvect.h" /// NEUT Generator Container to save extra particle status codes. class NEUTGeneratorInfo : public GeneratorInfoBase { public: NEUTGeneratorInfo() {}; virtual ~NEUTGeneratorInfo(); /// Assigns information to branches void AddBranchesToTree(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); /// Clear any dynamic arrays void DeallocateParticleStack(); /// Read extra NEUT information from the event void FillGeneratorInfo(NeutVect* nevent); /// Reset extra information to default/empty values void Reset(); int kMaxParticles; ///< Number of particles in stack int* fNEUTParticleStatusCode; ///= 541) + + // No need to vomit the contents of the card file all over my screen + StopTalking(); + +#ifdef T2KRW_OA2021_INTERFACE + if (!t2krew::T2KNEUTUtils::CardIsSet()) { + std::string neut_card = FitPar::Config().GetParS("NEUT_CARD"); + if (neut_card.size()) { + t2krew::T2KNEUTUtils::SetCardFile(neut_card); + } + } +#else std::string neut_card = FitPar::Config().GetParS("NEUT_CARD"); if (!neut_card.size()) { NUIS_ABORT( "[ERROR]: When using T2KReWeight must set NEUT_CARD config option."); } - // No need to vomit the contents of the card file all over my screen - StopTalking(); + NUIS_LOG(Fit, "Using NEUT card file: " << neut_card); t2krew::T2KNeutUtils::SetCardFile(neut_card); +#endif StartTalking(); #endif // Setup the NEUT Reweight engien fCalcName = name; NUIS_LOG(FIT, "Setting up T2K RW: " << fCalcName); // Create RW Engine suppressing cout StopTalking(); +#ifdef T2KRW_OA2021_INTERFACE + fT2KRW = t2krew::MakeT2KReWeightInstance(); +#else // Create Main RW Engine fT2KRW = new t2krew::T2KReWeight(); // Setup Sub RW Engines (Only activated for neut and niwg) fT2KNeutRW = new t2krew::T2KNeutReWeight(); fT2KNIWGRW = new t2krew::T2KNIWGReWeight(); fT2KRW->AdoptWghtEngine("fNeutRW", fT2KNeutRW); fT2KRW->AdoptWghtEngine("fNIWGRW", fT2KNIWGRW); fT2KRW->Reconfigure(); // allow cout again StartTalking(); // Set Abs Twk Config fIsAbsTwk = (FitPar::Config().GetParB("setabstwk")); - +#endif #else NUIS_ABORT("T2K RW NOT ENABLED"); #endif }; void T2KWeightEngine::IncludeDial(std::string name, double startval) { #ifdef __T2KREW_ENABLED__ - // Get First enum int nuisenum = Reweight::ConvDial(name, kT2K); // Setup Maps fEnumIndex[nuisenum]; // = std::vector(0); fNameIndex[name]; // = std::vector(0); // Split by commas std::vector allnames = GeneralUtils::ParseToStr(name, ","); for (uint i = 0; i < allnames.size(); i++) { std::string singlename = allnames[i]; - // Get RW +// Get RW +#ifdef T2KRW_OA2021_INTERFACE + int gensyst = t2krew::T2KSystToInt(fT2KRW->DialFromString(name)); +#else t2krew::T2KSyst_t gensyst = t2krew::T2KSyst::FromString(name); +#endif // Fill Maps int index = fValues.size(); fValues.push_back(0.0); fT2KSysts.push_back(gensyst); // Initialize dial NUIS_LOG(REC, "Registering " << singlename << " from " << name); +#ifndef T2KRW_OA2021_INTERFACE fT2KRW->Systematics().Include(gensyst); // If Absolute if (fIsAbsTwk) { fT2KRW->Systematics().SetAbsTwk(gensyst); } +#endif // Setup index fEnumIndex[nuisenum].push_back(index); fNameIndex[name].push_back(index); } // Set Value if given if (startval != _UNDEF_DIAL_VALUE_) { SetDialValue(nuisenum, startval); } #endif } void T2KWeightEngine::SetDialValue(int nuisenum, double val) { #ifdef __T2KREW_ENABLED__ std::vector indices = fEnumIndex[nuisenum]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; +#ifdef T2KRW_OA2021_INTERFACE + fT2KRW->SetDial_To_Value(t2krew::IntToT2KSyst(fT2KSysts[indices[i]]), val); +#else fT2KRW->Systematics().SetTwkDial(fT2KSysts[indices[i]], val); +#endif } #endif } void T2KWeightEngine::SetDialValue(std::string name, double val) { #ifdef __T2KREW_ENABLED__ std::vector indices = fNameIndex[name]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; +#ifdef T2KRW_OA2021_INTERFACE + fT2KRW->SetDial_To_Value(t2krew::IntToT2KSyst(fT2KSysts[indices[i]]), val); +#else fT2KRW->Systematics().SetTwkDial(fT2KSysts[indices[i]], val); +#endif } #endif } void T2KWeightEngine::Reconfigure(bool silent) { #ifdef __T2KREW_ENABLED__ // Hush now... if (silent) StopTalking(); // Reconf StopTalking(); fT2KRW->Reconfigure(); StartTalking(); // Shout again if (silent) StartTalking(); #endif } double T2KWeightEngine::CalcWeight(BaseFitEvt *evt) { double rw_weight = 1.0; #ifdef __T2KREW_ENABLED__ // Skip Non-NEUT - if (evt->fType != kNEUT){ - return 1.0; + if (evt->fType != kNEUT) { + return 1.0; } - + // Hush now StopTalking(); - // Get Weight For NEUT +// Get Weight For NEUT +#ifdef T2KRW_OA2021_INTERFACE + rw_weight = fT2KRW->CalcWeight(t2krew::Event::Make(evt->fNeutVect)); +#else rw_weight = fT2KRW->CalcWeight(evt->fNeutVect); - +#endif // Speak Now StartTalking(); #endif - if (!std::isnormal(rw_weight)){ + if (!std::isnormal(rw_weight)) { NUIS_ERR(WRN, "NEUT returned weight: " << rw_weight); rw_weight = 0; } // Return rw_weight return rw_weight; } diff --git a/src/Reweight/T2KWeightEngine.h b/src/Reweight/T2KWeightEngine.h index d7fe980..cb1b154 100644 --- a/src/Reweight/T2KWeightEngine.h +++ b/src/Reweight/T2KWeightEngine.h @@ -1,43 +1,52 @@ #ifndef WEIGHT_ENGINE_T2K_H #define WEIGHT_ENGINE_T2K_H #include "FitLogger.h" #ifdef __T2KREW_ENABLED__ +#ifdef T2KRW_OA2021_INTERFACE +#include "T2KReWeight/WeightEngines/T2KReWeightFactory.h" +#else #include "T2KNIWGReWeight.h" #include "T2KNIWGUtils.h" #include "T2KNeutReWeight.h" #include "T2KNeutUtils.h" #include "T2KReWeight.h" using namespace t2krew; #endif +#endif +#include "FitWeight.h" #include "GeneratorUtils.h" #include "WeightEngineBase.h" -#include "FitWeight.h" class T2KWeightEngine : public WeightEngineBase { public: - T2KWeightEngine(std::string name); - ~T2KWeightEngine() {}; + T2KWeightEngine(std::string name); + ~T2KWeightEngine(){}; - void IncludeDial(std::string name, double startval); + void IncludeDial(std::string name, double startval); - void SetDialValue(std::string name, double val); - void SetDialValue(int nuisenum, double val); + void SetDialValue(std::string name, double val); + void SetDialValue(int nuisenum, double val); - void Reconfigure(bool silent = false); + void Reconfigure(bool silent = false); - double CalcWeight(BaseFitEvt* evt); + double CalcWeight(BaseFitEvt *evt); + + inline bool NeedsEventReWeight() { return true; }; - inline bool NeedsEventReWeight() { return true; }; - #ifdef __T2KREW_ENABLED__ - std::vector fT2KSysts; - t2krew::T2KReWeight* fT2KRW; //!< T2K RW Object - t2krew::T2KNeutReWeight* fT2KNeutRW; - t2krew::T2KNIWGReWeight* fT2KNIWGRW; +#ifdef T2KRW_OA2021_INTERFACE + std::vector fT2KSysts; + std::unique_ptr fT2KRW; +#else + std::vector fT2KSysts; + t2krew::T2KReWeight *fT2KRW; //!< T2K RW Object + t2krew::T2KNeutReWeight *fT2KNeutRW; + t2krew::T2KNIWGReWeight *fT2KNIWGRW; +#endif #endif }; #endif diff --git a/src/Reweight/WeightUtils.cxx b/src/Reweight/WeightUtils.cxx index 59cfa90..b8e0335 100644 --- a/src/Reweight/WeightUtils.cxx +++ b/src/Reweight/WeightUtils.cxx @@ -1,613 +1,647 @@ #include "WeightUtils.h" #include "FitLogger.h" #ifndef __NO_REWEIGHT__ #ifdef __T2KREW_ENABLED__ +#ifdef T2KRW_OA2021_INTERFACE +#include "T2KReWeight/WeightEngines/T2KReWeightFactory.h" +#include "T2KReWeight/WeightEngines/NEUT/T2KNEUTUtils.h" +#else #include "T2KGenieReWeight.h" #include "T2KNIWGReWeight.h" #include "T2KNIWGUtils.h" #include "T2KNeutReWeight.h" #include "T2KNeutUtils.h" #include "T2KReWeight.h" using namespace t2krew; #endif +#endif #ifdef __NIWG_ENABLED__ #include "NIWGReWeight.h" #include "NIWGSyst.h" #endif #ifdef __NEUT_ENABLED__ #include "NReWeight.h" #include "NSyst.h" #endif #ifdef __NUWRO_REWEIGHT_ENABLED__ #include "NuwroReWeight.h" #include "NuwroSyst.h" #endif #ifdef __GENIE_ENABLED__ #ifdef GENIE_PRE_R3 #ifndef __NO_GENIE_REWEIGHT__ #include "ReWeight/GReWeight.h" #include "ReWeight/GSyst.h" #endif #else using namespace genie; #ifndef __NO_GENIE_REWEIGHT__ #include "RwFramework/GReWeight.h" #include "RwFramework/GSyst.h" using namespace genie::rew; #endif #endif #endif #ifdef __NOVA_ENABLED__ #include "NOvARwgtEngine.h" #endif #ifdef __NUSYST_ENABLED__ #include "nusystematicsWeightEngine.h" #endif #endif // end of no reweight #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"; 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; // Check whether this is a comment if (inputlist[0].c_str()[0] == '#') continue; // Check whether this is the correct parameter type if (inputlist[0].compare(parType) != 0) continue; // Check the parameter name 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]); } 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; // Check whether this is a comment if (inputlist[0].c_str()[0] == '#') continue; // Check whether this is the correct parameter type if (inputlist[0].compare(parType) != 0) continue; // Check the parameter name 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; NUIS_LOG(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; 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("nova_parameter")) return kNOvARWGT; else if (!type.compare("nusyst_parameter")) return kNuSystematics; 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"; } case kNOvARWGT: { return "nova_parameter"; } case kNuSystematics: { return "nusyst_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 NUIS_LOG(DEB, "Getting dial enum " << type << " " << name); // Select Types switch (type) { // NEUT DIAL TYPE case kNEUT: { #if defined(__NEUT_ENABLED__) && !defined(__NO_REWEIGHT__) 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 #endif break; } // NIWG DIAL TYPE case kNIWG: { #if defined(__NIWG_ENABLED__) && !defined(__NO_REWEIGHT__) int niwg_enum = (int)niwg::rew::NIWGSyst::FromString(name); if (niwg_enum != 0) { this_enum = niwg_enum + offset; } #else this_enum = Reweight::kNoTypeFound; #endif break; } // NUWRO DIAL TYPE case kNUWRO: { #if defined(__NUWRO_REWEIGHT_ENABLED__) && !defined(__NO_REWEIGHT__) int nuwro_enum = (int)nuwro::rew::NuwroSyst::FromString(name); if (nuwro_enum > 0) { this_enum = nuwro_enum + offset; } #else this_enum = Reweight::kNoTypeFound; #endif } // GENIE DIAL TYPE case kGENIE: { #if defined(__GENIE_ENABLED__) && !defined(__NO_REWEIGHT__) int genie_enum = (int)genie::rew::GSyst::FromString(name); if (genie_enum > 0) { this_enum = genie_enum + offset; } #else this_enum = Reweight::kNoTypeFound; #endif break; } case kCUSTOM: { int custom_enum = 0; // PLACEHOLDER this_enum = custom_enum + offset; break; } // T2K DIAL TYPE case kT2K: { #if defined(__T2KREW_ENABLED__) && !defined(__NO_REWEIGHT__) +#ifdef T2KRW_OA2021_INTERFACE + // This is possibly inefficient, this should probably not be called per fit + // step. + if (!t2krew::T2KNEUTUtils::CardIsSet()) { + std::string neut_card = FitPar::Config().GetParS("NEUT_CARD"); + if (neut_card.size()) { + t2krew::T2KNEUTUtils::SetCardFile(neut_card); + } + } + + int t2k_enum = t2krew::T2KSystToInt( + t2krew::MakeT2KReWeightInstance()->DialFromString(name)); +#else int t2k_enum = (int)t2krew::T2KSyst::FromString(name); +#endif if (t2k_enum > 0) { this_enum = t2k_enum + offset; } #else this_enum = Reweight::kNoTypeFound; #endif 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 kSPLINEPARAMETER: { if (gSplineParameterEnums.find(name) == gSplineParameterEnums.end()) { gSplineParameterEnums[name] = gSplineParameterEnums.size() + 1 + offset; } this_enum = gSplineParameterEnums[name]; break; } case kOSCILLATION: { #ifdef __PROB3PP_ENABLED__ int oscEnum = OscWeightEngine::SystEnumFromString(name); if (oscEnum != 0) { this_enum = oscEnum + offset; } #else 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()); NUIS_LOG(FTL, "Getting mode num " << mode_num); if (!mode_num) { NUIS_ABORT("Attempting to parse dial name: \"" << name << "\" as a mode norm dial but failed."); } this_enum = 60 + mode_num + offset; break; } } // If Not Enabled if (this_enum == Reweight::kNoTypeFound) { NUIS_ERR(FTL, "RW Engine not supported for " << FitBase::ConvDialType(type)); NUIS_ABORT("Check dial " << name); } // If Not Found if (this_enum == Reweight::kNoDialFound) { NUIS_ABORT("Dial " << name << " not found."); } 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 > kNuSystematics ? Reweight::kNoDialFound : t; } int Reweight::RemoveDialType(int type) { return (type % 1000); } int Reweight::NEUTEnumFromName(std::string const &name) { #if defined(__NEUT_ENABLED__) && !defined(__NO_REWEIGHT__) 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) { #if defined(__NIWG_ENABLED__) && !defined(__NO_REWEIGHT__) 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) { #if defined(__NUWRO_REWEIGHT_ENABLED__) && !defined(__NO_REWEIGHT__) 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) { #if defined(__GENIE_ENABLED__) && !defined(__NO_REWEIGHT__) 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) { #if defined(__T2KREW_ENABLED__) && !defined(__NO_REWEIGHT__) +#ifdef T2KRW_OA2021_INTERFACE + // This is possibly inefficient, this should probably not be called per fit + // step. + if (!t2krew::T2KNEUTUtils::CardIsSet()) { + std::string neut_card = FitPar::Config().GetParS("NEUT_CARD"); + if (neut_card.size()) { + t2krew::T2KNEUTUtils::SetCardFile(neut_card); + } + } + + int t2kenum = t2krew::T2KSystToInt( + t2krew::MakeT2KReWeightInstance()->DialFromString(name)); +#else int t2kenum = (int)t2krew::T2KSyst::FromString(name); +#endif 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 != kUnknownNUISANCEDial ? custenum : kNoDialFound); } -int Reweight::ConvDial(std::string const &name, std::string const &type, bool exceptions) { +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 // 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; #ifdef __NOVA_ENABLED__ case kNOvARWGT: genenum = NOvARwgtEngine::GetWeightGeneratorIndex(name); break; #endif #ifdef __NUSYST_ENABLED__ case kNuSystematics: { // Super inefficient... nusystematicsWeightEngine we; genenum = we.ConvDial(name); break; } #endif default: genenum = Reweight::kNoTypeFound; break; } // Throw if required. if (exceptions) { // If Not Enabled if (genenum == Reweight::kGeneratorNotBuilt) { NUIS_ERR(FTL, "RW Engine not supported for " << FitBase::ConvDialType(type)); NUIS_ABORT("Check dial " << name); } // If no type enabled if (genenum == Reweight::kNoTypeFound) { NUIS_ABORT("Type mismatch inside ConvDialEnum"); } // If Not Found if (genenum == Reweight::kNoDialFound) { NUIS_ABORT("Dial " << name << " not found."); } } // 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]; } } NUIS_LOG(FIT, "Cannot find dial with enum = " << nuisenum); return ""; }