diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt
index 6e87101..8a18777 100644
--- a/app/CMakeLists.txt
+++ b/app/CMakeLists.txt
@@ -1,195 +1,100 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
set(TARGETS_TO_BUILD)
if(USE_MINIMIZER)
add_executable(nuismin nuismin.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuismin)
- target_link_libraries(nuismin ${MODULETargets})
- target_link_libraries(nuismin ${CMAKE_DEPENDLIB_FLAGS})
- # target_link_libraries(nuismin ${ROOT_LIBS})
- if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
- set_target_properties(nuismin PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
- endif()
add_executable(nuissplines nuissplines.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissplines)
- target_link_libraries(nuissplines ${MODULETargets})
- target_link_libraries(nuissplines ${CMAKE_DEPENDLIB_FLAGS})
- # target_link_libraries(nuissplines ${ROOT_LIBS})
- if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
- set_target_properties(nuissplines PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
- endif()
endif()
include_directories(${RWENGINE_INCLUDE_DIRECTORIES})
include_directories(${CMAKE_SOURCE_DIR}/src/Routines)
include_directories(${CMAKE_SOURCE_DIR}/src/InputHandler)
include_directories(${CMAKE_SOURCE_DIR}/src/Genie)
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
include_directories(${CMAKE_SOURCE_DIR}/src/Statistical)
include_directories(${CMAKE_SOURCE_DIR}/src/Utils)
include_directories(${CMAKE_SOURCE_DIR}/src/Config)
include_directories(${CMAKE_SOURCE_DIR}/src/Logger)
include_directories(${CMAKE_SOURCE_DIR}/src/Splines)
include_directories(${CMAKE_SOURCE_DIR}/src/Reweight)
include_directories(${CMAKE_SOURCE_DIR}/src/FCN)
include_directories(${CMAKE_SOURCE_DIR}/src/MCStudies)
include_directories(${CMAKE_SOURCE_DIR}/src/Smearceptance)
include_directories(${EXP_INCLUDE_DIRECTORIES})
if (USE_NuWro AND NOT NUWRO_BUILT_FROM_FILE)
add_executable(nuwro_nuisance nuwro_NUISANCE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuwro_nuisance)
- target_link_libraries(nuwro_nuisance ${MODULETargets})
- target_link_libraries(nuwro_nuisance ${CMAKE_DEPENDLIB_FLAGS})
- # target_link_libraries(nuwro_nuisance ${ROOT_LIBS})
- include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
- if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
- set_target_properties(nuwro_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
- endif()
endif()
-# if (USE_NEUT)
-# add_executable(neut_nuisance neut_NUISANCE.cxx)
-# set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};neut_nuisance)
-# target_link_libraries(neut_nuisance ${MODULETargets})
-# target_link_libraries(neut_nuisance ${CMAKE_DEPENDLIB_FLAGS})
-# target_link_libraries(neut_nuisance ${ROOT_LIBS})
-# include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
-# if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
-# set_target_properties(neut_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
-# endif()
-# endif()
-
if (USE_GiBUU)
add_executable(DumpGiBUUEvents DumpGiBUUEvents.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};DumpGiBUUEvents)
- target_link_libraries(DumpGiBUUEvents ${MODULETargets})
- target_link_libraries(DumpGiBUUEvents ${CMAKE_DEPENDLIB_FLAGS})
- # target_link_libraries(DumpGiBUUEvents ${ROOT_LIBS})
- include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
- if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
- set_target_properties(DumpGiBUUEvents PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
- endif()
endif()
add_executable(nuiscomp nuiscomp.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuiscomp)
-target_link_libraries(nuiscomp ${MODULETargets})
-target_link_libraries(nuiscomp ${CMAKE_DEPENDLIB_FLAGS})
-# target_link_libraries(nuiscomp ${ROOT_LIBS})
-if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
- set_target_properties(nuiscomp PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
-endif()
add_executable(nuisflat nuisflat.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuisflat)
-target_link_libraries(nuisflat ${MODULETargets})
-target_link_libraries(nuisflat ${CMAKE_DEPENDLIB_FLAGS})
-# target_link_libraries(nuisflat ${ROOT_LIBS})
-if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
- set_target_properties(nuisflat PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
-endif()
add_executable(nuissmear nuissmear.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissmear)
-target_link_libraries(nuissmear ${MODULETargets})
-target_link_libraries(nuissmear ${CMAKE_DEPENDLIB_FLAGS})
-# target_link_libraries(nuissmear ${ROOT_LIBS})
-if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
- set_target_properties(nuissmear PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
-endif()
add_executable(nuissyst nuissyst.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissyst)
-target_link_libraries(nuissyst ${MODULETargets})
-target_link_libraries(nuissyst ${CMAKE_DEPENDLIB_FLAGS})
-# target_link_libraries(nuissyst ${ROOT_LIBS})
-if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
- set_target_properties(nuissyst PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
-endif()
add_executable(nuisbayes nuisbayes.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuisbayes)
-target_link_libraries(nuisbayes ${MODULETargets})
-target_link_libraries(nuisbayes ${CMAKE_DEPENDLIB_FLAGS})
-# target_link_libraries(nuisbayes ${ROOT_LIBS})
-if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
- set_target_properties(nuisbayes PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
-endif()
if(USE_GENIE)
add_executable(PrepareGENIE PrepareGENIE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareGENIE)
- target_link_libraries(PrepareGENIE ${MODULETargets})
- target_link_libraries(PrepareGENIE ${CMAKE_DEPENDLIB_FLAGS})
- # target_link_libraries(PrepareGENIE ${ROOT_LIBS})
- if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
- set_target_properties(PrepareGENIE PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
- endif()
endif()
if(USE_NEUT)
add_executable(PrepareNEUT PrepareNEUT.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareNEUT)
- target_link_libraries(PrepareNEUT ${MODULETargets})
- target_link_libraries(PrepareNEUT ${CMAKE_DEPENDLIB_FLAGS})
- # target_link_libraries(PrepareNEUT ${ROOT_LIBS})
- if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
- set_target_properties(PrepareNEUT PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
- endif()
endif()
# PREPARE NUWRO
# Commented out for the time being until it is finished..
if(USE_NuWro)
add_executable(PrepareNuwro PrepareNuwroEvents.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareNuwro)
- target_link_libraries(PrepareNuwro ${MODULETargets})
- target_link_libraries(PrepareNuwro ${CMAKE_DEPENDLIB_FLAGS})
- # target_link_libraries(PrepareNuwro ${ROOT_LIBS})
- if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
- set_target_properties(PrepareNuwro PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
- endif()
endif()
add_executable(nuisbac nuisbac.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuisbac)
-target_link_libraries(nuisbac ${MODULETargets})
-target_link_libraries(nuisbac ${CMAKE_DEPENDLIB_FLAGS})
-if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
- set_target_properties(nuisbac PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
-endif()
-install(TARGETS ${TARGETS_TO_BUILD} DESTINATION bin)
-
-#add_executable(DumpROOTClassesFromVector DumpROOTClassesFromVector.cxx)
-# #Strip out -lNuWro_event1
-# string(REPLACE "-lNuWro_event1" "" NWEVSTRIPPED_CDF ${CMAKE_DEPENDLIB_FLAGS})
-# cmessage(DEBUG "Attempted to strip out nuwro library: \"${CMAKE_DEPENDLIB_FLAGS}\" -> \"${NWEVSTRIPPED_CDF}\"")
-# add_executable(PrepareNEUT PrepareNEUT.cxx)
-# target_link_libraries(DumpROOTClassesFromVector ${MODULETargets})
-# target_link_libraries(DumpROOTClassesFromVector ${NWEVSTRIPPED_CDF})
-# if(NOT CMAKE_LINK_FLAGS STREQUAL "")
-# set_target_properties(DumpROOTClassesFromVector PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
-# endif()
-#install(TARGETS DumpROOTClassesFromVector DESTINATION bin)
+foreach(targ ${TARGETS_TO_BUILD})
+ if(NOT "${CMAKE_LINK_FLAGS} ${NUIS_EXE_FLAGS}" STREQUAL " ")
+ set_target_properties(${targ} PROPERTIES LINK_FLAGS "${CMAKE_LINK_FLAGS} ${NUIS_EXE_FLAGS}")
+ endif()
+ target_link_libraries(${targ} ${MODULETargets})
+ target_link_libraries(${targ} ${CMAKE_DEPENDLIB_FLAGS})
+endforeach()
+
+install(TARGETS ${TARGETS_TO_BUILD} DESTINATION bin)
\ No newline at end of file
diff --git a/cmake/NEUTSetup.cmake b/cmake/NEUTSetup.cmake
index 0bb4da7..4cd20ba 100644
--- a/cmake/NEUTSetup.cmake
+++ b/cmake/NEUTSetup.cmake
@@ -1,227 +1,232 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
include(cmake/parseConfigApp.cmake)
find_program(NEUTCONFIG NAMES neut-config)
LIST(APPEND EXTRA_CXX_FLAGS -DNEED_FILL_NEUT_COMMONS)
SET(HAVENEUTCONFIG FALSE)
# We are dealing with shiny NEUT
if(NOT "${NEUTCONFIG}" STREQUAL "NEUTCONFIG-NOTFOUND")
SET(HAVENEUTCONFIG TRUE)
cmessage(STATUS "Found neut-config, using it to determine configuration.")
else()
cmessage(STATUS "Failed to find neut-config, assuming older NEUT build.")
endif()
if(HAVENEUTCONFIG)
execute_process (COMMAND neut-config
--version OUTPUT_VARIABLE NEUT_VER
OUTPUT_STRIP_TRAILING_WHITESPACE)
execute_process (COMMAND neut-config
--incdir OUTPUT_VARIABLE NEUT_INCLUDE_DIRS
OUTPUT_STRIP_TRAILING_WHITESPACE)
execute_process (COMMAND neut-config
--libdir OUTPUT_VARIABLE NEUT_LINK_DIRS
OUTPUT_STRIP_TRAILING_WHITESPACE)
GETLIBDIRS(neut-config --cernflags CERN_LIB_DIR)
LIST(APPEND NEUT_LINK_DIRS ${CERN_LIB_DIR})
GETLIBS(neut-config --cernflags CERN_LIBS)
if(USE_REWEIGHT)
execute_process (COMMAND neut-config
--rwlibflags OUTPUT_VARIABLE NEUT_RWLIBS
OUTPUT_STRIP_TRAILING_WHITESPACE)
GETLIBS(neut-config --rwlibflags NEUT_RWLIBS)
LIST(APPEND NEUT_LIBS ${NEUT_RWLIBS})
else()
GETLIBS(neut-config --libflags NEUT_GENLIBS)
GETLIBS(neut-config --iolibflags NEUT_LIBS)
LIST(APPEND NEUT_LIBS ${NEUT_IOLIBS})
LIST(APPEND NEUT_LIBS ${NEUT_GENLIBS})
endif()
LIST(APPEND NEUT_LIBS ${CERN_LIBS};gfortran)
LIST(APPEND EXTRA_LIBS ${NEUT_LIBS})
string(REPLACE "." "" NEUT_VERSION ${NEUT_VER})
PrefixList(NEUT_INCLUDE_DIRS "-I" ${NEUT_INCLUDE_DIRS})
LIST(APPEND EXTRA_CXX_FLAGS ${NEUT_INCLUDE_DIRS} -D__NEUT_ENABLED__ -D__NEUT_VERSION__=${NEUT_VERSION})
LIST(APPEND EXTRA_LINK_DIRS ${NEUT_LINK_DIRS})
+ LIST(APPEND EXTRA_EXE_FLAGS
+ -fno-pie -fno-PIE -no-pie)
+
cmessage(STATUS "NEUT")
cmessage(STATUS " Version : ${NEUT_VER}")
cmessage(STATUS " Flags : ${NEUT_CXX_FLAGS}")
cmessage(STATUS " Includes : ${NEUT_INCLUDE_DIRS}")
cmessage(STATUS " Link Dirs : ${NEUT_LINK_DIRS}")
cmessage(STATUS " Libs : ${NEUT_LIBS}")
+ cmessage(STATUS " Exe Flags : -fno-pie -fno-PIE -no-pie")
else() # Everything better be set up already
if(NEUT_ROOT STREQUAL "")
cmessage(FATAL_ERROR "Variable NEUT_ROOT is not defined. Please export environment variable NEUT_ROOT or configure with -DNEUT_ROOT=/path/to/NEUT. This must be set to point to a prebuilt NEUT instance.")
endif()
if(CERN STREQUAL "")
cmessage(FATAL_ERROR "Variable CERN is not defined. Please export environment variable CERN or configure with -DCERN=/path/to/CERNLIB. This must be set to point to a prebuilt CERNLIB instance.")
endif()
if(CERN_LEVEL STREQUAL "")
cmessage(FATAL_ERROR "Variable CERN_LEVEL is not defined. Please export environment variable CERN_LEVEL or configure with -DCERN_LEVEL=XXXX (likely to be 2005).")
endif()
if(${NEUT_VERSION} VERSION_LESS 5.4.0)
set(NEUT_LIB_DIR ${NEUT_ROOT}/lib/Linux_pc)
else()
set(NEUT_LIB_DIR ${NEUT_ROOT}/lib)
endif()
set(NEUT_CLASS ${NEUT_ROOT}/src/neutclass)
LIST(APPEND EXTRA_CXX_FLAGS -D__NEUT_ENABLED__ -DNEUT_VERSION=${NEUT_VERSION})
LIST(APPEND EXTRA_CXX_FLAGS
-I${NEUT_ROOT}/include
-I${NEUT_ROOT}/src/neutclass)
+
LIST(APPEND EXTRA_LINK_DIRS
${NEUT_LIB_DIR}
${CERN}/${CERN_LEVEL}/lib)
if(USE_REWEIGHT)
LIST(APPEND EXTRA_CXX_FLAGS
-I${NEUT_ROOT}/src/reweight)
LIST(APPEND EXTRA_LINK_DIRS
${NEUT_ROOT}/src/reweight)
endif()
if(${NEUT_VERSION} VERSION_EQUAL 5.4.2)
LIST(APPEND EXTRA_LIBS
-Wl,--as-needed)
if(USE_REWEIGHT)
LIST(APPEND EXTRA_LIBS
NReWeight)
endif()
LIST(APPEND EXTRA_LIBS
-Wl,--start-group
neutcore_5.4.2
nuccorspl_5.4.2 #typo in NEUT, may hopefully disappear
nuceff_5.4.2
partnuck_5.4.2
skmcsvc_5.4.2
tauola_5.4.2
HT2p2h_5.4.0
N1p1h_5.4.0
-Wl,--end-group
jetset74
pdflib804
mathlib
packlib
pawlib)
LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_COMMON_QEAV)
elseif(${NEUT_VERSION} VERSION_EQUAL 5.4.0)
LIST(APPEND EXTRA_LIBS
-Wl,--as-needed)
if(USE_REWEIGHT)
LIST(APPEND EXTRA_LIBS
NReWeight)
endif()
LIST(APPEND EXTRA_LIBS
-Wl,--start-group
neutcore_5.4.0
nuccorspl_5.4.0 #typo in NEUT, may hopefully disappear
nuceff_5.4.0
partnuck_5.4.0
skmcsvc_5.4.0
tauola_5.4.0
HT2p2h_5.4.0
N1p1h_5.4.0
specfunc_5.4.0
radcorr_5.4.0
gfortran
-Wl,--end-group
jetset74
pdflib804
mathlib
packlib
pawlib)
else()
LIST(APPEND EXTRA_LIBS
-Wl,--as-needed)
if(USE_REWEIGHT)
LIST(APPEND EXTRA_LIBS
NReWeight)
endif()
LIST(APPEND EXTRA_LIBS
-Wl,--start-group
neutcore
nuccorrspl
nuceff
partnuck
skmcsvc
tauola
-Wl,--end-group
jetset74
pdflib804
mathlib
packlib
pawlib)
endif()
set(NEUT_ROOT_LIBS)
LIST(APPEND NEUT_ROOT_LIBS
${NEUT_CLASS}/neutctrl.so
${NEUT_CLASS}/neutfsivert.so)
# Check for new versions of NEUT with NUCLEON FSI
if(EXISTS "${NEUT_CLASS}/neutnucfsistep.so")
set(NEUT_NUCFSI 1)
LIST(APPEND EXTRA_CXX_FLAGS -DNEUT_NUCFSI_ENABLED)
LIST(APPEND NEUT_ROOT_LIBS
${NEUT_CLASS}/neutnucfsistep.so
${NEUT_CLASS}/neutnucfsivert.so
)
endif()
if(${NEUT_VERSION} VERSION_LESS 5.4.0)
LIST(APPEND NEUT_ROOT_LIBS
${NEUT_CLASS}/neutrootTreeSingleton.so)
endif()
LIST(APPEND NEUT_ROOT_LIBS
${NEUT_CLASS}/neutvtx.so
${NEUT_CLASS}/neutfsipart.so
${NEUT_CLASS}/neutpart.so
${NEUT_CLASS}/neutvect.so
)
foreach(OBJ ${NEUT_ROOT_LIBS})
LIST(APPEND EXTRA_SHAREDOBJS ${OBJ})
endforeach()
endif()
diff --git a/cmake/ROOTSetup.cmake b/cmake/ROOTSetup.cmake
index 622aea7..151ce2f 100644
--- a/cmake/ROOTSetup.cmake
+++ b/cmake/ROOTSetup.cmake
@@ -1,172 +1,172 @@
# 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 ( NOT DEFINED ENV{ROOTSYS} )
cmessage (FATAL_ERROR "$ROOTSYS is not defined, please set up ROOT first.")
else()
cmessage(STATUS "Using ROOT installed at $ENV{ROOTSYS}")
set(CMAKE_ROOTSYS $ENV{ROOTSYS})
endif()
# Get cflags from ROOT
execute_process (COMMAND root-config
--cflags OUTPUT_VARIABLE ROOT_CXX_FLAGS_RAW OUTPUT_STRIP_TRAILING_WHITESPACE)
string(REPLACE " " ";" ROOT_CXX_FLAGS "${ROOT_CXX_FLAGS_RAW}")
# Get libdir from ROOT
execute_process (COMMAND root-config
--libdir OUTPUT_VARIABLE ROOT_LIBDIR OUTPUT_STRIP_TRAILING_WHITESPACE)
# Get version from ROOT
execute_process (COMMAND root-config
--version OUTPUT_VARIABLE ROOT_VERSION OUTPUT_STRIP_TRAILING_WHITESPACE)
# Get features from ROOT
execute_process (COMMAND root-config
--features OUTPUT_VARIABLE ROOT_FEATURES OUTPUT_STRIP_TRAILING_WHITESPACE)
LIST(APPEND EXTRA_LINK_DIRS ${ROOT_LIBDIR})
LIST(APPEND ROOT_LIBS
Core
Cint
RIO
XMLIO
Net
Hist
Graf
Graf3d
Gpad
Tree
Rint
Postscript
Matrix
Physics
MathMore
MathCore
Thread
EG
Geom
GenVector)
cmessage(STATUS "Checking ROOT version: ${ROOT_VERSION}")
string(REGEX MATCH "^6.*" ROOTVERSIXMATCH ${ROOT_VERSION})
if(ROOTVERSIXMATCH)
cmessage(STATUS "Using ROOT6, We are essentially flying blind here.")
LIST(REMOVE_ITEM ROOT_LIBS Cint)
- LIST(APPEND EXTRA_CXX_FLAGS -DROOT6_USE_FIT_FITTER_INTERFACE)
+ LIST(APPEND EXTRA_CXX_FLAGS -DROOT6_USE_FIT_FITTER_INTERFACE -DROOT6)
set(USE_ROOT6 True)
else()
string(REGEX MATCH "5.34/([0-9]+)" ROOTVERSMATCH ${ROOT_VERSION})
if(NOT ROOTVERSMATCH OR ${CMAKE_MATCH_1} LESS "19")
cmessage(FATAL_ERROR "ROOT Version: ${ROOT_VERSION} has out of date minimizer interface, but minimizer functionality requested. Please configure with -DUSE_MINIMIZER=FALSE or update to 5.34/19 or greater to enable minimization features.")
endif()
endif()
if(USE_MINIMIZER)
if("${ROOT_FEATURES}" MATCHES "minuit2")
cmessage(STATUS "ROOT built with MINUIT2 support")
LIST(APPEND EXTRA_CXX_FLAGS -D__MINUIT2_ENABLED__)
else()
cmessage(FATAL_ERROR "ROOT built without MINUIT2 support but minimizer functionality requested. Either configure with -DUSE_MINIMIZER=FALSE or use a version of ROOT with MINUIT2 support.")
endif()
endif()
if("${ROOT_FEATURES}" MATCHES "opengl")
cmessage(STATUS "ROOT built with OpenGL support")
LIST(APPEND ROOT_LIBS RGL)
endif()
if(DEFINED NEED_ROOTPYTHIA6 AND NEED_ROOTPYTHIA6)
LIST(APPEND ROOT_LIBS EGPythia6 Pythia6)
endif()
cmessage ( STATUS "[ROOT]: root-config --version: ${ROOT_VERSION} ")
cmessage ( STATUS "[ROOT]: root-config --cflags : ${ROOT_CXX_FLAGS} ")
cmessage ( STATUS "[ROOT]: root-config --libs : ${ROOT_LD_FLAGS} ")
LIST(APPEND EXTRA_CXX_FLAGS ${ROOT_CXX_FLAGS})
#Helper functions for building dictionaries
function(GenROOTDictionary OutputDictName Header LinkDef)
get_directory_property(incdirs INCLUDE_DIRECTORIES)
string(REPLACE ";" ";-I" LISTDIRINCLUDES "-I${incdirs}")
string(REPLACE " " ";" LISTCPPFLAGS "${EXTRA_CXX_FLAGS}")
#ROOT5 CINT cannot handle it.
list(REMOVE_ITEM LISTCPPFLAGS "-std=c++11")
message(STATUS "LISTCPPFLAGS: ${LISTCPPFLAGS}")
message(STATUS "LISTINCLUDES: ${LISTDIRINCLUDES}")
#Learn how to generate the Dict.cxx and Dict.hxx
add_custom_command(
OUTPUT "${OutputDictName}.cxx" "${OutputDictName}.h"
COMMAND rootcint
ARGS -f ${OutputDictName}.cxx -c
-p ${LISTDIRINCLUDES} ${LISTCPPFLAGS} ${Header} ${LinkDef}
DEPENDS ${Header};${LinkDef})
endfunction()
function(BuildROOTProject ProjectName InputFile CommaSeparatedClassesToDump LIBLINKMODE)
string(REPLACE "," ";" HeadersToDump ${CommaSeparatedClassesToDump})
set(OUTPUTFILES ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectSource.cxx
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}LinkDef.h
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectHeaders.h
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectInstances.h)
cmessage(STATUS "As part of ROOT project: ${ProjectName}")
foreach (header ${HeadersToDump})
LIST(APPEND OUTPUTFILES "${CMAKE_BINARY_DIR}/${ProjectName}/${header}.h")
cmessage(STATUS "Will generate: ${CMAKE_BINARY_DIR}/${ProjectName}/${header}.h")
endforeach()
add_custom_command(
OUTPUT ${OUTPUTFILES}
COMMAND ${CMAKE_BINARY_DIR}/src/Utils/DumpROOTClassesFromVector
ARGS ${InputFile}
${CMAKE_BINARY_DIR}/${ProjectName}
${CommaSeparatedClassesToDump}
VERBATIM
DEPENDS DumpROOTClassesFromVector)
add_custom_target(${ProjectName}_sources
DEPENDS ${OUTPUTFILES})
GenROOTDictionary(
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectDict
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectHeaders.h
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}LinkDef.h
)
add_custom_target(${ProjectName}ProjectDict
DEPENDS
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectDict.cxx
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectDict.h )
# add_dependencies(${ProjectName}ProjectDict ${ProjectName}_sources)
#ProjectSource.cxx includes ProjectDict.cxx, so no need to add to compilation.
set(ROAA_SOURCEFILES
${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectSource.cxx)
add_library(${ProjectName} ${LIBLINKMODE} ${ROAA_SOURCEFILES})
add_dependencies(${ProjectName} ${ProjectName}ProjectDict)
endfunction()
diff --git a/cmake/c++CompilerSetup.cmake b/cmake/c++CompilerSetup.cmake
index bc6b473..11f6f33 100644
--- a/cmake/c++CompilerSetup.cmake
+++ b/cmake/c++CompilerSetup.cmake
@@ -1,121 +1,122 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
if(USE_OMP)
LIST(APPEND EXTRA_CXX_FLAGS -fopenmp)
endif()
if(USE_DYNSAMPLES)
LIST(APPEND EXTRA_LIBS dl)
LIST(APPEND EXTRA_CXX_FLAGS -D__USE_DYNSAMPLES__)
endif()
set(CXX_WARNINGS -Wall )
cmessage(DEBUG "EXTRA_CXX_FLAGS: ${EXTRA_CXX_FLAGS}")
string(REPLACE ";" " " STR_EXTRA_CXX_FLAGS "${EXTRA_CXX_FLAGS}")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${STR_EXTRA_CXX_FLAGS} ${CXX_WARNINGS}")
set(CMAKE_Fortran_FLAGS_RELEASE "-fPIC")
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0")
if(USE_DYNSAMPLES)
set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -fPIC")
set(CMAKE_Fortran_FLAGS_DEBUG "-fPIC")
endif()
set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -fPIC -O3")
if(CMAKE_BUILD_TYPE MATCHES DEBUG)
set(CURRENT_CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS_DEBUG})
elseif(CMAKE_BUILD_TYPE MATCHES RELEASE)
set(CURRENT_CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS_RELEASE})
else()
cmessage(FATAL_ERROR "[ERROR]: Unknown CMAKE_BUILD_TYPE (\"${CMAKE_BUILD_TYPE}\"): Should be \"DEBUG\" or \"RELEASE\".")
endif()
SET(STR_EXTRA_LINK_DIRS)
if(NOT EXTRA_LINK_DIRS STREQUAL "")
string(REPLACE ";" " -L" STR_EXTRA_LINK_DIRS "-L${EXTRA_LINK_DIRS}")
endif()
SET(STR_EXTRA_LIBS)
if(NOT EXTRA_LIBS STREQUAL "")
SET(STR_EXTRA_LIBS_NO_SCRUB_LINKOPTS)
string(REPLACE ";" " -l" STR_EXTRA_LIBS_NO_SCRUB_LINKOPTS "-l${EXTRA_LIBS}")
string(REPLACE "-l-" "-" STR_EXTRA_LIBS ${STR_EXTRA_LIBS_NO_SCRUB_LINKOPTS})
endif()
SET(STR_EXTRA_SHAREDOBJS)
if(NOT EXTRA_SHAREDOBJS STREQUAL "")
string(REPLACE ";" " " STR_EXTRA_SHAREDOBJS "${EXTRA_SHAREDOBJS}")
endif()
SET(STR_EXTRA_LINK_FLAGS)
if(NOT EXTRA_LINK_FLAGS STREQUAL "")
string(REPLACE ";" " " STR_EXTRA_LINK_FLAGS "${EXTRA_LINK_FLAGS}")
endif()
cmessage(DEBUG "EXTRA_LINK_DIRS: ${STR_EXTRA_LINK_DIRS}")
cmessage(DEBUG "EXTRA_LIBS: ${STR_EXTRA_LIBS}")
cmessage(DEBUG "EXTRA_SHAREDOBJS: ${STR_EXTRA_SHAREDOBJS}")
cmessage(DEBUG "EXTRA_LINK_FLAGS: ${STR_EXTRA_LINK_FLAGS}")
if(NOT STR_EXTRA_LINK_DIRS STREQUAL "" AND NOT STR_EXTRA_LIBS STREQUAL "")
SET(CMAKE_DEPENDLIB_FLAGS "${STR_EXTRA_LINK_DIRS} ${STR_EXTRA_LIBS}")
endif()
if(NOT EXTRA_SHAREDOBJS STREQUAL "")
if(NOT STR_EXTRA_LINK_FLAGS STREQUAL "")
SET(STR_EXTRA_LINK_FLAGS "${STR_EXTRA_SHAREDOBJS} ${STR_EXTRA_LINK_FLAGS}")
else()
SET(STR_EXTRA_LINK_FLAGS "${STR_EXTRA_SHAREDOBJS}")
endif()
endif()
if(NOT EXTRA_LINK_FLAGS STREQUAL "")
if(NOT CMAKE_LINK_FLAGS STREQUAL "")
SET(CMAKE_LINK_FLAGS "${CMAKE_LINK_FLAGS} ${STR_EXTRA_LINK_FLAGS}")
else()
SET(CMAKE_LINK_FLAGS "${STR_EXTRA_LINK_FLAGS}")
endif()
endif()
if(USE_OMP)
cmessage(FATAL_ERROR "No OMP features currently enabled so this is a FATAL_ERROR to let you know that you don't gain anything with this declaration.")
endif()
+string(REPLACE ";" " " NUIS_EXE_FLAGS "${EXTRA_EXE_FLAGS}")
string(STRIP ${CMAKE_CXX_FLAGS} CMAKE_CXX_FLAGS)
string(STRIP ${CMAKE_CXX_FLAGS_RELEASE} CMAKE_CXX_FLAGS_RELEASE)
string(STRIP ${CMAKE_CXX_FLAGS_DEBUG} CMAKE_CXX_FLAGS_DEBUG)
string(STRIP ${CMAKE_LINK_FLAGS} CMAKE_LINK_FLAGS)
string(STRIP ${CMAKE_DEPENDLIB_FLAGS} CMAKE_DEPENDLIB_FLAGS)
if (VERBOSE)
cmessage (STATUS "C++ Compiler : ${CXX_COMPILER_NAME}")
cmessage (STATUS " flags : ${CMAKE_CXX_FLAGS}")
cmessage (STATUS " Release flags : ${CMAKE_CXX_FLAGS_RELEASE}")
cmessage (STATUS " Debug flags : ${CMAKE_CXX_FLAGS_DEBUG}")
cmessage (STATUS " Link Flags : ${CMAKE_LINK_FLAGS}")
cmessage (STATUS " Lib Flags : ${CMAKE_DEPENDLIB_FLAGS}")
endif()
diff --git a/src/FitBase/StackBase.cxx b/src/FitBase/StackBase.cxx
index 1a721a1..bd6aee0 100644
--- a/src/FitBase/StackBase.cxx
+++ b/src/FitBase/StackBase.cxx
@@ -1,216 +1,250 @@
#include "StackBase.h"
void StackBase::AddMode(std::string name, std::string title, int linecolor,
int linewidth, int fillstyle) {
// int ncur = fAllLabels.size();
fAllLabels.push_back(name);
fAllTitles.push_back(title);
std::vector temp;
temp.push_back(linecolor);
temp.push_back(linewidth);
temp.push_back(fillstyle);
fAllStyles.push_back(temp);
}
void StackBase::FluxUnfold(TH1D *flux, TH1D *events, double scalefactor,
int nevents) {
for (size_t i = 0; i < fAllLabels.size(); i++) {
if (fNDim == 1) {
PlotUtils::FluxUnfoldedScaling((TH1D *)fAllHists[i], flux, events,
scalefactor, nevents);
} else if (fNDim == 2) {
PlotUtils::FluxUnfoldedScaling((TH2D *)fAllHists[i], flux, events,
scalefactor);
}
}
}
void StackBase::AddMode(int index, std::string name, std::string title,
int linecolor, int linewidth, int fillstyle) {
while (fAllLabels.size() <= (UInt_t)index) {
fAllLabels.push_back("");
fAllTitles.push_back("");
fAllStyles.push_back(std::vector(1, 1));
}
fAllLabels[index] = (name);
fAllTitles[index] = (title);
std::vector temp;
temp.push_back(linecolor);
temp.push_back(linewidth);
temp.push_back(fillstyle);
fAllStyles[index] = temp;
}
bool StackBase::IncludeInStack(TH1 *hist) {
if (!FitPar::Config().GetParB("includeemptystackhists") and
hist->Integral() == 0.0)
return false;
return true;
}
bool StackBase::IncludeInStack(int index) { return true; }
void StackBase::SetupStack(TH1 *hist) {
fTemplate = (TH1 *)hist->Clone(fName.c_str());
fTemplate->Reset();
// Determine template dim
fNDim = fTemplate->GetDimension();
for (size_t i = 0; i < fAllLabels.size(); i++) {
fAllHists.push_back(
(TH1 *)fTemplate->Clone((fName + "_" + fAllLabels[i]).c_str()));
}
};
void StackBase::Scale(double sf, std::string opt) {
for (size_t i = 0; i < fAllLabels.size(); i++) {
// std::cout << "Scaling Stack Hist " << i << " by " << sf << std::endl;
fAllHists[i]->Scale(sf, opt.c_str());
}
};
void StackBase::Reset() {
for (size_t i = 0; i < fAllLabels.size(); i++) {
fAllHists[i]->Reset();
}
};
void StackBase::FillStack(int index, double x, double y, double z,
double weight) {
if (index < 0 or (UInt_t) index >= fAllLabels.size()) {
NUIS_ERR(WRN, "Returning Stack Fill Because Range = " << index << " "
- << fAllLabels.size());
+ << fAllLabels.size());
return;
}
if (fNDim == 1)
fAllHists[index]->Fill(x, y);
else if (fNDim == 2) {
// std::cout << "Filling 2D Stack " << index << " " << x << " " << y << " "
// << z << std::endl;
((TH2 *)fAllHists[index])->Fill(x, y, z);
}
else if (fNDim == 3)
((TH3 *)fAllHists[index])->Fill(x, y, z, weight);
}
+void StackBase::SetBinContentStack(int index, int binx, int biny, int binz,
+ double content) {
+ if (index < 0 or (UInt_t) index >= fAllLabels.size()) {
+ NUIS_ERR(WRN, "Returning Stack Fill Because Range = " << index << " "
+ << fAllLabels.size());
+ return;
+ }
+
+ if (fNDim == 1) {
+ fAllHists[index]->SetBinContent(binx, content);
+ } else if (fNDim == 2) {
+ ((TH2 *)fAllHists[index])->SetBinContent(binx, biny, content);
+ } else if (fNDim == 3) {
+ ((TH3 *)fAllHists[index])->SetBinContent(binx, biny, binz, content);
+ }
+}
+
+void StackBase::SetBinErrorStack(int index, int binx, int biny, int binz,
+ double error) {
+ if (index < 0 or (UInt_t) index >= fAllLabels.size()) {
+ NUIS_ERR(WRN, "Returning Stack Fill Because Range = " << index << " "
+ << fAllLabels.size());
+ return;
+ }
+
+ if (fNDim == 1) {
+ fAllHists[index]->SetBinError(binx, error);
+ } else if (fNDim == 2) {
+ ((TH2 *)fAllHists[index])->SetBinError(binx, biny, error);
+ } else if (fNDim == 3) {
+ ((TH3 *)fAllHists[index])->SetBinError(binx, biny, binz, error);
+ }
+}
+
void StackBase::Write() {
THStack *st = new THStack();
// Loop and add all histograms
bool saveseperate = FitPar::Config().GetParB("WriteSeperateStacks");
for (size_t i = 0; i < fAllLabels.size(); i++) {
if (!IncludeInStack(fAllHists[i]))
continue;
if (!IncludeInStack(i))
continue;
fAllHists[i]->SetTitle(fAllTitles[i].c_str());
fAllHists[i]->GetXaxis()->SetTitle(fXTitle.c_str());
fAllHists[i]->GetYaxis()->SetTitle(fYTitle.c_str());
fAllHists[i]->GetZaxis()->SetTitle(fZTitle.c_str());
fAllHists[i]->SetLineColor(fAllStyles[i][0]);
fAllHists[i]->SetLineWidth(fAllStyles[i][1]);
fAllHists[i]->SetFillStyle(fAllStyles[i][2]);
fAllHists[i]->SetFillColor(fAllStyles[i][0]);
if (saveseperate)
fAllHists[i]->Write();
st->Add(fAllHists[i]);
}
st->SetTitle(fTitle.c_str());
st->SetName(fName.c_str());
st->Write();
delete st;
};
void StackBase::Multiply(TH1 *hist) {
for (size_t i = 0; i < fAllLabels.size(); i++) {
fAllHists[i]->Multiply(hist);
}
}
void StackBase::Divide(TH1 *hist) {
for (size_t i = 0; i < fAllLabels.size(); i++) {
fAllHists[i]->Divide(hist);
}
}
void StackBase::Add(TH1 *hist, double scale) {
for (size_t i = 0; i < fAllLabels.size(); i++) {
fAllHists[i]->Add(hist, scale);
}
}
void StackBase::Add(StackBase *hist, double scale) {
if (hist->GetType() != fType) {
NUIS_ERR(WRN, "Trying to add two StackBases of different types!");
NUIS_ERR(WRN, fType << " + " << hist->GetType() << " = Undefined.");
NUIS_ERR(WRN, "Doing nothing...");
return;
}
for (size_t i = 0; i < fAllLabels.size(); i++) {
fAllHists[i]->Add(hist->GetHist(i));
}
}
TH1 *StackBase::GetHist(int entry) { return fAllHists[entry]; }
TH1 *StackBase::GetHist(std::string label) {
TH1 *hist = NULL;
std::vector splitlabels = GeneralUtils::ParseToStr(label, "+");
for (size_t j = 0; j < splitlabels.size(); j++) {
std::string newlabel = splitlabels[j];
for (size_t i = 0; i < fAllLabels.size(); i++) {
if (newlabel == fAllLabels[i]) {
if (!hist)
hist = (TH1 *)fAllHists[i]->Clone();
else
hist->Add(fAllHists[i]);
}
}
}
return hist;
}
THStack StackBase::GetStack() {
THStack st = THStack();
for (size_t i = 0; i < fAllLabels.size(); i++) {
st.Add(fAllHists[i]);
}
return st;
}
void StackBase::AddNewHist(std::string name, TH1 *hist) {
AddMode(fAllLabels.size(), name, hist->GetTitle(), hist->GetLineColor());
fAllHists.push_back((TH1 *)hist->Clone());
}
void StackBase::AddToCategory(std::string name, TH1 *hist) {
for (size_t i = 0; i < fAllLabels.size(); i++) {
if (name == fAllLabels[i]) {
fAllHists[i]->Add(hist);
break;
}
}
}
void StackBase::AddToCategory(int index, TH1 *hist) {
fAllHists[index]->Add(hist);
}
diff --git a/src/FitBase/StackBase.h b/src/FitBase/StackBase.h
index c227813..a2bfe1e 100644
--- a/src/FitBase/StackBase.h
+++ b/src/FitBase/StackBase.h
@@ -1,126 +1,129 @@
#ifndef STACK_BASE_H
#define STACK_BASE_H
-#include "MeasurementVariableBox.h"
#include "FitLogger.h"
#include "GeneralUtils.h"
+#include "MeasurementVariableBox.h"
#include "TH1.h"
#include "TH1D.h"
-#include "TH2D.h"
-#include "THStack.h"
#include "TH2.h"
+#include "TH2D.h"
#include "TH3.h"
+#include "THStack.h"
#include "PlotUtils.h"
class StackBase {
public:
- StackBase() {};
- ~StackBase() {};
-
- virtual void AddMode(std::string name, std::string title,
- int linecolor = 1, int linewidth = 1, int fillstyle = 1001);
- virtual void AddMode(int index, std::string name, std::string title,
- int linecolor = 1, int linewidth = 1, int fillstyle = 1001);
-
- virtual bool IncludeInStack(TH1* hist);
- virtual bool IncludeInStack(int index);
-
- virtual void SetupStack(TH1* hist);
- virtual void Scale(double sf, std::string opt = "");
- virtual void FluxUnfold(TH1D* flux, TH1D* events, double scalefactor, int nevents);
- virtual void Reset();
- virtual void FillStack(int index, double x, double y = 1.0, double z = 1.0, double weight = 1.0);
- virtual void Write();
-
- virtual void Add(StackBase* hist, double scale);
- virtual void Add(TH1* hist, double scale);
-
- virtual void AddNewHist(std::string name, TH1* hist);
- virtual void AddToCategory(std::string name, TH1* hist);
- virtual void AddToCategory(int index, TH1* hist);
-
- virtual void Divide(TH1* hist);
- virtual void Multiply(TH1* hist);
- virtual TH1* GetHist(int entry);
- virtual TH1* GetHist(std::string label);
- virtual THStack GetStack();
-
- std::string GetType(){return fType;};
-
- std::string fName;
- std::string fTitle;
- std::string fXTitle;
- std::string fYTitle;
- std::string fZTitle;
- std::string fType;
-
- TH1* fTemplate;
- int fNDim;
-
- // Maps incase we want to be selective
- std::vector< std::vector > fAllStyles;
- std::vector fAllTitles;
- std::vector fAllLabels;
- std::vector fAllHists;
+ StackBase(){};
+ ~StackBase(){};
+
+ virtual void AddMode(std::string name, std::string title, int linecolor = 1,
+ int linewidth = 1, int fillstyle = 1001);
+ virtual void AddMode(int index, std::string name, std::string title,
+ int linecolor = 1, int linewidth = 1,
+ int fillstyle = 1001);
+
+ virtual bool IncludeInStack(TH1 *hist);
+ virtual bool IncludeInStack(int index);
+
+ virtual void SetupStack(TH1 *hist);
+ virtual void Scale(double sf, std::string opt = "");
+ virtual void FluxUnfold(TH1D *flux, TH1D *events, double scalefactor,
+ int nevents);
+ virtual void Reset();
+ virtual void FillStack(int index, double x, double y = 1.0, double z = 1.0,
+ double weight = 1.0);
+ virtual void SetBinContentStack(int index, int binx, int biny, int binz,
+ double content);
+ virtual void SetBinErrorStack(int index, int binx, int biny, int binz,
+ double error);
+ virtual void Write();
+
+ virtual void Add(StackBase *hist, double scale);
+ virtual void Add(TH1 *hist, double scale);
+
+ virtual void AddNewHist(std::string name, TH1 *hist);
+ virtual void AddToCategory(std::string name, TH1 *hist);
+ virtual void AddToCategory(int index, TH1 *hist);
+
+ virtual void Divide(TH1 *hist);
+ virtual void Multiply(TH1 *hist);
+ virtual TH1 *GetHist(int entry);
+ virtual TH1 *GetHist(std::string label);
+ virtual THStack GetStack();
+
+ std::string GetType() { return fType; };
+
+ std::string fName;
+ std::string fTitle;
+ std::string fXTitle;
+ std::string fYTitle;
+ std::string fZTitle;
+ std::string fType;
+
+ TH1 *fTemplate;
+ int fNDim;
+
+ // Maps incase we want to be selective
+ std::vector > fAllStyles;
+ std::vector fAllTitles;
+ std::vector fAllLabels;
+ std::vector fAllHists;
};
-
/*
class NuSpeciesStack : public StackBase {
public:
- SetupStack(TH1* hist) {
- AddMode("numu", "numu", kBlue, 2, 3004);
- AddMode("numubar", "numubar", kRed, 2, 3004 );
- AddMode("nue", "nue", kGreen, 2, 3004 );
- StackBase::SetupStack(hist);
- };
-
- void NuSpeciesStack::FillStack(int species, double x, double y = 1.0, double z = 1.0, double weight = 1.0) {
- Stackbase::FillStack(ConvertSpeciesToIndex(species), x, y, z, weight);
- }
-
- int ConvertSpeciesToIndex(int species) {
- switch (species) {
- case 14: return 0;
- case -14: return 1;
- case 11: return 2;
- default: return -1;
- }
- };
+ SetupStack(TH1* hist) {
+ AddMode("numu", "numu", kBlue, 2, 3004);
+ AddMode("numubar", "numubar", kRed, 2, 3004 );
+ AddMode("nue", "nue", kGreen, 2, 3004 );
+ StackBase::SetupStack(hist);
+ };
+
+ void NuSpeciesStack::FillStack(int species, double x, double y = 1.0,
+double z = 1.0, double weight = 1.0) {
+ Stackbase::FillStack(ConvertSpeciesToIndex(species), x, y, z,
+weight);
+ }
+
+ int ConvertSpeciesToIndex(int species) {
+ switch (species) {
+ case 14: return 0;
+ case -14: return 1;
+ case 11: return 2;
+ default: return -1;
+ }
+ };
};
class TargetStack : public StackBase {
public:
- SetupStack(TH1* hist) {
- AddMode("C", "C", kBlue, 2, 3004);
- AddMode("H", "H", kRed, 2, 3004 );
- AddMode("O", "O", kGreen, 2, 3004 );
- StackBase::SetupStack(hist);
- };
-
- void NuSpeciesStack::FillStack(int species, double x,
- double y = 1.0, double z = 1.0,
- double weight = 1.0) {
- Stackbase::FillStack(ConvertTargetToIndex(species), x, y, z, weight);
- }
-
- int ConvertTargetToIndex(int target) {
- switch (species) {
- case 1000010010: return 0;
- case 1000: return 1;
- case 1000: return 2;
- default: return -1;
- }
- };
+ SetupStack(TH1* hist) {
+ AddMode("C", "C", kBlue, 2, 3004);
+ AddMode("H", "H", kRed, 2, 3004 );
+ AddMode("O", "O", kGreen, 2, 3004 );
+ StackBase::SetupStack(hist);
+ };
+
+ void NuSpeciesStack::FillStack(int species, double x,
+ double y = 1.0, double z = 1.0,
+ double weight = 1.0) {
+ Stackbase::FillStack(ConvertTargetToIndex(species), x, y, z,
+weight);
+ }
+
+ int ConvertTargetToIndex(int target) {
+ switch (species) {
+ case 1000010010: return 0;
+ case 1000: return 1;
+ case 1000: return 2;
+ default: return -1;
+ }
+ };
}
*/
-
#endif
-
-
-
-
-
diff --git a/src/FitBase/StandardStacks.cxx b/src/FitBase/StandardStacks.cxx
index da2c562..6088382 100644
--- a/src/FitBase/StandardStacks.cxx
+++ b/src/FitBase/StandardStacks.cxx
@@ -1,428 +1,527 @@
#include "StandardStacks.h"
/// Fake stack functions
-FakeStack::FakeStack(TH1* hist) {
- fTemplate = (TH1*)hist;
+FakeStack::FakeStack(TH1 *hist) {
+ fTemplate = (TH1 *)hist;
fNDim = fTemplate->GetDimension();
fFakeType = "TH1";
fTGraphObject = NULL;
fTF1Object = NULL;
}
-
-FakeStack::FakeStack(TF1* f){
+FakeStack::FakeStack(TF1 *f) {
fTemplate = NULL;
fNDim = 1;
fTF1Object = f;
fFakeType = "TF1";
fTGraphObject = NULL;
}
-FakeStack::FakeStack(TGraph* gr){
+FakeStack::FakeStack(TGraph *gr) {
fTemplate = NULL;
fNDim = 1;
fTGraphObject = gr;
fFakeType = "TGRAPH";
fTF1Object = NULL;
}
FakeStack::~FakeStack() {
fTemplate = NULL;
fNDim = 0;
fFakeType = "NULL";
fTF1Object = NULL;
fTGraphObject = NULL;
}
void FakeStack::Fill(double x, double y, double z, double weight) {
- if (fTemplate){
- if (fNDim == 1) fTemplate->Fill(x, y);
- else if (fNDim == 2) ((TH2*)fTemplate)->Fill(x, y, z);
- else if (fNDim == 3) ((TH3*)fTemplate)->Fill(x, y, z, weight);
+ if (fTemplate) {
+ if (fNDim == 1)
+ fTemplate->Fill(x, y);
+ else if (fNDim == 2)
+ ((TH2 *)fTemplate)->Fill(x, y, z);
+ else if (fNDim == 3)
+ ((TH3 *)fTemplate)->Fill(x, y, z, weight);
}
- if (fTGraphObject){
- fTGraphObject->SetPoint( fTGraphObject->GetN(), x, y);
+ if (fTGraphObject) {
+ fTGraphObject->SetPoint(fTGraphObject->GetN(), x, y);
}
}
void FakeStack::Scale(double norm, std::string opt) {
- if (fTemplate){ fTemplate->Scale(norm, opt.c_str()); }
- if (fTGraphObject){
- for (int i = 0; i < fTGraphObject->GetN(); i++){
- const double* x = fTGraphObject->GetX();
- const double* y = fTGraphObject->GetY();
+ if (fTemplate) {
+ fTemplate->Scale(norm, opt.c_str());
+ }
+ if (fTGraphObject) {
+ for (int i = 0; i < fTGraphObject->GetN(); i++) {
+ const double *x = fTGraphObject->GetX();
+ const double *y = fTGraphObject->GetY();
fTGraphObject->SetPoint(i, x[i], y[i] * norm);
}
}
}
void FakeStack::Reset() {
- if (fTemplate){ fTemplate->Reset(); }
+ if (fTemplate) {
+ fTemplate->Reset();
+ }
}
void FakeStack::Write() {
- if (fTemplate){ fTemplate->Write(); }
- if (fTF1Object){ fTF1Object->Write(); }
- if (fTGraphObject){ fTGraphObject->Write(); }
+ if (fTemplate) {
+ fTemplate->Write();
+ }
+ if (fTF1Object) {
+ fTF1Object->Write();
+ }
+ if (fTGraphObject) {
+ fTGraphObject->Write();
+ }
}
-
-
/// TrueModeStack Functions
-TrueModeStack::TrueModeStack(std::string name, std::string title, TH1* hist) {
+TrueModeStack::TrueModeStack(std::string name, std::string title, TH1 *hist) {
fName = name;
fTitle = title;
// CC
- AddMode(0, "CCQE", "CCQE", kBlue, 2, 1001);
- AddMode(1, "CC2p2h", "2p2h", kRed, 2, 1001);
- AddMode(2, "CC1piponp", "CC1#pi^{+} on p", kGreen, 2, 1001);
- AddMode(3, "CC1pi0onn", "CC1#pi^{0} on n", kGreen + 3, 2, 1001);
- AddMode(4, "CC1piponn", "CC1#pi^{+} on n", kGreen - 2, 2, 1001);
- AddMode(5, "CCcoh", "CC coherent", kBlue, 2, 1001);
- AddMode(6, "CC1gamma", "CC1#gamma", kMagenta, 2, 1001);
+ AddMode(0, "CCQE", "CCQE", kBlue, 2, 1001);
+ AddMode(1, "CC2p2h", "2p2h", kRed, 2, 1001);
+ AddMode(2, "CC1piponp", "CC1#pi^{+} on p", kGreen, 2, 1001);
+ AddMode(3, "CC1pi0onn", "CC1#pi^{0} on n", kGreen + 3, 2, 1001);
+ AddMode(4, "CC1piponn", "CC1#pi^{+} on n", kGreen - 2, 2, 1001);
+ AddMode(5, "CCcoh", "CC coherent", kBlue, 2, 1001);
+ AddMode(6, "CC1gamma", "CC1#gamma", kMagenta, 2, 1001);
AddMode(7, "CCMultipi", "Multi #pi (1.3 < W < 2.0)", kYellow, 2, 1001);
- AddMode(8, "CC1eta", "CC1#eta^{0} on n", kYellow - 2, 2, 1001);
- AddMode(9, "CC1lamkp", "CC1#Lambda1K^{+}", kYellow - 6, 2, 1001);
- AddMode(10, "CCDIS", "DIS (W > 2.0)", kRed, 2, 1001);
+ AddMode(8, "CC1eta", "CC1#eta^{0} on n", kYellow - 2, 2, 1001);
+ AddMode(9, "CC1lamkp", "CC1#Lambda1K^{+}", kYellow - 6, 2, 1001);
+ AddMode(10, "CCDIS", "DIS (W > 2.0)", kRed, 2, 1001);
// NC
- AddMode(11, "NC1pi0onn", "NC1#pi^{0} on n", kBlue, 2, 3004);
+ AddMode(11, "NC1pi0onn", "NC1#pi^{0} on n", kBlue, 2, 3004);
AddMode(12, "NC1pi0onp", "NC1#pi^{0} on p", kBlue + 3, 2, 3004);
- AddMode(13, "NC1pimonn", "NC1#pi^{-} on n", kBlue - 2, 2, 3004);
+ AddMode(13, "NC1pimonn", "NC1#pi^{-} on n", kBlue - 2, 2, 3004);
AddMode(14, "NC1piponp", "NC1#pi^{+} on p", kBlue - 8, 2, 3004);
- AddMode(15, "NCcoh", "NC Coherent", kBlue + 8, 2, 3004);
- AddMode(16, "NC1gammaonn", "NC1#gamma on n", kMagenta, 2, 3004);
- AddMode(17, "NC1gammaonp", "NC1#gamma on p", kMagenta - 10, 2, 3004);
+ AddMode(15, "NCcoh", "NC Coherent", kBlue + 8, 2, 3004);
+ AddMode(16, "NC1gammaonn", "NC1#gamma on n", kMagenta, 2, 3004);
+ AddMode(17, "NC1gammaonp", "NC1#gamma on p", kMagenta - 10, 2, 3004);
AddMode(18, "NCMultipi", "Multi #pi (1.3 < W < 2.0)", kBlue - 10, 2, 3004);
AddMode(19, "NC1etaonn", "NC1#eta^{0} on n", kYellow - 2, 2, 3004);
AddMode(20, "NC1etaonp", "NC1#eta^{0} on p", kYellow - 4, 2, 3004);
- AddMode(21, "NC1kamk0", "NC1#Lambda1K^{0} on n", kYellow - 6, 2, 3004);
- AddMode(22, "NC1lamkp", "NC1#Lambda1K^{+}", kYellow - 10, 2, 3004);
+ AddMode(21, "NC1kamk0", "NC1#Lambda1K^{0} on n", kYellow - 6, 2, 3004);
+ AddMode(22, "NC1lamkp", "NC1#Lambda1K^{+}", kYellow - 10, 2, 3004);
AddMode(23, "NCDIS", "DIS (W > 2.0)", kRed, 2, 3004);
AddMode(24, "NCELonp", "NCEL on p", kBlack, 2, 3004);
AddMode(25, "NCELonn", "NCEL on n", kGray, 2, 3004);
- AddMode(26, "NC2p2h", "NC 2p2h", kRed+1, 2, 3004);
+ AddMode(26, "NC2p2h", "NC 2p2h", kRed + 1, 2, 3004);
// Undefined
AddMode(27, "UNDEFINED", "Undefined", kRed + 1, 2, 3000);
StackBase::SetupStack(hist);
};
int TrueModeStack::ConvertModeToIndex(int mode) {
// std::cout << "Converting Mode " << (mode) << std::endl;
switch (abs(mode)) {
- case 1: return 0; // CCQE
- case 2: return 1; // CC2p2h
- case 11: return 2; // CC1piponp
- case 12: return 3; // CC1pi0onn
- case 13: return 4; // CC1piponn
- case 16: return 5; // CCcoh
- case 17: return 6; // CC1gamma
- case 21: return 7; // CCMultipi
- case 22: return 8; // CC1eta
- case 23: return 9; // CC1lamkp
- case 26: return 10; // CCDIS
-
- case 31: return 11; // NC1pi0onn
- case 32: return 12; // NC1pi0onp
- case 33: return 13; // NC1pimonn
- case 34: return 14; // NC1piponp
- case 36: return 15; // NCcoh
- case 38: return 16; // NC1gammaonn
- case 39: return 17; // NC1gammaonp
- case 41: return 18; // NC Multipi
- case 42: return 19; // NC1etaonn
- case 43: return 20; // NC1etaonp
- case 44: return 21; // NC1kamk0
- case 45: return 22; // NC1lamkp
- case 46: return 23; // NCDIS
- case 51: return 24; // NCEL on p
- case 52: return 25; // NCEL on n
- case 53: return 26; // NC 2p2h
- default: return 27; // Undefined
+ case 1:
+ return 0; // CCQE
+ case 2:
+ return 1; // CC2p2h
+ case 11:
+ return 2; // CC1piponp
+ case 12:
+ return 3; // CC1pi0onn
+ case 13:
+ return 4; // CC1piponn
+ case 16:
+ return 5; // CCcoh
+ case 17:
+ return 6; // CC1gamma
+ case 21:
+ return 7; // CCMultipi
+ case 22:
+ return 8; // CC1eta
+ case 23:
+ return 9; // CC1lamkp
+ case 26:
+ return 10; // CCDIS
+
+ case 31:
+ return 11; // NC1pi0onn
+ case 32:
+ return 12; // NC1pi0onp
+ case 33:
+ return 13; // NC1pimonn
+ case 34:
+ return 14; // NC1piponp
+ case 36:
+ return 15; // NCcoh
+ case 38:
+ return 16; // NC1gammaonn
+ case 39:
+ return 17; // NC1gammaonp
+ case 41:
+ return 18; // NC Multipi
+ case 42:
+ return 19; // NC1etaonn
+ case 43:
+ return 20; // NC1etaonp
+ case 44:
+ return 21; // NC1kamk0
+ case 45:
+ return 22; // NC1lamkp
+ case 46:
+ return 23; // NCDIS
+ case 51:
+ return 24; // NCEL on p
+ case 52:
+ return 25; // NCEL on n
+ case 53:
+ return 26; // NC 2p2h
+ default:
+ return 27; // Undefined
}
};
-void TrueModeStack::Fill(int mode, double x, double y, double z, double weight) {
+void TrueModeStack::Fill(int mode, double x, double y, double z,
+ double weight) {
StackBase::FillStack(ConvertModeToIndex(mode), x, y, z, weight);
};
-void TrueModeStack::Fill(FitEvent* evt, double x, double y, double z, double weight) {
+void TrueModeStack::SetBinContent(int mode, int binx, int biny, int binz, double content) {
+ StackBase::SetBinContentStack(ConvertModeToIndex(mode), binx, biny, binz,
+ content);
+}
+void TrueModeStack::SetBinError(int mode, int binx, int biny, int binz, double error) {
+ StackBase::SetBinErrorStack(ConvertModeToIndex(mode), binx, biny, binz,
+ error);
+}
+
+void TrueModeStack::Fill(FitEvent *evt, double x, double y, double z,
+ double weight) {
StackBase::FillStack(ConvertModeToIndex(evt->Mode), x, y, z, weight);
};
-void TrueModeStack::Fill(BaseFitEvt* evt, double x, double y, double z, double weight) {
+void TrueModeStack::Fill(BaseFitEvt *evt, double x, double y, double z,
+ double weight) {
StackBase::FillStack(ConvertModeToIndex(evt->Mode), x, y, z, weight);
};
-
-
/// TrueModeStack Functions
-NuNuBarTrueModeStack::NuNuBarTrueModeStack(std::string name, std::string title, TH1* hist) {
+NuNuBarTrueModeStack::NuNuBarTrueModeStack(std::string name, std::string title,
+ TH1 *hist) {
fName = name;
fTitle = title;
// Neutrino
// CC
- AddMode(0, "NU_CCQE", "#nu CCQE", kBlue, 2, 1001);
- AddMode(1, "NU_CC2p2h", "#nu 2p2h", kRed, 2, 1001);
- AddMode(2, "NU_CC1piponp", "#nu CC1#pi^{+} on p", kGreen, 2, 1001);
- AddMode(3, "NU_CC1pi0onn", "#nu CC1#pi^{0} on n", kGreen + 3, 2, 1001);
- AddMode(4, "NU_CC1piponn", "#nu CC1#pi^{+} on n", kGreen - 2, 2, 1001);
- AddMode(5, "NU_CCcoh", "#nu CC coherent", kBlue, 2, 1001);
- AddMode(6, "NU_CC1gamma", "#nu CC1#gamma", kMagenta, 2, 1001);
- AddMode(7, "NU_CCMultipi", "#nu Multi #pi (1.3 < W < 2.0)", kYellow, 2, 1001);
- AddMode(8, "NU_CC1eta", "#nu CC1#eta^{0} on n", kYellow - 2, 2, 1001);
- AddMode(9, "NU_CC1lamkp", "#nu CC1#Lambda1K^{+}", kYellow - 6, 2, 1001);
- AddMode(10, "NU_CCDIS", "#nu DIS (W > 2.0)", kRed, 2, 1001);
+ AddMode(0, "NU_CCQE", "#nu CCQE", kBlue, 2, 1001);
+ AddMode(1, "NU_CC2p2h", "#nu 2p2h", kRed, 2, 1001);
+ AddMode(2, "NU_CC1piponp", "#nu CC1#pi^{+} on p", kGreen, 2, 1001);
+ AddMode(3, "NU_CC1pi0onn", "#nu CC1#pi^{0} on n", kGreen + 3, 2, 1001);
+ AddMode(4, "NU_CC1piponn", "#nu CC1#pi^{+} on n", kGreen - 2, 2, 1001);
+ AddMode(5, "NU_CCcoh", "#nu CC coherent", kBlue, 2, 1001);
+ AddMode(6, "NU_CC1gamma", "#nu CC1#gamma", kMagenta, 2, 1001);
+ AddMode(7, "NU_CCMultipi", "#nu Multi #pi (1.3 < W < 2.0)", kYellow, 2, 1001);
+ AddMode(8, "NU_CC1eta", "#nu CC1#eta^{0} on n", kYellow - 2, 2, 1001);
+ AddMode(9, "NU_CC1lamkp", "#nu CC1#Lambda1K^{+}", kYellow - 6, 2, 1001);
+ AddMode(10, "NU_CCDIS", "#nu DIS (W > 2.0)", kRed, 2, 1001);
// NC
- AddMode(11, "NU_NC1pi0onn", "#nu NC1#pi^{0} on n", kBlue, 2, 3004);
- AddMode(12, "NU_NC1pi0onp", "#nu NC1#pi^{0} on p", kBlue + 3, 2, 3004);
- AddMode(13, "NU_NC1pimonn", "#nu NC1#pi^{-} on n", kBlue - 2, 2, 3004);
- AddMode(14, "NU_NC1piponp", "#nu NC1#pi^{+} on p", kBlue - 8, 2, 3004);
- AddMode(15, "NU_NCcoh", "#nu NC Coherent", kBlue + 8, 2, 3004);
- AddMode(16, "NU_NC1gammaonn", "#nu NC1#gamma on n", kMagenta, 2, 3004);
- AddMode(17, "NU_NC1gammaonp", "#nu NC1#gamma on p", kMagenta - 10, 2, 3004);
- AddMode(18, "NU_NCMultipi", "#nu Multi #pi (1.3 < W < 2.0)", kBlue - 10, 2, 3004);
- AddMode(19, "NU_NC1etaonn", "#nu NC1#eta^{0} on n", kYellow - 2, 2, 3004);
- AddMode(20, "NU_NC1etaonp", "#nu NC1#eta^{0} on p", kYellow - 4, 2, 3004);
- AddMode(21, "NU_NC1kamk0", "#nu NC1#Lambda1K^{0} on n", kYellow - 6, 2, 3004);
- AddMode(22, "NU_NC1lamkp", "#nu NC1#Lambda1K^{+}", kYellow - 10, 2, 3004);
- AddMode(23, "NU_NCDIS", "#nu DIS (W > 2.0)", kRed, 2, 3004);
- AddMode(24, "NU_NCELonp", "#nu NCEL on p", kBlack, 2, 3004);
- AddMode(25, "NU_NCELonn", "#nu NCEL on n", kGray, 2, 3004);
+ AddMode(11, "NU_NC1pi0onn", "#nu NC1#pi^{0} on n", kBlue, 2, 3004);
+ AddMode(12, "NU_NC1pi0onp", "#nu NC1#pi^{0} on p", kBlue + 3, 2, 3004);
+ AddMode(13, "NU_NC1pimonn", "#nu NC1#pi^{-} on n", kBlue - 2, 2, 3004);
+ AddMode(14, "NU_NC1piponp", "#nu NC1#pi^{+} on p", kBlue - 8, 2, 3004);
+ AddMode(15, "NU_NCcoh", "#nu NC Coherent", kBlue + 8, 2, 3004);
+ AddMode(16, "NU_NC1gammaonn", "#nu NC1#gamma on n", kMagenta, 2, 3004);
+ AddMode(17, "NU_NC1gammaonp", "#nu NC1#gamma on p", kMagenta - 10, 2, 3004);
+ AddMode(18, "NU_NCMultipi", "#nu Multi #pi (1.3 < W < 2.0)", kBlue - 10, 2,
+ 3004);
+ AddMode(19, "NU_NC1etaonn", "#nu NC1#eta^{0} on n", kYellow - 2, 2, 3004);
+ AddMode(20, "NU_NC1etaonp", "#nu NC1#eta^{0} on p", kYellow - 4, 2, 3004);
+ AddMode(21, "NU_NC1kamk0", "#nu NC1#Lambda1K^{0} on n", kYellow - 6, 2, 3004);
+ AddMode(22, "NU_NC1lamkp", "#nu NC1#Lambda1K^{+}", kYellow - 10, 2, 3004);
+ AddMode(23, "NU_NCDIS", "#nu DIS (W > 2.0)", kRed, 2, 3004);
+ AddMode(24, "NU_NCELonp", "#nu NCEL on p", kBlack, 2, 3004);
+ AddMode(25, "NU_NCELonn", "#nu NCEL on n", kGray, 2, 3004);
// Undefined
- AddMode(26, "NU_UNDEFINED", "#nu Undefined", kRed + 2, 2, 3000);
+ AddMode(26, "NU_UNDEFINED", "#nu Undefined", kRed + 2, 2, 3000);
// CC
- AddMode(27, "ANTINU_CCQE", "#bar{#nu} CCQE", kBlue, 2, 1001);
- AddMode(28, "ANTINU_CC2p2h", "#bar{#nu} 2p2h", kRed, 2, 1001);
- AddMode(29, "ANTINU_CC1piponp", "#bar{#nu} CC1#pi^{+} on p", kGreen, 2, 1001);
- AddMode(30, "ANTINU_CC1pi0onn", "#bar{#nu} CC1#pi^{0} on n", kGreen + 3, 2, 1001);
- AddMode(31, "ANTINU_CC1piponn", "#bar{#nu} CC1#pi^{+} on n", kGreen - 2, 2, 1001);
- AddMode(32, "ANTINU_CCcoh", "#bar{#nu} CC coherent", kBlue, 2, 1001);
- AddMode(33, "ANTINU_CC1gamma", "#bar{#nu} CC1#gamma", kMagenta, 2, 1001);
- AddMode(34, "ANTINU_CCMultipi", "#bar{#nu} Multi #pi (1.3 < W < 2.0)", kYellow, 2, 1001);
- AddMode(35, "ANTINU_CC1eta", "#bar{#nu} CC1#eta^{0} on n", kYellow - 2, 2, 1001);
- AddMode(36, "ANTINU_CC1lamkp", "#bar{#nu} CC1#Lambda1K^{+}", kYellow - 6, 2, 1001);
- AddMode(37, "ANTINU_CCDIS", "#bar{#nu} DIS (W > 2.0)", kRed, 2, 1001);
+ AddMode(27, "ANTINU_CCQE", "#bar{#nu} CCQE", kBlue, 2, 1001);
+ AddMode(28, "ANTINU_CC2p2h", "#bar{#nu} 2p2h", kRed, 2, 1001);
+ AddMode(29, "ANTINU_CC1piponp", "#bar{#nu} CC1#pi^{+} on p", kGreen, 2, 1001);
+ AddMode(30, "ANTINU_CC1pi0onn", "#bar{#nu} CC1#pi^{0} on n", kGreen + 3, 2,
+ 1001);
+ AddMode(31, "ANTINU_CC1piponn", "#bar{#nu} CC1#pi^{+} on n", kGreen - 2, 2,
+ 1001);
+ AddMode(32, "ANTINU_CCcoh", "#bar{#nu} CC coherent", kBlue, 2, 1001);
+ AddMode(33, "ANTINU_CC1gamma", "#bar{#nu} CC1#gamma", kMagenta, 2, 1001);
+ AddMode(34, "ANTINU_CCMultipi", "#bar{#nu} Multi #pi (1.3 < W < 2.0)",
+ kYellow, 2, 1001);
+ AddMode(35, "ANTINU_CC1eta", "#bar{#nu} CC1#eta^{0} on n", kYellow - 2, 2,
+ 1001);
+ AddMode(36, "ANTINU_CC1lamkp", "#bar{#nu} CC1#Lambda1K^{+}", kYellow - 6, 2,
+ 1001);
+ AddMode(37, "ANTINU_CCDIS", "#bar{#nu} DIS (W > 2.0)", kRed, 2, 1001);
// NC
- AddMode(38, "ANTINU_NC1pi0onn", "#bar{#nu} NC1#pi^{0} on n", kBlue, 2, 3004);
- AddMode(39, "ANTINU_NC1pi0onp", "#bar{#nu} NC1#pi^{0} on p", kBlue + 3, 2, 3004);
- AddMode(40, "ANTINU_NC1pimonn", "#bar{#nu} NC1#pi^{-} on n", kBlue - 2, 2, 3004);
- AddMode(41, "ANTINU_NC1piponp", "#bar{#nu} NC1#pi^{+} on p", kBlue - 8, 2, 3004);
- AddMode(42, "ANTINU_NCcoh", "#bar{#nu} NC Coherent", kBlue + 8, 2, 3004);
- AddMode(43, "ANTINU_NC1gammaonn", "#bar{#nu} NC1#gamma on n", kMagenta, 2, 3004);
- AddMode(44, "ANTINU_NC1gammaonp", "#bar{#nu} NC1#gamma on p", kMagenta - 10, 2, 3004);
- AddMode(45, "ANTINU_NCMultipi", "#bar{#nu} Multi #pi (1.3 < W < 2.0)", kBlue - 10, 2, 3004);
- AddMode(46, "ANTINU_NC1etaonn", "#bar{#nu} NC1#eta^{0} on n", kYellow - 2, 2, 3004);
- AddMode(47, "ANTINU_NC1etaonp", "#bar{#nu} NC1#eta^{0} on p", kYellow - 4, 2, 3004);
- AddMode(48, "ANTINU_NC1kamk0", "#bar{#nu} NC1#Lambda1K^{0} on n", kYellow - 6, 2, 3004);
- AddMode(49, "ANTINU_NC1lamkp", "#bar{#nu} NC1#Lambda1K^{+}", kYellow - 10, 2, 3004);
- AddMode(50, "ANTINU_NCDIS", "#bar{#nu} DIS (W > 2.0)", kRed, 2, 3004);
- AddMode(51, "ANTINU_NCELonp", "#bar{#nu} NCEL on p", kBlack, 2, 3004);
- AddMode(52, "ANTINU_NCELonn", "#bar{#nu} NCEL on n", kGray, 2, 3004);
+ AddMode(38, "ANTINU_NC1pi0onn", "#bar{#nu} NC1#pi^{0} on n", kBlue, 2, 3004);
+ AddMode(39, "ANTINU_NC1pi0onp", "#bar{#nu} NC1#pi^{0} on p", kBlue + 3, 2,
+ 3004);
+ AddMode(40, "ANTINU_NC1pimonn", "#bar{#nu} NC1#pi^{-} on n", kBlue - 2, 2,
+ 3004);
+ AddMode(41, "ANTINU_NC1piponp", "#bar{#nu} NC1#pi^{+} on p", kBlue - 8, 2,
+ 3004);
+ AddMode(42, "ANTINU_NCcoh", "#bar{#nu} NC Coherent", kBlue + 8, 2, 3004);
+ AddMode(43, "ANTINU_NC1gammaonn", "#bar{#nu} NC1#gamma on n", kMagenta, 2,
+ 3004);
+ AddMode(44, "ANTINU_NC1gammaonp", "#bar{#nu} NC1#gamma on p", kMagenta - 10,
+ 2, 3004);
+ AddMode(45, "ANTINU_NCMultipi", "#bar{#nu} Multi #pi (1.3 < W < 2.0)",
+ kBlue - 10, 2, 3004);
+ AddMode(46, "ANTINU_NC1etaonn", "#bar{#nu} NC1#eta^{0} on n", kYellow - 2, 2,
+ 3004);
+ AddMode(47, "ANTINU_NC1etaonp", "#bar{#nu} NC1#eta^{0} on p", kYellow - 4, 2,
+ 3004);
+ AddMode(48, "ANTINU_NC1kamk0", "#bar{#nu} NC1#Lambda1K^{0} on n", kYellow - 6,
+ 2, 3004);
+ AddMode(49, "ANTINU_NC1lamkp", "#bar{#nu} NC1#Lambda1K^{+}", kYellow - 10, 2,
+ 3004);
+ AddMode(50, "ANTINU_NCDIS", "#bar{#nu} DIS (W > 2.0)", kRed, 2, 3004);
+ AddMode(51, "ANTINU_NCELonp", "#bar{#nu} NCEL on p", kBlack, 2, 3004);
+ AddMode(52, "ANTINU_NCELonn", "#bar{#nu} NCEL on n", kGray, 2, 3004);
// Undefined
- AddMode(53, "NU_UNDEFINED", "#bar{#nu} Undefined", kRed + 2, 2, 3000);
+ AddMode(53, "NU_UNDEFINED", "#bar{#nu} Undefined", kRed + 2, 2, 3000);
// Non Neutrino
AddMode(54, "UNDEFINED", "Non-#nu Undefined", kBlack, 2, 3000);
StackBase::SetupStack(hist);
};
int NuNuBarTrueModeStack::ConvertModeToIndex(int mode) {
switch (abs(mode)) {
- case 1: return 0; // CCQE
- case 2: return 1; // CC2p2h
- case 11: return 2; // CC1piponp
- case 12: return 3; // CC1pi0onn
- case 13: return 4; // CC1piponn
- case 16: return 5; // CCcoh
- case 17: return 6; // CC1gamma
- case 21: return 7; // CCMultipi
- case 22: return 8; // CC1eta
- case 23: return 9; // CC1lamkp
- case 26: return 10; // CCDIS
-
- case 31: return 11; // NC1pi0onn
- case 32: return 12; // NC1pi0onp
- case 33: return 13; // NC1pimonn
- case 34: return 14; // NC1piponp
- case 36: return 15; // NCcoh
- case 38: return 16; // NC1gammaonn
- case 39: return 17; // NC1gammaonp
- case 41: return 18; // NC Multipi
- case 42: return 19; // NC1etaonn
- case 43: return 20; // NC1etaonp
- case 44: return 21; // NC1kamk0
- case 45: return 22; // NC1lamkp
- case 46: return 23; // NCDIS
- case 51: return 24; // NCEL on p
- case 52: return 25; // NCEL on n
- default: return 26; // Undefined
+ case 1:
+ return 0; // CCQE
+ case 2:
+ return 1; // CC2p2h
+ case 11:
+ return 2; // CC1piponp
+ case 12:
+ return 3; // CC1pi0onn
+ case 13:
+ return 4; // CC1piponn
+ case 16:
+ return 5; // CCcoh
+ case 17:
+ return 6; // CC1gamma
+ case 21:
+ return 7; // CCMultipi
+ case 22:
+ return 8; // CC1eta
+ case 23:
+ return 9; // CC1lamkp
+ case 26:
+ return 10; // CCDIS
+
+ case 31:
+ return 11; // NC1pi0onn
+ case 32:
+ return 12; // NC1pi0onp
+ case 33:
+ return 13; // NC1pimonn
+ case 34:
+ return 14; // NC1piponp
+ case 36:
+ return 15; // NCcoh
+ case 38:
+ return 16; // NC1gammaonn
+ case 39:
+ return 17; // NC1gammaonp
+ case 41:
+ return 18; // NC Multipi
+ case 42:
+ return 19; // NC1etaonn
+ case 43:
+ return 20; // NC1etaonp
+ case 44:
+ return 21; // NC1kamk0
+ case 45:
+ return 22; // NC1lamkp
+ case 46:
+ return 23; // NCDIS
+ case 51:
+ return 24; // NCEL on p
+ case 52:
+ return 25; // NCEL on n
+ default:
+ return 26; // Undefined
}
};
-void NuNuBarTrueModeStack::Fill(int species, int mode, double x, double y, double z, double weight) {
+void NuNuBarTrueModeStack::Fill(int species, int mode, double x, double y,
+ double z, double weight) {
int modeindex = ConvertModeToIndex(mode);
int index = 54; // undefined
- if (species == 12 or species == 14 or species == 16) index = modeindex;
- else if (species == -12 or species == -14 or species == -16) index = modeindex + 27;
+ if (species == 12 or species == 14 or species == 16)
+ index = modeindex;
+ else if (species == -12 or species == -14 or species == -16)
+ index = modeindex + 27;
StackBase::FillStack(index, x, y, z, weight);
};
-
-
-
// Species Stack Functions
-BeamSpeciesStack::BeamSpeciesStack(std::string name, std::string title, TH1* hist) {
+BeamSpeciesStack::BeamSpeciesStack(std::string name, std::string title,
+ TH1 *hist) {
fName = name;
fTitle = title;
// charged eptons
AddMode(0, "electron", "e^{-}", kBlue, 2, 1001);
AddMode(1, "positron", "e^{+}", kBlue - 2, 2, 3004);
AddMode(2, "muon", "#mu^{-}", kRed, 2, 1001);
AddMode(3, "antimuon", "#mu^{+}", kRed - 2, 2, 3004);
AddMode(4, "tau", "#tau^{-}", kGreen, 2, 1001);
AddMode(5, "antitau", "#tau^{+}", kGreen - 2, 2, 3004);
// neutrinos
AddMode(6, "nue", "#nu_e", kBlue, 2, 1001);
AddMode(7, "antinue", "#bar{#nu}_e", kBlue - 2, 2, 3004);
AddMode(8, "numu", "#nu_#mu", kRed, 2, 1001);
AddMode(9, "antinumu", "#bar{#nu}_#mu", kRed - 2, 2, 3004);
AddMode(10, "nutau", "#nu_#tau", kGreen, 2, 1001);
AddMode(11, "antinutau", "#bar{#nu}_#tau", kGreen - 2, 2, 3004);
StackBase::SetupStack(hist);
};
int BeamSpeciesStack::ConvertSpeciesToIndex(int species) {
switch (species) {
- case 11: return 0; // e-
- case -11: return 1; // e+
- case 13: return 2; // mu-
- case -13: return 3; // mu+
- case 15: return 4; // tau-
- case -15: return 5; // tau+
- case 12: return 6; // nue
- case -12: return 7; // nuebar
- case 14: return 8; // numu
- case -14: return 9; // numubar
- case 16: return 10; // nutau
- case -16: return 11; //nutaubar
- default: return 12;
+ case 11:
+ return 0; // e-
+ case -11:
+ return 1; // e+
+ case 13:
+ return 2; // mu-
+ case -13:
+ return 3; // mu+
+ case 15:
+ return 4; // tau-
+ case -15:
+ return 5; // tau+
+ case 12:
+ return 6; // nue
+ case -12:
+ return 7; // nuebar
+ case 14:
+ return 8; // numu
+ case -14:
+ return 9; // numubar
+ case 16:
+ return 10; // nutau
+ case -16:
+ return 11; // nutaubar
+ default:
+ return 12;
}
};
-void BeamSpeciesStack::Fill(int species, double x, double y, double z, double weight) {
+void BeamSpeciesStack::Fill(int species, double x, double y, double z,
+ double weight) {
StackBase::FillStack(ConvertSpeciesToIndex(species), x, y, z, weight);
}
-
-
// Target Stack Functions
-TargetTypeStack::TargetTypeStack(std::string name, std::string title, TH1* hist){
+TargetTypeStack::TargetTypeStack(std::string name, std::string title,
+ TH1 *hist) {
fName = name;
fTitle = title;
AddMode(0, "H", "Hydrogen", kBlue, 2, 1001);
AddMode(1, "C", "Carbon", kRed, 2, 1001);
AddMode(2, "O", "Oxygen", kViolet, 2, 1001);
- AddMode(3, "UNDEFINED", "Undefined", kBlack, 2, 1001 );
+ AddMode(3, "UNDEFINED", "Undefined", kBlack, 2, 1001);
StackBase::SetupStack(hist);
}
-void TargetTypeStack::Fill(int pdg, double x, double y, double z, double weight){
+void TargetTypeStack::Fill(int pdg, double x, double y, double z,
+ double weight) {
int index = ConvertPDGToIndex(pdg);
StackBase::FillStack(index, x, y, z, weight);
}
-int TargetTypeStack::ConvertPDGToIndex(int pdg){
- switch(pdg){
- case 1000010010: return 0; // H
- case 1000060120: return 1; // C
- case 1000080160: return 2; // O
- default: return 3; // Undef
+int TargetTypeStack::ConvertPDGToIndex(int pdg) {
+ switch (pdg) {
+ case 1000010010:
+ return 0; // H
+ case 1000060120:
+ return 1; // C
+ case 1000080160:
+ return 2; // O
+ default:
+ return 3; // Undef
}
}
-bool TargetTypeStack::IncludeInStack(TH1* hist){
+bool TargetTypeStack::IncludeInStack(TH1 *hist) {
return (hist->Integral() > 0.0);
}
-
-
-
// CC Topology Stack Functions
-CCTopologyStack::CCTopologyStack(std::string name, std::string title, TH1* hist) {
+CCTopologyStack::CCTopologyStack(std::string name, std::string title,
+ TH1 *hist) {
fName = name;
fTitle = title;
AddMode(0, "CC0pi", "CC-0#pi", kBlue, 2, 1001);
AddMode(1, "CC1pip", "CC-1#pi^{+}", kRed, 2, 1001);
AddMode(2, "CC1pim", "CC-1#pi^{-}", kGreen, 2, 1001);
AddMode(3, "CC1pi0", "CC-1#pi^{0}", kYellow, 2, 1001);
AddMode(4, "CCNpi", "CC-N#pi", kGray, 2, 1001);
AddMode(5, "CCOther", "CC-Other", kViolet, 2, 1001);
- AddMode(6, "NC", "NC", kMagenta, 2, 1001);
+ AddMode(6, "NC", "NC", kMagenta, 2, 1001);
AddMode(7, "UNDEFINED", "Undefined", kBlack, 2, 1001);
-
-
}
-void CCTopologyStack::Fill(FitEvent* evt, double x, double y, double z, double weight) {
+void CCTopologyStack::Fill(FitEvent *evt, double x, double y, double z,
+ double weight) {
int index = GetIndexFromEventParticles(evt);
StackBase::FillStack(index, x, y, z, weight);
}
-int CCTopologyStack::GetIndexFromEventParticles(FitEvent* evt) {
+int CCTopologyStack::GetIndexFromEventParticles(FitEvent *evt) {
int nleptons = evt->NumFSLeptons();
- int npiplus = evt->NumFSParticle(211);
- int npineg = evt->NumFSParticle(-211);
- int npi0 = evt->NumFSParticle(111);
+ int npiplus = evt->NumFSParticle(211);
+ int npineg = evt->NumFSParticle(-211);
+ int npi0 = evt->NumFSParticle(111);
int npions = npiplus + npineg + npi0;
if (nleptons == 1) {
if (npions == 0) {
return 0; // CC0pi
} else if (npions == 1) {
- if (npiplus == 1) return 1; //CC1pi+
- else if (npineg == 1) return 2; //CC1pi-
- else if (npi0 == 1) return 3; //CC1pi0
+ if (npiplus == 1)
+ return 1; // CC1pi+
+ else if (npineg == 1)
+ return 2; // CC1pi-
+ else if (npi0 == 1)
+ return 3; // CC1pi0
} else if (npions > 1) {
return 4; // CCNpi
}
} else if (nleptons > 1) {
return 5; // CCOther
} else if (nleptons < 1) {
return 6;
}
return 7; // Undefined?
-
}
-
-
-
-
-
-
-
-
-
-
-
diff --git a/src/FitBase/StandardStacks.h b/src/FitBase/StandardStacks.h
index 4603608..322832b 100644
--- a/src/FitBase/StandardStacks.h
+++ b/src/FitBase/StandardStacks.h
@@ -1,133 +1,129 @@
#ifndef STANDARD_STACKS_H
#define STANDARD_STACKS_H
-#include "StackBase.h"
-#include "FitEvent.h"
#include "BaseFitEvt.h"
+#include "FitEvent.h"
+#include "StackBase.h"
-/// Single stack class, for consistent handling of TH1* and StackBase* histograms.
+/// Single stack class, for consistent handling of TH1* and StackBase*
+/// histograms.
class FakeStack : public StackBase {
public:
- /// Sets Template to exact pointer to hist without cloning.
- FakeStack(TH1* hist) ;
+ /// Sets Template to exact pointer to hist without cloning.
+ FakeStack(TH1 *hist);
- /// Sets NULL template and saves TF1 object instead
- FakeStack(TF1* f);
+ /// Sets NULL template and saves TF1 object instead
+ FakeStack(TF1 *f);
- /// Sets NULL template and saves TGraph object instead
- FakeStack(TGraph* gr);
+ /// Sets NULL template and saves TGraph object instead
+ FakeStack(TGraph *gr);
- /// Unlinks pointer to original histogram
- ~FakeStack();
+ /// Unlinks pointer to original histogram
+ ~FakeStack();
- /// Fills the normal fTemplate histogram
- void Fill(double x, double y = 1.0, double z = 1.0, double weight = 1.0);
+ /// Fills the normal fTemplate histogram
+ void Fill(double x, double y = 1.0, double z = 1.0, double weight = 1.0);
- /// Scales the normal fTemplate histogram
- void Scale(double norm, std::string opt);
+ /// Scales the normal fTemplate histogram
+ void Scale(double norm, std::string opt);
- /// Resets the normal fTemplate histogram
- void Reset();
+ /// Resets the normal fTemplate histogram
+ void Reset();
- /// Writes the normal fTemplate histogram
- void Write();
+ /// Writes the normal fTemplate histogram
+ void Write();
- std::string fFakeType;
- TGraph* fTGraphObject;
- TF1* fTF1Object;
+ std::string fFakeType;
+ TGraph *fTGraphObject;
+ TF1 *fTF1Object;
};
-
-
-
-
/// True Mode stack, an array of neut true interaction channels
class TrueModeStack : public StackBase {
public:
+ /// Main constructor listing true mode categories.
+ TrueModeStack(std::string name, std::string title, TH1 *hist);
- /// Main constructor listing true mode categories.
- TrueModeStack(std::string name, std::string title, TH1* hist);
-
- /// List to convert Modes to Index.
- /// Should be kept in sync with constructor.
- int ConvertModeToIndex(int mode);
+ /// List to convert Modes to Index.
+ /// Should be kept in sync with constructor.
+ static int ConvertModeToIndex(int mode);
- /// Fill fromgiven mode integer
- void Fill(int mode, double x, double y = 1.0, double z = 1.0, double weight = 1.0);
+ /// Fill fromgiven mode integer
+ void Fill(int mode, double x, double y = 1.0, double z = 1.0,
+ double weight = 1.0);
- /// Extracts Mode from FitEvent and fills
- void Fill(FitEvent* evt, double x, double y = 1.0, double z = 1.0, double weight = 1.0);
+ void SetBinContent(int mode, int binx, int biny, int binz, double content);
+ void SetBinError(int mode, int binx, int biny, int binz, double error);
- /// Extracts Mode from BaseFitEvt
- void Fill(BaseFitEvt* evt, double x, double y = 1.0, double z = 1.0, double weight = 1.0);
+ /// Extracts Mode from FitEvent and fills
+ void Fill(FitEvent *evt, double x, double y = 1.0, double z = 1.0,
+ double weight = 1.0);
+ /// Extracts Mode from BaseFitEvt
+ void Fill(BaseFitEvt *evt, double x, double y = 1.0, double z = 1.0,
+ double weight = 1.0);
};
/// True Mode NuNuBar stack, array of true channels split by nu/nubar
class NuNuBarTrueModeStack : public StackBase {
public:
+ /// Main constructor listing true mode categories.
+ NuNuBarTrueModeStack(std::string name, std::string title, TH1 *hist);
- /// Main constructor listing true mode categories.
- NuNuBarTrueModeStack(std::string name, std::string title, TH1* hist);
+ /// List to convert Modes to Index.
+ /// Should be kept in sync with constructor.
+ int ConvertModeToIndex(int mode);
- /// List to convert Modes to Index.
- /// Should be kept in sync with constructor.
- int ConvertModeToIndex(int mode);
-
- /// Fill fromgiven mode integer
- void Fill(int species, int mode, double x, double y = 1.0, double z = 1.0, double weight = 1.0);
+ /// Fill fromgiven mode integer
+ void Fill(int species, int mode, double x, double y = 1.0, double z = 1.0,
+ double weight = 1.0);
};
/// Species stack to look at contributions from multiple beam leptons
class BeamSpeciesStack : public StackBase {
public:
+ /// main constructor listing beam categories
+ BeamSpeciesStack(std::string name, std::string title, TH1 *hist);
- /// main constructor listing beam categories
- BeamSpeciesStack(std::string name, std::string title, TH1* hist);
-
- /// Fills stack using neutrino species
- void Fill(int species, double x, double y = 1.0, double z = 1.0, double weight = 1.0);
+ /// Fills stack using neutrino species
+ void Fill(int species, double x, double y = 1.0, double z = 1.0,
+ double weight = 1.0);
- /// List converts PDG Beam to index.
- /// Should be kept in sync with constructor.
- int ConvertSpeciesToIndex(int species);
+ /// List converts PDG Beam to index.
+ /// Should be kept in sync with constructor.
+ int ConvertSpeciesToIndex(int species);
};
-
/// Species stack to look at contributions from multiple beam leptons
class TargetTypeStack : public StackBase {
public:
+ /// main constructor listing beam categories
+ TargetTypeStack(std::string name, std::string title, TH1 *hist);
- /// main constructor listing beam categories
- TargetTypeStack(std::string name, std::string title, TH1* hist);
+ /// Fills stack using target pdg
+ void Fill(int pdg, double x, double y = 1.0, double z = 1.0,
+ double weight = 1.0);
- /// Fills stack using target pdg
- void Fill(int pdg, double x, double y = 1.0, double z = 1.0, double weight = 1.0);
+ /// List converts PDG Beam to index.
+ /// Should be kept in sync with constructor.
+ int ConvertPDGToIndex(int pdg);
- /// List converts PDG Beam to index.
- /// Should be kept in sync with constructor.
- int ConvertPDGToIndex(int pdg);
-
- /// Specie empty histograms to not be included in final stack
- bool IncludeInStack(TH1* hist);
+ /// Specie empty histograms to not be included in final stack
+ bool IncludeInStack(TH1 *hist);
};
-
-
/// CC Topology Stack, categories defined by final state particle counts.
class CCTopologyStack : public StackBase {
public:
+ /// main constructor listing beam categories
+ CCTopologyStack(std::string name, std::string title, TH1 *hist);
- /// main constructor listing beam categories
- CCTopologyStack(std::string name, std::string title, TH1* hist);
-
- /// Fills stack using FitEvent
- void Fill(FitEvent* evt, double x, double y = 1.0, double z = 1.0, double weight = 1.0);
-
- /// Extracts index from evt particle counts
- int GetIndexFromEventParticles(FitEvent* evt);
-
+ /// Fills stack using FitEvent
+ void Fill(FitEvent *evt, double x, double y = 1.0, double z = 1.0,
+ double weight = 1.0);
+ /// Extracts index from evt particle counts
+ int GetIndexFromEventParticles(FitEvent *evt);
};
#endif
diff --git a/src/InputHandler/NEUTInputHandler.cxx b/src/InputHandler/NEUTInputHandler.cxx
index 53cc2a9..a5b77dd 100644
--- a/src/InputHandler/NEUTInputHandler.cxx
+++ b/src/InputHandler/NEUTInputHandler.cxx
@@ -1,544 +1,547 @@
#ifdef __NEUT_ENABLED__
#include "NEUTInputHandler.h"
#include "InputUtils.h"
#include "PlotUtils.h"
#include "TTreePerfStats.h"
#include "fsihistC.h"
#include "necardC.h"
#include "nefillverC.h"
#include "neutcrsC.h"
#include "neutfsipart.h"
#include "neutfsivert.h"
#include "neutmodelC.h"
#include "neutparamsC.h"
#include "neutpart.h"
#include "neutrootTreeSingleton.h"
#include "neutvect.h"
#include "neworkC.h"
#include "vcworkC.h"
#include "posinnucC.h"
#ifdef __NEUT_NUCFSI_ENABLED__
#include "neutnucfsistep.h"
#include "neutnucfsivert.h"
#include "nucleonfsihistC.h"
#endif
NEUTGeneratorInfo::~NEUTGeneratorInfo() { DeallocateParticleStack(); }
void NEUTGeneratorInfo::AddBranchesToTree(TTree *tn) {
tn->Branch("NEUTParticleN", fNEUTParticleN, "NEUTParticleN/I");
tn->Branch("NEUTParticleStatusCode", fNEUTParticleStatusCode,
"NEUTParticleStatusCode[NEUTParticleN]/I");
tn->Branch("NEUTParticleAliveCode", fNEUTParticleAliveCode,
"NEUTParticleAliveCode[NEUTParticleN]/I");
}
void NEUTGeneratorInfo::SetBranchesFromTree(TTree *tn) {
tn->SetBranchAddress("NEUTParticleN", &fNEUTParticleN);
tn->SetBranchAddress("NEUTParticleStatusCode", &fNEUTParticleStatusCode);
tn->SetBranchAddress("NEUTParticleAliveCode", &fNEUTParticleAliveCode);
}
void NEUTGeneratorInfo::AllocateParticleStack(int stacksize) {
fNEUTParticleN = 0;
fNEUTParticleStatusCode = new int[stacksize];
fNEUTParticleStatusCode = new int[stacksize];
}
void NEUTGeneratorInfo::DeallocateParticleStack() {
delete fNEUTParticleStatusCode;
delete fNEUTParticleAliveCode;
}
void NEUTGeneratorInfo::FillGeneratorInfo(NeutVect *nevent) {
Reset();
for (int i = 0; i < nevent->Npart(); i++) {
fNEUTParticleStatusCode[i] = nevent->PartInfo(i)->fStatus;
fNEUTParticleAliveCode[i] = nevent->PartInfo(i)->fIsAlive;
fNEUTParticleN++;
}
}
void NEUTGeneratorInfo::Reset() {
for (int i = 0; i < fNEUTParticleN; i++) {
fNEUTParticleStatusCode[i] = -1;
fNEUTParticleAliveCode[i] = 9;
}
fNEUTParticleN = 0;
}
NEUTInputHandler::NEUTInputHandler(std::string const &handle,
std::string const &rawinputs) {
NUIS_LOG(SAM, "Creating NEUTInputHandler : " << handle);
// Run a joint input handling
fName = handle;
// Setup the TChain
fNEUTTree = new TChain("neuttree");
fSaveExtra = FitPar::Config().GetParB("SaveExtraNEUT");
fCacheSize = FitPar::Config().GetParI("CacheSize");
fMaxEvents = FitPar::Config().GetParI("MAXEVENTS");
// Loop over all inputs and grab flux, eventhist, and nevents
std::vector inputs = InputUtils::ParseInputFileList(rawinputs);
for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) {
// Open File for histogram access
TFile *inp_file = new TFile(inputs[inp_it].c_str(), "READ");
if (!inp_file or inp_file->IsZombie()) {
NUIS_ABORT(
"NEUT File IsZombie() at : '"
<< inputs[inp_it] << "'" << std::endl
<< "Check that your file paths are correct and the file exists!"
<< std::endl
<< "$ ls -lh " << inputs[inp_it]);
}
// Get Flux/Event hist
TH1D *fluxhist = (TH1D *)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "flux")).c_str());
TH1D *eventhist = (TH1D *)inp_file->Get(
(PlotUtils::GetObjectWithName(inp_file, "evt")).c_str());
if (!fluxhist or !eventhist) {
NUIS_ERR(FTL, "Input File Contents: " << inputs[inp_it]);
inp_file->ls();
NUIS_ABORT("NEUT FILE doesn't contain flux/xsec info. You may have to "
"regenerate your MC!");
}
// Get N Events
TTree *neuttree = (TTree *)inp_file->Get("neuttree");
if (!neuttree) {
NUIS_ERR(FTL, "neuttree not located in NEUT file: " << inputs[inp_it]);
NUIS_ABORT(
"Check your inputs, they may need to be completely regenerated!");
throw;
}
int nevents = neuttree->GetEntries();
if (nevents <= 0) {
NUIS_ABORT("Trying to a TTree with "
<< nevents << " to TChain from : " << inputs[inp_it]);
}
// Register input to form flux/event rate hists
RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist);
// Add To TChain
fNEUTTree->AddFile(inputs[inp_it].c_str());
}
// Registor all our file inputs
SetupJointInputs();
// Assign to tree
fEventType = kNEUT;
fNeutVect = NULL;
fNEUTTree->SetBranchAddress("vectorbranch", &fNeutVect);
+ #ifdef ROOT6
+ fNEUTTree->SetAutoDelete(true);
+ #endif
fNEUTTree->GetEntry(0);
// Create Fit Event
fNUISANCEEvent = new FitEvent();
fNUISANCEEvent->SetNeutVect(fNeutVect);
if (fSaveExtra) {
fNeutInfo = new NEUTGeneratorInfo();
fNUISANCEEvent->AddGeneratorInfo(fNeutInfo);
}
fNUISANCEEvent->HardReset();
};
NEUTInputHandler::~NEUTInputHandler(){
// if (fNEUTTree) delete fNEUTTree;
// if (fNeutVect) delete fNeutVect;
// if (fNeutInfo) delete fNeutInfo;
};
void NEUTInputHandler::CreateCache() {
if (fCacheSize > 0) {
// fNEUTTree->SetCacheEntryRange(0, fNEvents);
fNEUTTree->AddBranchToCache("vectorbranch", 1);
fNEUTTree->SetCacheSize(fCacheSize);
}
}
void NEUTInputHandler::RemoveCache() {
// fNEUTTree->SetCacheEntryRange(0, fNEvents);
fNEUTTree->AddBranchToCache("vectorbranch", 0);
fNEUTTree->SetCacheSize(0);
}
FitEvent *NEUTInputHandler::GetNuisanceEvent(const UInt_t ent,
const bool lightweight) {
UInt_t entry = ent + fSkip;
// Catch too large entries
if (entry >= (UInt_t)fNEvents)
return NULL;
// Read Entry from TTree to fill NEUT Vect in BaseFitEvt;
fNEUTTree->GetEntry(entry);
// Run NUISANCE Vector Filler
if (!lightweight) {
CalcNUISANCEKinematics();
}
#ifdef __PROB3PP_ENABLED__
else {
UInt_t npart = fNeutVect->Npart();
for (size_t i = 0; i < npart; i++) {
NeutPart *part = fNUISANCEEvent->fNeutVect->PartInfo(i);
if ((part->fIsAlive == false) && (part->fStatus == -1) &&
std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4,
part->fPID)) {
fNUISANCEEvent->probe_E = part->fP.T();
fNUISANCEEvent->probe_pdg = part->fPID;
break;
} else {
continue;
}
}
}
#endif
// Setup Input scaling for joint inputs
fNUISANCEEvent->InputWeight = GetInputWeight(entry);
// Return event pointer
return fNUISANCEEvent;
}
// From NEUT neutclass/neutpart.h
// Bool_t fIsAlive; // Particle should be tracked or not
// ( in the detector simulator )
//
// Int_t fStatus; // Status flag of this particle
// -2: Non existing particle
// -1: Initial state particle
// 0: Normal
// 1: Decayed to the other particle
// 2: Escaped from the detector
// 3: Absorped
// 4: Charge exchanged
// 5: Pauli blocked
// 6: N/A
// 7: Produced child particles
// 8: Inelastically scattered
//
int NEUTInputHandler::GetNeutParticleStatus(NeutPart *part) {
// State
int state = kUndefinedState;
// Remove Pauli blocked events, probably just single pion events
if (part->fStatus == 5) {
state = kFSIState;
// fStatus == -1 means initial state
} else if (part->fIsAlive == false && part->fStatus == -1) {
state = kInitialState;
// NEUT has a bit of a strange convention for fIsAlive and fStatus
// combinations
// for NC and neutrino particle isAlive true/false and status 2 means
// final state particle
// for other particles in NC status 2 means it's an FSI particle
// for CC it means it was an FSI particle
} else if (part->fStatus == 2) {
// NC case is a little strange... The outgoing neutrino might be alive or
// not alive. Remaining particles with status 2 are FSI particles that
// reinteracted
if (abs(fNeutVect->Mode) > 30 &&
(abs(part->fPID) == 16 || abs(part->fPID) == 14 ||
abs(part->fPID) == 12)) {
state = kFinalState;
// The usual CC case
} else if (part->fIsAlive == true) {
state = kFSIState;
}
} else if (part->fIsAlive == true && part->fStatus == 2 &&
(abs(part->fPID) == 16 || abs(part->fPID) == 14 ||
abs(part->fPID) == 12)) {
state = kFinalState;
} else if (part->fIsAlive == true && part->fStatus == 0) {
state = kFinalState;
} else if (!part->fIsAlive &&
(part->fStatus == 1 || part->fStatus == 3 || part->fStatus == 4 ||
part->fStatus == 7 || part->fStatus == 8)) {
state = kFSIState;
// There's one hyper weird case where fStatus = -3. This apparently
// corresponds to a nucleon being ejected via pion FSI when there is "data
// available"
} else if (!part->fIsAlive && (part->fStatus == -3)) {
state = kUndefinedState;
// NC neutrino outgoing
} else if (!part->fIsAlive && part->fStatus == 0 &&
(abs(part->fPID) == 16 || abs(part->fPID) == 14 ||
abs(part->fPID) == 12)) {
state = kFinalState;
// Warn if we still find alive particles without classifying them
} else if (part->fIsAlive == true) {
NUIS_ABORT("Undefined NEUT state "
<< " Alive: " << part->fIsAlive << " Status: " << part->fStatus
<< " PDG: " << part->fPID);
// Warn if we find dead particles that we haven't classified
} else {
NUIS_ABORT("Undefined NEUT state "
<< " Alive: " << part->fIsAlive << " Status: " << part->fStatus
<< " PDG: " << part->fPID);
}
return state;
}
void NEUTInputHandler::CalcNUISANCEKinematics() {
// Reset all variables
fNUISANCEEvent->ResetEvent();
// Fill Globals
fNUISANCEEvent->Mode = fNeutVect->Mode;
fNUISANCEEvent->fEventNo = fNeutVect->EventNo;
fNUISANCEEvent->fTargetA = fNeutVect->TargetA;
fNUISANCEEvent->fTargetZ = fNeutVect->TargetZ;
fNUISANCEEvent->fTargetH = fNeutVect->TargetH;
fNUISANCEEvent->fBound = bool(fNeutVect->Ibound);
if (fNUISANCEEvent->fBound) {
fNUISANCEEvent->fTargetPDG = TargetUtils::GetTargetPDGFromZA(
fNUISANCEEvent->fTargetZ, fNUISANCEEvent->fTargetA);
} else {
fNUISANCEEvent->fTargetPDG = 1000010010;
}
// Check Particle Stack
UInt_t npart = fNeutVect->Npart();
UInt_t kmax = fNUISANCEEvent->kMaxParticles;
if (npart > kmax) {
NUIS_ERR(WRN, "NEUT has too many particles. Expanding stack.");
fNUISANCEEvent->ExpandParticleStack(npart);
}
int nprimary = fNeutVect->Nprimary();
// Fill Particle Stack
for (size_t i = 0; i < npart; i++) {
// Get Current Count
int curpart = fNUISANCEEvent->fNParticles;
// Get NEUT Particle
NeutPart *part = fNeutVect->PartInfo(i);
// State
int state = GetNeutParticleStatus(part);
// Remove Undefined
if (kRemoveUndefParticles && state == kUndefinedState)
continue;
// Remove FSI
if (kRemoveFSIParticles && state == kFSIState)
continue;
// Remove Nuclear
if (kRemoveNuclearParticles &&
(state == kNuclearInitial || state == kNuclearRemnant))
continue;
// State
fNUISANCEEvent->fParticleState[curpart] = state;
// Is the paricle associated with the primary vertex?
bool primary = false;
// NEUT events are just popped onto the stack as primary, then continues to
// be non-primary
if (i < nprimary)
primary = true;
fNUISANCEEvent->fPrimaryVertex[curpart] = primary;
// Mom
fNUISANCEEvent->fParticleMom[curpart][0] = part->fP.X();
fNUISANCEEvent->fParticleMom[curpart][1] = part->fP.Y();
fNUISANCEEvent->fParticleMom[curpart][2] = part->fP.Z();
fNUISANCEEvent->fParticleMom[curpart][3] = part->fP.T();
// PDG
fNUISANCEEvent->fParticlePDG[curpart] = part->fPID;
// Add up particle count
fNUISANCEEvent->fNParticles++;
}
// Save Extra Generator Info
if (fSaveExtra) {
fNeutInfo->FillGeneratorInfo(fNeutVect);
}
// Run Initial, FSI, Final, Other ordering.
fNUISANCEEvent->OrderStack();
FitParticle *ISAnyLepton = fNUISANCEEvent->GetHMISAnyLeptons();
if (ISAnyLepton) {
fNUISANCEEvent->probe_E = ISAnyLepton->E();
fNUISANCEEvent->probe_pdg = ISAnyLepton->PDG();
}
return;
}
void NEUTUtils::FillNeutCommons(NeutVect *nvect) {
// WARNING: This has only been implemented for a neuttree and not GENIE
// This should be kept in sync with T2KNIWGUtils::GetNIWGEvent(TTree)
// NEUT version info. Can't get it to compile properly with this yet
// neutversion_.corev = nvect->COREVer;
// neutversion_.nucev = nvect->NUCEVer;
// neutversion_.nuccv = nvect->NUCCVer;
// Documentation: See nework.h
nework_.modene = nvect->Mode;
nework_.numne = nvect->Npart();
#ifdef NEUT_COMMON_QEAV
nemdls_.mdlqeaf = nvect->QEAVForm;
#else
nemdls_.mdlqeaf = nvect->QEVForm;
#endif
nemdls_.mdlqe = nvect->QEModel;
nemdls_.mdlspi = nvect->SPIModel;
nemdls_.mdldis = nvect->DISModel;
nemdls_.mdlcoh = nvect->COHModel;
neutcoh_.necohepi = nvect->COHModel;
nemdls_.xmaqe = nvect->QEMA;
nemdls_.xmvqe = nvect->QEMV;
nemdls_.kapp = nvect->KAPPA;
// nemdls_.sccfv = SCCFVdef;
// nemdls_.sccfa = SCCFAdef;
// nemdls_.fpqe = FPQEdef;
nemdls_.xmaspi = nvect->SPIMA;
nemdls_.xmvspi = nvect->SPIMV;
nemdls_.xmares = nvect->RESMA;
nemdls_.xmvres = nvect->RESMV;
neut1pi_.xmanffres = nvect->SPIMA;
neut1pi_.xmvnffres = nvect->SPIMV;
neut1pi_.xmarsres = nvect->RESMA;
neut1pi_.xmvrsres = nvect->RESMV;
neut1pi_.neiff = nvect->SPIForm;
neut1pi_.nenrtype = nvect->SPINRType;
neut1pi_.rneca5i = nvect->SPICA5I;
neut1pi_.rnebgscl = nvect->SPIBGScale;
nemdls_.xmacoh = nvect->COHMA;
nemdls_.rad0nu = nvect->COHR0;
// nemdls_.fa1coh = nvect->COHA1err;
// nemdls_.fb1coh = nvect->COHb1err;
// neutdis_.nepdf = NEPDFdef;
// neutdis_.nebodek = NEBODEKdef;
neutcard_.nefrmflg = nvect->FrmFlg;
neutcard_.nepauflg = nvect->PauFlg;
neutcard_.nenefo16 = nvect->NefO16;
neutcard_.nemodflg = nvect->ModFlg;
// neutcard_.nenefmodl = 1;
// neutcard_.nenefmodh = 1;
// neutcard_.nenefkinh = 1;
// neutpiabs_.neabspiemit = 1;
nenupr_.iformlen = nvect->FormLen;
neutpiless_.ipilessdcy = nvect->IPilessDcy;
neutpiless_.rpilessdcy = nvect->RPilessDcy;
neutpiless_.ipilessdcy = nvect->IPilessDcy;
neutpiless_.rpilessdcy = nvect->RPilessDcy;
neffpr_.fefqe = nvect->NuceffFactorPIQE;
neffpr_.fefqeh = nvect->NuceffFactorPIQEH;
neffpr_.fefinel = nvect->NuceffFactorPIInel;
neffpr_.fefabs = nvect->NuceffFactorPIAbs;
neffpr_.fefcx = nvect->NuceffFactorPICX;
neffpr_.fefcxh = nvect->NuceffFactorPICXH;
neffpr_.fefcoh = nvect->NuceffFactorPICoh;
neffpr_.fefqehf = nvect->NuceffFactorPIQEHKin;
neffpr_.fefcxhf = nvect->NuceffFactorPICXKin;
neffpr_.fefcohf = nvect->NuceffFactorPIQELKin;
for (int i = 0; i < nework_.numne; i++) {
nework_.ipne[i] = nvect->PartInfo(i)->fPID;
nework_.pne[i][0] =
(float)nvect->PartInfo(i)->fP.X() / 1000; // VC(NE)WORK in M(G)eV
nework_.pne[i][1] =
(float)nvect->PartInfo(i)->fP.Y() / 1000; // VC(NE)WORK in M(G)eV
nework_.pne[i][2] =
(float)nvect->PartInfo(i)->fP.Z() / 1000; // VC(NE)WORK in M(G)eV
}
// fsihist.h
// neutroot fills a dummy object for events with no FSI to prevent memory leak
// when
// reading the TTree, so check for it here
if ((int)nvect->NfsiVert() ==
1) { // An event with FSI must have at least two vertices
// if (nvect->NfsiPart()!=1 || nvect->Fsiprob!=-1)
// ERR(WRN) << "T2KNeutUtils::fill_neut_commons(TTree) NfsiPart!=1 or
// Fsiprob!=-1 when NfsiVert==1" << std::endl;
fsihist_.nvert = 0;
fsihist_.nvcvert = 0;
fsihist_.fsiprob = 1;
} else { // Real FSI event
fsihist_.nvert = (int)nvect->NfsiVert();
for (int ivert = 0; ivert < fsihist_.nvert; ivert++) {
fsihist_.iflgvert[ivert] = nvect->FsiVertInfo(ivert)->fVertID;
fsihist_.posvert[ivert][0] = (float)nvect->FsiVertInfo(ivert)->fPos.X();
fsihist_.posvert[ivert][1] = (float)nvect->FsiVertInfo(ivert)->fPos.Y();
fsihist_.posvert[ivert][2] = (float)nvect->FsiVertInfo(ivert)->fPos.Z();
}
fsihist_.nvcvert = nvect->NfsiPart();
for (int ip = 0; ip < fsihist_.nvcvert; ip++) {
fsihist_.abspvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomLab;
fsihist_.abstpvert[ip] = (float)nvect->FsiPartInfo(ip)->fMomNuc;
fsihist_.ipvert[ip] = nvect->FsiPartInfo(ip)->fPID;
fsihist_.iverti[ip] = nvect->FsiPartInfo(ip)->fVertStart;
fsihist_.ivertf[ip] = nvect->FsiPartInfo(ip)->fVertEnd;
fsihist_.dirvert[ip][0] = (float)nvect->FsiPartInfo(ip)->fDir.X();
fsihist_.dirvert[ip][1] = (float)nvect->FsiPartInfo(ip)->fDir.Y();
fsihist_.dirvert[ip][2] = (float)nvect->FsiPartInfo(ip)->fDir.Z();
}
fsihist_.fsiprob = nvect->Fsiprob;
}
neutcrscom_.crsx = nvect->Crsx;
neutcrscom_.crsy = nvect->Crsy;
neutcrscom_.crsz = nvect->Crsz;
neutcrscom_.crsphi = nvect->Crsphi;
neutcrscom_.crsq2 = nvect->Crsq2;
neuttarget_.numbndn = nvect->TargetA - nvect->TargetZ;
neuttarget_.numbndp = nvect->TargetZ;
neuttarget_.numfrep = nvect->TargetH;
neuttarget_.numatom = nvect->TargetA;
posinnuc_.ibound = nvect->Ibound;
// put empty nucleon FSI history (since it is not saved in the NeutVect
// format)
// Comment out as NEUT does not have the necessary proton FSI information yet
// nucleonfsihist_.nfnvert = 0;
// nucleonfsihist_.nfnstep = 0;
}
#endif
diff --git a/src/MCStudies/GenericFlux_Vectors.cxx b/src/MCStudies/GenericFlux_Vectors.cxx
index 13306b7..3eb4410 100644
--- a/src/MCStudies/GenericFlux_Vectors.cxx
+++ b/src/MCStudies/GenericFlux_Vectors.cxx
@@ -1,503 +1,503 @@
// 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 "GenericFlux_Vectors.h"
#ifndef __NO_MINERvA__
#include "MINERvA_SignalDef.h"
#endif
#ifndef __NO_T2K__
#include "T2K_SignalDef.h"
#endif
GenericFlux_Vectors::GenericFlux_Vectors(std::string name,
std::string inputfile, FitWeight *rw,
std::string type,
std::string fakeDataFile) {
// Measurement Details
fName = name;
eventVariables = NULL;
// Define our energy range for flux calcs
EnuMin = 0.;
EnuMax = 1E10; // Arbritrarily high energy limit
if (Config::HasPar("EnuMin")) {
EnuMin = Config::GetParD("EnuMin");
}
if (Config::HasPar("EnuMax")) {
EnuMax = Config::GetParD("EnuMax");
}
SavePreFSI = Config::Get().GetParB("nuisflat_SavePreFSI");
NUIS_LOG(SAM, "Running GenericFlux_Vectors saving pre-FSI particles? "
<< SavePreFSI);
// Set default fitter flags
fIsDiag = true;
fIsShape = false;
fIsRawEvents = false;
// This function will sort out the input files automatically and parse all the
// inputs,flags,etc.
// There may be complex cases where you have to do this by hand, but usually
// this will do.
Measurement1D::SetupMeasurement(inputfile, type, rw, fakeDataFile);
eventVariables = NULL;
// Setup fDataHist as a placeholder
this->fDataHist = new TH1D(("empty_data"), ("empty-data"), 1, 0, 1);
this->SetupDefaultHist();
fFullCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist);
covar = StatUtils::GetInvert(fFullCovar);
// 1. The generator is organised in SetupMeasurement so it gives the
// cross-section in "per nucleon" units.
// So some extra scaling for a specific measurement may be required. For
// Example to get a "per neutron" measurement on carbon
// which we do here, we have to multiple by the number of nucleons 12 and
// divide by the number of neutrons 6.
// N.B. MeasurementBase::PredictedEventRate includes the 1E-38 factor that is
// often included here in other classes that directly integrate the event
// histogram. This method is used here as it now respects EnuMin and EnuMax
// correctly.
this->fScaleFactor =
(this->PredictedEventRate("width", 0, EnuMax) / double(fNEvents)) /
this->TotalIntegratedFlux("width");
NUIS_LOG(SAM, " Generic Flux Scaling Factor = "
<< fScaleFactor << " [= "
<< (GetEventHistogram()->Integral("width") * 1E-38) << "/("
<< (fNEvents + 0.) << "*" << TotalIntegratedFlux("width")
<< ")]");
if (fScaleFactor <= 0.0) {
NUIS_ABORT("SCALE FACTOR TOO LOW");
}
// Setup our TTrees
this->AddEventVariablesToTree();
this->AddSignalFlagsToTree();
}
void GenericFlux_Vectors::AddEventVariablesToTree() {
// Setup the TTree to save everything
if (!eventVariables) {
Config::Get().out->cd();
eventVariables = new TTree((this->fName + "_VARS").c_str(),
(this->fName + "_VARS").c_str());
}
NUIS_LOG(SAM, "Adding Event Variables");
eventVariables->Branch("Mode", &Mode, "Mode/I");
eventVariables->Branch("cc", &cc, "cc/B");
eventVariables->Branch("PDGnu", &PDGnu, "PDGnu/I");
eventVariables->Branch("Enu_true", &Enu_true, "Enu_true/F");
eventVariables->Branch("tgt", &tgt, "tgt/I");
eventVariables->Branch("tgta", &tgta, "tgta/I");
eventVariables->Branch("tgtz", &tgtz, "tgtz/I");
eventVariables->Branch("PDGLep", &PDGLep, "PDGLep/I");
eventVariables->Branch("ELep", &ELep, "ELep/F");
eventVariables->Branch("CosLep", &CosLep, "CosLep/F");
// Basic interaction kinematics
eventVariables->Branch("Q2", &Q2, "Q2/F");
eventVariables->Branch("q0", &q0, "q0/F");
eventVariables->Branch("q3", &q3, "q3/F");
eventVariables->Branch("Enu_QE", &Enu_QE, "Enu_QE/F");
eventVariables->Branch("Q2_QE", &Q2_QE, "Q2_QE/F");
eventVariables->Branch("W_nuc_rest", &W_nuc_rest, "W_nuc_rest/F");
eventVariables->Branch("W", &W, "W/F");
eventVariables->Branch("W_genie", &W_genie, "W_genie/F");
eventVariables->Branch("x", &x, "x/F");
eventVariables->Branch("y", &y, "y/F");
eventVariables->Branch("Eav", &Eav, "Eav/F");
eventVariables->Branch("EavAlt", &EavAlt, "EavAlt/F");
eventVariables->Branch("CosThetaAdler", &CosThetaAdler, "CosThetaAdler/F");
eventVariables->Branch("PhiAdler", &PhiAdler, "PhiAdler/F");
eventVariables->Branch("dalphat", &dalphat, "dalphat/F");
eventVariables->Branch("dpt", &dpt, "dpt/F");
eventVariables->Branch("dphit", &dphit, "dphit/F");
eventVariables->Branch("pnreco_C", &pnreco_C, "pnreco_C/F");
// Save outgoing particle vectors
eventVariables->Branch("nfsp", &nfsp, "nfsp/I");
eventVariables->Branch("px", px, "px[nfsp]/F");
eventVariables->Branch("py", py, "py[nfsp]/F");
eventVariables->Branch("pz", pz, "pz[nfsp]/F");
eventVariables->Branch("E", E, "E[nfsp]/F");
eventVariables->Branch("pdg", pdg, "pdg[nfsp]/I");
eventVariables->Branch("pdg_rank", pdg_rank, "pdg_rank[nfsp]/I");
// Save init particle vectors
eventVariables->Branch("ninitp", &ninitp, "ninitp/I");
eventVariables->Branch("px_init", px_init, "px_init[ninitp]/F");
eventVariables->Branch("py_init", py_init, "py_init[ninitp]/F");
eventVariables->Branch("pz_init", pz_init, "pz_init[ninitp]/F");
eventVariables->Branch("E_init", E_init, "E_init[ninitp]/F");
eventVariables->Branch("pdg_init", pdg_init, "pdg_init[ninitp]/I");
// Save pre-FSI vectors
eventVariables->Branch("nvertp", &nvertp, "nvertp/I");
eventVariables->Branch("px_vert", px_vert, "px_vert[nvertp]/F");
eventVariables->Branch("py_vert", py_vert, "py_vert[nvertp]/F");
eventVariables->Branch("pz_vert", pz_vert, "pz_vert[nvertp]/F");
eventVariables->Branch("E_vert", E_vert, "E_vert[nvertp]/F");
eventVariables->Branch("pdg_vert", pdg_vert, "pdg_vert[nvertp]/I");
// Event Scaling Information
eventVariables->Branch("Weight", &Weight, "Weight/F");
eventVariables->Branch("InputWeight", &InputWeight, "InputWeight/F");
eventVariables->Branch("RWWeight", &RWWeight, "RWWeight/F");
// Should be a double because may be 1E-39 and less
eventVariables->Branch("fScaleFactor", &fScaleFactor, "fScaleFactor/D");
// The customs
eventVariables->Branch("CustomWeight", &CustomWeight, "CustomWeight/F");
eventVariables->Branch("CustomWeightArray", CustomWeightArray,
"CustomWeightArray[6]/F");
return;
}
void GenericFlux_Vectors::FillEventVariables(FitEvent *event) {
ResetVariables();
// Fill Signal Variables
FillSignalFlags(event);
NUIS_LOG(DEB, "Filling signal");
// Now fill the information
Mode = event->Mode;
cc = event->IsCC();
// Get the incoming neutrino and outgoing lepton
FitParticle *nu = event->GetBeamPart();
FitParticle *lep = event->GetHMFSAnyLepton();
PDGnu = nu->fPID;
Enu_true = nu->fP.E() / 1E3;
tgt = event->fTargetPDG;
tgta = event->fTargetA;
tgtz = event->fTargetZ;
TLorentzVector ISP4 = nu->fP;
if (lep != NULL) {
PDGLep = lep->fPID;
ELep = lep->fP.E() / 1E3;
CosLep = cos(nu->fP.Vect().Angle(lep->fP.Vect()));
// Basic interaction kinematics
Q2 = -1 * (nu->fP - lep->fP).Mag2() / 1E6;
q0 = (nu->fP - lep->fP).E() / 1E3;
q3 = (nu->fP - lep->fP).Vect().Mag() / 1E3;
// These assume C12 binding from MINERvA... not ideal
Enu_QE = FitUtils::EnuQErec(lep->fP, CosLep, 34., true);
Q2_QE = FitUtils::Q2QErec(lep->fP, CosLep, 34., true);
Eav = FitUtils::GetErecoil_MINERvA_LowRecoil(event) / 1.E3;
EavAlt = FitUtils::Eavailable(event) / 1.E3;
// Check if this is a 1pi+ or 1pi0 event
if ((SignalDef::isCC1pi(event, PDGnu, 211) ||
SignalDef::isCC1pi(event, PDGnu, -211) ||
SignalDef::isCC1pi(event, PDGnu, 111)) &&
event->NumFSNucleons() == 1) {
TLorentzVector Pnu = nu->fP;
TLorentzVector Pmu = lep->fP;
TLorentzVector Ppi = event->GetHMFSPions()->fP;
TLorentzVector Pprot = event->GetHMFSNucleons()->fP;
CosThetaAdler = FitUtils::CosThAdler(Pnu, Pmu, Ppi, Pprot);
PhiAdler = FitUtils::PhiAdler(Pnu, Pmu, Ppi, Pprot);
}
// Get W_true with assumption of initial state nucleon at rest
float m_n = (float)PhysConst::mass_proton;
// Q2 assuming nucleon at rest
W_nuc_rest = sqrt(-Q2 + 2 * m_n * q0 + m_n * m_n);
// True Q2
x = Q2 / (2 * m_n * q0);
y = 1 - ELep / Enu_true;
dalphat = FitUtils::Get_STV_dalphat_HMProton(event, PDGnu, true);
dpt = FitUtils::Get_STV_dpt_HMProton(event, PDGnu, true);
dphit = FitUtils::Get_STV_dphit_HMProton(event, PDGnu, true);
pnreco_C = FitUtils::Get_pn_reco_C_HMProton(event, PDGnu, true);
}
// Loop over the particles and store all the final state particles in a vector
for (UInt_t i = 0; i < event->Npart(); ++i) {
if (event->PartInfo(i)->fIsAlive &&
event->PartInfo(i)->Status() == kFinalState)
partList.push_back(event->PartInfo(i));
if (SavePreFSI && event->fPrimaryVertex[i])
vertList.push_back(event->PartInfo(i));
if (SavePreFSI && event->PartInfo(i)->IsInitialState())
initList.push_back(event->PartInfo(i));
if (event->PartInfo(i)->IsInitialState()) {
ISP4 += event->PartInfo(i)->fP;
}
}
// Save outgoing particle vectors
nfsp = (int)partList.size();
std::map > > pdgMap;
for (int i = 0; i < nfsp; ++i) {
px[i] = partList[i]->fP.X() / 1E3;
py[i] = partList[i]->fP.Y() / 1E3;
pz[i] = partList[i]->fP.Z() / 1E3;
E[i] = partList[i]->fP.E() / 1E3;
pdg[i] = partList[i]->fPID;
pdgMap[pdg[i]].push_back(std::make_pair(partList[i]->fP.Vect().Mag(), i));
}
for (std::map > >::iterator iter =
pdgMap.begin();
iter != pdgMap.end(); ++iter) {
std::vector > thisVect = iter->second;
std::sort(thisVect.begin(), thisVect.end());
// Now save the order... a bit funky to avoid inverting
int nPart = (int)thisVect.size() - 1;
for (int i = nPart; i >= 0; --i) {
pdg_rank[thisVect[i].second] = nPart - i;
}
}
// Save pre-FSI particles
nvertp = (int)vertList.size();
for (int i = 0; i < nvertp; ++i) {
px_vert[i] = vertList[i]->fP.X() / 1E3;
py_vert[i] = vertList[i]->fP.Y() / 1E3;
pz_vert[i] = vertList[i]->fP.Z() / 1E3;
E_vert[i] = vertList[i]->fP.E() / 1E3;
pdg_vert[i] = vertList[i]->fPID;
}
// Save init particles
ninitp = (int)initList.size();
for (int i = 0; i < ninitp; ++i) {
px_init[i] = initList[i]->fP.X() / 1E3;
py_init[i] = initList[i]->fP.Y() / 1E3;
pz_init[i] = initList[i]->fP.Z() / 1E3;
E_init[i] = initList[i]->fP.E() / 1E3;
pdg_init[i] = initList[i]->fPID;
}
#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();
W_genie = kine.W();
}
#endif
if (lep != NULL) {
W = (ISP4 - lep->fP).M();
} else {
W = 0;
}
// Fill event weights
Weight = event->RWWeight * event->InputWeight;
RWWeight = event->RWWeight;
InputWeight = event->InputWeight;
// And the Customs
CustomWeight = event->CustomWeight;
for (int i = 0; i < 6; ++i) {
CustomWeightArray[i] = event->CustomWeightArray[i];
}
// Fill the eventVariables Tree
eventVariables->Fill();
return;
};
//********************************************************************
void GenericFlux_Vectors::ResetVariables() {
//********************************************************************
cc = false;
// Reset all Function used to extract any variables of interest to the event
Mode = PDGnu = tgt = tgta = tgtz = PDGLep = 0;
Enu_true = ELep = CosLep = Q2 = q0 = q3 = Enu_QE = Q2_QE = W_nuc_rest = W =
x = y = Eav = EavAlt = CosThetaAdler = PhiAdler = -999.9;
W_genie = -999;
// Other fun variables
// MINERvA-like ones
dalphat = dpt = dphit = pnreco_C = -999.99;
nfsp = ninitp = nvertp = 0;
for (int i = 0; i < kMAX; ++i) {
px[i] = py[i] = pz[i] = E[i] = -999;
pdg[i] = pdg_rank[i] = 0;
px_init[i] = py_init[i] = pz_init[i] = E_init[i] = -999;
pdg_init[i] = 0;
px_vert[i] = py_vert[i] = pz_vert[i] = E_vert[i] = -999;
pdg_vert[i] = 0;
}
Weight = InputWeight = RWWeight = 0.0;
CustomWeight = 0.0;
for (int i = 0; i < 6; ++i)
CustomWeightArray[i] = 0.0;
partList.clear();
initList.clear();
vertList.clear();
flagCCINC = flagNCINC = flagCCQE = flagCC0pi = flagCCQELike = flagNCEL =
flagNC0pi = flagCCcoh = flagNCcoh = flagCC1pip = flagNC1pip = flagCC1pim =
flagNC1pim = flagCC1pi0 = flagNC1pi0 = false;
#ifndef __NO_MINERvA__
flagCC0piMINERvA = false;
#endif
#ifndef __NO_T2K__
flagCC0Pi_T2K_AnaI = false;
flagCC0Pi_T2K_AnaII = false;
#endif
}
//********************************************************************
void GenericFlux_Vectors::FillSignalFlags(FitEvent *event) {
//********************************************************************
// Some example flags are given from SignalDef.
// See src/Utils/SignalDef.cxx for more.
int nuPDG = event->PartInfo(0)->fPID;
// Generic signal flags
flagCCINC = SignalDef::isCCINC(event, nuPDG);
flagNCINC = SignalDef::isNCINC(event, nuPDG);
flagCCQE = SignalDef::isCCQE(event, nuPDG);
flagCCQELike = SignalDef::isCCQELike(event, nuPDG);
flagCC0pi = SignalDef::isCC0pi(event, nuPDG);
flagNCEL = SignalDef::isNCEL(event, nuPDG);
flagNC0pi = SignalDef::isNC0pi(event, nuPDG);
flagCCcoh = SignalDef::isCCCOH(event, nuPDG, 211);
flagNCcoh = SignalDef::isNCCOH(event, nuPDG, 111);
flagCC1pip = SignalDef::isCC1pi(event, nuPDG, 211);
flagNC1pip = SignalDef::isNC1pi(event, nuPDG, 211);
flagCC1pim = SignalDef::isCC1pi(event, nuPDG, -211);
flagNC1pim = SignalDef::isNC1pi(event, nuPDG, -211);
flagCC1pi0 = SignalDef::isCC1pi(event, nuPDG, 111);
flagNC1pi0 = SignalDef::isNC1pi(event, nuPDG, 111);
#ifndef __NO_MINERvA__
flagCC0piMINERvA = SignalDef::isCC0pi_MINERvAPTPZ(event, 14);
#endif
#ifndef __NO_T2K__
flagCC0Pi_T2K_AnaI =
SignalDef::isT2K_CC0pi(event, EnuMin, EnuMax, SignalDef::kAnalysis_I);
flagCC0Pi_T2K_AnaII =
SignalDef::isT2K_CC0pi(event, EnuMin, EnuMax, SignalDef::kAnalysis_II);
#endif
}
void GenericFlux_Vectors::AddSignalFlagsToTree() {
if (!eventVariables) {
Config::Get().out->cd();
eventVariables = new TTree((this->fName + "_VARS").c_str(),
(this->fName + "_VARS").c_str());
}
NUIS_LOG(SAM, "Adding signal flags");
// Signal Definitions from SignalDef.cxx
eventVariables->Branch("flagCCINC", &flagCCINC, "flagCCINC/O");
eventVariables->Branch("flagNCINC", &flagNCINC, "flagNCINC/O");
eventVariables->Branch("flagCCQE", &flagCCQE, "flagCCQE/O");
eventVariables->Branch("flagCC0pi", &flagCC0pi, "flagCC0pi/O");
eventVariables->Branch("flagCCQELike", &flagCCQELike, "flagCCQELike/O");
eventVariables->Branch("flagNCEL", &flagNCEL, "flagNCEL/O");
eventVariables->Branch("flagNC0pi", &flagNC0pi, "flagNC0pi/O");
eventVariables->Branch("flagCCcoh", &flagCCcoh, "flagCCcoh/O");
eventVariables->Branch("flagNCcoh", &flagNCcoh, "flagNCcoh/O");
eventVariables->Branch("flagCC1pip", &flagCC1pip, "flagCC1pip/O");
eventVariables->Branch("flagNC1pip", &flagNC1pip, "flagNC1pip/O");
eventVariables->Branch("flagCC1pim", &flagCC1pim, "flagCC1pim/O");
eventVariables->Branch("flagNC1pim", &flagNC1pim, "flagNC1pim/O");
eventVariables->Branch("flagCC1pi0", &flagCC1pi0, "flagCC1pi0/O");
eventVariables->Branch("flagNC1pi0", &flagNC1pi0, "flagNC1pi0/O");
#ifndef __NO_MINERvA__
eventVariables->Branch("flagCC0piMINERvA", &flagCC0piMINERvA,
"flagCC0piMINERvA/O");
#endif
#ifndef __NO_T2K__
eventVariables->Branch("flagCC0Pi_T2K_AnaI", &flagCC0Pi_T2K_AnaI,
"flagCC0Pi_T2K_AnaI/O");
eventVariables->Branch("flagCC0Pi_T2K_AnaII", &flagCC0Pi_T2K_AnaII,
"flagCC0Pi_T2K_AnaII/O");
#endif
};
void GenericFlux_Vectors::Write(std::string drawOpt) {
// First save the TTree
eventVariables->Write();
// Save Flux and Event Histograms too
GetInput()->GetFluxHistogram()->Write();
GetInput()->GetEventHistogram()->Write();
return;
}
// Override functions which aren't really necessary
bool GenericFlux_Vectors::isSignal(FitEvent *event) {
(void)event;
return true;
};
void GenericFlux_Vectors::ScaleEvents() { return; }
void GenericFlux_Vectors::ApplyNormScale(float norm) {
this->fCurrentNorm = norm;
return;
}
void GenericFlux_Vectors::FillHistograms() { return; }
void GenericFlux_Vectors::ResetAll() {
- eventVariables->Reset();
+ // eventVariables->Reset();
return;
}
float GenericFlux_Vectors::GetChi2() { return 0.0; }
diff --git a/src/Reweight/NOvARwgtEngine.cxx b/src/Reweight/NOvARwgtEngine.cxx
index 1277db9..7a60320 100644
--- a/src/Reweight/NOvARwgtEngine.cxx
+++ b/src/Reweight/NOvARwgtEngine.cxx
@@ -1,232 +1,317 @@
#include "NOvARwgtEngine.h"
#include "NOvARwgt/interfaces/GenieInterface.h"
#include "NOvARwgt/rwgt/tunes/Tunes2017.h"
#include "NOvARwgt/rwgt/tunes/Tunes2018.h"
-#include "NOvARwgt/rwgt/tunes/Tunes2019.h"
+#include "NOvARwgt/rwgt/tunes/Tunes2020.h"
#include "NOvARwgt/rwgt/tunes/TunesSA.h"
#include "NOvARwgt/rwgt/EventRecord.h"
-#include "NOvARwgt/rwgt/Tune.h"
#include "NOvARwgt/rwgt/ISystKnob.h"
+#include "NOvARwgt/rwgt/Tune.h"
+
+// GENIEv3
+#include "Framework/Utils/AppInit.h"
+#include "Framework/Utils/RunOpt.h"
static size_t const kCVTune2017 = 100;
static size_t const kCVTune2018 = 200;
-static size_t const kCVTune2019 = 300;
+static size_t const kCVTune2020 = 300;
static size_t const kCVTuneSA = 400;
static size_t const kNoSuchKnob = std::numeric_limits::max();
+novarwgt::Tune const *TuneFactory(size_t e) {
+ switch (e) {
+ case kCVTune2017: {
+ return &novarwgt::kCVTune2017;
+ }
+ case kCVTune2018: {
+ return &novarwgt::kCVTune2018;
+ }
+ case kCVTune2020: {
+ return &novarwgt::kCVTune2020;
+ }
+ case kCVTuneSA: {
+ return &novarwgt::kCVTuneSA;
+ }
+ default: {
+ return NULL;
+ }
+ }
+}
+
+void NOvARwgtEngine::InitializeKnobs() {
+
+ fTunes.push_back(TuneFactory(kCVTune2017));
+ fTuneEnums[kCVTune2017] = 0;
+ fTunes.push_back(TuneFactory(kCVTune2018));
+ fTuneEnums[kCVTune2018] = 1;
+ fTunes.push_back(TuneFactory(kCVTune2020));
+ fTuneEnums[kCVTune2020] = 2;
+ fTunes.push_back(TuneFactory(kCVTuneSA));
+ fTuneEnums[kCVTuneSA] = 3;
+ fTuneInUse = {false, false, false, false};
+ fTuneValues = {0, 0, 0, 0};
+
+ size_t ctr = 0;
+
+ NUIS_LOG(FIT, "NOvARwgt kCVTune2017 sub-knobs:");
+ size_t tune_ctr = 0;
+ for (auto &k : novarwgt::kCVTune2017.SystKnobs()) {
+ NUIS_LOG(FIT, "\t" << k.first);
+ fKnobs.push_back(k.second);
+ fKnobTunes.push_back(&novarwgt::kCVTune2017);
+ fKnobTuneidx.push_back(0);
+ fKnobEnums[kCVTune2017 + 1 + tune_ctr++] = ctr++;
+ fKnobInUse.push_back(false);
+ fKnobValues.push_back(0);
+ }
+
+ tune_ctr = 0;
+ NUIS_LOG(FIT, "NOvARwgt kCVTune2018 sub-knobs:");
+ for (auto &k : novarwgt::kCVTune2018.SystKnobs()) {
+ NUIS_LOG(FIT, "\t" << k.first);
+ fKnobs.push_back(k.second);
+ fKnobTunes.push_back(&novarwgt::kCVTune2018);
+ fKnobTuneidx.push_back(1);
+ fKnobEnums[kCVTune2018 + 1 + tune_ctr++] = ctr++;
+ fKnobInUse.push_back(false);
+ fKnobValues.push_back(0);
+ }
+
+ tune_ctr = 0;
+ NUIS_LOG(FIT, "NOvARwgt kCVTune2020 sub-knobs:");
+ for (auto &k : novarwgt::kCVTune2020.SystKnobs()) {
+ NUIS_LOG(FIT, "\t" << k.first);
+ fKnobs.push_back(k.second);
+ fKnobTunes.push_back(&novarwgt::kCVTune2020);
+ fKnobTuneidx.push_back(2);
+ fKnobEnums[kCVTune2020 + 1 + tune_ctr++] = ctr++;
+ fKnobInUse.push_back(false);
+ fKnobValues.push_back(0);
+ }
+
+ tune_ctr = 0;
+ NUIS_LOG(FIT, "NOvARwgt kCVTuneSA sub-knobs:");
+ for (auto &k : novarwgt::kCVTuneSA.SystKnobs()) {
+ NUIS_LOG(FIT, "\t" << k.first);
+ fKnobs.push_back(k.second);
+ fKnobTunes.push_back(&novarwgt::kCVTuneSA);
+ fKnobTuneidx.push_back(3);
+ fKnobEnums[kCVTuneSA + 1 + tune_ctr++] = ctr++;
+ fKnobInUse.push_back(false);
+ fKnobValues.push_back(0);
+ }
+}
+
+void NOvARwgtEngine::InitializeGENIE() {
+ genie::RunOpt *grunopt = genie::RunOpt::Instance();
+ grunopt->EnableBareXSecPreCalc(true);
+ grunopt->SetEventGeneratorList(Config::GetParS("GENIEEventGeneratorList"));
+ if (!Config::HasPar("GENIETune")) {
+ NUIS_ABORT(
+ "GENIE tune was not specified, this is required when reweighting GENIE "
+ "V3+ events. Add a config parameter like: to your nuisance card.");
+ }
+ grunopt->SetTuneName(Config::GetParS("GENIETune"));
+ grunopt->BuildTune();
+ std::string genv =
+ std::string(getenv("GENIE")) + "/config/Messenger_laconic.xml";
+ genie::utils::app_init::MesgThresholds(genv);
+}
+
size_t NOvARwgtEngine::GetWeightGeneratorIndex(std::string const &strname) {
- int upos = strname.find_first_of("_");
+ size_t upos = strname.find_first_of("_");
if (strname.find("CVTune2017") == 0) {
if (upos == std::string::npos) {
return kCVTune2017;
}
std::string knobname = strname.substr(upos + 1);
- if (novarwgt::kCVTune2017.SystKnobs().contains(knobname)) {
- auto loc = std::find(fKnobs.begin(), fKnobs.end(),
- novarwgt::kCVTune2017.SystKnobs()[knobname]);
- if (loc == fKnobs.end()) {
- size_t rtn = kCVTune2017 + 1 + fKnobs.size();
- fKnobs.push_back(novarwgt::kCVTune2017.SystKnobs()[knobname]);
- return rtn;
- } else {
- return kCVTune2017 + 1 + std::distance(fKnobs.begin(), loc);
- }
+ if (novarwgt::kCVTune2017.SystKnobs().count(knobname)) {
+ auto loc = novarwgt::kCVTune2017.SystKnobs().find(knobname);
+
+ return kCVTune2017 + 1 +
+ std::distance(novarwgt::kCVTune2017.SystKnobs().begin(), loc);
}
+
} else if (strname.find("CVTune2018") == 0) {
if (upos == std::string::npos) {
return kCVTune2018;
}
std::string knobname = strname.substr(upos + 1);
- if (novarwgt::kCVTune2018.SystKnobs().contains(knobname)) {
- auto loc = std::find(fKnobs.begin(), fKnobs.end(),
- novarwgt::kCVTune2018.SystKnobs()[knobname]);
- if (loc == fKnobs.end()) {
- size_t rtn = kCVTune2018 + 1 + fKnobs.size();
- fKnobs.push_back(novarwgt::kCVTune2018.SystKnobs()[knobname]);
- return rtn;
- } else {
- return kCVTune2018 + 1 + std::distance(fKnobs.begin(), loc);
- }
+ if (novarwgt::kCVTune2018.SystKnobs().count(knobname)) {
+ auto loc = novarwgt::kCVTune2018.SystKnobs().find(knobname);
+
+ return kCVTune2018 + 1 +
+ std::distance(novarwgt::kCVTune2018.SystKnobs().begin(), loc);
}
- } else if (strname.find("CVTune2019") == 0) {
+ } else if (strname.find("CVTune2020") == 0) {
if (upos == std::string::npos) {
- return kCVTune2019;
+ return kCVTune2020;
}
std::string knobname = strname.substr(upos + 1);
- if (novarwgt::kCVTune2019.SystKnobs().contains(knobname)) {
- auto loc = std::find(fKnobs.begin(), fKnobs.end(),
- novarwgt::kCVTune2019.SystKnobs()[knobname]);
- if (loc == fKnobs.end()) {
- size_t rtn = kCVTune2019 + 1 + fKnobs.size();
- fKnobs.push_back(novarwgt::kCVTune2019.SystKnobs()[knobname]);
- return rtn;
- } else {
- return kCVTune2019 + 1 + std::distance(fKnobs.begin(), loc);
- }
+ if (novarwgt::kCVTune2020.SystKnobs().count(knobname)) {
+ auto loc = novarwgt::kCVTune2020.SystKnobs().find(knobname);
+
+ return kCVTune2020 + 1 +
+ std::distance(novarwgt::kCVTune2020.SystKnobs().begin(), loc);
}
} else if (strname.find("CVTuneSA") == 0) {
if (upos == std::string::npos) {
return kCVTuneSA;
}
std::string knobname = strname.substr(upos + 1);
- if (novarwgt::kCVTuneSA.SystKnobs().contains(knobname)) {
- auto loc = std::find(fKnobs.begin(), fKnobs.end(),
- novarwgt::kCVTuneSA.SystKnobs()[knobname]);
- if (loc == fKnobs.end()) {
- size_t rtn = kCVTuneSA + 1 + fKnobs.size();
- fKnobs.push_back(novarwgt::kCVTuneSA.SystKnobs()[knobname]);
- return rtn;
- } else {
- return kCVTuneSA + 1 + std::distance(fKnobs.begin(), loc);
- }
+ if (novarwgt::kCVTuneSA.SystKnobs().count(knobname)) {
+ auto loc = novarwgt::kCVTuneSA.SystKnobs().find(knobname);
+
+ return kCVTuneSA + 1 +
+ std::distance(novarwgt::kCVTuneSA.SystKnobs().begin(), loc);
}
}
return kNoSuchKnob;
}
-novarwgt::Tune const *TuneFactory(size_t e) {
- switch (e) {
- case kTune2017: {
- return &novarwgt::kCVTune2017;
- }
- case kTune2018: {
- return &novarwgt::kCVTune2018;
- }
- case kTune2019: {
- return &novarwgt::kCVTune2019;
- }
- case kTuneSA: {
- return &novarwgt::kCVTuneSA;
- }
- default: { return NULL; }
- }
-}
-
void NOvARwgtEngine::IncludeDial(std::string name, double startval) {
- size_t we_e = GetWeightGeneratorIndex(name);
- if (we_e == kNoSuchKnob) {
+ size_t we_indx = GetWeightGeneratorIndex(name);
+ if (we_indx == kNoSuchKnob) {
NUIS_ABORT("[ERROR]: Invalid NOvARwgt Engine name: " << name);
}
- bool IsTune = !(we_e % 100);
+ bool IsTune = !(we_indx % 100);
if (IsTune) {
- if (fTuneEnums.find(we_e) != fTuneEnums.end()) {
- NUIS_ABORT("[ERROR]: NOvARwgt Tune name: " << name
- << " already included.");
- }
- fTuneEnums[we_e] = fTunes.size();
- fTunes.push_back(TuneFactory(we_e));
- fTuneValues.push_back(startval);
+ auto tune_idx = fTuneEnums[we_indx];
+ fTuneValues[tune_idx] = startval;
+ fTuneInUse[tune_idx] = true;
} else {
- // The GetWeightGeneratorIndex caches the Knob so we don't have to.
- fKnobEnums[we_e] = fKnobValues.size();
- fKnobValues.push_back(startval);
+ auto knob_idx = fKnobEnums[we_indx];
+ fKnobValues[knob_idx] = startval;
+ fKnobInUse[knob_idx] = true;
}
};
void NOvARwgtEngine::SetDialValue(int nuisenum, double val) {
size_t we_indx = (nuisenum % 1000);
- bool IsTune = !(we_e % 100);
+ bool IsTune = !(we_indx % 100);
if (IsTune) {
- if (!fTuneEnums.count(we_indx)) {
+ auto tune_idx = fTuneEnums[we_indx];
+ if (!fTuneInUse[tune_idx]) {
NUIS_ABORT("[ERROR]: SetDialValue for NOvARwgt dial: "
<< we_indx << " but that tune hasn't been included yet.");
}
- size_t engine_index = fTuneEnums[we_indx];
- fTuneValues[engine_index] = val;
+ fTuneValues[tune_idx] = val;
} else {
- if (!fKnobEnums.count(we_indx)) {
+ auto knob_idx = fKnobEnums[we_indx];
+ if (!fKnobInUse[knob_idx]) {
NUIS_ABORT("[ERROR]: SetDialValue for NOvARwgt dial: "
<< we_indx << " but that Knob hasn't been included yet.");
}
- size_t engine_index = fKnobEnums[we_indx];
- fKnobValues[engine_index] = val;
+ fKnobValues[knob_idx] = val;
}
}
void NOvARwgtEngine::SetDialValue(std::string name, double val) {
SetDialValue(GetWeightGeneratorIndex(name), val);
}
bool NOvARwgtEngine::IsDialIncluded(std::string name) {
return IsDialIncluded(GetWeightGeneratorIndex(name));
}
bool NOvARwgtEngine::IsDialIncluded(int nuisenum) {
size_t we_indx = (nuisenum % 1000);
- bool IsTune = !(we_e % 100);
+ bool IsTune = !(we_indx % 100);
if (IsTune) {
- return fTuneEnums.count(we_indx);
+ auto tune_idx = fTuneEnums[we_indx];
+ return fTuneInUse[tune_idx];
} else {
- return fKnobEnums.count(we_indx);
+ auto knob_idx = fKnobEnums[we_indx];
+ return fKnobInUse[knob_idx];
}
}
double NOvARwgtEngine::GetDialValue(std::string name) {
- return GetDialValueGetWeightGeneratorIndex(name));
+ return GetDialValue(GetWeightGeneratorIndex(name));
}
double NOvARwgtEngine::GetDialValue(int nuisenum) {
size_t we_indx = (nuisenum % 1000);
if (we_indx == kNoSuchKnob) {
- NUIS_ABORT("[ERROR]: Invalid NOvARwgt Engine name: " << name);
+ NUIS_ABORT("[ERROR]: Invalid NOvARwgt Engine enum: " << nuisenum);
}
- bool IsTune = !(we_e % 100);
+ bool IsTune = !(we_indx % 100);
if (IsTune) {
- if (!fTuneEnums.count(we_indx)) {
- NUIS_ABORT("[ERROR]: SetDialValue for NOvARwgt dial: "
+ auto tune_idx = fTuneEnums[we_indx];
+ if (!fTuneInUse[tune_idx]) {
+ NUIS_ABORT("[ERROR]: GetDialValue for NOvARwgt dial: "
<< we_indx << " but that tune hasn't been included yet.");
}
- return fTuneValues[fTuneEnums[we_indx]];
+ return fTuneValues[tune_idx];
} else {
-
- if (!fKnobEnums.count(we_indx)) {
- NUIS_ABORT("[ERROR]: SetDialValue for NOvARwgt dial: "
- << we_indx << " but that engine hasn't been included yet.");
+ auto knob_idx = fKnobEnums[we_indx];
+ if (!fKnobInUse[knob_idx]) {
+ NUIS_ABORT("[ERROR]: GetDialValue for NOvARwgt dial: "
+ << we_indx << " but that Knob hasn't been included yet.");
}
- return fKnobValues[fKnobEnums[we_indx]];
+ return fKnobValues[knob_idx];
}
}
double NOvARwgtEngine::CalcWeight(BaseFitEvt *evt) {
double rw_weight = 1.0;
// Make nom weight
if (!evt) {
NUIS_ABORT("evt not found : " << evt);
}
// Skip Non GENIE
if (evt->fType != kGENIE)
return 1.0;
if (!(evt->genie_event)) {
NUIS_ABORT("evt->genie_event not found!" << evt->genie_event);
}
if (!(evt->genie_event->event)) {
NUIS_ABORT("evt->genie_event->event GHepRecord not found!"
<< (evt->genie_event->event));
}
novarwgt::EventRecord rcd =
novarwgt::ConvertGenieEvent(evt->genie_event->event);
- for (size_t w_it = 0; w_it < fKnobs.size(); ++w_it) {
- rw_weight *= fKnobs[w_it]->CalcWeight(fKnobValues[w_it], rcd);
+ for (size_t k_it = 0; k_it < fKnobs.size(); ++k_it) {
+ if (!fKnobInUse[k_it]) {
+ continue;
+ }
+
+ double wght = fKnobTunes[k_it]->EventSystKnobWeight(
+ fKnobs[k_it]->Name(), fKnobValues[k_it], rcd, {},
+ fTuneInUse[fKnobTuneidx[k_it]] && fTuneValues[fKnobTuneidx[k_it]]);
+
+ // // have to divide out the CV weight for this, ugly hack because the last
+ // // parameter doesn't do what I want
+ // if (fTuneInUse[fKnobTuneidx[k_it]] && fTuneValues[fKnobTuneidx[k_it]]) {
+ // wght /= fKnobTunes[k_it]->EventSystKnobWeight(fKnobs[k_it]->Name(), 0,
+ // rcd, {}, false);
+ // }
+
+ rw_weight *= wght;
}
- for (size_t w_it = 0; w_it < fTunes.size(); ++w_it) {
- if (fTuneValues[w_it] ==
- 0) { // if a dial is set to 0, don't include its weight
+ for (size_t k_it = 0; k_it < fTunes.size(); ++k_it) {
+ if (!fTuneInUse[k_it]) {
+ continue;
+ }
+ if (!fTuneValues[k_it]) {
continue;
}
- rw_weight *= fTunes[w_it]->EventWeight(rcd);
+ double wght = fTunes[k_it]->EventWeight(rcd);
+ rw_weight *= wght;
}
return rw_weight;
}
diff --git a/src/Reweight/NOvARwgtEngine.h b/src/Reweight/NOvARwgtEngine.h
index 33d96d1..9f3dbeb 100644
--- a/src/Reweight/NOvARwgtEngine.h
+++ b/src/Reweight/NOvARwgtEngine.h
@@ -1,41 +1,50 @@
#pragma once
#include "WeightEngineBase.h"
namespace novarwgt {
class ISystKnob;
class Tune;
-}
+} // namespace novarwgt
class NOvARwgtEngine : public WeightEngineBase {
public:
- NOvARwgtEngine(){};
+ NOvARwgtEngine() {
+ InitializeKnobs();
+ InitializeGENIE();
+ };
virtual ~NOvARwgtEngine(){};
+ void InitializeKnobs();
+ void InitializeGENIE();
static size_t GetWeightGeneratorIndex(std::string const &);
// Functions requiring Override
void IncludeDial(std::string name, double startval);
void SetDialValue(int nuisenum, double val);
void SetDialValue(std::string name, double val);
bool IsDialIncluded(std::string name);
bool IsDialIncluded(int nuisenum);
double GetDialValue(std::string name);
double GetDialValue(int nuisenum);
void Reconfigure(bool silent){};
double CalcWeight(BaseFitEvt *evt);
bool NeedsEventReWeight() { return true; }
std::map fTuneEnums;
std::vector fTunes;
std::vector fTuneValues;
+ std::vector fTuneInUse;
std::map fKnobEnums;
std::vector fKnobs;
+ std::vector fKnobTunes;
+ std::vector fKnobTuneidx;
std::vector fKnobValues;
+ std::vector fKnobInUse;
};
diff --git a/src/Reweight/T2KWeightEngine.cxx b/src/Reweight/T2KWeightEngine.cxx
index 134700a..c4725e9 100644
--- a/src/Reweight/T2KWeightEngine.cxx
+++ b/src/Reweight/T2KWeightEngine.cxx
@@ -1,137 +1,148 @@
#include "T2KWeightEngine.h"
+#ifdef __T2KREW_ENABLED__
+#include "T2KNeutUtils.h"
+#endif
T2KWeightEngine::T2KWeightEngine(std::string name) {
#ifdef __T2KREW_ENABLED__
+ std::string neut_card = FitPar::Config().GetParS("NEUT_CARD");
+ if (!neut_card.size()) {
+ NUIS_ABORT(
+ "[ERROR]: When using T2KReWeight must set NEUT_CARD config option.");
+ }
+
+ t2krew::T2KNeutUtils::SetCardFile(neut_card);
+
// Setup the NEUT Reweight engien
fCalcName = name;
NUIS_LOG(FIT, "Setting up T2K RW : " << fCalcName);
// Create RW Engine suppressing cout
StopTalking();
// Create Main RW Engine
fT2KRW = new t2krew::T2KReWeight();
// Setup Sub RW Engines (Only activated for neut and niwg)
fT2KNeutRW = new t2krew::T2KNeutReWeight();
fT2KNIWGRW = new t2krew::T2KNIWGReWeight();
fT2KRW->AdoptWghtEngine("fNeutRW", fT2KNeutRW);
fT2KRW->AdoptWghtEngine("fNIWGRW", fT2KNIWGRW);
fT2KRW->Reconfigure();
// allow cout again
StartTalking();
// Set Abs Twk Config
fIsAbsTwk = (FitPar::Config().GetParB("setabstwk"));
#else
NUIS_ABORT("T2K RW NOT ENABLED");
#endif
};
void T2KWeightEngine::IncludeDial(std::string name, double startval) {
#ifdef __T2KREW_ENABLED__
// Get First enum
int nuisenum = Reweight::ConvDial(name, kT2K);
// Setup Maps
fEnumIndex[nuisenum]; // = std::vector(0);
fNameIndex[name]; // = std::vector(0);
// Split by commas
std::vector allnames = GeneralUtils::ParseToStr(name, ",");
for (uint i = 0; i < allnames.size(); i++) {
std::string singlename = allnames[i];
// Get RW
t2krew::T2KSyst_t gensyst = t2krew::T2KSyst::FromString(name);
// Fill Maps
int index = fValues.size();
fValues.push_back(0.0);
fT2KSysts.push_back(gensyst);
// Initialize dial
std::cout << "Registering " << singlename << " from " << name << std::endl;
fT2KRW->Systematics().Include(gensyst);
// If Absolute
if (fIsAbsTwk) {
fT2KRW->Systematics().SetAbsTwk(gensyst);
}
// Setup index
fEnumIndex[nuisenum].push_back(index);
fNameIndex[name].push_back(index);
}
// Set Value if given
if (startval != _UNDEF_DIAL_VALUE_) {
SetDialValue(nuisenum, startval);
}
#endif
}
void T2KWeightEngine::SetDialValue(int nuisenum, double val) {
#ifdef __T2KREW_ENABLED__
std::vector indices = fEnumIndex[nuisenum];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
fT2KRW->Systematics().SetTwkDial(fT2KSysts[indices[i]], val);
}
#endif
}
void T2KWeightEngine::SetDialValue(std::string name, double val) {
#ifdef __T2KREW_ENABLED__
std::vector indices = fNameIndex[name];
for (uint i = 0; i < indices.size(); i++) {
fValues[indices[i]] = val;
fT2KRW->Systematics().SetTwkDial(fT2KSysts[indices[i]], val);
}
#endif
}
void T2KWeightEngine::Reconfigure(bool silent) {
#ifdef __T2KREW_ENABLED__
// Hush now...
if (silent)
StopTalking();
// Reconf
StopTalking();
fT2KRW->Reconfigure();
StartTalking();
// Shout again
if (silent)
StartTalking();
#endif
}
double T2KWeightEngine::CalcWeight(BaseFitEvt *evt) {
double rw_weight = 1.0;
#ifdef __T2KREW_ENABLED__
// Skip Non GENIE
if (evt->fType != kNEUT)
return 1.0;
// Hush now
StopTalking();
// Get Weight For NEUT
rw_weight = fT2KRW->CalcWeight(evt->fNeutVect);
// Speak Now
StartTalking();
#endif
// Return rw_weight
return rw_weight;
}
diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx
index 0875981..be6e097 100644
--- a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx
+++ b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.cxx
@@ -1,330 +1,373 @@
// 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 "T2K_SignalDef.h"
#include "T2K_CC0pi_XSec_2DPcos_nu_I.h"
static size_t nangbins = 9;
static double angular_binning_costheta[] = {-1, 0, 0.6, 0.7, 0.8,
0.85, 0.9, 0.94, 0.98, 1};
//********************************************************************
T2K_CC0pi_XSec_2DPcos_nu_I::T2K_CC0pi_XSec_2DPcos_nu_I(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "T2K_CC0pi_XSec_2DPcos_nu_I sample. \n"
"Target: CH \n"
"Flux: T2K 2.5 degree off-axis (ND280) \n"
"Signal: CC0pi\n"
"https://journals.aps.org/prd/abstract/10.1103/"
"PhysRevD.93.112012 Analysis I";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("P_{#mu} (GeV)");
fSettings.SetYTitle("cos#theta_{#mu}");
fSettings.SetZTitle("d^{2}#sigma/dP_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)");
fSettings.SetAllowedTypes("FULL,DIAG/FREE,SHAPE,FIX/SYSTCOV/STATCOV",
"FIX/FULL");
fSettings.DefineAllowedTargets("C,H");
fSettings.SetEnuRangeFromFlux(fFluxHist);
// CCQELike plot information
fSettings.SetTitle("T2K_CC0pi_XSec_2DPcos_nu_I");
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
fMaskMomOverflow = false;
if (samplekey.Has("mask_mom_overflow")) {
fMaskMomOverflow = samplekey.GetB("mask_mom_overflow");
}
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = ((GetEventHistogram()->Integral("width") / (fNEvents + 0.)) *
1E-38 / (TotalIntegratedFlux()));
assert(std::isnormal(fScaleFactor));
// Setup Histograms
SetHistograms();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
bool T2K_CC0pi_XSec_2DPcos_nu_I::isSignal(FitEvent *event) {
return SignalDef::isT2K_CC0pi(event, EnuMin, EnuMax, SignalDef::kAnalysis_I);
};
void T2K_CC0pi_XSec_2DPcos_nu_I::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double pmu = Pmu.Vect().Mag() / 1000.;
double CosThetaMu = cos(Pnu.Vect().Angle(Pmu.Vect()));
fXVar = pmu;
fYVar = CosThetaMu;
+ fMode = event->Mode;
+
return;
};
void T2K_CC0pi_XSec_2DPcos_nu_I::FillHistograms() {
Measurement1D::FillHistograms();
if (Signal) {
fMCHist_Fine2D->Fill(fXVar, fYVar, Weight);
FillMCSlice(fXVar, fYVar, Weight);
}
}
// Modification is needed after the full reconfigure to move bins around
// Otherwise this would need to be replaced by a TH2Poly which is too awkward.
void T2K_CC0pi_XSec_2DPcos_nu_I::ConvertEventRates() {
// Do standard conversion.
Measurement1D::ConvertEventRates();
// Scale MC slices by their bin area
for (size_t i = 0; i < nangbins; ++i) {
fMCHist_Slices[i]->Scale(fScaleFactor / (angular_binning_costheta[i + 1] -
angular_binning_costheta[i]),
"width");
}
// Now Convert into 1D list
fMCHist->Reset();
int bincount = 0;
for (size_t i = 0; i < nangbins; i++) {
for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) {
fMCHist->SetBinContent(bincount + 1,
fMCHist_Slices[i]->GetBinContent(j + 1));
fMCHist->SetBinError(bincount + 1, fMCHist_Slices[i]->GetBinError(j + 1));
+
bincount++;
}
}
+ if (fMCHist_Modes) {
+ fMCHist_Modes->Reset();
+
+ for (std::map >::iterator it =
+ fMCModeHists_Slices.begin();
+ it != fMCModeHists_Slices.end(); ++it) {
+
+ bincount = 0;
+ for (size_t i = 0; i < nangbins; i++) {
+ it->second[i]->Scale(fScaleFactor / (angular_binning_costheta[i + 1] -
+ angular_binning_costheta[i]),
+ "width");
+
+ for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) {
+ fMCHist_Modes->SetBinContent(it->first, bincount + 1, 0, 0,
+ it->second[i]->GetBinContent(j + 1));
+ fMCHist_Modes->SetBinError(it->first, bincount + 1, 0, 0,
+ it->second[i]->GetBinError(j + 1));
+
+ bincount++;
+
+ }
+ }
+ }
+ }
+
return;
}
void T2K_CC0pi_XSec_2DPcos_nu_I::FillMCSlice(double x, double y, double w) {
+ if (fMCHist_Modes && !fMCModeHists_Slices.count(fMode)) {
+ for (size_t i = 0; i < nangbins; ++i) {
+ std::stringstream ss("");
+ ss << "T2K_CC0pi_XSec_2DPcos_nu_I_MODE_ " << fMode << "_slice" << i
+ << std::endl;
+ fMCModeHists_Slices[fMode].push_back(
+ static_cast(fMCHist_Slices[i]->Clone(ss.str().c_str())));
+ fMCModeHists_Slices[fMode].back()->Reset();
+ }
+ }
+
for (size_t i = 0; i < nangbins; ++i) {
if ((y >= angular_binning_costheta[i]) &&
(y < angular_binning_costheta[i + 1])) {
fMCHist_Slices[i]->Fill(x, w);
+ if (fMCHist_Modes) {
+ fMCModeHists_Slices[fMode][i]->Fill(x, w);
+ }
}
}
}
void T2K_CC0pi_XSec_2DPcos_nu_I::SetHistograms() {
// Read in 1D Data Histograms
TFile input(
(FitPar::GetDataBase() + "/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root")
.c_str(),
"READ");
fMCHist_Fine2D = new TH2D("T2K_CC0pi_XSec_2DPcos_nu_I_Fine2D",
"T2K_CC0pi_XSec_2DPcos_nu_I_Fine2D", 400, 0.0, 30.0,
100, -1.0, 1.0);
fMCHist_Fine2D->SetDirectory(NULL);
SetAutoProcessTH1(fMCHist_Fine2D);
TH2D *tempcov = (TH2D *)input.Get("analysis1_totcov");
fFullCovar = new TMatrixDSym(tempcov->GetNbinsX());
for (int i = 0; i < tempcov->GetNbinsX(); i++) {
for (int j = 0; j < tempcov->GetNbinsX(); j++) {
(*fFullCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1);
}
}
// Read in 2D Data Slices and Make MC Slices
int bincount = 0;
for (size_t i = 0; i < nangbins; i++) {
fDataHist_Slices.push_back(
(TH1D *)input.Get(Form("dataslice_%lu", i))->Clone());
fDataHist_Slices[i]->SetNameTitle(
Form("T2K_CC0pi_XSec_2DPcos_nu_I_data_Slice%lu", i),
(Form("T2K_CC0pi_XSec_2DPcos_nu_I_data_Slice%lu, cos(#theta) [%f,%f] ",
i, angular_binning_costheta[i],
angular_binning_costheta[i + 1])));
fDataHist_Slices.back()->SetDirectory(NULL);
// Loop over nbins and set errors from covar
for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) {
fDataHist_Slices[i]->SetBinError(
j + 1, sqrt((*fFullCovar)(bincount, bincount)) * 1E-38);
bincount++;
}
}
assert(bincount == tempcov->GetNbinsX());
if (fMaskMomOverflow) {
MaskMomOverflow();
bincount = fFullCovar->GetNcols();
}
std::vector > data_slice_bcbes;
for (size_t i = 0; i < nangbins; i++) {
for (int j = 0; j < fDataHist_Slices[i]->GetNbinsX(); j++) {
data_slice_bcbes.push_back(
std::make_pair(fDataHist_Slices[i]->GetBinContent(j + 1),
fDataHist_Slices[i]->GetBinError(j + 1)));
}
}
for (size_t i = 0; i < nangbins; i++) {
fMCHist_Slices.push_back((TH1D *)fDataHist_Slices[i]->Clone());
fMCHist_Slices.back()->SetDirectory(NULL);
fMCHist_Slices.back()->Reset();
fMCHist_Slices.back()->SetLineColor(kRed);
fMCHist_Slices[i]->SetNameTitle(
Form("T2K_CC0pi_XSec_2DPcos_nu_I_MC_Slice%lu", i),
(Form("T2K_CC0pi_XSec_2DPcos_nu_I_MC_Slice%lu, cos(#theta) [%f,%f] ", i,
angular_binning_costheta[i], angular_binning_costheta[i + 1])));
}
covar = StatUtils::GetInvert(fFullCovar);
fDecomp = StatUtils::GetDecomp(fFullCovar);
fDataHist =
new TH1D("T2K_CC0pi_XSec_2DPcos_nu_I_DATA_1D",
"T2K_CC0pi_XSec_2DPcos_nu_I_DATA_1D", bincount, 0, bincount);
fDataHist->SetDirectory(NULL);
for (size_t i = 0; i < data_slice_bcbes.size(); ++i) {
fDataHist->SetBinContent(i + 1, data_slice_bcbes[i].first);
fDataHist->SetBinError(i + 1, data_slice_bcbes[i].second);
}
SetShapeCovar();
fMCHist = (TH1D *)fDataHist->Clone("T2K_CC0pi_XSec_2DPcos_nu_I_MC_1D");
fMCHist->SetDirectory(NULL);
return;
}
void T2K_CC0pi_XSec_2DPcos_nu_I::MaskMomOverflow() {
std::vector bins_to_cut;
size_t nallbins = 0;
for (size_t i = 0; i < nangbins; i++) {
std::vector slice_bin_edges;
slice_bin_edges.push_back(
fDataHist_Slices[i]->GetXaxis()->GetBinLowEdge(1));
for (int j = 0; j < (fDataHist_Slices[i]->GetNbinsX() - 1); j++) {
slice_bin_edges.push_back(
fDataHist_Slices[i]->GetXaxis()->GetBinUpEdge(j + 1));
nallbins++;
}
bins_to_cut.push_back(nallbins++);
TH1D *tmp = new TH1D(fDataHist_Slices[i]->GetName(),
fDataHist_Slices[i]->GetTitle(),
slice_bin_edges.size() - 1, slice_bin_edges.data());
tmp->SetDirectory(NULL);
for (int j = 0; j < (fDataHist_Slices[i]->GetNbinsX() - 1); j++) {
tmp->SetBinContent(j + 1, fDataHist_Slices[i]->GetBinContent(j + 1));
tmp->SetBinError(j + 1, fDataHist_Slices[i]->GetBinError(j + 1));
}
delete fDataHist_Slices[i];
fDataHist_Slices[i] = tmp;
}
TMatrixDSym *tmpcovar = new TMatrixDSym(nallbins - bins_to_cut.size());
int icut = 0;
for (int ifull = 0; ifull < fFullCovar->GetNcols(); ifull++) {
if (std::find(bins_to_cut.begin(), bins_to_cut.end(), ifull) !=
bins_to_cut.end()) {
continue;
}
int jcut = 0;
for (int jfull = 0; jfull < fFullCovar->GetNcols(); jfull++) {
if (std::find(bins_to_cut.begin(), bins_to_cut.end(), jfull) !=
bins_to_cut.end()) {
continue;
}
(*tmpcovar)[icut][jcut] = (*fFullCovar)[ifull][jfull];
jcut++;
}
icut++;
}
delete fFullCovar;
fFullCovar = tmpcovar;
}
void T2K_CC0pi_XSec_2DPcos_nu_I::Write(std::string drawopt) {
this->Measurement1D::Write(drawopt);
for (size_t i = 0; i < nangbins; i++) {
fMCHist_Slices[i]->Write();
fDataHist_Slices[i]->Write();
}
if (fResidualHist) {
std::vector MCResidual_Slices;
size_t tb_it = 0;
for (size_t i = 0; i < fMCHist_Slices.size(); ++i) {
std::string name =
Form("T2K_CC0pi_XSec_2DPcos_nu_I_RESIDUAL_Slice%lu", i);
MCResidual_Slices.push_back(
static_cast(fMCHist_Slices[i]->Clone(name.c_str())));
MCResidual_Slices.back()->Reset();
for (int j = 0; j < fMCHist_Slices[i]->GetXaxis()->GetNbins(); ++j) {
double bc = fResidualHist->GetBinContent(tb_it + 1);
MCResidual_Slices.back()->SetBinContent(j + 1, bc);
tb_it++;
}
MCResidual_Slices.back()->Write();
}
}
if (fChi2LessBinHist) {
std::vector MCChi2LessBin_Slices;
size_t tb_it = 0;
for (size_t i = 0; i < fMCHist_Slices.size(); ++i) {
std::string name =
Form("T2K_CC0pi_XSec_2DPcos_nu_I_Chi2NMinusOne_Slice%lu", i);
MCChi2LessBin_Slices.push_back(
static_cast(fMCHist_Slices[i]->Clone(name.c_str())));
MCChi2LessBin_Slices.back()->Reset();
for (int j = 0; j < fMCHist_Slices[i]->GetXaxis()->GetNbins(); ++j) {
double bc = fChi2LessBinHist->GetBinContent(tb_it + 1);
MCChi2LessBin_Slices.back()->SetBinContent(j + 1, bc);
tb_it++;
}
MCChi2LessBin_Slices.back()->Write();
}
}
}
diff --git a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h
index 2570f1e..63f7a26 100644
--- a/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h
+++ b/src/T2K/T2K_CC0pi_XSec_2DPcos_nu_I.h
@@ -1,72 +1,74 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
-* This file is part of NUISANCE.
-*
-* NUISANCE is free software: you can redistribute it and/or modify
-* it under the terms of the GNU General Public License as published by
-* the Free Software Foundation, either version 3 of the License, or
-* (at your option) any later version.
-*
-* NUISANCE is distributed in the hope that it will be useful,
-* but WITHOUT ANY WARRANTY; without even the implied warranty of
-* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-* GNU General Public License for more details.
-*
-* You should have received a copy of the GNU General Public License
-* along with NUISANCE. If not, see .
-*******************************************************************************/
+ * This file is part of NUISANCE.
+ *
+ * NUISANCE is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * NUISANCE is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with NUISANCE. If not, see .
+ *******************************************************************************/
#ifndef T2K_CC0PI_2DPCOS_NU_NONUNIFORM_H_SEEN
#define T2K_CC0PI_2DPCOS_NU_NONUNIFORM_H_SEEN
#include "Measurement1D.h"
-#include "TH2Poly.h"
#include "MeasurementVariableBox2D.h"
+#include "TH2Poly.h"
class T2K_CC0pi_XSec_2DPcos_nu_I : public Measurement1D {
public:
-
T2K_CC0pi_XSec_2DPcos_nu_I(nuiskey samplekey);
/// Virtual Destructor
- ~T2K_CC0pi_XSec_2DPcos_nu_I() {};
+ ~T2K_CC0pi_XSec_2DPcos_nu_I(){};
/// Numu CC0PI Signal Definition
///
/// /item
bool isSignal(FitEvent *nvect);
/// Read histograms in a special way because format is different.
- /// Read from FitPar::GetDataBase()+"/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root"
+ /// Read from
+ /// FitPar::GetDataBase()+"/T2K/CC0pi/T2K_CC0PI_2DPmuCosmu_Data.root"
void SetHistograms();
/// Bin Tmu CosThetaMu
- void FillEventVariables(FitEvent* customEvent);
+ void FillEventVariables(FitEvent *customEvent);
// Fill Histograms
void FillHistograms();
/// Have to do a weird event scaling for analysis 1
void ConvertEventRates();
void Write(std::string drawopt);
/// \brief Create Q2 Box to save correction info
- inline MeasurementVariableBox* CreateBox(){ return new MeasurementVariableBox2D(); };
-
- private:
+ inline MeasurementVariableBox *CreateBox() {
+ return new MeasurementVariableBox2D();
+ };
+private:
bool fIsSystCov, fIsStatCov, fIsNormCov;
bool fMaskMomOverflow;
+ int fMode;
- TH2D* fMCHist_Fine2D;
+ TH2D *fMCHist_Fine2D;
- std::vector fMCHist_Slices;
- std::vector fDataHist_Slices;
+ std::vector fMCHist_Slices;
+ std::map > fMCModeHists_Slices;
+ std::vector fDataHist_Slices;
void FillMCSlice(double x, double y, double w);
void MaskMomOverflow();
-
};
#endif