diff --git a/CMakeLists.txt b/CMakeLists.txt index 38aae55..bea638b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,82 +1,83 @@ # 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) 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) ################################## FHICLCPP #################################### include(${CMAKE_SOURCE_DIR}/cmake/fhiclcppSetup.cmake) #Need this to be at the front LIST(REVERSE EXTRA_CXX_FLAGS) LIST(APPEND EXTRA_CXX_FLAGS -I${CMAKE_SOURCE_DIR}/src) LIST(REVERSE EXTRA_CXX_FLAGS) ################################## COMPILER #################################### include(${CMAKE_SOURCE_DIR}/cmake/c++CompilerSetup.cmake) ################################################################################ add_subdirectory(src) configure_file(cmake/setup.sh.in "${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.sh" @ONLY) install(FILES "${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/setup.sh" DESTINATION ${CMAKE_INSTALL_PREFIX}) add_subdirectory(config) +add_subdirectory(data) diff --git a/cmake/NuWroSetup.cmake b/cmake/NuWroSetup.cmake index 332c32e..971cdf7 100644 --- a/cmake/NuWroSetup.cmake +++ b/cmake/NuWroSetup.cmake @@ -1,48 +1,48 @@ # 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 . ################################################################################ if(NUWRO STREQUAL "") cmessage(FATAL_ERROR "Variable NUWRO is not defined. " "This must be set to point to a prebuilt NuWro instance.") endif() -LIST(APPEND EXTRA_CXX_FLAGS -DNUWRO_ENABLED) +LIST(APPEND EXTRA_CXX_FLAGS -DNUWRO_ENABLED -Wno-sign-compare -Wno-unused-variable -Wno-reorder) LIST(APPEND EXTRA_CXX_FLAGS -I${NUWRO}/src) if(NOT EXISTS ${NUWRO}/bin/event1.so) if(EXISTS ${NUWRO}/build/${CMAKE_SYSTEM_NAME}/lib) if(NUWRO_INC STREQUAL "") cmessage(FATAL_ERROR "Variable NUWRO_INC is not defined. " "This must be set to point to an installed NuWro instance.") endif() LIST(APPEND EXTRA_LINK_DIRS ${NUWRO}/build/${CMAKE_SYSTEM_NAME}/lib) LIST(APPEND EXTRA_LIBS event) else() cmessage(FATAL_ERROR "Expected to find the NuWro event library in: ${NUWRO}/bin/event1.so, or if using NuWro with reweight support: ${NUWRO}/build/${CMAKE_SYSTEM_NAME}/lib/libevent.a. Is NuWro built?") endif() else() LIST(APPEND EXTRA_SHAREDOBJS ${NUWRO}/bin/event1.so) endif() set(NEED_PYTHIA6 TRUE) set(NEED_ROOTPYTHIA6 TRUE) diff --git a/cmake/c++CompilerSetup.cmake b/cmake/c++CompilerSetup.cmake index 205fd19..0443100 100644 --- a/cmake/c++CompilerSetup.cmake +++ b/cmake/c++CompilerSetup.cmake @@ -1,83 +1,83 @@ # 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 . ################################################################################ set(CXX_WARNINGS -Wall ) -LIST(APPEND EXTRA_CXX_FLAGS -std=c++1y -fPIC) +LIST(APPEND EXTRA_CXX_FLAGS -std=c++1y -fPIC "-D__FILENAME__='\"$(subst ${CMAKE_SOURCE_DIR}/,,$(abspath $<))\"'") cmessage(DEBUG "EXTRA_CXX_FLAGS: ${EXTRA_CXX_FLAGS}") string(REPLACE ";" " " STR_EXTRA_CXX_FLAGS "${EXTRA_CXX_FLAGS}") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${STR_EXTRA_CXX_FLAGS} ${CXX_WARNINGS}") set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0") set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -O3") SET(STR_EXTRA_LINK_DIRS) if(NOT EXTRA_LINK_DIRS STREQUAL "") string(REPLACE ";" " -L" STR_EXTRA_LINK_DIRS "-L${EXTRA_LINK_DIRS}") endif() LIST(APPEND EXTRA_LIBS dl) SET(STR_EXTRA_LIBS) if(NOT EXTRA_LIBS STREQUAL "") SET(STR_EXTRA_LIBS_NO_SCRUB_LINKOPTS) string(REPLACE ";" " -l" STR_EXTRA_LIBS_NO_SCRUB_LINKOPTS "-l${EXTRA_LIBS}") string(REPLACE "-l-" "-" STR_EXTRA_LIBS ${STR_EXTRA_LIBS_NO_SCRUB_LINKOPTS}) endif() SET(STR_EXTRA_SHAREDOBJS) if(NOT EXTRA_SHAREDOBJS STREQUAL "") string(REPLACE ";" " " STR_EXTRA_SHAREDOBJS "${EXTRA_SHAREDOBJS}") endif() LIST(APPEND EXTRA_LINK_FLAGS -Wl,--no-as-needed) SET(STR_EXTRA_LINK_FLAGS) if(NOT EXTRA_LINK_FLAGS STREQUAL "") string(REPLACE ";" " " STR_EXTRA_LINK_FLAGS "${EXTRA_LINK_FLAGS}") endif() cmessage(DEBUG "EXTRA_LINK_DIRS: ${STR_EXTRA_LINK_DIRS}") cmessage(DEBUG "EXTRA_LIBS: ${STR_EXTRA_LIBS}") cmessage(DEBUG "EXTRA_SHAREDOBJS: ${STR_EXTRA_SHAREDOBJS}") cmessage(DEBUG "EXTRA_LINK_FLAGS: ${STR_EXTRA_LINK_FLAGS}") if(NOT STR_EXTRA_LINK_DIRS STREQUAL "" AND NOT STR_EXTRA_LIBS STREQUAL "") SET(CMAKE_DEPENDLIB_FLAGS "${STR_EXTRA_LINK_DIRS} ${STR_EXTRA_LIBS}") endif() if(NOT STR_EXTRA_SHAREDOBJS STREQUAL "") SET(CMAKE_DEPENDLIB_FLAGS "${CMAKE_DEPENDLIB_FLAGS} ${STR_EXTRA_SHAREDOBJS}") endif() if(NOT EXTRA_LINK_FLAGS STREQUAL "") if(NOT CMAKE_LINK_FLAGS STREQUAL "") SET(CMAKE_LINK_FLAGS "${CMAKE_LINK_FLAGS} ${STR_EXTRA_LINK_FLAGS}") else() SET(CMAKE_LINK_FLAGS "${STR_EXTRA_LINK_FLAGS}") endif() endif() if (VERBOSE) cmessage (STATUS "C++ Compiler : ${CXX_COMPILER_NAME}") cmessage (STATUS " flags : ${CMAKE_CXX_FLAGS}") cmessage (STATUS " Release flags : ${CMAKE_CXX_FLAGS_RELEASE}") cmessage (STATUS " Debug flags : ${CMAKE_CXX_FLAGS_DEBUG}") cmessage (STATUS " Link Flags : ${CMAKE_LINK_FLAGS}") cmessage (STATUS " Lib Flags : ${CMAKE_DEPENDLIB_FLAGS}") endif() diff --git a/config/global/CMakeLists.txt b/config/global/CMakeLists.txt index 307d738..caaac25 100644 --- a/config/global/CMakeLists.txt +++ b/config/global/CMakeLists.txt @@ -1,13 +1,18 @@ configure_file(nuis.global.config.fcl.in "${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/nuis.global.config.fcl" @ONLY) install(FILES "${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/nuis.global.config.fcl" DESTINATION ${CMAKE_INSTALL_PREFIX}/fcl) +SET(INuADataComparisons_FHiCL_Includes) +foreach(fcl ${INuADataComparisons_FHiCL}) + SET(INuADataComparisons_FHiCL_Includes "${INuADataComparisons_FHiCL_Includes}#include \"${fcl}\"\n") +endforeach() + configure_file(nuis.datacomparisons.fcl.in "${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/nuis.datacomparisons.fcl" @ONLY) install(FILES "${PROJECT_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/nuis.datacomparisons.fcl" DESTINATION ${CMAKE_INSTALL_PREFIX}/fcl) diff --git a/config/global/nuis.datacomparisons.fcl.in b/config/global/nuis.datacomparisons.fcl.in index 565e6e9..471d2a8 100644 --- a/config/global/nuis.datacomparisons.fcl.in +++ b/config/global/nuis.datacomparisons.fcl.in @@ -1,6 +1,8 @@ +@INuADataComparisons_FHiCL_Includes@ + data_comparisons: { nuA: [@INuADataComparisons_List@] eA: [] piA: [] NA: [] } diff --git a/config/global/nuis.global.config.fcl.in b/config/global/nuis.global.config.fcl.in index aeb5331..5b964a4 100644 --- a/config/global/nuis.global.config.fcl.in +++ b/config/global/nuis.global.config.fcl.in @@ -1,12 +1,18 @@ plugins: { search_paths: { #Default plugin install directory installed_plugins: [ @CMAKE_INSTALL_PREFIX@/plugins ] } } data_dir: @CMAKE_INSTALL_PREFIX@/data persistency: { default_output_file: nuis.out.root } + +global: { + sample: { + verbosity_default: Reticent + } +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5fd3a7e..b98e8d5 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,12 +1,13 @@ add_subdirectory(config) add_subdirectory(event) add_subdirectory(input) add_subdirectory(plugins) add_subdirectory(utility) add_subdirectory(generator) add_subdirectory(persistency) add_subdirectory(samples) SET(INuADataComparisons_List ${INuADataComparisons_List} PARENT_SCOPE) +SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE) add_subdirectory(app) diff --git a/src/app/CMakeLists.txt b/src/app/CMakeLists.txt index 6b6fc31..3f12000 100644 --- a/src/app/CMakeLists.txt +++ b/src/app/CMakeLists.txt @@ -1,13 +1,13 @@ -SET(APPS nuissamples DumpEventInfo) +SET(APPS nuissamples nuiscomp nuisevsum) foreach(a ${APPS}) add_executable(${a} ${a}.cxx) target_link_libraries(${a} nuis_event nuis_input nuis_utility nuis_config nuis_persistency) target_link_libraries(${a} ${CMAKE_DEPENDLIB_FLAGS}) 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/DumpEventInfo.cxx b/src/app/nuiscomp.cxx similarity index 77% rename from src/app/DumpEventInfo.cxx rename to src/app/nuiscomp.cxx index a602fec..cfb36c4 100644 --- a/src/app/DumpEventInfo.cxx +++ b/src/app/nuiscomp.cxx @@ -1,52 +1,51 @@ #include "config/GlobalConfiguration.hxx" #include "input/IInputHandler.hxx" #include "event/MinimalEvent.hxx" #include "samples/ISample.hxx" #include "plugins/Instantiate.hxx" #include "exception/exception.hxx" #include "fhiclcpp/make_ParameterSet.h" #include NEW_NUIS_EXCEPT(invalid_cli_arguments); int main(int argc, char const *argv[]) { nuis::config::EnsureConfigurationRead("nuis.global.config.fcl"); + nuis::config::EnsureConfigurationRead("nuis.datacomparisons.fcl"); if (argc != 2) { throw invalid_cli_arguments() << "[ERROR]: Expected to be passed a single FHiCL file name or " "absolute or relative path. N.B. Files in the local directory must " "be fully qualified like \"$ " << argv[0] << " ./myconf.fcl\"."; } - fhicl::ParameterSet ps = fhicl::make_ParameterSet(argv[1]); + nuis::config::EnsureConfigurationRead(argv[1]); - size_t NMax = std::numeric_limits::max(); - - if (ps.has_key("nmax")) { - NMax = ps.get("nmax"); - } + size_t NMax = nuis::config::GetDocument().get( + "nmax", std::numeric_limits::max()); for (fhicl::ParameterSet const &samp_config : - ps.get>("samples")) { + nuis::config::GetDocument().get>( + "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(); } } diff --git a/src/app/nuisevsum.cxx b/src/app/nuisevsum.cxx new file mode 100644 index 0000000..c45f677 --- /dev/null +++ b/src/app/nuisevsum.cxx @@ -0,0 +1,79 @@ +#include "config/GlobalConfiguration.hxx" + +#include "input/IInputHandler.hxx" + +#include "event/MinimalEvent.hxx" + +#include "samples/ISample.hxx" + +#include "plugins/Instantiate.hxx" + +#include "exception/exception.hxx" + +#include "fhiclcpp/make_ParameterSet.h" +#include "string_parsers/from_string.hxx" + +#include + +NEW_NUIS_EXCEPT(invalid_cli_arguments); + +size_t NMax = std::numeric_limits::max(); +std::string input_file; +std::string input_type; + +void SayUsage(char const *argv[]) { + std::cout << "[USAGE]: " << argv[0] + << "\n" + "\t-i : Input file passed to named " + "IInputHandler instance \n" + "\t-H : Name of IInputHandler subclass " + "capable of reading NUISANCE events from the argument of -i.\n" + "\t-n : Maximum number of events to " + "read. Will read entire input file by default.\n" + << std::endl; +} + +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]) == "-i") { + input_file = argv[++opt]; + } else if (std::string(argv[opt]) == "-H") { + input_type = argv[++opt]; + } else if (std::string(argv[opt]) == "-n") { + NMax = fhicl::string_parsers::str2T(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"); + + handleOpts(argc, argv); + + if (!input_type.length() || !input_file.length()) { + SayUsage(argv); + throw invalid_cli_arguments() + << "[ERROR]: Require both -i and -H cli options to be passed."; + } + + fhicl::ParameterSet sample_config; + + sample_config.put("input_type", input_type); + sample_config.put("file", input_file); + + nuis::plugins::plugin_traits::unique_ptr_t VerboseEventSummary = + nuis::plugins::Instantiate("VerboseEventSummary"); + + VerboseEventSummary->Initialize(sample_config); + VerboseEventSummary->ProcessSample(NMax); +} diff --git a/src/app/nuissamples.cxx b/src/app/nuissamples.cxx index 68a17bc..d61d1e3 100644 --- a/src/app/nuissamples.cxx +++ b/src/app/nuissamples.cxx @@ -1,36 +1,36 @@ #include "config/GlobalConfiguration.hxx" #include "input/IInputHandler.hxx" #include "plugins/Instantiate.hxx" #include "samples/IDataComparison.hxx" #include #include int main() { nuis::config::EnsureConfigurationRead("nuis.global.config.fcl"); nuis::config::EnsureConfigurationRead("nuis.datacomparisons.fcl"); for (std::string const &comparison_set_key : nuis::config::GetDocument() .get("data_comparisons") .get_names()) { for (std::string const &sample_name : nuis::config::GetDocument().get>( std::string("data_comparisons.") + comparison_set_key)) { nuis::plugins::plugin_traits::unique_ptr_t sample = nuis::plugins::Instantiate(sample_name); - // - // std::cout << sample->Name() << std::endl; - // std::cout << "\tJournal: " << sample->GetJournalReference() << std::endl; - // std::cout << "\tTarget: " << sample->GetTargetMaterial() << std::endl; - // std::cout << "\tFlux: " << sample->GetFluxDescription() << std::endl; - // std::cout << "\tSignal: " << sample->GetSignalDescription() << std::endl; - // std::cout << "\tDocs: " << sample->GetDocumentation() << std::endl; + + std::cout << sample->Name() << std::endl; + std::cout << "\tJournal: " << sample->GetJournalReference() << std::endl; + std::cout << "\tTarget: " << sample->GetTargetMaterial() << std::endl; + std::cout << "\tFlux: " << sample->GetFluxDescription() << std::endl; + std::cout << "\tSignal: " << sample->GetSignalDescription() << std::endl; + std::cout << "\tDocs: " << sample->GetDocumentation() << std::endl; } } } diff --git a/src/event/FullEvent.hxx b/src/event/FullEvent.hxx index 12cbfe9..74e4662 100644 --- a/src/event/FullEvent.hxx +++ b/src/event/FullEvent.hxx @@ -1,47 +1,70 @@ // 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_FULLEVENT_HXX_SEEN #define EVENT_FULLEVENT_HXX_SEEN #include "event/MinimalEvent.hxx" #include "event/Particle.hxx" +#include "string_parsers/to_string.hxx" + +#include + namespace nuis { namespace event { ///\brief The full, internal event format. class FullEvent : public MinimalEvent { public: struct StatusParticles { Particle::Status_t status; std::vector particles; }; FullEvent(); - FullEvent(FullEvent const&) = delete; - FullEvent(FullEvent&&); + FullEvent(FullEvent const &) = delete; + FullEvent(FullEvent &&); std::vector ParticleStack; void ClearParticleStack(); + + std::string to_string() const { + std::stringstream ss(""); + ss << "Event: Interaction mode = " << mode + << ", probe: { PDG: " << probe_pdg << ", Energy: " << probe_E + << " MeV }." << std::endl; + for (auto &status_stack : ParticleStack) { + ss << "\t[" << status_stack.status << "]" << std::endl; + + for (Particle const &part : status_stack.particles) { + ss << "\t\t{ PDG: " << part.pdg << ", P3: [ " << part.P4[0] << ", " + << part.P4[1] << ", " << part.P4[2] << "], E: " << part.P4[3] + << ", M: " << part.P4.M() << " }" << std::endl; + } + } + ss << std::endl; + return ss.str(); + } }; -} // namespace core +} // namespace event } // namespace nuis + #endif diff --git a/src/event/Particle.hxx b/src/event/Particle.hxx index accfe51..0bcad11 100644 --- a/src/event/Particle.hxx +++ b/src/event/Particle.hxx @@ -1,68 +1,77 @@ // 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_PARTICLE_HXX_SEEN #define EVENT_PARTICLE_HXX_SEEN #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!() { return (pdg == std::numeric_limits::max()); } + bool operator!() const { return (pdg == std::numeric_limits::max()); } + + double E() const { return P4.E(); } + double P() const { return P4.Vect().Mag(); } + TVector3 P3() const { return P4.Vect(); } + double M() const { return P4.M(); } + double Theta() const { return P4.Vect().Theta(); } + double CosTheta() const { return P4.Vect().Theta(); } }; } // namespace event } // namespace nuis #define X(A, B) \ - case nuis::event::Particle::Status_t::A: { \ + 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 #endif diff --git a/src/generator/input/NEUTInputHandler.cxx b/src/generator/input/NEUTInputHandler.cxx index f1b4920..a4133b6 100644 --- a/src/generator/input/NEUTInputHandler.cxx +++ b/src/generator/input/NEUTInputHandler.cxx @@ -1,106 +1,103 @@ #include "generator/input/NEUTInputHandler.hxx" -#include "utility/ROOTUtility.hxx" -#include "utility/PDGCodeUtility.hxx" #include "utility/InteractionChannelUtility.hxx" +#include "utility/PDGCodeUtility.hxx" +#include "utility/ROOTUtility.hxx" #include "generator/utility/NEUTUtility.hxx" #include "fhiclcpp/ParameterSet.h" using namespace nuis::event; using namespace nuis::utility; +using namespace nuis::neuttools; NEUTInputHandler::NEUTInputHandler() : fInputTree(nullptr) {} NEUTInputHandler::NEUTInputHandler(NEUTInputHandler &&other) : fInputTree(std::move(other.fInputTree)), fReaderEvent(std::move(other.fReaderEvent)) {} void NEUTInputHandler::Initialize(fhicl::ParameterSet const &ps) { fInputTree = CheckGetTTree(ps.get("file"), "neuttree"); fReaderEvent.fNeutVect = nullptr; fInputTree->tree->SetBranchAddress("vectorbranch", &fReaderEvent.fNeutVect); + + fKeepIntermediates = ps.get("keep_intermediates", false); } 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(); } fInputTree->tree->GetEntry(idx); fReaderEvent.mode = IntToChannel(fReaderEvent.fNeutVect->Mode); size_t NPart = fReaderEvent.fNeutVect->Npart(); for (size_t part_it = 0; part_it < NPart; part_it++) { NeutPart const &part = (*fReaderEvent.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; } } return fReaderEvent; } FullEvent const &NEUTInputHandler::GetFullEvent(ev_index_t idx) const { (void)GetMinimalEvent(idx); fReaderEvent.ClearParticleStack(); if (fReaderEvent.fNeutVect->Ibound) { fReaderEvent.target_pdg = MakeNuclearPDG(fReaderEvent.fNeutVect->TargetA, fReaderEvent.fNeutVect->TargetZ); } else { fReaderEvent.target_pdg = MakeNuclearPDG(1, 1); } size_t NPart = fReaderEvent.fNeutVect->Npart(); - bool FoundIntermediateStateParticle = false; + size_t NPrimary = fReaderEvent.fNeutVect->Nprimary(); for (size_t part_it = 0; part_it < NPart; part_it++) { NeutPart const &part = (*fReaderEvent.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); - size_t state_int = static_cast(state); - if ((!FoundIntermediateStateParticle) && - (state == Particle::Status_t::kIntermediate)) { - FoundIntermediateStateParticle = true; + 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 (!IsCoh(fReaderEvent.mode) && (part_it > 1) && (state_int == 0) && - (!FoundIntermediateStateParticle)) { + if ((part_it < NPrimary) && + (state != Particle::Status_t::kPrimaryInitialState)) { fReaderEvent .ParticleStack[static_cast( Particle::Status_t::kPrimaryFinalState)] .particles.push_back(nuis_part); } - // Intermediate particles should be pushed onto the primary final state - // stack for NEUT - if (state == Particle::Status_t::kIntermediate) { - state_int = static_cast(Particle::Status_t::kPrimaryFinalState); - } - fReaderEvent.ParticleStack[state_int].particles.push_back(nuis_part); } return fReaderEvent; } size_t NEUTInputHandler::GetNEvents() const { return fInputTree->tree->GetEntries(); } DECLARE_PLUGIN(IInputHandler, NEUTInputHandler); diff --git a/src/generator/input/NEUTInputHandler.hxx b/src/generator/input/NEUTInputHandler.hxx index 8216d9b..9f4d8f6 100644 --- a/src/generator/input/NEUTInputHandler.hxx +++ b/src/generator/input/NEUTInputHandler.hxx @@ -1,54 +1,56 @@ // 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 GENERATOR_INPUT_NEUTINPUTHANDLER_HXX_SEEN #define GENERATOR_INPUT_NEUTINPUTHANDLER_HXX_SEEN #include "event/FullEvent.hxx" #include "input/IInputHandler.hxx" #include namespace fhicl { class ParameterSet; } namespace nuis { namespace utility { class TreeFile; } } // namespace nuis class NEUTInputHandler : public IInputHandler { mutable std::unique_ptr fInputTree; mutable nuis::event::FullEvent fReaderEvent; + bool fKeepIntermediates; + public: NEUTInputHandler(); NEUTInputHandler(NEUTInputHandler const &) = delete; NEUTInputHandler(NEUTInputHandler &&); 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; size_t GetNEvents() const; }; #endif diff --git a/src/generator/input/NuWroInputHandler.cxx b/src/generator/input/NuWroInputHandler.cxx index 8d40dd8..2fd3d52 100644 --- a/src/generator/input/NuWroInputHandler.cxx +++ b/src/generator/input/NuWroInputHandler.cxx @@ -1,42 +1,115 @@ #include "generator/input/NuWroInputHandler.hxx" +#include "generator/utility/NuWroUtility.hxx" + #include "utility/ROOTUtility.hxx" #include "fhiclcpp/ParameterSet.h" +#include "particle.h" +typedef ::particle NuWroParticle; + using namespace nuis::event; using namespace nuis::utility; +using namespace nuis::nuwrotools; NuWroInputHandler::NuWroInputHandler() : fInputTree(nullptr) {} NuWroInputHandler::NuWroInputHandler(NuWroInputHandler &&other) : fInputTree(std::move(other.fInputTree)), fReaderEvent(std::move(other.fReaderEvent)) {} void NuWroInputHandler::Initialize(fhicl::ParameterSet const &ps) { fInputTree = CheckGetTTree(ps.get("file"), "treeout"); fReaderEvent.fNuWroEvent = nullptr; fInputTree->tree->SetBranchAddress("e", &fReaderEvent.fNuWroEvent); + + fKeepIntermediates = ps.get("keep_intermediates", false); } MinimalEvent const &NuWroInputHandler::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); + + fReaderEvent.mode = NuWroEventChannel(*fReaderEvent.fNuWroEvent); + fReaderEvent.probe_E = fReaderEvent.fNuWroEvent->in[0].E(); + fReaderEvent.probe_pdg = fReaderEvent.fNuWroEvent->in[0].pdg; + fReaderEvent.XSecWeight = + fReaderEvent.fNuWroEvent->weight / double(GetNEvents()); + + if (fWeightCache.size() <= idx) { + fWeightCache.push_back(fReaderEvent.XSecWeight); + } + return fReaderEvent; } FullEvent const &NuWroInputHandler::GetFullEvent(ev_index_t idx) const { (void)GetMinimalEvent(idx); - // Fill particle stack + + fReaderEvent.ClearParticleStack(); + + for (size_t p_it = 0; p_it < fReaderEvent.fNuWroEvent->in.size(); ++p_it) { + NuWroParticle &part = fReaderEvent.fNuWroEvent->in[p_it]; + + Particle nuis_part; + + nuis_part.pdg = part.pdg; + nuis_part.P4 = TLorentzVector(part[1], part[2], part[3], part[0]); + + fReaderEvent + .ParticleStack[static_cast( + Particle::Status_t::kPrimaryInitialState)] + .particles.push_back(nuis_part); + } + + for (size_t p_it = 0; + p_it < fKeepIntermediates && fReaderEvent.fNuWroEvent->out.size(); + ++p_it) { + NuWroParticle &part = fReaderEvent.fNuWroEvent->out[p_it]; + + Particle nuis_part; + + nuis_part.pdg = part.pdg; + nuis_part.P4 = TLorentzVector(part[1], part[2], part[3], part[0]); + + fReaderEvent + .ParticleStack[static_cast( + Particle::Status_t::kPrimaryFinalState)] + .particles.push_back(nuis_part); + } + for (size_t p_it = 0; (p_it < fReaderEvent.fNuWroEvent->post.size()); + ++p_it) { + NuWroParticle &part = fReaderEvent.fNuWroEvent->post[p_it]; + + Particle nuis_part; + + nuis_part.pdg = part.pdg; + nuis_part.P4 = TLorentzVector(part[1], part[2], part[3], part[0]); + + fReaderEvent + .ParticleStack[static_cast(Particle::Status_t::kNuclearLeaving)] + .particles.push_back(nuis_part); + } + return fReaderEvent; } +double NuWroInputHandler::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 NuWroInputHandler::GetNEvents() const { return fInputTree->tree->GetEntries(); } DECLARE_PLUGIN(IInputHandler, NuWroInputHandler); diff --git a/src/generator/input/NuWroInputHandler.hxx b/src/generator/input/NuWroInputHandler.hxx index a2567ac..bb050e3 100644 --- a/src/generator/input/NuWroInputHandler.hxx +++ b/src/generator/input/NuWroInputHandler.hxx @@ -1,54 +1,63 @@ // 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 GENERATOR_INPUT_NUWROINPUTHANDLER_HXX_SEEN #define GENERATOR_INPUT_NUWROINPUTHANDLER_HXX_SEEN #include "event/FullEvent.hxx" #include "input/IInputHandler.hxx" +#include "exception/exception.hxx" + #include namespace fhicl { class ParameterSet; } namespace nuis { namespace utility { class TreeFile; } } // namespace nuis class NuWroInputHandler : public IInputHandler { mutable std::unique_ptr fInputTree; mutable nuis::event::FullEvent fReaderEvent; + mutable std::vector fWeightCache; + + bool fKeepIntermediates; public: + + NEW_NUIS_EXCEPT(weight_cache_miss); + NuWroInputHandler(); NuWroInputHandler(NuWroInputHandler const &) = delete; NuWroInputHandler(NuWroInputHandler &&); 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; }; #endif diff --git a/src/generator/utility/CMakeLists.txt b/src/generator/utility/CMakeLists.txt index 6075970..7b60771 100644 --- a/src/generator/utility/CMakeLists.txt +++ b/src/generator/utility/CMakeLists.txt @@ -1,15 +1,20 @@ +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(GENERATOR_UTILS_IMPL) add_library(nuis_generator_utility SHARED ${GENERATOR_UTILS_IMPL}) target_link_libraries(nuis_generator_utility nuis_event) install(TARGETS nuis_generator_utility DESTINATION lib) endif() if(GENERATOR_UTILS_HDR) install(FILES ${GENERATOR_UTILS_HDR} DESTINATION include/generator/utility) endif() diff --git a/src/generator/utility/NEUTUtility.cxx b/src/generator/utility/NEUTUtility.cxx index d3260ca..56cd4a9 100644 --- a/src/generator/utility/NEUTUtility.cxx +++ b/src/generator/utility/NEUTUtility.cxx @@ -1,86 +1,90 @@ #include "generator/utility/NEUTUtility.hxx" #include "exception/exception.hxx" #include "utility/InteractionChannelUtility.hxx" #include "utility/PDGCodeUtility.hxx" #include "neutpart.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 << ", Part: { Status: " << part.fStatus << ", IsAlive: " << part.fIsAlive << ", PDG: " << 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; } // 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; } + +} // namespace neuttools +} // namespace nuis diff --git a/src/generator/utility/NEUTUtility.hxx b/src/generator/utility/NEUTUtility.hxx index ea2b3ab..c191e52 100644 --- a/src/generator/utility/NEUTUtility.hxx +++ b/src/generator/utility/NEUTUtility.hxx @@ -1,11 +1,15 @@ #ifndef GENERATOR_UTILITY_NEUTUTILITY_HXX_SEEN #define GENERATOR_UTILITY_NEUTUTILITY_HXX_SEEN #include "event/Particle.hxx" #include "event/types.hxx" class NeutPart; -nuis::event::Particle::Status_t -GetNeutParticleStatus(NeutPart const &, nuis::event::Channel_t); +namespace nuis { +namespace neuttools { +nuis::event::Particle::Status_t GetNeutParticleStatus(NeutPart const &, + nuis::event::Channel_t); +} +} // namespace nuis #endif diff --git a/src/generator/utility/NuWroUtility.cxx b/src/generator/utility/NuWroUtility.cxx new file mode 100644 index 0000000..d42456e --- /dev/null +++ b/src/generator/utility/NuWroUtility.cxx @@ -0,0 +1,400 @@ +#include "generator/utility/NuWroUtility.hxx" + +#include "exception/exception.hxx" + +using namespace nuis::event; + +namespace nuis { +namespace nuwrotools { + +NEW_NUIS_EXCEPT(invalid_channel_to_NuWro); +NEW_NUIS_EXCEPT(invalid_NuWro_dyn_found); + +std::pair GetFlagsDynEquivalent(Channel_t chan) { + + std::pair NuMode; + + NuMode.first.qel = false; + NuMode.first.res = false; + NuMode.first.dis = false; + NuMode.first.coh = false; + NuMode.first.mec = false; + NuMode.first.hip = false; + NuMode.first.nc = false; + NuMode.first.cc = false; + NuMode.first.anty = false; + NuMode.first.res_delta = false; + + switch (chan) { + case Channel_t::kCCQE: { + NuMode.first.qel = true; + NuMode.first.cc = true; + break; + } + case Channel_t::kCC2p2h: { + NuMode.first.mec = true; + NuMode.first.cc = true; + break; + } + case Channel_t::kCCSPP_PPip: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.cc = true; + break; + } + case Channel_t::kCCSPP_PPi0: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.cc = true; + break; + } + case Channel_t::kCCSPP_NPip: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.cc = true; + break; + } + case Channel_t::kCCCohPi: { + NuMode.first.coh = true; + NuMode.first.cc = true; + break; + } + case Channel_t::kCCResGamma: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.cc = true; + break; + } + case Channel_t::kCCTransitionMPi: { + NuMode.first.res = true; + NuMode.first.cc = true; + break; + } + case Channel_t::kCCResEta0: { + NuMode.first.res = true; + NuMode.first.cc = true; + break; + } + case Channel_t::kCCResK: { + NuMode.first.res = true; + NuMode.first.cc = true; + break; + } + case Channel_t::kCCDIS: { + NuMode.first.dis = true; + NuMode.first.cc = true; + break; + } + case Channel_t::kNCSPP_NPi0: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.nc = true; + break; + } + case Channel_t::kNCSPP_PPi0: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.nc = true; + break; + } + case Channel_t::kNCSPP_PPim: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.nc = true; + break; + } + case Channel_t::kNCSPP_NPip: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.nc = true; + break; + } + case Channel_t::kNCCohPi: { + NuMode.first.coh = true; + NuMode.first.nc = true; + break; + } + case Channel_t::kNCResNGamma: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.nc = true; + break; + } + case Channel_t::kNCResPGamma: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.nc = true; + break; + } + case Channel_t::kNCTransitionMPi: { + NuMode.first.res = true; + NuMode.first.nc = true; + break; + } + case Channel_t::kNCResNEta0: { + NuMode.first.res = true; + NuMode.first.nc = true; + break; + } + case Channel_t::kNCResPEta0: { + NuMode.first.res = true; + NuMode.first.nc = true; + break; + } + case Channel_t::kNCResK0: { + NuMode.first.res = true; + NuMode.first.nc = true; + break; + } + case Channel_t::kNCResKp: { + NuMode.first.res = true; + NuMode.first.nc = true; + break; + } + case Channel_t::kNCDIS: { + NuMode.first.dis = true; + NuMode.first.nc = true; + break; + } + case Channel_t::kNCELP: { + NuMode.first.qel = true; + NuMode.first.nc = true; + break; + } + case Channel_t::kNCELN: { + NuMode.first.qel = true; + NuMode.first.nc = true; + break; + } + case Channel_t::kCCQE_nub: { + NuMode.first.qel = true; + NuMode.first.cc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kCC2p2h_nub: { + NuMode.first.mec = true; + NuMode.first.cc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kCCSPP_NPim_nub: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.cc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kCCSPP_NPi0_nub: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.cc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kCCSPP_PPim_nub: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.cc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kCCCohPi_nub: { + NuMode.first.coh = true; + NuMode.first.cc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kCCResGamma_nub: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.cc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kCCTransitionMPi_nub: { + NuMode.first.res = true; + NuMode.first.cc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kCCResEta0_nub: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.cc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kCCResK_nub: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.cc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kCCDIS_nub: { + NuMode.first.dis = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kNCSPP_NPi0_nub: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.nc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kNCSPP_PPi0_nub: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.nc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kNCSPP_PPim_nub: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.nc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kNCSPP_NPip_nub: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.nc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kNCCohPi_nub: { + NuMode.first.coh = true; + NuMode.first.nc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kNCResNGamma_nub: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.nc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kNCResPGamma_nub: { + NuMode.first.res = true; + NuMode.first.res_delta = true; + NuMode.first.nc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kNCTransitionMPi_nub: { + NuMode.first.res = true; + NuMode.first.nc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kNCResNEta0_nub: { + NuMode.first.res = true; + NuMode.first.nc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kNCResPEta0_nub: { + NuMode.first.res = true; + NuMode.first.nc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kNCResK0_nub: { + NuMode.first.res = true; + NuMode.first.nc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kNCResKp_nub: { + NuMode.first.res = true; + NuMode.first.nc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kNCDIS_nub: { + NuMode.first.dis = true; + NuMode.first.nc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kNCELP_nub: { + NuMode.first.qel = true; + NuMode.first.nc = true; + NuMode.first.anty = true; + break; + } + case Channel_t::kNCELN_nub: { + NuMode.first.qel = true; + NuMode.first.nc = true; + NuMode.first.anty = true; + break; + } + + case Channel_t::kUndefined: + default: { + throw invalid_channel_to_NuWro() + << "[ERROR]: Encountered unhandled Channel_t: " << chan + << " when build Nuwro-equivalent dyn and flags."; + } + } + + NuMode.second = 0; + if (NuMode.first.qel) { + NuMode.second = NuMode.first.cc ? 0 : 1; + } + if (NuMode.first.res) { + NuMode.second = NuMode.first.cc ? 2 : 3; + } + if (NuMode.first.dis) { + NuMode.second = NuMode.first.cc ? 4 : 5; + } + if (NuMode.first.coh) { + NuMode.second = NuMode.first.cc ? 6 : 7; + } + if (NuMode.first.mec) { + NuMode.second = NuMode.first.cc ? 8 : 9; + } + + return NuMode; +} + +Channel_t NuWroEventChannel(NuWroEvent const &ev) { + switch (ev.dyn) { + case 0: { + return ev.flag.anty ? Channel_t::kCCQE_nub : Channel_t::kCCQE; + } + case 1: { + return ev.flag.anty ? Channel_t::kNCELN_nub : Channel_t::kNCELN; + } + case 2: { + return ev.flag.anty ? Channel_t::kCCSPP_NPi0_nub : Channel_t::kCCSPP_NPip; + } + case 3: { + return ev.flag.anty ? Channel_t::kNCSPP_NPi0_nub : Channel_t::kNCSPP_NPi0; + } + case 4: { + return ev.flag.anty ? Channel_t::kCCDIS_nub : Channel_t::kCCDIS; + } + case 5: { + return ev.flag.anty ? Channel_t::kNCDIS_nub : Channel_t::kNCDIS; + } + case 6: { + return ev.flag.anty ? Channel_t::kCCCohPi_nub : Channel_t::kCCCohPi; + } + case 7: { + return ev.flag.anty ? Channel_t::kNCCohPi_nub : Channel_t::kNCCohPi; + } + case 8: { + return ev.flag.anty ? Channel_t::kCC2p2h_nub : Channel_t::kCC2p2h; + } + default: + throw invalid_NuWro_dyn_found() << "[ERROR]: Found NuWro dyn: " << ev.dyn; + } +} + +} // namespace nuwrotools +} // namespace nuis diff --git a/src/generator/utility/NuWroUtility.hxx b/src/generator/utility/NuWroUtility.hxx new file mode 100644 index 0000000..141bcd9 --- /dev/null +++ b/src/generator/utility/NuWroUtility.hxx @@ -0,0 +1,19 @@ +#ifndef GENERATOR_UTILITY_NUWROUTILITY_HXX_SEEN +#define GENERATOR_UTILITY_NUWROUTILITY_HXX_SEEN + +#include "event/types.hxx" +#include "event/MinimalEvent.hxx" + +#include "event1.h" +typedef ::flags NuWroFlags; + +#include + +namespace nuis { +namespace nuwrotools { +std::pair GetFlagsDynEquivalent(nuis::event::Channel_t); +nuis::event::Channel_t NuWroEventChannel(NuWroEvent const &); +} +} // namespace nuis + +#endif diff --git a/src/input/IInputHandler.hxx b/src/input/IInputHandler.hxx index 8be40d3..f0a2890 100644 --- a/src/input/IInputHandler.hxx +++ b/src/input/IInputHandler.hxx @@ -1,111 +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 . *******************************************************************************/ #ifndef CORE_IINPUTHANDLER_HXX_SEEN #define CORE_IINPUTHANDLER_HXX_SEEN #include "plugins/traits.hxx" #include "exception/exception.hxx" #include namespace fhicl { class ParameterSet; } namespace nuis { namespace event { class MinimalEvent; class FullEvent; -} // namespace core +} // namespace event } // namespace nuis class IInputHandler { public: struct FullEvent_const_iterator : public std::iterator< std::input_iterator_tag, nuis::event::FullEvent const, size_t, nuis::event::FullEvent const *, nuis::event::FullEvent const &> { FullEvent_const_iterator(size_t _idx, IInputHandler const *_ih) { idx = _idx; ih = _ih; } FullEvent_const_iterator(FullEvent_const_iterator const &other) { idx = other.idx; ih = other.ih; } FullEvent_const_iterator operator=(FullEvent_const_iterator const &other) { idx = other.idx; ih = other.ih; return (*this); } bool operator==(FullEvent_const_iterator const &other) { return (idx == other.idx); } bool operator!=(FullEvent_const_iterator const &other) { return !(*this == other); } nuis::event::FullEvent const &operator*() { return ih->GetFullEvent(idx); } - nuis::event::FullEvent const *operator->() { return &ih->GetFullEvent(idx); } + nuis::event::FullEvent const *operator->() { + return &ih->GetFullEvent(idx); + } FullEvent_const_iterator operator++() { idx++; return *this; } FullEvent_const_iterator operator++(int) { FullEvent_const_iterator tmp(*this); idx++; return tmp; } private: size_t idx; IInputHandler const *ih; }; NEW_NUIS_EXCEPT(invalid_input_file); NEW_NUIS_EXCEPT(invalid_entry); typedef size_t ev_index_t; virtual void Initialize(fhicl::ParameterSet const &) = 0; virtual nuis::event::MinimalEvent const & GetMinimalEvent(ev_index_t idx) const = 0; virtual nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const = 0; virtual void RecalculateEventWeights(){}; - virtual double GetEventWeight(ev_index_t idx) const {}; + virtual double GetEventWeight(ev_index_t idx) const { return 1; }; virtual size_t GetNEvents() const = 0; FullEvent_const_iterator begin() const { return FullEvent_const_iterator(0, this); } FullEvent_const_iterator end() const { return FullEvent_const_iterator(GetNEvents(), this); } virtual ~IInputHandler() {} }; DECLARE_PLUGIN_INTERFACE(IInputHandler); #endif diff --git a/src/persistency/ROOTOutput.hxx b/src/persistency/ROOTOutput.hxx index a10c7f5..3b45a71 100644 --- a/src/persistency/ROOTOutput.hxx +++ b/src/persistency/ROOTOutput.hxx @@ -1,36 +1,66 @@ #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::shared_ptr GetOutputFile(std::string const &name = "default"); template -void WriteToOutputFile(T *object, std::string const &object_name, - std::string const &dir_name = "", - std::string const &file_name = "default") { +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::shared_ptr f = GetOutputFile(file_name); - f->cd(); - if (dir_name.length()) { - if (!f->cd(dir_name.c_str())) { - f->mkdir(dir_name.c_str())->cd(); + 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(); + } } } // namespace persistency } // namespace nuis #endif diff --git a/src/samples/CMakeLists.txt b/src/samples/CMakeLists.txt index 89b8852..ed8bfaf 100644 --- a/src/samples/CMakeLists.txt +++ b/src/samples/CMakeLists.txt @@ -1,17 +1,20 @@ set(samples_header_files -ISample.hxx) + ISample.hxx + IDataComparison.hxx + SimpleDataComparison.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() SET(INuADataComparisons_List ${INuADataComparisons_List} PARENT_SCOPE) +SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE) diff --git a/src/samples/IDataComparison.hxx b/src/samples/IDataComparison.hxx index ef9e78d..cfda1fc 100644 --- a/src/samples/IDataComparison.hxx +++ b/src/samples/IDataComparison.hxx @@ -1,72 +1,73 @@ // 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_IDATACOMPARISON_HXX_SEEN #define SAMPLES_IDATACOMPARISON_HXX_SEEN #include "samples/ISample.hxx" #include "fhiclcpp/ParameterSet.h" #include #include class IDataComparison : public ISample { public: virtual double GetGOF() = 0; + virtual double GetNDOGuess() = 0; virtual std::string GetJournalReference() { std::stringstream ss(""); ss << "Unknown Journal Ref. for IDataComparison: " << std::quoted(Name()); return ss.str(); } virtual std::string GetTargetMaterial() { std::stringstream ss(""); ss << "Unknown Target material for IDataComparison: " << std::quoted(Name()); return ss.str(); } virtual std::string GetFluxDescription() { std::stringstream ss(""); ss << "Unknown Flux description for IDataComparison: " << std::quoted(Name()); return ss.str(); } virtual std::string GetSignalDescription() { std::stringstream ss(""); ss << "Unknown Signal description for IDataComparison: " << std::quoted( Name()); return ss.str(); } virtual std::string GetDocumentation() { std::stringstream ss(""); ss << "No documentation provided for IDataComparison: " << std::quoted(Name()); return ss.str(); } virtual fhicl::ParameterSet GetExampleConfiguration() { return fhicl::ParameterSet(); } }; DECLARE_PLUGIN_INTERFACE(IDataComparison); #endif diff --git a/src/samples/ISample.hxx b/src/samples/ISample.hxx index a45b135..0b195b6 100644 --- a/src/samples/ISample.hxx +++ b/src/samples/ISample.hxx @@ -1,68 +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 . *******************************************************************************/ #ifndef SAMPLES_ISAMPLE_HXX_SEEN #define SAMPLES_ISAMPLE_HXX_SEEN #include "plugins/traits.hxx" #include "exception/exception.hxx" -#include -#include +#include "config/GlobalConfiguration.hxx" -namespace fhicl { -class ParameterSet; -} +#include "fhiclcpp/ParameterSet.h" + +#include +#include +#include namespace nuis { namespace event { class FullEvent; class MinimalEvent; -} // namespace core +} // namespace event } // namespace nuis +#define ISAMPLE_DEBUG(v) \ + if (verb > 2) { \ + std::cout << "[DEBUG]: " << v << std::endl; \ + } +#define ISAMPLE_INFO(v) \ + if (verb > 1) { \ + std::cout << "[INFO]: " << v << std::endl; \ + } +#define ISAMPLE_WARN(v) \ + if (verb > 0) { \ + std::cout << "[WARN] " << __FILENAME__ << ":" << __LINE__ << " : " << v \ + << std::endl; \ + } + class ISample { +protected: + NEW_NUIS_EXCEPT(unknown_ISample_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_ISample_verbosity() + << "[ERROR]: Failed to parse ISample verbosity setting from: " + << std::quoted(v); + } + } + public: NEW_NUIS_EXCEPT(uninitialized_ISample); NEW_NUIS_EXCEPT(unimplemented_ISample_optional_method); + ISample() { + SetSampleVerbosity(nuis::config::GetDocument().get( + "global.sample.verbosity_default", "Silent")); + } + virtual void Initialize(fhicl::ParameterSet const &) = 0; - //Interface for processing a single, external event + // Interface for processing a single, external event // // ISamples are not required to implement processing events from 'outside'. virtual void ProcessEvent(nuis::event::FullEvent const &) { throw unimplemented_ISample_optional_method() << "[ERROR]: ISample " << Name() << " does not implement ProcessEvent."; } virtual void ProcessSample(size_t nmax = std::numeric_limits::max()) = 0; virtual void Write() = 0; virtual std::string Name() = 0; virtual ~ISample() {} }; DECLARE_PLUGIN_INTERFACE(ISample); #endif diff --git a/src/samples/MCTools/CMakeLists.txt b/src/samples/MCTools/CMakeLists.txt index 86f316d..a7a2ee2 100644 --- a/src/samples/MCTools/CMakeLists.txt +++ b/src/samples/MCTools/CMakeLists.txt @@ -1,10 +1,13 @@ +LIST(APPEND SAMPLES_MCTOOLS_LIBS nuis_event nuis_input nuis_config) + if(USE_NuWro) LIST(APPEND MC_TOOL_IMPL NuisToNuWro.cxx) + LIST(APPEND SAMPLES_MCTOOLS_LIBS nuis_generator_utility) endif(USE_NuWro) LIST(APPEND MC_TOOL_IMPL VerboseEventSummary.cxx) add_library(MCTools SHARED ${MC_TOOL_IMPL}) -target_link_libraries(MCTools nuis_event nuis_input nuis_config) +target_link_libraries(MCTools ${SAMPLES_MCTOOLS_LIBS}) install(TARGETS MCTools DESTINATION plugins) diff --git a/src/samples/MCTools/NuisToNuWro.cxx b/src/samples/MCTools/NuisToNuWro.cxx index ebba240..b408c53 100644 --- a/src/samples/MCTools/NuisToNuWro.cxx +++ b/src/samples/MCTools/NuisToNuWro.cxx @@ -1,86 +1,103 @@ #include "samples/ISample.hxx" #include "event/FullEvent.hxx" #include "input/InputManager.hxx" #include "utility/FullEventUtility.hxx" #include "utility/ROOTUtility.hxx" #include "fhiclcpp/ParameterSet.h" +#include "generator/utility/NuWroUtility.hxx" + #include #include using namespace nuis::event; using namespace nuis::input; using namespace nuis::utility; +using namespace nuis::nuwrotools; class NuisToNuWro : public ISample { public: InputManager::Input_id_t fIH_id; std::unique_ptr fOutputTree; - event *fOutputEvent; + NuWroEvent *fOutputEvent; NuisToNuWro() : fIH_id(std::numeric_limits::max()), fOutputTree(nullptr), fOutputEvent(nullptr) {} void Initialize(fhicl::ParameterSet const &ps) { fIH_id = InputManager::Get().EnsureInputLoaded(ps); fOutputTree = MakeNewTTree(ps.get("output_file"), "treeout", "RECREATE"); fOutputTree->tree->Branch("e", &fOutputEvent); } void ProcessEvent(FullEvent const &ps) { + fOutputEvent->in.clear(); fOutputEvent->out.clear(); fOutputEvent->post.clear(); + std::pair NuMode = GetFlagsDynEquivalent(ps.mode); + + fOutputEvent->flag = NuMode.first; + fOutputEvent->dyn = NuMode.second; + for (Particle const &part : GetISParticles(ps)) { particle nuwro_part(part.pdg, part.P4.M()); nuwro_part.set_momentum(vect(part.P4.X(), part.P4.Y(), part.P4.Z())); fOutputEvent->in.push_back(nuwro_part); } for (Particle const &part : GetPrimaryFSParticles(ps)) { particle nuwro_part(part.pdg, part.P4.M()); nuwro_part.set_momentum(vect(part.P4.X(), part.P4.Y(), part.P4.Z())); fOutputEvent->out.push_back(nuwro_part); } for (Particle const &part : GetNuclearLeavingParticles(ps)) { particle nuwro_part(part.pdg, part.P4.M()); nuwro_part.set_momentum(vect(part.P4.X(), part.P4.Y(), part.P4.Z())); fOutputEvent->post.push_back(nuwro_part); } fOutputTree->tree->Fill(); } void ProcessSample(size_t nmax) { if (fIH_id == std::numeric_limits::max()) { throw uninitialized_ISample(); } IInputHandler const &IH = InputManager::Get().GetInputHandler(fIH_id); + + size_t NEvsToProcess = std::min(nmax, IH.GetNEvents()); + size_t NToShout = NEvsToProcess / 10; + std::cout << "[INFO]: Processing " << NEvsToProcess + << " input events to NuWro format." << std::endl; + size_t n = 0; for (auto const &fe : IH) { - if (++n > nmax) { + if (++n > NEvsToProcess) { break; } + if (NToShout && !(n % NToShout)) { + std::cout << "[INFO]: Processed " << n << "/" << NEvsToProcess + << " events." << std::endl; + } ProcessEvent(fe); } } - void Write() { - fOutputTree->file->Write(); - } + void Write() { fOutputTree->file->Write(); } std::string Name() { return "NuisToNuWro"; } }; DECLARE_PLUGIN(ISample, NuisToNuWro); diff --git a/src/samples/MCTools/VerboseEventSummary.cxx b/src/samples/MCTools/VerboseEventSummary.cxx index f56c7ee..726a5a5 100644 --- a/src/samples/MCTools/VerboseEventSummary.cxx +++ b/src/samples/MCTools/VerboseEventSummary.cxx @@ -1,60 +1,56 @@ #include "samples/ISample.hxx" #include "event/FullEvent.hxx" #include "input/InputManager.hxx" #include #include using namespace nuis::event; using namespace nuis::input; class VerboseEventSummary : public ISample { public: InputManager::Input_id_t fIH_id; VerboseEventSummary() : fIH_id(std::numeric_limits::max()) {} void Initialize(fhicl::ParameterSet const &ps) { fIH_id = InputManager::Get().EnsureInputLoaded(ps); } void ProcessEvent(FullEvent const &ps) { - std::cout << "Event: Interaction mode = " << ps.mode - << ", probe: { PDG: " << ps.probe_pdg - << ", Energy: " << ps.probe_E << " MeV }." << std::endl; - for (auto &status_stack : ps.ParticleStack) { - std::cout << "\t[" << status_stack.status << "]" << std::endl; - - for (Particle const &part : status_stack.particles) { - std::cout << "\t\t{ PDG: " << part.pdg << ", P3: [ " << part.P4[0] - << ", " << part.P4[1] << ", " << part.P4[2] - << "], E: " << part.P4[3] << ", M: " << part.P4.M() - << " }" << std::endl; - } - } - std::cout << std::endl; + std::cout << ps.to_string() << std::endl; } void ProcessSample(size_t nmax) { if (fIH_id == std::numeric_limits::max()) { throw uninitialized_ISample(); } IInputHandler const &IH = InputManager::Get().GetInputHandler(fIH_id); + size_t NEvsToProcess = std::min(nmax, IH.GetNEvents()); + size_t NToShout = NEvsToProcess / 10; + std::cout << "[INFO]: Processing " << NEvsToProcess << " input events." + << std::endl; + size_t n = 0; for (auto const &fe : IH) { if (++n > nmax) { break; } + if (NToShout && !(n % NToShout)) { + std::cout << "[INFO]: Processed " << n << "/" << NEvsToProcess + << " events." << std::endl; + } ProcessEvent(fe); } } void Write() {} std::string Name() { return "VerboseEventSummary"; } }; DECLARE_PLUGIN(ISample, VerboseEventSummary); diff --git a/src/samples/SimpleDataComparison.hxx b/src/samples/SimpleDataComparison.hxx new file mode 100644 index 0000000..15559bd --- /dev/null +++ b/src/samples/SimpleDataComparison.hxx @@ -0,0 +1,346 @@ +// 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_SIMPLEDATACOMPARISON_HXX_SEEN +#define SAMPLES_SIMPLEDATACOMPARISON_HXX_SEEN + +#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/ROOTUtility.hxx" + +#include +#include +#include +#include +#include + +template struct TH_dim_helper {}; +template <> struct TH_dim_helper<1> { typedef TH1D type; }; +template <> struct TH_dim_helper<2> { typedef TH2D type; }; + +template class SimpleDataComparison : public IDataComparison { + + NEW_NUIS_EXCEPT(invalid_SimpleDataComparison_initialization); + +protected: + 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; + + std::string fDataInputDescriptor; + std::unique_ptr::type> fData; + std::string fCovarianceInputDescriptor; + std::unique_ptr fCovariance; + std::unique_ptr::type> fPrediction; + std::unique_ptr::type> fPrediction_xsec; + std::unique_ptr::type> fPrediction_shape; + std::unique_ptr::type> fPrediction_comparison; + bool fComparisonFinalized; + + std::string fJournalReference; + std::string fTargetMaterial; + std::string fFluxDescription; + std::string fSignalDescription; + + std::pair EnuRange; + + std::function IsSigFunc; + std::function(nuis::event::FullEvent const &)> + CompProjFunc; + +public: + SimpleDataComparison() { + fIH_id = std::numeric_limits::max(); + write_directory = ""; + NMaxSample_override = std::numeric_limits::max(); + fDataInputDescriptor = ""; + fData = 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 (double &el : arr) { + el = 0; + } + return arr; + }; + + fJournalReference = ""; + fTargetMaterial = ""; + fFluxDescription = ""; + fSignalDescription = ""; + fIsShapeOnly = -1; + fIsFluxUnfolded = -1; + EnuRange = std::pair{std::numeric_limits::max(), + std::numeric_limits::max()}; + } + + void ReadGlobalConfigDefaults() { + fhicl::ParameterSet const &global_sample_configuration = + nuis::config::GetDocument().get( + std::string("global.sample_configuration.") + Name(), + fhicl::ParameterSet()); + + if (!fJournalReference.length()) { + fJournalReference = global_sample_configuration.get( + "journal_reference", fJournalReference); + } + if (!fTargetMaterial.length()) { + fTargetMaterial = global_sample_configuration.get( + "target_material", fTargetMaterial); + } + if (!fFluxDescription.length()) { + fFluxDescription = global_sample_configuration.get( + "flux_description", fFluxDescription); + } + if (!fSignalDescription.length()) { + fSignalDescription = global_sample_configuration.get( + "signal_description", fSignalDescription); + } + + if (fIsShapeOnly == -1) { + fIsShapeOnly = global_sample_configuration.get("shape_only", false); + } + if (fIsFluxUnfolded == -1) { + fIsFluxUnfolded = + global_sample_configuration.get("flux_unfolded", false); + } + + if ((EnuRange.first == std::numeric_limits::max()) && + (global_sample_configuration.has_key("enu_range"))) { + EnuRange = global_sample_configuration.get>( + "enu_range"); + } + } + + virtual std::string GetJournalReference() { return fJournalReference; } + 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 SetCovariance(std::string const &cov_descriptor) { + fCovarianceInputDescriptor = cov_descriptor; + } + + virtual void FillProjection(std::array const &proj, + double event_weight) { + nuis::utility::Fill(fPrediction.get(), proj, event_weight); + } + + virtual void FinalizeComparison() { + fPrediction_xsec = nuis::utility::Clone(fPrediction); + fPrediction_xsec->Scale(1.0, "width"); + fPrediction_shape = nuis::utility::Clone(fPrediction_xsec); + if (fData) { + fPrediction_shape->Scale(fData->Integral() / + fPrediction_shape->Integral()); + } else { + ISAMPLE_WARN("When Finalizing comparison, no Data histogram available."); + } + + if (fIsFluxUnfolded) { + // fPrediction_comparison + } else if (fIsShapeOnly) { + fPrediction_comparison = nuis::utility::Clone(fPrediction_shape); + } else { + fPrediction_comparison = nuis::utility::Clone(fPrediction_xsec); + } + fComparisonFinalized = true; + } + + void Initialize(fhicl::ParameterSet const &ps) { + + if (ps.has_key("verbosity")) { + SetSampleVerbosity(ps.get("verbosity")); + } else { + SetSampleVerbosity("Reticent"); + } + + ReadGlobalConfigDefaults(); + + fhicl::ParameterSet const &global_sample_configuration = + nuis::config::GetDocument().get( + std::string("global.sample_configuration.") + Name(), + fhicl::ParameterSet()); + + if (ps.has_key("fake_data")) { + fData = nuis::utility::GetHistogram::type>( + ps.get("fake_data_histogram")); + } else { + if (!fDataInputDescriptor.length()) { + if (!global_sample_configuration.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 = + global_sample_configuration.get("data_descriptor"); + } + fData = nuis::utility::GetHistogram::type>( + fDataInputDescriptor); + } + + fPrediction = nuis::utility::Clone(fData, true); + + if (fCovarianceInputDescriptor.length()) { + fCovariance = + nuis::utility::GetHistogram(fCovarianceInputDescriptor); + } else if (global_sample_configuration.has_key("covariance_descriptor")) { + fCovariance = nuis::utility::GetHistogram( + global_sample_configuration.get( + "covariance_descriptor")); + } + + if (ps.has_key("verbosity")) { + SetSampleVerbosity(ps.get("verbosity")); + } + + NMaxSample_override = + ps.get("nmax", std::numeric_limits::max()); + + write_directory = ps.get("write_directory", ""); + + fIH_id = nuis::input::InputManager::Get().EnsureInputLoaded(ps); + } + void ProcessSample(size_t nmax = std::numeric_limits::max()) { + if (fIH_id == + std::numeric_limits::max()) { + throw uninitialized_ISample(); + } + 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; + ISAMPLE_INFO("Sample " << std::quoted(Name()) << " processing " + << NEvsToProcess << " events."); + + IInputHandler::ev_index_t ev_idx = 0; + size_t NSigEvents = 0; + + bool DetermineSignalEvents = !fSignalCache.size(); + + nuis::utility::Clear(fPrediction.get()); + fComparisonFinalized = false; + + while (ev_idx < NEvsToProcess) { + if (DetermineSignalEvents) { + nuis::event::FullEvent const &fev = IH.GetFullEvent(ev_idx); + bool is_sig = IsSigFunc(fev); + fSignalCache.push_back(is_sig); + if (is_sig) { + fProjectionCache.push_back(CompProjFunc(fev)); + } + } + + if (NToShout && !(ev_idx % NToShout)) { + ISAMPLE_INFO("\tProcessed " << ev_idx << "/" << NEvsToProcess + << " events."); + } + + if (fSignalCache[ev_idx]) { + FillProjection(fProjectionCache[ev_idx], + IH.GetEventWeight(ev_idx) * nmax_scaling); + } + + ev_idx++; + } + ISAMPLE_INFO("\t" << fProjectionCache.size() << "/" << NEvsToProcess + << " events passed selection."); + } + void Write() { + if (!fComparisonFinalized) { + FinalizeComparison(); + } + + nuis::persistency::WriteToOutputFile::type>( + fPrediction_comparison.get(), "Prediction", write_directory); + nuis::persistency::WriteToOutputFile::type>( + fPrediction_xsec.get(), "Prediction_xsec", write_directory); + + if (fData) { + nuis::persistency::WriteToOutputFile::type>( + fData.get(), "Data", write_directory); + nuis::persistency::WriteToOutputFile::type>( + fPrediction_shape.get(), "Prediction_shape", write_directory); + } + } + + double GetGOF() { return 1; } + double GetNDOGuess() { + if (fData) { + return nuis::utility::TH_traits< + typename TH_dim_helper::type>::NbinsIncludeFlow(fData.get()); + } + return 0; + } + + fhicl::ParameterSet GetExampleConfiguration() { + fhicl::ParameterSet exps; + exps.put("name", Name()); + exps.put("input_type", "Generator"); + exps.put("file", "Events.root"); + exps.put("write_directory", Name()); + exps.put("fake_data_histogram", "/path/to/file.root;h_name"); + + return exps; + } +}; + +typedef SimpleDataComparison<1> SimpleDataComparison_1D; +typedef SimpleDataComparison<2> SimpleDataComparison_2D; + +#endif diff --git a/src/samples/nuA/BubbleChamber/ANL/ANL.sample.config.fcl b/src/samples/nuA/BubbleChamber/ANL/ANL.sample.config.fcl new file mode 100644 index 0000000..80c8e58 --- /dev/null +++ b/src/samples/nuA/BubbleChamber/ANL/ANL.sample.config.fcl @@ -0,0 +1,9 @@ +global.sample_configuration.ANL_CCQE_Evt_1DQ2_nu: { + journal_reference: "PRL 31 844 / PRD 16 3103 / PRD 26 537" + target_material: "D2" + flux_description: "ANL Muon Neutrino" + signal_description: "True CCQE" + shape_only: true + flux_unfolded: false + use_D2_correction: false +} 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 07504b1..6fa331f 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,236 +1,158 @@ // 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/IDataComparison.hxx" +#include "samples/SimpleDataComparison.hxx" -#include "event/FullEvent.hxx" - -#include "input/InputManager.hxx" - -#include "persistency/ROOTOutput.hxx" - -#include "utility/FileSystemUtility.hxx" #include "utility/FullEventUtility.hxx" #include "utility/KinematicUtility.hxx" #include "utility/PDGCodeUtility.hxx" -#include "utility/ROOTUtility.hxx" - -#include "TH1D.h" using namespace nuis::event; -using namespace nuis::input; using namespace nuis::utility; -using namespace nuis::persistency; -class ANL_CCQE_Evt_1DQ2_nu : public IDataComparison { +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; - InputManager::Input_id_t fIH_id; - std::string write_directory; - - std::unique_ptr fData; - std::unique_ptr fPrediction; - - std::pair EnuRange; - - ANL_CCQE_Evt_1DQ2_nu() - : Pub(kPRD26), Pub_str(""), UseD2Corr(false), - fIH_id(std::numeric_limits::max()), - write_directory("ANL_CCQE_Evt_1DQ2_nu"), fData(nullptr), - fPrediction(nullptr) {} - - std::string GetJournalReference() { - return "PRL 31 844 / PRD 16 3103 / PRD 26 537"; + ANL_CCQE_Evt_1DQ2_nu() : Pub(kPRD26), Pub_str(""), UseD2Corr(false) { + ReadGlobalConfigDefaults(); } - std::string GetTargetMaterial() { return "D2"; } - std::string GetFluxDescription() { return "ANL Muon Neutrino"; } - std::string GetSignalDescription() { return "True CCQE"; } + 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; + fhicl::ParameterSet exps = + SimpleDataComparison_1D::GetExampleConfiguration(); - exps.put("name", "ANL_CCQE_Evt_1DQ2_nu"); - exps.put("input_type", "Generator"); - exps.put("file", "ANL_Events.root"); - exps.put("write_directory", "ANL_CCQE_Evt_1DQ2_nu_Generator"); exps.put("publication", "PRD26"); - exps.put("use_D2_correction", true); - - fhicl::ParameterSet fd; - fd.put("file", "ANL_fake_data.root"); - fd.put("histogram_name", "fake_data"); - - exps.put("fake_data", fd); + exps.put("use_D2_correction", false); return exps; } void Initialize(fhicl::ParameterSet const &ps) { - if (!ps.has_key("publication")) { - std::string publication = ps.get("publication"); - 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 ]"; - } + if (ps.has_key("verbosity")) { + SetSampleVerbosity(ps.get("verbosity")); + } + + std::string publication = ps.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"; - EnuRange = std::pair{0.0, 3.0}; + EnuRange = std::pair{0, 3E3}; + ISAMPLE_INFO("Sample " << Name() + << " specialized for publication: " << Pub_str); break; } case kPRD16: { Pub_str = "PRD16_3103"; - EnuRange = std::pair{0.0, 6.0}; + EnuRange = std::pair{0, 6E3}; + ISAMPLE_INFO("Sample " << Name() + << " specialized for publication: " << Pub_str); break; } case kPRD26: { Pub_str = "PRD26_537"; - EnuRange = std::pair{0.0, 6.0}; + EnuRange = std::pair{0, 6E3}; + ISAMPLE_INFO("Sample " << Name() + << " specialized for publication: " << Pub_str); break; } } - if (ps.has_key("use_D2_correction")) { - UseD2Corr = ps.get("use_D2_correction"); - } - - if (ps.has_key("use_D2_correction")) { - UseD2Corr = ps.get("use_D2_correction"); - } - - if (ps.has_key("write_directory")) { - write_directory = ps.get("write_directory"); - } - - fIH_id = InputManager::Get().EnsureInputLoaded(ps); - - if (ps.has_key("fake_data")) { - fhicl::ParameterSet const &fd = ps.get("fake_data"); - fData = GetHistogramFromROOTFile( - fd.get("file"), fd.get("histogram_name")); - } else { - fData = GetHistogramFromROOTFile( - GetDataDir() + "nuA/BubbleChamber/ANL/ANL_CCQE_Data_" + Pub_str + - ".root", - "ANL_1DQ2_Data"); - } - fPrediction = CloneHistogram(fData, true); - } - - std::vector fIsSignal; - std::vector fQ2; + fhicl::ParameterSet const &global_sample_configuration = + nuis::config::GetDocument().get( + std::string("global.sample_configuration.") + Name(), + fhicl::ParameterSet()); - void ProcessEvent(FullEvent const &fev) { - fQ2.push_back(GetNeutrinoQ2QERec(fev)); - } + UseD2Corr = ps.get( + "use_D2_correction", + global_sample_configuration.get("use_D2_correction", false)); - bool IsSignal(FullEvent const &fev) { + SetData(GetDataDir() + "nuA/BubbleChamber/ANL/CCQE/ANL_CCQE_Data_" + + Pub_str + ".root;ANL_1DQ2_Data"); - if (fev.mode != Channel_t::kCCQE) { - return false; - } - - Particle ISNumu = GetISNeutralLepton(fev); + SimpleDataComparison_1D::Initialize(ps); - if (!ISNumu) { - return false; - } - - if (ISNumu.pdg != pdgcodes::kNuMu) { - return false; - } + // Signal selection function + IsSigFunc = [&](FullEvent const &fev) -> bool { + if (fev.mode != Channel_t::kCCQE) { + return false; + } - if ((ISNumu.P4.E() < EnuRange.first) || (ISNumu.P4.E() > EnuRange.second)) { - return false; - } + Particle ISNumu = GetHMISNeutralLepton(fev); - double Q2 = GetNeutrinoQ2QERec(fev); - if (Q2 <= 0) { - return false; - } + if (!ISNumu) { + return false; + } - return true; - } + if (ISNumu.pdg != pdgcodes::kNuMu) { + return false; + } - void ProcessSample(size_t nmax) { - if (fIH_id == std::numeric_limits::max()) { - throw uninitialized_ISample(); - } - IInputHandler const &IH = InputManager::Get().GetInputHandler(fIH_id); - - size_t NEvsToProcess = std::min(nmax, IH.GetNEvents()); - IInputHandler::ev_index_t ev_idx = 0; - size_t NSigEvents = 0; - - bool DetermineSignalEvents = !fIsSignal.size(); - - while (ev_idx < NEvsToProcess) { - if (DetermineSignalEvents) { - FullEvent const &fev = IH.GetFullEvent(ev_idx); - bool is_sig = IsSignal(fev); - fIsSignal.push_back(is_sig); - if (is_sig) { - ProcessEvent(fev); - } + if ((ISNumu.P4.E() < EnuRange.first) || + (ISNumu.P4.E() > EnuRange.second)) { + return false; } - if (fIsSignal[ev_idx]) { - fPrediction->Fill(fQ2[NSigEvents++], IH.GetEventWeight(ev_idx)); + double Q2 = GetNeutrinoQ2QERec(fev, 0); + if (Q2 <= 0) { + return false; } - ev_idx++; - } + return true; + }; + // 1D Projection function + CompProjFunc = [](FullEvent const &fev) -> std::array { + return {{GetNeutrinoQ2QERec(fev, 0)}}; + }; } - void Write() { - WriteToOutputFile(fData.get(), "Data", write_directory); - WriteToOutputFile(fPrediction.get(), "Prediction", write_directory); - } std::string Name() { return "ANL_CCQE_Evt_1DQ2_nu"; } - - double GetGOF() { return 0; /*CalcChi2(fData, fPrediction);*/ } }; DECLARE_PLUGIN(IDataComparison, ANL_CCQE_Evt_1DQ2_nu); +DECLARE_PLUGIN(ISample, ANL_CCQE_Evt_1DQ2_nu); diff --git a/src/samples/nuA/BubbleChamber/ANL/CMakeLists.txt b/src/samples/nuA/BubbleChamber/ANL/CMakeLists.txt index 0aa6e75..819d807 100644 --- a/src/samples/nuA/BubbleChamber/ANL/CMakeLists.txt +++ b/src/samples/nuA/BubbleChamber/ANL/CMakeLists.txt @@ -1,14 +1,17 @@ SET(SAMPLES ANL_CCQE_Evt_1DQ2_nu) -LIST(APPEND INuADataComparisons ANL_CCQE_Evt_1DQ2_nu) - foreach(S ${SAMPLES}) + LIST(APPEND INuADataComparisons ${S}) LIST(APPEND SAMPLES_src ${S}.cxx) endforeach() add_library(ANLDataComparisons SHARED ${SAMPLES_src}) target_link_libraries(ANLDataComparisons nuis_event nuis_input nuis_config nuis_persistency) install(TARGETS ANLDataComparisons DESTINATION plugins) +LIST(APPEND INuADataComparisons_FHiCL ANL.sample.config.fcl) +install(FILES ANL.sample.config.fcl DESTINATION fcl) + SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) +SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE) diff --git a/src/samples/nuA/BubbleChamber/CMakeLists.txt b/src/samples/nuA/BubbleChamber/CMakeLists.txt index 0c841fc..9c48c97 100644 --- a/src/samples/nuA/BubbleChamber/CMakeLists.txt +++ b/src/samples/nuA/BubbleChamber/CMakeLists.txt @@ -1,3 +1,4 @@ add_subdirectory(ANL) SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) +SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE) diff --git a/src/samples/nuA/CMakeLists.txt b/src/samples/nuA/CMakeLists.txt index 836e6d7..36dcced 100644 --- a/src/samples/nuA/CMakeLists.txt +++ b/src/samples/nuA/CMakeLists.txt @@ -1,3 +1,4 @@ add_subdirectory(BubbleChamber) SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE) +SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE) diff --git a/src/utility/CMakeLists.txt b/src/utility/CMakeLists.txt index 88ddaf6..7b26067 100644 --- a/src/utility/CMakeLists.txt +++ b/src/utility/CMakeLists.txt @@ -1,19 +1,21 @@ set(utility_implementation_files FileSystemUtility.cxx FullEventUtility.cxx - KinematicUtility.cxx) + KinematicUtility.cxx + PDGCodeUtility.cxx) set(utility_header_files FileSystemUtility.hxx FullEventUtility.hxx InteractionChannelUtility.hxx KinematicUtility.hxx PDGCodeUtility.hxx - ROOTUtility.hxx) + ROOTUtility.hxx + HistogramUtility.hxx) 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) diff --git a/src/utility/FullEventUtility.cxx b/src/utility/FullEventUtility.cxx index a6be49a..8188289 100644 --- a/src/utility/FullEventUtility.cxx +++ b/src/utility/FullEventUtility.cxx @@ -1,147 +1,147 @@ #include "utility/FullEventUtility.hxx" #include "event/FullEvent.hxx" #include "utility/PDGCodeUtility.hxx" using namespace nuis::event; namespace nuis { namespace utility { 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.P4.Vect().Mag() > HMParticle.P4.Vect().Mag()) { HMParticle = part; } } } return HMParticle; } 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 GetFSChargedLepton(FullEvent const &ev) { +Particle GetHMFSChargedLepton(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::ChargedLeptons); } -Particle GetFSNeutralLepton(FullEvent const &ev) { +Particle GetHMFSNeutralLepton(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::NeutralLeptons); } -Particle GetISNeutralLepton(FullEvent const &ev) { +Particle GetHMISNeutralLepton(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::NeutralLeptons, Particle::Status_t::kPrimaryInitialState); } -Particle GetFSChargedPion(FullEvent const &ev) { +Particle GetHMFSChargedPion(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::ChargedPions); } -Particle GetFSNeutralPion(FullEvent const &ev) { +Particle GetHMFSNeutralPion(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::NeutralPions); } -Particle GetFSPion(FullEvent const &ev) { +Particle GetHMFSPion(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Pions); } -Particle GetFSProton(FullEvent const &ev) { +Particle GetHMFSProton(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Protons); } -Particle GetFSNeutron(FullEvent const &ev) { +Particle GetHMFSNeutron(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Neutron); } -Particle GetFSNucleon(FullEvent const &ev) { +Particle GetHMFSNucleon(FullEvent const &ev) { return GetHMParticle(ev, pdgcodes::Nucleons); } -Particle GetFSOther(FullEvent const &ev) { +Particle GetHMFSOther(FullEvent const &ev) { return GetHMParticle(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 f497ea4..7d11a96 100644 --- a/src/utility/FullEventUtility.hxx +++ b/src/utility/FullEventUtility.hxx @@ -1,80 +1,80 @@ // 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_FULLEVENTUTILITY_HXX_SEEN #define UTILITY_FULLEVENTUTILITY_HXX_SEEN #include "event/types.hxx" #include "event/Particle.hxx" #include namespace nuis { namespace event { class FullEvent; } // namespace event } // namespace nuis namespace nuis { namespace utility { 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); 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 GetFSChargedLepton(event::FullEvent const &); -event::Particle GetFSNeutralLepton(event::FullEvent const &); -event::Particle GetISNeutralLepton(event::FullEvent const &); -event::Particle GetFSChargedPion(event::FullEvent const &); -event::Particle GetFSNeutralPion(event::FullEvent const &); -event::Particle GetFSPion(event::FullEvent const &); -event::Particle GetFSProton(event::FullEvent const &); -event::Particle GetFSNeutron(event::FullEvent const &); -event::Particle GetFSNucleon(event::FullEvent const &); -event::Particle GetFSOther(event::FullEvent const &); +event::Particle GetHMFSChargedLepton(event::FullEvent const &); +event::Particle GetHMFSNeutralLepton(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 &); } // namespace utility } // namespace nuis #endif diff --git a/src/utility/HistogramUtility.hxx b/src/utility/HistogramUtility.hxx new file mode 100644 index 0000000..34e99bc --- /dev/null +++ b/src/utility/HistogramUtility.hxx @@ -0,0 +1,61 @@ +// 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_HISTOGRAMUTILITY_HXX_SEEN +#define UTILITY_HISTOGRAMUTILITY_HXX_SEEN + +#include "utility/ROOTUtility.hxx" + +#include "exception/exception.hxx" + +#include "string_parsers/from_string.hxx" + +#include +#include +#include +#include + +namespace nuis { +namespace utility { + +NEW_NUIS_EXCEPT(unimplemented_GetHistogram_method); +NEW_NUIS_EXCEPT(invalid_histogram_descriptor); + +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)."; + } +} + +} // namespace utility +} // namespace nuis +#endif diff --git a/src/utility/InteractionChannelUtility.hxx b/src/utility/InteractionChannelUtility.hxx index 340c565..54b1c78 100644 --- a/src/utility/InteractionChannelUtility.hxx +++ b/src/utility/InteractionChannelUtility.hxx @@ -1,64 +1,66 @@ #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)); } } // namespace utility } // namespace nuis #endif diff --git a/src/utility/KinematicUtility.cxx b/src/utility/KinematicUtility.cxx index 711c65c..8337a7f 100644 --- a/src/utility/KinematicUtility.cxx +++ b/src/utility/KinematicUtility.cxx @@ -1,11 +1,86 @@ #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 GetNeutrinoQ2QERec(event::FullEvent const &fev){ - return 0; +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) { + 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 Theta_nu_mu = neutrino.P3().Angle(charged_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) { + 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 Theta_nu_mu = neutrino.P3().Angle(charged_lepton.P3()); + + return -ml_GeV * ml_GeV + 2.0 * GetNeutrinoEQERec(fev, SeparationEnergy_MeV) * + (el_GeV - pl_GeV * cos(Theta_nu_mu)); } } // namespace utility } // namespace nuis diff --git a/src/utility/KinematicUtility.hxx b/src/utility/KinematicUtility.hxx index 5b3bc73..f9e698d 100644 --- a/src/utility/KinematicUtility.hxx +++ b/src/utility/KinematicUtility.hxx @@ -1,34 +1,40 @@ // 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 "event/FullEvent.hxx" -#include "event/Particle.hxx" +namespace nuis { +namespace event { +class FullEvent; +} +} // namespace nuis namespace nuis { namespace utility { -double GetNeutrinoQ2QERec(event::FullEvent const &fev); +double GetNeutrinoEQERec(event::FullEvent const &fev, + double SeparationEnergy_MeV); +double GetNeutrinoQ2QERec(event::FullEvent const &fev, + double SeparationEnergy_MeV); } // namespace utility } // namespace nuis #endif diff --git a/src/utility/PDGCodeUtility.cxx b/src/utility/PDGCodeUtility.cxx new file mode 100644 index 0000000..c5fc138 --- /dev/null +++ b/src/utility/PDGCodeUtility.cxx @@ -0,0 +1,116 @@ +#include "utility/PDGCodeUtility.hxx" + +#include +#include + +using namespace nuis::event; +using namespace nuis::utility::pdgcodes; +using namespace nuis::utility::pdgmasses; + +namespace nuis { +namespace utility { + +static std::map const PDGMasses{ + {kNuMu, kNuMuMass_MeV}, {kNuMuBar, kNuMuBarMass_MeV}, + {kMu, kMuMass_MeV}, {kMuPlus, kMuPlusMass_MeV}, + {kNue, kNueMass_MeV}, {kNueBar, kNueBarMass_MeV}, + {kPiPlus, kPiPlusMass_MeV}, {kPiMinus, kPiMinusMass_MeV}, + {kPi0, kPi0Mass_MeV}, {kProton, kNeutronMass_MeV}, + {kNeutron, kProtonMass_MeV}}; + +double GetPDGMass(PDG_t pdg) { + if (PDGMasses.find(pdg) == PDGMasses.end()) { + throw unhandled_pdg_code() + << "[ERROR]: Unknown mass for particle with PDG code: " << pdg + << ", please add it to src/utility/PDGCodeUtility.cxx:PDGMasses"; + } + return PDGMasses.at(pdg); +} + +NEW_NUIS_EXCEPT(invalid_MatterType); + +bool IsInPDGList(PDG_t pdg, std::vector const &MatterList, + std::vector const &AntiMatterList, MatterType type) { + switch (type) { + case kMatter: { + return std::count(MatterList.begin(), MatterList.end(), pdg); + } + case kMatterAntimatter: { + return std::count(MatterList.begin(), MatterList.end(), pdg) || + std::count(AntiMatterList.begin(), AntiMatterList.end(), pdg); + } + case kAntimatter: { + return std::count(AntiMatterList.begin(), AntiMatterList.end(), pdg); + } + default: { throw invalid_MatterType(); } + } +} +bool IsNeutralLepton(PDG_t pdg, MatterType type) { + return IsInPDGList(pdg, NeutralLeptons_matter, NeutralLeptons_antimatter, + type); +} +bool IsChargedLepton(PDG_t pdg, MatterType type) { + return IsInPDGList(pdg, ChargedLeptons_matter, ChargedLeptons_antimatter, + type); +} +bool IsProton(PDG_t pdg, MatterType type) { + return IsInPDGList(pdg, Proton_matter, Proton_antimatter, type); +} +bool IsNeutron(PDG_t pdg) { return pdg == Neutron[0]; } +bool IsChargedPion(PDG_t pdg) { + return std::count(ChargedPions.begin(), ChargedPions.end(), pdg); +} +bool IsNeutralPion(PDG_t pdg) { + return std::count(NeutralPions.begin(), NeutralPions.end(), pdg); +} +bool IsPion(PDG_t pdg) { return std::count(Pions.begin(), Pions.end(), pdg); } +bool IsOther(PDG_t pdg) { + return !std::count(CommonParticles.begin(), CommonParticles.end(), pdg); +} + +bool IsMatter(PDG_t pdg) { + // Special cases + switch (pdg) { + case kPiPlus: { + return true; + } + case kPiMinus: { + return true; + } + case kPi0: { + return true; + } + case kNeutron: { + return true; + } + } + return (pdg > 0); +} +bool IsAntiMatter(PDG_t pdg) { + // Special cases + switch (pdg) { + case kPiPlus: { + return true; + } + case kPiMinus: { + return true; + } + case kPi0: { + return true; + } + case kNeutron: { + return true; + } + } + return (pdg < 0); +} + +PDG_t MakeNuclearPDG(size_t A, size_t Z) { + return 1000 * Z + 10 * A + 1000000000; +} + +size_t GetA(PDG_t pdg) { return ((pdg / 10) % 1000); } +size_t GetZ(PDG_t pdg) { return ((pdg / 1000) % 1000); } + +} // namespace utility +} // namespace nuis diff --git a/src/utility/PDGCodeUtility.hxx b/src/utility/PDGCodeUtility.hxx index 711f575..0cb2f30 100644 --- a/src/utility/PDGCodeUtility.hxx +++ b/src/utility/PDGCodeUtility.hxx @@ -1,135 +1,127 @@ // 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_PDGCODEUTILITY_HXX_SEEN #define UTILITY_PDGCODEUTILITY_HXX_SEEN -#include "core/types.hxx" +#include "event/types.hxx" + +#include "exception/exception.hxx" + +#include namespace nuis { namespace utility { namespace pdgcodes { +NEW_NUIS_EXCEPT(unhandled_pdg_code); + enum MatterType { kMatter = 1, kMatterAntimatter = 0, kAntimatter = -1 }; -static core::PDG_t const kNuMu = 14; -static core::PDG_t const kNuMuBar = -14; -static core::PDG_t const kMu = 13; -static core::PDG_t const kMuPlus = -13; +static event::PDG_t const kNuMu = 14; +static event::PDG_t const kNuMuBar = -14; + +static event::PDG_t const kMu = 13; +static event::PDG_t const kMuPlus = -13; -static core::PDG_t const kNue = 12; -static core::PDG_t const kNueBar = 12; +static event::PDG_t const kNue = 12; +static event::PDG_t const kNueBar = 12; -static core::PDG_t const kPiPlus = 211; -static core::PDG_t const kPiMinus = -211; -static core::PDG_t const kPi0 = 111; +static event::PDG_t const kPiPlus = 211; +static event::PDG_t const kPiMinus = -211; +static event::PDG_t const kPi0 = 111; -static core::PDG_t const kProton = 2212; -static core::PDG_t const kNeutron = 2112; +static event::PDG_t const kProton = 2212; +static event::PDG_t const kNeutron = 2112; -static std::vector const ChargedLeptons{11, kMu, 15, +static std::vector const ChargedLeptons{11, kMu, 15, -11, kMuPlus, -15}; -static std::vector const ChargedLeptons_matter{11, kMu, 15}; -static std::vector const ChargedLeptons_antimatter{-11, kMuPlus, +static std::vector const ChargedLeptons_matter{11, kMu, 15}; +static std::vector const ChargedLeptons_antimatter{-11, kMuPlus, -15}; -static std::vector const NeutralLeptons{kNue, kNuMu, 16, +static std::vector const NeutralLeptons{kNue, kNuMu, 16, kNueBar, kNuMuBar, -16}; -static std::vector const NeutralLeptons_matter{kNue, kNuMu, 16}; -static std::vector const NeutralLeptons_antimatter{kNueBar, +static std::vector const NeutralLeptons_matter{kNue, kNuMu, 16}; +static std::vector const NeutralLeptons_antimatter{kNueBar, kNuMuBar, -16}; -static std::vector const ChargedPions{kPiPlus, kPiMinus}; -static std::vector const NeutralPions{kPi0}; -static std::vector const Pions{kPiPlus, kPiMinus, kPi0}; -static std::vector const Protons{2212, -2122}; -static std::vector const Proton_matter{2212}; -static std::vector const Proton_antimatter{-2212}; -static std::vector const Neutron{kNeutron}; -static std::vector const Nucleons{kProton, 2112, -2212}; -static std::vector const Nucleons_matter{kProton, kNeutron}; -static std::vector const Nucleons_antimatter{-2212, 2112}; - -static std::vector const CommonParticles{ +static std::vector const ChargedPions{kPiPlus, kPiMinus}; +static std::vector const NeutralPions{kPi0}; +static std::vector const Pions{kPiPlus, kPiMinus, kPi0}; +static std::vector const Protons{kProton, -kProton}; +static std::vector const Proton_matter{kProton}; +static std::vector const Proton_antimatter{-kProton}; +static std::vector const Neutron{kNeutron}; +static std::vector const Nucleons{kProton, kNeutron, -kProton}; +static std::vector const Nucleons_matter{kProton, kNeutron}; +static std::vector const Nucleons_antimatter{-kProton, kNeutron}; + +static std::vector const CommonParticles{ 11, kMu, 15, -11, kMuPlus, -15, kNue, kNuMu, 16, kNueBar, kNuMuBar, -16, kPiPlus, kPiMinus, kPi0, kProton, kNeutron}; } // namespace pdgcodes -inline bool -IsInPDGList(core::PDG_t pdg, std::vector const &MatterList, - std::vector const &AntiMatterList, - pdgcodes::MatterType type = pdgcodes::kMatterAntimatter) { - switch (type) { - case pdgcodes::kMatter: { - return std::count(MatterList.begin(), MatterList.end(), pdg); - } - case pdgcodes::kMatterAntimatter: { - return std::count(MatterList.begin(), MatterList.end(), pdg) || - std::count(AntiMatterList.begin(), AntiMatterList.end(), pdg); - } - case pdgcodes::kAntimatter: { - return std::count(AntiMatterList.begin(), AntiMatterList.end(), pdg); - } - } -} -inline bool -IsNeutralLepton(core::PDG_t pdg, - pdgcodes::MatterType type = pdgcodes::kMatterAntimatter) { - return IsInPDGList(pdg, pdgcodes::NeutralLeptons_matter, - pdgcodes::NeutralLeptons_antimatter, type); -} -inline bool -IsChargedLepton(core::PDG_t pdg, - pdgcodes::MatterType type = pdgcodes::kMatterAntimatter) { - return IsInPDGList(pdg, pdgcodes::ChargedLeptons_matter, - pdgcodes::ChargedLeptons_antimatter, type); -} -inline bool IsProton(core::PDG_t pdg, - pdgcodes::MatterType type = pdgcodes::kMatterAntimatter) { - return IsInPDGList(pdg, pdgcodes::Proton_matter, pdgcodes::Proton_antimatter, - type); -} -inline bool IsNeutron(core::PDG_t pdg) { return pdg == pdgcodes::Neutron[0]; } -inline bool IsChargedPion(core::PDG_t pdg) { - return std::count(pdgcodes::ChargedPions.begin(), - pdgcodes::ChargedPions.end(), pdg); -} -inline bool IsNeutralPion(core::PDG_t pdg) { - return std::count(pdgcodes::NeutralPions.begin(), - pdgcodes::NeutralPions.end(), pdg); -} -inline bool IsPion(core::PDG_t pdg) { - return std::count(pdgcodes::Pions.begin(), pdgcodes::Pions.end(), pdg); -} -inline bool IsOther(core::PDG_t pdg) { - return !std::count(pdgcodes::CommonParticles.begin(), - pdgcodes::CommonParticles.end(), pdg); -} - -inline core::PDG_t MakeNuclearPDG(size_t A, size_t Z) { - return 1000 * Z + 10 * A + 1000000000; -} - -inline size_t GetA(core::PDG_t pdg) { return ((pdg / 10) % 1000); } -inline size_t GetZ(core::PDG_t pdg) { return ((pdg / 1000) % 1000); } +namespace pdgmasses { +static double const kNuMuMass_MeV = 0; +static double const kNuMuBarMass_MeV = 0; + +static double const kMuMass_MeV = 105.65; +static double const kMuPlusMass_MeV = 105.65; + +static double const kNueMass_MeV = 0; +static double const kNueBarMass_MeV = 0; + +static double const kPiPlusMass_MeV = 139.57; +static double const kPiMinusMass_MeV = 139.57; +static double const kPi0Mass_MeV = 134.97; + +static double const kNeutronMass_MeV = 939.56; +static double const kProtonMass_MeV = 938.27; +} // namespace pdgmasses + +double GetPDGMass(event::PDG_t); + +bool IsInPDGList(event::PDG_t pdg, std::vector const &MatterList, + std::vector const &AntiMatterList, + pdgcodes::MatterType type = pdgcodes::kMatterAntimatter); +bool IsNeutralLepton(event::PDG_t pdg, + pdgcodes::MatterType type = pdgcodes::kMatterAntimatter); +bool IsChargedLepton(event::PDG_t pdg, + pdgcodes::MatterType type = pdgcodes::kMatterAntimatter); +bool IsProton(event::PDG_t pdg, + pdgcodes::MatterType type = pdgcodes::kMatterAntimatter); +bool IsNeutron(event::PDG_t pdg); +bool IsChargedPion(event::PDG_t pdg); +bool IsNeutralPion(event::PDG_t pdg); +bool IsPion(event::PDG_t pdg); +bool IsOther(event::PDG_t pdg); + +bool IsMatter(event::PDG_t pdg); +bool IsAntiMatter(event::PDG_t pdg); + +event::PDG_t MakeNuclearPDG(size_t A, size_t Z); + +size_t GetA(event::PDG_t); +size_t GetZ(event::PDG_t); } // namespace utility } // namespace nuis #endif diff --git a/src/utility/ROOTUtility.hxx b/src/utility/ROOTUtility.hxx index f738c7d..34231f3 100644 --- a/src/utility/ROOTUtility.hxx +++ b/src/utility/ROOTUtility.hxx @@ -1,166 +1,211 @@ // 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_ROOTUTILITY_HXX_SEEN #define UTILITY_ROOTUTILITY_HXX_SEEN +#include "exception/exception.hxx" + #include "TFile.h" +#include "TH1D.h" +#include "TH1F.h" +#include "TH2D.h" +#include "TH2F.h" #include "TTree.h" -#include "exception/exception.hxx" - #include #include #include #include -class TH1; -class TH1D; -class TH1F; -class TH2; -class TH2D; -class TH2F; - namespace nuis { namespace utility { NEW_NUIS_EXCEPT(failed_to_open_TFile); NEW_NUIS_EXCEPT(failed_to_get_TTree); NEW_NUIS_EXCEPT(invalid_histogram_name); NEW_NUIS_EXCEPT(failed_to_clone); inline TFile *CheckOpenTFile(std::string const &fname, char const *opts = "") { TDirectory *ogDir = gDirectory; TFile *inpF = new TFile(fname.c_str(), opts); if (!inpF || !inpF->IsOpen()) { throw failed_to_open_TFile() << "[ERROR]: Couldn't open input file: " << std::quoted(fname); } if (ogDir) { ogDir->cd(); } return inpF; } struct TreeFile { TFile *file; TTree *tree; TreeFile() : file(nullptr), tree(nullptr) {} TreeFile(TreeFile const &other) = delete; TreeFile(TreeFile &&other) : file(other.file), tree(other.tree) { other.file = nullptr; other.tree = nullptr; } ~TreeFile() { if (file) { file->Close(); } } }; inline std::unique_ptr CheckGetTTree(std::string const &fname, std::string const &tname, char const *opts = "") { TreeFile tf; tf.file = CheckOpenTFile(fname, opts); tf.tree = dynamic_cast(tf.file->Get(tname.c_str())); if (!tf.tree) { throw failed_to_get_TTree() << "[ERROR]: Failed to get TTree named: " << std::quoted(tname) << " from TFile: " << std::quoted(fname); } return std::make_unique(std::move(tf)); } inline std::unique_ptr MakeNewTTree(std::string const &fname, std::string const &tname, char const *opts = "") { TreeFile tf; tf.file = CheckOpenTFile(fname, opts); tf.tree = new TTree(tname.c_str(), ""); tf.tree->SetDirectory(tf.file); return std::make_unique(std::move(tf)); } template struct TH_traits {}; template <> struct TH_traits { + static size_t const NDims = 1; static std::string name() { return "TH1"; } + static Int_t NbinsIncludeFlow(TH1 *const h) { + return h->GetXaxis()->GetNbins() + 2; + } }; template <> struct TH_traits { + static size_t const NDims = 1; static std::string name() { return "TH1D"; } + static Int_t NbinsIncludeFlow(TH1D *const h) { + return h->GetXaxis()->GetNbins() + 2; + } }; template <> struct TH_traits { + static size_t const NDims = 1; static std::string name() { return "TH1F"; } + static Int_t NbinsIncludeFlow(TH1F *const h) { + return h->GetXaxis()->GetNbins() + 2; + } }; template <> struct TH_traits { + static size_t const NDims = 2; static std::string name() { return "TH2"; } + static Int_t NbinsIncludeFlow(TH2 *const h) { + return (h->GetXaxis()->GetNbins() + 2) * (h->GetYaxis()->GetNbins() + 2); + } }; template <> struct TH_traits { + static size_t const NDims = 2; static std::string name() { return "TH2D"; } + static Int_t NbinsIncludeFlow(TH2D *const h) { + return (h->GetXaxis()->GetNbins() + 2) * (h->GetYaxis()->GetNbins() + 2); + } }; template <> struct TH_traits { + static size_t const NDims = 2; static std::string name() { return "TH2F"; } + static Int_t NbinsIncludeFlow(TH2F *const h) { + return (h->GetXaxis()->GetNbins() + 2) * (h->GetYaxis()->GetNbins() + 2); + } }; template inline std::unique_ptr GetHistogramFromROOTFile(std::string const &fname, std::string const &hname) { TFile *f = CheckOpenTFile(fname, "READ"); HT *h = dynamic_cast(f->Get(hname.c_str())); if (!h) { throw invalid_histogram_name() << "[ERROR]: Failed to get " << TH_traits::name() << " named " << std::quoted(hname) << " from input file " << std::quoted(fname); } std::unique_ptr clone(dynamic_cast(h->Clone())); clone->SetDirectory(nullptr); f->Close(); delete f; return clone; } +template void Clear(HT *h) { + for (Int_t bin_it = 0; bin_it < TH_traits::NbinsIncludeFlow(h); + ++bin_it) { + h->SetBinContent(bin_it, 0); + h->SetBinError(bin_it, 0); + } +} + +// Fill for 1D histograms uses type system to enforce correct number of +// dimensions +template +typename std::enable_if::NDims == 1, void>::type +Fill(HT *h, std::array::NDims> const &val, double weight = 1) { + h->Fill(val[0], weight); +} + +/// Fill for 2D histograms uses type system to enforce correct number of +/// dimensions +template +typename std::enable_if::NDims == 2, void>::type +Fill(HT *h, std::array::NDims> const &val, double weight = 1) { + h->Fill(val[0], val[1], weight); +} + template -inline std::unique_ptr CloneHistogram(std::unique_ptr const &source, - bool clear = false) { +inline std::unique_ptr Clone(std::unique_ptr const &source, + bool clear = false) { std::unique_ptr target(dynamic_cast(source->Clone())); if (!target) { throw failed_to_clone() << "[ERROR]: Failed to clone a " << TH_traits::name() << ", source = " << source.get(); } target->SetDirectory(nullptr); if (clear) { - target->Clear(); + Clear(target.get()); } return target; } } // namespace utility } // namespace nuis #endif