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