diff --git a/CMakeLists.txt b/CMakeLists.txt
index 552aa22..284f985 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,189 +1,190 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
cmake_minimum_required (VERSION 2.6 FATAL_ERROR)
project(NUISANCE)
include(ExternalProject)
enable_language(Fortran)
set (NUISANCE_VERSION_MAJOR 1)
set (NUISANCE_VERSION_MINOR 0)
set (NUISANCE_VERSION_REVISION 0)
set (NUISANCE_VERSION_STRING "v${NUISANCE_VERSION_MAJOR}r${NUISANCE_VERSION_MINOR}")
if(${NUISANCE_VERSION_REVISION} STRGREATER "0")
set (NUISANCE_VERSION_STRING "${NUISANCE_VERSION_STRING}p${NUISANCE_VERSION_REVISION}")
endif()
#Set this to TRUE to enable build debugging messages
set(BUILD_DEBUG_MSGS TRUE)
include(${CMAKE_SOURCE_DIR}/cmake/cmessage.cmake)
include(${CMAKE_SOURCE_DIR}/cmake/cacheVariables.cmake)
cmessage(STATUS "CMAKE_INSTALL_PREFIX: \"${CMAKE_INSTALL_PREFIX}\"")
cmessage(STATUS "CMAKE_BUILD_TYPE: \"${CMAKE_BUILD_TYPE}\"")
################################################################################
# Check Dependencies
################################################################################
################################## ROOT ######################################
include(${CMAKE_SOURCE_DIR}/cmake/ROOTSetup.cmake)
################################# HEPMC ######################################
include(${CMAKE_SOURCE_DIR}/cmake/HepMC.cmake)
############################ Reweight Engines ################################
include(${CMAKE_SOURCE_DIR}/cmake/ReweightEnginesSetup.cmake)
############################ Other Generators ################################
include(${CMAKE_SOURCE_DIR}/cmake/GiBUUSetup.cmake)
if(USE_NUANCE)
LIST(APPEND EXTRA_CXX_FLAGS -D__NUANCE_ENABLED__)
endif()
-################################# Pythia6 ####################################
+################################# Pythia6/8 ####################################
include(${CMAKE_SOURCE_DIR}/cmake/pythia6Setup.cmake)
+include(${CMAKE_SOURCE_DIR}/cmake/pythia8Setup.cmake)
################################# gperftools ###################################
include(${CMAKE_SOURCE_DIR}/cmake/gperfSetup.cmake)
if(NOT NOTEST)
enable_testing()
endif()
SET(GENERATOR_SUPPORT)
foreach(gen NEUT;NuWro;GENIE;GiBUU;NUANCE)
if(USE_${gen})
SET(GENERATOR_SUPPORT "${GENERATOR_SUPPORT}${gen} ")
endif()
endforeach(gen)
cmessage(STATUS "Generator Input Support: ${GENERATOR_SUPPORT}")
set(MINCODE
Routines
FCN)
set(CORE
MCStudies
Genie
FitBase
InputHandler
Splines
Reweight
Utils
#Devel
Smearceptance
)
LIST(APPEND ALLEXPERIMENTS
ANL
ArgoNeuT
BEBC
BNL
Electron
FNAL
GGM
K2K
MINERvA
MiniBooNE
SciBooNE
T2K)
foreach(exp ${ALLEXPERIMENTS})
if(NOT NO_${exp})
LIST(APPEND EXPERIMENTS_TO_BUILD ${exp})
else()
LIST(REVERSE EXTRA_CXX_FLAGS)
LIST(APPEND EXTRA_CXX_FLAGS -D__NO_${exp}__)
LIST(REVERSE EXTRA_CXX_FLAGS)
endif()
endforeach()
################################## COMPILER ####################################
include(${CMAKE_SOURCE_DIR}/cmake/c++CompilerSetup.cmake)
################################### doxygen ###################################
include(${CMAKE_SOURCE_DIR}/cmake/docsSetup.cmake)
################################################################################
set(MINIMUM_INCLUDE_DIRECTORIES)
LIST(APPEND MINIMUM_INCLUDE_DIRECTORIES
${RWENGINE_INCLUDE_DIRECTORIES}
${CMAKE_SOURCE_DIR}/src/FitBase
${CMAKE_SOURCE_DIR}/src/Reweight
${CMAKE_SOURCE_DIR}/src/InputHandler
${CMAKE_SOURCE_DIR}/src/Splines
${CMAKE_SOURCE_DIR}/src/Utils
${CMAKE_SOURCE_DIR}/src/Genie)
cmessage(DEBUG "Base include directories: ${MINIMUM_INCLUDE_DIRECTORIES}")
set(EXP_INCLUDE_DIRECTORIES)
foreach(edir ${EXPERIMENTS_TO_BUILD})
LIST(APPEND EXP_INCLUDE_DIRECTORIES ${CMAKE_SOURCE_DIR}/src/${edir})
endforeach()
cmessage(DEBUG "Included experiments: ${EXP_INCLUDE_DIRECTORIES}")
foreach(mdir ${MINCODE})
cmessage (DEBUG "Configuring directory: src/${mdir}")
add_subdirectory(src/${mdir})
endforeach()
foreach(edir ${EXPERIMENTS_TO_BUILD})
cmessage (DEBUG "Configuring directory: src/${edir}")
add_subdirectory(src/${edir})
endforeach()
foreach(cdir ${CORE})
cmessage (DEBUG "Configuring directory: src/${cdir}")
add_subdirectory(src/${cdir})
endforeach()
cmessage(DEBUG "Module targets: ${MODULETargets}")
add_subdirectory(app)
add_subdirectory(src/Tests)
configure_file(cmake/setup.sh.in
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.sh" @ONLY)
install(FILES
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.sh" DESTINATION
${CMAKE_INSTALL_PREFIX})
install(PROGRAMS
"${PROJECT_SOURCE_DIR}/scripts/nuiscardgen" DESTINATION
bin)
install(PROGRAMS
"${PROJECT_SOURCE_DIR}/scripts/nuissamples" DESTINATION
bin)
diff --git a/cmake/GENIESetup.cmake b/cmake/GENIESetup.cmake
index 1fdb3dd..148d685 100644
--- a/cmake/GENIESetup.cmake
+++ b/cmake/GENIESetup.cmake
@@ -1,150 +1,155 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
# TODO
# check system for libxml2
# check whether we need the includes
# check if we can use a subset of the GENIE libraries
################################################################################
# Check Dependencies
################################################################################
################################# GENIE ######################################
if(GENIE STREQUAL "")
cmessage(FATAL_ERROR "Variable GENIE is not defined. "
"The location of a pre-built GENIE install must be defined either as"
" $ cmake -DGENIE=/path/to/GENIE or as and environment vairable"
" $ export GENIE=/path/to/GENIE")
endif()
if (BUILD_GEVGEN)
cmessage(STATUS "Building custom gevgen")
LIST(APPEND EXTRA_CXX_FLAGS -D__GEVGEN_ENABLED__)
endif()
# Extract GENIE VERSION
execute_process (COMMAND ${CMAKE_SOURCE_DIR}/cmake/getgenieversion.sh ${GENIE}
OUTPUT_VARIABLE GENIE_VERSION OUTPUT_STRIP_TRAILING_WHITESPACE)
execute_process (COMMAND genie-config
--libs OUTPUT_VARIABLE GENIE_LD_FLAGS_STR OUTPUT_STRIP_TRAILING_WHITESPACE)
execute_process (COMMAND genie-config
--topsrcdir OUTPUT_VARIABLE GENIE_INCLUDES_DIR OUTPUT_STRIP_TRAILING_WHITESPACE)
string(REGEX MATCH "-L\([^ ]+\) \(.*\)$" PARSE_GENIE_LIBS_MATCH ${GENIE_LD_FLAGS_STR})
cmessage(DEBUG "genie-config --libs: ${GENIE_LD_FLAGS_STR}")
if(NOT PARSE_GENIE_LIBS_MATCH)
cmessage(FATAL_ERROR "Expected to be able to parse the result of genie-config --libs to a lib directory and a list of libraries to include, but got: \"${GENIE_LD_FLAGS_STR}\"")
endif()
set(GENIE_LIB_DIR ${CMAKE_MATCH_1})
set(GENIE_LIBS_RAW ${CMAKE_MATCH_2})
string(REPLACE "-l" "" GENIE_LIBS_STRIPED "${GENIE_LIBS_RAW}")
cmessage(STATUS "GENIE version : ${GENIE_VERSION}")
cmessage(STATUS "GENIE libdir : ${GENIE_LIB_DIR}")
cmessage(STATUS "GENIE libs : ${GENIE_LIBS_STRIPED}")
string(REGEX MATCH "ReinSeghal" WASMATCHED ${GENIE_LIBS_STRIPED})
if(WASMATCHED AND GENIE_VERSION STREQUAL "210")
set(GENIE_SEHGAL ${GENIE_LIBS_STRIPED})
STRING(REPLACE "ReinSeghal" "ReinSehgal" GENIE_LIBS_STRIPED ${GENIE_SEHGAL})
cmessage(DEBUG "Fixed inconsistency in library naming: ${GENIE_LIBS_STRIPED}")
endif()
string(REGEX MATCH "ReWeight" WASMATCHED ${GENIE_LIBS_STRIPED})
if(NOT WASMATCHED)
set(GENIE_LIBS_STRIPED "GReWeight ${GENIE_LIBS_STRIPED}")
cmessage(DEBUG "Force added ReWeight library: ${GENIE_LIBS_STRIPED}")
endif()
string(REPLACE " " ";" GENIE_LIBS_LIST "${GENIE_LIBS_STRIPED}")
cmessage(DEBUG "genie-config --libs -- MATCH1: ${CMAKE_MATCH_1}")
cmessage(DEBUG "genie-config --libs -- MATCH2: ${CMAKE_MATCH_2}")
cmessage(DEBUG "genie-config --libs -- libs stripped: ${GENIE_LIBS_STRIPED}")
cmessage(DEBUG "genie-config --libs -- libs list: ${GENIE_LIBS_LIST}")
################################ LHAPDF ######################################
if(LHAPDF_LIB STREQUAL "")
cmessage(FATAL_ERROR "Variable LHAPDF_LIB is not defined. The location of a pre-built lhapdf install must be defined either as $ cmake -DLHAPDF_LIB=/path/to/LHAPDF_libraries or as and environment vairable $ export LHAPDF_LIB=/path/to/LHAPDF_libraries")
endif()
if(LHAPDF_INC STREQUAL "")
cmessage(FATAL_ERROR "Variable LHAPDF_INC is not defined. The location of a pre-built lhapdf install must be defined either as $ cmake -DLHAPDF_INC=/path/to/LHAPDF_includes or as and environment vairable $ export LHAPDF_INC=/path/to/LHAPDF_includes")
endif()
if(LHAPATH STREQUAL "")
cmessage(FATAL_ERROR "Variable LHAPATH is not defined. The location of a the LHAPATH directory must be defined either as $ cmake -DLHAPATH=/path/to/LHAPATH or as and environment variable $ export LHAPATH=/path/to/LHAPATH")
endif()
################################ LIBXML ######################################
if(LIBXML2_LIB STREQUAL "")
cmessage(FATAL_ERROR "Variable LIBXML2_LIB is not defined. The location of a pre-built libxml2 install must be defined either as $ cmake -DLIBXML2_LIB=/path/to/LIBXML2_libraries or as and environment vairable $ export LIBXML2_LIB=/path/to/LIBXML2_libraries")
endif()
if(LIBXML2_INC STREQUAL "")
cmessage(FATAL_ERROR "Variable LIBXML2_INC is not defined. The location of a pre-built libxml2 install must be defined either as $ cmake -DLIBXML2_INC=/path/to/LIBXML2_includes or as and environment vairable $ export LIBXML2_INC=/path/to/LIBXML2_includes")
endif()
############################### log4cpp ######################################
if(LOG4CPP_LIB STREQUAL "")
cmessage(FATAL_ERROR "Variable LOG4CPP_LIB is not defined. The location of a pre-built log4cpp install must be defined either as $ cmake -DLOG4CPP_LIB=/path/to/LOG4CPP_libraries or as and environment vairable $ export LOG4CPP_LIB=/path/to/LOG4CPP_libraries")
endif()
if(LOG4CPP_INC STREQUAL "")
cmessage(FATAL_ERROR "Variable LOG4CPP_INC is not defined. The location of a pre-built log4cpp install must be defined either as $ cmake -DGENIE_LOG4CPP_INC=/path/to/LOG4CPP_includes or as and environment vairable $ export LOG4CPP_INC=/path/to/LOG4CPP_includes")
endif()
################################################################################
LIST(APPEND EXTRA_CXX_FLAGS -D__GENIE_ENABLED__ -D__GENIE_VERSION__=${GENIE_VERSION})
LIST(APPEND RWENGINE_INCLUDE_DIRECTORIES
${GENIE_INCLUDES_DIR}
${GENIE_INCLUDES_DIR}/GHEP
${GENIE_INCLUDES_DIR}/Ntuple
${GENIE_INCLUDES_DIR}/ReWeight
${GENIE_INCLUDES_DIR}/Apps
${GENIE_INCLUDES_DIR}/FluxDrivers
${GENIE_INCLUDES_DIR}/EVGDrivers
${LHAPDF_INC}
${LIBXML2_INC}
${LOG4CPP_INC})
SAYVARS()
LIST(APPEND EXTRA_LINK_DIRS
${GENIE_LIB_DIR}
${LHAPDF_LIB}
${LIBXML2_LIB}
${LOG4CPP_LIB})
#LIST(REVERSE EXTRA_LIBS)
#LIST(REVERSE GENIE_LIBS_LIST)
LIST(APPEND EXTRA_LIBS ${GENIE_LIBS_LIST})
#LIST(REVERSE EXTRA_LIBS)
LIST(APPEND EXTRA_LIBS LHAPDF xml2 log4cpp)
-set(NEED_PYTHIA6 TRUE)
-set(NEED_ROOTPYTHIA6 TRUE)
+if(USE_PYTHIA8)
+ set(NEED_PYTHIA8 TRUE)
+ set(NEED_ROOTPYTHIA8 TRUE)
+else()
+ set(NEED_PYTHIA6 TRUE)
+ set(NEED_ROOTPYTHIA6 TRUE)
+endif()
set(NEED_ROOTEVEGEN TRUE)
SET(USE_GENIE TRUE CACHE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. " FORCE)
diff --git a/cmake/cacheVariables.cmake b/cmake/cacheVariables.cmake
index 453a83e..5e88e9c 100644
--- a/cmake/cacheVariables.cmake
+++ b/cmake/cacheVariables.cmake
@@ -1,202 +1,206 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
function(CheckAndSetDefaultEnv VARNAME DEFAULT CACHETYPE DOCSTRING ENVNAME)
#cmessage(DEBUG "Trying to assign variable ${VARNAME} into the cache.")
if(NOT DEFINED ${VARNAME})
if(DEFINED ENV{${ENVNAME}} AND NOT $ENV{${ENVNAME}} STREQUAL "")
set(${VARNAME} $ENV{${ENVNAME}} CACHE ${CACHETYPE} ${DOCSTRING})
cmessage(DEBUG " Read ${VARNAME} from ENVVAR ${ENVNAME} as $ENV{${ENVNAME}}.")
else()
set(${VARNAME} ${DEFAULT} CACHE ${CACHETYPE} ${DOCSTRING})
endif()
else()
set(${VARNAME} ${${VARNAME}} CACHE ${CACHETYPE} ${DOCSTRING})
unset(${VARNAME})
endif()
cmessage(CACHE "--Set cache variable: \"${VARNAME}\" to \"${${VARNAME}}\", in cache ${CACHETYPE}.")
endfunction()
function(CheckAndSetDefaultCache VARNAME DEFAULT CACHETYPE DOCSTRING)
# cmessage(DEBUG "Trying to assign variable ${VARNAME} into the cache.")
if(NOT DEFINED ${VARNAME})
set(${VARNAME} ${DEFAULT} CACHE ${CACHETYPE} ${DOCSTRING})
else()
set(${VARNAME} ${${VARNAME}} CACHE ${CACHETYPE} ${DOCSTRING})
unset(${VARNAME})
endif()
cmessage(CACHE "--Set cache variable: \"${VARNAME}\" to \"${${VARNAME}}\", in cache ${CACHETYPE}.")
endfunction()
function(CheckAndSetDefault VARNAME DEFAULT)
# cmessage(DEBUG "Trying to assign variable ${VARNAME}.")
if(NOT DEFINED ${VARNAME})
set(${VARNAME} ${DEFAULT} PARENT_SCOPE)
set(${VARNAME} ${DEFAULT})
endif()
cmessage(CACHE "--Set variable: \"${VARNAME}\" to \"${${VARNAME}}\".")
endfunction()
CheckAndSetDefaultCache(VERBOSE TRUE BOOL "Whether to configure loudly.")
set (CMAKE_SKIP_BUILD_RPATH TRUE)
#Changes default install path to be a subdirectory of the build dir.
#Can set build dir at configure time with -DCMAKE_INSTALL_PREFIX=/install/path
if(CMAKE_INSTALL_PREFIX STREQUAL "" OR CMAKE_INSTALL_PREFIX STREQUAL
"/usr/local")
set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}")
elseif(NOT DEFINED CMAKE_INSTALL_PREFIX)
set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}")
endif()
if(CMAKE_BUILD_TYPE STREQUAL "")
set(CMAKE_BUILD_TYPE DEBUG)
elseif(NOT DEFINED CMAKE_BUILD_TYPE)
set(CMAKE_BUILD_TYPE DEBUG)
endif()
CheckAndSetDefaultCache(USE_MINIMIZER TRUE INTERNAL "Whether we are using the ROOT minimization libraries. ")
CheckAndSetDefaultCache(USE_HEPMC FALSE BOOL "Whether to enable HepMC input support. ")
CheckAndSetDefaultEnv(HEPMC "" PATH "Path to HepMC source tree root directory. Overrides environment variable \$HEPMC <>" HEPMC)
CheckAndSetDefaultCache(HEPMC_MOMUNIT "GEV" STRING "HepMC momentum units [MEV|GEV]. ")
CheckAndSetDefaultCache(HEPMC_LENUNIT "CM" STRING "HepMC momentum units [MM|CM]. ")
CheckAndSetDefaultCache(HEPMC_USED_EP FALSE INTERNAL "Whether we built HepMC or not. ")
CheckAndSetDefaultCache(USE_NEUT FALSE BOOL "Whether to enable NEUT (reweight) support. Requires external libraries. ")
CheckAndSetDefaultEnv(NEUT_ROOT "" PATH "Path to NEUT source tree root directory. Overrides environment variable \$NEUT_ROOT <>" NEUT_ROOT)
CheckAndSetDefaultEnv(CERN "" PATH "Path to CERNLIB source tree root directory that NEUT was built against. Overrides environment variable \$CERN <>" CERN)
CheckAndSetDefaultEnv(CERN_LEVEL "" STRING "CERNLIB Library version. Overrides environment variable \$CERN_LEVEL <>" CERN_LEVEL)
CheckAndSetDefaultCache(USE_NuWro FALSE BOOL "Whether to enable NuWro support. ")
CheckAndSetDefaultEnv(NUWRO "" PATH "Path to NuWro source tree root directory. Overrides environment variable \$NUWRO <>" NUWRO)
CheckAndSetDefaultEnv(NUWRO_INC "" PATH "Path to NuWro installed includes directory, needs to contain \"params_all.h\". Overrides environment variable \$NUWRO_INC <>" NUWRO_INC)
CheckAndSetDefaultCache(NUWRO_INPUT_FILE "" FILEPATH "Path to an input NuWro event vector, which can be used to build NuWro i/o libraries. <>")
CheckAndSetDefaultCache(NUWRO_BUILT_FROM_FILE FALSE INTERNAL "Whether the NuWro libraries were built by NUISANCE. ")
CheckAndSetDefaultCache(USE_NuWro_RW FALSE BOOL "Whether to try and build support for NuWro reweighting. ")
CheckAndSetDefaultCache(USE_NuWro_SRW_Event FALSE BOOL "Whether to use cut down NuWro reweight event format. Requires NuWro reweight. ")
CheckAndSetDefaultCache(USE_GENIE FALSE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. ")
CheckAndSetDefaultEnv(GENIE "" PATH "Path to GENIE source tree root directory. Overrides environment variable \$GENIE <>" GENIE)
CheckAndSetDefaultEnv(LHAPDF_LIB "" PATH "Path to pre-built LHAPDF libraries. Overrides environment variable \$LHAPDF_LIB. <>" LHAPDF_LIB)
CheckAndSetDefaultEnv(LHAPDF_INC "" PATH "Path to installed LHAPDF headers. Overrides environment variable \$LHAPDF_INC. <>" LHAPDF_INC)
CheckAndSetDefaultEnv(LHAPATH "" PATH "Path to LHA PDF inputs. Overrides environment variable \$LHAPATH. <>" LHAPATH)
CheckAndSetDefaultEnv(LIBXML2_LIB "" PATH "Path to pre-built LIBXML2 libraries. Overrides environment variable \$LIBXML2_LIB. <>" LIBXML2_LIB)
CheckAndSetDefaultEnv(LIBXML2_INC "" PATH "Path to installed LIBXML2 headers. Overrides environment variable \$LIBXML2_INC. <>" LIBXML2_INC)
CheckAndSetDefaultEnv(LOG4CPP_LIB "" PATH "Path to pre-built LOG4CPP libraries. Overrides environment variable \$LOG4CPP_LIB. <>" LOG4CPP_LIB)
CheckAndSetDefaultEnv(LOG4CPP_INC "" PATH "Path to installed LOG4CPP headers. Overrides environment variable \$LOG4CPP_INC. <>" LOG4CPP_INC)
CheckAndSetDefaultCache(BUILD_GEVGEN FALSE BOOL "Whether to build nuisance_gevgen app.")
CheckAndSetDefaultCache(USE_T2K FALSE BOOL "Whether to enable T2KReWeight support. Requires external libraries. ")
CheckAndSetDefaultEnv(T2KREWEIGHT "" PATH "Path to installed T2KREWEIGHTReWeight. Overrides environment variable \$T2KREWEIGHT. <>" T2KREWEIGHT)
CheckAndSetDefaultCache(USE_NIWG FALSE BOOL "Whether to enable (T2K) NIWG ReWeight support. Requires external libraries. ")
CheckAndSetDefaultEnv(NIWG_ROOT "" PATH "Path to installed NIWGReWeight. Overrides environment variable \$NIWG. <>" NIWG)
CheckAndSetDefaultEnv(PYTHIA6 "" PATH "Path to directory containing libPythia6.so. Overrides environment variable \$PYTHIA6 <>" PYTHIA6)
+CheckAndSetDefaultEnv(PYTHIA8 "" PATH "Path to directory containing libPythia8.so. Overrides environment variable \$PYTHIA8 <>" PYTHIA8)
+CheckAndSetDefaultCache(USE_PYTHIA8 FALSE BOOL "Whether to enable PYTHIA8 event support. ")
+
CheckAndSetDefaultCache(USE_GiBUU TRUE BOOL "Whether to enable GiBUU event support. ")
CheckAndSetDefaultCache(BUILD_GiBUU FALSE BOOL "Whether to build supporting GiBUU event tools along with a patched version of GiBUU. ")
CheckAndSetDefaultCache(USE_NUANCE TRUE BOOL "Whether to enable NUANCE event support. ")
CheckAndSetDefaultCache(USE_PROB3PP FALSE BOOL "Whether to download and compile in Prob3++ support. ")
CheckAndSetDefaultCache(NO_EXTERNAL_UPDATE TRUE BOOL "Whether to perform the update target for external dependencies. ")
CheckAndSetDefaultCache(USE_GPERFTOOLS FALSE BOOL "Whether to compile in google performance tools. ")
CheckAndSetDefault(NEED_PYTHIA6 FALSE)
+CheckAndSetDefault(NEED_PYTHIA8 FALSE)
CheckAndSetDefault(NEED_ROOTEVEGEN FALSE)
CheckAndSetDefault(NEED_ROOTPYTHIA6 FALSE)
CheckAndSetDefaultCache(USE_OMP FALSE BOOL "Whether to enable multicore features (there currently are none...). ")
CheckAndSetDefault(NO_EXPERIMENTS FALSE)
cmessage(STATUS "NO_EXPERIMENTS: ${NO_EXPERIMENTS}")
CheckAndSetDefaultCache(NO_ANL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build ANL samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_ArgoNeuT ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build ArgoNeuT samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_BEBC ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build BEBC samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_BNL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build BNL samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_FNAL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build FNAL samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_GGM ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build GGM samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_K2K ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build K2K samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_MINERvA ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build MINERvA samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_MiniBooNE ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build MiniBooNE samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_T2K ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build T2K samples. <-DNO_EXPERIMENTS=FALSE>")
CheckAndSetDefaultCache(NO_SciBooNE ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build SciBooNE samples. <-DNO_EXPERIMENTS=FALSE>")
function(SAYVARS)
LIST(APPEND VARS
USE_HEPMC
HEPMC
HEPMC_MOMUNIT
HEPMC_LENUNIT
HEPMC_USED_EP
USE_NEUT
NEUT_ROOT
CERN
CERN_LEVEL
USE_NuWro
NUWRO
NUWRO_INC
NUWRO_INPUT_FILE
NUWRO_BUILT_FROM_FILE
USE_GENIE
GENIE
LHAPDF_LIB
LHAPDF_INC
LIBXML2_LIB
LIBXML2_INC
LOG4CPP_LIB
GENIE_LOG4CPP_INC
BUILD_GEVGEN
USE_T2K
USE_NIWG
USE_GiBUU
BUILD_GiBUU
USE_NUANCE
NO_EXTERNAL_UPDATE
USE_GPERFTOOLS
NO_ANL
NO_ArgoNeuT
NO_BEBC
NO_BNL
NO_FNAL
NO_GGM
NO_K2K
NO_MINERvA
NO_MiniBooNE
NO_T2K
NO_SciBooNE)
foreach(v ${VARS})
if(DEFINED ${v})
cmessage(DEBUG "VARIABLE: \"${v}\" = \"${${v}}\"")
endif()
endforeach(v)
endfunction()
diff --git a/cmake/pythia8Setup.cmake b/cmake/pythia8Setup.cmake
new file mode 100644
index 0000000..d46002d
--- /dev/null
+++ b/cmake/pythia8Setup.cmake
@@ -0,0 +1,30 @@
+# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
+
+################################################################################
+# This file is part of NUISANCE.
+#
+# NUISANCE is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# NUISANCE is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with NUISANCE. If not, see .
+################################################################################
+
+if(NEED_PYTHIA8)
+ if(PYTHIA8 STREQUAL "")
+ cmessage(FATAL_ERROR "Variable PYTHIA8 is not defined. This must be set to point to a prebuilt PYTHIA8 instance, please set the \$PYTHIA8 environment variable or configure with -DPYTHIA8=/path/to/pythia8.")
+ endif()
+
+ LIST(APPEND EXTRA_LINK_DIRS ${PYTHIA8})
+
+ LIST(REVERSE EXTRA_LIBS)
+ LIST(APPEND EXTRA_LIBS Pythia8 gfortran)
+ LIST(REVERSE EXTRA_LIBS)
+endif()
diff --git a/data/MINERvA/CC0pi/MINERvA_CC0pi_nue_Q2Ratio_cov.txt b/data/MINERvA/CC0pi/MINERvA_CC0pi_nue_Q2Ratio_cov.txt
new file mode 100644
index 0000000..d0a6c85
--- /dev/null
+++ b/data/MINERvA/CC0pi/MINERvA_CC0pi_nue_Q2Ratio_cov.txt
@@ -0,0 +1,6 @@
+0.039 0.018 0.015 0.01 0.0059 -0.0053
+0.018 0.018 0.013 0.01 0.0091 0.0017
+0.015 0.013 0.02 0.014 0.018 0.013
+0.01 0.01 0.014 0.027 0.031 0.033
+0.0059 0.0091 0.018 0.031 0.1 0.084
+-0.0053 0.0017 0.013 0.033 0.084 0.15
diff --git a/data/MINERvA/CC0pi/MINERvA_CC0pi_nue_Q2Ratio_data.txt b/data/MINERvA/CC0pi/MINERvA_CC0pi_nue_Q2Ratio_data.txt
new file mode 100644
index 0000000..5fca0f9
--- /dev/null
+++ b/data/MINERvA/CC0pi/MINERvA_CC0pi_nue_Q2Ratio_data.txt
@@ -0,0 +1,7 @@
+0 0.72 0.197230829
+0.1 0.77 0.13892444
+0.2 0.93 0.13892444
+0.4 1.09 0.161245155
+0.8 1.27 0.317647603
+1.2 1.16 0.39661064
+2 0 0
diff --git a/data/MINERvA/CCcoh/Enu_nu_cov.csv b/data/MINERvA/CCcoh/Enu_nu_cov.csv
new file mode 100644
index 0000000..d3c67e4
--- /dev/null
+++ b/data/MINERvA/CCcoh/Enu_nu_cov.csv
@@ -0,0 +1,9 @@
+0.000966 0.0013858 -0.0000418 0.0000749 0.0000835 0.0000468 0.0000249 0.0000121 0.0000058
+0.0013858 0.0037297 0.0010112 0.0008641 0.0005755 0.0002781 0.0001472 0.0000721 0.0000347
+-0.0000418 0.0010112 0.0006594 0.0004378 0.0002607 0.0001148 0.0000578 0.0000279 0.000014
+0.0000749 0.0008641 0.0004378 0.0005827 0.0002926 0.0001208 0.0000568 0.000027 0.0000136
+0.0000835 0.0005755 0.0002607 0.0002926 0.0003194 0.0001052 0.0000456 0.0000194 0.0000092
+0.0000468 0.0002781 0.0001148 0.0001208 0.0001052 0.0001256 0.0000411 0.0000135 0.0000044
+0.0000249 0.0001472 0.0000578 0.0000568 0.0000456 0.0000411 0.0000588 0.0000138 0.0000041
+0.0000121 0.0000721 0.0000279 0.000027 0.0000194 0.0000135 0.0000138 0.0000207 0.0000051
+0.0000058 0.0000347 0.000014 0.0000136 0.0000092 0.0000044 0.0000041 0.0000051 0.0000079
diff --git a/data/MINERvA/CCcoh/Enu_nu_data.csv b/data/MINERvA/CCcoh/Enu_nu_data.csv
new file mode 100644
index 0000000..09b96cd
--- /dev/null
+++ b/data/MINERvA/CCcoh/Enu_nu_data.csv
@@ -0,0 +1,10 @@
+1.5 3.952E-39 1.02054E-39
+2 2.485E-39 3.97519E-40
+3 2.736E-39 3.94069E-40
+4 3.096E-39 6.23424E-40
+5 5.037E-39 9.30209E-40
+7 6.54E-39 1.15673E-39
+9 6.321E-39 1.3913E-39
+11 8.205E-39 1.60107E-39
+15 1.233E-38 2.33955E-39
+20 0 0
\ No newline at end of file
diff --git a/data/MINERvA/CCcoh/Enu_nubar_cov.csv b/data/MINERvA/CCcoh/Enu_nubar_cov.csv
new file mode 100644
index 0000000..9927991
--- /dev/null
+++ b/data/MINERvA/CCcoh/Enu_nubar_cov.csv
@@ -0,0 +1,9 @@
+0.00276 0.000561 -0.000195 0.000146 0.000396 0.000638 0.000844 0.001117 0.00086
+0.000561 0.002771 0.00159 0.000965 0.002011 0.003759 0.005241 0.00715 0.005339
+-0.000195 0.00159 0.002332 0.001524 0.002409 0.003629 0.005109 0.007435 0.005552
+0.000146 0.000965 0.001524 0.003053 0.002976 0.002489 0.003197 0.004467 0.003481
+0.000396 0.002011 0.002409 0.002976 0.010258 0.006499 0.007921 0.011924 0.006963
+0.000638 0.003759 0.003629 0.002489 0.006499 0.028482 0.019417 0.02304 0.01485
+0.000844 0.005241 0.005109 0.003197 0.007921 0.019417 0.06033 0.04348 0.021354
+0.001117 0.00715 0.007435 0.004467 0.011924 0.02304 0.04348 0.114784 0.039661
+0.00086 0.005339 0.005552 0.003481 0.006963 0.01485 0.021354 0.039661 0.152195
\ No newline at end of file
diff --git a/data/MINERvA/CCcoh/Enu_nubar_data.csv b/data/MINERvA/CCcoh/Enu_nubar_data.csv
new file mode 100644
index 0000000..0be5d2b
--- /dev/null
+++ b/data/MINERvA/CCcoh/Enu_nubar_data.csv
@@ -0,0 +1,10 @@
+1.5 2.54E-40 5.2543E-40
+2 2.049E-39 5.26972E-40
+3 2.673E-39 4.82637E-40
+4 1.725E-39 5.52784E-40
+5 3.928E-39 1.01278E-39
+7 6.84E-39 1.68804E-39
+9 9.917E-39 2.45641E-39
+11 1.5896E-38 3.38843E-39
+15 9.931E-39 3.90069E-39
+20 0 0
\ No newline at end of file
diff --git a/data/MINERvA/CCcoh/Epi_nu_cov.csv b/data/MINERvA/CCcoh/Epi_nu_cov.csv
new file mode 100644
index 0000000..69541d5
--- /dev/null
+++ b/data/MINERvA/CCcoh/Epi_nu_cov.csv
@@ -0,0 +1,9 @@
+0.000966 0.0013858 -0.0000418 0.0000749 0.0000835 0.0000468 0.0000249 0.0000121 0.0000058
+0.0013858 0.0037297 0.0010112 0.0008641 0.0005755 0.0002781 0.0001472 0.0000721 0.0000347
+-0.0000418 0.0010112 0.0006594 0.0004378 0.0002607 0.0001148 0.0000578 0.0000279 0.000014
+0.0000749 0.0008641 0.0004378 0.0005827 0.0002926 0.0001208 0.0000568 0.000027 0.0000136
+0.0000835 0.0005755 0.0002607 0.0002926 0.0003194 0.0001052 0.0000456 0.0000194 0.0000092
+0.0000468 0.0002781 0.0001148 0.0001208 0.0001052 0.0001256 0.0000411 0.0000135 0.0000044
+0.0000249 0.0001472 0.0000578 0.0000568 0.0000456 0.0000411 0.0000588 0.0000138 0.0000041
+0.0000121 0.0000721 0.0000279 0.000027 0.0000194 0.0000135 0.0000138 0.0000207 0.0000051
+0.0000058 0.0000347 0.000014 0.0000136 0.0000092 0.0000044 0.0000041 0.0000051 0.0000079
\ No newline at end of file
diff --git a/data/MINERvA/CCcoh/Epi_nu_data.csv b/data/MINERvA/CCcoh/Epi_nu_data.csv
new file mode 100644
index 0000000..7917269
--- /dev/null
+++ b/data/MINERvA/CCcoh/Epi_nu_data.csv
@@ -0,0 +1,10 @@
+0 5.67E-40 2.88E-40
+0.25 3.473E-39 5.68E-40
+0.5 1.426E-39 2.34E-40
+0.75 1.508E-39 2.22E-40
+1 1.213E-39 1.61E-40
+1.5 7.14E-40 9.4E-41
+2 4.29E-40 6E-41
+2.5 2.27E-40 3.4E-41
+3.5 1.16E-40 2E-41
+4 0 0
\ No newline at end of file
diff --git a/data/MINERvA/CCcoh/Epi_nubar_cov.csv b/data/MINERvA/CCcoh/Epi_nubar_cov.csv
new file mode 100644
index 0000000..0ffac41
--- /dev/null
+++ b/data/MINERvA/CCcoh/Epi_nubar_cov.csv
@@ -0,0 +1,9 @@
+0.0007803 0.0014033 0.0001656 0.0001398 0.0000948 0.0000417 0.0000181 0.000004 0.0000014
+0.0014033 0.0077102 0.0016934 0.0010592 0.0005185 0.0001968 0.0000902 0.0000224 0.0000075
+0.0001656 0.0016934 0.0011912 0.0006337 0.0002988 0.0001055 0.0000484 0.0000138 0.0000057
+0.0001398 0.0010592 0.0006337 0.0008985 0.0003508 0.0001187 0.0000451 0.0000109 0.0000052
+0.0000948 0.0005185 0.0002988 0.0003508 0.0004187 0.0001104 0.0000365 0.0000066 0.0000028
+0.0000417 0.0001968 0.0001055 0.0001187 0.0001104 0.0001744 0.000041 0.000009 0.0000016
+0.0000181 0.0000902 0.0000484 0.0000451 0.0000365 0.000041 0.0000898 0.000016 0.0000031
+0.000004 0.0000224 0.0000138 0.0000109 0.0000066 0.000009 0.000016 0.0000326 0.0000054
+0.0000014 0.0000075 0.0000057 0.0000052 0.0000028 0.0000016 0.0000031 0.0000054 0.000007
\ No newline at end of file
diff --git a/data/MINERvA/CCcoh/Epi_nubar_data.csv b/data/MINERvA/CCcoh/Epi_nubar_data.csv
new file mode 100644
index 0000000..4db6975
--- /dev/null
+++ b/data/MINERvA/CCcoh/Epi_nubar_data.csv
@@ -0,0 +1,10 @@
+0 4.23E-40 2.79213E-40
+0.25 2.245E-39 8.77707E-40
+0.5 1.342E-39 3.45527E-40
+0.75 1.467E-39 3.00188E-40
+1 1.129E-39 2.04841E-40
+1.5 6.59E-40 1.32019E-40
+2 3.93E-40 9.46414E-41
+2.5 1.82E-40 5.70088E-41
+3.5 7.3E-41 2.62488E-41
+4 0 0
\ No newline at end of file
diff --git a/data/MINERvA/CCcoh/Thpi_nu_cov.csv b/data/MINERvA/CCcoh/Thpi_nu_cov.csv
new file mode 100644
index 0000000..1a831bb
--- /dev/null
+++ b/data/MINERvA/CCcoh/Thpi_nu_cov.csv
@@ -0,0 +1,12 @@
+3.7820000E-06 2.7170000E-06 1.8670000E-06 1.2740000E-06 1.1640000E-06 8.7800000E-07 1.0400000E-06 9.2200000E-07 7.0900000E-07 4.0700000E-07 5.4600000E-07 5.4600000E-07
+2.7170000E-06 4.5150000E-06 2.4610000E-06 1.4840000E-06 1.3900000E-06 1.0720000E-06 1.2670000E-06 1.1370000E-06 8.7300000E-07 5.1500000E-07 6.9800000E-07 7.0000000E-07
+1.8670000E-06 2.4610000E-06 2.6340000E-06 1.3920000E-06 1.1070000E-06 8.7600000E-07 1.0050000E-06 8.8800000E-07 6.8900000E-07 4.0500000E-07 5.0700000E-07 4.8200000E-07
+1.2740000E-06 1.4840000E-06 1.3920000E-06 1.8430000E-06 1.0270000E-06 7.4600000E-07 7.2000000E-07 5.9600000E-07 4.7200000E-07 2.6100000E-07 2.3900000E-07 1.7200000E-07
+1.1640000E-06 1.3900000E-06 1.1070000E-06 1.0270000E-06 1.2880000E-06 7.0900000E-07 6.4300000E-07 5.5200000E-07 4.3400000E-07 2.5300000E-07 2.6300000E-07 2.1600000E-07
+8.7800000E-07 1.0720000E-06 8.7600000E-07 7.4600000E-07 7.0900000E-07 9.3400000E-07 5.7500000E-07 4.3200000E-07 3.2100000E-07 2.0000000E-07 2.0400000E-07 1.5600000E-07
+1.0400000E-06 1.2670000E-06 1.0050000E-06 7.2000000E-07 6.4300000E-07 5.7500000E-07 9.4600000E-07 5.9300000E-07 4.0300000E-07 2.4000000E-07 3.0700000E-07 2.7600000E-07
+9.2200000E-07 1.1370000E-06 8.8800000E-07 5.9600000E-07 5.5200000E-07 4.3200000E-07 5.9300000E-07 7.6900000E-07 4.6300000E-07 2.4300000E-07 2.9800000E-07 2.7800000E-07
+7.0900000E-07 8.7300000E-07 6.8900000E-07 4.7200000E-07 4.3400000E-07 3.2100000E-07 4.0300000E-07 4.6300000E-07 5.1800000E-07 2.6500000E-07 2.4900000E-07 2.2100000E-07
+4.0700000E-07 5.1500000E-07 4.0500000E-07 2.6100000E-07 2.5300000E-07 2.0000000E-07 2.4000000E-07 2.4300000E-07 2.6500000E-07 2.9000000E-07 2.0700000E-07 1.5900000E-07
+5.4600000E-07 6.9800000E-07 5.0700000E-07 2.3900000E-07 2.6300000E-07 2.0400000E-07 3.0700000E-07 2.9800000E-07 2.4900000E-07 2.0700000E-07 3.6400000E-07 2.6600000E-07
+5.4600000E-07 7.0000000E-07 4.8200000E-07 1.7200000E-07 2.1600000E-07 1.5600000E-07 2.7600000E-07 2.7800000E-07 2.2100000E-07 1.5900000E-07 2.6600000E-07 3.2900000E-07
\ No newline at end of file
diff --git a/data/MINERvA/CCcoh/Thpi_nu_data.csv b/data/MINERvA/CCcoh/Thpi_nu_data.csv
new file mode 100644
index 0000000..6c964dc
--- /dev/null
+++ b/data/MINERvA/CCcoh/Thpi_nu_data.csv
@@ -0,0 +1,13 @@
+0 1.14E-40 1.92094E-41
+5 1.33E-40 2.1095E-41
+10 1.15E-40 1.66433E-41
+15 1.09E-40 1.36015E-41
+20 7.3E-41 1.16619E-41
+25 5.4E-41 9.43398E-42
+30 4.2E-41 9.43398E-42
+35 2.7E-41 8.94427E-42
+40 1.8E-41 7.2111E-42
+45 5E-42 5E-42
+50 -5E-42 6.32456E-42
+60 -8E-42 5.38516E-42
+70 0 0
\ No newline at end of file
diff --git a/data/MINERvA/CCcoh/Thpi_nubar_cov.csv b/data/MINERvA/CCcoh/Thpi_nubar_cov.csv
new file mode 100644
index 0000000..fac4f43
--- /dev/null
+++ b/data/MINERvA/CCcoh/Thpi_nubar_cov.csv
@@ -0,0 +1,12 @@
+3.3460000E-06 1.7960000E-06 1.2590000E-06 8.5100000E-07 1.3010000E-06 7.7800000E-07 9.9000000E-07 7.7600000E-07 7.2300000E-07 1.8300000E-07 2.5100000E-07 9.0100000E-07
+1.7960000E-06 4.3450000E-06 2.1650000E-06 1.5990000E-06 2.5320000E-06 1.9160000E-06 2.4200000E-06 1.8270000E-06 9.3500000E-07 2.6600000E-07 3.1000000E-07 1.0490000E-06
+1.2590000E-06 2.1650000E-06 3.0030000E-06 1.4040000E-06 1.8650000E-06 1.4260000E-06 1.8250000E-06 1.4050000E-06 8.7600000E-07 3.1300000E-07 3.4900000E-07 9.6000000E-07
+8.5100000E-07 1.5990000E-06 1.4040000E-06 2.3250000E-06 1.8410000E-06 1.3210000E-06 1.5460000E-06 1.1950000E-06 5.5800000E-07 2.3500000E-07 1.7600000E-07 3.5700000E-07
+1.3010000E-06 2.5320000E-06 1.8650000E-06 1.8410000E-06 3.3490000E-06 2.1040000E-06 2.4450000E-06 1.8890000E-06 8.8800000E-07 3.5600000E-07 3.3700000E-07 8.5800000E-07
+7.7800000E-07 1.9160000E-06 1.4260000E-06 1.3210000E-06 2.1040000E-06 2.1650000E-06 2.0470000E-06 1.4970000E-06 6.0100000E-07 3.1800000E-07 2.7900000E-07 5.5400000E-07
+9.9000000E-07 2.4200000E-06 1.8250000E-06 1.5460000E-06 2.4450000E-06 2.0470000E-06 2.8370000E-06 1.9090000E-06 7.4600000E-07 3.1400000E-07 3.4500000E-07 8.3100000E-07
+7.7600000E-07 1.8270000E-06 1.4050000E-06 1.1950000E-06 1.8890000E-06 1.4970000E-06 1.9090000E-06 1.9270000E-06 8.5000000E-07 3.3500000E-07 2.9700000E-07 6.8900000E-07
+7.2300000E-07 9.3500000E-07 8.7600000E-07 5.5800000E-07 8.8800000E-07 6.0100000E-07 7.4600000E-07 8.5000000E-07 1.0810000E-06 4.9500000E-07 3.1000000E-07 6.6800000E-07
+1.8300000E-07 2.6600000E-07 3.1300000E-07 2.3500000E-07 3.5600000E-07 3.1800000E-07 3.1400000E-07 3.3500000E-07 4.9500000E-07 5.1000000E-07 2.4000000E-07 2.4500000E-07
+2.5100000E-07 3.1000000E-07 3.4900000E-07 1.7600000E-07 3.3700000E-07 2.7900000E-07 3.4500000E-07 2.9700000E-07 3.1000000E-07 2.4000000E-07 3.4000000E-07 3.5000000E-07
+9.0100000E-07 1.0490000E-06 9.6000000E-07 3.5700000E-07 8.5800000E-07 5.5400000E-07 8.3100000E-07 6.8900000E-07 6.6800000E-07 2.4500000E-07 3.5000000E-07 1.1260000E-06
\ No newline at end of file
diff --git a/data/MINERvA/CCcoh/Thpi_nubar_data.csv b/data/MINERvA/CCcoh/Thpi_nubar_data.csv
new file mode 100644
index 0000000..3cdfed1
--- /dev/null
+++ b/data/MINERvA/CCcoh/Thpi_nubar_data.csv
@@ -0,0 +1,13 @@
+0 6.8E-41 1.84391E-41
+5 1.2E-40 2.12132E-41
+10 8.8E-41 1.76918E-41
+15 9.1E-41 1.48661E-41
+20 7.2E-41 1.83848E-41
+25 3.9E-41 1.5E-41
+30 3.1E-41 1.64012E-41
+35 1.9E-41 1.42127E-41
+40 1.1E-41 1.08167E-41
+45 0.1E-62 7.2111E-42
+50 -4E-42 6.40312E-42
+60 -5E-42 1.04403E-41
+70 0 0
\ No newline at end of file
diff --git a/data/MINERvA/CCcoh/Enu_antinu_Covar_flux.csv b/data/MINERvA/CCcoh/callums/Enu_antinu_Covar_flux.csv
similarity index 100%
rename from data/MINERvA/CCcoh/Enu_antinu_Covar_flux.csv
rename to data/MINERvA/CCcoh/callums/Enu_antinu_Covar_flux.csv
diff --git a/data/MINERvA/CCcoh/Enu_antinu_Covar_stat.csv b/data/MINERvA/CCcoh/callums/Enu_antinu_Covar_stat.csv
similarity index 100%
rename from data/MINERvA/CCcoh/Enu_antinu_Covar_stat.csv
rename to data/MINERvA/CCcoh/callums/Enu_antinu_Covar_stat.csv
diff --git a/data/MINERvA/CCcoh/Enu_antinu_Covar_sys.csv b/data/MINERvA/CCcoh/callums/Enu_antinu_Covar_sys.csv
similarity index 100%
rename from data/MINERvA/CCcoh/Enu_antinu_Covar_sys.csv
rename to data/MINERvA/CCcoh/callums/Enu_antinu_Covar_sys.csv
diff --git a/data/MINERvA/CCcoh/Enu_antinu_XSec.csv b/data/MINERvA/CCcoh/callums/Enu_antinu_XSec.csv
similarity index 100%
rename from data/MINERvA/CCcoh/Enu_antinu_XSec.csv
rename to data/MINERvA/CCcoh/callums/Enu_antinu_XSec.csv
diff --git a/data/MINERvA/CCcoh/Enu_nu_Covar_flux.csv b/data/MINERvA/CCcoh/callums/Enu_nu_Covar_flux.csv
similarity index 100%
rename from data/MINERvA/CCcoh/Enu_nu_Covar_flux.csv
rename to data/MINERvA/CCcoh/callums/Enu_nu_Covar_flux.csv
diff --git a/data/MINERvA/CCcoh/Enu_nu_Covar_stat.csv b/data/MINERvA/CCcoh/callums/Enu_nu_Covar_stat.csv
similarity index 100%
rename from data/MINERvA/CCcoh/Enu_nu_Covar_stat.csv
rename to data/MINERvA/CCcoh/callums/Enu_nu_Covar_stat.csv
diff --git a/data/MINERvA/CCcoh/Enu_nu_Covar_sys.csv b/data/MINERvA/CCcoh/callums/Enu_nu_Covar_sys.csv
similarity index 100%
rename from data/MINERvA/CCcoh/Enu_nu_Covar_sys.csv
rename to data/MINERvA/CCcoh/callums/Enu_nu_Covar_sys.csv
diff --git a/data/MINERvA/CCcoh/Enu_nu_XSec.csv b/data/MINERvA/CCcoh/callums/Enu_nu_XSec.csv
similarity index 100%
rename from data/MINERvA/CCcoh/Enu_nu_XSec.csv
rename to data/MINERvA/CCcoh/callums/Enu_nu_XSec.csv
diff --git a/data/MINERvA/CCcoh/Epi_antinu_Covar_flux.csv b/data/MINERvA/CCcoh/callums/Epi_antinu_Covar_flux.csv
similarity index 100%
rename from data/MINERvA/CCcoh/Epi_antinu_Covar_flux.csv
rename to data/MINERvA/CCcoh/callums/Epi_antinu_Covar_flux.csv
diff --git a/data/MINERvA/CCcoh/Epi_antinu_Covar_stat.csv b/data/MINERvA/CCcoh/callums/Epi_antinu_Covar_stat.csv
similarity index 100%
rename from data/MINERvA/CCcoh/Epi_antinu_Covar_stat.csv
rename to data/MINERvA/CCcoh/callums/Epi_antinu_Covar_stat.csv
diff --git a/data/MINERvA/CCcoh/Epi_antinu_Covar_sys.csv b/data/MINERvA/CCcoh/callums/Epi_antinu_Covar_sys.csv
similarity index 100%
rename from data/MINERvA/CCcoh/Epi_antinu_Covar_sys.csv
rename to data/MINERvA/CCcoh/callums/Epi_antinu_Covar_sys.csv
diff --git a/data/MINERvA/CCcoh/callums/Epi_antinu_XSec.csv b/data/MINERvA/CCcoh/callums/Epi_antinu_XSec.csv
new file mode 100755
index 0000000..28e8693
--- /dev/null
+++ b/data/MINERvA/CCcoh/callums/Epi_antinu_XSec.csv
@@ -0,0 +1,10 @@
+0 0.423 0.146 0.238
+0.25 2.245 0.313 0.820
+0.5 1.342 0.142 0.315
+0.75 1.467 0.127 0.272
+1.0 1.129 0.094 0.182
+1.5 0.659 0.073 0.110
+2.0 0.393 0.059 0.074
+2.5 0.182 0.035 0.045
+3.5 0.073 0.017 0.020
+4.0 0.000 0.000 0.000
diff --git a/data/MINERvA/CCcoh/Epi_nu_Covar_flux.csv b/data/MINERvA/CCcoh/callums/Epi_nu_Covar_flux.csv
similarity index 100%
rename from data/MINERvA/CCcoh/Epi_nu_Covar_flux.csv
rename to data/MINERvA/CCcoh/callums/Epi_nu_Covar_flux.csv
diff --git a/data/MINERvA/CCcoh/Epi_nu_Covar_stat.csv b/data/MINERvA/CCcoh/callums/Epi_nu_Covar_stat.csv
similarity index 100%
rename from data/MINERvA/CCcoh/Epi_nu_Covar_stat.csv
rename to data/MINERvA/CCcoh/callums/Epi_nu_Covar_stat.csv
diff --git a/data/MINERvA/CCcoh/Epi_nu_Covar_sys.csv b/data/MINERvA/CCcoh/callums/Epi_nu_Covar_sys.csv
similarity index 100%
rename from data/MINERvA/CCcoh/Epi_nu_Covar_sys.csv
rename to data/MINERvA/CCcoh/callums/Epi_nu_Covar_sys.csv
diff --git a/data/MINERvA/CCcoh/Epi_nu_XSec.csv b/data/MINERvA/CCcoh/callums/Epi_nu_XSec.csv
similarity index 100%
rename from data/MINERvA/CCcoh/Epi_nu_XSec.csv
rename to data/MINERvA/CCcoh/callums/Epi_nu_XSec.csv
diff --git a/data/MINERvA/CCcoh/th_antinu_Covar_flux.csv b/data/MINERvA/CCcoh/callums/th_antinu_Covar_flux.csv
similarity index 100%
rename from data/MINERvA/CCcoh/th_antinu_Covar_flux.csv
rename to data/MINERvA/CCcoh/callums/th_antinu_Covar_flux.csv
diff --git a/data/MINERvA/CCcoh/th_antinu_Covar_stat.csv b/data/MINERvA/CCcoh/callums/th_antinu_Covar_stat.csv
similarity index 100%
rename from data/MINERvA/CCcoh/th_antinu_Covar_stat.csv
rename to data/MINERvA/CCcoh/callums/th_antinu_Covar_stat.csv
diff --git a/data/MINERvA/CCcoh/th_antinu_Covar_sys.csv b/data/MINERvA/CCcoh/callums/th_antinu_Covar_sys.csv
similarity index 100%
rename from data/MINERvA/CCcoh/th_antinu_Covar_sys.csv
rename to data/MINERvA/CCcoh/callums/th_antinu_Covar_sys.csv
diff --git a/data/MINERvA/CCcoh/th_antinu_XSec.csv b/data/MINERvA/CCcoh/callums/th_antinu_XSec.csv
similarity index 100%
rename from data/MINERvA/CCcoh/th_antinu_XSec.csv
rename to data/MINERvA/CCcoh/callums/th_antinu_XSec.csv
diff --git a/data/MINERvA/CCcoh/th_nu_Covar_flux.csv b/data/MINERvA/CCcoh/callums/th_nu_Covar_flux.csv
similarity index 100%
rename from data/MINERvA/CCcoh/th_nu_Covar_flux.csv
rename to data/MINERvA/CCcoh/callums/th_nu_Covar_flux.csv
diff --git a/data/MINERvA/CCcoh/th_nu_Covar_stat.csv b/data/MINERvA/CCcoh/callums/th_nu_Covar_stat.csv
similarity index 100%
rename from data/MINERvA/CCcoh/th_nu_Covar_stat.csv
rename to data/MINERvA/CCcoh/callums/th_nu_Covar_stat.csv
diff --git a/data/MINERvA/CCcoh/th_nu_Covar_sys.csv b/data/MINERvA/CCcoh/callums/th_nu_Covar_sys.csv
similarity index 100%
rename from data/MINERvA/CCcoh/th_nu_Covar_sys.csv
rename to data/MINERvA/CCcoh/callums/th_nu_Covar_sys.csv
diff --git a/data/MINERvA/CCcoh/th_nu_XSec.csv b/data/MINERvA/CCcoh/callums/th_nu_XSec.csv
similarity index 100%
rename from data/MINERvA/CCcoh/th_nu_XSec.csv
rename to data/MINERvA/CCcoh/callums/th_nu_XSec.csv
diff --git a/src/FCN/JointFCN.cxx b/src/FCN/JointFCN.cxx
index 32b2087..3d2e7c0 100755
--- a/src/FCN/JointFCN.cxx
+++ b/src/FCN/JointFCN.cxx
@@ -1,1035 +1,1052 @@
#include "JointFCN.h"
#include
#include "FitUtils.h"
// //***************************************************
// JointFCN::JointFCN(std::string cardfile, TFile* outfile) {
// //***************************************************
// fOutputDir = gDirectory;
// FitPar::Config().out = outfile;
// fCardFile = cardfile;
// LoadSamples(fCardFile);
// fCurIter = 0;
// fMCFilled = false;
// fOutputDir->cd();
// fIterationTree = NULL;
// fDialVals = NULL;
// fNDials = 0;
// fUsingEventManager = FitPar::Config().GetParB("EventManager");
// fOutputDir->cd();
// };
JointFCN::JointFCN(TFile* outfile) {
fOutputDir = gDirectory;
if (outfile) FitPar::Config().out = outfile;
std::vector samplekeys = Config::QueryKeys("sample");
LoadSamples(samplekeys);
std::vector covarkeys = Config::QueryKeys("covar");
LoadPulls(covarkeys);
fCurIter = 0;
fMCFilled = false;
fIterationTree = NULL;
fDialVals = NULL;
fNDials = 0;
fUsingEventManager = FitPar::Config().GetParB("EventManager");
fOutputDir->cd();
}
JointFCN::JointFCN(std::vector samplekeys, TFile* outfile) {
fOutputDir = gDirectory;
if (outfile) FitPar::Config().out = outfile;
LoadSamples(samplekeys);
fCurIter = 0;
fMCFilled = false;
fOutputDir->cd();
fIterationTree = NULL;
fDialVals = NULL;
fNDials = 0;
fUsingEventManager = FitPar::Config().GetParB("EventManager");
fOutputDir->cd();
}
//***************************************************
JointFCN::~JointFCN() {
//***************************************************
// Delete Samples
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
delete exp;
}
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
delete pull;
}
// Sort Tree
if (fIterationTree) DestroyIterationTree();
if (fDialVals) delete fDialVals;
if (fSampleLikes) delete fSampleLikes;
};
//***************************************************
void JointFCN::CreateIterationTree(std::string name, FitWeight* rw) {
//***************************************************
LOG(FIT) << " Creating new iteration tree! " << std::endl;
if (fIterationTree && !name.compare(fIterationTree->GetName())) {
fIterationTree->Reset();
return;
}
fIterationTree = new TTree(name.c_str(), name.c_str());
// Setup Main Branches
fIterationTree->Branch("total_likelihood", &fLikelihood,
"total_likelihood/D");
fIterationTree->Branch("total_ndof", &fNDOF, "total_ndof/I");
// Setup Sample Arrays
int ninputs = fSamples.size() + fPulls.size();
fSampleLikes = new double[ninputs];
fSampleNDOF = new int[ninputs];
fNDOF = GetNDOF();
// Setup Sample Branches
int count = 0;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
std::string name = exp->GetName();
std::string liketag = name + "_likelihood";
std::string ndoftag = name + "_ndof";
fIterationTree->Branch(liketag.c_str(), &fSampleLikes[count],
(liketag + "/D").c_str());
fIterationTree->Branch(ndoftag.c_str(), &fSampleNDOF[count],
(ndoftag + "/D").c_str());
count++;
}
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
std::string name = pull->GetName();
std::string liketag = name + "_likelihood";
std::string ndoftag = name + "_ndof";
fIterationTree->Branch(liketag.c_str(), &fSampleLikes[count],
(liketag + "/D").c_str());
fIterationTree->Branch(ndoftag.c_str(), &fSampleNDOF[count],
(ndoftag + "/D").c_str());
count++;
}
// Add Dial Branches
std::vector dials = rw->GetDialNames();
fNDials = dials.size();
fDialVals = new double[fNDials];
for (int i = 0; i < fNDials; i++) {
fIterationTree->Branch(dials[i].c_str(), &fDialVals[i],
(dials[i] + "/D").c_str());
}
}
//***************************************************
void JointFCN::DestroyIterationTree() {
//***************************************************
if (!fIterationTree) {
delete fIterationTree;
}
}
//***************************************************
void JointFCN::WriteIterationTree() {
//***************************************************
if (!fIterationTree) {
ERR(FTL) << "Can't save empty iteration tree!" << std::endl;
throw;
}
fIterationTree->Write();
}
//***************************************************
void JointFCN::FillIterationTree(FitWeight* rw) {
//***************************************************
if (!fIterationTree) {
ERR(FTL) << "Trying to fill iteration_tree when it is NULL!" << std::endl;
throw;
}
rw->GetAllDials(fDialVals, fNDials);
fIterationTree->Fill();
}
//***************************************************
double JointFCN::DoEval(const double* x) {
//***************************************************
// WEIGHT ENGINE
fDialChanged = FitBase::GetRW()->HasRWDialChanged(x);
FitBase::GetRW()->UpdateWeightEngine(x);
if (fDialChanged) {
FitBase::GetRW()->Reconfigure();
FitBase::EvtManager().ResetWeightFlags();
}
if (LOG_LEVEL(REC)) {
FitBase::GetRW()->Print();
}
// SORT SAMPLES
ReconfigureSamples();
// GET TEST STAT
fLikelihood = GetLikelihood();
// PRINT PROGRESS
LOG(FIT) << "Current Stat (iter. " << this->fCurIter << ") = " << fLikelihood
<< std::endl;
// UPDATE TREE
if (fIterationTree) FillIterationTree(FitBase::GetRW());
return fLikelihood;
}
//***************************************************
int JointFCN::GetNDOF() {
//***************************************************
int totaldof = 0;
int count = 0;
// Total number of Free bins in each MC prediction
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
int dof = exp->GetNDOF();
// Save Seperate DOF
if (fIterationTree) {
fSampleNDOF[count] = dof;
}
// Add to total
totaldof += dof;
count++;
}
// Loop over pulls
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
double dof = pull->GetLikelihood();
// Save seperate DOF
if (fIterationTree) {
fSampleNDOF[count] = dof;
}
// Add to total
totaldof += dof;
count++;
}
// Set Data Variable
fNDOF = totaldof;
return totaldof;
}
//***************************************************
double JointFCN::GetLikelihood() {
//***************************************************
LOG(MIN) << std::left << std::setw(43) << "Getting likelihoods..."
<< " : "
<< "-2logL" << std::endl;
// Loop and add up likelihoods in an uncorrelated way
double like = 0.0;
int count = 0;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
double newlike = exp->GetLikelihood();
int ndof = exp->GetNDOF();
// Save seperate likelihoods
if (fIterationTree) {
fSampleLikes[count] = newlike;
}
LOG(MIN) << "-> " << std::left << std::setw(40) << exp->GetName() << " : "
<< newlike << "/" << ndof << std::endl;
// Add Weight Scaling
// like *= FitBase::GetRW()->GetSampleLikelihoodWeight(exp->GetName());
// Add to total
like += newlike;
count++;
}
// Loop over pulls
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
double newlike = pull->GetLikelihood();
// Save seperate likelihoods
if (fIterationTree) {
fSampleLikes[count] = newlike;
}
// Add to total
like += newlike;
count++;
}
// Set Data Variable
fLikelihood = like;
return like;
};
void JointFCN::LoadSamples(std::vector samplekeys) {
LOG(MIN) << "Loading Samples : " << samplekeys.size() << std::endl;
for (size_t i = 0; i < samplekeys.size(); i++) {
nuiskey key = samplekeys[i];
// Get Sample Options
std::string samplename = key.GetS("name");
std::string samplefile = key.GetS("input");
std::string sampletype = key.GetS("type");
std::string fakeData = "";
LOG(MIN) << "Loading Sample : " << samplename << std::endl;
fOutputDir->cd();
MeasurementBase* NewLoadedSample = SampleUtils::CreateSample(key);
if (!NewLoadedSample) {
ERR(FTL) << "Could not load sample provided: " << samplename << std::endl;
ERR(FTL) << "Check spelling with that in src/FCN/SampleList.cxx"
<< std::endl;
throw;
} else {
fSamples.push_back(NewLoadedSample);
}
}
}
void JointFCN::LoadPulls(std::vector pullkeys) {
for (size_t i = 0; i < pullkeys.size(); i++) {
nuiskey key = pullkeys[i];
std::string pullname = key.GetS("name");
std::string pullfile = key.GetS("input");
std::string pulltype = key.GetS("type");
fOutputDir->cd();
std::cout << "Creating Pull Term : " << std::endl;
sleep(1);
fPulls.push_back(new ParamPull(pullname, pullfile, pulltype));
}
// // Sample Inputs
// if (!samplevect[0].compare("covar") || !samplevect[0].compare("pulls")
// ||
// !samplevect[0].compare("throws")) {
// // Get all inputs
// std::string name = samplevect[1];
// std::string files = samplevect[2];
// std::string type = "DEFAULT";
// if (samplevect.size() > 3) {
// type = samplevect[3];
// } else if (!samplevect[0].compare("pull")) {
// type = "GAUSPULL";
// } else if (!samplevect[0].compare("throws")) {
// type = "GAUSTHROWS";
// }
// // Create Pull Class
// LOG(MIN) << "Loading up pull term: " << name << " << " << files << "
// ("
// << type << ")" << std::endl;
// std::string fakeData = "";
// fOutputDir->cd();
// fPulls.push_back(new ParamPull(name, files, type));
// }
}
// //***************************************************
// void JointFCN::LoadSamples(std::string cardinput)
// //***************************************************
// {
// LOG(MIN) << "Initializing Samples" << std::endl;
// // Read the card file here and load objects
// std::string line;
// std::ifstream card(cardinput.c_str(), ifstream::in);
// // Make sure they are created in correct working DIR
// fOutputDir->cd();
// while (std::getline(card >> std::ws, line, '\n')) {
// // Skip Empties
// if (line.c_str()[0] == '#') continue;
// if (line.empty()) continue;
// // Parse line
// std::vector samplevect = GeneralUtils::ParseToStr(line, "
// ");
// // Sample Inputs
// if (!samplevect[0].compare("sample")) {
// // Get all inputs
// std::string name = samplevect[1];
// std::string files = samplevect[2];
// std::string type = "DEFAULT";
// if (samplevect.size() > 3) type = samplevect[3];
// // Create Sample Class
// LOG(MIN) << "Loading up sample: " << name << " << " << files << " ("
// << type << ")" << std::endl;
// std::string fakeData = "";
// fOutputDir->cd();
// MeasurementBase* NewLoadedSample = SampleUtils::CreateSample(name,
// files, type,
// fakeData, FitBase::GetRW());
// if (!NewLoadedSample) {
// ERR(FTL) << "Could not load sample provided: " << name << std::endl;
// ERR(FTL) << "Check spelling with that in src/FCN/SampleList.cxx"
// << std::endl;
// throw;
// } else {
// fSamples.push_back(NewLoadedSample);
// }
// }
// // Sample Inputs
// if (!samplevect[0].compare("covar") || !samplevect[0].compare("pulls") ||
// !samplevect[0].compare("throws")) {
// // Get all inputs
// std::string name = samplevect[1];
// std::string files = samplevect[2];
// std::string type = "DEFAULT";
// if (samplevect.size() > 3) {
// type = samplevect[3];
// } else if (!samplevect[0].compare("pull")) {
// type = "GAUSPULL";
// } else if (!samplevect[0].compare("throws")) {
// type = "GAUSTHROWS";
// }
// // Create Pull Class
// LOG(MIN) << "Loading up pull term: " << name << " << " << files << " ("
// << type << ")" << std::endl;
// std::string fakeData = "";
// fOutputDir->cd();
// fPulls.push_back(new ParamPull(name, files, type));
// }
// }
// card.close();
// };
//***************************************************
void JointFCN::ReconfigureSamples(bool fullconfig) {
//***************************************************
int starttime = time(NULL);
LOG(REC) << "Starting Reconfigure iter. " << this->fCurIter << std::endl;
// std::cout << fUsingEventManager << " " << fullconfig << " " << fMCFilled <<
// std::endl;
// Event Manager Reconf
if (fUsingEventManager) {
if (!fullconfig and fMCFilled)
ReconfigureFastUsingManager();
else
ReconfigureUsingManager();
} else {
// Loop over all Measurement Classes
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
// If RW Either do signal or full reconfigure.
if (fDialChanged or !fMCFilled or fullconfig) {
if (!fullconfig and fMCFilled)
exp->ReconfigureFast();
else
exp->Reconfigure();
// If RW Not needed just do normalisation
} else {
exp->Renormalise();
}
}
}
// Loop over pulls and update
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
pull->Reconfigure();
}
fMCFilled = true;
LOG(MIN) << "Finished Reconfigure iter. " << fCurIter << " in "
<< time(NULL) - starttime << "s" << std::endl;
fCurIter++;
}
//***************************************************
void JointFCN::ReconfigureSignal() {
//***************************************************
this->ReconfigureSamples(false);
}
//***************************************************
void JointFCN::ReconfigureAllEvents() {
//***************************************************
FitBase::GetRW()->Reconfigure();
FitBase::EvtManager().ResetWeightFlags();
ReconfigureSamples(true);
}
std::vector JointFCN::GetInputList() {
std::vector InputList;
fIsAllSplines = true;
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
std::vector subsamples = exp->GetSubSamples();
for (size_t i = 0; i < subsamples.size(); i++) {
InputHandlerBase* inp = subsamples[i]->GetInput();
if (std::find(InputList.begin(), InputList.end(), inp) ==
InputList.end()) {
if (subsamples[i]->GetInput()->GetType() != kSPLINEPARAMETER)
fIsAllSplines = false;
InputList.push_back(subsamples[i]->GetInput());
}
}
}
return InputList;
}
std::vector JointFCN::GetSubSampleList() {
std::vector SampleList;
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
std::vector subsamples = exp->GetSubSamples();
for (size_t i = 0; i < subsamples.size(); i++) {
SampleList.push_back(subsamples[i]);
}
}
return SampleList;
}
//***************************************************
void JointFCN::ReconfigureUsingManager() {
//***************************************************
// 'Slow' Event Manager Reconfigure
LOG(REC) << "Event Manager Reconfigure" << std::endl;
int timestart = time(NULL);
// Reset all samples
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ResetAll();
}
// If we are siving signal, reset all containers.
bool savesignal = (FitPar::Config().GetParB("SignalReconfigures"));
if (savesignal) {
// Reset all of our event signal vectors
fSignalEventBoxes.clear();
fSignalEventFlags.clear();
fSampleSignalFlags.clear();
fSignalEventSplines.clear();
}
// Make sure we have a list of inputs
if (fInputList.empty()) {
fInputList = GetInputList();
fSubSampleList = GetSubSampleList();
}
// If all inputs are splines make sure the readers are told
// they need to be reconfigured.
std::vector::iterator inp_iter = fInputList.begin();
if (fIsAllSplines) {
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Tell reader in each BaseEvent it needs a Reconfigure next weight calc.
BaseFitEvt* curevent = curinput->FirstBaseEvent();
if (curevent->fSplineRead) {
curevent->fSplineRead->SetNeedsReconfigure(true);
}
}
}
// MAIN INPUT LOOP ====================
int fillcount = 0;
int inputcount = 0;
inp_iter = fInputList.begin();
// Loop over each input in manager
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Get event information
FitEvent* curevent = curinput->FirstNuisanceEvent();
curinput->CreateCache();
int i = 0;
int nevents = curinput->GetNEvents();
int countwidth = nevents / 5;
// Start event loop iterating until we get a NULL pointer.
while (curevent) {
// Get Event Weight
curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
curevent->Weight = curevent->RWWeight * curevent->InputWeight;
double rwweight = curevent->Weight;
// std::cout << "RWWeight = " << curevent->RWWeight << " " <<
// curevent->InputWeight << std::endl;
// Logging
// std::cout << CHECKLOG(1) << std::endl;
if (LOGGING(REC)) {
if (i % countwidth == 0) {
QLOG(REC, curinput->GetName()
<< " : Processed " << i << " events. [M, W] = ["
<< curevent->Mode << ", " << rwweight << "]");
}
}
// Setup flag for if signal found in at least one sample
bool foundsignal = false;
// Create a new signal bitset for this event
std::vector signalbitset(fSubSampleList.size());
// Create a new signal box vector for this event
std::vector signalboxes;
// Start measurement iterator
size_t measitercount = 0;
std::vector::iterator meas_iter =
fSubSampleList.begin();
// Loop over all subsamples (sub in JointMeas)
for (; meas_iter != fSubSampleList.end(); meas_iter++) {
MeasurementBase* curmeas = (*meas_iter);
// Compare input pointers, to current input, skip if not.
// Pointer tells us if it matches without doing ID checks.
if (curinput != curmeas->GetInput()) {
if (savesignal) {
// Set bit to 0 as definitely not signal
signalbitset[measitercount] = 0;
}
// Count up what measurement we are on.
measitercount++;
// Skip sample as input not signal.
continue;
}
// Fill events for matching inputs.
MeasurementVariableBox* box = curmeas->FillVariableBox(curevent);
bool signal = curmeas->isSignal(curevent);
curmeas->SetSignal(signal);
curmeas->FillHistograms(curevent->Weight);
// If its Signal tally up fills
if (signal) {
fillcount++;
}
// If we are saving signal/splines fill the bitset
if (savesignal) {
signalbitset[measitercount] = signal;
}
// If signal save a clone of the event box for use later.
if (savesignal and signal) {
foundsignal = true;
signalboxes.push_back(box->CloneSignalBox());
}
// Keep track of Measurement we are on.
measitercount++;
}
// Once we've filled the measurements, if saving signal
// push back if any sample flagged this event as signal
if (savesignal) {
fSignalEventFlags.push_back(foundsignal);
}
// Save the vector of signal boxes for this event
if (savesignal and foundsignal) {
fSignalEventBoxes.push_back(signalboxes);
fSampleSignalFlags.push_back(signalbitset);
}
// If all inputs are splines we can save the spline coefficients
// for fast in memory reconfigures later.
if (fIsAllSplines and savesignal and foundsignal) {
// Make temp vector to push back with
std::vector coeff;
for (size_t l = 0; l < (UInt_t)curevent->fSplineRead->GetNPar(); l++) {
coeff.push_back(curevent->fSplineCoeff[l]);
}
// Push back to signal event splines. Kept in sync with
// fSignalEventBoxes size.
// int splinecount = fSignalEventSplines.size();
fSignalEventSplines.push_back(coeff);
// if (splinecount % 1000 == 0) {
// std::cout << "Pushed Back Coeff " << splinecount << " : ";
// for (size_t l = 0; l < fSignalEventSplines[splinecount].size(); l++)
// {
// std::cout << " " << fSignalEventSplines[splinecount][l];
// }
// std::cout << std::endl;
// }
}
// Clean up vectors once done with this event
signalboxes.clear();
signalbitset.clear();
// Iterate to the next event.
curevent = curinput->NextNuisanceEvent();
i++;
}
curinput->RemoveCache();
// Keep track of what input we are on.
inputcount++;
}
// End of Event Loop ===============================
// Now event loop is finished loop over all Measurements
// Converting Binned events to XSec Distributions
iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ConvertEventRates();
}
// Print out statements on approximate memory usage for profiling.
LOG(REC) << "Filled " << fillcount << " signal events." << std::endl;
if (savesignal) {
int mem =
( // sizeof(fSignalEventBoxes) +
// fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) +
sizeof(MeasurementVariableBox1D) * fillcount) *
1E-6;
LOG(REC) << " -> Saved " << fillcount
<< " signal boxes for faster access. (~" << mem << " MB)"
<< std::endl;
if (fIsAllSplines and !fSignalEventSplines.empty()) {
int splmem = sizeof(float) * fSignalEventSplines.size() *
fSignalEventSplines[0].size() * 1E-6;
LOG(REC) << " -> Saved " << fillcount << " " << fSignalEventSplines.size()
<< " spline sets into memory. (~" << splmem << " MB)"
<< std::endl;
}
}
LOG(REC) << "Time taken ReconfigureUsingManager() : "
<< time(NULL) - timestart << std::endl;
+ // Check SignalReconfigures works for all samples
+ if (savesignal){
+ double likefull = GetLikelihood();
+ ReconfigureFastUsingManager();
+ double likefast = GetLikelihood();
+
+ if (fabs(likefull - likefast) > 0.0001)
+ {
+ ERROR(FTL,"Fast and Full Likelihoods DIFFER! : " << likefull << " : " << likefast);
+ ERROR(FTL,"This means some samples you are using are not setup to use SignalReconfigures=1");
+ ERROR(FTL,"Please turn OFF signal reconfigures.");
+ throw;
+ } else {
+ LOG(FIT) << "Likelihoods for FULL and FAST match. Will use FAST next time." << std::endl;
+ }
+ }
+
// End of reconfigure
return;
};
//***************************************************
void JointFCN::ReconfigureFastUsingManager() {
//***************************************************
LOG(FIT) << " -> Doing FAST using manager" << std::endl;
// Get Start time for profilling
int timestart = time(NULL);
// Reset all samples
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ResetAll();
}
// Check for saved variables if not do a full reconfigure.
if (fSignalEventFlags.empty()) {
ERR(WRN) << "Signal Flags Empty! Using normal manager." << std::endl;
ReconfigureUsingManager();
return;
}
bool fFillNuisanceEvent =
FitPar::Config().GetParB("FullEventOnSignalReconfigure");
// Setup fast vector iterators.
std::vector::iterator inpsig_iter = fSignalEventFlags.begin();
std::vector >::iterator box_iter =
fSignalEventBoxes.begin();
std::vector >::iterator spline_iter =
fSignalEventSplines.begin();
std::vector >::iterator samsig_iter =
fSampleSignalFlags.begin();
int splinecount = 0;
// Setup stuff for logging
int fillcount = 0;
int nevents = fSignalEventFlags.size();
int countwidth = nevents / 20;
// If All Splines tell splines they need a reconfigure.
std::vector::iterator inp_iter = fInputList.begin();
if (fIsAllSplines) {
LOG(REC) << "All Spline Inputs so using fast spline loop." << std::endl;
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Tell each fSplineRead in BaseFitEvent to reconf next weight calc
BaseFitEvt* curevent = curinput->FirstBaseEvent();
if (curevent->fSplineRead)
curevent->fSplineRead->SetNeedsReconfigure(true);
}
}
// Loop over all possible spline inputs
double* coreeventweights = new double[fSignalEventBoxes.size()];
splinecount = 0;
inp_iter = fInputList.begin();
inpsig_iter = fSignalEventFlags.begin();
spline_iter = fSignalEventSplines.begin();
// Loop over all signal flags
// For each valid signal flag add one to splinecount
// Get Splines from that count and add to weight
// Add splinecount
int sigcount = 0;
splinecount = 0;
// #pragma omp parallel for shared(splinecount,sigcount)
for (uint iinput = 0; iinput < fInputList.size(); iinput++) {
InputHandlerBase* curinput = fInputList[iinput];
BaseFitEvt* curevent = curinput->FirstBaseEvent();
for (int i = 0; i < curinput->GetNEvents(); i++) {
double rwweight = 0.0;
if (fSignalEventFlags[sigcount]) {
// Get Event Info
if (!fIsAllSplines) {
if (fFillNuisanceEvent)
curinput->GetNuisanceEvent(i);
else
curevent = curinput->GetBaseEvent(i);
} else {
curevent->fSplineCoeff = &fSignalEventSplines[splinecount][0];
}
curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
curevent->Weight = curevent->RWWeight * curevent->InputWeight;
rwweight = curevent->Weight;
coreeventweights[splinecount] = rwweight;
if (splinecount % countwidth == 0) {
LOG(REC) << "Processed " << splinecount
<< " event weights. W = " << rwweight << std::endl;
}
// #pragma omp atomic
splinecount++;
}
// #pragma omp atomic
sigcount++;
}
}
LOG(SAM) << "Processed event weights." << std::endl;
// #pragma omp barrier
// Reset Iterators
inpsig_iter = fSignalEventFlags.begin();
spline_iter = fSignalEventSplines.begin();
box_iter = fSignalEventBoxes.begin();
samsig_iter = fSampleSignalFlags.begin();
int nsplineweights = splinecount;
splinecount = 0;
// Start of Fast Event Loop ============================
// Start input iterators
// Loop over number of inputs
for (int ispline = 0; ispline < nsplineweights; ispline++) {
double rwweight = coreeventweights[ispline];
// Get iterators for this event
std::vector::iterator subsamsig_iter = (*samsig_iter).begin();
std::vector::iterator subbox_iter =
(*box_iter).begin();
// Loop over all sub measurements.
std::vector::iterator meas_iter = fSubSampleList.begin();
for (; meas_iter != fSubSampleList.end(); meas_iter++, subsamsig_iter++) {
MeasurementBase* curmeas = (*meas_iter);
// If event flagged as signal for this sample fill from the box.
if (*subsamsig_iter) {
curmeas->SetSignal(true);
curmeas->FillHistogramsFromBox((*subbox_iter), rwweight);
// Move onto next box if there is one.
subbox_iter++;
fillcount++;
}
}
if (ispline % countwidth == 0) {
LOG(REC) << "Filled " << ispline << " sample weights." << std::endl;
}
// Iterate over the main signal event containers.
samsig_iter++;
box_iter++;
spline_iter++;
splinecount++;
}
// End of Fast Event Loop ===================
LOG(SAM) << "Filled sample distributions." << std::endl;
// Now loop over all Measurements
// Convert Binned events
iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ConvertEventRates();
}
// Cleanup coreeventweights
if (fIsAllSplines) {
delete coreeventweights;
}
// Print some reconfigure profiling.
LOG(REC) << "Filled " << fillcount << " signal events." << std::endl;
LOG(REC) << "Time taken ReconfigureFastUsingManager() : "
<< time(NULL) - timestart << std::endl;
}
//***************************************************
void JointFCN::Write() {
//***************************************************
// Loop over individual experiments and call Write
LOG(MIN) << "Writing each of the data classes..." << std::endl;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
exp->Write();
}
// Save Pull Terms
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
pull->Write();
}
if (FitPar::Config().GetParB("EventManager")) {
// Get list of inputs
std::map fInputs =
FitBase::EvtManager().GetInputs();
std::map::const_iterator iterInp;
for (iterInp = fInputs.begin(); iterInp != fInputs.end(); iterInp++) {
InputHandlerBase* input = (iterInp->second);
input->GetFluxHistogram()->Write();
input->GetXSecHistogram()->Write();
input->GetEventHistogram()->Write();
}
}
};
//***************************************************
void JointFCN::SetFakeData(std::string fakeinput) {
//***************************************************
LOG(MIN) << "Setting fake data from " << fakeinput << std::endl;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
exp->SetFakeDataValues(fakeinput);
}
return;
}
//***************************************************
void JointFCN::ThrowDataToy() {
//***************************************************
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
exp->ThrowDataToy();
}
return;
}
diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_1DEe_nue.cxx b/src/MINERvA/MINERvA_CC0pi_XSec_1DEe_nue.cxx
index 2a36858..944672c 100644
--- a/src/MINERvA/MINERvA_CC0pi_XSec_1DEe_nue.cxx
+++ b/src/MINERvA/MINERvA_CC0pi_XSec_1DEe_nue.cxx
@@ -1,102 +1,102 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CC0pi_XSec_1DEe_nue.h"
//********************************************************************
MINERvA_CC0pi_XSec_1DEe_nue::MINERvA_CC0pi_XSec_1DEe_nue(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA CC0pi nue Ee sample. \n" \
"Target: CH \n" \
"Flux: MINERvA Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("E_{e} (GeV)");
fSettings.SetYTitle("d#sigma/dE_{e} (cm^{2}/GeV)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(0.0, 10.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("MINERvA #nu_e CC0#pi");
fSettings.SetDataInput( FitPar::GetDataBase() + "/MINERvA/CC0pi/MINERvA_CC0pi_nue_Data_ARX1509_05729.root;Data_1DEe" );
fSettings.SetCovarInput( FitPar::GetDataBase() + "/MINERvA/CC0pi/MINERvA_CC0pi_nue_Data_ARX1509_05729.root;Covar_1DEe" );
fSettings.DefineAllowedSpecies("nue,nueb");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / TotalIntegratedFlux();
// Plot Setup -------------------------------------------------------
SetDataFromRootFile( fSettings.GetDataInput() );
SetCovarFromRootFile(fSettings.GetCovarInput() );
- ScaleCovar(1.0 / 100.0);
+ ScaleCovar(1.0 / 1000.0);
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MINERvA_CC0pi_XSec_1DEe_nue::FillEventVariables(FitEvent *event) {
//********************************************************************
int PDGnu = event->GetNeutrinoIn()->fPID;
int PDGe = 0;
if (PDGnu == 12) PDGe = 11;
else if (PDGnu == -12) PDGe = -11;
if (event->NumFSParticle(PDGe) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pe = event->GetHMFSParticle(PDGe)->fP;
Thetae = Pnu.Vect().Angle(Pe.Vect());
Q2QEe = FitUtils::Q2QErec(Pe, cos(Thetae), 34., true);
Ee = Pe.E() / 1000.0;
fXVar = Ee;
return;
}
//********************************************************************
bool MINERvA_CC0pi_XSec_1DEe_nue::isSignal(FitEvent *event) {
//*******************************************************************
// Check that this is a nue CC0pi event
if (!(SignalDef::isCC0pi(event, 12, EnuMin, EnuMax)) and
!(SignalDef::isCC0pi(event, -12, EnuMin, EnuMax))) return false;
// Electron Enrgy
if (Ee < 0.5) return false;
return true;
};
diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_1DQ2_nue.cxx b/src/MINERvA/MINERvA_CC0pi_XSec_1DQ2_nue.cxx
index 6846f22..0264ddc 100644
--- a/src/MINERvA/MINERvA_CC0pi_XSec_1DQ2_nue.cxx
+++ b/src/MINERvA/MINERvA_CC0pi_XSec_1DQ2_nue.cxx
@@ -1,103 +1,103 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CC0pi_XSec_1DQ2_nue.h"
//********************************************************************
MINERvA_CC0pi_XSec_1DQ2_nue::MINERvA_CC0pi_XSec_1DQ2_nue(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CC0pi_XSec_1DQ2_nue sample. \n" \
"Target: CH \n" \
"Flux: MINERvA Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("Q^{2}_{e} (GeV)");
fSettings.SetYTitle("d#sigma/dQ^{2}_{e} (cm^{2}/GeV)^{2}");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(0.0, 10.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("MINERvA_CC0pi_XSec_1DQ2_nue");
fSettings.SetDataInput( FitPar::GetDataBase() + "/MINERvA/CC0pi/MINERvA_CC0pi_nue_Data_ARX1509_05729.root;Data_1DQ2" );
fSettings.SetCovarInput( FitPar::GetDataBase() + "/MINERvA/CC0pi/MINERvA_CC0pi_nue_Data_ARX1509_05729.root;Covar_1DQ2" );
fSettings.DefineAllowedSpecies("nue,nueb");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / TotalIntegratedFlux();
// Plot Setup -------------------------------------------------------
SetDataFromRootFile( fSettings.GetDataInput() );
SetCovarFromRootFile(fSettings.GetCovarInput() );
- ScaleCovar(1.0 / 100.0);
+ ScaleCovar(1.0 / 1000.0);
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MINERvA_CC0pi_XSec_1DQ2_nue::FillEventVariables(FitEvent *event){
//********************************************************************
int PDGnu = event->GetNeutrinoIn()->fPID;
int PDGe = 0;
if (PDGnu == 12) PDGe= 11;
else if (PDGnu == -12) PDGe = -11;
if (event->NumFSParticle(PDGe) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pe = event->GetHMFSParticle(PDGe)->fP;
Thetae = Pnu.Vect().Angle(Pe.Vect());
Q2QEe = FitUtils::Q2QErec(Pe, cos(Thetae), 34., true);
Ee = Pe.E() / 1000.0;
fXVar = Q2QEe;
return;
}
//********************************************************************
bool MINERvA_CC0pi_XSec_1DQ2_nue::isSignal(FitEvent *event){
//*******************************************************************
// Check that this is a nue CC0pi event
if (!(SignalDef::isCC0pi(event, 12, EnuMin, EnuMax)) and
!(SignalDef::isCC0pi(event, -12, EnuMin, EnuMax))) return false;
// Restrct Ee
if (Ee < 0.5) return false;
return true;
};
diff --git a/src/MINERvA/MINERvA_CC0pi_XSec_1DThetae_nue.cxx b/src/MINERvA/MINERvA_CC0pi_XSec_1DThetae_nue.cxx
index 468fe7b..416273e 100644
--- a/src/MINERvA/MINERvA_CC0pi_XSec_1DThetae_nue.cxx
+++ b/src/MINERvA/MINERvA_CC0pi_XSec_1DThetae_nue.cxx
@@ -1,104 +1,104 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CC0pi_XSec_1DThetae_nue.h"
//********************************************************************
MINERvA_CC0pi_XSec_1DThetae_nue::MINERvA_CC0pi_XSec_1DThetae_nue(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CC0pi_XSec_1DThetae_nue sample. \n" \
"Target: CH \n" \
"Flux: MINERvA Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#theta_{e}");
fSettings.SetYTitle("d#sigma/d#theta_{e} (cm^{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(0.0, 10.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("MINERvA #nu_e CC0#pi");
fSettings.SetDataInput( FitPar::GetDataBase() + "/MINERvA/CC0pi/MINERvA_CC0pi_nue_Data_ARX1509_05729.root;Data_1DThetae" );
fSettings.SetCovarInput( FitPar::GetDataBase() + "/MINERvA/CC0pi/MINERvA_CC0pi_nue_Data_ARX1509_05729.root;Covar_1DThetae" );
fSettings.DefineAllowedSpecies("nue,nueb");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) / TotalIntegratedFlux();
// Plot Setup -------------------------------------------------------
SetDataFromRootFile( fSettings.GetDataInput() );
SetCovarFromRootFile(fSettings.GetCovarInput() );
- ScaleCovar(1.0 / 100.0);
+ ScaleCovar(1.0 / 1000.0);
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MINERvA_CC0pi_XSec_1DThetae_nue::FillEventVariables(FitEvent *event){
//********************************************************************
int PDGnu = event->GetNeutrinoIn()->fPID;
int PDGe = 0;
if (PDGnu == 12) PDGe= 11;
else if (PDGnu == -12) PDGe = -11;
if (event->NumFSParticle(PDGe) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pe = event->GetHMFSParticle(PDGe)->fP;
Thetae = Pnu.Vect().Angle(Pe.Vect());
Q2QEe = FitUtils::Q2QErec(Pe, cos(Thetae), 34., true);
Ee = Pe.E() / 1000.0;
fXVar = Thetae * 180. / TMath::Pi();
return;
}
//********************************************************************
bool MINERvA_CC0pi_XSec_1DThetae_nue::isSignal(FitEvent *event){
//*******************************************************************
// Check this is a nue CC0pi event
if (!(SignalDef::isCC0pi(event, 12, EnuMin, EnuMax)) and
!(SignalDef::isCC0pi(event, -12, EnuMin, EnuMax))) return false;
// Restrict EE
if (Ee < 0.5) return false;
return true;
};
diff --git a/src/MINERvA/MINERvA_CC1pip_XSec_1D_2017Update.cxx b/src/MINERvA/MINERvA_CC1pip_XSec_1D_2017Update.cxx
index c90aa64..ed7ca53 100644
--- a/src/MINERvA/MINERvA_CC1pip_XSec_1D_2017Update.cxx
+++ b/src/MINERvA/MINERvA_CC1pip_XSec_1D_2017Update.cxx
@@ -1,181 +1,169 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CC1pip_XSec_1D_2017Update.h"
//********************************************************************
void MINERvA_CC1pip_XSec_1D_2017Update::SetupDataSettings(){
//********************************************************************
// Set Distribution
std::string name = fSettings.GetS("name");
- std::cout << "Parsing NAME " << name << std::endl;
- sleep(5);
if (!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_2017")) fDist = kTpi;
else if (!name.compare("MINERvA_CC1pip_XSec_1Dth_nu_2017")) fDist= kth;
else if (!name.compare("MINERvA_CC1pip_XSec_1Dpmu_nu_2017")) fDist= kpmu;
else if (!name.compare("MINERvA_CC1pip_XSec_1Dthmu_nu_2017")) fDist= kthmu;
else if (!name.compare("MINERvA_CC1pip_XSec_1DQ2_nu_2017")) fDist= kQ2;
else if (!name.compare("MINERvA_CC1pip_XSec_1DEnu_nu_2017")) fDist= kEnu;
// Define what files to use from the dist
- /*
- ls [stowell@hep22 MINERvA]$ ls ../../data/MINERvA/CC1pip/070717/
-cc1pip_updated_1DEnu_data.txt cc1pip_updated_1Dpmu_shapecov.txt cc1pip_updated_1Dthmu_data.txt cc1pip_updated_1DTpi_data.txt
-cc1pip_updated_1DEnu_ratecov.txt cc1pip_updated_1DQ2_data.txt cc1pip_updated_1Dthmu_ratecov.txt cc1pip_updated_1DTpi_ratecov.txt
-cc1pip_updated_1DEnu_shapecov.txt cc1pip_updated_1DQ2_ratecov.txt cc1pip_updated_1Dthmu_shapecov.txt cc1pip_updated_1DTpi_shapecov.txt
-cc1pip_updated_1Dpmu_data.txt cc1pip_updated_1DQ2_shapecov.txt cc1pip_updated_1Dth_ratecov.txt README
-cc1pip_updated_1Dpmu_ratecov.txt cc1pip_updated_1Dth_data.txt cc1pip_updated_1Dth_shapecov.txt
- */
std::string datafile = "";
std::string covarfile = "";
std::string titles = "";
std::string distdescript = "";
switch(fDist){
case (kTpi):
datafile = "cc1pip_updated_1DTpi";
covarfile = "cc1pip_updated_1DTpi";
titles = "CC1#pi Updated;T_{#pi} (MeV);d#sigma/dT_{#pi} (cm^{2}/nucleon/MeV)";
break;
case (kth):
datafile = "cc1pip_updated_1Dth";
covarfile = "cc1pip_updated_1Dth";
titles = "CC1#pi Updated;#theta_{#pi};d#sigma/d#theta_{#pi} (cm^{2}/nucleon)";
break;
case (kpmu):
datafile = "cc1pip_updated_1Dpmu";
covarfile = "cc1pip_updated_1Dpmu";
titles = "CC1#pi Updated;p_{#mu} (GeV);d#sigma/dp_{#mu} (cm^{2}/nucleon/GeV)";
break;
case (kthmu):
datafile = "cc1pip_updated_1Dthmu";
covarfile = "cc1pip_updated_1Dthmu";
titles ="CC1#pi Updated;#theta_{#mu};d#sigma/d#theta_{#mu} (cm^{2}/nucleon)";
break;
case (kQ2):
datafile = "cc1pip_updated_1DQ2";
covarfile = "cc1pip_updated_1DQ2";
titles ="CC1#pi Updated;Q^{2} (GeV^{2});d#sigma/dQ^{2} (cm^{2}/nucleon/GeV^{2})";
break;
case (kEnu):
datafile = "cc1pip_updated_1DEnu";
covarfile = "cc1pip_updated_1DEnu";
titles ="CC1#pi Updated;E_{#nu} (GeV);#sigma(E_#nu) (cm^{2}/nucleon)";
break;
default:
THROW("Unknown Analysis Distribution : " << fDist);
}
// Choose shape or rate covariance
fIsShape = fSettings.Found("type","SHAPE");
std::string covid = fIsShape ? "_shapecov.txt" : "_ratecov.txt";
// Now setup each data distribution and description.
std::string descrip = distdescript + \
"Target: CH \n" \
"Flux: MINERvA Forward Horn Current numu ONLY \n" \
"Signal: Any event with 1 muon, and 1pi+ or 1pi- in FS. W < 1.4";
fSettings.SetDescription(descrip);
fSettings.SetDataInput( GeneralUtils::GetTopLevelDir()+"/data/MINERvA/CC1pip/070717/" + datafile + "_data.txt" );
fSettings.SetCovarInput( GeneralUtils::GetTopLevelDir()+"/data/MINERvA/CC1pip/070717/" + covarfile + covid );
fSettings.SetTitle( GeneralUtils::ParseToStr(titles,";")[0] );
fSettings.SetXTitle( GeneralUtils::ParseToStr(titles,";")[1] );
fSettings.SetYTitle( GeneralUtils::ParseToStr(titles,";")[2] );
return;
}
//********************************************************************
MINERvA_CC1pip_XSec_1D_2017Update::MINERvA_CC1pip_XSec_1D_2017Update(nuiskey samplekey) {
//********************************************************************
// Define Sample Settings common to all data distributions
fSettings = LoadSampleSettings(samplekey);
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(1.5, 10.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
SetupDataSettings();
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// If Enu setup scale factor for Enu Unfolded, otherwise use differential
if (fDist == kEnu) fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents);
else fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents) / TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetCorrelationFromTextFile( fSettings.GetCovarInput() );
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MINERvA_CC1pip_XSec_1D_2017Update::FillEventVariables(FitEvent *event) {
//********************************************************************
fXVar = -999.9;
if (event->NumFSParticle(PhysConst::pdg_charged_pions) == 0 ||
event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Ppip = event->GetHMFSParticle(PhysConst::pdg_charged_pions)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
- double Tpi = FitUtils::T(Ppip); // MeV
+ double Tpi = Ppip.E() - Ppip.Mag();
double th = (180./M_PI)*FitUtils::th(Pnu, Ppip);
double pmu = Pmu.Vect().Mag()/1.E3; // GeV
double thmu = (180.0/M_PI)*FitUtils::th(Pnu, Pmu);
double Q2 = fabs((Pmu - Pnu).Mag2()) / 1.E6; // Using true here?
double Enu = Pnu.E() / 1.E3;
- double hadMass = FitUtils::Wrec(Pnu, Pmu);
- if (hadMass < 1400){
- switch(fDist){
- case kTpi: fXVar = Tpi; break;
- case kth: fXVar = th; break;
- case kpmu: fXVar = pmu; break;
- case kthmu: fXVar = thmu; break;
- case kQ2: fXVar = Q2; break;
- case kEnu: fXVar = Enu; break;
- default:
- THROW("DIST NOT FOUND : " << fDist);
- }
+
+ switch(fDist){
+ case kTpi: fXVar = Tpi; break;
+ case kth: fXVar = th; break;
+ case kpmu: fXVar = pmu; break;
+ case kthmu: fXVar = thmu; break;
+ case kQ2: fXVar = Q2; break;
+ case kEnu: fXVar = Enu; break;
+ default:
+ THROW("DIST NOT FOUND : " << fDist);
}
return;
};
//********************************************************************
bool MINERvA_CC1pip_XSec_1D_2017Update::isSignal(FitEvent *event) {
//********************************************************************
// Only seem to release full phase space
- return SignalDef::isCC1pip_MINERvA(event, EnuMin, EnuMax, false);
+ return SignalDef::isCC1pip_MINERvA_2017(event, EnuMin, EnuMax);
}
diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEnu_antinu.cxx b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEnu_antinu.cxx
index 2ee44ec..d31528d 100644
--- a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEnu_antinu.cxx
+++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEnu_antinu.cxx
@@ -1,87 +1,77 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CCCOHPI_XSec_1DEnu_antinu.h"
//********************************************************************
MINERvA_CCCOHPI_XSec_1DEnu_antinu::MINERvA_CCCOHPI_XSec_1DEnu_antinu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CCCOHPI_XSec_1DEnu_antinu sample. \n" \
"Target: CH \n" \
"Flux: MINERvA Reverse Horn Current numu \n" \
"Signal: Any event with 1 mu+, 1pi-, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("E_{#nu} (GeV)");
fSettings.SetYTitle("d#sigma/dE_{#nu} (cm^{2}/GeV/C^{12})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(1.5, 20.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numub");
fSettings.SetTitle("MINERvA_CCCOHPI_XSec_1DEnu_antinu");
- fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Enu_antinu_XSec.csv");
- fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Enu_antinu_Covar_stat.csv;" +
- GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Enu_antinu_Covar_flux.csv;" +
- GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Enu_antinu_Covar_sys.csv");
+ fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Enu_nubar_data.csv");
+ fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Enu_nubar_cov.csv");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents)*13;
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetCovarFromMultipleTextFiles(fSettings.GetCovarInput());
- // Apply scalings based on the data release
- ScaleData(1E-39);
-
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void MINERvA_CCCOHPI_XSec_1DEnu_antinu::FillEventVariables(FitEvent *event) {
fXVar = event->GetNeutrinoIn()->fP.E()/1000.;
return;
};
//********************************************************************
bool MINERvA_CCCOHPI_XSec_1DEnu_antinu::isSignal(FitEvent *event) {
//********************************************************************
// Signal: numu + 12C --> mu- + pi+ + 12C'
return SignalDef::isCCCOH(event, -14, -211, EnuMin, EnuMax);
}
-
-double MINERvA_CCCOHPI_XSec_1DEnu_antinu::GetLikelihood(){
- double chi2 = StatUtils::GetChi2FromCov(this->fDataHist, this->fMCHist, this->covar, NULL, 1E39, 10);
- return chi2;
-}
diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEnu_antinu.h b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEnu_antinu.h
index b075690..d7e5d81 100644
--- a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEnu_antinu.h
+++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEnu_antinu.h
@@ -1,37 +1,36 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef MINERVA_CCCOHPI_XSEC_1DENU_ANTINU_H_SEEN
#define MINERVA_CCCOHPI_XSEC_1DENU_ANTINU_H_SEEN
#include "Measurement1D.h"
class MINERvA_CCCOHPI_XSec_1DEnu_antinu : public Measurement1D {
public:
MINERvA_CCCOHPI_XSec_1DEnu_antinu(nuiskey samplekey);
virtual ~MINERvA_CCCOHPI_XSec_1DEnu_antinu() {};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
- double GetLikelihood();
private:
};
#endif
diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEnu_nu.cxx b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEnu_nu.cxx
index f38eb52..4165ef6 100644
--- a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEnu_nu.cxx
+++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEnu_nu.cxx
@@ -1,87 +1,76 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CCCOHPI_XSec_1DEnu_nu.h"
//********************************************************************
MINERvA_CCCOHPI_XSec_1DEnu_nu::MINERvA_CCCOHPI_XSec_1DEnu_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CCCOHPI_XSec_1DEnu_nu sample. \n" \
"Target: CH \n" \
"Flux: MINERvA Forward Horn Current numu \n" \
"Signal: Any event with 1 mu-, 1pi+, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("E_{#nu} (GeV)");
fSettings.SetYTitle("d#sigma/dE_{#nu} (cm^{2}/GeV/C^{12})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(1.5, 20.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
fSettings.SetTitle("MINERvA_CCCOHPI_XSec_1DEnu_nu");
- fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Enu_nu_XSec.csv");
- fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Enu_nu_Covar_stat.csv;" +
- GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Enu_nu_Covar_flux.csv;" +
- GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Enu_nu_Covar_sys.csv");
+ fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Enu_nu_data.csv");
+ fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Enu_nu_cov.csv");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents)*13;
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetCovarFromMultipleTextFiles(fSettings.GetCovarInput());
- // Apply scalings based on the data release
- ScaleData(1E-39);
-
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void MINERvA_CCCOHPI_XSec_1DEnu_nu::FillEventVariables(FitEvent *event) {
-
fXVar = event->GetNeutrinoIn()->fP.E()/1000.;
return;
};
//********************************************************************
bool MINERvA_CCCOHPI_XSec_1DEnu_nu::isSignal(FitEvent *event) {
//********************************************************************
// Signal: numu + 12C --> mu- + pi+ + 12C'
return SignalDef::isCCCOH(event, 14, 211, EnuMin, EnuMax);
}
-
-double MINERvA_CCCOHPI_XSec_1DEnu_nu::GetLikelihood(){
- double chi2 = StatUtils::GetChi2FromCov(this->fDataHist, this->fMCHist, this->covar, NULL, 1E39, 10);
- return chi2;
-}
diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEnu_nu.h b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEnu_nu.h
index 647e5d5..8619992 100644
--- a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEnu_nu.h
+++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEnu_nu.h
@@ -1,37 +1,36 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef MINERVA_CCCOHPI_XSEC_1DENU_NU_H_SEEN
#define MINERVA_CCCOHPI_XSEC_1DENU_NU_H_SEEN
#include "Measurement1D.h"
class MINERvA_CCCOHPI_XSec_1DEnu_nu : public Measurement1D {
public:
MINERvA_CCCOHPI_XSec_1DEnu_nu(nuiskey samplekey);
virtual ~MINERvA_CCCOHPI_XSec_1DEnu_nu() {};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
- double GetLikelihood();
private:
};
#endif
diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEpi_antinu.cxx b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEpi_antinu.cxx
index 6baf075..bae166a 100644
--- a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEpi_antinu.cxx
+++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEpi_antinu.cxx
@@ -1,91 +1,81 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CCCOHPI_XSec_1DEpi_antinu.h"
//********************************************************************
MINERvA_CCCOHPI_XSec_1DEpi_antinu::MINERvA_CCCOHPI_XSec_1DEpi_antinu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CCCOHPI_XSec_1DEpi_antinu sample. \n" \
"Target: CH \n" \
"Flux: MINERvA Reverse Horn Current numu \n" \
"Signal: Any event with 1 mu+, 1pi-, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("E_{#pi} (GeV)");
fSettings.SetYTitle("d#sigma/dE_{#pi} (cm^{2}/GeV/C^{12})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(1.5, 20.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numub");
fSettings.SetTitle("MINERvA_CCCOHPI_XSec_1DEpi_antinu");
- fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Epi_antinu_XSec.csv");
- fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Epi_antinu_Covar_stat.csv;" +
- GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Epi_antinu_Covar_flux.csv;" +
- GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Epi_antinu_Covar_sys.csv");
+ fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Epi_nubar_data.csv");
+ fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Epi_nubar_cov.csv");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents)/TotalIntegratedFlux("width")*13;
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetCovarFromMultipleTextFiles(fSettings.GetCovarInput());
- // Apply scalings based on the data release
- ScaleData(1E-39);
-
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void MINERvA_CCCOHPI_XSec_1DEpi_antinu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(-211) == 0) return;
TLorentzVector Ppi = event->GetHMFSParticle(-211)->fP;
fXVar = Ppi.E()/1000.;
return;
};
//********************************************************************
bool MINERvA_CCCOHPI_XSec_1DEpi_antinu::isSignal(FitEvent *event) {
//********************************************************************
// Signal: numu + 12C --> mu- + pi+ + 12C'
return SignalDef::isCCCOH(event, -14, -211, EnuMin, EnuMax);
}
-
-double MINERvA_CCCOHPI_XSec_1DEpi_antinu::GetLikelihood(){
- double chi2 = StatUtils::GetChi2FromCov(this->fDataHist, this->fMCHist, this->covar, NULL, 1E39, 10);
- return chi2;
-}
diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEpi_antinu.h b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEpi_antinu.h
index fbcee07..6f45792 100644
--- a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEpi_antinu.h
+++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEpi_antinu.h
@@ -1,37 +1,36 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef MINERVA_CCCOHPI_XSEC_1DEPI_ANTINU_H_SEEN
#define MINERVA_CCCOHPI_XSEC_1DEPI_ANTINU_H_SEEN
#include "Measurement1D.h"
class MINERvA_CCCOHPI_XSec_1DEpi_antinu : public Measurement1D {
public:
MINERvA_CCCOHPI_XSec_1DEpi_antinu(nuiskey samplekey);
virtual ~MINERvA_CCCOHPI_XSec_1DEpi_antinu() {};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
- double GetLikelihood();
private:
};
#endif
diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEpi_nu.cxx b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEpi_nu.cxx
index 6701f30..575659b 100644
--- a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEpi_nu.cxx
+++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEpi_nu.cxx
@@ -1,91 +1,81 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CCCOHPI_XSec_1DEpi_nu.h"
//********************************************************************
MINERvA_CCCOHPI_XSec_1DEpi_nu::MINERvA_CCCOHPI_XSec_1DEpi_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CCCOHPI_XSec_1DEpi_nu sample. \n" \
"Target: CH \n" \
"Flux: MINERvA Forward Horn Current numu \n" \
"Signal: Any event with 1 mu-, 1pi+, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("E_{#pi} (GeV)");
fSettings.SetYTitle("d#sigma/dE_{#pi} (cm^{2}/GeV/C^{12})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(1.5, 20.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
fSettings.SetTitle("MINERvA_CCCOHPI_XSec_1DEpi_nu");
- fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Epi_nu_XSec.csv");
- fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Epi_nu_Covar_stat.csv;" +
- GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Epi_nu_Covar_flux.csv;" +
- GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Epi_nu_Covar_sys.csv");
+ fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Epi_nu_data.csv");
+ fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Epi_nu_cov.csv");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents)/TotalIntegratedFlux("width")*13;
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetCovarFromMultipleTextFiles(fSettings.GetCovarInput());
- // Apply scalings based on the data release
- ScaleData(1E-39);
-
// Final setup ---------------------------------------------------
FinaliseMeasurement();
-
+ std::cout << "MINERvA_CCCOHPI_XSec_1DEpi_nu.cxx : Data Integral = " << fDataHist->Integral() << std::endl;
};
void MINERvA_CCCOHPI_XSec_1DEpi_nu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(211) == 0) return;
TLorentzVector Ppi = event->GetHMFSParticle(211)->fP;
fXVar = Ppi.E()/1000.;
return;
};
//********************************************************************
bool MINERvA_CCCOHPI_XSec_1DEpi_nu::isSignal(FitEvent *event) {
//********************************************************************
// Signal: numu + 12C --> mu- + pi+ + 12C'
return SignalDef::isCCCOH(event, 14, 211, EnuMin, EnuMax);
}
-
-double MINERvA_CCCOHPI_XSec_1DEpi_nu::GetLikelihood(){
- double chi2 = StatUtils::GetChi2FromCov(this->fDataHist, this->fMCHist, this->covar, NULL, 1E39, 10);
- return chi2;
-}
diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEpi_nu.h b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEpi_nu.h
index 3233b49..53045e4 100644
--- a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEpi_nu.h
+++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DEpi_nu.h
@@ -1,37 +1,36 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef MINERVA_CCCOHPI_XSEC_1DEPI_NU_H_SEEN
#define MINERVA_CCCOHPI_XSEC_1DEPI_NU_H_SEEN
#include "Measurement1D.h"
class MINERvA_CCCOHPI_XSec_1DEpi_nu : public Measurement1D {
public:
MINERvA_CCCOHPI_XSec_1DEpi_nu(nuiskey samplekey);
virtual ~MINERvA_CCCOHPI_XSec_1DEpi_nu() {};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
- double GetLikelihood();
private:
};
#endif
diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_1Dth_antinu.cxx b/src/MINERvA/MINERvA_CCCOHPI_XSec_1Dth_antinu.cxx
index 1912a63..774c621 100644
--- a/src/MINERvA/MINERvA_CCCOHPI_XSec_1Dth_antinu.cxx
+++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_1Dth_antinu.cxx
@@ -1,92 +1,82 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CCCOHPI_XSec_1Dth_antinu.h"
//********************************************************************
MINERvA_CCCOHPI_XSec_1Dth_antinu::MINERvA_CCCOHPI_XSec_1Dth_antinu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CCCOHPI_XSec_1Dth_antinu sample. \n" \
"Target: CH \n" \
"Flux: MINERvA Reverse Horn Current numu \n" \
"Signal: Any event with 1 mu+, 1pi-, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#theta_{#pi} (deg.)");
fSettings.SetYTitle("d#sigma/d#theta_{#pi} (cm^{2}/deg./C^{12})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(1.5, 20.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numub");
fSettings.SetTitle("MINERvA_CCCOHPI_XSec_1Dth_antinu");
- fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/th_antinu_XSec.csv");
- fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/th_antinu_Covar_stat.csv;" +
- GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/th_antinu_Covar_flux.csv;" +
- GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/th_antinu_Covar_sys.csv");
+ fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Thpi_nubar_data.csv");
+ fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Thpi_nubar_cov.csv");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents)/TotalIntegratedFlux("width")*13;
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetCovarFromMultipleTextFiles(fSettings.GetCovarInput());
- // Apply scalings based on the data release
- ScaleData(1E-39);
-
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void MINERvA_CCCOHPI_XSec_1Dth_antinu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(-211) == 0) return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Ppi = event->GetHMFSParticle(-211)->fP;
fXVar = (180. / M_PI) * FitUtils::th(Pnu, Ppi);
return;
};
//********************************************************************
bool MINERvA_CCCOHPI_XSec_1Dth_antinu::isSignal(FitEvent *event) {
//********************************************************************
// Signal: numub + 12C --> mu+ + pi- + 12C'
return SignalDef::isCCCOH(event, -14, -211, EnuMin, EnuMax);
}
-
-double MINERvA_CCCOHPI_XSec_1Dth_antinu::GetLikelihood(){
- double chi2 = StatUtils::GetChi2FromCov(this->fDataHist, this->fMCHist, this->covar, NULL, 1E39, 10);
- return chi2;
-}
diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_1Dth_antinu.h b/src/MINERvA/MINERvA_CCCOHPI_XSec_1Dth_antinu.h
index 45b7f65..8626bb8 100644
--- a/src/MINERvA/MINERvA_CCCOHPI_XSec_1Dth_antinu.h
+++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_1Dth_antinu.h
@@ -1,37 +1,36 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef MINERVA_CCCOHPI_XSEC_1DTH_ANTINU_H_SEEN
#define MINERVA_CCCOHPI_XSEC_1DTH_ANTINU_H_SEEN
#include "Measurement1D.h"
class MINERvA_CCCOHPI_XSec_1Dth_antinu : public Measurement1D {
public:
MINERvA_CCCOHPI_XSec_1Dth_antinu(nuiskey samplekey);
virtual ~MINERvA_CCCOHPI_XSec_1Dth_antinu() {};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
- double GetLikelihood();
private:
};
#endif
diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_1Dth_nu.cxx b/src/MINERvA/MINERvA_CCCOHPI_XSec_1Dth_nu.cxx
index 95f2706..3219b80 100644
--- a/src/MINERvA/MINERvA_CCCOHPI_XSec_1Dth_nu.cxx
+++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_1Dth_nu.cxx
@@ -1,92 +1,82 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CCCOHPI_XSec_1Dth_nu.h"
//********************************************************************
MINERvA_CCCOHPI_XSec_1Dth_nu::MINERvA_CCCOHPI_XSec_1Dth_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CCCOHPI_XSec_1Dth_nu sample. \n" \
"Target: CH \n" \
"Flux: MINERvA Forward Horn Current numu \n" \
"Signal: Any event with 1 mu-, 1pi+, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#theta_{#pi} (deg.)");
fSettings.SetYTitle("d#sigma/d#theta_{#pi} (cm^{2}/deg./C^{12})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(1.5, 20.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
fSettings.SetTitle("MINERvA_CCCOHPI_XSec_1Dth_nu");
- fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/th_nu_XSec.csv");
- fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/th_nu_Covar_stat.csv;" +
- GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/th_nu_Covar_flux.csv;" +
- GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/th_nu_Covar_sys.csv");
+ fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Thpi_nu_data.csv");
+ fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Thpi_nu_cov.csv");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents)/TotalIntegratedFlux("width")*13;
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetCovarFromMultipleTextFiles(fSettings.GetCovarInput());
- // Apply scalings based on the data release
- ScaleData(1E-39);
-
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void MINERvA_CCCOHPI_XSec_1Dth_nu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(211) == 0) return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Ppi = event->GetHMFSParticle(211)->fP;
fXVar = (180. / M_PI) * FitUtils::th(Pnu, Ppi);
return;
};
//********************************************************************
bool MINERvA_CCCOHPI_XSec_1Dth_nu::isSignal(FitEvent *event) {
//********************************************************************
// Signal: numu + 12C --> mu- + pi+ + 12C'
return SignalDef::isCCCOH(event, 14, 211, EnuMin, EnuMax);
}
-
-double MINERvA_CCCOHPI_XSec_1Dth_nu::GetLikelihood(){
- double chi2 = StatUtils::GetChi2FromCov(this->fDataHist, this->fMCHist, this->covar, NULL, 1E39, 10);
- return chi2;
-}
diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_1Dth_nu.h b/src/MINERvA/MINERvA_CCCOHPI_XSec_1Dth_nu.h
index 1db7fc8..f0e7e4e 100644
--- a/src/MINERvA/MINERvA_CCCOHPI_XSec_1Dth_nu.h
+++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_1Dth_nu.h
@@ -1,37 +1,36 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef MINERVA_CCCOHPI_XSEC_1DTH_NU_H_SEEN
#define MINERVA_CCCOHPI_XSEC_1DTH_NU_H_SEEN
#include "Measurement1D.h"
class MINERvA_CCCOHPI_XSec_1Dth_nu : public Measurement1D {
public:
MINERvA_CCCOHPI_XSec_1Dth_nu(nuiskey samplekey);
virtual ~MINERvA_CCCOHPI_XSec_1Dth_nu() {};
void FillEventVariables(FitEvent *event);
bool isSignal(FitEvent *event);
- double GetLikelihood();
private:
};
#endif
diff --git a/src/MINERvA/MINERvA_CCNpip_XSec_1DEnu_nu.cxx b/src/MINERvA/MINERvA_CCNpip_XSec_1DEnu_nu.cxx
index 6e6703d..78fcb7a 100644
--- a/src/MINERvA/MINERvA_CCNpip_XSec_1DEnu_nu.cxx
+++ b/src/MINERvA/MINERvA_CCNpip_XSec_1DEnu_nu.cxx
@@ -1,90 +1,90 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CCNpip_XSec_1DEnu_nu.h"
//********************************************************************
MINERvA_CCNpip_XSec_1DEnu_nu::MINERvA_CCNpip_XSec_1DEnu_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CCNpip_XSec_1DEnu_nu sample. \n" \
"Target: CH \n" \
"Flux: MINERvA Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("E_{#nu} (GeV)");
fSettings.SetYTitle("d#sigma(E_{#nu}) (cm^{2}/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(1.5, 10.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("MINERvA_CCNpip_XSec_1DEnu_nu");
fSettings.SetDataInput( FitPar::GetDataBase() + "/MINERvA/CCNpip/2016/nu-ccNpi+-xsec-enu.csv" );
fSettings.SetCovarInput( FitPar::GetDataBase() + "/MINERvA/CCNpip/2016/nu-ccNpi+-correlation-enu.csv");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width")*double(1E-38)/double(fNEvents);
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
// MINERvA has the error quoted as a percentage of the cross-section
// Need to make this into an absolute error before we go from correlation matrix -> covariance matrix since it depends on the error in the ith bin
for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) {
fDataHist->SetBinError(i + 1, fDataHist->GetBinContent(i + 1) * (fDataHist->GetBinError(i + 1) / 100.));
}
SetCorrelationFromTextFile(fSettings.GetCovarInput() );
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MINERvA_CCNpip_XSec_1DEnu_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
if (event->NumFSParticle(13) == 0) return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double Enu = Pnu.E()/1000.;
fXVar = Enu;
};
//********************************************************************
// The false refers to that this cross-section uses the full phase space
bool MINERvA_CCNpip_XSec_1DEnu_nu::isSignal(FitEvent *event) {
//********************************************************************
- return SignalDef::isCCNpip_MINERvA(event, EnuMin, EnuMax, false);
+ return SignalDef::isCCNpip_MINERvA(event, EnuMin, EnuMax, false, true);
}
diff --git a/src/MINERvA/MINERvA_CCNpip_XSec_1DTpi_nu.cxx b/src/MINERvA/MINERvA_CCNpip_XSec_1DTpi_nu.cxx
index 23ffdf7..51fa133 100644
--- a/src/MINERvA/MINERvA_CCNpip_XSec_1DTpi_nu.cxx
+++ b/src/MINERvA/MINERvA_CCNpip_XSec_1DTpi_nu.cxx
@@ -1,267 +1,267 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CCNpip_XSec_1DTpi_nu.h"
//********************************************************************
MINERvA_CCNpip_XSec_1DTpi_nu::MINERvA_CCNpip_XSec_1DTpi_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CCNpip_XSec_1DTpi_nu sample. \n" \
"Target: CH \n" \
"Flux: MINERvA Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("T_{#pi} (MeV)");
fSettings.SetYTitle("(1/T#Phi) dN_{#pi}/dT_{#pi} (cm^{2}/MeV/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(1.5, 10.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
fFullPhaseSpace = !fSettings.Found("name", "_20deg");
fFluxCorrection = fSettings.Found("name", "fluxcorr");
fUpdatedData = !fSettings.Found("name", "2015");
fSettings.SetTitle("MINERvA_CCNpip_XSec_1DTpi_nu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents) / TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
// Full Phase Space
if (fFullPhaseSpace) {
// 2016 release
if (fUpdatedData) {
SetDataFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2016/nu-ccNpi+-xsec-pion-kinetic-energy.csv");
// MINERvA has the error quoted as a percentage of the cross-section
// Need to make this into an absolute error before we go from correlation
// matrix -> covariance matrix since it depends on the error in the ith
// bin
for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) {
fDataHist->SetBinError(i + 1, fDataHist->GetBinContent(i + 1) * (fDataHist->GetBinError(i + 1) / 100.));
}
// This is a correlation matrix, not covariance matrix, so needs to be
// converted
SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2016/nu-ccNpi+-correlation-pion-kinetic-energy.csv");
// 2015 release
} else {
// If we're doing shape only
if (fIsShape) {
SetDataFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_Tpi_shape.txt");
SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_Tpi_shape_cov.txt");
// If we're doing full cross-section
} else {
SetDataFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_Tpi.txt");
SetCorrelationFromTextFile( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_Tpi_cov.txt");
}
}
// Restricted Phase Space
} else {
// Only 2015 data released restricted muon phase space cross-section
// unfortunately
if (fUpdatedData) {
ERR(FTL) << fName << " has no updated 2016 data for restricted phase space! Using 2015 data." << std::endl;
throw;
}
// If we're using the shape only data
if (fIsShape) {
SetDataFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_Tpi_20deg_shape.txt");
SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_Tpi_20deg_shape_cov.txt");
// Or total cross-section
} else {
SetDataFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_Tpi_20deg.txt");
SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_Tpi_20deg_cov.txt");
}
}
// Scale the MINERvA data to account for the flux difference
// Adjust MINERvA data to flux correction; roughly a 11% normalisation increase in data
// Please change when MINERvA releases new data!
if (fFluxCorrection) {
for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) {
fDataHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) * 1.11);
}
}
// Make some auxillary helper plots
onePions = (TH1D*)(fDataHist->Clone());
onePions->SetNameTitle((fName + "_1pions").c_str(), (fName + "_1pions" + fPlotTitles).c_str());
SetAutoProcessTH1(onePions, kCMD_Reset, kCMD_Scale, kCMD_Norm);
twoPions = (TH1D*)(fDataHist->Clone());
twoPions->SetNameTitle((fName + "_2pions").c_str(), (fName + "_2pions;" + fPlotTitles).c_str());
SetAutoProcessTH1(twoPions, kCMD_Reset, kCMD_Scale, kCMD_Norm);
threePions = (TH1D*)(fDataHist->Clone());
threePions->SetNameTitle((fName + "_3pions").c_str(), (fName + "_3pions" + fPlotTitles).c_str());
SetAutoProcessTH1(threePions, kCMD_Reset, kCMD_Scale, kCMD_Norm);
morePions = (TH1D*)(fDataHist->Clone());
morePions->SetNameTitle((fName + "_4pions").c_str(), (fName + "_4pions" + fPlotTitles).c_str());
SetAutoProcessTH1(morePions, kCMD_Reset, kCMD_Scale, kCMD_Norm);
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
// Here we have to fill for every pion we find in the event
void MINERvA_CCNpip_XSec_1DTpi_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
if (event->NumFSParticle(211) == 0 && event->NumFSParticle(-211) == 0) return;
if (event->NumFSParticle(13) == 0) return;
// Need to make this use event boxes
// Clear out the vectors
GetPionBox()->Reset();
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
// Loop over the particle stack
for (unsigned int j = 2; j < event->Npart(); ++j) {
// Only include alive particles
if (event->GetParticleState(j) != kFinalState) continue;
int PID = (event->PartInfo(j))->fPID;
// Pick up the charged pions in the event
if (abs(PID) == 211) {
double ppi = FitUtils::T(event->PartInfo(j)->fP) * 1000.;
GetPionBox()->fTpiVect.push_back(ppi);
}
}
fXVar = 0;
return;
};
//********************************************************************
// The last bool refers to if we're using restricted phase space or not
bool MINERvA_CCNpip_XSec_1DTpi_nu::isSignal(FitEvent *event) {
//********************************************************************
// Last false refers to that this is NOT the restricted MINERvA phase space,
// in which only forward-going muons are accepted
- return SignalDef::isCCNpip_MINERvA(event, EnuMin, EnuMax, !fFullPhaseSpace);
+ return SignalDef::isCCNpip_MINERvA(event, EnuMin, EnuMax, !fFullPhaseSpace, !fUpdatedData);
}
//********************************************************************
// Need to override FillHistograms() here because we fill the histogram N_pion
// times
void MINERvA_CCNpip_XSec_1DTpi_nu::FillHistograms() {
//********************************************************************
if (Signal) {
unsigned int nPions = GetPionBox()->fTpiVect.size();
// Need to loop over all the pions in the sample
for (size_t k = 0; k < nPions; ++k) {
double tpi = GetPionBox()->fTpiVect[k];
this->fMCHist->Fill(tpi, Weight);
this->fMCFine->Fill(tpi, Weight);
this->fMCStat->Fill(tpi, 1.0);
if (nPions == 1) {
onePions->Fill(tpi, Weight);
} else if (nPions == 2) {
twoPions->Fill(tpi, Weight);
} else if (nPions == 3) {
threePions->Fill(tpi, Weight);
} else if (nPions > 3) {
morePions->Fill(tpi, Weight);
}
if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, tpi, Weight);
}
}
}
//********************************************************************
void MINERvA_CCNpip_XSec_1DTpi_nu::ScaleEvents() {
//********************************************************************
Measurement1D::ScaleEvents();
onePions->Scale(this->fScaleFactor, "width");
twoPions->Scale(this->fScaleFactor, "width");
threePions->Scale(this->fScaleFactor, "width");
morePions->Scale(this->fScaleFactor, "width");
return;
}
//********************************************************************
void MINERvA_CCNpip_XSec_1DTpi_nu::Write(std::string drawOpts) {
//********************************************************************
Measurement1D::Write(drawOpts);
// Make an auto processed pion stack
// Draw the npions stack
onePions->SetTitle("1#pi");
onePions->SetLineColor(kBlack);
// onePions->SetFillStyle(0);
onePions->SetFillColor(onePions->GetLineColor());
twoPions->SetTitle("2#pi");
twoPions->SetLineColor(kRed);
// twoPions->SetFillStyle(0);
twoPions->SetFillColor(twoPions->GetLineColor());
threePions->SetTitle("3#pi");
threePions->SetLineColor(kGreen);
// threePions->SetFillStyle(0);
threePions->SetFillColor(threePions->GetLineColor());
morePions->SetTitle(">3#pi");
morePions->SetLineColor(kBlue);
// morePions->SetFillStyle(0);
morePions->SetFillColor(morePions->GetLineColor());
THStack pionStack =
THStack((fName + "_pionStack").c_str(), (fName + "_pionStack").c_str());
pionStack.Add(onePions);
pionStack.Add(twoPions);
pionStack.Add(threePions);
pionStack.Add(morePions);
pionStack.Write();
return;
}
diff --git a/src/MINERvA/MINERvA_CCNpip_XSec_1Dth_nu.cxx b/src/MINERvA/MINERvA_CCNpip_XSec_1Dth_nu.cxx
index f753eba..f12a5fe 100644
--- a/src/MINERvA/MINERvA_CCNpip_XSec_1Dth_nu.cxx
+++ b/src/MINERvA/MINERvA_CCNpip_XSec_1Dth_nu.cxx
@@ -1,253 +1,253 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CCNpip_XSec_1Dth_nu.h"
//********************************************************************
MINERvA_CCNpip_XSec_1Dth_nu::MINERvA_CCNpip_XSec_1Dth_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CCNpip_XSec_1Dth_nu sample. \n" \
"Target: CH \n" \
"Flux: MINERvA Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#theta_{#pi} (degrees)");
fSettings.SetYTitle("(1/T#Phi) dN_{#pi}/d#theta_{#pi} (cm^{2}/degrees/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(1.5, 10.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
fFullPhaseSpace = !fSettings.Found("name", "_20deg");
fFluxCorrection = fSettings.Found("name", "fluxcorr");
fUpdatedData = !fSettings.Found("name", "2015");
fSettings.SetTitle("MINERvA_CCNpip_XSec_1Dth_nu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents) / TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
// Full Phase Space
if (fFullPhaseSpace) {
// 2016 release data
if (fUpdatedData) {
SetDataFromTextFile( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2016/nu-ccNpi+-xsec-pion-angle.csv");
// MINERvA has the error quoted as a percentage of the cross-section
// Need to make this into an absolute error before we go from correlation matrix -> covariance matrix since it depends on the error in the ith bin
for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) {
fDataHist->SetBinError(i + 1, fDataHist->GetBinContent(i + 1) * (fDataHist->GetBinError(i + 1) / 100.));
}
// This is a correlation matrix! but it's all fixed in SetCovarMatrixFromText
SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2016/nu-ccNpi+-correlation-pion-angle.csv");
// 2015 release data
} else {
// 2015 release allows for shape comparisons
if (fIsShape) {
SetDataFromTextFile( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_th_shape.txt");
SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_th_shape_cov.txt");
} else {
SetDataFromTextFile( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_th.txt");
SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_th_cov.txt");
}
}
// Restricted Phase Space Data
} else {
// 2016 release data unfortunately not released in 20degree forward-going, revert to 2015 data
if (fUpdatedData) {
ERR(FTL) << fName << " has no updated 2016 data for restricted phase space! Using 2015 data." << std::endl;
throw;
}
// Only 2015 20deg data
if (fIsShape) {
SetDataFromTextFile( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_th_20deg_shape.txt");
SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_th_20deg_shape_cov.txt");
} else {
SetDataFromTextFile( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_th_20deg.txt");
SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_th_20deg_cov.txt");
}
}
// Scale the MINERvA data to account for the flux difference
// Adjust MINERvA data to flux correction; roughly a 11% normalisation increase in data
// Please change when MINERvA releases new data!
if (fFluxCorrection) {
for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) {
fDataHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) * 1.11);
}
}
// Make some auxillary helper plots
onePions = (TH1D*)(fDataHist->Clone());
onePions->SetNameTitle((fName + "_1pions").c_str(), (fName + "_1pions" + fPlotTitles).c_str());
SetAutoProcessTH1(onePions, kCMD_Reset, kCMD_Scale, kCMD_Norm);
twoPions = (TH1D*)(fDataHist->Clone());
twoPions->SetNameTitle((fName + "_2pions").c_str(), (fName + "_2pions;" + fPlotTitles).c_str());
SetAutoProcessTH1(twoPions, kCMD_Reset, kCMD_Scale, kCMD_Norm);
threePions = (TH1D*)(fDataHist->Clone());
threePions->SetNameTitle((fName + "_3pions").c_str(), (fName + "_3pions" + fPlotTitles).c_str());
SetAutoProcessTH1(threePions, kCMD_Reset, kCMD_Scale, kCMD_Norm);
morePions = (TH1D*)(fDataHist->Clone());
morePions->SetNameTitle((fName + "_4pions").c_str(), (fName + "_4pions" + fPlotTitles).c_str());
SetAutoProcessTH1(morePions, kCMD_Reset, kCMD_Scale, kCMD_Norm);
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
// ********************************************
// Fill the event variables
// Here we want to fill the angle for every pion we can find in the event
void MINERvA_CCNpip_XSec_1Dth_nu::FillEventVariables(FitEvent *event) {
// ********************************************
if (event->NumFSParticle(211) == 0 && event->NumFSParticle(-211) == 0) return;
if (event->NumFSParticle(13) == 0) return;
GetBox()->Reset();
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
TLorentzVector Ppip;
// Loop over the particle stack
for (unsigned int j = 2; j < event->Npart(); ++j) {
// Only include alive particles
if ((event->PartInfo(j))->fIsAlive <= 0) continue;
if ((event->PartInfo(j))->fNEUTStatusCode != 0) continue;
int PID = (event->PartInfo(j))->fPID;
// Select highest momentum (energy) charged pion
if (abs(PID) == 211) {
Ppip = (event->PartInfo(j))->fP;
double th = (180. / M_PI) * FitUtils::th(Pnu, Ppip);
GetPionBox()->fthpiVect.push_back(th);
}
}
fXVar = 0;
return;
};
//********************************************************************
// The signal definition for MINERvA CCNpi+
// Last bool refers to if we're selecting on the full phase space or not
bool MINERvA_CCNpip_XSec_1Dth_nu::isSignal(FitEvent *event) {
//********************************************************************
- return SignalDef::isCCNpip_MINERvA(event, EnuMin, EnuMax, !fFullPhaseSpace);
+ return SignalDef::isCCNpip_MINERvA(event, EnuMin, EnuMax, !fFullPhaseSpace, !fUpdatedData);
}
//********************************************************************
// Need to override FillHistograms() here because we fill the histogram N_pion times
void MINERvA_CCNpip_XSec_1Dth_nu::FillHistograms() {
//********************************************************************
if (Signal) {
unsigned int nPions = GetPionBox()->fthpiVect.size();
// Need to loop over all the pions in the event
for (size_t k = 0; k < nPions; ++k) {
double th = GetPionBox()->fthpiVect[k];
this->fMCHist->Fill(th, Weight);
this->fMCFine->Fill(th, Weight);
this->fMCStat->Fill(th, 1.0);
if (nPions == 1) {
onePions->Fill(th, Weight);
} else if (nPions == 2) {
twoPions->Fill(th, Weight);
} else if (nPions == 3) {
threePions->Fill(th, Weight);
} else if (nPions > 3) {
morePions->Fill(th, Weight);
}
if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, th, Weight);
// PlotUtils::FillNeutModeArray(fMCHist_PDG, Mode, th, Weight);
}
}
return;
}
//********************************************************************
void MINERvA_CCNpip_XSec_1Dth_nu::Write(std::string drawOpts) {
//********************************************************************
Measurement1D::Write(drawOpts);
// Draw the npions stack
onePions->SetTitle("1#pi");
onePions->SetLineColor(kBlack);
//onePions->SetFillStyle(0);
onePions->SetFillColor(onePions->GetLineColor());
twoPions->SetTitle("2#pi");
twoPions->SetLineColor(kRed);
//twoPions->SetFillStyle(0);
twoPions->SetFillColor(twoPions->GetLineColor());
threePions->SetTitle("3#pi");
threePions->SetLineColor(kGreen);
//threePions->SetFillStyle(0);
threePions->SetFillColor(threePions->GetLineColor());
morePions->SetTitle(">3#pi");
morePions->SetLineColor(kBlue);
//morePions->SetFillStyle(0);
morePions->SetFillColor(morePions->GetLineColor());
THStack pionStack = THStack((fName + "_pionStack").c_str(), (fName + "_pionStack").c_str());
pionStack.Add(onePions);
pionStack.Add(twoPions);
pionStack.Add(threePions);
pionStack.Add(morePions);
pionStack.Write();
return;
}
diff --git a/src/MINERvA/MINERvA_SignalDef.cxx b/src/MINERvA/MINERvA_SignalDef.cxx
index fa7d34c..8c2be8c 100644
--- a/src/MINERvA/MINERvA_SignalDef.cxx
+++ b/src/MINERvA/MINERvA_SignalDef.cxx
@@ -1,271 +1,305 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "SignalDef.h"
#include "FitUtils.h"
#include "MINERvA_SignalDef.h"
namespace SignalDef {
// *********************************
// MINERvA CC1pi+/- signal definition (2015 release)
// Note: There is a 2016 release which is different to this (listed below), but
// it is CCNpi+ and has a different W cut
// Note2: The W cut is implemented in the class implementation in MINERvA/
// rather than here so we can draw events that don't pass the W cut (W cut is
// 1.4 GeV)
// Could possibly be changed for slight speed increase since less events
// would be used
//
// MINERvA signal is slightly different to MiniBooNE
//
// Exactly one negative muon
// Exactly one charged pion (both + and -); however, there is a Michel e-
// requirement but UNCLEAR IF UNFOLDED OR NOT (so don't know if should be
// applied)
// Exactly 1 charged pion exits (so + and - charge), however, has Michel
// electron requirement, so look for + only here?
// No restriction on neutral pions or other mesons
// MINERvA has unfolded and not unfolded muon phase space for 2015
//
// Possible issues with the data:
// 1) pi- is allowed in signal even when Michel cut included; most pi- is efficiency corrected in GENIE
// 2) There is a T_pi < 350 MeV cut coming from requiring a stopping pion; this is efficiency corrected in GENIE
// 3) There is a 1.5 < Enu < 10.0 cut in signal definition
// 4) There is an angular muon cut which is sometimes efficiency corrected (why we have bool isRestricted below)
//
// Nice things:
// Much data given: with and without muon angle cuts and with and without shape
// only data + covariance
//
bool isCC1pip_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool isRestricted) {
// *********************************
// Signal is both pi+ and pi-
// WARNING: PI- CONTAMINATION IS FULLY GENIE BECAUSE THE MICHEL TAG
// First, make sure it's CCINC
if (!isCCINC(event, 14, EnuMin, EnuMax)) return false;
// Allow pi+/pi-
int piPDG[] = {211, -211};
int nLeptons = event->NumFSLeptons();
int nPion = event->NumFSParticle(piPDG);
// Check that the desired pion exists and is the only meson
if (nPion != 1) return false;
// Check that there is only one final state lepton
if (nLeptons != 1) return false;
// MINERvA released another CC1pi+ xsec without muon unfolding!
// here the muon angle is < 20 degrees (seen in MINOS)
TLorentzVector pnu = event->GetHMISParticle(14)->fP;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
if (isRestricted) {
double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
if (th_nu_mu >= 20) return false;
}
// Extract Hadronic Mass
double hadMass = FitUtils::Wrec(pnu, pmu);
// Actual cut is True GENIE Ws! Arg.! Use gNtpcConv definition.
#ifdef __GENIE_ENABLED__
if (event->fType == kGENIE){
EventRecord * gevent = static_cast(event->genie_event->event);
const Interaction * interaction = gevent->Summary();
const Kinematics & kine = interaction->Kine();
double Ws = kine.W (true);
// std::cout << "Ws versus WRec = " << Ws << " vs " << hadMass << " " << kine.W(false) << std::endl;
hadMass = Ws * 1000.0;
}
#endif
if (hadMass > 1400.0) return false;
return true;
};
+
+ // Updated MINERvA 2017 Signal using Wexp and no restriction on angle
+bool isCC1pip_MINERvA_2017(FitEvent *event, double EnuMin, double EnuMax){
+
+ // Signal is both pi+ and pi-
+ // WARNING: PI- CONTAMINATION IS FULLY GENIE BECAUSE THE MICHEL TAG
+ // First, make sure it's CCINC
+ if (!isCCINC(event, 14, EnuMin, EnuMax)) return false;
+
+ // Allow pi+/pi-
+ int piPDG[] = {211, -211};
+ int nLeptons = event->NumFSLeptons();
+ int nPion = event->NumFSParticle(piPDG);
+
+ // Check that the desired pion exists and is the only meson
+ if (nPion != 1) return false;
+
+ // Check that there is only one final state lepton
+ if (nLeptons != 1) return false;
+
+ // Get Muon and Lepton Kinematics
+ TLorentzVector pnu = event->GetHMISParticle(14)->fP;
+ TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
+
+ // Extract Hadronic Mass
+ double hadMass = FitUtils::Wrec(pnu, pmu);
+ if (hadMass > 1400.0) return false;
+
+ return true;
+};
+
// *********************************
// MINERvA CCNpi+/- signal definition from 2016 publication
// Different to CC1pi+/- listed above; additional has W < 1.8 GeV
//
// For notes on strangeness of signal definition, see CC1pip_MINERvA
//
// One negative muon
// At least one charged pion
// 1.5 < Enu < 10
// No restrictions on pi0 or other mesons or baryons
// W_reconstructed (ignoring initial state motion) cut at 1.8 GeV
//
// Also writes number of pions (nPions) if studies on this want to be done...
bool isCCNpip_MINERvA(FitEvent *event, double EnuMin,
- double EnuMax, bool isRestricted) {
+ double EnuMax, bool isRestricted, bool isWtrue) {
// *********************************
// First, make sure it's CCINC
if (!isCCINC(event, 14, EnuMin, EnuMax)) return false;
int nLeptons = event->NumFSLeptons();
// Write the number of pions to the measurement class...
// Maybe better to just do that inside the class?
int nPions = event->NumFSParticle(PhysConst::pdg_charged_pions);
// Check that there is a pion!
if (nPions == 0) return false;
// Check that there is only one final state lepton
if (nLeptons != 1) return false;
// Need the muon and the neutrino to check angles and W
TLorentzVector pnu = event->GetNeutrinoIn()->fP;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
// MINERvA released some data with restricted muon angle
// Here the muon angle is < 20 degrees (seen in MINOS)
if (isRestricted) {
double th_nu_mu = FitUtils::th(pmu, pnu) * 180. / M_PI;
if (th_nu_mu >= 20.) return false;
}
// Lastly check the W cut (always at 1.8 GeV)
double Wrec = FitUtils::Wrec(pnu, pmu) + 0.;
// Actual cut is True GENIE Ws! Arg.! Use gNtpcConv definition.
+ if (isWtrue){
#ifdef __GENIE_ENABLED__
- if (event->fType == kGENIE){
- GHepRecord* ghep = static_cast(event->genie_event->event);
- const Interaction * interaction = ghep->Summary();
- const Kinematics & kine = interaction->Kine();
- double Ws = kine.W (true);
- Wrec = Ws * 1000.0; // Say Wrec is Ws
- }
+ if (event->fType == kGENIE){
+ GHepRecord* ghep = static_cast(event->genie_event->event);
+ const Interaction * interaction = ghep->Summary();
+ const Kinematics & kine = interaction->Kine();
+ double Ws = kine.W (true);
+ Wrec = Ws * 1000.0; // Say Wrec is Ws
+ }
#endif
+ }
+
if (Wrec > 1800. || Wrec < 0.0) return false;
return true;
};
//********************************************************************
bool isCCQEnumu_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool fullphasespace) {
//********************************************************************
if (!isCCQELike(event, 14, EnuMin, EnuMax)) return false;
TLorentzVector pnu = event->GetHMISParticle(14)->fP;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
double ThetaMu = pnu.Vect().Angle(pmu.Vect());
double Enu_rec = FitUtils::EnuQErec(pmu, cos(ThetaMu), 34., true);
// If Restricted phase space
if (!fullphasespace && ThetaMu > 0.34906585) return false;
// restrict energy range
if (Enu_rec < EnuMin || Enu_rec > EnuMax) return false;
return true;
};
//********************************************************************
bool isCCQEnumubar_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool fullphasespace) {
//********************************************************************
if (!isCCQELike(event, -14, EnuMin, EnuMax)) return false;
TLorentzVector pnu = event->GetHMISParticle(-14)->fP;
TLorentzVector pmu = event->GetHMFSParticle(-13)->fP;
double ThetaMu = pnu.Vect().Angle(pmu.Vect());
double Enu_rec = FitUtils::EnuQErec(pmu, cos(ThetaMu), 30., true);
// If Restricted phase space
if (!fullphasespace && ThetaMu > 0.34906585) return false;
// restrict energy range
if (Enu_rec < EnuMin || Enu_rec > EnuMax) return false;
return true;
}
//********************************************************************
bool isCCincLowRecoil_MINERvA(FitEvent *event, double EnuMin, double EnuMax) {
//********************************************************************
if (!isCCINC(event, 14, EnuMin, EnuMax)) return false;
// Need at least one muon
if (event->NumFSParticle(13) < 1) return false;
TLorentzVector pmu = event->GetHMFSParticle(13)->fP;
TLorentzVector pnu = event->GetHMISParticle(14)->fP;
// Cut on muon angle greated than 20deg
if (cos(pnu.Vect().Angle(pmu.Vect())) < 0.93969262078) return false;
// Cut on muon energy < 1.5 GeV
if (pmu.E()/1000.0 < 1.5) return false;
return true;
}
bool isCC0pi1p_MINERvA(FitEvent *event, double enumin, double enumax) {
// Require numu CC0pi event with a proton above threshold
bool signal = (isCC0pi(event, 14, enumin, enumax) &&
HasProtonKEAboveThreshold(event, 110.0));
return signal;
}
// 2015 analysis just asks for 1pi0 and no pi+/pi-
bool isCC1pi0_MINERvA_2015(FitEvent *event, double EnuMin, double EnuMax) {
bool CC1pi0_anu = SignalDef::isCC1pi(event, -14, 111, EnuMin, EnuMax);
return CC1pi0_anu;
}
// 2016 analysis just asks for 1pi0 and no other charged tracks
bool isCC1pi0_MINERvA_2016(FitEvent *event, double EnuMin, double EnuMax) {
bool CC1pi0_anu = SignalDef::isCC1pi(event, -14, 111, EnuMin, EnuMax);
/*
// Additionally look for charged proton track
bool HasProton = event->HasFSParticle(2212);
if (CC1pi0_anu) {
if (!HasProton) {
return true;
} else {
return false;
}
} else {
return false;
}
*/
return CC1pi0_anu;
}
}
diff --git a/src/MINERvA/MINERvA_SignalDef.h b/src/MINERvA/MINERvA_SignalDef.h
index 34e9ac2..2b1f03e 100644
--- a/src/MINERvA/MINERvA_SignalDef.h
+++ b/src/MINERvA/MINERvA_SignalDef.h
@@ -1,92 +1,94 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef MINERVA_SIGNALDEF_H_SEEN
#define MINERVA_SIGNALDEF_H_SEEN
#include "FitEvent.h"
namespace SignalDef {
// *********************************
/// MINERvA CC1pi+/- signal definition (2015 release)
/// Note: There is a 2016 release which is different to this (listed below),
/// but
/// it is CCNpi+ and has a different W cut
/// Note2: The W cut is implemented in the class implementation in MINERvA/
/// rather than here so we can draw events that don't pass the W cut (W cut is
/// 1.4 GeV)
/// Could possibly be changed for slight speed increase since less events
/// would be used
///
/// MINERvA signal is slightly different to MiniBooNE
///
/// Exactly one negative muon
/// Exactly one charged pion (both + and -); however, there is a Michel e-
/// requirement but UNCLEAR IF UNFOLDED OR NOT (so don't know if should be
/// applied)
/// Exactly 1 charged pion exits (so + and - charge), however, has Michel
/// electron requirement, so look for + only here?
/// No restriction on neutral pions or other mesons
/// MINERvA has unfolded and not unfolded muon phase space for 2015
///
/// Possible problems:
/// 1) Should there be a pi+ only cut implemented due to Michel requirement, or
/// is pi- events filled from MC?
/// 2) There is a T_pi < 350 MeV cut coming from requiring a stopping pion so
/// the
/// Michel e is seen, this is also unclear if it's unfolded so any pion is OK
///
/// Nice things:
/// Much data given: with and without muon angle cuts and with and without shape
/// only data + covariance
bool isCC1pip_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool isRestricted = false);
+bool isCC1pip_MINERvA_2017(FitEvent *event, double EnuMin, double EnuMax);
+
// *********************************
/// MINERvA CCNpi+/- signal definition from 2016 publication
/// Different to CC1pi+/- listed above; additional has W < 1.8 GeV
///
/// Still asks for a Michel e and still unclear if this is unfolded or not
/// Says stuff like "requirement that a Michel e isolates a subsample that is
/// more nearly a pi+ prodution", yet the signal definition is both pi+ and pi-?
///
/// One negative muon
/// At least one charged pion
/// 1.5 < Enu < 10
/// No restrictions on pi0 or other mesons or baryons
///
/// Also writes number of pions (nPions) if studies on this want to be done...
bool isCCNpip_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
- bool isRestricted = false);
+ bool isRestricted = false, bool isWtrue=false);
bool isCCQEnumu_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool fullphasespace = true);
bool isCCQEnumubar_MINERvA(FitEvent *event, double EnuMin, double EnuMax,
bool fullphasespace = true);
bool isCCincLowRecoil_MINERvA(FitEvent *event, double EnuMin, double EnuMax);
bool isCC0pi1p_MINERvA(FitEvent *event, double enumin, double enumax);
bool isCC1pi0_MINERvA_2015(FitEvent *event, double EnuMin, double EnuMax);
bool isCC1pi0_MINERvA_2016(FitEvent *event, double EnuMin, double EnuMax);
}
#endif