diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt index 3983c0a..c260ba6 100644 --- a/app/CMakeLists.txt +++ b/app/CMakeLists.txt @@ -1,178 +1,192 @@ # 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/Utils) 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() if(USE_GENIE) add_executable(PrepareGENIE PrepareGENIE.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareGENIE) target_link_libraries(PrepareGENIE ${MODULETargets}) target_link_libraries(PrepareGENIE ${CMAKE_DEPENDLIB_FLAGS}) # target_link_libraries(PrepareGENIE ${ROOT_LIBS}) if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "") set_target_properties(PrepareGENIE PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS}) endif() endif() if(USE_NEUT) add_executable(PrepareNEUT PrepareNEUT.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareNEUT) target_link_libraries(PrepareNEUT ${MODULETargets}) target_link_libraries(PrepareNEUT ${CMAKE_DEPENDLIB_FLAGS}) # target_link_libraries(PrepareNEUT ${ROOT_LIBS}) if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "") set_target_properties(PrepareNEUT PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS}) endif() endif() # PREPARE NUWRO # Commented out for the time being until it is finished.. if(USE_NuWro) add_executable(PrepareNuwro PrepareNuwroEvents.cxx) set(TARGETS_TO_BUILD ${TARGETS_TO_BUILD};PrepareNuwro) target_link_libraries(PrepareNuwro ${MODULETargets}) target_link_libraries(PrepareNuwro ${CMAKE_DEPENDLIB_FLAGS}) # target_link_libraries(PrepareNuwro ${ROOT_LIBS}) if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "") set_target_properties(PrepareNuwro PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS}) endif() endif() install(TARGETS ${TARGETS_TO_BUILD} DESTINATION bin) diff --git a/app/nuisflat.cxx b/app/nuisflat.cxx index 4be70c8..d0772bf 100644 --- a/app/nuisflat.cxx +++ b/app/nuisflat.cxx @@ -1,176 +1,177 @@ // 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 "ComparisonRoutines.h" #include "InputUtils.h" - +#include "MeasurementBase.h" +#include "GenericFlux_Tester.h" #include "Smearceptance_Tester.h" // Global Arguments std::string gOptInputFile = ""; std::string gOptFormat = ""; std::string gOptOutputFile = ""; std::string gOptType = "DEFAULT"; std::string gOptNumberEvents = "NULL"; std::string gOptCardInput = ""; std::string gOptOptions = ""; //******************************* void PrintSyntax() { //******************************* std::cout << "nuisflat -i input -f format [-o outfile] [-n nevents] [-t " "options] [-q con=val] \n"; std::cout << "\n Arguments : " << "\n\t -i input : Path to input vector of events to flatten" << "\n\t" << "\n\t This should be given in the same format a normal " "input file" << "\n\t is given to NUISANCE. {e.g. NUWRO:eventsout.root}." << "\n\t" << "\n\t -f format : FlatTree format to output:" << "\n\t\t GenericFlux : Standard event summary format." << "\n\t " << "\n\t[-c crd.xml]: Input card file to override configs or set dial values." << "\n\t " << "\n\t[-o outfile]: Optional output file path. " << "\n\t " << "\n\t If none given, input.format.root is chosen." << "\n\t" << "\n\t[-n nevents]: Optional choice of Nevents to run over. Default is " "all." << "\n\t" << "\n\t[-t options]: Pass OPTION to the FlatTree sample. " << "\n\t Similar to type field in comparison xml configs." << "\n\t" << "\n\t[-q con=val]: Configuration overrides." << std::endl; exit(-1); }; //____________________________________________________________________________ void GetCommandLineArgs(int argc, char** argv) { // Check for -h flag. for (int i = 0; i < argc; i++) { if ((!std::string(argv[i]).compare("-h")) || (!std::string(argv[i]).compare("-?")) || (!std::string(argv[i]).compare("--help"))) PrintSyntax(); } // Format is nuwro -r run_number -n n events std::vector args = GeneralUtils::LoadCharToVectStr(argc, argv); // Parse input file ParserUtils::ParseArgument(args, "-i", gOptInputFile, false); if (gOptInputFile == "") { THROW("Need to provide a valid input file to nuisflat using -i flag!"); } else { LOG(FIT) << "Reading Input File = " << gOptInputFile << std::endl; } // Get Output Format ParserUtils::ParseArgument(args, "-f", gOptFormat, false); if (gOptFormat == "") { THROW("Need to provide a valid output format to nuisflat!"); } else { LOG(FIT) << "Saving flattree in format = " << gOptFormat << std::endl; } // Get Output File ParserUtils::ParseArgument(args, "-o", gOptOutputFile, false); if (gOptOutputFile == "") { gOptOutputFile = gOptInputFile + "." + gOptFormat + ".root"; LOG(FIT) << "No output file given so saving nuisflat output to:" << gOptOutputFile << std::endl; } else { LOG(FIT) << "Saving nuisflat output to " << gOptOutputFile << std::endl; } // Get N Events and Configs nuisconfig configuration = Config::Get(); ParserUtils::ParseArgument(args, "-n", gOptNumberEvents, false); if (gOptNumberEvents.compare("NULL")) { configuration.OverrideConfig("MAXEVENTS=" + gOptNumberEvents); } std::vector configargs; ParserUtils::ParseArgument(args, "-q", configargs); for (size_t i = 0; i < configargs.size(); i++) { configuration.OverrideConfig(configargs[i]); } ParserUtils::ParseArgument(args, "-c", gOptCardInput, false); if (gOptCardInput != "") { QLOG(FIT, "Reading cardfile: " << gOptCardInput); configuration.LoadConfig(gOptCardInput, ""); } ParserUtils::ParseArgument(args, "-t", gOptOptions, false); if (gOptOptions != "") { QLOG(FIT, "Read options: \"" << gOptOptions << "\'"); } return; } //******************************* int main(int argc, char* argv[]) { //******************************* // Parse GetCommandLineArgs(argc, argv); // Make output file TFile* f = new TFile(gOptOutputFile.c_str(), "RECREATE"); if (f->IsZombie()) { THROW("Cannot create output file!"); } f->cd(); FitPar::Config().out = f; // Create a new measurementbase class depending on the Format MeasurementBase* flattreecreator = NULL; // Make a new sample key for the format of interest. nuiskey samplekey = Config::CreateKey("sample"); if (!gOptFormat.compare("GenericFlux")) { samplekey.AddS("name", "FlatTree"); samplekey.AddS("input", gOptInputFile); samplekey.AddS("type", gOptType); flattreecreator = new GenericFlux_Tester("FlatTree", gOptInputFile, FitBase::GetRW(), gOptType, ""); } else { ERR(FTL) << "Unknown FlatTree format!" << std::endl; } // Make the FlatTree reconfigure flattreecreator->Reconfigure(); f->cd(); flattreecreator->Write(); f->Close(); // Show Final Status LOG(FIT) << "-------------------------------------" << std::endl; LOG(FIT) << "Flattree Generation Complete." << std::endl; LOG(FIT) << "-------------------------------------" << std::endl; return 0; } diff --git a/cmake/MINERvASetup.cmake b/cmake/MINERvASetup.cmake new file mode 100644 index 0000000..0ccfcd9 --- /dev/null +++ b/cmake/MINERvASetup.cmake @@ -0,0 +1,19 @@ +# 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 . +################################################################################ +LIST(APPEND EXTRA_CXX_FLAGS -D__MINERVA_RW_ENABLED__) diff --git a/cmake/ReweightEnginesSetup.cmake b/cmake/ReweightEnginesSetup.cmake index 44c80ff..87f47d4 100644 --- a/cmake/ReweightEnginesSetup.cmake +++ b/cmake/ReweightEnginesSetup.cmake @@ -1,82 +1,87 @@ # 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 . ################################################################################ ################################## NEUT ###################################### if(USE_NEUT) include(${CMAKE_SOURCE_DIR}/cmake/NEUTSetup.cmake) cmessage(STATUS "Using NEUT Reweight engine.") set(USE_NEUT TRUE CACHE BOOL "Whether to enable NEUT (reweight) support. Requires external libraries. " FORCE) endif() ################################# NuWro ###################################### if(USE_NuWro) include(${CMAKE_SOURCE_DIR}/cmake/NuWroSetup.cmake) cmessage(STATUS "Using NuWro Reweight engine.") set(USE_NuWro TRUE CACHE BOOL "Whether to enable NuWro support. " FORCE) endif() ################################## GENIE ##################################### if(USE_GENIE) include(${CMAKE_SOURCE_DIR}/cmake/GENIESetup.cmake) cmessage(STATUS "Using GENIE Reweight engine.") set(USE_GENIE TRUE CACHE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. " FORCE) endif() ################################## NIWG ###################################### if(USE_NIWG) include(${CMAKE_SOURCE_DIR}/cmake/NIWGSetup.cmake) cmessage(STATUS "Using NIWG Reweight engine.") set(USE_NIWG TRUE CACHE BOOL "Whether to enable (T2K) NIWG ReWeight support. Requires external libraries. " FORCE) endif() ################################## T2K ###################################### if(USE_T2K) include(${CMAKE_SOURCE_DIR}/cmake/T2KSetup.cmake) cmessage(STATUS "Using T2K Reweight engine.") set(USE_T2K TRUE CACHE BOOL "Whether to enable T2KReWeight support. Requires external libraries. " FORCE) endif() - +################################## MINERvA ###################################### +if(USE_MINERvA_RW) + include(${CMAKE_SOURCE_DIR}/cmake/MINERvASetup.cmake) + cmessage(STATUS "Using MINERvA Reweight engine.") + set(USE_MINERvA_RW TRUE CACHE BOOL "Whether to enable MINERvA ReWeight support. " FORCE) +endif() ################################################################################ ################################ Prob3++ #################################### include(${CMAKE_SOURCE_DIR}/cmake/Prob3++Setup.cmake) ################################################################################ cmessage(STATUS "Reweight engine include directories: ${RWENGINE_INCLUDE_DIRECTORIES}") if(NEED_ROOTEVEGEN) cmessage(STATUS "Require ROOT eve generation libraries") LIST(REVERSE ROOT_LIBS) LIST(APPEND ROOT_LIBS Gui Ged Geom TreePlayer EG Eve) LIST(REVERSE ROOT_LIBS) endif() if(NEED_ROOTPYTHIA6) cmessage(STATUS "Require ROOT Pythia6 libraries") LIST(APPEND ROOT_LIBS EGPythia6 Pythia6) endif() LIST(APPEND EXTRA_LIBS ${ROOT_LIBS}) diff --git a/cmake/cacheVariables.cmake b/cmake/cacheVariables.cmake index 5e88e9c..1c75a96 100644 --- a/cmake/cacheVariables.cmake +++ b/cmake/cacheVariables.cmake @@ -1,206 +1,208 @@ # Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret ################################################################################ # This file is part of NUISANCE. # # NUISANCE is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # NUISANCE is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with NUISANCE. If not, see . ################################################################################ function(CheckAndSetDefaultEnv VARNAME DEFAULT CACHETYPE DOCSTRING ENVNAME) #cmessage(DEBUG "Trying to assign variable ${VARNAME} into the cache.") if(NOT DEFINED ${VARNAME}) if(DEFINED ENV{${ENVNAME}} AND NOT $ENV{${ENVNAME}} STREQUAL "") set(${VARNAME} $ENV{${ENVNAME}} CACHE ${CACHETYPE} ${DOCSTRING}) cmessage(DEBUG " Read ${VARNAME} from ENVVAR ${ENVNAME} as $ENV{${ENVNAME}}.") else() set(${VARNAME} ${DEFAULT} CACHE ${CACHETYPE} ${DOCSTRING}) endif() else() set(${VARNAME} ${${VARNAME}} CACHE ${CACHETYPE} ${DOCSTRING}) unset(${VARNAME}) endif() cmessage(CACHE "--Set cache variable: \"${VARNAME}\" to \"${${VARNAME}}\", in cache ${CACHETYPE}.") endfunction() function(CheckAndSetDefaultCache VARNAME DEFAULT CACHETYPE DOCSTRING) # cmessage(DEBUG "Trying to assign variable ${VARNAME} into the cache.") if(NOT DEFINED ${VARNAME}) set(${VARNAME} ${DEFAULT} CACHE ${CACHETYPE} ${DOCSTRING}) else() set(${VARNAME} ${${VARNAME}} CACHE ${CACHETYPE} ${DOCSTRING}) unset(${VARNAME}) endif() cmessage(CACHE "--Set cache variable: \"${VARNAME}\" to \"${${VARNAME}}\", in cache ${CACHETYPE}.") endfunction() function(CheckAndSetDefault VARNAME DEFAULT) # cmessage(DEBUG "Trying to assign variable ${VARNAME}.") if(NOT DEFINED ${VARNAME}) set(${VARNAME} ${DEFAULT} PARENT_SCOPE) set(${VARNAME} ${DEFAULT}) endif() cmessage(CACHE "--Set variable: \"${VARNAME}\" to \"${${VARNAME}}\".") endfunction() CheckAndSetDefaultCache(VERBOSE TRUE BOOL "Whether to configure loudly.") set (CMAKE_SKIP_BUILD_RPATH TRUE) #Changes default install path to be a subdirectory of the build dir. #Can set build dir at configure time with -DCMAKE_INSTALL_PREFIX=/install/path if(CMAKE_INSTALL_PREFIX STREQUAL "" OR CMAKE_INSTALL_PREFIX STREQUAL "/usr/local") set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}") elseif(NOT DEFINED CMAKE_INSTALL_PREFIX) set(CMAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/${CMAKE_SYSTEM_NAME}") endif() if(CMAKE_BUILD_TYPE STREQUAL "") set(CMAKE_BUILD_TYPE DEBUG) elseif(NOT DEFINED CMAKE_BUILD_TYPE) set(CMAKE_BUILD_TYPE DEBUG) endif() CheckAndSetDefaultCache(USE_MINIMIZER TRUE INTERNAL "Whether we are using the ROOT minimization libraries. ") CheckAndSetDefaultCache(USE_HEPMC FALSE BOOL "Whether to enable HepMC input support. ") CheckAndSetDefaultEnv(HEPMC "" PATH "Path to HepMC source tree root directory. Overrides environment variable \$HEPMC <>" HEPMC) CheckAndSetDefaultCache(HEPMC_MOMUNIT "GEV" STRING "HepMC momentum units [MEV|GEV]. ") CheckAndSetDefaultCache(HEPMC_LENUNIT "CM" STRING "HepMC momentum units [MM|CM]. ") CheckAndSetDefaultCache(HEPMC_USED_EP FALSE INTERNAL "Whether we built HepMC or not. ") CheckAndSetDefaultCache(USE_NEUT FALSE BOOL "Whether to enable NEUT (reweight) support. Requires external libraries. ") CheckAndSetDefaultEnv(NEUT_ROOT "" PATH "Path to NEUT source tree root directory. Overrides environment variable \$NEUT_ROOT <>" NEUT_ROOT) CheckAndSetDefaultEnv(CERN "" PATH "Path to CERNLIB source tree root directory that NEUT was built against. Overrides environment variable \$CERN <>" CERN) CheckAndSetDefaultEnv(CERN_LEVEL "" STRING "CERNLIB Library version. Overrides environment variable \$CERN_LEVEL <>" CERN_LEVEL) CheckAndSetDefaultCache(USE_NuWro FALSE BOOL "Whether to enable NuWro support. ") CheckAndSetDefaultEnv(NUWRO "" PATH "Path to NuWro source tree root directory. Overrides environment variable \$NUWRO <>" NUWRO) CheckAndSetDefaultEnv(NUWRO_INC "" PATH "Path to NuWro installed includes directory, needs to contain \"params_all.h\". Overrides environment variable \$NUWRO_INC <>" NUWRO_INC) CheckAndSetDefaultCache(NUWRO_INPUT_FILE "" FILEPATH "Path to an input NuWro event vector, which can be used to build NuWro i/o libraries. <>") CheckAndSetDefaultCache(NUWRO_BUILT_FROM_FILE FALSE INTERNAL "Whether the NuWro libraries were built by NUISANCE. ") CheckAndSetDefaultCache(USE_NuWro_RW FALSE BOOL "Whether to try and build support for NuWro reweighting. ") CheckAndSetDefaultCache(USE_NuWro_SRW_Event FALSE BOOL "Whether to use cut down NuWro reweight event format. Requires NuWro reweight. ") CheckAndSetDefaultCache(USE_GENIE FALSE BOOL "Whether to enable GENIE (reweight) support. Requires external libraries. ") CheckAndSetDefaultEnv(GENIE "" PATH "Path to GENIE source tree root directory. Overrides environment variable \$GENIE <>" GENIE) CheckAndSetDefaultEnv(LHAPDF_LIB "" PATH "Path to pre-built LHAPDF libraries. Overrides environment variable \$LHAPDF_LIB. <>" LHAPDF_LIB) CheckAndSetDefaultEnv(LHAPDF_INC "" PATH "Path to installed LHAPDF headers. Overrides environment variable \$LHAPDF_INC. <>" LHAPDF_INC) CheckAndSetDefaultEnv(LHAPATH "" PATH "Path to LHA PDF inputs. Overrides environment variable \$LHAPATH. <>" LHAPATH) CheckAndSetDefaultEnv(LIBXML2_LIB "" PATH "Path to pre-built LIBXML2 libraries. Overrides environment variable \$LIBXML2_LIB. <>" LIBXML2_LIB) CheckAndSetDefaultEnv(LIBXML2_INC "" PATH "Path to installed LIBXML2 headers. Overrides environment variable \$LIBXML2_INC. <>" LIBXML2_INC) CheckAndSetDefaultEnv(LOG4CPP_LIB "" PATH "Path to pre-built LOG4CPP libraries. Overrides environment variable \$LOG4CPP_LIB. <>" LOG4CPP_LIB) CheckAndSetDefaultEnv(LOG4CPP_INC "" PATH "Path to installed LOG4CPP headers. Overrides environment variable \$LOG4CPP_INC. <>" LOG4CPP_INC) CheckAndSetDefaultCache(BUILD_GEVGEN FALSE BOOL "Whether to build nuisance_gevgen app.") CheckAndSetDefaultCache(USE_T2K FALSE BOOL "Whether to enable T2KReWeight support. Requires external libraries. ") CheckAndSetDefaultEnv(T2KREWEIGHT "" PATH "Path to installed T2KREWEIGHTReWeight. Overrides environment variable \$T2KREWEIGHT. <>" T2KREWEIGHT) CheckAndSetDefaultCache(USE_NIWG FALSE BOOL "Whether to enable (T2K) NIWG ReWeight support. Requires external libraries. ") CheckAndSetDefaultEnv(NIWG_ROOT "" PATH "Path to installed NIWGReWeight. Overrides environment variable \$NIWG. <>" NIWG) +CheckAndSetDefaultCache(USE_MINERvA_RW FALSE BOOL "Whether to enable MINERvA ReWeight support. ") + CheckAndSetDefaultEnv(PYTHIA6 "" PATH "Path to directory containing libPythia6.so. Overrides environment variable \$PYTHIA6 <>" PYTHIA6) CheckAndSetDefaultEnv(PYTHIA8 "" PATH "Path to directory containing libPythia8.so. Overrides environment variable \$PYTHIA8 <>" PYTHIA8) CheckAndSetDefaultCache(USE_PYTHIA8 FALSE BOOL "Whether to enable PYTHIA8 event support. ") CheckAndSetDefaultCache(USE_GiBUU TRUE BOOL "Whether to enable GiBUU event support. ") CheckAndSetDefaultCache(BUILD_GiBUU FALSE BOOL "Whether to build supporting GiBUU event tools along with a patched version of GiBUU. ") CheckAndSetDefaultCache(USE_NUANCE TRUE BOOL "Whether to enable NUANCE event support. ") CheckAndSetDefaultCache(USE_PROB3PP FALSE BOOL "Whether to download and compile in Prob3++ support. ") CheckAndSetDefaultCache(NO_EXTERNAL_UPDATE TRUE BOOL "Whether to perform the update target for external dependencies. ") CheckAndSetDefaultCache(USE_GPERFTOOLS FALSE BOOL "Whether to compile in google performance tools. ") CheckAndSetDefault(NEED_PYTHIA6 FALSE) CheckAndSetDefault(NEED_PYTHIA8 FALSE) CheckAndSetDefault(NEED_ROOTEVEGEN FALSE) CheckAndSetDefault(NEED_ROOTPYTHIA6 FALSE) CheckAndSetDefaultCache(USE_OMP FALSE BOOL "Whether to enable multicore features (there currently are none...). ") CheckAndSetDefault(NO_EXPERIMENTS FALSE) cmessage(STATUS "NO_EXPERIMENTS: ${NO_EXPERIMENTS}") CheckAndSetDefaultCache(NO_ANL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build ANL samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_ArgoNeuT ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build ArgoNeuT samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_BEBC ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build BEBC samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_BNL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build BNL samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_FNAL ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build FNAL samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_GGM ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build GGM samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_K2K ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build K2K samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_MINERvA ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build MINERvA samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_MiniBooNE ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build MiniBooNE samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_T2K ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build T2K samples. <-DNO_EXPERIMENTS=FALSE>") CheckAndSetDefaultCache(NO_SciBooNE ${NO_EXPERIMENTS} BOOL "Whether to *NOT* build SciBooNE samples. <-DNO_EXPERIMENTS=FALSE>") function(SAYVARS) LIST(APPEND VARS USE_HEPMC HEPMC HEPMC_MOMUNIT HEPMC_LENUNIT HEPMC_USED_EP USE_NEUT NEUT_ROOT CERN CERN_LEVEL USE_NuWro NUWRO NUWRO_INC NUWRO_INPUT_FILE NUWRO_BUILT_FROM_FILE USE_GENIE GENIE LHAPDF_LIB LHAPDF_INC LIBXML2_LIB LIBXML2_INC LOG4CPP_LIB GENIE_LOG4CPP_INC BUILD_GEVGEN USE_T2K USE_NIWG USE_GiBUU BUILD_GiBUU USE_NUANCE NO_EXTERNAL_UPDATE USE_GPERFTOOLS NO_ANL NO_ArgoNeuT NO_BEBC NO_BNL NO_FNAL NO_GGM NO_K2K NO_MINERvA NO_MiniBooNE NO_T2K NO_SciBooNE) foreach(v ${VARS}) if(DEFINED ${v}) cmessage(DEBUG "VARIABLE: \"${v}\" = \"${${v}}\"") endif() endforeach(v) endfunction() diff --git a/src/FCN/JointFCN.cxx b/src/FCN/JointFCN.cxx index a4136e2..31023b2 100755 --- a/src/FCN/JointFCN.cxx +++ b/src/FCN/JointFCN.cxx @@ -1,1087 +1,1010 @@ #include "JointFCN.h" #include #include "FitUtils.h" -// //*************************************************** -// JointFCN::JointFCN(std::string cardfile, TFile* outfile) { -// //*************************************************** - -// fOutputDir = gDirectory; -// FitPar::Config().out = outfile; - -// fCardFile = cardfile; - -// LoadSamples(fCardFile); - -// fCurIter = 0; -// fMCFilled = false; - -// fOutputDir->cd(); - -// fIterationTree = NULL; -// fDialVals = NULL; -// fNDials = 0; - -// fUsingEventManager = FitPar::Config().GetParB("EventManager"); -// fOutputDir->cd(); -// }; +//*************************************************** JointFCN::JointFCN(TFile* outfile) { +//*************************************************** + fOutputDir = gDirectory; if (outfile) FitPar::Config().out = outfile; std::vector samplekeys = Config::QueryKeys("sample"); LoadSamples(samplekeys); std::vector covarkeys = Config::QueryKeys("covar"); LoadPulls(covarkeys); fCurIter = 0; fMCFilled = false; - fIterationTree = NULL; + fIterationTree = false; fDialVals = NULL; fNDials = 0; fUsingEventManager = FitPar::Config().GetParB("EventManager"); fOutputDir->cd(); } +//*************************************************** JointFCN::JointFCN(std::vector samplekeys, TFile* outfile) { +//*************************************************** + fOutputDir = gDirectory; if (outfile) FitPar::Config().out = outfile; LoadSamples(samplekeys); fCurIter = 0; fMCFilled = false; fOutputDir->cd(); - fIterationTree = NULL; + fIterationTree = false; fDialVals = NULL; fNDials = 0; fUsingEventManager = FitPar::Config().GetParB("EventManager"); fOutputDir->cd(); } //*************************************************** JointFCN::~JointFCN() { //*************************************************** // Delete Samples for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; delete exp; } for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull* pull = *iter; delete pull; } // Sort Tree if (fIterationTree) DestroyIterationTree(); if (fDialVals) delete fDialVals; if (fSampleLikes) delete fSampleLikes; }; //*************************************************** void JointFCN::CreateIterationTree(std::string name, FitWeight* rw) { - //*************************************************** - - LOG(FIT) << " Creating new iteration tree! " << std::endl; - if (fIterationTree && !name.compare(fIterationTree->GetName())) { - fIterationTree->Reset(); - return; - } - - fIterationTree = new TTree(name.c_str(), name.c_str()); - - // Setup Main Branches - fIterationTree->Branch("total_likelihood", &fLikelihood, - "total_likelihood/D"); - fIterationTree->Branch("total_ndof", &fNDOF, "total_ndof/I"); +//*************************************************** - // Setup Sample Arrays - int ninputs = fSamples.size() + fPulls.size(); - fSampleLikes = new double[ninputs]; - fSampleNDOF = new int[ninputs]; - fNDOF = GetNDOF(); + LOG(FIT) << " Creating new iteration container! " << std::endl; + DestroyIterationTree(); + fIterationTreeName = name; - // Setup Sample Branches - int count = 0; - for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); + // Add sample likelihoods and ndof + for (MeasListConstIter iter = fSamples.begin(); + iter != fSamples.end(); iter++) { - MeasurementBase* exp = *iter; + MeasurementBase* exp = *iter; std::string name = exp->GetName(); - std::string liketag = name + "_likelihood"; - std::string ndoftag = name + "_ndof"; - fIterationTree->Branch(liketag.c_str(), &fSampleLikes[count], - (liketag + "/D").c_str()); - fIterationTree->Branch(ndoftag.c_str(), &fSampleNDOF[count], - (ndoftag + "/D").c_str()); + std::string liketag = name + "_likelihood"; + fNameValues.push_back(liketag); + fCurrentValues.push_back(0.0); - count++; + std::string ndoftag = name + "_ndof"; + fNameValues.push_back(ndoftag); + fCurrentValues.push_back(0.0); } - for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { - ParamPull* pull = *iter; + // Add Pull terms + for (PullListConstIter iter = fPulls.begin(); + iter != fPulls.end(); iter++) { + ParamPull* pull = *iter; std::string name = pull->GetName(); + std::string liketag = name + "_likelihood"; + fNameValues.push_back(liketag); + fCurrentValues.push_back(0.0); + std::string ndoftag = name + "_ndof"; + fNameValues.push_back(ndoftag); + fCurrentValues.push_back(0.0); + } - fIterationTree->Branch(liketag.c_str(), &fSampleLikes[count], - (liketag + "/D").c_str()); - fIterationTree->Branch(ndoftag.c_str(), &fSampleNDOF[count], - (ndoftag + "/D").c_str()); + // Add Likelihoods + fNameValues.push_back("total_likelihood"); + fCurrentValues.push_back(0.0); - count++; - } + fNameValues.push_back("total_ndof"); + fCurrentValues.push_back(0.0); - // Add Dial Branches - std::vector dials = rw->GetDialNames(); - fNDials = dials.size(); - fDialVals = new double[fNDials]; + // Setup Containers + fSampleN = fSamples.size() + fPulls.size(); + fSampleLikes = new double[fSampleN]; + fSampleNDOF = new int[fSampleN]; - for (int i = 0; i < fNDials; i++) { - fIterationTree->Branch(dials[i].c_str(), &fDialVals[i], - (dials[i] + "/D").c_str()); + // Add Dials + std::vector dials = rw->GetDialNames(); + for (size_t i = 0; i < dials.size(); i++){ + fNameValues.push_back( dials[i] ); + fCurrentValues.push_back( 0.0 ); } + fNDials = dials.size(); + fDialVals = new double[fNDials]; + + // Set IterationTree Flag + fIterationTree = true; + } //*************************************************** void JointFCN::DestroyIterationTree() { - //*************************************************** +//*************************************************** + + fIterationCount.clear(); + fCurrentValues.clear(); + fNameValues.clear(); + fIterationValues.clear(); - if (!fIterationTree) { - delete fIterationTree; - } } //*************************************************** void JointFCN::WriteIterationTree() { - //*************************************************** +//*************************************************** + LOG(FIT) << "Writing iteration tree" << std::endl; + + // Make a new TTree + TTree* itree = new TTree(fIterationTreeName.c_str(), + fIterationTreeName.c_str()); + + double* vals = new double[fNameValues.size()]; + int count = 0; + + itree->Branch("iteration",&count,"Iteration/I"); + for (int i = 0; i < fNameValues.size(); i++) { + itree->Branch( fNameValues[i].c_str(), + &vals[i], + (fNameValues[i] + "/D").c_str() ); + } + + // Fill Iterations + for (size_t i = 0; i < fIterationValues.size(); i++){ + std::vector itervals = fIterationValues[i]; + + // Fill iteration state + count = fIterationCount[i]; + for (size_t j = 0; j < itervals.size(); j++){ + vals[j] = itervals[j]; + } - if (!fIterationTree) { - ERR(FTL) << "Can't save empty iteration tree!" << std::endl; - throw; + // Save to TTree + itree->Fill(); } - fIterationTree->Write(); + + // Write to file + itree->Write(); } //*************************************************** void JointFCN::FillIterationTree(FitWeight* rw) { - //*************************************************** +//*************************************************** - if (!fIterationTree) { - ERR(FTL) << "Trying to fill iteration_tree when it is NULL!" << std::endl; - throw; + // Loop over samples count + int count = 0; + for (int i = 0; i < fSampleN; i++){ + fCurrentValues[count++] = fSampleLikes[i]; + fCurrentValues[count++] = double(fSampleNDOF[i]); } + // Fill Totals + fCurrentValues[count++] = fLikelihood; + fCurrentValues[count++] = double(fNDOF); + + // Loop Over Parameter Counts rw->GetAllDials(fDialVals, fNDials); - fIterationTree->Fill(); + for (int i = 0; i < fNDials; i++){ + fCurrentValues[count++] = double(fDialVals[i]); + } + + // Push Back Into Container + fIterationCount.push_back( fCurIter ); + fIterationValues.push_back(fCurrentValues); } //*************************************************** double JointFCN::DoEval(const double* x) { //*************************************************** // WEIGHT ENGINE fDialChanged = FitBase::GetRW()->HasRWDialChanged(x); FitBase::GetRW()->UpdateWeightEngine(x); if (fDialChanged) { FitBase::GetRW()->Reconfigure(); FitBase::EvtManager().ResetWeightFlags(); } if (LOG_LEVEL(REC)) { FitBase::GetRW()->Print(); } // SORT SAMPLES ReconfigureSamples(); // GET TEST STAT fLikelihood = GetLikelihood(); - + fNDOF = GetNDOF(); + // PRINT PROGRESS LOG(FIT) << "Current Stat (iter. " << this->fCurIter << ") = " << fLikelihood << std::endl; // UPDATE TREE if (fIterationTree) FillIterationTree(FitBase::GetRW()); return fLikelihood; } //*************************************************** int JointFCN::GetNDOF() { //*************************************************** int totaldof = 0; int count = 0; // Total number of Free bins in each MC prediction for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; int dof = exp->GetNDOF(); // Save Seperate DOF if (fIterationTree) { fSampleNDOF[count] = dof; } // Add to total totaldof += dof; count++; } // Loop over pulls for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull* pull = *iter; double dof = pull->GetLikelihood(); // Save seperate DOF if (fIterationTree) { fSampleNDOF[count] = dof; } // Add to total totaldof += dof; count++; } // Set Data Variable - fNDOF = totaldof; + fSampleNDOF[count] = totaldof; return totaldof; } //*************************************************** double JointFCN::GetLikelihood() { //*************************************************** LOG(MIN) << std::left << std::setw(43) << "Getting likelihoods..." << " : " << "-2logL" << std::endl; // Loop and add up likelihoods in an uncorrelated way double like = 0.0; int count = 0; for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; double newlike = exp->GetLikelihood(); int ndof = exp->GetNDOF(); // Save seperate likelihoods if (fIterationTree) { fSampleLikes[count] = newlike; } LOG(MIN) << "-> " << std::left << std::setw(40) << exp->GetName() << " : " << newlike << "/" << ndof << std::endl; // Add Weight Scaling // like *= FitBase::GetRW()->GetSampleLikelihoodWeight(exp->GetName()); // Add to total like += newlike; count++; } // Loop over pulls for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull* pull = *iter; double newlike = pull->GetLikelihood(); // Save seperate likelihoods if (fIterationTree) { fSampleLikes[count] = newlike; } // Add to total like += newlike; count++; } // Set Data Variable fLikelihood = like; return like; }; void JointFCN::LoadSamples(std::vector samplekeys) { LOG(MIN) << "Loading Samples : " << samplekeys.size() << std::endl; for (size_t i = 0; i < samplekeys.size(); i++) { nuiskey key = samplekeys[i]; // Get Sample Options std::string samplename = key.GetS("name"); std::string samplefile = key.GetS("input"); std::string sampletype = key.GetS("type"); std::string fakeData = ""; LOG(MIN) << "Loading Sample : " << samplename << std::endl; fOutputDir->cd(); MeasurementBase* NewLoadedSample = SampleUtils::CreateSample(key); if (!NewLoadedSample) { ERR(FTL) << "Could not load sample provided: " << samplename << std::endl; ERR(FTL) << "Check spelling with that in src/FCN/SampleList.cxx" << std::endl; throw; } else { fSamples.push_back(NewLoadedSample); } } } +//*************************************************** void JointFCN::LoadPulls(std::vector pullkeys) { +//*************************************************** for (size_t i = 0; i < pullkeys.size(); i++) { nuiskey key = pullkeys[i]; std::string pullname = key.GetS("name"); std::string pullfile = key.GetS("input"); std::string pulltype = key.GetS("type"); fOutputDir->cd(); - std::cout << "Creating Pull Term : " << std::endl; - sleep(1); fPulls.push_back(new ParamPull(pullname, pullfile, pulltype)); } - - // // Sample Inputs - // if (!samplevect[0].compare("covar") || !samplevect[0].compare("pulls") - // || - // !samplevect[0].compare("throws")) { - // // Get all inputs - // std::string name = samplevect[1]; - // std::string files = samplevect[2]; - // std::string type = "DEFAULT"; - - // if (samplevect.size() > 3) { - // type = samplevect[3]; - // } else if (!samplevect[0].compare("pull")) { - // type = "GAUSPULL"; - // } else if (!samplevect[0].compare("throws")) { - // type = "GAUSTHROWS"; - // } - - // // Create Pull Class - // LOG(MIN) << "Loading up pull term: " << name << " << " << files << " - // (" - // << type << ")" << std::endl; - // std::string fakeData = ""; - // fOutputDir->cd(); - // fPulls.push_back(new ParamPull(name, files, type)); - // } } -// //*************************************************** -// void JointFCN::LoadSamples(std::string cardinput) -// //*************************************************** -// { -// LOG(MIN) << "Initializing Samples" << std::endl; - -// // Read the card file here and load objects -// std::string line; -// std::ifstream card(cardinput.c_str(), ifstream::in); - -// // Make sure they are created in correct working DIR -// fOutputDir->cd(); - -// while (std::getline(card >> std::ws, line, '\n')) { -// // Skip Empties -// if (line.c_str()[0] == '#') continue; -// if (line.empty()) continue; - -// // Parse line -// std::vector samplevect = GeneralUtils::ParseToStr(line, " -// "); - -// // Sample Inputs -// if (!samplevect[0].compare("sample")) { -// // Get all inputs -// std::string name = samplevect[1]; -// std::string files = samplevect[2]; -// std::string type = "DEFAULT"; -// if (samplevect.size() > 3) type = samplevect[3]; - -// // Create Sample Class -// LOG(MIN) << "Loading up sample: " << name << " << " << files << " (" -// << type << ")" << std::endl; -// std::string fakeData = ""; -// fOutputDir->cd(); -// MeasurementBase* NewLoadedSample = SampleUtils::CreateSample(name, -// files, type, -// fakeData, FitBase::GetRW()); - -// if (!NewLoadedSample) { -// ERR(FTL) << "Could not load sample provided: " << name << std::endl; -// ERR(FTL) << "Check spelling with that in src/FCN/SampleList.cxx" -// << std::endl; -// throw; -// } else { -// fSamples.push_back(NewLoadedSample); -// } -// } - -// // Sample Inputs -// if (!samplevect[0].compare("covar") || !samplevect[0].compare("pulls") || -// !samplevect[0].compare("throws")) { -// // Get all inputs -// std::string name = samplevect[1]; -// std::string files = samplevect[2]; -// std::string type = "DEFAULT"; - -// if (samplevect.size() > 3) { -// type = samplevect[3]; -// } else if (!samplevect[0].compare("pull")) { -// type = "GAUSPULL"; -// } else if (!samplevect[0].compare("throws")) { -// type = "GAUSTHROWS"; -// } - -// // Create Pull Class -// LOG(MIN) << "Loading up pull term: " << name << " << " << files << " (" -// << type << ")" << std::endl; -// std::string fakeData = ""; -// fOutputDir->cd(); -// fPulls.push_back(new ParamPull(name, files, type)); -// } -// } -// card.close(); -// }; - //*************************************************** void JointFCN::ReconfigureSamples(bool fullconfig) { - //*************************************************** +//*************************************************** int starttime = time(NULL); LOG(REC) << "Starting Reconfigure iter. " << this->fCurIter << std::endl; // std::cout << fUsingEventManager << " " << fullconfig << " " << fMCFilled << // std::endl; // Event Manager Reconf if (fUsingEventManager) { if (!fullconfig and fMCFilled) ReconfigureFastUsingManager(); else ReconfigureUsingManager(); } else { // Loop over all Measurement Classes for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; // If RW Either do signal or full reconfigure. if (fDialChanged or !fMCFilled or fullconfig) { if (!fullconfig and fMCFilled) exp->ReconfigureFast(); else exp->Reconfigure(); // If RW Not needed just do normalisation } else { exp->Renormalise(); } } } // Loop over pulls and update for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull* pull = *iter; pull->Reconfigure(); } fMCFilled = true; LOG(MIN) << "Finished Reconfigure iter. " << fCurIter << " in " << time(NULL) - starttime << "s" << std::endl; fCurIter++; } //*************************************************** void JointFCN::ReconfigureSignal() { - //*************************************************** - this->ReconfigureSamples(false); +//*************************************************** + ReconfigureSamples(false); } //*************************************************** void JointFCN::ReconfigureAllEvents() { //*************************************************** FitBase::GetRW()->Reconfigure(); FitBase::EvtManager().ResetWeightFlags(); ReconfigureSamples(true); } std::vector JointFCN::GetInputList() { std::vector InputList; fIsAllSplines = true; MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); std::vector subsamples = exp->GetSubSamples(); for (size_t i = 0; i < subsamples.size(); i++) { InputHandlerBase* inp = subsamples[i]->GetInput(); if (std::find(InputList.begin(), InputList.end(), inp) == InputList.end()) { if (subsamples[i]->GetInput()->GetType() != kSPLINEPARAMETER) fIsAllSplines = false; InputList.push_back(subsamples[i]->GetInput()); } } } return InputList; } std::vector JointFCN::GetSubSampleList() { std::vector SampleList; MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); std::vector subsamples = exp->GetSubSamples(); for (size_t i = 0; i < subsamples.size(); i++) { SampleList.push_back(subsamples[i]); } } return SampleList; } //*************************************************** void JointFCN::ReconfigureUsingManager() { - //*************************************************** +//*************************************************** // 'Slow' Event Manager Reconfigure LOG(REC) << "Event Manager Reconfigure" << std::endl; int timestart = time(NULL); // Reset all samples MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); exp->ResetAll(); } // If we are siving signal, reset all containers. bool savesignal = (FitPar::Config().GetParB("SignalReconfigures")); if (savesignal) { // Reset all of our event signal vectors fSignalEventBoxes.clear(); fSignalEventFlags.clear(); fSampleSignalFlags.clear(); fSignalEventSplines.clear(); } // Make sure we have a list of inputs if (fInputList.empty()) { fInputList = GetInputList(); fSubSampleList = GetSubSampleList(); } // If all inputs are splines make sure the readers are told // they need to be reconfigured. std::vector::iterator inp_iter = fInputList.begin(); if (fIsAllSplines) { for (; inp_iter != fInputList.end(); inp_iter++) { InputHandlerBase* curinput = (*inp_iter); // Tell reader in each BaseEvent it needs a Reconfigure next weight calc. BaseFitEvt* curevent = curinput->FirstBaseEvent(); if (curevent->fSplineRead) { curevent->fSplineRead->SetNeedsReconfigure(true); } } } // MAIN INPUT LOOP ==================== int fillcount = 0; int inputcount = 0; inp_iter = fInputList.begin(); // Loop over each input in manager for (; inp_iter != fInputList.end(); inp_iter++) { InputHandlerBase* curinput = (*inp_iter); // Get event information FitEvent* curevent = curinput->FirstNuisanceEvent(); curinput->CreateCache(); int i = 0; int nevents = curinput->GetNEvents(); int countwidth = nevents / 5; // Start event loop iterating until we get a NULL pointer. while (curevent) { // Get Event Weight curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent); curevent->Weight = curevent->RWWeight * curevent->InputWeight; double rwweight = curevent->Weight; // std::cout << "RWWeight = " << curevent->RWWeight << " " << // curevent->InputWeight << std::endl; // Logging // std::cout << CHECKLOG(1) << std::endl; if (LOGGING(REC)) { if (i % countwidth == 0) { QLOG(REC, curinput->GetName() - << " : Processed " << i << " events. [M, W] = [" - << curevent->Mode << ", " << rwweight << "]"); + << " : Processed " << i << " events. [M, W] = [" + << curevent->Mode << ", " << rwweight << "]"); } } // Setup flag for if signal found in at least one sample bool foundsignal = false; // Create a new signal bitset for this event std::vector signalbitset(fSubSampleList.size()); // Create a new signal box vector for this event std::vector signalboxes; // Start measurement iterator size_t measitercount = 0; std::vector::iterator meas_iter = - fSubSampleList.begin(); + fSubSampleList.begin(); // Loop over all subsamples (sub in JointMeas) for (; meas_iter != fSubSampleList.end(); meas_iter++) { MeasurementBase* curmeas = (*meas_iter); // Compare input pointers, to current input, skip if not. // Pointer tells us if it matches without doing ID checks. if (curinput != curmeas->GetInput()) { if (savesignal) { // Set bit to 0 as definitely not signal signalbitset[measitercount] = 0; } // Count up what measurement we are on. measitercount++; // Skip sample as input not signal. continue; } // Fill events for matching inputs. MeasurementVariableBox* box = curmeas->FillVariableBox(curevent); bool signal = curmeas->isSignal(curevent); curmeas->SetSignal(signal); curmeas->FillHistograms(curevent->Weight); // If its Signal tally up fills if (signal) { fillcount++; } // If we are saving signal/splines fill the bitset if (savesignal) { signalbitset[measitercount] = signal; } // If signal save a clone of the event box for use later. if (savesignal and signal) { foundsignal = true; signalboxes.push_back(box->CloneSignalBox()); } // Keep track of Measurement we are on. measitercount++; } // Once we've filled the measurements, if saving signal // push back if any sample flagged this event as signal if (savesignal) { fSignalEventFlags.push_back(foundsignal); } // Save the vector of signal boxes for this event if (savesignal and foundsignal) { fSignalEventBoxes.push_back(signalboxes); fSampleSignalFlags.push_back(signalbitset); } // If all inputs are splines we can save the spline coefficients // for fast in memory reconfigures later. if (fIsAllSplines and savesignal and foundsignal) { // Make temp vector to push back with std::vector coeff; for (size_t l = 0; l < (UInt_t)curevent->fSplineRead->GetNPar(); l++) { coeff.push_back(curevent->fSplineCoeff[l]); } // Push back to signal event splines. Kept in sync with // fSignalEventBoxes size. // int splinecount = fSignalEventSplines.size(); fSignalEventSplines.push_back(coeff); // if (splinecount % 1000 == 0) { // std::cout << "Pushed Back Coeff " << splinecount << " : "; // for (size_t l = 0; l < fSignalEventSplines[splinecount].size(); l++) // { // std::cout << " " << fSignalEventSplines[splinecount][l]; // } // std::cout << std::endl; // } } // Clean up vectors once done with this event signalboxes.clear(); signalbitset.clear(); // Iterate to the next event. curevent = curinput->NextNuisanceEvent(); i++; } curinput->RemoveCache(); // Keep track of what input we are on. inputcount++; } // End of Event Loop =============================== // Now event loop is finished loop over all Measurements // Converting Binned events to XSec Distributions iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); exp->ConvertEventRates(); } // Print out statements on approximate memory usage for profiling. LOG(REC) << "Filled " << fillcount << " signal events." << std::endl; if (savesignal) { int mem = - ( // sizeof(fSignalEventBoxes) + - // fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) + - sizeof(MeasurementVariableBox1D) * fillcount) * - 1E-6; + ( // sizeof(fSignalEventBoxes) + + // fSignalEventBoxes.size() * sizeof(fSignalEventBoxes.at(0)) + + sizeof(MeasurementVariableBox1D) * fillcount) * + 1E-6; LOG(REC) << " -> Saved " << fillcount << " signal boxes for faster access. (~" << mem << " MB)" << std::endl; if (fIsAllSplines and !fSignalEventSplines.empty()) { int splmem = sizeof(float) * fSignalEventSplines.size() * fSignalEventSplines[0].size() * 1E-6; LOG(REC) << " -> Saved " << fillcount << " " << fSignalEventSplines.size() << " spline sets into memory. (~" << splmem << " MB)" << std::endl; } } LOG(REC) << "Time taken ReconfigureUsingManager() : " << time(NULL) - timestart << std::endl; // Check SignalReconfigures works for all samples - if (savesignal){ + if (savesignal) { double likefull = GetLikelihood(); ReconfigureFastUsingManager(); double likefast = GetLikelihood(); - + if (fabs(likefull - likefast) > 0.0001) - { - ERROR(FTL,"Fast and Full Likelihoods DIFFER! : " << likefull << " : " << likefast); - ERROR(FTL,"This means some samples you are using are not setup to use SignalReconfigures=1"); - ERROR(FTL,"Please turn OFF signal reconfigures."); - throw; - } else { + { + ERROR(FTL, "Fast and Full Likelihoods DIFFER! : " << likefull << " : " << likefast); + ERROR(FTL, "This means some samples you are using are not setup to use SignalReconfigures=1"); + ERROR(FTL, "Please turn OFF signal reconfigures."); + throw; + } else { LOG(FIT) << "Likelihoods for FULL and FAST match. Will use FAST next time." << std::endl; } } // End of reconfigure return; }; //*************************************************** void JointFCN::ReconfigureFastUsingManager() { - //*************************************************** +//*************************************************** LOG(FIT) << " -> Doing FAST using manager" << std::endl; // Get Start time for profilling int timestart = time(NULL); // Reset all samples MeasListConstIter iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); exp->ResetAll(); } // Check for saved variables if not do a full reconfigure. if (fSignalEventFlags.empty()) { ERR(WRN) << "Signal Flags Empty! Using normal manager." << std::endl; ReconfigureUsingManager(); return; } bool fFillNuisanceEvent = - FitPar::Config().GetParB("FullEventOnSignalReconfigure"); + FitPar::Config().GetParB("FullEventOnSignalReconfigure"); // Setup fast vector iterators. std::vector::iterator inpsig_iter = fSignalEventFlags.begin(); std::vector >::iterator box_iter = - fSignalEventBoxes.begin(); + fSignalEventBoxes.begin(); std::vector >::iterator spline_iter = - fSignalEventSplines.begin(); + fSignalEventSplines.begin(); std::vector >::iterator samsig_iter = - fSampleSignalFlags.begin(); + fSampleSignalFlags.begin(); int splinecount = 0; // Setup stuff for logging int fillcount = 0; int nevents = fSignalEventFlags.size(); int countwidth = nevents / 20; // If All Splines tell splines they need a reconfigure. std::vector::iterator inp_iter = fInputList.begin(); if (fIsAllSplines) { LOG(REC) << "All Spline Inputs so using fast spline loop." << std::endl; for (; inp_iter != fInputList.end(); inp_iter++) { InputHandlerBase* curinput = (*inp_iter); // Tell each fSplineRead in BaseFitEvent to reconf next weight calc BaseFitEvt* curevent = curinput->FirstBaseEvent(); if (curevent->fSplineRead) curevent->fSplineRead->SetNeedsReconfigure(true); } } // Loop over all possible spline inputs double* coreeventweights = new double[fSignalEventBoxes.size()]; splinecount = 0; inp_iter = fInputList.begin(); inpsig_iter = fSignalEventFlags.begin(); spline_iter = fSignalEventSplines.begin(); // Loop over all signal flags // For each valid signal flag add one to splinecount // Get Splines from that count and add to weight // Add splinecount int sigcount = 0; splinecount = 0; // #pragma omp parallel for shared(splinecount,sigcount) for (uint iinput = 0; iinput < fInputList.size(); iinput++) { InputHandlerBase* curinput = fInputList[iinput]; BaseFitEvt* curevent = curinput->FirstBaseEvent(); for (int i = 0; i < curinput->GetNEvents(); i++) { double rwweight = 0.0; if (fSignalEventFlags[sigcount]) { // Get Event Info if (!fIsAllSplines) { if (fFillNuisanceEvent) curinput->GetNuisanceEvent(i); else curevent = curinput->GetBaseEvent(i); } else { curevent->fSplineCoeff = &fSignalEventSplines[splinecount][0]; } curevent->RWWeight = FitBase::GetRW()->CalcWeight(curevent); curevent->Weight = curevent->RWWeight * curevent->InputWeight; rwweight = curevent->Weight; coreeventweights[splinecount] = rwweight; if (splinecount % countwidth == 0) { LOG(REC) << "Processed " << splinecount << " event weights. W = " << rwweight << std::endl; } // #pragma omp atomic splinecount++; } // #pragma omp atomic sigcount++; } } LOG(SAM) << "Processed event weights." << std::endl; // #pragma omp barrier // Reset Iterators inpsig_iter = fSignalEventFlags.begin(); spline_iter = fSignalEventSplines.begin(); box_iter = fSignalEventBoxes.begin(); samsig_iter = fSampleSignalFlags.begin(); int nsplineweights = splinecount; splinecount = 0; // Start of Fast Event Loop ============================ // Start input iterators // Loop over number of inputs for (int ispline = 0; ispline < nsplineweights; ispline++) { double rwweight = coreeventweights[ispline]; // Get iterators for this event std::vector::iterator subsamsig_iter = (*samsig_iter).begin(); std::vector::iterator subbox_iter = - (*box_iter).begin(); + (*box_iter).begin(); // Loop over all sub measurements. std::vector::iterator meas_iter = fSubSampleList.begin(); for (; meas_iter != fSubSampleList.end(); meas_iter++, subsamsig_iter++) { MeasurementBase* curmeas = (*meas_iter); // If event flagged as signal for this sample fill from the box. if (*subsamsig_iter) { curmeas->SetSignal(true); curmeas->FillHistogramsFromBox((*subbox_iter), rwweight); // Move onto next box if there is one. subbox_iter++; fillcount++; } } if (ispline % countwidth == 0) { LOG(REC) << "Filled " << ispline << " sample weights." << std::endl; } // Iterate over the main signal event containers. samsig_iter++; box_iter++; spline_iter++; splinecount++; } // End of Fast Event Loop =================== LOG(SAM) << "Filled sample distributions." << std::endl; // Now loop over all Measurements // Convert Binned events iterSam = fSamples.begin(); for (; iterSam != fSamples.end(); iterSam++) { MeasurementBase* exp = (*iterSam); exp->ConvertEventRates(); } // Cleanup coreeventweights if (fIsAllSplines) { delete coreeventweights; } // Print some reconfigure profiling. LOG(REC) << "Filled " << fillcount << " signal events." << std::endl; LOG(REC) << "Time taken ReconfigureFastUsingManager() : " << time(NULL) - timestart << std::endl; } //*************************************************** void JointFCN::Write() { - //*************************************************** +//*************************************************** // Save a likelihood/ndof plot LOG(MIN) << "Writing likelihood plot.." << std::endl; std::vector likes; std::vector ndofs; std::vector names; for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); - iter++){ + iter++) { MeasurementBase* exp = *iter; double like = exp->GetLikelihood(); double ndof = exp->GetNDOF(); std::string name = exp->GetName(); likes.push_back(like); ndofs.push_back(ndof); names.push_back(name); } - TH1D likehist = TH1D("likelihood_hist","likelihood_hist", - likes.size(), 0.0, double(likes.size())); - TH1D ndofhist = TH1D("ndof_hist","ndof_hist", - ndofs.size(), 0.0, double(ndofs.size())); - TH1D divhist = TH1D("likedivndof_hist","likedivndof_hist", - likes.size(), 0.0, double(likes.size())); - for (size_t i = 0; i < likehist.GetNbinsX(); i++){ - likehist.SetBinContent(i+1, likes[i]); - ndofhist.SetBinContent(i+1, ndofs[i]); - if (ndofs[i] != 0.0){ - divhist.SetBinContent(i+1, likes[i]/ndofs[i]); + TH1D likehist = TH1D("likelihood_hist", "likelihood_hist", + likes.size(), 0.0, double(likes.size())); + TH1D ndofhist = TH1D("ndof_hist", "ndof_hist", + ndofs.size(), 0.0, double(ndofs.size())); + TH1D divhist = TH1D("likedivndof_hist", "likedivndof_hist", + likes.size(), 0.0, double(likes.size())); + for (size_t i = 0; i < likehist.GetNbinsX(); i++) { + likehist.SetBinContent(i + 1, likes[i]); + ndofhist.SetBinContent(i + 1, ndofs[i]); + if (ndofs[i] != 0.0) { + divhist.SetBinContent(i + 1, likes[i] / ndofs[i]); } - likehist.GetXaxis()->SetBinLabel(i+1, names[i].c_str()); - ndofhist.GetXaxis()->SetBinLabel(i+1, names[i].c_str()); - divhist.GetXaxis()->SetBinLabel(i+1, names[i].c_str()); + likehist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str()); + ndofhist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str()); + divhist.GetXaxis()->SetBinLabel(i + 1, names[i].c_str()); } likehist.Write(); ndofhist.Write(); divhist.Write(); // Loop over individual experiments and call Write LOG(MIN) << "Writing each of the data classes..." << std::endl; for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; exp->Write(); } // Save Pull Terms for (PullListConstIter iter = fPulls.begin(); iter != fPulls.end(); iter++) { ParamPull* pull = *iter; pull->Write(); } if (FitPar::Config().GetParB("EventManager")) { // Get list of inputs std::map fInputs = - FitBase::EvtManager().GetInputs(); + FitBase::EvtManager().GetInputs(); std::map::const_iterator iterInp; for (iterInp = fInputs.begin(); iterInp != fInputs.end(); iterInp++) { InputHandlerBase* input = (iterInp->second); input->GetFluxHistogram()->Write(); input->GetXSecHistogram()->Write(); input->GetEventHistogram()->Write(); } } }; //*************************************************** void JointFCN::SetFakeData(std::string fakeinput) { - //*************************************************** +//*************************************************** LOG(MIN) << "Setting fake data from " << fakeinput << std::endl; for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; exp->SetFakeDataValues(fakeinput); } return; } //*************************************************** void JointFCN::ThrowDataToy() { - //*************************************************** +//*************************************************** for (MeasListConstIter iter = fSamples.begin(); iter != fSamples.end(); iter++) { MeasurementBase* exp = *iter; exp->ThrowDataToy(); } return; } diff --git a/src/FCN/JointFCN.h b/src/FCN/JointFCN.h index 361a75a..bf87424 100755 --- a/src/FCN/JointFCN.h +++ b/src/FCN/JointFCN.h @@ -1,165 +1,176 @@ #ifndef _JOINT_FCN_H_ #define _JOINT_FCN_H_ /*! * \addtogroup FCN * @{ */ #include #include #include #include // ROOT headers #include "TTree.h" #include "TH1D.h" #include "TH2D.h" #include #include "TGraphErrors.h" #include #include #include "FitUtils.h" // External fitter headers #include "MeasurementBase.h" #include "SampleList.h" #define GetCurrentDir getcwd #include "EventManager.h" #include "Math/Functor.h" #include "ParamPull.h" #include "NuisConfig.h" #include "NuisKey.h" #include "MeasurementVariableBox.h" #include "MeasurementVariableBox1D.h" using namespace FitUtils; using namespace FitBase; //! Main FCN Class which ROOT's joint function needs to evaulate the chi2 at each stage of the fit. class JointFCN { public: //! Constructor //! cardfile = Path to input card file listing samples JointFCN(std::vector samplekeys, TFile* outfile=NULL); JointFCN(TFile* outfile=NULL); // Loads from global config //! Destructor ~JointFCN(); //! Create sample list from cardfile void LoadSamples(std::vector samplekeys); void LoadPulls(std::vector pullkeys); //! Main Likelihood evaluation FCN double DoEval(const double *x); //! Func Wrapper for ROOT inline double operator() (const std::vector & x) { double* x_array = new double[x.size()]; return this->DoEval(x_array); }; //! Func Wrapper for ROOT inline double operator() (const double *x) { return this->DoEval(x); }; //! Create a TTree to save all dial value iterations for this FCN void CreateIterationTree(std::string name, FitWeight* rw); //! Fills dial values and sample likelihoods into tree void FillIterationTree(FitWeight* rw); //! Writes TTree to fOutput directory void WriteIterationTree(); //! Deletes TTree void DestroyIterationTree(); //! Get Degrees of Freedom for samples (NBins) int GetNDOF(); //! Return NDOF wrapper inline unsigned int NDim() {return this->GetNDOF();}; //! Reconfigure samples where we force all events to be looped over. void ReconfigureAllEvents() ; //! Call Reconfigure on samples. //! Choice of reconfigure type depends on whether dials have changed //! and the MC has been filled. void ReconfigureSamples(bool fullconfig = false); //! Call reconfigure on only signal events (defaults to all events if CurIter=0) void ReconfigureSignal(); //! Gets likelihood for all samples in FCN (assuming uncorrelated) double GetLikelihood(); //! Returns list of pointers to the samples inline std::list GetSampleList(){ return fSamples; } //! Return list of pointers to all the pulls inline std::list GetPullList(){ return fPulls; }; //! Write all samples to output DIR void Write(); //! Set Fake data from file/MC void SetFakeData(std::string fakeinput); //! Reconfigure looping over duplicate inputs void ReconfigureUsingManager(); //! Reconfigure Fast looping over duplicate inputs void ReconfigureFastUsingManager(); /// Throws data according to current stats void ThrowDataToy(); std::vector GetSubSampleList(); std::vector GetInputList(); private: //! Append the experiments to include in the fit to this list std::list fSamples; //! Append parameter pull terms to include penalties in the fit to this list std::list fPulls; TDirectory *fOutputDir; //!< Directory to save contents std::string fCardFile; //!< Input Card text file bool fDialChanged; //!< Flag for if RW dials changed UInt_t fCurIter; //!< Counter for how many times reconfigure called bool fMCFilled; //!< Check MC has at least been filled once - TTree* fIterationTree; //!< Tree to save RW values on each function call + bool fIterationTree; //!< Tree to save RW values on each function call int fNDials; //!< Number of RW Dials in FitWeight double* fDialVals; //!< Current Values of RW Dials double fLikelihood; //!< Current likelihood for joint sample likelihood double fNDOF; //!< Total NDOF double* fSampleLikes; //!< Likelihoods for each individual measurement in list int * fSampleNDOF; //!< NDOF for each individual measurement in list bool fUsingEventManager; //!< Flag for doing joint comparisons std::vector< std::vector > fSignalEventSplines; std::vector< std::vector > fSignalEventBoxes; std::vector< bool > fSignalEventFlags; std::vector< std::vector > fSampleSignalFlags; std::vector fInputList; std::vector fSubSampleList; bool fIsAllSplines; + + std::vector< int > fIterationCount; + std::vector< double > fCurrentValues; + std::vector< std::string > fNameValues; + std::vector< std::vector > fIterationValues; + int fSampleN; + std::string fIterationTreeName; + + + + }; /*! @} */ #endif // _JOINT_FCN_H_ diff --git a/src/FCN/SampleList.cxx b/src/FCN/SampleList.cxx index a92ffa2..3db2359 100644 --- a/src/FCN/SampleList.cxx +++ b/src/FCN/SampleList.cxx @@ -1,725 +1,1016 @@ #include "SampleList.h" +#ifndef __NO_ANL__ +#include "ANL_CCQE_Evt_1DQ2_nu.h" +#include "ANL_CCQE_XSec_1DEnu_nu.h" + +// ANL CC1ppip +#include "ANL_CC1ppip_Evt_1DQ2_nu.h" +#include "ANL_CC1ppip_Evt_1DcosmuStar_nu.h" +#include "ANL_CC1ppip_Evt_1DcosmuStar_nu.h" +#include "ANL_CC1ppip_Evt_1DcosthAdler_nu.h" +#include "ANL_CC1ppip_Evt_1Dphi_nu.h" +#include "ANL_CC1ppip_Evt_1Dppi_nu.h" +#include "ANL_CC1ppip_Evt_1Dthpr_nu.h" +#include "ANL_CC1ppip_XSec_1DEnu_nu.h" +#include "ANL_CC1ppip_XSec_1DQ2_nu.h" +// ANL CC1npip +#include "ANL_CC1npip_Evt_1DQ2_nu.h" +#include "ANL_CC1npip_Evt_1DcosmuStar_nu.h" +#include "ANL_CC1npip_Evt_1Dppi_nu.h" +#include "ANL_CC1npip_XSec_1DEnu_nu.h" +// ANL CC1pi0 +#include "ANL_CC1pi0_Evt_1DQ2_nu.h" +#include "ANL_CC1pi0_Evt_1DcosmuStar_nu.h" +#include "ANL_CC1pi0_XSec_1DEnu_nu.h" +// ANL NC1npip (mm, exotic!) +#include "ANL_NC1npip_Evt_1Dppi_nu.h" +// ANL NC1ppim (mm, exotic!) +#include "ANL_NC1ppim_Evt_1DcosmuStar_nu.h" +#include "ANL_NC1ppim_XSec_1DEnu_nu.h" +// ANL CC2pi 1pim1pip (mm, even more exotic!) +#include "ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu.h" +#include "ANL_CC2pi_1pim1pip_Evt_1Dppim_nu.h" +#include "ANL_CC2pi_1pim1pip_Evt_1Dppip_nu.h" +#include "ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu.h" +#include "ANL_CC2pi_1pim1pip_XSec_1DEnu_nu.h" +// ANL CC2pi 1pip1pip (mm, even more exotic!) +#include "ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu.h" +#include "ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu.h" +#include "ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu.h" +#include "ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu.h" +#include "ANL_CC2pi_1pip1pip_XSec_1DEnu_nu.h" +// ANL CC2pi 1pip1pi0 (mm, even more exotic!) +#include "ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu.h" +#include "ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu.h" +#include "ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu.h" +#include "ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu.h" +#include "ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu.h" +#endif + +#ifndef __NO_ArgoNeuT__ +// ArgoNeuT CC-inclusive +#include "ArgoNeuT_CCInc_XSec_1Dpmu_antinu.h" +#include "ArgoNeuT_CCInc_XSec_1Dpmu_nu.h" +#include "ArgoNeuT_CCInc_XSec_1Dthetamu_antinu.h" +#include "ArgoNeuT_CCInc_XSec_1Dthetamu_nu.h" +#endif + +#ifndef __NO_BNL__ +// BNL CCQE +#include "BNL_CCQE_Evt_1DQ2_nu.h" +#include "BNL_CCQE_XSec_1DEnu_nu.h" +// BNL CC1ppip +#include "BNL_CC1ppip_Evt_1DQ2_nu.h" +#include "BNL_CC1ppip_Evt_1DQ2_nu.h" +#include "BNL_CC1ppip_Evt_1DcosthAdler_nu.h" +#include "BNL_CC1ppip_Evt_1Dphi_nu.h" +#include "BNL_CC1ppip_XSec_1DEnu_nu.h" +// BNL CC1npip +#include "BNL_CC1npip_Evt_1DQ2_nu.h" +#include "BNL_CC1npip_XSec_1DEnu_nu.h" +// BNL CC1pi0 +#include "BNL_CC1pi0_Evt_1DQ2_nu.h" +#include "BNL_CC1pi0_XSec_1DEnu_nu.h" +#endif + +#ifndef __NO_FNAL__ +// FNAL CCQE +#include "FNAL_CCQE_Evt_1DQ2_nu.h" +// FNAL CC1ppip +#include "FNAL_CC1ppip_Evt_1DQ2_nu.h" +#include "FNAL_CC1ppip_XSec_1DEnu_nu.h" +#include "FNAL_CC1ppip_XSec_1DQ2_nu.h" +// FNAL CC1ppim +#include "FNAL_CC1ppim_XSec_1DEnu_antinu.h" +#endif + +#ifndef __NO_BEBC__ +// BEBC CCQE +#include "BEBC_CCQE_XSec_1DQ2_nu.h" +// BEBC CC1ppip +#include "BEBC_CC1ppip_XSec_1DEnu_nu.h" +#include "BEBC_CC1ppip_XSec_1DQ2_nu.h" +// BEBC CC1npip +#include "BEBC_CC1npip_XSec_1DEnu_nu.h" +#include "BEBC_CC1npip_XSec_1DQ2_nu.h" +// BEBC CC1pi0 +#include "BEBC_CC1pi0_XSec_1DEnu_nu.h" +#include "BEBC_CC1pi0_XSec_1DQ2_nu.h" +// BEBC CC1npim +#include "BEBC_CC1npim_XSec_1DEnu_antinu.h" +#include "BEBC_CC1npim_XSec_1DQ2_antinu.h" +// BEBC CC1ppim +#include "BEBC_CC1ppim_XSec_1DEnu_antinu.h" +#include "BEBC_CC1ppim_XSec_1DQ2_antinu.h" +#endif + +#ifndef __NO_GGM__ +// GGM CC1ppip +#include "GGM_CC1ppip_Evt_1DQ2_nu.h" +#include "GGM_CC1ppip_XSec_1DEnu_nu.h" +#endif + +#ifndef __NO_MiniBooNE__ +// MiniBooNE CCQE +#include "MiniBooNE_CCQE_XSec_1DQ2_antinu.h" +#include "MiniBooNE_CCQE_XSec_1DQ2_nu.h" +#include "MiniBooNE_CCQE_XSec_2DTcos_antinu.h" +#include "MiniBooNE_CCQE_XSec_2DTcos_antinu.h" +#include "MiniBooNE_CCQE_XSec_2DTcos_nu.h" + +// MiniBooNE CC1pi+ 1D +#include "MiniBooNE_CC1pip_XSec_1DEnu_nu.h" +#include "MiniBooNE_CC1pip_XSec_1DQ2_nu.h" +#include "MiniBooNE_CC1pip_XSec_1DTpi_nu.h" +#include "MiniBooNE_CC1pip_XSec_1DTu_nu.h" +// MiniBooNE CC1pi+ 2D +#include "MiniBooNE_CC1pip_XSec_2DQ2Enu_nu.h" +#include "MiniBooNE_CC1pip_XSec_2DTpiCospi_nu.h" +#include "MiniBooNE_CC1pip_XSec_2DTpiEnu_nu.h" +#include "MiniBooNE_CC1pip_XSec_2DTuCosmu_nu.h" +#include "MiniBooNE_CC1pip_XSec_2DTuEnu_nu.h" + +// MiniBooNE CC1pi0 +#include "MiniBooNE_CC1pi0_XSec_1DEnu_nu.h" +#include "MiniBooNE_CC1pi0_XSec_1DQ2_nu.h" +#include "MiniBooNE_CC1pi0_XSec_1DTu_nu.h" +#include "MiniBooNE_CC1pi0_XSec_1Dcosmu_nu.h" +#include "MiniBooNE_CC1pi0_XSec_1Dcospi0_nu.h" +#include "MiniBooNE_CC1pi0_XSec_1Dppi0_nu.h" +#include "MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu.h" +#include "MiniBooNE_NC1pi0_XSec_1Dcospi0_nu.h" +#include "MiniBooNE_NC1pi0_XSec_1Dppi0_antinu.h" +#include "MiniBooNE_NC1pi0_XSec_1Dppi0_nu.h" + +// MiniBooNE NC1pi0 +//#include "MiniBooNE_NCpi0_XSec_1Dppi0_nu.h" + +// MiniBooNE NCEL +#include "MiniBooNE_NCEL_XSec_Treco_nu.h" +#endif + +#ifndef __NO_MINERvA__ +// MINERvA CCQE +#include "MINERvA_CCQE_XSec_1DQ2_antinu.h" +#include "MINERvA_CCQE_XSec_1DQ2_joint.h" +#include "MINERvA_CCQE_XSec_1DQ2_nu.h" + +// MINERvA CC0pi +#include "MINERvA_CC0pi_XSec_1DEe_nue.h" +#include "MINERvA_CC0pi_XSec_1DQ2_nu_proton.h" +#include "MINERvA_CC0pi_XSec_1DQ2_nue.h" +#include "MINERvA_CC0pi_XSec_1DThetae_nue.h" + +// MINERvA CC1pi+ +#include "MINERvA_CC1pip_XSec_1DTpi_20deg_nu.h" +#include "MINERvA_CC1pip_XSec_1DTpi_nu.h" +#include "MINERvA_CC1pip_XSec_1Dth_20deg_nu.h" +#include "MINERvA_CC1pip_XSec_1Dth_nu.h" + +// MINERvA CCNpi+ +#include "MINERvA_CCNpip_XSec_1DEnu_nu.h" +#include "MINERvA_CCNpip_XSec_1DQ2_nu.h" +#include "MINERvA_CCNpip_XSec_1DTpi_nu.h" +#include "MINERvA_CCNpip_XSec_1Dpmu_nu.h" +#include "MINERvA_CCNpip_XSec_1Dth_nu.h" +#include "MINERvA_CCNpip_XSec_1Dthmu_nu.h" + +// MINERvA CC1pi0 +#include "MINERvA_CC1pi0_XSec_1DEnu_antinu.h" +#include "MINERvA_CC1pi0_XSec_1DQ2_antinu.h" +#include "MINERvA_CC1pi0_XSec_1DTpi0_antinu.h" +#include "MINERvA_CC1pi0_XSec_1Dpmu_antinu.h" +#include "MINERvA_CC1pi0_XSec_1Dppi0_antinu.h" +#include "MINERvA_CC1pi0_XSec_1Dth_antinu.h" +#include "MINERvA_CC1pi0_XSec_1Dthmu_antinu.h" + +// MINERvA CCINC +#include "MINERvA_CCinc_XSec_1DEnu_ratio.h" +#include "MINERvA_CCinc_XSec_1Dx_ratio.h" +#include "MINERvA_CCinc_XSec_2DEavq3_nu.h" + +// MINERvA CCDIS +#include "MINERvA_CCDIS_XSec_1DEnu_ratio.h" +#include "MINERvA_CCDIS_XSec_1Dx_ratio.h" + +// MINERvA CCCOH pion +#include "MINERvA_CCCOHPI_XSec_1DEnu_antinu.h" +#include "MINERvA_CCCOHPI_XSec_1DEnu_nu.h" +#include "MINERvA_CCCOHPI_XSec_1DEpi_antinu.h" +#include "MINERvA_CCCOHPI_XSec_1DQ2_antinu.h" + +#include "MINERvA_CCCOHPI_XSec_1DEpi_nu.h" +#include "MINERvA_CCCOHPI_XSec_1Dth_antinu.h" +#include "MINERvA_CCCOHPI_XSec_1Dth_nu.h" +#include "MINERvA_CCCOHPI_XSec_1DQ2_nu.h" + +#include "MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.h" +#include "MINERvA_CC0pi_XSec_1DQ2_Tgt_nu.h" +#include "MINERvA_CC1pip_XSec_1D_2017Update.h" + + + +#endif + +#ifndef __NO_T2K__ +// T2K CC0pi +#include "T2K_CC0pi_XSec_2DPcos_nu.h" + +// T2K CC1pi+ on CH +#include "T2K_CC1pip_CH_XSec_1DQ2_nu.h" +#include "T2K_CC1pip_CH_XSec_1DWrec_nu.h" +#include "T2K_CC1pip_CH_XSec_1Dpmu_nu.h" +#include "T2K_CC1pip_CH_XSec_1Dppi_nu.h" +#include "T2K_CC1pip_CH_XSec_1Dq3_nu.h" +#include "T2K_CC1pip_CH_XSec_1Dthmupi_nu.h" +#include "T2K_CC1pip_CH_XSec_1Dthpi_nu.h" +#include "T2K_CC1pip_CH_XSec_1Dthq3pi_nu.h" + +// T2K CC1pi+ on H2O +#include "T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.h" +#include "T2K_CC1pip_H2O_XSec_1DEnuMB_nu.h" +#include "T2K_CC1pip_H2O_XSec_1Dcosmu_nu.h" +#include "T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.h" +#include "T2K_CC1pip_H2O_XSec_1Dcospi_nu.h" +#include "T2K_CC1pip_H2O_XSec_1Dpmu_nu.h" +#include "T2K_CC1pip_H2O_XSec_1Dppi_nu.h" + +// T2K STV CC0pi +#include "T2K_CC0pinp_STV_XSec_1Ddpt_nu.h" + +#include "T2K_CC0pi_XSec_2DPcos_nu_nonuniform.h" + +#endif + +#ifndef __NO_SciBooNE__ + +// SciBooNE COH studies +#include "SciBooNE_CCCOH_1TRK_1DQ2_nu.h" +#include "SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu.h" +#include "SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu.h" +#include "SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu.h" +#include "SciBooNE_CCCOH_MuPiVA_1DQ2_nu.h" +#include "SciBooNE_CCCOH_MuPr_1DQ2_nu.h" +#include "SciBooNE_CCCOH_STOPFINAL_1DQ2_nu.h" +#include "SciBooNE_CCCOH_STOP_NTrks_nu.h" +#endif + +#ifndef __NO_K2K__ +// K2K NC1pi0 +#include "K2K_NC1pi0_Evt_1Dppi0_nu.h" +#endif + +// MC Studies +#include "ExpMultDist_CCQE_XSec_1DVar_FakeStudy.h" +#include "ExpMultDist_CCQE_XSec_2DVar_FakeStudy.h" +#include "MCStudy_CCQEHistograms.h" + +#include "GenericFlux_Tester.h" +#include "GenericFlux_Vectors.h" + +#include "ElectronFlux_FlatTree.h" +#include "ElectronScattering_DurhamData.h" +#include "MCStudy_KaonPreSelection.h" +#include "MCStudy_MuonValidation.h" + +#include "OfficialNIWGPlots.h" +#include "T2K2017_FakeData.h" + +#include "Simple_Osc.h" + +#include "FitWeight.h" + +#include "NuisConfig.h" +#include "NuisKey.h" + + + //! Functions to make it easier for samples to be created and handled. namespace SampleUtils { //! Create a given sample given its name, file, type, fakdata(fkdt) file and the //! current rw engine and push it back into the list fChain. MeasurementBase* CreateSample(std::string name, std::string file, std::string type, std::string fkdt, FitWeight* rw) { nuiskey samplekey = Config::CreateKey("sample"); samplekey.AddS("name", name); samplekey.AddS("input", file); samplekey.AddS("type", type); return CreateSample(samplekey); } MeasurementBase* CreateSample(nuiskey samplekey) { FitWeight* rw = FitBase::GetRW(); std::string name = samplekey.GetS("name"); std::string file = samplekey.GetS("input"); std::string type = samplekey.GetS("type"); std::string fkdt = ""; /* ANL CCQE Samples */ #ifndef __NO_ANL__ if (!name.compare("ANL_CCQE_XSec_1DEnu_nu") || !name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD26") || !name.compare("ANL_CCQE_XSec_1DEnu_nu_PRL31") || !name.compare("ANL_CCQE_XSec_1DEnu_nu_PRD16")) { return (new ANL_CCQE_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CCQE_Evt_1DQ2_nu") || !name.compare("ANL_CCQE_Evt_1DQ2_nu_PRL31") || !name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD26") || !name.compare("ANL_CCQE_Evt_1DQ2_nu_PRD16")) { return (new ANL_CCQE_Evt_1DQ2_nu(samplekey)); /* ANL CC1ppip samples */ } else if (!name.compare("ANL_CC1ppip_XSec_1DEnu_nu") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W14Cut") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_Uncorr") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W14Cut_Uncorr") || !name.compare("ANL_CC1ppip_XSec_1DEnu_nu_W16Cut_Uncorr")) { return (new ANL_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_XSec_1DQ2_nu")) { return (new ANL_CC1ppip_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DQ2_nu") || !name.compare("ANL_CC1ppip_Evt_1DQ2_nu_W14Cut")) { return (new ANL_CC1ppip_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1Dppi_nu")) { return (new ANL_CC1ppip_Evt_1Dppi_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1Dthpr_nu")) { return (new ANL_CC1ppip_Evt_1Dthpr_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DcosmuStar_nu")) { return (new ANL_CC1ppip_Evt_1DcosmuStar_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1DcosthAdler_nu")) { return (new ANL_CC1ppip_Evt_1DcosthAdler_nu(samplekey)); } else if (!name.compare("ANL_CC1ppip_Evt_1Dphi_nu")) { return (new ANL_CC1ppip_Evt_1Dphi_nu(samplekey)); /* ANL CC1npip sample */ } else if (!name.compare("ANL_CC1npip_XSec_1DEnu_nu") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_W14Cut") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_Uncorr") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_W14Cut_Uncorr") || !name.compare("ANL_CC1npip_XSec_1DEnu_nu_W16Cut_Uncorr")) { return (new ANL_CC1npip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1DQ2_nu") || !name.compare("ANL_CC1npip_Evt_1DQ2_nu_W14Cut")) { return (new ANL_CC1npip_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1Dppi_nu")) { return (new ANL_CC1npip_Evt_1Dppi_nu(samplekey)); } else if (!name.compare("ANL_CC1npip_Evt_1DcosmuStar_nu")) { return (new ANL_CC1npip_Evt_1DcosmuStar_nu(samplekey)); /* ANL CC1pi0 sample */ } else if (!name.compare("ANL_CC1pi0_XSec_1DEnu_nu") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W14Cut") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_Uncorr") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W14Cut_Uncorr") || !name.compare("ANL_CC1pi0_XSec_1DEnu_nu_W16Cut_Uncorr")) { return (new ANL_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC1pi0_Evt_1DQ2_nu") || !name.compare("ANL_CC1pi0_Evt_1DQ2_nu_W14Cut")) { return (new ANL_CC1pi0_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("ANL_CC1pi0_Evt_1DcosmuStar_nu")) { return (new ANL_CC1pi0_Evt_1DcosmuStar_nu(samplekey)); /* ANL NC1npip sample */ } else if (!name.compare("ANL_NC1npip_Evt_1Dppi_nu")) { return (new ANL_NC1npip_Evt_1Dppi_nu(samplekey)); /* ANL NC1ppim sample */ } else if (!name.compare("ANL_NC1ppim_XSec_1DEnu_nu")) { return (new ANL_NC1ppim_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_NC1ppim_Evt_1DcosmuStar_nu")) { return (new ANL_NC1ppim_Evt_1DcosmuStar_nu(samplekey)); /* ANL CC2pi sample */ } else if (!name.compare("ANL_CC2pi_1pim1pip_XSec_1DEnu_nu")) { return (new ANL_CC2pi_1pim1pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dppip_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dppip_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dppim_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dppim_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu")) { return (new ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_XSec_1DEnu_nu")) { return (new ANL_CC2pi_1pip1pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu")) { return (new ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu")) { return (new ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu(samplekey)); } else if (!name.compare("ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu")) { return (new ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu(samplekey)); /* ArgoNeut Samples */ } else #endif #ifndef __NO_ArgoNeuT__ if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_antinu")) { return (new ArgoNeuT_CCInc_XSec_1Dpmu_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dpmu_nu")) { return (new ArgoNeuT_CCInc_XSec_1Dpmu_nu(samplekey)); } else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_antinu")) { return (new ArgoNeuT_CCInc_XSec_1Dthetamu_antinu(samplekey)); } else if (!name.compare("ArgoNeuT_CCInc_XSec_1Dthetamu_nu")) { return (new ArgoNeuT_CCInc_XSec_1Dthetamu_nu(samplekey)); /* BNL Samples */ } else #endif #ifndef __NO_BNL__ if (!name.compare("BNL_CCQE_XSec_1DEnu_nu")) { return (new BNL_CCQE_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CCQE_Evt_1DQ2_nu")) { return (new BNL_CCQE_Evt_1DQ2_nu(samplekey)); /* BNL CC1ppip samples */ } else if (!name.compare("BNL_CC1ppip_XSec_1DEnu_nu") || !name.compare("BNL_CC1ppip_XSec_1DEnu_nu_Uncorr") || !name.compare("BNL_CC1ppip_XSec_1DEnu_nu_W14Cut") || !name.compare("BNL_CC1ppip_XSec_1DEnu_nu_W14Cut_Uncorr")) { return (new BNL_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CC1ppip_Evt_1DQ2_nu") || !name.compare("BNL_CC1ppip_Evt_1DQ2_nu_W14Cut")) { return (new BNL_CC1ppip_Evt_1DQ2_nu(samplekey)); } else if (!name.compare("BNL_CC1ppip_Evt_1DcosthAdler_nu")) { return (new BNL_CC1ppip_Evt_1DcosthAdler_nu(samplekey)); } else if (!name.compare("BNL_CC1ppip_Evt_1Dphi_nu")) { return (new BNL_CC1ppip_Evt_1Dphi_nu(samplekey)); /* BNL CC1npip samples */ } else if (!name.compare("BNL_CC1npip_XSec_1DEnu_nu") || !name.compare("BNL_CC1npip_XSec_1DEnu_nu_Uncorr")) { return (new BNL_CC1npip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CC1npip_Evt_1DQ2_nu")) { return (new BNL_CC1npip_Evt_1DQ2_nu(samplekey)); /* BNL CC1pi0 samples */ } else if (!name.compare("BNL_CC1pi0_XSec_1DEnu_nu")) { return (new BNL_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BNL_CC1pi0_Evt_1DQ2_nu")) { return (new BNL_CC1pi0_Evt_1DQ2_nu(samplekey)); /* FNAL Samples */ } else #endif #ifndef __NO_FNAL__ if (!name.compare("FNAL_CCQE_Evt_1DQ2_nu")) { return (new FNAL_CCQE_Evt_1DQ2_nu(samplekey)); /* FNAL CC1ppip */ } else if (!name.compare("FNAL_CC1ppip_XSec_1DEnu_nu")) { return (new FNAL_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("FNAL_CC1ppip_XSec_1DQ2_nu")) { return (new FNAL_CC1ppip_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("FNAL_CC1ppip_Evt_1DQ2_nu")) { return (new FNAL_CC1ppip_Evt_1DQ2_nu(samplekey)); /* FNAL CC1ppim */ } else if (!name.compare("FNAL_CC1ppim_XSec_1DEnu_antinu")) { return (new FNAL_CC1ppim_XSec_1DEnu_antinu(samplekey)); /* BEBC Samples */ } else #endif #ifndef __NO_BEBC__ if (!name.compare("BEBC_CCQE_XSec_1DQ2_nu")) { return (new BEBC_CCQE_XSec_1DQ2_nu(samplekey)); /* BEBC CC1ppip samples */ } else if (!name.compare("BEBC_CC1ppip_XSec_1DEnu_nu")) { return (new BEBC_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BEBC_CC1ppip_XSec_1DQ2_nu")) { return (new BEBC_CC1ppip_XSec_1DQ2_nu(samplekey)); /* BEBC CC1npip samples */ } else if (!name.compare("BEBC_CC1npip_XSec_1DEnu_nu")) { return (new BEBC_CC1npip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BEBC_CC1npip_XSec_1DQ2_nu")) { return (new BEBC_CC1npip_XSec_1DQ2_nu(samplekey)); /* BEBC CC1pi0 samples */ } else if (!name.compare("BEBC_CC1pi0_XSec_1DEnu_nu")) { return (new BEBC_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("BEBC_CC1pi0_XSec_1DQ2_nu")) { return (new BEBC_CC1pi0_XSec_1DQ2_nu(samplekey)); /* BEBC CC1npim samples */ } else if (!name.compare("BEBC_CC1npim_XSec_1DEnu_antinu")) { return (new BEBC_CC1npim_XSec_1DEnu_antinu(samplekey)); } else if (!name.compare("BEBC_CC1npim_XSec_1DQ2_antinu")) { return (new BEBC_CC1npim_XSec_1DQ2_antinu(samplekey)); /* BEBC CC1ppim samples */ } else if (!name.compare("BEBC_CC1ppim_XSec_1DEnu_antinu")) { return (new BEBC_CC1ppim_XSec_1DEnu_antinu(samplekey)); } else if (!name.compare("BEBC_CC1ppim_XSec_1DQ2_antinu")) { return (new BEBC_CC1ppim_XSec_1DQ2_antinu(samplekey)); /* GGM CC1ppip samples */ } else #endif #ifndef __NO_GGM__ if (!name.compare("GGM_CC1ppip_XSec_1DEnu_nu")) { return (new GGM_CC1ppip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("GGM_CC1ppip_Evt_1DQ2_nu")) { return (new GGM_CC1ppip_Evt_1DQ2_nu(samplekey)); /* MiniBooNE Samples */ /* CCQE */ } else #endif #ifndef __NO_MiniBooNE__ if (!name.compare("MiniBooNE_CCQE_XSec_1DQ2_nu") || !name.compare("MiniBooNE_CCQELike_XSec_1DQ2_nu")) { return (new MiniBooNE_CCQE_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MiniBooNE_CCQE_XSec_1DQ2_antinu") || !name.compare("MiniBooNE_CCQELike_XSec_1DQ2_antinu") || !name.compare("MiniBooNE_CCQE_CTarg_XSec_1DQ2_antinu")) { return (new MiniBooNE_CCQE_XSec_1DQ2_antinu(samplekey)); } else if (!name.compare("MiniBooNE_CCQE_XSec_2DTcos_nu") || !name.compare("MiniBooNE_CCQELike_XSec_2DTcos_nu")) { return (new MiniBooNE_CCQE_XSec_2DTcos_nu(samplekey)); } else if (!name.compare("MiniBooNE_CCQE_XSec_2DTcos_antinu") || !name.compare("MiniBooNE_CCQELike_XSec_2DTcos_antinu")) { return (new MiniBooNE_CCQE_XSec_2DTcos_antinu(samplekey)); /* MiniBooNE CC1pi+ */ // 1D } else if (!name.compare("MiniBooNE_CC1pip_XSec_1DEnu_nu")) { return (new MiniBooNE_CC1pip_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_1DQ2_nu")) { return (new MiniBooNE_CC1pip_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_1DTpi_nu")) { return (new MiniBooNE_CC1pip_XSec_1DTpi_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_1DTu_nu")) { return (new MiniBooNE_CC1pip_XSec_1DTu_nu(samplekey)); // 2D } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DQ2Enu_nu")) { return (new MiniBooNE_CC1pip_XSec_2DQ2Enu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTpiCospi_nu")) { return (new MiniBooNE_CC1pip_XSec_2DTpiCospi_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTpiEnu_nu")) { return (new MiniBooNE_CC1pip_XSec_2DTpiEnu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTuCosmu_nu")) { return (new MiniBooNE_CC1pip_XSec_2DTuCosmu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pip_XSec_2DTuEnu_nu")) { return (new MiniBooNE_CC1pip_XSec_2DTuEnu_nu(samplekey)); /* MiniBooNE CC1pi0 */ } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DEnu_nu")) { return (new MiniBooNE_CC1pi0_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DQ2_nu")) { return (new MiniBooNE_CC1pi0_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1DTu_nu")) { return (new MiniBooNE_CC1pi0_XSec_1DTu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dcosmu_nu")) { return (new MiniBooNE_CC1pi0_XSec_1Dcosmu_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dcospi0_nu")) { return (new MiniBooNE_CC1pi0_XSec_1Dcospi0_nu(samplekey)); } else if (!name.compare("MiniBooNE_CC1pi0_XSec_1Dppi0_nu")) { return (new MiniBooNE_CC1pi0_XSec_1Dppi0_nu(samplekey)); } else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu") || !name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_rhc")) { return (new MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu(samplekey)); } else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_nu") || !name.compare("MiniBooNE_NC1pi0_XSec_1Dcospi0_fhc")) { return (new MiniBooNE_NC1pi0_XSec_1Dcospi0_nu(samplekey)); } else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_antinu") || !name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_rhc")) { return (new MiniBooNE_NC1pi0_XSec_1Dppi0_antinu(samplekey)); } else if (!name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_nu") || !name.compare("MiniBooNE_NC1pi0_XSec_1Dppi0_fhc")) { return (new MiniBooNE_NC1pi0_XSec_1Dppi0_nu(samplekey)); /* MiniBooNE NCEL */ } else if (!name.compare("MiniBooNE_NCEL_XSec_Treco_nu")) { return (new MiniBooNE_NCEL_XSec_Treco_nu(samplekey)); /* MINERvA Samples */ } else #endif #ifndef __NO_MINERvA__ if (!name.compare("MINERvA_CCQE_XSec_1DQ2_nu") || !name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg") || !name.compare("MINERvA_CCQE_XSec_1DQ2_nu_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_nu_20deg_oldflux")) { return (new MINERvA_CCQE_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_antinu") || !name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg") || !name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_antinu_20deg_oldflux")) { return (new MINERvA_CCQE_XSec_1DQ2_antinu(samplekey)); } else if (!name.compare("MINERvA_CCQE_XSec_1DQ2_joint_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg_oldflux") || !name.compare("MINERvA_CCQE_XSec_1DQ2_joint") || !name.compare("MINERvA_CCQE_XSec_1DQ2_joint_20deg")) { return (new MINERvA_CCQE_XSec_1DQ2_joint(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DEe_nue")) { return (new MINERvA_CC0pi_XSec_1DEe_nue(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nue")) { return (new MINERvA_CC0pi_XSec_1DQ2_nue(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DThetae_nue")) { return (new MINERvA_CC0pi_XSec_1DThetae_nue(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_nu_proton")) { return (new MINERvA_CC0pi_XSec_1DQ2_nu_proton(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtC_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtCH_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtFe_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtPb_nu")) { return (new MINERvA_CC0pi_XSec_1DQ2_Tgt_nu(samplekey)); } else if (!name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioC_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioFe_nu") || !name.compare("MINERvA_CC0pi_XSec_1DQ2_TgtRatioPb_nu")) { return (new MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu(samplekey)); /* CC1pi+ */ // DONE } else if (!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu") || !name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_20deg") || !name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_fluxcorr") || !name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_20deg_fluxcorr")) { return (new MINERvA_CC1pip_XSec_1DTpi_nu(samplekey)); // DONE } else if (!name.compare("MINERvA_CC1pip_XSec_1Dth_nu") || !name.compare("MINERvA_CC1pip_XSec_1Dth_nu_20deg") || !name.compare("MINERvA_CC1pip_XSec_1Dth_nu_fluxcorr") || !name.compare("MINERvA_CC1pip_XSec_1Dth_nu_20deg_fluxcorr")) { return (new MINERvA_CC1pip_XSec_1Dth_nu(samplekey)); } else if (!name.compare("MINERvA_CC1pip_XSec_1DTpi_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1Dth_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1Dpmu_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1Dthmu_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1DQ2_nu_2017") || !name.compare("MINERvA_CC1pip_XSec_1DEnu_nu_2017")) { return (new MINERvA_CC1pip_XSec_1D_2017Update(samplekey)); /* CCNpi+ */ } else if (!name.compare("MINERvA_CCNpip_XSec_1Dth_nu") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2016") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_20deg") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_fluxcorr") || !name.compare("MINERvA_CCNpip_XSec_1Dth_nu_2015_20deg_fluxcorr")) { return (new MINERvA_CCNpip_XSec_1Dth_nu(samplekey)); } else if (!name.compare("MINERvA_CCNpip_XSec_1DTpi_nu") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2016") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_20deg") || !name.compare("MINERvA_CCNpip_XSec_1DTpi_nu_2015_fluxcorr") || !name.compare( "MINERvA_CCNpip_XSec_1DTpi_nu_2015_20deg_fluxcorr")) { return (new MINERvA_CCNpip_XSec_1DTpi_nu(samplekey)); // Done } else if (!name.compare("MINERvA_CCNpip_XSec_1Dthmu_nu")) { return (new MINERvA_CCNpip_XSec_1Dthmu_nu(samplekey)); // Done } else if (!name.compare("MINERvA_CCNpip_XSec_1Dpmu_nu")) { return (new MINERvA_CCNpip_XSec_1Dpmu_nu(samplekey)); // Done } else if (!name.compare("MINERvA_CCNpip_XSec_1DQ2_nu")) { return (new MINERvA_CCNpip_XSec_1DQ2_nu(samplekey)); // Done } else if (!name.compare("MINERvA_CCNpip_XSec_1DEnu_nu")) { return (new MINERvA_CCNpip_XSec_1DEnu_nu(samplekey)); /* CC1pi0 */ // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2015") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2016") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_fluxcorr") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2015_fluxcorr") || !name.compare("MINERvA_CC1pi0_XSec_1Dth_antinu_2016_fluxcorr")) { return (new MINERvA_CC1pi0_XSec_1Dth_antinu(samplekey)); } else if (!name.compare("MINERvA_CC1pi0_XSec_1Dppi0_antinu") || !name.compare("MINERvA_CC1pi0_XSec_1Dppi0_antinu_fluxcorr")) { return (new MINERvA_CC1pi0_XSec_1Dppi0_antinu(samplekey)); } else if (!name.compare("MINERvA_CC1pi0_XSec_1DTpi0_antinu")) { return (new MINERvA_CC1pi0_XSec_1DTpi0_antinu(samplekey)); // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1DQ2_antinu")) { return (new MINERvA_CC1pi0_XSec_1DQ2_antinu(samplekey)); // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1Dthmu_antinu")) { return (new MINERvA_CC1pi0_XSec_1Dthmu_antinu(samplekey)); // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1Dpmu_antinu")) { return (new MINERvA_CC1pi0_XSec_1Dpmu_antinu(samplekey)); // Done } else if (!name.compare("MINERvA_CC1pi0_XSec_1DEnu_antinu")) { return (new MINERvA_CC1pi0_XSec_1DEnu_antinu(samplekey)); /* CCINC */ } else if (!name.compare("MINERvA_CCinc_XSec_2DEavq3_nu")) { return (new MINERvA_CCinc_XSec_2DEavq3_nu(samplekey)); } else if (!name.compare("MINERvA_CCinc_XSec_1Dx_ratio_C12_CH") || !name.compare("MINERvA_CCinc_XSec_1Dx_ratio_Fe56_CH") || !name.compare("MINERvA_CCinc_XSec_1Dx_ratio_Pb208_CH")) { return (new MINERvA_CCinc_XSec_1Dx_ratio(samplekey)); } else if (!name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_C12_CH") || !name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_Fe56_CH") || !name.compare("MINERvA_CCinc_XSec_1DEnu_ratio_Pb208_CH")) { return (new MINERvA_CCinc_XSec_1DEnu_ratio(samplekey)); /* CCDIS */ } else if (!name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_C12_CH") || !name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_Fe56_CH") || !name.compare("MINERvA_CCDIS_XSec_1Dx_ratio_Pb208_CH")) { return (new MINERvA_CCDIS_XSec_1Dx_ratio(samplekey)); } else if (!name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_C12_CH") || !name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_Fe56_CH") || !name.compare("MINERvA_CCDIS_XSec_1DEnu_ratio_Pb208_CH")) { return (new MINERvA_CCDIS_XSec_1DEnu_ratio(samplekey)); /* CC-COH */ } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_nu")) { return (new MINERvA_CCCOHPI_XSec_1DEnu_nu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_nu")) { return (new MINERvA_CCCOHPI_XSec_1DEpi_nu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_nu")) { return (new MINERvA_CCCOHPI_XSec_1Dth_nu(samplekey)); + } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_nu")) { + return (new MINERvA_CCCOHPI_XSec_1DQ2_nu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEnu_antinu")) { return (new MINERvA_CCCOHPI_XSec_1DEnu_antinu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DEpi_antinu")) { return (new MINERvA_CCCOHPI_XSec_1DEpi_antinu(samplekey)); } else if (!name.compare("MINERvA_CCCOHPI_XSec_1Dth_antinu")) { return (new MINERvA_CCCOHPI_XSec_1Dth_antinu(samplekey)); + } else if (!name.compare("MINERvA_CCCOHPI_XSec_1DQ2_antinu")) { + return (new MINERvA_CCCOHPI_XSec_1DQ2_antinu(samplekey)); + /* T2K Samples */ } else #endif #ifndef __NO_T2K__ if (!name.compare("T2K_CC0pi_XSec_2DPcos_nu") || !name.compare("T2K_CC0pi_XSec_2DPcos_nu_I") || !name.compare("T2K_CC0pi_XSec_2DPcos_nu_II")) { return (new T2K_CC0pi_XSec_2DPcos_nu(samplekey)); } else if (!name.compare("T2K_CC0pi_XSec_2DPcos_nu_nonuniform")) { return (new T2K_CC0pi_XSec_2DPcos_nu_nonuniform(samplekey)); /* T2K CC1pi+ CH samples */ // Comment these out for now because we don't have the proper data /* } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dpmu_nu")) { return (new T2K_CC1pip_CH_XSec_1Dpmu_nu(file, rw, type, fkdt)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dppi_nu")) { return (new T2K_CC1pip_CH_XSec_1Dppi_nu(file, rw, type, fkdt)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1DQ2_nu")) { return (new T2K_CC1pip_CH_XSec_1DQ2_nu(file, rw, type, fkdt)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dq3_nu")) { return (new T2K_CC1pip_CH_XSec_1Dq3_nu(file, rw, type, fkdt)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthmupi_nu")) { return (new T2K_CC1pip_CH_XSec_1Dthmupi_nu(file, rw, type, fkdt)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthpi_nu")) { return (new T2K_CC1pip_CH_XSec_1Dthpi_nu(file, rw, type, fkdt)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1Dthq3pi_nu")) { return (new T2K_CC1pip_CH_XSec_1Dthq3pi_nu(file, rw, type, fkdt)); } else if (!name.compare("T2K_CC1pip_CH_XSec_1DWrec_nu")) { return (new T2K_CC1pip_CH_XSec_1DWrec_nu(file, rw, type, fkdt)); */ /* T2K CC1pi+ H2O samples */ } else if (!name.compare("T2K_CC1pip_H2O_XSec_1DEnuDelta_nu")) { return (new T2K_CC1pip_H2O_XSec_1DEnuDelta_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1DEnuMB_nu")) { return (new T2K_CC1pip_H2O_XSec_1DEnuMB_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcosmu_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dcosmu_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcosmupi_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dcosmupi_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dcospi_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dcospi_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dpmu_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dpmu_nu(samplekey)); } else if (!name.compare("T2K_CC1pip_H2O_XSec_1Dppi_nu")) { return (new T2K_CC1pip_H2O_XSec_1Dppi_nu(samplekey)); /* T2K CC0pi + np CH samples */ } else if (!name.compare("T2K_CC0pinp_STV_XSec_1Ddpt_nu")) { return (new T2K_CC0pinp_STV_XSec_1Ddpt_nu(samplekey)); // SciBooNE COH studies } else #endif #ifndef __NO_SciBooNE__ if (!name.compare("SciBooNE_CCCOH_STOP_NTrks_nu")) { return (new SciBooNE_CCCOH_STOP_NTrks_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_1TRK_1DQ2_nu")) { return (new SciBooNE_CCCOH_1TRK_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPr_1DQ2_nu")) { return (new SciBooNE_CCCOH_MuPr_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiVA_1DQ2_nu")) { return (new SciBooNE_CCCOH_MuPiVA_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu")) { return (new SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu(samplekey)); } else if (!name.compare("SciBooNE_CCCOH_STOPFINAL_1DQ2_nu")) { return (new SciBooNE_CCCOH_STOPFINAL_1DQ2_nu(samplekey)); /* K2K Samples */ /* NC1pi0 */ } else #endif #ifndef __NO_K2K__ if (!name.compare("K2K_NC1pi0_Evt_1Dppi0_nu")) { return (new K2K_NC1pi0_Evt_1Dppi0_nu(samplekey)); /* Fake Studies */ } else #endif if (name.find("ExpMultDist_CCQE_XSec_1D") != std::string::npos && name.find("_FakeStudy") != std::string::npos) { return ( new ExpMultDist_CCQE_XSec_1DVar_FakeStudy(name, file, rw, type, fkdt)); } else if (name.find("ExpMultDist_CCQE_XSec_2D") != std::string::npos && name.find("_FakeStudy") != std::string::npos) { return ( new ExpMultDist_CCQE_XSec_2DVar_FakeStudy(name, file, rw, type, fkdt)); } else if (name.find("GenericFlux_") != std::string::npos) { return (new GenericFlux_Tester(name, file, rw, type, fkdt)); } else if (name.find("GenericVectors_") != std::string::npos) { return (new GenericFlux_Vectors(name, file, rw, type, fkdt)); } else if (!name.compare("T2K2017_FakeData")) { return (new T2K2017_FakeData(samplekey)); } else if (!name.compare("MCStudy_CCQE")) { return (new MCStudy_CCQEHistograms(name, file, rw, type, fkdt)); } else if (!name.compare("ElectronFlux_FlatTree")) { return (new ElectronFlux_FlatTree(name, file, rw, type, fkdt)); } else if (name.find("ElectronData_") != std::string::npos) { return new ElectronScattering_DurhamData(samplekey); } else if (name.find("MuonValidation_") != std::string::npos) { return (new MCStudy_MuonValidation(name, file, rw, type, fkdt)); } else if (!name.compare("NIWGOfficialPlots")) { return (new OfficialNIWGPlots(samplekey)); }else if (!name.compare("Simple_Osc")) { return (new Simple_Osc(samplekey)); } else { ERR(FTL) << "Error: No such sample: " << name << std::endl; exit(-1); return NULL; } // Return NULL if no sample loaded. return NULL; } } diff --git a/src/FCN/SampleList.h b/src/FCN/SampleList.h index fbf54a8..dbf0356 100644 --- a/src/FCN/SampleList.h +++ b/src/FCN/SampleList.h @@ -1,303 +1,29 @@ #ifndef _SAMPLE_LIST_H_ #define _SAMPLE_LIST_H_ /*! * \addtogroup FCN * @{ */ +#include -#ifndef __NO_ANL__ -#include "ANL_CCQE_Evt_1DQ2_nu.h" -#include "ANL_CCQE_XSec_1DEnu_nu.h" - -// ANL CC1ppip -#include "ANL_CC1ppip_Evt_1DQ2_nu.h" -#include "ANL_CC1ppip_Evt_1DcosmuStar_nu.h" -#include "ANL_CC1ppip_Evt_1DcosmuStar_nu.h" -#include "ANL_CC1ppip_Evt_1DcosthAdler_nu.h" -#include "ANL_CC1ppip_Evt_1Dphi_nu.h" -#include "ANL_CC1ppip_Evt_1Dppi_nu.h" -#include "ANL_CC1ppip_Evt_1Dthpr_nu.h" -#include "ANL_CC1ppip_XSec_1DEnu_nu.h" -#include "ANL_CC1ppip_XSec_1DQ2_nu.h" -// ANL CC1npip -#include "ANL_CC1npip_Evt_1DQ2_nu.h" -#include "ANL_CC1npip_Evt_1DcosmuStar_nu.h" -#include "ANL_CC1npip_Evt_1Dppi_nu.h" -#include "ANL_CC1npip_XSec_1DEnu_nu.h" -// ANL CC1pi0 -#include "ANL_CC1pi0_Evt_1DQ2_nu.h" -#include "ANL_CC1pi0_Evt_1DcosmuStar_nu.h" -#include "ANL_CC1pi0_XSec_1DEnu_nu.h" -// ANL NC1npip (mm, exotic!) -#include "ANL_NC1npip_Evt_1Dppi_nu.h" -// ANL NC1ppim (mm, exotic!) -#include "ANL_NC1ppim_Evt_1DcosmuStar_nu.h" -#include "ANL_NC1ppim_XSec_1DEnu_nu.h" -// ANL CC2pi 1pim1pip (mm, even more exotic!) -#include "ANL_CC2pi_1pim1pip_Evt_1Dpmu_nu.h" -#include "ANL_CC2pi_1pim1pip_Evt_1Dppim_nu.h" -#include "ANL_CC2pi_1pim1pip_Evt_1Dppip_nu.h" -#include "ANL_CC2pi_1pim1pip_Evt_1Dpprot_nu.h" -#include "ANL_CC2pi_1pim1pip_XSec_1DEnu_nu.h" -// ANL CC2pi 1pip1pip (mm, even more exotic!) -#include "ANL_CC2pi_1pip1pip_Evt_1Dpmu_nu.h" -#include "ANL_CC2pi_1pip1pip_Evt_1Dpneut_nu.h" -#include "ANL_CC2pi_1pip1pip_Evt_1DppipHigh_nu.h" -#include "ANL_CC2pi_1pip1pip_Evt_1DppipLow_nu.h" -#include "ANL_CC2pi_1pip1pip_XSec_1DEnu_nu.h" -// ANL CC2pi 1pip1pi0 (mm, even more exotic!) -#include "ANL_CC2pi_1pip1pi0_Evt_1Dpmu_nu.h" -#include "ANL_CC2pi_1pip1pi0_Evt_1Dppi0_nu.h" -#include "ANL_CC2pi_1pip1pi0_Evt_1Dppip_nu.h" -#include "ANL_CC2pi_1pip1pi0_Evt_1Dpprot_nu.h" -#include "ANL_CC2pi_1pip1pi0_XSec_1DEnu_nu.h" -#endif - -#ifndef __NO_ArgoNeuT__ -// ArgoNeuT CC-inclusive -#include "ArgoNeuT_CCInc_XSec_1Dpmu_antinu.h" -#include "ArgoNeuT_CCInc_XSec_1Dpmu_nu.h" -#include "ArgoNeuT_CCInc_XSec_1Dthetamu_antinu.h" -#include "ArgoNeuT_CCInc_XSec_1Dthetamu_nu.h" -#endif - -#ifndef __NO_BNL__ -// BNL CCQE -#include "BNL_CCQE_Evt_1DQ2_nu.h" -#include "BNL_CCQE_XSec_1DEnu_nu.h" -// BNL CC1ppip -#include "BNL_CC1ppip_Evt_1DQ2_nu.h" -#include "BNL_CC1ppip_Evt_1DQ2_nu.h" -#include "BNL_CC1ppip_Evt_1DcosthAdler_nu.h" -#include "BNL_CC1ppip_Evt_1Dphi_nu.h" -#include "BNL_CC1ppip_XSec_1DEnu_nu.h" -// BNL CC1npip -#include "BNL_CC1npip_Evt_1DQ2_nu.h" -#include "BNL_CC1npip_XSec_1DEnu_nu.h" -// BNL CC1pi0 -#include "BNL_CC1pi0_Evt_1DQ2_nu.h" -#include "BNL_CC1pi0_XSec_1DEnu_nu.h" -#endif - -#ifndef __NO_FNAL__ -// FNAL CCQE -#include "FNAL_CCQE_Evt_1DQ2_nu.h" -// FNAL CC1ppip -#include "FNAL_CC1ppip_Evt_1DQ2_nu.h" -#include "FNAL_CC1ppip_XSec_1DEnu_nu.h" -#include "FNAL_CC1ppip_XSec_1DQ2_nu.h" -// FNAL CC1ppim -#include "FNAL_CC1ppim_XSec_1DEnu_antinu.h" -#endif - -#ifndef __NO_BEBC__ -// BEBC CCQE -#include "BEBC_CCQE_XSec_1DQ2_nu.h" -// BEBC CC1ppip -#include "BEBC_CC1ppip_XSec_1DEnu_nu.h" -#include "BEBC_CC1ppip_XSec_1DQ2_nu.h" -// BEBC CC1npip -#include "BEBC_CC1npip_XSec_1DEnu_nu.h" -#include "BEBC_CC1npip_XSec_1DQ2_nu.h" -// BEBC CC1pi0 -#include "BEBC_CC1pi0_XSec_1DEnu_nu.h" -#include "BEBC_CC1pi0_XSec_1DQ2_nu.h" -// BEBC CC1npim -#include "BEBC_CC1npim_XSec_1DEnu_antinu.h" -#include "BEBC_CC1npim_XSec_1DQ2_antinu.h" -// BEBC CC1ppim -#include "BEBC_CC1ppim_XSec_1DEnu_antinu.h" -#include "BEBC_CC1ppim_XSec_1DQ2_antinu.h" -#endif - -#ifndef __NO_GGM__ -// GGM CC1ppip -#include "GGM_CC1ppip_Evt_1DQ2_nu.h" -#include "GGM_CC1ppip_XSec_1DEnu_nu.h" -#endif - -#ifndef __NO_MiniBooNE__ -// MiniBooNE CCQE -#include "MiniBooNE_CCQE_XSec_1DQ2_antinu.h" -#include "MiniBooNE_CCQE_XSec_1DQ2_nu.h" -#include "MiniBooNE_CCQE_XSec_2DTcos_antinu.h" -#include "MiniBooNE_CCQE_XSec_2DTcos_antinu.h" -#include "MiniBooNE_CCQE_XSec_2DTcos_nu.h" - -// MiniBooNE CC1pi+ 1D -#include "MiniBooNE_CC1pip_XSec_1DEnu_nu.h" -#include "MiniBooNE_CC1pip_XSec_1DQ2_nu.h" -#include "MiniBooNE_CC1pip_XSec_1DTpi_nu.h" -#include "MiniBooNE_CC1pip_XSec_1DTu_nu.h" -// MiniBooNE CC1pi+ 2D -#include "MiniBooNE_CC1pip_XSec_2DQ2Enu_nu.h" -#include "MiniBooNE_CC1pip_XSec_2DTpiCospi_nu.h" -#include "MiniBooNE_CC1pip_XSec_2DTpiEnu_nu.h" -#include "MiniBooNE_CC1pip_XSec_2DTuCosmu_nu.h" -#include "MiniBooNE_CC1pip_XSec_2DTuEnu_nu.h" - -// MiniBooNE CC1pi0 -#include "MiniBooNE_CC1pi0_XSec_1DEnu_nu.h" -#include "MiniBooNE_CC1pi0_XSec_1DQ2_nu.h" -#include "MiniBooNE_CC1pi0_XSec_1DTu_nu.h" -#include "MiniBooNE_CC1pi0_XSec_1Dcosmu_nu.h" -#include "MiniBooNE_CC1pi0_XSec_1Dcospi0_nu.h" -#include "MiniBooNE_CC1pi0_XSec_1Dppi0_nu.h" -#include "MiniBooNE_NC1pi0_XSec_1Dcospi0_antinu.h" -#include "MiniBooNE_NC1pi0_XSec_1Dcospi0_nu.h" -#include "MiniBooNE_NC1pi0_XSec_1Dppi0_antinu.h" -#include "MiniBooNE_NC1pi0_XSec_1Dppi0_nu.h" - -// MiniBooNE NC1pi0 -//#include "MiniBooNE_NCpi0_XSec_1Dppi0_nu.h" - -// MiniBooNE NCEL -#include "MiniBooNE_NCEL_XSec_Treco_nu.h" -#endif - -#ifndef __NO_MINERvA__ -// MINERvA CCQE -#include "MINERvA_CCQE_XSec_1DQ2_antinu.h" -#include "MINERvA_CCQE_XSec_1DQ2_joint.h" -#include "MINERvA_CCQE_XSec_1DQ2_nu.h" - -// MINERvA CC0pi -#include "MINERvA_CC0pi_XSec_1DEe_nue.h" -#include "MINERvA_CC0pi_XSec_1DQ2_nu_proton.h" -#include "MINERvA_CC0pi_XSec_1DQ2_nue.h" -#include "MINERvA_CC0pi_XSec_1DThetae_nue.h" - -// MINERvA CC1pi+ -#include "MINERvA_CC1pip_XSec_1DTpi_20deg_nu.h" -#include "MINERvA_CC1pip_XSec_1DTpi_nu.h" -#include "MINERvA_CC1pip_XSec_1Dth_20deg_nu.h" -#include "MINERvA_CC1pip_XSec_1Dth_nu.h" - -// MINERvA CCNpi+ -#include "MINERvA_CCNpip_XSec_1DEnu_nu.h" -#include "MINERvA_CCNpip_XSec_1DQ2_nu.h" -#include "MINERvA_CCNpip_XSec_1DTpi_nu.h" -#include "MINERvA_CCNpip_XSec_1Dpmu_nu.h" -#include "MINERvA_CCNpip_XSec_1Dth_nu.h" -#include "MINERvA_CCNpip_XSec_1Dthmu_nu.h" - -// MINERvA CC1pi0 -#include "MINERvA_CC1pi0_XSec_1DEnu_antinu.h" -#include "MINERvA_CC1pi0_XSec_1DQ2_antinu.h" -#include "MINERvA_CC1pi0_XSec_1DTpi0_antinu.h" -#include "MINERvA_CC1pi0_XSec_1Dpmu_antinu.h" -#include "MINERvA_CC1pi0_XSec_1Dppi0_antinu.h" -#include "MINERvA_CC1pi0_XSec_1Dth_antinu.h" -#include "MINERvA_CC1pi0_XSec_1Dthmu_antinu.h" - -// MINERvA CCINC -#include "MINERvA_CCinc_XSec_1DEnu_ratio.h" -#include "MINERvA_CCinc_XSec_1Dx_ratio.h" -#include "MINERvA_CCinc_XSec_2DEavq3_nu.h" - -// MINERvA CCDIS -#include "MINERvA_CCDIS_XSec_1DEnu_ratio.h" -#include "MINERvA_CCDIS_XSec_1Dx_ratio.h" - -// MINERvA CCCOH pion -#include "MINERvA_CCCOHPI_XSec_1DEnu_antinu.h" -#include "MINERvA_CCCOHPI_XSec_1DEnu_nu.h" -#include "MINERvA_CCCOHPI_XSec_1DEpi_antinu.h" -#include "MINERvA_CCCOHPI_XSec_1DEpi_nu.h" -#include "MINERvA_CCCOHPI_XSec_1Dth_antinu.h" -#include "MINERvA_CCCOHPI_XSec_1Dth_nu.h" - -#include "MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.h" -#include "MINERvA_CC0pi_XSec_1DQ2_Tgt_nu.h" -#include "MINERvA_CC1pip_XSec_1D_2017Update.h" - -#endif - -#ifndef __NO_T2K__ -// T2K CC0pi -#include "T2K_CC0pi_XSec_2DPcos_nu.h" - -// T2K CC1pi+ on CH -#include "T2K_CC1pip_CH_XSec_1DQ2_nu.h" -#include "T2K_CC1pip_CH_XSec_1DWrec_nu.h" -#include "T2K_CC1pip_CH_XSec_1Dpmu_nu.h" -#include "T2K_CC1pip_CH_XSec_1Dppi_nu.h" -#include "T2K_CC1pip_CH_XSec_1Dq3_nu.h" -#include "T2K_CC1pip_CH_XSec_1Dthmupi_nu.h" -#include "T2K_CC1pip_CH_XSec_1Dthpi_nu.h" -#include "T2K_CC1pip_CH_XSec_1Dthq3pi_nu.h" - -// T2K CC1pi+ on H2O -#include "T2K_CC1pip_H2O_XSec_1DEnuDelta_nu.h" -#include "T2K_CC1pip_H2O_XSec_1DEnuMB_nu.h" -#include "T2K_CC1pip_H2O_XSec_1Dcosmu_nu.h" -#include "T2K_CC1pip_H2O_XSec_1Dcosmupi_nu.h" -#include "T2K_CC1pip_H2O_XSec_1Dcospi_nu.h" -#include "T2K_CC1pip_H2O_XSec_1Dpmu_nu.h" -#include "T2K_CC1pip_H2O_XSec_1Dppi_nu.h" - -// T2K STV CC0pi -#include "T2K_CC0pinp_STV_XSec_1Ddpt_nu.h" - -#include "T2K_CC0pi_XSec_2DPcos_nu_nonuniform.h" - -#endif - -#ifndef __NO_SciBooNE__ - -// SciBooNE COH studies -#include "SciBooNE_CCCOH_1TRK_1DQ2_nu.h" -#include "SciBooNE_CCCOH_MuPiNoVA_1DQ2_nu.h" -#include "SciBooNE_CCCOH_MuPiNoVA_1Dthetapi_nu.h" -#include "SciBooNE_CCCOH_MuPiNoVA_1Dthetapr_nu.h" -#include "SciBooNE_CCCOH_MuPiVA_1DQ2_nu.h" -#include "SciBooNE_CCCOH_MuPr_1DQ2_nu.h" -#include "SciBooNE_CCCOH_STOPFINAL_1DQ2_nu.h" -#include "SciBooNE_CCCOH_STOP_NTrks_nu.h" -#endif - -#ifndef __NO_K2K__ -// K2K NC1pi0 -#include "K2K_NC1pi0_Evt_1Dppi0_nu.h" -#endif - -// MC Studies -#include "ExpMultDist_CCQE_XSec_1DVar_FakeStudy.h" -#include "ExpMultDist_CCQE_XSec_2DVar_FakeStudy.h" -#include "MCStudy_CCQEHistograms.h" - -#include "GenericFlux_Tester.h" -#include "GenericFlux_Vectors.h" - -#include "ElectronFlux_FlatTree.h" -#include "ElectronScattering_DurhamData.h" -#include "MCStudy_KaonPreSelection.h" -#include "MCStudy_MuonValidation.h" - -#include "OfficialNIWGPlots.h" -#include "T2K2017_FakeData.h" - -#include "Simple_Osc.h" - -#include "FitWeight.h" - -#include "NuisConfig.h" -#include "NuisKey.h" +class FitWeight; +class nuiskey; +class MeasurementBase; //! Functions to make it easier for samples to be created and handled. namespace SampleUtils { //! Create a given sample given its name, file, type, fakdata(fkdt) file and the //! current rw engine and push it back into the list fChain. /*bool LoadSample(std::list* fChain, std::string name, std::string file, std::string type, std::string fkdt, FitWeight* rw);*/ MeasurementBase* CreateSample(std::string name, std::string file, std::string type, std::string fkdt, FitWeight* rw); MeasurementBase* CreateSample(nuiskey samplekey); } /*! @} */ #endif diff --git a/src/FitBase/ParamPull.cxx b/src/FitBase/ParamPull.cxx index 4a0c6bc..8f9cb6e 100644 --- a/src/FitBase/ParamPull.cxx +++ b/src/FitBase/ParamPull.cxx @@ -1,810 +1,844 @@ // 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 "ParamPull.h" //******************************************************************************* ParamPull::ParamPull(std::string name, std::string inputfile, std::string type, std::string dials) { //******************************************************************************* fMinHist = NULL; fMaxHist = NULL; fTypeHist = NULL; fDialSelection = dials; + fLimitHist = NULL; fName = name; fInput = inputfile; fType = type; // Set the pull type SetType(fType); std::cout << fType << std::endl; // Setup Histograms from input file SetupHistograms(fInput); }; //******************************************************************************* void ParamPull::SetType(std::string type) { //******************************************************************************* fType = type; // Assume Default if empty if (type.empty() || type == "DEFAULT") { ERR(WRN) << "No type specified for ParmPull class " << fName << std::endl; ERR(WRN) << "Assuming GAUSTHROW/GAUSPULL" << std::endl; type = "GAUSTHROW/GAUSPULL"; } // Set Dial options if (type.find("FRAC") != std::string::npos) { fDialOptions = "FRAC"; fPlotTitles = ";; Fractional RW Value"; } else if (type.find("ABS") != std::string::npos) { fDialOptions = "ABS"; fPlotTitles = ";; ABS RW Value"; } else { fDialOptions = ""; fPlotTitles = ";; RW Value"; } // Parse throw types if (type.find("GAUSPULL") != std::string::npos) fCalcType = kGausPull; else fCalcType = kNoPull; if (type.find("GAUSTHROW") != std::string::npos) fThrowType = kGausThrow; + else if (type.find("FLATTHROW") != std::string::npos) fThrowType = kFlatThrow; else fThrowType = kNoThrow; // Extra check to see if throws or pulls are turned off if (type.find("NOPULL") != std::string::npos) fCalcType = kNoPull; if (type.find("NOTHROW") != std::string::npos) fThrowType = kNoThrow; } //******************************************************************************* void ParamPull::SetupHistograms(std::string input) { //******************************************************************************* // Extract Types from Input fFileType = ""; const int nfiletypes = 4; const std::string filetypes[nfiletypes] = {"FIT", "ROOT", "TXT", "DIAL"}; for (int i = 0; i < nfiletypes; i++) { std::string tempTypes = filetypes[i] + ":"; if (input.find(tempTypes) != std::string::npos) { fFileType = filetypes[i]; input.replace(input.find(tempTypes), tempTypes.size(), ""); break; } } // Read Files if (!fFileType.compare("FIT")) ReadFitFile(input); else if (!fFileType.compare("ROOT")) ReadRootFile(input); else if (!fFileType.compare("VECT")) ReadVectFile(input); else if (!fFileType.compare("DIAL")) ReadDialInput(input); else { ERR(FTL) << "Unknown ParamPull Type: " << input << std::endl; throw; } // Check Dials are all good CheckDialsValid(); // Setup MC Histogram fMCHist = (TH1D*) fDataHist->Clone(); fMCHist->Reset(); fMCHist->SetNameTitle( (fName + "_MC").c_str(), (fName + " MC" + fPlotTitles).c_str() ); // If no Covar input make an uncorrelated one if (!fCovar) { fCovar = StatUtils::MakeDiagonalCovarMatrix(fDataHist, 1.0); } // If no types or limits are provided give them a default option if (!fMinHist) { fMinHist = (TH1D*) fDataHist->Clone(); fMinHist->SetNameTitle( (fName + "_min").c_str(), (fName + " min" + fPlotTitles).c_str() ); for (int i = 0; i < fMinHist->GetNbinsX(); i++) { // TODO (P.Stowell) Change this to a NULL system where limits are actually free! fMinHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) - 1E6); } } if (!fMaxHist) { fMaxHist = (TH1D*) fDataHist->Clone(); fMaxHist->SetNameTitle( (fName + "_min").c_str(), (fName + " min" + fPlotTitles).c_str() ); for (int i = 0; i < fMaxHist->GetNbinsX(); i++) { fMaxHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) - 1E6); } } // Set types from state, or to unknown if (!fTypeHist) { int deftype = -1; if (fType.find("T2K") != std::string::npos) { deftype = kT2K; } else if (fType.find("NEUT") != std::string::npos) { deftype = kNEUT; } else if (fType.find("NIWG") != std::string::npos) { deftype = kNIWG; } else if (fType.find("GENIE") != std::string::npos) { deftype = kGENIE; } else if (fType.find("NORM") != std::string::npos) { deftype = kNORM; } else if (fType.find("NUWRO") != std::string::npos) { deftype = kNUWRO; } fTypeHist = new TH1I( (fName + "_type").c_str(), (fName + " type" + fPlotTitles).c_str(), fDataHist->GetNbinsX(), 0, fDataHist->GetNbinsX() ); for (int i = 0; i < fTypeHist->GetNbinsX(); i++) { fTypeHist->SetBinContent(i + 1, deftype ); } } // Sort Covariances fInvCovar = StatUtils::GetInvert(fCovar); fDecomp = StatUtils::GetDecomp(fCovar); // Create DataTrue for Throws fDataTrue = (TH1D*) fDataHist->Clone(); fDataTrue->SetNameTitle( (fName + "_truedata").c_str(), (fName + " truedata" + fPlotTitles).c_str() ); fDataOrig = (TH1D*) fDataHist->Clone(); fDataOrig->SetNameTitle( (fName + "_origdata").c_str(), (fName + " origdata" + fPlotTitles).c_str() ); // Select only dials we want if (!fDialSelection.empty()) { (*fDataHist) = RemoveBinsNotInString(*fDataHist, fDialSelection); } } //******************************************************************************* TH1D ParamPull::RemoveBinsNotInString(TH1D hist, std::string mystr) { //******************************************************************************* // Make list of allowed bins std::vector allowedbins; for (int i = 0; i < hist.GetNbinsX(); i++) { std::string syst = std::string(hist.GetXaxis()->GetBinLabel(i + 1)); if (mystr.find(syst) != std::string::npos) { allowedbins.push_back(syst); } } // Make new histogram UInt_t nbins = allowedbins.size(); TH1D newhist = TH1D( hist.GetName(), hist.GetTitle(), (Int_t)nbins, 0.0, (Double_t)nbins); // Setup bins for (UInt_t i = 0; i < nbins; i++) { // Set Labels newhist.GetXaxis()->SetBinLabel(i + 1, allowedbins[i].c_str()); // Copy Values for (Int_t j = 0; j < hist.GetNbinsX(); j++) { if (!allowedbins[i].compare(hist.GetXaxis()->GetBinLabel(j + 1))) { newhist.SetBinContent(i + 1, hist.GetBinContent(j + 1) ); newhist.SetBinError(i + 1, hist.GetBinError(j + 1) ); } } } return newhist; } //******************************************************************************* TH1I ParamPull::RemoveBinsNotInString(TH1I hist, std::string mystr) { //******************************************************************************* // Make list of allowed bins std::vector allowedbins; for (int i = 0; i < hist.GetNbinsX(); i++) { std::string syst = std::string(hist.GetXaxis()->GetBinLabel(i + 1)); if (mystr.find(syst) != std::string::npos) { allowedbins.push_back(syst); } } // Make new histogram UInt_t nbins = allowedbins.size(); TH1I newhist = TH1I( hist.GetName(), hist.GetTitle(), (Int_t)nbins, 0.0, (Int_t)nbins); // Setup bins for (UInt_t i = 0; i < nbins; i++) { // Set Labels newhist.GetXaxis()->SetBinLabel(i + 1, allowedbins[i].c_str()); // Copy Values for (Int_t j = 0; j < hist.GetNbinsX(); j++) { if (!allowedbins[i].compare(hist.GetXaxis()->GetBinLabel(j + 1))) { newhist.SetBinContent(i + 1, hist.GetBinContent(j + 1) ); newhist.SetBinError(i + 1, hist.GetBinError(j + 1) ); } } } return newhist; } //******************************************************************************* void ParamPull::ReadFitFile(std::string input) { //******************************************************************************* TFile* tempfile = new TFile(input.c_str(), "READ"); // Read Data fDataHist = (TH1D*) tempfile->Get("fit_dials_free"); if (!fDataHist) { ERR(FTL) << "Can't find TH1D hist fit_dials in " << fName << std::endl; ERR(FTL) << "File Entries:" << std::endl; tempfile->ls(); throw; } fDataHist->SetDirectory(0); fDataHist->SetNameTitle( (fName + "_data").c_str(), (fName + " data" + fPlotTitles).c_str() ); // Read Covar TH2D* tempcov = (TH2D*) tempfile->Get("covariance_free"); if (!tempcov) { ERR(FTL) << "Can't find TH2D covariance_free in " << fName << std::endl; ERR(FTL) << "File Entries:" << std::endl; tempfile->ls(); throw; } // Setup Covar int nbins = fDataHist->GetNbinsX(); fCovar = new TMatrixDSym( nbins ); for (int i = 0; i < nbins; i++) { for (int j = 0; j < nbins; j++) { (*fCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1); } } return; } //******************************************************************************* void ParamPull::ReadRootFile(std::string input) { //******************************************************************************* std::vector inputlist = GeneralUtils::ParseToStr(input, ";"); // Check all given if (inputlist.size() < 2) { ERR(FTL) << "Covar supplied in 'ROOT' format should have 3 semi-colon seperated entries!" << std::endl << "ROOT:filename;histname[;covarname]" << std::endl; ERR(FTL) << "histname = TH1D, covarname = TH2D" << std::endl; throw; } // Get Entries std::string filename = inputlist[0]; LOG(DEB) << filename << std::endl; std::string histname = inputlist[1]; LOG(DEB) << histname << std::endl; LOG(DEB) << input << std::endl; // Read File TFile* tempfile = new TFile(filename.c_str(), "READ"); if (tempfile->IsZombie()) { ERR(FTL) << "Can't find file in " << fName << std::endl; ERR(FTL) << "location = " << filename << std::endl; throw; } // Read Hist fDataHist = (TH1D*) tempfile->Get(histname.c_str()); if (!fDataHist) { ERR(FTL) << "Can't find TH1D hist " << histname << " in " << fName << std::endl; ERR(FTL) << "File Entries:" << std::endl; tempfile->ls(); throw; } fDataHist->SetDirectory(0); fDataHist->SetNameTitle( (fName + "_data").c_str(), (fName + " data" + fPlotTitles).c_str() ); LOG(DEB) << "READING COVAR" << std::endl; // Read Covar if (inputlist.size() > 2) { std::string covarname = inputlist[2]; LOG(DEB) << "COVARNAME = " << covarname << std::endl; TH2D* tempcov = (TH2D*) tempfile->Get(covarname.c_str()); if (!tempcov) { ERR(FTL) << "Can't find TH2D covar " << covarname << " in " << fName << std::endl; ERR(FTL) << "File Entries:" << std::endl; tempfile->ls(); throw; } // Setup Covar int nbins = fDataHist->GetNbinsX(); fCovar = new TMatrixDSym( nbins ); for (int i = 0; i < nbins; i++) { for (int j = 0; j < nbins; j++) { (*fCovar)(i, j) = tempcov->GetBinContent(i + 1, j + 1); } } // Uncorrelated } else { LOG(SAM) << "No Covar provided so using diagonal errors for " << fName << std::endl; fCovar = NULL; } } //******************************************************************************* void ParamPull::ReadVectFile(std::string input) { //******************************************************************************* std::vector inputlist = GeneralUtils::ParseToStr(input, ";"); if (inputlist.size() < 4) { ERR(FTL) << "Need 3 inputs for vector input in " << fName << std::endl; ERR(FTL) << "Inputs: " << input << std::endl; throw; } // Open File std::string rootname = inputlist[0]; TFile* tempfile = new TFile(rootname.c_str(), "READ"); if (tempfile->IsZombie()) { ERR(FTL) << "Can't find file in " << fName << std::endl; ERR(FTL) << "location = " << rootname << std::endl; throw; } // Get Name std::string tagname = inputlist[1]; // TVector dialtags = tempfile->Get(tagname.c_str()); // if (!dialtags){ // ERR(FTL) << "Can't find list of dial names!" << std::endl; // } // Get Values std::string valuename = inputlist[2]; TVectorD* dialvals = (TVectorD*)tempfile->Get(valuename.c_str()); if (!dialvals) { ERR(FTL) << "Can't find dial values" << std::endl; } // Get Matrix std::string matrixname = inputlist[3]; TMatrixD* matrixvals = (TMatrixD*)tempfile->Get(matrixname.c_str()); if (!matrixvals) { ERR(FTL) << "Can't find matirx values" << std::endl; } // Get Types if (inputlist.size() > 4) { std::string typesname = inputlist[3]; } // Get Minimum if (inputlist.size() > 5) { std::string minname = inputlist[4]; } // Get Maximum if (inputlist.size() > 6) { std::string maxname = inputlist[5]; } } //******************************************************************************* void ParamPull::ReadDialInput(std::string input) { //******************************************************************************* std::vector inputlist = GeneralUtils::ParseToStr(input, ";"); if (inputlist.size() < 3) { ERR(FTL) << "Need 3 inputs for dial input in " << fName << std::endl; ERR(FTL) << "Inputs: " << input << std::endl; throw; } std::vector inputvals = GeneralUtils::ParseToDbl(input, ";"); std::string dialname = inputlist[0]; double val = inputvals[1]; double err = inputvals[2]; fDataHist = new TH1D( (fName + "_data").c_str(), (fName + "_data" + fPlotTitles).c_str(), 1, 0, 1); fDataHist->SetBinContent(1, val); fDataHist->SetBinError(1, err); fDataHist->GetXaxis()->SetBinLabel(1, dialname.c_str()); + + fLimitHist = new TH1D( (fName + "_limits").c_str(), + (fName + "_limits" + fPlotTitles).c_str(), 1, 0, 1); + fLimitHist->Reset(); + if (inputvals.size() > 4) { + fLimitHist->SetBinContent(1, (inputvals[3] + inputvals[4]) / 2.0); + fLimitHist->SetBinError(1, (inputvals[4] - inputvals[3]) / 2.0); + } + fCovar = NULL; } //******************************************************************************* std::map ParamPull::GetAllDials() { //******************************************************************************* std::map dialtypemap; for (int i = 0; i < fDataHist->GetNbinsX(); i++) { std::string name = fDataHist->GetXaxis()->GetBinLabel(i + 1); int type = fTypeHist->GetBinContent(i + 1); dialtypemap[name] = type; } return dialtypemap; } //******************************************************************************* bool ParamPull::CheckDialsValid() { //******************************************************************************* return true; std::string helpstring = ""; for (int i = 0; i < fDataHist->GetNbinsX(); i++) { std::string name = std::string(fDataHist->GetXaxis()->GetBinLabel(i + 1)); // If dial exists its all good if (FitBase::GetRW()->DialIncluded(name)) continue; // If it doesn't but its a sample norm also continue if (name.find("_norm") != std::string::npos) { ERR(WRN) << "Norm dial included in covar but not set in FitWeight." << std::endl; ERR(WRN) << "Assuming its a sample norm and skipping..." << std::endl; } // Dial unknown so print a help statement std::ostringstream tempstr; tempstr << "unknown_parameter " << name << " " << fDataHist->GetBinContent(i + 1) << " " << fDataHist->GetBinContent(i + 1) - fDataHist->GetBinError(i + 1) << " " << fDataHist->GetBinContent(i + 1) + fDataHist->GetBinError(i + 1) << " " << fDataHist->GetBinError(i + 1) << " "; if (!fType.empty()) tempstr << fType << std::endl; else tempstr << "FREE" << std::endl; helpstring += tempstr.str(); } // Show statement before failing if (!helpstring.empty()) { ERR(WRN) << "Dial(s) included in covar but not set in FitWeight." << std::endl << "ParamPulls needs to know how you want it to be treated." << std::endl << "Include the following lines into your card:" << std::endl; ERR(WRN) << helpstring << std::endl; throw; return false; } else { return true; } } //******************************************************************************* void ParamPull::Reconfigure() { //******************************************************************************* FitWeight* rw = FitBase::GetRW(); // Get Dial Names that are valid std::vector namevec = rw->GetDialNames(); std::vector valuevec = rw->GetDialValues(); // Set Bin Values from RW for (UInt_t i = 0; i < namevec.size(); i++) { // Loop over bins and check name matches std::string syst = namevec.at(i); double systval = valuevec.at(i); std::vector allsyst = GeneralUtils::ParseToStr(syst, ","); // Proper Reconf using RW for (int j = 0; j < fMCHist->GetNbinsX(); j++) { // Search for the name of this bin in the corrent dial std::string binname = std::string(fMCHist->GetXaxis()->GetBinLabel(j + 1) ); // Check Full Name if (!syst.compare(binname.c_str())) { fMCHist->SetBinContent(j + 1, systval); break; } std::vector splitbinname = GeneralUtils::ParseToStr(binname, ","); for (size_t l = 0; l < splitbinname.size(); l++) { std::string singlebinname = splitbinname[l]; for (size_t k = 0; k < allsyst.size(); k++) { if (!allsyst[k].compare(singlebinname.c_str())) { fMCHist->SetBinContent(j + 1, systval); } } } } } return; }; //******************************************************************************* void ParamPull::ResetToy(void) { //******************************************************************************* if (fDataHist) delete fDataHist; LOG(DEB) << "Resetting toy" << std::endl; LOG(DEB) << fDataTrue << std::endl; fDataHist = (TH1D*)fDataTrue->Clone(); LOG(DEB) << "Setting name" << std::endl; fDataHist->SetNameTitle( (fName + "_data").c_str(), (fName + " data" + fPlotTitles).c_str() ); } //******************************************************************************* void ParamPull::SetFakeData(std::string fakeinput) { //******************************************************************************* // Set from MC Setting if (!fakeinput.compare("MC")) { // Copy MC into data if (fDataHist) delete fDataHist; fDataHist = (TH1D*)fMCHist->Clone(); fDataHist->SetNameTitle( (fName + "_data").c_str(), (fName + " fakedata" + fPlotTitles).c_str() ); // Copy original data errors for (int i = 0; i < fDataOrig->GetNbinsX(); i++) { fDataHist->SetBinError(i + 1, fDataOrig->GetBinError(i + 1) ); } // Make True Toy Central Value Hist fDataTrue = (TH1D*) fDataHist->Clone(); fDataTrue->SetNameTitle( (fName + "_truedata").c_str(), (fName + " truedata" + fPlotTitles).c_str() ); } else { ERR(FTL) << "Trying to set fake data for ParamPulls not from MC!" << std::endl; ERR(FTL) << "Not currently implemented.." << std::endl; throw; } } //******************************************************************************* void ParamPull::RemoveFakeData() { //******************************************************************************* delete fDataHist; fDataHist = (TH1D*)fDataOrig->Clone(); fDataHist->SetNameTitle( (fName + "_data").c_str(), (fName + " data" + fPlotTitles).c_str() ); fDataTrue = (TH1D*) fDataHist->Clone(); fDataTrue->SetNameTitle( (fName + "_truedata").c_str(), (fName + " truedata" + fPlotTitles).c_str() ); } //******************************************************************************* double ParamPull::GetLikelihood() { //******************************************************************************* double like = 0.0; switch (fCalcType) { // Gaussian Calculation with correlations case kGausPull: like = StatUtils::GetChi2FromCov(fDataHist, fMCHist, fInvCovar, NULL); like *= 1E-76; break; // Default says this has no pull case kNoThrow: default: like = 0.0; break; } LOG(DEB) << "Likelihood = " << like << " " << fCalcType << std::endl; return like; }; //******************************************************************************* int ParamPull::GetNDOF() { //******************************************************************************* int ndof = 0; if (fCalcType != kNoThrow) { - ndof += fDataHist->GetNbinsX(); + ndof = fDataHist->GetNbinsX(); } return ndof; }; //******************************************************************************* void ParamPull::ThrowCovariance() { //******************************************************************************* // Reset toy for throw ResetToy(); - LOG(DEB) << "Toy Reset " << std::endl; + LOG(FIT) << "Creating new toy dataset" << std::endl; // Generate random Gaussian throws std::vector randthrows; for (int i = 0; i < fDataHist->GetNbinsX(); i++) { double randtemp = 0.0; switch (fThrowType) { // Gaussian Throws case kGausThrow: randtemp = gRandom->Gaus(0.0, 1.0); break; + // Uniform Throws + case kFlatThrow: + randtemp = gRandom->Uniform(0.0, 1.0); + if (fLimitHist) { + randtemp = fLimitHist->GetBinContent(i + 1) + fLimitHist->GetBinError(i + 1) * ( randtemp * 2 - 1 ); + } + break; + // No Throws (DEFAULT) default: break; } randthrows.push_back(randtemp); } // Create Bin Modifications double totalres = 0.0; for (int i = 0; i < fDataHist->GetNbinsX(); i++) { // Calc Bin Mod double binmod = 0.0; - for (int j = 0; j < fDataHist->GetNbinsX(); j++) { - // std::cout << "DECOMP " << j << " " << i << " " << randthrows.at(j) << std::endl; - binmod += (*fDecomp)(j, i) * randthrows.at(j); + + if (fThrowType == kGausThrow) { + for (int j = 0; j < fDataHist->GetNbinsX(); j++) { + binmod += (*fDecomp)(j, i) * randthrows.at(j); + } + } else if (fThrowType == kFlatThrow) { + binmod = randthrows.at(i) - fDataHist->GetBinContent(i + 1); } + // Add up fraction dif totalres += binmod; // Add to current data fDataHist->SetBinContent(i + 1, fDataHist->GetBinContent(i + 1) + binmod); } // Rename fDataHist->SetNameTitle( (fName + "_data").c_str(), (fName + " toydata" + fPlotTitles).c_str() ); - // Print Status - LOG(REC) << "Created new toy histogram. Total Dif = " - << totalres << std::endl; + // Check Limits + if (fLimitHist) { + for (int i = 0; i < fLimitHist->GetNbinsX(); i++) { + if (fLimitHist->GetBinError(i + 1) == 0.0) continue; + if (fDataHist->GetBinContent(i + 1) > fLimitHist->GetBinContent(i + 1) + fLimitHist->GetBinError(i + 1) || + fDataHist->GetBinContent(i + 1) < fLimitHist->GetBinContent(i + 1) - fLimitHist->GetBinError(i + 1)) { + this->ThrowCovariance(); + } + } + } + return; }; //******************************************************************************* TH2D ParamPull::GetCovar() { //******************************************************************************* TH2D tempCov = TH2D(*fInvCovar); for (int i = 0; i < tempCov.GetNbinsX(); i++) { tempCov.GetXaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1)); tempCov.GetYaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1)); } tempCov.SetNameTitle((fName + "_INVCOV").c_str(), (fName + " InvertedCovariance;Dials;Dials").c_str()); return tempCov; } //******************************************************************************* TH2D ParamPull::GetFullCovar() { //******************************************************************************* TH2D tempCov = TH2D(*fCovar); for (int i = 0; i < tempCov.GetNbinsX(); i++) { tempCov.GetXaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1)); tempCov.GetYaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1)); } tempCov.SetNameTitle((fName + "_COV").c_str(), (fName + " Covariance;Dials;Dials").c_str()); return tempCov; } //******************************************************************************* TH2D ParamPull::GetDecompCovar() { //******************************************************************************* TH2D tempCov = TH2D(*fCovar); for (int i = 0; i < tempCov.GetNbinsX(); i++) { tempCov.GetXaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1)); tempCov.GetYaxis()->SetBinLabel(i + 1, fDataHist->GetXaxis()->GetBinLabel(i + 1)); } tempCov.SetNameTitle((fName + "_DEC").c_str(), (fName + " Decomposition;Dials;Dials").c_str()); return tempCov; } //******************************************************************************* void ParamPull::Write(std::string writeoptt) { //******************************************************************************* fDataHist->Write(); fMCHist->Write(); - + if (fLimitHist) { + fLimitHist->Write(); + } GetCovar().Write(); GetFullCovar().Write(); GetDecompCovar().Write(); return; }; diff --git a/src/FitBase/ParamPull.h b/src/FitBase/ParamPull.h index 22e707d..7da9eac 100644 --- a/src/FitBase/ParamPull.h +++ b/src/FitBase/ParamPull.h @@ -1,188 +1,193 @@ // Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #ifndef PARAM_PULL_H_SEEN #define PARAM_PULL_H_SEEN /*! * \addtogroup FitBase * @{ */ #include #include #include #include #include #include #include #include // ROOT includes #include #include #include #include #include #include #include // Fit Includes #include "PlotUtils.h" #include "StatUtils.h" #include "FitWeight.h" #include "FitLogger.h" #include "EventManager.h" #include "TVector.h" using namespace std; //! Enums to allow Non Guassian Pulls in Future enum FitPullTypes { kUnknownPull = -1, kNoPull = 0, kGausPull = 1 }; enum FitThrowTypes { kUnknownThrow = -1, kNoThrow = 0, - kGausThrow = 1 + kGausThrow = 1, + kFlatThrow = 2 }; //! Used to produce gaussian penalty terms in the fit. class ParamPull { public: //! Default Constructor ParamPull(std::string name, std::string inputfile, std::string type, std::string dials=""); //! Default destructor virtual ~ParamPull(void) {}; // Set dial types (DEFAULT,ABS,FRAC) void SetType(std::string type); // Setup Histogram inputs (from previous fit file, or ROOT file) void SetupHistograms(std::string input); // Read a previous NUISANCE file void ReadFitFile(std::string input); // Read a ROOT file with any histograms in void ReadRootFile(std::string input); // Read Text file void ReadVectFile(std::string input); // Read a single dial central value and error void ReadDialInput(std::string input); //! Reconfigure function reads in current RW engine dials and sets their value to MC void Reconfigure(void); //! Get likelihood given the current values double GetLikelihood(void); //! Get NDOF if used in likelihoods int GetNDOF(void); // Get Covariance Matrices as plots TH2D GetCovar(void); TH2D GetFullCovar(void); TH2D GetDecompCovar(void); // Get Covariance Matrices inline TMatrixDSym GetCovarMatrix (void) const { return *fInvCovar; }; inline TMatrixDSym GetFullCovarMatrix (void) const { return *fCovar; }; inline TMatrixDSym GetDecompCovarMatrix (void) const { return *fDecomp; }; //! Save the histograms void Write(std::string writeopt=""); //! Throw the dial values using the current covariance. Useful for parameter throws. void ThrowCovariance(void); //! Compare dials to RW bool CheckDialsValid(void); //! Reset toy data back to original data void ResetToy(void); //! Read fake data from MC void SetFakeData(std::string fakeinput); //! Reset fake data back to original data (before fake data or throws) void RemoveFakeData(); // Get Functions inline std::string GetName (void) const { return fName; }; inline std::string GetInput (void) const { return fInput; }; inline std::string GetType (void) const { return fType; }; inline std::string GetFileType (void) const { return fFileType; }; inline std::string GetDialOptions (void) const { return fDialOptions; }; std::map GetAllDials(); inline TH1D GetDataHist (void) const { return *fDataHist; }; inline TH1D GetDataTrue (void) const { return *fDataTrue; }; inline TH1D GetDataOrig (void) const { return *fDataOrig; }; inline TH1D GetMCHist (void) const { return *fMCHist; }; inline TH1D GetMaxHist (void) const { return *fMaxHist; }; inline TH1D GetMinHist (void) const { return *fMinHist; }; inline TH1I GetDialTypes (void) const { return *fTypeHist; }; + inline TH1D GetLimitHist (void) const { return *fLimitHist; }; + private: TH1D RemoveBinsNotInString(TH1D hist, std::string mystr); TH1I RemoveBinsNotInString(TH1I hist, std::string mystr); std::string fName; //!< Pull Name std::string fInput; //!< Pull input string std::string fType; //!< Pull options type std::string fFileType; //!< Pull input file types std::string fPlotTitles; //! Axis format std::string fDialOptions; //!< Dial handling options std::string fDialSelection; //!< Dial Selection TH1D* fMCHist; //!< Current MC Histogram TH1D* fDataHist; //!< Current data Histogram TH1D* fDataTrue; //!< True Data (before histogram throws) TH1D* fDataOrig; //!< Orig Data (without toys or fake data) TH1D* fMaxHist; //!< Maximum limit on MC/Data TH1D* fMinHist; //!< Maximum limit on MC/Data TH1I* fTypeHist; //!< Dial Types int fCalcType; //!< Method to calculate likelihood int fThrowType; //!< Method to calculate throws TMatrixDSym* fCovar; //!< Covariance TMatrixDSym* fInvCovar; //!< Inverted Covariance TMatrixDSym* fDecomp; //!< Decomposition + + TH1D* fLimitHist; }; // Class TypeDefs typedef std::list::const_iterator PullListConstIter; typedef std::list::iterator PullListIter; typedef std::vector::const_iterator PullVectConstIter; typedef std::vector::iterator PullVectIter; /*! @} */ #endif diff --git a/src/FitBase/SampleSettings.cxx b/src/FitBase/SampleSettings.cxx index 5388c0d..f574209 100644 --- a/src/FitBase/SampleSettings.cxx +++ b/src/FitBase/SampleSettings.cxx @@ -1,226 +1,231 @@ #include "SampleSettings.h" SampleSettings::SampleSettings() { }; SampleSettings::SampleSettings(nuiskey key) { fKeyValues = key; if (!key.Has("type")) key.AddS("type", "DEFAULT"); } std::string SampleSettings::GetName() { return GetS("name"); } void SampleSettings::SetS(std::string name, std::string val) { fKeyValues.SetS(name, val); }; bool SampleSettings::Found(std::string name, std::string substr) { std::string compstring = fKeyValues.GetS(name); return compstring.find(substr) != std::string::npos; } void SampleSettings::SetXTitle(std::string name) { SetDefault("xtitle", name); }; void SampleSettings::SetYTitle(std::string name) { SetDefault("ytitle", name); }; void SampleSettings::SetZTitle(std::string name) { SetDefault("ztitle", name); }; void SampleSettings::SetNormError(double norm) { SetDefault("norm_error", GeneralUtils::DblToStr(norm)); }; double SampleSettings::GetNormError() { return GetD("norm_error"); }; std::string SampleSettings::GetCovarInput() { return GetS("covar"); } void SampleSettings::SetAllowedTypes(std::string allowed, std::string defaulttype) { SetDefault("default_types", allowed); SetDefault("allowed_types", defaulttype); }; void SampleSettings::SetEnuRangeFromFlux(TH1D* fluxhist) { double enu_min = fluxhist->GetXaxis()->GetXmin(); double enu_max = fluxhist->GetXaxis()->GetXmax(); SetEnuRange(enu_min, enu_max); }; void SampleSettings::SetEnuRange(double min, double max) { SetDefault("enu_min", min); SetDefault("enu_max", max); }; void SampleSettings::Set(std::string name, double d) { SetDefault(name, d); } void SampleSettings::Set(std::string name, int i) { SetDefault(name, i); } void SampleSettings::Set(std::string name, std::string s) { SetDefault(name, s); } void SampleSettings::DefineAllowedTargets(std::string targ) { // fAllowedTargets = TargetUtils::ParseTargetsToIntVect(targ); }; void SampleSettings::FoundFill(std::string name, std::string substr, bool& cont, bool def) { std::string compstring = fKeyValues.GetS(name); if (compstring.find(substr) != std::string::npos) { cont = def; } else { cont = !def; } }; void SampleSettings::SetTitle(std::string val) { SetDefault("title", val); }; void SampleSettings::SetDataInput(std::string val) { SetDefault("data", val); }; std::string SampleSettings::GetMapInput() { return GetS("map"); } void SampleSettings::SetMapInput(std::string val) { SetDefault("map", val); } void SampleSettings::SetErrorInput(std::string val) { SetDefault("error", val); }; void SampleSettings::SetCovarInput(std::string val) { SetDefault("covar", val); }; +void SampleSettings::SetShapeCovarInput(std::string val) { + SetDefault("shapecovar", val); +}; + + void SampleSettings::SetDefault(std::string name, std::string val) { if (!fKeyValues.Has(name)) fKeyValues.AddS(name, val); }; void SampleSettings::SetDefault(std::string name, double val ) { if (!fKeyValues.Has(name)) fKeyValues.AddD(name, val); }; void SampleSettings::SetHasExtraHistograms(bool opt) { fHasExtraHistograms = opt; }; void SampleSettings::DefineAllowedSpecies(std::string species) { fAllowedTargets = BeamUtils::ParseSpeciesToIntVect(species); }; std::string SampleSettings::Title() { return GetS("title"); } std::string SampleSettings::GetDataInput() { return GetS("data"); }; std::string SampleSettings::GetErrorInput() { return GetS("error"); }; std::string SampleSettings::PlotTitles() { return ";" + GetS("xtitle") + ";" + GetS("ytitle"); }; std::string SampleSettings::GetFullTitles() { return Title() + PlotTitles(); } int SampleSettings::GetI(std::string name) { return fKeyValues.GetI(name); } bool SampleSettings::GetB(std::string name){ return fKeyValues.GetB(name); } double SampleSettings::GetD(std::string name) { return fKeyValues.GetD(name); } std::string SampleSettings::GetS(std::string name) { return fKeyValues.GetS(name); } void SampleSettings::SetSuggestedFlux(std::string str) { SetS("suggested_flux", str); } void SampleSettings::SetDescription(std::string str) { SetS("description", str); } void SampleSettings::Write(std::string name) { if (name.empty()) name = this->GetS("name") + "_settings"; // Make a new canvas TCanvas* c1 = new TCanvas( name.c_str(), name.c_str(), 800, 600 ); c1->cd(); // Make new TPaveText TPaveText pttitle = TPaveText(0.05, 0.85, 0.95, 0.95); pttitle.AddText( name.c_str() ); pttitle.SetTextAlign(11); pttitle.SetTextSize(15); pttitle.Draw(); TPaveText pt = TPaveText(0.05, 0.05, 0.95, 0.80); std::vector keys = fKeyValues.GetAllKeys(); for (int i = 0; i < keys.size(); i++) { std::string keyval = fKeyValues.GetS(keys[i]); std::vector lines = GeneralUtils::ParseToStr(keyval, "\n"); if (lines.size() == 1) { pt.AddText( ("#bullet #bf{" + keys[i] + "} : " + fKeyValues.GetS(keys[i])).c_str() ); } else { pt.AddText( ("#bullet #bf{" + keys[i] + "} : ").c_str() ); for (int j = 0; j < lines.size(); j++) { pt.AddText((" |--> " + lines[j]).c_str() ); } } } pt.SetTextAlign(11); pt.SetTextSize(14); pt.Draw("SAME"); c1->Write(); delete c1; } void SampleSettings::SetOnlyMC(bool state) { SetDefault("onlymc", state); } diff --git a/src/FitBase/SampleSettings.h b/src/FitBase/SampleSettings.h index e689eeb..358a584 100644 --- a/src/FitBase/SampleSettings.h +++ b/src/FitBase/SampleSettings.h @@ -1,74 +1,77 @@ #ifndef SAMPLESETTINGS_H #define SAMPLESETTINGS_H #include "NuisConfig.h" #include "NuisKey.h" #include "TH1D.h" #include "BeamUtils.h" #include "TPaveText.h" #include "TCanvas.h" class SampleSettings { public: SampleSettings(); SampleSettings(nuiskey key); inline std::string Name(){return GetName();}; std::string GetName(); void SetS(std::string name, std::string val); void SetXTitle(std::string name); void SetYTitle(std::string name); void SetZTitle(std::string name); void SetNormError(double norm); double GetNormError(); void SetAllowedTypes(std::string allowed, std::string defaulttype="FIX"); void SetEnuRangeFromFlux(TH1D* fluxhist); void SetEnuRange(double min, double max); std::string Title(); void DefineAllowedTargets(std::string targ); void FoundFill(std::string name, std::string substr, bool& cont, bool def); // void FoundFill(std::string name, std::string substr, double& val, double ) bool Found(std::string name, std::string substr); void SetTitle(std::string val); void SetDataInput(std::string val); void SetErrorInput(std::string val); std::string GetErrorInput(); std::string GetMapInput(); void SetMapInput(std::string val); void SetCovarInput(std::string val); + void SetShapeCovarInput(std::string val); void SetDefault(std::string name, std::string val); void SetDefault(std::string name, double val); void SetHasExtraHistograms(bool opt = true); void DefineAllowedSpecies(std::string species); void SetSuggestedFlux(std::string str); void SetDescription(std::string str); void Write(std::string name=""); std::string GetFullTitles(); bool Has(std::string name){return fKeyValues.Has(name);}; std::string GetDataInput(); std::string PlotTitles(); std::string GetS(std::string name); int GetI(std::string name); double GetD(std::string name); std::string GetCovarInput(); void SetOnlyMC(bool state=true); bool GetB(std::string name); void Set(std::string name, int i); void Set(std::string name, std::string s); void Set(std::string name, double d); + + std::vector fAllowedTargets; std::vector fAllowedSpecies; nuiskey fKeyValues; bool fHasExtraHistograms; }; #endif diff --git a/src/InputHandler/GENIEInputHandler.cxx b/src/InputHandler/GENIEInputHandler.cxx index b0efcf7..b03b8f7 100644 --- a/src/InputHandler/GENIEInputHandler.cxx +++ b/src/InputHandler/GENIEInputHandler.cxx @@ -1,390 +1,390 @@ #include "GENIEInputHandler.h" #ifdef __GENIE_ENABLED__ GENIEGeneratorInfo::~GENIEGeneratorInfo() { DeallocateParticleStack(); } void GENIEGeneratorInfo::AddBranchesToTree(TTree* tn) { tn->Branch("GenieParticlePDGs", &fGenieParticlePDGs, "GenieParticlePDGs/I"); } void GENIEGeneratorInfo::SetBranchesFromTree(TTree* tn) { tn->SetBranchAddress("GenieParticlePDGs", &fGenieParticlePDGs); } void GENIEGeneratorInfo::AllocateParticleStack(int stacksize) { fGenieParticlePDGs = new int[stacksize]; } void GENIEGeneratorInfo::DeallocateParticleStack() { delete fGenieParticlePDGs; } void GENIEGeneratorInfo::FillGeneratorInfo(NtpMCEventRecord* ntpl) { Reset(); // Check for GENIE Event if (!ntpl) return; if (!ntpl->event) return; // Cast Event Record GHepRecord* ghep = static_cast(ntpl->event); if (!ghep) return; // Fill Particle Stack GHepParticle* p = 0; TObjArrayIter iter(ghep); // Loop over all particles int i = 0; while ((p = (dynamic_cast((iter).Next())))) { if (!p) continue; // Get PDG fGenieParticlePDGs[i] = p->Pdg(); i++; } } void GENIEGeneratorInfo::Reset() { for (int i = 0; i < kMaxParticles; i++) { fGenieParticlePDGs[i] = 0; } } GENIEInputHandler::GENIEInputHandler(std::string const& handle, std::string const& rawinputs) { LOG(SAM) << "Creating GENIEInputHandler : " << handle << std::endl; // Run a joint input handling fName = handle; // Setup the TChain fGENIETree = new TChain("gtree"); fSaveExtra = FitPar::Config().GetParB("SaveExtraGenie"); fCacheSize = FitPar::Config().GetParI("CacheSize"); fMaxEvents = FitPar::Config().GetParI("MAXEVENTS"); // Loop over all inputs and grab flux, eventhist, and nevents std::vector inputs = InputUtils::ParseInputFileList(rawinputs); for (size_t inp_it = 0; inp_it < inputs.size(); ++inp_it) { // Open File for histogram access TFile* inp_file = new TFile(InputUtils::ExpandInputDirectories(inputs[inp_it]).c_str(), "READ"); if (!inp_file or inp_file->IsZombie()) { THROW( "GENIE File IsZombie() at : '" << inputs[inp_it] << "'" << std::endl << "Check that your file paths are correct and the file exists!" << std::endl << "$ ls -lh " << inputs[inp_it] ); } // Get Flux/Event hist TH1D* fluxhist = (TH1D*)inp_file->Get("nuisance_flux"); TH1D* eventhist = (TH1D*)inp_file->Get("nuisance_events"); if (!fluxhist or !eventhist) { ERROR(FTL, "Input File Contents: " << inputs[inp_it] ); inp_file->ls(); THROW( "GENIE FILE doesn't contain flux/xsec info." << std::endl << "Try running the app PrepareGENIE first on :" << inputs[inp_it] << std::endl << "$ PrepareGENIE -h" ); } // Get N Events TTree* genietree = (TTree*)inp_file->Get("gtree"); if (!genietree) { ERROR(FTL, "gtree not located in GENIE file: " << inputs[inp_it]); THROW("Check your inputs, they may need to be completely regenerated!"); throw; } int nevents = genietree->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 fGENIETree->AddFile( inputs[inp_it].c_str() ); } // Registor all our file inputs SetupJointInputs(); // Assign to tree fEventType = kGENIE; fGenieNtpl = NULL; fGENIETree->SetBranchAddress("gmcrec", &fGenieNtpl); fGENIETree->GetEntry(0); // Create Fit Event fNUISANCEEvent = new FitEvent(); fNUISANCEEvent->SetGenieEvent(fGenieNtpl); if (fSaveExtra) { fGenieInfo = new GENIEGeneratorInfo(); fNUISANCEEvent->AddGeneratorInfo(fGenieInfo); } fNUISANCEEvent->HardReset(); }; GENIEInputHandler::~GENIEInputHandler() { //if (fGenieGHep) delete fGenieGHep; - if (fGenieNtpl) delete fGenieNtpl; - if (fGENIETree) delete fGENIETree; - if (fGenieInfo) delete fGenieInfo; + //if (fGenieNtpl) delete fGenieNtpl; + //if (fGENIETree) delete fGENIETree; + //if (fGenieInfo) delete fGenieInfo; } void GENIEInputHandler::CreateCache() { if (fCacheSize > 0) { // fGENIETree->SetCacheEntryRange(0, fNEvents); fGENIETree->AddBranchToCache("*", 1); fGENIETree->SetCacheSize(fCacheSize); } } void GENIEInputHandler::RemoveCache() { // fGENIETree->SetCacheEntryRange(0, fNEvents); fGENIETree->AddBranchToCache("*", 0); fGENIETree->SetCacheSize(0); } FitEvent* GENIEInputHandler::GetNuisanceEvent(const UInt_t entry, const bool lightweight) { if (entry >= (UInt_t)fNEvents) return NULL; // Read Entry from TTree to fill NEUT Vect in BaseFitEvt; fGENIETree->GetEntry(entry); // Run NUISANCE Vector Filler if (!lightweight) { CalcNUISANCEKinematics(); } // Setup Input scaling for joint inputs fNUISANCEEvent->InputWeight = GetInputWeight(entry); return fNUISANCEEvent; } int GENIEInputHandler::GetGENIEParticleStatus(genie::GHepParticle* p, int mode) { /* kIStUndefined = -1, kIStInitialState = 0, / generator-level initial state / kIStStableFinalState = 1, / generator-level final state: particles to be tracked by detector-level MC / kIStIntermediateState = 2, kIStDecayedState = 3, kIStCorrelatedNucleon = 10, kIStNucleonTarget = 11, kIStDISPreFragmHadronicState = 12, kIStPreDecayResonantState = 13, kIStHadronInTheNucleus = 14, / hadrons inside the nucleus: marked for hadron transport modules to act on / kIStFinalStateNuclearRemnant = 15, / low energy nuclear fragments entering the record collectively as a 'hadronic blob' pseudo-particle / kIStNucleonClusterTarget = 16, // for composite nucleons before phase space decay */ int state = kUndefinedState; switch (p->Status()) { case genie::kIStNucleonTarget: case genie::kIStInitialState: case genie::kIStCorrelatedNucleon: case genie::kIStNucleonClusterTarget: state = kInitialState; break; case genie::kIStStableFinalState: state = kFinalState; break; case genie::kIStHadronInTheNucleus: if (abs(mode) == 2) state = kInitialState; else state = kFSIState; break; case genie::kIStPreDecayResonantState: case genie::kIStDISPreFragmHadronicState: case genie::kIStIntermediateState: state = kFSIState; break; case genie::kIStFinalStateNuclearRemnant: case genie::kIStUndefined: case genie::kIStDecayedState: default: break; } // Flag to remove nuclear part in genie if (p->Pdg() > 1000000) { if (state == kInitialState) state = kNuclearInitial; else if (state == kFinalState) state = kNuclearRemnant; } return state; } #endif #ifdef __GENIE_ENABLED__ int GENIEInputHandler::ConvertGENIEReactionCode(GHepRecord* gheprec) { // Electron Scattering if (gheprec->Summary()->ProcInfo().IsEM()) { if (gheprec->Summary()->InitState().ProbePdg() == 11) { if (gheprec->Summary()->ProcInfo().IsQuasiElastic()) return 1; else if (gheprec->Summary()->ProcInfo().IsMEC()) return 2; else if (gheprec->Summary()->ProcInfo().IsResonant()) return 13; else if (gheprec->Summary()->ProcInfo().IsDeepInelastic()) return 26; else { ERROR(WRN, "Unknown GENIE Electron Scattering Mode!" << std::endl << "ScatteringTypeId = " << gheprec->Summary()->ProcInfo().ScatteringTypeId() << " " << "InteractionTypeId = " << gheprec->Summary()->ProcInfo().InteractionTypeId() << std::endl << genie::ScatteringType::AsString(gheprec->Summary()->ProcInfo().ScatteringTypeId()) << " " << genie::InteractionType::AsString(gheprec->Summary()->ProcInfo().InteractionTypeId()) << " " << gheprec->Summary()->ProcInfo().IsMEC()); return 0; } } // Weak CC } else if (gheprec->Summary()->ProcInfo().IsWeakCC()) { // CC MEC if (gheprec->Summary()->ProcInfo().IsMEC()) { if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg())) return 2; else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg())) return -2; // CC OTHER } else { return utils::ghep::NeutReactionCode(gheprec); } // Weak NC } else if (gheprec->Summary()->ProcInfo().IsWeakNC()) { // NC MEC if (gheprec->Summary()->ProcInfo().IsMEC()) { if (pdg::IsNeutrino(gheprec->Summary()->InitState().ProbePdg())) return 32; else if (pdg::IsAntiNeutrino(gheprec->Summary()->InitState().ProbePdg())) return -32; // NC OTHER } else { return utils::ghep::NeutReactionCode(gheprec); } } return 0; } void GENIEInputHandler::CalcNUISANCEKinematics() { // Reset all variables fNUISANCEEvent->ResetEvent(); // Check for GENIE Event if (!fGenieNtpl) return; if (!fGenieNtpl->event) return; // Cast Event Record fGenieGHep = static_cast(fGenieNtpl->event); if (!fGenieGHep) return; // Convert GENIE Reaction Code fNUISANCEEvent->fMode = ConvertGENIEReactionCode(fGenieGHep); // Set Event Info fNUISANCEEvent->Mode = fNUISANCEEvent->fMode; fNUISANCEEvent->fEventNo = 0.0; fNUISANCEEvent->fTotCrs = fGenieGHep->XSec(); fNUISANCEEvent->fTargetA = 0.0; fNUISANCEEvent->fTargetZ = 0.0; fNUISANCEEvent->fTargetH = 0; fNUISANCEEvent->fBound = 0.0; fNUISANCEEvent->InputWeight = 1.0; //(1E+38 / genie::units::cm2) * fGenieGHep->XSec(); // Get N Particle Stack unsigned int npart = fGenieGHep->GetEntries(); unsigned int kmax = fNUISANCEEvent->kMaxParticles; if (npart > kmax) { ERR(WRN) << "GENIE has too many particles, expanding stack." << std::endl; fNUISANCEEvent->ExpandParticleStack(npart); } // Fill Particle Stack GHepParticle* p = 0; TObjArrayIter iter(fGenieGHep); fNUISANCEEvent->fNParticles = 0; // Loop over all particles while ((p = (dynamic_cast((iter).Next())))) { if (!p) continue; // Get Status int state = GetGENIEParticleStatus(p, fNUISANCEEvent->fMode); // Remove Undefined if (kRemoveUndefParticles && state == kUndefinedState) continue; // Remove FSI if (kRemoveFSIParticles && state == kFSIState) continue; if (kRemoveNuclearParticles && (state == kNuclearInitial || state == kNuclearRemnant)) continue; // Fill Vectors int curpart = fNUISANCEEvent->fNParticles; fNUISANCEEvent->fParticleState[curpart] = state; // Mom fNUISANCEEvent->fParticleMom[curpart][0] = p->Px() * 1.E3; fNUISANCEEvent->fParticleMom[curpart][1] = p->Py() * 1.E3; fNUISANCEEvent->fParticleMom[curpart][2] = p->Pz() * 1.E3; fNUISANCEEvent->fParticleMom[curpart][3] = p->E() * 1.E3; // PDG fNUISANCEEvent->fParticlePDG[curpart] = p->Pdg(); // Add to N particle count fNUISANCEEvent->fNParticles++; // Extra Check incase GENIE fails. if ((UInt_t)fNUISANCEEvent->fNParticles == kmax) { ERR(WRN) << "Number of GENIE Particles exceeds maximum!" << std::endl; ERR(WRN) << "Extend kMax, or run without including FSI particles!" << std::endl; break; } } // Fill Extra Stack if (fSaveExtra) fGenieInfo->FillGeneratorInfo(fGenieNtpl); // 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 GENIEInputHandler::Print() { } #endif diff --git a/src/MINERvA/CMakeLists.txt b/src/MINERvA/CMakeLists.txt index 91889ab..d41a492 100644 --- a/src/MINERvA/CMakeLists.txt +++ b/src/MINERvA/CMakeLists.txt @@ -1,153 +1,156 @@ # 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(IMPLFILES MINERvA_CCQE_XSec_1DQ2_antinu.cxx MINERvA_CCQE_XSec_1DQ2_joint.cxx MINERvA_CCQE_XSec_1DQ2_nu.cxx MINERvA_CC0pi_XSec_1DEe_nue.cxx MINERvA_CC0pi_XSec_1DQ2_nue.cxx MINERvA_CC0pi_XSec_1DQ2_nu_proton.cxx MINERvA_CC0pi_XSec_1DThetae_nue.cxx MINERvA_CC1pi0_XSec_1DEnu_antinu.cxx MINERvA_CC1pi0_XSec_1DQ2_antinu.cxx MINERvA_CC1pi0_XSec_1Dpmu_antinu.cxx MINERvA_CC1pi0_XSec_1Dppi0_antinu.cxx MINERvA_CC1pi0_XSec_1DTpi0_antinu.cxx MINERvA_CC1pi0_XSec_1Dth_antinu.cxx MINERvA_CC1pi0_XSec_1Dthmu_antinu.cxx MINERvA_CC1pip_XSec_1DTpi_20deg_nu.cxx MINERvA_CC1pip_XSec_1DTpi_nu.cxx MINERvA_CC1pip_XSec_1Dth_20deg_nu.cxx MINERvA_CC1pip_XSec_1Dth_nu.cxx MINERvA_CCNpip_XSec_1DEnu_nu.cxx MINERvA_CCNpip_XSec_1DQ2_nu.cxx MINERvA_CCNpip_XSec_1DTpi_nu.cxx MINERvA_CCNpip_XSec_1Dpmu_nu.cxx MINERvA_CCNpip_XSec_1Dth_nu.cxx MINERvA_CCNpip_XSec_1Dthmu_nu.cxx MINERvA_CCinc_XSec_2DEavq3_nu.cxx MINERvA_CCinc_XSec_1Dx_ratio.cxx MINERvA_CCinc_XSec_1DEnu_ratio.cxx MINERvA_CCinc_XSec_1Dx_nu.cxx MINERvA_CCinc_XSec_1DEnu_nu.cxx MINERvA_CCDIS_XSec_1Dx_ratio.cxx MINERvA_CCDIS_XSec_1DEnu_ratio.cxx MINERvA_CCDIS_XSec_1Dx_nu.cxx MINERvA_CCDIS_XSec_1DEnu_nu.cxx MINERvA_CC0pi_XSec_1DQ2_Tgt_nu.cxx MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.cxx MINERvA_CC1pip_XSec_1D_2017Update.cxx MINERvA_CCCOHPI_XSec_1DEnu_nu.cxx MINERvA_CCCOHPI_XSec_1DEpi_nu.cxx MINERvA_CCCOHPI_XSec_1Dth_nu.cxx +MINERvA_CCCOHPI_XSec_1DQ2_nu.cxx + MINERvA_CCCOHPI_XSec_1DEnu_antinu.cxx MINERvA_CCCOHPI_XSec_1DEpi_antinu.cxx MINERvA_CCCOHPI_XSec_1Dth_antinu.cxx +MINERvA_CCCOHPI_XSec_1DQ2_antinu.cxx MINERvAUtils.cxx MINERvA_SignalDef.cxx ) set(HEADERFILES MINERvA_CCQE_XSec_1DQ2_antinu.h MINERvA_CCQE_XSec_1DQ2_joint.h MINERvA_CCQE_XSec_1DQ2_nu.h MINERvA_CC0pi_XSec_1DEe_nue.h MINERvA_CC0pi_XSec_1DQ2_nue.h MINERvA_CC0pi_XSec_1DQ2_nu_proton.h MINERvA_CC0pi_XSec_1DThetae_nue.h MINERvA_CC1pi0_XSec_1DEnu_antinu.h MINERvA_CC1pi0_XSec_1DQ2_antinu.h MINERvA_CC1pi0_XSec_1Dpmu_antinu.h MINERvA_CC1pi0_XSec_1Dppi0_antinu.h MINERvA_CC1pi0_XSec_1DTpi0_antinu.h MINERvA_CC1pi0_XSec_1Dth_antinu.h MINERvA_CC1pi0_XSec_1Dthmu_antinu.h MINERvA_CC1pip_XSec_1DTpi_20deg_nu.h MINERvA_CC1pip_XSec_1DTpi_nu.h MINERvA_CC1pip_XSec_1Dth_20deg_nu.h MINERvA_CC1pip_XSec_1Dth_nu.h MINERvA_CCNpip_XSec_1DEnu_nu.h MINERvA_CCNpip_XSec_1DQ2_nu.h MINERvA_CCNpip_XSec_1DTpi_nu.h MINERvA_CCNpip_XSec_1Dpmu_nu.h MINERvA_CCNpip_XSec_1Dth_nu.h MINERvA_CCNpip_XSec_1Dthmu_nu.h MINERvA_CCinc_XSec_2DEavq3_nu.h MINERvA_CCinc_XSec_1Dx_ratio.h MINERvA_CCinc_XSec_1DEnu_ratio.h MINERvA_CCinc_XSec_1Dx_nu.h MINERvA_CCinc_XSec_1DEnu_nu.h MINERvA_CCDIS_XSec_1Dx_ratio.h MINERvA_CCDIS_XSec_1DEnu_ratio.h MINERvA_CCDIS_XSec_1Dx_nu.h MINERvA_CCDIS_XSec_1DEnu_nu.h MINERvA_CC0pi_XSec_1DQ2_Tgt_nu.h MINERvA_CC0pi_XSec_1DQ2_TgtRatio_nu.h MINERvA_CC1pip_XSec_1D_2017Update.h MINERvA_CCCOHPI_XSec_1DEnu_nu.h MINERvA_CCCOHPI_XSec_1DEpi_nu.h MINERvA_CCCOHPI_XSec_1Dth_nu.h MINERvA_CCCOHPI_XSec_1DEnu_antinu.h MINERvA_CCCOHPI_XSec_1DEpi_antinu.h MINERvA_CCCOHPI_XSec_1Dth_antinu.h MINERvAUtils.h MINERvA_SignalDef.h MINERvAVariableBoxes.h ) set(LIBNAME expMINERvA) if(CMAKE_BUILD_TYPE MATCHES DEBUG) add_library(${LIBNAME} STATIC ${IMPLFILES}) else(CMAKE_BUILD_TYPE MATCHES RELEASE) add_library(${LIBNAME} SHARED ${IMPLFILES}) endif() include_directories(${MINIMUM_INCLUDE_DIRECTORIES}) set_target_properties(${LIBNAME} PROPERTIES VERSION "${NUISANCE_VERSION_MAJOR}.${NUISANCE_VERSION_MINOR}.${NUISANCE_VERSION_REVISION}") #set_target_properties(${LIBNAME} PROPERTIES LINK_FLAGS ${ROOT_LD_FLAGS}) if(DEFINED PROJECTWIDE_EXTRA_DEPENDENCIES) add_dependencies(${LIBNAME} ${PROJECTWIDE_EXTRA_DEPENDENCIES}) endif() install(TARGETS ${LIBNAME} DESTINATION lib) #Can uncomment this to install the headers... but is it really neccessary? #install(FILES ${HEADERFILES} DESTINATION include) set(MODULETargets ${MODULETargets} ${LIBNAME} PARENT_SCOPE) diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DQ2_antinu.cxx b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DQ2_antinu.cxx new file mode 100644 index 0000000..0cc2838 --- /dev/null +++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DQ2_antinu.cxx @@ -0,0 +1,84 @@ +// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret + +/******************************************************************************* +* This file is part of NUISANCE. +* +* NUISANCE is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* NUISANCE is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with NUISANCE. If not, see . +*******************************************************************************/ +#include "MINERvA_SignalDef.h" +#include "MINERvA_CCCOHPI_XSec_1DQ2_antinu.h" + + +//******************************************************************** +MINERvA_CCCOHPI_XSec_1DQ2_antinu::MINERvA_CCCOHPI_XSec_1DQ2_antinu(nuiskey samplekey) { +//******************************************************************** + + // Sample overview --------------------------------------------------- + std::string descrip = "MINERvA_CCCOHPI_XSec_1DQ2_antinu sample. \n" \ + "Target: CH \n" \ + "Flux: MINERvA Reverse Horn Current numu \n" \ + "Signal: Any event with 1 mu+, 1pi-, and no other FS particles \n"; + + // Setup common settings + fSettings = LoadSampleSettings(samplekey); + fSettings.SetDescription(descrip); + fSettings.SetXTitle("E_{#pi} (GeV)"); + fSettings.SetYTitle("d#sigma/dE_{#pi} (cm^{2}/GeV/C^{12})"); + fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL"); + fSettings.SetEnuRange(1.5, 20.0); + fSettings.DefineAllowedTargets("C,H"); + fSettings.DefineAllowedSpecies("numub"); + fSettings.SetTitle("MINERvA_CCCOHPI_XSec_1DQ2_antinu"); + + fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Q2_nubar_data.csv"); + fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Q2_nubar_cov.csv"); + + FinaliseSampleSettings(); + + // Scaling Setup --------------------------------------------------- + // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon + fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents)/TotalIntegratedFlux("width")*13; + + // Plot Setup ------------------------------------------------------- + SetDataFromTextFile( fSettings.GetDataInput() ); + SetCovarFromMultipleTextFiles(fSettings.GetCovarInput()); + + // Final setup --------------------------------------------------- + FinaliseMeasurement(); + +}; + + +void MINERvA_CCCOHPI_XSec_1DQ2_antinu::FillEventVariables(FitEvent *event) { + + if (event->NumFSParticle(-211) == 0 || + event->NumISParticle(-14) == 0 || + event->NumFSParticle(-13) == 0) return; + + TLorentzVector Ppi = event->GetHMFSParticle(-211)->fP; + TLorentzVector Pnu = event->GetHMISParticle(-14)->fP; + TLorentzVector Pmu = event->GetHMFSParticle(-13)->fP; + double Q2 = (Pmu - Pnu).Mag2()/1.E6; + fXVar = Q2; + + return; +}; + +//******************************************************************** +bool MINERvA_CCCOHPI_XSec_1DQ2_antinu::isSignal(FitEvent *event) { +//******************************************************************** + // Signal: numu + 12C --> mu- + pi+ + 12C' + return SignalDef::isCCCOH(event, -14, -211, EnuMin, EnuMax); +} + diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DQ2_antinu.h b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DQ2_antinu.h new file mode 100644 index 0000000..dde6ce6 --- /dev/null +++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DQ2_antinu.h @@ -0,0 +1,36 @@ +// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret + +/******************************************************************************* +* This file is part of NUISANCE. +* +* NUISANCE is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* NUISANCE is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with NUISANCE. If not, see . +*******************************************************************************/ + +#ifndef MINERVA_CCCOHPI_XSEC_1DQ2_ANTINU_H_SEEN +#define MINERVA_CCCOHPI_XSEC_1DQ2_ANTINU_H_SEEN + +#include "Measurement1D.h" + +class MINERvA_CCCOHPI_XSec_1DQ2_antinu : public Measurement1D { +public: + MINERvA_CCCOHPI_XSec_1DQ2_antinu(nuiskey samplekey); + virtual ~MINERvA_CCCOHPI_XSec_1DQ2_antinu() {}; + + void FillEventVariables(FitEvent *event); + bool isSignal(FitEvent *event); + + private: +}; + +#endif diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DQ2_nu.cxx b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DQ2_nu.cxx new file mode 100644 index 0000000..fad6408 --- /dev/null +++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DQ2_nu.cxx @@ -0,0 +1,82 @@ +// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret + +/******************************************************************************* +* This file is part of NUISANCE. +* +* NUISANCE is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* NUISANCE is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with NUISANCE. If not, see . +*******************************************************************************/ +#include "MINERvA_SignalDef.h" +#include "MINERvA_CCCOHPI_XSec_1DQ2_nu.h" + +//******************************************************************** +MINERvA_CCCOHPI_XSec_1DQ2_nu::MINERvA_CCCOHPI_XSec_1DQ2_nu(nuiskey samplekey) { +//******************************************************************** + + // Sample overview --------------------------------------------------- + std::string descrip = "MINERvA_CCCOHPI_XSec_1DQ2_nu sample. \n" \ + "Target: CH \n" \ + "Flux: MINERvA Forward Horn Current numu \n" \ + "Signal: Any event with 1 mu-, 1pi+, and no other FS particles \n"; + + // Setup common settings + fSettings = LoadSampleSettings(samplekey); + fSettings.SetDescription(descrip); + fSettings.SetXTitle("E_{#pi} (GeV)"); + fSettings.SetYTitle("d#sigma/dE_{#pi} (cm^{2}/GeV/C^{12})"); + fSettings.SetAllowedTypes("FIX,FREE,SHAPE/DIAG,FULL/NORM/MASK", "FIX/FULL"); + fSettings.SetEnuRange(1.5, 20.0); + fSettings.DefineAllowedTargets("C,H"); + fSettings.DefineAllowedSpecies("numu"); + fSettings.SetTitle("MINERvA_CCCOHPI_XSec_1DQ2_nu"); + + fSettings.SetDataInput( GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Q2_nu_data.csv"); + fSettings.SetCovarInput(GeneralUtils::GetTopLevelDir() + "/data/MINERvA/CCcoh/Q2_nu_cov.csv"); + + FinaliseSampleSettings(); + + // Scaling Setup --------------------------------------------------- + // ScaleFactor automatically setup for DiffXSec/cm2/Nucleon + fScaleFactor = GetEventHistogram()->Integral("width") * double(1E-38) / double(fNEvents)/TotalIntegratedFlux("width")*13; + + // Plot Setup ------------------------------------------------------- + SetDataFromTextFile( fSettings.GetDataInput() ); + SetCovarFromMultipleTextFiles(fSettings.GetCovarInput()); + + // Final setup --------------------------------------------------- + FinaliseMeasurement(); +}; + + +void MINERvA_CCCOHPI_XSec_1DQ2_nu::FillEventVariables(FitEvent *event) { + + if (event->NumFSParticle(211) == 0 || + event->NumISParticle(14) == 0 || + event->NumFSParticle(13) == 0) return; + + TLorentzVector Ppi = event->GetHMFSParticle(211)->fP; + TLorentzVector Pnu = event->GetHMISParticle(14)->fP; + TLorentzVector Pmu = event->GetHMFSParticle(13)->fP; + double Q2 = (Pmu - Pnu).Mag2()/1.E6; + fXVar = Q2; + + return; +}; + +//******************************************************************** +bool MINERvA_CCCOHPI_XSec_1DQ2_nu::isSignal(FitEvent *event) { +//******************************************************************** + // Signal: numu + 12C --> mu- + pi+ + 12C' + return SignalDef::isCCCOH(event, 14, 211, EnuMin, EnuMax); +} + diff --git a/src/MINERvA/MINERvA_CCCOHPI_XSec_1DQ2_nu.h b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DQ2_nu.h new file mode 100644 index 0000000..d22b181 --- /dev/null +++ b/src/MINERvA/MINERvA_CCCOHPI_XSec_1DQ2_nu.h @@ -0,0 +1,36 @@ +// Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret + +/******************************************************************************* +* This file is part of NUISANCE. +* +* NUISANCE is free software: you can redistribute it and/or modify +* it under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* NUISANCE is distributed in the hope that it will be useful, +* but WITHOUT ANY WARRANTY; without even the implied warranty of +* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +* GNU General Public License for more details. +* +* You should have received a copy of the GNU General Public License +* along with NUISANCE. If not, see . +*******************************************************************************/ + +#ifndef MINERVA_CCCOHPI_XSEC_1DQ2_NU_H_SEEN +#define MINERVA_CCCOHPI_XSEC_1DQ2_NU_H_SEEN + +#include "Measurement1D.h" + +class MINERvA_CCCOHPI_XSec_1DQ2_nu : public Measurement1D { +public: + MINERvA_CCCOHPI_XSec_1DQ2_nu(nuiskey samplekey); + virtual ~MINERvA_CCCOHPI_XSec_1DQ2_nu() {}; + + void FillEventVariables(FitEvent *event); + bool isSignal(FitEvent *event); + + private: +}; + +#endif diff --git a/src/Reweight/CMakeLists.txt b/src/Reweight/CMakeLists.txt index 1f9270c..685133c 100644 --- a/src/Reweight/CMakeLists.txt +++ b/src/Reweight/CMakeLists.txt @@ -1,79 +1,81 @@ # Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret ################################################################################ # This file is part of NUISANCE. # # NUISANCE is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # NUISANCE is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with NUISANCE. If not, see . ################################################################################ set(IMPLFILES GlobalDialList.cxx FitWeight.cxx WeightEngineBase.cxx NEUTWeightEngine.cxx NuWroWeightEngine.cxx GENIEWeightEngine.cxx WeightUtils.cxx SampleNormEngine.cxx LikelihoodWeightEngine.cxx SplineWeightEngine.cxx NUISANCESyst.cxx T2KWeightEngine.cxx NUISANCEWeightEngine.cxx NUISANCEWeightCalcs.cxx NIWGWeightEngine.cxx OscWeightEngine.cxx +MINERvAWeightCalcs.cxx ) set(HEADERFILES GlobalDialList.h FitWeight.h WeightEngineBase.h NEUTWeightEngine.h NuWroWeightEngine.h GENIEWeightEngine.h WeightUtils.h SampleNormEngine.h LikelihoodWeightEngine.h SplineWeightEngine.h NUISANCESyst.h T2KWeightEngine.h NUISANCEWeightEngine.h NUISANCEWeightCalcs.h NIWGWeightEngine.h OscWeightEngine.h +MINERvAWeightCalcs.h ) set(LIBNAME Reweight) if(CMAKE_BUILD_TYPE MATCHES DEBUG) add_library(${LIBNAME} STATIC ${IMPLFILES}) else(CMAKE_BUILD_TYPE MATCHES RELEASE) add_library(${LIBNAME} SHARED ${IMPLFILES}) endif() include_directories(${MINIMUM_INCLUDE_DIRECTORIES}) set_target_properties(${LIBNAME} PROPERTIES VERSION "${ExtFit_VERSION_MAJOR}.${ExtFit_VERSION_MINOR}.${ExtFit_VERSION_REVISION}") #set_target_properties(${LIBNAME} PROPERTIES LINK_FLAGS ${ROOT_LD_FLAGS}) if(DEFINED PROJECTWIDE_EXTRA_DEPENDENCIES) add_dependencies(${LIBNAME} ${PROJECTWIDE_EXTRA_DEPENDENCIES}) endif() install(TARGETS ${LIBNAME} DESTINATION lib) #Can uncomment this to install the headers... but is it really neccessary? install(FILES ${HEADERFILES} DESTINATION include) set(MODULETargets ${MODULETargets} ${LIBNAME} PARENT_SCOPE) diff --git a/src/Reweight/GENIEWeightEngine.cxx b/src/Reweight/GENIEWeightEngine.cxx index 19d33f8..4b4ceb1 100644 --- a/src/Reweight/GENIEWeightEngine.cxx +++ b/src/Reweight/GENIEWeightEngine.cxx @@ -1,234 +1,231 @@ #include "GENIEWeightEngine.h" GENIEWeightEngine::GENIEWeightEngine(std::string name) { #ifdef __GENIE_ENABLED__ // Setup the NEUT Reweight engien fCalcName = name; LOG(FIT) << "Setting up GENIE RW : " << fCalcName << std::endl; // Create RW Engine suppressing cout StopTalking(); fGenieRW = new genie::rew::GReWeight(); // Get List of Vetos (Just for debugging) std::string rw_engine_list = FitPar::Config().GetParS("FitWeight.fGenieRW_veto"); bool xsec_ncel = rw_engine_list.find("xsec_ncel") == std::string::npos; bool xsec_ccqe = rw_engine_list.find("xsec_ccqe") == std::string::npos; bool xsec_coh = rw_engine_list.find("xsec_coh") == std::string::npos; bool xsec_nnres = rw_engine_list.find("xsec_nonresbkg") == std::string::npos; bool xsec_nudis = rw_engine_list.find("nuclear_dis") == std::string::npos; bool xsec_resdec = rw_engine_list.find("hadro_res_decay") == std::string::npos; bool xsec_fzone = rw_engine_list.find("hadro_intranuke") == std::string::npos; bool xsec_intra = rw_engine_list.find("hadro_fzone") == std::string::npos; bool xsec_agky = rw_engine_list.find("hadro_agky") == std::string::npos; bool xsec_qevec = rw_engine_list.find("xsec_ccqe_vec") == std::string::npos; bool xsec_dis = rw_engine_list.find("xsec_dis") == std::string::npos; bool xsec_nc = rw_engine_list.find("xsec_nc") == std::string::npos; bool xsec_ccres = rw_engine_list.find("xsec_ccres") == std::string::npos; bool xsec_ncres = rw_engine_list.find("xsec_ncres") == std::string::npos; bool xsec_nucqe = rw_engine_list.find("nuclear_qe") == std::string::npos; bool xsec_qeaxial = rw_engine_list.find("xsec_ccqe_axial") == std::string::npos; // Now actually add the RW Calcs if (xsec_ncel) fGenieRW->AdoptWghtCalc("xsec_ncel", new genie::rew::GReWeightNuXSecNCEL); if (xsec_ccqe) fGenieRW->AdoptWghtCalc("xsec_ccqe", new genie::rew::GReWeightNuXSecCCQE); if (xsec_coh) fGenieRW->AdoptWghtCalc("xsec_coh", new genie::rew::GReWeightNuXSecCOH); if (xsec_nnres) fGenieRW->AdoptWghtCalc("xsec_nonresbkg", new genie::rew::GReWeightNonResonanceBkg); if (xsec_nudis) fGenieRW->AdoptWghtCalc("nuclear_dis", new genie::rew::GReWeightDISNuclMod); if (xsec_resdec) fGenieRW->AdoptWghtCalc("hadro_res_decay", new genie::rew::GReWeightResonanceDecay); if (xsec_fzone) fGenieRW->AdoptWghtCalc("hadro_fzone", new genie::rew::GReWeightFZone); if (xsec_intra) fGenieRW->AdoptWghtCalc("hadro_intranuke", new genie::rew::GReWeightINuke); if (xsec_agky) fGenieRW->AdoptWghtCalc("hadro_agky", new genie::rew::GReWeightAGKY); if (xsec_qevec) fGenieRW->AdoptWghtCalc("xsec_ccqe_vec", new genie::rew::GReWeightNuXSecCCQEvec); #if __GENIE_VERSION__ >= 212 if (xsec_qeaxial) fGenieRW->AdoptWghtCalc("xsec_ccqe_axial", new genie::rew::GReWeightNuXSecCCQEaxial); #endif if (xsec_dis) fGenieRW->AdoptWghtCalc("xsec_dis", new genie::rew::GReWeightNuXSecDIS); if (xsec_nc) fGenieRW->AdoptWghtCalc("xsec_nc", new genie::rew::GReWeightNuXSecNC); if (xsec_ccres) fGenieRW->AdoptWghtCalc("xsec_ccres", new genie::rew::GReWeightNuXSecCCRES); if (xsec_ncres) fGenieRW->AdoptWghtCalc("xsec_ncres", new genie::rew::GReWeightNuXSecNCRES); if (xsec_nucqe) fGenieRW->AdoptWghtCalc("nuclear_qe", new genie::rew::GReWeightFGM); GReWeightNuXSecCCQE * rwccqe = dynamic_cast (fGenieRW->WghtCalc("xsec_ccqe")); rwccqe->SetMode(GReWeightNuXSecCCQE::kModeMa); // Default to include shape and normalization changes for CCRES (can be changed downstream if desired) GReWeightNuXSecCCRES * rwccres = dynamic_cast (fGenieRW->WghtCalc("xsec_ccres")); std::string marestype = FitPar::Config().GetParS("GENIEWeightEngine_CCRESMode"); if (!marestype.compare("kModeNormAndMaMvShape")){ rwccres->SetMode(GReWeightNuXSecCCRES::kModeNormAndMaMvShape); } else if (!marestype.compare("kModeMaMv")){ rwccres->SetMode(GReWeightNuXSecCCRES::kModeMaMv); } else { THROW("Unkown MARES Mode in GENIE Weight Engine : " << marestype ); } // Default to include shape and normalization changes for NCRES (can be changed downstream if desired) GReWeightNuXSecNCRES * rwncres = dynamic_cast (fGenieRW->WghtCalc("xsec_ncres")); rwncres->SetMode(GReWeightNuXSecNCRES::kModeMaMv); // Default to include shape and normalization changes for DIS (can be changed downstream if desired) GReWeightNuXSecDIS * rwdis = dynamic_cast (fGenieRW->WghtCalc("xsec_dis")); rwdis->SetMode(GReWeightNuXSecDIS::kModeABCV12u); - // Final Reconfigure - fGenieRW->Reconfigure(); - // Set Abs Twk Config fIsAbsTwk = (FitPar::Config().GetParB("setabstwk")); // allow cout again StartTalking(); #else ERR(FTL) << "GENIE RW NOT ENABLED" << std::endl; #endif }; void GENIEWeightEngine::IncludeDial(std::string name, double startval) { #ifdef __GENIE_ENABLED__ // Get First enum int nuisenum = Reweight::ConvDial(name, kGENIE); // Setup Maps fEnumIndex[nuisenum];// = std::vector(0); fNameIndex[name]; // = std::vector(0); // Split by commas std::vector allnames = GeneralUtils::ParseToStr(name, ","); for (uint i = 0; i < allnames.size(); i++) { std::string singlename = allnames[i]; // Get RW genie::rew::GSyst_t rwsyst = GSyst::FromString(singlename); // Fill Maps int index = fValues.size(); fValues.push_back(0.0); fGENIESysts.push_back(rwsyst); // Initialize dial std::cout << "Registering " << singlename << " from " << name << std::endl; fGenieRW->Systematics().Init( fGENIESysts[index] ); // If Absolute if (fIsAbsTwk) { GSystUncertainty::Instance()->SetUncertainty( rwsyst, 1.0, 1.0 ); } // Setup index fEnumIndex[nuisenum].push_back(index); fNameIndex[name].push_back(index); } // Set Value if given if (startval != -999.9) { SetDialValue(nuisenum, startval); } #endif }; void GENIEWeightEngine::SetDialValue(int nuisenum, double val) { #ifdef __GENIE_ENABLED__ std::vector indices = fEnumIndex[nuisenum]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; fGenieRW->Systematics().Set( fGENIESysts[indices[i]], val); } #endif } void GENIEWeightEngine::SetDialValue(std::string name, double val) { #ifdef __GENIE_ENABLED__ std::vector indices = fNameIndex[name]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; fGenieRW->Systematics().Set(fGENIESysts[indices[i]], val); } #endif } void GENIEWeightEngine::Reconfigure(bool silent) { #ifdef __GENIE_ENABLED__ // Hush now... if (silent) StopTalking(); // Reconf fGenieRW->Reconfigure(); fGenieRW->Print(); // Shout again if (silent) StartTalking(); #endif } double GENIEWeightEngine::CalcWeight(BaseFitEvt* evt) { double rw_weight = 1.0; #ifdef __GENIE_ENABLED__ // Skip Non GENIE if (evt->fType != kGENIE) return 1.0; // Make nom weight if (!evt) { THROW("evt not found : " << evt); } if (!(evt->genie_event)) { THROW("evt->genie_event not found!" << evt->genie_event); } if (!(evt->genie_event->event)) { THROW("evt->genie_event->event GHepRecord not found!" << (evt->genie_event->event)); } if (!fGenieRW) { THROW("GENIE RW Not Found!" << fGenieRW); } rw_weight = fGenieRW->CalcWeight(*(evt->genie_event->event)); // std::cout << "Returning GENIE Weight for electron scattering = " << rw_weight << std::endl; #endif // Return rw_weight return rw_weight; } diff --git a/src/Reweight/MINERvAWeightCalcs.cxx b/src/Reweight/MINERvAWeightCalcs.cxx new file mode 100644 index 0000000..ce70e83 --- /dev/null +++ b/src/Reweight/MINERvAWeightCalcs.cxx @@ -0,0 +1,403 @@ +#ifdef __MINERVA_RW_ENABLED__ +#ifdef __GENIE_ENABLED__ + +#include "MINERvAWeightCalcs.h" +#include "BaseFitEvt.h" + +namespace nuisance { +namespace reweight { + +//******************************************************* +MINERvAReWeight_MEC::MINERvAReWeight_MEC() { + //******************************************************* + fTwk_NormCCMEC = 0.0; + fDef_NormCCMEC = 1.0; + fCur_NormCCMEC = fDef_NormCCMEC; +} + +MINERvAReWeight_MEC::~MINERvAReWeight_MEC(){}; + +double MINERvAReWeight_MEC::CalcWeight(BaseFitEvt* evt) { + // Check GENIE + if (evt->fType != kGENIE) return 1.0; + + // Extract the GENIE Record + GHepRecord* ghep = static_cast(evt->genie_event->event); + const Interaction* interaction = ghep->Summary(); + const InitialState& init_state = interaction->InitState(); + const ProcessInfo& proc_info = interaction->ProcInfo(); + const Target& tgt = init_state.Tgt(); + + // If the event is not MEC this Calc doesn't handle it + if (!proc_info.IsMEC()) return 1.0; + + // WEIGHT CALCULATIONS ------------- + double w = 1.0; + + // CCMEC Dial + if (!proc_info.IsWeakCC()) w *= fCur_NormCCMEC; + + // Return Combined Weight + return w; +} + +void MINERvAReWeight_MEC::SetDialValue(std::string name, double val) { + SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); +} + +void MINERvAReWeight_MEC::SetDialValue(int rwenum, double val) { + // Check Handled + int curenum = rwenum % 1000; + if (!IsHandled(curenum)) return; + + // Set Values + if (curenum == kMINERvARW_NormCCMEC) { + fTwk_NormCCMEC = val; + fCur_NormCCMEC = fDef_NormCCMEC + fTwk_NormCCMEC; + } + + // Define Tweaked + fTweaked = ((fTwk_NormCCMEC != 0.0)); +} + +bool MINERvAReWeight_MEC::IsHandled(int rwenum) { + int curenum = rwenum % 1000; + + switch (curenum) { + case kMINERvARW_NormCCMEC: + return true; + default: + return false; + } +} + +//******************************************************* +MINERvAReWeight_RES::MINERvAReWeight_RES() { + //******************************************************* + fTwk_NormCCRES = 0.0; + fDef_NormCCRES = 1.0; + fCur_NormCCRES = fDef_NormCCRES; +} + +MINERvAReWeight_RES::~MINERvAReWeight_RES(){}; + +double MINERvAReWeight_RES::CalcWeight(BaseFitEvt* evt) { + + // std::cout << "Caculating RES" << std::endl; + // Check GENIE + if (evt->fType != kGENIE) return 1.0; + + // Extract the GENIE Record + GHepRecord* ghep = static_cast(evt->genie_event->event); + const Interaction* interaction = ghep->Summary(); + const InitialState& init_state = interaction->InitState(); + const ProcessInfo& proc_info = interaction->ProcInfo(); + const Target& tgt = init_state.Tgt(); + + // If the event is not RES this Calc doesn't handle it + if (!proc_info.IsResonant()) return 1.0; + + // WEIGHT CALCULATIONS ------------- + double w = 1.0; + + // CCRES Dial + if (proc_info.IsWeakCC()) w *= fCur_NormCCRES; + + // Return Combined Weight + return w; +} + +void MINERvAReWeight_RES::SetDialValue(std::string name, double val) { + SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); +} + +void MINERvAReWeight_RES::SetDialValue(int rwenum, double val) { + // Check Handled + int curenum = rwenum % 1000; + if (!IsHandled(curenum)) return; + + // Set Values + if (curenum == kMINERvARW_NormCCRES) { + fTwk_NormCCRES = val; + fCur_NormCCRES = fDef_NormCCRES + fTwk_NormCCRES; + } + + // Define Tweaked + fTweaked = ((fTwk_NormCCRES != 0.0)); +} + +bool MINERvAReWeight_RES::IsHandled(int rwenum) { + int curenum = rwenum % 1000; + + switch (curenum) { + case kMINERvARW_NormCCRES: + return true; + default: + return false; + } +} + +//******************************************************* +RikRPA::RikRPA() { + //******************************************************* + + // - Syst : kMINERvA_RikRPA_ApplyRPA + // - Type : Binary + // - Limits : 0.0 (false) -> 1.0 (true) + // - Default : 0.0 + fApplyDial_RPACorrection = false; + + // - Syst : kMINERvA_RikRPA_LowQ2 + // - Type : Absolute + // - Limits : 1.0 -> 1.0 + // - Default : 0.0 + // - Frac Error : 100% + fDefDial_RPALowQ2 = 0.0; + fCurDial_RPALowQ2 = fDefDial_RPALowQ2; + fErrDial_RPALowQ2 = 0.0; + + // - Syst : kMINERvA_RikRPA_HighQ2 + // - Type : Absolute + // - Limits : 1.0 -> 1.0 + // - Default : 0.0 + // - Frac Error : 100% + fDefDial_RPAHighQ2 = 0.0; + fCurDial_RPAHighQ2 = fDefDial_RPAHighQ2; + fErrDial_RPAHighQ2 = 1.0; + + fEventWeights = new double[5]; + + for (size_t i = 0; i < kMaxCalculators; i++) { + // fRPACalculators[i] = NULL; + } + + fTweaked = false; +} + +RikRPA::~RikRPA() { + delete fEventWeights; + + for (size_t i = 0; i < kMaxCalculators; i++) { + // if (fRPACalculators[i]) delete fRPACalculators[i]; + // fRPACalculators[i] = NULL; + } +} + +double RikRPA::CalcWeight(BaseFitEvt* evt) { + // LOG(FIT) << "Calculating RikRPA" << std::endl; + // Return 1.0 if not tweaked + if (!fTweaked) return 1.0; + + double w = 1.0; + + // Extract the GENIE Record + GHepRecord* ghep = static_cast(evt->genie_event->event); + const Interaction* interaction = ghep->Summary(); + const InitialState& init_state = interaction->InitState(); + const ProcessInfo& proc_info = interaction->ProcInfo(); + // const Kinematics & kine = interaction->Kine(); + // const XclsTag & xcls = interaction->ExclTag(); + const Target& tgt = init_state.Tgt(); + + // If not QE return 1.0 + // LOG(FIT) << "RikRPA : Event QE = " << proc_info.IsQuasiElastic() << + // std::endl; + if (!tgt.IsNucleus()) return 1.0; + if (!proc_info.IsQuasiElastic()) return 1.0; + + // Extract Beam and Target PDG + GHepParticle* neutrino = ghep->Probe(); + int bpdg = neutrino->Pdg(); + + GHepParticle* target = ghep->Particle(1); + assert(target); + int tpdg = target->Pdg(); + + // Find the enum we need + int calcenum = GetRPACalcEnum(bpdg, tpdg); + if (calcenum == -1) return 1.0; + + // Check we have the RPA Calc setup for this enum + // if not, set it up at that point + // if (!fRPACalculators[calcenum]) SetupRPACalculator(calcenum); + // weightRPA* rpacalc = fRPACalculators[calcenum]; + // if (!rpacalc) { + // THROW("Failed to grab the RPA Calculator : " << calcenum); + // } + + // Extract Q0-Q3 + GHepParticle* fsl = ghep->FinalStatePrimaryLepton(); + const TLorentzVector& k1 = *(neutrino->P4()); + const TLorentzVector& k2 = *(fsl->P4()); + double q0 = fabs((k1 - k2).E()); + double q3 = fabs((k1 - k2).Vect().Mag()); + + // Now use q0-q3 and RPA Calculator to fill fWeights + // LOG(FIT) << "Getting Weights = " << q0 << " " << q3 << std::endl; + // rpacalc->getWeight(q0, q3, fEventWeights); + + // Apply Interpolation (for the time being simple linear) + + // Syst Application : kMINERvA_RikRPA_ApplyRPA + if (fApplyDial_RPACorrection) { + w *= fEventWeights[0]; // CV + } + + LOG(FIT) << " fCurDial_RPALowQ2 = " << fCurDial_RPALowQ2 + << " fCurDial_RPAHighQ2 = " << fCurDial_RPAHighQ2 << " Weights " + << fEventWeights[0] << " " << fEventWeights[1] << " " + << fEventWeights[2] << " " << fEventWeights[3] << " " + << fEventWeights[4] << std::endl; + + // Syst Application : kMINERvA_RikRPA_LowQ2 + if (fabs(fCurDial_RPALowQ2) > 0.0) { + double interpw = fEventWeights[0]; + + if (fCurDial_RPALowQ2 > 0.0) { + interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[1]) * + fCurDial_RPALowQ2; // WLow+ } else if + } else if (fCurDial_RPALowQ2 < 0.0) { + interpw = fEventWeights[0] - (fEventWeights[2] - fEventWeights[0]) * + fCurDial_RPALowQ2; // WLow- + } + w *= interpw / fEventWeights[0]; // Div by CV again + } + + // Syst Application : kMINERvA_RikRPA_HighQ2 + if (fabs(fCurDial_RPAHighQ2) > 0.0) { + double interpw = fEventWeights[0]; + if (fCurDial_RPAHighQ2 > 0.0) { + interpw = fEventWeights[0] - (fEventWeights[0] - fEventWeights[3]) * + fCurDial_RPAHighQ2; // WHigh+ } else + if (fCurDial_RPAHighQ2 < 0.0) { + interpw = fEventWeights[0] - (fEventWeights[4] - fEventWeights[0]) * + fCurDial_RPAHighQ2; // WHigh- + } + w *= interpw / fEventWeights[0]; // Div by CV again + } + } + + // LOG(FIT) << "RPA Weight = " << w << std::endl; + return w; +} // namespace reweight + +void RikRPA::SetDialValue(std::string name, double val) { + SetDialValue(Reweight::ConvDial(name, kCUSTOM), val); +} + +void RikRPA::SetDialValue(int rwenum, double val) { + int curenum = rwenum % 1000; + + // Check Handled + if (!IsHandled(curenum)) return; + if (curenum == kMINERvARW_RikRPA_ApplyRPA) + fApplyDial_RPACorrection = (val > 0.5); + if (curenum == kMINERvARW_RikRPA_LowQ2) fCurDial_RPALowQ2 = val; + if (curenum == kMINERvARW_RikRPA_HighQ2) fCurDial_RPAHighQ2 = val; + + // Assign flag to say stuff has changed + fTweaked = (fApplyDial_RPACorrection || + fabs(fCurDial_RPAHighQ2 - fDefDial_RPAHighQ2) > 0.0 || + fabs(fCurDial_RPALowQ2 - fDefDial_RPALowQ2) > 0.0); +} + +bool RikRPA::IsHandled(int rwenum) { + int curenum = rwenum % 1000; + switch (curenum) { + case kMINERvARW_RikRPA_ApplyRPA: + return true; + case kMINERvARW_RikRPA_LowQ2: + return true; + case kMINERvARW_RikRPA_HighQ2: + return true; + default: + return false; + } +} + +void RikRPA::SetupRPACalculator(int calcenum) { + std::string rwdir = FitPar::GetDataBase() + "/reweight/MINERvA/RikRPA/"; + std::string fidir = ""; + switch (calcenum) { + case kNuMuC12: + fidir = "outNievesRPAratio-nu12C-20GeV-20170202.root"; + break; + + case kNuMuO16: + fidir = "outNievesRPAratio-nu16O-20GeV-20170202.root"; + break; + + case kNuMuAr40: + fidir = "outNievesRPAratio-nu40Ar-20GeV-20170202.root"; + break; + + case kNuMuCa40: + fidir = "outNievesRPAratio-nu40Ca-20GeV-20170202.root"; + break; + + case kNuMuFe56: + fidir = "outNievesRPAratio-nu56Fe-20GeV-20170202.root"; + break; + + case kNuMuBarC12: + fidir = "outNievesRPAratio-anu12C-20GeV-20170202.root"; + break; + + case kNuMuBarO16: + fidir = "outNievesRPAratio-anu16O-20GeV-20170202.root"; + break; + + case kNuMuBarAr40: + fidir = "outNievesRPAratio-anu40Ar-20GeV-20170202.root"; + break; + + case kNuMuBarCa40: + fidir = "outNievesRPAratio-anu40Ca-20GeV-20170202.root"; + break; + + case kNuMuBarFe56: + fidir = "outNievesRPAratio-anu56Fe-20GeV-20170202.root"; + break; + } + + LOG(FIT) << "Loading RPA CALC : " << fidir << std::endl; + TDirectory* olddir = gDirectory; + + // fRPACalculators[calcenum] = new weightRPA(rwdir + "/" + fidir); + olddir->cd(); + return; +} + +int RikRPA::GetRPACalcEnum(int bpdg, int tpdg) { + if (bpdg == 14 && tpdg == 1000060120) + return kNuMuC12; + else if (bpdg == 14 && tpdg == 1000080160) + return kNuMuO16; + else if (bpdg == 14 && tpdg == 1000180400) + return kNuMuAr40; + else if (bpdg == 14 && tpdg == 1000200400) + return kNuMuCa40; + else if (bpdg == 14 && tpdg == 1000280560) + return kNuMuFe56; + else if (bpdg == -14 && tpdg == 1000060120) + return kNuMuBarC12; + else if (bpdg == -14 && tpdg == 1000080160) + return kNuMuBarO16; + else if (bpdg == -14 && tpdg == 1000180400) + return kNuMuBarAr40; + else if (bpdg == -14 && tpdg == 1000200400) + return kNuMuBarCa40; + else if (bpdg == -14 && tpdg == 1000280560) + return kNuMuBarFe56; + else { + ERROR(WRN, "Unknown beam and target combination for RPA Calcs! " + << bpdg << " " << tpdg); + } + + return -1; +} + +} // namespace reweight +} // namespace nuisance +#endif +#endif diff --git a/src/Reweight/MINERvAWeightCalcs.h b/src/Reweight/MINERvAWeightCalcs.h new file mode 100644 index 0000000..8bcba23 --- /dev/null +++ b/src/Reweight/MINERvAWeightCalcs.h @@ -0,0 +1,117 @@ +#ifndef MINERVA_WEIGHT_CALCS +#define MINERVA_WEIGHT_CALCS + +#include + +#ifdef __MINERVA_RW_ENABLED__ +#ifdef __GENIE_ENABLED__ +#include "Conventions/Units.h" +#include "EVGCore/EventRecord.h" +#include "FitEvent.h" +#include "FitParameters.h" +#include "GHEP/GHepParticle.h" +#include "GHEP/GHepRecord.h" +#include "GHEP/GHepUtils.h" +#include "GeneralUtils.h" +#include "NUISANCESyst.h" +#include "NUISANCEWeightCalcs.h" +#include "Ntuple/NtpMCEventRecord.h" +#include "PDG/PDGUtils.h" +#include "WeightUtils.h" +// #include "weightRPA.h" +using namespace genie; +class BaseFitEvt; + + +namespace nuisance { +namespace reweight { + +// MEC Dials +class MINERvAReWeight_MEC : public NUISANCEWeightCalc { + public: + MINERvAReWeight_MEC(); + virtual ~MINERvAReWeight_MEC(); + + double CalcWeight(BaseFitEvt* evt); + void SetDialValue(std::string name, double val); + void SetDialValue(int rwenum, double val); + bool IsHandled(int rwenum); + + double fTwk_NormCCMEC; + double fCur_NormCCMEC; + double fDef_NormCCMEC; + bool fTweaked; + +}; + +// RES Dials +class MINERvAReWeight_RES : public NUISANCEWeightCalc { + public: + MINERvAReWeight_RES(); + virtual ~MINERvAReWeight_RES(); + + double CalcWeight(BaseFitEvt* evt); + void SetDialValue(std::string name, double val); + void SetDialValue(int rwenum, double val); + bool IsHandled(int rwenum); + + double fTwk_NormCCRES; + double fCur_NormCCRES; + double fDef_NormCCRES; + bool fTweaked; + + }; + +/// RPA Weight Calculator that applies RPA systematics +/// to GENIE events. GENIE EVENTS ONLY! +class RikRPA : public NUISANCEWeightCalc { +public: + RikRPA(); + ~RikRPA(); + + double CalcWeight(BaseFitEvt* evt); + void SetDialValue(std::string name, double val); + void SetDialValue(int rwenum, double val); + bool IsHandled(int rwenum); + + void SetupRPACalculator(int calcenum); + int GetRPACalcEnum(int bpdg, int tpdg); + + bool fApplyDial_RPACorrection; + + double fTwkDial_RPALowQ2; + double fDefDial_RPALowQ2; + double fCurDial_RPALowQ2; + double fErrDial_RPALowQ2; + + double fTwkDial_RPAHighQ2; + double fDefDial_RPAHighQ2; + double fCurDial_RPAHighQ2; + double fErrDial_RPAHighQ2; + + double* fEventWeights; + bool fTweaked; + + const static int kMaxCalculators = 10; + enum rpacalcenums { + kNuMuC12, + kNuMuO16, + kNuMuAr40, + kNuMuCa40, + kNuMuFe56, + kNuMuBarC12, + kNuMuBarO16, + kNuMuBarAr40, + kNuMuBarCa40, + kNuMuBarFe56 + }; + // weightRPA* fRPACalculators[kMaxCalculators]; +}; + +}; // namespace reweight +}; // namespace nuisance + +#endif // __GENIE_ENABLED__ +#endif //__MINERVA_RW_ENABLED__ + +#endif diff --git a/src/Reweight/NUISANCESyst.cxx b/src/Reweight/NUISANCESyst.cxx index 44366e1..40aef5e 100644 --- a/src/Reweight/NUISANCESyst.cxx +++ b/src/Reweight/NUISANCESyst.cxx @@ -1,52 +1,60 @@ #include "NUISANCESyst.h" int Reweight::ConvertNUISANCEDial(std::string type) { for (int i = kUnkownNUISANCEDial + 1; i < kNUISANCEDial_LAST; i++) { if (!type.compare(ConvNUISANCEDial(i).c_str())) { return i; } } return kUnkownNUISANCEDial; }; std::string Reweight::ConvNUISANCEDial(int type) { switch (type) { case kGaussianCorr_CCQE_norm: { return "GaussianCorr_CCQE_norm"; } case kGaussianCorr_CCQE_tilt: { return "GaussianCorr_CCQE_tilt"; } case kGaussianCorr_CCQE_Pq0: { return "GaussianCorr_CCQE_Pq0"; } case kGaussianCorr_CCQE_Wq0: { return "GaussianCorr_CCQE_Wq0"; } case kGaussianCorr_CCQE_Pq3: { return "GaussianCorr_CCQE_Pq3"; } case kGaussianCorr_CCQE_Wq3: { return "GaussianCorr_CCQE_Wq3"; } case kGaussianCorr_2p2h_norm: { return "GaussianCorr_2p2h_norm"; } case kGaussianCorr_2p2h_tilt: { return "GaussianCorr_2p2h_tilt"; } case kGaussianCorr_2p2h_Pq0: { return "GaussianCorr_2p2h_Pq0"; } case kGaussianCorr_2p2h_Wq0: { return "GaussianCorr_2p2h_Wq0"; } case kGaussianCorr_2p2h_Pq3: { return "GaussianCorr_2p2h_Pq3"; } case kGaussianCorr_2p2h_Wq3: { return "GaussianCorr_2p2h_Wq3"; } case kGaussianCorr_2p2h_PPandNN_norm: { return "GaussianCorr_2p2h_PPandNN_norm"; } case kGaussianCorr_2p2h_PPandNN_tilt: { return "GaussianCorr_2p2h_PPandNN_tilt"; } case kGaussianCorr_2p2h_PPandNN_Pq0: { return "GaussianCorr_2p2h_PPandNN_Pq0"; } case kGaussianCorr_2p2h_PPandNN_Wq0: { return "GaussianCorr_2p2h_PPandNN_Wq0"; } case kGaussianCorr_2p2h_PPandNN_Pq3: { return "GaussianCorr_2p2h_PPandNN_Pq3"; } case kGaussianCorr_2p2h_PPandNN_Wq3: { return "GaussianCorr_2p2h_PPandNN_Wq3"; } case kGaussianCorr_2p2h_NP_norm: { return "GaussianCorr_2p2h_NP_norm"; } case kGaussianCorr_2p2h_NP_tilt: { return "GaussianCorr_2p2h_NP_tilt"; } case kGaussianCorr_2p2h_NP_Pq0: { return "GaussianCorr_2p2h_NP_Pq0"; } case kGaussianCorr_2p2h_NP_Wq0: { return "GaussianCorr_2p2h_NP_Wq0"; } case kGaussianCorr_2p2h_NP_Pq3: { return "GaussianCorr_2p2h_NP_Pq3"; } case kGaussianCorr_2p2h_NP_Wq3: { return "GaussianCorr_2p2h_NP_Wq3"; } case kGaussianCorr_CC1pi_norm: { return "GaussianCorr_CC1pi_norm"; } case kGaussianCorr_CC1pi_tilt: { return "GaussianCorr_CC1pi_tilt"; } case kGaussianCorr_CC1pi_Pq0: { return "GaussianCorr_CC1pi_Pq0"; } case kGaussianCorr_CC1pi_Wq0: { return "GaussianCorr_CC1pi_Wq0"; } case kGaussianCorr_CC1pi_Pq3: { return "GaussianCorr_CC1pi_Pq3"; } case kGaussianCorr_CC1pi_Wq3: { return "GaussianCorr_CC1pi_Wq3"; } case kGaussianCorr_AllowSuppression: { return "GaussianCorr_AllowSuppression"; } + case kModeNorm_NormRES: { return "NormRES"; } + case kMINERvARW_NormCCMEC: { return "MINERvARW_NormCCMEC"; } + case kMINERvARW_NormCCRES: { return "MINERvARW_NormCCRES"; } + + case kMINERvARW_RikRPA_ApplyRPA: { return "MINERvARW_RikRPA_ApplyRPA"; } + case kMINERvARW_RikRPA_LowQ2: { return "MINERvARW_RikRPA_LowQ2"; } + case kMINERvARW_RikRPA_HighQ2: { return "MINERvARW_RikRPA_HighQ2"; } + default: return "unknown_nuisance_dial"; } }; diff --git a/src/Reweight/NUISANCESyst.h b/src/Reweight/NUISANCESyst.h index 296c956..a3cfb74 100644 --- a/src/Reweight/NUISANCESyst.h +++ b/src/Reweight/NUISANCESyst.h @@ -1,50 +1,57 @@ #ifndef NUISANCESyst_H #define NUISANCESyst_H #include "GeneralUtils.h" namespace Reweight { enum NUISANCESyst { kUnkownNUISANCEDial = 0, kGaussianCorr_CCQE_norm, kGaussianCorr_CCQE_tilt, kGaussianCorr_CCQE_Pq0, kGaussianCorr_CCQE_Wq0, kGaussianCorr_CCQE_Pq3, kGaussianCorr_CCQE_Wq3, kGaussianCorr_2p2h_norm, kGaussianCorr_2p2h_tilt, kGaussianCorr_2p2h_Pq0, kGaussianCorr_2p2h_Wq0, kGaussianCorr_2p2h_Pq3, kGaussianCorr_2p2h_Wq3, kGaussianCorr_2p2h_PPandNN_norm, kGaussianCorr_2p2h_PPandNN_tilt, kGaussianCorr_2p2h_PPandNN_Pq0, kGaussianCorr_2p2h_PPandNN_Wq0, kGaussianCorr_2p2h_PPandNN_Pq3, kGaussianCorr_2p2h_PPandNN_Wq3, kGaussianCorr_2p2h_NP_norm, kGaussianCorr_2p2h_NP_tilt, kGaussianCorr_2p2h_NP_Pq0, kGaussianCorr_2p2h_NP_Wq0, kGaussianCorr_2p2h_NP_Pq3, kGaussianCorr_2p2h_NP_Wq3, kGaussianCorr_CC1pi_norm, kGaussianCorr_CC1pi_tilt, kGaussianCorr_CC1pi_Pq0, kGaussianCorr_CC1pi_Wq0, kGaussianCorr_CC1pi_Pq3, kGaussianCorr_CC1pi_Wq3, kGaussianCorr_AllowSuppression, kModeNorm_NormRES, + kMINERvARW_NormCCMEC, + kMINERvARW_NormCCRES, + + kMINERvARW_RikRPA_ApplyRPA, + kMINERvARW_RikRPA_LowQ2, + kMINERvARW_RikRPA_HighQ2, + kNUISANCEDial_LAST }; int ConvertNUISANCEDial(std::string type); std::string ConvNUISANCEDial(int type); }; #endif diff --git a/src/Reweight/NUISANCEWeightEngine.cxx b/src/Reweight/NUISANCEWeightEngine.cxx index e7344ca..116385c 100644 --- a/src/Reweight/NUISANCEWeightEngine.cxx +++ b/src/Reweight/NUISANCEWeightEngine.cxx @@ -1,129 +1,110 @@ #include "NUISANCEWeightEngine.h" #include "NUISANCEWeightCalcs.h" -NUISANCEWeightEngine::NUISANCEWeightEngine(std::string name) { +#ifdef __MINERVA_RW_ENABLED__ +#ifdef __GENIE_ENABLED__ +#include "MINERvAWeightCalcs.h" +#endif +#endif + +NUISANCEWeightEngine::NUISANCEWeightEngine(std::string name) { // Setup the NUISANCE Reweight engine fCalcName = name; LOG(FIT) << "Setting up NUISANCE Custom RW : " << fCalcName << std::endl; // Load in all Weight Calculations - fWeightCalculators.push_back( new GaussianModeCorr() ); - fWeightCalculators.push_back( new ModeNormCalc() ); + fWeightCalculators.push_back(new GaussianModeCorr()); + fWeightCalculators.push_back(new ModeNormCalc()); + +#ifdef __MINERVA_RW_ENABLED__ +#ifdef __GENIE_ENABLED__ + fWeightCalculators.push_back(new nuisance::reweight::MINERvAReWeight_MEC()); + fWeightCalculators.push_back(new nuisance::reweight::MINERvAReWeight_RES()); + fWeightCalculators.push_back(new nuisance::reweight::RikRPA()); +#endif +#endif + // Set Abs Twk Config fIsAbsTwk = true; - }; - void NUISANCEWeightEngine::IncludeDial(std::string name, double startval) { - // Get First enum int nuisenum = Reweight::ConvDial(name, kCUSTOM); // Setup Maps - fEnumIndex[nuisenum];// = std::vector(0); - fNameIndex[name]; // = std::vector(0); + fEnumIndex[nuisenum]; // = std::vector(0); + fNameIndex[name]; // = std::vector(0); // Split by commas std::vector allnames = GeneralUtils::ParseToStr(name, ","); for (uint i = 0; i < allnames.size(); i++) { std::string singlename = allnames[i]; // Get RW - int singleenum = Reweight::ConvDial(singlename, kCUSTOM); + int singleenum = Reweight::ConvDial(singlename, kCUSTOM); // Fill Maps int index = fValues.size(); fValues.push_back(0.0); fNUISANCEEnums.push_back(singleenum); // Initialize dial std::cout << "Registering " << singlename << " from " << name << std::endl; // Setup index fEnumIndex[nuisenum].push_back(index); fNameIndex[name].push_back(index); } // Set Value if given if (startval != -999.9) { SetDialValue(nuisenum, startval); } }; - void NUISANCEWeightEngine::SetDialValue(int nuisenum, double val) { std::vector indices = fEnumIndex[nuisenum]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; } } void NUISANCEWeightEngine::SetDialValue(std::string name, double val) { std::vector indices = fNameIndex[name]; for (uint i = 0; i < indices.size(); i++) { fValues[indices[i]] = val; } } void NUISANCEWeightEngine::Reconfigure(bool silent) { - for (size_t i = 0; i < fNUISANCEEnums.size(); i++) { - - for (std::vector::iterator calciter = fWeightCalculators.begin(); + for (std::vector::iterator calciter = + fWeightCalculators.begin(); calciter != fWeightCalculators.end(); calciter++) { - - NUISANCEWeightCalc* nuiscalc = static_cast(*calciter); + NUISANCEWeightCalc* nuiscalc = + static_cast(*calciter); if (nuiscalc->IsHandled(fNUISANCEEnums[i])) { nuiscalc->SetDialValue(fNUISANCEEnums[i], fValues[i]); } } } } - - -double NUISANCEWeightEngine::CalcWeight(BaseFitEvt * evt) { +double NUISANCEWeightEngine::CalcWeight(BaseFitEvt* evt) { double rw_weight = 1.0; // Cast as usable class - for (std::vector::iterator iter = fWeightCalculators.begin(); + for (std::vector::iterator iter = + fWeightCalculators.begin(); iter != fWeightCalculators.end(); iter++) { NUISANCEWeightCalc* nuiscalc = static_cast(*iter); rw_weight *= nuiscalc->CalcWeight(evt); } // Return rw_weight return rw_weight; } - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/src/Routines/SystematicRoutines.cxx b/src/Routines/SystematicRoutines.cxx index 6654f71..d33b35e 100755 --- a/src/Routines/SystematicRoutines.cxx +++ b/src/Routines/SystematicRoutines.cxx @@ -1,1515 +1,1517 @@ // 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 "SystematicRoutines.h" void SystematicRoutines::Init(){ fInputFile = ""; fInputRootFile = NULL; fOutputFile = ""; fOutputRootFile = NULL; fCovar = fCovarFree = NULL; fCorrel = fCorrelFree = NULL; fDecomp = fDecompFree = NULL; fStrategy = "ErrorBands"; fRoutines.clear(); fRoutines.push_back("ErrorBands"); fCardFile = ""; fFakeDataInput = ""; fSampleFCN = NULL; fAllowedRoutines = ("ErrorBands,PlotLimits"); }; SystematicRoutines::~SystematicRoutines(){ }; SystematicRoutines::SystematicRoutines(int argc, char* argv[]){ // Initialise Defaults Init(); nuisconfig configuration = Config::Get(); // Default containers std::string cardfile = ""; std::string maxevents = "-1"; int errorcount = 0; int verbocount = 0; std::vector xmlcmds; std::vector configargs; fNThrows = 250; fStartThrows = 0; fThrowString = ""; // 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, "-s", fStartThrows, false, false); ParserUtils::ParseArgument(args, "-t", fNThrows, false, false); ParserUtils::ParseArgument(args, "-p", fThrowString, false, false); ParserUtils::ParseArgument(args, "-i", xmlcmds); ParserUtils::ParseArgument(args, "-q", configargs); ParserUtils::ParseCounter(args, "e", errorcount); ParserUtils::ParseCounter(args, "v", verbocount); ParserUtils::CheckBadArguments(args); // Add extra defaults if none given if (fCardFile.empty() and xmlcmds.empty()) { ERR(FTL) << "No input supplied!" << std::endl; throw; } if (fOutputFile.empty() and !fCardFile.empty()) { fOutputFile = fCardFile + ".root"; ERR(WRN) << "No output supplied so saving it to: " << fOutputFile << std::endl; } else if (fOutputFile.empty()) { ERR(FTL) << "No output file or cardfile supplied!" << std::endl; throw; } // Configuration Setup ============================= // Check no comp key is available if (Config::Get().GetNodes("nuiscomp").empty()) { fCompKey = Config::Get().CreateNode("nuiscomp"); } else { fCompKey = Config::Get().GetNodes("nuiscomp")[0]; } if (!fCardFile.empty()) fCompKey.AddS("cardfile", fCardFile); if (!fOutputFile.empty()) fCompKey.AddS("outputfile", fOutputFile); if (!fStrategy.empty()) fCompKey.AddS("strategy", fStrategy); // Load XML Cardfile configuration.LoadConfig( fCompKey.GetS("cardfile"), ""); // Add CMD XML Structs for (size_t i = 0; i < xmlcmds.size(); i++) { configuration.AddXMLLine(xmlcmds[i]); } // Add Config Args for (size_t i = 0; i < configargs.size(); i++) { configuration.OverrideConfig(configargs[i]); } if (maxevents.compare("-1")){ configuration.OverrideConfig("MAXEVENTS=" + maxevents); } // Finish configuration XML configuration.FinaliseConfig(fCompKey.GetS("outputfile") + ".xml"); // Add Error Verbo Lines verbocount += Config::Get().GetParI("VERBOSITY"); errorcount += Config::Get().GetParI("ERROR"); std::cout << "[ NUISANCE ]: Setting VERBOSITY=" << verbocount << std::endl; std::cout << "[ NUISANCE ]: Setting ERROR=" << errorcount << std::endl; SETVERBOSITY(verbocount); // Proper Setup if (fStrategy.find("ErrorBands") != std::string::npos || fStrategy.find("MergeErrors") != std::string::npos){ fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); } // fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); SetupSystematicsFromXML(); SetupCovariance(); SetupRWEngine(); SetupFCN(); GetCovarFromFCN(); // Run(); return; }; void SystematicRoutines::SetupSystematicsFromXML(){ LOG(FIT) << "Setting up nuismin" << std::endl; // Setup Parameters ------------------------------------------ std::vector 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 if state not 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) { parnom = FitBase::RWAbsToSigma( partype, parname, parnom ); parlow = FitBase::RWAbsToSigma( partype, parname, parlow ); parhigh = FitBase::RWAbsToSigma( partype, parname, parhigh ); parstep = FitBase::RWAbsToSigma( partype, parname, parstep ); } else if (parstate.find("FRAC") != std::string::npos) { parnom = FitBase::RWFracToSigma( partype, parname, parnom ); parlow = FitBase::RWFracToSigma( partype, parname, parlow ); parhigh = FitBase::RWFracToSigma( partype, parname, parhigh ); parstep = FitBase::RWFracToSigma( partype, parname, parstep ); } // Push into vectors fParams.push_back(parname); fTypeVals[parname] = FitBase::ConvDialType(partype);; fStartVals[parname] = parnom; fCurVals[parname] = parnom; fErrorVals[parname] = 0.0; fStateVals[parname] = parstate; bool fixstate = parstate.find("FIX") != std::string::npos; fFixVals[parname] = fixstate; fStartFixVals[parname] = fFixVals[parname]; fMinVals[parname] = parlow; fMaxVals[parname] = parhigh; fStepVals[parname] = parstep; } // Setup Samples ---------------------------------------------- std::vector 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; 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; } } void SystematicRoutines::ReadCard(std::string cardfile){ // Read cardlines into vector std::vector cardlines = GeneralUtils::ParseFileToStr(cardfile,"\n"); FitPar::Config().cardLines = cardlines; // Read Samples first (norm params can be overridden) int linecount = 0; for (std::vector::iterator iter = cardlines.begin(); iter != cardlines.end(); iter++){ std::string line = (*iter); linecount++; // Skip Empties if (line.empty()) continue; if (line.c_str()[0] == '#') continue; // Read Valid Samples int samstatus = ReadSamples(line); // Show line if bad to help user if (samstatus == kErrorStatus) { ERR(FTL) << "Bad Input in cardfile " << fCardFile << " at line " << linecount << "!" << std::endl; LOG(FIT) << line << std::endl; throw; } } // Read Parameters second linecount = 0; for (std::vector::iterator iter = cardlines.begin(); iter != cardlines.end(); iter++){ std::string line = (*iter); linecount++; // Skip Empties if (line.empty()) continue; if (line.c_str()[0] == '#') continue; // Try Parameter Reads int parstatus = ReadParameters(line); int fakstatus = ReadFakeDataPars(line); // Show line if bad to help user if (parstatus == kErrorStatus || fakstatus == kErrorStatus ){ ERR(FTL) << "Bad Parameter Input in cardfile " << fCardFile << " at line " << linecount << "!" << std::endl; LOG(FIT) << line << std::endl; throw; } } return; }; int SystematicRoutines::ReadParameters(std::string parstring){ std::string inputspec = "RW Dial Inputs Syntax \n" "free input w/ limits: TYPE NAME START MIN MAX STEP [STATE] \n" "fix input: TYPE NAME VALUE [STATE] \n" "free input w/o limits: TYPE NAME START FREE,[STATE] \n" "Allowed Types: \n" "neut_parameter,niwg_parameter,t2k_parameter," "nuwro_parameter,gibuu_parameter"; // Check sample input if (parstring.find("parameter") == std::string::npos) return kGoodStatus; // Parse inputs std::vector strvct = GeneralUtils::ParseToStr(parstring, " "); // Skip if comment or parameter somewhere later in line if (strvct[0].c_str()[0] == '#' || strvct[0].find("parameter") == std::string::npos){ return kGoodStatus; } // Check length if (strvct.size() < 3){ ERR(FTL) << "Input rw dials need to provide at least 3 inputs." << std::endl; std::cout << inputspec << std::endl; return kErrorStatus; } // Setup default inputs std::string partype = strvct[0]; std::string parname = strvct[1]; double parval = GeneralUtils::StrToDbl(strvct[2]); double minval = parval - 1.0; double maxval = parval + 1.0; double stepval = 1.0; std::string state = "FIX"; //[DEFAULT] // Check Type if (FitBase::ConvDialType(partype) == kUNKNOWN){ ERR(FTL) << "Unknown parameter type! " << partype << std::endl; std::cout << inputspec << std::endl; return kErrorStatus; } // Check Parameter Name if (FitBase::GetDialEnum(partype, parname) == -1){ ERR(FTL) << "Bad RW parameter name! " << partype << " " << parname << std::endl; std::cout << inputspec << std::endl; return kErrorStatus; } // Option Extra (No Limits) if (strvct.size() == 4){ state = strvct[3]; } // Check for weirder inputs if (strvct.size() > 4 && strvct.size() < 6){ ERR(FTL) << "Provided incomplete limits for " << parname << std::endl; std::cout << inputspec << std::endl; return kErrorStatus; } // Option Extra (With limits and steps) if (strvct.size() >= 6){ minval = GeneralUtils::StrToDbl(strvct[3]); maxval = GeneralUtils::StrToDbl(strvct[4]); stepval = GeneralUtils::StrToDbl(strvct[5]); } // Option Extra (dial state after limits) if (strvct.size() == 7){ state = strvct[6]; } // Run Parameter Conversion if needed if (state.find("ABS") != std::string::npos){ parval = FitBase::RWAbsToSigma( partype, parname, parval ); minval = FitBase::RWAbsToSigma( partype, parname, minval ); maxval = FitBase::RWAbsToSigma( partype, parname, maxval ); stepval = FitBase::RWAbsToSigma( partype, parname, stepval ); } else if (state.find("FRAC") != std::string::npos){ parval = FitBase::RWFracToSigma( partype, parname, parval ); minval = FitBase::RWFracToSigma( partype, parname, minval ); maxval = FitBase::RWFracToSigma( partype, parname, maxval ); stepval = FitBase::RWFracToSigma( partype, parname, stepval ); } // Check no repeat params if (std::find(fParams.begin(), fParams.end(), parname) != fParams.end()){ ERR(FTL) << "Duplicate parameter names given for " << parname << std::endl; throw; } // Setup Containers fParams.push_back(parname); fTypeVals[parname] = FitBase::ConvDialType(partype); fStartVals[parname] = parval; fCurVals[parname] = fStartVals[parname]; fErrorVals[parname] = 0.0; fStateVals[parname] = state; bool fixstate = state.find("FIX") != std::string::npos; fFixVals[parname] = fixstate; fStartFixVals[parname] = fFixVals[parname]; fMinVals[parname] = minval; fMaxVals[parname] = maxval; fStepVals[parname] = stepval; // Print the parameter LOG(MIN) << "Read Parameter " << parname << " " << parval << " " << minval << " " << maxval << " " << stepval << " " << state << std::endl; // Tell reader its all good return kGoodStatus; } //******************************************* int SystematicRoutines::ReadFakeDataPars(std::string parstring){ //****************************************** std::string inputspec = "Fake Data Dial Inputs Syntax \n" "fake value: fake_parameter NAME VALUE \n" "Name should match dialnames given in actual dial specification."; // Check sample input if (parstring.find("fake_parameter") == std::string::npos) return kGoodStatus; // Parse inputs std::vector strvct = GeneralUtils::ParseToStr(parstring, " "); // Skip if comment or parameter somewhere later in line if (strvct[0].c_str()[0] == '#' || strvct[0] == "fake_parameter"){ return kGoodStatus; } // Check length if (strvct.size() < 3){ ERR(FTL) << "Fake dials need to provide at least 3 inputs." << std::endl; std::cout << inputspec << std::endl; return kErrorStatus; } // Read Inputs std::string parname = strvct[1]; double parval = GeneralUtils::StrToDbl(strvct[2]); // Setup Container fFakeVals[parname] = parval; // Print the fake parameter LOG(MIN) << "Read Fake Parameter " << parname << " " << parval << std::endl; // Tell reader its all good return kGoodStatus; } //****************************************** int SystematicRoutines::ReadSamples(std::string samstring){ //****************************************** const static std::string inputspec = "\tsample :inputfile.root [OPTS] " "[norm]\nsample_name: Name " "of sample to include. e.g. MiniBooNE_CCQE_XSec_1DQ2_nu\ninput_type: The " "input event format. e.g. NEUT, GENIE, EVSPLN, ...\nOPTS: Additional, " "optional sample options.\nnorm: Additional, optional sample " "normalisation factor."; // Check sample input if (samstring.find("sample") == std::string::npos) return kGoodStatus; // Parse inputs std::vector strvct = GeneralUtils::ParseToStr(samstring, " "); // Skip if comment or parameter somewhere later in line if (strvct[0].c_str()[0] == '#' || strvct[0] != "sample"){ return kGoodStatus; } // Check length if (strvct.size() < 3){ ERR(FTL) << "Sample need to provide at least 3 inputs." << std::endl; return kErrorStatus; } // Setup default inputs std::string samname = strvct[1]; std::string samfile = strvct[2]; if (samfile == "FIX") { ERR(FTL) << "Input filename was \"FIX\", this line is probably malformed " "in the input card file. Line:\'" << samstring << "\'" << std::endl; ERR(FTL) << "Expect sample lines to look like:\n\t" << inputspec << std::endl; throw; } std::string samtype = "DEFAULT"; double samnorm = 1.0; // Optional Type if (strvct.size() > 3) { samtype = strvct[3]; // samname += "_"+samtype; // Also get rid of the / and replace it with underscore because it might not be supported character // while (samname.find("/") != std::string::npos) { // samname.replace(samname.find("/"), 1, std::string("_")); // } } // Optional Norm if (strvct.size() > 4) samnorm = GeneralUtils::StrToDbl(strvct[4]); // Add Sample Names as Norm Dials std::string normname = samname + "_norm"; // Check no repeat params if (std::find(fParams.begin(), fParams.end(), normname) != fParams.end()){ ERR(FTL) << "Duplicate samples given for " << samname << std::endl; throw; } fParams.push_back(normname); fTypeVals[normname] = kNORM; fStartVals[normname] = samnorm; fCurVals[normname] = fStartVals[normname]; fErrorVals[normname] = 0.0; fMinVals[normname] = 0.1; fMaxVals[normname] = 10.0; fStepVals[normname] = 0.5; bool state = samtype.find("FREE") == std::string::npos; fFixVals[normname] = state; fStartFixVals[normname] = state; // Print read in LOG(MIN) << "Read sample " << samname << " " << samfile << " " << samtype << " " << samnorm << std::endl; // Tell reader its all good return kGoodStatus; } /* Setup Functions */ //************************************* void SystematicRoutines::SetupRWEngine(){ //************************************* for (UInt_t i = 0; i < fParams.size(); i++){ std::string name = fParams[i]; FitBase::GetRW() -> IncludeDial(name, fTypeVals.at(name) ); } UpdateRWEngine(fStartVals); return; } //************************************* void SystematicRoutines::SetupFCN(){ //************************************* LOG(FIT)<<"Making the jointFCN"<Reconfigure(); fSampleFCN->ReconfigureAllEvents(); fSampleFCN->SetFakeData("MC"); UpdateRWEngine(fCurVals); LOG(FIT)<<"Set all data to fake MC predictions."<SetFakeData(fFakeDataInput); } return; } //***************************************** void SystematicRoutines::GetCovarFromFCN(){ //***************************************** LOG(FIT) << "Loading ParamPull objects from FCN to build covar" << std::endl; // Make helperstring std::ostringstream helperstr; // Keep track of what is being thrown std::map dialthrowhandle; // Get Covariance Objects from FCN std::list inputpulls = fSampleFCN->GetPullList(); for (PullListConstIter iter = inputpulls.begin(); iter != inputpulls.end(); iter++){ ParamPull* pull = (*iter); if (pull->GetType().find("THROW")){ fInputThrows.push_back(pull); fInputCovar.push_back(pull->GetFullCovarMatrix()); fInputDials.push_back(pull->GetDataHist()); LOG(FIT) << "Read ParamPull: " << pull->GetName() << " " << pull->GetType() << std::endl; } TH1D dialhist = pull->GetDataHist(); TH1D minhist = pull->GetMinHist(); TH1D maxhist = pull->GetMaxHist(); TH1I typehist = pull->GetDialTypes(); for (int i = 0; i < dialhist.GetNbinsX(); i++){ std::string name = std::string(dialhist.GetXaxis()->GetBinLabel(i+1)); dialthrowhandle[name] = pull->GetName(); if (fCurVals.find(name) == fCurVals.end()){ // Add to Containers fParams.push_back(name); fCurVals[name] = dialhist.GetBinContent(i+1); fStartVals[name] = dialhist.GetBinContent(i+1); fMinVals[name] = minhist.GetBinContent(i+1); fMaxVals[name] = maxhist.GetBinContent(i+1); fStepVals[name] = 1.0; fFixVals[name] = false; fStartFixVals[name] = false; fTypeVals[name] = typehist.GetBinContent(i+1); fStateVals[name] = "FREE" + pull->GetType(); // Maker Helper helperstr << std::string(16, ' ' ) << FitBase::ConvDialType(fTypeVals[name]) << " " << name << " " << fMinVals[name] << " " << fMaxVals[name] << " " << fStepVals[name] << " " << fStateVals[name] << std::endl; } } } // Check if no throws given if (fInputThrows.empty()){ ERR(WRN) << "No covariances given to nuissyst" << std::endl; ERR(WRN) << "Pushing back an uncorrelated gaussian throw error for each free parameter using step size" << std::endl; for (UInt_t i = 0; i < fParams.size(); i++){ std::string syst = fParams[i]; if (fFixVals[syst]) continue; // Make Terms std::string name = syst + "_pull"; std::ostringstream pullterm; pullterm << "DIAL:" << syst << ";" << fStartVals[syst] << ";" << fStepVals[syst]; std::string type = "GAUSTHROW/NEUT"; // Push Back Pulls ParamPull* pull = new ParamPull( name, pullterm.str(), type ); fInputThrows.push_back(pull); fInputCovar.push_back(pull->GetFullCovarMatrix()); fInputDials.push_back(pull->GetDataHist()); // Print Whats added ERR(WRN) << "Added ParamPull : " << name << " " << pullterm.str() << " " << type << std::endl; // Add helper string for future fits helperstr << std::string(16, ' ' ) << "covar " << name << " " << pullterm.str() << " " << type << std::endl; // Keep Track of Throws dialthrowhandle[syst] = pull->GetName(); } } // Print Helper String if (!helperstr.str().empty()){ LOG(FIT) << "To remove these statements in future studies, add the lines below to your card:" << std::endl; // Can't use the logger properly because this can be multi-line. Use cout and added spaces to look better! std::cout << helperstr.str(); sleep(2); } // Print Throw State for (UInt_t i = 0; i < fParams.size(); i++){ std::string syst = fParams[i]; if (dialthrowhandle.find(syst) != dialthrowhandle.end()){ LOG(FIT) << "Dial " << i << ". " << setw(40) << syst << " = THROWING with " << dialthrowhandle[syst] << std::endl; } else { LOG(FIT) << "Dial " << i << ". " << setw(40) << syst << " = FIXED" << std::endl; } } // Pause anyway sleep(1); return; } /* Fitting Functions */ //************************************* void SystematicRoutines::UpdateRWEngine(std::map& updateVals){ //************************************* for (UInt_t i = 0; i < fParams.size(); i++){ std::string name = fParams[i]; if (updateVals.find(name) == updateVals.end()) continue; FitBase::GetRW()->SetDialValue(name,updateVals.at(name)); } FitBase::GetRW()->Reconfigure(); return; } //************************************* void SystematicRoutines::PrintState(){ //************************************* LOG(FIT)<<"------------"<GetLikelihood(); LOG(FIT) << std::left << std::setw(46) << "Likelihood for JointFCN: " << like << std::endl; LOG(FIT)<<"------------"<cd(); SaveCurrentState(); } //************************************* void SystematicRoutines::SaveCurrentState(std::string subdir){ //************************************* LOG(FIT)<<"Saving current full FCN predictions" <mkdir(subdir.c_str()); newdir->cd(); } FitBase::GetRW()->Reconfigure(); fSampleFCN->ReconfigureAllEvents(); fSampleFCN->Write(); // Change back to current DIR curdir->cd(); return; } //************************************* void SystematicRoutines::SaveNominal(){ //************************************* if (!fOutputRootFile) fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); fOutputRootFile->cd(); LOG(FIT)<<"Saving Nominal Predictions (be cautious with this)" <Reconfigure(); SaveCurrentState("nominal"); }; //************************************* void SystematicRoutines::SavePrefit(){ //************************************* if (!fOutputRootFile) fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); fOutputRootFile->cd(); LOG(FIT)<<"Saving Prefit Predictions"< 0){ fCovarFree = new TH2D("covariance_free", "covariance_free", NFREE,0,NFREE, NFREE,0,NFREE); } // Set Bin Labels int countall = 0; int countfree = 0; for (UInt_t i = 0; i < fParams.size(); i++){ fCovar->GetXaxis()->SetBinLabel(countall+1,fParams[i].c_str()); fCovar->GetYaxis()->SetBinLabel(countall+1,fParams[i].c_str()); countall++; if (!fFixVals[fParams[i]] and NFREE > 0){ fCovarFree->GetXaxis()->SetBinLabel(countfree+1,fParams[i].c_str()); fCovarFree->GetYaxis()->SetBinLabel(countfree+1,fParams[i].c_str()); countfree++; } } fCorrel = PlotUtils::GetCorrelationPlot(fCovar,"correlation"); fDecomp = PlotUtils::GetDecompPlot(fCovar,"decomposition"); if (NFREE > 0)fCorrelFree = PlotUtils::GetCorrelationPlot(fCovarFree, "correlation_free"); if (NFREE > 0)fDecompFree = PlotUtils::GetDecompPlot(fCovarFree,"decomposition_free"); return; }; //************************************* void SystematicRoutines::ThrowCovariance(bool uniformly){ //************************************* // Set fThrownVals to all values in currentVals for (UInt_t i = 0; i < fParams.size(); i++){ std::string name = fParams.at(i); fThrownVals[name] = fCurVals[name]; } for (PullListConstIter iter = fInputThrows.begin(); iter != fInputThrows.end(); iter++){ ParamPull* pull = *iter; pull->ThrowCovariance(); TH1D dialhist = pull->GetDataHist(); for (int i = 0; i < dialhist.GetNbinsX(); i++){ std::string name = std::string(dialhist.GetXaxis()->GetBinLabel(i+1)); if (fCurVals.find(name) != fCurVals.end()){ fThrownVals[name] = dialhist.GetBinContent(i+1); } } // Reset throw incase pulls are calculated. pull->ResetToy(); } return; }; //************************************* void SystematicRoutines::PlotLimits(){ //************************************* std::cout << "Plotting Limits" << std::endl; if (!fOutputRootFile) fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); TDirectory* limfolder = (TDirectory*) fOutputRootFile->mkdir("Limits"); limfolder->cd(); // Set all parameters at their starting values for (UInt_t i = 0; i < fParams.size(); i++){ fCurVals[fParams[i]] = fStartVals[fParams[i]]; } TDirectory* nomfolder = (TDirectory*) limfolder->mkdir("nominal"); nomfolder->cd(); UpdateRWEngine(fCurVals); fSampleFCN->ReconfigureAllEvents(); fSampleFCN->Write(); limfolder->cd(); std::vector allfolders; // Loop through each parameter for (UInt_t i = 0; i < fParams.size(); i++){ std::string syst = fParams[i]; std::cout << "Starting Param " << syst << std::endl; if (fFixVals[syst]) continue; // Loop Downwards while (fCurVals[syst] > fMinVals[syst]){ fCurVals[syst] = fCurVals[syst] - fStepVals[syst]; // Check Limit if (fCurVals[syst] < fMinVals[syst]) fCurVals[syst] = fMinVals[syst]; // Check folder exists std::string curvalstring = std::string( Form( (syst + "_%f").c_str(), fCurVals[syst] ) ); if (std::find(allfolders.begin(), allfolders.end(), curvalstring) != allfolders.end()) break; // Make new folder for variation TDirectory* minfolder = (TDirectory*) limfolder->mkdir(Form( (syst + "_%f").c_str(), fCurVals[syst] ) ); minfolder->cd(); allfolders.push_back(curvalstring); // Update Iterations double *vals = FitUtils::GetArrayFromMap( fParams, fCurVals ); fSampleFCN->DoEval( vals ); delete vals; // Save to folder fSampleFCN->Write(); } // Reset before next loop fCurVals[syst] = fStartVals[syst]; // Loop Upwards now while (fCurVals[syst] < fMaxVals[syst]){ fCurVals[syst] = fCurVals[syst] + fStepVals[syst]; // Check Limit if (fCurVals[syst] > fMaxVals[syst]) fCurVals[syst] = fMaxVals[syst]; // Check folder exists std::string curvalstring = std::string( Form( (syst + "_%f").c_str(), fCurVals[syst] ) ); if (std::find(allfolders.begin(), allfolders.end(), curvalstring) != allfolders.end()) break; // Make new folder TDirectory* maxfolder = (TDirectory*) limfolder->mkdir(Form( (syst + "_%f").c_str(), fCurVals[syst] ) ); maxfolder->cd(); allfolders.push_back(curvalstring); // Update Iterations double *vals = FitUtils::GetArrayFromMap( fParams, fCurVals ); fSampleFCN->DoEval( vals ); delete vals; // Save to file fSampleFCN->Write(); } // Reset before leaving fCurVals[syst] = fStartVals[syst]; UpdateRWEngine(fCurVals); } return; } //************************************* void SystematicRoutines::Run(){ //************************************* std::cout << "Running routines "<< std::endl; fRoutines = GeneralUtils::ParseToStr(fStrategy,","); for (UInt_t i = 0; i < fRoutines.size(); i++){ std::string routine = fRoutines.at(i); int fitstate = kFitUnfinished; LOG(FIT)<<"Running Routine: "<cd(); int nthrows = fNThrows; int startthrows = fStartThrows; int endthrows = startthrows + nthrows; if (nthrows < 0) nthrows = endthrows; if (startthrows < 0) startthrows = 0; if (endthrows < 0) endthrows = startthrows + nthrows; int seed = (gRandom->Uniform(0.0,1.0)*100000 + 100000000*(startthrows + endthrows) + time(NULL))/35; gRandom->SetSeed(seed); LOG(FIT) << "Using Seed : " << seed << std::endl; LOG(FIT) << "nthrows = " << nthrows << std::endl; LOG(FIT) << "startthrows = " << startthrows << std::endl; LOG(FIT) << "endthrows = " << endthrows << std::endl; UpdateRWEngine(fCurVals); fSampleFCN->ReconfigureAllEvents(); if (startthrows == 0){ LOG(FIT) << "Making nominal " << std::endl; TDirectory* nominal = (TDirectory*) tempfile->mkdir("nominal"); nominal->cd(); fSampleFCN->Write(); } LOG(SAM) << "nthrows = " << nthrows << std::endl; LOG(SAM) << "startthrows = " << startthrows << std::endl; LOG(SAM) << "endthrows = " << endthrows << std::endl; TTree* parameterTree = new TTree("throws","throws"); double chi2; for (UInt_t i = 0; i < fParams.size(); i++) parameterTree->Branch(fParams[i].c_str(), &fThrownVals[fParams[i]], (fParams[i] + "/D").c_str()); parameterTree->Branch("chi2",&chi2,"chi2/D"); fSampleFCN->CreateIterationTree("error_iterations", FitBase::GetRW()); // Would anybody actually want to do uniform throws of any parameter?? bool uniformly = FitPar::Config().GetParB("error_uniform"); // Run Throws and save for (Int_t i = 0; i < endthrows+1; i++){ LOG(FIT) << "Loop " << i << std::endl; ThrowCovariance(uniformly); if (i < startthrows) continue; if (i == 0) continue; LOG(FIT) << "Throw " << i << " ================================" << std::endl; // Generate Random Parameter Throw // ThrowCovariance(uniformly); TDirectory* throwfolder = (TDirectory*)tempfile->mkdir(Form("throw_%i",i)); throwfolder->cd(); // Run Eval double *vals = FitUtils::GetArrayFromMap( fParams, fThrownVals ); chi2 = fSampleFCN->DoEval( vals ); delete vals; // Save the FCN fSampleFCN->Write(); parameterTree->Fill(); } + tempfile->cd(); fSampleFCN->WriteIterationTree(); tempfile->Close(); } void SystematicRoutines::MergeThrows(){ fOutputRootFile = new TFile(fCompKey.GetS("outputfile").c_str(), "RECREATE"); fOutputRootFile->cd(); // Make a container folder TDirectory* errorDIR = (TDirectory*) fOutputRootFile->mkdir("error_bands"); errorDIR->cd(); TDirectory* outnominal = (TDirectory*) fOutputRootFile->mkdir("nominal_throw"); outnominal->cd(); // Split Input Files if (!fThrowString.empty()) fThrowList = GeneralUtils::ParseToStr(fThrowString,","); // Add default if no throwlist given if (fThrowList.size() < 1) fThrowList.push_back( fOutputFile + ".throws.root" ); /// Save location of file containing nominal std::string nominalfile; bool nominalfound; // Loop over files and check they exist. for (uint i = 0; i < fThrowList.size(); i++){ std::string file = fThrowList[i]; bool found = false; // normal std::string newfile = file; TFile* throwfile = new TFile(file.c_str(),"READ"); if (throwfile and !throwfile->IsZombie()){ found = true; } // normal.throws.root if (!found){ newfile = file + ".throws.root"; throwfile = new TFile((file + ".throws.root").c_str(),"READ"); if (throwfile and !throwfile->IsZombie()) { found = true; } } // If its found save to throwlist, else save empty. // Also search for nominal if (found){ fThrowList[i] = newfile; LOG(FIT) << "Throws File :" << newfile << std::endl; // Find input which contains nominal if (throwfile->Get("nominal")){ nominalfound = true; nominalfile = newfile; } throwfile->Close(); } else { fThrowList[i] = ""; } delete throwfile; } // Make sure we have a nominal file if (!nominalfound or nominalfile.empty()){ ERR(FTL) << "No nominal found when mergining! Exiting!" << std::endl; throw; } // Get the nominal throws file TFile* tempfile = new TFile((nominalfile).c_str(),"READ"); tempfile->cd(); TDirectory* nominal = (TDirectory*)tempfile->Get("nominal"); // int nthrows = FitPar::Config().GetParI("error_throws"); bool uniformly = FitPar::Config().GetParB("error_uniform"); // Check percentage of bad files is okay. int badfilecount = 0; for (uint i = 0; i < fThrowList.size(); i++){ if (!fThrowList[i].empty()){ LOG(FIT) << "Loading Throws From File " << i << " : " << fThrowList[i] << std::endl; } else { badfilecount++; } } // Check we have at least one good file if ((uint)badfilecount == fThrowList.size()){ ERR(FTL) << "Found no good throw files for MergeThrows" << std::endl; throw; } else if (badfilecount > fThrowList.size()*0.25){ ERR(WRN) << "Over 25% of your throw files are dodgy. Please check this is okay!" << std::endl; ERR(WRN) << "Will continue for the time being..." << std::endl; sleep(5); } // Now go through the keys in the temporary file and look for TH1D, and TH2D plots TIter next(nominal->GetListOfKeys()); TKey *key; while ((key = (TKey*)next())) { TClass *cl = gROOT->GetClass(key->GetClassName()); if (!cl->InheritsFrom("TH1D") and !cl->InheritsFrom("TH2D")) continue; TH1* baseplot = (TH1D*)key->ReadObj(); std::string plotname = std::string(baseplot->GetName()); LOG(FIT) << "Creating error bands for " << plotname; if (LOG_LEVEL(FIT)){ if (!uniformly) std::cout << " : Using COVARIANCE Throws! " << std::endl; else std::cout << " : Using UNIFORM THROWS!!! " << std::endl; } int nbins = 0; if (cl->InheritsFrom("TH1D")) nbins = ((TH1D*)baseplot)->GetNbinsX(); else nbins = ((TH1D*)baseplot)->GetNbinsX()* ((TH1D*)baseplot)->GetNbinsY(); // Setup TProfile with RMS option TProfile* tprof = new TProfile((plotname + "_prof").c_str(),(plotname + "_prof").c_str(),nbins, 0, nbins, "S"); // Setup The TTREE double* bincontents; bincontents = new double[nbins]; double* binlowest; binlowest = new double[nbins]; double* binhighest; binhighest = new double[nbins]; errorDIR->cd(); TTree* bintree = new TTree((plotname + "_tree").c_str(), (plotname + "_tree").c_str()); for (Int_t i = 0; i < nbins; i++){ bincontents[i] = 0.0; binhighest[i] = 0.0; binlowest[i] = 0.0; bintree->Branch(Form("content_%i",i),&bincontents[i],Form("content_%i/D",i)); } // Make new throw plot TH1* newplot; // Run Throw Merging. for (UInt_t i = 0; i < fThrowList.size(); i++){ TFile* throwfile = new TFile(fThrowList[i].c_str(), "READ"); // Loop over all throws in a folder TIter nextthrow(throwfile->GetListOfKeys()); TKey *throwkey; while ((throwkey = (TKey*)nextthrow())) { // Skip non throw folders if (std::string(throwkey->GetName()).find("throw_") == std::string::npos) continue; // Get Throw DIR TDirectory* throwdir = (TDirectory*)throwkey->ReadObj(); // Get Plot From Throw newplot = (TH1*)throwdir->Get(plotname.c_str()); if (!newplot) continue; // Loop Over Plot for (Int_t j = 0; j < nbins; j++){ tprof->Fill(j+0.5, newplot->GetBinContent(j+1)); bincontents[j] = newplot->GetBinContent(j+1); if (bincontents[j] < binlowest[j] or i == 0) binlowest[j] = bincontents[j]; if (bincontents[j] > binhighest[j] or i == 0) binhighest[j] = bincontents[j]; } errorDIR->cd(); bintree->Fill(); } } errorDIR->cd(); if (uniformly){ LOG(FIT) << "Uniformly Calculating Plot Errors!" << std::endl; } TH1* statplot = (TH1*) baseplot->Clone(); for (Int_t j = 0; j < nbins; j++){ if (!uniformly){ // if ((baseplot->GetBinError(j+1)/baseplot->GetBinContent(j+1)) < 1.0) { // baseplot->SetBinError(j+1,sqrt(pow(tprof->GetBinError(j+1),2) + pow(baseplot->GetBinError(j+1),2))); // } else { - baseplot->SetBinError(j+1,tprof->GetBinError(j+1)); + baseplot->SetBinContent(j+1,tprof->GetBinContent(j+1)); + baseplot->SetBinError(j+1,tprof->GetBinError(j+1)); // } } else { baseplot->SetBinContent(j+1, 0.0);//(binlowest[j] + binhighest[j]) / 2.0); baseplot->SetBinError(j+1, 0.0); //(binhighest[j] - binlowest[j])/2.0); } } errorDIR->cd(); baseplot->Write(); tprof->Write(); bintree->Write(); outnominal->cd(); for (int i = 0; i < nbins; i++){ baseplot->SetBinError(i+1, sqrt(pow(statplot->GetBinError(i+1),2) + pow(baseplot->GetBinError(i+1),2))); } baseplot->Write(); delete statplot; delete baseplot; delete tprof; delete bintree; delete [] bincontents; } return; };