diff --git a/cmake/NEUTSetup.cmake b/cmake/NEUTSetup.cmake
index ac6fc40..5c8c2ce 100644
--- a/cmake/NEUTSetup.cmake
+++ b/cmake/NEUTSetup.cmake
@@ -1,252 +1,261 @@
# Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
include(cmake/parseConfigApp.cmake)
find_program(NEUTCONFIG NAMES neut-config)
LIST(APPEND EXTRA_CXX_FLAGS -DNEED_FILL_NEUT_COMMONS)
SET(HAVENEUTCONFIG FALSE)
# We are dealing with shiny NEUT
if(NOT "${NEUTCONFIG}" STREQUAL "NEUTCONFIG-NOTFOUND")
SET(HAVENEUTCONFIG TRUE)
cmessage(STATUS "Found neut-config, using it to determine configuration.")
else()
cmessage(STATUS "Failed to find neut-config, assuming older NEUT build.")
endif()
if(HAVENEUTCONFIG)
execute_process (COMMAND neut-config
--version OUTPUT_VARIABLE NEUT_VER
OUTPUT_STRIP_TRAILING_WHITESPACE)
execute_process (COMMAND neut-config
--incdir OUTPUT_VARIABLE NEUT_INCLUDE_DIRS
OUTPUT_STRIP_TRAILING_WHITESPACE)
execute_process (COMMAND neut-config
--libdir OUTPUT_VARIABLE NEUT_LINK_DIRS
OUTPUT_STRIP_TRAILING_WHITESPACE)
GETLIBDIRS(neut-config --cernflags CERN_LIB_DIR)
LIST(APPEND NEUT_LINK_DIRS ${CERN_LIB_DIR})
GETLIBS(neut-config --cernflags CERN_LIBS)
- if(USE_REWEIGHT)
+ if(NOT DEFINED USE_NEUT_REWEIGHT)
+ if(USE_REWEIGHT)
+ SET(USE_NEUT_REWEIGHT ON)
+ else()
+ SET(USE_NEUT_REWEIGHT OFF)
+ endif()
+ endif()
+
+ if(USE_REWEIGHT AND USE_NEUT_REWEIGHT)
execute_process (COMMAND neut-config
--rwlibflags OUTPUT_VARIABLE NEUT_RWLIBS
OUTPUT_STRIP_TRAILING_WHITESPACE)
GETLIBS(neut-config --rwlibflags NEUT_RWLIBS)
LIST(APPEND NEUT_LIBS ${NEUT_RWLIBS})
+ LIST(APPEND EXTRA_CXX_FLAGS ${NEUT_INCLUDE_DIRS} -D__USE_NEUT_REWEIGHT__)
else()
GETLIBS(neut-config --libflags NEUT_GENLIBS)
GETLIBS(neut-config --iolibflags NEUT_LIBS)
LIST(APPEND NEUT_LIBS ${NEUT_IOLIBS})
LIST(APPEND NEUT_LIBS ${NEUT_GENLIBS})
endif()
LIST(APPEND NEUT_LIBS ${CERN_LIBS};gfortran)
LIST(APPEND EXTRA_LIBS ${NEUT_LIBS})
string(REPLACE "." "" NEUT_VERSION ${NEUT_VER})
PrefixList(NEUT_INCLUDE_DIRS "-I" ${NEUT_INCLUDE_DIRS})
LIST(APPEND EXTRA_CXX_FLAGS ${NEUT_INCLUDE_DIRS} -D__NEUT_ENABLED__ -D__NEUT_VERSION__=${NEUT_VERSION})
LIST(APPEND EXTRA_LINK_DIRS ${NEUT_LINK_DIRS})
include(CheckCXXCompilerFlag)
CHECK_CXX_COMPILER_FLAG(-Wl,--allow-multiple-definition COMPILER_SUPPORTS_ALLOW_MULTIPLE_DEFINITION)
IF(COMPILER_SUPPORTS_ALLOW_MULTIPLE_DEFINITION)
LIST(APPEND EXTRA_LINK_FLAGS -Wl,--allow-multiple-definition)
ENDIF()
CHECK_CXX_COMPILER_FLAG(-no-pie COMPILER_SUPPORTS_NO_PIE)
CHECK_CXX_COMPILER_FLAG(-fno-pie COMPILER_SUPPORTS_FNO_PIE)
CHECK_CXX_COMPILER_FLAG(-fno-PIE COMPILER_SUPPORTS_FNO_PIE_CAP)
if(COMPILER_SUPPORTS_NO_PIE)
set(PIE_FLAGS "-no-pie")
elseif(COMPILER_SUPPORTS_FNO_PIE)
set(PIE_FLAGS "-fno-pie")
elseif(COMPILER_SUPPOERTS_FNO_PIE_CAP)
set(PIE_FLAGS "-fno-PIE")
else()
message(STATUS "The compiler ${CMAKE_CXX_COMPILER} has no C++11 support. Please use a different C++ compiler.")
endif()
LIST(APPEND EXTRA_EXE_FLAGS
${PIE_FLAGS})
cmessage(STATUS "NEUT")
cmessage(STATUS " Version : ${NEUT_VER}")
cmessage(STATUS " Flags : ${NEUT_CXX_FLAGS}")
cmessage(STATUS " Includes : ${NEUT_INCLUDE_DIRS}")
cmessage(STATUS " Link Dirs : ${NEUT_LINK_DIRS}")
cmessage(STATUS " Libs : ${NEUT_LIBS}")
cmessage(STATUS " Exe Flags : ${EXTRA_EXE_FLAGS}")
else() # Everything better be set up already
if(NEUT_ROOT STREQUAL "")
cmessage(FATAL_ERROR "Variable NEUT_ROOT is not defined. Please export environment variable NEUT_ROOT or configure with -DNEUT_ROOT=/path/to/NEUT. This must be set to point to a prebuilt NEUT instance.")
endif()
if(CERN STREQUAL "")
cmessage(FATAL_ERROR "Variable CERN is not defined. Please export environment variable CERN or configure with -DCERN=/path/to/CERNLIB. This must be set to point to a prebuilt CERNLIB instance.")
endif()
if(CERN_LEVEL STREQUAL "")
cmessage(FATAL_ERROR "Variable CERN_LEVEL is not defined. Please export environment variable CERN_LEVEL or configure with -DCERN_LEVEL=XXXX (likely to be 2005).")
endif()
if(${NEUT_VERSION} VERSION_LESS 5.4.0)
set(NEUT_LIB_DIR ${NEUT_ROOT}/lib/Linux_pc)
else()
set(NEUT_LIB_DIR ${NEUT_ROOT}/lib)
endif()
set(NEUT_CLASS ${NEUT_ROOT}/src/neutclass)
LIST(APPEND EXTRA_CXX_FLAGS -D__NEUT_ENABLED__ -DNEUT_VERSION=${NEUT_VERSION})
LIST(APPEND EXTRA_CXX_FLAGS
-I${NEUT_ROOT}/include
-I${NEUT_ROOT}/src/neutclass)
LIST(APPEND EXTRA_LINK_DIRS
${NEUT_LIB_DIR}
${CERN}/${CERN_LEVEL}/lib)
if(USE_REWEIGHT)
LIST(APPEND EXTRA_CXX_FLAGS
-I${NEUT_ROOT}/src/reweight)
LIST(APPEND EXTRA_LINK_DIRS
${NEUT_ROOT}/src/reweight)
endif()
if(${NEUT_VERSION} VERSION_EQUAL 5.4.2)
LIST(APPEND EXTRA_LIBS
-Wl,--as-needed)
if(USE_REWEIGHT)
LIST(APPEND EXTRA_LIBS
NReWeight)
endif()
LIST(APPEND EXTRA_LIBS
-Wl,--start-group
neutcore_5.4.2
nuccorspl_5.4.2 #typo in NEUT, may hopefully disappear
nuceff_5.4.2
partnuck_5.4.2
skmcsvc_5.4.2
tauola_5.4.2
HT2p2h_5.4.0
N1p1h_5.4.0
-Wl,--end-group
jetset74
pdflib804
mathlib
packlib
pawlib)
LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_COMMON_QEAV)
elseif(${NEUT_VERSION} VERSION_EQUAL 5.4.0)
LIST(APPEND EXTRA_LIBS
-Wl,--as-needed)
if(USE_REWEIGHT)
LIST(APPEND EXTRA_LIBS
NReWeight)
endif()
LIST(APPEND EXTRA_LIBS
-Wl,--start-group
neutcore_5.4.0
nuccorspl_5.4.0 #typo in NEUT, may hopefully disappear
nuceff_5.4.0
partnuck_5.4.0
skmcsvc_5.4.0
tauola_5.4.0
HT2p2h_5.4.0
N1p1h_5.4.0
specfunc_5.4.0
radcorr_5.4.0
gfortran
-Wl,--end-group
jetset74
pdflib804
mathlib
packlib
pawlib)
else()
LIST(APPEND EXTRA_LIBS
-Wl,--as-needed)
if(USE_REWEIGHT)
LIST(APPEND EXTRA_LIBS
NReWeight)
endif()
LIST(APPEND EXTRA_LIBS
-Wl,--start-group
neutcore
nuccorrspl
nuceff
partnuck
skmcsvc
tauola
-Wl,--end-group
jetset74
pdflib804
mathlib
packlib
pawlib)
endif()
set(NEUT_ROOT_LIBS)
LIST(APPEND NEUT_ROOT_LIBS
${NEUT_CLASS}/neutctrl.so
${NEUT_CLASS}/neutfsivert.so)
# Check for new versions of NEUT with NUCLEON FSI
if(EXISTS "${NEUT_CLASS}/neutnucfsistep.so")
set(NEUT_NUCFSI 1)
LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_NUCFSI_ENABLED)
LIST(APPEND NEUT_ROOT_LIBS
${NEUT_CLASS}/neutnucfsistep.so
${NEUT_CLASS}/neutnucfsivert.so
)
endif()
if(${NEUT_VERSION} VERSION_LESS 5.4.0)
LIST(APPEND NEUT_ROOT_LIBS
${NEUT_CLASS}/neutrootTreeSingleton.so)
endif()
LIST(APPEND NEUT_ROOT_LIBS
${NEUT_CLASS}/neutvtx.so
${NEUT_CLASS}/neutfsipart.so
${NEUT_CLASS}/neutpart.so
${NEUT_CLASS}/neutvect.so
)
foreach(OBJ ${NEUT_ROOT_LIBS})
LIST(APPEND EXTRA_SHAREDOBJS ${OBJ})
endforeach()
endif()
diff --git a/cmake/T2KSetup.cmake b/cmake/T2KSetup.cmake
index a4d3028..19c1582 100644
--- a/cmake/T2KSetup.cmake
+++ b/cmake/T2KSetup.cmake
@@ -1,160 +1,160 @@
# Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
include(cmake/parseConfigApp.cmake)
find_program(T2KRWCONFIG NAMES t2kreweight-config)
SET(HAVET2KRWCONFIG FALSE)
# We are dealing with shiny NEUT
if(NOT "${T2KRWCONFIG}" STREQUAL "T2KRWCONFIG-NOTFOUND")
SET(HAVET2KRWCONFIG TRUE)
- cmessage(STATUS "Found neut-config, using it to determine configuration.")
+ cmessage(STATUS "Found t2kreweight-config, using it to determine configuration.")
else()
- cmessage(STATUS "Failed to find neut-config, assuming older NEUT build.")
+ cmessage(STATUS "Failed to find t2kreweight-config, assuming older NEUT build.")
endif()
if(HAVET2KRWCONFIG)
if(NOT DEFINED CMAKE_CXX_STANDARD OR "${CMAKE_CXX_STANDARD} " STREQUAL " ")
SET(CMAKE_CXX_STANDARD 11)
endif()
LIST(APPEND EXTRA_CXX_FLAGS -DT2KRW_OA2021_INTERFACE -D__T2KREW_ENABLED__)
GETLIBDIRS(t2kreweight-config --linkflags T2KRW_LINK_DIRS)
GETINCDIRS(t2kreweight-config --cflags T2KRW_INC_DIRS)
GETLIBS(t2kreweight-config --linkflags T2KRW_LIBS)
cmessage(STATUS "T2KReWeight:")
cmessage(STATUS " LINK DIRS: ${T2KRW_LINK_DIRS}")
cmessage(STATUS " INC DIRS: ${T2KRW_INC_DIRS}")
cmessage(STATUS " LIBS: ${T2KRW_LIBS}")
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${T2KRW_INC_DIRS})
LIST(APPEND EXTRA_LINK_DIRS ${T2KRW_LINK_DIRS})
LIST(APPEND EXTRA_LIBS ${T2KRW_LIBS})
execute_process (COMMAND t2kreweight-config
--has-feature NIWG RESULT_VARIABLE T2KRW_HAS_NIWG)
if("${T2KRW_HAS_NIWG} " STREQUAL "0 ")
GETLIBDIRS(t2kreweight-config --niwgflags NIWG_LINK_DIRS)
GETINCDIRS(t2kreweight-config --niwgflags NIWG_INC_DIRS)
GETLIBS(t2kreweight-config --niwgflags NIWG_LIBS)
cmessage(STATUS "NIWG:")
cmessage(STATUS " LINK DIRS: ${NIWG_LINK_DIRS}")
cmessage(STATUS " INC DIRS: ${NIWG_INC_DIRS}")
cmessage(STATUS " LIBS: ${NIWG_LIBS}")
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${NIWG_INC_DIRS})
LIST(APPEND EXTRA_LINK_DIRS ${NIWG_LINK_DIRS})
LIST(APPEND EXTRA_LIBS ${NIWG_LIBS})
endif()
execute_process (COMMAND t2kreweight-config
--has-feature NEUT RESULT_VARIABLE T2KRW_HAS_NEUT)
if("${T2KRW_HAS_NEUT} " STREQUAL "0 ")
LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_EVENT_ENABLED)
set(USE_NEUT_EVENT TRUE)
set(USE_NEUT_EVENT TRUE CACHE BOOL "Whether to enable NEUT (event i/o) support. Requires external libraries. " FORCE)
GETLIBDIRS(t2kreweight-config --neutflags NEUT_LINK_DIRS)
GETINCDIRS(t2kreweight-config --neutflags NEUT_INC_DIRS)
GETLIBS(t2kreweight-config --neutflags NEUT_LIBS)
cmessage(STATUS "NEUT:")
cmessage(STATUS " LINK DIRS: ${NEUT_LINK_DIRS}")
cmessage(STATUS " INC DIRS: ${NEUT_INC_DIRS}")
cmessage(STATUS " LIBS: ${NEUT_LIBS}")
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${NEUT_INC_DIRS})
LIST(APPEND EXTRA_LINK_DIRS ${NEUT_LINK_DIRS})
LIST(APPEND EXTRA_LIBS ${NEUT_LIBS})
endif()
execute_process (COMMAND t2kreweight-config
--has-feature GEANT RESULT_VARIABLE T2KRW_HAS_GEANT)
if("${T2KRW_HAS_GEANT} " STREQUAL "0 ")
GETLIBDIRS(t2kreweight-config --geantflags GEANT_LINK_DIRS)
GETINCDIRS(t2kreweight-config --geantflags GEANT_INC_DIRS)
GETLIBS(t2kreweight-config --geantflags GEANT_LIBS)
cmessage(STATUS "GEANT:")
cmessage(STATUS " LINK DIRS: ${GEANT_LINK_DIRS}")
cmessage(STATUS " INC DIRS: ${GEANT_INC_DIRS}")
cmessage(STATUS " LIBS: ${GEANT_LIBS}")
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${GEANT_INC_DIRS})
LIST(APPEND EXTRA_LINK_DIRS ${GEANT_LINK_DIRS})
LIST(APPEND EXTRA_LIBS ${GEANT_LIBS})
endif()
execute_process (COMMAND t2kreweight-config
--has-feature ND280 RESULT_VARIABLE T2KRW_HAS_ND280)
if("${T2KRW_HAS_ND280} " STREQUAL "0 ")
GETLIBDIRS(t2kreweight-config --nd280flags ND280_LINK_DIRS)
GETINCDIRS(t2kreweight-config --nd280flags ND280_INC_DIRS)
GETLIBS(t2kreweight-config --nd280flags ND280_LIBS)
cmessage(STATUS "ND280:")
cmessage(STATUS " LINK DIRS: ${ND280_LINK_DIRS}")
cmessage(STATUS " INC DIRS: ${ND280_INC_DIRS}")
cmessage(STATUS " LIBS: ${ND280_LIBS}")
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${ND280_INC_DIRS})
LIST(APPEND EXTRA_LINK_DIRS ${ND280_LINK_DIRS})
LIST(APPEND EXTRA_LIBS ${ND280_LIBS})
endif()
else()
if(T2KREWEIGHT STREQUAL "")
cmessage(FATAL_ERROR "Variable T2KREWEIGHT is not defined. Either configure with -DT2KREWEIGHT or \"\$ export T2KREWEIGHT=/path/to/T2KReWeight\". This must be set to point to a prebuilt T2KReWeight instance.")
endif()
LIST(APPEND EXTRA_CXX_FLAGS -D__T2KREW_ENABLED__ )
# First check the OAANALYSIS libs (need to grab some headers for T2KReWeight linking if compiled with oaAnalysisReader)
IF(NOT $ENV{OAANALYSISREADERROOT} STREQUAL "")
cmessage(STATUS "Found OAANALYSISREADERROOT $ENV{OAANALYSISREADERROOT}, appending...")
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES $ENV{OAANALYSISREADERROOT}/$ENV{OAANALYSISREADERCONFIG})
LIST(APPEND EXTRA_LINK_DIRS $ENV{OAANALYSISREADERROOT}/$ENV{OAANALYSISREADERCONFIG})
LIST(APPEND EXTRA_LIBS oaAnalysisReader)
# Don't have to append this; should be appeneded in ${T2KReWeight/src/T2KBuild.h}
#LIST(APPEND EXTRA_CXX_FLAGS -D__T2KRW_OAANALYSIS_ENABLED__)
endif()
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES ${T2KREWEIGHT}/src/)
LIST(APPEND EXTRA_LINK_DIRS ${T2KREWEIGHT}/lib)
LIST(APPEND EXTRA_LIBS T2KReWeight)
endif()
\ No newline at end of file
diff --git a/src/FCN/SampleList.cxx b/src/FCN/SampleList.cxx
index 12d065e..6b28a5b 100644
--- a/src/FCN/SampleList.cxx
+++ b/src/FCN/SampleList.cxx
@@ -1,1579 +1,1576 @@
#include "SampleList.h"
#ifndef __NO_ANL__
#include "ANL_CCQE_Evt_1DQ2_nu.h"
#include "ANL_CCQE_XSec_1DEnu_nu.h"
// ANL CC1ppip
#include "ANL_CC1ppip_Evt_1DQ2_nu.h"
#include "ANL_CC1ppip_Evt_1DcosmuStar_nu.h"
#include "ANL_CC1ppip_Evt_1DcosthAdler_nu.h"
#include "ANL_CC1ppip_Evt_1Dphi_nu.h"
#include "ANL_CC1ppip_Evt_1Dppi_nu.h"
#include "ANL_CC1ppip_Evt_1Dthpr_nu.h"
#include "ANL_CC1ppip_XSec_1DEnu_nu.h"
#include "ANL_CC1ppip_XSec_1DQ2_nu.h"
#include "ANL_CC1ppip_Evt_1DWNpi_nu.h"
#include "ANL_CC1ppip_Evt_1DWNmu_nu.h"
#include "ANL_CC1ppip_Evt_1DWmupi_nu.h"
// ANL CC1npip
#include "ANL_CC1npip_Evt_1DQ2_nu.h"
#include "ANL_CC1npip_Evt_1DcosmuStar_nu.h"
#include "ANL_CC1npip_Evt_1Dppi_nu.h"
#include "ANL_CC1npip_XSec_1DEnu_nu.h"
#include "ANL_CC1npip_Evt_1DWNpi_nu.h"
#include "ANL_CC1npip_Evt_1DWNmu_nu.h"
#include "ANL_CC1npip_Evt_1DWmupi_nu.h"
// ANL CC1pi0
#include "ANL_CC1pi0_Evt_1DQ2_nu.h"
#include "ANL_CC1pi0_Evt_1DcosmuStar_nu.h"
#include "ANL_CC1pi0_XSec_1DEnu_nu.h"
#include "ANL_CC1pi0_Evt_1DWNpi_nu.h"
#include "ANL_CC1pi0_Evt_1DWNmu_nu.h"
#include "ANL_CC1pi0_Evt_1DWmupi_nu.h"
// ANL NC1npip (mm, exotic!)
#include "ANL_NC1npip_Evt_1Dppi_nu.h"
// ANL NC1ppim (mm, exotic!)
#include "ANL_NC1ppim_Evt_1DcosmuStar_nu.h"
#include "ANL_NC1ppim_XSec_1DEnu_nu.h"
// ANL CC2pi 1pim1pip (mm, even more exotic!)
#include "ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu.h"
#include "ANL_CC2pi_1pim1pip_Evt_1Dppim_nu.h"
#include "ANL_CC2pi_1pim1pip_Evt_1Dppip_nu.h"
#include "ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu.h"
#include "ANL_CC2pi_1pim1pip_XSec_1DEnu_nu.h"
// ANL CC2pi 1pip1pip (mm, even more exotic!)
#include "ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu.h"
#include "ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu.h"
#include "ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu.h"
#include "ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu.h"
#include "ANL_CC2pi_1pip1pip_XSec_1DEnu_nu.h"
// ANL CC2pi 1pip1pi0 (mm, even more exotic!)
#include "ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu.h"
#include "ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu.h"
#include "ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu.h"
#include "ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu.h"
#include "ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu.h"
#endif
#ifndef __NO_ArgoNeuT__
// ArgoNeuT CC1Pi
#include "ArgoNeuT_CC1Pi_XSec_1Dpmu_antinu.h"
#include "ArgoNeuT_CC1Pi_XSec_1Dpmu_nu.h"
#include "ArgoNeuT_CC1Pi_XSec_1Dthetamu_antinu.h"
#include "ArgoNeuT_CC1Pi_XSec_1Dthetamu_nu.h"
#include "ArgoNeuT_CC1Pi_XSec_1Dthetamupi_antinu.h"
#include "ArgoNeuT_CC1Pi_XSec_1Dthetamupi_nu.h"
#include "ArgoNeuT_CC1Pi_XSec_1Dthetapi_antinu.h"
#include "ArgoNeuT_CC1Pi_XSec_1Dthetapi_nu.h"
// ArgoNeuT CC-inclusive
#include "ArgoNeuT_CCInc_XSec_1Dpmu_antinu.h"
#include "ArgoNeuT_CCInc_XSec_1Dpmu_nu.h"
#include "ArgoNeuT_CCInc_XSec_1Dthetamu_antinu.h"
#include "ArgoNeuT_CCInc_XSec_1Dthetamu_nu.h"
#endif
#ifndef __NO_BNL__
// BNL CCQE
#include "BNL_CCQE_Evt_1DQ2_nu.h"
#include "BNL_CCQE_XSec_1DEnu_nu.h"
// BNL CC1ppip
#include "BNL_CC1ppip_Evt_1DQ2_nu.h"
#include "BNL_CC1ppip_Evt_1DcosthAdler_nu.h"
#include "BNL_CC1ppip_Evt_1Dphi_nu.h"
#include "BNL_CC1ppip_XSec_1DEnu_nu.h"
#include "BNL_CC1ppip_Evt_1DWNpi_nu.h"
#include "BNL_CC1ppip_Evt_1DWNmu_nu.h"
#include "BNL_CC1ppip_Evt_1DWmupi_nu.h"
// BNL CC1npip
#include "BNL_CC1npip_Evt_1DQ2_nu.h"
#include "BNL_CC1npip_XSec_1DEnu_nu.h"
#include "BNL_CC1npip_Evt_1DWNpi_nu.h"
#include "BNL_CC1npip_Evt_1DWNmu_nu.h"
#include "BNL_CC1npip_Evt_1DWmupi_nu.h"
// BNL CC1pi0
#include "BNL_CC1pi0_Evt_1DQ2_nu.h"
#include "BNL_CC1pi0_XSec_1DEnu_nu.h"
#include "BNL_CC1pi0_Evt_1DWNpi_nu.h"
#include "BNL_CC1pi0_Evt_1DWNmu_nu.h"
#include "BNL_CC1pi0_Evt_1DWmupi_nu.h"
// BNL multipi
#include "BNL_CC2pi_1pim1pip_XSec_1DEnu_nu.cxx"
#include "BNL_CC3pi_1pim2pip_XSec_1DEnu_nu.cxx"
#include "BNL_CC4pi_2pim2pip_XSec_1DEnu_nu.cxx"
#include "BNL_CC2pi_1pim1pip_Evt_1DWpippim_nu.cxx"
#include "BNL_CC2pi_1pim1pip_Evt_1DWpippr_nu.cxx"
#endif
#ifndef __NO_FNAL__
// FNAL CCQE
#include "FNAL_CCQE_Evt_1DQ2_nu.h"
// FNAL CC1ppip
#include "FNAL_CC1ppip_Evt_1DQ2_nu.h"
#include "FNAL_CC1ppip_XSec_1DEnu_nu.h"
#include "FNAL_CC1ppip_XSec_1DQ2_nu.h"
// FNAL CC1ppim
#include "FNAL_CC1ppim_XSec_1DEnu_antinu.h"
#endif
#ifndef __NO_BEBC__
// BEBC CCQE
#include "BEBC_CCQE_XSec_1DQ2_nu.h"
// BEBC CC1ppip
#include "BEBC_CC1ppip_XSec_1DEnu_nu.h"
#include "BEBC_CC1ppip_XSec_1DQ2_nu.h"
// BEBC CC1npip
#include "BEBC_CC1npip_XSec_1DEnu_nu.h"
#include "BEBC_CC1npip_XSec_1DQ2_nu.h"
// BEBC CC1pi0
#include "BEBC_CC1pi0_XSec_1DEnu_nu.h"
#include "BEBC_CC1pi0_XSec_1DQ2_nu.h"
// BEBC CC1npim
#include "BEBC_CC1npim_XSec_1DEnu_antinu.h"
#include "BEBC_CC1npim_XSec_1DQ2_antinu.h"
// BEBC CC1ppim
#include "BEBC_CC1ppim_XSec_1DEnu_antinu.h"
#include "BEBC_CC1ppim_XSec_1DQ2_antinu.h"
#endif
#ifndef __NO_GGM__
// GGM CC1ppip
#include "GGM_CC1ppip_Evt_1DQ2_nu.h"
#include "GGM_CC1ppip_XSec_1DEnu_nu.h"
#endif
#ifndef __NO_MiniBooNE__
// MiniBooNE CCQE
#include "MiniBooNE_CCQE_XSec_1DEnu_nu.h"
#include "MiniBooNE_CCQE_XSec_1DQ2_antinu.h"
#include "MiniBooNE_CCQE_XSec_1DQ2_nu.h"
#include "MiniBooNE_CCQE_XSec_2DTcos_antinu.h"
#include "MiniBooNE_CCQE_XSec_2DTcos_nu.h"
// MiniBooNE CC1pi+ 1D
#include "MiniBooNE_CC1pip_XSec_1DEnu_nu.h"
#include "MiniBooNE_CC1pip_XSec_1DQ2_nu.h"
#include "MiniBooNE_CC1pip_XSec_1DTpi_nu.h"
#include "MiniBooNE_CC1pip_XSec_1DTu_nu.h"
// MiniBooNE CC1pi+ 2D
#include "MiniBooNE_CC1pip_XSec_2DQ2Enu_nu.h"
#include "MiniBooNE_CC1pip_XSec_2DTpiCospi_nu.h"
#include "MiniBooNE_CC1pip_XSec_2DTpiEnu_nu.h"
#include "MiniBooNE_CC1pip_XSec_2DTuCosmu_nu.h"
#include "MiniBooNE_CC1pip_XSec_2DTuEnu_nu.h"
// MiniBooNE CC1pi0
#include "MiniBooNE_CC1pi0_XSec_1DEnu_nu.h"
#include "MiniBooNE_CC1pi0_XSec_1DQ2_nu.h"
#include "MiniBooNE_CC1pi0_XSec_1DTu_nu.h"
#include "MiniBooNE_CC1pi0_XSec_1Dcosmu_nu.h"
#include "MiniBooNE_CC1pi0_XSec_1Dcospi0_nu.h"
#include "MiniBooNE_CC1pi0_XSec_1Dppi0_nu.h"
#include "MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu.h"
#include "MiniBooNE_NC1pi0_XSec_1Dcospi0_nu.h"
#include "MiniBooNE_NC1pi0_XSec_1Dppi0_antinu.h"
#include "MiniBooNE_NC1pi0_XSec_1Dppi0_nu.h"
-// MiniBooNE NC1pi0
-//#include "MiniBooNE_NCpi0_XSec_1Dppi0_nu.h"
-
// MiniBooNE NCEL
#include "MiniBooNE_NCEL_XSec_Treco_nu.h"
#endif
#ifndef __NO_MicroBooNE__
#include "MicroBooNE_CCInc_XSec_2DPcos_nu.h"
#include "MicroBooNE_CC1MuNp_XSec_1D_nu.h"
#endif
#ifndef __NO_MINERvA__
// MINERvA CCQE
#include "MINERvA_CCQE_XSec_1DQ2_antinu.h"
#include "MINERvA_CCQE_XSec_1DQ2_joint.h"
#include "MINERvA_CCQE_XSec_1DQ2_nu.h"
// MINERvA CC0pi
#include "MINERvA_CC0pi_XSec_1DEe_nue.h"
#include "MINERvA_CC0pi_XSec_1DQ2_nu_proton.h"
#include "MINERvA_CC0pi_XSec_1DQ2_nue.h"
#include "MINERvA_CC0pi_XSec_1DThetae_nue.h"
// 2018 MINERvA CC0pi STV
#include "MINERvA_CC0pinp_STV_XSec_1D_nu.h"
// 2018 MINERvA CC0pi 2D
#include "MINERvA_CC0pi_XSec_1D_2018_nu.h"
#include "MINERvA_CC0pi_XSec_2D_nu.h"
// #include "MINERvA_CC0pi_XSec_3DptpzTp_nu.h"
// 2018 MINERvA CC0pi 2D antinu
#include "MINERvA_CC0pi_XSec_2D_antinu.h"
// MINERvA CC1pi+
#include "MINERvA_CC1pip_XSec_1DTpi_20deg_nu.h"
#include "MINERvA_CC1pip_XSec_1DTpi_nu.h"
#include "MINERvA_CC1pip_XSec_1Dth_20deg_nu.h"
#include "MINERvA_CC1pip_XSec_1Dth_nu.h"
// 2017 data update
#include "MINERvA_CC1pip_XSec_1D_2017Update.h"
// MINERvA CC1pi-
#include "MINERvA_CC1pim_XSec_1DEnu_antinu.h"
#include "MINERvA_CC1pim_XSec_1DQ2_antinu.h"
#include "MINERvA_CC1pim_XSec_1DTpi_antinu.h"
#include "MINERvA_CC1pim_XSec_1Dpmu_antinu.h"
#include "MINERvA_CC1pim_XSec_1Dth_antinu.h"
#include "MINERvA_CC1pim_XSec_1Dthmu_antinu.h"
// MINERvA CCNpi+
#include "MINERvA_CCNpip_XSec_1DEnu_nu.h"
#include "MINERvA_CCNpip_XSec_1DQ2_nu.h"
#include "MINERvA_CCNpip_XSec_1DTpi_nu.h"
#include "MINERvA_CCNpip_XSec_1Dpmu_nu.h"
#include "MINERvA_CCNpip_XSec_1Dth_nu.h"
#include "MINERvA_CCNpip_XSec_1Dthmu_nu.h"
// MINERvA CC1pi0
#include "MINERvA_CC1pi0_XSec_1DEnu_antinu.h"
#include "MINERvA_CC1pi0_XSec_1DQ2_antinu.h"
#include "MINERvA_CC1pi0_XSec_1DTpi0_antinu.h"
#include "MINERvA_CC1pi0_XSec_1Dpmu_antinu.h"
#include "MINERvA_CC1pi0_XSec_1Dppi0_antinu.h"
#include "MINERvA_CC1pi0_XSec_1Dth_antinu.h"
#include "MINERvA_CC1pi0_XSec_1Dthmu_antinu.h"
// MINERvA CC1pi0 neutrino
#include "MINERvA_CC1pi0_XSec_1D_nu.h"
// MINERvA CCINC
#include "MINERvA_CCinc_XSec_1DEnu_ratio.h"
#include "MINERvA_CCinc_XSec_1Dx_ratio.h"
#include "MINERvA_CCinc_XSec_2DEavq3_nu.h"
// MINERvA CCDIS
#include "MINERvA_CCDIS_XSec_1DEnu_ratio.h"
#include "MINERvA_CCDIS_XSec_1Dx_ratio.h"
// MINERvA CCCOH pion
#include "MINERvA_CCCOHPI_XSec_1DEnu_antinu.h"
#include "MINERvA_CCCOHPI_XSec_1DEpi_antinu.h"
#include "MINERvA_CCCOHPI_XSec_1DQ2_antinu.h"
#include "MINERvA_CCCOHPI_XSec_1DEpi_nu.h"
#include "MINERvA_CCCOHPI_XSec_1DQ2_nu.h"
#include "MINERvA_CCCOHPI_XSec_1Dth_nu.h"
#include "MINERvA_CCCOHPI_XSec_joint.h"
#include "MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.h"
#include "MINERvA_CC0pi_XSec_1DQ2_Tgt_nu.h"
#endif
#ifndef __NO_T2K__
// T2K CC0pi 2016
#include "T2K_CC0pi_XSec_2DPcos_nu_I.h"
#include "T2K_CC0pi_XSec_2DPcos_nu_II.h"
// T2K CC0pi 2020 arXiv:1908.10249
#include "T2K_CC0pi_XSec_H2O_2DPcos_anu.h"
// T2K CC0pi 2020 arXiv:2004.05434
#include "T2K_NuMu_CC0pi_OC_XSec_2DPcos.h"
#include "T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint.h"
// T2K CC0pi 2020 arXiv:2002.09323
#include "T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos.h"
#include "T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint.h"
// T2K CC-inclusive with full acceptance 2018
#include "T2K_CCinc_XSec_2DPcos_nu_nonuniform.h"
// T2K nue CC-inclusive 2019
#include "T2K_nueCCinc_XSec_1Dpe.h"
#include "T2K_nueCCinc_XSec_1Dthe.h"
#include "T2K_nueCCinc_XSec_1Dpe_joint.h"
#include "T2K_nueCCinc_XSec_1Dthe_joint.h"
#include "T2K_nueCCinc_XSec_joint.h"
// T2K STV CC0pi 2018
#include "T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p_Np.h"
#include "T2K_CC0pinp_STV_XSec_1Ddat_nu.h"
#include "T2K_CC0pinp_STV_XSec_1Ddphit_nu.h"
#include "T2K_CC0pinp_STV_XSec_1Ddpt_nu.h"
#include "T2K_CC0pinp_ifk_XSec_3Dinfa_nu.h"
#include "T2K_CC0pinp_ifk_XSec_3Dinfip_nu.h"
#include "T2K_CC0pinp_ifk_XSec_3Dinfp_nu.h"
// T2K CC1pi+ on CH
#include "T2K_CC1pip_CH_XSec_1DAdlerPhi_nu.h"
#include "T2K_CC1pip_CH_XSec_1DCosThAdler_nu.h"
#include "T2K_CC1pip_CH_XSec_1DQ2_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dppi_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dthmupi_nu.h"
#include "T2K_CC1pip_CH_XSec_1Dthpi_nu.h"
#include "T2K_CC1pip_CH_XSec_2Dpmucosmu_nu.h"
// T2K CC1pi+ on H2O
#include "T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.h"
#include "T2K_CC1pip_H2O_XSec_1DEnuMB_nu.h"
#include "T2K_CC1pip_H2O_XSec_1Dcosmu_nu.h"
#include "T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.h"
#include "T2K_CC1pip_H2O_XSec_1Dcospi_nu.h"
#include "T2K_CC1pip_H2O_XSec_1Dpmu_nu.h"
#include "T2K_CC1pip_H2O_XSec_1Dppi_nu.h"
// add header here
#endif
#ifndef __NO_SciBooNE__
// SciBooNE COH studies
#include "SciBooNE_CCCOH_1TRK_1DQ2_nu.h"
#include "SciBooNE_CCCOH_1TRK_1Dpmu_nu.h"
#include "SciBooNE_CCCOH_1TRK_1Dthetamu_nu.h"
#include "SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu.h"
#include "SciBooNE_CCCOH_MuPiNoVA_1Dpmu_nu.h"
#include "SciBooNE_CCCOH_MuPiNoVA_1Dthetamu_nu.h"
#include "SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu.h"
#include "SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu.h"
#include "SciBooNE_CCCOH_MuPiVA_1DQ2_nu.h"
#include "SciBooNE_CCCOH_MuPiVA_1Dpmu_nu.h"
#include "SciBooNE_CCCOH_MuPiVA_1Dthetamu_nu.h"
#include "SciBooNE_CCCOH_MuPr_1DQ2_nu.h"
#include "SciBooNE_CCCOH_MuPr_1Dpmu_nu.h"
#include "SciBooNE_CCCOH_MuPr_1Dthetamu_nu.h"
#include "SciBooNE_CCCOH_STOPFINAL_1DQ2_nu.h"
#include "SciBooNE_CCCOH_STOP_NTrks_nu.h"
#include "SciBooNE_CCInc_XSec_1DEnu_nu.h"
#endif
#ifndef __NO_K2K__
// K2K NC1pi0
#include "K2K_NC1pi0_Evt_1Dppi0_nu.h"
#endif
// MC Studies
#include "ExpMultDist_CCQE_XSec_1DVar_FakeStudy.h"
#include "ExpMultDist_CCQE_XSec_2DVar_FakeStudy.h"
#include "MCStudy_CCQEHistograms.h"
#include "GenericFlux_Tester.h"
#include "GenericFlux_Vectors.h"
#include "ElectronFlux_FlatTree.h"
#include "ElectronScattering_DurhamData.h"
#include "MCStudy_KaonPreSelection.h"
#include "MCStudy_MuonValidation.h"
#include "OfficialNIWGPlots.h"
#include "T2K2017_FakeData.h"
#include "SigmaEnuHists.h"
#include "Simple_Osc.h"
#include "Smear_SVDUnfold_Propagation_Osc.h"
#include "FitWeight.h"
#include "NuisConfig.h"
#include "NuisKey.h"
#ifdef __USE_DYNSAMPLES__
#include "TRegexp.h"
#include
// linux
#include
DynamicSampleFactory::DynamicSampleFactory() : NSamples(0), NManifests(0) {
LoadPlugins();
NUIS_LOG(FIT, "Loaded " << NSamples << " from " << NManifests
<< " shared object libraries.");
}
DynamicSampleFactory *DynamicSampleFactory::glblDSF = NULL;
DynamicSampleFactory::PluginManifest::~PluginManifest() {
for (size_t i_it = 0; i_it < Instances.size(); ++i_it) {
(*(DSF_DestroySample))(Instances[i_it]);
}
}
std::string EnsureTrailingSlash(std::string const &inp) {
if (!inp.length()) {
return "/";
}
if (inp[inp.length() - 1] == '/') {
return inp;
}
return inp + "/";
}
void DynamicSampleFactory::LoadPlugins() {
std::vector SearchDirectories;
if (Config::HasPar("dynamic_sample.path")) {
SearchDirectories =
GeneralUtils::ParseToStr(Config::GetParS("dynamic_sample.path"), ":");
}
char const *envPath = getenv("NUISANCE_DS_PATH");
if (envPath) {
std::vector envPaths = GeneralUtils::ParseToStr(envPath, ":");
for (size_t ep_it = 0; ep_it < envPaths.size(); ++ep_it) {
SearchDirectories.push_back(envPaths[ep_it]);
}
}
if (!SearchDirectories.size()) {
char const *pwdPath = getenv("PWD");
if (pwdPath) {
SearchDirectories.push_back(pwdPath);
}
}
for (size_t sp_it = 0; sp_it < SearchDirectories.size(); ++sp_it) {
std::string dirpath = EnsureTrailingSlash(SearchDirectories[sp_it]);
NUIS_LOG(FIT, "Searching for dynamic sample manifests in: " << dirpath);
Ssiz_t len = 0;
DIR *dir;
struct dirent *ent;
dir = opendir(dirpath.c_str());
if (dir != NULL) {
TRegexp matchExp("*.so", true);
while ((ent = readdir(dir)) != NULL) {
if (matchExp.Index(TString(ent->d_name), &len) != Ssiz_t(-1)) {
NUIS_LOG(FIT, "\tFound shared object: "
<< ent->d_name
<< " checking for relevant methods...");
void *dlobj =
dlopen((dirpath + ent->d_name).c_str(), RTLD_NOW | RTLD_GLOBAL);
char const *dlerr_cstr = dlerror();
std::string dlerr;
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr.length()) {
NUIS_ERR(WRN, "\tDL Load Error: " << dlerr);
continue;
}
PluginManifest plgManif;
plgManif.dllib = dlobj;
plgManif.soloc = (dirpath + ent->d_name);
plgManif.DSF_NSamples =
reinterpret_cast(dlsym(dlobj, "DSF_NSamples"));
dlerr = "";
dlerr_cstr = dlerror();
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr.length()) {
NUIS_ERR(WRN, "\tFailed to load symbol \"DSF_NSamples\" from "
<< (dirpath + ent->d_name) << ": " << dlerr);
dlclose(dlobj);
continue;
}
plgManif.DSF_GetSampleName = reinterpret_cast(
dlsym(dlobj, "DSF_GetSampleName"));
dlerr = "";
dlerr_cstr = dlerror();
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr.length()) {
NUIS_ERR(WRN, "\tFailed to load symbol \"DSF_GetSampleName\" from "
<< (dirpath + ent->d_name) << ": " << dlerr);
dlclose(dlobj);
continue;
}
plgManif.DSF_GetSample = reinterpret_cast(
dlsym(dlobj, "DSF_GetSample"));
dlerr = "";
dlerr_cstr = dlerror();
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr.length()) {
NUIS_ERR(WRN, "\tFailed to load symbol \"DSF_GetSample\" from "
<< (dirpath + ent->d_name) << ": " << dlerr);
dlclose(dlobj);
continue;
}
plgManif.DSF_DestroySample = reinterpret_cast(
dlsym(dlobj, "DSF_DestroySample"));
dlerr = "";
dlerr_cstr = dlerror();
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr.length()) {
NUIS_ERR(WRN, "Failed to load symbol \"DSF_DestroySample\" from "
<< (dirpath + ent->d_name) << ": " << dlerr);
dlclose(dlobj);
continue;
}
plgManif.NSamples = (*(plgManif.DSF_NSamples))();
NUIS_LOG(FIT, "\tSuccessfully loaded dynamic sample manifest: "
<< plgManif.soloc << ". Contains "
<< plgManif.NSamples << " samples.");
for (size_t smp_it = 0; smp_it < plgManif.NSamples; ++smp_it) {
char const *smp_name = (*(plgManif.DSF_GetSampleName))(smp_it);
if (!smp_name) {
NUIS_ABORT("Could not load sample "
<< smp_it << " / " << plgManif.NSamples << " from "
<< plgManif.soloc);
}
if (Samples.count(smp_name)) {
NUIS_ERR(WRN, "Already loaded a sample named: \""
<< smp_name
<< "\". cannot load duplciates. This "
"sample will be skipped.");
continue;
}
plgManif.SamplesProvided.push_back(smp_name);
Samples[smp_name] = std::make_pair(plgManif.soloc, smp_it);
NUIS_LOG(FIT, "\t\t" << smp_name);
}
if (plgManif.SamplesProvided.size()) {
Manifests[plgManif.soloc] = plgManif;
NSamples += plgManif.SamplesProvided.size();
NManifests++;
} else {
dlclose(dlobj);
}
}
}
closedir(dir);
} else {
NUIS_ERR(WRN, "Tried to open non-existant directory.");
}
}
}
DynamicSampleFactory &DynamicSampleFactory::Get() {
if (!glblDSF) {
glblDSF = new DynamicSampleFactory();
}
return *glblDSF;
}
void DynamicSampleFactory::Print() {
std::map > ManifestSamples;
for (std::map >::iterator smp_it =
Samples.begin();
smp_it != Samples.end(); ++smp_it) {
if (!ManifestSamples.count(smp_it->second.first)) {
ManifestSamples[smp_it->second.first] = std::vector();
}
ManifestSamples[smp_it->second.first].push_back(smp_it->first);
}
NUIS_LOG(FIT, "Dynamic sample manifest: ");
for (std::map >::iterator m_it =
ManifestSamples.begin();
m_it != ManifestSamples.end(); ++m_it) {
NUIS_LOG(FIT, "\tLibrary " << m_it->first << " contains: ");
for (size_t s_it = 0; s_it < m_it->second.size(); ++s_it) {
NUIS_LOG(FIT, "\t\t" << m_it->second[s_it]);
}
}
}
bool DynamicSampleFactory::HasSample(std::string const &name) {
return Samples.count(name);
}
bool DynamicSampleFactory::HasSample(nuiskey &samplekey) {
return HasSample(samplekey.GetS("name"));
}
MeasurementBase *DynamicSampleFactory::CreateSample(nuiskey &samplekey) {
if (!HasSample(samplekey)) {
NUIS_ERR(WRN, "Asked to load unknown sample: \"" << samplekey.GetS("name")
<< "\".");
return NULL;
}
std::pair sample = Samples[samplekey.GetS("name")];
NUIS_LOG(SAM,
"\tLoading sample " << sample.second << " from " << sample.first);
return (*(Manifests[sample.first].DSF_GetSample))(sample.second, &samplekey);
}
DynamicSampleFactory::~DynamicSampleFactory() { Manifests.clear(); }
#endif
//! Functions to make it easier for samples to be created and handled.
namespace SampleUtils {
//! Create a given sample given its name, file, type, fakdata(fkdt) file and the
//! current rw engine and push it back into the list fChain.
MeasurementBase *CreateSample(std::string name, std::string file,
std::string type, std::string fkdt,
FitWeight *rw) {
nuiskey samplekey = Config::CreateKey("sample");
samplekey.Set("name", name);
samplekey.Set("input", file);
samplekey.Set("type", type);
return CreateSample(samplekey);
}
MeasurementBase *CreateSample(nuiskey samplekey) {
#ifdef __USE_DYNSAMPLES__
if (DynamicSampleFactory::Get().HasSample(samplekey)) {
NUIS_LOG(SAM, "Instantiating dynamic sample...");
MeasurementBase *ds = DynamicSampleFactory::Get().CreateSample(samplekey);
if (ds) {
NUIS_LOG(SAM, "Done.");
return ds;
}
NUIS_ABORT("Failed to instantiate dynamic sample.");
}
#endif
FitWeight *rw = FitBase::GetRW();
std::string name = samplekey.GetS("name");
std::string file = samplekey.GetS("input");
std::string type = samplekey.GetS("type");
std::string fkdt = "";
/*
ANL CCQE Samples
*/
#ifndef __NO_ANL__
if (!name.compare("ANL_CCQE_XSec_1DEnu_nu") ||
!name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD26") ||
!name.compare("ANL_CCQE_XSec_1DEnu_nu_PRL31") ||
!name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD16")) {
return (new ANL_CCQE_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CCQE_Evt_1DQ2_nu") ||
!name.compare("ANL_CCQE_Evt_1DQ2_nu_PRL31") ||
!name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD26") ||
!name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD16")) {
return (new ANL_CCQE_Evt_1DQ2_nu(samplekey));
/*
ANL CC1ppip samples
*/
} else if (!name.compare("ANL_CC1ppip_XSec_1DEnu_nu") ||
!name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W14Cut") ||
!name.compare("ANL_CC1ppip_XSec_1DEnu_nu_Uncorr") ||
!name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W14Cut_Uncorr") ||
!name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W16Cut_Uncorr")) {
return (new ANL_CC1ppip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_XSec_1DQ2_nu")) {
return (new ANL_CC1ppip_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1DQ2_nu") ||
!name.compare("ANL_CC1ppip_Evt_1DQ2_nu_W14Cut")) {
return (new ANL_CC1ppip_Evt_1DQ2_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1Dppi_nu")) {
return (new ANL_CC1ppip_Evt_1Dppi_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1Dthpr_nu")) {
return (new ANL_CC1ppip_Evt_1Dthpr_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1DcosmuStar_nu")) {
return (new ANL_CC1ppip_Evt_1DcosmuStar_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1DcosthAdler_nu")) {
return (new ANL_CC1ppip_Evt_1DcosthAdler_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1Dphi_nu")) {
return (new ANL_CC1ppip_Evt_1Dphi_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1DWNpi_nu")) {
return (new ANL_CC1ppip_Evt_1DWNpi_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1DWNmu_nu")) {
return (new ANL_CC1ppip_Evt_1DWNmu_nu(samplekey));
} else if (!name.compare("ANL_CC1ppip_Evt_1DWmupi_nu")) {
return (new ANL_CC1ppip_Evt_1DWmupi_nu(samplekey));
/*
ANL CC1npip sample
*/
} else if (!name.compare("ANL_CC1npip_XSec_1DEnu_nu") ||
!name.compare("ANL_CC1npip_XSec_1DEnu_nu_W14Cut") ||
!name.compare("ANL_CC1npip_XSec_1DEnu_nu_Uncorr") ||
!name.compare("ANL_CC1npip_XSec_1DEnu_nu_W14Cut_Uncorr") ||
!name.compare("ANL_CC1npip_XSec_1DEnu_nu_W16Cut_Uncorr")) {
return (new ANL_CC1npip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC1npip_Evt_1DQ2_nu") ||
!name.compare("ANL_CC1npip_Evt_1DQ2_nu_W14Cut")) {
return (new ANL_CC1npip_Evt_1DQ2_nu(samplekey));
} else if (!name.compare("ANL_CC1npip_Evt_1Dppi_nu")) {
return (new ANL_CC1npip_Evt_1Dppi_nu(samplekey));
} else if (!name.compare("ANL_CC1npip_Evt_1DcosmuStar_nu")) {
return (new ANL_CC1npip_Evt_1DcosmuStar_nu(samplekey));
} else if (!name.compare("ANL_CC1npip_Evt_1DWNpi_nu")) {
return (new ANL_CC1npip_Evt_1DWNpi_nu(samplekey));
} else if (!name.compare("ANL_CC1npip_Evt_1DWNmu_nu")) {
return (new ANL_CC1npip_Evt_1DWNmu_nu(samplekey));
} else if (!name.compare("ANL_CC1npip_Evt_1DWmupi_nu")) {
return (new ANL_CC1npip_Evt_1DWmupi_nu(samplekey));
/*
ANL CC1pi0 sample
*/
} else if (!name.compare("ANL_CC1pi0_XSec_1DEnu_nu") ||
!name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W14Cut") ||
!name.compare("ANL_CC1pi0_XSec_1DEnu_nu_Uncorr") ||
!name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W14Cut_Uncorr") ||
!name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W16Cut_Uncorr")) {
return (new ANL_CC1pi0_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC1pi0_Evt_1DQ2_nu") ||
!name.compare("ANL_CC1pi0_Evt_1DQ2_nu_W14Cut")) {
return (new ANL_CC1pi0_Evt_1DQ2_nu(samplekey));
} else if (!name.compare("ANL_CC1pi0_Evt_1DcosmuStar_nu")) {
return (new ANL_CC1pi0_Evt_1DcosmuStar_nu(samplekey));
} else if (!name.compare("ANL_CC1pi0_Evt_1DWNpi_nu")) {
return (new ANL_CC1pi0_Evt_1DWNpi_nu(samplekey));
} else if (!name.compare("ANL_CC1pi0_Evt_1DWNmu_nu")) {
return (new ANL_CC1pi0_Evt_1DWNmu_nu(samplekey));
} else if (!name.compare("ANL_CC1pi0_Evt_1DWmupi_nu")) {
return (new ANL_CC1pi0_Evt_1DWmupi_nu(samplekey));
/*
ANL NC1npip sample
*/
} else if (!name.compare("ANL_NC1npip_Evt_1Dppi_nu")) {
return (new ANL_NC1npip_Evt_1Dppi_nu(samplekey));
/*
ANL NC1ppim sample
*/
} else if (!name.compare("ANL_NC1ppim_XSec_1DEnu_nu")) {
return (new ANL_NC1ppim_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_NC1ppim_Evt_1DcosmuStar_nu")) {
return (new ANL_NC1ppim_Evt_1DcosmuStar_nu(samplekey));
/*
ANL CC2pi sample
*/
} else if (!name.compare("ANL_CC2pi_1pim1pip_XSec_1DEnu_nu")) {
return (new ANL_CC2pi_1pim1pip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu")) {
return (new ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dppip_nu")) {
return (new ANL_CC2pi_1pim1pip_Evt_1Dppip_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dppim_nu")) {
return (new ANL_CC2pi_1pim1pip_Evt_1Dppim_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu")) {
return (new ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pip_XSec_1DEnu_nu")) {
return (new ANL_CC2pi_1pip1pip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu")) {
return (new ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu")) {
return (new ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu")) {
return (new ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu")) {
return (new ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu")) {
return (new ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu")) {
return (new ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu")) {
return (new ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu")) {
return (new ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu(samplekey));
} else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu")) {
return (new ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu(samplekey));
/*
ArgoNeut Samples
*/
} else
#endif
#ifndef __NO_ArgoNeuT__
if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_antinu")) {
return (new ArgoNeuT_CCInc_XSec_1Dpmu_antinu(samplekey));
} else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_nu")) {
return (new ArgoNeuT_CCInc_XSec_1Dpmu_nu(samplekey));
} else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_antinu")) {
return (new ArgoNeuT_CCInc_XSec_1Dthetamu_antinu(samplekey));
} else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_nu")) {
return (new ArgoNeuT_CCInc_XSec_1Dthetamu_nu(samplekey));
} else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dpmu_nu")) {
return (new ArgoNeuT_CC1Pi_XSec_1Dpmu_nu(samplekey));
} else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamu_nu")) {
return (new ArgoNeuT_CC1Pi_XSec_1Dthetamu_nu(samplekey));
} else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetapi_nu")) {
return (new ArgoNeuT_CC1Pi_XSec_1Dthetapi_nu(samplekey));
} else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamupi_nu")) {
return (new ArgoNeuT_CC1Pi_XSec_1Dthetamupi_nu(samplekey));
} else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dpmu_antinu")) {
return (new ArgoNeuT_CC1Pi_XSec_1Dpmu_antinu(samplekey));
} else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamu_antinu")) {
return (new ArgoNeuT_CC1Pi_XSec_1Dthetamu_antinu(samplekey));
} else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetapi_antinu")) {
return (new ArgoNeuT_CC1Pi_XSec_1Dthetapi_antinu(samplekey));
} else if (!name.compare("ArgoNeuT_CC1Pi_XSec_1Dthetamupi_antinu")) {
return (new ArgoNeuT_CC1Pi_XSec_1Dthetamupi_antinu(samplekey));
/*
BNL Samples
*/
} else
#endif
#ifndef __NO_BNL__
if (!name.compare("BNL_CCQE_XSec_1DEnu_nu")) {
return (new BNL_CCQE_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BNL_CCQE_Evt_1DQ2_nu")) {
return (new BNL_CCQE_Evt_1DQ2_nu(samplekey));
/*
BNL CC1ppip samples
*/
} else if (!name.compare("BNL_CC1ppip_XSec_1DEnu_nu") ||
!name.compare("BNL_CC1ppip_XSec_1DEnu_nu_Uncorr") ||
!name.compare("BNL_CC1ppip_XSec_1DEnu_nu_W14Cut") ||
!name.compare("BNL_CC1ppip_XSec_1DEnu_nu_W14Cut_Uncorr")) {
return (new BNL_CC1ppip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BNL_CC1ppip_Evt_1DQ2_nu") ||
!name.compare("BNL_CC1ppip_Evt_1DQ2_nu_W14Cut")) {
return (new BNL_CC1ppip_Evt_1DQ2_nu(samplekey));
} else if (!name.compare("BNL_CC1ppip_Evt_1DcosthAdler_nu")) {
return (new BNL_CC1ppip_Evt_1DcosthAdler_nu(samplekey));
} else if (!name.compare("BNL_CC1ppip_Evt_1Dphi_nu")) {
return (new BNL_CC1ppip_Evt_1Dphi_nu(samplekey));
} else if (!name.compare("BNL_CC1ppip_Evt_1DWNpi_nu")) {
return (new BNL_CC1ppip_Evt_1DWNpi_nu(samplekey));
} else if (!name.compare("BNL_CC1ppip_Evt_1DWNmu_nu")) {
return (new BNL_CC1ppip_Evt_1DWNmu_nu(samplekey));
} else if (!name.compare("BNL_CC1ppip_Evt_1DWmupi_nu")) {
return (new BNL_CC1ppip_Evt_1DWmupi_nu(samplekey));
/*
BNL CC1npip samples
*/
} else if (!name.compare("BNL_CC1npip_XSec_1DEnu_nu") ||
!name.compare("BNL_CC1npip_XSec_1DEnu_nu_Uncorr")) {
return (new BNL_CC1npip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BNL_CC1npip_Evt_1DQ2_nu")) {
return (new BNL_CC1npip_Evt_1DQ2_nu(samplekey));
} else if (!name.compare("BNL_CC1npip_Evt_1DWNpi_nu")) {
return (new BNL_CC1npip_Evt_1DWNpi_nu(samplekey));
} else if (!name.compare("BNL_CC1npip_Evt_1DWNmu_nu")) {
return (new BNL_CC1npip_Evt_1DWNmu_nu(samplekey));
} else if (!name.compare("BNL_CC1npip_Evt_1DWmupi_nu")) {
return (new BNL_CC1npip_Evt_1DWmupi_nu(samplekey));
/*
BNL CC1pi0 samples
*/
} else if (!name.compare("BNL_CC1pi0_XSec_1DEnu_nu")) {
return (new BNL_CC1pi0_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BNL_CC1pi0_Evt_1DQ2_nu")) {
return (new BNL_CC1pi0_Evt_1DQ2_nu(samplekey));
} else if (!name.compare("BNL_CC1pi0_Evt_1DWNpi_nu")) {
return (new BNL_CC1pi0_Evt_1DWNpi_nu(samplekey));
} else if (!name.compare("BNL_CC1pi0_Evt_1DWNmu_nu")) {
return (new BNL_CC1pi0_Evt_1DWNmu_nu(samplekey));
} else if (!name.compare("BNL_CC1pi0_Evt_1DWmupi_nu")) {
return (new BNL_CC1pi0_Evt_1DWmupi_nu(samplekey));
/*
BNL multi-pi
*/
} else if (!name.compare("BNL_CC2pi_1pim1pip_XSec_1DEnu_nu")) {
return (new BNL_CC2pi_1pim1pip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BNL_CC3pi_1pim2pip_XSec_1DEnu_nu")) {
return (new BNL_CC3pi_1pim2pip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BNL_CC4pi_2pim2pip_XSec_1DEnu_nu")) {
return (new BNL_CC4pi_2pim2pip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BNL_CC2pi_1pim1pip_Evt_1DWpippim_nu")) {
return (new BNL_CC2pi_1pim1pip_Evt_1DWpippim_nu(samplekey));
} else if (!name.compare("BNL_CC2pi_1pim1pip_Evt_1DWpippr_nu")) {
return (new BNL_CC2pi_1pim1pip_Evt_1DWpippr_nu(samplekey));
/*
FNAL Samples
*/
} else
#endif
#ifndef __NO_FNAL__
if (!name.compare("FNAL_CCQE_Evt_1DQ2_nu")) {
return (new FNAL_CCQE_Evt_1DQ2_nu(samplekey));
/*
FNAL CC1ppip
*/
} else if (!name.compare("FNAL_CC1ppip_XSec_1DEnu_nu")) {
return (new FNAL_CC1ppip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("FNAL_CC1ppip_XSec_1DQ2_nu")) {
return (new FNAL_CC1ppip_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("FNAL_CC1ppip_Evt_1DQ2_nu")) {
return (new FNAL_CC1ppip_Evt_1DQ2_nu(samplekey));
/*
FNAL CC1ppim
*/
} else if (!name.compare("FNAL_CC1ppim_XSec_1DEnu_antinu")) {
return (new FNAL_CC1ppim_XSec_1DEnu_antinu(samplekey));
/*
BEBC Samples
*/
} else
#endif
#ifndef __NO_BEBC__
if (!name.compare("BEBC_CCQE_XSec_1DQ2_nu")) {
return (new BEBC_CCQE_XSec_1DQ2_nu(samplekey));
/*
BEBC CC1ppip samples
*/
} else if (!name.compare("BEBC_CC1ppip_XSec_1DEnu_nu")) {
return (new BEBC_CC1ppip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BEBC_CC1ppip_XSec_1DQ2_nu")) {
return (new BEBC_CC1ppip_XSec_1DQ2_nu(samplekey));
/*
BEBC CC1npip samples
*/
} else if (!name.compare("BEBC_CC1npip_XSec_1DEnu_nu")) {
return (new BEBC_CC1npip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BEBC_CC1npip_XSec_1DQ2_nu")) {
return (new BEBC_CC1npip_XSec_1DQ2_nu(samplekey));
/*
BEBC CC1pi0 samples
*/
} else if (!name.compare("BEBC_CC1pi0_XSec_1DEnu_nu")) {
return (new BEBC_CC1pi0_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("BEBC_CC1pi0_XSec_1DQ2_nu")) {
return (new BEBC_CC1pi0_XSec_1DQ2_nu(samplekey));
/*
BEBC CC1npim samples
*/
} else if (!name.compare("BEBC_CC1npim_XSec_1DEnu_antinu")) {
return (new BEBC_CC1npim_XSec_1DEnu_antinu(samplekey));
} else if (!name.compare("BEBC_CC1npim_XSec_1DQ2_antinu")) {
return (new BEBC_CC1npim_XSec_1DQ2_antinu(samplekey));
/*
BEBC CC1ppim samples
*/
} else if (!name.compare("BEBC_CC1ppim_XSec_1DEnu_antinu")) {
return (new BEBC_CC1ppim_XSec_1DEnu_antinu(samplekey));
} else if (!name.compare("BEBC_CC1ppim_XSec_1DQ2_antinu")) {
return (new BEBC_CC1ppim_XSec_1DQ2_antinu(samplekey));
/*
GGM CC1ppip samples
*/
} else
#endif
#ifndef __NO_GGM__
if (!name.compare("GGM_CC1ppip_XSec_1DEnu_nu")) {
return (new GGM_CC1ppip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("GGM_CC1ppip_Evt_1DQ2_nu")) {
return (new GGM_CC1ppip_Evt_1DQ2_nu(samplekey));
/*
MiniBooNE Samples
*/
/*
CCQE
*/
} else
#endif
#ifndef __NO_MiniBooNE__
if (!name.compare("MiniBooNE_CCQE_XSec_1DQ2_nu") ||
!name.compare("MiniBooNE_CCQELike_XSec_1DQ2_nu")) {
return (new MiniBooNE_CCQE_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("MiniBooNE_CCQE_XSec_1DEnu_nu") ||
!name.compare("MiniBooNE_CCQELike_XSec_1DEnu_nu")) {
return (new MiniBooNE_CCQE_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CCQE_XSec_1DQ2_antinu") ||
!name.compare("MiniBooNE_CCQELike_XSec_1DQ2_antinu") ||
!name.compare("MiniBooNE_CCQE_CTarg_XSec_1DQ2_antinu")) {
return (new MiniBooNE_CCQE_XSec_1DQ2_antinu(samplekey));
} else if (!name.compare("MiniBooNE_CCQE_XSec_2DTcos_nu") ||
!name.compare("MiniBooNE_CCQELike_XSec_2DTcos_nu")) {
return (new MiniBooNE_CCQE_XSec_2DTcos_nu(samplekey));
} else if (!name.compare("MiniBooNE_CCQE_XSec_2DTcos_antinu") ||
!name.compare("MiniBooNE_CCQELike_XSec_2DTcos_antinu")) {
return (new MiniBooNE_CCQE_XSec_2DTcos_antinu(samplekey));
/*
MiniBooNE CC1pi+
*/
// 1D
} else if (!name.compare("MiniBooNE_CC1pip_XSec_1DEnu_nu")) {
return (new MiniBooNE_CC1pip_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_1DQ2_nu")) {
return (new MiniBooNE_CC1pip_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_1DTpi_nu")) {
return (new MiniBooNE_CC1pip_XSec_1DTpi_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_1DTu_nu")) {
return (new MiniBooNE_CC1pip_XSec_1DTu_nu(samplekey));
// 2D
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DQ2Enu_nu")) {
return (new MiniBooNE_CC1pip_XSec_2DQ2Enu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTpiCospi_nu")) {
return (new MiniBooNE_CC1pip_XSec_2DTpiCospi_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTpiEnu_nu")) {
return (new MiniBooNE_CC1pip_XSec_2DTpiEnu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTuCosmu_nu")) {
return (new MiniBooNE_CC1pip_XSec_2DTuCosmu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTuEnu_nu")) {
return (new MiniBooNE_CC1pip_XSec_2DTuEnu_nu(samplekey));
/*
MiniBooNE CC1pi0
*/
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DEnu_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DQ2_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DTu_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1DTu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dcosmu_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1Dcosmu_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dcospi0_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1Dcospi0_nu(samplekey));
} else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dppi0_nu")) {
return (new MiniBooNE_CC1pi0_XSec_1Dppi0_nu(samplekey));
} else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu") ||
!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_rhc")) {
return (new MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu(samplekey));
} else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_nu") ||
!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_fhc")) {
return (new MiniBooNE_NC1pi0_XSec_1Dcospi0_nu(samplekey));
} else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_antinu") ||
!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_rhc")) {
return (new MiniBooNE_NC1pi0_XSec_1Dppi0_antinu(samplekey));
} else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_nu") ||
!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_fhc")) {
return (new MiniBooNE_NC1pi0_XSec_1Dppi0_nu(samplekey));
/*
MiniBooNE NCEL
*/
} else if (!name.compare("MiniBooNE_NCEL_XSec_Treco_nu")) {
return (new MiniBooNE_NCEL_XSec_Treco_nu(samplekey));
} else
#endif
#ifndef __NO_MicroBooNE__
/*
MicroBooNE Samples
*/
if (!name.compare("MicroBooNE_CCInc_XSec_2DPcos_nu")) {
return (new MicroBooNE_CCInc_XSec_2DPcos_nu(samplekey));
} else if (!name.compare("MicroBooNE_CC1MuNp_XSec_1DPmu_nu") ||
!name.compare("MicroBooNE_CC1MuNp_XSec_1Dcosmu_nu") ||
!name.compare("MicroBooNE_CC1MuNp_XSec_1DPp_nu") ||
!name.compare("MicroBooNE_CC1MuNp_XSec_1Dcosp_nu") ||
!name.compare("MicroBooNE_CC1MuNp_XSec_1Dthetamup_nu")) {
return (new MicroBooNE_CC1MuNp_XSec_1D_nu(samplekey));
} else
#endif
#ifndef __NO_MINERvA__
/*
MINERvA Samples
*/
if (!name.compare("MINERvA_CCQE_XSec_1DQ2_nu") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_nu_oldflux") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg_oldflux")) {
return (new MINERvA_CCQE_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_oldflux") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg_oldflux")) {
return (new MINERvA_CCQE_XSec_1DQ2_antinu(samplekey));
} else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_joint_oldflux") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg_oldflux") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_joint") ||
!name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg")) {
return (new MINERvA_CCQE_XSec_1DQ2_joint(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DEe_nue")) {
return (new MINERvA_CC0pi_XSec_1DEe_nue(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nue")) {
return (new MINERvA_CC0pi_XSec_1DQ2_nue(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DThetae_nue")) {
return (new MINERvA_CC0pi_XSec_1DThetae_nue(samplekey));
} else if (!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpmu_nu") ||
!name.compare("MINERvA_CC0pinp_STV_XSec_1Dthmu_nu") ||
!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpprot_nu") ||
!name.compare("MINERvA_CC0pinp_STV_XSec_1Dthprot_nu") ||
!name.compare("MINERvA_CC0pinp_STV_XSec_1Dpnreco_nu") ||
!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddalphat_nu") ||
!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddpt_nu") ||
!name.compare("MINERvA_CC0pinp_STV_XSec_1Ddphit_nu")) {
return (new MINERvA_CC0pinp_STV_XSec_1D_nu(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nu_proton")) {
return (new MINERvA_CC0pi_XSec_1DQ2_nu_proton(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtC_nu") ||
!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtCH_nu") ||
!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtFe_nu") ||
!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtPb_nu")) {
return (new MINERvA_CC0pi_XSec_1DQ2_Tgt_nu(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioC_nu") ||
!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioFe_nu") ||
!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioPb_nu")) {
return (new MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu(samplekey));
// Dan Ruterbories measurements of late 2018
} else if (!name.compare("MINERvA_CC0pi_XSec_2Dptpz_nu")) {
return (new MINERvA_CC0pi_XSec_2D_nu(samplekey));
// } else if (!name.compare("MINERvA_CC0pi_XSec_3DptpzTp_nu")) {
// return (new MINERvA_CC0pi_XSec_3DptpzTp_nu(samplekey));
} else if (!name.compare("MINERvA_CC0pi_XSec_1Dpt_nu") ||
!name.compare("MINERvA_CC0pi_XSec_1Dpz_nu") ||
!name.compare("MINERvA_CC0pi_XSec_1DQ2QE_nu") ||
!name.compare("MINERvA_CC0pi_XSec_1DEnuQE_nu")) {
return (new MINERvA_CC0pi_XSec_1D_2018_nu(samplekey));
// C. Patrick's early 2018 measurements
} else if (!name.compare("MINERvA_CC0pi_XSec_2Dptpz_antinu") ||
!name.compare("MINERvA_CC0pi_XSec_2DQ2QEEnuQE_antinu") ||
!name.compare("MINERvA_CC0pi_XSec_2DQ2QEEnuTrue_antinu")) {
return (new MINERvA_CC0pi_XSec_2D_antinu(samplekey));
/*
CC1pi+
*/
// DONE
} else if (!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu") ||
!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_20deg") ||
!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_fluxcorr") ||
!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_20deg_fluxcorr")) {
return (new MINERvA_CC1pip_XSec_1DTpi_nu(samplekey));
// DONE
} else if (!name.compare("MINERvA_CC1pip_XSec_1Dth_nu") ||
!name.compare("MINERvA_CC1pip_XSec_1Dth_nu_20deg") ||
!name.compare("MINERvA_CC1pip_XSec_1Dth_nu_fluxcorr") ||
!name.compare("MINERvA_CC1pip_XSec_1Dth_nu_20deg_fluxcorr")) {
return (new MINERvA_CC1pip_XSec_1Dth_nu(samplekey));
} else if (!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_2017") ||
!name.compare("MINERvA_CC1pip_XSec_1Dth_nu_2017") ||
!name.compare("MINERvA_CC1pip_XSec_1Dpmu_nu_2017") ||
!name.compare("MINERvA_CC1pip_XSec_1Dthmu_nu_2017") ||
!name.compare("MINERvA_CC1pip_XSec_1DQ2_nu_2017") ||
!name.compare("MINERvA_CC1pip_XSec_1DEnu_nu_2017")) {
return (new MINERvA_CC1pip_XSec_1D_2017Update(samplekey));
/*
CC1pi-
*/
} else if (!name.compare("MINERvA_CC1pim_XSec_1DEnu_antinu")) {
return (new MINERvA_CC1pim_XSec_1DEnu_antinu(samplekey));
} else if (!name.compare("MINERvA_CC1pim_XSec_1DQ2_antinu")) {
return (new MINERvA_CC1pim_XSec_1DQ2_antinu(samplekey));
} else if (!name.compare("MINERvA_CC1pim_XSec_1DTpi_antinu")) {
return (new MINERvA_CC1pim_XSec_1DTpi_antinu(samplekey));
} else if (!name.compare("MINERvA_CC1pim_XSec_1Dpmu_antinu")) {
return (new MINERvA_CC1pim_XSec_1Dpmu_antinu(samplekey));
} else if (!name.compare("MINERvA_CC1pim_XSec_1Dth_antinu")) {
return (new MINERvA_CC1pim_XSec_1Dth_antinu(samplekey));
} else if (!name.compare("MINERvA_CC1pim_XSec_1Dthmu_antinu")) {
return (new MINERvA_CC1pim_XSec_1Dthmu_antinu(samplekey));
/*
CCNpi+
*/
} else if (!name.compare("MINERvA_CCNpip_XSec_1Dth_nu") ||
!name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015") ||
!name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2016") ||
!name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_20deg") ||
!name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_fluxcorr") ||
!name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_20deg_fluxcorr")) {
return (new MINERvA_CCNpip_XSec_1Dth_nu(samplekey));
} else if (!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu") ||
!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015") ||
!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2016") ||
!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_20deg") ||
!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_fluxcorr") ||
!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_20deg_fluxcorr")) {
return (new MINERvA_CCNpip_XSec_1DTpi_nu(samplekey));
} else if (!name.compare("MINERvA_CCNpip_XSec_1Dthmu_nu")) {
return (new MINERvA_CCNpip_XSec_1Dthmu_nu(samplekey));
} else if (!name.compare("MINERvA_CCNpip_XSec_1Dpmu_nu")) {
return (new MINERvA_CCNpip_XSec_1Dpmu_nu(samplekey));
} else if (!name.compare("MINERvA_CCNpip_XSec_1DQ2_nu")) {
return (new MINERvA_CCNpip_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("MINERvA_CCNpip_XSec_1DEnu_nu")) {
return (new MINERvA_CCNpip_XSec_1DEnu_nu(samplekey));
/*
MINERvA CC1pi0 anti-nu
*/
// Done
} else if (!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2015") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2016") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_fluxcorr") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2015_fluxcorr") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2016_fluxcorr")) {
return (new MINERvA_CC1pi0_XSec_1Dth_antinu(samplekey));
} else if (!name.compare("MINERvA_CC1pi0_XSec_1Dppi0_antinu") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dppi0_antinu_fluxcorr")) {
return (new MINERvA_CC1pi0_XSec_1Dppi0_antinu(samplekey));
} else if (!name.compare("MINERvA_CC1pi0_XSec_1DTpi0_antinu")) {
return (new MINERvA_CC1pi0_XSec_1DTpi0_antinu(samplekey));
// Done
} else if (!name.compare("MINERvA_CC1pi0_XSec_1DQ2_antinu")) {
return (new MINERvA_CC1pi0_XSec_1DQ2_antinu(samplekey));
// Done
} else if (!name.compare("MINERvA_CC1pi0_XSec_1Dthmu_antinu")) {
return (new MINERvA_CC1pi0_XSec_1Dthmu_antinu(samplekey));
// Done
} else if (!name.compare("MINERvA_CC1pi0_XSec_1Dpmu_antinu")) {
return (new MINERvA_CC1pi0_XSec_1Dpmu_antinu(samplekey));
// Done
} else if (!name.compare("MINERvA_CC1pi0_XSec_1DEnu_antinu")) {
return (new MINERvA_CC1pi0_XSec_1DEnu_antinu(samplekey));
// MINERvA CC1pi0 nu
} else if (!name.compare("MINERvA_CC1pi0_XSec_1DTpi_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dth_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dpmu_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1Dthmu_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1DQ2_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1DEnu_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1DWexp_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1DPPi0Mass_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1DPPi0MassDelta_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1DCosAdler_nu") ||
!name.compare("MINERvA_CC1pi0_XSec_1DPhiAdler_nu")) {
return (new MINERvA_CC1pi0_XSec_1D_nu(samplekey));
/*
CCINC
*/
} else if (!name.compare("MINERvA_CCinc_XSec_2DEavq3_nu")) {
return (new MINERvA_CCinc_XSec_2DEavq3_nu(samplekey));
} else if (!name.compare("MINERvA_CCinc_XSec_1Dx_ratio_C12_CH") ||
!name.compare("MINERvA_CCinc_XSec_1Dx_ratio_Fe56_CH") ||
!name.compare("MINERvA_CCinc_XSec_1Dx_ratio_Pb208_CH")) {
return (new MINERvA_CCinc_XSec_1Dx_ratio(samplekey));
} else if (!name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_C12_CH") ||
!name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_Fe56_CH") ||
!name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_Pb208_CH")) {
return (new MINERvA_CCinc_XSec_1DEnu_ratio(samplekey));
/*
CCDIS
*/
} else if (!name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_C12_CH") ||
!name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_Fe56_CH") ||
!name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_Pb208_CH")) {
return (new MINERvA_CCDIS_XSec_1Dx_ratio(samplekey));
} else if (!name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_C12_CH") ||
!name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_Fe56_CH") ||
!name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_Pb208_CH")) {
return (new MINERvA_CCDIS_XSec_1DEnu_ratio(samplekey));
/*
CC-COH
*/
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_nu")) {
return (new MINERvA_CCCOHPI_XSec_1DEnu_nu(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_nu")) {
return (new MINERvA_CCCOHPI_XSec_1DEpi_nu(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_nu")) {
return (new MINERvA_CCCOHPI_XSec_1Dth_nu(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_nu")) {
return (new MINERvA_CCCOHPI_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_antinu")) {
return (new MINERvA_CCCOHPI_XSec_1DEnu_antinu(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_antinu")) {
return (new MINERvA_CCCOHPI_XSec_1DEpi_antinu(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_antinu")) {
return (new MINERvA_CCCOHPI_XSec_1Dth_antinu(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_antinu")) {
return (new MINERvA_CCCOHPI_XSec_1DQ2_antinu(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_joint")) {
return (new MINERvA_CCCOHPI_XSec_joint(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_joint")) {
return (new MINERvA_CCCOHPI_XSec_joint(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_joint")) {
return (new MINERvA_CCCOHPI_XSec_joint(samplekey));
} else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_joint")) {
return (new MINERvA_CCCOHPI_XSec_joint(samplekey));
/*
T2K Samples
*/
} else
#endif
#ifndef __NO_T2K__
if (!name.compare("T2K_CC0pi_XSec_2DPcos_nu_I")) {
return (new T2K_CC0pi_XSec_2DPcos_nu_I(samplekey));
} else if (!name.compare("T2K_CC0pi_XSec_2DPcos_nu_II")) {
return (new T2K_CC0pi_XSec_2DPcos_nu_II(samplekey));
} else if (!name.compare("T2K_CCinc_XSec_2DPcos_nu_nonuniform")) {
return (new T2K_CCinc_XSec_2DPcos_nu_nonuniform(samplekey));
} else if (!name.compare("T2K_CC0pi_XSec_H2O_2DPcos_anu")) {
return (new T2K_CC0pi_XSec_H2O_2DPcos_anu(samplekey));
} else if (!name.compare("T2K_NuMu_CC0pi_O_XSec_2DPcos") ||
!name.compare("T2K_NuMu_CC0pi_C_XSec_2DPcos")) {
return (new T2K_NuMu_CC0pi_OC_XSec_2DPcos(samplekey));
} else if (!name.compare("T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint")) {
return (new T2K_NuMu_CC0pi_OC_XSec_2DPcos_joint(samplekey));
} else if (!name.compare("T2K_NuMu_CC0pi_CH_XSec_2DPcos") ||
!name.compare("T2K_AntiNuMu_CC0pi_CH_XSec_2DPcos")) {
return (new T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos(samplekey));
} else if (!name.compare("T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint")) {
return (new T2K_NuMuAntiNuMu_CC0pi_CH_XSec_2DPcos_joint(samplekey));
} else if (!name.compare("T2K_nueCCinc_XSec_1Dpe_FHC") ||
!name.compare("T2K_nueCCinc_XSec_1Dpe_RHC") ||
!name.compare("T2K_nuebarCCinc_XSec_1Dpe_RHC")) {
return (new T2K_nueCCinc_XSec_1Dpe(samplekey));
} else if (!name.compare("T2K_nueCCinc_XSec_1Dthe_FHC") ||
!name.compare("T2K_nueCCinc_XSec_1Dthe_RHC") ||
!name.compare("T2K_nuebarCCinc_XSec_1Dthe_RHC")) {
return (new T2K_nueCCinc_XSec_1Dthe(samplekey));
} else if (!name.compare("T2K_nueCCinc_XSec_1Dpe_joint")) {
return (new T2K_nueCCinc_XSec_1Dpe_joint(samplekey));
} else if (!name.compare("T2K_nueCCinc_XSec_1Dthe_joint")) {
return (new T2K_nueCCinc_XSec_1Dthe_joint(samplekey));
} else if (!name.compare("T2K_nueCCinc_XSec_joint")) {
return (new T2K_nueCCinc_XSec_joint(samplekey));
/*
T2K CC1pi+ CH samples
*/
// Comment these out for now because we don't have the proper data
} else if (!name.compare("T2K_CC1pip_CH_XSec_2Dpmucosmu_nu")) {
return (new T2K_CC1pip_CH_XSec_2Dpmucosmu_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dppi_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dppi_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthpi_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dthpi_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthmupi_nu")) {
return (new T2K_CC1pip_CH_XSec_1Dthmupi_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1DQ2_nu")) {
return (new T2K_CC1pip_CH_XSec_1DQ2_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1DAdlerPhi_nu")) {
return (new T2K_CC1pip_CH_XSec_1DAdlerPhi_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_CH_XSec_1DCosThAdler_nu")) {
return (new T2K_CC1pip_CH_XSec_1DCosThAdler_nu(samplekey));
/*
T2K CC1pi+ H2O samples
*/
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1DEnuDelta_nu")) {
return (new T2K_CC1pip_H2O_XSec_1DEnuDelta_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1DEnuMB_nu")) {
return (new T2K_CC1pip_H2O_XSec_1DEnuMB_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcosmu_nu")) {
return (new T2K_CC1pip_H2O_XSec_1Dcosmu_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcosmupi_nu")) {
return (new T2K_CC1pip_H2O_XSec_1Dcosmupi_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcospi_nu")) {
return (new T2K_CC1pip_H2O_XSec_1Dcospi_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dpmu_nu")) {
return (new T2K_CC1pip_H2O_XSec_1Dpmu_nu(samplekey));
} else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dppi_nu")) {
return (new T2K_CC1pip_H2O_XSec_1Dppi_nu(samplekey));
/*
T2K CC0pi + np CH samples
*/
} else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddpt_nu")) {
return (new T2K_CC0pinp_STV_XSec_1Ddpt_nu(samplekey));
} else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddphit_nu")) {
return (new T2K_CC0pinp_STV_XSec_1Ddphit_nu(samplekey));
} else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddat_nu")) {
return (new T2K_CC0pinp_STV_XSec_1Ddat_nu(samplekey));
} else if (!name.compare("T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p_Np") || !name.compare("T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p") || !name.compare("T2K_CC0piWithProtons_XSec_2018_multidif_0p") || !name.compare("T2K_CC0piWithProtons_XSec_2018_multidif_1p")) {
return (new T2K_CC0piWithProtons_XSec_2018_multidif_0p_1p_Np(samplekey));
} else if (!name.compare("T2K_CC0pinp_ifk_XSec_3Dinfp_nu")) {
return (new T2K_CC0pinp_ifk_XSec_3Dinfp_nu(samplekey));
} else if (!name.compare("T2K_CC0pinp_ifk_XSec_3Dinfa_nu")) {
return (new T2K_CC0pinp_ifk_XSec_3Dinfa_nu(samplekey));
} else if (!name.compare("T2K_CC0pinp_ifk_XSec_3Dinfip_nu")) {
return (new T2K_CC0pinp_ifk_XSec_3Dinfip_nu(samplekey));
// SciBooNE COH studies
} else
#endif
#ifndef __NO_SciBooNE__
if (!name.compare("SciBooNE_CCCOH_STOP_NTrks_nu")) {
return (new SciBooNE_CCCOH_STOP_NTrks_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_1TRK_1DQ2_nu")) {
return (new SciBooNE_CCCOH_1TRK_1DQ2_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_1TRK_1Dpmu_nu")) {
return (new SciBooNE_CCCOH_1TRK_1Dpmu_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_1TRK_1Dthetamu_nu")) {
return (new SciBooNE_CCCOH_1TRK_1Dthetamu_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPr_1DQ2_nu")) {
return (new SciBooNE_CCCOH_MuPr_1DQ2_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPr_1Dpmu_nu")) {
return (new SciBooNE_CCCOH_MuPr_1Dpmu_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPr_1Dthetamu_nu")) {
return (new SciBooNE_CCCOH_MuPr_1Dthetamu_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPiVA_1DQ2_nu")) {
return (new SciBooNE_CCCOH_MuPiVA_1DQ2_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPiVA_1Dpmu_nu")) {
return (new SciBooNE_CCCOH_MuPiVA_1Dpmu_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPiVA_1Dthetamu_nu")) {
return (new SciBooNE_CCCOH_MuPiVA_1Dthetamu_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu")) {
return (new SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu")) {
return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu")) {
return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dthetamu_nu")) {
return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetamu_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dpmu_nu")) {
return (new SciBooNE_CCCOH_MuPiNoVA_1Dpmu_nu(samplekey));
} else if (!name.compare("SciBooNE_CCCOH_STOPFINAL_1DQ2_nu")) {
return (new SciBooNE_CCCOH_STOPFINAL_1DQ2_nu(samplekey));
} else if (!name.compare("SciBooNE_CCInc_XSec_1DEnu_nu") ||
!name.compare("SciBooNE_CCInc_XSec_1DEnu_nu_NEUT") ||
!name.compare("SciBooNE_CCInc_XSec_1DEnu_nu_NUANCE")) {
return (new SciBooNE_CCInc_XSec_1DEnu_nu(samplekey));
/*
K2K Samples
*/
/*
NC1pi0
*/
} else
#endif
#ifndef __NO_K2K__
if (!name.compare("K2K_NC1pi0_Evt_1Dppi0_nu")) {
return (new K2K_NC1pi0_Evt_1Dppi0_nu(samplekey));
/*
Fake Studies
*/
} else
#endif
if (name.find("ExpMultDist_CCQE_XSec_1D") != std::string::npos &&
name.find("_FakeStudy") != std::string::npos) {
return (
new ExpMultDist_CCQE_XSec_1DVar_FakeStudy(name, file, rw, type, fkdt));
} else if (name.find("ExpMultDist_CCQE_XSec_2D") != std::string::npos &&
name.find("_FakeStudy") != std::string::npos) {
return (
new ExpMultDist_CCQE_XSec_2DVar_FakeStudy(name, file, rw, type, fkdt));
} else if (name.find("GenericFlux") != std::string::npos) {
return (new GenericFlux_Tester(name, file, rw, type, fkdt));
} else if (name.find("GenericVectors") != std::string::npos) {
return (new GenericFlux_Vectors(name, file, rw, type, fkdt));
} else if (!name.compare("T2K2017_FakeData")) {
return (new T2K2017_FakeData(samplekey));
} else if (!name.compare("MCStudy_CCQE")) {
return (new MCStudy_CCQEHistograms(name, file, rw, type, fkdt));
} else if (!name.compare("ElectronFlux_FlatTree")) {
return (new ElectronFlux_FlatTree(name, file, rw, type, fkdt));
} else if (name.find("ElectronData_") != std::string::npos) {
return new ElectronScattering_DurhamData(samplekey);
} else if (name.find("MuonValidation_") != std::string::npos) {
return (new MCStudy_MuonValidation(name, file, rw, type, fkdt));
} else if (!name.compare("NIWGOfficialPlots")) {
return (new OfficialNIWGPlots(samplekey));
} else if ((name.find("SigmaEnuHists") != std::string::npos) ||
(name.find("SigmaEnuPerEHists") != std::string::npos)) {
return (new SigmaEnuHists(samplekey));
} else if (!name.compare("Simple_Osc")) {
return (new Simple_Osc(samplekey));
} else if (!name.compare("Smear_SVDUnfold_Propagation_Osc")) {
return (new Smear_SVDUnfold_Propagation_Osc(samplekey));
} else {
NUIS_ABORT("Error: No such sample: " << name << std::endl);
}
// Return NULL if no sample loaded.
return NULL;
} // namespace SampleUtils
} // namespace SampleUtils
diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx
index 36b63d2..60ef26f 100644
--- a/src/InputHandler/GENIEInputHandler.cxx
+++ b/src/InputHandler/GENIEInputHandler.cxx
@@ -1,635 +1,636 @@
// Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifdef __GENIE_ENABLED__
#include "GENIEInputHandler.h"
#ifdef GENIE_PRE_R3
#include "Messenger/Messenger.h"
#else
#include "Framework/Messenger/Messenger.h"
#endif
#include "InputUtils.h"
GENIEGeneratorInfo::~GENIEGeneratorInfo() { DeallocateParticleStack(); }
void GENIEGeneratorInfo::AddBranchesToTree(TTree *tn) {
tn->Branch("GenieParticlePDGs", &fGenieParticlePDGs, "GenieParticlePDGs/I");
}
void GENIEGeneratorInfo::SetBranchesFromTree(TTree *tn) {
tn->SetBranchAddress("GenieParticlePDGs", &fGenieParticlePDGs);
}
void GENIEGeneratorInfo::AllocateParticleStack(int stacksize) {
fGenieParticlePDGs = new int[stacksize];
}
void GENIEGeneratorInfo::DeallocateParticleStack() {
delete fGenieParticlePDGs;
}
void GENIEGeneratorInfo::FillGeneratorInfo(NtpMCEventRecord *ntpl) {
Reset();
// Check for GENIE Event
if (!ntpl)
return;
if (!ntpl->event)
return;
// Cast Event Record
GHepRecord *ghep = static_cast(ntpl->event);
if (!ghep)
return;
// Fill Particle Stack
GHepParticle *p = 0;
TObjArrayIter iter(ghep);
// Loop over all particles
int i = 0;
while ((p = (dynamic_cast((iter).Next())))) {
if (!p)
continue;
// Get PDG
fGenieParticlePDGs[i] = p->Pdg();
i++;
}
}
void GENIEGeneratorInfo::Reset() {
for (int i = 0; i < kMaxParticles; i++) {
fGenieParticlePDGs[i] = 0;
}
}
+bool GENIEInputHandler::IsPrimary(GHepParticle *p) {
+
+ // If initial state nucleon, or nucleon target, it is definintely primary!
+ if (p->Status() == genie::kIStInitialState ||
+ p->Status() == genie::kIStNucleonTarget) return true;
+
+ // Reject intermediate resonant state or pre-DIS state
+ if (p->Status() == genie::kIStDISPreFragmHadronicState || // DIS before fragmentation
+ p->Status() == genie::kIStPreDecayResonantState || // pre decay resonance state
+ p->Status() == genie::kIStIntermediateState || // intermediate state
+ p->Status() == genie::kIStDecayedState || // Decayed state
+ p->Status() == genie::kIStUndefined) return false; // undefined state
+
+ // Check if the mother is the neutrino or IS nucleon
+ if (p->FirstMother() < 2) return true;
+
+ // We've now filtered off intermediate states and obvious initial states
+
+ // Loop over this particle's mothers, grandmothers, and so on cleaning out intermediate particles
+ // Set the starting particle to be our particle
+ GHepParticle *mother = fGenieGHep->Particle(p->FirstMother());
+ while (mother->FirstMother() > 1) {
+ // It could be that the mother's status is actually a decayed state linked back to the vertex
+ if (mother->Status() == genie::kIStDecayedState || // A decayed state
+ mother->Status() == genie::kIStDISPreFragmHadronicState || // A DIS state before fragmentation
+ mother->Status() == genie::kIStPreDecayResonantState ) { // A pre-decay resonant state
+ mother = fGenieGHep->Particle(mother->FirstMother());
+ } else { // If not, move out of the loop
+ break;
+ }
+ }
+
+ // Then do a simple check of mother is associated with the primary
+ int MotherID = mother->FirstMother();
+ if (MotherID > 2) return false;
+
+ // Finally, this should mean that our partcile is marked for transport through the nucleus
+ // Could also be interactions of free proton
+ if (p->Status() == genie::kIStHadronInTheNucleus || // Then require the particle to be paseed to FSI
+ (p->Status() == genie::kIStStableFinalState && // Can also have interaction on free proton
+ fGenieGHep->Summary()->InitState().TgtPtr()->A() == 1 &&
+ fGenieGHep->Summary()->InitState().TgtPtr()->Z() == 1) ) {
+ return true;
+ }
+
+ return false;
+}
+
GENIEInputHandler::GENIEInputHandler(std::string const &handle,
- std::string const &rawinputs) {
+ std::string const &rawinputs) {
NUIS_LOG(SAM, "Creating GENIEInputHandler : " << handle);
// Plz no shouting
StopTalking();
genie::Messenger::Instance()->SetPriorityLevel("GHepUtils", pFATAL);
StartTalking();
// Shout all you want
// Run a joint input handling
fName = handle;
// Setup the TChain
fGENIETree = new TChain("gtree");
fSaveExtra = FitPar::Config().GetParB("SaveExtraGenie");
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
- // Are we running with NOvA weights
- fNOvAWeights = FitPar::Config().GetParB("NOvA_Weights");
- MAQEw = 1.0;
- NonResw = 1.0;
- RPAQEw = 1.0;
- RPARESw = 1.0;
- MECw = 1.0;
- DISw = 1.0;
- NOVAw = 1.0;
-
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
TFile *inp_file = new TFile(
InputUtils::ExpandInputDirectories(inputs[inp_it]).c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
NUIS_ABORT(
"GENIE File IsZombie() at : '"
<< inputs[inp_it] << "'" << std::endl
<< "Check that your file paths are correct and the file exists!"
<< std::endl
<< "$ ls -lh " << inputs[inp_it]);
}
// Get Flux/Event hist
TH1D *fluxhist = (TH1D *)inp_file->Get("nuisance_flux");
TH1D *eventhist = (TH1D *)inp_file->Get("nuisance_events");
if (!fluxhist or !eventhist) {
NUIS_ERR(FTL, "Input File Contents: " << inputs[inp_it]);
inp_file->ls();
NUIS_ABORT("GENIE FILE doesn't contain flux/xsec info."
- << std::endl
- << "Try running the app PrepareGENIE first on :"
- << inputs[inp_it] << std::endl
- << "$ PrepareGENIE -h");
+ << std::endl
+ << "Try running the app PrepareGENIE first on :"
+ << inputs[inp_it] << std::endl
+ << "$ PrepareGENIE -h");
}
// Get N Events
TTree *genietree = (TTree *)inp_file->Get("gtree");
if (!genietree) {
NUIS_ERR(FTL, "gtree not located in GENIE file: " << inputs[inp_it]);
NUIS_ABORT(
"Check your inputs, they may need to be completely regenerated!");
}
int nevents = genietree->GetEntries();
if (nevents <= 0) {
NUIS_ABORT("Trying to a TTree with "
- << nevents << " to TChain from : " << inputs[inp_it]);
- }
-
- // Check for precomputed weights
- TTree *weighttree = (TTree *)inp_file->Get("nova_wgts");
- if (fNOvAWeights) {
- if (!weighttree) {
- NUIS_ABORT("Did not find nova_wgts tree in file "
- << inputs[inp_it] << " but you specified it" << std::endl);
- } else {
- NUIS_LOG(FIT, "Found nova_wgts tree in file " << inputs[inp_it]);
- }
+ << nevents << " to TChain from : " << inputs[inp_it]);
}
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add To TChain
fGENIETree->AddFile(inputs[inp_it].c_str());
- if (weighttree != NULL)
- fGENIETree->AddFriend(weighttree);
}
// Registor all our file inputs
SetupJointInputs();
// Assign to tree
fEventType = kGENIE;
fGenieNtpl = NULL;
fGENIETree->SetBranchAddress("gmcrec", &fGenieNtpl);
- // Set up the custom weights
- if (fNOvAWeights) {
- fGENIETree->SetBranchAddress("MAQEwgt", &MAQEw);
- fGENIETree->SetBranchAddress("nonResNormWgt", &NonResw);
- fGENIETree->SetBranchAddress("RPAQEWgt", &RPAQEw);
- fGENIETree->SetBranchAddress("RPARESWgt", &RPARESw);
- fGENIETree->SetBranchAddress("MECWgt", &MECw);
- fGENIETree->SetBranchAddress("DISWgt", &DISw);
- fGENIETree->SetBranchAddress("nova2018CVWgt", &NOVAw);
- }
-
// Libraries should be seen but not heard...
StopTalking();
fGENIETree->GetEntry(0);
StartTalking();
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetGenieEvent(fGenieNtpl);
if (fSaveExtra) {
fGenieInfo = new GENIEGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fGenieInfo);
}
fNUISANCEEvent->HardReset();
};
GENIEInputHandler::~GENIEInputHandler() {
// if (fGenieGHep) delete fGenieGHep;
// if (fGenieNtpl) delete fGenieNtpl;
// if (fGENIETree) delete fGENIETree;
// if (fGenieInfo) delete fGenieInfo;
}
void GENIEInputHandler::CreateCache() {
if (fCacheSize > 0) {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 1);
fGENIETree->SetCacheSize(fCacheSize);
}
}
void GENIEInputHandler::RemoveCache() {
// fGENIETree->SetCacheEntryRange(0, fNEvents);
fGENIETree->AddBranchToCache("*", 0);
fGENIETree->SetCacheSize(0);
}
FitEvent *GENIEInputHandler::GetNuisanceEvent(const UInt_t ent,
- const bool lightweight) {
+ const bool lightweight) {
UInt_t entry = ent + fSkip;
if (entry >= (UInt_t)fNEvents)
return NULL;
// Clear the previous event (See Note 1 in ROOT TClonesArray documentation)
if (fGenieNtpl) {
fGenieNtpl->Clear();
}
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fGENIETree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
else {
// Check for GENIE Event
if (!fGenieNtpl)
return NULL;
if (!fGenieNtpl->event)
return NULL;
// Cast Event Record
fGenieGHep = static_cast(fGenieNtpl->event);
if (!fGenieGHep)
return NULL;
TObjArrayIter iter(fGenieGHep);
genie::GHepParticle *p;
while ((p = (dynamic_cast((iter).Next())))) {
if (!p) {
continue;
}
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
if (state != genie::kIStInitialState) {
continue;
}
fNUISANCEEvent->probe_E = p->E() * 1.E3;
fNUISANCEEvent->probe_pdg = p->Pdg();
break;
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
return fNUISANCEEvent;
}
int GENIEInputHandler::GetGENIEParticleStatus(genie::GHepParticle *p,
- int mode) {
+ int mode) {
/*
kIStUndefined = -1,
kIStInitialState = 0, / generator-level initial state /
kIStStableFinalState = 1, / generator-level final state:
particles to be tracked by detector-level MC /
kIStIntermediateState = 2,
kIStDecayedState = 3,
kIStCorrelatedNucleon = 10,
kIStNucleonTarget = 11,
kIStDISPreFragmHadronicState = 12,
kIStPreDecayResonantState = 13,
kIStHadronInTheNucleus = 14, / hadrons inside the nucleus: marked
for hadron transport modules to act on /
kIStFinalStateNuclearRemnant = 15, / low energy nuclear fragments
entering the record collectively as a 'hadronic blob' pseudo-particle /
kIStNucleonClusterTarget = 16, // for composite nucleons before
phase space decay
*/
int state = kUndefinedState;
switch (p->Status()) {
- case genie::kIStNucleonTarget:
- case genie::kIStInitialState:
- case genie::kIStCorrelatedNucleon:
- case genie::kIStNucleonClusterTarget:
- state = kInitialState;
- break;
-
- case genie::kIStStableFinalState:
- state = kFinalState;
- break;
-
- case genie::kIStHadronInTheNucleus:
- if (abs(mode) == 2)
+ case genie::kIStNucleonTarget:
+ case genie::kIStInitialState:
+ case genie::kIStCorrelatedNucleon:
+ case genie::kIStNucleonClusterTarget:
state = kInitialState;
- else
+ break;
+
+ case genie::kIStStableFinalState:
+ state = kFinalState;
+ break;
+
+ case genie::kIStHadronInTheNucleus:
+ if (abs(mode) == 2)
+ state = kInitialState;
+ else
+ state = kFSIState;
+ break;
+
+ case genie::kIStPreDecayResonantState:
+ case genie::kIStDISPreFragmHadronicState:
+ case genie::kIStIntermediateState:
state = kFSIState;
- break;
-
- case genie::kIStPreDecayResonantState:
- case genie::kIStDISPreFragmHadronicState:
- case genie::kIStIntermediateState:
- state = kFSIState;
- break;
-
- case genie::kIStFinalStateNuclearRemnant:
- case genie::kIStUndefined:
- case genie::kIStDecayedState:
- default:
- break;
+ break;
+
+ case genie::kIStFinalStateNuclearRemnant:
+ case genie::kIStUndefined:
+ case genie::kIStDecayedState:
+ default:
+ break;
}
// Flag to remove nuclear part in genie
- if (p->Pdg() > 1000000) {
- if (state == kInitialState)
- state = kNuclearInitial;
- else if (state == kFinalState)
- state = kNuclearRemnant;
+ if (kRemoveNuclearParticles) {
+ if (p->Pdg() > 1000000 &&
+ pdg::IsNeutronOrProton(
+ fGenieGHep->Summary()->InitState().TgtPtr()->HitNucPdg())) { // Check that it isn't coherent! Want to keep initial state if so
+ if (state == kInitialState)
+ state = kNuclearInitial;
+ else if (state == kFinalState)
+ state = kNuclearRemnant;
+ }
}
return state;
}
#endif
#ifdef __GENIE_ENABLED__
int GENIEInputHandler::ConvertGENIEReactionCode(GHepRecord *gheprec) {
// I randomly picked 53 here because NEUT doesn't have an appropriate mode...
if (gheprec->Summary()->ProcInfo().IsNuElectronElastic()){
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg())) return 53;
else return -53;
}
// And the same story for 54
if (gheprec->Summary()->ProcInfo().IsIMDAnnihilation()){
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg())) return 54;
else return -54;
}
// Electron Scattering
if (gheprec->Summary()->ProcInfo().IsEM()) {
if (gheprec->Summary()->InitState().ProbePdg() == 11) {
if (gheprec->Summary()->ProcInfo().IsQuasiElastic())
return 1;
else if (gheprec->Summary()->ProcInfo().IsMEC())
return 2;
else if (gheprec->Summary()->ProcInfo().IsResonant())
return 13;
else if (gheprec->Summary()->ProcInfo().IsDeepInelastic())
return 26;
else {
NUIS_ERR(WRN,
- "Unknown GENIE Electron Scattering Mode!"
- << std::endl
- << "ScatteringTypeId = "
- << gheprec->Summary()->ProcInfo().ScatteringTypeId() << " "
- << "InteractionTypeId = "
- << gheprec->Summary()->ProcInfo().InteractionTypeId()
- << std::endl
- << genie::ScatteringType::AsString(
- gheprec->Summary()->ProcInfo().ScatteringTypeId())
- << " "
- << genie::InteractionType::AsString(
- gheprec->Summary()->ProcInfo().InteractionTypeId())
- << " " << gheprec->Summary()->ProcInfo().IsMEC());
+ "Unknown GENIE Electron Scattering Mode!"
+ << std::endl
+ << "ScatteringTypeId = "
+ << gheprec->Summary()->ProcInfo().ScatteringTypeId() << " "
+ << "InteractionTypeId = "
+ << gheprec->Summary()->ProcInfo().InteractionTypeId()
+ << std::endl
+ << genie::ScatteringType::AsString(
+ gheprec->Summary()->ProcInfo().ScatteringTypeId())
+ << " "
+ << genie::InteractionType::AsString(
+ gheprec->Summary()->ProcInfo().InteractionTypeId())
+ << " " << gheprec->Summary()->ProcInfo().IsMEC());
return 0;
}
}
// Weak CC
} else if (gheprec->Summary()->ProcInfo().IsWeakCC()) {
// CC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 2;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -2;
#ifndef GENIE_PRE_R3
} else if (gheprec->Summary()->ProcInfo().IsDiffractive()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 15;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -15;
#endif
// CC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
// Weak NC
} else if (gheprec->Summary()->ProcInfo().IsWeakNC()) {
// NC MEC
if (gheprec->Summary()->ProcInfo().IsMEC()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 32;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -32;
#ifndef GENIE_PRE_R3
} else if (gheprec->Summary()->ProcInfo().IsDiffractive()) {
if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return 35;
else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg()))
return -35;
#endif
// NC OTHER
} else {
return utils::ghep::NeutReactionCode(gheprec);
}
}
return 0;
}
void GENIEInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Check for GENIE Event
if (!fGenieNtpl)
return;
if (!fGenieNtpl->event)
return;
// Cast Event Record
fGenieGHep = static_cast(fGenieNtpl->event);
if (!fGenieGHep)
return;
// Convert GENIE Reaction Code
fNUISANCEEvent->Mode = ConvertGENIEReactionCode(fGenieGHep);
if (!fNUISANCEEvent->Mode) {
NUIS_ERR(WRN, "Failed to determine mode for GENIE event: ");
std::cout << *fGenieGHep << std::endl;
}
// Set Event Info
fNUISANCEEvent->fEventNo = 0.0;
fNUISANCEEvent->fTotCrs = fGenieGHep->XSec();
+
// Have a bool storing if interaction happened on free or bound nucleon
bool IsFree = false;
+
// Set the TargetPDG
if (fGenieGHep->TargetNucleus() != NULL) {
fNUISANCEEvent->fTargetPDG = fGenieGHep->TargetNucleus()->Pdg();
IsFree = false;
// Sometimes GENIE scatters off free nucleons, electrons, photons
// In which TargetNucleus is NULL and we need to find the initial state
// particle
} else {
// Check the particle is an initial state particle
// Follows GHepRecord::TargetNucleusPosition but doesn't do check on
// pdg::IsIon
GHepParticle *p = fGenieGHep->Particle(1);
// Check that particle 1 actually exists
if (!p) {
NUIS_ABORT("Can't find particle 1 for GHepRecord");
}
// If not an ion but is an initial state particle
if (!pdg::IsIon(p->Pdg()) && p->Status() == kIStInitialState) {
IsFree = true;
fNUISANCEEvent->fTargetPDG = p->Pdg();
// Catch if something strange happens:
// Here particle 1 is not an initial state particle OR
// particle 1 is an ion OR
// both
} else {
if (pdg::IsIon(p->Pdg())) {
NUIS_ABORT(
"Particle 1 in GHepRecord stack is an ion but isn't an initial "
"state particle");
} else {
NUIS_ABORT(
"Particle 1 in GHepRecord stack is not an ion but is an initial "
"state particle");
}
}
}
// Set the A and Z and H from the target PDG
// Depends on if we scattered off a free or bound nucleon
if (!IsFree) {
fNUISANCEEvent->fTargetA =
- TargetUtils::GetTargetAFromPDG(fNUISANCEEvent->fTargetPDG);
+ TargetUtils::GetTargetAFromPDG(fNUISANCEEvent->fTargetPDG);
fNUISANCEEvent->fTargetZ =
- TargetUtils::GetTargetZFromPDG(fNUISANCEEvent->fTargetPDG);
+ TargetUtils::GetTargetZFromPDG(fNUISANCEEvent->fTargetPDG);
fNUISANCEEvent->fTargetH = 0;
} else {
// If free proton scattering
if (fNUISANCEEvent->fTargetPDG == 2212) {
fNUISANCEEvent->fTargetA = 1;
fNUISANCEEvent->fTargetZ = 1;
fNUISANCEEvent->fTargetH = 1;
// If free neutron scattering
} else if (fNUISANCEEvent->fTargetPDG == 2112) {
fNUISANCEEvent->fTargetA = 0;
fNUISANCEEvent->fTargetZ = 1;
fNUISANCEEvent->fTargetH = 0;
// If neither
} else {
fNUISANCEEvent->fTargetA = 0;
fNUISANCEEvent->fTargetZ = 0;
fNUISANCEEvent->fTargetH = 0;
}
}
fNUISANCEEvent->fBound = !IsFree;
fNUISANCEEvent->InputWeight =
- 1.0; //(1E+38 / genie::units::cm2) * fGenieGHep->XSec();
-
- // And the custom weights
- if (fNOvAWeights) {
- fNUISANCEEvent->CustomWeight = NOVAw;
- fNUISANCEEvent->CustomWeightArray[0] = MAQEw;
- fNUISANCEEvent->CustomWeightArray[1] = NonResw;
- fNUISANCEEvent->CustomWeightArray[2] = RPAQEw;
- fNUISANCEEvent->CustomWeightArray[3] = RPARESw;
- fNUISANCEEvent->CustomWeightArray[4] = MECw;
- fNUISANCEEvent->CustomWeightArray[5] = NOVAw;
- } else {
- fNUISANCEEvent->CustomWeight = 1.0;
- fNUISANCEEvent->CustomWeightArray[0] = 1.0;
- fNUISANCEEvent->CustomWeightArray[1] = 1.0;
- fNUISANCEEvent->CustomWeightArray[2] = 1.0;
- fNUISANCEEvent->CustomWeightArray[3] = 1.0;
- fNUISANCEEvent->CustomWeightArray[4] = 1.0;
- fNUISANCEEvent->CustomWeightArray[5] = 1.0;
- }
+ 1.0; //(1E+38 / genie::units::cm2) * fGenieGHep->XSec();
// Get N Particle Stack
unsigned int npart = fGenieGHep->GetEntries();
unsigned int kmax = fNUISANCEEvent->kMaxParticles;
if (npart > kmax) {
NUIS_ERR(WRN, "GENIE has too many particles, expanding stack.");
fNUISANCEEvent->ExpandParticleStack(npart);
}
// Fill Particle Stack
GHepParticle *p = 0;
TObjArrayIter iter(fGenieGHep);
fNUISANCEEvent->fNParticles = 0;
// Loop over all particles
while ((p = (dynamic_cast((iter).Next())))) {
if (!p)
continue;
// Get Status
int state = GetGENIEParticleStatus(p, fNUISANCEEvent->Mode);
// Remove Undefined
if (kRemoveUndefParticles && state == kUndefinedState)
continue;
// Remove FSI
if (kRemoveFSIParticles && state == kFSIState)
continue;
if (kRemoveNuclearParticles &&
(state == kNuclearInitial || state == kNuclearRemnant))
continue;
// Fill Vectors
int curpart = fNUISANCEEvent->fNParticles;
fNUISANCEEvent->fParticleState[curpart] = state;
// Mom
fNUISANCEEvent->fParticleMom[curpart][0] = p->Px() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][1] = p->Py() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][2] = p->Pz() * 1.E3;
fNUISANCEEvent->fParticleMom[curpart][3] = p->E() * 1.E3;
// PDG
fNUISANCEEvent->fParticlePDG[curpart] = p->Pdg();
// Set if the particle was on the fundamental vertex
- fNUISANCEEvent->fPrimaryVertex[curpart] = (p->FirstMother() < 2);
+ fNUISANCEEvent->fPrimaryVertex[curpart] = IsPrimary(p);
// Add to N particle count
fNUISANCEEvent->fNParticles++;
// Extra Check incase GENIE fails.
if ((UInt_t)fNUISANCEEvent->fNParticles == kmax) {
NUIS_ERR(WRN, "Number of GENIE Particles exceeds maximum!");
NUIS_ERR(WRN, "Extend kMax, or run without including FSI particles!");
break;
}
}
// Fill Extra Stack
if (fSaveExtra)
fGenieInfo->FillGeneratorInfo(fGenieNtpl);
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
FitParticle *ISAnyLepton = fNUISANCEEvent->GetHMISAnyLeptons();
if (ISAnyLepton) {
fNUISANCEEvent->probe_E = ISAnyLepton->E();
fNUISANCEEvent->probe_pdg = ISAnyLepton->PDG();
}
return;
}
void GENIEInputHandler::Print() {}
#endif
diff --git a/src/InputHandler/GENIEInputHandler.h b/src/InputHandler/GENIEInputHandler.h
index a09cda0..44d5d1a 100644
--- a/src/InputHandler/GENIEInputHandler.h
+++ b/src/InputHandler/GENIEInputHandler.h
@@ -1,130 +1,121 @@
// Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef GENIEINPUTHANDLER_H
#define GENIEINPUTHANDLER_H
/*!
* \addtogroup InputHandler
* @{
*/
#ifdef __GENIE_ENABLED__
#include "InputHandler.h"
#include "InputUtils.h"
#include "PlotUtils.h"
#ifdef GENIE_PRE_R3
#include "GHEP/GHepParticle.h"
#include "PDG/PDGUtils.h"
#include "GHEP/GHepUtils.h"
#include "Conventions/Units.h"
#include "EVGCore/EventRecord.h"
#include "GHEP/GHepRecord.h"
#include "Ntuple/NtpMCEventRecord.h"
#else
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Framework/GHEP/GHepUtils.h"
#include "Framework/Conventions/Units.h"
#include "Framework/EventGen/EventRecord.h"
#include "Framework/GHEP/GHepRecord.h"
#include "Framework/Ntuple/NtpMCEventRecord.h"
#endif
using namespace genie;
/// GENIE Generator Container to save extra particle status codes.
class GENIEGeneratorInfo : public GeneratorInfoBase {
-public:
- GENIEGeneratorInfo() {};
- virtual ~GENIEGeneratorInfo();
+ public:
+ GENIEGeneratorInfo() {};
+ virtual ~GENIEGeneratorInfo();
- /// Assigns information to branches
- void AddBranchesToTree(TTree* tn);
+ /// Assigns information to branches
+ void AddBranchesToTree(TTree* tn);
- /// Setup reading information from branches
- void SetBranchesFromTree(TTree* tn);
+ /// Setup reading information from branches
+ void SetBranchesFromTree(TTree* tn);
- /// Allocate any dynamic arrays for a new particle stack size
- void AllocateParticleStack(int stacksize);
+ /// Allocate any dynamic arrays for a new particle stack size
+ void AllocateParticleStack(int stacksize);
- /// Clear any dynamic arrays
- void DeallocateParticleStack();
+ /// Clear any dynamic arrays
+ void DeallocateParticleStack();
- /// Read extra genie information from the event
- void FillGeneratorInfo(NtpMCEventRecord* ntpl);
+ /// Read extra genie information from the event
+ void FillGeneratorInfo(NtpMCEventRecord* ntpl);
- /// Reset extra information to default/empty values
- void Reset();
+ /// Reset extra information to default/empty values
+ void Reset();
- int kMaxParticles; ///< Number of particles in stack
- int* fGenieParticlePDGs; ///< GENIE Particle PDGs (example)
+ int kMaxParticles; ///< Number of particles in stack
+ int* fGenieParticlePDGs; ///< GENIE Particle PDGs (example)
};
/// Main GENIE InputHandler
class GENIEInputHandler : public InputHandlerBase {
-public:
+ public:
- /// Standard constructor given a name and input files
- GENIEInputHandler(std::string const& handle, std::string const& rawinputs);
- virtual ~GENIEInputHandler();
+ /// Standard constructor given a name and input files
+ GENIEInputHandler(std::string const& handle, std::string const& rawinputs);
+ virtual ~GENIEInputHandler();
- /// Create a TTree Cache to speed up file read
- void CreateCache();
+ /// Create a TTree Cache to speed up file read
+ void CreateCache();
- /// Remove TTree Cache to save memory
- void RemoveCache();
-
- /// Returns a NUISANCE format event from the GENIE TTree. If !lightweight
- /// then CalcNUISANCEKinematics() is called to convert the GENIE event into
- /// a standard NUISANCE format.
- FitEvent* GetNuisanceEvent(const UInt_t entry, const bool lightweight = false);
+ /// Remove TTree Cache to save memory
+ void RemoveCache();
- /// Converts GENIE event into standard NUISANCE FitEvent by looping over all
- /// particles in the event and adding them to stack in fNUISANCEEvent.
- void CalcNUISANCEKinematics();
+ /// Returns a NUISANCE format event from the GENIE TTree. If !lightweight
+ /// then CalcNUISANCEKinematics() is called to convert the GENIE event into
+ /// a standard NUISANCE format.
+ FitEvent* GetNuisanceEvent(const UInt_t entry, const bool lightweight = false);
- /// Placeholder for GENIE related event printing.
- void Print();
+ /// Converts GENIE event into standard NUISANCE FitEvent by looping over all
+ /// particles in the event and adding them to stack in fNUISANCEEvent.
+ void CalcNUISANCEKinematics();
- /// Converts GENIE particle status codes into NUISANCE status codes.
- int GetGENIEParticleStatus(genie::GHepParticle* part, int mode = 0);
+ /// Placeholder for GENIE related event printing.
+ void Print();
- /// Converts GENIE event reaction codes into NUISANCE reaction codes.
- int ConvertGENIEReactionCode(GHepRecord* gheprec);
+ /// Converts GENIE particle status codes into NUISANCE status codes.
+ int GetGENIEParticleStatus(genie::GHepParticle* part, int mode = 0);
- GHepRecord* fGenieGHep; ///< Pointer to actual event record
- NtpMCEventRecord* fGenieNtpl; ///< Ntpl Wrapper Class
+ /// Converts GENIE event reaction codes into NUISANCE reaction codes.
+ int ConvertGENIEReactionCode(GHepRecord* gheprec);
- TChain* fGENIETree; ///< Main GENIE Event TTree
- bool fSaveExtra; ///< Flag to save Extra GENIE info into Nuisance Event
- GENIEGeneratorInfo* fGenieInfo; ///< Extra GENIE Generator Info Writer
+ GHepRecord* fGenieGHep; ///< Pointer to actual event record
+ NtpMCEventRecord* fGenieNtpl; ///< Ntpl Wrapper Class
- bool fNOvAWeights; ///< Flag to save nova weights or not
+ TChain* fGENIETree; ///< Main GENIE Event TTree
+ bool fSaveExtra; ///< Flag to save Extra GENIE info into Nuisance Event
+ GENIEGeneratorInfo* fGenieInfo; ///< Extra GENIE Generator Info Writer
- // Extra weights from Jeremy for NOvA weights
- double MAQEw;
- double NonResw;
- double RPAQEw;
- double RPARESw;
- double MECw;
- double DISw;
- double NOVAw;
+ bool IsPrimary(GHepParticle *p);
};
/*! @} */
#endif
#endif
diff --git a/src/Reweight/FitWeight.cxx b/src/Reweight/FitWeight.cxx
index be3b57d..dfbc55d 100644
--- a/src/Reweight/FitWeight.cxx
+++ b/src/Reweight/FitWeight.cxx
@@ -1,305 +1,314 @@
#include "FitWeight.h"
#include "GENIEWeightEngine.h"
#include "LikelihoodWeightEngine.h"
#include "ModeNormEngine.h"
#include "NEUTWeightEngine.h"
#include "NIWGWeightEngine.h"
#include "NUISANCEWeightEngine.h"
#include "NuWroWeightEngine.h"
#include "OscWeightEngine.h"
#include "SampleNormEngine.h"
#include "SplineWeightEngine.h"
#include "T2KWeightEngine.h"
#ifdef __NOVA_ENABLED__
#include "NOvARwgtEngine.h"
#endif
#ifdef __NUSYST_ENABLED__
#include "nusystematicsWeightEngine.h"
#endif
void FitWeight::AddRWEngine(int type) {
NUIS_LOG(FIT, "Adding reweight engine " << type);
switch (type) {
case kNEUT:
fAllRW[type] = new NEUTWeightEngine("neutrw");
break;
case kNUWRO:
fAllRW[type] = new NuWroWeightEngine("nuwrorw");
break;
case kGENIE:
fAllRW[type] = new GENIEWeightEngine("genierw");
break;
case kNORM:
fAllRW[type] = new SampleNormEngine("normrw");
break;
case kLIKEWEIGHT:
fAllRW[type] = new LikelihoodWeightEngine("likerw");
break;
case kT2K:
fAllRW[type] = new T2KWeightEngine("t2krw");
break;
case kCUSTOM:
fAllRW[type] = new NUISANCEWeightEngine("nuisrw");
break;
case kSPLINEPARAMETER:
fAllRW[type] = new SplineWeightEngine("splinerw");
break;
case kNIWG:
fAllRW[type] = new NIWGWeightEngine("niwgrw");
break;
case kOSCILLATION:
fAllRW[type] = new OscWeightEngine();
break;
case kMODENORM:
fAllRW[type] = new ModeNormEngine();
break;
#ifdef __NOVA_ENABLED__
case kNOvARWGT:
fAllRW[type] = new NOvARwgtEngine();
break;
#endif
#ifdef __NUSYST_ENABLED__
case kNuSystematics:
fAllRW[type] = new nusystematicsWeightEngine();
break;
#endif
default:
NUIS_ABORT("CANNOT ADD RW Engine for unknown dial type: " << type);
break;
}
}
WeightEngineBase *FitWeight::GetRWEngine(int type) {
if (HasRWEngine(type)) {
return fAllRW[type];
}
NUIS_ABORT("CANNOT get RW Engine for dial type: " << type);
}
bool FitWeight::HasRWEngine(int type) {
switch (type) {
case kNEUT:
case kNUWRO:
case kGENIE:
case kNORM:
case kLIKEWEIGHT:
case kT2K:
case kCUSTOM:
case kSPLINEPARAMETER:
case kNIWG:
case kOSCILLATION:
#ifdef __NOVA_ENABLED__
case kNOvARWGT:
#endif
#ifdef __NUSYST_ENABLED__
case kNuSystematics:
#endif
{
return fAllRW.count(type);
}
default: { NUIS_ABORT("CANNOT get RW Engine for dial type: " << type); }
}
}
void FitWeight::IncludeDial(std::string name, std::string type, double val) {
// Should register the dial here.
int typeenum = Reweight::ConvDialType(type);
IncludeDial(name, typeenum, val);
}
void FitWeight::IncludeDial(std::string name, int dialtype, double val) {
// Get the dial enum
int nuisenum = Reweight::ConvDial(name, dialtype);
if (nuisenum == -1) {
NUIS_ERR(FTL, "Could not include dial " << name);
NUIS_ERR(FTL, "With dialtype: " << dialtype);
NUIS_ERR(FTL, "With value: " << val);
NUIS_ABORT("With nuisenum: " << nuisenum);
}
// Setup RW Engine Pointer
if (fAllRW.find(dialtype) == fAllRW.end()) {
AddRWEngine(dialtype);
}
WeightEngineBase *rw = fAllRW[dialtype];
// Include the dial
rw->IncludeDial(name, val);
// Set Dial Value
if (val != -9999.9) {
rw->SetDialValue(name, val);
}
// Sort Maps
fAllEnums[name] = nuisenum;
fAllValues[nuisenum] = val;
// Sort Lists
fNameList.push_back(name);
fEnumList.push_back(nuisenum);
fValueList.push_back(val);
}
void FitWeight::Reconfigure(bool silent) {
// Reconfigure all added RW engines
for (std::map::iterator iter = fAllRW.begin();
iter != fAllRW.end(); iter++) {
(*iter).second->Reconfigure(silent);
}
}
void FitWeight::SetDialValue(std::string name, double val) {
// Add extra check, if name not found look for one with name in it.
int nuisenum = fAllEnums[name];
SetDialValue(nuisenum, val);
}
// Allow for name aswell using GlobalList to determine sample name.
void FitWeight::SetDialValue(int nuisenum, double val) {
// Conv dial type
int dialtype = Reweight::GetDialType(nuisenum);
if (fAllRW.find(dialtype) == fAllRW.end()) {
- NUIS_ERR(FTL, "Can't find RW engine for parameter " << fNameList[dialtype]);
+
+ std::string name = "";
+ for(size_t i = 0; i < fEnumList.size(); ++i){
+ if(fEnumList[i] == nuisenum){
+ name = fNameList[i];
+ break;
+ }
+ }
+
+ NUIS_ERR(FTL, "Can't find RW engine for parameter " << name);
NUIS_ERR(FTL, "With dialtype " << dialtype << ", "
<< Reweight::RemoveDialType(nuisenum));
NUIS_ABORT("Are you sure you enabled the right engines?");
}
// Get RW Engine for this dial
fAllRW[dialtype]->SetDialValue(nuisenum, val);
fAllValues[nuisenum] = val;
// Update ValueList
for (size_t i = 0; i < fEnumList.size(); i++) {
if (fEnumList[i] == nuisenum) {
fValueList[i] = val;
}
}
}
void FitWeight::SetAllDials(const double *x, int n) {
for (size_t i = 0; i < (UInt_t)n; i++) {
int rwenum = fEnumList[i];
SetDialValue(rwenum, x[i]);
}
Reconfigure();
}
double FitWeight::GetDialValue(std::string name) {
// Add extra check, if name not found look for one with name in it.
int nuisenum = fAllEnums[name];
return GetDialValue(nuisenum);
}
double FitWeight::GetDialValue(int nuisenum) { return fAllValues[nuisenum]; }
int FitWeight::GetDialPos(std::string name) {
int rwenum = fAllEnums[name];
return GetDialPos(rwenum);
}
int FitWeight::GetDialPos(int nuisenum) {
for (size_t i = 0; i < fEnumList.size(); i++) {
if (fEnumList[i] == nuisenum) {
return i;
}
}
NUIS_ABORT("No Dial Found! (enum = " << nuisenum << ") ");
}
bool FitWeight::DialIncluded(std::string name) {
return (fAllEnums.find(name) != fAllEnums.end());
}
bool FitWeight::DialIncluded(int rwenum) {
return (fAllValues.find(rwenum) != fAllValues.end());
}
double FitWeight::CalcWeight(BaseFitEvt *evt) {
double rwweight = 1.0;
for (std::map::iterator iter = fAllRW.begin();
iter != fAllRW.end(); iter++) {
double w = (*iter).second->CalcWeight(evt);
rwweight *= w;
}
return rwweight;
}
void FitWeight::UpdateWeightEngine(const double *x) {
size_t count = 0;
for (std::vector::iterator iter = fEnumList.begin();
iter != fEnumList.end(); iter++) {
SetDialValue((*iter), x[count]);
count++;
}
}
void FitWeight::GetAllDials(double *x, int n) {
for (int i = 0; i < n; i++) {
x[i] = GetDialValue(fEnumList[i]);
}
}
// bool FitWeight::NeedsEventReWeight(const double* x) {
// bool haschange = false;
// size_t count = 0;
// // Compare old to new and decide if RW needed.
// for (std::vector::iterator iter = fEnumList.begin();
// iter != fEnumList.end(); iter++) {
// int nuisenum = (*iter);
-// int type = (nuisenum / 1000) - (nuisenum % 1000);
+// int type = (nuisenum / NUIS_DIAL_OFFSET) - (nuisenum % NUIS_DIAL_OFFSET);
// // Compare old to new
// double oldval = GetDialValue(nuisenum);
// double newval = x[count];
// if (oldval != newval) {
// if (fAllRW[type]->NeedsEventReWeight()) {
// haschange = true;
// }
// }
// count++;
// }
// return haschange;
// }
double FitWeight::GetSampleNorm(std::string name) {
if (name.empty())
return 1.0;
// Find norm dial
if (fAllEnums.find(name + "_norm") != fAllEnums.end()) {
if (fAllValues.find(fAllEnums[name + "_norm"]) != fAllValues.end()) {
return fAllValues[fAllEnums[name + "_norm"]];
} else {
return 1.0;
}
} else {
return 1.0;
}
}
void FitWeight::Print() {
NUIS_LOG(REC, "Fit Weight State: ");
for (size_t i = 0; i < fNameList.size(); i++) {
NUIS_LOG(REC,
"|-> Par " << i << ". " << fNameList[i] << " " << fValueList[i]);
}
}
diff --git a/src/Reweight/MINERvAWeightCalcs.cxx b/src/Reweight/MINERvAWeightCalcs.cxx
index 0f555b4..68f7aa7 100644
--- a/src/Reweight/MINERvAWeightCalcs.cxx
+++ b/src/Reweight/MINERvAWeightCalcs.cxx
@@ -1,871 +1,871 @@
#ifdef __MINERVA_RW_ENABLED__
#ifdef __GENIE_ENABLED__
#include "MINERvAWeightCalcs.h"
#include "BaseFitEvt.h"
namespace nuisance {
namespace reweight {
//*******************************************************
MINERvAReWeight_QE::MINERvAReWeight_QE() {
//*******************************************************
fTwk_NormCCQE = 0.0;
fDef_NormCCQE = 1.0;
fCur_NormCCQE = fDef_NormCCQE;
}
//*******************************************************
MINERvAReWeight_QE::~MINERvAReWeight_QE(){};
//*******************************************************
//*******************************************************
double MINERvAReWeight_QE::CalcWeight(BaseFitEvt *evt) {
//*******************************************************
// Check GENIE
if (evt->fType != kGENIE) return 1.0;
// Extract the GENIE Record
GHepRecord *ghep = static_cast(evt->genie_event->event);
const Interaction *interaction = ghep->Summary();
// const InitialState& init_state = interaction->InitState();
const ProcessInfo &proc_info = interaction->ProcInfo();
// const Target& tgt = init_state.Tgt();
// If the event is not QE this Calc doesn't handle it
if (!proc_info.IsQuasiElastic()) return 1.0;
// WEIGHT CALCULATIONS -------------
double w = 1.0;
// CCQE Dial
if (!proc_info.IsWeakCC()) w *= fCur_NormCCQE;
// Return Combined Weight
return w;
}
//*******************************************************
void MINERvAReWeight_QE::SetDialValue(std::string name, double val) {
//*******************************************************
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//*******************************************************
void MINERvAReWeight_QE::SetDialValue(int rwenum, double val) {
//*******************************************************
// Check Handled
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
if (!IsHandled(curenum)) return;
// Set Values
if (curenum == Reweight::kMINERvARW_NormCCQE) {
fTwk_NormCCQE = val;
fCur_NormCCQE = fDef_NormCCQE + fTwk_NormCCQE;
}
// Define Tweaked
fTweaked = ((fTwk_NormCCQE != 0.0));
}
//*******************************************************
bool MINERvAReWeight_QE::IsHandled(int rwenum) {
//*******************************************************
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::kMINERvARW_NormCCQE:
return true;
default:
return false;
}
}
//*******************************************************
MINERvAReWeight_MEC::MINERvAReWeight_MEC() {
//*******************************************************
fTwk_NormCCMEC = 0.0;
fDef_NormCCMEC = 1.0;
fCur_NormCCMEC = fDef_NormCCMEC;
}
//*******************************************************
MINERvAReWeight_MEC::~MINERvAReWeight_MEC(){};
//*******************************************************
//*******************************************************
double MINERvAReWeight_MEC::CalcWeight(BaseFitEvt *evt) {
//*******************************************************
// Check GENIE
if (evt->fType != kGENIE)
return 1.0;
// Extract the GENIE Record
GHepRecord *ghep = static_cast(evt->genie_event->event);
const Interaction *interaction = ghep->Summary();
// const InitialState& init_state = interaction->InitState();
const ProcessInfo &proc_info = interaction->ProcInfo();
// const Target& tgt = init_state.Tgt();
// If the event is not MEC this Calc doesn't handle it
if (!proc_info.IsMEC()) return 1.0;
// WEIGHT CALCULATIONS -------------
double w = 1.0;
// CCMEC Dial
if (!proc_info.IsWeakCC()) w *= fCur_NormCCMEC;
// Return Combined Weight
return w;
}
//*******************************************************
void MINERvAReWeight_MEC::SetDialValue(std::string name, double val) {
//*******************************************************
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//*******************************************************
void MINERvAReWeight_MEC::SetDialValue(int rwenum, double val) {
//*******************************************************
// Check Handled
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
if (!IsHandled(curenum))
return;
// Set Values
if (curenum == Reweight::kMINERvARW_NormCCMEC) {
fTwk_NormCCMEC = val;
fCur_NormCCMEC = fDef_NormCCMEC + fTwk_NormCCMEC;
}
// Define Tweaked
fTweaked = ((fTwk_NormCCMEC != 0.0));
}
//*******************************************************
bool MINERvAReWeight_MEC::IsHandled(int rwenum) {
//*******************************************************
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::kMINERvARW_NormCCMEC:
return true;
default:
return false;
}
}
//*******************************************************
MINERvAReWeight_RES::MINERvAReWeight_RES() {
//*******************************************************
fTwk_NormCCRES = 0.0;
fDef_NormCCRES = 1.0;
fCur_NormCCRES = fDef_NormCCRES;
}
//*******************************************************
MINERvAReWeight_RES::~MINERvAReWeight_RES(){};
//*******************************************************
//*******************************************************
double MINERvAReWeight_RES::CalcWeight(BaseFitEvt *evt) {
//*******************************************************
// std::cout << "Caculating RES" << std::endl;
// Check GENIE
if (evt->fType != kGENIE) return 1.0;
// Extract the GENIE Record
GHepRecord *ghep = static_cast(evt->genie_event->event);
const Interaction *interaction = ghep->Summary();
// const InitialState& init_state = interaction->InitState();
const ProcessInfo &proc_info = interaction->ProcInfo();
// const Target& tgt = init_state.Tgt();
// If the event is not RES this Calc doesn't handle it
if (!proc_info.IsResonant()) return 1.0;
// WEIGHT CALCULATIONS -------------
double w = 1.0;
// CCRES Dial
if (proc_info.IsWeakCC()) w *= fCur_NormCCRES;
// Return Combined Weight
return w;
}
//*******************************************************
void MINERvAReWeight_RES::SetDialValue(std::string name, double val) {
//*******************************************************
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//*******************************************************
void MINERvAReWeight_RES::SetDialValue(int rwenum, double val) {
//*******************************************************
// Check Handled
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
if (!IsHandled(curenum)) return;
// Set Values
if (curenum == Reweight::kMINERvARW_NormCCRES) {
fTwk_NormCCRES = val;
fCur_NormCCRES = fDef_NormCCRES + fTwk_NormCCRES;
}
// Define Tweaked
fTweaked = ((fTwk_NormCCRES != 0.0));
}
//*******************************************************
bool MINERvAReWeight_RES::IsHandled(int rwenum) {
//*******************************************************
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::kMINERvARW_NormCCRES:
return true;
default:
return false;
}
}
//*******************************************************
RikRPA::RikRPA() {
//*******************************************************
// - Syst : kMINERvA_RikRPA_ApplyRPA
// - Type : Binary
// - Limits : 0.0 (false) -> 1.0 (true)
// - Default : 0.0
fApplyDial_RPACorrection = false;
// - Syst : kMINERvA_RikRPA_LowQ2
// - Type : Absolute
// - Limits : 1.0 -> 1.0
// - Default : 0.0
// - Frac Error : 100%
fDefDial_RPALowQ2 = 0.0;
fCurDial_RPALowQ2 = fDefDial_RPALowQ2;
fErrDial_RPALowQ2 = 0.0;
// - Syst : kMINERvA_RikRPA_HighQ2
// - Type : Absolute
// - Limits : 1.0 -> 1.0
// - Default : 0.0
// - Frac Error : 100%
fDefDial_RPAHighQ2 = 0.0;
fCurDial_RPAHighQ2 = fDefDial_RPAHighQ2;
fErrDial_RPAHighQ2 = 1.0;
// - Syst : kMINERvA_RikRESRPA_ApplyRPA
// - Type : Binary
// - Limits : 0.0 (false) -> 1.0 (true)
// - Default : 0.0
fApplyDial_RESRPACorrection = false;
// - Syst : kMINERvA_RikRESRPA_LowQ2
// - Type : Absolute
// - Limits : 1.0 -> 1.0
// - Default : 0.0
// - Frac Error : 100%
fDefDial_RESRPALowQ2 = 0.0;
fCurDial_RESRPALowQ2 = fDefDial_RESRPALowQ2;
fErrDial_RESRPALowQ2 = 0.0;
// - Syst : kMINERvA_RikRESRPA_HighQ2
// - Type : Absolute
// - Limits : 1.0 -> 1.0
// - Default : 0.0
// - Frac Error : 100%
fDefDial_RESRPAHighQ2 = 0.0;
fCurDial_RESRPAHighQ2 = fDefDial_RESRPAHighQ2;
fErrDial_RESRPAHighQ2 = 1.0;
// Setup Calculators
fEventWeights = new double[5];
for (int i = 0; i < kMaxCalculators; i++) {
fRPACalculators[i] = NULL;
}
fTweaked = false;
}
//*******************************************************
RikRPA::~RikRPA() {
//*******************************************************
// delete fEventWeights;
// for (size_t i = 0; i < kMaxCalculators; i++) {
// if (fRPACalculators[i]) delete fRPACalculators[i];
// fRPACalculators[i] = NULL;
// }
}
//*******************************************************
double RikRPA::CalcWeight(BaseFitEvt *evt) {
//*******************************************************
// LOG(FIT) << "Calculating RikRPA" << std::endl;
// Return 1.0 if not tweaked
if (!fTweaked)
return 1.0;
double w = 1.0;
// Extract the GENIE Record
GHepRecord *ghep = static_cast(evt->genie_event->event);
const Interaction *interaction = ghep->Summary();
const InitialState &init_state = interaction->InitState();
const ProcessInfo &proc_info = interaction->ProcInfo();
// const Kinematics & kine = interaction->Kine();
// const XclsTag & xcls = interaction->ExclTag();
const Target &tgt = init_state.Tgt();
// If not QE return 1.0
// LOG(FIT) << "RikRPA : Event QE = " << proc_info.IsQuasiElastic() <<
// std::endl;
if (!tgt.IsNucleus()) {
return 1.0;
}
if (!proc_info.IsQuasiElastic() && !proc_info.IsResonant())
return 1.0;
// Extract Beam and Target PDG
GHepParticle *neutrino = ghep->Probe();
int bpdg = neutrino->Pdg();
GHepParticle *target = ghep->TargetNucleus();
assert(target);
int tpdg = target->Pdg();
// Find the enum we need
int calcenum = GetRPACalcEnum(bpdg, tpdg);
if (calcenum == -1)
return 1.0;
// Check we have the RPA Calc setup for this enum
// if not, set it up at that point
if (!fRPACalculators[calcenum])
SetupRPACalculator(calcenum);
weightRPA *rpacalc = fRPACalculators[calcenum];
if (!rpacalc) {
NUIS_ABORT("Failed to grab the RPA Calculator : " << calcenum);
}
// Extract Q0-Q3
GHepParticle *fsl = ghep->FinalStatePrimaryLepton();
const TLorentzVector &k1 = *(neutrino->P4());
const TLorentzVector &k2 = *(fsl->P4());
double q0 = fabs((k1 - k2).E());
double q3 = fabs((k1 - k2).Vect().Mag());
double Q2 = fabs((k1 - k2).Mag2());
// Quasielastic
if (proc_info.IsQuasiElastic()) {
// Now use q0-q3 and RPA Calculator to fill fWeights
rpacalc->getWeight(q0, q3, fEventWeights);
if (fApplyDial_RPACorrection) {
w *= fEventWeights[0]; // CV
}
// Syst Application : kMINERvA_RikRPA_LowQ2
if (fabs(fCurDial_RPALowQ2) > 0.0) {
double interpw = fEventWeights[0];
if (fCurDial_RPALowQ2 > 0.0 && Q2 < 2.0) {
interpw =
fEventWeights[0] - (fEventWeights[0] - fEventWeights[1]) *
fCurDial_RPALowQ2; // WLow+
} else if (fCurDial_RPALowQ2 < 0.0 && Q2 < 2.0) {
interpw = fEventWeights[0] - (fEventWeights[2] - fEventWeights[0]) *
fCurDial_RPALowQ2; // WLow-
}
w *= interpw / fEventWeights[0]; // Div by CV again
}
// Syst Application : kMINERvA_RikRPA_HighQ2
if (fabs(fCurDial_RPAHighQ2) > 0.0) {
double interpw = fEventWeights[0];
if (fCurDial_RPAHighQ2 > 0.0) {
interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) *
fCurDial_RPAHighQ2; // WHigh+
} else if (fCurDial_RPAHighQ2 < 0.0) {
interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) *
fCurDial_RPAHighQ2; // WHigh-
}
w *= interpw / fEventWeights[0]; // Div by CV again
}
}
// Resonant Events
if (proc_info.IsResonant()) {
// Now use Q2 and RESRPA Calculator to fill fWeights
double CV = rpacalc->getWeight(Q2);
if (fApplyDial_RESRPACorrection) {
w *= CV; // fEventWeights[0]; // CVa
}
/*
// Syst Application : kMINERvA_RikRESRPA_LowQ2
if (fabs(fCurDial_RESRPAHighQ2) > 0.0) {
double interpw = fEventWeights[0];
if (fCurDial_RESRPAHighQ2 > 0.0) {
interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) *
fCurDial_RESRPAHighQ2; // WHigh+
} else if (fCurDial_RESRPAHighQ2 < 0.0) {
interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) *
fCurDial_RESRPAHighQ2; // WHigh-
}
w *= interpw / fEventWeights[0]; // Div by CV again
}
// Syst Application : kMINERvA_RikRESRPA_HighQ2
if (fabs(fCurDial_RESRPAHighQ2) > 0.0) {
double interpw = fEventWeights[0];
if (fCurDial_RESRPAHighQ2 > 0.0) {
interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) *
fCurDial_RESRPAHighQ2; // WHigh+
} else if (fCurDial_RESRPAHighQ2 < 0.0) {
interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) *
fCurDial_RESRPAHighQ2; // WHigh-
}
w *= interpw / fEventWeights[0]; // Div by CV again
}
*/
}
// LOG(FIT) << "RPA Weight = " << w << std::endl;
return w;
} // namespace reweight
//*******************************************************
void RikRPA::SetDialValue(std::string name, double val) {
//*******************************************************
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//*******************************************************
void RikRPA::SetDialValue(int rwenum, double val) {
//*******************************************************
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
// Check Handled
if (!IsHandled(curenum))
return;
if (curenum == Reweight::kMINERvARW_RikRPA_ApplyRPA)
fApplyDial_RPACorrection = (val > 0.5);
if (curenum == Reweight::kMINERvARW_RikRPA_LowQ2)
fCurDial_RPALowQ2 = val;
if (curenum == Reweight::kMINERvARW_RikRPA_HighQ2)
fCurDial_RPAHighQ2 = val;
if (curenum == Reweight::kMINERvARW_RikRESRPA_ApplyRPA)
fApplyDial_RESRPACorrection = (val > 0.5);
if (curenum == Reweight::kMINERvARW_RikRESRPA_LowQ2)
fCurDial_RESRPALowQ2 = val;
if (curenum == Reweight::kMINERvARW_RikRESRPA_HighQ2)
fCurDial_RESRPAHighQ2 = val;
// Assign flag to say stuff has changed
fTweaked = (fApplyDial_RPACorrection ||
fabs(fCurDial_RPAHighQ2 - fDefDial_RPAHighQ2) > 0.0 ||
fabs(fCurDial_RPALowQ2 - fDefDial_RPALowQ2) > 0.0 ||
fApplyDial_RESRPACorrection ||
fabs(fCurDial_RESRPAHighQ2 - fDefDial_RESRPAHighQ2) > 0.0 ||
fabs(fCurDial_RESRPALowQ2 - fDefDial_RESRPALowQ2) > 0.0);
}
//*******************************************************
bool RikRPA::IsHandled(int rwenum) {
//*******************************************************
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::kMINERvARW_RikRESRPA_ApplyRPA:
return true;
case Reweight::kMINERvARW_RikRESRPA_LowQ2:
return true;
case Reweight::kMINERvARW_RikRESRPA_HighQ2:
return true;
case Reweight::kMINERvARW_RikRPA_ApplyRPA:
return true;
case Reweight::kMINERvARW_RikRPA_LowQ2:
return true;
case Reweight::kMINERvARW_RikRPA_HighQ2:
return true;
default:
return false;
}
}
//*******************************************************
void RikRPA::SetupRPACalculator(int calcenum) {
//*******************************************************
std::string rwdir = FitPar::GetDataBase() + "reweight/MINERvA/RikRPA/";
std::string fidir = "";
switch (calcenum) {
case kNuMuC12:
fidir = "outNievesRPAratio-nu12C-20GeV-20170202.root";
break;
case kNuMuO16:
fidir = "outNievesRPAratio-nu16O-20GeV-20170202.root";
break;
case kNuMuAr40:
fidir = "outNievesRPAratio-nu40Ar-20GeV-20170202.root";
break;
case kNuMuCa40:
fidir = "outNievesRPAratio-nu40Ca-20GeV-20170202.root";
break;
case kNuMuFe56:
fidir = "outNievesRPAratio-nu56Fe-20GeV-20170202.root";
break;
case kNuMuBarC12:
fidir = "outNievesRPAratio-anu12C-20GeV-20170202.root";
break;
case kNuMuBarO16:
fidir = "outNievesRPAratio-anu16O-20GeV-20170202.root";
break;
case kNuMuBarAr40:
fidir = "outNievesRPAratio-anu40Ar-20GeV-20170202.root";
break;
case kNuMuBarCa40:
fidir = "outNievesRPAratio-anu40Ca-20GeV-20170202.root";
break;
case kNuMuBarFe56:
fidir = "outNievesRPAratio-anu56Fe-20GeV-20170202.root";
break;
}
NUIS_LOG(FIT, "Loading RPA CALC : " << fidir);
TDirectory *olddir = gDirectory;
NUIS_LOG(FIT, "***********************************************");
NUIS_LOG(FIT, "Loading a new weightRPA calculator");
NUIS_LOG(FIT, "Authors: Rik Gran, Heidi Schellman");
NUIS_LOG(FIT, "Citation: arXiv:1705.02932 [hep-ex]");
NUIS_LOG(FIT, "***********************************************");
// Test the file exists
std::ifstream infile((rwdir + fidir).c_str());
if (!infile.good()) {
NUIS_ERR(FTL, "*** ERROR ***");
NUIS_ERR(FTL, "RikRPA file " << rwdir + fidir << " does not exist!");
NUIS_ERR(FTL,
"These can be found at https://nuisance.hepforge.org/files/RikRPA/");
NUIS_ERR(FTL,
"Please run: wget -r -nH --cut-dirs=2 -np -e robots=off -R "
"\"index.html*\" https://nuisance.hepforge.org/files/RikRPA/ -P "
<< rwdir);
NUIS_ERR(FTL, "And try again");
NUIS_ABORT("*************");
}
fRPACalculators[calcenum] = new weightRPA(rwdir + fidir);
olddir->cd();
return;
}
//*******************************************************
int RikRPA::GetRPACalcEnum(int bpdg, int tpdg) {
//*******************************************************
if (bpdg == 14 && tpdg == 1000060120)
return kNuMuC12;
else if (bpdg == 14 && tpdg == 1000080160)
return kNuMuO16;
else if (bpdg == 14 && tpdg == 1000180400)
return kNuMuAr40;
else if (bpdg == 14 && tpdg == 1000200400)
return kNuMuCa40;
else if (bpdg == 14 && tpdg == 1000280560)
return kNuMuFe56;
else if (bpdg == -14 && tpdg == 1000060120)
return kNuMuBarC12;
else if (bpdg == -14 && tpdg == 1000080160)
return kNuMuBarO16;
else if (bpdg == -14 && tpdg == 1000180400)
return kNuMuBarAr40;
else if (bpdg == -14 && tpdg == 1000200400)
return kNuMuBarCa40;
else if (bpdg == -14 && tpdg == 1000280560)
return kNuMuBarFe56;
else {
// NUIS_ERR(WRN, "Unknown beam and target combination for RPA Calcs! "
//<< bpdg << " " << tpdg);
}
return -1;
}
//*****************************************************************************
COHBrandon::COHBrandon() {
fApply_COHNorm = false;
fDef_COHNorm = 0.5;
fCur_COHNorm = fDef_COHNorm;
fTwk_COHNorm = 0.0;
fDef_COHCut = 0.450;
fCur_COHCut = fDef_COHCut;
fTwk_COHNorm = 0.0;
}
COHBrandon::~COHBrandon() {};
double COHBrandon::CalcWeight(BaseFitEvt* evt) {
// Check GENIE
if (evt->fType != kGENIE) return 1.0;
// Extract the GENIE Record
GHepRecord* ghep = static_cast(evt->genie_event->event);
const Interaction* interaction = ghep->Summary();
//const InitialState& init_state = interaction->InitState();
const ProcessInfo& proc_info = interaction->ProcInfo();
//const Target& tgt = init_state.Tgt();
// If the event is not QE this Calc doesn't handle it
if (!proc_info.IsCoherent()) return 1.0;
// WEIGHT CALCULATIONS -------------
double w = 1.0;
double pionE = -999.9;
TObjArrayIter piter(ghep);
GHepParticle * p = 0;
//int ip = -1;
while ( (p = (GHepParticle *) piter.Next()) ) {
int pdgc = p->Pdg();
int ist = p->Status();
// if (ghep->Particle(ip)->FirstMother()==0) continue;
if (pdg::IsPseudoParticle(p->Pdg())) continue;
if (ist == kIStStableFinalState) {
if (pdgc == 211 || pdgc == -211 || pdgc == 111) {
pionE = p->Energy();
break;
}
}
}
// std::cout << "Coherent Pion Energy : " << pionE << std::endl;
// Apply Weight
if (fApply_COHNorm && pionE > 0.0 && pionE < fCur_COHCut) {
w *= fCur_COHNorm;
}
// Return Combined Weight
return w;
}
void COHBrandon::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void COHBrandon::SetDialValue(int rwenum, double val) {
// Check Handled
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
if (!IsHandled(curenum)) return;
// Set Values
if (curenum == Reweight::kMINERvARW_NormCOH) {
fTwk_COHNorm = val;
fCur_COHNorm = fDef_COHNorm * (1.0 + 0.1 * fTwk_COHNorm);
}
if (curenum == Reweight::kMINERvARW_CutCOH) {
fTwk_COHCut = val;
fCur_COHCut = fDef_COHCut * (1.0 + 0.1 * fTwk_COHCut);
}
if (curenum == Reweight::kMINERvARW_ApplyCOH) {
fApply_COHNorm = (val > 0.5);
}
// Define Tweaked
fTweaked = (fApply_COHNorm);
}
bool COHBrandon::IsHandled(int rwenum) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::kMINERvARW_NormCOH:
case Reweight::kMINERvARW_CutCOH:
case Reweight::kMINERvARW_ApplyCOH:
return true;
default:
return false;
}
}
//*****************************************************************************
WEnhancement::WEnhancement() {
fApply_Enhancement = false;
fDef_WNorm = 1.0;
fCur_WNorm = fDef_WNorm;
fTwk_WNorm = 0.0;
fDef_WMean = 0.95;
fCur_WMean = fDef_WMean;
fTwk_WMean = 0.0;
fDef_WSigma = 0.225;
fCur_WSigma = fDef_WSigma;
fTwk_WSigma = 0.0;
}
WEnhancement::~WEnhancement() {};
double WEnhancement::CalcWeight(BaseFitEvt* evt) {
// Check GENIE
if (evt->fType != kGENIE) return 1.0;
// Extract the GENIE Record
GHepRecord* ghep = static_cast(evt->genie_event->event);
const Interaction* interaction = ghep->Summary();
//const Kinematics & kine = interaction->Kine();
//const InitialState& init_state = interaction->InitState();
const ProcessInfo& proc_info = interaction->ProcInfo();
//const Target& tgt = init_state.Tgt();
// If the event is not QE this Calc doesn't handle it
if (!proc_info.IsWeakCC()) return 1.0;
if (!proc_info.IsResonant()) return 1.0;
if (!fApply_Enhancement) return 1.0;
// WEIGHT CALCULATIONS -------------
double w = 1.0;
bool isCC1pi0AtVertex = false;
// Get W
GHepParticle* neutrino = ghep->Probe();
GHepParticle* fsl = ghep->FinalStatePrimaryLepton();
const TLorentzVector& k1 = *(neutrino->P4());
const TLorentzVector& k2 = *(fsl->P4());
double E_mu = k2.E();
double p_mu = k2.Vect().Mag();
double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu);
double th_nu_mu = k1.Vect().Angle(k2.Vect());
// The factor of 1000 is necessary for downstream functions
const double m_p = PhysConst::mass_proton;
double E_nu = k1.E();
// double hadMass = kine.W (true);
double hadMass = sqrt(m_p * m_p + m_mu * m_mu - 2 * m_p * E_mu + \
2 * E_nu * (m_p - E_mu + p_mu * cos(th_nu_mu)));
// Determine if event is CC1pi0 at vertex
TObjArrayIter piter(ghep);
GHepParticle * p = 0;
int pi0 = 0;
int piother = 0;
while ( (p = (GHepParticle *) piter.Next()) ) {
int pdgc = p->Pdg();
int ist = p->Status();
if (pdg::IsPseudoParticle(p->Pdg())) continue;
if (ist == genie::kIStStableFinalState) {
if (pdgc == 111) {
pi0++;
} else if (pdgc == 211 || pdgc == -211) {
piother++;
}
}
}
// std::cout << "Npi0 = " << pi0 << std::endl;
isCC1pi0AtVertex = (pi0 == 1 && piother == 0);
// Apply Weight
if (fApply_Enhancement && isCC1pi0AtVertex) {
double enhancement = 1.0 + fCur_WNorm * (exp( -0.5 * pow((fCur_WMean - hadMass) / (fCur_WSigma), 2) ));
w *= enhancement;
}
// Return Combined Weight
if (isnan(w) || w < 0.0 || w > 400.0) {
w = 1.0;
}
return w;
}
void WEnhancement::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void WEnhancement::SetDialValue(int rwenum, double val) {
// Check Handled
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
if (!IsHandled(curenum)) return;
// Set Values
if (curenum == Reweight::kMINERvARW_ApplyWTune) {
fApply_Enhancement = (val > 0.5);
}
if (curenum == Reweight::kMINERvARW_NormWTune) {
fTwk_WNorm = val;
fCur_WNorm = fDef_WNorm * (1.0 + 0.1 * fTwk_WNorm);
}
if (curenum == Reweight::kMINERvARW_MeanWTune) {
fTwk_WMean = val;
fCur_WMean = fDef_WMean * (1.0 + 0.1 * fTwk_WMean);
}
if (curenum == Reweight::kMINERvARW_SigmaWTune) {
fTwk_WSigma = val;
fCur_WSigma = fDef_WSigma * (1.0 + 0.1 * fTwk_WSigma);
}
// Define Tweaked
fTweaked = (fApply_Enhancement);
}
bool WEnhancement::IsHandled(int rwenum) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::kMINERvARW_ApplyWTune:
case Reweight::kMINERvARW_NormWTune:
case Reweight::kMINERvARW_MeanWTune:
case Reweight::kMINERvARW_SigmaWTune:
return true;
default:
return false;
}
}
} // namespace reweight
} // namespace nuisance
#endif
#endif
diff --git a/src/Reweight/NEUTWeightEngine.cxx b/src/Reweight/NEUTWeightEngine.cxx
index 66aed33..c4542b8 100644
--- a/src/Reweight/NEUTWeightEngine.cxx
+++ b/src/Reweight/NEUTWeightEngine.cxx
@@ -1,208 +1,208 @@
#include "NEUTWeightEngine.h"
NEUTWeightEngine::NEUTWeightEngine(std::string name) {
-#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
+#if defined(__NEUT_ENABLED__) and defined(__USE_NEUT_REWEIGHT__)
#if defined(__NEUT_VERSION__) && (__NEUT_VERSION__ >= 541)
std::string neut_card = FitPar::Config().GetParS("NEUT_CARD");
if (!neut_card.size()) {
NUIS_ABORT(
"[ERROR]: When using NEUTReWeight must set NEUT_CARD config option.");
}
// No need to vomit the contents of the card file all over my screen
StopTalking();
neut::CommonBlockIFace::Initialize(neut_card);
StartTalking();
#endif
// Setup the NEUT Reweight engien
fCalcName = name;
NUIS_LOG(FIT, "Setting up NEUT RW : " << fCalcName);
// Create RW Engine suppressing cout
StopTalking();
fNeutRW = new neut::rew::NReWeight();
TDirectory *olddir = gDirectory;
// get list of vetoed calc engines (just for debug really)
std::string rw_engine_list =
FitPar::Config().GetParS("FitWeight_fNeutRW_veto");
bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos;
bool xsec_res = rw_engine_list.find("xsec_res") == std::string::npos;
// Activate each calc engine
if (xsec_ccqe)
fNeutRW->AdoptWghtCalc("xsec_ccqe", new neut::rew::NReWeightNuXSecCCQE);
if (xsec_res)
fNeutRW->AdoptWghtCalc("xsec_res", new neut::rew::NReWeightNuXSecRES);
// Dials removed in NEUT 5.4.1
#if __NEUT_VERSION__ < 541
bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos;
bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos;
bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos;
bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos;
bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos;
bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos;
bool nucl_casc = rw_engine_list.find("nucl_casc") == std::string::npos;
bool nucl_piless = rw_engine_list.find("nucl_piless") == std::string::npos;
if (nucl_casc)
fNeutRW->AdoptWghtCalc("nucl_casc", new neut::rew::NReWeightCasc);
if (xsec_coh)
fNeutRW->AdoptWghtCalc("xsec_coh", new neut::rew::NReWeightNuXSecCOH);
if (xsec_nc)
fNeutRW->AdoptWghtCalc("xsec_nc", new neut::rew::NReWeightNuXSecNC);
if (nucl_piless)
fNeutRW->AdoptWghtCalc("nucl_piless", new neut::rew::NReWeightNuclPiless);
if (xsec_ncres)
fNeutRW->AdoptWghtCalc("xsec_ncres", new neut::rew::NReWeightNuXSecNCRES);
if (xsec_ccres)
fNeutRW->AdoptWghtCalc("xsec_ccres", new neut::rew::NReWeightNuXSecCCRES);
if (xsec_ncel)
fNeutRW->AdoptWghtCalc("xsec_ncel", new neut::rew::NReWeightNuXSecNCEL);
if (xsec_dis)
fNeutRW->AdoptWghtCalc("xsec_dis", new neut::rew::NReWeightNuXSecDIS);
#endif
fNeutRW->Reconfigure();
olddir->cd();
// Set Abs Twk Config
fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
// allow cout again
StartTalking();
#else
NUIS_ABORT("NEUT RW NOT ENABLED!");
#endif
};
void NEUTWeightEngine::IncludeDial(std::string name, double startval) {
-#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
+#if defined(__NEUT_ENABLED__) and defined(__USE_NEUT_REWEIGHT__)
// Get First enum
int nuisenum = Reweight::ConvDial(name, kNEUT);
// Setup Maps
fEnumIndex[nuisenum]; // = std::vector(0);
fNameIndex[name]; // = std::vector(0);
// Split by commas
std::vector allnames = GeneralUtils::ParseToStr(name, ",");
for (uint i = 0; i < allnames.size(); i++) {
std::string singlename = allnames[i];
// Get Syst
#if __NEUT_VERSION__ < 541
neut::rew::NSyst_t gensyst = NSyst::FromString(singlename);
#else
neut::rew::NSyst_t gensyst = neut::rew::NSyst::FromString(singlename);
#endif
// Fill Maps
int index = fValues.size();
fValues.push_back(0.0);
fNEUTSysts.push_back(gensyst);
// Initialize dial
NUIS_LOG(FIT, "Registering " << singlename << " dial.");
fNeutRW->Systematics().Init(fNEUTSysts[index]);
// If Absolute
if (fIsAbsTwk) {
#if __NEUT_VERSION__ < 541
NSystUncertainty::Instance()->SetUncertainty(fNEUTSysts[index], 1.0, 1.0);
#else
neut::rew::NSystUncertainty::Instance()->SetUncertainty(fNEUTSysts[index], 1.0, 1.0);
#endif
}
// Setup index
fEnumIndex[nuisenum].push_back(index);
fNameIndex[name].push_back(index);
}
// Set Value if given
if (startval != _UNDEF_DIAL_VALUE_) {
SetDialValue(nuisenum, startval);
}
#endif
}
void NEUTWeightEngine::SetDialValue(int nuisenum, double val) {
-#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
+#if defined(__NEUT_ENABLED__) and defined(__USE_NEUT_REWEIGHT__)
std::vector indices = fEnumIndex[nuisenum];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
NUIS_LOG(FIT, "Setting Dial Value for "
<< nuisenum << " " << i << " " << indices[i] << " "
<< fValues[indices[i]]
<< " Enum: " << fNEUTSysts[indices[i]]);
fNeutRW->Systematics().Set(fNEUTSysts[indices[i]], val);
}
#endif
}
void NEUTWeightEngine::SetDialValue(std::string name, double val) {
-#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
+#if defined(__NEUT_ENABLED__) and defined(__USE_NEUT_REWEIGHT__)
std::vector indices = fNameIndex[name];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
NUIS_LOG(FIT, "Setting Dial Value for "
<< name << " = " << i << " " << indices[i] << " "
<< fValues[indices[i]]
<< " Enum: " << fNEUTSysts[indices[i]]);
fNeutRW->Systematics().Set(fNEUTSysts[indices[i]], val);
}
#endif
}
void NEUTWeightEngine::Reconfigure(bool silent) {
-#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
+#if defined(__NEUT_ENABLED__) and defined(__USE_NEUT_REWEIGHT__)
// Hush now...
if (silent)
StopTalking();
// Reconf
fNeutRW->Reconfigure();
// if (LOG_LEVEL(DEB)){
fNeutRW->Print();
// }
// Shout again
if (silent)
StartTalking();
#endif
}
double NEUTWeightEngine::CalcWeight(BaseFitEvt *evt) {
double rw_weight = 1.0;
-#if defined(__NEUT_ENABLED__) and !defined(__NO_REWEIGHT__)
+#if defined(__NEUT_ENABLED__) and defined(__USE_NEUT_REWEIGHT__)
// Skip Non NEUT
if (evt->fType != kNEUT)
return 1.0;
// Hush now
StopTalking();
// Fill NEUT Common blocks
NEUTUtils::FillNeutCommons(evt->fNeutVect);
// Call Weight calculation
rw_weight = fNeutRW->CalcWeight();
// Speak Now
StartTalking();
#endif
if (!std::isnormal(rw_weight)){
NUIS_ERR(WRN, "NEUT returned weight: " << rw_weight);
rw_weight = 0;
}
// Return rw_weight
return rw_weight;
}
diff --git a/src/Reweight/NEUTWeightEngine.h b/src/Reweight/NEUTWeightEngine.h
index 17ff6ca..009cd9f 100644
--- a/src/Reweight/NEUTWeightEngine.h
+++ b/src/Reweight/NEUTWeightEngine.h
@@ -1,61 +1,57 @@
#ifndef WEIGHT_ENGINE_NEUT_H
#define WEIGHT_ENGINE_NEUT_H
#include "FitLogger.h"
-#ifdef __NEUT_ENABLED__
-#ifndef __NO_REWEIGHT__
+#if defined(__NEUT_ENABLED__) and defined(__USE_NEUT_REWEIGHT__)
#include "NEUTInputHandler.h"
#include "NReWeight.h"
#include "NReWeightNuXSecCCQE.h"
#include "NReWeightNuXSecRES.h"
// Dials removed in NEUT 5.4.1
#if __NEUT_VERSION__ < 541
#include "NReWeightCasc.h"
#include "NReWeightNuclPiless.h"
#include "NReWeightNuXSecNCRES.h"
#include "NReWeightNuXSecCCRES.h"
#include "NReWeightNuXSecNC.h"
#include "NReWeightNuXSecCOH.h"
#include "NReWeightNuXSecNCEL.h"
#include "NReWeightNuXSecDIS.h"
#endif
#if __NEUT_VERSION__ >= 541
#include "CommonBlockIFace.h"
#endif
#include "NSyst.h"
#include "NSystUncertainty.h"
#include "neutpart.h"
#include "neutvect.h"
#endif
-#endif
#include "FitWeight.h"
#include "GeneratorUtils.h"
#include "WeightEngineBase.h"
class NEUTWeightEngine : public WeightEngineBase {
public:
NEUTWeightEngine(std::string name);
~NEUTWeightEngine(){};
void IncludeDial(std::string name, double startval);
void SetDialValue(std::string name, double val);
void SetDialValue(int nuisenum, double val);
void Reconfigure(bool silent = false);
double CalcWeight(BaseFitEvt *evt);
inline bool NeedsEventReWeight() { return true; };
-#ifdef __NEUT_ENABLED__
-#ifndef __NO_REWEIGHT__
+#if defined(__NEUT_ENABLED__) and defined(__USE_NEUT_REWEIGHT__)
std::vector fNEUTSysts;
neut::rew::NReWeight *fNeutRW;
#endif
-#endif
};
#endif
diff --git a/src/Reweight/NOvARwgtEngine.cxx b/src/Reweight/NOvARwgtEngine.cxx
index 7a60320..e184834 100644
--- a/src/Reweight/NOvARwgtEngine.cxx
+++ b/src/Reweight/NOvARwgtEngine.cxx
@@ -1,317 +1,317 @@
#include "NOvARwgtEngine.h"
#include "NOvARwgt/interfaces/GenieInterface.h"
#include "NOvARwgt/rwgt/tunes/Tunes2017.h"
#include "NOvARwgt/rwgt/tunes/Tunes2018.h"
#include "NOvARwgt/rwgt/tunes/Tunes2020.h"
#include "NOvARwgt/rwgt/tunes/TunesSA.h"
#include "NOvARwgt/rwgt/EventRecord.h"
#include "NOvARwgt/rwgt/ISystKnob.h"
#include "NOvARwgt/rwgt/Tune.h"
// GENIEv3
#include "Framework/Utils/AppInit.h"
#include "Framework/Utils/RunOpt.h"
static size_t const kCVTune2017 = 100;
static size_t const kCVTune2018 = 200;
static size_t const kCVTune2020 = 300;
static size_t const kCVTuneSA = 400;
static size_t const kNoSuchKnob = std::numeric_limits::max();
novarwgt::Tune const *TuneFactory(size_t e) {
switch (e) {
case kCVTune2017: {
return &novarwgt::kCVTune2017;
}
case kCVTune2018: {
return &novarwgt::kCVTune2018;
}
case kCVTune2020: {
return &novarwgt::kCVTune2020;
}
case kCVTuneSA: {
return &novarwgt::kCVTuneSA;
}
default: {
return NULL;
}
}
}
void NOvARwgtEngine::InitializeKnobs() {
fTunes.push_back(TuneFactory(kCVTune2017));
fTuneEnums[kCVTune2017] = 0;
fTunes.push_back(TuneFactory(kCVTune2018));
fTuneEnums[kCVTune2018] = 1;
fTunes.push_back(TuneFactory(kCVTune2020));
fTuneEnums[kCVTune2020] = 2;
fTunes.push_back(TuneFactory(kCVTuneSA));
fTuneEnums[kCVTuneSA] = 3;
fTuneInUse = {false, false, false, false};
fTuneValues = {0, 0, 0, 0};
size_t ctr = 0;
NUIS_LOG(FIT, "NOvARwgt kCVTune2017 sub-knobs:");
size_t tune_ctr = 0;
for (auto &k : novarwgt::kCVTune2017.SystKnobs()) {
NUIS_LOG(FIT, "\t" << k.first);
fKnobs.push_back(k.second);
fKnobTunes.push_back(&novarwgt::kCVTune2017);
fKnobTuneidx.push_back(0);
fKnobEnums[kCVTune2017 + 1 + tune_ctr++] = ctr++;
fKnobInUse.push_back(false);
fKnobValues.push_back(0);
}
tune_ctr = 0;
NUIS_LOG(FIT, "NOvARwgt kCVTune2018 sub-knobs:");
for (auto &k : novarwgt::kCVTune2018.SystKnobs()) {
NUIS_LOG(FIT, "\t" << k.first);
fKnobs.push_back(k.second);
fKnobTunes.push_back(&novarwgt::kCVTune2018);
fKnobTuneidx.push_back(1);
fKnobEnums[kCVTune2018 + 1 + tune_ctr++] = ctr++;
fKnobInUse.push_back(false);
fKnobValues.push_back(0);
}
tune_ctr = 0;
NUIS_LOG(FIT, "NOvARwgt kCVTune2020 sub-knobs:");
for (auto &k : novarwgt::kCVTune2020.SystKnobs()) {
NUIS_LOG(FIT, "\t" << k.first);
fKnobs.push_back(k.second);
fKnobTunes.push_back(&novarwgt::kCVTune2020);
fKnobTuneidx.push_back(2);
fKnobEnums[kCVTune2020 + 1 + tune_ctr++] = ctr++;
fKnobInUse.push_back(false);
fKnobValues.push_back(0);
}
tune_ctr = 0;
NUIS_LOG(FIT, "NOvARwgt kCVTuneSA sub-knobs:");
for (auto &k : novarwgt::kCVTuneSA.SystKnobs()) {
NUIS_LOG(FIT, "\t" << k.first);
fKnobs.push_back(k.second);
fKnobTunes.push_back(&novarwgt::kCVTuneSA);
fKnobTuneidx.push_back(3);
fKnobEnums[kCVTuneSA + 1 + tune_ctr++] = ctr++;
fKnobInUse.push_back(false);
fKnobValues.push_back(0);
}
}
void NOvARwgtEngine::InitializeGENIE() {
genie::RunOpt *grunopt = genie::RunOpt::Instance();
grunopt->EnableBareXSecPreCalc(true);
grunopt->SetEventGeneratorList(Config::GetParS("GENIEEventGeneratorList"));
if (!Config::HasPar("GENIETune")) {
NUIS_ABORT(
"GENIE tune was not specified, this is required when reweighting GENIE "
"V3+ events. Add a config parameter like: to your nuisance card.");
}
grunopt->SetTuneName(Config::GetParS("GENIETune"));
grunopt->BuildTune();
std::string genv =
std::string(getenv("GENIE")) + "/config/Messenger_laconic.xml";
genie::utils::app_init::MesgThresholds(genv);
}
size_t NOvARwgtEngine::GetWeightGeneratorIndex(std::string const &strname) {
size_t upos = strname.find_first_of("_");
if (strname.find("CVTune2017") == 0) {
if (upos == std::string::npos) {
return kCVTune2017;
}
std::string knobname = strname.substr(upos + 1);
if (novarwgt::kCVTune2017.SystKnobs().count(knobname)) {
auto loc = novarwgt::kCVTune2017.SystKnobs().find(knobname);
return kCVTune2017 + 1 +
std::distance(novarwgt::kCVTune2017.SystKnobs().begin(), loc);
}
} else if (strname.find("CVTune2018") == 0) {
if (upos == std::string::npos) {
return kCVTune2018;
}
std::string knobname = strname.substr(upos + 1);
if (novarwgt::kCVTune2018.SystKnobs().count(knobname)) {
auto loc = novarwgt::kCVTune2018.SystKnobs().find(knobname);
return kCVTune2018 + 1 +
std::distance(novarwgt::kCVTune2018.SystKnobs().begin(), loc);
}
} else if (strname.find("CVTune2020") == 0) {
if (upos == std::string::npos) {
return kCVTune2020;
}
std::string knobname = strname.substr(upos + 1);
if (novarwgt::kCVTune2020.SystKnobs().count(knobname)) {
auto loc = novarwgt::kCVTune2020.SystKnobs().find(knobname);
return kCVTune2020 + 1 +
std::distance(novarwgt::kCVTune2020.SystKnobs().begin(), loc);
}
} else if (strname.find("CVTuneSA") == 0) {
if (upos == std::string::npos) {
return kCVTuneSA;
}
std::string knobname = strname.substr(upos + 1);
if (novarwgt::kCVTuneSA.SystKnobs().count(knobname)) {
auto loc = novarwgt::kCVTuneSA.SystKnobs().find(knobname);
return kCVTuneSA + 1 +
std::distance(novarwgt::kCVTuneSA.SystKnobs().begin(), loc);
}
}
return kNoSuchKnob;
}
void NOvARwgtEngine::IncludeDial(std::string name, double startval) {
size_t we_indx = GetWeightGeneratorIndex(name);
if (we_indx == kNoSuchKnob) {
NUIS_ABORT("[ERROR]: Invalid NOvARwgt Engine name: " << name);
}
- bool IsTune = !(we_indx % 100);
+ bool IsTune = !(we_indx % (NUIS_DIAL_OFFSET /10));
if (IsTune) {
auto tune_idx = fTuneEnums[we_indx];
fTuneValues[tune_idx] = startval;
fTuneInUse[tune_idx] = true;
} else {
auto knob_idx = fKnobEnums[we_indx];
fKnobValues[knob_idx] = startval;
fKnobInUse[knob_idx] = true;
}
};
void NOvARwgtEngine::SetDialValue(int nuisenum, double val) {
- size_t we_indx = (nuisenum % 1000);
- bool IsTune = !(we_indx % 100);
+ size_t we_indx = (nuisenum % NUIS_DIAL_OFFSET);
+ bool IsTune = !(we_indx % (NUIS_DIAL_OFFSET /10));
if (IsTune) {
auto tune_idx = fTuneEnums[we_indx];
if (!fTuneInUse[tune_idx]) {
NUIS_ABORT("[ERROR]: SetDialValue for NOvARwgt dial: "
<< we_indx << " but that tune hasn't been included yet.");
}
fTuneValues[tune_idx] = val;
} else {
auto knob_idx = fKnobEnums[we_indx];
if (!fKnobInUse[knob_idx]) {
NUIS_ABORT("[ERROR]: SetDialValue for NOvARwgt dial: "
<< we_indx << " but that Knob hasn't been included yet.");
}
fKnobValues[knob_idx] = val;
}
}
void NOvARwgtEngine::SetDialValue(std::string name, double val) {
SetDialValue(GetWeightGeneratorIndex(name), val);
}
bool NOvARwgtEngine::IsDialIncluded(std::string name) {
return IsDialIncluded(GetWeightGeneratorIndex(name));
}
bool NOvARwgtEngine::IsDialIncluded(int nuisenum) {
- size_t we_indx = (nuisenum % 1000);
- bool IsTune = !(we_indx % 100);
+ size_t we_indx = (nuisenum % NUIS_DIAL_OFFSET);
+ bool IsTune = !(we_indx % (NUIS_DIAL_OFFSET /10));
if (IsTune) {
auto tune_idx = fTuneEnums[we_indx];
return fTuneInUse[tune_idx];
} else {
auto knob_idx = fKnobEnums[we_indx];
return fKnobInUse[knob_idx];
}
}
double NOvARwgtEngine::GetDialValue(std::string name) {
return GetDialValue(GetWeightGeneratorIndex(name));
}
double NOvARwgtEngine::GetDialValue(int nuisenum) {
- size_t we_indx = (nuisenum % 1000);
+ size_t we_indx = (nuisenum % NUIS_DIAL_OFFSET);
if (we_indx == kNoSuchKnob) {
NUIS_ABORT("[ERROR]: Invalid NOvARwgt Engine enum: " << nuisenum);
}
- bool IsTune = !(we_indx % 100);
+ bool IsTune = !(we_indx % (NUIS_DIAL_OFFSET /10));
if (IsTune) {
auto tune_idx = fTuneEnums[we_indx];
if (!fTuneInUse[tune_idx]) {
NUIS_ABORT("[ERROR]: GetDialValue for NOvARwgt dial: "
<< we_indx << " but that tune hasn't been included yet.");
}
return fTuneValues[tune_idx];
} else {
auto knob_idx = fKnobEnums[we_indx];
if (!fKnobInUse[knob_idx]) {
NUIS_ABORT("[ERROR]: GetDialValue for NOvARwgt dial: "
<< we_indx << " but that Knob hasn't been included yet.");
}
return fKnobValues[knob_idx];
}
}
double NOvARwgtEngine::CalcWeight(BaseFitEvt *evt) {
double rw_weight = 1.0;
// Make nom weight
if (!evt) {
NUIS_ABORT("evt not found : " << evt);
}
// Skip Non GENIE
if (evt->fType != kGENIE)
return 1.0;
if (!(evt->genie_event)) {
NUIS_ABORT("evt->genie_event not found!" << evt->genie_event);
}
if (!(evt->genie_event->event)) {
NUIS_ABORT("evt->genie_event->event GHepRecord not found!"
<< (evt->genie_event->event));
}
novarwgt::EventRecord rcd =
novarwgt::ConvertGenieEvent(evt->genie_event->event);
for (size_t k_it = 0; k_it < fKnobs.size(); ++k_it) {
if (!fKnobInUse[k_it]) {
continue;
}
double wght = fKnobTunes[k_it]->EventSystKnobWeight(
fKnobs[k_it]->Name(), fKnobValues[k_it], rcd, {},
fTuneInUse[fKnobTuneidx[k_it]] && fTuneValues[fKnobTuneidx[k_it]]);
// // have to divide out the CV weight for this, ugly hack because the last
// // parameter doesn't do what I want
// if (fTuneInUse[fKnobTuneidx[k_it]] && fTuneValues[fKnobTuneidx[k_it]]) {
// wght /= fKnobTunes[k_it]->EventSystKnobWeight(fKnobs[k_it]->Name(), 0,
// rcd, {}, false);
// }
rw_weight *= wght;
}
for (size_t k_it = 0; k_it < fTunes.size(); ++k_it) {
if (!fTuneInUse[k_it]) {
continue;
}
if (!fTuneValues[k_it]) {
continue;
}
double wght = fTunes[k_it]->EventWeight(rcd);
rw_weight *= wght;
}
return rw_weight;
}
diff --git a/src/Reweight/NUISANCEWeightCalcs.cxx b/src/Reweight/NUISANCEWeightCalcs.cxx
index 426f582..7f9d8d2 100644
--- a/src/Reweight/NUISANCEWeightCalcs.cxx
+++ b/src/Reweight/NUISANCEWeightCalcs.cxx
@@ -1,792 +1,794 @@
#include "NUISANCEWeightCalcs.h"
#include "FitEvent.h"
#include "GeneralUtils.h"
#include "NUISANCESyst.h"
#include "WeightUtils.h"
+#include "WeightEngineBase.h"
+
using namespace Reweight;
ModeNormCalc::ModeNormCalc() { fNormRES = 1.0; }
double ModeNormCalc::CalcWeight(BaseFitEvt *evt) {
int mode = abs(evt->Mode);
double w = 1.0;
if (mode == 11 or mode == 12 or mode == 13) {
w *= fNormRES;
}
return w;
}
void ModeNormCalc::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void ModeNormCalc::SetDialValue(int rwenum, double val) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
// Check Handled
if (!IsHandled(curenum))
return;
if (curenum == kModeNorm_NormRES)
fNormRES = val;
}
bool ModeNormCalc::IsHandled(int rwenum) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case kModeNorm_NormRES:
return true;
default:
return false;
}
}
//*****************************************************************************
MINOSRPA::MINOSRPA() {
fApply_MINOSRPA = false;
fDef_MINOSRPA_A = 1.010;
fTwk_MINOSRPA_A = 0.000;
fDef_MINOSRPA_B = 0.156;
fTwk_MINOSRPA_B = 0.000;
}
//
double MINOSRPA::CalcWeight(BaseFitEvt* evt) {
if (!fTweaked || !fApply_MINOSRPA) return 1.0;
double w = 1.0;
// If GENIE is enabled, use old code
#ifdef __GENIE_ENABLED__
// Extract the GENIE Record
GHepRecord* ghep = static_cast(evt->genie_event->event);
const Interaction* interaction = ghep->Summary();
const InitialState& init_state = interaction->InitState();
const ProcessInfo& proc_info = interaction->ProcInfo();
const Target& tgt = init_state.Tgt();
// If not on nucleus, not resonant, or NC
if (!tgt.IsNucleus() || !proc_info.IsResonant() || proc_info.IsWeakNC()) return 1.0;
// Extract Beam and Target PDG
GHepParticle* neutrino = ghep->Probe();
// int bpdg = neutrino->Pdg();
GHepParticle* target = ghep->Particle(1);
assert(target);
// int tpdg = target->Pdg();
// Extract Q0-Q3
GHepParticle* fsl = ghep->FinalStatePrimaryLepton();
const TLorentzVector& k1 = *(neutrino->P4());
const TLorentzVector& k2 = *(fsl->P4());
//double q0 = fabs((k1 - k2).E());
//double q3 = fabs((k1 - k2).Vect().Mag());
double Q2 = fabs((k1 - k2).Mag2());
w *= GetRPAWeight(Q2);
#else
// Get the Q2 from NUISANCE if not GENIE
FitEvent *fevt = static_cast(evt);
// Check the event is resonant
if (!fevt->IsResonant() || fevt->IsNC()) return 1.0;
int targeta = fevt->GetTargetA();
int targetz = fevt->GetTargetZ();
// Apply only to nuclear targets, ignore free protons
if (targeta == 1 || targetz == 1) return 1.0;
// Q2 in GeV2
double Q2 = fevt->GetQ2();
w *= GetRPAWeight(Q2);
#endif
return w;
}
// Do the actual weight calculation
double MINOSRPA::GetRPAWeight(double Q2) {
if (Q2 > 0.7) return 1.0;
double w = fCur_MINOSRPA_A / (1.0 + TMath::Exp(1.0 - TMath::Sqrt(Q2) / fCur_MINOSRPA_B));
return w;
}
bool MINOSRPA::IsHandled(int rwenum) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::kMINERvARW_MINOSRPA_Apply:
case Reweight::kMINERvARW_MINOSRPA_A:
case Reweight::kMINERvARW_MINOSRPA_B:
return true;
default:
return false;
}
}
//
void MINOSRPA::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//
void MINOSRPA::SetDialValue(int rwenum, double val) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
// Check Handled
if (!IsHandled(curenum)) return;
if (curenum == Reweight::kMINERvARW_MINOSRPA_Apply) fApply_MINOSRPA = (val > 0.5);
if (curenum == Reweight::kMINERvARW_MINOSRPA_A) fTwk_MINOSRPA_A = val;
if (curenum == Reweight::kMINERvARW_MINOSRPA_B) fTwk_MINOSRPA_B = val;
// Check for changes
fTweaked = (fApply_MINOSRPA ||
fabs(fTwk_MINOSRPA_A) > 0.0 ||
fabs(fTwk_MINOSRPA_B) > 0.0);
// Update Values
fCur_MINOSRPA_A = fDef_MINOSRPA_A * (1.0 + 0.1 * fTwk_MINOSRPA_A);
fCur_MINOSRPA_B = fDef_MINOSRPA_B * (1.0 + 0.1 * fTwk_MINOSRPA_B);
}
//
//
//*****************************************************************************
LagrangeRPA::LagrangeRPA() {
fApplyRPA = false;
/*
fI1_Def = 4.18962e-01;
fI1 = fI1_Def;
fI2_Def = 7.39927e-01;
fI2 = fI2_Def;
*/
// Table VIII https://arxiv.org/pdf/1903.01558.pdf
fR1_Def = 0.37;
fR1 = fR1_Def;
fR2_Def = 0.60;
fR2 = fR2_Def;
}
//
double LagrangeRPA::CalcWeight(BaseFitEvt* evt) {
if (!fTweaked || !fApplyRPA) return 1.0;
double w = 1.0;
// If GENIE is enabled, use old code
#ifdef __GENIE_ENABLED__
// Extract the GENIE Record
GHepRecord* ghep = static_cast(evt->genie_event->event);
const Interaction* interaction = ghep->Summary();
const InitialState& init_state = interaction->InitState();
const ProcessInfo& proc_info = interaction->ProcInfo();
const Target& tgt = init_state.Tgt();
// If not on nucleus, not resonant, or NC
if (!tgt.IsNucleus() || !proc_info.IsResonant() || proc_info.IsWeakNC()) return 1.0;
// Extract Beam and Target PDG
GHepParticle* neutrino = ghep->Probe();
// int bpdg = neutrino->Pdg();
GHepParticle* target = ghep->Particle(1);
assert(target);
// int tpdg = target->Pdg();
// Extract Q0-Q3
GHepParticle* fsl = ghep->FinalStatePrimaryLepton();
const TLorentzVector& k1 = *(neutrino->P4());
const TLorentzVector& k2 = *(fsl->P4());
//double q0 = fabs((k1 - k2).E());
//double q3 = fabs((k1 - k2).Vect().Mag());
double Q2 = fabs((k1 - k2).Mag2());
w *= GetRPAWeight(Q2);
#else
// Get the Q2 from NUISANCE if not GENIE
FitEvent *fevt = static_cast(evt);
// Check the event is resonant
if (!fevt->IsResonant() || fevt->IsNC()) return 1.0;
int targeta = fevt->GetTargetA();
int targetz = fevt->GetTargetZ();
// Apply only to nuclear targets, ignore free protons
if (targeta == 1 || targetz == 1) return 1.0;
// Q2 in GeV2
double Q2 = fevt->GetQ2();
w *= GetRPAWeight(Q2);
#endif
return w;
}
//
//
double LagrangeRPA::GetRPAWeight(double Q2) {
//std::cout << "Getting RPA Weight : " << Q2 << std::endl;
if (Q2 > 0.7) return 1.0;
// Keep original Lagrange RPA for documentation
/*
double x1 = 0.00;
double x2 = 0.30;
double x3 = 0.70;
double y1 = 0.00;
double y2 = fI2;
double y3 = 1.00;
double xv = Q2;
// Automatically 0 because y1 choice
double W1 = y1 * (xv-x2)*(xv-x3)/((x1-x2)*(x1-x3));
double W2 = y2 * (xv-x1)*(xv-x3)/((x2-x1)*(x2-x3));
double W3 = y3 * (xv-x1)*(xv-x2)/((x3-x1)*(x3-x2));
double P = W1 + W2 + W3;
double A1 = (1.0 - sqrt(1.0 - fI1));
double R = P * (1.0 - A1) + A1;
return 1.0 - (1.0-R)*(1.0-R);
*/
// Equation 7 https://arxiv.org/pdf/1903.01558.pdf
const double x1 = 0.00;
const double x2 = 0.35;
const double x3 = 0.70;
// Equation 6 https://arxiv.org/pdf/1903.01558.pdf
double RQ2 = fR2 *( (Q2-x1)*(Q2-x3)/((x2-x1)*(x2-x3)) )
+ (Q2-x1)*(Q2-x2)/((x3-x1)*(x3-x2));
double weight = 1-(1-fR1)*(1-RQ2)*(1-RQ2);
// Check range of R1 and R2
// Commented out because this is implementation dependent: user may want strange double peaks
/*
if (fR1 > 1) return 1;
if (fR2 > 1 || fR2 < 0.5) return 1;
*/
return weight;
}
//
bool LagrangeRPA::IsHandled(int rwenum) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case Reweight::kMINERvARW_LagrangeRPA_Apply:
case Reweight::kMINERvARW_LagrangeRPA_R1:
case Reweight::kMINERvARW_LagrangeRPA_R2:
return true;
default:
return false;
}
}
//
void LagrangeRPA::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
//
void LagrangeRPA::SetDialValue(int rwenum, double val) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
// Check Handled
if (!IsHandled(curenum)) return;
if (curenum == Reweight::kMINERvARW_LagrangeRPA_Apply) fApplyRPA = (val > 0.5);
if (curenum == Reweight::kMINERvARW_LagrangeRPA_R1) fR1 = val;
if (curenum == Reweight::kMINERvARW_LagrangeRPA_R2) fR2 = val;
// Check for changes
fTweaked = (fApplyRPA);
}
//
BeRPACalc::BeRPACalc()
: fBeRPA_A(0.59), fBeRPA_B(1.05), fBeRPA_D(1.13), fBeRPA_E(0.88),
fBeRPA_U(1.2), nParams(0) {
// A = 0.59 +/- 20%
// B = 1.05 +/- 20%
// D = 1.13 +/- 15%
// E = 0.88 +/- 40%
// U = 1.2
}
double BeRPACalc::CalcWeight(BaseFitEvt *evt) {
int mode = abs(evt->Mode);
double w = 1.0;
if (nParams == 0) {
return w;
}
// Get Q2
// Get final state lepton
if (mode == 1) {
FitEvent *fevt = static_cast(evt);
double Q2 = fevt->GetQ2();
// Only CCQE events
w *= calcRPA(Q2, fBeRPA_A, fBeRPA_B, fBeRPA_D, fBeRPA_E, fBeRPA_U);
}
return w;
}
void BeRPACalc::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void BeRPACalc::SetDialValue(int rwenum, double val) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
// Check Handled
if (!IsHandled(curenum))
return;
// Need 4 or 5 reconfigures
if (curenum == kBeRPA_A)
fBeRPA_A = val;
else if (curenum == kBeRPA_B)
fBeRPA_B = val;
else if (curenum == kBeRPA_D)
fBeRPA_D = val;
else if (curenum == kBeRPA_E)
fBeRPA_E = val;
else if (curenum == kBeRPA_U)
fBeRPA_U = val;
nParams++;
}
bool BeRPACalc::IsHandled(int rwenum) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case kBeRPA_A:
case kBeRPA_B:
case kBeRPA_D:
case kBeRPA_E:
case kBeRPA_U:
return true;
default:
return false;
}
}
SBLOscWeightCalc::SBLOscWeightCalc() {
fDistance = 0.0;
fMassSplitting = 0.0;
fSin2Theta = 0.0;
}
double SBLOscWeightCalc::CalcWeight(BaseFitEvt *evt) {
FitEvent *fevt = static_cast(evt);
FitParticle *pnu = fevt->PartInfo(0);
double E = pnu->fP.E() / 1.E3;
// Extract energy
return GetSBLOscWeight(E);
}
void SBLOscWeightCalc::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void SBLOscWeightCalc::SetDialValue(int rwenum, double val) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
if (!IsHandled(curenum))
return;
if (curenum == kSBLOsc_Distance)
fDistance = val;
if (curenum == kSBLOsc_MassSplitting)
fMassSplitting = val;
if (curenum == kSBLOsc_Sin2Theta)
fSin2Theta = val;
}
bool SBLOscWeightCalc::IsHandled(int rwenum) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case kSBLOsc_Distance:
return true;
case kSBLOsc_MassSplitting:
return true;
case kSBLOsc_Sin2Theta:
return true;
default:
return false;
}
}
double SBLOscWeightCalc::GetSBLOscWeight(double E) {
if (E <= 0.0)
return 1.0 - 0.5 * fSin2Theta;
return 1.0 - fSin2Theta * pow(sin(1.267 * fMassSplitting * fDistance / E), 2);
}
GaussianModeCorr::GaussianModeCorr() {
// Apply the tilt-shift Gauss by Patrick
// Alternatively set in config
fMethod = true;
// Init
fApply_CCQE = false;
fGausVal_CCQE[kPosNorm] = 0.0;
fGausVal_CCQE[kPosTilt] = 0.0;
fGausVal_CCQE[kPosPq0] = 1.0;
fGausVal_CCQE[kPosWq0] = 1.0;
fGausVal_CCQE[kPosPq3] = 1.0;
fGausVal_CCQE[kPosWq3] = 1.0;
fApply_2p2h = false;
fGausVal_2p2h[kPosNorm] = 0.0;
fGausVal_2p2h[kPosTilt] = 0.0;
fGausVal_2p2h[kPosPq0] = 1.0;
fGausVal_2p2h[kPosWq0] = 1.0;
fGausVal_2p2h[kPosPq3] = 1.0;
fGausVal_2p2h[kPosWq3] = 1.0;
fApply_2p2h_PPandNN = false;
fGausVal_2p2h_PPandNN[kPosNorm] = 0.0;
fGausVal_2p2h_PPandNN[kPosTilt] = 0.0;
fGausVal_2p2h_PPandNN[kPosPq0] = 1.0;
fGausVal_2p2h_PPandNN[kPosWq0] = 1.0;
fGausVal_2p2h_PPandNN[kPosPq3] = 1.0;
fGausVal_2p2h_PPandNN[kPosWq3] = 1.0;
fApply_2p2h_NP = false;
fGausVal_2p2h_NP[kPosNorm] = 0.0;
fGausVal_2p2h_NP[kPosTilt] = 0.0;
fGausVal_2p2h_NP[kPosPq0] = 1.0;
fGausVal_2p2h_NP[kPosWq0] = 1.0;
fGausVal_2p2h_NP[kPosPq3] = 1.0;
fGausVal_2p2h_NP[kPosWq3] = 1.0;
fApply_CC1pi = false;
fGausVal_CC1pi[kPosNorm] = 0.0;
fGausVal_CC1pi[kPosTilt] = 0.0;
fGausVal_CC1pi[kPosPq0] = 1.0;
fGausVal_CC1pi[kPosWq0] = 1.0;
fGausVal_CC1pi[kPosPq3] = 1.0;
fGausVal_CC1pi[kPosWq3] = 1.0;
fAllowSuppression = false;
fDebugStatements = FitPar::Config().GetParB("GaussianModeCorr_DEBUG");
// fDebugStatements = true;
}
double GaussianModeCorr::CalcWeight(BaseFitEvt *evt) {
FitEvent *fevt = static_cast(evt);
double rw_weight = 1.0;
// Get Neutrino
if (!fevt->Npart()) {
NUIS_ABORT("NO particles found in stack!");
}
FitParticle *pnu = fevt->GetNeutrinoIn();
if (!pnu) {
NUIS_ABORT("NO Starting particle found in stack!");
}
FitParticle *plep = fevt->GetLeptonOut();
if (!plep) return 1.0;
TLorentzVector q = pnu->fP - plep->fP;
// Extra q0,q3
double q0 = fabs(q.E()) / 1.E3;
double q3 = fabs(q.Vect().Mag()) / 1.E3;
int initialstate = -1; // Undef
if (abs(fevt->Mode) == 2) {
int npr = 0;
int nne = 0;
for (UInt_t j = 0; j < fevt->Npart(); j++) {
if ((fevt->PartInfo(j))->fIsAlive)
continue;
if (fevt->PartInfo(j)->fPID == 2212)
npr++;
else if (fevt->PartInfo(j)->fPID == 2112)
nne++;
}
if (fevt->Mode == 2 && npr == 1 && nne == 1) {
initialstate = 2;
} else if (fevt->Mode == 2 &&
((npr == 0 && nne == 2) || (npr == 2 && nne == 0))) {
initialstate = 1;
}
}
// Apply weighting
if (fApply_CCQE && abs(fevt->Mode) == 1) {
if (fDebugStatements)
std::cout << "Getting CCQE Weight" << std::endl;
double g = GetGausWeight(q0, q3, fGausVal_CCQE);
if (g < 1.0)
g = 1.0;
rw_weight *= g;
}
if (fApply_2p2h && abs(fevt->Mode) == 2) {
if (fDebugStatements) std::cout << "Getting 2p2h Weight" << std::endl;
if (fDebugStatements) std::cout << "Got q0 q3 = " << q0 << " " << q3 << " mode = " << fevt->Mode << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h);
if (fDebugStatements) std::cout << "Returning Weight " << rw_weight << std::endl;
}
if (fApply_2p2h_PPandNN && abs(fevt->Mode) == 2 && initialstate == 1) {
if (fDebugStatements) std::cout << "Getting 2p2h PPandNN Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_PPandNN);
}
if (fApply_2p2h_NP && abs(fevt->Mode) == 2 && initialstate == 2) {
if (fDebugStatements) std::cout << "Getting 2p2h NP Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_NP);
}
if (fApply_CC1pi && abs(fevt->Mode) >= 11 && abs(fevt->Mode) <= 13) {
if (fDebugStatements)
std::cout << "Getting CC1pi Weight" << std::endl;
rw_weight *= GetGausWeight(q0, q3, fGausVal_CC1pi);
}
return rw_weight;
}
void GaussianModeCorr::SetMethod(bool method) {
fMethod = method;
if (fMethod == true) {
NUIS_LOG(FIT,
" Using tilt-shift Gaussian parameters for Gaussian enhancement...");
} else {
NUIS_LOG(FIT, " Using Normal Gaussian parameters for Gaussian enhancement...");
}
};
double GaussianModeCorr::GetGausWeight(double q0, double q3, double vals[]) {
// The weight
double w = 1.0;
// Use tilt-shift method by Patrick
if (fMethod) {
if (fDebugStatements) {
std::cout << "Using Patrick gaussian" << std::endl;
}
// // CCQE Without Suppression
// double Norm = 4.82788679036;
// double Tilt = 2.3501416116;
// double Pq0 = 0.363964889702;
// double Wq0 = 0.133976806938;
// double Pq3 = 0.431769740224;
// double Wq3 = 0.207666663434;
// // Also add for CCQE at the end
// return (w > 1.0) ? w : 1.0;
// // 2p2h with suppression
// double Norm = 15.967;
// double Tilt = -0.455655;
// double Pq0 = 0.214598;
// double Wq0 = 0.0291061;
// double Pq3 = 0.480194;
// double Wq3 = 0.134588;
double Norm = vals[kPosNorm];
double Tilt = vals[kPosTilt];
double Pq0 = vals[kPosPq0];
double Wq0 = vals[kPosWq0];
double Pq3 = vals[kPosPq3];
double Wq3 = vals[kPosWq3];
double a = cos(Tilt) * cos(Tilt) / (2 * Wq0 * Wq0);
a += sin(Tilt) * sin(Tilt) / (2 * Wq3 * Wq3);
double b = -sin(2 * Tilt) / (4 * Wq0 * Wq0);
b += sin(2 * Tilt) / (4 * Wq3 * Wq3);
double c = sin(Tilt) * sin(Tilt) / (2 * Wq0 * Wq0);
c += cos(Tilt) * cos(Tilt) / (2 * Wq3 * Wq3);
w = Norm;
w *= exp(-a * (q0 - Pq0) * (q0 - Pq0));
w *= exp(+2.0 * b * (q0 - Pq0) * (q3 - Pq3));
w *= exp(-c * (q3 - Pq3) * (q3 - Pq3));
if (fDebugStatements) {
std::cout << "Applied Tilt " << Tilt << " " << cos(Tilt) << " "
<< sin(Tilt) << std::endl;
std::cout << "abc = " << a << " " << b << " " << c << std::endl;
std::cout << "Returning " << Norm << " " << Pq0 << " " << Wq0 << " "
<< Pq3 << " " << Wq3 << " " << w << std::endl;
}
if (w != w || std::isnan(w) || w < 0.0) {
w = 0.0;
}
if (w < 1.0 and !fAllowSuppression) {
w = 1.0;
}
return w;
// Use the MINERvA Gaussian method
} else {
/*
* From MINERvA and Daniel Ruterbories:
* Old notes here: *
* http://cdcvs.fnal.gov/cgi-bin/public-cvs/cvsweb-public.cgi/AnalysisFramework/Ana/CCQENu/ana_common/data/?cvsroot=mnvsoft
* These parameters are slightly altered
*
* FRESH:
* 10.5798
* 0.254032
* 0.50834
* 0.0571035
* 0.129051
* 0.875287
*/
if (fDebugStatements) {
std::cout << "Using MINERvA Gaussian" << std::endl;
}
double norm = vals[kPosNorm];
double meanq0 = vals[kPosTilt];
double meanq3 = vals[kPosPq0];
double sigmaq0 = vals[kPosWq0];
double sigmaq3 = vals[kPosPq3];
double corr = vals[kPosWq3];
double z = (q0 - meanq0) * (q0 - meanq0) / (sigmaq0 * sigmaq0) +
(q3 - meanq3) * (q3 - meanq3) / (sigmaq3 * sigmaq3) -
2 * corr * (q0 - meanq0) * (q3 - meanq3) / (sigmaq0 * sigmaq3);
double ret = 1.0;
if ( fabs(1 - corr*corr) < 1.E-5 ) {
return 1.0;
}
if ( (-0.5 * z / (1 - corr*corr)) > 200 or (-0.5 * z / (1 - corr*corr)) < -200 ) {
return 1.0;
} else {
ret = norm * exp( -0.5 * z / (1 - corr*corr) );
}
if (ret != ret or ret < 0.0 or isnan(ret)) {
return 1.0;
}
if (fAllowSuppression) return ret;
return ret + 1.0;
// Need to add 1 to the results
}
return w;
}
void GaussianModeCorr::SetDialValue(std::string name, double val) {
SetDialValue(Reweight::ConvDial(name, kCUSTOM), val);
}
void GaussianModeCorr::SetDialValue(int rwenum, double val) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
// Check Handled
if (!IsHandled(curenum)) return;
// CCQE Setting
for (int i = kGaussianCorr_CCQE_norm; i <= kGaussianCorr_CCQE_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_CCQE_norm;
fGausVal_CCQE[index] = val;
fApply_CCQE = true;
}
}
// 2p2h Setting
for (int i = kGaussianCorr_2p2h_norm; i <= kGaussianCorr_2p2h_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_2p2h_norm;
fGausVal_2p2h[index] = val;
fApply_2p2h = true;
}
}
// 2p2h_PPandNN Setting
for (int i = kGaussianCorr_2p2h_PPandNN_norm; i <= kGaussianCorr_2p2h_PPandNN_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_2p2h_PPandNN_norm;
fGausVal_2p2h_PPandNN[index] = val;
fApply_2p2h_PPandNN = true;
}
}
// 2p2h_NP Setting
for (int i = kGaussianCorr_2p2h_NP_norm; i <= kGaussianCorr_2p2h_NP_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_2p2h_NP_norm;
fGausVal_2p2h_NP[index] = val;
fApply_2p2h_NP = true;
}
}
// CC1pi Setting
for (int i = kGaussianCorr_CC1pi_norm; i <= kGaussianCorr_CC1pi_Wq3; i++) {
if (i == curenum) {
int index = i - kGaussianCorr_CC1pi_norm;
fGausVal_CC1pi[index] = val;
fApply_CC1pi = true;
}
}
if (curenum == kGaussianCorr_AllowSuppression) {
fAllowSuppression = (val > 0.5);
}
}
bool GaussianModeCorr::IsHandled(int rwenum) {
- int curenum = rwenum % 1000;
+ int curenum = rwenum % NUIS_DIAL_OFFSET;
switch (curenum) {
case kGaussianCorr_CCQE_norm:
case kGaussianCorr_CCQE_tilt:
case kGaussianCorr_CCQE_Pq0:
case kGaussianCorr_CCQE_Wq0:
case kGaussianCorr_CCQE_Pq3:
case kGaussianCorr_CCQE_Wq3:
case kGaussianCorr_2p2h_norm:
case kGaussianCorr_2p2h_tilt:
case kGaussianCorr_2p2h_Pq0:
case kGaussianCorr_2p2h_Wq0:
case kGaussianCorr_2p2h_Pq3:
case kGaussianCorr_2p2h_Wq3:
case kGaussianCorr_2p2h_PPandNN_norm:
case kGaussianCorr_2p2h_PPandNN_tilt:
case kGaussianCorr_2p2h_PPandNN_Pq0:
case kGaussianCorr_2p2h_PPandNN_Wq0:
case kGaussianCorr_2p2h_PPandNN_Pq3:
case kGaussianCorr_2p2h_PPandNN_Wq3:
case kGaussianCorr_2p2h_NP_norm:
case kGaussianCorr_2p2h_NP_tilt:
case kGaussianCorr_2p2h_NP_Pq0:
case kGaussianCorr_2p2h_NP_Wq0:
case kGaussianCorr_2p2h_NP_Pq3:
case kGaussianCorr_2p2h_NP_Wq3:
case kGaussianCorr_CC1pi_norm:
case kGaussianCorr_CC1pi_tilt:
case kGaussianCorr_CC1pi_Pq0:
case kGaussianCorr_CC1pi_Wq0:
case kGaussianCorr_CC1pi_Pq3:
case kGaussianCorr_CC1pi_Wq3:
case kGaussianCorr_AllowSuppression:
return true;
default:
return false;
}
}
diff --git a/src/Reweight/OscWeightEngine.cxx b/src/Reweight/OscWeightEngine.cxx
index 1224e02..43d9c91 100644
--- a/src/Reweight/OscWeightEngine.cxx
+++ b/src/Reweight/OscWeightEngine.cxx
@@ -1,331 +1,331 @@
// Copyright 2016-2021 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
//#define DEBUG_OSC_WE
#include "OscWeightEngine.h"
#include
enum nuTypes {
kNuebarType = -1,
kNumubarType = -2,
kNutaubarType = -3,
kNueType = 1,
kNumuType = 2,
kNutauType = 3,
};
nuTypes GetNuType(int pdg) {
switch (pdg) {
case 16:
return kNutauType;
case 14:
return kNumuType;
case 12:
return kNueType;
case -16:
return kNutaubarType;
case -14:
return kNumubarType;
case -12:
return kNuebarType;
default: { NUIS_ABORT("Attempting to convert \"neutrino pdg\": " << pdg); }
}
}
OscWeightEngine::OscWeightEngine()
:
#ifdef __PROB3PP_ENABLED__
bp(),
#endif
theta12(0.825),
theta13(0.10),
theta23(1.0),
dm12(7.9e-5),
dm23(2.5e-3),
dcp(0.0),
LengthParam(0xdeadbeef),
TargetNuType(0),
ForceFromNuPDG(0) {
Config();
}
void OscWeightEngine::Config() {
std::vector OscParam = Config::QueryKeys("OscParam");
if (OscParam.size() < 1) {
NUIS_ERR(WRN,
"Oscillation parameters specified but no OscParam element "
"configuring the experimental characteristics found.\nExpect at "
"least . Pausing for "
"10...");
sleep(10);
return;
}
if (OscParam[0].Has("baseline_km")) {
LengthParamIsZenith = false;
LengthParam = OscParam[0].GetD("baseline_km");
constant_density = OscParam[0].Has("matter_density")
? OscParam[0].GetD("matter_density")
: 0xdeadbeef;
} else if (OscParam[0].Has("detection_zenith_deg")) {
LengthParamIsZenith = true;
static const double deg2rad = asin(1) / 90.0;
LengthParam = cos(OscParam[0].GetD("detection_zenith_deg") * deg2rad);
} else {
NUIS_ERR(WRN,
"It appeared that you wanted to set up an oscillation weight "
"branch, but it was not correctly configured. You need to specify "
"either: detection_zenith_deg or baseline_km attributes on the "
"OscParam element, and if baseline_km is specified, you can "
"optionally also set matter_density for oscillations through a "
"constant matter density. Pausing for 10...");
sleep(10);
return;
}
dm23 = OscParam[0].Has("dm23") ? OscParam[0].GetD("dm23") : dm23;
theta23 = OscParam[0].Has("sinsq_theta23") ? OscParam[0].GetD("sinsq_theta23")
: theta23;
theta13 = OscParam[0].Has("sinsq_theta13") ? OscParam[0].GetD("sinsq_theta13")
: theta13;
dm12 = OscParam[0].Has("dm12") ? OscParam[0].GetD("dm12") : dm12;
theta12 = OscParam[0].Has("sinsq_theta12") ? OscParam[0].GetD("sinsq_theta12")
: theta12;
dcp = OscParam[0].Has("dcp") ? OscParam[0].GetD("dcp") : dcp;
TargetNuType = OscParam[0].Has("TargetNuPDG")
? GetNuType(OscParam[0].GetI("TargetNuPDG"))
: 0;
ForceFromNuPDG = OscParam[0].Has("ForceFromNuPDG")
? GetNuType(OscParam[0].GetI("ForceFromNuPDG"))
: 0;
NUIS_LOG(FIT, "Configured oscillation weighter:");
if (LengthParamIsZenith) {
NUIS_LOG(FIT,
"Earth density profile with detection cos(zenith) = " << LengthParam);
} else {
if (constant_density != 0xdeadbeef) {
NUIS_LOG(FIT,
"Constant density with experimental baseline = " << LengthParam);
} else {
NUIS_LOG(FIT,
"Vacuum oscillations with experimental baseline = " << LengthParam);
}
}
params[0] = dm23;
params[1] = theta23;
params[2] = theta13;
params[3] = dm12;
params[4] = theta12;
params[5] = dcp;
NUIS_LOG(FIT, "\tdm23 : " << params[0]);
NUIS_LOG(FIT, "\tsinsq_theta23: " << params[1]);
NUIS_LOG(FIT, "\tsinsq_theta13: " << params[2]);
NUIS_LOG(FIT, "\tdm12 : " << params[3]);
NUIS_LOG(FIT, "\tsinsq_theta12: " << params[4]);
NUIS_LOG(FIT, "\tdcp : " << params[5]);
if (TargetNuType) {
NUIS_LOG(FIT, "\tTargetNuType: " << TargetNuType);
}
if (ForceFromNuPDG) {
NUIS_LOG(FIT, "\tForceFromNuPDG: " << ForceFromNuPDG);
}
#ifdef __PROB3PP_ENABLED__
bp.SetMNS(params[theta12_idx], params[theta13_idx], params[theta23_idx],
params[dm12_idx], params[dm23_idx], params[dcp_idx], 1, true, 2);
bp.DefinePath(LengthParam, 0);
if (LengthParamIsZenith) {
NUIS_LOG(FIT, "\tBaseline : " << (bp.GetBaseline() / 100.0) << " km.");
}
#endif
}
void OscWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef DEBUG_OSC_WE
std::cout << "IncludeDial: " << name << " at " << startval << std::endl;
#endif
int dial = SystEnumFromString(name);
if (!dial) {
NUIS_ABORT("OscWeightEngine passed dial: " << name
<< " that it does not understand.");
}
params[dial - 1] = startval;
}
void OscWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef DEBUG_OSC_WE
- std::cout << "SetDial: " << (nuisenum % 1000) << " at " << val << std::endl;
+ std::cout << "SetDial: " << (nuisenum % NUIS_DIAL_OFFSET) << " at " << val << std::endl;
#endif
- fHasChanged = (params[(nuisenum % 1000) - 1] - val) >
+ fHasChanged = (params[(nuisenum % NUIS_DIAL_OFFSET) - 1] - val) >
std::numeric_limits::epsilon();
- params[(nuisenum % 1000) - 1] = val;
+ params[(nuisenum % NUIS_DIAL_OFFSET) - 1] = val;
}
void OscWeightEngine::SetDialValue(std::string name, double val) {
#ifdef DEBUG_OSC_WE
std::cout << "SetDial: " << name << " at " << val << std::endl;
#endif
int dial = SystEnumFromString(name);
if (!dial) {
NUIS_ABORT("OscWeightEngine passed dial: " << name
<< " that it does not understand.");
}
fHasChanged =
(params[dial - 1] - val) > std::numeric_limits::epsilon();
params[dial - 1] = val;
}
bool OscWeightEngine::IsDialIncluded(std::string name) {
return SystEnumFromString(name);
}
bool OscWeightEngine::IsDialIncluded(int nuisenum) {
- return ((nuisenum % 1000) > 0) && ((nuisenum % 1000) < 6);
+ return ((nuisenum % NUIS_DIAL_OFFSET) > 0) && ((nuisenum % NUIS_DIAL_OFFSET) < 6);
}
double OscWeightEngine::GetDialValue(std::string name) {
int dial = SystEnumFromString(name);
if (!dial) {
NUIS_ABORT("OscWeightEngine passed dial: " << name
<< " that it does not understand.");
}
return params[dial - 1];
}
double OscWeightEngine::GetDialValue(int nuisenum) {
- if (!(nuisenum % 1000) || (nuisenum % 1000) > 6) {
+ if (!(nuisenum % NUIS_DIAL_OFFSET) || (nuisenum % NUIS_DIAL_OFFSET) > 6) {
NUIS_ABORT("OscWeightEngine passed dial enum: "
- << (nuisenum % 1000)
+ << (nuisenum % NUIS_DIAL_OFFSET)
<< " that it does not understand, expected [1,6].");
}
- return params[(nuisenum % 1000) - 1];
+ return params[(nuisenum % NUIS_DIAL_OFFSET) - 1];
}
void OscWeightEngine::Reconfigure(bool silent) { fHasChanged = false; };
bool OscWeightEngine::NeedsEventReWeight() {
if (fHasChanged) {
return true;
}
return false;
}
double OscWeightEngine::CalcWeight(BaseFitEvt* evt) {
static bool Warned = false;
if (evt->probe_E == 0xdeadbeef) {
if (!Warned) {
NUIS_ERR(WRN,
"Oscillation weights asked for but using 'litemode' or "
"unsupported generator input. Pasuing for 10...");
sleep(10);
Warned = true;
}
return 1;
}
return CalcWeight(evt->probe_E * 1E-3, evt->probe_pdg);
}
double OscWeightEngine::CalcWeight(double ENu, int PDGNu, int TargetPDGNu) {
if (LengthParam == 0xdeadbeef) { // not configured.
return 1;
}
#ifdef __PROB3PP_ENABLED__
int NuType = (ForceFromNuPDG != 0) ? ForceFromNuPDG : GetNuType(PDGNu);
bp.SetMNS(params[theta12_idx], params[theta13_idx], params[theta23_idx],
params[dm12_idx], params[dm23_idx], params[dcp_idx], ENu, true,
NuType);
int pmt = 0;
double prob_weight = 1;
TargetPDGNu = (TargetPDGNu == -1) ? (TargetNuType ? TargetNuType : NuType)
: GetNuType(TargetPDGNu);
if (LengthParamIsZenith) { // Use earth density
bp.DefinePath(LengthParam, 0);
bp.propagate(NuType);
pmt = 0;
prob_weight = bp.GetProb(NuType, TargetPDGNu);
} else {
if (constant_density != 0xdeadbeef) {
bp.propagateLinear(NuType, LengthParam, constant_density);
pmt = 1;
prob_weight = bp.GetProb(NuType, TargetPDGNu);
} else {
pmt = 2;
prob_weight =
bp.GetVacuumProb(NuType, TargetPDGNu, ENu * 1E-3, LengthParam);
}
}
#ifdef DEBUG_OSC_WE
if (prob_weight != prob_weight) {
NUIS_ABORT("Calculated bad prob weight: " << prob_weight << "(Osc Type: " << pmt
<< " -- " << NuType << " -> "
<< TargetPDGNu << ")");
}
if (prob_weight > 1) {
NUIS_ABORT("Calculated bad prob weight: " << prob_weight << "(Osc Type: " << pmt
<< " -- " << NuType << " -> "
<< TargetPDGNu << ")");
}
std::cout << NuType << " -> " << TargetPDGNu << ": " << ENu << " = "
<< prob_weight << "%%." << std::endl;
#endif
return prob_weight;
#else
return 1;
#endif
}
int OscWeightEngine::SystEnumFromString(std::string const& name) {
if (name == "dm23") {
return 1;
} else if (name == "sinsq_theta23") {
return 2;
} else if (name == "sinsq_theta13") {
return 3;
} else if (name == "dm12") {
return 4;
} else if (name == "sinsq_theta12") {
return 5;
} else if (name == "dcp") {
return 6;
} else {
return 0;
}
}
void OscWeightEngine::Print() {
std::cout << "OscWeightEngine: " << std::endl;
std::cout << "\t theta12: " << params[theta12_idx] << std::endl;
std::cout << "\t theta13: " << params[theta13_idx] << std::endl;
std::cout << "\t theta23: " << params[theta23_idx] << std::endl;
std::cout << "\t dm12: " << params[dm12_idx] << std::endl;
std::cout << "\t dm23: " << params[dm23_idx] << std::endl;
std::cout << "\t dcp: " << params[dcp_idx] << std::endl;
}
diff --git a/src/Reweight/WeightEngineBase.h b/src/Reweight/WeightEngineBase.h
index 5a6e6cf..df326bd 100644
--- a/src/Reweight/WeightEngineBase.h
+++ b/src/Reweight/WeightEngineBase.h
@@ -1,56 +1,57 @@
#ifndef WEIGHTENGINE_BASE_H
#define WEIGHTENGINE_BASE_H
#include "BaseFitEvt.h"
#include "FitLogger.h"
#include "FitUtils.h"
#include
#include
#include
#include
#include
#include
#include