diff --git a/CMakeLists.txt b/CMakeLists.txt
index 068cf84..6f53a19 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,196 +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 .
################################################################################
cmake_minimum_required (VERSION 2.6 FATAL_ERROR)
project(NUISANCE)
include(ExternalProject)
enable_language(Fortran)
set (NUISANCE_VERSION_MAJOR 1)
set (NUISANCE_VERSION_MINOR 0)
set (NUISANCE_VERSION_REVISION 0)
set (NUISANCE_VERSION_STRING "v${NUISANCE_VERSION_MAJOR}r${NUISANCE_VERSION_MINOR}")
if(${NUISANCE_VERSION_REVISION} STRGREATER "0")
set (NUISANCE_VERSION_STRING "${NUISANCE_VERSION_STRING}p${NUISANCE_VERSION_REVISION}")
endif()
#Set this to TRUE to enable build debugging messages
set(BUILD_DEBUG_MSGS TRUE)
include(${CMAKE_SOURCE_DIR}/cmake/cmessage.cmake)
include(${CMAKE_SOURCE_DIR}/cmake/cacheVariables.cmake)
cmessage(STATUS "CMAKE_INSTALL_PREFIX: \"${CMAKE_INSTALL_PREFIX}\"")
cmessage(STATUS "CMAKE_BUILD_TYPE: \"${CMAKE_BUILD_TYPE}\"")
################################################################################
# Check Dependencies
################################################################################
################################## ROOT ######################################
include(${CMAKE_SOURCE_DIR}/cmake/ROOTSetup.cmake)
################################# HEPMC ######################################
include(${CMAKE_SOURCE_DIR}/cmake/HepMC.cmake)
############################ Reweight Engines ################################
include(${CMAKE_SOURCE_DIR}/cmake/ReweightEnginesSetup.cmake)
############################ Other Generators ################################
include(${CMAKE_SOURCE_DIR}/cmake/GiBUUSetup.cmake)
if(USE_NUANCE)
LIST(APPEND EXTRA_CXX_FLAGS -D__NUANCE_ENABLED__)
endif()
################################# Pythia6/8 ####################################
include(${CMAKE_SOURCE_DIR}/cmake/pythia6Setup.cmake)
include(${CMAKE_SOURCE_DIR}/cmake/pythia8Setup.cmake)
################################# gperftools ###################################
include(${CMAKE_SOURCE_DIR}/cmake/gperfSetup.cmake)
if(NOT NOTEST)
enable_testing()
endif()
SET(GENERATOR_SUPPORT)
foreach(gen NEUT;NuWro;GENIE;GiBUU;NUANCE)
if(USE_${gen})
SET(GENERATOR_SUPPORT "${GENERATOR_SUPPORT}${gen} ")
endif()
endforeach(gen)
cmessage(STATUS "Generator Input Support: ${GENERATOR_SUPPORT}")
set(MINCODE
Routines
FCN)
set(CORE
MCStudies
Genie
FitBase
Config
Logger
- Statistical
InputHandler
Splines
Reweight
Utils
+ Statistical
#Devel
Smearceptance
)
LIST(APPEND ALLEXPERIMENTS
ANL
ArgoNeuT
BEBC
BNL
Electron
FNAL
GGM
K2K
MINERvA
MiniBooNE
SciBooNE
T2K)
foreach(exp ${ALLEXPERIMENTS})
if(NOT NO_${exp})
LIST(APPEND EXPERIMENTS_TO_BUILD ${exp})
else()
LIST(REVERSE EXTRA_CXX_FLAGS)
LIST(APPEND EXTRA_CXX_FLAGS -D__NO_${exp}__)
LIST(REVERSE EXTRA_CXX_FLAGS)
endif()
endforeach()
################################## COMPILER ####################################
include(${CMAKE_SOURCE_DIR}/cmake/c++CompilerSetup.cmake)
################################### doxygen ###################################
include(${CMAKE_SOURCE_DIR}/cmake/docsSetup.cmake)
################################################################################
set(MINIMUM_INCLUDE_DIRECTORIES)
LIST(APPEND MINIMUM_INCLUDE_DIRECTORIES
${RWENGINE_INCLUDE_DIRECTORIES}
${CMAKE_SOURCE_DIR}/src/FitBase
${CMAKE_SOURCE_DIR}/src/Reweight
${CMAKE_SOURCE_DIR}/src/InputHandler
${CMAKE_SOURCE_DIR}/src/Config
${CMAKE_SOURCE_DIR}/src/Logger
${CMAKE_SOURCE_DIR}/src/Statistical
${CMAKE_SOURCE_DIR}/src/Splines
${CMAKE_SOURCE_DIR}/src/Utils
${CMAKE_SOURCE_DIR}/src/Genie)
cmessage(DEBUG "Base include directories: ${MINIMUM_INCLUDE_DIRECTORIES}")
set(EXP_INCLUDE_DIRECTORIES)
foreach(edir ${EXPERIMENTS_TO_BUILD})
LIST(APPEND EXP_INCLUDE_DIRECTORIES ${CMAKE_SOURCE_DIR}/src/${edir})
endforeach()
cmessage(DEBUG "Included experiments: ${EXP_INCLUDE_DIRECTORIES}")
foreach(mdir ${MINCODE})
cmessage (DEBUG "Configuring directory: src/${mdir}")
add_subdirectory(src/${mdir})
endforeach()
foreach(edir ${EXPERIMENTS_TO_BUILD})
cmessage (DEBUG "Configuring directory: src/${edir}")
add_subdirectory(src/${edir})
endforeach()
foreach(cdir ${CORE})
cmessage (DEBUG "Configuring directory: src/${cdir}")
add_subdirectory(src/${cdir})
endforeach()
cmessage(DEBUG "Module targets: ${MODULETargets}")
add_subdirectory(app)
add_subdirectory(src/Tests)
configure_file(cmake/setup.sh.in
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.sh" @ONLY)
install(FILES
"${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.sh" DESTINATION
${CMAKE_INSTALL_PREFIX})
install(PROGRAMS
"${PROJECT_SOURCE_DIR}/scripts/nuiscardgen" DESTINATION
bin)
install(PROGRAMS
"${PROJECT_SOURCE_DIR}/scripts/nuissamples" DESTINATION
bin)
diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt
index e4e2350..1d488e3 100644
--- a/app/CMakeLists.txt
+++ b/app/CMakeLists.txt
@@ -1,217 +1,217 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
set(TARGETS_TO_BUILD)
if(USE_MINIMIZER)
add_executable(nuismin nuismin.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuismin)
target_link_libraries(nuismin ${MODULETargets})
target_link_libraries(nuismin ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuismin ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuismin PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuissplines nuissplines.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissplines)
target_link_libraries(nuissplines ${MODULETargets})
target_link_libraries(nuissplines ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuissplines ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuissplines PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
include_directories(${RWENGINE_INCLUDE_DIRECTORIES})
include_directories(${CMAKE_SOURCE_DIR}/src/Routines)
include_directories(${CMAKE_SOURCE_DIR}/src/InputHandler)
include_directories(${CMAKE_SOURCE_DIR}/src/Genie)
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
+include_directories(${CMAKE_SOURCE_DIR}/src/Statistical)
include_directories(${CMAKE_SOURCE_DIR}/src/Utils)
include_directories(${CMAKE_SOURCE_DIR}/src/Config)
include_directories(${CMAKE_SOURCE_DIR}/src/Logger)
-include_directories(${CMAKE_SOURCE_DIR}/src/Statistical)
include_directories(${CMAKE_SOURCE_DIR}/src/Splines)
include_directories(${CMAKE_SOURCE_DIR}/src/Reweight)
include_directories(${CMAKE_SOURCE_DIR}/src/FCN)
include_directories(${CMAKE_SOURCE_DIR}/src/MCStudies)
include_directories(${CMAKE_SOURCE_DIR}/src/Smearceptance)
include_directories(${EXP_INCLUDE_DIRECTORIES})
if (USE_NuWro AND NOT NUWRO_BUILT_FROM_FILE)
add_executable(nuwro_nuisance nuwro_NUISANCE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuwro_nuisance)
target_link_libraries(nuwro_nuisance ${MODULETargets})
target_link_libraries(nuwro_nuisance ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuwro_nuisance ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuwro_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
if (USE_NEUT)
add_executable(neut_nuisance neut_NUISANCE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};neut_nuisance)
target_link_libraries(neut_nuisance ${MODULETargets})
target_link_libraries(neut_nuisance ${CMAKE_DEPENDLIB_FLAGS})
target_link_libraries(neut_nuisance ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(neut_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
if (BUILD_GEVGEN)
add_executable(gevgen_nuisance gEvGen_NUISANCE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};gevgen_nuisance)
target_link_libraries(gevgen_nuisance ${MODULETargets})
target_link_libraries(gevgen_nuisance ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(gevgen_nuisance ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
include_directories(${GENIE_INCLUDES}/Apps)
include_directories(${GENIE_INCLUDES}/FluxDrivers)
include_directories(${GENIE_INCLUDES}/EVGDrivers)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(gevgen_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(gevgen_nuisance_mixed gEvGen_NUISANCE_MIXED.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};gevgen_nuisance_mixed)
target_link_libraries(gevgen_nuisance_mixed ${MODULETargets})
target_link_libraries(gevgen_nuisance_mixed ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(gevgen_nuisance_mixed ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
include_directories(${GENIE_INCLUDES}/Apps)
include_directories(${GENIE_INCLUDES}/FluxDrivers)
include_directories(${GENIE_INCLUDES}/EVGDrivers)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(gevgen_nuisance_mixed PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
if (USE_GiBUU)
add_executable(DumpGiBUUEvents DumpGiBUUEvents.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};DumpGiBUUEvents)
target_link_libraries(DumpGiBUUEvents ${MODULETargets})
target_link_libraries(DumpGiBUUEvents ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(DumpGiBUUEvents ${ROOT_LIBS})
include_directories(${CMAKE_SOURCE_DIR}/src/FitBase)
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(DumpGiBUUEvents PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
add_executable(nuiscomp nuiscomp.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuiscomp)
target_link_libraries(nuiscomp ${MODULETargets})
target_link_libraries(nuiscomp ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuiscomp ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuiscomp PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuisflat nuisflat.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuisflat)
target_link_libraries(nuisflat ${MODULETargets})
target_link_libraries(nuisflat ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuisflat ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuisflat PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuissmear nuissmear.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissmear)
target_link_libraries(nuissmear ${MODULETargets})
target_link_libraries(nuissmear ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuissmear ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuissmear PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuissyst nuissyst.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissyst)
target_link_libraries(nuissyst ${MODULETargets})
target_link_libraries(nuissyst ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuissyst ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuissyst PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
add_executable(nuisbayes nuisbayes.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuisbayes)
target_link_libraries(nuisbayes ${MODULETargets})
target_link_libraries(nuisbayes ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(nuisbayes ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(nuisbayes PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
if(USE_GENIE)
add_executable(PrepareGENIE PrepareGENIE.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareGENIE)
target_link_libraries(PrepareGENIE ${MODULETargets})
target_link_libraries(PrepareGENIE ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(PrepareGENIE ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(PrepareGENIE PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
if(USE_NEUT)
add_executable(PrepareNEUT PrepareNEUT.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareNEUT)
target_link_libraries(PrepareNEUT ${MODULETargets})
target_link_libraries(PrepareNEUT ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(PrepareNEUT ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(PrepareNEUT PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
# PREPARE NUWRO
# Commented out for the time being until it is finished..
if(USE_NuWro)
add_executable(PrepareNuwro PrepareNuwroEvents.cxx)
set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareNuwro)
target_link_libraries(PrepareNuwro ${MODULETargets})
target_link_libraries(PrepareNuwro ${CMAKE_DEPENDLIB_FLAGS})
# target_link_libraries(PrepareNuwro ${ROOT_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(PrepareNuwro PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
endif()
install(TARGETS ${TARGETS_TO_BUILD} DESTINATION bin)
#add_executable(DumpROOTClassesFromVector DumpROOTClassesFromVector.cxx)
# #Strip out -lNuWro_event1
# string(REPLACE "-lNuWro_event1" "" NWEVSTRIPPED_CDF ${CMAKE_DEPENDLIB_FLAGS})
# cmessage(DEBUG "Attempted to strip out nuwro library: \"${CMAKE_DEPENDLIB_FLAGS}\" -> \"${NWEVSTRIPPED_CDF}\"")
# add_executable(PrepareNEUT PrepareNEUT.cxx)
# target_link_libraries(DumpROOTClassesFromVector ${MODULETargets})
# target_link_libraries(DumpROOTClassesFromVector ${NWEVSTRIPPED_CDF})
# if(NOT CMAKE_LINK_FLAGS STREQUAL "")
# set_target_properties(DumpROOTClassesFromVector PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
# endif()
#install(TARGETS DumpROOTClassesFromVector DESTINATION bin)
diff --git a/app/PrepareNEUT.cxx b/app/PrepareNEUT.cxx
index a4c31ed..3b37e40 100644
--- a/app/PrepareNEUT.cxx
+++ b/app/PrepareNEUT.cxx
@@ -1,237 +1,238 @@
#include
#include
#include "TFile.h"
#include "TH1D.h"
#include "TTree.h"
#include "PlotUtils.h"
+#include "StatUtils.h"
#include "FitLogger.h"
// If you don't have NEUT enabled, you shouldn't compile this...
#include "neutpart.h"
#include "neutvect.h"
std::string fInputFiles = "";
std::string fOutputFile = "";
std::string fFluxFile = "";
bool fFluxInGeV = false;
void PrintOptions();
void ParseOptions(int argc, char* argv[]);
void CreateRateHistogram(std::string inputList, std::string flux, std::string output);
//*******************************
int main(int argc, char* argv[]){
//*******************************
LOG_VERB(FitPar::Config().GetParI("VERBOSITY"));
ERR_VERB(FitPar::Config().GetParI("ERROR"));
ParseOptions(argc, argv);
LOG(FIT) << "Running PrepareNEUT" << std::endl;
CreateRateHistogram(fInputFiles, fFluxFile, fOutputFile);
};
//*******************************
void CreateRateHistogram(std::string inputList, std::string flux, std::string output){
//*******************************
// Need to allow for more than one file... will do soon
TChain* tn = new TChain("neuttree");
std::vector inputs = GeneralUtils::ParseToStr(inputList,",");
for (std::vector::iterator it = inputs.begin();
it != inputs.end(); ++it){
LOG(FIT) << "Adding " << *it << " to the output" << std::endl;
tn->AddFile((*it).c_str());
}
if (inputs.size() > 1 && output.empty()){
ERR(FTL) << "You must provide a new output file name if you want to have more than 1 input file!" << std::endl;
throw;
}
int nevts = tn->GetEntries();
if (!nevts){
ERR(FTL) << "Either the input file is not from NEUT, or it's empty..." <SetBranchAddress("vectorbranch", &fNeutVect);
// Get Flux Hist
std::vector fluxvect = GeneralUtils::ParseToStr(flux,",");
TH1D* fluxHist = NULL;
if (fluxvect.size() > 1){
TFile* fluxfile = new TFile(fluxvect[0].c_str(),"READ");
fluxHist = (TH1D*) fluxfile->Get(fluxvect[1].c_str());
fluxHist->SetDirectory(0);
} else {
ERR(FTL) << "NO FLUX SPECIFIED" << std::endl;
throw;
}
// Decide what type of flux was given
if (fFluxInGeV) LOG(FIT) << "Assuming flux histogram is in GeV" << std::endl;
else LOG(FIT) << "Assuming flux histogram is in MeV" << std::endl;
// Make Event Hist
TH1D* xsecHist = (TH1D*)fluxHist->Clone();
xsecHist->Reset();
// Make a total cross section hist for shits and giggles
TH1D* entryHist = (TH1D*) xsecHist->Clone();
for (int i = 0; i < nevts; ++i){
tn->GetEntry(i);
NeutPart* part = fNeutVect->PartInfo(0);
double E = part->fP.E();
double xsec = fNeutVect->Totcrs;
// Unit conversion
if (fFluxInGeV) E*=1E-3;
xsecHist ->Fill(E, xsec);
entryHist->Fill(E);
if (i % (nevts/20) == 0){
LOG(FIT) << "Processed " << i << "/" << nevts << " NEUT events." << std::endl;
}
}
LOG(FIT) << "Processed all events" << std::endl;
xsecHist ->Divide(entryHist);
// This will be the evtrt histogram
TH1D* evtHist = NULL;
// If the integral of xsecHist is 0 the input file used a really old version of NEUT without Totcrs
if (!xsecHist->Integral(0, -1)){
ERR(WRN) << "Old NEUT input file: events will not be correctly normalized" << std::endl;
evtHist = (TH1D*)entryHist->Clone();
if (evtHist->Integral() != 0)
evtHist ->Scale(fluxHist->Integral()/float(evtHist->Integral()));
} else {
evtHist = (TH1D*)xsecHist->Clone();
evtHist ->Multiply(fluxHist);
}
// Check whether the overflow is empty. If not, advise that either the wrong flux histogram or units were used...
// If the events were generated with a limited range of the flux histogram, this may be benign
if (evtHist->Integral(0, -1) != evtHist->Integral() || evtHist->Integral(0, -1) == 0){
ERR(WRN) << "The input file and flux histogram provided do not match... " << std::endl;
ERR(WRN) << "Are the units correct? Did you provide the correct flux file?" << std::endl;
ERR(WRN) << "Use output with caution..." << std::endl;
}
// Pick where the output should go
TFile *outFile = NULL;
if (!output.empty()){
LOG(FIT) << "Saving histograms in " << output << std::endl;
outFile = new TFile(output.c_str(), "RECREATE");
} else {
LOG(FIT) << "Saving histograms in " << inputs[0] << std::endl;
outFile = new TFile(inputs[0].c_str(), "UPDATE");
}
outFile->cd();
std::string xsec_name = "xsec_PrepareNeut";
std::string flux_name = "flux_PrepareNeut";
std::string rate_name = "evtrt_PrepareNeut";
if (output.empty()){
// Check whether we should overwrite existing histograms
std::string input_xsec = PlotUtils::GetObjectWithName(outFile, "xsec");
std::string input_flux = PlotUtils::GetObjectWithName(outFile, "flux");
std::string input_rate = PlotUtils::GetObjectWithName(outFile, "evtrt");
if (!input_xsec.empty()) {
LOG(FIT) << "Updating histogram: " << input_xsec << std::endl;
xsec_name = input_xsec;
}
if (!input_flux.empty()) {
LOG(FIT) << "Updating histogram: " << input_flux << std::endl;
flux_name = input_flux;
}
if (!input_rate.empty()) {
LOG(FIT) << "Updating histogram: " << input_rate << std::endl;
rate_name = input_rate;
}
} else {
LOG(FIT) << "Cloning neuttree into output file." << std::endl;
StopTalking();
TTree* newtree = (TTree*)tn->CloneTree(-1, "fast");
StartTalking();
newtree->Write();
}
xsecHist->Write(xsec_name.c_str(), TObject::kOverwrite);
fluxHist->Write(flux_name.c_str(), TObject::kOverwrite);
evtHist ->Write(rate_name.c_str(), TObject::kOverwrite);
outFile ->Close();
return;
}
void PrintOptions(){
std::cout << "PrepareNEUT NUISANCE app. " << std::endl
<< "Produces or recalculates evtrt and flux histograms necessary for NUISANCE normalization." << std::endl;
std::cout << "PrepareNEUT: " << std::endl;
std::cout << " [-h,-help,--h,--help]" << std::endl;
std::cout << " -i inputfile1.root,inputfile2.root,inputfile3.root,..." << std::endl;
std::cout << " Takes any number of files, but assumes all are produced with a single flux" << std::endl;
std::cout << " -f flux_root_file.root,flux_hist_name" << std::endl;
std::cout << " Path to root file containing the flux histogram used when generating the NEUT files" << std::endl;
std::cout << " [-o outputfile.root] " << std::endl;
std::cout << " If an output file is not given, the input file will be used" << std::endl;
std::cout << " If more than one input file is given, an output file must be given" << std::endl;
std::cout << " [-G]" << std::endl;
std::cout << " Flux is assumed to be in MeV. This switch indicates the input flux is in GeV" << std::endl << std::endl;
}
void ParseOptions(int argc, char* argv[]){
bool flagopt = false;
// If No Arguments print commands
for (int i = 1; i< argc; ++i){
if (!std::strcmp(argv[i], "-h")) { flagopt = true; break; }
else if (!std::strcmp(argv[i], "-G")) { fFluxInGeV = true; }
if (i+1 != argc){
// Cardfile
if (!std::strcmp(argv[i], "-h")) { flagopt = true; break; }
else if (!std::strcmp(argv[i], "-i")) { fInputFiles = argv[i+1]; ++i; }
else if (!std::strcmp(argv[i], "-o")) { fOutputFile = argv[i+1]; ++i; }
else if (!std::strcmp(argv[i], "-f")) { fFluxFile = argv[i+1]; ++i; }
else {
ERR(FTL) << "ERROR: unknown command line option given! - '"
<.
################################################################################
set(HEADERFILES
FitUtils.h
GeneralUtils.h
ParserUtils.h
PlotUtils.h
SignalDef.h
BeamUtils.h
TargetUtils.h
OpenMPWrapper.h
PhysConst.h
)
set(IMPLFILES
FitUtils.cxx
GeneralUtils.cxx
PlotUtils.cxx
SignalDef.cxx
BeamUtils.cxx
TargetUtils.cxx
ParserUtils.cxx
)
set(LIBNAME Utils)
if(CMAKE_BUILD_TYPE MATCHES DEBUG)
add_library(${LIBNAME} STATIC ${IMPLFILES})
else(CMAKE_BUILD_TYPE MATCHES RELEASE)
add_library(${LIBNAME} SHARED ${IMPLFILES})
endif()
include_directories(${MINIMUM_INCLUDE_DIRECTORIES})
include_directories(${CMAKE_SOURCE_DIR}/src/Logger)
+include_directories(${CMAKE_SOURCE_DIR}/src/Statistical)
set_target_properties(${LIBNAME} PROPERTIES VERSION
"${NUISANCE_VERSION_MAJOR}.${NUISANCE_VERSION_MINOR}.${NUISANCE_VERSION_REVISION}")
#set_target_properties(${LIBNAME} PROPERTIES LINK_FLAGS ${ROOT_LD_FLAGS})
if(DEFINED PROJECTWIDE_EXTRA_DEPENDENCIES)
add_dependencies(${LIBNAME} ${PROJECTWIDE_EXTRA_DEPENDENCIES})
endif()
install(TARGETS ${LIBNAME} DESTINATION lib)
#Can uncomment this to install the headers... but is it really neccessary?
#install(FILES ${HEADERFILES} DESTINATION include)
set(MODULETargets ${MODULETargets} ${LIBNAME} PARENT_SCOPE)
#add_executable(DumpROOTClassesFromVector DumpROOTClassesFromVector.cxx GeneralUtils.cxx FitLogger.cxx PythiaQuiet.f)
#
# #Strip out -lNuWro_event1
# string(REPLACE "-lNuWro_event1" "" NWEVSTRIPPED_CDF ${CMAKE_DEPENDLIB_FLAGS})
# cmessage(DEBUG "Attempted to strip out nuwro library: \"${CMAKE_DEPENDLIB_FLAGS}\" -> \"${NWEVSTRIPPED_CDF}\"")
# target_link_libraries(DumpROOTClassesFromVector ${NWEVSTRIPPED_CDF})
# if(NOT CMAKE_LINK_FLAGS STREQUAL "")
# set_target_properties(DumpROOTClassesFromVector PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
# endif()
#install(TARGETS DumpROOTClassesFromVector DESTINATION bin)
diff --git a/src/Utils/PlotUtils.cxx b/src/Utils/PlotUtils.cxx
index c7d17fd..0a86acc 100644
--- a/src/Utils/PlotUtils.cxx
+++ b/src/Utils/PlotUtils.cxx
@@ -1,927 +1,927 @@
// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "PlotUtils.h"
#include "FitEvent.h"
-
+#include "StatUtils.h"
// MOVE TO MB UTILS!
// This function is intended to be modified to enforce a consistent masking for all models.
TH2D* PlotUtils::SetMaskHist(std::string type, TH2D* data) {
TH2D *fMaskHist = (TH2D*)data->Clone("fMaskHist");
for (int xBin = 0; xBin < fMaskHist->GetNbinsX(); ++xBin) {
for (int yBin = 0; yBin < fMaskHist->GetNbinsY(); ++yBin) {
if (data->GetBinContent(xBin + 1, yBin + 1) == 0) fMaskHist->SetBinContent(xBin + 1, yBin + 1, 0);
else fMaskHist->SetBinContent(xBin + 1, yBin + 1, 0.5);
if (!type.compare("MB_numu_2D")) {
if (yBin == 19 && xBin < 8) fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0);
} else {
if (yBin == 19 && xBin < 11) fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0);
}
if (yBin == 18 && xBin < 3) fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0);
if (yBin == 17 && xBin < 2) fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0);
if (yBin == 16 && xBin < 1) fMaskHist->SetBinContent(xBin + 1, yBin + 1, 1.0);
}
}
return fMaskHist;
};
// MOVE TO GENERAL UTILS?
bool PlotUtils::CheckObjectWithName(TFile *inFile, std::string substring) {
TIter nextkey(inFile->GetListOfKeys());
TKey *key;
while ( (key = (TKey*)nextkey()) ) {
std::string test(key->GetName());
if (test.find(substring) != std::string::npos) return true;
}
return false;
};
// MOVE TO GENERAL UTILS?
std::string PlotUtils::GetObjectWithName(TFile *inFile, std::string substring) {
TIter nextkey(inFile->GetListOfKeys());
TKey *key;
std::string output="";
while ( (key = (TKey*)nextkey()) ) {
std::string test(key->GetName());
if (test.find(substring) != std::string::npos)
output = test;
}
return output;
};
void PlotUtils::CreateNeutModeArray(TH1* hist, TH1* neutarray[]) {
for (int i = 0; i < 60; i++) {
neutarray[i] = (TH1*)hist->Clone(Form("%s_NMODE_%i", hist->GetName(), i));
}
return;
};
void PlotUtils::DeleteNeutModeArray(TH1* neutarray[]) {
for (int i = 0; i < 60; i++) {
delete neutarray[i];
}
return;
};
void PlotUtils::FillNeutModeArray(TH1D* hist[], int mode, double xval, double weight) {
if (abs(mode) > 60) return;
hist[abs(mode)]->Fill(xval, weight);
return;
};
void PlotUtils::FillNeutModeArray(TH2D* hist[], int mode, double xval, double yval, double weight) {
if (abs(mode) > 60) return;
hist[abs(mode)]->Fill(xval, yval, weight);
return;
};
THStack PlotUtils::GetNeutModeStack(std::string title, TH1* ModeStack[], int option) {
(void) option;
THStack allmodes = THStack(title.c_str(), title.c_str());
for (int i = 0; i < 60; i++) {
allmodes.Add(ModeStack[i]);
}
// Credit to Clarence for copying all this out.
// CC
ModeStack[1]->SetTitle("CCQE");
ModeStack[1]->SetFillColor(kBlue);
// ModeStack[1]->SetFillStyle(3444);
ModeStack[1]->SetLineColor(kBlue);
ModeStack[2]->SetTitle("2p/2h Nieves");
ModeStack[2]->SetFillColor(kRed);
//ModeStack[2]->SetFillStyle(3344);
ModeStack[2]->SetLineColor(kRed);
//ModeStack[11]->SetTitle("#it{#nu + p #rightarrow l^{-} + p + #pi^{+}}");
ModeStack[11]->SetTitle("CC1#pi^{+} on p");
ModeStack[11]->SetFillColor(kGreen);
//ModeStack[11]->SetFillStyle(3004);
ModeStack[11]->SetLineColor(kGreen);
//ModeStack[12]->SetTitle("#it{#nu + n #rightarrow l^{-} + p + #pi^{0}}");
ModeStack[12]->SetTitle("CC1#pi^{0} on n");
ModeStack[12]->SetFillColor(kGreen + 3);
//ModeStack[12]->SetFillStyle(3005);
ModeStack[12]->SetLineColor(kGreen);
//ModeStack[13]->SetTitle("#it{#nu + n #rightarrow l^{-} + n + #pi^{+}}");
ModeStack[13]->SetTitle("CC1#pi^{+} on n");
ModeStack[13]->SetFillColor(kGreen - 2);
//ModeStack[13]->SetFillStyle(3004);
ModeStack[13]->SetLineColor(kGreen);
ModeStack[16]->SetTitle("CC coherent");
ModeStack[16]->SetFillColor(kBlue);
//ModeStack[16]->SetFillStyle(3644);
ModeStack[16]->SetLineColor(kBlue);
//ModeStack[17]->SetTitle("#it{#nu + n #rightarrow l^{-} + p + #gamma}");
ModeStack[17]->SetTitle("CC1#gamma");
ModeStack[17]->SetFillColor(kMagenta);
//ModeStack[17]->SetFillStyle(3001);
ModeStack[17]->SetLineColor(kMagenta);
ModeStack[21]->SetTitle("Multi #pi (1.3 < W < 2.0)");
ModeStack[21]->SetFillColor(kYellow);
//ModeStack[21]->SetFillStyle(3005);
ModeStack[21]->SetLineColor(kYellow);
//ModeStack[22]->SetTitle("#it{#nu + n #rightarrow l^{-} + p + #eta^{0}}");
ModeStack[22]->SetTitle("CC1#eta^{0} on n");
ModeStack[22]->SetFillColor(kYellow - 2);
//ModeStack[22]->SetFillStyle(3013);
ModeStack[22]->SetLineColor(kYellow - 2);
//ModeStack[23]->SetTitle("#it{#nu + n #rightarrow l^{-} + #Lambda + K^{+}}");
ModeStack[23]->SetTitle("CC1#Labda1K^{+}");
ModeStack[23]->SetFillColor(kYellow - 6);
//ModeStack[23]->SetFillStyle(3013);
ModeStack[23]->SetLineColor(kYellow - 6);
ModeStack[26]->SetTitle("DIS (W > 2.0)");
ModeStack[26]->SetFillColor(kRed);
//ModeStack[26]->SetFillStyle(3006);
ModeStack[26]->SetLineColor(kRed);
// NC
//ModeStack[31]->SetTitle("#it{#nu + n #rightarrow #nu + n + #pi^{0}}");
ModeStack[31]->SetTitle("NC1#pi^{0} on n");
ModeStack[31]->SetFillColor(kBlue);
//ModeStack[31]->SetFillStyle(3004);
ModeStack[31]->SetLineColor(kBlue);
//ModeStack[32]->SetTitle("#it{#nu + p #rightarrow #nu + p + #pi^{0}}");
ModeStack[32]->SetTitle("NC1#pi^{0} on p");
ModeStack[32]->SetFillColor(kBlue + 3);
//ModeStack[32]->SetFillStyle(3004);
ModeStack[32]->SetLineColor(kBlue + 3);
//ModeStack[33]->SetTitle("#it{#nu + n #rightarrow #nu + p + #pi^{-}}");
ModeStack[33]->SetTitle("NC1#pi^{-} on n");
ModeStack[33]->SetFillColor(kBlue - 2);
//ModeStack[33]->SetFillStyle(3005);
ModeStack[33]->SetLineColor(kBlue - 2);
//ModeStack[34]->SetTitle("#it{#nu + p #rightarrow #nu + n + #pi^{+}}");
ModeStack[34]->SetTitle("NC1#pi^{+} on p");
ModeStack[34]->SetFillColor(kBlue - 8);
//ModeStack[34]->SetFillStyle(3005);
ModeStack[34]->SetLineColor(kBlue - 8);
ModeStack[36]->SetTitle("NC Coherent");
ModeStack[36]->SetFillColor(kBlue + 8);
//ModeStack[36]->SetFillStyle(3644);
ModeStack[36]->SetLineColor(kBlue + 8);
//ModeStack[38]->SetTitle("#it{#nu + n #rightarrow #nu + n + #gamma}");
ModeStack[38]->SetTitle("NC1#gamma on n");
ModeStack[38]->SetFillColor(kMagenta);
//ModeStack[38]->SetFillStyle(3001);
ModeStack[38]->SetLineColor(kMagenta);
//ModeStack[39]->SetTitle("#it{#nu + p #rightarrow #nu + p + #gamma}");
ModeStack[39]->SetTitle("NC1#gamma on p");
ModeStack[39]->SetFillColor(kMagenta - 10);
//ModeStack[39]->SetFillStyle(3001);
ModeStack[39]->SetLineColor(kMagenta - 10);
ModeStack[41]->SetTitle("Multi #pi (1.3 < W < 2.0)");
ModeStack[41]->SetFillColor(kBlue - 10);
//ModeStack[41]->SetFillStyle(3005);
ModeStack[41]->SetLineColor(kBlue - 10);
//ModeStack[42]->SetTitle("#it{#nu + n #rightarrow #nu + n + #eta^{0}}");
ModeStack[42]->SetTitle("NC1#eta^{0} on n");
ModeStack[42]->SetFillColor(kYellow - 2);
//ModeStack[42]->SetFillStyle(3013);
ModeStack[42]->SetLineColor(kYellow - 2);
//ModeStack[43]->SetTitle("#it{#nu + p #rightarrow #nu + p + #eta^{0}}");
ModeStack[43]->SetTitle("NC1#eta^{0} on p");
ModeStack[43]->SetFillColor(kYellow - 4);
//ModeStack[43]->SetFillStyle(3013);
ModeStack[43]->SetLineColor(kYellow - 4);
//ModeStack[44]->SetTitle("#it{#nu + n #rightarrow #nu + #Lambda + K^{0}}");
ModeStack[44]->SetTitle("NC1#Lambda1K^{0} on n");
ModeStack[44]->SetFillColor(kYellow - 6);
//ModeStack[44]->SetFillStyle(3014);
ModeStack[44]->SetLineColor(kYellow - 6);
//ModeStack[45]->SetTitle("#it{#nu + p #rightarrow #nu + #Lambda + K^{+}}");
ModeStack[45]->SetTitle("NC1#Lambda1K^{+}");
ModeStack[45]->SetFillColor(kYellow - 10);
//ModeStack[45]->SetFillStyle(3014);
ModeStack[45]->SetLineColor(kYellow - 10);
ModeStack[46]->SetTitle("DIS (W > 2.0)");
ModeStack[46]->SetFillColor(kRed);
//ModeStack[46]->SetFillStyle(3006);
ModeStack[46]->SetLineColor(kRed);
//ModeStack[51]->SetTitle("#it{#nu + p #rightarrow #nu + p}");
ModeStack[51]->SetTitle("NC on p");
ModeStack[51]->SetFillColor(kBlack);
//ModeStack[51]->SetFillStyle(3444);
ModeStack[51]->SetLineColor(kBlack);
//ModeStack[52]->SetTitle("#it{#nu + n #rightarrow #nu + n}");
ModeStack[52]->SetTitle("NC on n");
ModeStack[52]->SetFillColor(kGray);
//ModeStack[52]->SetFillStyle(3444);
ModeStack[52]->SetLineColor(kGray);
return allmodes;
};
TLegend PlotUtils::GenerateStackLegend(THStack stack, int xlow, int ylow, int xhigh, int yhigh) {
TLegend leg = TLegend(xlow, ylow, xhigh, yhigh);
TObjArray* histarray = stack.GetStack();
int nhist = histarray->GetEntries();
for (int i = 0; i < nhist; i++) {
TH1* hist = (TH1*)(histarray->At(i));
leg.AddEntry( (hist), ((TH1*)histarray->At(i))->GetTitle(), "fl");
}
leg.SetName(Form("%s_LEG", stack.GetName()));
return leg;
};
void PlotUtils::ScaleNeutModeArray(TH1* hist[], double factor, std::string option) {
for (int i = 0; i < 60; i++) {
if (hist[i])
hist[i]->Scale(factor, option.c_str());
}
return;
};
void PlotUtils::ResetNeutModeArray(TH1* hist[]) {
for (int i = 0; i < 60; i++) {
if (hist[i])
hist[i]->Reset();
}
return;
};
//********************************************************************
// This assumes the Enu axis is the x axis, as is the case for MiniBooNE 2D distributions
void PlotUtils::FluxUnfoldedScaling(TH2D* fMCHist, TH1D* fhist, TH1D* ehist, double scalefactor) {
//********************************************************************
// Make clones to avoid changing stuff
TH1D* eventhist = (TH1D*)ehist->Clone();
TH1D* fFluxHist = (TH1D*)fhist->Clone();
// Undo width integral in SF
fMCHist->Scale( scalefactor / eventhist->Integral(1, eventhist->GetNbinsX() + 1, "width"));
// Standardise The Flux
eventhist->Scale(1.0 / fFluxHist->Integral());
fFluxHist->Scale(1.0 / fFluxHist->Integral());
// Do interpolation for 2D plots?
// fFluxHist = PlotUtils::InterpolateFineHistogram(fFluxHist,100,"width");
// eventhist = PlotUtils::InterpolateFineHistogram(eventhist,100,"width");
// eventhist->Scale(1.0/fFluxHist->Integral());
// fFluxHist->Scale(1.0/fFluxHist->Integral());
// Scale fMCHist by eventhist integral
fMCHist->Scale( eventhist->Integral(1, eventhist->GetNbinsX() + 1) );
// Now Get a flux PDF assuming X axis is Enu
TH1D* pdfflux = (TH1D*) fMCHist->ProjectionX()->Clone();
// pdfflux->Write( (std::string(fMCHist->GetName()) + "_PROJX").c_str());
pdfflux->Reset();
// Awful MiniBooNE Check for the time being
bool ismb = std::string(fMCHist->GetName()).find("MiniBooNE") != std::string::npos;
for (int i = 0; i < pdfflux->GetNbinsX(); i++) {
double Ml = pdfflux->GetXaxis()->GetBinLowEdge(i + 1);
double Mh = pdfflux->GetXaxis()->GetBinLowEdge(i + 2);
//double Mc = pdfflux->GetXaxis()->GetBinCenter(i+1);
//double Mw = pdfflux->GetBinWidth(i+1);
double fluxint = 0.0;
// Scaling to match flux for MB
if (ismb) {
Ml /= 1.E3;
Mh /= 1.E3;
// Mc /= 1.E3;
// Mw /= 1.E3;
}
for (int j = 0; j < fFluxHist->GetNbinsX(); j++) {
//double Fc = fFluxHist->GetXaxis()->GetBinCenter(j+1);
double Fl = fFluxHist->GetXaxis()->GetBinLowEdge(j + 1);
double Fh = fFluxHist->GetXaxis()->GetBinLowEdge(j + 2);
double Fe = fFluxHist->GetBinContent(j + 1);
double Fw = fFluxHist->GetXaxis()->GetBinWidth(j + 1);
if (Fl >= Ml and Fh <= Mh) { fluxint += Fe; }
else if (Fl < Ml and Fl < Mh and Fh > Ml and Fh < Mh) { fluxint += Fe * (Fh - Ml) / Fw; }
else if (Fh > Mh and Fl < Mh and Fh > Ml and Fl > Ml) { fluxint += Fe * (Mh - Fl) / Fw; }
else if (Ml >= Fl and Mh <= Fh) { fluxint += Fe * (Mh - Ml) / Fw; }
else { continue; }
}
pdfflux->SetBinContent(i + 1, fluxint);
}
for (int i = 0; i < fMCHist->GetNbinsX(); i++) {
for (int j = 0; j < fMCHist->GetNbinsY(); j++) {
if (pdfflux->GetBinContent(i + 1) == 0.0) continue;
double binWidth = fMCHist->GetYaxis()->GetBinLowEdge(j + 2) - fMCHist->GetYaxis()->GetBinLowEdge(j + 1);
fMCHist->SetBinContent(i + 1, j + 1, fMCHist->GetBinContent(i + 1, j + 1) / pdfflux->GetBinContent(i + 1) / binWidth);
fMCHist->SetBinError(i + 1, j + 1, fMCHist->GetBinError(i + 1, j + 1) / pdfflux->GetBinContent(i + 1) / binWidth);
}
}
delete eventhist;
delete fFluxHist;
return;
};
TH1D* PlotUtils::InterpolateFineHistogram(TH1D* hist, int res, std::string opt) {
int nbins = hist->GetNbinsX();
double elow = hist->GetXaxis()->GetBinLowEdge(1);
double ehigh = hist->GetXaxis()->GetBinLowEdge(nbins + 1);
bool width = true; //opt.find("width") != std::string::npos;
TH1D* fine = new TH1D("fine", "fine", nbins * res, elow, ehigh);
TGraph* temp = new TGraph();
for (int i = 0; i < nbins; i++) {
double E = hist->GetXaxis()->GetBinCenter(i + 1);
double C = hist->GetBinContent(i + 1);
double W = hist->GetXaxis()->GetBinWidth(i + 1);
if (!width) W = 1.0;
if (W != 0.0) temp->SetPoint( temp->GetN(), E, C / W );
}
for (int i = 0; i < fine->GetNbinsX(); i++) {
double E = fine->GetXaxis()->GetBinCenter(i + 1);
double W = fine->GetBinWidth(i + 1);
if (!width) W = 1.0;
fine->SetBinContent(i + 1, temp->Eval(E, 0, "S") * W);
}
fine->Scale(hist->Integral(1, hist->GetNbinsX() + 1) / fine->Integral(1, fine->GetNbinsX() + 1));
std::cout << "Interpolation Difference = " << fine->Integral(1, fine->GetNbinsX() + 1) << "/" << hist->Integral(1, hist->GetNbinsX() + 1) << std::endl;
return fine;
}
//********************************************************************
// This interpolates the flux by a TGraph instead of requiring the flux and MC flux to have the same binning
void PlotUtils::FluxUnfoldedScaling(TH1D* mcHist, TH1D* fhist, TH1D* ehist, double scalefactor, int nevents) {
//********************************************************************
TH1D* eventhist = (TH1D*)ehist->Clone();
TH1D* fFluxHist = (TH1D*)fhist->Clone();
if (FitPar::Config().GetParB("save_flux_debug")) {
std::string name = std::string(mcHist->GetName());
mcHist->Write((name + "_UNF_MC").c_str());
fFluxHist->Write((name + "_UNF_FLUX").c_str());
eventhist->Write((name + "_UNF_EVT").c_str());
TH1D* scalehist = new TH1D("scalehist", "scalehist", 1, 0.0, 1.0);
scalehist->SetBinContent(1, scalefactor);
scalehist->SetBinContent(2, nevents);
scalehist->Write((name + "_UNF_SCALE").c_str());
}
// Undo width integral in SF
mcHist->Scale( scalefactor / eventhist->Integral(1, eventhist->GetNbinsX() + 1, "width"));
// Standardise The Flux
eventhist->Scale(1.0 / fFluxHist->Integral());
fFluxHist->Scale(1.0 / fFluxHist->Integral());
// Scale mcHist by eventhist integral
mcHist->Scale( eventhist->Integral(1, eventhist->GetNbinsX() + 1) );
// Now Get a flux PDF
TH1D* pdfflux = (TH1D*) mcHist->Clone();
pdfflux->Reset();
for (int i = 0; i < mcHist->GetNbinsX(); i++) {
double Ml = mcHist->GetXaxis()->GetBinLowEdge(i + 1);
double Mh = mcHist->GetXaxis()->GetBinLowEdge(i + 2);
//double Mc = mcHist->GetXaxis()->GetBinCenter(i+1);
//double Me = mcHist->GetBinContent(i+1);
//double Mw = mcHist->GetBinWidth(i+1);
double fluxint = 0.0;
for (int j = 0; j < fFluxHist->GetNbinsX(); j++) {
//double Fc = fFluxHist->GetXaxis()->GetBinCenter(j+1);
double Fl = fFluxHist->GetXaxis()->GetBinLowEdge(j + 1);
double Fh = fFluxHist->GetXaxis()->GetBinLowEdge(j + 2);
double Fe = fFluxHist->GetBinContent(j + 1);
double Fw = fFluxHist->GetXaxis()->GetBinWidth(j + 1);
if (Fl >= Ml and Fh <= Mh) { fluxint += Fe; }
else if (Fl < Ml and Fl < Mh and Fh > Ml and Fh < Mh) { fluxint += Fe * (Fh - Ml) / Fw; }
else if (Fh > Mh and Fl < Mh and Fh > Ml and Fl > Ml) { fluxint += Fe * (Mh - Fl) / Fw; }
else if (Ml >= Fl and Mh <= Fh) { fluxint += Fe * (Mh - Ml) / Fw; }
else { continue; }
}
pdfflux->SetBinContent(i + 1, fluxint);
}
// Scale MC hist by pdfflux
for (int i = 0; i < mcHist->GetNbinsX(); i++) {
if (pdfflux->GetBinContent(i + 1) == 0.0) continue;
mcHist->SetBinContent(i + 1, mcHist->GetBinContent(i + 1) / pdfflux->GetBinContent(i + 1));
mcHist->SetBinError(i + 1, mcHist->GetBinError(i + 1) / pdfflux->GetBinContent(i + 1));
}
delete eventhist;
delete fFluxHist;
return;
};
// MOVE TO GENERAL UTILS
//********************************************************************
void PlotUtils::Set2DHistFromText(std::string dataFile, TH2* hist, double norm, bool skipbins) {
//********************************************************************
std::string line;
std::ifstream data(dataFile.c_str(), ifstream::in);
int yBin = 0;
while (std::getline(data >> std::ws, line, '\n')) {
std::vector entries = GeneralUtils::ParseToDbl(line, " ");
// Loop over entries and insert them into the histogram
for (uint xBin = 0; xBin < entries.size(); xBin++) {
if (!skipbins or entries[xBin] != -1.0)
hist->SetBinContent(xBin + 1, yBin + 1, entries[xBin]*norm);
}
yBin++;
}
return;
}
// MOVE TO GENERAL UTILS
TH1D* PlotUtils::GetTH1DFromFile(std::string dataFile, std::string title, std::string fPlotTitles, std::string alt_name) {
TH1D* tempPlot;
// If format is a root file
if (dataFile.find(".root") != std::string::npos) {
TFile* temp_infile = new TFile(dataFile.c_str(), "READ");
tempPlot = (TH1D*)temp_infile->Get(title.c_str());
tempPlot->SetDirectory(0);
temp_infile->Close();
delete temp_infile;
// Else its a space seperated txt file
} else {
// Make a TGraph Errors
TGraphErrors *gr = new TGraphErrors(dataFile.c_str(), "%lg %lg %lg");
if (gr->IsZombie()) {
exit(-1);
}
double* bins = gr->GetX();
double* values = gr->GetY();
double* errors = gr->GetEY();
int npoints = gr->GetN();
// Fill the histogram from it
tempPlot = new TH1D(title.c_str(), title.c_str(), npoints - 1, bins);
for (int i = 0; i < npoints; ++i) {
tempPlot->SetBinContent(i+1, values[i]);
// If only two columns are present in the input file, use the sqrt(values) as the error
// equivalent to assuming that the error is statistical
if (!errors[i]) tempPlot->SetBinError(i+1, sqrt(values[i]));
else tempPlot->SetBinError(i+1, errors[i]);
}
delete gr;
}
// Allow alternate naming for root files
if (!alt_name.empty()) {
tempPlot->SetNameTitle(alt_name.c_str(), alt_name.c_str());
}
// Allow alternate axis titles
if (!fPlotTitles.empty()) {
tempPlot->SetNameTitle( tempPlot->GetName(), (std::string(tempPlot->GetTitle()) + fPlotTitles).c_str() );
}
return tempPlot;
};
TH1D* PlotUtils::GetRatioPlot(TH1D* hist1, TH1D* hist2) {
// make copy of first hist
TH1D* new_hist = (TH1D*) hist1->Clone();
// Do bins and errors ourselves as scales can go awkward
for (int i = 0; i < new_hist->GetNbinsX(); i++) {
if (hist2->GetBinContent(i + 1) == 0.0) {
new_hist->SetBinContent(i + 1, 0.0);
}
new_hist->SetBinContent(i + 1, hist1->GetBinContent(i + 1) / hist2->GetBinContent(i + 1) );
new_hist->SetBinError(i + 1, hist1->GetBinError(i + 1) / hist2->GetBinContent(i + 1) );
}
return new_hist;
};
TH1D* PlotUtils::GetRenormalisedPlot(TH1D* hist1, TH1D* hist2) {
// make copy of first hist
TH1D* new_hist = (TH1D*) hist1->Clone();
if (hist1->Integral("width") == 0 or hist2->Integral("width") == 0) {
new_hist->Reset();
return new_hist;
}
Double_t scaleF = hist2->Integral("width") / hist1->Integral("width");
new_hist->Scale(scaleF);
return new_hist;
};
TH1D* PlotUtils::GetShapePlot(TH1D* hist1) {
// make copy of first hist
TH1D* new_hist = (TH1D*) hist1->Clone();
if (hist1->Integral("width") == 0) {
new_hist->Reset();
return new_hist;
}
Double_t scaleF1 = 1.0 / hist1->Integral("width");
new_hist->Scale(scaleF1);
return new_hist;
};
TH1D* PlotUtils::GetShapeRatio(TH1D* hist1, TH1D* hist2) {
TH1D* new_hist1 = GetShapePlot(hist1);
TH1D* new_hist2 = GetShapePlot(hist2);
// Do bins and errors ourselves as scales can go awkward
for (int i = 0; i < new_hist1->GetNbinsX(); i++) {
if (hist2->GetBinContent(i + 1) == 0) {
new_hist1->SetBinContent(i + 1, 0.0);
}
new_hist1->SetBinContent(i + 1, new_hist1->GetBinContent(i + 1) / new_hist2->GetBinContent(i + 1) );
new_hist1->SetBinError(i + 1, new_hist1->GetBinError(i + 1) / new_hist2->GetBinContent(i + 1) );
}
delete new_hist2;
return new_hist1;
};
TH2D* PlotUtils::GetCovarPlot(TMatrixDSym* cov, std::string name, std::string title) {
TH2D* CovarPlot;
if (cov) CovarPlot = new TH2D((*cov));
else CovarPlot = new TH2D(name.c_str(), title.c_str(), 1, 0, 1, 1, 0, 1);
CovarPlot->SetName(name.c_str());
CovarPlot->SetTitle(title.c_str());
return CovarPlot;
}
TH2D* PlotUtils::GetFullCovarPlot(TMatrixDSym* cov, std::string name) {
return PlotUtils::GetCovarPlot(cov, name + "_COV", name + "_COV;Bins;Bins;Covariance (#times10^{-76})");
}
TH2D* PlotUtils::GetInvCovarPlot(TMatrixDSym* cov, std::string name) {
return PlotUtils::GetCovarPlot(cov, name + "_INVCOV", name + "_INVCOV;Bins;Bins;Inv. Covariance (#times10^{-76})");
}
TH2D* PlotUtils::GetDecompCovarPlot(TMatrixDSym* cov, std::string name) {
return PlotUtils::GetCovarPlot(cov, name + "_DECCOV", name + "_DECCOV;Bins;Bins;Decomp Covariance (#times10^{-76})");
}
TH1D* PlotUtils::GetTH1DFromRootFile(std::string file, std::string name) {
if (name.empty()) {
std::vector tempfile = GeneralUtils::ParseToStr(file, ";");
file = tempfile[0];
name = tempfile[1];
}
TFile* rootHistFile = new TFile(file.c_str(), "READ");
TH1D* tempHist = (TH1D*) rootHistFile->Get(name.c_str())->Clone();
tempHist->SetDirectory(0);
rootHistFile->Close();
return tempHist;
}
TH2D* PlotUtils::GetTH2DFromRootFile(std::string file, std::string name) {
if (name.empty()) {
std::vector tempfile = GeneralUtils::ParseToStr(file, ";");
file = tempfile[0];
name = tempfile[1];
}
TFile* rootHistFile = new TFile(file.c_str(), "READ");
TH2D* tempHist = (TH2D*) rootHistFile->Get(name.c_str())->Clone();
tempHist->SetDirectory(0);
rootHistFile->Close();
return tempHist;
}
TH2D* PlotUtils::GetTH2DFromTextFile(std::string file) {
/// Contents should be
/// Low Edfe
return NULL;
}
void PlotUtils::AddNeutModeArray(TH1D* hist1[], TH1D* hist2[], double scaling) {
for (int i = 0; i < 60; i++) {
if (!hist2[i]) continue;
if (!hist1[i]) continue;
hist1[i]->Add(hist2[i], scaling);
}
return;
}
void PlotUtils::ScaleToData(TH1D* data, TH1D* mc, TH1I* mask) {
double scaleF = GetDataMCRatio(data, mc, mask);
mc->Scale(scaleF);
return;
}
void PlotUtils::MaskBins(TH1D* hist, TH1I* mask) {
for (int i = 0; i < hist->GetNbinsX(); i++) {
if (mask->GetBinContent(i + 1) <= 0.5) continue;
hist->SetBinContent(i + 1, 0.0);
hist->SetBinError(i + 1, 0.0);
LOG(REC) << "MaskBins: Set " << hist->GetName() << " Bin " << i + 1 << " to 0.0 +- 0.0" << std::endl;
}
return;
}
void PlotUtils::MaskBins(TH2D* hist, TH2I* mask) {
for (int i = 0; i < hist->GetNbinsX(); i++) {
for (int j = 0; j < hist->GetNbinsY(); j++) {
if (mask->GetBinContent(i + 1, j + 1) <= 0.5) continue;
hist->SetBinContent(i + 1, j + 1, 0.0);
hist->SetBinError(i + 1, j + 1, 0.0);
LOG(REC) << "MaskBins: Set " << hist->GetName() << " Bin " << i + 1 << " " << j + 1 << " to 0.0 +- 0.0" << std::endl;
}
}
return;
}
double PlotUtils::GetDataMCRatio(TH1D* data, TH1D* mc, TH1I* mask) {
double rat = 1.0;
TH1D* newmc = (TH1D*)mc->Clone();
TH1D* newdt = (TH1D*)data->Clone();
if (mask) {
MaskBins(newmc, mask);
MaskBins(newdt, mask);
}
rat = newdt->Integral() / newmc->Integral();
return rat;
}
TH2D* PlotUtils::GetCorrelationPlot(TH2D* cov, std::string name) {
TH2D* cor = (TH2D*)cov->Clone();
cor->Reset();
for (int i = 0; i < cov->GetNbinsX(); i++) {
for (int j = 0; j < cov->GetNbinsY(); j++) {
if (cov->GetBinContent(i + 1, i + 1) != 0.0 and cov->GetBinContent(j + 1, j + 1) != 0.0)
cor->SetBinContent(i + 1, j + 1, cov->GetBinContent(i + 1, j + 1) / (sqrt(cov->GetBinContent(i + 1, i + 1) * cov->GetBinContent(j + 1, j + 1))));
}
}
if (!name.empty()) {
cor->SetNameTitle(name.c_str(), (name + ";;correlation").c_str());
}
cor->SetMinimum(-1);
cor->SetMaximum(1);
return cor;
}
TH2D* PlotUtils::GetDecompPlot(TH2D* cov, std::string name) {
TMatrixDSym* covarmat = new TMatrixDSym(cov->GetNbinsX());
for (int i = 0; i < cov->GetNbinsX(); i++)
for (int j = 0; j < cov->GetNbinsY(); j++)
(*covarmat)(i, j) = cov->GetBinContent(i + 1, j + 1);
TMatrixDSym* decompmat = StatUtils::GetDecomp(covarmat);
TH2D* dec = (TH2D*) cov->Clone();
for (int i = 0; i < cov->GetNbinsX(); i++)
for (int j = 0; j < cov->GetNbinsY(); j++)
dec->SetBinContent(i + 1, j + 1, (*decompmat)(i, j));
delete covarmat;
delete decompmat;
dec->SetNameTitle(name.c_str(), (name + ";;;decomposition").c_str());
return dec;
}
TH2D* PlotUtils::MergeIntoTH2D(TH1D* xhist, TH1D* yhist, std::string zname) {
std::vector xedges, yedges;
for (int i = 0; i < xhist->GetNbinsX() + 2; i++) {
xedges.push_back(xhist->GetXaxis()->GetBinLowEdge(i + 1));
}
for (int i = 0; i < yhist->GetNbinsX() + 2; i++) {
yedges.push_back(yhist->GetXaxis()->GetBinLowEdge(i + 1));
}
int nbinsx = xhist->GetNbinsX();
int nbinsy = yhist->GetNbinsX();
std::string name = std::string(xhist->GetName()) + "_vs_" + std::string(yhist->GetName());
std::string titles = ";" + std::string(xhist->GetXaxis()->GetTitle()) + ";" + std::string(yhist->GetXaxis()->GetTitle()) + ";" + zname;
TH2D* newplot = new TH2D(name.c_str(), (name + titles).c_str(), nbinsx, &xedges[0], nbinsy, &yedges[0]);
return newplot;
}
//***************************************************
void PlotUtils::MatchEmptyBins(TH1D* data, TH1D* mc) {
//**************************************************
for (int i = 0; i < data->GetNbinsX(); i++) {
if (data->GetBinContent(i + 1) == 0.0 or data->GetBinError(i + 1) == 0.0)
mc->SetBinContent(i + 1, 0.0);
}
return;
}
//***************************************************
void PlotUtils::MatchEmptyBins(TH2D* data, TH2D* mc) {
//**************************************************
for (int i = 0; i < data->GetNbinsX(); i++) {
for (int j = 0; j < data->GetNbinsY(); j++) {
if (data->GetBinContent(i + 1, j + 1) == 0.0 or data->GetBinError(i + 1, j + 1) == 0.0)
mc->SetBinContent(i + 1, j + 1, 0.0);
}
}
return;
}
//***************************************************
TH1D* PlotUtils::GetProjectionX(TH2D* hist, TH2I* mask) {
//***************************************************
TH2D* maskedhist = StatUtils::ApplyHistogramMasking(hist, mask);
TH1D* hist_X = maskedhist->ProjectionX();
delete maskedhist;
return hist_X;
}
//***************************************************
TH1D* PlotUtils::GetProjectionY(TH2D* hist, TH2I* mask) {
//***************************************************
TH2D* maskedhist = StatUtils::ApplyHistogramMasking(hist, mask);
TH1D* hist_Y = maskedhist->ProjectionY();
delete maskedhist;
return hist_Y;
}