diff --git a/CMakeLists.txt b/CMakeLists.txt index eb7a6d1..2b088f9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,78 +1,78 @@ # Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret ################################################################################ # This file is part of NUISANCE. # # NUISANCE is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # NUISANCE is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with NUISANCE. If not, see . ################################################################################ cmake_minimum_required (VERSION 2.8 FATAL_ERROR) #Use the compilers found in the path find_program(CMAKE_C_COMPILER NAMES $ENV{CC} gcc PATHS ENV PATH NO_DEFAULT_PATH) find_program(CMAKE_CXX_COMPILER NAMES $ENV{CXX} g++ PATHS ENV PATH NO_DEFAULT_PATH) project(NUISANCE) include(ExternalProject) set (NUISANCE_VERSION_MAJOR 3) set (NUISANCE_VERSION_MINOR 0) set (NUISANCE_VERSION_REVISION 0) set (NUISANCE_VERSION_STRING "v${NUISANCE_VERSION_MAJOR}r${NUISANCE_VERSION_MINOR}") if(${NUISANCE_VERSION_REVISION} STRGREATER "0") set (NUISANCE_VERSION_STRING "${NUISANCE_VERSION_STRING}p${NUISANCE_VERSION_REVISION}") endif() #Set this to TRUE to enable build debugging messages set(BUILD_DEBUG_MSGS TRUE) include(${CMAKE_SOURCE_DIR}/cmake/cmessage.cmake) include(${CMAKE_SOURCE_DIR}/cmake/cacheVariables.cmake) include(${CMAKE_SOURCE_DIR}/cmake/parseConfigApp.cmake) include(${CMAKE_SOURCE_DIR}/cmake/StringAndListUtils.cmake) cmessage(STATUS "CMAKE_INSTALL_PREFIX: \"${CMAKE_INSTALL_PREFIX}\"") cmessage(STATUS "CMAKE_BUILD_TYPE: \"${CMAKE_BUILD_TYPE}\"") ################################################################################ # Check Dependencies ################################################################################ ################################## ROOT ###################################### include(${CMAKE_SOURCE_DIR}/cmake/ROOTSetup.cmake) ################################# InputHandler ################################# include(${CMAKE_SOURCE_DIR}/cmake/InputHandlerSetup.cmake) ################################# Pythia6/8 ################################### include(${CMAKE_SOURCE_DIR}/cmake/pythia6Setup.cmake) +#Want this before fhiclcpp which will add the install directory +include_directories(${CMAKE_SOURCE_DIR}/src) + ################################## FHICLCPP #################################### include(${CMAKE_SOURCE_DIR}/cmake/fhiclcppSetup.cmake) -#Want this after fhiclcpp which will add the install directory -include_directories(${CMAKE_SOURCE_DIR}/src) - ################################## COMPILER #################################### include(${CMAKE_SOURCE_DIR}/cmake/c++CompilerSetup.cmake) ################################################################################ add_subdirectory(src) add_subdirectory(scripts) add_subdirectory(config) add_subdirectory(data) diff --git a/data/nuA/Nuclear/CMakeLists.txt b/data/nuA/Nuclear/CMakeLists.txt index b7a435f..ca939a6 100644 --- a/data/nuA/Nuclear/CMakeLists.txt +++ b/data/nuA/Nuclear/CMakeLists.txt @@ -1,3 +1,4 @@ SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/Nuclear) add_subdirectory(T2K) +add_subdirectory(MINERvA) diff --git a/data/nuA/Nuclear/MINERvA/CC0Pi/CMakeLists.txt b/data/nuA/Nuclear/MINERvA/CC0Pi/CMakeLists.txt new file mode 100644 index 0000000..f3c91ea --- /dev/null +++ b/data/nuA/Nuclear/MINERvA/CC0Pi/CMakeLists.txt @@ -0,0 +1,3 @@ +SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/CC0Pi) + +add_subdirectory(NProt_CH_xsec_STV_nu) diff --git a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt b/data/nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/CMakeLists.txt similarity index 59% copy from data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt copy to data/nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/CMakeLists.txt index 18cc6d1..2cb89d2 100644 --- a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt +++ b/data/nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/CMakeLists.txt @@ -1,5 +1,5 @@ -SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/CC0Pi) +SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/NProt_CH_xsec_STV_nu) FILE(GLOB DATA_FILES *.root) install(FILES ${DATA_FILES} DESTINATION ${DATA_INSTALL_DIR}) diff --git a/data/nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/MINERvA_1805.05486.root b/data/nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/MINERvA_1805.05486.root new file mode 100644 index 0000000..7d51761 Binary files /dev/null and b/data/nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/MINERvA_1805.05486.root differ diff --git a/data/nuA/Nuclear/MINERvA/CMakeLists.txt b/data/nuA/Nuclear/MINERvA/CMakeLists.txt new file mode 100644 index 0000000..4a70bbe --- /dev/null +++ b/data/nuA/Nuclear/MINERvA/CMakeLists.txt @@ -0,0 +1,3 @@ +SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/MINERvA) + +add_subdirectory(CC0Pi) diff --git a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt b/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt index 18cc6d1..e447144 100644 --- a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt +++ b/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt @@ -1,5 +1,4 @@ SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/CC0Pi) -FILE(GLOB DATA_FILES *.root) - -install(FILES ${DATA_FILES} DESTINATION ${DATA_INSTALL_DIR}) +add_subdirectory(NProt_CH_xsec_STV_nu) +add_subdirectory(H20_xsec_2Dpcthetamu_numubar) diff --git a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt b/data/nuA/Nuclear/T2K/CC0Pi/H20_xsec_2Dpcthetamu_numubar/CMakeLists.txt similarity index 56% copy from data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt copy to data/nuA/Nuclear/T2K/CC0Pi/H20_xsec_2Dpcthetamu_numubar/CMakeLists.txt index 18cc6d1..5ca9b21 100644 --- a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt +++ b/data/nuA/Nuclear/T2K/CC0Pi/H20_xsec_2Dpcthetamu_numubar/CMakeLists.txt @@ -1,5 +1,5 @@ -SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/CC0Pi) +SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/H20_xsec_2Dpcthetamu_numubar) FILE(GLOB DATA_FILES *.root) install(FILES ${DATA_FILES} DESTINATION ${DATA_INSTALL_DIR}) diff --git a/data/nuA/Nuclear/T2K/CC0Pi/H2O_xsec_2Dpmuthetamu_numubar.root b/data/nuA/Nuclear/T2K/CC0Pi/H20_xsec_2Dpcthetamu_numubar/H2O_xsec_2Dpmuthetamu_numubar.root similarity index 100% rename from data/nuA/Nuclear/T2K/CC0Pi/H2O_xsec_2Dpmuthetamu_numubar.root rename to data/nuA/Nuclear/T2K/CC0Pi/H20_xsec_2Dpcthetamu_numubar/H2O_xsec_2Dpmuthetamu_numubar.root diff --git a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt b/data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/CMakeLists.txt similarity index 59% copy from data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt copy to data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/CMakeLists.txt index 18cc6d1..2cb89d2 100644 --- a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt +++ b/data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/CMakeLists.txt @@ -1,5 +1,5 @@ -SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/CC0Pi) +SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/NProt_CH_xsec_STV_nu) FILE(GLOB DATA_FILES *.root) install(FILES ${DATA_FILES} DESTINATION ${DATA_INSTALL_DIR}) diff --git a/data/nuA/Nuclear/T2K/CC0Pi/datResults.root b/data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/datResults.root similarity index 100% rename from data/nuA/Nuclear/T2K/CC0Pi/datResults.root rename to data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/datResults.root diff --git a/data/nuA/Nuclear/T2K/CC0Pi/dphitResults.root b/data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/dphitResults.root similarity index 100% rename from data/nuA/Nuclear/T2K/CC0Pi/dphitResults.root rename to data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/dphitResults.root diff --git a/data/nuA/Nuclear/T2K/CC0Pi/dptResults.root b/data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/dptResults.root similarity index 100% rename from data/nuA/Nuclear/T2K/CC0Pi/dptResults.root rename to data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/dptResults.root diff --git a/data/nuA/Nuclear/T2K/CMakeLists.txt b/data/nuA/Nuclear/T2K/CMakeLists.txt index 25096b5..da3e52c 100644 --- a/data/nuA/Nuclear/T2K/CMakeLists.txt +++ b/data/nuA/Nuclear/T2K/CMakeLists.txt @@ -1,3 +1,4 @@ SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/T2K) add_subdirectory(CC0Pi) +add_subdirectory(CC1Pi) diff --git a/src/app/CMakeLists.txt b/src/app/CMakeLists.txt index bb1876d..c84c5fe 100644 --- a/src/app/CMakeLists.txt +++ b/src/app/CMakeLists.txt @@ -1,15 +1,15 @@ SET(APPS nuissamples nuiscomp nuisevsum nuisstudy) foreach(a ${APPS}) add_executable(${a} ${a}.cxx) - target_link_libraries(${a} -Wl,--no-as-needed nuis_event nuis_input nuis_utility nuis_config nuis_persistency nuis_plugins nuis_params nuis_variation) + target_link_libraries(${a} -Wl,--no-as-needed nuis_event nuis_input nuis_utility_experimental nuis_utility nuis_config nuis_persistency nuis_plugins nuis_params nuis_variation) target_link_libraries(${a} ${GENERATOR_DEPENDENT_TARGETS} -Wl,--as-needed) target_link_libraries(${a} ${NUISANCE_LINK_DIRS}) target_link_libraries(${a} ${NUISANCE_DEPEND_LIBS}) if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "") set_target_properties(${a} PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS}) endif() install(TARGETS ${a} DESTINATION bin) endforeach() diff --git a/src/app/nuisstudy.cxx b/src/app/nuisstudy.cxx index 4098cd8..5d1d199 100644 --- a/src/app/nuisstudy.cxx +++ b/src/app/nuisstudy.cxx @@ -1,97 +1,112 @@ #include "config/GlobalConfiguration.hxx" #include "input/InputManager.hxx" #include "event/MinimalEvent.hxx" #include "samples/IEventProcessor.hxx" #include "plugins/Instantiate.hxx" #include "exception/exception.hxx" #include "persistency/ROOTOutput.hxx" +#include "utility/StringUtility.hxx" + #include "fhiclcpp/make_ParameterSet.h" #include NEW_NUIS_EXCEPT(invalid_cli_arguments); void SayUsage(char const *argv[]) { std::cout << "[USAGE]: " << argv[0] << "\n" "\t-c : FHiCL file containing study " "configuration. \n" "\t-s : FHiCL key of a single sample " "to run from the -c argument. \n" + "\t-o : Default output stream name " + "override, if unspecified mode will be 'CREATE'. \n" << std::endl; } std::string fhicl_file = ""; std::string named_sample = ""; +std::string default_stream_override = ""; void handleOpts(int argc, char const *argv[]) { int opt = 1; while (opt < argc) { if ((std::string(argv[opt]) == "-?") || (std::string(argv[opt]) == "--help")) { SayUsage(argv); exit(0); } else if (std::string(argv[opt]) == "-c") { fhicl_file = argv[++opt]; } else if (std::string(argv[opt]) == "-s") { named_sample = argv[++opt]; + } else if (std::string(argv[opt]) == "-o") { + default_stream_override = argv[++opt]; } else { std::cout << "[ERROR]: Unknown option: " << argv[opt] << std::endl; SayUsage(argv); exit(1); } opt++; } } int main(int argc, char const *argv[]) { nuis::config::EnsureConfigurationRead("nuis.global.config.fcl"); nuis::config::EnsureConfigurationRead("nuis.datacomparisons.fcl"); handleOpts(argc, argv); if (!fhicl_file.size()) { SayUsage(argv); throw invalid_cli_arguments(); } nuis::config::EnsureConfigurationRead(fhicl_file); + if (default_stream_override.size()) { + std::vector split = + nuis::utility::split(default_stream_override, ":"); + std::string open_mode = (split.size() > 1) ? split[1] : "CREATE"; + + nuis::persistency::NewStream("default", split[0], open_mode); + } + size_t NMax = nuis::config::GetDocument().get( "nmax", std::numeric_limits::max()); std::vector samples; if (named_sample.size()) { samples.push_back( nuis::config::GetDocument().get(named_sample)); } else { samples = nuis::config::GetDocument().get>( "samples"); } for (fhicl::ParameterSet const &samp_config : samples) { std::cout << "[INFO]: Reading sample: " << samp_config.get("name") << std::endl; nuis::plugins::plugin_traits::unique_ptr_t sample = nuis::plugins::Instantiate( samp_config.get("name")); sample->Initialize(samp_config); sample->ProcessSample(NMax); sample->Write(); // Ensures no re-use of samples but cleans up the memory. nuis::input::InputManager::Get().Clear(); } nuis::persistency::CloseOpenTFiles(); } diff --git a/src/event/Particle.hxx b/src/event/Particle.hxx index f71d2e9..67a7863 100644 --- a/src/event/Particle.hxx +++ b/src/event/Particle.hxx @@ -1,89 +1,90 @@ // Copyright 2018 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 . *******************************************************************************/ #pragma once #include "event/types.hxx" #include "TLorentzVector.h" namespace nuis { namespace event { class Particle { public: NEW_NUIS_EXCEPT(invalid_particle); #define STATUS_LIST \ X(kNuclearLeaving, 0) \ X(kPrimaryInitialState, 1) \ X(kPrimaryFinalState, 2) \ X(kIntermediate, 3) \ X(kUnknown, 4) \ X(kBlocked, 5) \ X(kNParticleStatus, 6) #define X(A, B) A = B, enum class Status_t { STATUS_LIST }; #undef X Particle(); Particle(Particle const &); PDG_t pdg; TLorentzVector P4; bool operator!() const { return (pdg == std::numeric_limits::max()); } double E() const { return P4.E(); } double KE() const { return P4.E() - P4.M(); } double P() const { return P4.Vect().Mag(); } TVector3 P3() const { return P4.Vect(); } TVector3 Dir() const { return P3().Unit(); } double M() const { return P4.M(); } double Theta() const { return P4.Vect().Theta(); } + double Theta_deg() const { return P4.Vect().Theta() * (180.0 / M_PI); } double CosTheta() const { return P4.Vect().CosTheta(); } }; } // namespace event } // namespace nuis #define X(A, B) \ case nuis::event::Particle::Status_t::A: { \ return os << #A; \ } inline std::ostream &operator<<(std::ostream &os, nuis::event::Particle::Status_t te) { switch (te) { STATUS_LIST } return os; } #undef X #undef STATUS_LIST inline bool operator==(nuis::event::Particle const &l, nuis::event::Particle const &r) { if (l.pdg != r.pdg) { return false; } for (size_t i = 0; i < 4; ++i) { if (l.P4[i] != r.P4[i]) { return false; } } return true; } diff --git a/src/event/types.hxx b/src/event/types.hxx index 8850d33..cfcc983 100644 --- a/src/event/types.hxx +++ b/src/event/types.hxx @@ -1,129 +1,124 @@ // Copyright 2018 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 EVENT_TYPES_HXX_SEEN #define EVENT_TYPES_HXX_SEEN #include "exception/exception.hxx" namespace nuis { namespace event { // Modelled off the NEUT interaction codes. #define NUIS_INTERACTION_CHANNEL_LIST \ X(kCCQE, 1) \ X(kCC2p2h, 2) \ X(kCCSPP_PPip, 11) \ X(kCCSPP_PPi0, 12) \ X(kCCSPP_NPip, 13) \ X(kCCDiffPP, 15) \ X(kCCCohPi, 16) \ X(kCCResGamma, 17) \ X(kCCTransitionMPi, 21) \ X(kCCResEta0, 22) \ X(kCCResK, 23) \ X(kCCDIS, 26) \ \ X(kNCSPP_NPi0, 31) \ X(kNCSPP_PPi0, 32) \ X(kNCSPP_PPim, 33) \ X(kNCSPP_NPip, 34) \ X(kNCDiffPP, 35) \ X(kNCCohPi, 36) \ X(kNCResNGamma, 38) \ X(kNCResPGamma, 39) \ X(kNCTransitionMPi, 41) \ X(kNCResNEta0, 42) \ X(kNCResPEta0, 43) \ X(kNCResK0, 44) \ X(kNCResKp, 45) \ X(kNCDIS, 46) \ X(kNCELP, 51) \ X(kNCELN, 52) \ X(kNC2p2h, 53) \ \ X(kCCQE_nub, -1) \ X(kCC2p2h_nub, -2) \ X(kCCSPP_NPim_nub, -11) \ X(kCCSPP_NPi0_nub, -12) \ X(kCCSPP_PPim_nub, -13) \ X(kCCDiffPP_nub, -15) \ X(kCCCohPi_nub, -16) \ X(kCCResGamma_nub, -17) \ X(kCCTransitionMPi_nub, -21) \ X(kCCResEta0_nub, -22) \ X(kCCResK_nub, -23) \ X(kCCDIS_nub, -26) \ \ X(kNCSPP_NPi0_nub, -31) \ X(kNCSPP_PPi0_nub, -32) \ X(kNCSPP_PPim_nub, -33) \ X(kNCSPP_NPip_nub, -34) \ X(kNCDiffPP_nub, -35) \ X(kNCCohPi_nub, -36) \ X(kNCResNGamma_nub, -38) \ X(kNCResPGamma_nub, -39) \ X(kNCTransitionMPi_nub, -41) \ X(kNCResNEta0_nub, -42) \ X(kNCResPEta0_nub, -43) \ X(kNCResK0_nub, -44) \ X(kNCResKp_nub, -45) \ X(kNCDIS_nub, -46) \ X(kNCELP_nub, -51) \ X(kNCELN_nub, -52) \ X(kNC2p2h_nub, -53) \ \ X(kUndefined, 0) #define X(A, B) A = B, enum class Channel_t { NUIS_INTERACTION_CHANNEL_LIST }; #undef X #define X(A, B) \ case B: { \ return Channel_t::A; \ } inline Channel_t FromNEUTCode(int nc) { switch (nc) { NUIS_INTERACTION_CHANNEL_LIST default: { return Channel_t::kUndefined; } } } #undef X using PDG_t = long; -inline bool Is2p2h(Channel_t chan) { - return ((chan == Channel_t::kCC2p2h) || (chan == Channel_t::kNC2p2h) || - (chan == Channel_t::kCC2p2h_nub) || (chan == Channel_t::kNC2p2h_nub)); -} - } // namespace event } // namespace nuis #define X(A, B) \ case nuis::event::Channel_t::A: { \ return os << #A; \ } inline std::ostream &operator<<(std::ostream &os, nuis::event::Channel_t te) { switch (te) { NUIS_INTERACTION_CHANNEL_LIST } return os; } #undef X #endif diff --git a/src/exception/exception.hxx b/src/exception/exception.hxx index 294ee46..ec72518 100644 --- a/src/exception/exception.hxx +++ b/src/exception/exception.hxx @@ -1,40 +1,37 @@ -#ifndef EXCEPTION_EXCEPTION_SEEN -#define EXCEPTION_EXCEPTION_SEEN +#pragma once #include #include #include namespace nuis { struct nuis_except : public std::exception { std::stringstream msgstrm; std::string msg; nuis_except() : msgstrm(), msg() {} nuis_except(nuis_except const &other) : msgstrm(), msg() { msgstrm << other.msg; msg = other.msg; } const char *what() const noexcept { return msg.c_str(); } template nuis_except &operator<<(T const &obj) { msgstrm << obj; msg = msgstrm.str(); return (*this); } }; } // namespace nuis #define NEW_NUIS_EXCEPT(EXCEPT_NAME) \ struct EXCEPT_NAME : public nuis::nuis_except { \ EXCEPT_NAME() : nuis::nuis_except() {} \ EXCEPT_NAME(EXCEPT_NAME const &other) : nuis::nuis_except(other) {} \ template EXCEPT_NAME &operator<<(T const &obj) { \ msgstrm << obj; \ msg = msgstrm.str(); \ return (*this); \ } \ } - -#endif diff --git a/src/generator/input/GENIEInputHandler.cxx b/src/generator/input/GENIEInputHandler.cxx index baa784a..cc5c146 100644 --- a/src/generator/input/GENIEInputHandler.cxx +++ b/src/generator/input/GENIEInputHandler.cxx @@ -1,149 +1,165 @@ #include "generator/input/GENIEInputHandler.hxx" #include "generator/utility/GENIEUtility.hxx" #include "utility/PDGCodeUtility.hxx" #ifdef GENIE_V3_INTERFACE #include "Framework/EventGen/EventRecord.h" #include "Framework/GHEP/GHepParticle.h" #include "Framework/GHEP/GHepRecord.h" #else #include "EVGCore/EventRecord.h" #include "GHEP/GHepParticle.h" #include "GHEP/GHepRecord.h" #endif #include "fhiclcpp/ParameterSet.h" using namespace nuis; using namespace nuis::event; using namespace nuis::utility; using namespace nuis::genietools; -GENIEInputHandler::GENIEInputHandler() : fInputTree(), fGenieNtpl(nullptr) {} +GENIEInputHandler::GENIEInputHandler() + : fInputTree(), fGenieNtpl(nullptr), fFileWeight(1) {} GENIEInputHandler::GENIEInputHandler(GENIEInputHandler &&other) : fInputTree(std::move(other.fInputTree)), - fReaderEvent(std::move(other.fReaderEvent)), fGenieNtpl(nullptr) {} + fReaderEvent(std::move(other.fReaderEvent)), fGenieNtpl(nullptr), + fFileWeight(1) {} void GENIEInputHandler::Initialize(fhicl::ParameterSet const &ps) { fInputTree = CheckGetTTree(ps.get("file"), "gtree"); fInputTree.tree->SetBranchAddress("gmcrec", &fGenieNtpl); fKeepIntermediates = ps.get("keep_intermediates", false); fKeepNuclearParticles = ps.get("keep_nuclear_particles", false); + + fhicl::ParameterSet XSecInfo = ps.get("xsec_info", {}); + if (XSecInfo.has_key("weight")) { + fFileWeight = XSecInfo.get("weight"); + std::cout << "[INFO]: Average GENIE XSecWeight = " << fFileWeight << "" + << std::endl; + } else if (XSecInfo.has_key("flux") || XSecInfo.has_key("target") || + XSecInfo.has_key("spline_file")) { + std::cout + << "[INFO]: Attempting to build GENIE file weight with input splines." + << std::endl; + fFileWeight = genietools::GetFileWeight(XSecInfo) / double(GetNEvents()); + std::cout << "[INFO]: Average GENIE XSecWeight = " << fFileWeight << "" + << std::endl; + } } MinimalEvent const &GENIEInputHandler::GetMinimalEvent(ev_index_t idx) const { if (idx >= GetNEvents()) { throw IInputHandler::invalid_entry() << "[ERROR]: Attempted to get entry " << idx << " from an InputHandler with only " << GetNEvents(); } fInputTree.tree->GetEntry(idx); genie::GHepRecord *GHep = static_cast(fGenieNtpl->event); if (!GHep) { throw invalid_GENIE_event() << "[ERROR]: GENIE event " << idx << " failed to contain a GHepRecord"; } fReaderEvent.mode = GetEventChannel(*GHep); TObjArrayIter iter(GHep); genie::GHepParticle *p; while ((p = (dynamic_cast((iter).Next())))) { if (!p) { continue; } Particle::Status_t state = GetParticleStatus(*p, fReaderEvent.mode); if (state != Particle::Status_t::kPrimaryInitialState) { continue; } if (!IsNeutralLepton(p->Pdg(), pdgcodes::kMatterAntimatter) && !IsChargedLepton(p->Pdg(), pdgcodes::kMatterAntimatter)) { continue; } fReaderEvent.probe_E = p->E() * 1.E3; fReaderEvent.probe_pdg = p->Pdg(); break; } - fReaderEvent.XSecWeight = 1; + fReaderEvent.XSecWeight = fFileWeight; if (fWeightCache.size() <= idx) { fWeightCache.push_back(fReaderEvent.XSecWeight); } return fReaderEvent; } FullEvent const &GENIEInputHandler::GetFullEvent(ev_index_t idx) const { (void)GetMinimalEvent(idx); fReaderEvent.ClearParticleStack(); genie::GHepRecord *GHep = static_cast(fGenieNtpl->event); unsigned int npart = GHep->GetEntries(); // Fill Particle Stack genie::GHepParticle *p = 0; TObjArrayIter iter(GHep); // Loop over all particles while ((p = (dynamic_cast((iter).Next())))) { if (!p) { continue; } // Get Status Particle::Status_t state = GetParticleStatus(*p, fReaderEvent.mode); - if (!fKeepIntermediates) { + if (!fKeepIntermediates && (state == Particle::Status_t::kIntermediate)) { continue; } if (!fKeepNuclearParticles && IsNuclearPDG(p->Pdg())) { continue; } Particle nuis_part; - nuis_part.pdg = p->Pdg(); nuis_part.P4 = TLorentzVector(p->Px(), p->Py(), p->Pz(), p->E()) * 1E3; fReaderEvent.ParticleStack[static_cast(state)].particles.push_back( nuis_part); } return fReaderEvent; } double GENIEInputHandler::GetEventWeight(ev_index_t idx) const { if (idx > fWeightCache.size()) { throw weight_cache_miss() << "[ERROR]: Failed to get cached weight for event index: " << idx; } return fWeightCache[idx]; } double GENIEInputHandler::GetXSecScaleFactor( std::pair const &EnuRange) const { throw input_handler_feature_unimplemented() << "[ERROR]: Flux cuts not yet implemented for GENIE input handler."; } size_t GENIEInputHandler::GetNEvents() const { return fInputTree.tree->GetEntries(); } GeneratorManager::Generator_id_t GENIEInputHandler::GetGeneratorId() const { return GeneratorManager::Get().EnsureGeneratorRegistered("GENIE"); } DECLARE_PLUGIN(IInputHandler, GENIEInputHandler); diff --git a/src/generator/input/GENIEInputHandler.hxx b/src/generator/input/GENIEInputHandler.hxx index 7c87f17..312b6b7 100644 --- a/src/generator/input/GENIEInputHandler.hxx +++ b/src/generator/input/GENIEInputHandler.hxx @@ -1,67 +1,69 @@ // Copyright 2018 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 . *******************************************************************************/ #pragma once #include "event/FullEvent.hxx" #include "input/IInputHandler.hxx" #include "exception/exception.hxx" #include "utility/ROOTUtility.hxx" #ifdef GENIE_V3_INTERFACE #include "Framework/Ntuple/NtpMCEventRecord.h" #else #include "Ntuple/NtpMCEventRecord.h" #endif -#include - namespace fhicl { -class ParameterSet; + class ParameterSet; } +#include + class GENIEInputHandler : public IInputHandler { mutable nuis::utility::TreeFile fInputTree; mutable nuis::event::FullEvent fReaderEvent; mutable std::vector fWeightCache; mutable genie::NtpMCEventRecord *fGenieNtpl; bool fKeepIntermediates; bool fKeepNuclearParticles; + double fFileWeight; + public: NEW_NUIS_EXCEPT(weight_cache_miss); GENIEInputHandler(); GENIEInputHandler(GENIEInputHandler const &) = delete; GENIEInputHandler(GENIEInputHandler &&); void Initialize(fhicl::ParameterSet const &); nuis::event::MinimalEvent const &GetMinimalEvent(ev_index_t idx) const; nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const; double GetEventWeight(ev_index_t idx) const; size_t GetNEvents() const; double GetXSecScaleFactor(std::pair const &EnuRange) const; nuis::GeneratorManager::Generator_id_t GetGeneratorId() const; }; diff --git a/src/generator/input/NEUTInputHandler.cxx b/src/generator/input/NEUTInputHandler.cxx index ea62b19..9fc688d 100644 --- a/src/generator/input/NEUTInputHandler.cxx +++ b/src/generator/input/NEUTInputHandler.cxx @@ -1,257 +1,260 @@ #include "generator/input/NEUTInputHandler.hxx" #include "persistency/ROOTOutput.hxx" #include "utility/HistogramUtility.hxx" #include "utility/InteractionChannelUtility.hxx" #include "utility/PDGCodeUtility.hxx" #include "generator/utility/NEUTUtility.hxx" #include "neutpart.h" #include "neutvect.h" #include "fhiclcpp/ParameterSet.h" using namespace nuis; using namespace nuis::event; using namespace nuis::utility; using namespace nuis::neuttools; NEUTInputHandler::NEUTInputHandler() : fInputTreeFile() {} NEUTInputHandler::NEUTInputHandler(NEUTInputHandler &&other) : fInputTreeFile(std::move(other.fInputTreeFile)), fReaderEvent(std::move(other.fReaderEvent)) {} void NEUTInputHandler::Initialize(fhicl::ParameterSet const &ps) { fInputTreeFile = CheckGetTTree(ps.get("file"), "neuttree"); fNeutVect = nullptr; fInputTreeFile.tree->SetBranchAddress("vectorbranch", &fNeutVect); fKeepIntermediates = ps.get("keep_intermediates", false); std::string flux_name = ps.get("flux_hist_name", "flux_numu"); std::string evtrt_name = ps.get("evtrt_hist_name", "evtrt_numu"); fXSecRescaleFactor = ps.get("xsec_rescale_factor",1); std::string override_flux_file = ps.get("override_flux_file", ""); if (override_flux_file.size()) { fFlux = GetHistogramFromROOTFile(override_flux_file, flux_name); RebuildEventRate(!ps.get("flux_hist_in_MeV", false)); } else { fFlux = GetHistogramFromROOTFile(fInputTreeFile.file, flux_name, false); } if (fFlux) { if (!fEvtrt) { fEvtrt = GetHistogramFromROOTFile(fInputTreeFile.file, evtrt_name); } fFileWeight = fEvtrt->Integral() * 1.0E-38 / (fFlux->Integral() * double(GetNEvents())); std::cout << "[INFO]: Average NEUT XSecWeight = " << fFileWeight << "" << std::endl; } else { std::cout << "[INFO]: Input NEUT file doesn't have flux information, assuming " "mono energetic and calculating average cross-section..." << std::endl; fFileWeight = GetMonoEXSecWeight() / double(GetNEvents()); std::cout << "[INFO]: Done (Average NEUT XSecWeight = " << fFileWeight << ")!" << std::endl; } fFileWeight *= fXSecRescaleFactor; } NEW_NUIS_EXCEPT(non_mono_energetic_input_file); double NEUTInputHandler::GetMonoEXSecWeight() { double xsec = 0; double count = 0; double E_first = std::numeric_limits::max(); size_t NEvents = GetNEvents(); size_t ShoutEvery = NEvents / 100; std::cout << "[INFO]: Read " << 0 << "/" << NEvents << " NEUT events." << std::flush; for (size_t nev_it = 0; nev_it < NEvents; ++nev_it) { if (ShoutEvery && !(nev_it % ShoutEvery)) { std::cout << "\r[INFO]: Read " << nev_it << "/" << NEvents << " NEUT events." << std::flush; } fInputTreeFile.tree->GetEntry(nev_it); xsec += fNeutVect->Totcrs; + //Assume monoE. + return xsec * 1E-38; + count++; if (E_first == std::numeric_limits::max()) { NeutPart *part = fNeutVect->PartInfo(0); E_first = part->fP.E(); } else { NeutPart *part = fNeutVect->PartInfo(0); if (E_first != part->fP.E()) { throw non_mono_energetic_input_file() << "[ERROR]: When calculating NEUT xsec weight for a mono " "energetic beam, found different energies: " << E_first << " and " << part->fP.E() << ". Cannot continue with this input file."; } } } std::cout << "\r[INFO]: Read " << NEvents << "/" << NEvents << " NEUT events." << std::endl; return (xsec / count) * 1E-38; } void NEUTInputHandler::RebuildEventRate(bool FluxInGeV) { auto XSec = Clone(fFlux, true, "xsec"); auto Entry = Clone(fFlux, true, "entry"); std::cout << "[INFO]: Rebuilding total cross-section prediction..." << std::endl; size_t NEvents = GetNEvents(); size_t ShoutEvery = NEvents / 100; std::cout << "[INFO]: Read " << 0 << "/" << NEvents << " NEUT events." << std::flush; for (size_t nev_it = 0; nev_it < NEvents; ++nev_it) { if (ShoutEvery && !(nev_it % ShoutEvery)) { std::cout << "\r[INFO]: Read " << nev_it << "/" << NEvents << " NEUT events." << std::flush; } fInputTreeFile.tree->GetEntry(nev_it); NeutPart *part = fNeutVect->PartInfo(0); double E = part->fP.E() * (FluxInGeV ? 1E-3 : 1); XSec->Fill(E, fNeutVect->Totcrs); Entry->Fill(E); } std::cout << std::endl << "[INFO]: Done!" << std::endl; // Binned total xsec XSec->Divide(Entry.get()); fEvtrt = Clone(XSec, false, "evtrt"); // Event rate prediction fEvtrt->Multiply(fFlux.get()); } double NEUTInputHandler::GetXSecScaleFactor(std::pair const &) const { throw input_handler_feature_unimplemented() << "[ERROR]: Flux cuts not yet implemented for NEUT input handler."; } NeutVect *NEUTInputHandler::GetNeutEvent(ev_index_t idx) const { if (idx >= GetNEvents()) { throw IInputHandler::invalid_entry() << "[ERROR]: Attempted to get entry " << idx << " from an InputHandler with only " << GetNEvents(); } fInputTreeFile.tree->GetEntry(idx); return fNeutVect; } MinimalEvent const &NEUTInputHandler::GetMinimalEvent(ev_index_t idx) const { if (idx >= GetNEvents()) { throw IInputHandler::invalid_entry() << "[ERROR]: Attempted to get entry " << idx << " from an InputHandler with only " << GetNEvents(); } fInputTreeFile.tree->GetEntry(idx); fReaderEvent.mode = IntToChannel(fNeutVect->Mode); size_t NPart = fNeutVect->Npart(); for (size_t part_it = 0; part_it < NPart; part_it++) { NeutPart const &part = (*fNeutVect->PartInfo(part_it)); if ((part.fIsAlive == false) && (part.fStatus == -1) && IsNeutralLepton(part.fPID)) { fReaderEvent.probe_E = part.fP.T(); fReaderEvent.probe_pdg = part.fPID; break; } } fReaderEvent.XSecWeight = fFileWeight; if (fWeightCache.size() <= idx) { fWeightCache.push_back(fReaderEvent.XSecWeight); } fReaderEvent.fGenEvent = fNeutVect; return fReaderEvent; } FullEvent const &NEUTInputHandler::GetFullEvent(ev_index_t idx) const { (void)GetMinimalEvent(idx); fReaderEvent.ClearParticleStack(); if (fNeutVect->Ibound) { fReaderEvent.target_pdg = MakeNuclearPDG(fNeutVect->TargetA, fNeutVect->TargetZ); } else { fReaderEvent.target_pdg = MakeNuclearPDG(1, 1); } size_t NPart = fNeutVect->Npart(); size_t NPrimary = fNeutVect->Nprimary(); for (size_t part_it = 0; part_it < NPart; part_it++) { NeutPart const &part = (*fNeutVect->PartInfo(part_it)); Particle nuis_part; nuis_part.pdg = part.fPID; nuis_part.P4 = part.fP; Particle::Status_t state = GetNeutParticleStatus(part, fReaderEvent.mode); if (!fKeepIntermediates && (state == Particle::Status_t::kIntermediate)) { continue; } size_t state_int = static_cast(state); // Add status == 0 particles as pre-FSI particles until we find an // intermediate state particle if ((part_it < NPrimary) && (state != Particle::Status_t::kPrimaryInitialState)) { fReaderEvent .ParticleStack[static_cast( Particle::Status_t::kPrimaryFinalState)] .particles.push_back(nuis_part); } fReaderEvent.ParticleStack[state_int].particles.push_back(nuis_part); } return fReaderEvent; } double NEUTInputHandler::GetEventWeight(ev_index_t idx) const { if (idx > fWeightCache.size()) { throw weight_cache_miss() << "[ERROR]: Failed to get cached weight for event index: " << idx; } return fWeightCache[idx]; } size_t NEUTInputHandler::GetNEvents() const { return fInputTreeFile.tree->GetEntries(); } GeneratorManager::Generator_id_t NEUTInputHandler::GetGeneratorId() const { return GeneratorManager::Get().EnsureGeneratorRegistered("NEUT"); } DECLARE_PLUGIN(IInputHandler, NEUTInputHandler); diff --git a/src/generator/utility/CMakeLists.txt b/src/generator/utility/CMakeLists.txt index 00e5686..bca56b8 100644 --- a/src/generator/utility/CMakeLists.txt +++ b/src/generator/utility/CMakeLists.txt @@ -1,57 +1,63 @@ if(USE_NUWRO) LIST(APPEND GENERATOR_UTILS_IMPL NuWroUtility.cxx) LIST(APPEND GENERATOR_UTILS_HDR NuWroUtility.hxx) endif(USE_NUWRO) if(USE_NEUT) LIST(APPEND GENERATOR_UTILS_IMPL NEUTUtility.cxx) LIST(APPEND GENERATOR_UTILS_HDR NEUTUtility.hxx) endif(USE_NEUT) if(USE_GENIE) - LIST(APPEND GENERATOR_UTILS_IMPL GENIEUtility.cxx) - LIST(APPEND GENERATOR_UTILS_HDR GENIEUtility.hxx) + LIST(APPEND GENERATOR_UTILS_IMPL + GENIEUtility.cxx + GENIESplineReader.cxx + ) + LIST(APPEND GENERATOR_UTILS_HDR + GENIEUtility.hxx + GENIESplineReader.hxx + ) endif(USE_GENIE) if(GENERATOR_UTILS_IMPL) add_library(nuis_generator_utility SHARED ${GENERATOR_UTILS_IMPL}) if(USE_GENIE) target_compile_options(nuis_generator_utility PUBLIC ${GENIE_CXX_FLAGS}) include_directories(${GENIE_INCLUDE_DIRS}) target_link_libraries(nuis_generator_utility ${GENIE_LINK_DIRS}) target_link_libraries(nuis_generator_utility -Wl,--no-as-needed ${GENIE_LIBS}) endif() if(USE_NEUT) target_compile_options(nuis_generator_utility PUBLIC ${NEUT_CXX_FLAGS}) include_directories(${NEUT_INCLUDE_DIRS}) target_link_libraries(nuis_generator_utility ${NEUT_LINK_DIRS}) target_link_libraries(nuis_generator_utility -Wl,--no-as-needed ${NEUT_LIBS}) target_link_libraries(nuis_generator_utility ${NEUT_IMPORTED_TARGETS}) endif() if(USE_NUWRO) target_compile_options(nuis_generator_utility PUBLIC ${NUWRO_CXX_FLAGS}) include_directories(${NUWRO_INCLUDE_DIRS}) target_link_libraries(nuis_generator_utility -Wl,--no-as-needed ${NUWRO_IMPORTED_TARGETS}) endif() if(${NEED_PYTHIA6}) target_link_libraries(nuis_generator_utility -L${PYTHIA6}) endif() target_link_libraries(nuis_generator_utility -L${ROOT_LIBDIR} ${ROOT_LIBS}) install(TARGETS nuis_generator_utility DESTINATION lib) endif() if(GENERATOR_UTILS_HDR) install(FILES ${GENERATOR_UTILS_HDR} DESTINATION include/generator/utility) endif() LIST(APPEND GENERATOR_DEPENDENT_TARGETS nuis_generator_utility) SET(GENERATOR_DEPENDENT_TARGETS ${GENERATOR_DEPENDENT_TARGETS} PARENT_SCOPE) diff --git a/src/generator/utility/GENIESplineReader.cxx b/src/generator/utility/GENIESplineReader.cxx new file mode 100644 index 0000000..9f9657c --- /dev/null +++ b/src/generator/utility/GENIESplineReader.cxx @@ -0,0 +1,483 @@ +#include "GENIESplineReader.hxx" + +#include "exception/exception.hxx" + +#include "utility/StringUtility.hxx" + +#include "string_parsers/from_string.hxx" +#include "string_parsers/utility.hxx" + +#include + +#include +#include + +// #define SPLINE_NAME_MATCHER_VERBOSE + +namespace nuis { +namespace genietools { + +SSAX::SSAX() + : EventOnAllTags(true), LeavingAlready(false), InterestingTags(), + buf_size(1), buf_usage(0), buf(NULL), depth(0) { + buf = new char[buf_size]; +}; + +SSAX::~SSAX() { delete buf; } + +void SSAX::OnStartDocument() {} +void SSAX::OnEndDocument() {} +void SSAX::OnStartElement( + const char *name, + std::vector> &attributes) { + std::cout << "[SSAX]: OnStartElement(\"" << name << "\", " << std::endl; + for (size_t a_it = 0; a_it < attributes.size(); ++a_it) { + std::cout << "\t\"" << attributes[a_it].first << "\" = \"" + << attributes[a_it].second << "\"" << std::endl; + } + std::cout << "\t)" << std::endl; +} +void SSAX::OnCharacters(const char *characters) { + std::cout << "[SSAX]: OnCharacters(\"" << characters << "\")" << std::endl; +} +void SSAX::OnEndElement(const char *name) { + std::cout << "[SSAX]: OnEndElement(\"" << name << "\")" << std::endl; +} +void SSAX::OnComment(const char *) {} + +void SSAX::NotifyInterestingTags(std::vector it) { + InterestingTags = it; + EventOnAllTags = !InterestingTags.size(); +} +void SSAX::NotifyEarlyExit() { LeavingAlready = true; } + +void SSAX::ExpandBuf() { + size_t new_buf_size = (buf_size <= 1E6) ? buf_size * 10 : (buf_size + 1E6); + + std::cout << "[SSAX]: Expanding internal buffer to: " << new_buf_size + << " bytes." << std::endl; + + char *new_buf = new char[new_buf_size]; + memcpy(new_buf, buf, buf_size); + delete[] buf; + + buf = new_buf; + buf_size = new_buf_size; +} + +void SSAX::AddToBuf(char character) { + if (buf_usage == buf_size) { + ExpandBuf(); + } + buf[buf_usage++] = character; +} + +void SSAX::NullTermBuf() { AddToBuf('\0'); } + +bool SSAX::IsNameStartChar(char c) { + return (isalpha(c) || (c == ':') || (c == '_')); +} + +bool SSAX::IsNameChar(char c) { + return (IsNameStartChar(c) || isdigit(c) || (c == '-') || (c == '.')); +} + +void SSAX::IgnoreAllToRightBracket(std::ifstream &ifs) { + char prev, next; + while (EOF != (next = ifs.get())) { + if (next == '>') { + if (prev == '/') { // Empty element + depth--; + } + + return; + } + prev = next; + } + std::cerr << "[SSAX - ERROR]: encountered EOF before ending '>' at depth: " + << depth << std::endl; + throw; +} + +char SSAX::IgnoreSpace(std::ifstream &ifs) { + char next; + while (EOF != (next = ifs.get())) { + if (isspace(next)) { + continue; + } else { + return next; + } + } + std::cerr << "[SSAX - ERROR]: encountered EOF early at depth: " << depth + << std::endl; + throw; +} + +char SSAX::BufNameToBoundary(std::ifstream &ifs) { + char next; + while (EOF != (next = ifs.get())) { + if (isspace(next) || (next == '>') || (next == '/') || (next == '=')) { + return next; + } + if (IsNameChar(next)) { + AddToBuf(next); + } else { + std::cerr << "[SSAX - ERROR]: encountered bad character " << next + << std::endl; + throw; + } + } + std::cerr << "[SSAX - ERROR]: encountered EOF early at depth: " << depth + << std::endl; + throw; +} + +char SSAX::BufAttrValueToMatchingQuote(std::ifstream &ifs, char quot) { + char next; + while (EOF != (next = ifs.get())) { + if (next == quot) { + return next; + } + AddToBuf(next); + } + std::cerr << "[SSAX - ERROR]: encountered EOF early at depth: " << depth + << std::endl; + throw; +} + +std::pair +SSAX::HandleAttribute(std::ifstream &ifs) { + return std::pair(nullptr, nullptr); +} + +void SSAX::HandleOpenTag(std::ifstream &ifs) { + char next = ifs.get(); + size_t AttrNameStart = 0; + size_t AttrValueStart = 0; + attrs.clear(); + bool IsInteresting = true; + + if (next == EOF) { + std::cerr << "[SSAX - ERROR]: encountered EOF directly after element " + "opening at depth: " + << depth << std::endl; + throw; + } + + if (!IsNameStartChar(next)) { + std::cerr << "[SSAX - ERROR]: Bad element name first character: " << next + << std::endl; + throw; + } + AddToBuf(next); + + next = BufNameToBoundary(ifs); + NullTermBuf(); + AttrNameStart = buf_usage; + + if (!EventOnAllTags) { + IsInteresting = false; + for (size_t it_it = 0; it_it < InterestingTags.size(); ++it_it) { + if (InterestingTags[it_it] == buf) { + IsInteresting = true; + break; + } + } + } + + if ((next == '>')) { // Handles elements like + if (IsInteresting) { + OnStartElement(buf, attrs); + } + depth++; + return; + } + if ((next == '/')) { // Handles elements like + next = ifs.get(); + if (EOF == next) { + std::cerr << "[SSAX - ERROR]: encountered EOF early at depth: " << depth + << std::endl; + throw; + } + + if ((next == '>')) { + if (IsInteresting) { + OnStartElement(buf, attrs); + } + return; + } else { + std::cerr << "[SSAX - ERROR]: encountered unexpected character " << next + << " after /. Expected >. Exiting early at depth: " << depth + << std::endl; + throw; + } + } + + while (EOF != (next = ifs.get())) { + if (isspace(next)) { + next = IgnoreSpace(ifs); + } + + if (IsNameStartChar(next)) { + AddToBuf(next); + next = BufNameToBoundary(ifs); + if (next != '=') { + std::cerr << "[SSAX - ERROR]: encountered unexpected character \"" + << next + << "\". Expected \"=\". Exiting early at depth: " << depth + << std::endl; + throw; + } + NullTermBuf(); + AttrValueStart = buf_usage; + + next = ifs.get(); + if ((next != '\"') && (next != '\'')) { + std::cerr << "[SSAX - ERROR]: encountered unexpected character \"" + << next << "\" (" << int(next) << "). Expected \"\'\"(" + << int('\'') << ") or \"\"\"(" << int('\"') + << "). Exiting early at depth: " << depth << std::endl; + throw; + } + next = BufAttrValueToMatchingQuote(ifs, next); + NullTermBuf(); + attrs.push_back( + std::make_pair(&buf[AttrNameStart], &buf[AttrValueStart])); + AttrNameStart = buf_usage; + continue; + } + + if ((next == '>')) { // Handles elements like + if (IsInteresting) { + OnStartElement(buf, attrs); + } + depth++; + return; + } + if ((next == '/')) { // Handles elements like + next = ifs.get(); + if (EOF == next) { + std::cerr << "[SSAX - ERROR]: encountered EOF early at depth: " << depth + << std::endl; + throw; + } + + if ((next == '>')) { + if (IsInteresting) { + OnStartElement(buf, attrs); + } + // No depth ++; + return; + } else { + std::cerr << "[SSAX - ERROR]: encountered unexpected character " << next + << " after /. Expected >. Exiting early at depth: " << depth + << std::endl; + throw; + } + } + + std::cerr << "[SSAX - ERROR]: encountered unexpected character " << next + << " after /. Expected >. Exiting early at depth: " << depth + << std::endl; + throw; + } + std::cerr << "[SSAX - ERROR]: encountered EOF early at depth: " << depth + << std::endl; + throw; +} + +void SSAX::HandleCloseTag(std::ifstream &ifs) { + ifs.ignore(); + char next = ifs.get(); + bool IsInteresting = true; + + if (!IsNameStartChar(next)) { + std::cerr << "[SSAX - ERROR]: Bad element name first character: " << next + << std::endl; + throw; + } + AddToBuf(next); + + next = BufNameToBoundary(ifs); + NullTermBuf(); + + if (!EventOnAllTags) { + IsInteresting = false; + for (size_t it_it = 0; it_it < InterestingTags.size(); ++it_it) { + if (InterestingTags[it_it] == buf) { + IsInteresting = true; + break; + } + } + } + + if (IsInteresting) { + OnEndElement(buf); + } +} + +void SSAX::HandleCommentTag(std::ifstream &ifs) { + IgnoreAllToRightBracket(ifs); +} + +void SSAX::FlushBufferToCharacterHandler() { + NullTermBuf(); + OnCharacters(buf); +} +void SSAX::ClearBuf() { buf_usage = 0; } + +void SSAX::ParseFile(const char *filename) { + std::ifstream ifs(filename); + + if (!ifs.good()) { + std::cerr << "[SSAX]: Failed to open " << filename << " for reading." + << std::endl; + throw; + } + ifs.clear(); + ifs.seekg(0); + + OnStartDocument(); + + char next; + while (EOF != (next = ifs.get())) { + if (LeavingAlready) { + return; + } + if (next == '<') { + FlushBufferToCharacterHandler(); + ClearBuf(); + switch (ifs.peek()) { + case '!': + case '?': { + HandleCommentTag(ifs); + ClearBuf(); + break; + } + case '/': { + HandleCloseTag(ifs); + ClearBuf(); + break; + } + default: { + HandleOpenTag(ifs); + ClearBuf(); + } + } + continue; + } + AddToBuf(next); + } + OnEndDocument(); +} + +GENIESplineGetter::GENIESplineGetter(std::vector sps) : SSAX() { + SplinePatternStrings = sps; + for (auto const &sp : SplinePatternStrings) { + SplinePatterns.emplace_back(nuis::utility::DeGlobPattern(sp)); + Splines.push_back({}); + } + + ReadingDocument = false; + NSplines = 0; + + InSplineElement = false; + InEElement = false; + InxsecElement = false; +} + +void GENIESplineGetter::OnStartDocument() { + ReadingDocument = true; + std::cout << "[GENIESplineGetter]: Begin document..." << std::endl; +} +void GENIESplineGetter::OnEndDocument() { + ReadingDocument = false; + std::cout << "[GENIESplineGetter]: Finished document, read " << NSplines + << " splines." << std::endl; +} + +void GENIESplineGetter::OnStartElement( + const char *name, + std::vector> &attributes) { + std::string name_str = std::string(name); + + if (InSplineElement) { + if (name_str == "E") { + InEElement = true; + } + if (name_str == "xsec") { + InxsecElement = true; + } + } + + if (name_str == "spline") { + for (size_t attr_it = 0; attr_it < attributes.size(); ++attr_it) { + if (std::string(attributes[attr_it].first) == "name") { + std::string attr_val = std::string(attributes[attr_it].second); +#ifdef SPLINE_NAME_MATCHER_VERBOSE + std::cout << "[DEBUG]: Found spline named: " << attr_val << std::endl; +#endif + for (sp_it = 0; sp_it < SplinePatterns.size(); ++sp_it) { +#ifdef SPLINE_NAME_MATCHER_VERBOSE + std::cout << "[DEBUG]: Checking against " + << SplinePatternStrings[sp_it] << std::endl; +#endif + if (std::regex_match(attr_val, SplinePatterns[sp_it])) { // match +#ifdef SPLINE_NAME_MATCHER_VERBOSE + std::cout << "[DEBUG]: Matched! " << std::endl; +#endif + InSplineElement = true; + Splines[sp_it].push_back({{}, {}, attr_val}); + break; + } + } + return; + } + } + } +} + +void GENIESplineGetter::OnCharacters(const char *characters) { + if (InEElement) { + std::string ch(characters); + fhicl::string_parsers::trim(ch); + Splines[sp_it].back().EValues.push_back( + fhicl::string_parsers::str2T(ch)); + InEElement = false; + } + if (InxsecElement) { + std::string ch(characters); + fhicl::string_parsers::trim(ch); + Splines[sp_it].back().XSecValues.push_back( + fhicl::string_parsers::str2T(ch)); + InxsecElement = false; + } +} + +void GENIESplineGetter::OnEndElement(const char *name) { + if (InSplineElement && (std::string(name) == "spline")) { + std::cout << "[GENIESplineGetter]: Finished reading spline element: " + << Splines[sp_it].back().SplineName << std::endl; + InSplineElement = false; + NSplines++; + } +} + +void GENIESplineGetter::OnComment(const char *) {} + +std::vector> GENIESplineGetter::GetTGraphs() { + std::vector> graphs; + for (sp_it = 0; sp_it < SplinePatterns.size(); ++sp_it) { + graphs.push_back({}); + for (size_t spl_it = 0; spl_it < Splines[sp_it].size(); ++spl_it) { + graphs.back().emplace_back(Splines[sp_it][spl_it].EValues.size(), + Splines[sp_it][spl_it].EValues.data(), + Splines[sp_it][spl_it].XSecValues.data()); + graphs.back().back().SetName(Splines[sp_it][spl_it].SplineName.c_str()); + } + } + return graphs; +} + +} // namespace genietools +} // namespace nuis diff --git a/src/generator/utility/GENIESplineReader.hxx b/src/generator/utility/GENIESplineReader.hxx new file mode 100644 index 0000000..a3d798e --- /dev/null +++ b/src/generator/utility/GENIESplineReader.hxx @@ -0,0 +1,110 @@ +#pragma once + +#include "TGraph.h" + +#include +#include +#include +#include +#include + +namespace nuis { +namespace genietools { + +class SSAX { + bool EventOnAllTags; + bool LeavingAlready; + std::vector InterestingTags; + + size_t buf_size; + size_t buf_usage; + char *buf; + std::vector> attrs; + + size_t depth; + +public: + SSAX(); + + virtual ~SSAX(); + + virtual void OnStartDocument(); + virtual void OnEndDocument(); + virtual void OnStartElement( + const char *name, + std::vector> &attributes); + virtual void OnCharacters(const char *characters); + virtual void OnEndElement(const char *name); + virtual void OnComment(const char *); + + void NotifyInterestingTags(std::vector it); + void NotifyEarlyExit(); + +private: + void ExpandBuf(); + + void AddToBuf(char character); + void NullTermBuf(); + + bool IsNameStartChar(char c); + bool IsNameChar(char c); + + void IgnoreAllToRightBracket(std::ifstream &ifs); + + char IgnoreSpace(std::ifstream &ifs); + + char BufNameToBoundary(std::ifstream &ifs); + + char BufAttrValueToMatchingQuote(std::ifstream &ifs, char quot); + + std::pair HandleAttribute(std::ifstream &ifs); + + void HandleOpenTag(std::ifstream &ifs); + void HandleCloseTag(std::ifstream &ifs); + void HandleCommentTag(std::ifstream &ifs); + + void FlushBufferToCharacterHandler(); + + void ClearBuf(); + +public: + void ParseFile(const char *filename); +}; + +struct spline { + std::vector EValues; + std::vector XSecValues; + std::string SplineName; +}; + +class GENIESplineGetter : public SSAX { +public: + std::vector SplinePatternStrings; + std::vector SplinePatterns; + std::vector> Splines; + + bool InEElement; + bool InxsecElement; + bool InSplineElement; + + bool ReadingDocument; + + size_t sp_it; + size_t NSplines; + + GENIESplineGetter(std::vector sn); + + void OnStartDocument(); + void OnEndDocument(); + void OnStartElement( + const char *name, + std::vector> &attributes); + void OnCharacters(const char *characters); + void OnEndElement(const char *name); + void OnComment(const char *); + + std::vector> GetTGraphs(); +}; + +} // namespace genietools +} // namespace nuis diff --git a/src/generator/utility/GENIEUtility.cxx b/src/generator/utility/GENIEUtility.cxx index 2d93aea..b3be697 100644 --- a/src/generator/utility/GENIEUtility.cxx +++ b/src/generator/utility/GENIEUtility.cxx @@ -1,298 +1,465 @@ // Copyright 2018 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 "generator/utility/GENIEUtility.hxx" +#include "generator/utility/GENIESplineReader.hxx" +#include "utility/HistogramUtility.hxx" +#include "utility/StringUtility.hxx" +#include "utility/InteractionChannelUtility.hxx" #include "utility/PDGCodeUtility.hxx" #ifdef GENIE_V3_INTERFACE #include "Framework/GHEP/GHepParticle.h" #include "Framework/GHEP/GHepRecord.h" #include "Framework/GHEP/GHepUtils.h" #include "Framework/ParticleData/PDGCodes.h" #else #include "GHEP/GHepParticle.h" #include "GHEP/GHepRecord.h" #include "GHEP/GHepUtils.h" #include "PDG/PDGCodes.h" #endif +#include "fhiclcpp/ParameterSet.h" + +#include "string_parsers/from_string.hxx" + #include "TGraph.h" namespace nuis { using namespace event; namespace genietools { static std::map SplineCache; TGraph dum; TGraph const &GetGENIESpline(std::string const &SplineFile, std::string const &SplineIdentifier) { if (SplineCache.find(SplineFile + SplineIdentifier) != SplineCache.end()) { return SplineCache.find(SplineFile + SplineIdentifier)->second; } return dum; } struct NFSParticleCount { size_t NProton; size_t NNeutron; size_t NPip; size_t NPi0; size_t NPim; size_t NOther; }; NFSParticleCount CountPreFSIParticles(genie::GHepRecord const &ev) { // This code in this method is adapted from the GENIE source code found in // GHep/GHepUtils.cxx This method therefore carries the GENIE copyright // licence as copied below: // /// Copyright (c) 2003-2017, GENIE Neutrino MC Generator Collaboration /// For the full text of the license visit http://copyright.genie-mc.org /// or see $GENIE/LICENSE // genie::Target const &tgt = ev.Summary()->InitState().Tgt(); if (!tgt.HitNucIsSet()) { throw invalid_GENIE_event() << "[ERROR]: Failed to get hit nucleon kinematics as it was not " "included in this GHep event. This is a fatal error."; } genie::GHepParticle *FSLep = ev.FinalStatePrimaryLepton(); genie::GHepParticle *ISLep = ev.Probe(); if (!FSLep || !ISLep) { throw invalid_GENIE_event() << "[ERROR]: Failed to find IS and FS lepton in event: " << ev.Summary()->AsString(); } size_t NPi0 = 0, NPip = 0, NPim = 0, NProton = 0, NNeutron = 0, NOther = 0; bool nuclear_target = tgt.IsNucleus(); TIter event_iter(&ev); genie::GHepParticle *p = 0; while ((p = dynamic_cast(event_iter.Next()))) { genie::GHepStatus_t ghep_ist = (genie::GHepStatus_t)p->Status(); int ghep_pdgc = p->Pdg(); int ghep_fm = p->FirstMother(); int ghep_fmpdgc = (ghep_fm == -1) ? 0 : ev.Particle(ghep_fm)->Pdg(); // For nuclear targets use hadrons marked as 'hadron in the nucleus' // which are the ones passed in the intranuclear rescattering // For free nucleon targets use particles marked as 'final state' // but make an exception for decayed pi0's,eta's (count them and not their // daughters) bool decayed = (ghep_ist == genie::kIStDecayedState && (ghep_pdgc == genie::kPdgPi0 || ghep_pdgc == genie::kPdgEta)); bool parent_included = (ghep_fmpdgc == genie::kPdgPi0 || ghep_fmpdgc == genie::kPdgEta); bool count_it = (nuclear_target && ghep_ist == genie::kIStHadronInTheNucleus) || (!nuclear_target && decayed) || (!nuclear_target && ghep_ist == genie::kIStStableFinalState && !parent_included); if (!count_it) { continue; } if (ghep_pdgc == genie::kPdgPiP) { NPip++; } else if (ghep_pdgc == genie::kPdgPiM) { NPim++; } else if (ghep_pdgc == genie::kPdgPi0) { NPi0++; } else if (ghep_pdgc == genie::kPdgProton) { NProton++; } else if (ghep_pdgc == genie::kPdgNeutron) { NNeutron++; } else if (!utility::IsNeutralLepton( ghep_pdgc, utility::pdgcodes::kMatterAntimatter) && !utility::IsChargedLepton( ghep_pdgc, utility::pdgcodes::kMatterAntimatter)) { NOther++; } } return NFSParticleCount{NProton, NNeutron, NPip, NPi0, NPim, NOther}; } Channel_t GetEventChannel(genie::GHepRecord const &gev) { // Electron Scattering if (gev.Summary()->ProcInfo().IsEM()) { if (gev.Summary()->InitState().ProbePdg() == utility::pdgcodes::kElectron) { if (gev.Summary()->ProcInfo().IsQuasiElastic()) { NFSParticleCount fsparts = CountPreFSIParticles(gev); if (fsparts.NProton) { return Channel_t::kNCELP; } else { return Channel_t::kNCELN; } } else if (gev.Summary()->ProcInfo().IsMEC()) { return Channel_t::kNC2p2h; } else if (gev.Summary()->ProcInfo().IsResonant()) { NFSParticleCount fsparts = CountPreFSIParticles(gev); if (fsparts.NOther || ((fsparts.NPip + fsparts.NPi0 + fsparts.NPim) > 1)) { return Channel_t::kNCTransitionMPi; } else if (fsparts.NPip) { return Channel_t::kNCSPP_NPip; } else if (fsparts.NPi0) { return fsparts.NProton ? Channel_t::kNCSPP_PPi0 : Channel_t::kNCSPP_NPi0; } else if (fsparts.NPim) { return Channel_t::kNCSPP_PPim; } return Channel_t::kNCTransitionMPi; } else if (gev.Summary()->ProcInfo().IsDeepInelastic()) { return Channel_t::kNCDIS; } else { std::cout << "Unknown GENIE Electron Scattering Mode!" << std::endl << "ScatteringTypeId = " << gev.Summary()->ProcInfo().ScatteringTypeId() << " " << "InteractionTypeId = " << gev.Summary()->ProcInfo().InteractionTypeId() << std::endl << genie::ScatteringType::AsString( gev.Summary()->ProcInfo().ScatteringTypeId()) << " " << genie::InteractionType::AsString( gev.Summary()->ProcInfo().InteractionTypeId()) << " " << gev.Summary()->ProcInfo().IsMEC() << std::endl; return Channel_t::kUndefined; } } // Weak CC } else if (gev.Summary()->ProcInfo().IsWeakCC()) { // CC MEC if (gev.Summary()->ProcInfo().IsMEC()) { if (utility::IsNeutralLepton(gev.Summary()->InitState().ProbePdg(), utility::pdgcodes::kMatter)) { return Channel_t::kCC2p2h; } else if (utility::IsNeutralLepton(gev.Summary()->InitState().ProbePdg(), utility::pdgcodes::kAntimatter)) { return Channel_t::kCC2p2h_nub; } // CC OTHER } else { return FromNEUTCode(genie::utils::ghep::NeutReactionCode(&gev)); } // Weak NC } else if (gev.Summary()->ProcInfo().IsWeakNC()) { // NC MEC if (gev.Summary()->ProcInfo().IsMEC()) { if (utility::IsNeutralLepton(gev.Summary()->InitState().ProbePdg(), utility::pdgcodes::kMatter)) { return Channel_t::kNC2p2h; } else if (utility::IsNeutralLepton(gev.Summary()->InitState().ProbePdg(), utility::pdgcodes::kAntimatter)) { return Channel_t::kNC2p2h_nub; } // NC OTHER } else { return FromNEUTCode(genie::utils::ghep::NeutReactionCode(&gev)); } } return Channel_t::kUndefined; } Particle::Status_t GetParticleStatus(genie::GHepParticle const &p, - Channel_t chan) { + Channel_t chan) { /* 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 */ Particle::Status_t state = Particle::Status_t::kUnknown; switch (p.Status()) { case genie::kIStNucleonTarget: case genie::kIStInitialState: case genie::kIStCorrelatedNucleon: case genie::kIStNucleonClusterTarget: { state = Particle::Status_t::kPrimaryInitialState; break; } case genie::kIStStableFinalState: { state = Particle::Status_t::kNuclearLeaving; break; } case genie::kIStHadronInTheNucleus: { - state = Is2p2h(chan) ? Particle::Status_t::kPrimaryInitialState - : Particle::Status_t::kIntermediate; + state = utility::Is2p2h(chan) ? Particle::Status_t::kPrimaryInitialState + : Particle::Status_t::kIntermediate; break; } case genie::kIStPreDecayResonantState: case genie::kIStDISPreFragmHadronicState: case genie::kIStIntermediateState: { state = Particle::Status_t::kIntermediate; break; } case genie::kIStFinalStateNuclearRemnant: case genie::kIStUndefined: case genie::kIStDecayedState: default: { state = Particle::Status_t::kUnknown; } } if (utility::IsNuclearPDG(p.Pdg())) { if (state == Particle::Status_t::kPrimaryInitialState) { state = Particle::Status_t::kPrimaryInitialState; } else if (state == Particle::Status_t::kNuclearLeaving) { state = Particle::Status_t::kPrimaryFinalState; } } return state; } +struct TargetSplineBlob { + double MolecularWeight; + size_t NNucleons; + std::string search_term; + TGraph TotSpline; +}; + +NEW_NUIS_EXCEPT(Invalid_target_specifier); + +double GetFileWeight(fhicl::ParameterSet const &xsecinfo) { + double mec_scale = xsecinfo.get("mec_scale", 1); + double EMin = xsecinfo.get("EMin", 0); + double EMax = xsecinfo.get("EMax", 10); + size_t NKnots = xsecinfo.get("NKnots", 100); + + std::vector TargetSplines; + int probe = xsecinfo.get("nu_pdg"); + + for (std::string const &target : + xsecinfo.get>("target")) { + std::vector splits = nuis::utility::split(target, ";"); + if (splits.size() != 1 && splits.size() != 2) { + throw Invalid_target_specifier() + << "Expected to find 100ZZZAAA0;, e.g. " + "1000060120;1, " + "but found: \"" + << target << "\""; + } + std::stringstream ss; + ss << ".*nu:" << probe << ";tgt:" << splits[0] << ".*"; + + size_t nuc_pdg = fhicl::string_parsers::str2T(splits[0]); + + double molwght = splits.size() == 2 + ? fhicl::string_parsers::str2T(splits[1]) + : 1; + + TargetSplines.push_back( + {molwght, (nuc_pdg % 1000) / 10, ss.str(), TGraph()}); + + std::cout << "[INFO]:\tGetting splines that match: \"" << ss.str() + << "\" with A = " << TargetSplines.back().NNucleons + << " and multiplied by " + << molwght / double(TargetSplines.back().NNucleons) << std::endl; + } + + size_t NTotNucleons = 0; + std::vector spline_patterns; + for (auto const &spb : TargetSplines) { + spline_patterns.push_back(spb.search_term); + NTotNucleons += spb.MolecularWeight * spb.NNucleons; + } + double WeightToPerNucleon = 1.0 / double(NTotNucleons); + + // Read the spline file + std::cout << "[INFO]: Reading input file: " + << xsecinfo.get("spline_file") + << ", this may take a while..." << std::endl; + GENIESplineGetter spg(spline_patterns); + spg.ParseFile(xsecinfo.get("spline_file").c_str()); + std::vector> all_splines = spg.GetTGraphs(); + + double step = (EMax - EMin) / double(NKnots); + + // Sum each targets-worth of interaction splines together, weighting for + // molecular fraction and possibly fixing the MEC spline oddity + for (size_t ts_it = 0; ts_it < TargetSplines.size(); ++ts_it) { + + TargetSplines[ts_it].TotSpline = TGraph(NKnots); + + for (Int_t p_it = 0; p_it < NKnots; ++p_it) { + TargetSplines[ts_it].TotSpline.SetPoint(p_it, EMin + p_it * step, 0); + } + + double target_weight = (TargetSplines[ts_it].MolecularWeight / + double(TargetSplines[ts_it].NNucleons)); + for (size_t c_it = 0; c_it < all_splines[ts_it].size(); ++c_it) { + double weight = target_weight; + if (std::string(all_splines[ts_it][c_it].GetName()).find("MEC") != + std::string::npos) { + if (mec_scale != 1) { + std::cout << "[INFO]: Weighting MEC splines by " << mec_scale + << std::endl; + weight *= mec_scale; + } + } + for (Int_t p_it = 0; p_it < NKnots; ++p_it) { + double E, XSec; + TargetSplines[ts_it].TotSpline.GetPoint(p_it, E, XSec); + + // Copied from GENIE/Convents/Units.h + static const double gigaelectronvolt = 1.; + static const double GeV = gigaelectronvolt; + static const double meter = 5.07e+15 / GeV; + static const double centimeter = 0.01 * meter; + static const double centimeter2 = centimeter * centimeter; + + XSec += all_splines[ts_it][c_it].Eval(E) * weight / centimeter2; + + TargetSplines[ts_it].TotSpline.SetPoint(p_it, E, XSec); + } + } + } + + // Sum all the correctly weighted per-target splines + TGraph MasterSpline(NKnots); + + for (Int_t p_it = 0; p_it < NKnots; ++p_it) { + MasterSpline.SetPoint(p_it, EMin + p_it * step, 0); + } + for (size_t c_it = 0; c_it < TargetSplines.size(); ++c_it) { + for (Int_t p_it = 0; p_it < NKnots; ++p_it) { + double E, XSec; + MasterSpline.GetPoint(p_it, E, XSec); + XSec += TargetSplines[c_it].TotSpline.Eval(E) * WeightToPerNucleon; + + MasterSpline.SetPoint(p_it, E, XSec); + } + } + + fhicl::ParameterSet fluxps = xsecinfo.get("flux"); + // If MonoE + if (fluxps.has_key("energy_GeV")) { + double monoE = fluxps.get("energy_GeV"); + + return MasterSpline.Eval(monoE); + } else { + // Use the master spline to get flux*xsec (which the events are thrown + // according to) and return the flux-averaged total xsec which is used for + // weighting + std::unique_ptr Flux = + nuis::utility::GetHistogram(fluxps.get("histogram")); + bool per_width = fluxps.get("is_divided_by_bin_width", true); + double FluxIntegral = Flux->Integral(per_width ? "width" : ""); + + for (int bi_it = 0; bi_it < Flux->GetXaxis()->GetNbins(); ++bi_it) { + double bc = Flux->GetBinContent(bi_it + 1); + double low_e = Flux->GetXaxis()->GetBinLowEdge(bi_it + 1); + double up_e = Flux->GetXaxis()->GetBinUpEdge(bi_it + 1); + + size_t NIntSteps = 100; + double step = (up_e - low_e) / double(NIntSteps); + + double avg_xsec = 0; + for (size_t i = 0; i < NIntSteps; ++i) { + double e = low_e + (double(i) + 0.5) * step; + avg_xsec += MasterSpline.Eval(e); + } + avg_xsec /= double(NIntSteps); + + Flux->SetBinContent(bi_it + 1, bc * avg_xsec); + } + + double EvRateIntegral = Flux->Integral(per_width ? "width" : ""); + + return EvRateIntegral / FluxIntegral; + } +} + } // namespace genietools } // namespace nuis diff --git a/src/generator/utility/GENIEUtility.hxx b/src/generator/utility/GENIEUtility.hxx index b95f450..3bc447a 100644 --- a/src/generator/utility/GENIEUtility.hxx +++ b/src/generator/utility/GENIEUtility.hxx @@ -1,47 +1,53 @@ // Copyright 2018 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 . *******************************************************************************/ #pragma once #include "event/Particle.hxx" #include "event/types.hxx" class TGraph; namespace genie { class GHepRecord; class GHepParticle; } // namespace genie +namespace fhicl { + class ParameterSet; +} + namespace nuis { namespace genietools { NEW_NUIS_EXCEPT(invalid_GENIE_event); TGraph const &GetGENIESpline(std::string const &SplineFile, std::string const &SplineIdentifier); event::Channel_t GetEventChannel(genie::GHepRecord const &); event::Particle::Status_t GetParticleStatus(genie::GHepParticle const &p, - event::Channel_t chan); + event::Channel_t chan); + +double GetFileWeight(fhicl::ParameterSet const &xsecinfo); } // namespace genietools } // namespace nuis diff --git a/src/generator/utility/NEUTUtility.cxx b/src/generator/utility/NEUTUtility.cxx index 5ff5b77..0b0e04a 100644 --- a/src/generator/utility/NEUTUtility.cxx +++ b/src/generator/utility/NEUTUtility.cxx @@ -1,95 +1,95 @@ #include "generator/utility/NEUTUtility.hxx" #include "generator/input/NEUTInputHandler.hxx" #include "exception/exception.hxx" #include "utility/InteractionChannelUtility.hxx" #include "utility/PDGCodeUtility.hxx" #include "neutpart.h" #include "neutvect.h" using namespace nuis::event; using namespace nuis::utility; namespace nuis { namespace neuttools { NEW_NUIS_EXCEPT(unexpected_NEUT_particle_state); Particle::Status_t GetNeutParticleStatus(NeutPart const &part, Channel_t mode) { #ifdef DEBUG_NEUT_UTILITY std::cout << "[DEBUG]: Mode: " << mode << ", Part: { Status: " << part.fStatus << ", IsAlive: " << part.fIsAlive << ", PDG: " << part.fPID << " }." << std::endl; #endif // Remove Pauli blocked events, probably just single pion events if (part.fStatus == 5) { return Particle::Status_t::kBlocked; // fStatus == -1 means initial state } else if (part.fIsAlive == false && part.fStatus == -1) { return Particle::Status_t::kPrimaryInitialState; // NEUT has a bit of a strange convention for fIsAlive and fStatus // combinations // for NC and neutrino particle isAlive true/false and status 2 means // final state particle // for other particles in NC status 2 means it's an FSI particle // for CC it means it was an FSI particle } else if (part.fStatus == 2) { // NC case is a little strange... The outgoing neutrino might be alive or // not alive. Remaining particles with status 2 are FSI particles that // reinteracted if (IsNC(mode) && IsNeutralLepton(part.fPID)) { return Particle::Status_t::kNuclearLeaving; // The usual CC case } else if (part.fIsAlive == true) { // return Particle::Status_t::kIntermediate; throw unexpected_NEUT_particle_state() << "[ERROR] Found unexpected NEUT particle in neutvect stack: Mode: " << mode << " (IsNC: " << IsNC(mode) << "), Part: { Status: " << part.fStatus << ", IsAlive: " << part.fIsAlive << ", PDG: " << part.fPID << ", IsNeutralLepton: " << IsNeutralLepton(part.fPID) << " }."; } } else if ((part.fIsAlive == true) && (part.fStatus == 2) && IsNeutralLepton(part.fPID)) { return Particle::Status_t::kNuclearLeaving; } else if ((part.fIsAlive == true) && (part.fStatus == 0)) { return Particle::Status_t::kNuclearLeaving; } else if (!part.fIsAlive && ((part.fStatus == 1) || (part.fStatus == 3) || (part.fStatus == 4) || (part.fStatus == 7) || (part.fStatus == 8))) { return Particle::Status_t::kIntermediate; // There's one hyper weird case where fStatus = -3. This apparently // corresponds to a nucleon being ejected via pion FSI when there is "data // available" } else if (!part.fIsAlive && (part.fStatus == -3)) { return Particle::Status_t::kUnknown; // NC neutrino outgoing } else if (!part.fIsAlive && part.fStatus == 0 && IsNeutralLepton(part.fPID)) { return Particle::Status_t::kNuclearLeaving; // Warn if we still find alive particles without classifying them } else if (part.fIsAlive == true) { - throw unexpected_NEUT_particle_state() - << "[ERROR]: Undefined NEUT state " - << " Alive: " << part.fIsAlive << " Status: " << part.fStatus - << " PDG: " << part.fPID; + std::cout << "[WARN]: Undefined NEUT state " + << " Alive: " << part.fIsAlive << " Status: " << part.fStatus + << " PDG: " << part.fPID << std::endl; + return Particle::Status_t::kUnknown; } // Warn if we find dead particles that we haven't classified - throw unexpected_NEUT_particle_state() - << "[ERROR]: Undefined NEUT state " - << " Alive: " << part.fIsAlive << " Status: " << part.fStatus - << " PDG: " << part.fPID; + std::cout << "[WARN]: Undefined NEUT state " + << " Alive: " << part.fIsAlive << " Status: " << part.fStatus + << " PDG: " << part.fPID << std::endl; + return Particle::Status_t::kUnknown; } } // namespace neuttools } // namespace nuis diff --git a/src/persistency/ROOTOutput.cxx b/src/persistency/ROOTOutput.cxx index 2de36b1..34ac610 100644 --- a/src/persistency/ROOTOutput.cxx +++ b/src/persistency/ROOTOutput.cxx @@ -1,81 +1,88 @@ #include "persistency/ROOTOutput.hxx" #include "utility/ROOTUtility.hxx" #include "config/GlobalConfiguration.hxx" #include "fhiclcpp/ParameterSet.h" #include "TFile.h" namespace nuis { namespace persistency { struct NamedTFile { std::string name; std::unique_ptr file; NamedTFile() : name(""), file(nullptr) {} NamedTFile(NamedTFile &&other) : name(std::move(other.name)), file(std::move(other.file)) {} ~NamedTFile() { if (file) { file->Write(); file->Close(); } } }; static std::vector Files; -std::unique_ptr &GetOutputFile(std::string const &name) { - for (NamedTFile &file : Files) { - if (file.name == name) { - return file.file; - } - } - - fhicl::ParameterSet const &persistency = - config::GetDocument().get("persistency"); - std::string file_name = persistency.get(name + ".output_file"); - std::string open_opts = - persistency.get(name + ".open_mode", "CREATE"); - +void NewStream(std::string const &name, + std::string file_name, + std::string opts) { NamedTFile ntf; ntf.name = name; - ntf.file = std::make_unique(file_name.c_str(), open_opts.c_str()); + ntf.file = std::make_unique(file_name.c_str(), opts.c_str()); if (!ntf.file || !ntf.file->IsOpen()) { - if ((name == "default") && (open_opts == "CREATE")) { + if ((name == "default") && (opts == "CREATE")) { std::cout << "[WARN]: It appears the default file cannot be opened because it " "exists already, to stop you wasting processing time, the default " "output stream will be written to nuis.default.tmp.root in the " "current directory in RECREATE mode." << std::endl; file_name = "nuis.default.tmp.root"; - open_opts = "RECREATE"; - ntf.file = std::make_unique(file_name.c_str(), open_opts.c_str()); + opts = "RECREATE"; + ntf.file = std::make_unique(file_name.c_str(), opts.c_str()); if (!ntf.file || !ntf.file->IsOpen()) { throw utility::failed_to_open_TFile() << "[ERROR]: Failed to open output file: " << std::quoted(file_name) - << " in write mode with opts = " << std::quoted(open_opts); + << " in write mode with opts = " << std::quoted(opts); } } else { throw utility::failed_to_open_TFile() << "[ERROR]: Failed to open output file: " << std::quoted(file_name) - << " in write mode with opts = " << std::quoted(open_opts); + << " in write mode with opts = " << std::quoted(opts); } } Files.push_back(std::move(ntf)); - return Files.back().file; +} + +std::unique_ptr &GetOutputFile(std::string const &name) { + for (NamedTFile &file : Files) { + if (file.name == name) { + return file.file; + } + } + + fhicl::ParameterSet const &persistency = + config::GetDocument().get("persistency"); + std::string file_name = persistency.get(name + ".output_file"); + std::string opts = + persistency.get(name + ".open_mode", "CREATE"); + + NewStream(name, file_name, opts); + + return GetOutputFile(name); } void CloseOpenTFiles() { for (NamedTFile &f : Files) { std::cout << "[INFO]: Closing open TFile: " << f.name << " " << f.file->GetName() << std::endl; } Files.clear(); } } // namespace persistency } // namespace nuis diff --git a/src/persistency/ROOTOutput.hxx b/src/persistency/ROOTOutput.hxx index 1aa5bb6..f88da22 100644 --- a/src/persistency/ROOTOutput.hxx +++ b/src/persistency/ROOTOutput.hxx @@ -1,76 +1,80 @@ #ifndef PERSITENCY_ROOTOUTPUT_HXX_SEEN #define PERSITENCY_ROOTOUTPUT_HXX_SEEN #include "exception/exception.hxx" #include "TFile.h" #include #include #include namespace nuis { namespace persistency { NEW_NUIS_EXCEPT(WriteToOutputFile_nullptr); /// Will get/open a TFile that is described in the global config /// /// The named streams will be used to configure the file name and open mode from /// the global config element persistency.: {file: output.root opts: /// CREATE} std::unique_ptr &GetOutputFile(std::string const &name = "default"); +void NewStream(std::string const &name = "default", + std::string file_name = "default.nuis.root", + std::string opts = "CREATE"); + template inline void WriteToOutputFile(T *object, std::string const &object_name, std::string dir_name = "", std::string const &file_name = "default") { if (!object) { throw WriteToOutputFile_nullptr(); } TDirectory *ogdir = gDirectory; std::unique_ptr &f = GetOutputFile(file_name); TDirectory *d = f.get(); while (dir_name.length()) { size_t next_slash = dir_name.find_first_of('/'); std::string next_dir = dir_name.substr(0, next_slash); if (next_slash != std::string::npos) { dir_name = dir_name.substr(next_slash + 1); } else { dir_name = ""; } TDirectory *nd = d->GetDirectory(next_dir.c_str()); if (!nd) { nd = d->mkdir(next_dir.c_str()); } nd->cd(); d = nd; } d->cd(); object->Write(object_name.c_str(), TObject::kOverwrite); if (ogdir) { ogdir->cd(); } } template inline void WriteToOutputFile(std::unique_ptr &object, std::string const &object_name, std::string dir_name = "", std::string const &file_name = "default") { return WriteToOutputFile(object.get(), object_name, dir_name, file_name); } void CloseOpenTFiles(); } // namespace persistency } // namespace nuis #endif diff --git a/src/plugins/Instantiate.hxx b/src/plugins/Instantiate.hxx index 818d778..f4ce598 100644 --- a/src/plugins/Instantiate.hxx +++ b/src/plugins/Instantiate.hxx @@ -1,252 +1,253 @@ // Copyright 2018 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 . *******************************************************************************/ #pragma once #include "plugins/NamedSO.hxx" #include "plugins/traits.hxx" #include "config/GlobalConfiguration.hxx" #include "utility/FileSystemUtility.hxx" +#include "utility/StringUtility.hxx" #include "exception/exception.hxx" #include "fhiclcpp/ParameterSet.h" #include "string_parsers/to_string.hxx" // linux #include #include #include #include #include #include #include // #define DEBUG_INSTANTIATE namespace nuis { namespace plugins { extern std::vector LoadedSharedObjects; NEW_NUIS_EXCEPT(failed_to_find_instantiator); NEW_NUIS_EXCEPT(malformed_plugin_interface); NEW_NUIS_EXCEPT(failed_to_load_so); typedef void *(*inst_fcn)(); typedef void (*dltr_fcn)(void *); template struct PluginInstantiator { std::string FQ_so_path; std::string Base_classname; std::string Classname; void *dllib; inst_fcn Instantiator; dltr_fcn Deleter; PluginInstantiator() : FQ_so_path(""), Base_classname(""), Classname(""), dllib(nullptr), Instantiator(nullptr), Deleter(nullptr) {} PluginInstantiator(PluginInstantiator const &) = delete; PluginInstantiator(PluginInstantiator &&other) { FQ_so_path = std::move(other.FQ_so_path); Base_classname = std::move(other.Base_classname); Classname = std::move(other.Classname); dllib = other.dllib; Instantiator = other.Instantiator; Deleter = other.Deleter; other.FQ_so_path = ""; other.Base_classname = ""; other.Classname = ""; other.dllib = nullptr; other.Instantiator = nullptr; other.Deleter = nullptr; } typename plugin_traits::unique_ptr_t Instantiate() { T *inst = reinterpret_cast((*Instantiator)()); dltr_fcn dltr = Deleter; std::string cln = Classname; std::function deleter = [=](T *inst) { #ifdef DEBUG_INSTANTIATE std::cout << "[INFO]: Deleting instance of " << cln << " with " << (void *)dltr << std::endl; #endif (*dltr)(inst); }; return typename plugin_traits::unique_ptr_t(inst, deleter); } }; NamedSO &GetSharedObject(std::string const &FQPath) { for (NamedSO &so : LoadedSharedObjects) { if (so.name == FQPath) { return so; } } NamedSO so; so.name = FQPath; so.dllib = dlopen(FQPath.c_str(), RTLD_NOW | RTLD_GLOBAL); char const *dlerr_cstr = dlerror(); std::string dlerr; if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr.length()) { throw failed_to_load_so() << "[INFO]: Failed to load shared object: " << FQPath << " with dlerror: " << dlerr; } else { #ifdef DEBUG_INSTANTIATE std::cout << "[INFO]: Loaded shared object " << FQPath << std::endl; #endif } LoadedSharedObjects.push_back(std::move(so)); return LoadedSharedObjects.back(); } template typename plugin_traits::unique_ptr_t Instantiate(std::string const &classname) { static std::vector> LoadedPlugins; fhicl::ParameterSet const &plugins = config::GetDocument().get("plugins"); fhicl::ParameterSet const &search_paths = plugins.get("search_paths"); std::vector plugin_search_dirs; // Look for plugin search paths in sequence elements of the // plugins.search_paths table for (std::string const &key : search_paths.get_names()) { if (!search_paths.is_key_to_sequence(key)) { continue; } for (std::string const &path : search_paths.get>(key)) { plugin_search_dirs.push_back(path); } } for (std::string path : plugin_search_dirs) { path = utility::EnsureTrailingSlash(path); for (std::string const &so_name : utility::GetMatchingFiles(path, ".*\\.so")) { for (PluginInstantiator &plugin : LoadedPlugins) { if (plugin.FQ_so_path == (path + so_name) && (plugin.Base_classname == plugin_traits::interface_name()) && (plugin.Classname == classname)) { #ifdef DEBUG_INSTANTIATE std::cout << "[INFO]: Using already loaded PluginInstantiator" << std::endl; #endif return plugin.Instantiate(); } } PluginInstantiator plugin; plugin.FQ_so_path = path + so_name; plugin.Base_classname = plugin_traits::interface_name(); plugin.Classname = classname; plugin.dllib = GetSharedObject(plugin.FQ_so_path).dllib; char const *dlerr_cstr = nullptr; std::string dlerr(""); plugin.Instantiator = reinterpret_cast(dlsym( plugin.dllib, plugin_traits::instantiator_function_name(classname).c_str())); dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr_cstr) { #ifdef DEBUG_INSTANTIATE std::cout << "[INFO]: Failed to load appropriate instantiator method: " << plugin_traits::instantiator_function_name(classname) << " from shared object " << plugin.FQ_so_path << std::endl; #endif continue; } else { #ifdef DEBUG_INSTANTIATE std::cout << "[INFO]: Loaded instantiator method: " << plugin_traits::instantiator_function_name(classname) << " from shared object " << plugin.FQ_so_path << std::endl; #endif } plugin.Deleter = reinterpret_cast( dlsym(plugin.dllib, plugin_traits::deleter_function_name(classname).c_str())); dlerr_cstr = dlerror(); if (dlerr_cstr) { dlerr = dlerr_cstr; } if (dlerr_cstr) { throw malformed_plugin_interface() << "[ERROR]: Failed to load appropriate deleter method: " << plugin_traits::deleter_function_name(classname) << " from shared object " << plugin.FQ_so_path << " with error: " << std::quoted(dlerr); } else { #ifdef DEBUG_INSTANTIATE std::cout << "[INFO]: Loaded deleter method: " << plugin_traits::deleter_function_name(classname) << " from shared object " << plugin.FQ_so_path << std::endl; #endif } #ifdef DEBUG_INSTANTIATE std::cout << "[INFO]: Checking if shared object " << std::quoted(plugin.FQ_so_path) << " knows how to instantiate class " << std::quoted(classname) << " via interface " << std::quoted(plugin_traits::interface_name()) << std::endl; #endif LoadedPlugins.push_back(std::move(plugin)); return LoadedPlugins.back().Instantiate(); } } throw failed_to_find_instantiator() << "[ERROR]: Failed to find instantiator for classname: " << std::quoted(classname) << " using interface " << std::quoted(plugin_traits::interface_name()) << " from configured search paths: " << fhicl::string_parsers::T2Str>( plugin_search_dirs); } } // namespace plugins } // namespace nuis diff --git a/src/samples/CMakeLists.txt b/src/samples/CMakeLists.txt index 6313da0..c1a9654 100644 --- a/src/samples/CMakeLists.txt +++ b/src/samples/CMakeLists.txt @@ -1,24 +1,25 @@ set(samples_header_files IEventProcessor.hxx IDataComparison.hxx SimpleDataComparison.hxx + MultiDataComparison.hxx SimpleMCStudy.hxx) install(FILES ${samples_header_files} DESTINATION include/samples) add_subdirectory(MCTools) add_subdirectory(nuA) cmessage(DEBUG "INuADataComparisons: ${INuADataComparisons}") SET(INuADataComparisons_List) if(NOT IDataComparisons STREQUAL "") string(REPLACE ";" ", " INuADataComparisons_List "${INuADataComparisons}") endif() cmessage(DEBUG "INuADataComparisons_List: ${INuADataComparisons_List}") SET(INuADataComparisons_List ${INuADataComparisons_List} PARENT_SCOPE) SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE) diff --git a/src/samples/IEventProcessor.hxx b/src/samples/IEventProcessor.hxx index 90eff1a..a3693ce 100644 --- a/src/samples/IEventProcessor.hxx +++ b/src/samples/IEventProcessor.hxx @@ -1,119 +1,130 @@ // Copyright 2018 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 SAMPLES_IEventProcessor_HXX_SEEN #define SAMPLES_IEventProcessor_HXX_SEEN #include "plugins/traits.hxx" #include "exception/exception.hxx" #include "config/GlobalConfiguration.hxx" #include "fhiclcpp/ParameterSet.h" #include "TH1.h" #include #include #include namespace nuis { namespace event { class FullEvent; class MinimalEvent; } // namespace event } // namespace nuis #define IEventProcessor_DEBUG(v) \ if (verb > 2) { \ std::cout << "[DEBUG]: " << v << std::endl; \ } #define IEventProcessor_INFO(v) \ if (verb > 1) { \ std::cout << "[INFO]: " << v << std::endl; \ } #define IEventProcessor_WARN(v) \ if (verb > 0) { \ std::cout << "[WARN] " << __FILENAME__ << ":" << __LINE__ << " : " << v \ << std::endl; \ } class IEventProcessor { protected: NEW_NUIS_EXCEPT(unknown_IEventProcessor_verbosity); enum sample_verbosity { kSilent = 0, kWarn = 1, kReticent = 2, kVerbose = 3 }; sample_verbosity verb; void SetSampleVerbosity(std::string v) { v = nuis::config::GetDocument().get( "global.sample.verbosity_override", v); if (v == "Silent") { verb = kSilent; } else if (v == "Warn") { verb = kWarn; } else if (v == "Reticent") { verb = kReticent; } else if (v == "Verbose") { verb = kVerbose; } else { throw unknown_IEventProcessor_verbosity() << "[ERROR]: Failed to parse IEventProcessor verbosity setting from: " << std::quoted(v); } } public: NEW_NUIS_EXCEPT(uninitialized_IEventProcessor); NEW_NUIS_EXCEPT(unimplemented_IEventProcessor_optional_method); IEventProcessor() { SetSampleVerbosity(nuis::config::GetDocument().get( "global.sample.verbosity_default", "Silent")); TH1::SetDefaultSumw2(); } virtual void Initialize(fhicl::ParameterSet const &) = 0; // Interface for processing a single, external event // // IEventProcessors are not required to implement processing events from // 'outside'. - virtual void ProcessEvent(nuis::event::FullEvent const &) { + virtual void ProcessEvent(nuis::event::FullEvent const &, double weight = 1) { throw unimplemented_IEventProcessor_optional_method() << "[ERROR]: IEventProcessor " << Name() << " does not implement ProcessEvent."; } + // Interface for processing a pre-selected single, external event + // + // IEventProcessors are not required to implement processing events from + // 'outside'. + virtual void ProcessSignalEvent(nuis::event::FullEvent const &, + double weight = 1) { + throw unimplemented_IEventProcessor_optional_method() + << "[ERROR]: IEventProcessor " << Name() + << " does not implement ProcessSignalEvent."; + } + virtual void ProcessSample(size_t nmax = std::numeric_limits::max()) = 0; virtual void Write() = 0; virtual std::string Name() = 0; virtual ~IEventProcessor() {} }; DECLARE_PLUGIN_INTERFACE(IEventProcessor); #endif diff --git a/src/samples/MCTools/CMakeLists.txt b/src/samples/MCTools/CMakeLists.txt index 6d1455b..b325ae9 100644 --- a/src/samples/MCTools/CMakeLists.txt +++ b/src/samples/MCTools/CMakeLists.txt @@ -1,16 +1,16 @@ if(USE_NUWRO) LIST(APPEND MC_TOOL_IMPL NuisToNuWro.cxx) endif(USE_NUWRO) -LIST(APPEND MC_TOOL_IMPL VerboseEventSummary.cxx NuisToCAFTruth) +LIST(APPEND MC_TOOL_IMPL VerboseEventSummary.cxx NuisToCAFTruth.cxx EventSummary_ECTJune2019.cxx) add_library(MCTools SHARED ${MC_TOOL_IMPL}) target_link_libraries(MCTools ${SAMPLES_MCTOOLS_LIBS}) if(USE_NUWRO) target_compile_options(MCTools PUBLIC ${NUWRO_CXX_FLAGS}) include_directories(${NUWRO_INCLUDE_DIRS}) target_link_libraries(MCTools -Wl,--no-as-needed ${NUWRO_IMPORTED_TARGETS}) endif() install(TARGETS MCTools DESTINATION plugins) diff --git a/src/samples/MCTools/EventSummary_ECTJune2019.cxx b/src/samples/MCTools/EventSummary_ECTJune2019.cxx new file mode 100644 index 0000000..d029fc5 --- /dev/null +++ b/src/samples/MCTools/EventSummary_ECTJune2019.cxx @@ -0,0 +1,311 @@ +// Copyright 2018 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 "samples/SimpleMCStudy.hxx" + +#include "utility/EventTopologyUtility.hxx" +#include "utility/FullEventUtility.hxx" +#include "utility/InteractionChannelUtility.hxx" +#include "utility/KinematicUtility.hxx" +#include "utility/PDGCodeUtility.hxx" +#include "utility/ROOTUtility.hxx" + +#include "utility/experimental/MINERvAUtility.hxx" +#include "utility/experimental/T2KUtility.hxx" + +#include "persistency/ROOTOutput.hxx" + +using namespace nuis::event; +using namespace nuis::utility; + +class EventSummary_ECTJune2019 : public SimpleMCStudy { + + TreeFile EvSumTree; + + struct EvSum { + double e_nu; + int pdg_nu; + + double e_fslep; + double costheta_fslep; + int pdg_fslep; + + double e_hmfsproton; + double costheta_hmfsproton; + + double e_hmfspi; + double costheta_hmfspi; + int pdg_hmfspi; + + double dpt; + double dat; + double dphit; + double dpt_mnv; + double dat_mnv; + double dphit_mnv; + double pn_mnv; + double q0; + double q3; + double dpn; + double EAvailProxy; + double Q2; + double W_rec; + + int mode; + + int NFSPip; + int NFSPim; + int NFSPi0; + int NFSProton; + int NFSNeutron; + int NFSGamma; + int NFSKaon; + int NFSOther; + + bool isCC; + + bool isQE; + bool isQEL; + bool is2p2h; + bool isNucleonSPP; + bool isCoh; + bool isMultiPi; + bool isDIS; + bool isOtherChannel; + + bool is0Pi; + bool is1Pi; + bool isNPi; + bool isOtherTopo; + + bool isT2K0Pi; + bool isT2K1Pip_CH; + bool isT2K_STV; + + bool isMINERvA0Pi; + bool isMINERvA1CPi; + bool isMINERvA1Pi0; + bool isMINERvALowRecoil; + bool isMINERvASTV; + + double xsweight; + }; + + EvSum es; + +public: + EventSummary_ECTJune2019() { ReadGlobalConfigDefaults(); } + + void SetBranchAddresses() { + EvSumTree->Branch("e_nu", &es.e_nu, "e_nu/D"); + EvSumTree->Branch("pdg_nu", &es.pdg_nu, "pdg_nu/I"); + + EvSumTree->Branch("e_fslep", &es.e_fslep, "e_fslep/D"); + EvSumTree->Branch("costheta_fslep", &es.costheta_fslep, "costheta_fslep/D"); + EvSumTree->Branch("pdg_fslep", &es.pdg_fslep, "pdg_fslep/I"); + + EvSumTree->Branch("e_hmfsproton", &es.e_hmfsproton, "e_hmfsproton/D"); + EvSumTree->Branch("costheta_hmfsproton", &es.costheta_hmfsproton, + "costheta_hmfsproton/D"); + + EvSumTree->Branch("e_hmfspi", &es.e_hmfspi, "e_hmfspi/D"); + EvSumTree->Branch("costheta_hmfspi", &es.costheta_hmfspi, + "costheta_hmfspi/D"); + EvSumTree->Branch("pdg_hmfspi", &es.pdg_hmfspi, "pdg_hmfspi/I"); + + EvSumTree->Branch("dpt", &es.dpt, "dpt/D"); + EvSumTree->Branch("dat", &es.dat, "dat/D"); + EvSumTree->Branch("dphit", &es.dphit, "dphit/D"); + EvSumTree->Branch("dpt_mnv", &es.dpt_mnv, "dpt_mnv/D"); + EvSumTree->Branch("dat_mnv", &es.dat_mnv, "dat_mnv/D"); + EvSumTree->Branch("dphit_mnv", &es.dphit_mnv, "dphit_mnv/D"); + EvSumTree->Branch("pn_mnv", &es.pn_mnv, "pn_mnv/D"); + EvSumTree->Branch("q0", &es.q0, "q0/D"); + EvSumTree->Branch("q3", &es.q3, "q3/D"); + EvSumTree->Branch("EAvailProxy", &es.EAvailProxy, "EAvailProxy/D"); + EvSumTree->Branch("Q2", &es.Q2, "Q2/D"); + EvSumTree->Branch("W_rec", &es.W_rec, "W_rec/D"); + + EvSumTree->Branch("mode", &es.mode, "mode/I"); + + EvSumTree->Branch("NFSPip", &es.NFSPip, "NFSPip/I"); + EvSumTree->Branch("NFSPim", &es.NFSPim, "NFSPim/I"); + EvSumTree->Branch("NFSPi0", &es.NFSPi0, "NFSPi0/I"); + EvSumTree->Branch("NFSProton", &es.NFSProton, "NFSProton/I"); + EvSumTree->Branch("NFSNeutron", &es.NFSNeutron, "NFSNeutron/I"); + EvSumTree->Branch("NFSGamma", &es.NFSGamma, "NFSGamma/I"); + EvSumTree->Branch("NFSKaon", &es.NFSKaon, "NFSKaon/I"); + EvSumTree->Branch("NFSOther", &es.NFSOther, "NFSOther/I"); + + EvSumTree->Branch("isCC", &es.isCC, "isCC/O"); + + EvSumTree->Branch("isQE", &es.isQE, "isQE/O"); + EvSumTree->Branch("isQEL", &es.isQEL, "isQEL/O"); + EvSumTree->Branch("is2p2h", &es.is2p2h, "is2p2h/O"); + EvSumTree->Branch("isNucleonSPP", &es.isNucleonSPP, "isNucleonSPP/O"); + EvSumTree->Branch("isCoh", &es.isCoh, "isCoh/O"); + EvSumTree->Branch("isMultiPi", &es.isMultiPi, "isMultiPi/O"); + EvSumTree->Branch("isDIS", &es.isDIS, "isDIS/O"); + EvSumTree->Branch("isOtherChannel", &es.isOtherChannel, "isOtherChannel/O"); + + EvSumTree->Branch("is0Pi", &es.is0Pi, "is0Pi/O"); + EvSumTree->Branch("is1Pi", &es.is1Pi, "is1Pi/O"); + EvSumTree->Branch("isNPi", &es.isNPi, "isNPi/O"); + EvSumTree->Branch("isOtherTopo", &es.isOtherTopo, "isOtherTopo/O"); + + EvSumTree->Branch("isT2K0Pi", &es.isT2K0Pi, "isT2K0Pi/O"); + EvSumTree->Branch("isT2K1Pip_CH", &es.isT2K1Pip_CH, "isT2K1Pip_CH/O"); + EvSumTree->Branch("isT2K_STV", &es.isT2K_STV, "isT2K_STV/O"); + + EvSumTree->Branch("isMINERvA0Pi", &es.isMINERvA0Pi, "isMINERvA0Pi/O"); + EvSumTree->Branch("isMINERvA1CPi", &es.isMINERvA1CPi, "isMINERvA1CPi/O"); + EvSumTree->Branch("isMINERvA1Pi0", &es.isMINERvA1Pi0, "isMINERvA1Pi0/O"); + EvSumTree->Branch("isMINERvALowRecoil", &es.isMINERvALowRecoil, + "isMINERvALowRecoil/O"); + EvSumTree->Branch("isMINERvASTV", &es.isMINERvASTV, "isMINERvASTV/O"); + + EvSumTree->Branch("xsweight", &es.xsweight, "xsweight/D"); + } + + void Initialize(fhicl::ParameterSet const &instance_sample_configuration) { + + // Perform any per-sample configuration in the base class + SimpleMCStudy::Initialize(instance_sample_configuration); + + EvSumTree = AddNewTTreeToFile(nuis::persistency::GetOutputFile().get(), + "ECTEvSumTree", write_directory); + SetBranchAddresses(); + + //! Here you do whatever you want to do on a per-event basis. + //! This method will be run for every event given by the input handler, the + //! weight allows differential xsecs to be plotted easily and also contains + //! any requested reweightable parameter variations. + ProcessEventFunction = [&](nuis::event::FullEvent const &ev, + double weight) -> void { + Particle ISNu = GetHMISNeutralLepton(ev); + if (!ISNu) { + return; + } + + Particle FSLep = GetHMFSLepton(ev); + if (!FSLep) { + return; + } + + es.e_nu = ISNu.E(); + es.pdg_nu = ISNu.pdg; + + es.e_fslep = FSLep.E(); + es.costheta_fslep = FSLep.CosTheta(); + es.pdg_fslep = FSLep.pdg; + + Particle FSPi = GetHMFSPion(ev); + if (!FSPi) { + es.e_hmfspi = 0; + es.costheta_hmfspi = 0; + es.pdg_hmfspi = 0; + } else { + es.e_hmfspi = FSPi.E(); + es.costheta_hmfspi = FSPi.CosTheta(); + es.pdg_hmfspi = FSPi.pdg; + } + + Particle FSProton = GetHMFSProton(ev); + if (!FSProton) { + es.e_hmfsproton = 0; + es.costheta_hmfsproton = 0; + } else { + es.e_hmfsproton = FSProton.E(); + es.costheta_hmfsproton = FSProton.CosTheta(); + } + + es.dpt = GetDeltaPT_CC0PiN(ev, ISNu.pdg).Mag(); + es.dat = GetDeltaPhiT_CC0PiN(ev, ISNu.pdg); + es.dphit = GetDeltaAlphaT_CC0PiN(ev, ISNu.pdg); + es.dpt_mnv = mnv::GetDeltaPT_CC0PiN_mnv(ev).Mag(); + es.dat_mnv = mnv::GetDeltaAlphaT_CC0PiN_mnv(ev); + es.dphit_mnv = mnv::GetDeltaPhiT_CC0PiN_mnv(ev); + es.pn_mnv = mnv::GetNeutronMomentumReco_CC0PiN_mnv(ev); + es.EAvailProxy = GetEAvailProxy(ev); + TLorentzVector EMTransfer_lep = GetEnergyMomentumTransfer(ev); + es.Q2 = -EMTransfer_lep.Mag2(); + es.q0 = EMTransfer_lep.E(); + es.q3 = EMTransfer_lep.Vect().Mag(); + es.W_rec = GetNeutrinoWRec(ev); + + es.mode = ChannelToInt(ev.mode); + + es.NFSPip = GetNParticles(ev, {pdgcodes::kPiPlus}); + es.NFSPim = GetNParticles(ev, {pdgcodes::kPiMinus}); + es.NFSPi0 = GetNParticles(ev, {pdgcodes::kPi0}); + es.NFSProton = GetNParticles(ev, {pdgcodes::kProton}); + es.NFSNeutron = GetNParticles(ev, {pdgcodes::kNeutron}); + es.NFSGamma = GetNParticles(ev, {pdgcodes::kGamma}); + es.NFSKaon = GetNParticles(ev, pdgcodes::Kaons); + es.NFSOther = GetNFSOthers(ev) - (es.NFSGamma + es.NFSKaon); + + es.isCC = IsCC(ev.mode); + + es.isQE = IsQE(ev.mode); + es.isQEL = IsQEL(ev.mode); + es.is2p2h = Is2p2h(ev.mode); + es.isNucleonSPP = IsNucleonSPP(ev.mode); + es.isCoh = IsCoh(ev.mode); + es.isMultiPi = IsMultiPi(ev.mode); + es.isDIS = IsDIS(ev.mode); + es.isOtherChannel = + !(es.isQE || es.isQEL || es.is2p2h || es.isNucleonSPP || es.isCoh || + es.isMultiPi || es.isDIS); + + size_t NPi = GetNFSPions(ev); + size_t NOthers = GetNFSOthers(ev); + es.is0Pi = (NPi == 0) && (!NOthers); + es.is1Pi = (NPi == 1) && (!NOthers); + es.isNPi = (NPi > 1) && (!NOthers); + es.isOtherTopo = NOthers; + + es.isT2K0Pi = + t2k::IsCC0Pi_NumProtons(ev).first; // Don't care about n protons + es.isT2K1Pip_CH = t2k::IsCC1Pip_CH_RecPi(ev); + es.isT2K_STV = t2k::IsCC0Pi_STV(ev); + + es.isMINERvA0Pi = + mnv::IsCC0Pi_NumProtons(ev).first; // Don't care about n protons + es.isMINERvA1CPi = mnv::IsCC1CPi_2017(ev); + es.isMINERvA1Pi0 = mnv::IsCC1Pi0_2016(ev); + es.isMINERvALowRecoil = mnv::IsCCIncLowRecoil(ev); + es.isMINERvASTV = mnv::IsCC0PiNp_STV(ev); + + es.xsweight = weight; + + EvSumTree->Fill(); + }; + } + + std::string Name() { return "EventSummary_ECTJune2019"; } + + //! Here you can write any custom histograms to TTrees that your sample has + //! been handling. + void Write() {} +}; + +//! These declarations allow your class to be loaded dynamically by NUISANCE +DECLARE_PLUGIN(IEventProcessor, EventSummary_ECTJune2019); diff --git a/src/samples/MultiDataComparison.hxx b/src/samples/MultiDataComparison.hxx new file mode 100644 index 0000000..cb9472c --- /dev/null +++ b/src/samples/MultiDataComparison.hxx @@ -0,0 +1,231 @@ +// Copyright 2018 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 . + *******************************************************************************/ + +#pragma once + +#include "samples/IDataComparison.hxx" + +#include "event/FullEvent.hxx" + +#include "input/InputManager.hxx" + +#include "persistency/ROOTOutput.hxx" + +#include "utility/FileSystemUtility.hxx" +#include "utility/HistogramUtility.hxx" +#include "utility/KinematicUtility.hxx" +#include "utility/ROOTUtility.hxx" +#include "utility/StatsUtility.hxx" + +#include +#include +#include +#include +#include + +///Wraps multiple projections of the same signal selection +class MultiDataComparison : public IDataComparison { + + NEW_NUIS_EXCEPT(invalid_MultiDataComparison_initialization); + NEW_NUIS_EXCEPT(MultiDataComparison_already_finalized); + +protected: + std::string fName; + + nuis::input::InputManager::Input_id_t fIH_id; + std::string write_directory; + size_t NMaxSample_override; + + std::string fJournalReference; + std::string fDOI; + std::string fYear; + std::string fTargetMaterial; + std::string fFluxDescription; + std::string fSignalDescription; + + std::function IsSigFunc; + + std::vector>> + Comparisons; + + nuis::utility::KinematicRange energy_cut; + +public: + MultiDataComparison(std::string name) { + fName = std::move(name); + fIH_id = std::numeric_limits::max(); + write_directory = ""; + NMaxSample_override = std::numeric_limits::max(); + + fJournalReference = ""; + fDOI = ""; + fYear = ""; + fTargetMaterial = ""; + fFluxDescription = ""; + fSignalDescription = ""; + + energy_cut = nuis::utility::KinematicRange{ + std::numeric_limits::max(), std::numeric_limits::max()}; + } + + std::string Name() { return fName; } + + fhicl::ParameterSet fGlobalConfig; + fhicl::ParameterSet fInstanceConfig; + + void ReadGlobalConfigDefaults() { + fGlobalConfig = nuis::config::GetDocument().get( + std::string("global.sample_configuration.") + Name(), + fhicl::ParameterSet()); + + if (!fJournalReference.length()) { + fJournalReference = fGlobalConfig.get("journal_reference", + "Not yet published"); + } + if (!fDOI.length()) { + fDOI = fGlobalConfig.get("DOI", "Unknown DOI"); + } + if (!fYear.length()) { + fYear = fGlobalConfig.get("year", "Unknown year"); + } + + if (!fTargetMaterial.length()) { + fTargetMaterial = fGlobalConfig.get( + "target_material", "Unknown target material"); + } + if (!fFluxDescription.length()) { + fFluxDescription = + fGlobalConfig.get("flux_description", "Unknown flux"); + } + if (!fSignalDescription.length()) { + fSignalDescription = fGlobalConfig.get("signal_description", + "Unknown Signal"); + } + + if ((energy_cut.first == std::numeric_limits::max()) && + (fGlobalConfig.has_key("enu_range"))) { + energy_cut = fGlobalConfig.get>("enu_range"); + } + } + + virtual std::string GetJournalReference() { return fJournalReference; } + virtual std::string GetDOI() { return fDOI; } + virtual std::string GetYear() { return fYear; } + virtual std::string GetTargetMaterial() { return fTargetMaterial; } + virtual std::string GetFluxDescription() { return fFluxDescription; } + virtual std::string GetSignalDescription() { return fSignalDescription; } + + void Initialize(fhicl::ParameterSet const &instance_sample_configuration) { + + fInstanceConfig = instance_sample_configuration; + + SetSampleVerbosity(fInstanceConfig.get( + "verbosity", nuis::config::GetDocument().get( + "global.sample.verbosity_default", "Reticent"))); + + ReadGlobalConfigDefaults(); + + if (fInstanceConfig.has_key("verbosity")) { + SetSampleVerbosity(fInstanceConfig.get("verbosity")); + } + + NMaxSample_override = + fInstanceConfig.get("nmax", std::numeric_limits::max()); + + write_directory = + fInstanceConfig.get("write_directory", Name()); + } + + void ProcessSample(size_t nmax = std::numeric_limits::max()) { + if (fIH_id == + std::numeric_limits::max()) { + fIH_id = + nuis::input::InputManager::Get().EnsureInputLoaded(fInstanceConfig); + } + IInputHandler const &IH = + nuis::input::InputManager::Get().GetInputHandler(fIH_id); + + nmax = std::min(NMaxSample_override, nmax); + + size_t NEvsToProcess = std::min(nmax, IH.GetNEvents()); + + double nmax_scaling = double(IH.GetNEvents()) / double(NEvsToProcess); + + size_t NToShout = NEvsToProcess / 10; + IEventProcessor_INFO("Sample " << std::quoted(Name()) << " processing " + << NEvsToProcess << " events."); + + IInputHandler::ev_index_t ev_idx = 0; + + size_t NSig = 0; + while (ev_idx < NEvsToProcess) { + bool is_sig = false; + + nuis::event::FullEvent const &fev = IH.GetFullEvent(ev_idx); + + if (NToShout && !(ev_idx % NToShout)) { + IEventProcessor_INFO("\tProcessed " << ev_idx << "/" << NEvsToProcess + << " events."); + } + + if (IsSigFunc(fev)) { + NSig++; + for (auto &comp : Comparisons) { + comp.second->ProcessSignalEvent(fev, IH.GetEventWeight(ev_idx) * + nmax_scaling); + } + } + + ev_idx++; + } + IEventProcessor_INFO("\t" << NSig << "/" << NEvsToProcess + << " events passed selection."); + } + void Write() {} + + double GetGOF() { + double totGOF = 0; + for (auto const &comp : Comparisons) { + double sGOF = comp.second->GetGOF(); + if (sGOF != std::numeric_limits::max()) { + totGOF += sGOF; + } else { + std::cout << "[WARN]: Projection \"" << comp.first << "\" of sample " + << Name() << " produced bad GOF." << std::endl; + } + } + return totGOF; + } + double GetNDOGuess() { + size_t totNDOGuess = 0; + for (auto const &comp : Comparisons) { + totNDOGuess += comp.second->GetNDOGuess(); + } + return totNDOGuess; + } + + fhicl::ParameterSet GetExampleConfiguration() { + fhicl::ParameterSet exps; + exps.put("name", Name()); + exps.put("input_type", "NuWro/NEUT/GENIE"); + exps.put("file", "input_events.root"); + exps.put("write_directory", Name()); + return exps; + } +}; diff --git a/src/samples/SimpleDataComparison.hxx b/src/samples/SimpleDataComparison.hxx index b8f91e7..fbfb548 100644 --- a/src/samples/SimpleDataComparison.hxx +++ b/src/samples/SimpleDataComparison.hxx @@ -1,451 +1,493 @@ // Copyright 2018 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 . *******************************************************************************/ #pragma once #include "samples/IDataComparison.hxx" #include "event/FullEvent.hxx" #include "input/InputManager.hxx" #include "persistency/ROOTOutput.hxx" #include "utility/FileSystemUtility.hxx" #include "utility/HistogramUtility.hxx" #include "utility/KinematicUtility.hxx" #include "utility/ROOTUtility.hxx" #include "utility/StatsUtility.hxx" #include #include #include #include #include template ::type> class SimpleDataComparison : public IDataComparison { NEW_NUIS_EXCEPT(invalid_SimpleDataComparison_initialization); NEW_NUIS_EXCEPT(SimpleDataComparison_already_finalized); +public: + static size_t const NDim = nd; + std::string write_directory; + protected: using HistType = HT; using TH_Help = typename nuis::utility::TH_Helper; + nuis::input::InputManager::Input_id_t fIH_id; - static size_t const NDim = nd; + std::string fName; - nuis::input::InputManager::Input_id_t fIH_id; - std::string write_directory; size_t NMaxSample_override; int fIsShapeOnly; int fIsFluxUnfolded; std::vector fSignalCache; std::vector> fProjectionCache; bool fUseCache; + bool fModeHists; std::string fDataInputDescriptor; std::unique_ptr fData; std::string fMaskInputDescriptor; std::unique_ptr fMask; std::string fCovarianceInputDescriptor; std::unique_ptr fCovariance; std::unique_ptr fPrediction; + std::map> fPrediction_modes; std::unique_ptr fPrediction_xsec; std::unique_ptr fPrediction_shape; std::unique_ptr fPrediction_comparison; bool fComparisonFinalized; std::string fJournalReference; std::string fDOI; std::string fYear; std::string fTargetMaterial; std::string fFluxDescription; std::string fSignalDescription; - nuis::utility::ENuRange energy_cut; + nuis::utility::KinematicRange energy_cut; std::function IsSigFunc; - std::function(nuis::event::FullEvent const &)> - CompProjFunc; // If assigned by subclass will be called on for all events, bool signifies // whether the event was selected. std::function ProcessExtraFunc; public: - SimpleDataComparison() { + std::function(nuis::event::FullEvent const &)> + CompProjFunc; + + SimpleDataComparison(std::string name) { + fName = std::move(name); fIH_id = std::numeric_limits::max(); fUseCache = false; write_directory = ""; NMaxSample_override = std::numeric_limits::max(); fDataInputDescriptor = ""; fData = nullptr; fMaskInputDescriptor = ""; fMask = nullptr; fCovarianceInputDescriptor = ""; fCovariance = nullptr; fPrediction = nullptr; fPrediction_xsec = nullptr; fPrediction_shape = nullptr; fPrediction_comparison = nullptr; fComparisonFinalized = false; IsSigFunc = [](nuis::event::FullEvent const &) -> bool { return true; }; CompProjFunc = [](nuis::event::FullEvent const &) -> std::array { std::array arr; for (NumericT &el : arr) { el = 0; } return arr; }; ProcessExtraFunc = std::function(); fJournalReference = ""; fDOI = ""; fYear = ""; fTargetMaterial = ""; fFluxDescription = ""; fSignalDescription = ""; fIsShapeOnly = -1; fIsFluxUnfolded = -1; - energy_cut = nuis::utility::ENuRange{std::numeric_limits::max(), - std::numeric_limits::max()}; + energy_cut = nuis::utility::KinematicRange{ + std::numeric_limits::max(), std::numeric_limits::max()}; } + std::string Name() { return fName; } + void SetUseCache(bool uc = true) { fUseCache = uc; if (!uc) { fSignalCache.resize(0); fProjectionCache.resize(0); } } fhicl::ParameterSet fGlobalConfig; fhicl::ParameterSet fInstanceConfig; void ReadGlobalConfigDefaults() { fGlobalConfig = nuis::config::GetDocument().get( std::string("global.sample_configuration.") + Name(), fhicl::ParameterSet()); if (!fJournalReference.length()) { fJournalReference = fGlobalConfig.get("journal_reference", "Not yet published"); } if (!fDOI.length()) { fDOI = fGlobalConfig.get("DOI", "Unknown DOI"); } if (!fYear.length()) { fYear = fGlobalConfig.get("year", "Unknown year"); } if (!fTargetMaterial.length()) { fTargetMaterial = fGlobalConfig.get( "target_material", "Unknown target material"); } if (!fFluxDescription.length()) { fFluxDescription = fGlobalConfig.get("flux_description", "Unknown flux"); } if (!fSignalDescription.length()) { fSignalDescription = fGlobalConfig.get("signal_description", "Unknown Signal"); } if (fIsShapeOnly == -1) { fIsShapeOnly = fGlobalConfig.get("shape_only", false); } if (fIsFluxUnfolded == -1) { fIsFluxUnfolded = fGlobalConfig.get("flux_unfolded", false); } if ((energy_cut.first == std::numeric_limits::max()) && (fGlobalConfig.has_key("enu_range"))) { energy_cut = fGlobalConfig.get>("enu_range"); } } virtual std::string GetJournalReference() { return fJournalReference; } virtual std::string GetDOI() { return fDOI; } virtual std::string GetYear() { return fYear; } virtual std::string GetTargetMaterial() { return fTargetMaterial; } virtual std::string GetFluxDescription() { return fFluxDescription; } virtual std::string GetSignalDescription() { return fSignalDescription; } void SetShapeOnly(bool iso) { fIsShapeOnly = iso; } void SetFluxUnfolded(bool ifo) { fIsFluxUnfolded = ifo; } void SetData(std::string const &data_descriptor) { fDataInputDescriptor = data_descriptor; } void SetMask(std::string const &mask_descriptor) { fMaskInputDescriptor = mask_descriptor; } void SetCovariance(std::string const &cov_descriptor) { fCovarianceInputDescriptor = cov_descriptor; } virtual void FillProjection(std::array const &proj, NumericT event_weight) { TH_Help::Fill(fPrediction, proj, event_weight); } + virtual void FillModeProjection(std::array const &proj, + NumericT event_weight, + nuis::event::Channel_t mode) { + if (!fPrediction_modes.count(mode)) { + std::stringstream ss; + ss << "Prediction_" << mode; + fPrediction_modes[mode] = + nuis::utility::Clone(fPrediction, true, ss.str()); + } + TH_Help::Fill(fPrediction_modes[mode], proj, event_weight); + } + + virtual void ProcessSignalEvent(nuis::event::FullEvent const &fev, + double weight = 1) { + auto const &proj = CompProjFunc(fev); + FillProjection(proj, weight); + if (fModeHists) { + FillModeProjection(proj, weight, fev.mode); + } + } virtual void FinalizeComparison() { if (fComparisonFinalized) { throw SimpleDataComparison_already_finalized() << "[ERROR]: Attempted to re-finalize a comparison for " "SimpleDataComparison: " << std::quoted(Name()); } fPrediction_xsec = nuis::utility::Clone(fPrediction, false, "Prediction_xsec"); - IInputHandler const &IH = - nuis::input::InputManager::Get().GetInputHandler(fIH_id); - TH_Help::Scale(fPrediction_xsec, 1.0, "width"); - + if (fModeHists) { + for (auto &mh : fPrediction_modes) { + TH_Help::Scale(mh.second, 1.0, "width"); + } + } // If we have a flux cut if (energy_cut.first != std::numeric_limits::max()) { + IInputHandler const &IH = + nuis::input::InputManager::Get().GetInputHandler(fIH_id); TH_Help::Scale(fPrediction_xsec, IH.GetXSecScaleFactor(energy_cut)); } fPrediction_shape = nuis::utility::Clone(fPrediction_xsec, false, "Prediction_shape"); if (fData) { TH_Help::Scale(fPrediction_shape, TH_Help::Integral(fData, "width") / TH_Help::Integral(fPrediction_shape, "width")); } else { IEventProcessor_WARN( "When Finalizing comparison, no Data histogram available."); } if (fIsFluxUnfolded) { // fPrediction_comparison } else if (fIsShapeOnly) { fPrediction_comparison = nuis::utility::Clone(fPrediction_shape, false, "Prediction_comparison"); } else { fPrediction_comparison = nuis::utility::Clone(fPrediction_xsec, false, "Prediction_comparison"); } fComparisonFinalized = true; } void Initialize(fhicl::ParameterSet const &instance_sample_configuration) { fInstanceConfig = instance_sample_configuration; SetSampleVerbosity(fInstanceConfig.get( "verbosity", nuis::config::GetDocument().get( "global.sample.verbosity_default", "Reticent"))); ReadGlobalConfigDefaults(); if (fInstanceConfig.has_key("fake_data")) { fData = nuis::utility::GetHistogram( fInstanceConfig.get("fake_data_histogram")); } else if (!fGlobalConfig.get("has_data", true) || !fInstanceConfig.get("has_data", true)) { // Explicitly not expecting data } else { if (!fDataInputDescriptor.length()) { if (!fGlobalConfig.has_key("data_descriptor")) { throw invalid_SimpleDataComparison_initialization() << "[ERROR]: SimpleDataComparison::Initialize for " "IDataComparison: " << std::quoted(Name()) << " failed as no input data was set by a call to " "SimpleDataComparison::SetData and no data_descriptor for " "this SimpleDataComparison could be found in the global " "configuration."; } fDataInputDescriptor = fGlobalConfig.get("data_descriptor"); } fData = nuis::utility::GetHistogram(fDataInputDescriptor); } if (!fPrediction) { if (!fData) { throw invalid_SimpleDataComparison_initialization() << "[ERROR]: SimpleDataComparison::Initialize for " "IDataComparison: " << std::quoted(Name()) << " failed. As `has_data: false` was set in the configuration " "(global: " << (!fGlobalConfig.get("has_data", true)) << ", instance: " << !fInstanceConfig.get("has_data", true) << "), the instance constructor must supply the fPrediction " "binning, and it wasn't."; } fPrediction = nuis::utility::Clone(fData, true, "Prediction"); } if (fCovarianceInputDescriptor.length()) { fCovariance = nuis::utility::GetHistogram(fCovarianceInputDescriptor); } else if (fGlobalConfig.has_key("covariance_descriptor")) { fCovariance = nuis::utility::GetHistogram( fGlobalConfig.get("covariance_descriptor")); } if (fMaskInputDescriptor.length()) { fMask = nuis::utility::GetHistogram(fMaskInputDescriptor); } else if (fGlobalConfig.has_key("mask_descriptor")) { fMask = nuis::utility::GetHistogram( fGlobalConfig.get("mask_descriptor")); } if (fInstanceConfig.has_key("verbosity")) { SetSampleVerbosity(fInstanceConfig.get("verbosity")); } + fModeHists = fInstanceConfig.get("write_mode_hists", false); + NMaxSample_override = fInstanceConfig.get("nmax", std::numeric_limits::max()); write_directory = fInstanceConfig.get("write_directory", Name()); - - fIH_id = - nuis::input::InputManager::Get().EnsureInputLoaded(fInstanceConfig); } void ProcessSample(size_t nmax = std::numeric_limits::max()) { if (fIH_id == std::numeric_limits::max()) { - throw uninitialized_IEventProcessor(); + fIH_id = + nuis::input::InputManager::Get().EnsureInputLoaded(fInstanceConfig); } IInputHandler const &IH = nuis::input::InputManager::Get().GetInputHandler(fIH_id); nmax = std::min(NMaxSample_override, nmax); size_t NEvsToProcess = std::min(nmax, IH.GetNEvents()); double nmax_scaling = double(IH.GetNEvents()) / double(NEvsToProcess); size_t NToShout = NEvsToProcess / 10; IEventProcessor_INFO("Sample " << std::quoted(Name()) << " processing " << NEvsToProcess << " events."); IInputHandler::ev_index_t ev_idx = 0; bool DetermineSignalEvents = !fSignalCache.size() || !fUseCache; nuis::utility::Clear(*fPrediction); fComparisonFinalized = false; size_t cache_ctr = 0; while (ev_idx < NEvsToProcess) { bool is_sig = false; std::array proj; if (DetermineSignalEvents) { nuis::event::FullEvent const &fev = IH.GetFullEvent(ev_idx); is_sig = IsSigFunc(fev); fSignalCache.push_back(is_sig); if (is_sig) { proj = CompProjFunc(fev); fProjectionCache.push_back(proj); } if (ProcessExtraFunc) { ProcessExtraFunc(fev, is_sig, IH.GetEventWeight(ev_idx) * nmax_scaling); } + if (fModeHists && is_sig) { + FillModeProjection(proj, IH.GetEventWeight(ev_idx) * nmax_scaling, + fev.mode); + } } else { is_sig = fSignalCache[ev_idx]; proj = fProjectionCache[cache_ctr++]; } if (NToShout && !(ev_idx % NToShout)) { IEventProcessor_INFO("\tProcessed " << ev_idx << "/" << NEvsToProcess << " events."); } if (is_sig) { FillProjection(proj, IH.GetEventWeight(ev_idx) * nmax_scaling); } ev_idx++; } IEventProcessor_INFO("\t" << fProjectionCache.size() << "/" << NEvsToProcess << " events passed selection."); } void Write() { if (!fComparisonFinalized) { FinalizeComparison(); } nuis::persistency::WriteToOutputFile( fPrediction_comparison, "Prediction", write_directory); + if (fModeHists) { + for (auto &mh : fPrediction_modes) { + nuis::persistency::WriteToOutputFile( + mh.second, mh.second->GetName(), write_directory); + } + } nuis::persistency::WriteToOutputFile( fPrediction_xsec, "Prediction_xsec", write_directory); if (fData) { nuis::persistency::WriteToOutputFile(fData, "Data", write_directory); nuis::persistency::WriteToOutputFile( fPrediction_shape, "Prediction_shape", write_directory); } } double GetGOF() { if (!fComparisonFinalized) { FinalizeComparison(); } if (fData && fPrediction_comparison) { return nuis::utility::GetChi2(fData, fPrediction_comparison); } else return std::numeric_limits::max(); } double GetNDOGuess() { if (fData) { return TH_Help::Nbins(fData); } return 0; } fhicl::ParameterSet GetExampleConfiguration() { fhicl::ParameterSet exps; exps.put("name", Name()); exps.put("input_type", "NuWro/NEUT/GENIE"); exps.put("file", "input_events.root"); exps.put("write_directory", Name()); return exps; } }; typedef SimpleDataComparison<1, double, TH1D> SimpleDataComparison_1D; typedef SimpleDataComparison<2, double, TH2D> SimpleDataComparison_2D; typedef SimpleDataComparison<1, float, TH1F> SimpleDataComparison_1F; typedef SimpleDataComparison<2, float, TH2F> SimpleDataComparison_2F; typedef SimpleDataComparison<2, double, TH2Poly> SimpleDataComparison_2DPoly; diff --git a/src/samples/nuA/BubbleChamber/ANL/ANL_CCQE_Evt_1DQ2_nu.cxx b/src/samples/nuA/BubbleChamber/ANL/ANL_CCQE_Evt_1DQ2_nu.cxx index 53e7106..86d1a6d 100644 --- a/src/samples/nuA/BubbleChamber/ANL/ANL_CCQE_Evt_1DQ2_nu.cxx +++ b/src/samples/nuA/BubbleChamber/ANL/ANL_CCQE_Evt_1DQ2_nu.cxx @@ -1,190 +1,189 @@ // Copyright 2018 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 "samples/SimpleDataComparison.hxx" #include "utility/FullEventUtility.hxx" #include "utility/KinematicUtility.hxx" #include "utility/PDGCodeUtility.hxx" using namespace nuis::event; using namespace nuis::utility; class ANL_CCQE_Evt_1DQ2_nu : public SimpleDataComparison_1D { public: NEW_NUIS_EXCEPT(invalid_publication_specifier); enum Publication { kPRL31, kPRD16, kPRD26 }; Publication Pub; std::string Pub_str; bool UseD2Corr; std::unique_ptr fD2CorrHist; std::unique_ptr fPrediction_Uncorr; ANL_CCQE_Evt_1DQ2_nu() - : Pub(kPRD26), Pub_str(""), UseD2Corr(false), fD2CorrHist(nullptr) { + : SimpleDataComparison_1D("ANL_CCQE_Evt_1DQ2_nu"), Pub(kPRD26), + Pub_str(""), UseD2Corr(false), fD2CorrHist(nullptr) { ReadGlobalConfigDefaults(); } std::string GetDocumentation() { return "Can specify \"publication: \", where is one of [ PRL31, " "PRD16, PRD26 ] to clarify a publication for comparison. Defaults " "to PRD26.\n" "Can enable deuterium Q2 correction by specifying " "\"use_D2_correction: true\""; } fhicl::ParameterSet GetExampleConfiguration() { fhicl::ParameterSet exps = SimpleDataComparison_1D::GetExampleConfiguration(); exps.put("publication", "PRD26"); exps.put("use_D2_correction", false); return exps; } void Initialize(fhicl::ParameterSet const &instance_sample_configuration) { std::string publication = instance_sample_configuration.get("publication", "PRD26"); if (publication == "PRL31") { Pub = kPRL31; } else if (publication == "PRD16") { Pub = kPRD16; } else if (publication == "PRD26") { Pub = kPRD26; } else { throw invalid_publication_specifier() << "[ERROR]: Found unexpected publication specifier " << std::quoted(publication) << ". Expected one of [ PRL31, PRD16, PRD26 ]"; } switch (Pub) { case kPRL31: { Pub_str = "PRL31_844"; energy_cut = std::pair{0, 3E3}; IEventProcessor_INFO( "Sample " << Name() << " specialized for publication: " << Pub_str); break; } case kPRD16: { Pub_str = "PRD16_3103"; energy_cut = std::pair{0, 6E3}; IEventProcessor_INFO( "Sample " << Name() << " specialized for publication: " << Pub_str); break; } case kPRD26: { Pub_str = "PRD26_537"; energy_cut = std::pair{0, 6E3}; IEventProcessor_INFO( "Sample " << Name() << " specialized for publication: " << Pub_str); break; } } fhicl::ParameterSet const &global_sample_configuration = nuis::config::GetDocument().get( std::string("global.sample_configuration.") + Name(), fhicl::ParameterSet()); SetData(GetDataDir() + "nuA/BubbleChamber/ANL/CCQE/ANL_CCQE_Data_" + Pub_str + ".root;ANL_1DQ2_Data"); SimpleDataComparison_1D::Initialize(instance_sample_configuration); UseD2Corr = instance_sample_configuration.get( "use_D2_correction", global_sample_configuration.get("use_D2_correction", false)); if (UseD2Corr) { fD2CorrHist = nuis::utility::GetHistogram( GetDataDir() + "nuA/BubbleChamber/ANL/CCQE/" "ANL_CCQE_Data_PRL31_844.root;ANL_1DQ2_Correction"); fPrediction_Uncorr = Clone(fPrediction, true); } // Signal selection function IsSigFunc = [&](FullEvent const &fev) -> bool { if (fev.mode != Channel_t::kCCQE) { return false; } Particle ISNumu = GetHMISNeutralLepton(fev); if (!ISNumu) { return false; } if (ISNumu.pdg != pdgcodes::kNuMu) { return false; } if (!energy_cut.IsInRange(ISNumu.P4.E())) { return false; } double Q2 = GetNeutrinoQ2QERec(fev, 0); if (Q2 <= 0) { return false; } return true; }; // 1D Projection function CompProjFunc = [](FullEvent const &fev) -> std::array { return {GetNeutrinoQ2QERec(fev, 0)}; }; } // Used to apply D2 correction if requested virtual void FillProjection(std::array const &proj, double event_weight) { if (UseD2Corr) { TH_Help::Fill(fPrediction_Uncorr, proj, event_weight); event_weight *= fD2CorrHist->Interpolate(proj[0]); } TH_Help::Fill(fPrediction, proj, event_weight); } void FinalizeComparison() { SimpleDataComparison_1D::FinalizeComparison(); if (UseD2Corr) { fPrediction_Uncorr->Scale(1.0, "width"); } } void Write() { SimpleDataComparison_1D::Write(); if (UseD2Corr) { nuis::persistency::WriteToOutputFile( fPrediction_Uncorr, "Prediction_Uncorr", write_directory); } } - - std::string Name() { return "ANL_CCQE_Evt_1DQ2_nu"; } }; DECLARE_PLUGIN(IDataComparison, ANL_CCQE_Evt_1DQ2_nu); DECLARE_PLUGIN(IEventProcessor, ANL_CCQE_Evt_1DQ2_nu); diff --git a/src/samples/nuA/Nuclear/CMakeLists.txt b/src/samples/nuA/Nuclear/CMakeLists.txt index 39d211c..e8b1c7e 100644 --- a/src/samples/nuA/Nuclear/CMakeLists.txt +++ b/src/samples/nuA/Nuclear/CMakeLists.txt @@ -1,4 +1,5 @@ add_subdirectory(T2K) +add_subdirectory(MINERvA) SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE) diff --git a/src/samples/nuA/Nuclear/MINERvA/CC0Pi/CMakeLists.txt b/src/samples/nuA/Nuclear/MINERvA/CC0Pi/CMakeLists.txt new file mode 100644 index 0000000..44a69c1 --- /dev/null +++ b/src/samples/nuA/Nuclear/MINERvA/CC0Pi/CMakeLists.txt @@ -0,0 +1,17 @@ +SET(SAMPLES + MINERvA_CC0PiNProt_CH_xsec_STV_nu) + +foreach(S ${SAMPLES}) + LIST(APPEND INuADataComparisons ${S}) + LIST(APPEND SAMPLES_src ${S}.cxx) +endforeach() + +add_library(MINERvACC0PiDataComparisons SHARED ${SAMPLES_src}) +target_link_libraries(MINERvACC0PiDataComparisons + nuis_event nuis_input nuis_config nuis_persistency) + +install(TARGETS MINERvACC0PiDataComparisons DESTINATION plugins) + +install(FILES MINERvA.CC0Pi.sample.config.fcl DESTINATION fcl) + +SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) diff --git a/src/samples/nuA/Nuclear/MINERvA/CC0Pi/MINERvA.CC0Pi.sample.config.fcl b/src/samples/nuA/Nuclear/MINERvA/CC0Pi/MINERvA.CC0Pi.sample.config.fcl new file mode 100644 index 0000000..25697d8 --- /dev/null +++ b/src/samples/nuA/Nuclear/MINERvA/CC0Pi/MINERvA.CC0Pi.sample.config.fcl @@ -0,0 +1,9 @@ +global.sample_configuration.MINERvA_CC0PiNProt_CH_xsec_STV_nu: { + journal_reference: "" + doi: "" + target_material: "CH" + flux_description: "MINERvA Medium Energy FHC muon neutrino" + signal_description: "CC0Pi + N protons" + shape_only: false + flux_unfolded: false +} diff --git a/src/samples/nuA/Nuclear/MINERvA/CC0Pi/MINERvA_CC0PiNProt_CH_xsec_STV_nu.cxx b/src/samples/nuA/Nuclear/MINERvA/CC0Pi/MINERvA_CC0PiNProt_CH_xsec_STV_nu.cxx new file mode 100644 index 0000000..c47538c --- /dev/null +++ b/src/samples/nuA/Nuclear/MINERvA/CC0Pi/MINERvA_CC0PiNProt_CH_xsec_STV_nu.cxx @@ -0,0 +1,190 @@ +// Copyright 2018 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 "samples/MultiDataComparison.hxx" +#include "samples/SimpleDataComparison.hxx" + +#include "utility/EventTopologyUtility.hxx" +#include "utility/FullEventUtility.hxx" +#include "utility/KinematicUtility.hxx" +#include "utility/PDGCodeUtility.hxx" +#include "utility/experimental/MINERvAUtility.hxx" + +using namespace nuis::event; +using namespace nuis::utility; + +class MINERvA_CC0PiNProt_CH_xsec_STV_nu : public MultiDataComparison { + +public: + MINERvA_CC0PiNProt_CH_xsec_STV_nu() + : MultiDataComparison("MINERvA_CC0PiNProt_CH_xsec_STV_nu") { + ReadGlobalConfigDefaults(); + } + + void Initialize(fhicl::ParameterSet const &instance_sample_configuration) { + + // Perform any per-sample configuration in the base class + MultiDataComparison::Initialize(instance_sample_configuration); + + std::unique_ptr DPT( + new SimpleDataComparison_1D("MINERvA_CC0PiNProt_CH_xsec_STV_nu:dpt")); + DPT->SetData(GetDataDir() + + "nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/" + "MINERvA_1805.05486.root;dpt"); + DPT->Initialize( + instance_sample_configuration.get("dpt", {})); + DPT->CompProjFunc = [](FullEvent const &fev) + -> std::array { + return {mnv::GetDeltaPT_CC0PiN_mnv(fev).Mag() * 1E-3}; + }; + + std::unique_ptr DAT(new SimpleDataComparison_1D( + "MINERvA_CC0PiNProt_CH_xsec_STV_nu:dalphat")); + DAT->SetData(GetDataDir() + + "nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/" + "MINERvA_1805.05486.root;dalphat"); + DAT->Initialize( + instance_sample_configuration.get("dat", {})); + DAT->CompProjFunc = [](FullEvent const &fev) + -> std::array { + return {mnv::GetDeltaAlphaT_CC0PiN_mnv(fev)}; + }; + + std::unique_ptr DPhiT( + new SimpleDataComparison_1D("MINERvA_CC0PiNProt_CH_xsec_STV_nu:dphit")); + DPhiT->SetData(GetDataDir() + + "nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/" + "MINERvA_1805.05486.root;dphit"); + DPhiT->Initialize( + instance_sample_configuration.get("dphit", {})); + DPhiT->CompProjFunc = [](FullEvent const &fev) + -> std::array { + return {mnv::GetDeltaPhiT_CC0PiN_mnv(fev)}; + }; + + std::unique_ptr MuMom(new SimpleDataComparison_1D( + "MINERvA_CC0PiNProt_CH_xsec_STV_nu:muonmomentum")); + MuMom->SetData(GetDataDir() + + "nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/" + "MINERvA_1805.05486.root;muonmomentum"); + MuMom->Initialize(instance_sample_configuration.get( + "muonmomentum", {})); + MuMom->CompProjFunc = [](FullEvent const &fev) + -> std::array { + return {GetHMFSChargedLepton(fev).P() * 1E-3}; + }; + std::unique_ptr MuTheta( + new SimpleDataComparison_1D( + "MINERvA_CC0PiNProt_CH_xsec_STV_nu:muontheta")); + MuTheta->SetData(GetDataDir() + + "nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/" + "MINERvA_1805.05486.root;muontheta"); + MuTheta->Initialize(instance_sample_configuration.get( + "muontheta", {})); + MuTheta->CompProjFunc = [](FullEvent const &fev) + -> std::array { + return {GetHMFSChargedLepton(fev).Theta_deg()}; + }; + + std::unique_ptr ProtonMom( + new SimpleDataComparison_1D( + "MINERvA_CC0PiNProt_CH_xsec_STV_nu:protonmomentum")); + ProtonMom->SetData(GetDataDir() + + "nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/" + "MINERvA_1805.05486.root;protonmomentum"); + ProtonMom->Initialize( + instance_sample_configuration.get("protonmomentum", + {})); + ProtonMom->CompProjFunc = [](FullEvent const &fev) + -> std::array { + return {GetHMFSProtonInPhaseSpace(fev, mnv::CC0PiNProt_ProtonPS).P() * + 1E-3}; + }; + std::unique_ptr ProtonTheta( + new SimpleDataComparison_1D( + "MINERvA_CC0PiNProt_CH_xsec_STV_nu:protontheta")); + ProtonTheta->SetData(GetDataDir() + + "nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/" + "MINERvA_1805.05486.root;protontheta"); + ProtonTheta->Initialize( + instance_sample_configuration.get("protontheta", + {})); + ProtonTheta->CompProjFunc = [](FullEvent const &fev) + -> std::array { + return { + GetHMFSProtonInPhaseSpace(fev, mnv::CC0PiNProt_ProtonPS).Theta_deg()}; + }; + + std::unique_ptr NeutronMomReco( + new SimpleDataComparison_1D( + "MINERvA_CC0PiNProt_CH_xsec_STV_nu:neutronmomentum")); + NeutronMomReco->SetData(GetDataDir() + + "nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/" + "MINERvA_1805.05486.root;neutronmomentum"); + NeutronMomReco->Initialize( + instance_sample_configuration.get( + "neutronmomentum", {})); + NeutronMomReco->CompProjFunc = [](FullEvent const &fev) + -> std::array { + return {mnv::GetNeutronMomentumReco_CC0PiN_mnv(fev)*1E-3}; + }; + + std::string wd_stub = write_directory.size() ? write_directory + "/" : ""; + DPT->write_directory = wd_stub + "dpt"; + DAT->write_directory = wd_stub + "dat"; + DPhiT->write_directory = wd_stub + "dphit"; + MuMom->write_directory = wd_stub + "muonmomentum"; + MuTheta->write_directory = wd_stub + "muontheta"; + ProtonMom->write_directory = wd_stub + "protonmomentum"; + ProtonTheta->write_directory = wd_stub + "protontheta"; + NeutronMomReco->write_directory = wd_stub + "neutronmomentum"; + + Comparisons.emplace_back("dpt", std::move(DPT)); + Comparisons.emplace_back("dat", std::move(DAT)); + Comparisons.emplace_back("dphit", std::move(DPhiT)); + Comparisons.emplace_back("muonmomentum", std::move(MuMom)); + Comparisons.emplace_back("muontheta", std::move(MuTheta)); + Comparisons.emplace_back("protonmomentum", std::move(ProtonMom)); + Comparisons.emplace_back("protontheta", std::move(ProtonTheta)); + Comparisons.emplace_back("neutronmomentum", std::move(NeutronMomReco)); + + //! Define your event signal here. + IsSigFunc = [](FullEvent const &fev) -> bool { + // Is numuCC + if (GetNParticles(fev, {pdgcodes::kMu}) != 1) { + return false; + } + return mnv::IsCC0PiNp_STV(fev); + }; + } + + //! Here you can write any custom histograms to TTrees that your sample has + //! been handling. + void Write() { + for (auto &comp : Comparisons) { + comp.second->Write(); + } + } +}; + +//! These declarations allow your class to be loaded dynamically by NUISANCE +DECLARE_PLUGIN(IDataComparison, MINERvA_CC0PiNProt_CH_xsec_STV_nu); +DECLARE_PLUGIN(IEventProcessor, MINERvA_CC0PiNProt_CH_xsec_STV_nu); diff --git a/src/samples/nuA/Nuclear/T2K/CMakeLists.txt b/src/samples/nuA/Nuclear/MINERvA/CMakeLists.txt similarity index 56% copy from src/samples/nuA/Nuclear/T2K/CMakeLists.txt copy to src/samples/nuA/Nuclear/MINERvA/CMakeLists.txt index b667c78..c611e17 100644 --- a/src/samples/nuA/Nuclear/T2K/CMakeLists.txt +++ b/src/samples/nuA/Nuclear/MINERvA/CMakeLists.txt @@ -1,7 +1,7 @@ add_subdirectory(CC0Pi) -LIST(APPEND INuADataComparisons_FHiCL T2K.sample.config.fcl) -install(FILES T2K.sample.config.fcl DESTINATION fcl) +LIST(APPEND INuADataComparisons_FHiCL MINERvA.sample.config.fcl) +install(FILES MINERvA.sample.config.fcl DESTINATION fcl) SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE) diff --git a/src/samples/nuA/Nuclear/MINERvA/MINERvA.sample.config.fcl b/src/samples/nuA/Nuclear/MINERvA/MINERvA.sample.config.fcl new file mode 100644 index 0000000..09200c4 --- /dev/null +++ b/src/samples/nuA/Nuclear/MINERvA/MINERvA.sample.config.fcl @@ -0,0 +1 @@ +#include "MINERvA.CC0Pi.sample.config.fcl" diff --git a/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0PiNProt_CH_xsec_STV_nu.cxx b/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0PiNProt_CH_xsec_STV_nu.cxx index 8f7ebd6..968bbd9 100644 --- a/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0PiNProt_CH_xsec_STV_nu.cxx +++ b/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0PiNProt_CH_xsec_STV_nu.cxx @@ -1,190 +1,113 @@ // Copyright 2018 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 "samples/MultiDataComparison.hxx" #include "samples/SimpleDataComparison.hxx" #include "utility/EventTopologyUtility.hxx" #include "utility/FullEventUtility.hxx" #include "utility/KinematicUtility.hxx" #include "utility/PDGCodeUtility.hxx" +#include "utility/experimental/T2KUtility.hxx" using namespace nuis::event; using namespace nuis::utility; -class T2K_CC0PiNProt_CH_xsec_STV_nu : public SimpleDataComparison_1D { - - enum STVProjection { kDeltaPT, kDeltaPhiT, kDeltaAlphaT }; +class T2K_CC0PiNProt_CH_xsec_STV_nu : public MultiDataComparison { public: - T2K_CC0PiNProt_CH_xsec_STV_nu() { ReadGlobalConfigDefaults(); } - - std::string GetDocumentation() { - return "Can specify \"projection: \", where is one of [ " - "DeltaPT, DeltaPhiT, DeltaAlphaT ] to clarify a projection for " - "comparison. Defaults to DeltaPT.\n"; - } - - fhicl::ParameterSet GetExampleConfiguration() { - fhicl::ParameterSet exps = - SimpleDataComparison_1D::GetExampleConfiguration(); - - exps.put("projection", "DeltaPT"); - - return exps; - } - - NEW_NUIS_EXCEPT(invalid_projection_specifier); - - STVProjection GetProjection(std::string const &pstring) { - if ((pstring == "DeltaPT") || (pstring == "deltapt") || - (pstring == "dpt")) { - return kDeltaPT; - } else if ((pstring == "DeltaPhiT") || (pstring == "deltaphit") || - (pstring == "dphit")) { - return kDeltaPhiT; - } else if ((pstring == "DeltaAlphaT") || (pstring == "deltaalphat") || - (pstring == "dat")) { - return kDeltaAlphaT; - } - throw invalid_projection_specifier() - << "Unknown projection " << std::quoted(pstring) - << " for comparison: " << std::quoted(Name()); + T2K_CC0PiNProt_CH_xsec_STV_nu() + : MultiDataComparison("T2K_CC0PiNProt_CH_xsec_STV_nu") { + ReadGlobalConfigDefaults(); } void Initialize(fhicl::ParameterSet const &instance_sample_configuration) { - //! Set the verbosity of the sample logging macros. - if (instance_sample_configuration.has_key("verbosity")) { - SetSampleVerbosity( - instance_sample_configuration.get("verbosity")); - } - - STVProjection proj = - GetProjection(instance_sample_configuration.get( - "projection", "DeltaPT")); + // Perform any per-sample configuration in the base class + MultiDataComparison::Initialize(instance_sample_configuration); + + std::unique_ptr DPT( + new SimpleDataComparison_1D("T2K_CC0PiNProt_CH_xsec_STV_nu:dpt")); + DPT->SetData( + GetDataDir() + + "nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/dptResults.root;Result"); + DPT->Initialize( + instance_sample_configuration.get("dpt", {})); + DPT->CompProjFunc = [](FullEvent const &fev) + -> std::array { + return {GetDeltaPT_CC0PiN(fev).Mag() * 1E-3}; + }; - //! This will automatically set the data histogram to be loaded. - SetData(instance_sample_configuration.get("data_specifier")); + std::unique_ptr DAT( + new SimpleDataComparison_1D("T2K_CC0PiNProt_CH_xsec_STV_nu:dat")); + DAT->SetData( + GetDataDir() + + "nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/datResults.root;Result"); + DAT->Initialize( + instance_sample_configuration.get("dat", {})); + DAT->CompProjFunc = [](FullEvent const &fev) + -> std::array { + return {GetDeltaAlphaT_CC0PiN(fev)}; + }; - std::string projstr = ""; - switch (proj) { - case kDeltaPT: { - projstr = "dpt"; - break; - } - case kDeltaPhiT: { - projstr = "dphit"; - break; - } - case kDeltaAlphaT: { - projstr = "dat"; - break; - } - } + std::unique_ptr DPhiT( + new SimpleDataComparison_1D("T2K_CC0PiNProt_CH_xsec_STV_nu:dphit")); + DPhiT->SetData( + GetDataDir() + + "nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/dphitResults.root;Result"); + DPhiT->Initialize( + instance_sample_configuration.get("dphit", {})); + DPhiT->CompProjFunc = [](FullEvent const &fev) + -> std::array { + return {GetDeltaPhiT_CC0PiN(fev)}; + }; - SetData(GetDataDir() + "nuA/Nuclear/T2K/CC0Pi/" + projstr + - "Results.root;Result"); + std::string wd_stub = write_directory.size() ? write_directory + "/" : ""; + DPT->write_directory = wd_stub + "dpt"; + DAT->write_directory = wd_stub + "dat"; + DPhiT->write_directory = wd_stub + "dphit"; - // Perform any per-sample configuration in the base class - SimpleDataComparison_1D::Initialize(instance_sample_configuration); + Comparisons.emplace_back("dpt", std::move(DPT)); + Comparisons.emplace_back("dat", std::move(DAT)); + Comparisons.emplace_back("dphit", std::move(DPhiT)); //! Define your event signal here. IsSigFunc = [](FullEvent const &fev) -> bool { - //! See src/utility/EventTopologyUtility.hxx for more pre-defined - //! topological signals. - if (!IsCC0Pi(fev)) { + // Is numuCC + if (GetNParticles(fev, {pdgcodes::kMu}) != 1) { return false; } - - //! See src/event/FullEvent.hxx for the full event class definition. - //! See src/utility/FullEventUtility.hxx for more helper methods for - //! interacting with the event class. - - //! Get the initial state muon neutrino - Particle ISNumu = GetHMISParticle(fev, {pdgcodes::kNuMu}); - //! An nuis::event::Particle that return true for !part is invalid and - //! does not exist on the particle stack. i.e. here, the selection fails - //! if the event didn't have an initial state numubar. - if (!ISNumu) { - return false; - } - - //! Get the final state muon - Particle FSMu = GetHMFSParticle(fev, {pdgcodes::kMu}); - if (!FSMu) { - return false; - } - - //! Get the highest momentum final state proton - Particle FSProton = GetHMFSParticle(fev, {pdgcodes::kProton}); - if (!FSProton) { - return false; - } - - //! Cut on kinematic properties of the true final state. - - // Muon phase space - // Pmu > 250 MeV, cos(theta_mu) > -0.6 (Sweet phase space!) - if ((FSMu.P() < 250) || (cos(ISNumu.P3().Angle(FSMu.P3())) < -0.6)) { - return false; - } - - // Proton phase space - // Pprot > 450 MeV, cos(theta_proton) > 0.4 - if ((FSProton.P() < 450) || (FSProton.P() > 1E3) || - (cos(ISNumu.P3().Angle(FSProton.P3())) < 0.4)) { - return false; - } - - //! Select the event! - return true; - }; - - //! 1D Projection function - //! This function takes selected events and returns an array of size, the - //! dimensionality of the comparisons (SimpleDataComparison_1D::NDim) - CompProjFunc = [=](FullEvent const &fev) - -> std::array { - switch (proj) { - case kDeltaPT: { - return {GetDeltaPT_CC0PiN(fev).Mag() * 1E-3}; - } - case kDeltaPhiT: { - return {GetDeltaPhiT_CC0PiN(fev)}; - } - case kDeltaAlphaT: { - return {GetDeltaAlphaT_CC0PiN(fev)}; - } - } - return {-std::numeric_limits::max()}; + return t2k::IsCC0Pi_STV(fev); }; } - std::string Name() { return "T2K_CC0PiNProt_CH_xsec_STV_nu"; } - //! Here you can write any custom histograms to TTrees that your sample has //! been handling. - void Write() { SimpleDataComparison_1D::Write(); } + void Write() { + for (auto &comp : Comparisons) { + comp.second->Write(); + } + } }; //! These declarations allow your class to be loaded dynamically by NUISANCE DECLARE_PLUGIN(IDataComparison, T2K_CC0PiNProt_CH_xsec_STV_nu); DECLARE_PLUGIN(IEventProcessor, T2K_CC0PiNProt_CH_xsec_STV_nu); diff --git a/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar.cxx b/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar.cxx index a225472..1aab8d6 100644 --- a/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar.cxx +++ b/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar.cxx @@ -1,182 +1,183 @@ // Copyright 2018 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 "samples/SimpleDataComparison.hxx" #include "utility/EventTopologyUtility.hxx" #include "utility/FullEventUtility.hxx" #include "utility/KinematicUtility.hxx" #include "utility/PDGCodeUtility.hxx" using namespace nuis::event; using namespace nuis::utility; class T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar : public SimpleDataComparison_2DPoly { std::unique_ptr fPrediction_pmu; std::unique_ptr fPrediction_ctheta; std::unique_ptr fPrediction_fine; std::unique_ptr fPrediction_unscale; public: - T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar() { ReadGlobalConfigDefaults(); } + T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar() + : SimpleDataComparison_2DPoly("T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar") { + ReadGlobalConfigDefaults(); + } std::string GetDocumentation() { return ""; } fhicl::ParameterSet GetExampleConfiguration() { fhicl::ParameterSet exps = SimpleDataComparison_2DPoly::GetExampleConfiguration(); return exps; } void Initialize(fhicl::ParameterSet const &instance_sample_configuration) { SetData( GetDataDir() + "nuA/Nuclear/T2K/CC0Pi/H2O_xsec_2Dpmuthetamu_numubar.root;datapoly"); fPrediction_fine = std::unique_ptr( new TH2D("fPrediction_fine", ";#it{p}_{#mu} " "(MeV/#it{c});cos(#theta_{#mu});d^{2}#sigma/" "d#it{p}_{#mu}dcos(#theta_{#mu}) (cm^{2} MeV^{-1} A^{-1})", 100, 0, 3000, 50, 0.84, 1)); fPrediction_pmu = std::unique_ptr(new TH1D("fPrediction_pmu", ";#it{p}_{#mu} " "(MeV/#it{c});d#sigma/" "d#it{p}_{#mu} (cm^{2} MeV^{-1} A^{-1})", 100, 0, 3000)); fPrediction_ctheta = std::unique_ptr(new TH1D("fPrediction_ctheta", ";cos(#theta_{#mu});d#sigma/" "ddcos(#theta_{#mu}) (cm^{2} A^{-1})", 50, 0.84, 1)); SimpleDataComparison_2DPoly::Initialize(instance_sample_configuration); fPrediction_unscale = nuis::utility::Clone(fPrediction, false, "Prediction_unscale"); // Signal selection function IsSigFunc = [](FullEvent const &fev) -> bool { if (!IsCC0Pi(fev)) { return false; } Particle ISNumuBar = GetHMISParticle(fev, {pdgcodes::kNuMuBar}); if (!ISNumuBar) { return false; } Particle FSMuPlus = GetHMFSParticle(fev, {pdgcodes::kMuPlus}); if (!FSMuPlus) { return false; } if (FSMuPlus.CosTheta() < 0.84) { return false; } return true; }; // 1D Projection function CompProjFunc = [](FullEvent const &fev) -> std::array { Particle FSMuPlus = GetHMFSParticle(fev, {pdgcodes::kMuPlus}); return {FSMuPlus.P(), FSMuPlus.CosTheta()}; }; ProcessExtraFunc = [&](FullEvent const &fev, bool isSel, double weight) { if (isSel) { Particle FSMuPlus = GetHMFSParticle(fev, {pdgcodes::kMuPlus}); TH_Helper::Fill(fPrediction_fine, {FSMuPlus.P(), FSMuPlus.CosTheta()}, weight); TH_Help::Fill(fPrediction_unscale, {FSMuPlus.P(), FSMuPlus.CosTheta()}); } Particle ISNumuBar = GetHMISParticle(fev, {pdgcodes::kNuMuBar}); if (!ISNumuBar) { return; } Particle FSMuPlus = GetHMFSParticle(fev, {pdgcodes::kMuPlus}); if (!FSMuPlus) { return; } TH_Helper::Fill(fPrediction_pmu, {FSMuPlus.P()}, weight); TH_Helper::Fill(fPrediction_ctheta, {FSMuPlus.CosTheta()}, weight); }; } - std::string Name() { return "T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar"; } - void Write() { SimpleDataComparison_2DPoly::Write(); nuis::persistency::WriteToOutputFile( fPrediction_fine, "Prediction_fine", write_directory); nuis::persistency::WriteToOutputFile( fPrediction_pmu, "fPrediction_pmu", write_directory); nuis::persistency::WriteToOutputFile( fPrediction_ctheta, "fPrediction_ctheta", write_directory); nuis::persistency::WriteToOutputFile( fPrediction_unscale, "Prediction_unscale", write_directory); static std::vector> const binning = { {// Slice_0 YPolyBinSpec(450, 0.85), YPolyBinSpec(450, 0.95)}, {// Slice_1 YPolyBinSpec(550, 0.86), YPolyBinSpec(550, 0.93), YPolyBinSpec(550, 0.97)}, {// Slice_2 YPolyBinSpec(700, 0.89), YPolyBinSpec(700, 0.94), YPolyBinSpec(700, 0.98)}, {// Slice_3 YPolyBinSpec(900, 0.91), YPolyBinSpec(900, 0.95), YPolyBinSpec(900, 0.98)}, {// Slice_4 YPolyBinSpec(1200, 0.92), YPolyBinSpec(1200, 0.96), YPolyBinSpec(1200, 0.98)}, {// Slice_4 YPolyBinSpec(1600, 0.93), YPolyBinSpec(1600, 0.97), YPolyBinSpec(1600, 0.99)}, {// Slice_5 YPolyBinSpec(2500, 0.96), YPolyBinSpec(2500, 0.99)}, }; for (auto &slice : GetTH2PolySlices(fData, binning)) { nuis::persistency::WriteToOutputFile(slice, slice->GetName(), write_directory); } for (auto &slice : GetTH2PolySlices(fPrediction_xsec, binning)) { nuis::persistency::WriteToOutputFile(slice, slice->GetName(), write_directory); } } }; DECLARE_PLUGIN(IDataComparison, T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar); DECLARE_PLUGIN(IEventProcessor, T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar); diff --git a/src/samples/nuA/Nuclear/T2K/CC1Pi/CMakeLists.txt b/src/samples/nuA/Nuclear/T2K/CC1Pi/CMakeLists.txt new file mode 100644 index 0000000..b023764 --- /dev/null +++ b/src/samples/nuA/Nuclear/T2K/CC1Pi/CMakeLists.txt @@ -0,0 +1,18 @@ +SET(SAMPLES +T2K_CC1Pip_CH_xsec_muproj_nu +T2K_CC1Pip_CH_xsec_piproj_nu) + +foreach(S ${SAMPLES}) + LIST(APPEND INuADataComparisons ${S}) + LIST(APPEND SAMPLES_src ${S}.cxx) +endforeach() + +add_library(T2KCC1PiDataComparisons SHARED ${SAMPLES_src}) +target_link_libraries(T2KCC1PiDataComparisons + nuis_event nuis_input nuis_config nuis_persistency) + +install(TARGETS T2KCC1PiDataComparisons DESTINATION plugins) + +install(FILES T2K.CC1Pi.sample.config.fcl DESTINATION fcl) + +SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) diff --git a/src/samples/nuA/Nuclear/T2K/CC1Pi/T2K.CC1Pi.sample.config.fcl b/src/samples/nuA/Nuclear/T2K/CC1Pi/T2K.CC1Pi.sample.config.fcl new file mode 100644 index 0000000..9796eeb --- /dev/null +++ b/src/samples/nuA/Nuclear/T2K/CC1Pi/T2K.CC1Pi.sample.config.fcl @@ -0,0 +1,7 @@ +global.sample_configuration.T2K_CC1Pip_CH_xsec: { + journal_reference: "Not yet published" + doi: "https://doi.org/10.1103/PhysRevD.98.032003" + target_material: "CH" + flux_description: "ND280 FHC muon neutrino" + signal_description: "CC1Pi" +} diff --git a/src/samples/nuA/Nuclear/T2K/CC1Pi/T2K_CC1Pip_CH_xsec_muproj_nu.cxx b/src/samples/nuA/Nuclear/T2K/CC1Pi/T2K_CC1Pip_CH_xsec_muproj_nu.cxx new file mode 100644 index 0000000..f9fa284 --- /dev/null +++ b/src/samples/nuA/Nuclear/T2K/CC1Pi/T2K_CC1Pip_CH_xsec_muproj_nu.cxx @@ -0,0 +1,95 @@ +// Copyright 2018 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 "samples/MultiDataComparison.hxx" +#include "samples/SimpleDataComparison.hxx" + +#include "utility/EventTopologyUtility.hxx" +#include "utility/FullEventUtility.hxx" +#include "utility/KinematicUtility.hxx" +#include "utility/PDGCodeUtility.hxx" +#include "utility/experimental/T2KUtility.hxx" + +using namespace nuis::event; +using namespace nuis::utility; + +class T2K_CC1Pip_CH_xsec_muproj_nu : public MultiDataComparison { + +public: + T2K_CC1Pip_CH_xsec_muproj_nu() + : MultiDataComparison("T2K_CC1Pip_CH_xsec_muproj_nu") { + ReadGlobalConfigDefaults(); + } + + enum distribution { + kMomentumPion, kPmuThetamu, kPhi_adler, kQ2, kTheta_adler, kThetapimu, kThetapion + }; + + void Initialize(fhicl::ParameterSet const &instance_sample_configuration) { + + // Perform any per-sample configuration in the base class + MultiDataComparison::Initialize(instance_sample_configuration); + +// PmuThetamu.root +// MomentumPion.root Thetapimu.root +// MomentumPion.txt Q2.root Thetapimu.txt +// phi_adler.root Q2.txt Thetapion.root +// phi_adler.txt theta_adler.root Thetapion.txt + + std::unique_ptr MomentumPion( + new SimpleDataComparison_1D("T2K_CC1Pip_CH_xsec_muproj_nu:MomentumPion")); + MomentumPion->SetData( + GetDataDir() + + "nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/MomentumPion.root;Momentum_pion (GeV)"); + MomentumPion->Initialize( + instance_sample_configuration.get("MomentumPion", {})); + MomentumPion->CompProjFunc = [](FullEvent const &fev) + -> std::array { + return {GetHMFSParticle(fev,{pdgcodes::kPiPlus}).P() * 1E-3}; + }; + + std::string wd_stub = write_directory.size() ? write_directory + "/" : ""; + MomentumPion->write_directory = wd_stub + "MomentumPion"; + + Comparisons.emplace_back("MomentumPion", std::move(MomentumPion)); + + //! Define your event signal here. + IsSigFunc = [](FullEvent const &fev) -> bool { + // Is numuCC + if (GetNParticles(fev, {pdgcodes::kMu}) != 1) { + return false; + } + return t2k::IsCC1Pip_CH_MichTag(fev); + }; + } + + //! Here you can write any custom histograms to TTrees that your sample has + //! been handling. + void Write() { + for (auto &comp : Comparisons) { + comp.second->Write(); + } + } +}; + +//! These declarations allow your class to be loaded dynamically by NUISANCE +DECLARE_PLUGIN(IDataComparison, T2K_CC1Pip_CH_xsec_muproj_nu); +DECLARE_PLUGIN(IEventProcessor, T2K_CC1Pip_CH_xsec_muproj_nu); diff --git a/src/samples/nuA/Nuclear/T2K/CC1Pi/T2K_CC1Pip_CH_xsec_piproj_nu.cxx b/src/samples/nuA/Nuclear/T2K/CC1Pi/T2K_CC1Pip_CH_xsec_piproj_nu.cxx new file mode 100644 index 0000000..a0b5121 --- /dev/null +++ b/src/samples/nuA/Nuclear/T2K/CC1Pi/T2K_CC1Pip_CH_xsec_piproj_nu.cxx @@ -0,0 +1,95 @@ +// Copyright 2018 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 "samples/MultiDataComparison.hxx" +#include "samples/SimpleDataComparison.hxx" + +#include "utility/EventTopologyUtility.hxx" +#include "utility/FullEventUtility.hxx" +#include "utility/KinematicUtility.hxx" +#include "utility/PDGCodeUtility.hxx" +#include "utility/experimental/T2KUtility.hxx" + +using namespace nuis::event; +using namespace nuis::utility; + +class T2K_CC1Pip_CH_xsec_piproj_nu : public MultiDataComparison { + +public: + T2K_CC1Pip_CH_xsec_piproj_nu() + : MultiDataComparison("T2K_CC1Pip_CH_xsec_piproj_nu") { + ReadGlobalConfigDefaults(); + } + + enum distribution { + kMomentumPion, kPmuThetamu, kPhi_adler, kQ2, kTheta_adler, kThetapimu, kThetapion + }; + + void Initialize(fhicl::ParameterSet const &instance_sample_configuration) { + + // Perform any per-sample configuration in the base class + MultiDataComparison::Initialize(instance_sample_configuration); + +// PmuThetamu.root +// MomentumPion.root Thetapimu.root +// MomentumPion.txt Q2.root Thetapimu.txt +// phi_adler.root Q2.txt Thetapion.root +// phi_adler.txt theta_adler.root Thetapion.txt + + std::unique_ptr MomentumPion( + new SimpleDataComparison_1D("T2K_CC1Pip_CH_xsec_piproj_nu:MomentumPion")); + MomentumPion->SetData( + GetDataDir() + + "nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/MomentumPion.root;Momentum_pion (GeV)"); + MomentumPion->Initialize( + instance_sample_configuration.get("MomentumPion", {})); + MomentumPion->CompProjFunc = [](FullEvent const &fev) + -> std::array { + return {GetHMFSParticle(fev,{pdgcodes::kPiPlus}).P() * 1E-3}; + }; + + std::string wd_stub = write_directory.size() ? write_directory + "/" : ""; + MomentumPion->write_directory = wd_stub + "MomentumPion"; + + Comparisons.emplace_back("MomentumPion", std::move(MomentumPion)); + + //! Define your event signal here. + IsSigFunc = [](FullEvent const &fev) -> bool { + // Is numuCC + if (GetNParticles(fev, {pdgcodes::kMu}) != 1) { + return false; + } + return t2k::IsCC1Pip_CH_RecPi(fev); + }; + } + + //! Here you can write any custom histograms to TTrees that your sample has + //! been handling. + void Write() { + for (auto &comp : Comparisons) { + comp.second->Write(); + } + } +}; + +//! These declarations allow your class to be loaded dynamically by NUISANCE +DECLARE_PLUGIN(IDataComparison, T2K_CC1Pip_CH_xsec_piproj_nu); +DECLARE_PLUGIN(IEventProcessor, T2K_CC1Pip_CH_xsec_piproj_nu); diff --git a/src/samples/nuA/Nuclear/T2K/CMakeLists.txt b/src/samples/nuA/Nuclear/T2K/CMakeLists.txt index b667c78..d91f6b9 100644 --- a/src/samples/nuA/Nuclear/T2K/CMakeLists.txt +++ b/src/samples/nuA/Nuclear/T2K/CMakeLists.txt @@ -1,7 +1,8 @@ add_subdirectory(CC0Pi) +add_subdirectory(CC1Pi) LIST(APPEND INuADataComparisons_FHiCL T2K.sample.config.fcl) install(FILES T2K.sample.config.fcl DESTINATION fcl) SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE) diff --git a/src/samples/nuA/Nuclear/T2K/T2K.sample.config.fcl b/src/samples/nuA/Nuclear/T2K/T2K.sample.config.fcl index 18a0fd0..2ee611f 100644 --- a/src/samples/nuA/Nuclear/T2K/T2K.sample.config.fcl +++ b/src/samples/nuA/Nuclear/T2K/T2K.sample.config.fcl @@ -1 +1,2 @@ #include "T2K.CC0Pi.sample.config.fcl" +#include "T2K.CC1Pi.sample.config.fcl" diff --git a/src/utility/CMakeLists.txt b/src/utility/CMakeLists.txt index 5f8c5d6..8024eba 100644 --- a/src/utility/CMakeLists.txt +++ b/src/utility/CMakeLists.txt @@ -1,28 +1,32 @@ set(utility_implementation_files FileSystemUtility.cxx FullEventUtility.cxx KinematicUtility.cxx PDGCodeUtility.cxx StringUtility.cxx EventTopologyUtility.cxx - HistogramUtility.cxx) + HistogramUtility.cxx + PhaseSpaceRestriction.cxx) set(utility_header_files FileSystemUtility.hxx FullEventUtility.hxx InteractionChannelUtility.hxx KinematicUtility.hxx PDGCodeUtility.hxx ROOTUtility.hxx HistogramUtility.hxx StringUtility.hxx TerminalUtility.hxx StatsUtility.hxx - EventTopologyUtility.hxx) + EventTopologyUtility.hxx + PhaseSpaceRestriction.cxx) add_library(nuis_utility SHARED ${utility_implementation_files}) target_link_libraries(nuis_utility) install(TARGETS nuis_utility DESTINATION lib) install(FILES ${utility_header_files} DESTINATION include/utility) + +add_subdirectory(experimental) diff --git a/src/utility/EventTopologyUtility.cxx b/src/utility/EventTopologyUtility.cxx index 6ba1869..b68d4e3 100644 --- a/src/utility/EventTopologyUtility.cxx +++ b/src/utility/EventTopologyUtility.cxx @@ -1,33 +1,51 @@ #include "utility/EventTopologyUtility.hxx" #include "event/FullEvent.hxx" #include "utility/FullEventUtility.hxx" namespace nuis { namespace utility { bool IsCC0Pi(event::FullEvent const &ev) { + if (!IsCCInc(ev)) { + return false; + } - event::Particle ISNu = GetHMISNeutralLepton(ev); - if (!ISNu) { + if (GetNFSPions(ev) || GetNFSOthers(ev)) { return false; } - event::PDG_t expected_fslep_pdg = ISNu.pdg > 0 ? ISNu.pdg - 1 : ISNu.pdg + 1; + return true; +} +bool IsCC1Pi(event::FullEvent const &ev, std::vector PionPDGs) { + if (!IsCCInc(ev)) { + return false; + } - event::Particle FSLep = GetHMFSParticle(ev, {expected_fslep_pdg}); - if (!FSLep) { + size_t NChosenPi = GetNParticles(ev, PionPDGs); + if ((NChosenPi != 1) || (NChosenPi != GetNFSPions(ev)) || GetNFSOthers(ev)) { return false; } - if (GetFSPions(ev).size() || GetFSOthers(ev).size()) { + return true; +} + +bool IsCCInc(event::FullEvent const &ev) { + event::Particle ISNu = GetHMISNeutralLepton(ev); + if (!ISNu) { + return false; + } + + event::PDG_t expected_fslep_pdg = ISNu.pdg > 0 ? ISNu.pdg - 1 : ISNu.pdg + 1; + + if (GetNParticles(ev, {expected_fslep_pdg}) != GetNFSLeptons(ev)) { return false; } return true; } } // namespace utility } // namespace nuis diff --git a/src/utility/EventTopologyUtility.hxx b/src/utility/EventTopologyUtility.hxx index 2043e9d..1f8fee1 100644 --- a/src/utility/EventTopologyUtility.hxx +++ b/src/utility/EventTopologyUtility.hxx @@ -1,15 +1,19 @@ #pragma once +#include "utility/PDGCodeUtility.hxx" + namespace nuis { namespace event { class FullEvent; } // namespace event } // namespace nuis namespace nuis { namespace utility { bool IsCC0Pi(event::FullEvent const &); + bool IsCC1Pi(event::FullEvent const &, std::vector PionPDGs = pdgcodes::Pions); + bool IsCCInc(event::FullEvent const &); } } diff --git a/src/utility/FileSystemUtility.cxx b/src/utility/FileSystemUtility.cxx index a9acbb4..8f04cfb 100644 --- a/src/utility/FileSystemUtility.cxx +++ b/src/utility/FileSystemUtility.cxx @@ -1,53 +1,46 @@ #include "utility/FileSystemUtility.hxx" +#include "utility/StringUtility.hxx" #include "exception/exception.hxx" #include "config/GlobalConfiguration.hxx" #include "fhiclcpp/ParameterSet.h" #include -#include - namespace nuis { namespace utility { -NEW_NUIS_EXCEPT(unexpected_empty_string); - -std::string EnsureTrailingSlash(std::string str) { - if (!str.size()) { - throw unexpected_empty_string(); - } - if (str.back() != '/') { - return str + '/'; - } - return str; -} +std::vector GetMatchingFiles(std::string directory, + std::string pattern) { -std::vector GetMatchingFiles(std::string directory, std::string const &pattern) { + std::cout << "[INFO]: Looking for files matching: \"" << pattern + << "\" in directory: " << directory << std::endl; directory = EnsureTrailingSlash(directory); + pattern = DeGlobPattern(pattern); + std::regex rpattern(pattern); std::vector matches; DIR *dir; struct dirent *ent; dir = opendir(directory.c_str()); if (dir != NULL) { while ((ent = readdir(dir)) != NULL) { if (std::regex_match(ent->d_name, rpattern)) { matches.push_back(ent->d_name); } } } return matches; } std::string GetDataDir() { return EnsureTrailingSlash( config::GetDocument().get("data_dir")); } } // namespace utility } // namespace nuis diff --git a/src/utility/FileSystemUtility.hxx b/src/utility/FileSystemUtility.hxx index eb7f7a0..df7b270 100644 --- a/src/utility/FileSystemUtility.hxx +++ b/src/utility/FileSystemUtility.hxx @@ -1,37 +1,32 @@ // Copyright 2018 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 UTILITY_FILESYSTEMUTILITY_HXX_SEEN -#define UTILITY_FILESYSTEMUTILITY_HXX_SEEN - +#pragma once #include #include namespace nuis { namespace utility { -std::string EnsureTrailingSlash(std::string); - -std::vector GetMatchingFiles(std::string, std::string const &); +std::vector GetMatchingFiles(std::string, std::string); std::string GetDataDir(); } // namespace utility } // namespace nuis -#endif diff --git a/src/utility/FullEventUtility.cxx b/src/utility/FullEventUtility.cxx index 7046b2e..2dc3c74 100644 --- a/src/utility/FullEventUtility.cxx +++ b/src/utility/FullEventUtility.cxx @@ -1,172 +1,297 @@ #include "utility/FullEventUtility.hxx" #include "event/FullEvent.hxx" #include "utility/PDGCodeUtility.hxx" +#include "utility/PhaseSpaceRestriction.hxx" using namespace nuis::event; namespace nuis { namespace utility { event::Particle GetHMParticle(std::vector const &parts) { if (parts.size()) { return event::Particle(); } return *std::max_element( parts.begin(), parts.end(), [](event::Particle const &l, event::Particle const &r) { return l.P() < r.P(); }); } std::vector GetParticles(FullEvent const &ev, std::vector const &pdgs, Particle::Status_t status, bool include_matching_pdg) { std::vector selected_particles; for (auto const &parts : ev.ParticleStack) { if (parts.status != status) { continue; } for (Particle const &part : parts.particles) { bool matched_pdg = false; for (PDG_t pdg : pdgs) { matched_pdg = matched_pdg || (part.pdg == pdg); } bool keep = ((include_matching_pdg && matched_pdg) || (!include_matching_pdg && !matched_pdg)); if (!keep) { continue; } selected_particles.push_back(part); } } return selected_particles; } std::vector const &GetISParticles(FullEvent const &ev) { return ev .ParticleStack[static_cast( Particle::Status_t::kPrimaryInitialState)] .particles; } std::vector const &GetPrimaryFSParticles(FullEvent const &ev) { return ev .ParticleStack[static_cast( Particle::Status_t::kPrimaryFinalState)] .particles; } std::vector const &GetNuclearLeavingParticles(FullEvent const &ev) { return ev .ParticleStack[static_cast(Particle::Status_t::kNuclearLeaving)] .particles; } Particle GetHMParticle(FullEvent const &ev, std::vector const &pdgs, Particle::Status_t status, bool include_matching_pdg) { Particle HMParticle; for (auto const &parts : ev.ParticleStack) { if (parts.status != status) { continue; } for (Particle const &part : parts.particles) { bool matched_pdg = false; for (PDG_t pdg : pdgs) { matched_pdg = matched_pdg || (part.pdg == pdg); } bool keep = ((include_matching_pdg && matched_pdg) || (!include_matching_pdg && !matched_pdg)); if (!keep) { continue; } if (part.P() > HMParticle.P()) { HMParticle = part; } } } return HMParticle; } event::Particle GetHMISParticle(event::FullEvent const &ev, std::vector const &pdgs) { return GetHMParticle(ev, pdgs, Particle::Status_t::kPrimaryInitialState); } event::Particle GetHMFSParticle(event::FullEvent const &ev, std::vector const &pdgs) { return GetHMParticle(ev, pdgs); } std::vector GetFSChargedLeptons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::ChargedLeptons); } std::vector GetFSNeutralLeptons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::NeutralLeptons); } std::vector GetISNeutralLeptons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::NeutralLeptons, Particle::Status_t::kPrimaryInitialState); } std::vector GetFSChargedPions(FullEvent const &ev) { return GetParticles(ev, pdgcodes::ChargedPions); } std::vector GetFSNeutralPions(FullEvent const &ev) { return GetParticles(ev, pdgcodes::NeutralPions); } std::vector GetFSPions(FullEvent const &ev) { return GetParticles(ev, pdgcodes::Pions); } std::vector GetFSProtons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::Protons); } std::vector GetFSNeutrons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::Neutron); } std::vector GetFSNucleons(FullEvent const &ev) { return GetParticles(ev, pdgcodes::Nucleons); } std::vector GetFSOthers(FullEvent const &ev) { return GetParticles(ev, pdgcodes::CommonParticles, Particle::Status_t::kNuclearLeaving, false); } Particle GetHMFSChargedLepton(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::ChargedLeptons); } Particle GetHMFSNeutralLepton(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::NeutralLeptons); } Particle GetHMFSLepton(event::FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Leptons); } Particle GetHMISNeutralLepton(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::NeutralLeptons, Particle::Status_t::kPrimaryInitialState); } Particle GetHMFSChargedPion(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::ChargedPions); } Particle GetHMFSNeutralPion(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::NeutralPions); } Particle GetHMFSPion(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Pions); } Particle GetHMFSProton(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Protons); } Particle GetHMFSNeutron(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Neutron); } Particle GetHMFSNucleon(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Nucleons); } Particle GetHMFSOther(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::CommonParticles, Particle::Status_t::kNuclearLeaving, false); } +Particle GetHMParticleInPhaseSpace(FullEvent const &ev, + std::vector const &pdgs, + IParticlePhaseSpaceRestriction const &ps, + Particle::Status_t status, + bool include_matching_pdg) { + Particle HMParticle; + for (auto const &parts : ev.ParticleStack) { + if (parts.status != status) { + continue; + } + for (Particle const &part : parts.particles) { + if (!ps.Inside(part)) { + continue; + } + bool matched_pdg = false; + for (PDG_t pdg : pdgs) { + matched_pdg = matched_pdg || (part.pdg == pdg); + } + bool keep = ((include_matching_pdg && matched_pdg) || + (!include_matching_pdg && !matched_pdg)); + if (!keep) { + continue; + } + if (part.P() > HMParticle.P()) { + HMParticle = part; + } + } + } + return HMParticle; +} + +Particle GetHMFSProtonInPhaseSpace(FullEvent const &ev, + IParticlePhaseSpaceRestriction const &ps) { + return GetHMParticleInPhaseSpace(ev, pdgcodes::Protons, ps); +} + +size_t GetNParticles(FullEvent const &ev, std::vector const &pdgs, + Particle::Status_t status, bool include_matching_pdg) { + size_t N = 0; + for (auto const &parts : ev.ParticleStack) { + if (parts.status != status) { + continue; + } + for (Particle const &part : parts.particles) { + bool matched_pdg = false; + for (PDG_t pdg : pdgs) { + matched_pdg = matched_pdg || (part.pdg == pdg); + } + bool keep = ((include_matching_pdg && matched_pdg) || + (!include_matching_pdg && !matched_pdg)); + if (!keep) { + continue; + } + N++; + } + } + return N; +} + +size_t GetNParticlesInPhaseSpace(FullEvent const &ev, + std::vector const &pdgs, + IParticlePhaseSpaceRestriction const &ps, + Particle::Status_t status, + bool include_matching_pdg) { + size_t N = 0; + for (auto const &parts : ev.ParticleStack) { + if (parts.status != status) { + continue; + } + for (Particle const &part : parts.particles) { + if (!ps.Inside(part)) { + continue; + } + bool matched_pdg = false; + for (PDG_t pdg : pdgs) { + matched_pdg = matched_pdg || (part.pdg == pdg); + } + bool keep = ((include_matching_pdg && matched_pdg) || + (!include_matching_pdg && !matched_pdg)); + if (!keep) { + continue; + } + N++; + } + } + return N; +} + +size_t GetNFSChargedLeptons(FullEvent const &ev) { + return GetNParticles(ev, pdgcodes::ChargedLeptons); +} +size_t GetNFSNeutralLeptons(FullEvent const &ev) { + return GetNParticles(ev, pdgcodes::NeutralLeptons); +} +size_t GetNFSLeptons(FullEvent const &ev) { + return GetNParticles(ev, pdgcodes::Leptons); +} +size_t GetNISNeutralLeptons(FullEvent const &ev) { + return GetNParticles(ev, pdgcodes::NeutralLeptons, + Particle::Status_t::kPrimaryInitialState); +} +size_t GetNFSChargedPions(FullEvent const &ev) { + return GetNParticles(ev, pdgcodes::ChargedPions); +} +size_t GetNFSNeutralPions(FullEvent const &ev) { + return GetNParticles(ev, pdgcodes::NeutralPions); +} +size_t GetNFSPions(FullEvent const &ev) { + return GetNParticles(ev, pdgcodes::Pions); +} +size_t GetNFSProtons(FullEvent const &ev) { + return GetNParticles(ev, pdgcodes::Protons); +} +size_t GetNFSNeutrons(FullEvent const &ev) { + return GetNParticles(ev, pdgcodes::Neutron); +} +size_t GetNFSNucleons(FullEvent const &ev) { + return GetNParticles(ev, pdgcodes::Nucleons); +} +size_t GetNFSOthers(FullEvent const &ev) { + return GetNParticles(ev, pdgcodes::CommonParticles, + Particle::Status_t::kNuclearLeaving, false); +} + } // namespace utility } // namespace nuis diff --git a/src/utility/FullEventUtility.hxx b/src/utility/FullEventUtility.hxx index 22c7ac9..3b32390 100644 --- a/src/utility/FullEventUtility.hxx +++ b/src/utility/FullEventUtility.hxx @@ -1,85 +1,124 @@ // Copyright 2018 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 . *******************************************************************************/ #pragma once -#include "event/types.hxx" #include "event/Particle.hxx" +#include "event/types.hxx" #include namespace nuis { namespace event { class FullEvent; } // namespace event } // namespace nuis namespace nuis { namespace utility { +class IParticlePhaseSpaceRestriction; event::Particle GetHMParticle(std::vector); std::vector GetParticles(event::FullEvent const &, std::vector const &, event::Particle::Status_t status = event::Particle::Status_t::kNuclearLeaving, bool include_matching_pdg = true); std::vector const &GetISParticles(event::FullEvent const &); std::vector const & GetPrimaryFSParticles(event::FullEvent const &); std::vector const & GetNuclearLeavingParticles(event::FullEvent const &); event::Particle GetHMParticle(event::FullEvent const &, std::vector const &, event::Particle::Status_t status = event::Particle::Status_t::kNuclearLeaving, bool include_matching_pdg = true); event::Particle GetHMISParticle(event::FullEvent const &, - std::vector const &); + std::vector const &); event::Particle GetHMFSParticle(event::FullEvent const &, - std::vector const &); + std::vector const &); + +///***********Getters that include copies std::vector GetFSChargedLeptons(event::FullEvent const &); std::vector GetFSNeutralLeptons(event::FullEvent const &); std::vector GetISNeutralLeptons(event::FullEvent const &); std::vector GetFSChargedPions(event::FullEvent const &); std::vector GetFSNeutralPions(event::FullEvent const &); std::vector GetFSPions(event::FullEvent const &); std::vector GetFSProtons(event::FullEvent const &); std::vector GetFSNeutrons(event::FullEvent const &); std::vector GetFSNucleons(event::FullEvent const &); std::vector GetFSOthers(event::FullEvent const &); event::Particle GetHMFSChargedLepton(event::FullEvent const &); event::Particle GetHMFSNeutralLepton(event::FullEvent const &); event::Particle GetHMFSLepton(event::FullEvent const &); event::Particle GetHMISNeutralLepton(event::FullEvent const &); event::Particle GetHMFSChargedPion(event::FullEvent const &); event::Particle GetHMFSNeutralPion(event::FullEvent const &); event::Particle GetHMFSPion(event::FullEvent const &); event::Particle GetHMFSProton(event::FullEvent const &); event::Particle GetHMFSNeutron(event::FullEvent const &); event::Particle GetHMFSNucleon(event::FullEvent const &); event::Particle GetHMFSOther(event::FullEvent const &); +event::Particle +GetHMParticleInPhaseSpace(event::FullEvent const &ev, + std::vector const &pdgs, + IParticlePhaseSpaceRestriction const &ps, + event::Particle::Status_t status = + event::Particle::Status_t::kNuclearLeaving, + bool include_matching_pdg = true); +event::Particle +GetHMFSProtonInPhaseSpace(event::FullEvent const &ev, + IParticlePhaseSpaceRestriction const &ps); + +///***********Counters requiring no copies +size_t GetNParticles(event::FullEvent const &ev, + std::vector const &pdgs, + event::Particle::Status_t status = + event::Particle::Status_t::kNuclearLeaving, + bool include_matching_pdg = true); +size_t GetNParticlesInPhaseSpace(event::FullEvent const &ev, + std::vector const &pdgs, + IParticlePhaseSpaceRestriction const &ps, + event::Particle::Status_t status = + event::Particle::Status_t::kNuclearLeaving, + bool include_matching_pdg = true); + +size_t GetNFSChargedLeptons(event::FullEvent const &ev); +size_t GetNFSNeutralLeptons(event::FullEvent const &ev); +size_t GetNFSLeptons(event::FullEvent const &ev); +size_t GetNISNeutralLeptons(event::FullEvent const &ev); +size_t GetNFSChargedPions(event::FullEvent const &ev); +size_t GetNFSNeutralPions(event::FullEvent const &ev); +size_t GetNFSPions(event::FullEvent const &ev); +size_t GetNFSProtons(event::FullEvent const &ev); +size_t GetNFSNeutrons(event::FullEvent const &ev); +size_t GetNFSNucleons(event::FullEvent const &ev); +size_t GetNFSOthers(event::FullEvent const &ev); + } // namespace utility } // namespace nuis diff --git a/src/utility/HistogramUtility.hxx b/src/utility/HistogramUtility.hxx index c84d620..d9bb5d8 100644 --- a/src/utility/HistogramUtility.hxx +++ b/src/utility/HistogramUtility.hxx @@ -1,467 +1,471 @@ // Copyright 2018 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 . *******************************************************************************/ #pragma once #include "utility/ROOTUtility.hxx" #include "exception/exception.hxx" #include "string_parsers/from_string.hxx" #include "fhiclcpp/ParameterSet.h" #include "TAxis.h" #include "TH1D.h" #include "TH1F.h" #include "TH2D.h" #include "TH2F.h" #include "TH2Poly.h" #include #include #include #include namespace nuis { namespace utility { NEW_NUIS_EXCEPT(unimplemented_GetHistogram_method); NEW_NUIS_EXCEPT(invalid_histogram_descriptor); NEW_NUIS_EXCEPT(invalid_histogram_name); NEW_NUIS_EXCEPT(failed_to_clone); bool IsFlowBin(TAxis const &ax, Int_t bin_it); bool IsInHistogramRange(TAxis const &ax, double v); template struct HType_traits {}; template <> struct HType_traits { using type = TH1; static size_t const NDim = 1; using NumericT = double; static std::string name() { return "TH1"; } }; template <> struct HType_traits { using type = TH1D; static size_t const NDim = 1; using NumericT = double; static std::string name() { return "TH1D"; } }; template <> struct HType_traits { using type = TH1F; static size_t const NDim = 1; using NumericT = float; static std::string name() { return "TH1F"; } }; template <> struct HType_traits { using type = TH2; static size_t const NDim = 2; using NumericT = double; static std::string name() { return "TH2"; } }; template <> struct HType_traits { using type = TH2D; static size_t const NDim = 2; using NumericT = double; static std::string name() { return "TH2D"; } }; template <> struct HType_traits { using type = TH2F; static size_t const NDim = 2; using NumericT = float; static std::string name() { return "TH2F"; } }; template <> struct HType_traits { using type = TH2Poly; static size_t const NDim = 0; using NumericT = double; static std::string name() { return "TH2Poly"; } }; template struct HType_Helper {}; template <> struct HType_Helper<1, void> { using type = TH1; static size_t const NDim = HType_traits::NDim; using NumericT = HType_traits::NumericT; static std::string name() { return HType_traits::name(); } }; template <> struct HType_Helper<1, double> { using type = TH1D; static size_t const NDim = HType_traits::NDim; using NumericT = HType_traits::NumericT; static std::string name() { return HType_traits::name(); } }; template <> struct HType_Helper<1, float> { using type = TH1F; static size_t const NDim = HType_traits::NDim; using NumericT = HType_traits::NumericT; static std::string name() { return HType_traits::name(); } }; template <> struct HType_Helper<2, void> { using type = TH2; static size_t const NDim = HType_traits::NDim; using NumericT = HType_traits::NumericT; static std::string name() { return HType_traits::name(); } }; template <> struct HType_Helper<2, double> { using type = TH2D; static size_t const NDim = HType_traits::NDim; using NumericT = HType_traits::NumericT; static std::string name() { return HType_traits::name(); } }; template <> struct HType_Helper<2, float> { using type = TH2F; static size_t const NDim = HType_traits::NDim; using NumericT = HType_traits::NumericT; static std::string name() { return HType_traits::name(); } }; template struct TH_Helper {}; template struct TH_Helper::NDim == 1>::type> { static size_t const NDim = HType_traits::NDim; using NumericT = typename HType_traits::NumericT; static std::string name() { return HType_traits::name(); } static Int_t NbinsIncludeFlow(HT const &h) { return h.GetXaxis()->GetNbins() + 2; } static Int_t Nbins(HT const &h) { return h.GetXaxis()->GetNbins(); } static bool IsFlowBin(HT const &h, Int_t bin_it) { return nuis::utility::IsFlowBin(*h.GetXaxis(), bin_it); } static void Fill(HT &h, std::array const &v, double w = 1) { h.Fill(v[0], w); } static void Scale(HT &h, NumericT SF, char const *opt = "") { h.Scale(SF, opt); } static double Integral(HT const &h, char const *opt = "") { return h.Integral(opt); } static Int_t NbinsIncludeFlow(std::unique_ptr const &h) { return NbinsIncludeFlow(*h); } static Int_t Nbins(std::unique_ptr const &h) { return Nbins(*h); } static bool IsFlowBin(std::unique_ptr const &h, Int_t bin_it) { return IsFlowBin(*h, bin_it); } static void Fill(std::unique_ptr &h, std::array const &v, double w = 1) { Fill(*h, v, w); } static void Scale(std::unique_ptr &h, NumericT SF, char const *opt = "") { Scale(*h, SF, opt); } static double Integral(std::unique_ptr const &h, char const *opt = "") { return Integral(*h, opt); } }; template struct TH_Helper::NDim == 2>::type> { static size_t const NDim = HType_traits::NDim; using NumericT = typename HType_traits::NumericT; static std::string name() { return HType_traits::name(); } // TH2 *************************************************************** static Int_t NbinsIncludeFlow(HT const &h) { return (h.GetXaxis()->GetNbins() + 2) * (h.GetYaxis()->GetNbins() + 2); } static Int_t Nbins(HT const &h) { return (h.GetXaxis()->GetNbins()) * (h.GetYaxis()->GetNbins()); } static bool IsFlowBin(HT const &h, Int_t xbin_it, Int_t ybin_it) { return nuis::utility::IsFlowBin(*h.GetXaxis(), xbin_it) || nuis::utility::IsFlowBin(*h.GetYaxis(), ybin_it); } static void Fill(HT &h, std::array const &v, double w = 1) { h.Fill(v[0], v[1], w); } static void Scale(HT &h, NumericT SF, char const *opt = "") { h.Scale(SF, opt); } static double Integral(HT const &h, char const *opt = "") { return h.Integral(opt); } static Int_t NbinsIncludeFlow(std::unique_ptr const &h) { return NbinsIncludeFlow(*h); } static Int_t Nbins(std::unique_ptr const &h) { return Nbins(*h); } static bool IsFlowBin(std::unique_ptr const &h, Int_t xbin_it, Int_t ybin_it) { return IsFlowBin(*h, xbin_it, ybin_it); } static void Fill(std::unique_ptr &h, std::array const &v, double w = 1) { Fill(*h, v, w); } static void Scale(std::unique_ptr &h, NumericT SF, char const *opt = "") { Scale(*h, SF, opt); } static double Integral(std::unique_ptr const &h, char const *opt = "") { return Integral(*h, opt); } }; template struct TH_Helper< HT, typename std::enable_if::value>::type> { static size_t const NDim = 2; using NumericT = typename HType_traits::NumericT; static std::string name() { return HType_traits::name(); } // TH2Poly *************************************************************** static Int_t NbinsIncludeFlow(HT const &h) { return h.GetNumberOfBins() + 9; } static Int_t Nbins(HT const &h) { return h.GetNumberOfBins(); } static bool IsFlowBin(HT const &h, Int_t bin_it) { return (bin_it < 0); } static void Fill(HT &h, std::array const &v, double w = 1) { h.Fill(v[0], v[1], w); } static void Scale(HT &h, NumericT SF, char const *opt = "") { bool width = (std::string(opt).find("width") != std::string::npos); size_t nbins = Nbins(h); for (size_t bin_it = 0; bin_it < nbins; ++bin_it) { double bin_area = 1; if (width) { TH2PolyBin *poly_bin = dynamic_cast(h.GetBins()->At(bin_it)); bin_area = poly_bin->GetArea(); } h.SetBinContent(bin_it + 1, h.GetBinContent(bin_it + 1) * (SF / bin_area)); h.SetBinError(bin_it + 1, h.GetBinError(bin_it + 1) * (SF / bin_area)); } } static double Integral(HT &h, char const *opt = "") { bool width = (std::string(opt).find("width") != std::string::npos); size_t nbins = Nbins(h); double integral = 0; for (size_t bin_it = 0; bin_it < nbins; ++bin_it) { double bin_area = 1; if (width) { TH2PolyBin *poly_bin = dynamic_cast(h.GetBins()->At(bin_it)); bin_area = poly_bin->GetArea(); } integral += h.GetBinContent(bin_it + 1) * bin_area; } return integral; } static Int_t NbinsIncludeFlow(std::unique_ptr const &h) { return NbinsIncludeFlow(*h); } static Int_t Nbins(std::unique_ptr const &h) { return Nbins(*h); } static bool IsFlowBin(std::unique_ptr const &h, Int_t bin_it) { return IsFlowBin(*h, bin_it); } static void Fill(std::unique_ptr &h, std::array const &v, double w = 1) { Fill(*h, v, w); } static void Scale(std::unique_ptr &h, NumericT SF, char const *opt = "") { Scale(*h, SF, opt); } static double Integral(std::unique_ptr &h, char const *opt = "") { return Integral(*h, opt); } }; template void Clear(typename std::enable_if::NDim != 0, HT>::type &h) { for (Int_t bin_it = 0; bin_it < TH_Helper::NbinsIncludeFlow(h); ++bin_it) { h.SetBinContent(bin_it, 0); h.SetBinError(bin_it, 0); } } template void Clear( typename std::enable_if::value, HT>::type &h) { h.ClearBinContents(); } template inline std::unique_ptr Clone(HT const &source, bool clear = false, std::string const &clone_name = "") { std::unique_ptr target(dynamic_cast( source.Clone(clone_name.size() ? clone_name.c_str() : ""))); if (!target) { throw failed_to_clone() << "[ERROR]: Failed to clone a " << TH_Helper::name() << "."; } target->SetDirectory(nullptr); if (clear) { Clear(*target); } return target; } template inline std::unique_ptr Clone(std::unique_ptr const &source, bool clear = false, std::string const &clone_name = "") { return Clone(*source, clear, clone_name); } template inline std::unique_ptr GetHistogramFromROOTFile(TFile_ptr &f, std::string const &hname, bool ThrowIfMissing = true) { - f->Get(hname.c_str()); - HT *h = dynamic_cast(f->Get(hname.c_str())); + TObject *obj = f->Get(hname.c_str()); + if(std::string(obj->ClassName()) == "TList"){ + obj = static_cast(obj)->At(0); + } + + HT *h = dynamic_cast(obj); if (!h) { if (ThrowIfMissing) { throw invalid_histogram_name() << "[ERROR]: Failed to get " << TH_Helper::name() << " named " << std::quoted(hname) << " from input file " << std::quoted(f->GetName()); } else { return nullptr; } } std::unique_ptr clone = Clone(*h); return clone; } template inline std::unique_ptr GetHistogramFromROOTFile(std::string const &fname, std::string const &hname, bool ThrowIfMissing = true) { TFile_ptr temp = CheckOpenTFile(fname, "READ"); return GetHistogramFromROOTFile(temp, hname, ThrowIfMissing); } template std::unique_ptr GetHistogram(std::string const &input_descriptor) { std::vector split_descriptor = fhicl::string_parsers::ParseToVect(input_descriptor, ";", true, true); if (split_descriptor.size() == 1) { // assume text throw unimplemented_GetHistogram_method(); } else if (split_descriptor.size() == 2) { return GetHistogramFromROOTFile(split_descriptor[0], split_descriptor[1]); } else { throw invalid_histogram_descriptor() << "[ERROR]: Failed to parse histogram descriptor: " << std::quoted(input_descriptor) << " as an input histogram (Text/ROOT)."; } } struct PolyBinSpecifier { double X, Y; bool UseXAxis; }; constexpr PolyBinSpecifier XPolyBinSpec(double X, double Y) { return PolyBinSpecifier{X + std::numeric_limits::epsilon() * 1E2, Y + std::numeric_limits::epsilon() * 1E2, true}; } constexpr PolyBinSpecifier YPolyBinSpec(double X, double Y) { return PolyBinSpecifier{X + std::numeric_limits::epsilon() * 1E2, Y + std::numeric_limits::epsilon() * 1E2, false}; } std::vector> GetTH2PolySlices( std::unique_ptr &hinp, std::vector> const &BinsSpecifiers); template inline typename std::enable_if::NDim == 1, std::unique_ptr>::type BuildHistFromFHiCL(fhicl::ParameterSet const &ps) { std::array xaxis = ps.get>("xaxis_descriptor"); std::unique_ptr rtn = std::make_unique( ps.get("name", "").c_str(), ps.get("title", "").c_str(), xaxis[0], xaxis[1], xaxis[2]); rtn->SetDirectory(nullptr); return rtn; } template inline typename std::enable_if::NDim == 2, std::unique_ptr>::type BuildHistFromFHiCL(fhicl::ParameterSet const &ps) { std::array xaxis = ps.get>("xaxis_descriptor"); std::array yaxis = ps.get>("yaxis_descriptor"); std::unique_ptr rtn = std::make_unique(ps.get("name", "").c_str(), ps.get("title", "").c_str(), xaxis[0], xaxis[1], xaxis[2], yaxis[0], yaxis[1], yaxis[2]); rtn->SetDirectory(nullptr); return rtn; } static bool const kYSlice = true; static bool const kXSlice = false; void SliceNorm(std::unique_ptr &hist, bool AlongY = kYSlice, char const *opt = ""); } // namespace utility } // namespace nuis diff --git a/src/utility/InteractionChannelUtility.hxx b/src/utility/InteractionChannelUtility.hxx index 00b8a8c..78abab2 100644 --- a/src/utility/InteractionChannelUtility.hxx +++ b/src/utility/InteractionChannelUtility.hxx @@ -1,74 +1,109 @@ #ifndef UTILITY_CHANNELUTILITY_HXX_SEEN #define UTILITY_CHANNELUTILITY_HXX_SEEN #include "event/types.hxx" #include "exception/exception.hxx" #include namespace nuis { namespace utility { NEW_NUIS_EXCEPT(invalid_channel_conversion); #define X(A, B) \ case B: { \ return nuis::event::Channel_t::A; \ } inline event::Channel_t IntToChannel(int mode) { switch (mode) { NUIS_INTERACTION_CHANNEL_LIST default: { throw invalid_channel_conversion() << "[ERROR]: Failed to parse mode integer " << mode << " as a nuis::event::Channel_t."; } } } #undef X #define X(A, B) \ case event::Channel_t::A: { \ return B; \ } inline int ChannelToInt(event::Channel_t mode) { switch (mode) { NUIS_INTERACTION_CHANNEL_LIST default: { throw invalid_channel_conversion() << "[ERROR]: Attempting to convert " "undefined nuis::event::Channel_t to an " "integer."; } } } #undef X inline bool IsNC(event::Channel_t mode) { return abs(ChannelToInt(mode)) > 30; } inline bool IsCC(event::Channel_t mode) { return !IsNC(mode); } inline bool IsNu(event::Channel_t mode) { return ChannelToInt(mode) > 0; } inline bool IsNub(event::Channel_t mode) { return !IsNu(mode); } -inline bool IsCoh(event::Channel_t mode) { - return ((mode == event::Channel_t::kCCCohPi) || - (mode == event::Channel_t::kNCCohPi) || - (mode == event::Channel_t::kCCCohPi_nub) || - (mode == event::Channel_t::kNCCohPi_nub)); -} inline bool IsQE(event::Channel_t mode) { return ((mode == event::Channel_t::kCCQE) || (mode == event::Channel_t::kNCELP) || (mode == event::Channel_t::kNCELN) || (mode == event::Channel_t::kCCQE_nub) || (mode == event::Channel_t::kNCELP_nub) || (mode == event::Channel_t::kNCELN_nub)); } - +inline bool Is2p2h(event::Channel_t mode) { + return ((mode == event::Channel_t::kCC2p2h) || + (mode == event::Channel_t::kNC2p2h) || + (mode == event::Channel_t::kCC2p2h_nub) || + (mode == event::Channel_t::kNC2p2h_nub)); +} +inline bool IsQEL(event::Channel_t mode) { + return (IsQE(mode) || Is2p2h(mode)); +} +inline bool IsNucleonSPP(event::Channel_t mode) { + return ((mode == event::Channel_t::kCCSPP_PPip) || + (mode == event::Channel_t::kCCSPP_PPi0) || + (mode == event::Channel_t::kCCSPP_NPip) || + (mode == event::Channel_t::kNCSPP_NPi0) || + (mode == event::Channel_t::kNCSPP_PPi0) || + (mode == event::Channel_t::kNCSPP_PPim) || + (mode == event::Channel_t::kCCSPP_NPim_nub) || + (mode == event::Channel_t::kCCSPP_NPi0_nub) || + (mode == event::Channel_t::kCCSPP_PPim_nub) || + (mode == event::Channel_t::kNCSPP_NPi0_nub) || + (mode == event::Channel_t::kNCSPP_PPi0_nub) || + (mode == event::Channel_t::kNCSPP_PPim_nub) || + (mode == event::Channel_t::kNCSPP_NPip_nub)); +} +inline bool IsMultiPi(event::Channel_t mode) { + return ((mode == event::Channel_t::kCCTransitionMPi) || + (mode == event::Channel_t::kNCTransitionMPi) || + (mode == event::Channel_t::kCCTransitionMPi_nub) || + (mode == event::Channel_t::kNCTransitionMPi_nub)); +} +inline bool IsDIS(event::Channel_t mode) { + return ((mode == event::Channel_t::kCCDIS) || + (mode == event::Channel_t::kNCDIS) || + (mode == event::Channel_t::kCCDIS_nub) || + (mode == event::Channel_t::kNCDIS_nub)); +} +inline bool IsCoh(event::Channel_t mode) { + return ((mode == event::Channel_t::kCCCohPi) || + (mode == event::Channel_t::kNCCohPi) || + (mode == event::Channel_t::kCCCohPi_nub) || + (mode == event::Channel_t::kNCCohPi_nub)); +} } // namespace utility } // namespace nuis #endif diff --git a/src/utility/KinematicUtility.cxx b/src/utility/KinematicUtility.cxx index be2a305..369a3c7 100644 --- a/src/utility/KinematicUtility.cxx +++ b/src/utility/KinematicUtility.cxx @@ -1,228 +1,290 @@ #include "utility/KinematicUtility.hxx" #include "event/FullEvent.hxx" #include "event/Particle.hxx" #include "utility/FullEventUtility.hxx" #include "utility/PDGCodeUtility.hxx" using namespace nuis::event; namespace nuis { namespace utility { double GetNeutrinoEQERec(FullEvent const &fev, double SeparationEnergy_MeV) { Particle const &neutrino = GetHMISNeutralLepton(fev); if (!neutrino) { throw Particle::invalid_particle() << "[ERROR]: In GetNeutrinoEQERec, expected to be able to get IS " "neutral lepton, but found none:" << "\n" << fev.to_string(); } - Particle const &charged_lepton = GetHMFSChargedLepton(fev); - if (!charged_lepton) { + Particle const &lepton = GetHMFSLepton(fev); + if (!lepton) { throw Particle::invalid_particle() << "[ERROR]: In GetNeutrinoEQERec, expected to be able to get FS " "charged lepton, but found none:" << "\n" << fev.to_string(); } double mass_ISN = (IsMatter(neutrino.pdg) ? (pdgmasses::kNeutronMass_MeV - SeparationEnergy_MeV) : (pdgmasses::kProtonMass_MeV - SeparationEnergy_MeV)) / 1000.0; double mass_FSN = (IsMatter(neutrino.pdg) ? pdgmasses::kProtonMass_MeV : pdgmasses::kNeutronMass_MeV) / 1000.0; - double el_GeV = charged_lepton.E() / 1000.0; - double pl_GeV = charged_lepton.P() / 1000.0; // momentum of lepton - double ml_GeV = charged_lepton.M() / 1000.0; // lepton mass + double el_GeV = lepton.E() / 1000.0; + double pl_GeV = lepton.P() / 1000.0; // momentum of lepton + double ml_GeV = lepton.M() / 1000.0; // lepton mass - double Theta_nu_mu = neutrino.P3().Angle(charged_lepton.P3()); + double Theta_nu_mu = neutrino.P3().Angle(lepton.P3()); return (2.0 * mass_ISN * el_GeV - ml_GeV * ml_GeV + mass_FSN * mass_FSN - mass_ISN * mass_ISN) / (2.0 * (mass_ISN - el_GeV + pl_GeV * cos(Theta_nu_mu))); } double GetNeutrinoQ2QERec(FullEvent const &fev, double SeparationEnergy_MeV) { Particle const &neutrino = GetHMISNeutralLepton(fev); if (!neutrino) { throw Particle::invalid_particle() << "[ERROR]: In GetNeutrinoQ2QERec, expected to be able to get IS " "neutral lepton, but found none:" << "\n" << fev.to_string(); } - Particle const &charged_lepton = GetHMFSChargedLepton(fev); - if (!charged_lepton) { + Particle const &lepton = GetHMFSLepton(fev); + if (!lepton) { throw Particle::invalid_particle() << "[ERROR]: In GetNeutrinoQ2QERec, expected to be able to get FS " "charged lepton, but found none:" << "\n" << fev.to_string(); } - double el_GeV = charged_lepton.E() / 1000.0; - double pl_GeV = charged_lepton.P() / 1000.0; // momentum of lepton - double ml_GeV = charged_lepton.M() / 1000.0; // lepton mass + double el_GeV = lepton.E() / 1000.0; + double pl_GeV = lepton.P() / 1000.0; // momentum of lepton + double ml_GeV = lepton.M() / 1000.0; // lepton mass - double Theta_nu_mu = neutrino.P3().Angle(charged_lepton.P3()); + double Theta_nu_mu = neutrino.P3().Angle(lepton.P3()); return -ml_GeV * ml_GeV + 2.0 * GetNeutrinoEQERec(fev, SeparationEnergy_MeV) * (el_GeV - pl_GeV * cos(Theta_nu_mu)); } +double GetNeutrinoWRec(event::FullEvent const &fev) { + Particle const &neutrino = GetHMISNeutralLepton(fev); + if (!neutrino) { + throw Particle::invalid_particle() + << "[ERROR]: In GetNeutrinoWRec, expected to be able to get IS " + "neutral lepton, but found none:" + << "\n" + << fev.to_string(); + } + + Particle const &lepton = GetHMFSLepton(fev); + if (!lepton) { + throw Particle::invalid_particle() + << "[ERROR]: In GetNeutrinoWRec, expected to be able to get FS " + "charged lepton, but found none:" + << "\n" + << fev.to_string(); + } + + double E_mu = lepton.E(); + double p_mu = lepton.P(); + double m_mu = sqrt(E_mu * E_mu - p_mu * p_mu); + double th_nu_mu = neutrino.P3().Angle(lepton.P3()); + + // The factor of 1000 is necessary for downstream functions + double m_p = pdgmasses::kProtonMass_MeV; + + // MINERvA cut on W_exp which is tuned to W_true; so use true Enu from + // generators + double E_nu = neutrino.E(); + + double w_rec = sqrt(m_p * m_p + m_mu * m_mu - 2 * m_p * E_mu + + 2 * E_nu * (m_p - E_mu + p_mu * cos(th_nu_mu))); + + return w_rec; +} + +double GetEAvailProxy(event::FullEvent const &fev) { + Particle ISNu = GetHMISNeutralLepton(fev); + if (!ISNu) { + return 0; + } + + Particle FSLep = GetHMFSLepton(fev); + if (!FSLep) { + return 0; + } + + TLorentzVector EMTransfer_lep = ISNu.P4 - FSLep.P4; + double eav = EMTransfer_lep.E(); + + for (auto const &neutron : GetFSNeutrons(fev)) { + eav -= neutron.KE(); + } + + for (auto const &pion : GetFSPions(fev)) { + eav -= pion.M(); + } + + return eav; +} + TVector3 GetVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal) { TVector3 pnUnit = planarNormal.Unit(); double inpProjectPN = inp.Dot(pnUnit); TVector3 InPlanarInput = inp - (inpProjectPN * pnUnit); return InPlanarInput; } TVector3 GetUnitVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal) { return GetVectorInTPlane(inp, planarNormal).Unit(); } double GetDeltaPhiT(TVector3 const &V_lepton, TVector3 const &V_other, TVector3 const &Normal, bool PiMinus) { TVector3 V_lepton_T = GetUnitVectorInTPlane(V_lepton, Normal); TVector3 V_other_T = GetUnitVectorInTPlane(V_other, Normal); return PiMinus ? acos(V_lepton_T.Dot(V_other_T)) : (M_PI - acos(V_lepton_T.Dot(V_other_T))); } TVector3 GetDeltaPT(TVector3 const &V_lepton, TVector3 const &V_other, TVector3 const &Normal) { TVector3 V_lepton_T = GetVectorInTPlane(V_lepton, Normal); TVector3 V_other_T = GetVectorInTPlane(V_other, Normal); return V_lepton_T + V_other_T; } double GetDeltaAlphaT(TVector3 const &V_lepton, TVector3 const &V_other, TVector3 const &Normal, bool PiMinus) { TVector3 DeltaPT = GetDeltaPT(V_lepton, V_other, Normal); return GetDeltaPhiT(V_lepton, DeltaPT, Normal, PiMinus); } TVector3 GetDeltaPT_CC0PiN(event::FullEvent const &fev, event::PDG_t islep_pdg, event::PDG_t fslep_pdg, event::PDG_t fsnuc_pdg) { if (fslep_pdg == nuis::utility::pdgcodes::kDefault) { // Assume neutrino CC fslep_pdg = islep_pdg > 0 ? islep_pdg - 1 : islep_pdg + 1; } if (fsnuc_pdg == nuis::utility::pdgcodes::kDefault) { fsnuc_pdg = islep_pdg > 0 ? nuis::utility::pdgcodes::kProton : nuis::utility::pdgcodes::kNeutron; } Particle ISSLep = GetHMISParticle(fev, {islep_pdg}); if (!ISSLep) { return {-std::numeric_limits::max(), -std::numeric_limits::max(), -std::numeric_limits::max()}; } Particle FSLep = GetHMFSParticle(fev, {fslep_pdg}); if (!FSLep) { return {-std::numeric_limits::max(), -std::numeric_limits::max(), -std::numeric_limits::max()}; } Particle FSNuc = GetHMFSParticle(fev, {fsnuc_pdg}); if (!FSNuc) { return {-std::numeric_limits::max(), -std::numeric_limits::max(), -std::numeric_limits::max()}; } return GetDeltaPT(FSLep.P3(), FSNuc.P3(), ISSLep.P3()); } double GetDeltaPhiT_CC0PiN(event::FullEvent const &fev, event::PDG_t islep_pdg, event::PDG_t fslep_pdg, event::PDG_t fsnuc_pdg) { if (fslep_pdg == nuis::utility::pdgcodes::kDefault) { // Assume neutrino CC fslep_pdg = islep_pdg > 0 ? islep_pdg - 1 : islep_pdg + 1; } if (fsnuc_pdg == nuis::utility::pdgcodes::kDefault) { fsnuc_pdg = islep_pdg > 0 ? nuis::utility::pdgcodes::kProton : nuis::utility::pdgcodes::kNeutron; } Particle ISSLep = GetHMISParticle(fev, {islep_pdg}); if (!ISSLep) { return -std::numeric_limits::max(); } Particle FSLep = GetHMFSParticle(fev, {fslep_pdg}); if (!FSLep) { return -std::numeric_limits::max(); } Particle FSNuc = GetHMFSParticle(fev, {fsnuc_pdg}); if (!FSNuc) { return -std::numeric_limits::max(); } return GetDeltaPhiT(FSLep.P3(), FSNuc.P3(), ISSLep.P3()); } double GetDeltaAlphaT_CC0PiN(event::FullEvent const &fev, event::PDG_t islep_pdg, event::PDG_t fslep_pdg, event::PDG_t fsnuc_pdg) { if (fslep_pdg == nuis::utility::pdgcodes::kDefault) { // Assume neutrino CC fslep_pdg = islep_pdg > 0 ? islep_pdg - 1 : islep_pdg + 1; } if (fsnuc_pdg == nuis::utility::pdgcodes::kDefault) { fsnuc_pdg = islep_pdg > 0 ? nuis::utility::pdgcodes::kProton : nuis::utility::pdgcodes::kNeutron; } Particle ISSLep = GetHMISParticle(fev, {islep_pdg}); if (!ISSLep) { return -std::numeric_limits::max(); } Particle FSLep = GetHMFSParticle(fev, {fslep_pdg}); if (!FSLep) { return -std::numeric_limits::max(); } Particle FSNuc = GetHMFSParticle(fev, {fsnuc_pdg}); if (!FSNuc) { return -std::numeric_limits::max(); } return GetDeltaAlphaT(FSLep.P3(), FSNuc.P3(), ISSLep.P3()); } TLorentzVector GetEnergyMomentumTransfer(event::FullEvent const &fev) { Particle const &neutrino = GetHMISNeutralLepton(fev); if (!neutrino) { throw Particle::invalid_particle() << "[ERROR]: In GetEnergyMomentumTransfer, expected to be able to get " "IS neutral lepton, but found none: \n" << fev.to_string(); } Particle const &fs_lepton = GetHMFSLepton(fev); if (!fs_lepton) { throw Particle::invalid_particle() << "[ERROR]: In GetEnergyMomentumTransfer, expected to be able to get " "FS lepton, but found none: \n" << fev.to_string(); } return (neutrino.P4 - fs_lepton.P4); } } // namespace utility } // namespace nuis diff --git a/src/utility/KinematicUtility.hxx b/src/utility/KinematicUtility.hxx index 52c0536..d2b445a 100644 --- a/src/utility/KinematicUtility.hxx +++ b/src/utility/KinematicUtility.hxx @@ -1,86 +1,95 @@ // Copyright 2018 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 UTILITY_KINEMATICUTILITY_HXX_SEEN #define UTILITY_KINEMATICUTILITY_HXX_SEEN #include "utility/PDGCodeUtility.hxx" #include "TLorentzVector.h" #include "TVector3.h" #include #include namespace nuis { namespace event { class FullEvent; } } // namespace nuis namespace nuis { namespace utility { double GetNeutrinoEQERec(event::FullEvent const &fev, double SeparationEnergy_MeV); double GetNeutrinoQ2QERec(event::FullEvent const &fev, double SeparationEnergy_MeV); +double GetNeutrinoWRec(event::FullEvent const &fev); +double GetEAvailProxy(event::FullEvent const &fev); TVector3 GetVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal); TVector3 GetUnitVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal); double GetDeltaPhiT(TVector3 const &V_lepton, TVector3 const &V_other, TVector3 const &Normal, bool PiMinus = false); TVector3 GetDeltaPT(TVector3 const &V_lepton, TVector3 const &V_other, TVector3 const &Normal); double GetDeltaAlphaT(TVector3 const &V_lepton, TVector3 const &V_other, TVector3 const &Normal, bool PiMinus = false); TVector3 GetDeltaPT_CC0PiN(event::FullEvent const &fev, event::PDG_t islep_pdg = nuis::utility::pdgcodes::kNuMu, event::PDG_t fslep_pdg = nuis::utility::pdgcodes::kDefault, event::PDG_t fsnuc_pdg = nuis::utility::pdgcodes::kDefault); double GetDeltaPhiT_CC0PiN(event::FullEvent const &fev, event::PDG_t islep_pdg = nuis::utility::pdgcodes::kNuMu, event::PDG_t fslep_pdg = nuis::utility::pdgcodes::kDefault, event::PDG_t fsnuc_pdg = nuis::utility::pdgcodes::kDefault); double GetDeltaAlphaT_CC0PiN( event::FullEvent const &fev, event::PDG_t islep_pdg = nuis::utility::pdgcodes::kNuMu, event::PDG_t fslep_pdg = nuis::utility::pdgcodes::kDefault, event::PDG_t fsnuc_pdg = nuis::utility::pdgcodes::kDefault); TLorentzVector GetEnergyMomentumTransfer(event::FullEvent const &fev); -struct ENuRange : public std::pair { +struct KinematicRange : public std::pair { using std::pair::pair; using std::pair::operator=; - ENuRange() + KinematicRange() : std::pair(0, std::numeric_limits::max()) {} - bool IsInRange(double enu) { return (enu > second) || (enu < first); } + // Not sure why something like this isn't being inherited... + KinematicRange(std::pair const &o) + : std::pair(std::min(o.first, o.second), + std::max(o.first, o.second)) {} + + bool IsInRange(double val) const { + return !((val > second) || (val < first)); + } }; } // namespace utility } // namespace nuis #endif diff --git a/src/utility/FileSystemUtility.hxx b/src/utility/PhaseSpaceRestriction.cxx similarity index 74% copy from src/utility/FileSystemUtility.hxx copy to src/utility/PhaseSpaceRestriction.cxx index eb7f7a0..fd3571b 100644 --- a/src/utility/FileSystemUtility.hxx +++ b/src/utility/PhaseSpaceRestriction.cxx @@ -1,37 +1,37 @@ // Copyright 2018 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 UTILITY_FILESYSTEMUTILITY_HXX_SEEN -#define UTILITY_FILESYSTEMUTILITY_HXX_SEEN +#include "utility/PhaseSpaceRestriction.hxx" -#include -#include +#include "event/Particle.hxx" namespace nuis { namespace utility { - -std::string EnsureTrailingSlash(std::string); - -std::vector GetMatchingFiles(std::string, std::string const &); - -std::string GetDataDir(); - +bool SimpleParticlePhaseSpaceRestriction::Inside( + nuis::event::Particle const &part) const { + double ep = IsERange ? part.E() : part.P(); + double ct = part.CosTheta(); + + if (!EMomRange.IsInRange(ep) || !CosThetaRange.IsInRange(ct)) { + return false; + } + return true; +} } // namespace utility } // namespace nuis -#endif diff --git a/src/utility/PhaseSpaceRestriction.hxx b/src/utility/PhaseSpaceRestriction.hxx new file mode 100644 index 0000000..06893fa --- /dev/null +++ b/src/utility/PhaseSpaceRestriction.hxx @@ -0,0 +1,147 @@ +// Copyright 2018 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 "utility/KinematicUtility.hxx" + +#include +#include + +#include + +#pragma once +namespace nuis { +namespace event { +class Particle; +} // namespace event +} // namespace nuis + +namespace nuis { +namespace utility { + +class IParticlePhaseSpaceRestriction { +public: + virtual bool Inside(nuis::event::Particle const &) const = 0; +}; + +class FullPhaseSpace : public IParticlePhaseSpaceRestriction { +public: + bool Inside(nuis::event::Particle const &) const { return true; } +}; + +class SimpleParticlePhaseSpaceRestriction + : public IParticlePhaseSpaceRestriction { + KinematicRange EMomRange; + KinematicRange CosThetaRange; + bool IsERange; + +public: + ///\brief Simple E/P and angle phase space restriction constructor + /// + ///\detail Best to use one of the named static methods to avoid E/P or + /// theta/costheta confusion. + SimpleParticlePhaseSpaceRestriction( + std::pair eprange = {0, + std::numeric_limits::max()}, + bool iserange = true, std::pair costhetarange = {0, M_PI}) + : EMomRange(eprange), CosThetaRange(costhetarange), IsERange(iserange) {} + + bool Inside(nuis::event::Particle const &) const; + + static SimpleParticlePhaseSpaceRestriction + Energy(double min = 0, double max = std::numeric_limits::max()) { + return SimpleParticlePhaseSpaceRestriction({min, max}); + } + + static SimpleParticlePhaseSpaceRestriction Theta_deg(double min = 0, + double max = 180) { + return SimpleParticlePhaseSpaceRestriction( + {0, std::numeric_limits::max()}, true, + {cos(min * (M_PI / 180.0)), cos(max * (M_PI / 180.0))}); + } + + static SimpleParticlePhaseSpaceRestriction Theta_rad(double min = 0, + double max = M_PI) { + return SimpleParticlePhaseSpaceRestriction( + {0, std::numeric_limits::max()}, true, {cos(min), cos(max)}); + } + + static SimpleParticlePhaseSpaceRestriction CosTheta(double min = -1, + double max = 1) { + return SimpleParticlePhaseSpaceRestriction( + {0, std::numeric_limits::max()}, true, {min, max}); + } + + static SimpleParticlePhaseSpaceRestriction + Momentum(double min = 0, double max = std::numeric_limits::max()) { + return SimpleParticlePhaseSpaceRestriction({min, max}, false); + } + + static SimpleParticlePhaseSpaceRestriction EnergyTheta_deg( + std::pair erange = {0, + std::numeric_limits::max()}, + std::pair thetarange = {0, 180}) { + return SimpleParticlePhaseSpaceRestriction( + erange, true, + {cos(thetarange.first * (M_PI / 180.0)), + cos(thetarange.second * (M_PI / 180.0))}); + } + + static SimpleParticlePhaseSpaceRestriction EnergyTheta_rad( + std::pair erange = {0, + std::numeric_limits::max()}, + std::pair thetarange = {0, M_PI}) { + return SimpleParticlePhaseSpaceRestriction( + erange, true, {cos(thetarange.first), cos(thetarange.second)}); + } + + static SimpleParticlePhaseSpaceRestriction EnergyCosTheta( + std::pair erange = {0, + std::numeric_limits::max()}, + std::pair costhetarange = {-1, 1}) { + return SimpleParticlePhaseSpaceRestriction(erange, true, costhetarange); + } + + static SimpleParticlePhaseSpaceRestriction MomentumTheta_deg( + std::pair prange = {0, + std::numeric_limits::max()}, + std::pair thetarange = {0, 180}) { + return SimpleParticlePhaseSpaceRestriction( + prange, false, + {cos(thetarange.first * (M_PI / 180.0)), + cos(thetarange.second * (M_PI / 180.0))}); + } + + static SimpleParticlePhaseSpaceRestriction MomentumTheta_rad( + std::pair prange = {0, + std::numeric_limits::max()}, + std::pair thetarange = {0, M_PI}) { + return SimpleParticlePhaseSpaceRestriction( + prange, false, {cos(thetarange.first), cos(thetarange.second)}); + } + + static SimpleParticlePhaseSpaceRestriction MomentumCosTheta( + std::pair prange = {0, + std::numeric_limits::max()}, + std::pair costhetarange = {-1, 1}) { + return SimpleParticlePhaseSpaceRestriction(prange, false, costhetarange); + } +}; + +} // namespace utility +} // namespace nuis diff --git a/src/utility/StringUtility.cxx b/src/utility/StringUtility.cxx index c0ceec6..5b4161d 100644 --- a/src/utility/StringUtility.cxx +++ b/src/utility/StringUtility.cxx @@ -1,71 +1,165 @@ +#include "exception/exception.hxx" + #include "utility/StringUtility.hxx" +#include "string_parsers/from_string.hxx" #include "string_parsers/utility.hxx" #include #include #include namespace nuis { namespace utility { +NEW_NUIS_EXCEPT(unexpected_empty_string); + +std::string EnsureTrailingSlash(std::string str) { + if (!str.size()) { + throw unexpected_empty_string(); + } + if (str.back() != '/') { + return str + '/'; + } + return str; +} + +std::string parseCode(std::regex_constants::error_type etype) { + switch (etype) { + case std::regex_constants::error_collate: + return "error_collate: invalid collating element request"; + case std::regex_constants::error_ctype: + return "error_ctype: invalid character class"; + case std::regex_constants::error_escape: + return "error_escape: invalid escape character or trailing escape"; + case std::regex_constants::error_backref: + return "error_backref: invalid back reference"; + case std::regex_constants::error_brack: + return "error_brack: mismatched bracket([ or ])"; + case std::regex_constants::error_paren: + return "error_paren: mismatched parentheses(( or ))"; + case std::regex_constants::error_brace: + return "error_brace: mismatched brace({ or })"; + case std::regex_constants::error_badbrace: + return "error_badbrace: invalid range inside a { }"; + case std::regex_constants::error_range: + return "erro_range: invalid character range(e.g., [z-a])"; + case std::regex_constants::error_space: + return "error_space: insufficient memory to handle this regular expression"; + case std::regex_constants::error_badrepeat: + return "error_badrepeat: a repetition character (*, ?, +, or {) was not " + "preceded by a valid regular expression"; + case std::regex_constants::error_complexity: + return "error_complexity: the requested match is too complex"; + case std::regex_constants::error_stack: + return "error_stack: insufficient memory to evaluate a match"; + default: + return ""; + } +} + +std::string DeGlobPattern(std::string const &pattern) { + std::stringstream ss(""); + size_t next_asterisk = pattern.find_first_of('*'); + size_t next_to_add = 0; + bool modified = false; + while (next_asterisk != std::string::npos) { + if ((pattern[next_asterisk - 1] != + ']') && // Try to allow valid uses of an asterisk without a preceding. + (pattern[next_asterisk - 1] != '.')) { + modified = true; + // Add a . + ss << pattern.substr(next_to_add, next_asterisk - next_to_add) << ".*"; + next_to_add = next_asterisk + 1; + if (next_to_add >= pattern.size()) { + next_to_add = std::string::npos; + } + } + next_asterisk = pattern.find_first_of('*', next_asterisk + 1); + } + + if (next_to_add != std::string::npos) { + ss << pattern.substr(next_to_add); + } + + if (modified) { + std::cout << "[INFO]: DeGlobified input pattern: " << pattern + << " to std::regex_friendly: " << ss.str() << std::endl; + } + + return ss.str(); +} + // #define DEBUG_INDENT_APPLY_WIDTH std::string indent_apply_width(std::string input, size_t indent, size_t width) { std::stringstream output_stream; std::string indent_string(indent, ' '); if (width) { width--; } if (width <= indent) { return input; } #ifdef DEBUG_INDENT_APPLY_WIDTH std::cout << "[IDEBUG]: Indenting and forcing width on: " << std::quoted(input) << ", width = " << width << ", indent = " << indent << std::endl; #endif while (input.length()) { #ifdef DEBUG_INDENT_APPLY_WIDTH std::cout << "[DEBUG]\tRemaining string: " << std::quoted(input) << ", width = " << width << ", indent = " << indent << std::endl; #endif if (input.length() < width - indent) { #ifdef DEBUG_INDENT_APPLY_WIDTH std::cout << "[DEBUG]\tRemaining string is less than allowed width " << std::quoted(input) << std::endl; #endif output_stream << indent_string << input << std::endl; break; } size_t last_ws = input.find_last_of(" ", width - indent); size_t first_nl = input.find_first_of("\n"); #ifdef DEBUG_INDENT_APPLY_WIDTH std::cout << "[DEBUG]: last_ws = " << last_ws << ", first_nl = " << first_nl << std::endl; #endif size_t next_split = std::min(last_ws, first_nl); #ifdef DEBUG_INDENT_APPLY_WIDTH std::cout << "[DEBUG]: Found split at " << next_split << std::endl; #endif if (next_split == std::string::npos) { next_split = std::min(width - indent, input.length()); } std::string next_line = input.substr(0, next_split); fhicl::string_parsers::trim(next_line); #ifdef DEBUG_INDENT_APPLY_WIDTH std::cout << "[DEBUG]: Streaming " << std::quoted(next_line) << std::endl; #endif output_stream << indent_string << next_line << std::endl; input = input.substr(next_split + 1); } return output_stream.str(); } + +std::vector split(std::string const &str, + std::string const &delim) { + return fhicl::string_parsers::ParseToVect(str, delim, false, + true); +} + +std::vector split(std::string const &str, + std::vector const &delims) { + return fhicl::string_parsers::ParseToVect(str, delims, false, + true); +} + } // namespace utility } // namespace nuis diff --git a/src/utility/StringUtility.hxx b/src/utility/StringUtility.hxx index 541393e..7249c6e 100644 --- a/src/utility/StringUtility.hxx +++ b/src/utility/StringUtility.hxx @@ -1,34 +1,43 @@ // Copyright 2018 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 UTILITY_STRINGUTILITY_HXX_SEEN -#define UTILITY_STRINGUTILITY_HXX_SEEN +#pragma once #include "utility/TerminalUtility.hxx" #include +#include +#include namespace nuis { namespace utility { -std::string indent_apply_width(std::string, - size_t indent = 0, size_t width = GetWindowWidth()); -} -} // namespace nuis -#endif +std::string EnsureTrailingSlash(std::string str); +std::string parseCode(std::regex_constants::error_type etype); +std::string DeGlobPattern(std::string const &pattern); + +std::string indent_apply_width(std::string, size_t indent = 0, + size_t width = GetWindowWidth()); + +std::vector split(std::string const &str, + std::string const &delim); +std::vector split(std::string const &str, + std::vector const &delims); +} // namespace utility +} // namespace nuis diff --git a/src/utility/experimental/CMakeLists.txt b/src/utility/experimental/CMakeLists.txt new file mode 100644 index 0000000..248876f --- /dev/null +++ b/src/utility/experimental/CMakeLists.txt @@ -0,0 +1,16 @@ +set(utility_experimental_implementation_files + T2KUtility.cxx + MINERvAUtility.cxx + ) + +set(utility_experimental_header_files + T2KUtility.hxx + MINERvAUtility.hxx + ) + +add_library(nuis_utility_experimental SHARED ${utility_experimental_implementation_files}) +target_link_libraries(nuis_utility_experimental) + +install(TARGETS nuis_utility_experimental DESTINATION lib) + +install(FILES ${utility_experimental_header_files} DESTINATION include/utility/experimental) diff --git a/src/utility/experimental/MINERvAUtility.cxx b/src/utility/experimental/MINERvAUtility.cxx new file mode 100644 index 0000000..b557426 --- /dev/null +++ b/src/utility/experimental/MINERvAUtility.cxx @@ -0,0 +1,192 @@ +#include "utility/experimental/MINERvAUtility.hxx" + +#include "event/FullEvent.hxx" + +#include "utility/EventTopologyUtility.hxx" +#include "utility/FullEventUtility.hxx" +#include "utility/KinematicUtility.hxx" + +using namespace nuis::event; + +namespace nuis { +namespace utility { +namespace mnv { + +SimpleParticlePhaseSpaceRestriction CC0PiNProt_ProtonPS = + SimpleParticlePhaseSpaceRestriction::MomentumTheta_deg({450, 1200}, {0, 70}); + +bool IsCCIncLowRecoil(event::FullEvent const &ev) { + if (!nuis::utility::IsCCInc(ev)) { + return false; + } + + static SimpleParticlePhaseSpaceRestriction muon_ps = + SimpleParticlePhaseSpaceRestriction::EnergyTheta_deg( + {1500, std::numeric_limits::max()}, {0, 20}); + + if (GetNParticlesInPhaseSpace(ev, {pdgcodes::kMu}, muon_ps) != 1) { + return false; + } + + return true; +} + +std::pair IsCC0Pi_NumProtons(event::FullEvent const &ev) { + if (!nuis::utility::IsCC0Pi(ev)) { + return {false, 0}; + } + + static SimpleParticlePhaseSpaceRestriction proton_ps = + SimpleParticlePhaseSpaceRestriction::Momentum(450); + + return {true, GetNParticlesInPhaseSpace(ev, {pdgcodes::kProton}, proton_ps)}; +} + +bool IsCC0PiNp(event::FullEvent const &ev) { + auto res = IsCC0Pi_NumProtons(ev); + return res.first && (res.second >= 1); +} + +bool IsCC0PiNp_STV(event::FullEvent const &ev) { + if (!nuis::utility::IsCC0Pi(ev)) { + return false; + } + + static SimpleParticlePhaseSpaceRestriction muon_ps = + SimpleParticlePhaseSpaceRestriction::EnergyTheta_deg({1500, 10E3}, + {0, 20}); + + if (GetNParticlesInPhaseSpace(ev, {pdgcodes::kMu}, muon_ps) != 1) { + return false; + } + + if (GetNParticlesInPhaseSpace(ev, {pdgcodes::kProton}, CC0PiNProt_ProtonPS) == + 0) { + return false; + } + + return true; +} + +bool IsCC1Pi0_2016(event::FullEvent const &ev) { + if (!nuis::utility::IsCC1Pi(ev, {pdgcodes::kPi0})) { + return false; + } + return true; +} +bool IsCC1CPi_2017(event::FullEvent const &ev) { + if (!nuis::utility::IsCC1Pi(ev, {pdgcodes::kPiPlus, pdgcodes::kPiMinus})) { + return false; + } + + if (nuis::utility::GetNeutrinoWRec(ev) > 1400) { + return false; + } + + return true; +} + +TVector3 GetDeltaPT_CC0PiN_mnv(event::FullEvent const &fev) { + + Particle ISSLep = GetHMISParticle(fev, {pdgcodes::kNuMu}); + if (!ISSLep) { + return {-std::numeric_limits::max(), + -std::numeric_limits::max(), + -std::numeric_limits::max()}; + } + Particle FSLep = GetHMFSParticle(fev, {pdgcodes::kMu}); + if (!FSLep) { + return {-std::numeric_limits::max(), + -std::numeric_limits::max(), + -std::numeric_limits::max()}; + } + + Particle FSNuc = GetHMFSProtonInPhaseSpace(fev, CC0PiNProt_ProtonPS); + if (!FSNuc) { + return {-std::numeric_limits::max(), + -std::numeric_limits::max(), + -std::numeric_limits::max()}; + } + + return GetDeltaPT(FSLep.P3(), FSNuc.P3(), ISSLep.P3()); +} +double GetDeltaPhiT_CC0PiN_mnv(event::FullEvent const &fev) { + + Particle ISSLep = GetHMISParticle(fev, {pdgcodes::kNuMu}); + if (!ISSLep) { + return -std::numeric_limits::max(); + } + Particle FSLep = GetHMFSParticle(fev, {pdgcodes::kMu}); + if (!FSLep) { + return -std::numeric_limits::max(); + } + + Particle FSNuc = GetHMFSProtonInPhaseSpace(fev, CC0PiNProt_ProtonPS); + if (!FSNuc) { + return -std::numeric_limits::max(); + } + + return GetDeltaPhiT(FSLep.P3(), FSNuc.P3(), ISSLep.P3()); +} +double GetDeltaAlphaT_CC0PiN_mnv(event::FullEvent const &fev) { + + Particle ISSLep = GetHMISParticle(fev, {pdgcodes::kNuMu}); + if (!ISSLep) { + return -std::numeric_limits::max(); + } + Particle FSLep = GetHMFSParticle(fev, {pdgcodes::kMu}); + if (!FSLep) { + return -std::numeric_limits::max(); + } + + Particle FSNuc = GetHMFSProtonInPhaseSpace(fev, CC0PiNProt_ProtonPS); + if (!FSNuc) { + return -std::numeric_limits::max(); + } + + return GetDeltaAlphaT(FSLep.P3(), FSNuc.P3(), ISSLep.P3()); +} + +double GetNeutronMomentumReco_CC0PiN_mnv(event::FullEvent const &fev) { + + Particle ISSLep = GetHMISParticle(fev, {pdgcodes::kNuMu}); + if (!ISSLep) { + return -std::numeric_limits::max(); + } + Particle FSLep = GetHMFSParticle(fev, {pdgcodes::kMu}); + if (!FSLep) { + return -std::numeric_limits::max(); + } + + Particle FSNuc = GetHMFSProtonInPhaseSpace(fev, CC0PiNProt_ProtonPS); + if (!FSNuc) { + return -std::numeric_limits::max(); + } + + double const el = FSLep.E(); + double const eh = FSNuc.E(); + + TVector3 dpt = GetDeltaPT(FSLep.P3(), FSNuc.P3(), ISSLep.P3()); + double dptMag = dpt.Mag(); + + double ma = 6 * pdgmasses::kNeutronMass_MeV + 6 * pdgmasses::kProtonMass_MeV - + 92.16; // target mass (E is from PhysRevC.95.065501) + double map = ma - pdgmasses::kNeutronMass_MeV + 27.13; // remnant mass + + double pmul = FSLep.P3().Dot(ISSLep.P3().Unit()); + double phl = FSNuc.P3().Dot(ISSLep.P3().Unit()); + + double R = ma + pmul + phl - el - eh; + + double dpl = 0.5 * R - (map * map + dptMag * dptMag) / (2 * R); + // double dpl = ((R*R)-(dptMag*dptMag)-(map*map))/(2*R); // as in in + // PhysRevC.95.065501 - gives same result + + double pn_reco = sqrt((dptMag * dptMag) + (dpl * dpl)); + + return pn_reco; +} + +} // namespace mnv +} // namespace utility +} // namespace nuis diff --git a/src/utility/experimental/MINERvAUtility.hxx b/src/utility/experimental/MINERvAUtility.hxx new file mode 100644 index 0000000..ce6b6dd --- /dev/null +++ b/src/utility/experimental/MINERvAUtility.hxx @@ -0,0 +1,37 @@ +#pragma once + +#include "utility/PhaseSpaceRestriction.hxx" + +#include //for size_t +#include + +namespace nuis { +namespace event { +class FullEvent; +} // namespace event +} // namespace nuis + +namespace nuis { +namespace utility { +namespace mnv { + +std::pair IsCC0Pi_NumProtons(event::FullEvent const &); + +bool IsCCIncLowRecoil(event::FullEvent const &); +bool IsCC0PiNp(event::FullEvent const &); +bool IsCC0PiNp_STV(event::FullEvent const &); + +bool IsCC1Pi0_2016(event::FullEvent const &); +bool IsCC1CPi_2017(event::FullEvent const &); + +TVector3 GetDeltaPT_CC0PiN_mnv(event::FullEvent const &fev); +double GetDeltaPhiT_CC0PiN_mnv(event::FullEvent const &fev); +double GetDeltaAlphaT_CC0PiN_mnv(event::FullEvent const &fev); + +double GetNeutronMomentumReco_CC0PiN_mnv(event::FullEvent const &fev); + +extern SimpleParticlePhaseSpaceRestriction CC0PiNProt_ProtonPS; + +} // namespace mnv +} // namespace utility +} // namespace nuis diff --git a/src/utility/experimental/T2KUtility.cxx b/src/utility/experimental/T2KUtility.cxx new file mode 100644 index 0000000..f0301fc --- /dev/null +++ b/src/utility/experimental/T2KUtility.cxx @@ -0,0 +1,113 @@ +#include "utility/experimental/T2KUtility.hxx" + +#include "event/FullEvent.hxx" + +#include "utility/EventTopologyUtility.hxx" +#include "utility/FullEventUtility.hxx" +#include "utility/PhaseSpaceRestriction.hxx" + +namespace nuis { +namespace utility { +namespace t2k { +std::pair IsCC0Pi_NumProtons(event::FullEvent const &ev) { + if (!nuis::utility::IsCC0Pi(ev)) { + return {false, 0}; + } + + static SimpleParticlePhaseSpaceRestriction prot_mom_cut = + SimpleParticlePhaseSpaceRestriction::Momentum(500); + + return {true, + GetNParticlesInPhaseSpace(ev, {pdgcodes::kProton}, prot_mom_cut)}; +} + +bool IsCC0PiNp(event::FullEvent const &ev) { + auto res = IsCC0Pi_NumProtons(ev); + return res.first && (res.second > 1); +} +bool IsCC0Pi1p(event::FullEvent const &ev) { + auto res = IsCC0Pi_NumProtons(ev); + return res.first && (res.second == 1); +} +bool IsCC0Pi0p(event::FullEvent const &ev) { + auto res = IsCC0Pi_NumProtons(ev); + return res.first && (res.second == 0); +} +bool IsCC0Pi_STV(event::FullEvent const &ev) { + if (!nuis::utility::IsCC0Pi(ev)) { + return false; + } + + //! Cut on kinematic properties of the true final state. + static SimpleParticlePhaseSpaceRestriction muon_PS = + SimpleParticlePhaseSpaceRestriction::MomentumCosTheta( + {250, std::numeric_limits::max()}, {-0.6, 1}); + // Muon phase space + // Pmu > 250 MeV, cos(theta_mu) > -0.6 (Sweet phase space!) + if (GetNParticlesInPhaseSpace(ev, {pdgcodes::kMu}, muon_PS) != 1) { + return false; + } + + static SimpleParticlePhaseSpaceRestriction prot_PS = + SimpleParticlePhaseSpaceRestriction::MomentumCosTheta({450, 1E3}, + {0.4, 1}); + + // Proton phase space + if (GetNParticlesInPhaseSpace(ev, {pdgcodes::kProton}, prot_PS) < 1) { + return false; + } + + return true; +} + +bool IsCC1Pip_CH_MichTag(event::FullEvent const &ev) { + if (!nuis::utility::IsCC1Pi(ev, {pdgcodes::kPiPlus})) { + return false; + } + + static SimpleParticlePhaseSpaceRestriction muon_ps = + SimpleParticlePhaseSpaceRestriction::MomentumCosTheta( + {200, std::numeric_limits::max()}, {0.2, 1}); + + if (GetNParticlesInPhaseSpace(ev, {pdgcodes::kMu}, muon_ps) != 1) { + return false; + } + return true; +} + +bool IsCC1Pip_CH_RecPi(event::FullEvent const &ev) { + if (!IsCC1Pip_CH_MichTag(ev)) { + return false; + } + + static SimpleParticlePhaseSpaceRestriction pi_ps = + SimpleParticlePhaseSpaceRestriction::MomentumCosTheta( + {200, std::numeric_limits::max()}, {0.2, 1}); + + if (GetNParticlesInPhaseSpace(ev, {pdgcodes::kPiPlus}, pi_ps) != 1) { + return false; + } + return true; +} + +bool IsCC1Pip_H20(event::FullEvent const &ev) { + if (!nuis::utility::IsCC1Pi(ev, {pdgcodes::kPiPlus})) { + return false; + } + + static SimpleParticlePhaseSpaceRestriction muon_pi_ps = + SimpleParticlePhaseSpaceRestriction::MomentumCosTheta( + {200, std::numeric_limits::max()}, {0.3, 1}); + + if (GetNParticlesInPhaseSpace(ev, {pdgcodes::kMu}, muon_pi_ps) != 1) { + return false; + } + if (GetNParticlesInPhaseSpace(ev, {pdgcodes::kPiPlus}, muon_pi_ps) != 1) { + return false; + } + return true; +} + +} // namespace t2k +} // namespace utility +} // namespace nuis diff --git a/src/utility/experimental/T2KUtility.hxx b/src/utility/experimental/T2KUtility.hxx new file mode 100644 index 0000000..8bbed3e --- /dev/null +++ b/src/utility/experimental/T2KUtility.hxx @@ -0,0 +1,44 @@ +#pragma once + +#include +#include //for size_t + +namespace nuis { +namespace event { +class FullEvent; +} // namespace event +} // namespace nuis + +namespace nuis { +namespace utility { +namespace t2k { + +std::pair IsCC0Pi_NumProtons(event::FullEvent const &ev); +bool IsCC0PiNp(event::FullEvent const &); +bool IsCC0Pi1p(event::FullEvent const &); +bool IsCC0Pi0p(event::FullEvent const &); +bool IsCC0Pi_STV(event::FullEvent const &); + +/// \brief T2K CC1pi+ CH analysis (Raquel's thesis) +/// Has different phase space cuts depending on if using Michel tag or not +/// +/// Essentially consists of two samples: one sample which has Michel e (which we +/// can't get pion direction from); this covers backwards angles quite well. +/// Measurements including this sample does not have include pion kinematics +/// cuts one sample which has PID in FGD and TPC and no Michel e. These are +/// mostly forward-going so require a pion kinematics cut +/// +/// Essentially, cuts are: +/// 1 negative muon +/// 1 positive pion +/// Any number of nucleons +/// No other particles in the final state +bool IsCC1Pip_CH_MichTag(event::FullEvent const &); +bool IsCC1Pip_CH_RecPi(event::FullEvent const &); + + + + +} // namespace t2k +} // namespace utility +} // namespace nuis