diff --git a/CMakeLists.txt b/CMakeLists.txt index bea638b..78fe008 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,83 +1,80 @@ # Copyright 2018 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) set (NUISANCE_VERSION_MAJOR 3) set (NUISANCE_VERSION_MINOR 0) 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) ################################# InputHandler ################################# include(${CMAKE_SOURCE_DIR}/cmake/InputHandlerSetup.cmake) ################################# Pythia6/8 ################################### include(${CMAKE_SOURCE_DIR}/cmake/pythia6Setup.cmake) -################################## FHICLCPP #################################### -include(${CMAKE_SOURCE_DIR}/cmake/fhiclcppSetup.cmake) - -#Need this to be at the front -LIST(REVERSE EXTRA_CXX_FLAGS) +#Want this before fhiclcpp which will add the install directory LIST(APPEND EXTRA_CXX_FLAGS -I${CMAKE_SOURCE_DIR}/src) -LIST(REVERSE EXTRA_CXX_FLAGS) +################################## FHICLCPP #################################### +include(${CMAKE_SOURCE_DIR}/cmake/fhiclcppSetup.cmake) ################################## COMPILER #################################### include(${CMAKE_SOURCE_DIR}/cmake/c++CompilerSetup.cmake) ################################################################################ add_subdirectory(src) 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}) add_subdirectory(config) add_subdirectory(data) diff --git a/checksamples.sh b/checksamples.sh index a0356af..edb23e5 100644 --- a/checksamples.sh +++ b/checksamples.sh @@ -1,29 +1,29 @@ for line in $(cat ./src/FCN/SampleList.cxx); #$(grep compare ../src/FCN/SampleList.cxx); do if [[ "$line" == *"compare"* ]]; then parsed=${line/'!name.compare'/} parsed=${parsed/'('/} parsed=${parsed/')'/} parsed=${parsed/'"'/} parsed=${parsed/'("'/} parsed=${parsed/'")'/} echo $parsed fi if [[ "$line" == *"find("* ]]; then parsed=${line/'!name.find'/} parsed=${parsed/'('/} parsed=${parsed/')'/} parsed=${parsed/'"'/} parsed=${parsed/'("'/} parsed=${parsed/'")'/} - echo "MultiSample: $line" + echo "MultIEventProcessor: $line" fi done; diff --git a/cmake/NEUTSetup.cmake b/cmake/NEUTSetup.cmake index 619d09f..2f7c60b 100644 --- a/cmake/NEUTSetup.cmake +++ b/cmake/NEUTSetup.cmake @@ -1,114 +1,130 @@ # 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(NEUT_ROOT STREQUAL "") cmessage(FATAL_ERROR "Variable NEUT_ROOT is not defined. Please export environment variable NEUT_ROOT or configure with -DNEUT_ROOT=/path/to/NEUT. This must be set to point to a prebuilt NEUT instance.") endif() if(CERN STREQUAL "") cmessage(FATAL_ERROR "Variable CERN is not defined. Please export environment variable CERN or configure with -DCERN=/path/to/CERNLIB. This must be set to point to a prebuilt CERNLIB instance.") endif() if(CERN_LEVEL STREQUAL "") cmessage(FATAL_ERROR "Variable CERN_LEVEL is not defined. Please export environment variable CERN_LEVEL or configure with -DCERN_LEVEL=XXXX (likely to be 2005).") endif() -if(NOT IS_NEUT_54) +if(${NEUT_VERSION} VERSION_LESS 5.4.0) set(NEUT_LIB_DIR ${NEUT_ROOT}/lib/Linux_pc) else() set(NEUT_LIB_DIR ${NEUT_ROOT}/lib) endif() + set(NEUT_CLASS ${NEUT_ROOT}/src/neutclass) -LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_ENABLED ) +LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_ENABLED -DNEUT_VERSION=${NEUT_VERSION}) + +if(${NEUT_VERSION} VERSION_GREATER 5.4.0) + LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_COMMON_QEAV) +endif() LIST(APPEND EXTRA_CXX_FLAGS -I${NEUT_ROOT}/include -I${NEUT_ROOT}/src/neutclass -I${NEUT_ROOT}/src/reweight) LIST(APPEND EXTRA_LINK_DIRS ${NEUT_LIB_DIR} ${CERN}/${CERN_LEVEL}/lib ${NEUT_ROOT}/src/reweight) -if(IS_NEUT_54) +LIST(APPEND EXTRA_LIBS NReWeight) + +if(${NEUT_VERSION} VERSION_EQUAL 5.4.2) + LIST(APPEND EXTRA_LIBS + neutcore_5.4.2 + nuccorspl_5.4.2 #typo in NEUT, may hopefully disappear + nuceff_5.4.2 + partnuck_5.4.2 + skmcsvc_5.4.2 + tauola_5.4.2 + HT2p2h_5.4.0 + N1p1h_5.4.0) + LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_COMMON_QEAV) +elseif(${NEUT_VERSION} VERSION_EQUAL 5.4.0) LIST(APPEND EXTRA_LIBS NReWeight neutcore_5.4.0 nuccorspl_5.4.0 #typo in NEUT, may hopefully disappear nuceff_5.4.0 partnuck_5.4.0 skmcsvc_5.4.0 tauola_5.4.0 HT2p2h_5.4.0 - N1p1h_5.4.0 - jetset74 - pdflib804 - mathlib - packlib - pawlib) + N1p1h_5.4.0) else() LIST(APPEND EXTRA_LIBS NReWeight neutcore nuccorrspl nuceff partnuck skmcsvc - tauola - jetset74 - pdflib804 - mathlib - packlib - pawlib) + tauola) endif() +LIST(APPEND EXTRA_LIBS + jetset74 + pdflib804 + mathlib + packlib + pawlib + gfortran) + set(NEUT_ROOT_LIBS) LIST(APPEND NEUT_ROOT_LIBS ${NEUT_CLASS}/neutctrl.so ${NEUT_CLASS}/neutfsivert.so) # Check for new versions of NEUT with NUCLEON FSI if(EXISTS "${NEUT_CLASS}/neutnucfsistep.so") set(NEUT_NUCFSI 1) LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_NUCFSI_ENABLED) LIST(APPEND NEUT_ROOT_LIBS ${NEUT_CLASS}/neutnucfsistep.so ${NEUT_CLASS}/neutnucfsivert.so ) endif() -if(NOT IS_NEUT_54) +if(${NEUT_VERSION} VERSION_LESS 5.4.0) LIST(APPEND NEUT_ROOT_LIBS ${NEUT_CLASS}/neutrootTreeSingleton.so) endif() LIST(APPEND NEUT_ROOT_LIBS ${NEUT_CLASS}/neutvtx.so ${NEUT_CLASS}/neutfsipart.so ${NEUT_CLASS}/neutpart.so ${NEUT_CLASS}/neutvect.so ) foreach(OBJ ${NEUT_ROOT_LIBS}) LIST(APPEND EXTRA_SHAREDOBJS ${OBJ}) endforeach() diff --git a/cmake/c++CompilerSetup.cmake b/cmake/c++CompilerSetup.cmake index 0443100..a6dd6cc 100644 --- a/cmake/c++CompilerSetup.cmake +++ b/cmake/c++CompilerSetup.cmake @@ -1,83 +1,83 @@ # Copyright 2018 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(CXX_WARNINGS -Wall ) -LIST(APPEND EXTRA_CXX_FLAGS -std=c++1y -fPIC "-D__FILENAME__='\"$(subst ${CMAKE_SOURCE_DIR}/,,$(abspath $<))\"'") +LIST(APPEND EXTRA_CXX_FLAGS -Wall -Wextra -Werror -Wno-delete-non-virtual-dtor -Wno-unused "-D__FILENAME__='\"$(subst ${CMAKE_SOURCE_DIR}/,,$(abspath $<))\"'") 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_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0") set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O3") SET(STR_EXTRA_LINK_DIRS) if(NOT EXTRA_LINK_DIRS STREQUAL "") string(REPLACE ";" " -L" STR_EXTRA_LINK_DIRS "-L${EXTRA_LINK_DIRS}") endif() LIST(APPEND EXTRA_LIBS dl) 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() LIST(APPEND EXTRA_LINK_FLAGS -Wl,--no-as-needed) 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(NOT STR_EXTRA_SHAREDOBJS STREQUAL "") SET(CMAKE_DEPENDLIB_FLAGS "${CMAKE_DEPENDLIB_FLAGS} ${STR_EXTRA_SHAREDOBJS}") 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 (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 11f3ae4..0c45699 100644 --- a/cmake/cacheVariables.cmake +++ b/cmake/cacheVariables.cmake @@ -1,99 +1,99 @@ # Copyright 2018 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() # Misc 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. ") # NuWro 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) # NEUT 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_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) # Pythia CheckAndSetDefaultEnv(PYTHIA6 "" PATH "Path to directory containing libPythia6.so. Overrides environment variable \$PYTHIA6 <>" PYTHIA6) CheckAndSetDefault(NEED_PYTHIA6 FALSE) CheckAndSetDefault(NEED_ROOTEVEGEN FALSE) CheckAndSetDefault(NEED_ROOTPYTHIA6 FALSE) diff --git a/config/global/nuis.global.config.fcl.in b/config/global/nuis.global.config.fcl.in index 5b964a4..2682d54 100644 --- a/config/global/nuis.global.config.fcl.in +++ b/config/global/nuis.global.config.fcl.in @@ -1,18 +1,20 @@ plugins: { search_paths: { #Default plugin install directory installed_plugins: [ @CMAKE_INSTALL_PREFIX@/plugins ] } } data_dir: @CMAKE_INSTALL_PREFIX@/data persistency: { - default_output_file: nuis.out.root + default: { + output_file: nuis.out.root + } } global: { sample: { verbosity_default: Reticent } } diff --git a/data/CMakeLists.txt b/data/CMakeLists.txt new file mode 100644 index 0000000..aafcffb --- /dev/null +++ b/data/CMakeLists.txt @@ -0,0 +1,2 @@ +SET(DATA_INSTALL_DIR data) +add_subdirectory(nuA) diff --git a/data/nuA/BubbleChamber/ANL/CCQE/ANL_CCQE_Data_PRD16_3103.root b/data/nuA/BubbleChamber/ANL/CCQE/ANL_CCQE_Data_PRD16_3103.root new file mode 100644 index 0000000..e9af236 Binary files /dev/null and b/data/nuA/BubbleChamber/ANL/CCQE/ANL_CCQE_Data_PRD16_3103.root differ diff --git a/data/nuA/BubbleChamber/ANL/CCQE/ANL_CCQE_Data_PRD26_537.root b/data/nuA/BubbleChamber/ANL/CCQE/ANL_CCQE_Data_PRD26_537.root new file mode 100644 index 0000000..965ec79 Binary files /dev/null and b/data/nuA/BubbleChamber/ANL/CCQE/ANL_CCQE_Data_PRD26_537.root differ diff --git a/data/nuA/BubbleChamber/ANL/CCQE/ANL_CCQE_Data_PRL31_844.root b/data/nuA/BubbleChamber/ANL/CCQE/ANL_CCQE_Data_PRL31_844.root new file mode 100644 index 0000000..11e20a0 Binary files /dev/null and b/data/nuA/BubbleChamber/ANL/CCQE/ANL_CCQE_Data_PRL31_844.root differ diff --git a/data/nuA/BubbleChamber/ANL/CCQE/CMakeLists.txt b/data/nuA/BubbleChamber/ANL/CCQE/CMakeLists.txt new file mode 100644 index 0000000..b71f174 --- /dev/null +++ b/data/nuA/BubbleChamber/ANL/CCQE/CMakeLists.txt @@ -0,0 +1,5 @@ +SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/CCQE) + +FILE(GLOB DATA_FILES *.root) + +install(FILES ${DATA_FILES} DESTINATION ${DATA_INSTALL_DIR}) diff --git a/data/nuA/BubbleChamber/ANL/CMakeLists.txt b/data/nuA/BubbleChamber/ANL/CMakeLists.txt new file mode 100644 index 0000000..0d377cb --- /dev/null +++ b/data/nuA/BubbleChamber/ANL/CMakeLists.txt @@ -0,0 +1,3 @@ +SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/ANL) + +add_subdirectory(CCQE) diff --git a/data/nuA/BubbleChamber/CMakeLists.txt b/data/nuA/BubbleChamber/CMakeLists.txt new file mode 100644 index 0000000..e4bba08 --- /dev/null +++ b/data/nuA/BubbleChamber/CMakeLists.txt @@ -0,0 +1,3 @@ +SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/BubbleChamber) + +add_subdirectory(ANL) diff --git a/data/nuA/CMakeLists.txt b/data/nuA/CMakeLists.txt new file mode 100644 index 0000000..36a3bdc --- /dev/null +++ b/data/nuA/CMakeLists.txt @@ -0,0 +1,4 @@ +SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/nuA) + +add_subdirectory(BubbleChamber) +add_subdirectory(Nuclear) diff --git a/data/nuA/Nuclear/CMakeLists.txt b/data/nuA/Nuclear/CMakeLists.txt new file mode 100644 index 0000000..b7a435f --- /dev/null +++ b/data/nuA/Nuclear/CMakeLists.txt @@ -0,0 +1,3 @@ +SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/Nuclear) + +add_subdirectory(T2K) diff --git a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt b/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt new file mode 100644 index 0000000..18cc6d1 --- /dev/null +++ b/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt @@ -0,0 +1,5 @@ +SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/CC0Pi) + +FILE(GLOB DATA_FILES *.root) + +install(FILES ${DATA_FILES} DESTINATION ${DATA_INSTALL_DIR}) diff --git a/data/nuA/Nuclear/T2K/CC0Pi/H2O_xsec_2Dpmuthetamu_numubar.root b/data/nuA/Nuclear/T2K/CC0Pi/H2O_xsec_2Dpmuthetamu_numubar.root new file mode 100644 index 0000000..cae42cc Binary files /dev/null and b/data/nuA/Nuclear/T2K/CC0Pi/H2O_xsec_2Dpmuthetamu_numubar.root differ diff --git a/data/nuA/Nuclear/T2K/CMakeLists.txt b/data/nuA/Nuclear/T2K/CMakeLists.txt new file mode 100644 index 0000000..25096b5 --- /dev/null +++ b/data/nuA/Nuclear/T2K/CMakeLists.txt @@ -0,0 +1,3 @@ +SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/T2K) + +add_subdirectory(CC0Pi) diff --git a/doc/design_choices.md b/doc/design_choices.md index b688f6a..86ad1f7 100644 --- a/doc/design_choices.md +++ b/doc/design_choices.md @@ -1,30 +1,30 @@ # Design Choices This section describes the guiding design choices that are attempted to be implemented in the NUISANCE v3 rewrite. Any perceived failure to adhere to these design choices should be regarded as a bug and submitted to the developers. ## Modularity of components The various moving parts that go into a NUISANCE analysis should be as independent as possible (and no more). For example, as in previous versions, 'input handlers' are implementations of an `InputHandler` ABC, they are responsible for converting a specific event format into the NUISANCE event format. This event format should be generator agnostic as far as possible, however, for the generator-based-event reweighting to perform efficiently, the original generator event must be associated to the NUISANCE event. Instead of having some per-generator event manager, or reweight engines being allowed specific access to known InputHandler subclasses, it was decided that a pointer to the generator event would be carried around by the NUISANCE event base-class. This design choice has been kept. -In NUISANCE v3, the design of the 'sample', the basic implementation unit of a comparison, has been further modularized. The previous Measurement{Base,1D,2D} base classes have been surplanted by the ISample and IDataComparison interfaces. IDataComparison being a subclass of ISample. These specify the interface that any NUISANCE sample must implement. The implementation of any subclass is left up to the user, however, a helper base class that should be used for almost all simple data comparisons can be found in SimpleDataComparison. It attempts to provide a fully functional base class, for data comparisons that previously subclassed from Measurement1D or 2D, that needs minimal specializing for specific data samples---most importantly the signal definition, the kinematic projection, and optional-but-encouraged data set metadata. As the interface suggests, the responsibility of looping over events and determining event variations is now fully that of the ISample, rather than a finely tuned collaboration between the sample and some calling method. However, subclasses of SimpleDataComparison can still be fully implemented by providing only per-event methods as the base class handles the event looping. +In NUISANCE v3, the design of the 'sample', the basic implementation unit of a comparison, has been further modularized. The previous Measurement{Base,1D,2D} base classes have been surplanted by the IEventProcessor and IDataComparison interfaces. IDataComparison being a subclass of IEventProcessor. These specify the interface that any NUISANCE sample must implement. The implementation of any subclass is left up to the user, however, a helper base class that should be used for almost all simple data comparisons can be found in SimpleDataComparison. It attempts to provide a fully functional base class, for data comparisons that previously subclassed from Measurement1D or 2D, that needs minimal specializing for specific data samples---most importantly the signal definition, the kinematic projection, and optional-but-encouraged data set metadata. As the interface suggests, the responsibility of looping over events and determining event variations is now fully that of the IEventProcessor, rather than a finely tuned collaboration between the sample and some calling method. However, subclasses of SimpleDataComparison can still be fully implemented by providing only per-event methods as the base class handles the event looping. Related to this, it may be useful to extend the 'reweight' concept, which was assumed to be the only form of event response in previous versions to a more general event 'variation'. However, since the majority of variations will still be weight-based, and as weight-based variations can often be calculated and applied fully with just the `nuis::event::MinimalEvent` format, the InputHandler will expose convenience methods that allow the current event weight to be without direct interrogation of any IWeightProvider instances. However, because general event variations can be arbitrarily complex and can require signal definitions to be rechecked (requiring a `nuis::event::FullEvent`), no such convenience method exists and samples looking to make use of such variations must explicitly call the VariationManager. ## Extensibility without recompilation In previous versions of NUISANCE, the addition of new studies by non-experts was somewhat involved. Even though well-documented in multiple tutorials the addition of new studies required multiple source code changes, for opaque reasons outside of the study implementation itself. It also required recompilation of the entire set of NUISANCE binaries, this makes sharing a NUISANCE release on a cluster problematic and limits the entry barrier for non-expert users. V3 of NUISANCE is built from the ground up to allow for use and extensibility by non-experts, this is primarily facilitated by the `Instantiate` template method that implements a plugin factory, somewhat inspired by the `art::make_tool` utility from the ART framework. While the architecture of c++ means that plugins compiled with different toolsets than the main NUISANCE binaries are unlikely to work (and thus checked and disallowed before attempting any c++ object instantiation), multiple users working on a cluster should now be able to have local analysis implementations that are dynamically instantiated at runtime and require no recompilation of a central NUISANCE install. ## Flexibility of configuration Originally NUISANCE used an inflexible DSL for configuration, in v2, a more-flexible XML-based configuration file was used. The change increased the scope of specific extensions to the core configuration, but the configuration API was clunky. In v3, NUISANCE uses an implementation of the FHiCL-c++ language bindings to provide both configuration parsing and validation, as well as programatic access to the runtime configuration. FHiCL is particularly adept at providing a naturally hierarchial configuration, whereby global defaults can be specified but overridden at runtime by user-supplied configuration files without special handling of the in-memory global configuration document. The FHiCL language bindings also specify how to retrieve configuration elements with specified types, reducing the need for comprehensive runtime checks for missing or incomplete configuration. ## Exceptional circumstances -Previously, when NUISANCE encountered exceptional circumstances an exception was thrown. This exception was usually unqualified and sometimes preceeded by an error message. This made debugging and reporting of runtime errors fiddly for non-experts. In v3, all exceptions should be sensibly named, derived from `nuis::nuis_except` and be supplied with a verbose message diagnosing the problem encountered. +Previously, when NUISANCE encountered exceptional circumstances an exception was thrown. This exception was usually unqualified and sometimes preceded by an error message. This made debugging and reporting of runtime errors fiddly for non-experts. In v3, all exceptions should be sensibly named, derived from `nuis::nuis_except` and be supplied with a verbose message diagnosing the problem encountered. ## Documentation -The inline documentation of the code in previous versions of NUISANCE has been acceptable. The barrier to entry for new users was significantly lowered by external tutorials and documents. For v3, the hope is that the interface documentation can be fully comprehensive, and a specific set of instructions for common and less-common tasks can be distributed with the code. In addition to this, applications and scripts that verbosely aid users in standard NUISANCE workflows should be provided, and all usage prompts should be kept informative and fully up-to-date. +The inline documentation of the code in previous versions of NUISANCE has been acceptable. The barrier to entry for new users was significantly lowered by external tutorials and documents. For v3, the hope is that the interface documentation can be fully comprehensive, and a specific set of instructions for common and less-common tasks can be distributed with the code. In addition to this, applications and scripts that verbosely aid users in standard NUISANCE workflows should be provided, and all usage prompts should be kept informative and fully up-to-date. Plugins should be queryable for example configurations. diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b98e8d5..b615e75 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,13 +1,15 @@ add_subdirectory(config) add_subdirectory(event) add_subdirectory(input) add_subdirectory(plugins) add_subdirectory(utility) add_subdirectory(generator) add_subdirectory(persistency) +add_subdirectory(parameters) +add_subdirectory(variation) add_subdirectory(samples) SET(INuADataComparisons_List ${INuADataComparisons_List} PARENT_SCOPE) SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE) add_subdirectory(app) diff --git a/src/app/CMakeLists.txt b/src/app/CMakeLists.txt index 3f12000..f11b727 100644 --- a/src/app/CMakeLists.txt +++ b/src/app/CMakeLists.txt @@ -1,13 +1,13 @@ SET(APPS nuissamples nuiscomp nuisevsum) foreach(a ${APPS}) add_executable(${a} ${a}.cxx) - target_link_libraries(${a} nuis_event nuis_input nuis_utility nuis_config nuis_persistency) + target_link_libraries(${a} nuis_event nuis_input nuis_utility nuis_config nuis_persistency nuis_plugins nuis_params nuis_variation) target_link_libraries(${a} ${CMAKE_DEPENDLIB_FLAGS}) if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "") set_target_properties(${a} PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS}) endif() install(TARGETS ${a} DESTINATION bin) endforeach() diff --git a/src/app/nuiscomp.cxx b/src/app/nuiscomp.cxx index da163f3..c69028f 100644 --- a/src/app/nuiscomp.cxx +++ b/src/app/nuiscomp.cxx @@ -1,53 +1,57 @@ #include "config/GlobalConfiguration.hxx" #include "input/IInputHandler.hxx" #include "event/MinimalEvent.hxx" #include "samples/IDataComparison.hxx" #include "plugins/Instantiate.hxx" #include "exception/exception.hxx" +#include "persistency/ROOTOutput.hxx" + #include "fhiclcpp/make_ParameterSet.h" #include NEW_NUIS_EXCEPT(invalid_cli_arguments); int main(int argc, char const *argv[]) { nuis::config::EnsureConfigurationRead("nuis.global.config.fcl"); nuis::config::EnsureConfigurationRead("nuis.datacomparisons.fcl"); if (argc != 2) { throw invalid_cli_arguments() << "[ERROR]: Expected to be passed a single FHiCL file name or " "absolute or relative path. N.B. Files in the local directory must " "be fully qualified like \"$ " << argv[0] << " ./myconf.fcl\"."; } nuis::config::EnsureConfigurationRead(argv[1]); size_t NMax = nuis::config::GetDocument().get( "nmax", std::numeric_limits::max()); for (fhicl::ParameterSet const &samp_config : nuis::config::GetDocument().get>( "samples")) { std::cout << "[INFO]: Reading sample: " << samp_config.get("name") << std::endl; nuis::plugins::plugin_traits::unique_ptr_t sample = nuis::plugins::Instantiate( samp_config.get("name")); sample->Initialize(samp_config); sample->ProcessSample(NMax); std::cout << "[INFO]:\t Sample GOF = " << sample->GetGOF() << " / " << sample->GetNDOGuess() << std::endl; sample->Write(); } + + nuis::persistency::CloseOpenTFiles(); } diff --git a/src/app/nuisevsum.cxx b/src/app/nuisevsum.cxx index c45f677..5461e22 100644 --- a/src/app/nuisevsum.cxx +++ b/src/app/nuisevsum.cxx @@ -1,79 +1,79 @@ #include "config/GlobalConfiguration.hxx" #include "input/IInputHandler.hxx" #include "event/MinimalEvent.hxx" -#include "samples/ISample.hxx" +#include "samples/IEventProcessor.hxx" #include "plugins/Instantiate.hxx" #include "exception/exception.hxx" #include "fhiclcpp/make_ParameterSet.h" #include "string_parsers/from_string.hxx" #include NEW_NUIS_EXCEPT(invalid_cli_arguments); size_t NMax = std::numeric_limits::max(); std::string input_file; std::string input_type; void SayUsage(char const *argv[]) { std::cout << "[USAGE]: " << argv[0] << "\n" "\t-i : Input file passed to named " "IInputHandler instance \n" "\t-H : Name of IInputHandler subclass " "capable of reading NUISANCE events from the argument of -i.\n" "\t-n : Maximum number of events to " "read. Will read entire input file by default.\n" << std::endl; } void handleOpts(int argc, char const *argv[]) { int opt = 1; while (opt < argc) { if ((std::string(argv[opt]) == "-?") || (std::string(argv[opt]) == "--help")) { SayUsage(argv); exit(0); } else if (std::string(argv[opt]) == "-i") { input_file = argv[++opt]; } else if (std::string(argv[opt]) == "-H") { input_type = argv[++opt]; } else if (std::string(argv[opt]) == "-n") { NMax = fhicl::string_parsers::str2T(argv[++opt]); } else { std::cout << "[ERROR]: Unknown option: " << argv[opt] << std::endl; SayUsage(argv); exit(1); } opt++; } } int main(int argc, char const *argv[]) { nuis::config::EnsureConfigurationRead("nuis.global.config.fcl"); handleOpts(argc, argv); if (!input_type.length() || !input_file.length()) { SayUsage(argv); throw invalid_cli_arguments() << "[ERROR]: Require both -i and -H cli options to be passed."; } fhicl::ParameterSet sample_config; sample_config.put("input_type", input_type); sample_config.put("file", input_file); - nuis::plugins::plugin_traits::unique_ptr_t VerboseEventSummary = - nuis::plugins::Instantiate("VerboseEventSummary"); + nuis::plugins::plugin_traits::unique_ptr_t VerboseEventSummary = + nuis::plugins::Instantiate("VerboseEventSummary"); VerboseEventSummary->Initialize(sample_config); VerboseEventSummary->ProcessSample(NMax); } diff --git a/src/app/nuissamples.cxx b/src/app/nuissamples.cxx index 3fdf60d..d3beeb5 100644 --- a/src/app/nuissamples.cxx +++ b/src/app/nuissamples.cxx @@ -1,111 +1,115 @@ #include "config/GlobalConfiguration.hxx" #include "input/IInputHandler.hxx" #include "plugins/Instantiate.hxx" #include "samples/IDataComparison.hxx" #include "utility/StringUtility.hxx" #include #include #include #include std::string search_term = ".*"; bool fStrictRegex = false; std::string config_out_filename = ""; void SayUsage(char const *argv[]) { std::cout << "[USAGE]: " << argv[0] << "\n" "\t-s : Used to filter known IDataComparisons. " "\n\t\tWild cards are added to either side of the search term.\n" "\t-S : Used to filter known IDataComparisons. " "\n\t\tThe exact passed term is used.\n" "\t-o : Dump example sample configuration file\n" "\t\tfor matching samples.\n" << std::endl; } void handleOpts(int argc, char const *argv[]) { int opt = 1; while (opt < argc) { if ((std::string(argv[opt]) == "-?") || (std::string(argv[opt]) == "--help")) { SayUsage(argv); exit(0); } else if (std::string(argv[opt]) == "-s") { search_term = argv[++opt]; } else if (std::string(argv[opt]) == "-S") { fStrictRegex = true; search_term = argv[++opt]; } else if (std::string(argv[opt]) == "-o") { config_out_filename = argv[++opt]; } else { std::cout << "[ERROR]: Unknown option: " << argv[opt] << std::endl; SayUsage(argv); exit(1); } opt++; } } NEW_NUIS_EXCEPT(invalid_output_file); int main(int argc, char const *argv[]) { nuis::config::EnsureConfigurationRead("nuis.global.config.fcl"); nuis::config::EnsureConfigurationRead("nuis.datacomparisons.fcl"); handleOpts(argc, argv); std::regex rpattern(fStrictRegex ? search_term : std::string(".*") + search_term + ".*"); std::vector example_sample_configs; for (std::string const &comparison_set_key : nuis::config::GetDocument() .get("data_comparisons") .get_names()) { for (std::string const &sample_name : nuis::config::GetDocument().get>( std::string("data_comparisons.") + comparison_set_key)) { if (!std::regex_match(sample_name, rpattern)) { continue; } nuis::plugins::plugin_traits::unique_ptr_t sample = nuis::plugins::Instantiate(sample_name); std::cout << sample->Name() << std::endl; std::cout << "\tJournal: " << sample->GetJournalReference() << std::endl; std::cout << "\tTarget: " << sample->GetTargetMaterial() << std::endl; std::cout << "\tFlux: " << sample->GetFluxDescription() << std::endl; std::cout << "\tSignal: " << sample->GetSignalDescription() << std::endl; std::cout << "\tDocs: \n" << nuis::utility::indent_apply_width(sample->GetDocumentation(), 10) << std::endl; - std::cout << "Config: " - << sample->GetExampleConfiguration().to_indented_string() + std::cout << "\tExample_Config: {\n" + << nuis::utility::indent_apply_width( + sample->GetExampleConfiguration().to_indented_string(), + 12) + << "\n\t}\n" << std::endl; + ; example_sample_configs.push_back(sample->GetExampleConfiguration()); } } if (config_out_filename.length()) { std::ofstream out_file(config_out_filename); if (!out_file) { throw invalid_output_file() << "[ERROR]: Failed to open output file: " << std::quoted(config_out_filename); } fhicl::ParameterSet example_config; example_config.put("samples", example_sample_configs); out_file << example_config.to_indented_string(); } } diff --git a/src/config/GlobalConfiguration.cxx b/src/config/GlobalConfiguration.cxx index ec35a15..772b429 100644 --- a/src/config/GlobalConfiguration.cxx +++ b/src/config/GlobalConfiguration.cxx @@ -1,41 +1,41 @@ #include "config/GlobalConfiguration.hxx" #include "fhiclcpp/ParameterSet.h" #include "fhiclcpp/make_ParameterSet.h" namespace nuis { namespace config { struct NamedDoc { std::string name; fhicl::ParameterSet document; std::vector files_read; }; std::vector Configurations; NamedDoc &GetDocument_modifyable(std::string const &name) { for (NamedDoc &doc : Configurations) { if (doc.name == name) { return doc; } } - Configurations.push_back(NamedDoc{name, fhicl::ParameterSet()}); + Configurations.push_back(NamedDoc{name, fhicl::ParameterSet(), {}}); return GetDocument_modifyable(name); } fhicl::ParameterSet const &GetDocument(std::string const &name) { return GetDocument_modifyable(name).document; } bool EnsureConfigurationRead(std::string const &fhicl_file, std::string const &name) { NamedDoc &doc = GetDocument_modifyable(name); if (std::count(doc.files_read.begin(), doc.files_read.end(), fhicl_file)) { return false; } doc.document.splice(fhicl::make_ParameterSet(fhicl_file)); return true; } } // namespace config } // namespace nuis diff --git a/src/event/FullEvent.cxx b/src/event/FullEvent.cxx index 8d77987..c058989 100644 --- a/src/event/FullEvent.cxx +++ b/src/event/FullEvent.cxx @@ -1,24 +1,54 @@ #include "event/FullEvent.hxx" namespace nuis { namespace event { FullEvent::FullEvent() : MinimalEvent() { for (size_t status_it = 0; status_it < static_cast(Particle::Status_t::kNParticleStatus); ++status_it) { ParticleStack.push_back({static_cast(status_it), {}}); } } FullEvent::FullEvent(FullEvent &&other) : MinimalEvent(std::move(other)), ParticleStack(std::move(other.ParticleStack)) {} +FullEvent &FullEvent::operator=(FullEvent &&other) { + MinimalEvent::operator=(std::move(other)); + ParticleStack = std::move(other.ParticleStack); + return *this; +} + +FullEvent FullEvent::Clone() const { + FullEvent clone; + clone.MinimalEvent::operator=(MinimalEvent::Clone()); + clone.ParticleStack = ParticleStack; + return clone; +} + void FullEvent::ClearParticleStack() { for (auto &status_stack : ParticleStack) { status_stack.particles.clear(); } } -} // namespace core +std::string FullEvent::to_string() const { + std::stringstream ss(""); + ss << "Event: Interaction mode = " << mode << ", probe: { PDG: " << probe_pdg + << ", Energy: " << probe_E << " MeV }." << std::endl; + for (auto &status_stack : ParticleStack) { + ss << "\t[" << status_stack.status << "]" << std::endl; + + for (Particle const &part : status_stack.particles) { + ss << "\t\t{ PDG: " << part.pdg << ", P3: [ " << part.P4[0] << ", " + << part.P4[1] << ", " << part.P4[2] << "], E: " << part.P4[3] + << ", M: " << part.P4.M() << " }" << std::endl; + } + } + ss << std::endl; + return ss.str(); +} + +} // namespace event } // namespace nuis diff --git a/src/event/FullEvent.hxx b/src/event/FullEvent.hxx index 74e4662..6e18f1a 100644 --- a/src/event/FullEvent.hxx +++ b/src/event/FullEvent.hxx @@ -1,70 +1,58 @@ // Copyright 2018 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 EVENT_FULLEVENT_HXX_SEEN #define EVENT_FULLEVENT_HXX_SEEN #include "event/MinimalEvent.hxx" #include "event/Particle.hxx" #include "string_parsers/to_string.hxx" #include namespace nuis { namespace event { ///\brief The full, internal event format. class FullEvent : public MinimalEvent { public: struct StatusParticles { Particle::Status_t status; std::vector particles; }; FullEvent(); FullEvent(FullEvent const &) = delete; FullEvent(FullEvent &&); + FullEvent& operator=(FullEvent &&); + + FullEvent Clone() const; + std::vector ParticleStack; void ClearParticleStack(); - std::string to_string() const { - std::stringstream ss(""); - ss << "Event: Interaction mode = " << mode - << ", probe: { PDG: " << probe_pdg << ", Energy: " << probe_E - << " MeV }." << std::endl; - for (auto &status_stack : ParticleStack) { - ss << "\t[" << status_stack.status << "]" << std::endl; - - for (Particle const &part : status_stack.particles) { - ss << "\t\t{ PDG: " << part.pdg << ", P3: [ " << part.P4[0] << ", " - << part.P4[1] << ", " << part.P4[2] << "], E: " << part.P4[3] - << ", M: " << part.P4.M() << " }" << std::endl; - } - } - ss << std::endl; - return ss.str(); - } + std::string to_string() const; }; } // namespace event } // namespace nuis #endif diff --git a/src/event/MinimalEvent.cxx b/src/event/MinimalEvent.cxx index 1e78ed6..875b99f 100644 --- a/src/event/MinimalEvent.cxx +++ b/src/event/MinimalEvent.cxx @@ -1,22 +1,49 @@ #include "event/MinimalEvent.hxx" namespace nuis { namespace event { MinimalEvent::MinimalEvent() : mode(Channel_t::kUndefined), probe_E(0), probe_pdg(0), XSecWeight(1), RWWeight(1) { #ifdef __NUWRO_ENABLED__ fNuWroEvent = nullptr; #endif } MinimalEvent::MinimalEvent(MinimalEvent &&other) : mode(other.mode), probe_E(other.probe_E), probe_pdg(other.probe_pdg), XSecWeight(other.XSecWeight), RWWeight(other.RWWeight) { #ifdef __NUWRO_ENABLED__ fNuWroEvent = other.fNuWroEvent; other.fNuWroEvent = nullptr; #endif } -} // namespace core + +MinimalEvent &MinimalEvent::operator=(MinimalEvent &&other) { + mode = other.mode; + probe_E = other.probe_E; + probe_pdg = other.probe_pdg; + XSecWeight = other.XSecWeight; + RWWeight = other.RWWeight; +#ifdef __NUWRO_ENABLED__ + fNuWroEvent = other.fNuWroEvent; + other.fNuWroEvent = nullptr; +#endif + return *this; +} + +MinimalEvent MinimalEvent::Clone() const { + MinimalEvent clone; + clone.mode = mode; + clone.probe_E = probe_E; + clone.probe_pdg = probe_pdg; + clone.XSecWeight = XSecWeight; + clone.RWWeight = RWWeight; +#ifdef __NUWRO_ENABLED__ + clone.fNuWroEvent = fNuWroEvent; +#endif + + return clone; +} +} // namespace event } // namespace nuis diff --git a/src/event/MinimalEvent.hxx b/src/event/MinimalEvent.hxx index 328965d..83f36f1 100644 --- a/src/event/MinimalEvent.hxx +++ b/src/event/MinimalEvent.hxx @@ -1,73 +1,78 @@ // Copyright 2018 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 EVENT_MINIMALEVENT_HXX_SEEN #define EVENT_MINIMALEVENT_HXX_SEEN #ifdef NUWRO_ENABLED #include "event1.h" typedef ::event NuWroEvent; #endif #ifdef NEUT_ENABLED #include "neutpart.h" #include "neutvect.h" #endif #include "event/types.hxx" namespace nuis { namespace event { ///\brief The minimal event information needed to perform reweights. /// /// Most often, event selections cannot be applied using this reduced format. class MinimalEvent { public: MinimalEvent(); MinimalEvent(MinimalEvent const &) = delete; MinimalEvent(MinimalEvent &&); + MinimalEvent &operator=(MinimalEvent &&); + + /// Make a clone of this MinimalEvent + MinimalEvent Clone() const; + /// True interaction mode Channel_t mode; /// True probe energy double probe_E; /// True probe particle code PDG_t probe_pdg; PDG_t target_pdg; /// Event-weight that can be used to scale to a cross-section prediction double XSecWeight; /// Event weight incurred from current reweight engine state. double RWWeight; #ifdef NUWRO_ENABLED ///\brief Pointer to Nuwro event /// /// This will usually be tied to a TTree and so we are not responsible for /// deleting it NuWroEvent *fNuWroEvent; #endif #ifdef NEUT_ENABLED NeutVect *fNeutVect; #endif }; } // namespace core } // namespace nuis #endif diff --git a/src/event/Particle.hxx b/src/event/Particle.hxx index 0bcad11..dc6f4c0 100644 --- a/src/event/Particle.hxx +++ b/src/event/Particle.hxx @@ -1,77 +1,77 @@ // Copyright 2018 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 EVENT_PARTICLE_HXX_SEEN #define EVENT_PARTICLE_HXX_SEEN #include "event/types.hxx" #include "TLorentzVector.h" namespace nuis { namespace event { class Particle { public: NEW_NUIS_EXCEPT(invalid_particle); #define STATUS_LIST \ X(kNuclearLeaving, 0) \ X(kPrimaryInitialState, 1) \ X(kPrimaryFinalState, 2) \ X(kIntermediate, 3) \ X(kUnknown, 4) \ X(kBlocked, 5) \ X(kNParticleStatus, 6) #define X(A, B) A = B, enum class Status_t { STATUS_LIST }; #undef X Particle(); Particle(Particle const &); PDG_t pdg; TLorentzVector P4; bool operator!() const { return (pdg == std::numeric_limits::max()); } double E() const { return P4.E(); } double P() const { return P4.Vect().Mag(); } TVector3 P3() const { return P4.Vect(); } double M() const { return P4.M(); } double Theta() const { return P4.Vect().Theta(); } - double CosTheta() const { return P4.Vect().Theta(); } + double CosTheta() const { return P4.Vect().CosTheta(); } }; } // namespace event } // namespace nuis #define X(A, B) \ case nuis::event::Particle::Status_t::A: { \ return os << #A; \ } inline std::ostream &operator<<(std::ostream &os, nuis::event::Particle::Status_t te) { switch (te) { STATUS_LIST } return os; } #undef X #undef STATUS_LIST #endif diff --git a/src/generator/CMakeLists.txt b/src/generator/CMakeLists.txt index 84e4aea..eda12fb 100644 --- a/src/generator/CMakeLists.txt +++ b/src/generator/CMakeLists.txt @@ -1,2 +1,3 @@ add_subdirectory(input) add_subdirectory(utility) +add_subdirectory(variation) diff --git a/src/generator/input/NEUTInputHandler.cxx b/src/generator/input/NEUTInputHandler.cxx index a4133b6..798299e 100644 --- a/src/generator/input/NEUTInputHandler.cxx +++ b/src/generator/input/NEUTInputHandler.cxx @@ -1,103 +1,119 @@ #include "generator/input/NEUTInputHandler.hxx" #include "utility/InteractionChannelUtility.hxx" #include "utility/PDGCodeUtility.hxx" #include "utility/ROOTUtility.hxx" #include "generator/utility/NEUTUtility.hxx" #include "fhiclcpp/ParameterSet.h" using namespace nuis::event; using namespace nuis::utility; using namespace nuis::neuttools; -NEUTInputHandler::NEUTInputHandler() : fInputTree(nullptr) {} +NEUTInputHandler::NEUTInputHandler() : fInputTreeFile(nullptr) {} NEUTInputHandler::NEUTInputHandler(NEUTInputHandler &&other) - : fInputTree(std::move(other.fInputTree)), + : fInputTreeFile(std::move(other.fInputTreeFile)), fReaderEvent(std::move(other.fReaderEvent)) {} void NEUTInputHandler::Initialize(fhicl::ParameterSet const &ps) { - fInputTree = CheckGetTTree(ps.get("file"), "neuttree"); + fInputTreeFile = CheckGetTTree(ps.get("file"), "neuttree"); fReaderEvent.fNeutVect = nullptr; - fInputTree->tree->SetBranchAddress("vectorbranch", &fReaderEvent.fNeutVect); + fInputTreeFile->tree->SetBranchAddress("vectorbranch", + &fReaderEvent.fNeutVect); fKeepIntermediates = ps.get("keep_intermediates", false); } +double NEUTInputHandler::GetXSecScaleFactor(std::pair const &) const { + // check input file for histograms + std::unique_ptr flux = + GetHistogramFromROOTFile(fInputTreeFile->file, "flux_numu"); + std::unique_ptr evtrt = + GetHistogramFromROOTFile(fInputTreeFile->file, "evtrt_numu"); + + if (flux && evtrt) { + return evtrt->Integral() / (flux->Integral() * double(GetNEvents())); + } + + // try and build them + throw; +} + MinimalEvent const &NEUTInputHandler::GetMinimalEvent(ev_index_t idx) const { if (idx >= GetNEvents()) { throw IInputHandler::invalid_entry() << "[ERROR]: Attempted to get entry " << idx << " from an InputHandler with only " << GetNEvents(); } - fInputTree->tree->GetEntry(idx); + fInputTreeFile->tree->GetEntry(idx); fReaderEvent.mode = IntToChannel(fReaderEvent.fNeutVect->Mode); size_t NPart = fReaderEvent.fNeutVect->Npart(); for (size_t part_it = 0; part_it < NPart; part_it++) { NeutPart const &part = (*fReaderEvent.fNeutVect->PartInfo(part_it)); if ((part.fIsAlive == false) && (part.fStatus == -1) && IsNeutralLepton(part.fPID)) { fReaderEvent.probe_E = part.fP.T(); fReaderEvent.probe_pdg = part.fPID; break; } } return fReaderEvent; } FullEvent const &NEUTInputHandler::GetFullEvent(ev_index_t idx) const { (void)GetMinimalEvent(idx); fReaderEvent.ClearParticleStack(); if (fReaderEvent.fNeutVect->Ibound) { fReaderEvent.target_pdg = MakeNuclearPDG(fReaderEvent.fNeutVect->TargetA, fReaderEvent.fNeutVect->TargetZ); } else { fReaderEvent.target_pdg = MakeNuclearPDG(1, 1); } size_t NPart = fReaderEvent.fNeutVect->Npart(); size_t NPrimary = fReaderEvent.fNeutVect->Nprimary(); for (size_t part_it = 0; part_it < NPart; part_it++) { NeutPart const &part = (*fReaderEvent.fNeutVect->PartInfo(part_it)); Particle nuis_part; nuis_part.pdg = part.fPID; nuis_part.P4 = part.fP; Particle::Status_t state = GetNeutParticleStatus(part, fReaderEvent.mode); if (!fKeepIntermediates && (state == Particle::Status_t::kIntermediate)) { continue; } size_t state_int = static_cast(state); // Add status == 0 particles as pre-FSI particles until we find an // intermediate state particle if ((part_it < NPrimary) && (state != Particle::Status_t::kPrimaryInitialState)) { fReaderEvent .ParticleStack[static_cast( Particle::Status_t::kPrimaryFinalState)] .particles.push_back(nuis_part); } fReaderEvent.ParticleStack[state_int].particles.push_back(nuis_part); } return fReaderEvent; } size_t NEUTInputHandler::GetNEvents() const { - return fInputTree->tree->GetEntries(); + return fInputTreeFile->tree->GetEntries(); } DECLARE_PLUGIN(IInputHandler, NEUTInputHandler); diff --git a/src/generator/input/NEUTInputHandler.hxx b/src/generator/input/NEUTInputHandler.hxx index 9f4d8f6..4d063ea 100644 --- a/src/generator/input/NEUTInputHandler.hxx +++ b/src/generator/input/NEUTInputHandler.hxx @@ -1,56 +1,58 @@ // Copyright 2018 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 GENERATOR_INPUT_NEUTINPUTHANDLER_HXX_SEEN #define GENERATOR_INPUT_NEUTINPUTHANDLER_HXX_SEEN #include "event/FullEvent.hxx" #include "input/IInputHandler.hxx" #include namespace fhicl { class ParameterSet; } namespace nuis { namespace utility { class TreeFile; } } // namespace nuis class NEUTInputHandler : public IInputHandler { - mutable std::unique_ptr fInputTree; + mutable std::unique_ptr fInputTreeFile; mutable nuis::event::FullEvent fReaderEvent; bool fKeepIntermediates; public: NEUTInputHandler(); NEUTInputHandler(NEUTInputHandler const &) = delete; NEUTInputHandler(NEUTInputHandler &&); void Initialize(fhicl::ParameterSet const &); nuis::event::MinimalEvent const &GetMinimalEvent(ev_index_t idx) const; nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const; size_t GetNEvents() const; + double GetXSecScaleFactor( + std::pair const &EnuRange) const; }; #endif diff --git a/src/generator/input/NuWroInputHandler.cxx b/src/generator/input/NuWroInputHandler.cxx index 2fd3d52..52f2a7e 100644 --- a/src/generator/input/NuWroInputHandler.cxx +++ b/src/generator/input/NuWroInputHandler.cxx @@ -1,115 +1,121 @@ #include "generator/input/NuWroInputHandler.hxx" #include "generator/utility/NuWroUtility.hxx" #include "utility/ROOTUtility.hxx" #include "fhiclcpp/ParameterSet.h" #include "particle.h" -typedef ::particle NuWroParticle; +using NuWroParticle = ::particle; using namespace nuis::event; using namespace nuis::utility; using namespace nuis::nuwrotools; -NuWroInputHandler::NuWroInputHandler() : fInputTree(nullptr) {} +NuWroInputHandler::NuWroInputHandler() : fInputTree() {} NuWroInputHandler::NuWroInputHandler(NuWroInputHandler &&other) : fInputTree(std::move(other.fInputTree)), fReaderEvent(std::move(other.fReaderEvent)) {} void NuWroInputHandler::Initialize(fhicl::ParameterSet const &ps) { fInputTree = CheckGetTTree(ps.get("file"), "treeout"); fReaderEvent.fNuWroEvent = nullptr; - fInputTree->tree->SetBranchAddress("e", &fReaderEvent.fNuWroEvent); + fInputTree.tree->SetBranchAddress("e", &fReaderEvent.fNuWroEvent); fKeepIntermediates = ps.get("keep_intermediates", false); } MinimalEvent const &NuWroInputHandler::GetMinimalEvent(ev_index_t idx) const { if (idx >= GetNEvents()) { throw IInputHandler::invalid_entry() << "[ERROR]: Attempted to get entry " << idx << " from an InputHandler with only " << GetNEvents(); } - fInputTree->tree->GetEntry(idx); + fInputTree.tree->GetEntry(idx); fReaderEvent.mode = NuWroEventChannel(*fReaderEvent.fNuWroEvent); fReaderEvent.probe_E = fReaderEvent.fNuWroEvent->in[0].E(); fReaderEvent.probe_pdg = fReaderEvent.fNuWroEvent->in[0].pdg; fReaderEvent.XSecWeight = fReaderEvent.fNuWroEvent->weight / double(GetNEvents()); if (fWeightCache.size() <= idx) { fWeightCache.push_back(fReaderEvent.XSecWeight); } return fReaderEvent; } FullEvent const &NuWroInputHandler::GetFullEvent(ev_index_t idx) const { (void)GetMinimalEvent(idx); fReaderEvent.ClearParticleStack(); for (size_t p_it = 0; p_it < fReaderEvent.fNuWroEvent->in.size(); ++p_it) { NuWroParticle &part = fReaderEvent.fNuWroEvent->in[p_it]; Particle nuis_part; nuis_part.pdg = part.pdg; nuis_part.P4 = TLorentzVector(part[1], part[2], part[3], part[0]); fReaderEvent .ParticleStack[static_cast( Particle::Status_t::kPrimaryInitialState)] .particles.push_back(nuis_part); } for (size_t p_it = 0; p_it < fKeepIntermediates && fReaderEvent.fNuWroEvent->out.size(); ++p_it) { NuWroParticle &part = fReaderEvent.fNuWroEvent->out[p_it]; Particle nuis_part; nuis_part.pdg = part.pdg; nuis_part.P4 = TLorentzVector(part[1], part[2], part[3], part[0]); fReaderEvent .ParticleStack[static_cast( Particle::Status_t::kPrimaryFinalState)] .particles.push_back(nuis_part); } for (size_t p_it = 0; (p_it < fReaderEvent.fNuWroEvent->post.size()); ++p_it) { NuWroParticle &part = fReaderEvent.fNuWroEvent->post[p_it]; Particle nuis_part; nuis_part.pdg = part.pdg; nuis_part.P4 = TLorentzVector(part[1], part[2], part[3], part[0]); fReaderEvent .ParticleStack[static_cast(Particle::Status_t::kNuclearLeaving)] .particles.push_back(nuis_part); } return fReaderEvent; } double NuWroInputHandler::GetEventWeight(ev_index_t idx) const { if (idx > fWeightCache.size()) { throw weight_cache_miss() << "[ERROR]: Failed to get cached weight for event index: " << idx; } return fWeightCache[idx]; } +double NuWroInputHandler::GetXSecScaleFactor( + std::pair const &EnuRange) const { + throw input_handler_feature_unimplemented() + << "[ERROR]: Flux cuts not yet implemented for NuWro input handler."; +} + size_t NuWroInputHandler::GetNEvents() const { - return fInputTree->tree->GetEntries(); + return fInputTree.tree->GetEntries(); } DECLARE_PLUGIN(IInputHandler, NuWroInputHandler); diff --git a/src/generator/input/NuWroInputHandler.hxx b/src/generator/input/NuWroInputHandler.hxx index bb050e3..a7fc22c 100644 --- a/src/generator/input/NuWroInputHandler.hxx +++ b/src/generator/input/NuWroInputHandler.hxx @@ -1,63 +1,59 @@ // Copyright 2018 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 GENERATOR_INPUT_NUWROINPUTHANDLER_HXX_SEEN -#define GENERATOR_INPUT_NUWROINPUTHANDLER_HXX_SEEN +#pragma once #include "event/FullEvent.hxx" #include "input/IInputHandler.hxx" #include "exception/exception.hxx" +#include "utility/ROOTUtility.hxx" + #include namespace fhicl { class ParameterSet; } -namespace nuis { -namespace utility { -class TreeFile; -} -} // namespace nuis - class NuWroInputHandler : public IInputHandler { - mutable std::unique_ptr fInputTree; + mutable nuis::utility::TreeFile fInputTree; mutable nuis::event::FullEvent fReaderEvent; mutable std::vector fWeightCache; bool fKeepIntermediates; public: NEW_NUIS_EXCEPT(weight_cache_miss); NuWroInputHandler(); NuWroInputHandler(NuWroInputHandler const &) = delete; NuWroInputHandler(NuWroInputHandler &&); void Initialize(fhicl::ParameterSet const &); nuis::event::MinimalEvent const &GetMinimalEvent(ev_index_t idx) const; nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const; double GetEventWeight(ev_index_t idx) const; size_t GetNEvents() const; -}; -#endif + double GetXSecScaleFactor( + std::pair const &EnuRange) const; +}; diff --git a/src/generator/utility/CMakeLists.txt b/src/generator/utility/CMakeLists.txt index 7b60771..437ae87 100644 --- a/src/generator/utility/CMakeLists.txt +++ b/src/generator/utility/CMakeLists.txt @@ -1,20 +1,20 @@ if(USE_NuWro) LIST(APPEND GENERATOR_UTILS_IMPL NuWroUtility.cxx) LIST(APPEND GENERATOR_UTILS_HDR NuWroUtility.hxx) endif(USE_NuWro) if(USE_NEUT) LIST(APPEND GENERATOR_UTILS_IMPL NEUTUtility.cxx) LIST(APPEND GENERATOR_UTILS_HDR NEUTUtility.hxx) endif(USE_NEUT) if(GENERATOR_UTILS_IMPL) add_library(nuis_generator_utility SHARED ${GENERATOR_UTILS_IMPL}) target_link_libraries(nuis_generator_utility nuis_event) - install(TARGETS nuis_generator_utility DESTINATION lib) + install(TARGETS nuis_generator_utility DESTINATION plugins) endif() if(GENERATOR_UTILS_HDR) install(FILES ${GENERATOR_UTILS_HDR} DESTINATION include/generator/utility) endif() diff --git a/src/generator/variation/CMakeLists.txt b/src/generator/variation/CMakeLists.txt index e69de29..6cd1d8d 100644 --- a/src/generator/variation/CMakeLists.txt +++ b/src/generator/variation/CMakeLists.txt @@ -0,0 +1,13 @@ +LIST(APPEND WEIGHT_ENGINES_LINK_LIBS nuis_event nuis_config nuis_params) + +if(USE_NEUT) + LIST(APPEND WEIGHT_ENGINES_IMPL NEUTWeightEngine.cxx) + LIST(APPEND WEIGHT_ENGINES_LINK_LIBS nuis_generator_utility) +endif(USE_NEUT) + +if(WEIGHT_ENGINES_IMPL) + add_library(generator_weight_engines SHARED ${WEIGHT_ENGINES_IMPL}) + target_link_libraries(generator_weight_engines ${WEIGHT_ENGINES_LINK_LIBS}) + + install(TARGETS generator_weight_engines DESTINATION plugins) +endif() diff --git a/src/generator/variation/FillNEUTCommons.hxx b/src/generator/variation/FillNEUTCommons.hxx new file mode 100644 index 0000000..8240d31 --- /dev/null +++ b/src/generator/variation/FillNEUTCommons.hxx @@ -0,0 +1,159 @@ +#ifndef GENERATOR_VARIATION_FILLENEUTCOMMONS_HXX_SEEN +#define GENERATOR_VARIATION_FILLENEUTCOMMONS_HXX_SEEN + +#include "NFortFns.h" // Contains all the NEUT common blocks + +namespace NEUTUtils { + +inline void 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; +} + +} // namespace NEUTUtils + +#endif diff --git a/src/generator/variation/NEUTReWeight.cxx b/src/generator/variation/NEUTReWeight.cxx deleted file mode 100644 index 4e61a9c..0000000 --- a/src/generator/variation/NEUTReWeight.cxx +++ /dev/null @@ -1,3 +0,0 @@ -#include "generator/variation/NEUTReWeight.hxx" - -DECLARE_PLUGIN(IWeightProvider, NEUTReWeight); diff --git a/src/generator/variation/NEUTWeightEngine.cxx b/src/generator/variation/NEUTWeightEngine.cxx new file mode 100644 index 0000000..51c1977 --- /dev/null +++ b/src/generator/variation/NEUTWeightEngine.cxx @@ -0,0 +1,114 @@ +#include "generator/variation/NEUTWeightEngine.hxx" +#include "generator/variation/FillNEUTCommons.hxx" + +#include "fhiclcpp/ParameterSet.h" + +#include "event/MinimalEvent.hxx" + +// NEUT Engine includes +#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 "NSystUncertainty.h" + +#include + +using namespace nuis; +using namespace nuis::params; +using namespace nuis::event; + +void NEUTWeightEngine::Initialize(fhicl::ParameterSet const &ps) { + + fNeutRW = std::make_unique(); + TDirectory *dir = gDirectory; + fNeutRW->AdoptWghtCalc("xsec_ccqe", new neut::rew::NReWeightNuXSecCCQE); + fNeutRW->AdoptWghtCalc("xsec_res", new neut::rew::NReWeightNuXSecRES); + fNeutRW->AdoptWghtCalc("xsec_ccres", new neut::rew::NReWeightNuXSecCCRES); + fNeutRW->AdoptWghtCalc("xsec_coh", new neut::rew::NReWeightNuXSecCOH); + fNeutRW->AdoptWghtCalc("xsec_dis", new neut::rew::NReWeightNuXSecDIS); + fNeutRW->AdoptWghtCalc("xsec_ncel", new neut::rew::NReWeightNuXSecNCEL); + fNeutRW->AdoptWghtCalc("xsec_nc", new neut::rew::NReWeightNuXSecNC); + fNeutRW->AdoptWghtCalc("xsec_ncres", new neut::rew::NReWeightNuXSecNCRES); + fNeutRW->AdoptWghtCalc("nucl_casc", new neut::rew::NReWeightCasc); + fNeutRW->AdoptWghtCalc("nucl_piless", new neut::rew::NReWeightNuclPiless); + dir->cd(); + + for (fhicl::ParameterSet param_ps : + ps.get>("parameters")) { + + param_ps.put("type", GetName()); + + std::string const ¶m_name = param_ps.get("name"); + + NEUTSystParam nsp; + + nsp.nsyst = neut::rew::NSyst::FromString(param_name); + + if (nsp.nsyst == neut::rew::kNullSystematic) { + throw invalid_NEUT_syst_name() << "[ERROR]: NReWeight failed to parse " + << std::quoted(param_name) << " as a NEUT dial."; + } + fNeutRW->Systematics().Init(nsp.nsyst); + + nsp.pid = ParameterManager::Get().EnsureParameterRegistered(param_ps); + + fNEUTSysts.push_back(nsp); + } + + Reconfigure(); +} +void NEUTWeightEngine::Reconfigure() { + + for (NEUTSystParam const &nsp : fNEUTSysts) { + double val = ParameterManager::Get().GetParameterValue(nsp.pid); + fNeutRW->Systematics().Set(nsp.nsyst, val); + } + + fNeutRW->Reconfigure(); +} +double NEUTWeightEngine::GetEventWeight(nuis::event::MinimalEvent const &ev) { + if (!ev.fNeutVect) { + return 1.0; + } + NEUTUtils::FillNeutCommons(ev.fNeutVect); + return fNeutRW->CalcWeight(); +} + +std::string NEUTWeightEngine::GetName() { return "NEUTWeightEngine"; } +std::string NEUTWeightEngine::GetDocumentation() { return ""; } +fhicl::ParameterSet NEUTWeightEngine::GetExampleConfiguration() { + fhicl::ParameterSet ps; + + ps.put("weight_engine_name", GetName()); + fhicl::ParameterSet dial_maqe; + dial_maqe.put("name", "MAQE"); + + dial_maqe.put("start", 0); + dial_maqe.put("min", -3); + dial_maqe.put("max", 3); + dial_maqe.put("step", 0.1); + + fhicl::ParameterSet dial_mares; + dial_mares.put("name", "MARES"); + + dial_mares.put("start", 0); + dial_mares.put("min", -3); + dial_mares.put("max", 3); + dial_mares.put("step", 0.1); + + ps.put>( + "parameters", std::vector{{dial_maqe, dial_mares}}); + + return ps; +} + +DECLARE_PLUGIN(IWeightProvider, NEUTWeightEngine); diff --git a/src/generator/variation/NEUTReWeight.hxx b/src/generator/variation/NEUTWeightEngine.hxx similarity index 67% rename from src/generator/variation/NEUTReWeight.hxx rename to src/generator/variation/NEUTWeightEngine.hxx index fd2c330..b7581c8 100644 --- a/src/generator/variation/NEUTReWeight.hxx +++ b/src/generator/variation/NEUTWeightEngine.hxx @@ -1,35 +1,61 @@ // Copyright 2018 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 GENERATOR_VARIATION_NEUTREWEIGHT_HXX_SEEN #define GENERATOR_VARIATION_NEUTREWEIGHT_HXX_SEEN #include "variation/IWeightProvider.hxx" -class NEUTReWeight : public IWeightProvider { +#include "exception/exception.hxx" + +#include "parameters/ParameterManager.hxx" + +#include "NSyst.h" + +#include +#include + +namespace neut { +namespace rew { +class NReWeight; +} +} // namespace neut + +class NEUTWeightEngine : public IWeightProvider { + + struct NEUTSystParam { + nuis::params::paramId_t pid; + neut::rew::NSyst_t nsyst; + }; + + std::vector fNEUTSysts; + std::unique_ptr fNeutRW; + public: - double GetEventWeight(nuis::event::MinimalEvent const &); + NEW_NUIS_EXCEPT(invalid_NEUT_syst_name); + void Initialize(fhicl::ParameterSet const &); - paramId_t GetParameterId(std::string const &); - void SetParameterValue(paramId_t, double); - bool ParametersVaried(); void Reconfigure(); + double GetEventWeight(nuis::event::MinimalEvent const &); + std::string GetName(); + std::string GetDocumentation(); + fhicl::ParameterSet GetExampleConfiguration(); }; #endif diff --git a/src/input/CMakeLists.txt b/src/input/CMakeLists.txt index 5aad051..6036777 100644 --- a/src/input/CMakeLists.txt +++ b/src/input/CMakeLists.txt @@ -1,12 +1,12 @@ SET(input_implementation_files InputManager.cxx) SET(input_header_files InputManager.hxx IInputHandler.hxx) add_library(nuis_input SHARED ${input_implementation_files}) -target_link_libraries(nuis_input) +target_link_libraries(nuis_input nuis_plugins) install(TARGETS nuis_input DESTINATION lib) install(FILES ${input_header_files} DESTINATION include/input) diff --git a/src/input/IInputHandler.hxx b/src/input/IInputHandler.hxx index e558904..58ea6c0 100644 --- a/src/input/IInputHandler.hxx +++ b/src/input/IInputHandler.hxx @@ -1,113 +1,118 @@ // Copyright 2018 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 INPUT_IINPUTHANDLER_HXX_SEEN -#define INPUT_IINPUTHANDLER_HXX_SEEN +#pragma once #include "plugins/traits.hxx" #include "exception/exception.hxx" #include +#include namespace fhicl { class ParameterSet; } namespace nuis { namespace event { class MinimalEvent; class FullEvent; } // namespace event } // namespace nuis class IInputHandler { public: struct FullEvent_const_iterator : public std::iterator< std::input_iterator_tag, nuis::event::FullEvent const, size_t, nuis::event::FullEvent const *, nuis::event::FullEvent const &> { FullEvent_const_iterator(size_t _idx, IInputHandler const *_ih) { idx = _idx; ih = _ih; } FullEvent_const_iterator(FullEvent_const_iterator const &other) { idx = other.idx; ih = other.ih; } FullEvent_const_iterator operator=(FullEvent_const_iterator const &other) { idx = other.idx; ih = other.ih; return (*this); } bool operator==(FullEvent_const_iterator const &other) { return (idx == other.idx); } bool operator!=(FullEvent_const_iterator const &other) { return !(*this == other); } nuis::event::FullEvent const &operator*() { return ih->GetFullEvent(idx); } nuis::event::FullEvent const *operator->() { return &ih->GetFullEvent(idx); } FullEvent_const_iterator operator++() { idx++; return *this; } FullEvent_const_iterator operator++(int) { FullEvent_const_iterator tmp(*this); idx++; return tmp; } private: size_t idx; IInputHandler const *ih; }; NEW_NUIS_EXCEPT(invalid_input_file); NEW_NUIS_EXCEPT(invalid_entry); + NEW_NUIS_EXCEPT(input_handler_feature_unimplemented); typedef size_t ev_index_t; virtual void Initialize(fhicl::ParameterSet const &) = 0; virtual nuis::event::MinimalEvent const & GetMinimalEvent(ev_index_t idx) const = 0; virtual nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const = 0; virtual void RecalculateEventWeights(){}; - virtual double GetEventWeight(ev_index_t idx) const { return 1; }; + virtual double GetEventWeight(ev_index_t) const { return 1; }; + + /// Allows samples that implement flux cuts to request an appropriate + virtual double GetXSecScaleFactor( + std::pair const &EnuRange = std::pair{ + std::numeric_limits::max(), + std::numeric_limits::max()}) const = 0; virtual size_t GetNEvents() const = 0; FullEvent_const_iterator begin() const { return FullEvent_const_iterator(0, this); } FullEvent_const_iterator end() const { return FullEvent_const_iterator(GetNEvents(), this); } virtual ~IInputHandler() {} }; DECLARE_PLUGIN_INTERFACE(IInputHandler); - -#endif diff --git a/src/input/InputManager.cxx b/src/input/InputManager.cxx index 8928d84..6e2dbca 100644 --- a/src/input/InputManager.cxx +++ b/src/input/InputManager.cxx @@ -1,65 +1,66 @@ #include "input/InputManager.hxx" #include "plugins/Instantiate.hxx" #include "fhiclcpp/ParameterSet.h" #include namespace nuis { namespace input { InputManager *InputManager::_global_inst = nullptr; InputManager::NamedInputHandler::NamedInputHandler( std::string const &file, plugins::plugin_traits::unique_ptr_t &&IH) { name = file; handler = std::move(IH); } InputManager::InputManager() : Inputs() {} InputManager &InputManager::Get() { if (!_global_inst) { _global_inst = new InputManager(); } return *_global_inst; } InputManager::Input_id_t InputManager::EnsureInputLoaded(fhicl::ParameterSet const &ps) { std::string const &file_name = ps.get("file"); for (size_t i = 0; i < Inputs.size(); ++i) { if (Inputs[i].name == file_name) { return i; } } + Input_id_t iid = Inputs.size(); Inputs.emplace_back(file_name, nuis::plugins::Instantiate( ps.get("input_type") + "InputHandler")); Inputs.back().handler->Initialize(ps); - return (Inputs.size() - 1); + return iid; } InputManager::Input_id_t InputManager::GetInputId(std::string const &file_name) const { for (size_t i = 0; i < Inputs.size(); ++i) { if (Inputs[i].name == file_name) { return i; } } throw unknown_input() << "[ERROR]: Input file " << std::quoted(file_name) << " has not been loaded."; } IInputHandler const & InputManager::GetInputHandler(InputManager::Input_id_t id) const { if (id >= Inputs.size()) { throw unknown_input() << "[ERROR]: Attempted to get input with id " << id << ", but only have " << Inputs.size() << " loaded inputs."; } return *Inputs[id].handler.get(); } } // namespace input } // namespace nuis diff --git a/src/parameters/ParameterManager.cxx b/src/parameters/ParameterManager.cxx index 76b17d6..4f39e47 100644 --- a/src/parameters/ParameterManager.cxx +++ b/src/parameters/ParameterManager.cxx @@ -1,133 +1,146 @@ #include "parameters/ParameterManager.hxx" +#include "fhiclcpp/ParameterSet.h" + namespace nuis { namespace params { +ParameterManager::NamedParameter::NamedParameter() + : name(""), type(""), value(kDefaultValue), start(kDefaultValue), + min(kDefaultValue), max(kDefaultValue), step(kDefaultValue), + Penalty([](double) -> double { return 0; }) {} + +ParameterManager::NamedParameter::NamedParameter(NamedParameter &&other) + : name(std::move(other.name)), type(std::move(other.type)), + value(other.value), start(other.start), min(other.min), max(other.max), + step(other.step), Penalty(std::move(other.Penalty)) {} + ParameterManager *ParameterManager::_global_inst = nullptr; ParameterManager::ParameterManager() : locked(false) {} ParameterManager &ParameterManager::Get() { if (!_global_inst) { _global_inst = new ParameterManager(); } return *_global_inst; } void ParameterManager::ValidateParamId(paramId_t pid) { if (pid >= Parameters.size()) { throw invalid_parameter_id() << "[ERROR]: Passed parameter id " << pid << ", but the ParameterManager only knows about " << Parameters.size() << " parameters."; } } void ParameterManager::LockParameterList() { locked = true; } void ParameterManager::UnlockParameterList() { locked = false; } paramId_t ParameterManager::EnsureParameterRegistered(fhicl::ParameterSet const &ps) { if (locked) { throw parameter_list_is_locked() << "[ERROR]: Attempted to register parameter: " << ps.to_string() << " when global ParameterManager was locked in state: " << StateString(); } NamedParameter np; np.name = ps.get("name"); np.type = ps.get("type"); paramId_t pid = GetParameterId(np.name, np.type); if (pid != kParamUnhandled) { return pid; } np.start = ps.get("start"); np.value = np.start; np.min = ps.get("min", kDefaultLimit); np.max = ps.get("max", kDefaultLimit); np.step = ps.get("step"); pid = Parameters.size(); - Parameters.push_back(np); + Parameters.emplace_back(std::move(np)); return pid; } paramId_t ParameterManager::GetParameterId(std::string const &name, std::string const &type) { paramId_t pid = kParamUnhandled; for (size_t p_it = 0; p_it < Parameters.size(); ++p_it) { if (name != Parameters[p_it].name) { continue; } if (type.size() && (type != Parameters[p_it].type)) { continue; } // matches search, check if it is the first to match the search. if (pid != kParamUnhandled) { throw ambiguous_parameter_specified() << "[ERROR]: When searching for parameter by name-only, found at " "least two matching parameters: { PID: " << pid << ", name: " << Parameters[pid].name << ", type: " << Parameters[pid].type << " } and { PID: " << p_it << ", name: " << Parameters[p_it].name << ", type: " << Parameters[p_it].type << " }"; } pid = p_it; } return pid; } void ParameterManager::SetParameterValue(paramId_t pid, double val) { ValidateParamId(pid); if (!IsValidParameterValue(pid, val)) { throw param_value_out_of_bounds() - << "[ERROR]: Attempting to set parameter { PID: " << pid ", name: " - << Parameters[pid].name << ", type: " << Parameters[pid].type - << " } to " << val << ", but this is out of the allowed range [" + << "[ERROR]: Attempting to set parameter { PID: " << pid + << ", name: " << Parameters[pid].name + << ", type: " << Parameters[pid].type << " } to " << val + << ", but this is out of the allowed range [" << ((Parameters[pid].min == kDefaultLimit) ? "unbounded" : std::to_string(Parameters[pid].min)) << "," << ((Parameters[pid].max == kDefaultLimit) ? "unbounded" : std::to_string(Parameters[pid].max)) << "]"; } Parameters[pid].value = val; } double ParameterManager::GetParameterValue(paramId_t pid) { ValidateParamId(pid); return Parameters[pid].value; } double ParameterManager::GetParameterStep(paramId_t pid) { ValidateParamId(pid); return Parameters[pid].step; } double ParameterManager::GetParameterStart(paramId_t pid) { ValidateParamId(pid); return Parameters[pid].start; } double ParameterManager::GetParameterMin(paramId_t pid) { ValidateParamId(pid); return Parameters[pid].min; } double ParameterManager::GetParameterMax(paramId_t pid) { ValidateParamId(pid); return Parameters[pid].max; } bool ParameterManager::IsValidParameterValue(paramId_t pid, double val) { ValidateParamId(pid); return ( ((Parameters[pid].min == kDefaultLimit) || (Parameters[pid].min < val)) && ((Parameters[pid].max == kDefaultLimit) || (Parameters[pid].max > val))); } std::string ParameterManager::StateString() { return "Parameter Manager state:"; } } // namespace params } // namespace nuis diff --git a/src/parameters/ParameterManager.hxx b/src/parameters/ParameterManager.hxx index 07f8651..ae1de88 100644 --- a/src/parameters/ParameterManager.hxx +++ b/src/parameters/ParameterManager.hxx @@ -1,91 +1,97 @@ // Copyright 2018 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 PARAMETERS_PARAMETERMANAGER_HXX_SEEN #define PARAMETERS_PARAMETERMANAGER_HXX_SEEN #include "exception/exception.hxx" #include #include #include +#include namespace fhicl { class ParameterSet; } namespace nuis { namespace params { typedef size_t paramId_t; static paramId_t const kParamUnhandled = std::numeric_limits::max(); static double const kDefaultLimit = 0xdeadbeef; +static double const kDefaultValue = 0xdeadbeef; class ParameterManager { struct NamedParameter { + NamedParameter(); + NamedParameter(NamedParameter const &) = delete; + NamedParameter(NamedParameter &&); std::string name; std::string type; double value; double start; double min; double max; double step; + std::function Penalty; }; std::vector Parameters; // TMatrixD describing parameter covariance. ParameterManager(); static ParameterManager *_global_inst; void ValidateParamId(paramId_t); bool locked; public: static ParameterManager &Get(); NEW_NUIS_EXCEPT(invalid_parameter_id); NEW_NUIS_EXCEPT(param_value_out_of_bounds); NEW_NUIS_EXCEPT(ambiguous_parameter_specified); NEW_NUIS_EXCEPT(parameter_list_is_locked); ///\brief Lock the parameter list so that subsequent calls to ///EnsureParameterRegistered cause an exception to be thrown. /// /// Useful for ensuring that plugins do not attempt to add parameters mid-fit. void LockParameterList(); void UnlockParameterList(); paramId_t EnsureParameterRegistered(fhicl::ParameterSet const &); paramId_t GetParameterId(std::string const &, std::string const &type = ""); void SetParameterValue(paramId_t, double); double GetParameterValue(paramId_t); double GetParameterStep(paramId_t); double GetParameterStart(paramId_t); double GetParameterMin(paramId_t); double GetParameterMax(paramId_t); bool IsValidParameterValue(paramId_t, double); std::string StateString(); }; } // namespace params } // namespace nuis #endif diff --git a/src/persistency/ROOTOutput.cxx b/src/persistency/ROOTOutput.cxx index 7d83df5..8103ee8 100644 --- a/src/persistency/ROOTOutput.cxx +++ b/src/persistency/ROOTOutput.cxx @@ -1,45 +1,64 @@ #include "persistency/ROOTOutput.hxx" #include "utility/ROOTUtility.hxx" #include "config/GlobalConfiguration.hxx" #include "fhiclcpp/ParameterSet.h" #include "TFile.h" namespace nuis { namespace persistency { struct NamedTFile { std::string name; - std::shared_ptr file; + std::unique_ptr file; + + NamedTFile() : name(""), file(nullptr) {} + NamedTFile(NamedTFile &&other) + : name(std::move(other.name)), file(std::move(other.file)) {} + ~NamedTFile() { + if (file) { + file->Write(); + file->Close(); + } + } }; + static std::vector Files; -std::shared_ptr GetOutputFile(std::string const &name) { +std::unique_ptr &GetOutputFile(std::string const &name) { for (NamedTFile &file : Files) { if (file.name == name) { return file.file; } } fhicl::ParameterSet const &persistency = config::GetDocument().get("persistency"); std::string file_name = persistency.get(name + ".output_file"); std::string open_opts = persistency.get(name + ".open_mode", "CREATE"); NamedTFile ntf; ntf.name = name; - ntf.file = std::make_shared(file_name.c_str(), open_opts.c_str()); + ntf.file = std::make_unique(file_name.c_str(), open_opts.c_str()); if (!ntf.file || !ntf.file->IsOpen()) { throw utility::failed_to_open_TFile() << "[ERROR]: Failed to open output file: " << std::quoted(file_name) << " in write mode with opts = " << std::quoted(open_opts); } - Files.push_back(ntf); + Files.push_back(std::move(ntf)); return Files.back().file; } + +void CloseOpenTFiles() { + for (NamedTFile &f : Files) { + std::cout << "[INFO]: Closing open TFile: " << f.name << " " + << f.file->GetName() << std::endl; + } + Files.clear(); +} } // namespace persistency } // namespace nuis diff --git a/src/persistency/ROOTOutput.hxx b/src/persistency/ROOTOutput.hxx index 3b45a71..1aa5bb6 100644 --- a/src/persistency/ROOTOutput.hxx +++ b/src/persistency/ROOTOutput.hxx @@ -1,66 +1,76 @@ #ifndef PERSITENCY_ROOTOUTPUT_HXX_SEEN #define PERSITENCY_ROOTOUTPUT_HXX_SEEN #include "exception/exception.hxx" #include "TFile.h" #include #include #include namespace nuis { namespace persistency { NEW_NUIS_EXCEPT(WriteToOutputFile_nullptr); /// Will get/open a TFile that is described in the global config /// /// The named streams will be used to configure the file name and open mode from /// the global config element persistency.: {file: output.root opts: /// CREATE} -std::shared_ptr GetOutputFile(std::string const &name = "default"); +std::unique_ptr &GetOutputFile(std::string const &name = "default"); template inline void WriteToOutputFile(T *object, std::string const &object_name, - std::string dir_name = "", - std::string const &file_name = "default") { + std::string dir_name = "", + std::string const &file_name = "default") { -if(!object){ - throw WriteToOutputFile_nullptr(); -} + if (!object) { + throw WriteToOutputFile_nullptr(); + } TDirectory *ogdir = gDirectory; - std::shared_ptr f = GetOutputFile(file_name); + std::unique_ptr &f = GetOutputFile(file_name); TDirectory *d = f.get(); while (dir_name.length()) { size_t next_slash = dir_name.find_first_of('/'); std::string next_dir = dir_name.substr(0, next_slash); if (next_slash != std::string::npos) { dir_name = dir_name.substr(next_slash + 1); } else { dir_name = ""; } TDirectory *nd = d->GetDirectory(next_dir.c_str()); if (!nd) { nd = d->mkdir(next_dir.c_str()); } nd->cd(); d = nd; } d->cd(); object->Write(object_name.c_str(), TObject::kOverwrite); if (ogdir) { ogdir->cd(); } } + +template +inline void WriteToOutputFile(std::unique_ptr &object, + std::string const &object_name, + std::string dir_name = "", + std::string const &file_name = "default") { + return WriteToOutputFile(object.get(), object_name, dir_name, file_name); +} + +void CloseOpenTFiles(); } // namespace persistency } // namespace nuis #endif diff --git a/src/plugins/CMakeLists.txt b/src/plugins/CMakeLists.txt index 5c98f59..b2ab67f 100644 --- a/src/plugins/CMakeLists.txt +++ b/src/plugins/CMakeLists.txt @@ -1,5 +1,13 @@ +SET(plugins_implementation_files + LoadedPlugins.cxx) + set(plugins_header_files -Instantiate.hxx -traits.hxx) + Instantiate.hxx + NamedSO.hxx + traits.hxx) + +add_library(nuis_plugins SHARED ${plugins_implementation_files}) +target_link_libraries(nuis_plugins) +install(TARGETS nuis_plugins DESTINATION lib) install(FILES ${plugins_header_files} DESTINATION include/plugins) diff --git a/src/plugins/Instantiate.hxx b/src/plugins/Instantiate.hxx index 05c2bbd..8e7bf57 100644 --- a/src/plugins/Instantiate.hxx +++ b/src/plugins/Instantiate.hxx @@ -1,273 +1,254 @@ // Copyright 2018 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 PLUGINS_PLUGINMANAGER_HXX_SEEN #define PLUGINS_PLUGINMANAGER_HXX_SEEN +#include "plugins/NamedSO.hxx" #include "plugins/traits.hxx" #include "config/GlobalConfiguration.hxx" #include "utility/FileSystemUtility.hxx" #include "exception/exception.hxx" #include "fhiclcpp/ParameterSet.h" #include "string_parsers/to_string.hxx" // linux #include #include #include #include #include #include #include // #define DEBUG_INSTANTIATE namespace nuis { namespace plugins { +extern std::vector LoadedSharedObjects; + NEW_NUIS_EXCEPT(failed_to_find_instantiator); NEW_NUIS_EXCEPT(malformed_plugin_interface); NEW_NUIS_EXCEPT(failed_to_load_so); typedef void *(*inst_fcn)(); typedef void (*dltr_fcn)(void *); template struct PluginInstantiator { std::string FQ_so_path; std::string Base_classname; std::string Classname; void *dllib; inst_fcn Instantiator; dltr_fcn Deleter; PluginInstantiator() : FQ_so_path(""), Base_classname(""), Classname(""), dllib(nullptr), Instantiator(nullptr), Deleter(nullptr) {} PluginInstantiator(PluginInstantiator const &) = delete; PluginInstantiator(PluginInstantiator &&other) { FQ_so_path = std::move(other.FQ_so_path); Base_classname = std::move(other.Base_classname); Classname = std::move(other.Classname); dllib = other.dllib; Instantiator = other.Instantiator; Deleter = other.Deleter; other.FQ_so_path = ""; other.Base_classname = ""; other.Classname = ""; other.dllib = nullptr; other.Instantiator = nullptr; other.Deleter = nullptr; } typename plugin_traits::unique_ptr_t Instantiate() { T *inst = reinterpret_cast((*Instantiator)()); dltr_fcn dltr = Deleter; std::string cln = Classname; std::function deleter = [=](T *inst) { #ifdef DEBUG_INSTANTIATE std::cout << "[INFO]: Deleting instance of " << cln << " with " << (void *)dltr << std::endl; #endif (*dltr)(inst); }; return typename plugin_traits::unique_ptr_t(inst, deleter); } }; -struct NamedSO { - std::string name; - void *dllib; - - NamedSO() : name(""), dllib(nullptr) {} - NamedSO(NamedSO const &) = delete; - NamedSO(NamedSO &&other) : name(std::move(other.name)), dllib(other.dllib) { - other.dllib = nullptr; - } - - ~NamedSO() { - if (dllib) { -#ifdef DEBUG_INSTANTIATE - std::cout << "[INFO]: dlclose on shared object: " << std::quoted(name) - << std::endl; -#endif - dlclose(dllib); - } - } -}; - NamedSO &GetSharedObject(std::string const &FQPath) { - static std::vector LoadedSharedObjects; for (NamedSO &so : LoadedSharedObjects) { if (so.name == FQPath) { return so; } } NamedSO so; so.name = FQPath; so.dllib = dlopen(FQPath.c_str(), RTLD_NOW | RTLD_GLOBAL); char const *dlerr_cstr = dlerror(); std::string dlerr; if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { throw failed_to_load_so() << "[INFO]: Failed to load shared object: " << FQPath << " with dlerror: " << dlerr; } else { #ifdef DEBUG_INSTANTIATE std::cout << "[INFO]: Loaded shared object " << FQPath << std::endl; #endif } LoadedSharedObjects.push_back(std::move(so)); return LoadedSharedObjects.back(); } template typename plugin_traits::unique_ptr_t Instantiate(std::string const &classname) { static std::vector> LoadedPlugins; fhicl::ParameterSet const &plugins = config::GetDocument().get("plugins"); fhicl::ParameterSet const &search_paths = plugins.get("search_paths"); std::vector plugin_search_dirs; // Look for plugin search paths in sequence elements of the // plugins.search_paths table for (std::string const &key : search_paths.get_names()) { if (!search_paths.is_key_to_sequence(key)) { continue; } for (std::string const &path : search_paths.get>(key)) { plugin_search_dirs.push_back(path); } } for (std::string path : plugin_search_dirs) { path = utility::EnsureTrailingSlash(path); for (std::string const &so_name : utility::GetMatchingFiles(path, ".*\\.so")) { for (PluginInstantiator &plugin : LoadedPlugins) { if (plugin.FQ_so_path == (path + so_name) && (plugin.Base_classname == plugin_traits::interface_name()) && (plugin.Classname == classname)) { #ifdef DEBUG_INSTANTIATE std::cout << "[INFO]: Using already loaded PluginInstantiator" << std::endl; #endif return plugin.Instantiate(); } } PluginInstantiator plugin; plugin.FQ_so_path = path + so_name; plugin.Base_classname = plugin_traits::interface_name(); plugin.Classname = classname; plugin.dllib = GetSharedObject(plugin.FQ_so_path).dllib; char const *dlerr_cstr = nullptr; std::string dlerr(""); plugin.Instantiator = reinterpret_cast(dlsym( plugin.dllib, plugin_traits::instantiator_function_name(classname).c_str())); dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr_cstr) { #ifdef DEBUG_INSTANTIATE std::cout << "[INFO]: Failed to load appropriate instantiator method: " << plugin_traits::instantiator_function_name(classname) << " from shared object " << plugin.FQ_so_path; #endif continue; } else { #ifdef DEBUG_INSTANTIATE std::cout << "[INFO]: Loaded instantiator method: " << plugin_traits::instantiator_function_name(classname) << " from shared object " << plugin.FQ_so_path << std::endl; #endif } plugin.Deleter = reinterpret_cast( dlsym(plugin.dllib, plugin_traits::deleter_function_name(classname).c_str())); dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr_cstr) { throw malformed_plugin_interface() << "[ERROR]: Failed to load appropriate deleter method: " << plugin_traits::deleter_function_name(classname) << " from shared object " << plugin.FQ_so_path << " with error: " << std::quoted(dlerr); } else { #ifdef DEBUG_INSTANTIATE std::cout << "[INFO]: Loaded deleter method: " << plugin_traits::deleter_function_name(classname) << " from shared object " << plugin.FQ_so_path << std::endl; #endif } #ifdef DEBUG_INSTANTIATE std::cout << "[INFO]: Checking if shared object " << std::quoted(plugin.FQ_so_path) << " knows how to instantiate class " << std::quoted(classname) << " via interface " << std::quoted(plugin_traits::interface_name()) << std::endl; #endif LoadedPlugins.push_back(std::move(plugin)); return LoadedPlugins.back().Instantiate(); } } throw failed_to_find_instantiator() << "[ERROR]: Failed to find instantiator for classname: " << std::quoted(classname) << " using interface " << std::quoted(plugin_traits::interface_name()) << " from configured search paths: " << fhicl::string_parsers::T2Str>( plugin_search_dirs); } } // namespace plugins } // namespace nuis #endif diff --git a/src/plugins/LoadedPlugins.cxx b/src/plugins/LoadedPlugins.cxx new file mode 100644 index 0000000..86e1edf --- /dev/null +++ b/src/plugins/LoadedPlugins.cxx @@ -0,0 +1,9 @@ +#include "plugins/NamedSO.hxx" + +#include + +namespace nuis { +namespace plugins { +std::vector LoadedSharedObjects; +} +} // namespace nuis diff --git a/src/plugins/NamedSO.hxx b/src/plugins/NamedSO.hxx new file mode 100644 index 0000000..5137700 --- /dev/null +++ b/src/plugins/NamedSO.hxx @@ -0,0 +1,32 @@ +#ifndef PLUGINS_NAMEDSO_HXX_SEEN +#define PLUGINS_NAMEDSO_HXX_SEEN + +// linux +#include + +#include +#include +#include + +struct NamedSO { + std::string name; + void *dllib; + + NamedSO() : name(""), dllib(nullptr) {} + NamedSO(NamedSO const &) = delete; + NamedSO(NamedSO &&other) : name(std::move(other.name)), dllib(other.dllib) { + other.dllib = nullptr; + } + + ~NamedSO() { + if (dllib) { +#ifdef DEBUG_INSTANTIATE + std::cout << "[INFO]: dlclose on shared object: " << std::quoted(name) + << std::endl; +#endif + dlclose(dllib); + } + } +}; + +#endif diff --git a/src/samples/CMakeLists.txt b/src/samples/CMakeLists.txt index ed8bfaf..a31710c 100644 --- a/src/samples/CMakeLists.txt +++ b/src/samples/CMakeLists.txt @@ -1,20 +1,23 @@ set(samples_header_files - ISample.hxx + IEventProcessor.hxx IDataComparison.hxx SimpleDataComparison.hxx) install(FILES ${samples_header_files} DESTINATION include/samples) add_subdirectory(MCTools) add_subdirectory(nuA) cmessage(DEBUG "INuADataComparisons: ${INuADataComparisons}") SET(INuADataComparisons_List) if(NOT IDataComparisons STREQUAL "") - string(REPLACE ";" "," INuADataComparisons_List ${INuADataComparisons}) + string(REPLACE ";" ", " INuADataComparisons_List "${INuADataComparisons}") endif() +cmessage(DEBUG "INuADataComparisons_List: ${INuADataComparisons_List}") + + SET(INuADataComparisons_List ${INuADataComparisons_List} PARENT_SCOPE) SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE) diff --git a/src/samples/IDataComparison.hxx b/src/samples/IDataComparison.hxx index cfda1fc..1fdc67d 100644 --- a/src/samples/IDataComparison.hxx +++ b/src/samples/IDataComparison.hxx @@ -1,73 +1,73 @@ // Copyright 2018 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 SAMPLES_IDATACOMPARISON_HXX_SEEN #define SAMPLES_IDATACOMPARISON_HXX_SEEN -#include "samples/ISample.hxx" +#include "samples/IEventProcessor.hxx" #include "fhiclcpp/ParameterSet.h" #include #include -class IDataComparison : public ISample { +class IDataComparison : public IEventProcessor { public: virtual double GetGOF() = 0; virtual double GetNDOGuess() = 0; virtual std::string GetJournalReference() { std::stringstream ss(""); ss << "Unknown Journal Ref. for IDataComparison: " << std::quoted(Name()); return ss.str(); } virtual std::string GetTargetMaterial() { std::stringstream ss(""); ss << "Unknown Target material for IDataComparison: " << std::quoted(Name()); return ss.str(); } virtual std::string GetFluxDescription() { std::stringstream ss(""); ss << "Unknown Flux description for IDataComparison: " << std::quoted(Name()); return ss.str(); } virtual std::string GetSignalDescription() { std::stringstream ss(""); ss << "Unknown Signal description for IDataComparison: " << std::quoted( Name()); return ss.str(); } virtual std::string GetDocumentation() { std::stringstream ss(""); ss << "No documentation provided for IDataComparison: " << std::quoted(Name()); return ss.str(); } virtual fhicl::ParameterSet GetExampleConfiguration() { return fhicl::ParameterSet(); } }; DECLARE_PLUGIN_INTERFACE(IDataComparison); #endif diff --git a/src/samples/ISample.hxx b/src/samples/IEventProcessor.hxx similarity index 74% rename from src/samples/ISample.hxx rename to src/samples/IEventProcessor.hxx index 0b195b6..90eff1a 100644 --- a/src/samples/ISample.hxx +++ b/src/samples/IEventProcessor.hxx @@ -1,113 +1,119 @@ // Copyright 2018 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 SAMPLES_ISAMPLE_HXX_SEEN -#define SAMPLES_ISAMPLE_HXX_SEEN +#ifndef SAMPLES_IEventProcessor_HXX_SEEN +#define SAMPLES_IEventProcessor_HXX_SEEN #include "plugins/traits.hxx" #include "exception/exception.hxx" #include "config/GlobalConfiguration.hxx" #include "fhiclcpp/ParameterSet.h" +#include "TH1.h" + #include #include #include namespace nuis { namespace event { class FullEvent; class MinimalEvent; } // namespace event } // namespace nuis -#define ISAMPLE_DEBUG(v) \ +#define IEventProcessor_DEBUG(v) \ if (verb > 2) { \ std::cout << "[DEBUG]: " << v << std::endl; \ } -#define ISAMPLE_INFO(v) \ +#define IEventProcessor_INFO(v) \ if (verb > 1) { \ std::cout << "[INFO]: " << v << std::endl; \ } -#define ISAMPLE_WARN(v) \ +#define IEventProcessor_WARN(v) \ if (verb > 0) { \ std::cout << "[WARN] " << __FILENAME__ << ":" << __LINE__ << " : " << v \ << std::endl; \ } -class ISample { +class IEventProcessor { protected: - NEW_NUIS_EXCEPT(unknown_ISample_verbosity); + NEW_NUIS_EXCEPT(unknown_IEventProcessor_verbosity); enum sample_verbosity { kSilent = 0, kWarn = 1, kReticent = 2, kVerbose = 3 }; sample_verbosity verb; void SetSampleVerbosity(std::string v) { v = nuis::config::GetDocument().get( "global.sample.verbosity_override", v); if (v == "Silent") { verb = kSilent; } else if (v == "Warn") { verb = kWarn; } else if (v == "Reticent") { verb = kReticent; } else if (v == "Verbose") { verb = kVerbose; } else { - throw unknown_ISample_verbosity() - << "[ERROR]: Failed to parse ISample verbosity setting from: " + throw unknown_IEventProcessor_verbosity() + << "[ERROR]: Failed to parse IEventProcessor verbosity setting from: " << std::quoted(v); } } public: - NEW_NUIS_EXCEPT(uninitialized_ISample); - NEW_NUIS_EXCEPT(unimplemented_ISample_optional_method); + NEW_NUIS_EXCEPT(uninitialized_IEventProcessor); + NEW_NUIS_EXCEPT(unimplemented_IEventProcessor_optional_method); - ISample() { + IEventProcessor() { SetSampleVerbosity(nuis::config::GetDocument().get( "global.sample.verbosity_default", "Silent")); + + TH1::SetDefaultSumw2(); } virtual void Initialize(fhicl::ParameterSet const &) = 0; // Interface for processing a single, external event // - // ISamples are not required to implement processing events from 'outside'. + // IEventProcessors are not required to implement processing events from + // 'outside'. virtual void ProcessEvent(nuis::event::FullEvent const &) { - throw unimplemented_ISample_optional_method() - << "[ERROR]: ISample " << Name() << " does not implement ProcessEvent."; + throw unimplemented_IEventProcessor_optional_method() + << "[ERROR]: IEventProcessor " << Name() + << " does not implement ProcessEvent."; } virtual void ProcessSample(size_t nmax = std::numeric_limits::max()) = 0; virtual void Write() = 0; virtual std::string Name() = 0; - virtual ~ISample() {} + virtual ~IEventProcessor() {} }; -DECLARE_PLUGIN_INTERFACE(ISample); +DECLARE_PLUGIN_INTERFACE(IEventProcessor); #endif diff --git a/src/samples/MCTools/NuisToNuWro.cxx b/src/samples/MCTools/NuisToNuWro.cxx index b408c53..dd1e6c4 100644 --- a/src/samples/MCTools/NuisToNuWro.cxx +++ b/src/samples/MCTools/NuisToNuWro.cxx @@ -1,103 +1,103 @@ -#include "samples/ISample.hxx" +#include "samples/IEventProcessor.hxx" #include "event/FullEvent.hxx" #include "input/InputManager.hxx" #include "utility/FullEventUtility.hxx" #include "utility/ROOTUtility.hxx" #include "fhiclcpp/ParameterSet.h" #include "generator/utility/NuWroUtility.hxx" #include #include using namespace nuis::event; using namespace nuis::input; using namespace nuis::utility; using namespace nuis::nuwrotools; -class NuisToNuWro : public ISample { +class NuisToNuWro : public IEventProcessor { public: InputManager::Input_id_t fIH_id; - std::unique_ptr fOutputTree; + TreeFile fOutputTree; NuWroEvent *fOutputEvent; NuisToNuWro() : fIH_id(std::numeric_limits::max()), - fOutputTree(nullptr), fOutputEvent(nullptr) {} + fOutputTree(), fOutputEvent(nullptr) {} void Initialize(fhicl::ParameterSet const &ps) { fIH_id = InputManager::Get().EnsureInputLoaded(ps); fOutputTree = MakeNewTTree(ps.get("output_file"), "treeout", "RECREATE"); - fOutputTree->tree->Branch("e", &fOutputEvent); + fOutputTree.tree->Branch("e", &fOutputEvent); } void ProcessEvent(FullEvent const &ps) { fOutputEvent->in.clear(); fOutputEvent->out.clear(); fOutputEvent->post.clear(); std::pair NuMode = GetFlagsDynEquivalent(ps.mode); fOutputEvent->flag = NuMode.first; fOutputEvent->dyn = NuMode.second; for (Particle const &part : GetISParticles(ps)) { particle nuwro_part(part.pdg, part.P4.M()); nuwro_part.set_momentum(vect(part.P4.X(), part.P4.Y(), part.P4.Z())); fOutputEvent->in.push_back(nuwro_part); } for (Particle const &part : GetPrimaryFSParticles(ps)) { particle nuwro_part(part.pdg, part.P4.M()); nuwro_part.set_momentum(vect(part.P4.X(), part.P4.Y(), part.P4.Z())); fOutputEvent->out.push_back(nuwro_part); } for (Particle const &part : GetNuclearLeavingParticles(ps)) { particle nuwro_part(part.pdg, part.P4.M()); nuwro_part.set_momentum(vect(part.P4.X(), part.P4.Y(), part.P4.Z())); fOutputEvent->post.push_back(nuwro_part); } - fOutputTree->tree->Fill(); + fOutputTree.tree->Fill(); } void ProcessSample(size_t nmax) { if (fIH_id == std::numeric_limits::max()) { - throw uninitialized_ISample(); + throw uninitialized_IEventProcessor(); } IInputHandler const &IH = InputManager::Get().GetInputHandler(fIH_id); size_t NEvsToProcess = std::min(nmax, IH.GetNEvents()); size_t NToShout = NEvsToProcess / 10; std::cout << "[INFO]: Processing " << NEvsToProcess << " input events to NuWro format." << std::endl; size_t n = 0; for (auto const &fe : IH) { if (++n > NEvsToProcess) { break; } if (NToShout && !(n % NToShout)) { std::cout << "[INFO]: Processed " << n << "/" << NEvsToProcess << " events." << std::endl; } ProcessEvent(fe); } } - void Write() { fOutputTree->file->Write(); } + void Write() { fOutputTree.file->Write(); } std::string Name() { return "NuisToNuWro"; } }; -DECLARE_PLUGIN(ISample, NuisToNuWro); +DECLARE_PLUGIN(IEventProcessor, NuisToNuWro); diff --git a/src/samples/MCTools/VerboseEventSummary.cxx b/src/samples/MCTools/VerboseEventSummary.cxx index 726a5a5..a007433 100644 --- a/src/samples/MCTools/VerboseEventSummary.cxx +++ b/src/samples/MCTools/VerboseEventSummary.cxx @@ -1,56 +1,56 @@ -#include "samples/ISample.hxx" +#include "samples/IEventProcessor.hxx" #include "event/FullEvent.hxx" #include "input/InputManager.hxx" #include #include using namespace nuis::event; using namespace nuis::input; -class VerboseEventSummary : public ISample { +class VerboseEventSummary : public IEventProcessor { public: InputManager::Input_id_t fIH_id; VerboseEventSummary() : fIH_id(std::numeric_limits::max()) {} void Initialize(fhicl::ParameterSet const &ps) { fIH_id = InputManager::Get().EnsureInputLoaded(ps); } void ProcessEvent(FullEvent const &ps) { std::cout << ps.to_string() << std::endl; } void ProcessSample(size_t nmax) { if (fIH_id == std::numeric_limits::max()) { - throw uninitialized_ISample(); + throw uninitialized_IEventProcessor(); } IInputHandler const &IH = InputManager::Get().GetInputHandler(fIH_id); size_t NEvsToProcess = std::min(nmax, IH.GetNEvents()); size_t NToShout = NEvsToProcess / 10; std::cout << "[INFO]: Processing " << NEvsToProcess << " input events." << std::endl; size_t n = 0; for (auto const &fe : IH) { if (++n > nmax) { break; } if (NToShout && !(n % NToShout)) { std::cout << "[INFO]: Processed " << n << "/" << NEvsToProcess << " events." << std::endl; } ProcessEvent(fe); } } void Write() {} std::string Name() { return "VerboseEventSummary"; } }; -DECLARE_PLUGIN(ISample, VerboseEventSummary); +DECLARE_PLUGIN(IEventProcessor, VerboseEventSummary); diff --git a/src/samples/SimpleDataComparison.hxx b/src/samples/SimpleDataComparison.hxx index bdce072..36f8e7a 100644 --- a/src/samples/SimpleDataComparison.hxx +++ b/src/samples/SimpleDataComparison.hxx @@ -1,375 +1,424 @@ // Copyright 2018 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 SAMPLES_SIMPLEDATACOMPARISON_HXX_SEEN #define SAMPLES_SIMPLEDATACOMPARISON_HXX_SEEN #include "samples/IDataComparison.hxx" #include "event/FullEvent.hxx" #include "input/InputManager.hxx" #include "persistency/ROOTOutput.hxx" #include "utility/FileSystemUtility.hxx" #include "utility/HistogramUtility.hxx" +#include "utility/KinematicUtility.hxx" #include "utility/ROOTUtility.hxx" #include "utility/StatsUtility.hxx" #include #include #include #include #include -template struct TH_dim_helper {}; -template <> struct TH_dim_helper<1> { typedef TH1D type; }; -template <> struct TH_dim_helper<2> { typedef TH2D type; }; - -template class SimpleDataComparison : public IDataComparison { +template ::type> +class SimpleDataComparison : public IDataComparison { NEW_NUIS_EXCEPT(invalid_SimpleDataComparison_initialization); NEW_NUIS_EXCEPT(SimpleDataComparison_already_finalized); protected: + using HistType = HT; + using TH_Help = typename nuis::utility::TH_Helper; + nuis::input::InputManager::Input_id_t fIH_id; std::string write_directory; size_t NMaxSample_override; int fIsShapeOnly; int fIsFluxUnfolded; std::vector fSignalCache; - std::vector> fProjectionCache; + std::vector> fProjectionCache; std::string fDataInputDescriptor; - std::unique_ptr::type> fData; + std::unique_ptr fData; std::string fMaskInputDescriptor; - std::unique_ptr::type> fMask; + std::unique_ptr fMask; std::string fCovarianceInputDescriptor; - std::unique_ptr fCovariance; - std::unique_ptr::type> fPrediction; - std::unique_ptr::type> fPrediction_xsec; - std::unique_ptr::type> fPrediction_shape; - std::unique_ptr::type> fPrediction_comparison; + std::unique_ptr fCovariance; + std::unique_ptr fPrediction; + std::unique_ptr fPrediction_xsec; + std::unique_ptr fPrediction_shape; + std::unique_ptr fPrediction_comparison; bool fComparisonFinalized; std::string fJournalReference; std::string fTargetMaterial; std::string fFluxDescription; std::string fSignalDescription; - std::pair EnuRange; + nuis::utility::ENuRange energy_cut; std::function IsSigFunc; - std::function(nuis::event::FullEvent const &)> + std::function(nuis::event::FullEvent const &)> CompProjFunc; + // If assigned by subclass will be called on for all events, bool signifies + // whether the event was selected. + std::function + ProcessExtraFunc; + public: SimpleDataComparison() { fIH_id = std::numeric_limits::max(); write_directory = ""; NMaxSample_override = std::numeric_limits::max(); fDataInputDescriptor = ""; fData = nullptr; fMaskInputDescriptor = ""; fMask = nullptr; fCovarianceInputDescriptor = ""; fCovariance = nullptr; fPrediction = nullptr; fPrediction_xsec = nullptr; fPrediction_shape = nullptr; fPrediction_comparison = nullptr; fComparisonFinalized = false; IsSigFunc = [](nuis::event::FullEvent const &) -> bool { return true; }; CompProjFunc = - [](nuis::event::FullEvent const &) -> std::array { - std::array arr; - for (double &el : arr) { + [](nuis::event::FullEvent const &) -> std::array { + std::array arr; + for (NumericT &el : arr) { el = 0; } return arr; }; + ProcessExtraFunc = + std::function(); + fJournalReference = ""; fTargetMaterial = ""; fFluxDescription = ""; fSignalDescription = ""; fIsShapeOnly = -1; fIsFluxUnfolded = -1; - EnuRange = std::pair{std::numeric_limits::max(), + + energy_cut = nuis::utility::ENuRange{std::numeric_limits::max(), std::numeric_limits::max()}; } + fhicl::ParameterSet fGlobalConfig; + fhicl::ParameterSet fInstanceConfig; + void ReadGlobalConfigDefaults() { - fhicl::ParameterSet const &global_sample_configuration = - nuis::config::GetDocument().get( - std::string("global.sample_configuration.") + Name(), - fhicl::ParameterSet()); + fGlobalConfig = nuis::config::GetDocument().get( + std::string("global.sample_configuration.") + Name(), + fhicl::ParameterSet()); if (!fJournalReference.length()) { - fJournalReference = global_sample_configuration.get( - "journal_reference", fJournalReference); + fJournalReference = fGlobalConfig.get("journal_reference", + fJournalReference); } if (!fTargetMaterial.length()) { - fTargetMaterial = global_sample_configuration.get( - "target_material", fTargetMaterial); + fTargetMaterial = + fGlobalConfig.get("target_material", fTargetMaterial); } if (!fFluxDescription.length()) { - fFluxDescription = global_sample_configuration.get( - "flux_description", fFluxDescription); + fFluxDescription = + fGlobalConfig.get("flux_description", fFluxDescription); } if (!fSignalDescription.length()) { - fSignalDescription = global_sample_configuration.get( - "signal_description", fSignalDescription); + fSignalDescription = fGlobalConfig.get("signal_description", + fSignalDescription); } if (fIsShapeOnly == -1) { - fIsShapeOnly = global_sample_configuration.get("shape_only", false); + fIsShapeOnly = fGlobalConfig.get("shape_only", false); } if (fIsFluxUnfolded == -1) { - fIsFluxUnfolded = - global_sample_configuration.get("flux_unfolded", false); + fIsFluxUnfolded = fGlobalConfig.get("flux_unfolded", false); } - if ((EnuRange.first == std::numeric_limits::max()) && - (global_sample_configuration.has_key("enu_range"))) { - EnuRange = global_sample_configuration.get>( - "enu_range"); + if ((energy_cut.first == std::numeric_limits::max()) && + (fGlobalConfig.has_key("enu_range"))) { + energy_cut = fGlobalConfig.get>("enu_range"); } } virtual std::string GetJournalReference() { return fJournalReference; } virtual std::string GetTargetMaterial() { return fTargetMaterial; } virtual std::string GetFluxDescription() { return fFluxDescription; } virtual std::string GetSignalDescription() { return fSignalDescription; } void SetShapeOnly(bool iso) { fIsShapeOnly = iso; } void SetFluxUnfolded(bool ifo) { fIsFluxUnfolded = ifo; } void SetData(std::string const &data_descriptor) { fDataInputDescriptor = data_descriptor; } void SetMask(std::string const &mask_descriptor) { fMaskInputDescriptor = mask_descriptor; } void SetCovariance(std::string const &cov_descriptor) { fCovarianceInputDescriptor = cov_descriptor; } - virtual void FillProjection(std::array const &proj, - double event_weight) { - nuis::utility::Fill(fPrediction.get(), proj, event_weight); + virtual void FillProjection(std::array const &proj, + NumericT event_weight) { + TH_Help::Fill(fPrediction, proj, event_weight); } virtual void FinalizeComparison() { if (fComparisonFinalized) { throw SimpleDataComparison_already_finalized() << "[ERROR]: Attempted to re-finalize a comparison for " - "IDataComparison: " + "SimpleDataComparison: " << std::quoted(Name()); } - fPrediction_xsec = nuis::utility::Clone(fPrediction); - fPrediction_xsec->Scale(1.0, "width"); - fPrediction_shape = nuis::utility::Clone(fPrediction_xsec); + fPrediction_xsec = + nuis::utility::Clone(fPrediction, false, "Prediction_xsec"); + + IInputHandler const &IH = + nuis::input::InputManager::Get().GetInputHandler(fIH_id); + + TH_Help::Scale(fPrediction_xsec, 1.0, "width"); + + // If we have a flux cut + if (energy_cut.first != std::numeric_limits::max()) { + TH_Help::Scale(fPrediction_xsec, IH.GetXSecScaleFactor(energy_cut)); + } + + fPrediction_shape = + nuis::utility::Clone(fPrediction_xsec, false, "Prediction_shape"); if (fData) { - fPrediction_shape->Scale(fData->Integral() / - fPrediction_shape->Integral()); + TH_Help::Scale(fPrediction_shape, + TH_Help::Integral(fData, "width") / + TH_Help::Integral(fPrediction_shape, "width")); + } else { - ISAMPLE_WARN("When Finalizing comparison, no Data histogram available."); + IEventProcessor_WARN( + "When Finalizing comparison, no Data histogram available."); } if (fIsFluxUnfolded) { // fPrediction_comparison } else if (fIsShapeOnly) { - fPrediction_comparison = nuis::utility::Clone(fPrediction_shape); + fPrediction_comparison = nuis::utility::Clone(fPrediction_shape, false, + "Prediction_comparison"); } else { - fPrediction_comparison = nuis::utility::Clone(fPrediction_xsec); + fPrediction_comparison = nuis::utility::Clone(fPrediction_xsec, false, + "Prediction_comparison"); } fComparisonFinalized = true; } - void Initialize(fhicl::ParameterSet const &ps) { + void Initialize(fhicl::ParameterSet const &instance_sample_configuration) { - if (ps.has_key("verbosity")) { - SetSampleVerbosity(ps.get("verbosity")); + fInstanceConfig = instance_sample_configuration; + + if (fInstanceConfig.has_key("verbosity")) { + SetSampleVerbosity(fInstanceConfig.get("verbosity")); } else { SetSampleVerbosity("Reticent"); } ReadGlobalConfigDefaults(); - fhicl::ParameterSet const &global_sample_configuration = - nuis::config::GetDocument().get( - std::string("global.sample_configuration.") + Name(), - fhicl::ParameterSet()); - - if (ps.has_key("fake_data")) { - fData = nuis::utility::GetHistogram::type>( - ps.get("fake_data_histogram")); + if (fInstanceConfig.has_key("fake_data")) { + fData = nuis::utility::GetHistogram( + fInstanceConfig.get("fake_data_histogram")); + } else if (!fGlobalConfig.get("has_data", true) || + !fInstanceConfig.get("has_data", true)) { + // Explicitly not expecting data } else { if (!fDataInputDescriptor.length()) { - if (!global_sample_configuration.has_key("data_descriptor")) { + if (!fGlobalConfig.has_key("data_descriptor")) { throw invalid_SimpleDataComparison_initialization() << "[ERROR]: SimpleDataComparison::Initialize for " "IDataComparison: " << std::quoted(Name()) << " failed as no input data was set by a call to " "SimpleDataComparison::SetData and no data_descriptor for " "this SimpleDataComparison could be found in the global " "configuration."; } fDataInputDescriptor = - global_sample_configuration.get("data_descriptor"); + fGlobalConfig.get("data_descriptor"); } - fData = nuis::utility::GetHistogram::type>( - fDataInputDescriptor); + fData = nuis::utility::GetHistogram(fDataInputDescriptor); } - fPrediction = nuis::utility::Clone(fData, true); + if (!fPrediction) { + if (!fData) { + throw invalid_SimpleDataComparison_initialization() + << "[ERROR]: SimpleDataComparison::Initialize for " + "IDataComparison: " + << std::quoted(Name()) + << " failed. As `has_data: false` was set in the configuration " + "(global: " + << (!fGlobalConfig.get("has_data", true)) + << ", instance: " << !fInstanceConfig.get("has_data", true) + << "), the instance constructor must supply the fPrediction " + "binning, and it wasn't."; + } + fPrediction = nuis::utility::Clone(fData, true, "Prediction"); + } if (fCovarianceInputDescriptor.length()) { fCovariance = nuis::utility::GetHistogram(fCovarianceInputDescriptor); - } else if (global_sample_configuration.has_key("covariance_descriptor")) { + } else if (fGlobalConfig.has_key("covariance_descriptor")) { fCovariance = nuis::utility::GetHistogram( - global_sample_configuration.get( - "covariance_descriptor")); + fGlobalConfig.get("covariance_descriptor")); } if (fMaskInputDescriptor.length()) { - fMask = nuis::utility::GetHistogram::type>( - fMaskInputDescriptor); - } else if (global_sample_configuration.has_key("mask_descriptor")) { - fMask = nuis::utility::GetHistogram::type>( - global_sample_configuration.get("mask_descriptor")); + fMask = nuis::utility::GetHistogram(fMaskInputDescriptor); + } else if (fGlobalConfig.has_key("mask_descriptor")) { + fMask = nuis::utility::GetHistogram( + fGlobalConfig.get("mask_descriptor")); } - if (ps.has_key("verbosity")) { - SetSampleVerbosity(ps.get("verbosity")); + if (fInstanceConfig.has_key("verbosity")) { + SetSampleVerbosity(fInstanceConfig.get("verbosity")); } NMaxSample_override = - ps.get("nmax", std::numeric_limits::max()); + fInstanceConfig.get("nmax", std::numeric_limits::max()); - write_directory = ps.get("write_directory", ""); + write_directory = + fInstanceConfig.get("write_directory", Name()); - fIH_id = nuis::input::InputManager::Get().EnsureInputLoaded(ps); + fIH_id = + nuis::input::InputManager::Get().EnsureInputLoaded(fInstanceConfig); } + void ProcessSample(size_t nmax = std::numeric_limits::max()) { if (fIH_id == std::numeric_limits::max()) { - throw uninitialized_ISample(); + throw uninitialized_IEventProcessor(); } IInputHandler const &IH = nuis::input::InputManager::Get().GetInputHandler(fIH_id); nmax = std::min(NMaxSample_override, nmax); size_t NEvsToProcess = std::min(nmax, IH.GetNEvents()); double nmax_scaling = double(IH.GetNEvents()) / double(NEvsToProcess); size_t NToShout = NEvsToProcess / 10; - ISAMPLE_INFO("Sample " << std::quoted(Name()) << " processing " - << NEvsToProcess << " events."); + IEventProcessor_INFO("Sample " << std::quoted(Name()) << " processing " + << NEvsToProcess << " events."); IInputHandler::ev_index_t ev_idx = 0; - size_t NSigEvents = 0; bool DetermineSignalEvents = !fSignalCache.size(); - nuis::utility::Clear(fPrediction.get()); + nuis::utility::Clear(*fPrediction); fComparisonFinalized = false; + size_t cache_ctr = 0; + while (ev_idx < NEvsToProcess) { if (DetermineSignalEvents) { nuis::event::FullEvent const &fev = IH.GetFullEvent(ev_idx); bool is_sig = IsSigFunc(fev); fSignalCache.push_back(is_sig); if (is_sig) { fProjectionCache.push_back(CompProjFunc(fev)); } + if (ProcessExtraFunc) { + ProcessExtraFunc(fev, is_sig, + IH.GetEventWeight(ev_idx) * nmax_scaling); + } } if (NToShout && !(ev_idx % NToShout)) { - ISAMPLE_INFO("\tProcessed " << ev_idx << "/" << NEvsToProcess - << " events."); + IEventProcessor_INFO("\tProcessed " << ev_idx << "/" << NEvsToProcess + << " events."); } if (fSignalCache[ev_idx]) { - FillProjection(fProjectionCache[ev_idx], + FillProjection(fProjectionCache[cache_ctr++], IH.GetEventWeight(ev_idx) * nmax_scaling); } ev_idx++; } - ISAMPLE_INFO("\t" << fProjectionCache.size() << "/" << NEvsToProcess - << " events passed selection."); + IEventProcessor_INFO("\t" << fProjectionCache.size() << "/" << NEvsToProcess + << " events passed selection."); } void Write() { if (!fComparisonFinalized) { FinalizeComparison(); } - nuis::persistency::WriteToOutputFile::type>( - fPrediction_comparison.get(), "Prediction", write_directory); - nuis::persistency::WriteToOutputFile::type>( - fPrediction_xsec.get(), "Prediction_xsec", write_directory); + nuis::persistency::WriteToOutputFile( + fPrediction_comparison, "Prediction", write_directory); + nuis::persistency::WriteToOutputFile( + fPrediction_xsec, "Prediction_xsec", write_directory); if (fData) { - nuis::persistency::WriteToOutputFile::type>( - fData.get(), "Data", write_directory); - nuis::persistency::WriteToOutputFile::type>( - fPrediction_shape.get(), "Prediction_shape", write_directory); + nuis::persistency::WriteToOutputFile(fData, "Data", + write_directory); + nuis::persistency::WriteToOutputFile( + fPrediction_shape, "Prediction_shape", write_directory); } } double GetGOF() { if (!fComparisonFinalized) { FinalizeComparison(); } - return nuis::utility::GetChi2(fData.get(), fPrediction_comparison.get()); + if (fData && fPrediction_comparison) { + return nuis::utility::GetChi2(fData, fPrediction_comparison); + } else + return std::numeric_limits::max(); } double GetNDOGuess() { if (fData) { - return nuis::utility::TH_traits< - typename TH_dim_helper::type>::NbinsIncludeFlow(fData.get()); + return TH_Help::Nbins(fData); } return 0; } fhicl::ParameterSet GetExampleConfiguration() { fhicl::ParameterSet exps; exps.put("name", Name()); exps.put("input_type", "Generator"); exps.put("file", "Events.root"); exps.put("write_directory", Name()); exps.put("fake_data_histogram", "/path/to/file.root;h_name"); return exps; } }; typedef SimpleDataComparison<1> SimpleDataComparison_1D; typedef SimpleDataComparison<2> SimpleDataComparison_2D; #endif diff --git a/src/samples/nuA/BubbleChamber/ANL/ANL_CCQE_Evt_1DQ2_nu.cxx b/src/samples/nuA/BubbleChamber/ANL/ANL_CCQE_Evt_1DQ2_nu.cxx index fa99d28..ebfd897 100644 --- a/src/samples/nuA/BubbleChamber/ANL/ANL_CCQE_Evt_1DQ2_nu.cxx +++ b/src/samples/nuA/BubbleChamber/ANL/ANL_CCQE_Evt_1DQ2_nu.cxx @@ -1,192 +1,195 @@ // Copyright 2018 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 "samples/SimpleDataComparison.hxx" #include "utility/FullEventUtility.hxx" #include "utility/KinematicUtility.hxx" #include "utility/PDGCodeUtility.hxx" using namespace nuis::event; using namespace nuis::utility; class ANL_CCQE_Evt_1DQ2_nu : public SimpleDataComparison_1D { public: NEW_NUIS_EXCEPT(invalid_publication_specifier); enum Publication { kPRL31, kPRD16, kPRD26 }; Publication Pub; std::string Pub_str; bool UseD2Corr; - std::unique_ptr fD2CorrHist; - std::unique_ptr fPrediction_Uncorr; + std::unique_ptr fD2CorrHist; + std::unique_ptr fPrediction_Uncorr; ANL_CCQE_Evt_1DQ2_nu() : Pub(kPRD26), Pub_str(""), UseD2Corr(false), fD2CorrHist(nullptr) { ReadGlobalConfigDefaults(); } std::string GetDocumentation() { return "Can specify \"publication: \", where is one of [ PRL31, " "PRD16, PRD26 ] to clarify a publication for comparison. Defaults " "to PRD26.\n" "Can enable deuterium Q2 correction by specifying " "\"use_D2_correction: true\""; } fhicl::ParameterSet GetExampleConfiguration() { fhicl::ParameterSet exps = SimpleDataComparison_1D::GetExampleConfiguration(); exps.put("publication", "PRD26"); exps.put("use_D2_correction", false); return exps; } - void Initialize(fhicl::ParameterSet const &ps) { + void Initialize(fhicl::ParameterSet const &instance_sample_configuration) { - if (ps.has_key("verbosity")) { - SetSampleVerbosity(ps.get("verbosity")); + if (instance_sample_configuration.has_key("verbosity")) { + SetSampleVerbosity( + instance_sample_configuration.get("verbosity")); } - std::string publication = ps.get("publication", "PRD26"); + std::string publication = + instance_sample_configuration.get("publication", "PRD26"); if (publication == "PRL31") { Pub = kPRL31; } else if (publication == "PRD16") { Pub = kPRD16; } else if (publication == "PRD26") { Pub = kPRD26; } else { throw invalid_publication_specifier() << "[ERROR]: Found unexpected publication specifier " << std::quoted(publication) << ". Expected one of [ PRL31, PRD16, PRD26 ]"; } switch (Pub) { case kPRL31: { Pub_str = "PRL31_844"; - EnuRange = std::pair{0, 3E3}; - ISAMPLE_INFO("Sample " << Name() - << " specialized for publication: " << Pub_str); + energy_cut = std::pair{0, 3E3}; + IEventProcessor_INFO( + "Sample " << Name() << " specialized for publication: " << Pub_str); break; } case kPRD16: { Pub_str = "PRD16_3103"; - EnuRange = std::pair{0, 6E3}; - ISAMPLE_INFO("Sample " << Name() - << " specialized for publication: " << Pub_str); + energy_cut = std::pair{0, 6E3}; + IEventProcessor_INFO( + "Sample " << Name() << " specialized for publication: " << Pub_str); break; } case kPRD26: { Pub_str = "PRD26_537"; - EnuRange = std::pair{0, 6E3}; - ISAMPLE_INFO("Sample " << Name() - << " specialized for publication: " << Pub_str); + energy_cut = std::pair{0, 6E3}; + IEventProcessor_INFO( + "Sample " << Name() << " specialized for publication: " << Pub_str); break; } } fhicl::ParameterSet const &global_sample_configuration = nuis::config::GetDocument().get( std::string("global.sample_configuration.") + Name(), fhicl::ParameterSet()); SetData(GetDataDir() + "nuA/BubbleChamber/ANL/CCQE/ANL_CCQE_Data_" + Pub_str + ".root;ANL_1DQ2_Data"); - SimpleDataComparison_1D::Initialize(ps); + SimpleDataComparison_1D::Initialize(instance_sample_configuration); - UseD2Corr = ps.get( + UseD2Corr = instance_sample_configuration.get( "use_D2_correction", global_sample_configuration.get("use_D2_correction", false)); if (UseD2Corr) { - fD2CorrHist = nuis::utility::GetHistogram( + fD2CorrHist = nuis::utility::GetHistogram( GetDataDir() + "nuA/BubbleChamber/ANL/CCQE/" "ANL_CCQE_Data_PRL31_844.root;ANL_1DQ2_Correction"); fPrediction_Uncorr = Clone(fPrediction, true); } // Signal selection function IsSigFunc = [&](FullEvent const &fev) -> bool { if (fev.mode != Channel_t::kCCQE) { return false; } Particle ISNumu = GetHMISNeutralLepton(fev); if (!ISNumu) { return false; } if (ISNumu.pdg != pdgcodes::kNuMu) { return false; } - if ((ISNumu.P4.E() < EnuRange.first) || - (ISNumu.P4.E() > EnuRange.second)) { + if (!energy_cut.IsInRange(ISNumu.P4.E())) { return false; } double Q2 = GetNeutrinoQ2QERec(fev, 0); if (Q2 <= 0) { return false; } return true; }; // 1D Projection function CompProjFunc = [](FullEvent const &fev) -> std::array { - return {{GetNeutrinoQ2QERec(fev, 0)}}; + return {GetNeutrinoQ2QERec(fev, 0)}; }; } // Used to apply D2 correction if requested virtual void FillProjection(std::array const &proj, double event_weight) { if (UseD2Corr) { - nuis::utility::Fill(fPrediction_Uncorr.get(), proj, event_weight); + TH_Help::Fill(fPrediction_Uncorr, proj, event_weight); event_weight *= fD2CorrHist->Interpolate(proj[0]); } - nuis::utility::Fill(fPrediction.get(), proj, event_weight); + TH_Help::Fill(fPrediction, proj, event_weight); } void FinalizeComparison() { SimpleDataComparison_1D::FinalizeComparison(); - fPrediction_Uncorr->Scale(1.0, "width"); + if (UseD2Corr) { + fPrediction_Uncorr->Scale(1.0, "width"); + } } void Write() { SimpleDataComparison_1D::Write(); if (UseD2Corr) { - nuis::persistency::WriteToOutputFile( - fPrediction_Uncorr.get(), "Prediction_Uncorr", write_directory); + nuis::persistency::WriteToOutputFile( + fPrediction_Uncorr, "Prediction_Uncorr", write_directory); } } std::string Name() { return "ANL_CCQE_Evt_1DQ2_nu"; } }; DECLARE_PLUGIN(IDataComparison, ANL_CCQE_Evt_1DQ2_nu); -DECLARE_PLUGIN(ISample, ANL_CCQE_Evt_1DQ2_nu); +DECLARE_PLUGIN(IEventProcessor, ANL_CCQE_Evt_1DQ2_nu); diff --git a/src/samples/nuA/CMakeLists.txt b/src/samples/nuA/CMakeLists.txt index 36dcced..9af875c 100644 --- a/src/samples/nuA/CMakeLists.txt +++ b/src/samples/nuA/CMakeLists.txt @@ -1,4 +1,5 @@ add_subdirectory(BubbleChamber) +add_subdirectory(Nuclear) SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE) diff --git a/src/samples/nuA/CMakeLists.txt b/src/samples/nuA/Nuclear/CMakeLists.txt similarity index 80% copy from src/samples/nuA/CMakeLists.txt copy to src/samples/nuA/Nuclear/CMakeLists.txt index 36dcced..39d211c 100644 --- a/src/samples/nuA/CMakeLists.txt +++ b/src/samples/nuA/Nuclear/CMakeLists.txt @@ -1,4 +1,4 @@ -add_subdirectory(BubbleChamber) +add_subdirectory(T2K) SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE) diff --git a/src/samples/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt b/src/samples/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt new file mode 100644 index 0000000..d62dda3 --- /dev/null +++ b/src/samples/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt @@ -0,0 +1,15 @@ +SET(SAMPLES T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar) + +foreach(S ${SAMPLES}) + LIST(APPEND INuADataComparisons ${S}) + LIST(APPEND SAMPLES_src ${S}.cxx) +endforeach() + +add_library(T2KCC0PiDataComparisons SHARED ${SAMPLES_src}) +target_link_libraries(T2KCC0PiDataComparisons nuis_event nuis_input nuis_config nuis_persistency) + +install(TARGETS T2KCC0PiDataComparisons DESTINATION plugins) + +install(FILES T2K.CC0Pi.sample.config.fcl DESTINATION fcl) + +SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) diff --git a/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K.CC0Pi.sample.config.fcl b/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K.CC0Pi.sample.config.fcl new file mode 100644 index 0000000..fb2a079 --- /dev/null +++ b/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K.CC0Pi.sample.config.fcl @@ -0,0 +1,9 @@ +global.sample_configuration.T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar: { + journal_reference: "Not yet published" + target_material: "H20" + flux_description: "ND280 FHC anti-muon neutrino" + signal_description: "CC0Pi" + shape_only: false + flux_unfolded: false + #has_data: false +} diff --git a/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar.cxx b/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar.cxx new file mode 100644 index 0000000..7576869 --- /dev/null +++ b/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar.cxx @@ -0,0 +1,194 @@ +// Copyright 2018 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 "samples/SimpleDataComparison.hxx" + +#include "utility/EventTopologyUtility.hxx" +#include "utility/FullEventUtility.hxx" +#include "utility/KinematicUtility.hxx" +#include "utility/PDGCodeUtility.hxx" + +using namespace nuis::event; +using namespace nuis::utility; + +using SimpleDataComparison_2DPoly = SimpleDataComparison<2, double, TH2Poly>; + +class T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar + : public SimpleDataComparison_2DPoly { + + std::unique_ptr fPrediction_pmu; + std::unique_ptr fPrediction_ctheta; + + std::unique_ptr fPrediction_fine; + std::unique_ptr fPrediction_unscale; + +public: + T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar() { ReadGlobalConfigDefaults(); } + + std::string GetDocumentation() { return ""; } + fhicl::ParameterSet GetExampleConfiguration() { + fhicl::ParameterSet exps = + SimpleDataComparison_2DPoly::GetExampleConfiguration(); + + return exps; + } + + void Initialize(fhicl::ParameterSet const &instance_sample_configuration) { + + if (instance_sample_configuration.has_key("verbosity")) { + SetSampleVerbosity( + instance_sample_configuration.get("verbosity")); + } + + fhicl::ParameterSet const &global_sample_configuration = + nuis::config::GetDocument().get( + std::string("global.sample_configuration.") + Name(), + fhicl::ParameterSet()); + + SetData( + GetDataDir() + + "nuA/Nuclear/T2K/CC0Pi/H2O_xsec_2Dpmuthetamu_numubar.root;datapoly"); + + fPrediction_fine = std::unique_ptr( + new TH2D("fPrediction_fine", + ";#it{p}_{#mu} " + "(MeV/#it{c});cos(#theta_{#mu});d^{2}#sigma/" + "d#it{p}_{#mu}dcos(#theta_{#mu}) (cm^{2} MeV^{-1} A^{-1})", + 100, 0, 3000, 50, 0.84, 1)); + + fPrediction_pmu = + std::unique_ptr(new TH1D("fPrediction_pmu", + ";#it{p}_{#mu} " + "(MeV/#it{c});d#sigma/" + "d#it{p}_{#mu} (cm^{2} MeV^{-1} A^{-1})", + 100, 0, 3000)); + fPrediction_ctheta = + std::unique_ptr(new TH1D("fPrediction_ctheta", + ";cos(#theta_{#mu});d#sigma/" + "ddcos(#theta_{#mu}) (cm^{2} A^{-1})", + 50, 0.84, 1)); + + SimpleDataComparison_2DPoly::Initialize(instance_sample_configuration); + + fPrediction_unscale = + nuis::utility::Clone(fPrediction, false, "Prediction_unscale"); + + // Signal selection function + IsSigFunc = [](FullEvent const &fev) -> bool { + if (!IsCC0Pi(fev)) { + return false; + } + + Particle ISNumuBar = GetHMISParticle(fev, {pdgcodes::kNuMuBar}); + if (!ISNumuBar) { + return false; + } + + Particle FSMuPlus = GetHMFSParticle(fev, {pdgcodes::kMuPlus}); + if (!FSMuPlus) { + return false; + } + + if (FSMuPlus.CosTheta() < 0.84) { + return false; + } + + return true; + }; + // 1D Projection function + CompProjFunc = [](FullEvent const &fev) -> std::array { + Particle FSMuPlus = GetHMFSParticle(fev, {pdgcodes::kMuPlus}); + return {FSMuPlus.P(), FSMuPlus.CosTheta()}; + }; + + ProcessExtraFunc = [&](FullEvent const &fev, bool isSel, double weight) { + if (isSel) { + Particle FSMuPlus = GetHMFSParticle(fev, {pdgcodes::kMuPlus}); + TH_Helper::Fill(fPrediction_fine, + {FSMuPlus.P(), FSMuPlus.CosTheta()}, weight); + TH_Help::Fill(fPrediction_unscale, {FSMuPlus.P(), FSMuPlus.CosTheta()}); + } + + Particle ISNumuBar = GetHMISParticle(fev, {pdgcodes::kNuMuBar}); + if (!ISNumuBar) { + return; + } + + Particle FSMuPlus = GetHMFSParticle(fev, {pdgcodes::kMuPlus}); + if (!FSMuPlus) { + return; + } + + TH_Helper::Fill(fPrediction_pmu, {FSMuPlus.P()}, weight); + TH_Helper::Fill(fPrediction_ctheta, {FSMuPlus.CosTheta()}, weight); + }; + } + + std::string Name() { return "T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar"; } + + void Write() { + SimpleDataComparison_2DPoly::Write(); + nuis::persistency::WriteToOutputFile( + fPrediction_fine, "Prediction_fine", write_directory); + nuis::persistency::WriteToOutputFile( + fPrediction_pmu, "fPrediction_pmu", write_directory); + nuis::persistency::WriteToOutputFile( + fPrediction_ctheta, "fPrediction_ctheta", write_directory); + + nuis::persistency::WriteToOutputFile( + fPrediction_unscale, "Prediction_unscale", write_directory); + + static std::vector> const + binning = { + {// Slice_0 + YPolyBinSpec(450, 0.85), YPolyBinSpec(450, 0.95)}, + {// Slice_1 + YPolyBinSpec(550, 0.86), YPolyBinSpec(550, 0.93), + YPolyBinSpec(550, 0.97)}, + {// Slice_2 + YPolyBinSpec(700, 0.89), YPolyBinSpec(700, 0.94), + YPolyBinSpec(700, 0.98)}, + {// Slice_3 + YPolyBinSpec(900, 0.91), YPolyBinSpec(900, 0.95), + YPolyBinSpec(900, 0.98)}, + {// Slice_4 + YPolyBinSpec(1200, 0.92), YPolyBinSpec(1200, 0.96), + YPolyBinSpec(1200, 0.98)}, + {// Slice_4 + YPolyBinSpec(1600, 0.93), YPolyBinSpec(1600, 0.97), + YPolyBinSpec(1600, 0.99)}, + {// Slice_5 + YPolyBinSpec(2500, 0.96), YPolyBinSpec(2500, 0.99)}, + }; + + for (auto &slice : GetTH2PolySlices(fData, binning)) { + nuis::persistency::WriteToOutputFile(slice, slice->GetName(), + write_directory); + } + for (auto &slice : GetTH2PolySlices(fPrediction_xsec, binning)) { + nuis::persistency::WriteToOutputFile(slice, slice->GetName(), + write_directory); + } + } +}; + +DECLARE_PLUGIN(IDataComparison, T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar); +DECLARE_PLUGIN(IEventProcessor, T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar); diff --git a/src/samples/nuA/Nuclear/T2K/CMakeLists.txt b/src/samples/nuA/Nuclear/T2K/CMakeLists.txt new file mode 100644 index 0000000..b667c78 --- /dev/null +++ b/src/samples/nuA/Nuclear/T2K/CMakeLists.txt @@ -0,0 +1,7 @@ +add_subdirectory(CC0Pi) + +LIST(APPEND INuADataComparisons_FHiCL T2K.sample.config.fcl) +install(FILES T2K.sample.config.fcl DESTINATION fcl) + +SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) +SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE) diff --git a/src/samples/nuA/Nuclear/T2K/T2K.sample.config.fcl b/src/samples/nuA/Nuclear/T2K/T2K.sample.config.fcl new file mode 100644 index 0000000..18a0fd0 --- /dev/null +++ b/src/samples/nuA/Nuclear/T2K/T2K.sample.config.fcl @@ -0,0 +1 @@ +#include "T2K.CC0Pi.sample.config.fcl" diff --git a/src/utility/CMakeLists.txt b/src/utility/CMakeLists.txt index abde72b..348e8e9 100644 --- a/src/utility/CMakeLists.txt +++ b/src/utility/CMakeLists.txt @@ -1,25 +1,27 @@ set(utility_implementation_files FileSystemUtility.cxx FullEventUtility.cxx KinematicUtility.cxx PDGCodeUtility.cxx - StringUtility.cxx) + StringUtility.cxx + EventTopologyUtility.cxx) set(utility_header_files FileSystemUtility.hxx FullEventUtility.hxx InteractionChannelUtility.hxx KinematicUtility.hxx PDGCodeUtility.hxx ROOTUtility.hxx HistogramUtility.hxx StringUtility.hxx TerminalUtility.hxx - StatsUtility.hxx) + StatsUtility.hxx + EventTopologyUtility.hxx) add_library(nuis_utility SHARED ${utility_implementation_files}) target_link_libraries(nuis_utility) install(TARGETS nuis_utility DESTINATION lib) install(FILES ${utility_header_files} DESTINATION include/utility) diff --git a/src/utility/EventTopologyUtility.cxx b/src/utility/EventTopologyUtility.cxx new file mode 100644 index 0000000..6ba1869 --- /dev/null +++ b/src/utility/EventTopologyUtility.cxx @@ -0,0 +1,33 @@ +#include "utility/EventTopologyUtility.hxx" + +#include "event/FullEvent.hxx" + +#include "utility/FullEventUtility.hxx" + +namespace nuis { +namespace utility { + +bool IsCC0Pi(event::FullEvent const &ev) { + + event::Particle ISNu = GetHMISNeutralLepton(ev); + if (!ISNu) { + return false; + } + + event::PDG_t expected_fslep_pdg = ISNu.pdg > 0 ? ISNu.pdg - 1 : ISNu.pdg + 1; + + + event::Particle FSLep = GetHMFSParticle(ev, {expected_fslep_pdg}); + if (!FSLep) { + return false; + } + + if (GetFSPions(ev).size() || GetFSOthers(ev).size()) { + return false; + } + + return true; +} + +} // namespace utility +} // namespace nuis diff --git a/src/utility/EventTopologyUtility.hxx b/src/utility/EventTopologyUtility.hxx new file mode 100644 index 0000000..2043e9d --- /dev/null +++ b/src/utility/EventTopologyUtility.hxx @@ -0,0 +1,15 @@ +#pragma once + +namespace nuis { +namespace event { +class FullEvent; +} // namespace event +} // namespace nuis + +namespace nuis { +namespace utility { + + bool IsCC0Pi(event::FullEvent const &); + +} +} diff --git a/src/utility/FullEventUtility.cxx b/src/utility/FullEventUtility.cxx index 8188289..1406c6b 100644 --- a/src/utility/FullEventUtility.cxx +++ b/src/utility/FullEventUtility.cxx @@ -1,147 +1,169 @@ #include "utility/FullEventUtility.hxx" #include "event/FullEvent.hxx" #include "utility/PDGCodeUtility.hxx" using namespace nuis::event; namespace nuis { namespace utility { +event::Particle GetHMParticle(std::vector const &parts) { + if (parts.size()) { + return event::Particle(); + } + + return *std::max_element( + parts.begin(), parts.end(), + [](event::Particle const &l, event::Particle const &r) { + return l.P() < r.P(); + }); +} + std::vector GetParticles(FullEvent const &ev, std::vector const &pdgs, Particle::Status_t status, bool include_matching_pdg) { std::vector selected_particles; for (auto const &parts : ev.ParticleStack) { if (parts.status != status) { continue; } for (Particle const &part : parts.particles) { bool matched_pdg = false; for (PDG_t pdg : pdgs) { matched_pdg = matched_pdg || (part.pdg == pdg); } bool keep = ((include_matching_pdg && matched_pdg) || (!include_matching_pdg && !matched_pdg)); if (!keep) { continue; } selected_particles.push_back(part); } } return selected_particles; } std::vector const &GetISParticles(FullEvent const &ev) { return ev .ParticleStack[static_cast( Particle::Status_t::kPrimaryInitialState)] .particles; } std::vector const &GetPrimaryFSParticles(FullEvent const &ev) { return ev .ParticleStack[static_cast( Particle::Status_t::kPrimaryFinalState)] .particles; } std::vector const &GetNuclearLeavingParticles(FullEvent const &ev) { return ev .ParticleStack[static_cast(Particle::Status_t::kNuclearLeaving)] .particles; } Particle GetHMParticle(FullEvent const &ev, std::vector const &pdgs, Particle::Status_t status, bool include_matching_pdg) { Particle HMParticle; for (auto const &parts : ev.ParticleStack) { if (parts.status != status) { continue; } for (Particle const &part : parts.particles) { bool matched_pdg = false; for (PDG_t pdg : pdgs) { matched_pdg = matched_pdg || (part.pdg == pdg); } bool keep = ((include_matching_pdg && matched_pdg) || (!include_matching_pdg && !matched_pdg)); if (!keep) { continue; } - if (part.P4.Vect().Mag() > HMParticle.P4.Vect().Mag()) { + if (part.P() > HMParticle.P()) { HMParticle = part; } } } return HMParticle; } +event::Particle GetHMISParticle(event::FullEvent const &ev, + std::vector const &pdgs) { + return GetHMParticle(ev, pdgs, Particle::Status_t::kPrimaryInitialState); +} + +event::Particle GetHMFSParticle(event::FullEvent const &ev, + std::vector const &pdgs) { + return GetHMParticle(ev, pdgs); +} + std::vector GetFSChargedLeptons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::ChargedLeptons); } std::vector GetFSNeutralLeptons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::NeutralLeptons); } std::vector GetISNeutralLeptons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::NeutralLeptons, Particle::Status_t::kPrimaryInitialState); } std::vector GetFSChargedPions(FullEvent const &ev) { return GetParticles(ev, pdgcodes::ChargedPions); } std::vector GetFSNeutralPions(FullEvent const &ev) { return GetParticles(ev, pdgcodes::NeutralPions); } std::vector GetFSPions(FullEvent const &ev) { return GetParticles(ev, pdgcodes::Pions); } std::vector GetFSProtons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::Protons); } std::vector GetFSNeutrons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::Neutron); } std::vector GetFSNucleons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::Nucleons); } std::vector GetFSOthers(FullEvent const &ev) { return GetParticles(ev, pdgcodes::CommonParticles, Particle::Status_t::kNuclearLeaving, false); } Particle GetHMFSChargedLepton(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::ChargedLeptons); } Particle GetHMFSNeutralLepton(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::NeutralLeptons); } Particle GetHMISNeutralLepton(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::NeutralLeptons, Particle::Status_t::kPrimaryInitialState); } Particle GetHMFSChargedPion(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::ChargedPions); } Particle GetHMFSNeutralPion(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::NeutralPions); } Particle GetHMFSPion(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Pions); } Particle GetHMFSProton(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Protons); } Particle GetHMFSNeutron(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Neutron); } Particle GetHMFSNucleon(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Nucleons); } Particle GetHMFSOther(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::CommonParticles, Particle::Status_t::kNuclearLeaving, false); } } // namespace utility } // namespace nuis diff --git a/src/utility/FullEventUtility.hxx b/src/utility/FullEventUtility.hxx index 7d11a96..aa25bea 100644 --- a/src/utility/FullEventUtility.hxx +++ b/src/utility/FullEventUtility.hxx @@ -1,80 +1,84 @@ // Copyright 2018 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 UTILITY_FULLEVENTUTILITY_HXX_SEEN -#define UTILITY_FULLEVENTUTILITY_HXX_SEEN +#pragma once #include "event/types.hxx" #include "event/Particle.hxx" #include namespace nuis { namespace event { class FullEvent; } // namespace event } // namespace nuis namespace nuis { namespace utility { +event::Particle GetHMParticle(std::vector); + std::vector GetParticles(event::FullEvent const &, std::vector const &, event::Particle::Status_t status = event::Particle::Status_t::kNuclearLeaving, bool include_matching_pdg = true); std::vector const &GetISParticles(event::FullEvent const &); std::vector const & GetPrimaryFSParticles(event::FullEvent const &); std::vector const & GetNuclearLeavingParticles(event::FullEvent const &); event::Particle GetHMParticle(event::FullEvent const &, std::vector const &, event::Particle::Status_t status = event::Particle::Status_t::kNuclearLeaving, bool include_matching_pdg = true); +event::Particle GetHMISParticle(event::FullEvent const &, + std::vector const &); +event::Particle GetHMFSParticle(event::FullEvent const &, + std::vector const &); + std::vector GetFSChargedLeptons(event::FullEvent const &); std::vector GetFSNeutralLeptons(event::FullEvent const &); std::vector GetISNeutralLeptons(event::FullEvent const &); std::vector GetFSChargedPions(event::FullEvent const &); std::vector GetFSNeutralPions(event::FullEvent const &); std::vector GetFSPions(event::FullEvent const &); std::vector GetFSProtons(event::FullEvent const &); std::vector GetFSNeutrons(event::FullEvent const &); std::vector GetFSNucleons(event::FullEvent const &); std::vector GetFSOthers(event::FullEvent const &); event::Particle GetHMFSChargedLepton(event::FullEvent const &); event::Particle GetHMFSNeutralLepton(event::FullEvent const &); event::Particle GetHMISNeutralLepton(event::FullEvent const &); event::Particle GetHMFSChargedPion(event::FullEvent const &); event::Particle GetHMFSNeutralPion(event::FullEvent const &); event::Particle GetHMFSPion(event::FullEvent const &); event::Particle GetHMFSProton(event::FullEvent const &); event::Particle GetHMFSNeutron(event::FullEvent const &); event::Particle GetHMFSNucleon(event::FullEvent const &); event::Particle GetHMFSOther(event::FullEvent const &); } // namespace utility } // namespace nuis - -#endif diff --git a/src/utility/HistogramUtility.hxx b/src/utility/HistogramUtility.hxx index 34e99bc..b3f8ef4 100644 --- a/src/utility/HistogramUtility.hxx +++ b/src/utility/HistogramUtility.hxx @@ -1,61 +1,519 @@ // Copyright 2018 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 UTILITY_HISTOGRAMUTILITY_HXX_SEEN #define UTILITY_HISTOGRAMUTILITY_HXX_SEEN #include "utility/ROOTUtility.hxx" #include "exception/exception.hxx" #include "string_parsers/from_string.hxx" +#include "TAxis.h" +#include "TH1D.h" +#include "TH1F.h" +#include "TH2D.h" +#include "TH2F.h" +#include "TH2Poly.h" + +#include #include #include #include -#include namespace nuis { namespace utility { NEW_NUIS_EXCEPT(unimplemented_GetHistogram_method); NEW_NUIS_EXCEPT(invalid_histogram_descriptor); +NEW_NUIS_EXCEPT(invalid_histogram_name); +NEW_NUIS_EXCEPT(failed_to_clone); + +inline bool IsFlowBin(TAxis const &ax, Int_t bin_it) { + return ((bin_it <= 0) || (bin_it >= (ax.GetNbins() + 1))); +} + +inline bool IsInHistogramRange(TAxis const &ax, double v) { + Int_t bin_it = ax.FindFixBin(v); + return !IsFlowBin(ax, bin_it); +} + +template struct HType_traits {}; +template <> struct HType_traits { + using type = TH1; + static size_t const NDim = 1; + using NumericT = double; + static std::string name() { return "TH1"; } +}; +template <> struct HType_traits { + using type = TH1D; + static size_t const NDim = 1; + using NumericT = double; + static std::string name() { return "TH1D"; } +}; +template <> struct HType_traits { + using type = TH1F; + static size_t const NDim = 1; + using NumericT = float; + static std::string name() { return "TH1F"; } +}; +template <> struct HType_traits { + using type = TH2; + static size_t const NDim = 2; + using NumericT = double; + static std::string name() { return "TH2"; } +}; +template <> struct HType_traits { + using type = TH2D; + static size_t const NDim = 2; + using NumericT = double; + static std::string name() { return "TH2D"; } +}; +template <> struct HType_traits { + using type = TH2F; + static size_t const NDim = 2; + using NumericT = float; + static std::string name() { return "TH2F"; } +}; +template <> struct HType_traits { + using type = TH2Poly; + static size_t const NDim = 0; + using NumericT = double; + static std::string name() { return "TH2Poly"; } +}; + +template struct HType_Helper {}; +template <> struct HType_Helper<1, void> { + using type = TH1; + static size_t const NDim = HType_traits::NDim; + using NumericT = HType_traits::NumericT; + static std::string name() { return HType_traits::name(); } +}; +template <> struct HType_Helper<1, double> { + using type = TH1D; + static size_t const NDim = HType_traits::NDim; + using NumericT = HType_traits::NumericT; + static std::string name() { return HType_traits::name(); } +}; +template <> struct HType_Helper<1, float> { + using type = TH1F; + static size_t const NDim = HType_traits::NDim; + using NumericT = HType_traits::NumericT; + static std::string name() { return HType_traits::name(); } +}; +template <> struct HType_Helper<2, void> { + using type = TH2; + static size_t const NDim = HType_traits::NDim; + using NumericT = HType_traits::NumericT; + static std::string name() { return HType_traits::name(); } +}; +template <> struct HType_Helper<2, double> { + using type = TH2D; + static size_t const NDim = HType_traits::NDim; + using NumericT = HType_traits::NumericT; + static std::string name() { return HType_traits::name(); } +}; +template <> struct HType_Helper<2, float> { + using type = TH2F; + static size_t const NDim = HType_traits::NDim; + using NumericT = HType_traits::NumericT; + static std::string name() { return HType_traits::name(); } +}; + +template struct TH_Helper {}; + +template +struct TH_Helper::NDim == 1>::type> { + static size_t const NDim = HType_traits::NDim; + using NumericT = typename HType_traits::NumericT; + static std::string name() { return HType_traits::name(); } + + static Int_t NbinsIncludeFlow(HT const &h) { + return h.GetXaxis()->GetNbins() + 2; + } + static Int_t Nbins(HT const &h) { return h.GetXaxis()->GetNbins(); } + + static bool IsFlowBin(HT const &h, Int_t bin_it) { + return nuis::utility::IsFlowBin(*h.GetXaxis(), bin_it); + } + + static void Fill(HT &h, std::array const &v, + double w = 1) { + h.Fill(v[0], w); + } + + static void Scale(HT &h, NumericT SF, char const *opt = "") { + h.Scale(SF, opt); + } + + static double Integral(HT const &h, char const *opt = "") { + return h.Integral(opt); + } + + static Int_t NbinsIncludeFlow(std::unique_ptr const &h) { + return NbinsIncludeFlow(*h); + } + static Int_t Nbins(std::unique_ptr const &h) { return Nbins(*h); } + + static bool IsFlowBin(std::unique_ptr const &h, Int_t bin_it) { + return IsFlowBin(*h, bin_it); + } + + static void Fill(std::unique_ptr &h, + std::array const &v, double w = 1) { + Fill(*h, v, w); + } + + static void Scale(std::unique_ptr &h, NumericT SF, char const *opt = "") { + Scale(*h, SF, opt); + } + + static double Integral(std::unique_ptr const &h, char const *opt = "") { + return Integral(*h, opt); + } +}; + +template +struct TH_Helper::NDim == 2>::type> { + static size_t const NDim = HType_traits::NDim; + using NumericT = typename HType_traits::NumericT; + static std::string name() { return HType_traits::name(); } + + // TH2 *************************************************************** + static Int_t NbinsIncludeFlow(HT const &h) { + return (h.GetXaxis()->GetNbins() + 2) * (h.GetYaxis()->GetNbins() + 2); + } + static Int_t Nbins(HT const &h) { + return (h.GetXaxis()->GetNbins()) * (h.GetYaxis()->GetNbins()); + } + + static bool IsFlowBin(HT const &h, Int_t xbin_it, Int_t ybin_it) { + return nuis::utility::IsFlowBin(*h.GetXaxis(), xbin_it) || + nuis::utility::IsFlowBin(*h.GetYaxis(), ybin_it); + } + + static void Fill(HT &h, std::array const &v, + double w = 1) { + h.Fill(v[0], v[1], w); + } + + static void Scale(HT &h, NumericT SF, char const *opt = "") { + h.Scale(SF, opt); + } + + static double Integral(HT const &h, char const *opt = "") { + return h.Integral(opt); + } + + static Int_t NbinsIncludeFlow(std::unique_ptr const &h) { + return NbinsIncludeFlow(*h); + } + static Int_t Nbins(std::unique_ptr const &h) { return Nbins(*h); } + + static bool IsFlowBin(std::unique_ptr const &h, Int_t xbin_it, + Int_t ybin_it) { + return IsFlowBin(*h, xbin_it, ybin_it); + } + + static void Fill(std::unique_ptr &h, + std::array const &v, double w = 1) { + Fill(*h, v, w); + } + + static void Scale(std::unique_ptr &h, NumericT SF, char const *opt = "") { + Scale(*h, SF, opt); + } + + static double Integral(std::unique_ptr const &h, char const *opt = "") { + return Integral(*h, opt); + } +}; + +template +struct TH_Helper< + HT, typename std::enable_if::value>::type> { + static size_t const NDim = 2; + using NumericT = typename HType_traits::NumericT; + static std::string name() { return HType_traits::name(); } + + // TH2Poly *************************************************************** + static Int_t NbinsIncludeFlow(HT const &h) { return h.GetNumberOfBins() + 9; } + static Int_t Nbins(HT const &h) { return h.GetNumberOfBins(); } + + static bool IsFlowBin(HT const &h, Int_t bin_it) { return (bin_it < 0); } + + static void Fill(HT &h, std::array const &v, + double w = 1) { + h.Fill(v[0], v[1], w); + } + + static void Scale(HT &h, NumericT SF, char const *opt = "") { + + bool width = (std::string(opt).find("width") != std::string::npos); + size_t nbins = Nbins(h); + for (size_t bin_it = 0; bin_it < nbins; ++bin_it) { + double bin_area = 1; + + if (width) { + TH2PolyBin *poly_bin = + dynamic_cast(h.GetBins()->At(bin_it)); + + bin_area = poly_bin->GetArea(); + } + + h.SetBinContent(bin_it + 1, + h.GetBinContent(bin_it + 1) * (SF / bin_area)); + h.SetBinError(bin_it + 1, h.GetBinError(bin_it + 1) * (SF / bin_area)); + } + } + + static double Integral(HT &h, char const *opt = "") { + bool width = (std::string(opt).find("width") != std::string::npos); + size_t nbins = Nbins(h); + double integral = 0; + for (size_t bin_it = 0; bin_it < nbins; ++bin_it) { + double bin_area = 1; + + if (width) { + TH2PolyBin *poly_bin = + dynamic_cast(h.GetBins()->At(bin_it)); + + bin_area = poly_bin->GetArea(); + } + + integral += h.GetBinContent(bin_it + 1) * bin_area; + } + return integral; + } + + static Int_t NbinsIncludeFlow(std::unique_ptr const &h) { + return NbinsIncludeFlow(*h); + } + static Int_t Nbins(std::unique_ptr const &h) { return Nbins(*h); } + + static bool IsFlowBin(std::unique_ptr const &h, Int_t bin_it) { + return IsFlowBin(*h, bin_it); + } + + static void Fill(std::unique_ptr &h, + std::array const &v, double w = 1) { + Fill(*h, v, w); + } + + static void Scale(std::unique_ptr &h, NumericT SF, char const *opt = "") { + Scale(*h, SF, opt); + } + + static double Integral(std::unique_ptr &h, char const *opt = "") { + return Integral(*h, opt); + } +}; + +template +void Clear(typename std::enable_if::NDim != 0, HT>::type &h) { + for (Int_t bin_it = 0; bin_it < TH_Helper::NbinsIncludeFlow(h); + ++bin_it) { + h.SetBinContent(bin_it, 0); + h.SetBinError(bin_it, 0); + } +} + +template +void Clear( + typename std::enable_if::value, HT>::type &h) { + h.ClearBinContents(); +} + +template +inline std::unique_ptr Clone(HT const &source, bool clear = false, + std::string const &clone_name = "") { + std::unique_ptr target(dynamic_cast( + source.Clone(clone_name.size() ? clone_name.c_str() : ""))); + if (!target) { + throw failed_to_clone() + << "[ERROR]: Failed to clone a " << TH_Helper::name() << "."; + } + target->SetDirectory(nullptr); + + if (clear) { + Clear(*target); + } + + return target; +} +template +inline std::unique_ptr Clone(std::unique_ptr const &source, + bool clear = false, + std::string const &clone_name = "") { + return Clone(*source, clear, clone_name); +} + +template +inline std::unique_ptr GetHistogramFromROOTFile(TFile_ptr &f, + std::string const &hname) { + + f->Get(hname.c_str()); + HT *h = dynamic_cast(f->Get(hname.c_str())); + if (!h) { + throw invalid_histogram_name() + << "[ERROR]: Failed to get " << TH_Helper::name() << " named " + << std::quoted(hname) << " from input file " + << std::quoted(f->GetName()); + } + std::unique_ptr clone = Clone(*h); + return clone; +} + +template +inline std::unique_ptr GetHistogramFromROOTFile(std::string const &fname, + std::string const &hname) { + TFile_ptr temp = CheckOpenTFile(fname, "READ"); + return GetHistogramFromROOTFile(temp, hname); +} template std::unique_ptr GetHistogram(std::string const &input_descriptor) { std::vector split_descriptor = fhicl::string_parsers::ParseToVect(input_descriptor, ";", true, true); if (split_descriptor.size() == 1) { // assume text throw unimplemented_GetHistogram_method(); } else if (split_descriptor.size() == 2) { return GetHistogramFromROOTFile(split_descriptor[0], split_descriptor[1]); } else { throw invalid_histogram_descriptor() << "[ERROR]: Failed to parse histogram descriptor: " << std::quoted(input_descriptor) << " as an input histogram (Text/ROOT)."; } } +NEW_NUIS_EXCEPT(invalid_TH2Poly); +NEW_NUIS_EXCEPT(invalid_PolyBinSpecifierList); + +// List of bins to put in a row, X/Y +struct PolyBinSpecifier { + double X, Y; + bool UseXAxis; +}; + +constexpr PolyBinSpecifier XPolyBinSpec(double X, double Y) { + return PolyBinSpecifier{X + std::numeric_limits::epsilon() * 1E2, + Y + std::numeric_limits::epsilon() * 1E2, + true}; +} + +constexpr PolyBinSpecifier YPolyBinSpec(double X, double Y) { + return PolyBinSpecifier{X + std::numeric_limits::epsilon() * 1E2, + Y + std::numeric_limits::epsilon() * 1E2, + false}; +} + +std::vector> GetTH2PolySlices( + std::unique_ptr &hinp, + std::vector> const &BinsSpecifiers) { + + std::vector> slices; + + size_t sl_it = 0; + for (auto &slice_spec : BinsSpecifiers) { + std::vector Binning; + std::vector BinContent; + std::vector BinError; + bool UseXAxis = false; + size_t bin_ctr = 0; + for (auto poly_bin_spec : slice_spec) { + Int_t bin_it = hinp->FindBin(poly_bin_spec.X, poly_bin_spec.Y); + if (bin_it < 1) { + std::cout << "[WARN]: When searching for matching bin: { X: " + << poly_bin_spec.X << ", Y: " << poly_bin_spec.Y + << "} got flow bin: " << bin_it << std::endl; + continue; + } + TH2PolyBin *poly_bin = + dynamic_cast(hinp->GetBins()->At(bin_it - 1)); + + if (!bin_ctr) { + UseXAxis = poly_bin_spec.UseXAxis; + } else if (UseXAxis != poly_bin_spec.UseXAxis) { + throw invalid_PolyBinSpecifierList() + << "[ERROR]: For slice: " << sl_it + << " of TH2Poly: " << std::quoted(hinp->GetName()) + << " bin specifier: " << bin_ctr << " was set to use the " + << (poly_bin_spec.UseXAxis ? "X" : "Y") + << " axis as the dependent axis of the slice, but previous bins " + "were set to use the " + << (poly_bin_spec.UseXAxis ? "X" : "Y") << " axis."; + } + + if ((poly_bin_spec.X < poly_bin->GetXMin()) || + (poly_bin_spec.X >= poly_bin->GetXMax()) || + (poly_bin_spec.Y < poly_bin->GetYMin()) || + (poly_bin_spec.Y >= poly_bin->GetYMax())) { + std::cout + << "[WARN]: Found bin doesn't seem to contain expected point: { X: " + << poly_bin_spec.X << ", Y: " << poly_bin_spec.Y + << "}, got bin_it = " << bin_it + << " which had x_low: " << poly_bin->GetXMin() + << ", and x_up: " << poly_bin->GetXMax() + << ", y_low: " << poly_bin->GetYMin() + << ", y_up: " << poly_bin->GetYMax() << std::endl; + } + + double low = UseXAxis ? poly_bin->GetXMin() : poly_bin->GetYMin(); + double up = UseXAxis ? poly_bin->GetXMax() : poly_bin->GetYMax(); + + if (!Binning.size()) { // Add low edge + Binning.push_back(low); + } else if (std::abs(Binning.back() - low) > + (std::numeric_limits::epsilon() * 1E2)) { + BinContent.push_back(0); + BinError.push_back(0); + Binning.push_back(low); + } + + BinContent.push_back(hinp->GetBinContent(bin_it)); + BinError.push_back(hinp->GetBinError(bin_it)); + Binning.push_back(up); + } + slices.emplace_back(new TH1D( + (std::string(hinp->GetName()) + "_slice" + std::to_string(sl_it++)) + .c_str(), + (std::string(";") + + (UseXAxis ? hinp->GetXaxis() : hinp->GetYaxis())->GetTitle() + ";" + + hinp->GetZaxis()->GetTitle()) + .c_str(), + Binning.size() - 1, Binning.data())); + + for (size_t bin_it = 0; bin_it < BinContent.size(); ++bin_it) { + slices.back()->SetBinContent(bin_it + 1, BinContent[bin_it]); + slices.back()->SetBinError(bin_it + 1, BinError[bin_it]); + } + } + return slices; +} + } // namespace utility } // namespace nuis #endif diff --git a/src/utility/KinematicUtility.hxx b/src/utility/KinematicUtility.hxx index f9e698d..5143d44 100644 --- a/src/utility/KinematicUtility.hxx +++ b/src/utility/KinematicUtility.hxx @@ -1,40 +1,53 @@ // Copyright 2018 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 UTILITY_KINEMATICUTILITY_HXX_SEEN #define UTILITY_KINEMATICUTILITY_HXX_SEEN +#include +#include + namespace nuis { namespace event { class FullEvent; } } // namespace nuis namespace nuis { namespace utility { double GetNeutrinoEQERec(event::FullEvent const &fev, double SeparationEnergy_MeV); double GetNeutrinoQ2QERec(event::FullEvent const &fev, double SeparationEnergy_MeV); +struct ENuRange : public std::pair { + using std::pair::pair; + using std::pair::operator=; + + ENuRange() + : std::pair(0, std::numeric_limits::max()) {} + + bool IsInRange(double enu) { return (enu > second) || (enu < first); } +}; + } // namespace utility } // namespace nuis #endif diff --git a/src/utility/ROOTUtility.hxx b/src/utility/ROOTUtility.hxx index 34231f3..1eef043 100644 --- a/src/utility/ROOTUtility.hxx +++ b/src/utility/ROOTUtility.hxx @@ -1,211 +1,142 @@ // Copyright 2018 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 UTILITY_ROOTUTILITY_HXX_SEEN #define UTILITY_ROOTUTILITY_HXX_SEEN #include "exception/exception.hxx" #include "TFile.h" -#include "TH1D.h" -#include "TH1F.h" -#include "TH2D.h" -#include "TH2F.h" #include "TTree.h" #include #include #include #include namespace nuis { namespace utility { NEW_NUIS_EXCEPT(failed_to_open_TFile); NEW_NUIS_EXCEPT(failed_to_get_TTree); -NEW_NUIS_EXCEPT(invalid_histogram_name); -NEW_NUIS_EXCEPT(failed_to_clone); -inline TFile *CheckOpenTFile(std::string const &fname, char const *opts = "") { +void CloseTFile(TFile *f = nullptr) { + if (f) { + std::cout << "[INFO]: Shutting TFile: " << f->GetName() << ", " << f + << std::endl; + f->Close(); + } + delete f; +} + +using TFile_ptr = std::unique_ptr; +inline TFile_ptr make_unique_TFile(TFile *f) { + return TFile_ptr(f, &CloseTFile); +} + +inline TFile_ptr CheckOpenTFile(std::string const &fname, + char const *opts = "") { TDirectory *ogDir = gDirectory; TFile *inpF = new TFile(fname.c_str(), opts); if (!inpF || !inpF->IsOpen()) { throw failed_to_open_TFile() << "[ERROR]: Couldn't open input file: " << std::quoted(fname); } if (ogDir) { ogDir->cd(); } - return inpF; + return make_unique_TFile(inpF); } struct TreeFile { - TFile *file; + TFile_ptr file; TTree *tree; - TreeFile() : file(nullptr), tree(nullptr) {} + bool file_owned; + TreeFile() + : file(make_unique_TFile(nullptr)), tree(nullptr), file_owned(false) {} TreeFile(TreeFile const &other) = delete; - TreeFile(TreeFile &&other) : file(other.file), tree(other.tree) { + TreeFile &operator=(TreeFile &&other) { + file = std::move(other.file); + tree = other.tree; + file_owned = other.file_owned; + + other.file = nullptr; + other.tree = nullptr; + other.file_owned = false; + + return *this; + } + + TreeFile(TreeFile &&other) + : file(std::move(other.file)), tree(other.tree), + file_owned(other.file_owned) { other.file = nullptr; other.tree = nullptr; + other.file_owned = false; } ~TreeFile() { - if (file) { - file->Close(); + if (!file_owned) { + file.release(); } } }; -inline std::unique_ptr CheckGetTTree(std::string const &fname, - std::string const &tname, - char const *opts = "") { +inline TreeFile CheckGetTTree(TFile *file, std::string const &tname) { TreeFile tf; - tf.file = CheckOpenTFile(fname, opts); + tf.file = make_unique_TFile(file); tf.tree = dynamic_cast(tf.file->Get(tname.c_str())); + tf.file_owned = false; + if (!tf.tree) { throw failed_to_get_TTree() << "[ERROR]: Failed to get TTree named: " << std::quoted(tname) - << " from TFile: " << std::quoted(fname); + << " from TFile: " << std::quoted(file->GetName()); } - return std::make_unique(std::move(tf)); + std::cout << "[INFO]: Opened TFile: " << file->GetName() << ", " << file + << std::endl; + return tf; +} + +inline TreeFile CheckGetTTree(std::string const &fname, + std::string const &tname, char const *opts = "") { + TreeFile tf = CheckGetTTree(CheckOpenTFile(fname, opts).release(), tname); + tf.file_owned = true; + return tf; } -inline std::unique_ptr MakeNewTTree(std::string const &fname, - std::string const &tname, - char const *opts = "") { +inline TreeFile MakeNewTTree(std::string const &fname, std::string const &tname, + char const *opts = "") { TreeFile tf; tf.file = CheckOpenTFile(fname, opts); tf.tree = new TTree(tname.c_str(), ""); - tf.tree->SetDirectory(tf.file); - return std::make_unique(std::move(tf)); -} - -template struct TH_traits {}; - -template <> struct TH_traits { - static size_t const NDims = 1; - static std::string name() { return "TH1"; } - static Int_t NbinsIncludeFlow(TH1 *const h) { - return h->GetXaxis()->GetNbins() + 2; - } -}; -template <> struct TH_traits { - static size_t const NDims = 1; - static std::string name() { return "TH1D"; } - static Int_t NbinsIncludeFlow(TH1D *const h) { - return h->GetXaxis()->GetNbins() + 2; - } -}; -template <> struct TH_traits { - static size_t const NDims = 1; - static std::string name() { return "TH1F"; } - static Int_t NbinsIncludeFlow(TH1F *const h) { - return h->GetXaxis()->GetNbins() + 2; - } -}; -template <> struct TH_traits { - static size_t const NDims = 2; - static std::string name() { return "TH2"; } - static Int_t NbinsIncludeFlow(TH2 *const h) { - return (h->GetXaxis()->GetNbins() + 2) * (h->GetYaxis()->GetNbins() + 2); - } -}; -template <> struct TH_traits { - static size_t const NDims = 2; - static std::string name() { return "TH2D"; } - static Int_t NbinsIncludeFlow(TH2D *const h) { - return (h->GetXaxis()->GetNbins() + 2) * (h->GetYaxis()->GetNbins() + 2); - } -}; -template <> struct TH_traits { - static size_t const NDims = 2; - static std::string name() { return "TH2F"; } - static Int_t NbinsIncludeFlow(TH2F *const h) { - return (h->GetXaxis()->GetNbins() + 2) * (h->GetYaxis()->GetNbins() + 2); - } -}; - -template -inline std::unique_ptr GetHistogramFromROOTFile(std::string const &fname, - std::string const &hname) { - TFile *f = CheckOpenTFile(fname, "READ"); - HT *h = dynamic_cast(f->Get(hname.c_str())); - if (!h) { - throw invalid_histogram_name() - << "[ERROR]: Failed to get " << TH_traits::name() << " named " - << std::quoted(hname) << " from input file " << std::quoted(fname); - } - std::unique_ptr clone(dynamic_cast(h->Clone())); - clone->SetDirectory(nullptr); - - f->Close(); - delete f; - - return clone; + tf.tree->SetDirectory(tf.file.get()); + tf.file_owned = true; + return tf; } -template void Clear(HT *h) { - for (Int_t bin_it = 0; bin_it < TH_traits::NbinsIncludeFlow(h); - ++bin_it) { - h->SetBinContent(bin_it, 0); - h->SetBinError(bin_it, 0); - } -} - -// Fill for 1D histograms uses type system to enforce correct number of -// dimensions -template -typename std::enable_if::NDims == 1, void>::type -Fill(HT *h, std::array::NDims> const &val, double weight = 1) { - h->Fill(val[0], weight); -} - -/// Fill for 2D histograms uses type system to enforce correct number of -/// dimensions -template -typename std::enable_if::NDims == 2, void>::type -Fill(HT *h, std::array::NDims> const &val, double weight = 1) { - h->Fill(val[0], val[1], weight); -} - -template -inline std::unique_ptr Clone(std::unique_ptr const &source, - bool clear = false) { - std::unique_ptr target(dynamic_cast(source->Clone())); - if (!target) { - throw failed_to_clone() - << "[ERROR]: Failed to clone a " << TH_traits::name() - << ", source = " << source.get(); - } - target->SetDirectory(nullptr); - - if (clear) { - Clear(target.get()); - } - - return target; -} +std::string SanitizeROOTObjectName(std::string name) { return name; } } // namespace utility } // namespace nuis #endif diff --git a/src/utility/StatsUtility.hxx b/src/utility/StatsUtility.hxx index 4e207c4..8c559d9 100644 --- a/src/utility/StatsUtility.hxx +++ b/src/utility/StatsUtility.hxx @@ -1,45 +1,66 @@ #ifndef UTILITY_STATSUTILITY_HXX_SEEN #define UTILITY_STATSUTILITY_HXX_SEEN +#include "utility/HistogramUtility.hxx" + namespace nuis { namespace utility { NEW_NUIS_EXCEPT(unimplemented_covariance_usage); NEW_NUIS_EXCEPT(Mismatched_NBins); template -typename std::enable_if::NDims == 1, double>::type -GetChi2(HT const *a, HT const *b, TH2D const *Covariance = nullptr) { +double GetChi2(typename std::enable_if::NDim == 1, + std::unique_ptr>::type const &a, + std::unique_ptr const &b, + std::unique_ptr const &Covariance = nullptr) { if (Covariance) { throw unimplemented_covariance_usage() << "[ERROR]: Using a covariance in the Chi2 evaluation is not yet " "implemented."; } if (a->GetXaxis()->GetNbins() != b->GetXaxis()->GetNbins()) { Mismatched_NBins() << "[ERROR]: Attempted to evaluate Chi2 between two " "histograms with differing bin contents: NBins(a) = " << a->GetXaxis()->GetNbins() << ", NBins(b) = " << b->GetXaxis()->GetNbins(); } double chi2 = 0; for (Int_t bin_it_i = 0; bin_it_i < a->GetXaxis()->GetNbins(); ++bin_it_i) { for (Int_t bin_it_j = 0; bin_it_j < a->GetXaxis()->GetNbins(); ++bin_it_j) { double err = 1.0 / (b->GetBinError(bin_it_i + 1) * b->GetBinError(bin_it_j + 1)); double contrib = (a->GetBinContent(bin_it_i + 1) - b->GetBinContent(bin_it_i + 1)) * (err) * (a->GetBinContent(bin_it_j + 1) - b->GetBinContent(bin_it_j + 1)); chi2 += contrib; } } return chi2; } + +template +double GetChi2(typename std::enable_if::NDim == 2, + std::unique_ptr>::type const &a, + std::unique_ptr const &b, + std::unique_ptr const &Covariance = nullptr) { + return std::numeric_limits::max(); +} + +template +double GetChi2(typename std::enable_if::value, + std::unique_ptr>::type const &a, + std::unique_ptr const &b, + std::unique_ptr const &Covariance = nullptr) { + return std::numeric_limits::max(); +} + } // namespace utility } // namespace nuis #endif diff --git a/src/variation/CMakeLists.txt b/src/variation/CMakeLists.txt index e69de29..fe425d0 100644 --- a/src/variation/CMakeLists.txt +++ b/src/variation/CMakeLists.txt @@ -0,0 +1,14 @@ +SET(variation_implementation_files + IWeightProvider.cxx + WeightManager.cxx) + +SET(variation_header_files + WeightManager.hxx + IVariationProvider.hxx + IWeightProvider.hxx) + +add_library(nuis_variation SHARED ${variation_implementation_files}) +target_link_libraries(nuis_variation nuis_event nuis_plugins) + +install(TARGETS nuis_variation DESTINATION lib) +install(FILES ${variation_header_files} DESTINATION include/variation) diff --git a/src/variation/IVariationProvider.hxx b/src/variation/IVariationProvider.hxx index a7ccbb4..3c4a577 100644 --- a/src/variation/IVariationProvider.hxx +++ b/src/variation/IVariationProvider.hxx @@ -1,60 +1,55 @@ // Copyright 2018 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 VARIATION_IVARIATIONPROVIDER_HXX_SEEN -#define VARIATION_IVARIATIONPROVIDER_HXX_SEEN +#pragma once #include "plugins/traits.hxx" -#include "parameters/ParameterManager.hxx" - namespace fhicl { class ParameterSet; } namespace nuis { namespace event { class MinimalEvent; class FullEvent; } // namespace event } // namespace nuis class IVariationProvider { public: virtual void Initialize(fhicl::ParameterSet const &) = 0; - virtual nuis::params::paramId_t GetParameterId(std::string const &) = 0; - - bool HandlesParameter(std::string const ¶m_name) { - return (GetParameterId(param_name) != nuis::params::kParamUnhandled); - } - - virtual void SetParameterValue(nuis::params::paramId_t, double) = 0; - virtual bool ParametersVaried() = 0; virtual void Reconfigure() = 0; virtual nuis::event::FullEvent VaryFullEvent(nuis::event::FullEvent const &) = 0; + virtual std::string GetName() = 0; + + virtual std::string GetDocumentation() { + return "No documentation for IVariationProvider: " + GetName(); + } + + virtual fhicl::ParameterSet GetExampleConfiguration() = 0; + virtual ~IVariationProvider() {} }; DECLARE_PLUGIN_INTERFACE(IVariationProvider); - -#endif diff --git a/src/variation/IWeightProvider.cxx b/src/variation/IWeightProvider.cxx new file mode 100644 index 0000000..473f5e4 --- /dev/null +++ b/src/variation/IWeightProvider.cxx @@ -0,0 +1,10 @@ +#include "variation/IWeightProvider.hxx" + +#include "event/FullEvent.hxx" + +nuis::event::FullEvent +IWeightProvider::VaryFullEvent(nuis::event::FullEvent const &fe) { + nuis::event::FullEvent fe_clone = fe.Clone(); + fe_clone.RWWeight = GetEventWeight(fe_clone); + return fe_clone; +} diff --git a/src/variation/IWeightProvider.hxx b/src/variation/IWeightProvider.hxx index 2feabbb..0c89690 100644 --- a/src/variation/IWeightProvider.hxx +++ b/src/variation/IWeightProvider.hxx @@ -1,41 +1,34 @@ // Copyright 2018 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 VARIATION_IWEIGHTPROVIDER_HXX_SEEN -#define VARIATION_IWEIGHTPROVIDER_HXX_SEEN +#pragma once #include "variation/IVariationProvider.hxx" class IWeightProvider : public IVariationProvider { public: virtual double GetEventWeight(nuis::event::MinimalEvent const &) = 0; /// For weight providers, the full variation is just the application of the reweight_weight. - nuis::event::FullEvent VaryFullEvent(nuis::event::FullEvent const &fe) { - nuis::event::FullEvent fe_clone = fe.clone(); - fe_clone.RWWeight = GetEventWeight(fe_clone); - return fe_clone; - } + nuis::event::FullEvent VaryFullEvent(nuis::event::FullEvent const &); virtual ~IWeightProvider() {} }; DECLARE_PLUGIN_INTERFACE(IWeightProvider); - -#endif diff --git a/src/variation/VariationManager.cxx b/src/variation/VariationManager.cxx deleted file mode 100644 index e69de29..0000000 diff --git a/src/variation/WeightManager.cxx b/src/variation/WeightManager.cxx new file mode 100644 index 0000000..e64e259 --- /dev/null +++ b/src/variation/WeightManager.cxx @@ -0,0 +1,50 @@ +#include "variation/WeightManager.hxx" + +#include "fhiclcpp/ParameterSet.h" + +#include "plugins/Instantiate.hxx" + +namespace nuis { +namespace variation { + +WeightManager::NamedWeightProvider::NamedWeightProvider( + std::string const &_name, + plugins::plugin_traits::unique_ptr_t &&_handler) + : name(_name), handler(std::move(_handler)) {} + +WeightManager::NamedWeightProvider::NamedWeightProvider( + NamedWeightProvider &&other) + : name(std::move(other.name)), handler(std::move(other.handler)) {} + +WeightManager *WeightManager::_global_inst = nullptr; + +WeightManager::WeightManager() {} + +WeightManager &WeightManager::Get() { + if (!_global_inst) { + _global_inst = new WeightManager(); + } + return *_global_inst; +} + +WeightManager::WeightProv_id_t +WeightManager::EnsureWeightProviderLoaded(fhicl::ParameterSet const &ps) { + std::string const &prov_name = ps.get("weight_engine_name"); + for (size_t i = 0; i < WeightEngines.size(); ++i) { + if (WeightEngines[i].name == prov_name) { + return i; + } + } + + WeightProv_id_t wid = WeightEngines.size(); + WeightEngines.emplace_back( + prov_name, nuis::plugins::Instantiate(prov_name)); + WeightEngines.back().handler->Initialize(ps); + + return wid; +} +double WeightManager::GetEventWeight(nuis::event::MinimalEvent const &) { + return 1; +} +} // namespace variation +} // namespace nuis diff --git a/src/variation/VariationManager.hxx b/src/variation/WeightManager.hxx similarity index 75% rename from src/variation/VariationManager.hxx rename to src/variation/WeightManager.hxx index 79fac54..4221f06 100644 --- a/src/variation/VariationManager.hxx +++ b/src/variation/WeightManager.hxx @@ -1,69 +1,68 @@ // Copyright 2018 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 VARIATION_VARIATIONMANAGER_HXX_SEEN -#define VARIATION_VARIATIONMANAGER_HXX_SEEN +#pragma once #include "variation/IWeightProvider.hxx" #include "plugins/traits.hxx" #include "exception/exception.hxx" #include #include namespace fhicl { class ParameterSet; } namespace nuis { namespace variation { -class VariationManager { + +///\brief Provides event-weight responses for reweighting tools +class WeightManager { public: - typedef size_t VarProv_id_t; + typedef size_t WeightProv_id_t; private: struct NamedWeightProvider { NamedWeightProvider( std::string const &, plugins::plugin_traits::unique_ptr_t &&); + NamedWeightProvider(NamedWeightProvider const &) = delete; + NamedWeightProvider(NamedWeightProvider &&); std::string name; plugins::plugin_traits::unique_ptr_t handler; }; - std::vector VarProvs; + std::vector WeightEngines; - VariationManager(); + WeightManager(); - static VariationManager *_global_inst; + static WeightManager *_global_inst; public: + static WeightManager &Get(); - static VariationManager &Get(); - - paramId_t EnsureParameterHandled(fhicl::ParameterSet const &); - void SetParameterValue(paramId_t, double); - void GetParameterPull(paramId_t); - + WeightProv_id_t EnsureWeightProviderLoaded(fhicl::ParameterSet const &); double GetEventWeight(nuis::event::MinimalEvent const &); + + void ReconfigureWeightEngines(); }; } // namespace variation } // namespace nuis - -#endif