Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F10882082
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Flag For Later
Size
248 KB
Subscribers
None
View Options
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b6c6bff..4a23728 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,183 +1,187 @@
cmake_minimum_required (VERSION 2.6 FATAL_ERROR)
project(ExternalDataFitter)
enable_language(Fortran)
set (ExtFit_VERSION_MAJOR 1)
set (ExtFit_VERSION_MINOR 0) #The q+1'th letter of the alphabet
set (ExtFit_VERSION_REVISION 0)
set (ExtFit_VERSION_STRING "v${ExtFit_VERSION_MAJOR}r${ExtFit_VERSION_MINOR}")
if(${ExtFit_VERSION_REVISION} STRGREATER "0")
set (ExtFit_VERSION_STRING "${ExtFit_VERSION_STRING}p${ExtFit_VERSION_REVISION}")
endif()
set (VERBOSE TRUE)
set (CMAKE_SKIP_BUILD_RPATH TRUE)
if(NOT DEFINED NOTEST OR NOT NOTEST)
enable_testing()
endif()
include(${CMAKE_SOURCE_DIR}/cmake/cmessage.cmake)
if(NOT DEFINED USE_NEUT AND
NOT DEFINED USE_NuWro AND
NOT DEFINED USE_GENIE AND
NOT DEFINED USE_T2K AND
NOT DEFINED USE_NIWG AND
NOT DEFINED USE_GiBUU AND
NOT DEFINED USE_NUANCE)
cmessage(FATAL_ERROR "No reweight engines requested. Configure with at least "
"one of -DUSE_{NEUT,NuWro,GENIE,NIWG,T2K,GiBUU,NUANCE}.")
else()
cmessage(STATUS "Generator Input Support:
NEUT:${USE_NEUT},
NuWro:${USE_NuWro},
GENIE:${USE_GENIE},
NIWG:${USE_NIWG},
GiBUU:${USE_GiBUU},
NUANCE:${USE_NUANCE}")
endif()
#Changes default install path to be a subdirectory of the build dir.
#Can set build dir at configure time with -DCMAKE_INSTALL_PREFIX=/install/path
if(CMAKE_INSTALL_PREFIX STREQUAL "" OR CMAKE_INSTALL_PREFIX STREQUAL
"/usr/local")
set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}")
elseif(NOT DEFINED CMAKE_INSTALL_PREFIX)
set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}")
endif()
cmessage(STATUS "CMAKE_INSTALL_PREFIX: \"${CMAKE_INSTALL_PREFIX}\"")
if(CMAKE_BUILD_TYPE STREQUAL "")
set(CMAKE_BUILD_TYPE DEBUG)
elseif(NOT DEFINED CMAKE_BUILD_TYPE)
set(CMAKE_BUILD_TYPE DEBUG)
endif()
cmessage(STATUS "CMAKE_BUILD_TYPE: \"${CMAKE_BUILD_TYPE}\"")
################################################################################
# Check Dependencies
################################################################################
################################## ROOT ######################################
include(${CMAKE_SOURCE_DIR}/cmake/ROOTSetup.cmake)
############################ Reweight Engines ################################
include(${CMAKE_SOURCE_DIR}/cmake/ReweightEnginesSetup.cmake)
############################### GiBUU + NUANCE ####################################
if(DEFINED USE_EXP AND USE_EXP)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DINMEMORYEVENTCLASS")
endif()
if(DEFINED USE_GiBUU AND USE_GiBUU)
cmessage(STATUS "Included GiBUU")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -D__GiBUU_ENABLED__")
include(${CMAKE_SOURCE_DIR}/cmake/GiBUUSetup.cmake)
endif()
+if(NOT DEFINED BUILD_GiBUU)
+ set(BUILD_GiBUU 0)
+endif()
+
if(DEFINED USE_NUANCE AND USE_NUANCE)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -D__NUANCE_ENABLED__")
endif()
################################# Pythia6 ####################################
include(${CMAKE_SOURCE_DIR}/cmake/pythia6Setup.cmake)
################################## COMPILER ####################################
include(${CMAKE_SOURCE_DIR}/cmake/c++CompilerSetup.cmake)
################################################################################
################################# gperftools ###################################
include(${CMAKE_SOURCE_DIR}/cmake/gperfSetup.cmake)
################################### doxygen ###################################
include(${CMAKE_SOURCE_DIR}/cmake/docsSetup.cmake)
###############################################################################
set(MINCODE
Routines
FCN)
set(CORE
MCStudies
Genie
FitBase
InputHandler
Splines
Reweight
Utils
#Devel
)
###############
# Allow compilation against single experiment folder
# Add later..
##############
set(EXPERIMENTS
ANL
ArgoNeuT
BEBC
BNL
Electron
FNAL
GGM
K2K
MINERvA
MiniBooNE
SciBooNE
T2K)
set(EXP_INCLUDE_DIRECTORIES)
foreach(edir ${EXPERIMENTS})
set(EXP_INCLUDE_DIRECTORIES ${EXP_INCLUDE_DIRECTORIES};${CMAKE_SOURCE_DIR}/src/${edir})
endforeach()
cmessage(STATUS "Included experiments: ${EXP_INCLUDE_DIRECTORIES}")
foreach(mdir ${MINCODE})
cmessage (STATUS "Configuring directory: src/${mdir}")
add_subdirectory(src/${mdir})
endforeach()
foreach(edir ${EXPERIMENTS})
cmessage (STATUS "Configuring directory: src/${edir}")
add_subdirectory(src/${edir})
endforeach()
foreach(cdir ${CORE})
cmessage (STATUS "Configuring directory: src/${cdir}")
add_subdirectory(src/${cdir})
endforeach()
cmessage(STATUS "Module targets: ${MODULETargets}")
add_subdirectory(app)
add_subdirectory(src/Tests)
configure_file(cmake/setup.sh.in
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.sh" @ONLY)
install(FILES
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.sh" DESTINATION
${CMAKE_INSTALL_PREFIX})
install(PROGRAMS
"${PROJECT_SOURCE_DIR}/scripts/nuiscardgen" DESTINATION
bin)
install(PROGRAMS
"${PROJECT_SOURCE_DIR}/scripts/nuissamples" DESTINATION
bin)
diff --git a/data/MiniBooNE/anti-ccqe/aski_bkg_ccpim.txt b/data/MiniBooNE/anti-ccqe/aski_bkg_ccpim.txt
index 76e47c4..f068322 100755
--- a/data/MiniBooNE/anti-ccqe/aski_bkg_ccpim.txt
+++ b/data/MiniBooNE/anti-ccqe/aski_bkg_ccpim.txt
@@ -1,20 +1,20 @@
-114.2 235.2 315.1 362.5 397.5 410.2 405.5 375.6 333.6 283.6 226.5 176.2 130.3 90.65 0.0 0.0 0.0 0.0
-91.75 170.2 194.0 190.6 177.6 148.7 119.5 91.94 61.22 40.46 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-67.57 110.1 110.0 91.36 67.25 44.58 28.53 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-48.98 70.90 60.55 40.69 24.32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-35.74 43.29 31.14 17.10 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-25.59 27.80 15.67 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-18.96 17.05 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-12.54 9.613 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+85.72 183.4 250.1 288.5 316.7 321.5 311.6 282.6 242.8 199.2 153.7 114.4 57.10 6.809 0.0 0.0 0.0 0.0
+69.90 133.9 153.9 150.5 138.5 113.4 88.08 64.12 39.46 11.47 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+49.96 84.65 85.19 70.93 50.74 31.92 17.41 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+35.68 53.85 45.96 30.16 15.68 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+25.57 31.62 22.11 10.67 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+16.70 19.44 9.089 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+3.479 9.247 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.222 0.567 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
diff --git a/data/MiniBooNE/anti-ccqe/aski_bkg_ccqe.txt b/data/MiniBooNE/anti-ccqe/aski_bkg_ccqe.txt
index f068322..76e47c4 100755
--- a/data/MiniBooNE/anti-ccqe/aski_bkg_ccqe.txt
+++ b/data/MiniBooNE/anti-ccqe/aski_bkg_ccqe.txt
@@ -1,20 +1,20 @@
-85.72 183.4 250.1 288.5 316.7 321.5 311.6 282.6 242.8 199.2 153.7 114.4 57.10 6.809 0.0 0.0 0.0 0.0
-69.90 133.9 153.9 150.5 138.5 113.4 88.08 64.12 39.46 11.47 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-49.96 84.65 85.19 70.93 50.74 31.92 17.41 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-35.68 53.85 45.96 30.16 15.68 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-25.57 31.62 22.11 10.67 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-16.70 19.44 9.089 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-3.479 9.247 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.222 0.567 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
-0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+114.2 235.2 315.1 362.5 397.5 410.2 405.5 375.6 333.6 283.6 226.5 176.2 130.3 90.65 0.0 0.0 0.0 0.0
+91.75 170.2 194.0 190.6 177.6 148.7 119.5 91.94 61.22 40.46 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+67.57 110.1 110.0 91.36 67.25 44.58 28.53 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+48.98 70.90 60.55 40.69 24.32 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+35.74 43.29 31.14 17.10 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+25.59 27.80 15.67 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+18.96 17.05 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+12.54 9.613 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
+0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
diff --git a/src/BNL/BNL_CC1ppip_Evt_1Dphi_nu.cxx b/src/BNL/BNL_CC1ppip_Evt_1Dphi_nu.cxx
index e33d926..977e295 100644
--- a/src/BNL/BNL_CC1ppip_Evt_1Dphi_nu.cxx
+++ b/src/BNL/BNL_CC1ppip_Evt_1Dphi_nu.cxx
@@ -1,159 +1,159 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "BNL_CC1ppip_Evt_1Dphi_nu.h"
//********************************************************************
BNL_CC1ppip_Evt_1Dphi_nu::BNL_CC1ppip_Evt_1Dphi_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "BNL_CC1ppip_Evt_1Dphi_nu sample. \n" \
"Target: D2 \n" \
"Flux: \n" \
"Signal: \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#phi_{Adler}");
fSettings.SetYTitle("Number of events");
fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG");
- fSettings.SetEnuRange(0.5, 6.0);
+ fSettings.SetEnuRange(0.0, 6.0);
fSettings.DefineAllowedTargets("D,H");
// CCQELike plot information
fSettings.SetTitle("BNL_CC1ppip_Evt_1Dphi_nu");
fSettings.SetDataInput( FitPar::GetDataBase() + "/BNL/CC1pip_on_p/BNL_CC1ppip_W14_phiAdler.csv" );
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor =GetEventHistogram()->Integral("width")/(fNEvents+0.)*2./1.;
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetPoissonErrors();
SetCovarFromDiagonal();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void BNL_CC1ppip_Evt_1Dphi_nu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(211) == 0 ||
event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pp = event->GetHMFSParticle(2212)->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
// Get the hadronic mass
double hadMass = FitUtils::MpPi(Pp, Ppip);
// Need to boost pion and muon into resonance rest-frame to get phi (e.g. see F. Sanchez arxiv 1511.00501v2)
//
// Get the resonance 4-vector
TLorentzVector Pres = Ppip + Pp;
// Boost the pion 4-vector into the resonance 4-vector rest-frame
Ppip.Boost(Pres.BoostVector());
Pmu.Boost(-Pres.BoostVector());
Pp.Boost(-Pres.BoostVector());
// Get the vectors from the 4-vector
TVector3 PmuVect = Pmu.Vect();
TVector3 PnuVect = Pnu.Vect();
TVector3 PresVect = Pres.Vect();
TVector3 PpipVect = Ppip.Vect();
// Define y direction as being z (resonance direction) x pmu*
TVector3 zVect = (PnuVect-PmuVect);
zVect *= 1/double(zVect.Mag());
TVector3 yVect = zVect.Cross(PmuVect);
// Normalise yVector
yVect *= 1/double(yVect.Mag());
// define x direction as being y X z
TVector3 xVect = yVect.Cross(zVect);
// Normalise zVector
xVect *= 1/double(xVect.Mag());
// Project pion onto z axis
TVector3 PpipVectZ = zVect * PpipVect.Dot(zVect);
// Then subtract this vector off the pion vector
TVector3 PpipVectPlane = PpipVect - PpipVectZ;
// Then finally construct phi as the angle between pion projection and x axis
double phi = -999;
// BNL has a M(pi, p) < 1.4 GeV cut imposed
if (hadMass < 1400) {
if (PpipVectPlane.Y() > 0) {
phi = (180./M_PI)*PpipVectPlane.Angle(xVect);
} else if (PpipVectPlane.Y() < 0) {
phi = (180./M_PI)*(2*M_PI-PpipVectPlane.Angle(xVect));
} else if (PpipVectPlane.Y() == 0) {
double randNo = rand.Rndm();
if (randNo > 0.5) {
phi = (180./M_PI)*PpipVectPlane.Angle(xVect);
} else {
phi = (180./M_PI)*(2*M_PI-PpipVectPlane.Angle(xVect));
}
}
}
fXVar = phi;
return;
};
bool BNL_CC1ppip_Evt_1Dphi_nu::isSignal(FitEvent *event) {
return SignalDef::isCC1pi3Prong(event, 14, 211, 2212,EnuMin,EnuMax);
}
/*
void BNL_CC1ppip_Evt_1Dphi_nu::FillHistograms() {
if (makeHadronicMassHist) {
hadMassHist->Fill(hadMass);
}
Measurement1D::FillHistograms();
}
void BNL_CC1ppip_Evt_1Dphi_nu::ScaleEvents() {
PlotUtils::FluxUnfoldedScaling(fMCHist, GetFluxHistogram());
PlotUtils::FluxUnfoldedScaling(fMCFine, GetFluxHistogram());
fMCHist->Scale(fScaleFactor);
fMCFine->Scale(fScaleFactor);
return;
}
*/
diff --git a/src/FCN/JointFCN.cxx b/src/FCN/JointFCN.cxx
index 088ecd1..0d2d9b0 100755
--- a/src/FCN/JointFCN.cxx
+++ b/src/FCN/JointFCN.cxx
@@ -1,1036 +1,1032 @@
#include "JointFCN.h"
#include <stdio.h>
#include "FitUtils.h"
// //***************************************************
// JointFCN::JointFCN(std::string cardfile, TFile* outfile) {
// //***************************************************
// fOutputDir = gDirectory;
// FitPar::Config().out = outfile;
// fCardFile = cardfile;
// LoadSamples(fCardFile);
// fCurIter = 0;
// fMCFilled = false;
// fOutputDir->cd();
// fIterationTree = NULL;
// fDialVals = NULL;
// fNDials = 0;
// fUsingEventManager = FitPar::Config().GetParB("EventManager");
// fOutputDir->cd();
// };
JointFCN::JointFCN(TFile* outfile) {
fOutputDir = gDirectory;
if (outfile)
FitPar::Config().out = outfile;
std::vector<nuiskey> samplekeys = Config::QueryKeys("sample");
LoadSamples(samplekeys);
std::vector<nuiskey> covarkeys = Config::QueryKeys("covar");
LoadPulls(covarkeys);
fCurIter = 0;
fMCFilled = false;
fIterationTree = NULL;
fDialVals = NULL;
fNDials = 0;
fUsingEventManager = FitPar::Config().GetParB("EventManager");
fOutputDir->cd();
}
JointFCN::JointFCN(std::vector<nuiskey> samplekeys, TFile* outfile) {
fOutputDir = gDirectory;
if (outfile)
FitPar::Config().out = outfile;
LoadSamples(samplekeys);
fCurIter = 0;
fMCFilled = false;
fOutputDir->cd();
fIterationTree = NULL;
fDialVals = NULL;
fNDials = 0;
fUsingEventManager = FitPar::Config().GetParB("EventManager");
fOutputDir->cd();
}
//***************************************************
JointFCN::~JointFCN() {
//***************************************************
// Delete Samples
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
delete exp;
}
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
delete pull;
}
// Sort Tree
if (fIterationTree) DestroyIterationTree();
if (fDialVals) delete fDialVals;
if (fSampleLikes) delete fSampleLikes;
};
//***************************************************
void JointFCN::CreateIterationTree(std::string name, FitWeight* rw) {
//***************************************************
LOG(FIT) << " Creating new iteration tree! " << std::endl;
if (fIterationTree && !name.compare(fIterationTree->GetName())) {
fIterationTree->Reset();
return;
}
fIterationTree = new TTree(name.c_str(), name.c_str());
// Setup Main Branches
fIterationTree->Branch("total_likelihood", &fLikelihood,
"total_likelihood/D");
fIterationTree->Branch("total_ndof", &fNDOF, "total_ndof/I");
// Setup Sample Arrays
int ninputs = fSamples.size() + fPulls.size();
fSampleLikes = new double[ninputs];
fSampleNDOF = new int[ninputs];
fNDOF = GetNDOF();
// Setup Sample Branches
int count = 0;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
std::string name = exp->GetName();
std::string liketag = name + "_likelihood";
std::string ndoftag = name + "_ndof";
fIterationTree->Branch(liketag.c_str(), &fSampleLikes[count],
(liketag + "/D").c_str());
fIterationTree->Branch(ndoftag.c_str(), &fSampleNDOF[count],
(ndoftag + "/D").c_str());
count++;
}
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
std::string name = pull->GetName();
std::string liketag = name + "_likelihood";
std::string ndoftag = name + "_ndof";
fIterationTree->Branch(liketag.c_str(), &fSampleLikes[count],
(liketag + "/D").c_str());
fIterationTree->Branch(ndoftag.c_str(), &fSampleNDOF[count],
(ndoftag + "/D").c_str());
count++;
}
// Add Dial Branches
std::vector<std::string> dials = rw->GetDialNames();
fNDials = dials.size();
fDialVals = new double[fNDials];
for (int i = 0; i < fNDials; i++) {
fIterationTree->Branch(dials[i].c_str(), &fDialVals[i],
(dials[i] + "/D").c_str());
}
}
//***************************************************
void JointFCN::DestroyIterationTree() {
//***************************************************
if (!fIterationTree) {
delete fIterationTree;
}
}
//***************************************************
void JointFCN::WriteIterationTree() {
//***************************************************
if (!fIterationTree) {
ERR(FTL) << "Can't save empty iteration tree!" << std::endl;
throw;
}
fIterationTree->Write();
}
//***************************************************
void JointFCN::FillIterationTree(FitWeight* rw) {
//***************************************************
if (!fIterationTree) {
ERR(FTL) << "Trying to fill iteration_tree when it is NULL!" << std::endl;
throw;
}
rw->GetAllDials(fDialVals, fNDials);
fIterationTree->Fill();
}
//***************************************************
double JointFCN::DoEval(const double* x) {
//***************************************************
// WEIGHT ENGINE
fDialChanged = FitBase::GetRW()->HasRWDialChanged(x);
FitBase::GetRW()->UpdateWeightEngine(x);
if (fDialChanged) {
FitBase::GetRW()->Reconfigure();
FitBase::EvtManager().ResetWeightFlags();
}
if (LOG_LEVEL(REC)) {
FitBase::GetRW()->Print();
}
// SORT SAMPLES
ReconfigureSamples();
// GET TEST STAT
fLikelihood = GetLikelihood();
// PRINT PROGRESS
LOG(FIT) << "Current Stat (iter. " << this->fCurIter << ") = " << fLikelihood
<< std::endl;
// UPDATE TREE
if (fIterationTree) FillIterationTree(FitBase::GetRW());
return fLikelihood;
}
//***************************************************
int JointFCN::GetNDOF() {
//***************************************************
int totaldof = 0;
int count = 0;
// Total number of Free bins in each MC prediction
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
int dof = exp->GetNDOF();
// Save Seperate DOF
if (fIterationTree) {
fSampleNDOF[count] = dof;
}
// Add to total
totaldof += dof;
count++;
}
// Loop over pulls
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
double dof = pull->GetLikelihood();
// Save seperate DOF
if (fIterationTree) {
fSampleNDOF[count] = dof;
}
// Add to total
totaldof += dof;
count++;
}
// Set Data Variable
fNDOF = totaldof;
return totaldof;
}
//***************************************************
double JointFCN::GetLikelihood() {
//***************************************************
LOG(MIN) << std::left << std::setw(43) << "Getting likelihoods..." << " : " << "-2logL" << std::endl;
// Loop and add up likelihoods in an uncorrelated way
double like = 0.0;
int count = 0;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
double newlike = exp->GetLikelihood();
// Save seperate likelihoods
if (fIterationTree) {
fSampleLikes[count] = newlike;
}
LOG(MIN) << "-> " << std::left << std::setw(40) << exp->GetName() << " : " << newlike << std::endl;
// Add Weight Scaling
// like *= FitBase::GetRW()->GetSampleLikelihoodWeight(exp->GetName());
// Add to total
like += newlike;
count++;
}
// Loop over pulls
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
double newlike = pull->GetLikelihood();
// Save seperate likelihoods
if (fIterationTree) {
fSampleLikes[count] = newlike;
}
// Add to total
like += newlike;
count++;
}
// Set Data Variable
fLikelihood = like;
return like;
};
void JointFCN::LoadSamples(std::vector<nuiskey> samplekeys) {
LOG(MIN) << "Loading Samples : " << samplekeys.size() << std::endl;
for (size_t i = 0; i < samplekeys.size(); i++) {
nuiskey key = samplekeys[i];
// Get Sample Options
std::string samplename = key.GetS("name");
std::string samplefile = key.GetS("input");
std::string sampletype = key.GetS("type");
std::string fakeData = "";
LOG(MIN) << "Loading Sample : " << samplename << std::endl;
fOutputDir->cd();
MeasurementBase* NewLoadedSample
= SampleUtils::CreateSample(key);
if (!NewLoadedSample) {
ERR(FTL) << "Could not load sample provided: " << samplename << std::endl;
ERR(FTL) << "Check spelling with that in src/FCN/SampleList.cxx"
<< std::endl;
throw;
} else {
fSamples.push_back(NewLoadedSample);
}
}
}
void JointFCN::LoadPulls(std::vector<nuiskey> pullkeys) {
for (size_t i = 0; i < pullkeys.size(); i++) {
nuiskey key = pullkeys[i];
std::string pullname = key.GetS("name");
std::string pullfile = key.GetS("input");
std::string pulltype = key.GetS("type");
fOutputDir->cd();
fPulls.push_back(new ParamPull(pullname, pullfile, pulltype));
}
// // Sample Inputs
// if (!samplevect[0].compare("covar") || !samplevect[0].compare("pulls") ||
// !samplevect[0].compare("throws")) {
// // Get all inputs
// std::string name = samplevect[1];
// std::string files = samplevect[2];
// std::string type = "DEFAULT";
// if (samplevect.size() > 3) {
// type = samplevect[3];
// } else if (!samplevect[0].compare("pull")) {
// type = "GAUSPULL";
// } else if (!samplevect[0].compare("throws")) {
// type = "GAUSTHROWS";
// }
// // Create Pull Class
// LOG(MIN) << "Loading up pull term: " << name << " << " << files << " ("
// << type << ")" << std::endl;
// std::string fakeData = "";
// fOutputDir->cd();
// fPulls.push_back(new ParamPull(name, files, type));
// }
}
// //***************************************************
// void JointFCN::LoadSamples(std::string cardinput)
// //***************************************************
// {
// LOG(MIN) << "Initializing Samples" << std::endl;
// // Read the card file here and load objects
// std::string line;
// std::ifstream card(cardinput.c_str(), ifstream::in);
// // Make sure they are created in correct working DIR
// fOutputDir->cd();
// while (std::getline(card >> std::ws, line, '\n')) {
// // Skip Empties
// if (line.c_str()[0] == '#') continue;
// if (line.empty()) continue;
// // Parse line
// std::vector<std::string> samplevect = GeneralUtils::ParseToStr(line, " ");
// // Sample Inputs
// if (!samplevect[0].compare("sample")) {
// // Get all inputs
// std::string name = samplevect[1];
// std::string files = samplevect[2];
// std::string type = "DEFAULT";
// if (samplevect.size() > 3) type = samplevect[3];
// // Create Sample Class
// LOG(MIN) << "Loading up sample: " << name << " << " << files << " ("
// << type << ")" << std::endl;
// std::string fakeData = "";
// fOutputDir->cd();
// MeasurementBase* NewLoadedSample = SampleUtils::CreateSample(name, files, type,
// fakeData, FitBase::GetRW());
// if (!NewLoadedSample) {
// ERR(FTL) << "Could not load sample provided: " << name << std::endl;
// ERR(FTL) << "Check spelling with that in src/FCN/SampleList.cxx"
// << std::endl;
// throw;
// } else {
// fSamples.push_back(NewLoadedSample);
// }
// }
// // Sample Inputs
// if (!samplevect[0].compare("covar") || !samplevect[0].compare("pulls") ||
// !samplevect[0].compare("throws")) {
// // Get all inputs
// std::string name = samplevect[1];
// std::string files = samplevect[2];
// std::string type = "DEFAULT";
// if (samplevect.size() > 3) {
// type = samplevect[3];
// } else if (!samplevect[0].compare("pull")) {
// type = "GAUSPULL";
// } else if (!samplevect[0].compare("throws")) {
// type = "GAUSTHROWS";
// }
// // Create Pull Class
// LOG(MIN) << "Loading up pull term: " << name << " << " << files << " ("
// << type << ")" << std::endl;
// std::string fakeData = "";
// fOutputDir->cd();
// fPulls.push_back(new ParamPull(name, files, type));
// }
// }
// card.close();
// };
//***************************************************
void JointFCN::ReconfigureSamples(bool fullconfig) {
//***************************************************
int starttime = time(NULL);
LOG(REC) << "Starting Reconfigure iter. " << this->fCurIter << std::endl;
// std::cout << fUsingEventManager << " " << fullconfig << " " << fMCFilled << std::endl;
// Event Manager Reconf
if (fUsingEventManager) {
if (!fullconfig and fMCFilled)
ReconfigureFastUsingManager();
else
ReconfigureUsingManager();
} else {
// Loop over all Measurement Classes
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
// If RW Either do signal or full reconfigure.
if (fDialChanged or !fMCFilled or fullconfig) {
if (!fullconfig and fMCFilled)
exp->ReconfigureFast();
else
exp->Reconfigure();
// If RW Not needed just do normalisation
} else {
exp->Renormalise();
}
}
}
// Loop over pulls and update
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
pull->Reconfigure();
}
fMCFilled = true;
LOG(MIN) << "Finished Reconfigure iter. " << fCurIter << " in "
<< time(NULL) - starttime << "s" << std::endl;
fCurIter++;
}
//***************************************************
void JointFCN::ReconfigureSignal() {
//***************************************************
this->ReconfigureSamples(false);
}
//***************************************************
void JointFCN::ReconfigureAllEvents() {
//***************************************************
FitBase::GetRW()->Reconfigure();
FitBase::EvtManager().ResetWeightFlags();
ReconfigureSamples(true);
}
std::vector<InputHandlerBase*> JointFCN::GetInputList() {
std::vector<InputHandlerBase*> InputList;
fIsAllSplines = true;
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
std::vector<MeasurementBase*> subsamples = exp->GetSubSamples();
for (size_t i = 0; i < subsamples.size(); i++) {
InputHandlerBase* inp = subsamples[i]->GetInput();
if (std::find(InputList.begin(), InputList.end(), inp) == InputList.end()) {
if (subsamples[i]->GetInput()->GetType() != kSPLINEPARAMETER) fIsAllSplines = false;
InputList.push_back(subsamples[i]->GetInput());
}
}
}
return InputList;
}
std::vector<MeasurementBase*> JointFCN::GetSubSampleList() {
std::vector<MeasurementBase*> SampleList;
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
std::vector<MeasurementBase*> subsamples = exp->GetSubSamples();
for (size_t i = 0; i < subsamples.size(); i++) {
SampleList.push_back(subsamples[i]);
}
}
return SampleList;
}
//***************************************************
void JointFCN::ReconfigureUsingManager() {
//***************************************************
// 'Slow' Event Manager Reconfigure
LOG(REC) << "Event Manager Reconfigure" << std::endl;
int timestart = time(NULL);
// Reset all samples
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ResetAll();
}
// If we are siving signal, reset all containers.
bool savesignal = (FitPar::Config().GetParB("SignalReconfigures"));
if (savesignal) {
// Reset all of our event signal vectors
fSignalEventBoxes.clear();
fSignalEventFlags.clear();
fSampleSignalFlags.clear();
fSignalEventSplines.clear();
}
// Make sure we have a list of inputs
if (fInputList.empty()) {
fInputList = GetInputList();
fSubSampleList = GetSubSampleList();
}
// If all inputs are splines make sure the readers are told
// they need to be reconfigured.
std::vector<InputHandlerBase*>::iterator inp_iter = fInputList.begin();
if (fIsAllSplines) {
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Tell reader in each BaseEvent it needs a Reconfigure next weight calc.
BaseFitEvt* curevent = curinput->FirstBaseEvent();
if (curevent->fSplineRead) {
curevent->fSplineRead->SetNeedsReconfigure(true);
}
}
}
// MAIN INPUT LOOP ====================
int fillcount = 0;
int inputcount = 0;
inp_iter = fInputList.begin();
// Loop over each input in manager
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Get event information
FitEvent* curevent = curinput->FirstNuisanceEvent();
curinput->CreateCache();
int i = 0;
int nevents = curinput->GetNEvents();
- int countwidth = nevents / 20;
+ int countwidth = nevents / 5;
// Start event loop iterating until we get a NULL pointer.
while (curevent) {
// Get Event Weight
curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
curevent->Weight = curevent->RWWeight * curevent->InputWeight;
double rwweight = curevent->Weight;
// std::cout << "RWWeight = " << curevent->RWWeight << " " << curevent->InputWeight << std::endl;
// Logging
// std::cout << CHECKLOG(1) << std::endl;
if (LOGGING(REC)) {
if (i % countwidth == 0) {
QLOG(REC, curinput->GetName() << " : Processed " << i
<< " events. [M, W] = ["
<< curevent->Mode << ", "
<< rwweight << "]");
}
}
// Setup flag for if signal found in at least one sample
bool foundsignal = false;
// Create a new signal bitset for this event
std::vector<bool> signalbitset(fSubSampleList.size());
// Create a new signal box vector for this event
std::vector<MeasurementVariableBox*> signalboxes;
// Start measurement iterator
size_t measitercount = 0;
std::vector<MeasurementBase*>::iterator meas_iter = fSubSampleList.begin();
// Loop over all subsamples (sub in JointMeas)
for (; meas_iter != fSubSampleList.end(); meas_iter++) {
MeasurementBase* curmeas = (*meas_iter);
// Compare input pointers, to current input, skip if not.
// Pointer tells us if it matches without doing ID checks.
if (curinput != curmeas->GetInput()) {
if (savesignal) {
// Set bit to 0 as definitely not signal
signalbitset[measitercount] = 0;
}
// Count up what measurement we are on.
measitercount++;
// Skip sample as input not signal.
continue;
}
// Fill events for matching inputs.
MeasurementVariableBox* box = curmeas->FillVariableBox(curevent);
bool signal = curmeas->isSignal(curevent);
curmeas->SetSignal(signal);
curmeas->FillHistograms(curevent->Weight);
// If its Signal tally up fills
if (signal) {
fillcount++;
}
// If we are saving signal/splines fill the bitset
if (savesignal) {
signalbitset[measitercount] = signal;
}
// If signal save a clone of the event box for use later.
if (savesignal and signal) {
foundsignal = true;
signalboxes.push_back( box->CloneSignalBox() );
}
// Keep track of Measurement we are on.
measitercount++;
}
// Once we've filled the measurements, if saving signal
// push back if any sample flagged this event as signal
if (savesignal) {
fSignalEventFlags.push_back(foundsignal);
}
// Save the vector of signal boxes for this event
if (savesignal and foundsignal) {
fSignalEventBoxes.push_back(signalboxes);
fSampleSignalFlags.push_back(signalbitset);
}
// If all inputs are splines we can save the spline coefficients
// for fast in memory reconfigures later.
if (fIsAllSplines and savesignal and foundsignal) {
// Make temp vector to push back with
std::vector<float> coeff;
for (size_t l = 0; l < (UInt_t)curevent->fSplineRead->GetNPar(); l++) {
coeff.push_back( curevent->fSplineCoeff[l] );
}
// Push back to signal event splines. Kept in sync with fSignalEventBoxes size.
// int splinecount = fSignalEventSplines.size();
fSignalEventSplines.push_back(coeff);
// if (splinecount % 1000 == 0) {
// std::cout << "Pushed Back Coeff " << splinecount << " : ";
// for (size_t l = 0; l < fSignalEventSplines[splinecount].size(); l++) {
// std::cout << " " << fSignalEventSplines[splinecount][l];
// }
// std::cout << std::endl;
// }
}
// Clean up vectors once done with this event
signalboxes.clear();
signalbitset.clear();
// Iterate to the next event.
curevent = curinput->NextNuisanceEvent();
i++;
}
curinput->RemoveCache();
// Keep track of what input we are on.
inputcount++;
}
// End of Event Loop ===============================
// Now event loop is finished loop over all Measurements
// Converting Binned events to XSec Distributions
iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ConvertEventRates();
}
// Print out statements on approximate memory usage for profiling.
LOG(REC) << "Filled " << fillcount << " signal events." << std::endl;
if (savesignal) {
int mem = ( //sizeof(fSignalEventBoxes) +
// fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) +
sizeof(MeasurementVariableBox1D) * fillcount) * 1E-6;
LOG(REC) << " -> Saved " << fillcount << " signal boxes for faster access. (~" << mem << " MB)" << std::endl;
if (fIsAllSplines and !fSignalEventSplines.empty()) {
int splmem = sizeof(float) * fSignalEventSplines.size() * fSignalEventSplines[0].size() * 1E-6;
LOG(REC) << " -> Saved " << fillcount << " " << fSignalEventSplines.size() << " spline sets into memory. (~" << splmem << " MB)" << std::endl;
}
}
LOG(REC) << "Time taken ReconfigureUsingManager() : " << time(NULL) - timestart << std::endl;
// End of reconfigure
return;
};
//***************************************************
void JointFCN::ReconfigureFastUsingManager() {
//***************************************************
LOG(FIT) << " -> Doing FAST using manager" << std::endl;
// Get Start time for profilling
int timestart = time(NULL);
// Reset all samples
MeasListConstIter iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ResetAll();
}
// Check for saved variables if not do a full reconfigure.
if (fSignalEventFlags.empty()) {
ERR(WRN) << "Signal Flags Empty! Using normal manager." << std::endl;
ReconfigureUsingManager();
return;
}
bool fFillNuisanceEvent = FitPar::Config().GetParB("FullEventOnSignalReconfigure");
// Setup fast vector iterators.
std::vector<bool>::iterator inpsig_iter = fSignalEventFlags.begin();
std::vector< std::vector<MeasurementVariableBox*> >::iterator box_iter = fSignalEventBoxes.begin();
std::vector< std::vector<float> >::iterator spline_iter = fSignalEventSplines.begin();
std::vector< std::vector<bool> >::iterator samsig_iter = fSampleSignalFlags.begin();
int splinecount = 0;
// Setup stuff for logging
int fillcount = 0;
int nevents = fSignalEventFlags.size();
- int countwidth = nevents / 500;
+ int countwidth = nevents / 20;
// If All Splines tell splines they need a reconfigure.
std::vector<InputHandlerBase*>::iterator inp_iter = fInputList.begin();
if (fIsAllSplines) {
LOG(REC) << "All Spline Inputs so using fast spline loop." << std::endl;
for (; inp_iter != fInputList.end(); inp_iter++) {
InputHandlerBase* curinput = (*inp_iter);
// Tell each fSplineRead in BaseFitEvent to reconf next weight calc
BaseFitEvt* curevent = curinput->FirstBaseEvent();
if (curevent->fSplineRead) curevent->fSplineRead->SetNeedsReconfigure(true);
}
}
// Loop over all possible spline inputs
double* coreeventweights = new double[fSignalEventBoxes.size()];
splinecount = 0;
inp_iter = fInputList.begin();
inpsig_iter = fSignalEventFlags.begin();
spline_iter = fSignalEventSplines.begin();
// Loop over all signal flags
// For each valid signal flag add one to splinecount
// Get Splines from that count and add to weight
// Add splinecount
int sigcount = 0;
splinecount = 0;
// #pragma omp parallel for shared(splinecount,sigcount)
for (uint iinput = 0; iinput < fInputList.size(); iinput++) {
InputHandlerBase* curinput = fInputList[iinput];
BaseFitEvt* curevent = curinput->FirstBaseEvent();
for (int i = 0; i < curinput->GetNEvents(); i++) {
double rwweight = 0.0;
if (fSignalEventFlags[sigcount]) {
// Get Event Info
if (!fIsAllSplines) {
if (fFillNuisanceEvent) curinput->GetNuisanceEvent(i);
else curevent = curinput->GetBaseEvent(i);
} else {
curevent->fSplineCoeff = &fSignalEventSplines[splinecount][0];
}
curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent);
curevent->Weight = curevent->RWWeight * curevent->InputWeight;
rwweight = curevent->Weight;
- // #pragma omp atomic
- if (fIsAllSplines) {
- coreeventweights[splinecount] = rwweight;
- }
+ coreeventweights[splinecount] = rwweight;
if (splinecount % countwidth == 0) {
- LOG(REC) << "Processed " << splinecount << " event weights." << std::endl;
+ LOG(REC) << "Processed " << splinecount << " event weights. W = " << rwweight << std::endl;
}
-
// #pragma omp atomic
splinecount++;
}
// #pragma omp atomic
sigcount++;
}
}
LOG(SAM) << "Processed event weights." << std::endl;
// #pragma omp barrier
// Reset Iterators
inpsig_iter = fSignalEventFlags.begin();
spline_iter = fSignalEventSplines.begin();
box_iter = fSignalEventBoxes.begin();
samsig_iter = fSampleSignalFlags.begin();
int nsplineweights = splinecount;
splinecount = 0;
// Start of Fast Event Loop ============================
// Start input iterators
// Loop over number of inputs
for (int ispline = 0; ispline < nsplineweights; ispline++) {
double rwweight = coreeventweights[ispline];
// Get iterators for this event
std::vector<bool>::iterator subsamsig_iter = (*samsig_iter).begin();
std::vector<MeasurementVariableBox*>::iterator subbox_iter = (*box_iter).begin();
// Loop over all sub measurements.
std::vector<MeasurementBase*>::iterator meas_iter = fSubSampleList.begin();
for (; meas_iter != fSubSampleList.end(); meas_iter++, subsamsig_iter++) {
MeasurementBase* curmeas = (*meas_iter);
// If event flagged as signal for this sample fill from the box.
if (*subsamsig_iter) {
curmeas->SetSignal(true);
curmeas->FillHistogramsFromBox((*subbox_iter), rwweight);
// Move onto next box if there is one.
subbox_iter++;
fillcount++;
}
}
if (ispline % countwidth == 0) {
LOG(REC) << "Filled " << ispline << " sample weights." << std::endl;
}
// Iterate over the main signal event containers.
samsig_iter++;
box_iter++;
spline_iter++;
splinecount++;
}
// End of Fast Event Loop ===================
LOG(SAM) << "Filled sample distributions." << std::endl;
// Now loop over all Measurements
// Convert Binned events
iterSam = fSamples.begin();
for (; iterSam != fSamples.end(); iterSam++) {
MeasurementBase* exp = (*iterSam);
exp->ConvertEventRates();
}
// Cleanup coreeventweights
if (fIsAllSplines) {
delete coreeventweights;
}
// Print some reconfigure profiling.
LOG(REC) << "Filled " << fillcount << " signal events." << std::endl;
LOG(REC) << "Time taken ReconfigureFastUsingManager() : " << time(NULL) - timestart << std::endl;
}
//***************************************************
void JointFCN::Write() {
//***************************************************
// Loop over individual experiments and call Write
LOG(MIN) << "Writing each of the data classes..." << std::endl;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
exp->Write();
}
// Save Pull Terms
for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) {
ParamPull* pull = *iter;
pull->Write();
}
if (FitPar::Config().GetParB("EventManager")) {
// Get list of inputs
std::map<int, InputHandlerBase*> fInputs = FitBase::EvtManager().GetInputs();
std::map<int, InputHandlerBase*>::const_iterator iterInp;
for (iterInp = fInputs.begin(); iterInp != fInputs.end(); iterInp++) {
InputHandlerBase* input = (iterInp->second);
input->GetFluxHistogram()->Write();
input->GetXSecHistogram()->Write();
input->GetEventHistogram()->Write();
}
}
};
//***************************************************
void JointFCN::SetFakeData(std::string fakeinput) {
//***************************************************
LOG(MIN) << "Setting fake data from " << fakeinput << std::endl;
for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end();
iter++) {
MeasurementBase* exp = *iter;
exp->SetFakeDataValues(fakeinput);
}
return;
}
diff --git a/src/FNAL/FNAL_CC1ppip_Evt_1DQ2_nu.cxx b/src/FNAL/FNAL_CC1ppip_Evt_1DQ2_nu.cxx
index 1607860..2e5edd6 100644
--- a/src/FNAL/FNAL_CC1ppip_Evt_1DQ2_nu.cxx
+++ b/src/FNAL/FNAL_CC1ppip_Evt_1DQ2_nu.cxx
@@ -1,88 +1,90 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "FNAL_CC1ppip_Evt_1DQ2_nu.h"
//********************************************************************
FNAL_CC1ppip_Evt_1DQ2_nu::FNAL_CC1ppip_Evt_1DQ2_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "FNAL_CC1ppip_Evt_1DQ2_nu sample. \n" \
"Target: D2 \n" \
"Flux: \n" \
"Signal: \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("Q^{2}_{CC1#pi} (GeV^{2})");
fSettings.SetYTitle("Number of events");
- fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG");
+ //fSettings.SetAllowedTypes("EVT/DIAG", "EVT/DIAG");
fSettings.SetEnuRange(10.0, 100.0);
- fSettings.DefineAllowedTargets("D,H");
+ fSettings.DefineAllowedTargets("D,H");
+ fIsDiag = true;
+ fIsRawEvents = true;
// plot information
fSettings.SetTitle("FNAL_CC1ppip_Evt_1DQ2_nu");
fSettings.DefineAllowedSpecies("numu");
fSettings.SetDataInput( FitPar::GetDataBase() + "FNAL/CC1pip_on_p/FNAL_cc1ppip_nEvent_Q2_w14_bin_edit.txt");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor for shape
fScaleFactor = GetEventHistogram()->Integral("width")/(fNEvents+0.)*2./1.;
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetPoissonErrors();
SetCovarFromDiagonal();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
}
void FNAL_CC1ppip_Evt_1DQ2_nu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(211) == 0 ||
event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pp = event->GetHMFSParticle(2212)->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double hadMass = FitUtils::MpPi(Pp, Ppip);
double q2CCpip = -1.0;
// FNAL has a M(pi, p) < 1.4 GeV cut imposed only on this channel
if (hadMass < 1400) q2CCpip = -1*(Pnu-Pmu).Mag2()/1.E6;
fXVar = q2CCpip;
return;
};
bool FNAL_CC1ppip_Evt_1DQ2_nu::isSignal(FitEvent *event) {
return SignalDef::isCC1pi3Prong(event, 14, 211, 2212, EnuMin, EnuMax);
}
diff --git a/src/GGM/GGM_CC1ppip_Evt_1DQ2_nu.cxx b/src/GGM/GGM_CC1ppip_Evt_1DQ2_nu.cxx
index 90fb1b0..82329a0 100644
--- a/src/GGM/GGM_CC1ppip_Evt_1DQ2_nu.cxx
+++ b/src/GGM/GGM_CC1ppip_Evt_1DQ2_nu.cxx
@@ -1,88 +1,90 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "GGM_CC1ppip_Evt_1DQ2_nu.h"
//********************************************************************
GGM_CC1ppip_Evt_1DQ2_nu::GGM_CC1ppip_Evt_1DQ2_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "GGM_CC1ppip_Evt_1DQ2_nu sample. \n" \
"Target: D2 \n" \
"Flux: \n" \
"Signal: \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("Q^{2}_{CC1#pi} (GeV^{2})");
fSettings.SetYTitle("Number of events");
- fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG");
+ // fSettings.SetAllowedTypes("EVT/SHAPE/DIAG", "EVT/SHAPE/DIAG");
fSettings.SetEnuRange(1.0, 10.0);
fSettings.DefineAllowedTargets("D,H");
+ fIsDiag = true;
+ fIsRawEvents = true;
// plot information
fSettings.SetTitle("GGM_CC1ppip_Evt_1DQ2_nu");
fSettings.DefineAllowedSpecies("numu");
fSettings.SetDataInput( FitPar::GetDataBase() + "GGM/CC1pip_on_p/GGM_CC1ppip_Q2_events_bin_edit.txt");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor for shape
fScaleFactor = GetEventHistogram()->Integral("width")*double(1E-38)/double(fNEvents+0.)*(2./1.);
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetPoissonErrors();
SetCovarFromDiagonal();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
}
void GGM_CC1ppip_Evt_1DQ2_nu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(2212) == 0 ||
event->NumFSParticle(211) == 0 ||
event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pp = event->GetHMFSParticle(2212)->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double hadMass = FitUtils::MpPi(Pp, Ppip);
double q2CCpip = -1.0;
// GGM has a M(pi, p) < 1.4 GeV cut imposed only on this channel
if (hadMass < 1400) q2CCpip = -1*(Pnu-Pmu).Mag2()/1.E6;
fXVar = q2CCpip;
return;
}
bool GGM_CC1ppip_Evt_1DQ2_nu::isSignal(FitEvent *event) {
return SignalDef::isCC1pi3Prong(event, 14, 211, 2212, EnuMin, EnuMax);
}
diff --git a/src/MINERvA/MINERvA_CCNpip_XSec_1Dth_nu.cxx b/src/MINERvA/MINERvA_CCNpip_XSec_1Dth_nu.cxx
index 71a712e..9007cdf 100644
--- a/src/MINERvA/MINERvA_CCNpip_XSec_1Dth_nu.cxx
+++ b/src/MINERvA/MINERvA_CCNpip_XSec_1Dth_nu.cxx
@@ -1,262 +1,263 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CCNpip_XSec_1Dth_nu.h"
//********************************************************************
MINERvA_CCNpip_XSec_1Dth_nu::MINERvA_CCNpip_XSec_1Dth_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CCNpip_XSec_1Dth_nu sample. \n" \
"Target: CH \n" \
"Flux: MINERvA Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("#theta_{#pi} (degrees)");
fSettings.SetYTitle("(1/T#Phi) dN_{#pi}/d#theta_{#pi} (cm^{2}/degrees/nucleon)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(1.5, 10.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.DefineAllowedSpecies("numu");
fFullPhaseSpace = !fSettings.Found("name", "_20deg");
fFluxCorrection = fSettings.Found("name", "fluxcorr");
fUpdatedData = !fSettings.Found("name", "2015");
fSettings.SetTitle("MINERvA_CCNpip_XSec_1Dth_nu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents) / TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
// Full Phase Space
if (fFullPhaseSpace) {
// 2016 release data
if (fUpdatedData) {
SetDataFromTextFile( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2016/nu-ccNpi+-xsec-pion-angle.csv");
// MINERvA has the error quoted as a percentage of the cross-section
// Need to make this into an absolute error before we go from correlation matrix -> covariance matrix since it depends on the error in the ith bin
for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) {
fDataHist->SetBinError(i + 1, fDataHist->GetBinContent(i + 1) * (fDataHist->GetBinError(i + 1) / 100.));
}
// This is a correlation matrix! but it's all fixed in SetCovarMatrixFromText
SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2016/nu-ccNpi+-correlation-pion-angle.csv");
// 2015 release data
} else {
// 2015 release allows for shape comparisons
if (fIsShape) {
SetDataFromTextFile( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_th_shape.txt");
SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_th_shape_cov.txt");
} else {
SetDataFromTextFile( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_th.txt");
SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_th_cov.txt");
}
}
// Restricted Phase Space Data
} else {
// 2016 release data unfortunately not released in 20degree forward-going, revert to 2015 data
if (fUpdatedData) {
ERR(FTL) << fName << " has no updated 2016 data for restricted phase space! Using 2015 data." << std::endl;
throw;
}
// Only 2015 20deg data
if (fIsShape) {
SetDataFromTextFile( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_th_20deg_shape.txt");
SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_th_20deg_shape_cov.txt");
} else {
SetDataFromTextFile( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_th_20deg.txt");
SetCorrelationFromTextFile(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCNpip/2015/MINERvA_CCNpi_th_20deg_cov.txt");
}
}
// Scale the MINERvA data to account for the flux difference
// Adjust MINERvA data to flux correction; roughly a 11% normalisation increase in data
// Please change when MINERvA releases new data!
if (fFluxCorrection) {
for (int i = 0; i < fDataHist->GetNbinsX() + 1; i++) {
fDataHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) * 1.11);
}
}
// Make some auxillary helper plots
onePions = (TH1D*)(fDataHist->Clone());
onePions->SetNameTitle((fName + "_1pions").c_str(), (fName + "_1pions" + fPlotTitles).c_str());
SetAutoProcessTH1(onePions, kCMD_Reset, kCMD_Scale, kCMD_Norm);
twoPions = (TH1D*)(fDataHist->Clone());
twoPions->SetNameTitle((fName + "_2pions").c_str(), (fName + "_2pions;" + fPlotTitles).c_str());
SetAutoProcessTH1(twoPions, kCMD_Reset, kCMD_Scale, kCMD_Norm);
threePions = (TH1D*)(fDataHist->Clone());
threePions->SetNameTitle((fName + "_3pions").c_str(), (fName + "_3pions" + fPlotTitles).c_str());
SetAutoProcessTH1(threePions, kCMD_Reset, kCMD_Scale, kCMD_Norm);
morePions = (TH1D*)(fDataHist->Clone());
morePions->SetNameTitle((fName + "_4pions").c_str(), (fName + "_4pions" + fPlotTitles).c_str());
SetAutoProcessTH1(morePions, kCMD_Reset, kCMD_Scale, kCMD_Norm);
// Reserve length 3 for the number of pions
// We fill once per pion found in the event, so can fill multiple times for one event
thVect.reserve(3);
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
// ********************************************
// Fill the event variables
// Here we want to fill the angle for every pion we can find in the event
void MINERvA_CCNpip_XSec_1Dth_nu::FillEventVariables(FitEvent *event) {
// ********************************************
thVect.clear();
if (event->NumFSParticle(211) == 0 && event->NumFSParticle(-211) == 0) return;
if (event->NumFSParticle(13) == 0) return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double hadMass = FitUtils::Wrec(Pnu, Pmu);
if (hadMass < 1800) {
TLorentzVector Ppip;
// Loop over the particle stack
for (unsigned int j = 2; j < event->Npart(); ++j) {
// Only include alive particles
- if (!(event->PartInfo(j))->fIsAlive && (event->PartInfo(j))->fNEUTStatusCode != 0) continue;
+ if ((event->PartInfo(j))->fIsAlive <= 0) continue;
+ if ((event->PartInfo(j))->fNEUTStatusCode != 0) continue;
int PID = (event->PartInfo(j))->fPID;
// Select highest momentum (energy) charged pion
if (abs(PID) == 211) {
Ppip = (event->PartInfo(j))->fP;
double th = (180. / M_PI) * FitUtils::th(Pnu, Ppip);
thVect.push_back(th);
}
}
}
fXVar = 0;
return;
};
//********************************************************************
// The signal definition for MINERvA CCNpi+
// Last bool refers to if we're selecting on the full phase space or not
bool MINERvA_CCNpip_XSec_1Dth_nu::isSignal(FitEvent *event) {
//********************************************************************
return SignalDef::isCCNpip_MINERvA(event, EnuMin, EnuMax, !fFullPhaseSpace);
}
//********************************************************************
// Need to override FillHistograms() here because we fill the histogram N_pion times
void MINERvA_CCNpip_XSec_1Dth_nu::FillHistograms() {
//********************************************************************
if (Signal) {
unsigned int nPions = thVect.size();
// Need to loop over all the pions in the event
for (size_t k = 0; k < nPions; ++k) {
double th = thVect[k];
this->fMCHist->Fill(th, Weight);
this->fMCFine->Fill(th, Weight);
this->fMCStat->Fill(th, 1.0);
if (nPions == 1) {
onePions->Fill(th, Weight);
} else if (nPions == 2) {
twoPions->Fill(th, Weight);
} else if (nPions == 3) {
threePions->Fill(th, Weight);
} else if (nPions > 3) {
morePions->Fill(th, Weight);
}
if (fMCHist_Modes) fMCHist_Modes->Fill(Mode, th, Weight);
// PlotUtils::FillNeutModeArray(fMCHist_PDG, Mode, th, Weight);
}
}
return;
}
//********************************************************************
void MINERvA_CCNpip_XSec_1Dth_nu::Write(std::string drawOpts) {
//********************************************************************
Measurement1D::Write(drawOpts);
// Draw the npions stack
onePions->SetTitle("1#pi");
onePions->SetLineColor(kBlack);
//onePions->SetFillStyle(0);
onePions->SetFillColor(onePions->GetLineColor());
twoPions->SetTitle("2#pi");
twoPions->SetLineColor(kRed);
//twoPions->SetFillStyle(0);
twoPions->SetFillColor(twoPions->GetLineColor());
threePions->SetTitle("3#pi");
threePions->SetLineColor(kGreen);
//threePions->SetFillStyle(0);
threePions->SetFillColor(threePions->GetLineColor());
morePions->SetTitle(">3#pi");
morePions->SetLineColor(kBlue);
//morePions->SetFillStyle(0);
morePions->SetFillColor(morePions->GetLineColor());
THStack pionStack = THStack((fName + "_pionStack").c_str(), (fName + "_pionStack").c_str());
pionStack.Add(onePions);
pionStack.Add(twoPions);
pionStack.Add(threePions);
pionStack.Add(morePions);
pionStack.Write();
return;
}
diff --git a/src/MINERvA/MINERvA_CCQE_XSec_1DQ2_antinu.cxx b/src/MINERvA/MINERvA_CCQE_XSec_1DQ2_antinu.cxx
index 7adf843..8ce4ee0 100644
--- a/src/MINERvA/MINERvA_CCQE_XSec_1DQ2_antinu.cxx
+++ b/src/MINERvA/MINERvA_CCQE_XSec_1DQ2_antinu.cxx
@@ -1,146 +1,140 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "MINERvA_SignalDef.h"
#include "MINERvA_CCQE_XSec_1DQ2_antinu.h"
//********************************************************************
MINERvA_CCQE_XSec_1DQ2_antinu::MINERvA_CCQE_XSec_1DQ2_antinu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MINERvA_CCQE_XSec_1DQ2_antinu sample. \n" \
"Target: CH \n" \
"Flux: MINERvA Forward Horn Current Numubar \n" \
"Signal: True CCQE/2p2h defined at the vertex level \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("Q^{2}_{QE} (GeV^{2})");
fSettings.SetYTitle("d#sigma/dQ_{QE}^{2} (cm^{2}/GeV^{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(1.5, 10.0);
fSettings.DefineAllowedTargets("C,H");
isFluxFix = !fSettings.Found("name", "_oldflux");
fullphasespace = !fSettings.Found("name", "_20deg");
// CCQELike plot information
fSettings.SetTitle("MINERvA_CCQE_XSec_1DQ2_antinu");
std::string basedir = FitPar::GetDataBase() + "/MINERvA/CCQE/";
std::string datafilename = "";
std::string covarfilename = "";
// Full Phase Space
if (fullphasespace) {
if (isFluxFix) {
if (fIsShape) {
ERR(WRN) << "SHAPE likelihood comparison not available for MINERvA "
<< "datasets with fixed flux information. NUISANCE will scale MC to match "
<< "data normalization but full covariance will be used. " << std::endl;
}
datafilename = "Q2QE_numubar_data_fluxfix.txt";
- covarfilename = "Q2QE_numubar_covar_fluxfix.txt";
+ covarfilename = "Q2QE_numubar_covar_fluxfix.txt"; // Correlation Matrix
} else {
if (fIsShape) {
datafilename = "Q2QE_numubar_data_SHAPE-extracted.txt";
- covarfilename = "Q2QE_numubar_covar_SHAPE-extracted.txt";
+ covarfilename = "Q2QE_numubar_covar_SHAPE-extracted.txt"; // correlation
} else {
datafilename = "Q2QE_numubar_data.txt";
- covarfilename = "Q2QE_numubar_covar.txt";
+ covarfilename = "Q2QE_numubar_covar.txt"; // Correlation
}
}
// Restricted Phase Space
} else {
if (isFluxFix) {
if (fIsShape) {
ERR(WRN) << "SHAPE likelihood comparison not available for MINERvA "
<< "datasets with fixed flux information. NUISANCE will scale MC to match "
<< "data normalization but full covariance will be used. " << std::endl;
}
datafilename = "20deg_Q2QE_numubar_data_fluxfix.txt";
- covarfilename = "20deg_Q2QE_numubar_covar_fluxfix.txt";
+ covarfilename = "20deg_Q2QE_numubar_covar_fluxfix.txt"; // Correlation
} else {
if (fIsShape) {
datafilename = "20deg_Q2QE_numubar_data_SHAPE-extracted.txt";
- covarfilename = "20deg_Q2QE_numubar_covar_SHAPE-extracted.txt";
+ covarfilename = "20deg_Q2QE_numubar_covar_SHAPE-extracted.txt"; // Correlation
} else {
datafilename = "20deg_Q2QE_numubar_data.txt";
- covarfilename = "20deg_Q2QE_numubar_covar.txt";
+ covarfilename = "20deg_Q2QE_numubar_covar.txt"; // Correlation
}
}
}
fSettings.SetDataInput( basedir + datafilename );
fSettings.SetCovarInput( basedir + covarfilename );
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = (GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) * 13. / 7. / TotalIntegratedFlux();
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
-
- if (!isFluxFix or !fullphasespace) {
- SetCorrelationFromTextFile( fSettings.GetCovarInput() );
- ScaleCovar(1E76);
- } else {
- SetCovarFromTextFile( fSettings.GetCovarInput() );
- }
-
+ SetCorrelationFromTextFile( fSettings.GetCovarInput() );
+
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MINERvA_CCQE_XSec_1DQ2_antinu::FillEventVariables(FitEvent *event) {
//********************************************************************
if (event->NumFSParticle(-13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Pmu = event->GetHMFSParticle(-13)->fP;
double ThetaMu = Pnu.Vect().Angle(Pmu.Vect());
double q2qe = FitUtils::Q2QErec(Pmu, cos(ThetaMu), 30., true);
fXVar = q2qe;
return;
}
//********************************************************************
bool MINERvA_CCQE_XSec_1DQ2_antinu::isSignal(FitEvent *event) {
//*******************************************************************
return SignalDef::isCCQEnumubar_MINERvA(event, EnuMin, EnuMax, fullphasespace);
}
diff --git a/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1DEnu_nu.cxx b/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1DEnu_nu.cxx
index e55e2f6..dfdd241 100644
--- a/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1DEnu_nu.cxx
+++ b/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1DEnu_nu.cxx
@@ -1,94 +1,94 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "MiniBooNE_CC1pi0_XSec_1DEnu_nu.h"
//********************************************************************
MiniBooNE_CC1pi0_XSec_1DEnu_nu::MiniBooNE_CC1pi0_XSec_1DEnu_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MiniBooNE_CC1pi0_XSec_1DEnu_nu sample. \n" \
"Target: CH \n" \
"Flux: MiniBooNE Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 muon, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("E_{#nu} (GeV)");
fSettings.SetYTitle("#sigma(E_{#nu}) (cm^{2}/CH_{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(0.5, 2.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("MiniBooNE_CC1pi0_XSec_1DEnu_nu");
fSettings.SetDataInput( FitPar::GetDataBase() + "MiniBooNE/CC1pi0/totalxsec_edit.txt" );
fSettings.SetCovarInput( FitPar::GetDataBase() + "MiniBooNE/CC1pi0/totalxsec_covar.txt" );
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width")*double(1E-38)/double(fNEvents)*(14.08);
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
- SetCovarFromTextFile( fSettings.GetCovarInput() );
+ SetCorrelationFromTextFile( fSettings.GetCovarInput() );
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void MiniBooNE_CC1pi0_XSec_1DEnu_nu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(111) == 0 ||
event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu =event->GetNeutrinoIn()->fP;
TLorentzVector Ppi0 = event->GetHMFSParticle(111)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double Enu = FitUtils::EnuCC1pi0rec(Pnu, Pmu, Ppi0);
fXVar = Enu;
return;
};
// **************************************************
// MiniBooNE CC1pi0 signal definition
//
// The signal definition is:
// Exactly one negative muon
// Exactly one pi0
// No additional mesons
// Any number of nucleons
//
// Does a few clever cuts on the likelihood to reduce CCQE contamination by
// looking at "fuzziness" of the ring; CCQE events are sharp, CC1pi0 are fuzzy
// (because of the pi0->2 gamma collinearity)
bool MiniBooNE_CC1pi0_XSec_1DEnu_nu::isSignal(FitEvent *event) {
return SignalDef::isCC1pi(event, 14, 111, EnuMin, EnuMax);
}
diff --git a/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1DQ2_nu.cxx b/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1DQ2_nu.cxx
index 17d43dc..4dfaf38 100644
--- a/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1DQ2_nu.cxx
+++ b/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1DQ2_nu.cxx
@@ -1,84 +1,84 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "MiniBooNE_CC1pi0_XSec_1DQ2_nu.h"
//********************************************************************
MiniBooNE_CC1pi0_XSec_1DQ2_nu::MiniBooNE_CC1pi0_XSec_1DQ2_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MiniBooNE_CC1pi0_XSec_1DQ2_nu sample. \n" \
"Target: CH \n" \
"Flux: MiniBooNE Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 muon, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("Q^{2}_{CC#pi} (GeV^{2})");
fSettings.SetYTitle("d#sigma/dQ_{CC#pi}^{2} (cm^{2}/GeV^{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(0.5, 2.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("MiniBooNE_CC1pi0_XSec_1DQ2_nu");
fSettings.SetDataInput( FitPar::GetDataBase() + "MiniBooNE/CC1pi0/dxsecdq2_edit.txt" );
fSettings.SetCovarInput( FitPar::GetDataBase() + "MiniBooNE/CC1pi0/dxsecdq2_covar.txt" );
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width")*double(1E-38)/double(fNEvents)*(14.08)/TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
- SetCovarFromTextFile( fSettings.GetCovarInput() );
+ SetCorrelationFromTextFile( fSettings.GetCovarInput() );
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void MiniBooNE_CC1pi0_XSec_1DQ2_nu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(111) == 0 ||
event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Ppi0 = event->GetHMFSParticle(111)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double q2CCpi0 = FitUtils::Q2CC1pi0rec(Pnu, Pmu, Ppi0);
fXVar = q2CCpi0;
return;
};
//********************************************************************
bool MiniBooNE_CC1pi0_XSec_1DQ2_nu::isSignal(FitEvent *event) {
//********************************************************************
return SignalDef::isCC1pi(event, 14, 111, EnuMin, EnuMax);
}
diff --git a/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1DTu_nu.cxx b/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1DTu_nu.cxx
index bac106d..3245928 100644
--- a/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1DTu_nu.cxx
+++ b/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1DTu_nu.cxx
@@ -1,83 +1,83 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "MiniBooNE_CC1pi0_XSec_1DTu_nu.h"
//********************************************************************
MiniBooNE_CC1pi0_XSec_1DTu_nu::MiniBooNE_CC1pi0_XSec_1DTu_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MiniBooNE_CC1pi0_XSec_1DTu_nu sample. \n" \
"Target: CH \n" \
"Flux: MiniBooNE Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 muon, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("T_{#mu} (GeV)");
fSettings.SetYTitle("d#sigma/dE_{#mu} (cm^{2}/GeV^{2}/CH_{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(0.5, 2.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("MiniBooNE_CC1pi0_XSec_1DTu_nu");
fSettings.SetDataInput( FitPar::GetDataBase() + "MiniBooNE/CC1pi0/dxsecdemu_edit.txt" );
fSettings.SetCovarInput( FitPar::GetDataBase() + "MiniBooNE/CC1pi0/dxsecdemu_covar.txt" );
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width")*double(1E-38)/double(fNEvents)*(14.08)/TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
- SetCovarFromTextFile( fSettings.GetCovarInput() );
+ SetCorrelationFromTextFile( fSettings.GetCovarInput() );
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void MiniBooNE_CC1pi0_XSec_1DTu_nu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(111) == 0 ||
event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Ppi0 = event->GetHMFSParticle(111)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double TuCCpi0 = FitUtils::T(Pmu);
fXVar = TuCCpi0;
return;
};
//********************************************************************
bool MiniBooNE_CC1pi0_XSec_1DTu_nu::isSignal(FitEvent *event) {
//********************************************************************
return SignalDef::isCC1pi(event, 14, 111, EnuMin, EnuMax);
}
diff --git a/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1Dcosmu_nu.cxx b/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1Dcosmu_nu.cxx
index ca3e1a6..7ed184e 100644
--- a/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1Dcosmu_nu.cxx
+++ b/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1Dcosmu_nu.cxx
@@ -1,85 +1,85 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "MiniBooNE_CC1pi0_XSec_1Dcosmu_nu.h"
//********************************************************************
MiniBooNE_CC1pi0_XSec_1Dcosmu_nu::MiniBooNE_CC1pi0_XSec_1Dcosmu_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MiniBooNE_CC1pi0_XSec_1Dcosmu_nu sample. \n" \
"Target: CH \n" \
"Flux: MiniBooNE Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 muon, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("cos#theta_{#mu}");
fSettings.SetYTitle("d#sigma/dcos#theta_{#mu} (cm^{2}/CH_{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(0.5, 2.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("MiniBooNE_CC1pi0_XSec_1Dcosmu_nu");
fSettings.SetDataInput( FitPar::GetDataBase() + "MiniBooNE/CC1pi0/dxsecdcosmu_edit.txt" );
fSettings.SetCovarInput( FitPar::GetDataBase() + "MiniBooNE/CC1pi0/dxsecdcosmu_covar.txt" );
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width")*double(1E-38)/double(fNEvents)*(14.08)/TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
- SetCovarFromTextFile( fSettings.GetCovarInput() );
+ SetCorrelationFromTextFile( fSettings.GetCovarInput() );
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void MiniBooNE_CC1pi0_XSec_1Dcosmu_nu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(111) == 0 ||
event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Ppi0 = event->GetHMFSParticle(111)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
// No W cut on MiniBooNE CC1pi+
double CosMu = cos(FitUtils::th(Pnu, Pmu));
fXVar = CosMu;
return;
};
//********************************************************************
bool MiniBooNE_CC1pi0_XSec_1Dcosmu_nu::isSignal(FitEvent *event) {
//********************************************************************
return SignalDef::isCC1pi(event, 14, 111, EnuMin, EnuMax);
}
diff --git a/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1Dcospi0_nu.cxx b/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1Dcospi0_nu.cxx
index 81d61ca..0758791 100644
--- a/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1Dcospi0_nu.cxx
+++ b/src/MiniBooNE/MiniBooNE_CC1pi0_XSec_1Dcospi0_nu.cxx
@@ -1,84 +1,84 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "MiniBooNE_CC1pi0_XSec_1Dcospi0_nu.h"
//********************************************************************
MiniBooNE_CC1pi0_XSec_1Dcospi0_nu::MiniBooNE_CC1pi0_XSec_1Dcospi0_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MiniBooNE_CC1pi0_XSec_1Dcospi0_nu sample. \n" \
"Target: CH \n" \
"Flux: MiniBooNE Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 muon, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("cos#theta_{#pi^{0}}");
fSettings.SetYTitle("d#sigma/dcos#theta_{#pi^{0}} (cm^{2}/CH_{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/FULL,DIAG/NORM/MASK", "FIX/FULL");
fSettings.SetEnuRange(0.5, 2.0);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("MiniBooNE_CC1pi0_XSec_1Dcospi0_nu");
fSettings.SetDataInput( FitPar::GetDataBase() + "MiniBooNE/CC1pi0/dxsecdcospi_edit.txt" );
fSettings.SetCovarInput( FitPar::GetDataBase() + "MiniBooNE/CC1pi0/dxsecdcospi_covar.txt" );
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width")*double(1E-38)/double(fNEvents)*(14.08)/TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
- SetCovarFromTextFile( fSettings.GetCovarInput() );
+ SetCorrelationFromTextFile( fSettings.GetCovarInput() );
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
void MiniBooNE_CC1pi0_XSec_1Dcospi0_nu::FillEventVariables(FitEvent *event) {
if (event->NumFSParticle(111) == 0 ||
event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu =event->GetNeutrinoIn()->fP;
TLorentzVector Ppi0 = event->GetHMFSParticle(111)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
double CosPi0 = cos(FitUtils::th(Pnu, Ppi0));
fXVar = CosPi0;
return;
};
//********************************************************************
bool MiniBooNE_CC1pi0_XSec_1Dcospi0_nu::isSignal(FitEvent *event) {
//********************************************************************
return SignalDef::isCC1pi(event, 14, 111, EnuMin, EnuMax);
}
diff --git a/src/MiniBooNE/MiniBooNE_CC1pip_XSec_1DTpi_nu.cxx b/src/MiniBooNE/MiniBooNE_CC1pip_XSec_1DTpi_nu.cxx
index 8524b8d..1becf55 100644
--- a/src/MiniBooNE/MiniBooNE_CC1pip_XSec_1DTpi_nu.cxx
+++ b/src/MiniBooNE/MiniBooNE_CC1pip_XSec_1DTpi_nu.cxx
@@ -1,90 +1,90 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "MiniBooNE_CC1pip_XSec_1DTpi_nu.h"
//********************************************************************
MiniBooNE_CC1pip_XSec_1DTpi_nu::MiniBooNE_CC1pip_XSec_1DTpi_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MiniBooNE_CC1pip_XSec_1DTpi_nu sample. \n" \
"Target: CH \n" \
"Flux: MiniBooNE Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 muon, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("T_{#pi} (MeV)");
fSettings.SetYTitle("d#sigma/dT_{#pi^{+}}} (cm^{2}/MeV/CH_{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG/NORM/MASK", "FIX/DIAG");
- fSettings.SetEnuRange(0.0, 3.0);
+ fSettings.SetEnuRange(0.0, 100.0);
fSettings.SetNormError(0.107);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("MiniBooNE_CC1pip_XSec_1DTpi_nu");
fSettings.SetDataInput( FitPar::GetDataBase() + "MiniBooNE/CC1pip/ccpipXSec_KEpi.txt" );
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width")*double(1E-38)/double(fNEvents)*(14.08)/TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetCovarFromDiagonal();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MiniBooNE_CC1pip_XSec_1DTpi_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
if (event->NumFSParticle(211) == 0 ||
event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
// No W cut on MiniBooNE data from publication
// WARNING: DRAGONS LAY HERE! Mike Wilking's thesis might not have this. Beware that the publication says W < 1.35 GeV, but this is "efficiency corrected", i.e. FILLED WITH MONTE-CARLO!!!!!!!! AAAAAH
//double hadMass = FitUtils::Wrec(Pnu, Pmu, Ppip);
double Tpi = FitUtils::T(Ppip)*1000.;
fXVar = Tpi;
return;
};
bool MiniBooNE_CC1pip_XSec_1DTpi_nu::isSignal(FitEvent *event) {
return SignalDef::isCC1pi(event, 14, 211, EnuMin, EnuMax);
}
diff --git a/src/MiniBooNE/MiniBooNE_CC1pip_XSec_1DTu_nu.cxx b/src/MiniBooNE/MiniBooNE_CC1pip_XSec_1DTu_nu.cxx
index d9fb93a..2979051 100644
--- a/src/MiniBooNE/MiniBooNE_CC1pip_XSec_1DTu_nu.cxx
+++ b/src/MiniBooNE/MiniBooNE_CC1pip_XSec_1DTu_nu.cxx
@@ -1,84 +1,84 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "MiniBooNE_CC1pip_XSec_1DTu_nu.h"
//********************************************************************
MiniBooNE_CC1pip_XSec_1DTu_nu::MiniBooNE_CC1pip_XSec_1DTu_nu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MiniBooNE_CC1pip_XSec_1DTu_nu sample. \n" \
"Target: CH \n" \
"Flux: MiniBooNE Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 muon, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("T_{#mu} (MeV)");
fSettings.SetYTitle("d#sigma/dT_{#mu} (cm^{2}/MeV/CH_{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG/NORM/MASK", "FIX/DIAG");
- fSettings.SetEnuRange(0.0, 3.0);
+ fSettings.SetEnuRange(0.0, 100.0);
fSettings.SetNormError(0.107);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("MiniBooNE_CC1pip_XSec_1DTu_nu");
fSettings.SetDataInput( FitPar::GetDataBase() + "MiniBooNE/CC1pip/ccpipXSec_KEmu.txt" );
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = GetEventHistogram()->Integral("width")*double(1E-38)/double(fNEvents)*(14.08)/TotalIntegratedFlux("width");
// Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetCovarFromDiagonal();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MiniBooNE_CC1pip_XSec_1DTu_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
if (event->NumFSParticle(211) == 0 ||
event->NumFSParticle(13) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
TLorentzVector Ppip = event->GetHMFSParticle(211)->fP;
TLorentzVector Pmu = event->GetHMFSParticle(13)->fP;
// No W cut on MiniBooNE CC1pi+
double Tmu = FitUtils::T(Pmu)*1000.;
fXVar = Tmu;
return;
};
bool MiniBooNE_CC1pip_XSec_1DTu_nu::isSignal(FitEvent *event) {
return SignalDef::isCC1pi(event, 14, 211, EnuMin, EnuMax);
}
diff --git a/src/MiniBooNE/MiniBooNE_CCQE_XSec_1DQ2_antinu.cxx b/src/MiniBooNE/MiniBooNE_CCQE_XSec_1DQ2_antinu.cxx
index a1cfcae..8483e70 100644
--- a/src/MiniBooNE/MiniBooNE_CCQE_XSec_1DQ2_antinu.cxx
+++ b/src/MiniBooNE/MiniBooNE_CCQE_XSec_1DQ2_antinu.cxx
@@ -1,195 +1,196 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "MiniBooNE_CCQE_XSec_1DQ2_antinu.h"
#include <csignal>
MiniBooNE_CCQE_XSec_1DQ2_antinu::MiniBooNE_CCQE_XSec_1DQ2_antinu(nuiskey samplekey) {
// Sample overview
std::string descrip = "MiniBooNE CCQE/CC0pi sample. \n" \
"Target: CH2.08 \n" \
"Flux: CCQE = Forward Horn Current numu \n" \
" CC0pi = Forward Horn Current numu+numub \n" \
"Signal: CCQE = True CCQE + True 2p2h (Mode == 1 or 2) \n" \
" CC0pi = Events with 1 mu+/mu-, N nucleons, 0 other";
// 1. Initalise sample Settings ---------------------------------------
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("Q^{2}_{QE} (GeV^{2})");
fSettings.SetYTitle("d#sigma/dQ_{QE}^{2} (cm^{2}/GeV^{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 3.0);
fSettings.SetSuggestedFlux( FitPar::GetDataBase() + "/MiniBooNE/ccqe/mb_ccqe_flux.root");
// Define input data information
fSettings.FoundFill("name", "CCQELike", fCCQElike, true);
fSettings.FoundFill("name", "CTarg", fUseCorrectedCTarget, true);
if (fCCQElike && fUseCorrectedCTarget) {
ERR(FTL) << "Sample: MiniBooNE_CCQE_XSec_1DQ2_antinu cannot run in both "
"QELike and C-Target mode. You're welcome to add the data set."
<< std::endl;
throw;
}
if (fCCQElike) {
// CCQELike plot information
fSettings.SetTitle("MiniBooNE #nu_#mu CCQE on CH");
fSettings.SetDataInput( FitPar::GetDataBase() + "/MiniBooNE/anti-ccqe/asqq_like.txt" );
fSettings.SetCovarInput( FitPar::GetDataBase() + "/MiniBooNE/anti-ccqe/asqq_diagcovar" );
fSettings.SetDefault( "ccqelikebkg_input", FitPar::GetDataBase() + "/MiniBooNE/anti-ccqe/asqq_bkg_ccqe.txt" );
fSettings.SetDefault( "ccpimbkg_input", FitPar::GetDataBase() + "/MiniBooNE/anti-ccqe/asqq_bkg_ccpim.txt" );
fSettings.SetHasExtraHistograms(true);
fSettings.DefineAllowedSpecies("numu,numub");
fSettings.DefineAllowedTargets("C,H");
} else if (!fUseCorrectedCTarget) {
// CCQE Plot Information
fSettings.SetTitle("MiniBooNE #nu_#mu CC0#pi on CH");
fSettings.SetDataInput( FitPar::GetDataBase() + "/MiniBooNE/anti-ccqe/asqq_con.txt" );
fSettings.SetCovarInput( FitPar::GetDataBase() + "/MiniBooNE/anti-ccqe/asqq_diagcovar" );
fSettings.DefineAllowedSpecies("numu");
fSettings.DefineAllowedTargets("C,H");
} else {
// CCQE Corrected Target Plot Information
fSettings.SetTitle("MiniBooNE #nu_#mu CC0#pi on C");
fSettings.SetDataInput( FitPar::GetDataBase() + "/MiniBooNE/anti-ccqe/asqq_con_ctarget.txt" );
fSettings.SetCovarInput( FitPar::GetDataBase() + "/MiniBooNE/anti-ccqe/asqq_diagcovar" );
fSettings.DefineAllowedSpecies("numu");
fSettings.DefineAllowedTargets("C");
}
FinaliseSampleSettings();
// 2. Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
// Multiply by 14.08/6.0 to get per neutron
double NNucPerNTarg = fUseCorrectedCTarget ? 12.0 / 6.0 : 14.08 / 8.0;
fScaleFactor = ((GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) *
NNucPerNTarg / TotalIntegratedFlux());
// 3. Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetCovarFromDiagonal();
///
/// If CCQELike is used an additional the CCQELike BKG is used and a PDG
/// Histogram is saved
if (fCCQElike) {
// CCQELike Data
fDataHist_CCQELIKE = PlotUtils::GetTH1DFromFile( fSettings.GetS("ccqelikebkg_input"),
fSettings.GetName() + "_CCQELIKE_data" );
fDataHist_CCQELIKE->SetNameTitle( (fSettings.Name() + "_CCQELIKE_BKG").c_str(),
("MiniBooNE #nu_#mu CCQE-Like Backgrounds" + fSettings.PlotTitles()).c_str() );
+ fDataHist->Add(fDataHist_CCQELIKE);
SetAutoProcessTH1(fDataHist_CCQELIKE, kCMD_Write);
// CCQELike MC
fMCHist_CCQELIKE = new NuNuBarTrueModeStack( fSettings.Name() + "_CCQELIKE_MC",
"CCQE-like MC" + fSettings.PlotTitles(),
fDataHist_CCQELIKE );
SetAutoProcessTH1(fMCHist_CCQELIKE);
// Data CCRES
fDataHist_CCPIM = PlotUtils::GetTH1DFromFile( fSettings.GetS("ccpimbkg_input"),
fSettings.GetName() + "_CCPIM_BKG_data" );
fDataHist_CCPIM->SetNameTitle( (fSettings.Name() + "_CCPIM_data").c_str(),
("MiniBooNE #nu_#mu CCQE-Like Backgrounds" + fSettings.PlotTitles()).c_str() );
SetAutoProcessTH1(fDataHist_CCQELIKE, kCMD_Write);
// MC CCRES
fMCHist_CCPIM = new NuNuBarTrueModeStack(fSettings.Name() + "_CCPIM_BKG_MC", "CCQE-like BKG CC-RES" + fSettings.PlotTitles(), fDataHist_CCPIM);
SetAutoProcessTH1(fMCHist_CCPIM);
// Make NON CCPIM
fDataHist_NONCCPIM = (TH1D *)fDataHist_CCQELIKE->Clone();
fDataHist_NONCCPIM->SetNameTitle((fName + "_data_NONCCPIM").c_str(),
(fName + "_data_NONCCPIM").c_str());
fDataHist_NONCCPIM->Add( fDataHist_CCPIM, -1.0 );
SetAutoProcessTH1(fDataHist_NONCCPIM, kCMD_Write);
fMCHist_NONCCPIM = new NuNuBarTrueModeStack(
fSettings.Name() + "_NONCCPIM_BKG",
"CCQE-like BKG CC-NonRES" + fSettings.PlotTitles(),
fDataHist_NONCCPIM);
SetAutoProcessTH1(fMCHist_NONCCPIM);
}
FinaliseMeasurement();
};
void MiniBooNE_CCQE_XSec_1DQ2_antinu::FillEventVariables(FitEvent * event) {
if (event->NumFSParticle(PhysConst::pdg_muons) == 0) return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
// The highest momentum mu+/mu-. The isSignal definition should make sure we
// only
// accept events we want, so no need to do an additional check here.
TLorentzVector Pmu = event->GetHMFSParticle(PhysConst::pdg_muons)->fP;
// Set X Variables
fXVar = FitUtils::Q2QErec(Pmu, cos(Pnu.Vect().Angle(Pmu.Vect())), 30., false);
fPDGnu = event->PDGnu();
return;
};
bool MiniBooNE_CCQE_XSec_1DQ2_antinu::isSignal(FitEvent * event) {
// If CC0pi, include both charges
if (fCCQElike) {
if (SignalDef::isCC0pi(event, 14, EnuMin, EnuMax) ||
SignalDef::isCC0pi(event, -14, EnuMin, EnuMax)) {
return true;
}
} else {
if (SignalDef::isCCQELike(event, -14, EnuMin, EnuMax)) return true;
}
return false;
};
void MiniBooNE_CCQE_XSec_1DQ2_antinu::FillExtraHistograms(MeasurementVariableBox* vars, double weight) {
// No Extra Hists if not ccqelike
if (!fCCQElike or !Signal) return;
// Fill Stacks
if (Mode != -1 and Mode != -2) {
if (fabs(Mode) == 11 or fabs(Mode) == 12 or fabs(Mode == 13)) {
fMCHist_CCPIM->Fill(fPDGnu, Mode, fXVar, weight);
} else {
fMCHist_NONCCPIM->Fill(fPDGnu, Mode, fXVar, weight);
}
}
fMCHist_CCQELIKE->Fill(fPDGnu, Mode, fXVar, weight);
}
diff --git a/src/MiniBooNE/MiniBooNE_CCQE_XSec_1DQ2_nu.cxx b/src/MiniBooNE/MiniBooNE_CCQE_XSec_1DQ2_nu.cxx
index 052ce2b..ecb4edf 100644
--- a/src/MiniBooNE/MiniBooNE_CCQE_XSec_1DQ2_nu.cxx
+++ b/src/MiniBooNE/MiniBooNE_CCQE_XSec_1DQ2_nu.cxx
@@ -1,182 +1,183 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "MiniBooNE_CCQE_XSec_1DQ2_nu.h"
//********************************************************************
MiniBooNE_CCQE_XSec_1DQ2_nu::MiniBooNE_CCQE_XSec_1DQ2_nu(nuiskey samplekey) {
//********************************************************************
// Initalise sample Settings ---------------------------------------
fSettings = LoadSampleSettings(samplekey);
fSettings.FoundFill("name", "CCQELike", ccqelike, true);
// Multiple constructors for similar samples
if (ccqelike) Setup_MiniBooNE_CCQELike_XSec_1DQ2_nu();
else Setup_MiniBooNE_CCQE_XSec_1DQ2_nu();
// Final Check for all requirements, necessary to setup all extra plots.
FinaliseMeasurement();
}
//********************************************************************
void MiniBooNE_CCQE_XSec_1DQ2_nu::Setup_MiniBooNE_CCQE_XSec_1DQ2_nu() {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MiniBooNE CCQE sample. \n" \
"Target: CH2.08 \n" \
"Flux: CCQE = Forward Horn Current numu \n" \
"Signal: CCQE = True CCQE + True 2p2h (Mode == 1 or 2) \n";
// Setup common settings
fSettings.SetDescription(descrip);
fSettings.SetXTitle("Q^{2}_{QE} (GeV^{2})");
fSettings.SetYTitle("d#sigma/dQ_{QE}^{2} (cm^{2}/GeV^{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 3.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.SetSuggestedFlux( FitPar::GetDataBase() + "/MiniBooNE/ccqe/mb_ccqe_flux.root");
// CCQE Plot Information
fSettings.SetTitle("MiniBooNE #nu_#mu CC0#pi");
fSettings.SetDataInput( FitPar::GetDataBase() + "/MiniBooNE/ccqe/asqq_con.txt" );
fSettings.DefineAllowedSpecies("numu");
FinaliseSampleSettings();
// 2. Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
// Multiply by 14.08/6.0 to get per neutron
fScaleFactor = ((GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) *
(14.08 / 6.0) / TotalIntegratedFlux());
// 3. Plot Setup -------------------------------------------------------
SetDataFromTextFile( fSettings.GetDataInput() );
SetCovarFromDiagonal();
};
//********************************************************************
void MiniBooNE_CCQE_XSec_1DQ2_nu::Setup_MiniBooNE_CCQELike_XSec_1DQ2_nu() {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MiniBooNE CC0pi sample. \n" \
"Target: CH2.08 \n" \
"Flux: MiniBooNE Forward Horn Current numu \n" \
"Signal: Any event with 1 muon, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings.SetDescription(descrip);
fSettings.SetXTitle("Q^{2}_{QE} (GeV^{2})");
fSettings.SetYTitle("d#sigma/dQ_{QE}^{2} (cm^{2}/GeV^{2})");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 3.0);
fSettings.DefineAllowedTargets("C,H");
fSettings.SetSuggestedFlux( FitPar::GetDataBase() + "/MiniBooNE/ccqe/mb_ccqe_flux.root");
// CCQELike plot information
fSettings.SetTitle("MiniBooNE #nu_#mu CCQE");
fSettings.SetDataInput( FitPar::GetDataBase() + "/MiniBooNE/ccqe/asqq_like.txt" );
fSettings.SetCovarInput( FitPar::GetDataBase() + "/MiniBooNE/ccqe/asqq_diagcovar" );
fSettings.SetDefault( "ccqelikebkg_input", FitPar::GetDataBase() + "/MiniBooNE/ccqe/asqq_bkg.txt" );
fSettings.SetHasExtraHistograms(true);
fSettings.DefineAllowedSpecies("numu,numub");
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
// Multiply by 14.08/6.0 to get per neutron
fScaleFactor = ((GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) *
(14.08 / 6.0) / TotalIntegratedFlux());
// Plot Setup -------------------------------------------------------
fDataHist = PlotUtils::GetTH1DFromFile( fSettings.GetDataInput(), fSettings.GetName() );
fDataHist->SetTitle( (fSettings.Title() + fSettings.PlotTitles()).c_str() );
// Make Data CCQELike BKG
fDataHist_CCQELIKE = PlotUtils::GetTH1DFromFile( fSettings.GetS("ccqelikebkg_input"),
fSettings.GetName() + "_CCQELIKEBKG_data" );
fDataHist_CCQELIKE->SetNameTitle( (fSettings.Name() + "_CCQELIKE_BKG_data").c_str(),
("MiniBooNE #nu_#mu CCQE-Like Backgrounds" + fSettings.PlotTitles()).c_str() );
+ // fDataHist->Add(fDataHist_CCQELIKE);
SetAutoProcessTH1(fDataHist_CCQELIKE, kCMD_Write);
// Make MC Clone
fMCHist_CCQELIKE = new TrueModeStack( fSettings.Name() + "_CCQELIKE_BKG_MC",
"CCQE-like BKG MC" + fSettings.PlotTitles(),
fDataHist_CCQELIKE);
SetAutoProcessTH1(fMCHist_CCQELIKE);
};
//********************************************************************
void MiniBooNE_CCQE_XSec_1DQ2_nu::FillEventVariables(FitEvent *event) {
//********************************************************************
if (event->NumFSParticle(PhysConst::pdg_muons) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
// The highest momentum mu+/mu-. The isSignal definition should make sure we only
// accept events we want, so no need to do an additional check here.
TLorentzVector Pmu = event->GetHMFSParticle(PhysConst::pdg_muons)->fP;
q2qe = FitUtils::Q2QErec(Pmu, cos(Pnu.Vect().Angle(Pmu.Vect())), 34., false);
// Set X Variables
fXVar = q2qe;
return;
};
//********************************************************************
bool MiniBooNE_CCQE_XSec_1DQ2_nu::isSignal(FitEvent *event) {
//********************************************************************
// If CC0pi, include both charges
if (ccqelike) {
if (SignalDef::isCC0pi(event, 14, EnuMin, EnuMax) ||
SignalDef::isCC0pi(event, -14, EnuMin, EnuMax))
return true;
} else {
if (SignalDef::isCCQELike(event, 14, EnuMin, EnuMax))
return true;
}
return false;
};
void MiniBooNE_CCQE_XSec_1DQ2_nu::FillExtraHistograms(MeasurementVariableBox* vars, double weight) {
// No Extra Hists if not ccqelike
if (!ccqelike) return;
if ((Mode != 1 and Mode != 2) and (Signal)) {
fMCHist_CCQELIKE->Fill(Mode, fXVar, weight);
}
return;
}
diff --git a/src/MiniBooNE/MiniBooNE_CCQE_XSec_2DTcos_antinu.cxx b/src/MiniBooNE/MiniBooNE_CCQE_XSec_2DTcos_antinu.cxx
index d7fc791..4701e28 100644
--- a/src/MiniBooNE/MiniBooNE_CCQE_XSec_2DTcos_antinu.cxx
+++ b/src/MiniBooNE/MiniBooNE_CCQE_XSec_2DTcos_antinu.cxx
@@ -1,126 +1,126 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "MiniBooNE_CCQE_XSec_2DTcos_antinu.h"
//********************************************************************
MiniBooNE_CCQE_XSec_2DTcos_antinu::MiniBooNE_CCQE_XSec_2DTcos_antinu(nuiskey samplekey) {
//********************************************************************
// Sample overview ---------------------------------------------------
std::string descrip = "MiniBooNE_CCQE_XSec_2DTcos_antinu sample. \n" \
"Target: CH \n" \
"Flux: MiniBooNE Forward Horn Current nue + nuebar \n" \
"Signal: Any event with 1 electron, any nucleons, and no other FS particles \n";
// Setup common settings
fSettings = LoadSampleSettings(samplekey);
fSettings.SetDescription(descrip);
fSettings.SetXTitle("T_{#mu} (GeV)");
fSettings.SetYTitle("cos#theta_{#mu}");
fSettings.SetZTitle("d^{2}#sigma/dT_{#mu}dcos#theta_{#mu} (cm^{2}/GeV)");
fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG/NORM/MASK", "FIX/DIAG");
fSettings.SetEnuRange(0.0, 3.0);
fSettings.SetNormError(0.130);
fSettings.DefineAllowedTargets("C,H");
// CCQELike plot information
fSettings.SetTitle("MiniBooNE_CCQE_XSec_2DTcos_antinu");
ccqelike = fSettings.Found("name", "CCQELike");
if (!ccqelike) {
- fSettings.SetDataInput( FitPar::GetDataBase() + "/MiniBooNE/anti-ccqe/aski_con.txt" );
+ fSettings.SetDataInput( FitPar::GetDataBase() + "/MiniBooNE/anti-ccqe/aski_like.txt" );
fSettings.SetErrorInput( FitPar::GetDataBase() + "/MiniBooNE/anti-ccqe/aski_err.txt" );
fSettings.DefineAllowedSpecies("numub");
} else {
fSettings.SetDataInput( FitPar::GetDataBase() + "/MiniBooNE/anti-ccqe/aski_like.txt" );
fSettings.SetErrorInput( FitPar::GetDataBase() + "/MiniBooNE/anti-ccqe/aski_err.txt" );
fMeasurementSpeciesType = kNumuWithWrongSignMeasurement;
fSettings.DefineAllowedSpecies("numu,numub");
}
FinaliseSampleSettings();
// Scaling Setup ---------------------------------------------------
// ScaleFactor automatically setup for DiffXSec/cm2/Nucleon
fScaleFactor = ((GetEventHistogram()->Integral("width") * 1E-38 / (fNEvents + 0.)) * (14.08 / 8.) / TotalIntegratedFlux());
// Plot Setup -------------------------------------------------------
Double_t binx[19] = { 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0};
Double_t biny[21] = { -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0};
CreateDataHistogram(19, binx, 21, biny);
SetDataValuesFromTextFile( fSettings.GetDataInput() );
ScaleData(1E-41);
SetDataErrorsFromTextFile( fSettings.GetErrorInput() );
ScaleDataErrors(1E-42);
SetCovarFromDiagonal();
// Final setup ---------------------------------------------------
FinaliseMeasurement();
};
//********************************************************************
void MiniBooNE_CCQE_XSec_2DTcos_antinu::FillEventVariables(FitEvent *event) {
//********************************************************************
if (event->NumFSParticle(PhysConst::pdg_muons) == 0)
return;
TLorentzVector Pnu = event->GetNeutrinoIn()->fP;
// The highest momentum mu+/mu-. The isSignal definition should make sure we only
// accept events we want, so no need to do an additional check here.
TLorentzVector Pmu = event->GetHMFSParticle(PhysConst::pdg_muons)->fP;
// Now find the kinematic values and fill the histogram
Ekmu = Pmu.E() / 1000.0 - PhysConst::mass_muon;
costheta = cos(Pnu.Vect().Angle(Pmu.Vect()));
// Set X Variables
fXVar = Ekmu;
fYVar = costheta;
return;
};
//********************************************************************
bool MiniBooNE_CCQE_XSec_2DTcos_antinu::isSignal(FitEvent *event) {
//********************************************************************
// If CC0pi, include both charges
if (ccqelike) {
if (SignalDef::isCC0pi(event, 14, EnuMin, EnuMax) ||
SignalDef::isCC0pi(event, -14, EnuMin, EnuMax))
return true;
} else {
if (SignalDef::isCCQELike(event, -14, EnuMin, EnuMax))
return true;
}
return false;
};
diff --git a/src/Reweight/FitWeight.cxx b/src/Reweight/FitWeight.cxx
index ffbac6f..50714a6 100644
--- a/src/Reweight/FitWeight.cxx
+++ b/src/Reweight/FitWeight.cxx
@@ -1,242 +1,242 @@
#include "FitWeight.h"
void FitWeight::AddRWEngine(int type) {
switch (type) {
case kNEUT:
fAllRW[type] = new NEUTWeightEngine("neutrw");
break;
case kNUWRO:
fAllRW[type] = new NuWroWeightEngine("nuwrorw");
break;
case kGENIE:
fAllRW[type] = new GENIEWeightEngine("genierw");
break;
case kNORM:
fAllRW[type] = new SampleNormEngine("normrw");
break;
case kLIKEWEIGHT:
fAllRW[type] = new LikelihoodWeightEngine("likerw");
break;
case kT2K:
fAllRW[type] = new T2KWeightEngine("t2krw");
break;
case kCUSTOM:
fAllRW[type] = new NUISANCEWeightEngine("nuisrw");
break;
case kSPLINEPARAMETER:
fAllRW[type] = new SplineWeightEngine("splinerw");
break;
case kNIWG:
fAllRW[type] = new NIWGWeightEngine("niwgrw");
break;
default:
THROW("CANNOT ADD RW Engine for unknown dial type: " << type);
break;
}
}
void FitWeight::IncludeDial(std::string name, std::string type, double val) {
// Should register the dial here.
int typeenum = Reweight::ConvDialType(type);
IncludeDial(name, typeenum, val);
}
void FitWeight::IncludeDial(std::string name, int dialtype, double val) {
// Get the dial enum
int nuisenum = Reweight::ConvDial(name, dialtype);
if (nuisenum == -1){
THROW("NUISENUM == " << nuisenum << " for " << name);
}
// Setup RW Engine Pointer
if (fAllRW.find(dialtype) == fAllRW.end()) {
AddRWEngine(dialtype);
}
WeightEngineBase* rw = fAllRW[dialtype];
// Include the dial
rw->IncludeDial(name, val);
// Set Dial Value
if (val != -9999.9) {
rw->SetDialValue(name, val);
}
// Sort Maps
fAllEnums[name] = nuisenum;
fAllValues[nuisenum] = val;
// Sort Lists
fNameList.push_back(name);
fEnumList.push_back(nuisenum);
fValueList.push_back(val);
}
void FitWeight::Reconfigure(bool silent) {
// Reconfigure all added RW engines
for (std::map<int, WeightEngineBase*>::iterator iter = fAllRW.begin();
iter != fAllRW.end(); iter++) {
(*iter).second->Reconfigure(silent);
}
}
void FitWeight::SetDialValue(std::string name, double val) {
int nuisenum = fAllEnums[name];
SetDialValue(nuisenum, val);
}
// Allow for name aswell using GlobalList to determine sample name.
void FitWeight::SetDialValue(int nuisenum, double val) {
// Conv dial type
int dialtype = int(nuisenum - (nuisenum % 1000)) / 1000;
if (fAllRW.find(dialtype) == fAllRW.end()){
THROW("Cannot find RW Engine for dialtype = " << dialtype << " " << nuisenum << " " << (nuisenum - (nuisenum % 1000)) / 1000);
}
// Get RW Engine for this dial
fAllRW[dialtype]->SetDialValue(nuisenum, val);
fAllValues[nuisenum] = val;
// Update ValueList
for (size_t i = 0; i < fEnumList.size(); i++) {
if (fEnumList[i] == nuisenum) {
fValueList[i] = val;
}
}
}
void FitWeight::SetAllDials(const double* x, int n) {
for (size_t i = 0; i < (UInt_t)n; i++) {
int rwenum = fEnumList[i];
SetDialValue(rwenum, x[i]);
}
Reconfigure();
}
double FitWeight::GetDialValue(std::string name) {
int nuisenum = fAllEnums[name];
return GetDialValue(nuisenum);
}
double FitWeight::GetDialValue(int nuisenum) {
return fAllValues[nuisenum];
}
int FitWeight::GetDialPos(std::string name) {
int rwenum = fAllEnums[name];
return GetDialPos(rwenum);
}
int FitWeight::GetDialPos(int nuisenum) {
for (size_t i = 0; i < fEnumList.size(); i++) {
if (fEnumList[i] == nuisenum) {
return i;
}
}
ERR(FTL) << "No Dial Found! " << std::endl;
throw;
return -1;
}
bool FitWeight::DialIncluded(std::string name) {
return (fAllEnums.find(name) != fAllEnums.end());
}
bool FitWeight::DialIncluded(int rwenum) {
return (fAllValues.find(rwenum) != fAllValues.end());
}
double FitWeight::CalcWeight(BaseFitEvt* evt) {
double rwweight = 1.0;
for (std::map<int, WeightEngineBase*>::iterator iter = fAllRW.begin();
iter != fAllRW.end(); iter++) {
double w = (*iter).second->CalcWeight(evt);
// LOG(FIT) << "Iter " << (*iter).second->fCalcName << " = " << w << std::endl;
rwweight *= w;
}
return rwweight;
}
void FitWeight::UpdateWeightEngine(const double* x) {
size_t count = 0;
for (std::vector<int>::iterator iter = fEnumList.begin();
iter != fEnumList.end(); iter++) {
SetDialValue( (*iter), x[count] );
count++;
}
}
void FitWeight::GetAllDials(double* x, int n) {
for (int i = 0; i < n; i++) {
x[i] = GetDialValue( fEnumList[i] );
}
}
bool FitWeight::NeedsEventReWeight(const double* x) {
bool haschange = false;
size_t count = 0;
// Compare old to new and decide if RW needed.
for (std::vector<int>::iterator iter = fEnumList.begin();
iter != fEnumList.end(); iter++) {
int nuisenum = (*iter);
int type = (nuisenum / 1000) - (nuisenum % 1000);
// Compare old to new
double oldval = GetDialValue(nuisenum);
double newval = x[count];
if (oldval != newval) {
if (fAllRW[type]->NeedsEventReWeight()) {
haschange = true;
}
}
count++;
}
return haschange;
}
double FitWeight::GetSampleNorm(std::string name) {
if (name.empty()) return 1.0;
// Find norm dial
- if (fAllEnums.find(name) != fAllEnums.end()) {
- if (fAllValues.find(fAllEnums[name]) != fAllValues.end()){
- return fAllValues[ fAllEnums[name] ];
+ if (fAllEnums.find(name + "_norm") != fAllEnums.end()) {
+ if (fAllValues.find(fAllEnums[name+"_norm"]) != fAllValues.end()){
+ return fAllValues[ fAllEnums[name+"_norm"] ];
} else {
return 1.0;
}
} else {
return 1.0;
}
}
void FitWeight::Print() {
LOG(REC) << "Fit Weight State: " << std::endl;
for (size_t i = 0; i < fNameList.size(); i++) {
LOG(REC) << " -> Par " << i << ". " << fNameList[i] << " " << fValueList[i] << std::endl;
}
}
diff --git a/src/Routines/ComparisonRoutines.cxx b/src/Routines/ComparisonRoutines.cxx
index b0e3b21..9ccfff6 100755
--- a/src/Routines/ComparisonRoutines.cxx
+++ b/src/Routines/ComparisonRoutines.cxx
@@ -1,531 +1,531 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "ComparisonRoutines.h"
/*
Constructor/Destructor
*/
//************************
void ComparisonRoutines::Init() {
//************************
fOutputFile = "";
fOutputRootFile = NULL;
fStrategy = "Compare";
fRoutines.clear();
fCardFile = "";
fFakeDataInput = "";
fSampleFCN = NULL;
fAllowedRoutines = ("Compare");
};
//*************************************
ComparisonRoutines::~ComparisonRoutines() {
//*************************************
};
/*
Input Functions
*/
//*************************************
ComparisonRoutines::ComparisonRoutines(int argc, char* argv[]) {
//*************************************
// Initialise Defaults
Init();
nuisconfig configuration = Config::Get();
// Default containers
std::string cardfile = "";
std::string maxevents = "-1";
int errorcount = 0;
int verbocount = 0;
std::vector<std::string> xmlcmds;
std::vector<std::string> configargs;
// Make easier to handle arguments.
std::vector<std::string> args = GeneralUtils::LoadCharToVectStr(argc, argv);
ParserUtils::ParseArgument(args, "-c", fCardFile, true);
ParserUtils::ParseArgument(args, "-o", fOutputFile, false, false);
ParserUtils::ParseArgument(args, "-n", maxevents, false, false);
ParserUtils::ParseArgument(args, "-f", fStrategy, false, false);
ParserUtils::ParseArgument(args, "-d", fFakeDataInput, false, false);
ParserUtils::ParseArgument(args, "-i", xmlcmds);
ParserUtils::ParseArgument(args, "-q", configargs);
ParserUtils::ParseCounter(args, "e", errorcount);
ParserUtils::ParseCounter(args, "v", verbocount);
ParserUtils::CheckBadArguments(args);
// Add extra defaults if none given
if (fCardFile.empty() and xmlcmds.empty()) {
ERR(FTL) << "No input supplied!" << std::endl;
throw;
}
if (fOutputFile.empty() and !fCardFile.empty()) {
fOutputFile = fCardFile + ".root";
ERR(WRN) << "No output supplied so saving it to: " << fOutputFile << std::endl;
} else if (fOutputFile.empty()) {
ERR(FTL) << "No output file or cardfile supplied!" << std::endl;
throw;
}
// Configuration Setup =============================
// Check no comp key is available
nuiskey fCompKey;
if (Config::Get().GetNodes("nuiscomp").empty()) {
fCompKey = Config::Get().CreateNode("nuiscomp");
} else {
fCompKey = Config::Get().GetNodes("nuiscomp")[0];
}
if (!fCardFile.empty()) fCompKey.AddS("cardfile", fCardFile);
if (!fOutputFile.empty()) fCompKey.AddS("outputfile", fOutputFile);
if (!fStrategy.empty()) fCompKey.AddS("strategy", fStrategy);
// Load XML Cardfile
configuration.LoadConfig( fCompKey.GetS("cardfile"), "");
// Add CMD XML Structs
for (size_t i = 0; i < xmlcmds.size(); i++) {
configuration.AddXMLLine(xmlcmds[i]);
}
// Add Config Args
for (size_t i = 0; i < configargs.size(); i++) {
configuration.OverrideConfig(configargs[i]);
}
if (maxevents.compare("-1")){
configuration.OverrideConfig("MAXEVENTS=" + maxevents);
}
// Finish configuration XML
configuration.FinaliseConfig(fCompKey.GetS("outputfile") + ".xml");
// Add Error Verbo Lines
verbocount += Config::Get().GetParI("VERBOSITY");
errorcount += Config::Get().GetParI("ERROR");
bool trace = Config::Get().GetParB("TRACE");
std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl;
std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl;
SETVERBOSITY(verbocount);
SETTRACE(trace);
// Comparison Setup ========================================
// Proper Setup
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
SetupComparisonsFromXML();
SetupRWEngine();
SetupFCN();
return;
};
//*************************************
void ComparisonRoutines::SetupComparisonsFromXML() {
//*************************************
LOG(FIT) << "Setting up nuiscomp" << std::endl;
// Setup Parameters ------------------------------------------
std::vector<nuiskey> parkeys = Config::QueryKeys("parameter");
if (!parkeys.empty()) {
LOG(FIT) << "Number of parameters : " << parkeys.size() << std::endl;
}
for (size_t i = 0; i < parkeys.size(); i++) {
nuiskey key = parkeys.at(i);
// Check for type,name,nom
if (!key.Has("type")) {
ERR(FTL) << "No type given for parameter " << i << std::endl;
ERR(FTL) << "type='PARAMETER_TYPE'" << std::endl;
throw;
} else if (!key.Has("name")) {
ERR(FTL) << "No name given for parameter " << i << std::endl;
ERR(FTL) << "name='SAMPLE_NAME'" << std::endl;
throw;
} else if (!key.Has("nominal")) {
ERR(FTL) << "No nominal given for parameter " << i << std::endl;
ERR(FTL) << "nominal='NOMINAL_VALUE'" << std::endl;
throw;
}
// Get Inputs
std::string partype = key.GetS("type");
std::string parname = key.GetS("name");
double parnom = key.GetD("nominal");
double parlow = parnom - 1;
double parhigh = parnom + 1;
double parstep = 1;
std::string parstate = key.GetS("state");
// Check for incomplete limtis
int limdef = ((int)key.Has("low") +
(int)key.Has("high") +
(int)key.Has("step"));
if (limdef > 0 and limdef < 3){
ERR(FTL) << "Incomplete limit set given for parameter : " << parname << std::endl;
ERR(FTL) << "Requires: low='LOWER_LIMIT' high='UPPER_LIMIT' step='STEP_SIZE' " << std::endl;
throw;
}
// Extra limits
if (key.Has("low")) {
parlow = key.GetD("low");
parhigh = key.GetD("high");
parstep = key.GetD("step");
LOG(FIT) << "Read " << partype << " : "
<< parname << " = "
<< parnom << " : "
<< parlow << " < p < " << parhigh
<< " : " << parstate << std::endl;
} else {
LOG(FIT) << "Read " << partype << " : "
<< parname << " = "
<< parnom << " : "
<< parstate << std::endl;
}
// Convert if required
if (parstate.find("ABS") != std::string::npos) {
parnom = FitBase::RWAbsToSigma( partype, parname, parnom );
parlow = FitBase::RWAbsToSigma( partype, parname, parlow );
parhigh = FitBase::RWAbsToSigma( partype, parname, parhigh );
parstep = FitBase::RWAbsToSigma( partype, parname, parstep );
} else if (parstate.find("FRAC") != std::string::npos) {
parnom = FitBase::RWFracToSigma( partype, parname, parnom );
parlow = FitBase::RWFracToSigma( partype, parname, parlow );
parhigh = FitBase::RWFracToSigma( partype, parname, parhigh );
parstep = FitBase::RWFracToSigma( partype, parname, parstep );
}
// Push into vectors
fParams.push_back(parname);
fTypeVals[parname] = FitBase::ConvDialType(partype);;
fCurVals[parname] = parnom;
fStateVals[parname] = parstate;
}
// Setup Samples ----------------------------------------------
std::vector<nuiskey> samplekeys = Config::QueryKeys("sample");
if (!samplekeys.empty()) {
LOG(FIT) << "Number of samples : " << samplekeys.size() << std::endl;
}
for (size_t i = 0; i < samplekeys.size(); i++) {
nuiskey key = samplekeys.at(i);
// Get Sample Options
std::string samplename = key.GetS("name");
std::string samplefile = key.GetS("input");
std::string sampletype =
key.Has("type") ? key.GetS("type") : "DEFAULT";
double samplenorm =
key.Has("norm") ? key.GetD("norm") : 1.0;
// Print out
LOG(FIT) << "Read Sample " << i << ". : "
<< samplename << " (" << sampletype << ") [Norm=" << samplenorm<<"]"<< std::endl
<< " -> input='" << samplefile << "'" << std::endl;
// If FREE add to parameters otherwise continue
if (sampletype.find("FREE") == std::string::npos) {
continue;
}
// Form norm dial from samplename + sampletype + "_norm";
- std::string normname = samplename + sampletype + "_norm";
+ std::string normname = samplename + "_norm";
// Check normname not already present
if (fTypeVals.find("normname") != fTypeVals.end()) {
continue;
}
// Add new norm dial to list if its passed above checks
fParams.push_back(normname);
fTypeVals[normname] = kNORM;
fStateVals[normname] = sampletype;
fCurVals[normname] = samplenorm;
}
// Setup Fake Parameters -----------------------------
std::vector<nuiskey> fakekeys = Config::QueryKeys("fakeparameter");
if (!fakekeys.empty()) {
LOG(FIT) << "Number of fake parameters : " << fakekeys.size() << std::endl;
}
for (size_t i = 0; i < fakekeys.size(); i++) {
nuiskey key = fakekeys.at(i);
// Check for type,name,nom
if (!key.Has("name")) {
ERR(FTL) << "No name given for fakeparameter " << i << std::endl;
throw;
} else if (!key.Has("nominal")) {
ERR(FTL) << "No nominal given for fakeparameter " << i << std::endl;
throw;
}
// Get Inputs
std::string parname = key.GetS("name");
double parnom = key.GetD("nominal");
// Push into vectors
fFakeVals[parname] = parnom;
}
}
//*************************************
void ComparisonRoutines::SetupRWEngine() {
//*************************************
LOG(FIT) << "Setting up FitWeight Engine" << std::endl;
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string name = fParams[i];
FitBase::GetRW()->IncludeDial(name, fTypeVals.at(name));
}
return;
}
//*************************************
void ComparisonRoutines::SetupFCN() {
//*************************************
LOG(FIT) << "Building the SampleFCN" << std::endl;
if (fSampleFCN) delete fSampleFCN;
FitPar::Config().out = fOutputRootFile;
fOutputRootFile->cd();
fSampleFCN = new JointFCN(fOutputRootFile);
SetFakeData();
return;
}
//*************************************
void ComparisonRoutines::SetFakeData() {
//*************************************
if (fFakeDataInput.empty()) return;
if (fFakeDataInput.compare("MC") == 0) {
LOG(FIT) << "Setting fake data from MC starting prediction." << std::endl;
UpdateRWEngine(fFakeVals);
FitBase::GetRW()->Reconfigure();
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->SetFakeData("MC");
UpdateRWEngine(fCurVals);
LOG(FIT) << "Set all data to fake MC predictions." << std::endl;
} else {
LOG(FIT) << "Setting fake data from: " << fFakeDataInput << std::endl;
fSampleFCN->SetFakeData(fFakeDataInput);
}
return;
}
/*
Fitting Functions
*/
//*************************************
void ComparisonRoutines::UpdateRWEngine(
std::map<std::string, double>& updateVals) {
//*************************************
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string name = fParams[i];
if (updateVals.find(name) == updateVals.end()) continue;
FitBase::GetRW()->SetDialValue(name, updateVals.at(name));
}
FitBase::GetRW()->Reconfigure();
return;
}
//*************************************
void ComparisonRoutines::Run() {
//*************************************
LOG(FIT) << "Running ComparisonRoutines : " << fStrategy << std::endl;
if (FitPar::Config().GetParB("save_nominal")) {
SaveNominal();
}
// Parse given routines
fRoutines = GeneralUtils::ParseToStr(fStrategy, ",");
if (fRoutines.empty()) {
ERR(FTL) << "Trying to run ComparisonRoutines with no routines given!" << std::endl;
throw;
}
for (UInt_t i = 0; i < fRoutines.size(); i++) {
std::string routine = fRoutines.at(i);
LOG(FIT) << "Routine: " << routine << std::endl;
if (!routine.compare("Compare")) {
UpdateRWEngine(fCurVals);
GenerateComparison();
PrintState();
SaveCurrentState();
}
}
return;
}
//*************************************
void ComparisonRoutines::GenerateComparison() {
//*************************************
LOG(FIT) << "Generating Comparison." << std::endl;
// Main Event Loop from event Manager
fSampleFCN->ReconfigureAllEvents();
return;
}
//*************************************
void ComparisonRoutines::PrintState() {
//*************************************
LOG(FIT) << "------------" << std::endl;
// Count max size
int maxcount = 0;
for (UInt_t i = 0; i < fParams.size(); i++) {
maxcount = max(int(fParams[i].size()), maxcount);
}
// Header
LOG(FIT) << " # " << left << setw(maxcount) << "Parameter "
<< " = " << setw(10) << "Value"
<< " +- " << setw(10) << "Error"
<< " " << setw(8) << "(Units)"
<< " " << setw(10) << "Conv. Val"
<< " +- " << setw(10) << "Conv. Err"
<< " " << setw(8) << "(Units)" << std::endl;
// Parameters
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string syst = fParams.at(i);
std::string typestr = FitBase::ConvDialType(fTypeVals[syst]);
std::string curunits = "(sig.)";
double curval = fCurVals[syst];
double curerr = 0.0;
if (fStateVals[syst].find("ABS") != std::string::npos) {
curval = FitBase::RWSigmaToAbs(typestr, syst, curval);
curerr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
FitBase::RWSigmaToAbs(typestr, syst, 0.0));
curunits = "(Abs.)";
} else if (fStateVals[syst].find("FRAC") != std::string::npos) {
curval = FitBase::RWSigmaToFrac(typestr, syst, curval);
curerr = (FitBase::RWSigmaToFrac(typestr, syst, curerr) -
FitBase::RWSigmaToFrac(typestr, syst, 0.0));
curunits = "(Frac)";
}
std::string convunits = "(" + FitBase::GetRWUnits(typestr, syst) + ")";
double convval = FitBase::RWSigmaToAbs(typestr, syst, curval);
double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
FitBase::RWSigmaToAbs(typestr, syst, 0.0));
std::ostringstream curparstring;
curparstring << " " << setw(3) << left << i << ". " << setw(maxcount)
<< syst << " = " << setw(10) << curval << " +- " << setw(10)
<< curerr << " " << setw(8) << curunits << " " << setw(10)
<< convval << " +- " << setw(10) << converr << " " << setw(8)
<< convunits;
LOG(FIT) << curparstring.str() << std::endl;
}
LOG(FIT) << "------------" << std::endl;
double like = fSampleFCN->GetLikelihood();
LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like << std::endl;
LOG(FIT) << "------------" << std::endl;
}
/*
Write Functions
*/
//*************************************
void ComparisonRoutines::SaveCurrentState(std::string subdir) {
//*************************************
LOG(FIT) << "Saving current full FCN predictions" << std::endl;
// Setup DIRS
TDirectory* curdir = gDirectory;
if (!subdir.empty()) {
TDirectory* newdir = (TDirectory*)gDirectory->mkdir(subdir.c_str());
newdir->cd();
}
fSampleFCN->Write();
// Change back to current DIR
curdir->cd();
return;
}
//*************************************
void ComparisonRoutines::SaveNominal() {
//*************************************
fOutputRootFile->cd();
LOG(FIT) << "Saving Nominal Predictions (be cautious with this)" << std::endl;
FitBase::GetRW()->Reconfigure();
GenerateComparison();
SaveCurrentState("nominal");
};
diff --git a/src/Routines/MinimizerRoutines.cxx b/src/Routines/MinimizerRoutines.cxx
index 61cd91b..de1cd1a 100755
--- a/src/Routines/MinimizerRoutines.cxx
+++ b/src/Routines/MinimizerRoutines.cxx
@@ -1,1450 +1,1470 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "MinimizerRoutines.h"
/*
Constructor/Destructor
*/
//************************
void MinimizerRoutines::Init() {
//************************
fInputFile = "";
fInputRootFile = NULL;
fOutputFile = "";
fOutputRootFile = NULL;
fCovar = NULL;
fCovFree = NULL;
fCorrel = NULL;
fCorFree = NULL;
fDecomp = NULL;
fDecFree = NULL;
fStrategy = "Migrad,FixAtLimBreak,Migrad";
fRoutines.clear();
fCardFile = "";
fFakeDataInput = "";
fSampleFCN = NULL;
fMinimizer = NULL;
fMinimizerFCN = NULL;
fCallFunctor = NULL;
fAllowedRoutines = ("Migrad,Simplex,Combined,"
"Brute,Fumili,ConjugateFR,"
"ConjugatePR,BFGS,BFGS2,"
"SteepDesc,GSLSimAn,FixAtLim,FixAtLimBreak"
"Chi2Scan1D,Chi2Scan2D,Contours,ErrorBands");
};
//*************************************
MinimizerRoutines::~MinimizerRoutines() {
//*************************************
};
/*
Input Functions
*/
//*************************************
MinimizerRoutines::MinimizerRoutines(int argc, char* argv[]) {
//*************************************
// Initialise Defaults
Init();
nuisconfig configuration = Config::Get();
// Default containers
std::string cardfile = "";
std::string maxevents = "-1";
int errorcount = 0;
int verbocount = 0;
std::vector<std::string> xmlcmds;
std::vector<std::string> configargs;
// Make easier to handle arguments.
std::vector<std::string> args = GeneralUtils::LoadCharToVectStr(argc, argv);
ParserUtils::ParseArgument(args, "-c", fCardFile, true);
ParserUtils::ParseArgument(args, "-o", fOutputFile, false, false);
ParserUtils::ParseArgument(args, "-n", maxevents, false, false);
ParserUtils::ParseArgument(args, "-f", fStrategy, false, false);
ParserUtils::ParseArgument(args, "-d", fFakeDataInput, false, false);
ParserUtils::ParseArgument(args, "-i", xmlcmds);
ParserUtils::ParseArgument(args, "-q", configargs);
ParserUtils::ParseCounter(args, "e", errorcount);
ParserUtils::ParseCounter(args, "v", verbocount);
ParserUtils::CheckBadArguments(args);
// Add extra defaults if none given
if (fCardFile.empty() and xmlcmds.empty()) {
ERR(FTL) << "No input supplied!" << std::endl;
throw;
}
if (fOutputFile.empty() and !fCardFile.empty()) {
fOutputFile = fCardFile + ".root";
ERR(WRN) << "No output supplied so saving it to: " << fOutputFile << std::endl;
} else if (fOutputFile.empty()) {
ERR(FTL) << "No output file or cardfile supplied!" << std::endl;
throw;
}
// Configuration Setup =============================
// Check no comp key is available
nuiskey fCompKey;
if (Config::Get().GetNodes("nuiscomp").empty()) {
fCompKey = Config::Get().CreateNode("nuiscomp");
} else {
fCompKey = Config::Get().GetNodes("nuiscomp")[0];
}
if (!fCardFile.empty()) fCompKey.AddS("cardfile", fCardFile);
if (!fOutputFile.empty()) fCompKey.AddS("outputfile", fOutputFile);
if (!fStrategy.empty()) fCompKey.AddS("strategy", fStrategy);
// Load XML Cardfile
configuration.LoadConfig( fCompKey.GetS("cardfile"), "");
// Add CMD XML Structs
for (size_t i = 0; i < xmlcmds.size(); i++) {
configuration.AddXMLLine(xmlcmds[i]);
}
// Add Config Args
for (size_t i = 0; i < configargs.size(); i++) {
configuration.OverrideConfig(configargs[i]);
}
if (maxevents.compare("-1")){
configuration.OverrideConfig("MAXEVENTS=" + maxevents);
}
// Finish configuration XML
configuration.FinaliseConfig(fCompKey.GetS("outputfile") + ".xml");
// Add Error Verbo Lines
verbocount += Config::Get().GetParI("VERBOSITY");
errorcount += Config::Get().GetParI("ERROR");
std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl;
std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl;
// FitPar::log_verb = verbocount;
SETVERBOSITY(verbocount);
// ERR_VERB(errorcount);
// Minimizer Setup ========================================
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
SetupMinimizerFromXML();
SetupCovariance();
SetupRWEngine();
SetupFCN();
return;
};
//*************************************
void MinimizerRoutines::SetupMinimizerFromXML() {
//*************************************
LOG(FIT) << "Setting up nuismin" << std::endl;
// Setup Parameters ------------------------------------------
std::vector<nuiskey> parkeys = Config::QueryKeys("parameter");
if (!parkeys.empty()) {
LOG(FIT) << "Number of parameters : " << parkeys.size() << std::endl;
}
for (size_t i = 0; i < parkeys.size(); i++) {
nuiskey key = parkeys.at(i);
// Check for type,name,nom
if (!key.Has("type")) {
ERR(FTL) << "No type given for parameter " << i << std::endl;
throw;
} else if (!key.Has("name")) {
ERR(FTL) << "No name given for parameter " << i << std::endl;
throw;
} else if (!key.Has("nominal")) {
ERR(FTL) << "No nominal given for parameter " << i << std::endl;
throw;
}
// Get Inputs
std::string partype = key.GetS("type");
std::string parname = key.GetS("name");
double parnom = key.GetD("nominal");
double parlow = parnom - 1;
double parhigh = parnom + 1;
double parstep = 1;
std::string parstate = key.GetS("state");
// Extra limits
if (key.Has("low")) {
parlow = key.GetD("low");
parhigh = key.GetD("high");
parstep = key.GetD("step");
LOG(FIT) << "Read " << partype << " : "
<< parname << " = "
<< parnom << " : "
<< parlow << " < p < " << parhigh
<< " : " << parstate << std::endl;
} else {
LOG(FIT) << "Read " << partype << " : "
<< parname << " = "
<< parnom << " : "
<< parstate << std::endl;
}
// Run Parameter Conversion if needed
if (parstate.find("ABS") != std::string::npos) {
parnom = FitBase::RWAbsToSigma( partype, parname, parnom );
parlow = FitBase::RWAbsToSigma( partype, parname, parlow );
parhigh = FitBase::RWAbsToSigma( partype, parname, parhigh );
parstep = FitBase::RWAbsToSigma( partype, parname, parstep );
} else if (parstate.find("FRAC") != std::string::npos) {
parnom = FitBase::RWFracToSigma( partype, parname, parnom );
parlow = FitBase::RWFracToSigma( partype, parname, parlow );
parhigh = FitBase::RWFracToSigma( partype, parname, parhigh );
parstep = FitBase::RWFracToSigma( partype, parname, parstep );
}
// Push into vectors
fParams.push_back(parname);
fTypeVals[parname] = FitBase::ConvDialType(partype);;
fStartVals[parname] = parnom;
fCurVals[parname] = parnom;
fErrorVals[parname] = 0.0;
fStateVals[parname] = parstate;
bool fixstate = parstate.find("FIX") != std::string::npos;
fFixVals[parname] = fixstate;
fStartFixVals[parname] = fFixVals[parname];
fMinVals[parname] = parlow;
fMaxVals[parname] = parhigh;
fStepVals[parname] = parstep;
}
// Setup Samples ----------------------------------------------
std::vector<nuiskey> samplekeys = Config::QueryKeys("sample");
if (!samplekeys.empty()) {
LOG(FIT) << "Number of samples : " << samplekeys.size() << std::endl;
}
for (size_t i = 0; i < samplekeys.size(); i++) {
nuiskey key = samplekeys.at(i);
// Get Sample Options
std::string samplename = key.GetS("name");
std::string samplefile = key.GetS("input");
std::string sampletype =
key.Has("type") ? key.GetS("type") : "DEFAULT";
double samplenorm =
key.Has("norm") ? key.GetD("norm") : 1.0;
// Print out
LOG(FIT) << "Read sample info " << i << " : "
<< samplename << std::endl
<< "\t\t input -> " << samplefile << std::endl
<< "\t\t state -> " << sampletype << std::endl
<< "\t\t norm -> " << samplenorm << std::endl;
// If FREE add to parameters otherwise continue
if (sampletype.find("FREE") == std::string::npos) {
continue;
}
// Form norm dial from samplename + sampletype + "_norm";
- std::string normname = samplename + sampletype + "_norm";
+ std::string normname = samplename + "_norm";
// Check normname not already present
if (fTypeVals.find(normname) != fTypeVals.end()) {
continue;
}
// Add new norm dial to list if its passed above checks
fParams.push_back(normname);
fTypeVals[normname] = kNORM;
fStateVals[normname] = sampletype;
+ fStartVals[normname] = samplenorm;
fCurVals[normname] = samplenorm;
fErrorVals[normname] = 0.0;
fMinVals[normname] = 0.1;
fMaxVals[normname] = 10.0;
fStepVals[normname] = 0.5;
bool state = sampletype.find("FREE") == std::string::npos;
fFixVals[normname] = state;
fStartFixVals[normname] = state;
+
+
+
}
// Setup Fake Parameters -----------------------------
std::vector<nuiskey> fakekeys = Config::QueryKeys("fakeparameter");
if (!fakekeys.empty()) {
LOG(FIT) << "Number of fake parameters : " << fakekeys.size() << std::endl;
}
for (size_t i = 0; i < fakekeys.size(); i++) {
nuiskey key = fakekeys.at(i);
// Check for type,name,nom
if (!key.Has("name")) {
ERR(FTL) << "No name given for fakeparameter " << i << std::endl;
throw;
} else if (!key.Has("nom")) {
ERR(FTL) << "No nominal given for fakeparameter " << i << std::endl;
throw;
}
// Get Inputs
std::string parname = key.GetS("name");
double parnom = key.GetD("nom");
// Push into vectors
fFakeVals[parname] = parnom;
}
}
/*
Setup Functions
*/
//*************************************
void MinimizerRoutines::SetupRWEngine() {
//*************************************
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string name = fParams[i];
FitBase::GetRW() -> IncludeDial(name, fTypeVals.at(name) );
}
UpdateRWEngine(fStartVals);
return;
}
//*************************************
void MinimizerRoutines::SetupFCN() {
//*************************************
LOG(FIT) << "Making the jointFCN" << std::endl;
if (fSampleFCN) delete fSampleFCN;
// fSampleFCN = new JointFCN(fCardFile, fOutputRootFile);
fSampleFCN = new JointFCN(fOutputRootFile);
SetFakeData();
fMinimizerFCN = new MinimizerFCN( fSampleFCN );
fCallFunctor = new ROOT::Math::Functor( *fMinimizerFCN, fParams.size() );
fSampleFCN->CreateIterationTree( "fit_iterations", FitBase::GetRW() );
return;
}
//******************************************
void MinimizerRoutines::SetupFitter(std::string routine) {
//******************************************
// Make the fitter
std::string fitclass = "";
std::string fittype = "";
// Get correct types
if (!routine.compare("Migrad")) {
fitclass = "Minuit2"; fittype = "Migrad";
} else if (!routine.compare("Simplex")) {
fitclass = "Minuit2"; fittype = "Simplex";
} else if (!routine.compare("Combined")) {
fitclass = "Minuit2"; fittype = "Combined";
} else if (!routine.compare("Brute")) {
fitclass = "Minuit2"; fittype = "Scan";
} else if (!routine.compare("Fumili")) {
fitclass = "Minuit2"; fittype = "Fumili";
} else if (!routine.compare("ConjugateFR")) {
fitclass = "GSLMultiMin"; fittype = "ConjugateFR";
} else if (!routine.compare("ConjugatePR")) {
fitclass = "GSLMultiMin"; fittype = "ConjugatePR";
} else if (!routine.compare("BFGS")) {
fitclass = "GSLMultiMin"; fittype = "BFGS";
} else if (!routine.compare("BFGS2")) {
fitclass = "GSLMultiMin"; fittype = "BFGS2";
} else if (!routine.compare("SteepDesc")) {
fitclass = "GSLMultiMin"; fittype = "SteepestDescent";
// } else if (!routine.compare("GSLMulti")) { fitclass = "GSLMultiFit"; fittype = ""; // Doesn't work out of the box
} else if (!routine.compare("GSLSimAn")) { fitclass = "GSLSimAn"; fittype = ""; }
// make minimizer
if (fMinimizer) delete fMinimizer;
fMinimizer = ROOT::Math::Factory::CreateMinimizer(fitclass, fittype);
fMinimizer->SetMaxFunctionCalls(FitPar::Config().GetParI("minimizer.maxcalls"));
if (!routine.compare("Brute")) {
fMinimizer->SetMaxFunctionCalls(fParams.size() * fParams.size() * 4);
fMinimizer->SetMaxIterations(fParams.size() * fParams.size() * 4);
}
fMinimizer->SetMaxIterations(FitPar::Config().GetParI("minimizer.maxiterations"));
fMinimizer->SetTolerance(FitPar::Config().GetParD("minimizer.tolerance"));
fMinimizer->SetStrategy(FitPar::Config().GetParI("minimizer.strategy"));
fMinimizer->SetFunction(*fCallFunctor);
int ipar = 0;
//Add Fit Parameters
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string syst = fParams.at(i);
bool fixed = true;
double vstart, vstep, vlow, vhigh;
vstart = vstep = vlow = vhigh = 0.0;
if (fCurVals.find(syst) != fCurVals.end() ) vstart = fCurVals.at(syst);
if (fMinVals.find(syst) != fMinVals.end() ) vlow = fMinVals.at(syst);
if (fMaxVals.find(syst) != fMaxVals.end() ) vhigh = fMaxVals.at(syst);
if (fStepVals.find(syst) != fStepVals.end()) vstep = fStepVals.at(syst);
if (fFixVals.find(syst) != fFixVals.end() ) fixed = fFixVals.at(syst);
// fix for errors
if (vhigh == vlow) vhigh += 1.0;
fMinimizer->SetVariable(ipar, syst, vstart, vstep);
fMinimizer->SetVariableLimits(ipar, vlow, vhigh);
if (fixed) {
fMinimizer->FixVariable(ipar);
LOG(FIT) << "Fixed Param: " << syst << std::endl;
} else {
LOG(FIT) << "Free Param: " << syst
<< " Start:" << vstart
<< " Range:" << vlow << " to " << vhigh
<< " Step:" << vstep << std::endl;
}
ipar++;
}
LOG(FIT) << "Setup Minimizer: " << fMinimizer->NDim() << "(NDim) " << fMinimizer->NFree() << "(NFree)" << std::endl;
return;
}
//*************************************
// Set fake data from user input
void MinimizerRoutines::SetFakeData() {
//*************************************
// If the fake data input field (-d) isn't provided, return to caller
if (fFakeDataInput.empty()) return;
// If user specifies -d MC we set the data to the MC
// User can also specify fake data parameters to reweight by doing "fake_parameter" in input card file
// "fake_parameter" gets read in ReadCard function (reads to fFakeVals)
if (fFakeDataInput.compare("MC") == 0) {
LOG(FIT) << "Setting fake data from MC starting prediction." << std::endl;
// fFakeVals get read in in ReadCard
UpdateRWEngine(fFakeVals);
// Reconfigure the reweight engine
FitBase::GetRW()->Reconfigure();
// Reconfigure all the samples to the new reweight
fSampleFCN->ReconfigureAllEvents();
// Feed on and set the fake-data in each measurement class
fSampleFCN->SetFakeData("MC");
// Changed the reweight engine values back to the current values
// So we start the fit at a different value than what we set the fake-data to
UpdateRWEngine(fCurVals);
LOG(FIT) << "Set all data to fake MC predictions." << std::endl;
} else {
fSampleFCN->SetFakeData(fFakeDataInput);
}
return;
}
/*
Fitting Functions
*/
//*************************************
void MinimizerRoutines::UpdateRWEngine(std::map<std::string, double>& updateVals) {
//*************************************
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string name = fParams[i];
if (updateVals.find(name) == updateVals.end()) continue;
FitBase::GetRW()->SetDialValue(name, updateVals.at(name));
}
FitBase::GetRW()->Reconfigure();
return;
}
//*************************************
void MinimizerRoutines::Run() {
//*************************************
LOG(FIT) << "Running MinimizerRoutines : " << fStrategy << std::endl;
if (FitPar::Config().GetParB("save_nominal")) {
SaveNominal();
}
// Parse given routines
fRoutines = GeneralUtils::ParseToStr(fStrategy,",");
if (fRoutines.empty()){
ERR(FTL) << "Trying to run MinimizerRoutines with no routines given!" << std::endl;
throw;
}
for (UInt_t i = 0; i < fRoutines.size(); i++) {
std::string routine = fRoutines.at(i);
int fitstate = kFitUnfinished;
LOG(FIT) << "Running Routine: " << routine << std::endl;
// Try Routines
if (routine.find("LowStat") != std::string::npos) LowStatRoutine(routine);
else if (routine == "FixAtLim") FixAtLimit();
else if (routine == "FixAtLimBreak") fitstate = FixAtLimit();
else if (routine.find("ErrorBands") != std::string::npos) GenerateErrorBands();
else if (!routine.compare("Chi2Scan1D")) Create1DScans();
else if (!routine.compare("Chi2Scan2D")) Chi2Scan2D();
else fitstate = RunFitRoutine(routine);
// If ending early break here
if (fitstate == kFitFinished || fitstate == kNoChange) {
LOG(FIT) << "Ending fit routines loop." << std::endl;
break;
}
}
return;
}
//*************************************
int MinimizerRoutines::RunFitRoutine(std::string routine) {
//*************************************
int endfits = kFitUnfinished;
// set fitter at the current start values
fOutputRootFile->cd();
SetupFitter(routine);
// choose what to do with the minimizer depending on routine.
if (!routine.compare("Migrad") or
!routine.compare("Simplex") or
!routine.compare("Combined") or
!routine.compare("Brute") or
!routine.compare("Fumili") or
!routine.compare("ConjugateFR") or
!routine.compare("ConjugatePR") or
!routine.compare("BFGS") or
!routine.compare("BFGS2") or
!routine.compare("SteepDesc") or
// !routine.compare("GSLMulti") or
!routine.compare("GSLSimAn")) {
if (fMinimizer->NFree() > 0) {
LOG(FIT) << fMinimizer->Minimize() << std::endl;
GetMinimizerState();
}
}
// other otptions
else if (!routine.compare("Contour")) {
CreateContours();
}
return endfits;
}
//*************************************
void MinimizerRoutines::PrintState() {
//*************************************
LOG(FIT) << "------------" << std::endl;
// Count max size
int maxcount = 0;
for (UInt_t i = 0; i < fParams.size(); i++) {
maxcount = max(int(fParams[i].size()), maxcount);
}
// Header
LOG(FIT) << " # " << left << setw(maxcount) << "Parameter "
<< " = "
<< setw(10) << "Value" << " +- "
<< setw(10) << "Error" << " "
<< setw(8) << "(Units)" << " "
<< setw(10) << "Conv. Val" << " +- "
<< setw(10) << "Conv. Err" << " "
<< setw(8) << "(Units)" << std::endl;
// Parameters
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string syst = fParams.at(i);
std::string typestr = FitBase::ConvDialType(fTypeVals[syst]);
std::string curunits = "(sig.)";
double curval = fCurVals[syst];
double curerr = fErrorVals[syst];
if (fStateVals[syst].find("ABS") != std::string::npos) {
curval = FitBase::RWSigmaToAbs(typestr, syst, curval);
curerr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
FitBase::RWSigmaToAbs(typestr, syst, 0.0));
curunits = "(Abs.)";
} else if (fStateVals[syst].find("FRAC") != std::string::npos) {
curval = FitBase::RWSigmaToFrac(typestr, syst, curval);
curerr = (FitBase::RWSigmaToFrac(typestr, syst, curerr) -
FitBase::RWSigmaToFrac(typestr, syst, 0.0));
curunits = "(Frac)";
}
std::string convunits = "(" + FitBase::GetRWUnits(typestr, syst) + ")";
double convval = FitBase::RWSigmaToAbs(typestr, syst, curval);
double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
FitBase::RWSigmaToAbs(typestr, syst, 0.0));
std::ostringstream curparstring;
curparstring << " " << setw(3) << left
<< i << ". "
<< setw(maxcount) << syst << " = "
<< setw(10) << curval << " +- "
<< setw(10) << curerr << " "
<< setw(8) << curunits << " "
<< setw(10) << convval << " +- "
<< setw(10) << converr << " "
<< setw(8) << convunits;
LOG(FIT) << curparstring.str() << std::endl;
}
LOG(FIT) << "------------" << std::endl;
double like = fSampleFCN->GetLikelihood();
LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like << std::endl;
LOG(FIT) << "------------" << std::endl;
}
//*************************************
void MinimizerRoutines::GetMinimizerState() {
//*************************************
LOG(FIT) << "Minimizer State: " << std::endl;
// Get X and Err
const double *values = fMinimizer->X();
const double *errors = fMinimizer->Errors();
// int ipar = 0;
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string syst = fParams.at(i);
fCurVals[syst] = values[i];
fErrorVals[syst] = errors[i];
}
PrintState();
// Covar
SetupCovariance();
if (fMinimizer->CovMatrixStatus() > 0) {
// Fill Full Covar
+ std::cout << "Filling covariance" << std::endl;
for (int i = 0; i < fCovar->GetNbinsX(); i++) {
for (int j = 0; j < fCovar->GetNbinsY(); j++) {
fCovar->SetBinContent(i + 1, j + 1, fMinimizer->CovMatrix(i, j));
}
}
int freex = 0;
int freey = 0;
for (int i = 0; i < fCovar->GetNbinsX(); i++) {
if (fMinimizer->IsFixedVariable(i)) continue;
freey = 0;
for (int j = 0; j < fCovar->GetNbinsY(); j++) {
if (fMinimizer->IsFixedVariable(j)) continue;
fCovFree->SetBinContent(freex + 1, freey + 1, fMinimizer->CovMatrix(i, j));
freey++;
}
freex++;
}
fCorrel = PlotUtils::GetCorrelationPlot(fCovar, "correlation");
fDecomp = PlotUtils::GetDecompPlot(fCovar, "decomposition");
if (fMinimizer->NFree() > 0) {
fCorFree = PlotUtils::GetCorrelationPlot(fCovFree, "correlation_free");
fDecFree = PlotUtils::GetDecompPlot(fCovFree, "decomposition_free");
}
}
+ std::cout << "Got STATE" << std::endl;
return;
};
//*************************************
void MinimizerRoutines::LowStatRoutine(std::string routine) {
//*************************************
LOG(FIT) << "Running Low Statistics Routine: " << routine << std::endl;
int lowstatsevents = FitPar::Config().GetParI("minimizer.lowstatevents");
int maxevents = FitPar::Config().GetParI("input.maxevents");
int verbosity = FitPar::Config().GetParI("VERBOSITY");
std::string trueroutine = routine;
std::string substring = "LowStat";
trueroutine.erase( trueroutine.find(substring),
substring.length() );
// Set MAX EVENTS=1000
FitPar::Config().SetParI("input.maxevents", lowstatsevents);
FitPar::Config().SetParI("VERBOSITY", 3);
SetupFCN();
RunFitRoutine(trueroutine);
FitPar::Config().SetParI("input.maxevents", maxevents);
SetupFCN();
FitPar::Config().SetParI("VERBOSITY", verbosity);
return;
}
//*************************************
void MinimizerRoutines::Create1DScans() {
//*************************************
// 1D Scan Routine
// Steps through all free parameters about nominal using the step size
// Creates a graph for each free parameter
// At the current point create a 1D Scan for all parametes (Uncorrelated)
for (UInt_t i = 0; i < fParams.size(); i++) {
if (fFixVals[fParams[i]]) continue;
LOG(FIT) << "Running 1D Scan for " << fParams[i] << std::endl;
fSampleFCN->CreateIterationTree(fParams[i] +
"_scan1D_iterations",
FitBase::GetRW());
double scanmiddlepoint = fCurVals[fParams[i]];
// Determine N points needed
double limlow = fMinVals[fParams[i]];
double limhigh = fMaxVals[fParams[i]];
double step = fStepVals[fParams[i]];
int npoints = int( fabs(limhigh - limlow) / (step + 0.) );
TH1D* contour = new TH1D(("Chi2Scan1D_" + fParams[i]).c_str(),
("Chi2Scan1D_" + fParams[i] +
";" + fParams[i]).c_str(),
npoints, limlow, limhigh);
// Fill bins
for (int x = 0; x < contour->GetNbinsX(); x++) {
// Set X Val
fCurVals[fParams[i]] = contour->GetXaxis()->GetBinCenter(x + 1);
// Run Eval
double *vals = FitUtils::GetArrayFromMap( fParams, fCurVals );
double chi2 = fSampleFCN->DoEval( vals );
delete vals;
// Fill Contour
contour->SetBinContent(x + 1, chi2);
}
// Save contour
contour->Write();
// Reset Parameter
fCurVals[fParams[i]] = scanmiddlepoint;
// Save TTree
fSampleFCN->WriteIterationTree();
}
return;
}
//*************************************
void MinimizerRoutines::Chi2Scan2D() {
//*************************************
// Chi2 Scan 2D
// Creates a 2D chi2 scan by stepping through all free parameters
// Works for all pairwise combos of free parameters
// Scan I
for (UInt_t i = 0; i < fParams.size(); i++) {
if (fFixVals[fParams[i]]) continue;
// Scan J
for (UInt_t j = 0; j < i; j++) {
if (fFixVals[fParams[j]]) continue;
fSampleFCN->CreateIterationTree( fParams[i] + "_" +
fParams[j] + "_" +
"scan2D_iterations",
FitBase::GetRW() );
double scanmid_i = fCurVals[fParams[i]];
double scanmid_j = fCurVals[fParams[j]];
double limlow_i = fMinVals[fParams[i]];
double limhigh_i = fMaxVals[fParams[i]];
double step_i = fStepVals[fParams[i]];
double limlow_j = fMinVals[fParams[j]];
double limhigh_j = fMaxVals[fParams[j]];
double step_j = fStepVals[fParams[j]];
int npoints_i = int( fabs(limhigh_i - limlow_i) / (step_i + 0.) ) + 1;
int npoints_j = int( fabs(limhigh_j - limlow_j) / (step_j + 0.) ) + 1;
TH2D* contour = new TH2D(("Chi2Scan2D_" + fParams[i] + "_" + fParams[j]).c_str(),
("Chi2Scan2D_" + fParams[i] + "_" + fParams[j] +
";" + fParams[i] + ";" + fParams[j]).c_str(),
npoints_i, limlow_i, limhigh_i,
npoints_j, limlow_j, limhigh_j );
// Begin Scan
LOG(FIT) << "Running scan for " << fParams[i] << " " << fParams[j] << std::endl;
// Fill bins
for (int x = 0; x < contour->GetNbinsX(); x++) {
// Set X Val
fCurVals[fParams[i]] = contour->GetXaxis()->GetBinCenter(x + 1);
// Loop Y
for (int y = 0; y < contour->GetNbinsY(); y++) {
// Set Y Val
fCurVals[fParams[j]] = contour->GetYaxis()->GetBinCenter(y + 1);
// Run Eval
double *vals = FitUtils::GetArrayFromMap( fParams, fCurVals );
double chi2 = fSampleFCN->DoEval( vals );
delete vals;
// Fill Contour
contour->SetBinContent(x + 1, y + 1, chi2);
fCurVals[fParams[j]] = scanmid_j;
}
fCurVals[fParams[i]] = scanmid_i;
fCurVals[fParams[j]] = scanmid_j;
}
// Save contour
contour->Write();
// Save Iterations
fSampleFCN->WriteIterationTree();
}
}
return;
}
//*************************************
void MinimizerRoutines::CreateContours() {
//*************************************
// Use MINUIT for this if possible
ERR(FTL) << " Contours not yet implemented as it is really slow!" << std::endl;
throw;
return;
}
//*************************************
int MinimizerRoutines::FixAtLimit() {
//*************************************
bool fixedparam = false;
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string syst = fParams.at(i);
if (fFixVals[syst]) continue;
double curVal = fCurVals.at(syst);
double minVal = fMinVals.at(syst);
double maxVal = fMinVals.at(syst);
if (fabs(curVal - minVal) < 0.0001) {
fCurVals[syst] = minVal;
fFixVals[syst] = true;
fixedparam = true;
}
if (fabs(maxVal - curVal) < 0.0001) {
fCurVals[syst] = maxVal;
fFixVals[syst] = true;
fixedparam = true;
}
}
if (!fixedparam) {
LOG(FIT) << "No dials needed fixing!" << std::endl;
return kNoChange;
} else return kStateChange;
}
/*
Write Functions
*/
//*************************************
void MinimizerRoutines::SaveResults() {
//*************************************
fOutputRootFile->cd();
if (fMinimizer) {
SetupCovariance();
SaveMinimizerState();
}
SaveCurrentState();
}
//*************************************
void MinimizerRoutines::SaveMinimizerState() {
//*************************************
+ std::cout << "Saving Minimizer State" << std::endl;
if (!fMinimizer) {
ERR(FTL) << "Can't save minimizer state without min object" << std::endl;
throw;
}
// Save main fit tree
fSampleFCN->WriteIterationTree();
-
+
// Get Vals and Errors
GetMinimizerState();
// Save tree with fit status
std::vector<std::string> nameVect;
std::vector<double> valVect;
std::vector<double> errVect;
std::vector<double> minVect;
std::vector<double> maxVect;
std::vector<double> startVect;
std::vector<int> endfixVect;
std::vector<int> startfixVect;
// int NFREEPARS = fMinimizer->NFree();
int NPARS = fMinimizer->NDim();
int ipar = 0;
// Dial Vals
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string name = fParams.at(i);
-
nameVect .push_back( name );
valVect .push_back( fCurVals.at(name) );
+
errVect .push_back( fErrorVals.at(name) );
+
minVect .push_back( fMinVals.at(name) );
+
maxVect .push_back( fMaxVals.at(name) );
startVect .push_back( fStartVals.at(name) );
+
endfixVect .push_back( fFixVals.at(name) );
+
startfixVect.push_back( fStartFixVals.at(name) );
ipar++;
}
int NFREE = fMinimizer->NFree();
int NDIM = fMinimizer->NDim();
double CHI2 = fSampleFCN->GetLikelihood();
int NBINS = fSampleFCN->GetNDOF();
int NDOF = NBINS - NFREE;
// Write fit results
TTree* fit_tree = new TTree("fit_result", "fit_result");
fit_tree->Branch("parameter_names", &nameVect);
fit_tree->Branch("parameter_values", &valVect);
fit_tree->Branch("parameter_errors", &errVect);
fit_tree->Branch("parameter_min", &minVect);
fit_tree->Branch("parameter_max", &maxVect);
fit_tree->Branch("parameter_start", &startVect);
fit_tree->Branch("parameter_fix", &endfixVect);
fit_tree->Branch("parameter_startfix", &startfixVect);
fit_tree->Branch("CHI2", &CHI2, "CHI2/D");
fit_tree->Branch("NDOF", &NDOF, "NDOF/I");
fit_tree->Branch("NBINS", &NBINS, "NBINS/I");
fit_tree->Branch("NDIM", &NDIM, "NDIM/I");
fit_tree->Branch("NFREE", &NFREE, "NFREE/I");
fit_tree->Fill();
fit_tree->Write();
// Make dial variables
TH1D dialvar = TH1D("fit_dials", "fit_dials", NPARS, 0, NPARS);
TH1D startvar = TH1D("start_dials", "start_dials", NPARS, 0, NPARS);
TH1D minvar = TH1D("min_dials", "min_dials", NPARS, 0, NPARS);
TH1D maxvar = TH1D("max_dials", "max_dials", NPARS, 0, NPARS);
TH1D dialvarfree = TH1D("fit_dials_free", "fit_dials_free", NFREE, 0, NFREE);
TH1D startvarfree = TH1D("start_dials_free", "start_dials_free", NFREE, 0, NFREE);
TH1D minvarfree = TH1D("min_dials_free", "min_dials_free", NFREE, 0, NFREE);
TH1D maxvarfree = TH1D("max_dials_free", "max_dials_free", NFREE, 0, NFREE);
int freecount = 0;
for (UInt_t i = 0; i < nameVect.size(); i++) {
std::string name = nameVect.at(i);
dialvar.SetBinContent(i + 1, valVect.at(i));
dialvar.SetBinError(i + 1, errVect.at(i));
dialvar.GetXaxis()->SetBinLabel(i + 1, name.c_str());
startvar.SetBinContent(i + 1, startVect.at(i));
startvar.GetXaxis()->SetBinLabel(i + 1, name.c_str());
minvar.SetBinContent(i + 1, minVect.at(i));
minvar.GetXaxis()->SetBinLabel(i + 1, name.c_str());
maxvar.SetBinContent(i + 1, maxVect.at(i));
maxvar.GetXaxis()->SetBinLabel(i + 1, name.c_str());
if (NFREE > 0) {
if (!startfixVect.at(i)) {
freecount++;
dialvarfree.SetBinContent(freecount, valVect.at(i));
dialvarfree.SetBinError(freecount, errVect.at(i));
dialvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str());
startvarfree.SetBinContent(freecount, startVect.at(i));
startvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str());
minvarfree.SetBinContent(freecount, minVect.at(i));
minvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str());
maxvarfree.SetBinContent(freecount, maxVect.at(i));
maxvarfree.GetXaxis()->SetBinLabel(freecount, name.c_str());
}
}
}
// Save Dial Plots
dialvar.Write();
startvar.Write();
minvar.Write();
maxvar.Write();
if (NFREE > 0) {
dialvarfree.Write();
startvarfree.Write();
minvarfree.Write();
maxvarfree.Write();
}
// Save fit_status plot
TH1D statusplot = TH1D("fit_status", "fit_status", 8, 0, 8);
std::string fit_labels[8] = {"status", "cov_status", \
"maxiter", "maxfunc", \
"iter", "func", \
"precision", "tolerance"
};
double fit_vals[8];
fit_vals[0] = fMinimizer->Status() + 0.;
fit_vals[1] = fMinimizer->CovMatrixStatus() + 0.;
fit_vals[2] = fMinimizer->MaxIterations() + 0.;
fit_vals[3] = fMinimizer->MaxFunctionCalls() + 0.;
fit_vals[4] = fMinimizer->NIterations() + 0.;
fit_vals[5] = fMinimizer->NCalls() + 0.;
fit_vals[6] = fMinimizer->Precision() + 0.;
fit_vals[7] = fMinimizer->Tolerance() + 0.;
for (int i = 0; i < 8; i++) {
statusplot.SetBinContent(i + 1, fit_vals[i]);
statusplot.GetXaxis()->SetBinLabel(i + 1, fit_labels[i].c_str());
}
statusplot.Write();
// Save Covars
if (fCovar) fCovar->Write();
if (fCovFree) fCovFree->Write();
if (fCorrel) fCorrel->Write();
if (fCorFree) fCorFree->Write();
if (fDecomp) fDecomp->Write();
if (fDecFree) fDecFree->Write();
return;
}
//*************************************
void MinimizerRoutines::SaveCurrentState(std::string subdir) {
//*************************************
LOG(FIT) << "Saving current full FCN predictions" << std::endl;
// Setup DIRS
TDirectory* curdir = gDirectory;
if (!subdir.empty()) {
TDirectory* newdir = (TDirectory*) gDirectory->mkdir(subdir.c_str());
newdir->cd();
}
FitBase::GetRW()->Reconfigure();
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->Write();
// Change back to current DIR
curdir->cd();
return;
}
//*************************************
void MinimizerRoutines::SaveNominal() {
//*************************************
fOutputRootFile->cd();
LOG(FIT) << "Saving Nominal Predictions (be cautious with this)" << std::endl;
FitBase::GetRW()->Reconfigure();
SaveCurrentState("nominal");
};
//*************************************
void MinimizerRoutines::SavePrefit() {
//*************************************
fOutputRootFile->cd();
LOG(FIT) << "Saving Prefit Predictions" << std::endl;
UpdateRWEngine(fStartVals);
SaveCurrentState("prefit");
UpdateRWEngine(fCurVals);
};
/*
MISC Functions
*/
//*************************************
int MinimizerRoutines::GetStatus() {
//*************************************
return 0;
}
//*************************************
void MinimizerRoutines::SetupCovariance() {
//*************************************
// Remove covares if they exist
if (fCovar) delete fCovar;
if (fCovFree) delete fCovFree;
if (fCorrel) delete fCorrel;
if (fCorFree) delete fCorFree;
if (fDecomp) delete fDecomp;
if (fDecFree) delete fDecFree;
LOG(FIT) << "Building covariance matrix.." << std::endl;
int NFREE = 0;
int NDIM = 0;
// Get NFREE from min or from vals (for cases when doing throws)
if (fMinimizer) {
+ std::cout << "NFREE FROM MINIMIZER" << std::endl;
NFREE = fMinimizer->NFree();
NDIM = fMinimizer->NDim();
} else {
NDIM = fParams.size();
for (UInt_t i = 0; i < fParams.size(); i++) {
+ std::cout << "Getting Param " << fParams[i] << std::endl;
+
if (!fFixVals[fParams[i]]) NFREE++;
}
}
if (NDIM == 0) return;
LOG(FIT) << "NFREE == " << NFREE << std::endl;
fCovar = new TH2D("covariance", "covariance", NDIM, 0, NDIM, NDIM, 0, NDIM);
if (NFREE > 0) {
fCovFree = new TH2D("covariance_free",
"covariance_free",
NFREE, 0, NFREE,
NFREE, 0, NFREE);
} else {
fCovFree = NULL;
}
// Set Bin Labels
int countall = 0;
int countfree = 0;
for (UInt_t i = 0; i < fParams.size(); i++) {
+ std::cout << "Getting Param " << i << std::endl;
+ std::cout << "ParamI = " << fParams[i] << std::endl;
+
fCovar->GetXaxis()->SetBinLabel(countall + 1, fParams[i].c_str());
fCovar->GetYaxis()->SetBinLabel(countall + 1, fParams[i].c_str());
countall++;
if (!fFixVals[fParams[i]] and NFREE > 0) {
fCovFree->GetXaxis()->SetBinLabel(countfree + 1, fParams[i].c_str());
fCovFree->GetYaxis()->SetBinLabel(countfree + 1, fParams[i].c_str());
countfree++;
}
}
+ std::cout << "Filling Matrices" << std::endl;
+
fCorrel = PlotUtils::GetCorrelationPlot(fCovar, "correlation");
fDecomp = PlotUtils::GetDecompPlot(fCovar, "decomposition");
if (NFREE > 0) {
fCorFree = PlotUtils::GetCorrelationPlot(fCovFree, "correlation_free");
fDecFree = PlotUtils::GetDecompPlot(fCovFree, "decomposition_free");
} else {
fCorFree = NULL;
fDecFree = NULL;
}
+ std::cout << " Set the covariance" << std::endl;
return;
};
//*************************************
void MinimizerRoutines::ThrowCovariance(bool uniformly) {
//*************************************
std::vector<double> rands;
if (!fDecFree) {
ERR(WRN) << "Trying to throw 0 free parameters" << std::endl;
return;
}
// Generate Random Gaussians
for (Int_t i = 0; i < fDecFree->GetNbinsX(); i++) {
rands.push_back(gRandom->Gaus(0.0, 1.0));
}
// Reset Thrown Values
for (UInt_t i = 0; i < fParams.size(); i++) {
fThrownVals[fParams[i]] = fCurVals[fParams[i]];
}
// Loop and get decomp
for (Int_t i = 0; i < fDecFree->GetNbinsX(); i++) {
std::string parname = std::string(fDecFree->GetXaxis()->GetBinLabel(i + 1));
double mod = 0.0;
if (!uniformly) {
for (Int_t j = 0; j < fDecFree->GetNbinsY(); j++) {
mod += rands[j] * fDecFree->GetBinContent(j + 1, i + 1);
}
}
if (fCurVals.find(parname) != fCurVals.end()) {
if (uniformly) fThrownVals[parname] = gRandom->Uniform(fMinVals[parname], fMaxVals[parname]);
else { fThrownVals[parname] = fCurVals[parname] + mod; }
}
}
// Check Limits
for (UInt_t i = 0; i < fParams.size(); i++) {
std::string syst = fParams[i];
if (fFixVals[syst]) continue;
if (fThrownVals[syst] < fMinVals[syst]) fThrownVals[syst] = fMinVals[syst];
if (fThrownVals[syst] > fMaxVals[syst]) fThrownVals[syst] = fMaxVals[syst];
}
return;
};
//*************************************
void MinimizerRoutines::GenerateErrorBands() {
//*************************************
TDirectory* errorDIR = (TDirectory*) fOutputRootFile->mkdir("error_bands");
errorDIR->cd();
// Make a second file to store throws
std::string tempFileName = fOutputFile;
if (tempFileName.find(".root") != std::string::npos) tempFileName.erase(tempFileName.find(".root"), 5);
tempFileName += ".throws.root";
TFile* tempfile = new TFile(tempFileName.c_str(),"RECREATE");
tempfile->cd();
int nthrows = FitPar::Config().GetParI("error_throws");
UpdateRWEngine(fCurVals);
fSampleFCN->ReconfigureAllEvents();
TDirectory* nominal = (TDirectory*) tempfile->mkdir("nominal");
nominal->cd();
fSampleFCN->Write();
TDirectory* outnominal = (TDirectory*) fOutputRootFile->mkdir("nominal_throw");
outnominal->cd();
fSampleFCN->Write();
errorDIR->cd();
TTree* parameterTree = new TTree("throws", "throws");
double chi2;
for (UInt_t i = 0; i < fParams.size(); i++)
parameterTree->Branch(fParams[i].c_str(), &fThrownVals[fParams[i]], (fParams[i] + "/D").c_str());
parameterTree->Branch("chi2", &chi2, "chi2/D");
bool uniformly = FitPar::Config().GetParB("error_uniform");
// Run Throws and save
for (Int_t i = 0; i < nthrows; i++) {
TDirectory* throwfolder = (TDirectory*)tempfile->mkdir(Form("throw_%i", i));
throwfolder->cd();
// Generate Random Parameter Throw
ThrowCovariance(uniformly);
// Run Eval
double *vals = FitUtils::GetArrayFromMap( fParams, fThrownVals );
chi2 = fSampleFCN->DoEval( vals );
delete vals;
// Save the FCN
fSampleFCN->Write();
parameterTree->Fill();
}
errorDIR->cd();
fDecFree->Write();
fCovFree->Write();
parameterTree->Write();
delete parameterTree;
// Now go through the keys in the temporary file and look for TH1D, and TH2D plots
TIter next(nominal->GetListOfKeys());
TKey *key;
while ((key = (TKey*)next())) {
TClass *cl = gROOT->GetClass(key->GetClassName());
if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue;
TH1D *baseplot = (TH1D*)key->ReadObj();
std::string plotname = std::string(baseplot->GetName());
int nbins = baseplot->GetNbinsX() * baseplot->GetNbinsY();
// Setup TProfile with RMS option
TProfile* tprof = new TProfile((plotname + "_prof").c_str(), (plotname + "_prof").c_str(), nbins, 0, nbins, "S");
// Setup The TTREE
double* bincontents;
bincontents = new double[nbins];
double* binlowest;
binlowest = new double[nbins];
double* binhighest;
binhighest = new double[nbins];
errorDIR->cd();
TTree* bintree = new TTree((plotname + "_tree").c_str(), (plotname + "_tree").c_str());
for (Int_t i = 0; i < nbins; i++) {
bincontents[i] = 0.0;
binhighest[i] = 0.0;
binlowest[i] = 0.0;
bintree->Branch(Form("content_%i", i), &bincontents[i], Form("content_%i/D", i));
}
for (Int_t i = 0; i < nthrows; i++) {
TH1* newplot = (TH1*)tempfile->Get(Form(("throw_%i/" + plotname).c_str(), i));
for (Int_t j = 0; j < nbins; j++) {
tprof->Fill(j + 0.5, newplot->GetBinContent(j + 1));
bincontents[j] = newplot->GetBinContent(j + 1);
if (bincontents[j] < binlowest[j] or i == 0) binlowest[j] = bincontents[j];
if (bincontents[j] > binhighest[j] or i == 0) binhighest[j] = bincontents[j];
}
errorDIR->cd();
bintree->Fill();
delete newplot;
}
errorDIR->cd();
for (Int_t j = 0; j < nbins; j++) {
if (!uniformly) {
baseplot->SetBinError(j + 1, tprof->GetBinError(j + 1));
} else {
baseplot->SetBinContent(j + 1, (binlowest[j] + binhighest[j]) / 2.0);
baseplot->SetBinError(j + 1, (binhighest[j] - binlowest[j]) / 2.0);
}
}
errorDIR->cd();
baseplot->Write();
tprof->Write();
bintree->Write();
delete baseplot;
delete tprof;
delete bintree;
delete [] bincontents;
}
return;
};
diff --git a/src/Routines/SystematicRoutines.cxx b/src/Routines/SystematicRoutines.cxx
index 124eb98..a09213a 100755
--- a/src/Routines/SystematicRoutines.cxx
+++ b/src/Routines/SystematicRoutines.cxx
@@ -1,1508 +1,1508 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "SystematicRoutines.h"
void SystematicRoutines::Init(){
fInputFile = "";
fInputRootFile = NULL;
fOutputFile = "";
fOutputRootFile = NULL;
fCovar = fCovarFree = NULL;
fCorrel = fCorrelFree = NULL;
fDecomp = fDecompFree = NULL;
fStrategy = "ErrorBands";
fRoutines.clear();
fRoutines.push_back("ErrorBands");
fCardFile = "";
fFakeDataInput = "";
fSampleFCN = NULL;
fAllowedRoutines = ("ErrorBands,PlotLimits");
};
SystematicRoutines::~SystematicRoutines(){
};
SystematicRoutines::SystematicRoutines(int argc, char* argv[]){
// Initialise Defaults
Init();
nuisconfig configuration = Config::Get();
// Default containers
std::string cardfile = "";
std::string maxevents = "-1";
int errorcount = 0;
int verbocount = 0;
std::vector<std::string> xmlcmds;
std::vector<std::string> configargs;
fNThrows = 250;
fStartThrows = 0;
fThrowString = "";
// Make easier to handle arguments.
std::vector<std::string> args = GeneralUtils::LoadCharToVectStr(argc, argv);
ParserUtils::ParseArgument(args, "-c", fCardFile, true);
ParserUtils::ParseArgument(args, "-o", fOutputFile, false, false);
ParserUtils::ParseArgument(args, "-n", maxevents, false, false);
ParserUtils::ParseArgument(args, "-f", fStrategy, false, false);
ParserUtils::ParseArgument(args, "-d", fFakeDataInput, false, false);
ParserUtils::ParseArgument(args, "-s", fStartThrows, false, false);
ParserUtils::ParseArgument(args, "-t", fNThrows, false, false);
ParserUtils::ParseArgument(args, "-p", fThrowString, false, false);
ParserUtils::ParseArgument(args, "-i", xmlcmds);
ParserUtils::ParseArgument(args, "-q", configargs);
ParserUtils::ParseCounter(args, "e", errorcount);
ParserUtils::ParseCounter(args, "v", verbocount);
ParserUtils::CheckBadArguments(args);
// Add extra defaults if none given
if (fCardFile.empty() and xmlcmds.empty()) {
ERR(FTL) << "No input supplied!" << std::endl;
throw;
}
if (fOutputFile.empty() and !fCardFile.empty()) {
fOutputFile = fCardFile + ".root";
ERR(WRN) << "No output supplied so saving it to: " << fOutputFile << std::endl;
} else if (fOutputFile.empty()) {
ERR(FTL) << "No output file or cardfile supplied!" << std::endl;
throw;
}
// Configuration Setup =============================
// Check no comp key is available
if (Config::Get().GetNodes("nuiscomp").empty()) {
fCompKey = Config::Get().CreateNode("nuiscomp");
} else {
fCompKey = Config::Get().GetNodes("nuiscomp")[0];
}
if (!fCardFile.empty()) fCompKey.AddS("cardfile", fCardFile);
if (!fOutputFile.empty()) fCompKey.AddS("outputfile", fOutputFile);
if (!fStrategy.empty()) fCompKey.AddS("strategy", fStrategy);
// Load XML Cardfile
configuration.LoadConfig( fCompKey.GetS("cardfile"), "");
// Add CMD XML Structs
for (size_t i = 0; i < xmlcmds.size(); i++) {
configuration.AddXMLLine(xmlcmds[i]);
}
// Add Config Args
for (size_t i = 0; i < configargs.size(); i++) {
configuration.OverrideConfig(configargs[i]);
}
if (maxevents.compare("-1")){
configuration.OverrideConfig("MAXEVENTS=" + maxevents);
}
// Finish configuration XML
configuration.FinaliseConfig(fCompKey.GetS("outputfile") + ".xml");
// Add Error Verbo Lines
verbocount += Config::Get().GetParI("VERBOSITY");
errorcount += Config::Get().GetParI("ERROR");
std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl;
std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl;
SETVERBOSITY(verbocount);
// Proper Setup
if (fStrategy.find("ErrorBands") != std::string::npos ||
fStrategy.find("MergeErrors") != std::string::npos){
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
}
// fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
SetupSystematicsFromXML();
SetupCovariance();
SetupRWEngine();
SetupFCN();
GetCovarFromFCN();
// Run();
return;
};
void SystematicRoutines::SetupSystematicsFromXML(){
LOG(FIT) << "Setting up nuismin" << std::endl;
// Setup Parameters ------------------------------------------
std::vector<nuiskey> parkeys = Config::QueryKeys("parameter");
if (!parkeys.empty()) {
LOG(FIT) << "Number of parameters : " << parkeys.size() << std::endl;
}
for (size_t i = 0; i < parkeys.size(); i++) {
nuiskey key = parkeys.at(i);
// Check for type,name,nom
if (!key.Has("type")) {
ERR(FTL) << "No type given for parameter " << i << std::endl;
throw;
} else if (!key.Has("name")) {
ERR(FTL) << "No name given for parameter " << i << std::endl;
throw;
} else if (!key.Has("nominal")) {
ERR(FTL) << "No nominal given for parameter " << i << std::endl;
throw;
}
// Get Inputs
std::string partype = key.GetS("type");
std::string parname = key.GetS("name");
double parnom = key.GetD("nominal");
double parlow = parnom - 1;
double parhigh = parnom + 1;
double parstep = 1;
std::string parstate = key.GetS("state");
// Extra limits
if (key.Has("low")) {
parlow = key.GetD("low");
parhigh = key.GetD("high");
parstep = key.GetD("step");
LOG(FIT) << "Read " << partype << " : "
<< parname << " = "
<< parnom << " : "
<< parlow << " < p < " << parhigh
<< " : " << parstate << std::endl;
} else {
LOG(FIT) << "Read " << partype << " : "
<< parname << " = "
<< parnom << " : "
<< parstate << std::endl;
}
// Run Parameter Conversion if needed
if (parstate.find("ABS") != std::string::npos) {
parnom = FitBase::RWAbsToSigma( partype, parname, parnom );
parlow = FitBase::RWAbsToSigma( partype, parname, parlow );
parhigh = FitBase::RWAbsToSigma( partype, parname, parhigh );
parstep = FitBase::RWAbsToSigma( partype, parname, parstep );
} else if (parstate.find("FRAC") != std::string::npos) {
parnom = FitBase::RWFracToSigma( partype, parname, parnom );
parlow = FitBase::RWFracToSigma( partype, parname, parlow );
parhigh = FitBase::RWFracToSigma( partype, parname, parhigh );
parstep = FitBase::RWFracToSigma( partype, parname, parstep );
}
// Push into vectors
fParams.push_back(parname);
fTypeVals[parname] = FitBase::ConvDialType(partype);;
fStartVals[parname] = parnom;
fCurVals[parname] = parnom;
fErrorVals[parname] = 0.0;
fStateVals[parname] = parstate;
bool fixstate = parstate.find("FIX") != std::string::npos;
fFixVals[parname] = fixstate;
fStartFixVals[parname] = fFixVals[parname];
fMinVals[parname] = parlow;
fMaxVals[parname] = parhigh;
fStepVals[parname] = parstep;
}
// Setup Samples ----------------------------------------------
std::vector<nuiskey> samplekeys = Config::QueryKeys("sample");
if (!samplekeys.empty()) {
LOG(FIT) << "Number of samples : " << samplekeys.size() << std::endl;
}
for (size_t i = 0; i < samplekeys.size(); i++) {
nuiskey key = samplekeys.at(i);
// Get Sample Options
std::string samplename = key.GetS("name");
std::string samplefile = key.GetS("input");
std::string sampletype =
key.Has("type") ? key.GetS("type") : "DEFAULT";
double samplenorm =
key.Has("norm") ? key.GetD("norm") : 1.0;
// Print out
LOG(FIT) << "Read sample info " << i << " : "
<< samplename << std::endl
<< "\t\t input -> " << samplefile << std::endl
<< "\t\t state -> " << sampletype << std::endl
<< "\t\t norm -> " << samplenorm << std::endl;
// If FREE add to parameters otherwise continue
if (sampletype.find("FREE") == std::string::npos) {
continue;
}
// Form norm dial from samplename + sampletype + "_norm";
- std::string normname = samplename + sampletype + "_norm";
+ std::string normname = samplename + "_norm";
// Check normname not already present
if (fTypeVals.find(normname) != fTypeVals.end()) {
continue;
}
// Add new norm dial to list if its passed above checks
fParams.push_back(normname);
fTypeVals[normname] = kNORM;
fStateVals[normname] = sampletype;
fCurVals[normname] = samplenorm;
fErrorVals[normname] = 0.0;
fMinVals[normname] = 0.1;
fMaxVals[normname] = 10.0;
fStepVals[normname] = 0.5;
bool state = sampletype.find("FREE") == std::string::npos;
fFixVals[normname] = state;
fStartFixVals[normname] = state;
}
// Setup Fake Parameters -----------------------------
std::vector<nuiskey> fakekeys = Config::QueryKeys("fakeparameter");
if (!fakekeys.empty()) {
LOG(FIT) << "Number of fake parameters : " << fakekeys.size() << std::endl;
}
for (size_t i = 0; i < fakekeys.size(); i++) {
nuiskey key = fakekeys.at(i);
// Check for type,name,nom
if (!key.Has("name")) {
ERR(FTL) << "No name given for fakeparameter " << i << std::endl;
throw;
} else if (!key.Has("nom")) {
ERR(FTL) << "No nominal given for fakeparameter " << i << std::endl;
throw;
}
// Get Inputs
std::string parname = key.GetS("name");
double parnom = key.GetD("nom");
// Push into vectors
fFakeVals[parname] = parnom;
}
}
void SystematicRoutines::ReadCard(std::string cardfile){
// Read cardlines into vector
std::vector<std::string> cardlines = GeneralUtils::ParseFileToStr(cardfile,"\n");
FitPar::Config().cardLines = cardlines;
// Read Samples first (norm params can be overridden)
int linecount = 0;
for (std::vector<std::string>::iterator iter = cardlines.begin();
iter != cardlines.end(); iter++){
std::string line = (*iter);
linecount++;
// Skip Empties
if (line.empty()) continue;
if (line.c_str()[0] == '#') continue;
// Read Valid Samples
int samstatus = ReadSamples(line);
// Show line if bad to help user
if (samstatus == kErrorStatus) {
ERR(FTL) << "Bad Input in cardfile " << fCardFile
<< " at line " << linecount << "!" << std::endl;
LOG(FIT) << line << std::endl;
throw;
}
}
// Read Parameters second
linecount = 0;
for (std::vector<std::string>::iterator iter = cardlines.begin();
iter != cardlines.end(); iter++){
std::string line = (*iter);
linecount++;
// Skip Empties
if (line.empty()) continue;
if (line.c_str()[0] == '#') continue;
// Try Parameter Reads
int parstatus = ReadParameters(line);
int fakstatus = ReadFakeDataPars(line);
// Show line if bad to help user
if (parstatus == kErrorStatus ||
fakstatus == kErrorStatus ){
ERR(FTL) << "Bad Parameter Input in cardfile " << fCardFile
<< " at line " << linecount << "!" << std::endl;
LOG(FIT) << line << std::endl;
throw;
}
}
return;
};
int SystematicRoutines::ReadParameters(std::string parstring){
std::string inputspec = "RW Dial Inputs Syntax \n"
"free input w/ limits: TYPE NAME START MIN MAX STEP [STATE] \n"
"fix input: TYPE NAME VALUE [STATE] \n"
"free input w/o limits: TYPE NAME START FREE,[STATE] \n"
"Allowed Types: \n"
"neut_parameter,niwg_parameter,t2k_parameter,"
"nuwro_parameter,gibuu_parameter";
// Check sample input
if (parstring.find("parameter") == std::string::npos) return kGoodStatus;
// Parse inputs
std::vector<std::string> strvct = GeneralUtils::ParseToStr(parstring, " ");
// Skip if comment or parameter somewhere later in line
if (strvct[0].c_str()[0] == '#' ||
strvct[0].find("parameter") == std::string::npos){
return kGoodStatus;
}
// Check length
if (strvct.size() < 3){
ERR(FTL) << "Input rw dials need to provide at least 3 inputs." << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Setup default inputs
std::string partype = strvct[0];
std::string parname = strvct[1];
double parval = GeneralUtils::StrToDbl(strvct[2]);
double minval = parval - 1.0;
double maxval = parval + 1.0;
double stepval = 1.0;
std::string state = "FIX"; //[DEFAULT]
// Check Type
if (FitBase::ConvDialType(partype) == kUNKNOWN){
ERR(FTL) << "Unknown parameter type! " << partype << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Check Parameter Name
if (FitBase::GetDialEnum(partype, parname) == -1){
ERR(FTL) << "Bad RW parameter name! " << partype << " " << parname << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Option Extra (No Limits)
if (strvct.size() == 4){
state = strvct[3];
}
// Check for weirder inputs
if (strvct.size() > 4 && strvct.size() < 6){
ERR(FTL) << "Provided incomplete limits for " << parname << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Option Extra (With limits and steps)
if (strvct.size() >= 6){
minval = GeneralUtils::StrToDbl(strvct[3]);
maxval = GeneralUtils::StrToDbl(strvct[4]);
stepval = GeneralUtils::StrToDbl(strvct[5]);
}
// Option Extra (dial state after limits)
if (strvct.size() == 7){
state = strvct[6];
}
// Run Parameter Conversion if needed
if (state.find("ABS") != std::string::npos){
parval = FitBase::RWAbsToSigma( partype, parname, parval );
minval = FitBase::RWAbsToSigma( partype, parname, minval );
maxval = FitBase::RWAbsToSigma( partype, parname, maxval );
stepval = FitBase::RWAbsToSigma( partype, parname, stepval );
} else if (state.find("FRAC") != std::string::npos){
parval = FitBase::RWFracToSigma( partype, parname, parval );
minval = FitBase::RWFracToSigma( partype, parname, minval );
maxval = FitBase::RWFracToSigma( partype, parname, maxval );
stepval = FitBase::RWFracToSigma( partype, parname, stepval );
}
// Check no repeat params
if (std::find(fParams.begin(), fParams.end(), parname) != fParams.end()){
ERR(FTL) << "Duplicate parameter names given for " << parname << std::endl;
throw;
}
// Setup Containers
fParams.push_back(parname);
fTypeVals[parname] = FitBase::ConvDialType(partype);
fStartVals[parname] = parval;
fCurVals[parname] = fStartVals[parname];
fErrorVals[parname] = 0.0;
fStateVals[parname] = state;
bool fixstate = state.find("FIX") != std::string::npos;
fFixVals[parname] = fixstate;
fStartFixVals[parname] = fFixVals[parname];
fMinVals[parname] = minval;
fMaxVals[parname] = maxval;
fStepVals[parname] = stepval;
// Print the parameter
LOG(MIN) << "Read Parameter " << parname << " " << parval << " "
<< minval << " " << maxval << " "
<< stepval << " " << state << std::endl;
// Tell reader its all good
return kGoodStatus;
}
//*******************************************
int SystematicRoutines::ReadFakeDataPars(std::string parstring){
//******************************************
std::string inputspec = "Fake Data Dial Inputs Syntax \n"
"fake value: fake_parameter NAME VALUE \n"
"Name should match dialnames given in actual dial specification.";
// Check sample input
if (parstring.find("fake_parameter") == std::string::npos)
return kGoodStatus;
// Parse inputs
std::vector<std::string> strvct = GeneralUtils::ParseToStr(parstring, " ");
// Skip if comment or parameter somewhere later in line
if (strvct[0].c_str()[0] == '#' ||
strvct[0] == "fake_parameter"){
return kGoodStatus;
}
// Check length
if (strvct.size() < 3){
ERR(FTL) << "Fake dials need to provide at least 3 inputs." << std::endl;
std::cout << inputspec << std::endl;
return kErrorStatus;
}
// Read Inputs
std::string parname = strvct[1];
double parval = GeneralUtils::StrToDbl(strvct[2]);
// Setup Container
fFakeVals[parname] = parval;
// Print the fake parameter
LOG(MIN) << "Read Fake Parameter " << parname << " " << parval << std::endl;
// Tell reader its all good
return kGoodStatus;
}
//******************************************
int SystematicRoutines::ReadSamples(std::string samstring){
//******************************************
const static std::string inputspec =
"\tsample <sample_name> <input_type>:inputfile.root [OPTS] "
"[norm]\nsample_name: Name "
"of sample to include. e.g. MiniBooNE_CCQE_XSec_1DQ2_nu\ninput_type: The "
"input event format. e.g. NEUT, GENIE, EVSPLN, ...\nOPTS: Additional, "
"optional sample options.\nnorm: Additional, optional sample "
"normalisation factor.";
// Check sample input
if (samstring.find("sample") == std::string::npos)
return kGoodStatus;
// Parse inputs
std::vector<std::string> strvct = GeneralUtils::ParseToStr(samstring, " ");
// Skip if comment or parameter somewhere later in line
if (strvct[0].c_str()[0] == '#' ||
strvct[0] != "sample"){
return kGoodStatus;
}
// Check length
if (strvct.size() < 3){
ERR(FTL) << "Sample need to provide at least 3 inputs." << std::endl;
return kErrorStatus;
}
// Setup default inputs
std::string samname = strvct[1];
std::string samfile = strvct[2];
if (samfile == "FIX") {
ERR(FTL) << "Input filename was \"FIX\", this line is probably malformed "
"in the input card file. Line:\'"
<< samstring << "\'" << std::endl;
ERR(FTL) << "Expect sample lines to look like:\n\t" << inputspec
<< std::endl;
throw;
}
std::string samtype = "DEFAULT";
double samnorm = 1.0;
// Optional Type
if (strvct.size() > 3) {
samtype = strvct[3];
// samname += "_"+samtype;
// Also get rid of the / and replace it with underscore because it might not be supported character
// while (samname.find("/") != std::string::npos) {
// samname.replace(samname.find("/"), 1, std::string("_"));
// }
}
// Optional Norm
if (strvct.size() > 4) samnorm = GeneralUtils::StrToDbl(strvct[4]);
// Add Sample Names as Norm Dials
std::string normname = samname + "_norm";
// Check no repeat params
if (std::find(fParams.begin(), fParams.end(), normname) != fParams.end()){
ERR(FTL) << "Duplicate samples given for " << samname << std::endl;
throw;
}
fParams.push_back(normname);
fTypeVals[normname] = kNORM;
fStartVals[normname] = samnorm;
fCurVals[normname] = fStartVals[normname];
fErrorVals[normname] = 0.0;
fMinVals[normname] = 0.1;
fMaxVals[normname] = 10.0;
fStepVals[normname] = 0.5;
bool state = samtype.find("FREE") == std::string::npos;
fFixVals[normname] = state;
fStartFixVals[normname] = state;
// Print read in
LOG(MIN) << "Read sample " << samname << " "
<< samfile << " " << samtype << " "
<< samnorm << std::endl;
// Tell reader its all good
return kGoodStatus;
}
/*
Setup Functions
*/
//*************************************
void SystematicRoutines::SetupRWEngine(){
//*************************************
for (UInt_t i = 0; i < fParams.size(); i++){
std::string name = fParams[i];
FitBase::GetRW() -> IncludeDial(name, fTypeVals.at(name) );
}
UpdateRWEngine(fStartVals);
return;
}
//*************************************
void SystematicRoutines::SetupFCN(){
//*************************************
LOG(FIT)<<"Making the jointFCN"<<std::endl;
if (fSampleFCN) delete fSampleFCN;
fSampleFCN = new JointFCN(fOutputRootFile);
SetFakeData();
return;
}
//*************************************
void SystematicRoutines::SetFakeData(){
//*************************************
if (fFakeDataInput.empty()) return;
if (fFakeDataInput.compare("MC") == 0){
LOG(FIT)<<"Setting fake data from MC starting prediction." <<std::endl;
UpdateRWEngine(fFakeVals);
FitBase::GetRW()->Reconfigure();
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->SetFakeData("MC");
UpdateRWEngine(fCurVals);
LOG(FIT)<<"Set all data to fake MC predictions."<<std::endl;
} else {
fSampleFCN->SetFakeData(fFakeDataInput);
}
return;
}
//*****************************************
void SystematicRoutines::GetCovarFromFCN(){
//*****************************************
LOG(FIT) << "Loading ParamPull objects from FCN to build covar" << std::endl;
// Make helperstring
std::ostringstream helperstr;
// Keep track of what is being thrown
std::map<std::string, std::string> dialthrowhandle;
// Get Covariance Objects from FCN
std::list<ParamPull*> inputpulls = fSampleFCN->GetPullList();
for (PullListConstIter iter = inputpulls.begin();
iter != inputpulls.end(); iter++){
ParamPull* pull = (*iter);
if (pull->GetType().find("THROW")){
fInputThrows.push_back(pull);
fInputCovar.push_back(pull->GetFullCovarMatrix());
fInputDials.push_back(pull->GetDataHist());
LOG(FIT) << "Read ParamPull: " << pull->GetName() << " " << pull->GetType() << std::endl;
}
TH1D dialhist = pull->GetDataHist();
TH1D minhist = pull->GetMinHist();
TH1D maxhist = pull->GetMaxHist();
TH1I typehist = pull->GetDialTypes();
for (int i = 0; i < dialhist.GetNbinsX(); i++){
std::string name = std::string(dialhist.GetXaxis()->GetBinLabel(i+1));
dialthrowhandle[name] = pull->GetName();
if (fCurVals.find(name) == fCurVals.end()){
// Add to Containers
fParams.push_back(name);
fCurVals[name] = dialhist.GetBinContent(i+1);
fStartVals[name] = dialhist.GetBinContent(i+1);
fMinVals[name] = minhist.GetBinContent(i+1);
fMaxVals[name] = maxhist.GetBinContent(i+1);
fStepVals[name] = 1.0;
fFixVals[name] = false;
fStartFixVals[name] = false;
fTypeVals[name] = typehist.GetBinContent(i+1);
fStateVals[name] = "FREE" + pull->GetType();
// Maker Helper
helperstr << std::string(16, ' ' ) << FitBase::ConvDialType(fTypeVals[name]) << " "
<< name << " " << fMinVals[name] << " "
<< fMaxVals[name] << " " << fStepVals[name] << " " << fStateVals[name]
<< std::endl;
}
}
}
// Check if no throws given
if (fInputThrows.empty()){
ERR(WRN) << "No covariances given to nuissyst" << std::endl;
ERR(WRN) << "Pushing back an uncorrelated gaussian throw error for each free parameter using step size" << std::endl;
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams[i];
if (fFixVals[syst]) continue;
// Make Terms
std::string name = syst + "_pull";
std::ostringstream pullterm;
pullterm << "DIAL:" << syst << ";"
<< fStartVals[syst] << ";"
<< fStepVals[syst];
std::string type = "GAUSTHROW/NEUT";
// Push Back Pulls
ParamPull* pull = new ParamPull( name, pullterm.str(), type );
fInputThrows.push_back(pull);
fInputCovar.push_back(pull->GetFullCovarMatrix());
fInputDials.push_back(pull->GetDataHist());
// Print Whats added
ERR(WRN) << "Added ParamPull : " << name << " " << pullterm.str() << " " << type << std::endl;
// Add helper string for future fits
helperstr << std::string(16, ' ' ) << "covar " << name << " " << pullterm.str() << " " << type << std::endl;
// Keep Track of Throws
dialthrowhandle[syst] = pull->GetName();
}
}
// Print Helper String
if (!helperstr.str().empty()){
LOG(FIT) << "To remove these statements in future studies, add the lines below to your card:" << std::endl;
// Can't use the logger properly because this can be multi-line. Use cout and added spaces to look better!
std::cout << helperstr.str();
sleep(2);
}
// Print Throw State
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams[i];
if (dialthrowhandle.find(syst) != dialthrowhandle.end()){
LOG(FIT) << "Dial " << i << ". " << setw(40) << syst << " = THROWING with " << dialthrowhandle[syst] << std::endl;
} else {
LOG(FIT) << "Dial " << i << ". " << setw(40) << syst << " = FIXED" << std::endl;
}
}
// Pause anyway
sleep(1);
return;
}
/*
Fitting Functions
*/
//*************************************
void SystematicRoutines::UpdateRWEngine(std::map<std::string,double>& updateVals){
//*************************************
for (UInt_t i = 0; i < fParams.size(); i++){
std::string name = fParams[i];
if (updateVals.find(name) == updateVals.end()) continue;
FitBase::GetRW()->SetDialValue(name,updateVals.at(name));
}
FitBase::GetRW()->Reconfigure();
return;
}
//*************************************
void SystematicRoutines::PrintState(){
//*************************************
LOG(FIT)<<"------------"<<std::endl;
// Count max size
int maxcount = 0;
for (UInt_t i = 0; i < fParams.size(); i++){
maxcount = max(int(fParams[i].size()), maxcount);
}
// Header
LOG(FIT) << " # " << left << setw(maxcount) << "Parameter "
<< " = "
<< setw(10) << "Value" << " +- "
<< setw(10) << "Error" << " "
<< setw(8) << "(Units)" << " "
<< setw(10) << "Conv. Val" << " +- "
<< setw(10) << "Conv. Err" << " "
<< setw(8) << "(Units)" << std::endl;
// Parameters
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams.at(i);
std::string typestr = FitBase::ConvDialType(fTypeVals[syst]);
std::string curunits = "(sig.)";
double curval = fCurVals[syst];
double curerr = fErrorVals[syst];
if (fStateVals[syst].find("ABS") != std::string::npos){
curval = FitBase::RWSigmaToAbs(typestr, syst, curval);
curerr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
FitBase::RWSigmaToAbs(typestr, syst, 0.0));
curunits = "(Abs.)";
} else if (fStateVals[syst].find("FRAC") != std::string::npos){
curval = FitBase::RWSigmaToFrac(typestr, syst, curval);
curerr = (FitBase::RWSigmaToFrac(typestr, syst, curerr) -
FitBase::RWSigmaToFrac(typestr, syst, 0.0));
curunits = "(Frac)";
}
std::string convunits = "(" + FitBase::GetRWUnits(typestr, syst) + ")";
double convval = FitBase::RWSigmaToAbs(typestr, syst, curval);
double converr = (FitBase::RWSigmaToAbs(typestr, syst, curerr) -
FitBase::RWSigmaToAbs(typestr, syst, 0.0));
std::ostringstream curparstring;
curparstring << " " << setw(3) << left
<< i << ". "
<< setw(maxcount) << syst << " = "
<< setw(10) << curval << " +- "
<< setw(10) << curerr << " "
<< setw(8) << curunits << " "
<< setw(10) << convval << " +- "
<< setw(10) << converr << " "
<< setw(8) << convunits;
LOG(FIT) << curparstring.str() << std::endl;
}
LOG(FIT)<<"------------"<<std::endl;
double like = fSampleFCN->GetLikelihood();
LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like << std::endl;
LOG(FIT)<<"------------"<<std::endl;
}
/*
Write Functions
*/
//*************************************
void SystematicRoutines::SaveResults(){
//*************************************
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
SaveCurrentState();
}
//*************************************
void SystematicRoutines::SaveCurrentState(std::string subdir){
//*************************************
LOG(FIT)<<"Saving current full FCN predictions" <<std::endl;
// Setup DIRS
TDirectory* curdir = gDirectory;
if (!subdir.empty()){
TDirectory* newdir =(TDirectory*) gDirectory->mkdir(subdir.c_str());
newdir->cd();
}
FitBase::GetRW()->Reconfigure();
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->Write();
// Change back to current DIR
curdir->cd();
return;
}
//*************************************
void SystematicRoutines::SaveNominal(){
//*************************************
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
LOG(FIT)<<"Saving Nominal Predictions (be cautious with this)" <<std::endl;
FitBase::GetRW()->Reconfigure();
SaveCurrentState("nominal");
};
//*************************************
void SystematicRoutines::SavePrefit(){
//*************************************
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
LOG(FIT)<<"Saving Prefit Predictions"<<std::endl;
UpdateRWEngine(fStartVals);
SaveCurrentState("prefit");
UpdateRWEngine(fCurVals);
};
/*
MISC Functions
*/
//*************************************
int SystematicRoutines::GetStatus(){
//*************************************
return 0;
}
//*************************************
void SystematicRoutines::SetupCovariance(){
//*************************************
// Remove covares if they exist
if (fCovar) delete fCovar;
if (fCovarFree) delete fCovarFree;
if (fCorrel) delete fCorrel;
if (fCorrelFree) delete fCorrelFree;
if (fDecomp) delete fDecomp;
if (fDecompFree) delete fDecompFree;
int NFREE = 0;
int NDIM = 0;
// Get NFREE from min or from vals (for cases when doing throws)
NDIM = fParams.size();
for (UInt_t i = 0; i < fParams.size(); i++){
if (!fFixVals[fParams[i]]) NFREE++;
}
if (NDIM == 0) return;
fCovar = new TH2D("covariance","covariance",NDIM,0,NDIM,NDIM,0,NDIM);
if (NFREE > 0){
fCovarFree = new TH2D("covariance_free",
"covariance_free",
NFREE,0,NFREE,
NFREE,0,NFREE);
}
// Set Bin Labels
int countall = 0;
int countfree = 0;
for (UInt_t i = 0; i < fParams.size(); i++){
fCovar->GetXaxis()->SetBinLabel(countall+1,fParams[i].c_str());
fCovar->GetYaxis()->SetBinLabel(countall+1,fParams[i].c_str());
countall++;
if (!fFixVals[fParams[i]] and NFREE > 0){
fCovarFree->GetXaxis()->SetBinLabel(countfree+1,fParams[i].c_str());
fCovarFree->GetYaxis()->SetBinLabel(countfree+1,fParams[i].c_str());
countfree++;
}
}
fCorrel = PlotUtils::GetCorrelationPlot(fCovar,"correlation");
fDecomp = PlotUtils::GetDecompPlot(fCovar,"decomposition");
if (NFREE > 0)fCorrelFree = PlotUtils::GetCorrelationPlot(fCovarFree, "correlation_free");
if (NFREE > 0)fDecompFree = PlotUtils::GetDecompPlot(fCovarFree,"decomposition_free");
return;
};
//*************************************
void SystematicRoutines::ThrowCovariance(bool uniformly){
//*************************************
// Set fThrownVals to all values in currentVals
for (UInt_t i = 0; i < fParams.size(); i++){
std::string name = fParams.at(i);
fThrownVals[name] = fCurVals[name];
}
for (PullListConstIter iter = fInputThrows.begin();
iter != fInputThrows.end(); iter++){
ParamPull* pull = *iter;
pull->ThrowCovariance();
TH1D dialhist = pull->GetDataHist();
for (int i = 0; i < dialhist.GetNbinsX(); i++){
std::string name = std::string(dialhist.GetXaxis()->GetBinLabel(i+1));
if (fCurVals.find(name) != fCurVals.end()){
fThrownVals[name] = dialhist.GetBinContent(i+1);
}
}
// Reset throw incase pulls are calculated.
pull->ResetToy();
}
return;
};
//*************************************
void SystematicRoutines::PlotLimits(){
//*************************************
std::cout << "Plotting Limits" << std::endl;
if (!fOutputRootFile)
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
TDirectory* limfolder = (TDirectory*) fOutputRootFile->mkdir("Limits");
limfolder->cd();
// Set all parameters at their starting values
for (UInt_t i = 0; i < fParams.size(); i++){
fCurVals[fParams[i]] = fStartVals[fParams[i]];
}
TDirectory* nomfolder = (TDirectory*) limfolder->mkdir("nominal");
nomfolder->cd();
UpdateRWEngine(fCurVals);
fSampleFCN->ReconfigureAllEvents();
fSampleFCN->Write();
limfolder->cd();
std::vector<std::string> allfolders;
// Loop through each parameter
for (UInt_t i = 0; i < fParams.size(); i++){
std::string syst = fParams[i];
std::cout << "Starting Param " << syst << std::endl;
if (fFixVals[syst]) continue;
// Loop Downwards
while (fCurVals[syst] > fMinVals[syst]){
fCurVals[syst] = fCurVals[syst] - fStepVals[syst];
// Check Limit
if (fCurVals[syst] < fMinVals[syst])
fCurVals[syst] = fMinVals[syst];
// Check folder exists
std::string curvalstring = std::string( Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
if (std::find(allfolders.begin(), allfolders.end(), curvalstring) != allfolders.end())
break;
// Make new folder for variation
TDirectory* minfolder = (TDirectory*) limfolder->mkdir(Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
minfolder->cd();
allfolders.push_back(curvalstring);
// Update Iterations
double *vals = FitUtils::GetArrayFromMap( fParams, fCurVals );
fSampleFCN->DoEval( vals );
delete vals;
// Save to folder
fSampleFCN->Write();
}
// Reset before next loop
fCurVals[syst] = fStartVals[syst];
// Loop Upwards now
while (fCurVals[syst] < fMaxVals[syst]){
fCurVals[syst] = fCurVals[syst] + fStepVals[syst];
// Check Limit
if (fCurVals[syst] > fMaxVals[syst])
fCurVals[syst] = fMaxVals[syst];
// Check folder exists
std::string curvalstring = std::string( Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
if (std::find(allfolders.begin(), allfolders.end(), curvalstring) != allfolders.end())
break;
// Make new folder
TDirectory* maxfolder = (TDirectory*) limfolder->mkdir(Form( (syst + "_%f").c_str(), fCurVals[syst] ) );
maxfolder->cd();
allfolders.push_back(curvalstring);
// Update Iterations
double *vals = FitUtils::GetArrayFromMap( fParams, fCurVals );
fSampleFCN->DoEval( vals );
delete vals;
// Save to file
fSampleFCN->Write();
}
// Reset before leaving
fCurVals[syst] = fStartVals[syst];
UpdateRWEngine(fCurVals);
}
return;
}
//*************************************
void SystematicRoutines::Run(){
//*************************************
std::cout << "Running routines "<< std::endl;
fRoutines = GeneralUtils::ParseToStr(fStrategy,",");
for (UInt_t i = 0; i < fRoutines.size(); i++){
std::string routine = fRoutines.at(i);
int fitstate = kFitUnfinished;
LOG(FIT)<<"Running Routine: "<<routine<<std::endl;
if (routine.compare("PlotLimits") == 0) PlotLimits();
else if (routine.compare("ErrorBands") == 0) GenerateErrorBands();
else if (routine.compare("ThrowErrors") == 0) GenerateThrows();
else if (routine.compare("MergeErrors") == 0) MergeThrows();
else {
std::cout << "Unknown ROUTINE : " << routine << std::endl;
}
// If ending early break here
if (fitstate == kFitFinished || fitstate == kNoChange){
LOG(FIT) << "Ending fit routines loop." << std::endl;
break;
}
}
return;
}
void SystematicRoutines::GenerateErrorBands(){
GenerateThrows();
MergeThrows();
}
//*************************************
void SystematicRoutines::GenerateThrows(){
//*************************************
TFile* tempfile = new TFile((fOutputFile + ".throws.root").c_str(),"RECREATE");
tempfile->cd();
int nthrows = fNThrows;
int startthrows = fStartThrows;
int endthrows = startthrows + nthrows;
if (nthrows < 0) nthrows = endthrows;
if (startthrows < 0) startthrows = 0;
if (endthrows < 0) endthrows = startthrows + nthrows;
int seed = (gRandom->Uniform(0.0,1.0)*100000 + 100000000*(startthrows + endthrows) + time(NULL))/35;
gRandom->SetSeed(seed);
LOG(FIT) << "Using Seed : " << seed << std::endl;
LOG(FIT) << "nthrows = " << nthrows << std::endl;
LOG(FIT) << "startthrows = " << startthrows << std::endl;
LOG(FIT) << "endthrows = " << endthrows << std::endl;
UpdateRWEngine(fCurVals);
fSampleFCN->ReconfigureAllEvents();
if (startthrows == 0){
LOG(FIT) << "Making nominal " << std::endl;
TDirectory* nominal = (TDirectory*) tempfile->mkdir("nominal");
nominal->cd();
fSampleFCN->Write();
}
LOG(SAM) << "nthrows = " << nthrows << std::endl;
LOG(SAM) << "startthrows = " << startthrows << std::endl;
LOG(SAM) << "endthrows = " << endthrows << std::endl;
TTree* parameterTree = new TTree("throws","throws");
double chi2;
for (UInt_t i = 0; i < fParams.size(); i++)
parameterTree->Branch(fParams[i].c_str(), &fThrownVals[fParams[i]], (fParams[i] + "/D").c_str());
parameterTree->Branch("chi2",&chi2,"chi2/D");
fSampleFCN->CreateIterationTree("error_iterations", FitBase::GetRW());
// Would anybody actually want to do uniform throws of any parameter??
bool uniformly = FitPar::Config().GetParB("error_uniform");
// Run Throws and save
for (Int_t i = 0; i < endthrows+1; i++){
LOG(FIT) << "Loop " << i << std::endl;
ThrowCovariance(uniformly);
if (i < startthrows) continue;
if (i == 0) continue;
LOG(FIT) << "Throw " << i << " ================================" << std::endl;
// Generate Random Parameter Throw
// ThrowCovariance(uniformly);
TDirectory* throwfolder = (TDirectory*)tempfile->mkdir(Form("throw_%i",i));
throwfolder->cd();
// Run Eval
double *vals = FitUtils::GetArrayFromMap( fParams, fThrownVals );
chi2 = fSampleFCN->DoEval( vals );
delete vals;
// Save the FCN
fSampleFCN->Write();
parameterTree->Fill();
}
fSampleFCN->WriteIterationTree();
tempfile->Close();
}
void SystematicRoutines::MergeThrows(){
fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE");
fOutputRootFile->cd();
// Make a container folder
TDirectory* errorDIR = (TDirectory*) fOutputRootFile->mkdir("error_bands");
errorDIR->cd();
TDirectory* outnominal = (TDirectory*) fOutputRootFile->mkdir("nominal_throw");
outnominal->cd();
// Split Input Files
if (!fThrowString.empty()) fThrowList = GeneralUtils::ParseToStr(fThrowString,",");
// Add default if no throwlist given
if (fThrowList.size() < 1) fThrowList.push_back( fOutputFile + ".throws.root" );
/// Save location of file containing nominal
std::string nominalfile;
bool nominalfound;
// Loop over files and check they exist.
for (uint i = 0; i < fThrowList.size(); i++){
std::string file = fThrowList[i];
bool found = false;
// normal
std::string newfile = file;
TFile* throwfile = new TFile(file.c_str(),"READ");
if (throwfile and !throwfile->IsZombie()){
found = true;
}
// normal.throws.root
if (!found){
newfile = file + ".throws.root";
throwfile = new TFile((file + ".throws.root").c_str(),"READ");
if (throwfile and !throwfile->IsZombie()) {
found = true;
}
}
// If its found save to throwlist, else save empty.
// Also search for nominal
if (found){
fThrowList[i] = newfile;
LOG(FIT) << "Throws File :" << newfile << std::endl;
// Find input which contains nominal
if (throwfile->Get("nominal")){
nominalfound = true;
nominalfile = newfile;
}
throwfile->Close();
} else {
fThrowList[i] = "";
}
delete throwfile;
}
// Make sure we have a nominal file
if (!nominalfound or nominalfile.empty()){
ERR(FTL) << "No nominal found when mergining! Exiting!" << std::endl;
throw;
}
// Get the nominal throws file
TFile* tempfile = new TFile((nominalfile).c_str(),"READ");
tempfile->cd();
TDirectory* nominal = (TDirectory*)tempfile->Get("nominal");
// int nthrows = FitPar::Config().GetParI("error_throws");
bool uniformly = FitPar::Config().GetParB("error_uniform");
// Check percentage of bad files is okay.
int badfilecount = 0;
for (uint i = 0; i < fThrowList.size(); i++){
if (!fThrowList[i].empty()){
LOG(FIT) << "Loading Throws From File " << i << " : "
<< fThrowList[i] << std::endl;
} else {
badfilecount++;
}
}
// Check we have at least one good file
if ((uint)badfilecount == fThrowList.size()){
ERR(FTL) << "Found no good throw files for MergeThrows" << std::endl;
throw;
} else if (badfilecount > fThrowList.size()*0.25){
ERR(WRN) << "Over 25% of your throw files are dodgy. Please check this is okay!" << std::endl;
ERR(WRN) << "Will continue for the time being..." << std::endl;
sleep(5);
}
// Now go through the keys in the temporary file and look for TH1D, and TH2D plots
TIter next(nominal->GetListOfKeys());
TKey *key;
while ((key = (TKey*)next())) {
TClass *cl = gROOT->GetClass(key->GetClassName());
if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue;
TH1* baseplot = (TH1D*)key->ReadObj();
std::string plotname = std::string(baseplot->GetName());
LOG(FIT) << "Creating error bands for " << plotname;
if (LOG_LEVEL(FIT)){
if (!uniformly) std::cout << " : Using COVARIANCE Throws! " << std::endl;
else std::cout << " : Using UNIFORM THROWS!!! " << std::endl;
}
int nbins = 0;
if (cl->InheritsFrom("TH1D")) nbins = ((TH1D*)baseplot)->GetNbinsX();
else nbins = ((TH1D*)baseplot)->GetNbinsX()* ((TH1D*)baseplot)->GetNbinsY();
// Setup TProfile with RMS option
TProfile* tprof = new TProfile((plotname + "_prof").c_str(),(plotname + "_prof").c_str(),nbins, 0, nbins, "S");
// Setup The TTREE
double* bincontents;
bincontents = new double[nbins];
double* binlowest;
binlowest = new double[nbins];
double* binhighest;
binhighest = new double[nbins];
errorDIR->cd();
TTree* bintree = new TTree((plotname + "_tree").c_str(), (plotname + "_tree").c_str());
for (Int_t i = 0; i < nbins; i++){
bincontents[i] = 0.0;
binhighest[i] = 0.0;
binlowest[i] = 0.0;
bintree->Branch(Form("content_%i",i),&bincontents[i],Form("content_%i/D",i));
}
// Make new throw plot
TH1* newplot;
// Run Throw Merging.
for (UInt_t i = 0; i < fThrowList.size(); i++){
TFile* throwfile = new TFile(fThrowList[i].c_str(), "READ");
// Loop over all throws in a folder
TIter nextthrow(throwfile->GetListOfKeys());
TKey *throwkey;
while ((throwkey = (TKey*)nextthrow())) {
// Skip non throw folders
if (std::string(throwkey->GetName()).find("throw_") == std::string::npos) continue;
// Get Throw DIR
TDirectory* throwdir = (TDirectory*)throwkey->ReadObj();
// Get Plot From Throw
newplot = (TH1*)throwdir->Get(plotname.c_str());
if (!newplot) continue;
// Loop Over Plot
for (Int_t j = 0; j < nbins; j++){
tprof->Fill(j+0.5, newplot->GetBinContent(j+1));
bincontents[j] = newplot->GetBinContent(j+1);
if (bincontents[j] < binlowest[j] or i == 0) binlowest[j] = bincontents[j];
if (bincontents[j] > binhighest[j] or i == 0) binhighest[j] = bincontents[j];
}
errorDIR->cd();
bintree->Fill();
}
}
errorDIR->cd();
if (uniformly){
LOG(FIT) << "Uniformly Calculating Plot Errors!" << std::endl;
}
TH1* statplot = (TH1*) baseplot->Clone();
for (Int_t j = 0; j < nbins; j++){
if (!uniformly){
// if ((baseplot->GetBinError(j+1)/baseplot->GetBinContent(j+1)) < 1.0) {
// baseplot->SetBinError(j+1,sqrt(pow(tprof->GetBinError(j+1),2) + pow(baseplot->GetBinError(j+1),2)));
// } else {
baseplot->SetBinError(j+1,tprof->GetBinError(j+1));
// }
} else {
baseplot->SetBinContent(j+1, 0.0);//(binlowest[j] + binhighest[j]) / 2.0);
baseplot->SetBinError(j+1, 0.0); //(binhighest[j] - binlowest[j])/2.0);
}
}
errorDIR->cd();
baseplot->Write();
tprof->Write();
bintree->Write();
outnominal->cd();
for (int i = 0; i < nbins; i++){
baseplot->SetBinError(i+1, sqrt(pow(statplot->GetBinError(i+1),2) + pow(baseplot->GetBinError(i+1),2)));
}
baseplot->Write();
delete statplot;
delete baseplot;
delete tprof;
delete bintree;
delete [] bincontents;
}
return;
};
diff --git a/src/Utils/SignalDef.cxx b/src/Utils/SignalDef.cxx
index 5fb0923..f85b4f0 100644
--- a/src/Utils/SignalDef.cxx
+++ b/src/Utils/SignalDef.cxx
@@ -1,289 +1,289 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#include "FitUtils.h"
#include "SignalDef.h"
bool SignalDef::isCCINC(FitEvent *event, int nuPDG, double EnuMin, double EnuMax) {
// Check for the desired PDG code
if (!event->HasISParticle(nuPDG)) return false;
// Check that it's within the allowed range if set
if (EnuMin != EnuMax) {
if (!SignalDef::IsEnuInRange(event, EnuMin*1000, EnuMax*1000)) {
return false;
}
}
// Check that the charged lepton we expect has been produced
if (!event->HasFSParticle(nuPDG > 0 ? nuPDG-1 : nuPDG+1)) return false;
return true;
}
bool SignalDef::isNCINC(FitEvent *event, int nuPDG, double EnuMin, double EnuMax) {
// Check for the desired PDG code before and after the interaction
if (!event->HasISParticle(nuPDG) ||
!event->HasFSParticle(nuPDG)) return false;
// Check that it's within the allowed range if set
if (EnuMin != EnuMax)
if (!SignalDef::IsEnuInRange(event, EnuMin*1000, EnuMax*1000))
return false;
return true;
}
bool SignalDef::isCC0pi(FitEvent *event, int nuPDG, double EnuMin, double EnuMax){
// Check it's CCINC
if (!SignalDef::isCCINC(event, nuPDG, EnuMin, EnuMax)) return false;
// Veto event with mesons
if (event->NumFSMesons() != 0) return false;
// Veto events which don't have exactly 1 outgoing charged lepton
if (event->NumFSLeptons() != 1) return false;
return true;
}
bool SignalDef::isNC0pi(FitEvent *event, int nuPDG, double EnuMin, double EnuMax){
// Check it's NCINC
if (!SignalDef::isNCINC(event, nuPDG, EnuMin, EnuMax)) return false;
// Veto event with mesons
if (event->NumFSMesons() != 0) return false;
// Veto events with a charged lepton
if (event->NumFSLeptons() != 0) return false;
return true;
}
bool SignalDef::isCCQE(FitEvent *event, int nuPDG, double EnuMin, double EnuMax){
// Check if it's CCINC
if (!SignalDef::isCCINC(event, nuPDG, EnuMin, EnuMax)) return false;
// Require modes 1 (CCQE)
if (abs(event->Mode) != 1) return false;
return true;
}
bool SignalDef::isNCEL(FitEvent *event, int nuPDG, double EnuMin, double EnuMax){
// Check if it's NCINC
if (!SignalDef::isNCINC(event, nuPDG, EnuMin, EnuMax)) return false;
// Require modes 51/52 (NCEL)
if (abs(event->Mode) != 51 && abs(event->Mode) != 52) return false;
return true;
}
bool SignalDef::isCCQELike(FitEvent *event, int nuPDG, double EnuMin, double EnuMax){
// Check if it's CCINC
if (!SignalDef::isCCINC(event, nuPDG, EnuMin, EnuMax)) return false;
// Require modes 1/2 (CCQE and MEC)
if (abs(event->Mode) != 1 && abs(event->Mode) != 2) return false;
return true;
}
// Require one meson, one charged lepton. types specified in the arguments
bool SignalDef::isCC1pi(FitEvent *event, int nuPDG, int piPDG,
double EnuMin, double EnuMax){
// First, make sure it's CCINC
if (!SignalDef::isCCINC(event, nuPDG, EnuMin, EnuMax)) return false;
int nMesons = event->NumFSMesons();
int nLeptons = event->NumFSLeptons();
int nPion = event->NumFSParticle(piPDG);
// Check that the desired pion exists and is the only meson
if (nPion != 1 || nMesons != 1) return false;
// Check that there is only one final state lepton
if (nLeptons != 1) return false;
return true;
}
// Require one meson, one neutrino. Types specified as arguments
bool SignalDef::isNC1pi(FitEvent *event, int nuPDG, int piPDG,
double EnuMin, double EnuMax){
// First, make sure it's NCINC
if (!SignalDef::isNCINC(event, nuPDG, EnuMin, EnuMax)) return false;
int nMesons = event->NumFSMesons();
int nLeptons = event->NumFSLeptons();
int nPion = event->NumFSParticle(piPDG);
// Check that the desired pion exists and is the only meson
if (nPion != 1 || nMesons != 1) return false;
// Check that there are no charged leptons
if (nLeptons != 0) return false;
return true;
}
// A slightly ugly function to replace the BC 2pi channels.
// All particles which are allowed in the final state are specified
bool SignalDef::isCCWithFS(FitEvent *event, int nuPDG, std::vector<int> pdgs,
double EnuMin, double EnuMax){
// Check it's CCINC
if (!SignalDef::isCCINC(event, nuPDG, EnuMin, EnuMax)) return false;
// Remove events where the number of final state particles
// do not match the number specified in the signal definition
if ((int)pdgs.size() != event->NumFSParticle()) return false;
// For every particle in the list, check the number in the FS
for (std::vector<int>::iterator it = pdgs.begin(); it != pdgs.end(); ++it){
// Check how many times this pdg is in the vector
int nEntries = std::count (pdgs.begin(), pdgs.end(), *it);
if (event->NumFSParticle(*it) != nEntries)
return false;
}
return true;
}
// Require one meson, one charged lepton, AND specify the only other final state particle
// This is only suitable for bubble chambers. Types specified in the arguments
bool SignalDef::isCC1pi3Prong(FitEvent *event, int nuPDG, int piPDG,
int thirdPDG, double EnuMin, double EnuMax){
// First, make sure it's CC1pi
if (!SignalDef::isCC1pi(event, nuPDG, piPDG, EnuMin, EnuMax)) return false;
// Check we have the third prong
if (event->NumFSParticle(thirdPDG) == 0) return false;
// Check that there are only three FS particles
if (event->NumFSParticle() != 3) return false;
return true;
}
// Require one meson, one neutrino, AND specify the only other final state particle
// This is only suitable for bubble chambers. Types specified in the arguments
bool SignalDef::isNC1pi3Prong(FitEvent *event, int nuPDG, int piPDG,
int thirdPDG, double EnuMin, double EnuMax){
// First, make sure it's NC1pi
if (!SignalDef::isNC1pi(event, nuPDG, piPDG, EnuMin, EnuMax)) return false;
// Check we have the third prong
if (event->NumFSParticle(thirdPDG) == 0) return false;
// Check that there are only three FS particles
if (event->NumFSParticle() != 3) return false;
return true;
}
bool SignalDef::isCCCOH(FitEvent *event, int nuPDG, int piPDG, double EnuMin, double EnuMax){
// Check this is a CCINC event
if (!SignalDef::isCCINC(event, nuPDG, EnuMin, EnuMax)) return false;
int nLepton = event->NumFSParticle(nuPDG > 0 ? nuPDG-1 : nuPDG+1);
int nPion = event->NumFSParticle(piPDG);
int nFS = event->NumFSParticle();
if (nLepton != 1 || nPion != 1) return false;
if (nFS != 2) return false;
return true;
}
bool SignalDef::isNCCOH(FitEvent *event, int nuPDG, int piPDG, double EnuMin, double EnuMax){
// Check this is an NCINC event
if (!SignalDef::isNCINC(event, nuPDG, EnuMin, EnuMax)) return false;
int nLepton = event->NumFSParticle(nuPDG);
int nPion = event->NumFSParticle(piPDG);
int nFS = event->NumFSParticle();
if (nLepton != 1 || nPion != 1) return false;
if (nFS != 2) return false;
return true;
}
bool SignalDef::HasProtonKEAboveThreshold(FitEvent* event, double threshold){
for (uint i = 0; i < event->Npart(); i++){
FitParticle* p = event->PartInfo(i);
if (!p->IsFinalState()) continue;
if (p->fPID != 2212) continue;
if (FitUtils::T(p->fP) > threshold / 1000.0) return true;
}
return false;
}
bool SignalDef::HasProtonMomAboveThreshold(FitEvent* event, double threshold){
for (uint i = 0; i < event->Npart(); i++){
FitParticle* p = event->PartInfo(i);
if (!p->IsFinalState()) continue;
if (p->fPID != 2212) continue;
if (p->fP.Vect().Mag() > threshold) return true;
}
return false;
}
// Calculate the angle between the neutrino and an outgoing particle, apply a cut
bool SignalDef::IsRestrictedAngle(FitEvent* event, int nuPDG, int otherPDG, double angle){
// If the particles don't exist, return false
if (!event->HasISParticle(nuPDG) || !event->HasFSParticle(otherPDG)) return false;
// Get Mom
TVector3 pnu = event->GetHMISParticle(nuPDG)->fP.Vect();
TVector3 p2 = event->GetHMFSParticle(otherPDG)->fP.Vect();
double theta = pnu.Angle(p2) * 180. / TMath::Pi();
return (theta < angle);
}
bool SignalDef::IsEnuInRange(FitEvent* event, double emin, double emax){
- return (event->PartInfo(0)->fP.E() > emin &&
- event->PartInfo(0)->fP.E() < emax);
+ return (event->GetNeutrinoIn()->fP.E() > emin &&
+ event->GetNeutrinoIn()->fP.E() < emax);
}
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Sat, May 3, 7:04 AM (6 h, 15 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4983270
Default Alt Text
(248 KB)
Attached To
rNUISANCEGIT nuisancegit
Event Timeline
Log In to Comment