diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt index 1d488e3..1afec7d 100644 --- a/app/CMakeLists.txt +++ b/app/CMakeLists.txt @@ -1,217 +1,225 @@ # Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret ################################################################################ # This file is part of NUISANCE. # # NUISANCE is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # NUISANCE is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with NUISANCE. If not, see . ################################################################################ set(TARGETS_TO_BUILD) if(USE_MINIMIZER) add_executable(nuismin nuismin.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuismin) target_link_libraries(nuismin ${MODULETargets}) target_link_libraries(nuismin ${CMAKE_DEPENDLIB_FLAGS}) # target_link_libraries(nuismin ${ROOT_LIBS}) if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "") set_target_properties(nuismin PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS}) endif() add_executable(nuissplines nuissplines.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuissplines) target_link_libraries(nuissplines ${MODULETargets}) target_link_libraries(nuissplines ${CMAKE_DEPENDLIB_FLAGS}) # target_link_libraries(nuissplines ${ROOT_LIBS}) if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "") set_target_properties(nuissplines PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS}) endif() endif() include_directories(${RWENGINE_INCLUDE_DIRECTORIES}) include_directories(${CMAKE_SOURCE_DIR}/src/Routines) include_directories(${CMAKE_SOURCE_DIR}/src/InputHandler) include_directories(${CMAKE_SOURCE_DIR}/src/Genie) include_directories(${CMAKE_SOURCE_DIR}/src/FitBase) include_directories(${CMAKE_SOURCE_DIR}/src/Statistical) include_directories(${CMAKE_SOURCE_DIR}/src/Utils) include_directories(${CMAKE_SOURCE_DIR}/src/Config) include_directories(${CMAKE_SOURCE_DIR}/src/Logger) include_directories(${CMAKE_SOURCE_DIR}/src/Splines) include_directories(${CMAKE_SOURCE_DIR}/src/Reweight) include_directories(${CMAKE_SOURCE_DIR}/src/FCN) include_directories(${CMAKE_SOURCE_DIR}/src/MCStudies) include_directories(${CMAKE_SOURCE_DIR}/src/Smearceptance) include_directories(${EXP_INCLUDE_DIRECTORIES}) if (USE_NuWro AND NOT NUWRO_BUILT_FROM_FILE) add_executable(nuwro_nuisance nuwro_NUISANCE.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuwro_nuisance) target_link_libraries(nuwro_nuisance ${MODULETargets}) target_link_libraries(nuwro_nuisance ${CMAKE_DEPENDLIB_FLAGS}) # target_link_libraries(nuwro_nuisance ${ROOT_LIBS}) include_directories(${CMAKE_SOURCE_DIR}/src/FitBase) if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "") set_target_properties(nuwro_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS}) endif() endif() if (USE_NEUT) add_executable(neut_nuisance neut_NUISANCE.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};neut_nuisance) target_link_libraries(neut_nuisance ${MODULETargets}) target_link_libraries(neut_nuisance ${CMAKE_DEPENDLIB_FLAGS}) target_link_libraries(neut_nuisance ${ROOT_LIBS}) include_directories(${CMAKE_SOURCE_DIR}/src/FitBase) if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "") set_target_properties(neut_nuisance PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS}) endif() endif() if (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() +add_executable(nuisbac nuisbac.cxx) +set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};nuisbac) +target_link_libraries(nuisbac ${MODULETargets}) +target_link_libraries(nuisbac ${CMAKE_DEPENDLIB_FLAGS}) +if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "") + set_target_properties(nuisbac PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS}) +endif() + install(TARGETS ${TARGETS_TO_BUILD} DESTINATION bin) #add_executable(DumpROOTClassesFromVector DumpROOTClassesFromVector.cxx) # #Strip out -lNuWro_event1 # string(REPLACE "-lNuWro_event1" "" NWEVSTRIPPED_CDF ${CMAKE_DEPENDLIB_FLAGS}) # cmessage(DEBUG "Attempted to strip out nuwro library: \"${CMAKE_DEPENDLIB_FLAGS}\" -> \"${NWEVSTRIPPED_CDF}\"") # add_executable(PrepareNEUT PrepareNEUT.cxx) # target_link_libraries(DumpROOTClassesFromVector ${MODULETargets}) # target_link_libraries(DumpROOTClassesFromVector ${NWEVSTRIPPED_CDF}) # if(NOT CMAKE_LINK_FLAGS STREQUAL "") # set_target_properties(DumpROOTClassesFromVector PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS}) # endif() #install(TARGETS DumpROOTClassesFromVector DESTINATION bin) diff --git a/app/nuisbac.cxx b/app/nuisbac.cxx new file mode 100644 index 0000000..0b6df2f --- /dev/null +++ b/app/nuisbac.cxx @@ -0,0 +1,9 @@ +#include "Initialiser.h" + +int main(){ + + for(size_t i = 0; i < 100; ++i){ + nuisance_init(); + } + +} diff --git a/cmake/ROOTSetup.cmake b/cmake/ROOTSetup.cmake index ebf9dee..37f6157 100644 --- a/cmake/ROOTSetup.cmake +++ b/cmake/ROOTSetup.cmake @@ -1,162 +1,167 @@ # Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret ################################################################################ # This file is part of NUISANCE. # # NUISANCE is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # NUISANCE is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with NUISANCE. If not, see . ################################################################################ if ( NOT DEFINED ENV{ROOTSYS} ) cmessage (FATAL_ERROR "$ROOTSYS is not defined, please set up ROOT first.") else() cmessage(STATUS "Using ROOT installed at $ENV{ROOTSYS}") set(CMAKE_ROOTSYS $ENV{ROOTSYS}) endif() # Get cflags from ROOT execute_process (COMMAND root-config --cflags OUTPUT_VARIABLE ROOT_CXX_FLAGS_RAW OUTPUT_STRIP_TRAILING_WHITESPACE) string(REPLACE " " ";" ROOT_CXX_FLAGS "${ROOT_CXX_FLAGS_RAW}") # Get libdir from ROOT execute_process (COMMAND root-config --libdir OUTPUT_VARIABLE ROOT_LIBDIR OUTPUT_STRIP_TRAILING_WHITESPACE) # Get version from ROOT execute_process (COMMAND root-config --version OUTPUT_VARIABLE ROOT_VERSION OUTPUT_STRIP_TRAILING_WHITESPACE) # Get features from ROOT execute_process (COMMAND root-config --features OUTPUT_VARIABLE ROOT_FEATURES OUTPUT_STRIP_TRAILING_WHITESPACE) LIST(APPEND EXTRA_LINK_DIRS ${ROOT_LIBDIR}) LIST(APPEND ROOT_LIBS Core Cint RIO XMLIO Net Hist Graf Graf3d Gpad Tree Rint Postscript Matrix Physics MathCore Thread EG Geom GenVector) if(USE_MINIMIZER) if("${ROOT_FEATURES}" MATCHES "minuit2") cmessage(STATUS "ROOT built with MINUIT2 support") LIST(APPEND EXTRA_CXX_FLAGS -D__MINUIT2_ENABLED__) else() cmessage(FATAL_ERROR "ROOT built without MINUIT2 support but minimizer functionality requested. Either configure with -DUSE_MINIMIZER=FALSE or use a version of ROOT with MINUIT2 support.") endif() - string(REGEX MATCH "5.34/([0-9]+)" ROOTVERSMATCH ${ROOT_VERSION}) - if(NOT ROOTVERSMATCH OR ${CMAKE_MATCH_1} LESS "19") - cmessage(FATAL_ERROR "ROOT Version: ${ROOT_VERSION} has out of date minimizer interface, but minimizer functionality requested. Please configure with -DUSE_MINIMIZER=FALSE or update to 5.34/19 or greater to enable minimization features.") + string(REGEX MATCH "6\..*+)" ROOTVERSIXMATCH ${ROOT_VERSION}) + if(ROOTVERSIXMATCH) + cmessage(STATUS "Using ROOT6, We are essentiall flying blind here.") + else() + string(REGEX MATCH "5.34/([0-9]+)" ROOTVERSMATCH ${ROOT_VERSION}) + if(NOT ROOTVERSMATCH OR ${CMAKE_MATCH_1} LESS "19") + cmessage(FATAL_ERROR "ROOT Version: ${ROOT_VERSION} has out of date minimizer interface, but minimizer functionality requested. Please configure with -DUSE_MINIMIZER=FALSE or update to 5.34/19 or greater to enable minimization features.") + endif() endif() endif() if("${ROOT_FEATURES}" MATCHES "opengl") cmessage(STATUS "ROOT built with OpenGL support") LIST(APPEND ROOT_LIBS RGL) endif() if(DEFINED NEED_ROOTPYTHIA6 AND NEED_ROOTPYTHIA6) LIST(APPEND ROOT_LIBS EGPythia6 Pythia6) endif() cmessage ( STATUS "[ROOT]: root-config --version: ${ROOT_VERSION} ") cmessage ( STATUS "[ROOT]: root-config --cflags : ${ROOT_CXX_FLAGS} ") cmessage ( STATUS "[ROOT]: root-config --libs : ${ROOT_LD_FLAGS} ") LIST(APPEND EXTRA_CXX_FLAGS ${ROOT_CXX_FLAGS}) #Helper functions for building dictionaries function(GenROOTDictionary OutputDictName Header LinkDef) get_directory_property(incdirs INCLUDE_DIRECTORIES) string(REPLACE ";" ";-I" LISTDIRINCLUDES "-I${incdirs}") string(REPLACE " " ";" LISTCPPFLAGS "${EXTRA_CXX_FLAGS}") #ROOT5 CINT cannot handle it. list(REMOVE_ITEM LISTCPPFLAGS "-std=c++11") message(STATUS "LISTCPPFLAGS: ${LISTCPPFLAGS}") message(STATUS "LISTINCLUDES: ${LISTDIRINCLUDES}") #Learn how to generate the Dict.cxx and Dict.hxx add_custom_command( OUTPUT "${OutputDictName}.cxx" "${OutputDictName}.h" COMMAND rootcint ARGS -f ${OutputDictName}.cxx -c -p ${LISTDIRINCLUDES} ${LISTCPPFLAGS} ${Header} ${LinkDef} DEPENDS ${Header};${LinkDef}) endfunction() function(BuildROOTProject ProjectName InputFile CommaSeparatedClassesToDump LIBLINKMODE) string(REPLACE "," ";" HeadersToDump ${CommaSeparatedClassesToDump}) set(OUTPUTFILES ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectSource.cxx ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}LinkDef.h ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectHeaders.h ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectInstances.h) cmessage(STATUS "As part of ROOT project: ${ProjectName}") foreach (header ${HeadersToDump}) LIST(APPEND OUTPUTFILES "${CMAKE_BINARY_DIR}/${ProjectName}/${header}.h") cmessage(STATUS "Will generate: ${CMAKE_BINARY_DIR}/${ProjectName}/${header}.h") endforeach() add_custom_command( OUTPUT ${OUTPUTFILES} COMMAND ${CMAKE_BINARY_DIR}/src/Utils/DumpROOTClassesFromVector ARGS ${InputFile} ${CMAKE_BINARY_DIR}/${ProjectName} ${CommaSeparatedClassesToDump} VERBATIM DEPENDS DumpROOTClassesFromVector) add_custom_target(${ProjectName}_sources DEPENDS ${OUTPUTFILES}) GenROOTDictionary( ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectDict ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectHeaders.h ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}LinkDef.h ) add_custom_target(${ProjectName}ProjectDict DEPENDS ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectDict.cxx ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectDict.h ) # add_dependencies(${ProjectName}ProjectDict ${ProjectName}_sources) #ProjectSource.cxx includes ProjectDict.cxx, so no need to add to compilation. set(ROAA_SOURCEFILES ${CMAKE_BINARY_DIR}/${ProjectName}/${ProjectName}ProjectSource.cxx) add_library(${ProjectName} ${LIBLINKMODE} ${ROAA_SOURCEFILES}) add_dependencies(${ProjectName} ${ProjectName}ProjectDict) endfunction() diff --git a/parameters/dial_conversion.card b/parameters/dial_conversion.card index b34a7d9..0227a07 100644 --- a/parameters/dial_conversion.card +++ b/parameters/dial_conversion.card @@ -1,23 +1,38 @@ # par Name Units Nominal FracErr ConvFunc neut_parameter MaCCQE GeV 1.21*(1.0+x*0.16) +neut_parameter MaNFFRES GeV 0.95*(1.0+x*0.157894737) +neut_parameter CA5RES x 1.01*(1.0+x*0.247524752) +neut_parameter BgSclRES x 1.30*(1.0+x*0.153846154) + +neut_parameter FrAbs_pi x 1.1*(1.0+x*0.5) +neut_parameter FrInelLow_pi x 1*(1.0+x*0.5) +neut_parameter FrInelHigh_pi x 1.8*(1.0+x*0.3) +neut_parameter FrPiProd_pi x 1*(1.0+x*0.5) +neut_parameter FrCExLow_pi x 1*(1.0+x*0.5) +neut_parameter FrCExHigh_pi x 1.8*(1.0+x*0.3) + niwg_parameter NIWGMEC_Norm_C12 % 100.0*(1.0+x) niwg_parameter VecFFCCQE MDLQE x niwg_parameter CCQEFermiSurfMom MeV 217*(1.0+x*0.15) - t2k_parameter NIWG2014a_pF_C12 MeV 217*(1.0+x*0.073733) t2k_parameter NIWG2014a_pF_O16 MeV 225*(1.0+x*0.071111) t2k_parameter NIWGMEC_PDDWeight_C12 x 0.5*(1.0+x*1.0) t2k_parameter NIWGMEC_PDDWeight_O16 x 0.5*(1.0+x*1.0) t2k_parameter NXSec_MaCCQE GeV 1.21*(1.0+x*0.165289256) t2k_parameter NXSec_MaNFFRES GeV 0.95*(1.0+x*0.157894737) t2k_parameter NXSec_CA5RES x 1.01*(1.0+x*0.247524752) t2k_parameter NXSec_BgSclCCRES x 1.30*(1.0+x*0.153846154) t2k_parameter NIWG_Effective_rpaCCQE_A x 1.0*(1.0 + x*0.1) t2k_parameter NIWG_Effective_rpaCCQE_B x 1.0*(1.0 + x*0.1) t2k_parameter NIWG_Effective_rpaCCQE_C x 1.0*(1.0 + x*0.1) t2k_parameter NIWG_Effective_rpaCCQE_D x 1.0*(1.0 + x*0.1) t2k_parameter NIWG_Effective_rpaCCQE_U x 1.0*(1.0 + x*0.1) + +nuwro_parameter kNuwro_Ma_CCQE MeV 1200*(1.0+x*0.160) +nuwro_parameter kNuwro_MaRES GeV 0.940*(1.0+x*0.1) +nuwro_parameter kNuwro_CA5 x 1.19*(1.0+x*0.1) +nuwro_parameter kNuwro_SPPBkgScale x 1.0*(1.0+x*0.26) diff --git a/src/Electron/CLAS6-EG2_Accepter.cxx b/src/Electron/CLAS6-EG2_Accepter.cxx index d2c4307..3b910ee 100644 --- a/src/Electron/CLAS6-EG2_Accepter.cxx +++ b/src/Electron/CLAS6-EG2_Accepter.cxx @@ -1,177 +1,252 @@ #include "ISmearcepter.h" #include "TH3D.h" #include "TRandom3.h" #include -#define DEBUG_CLASACCEPT +// #define DEBUG_CLASACCEPT 1 + +struct EffMap { + TH3D *Generated; + TH3D *Accepted; + + void Build(TFile *inpF, std::string const &GenName, + std::string const &AccName) { + Generated = dynamic_cast(inpF->Get(GenName.c_str())); + Accepted = dynamic_cast(inpF->Get(AccName.c_str())); + + if (!Generated) { + std::cout << "[ERROR]: Could not retrieve \"accepted\" histogram: \"" + << AccName << "\" from file: \"" << inpF->GetName() << "\"." + << std::endl; + exit(1); + } + if (!Accepted) { + std::cout << "[ERROR]: Could not retrieve \"generated\" histogram: \"" + << GenName << "\" from file: \"" << inpF->GetName() << "\"." + << std::endl; + exit(1); + } + + Generated = static_cast(Generated->Clone()); + Generated->SetDirectory(NULL); + Accepted = static_cast(Accepted->Clone()); + Accepted->SetDirectory(NULL); + } + + double GetAccRatio(double p_GeV, double costheta, double phi_deg) { + // For a bin in phase space defined by p, cost, phi: + // Find number of generated events + Int_t pbin = Generated->GetXaxis()->FindBin(p_GeV); + + Int_t tbin = Generated->GetYaxis()->FindBin(costheta); + Int_t phibin = Generated->GetZaxis()->FindBin(phi_deg); + + if (((pbin == 0) || (pbin == (Generated->GetXaxis()->GetNbins() + 1))) || + ((tbin == 0) || (tbin == (Generated->GetYaxis()->GetNbins() + 1))) || + ((phibin == 0) || + (phibin == (Generated->GetZaxis()->GetNbins() + 1)))) { + return 0; + } + + double num_gen = Generated->GetBinContent(pbin, tbin, phibin); + // Find number of accepted events + pbin = Accepted->GetXaxis()->FindBin(p_GeV); + tbin = Accepted->GetYaxis()->FindBin(costheta); + phibin = Accepted->GetZaxis()->FindBin(phi_deg); + double num_acc = Accepted->GetBinContent(pbin, tbin, phibin); + double acc_ratio = double(num_acc) / double(num_gen); + + if (((pbin == 0) || (pbin == (Accepted->GetXaxis()->GetNbins() + 1))) || + ((tbin == 0) || (tbin == (Accepted->GetYaxis()->GetNbins() + 1))) || + ((phibin == 0) || (phibin == (Accepted->GetZaxis()->GetNbins() + 1)))) { + return 0; + } + + if ((acc_ratio != 0 && !std::isnormal(acc_ratio)) || (acc_ratio > 1)) { + std::cout << "[BINS]: p " << Generated->GetXaxis()->GetBinLowEdge(1) + << " -- " + << Generated->GetXaxis()->GetBinUpEdge( + Generated->GetXaxis()->GetNbins()) + << ", cost " << Generated->GetYaxis()->GetBinLowEdge(1) + << " -- " + << Generated->GetYaxis()->GetBinUpEdge( + Generated->GetYaxis()->GetNbins()) + << ", phi " << Generated->GetZaxis()->GetBinLowEdge(1) << " -- " + << Generated->GetZaxis()->GetBinUpEdge( + Generated->GetZaxis()->GetNbins()) + << ". " << std::endl + << "[ERROR]: Bad acceptance ratio: " << acc_ratio << " = " + << num_acc << " / " << num_gen << ". (" << p_GeV << ", " + << costheta << ", " << phi_deg << ")." << std::endl; + exit(1); + } + return acc_ratio; + } +}; class CLASAccepter : public ISmearcepter { TRandom3 rand; // Maps a particle PDG to the relevant generated and accepted histograms from // the input map. - std::map > Acceptance; + std::map Acceptance; public: CLASAccepter() { ElementName = "CLASAccepter"; } void SpecifcSetup(nuiskey &nk) { rand.~TRandom3(); new (&rand) TRandom3(); InstanceName = nk.GetS("name"); std::string const &mapfile = nk.GetS("map"); if (!mapfile.length()) { std::cout << "[ERROR]: No input file specified by \"map\" attribute." << std::endl; exit(1); } TFile *f = new TFile(mapfile.c_str()); if (!f || !f->IsOpen()) { std::cout << "[ERROR]: Could not open root file specified by \"map\" " "attribute: \"" << mapfile << "\"" << std::endl; exit(1); } std::vector accepts = nk.GetListOfChildNodes("accept"); for (auto &acc : accepts) { std::string const &genStr = acc.GetS("generated"); std::string const &accStr = acc.GetS("accepted"); if (!genStr.length() || !accStr.length()) { std::cout << "[ERROR]: expected accept node to contain both " "\"generated\" and \"accepted\" attributes." << std::endl; exit(1); } std::string const &pdgs_s = acc.GetS("PDG"); std::vector pdgs_i = GeneralUtils::ParseToInt(pdgs_s, ","); if (!pdgs_i.size()) { std::cout << "[ERROR]: Could not find any applicable particle PDG codes." << std::endl; exit(1); } - std::pair genacc; - - genacc.first = dynamic_cast(f->Get(accStr.c_str())); - genacc.second = dynamic_cast(f->Get(genStr.c_str())); - - if (!genacc.first) { - std::cout << "[ERROR]: Could not retrieve \"accepted\" histogram: \"" - << accStr << "\" from file: \"" << mapfile << "\"." - << std::endl; - exit(1); - } - if (!genacc.second) { - std::cout << "[ERROR]: Could not retrieve \"accepted\" histogram: \"" - << genStr << "\" from file: \"" << mapfile << "\"." - << std::endl; - exit(1); - } + EffMap ef; - genacc.first = static_cast(genacc.first->Clone()); - genacc.first->SetDirectory(NULL); - genacc.second = static_cast(genacc.second->Clone()); - genacc.second->SetDirectory(NULL); + ef.Build(f, genStr, accStr); for (size_t pdg_it = 0; pdg_it < pdgs_i.size(); ++pdg_it) { if (Acceptance.count(pdgs_i[pdg_it])) { std::cout << "[WARN]: Acceptance map already contains acceptance for PDG: " << pdgs_i[pdg_it] << ". Overwriting..." << std::endl; } - Acceptance[pdgs_i[pdg_it]] = genacc; + Acceptance[pdgs_i[pdg_it]] = ef; } } std::cout << "Loaded " << Acceptance.size() << " particle acceptance definitions." << std::endl; f->Close(); delete f; } RecoInfo *Smearcept(FitEvent *fe) { RecoInfo *ri = new RecoInfo(); for (size_t p_it = 0; p_it < fe->NParticles(); ++p_it) { FitParticle *fp = fe->GetParticle(p_it); int PDG = fp->PDG(); - double p = fp->P3().Mag(); + double p = fp->P3().Mag() * 1E-3; double cost = fp->P3().CosTheta(); - double phi = fp->P3().Phi(); + double phi = fp->P3().Phi() * (180.0 / TMath::Pi()) + 150.0; -#ifdef DEBUG_CLASACCEPT +#if DEBUG_CLASACCEPT > 1 std::cout << std::endl; std::cout << "[" << p_it << "]: " << PDG << ", " << fp->Status() << ", " << fp->E() << " -- KE:" << fp->KE() << " Mom: " << p << std::flush; #endif if (fp->Status() != kFinalState) { -#ifdef DEBUG_CLASACCEPT +#if DEBUG_CLASACCEPT > 1 std::cout << " -- Not final state." << std::flush; #endif continue; } if (!Acceptance.count(PDG)) { -#ifdef DEBUG_CLASACCEPT +#if DEBUG_CLASACCEPT > 1 std::cout << " -- Unknown acceptance." << std::flush; #endif continue; } - std::pair acc = Acceptance[PDG]; - - // For a bin in phase space defined by p, cost, phi: - // Find number of generated events - Int_t pbin = acc.second->GetXaxis()->FindBin(p); - Int_t tbin = acc.second->GetYaxis()->FindBin(cost); - Int_t phibin = acc.second->GetZaxis()->FindBin(phi); - double num_gen = acc.second->GetBinContent(pbin, tbin, phibin); - // Find number of accepted events - pbin = acc.first->GetXaxis()->FindBin(p); - tbin = acc.first->GetYaxis()->FindBin(cost); - phibin = acc.first->GetZaxis()->FindBin(phi); - double num_acc = acc.first->GetBinContent(pbin, tbin, phibin); - double acc_ratio = double(num_acc) / double(num_gen); + EffMap eff = Acceptance[PDG]; + + double acc_ratio = eff.GetAccRatio(p, cost, phi); bool accepted = (rand.Uniform() < acc_ratio); if (accepted) { #ifdef DEBUG_CLASACCEPT + std::cout << "(" << p << ", " << cost << ", " << phi << ")." + << std::endl; std::cout << " -- Reconstructed with probability: " << acc_ratio << std::flush; #endif ri->RecObjMom.push_back(fp->P3()); ri->RecObjClass.push_back(fp->PDG()); continue; } #ifdef DEBUG_CLASACCEPT + std::cout << "(" << p << ", " << cost << ", " << phi << ")." << std::endl; std::cout << " -- Rejected with probability: " << acc_ratio << std::flush; #endif -#ifdef DEBUG_CLASACCEPT - std::cout << std::endl; -#endif } #ifdef DEBUG_CLASACCEPT - std::cout << "Reconstructed " << ri->RecObjMom.size() << " particles. " - << std::endl; + std::cout << std::endl; + + if (ri->RecObjMom.size()) { + std::cout << "Reconstructed " << ri->RecObjMom.size() << " particles. " + << std::endl; + } #endif return ri; } + + double GetEfficiency(FitEvent *fe) { + double effweight = 1; + for (size_t p_it = 0; p_it < fe->NParticles(); ++p_it) { + FitParticle *fp = fe->GetParticle(p_it); + + int PDG = fp->PDG(); + double p = fp->P3().Mag() * 1E-3; + double cost = fp->P3().CosTheta(); + double phi = fp->P3().Phi() * (180.0 / TMath::Pi()) + 150.0; + if (fp->Status() != kFinalState) { + continue; + } + if (!Acceptance.count(PDG)) { + continue; + } + EffMap eff = Acceptance[PDG]; + + effweight *= eff.GetAccRatio(p, cost, phi); + } + return effweight; + } }; diff --git a/src/FitBase/CustomVariableBoxes.h b/src/FitBase/CustomVariableBoxes.h index 8bb8b41..e900871 100644 --- a/src/FitBase/CustomVariableBoxes.h +++ b/src/FitBase/CustomVariableBoxes.h @@ -1,21 +1,28 @@ #ifndef CUSTOM_VARIABLES_BOX_H #define CUSTOM_VARIABLES_BOX_H #include "MeasurementVariableBox.h" #include "MeasurementVariableBox1D.h" #include "MeasurementVariableBox2D.h" /*! * \addtogroup FitBase * @{ */ /// Custom box used to also save Q2 for each event. class Q2VariableBox1D : public MeasurementVariableBox1D { -public: - inline Q2VariableBox1D() { Reset(); }; - inline void Reset() { fQ2 = -999.9; } - double fQ2; + public: + inline Q2VariableBox1D() { Reset(); }; + inline MeasurementVariableBox* CloneSignalBox() { + Q2VariableBox1D* box = new Q2VariableBox1D(); + box->fX = this->fX; + box->fSampleWeight = this->fSampleWeight; + box->fQ2 = this->fQ2; + return box; + }; + inline void Reset() { fQ2 = -999.9; } + double fQ2; }; #endif diff --git a/src/InputHandler/GIBUUInputHandler.cxx b/src/InputHandler/GIBUUInputHandler.cxx index c0e8724..f61559f 100644 --- a/src/InputHandler/GIBUUInputHandler.cxx +++ b/src/InputHandler/GIBUUInputHandler.cxx @@ -1,301 +1,302 @@ #ifdef __GiBUU_ENABLED__ #include "GIBUUInputHandler.h" #include "InputUtils.h" GIBUUGeneratorInfo::~GIBUUGeneratorInfo() { DeallocateParticleStack(); } void GIBUUGeneratorInfo::AddBranchesToTree(TTree* tn) { // tn->Branch("NEUTParticleN", fNEUTParticleN, "NEUTParticleN/I"); // tn->Branch("NEUTParticleStatusCode", fNEUTParticleStatusCode, // "NEUTParticleStatusCode[NEUTParticleN]/I"); // tn->Branch("NEUTParticleAliveCode", fNEUTParticleAliveCode, // "NEUTParticleAliveCode[NEUTParticleN]/I"); } void GIBUUGeneratorInfo::SetBranchesFromTree(TTree* tn) { // tn->SetBranchAddress("NEUTParticleN", &fNEUTParticleN ); // tn->SetBranchAddress("NEUTParticleStatusCode", &fNEUTParticleStatusCode ); // tn->SetBranchAddress("NEUTParticleAliveCode", &fNEUTParticleAliveCode ); } void GIBUUGeneratorInfo::AllocateParticleStack(int stacksize) { // fNEUTParticleN = 0; // fNEUTParticleStatusCode = new int[stacksize]; // fNEUTParticleStatusCode = new int[stacksize]; } void GIBUUGeneratorInfo::DeallocateParticleStack() { // delete fNEUTParticleStatusCode; // delete fNEUTParticleAliveCode; } void GIBUUGeneratorInfo::FillGeneratorInfo(GiBUUStdHepReader* nevent) { Reset(); // for (int i = 0; i < nevent->Npart(); i++) { // fNEUTParticleStatusCode[i] = nevent->PartInfo(i)->fStatus; // fNEUTParticleAliveCode[i] = nevent->PartInfo(i)->fIsAlive; // fNEUTParticleN++; // } } void GIBUUGeneratorInfo::Reset() { // for (int i = 0; i < fNEUTParticleN; i++) { // fNEUTParticleStatusCode[i] = -1; // fNEUTParticleAliveCode[i] = 9; // } // fNEUTParticleN = 0; } GIBUUInputHandler::GIBUUInputHandler(std::string const& handle, std::string const& rawinputs) { LOG(SAM) << "Creating GiBUUInputHandler : " << handle << std::endl; // Run a joint input handling fName = handle; fEventType = kGiBUU; fGIBUUTree = new TChain("giRooTracker"); // Loop over all inputs and grab flux, eventhist, and nevents std::vector inputs = InputUtils::ParseInputFileList(rawinputs); for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) { // Open File for histogram access LOG(SAM) << "Opening event file " << inputs[inp_it] << std::endl; TFile* inp_file = new TFile(inputs[inp_it].c_str(), "READ"); if ((!inp_file) || (!inp_file->IsOpen())) { THROW("GiBUU file !IsOpen() at : '" << inputs[inp_it] << "'" << std::endl << "Check that your file paths are correct and the file exists!"); } int NFluxes = bool(dynamic_cast(inp_file->Get("numu_flux"))) + bool(dynamic_cast(inp_file->Get("numub_flux"))) + bool(dynamic_cast(inp_file->Get("nue_flux"))) + bool(dynamic_cast(inp_file->Get("nueb_flux"))) + bool(dynamic_cast(inp_file->Get("e_flux"))); if (NFluxes != 1) { THROW("Found " << NFluxes << " input fluxes in " << inputs[inp_it] << ". The NUISANCE GiBUU interface expects to be " "passed multiple species vectors as separate " "input files like: " "\"GiBUU:(MINERVA_FHC_numu_evts.root,MINERVA_FHC_" "numubar_evts.root,[...])\""); } // Get Flux/Event hist TH1D* fluxhist = dynamic_cast(inp_file->Get("flux")); TH1D* eventhist = dynamic_cast(inp_file->Get("evt")); if (!fluxhist || !eventhist) { ERROR(FTL, "Input File Contents: " << inputs[inp_it]); inp_file->ls(); THROW( "GiBUU FILE doesn't contain flux/xsec info. You may have to " "regenerate your MC!"); } // Get N Events TTree* giRooTracker = dynamic_cast(inp_file->Get("giRooTracker")); if (!giRooTracker) { ERROR(FTL, "giRooTracker Tree not located in NEUT file: " << inputs[inp_it]); THROW("Check your inputs, they may need to be completely regenerated!"); throw; } int nevents = giRooTracker->GetEntries(); if (nevents <= 0) { THROW("Trying to a TTree with " << nevents << " to TChain from : " << inputs[inp_it]); } // Register input to form flux/event rate hists RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist); // Add To TChain fGIBUUTree->AddFile(inputs[inp_it].c_str()); } // Registor all our file inputs SetupJointInputs(); // Create Fit Event fNUISANCEEvent = new FitEvent(); fGiReader = new GiBUUStdHepReader(); fGiReader->SetBranchAddresses(fGIBUUTree); fNUISANCEEvent->HardReset(); }; FitEvent* GIBUUInputHandler::GetNuisanceEvent(const UInt_t entry, const bool lightweight) { // Check out of bounds if (entry >= (UInt_t)fNEvents) return NULL; // Read Entry from TTree to fill NEUT Vect in BaseFitEvt; fGIBUUTree->GetEntry(entry); // Run NUISANCE Vector Filler if (!lightweight) { CalcNUISANCEKinematics(); } #ifdef __PROB3PP_ENABLED__ else { for (int i = 0; i < fGiReader->StdHepN; i++) { int state = GetGIBUUParticleStatus(fGiReader->StdHepStatus[i], fGiReader->StdHepPdg[i]); if (state != kInitialState) { continue; } if (std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4, fGiReader->StdHepPdg[i])) { fNUISANCEEvent->probe_E = fGiReader->StdHepP4[i][3] * 1.E3; fNUISANCEEvent->probe_pdg = fGiReader->StdHepPdg[i]; break; } } } #endif fNUISANCEEvent->InputWeight *= GetInputWeight(entry); + fNUISANCEEvent->GiRead = fGiReader; return fNUISANCEEvent; } int GetGIBUUParticleStatus(int status, int pdg) { int state = kUndefinedState; switch (status) { case 0: // Incoming case 11: // Struck nucleon state = kInitialState; break; case 1: // Good Final State state = kFinalState; break; default: // Other break; } // Set Nuclear States Flag if (pdg > 1000000) { if (state == kInitialState) state = kNuclearInitial; else if (state == kFinalState) state = kNuclearRemnant; else state = kUndefinedState; } return state; } void GIBUUInputHandler::CalcNUISANCEKinematics() { // Reset all variables fNUISANCEEvent->ResetEvent(); FitEvent* evt = fNUISANCEEvent; evt->Mode = fGiReader->GiBUU2NeutCode; evt->fEventNo = 0.0; evt->fTotCrs = 0; evt->fTargetA = 0.0; // Change to get these from nuclear remnant. evt->fTargetZ = 0.0; evt->fTargetH = 0; evt->fBound = 0.0; // Extra GiBUU Input Weight evt->InputWeight = fGiReader->EvtWght; // Check Stack N int npart = fGiReader->StdHepN; int kmax = evt->kMaxParticles; if ((UInt_t)npart > (UInt_t)kmax) { ERROR(WRN, "GiBUU has too many particles. Expanding Stack."); fNUISANCEEvent->ExpandParticleStack(npart); } // Create Stack evt->fNParticles = 0; for (int i = 0; i < npart; i++) { // State int state = GetGIBUUParticleStatus(fGiReader->StdHepStatus[i], fGiReader->StdHepPdg[i]); int curpart = evt->fNParticles; // Set State evt->fParticleState[evt->fNParticles] = state; // Mom evt->fParticleMom[curpart][0] = fGiReader->StdHepP4[i][0] * 1.E3; evt->fParticleMom[curpart][1] = fGiReader->StdHepP4[i][1] * 1.E3; evt->fParticleMom[curpart][2] = fGiReader->StdHepP4[i][2] * 1.E3; evt->fParticleMom[curpart][3] = fGiReader->StdHepP4[i][3] * 1.E3; // PDG evt->fParticlePDG[curpart] = fGiReader->StdHepPdg[i]; // Add to total particles evt->fNParticles++; } // Run Initial, FSI, Final, Other ordering. fNUISANCEEvent->OrderStack(); FitParticle* ISNeutralLepton = fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos); if (ISNeutralLepton) { fNUISANCEEvent->probe_E = ISNeutralLepton->E(); fNUISANCEEvent->probe_pdg = ISNeutralLepton->PDG(); } return; } void GIBUUInputHandler::Print() {} void GIBUUInputHandler::SetupJointInputs() { if (jointeventinputs.size() <= 1) { jointinput = false; } else if (jointeventinputs.size() > 1) { jointinput = true; jointindexswitch = 0; } fMaxEvents = FitPar::Config().GetParI("MAXEVENTS"); if (fMaxEvents != -1 and jointeventinputs.size() > 1) { THROW("Can only handle joint inputs when config MAXEVENTS = -1!"); } for (size_t i = 0; i < jointeventinputs.size(); i++) { double scale = double(fNEvents) / fEventHist->Integral("width"); scale *= jointfluxinputs.at(i)->Integral("width"); jointindexscale.push_back(scale); } fEventHist->SetNameTitle((fName + "_EVT").c_str(), (fName + "_EVT").c_str()); fFluxHist->SetNameTitle((fName + "_FLUX").c_str(), (fName + "_FLUX").c_str()); // Setup Max Events if (fMaxEvents > 1 && fMaxEvents < fNEvents) { if (LOG_LEVEL(SAM)) { std::cout << "\t\t|-> Read Max Entries : " << fMaxEvents << std::endl; } fNEvents = fMaxEvents; } // Print out Status if (LOG_LEVEL(SAM)) { std::cout << "\t\t|-> Total Entries : " << fNEvents << std::endl << "\t\t|-> Event Integral : " << fEventHist->Integral("width") * 1.E-38 << " events/nucleon" << std::endl << "\t\t|-> Flux Integral : " << fFluxHist->Integral("width") << " /cm2" << std::endl << "\t\t|-> Event/Flux : " << fEventHist->Integral("width") * 1.E-38 / fFluxHist->Integral("width") << " cm2/nucleon" << std::endl; } } #endif diff --git a/src/InputHandler/NuWroInputHandler.cxx b/src/InputHandler/NuWroInputHandler.cxx index f4c62ee..384a326 100644 --- a/src/InputHandler/NuWroInputHandler.cxx +++ b/src/InputHandler/NuWroInputHandler.cxx @@ -1,473 +1,479 @@ #ifdef __NUWRO_ENABLED__ #include "NuWroInputHandler.h" #include "InputUtils.h" NuWroGeneratorInfo::~NuWroGeneratorInfo() { delete fNuWroParticlePDGs; } void NuWroGeneratorInfo::AddBranchesToTree(TTree* tn) { tn->Branch("NuWroParticlePDGs", &fNuWroParticlePDGs, "NuWroParticlePDGs/I"); } void NuWroGeneratorInfo::SetBranchesFromTree(TTree* tn) { tn->SetBranchAddress("NuWroParticlePDGs", &fNuWroParticlePDGs); } void NuWroGeneratorInfo::AllocateParticleStack(int stacksize) { fNuWroParticlePDGs = new int[stacksize]; } void NuWroGeneratorInfo::DeallocateParticleStack() { delete fNuWroParticlePDGs; } void NuWroGeneratorInfo::FillGeneratorInfo(event* e) { Reset(); } void NuWroGeneratorInfo::Reset() { for (int i = 0; i < kMaxParticles; i++) { fNuWroParticlePDGs[i] = 0; } } int event1_nof(event* e, int pdg) { int c = 0; for (size_t i = 0; i < e->out.size(); i++) if (e->out[i].pdg == pdg) c++; return c; } NuWroInputHandler::NuWroInputHandler(std::string const& handle, std::string const& rawinputs) { LOG(SAM) << "Creating NuWroInputHandler : " << handle << std::endl; // Run a joint input handling fName = handle; fMaxEvents = FitPar::Config().GetParI("MAXEVENTS"); fSaveExtra = false; // FitPar::Config().GetParB("NuWroSaveExtra"); // Setup the TChain fNuWroTree = new TChain("treeout"); // Loop over all inputs and grab flux, eventhist, and nevents std::vector inputs = InputUtils::ParseInputFileList(rawinputs); for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) { // Open File for histogram access TFile* inp_file = new TFile(inputs[inp_it].c_str(), "READ"); if (!inp_file or inp_file->IsZombie()) { ERR(FTL) << "nuwro File IsZombie() at " << inputs[inp_it] << std::endl; throw; } // Get Flux/Event hist TH1D* fluxhist = (TH1D*)inp_file->Get( (PlotUtils::GetObjectWithName(inp_file, "FluxHist")).c_str()); TH1D* eventhist = (TH1D*)inp_file->Get( (PlotUtils::GetObjectWithName(inp_file, "EvtHist")).c_str()); if (!fluxhist or !eventhist) { ERR(FTL) << "nuwro FILE doesn't contain flux/xsec info" << std::endl; if (FitPar::Config().GetParB("regennuwro")) { ERR(FTL) << "Regen NuWro has not been added yet. Email the developers!" << std::endl; // ProcessNuWroInputFlux(inputs[inp_it]); throw; } else { ERR(FTL) << "If you would like NUISANCE to generate these for you " << "please set parameter regennuwro=1 and re-run." << std::endl; throw; } } // Get N Events TTree* nuwrotree = (TTree*)inp_file->Get("treeout"); if (!nuwrotree) { ERR(FTL) << "treeout not located in nuwro file! " << inputs[inp_it] << std::endl; throw; } int nevents = nuwrotree->GetEntries(); // Register input to form flux/event rate hists RegisterJointInput(inputs[inp_it], nevents, fluxhist, eventhist); // Add to TChain fNuWroTree->Add(inputs[inp_it].c_str()); } // Registor all our file inputs SetupJointInputs(); // Setup Events fNuWroEvent = NULL; fNuWroTree->SetBranchAddress("e", &fNuWroEvent); fNuWroTree->GetEntry(0); fNUISANCEEvent = new FitEvent(); fNUISANCEEvent->fType = kNUWRO; fNUISANCEEvent->fNuwroEvent = fNuWroEvent; fNUISANCEEvent->HardReset(); if (fSaveExtra) { fNuWroInfo = new NuWroGeneratorInfo(); fNUISANCEEvent->AddGeneratorInfo(fNuWroInfo); } }; NuWroInputHandler::~NuWroInputHandler() { if (fNuWroTree) delete fNuWroTree; } void NuWroInputHandler::CreateCache() { // fNuWroTree->SetCacheEntryRange(0, fNEvents); // fNuWroTree->AddBranchToCache("*", 1); // fNuWroTree->SetCacheSize(fCacheSize); } void NuWroInputHandler::RemoveCache() { // fNuWroTree->SetCacheEntryRange(0, fNEvents); // fNuWroTree->AddBranchToCache("*", 0); // fNuWroTree->SetCacheSize(0); } void NuWroInputHandler::ProcessNuWroInputFlux(const std::string file) {} FitEvent* NuWroInputHandler::GetNuisanceEvent(const UInt_t entry, const bool lightweight) { // Catch too large entries if (entry >= (UInt_t)fNEvents) return NULL; // Read Entry from TTree to fill NEUT Vect in BaseFitEvt; fNuWroTree->GetEntry(entry); // Run NUISANCE Vector Filler if (!lightweight) { CalcNUISANCEKinematics(); } #ifdef __PROB3PP_ENABLED__ for (size_t i = 0; i < fNUISANCEEvent->fNuwroEvent->in.size(); i++) { if (std::count(PhysConst::pdg_neutrinos, PhysConst::pdg_neutrinos + 4, fNUISANCEEvent->fNuwroEvent->in[i].pdg)) { fNUISANCEEvent->probe_E = fNUISANCEEvent->fNuwroEvent->in[i].t; fNUISANCEEvent->probe_pdg = fNUISANCEEvent->fNuwroEvent->in[i].pdg; break; } } #endif // Setup Input scaling for joint inputs fNUISANCEEvent->InputWeight = GetInputWeight(entry); #ifdef __USE_NUWRO_SRW_EVENTS__ if (!rwEvs.size()) { fNuwroParams = fNuWroEvent->par; } - if (entry >= rwEvs.size()) { rwEvs.push_back(BaseFitEvt()); + rwEvs.back().fType = kNUWRO; + rwEvs.back().Mode = fNUISANCEEvent->Mode; rwEvs.back().fNuwroSRWEvent = SRW::SRWEvent(*fNuWroEvent); rwEvs.back().fNuwroEvent = NULL; rwEvs.back().fNuwroParams = &fNuwroParams; rwEvs.back().probe_E = rwEvs.back().fNuwroSRWEvent.NeutrinoEnergy; rwEvs.back().probe_pdg = rwEvs.back().fNuwroSRWEvent.NeutrinoPDG; } + + fNUISANCEEvent->fNuwroSRWEvent = SRW::SRWEvent(*fNuWroEvent); + fNUISANCEEvent->fNuwroParams = &fNuwroParams; + fNUISANCEEvent->probe_E = fNUISANCEEvent->fNuwroSRWEvent.NeutrinoEnergy; + fNUISANCEEvent->probe_pdg = fNUISANCEEvent->fNuwroSRWEvent.NeutrinoPDG; #endif return fNUISANCEEvent; } int NuWroInputHandler::ConvertNuwroMode(event* e) { Int_t proton_pdg, neutron_pdg, pion_pdg, pion_plus_pdg, pion_minus_pdg, lambda_pdg, eta_pdg, kaon_pdg, kaon_plus_pdg; proton_pdg = 2212; eta_pdg = 221; neutron_pdg = 2112; pion_pdg = 111; pion_plus_pdg = 211; pion_minus_pdg = -211; // O_16_pdg = 100069; // oznacznie z Neuta lambda_pdg = 3122; kaon_pdg = 311; kaon_plus_pdg = 321; if (e->flag.qel) // kwiazielastyczne oddziaływanie { if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem { if (e->flag.cc) return -1; else { if (event1_nof(e, proton_pdg)) return -51; else if (event1_nof(e, neutron_pdg)) return -52; // sprawdzam dodatkowo ? } } else // oddziaływanie z neutrinem { if (e->flag.cc) return 1; else { if (event1_nof(e, proton_pdg)) return 51; else if (event1_nof(e, neutron_pdg)) return 52; } } } if (e->flag.mec) { if (e->flag.anty) return -2; else return 2; } if (e->flag.res) // rezonansowa produkcja: pojedynczy pion, pojed.eta, kaon, // multipiony { Int_t liczba_pionow, liczba_kaonow; liczba_pionow = event1_nof(e, pion_pdg) + event1_nof(e, pion_plus_pdg) + event1_nof(e, pion_minus_pdg); liczba_kaonow = event1_nof(e, kaon_pdg) + event1_nof(e, kaon_pdg); if (liczba_pionow > 1 || liczba_pionow == 0) // multipiony { if (e->flag.anty) { if (e->flag.cc) return -21; else return -41; } else { if (e->flag.cc) return 21; else return 41; } } if (liczba_pionow == 1) { if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem { if (e->flag.cc) { if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_minus_pdg)) return -11; if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_pdg)) return -12; if (event1_nof(e, proton_pdg) && event1_nof(e, pion_minus_pdg)) return -13; } else { if (event1_nof(e, proton_pdg)) { if (event1_nof(e, pion_minus_pdg)) return -33; else if (event1_nof(e, pion_pdg)) return -32; } else if (event1_nof(e, neutron_pdg)) { if (event1_nof(e, pion_plus_pdg)) return -34; else if (event1_nof(e, pion_pdg)) return -31; } } } else // oddziaływanie z neutrinem { if (e->flag.cc) { if (event1_nof(e, proton_pdg) && event1_nof(e, pion_plus_pdg)) return 11; if (event1_nof(e, proton_pdg) && event1_nof(e, pion_pdg)) return 12; if (event1_nof(e, neutron_pdg) && event1_nof(e, pion_plus_pdg)) return 13; } else { if (event1_nof(e, proton_pdg)) { if (event1_nof(e, pion_minus_pdg)) return 33; else if (event1_nof(e, pion_pdg)) return 32; } else if (event1_nof(e, neutron_pdg)) { if (event1_nof(e, pion_plus_pdg)) return 34; else if (event1_nof(e, pion_pdg)) return 31; } } } } if (event1_nof(e, eta_pdg)) // produkcja rezonansowa ety { if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem { if (e->flag.cc) return -22; else { if (event1_nof(e, neutron_pdg)) return -42; else if (event1_nof(e, proton_pdg)) return -43; // sprawdzam dodatkowo ? } } else // oddziaływanie z neutrinem { if (e->flag.cc) return 22; else { if (event1_nof(e, neutron_pdg)) return 42; else if (event1_nof(e, proton_pdg)) return 43; } } } if (event1_nof(e, lambda_pdg) == 1 && liczba_kaonow == 1) // produkcja rezonansowa kaonu { if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem { if (e->flag.cc && event1_nof(e, kaon_pdg)) return -23; else { if (event1_nof(e, kaon_pdg)) return -44; else if (event1_nof(e, kaon_plus_pdg)) return -45; } } else // oddziaływanie z neutrinem { if (e->flag.cc && event1_nof(e, kaon_plus_pdg)) return 23; else { if (event1_nof(e, kaon_pdg)) return 44; else if (event1_nof(e, kaon_plus_pdg)) return 45; } } } } if (e->flag.coh) // koherentne oddziaływanie tylko na O(16) { Int_t _target; _target = e->par.nucleus_p + e->par.nucleus_n; // liczba masowa O(16) if (_target == 16) { if (e->flag.anty) // jeśli jest to oddziaływanie z antyneutrinem { if (e->flag.cc && event1_nof(e, pion_minus_pdg)) return -16; else if (event1_nof(e, pion_pdg)) return -36; } else // oddziaływanie z neutrinem { if (e->flag.cc && event1_nof(e, pion_plus_pdg)) return 16; else if (event1_nof(e, pion_pdg)) return 36; } } } // gleboko nieelastyczne rozpraszanie if (e->flag.dis) { if (e->flag.anty) { if (e->flag.cc) return -26; else return -46; } else { if (e->flag.cc) return 26; else return 46; } } return 9999; } void NuWroInputHandler::CalcNUISANCEKinematics() { // std::cout << "NuWro Event Address " << fNuWroEvent << std::endl; // Reset all variables fNUISANCEEvent->ResetEvent(); FitEvent* evt = fNUISANCEEvent; // Sort Event Info evt->Mode = ConvertNuwroMode(fNuWroEvent); if (abs(evt->Mode) > 60) { evt->Mode = 0; } evt->fEventNo = 0.0; evt->fTotCrs = 0.0; evt->fTargetA = fNuWroEvent->par.nucleus_p + fNuWroEvent->par.nucleus_n; evt->fTargetZ = fNuWroEvent->par.nucleus_p; evt->fTargetH = 0; evt->fBound = (evt->fTargetA) == 1; // Check Particle Stack UInt_t npart_in = fNuWroEvent->in.size(); UInt_t npart_out = fNuWroEvent->out.size(); UInt_t npart_post = fNuWroEvent->post.size(); UInt_t npart = npart_in + npart_out + npart_post; UInt_t kmax = evt->kMaxParticles; if (npart > kmax) { ERR(WRN) << "NUWRO has too many particles. Expanding stack." << std::endl; fNUISANCEEvent->ExpandParticleStack(npart); } // Sort Particles evt->fNParticles = 0; std::vector::iterator p_iter; // Initial State for (p_iter = fNuWroEvent->in.begin(); p_iter != fNuWroEvent->in.end(); p_iter++) { AddNuWroParticle(fNUISANCEEvent, (*p_iter), kInitialState); } // FSI State // for (size_t i = 0; i < npart_in; i++ ) { // AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFSIState); // } // Final State for (p_iter = fNuWroEvent->post.begin(); p_iter != fNuWroEvent->post.end(); p_iter++) { AddNuWroParticle(fNUISANCEEvent, (*p_iter), kFinalState); } // Fill Generator Info if (fSaveExtra) fNuWroInfo->FillGeneratorInfo(fNuWroEvent); // Run Initial, FSI, Final, Other ordering. fNUISANCEEvent->OrderStack(); FitParticle* ISNeutralLepton = fNUISANCEEvent->GetHMISParticle(PhysConst::pdg_neutrinos); if (ISNeutralLepton) { fNUISANCEEvent->probe_E = ISNeutralLepton->E(); fNUISANCEEvent->probe_pdg = ISNeutralLepton->PDG(); } return; } void NuWroInputHandler::AddNuWroParticle(FitEvent* evt, particle& p, int state) { // Add Mom evt->fParticleMom[evt->fNParticles][0] = static_cast(p).x; evt->fParticleMom[evt->fNParticles][1] = static_cast(p).y; evt->fParticleMom[evt->fNParticles][2] = static_cast(p).z; evt->fParticleMom[evt->fNParticles][3] = static_cast(p).t; // Status/PDG evt->fParticleState[evt->fNParticles] = state; evt->fParticlePDG[evt->fNParticles] = p.pdg; // Add to particle count evt->fNParticles++; } void NuWroInputHandler::Print() {} #endif diff --git a/src/InputHandler/StdHepEvt.cxx b/src/InputHandler/StdHepEvt.cxx index dfa93f8..02e821a 100644 --- a/src/InputHandler/StdHepEvt.cxx +++ b/src/InputHandler/StdHepEvt.cxx @@ -1,135 +1,141 @@ // 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 #include #include #include "StdHepEvt.h" // Include logging #include "FitLogger.h" StdHepReader::StdHepReader(){}; bool StdHepReader::SetBranchAddresses(TChain *chain) { bool ok = true; int SBAStatus = 0; SBAStatus = chain->SetBranchAddress("StdHepN", &StdHepN); ok = ok && (SBAStatus || SBAStatus == 5); if (!(!SBAStatus || SBAStatus == 5)) { ERR(WRN) << "Failed to set branch address for \"StdHepN\": " << SBAStatus << std::endl; } SBAStatus = chain->SetBranchAddress("StdHepPdg", StdHepPdg); ok = ok && (SBAStatus || SBAStatus == 5); if (!(!SBAStatus || SBAStatus == 5)) { ERR(WRN) << "Failed to set branch address for \"StdHepPdg\": " << SBAStatus << std::endl; } SBAStatus = chain->SetBranchAddress("StdHepStatus", StdHepStatus); ok = ok && (SBAStatus || SBAStatus == 5); if (!(!SBAStatus || SBAStatus == 5)) { ERR(WRN) << "Failed to set branch address for \"StdHepStatus\": " << SBAStatus << std::endl; } SBAStatus = chain->SetBranchAddress("StdHepP4", StdHepP4); ok = ok && (SBAStatus || SBAStatus == 5); if (!(!SBAStatus || SBAStatus == 5)) { ERR(WRN) << "Failed to set branch address for \"StdHepP4\": " << SBAStatus << std::endl; } return ok; } bool GiBUUStdHepReader::SetBranchAddresses(TChain *chain) { bool ok = true; int SBAStatus = 0; ok = ok && StdHepReader::SetBranchAddresses(chain); SBAStatus = chain->SetBranchAddress("GiBUU2NeutCode", &GiBUU2NeutCode); ok = ok && (SBAStatus || SBAStatus == 5); if (!(!SBAStatus || SBAStatus == 5)) { ERR(WRN) << "Failed to set branch address for \"GiBUU2NeutCode\": " << SBAStatus << std::endl; } + SBAStatus = chain->SetBranchAddress("GiBUUReactionCode", &GiBUUReactionCode); + ok = ok && (SBAStatus || SBAStatus == 5); + if (!(!SBAStatus || SBAStatus == 5)) { + ERR(WRN) << "Failed to set branch address for \"GiBUUReactionCode\": " + << SBAStatus << std::endl; + } SBAStatus = chain->SetBranchAddress("EvtWght", &EvtWght); ok = ok && (SBAStatus || SBAStatus == 5); if (!(!SBAStatus || SBAStatus == 5)) { ERR(WRN) << "Failed to set branch address for \"EvtWght\": " << SBAStatus << std::endl; } return ok; } std::string NegSpacer(double const &num) { return (num >= 0) ? " " : ""; } std::ostream &operator<<(std::ostream &os, TLorentzVector const &tlv) { std::streamsize prec = os.precision(); std::ios_base::fmtflags flags = os.flags(); os.precision(2); os.flags(std::ios::scientific); os << "[" << NegSpacer(tlv[0]) << tlv[0] << "," << NegSpacer(tlv[1]) << tlv[1] << "," << NegSpacer(tlv[2]) << tlv[2] << "," << NegSpacer(tlv[3]) << tlv[3] << ":M(" << tlv.M() << ")]"; os.precision(prec); os.flags(flags); return os; } std::string WriteGiBUUEvent(GiBUUStdHepReader const &gi) { std::stringstream ss(""); ss << "[INFO]: contained " << gi.StdHepN << ", Event Weight: " << std::setprecision(3) << gi.EvtWght << ", NeutConventionReactionCode: " << gi.GiBUU2NeutCode << "\n\t[Lep In](" << std::setw(3) << gi.StdHepPdg[0] << ") " << TLorentzVector(gi.StdHepP4[0][StdHepReader::kStdHepIdxPx], gi.StdHepP4[0][StdHepReader::kStdHepIdxPy], gi.StdHepP4[0][StdHepReader::kStdHepIdxPz], gi.StdHepP4[0][StdHepReader::kStdHepIdxE]) << std::endl; ss << "\t[Target] : " << gi.StdHepPdg[1] << std::endl; ss << "\t[Nuc In] : " << TLorentzVector(gi.StdHepP4[3][StdHepReader::kStdHepIdxPx], gi.StdHepP4[3][StdHepReader::kStdHepIdxPy], gi.StdHepP4[3][StdHepReader::kStdHepIdxPz], gi.StdHepP4[3][StdHepReader::kStdHepIdxE]) << " (" << std::setw(4) << gi.StdHepPdg[3] << ")" << std::endl; for (Int_t stdHepInd = 4; stdHepInd < gi.StdHepN; ++stdHepInd) { ss << "\t[" << std::setw(2) << (stdHepInd - (4)) << "](" << std::setw(5) << gi.StdHepPdg[stdHepInd] << ") " << TLorentzVector(gi.StdHepP4[stdHepInd][StdHepReader::kStdHepIdxPx], gi.StdHepP4[stdHepInd][StdHepReader::kStdHepIdxPy], gi.StdHepP4[stdHepInd][StdHepReader::kStdHepIdxPz], gi.StdHepP4[stdHepInd][StdHepReader::kStdHepIdxE]) << std::endl; } ss << "\t[Lep Out](" << std::setw(3) << gi.StdHepPdg[2] << ") " << TLorentzVector(gi.StdHepP4[2][StdHepReader::kStdHepIdxPx], gi.StdHepP4[2][StdHepReader::kStdHepIdxPy], gi.StdHepP4[2][StdHepReader::kStdHepIdxPz], gi.StdHepP4[2][StdHepReader::kStdHepIdxE]) << std::endl; return ss.str(); } diff --git a/src/InputHandler/StdHepEvt.h b/src/InputHandler/StdHepEvt.h index 198279b..96f280c 100644 --- a/src/InputHandler/StdHepEvt.h +++ b/src/InputHandler/StdHepEvt.h @@ -1,111 +1,116 @@ #ifndef __STDHEPEVT_SEEN__ #define __STDHEPEVT_SEEN__ // 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 "TChain.h" #include "TLorentzVector.h" struct StdHepReader { public: const static int kStdHepIdxPx = 0; const static int kStdHepIdxPy = 1; const static int kStdHepIdxPz = 2; const static int kStdHepIdxE = 3; const static int kStdHepNPmax = 100; StdHepReader(); ///\brief The number of StdHep particles in this event. Int_t StdHepN; ///\brief The PDG codes of particles in this event. /// /// This is determined from the GiBUU particle number by /// GiBUUUtils::GiBUUToPDG. ///\warning This is not a one-to-one mapping, e.g. resonances are not uniquely /// determined by the GiBUU scheme. Int_t StdHepPdg[kStdHepNPmax]; //[StdHepN] ///\brief The StdHep Status of particles in this event. /// /// Status Codes in use: /// - -1: Initial state real particle. /// - 1: Final state real particle. Int_t StdHepStatus[kStdHepNPmax]; //[StdHepN] ///\brief Four momentum for particles in this event. Double_t StdHepP4[kStdHepNPmax][4]; bool SetBranchAddresses(TChain*); }; struct GiBUUStdHepReader : public StdHepReader { GiBUUStdHepReader() : StdHepReader(){}; + ///\brief GiBUU interaction code + /// + /// See https://gibuu.hepforge.org/trac/wiki/LesHouches for details + Int_t GiBUUReactionCode; + ///\brief NEUT equivalent reaction code. /// CC: /// * 1 : QE /// * 2 : 2p2h /// * 10 : Single pion background (non-resonant) /// * 11 : Delta++ ( -11 : Delta- for nubar) /// * 12 : Delta+ (-12 : Delta0 for nubar) /// * 21 : Multi pion production /// * 26 : DIS /// * 4 : Higher resonance, charge: -1 /// * 5 : Higher resonance, charge: 0 /// * 6 : Higher resonance, charge: +1 /// * 7 : Higher resonance, charge: +2 /// /// NC: /// * 30 : Single pion background (non-resonant) /// * 31 : Delta0 /// * 32 : Delta+ /// * 41 : Multi pion production /// * 42 : 2p2h /// * 46 : DIS /// * 47 : Higher resonance, charge: -1 /// * 48 : Higher resonance, charge: 0 /// * 49 : Higher resonance, charge: +1 /// * 50 : Higher resonance, charge: +2 /// * 51 : NCEL proton-target /// * 52 : NCEL neutron-target /// Int_t GiBUU2NeutCode; ///\brief The total XSec weighting that should be applied to this event. Double_t EvtWght; ///\brief Weighting which takes account of multiple input numu species. /// /// Defined such that W_numu + W_numubar = 1 Double_t SpeciesWght_numu; ///\brief Weighting which takes account of multiple input nue species. /// /// Defined such that W_nue + W_nuebar = 1 Double_t SpeciesWght_nue; ///\brief Weighting which takes account of multiple input neutrino species. /// /// Defined such that \Sum_species W_species = 1 Double_t SpeciesWght; bool SetBranchAddresses(TChain*); }; std::string WriteGiBUUEvent(GiBUUStdHepReader const& gi); #endif diff --git a/src/Logger/Initialiser.cxx b/src/Logger/Initialiser.cxx index 923805a..8cb09f2 100644 --- a/src/Logger/Initialiser.cxx +++ b/src/Logger/Initialiser.cxx @@ -1,98 +1,101 @@ #include "Initialiser.h" -void RunNuisance(){ - std::cout << "Starting NUISANCE" << std::endl; -} +void RunNuisance() { std::cout << "Starting NUISANCE" << std::endl; } struct LetterBackronym { LetterBackronym(size_t n, std::string const &b, float p = 1.0, std::string const &t = "") { NUsed = n; Backkie = b; ProbAccept = p; TagLine = t; }; size_t NUsed; float ProbAccept; std::string Backkie; std::string TagLine; }; -__attribute__((constructor)) void nuisance_init(void) { +__attribute__((constructor)) void constructor(void) { nuisance_init(); } + +void nuisance_init(void) { std::vector > Letters; for (size_t i = 0; i < 8; ++i) { Letters.push_back(std::vector()); } Letters[0].push_back(LetterBackronym(2, "Neutrino")); Letters[0].push_back(LetterBackronym(3, "NUIsance", 0.2)); Letters[2].push_back(LetterBackronym(1, "Interaction")); Letters[3].push_back(LetterBackronym(1, "Systematics")); Letters[3].push_back(LetterBackronym( - 1, "Synthesiser", 0.2, "Playing on the comparisons you want to see")); + 1, "Synthesiser", 0.2, "Playing on the comparisons you want to see")); Letters[4].push_back(LetterBackronym(2, "ANalyser")); Letters[4].push_back(LetterBackronym(1, "Aggregating", 0.5)); Letters[4].push_back(LetterBackronym(3, "from A-Neutrino sCattering", 1, "You can always find a frame")); Letters[5].push_back( - LetterBackronym(1, "New", 1, "The freshest comparisons")); + LetterBackronym(1, "New", 1, "The freshest comparisons")); Letters[6].push_back(LetterBackronym(1, "by Comparing")); Letters[6].push_back(LetterBackronym(1, "Constraints from")); Letters[7].push_back(LetterBackronym(1, "Experiments")); std::vector TagLines; TagLines.push_back("Fit and compare."); std::stringstream back(""); TRandom3 tr; tr.SetSeed(); for (size_t i = 0; i < 8;) { LetterBackronym const &let = Letters[i][tr.Integer(Letters[i].size())]; if (tr.Uniform() > let.ProbAccept) { continue; } back << let.Backkie << " "; i += let.NUsed; if (let.TagLine.length()) { TagLines.push_back(let.TagLine); } } std::string Name = "Nuisance"; std::string TagL = TagLines[tr.Integer(TagLines.size())]; std::vector > > - OneBlob; + OneBlob; OneBlob.push_back( - std::make_pair("NUISANCE", std::make_pair("", "FiXing your Neutrinos"))); + std::make_pair("NUISANCE", std::make_pair("", "FiXing your Neutrinos"))); if (tr.Uniform() < 0.01) { std::pair > const &blob = - OneBlob[tr.Integer(OneBlob.size())]; + OneBlob[tr.Integer(OneBlob.size())]; Name = blob.first; back.str(""); back << blob.second.first; TagL = blob.second.second; } - std::cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" - "%%%%%%%%%%%%%%%" - "%%" - << std::endl - << "%% Welcome to " << Name << ": \033[5m" << back.str() - << "\033[0m-- " << TagL << std::endl - << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" - "%%%%%%%%%%%%%%%" - "%%" - << std::endl; -} \ No newline at end of file + // std::cout << + // "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + // "%%%%%%%%%%%%%%%" + // "%%" + // << std::endl + // << "%% Welcome to " << Name << ": \033[5m" << back.str() + // << "\033[0m-- " << TagL << std::endl + // << + // "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" + // "%%%%%%%%%%%%%%%" + // "%%" + // << std::endl; + std::cout << Name << ": " << back.str() << " -- " << TagL << std::endl; +} diff --git a/src/Logger/Initialiser.h b/src/Logger/Initialiser.h index 104ead9..4386235 100644 --- a/src/Logger/Initialiser.h +++ b/src/Logger/Initialiser.h @@ -1,20 +1,21 @@ #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include "TRandom3.h" void RunNuisance(); +void nuisance_init(void) ; diff --git a/src/Reweight/NUISANCEWeightCalcs.cxx b/src/Reweight/NUISANCEWeightCalcs.cxx index d21637b..b3d6d09 100644 --- a/src/Reweight/NUISANCEWeightCalcs.cxx +++ b/src/Reweight/NUISANCEWeightCalcs.cxx @@ -1,394 +1,391 @@ #include "NUISANCEWeightCalcs.h" -#include "GeneralUtils.h" #include "FitEvent.h" -#include "WeightUtils.h" +#include "GeneralUtils.h" #include "NUISANCESyst.h" +#include "WeightUtils.h" using namespace Reweight; -ModeNormCalc::ModeNormCalc(){ - fNormRES = 1.0; -} +ModeNormCalc::ModeNormCalc() { fNormRES = 1.0; } double ModeNormCalc::CalcWeight(BaseFitEvt* evt) { int mode = abs(evt->Mode); double w = 1.0; - if (mode == 11 or mode == 12 or mode == 13){ + if (mode == 11 or mode == 12 or mode == 13) { w *= fNormRES; } return w; } void ModeNormCalc::SetDialValue(std::string name, double val) { SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); } void ModeNormCalc::SetDialValue(int rwenum, double val) { - int curenum = rwenum % 1000; // Check Handled if (!IsHandled(curenum)) return; if (curenum == kModeNorm_NormRES) fNormRES = val; } bool ModeNormCalc::IsHandled(int rwenum) { - int curenum = rwenum % 1000; switch (curenum) { - case kModeNorm_NormRES: return true; - default: return false; + case kModeNorm_NormRES: + return true; + default: + return false; } } - - -SBLOscWeightCalc::SBLOscWeightCalc(){ +SBLOscWeightCalc::SBLOscWeightCalc() { fDistance = 0.0; fMassSplitting = 0.0; fSin2Theta = 0.0; } double SBLOscWeightCalc::CalcWeight(BaseFitEvt* evt) { FitEvent* fevt = static_cast(evt); FitParticle* pnu = fevt->PartInfo(0); double E = pnu->fP.E() / 1.E3; // Extract energy return GetSBLOscWeight(E); } void SBLOscWeightCalc::SetDialValue(std::string name, double val) { SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); } void SBLOscWeightCalc::SetDialValue(int rwenum, double val) { - int curenum = rwenum % 1000; if (!IsHandled(curenum)) return; - if (curenum == kSBLOsc_Distance) fDistance = val; + if (curenum == kSBLOsc_Distance) fDistance = val; if (curenum == kSBLOsc_MassSplitting) fMassSplitting = val; - if (curenum == kSBLOsc_Sin2Theta) fSin2Theta = val; + if (curenum == kSBLOsc_Sin2Theta) fSin2Theta = val; } bool SBLOscWeightCalc::IsHandled(int rwenum) { - int curenum = rwenum % 1000; switch (curenum) { - case kSBLOsc_Distance: return true; - case kSBLOsc_MassSplitting: return true; - case kSBLOsc_Sin2Theta: return true; - default: return false; + case kSBLOsc_Distance: + return true; + case kSBLOsc_MassSplitting: + return true; + case kSBLOsc_Sin2Theta: + return true; + default: + return false; } } -double SBLOscWeightCalc::GetSBLOscWeight(double E){ +double SBLOscWeightCalc::GetSBLOscWeight(double E) { if (E <= 0.0) return 1.0 - 0.5 * fSin2Theta; - return 1.0 - fSin2Theta * pow( sin(1.267 * fMassSplitting * fDistance / E ), 2); + return 1.0 - fSin2Theta * pow(sin(1.267 * fMassSplitting * fDistance / E), 2); } - GaussianModeCorr::GaussianModeCorr() { - - // Init - fApply_CCQE = false; - fGausVal_CCQE[kPosNorm] = 0.0; - fGausVal_CCQE[kPosTilt] = 0.0; - fGausVal_CCQE[kPosPq0] = 1.0; - fGausVal_CCQE[kPosWq0] = 1.0; - fGausVal_CCQE[kPosPq3] = 1.0; - fGausVal_CCQE[kPosWq3] = 1.0; - - fApply_2p2h = false; - fGausVal_2p2h[kPosNorm] = 0.0; - fGausVal_2p2h[kPosTilt] = 0.0; - fGausVal_2p2h[kPosPq0] = 1.0; - fGausVal_2p2h[kPosWq0] = 1.0; - fGausVal_2p2h[kPosPq3] = 1.0; - fGausVal_2p2h[kPosWq3] = 1.0; - - fApply_2p2h_PPandNN = false; - fGausVal_2p2h_PPandNN[kPosNorm] = 0.0; - fGausVal_2p2h_PPandNN[kPosTilt] = 0.0; - fGausVal_2p2h_PPandNN[kPosPq0] = 1.0; - fGausVal_2p2h_PPandNN[kPosWq0] = 1.0; - fGausVal_2p2h_PPandNN[kPosPq3] = 1.0; - fGausVal_2p2h_PPandNN[kPosWq3] = 1.0; - - fApply_2p2h_NP = false; - fGausVal_2p2h_NP[kPosNorm] = 0.0; - fGausVal_2p2h_NP[kPosTilt] = 0.0; - fGausVal_2p2h_NP[kPosPq0] = 1.0; - fGausVal_2p2h_NP[kPosWq0] = 1.0; - fGausVal_2p2h_NP[kPosPq3] = 1.0; - fGausVal_2p2h_NP[kPosWq3] = 1.0; - - fApply_CC1pi = false; - fGausVal_CC1pi[kPosNorm] = 0.0; - fGausVal_CC1pi[kPosTilt] = 0.0; - fGausVal_CC1pi[kPosPq0] = 1.0; - fGausVal_CC1pi[kPosWq0] = 1.0; - fGausVal_CC1pi[kPosPq3] = 1.0; - fGausVal_CC1pi[kPosWq3] = 1.0; - - fAllowSuppression = false; - - fDebugStatements = FitPar::Config().GetParB("GaussianModeCorr_DEBUG"); - + // Init + fApply_CCQE = false; + fGausVal_CCQE[kPosNorm] = 0.0; + fGausVal_CCQE[kPosTilt] = 0.0; + fGausVal_CCQE[kPosPq0] = 1.0; + fGausVal_CCQE[kPosWq0] = 1.0; + fGausVal_CCQE[kPosPq3] = 1.0; + fGausVal_CCQE[kPosWq3] = 1.0; + + fApply_2p2h = false; + fGausVal_2p2h[kPosNorm] = 0.0; + fGausVal_2p2h[kPosTilt] = 0.0; + fGausVal_2p2h[kPosPq0] = 1.0; + fGausVal_2p2h[kPosWq0] = 1.0; + fGausVal_2p2h[kPosPq3] = 1.0; + fGausVal_2p2h[kPosWq3] = 1.0; + + fApply_2p2h_PPandNN = false; + fGausVal_2p2h_PPandNN[kPosNorm] = 0.0; + fGausVal_2p2h_PPandNN[kPosTilt] = 0.0; + fGausVal_2p2h_PPandNN[kPosPq0] = 1.0; + fGausVal_2p2h_PPandNN[kPosWq0] = 1.0; + fGausVal_2p2h_PPandNN[kPosPq3] = 1.0; + fGausVal_2p2h_PPandNN[kPosWq3] = 1.0; + + fApply_2p2h_NP = false; + fGausVal_2p2h_NP[kPosNorm] = 0.0; + fGausVal_2p2h_NP[kPosTilt] = 0.0; + fGausVal_2p2h_NP[kPosPq0] = 1.0; + fGausVal_2p2h_NP[kPosWq0] = 1.0; + fGausVal_2p2h_NP[kPosPq3] = 1.0; + fGausVal_2p2h_NP[kPosWq3] = 1.0; + + fApply_CC1pi = false; + fGausVal_CC1pi[kPosNorm] = 0.0; + fGausVal_CC1pi[kPosTilt] = 0.0; + fGausVal_CC1pi[kPosPq0] = 1.0; + fGausVal_CC1pi[kPosWq0] = 1.0; + fGausVal_CC1pi[kPosPq3] = 1.0; + fGausVal_CC1pi[kPosWq3] = 1.0; + + fAllowSuppression = false; + + fDebugStatements = FitPar::Config().GetParB("GaussianModeCorr_DEBUG"); } double GaussianModeCorr::CalcWeight(BaseFitEvt* evt) { + FitEvent* fevt = static_cast(evt); + double rw_weight = 1.0; - FitEvent* fevt = static_cast(evt); - double rw_weight = 1.0; - - // Get Neutrino - if (!fevt->Npart()){ - THROW("NO particles found in stack!"); - } - FitParticle* pnu = fevt->PartInfo(0); - if (!pnu){ - THROW("NO Starting particle found in stack!"); - } - int pdgnu = pnu->fPID; - - FitParticle* plep = fevt->GetHMFSParticle(abs(pdgnu) - 1); - if (!plep) return 1.0; - - TLorentzVector q = pnu->fP - plep->fP; + // Get Neutrino + if (!fevt->Npart()) { + THROW("NO particles found in stack!"); + } + FitParticle* pnu = fevt->GetHMISAnyLeptons(); + if (!pnu) { + THROW("NO Starting particle found in stack!"); + } + int pdgnu = pnu->fPID; - // Extra q0,q3 - double q0 = fabs(q.E()) / 1.E3; - double q3 = fabs(q.Vect().Mag()) / 1.E3; + int expect_fsleppdg = 0; - int initialstate = -1; // Undef - if (abs(fevt->Mode) == 2) { + if (pdgnu & 1) { + expect_fsleppdg = pdgnu; + } else { + expect_fsleppdg = abs(pdgnu) - 1; + } - int npr = 0; - int nne = 0; + FitParticle* plep = fevt->GetHMFSParticle(expect_fsleppdg); + if (!plep) return 1.0; - for (UInt_t j = 0; j < fevt->Npart(); j++) { - if ((fevt->PartInfo(j))->fIsAlive) continue; + TLorentzVector q = pnu->fP - plep->fP; - if (fevt->PartInfo(j)->fPID == 2212) npr++; - else if (fevt->PartInfo(j)->fPID == 2112) nne++; - } - // std::cout << "PN State = " << npr << " " << nne << std::endl; + // Extra q0,q3 + double q0 = fabs(q.E()) / 1.E3; + double q3 = fabs(q.Vect().Mag()) / 1.E3; - if (fevt->Mode == 2 and npr == 1 and nne == 1) { - initialstate = 2; + int initialstate = -1; // Undef + if (abs(fevt->Mode) == 2) { + int npr = 0; + int nne = 0; - } else if (fevt->Mode == 2 and ((npr == 0 and nne == 2) or (npr == 2 and nne == 0))) { - initialstate = 1; - } - } + for (UInt_t j = 0; j < fevt->Npart(); j++) { + if ((fevt->PartInfo(j))->fIsAlive) continue; -// std::cout << "Got q0 q3 = " << q0 << " " << q3 << std::endl; + if (fevt->PartInfo(j)->fPID == 2212) + npr++; + else if (fevt->PartInfo(j)->fPID == 2112) + nne++; + } + // std::cout << "PN State = " << npr << " " << nne << std::endl; -// Apply weighting - if (fApply_CCQE and abs(fevt->Mode) == 1) { - if (fDebugStatements) std::cout << "Getting CCQE Weight" << std::endl; - double g = GetGausWeight(q0, q3, fGausVal_CCQE); - if (g < 1.0) g = 1.0; - rw_weight *= g; - } + if (fevt->Mode == 2 and npr == 1 and nne == 1) { + initialstate = 2; - if (fApply_2p2h and abs(fevt->Mode) == 2) { - if (fDebugStatements) std::cout << "Getting 2p2h Weight" << std::endl; - rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h); - } + } else if (fevt->Mode == 2 and + ((npr == 0 and nne == 2) or (npr == 2 and nne == 0))) { + initialstate = 1; + } + } - if (fApply_2p2h_PPandNN and abs(fevt->Mode) == 2 and initialstate == 1) { - if (fDebugStatements) std::cout << "Getting 2p2h PPandNN Weight" << std::endl; - rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_PPandNN); - } + // std::cout << "Got q0 q3 = " << q0 << " " << q3 << std::endl; - if (fApply_2p2h_NP and abs(fevt->Mode) == 2 and initialstate == 2) { - if (fDebugStatements) std::cout << "Getting 2p2h NP Weight" << std::endl; - rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_NP); - } + // Apply weighting + if (fApply_CCQE and abs(fevt->Mode) == 1) { + if (fDebugStatements) std::cout << "Getting CCQE Weight" << std::endl; + double g = GetGausWeight(q0, q3, fGausVal_CCQE); + if (g < 1.0) g = 1.0; + rw_weight *= g; + } - if (fApply_CC1pi and abs(fevt->Mode) >= 11 and abs(fevt->Mode) <= 13) { - if (fDebugStatements) std::cout << "Getting CC1pi Weight" << std::endl; - rw_weight *= GetGausWeight(q0, q3, fGausVal_CC1pi); - } + if (fApply_2p2h and abs(fevt->Mode) == 2) { + if (fDebugStatements) std::cout << "Getting 2p2h Weight" << std::endl; + rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h); + } + if (fApply_2p2h_PPandNN and abs(fevt->Mode) == 2 and initialstate == 1) { + if (fDebugStatements) + std::cout << "Getting 2p2h PPandNN Weight" << std::endl; + rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_PPandNN); + } + if (fApply_2p2h_NP and abs(fevt->Mode) == 2 and initialstate == 2) { + if (fDebugStatements) std::cout << "Getting 2p2h NP Weight" << std::endl; + rw_weight *= GetGausWeight(q0, q3, fGausVal_2p2h_NP); + } + if (fApply_CC1pi and abs(fevt->Mode) >= 11 and abs(fevt->Mode) <= 13) { + if (fDebugStatements) std::cout << "Getting CC1pi Weight" << std::endl; + rw_weight *= GetGausWeight(q0, q3, fGausVal_CC1pi); + } - // if (fDebugStatements) std::cout << "Returning Weight " << rw_weight << std::endl; - return rw_weight; + // if (fDebugStatements) std::cout << "Returning Weight " << rw_weight << + // std::endl; + return rw_weight; } double GaussianModeCorr::GetGausWeight(double q0, double q3, double vals[]) { + // // CCQE Without Suppression + // double Norm = 4.82788679036; + // double Tilt = 2.3501416116; + // double Pq0 = 0.363964889702; + // double Wq0 = 0.133976806938; + // double Pq3 = 0.431769740224; + // double Wq3 = 0.207666663434; + + // // Also add for CCQE at the end + // return (w > 1.0) ? w : 1.0; + + // // 2p2h with suppression + // double Norm = 15.967; + // double Tilt = -0.455655; + // double Pq0 = 0.214598; + // double Wq0 = 0.0291061; + // double Pq3 = 0.480194; + // double Wq3 = 0.134588; + + double Norm = vals[kPosNorm]; + double Tilt = vals[kPosTilt]; + double Pq0 = vals[kPosPq0]; + double Wq0 = vals[kPosWq0]; + double Pq3 = vals[kPosPq3]; + double Wq3 = vals[kPosWq3]; + + double a = cos(Tilt) * cos(Tilt) / (2 * Wq0 * Wq0); + a += sin(Tilt) * sin(Tilt) / (2 * Wq3 * Wq3); + + double b = -sin(2 * Tilt) / (4 * Wq0 * Wq0); + b += sin(2 * Tilt) / (4 * Wq3 * Wq3); + + double c = sin(Tilt) * sin(Tilt) / (2 * Wq0 * Wq0); + c += cos(Tilt) * cos(Tilt) / (2 * Wq3 * Wq3); + + double w = Norm; + w *= exp(-a * (q0 - Pq0) * (q0 - Pq0)); + w *= exp(+2.0 * b * (q0 - Pq0) * (q3 - Pq3)); + w *= exp(-c * (q3 - Pq3) * (q3 - Pq3)); + + if (fDebugStatements) { + std::cout << "Applied Tilt " << Tilt << " " << cos(Tilt) << " " << sin(Tilt) + << std::endl; + std::cout << "abc = " << a << " " << b << " " << c << std::endl; + std::cout << "Returning " << Norm << " " << Pq0 << " " << Wq0 << " " << Pq3 + << " " << Wq3 << " " << w << std::endl; + } - // // CCQE Without Suppression - // double Norm = 4.82788679036; - // double Tilt = 2.3501416116; - // double Pq0 = 0.363964889702; - // double Wq0 = 0.133976806938; - // double Pq3 = 0.431769740224; - // double Wq3 = 0.207666663434; - - // // Also add for CCQE at the end - // return (w > 1.0) ? w : 1.0; - - // // 2p2h with suppression - // double Norm = 15.967; - // double Tilt = -0.455655; - // double Pq0 = 0.214598; - // double Wq0 = 0.0291061; - // double Pq3 = 0.480194; - // double Wq3 = 0.134588; - - - double Norm = vals[kPosNorm]; - double Tilt = vals[kPosTilt]; - double Pq0 = vals[kPosPq0]; - double Wq0 = vals[kPosWq0]; - double Pq3 = vals[kPosPq3]; - double Wq3 = vals[kPosWq3]; - - double a = cos(Tilt) * cos(Tilt) / (2 * Wq0 * Wq0); - a += sin(Tilt) * sin(Tilt) / (2 * Wq3 * Wq3); - - double b = - sin(2 * Tilt) / (4 * Wq0 * Wq0); - b += sin(2 * Tilt) / (4 * Wq3 * Wq3); - - double c = sin(Tilt) * sin(Tilt) / (2 * Wq0 * Wq0); - c += cos(Tilt) * cos(Tilt) / (2 * Wq3 * Wq3); - - double w = Norm; - w *= exp(-a * (q0 - Pq0) * (q0 - Pq0)); - w *= exp(+2.0 * b * (q0 - Pq0) * (q3 - Pq3)); - w *= exp(-c * (q3 - Pq3) * (q3 - Pq3)); - - if (fDebugStatements) { + if (w != w || isnan(w) || w < 0.0) { + w = 0.0; + } - std::cout << "Applied Tilt " << Tilt << " " << cos(Tilt) << " " << sin(Tilt) << std::endl; - std::cout << "abc = " << a << " " << b << " " << c << std::endl; - std::cout << "Returning " << Norm << " " << Pq0 << " " << Wq0 << " " << Pq3 << " " << Wq3 << " " << w << std::endl; + if (w < 1.0 and !fAllowSuppression) { + w = 1.0; + } - } + return w; +} +void GaussianModeCorr::SetDialValue(std::string name, double val) { + SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); +} - if (w != w || isnan(w) || w < 0.0){ - w = 0.0; - } +void GaussianModeCorr::SetDialValue(int rwenum, double val) { + int curenum = rwenum % 1000; - if (w < 1.0 and !fAllowSuppression){ - w = 1.0; - } + // Check Handled + if (!IsHandled(curenum)) return; - return w; -} + // CCQE Setting + for (int i = kGaussianCorr_CCQE_norm; i <= kGaussianCorr_CCQE_Wq3; i++) { + if (i == curenum) { + int index = i - kGaussianCorr_CCQE_norm; + fGausVal_CCQE[index] = val; + fApply_CCQE = true; + } + } + // 2p2h Setting + for (int i = kGaussianCorr_2p2h_norm; i <= kGaussianCorr_2p2h_Wq3; i++) { + if (i == curenum) { + int index = i - kGaussianCorr_2p2h_norm; + fGausVal_2p2h[index] = val; + fApply_2p2h = true; + } + } -void GaussianModeCorr::SetDialValue(std::string name, double val) { - SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); -} + // 2p2h_PPandNN Setting + for (int i = kGaussianCorr_2p2h_PPandNN_norm; + i <= kGaussianCorr_2p2h_PPandNN_Wq3; i++) { + if (i == curenum) { + int index = i - kGaussianCorr_2p2h_PPandNN_norm; + fGausVal_2p2h_PPandNN[index] = val; + fApply_2p2h_PPandNN = true; + } + } -void GaussianModeCorr::SetDialValue(int rwenum, double val) { + // 2p2h_NP Setting + for (int i = kGaussianCorr_2p2h_NP_norm; i <= kGaussianCorr_2p2h_NP_Wq3; + i++) { + if (i == curenum) { + int index = i - kGaussianCorr_2p2h_NP_norm; + fGausVal_2p2h_NP[index] = val; + fApply_2p2h_NP = true; + } + } - int curenum = rwenum % 1000; - - // Check Handled - if (!IsHandled(curenum)) return; - - - // CCQE Setting - for (int i = kGaussianCorr_CCQE_norm; i <= kGaussianCorr_CCQE_Wq3; i++) { - if (i == curenum) { - int index = i - kGaussianCorr_CCQE_norm; - fGausVal_CCQE[index] = val; - fApply_CCQE = true; - } - } - - // 2p2h Setting - for (int i = kGaussianCorr_2p2h_norm; i <= kGaussianCorr_2p2h_Wq3; i++) { - if (i == curenum) { - int index = i - kGaussianCorr_2p2h_norm; - fGausVal_2p2h[index] = val; - fApply_2p2h = true; - } - } - - // 2p2h_PPandNN Setting - for (int i = kGaussianCorr_2p2h_PPandNN_norm; i <= kGaussianCorr_2p2h_PPandNN_Wq3; i++) { - if (i == curenum) { - int index = i - kGaussianCorr_2p2h_PPandNN_norm; - fGausVal_2p2h_PPandNN[index] = val; - fApply_2p2h_PPandNN = true; - } - } - - // 2p2h_NP Setting - for (int i = kGaussianCorr_2p2h_NP_norm; i <= kGaussianCorr_2p2h_NP_Wq3; i++) { - if (i == curenum) { - int index = i - kGaussianCorr_2p2h_NP_norm; - fGausVal_2p2h_NP[index] = val; - fApply_2p2h_NP = true; - } - } - - // CC1pi Setting - for (int i = kGaussianCorr_CC1pi_norm; i <= kGaussianCorr_CC1pi_Wq3; i++) { - if (i == curenum) { - int index = i - kGaussianCorr_CC1pi_norm; - fGausVal_CC1pi[index] = val; - fApply_CC1pi = true; - } - } - - if (curenum == kGaussianCorr_AllowSuppression){ - fAllowSuppression = (val > 0.5); - } + // CC1pi Setting + for (int i = kGaussianCorr_CC1pi_norm; i <= kGaussianCorr_CC1pi_Wq3; i++) { + if (i == curenum) { + int index = i - kGaussianCorr_CC1pi_norm; + fGausVal_CC1pi[index] = val; + fApply_CC1pi = true; + } + } + if (curenum == kGaussianCorr_AllowSuppression) { + fAllowSuppression = (val > 0.5); + } } bool GaussianModeCorr::IsHandled(int rwenum) { - - int curenum = rwenum % 1000; - switch (curenum) { - case kGaussianCorr_CCQE_norm: - case kGaussianCorr_CCQE_tilt: - case kGaussianCorr_CCQE_Pq0: - case kGaussianCorr_CCQE_Wq0: - case kGaussianCorr_CCQE_Pq3: - case kGaussianCorr_CCQE_Wq3: - case kGaussianCorr_2p2h_norm: - case kGaussianCorr_2p2h_tilt: - case kGaussianCorr_2p2h_Pq0: - case kGaussianCorr_2p2h_Wq0: - case kGaussianCorr_2p2h_Pq3: - case kGaussianCorr_2p2h_Wq3: - case kGaussianCorr_2p2h_PPandNN_norm: - case kGaussianCorr_2p2h_PPandNN_tilt: - case kGaussianCorr_2p2h_PPandNN_Pq0: - case kGaussianCorr_2p2h_PPandNN_Wq0: - case kGaussianCorr_2p2h_PPandNN_Pq3: - case kGaussianCorr_2p2h_PPandNN_Wq3: - case kGaussianCorr_2p2h_NP_norm: - case kGaussianCorr_2p2h_NP_tilt: - case kGaussianCorr_2p2h_NP_Pq0: - case kGaussianCorr_2p2h_NP_Wq0: - case kGaussianCorr_2p2h_NP_Pq3: - case kGaussianCorr_2p2h_NP_Wq3: - case kGaussianCorr_CC1pi_norm: - case kGaussianCorr_CC1pi_tilt: - case kGaussianCorr_CC1pi_Pq0: - case kGaussianCorr_CC1pi_Wq0: - case kGaussianCorr_CC1pi_Pq3: - case kGaussianCorr_CC1pi_Wq3: - case kGaussianCorr_AllowSuppression: - - return true; - default: - return false; - } + int curenum = rwenum % 1000; + switch (curenum) { + case kGaussianCorr_CCQE_norm: + case kGaussianCorr_CCQE_tilt: + case kGaussianCorr_CCQE_Pq0: + case kGaussianCorr_CCQE_Wq0: + case kGaussianCorr_CCQE_Pq3: + case kGaussianCorr_CCQE_Wq3: + case kGaussianCorr_2p2h_norm: + case kGaussianCorr_2p2h_tilt: + case kGaussianCorr_2p2h_Pq0: + case kGaussianCorr_2p2h_Wq0: + case kGaussianCorr_2p2h_Pq3: + case kGaussianCorr_2p2h_Wq3: + case kGaussianCorr_2p2h_PPandNN_norm: + case kGaussianCorr_2p2h_PPandNN_tilt: + case kGaussianCorr_2p2h_PPandNN_Pq0: + case kGaussianCorr_2p2h_PPandNN_Wq0: + case kGaussianCorr_2p2h_PPandNN_Pq3: + case kGaussianCorr_2p2h_PPandNN_Wq3: + case kGaussianCorr_2p2h_NP_norm: + case kGaussianCorr_2p2h_NP_tilt: + case kGaussianCorr_2p2h_NP_Pq0: + case kGaussianCorr_2p2h_NP_Wq0: + case kGaussianCorr_2p2h_NP_Pq3: + case kGaussianCorr_2p2h_NP_Wq3: + case kGaussianCorr_CC1pi_norm: + case kGaussianCorr_CC1pi_tilt: + case kGaussianCorr_CC1pi_Pq0: + case kGaussianCorr_CC1pi_Wq0: + case kGaussianCorr_CC1pi_Pq3: + case kGaussianCorr_CC1pi_Wq3: + case kGaussianCorr_AllowSuppression: + + return true; + default: + return false; + } } diff --git a/src/Reweight/OscWeightEngine.cxx b/src/Reweight/OscWeightEngine.cxx index bcb9292..1577175 100644 --- a/src/Reweight/OscWeightEngine.cxx +++ b/src/Reweight/OscWeightEngine.cxx @@ -1,329 +1,331 @@ // 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 . *******************************************************************************/ //#define DEBUG_OSC_WE #include "OscWeightEngine.h" #include enum nuTypes { kNuebarType = -1, kNumubarType = -2, kNutaubarType = -3, kNueType = 1, kNumuType = 2, kNutauType = 3, }; nuTypes GetNuType(int pdg) { switch (pdg) { case 16: return kNutauType; case 14: return kNumuType; case 12: return kNueType; case -16: return kNutaubarType; case -14: return kNumubarType; case -12: return kNuebarType; default: { THROW("Attempting to convert \"neutrino pdg\": " << pdg); } } } OscWeightEngine::OscWeightEngine() : #ifdef __PROB3PP_ENABLED__ bp(), #endif theta12(0.825), theta13(0.10), theta23(1.0), dm12(7.9e-5), dm23(2.5e-3), dcp(0.0), LengthParam(0xdeadbeef), TargetNuType(0), ForceFromNuPDG(0) { Config(); } void OscWeightEngine::Config() { std::vector OscParam = Config::QueryKeys("OscParam"); if (OscParam.size() < 1) { ERROR(WRN, "Oscillation parameters specified but no OscParam element " "configuring the experimental characteristics found.\nExpect at " "least . Pausing for " "10..."); sleep(10); return; } if (OscParam[0].Has("baseline_km")) { LengthParamIsZenith = false; LengthParam = OscParam[0].GetD("baseline_km"); constant_density = OscParam[0].Has("matter_density") ? OscParam[0].GetD("matter_density") : 0xdeadbeef; } else if (OscParam[0].Has("detection_zenith_deg")) { LengthParamIsZenith = true; static const double deg2rad = asin(1) / 90.0; LengthParam = cos(OscParam[0].GetD("detection_zenith_deg") * deg2rad); } else { ERROR(WRN, "It appeared that you wanted to set up an oscillation weight " "branch, but it was not correctly configured. You need to specify " "either: detection_zenith_deg or baseline_km attributes on the " "OscParam element, and if baseline_km is specified, you can " "optionally also set matter_density for oscillations through a " "constant matter density. Pausing for 10..."); sleep(10); return; } dm23 = OscParam[0].Has("dm23") ? OscParam[0].GetD("dm23") : dm23; theta23 = OscParam[0].Has("sinsq_theta23") ? OscParam[0].GetD("sinsq_theta23") : theta23; theta13 = OscParam[0].Has("sinsq_theta13") ? OscParam[0].GetD("sinsq_theta13") : theta13; dm12 = OscParam[0].Has("dm12") ? OscParam[0].GetD("dm12") : dm12; theta12 = OscParam[0].Has("sinsq_theta12") ? OscParam[0].GetD("sinsq_theta12") : theta12; dcp = OscParam[0].Has("dcp") ? OscParam[0].GetD("dcp") : dcp; TargetNuType = OscParam[0].Has("TargetNuPDG") ? GetNuType(OscParam[0].GetI("TargetNuPDG")) : 0; ForceFromNuPDG = OscParam[0].Has("ForceFromNuPDG") ? GetNuType(OscParam[0].GetI("ForceFromNuPDG")) : 0; QLOG(FIT, "Configured oscillation weighter:"); if (LengthParamIsZenith) { QLOG(FIT, "Earth density profile with detection cos(zenith) = " << LengthParam); } else { if (constant_density != 0xdeadbeef) { QLOG(FIT, "Constant density with experimental baseline = " << LengthParam); } else { QLOG(FIT, "Vacuum oscillations with experimental baseline = " << LengthParam); } } params[0] = dm23; params[1] = theta23; params[2] = theta13; params[3] = dm12; params[4] = theta12; params[5] = dcp; QLOG(FIT, "\tdm23 : " << params[0]); QLOG(FIT, "\tsinsq_theta23: " << params[1]); QLOG(FIT, "\tsinsq_theta13: " << params[2]); QLOG(FIT, "\tdm12 : " << params[3]); QLOG(FIT, "\tsinsq_theta12: " << params[4]); QLOG(FIT, "\tdcp : " << params[5]); if (TargetNuType) { QLOG(FIT, "\tTargetNuType: " << TargetNuType); } if (ForceFromNuPDG) { QLOG(FIT, "\tForceFromNuPDG: " << ForceFromNuPDG); } #ifdef __PROB3PP_ENABLED__ bp.SetMNS(params[theta12_idx], params[theta13_idx], params[theta23_idx], params[dm12_idx], params[dm23_idx], params[dcp_idx], 1, true, 2); bp.DefinePath(LengthParam, 0); - QLOG(FIT, "\tBaseline : " << (bp.GetBaseline() / 100.0) << " km."); + if (LengthParamIsZenith) { + QLOG(FIT, "\tBaseline : " << (bp.GetBaseline() / 100.0) << " km."); + } #endif } void OscWeightEngine::IncludeDial(std::string name, double startval) { #ifdef DEBUG_OSC_WE std::cout << "IncludeDial: " << name << " at " << startval << std::endl; #endif int dial = SystEnumFromString(name); if (!dial) { THROW("OscWeightEngine passed dial: " << name << " that it does not understand."); } params[dial - 1] = startval; } void OscWeightEngine::SetDialValue(int nuisenum, double val) { #ifdef DEBUG_OSC_WE std::cout << "SetDial: " << (nuisenum % 1000) << " at " << val << std::endl; #endif fHasChanged = (params[(nuisenum % 1000) - 1] - val) > std::numeric_limits::epsilon(); params[(nuisenum % 1000) - 1] = val; } void OscWeightEngine::SetDialValue(std::string name, double val) { #ifdef DEBUG_OSC_WE std::cout << "SetDial: " << name << " at " << val << std::endl; #endif int dial = SystEnumFromString(name); if (!dial) { THROW("OscWeightEngine passed dial: " << name << " that it does not understand."); } fHasChanged = (params[dial - 1] - val) > std::numeric_limits::epsilon(); params[dial - 1] = val; } bool OscWeightEngine::IsDialIncluded(std::string name) { return SystEnumFromString(name); } bool OscWeightEngine::IsDialIncluded(int nuisenum) { return ((nuisenum % 1000) > 0) && ((nuisenum % 1000) < 6); } double OscWeightEngine::GetDialValue(std::string name) { int dial = SystEnumFromString(name); if (!dial) { THROW("OscWeightEngine passed dial: " << name << " that it does not understand."); } return params[dial - 1]; } double OscWeightEngine::GetDialValue(int nuisenum) { if (!(nuisenum % 1000) || (nuisenum % 1000) > 6) { THROW("OscWeightEngine passed dial enum: " << (nuisenum % 1000) << " that it does not understand, expected [1,6]."); } return params[(nuisenum % 1000) - 1]; } void OscWeightEngine::Reconfigure(bool silent) { fHasChanged = false; }; bool OscWeightEngine::NeedsEventReWeight() { if (fHasChanged) { return true; } return false; } double OscWeightEngine::CalcWeight(BaseFitEvt* evt) { static bool Warned = false; if (evt->probe_E == 0xdeadbeef) { if (!Warned) { ERROR(WRN, "Oscillation weights asked for but using 'litemode' or " "unsupported generator input. Pasuing for 10..."); sleep(10); Warned = true; } return 1; } return CalcWeight(evt->probe_E * 1E-3, evt->probe_pdg); } double OscWeightEngine::CalcWeight(double ENu, int PDGNu, int TargetPDGNu) { if (LengthParam == 0xdeadbeef) { // not configured. return 1; } #ifdef __PROB3PP_ENABLED__ int NuType = (ForceFromNuPDG != 0) ? ForceFromNuPDG : GetNuType(PDGNu); bp.SetMNS(params[theta12_idx], params[theta13_idx], params[theta23_idx], params[dm12_idx], params[dm23_idx], params[dcp_idx], ENu, true, NuType); int pmt = 0; double prob_weight = 1; TargetPDGNu = (TargetPDGNu == -1) ? (TargetNuType ? TargetNuType : NuType) : GetNuType(TargetPDGNu); if (LengthParamIsZenith) { // Use earth density bp.DefinePath(LengthParam, 0); bp.propagate(NuType); pmt = 0; prob_weight = bp.GetProb(NuType, TargetPDGNu); } else { if (constant_density != 0xdeadbeef) { bp.propagateLinear(NuType, LengthParam, constant_density); pmt = 1; prob_weight = bp.GetProb(NuType, TargetPDGNu); } else { pmt = 2; prob_weight = bp.GetVacuumProb(NuType, TargetPDGNu, ENu * 1E-3, LengthParam); } } #ifdef DEBUG_OSC_WE if (prob_weight != prob_weight) { THROW("Calculated bad prob weight: " << prob_weight << "(Osc Type: " << pmt << " -- " << NuType << " -> " << TargetPDGNu << ")"); } if (prob_weight > 1) { THROW("Calculated bad prob weight: " << prob_weight << "(Osc Type: " << pmt << " -- " << NuType << " -> " << TargetPDGNu << ")"); } std::cout << NuType << " -> " << TargetPDGNu << ": " << ENu << " = " << prob_weight << "%%." << std::endl; #endif return prob_weight; #else return 1; #endif } int OscWeightEngine::SystEnumFromString(std::string const& name) { if (name == "dm23") { return 1; } else if (name == "sinsq_theta23") { return 2; } else if (name == "sinsq_theta13") { return 3; } else if (name == "dm12") { return 4; } else if (name == "sinsq_theta12") { return 5; } else if (name == "dcp") { return 6; } else { return 0; } } void OscWeightEngine::Print() { std::cout << "OscWeightEngine: " << std::endl; std::cout << "\t theta12: " << params[theta12_idx] << std::endl; std::cout << "\t theta13: " << params[theta13_idx] << std::endl; std::cout << "\t theta23: " << params[theta23_idx] << std::endl; std::cout << "\t dm12: " << params[dm12_idx] << std::endl; std::cout << "\t dm23: " << params[dm23_idx] << std::endl; std::cout << "\t dcp: " << params[dcp_idx] << std::endl; } diff --git a/src/Reweight/OscWeightEngine.h b/src/Reweight/OscWeightEngine.h index 8cacff7..774fb99 100644 --- a/src/Reweight/OscWeightEngine.h +++ b/src/Reweight/OscWeightEngine.h @@ -1,143 +1,144 @@ // 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 "FitLogger.h" #include "FitEvent.h" #include "PhysConst.h" #include "WeightEngineBase.h" #ifdef __PROB3PP_ENABLED__ #include "BargerPropagator.h" #endif #include #ifdef __PROB3PP_ENABLED__ class BG : public BargerPropagator { public: BG() : BargerPropagator(){}; double GetBaseline() { return Earth->get_Pathlength(); } }; #endif class OscWeightEngine : public WeightEngineBase { enum params { dm23_idx = 0, theta23_idx, theta13_idx, dm12_idx, theta12_idx, dcp_idx, }; #ifdef __PROB3PP_ENABLED__ BG bp; #endif //******************************* Osc params ****************************** double theta12; double theta13; double theta23; /// The 1-2 mass squared splitting (small) [eV] double dm12; /// The 2-3 mass squared splitting (large) [eV] double dm23; /// The PMNS CP-violating phase double dcp; ///\brief The constant matter density used for simple given baseline /// oscillation [g/cm^3] double constant_density; /// Whether LengthParam corresponds to a Zenith or a baseline. /// /// If we just want to calculate the osc. prob. with a constant matter density /// then this should be false and constant_density should be set /// (or 0 for vacuum prob). bool LengthParamIsZenith; /// Either a path length or a post oscillation zenith angle /// /// N.B. For a beamline that has a dip angle of X degrees, the post /// oscillation zenith angle will be 90+X degrees. double LengthParam; /// Holds current value of oscillation parameters. double params[6]; /// The oscillation target type /// /// If unspecified in the element, it will default to /// disappearance probability. int TargetNuType; /// The initial neutrino species /// /// If unspecified in the element, it will be determined by /// the incoming events. int ForceFromNuPDG; public: OscWeightEngine(); /// Configures oscillation parameters from input xml file. /// /// Osc parameters configured from OscParam XML element as: /// /// /// /// If matter_density and baseline are present, then oscillation probability /// is calculated for a constant matter density. /// If detection_zenith_deg is present, then the baseline and density are /// calculated from the density profile and radius of the earth. /// If none are present, a vacuum oscillation is calculated. /// If TargetNuPDG is unspecified, oscillation will default to /// disappearance probability. void Config(); // Functions requiring Override void IncludeDial(std::string name, double startval); void SetDialValue(int nuisenum, double val); void SetDialValue(std::string name, double val); bool IsDialIncluded(std::string name); bool IsDialIncluded(int nuisenum); double GetDialValue(std::string name); double GetDialValue(int nuisenum); void Reconfigure(bool silent); bool NeedsEventReWeight(); double CalcWeight(BaseFitEvt* evt); + /// ENu [GeV] double CalcWeight(double ENu, int PDGNu, int TargetPDGNu = -1); static int SystEnumFromString(std::string const& name); void Print(); }; diff --git a/src/Routines/MinimizerRoutines.cxx b/src/Routines/MinimizerRoutines.cxx index ee3e995..369cf45 100755 --- a/src/Routines/MinimizerRoutines.cxx +++ b/src/Routines/MinimizerRoutines.cxx @@ -1,1506 +1,1511 @@ // 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 "MinimizerRoutines.h" #include "Simple_MH_Sampler.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," "DataToys,MCMC"); }; //************************************* 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 xmlcmds; std::vector configargs; // Make easier to handle arguments. std::vector 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.Set("cardfile", fCardFile); + if (!fCardFile.empty()) fCompKey.Set("cardfile", fCardFile); if (!fOutputFile.empty()) fCompKey.Set("outputfile", fOutputFile); - if (!fStrategy.empty()) fCompKey.Set("strategy", fStrategy); + if (!fStrategy.empty()) fCompKey.Set("strategy", fStrategy); // Load XML Cardfile - configuration.LoadSettings( fCompKey.GetS("cardfile"), ""); + configuration.LoadSettings(fCompKey.GetS("cardfile"), ""); // 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.FinaliseSettings(fCompKey.GetS("outputfile") + ".xml"); // Add Error Verbo Lines verbocount += Config::GetParI("VERBOSITY"); errorcount += Config::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 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; // Override state if none given if (!key.Has("state")) { key.SetS("state", "FIX"); } 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) { + double opnom = parnom; + double oparstep = parstep; parnom = FitBase::RWAbsToSigma(partype, parname, parnom); parlow = FitBase::RWAbsToSigma(partype, parname, parlow); parhigh = FitBase::RWAbsToSigma(partype, parname, parhigh); - parstep = FitBase::RWAbsToSigma(partype, parname, parstep); + parstep = + FitBase::RWAbsToSigma(partype, parname, opnom + parstep) - parnom; + std::cout << "ParStep: " << parstep << " (" << oparstep << ")." + << std::endl; } 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 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 + "_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 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 = ""; bool UseMCMC = false; // 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 = ""; } else if (!routine.compare("MCMC")) { UseMCMC = true; } // make minimizer if (fMinimizer) delete fMinimizer; if (UseMCMC) { fMinimizer = new Simple_MH_Sampler(); } else { 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& 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.find("DataToys") != std::string::npos) ThrowDataToys(); 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") or !routine.compare("MCMC")) { 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 Config::SetPar("input.maxevents", lowstatsevents); Config::SetPar("VERBOSITY", 3); SetupFCN(); RunFitRoutine(trueroutine); Config::SetPar("input.maxevents", maxevents); SetupFCN(); Config::SetPar("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 nameVect; std::vector valVect; std::vector errVect; std::vector minVect; std::vector maxVect; std::vector startVect; std::vector endfixVect; std::vector 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 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; }; void MinimizerRoutines::ThrowDataToys() { LOG(FIT) << "Generating Toy Data Throws" << std::endl; int verb = Config::GetParI("VERBOSITY"); SETVERBOSITY(FIT); int nthrows = FitPar::Config().GetParI("NToyThrows"); double maxlike = -1.0; double minlike = -1.0; std::vector values; for (int i = 0; i < 1.E4; i++) { fSampleFCN->ThrowDataToy(); double like = fSampleFCN->GetLikelihood(); values.push_back(like); if (maxlike == -1.0 or like > maxlike) maxlike = like; if (minlike == -1.0 or like < minlike) minlike = like; } SETVERBOSITY(verb); // Fill Histogram TH1D* likes = new TH1D("toydatalikelihood", "toydatalikelihood", int(sqrt(nthrows)), minlike, maxlike); for (size_t i = 0; i < values.size(); i++) { likes->Fill(values[i]); } // Save to file LOG(FIT) << "Writing toy data throws" << std::endl; fOutputRootFile->cd(); likes->Write(); } diff --git a/src/Utils/PlotUtils.cxx b/src/Utils/PlotUtils.cxx index 4eabf23..ec742e3 100644 --- a/src/Utils/PlotUtils.cxx +++ b/src/Utils/PlotUtils.cxx @@ -1,999 +1,1001 @@ // 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(); + delete rootHistFile; return tempHist; } TH1* PlotUtils::GetTH1FromRootFile(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"); if (!rootHistFile || rootHistFile->IsZombie()) { THROW("Couldn't open root file: \"" << file << "\"."); } TH1* tempHist = dynamic_cast(rootHistFile->Get(name.c_str())->Clone()); if (!tempHist) { THROW("Couldn't retrieve: \"" << name << "\" from root file: \"" << file << "\"."); } tempHist->SetDirectory(0); rootHistFile->Close(); + delete rootHistFile; return tempHist; } /// Returns a vector of named TH1*s found in a single input file. /// /// Expects a descriptor like: file.root[hist1|hist2|...] std::vector PlotUtils::GetTH1sFromRootFile( std::string const& descriptor) { std::vector descriptors = GeneralUtils::ParseToStr(descriptor, ","); std::vector hists; for (size_t d_it = 0; d_it < descriptors.size(); ++d_it) { std::string& d = descriptors[d_it]; std::vector fname = GeneralUtils::ParseToStr(d, "["); if (!fname.size() || !fname[0].length()) { THROW("Couldn't find input file when attempting to parse : \"" << d << "\". Expected input.root[hist1|hist2|...]."); } if (fname[1][fname[1].length() - 1] == ']') { fname[1] = fname[1].substr(0, fname[1].length() - 1); } std::vector histnames = GeneralUtils::ParseToStr(fname[1], "|"); if (!histnames.size()) { THROW( "Couldn't find any histogram name specifiers when attempting to " "parse " ": \"" << fname[1] << "\". Expected hist1|hist2|..."); } TFile* rootHistFile = new TFile(fname[0].c_str(), "READ"); if (!rootHistFile || rootHistFile->IsZombie()) { THROW("Couldn't open root file: \"" << fname[0] << "\"."); } for (size_t i = 0; i < histnames.size(); ++i) { TH1* tempHist = dynamic_cast(rootHistFile->Get(histnames[i].c_str())->Clone()); if (!tempHist) { THROW("Couldn't retrieve: \"" << histnames[i] << "\" from root file: \"" << fname[0] << "\"."); } tempHist->SetDirectory(0); hists.push_back(tempHist); } rootHistFile->Close(); } return hists; } 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; }