diff --git a/CMakeLists.txt b/CMakeLists.txt
index eb7a6d1..2b088f9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,78 +1,78 @@
# Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
cmake_minimum_required (VERSION 2.8 FATAL_ERROR)
#Use the compilers found in the path
find_program(CMAKE_C_COMPILER NAMES $ENV{CC} gcc PATHS ENV PATH NO_DEFAULT_PATH)
find_program(CMAKE_CXX_COMPILER NAMES $ENV{CXX} g++ PATHS ENV PATH NO_DEFAULT_PATH)
project(NUISANCE)
include(ExternalProject)
set (NUISANCE_VERSION_MAJOR 3)
set (NUISANCE_VERSION_MINOR 0)
set (NUISANCE_VERSION_REVISION 0)
set (NUISANCE_VERSION_STRING "v${NUISANCE_VERSION_MAJOR}r${NUISANCE_VERSION_MINOR}")
if(${NUISANCE_VERSION_REVISION} STRGREATER "0")
set (NUISANCE_VERSION_STRING "${NUISANCE_VERSION_STRING}p${NUISANCE_VERSION_REVISION}")
endif()
#Set this to TRUE to enable build debugging messages
set(BUILD_DEBUG_MSGS TRUE)
include(${CMAKE_SOURCE_DIR}/cmake/cmessage.cmake)
include(${CMAKE_SOURCE_DIR}/cmake/cacheVariables.cmake)
include(${CMAKE_SOURCE_DIR}/cmake/parseConfigApp.cmake)
include(${CMAKE_SOURCE_DIR}/cmake/StringAndListUtils.cmake)
cmessage(STATUS "CMAKE_INSTALL_PREFIX: \"${CMAKE_INSTALL_PREFIX}\"")
cmessage(STATUS "CMAKE_BUILD_TYPE: \"${CMAKE_BUILD_TYPE}\"")
################################################################################
# Check Dependencies
################################################################################
################################## ROOT ######################################
include(${CMAKE_SOURCE_DIR}/cmake/ROOTSetup.cmake)
################################# InputHandler #################################
include(${CMAKE_SOURCE_DIR}/cmake/InputHandlerSetup.cmake)
################################# Pythia6/8 ###################################
include(${CMAKE_SOURCE_DIR}/cmake/pythia6Setup.cmake)
+#Want this before fhiclcpp which will add the install directory
+include_directories(${CMAKE_SOURCE_DIR}/src)
+
################################## FHICLCPP ####################################
include(${CMAKE_SOURCE_DIR}/cmake/fhiclcppSetup.cmake)
-#Want this after fhiclcpp which will add the install directory
-include_directories(${CMAKE_SOURCE_DIR}/src)
-
################################## COMPILER ####################################
include(${CMAKE_SOURCE_DIR}/cmake/c++CompilerSetup.cmake)
################################################################################
add_subdirectory(src)
add_subdirectory(scripts)
add_subdirectory(config)
add_subdirectory(data)
diff --git a/data/nuA/Nuclear/CMakeLists.txt b/data/nuA/Nuclear/CMakeLists.txt
index b7a435f..ca939a6 100644
--- a/data/nuA/Nuclear/CMakeLists.txt
+++ b/data/nuA/Nuclear/CMakeLists.txt
@@ -1,3 +1,4 @@
SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/Nuclear)
add_subdirectory(T2K)
+add_subdirectory(MINERvA)
diff --git a/data/nuA/Nuclear/MINERvA/CC0Pi/CMakeLists.txt b/data/nuA/Nuclear/MINERvA/CC0Pi/CMakeLists.txt
new file mode 100644
index 0000000..f3c91ea
--- /dev/null
+++ b/data/nuA/Nuclear/MINERvA/CC0Pi/CMakeLists.txt
@@ -0,0 +1,3 @@
+SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/CC0Pi)
+
+add_subdirectory(NProt_CH_xsec_STV_nu)
diff --git a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt b/data/nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/CMakeLists.txt
similarity index 59%
copy from data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt
copy to data/nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/CMakeLists.txt
index 18cc6d1..2cb89d2 100644
--- a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt
+++ b/data/nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/CMakeLists.txt
@@ -1,5 +1,5 @@
-SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/CC0Pi)
+SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/NProt_CH_xsec_STV_nu)
FILE(GLOB DATA_FILES *.root)
install(FILES ${DATA_FILES} DESTINATION ${DATA_INSTALL_DIR})
diff --git a/data/nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/MINERvA_1805.05486.root b/data/nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/MINERvA_1805.05486.root
new file mode 100644
index 0000000..7d51761
Binary files /dev/null and b/data/nuA/Nuclear/MINERvA/CC0Pi/NProt_CH_xsec_STV_nu/MINERvA_1805.05486.root differ
diff --git a/data/nuA/Nuclear/MINERvA/CMakeLists.txt b/data/nuA/Nuclear/MINERvA/CMakeLists.txt
new file mode 100644
index 0000000..4a70bbe
--- /dev/null
+++ b/data/nuA/Nuclear/MINERvA/CMakeLists.txt
@@ -0,0 +1,3 @@
+SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/MINERvA)
+
+add_subdirectory(CC0Pi)
diff --git a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt b/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt
index 18cc6d1..e447144 100644
--- a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt
+++ b/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt
@@ -1,5 +1,4 @@
SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/CC0Pi)
-FILE(GLOB DATA_FILES *.root)
-
-install(FILES ${DATA_FILES} DESTINATION ${DATA_INSTALL_DIR})
+add_subdirectory(NProt_CH_xsec_STV_nu)
+add_subdirectory(H20_xsec_2Dpcthetamu_numubar)
diff --git a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt b/data/nuA/Nuclear/T2K/CC0Pi/H20_xsec_2Dpcthetamu_numubar/CMakeLists.txt
similarity index 56%
copy from data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt
copy to data/nuA/Nuclear/T2K/CC0Pi/H20_xsec_2Dpcthetamu_numubar/CMakeLists.txt
index 18cc6d1..5ca9b21 100644
--- a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt
+++ b/data/nuA/Nuclear/T2K/CC0Pi/H20_xsec_2Dpcthetamu_numubar/CMakeLists.txt
@@ -1,5 +1,5 @@
-SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/CC0Pi)
+SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/H20_xsec_2Dpcthetamu_numubar)
FILE(GLOB DATA_FILES *.root)
install(FILES ${DATA_FILES} DESTINATION ${DATA_INSTALL_DIR})
diff --git a/data/nuA/Nuclear/T2K/CC0Pi/H2O_xsec_2Dpmuthetamu_numubar.root b/data/nuA/Nuclear/T2K/CC0Pi/H20_xsec_2Dpcthetamu_numubar/H2O_xsec_2Dpmuthetamu_numubar.root
similarity index 100%
rename from data/nuA/Nuclear/T2K/CC0Pi/H2O_xsec_2Dpmuthetamu_numubar.root
rename to data/nuA/Nuclear/T2K/CC0Pi/H20_xsec_2Dpcthetamu_numubar/H2O_xsec_2Dpmuthetamu_numubar.root
diff --git a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt b/data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/CMakeLists.txt
similarity index 59%
copy from data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt
copy to data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/CMakeLists.txt
index 18cc6d1..2cb89d2 100644
--- a/data/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt
+++ b/data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/CMakeLists.txt
@@ -1,5 +1,5 @@
-SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/CC0Pi)
+SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/NProt_CH_xsec_STV_nu)
FILE(GLOB DATA_FILES *.root)
install(FILES ${DATA_FILES} DESTINATION ${DATA_INSTALL_DIR})
diff --git a/data/nuA/Nuclear/T2K/CC0Pi/datResults.root b/data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/datResults.root
similarity index 100%
rename from data/nuA/Nuclear/T2K/CC0Pi/datResults.root
rename to data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/datResults.root
diff --git a/data/nuA/Nuclear/T2K/CC0Pi/dphitResults.root b/data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/dphitResults.root
similarity index 100%
rename from data/nuA/Nuclear/T2K/CC0Pi/dphitResults.root
rename to data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/dphitResults.root
diff --git a/data/nuA/Nuclear/T2K/CC0Pi/dptResults.root b/data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/dptResults.root
similarity index 100%
rename from data/nuA/Nuclear/T2K/CC0Pi/dptResults.root
rename to data/nuA/Nuclear/T2K/CC0Pi/NProt_CH_xsec_STV_nu/dptResults.root
diff --git a/data/nuA/Nuclear/T2K/CMakeLists.txt b/data/nuA/Nuclear/T2K/CMakeLists.txt
index 25096b5..da3e52c 100644
--- a/data/nuA/Nuclear/T2K/CMakeLists.txt
+++ b/data/nuA/Nuclear/T2K/CMakeLists.txt
@@ -1,3 +1,4 @@
SET(DATA_INSTALL_DIR ${DATA_INSTALL_DIR}/T2K)
add_subdirectory(CC0Pi)
+add_subdirectory(CC1Pi)
diff --git a/src/app/CMakeLists.txt b/src/app/CMakeLists.txt
index bb1876d..c84c5fe 100644
--- a/src/app/CMakeLists.txt
+++ b/src/app/CMakeLists.txt
@@ -1,15 +1,15 @@
SET(APPS nuissamples nuiscomp nuisevsum nuisstudy)
foreach(a ${APPS})
add_executable(${a} ${a}.cxx)
- target_link_libraries(${a} -Wl,--no-as-needed nuis_event nuis_input nuis_utility nuis_config nuis_persistency nuis_plugins nuis_params nuis_variation)
+ target_link_libraries(${a} -Wl,--no-as-needed nuis_event nuis_input nuis_utility_experimental nuis_utility nuis_config nuis_persistency nuis_plugins nuis_params nuis_variation)
target_link_libraries(${a} ${GENERATOR_DEPENDENT_TARGETS} -Wl,--as-needed)
target_link_libraries(${a} ${NUISANCE_LINK_DIRS})
target_link_libraries(${a} ${NUISANCE_DEPEND_LIBS})
if(NOT "${CMAKE_LINK_FLAGS}" STREQUAL "")
set_target_properties(${a} PROPERTIES LINK_FLAGS ${CMAKE_LINK_FLAGS})
endif()
install(TARGETS ${a} DESTINATION bin)
endforeach()
diff --git a/src/app/nuisstudy.cxx b/src/app/nuisstudy.cxx
index 4098cd8..5d1d199 100644
--- a/src/app/nuisstudy.cxx
+++ b/src/app/nuisstudy.cxx
@@ -1,97 +1,112 @@
#include "config/GlobalConfiguration.hxx"
#include "input/InputManager.hxx"
#include "event/MinimalEvent.hxx"
#include "samples/IEventProcessor.hxx"
#include "plugins/Instantiate.hxx"
#include "exception/exception.hxx"
#include "persistency/ROOTOutput.hxx"
+#include "utility/StringUtility.hxx"
+
#include "fhiclcpp/make_ParameterSet.h"
#include
NEW_NUIS_EXCEPT(invalid_cli_arguments);
void SayUsage(char const *argv[]) {
std::cout << "[USAGE]: " << argv[0]
<< "\n"
"\t-c : FHiCL file containing study "
"configuration. \n"
"\t-s : FHiCL key of a single sample "
"to run from the -c argument. \n"
+ "\t-o : Default output stream name "
+ "override, if unspecified mode will be 'CREATE'. \n"
<< std::endl;
}
std::string fhicl_file = "";
std::string named_sample = "";
+std::string default_stream_override = "";
void handleOpts(int argc, char const *argv[]) {
int opt = 1;
while (opt < argc) {
if ((std::string(argv[opt]) == "-?") ||
(std::string(argv[opt]) == "--help")) {
SayUsage(argv);
exit(0);
} else if (std::string(argv[opt]) == "-c") {
fhicl_file = argv[++opt];
} else if (std::string(argv[opt]) == "-s") {
named_sample = argv[++opt];
+ } else if (std::string(argv[opt]) == "-o") {
+ default_stream_override = argv[++opt];
} else {
std::cout << "[ERROR]: Unknown option: " << argv[opt] << std::endl;
SayUsage(argv);
exit(1);
}
opt++;
}
}
int main(int argc, char const *argv[]) {
nuis::config::EnsureConfigurationRead("nuis.global.config.fcl");
nuis::config::EnsureConfigurationRead("nuis.datacomparisons.fcl");
handleOpts(argc, argv);
if (!fhicl_file.size()) {
SayUsage(argv);
throw invalid_cli_arguments();
}
nuis::config::EnsureConfigurationRead(fhicl_file);
+ if (default_stream_override.size()) {
+ std::vector split =
+ nuis::utility::split(default_stream_override, ":");
+ std::string open_mode = (split.size() > 1) ? split[1] : "CREATE";
+
+ nuis::persistency::NewStream("default", split[0], open_mode);
+ }
+
size_t NMax = nuis::config::GetDocument().get(
"nmax", std::numeric_limits::max());
std::vector samples;
if (named_sample.size()) {
samples.push_back(
nuis::config::GetDocument().get(named_sample));
} else {
samples = nuis::config::GetDocument().get>(
"samples");
}
for (fhicl::ParameterSet const &samp_config : samples) {
std::cout << "[INFO]: Reading sample: "
<< samp_config.get("name") << std::endl;
nuis::plugins::plugin_traits::unique_ptr_t sample =
nuis::plugins::Instantiate(
samp_config.get("name"));
sample->Initialize(samp_config);
sample->ProcessSample(NMax);
sample->Write();
// Ensures no re-use of samples but cleans up the memory.
nuis::input::InputManager::Get().Clear();
}
nuis::persistency::CloseOpenTFiles();
}
diff --git a/src/event/Particle.hxx b/src/event/Particle.hxx
index f71d2e9..67a7863 100644
--- a/src/event/Particle.hxx
+++ b/src/event/Particle.hxx
@@ -1,89 +1,90 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#pragma once
#include "event/types.hxx"
#include "TLorentzVector.h"
namespace nuis {
namespace event {
class Particle {
public:
NEW_NUIS_EXCEPT(invalid_particle);
#define STATUS_LIST \
X(kNuclearLeaving, 0) \
X(kPrimaryInitialState, 1) \
X(kPrimaryFinalState, 2) \
X(kIntermediate, 3) \
X(kUnknown, 4) \
X(kBlocked, 5) \
X(kNParticleStatus, 6)
#define X(A, B) A = B,
enum class Status_t { STATUS_LIST };
#undef X
Particle();
Particle(Particle const &);
PDG_t pdg;
TLorentzVector P4;
bool operator!() const { return (pdg == std::numeric_limits::max()); }
double E() const { return P4.E(); }
double KE() const { return P4.E() - P4.M(); }
double P() const { return P4.Vect().Mag(); }
TVector3 P3() const { return P4.Vect(); }
TVector3 Dir() const { return P3().Unit(); }
double M() const { return P4.M(); }
double Theta() const { return P4.Vect().Theta(); }
+ double Theta_deg() const { return P4.Vect().Theta() * (180.0 / M_PI); }
double CosTheta() const { return P4.Vect().CosTheta(); }
};
} // namespace event
} // namespace nuis
#define X(A, B) \
case nuis::event::Particle::Status_t::A: { \
return os << #A; \
}
inline std::ostream &operator<<(std::ostream &os,
nuis::event::Particle::Status_t te) {
switch (te) { STATUS_LIST }
return os;
}
#undef X
#undef STATUS_LIST
inline bool operator==(nuis::event::Particle const &l,
nuis::event::Particle const &r) {
if (l.pdg != r.pdg) {
return false;
}
for (size_t i = 0; i < 4; ++i) {
if (l.P4[i] != r.P4[i]) {
return false;
}
}
return true;
}
diff --git a/src/event/types.hxx b/src/event/types.hxx
index 8850d33..cfcc983 100644
--- a/src/event/types.hxx
+++ b/src/event/types.hxx
@@ -1,129 +1,124 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef EVENT_TYPES_HXX_SEEN
#define EVENT_TYPES_HXX_SEEN
#include "exception/exception.hxx"
namespace nuis {
namespace event {
// Modelled off the NEUT interaction codes.
#define NUIS_INTERACTION_CHANNEL_LIST \
X(kCCQE, 1) \
X(kCC2p2h, 2) \
X(kCCSPP_PPip, 11) \
X(kCCSPP_PPi0, 12) \
X(kCCSPP_NPip, 13) \
X(kCCDiffPP, 15) \
X(kCCCohPi, 16) \
X(kCCResGamma, 17) \
X(kCCTransitionMPi, 21) \
X(kCCResEta0, 22) \
X(kCCResK, 23) \
X(kCCDIS, 26) \
\
X(kNCSPP_NPi0, 31) \
X(kNCSPP_PPi0, 32) \
X(kNCSPP_PPim, 33) \
X(kNCSPP_NPip, 34) \
X(kNCDiffPP, 35) \
X(kNCCohPi, 36) \
X(kNCResNGamma, 38) \
X(kNCResPGamma, 39) \
X(kNCTransitionMPi, 41) \
X(kNCResNEta0, 42) \
X(kNCResPEta0, 43) \
X(kNCResK0, 44) \
X(kNCResKp, 45) \
X(kNCDIS, 46) \
X(kNCELP, 51) \
X(kNCELN, 52) \
X(kNC2p2h, 53) \
\
X(kCCQE_nub, -1) \
X(kCC2p2h_nub, -2) \
X(kCCSPP_NPim_nub, -11) \
X(kCCSPP_NPi0_nub, -12) \
X(kCCSPP_PPim_nub, -13) \
X(kCCDiffPP_nub, -15) \
X(kCCCohPi_nub, -16) \
X(kCCResGamma_nub, -17) \
X(kCCTransitionMPi_nub, -21) \
X(kCCResEta0_nub, -22) \
X(kCCResK_nub, -23) \
X(kCCDIS_nub, -26) \
\
X(kNCSPP_NPi0_nub, -31) \
X(kNCSPP_PPi0_nub, -32) \
X(kNCSPP_PPim_nub, -33) \
X(kNCSPP_NPip_nub, -34) \
X(kNCDiffPP_nub, -35) \
X(kNCCohPi_nub, -36) \
X(kNCResNGamma_nub, -38) \
X(kNCResPGamma_nub, -39) \
X(kNCTransitionMPi_nub, -41) \
X(kNCResNEta0_nub, -42) \
X(kNCResPEta0_nub, -43) \
X(kNCResK0_nub, -44) \
X(kNCResKp_nub, -45) \
X(kNCDIS_nub, -46) \
X(kNCELP_nub, -51) \
X(kNCELN_nub, -52) \
X(kNC2p2h_nub, -53) \
\
X(kUndefined, 0)
#define X(A, B) A = B,
enum class Channel_t { NUIS_INTERACTION_CHANNEL_LIST };
#undef X
#define X(A, B) \
case B: { \
return Channel_t::A; \
}
inline Channel_t FromNEUTCode(int nc) {
switch (nc) {
NUIS_INTERACTION_CHANNEL_LIST
default: { return Channel_t::kUndefined; }
}
}
#undef X
using PDG_t = long;
-inline bool Is2p2h(Channel_t chan) {
- return ((chan == Channel_t::kCC2p2h) || (chan == Channel_t::kNC2p2h) ||
- (chan == Channel_t::kCC2p2h_nub) || (chan == Channel_t::kNC2p2h_nub));
-}
-
} // namespace event
} // namespace nuis
#define X(A, B) \
case nuis::event::Channel_t::A: { \
return os << #A; \
}
inline std::ostream &operator<<(std::ostream &os, nuis::event::Channel_t te) {
switch (te) { NUIS_INTERACTION_CHANNEL_LIST }
return os;
}
#undef X
#endif
diff --git a/src/exception/exception.hxx b/src/exception/exception.hxx
index 294ee46..ec72518 100644
--- a/src/exception/exception.hxx
+++ b/src/exception/exception.hxx
@@ -1,40 +1,37 @@
-#ifndef EXCEPTION_EXCEPTION_SEEN
-#define EXCEPTION_EXCEPTION_SEEN
+#pragma once
#include
#include
#include
namespace nuis {
struct nuis_except : public std::exception {
std::stringstream msgstrm;
std::string msg;
nuis_except() : msgstrm(), msg() {}
nuis_except(nuis_except const &other)
: msgstrm(), msg() {
msgstrm << other.msg;
msg = other.msg;
}
const char *what() const noexcept { return msg.c_str(); }
template nuis_except &operator<<(T const &obj) {
msgstrm << obj;
msg = msgstrm.str();
return (*this);
}
};
} // namespace nuis
#define NEW_NUIS_EXCEPT(EXCEPT_NAME) \
struct EXCEPT_NAME : public nuis::nuis_except { \
EXCEPT_NAME() : nuis::nuis_except() {} \
EXCEPT_NAME(EXCEPT_NAME const &other) : nuis::nuis_except(other) {} \
template EXCEPT_NAME &operator<<(T const &obj) { \
msgstrm << obj; \
msg = msgstrm.str(); \
return (*this); \
} \
}
-
-#endif
diff --git a/src/generator/input/GENIEInputHandler.cxx b/src/generator/input/GENIEInputHandler.cxx
index baa784a..cc5c146 100644
--- a/src/generator/input/GENIEInputHandler.cxx
+++ b/src/generator/input/GENIEInputHandler.cxx
@@ -1,149 +1,165 @@
#include "generator/input/GENIEInputHandler.hxx"
#include "generator/utility/GENIEUtility.hxx"
#include "utility/PDGCodeUtility.hxx"
#ifdef GENIE_V3_INTERFACE
#include "Framework/EventGen/EventRecord.h"
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/GHEP/GHepRecord.h"
#else
#include "EVGCore/EventRecord.h"
#include "GHEP/GHepParticle.h"
#include "GHEP/GHepRecord.h"
#endif
#include "fhiclcpp/ParameterSet.h"
using namespace nuis;
using namespace nuis::event;
using namespace nuis::utility;
using namespace nuis::genietools;
-GENIEInputHandler::GENIEInputHandler() : fInputTree(), fGenieNtpl(nullptr) {}
+GENIEInputHandler::GENIEInputHandler()
+ : fInputTree(), fGenieNtpl(nullptr), fFileWeight(1) {}
GENIEInputHandler::GENIEInputHandler(GENIEInputHandler &&other)
: fInputTree(std::move(other.fInputTree)),
- fReaderEvent(std::move(other.fReaderEvent)), fGenieNtpl(nullptr) {}
+ fReaderEvent(std::move(other.fReaderEvent)), fGenieNtpl(nullptr),
+ fFileWeight(1) {}
void GENIEInputHandler::Initialize(fhicl::ParameterSet const &ps) {
fInputTree = CheckGetTTree(ps.get("file"), "gtree");
fInputTree.tree->SetBranchAddress("gmcrec", &fGenieNtpl);
fKeepIntermediates = ps.get("keep_intermediates", false);
fKeepNuclearParticles = ps.get("keep_nuclear_particles", false);
+
+ fhicl::ParameterSet XSecInfo = ps.get("xsec_info", {});
+ if (XSecInfo.has_key("weight")) {
+ fFileWeight = XSecInfo.get("weight");
+ std::cout << "[INFO]: Average GENIE XSecWeight = " << fFileWeight << ""
+ << std::endl;
+ } else if (XSecInfo.has_key("flux") || XSecInfo.has_key("target") ||
+ XSecInfo.has_key("spline_file")) {
+ std::cout
+ << "[INFO]: Attempting to build GENIE file weight with input splines."
+ << std::endl;
+ fFileWeight = genietools::GetFileWeight(XSecInfo) / double(GetNEvents());
+ std::cout << "[INFO]: Average GENIE XSecWeight = " << fFileWeight << ""
+ << std::endl;
+ }
}
MinimalEvent const &GENIEInputHandler::GetMinimalEvent(ev_index_t idx) const {
if (idx >= GetNEvents()) {
throw IInputHandler::invalid_entry()
<< "[ERROR]: Attempted to get entry " << idx
<< " from an InputHandler with only " << GetNEvents();
}
fInputTree.tree->GetEntry(idx);
genie::GHepRecord *GHep = static_cast(fGenieNtpl->event);
if (!GHep) {
throw invalid_GENIE_event()
<< "[ERROR]: GENIE event " << idx << " failed to contain a GHepRecord";
}
fReaderEvent.mode = GetEventChannel(*GHep);
TObjArrayIter iter(GHep);
genie::GHepParticle *p;
while ((p = (dynamic_cast((iter).Next())))) {
if (!p) {
continue;
}
Particle::Status_t state = GetParticleStatus(*p, fReaderEvent.mode);
if (state != Particle::Status_t::kPrimaryInitialState) {
continue;
}
if (!IsNeutralLepton(p->Pdg(), pdgcodes::kMatterAntimatter) &&
!IsChargedLepton(p->Pdg(), pdgcodes::kMatterAntimatter)) {
continue;
}
fReaderEvent.probe_E = p->E() * 1.E3;
fReaderEvent.probe_pdg = p->Pdg();
break;
}
- fReaderEvent.XSecWeight = 1;
+ fReaderEvent.XSecWeight = fFileWeight;
if (fWeightCache.size() <= idx) {
fWeightCache.push_back(fReaderEvent.XSecWeight);
}
return fReaderEvent;
}
FullEvent const &GENIEInputHandler::GetFullEvent(ev_index_t idx) const {
(void)GetMinimalEvent(idx);
fReaderEvent.ClearParticleStack();
genie::GHepRecord *GHep = static_cast(fGenieNtpl->event);
unsigned int npart = GHep->GetEntries();
// Fill Particle Stack
genie::GHepParticle *p = 0;
TObjArrayIter iter(GHep);
// Loop over all particles
while ((p = (dynamic_cast((iter).Next())))) {
if (!p) {
continue;
}
// Get Status
Particle::Status_t state = GetParticleStatus(*p, fReaderEvent.mode);
- if (!fKeepIntermediates) {
+ if (!fKeepIntermediates && (state == Particle::Status_t::kIntermediate)) {
continue;
}
if (!fKeepNuclearParticles && IsNuclearPDG(p->Pdg())) {
continue;
}
Particle nuis_part;
-
nuis_part.pdg = p->Pdg();
nuis_part.P4 = TLorentzVector(p->Px(), p->Py(), p->Pz(), p->E()) * 1E3;
fReaderEvent.ParticleStack[static_cast(state)].particles.push_back(
nuis_part);
}
return fReaderEvent;
}
double GENIEInputHandler::GetEventWeight(ev_index_t idx) const {
if (idx > fWeightCache.size()) {
throw weight_cache_miss()
<< "[ERROR]: Failed to get cached weight for event index: " << idx;
}
return fWeightCache[idx];
}
double GENIEInputHandler::GetXSecScaleFactor(
std::pair const &EnuRange) const {
throw input_handler_feature_unimplemented()
<< "[ERROR]: Flux cuts not yet implemented for GENIE input handler.";
}
size_t GENIEInputHandler::GetNEvents() const {
return fInputTree.tree->GetEntries();
}
GeneratorManager::Generator_id_t GENIEInputHandler::GetGeneratorId() const {
return GeneratorManager::Get().EnsureGeneratorRegistered("GENIE");
}
DECLARE_PLUGIN(IInputHandler, GENIEInputHandler);
diff --git a/src/generator/input/GENIEInputHandler.hxx b/src/generator/input/GENIEInputHandler.hxx
index 7c87f17..312b6b7 100644
--- a/src/generator/input/GENIEInputHandler.hxx
+++ b/src/generator/input/GENIEInputHandler.hxx
@@ -1,67 +1,69 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#pragma once
#include "event/FullEvent.hxx"
#include "input/IInputHandler.hxx"
#include "exception/exception.hxx"
#include "utility/ROOTUtility.hxx"
#ifdef GENIE_V3_INTERFACE
#include "Framework/Ntuple/NtpMCEventRecord.h"
#else
#include "Ntuple/NtpMCEventRecord.h"
#endif
-#include
-
namespace fhicl {
-class ParameterSet;
+ class ParameterSet;
}
+#include
+
class GENIEInputHandler : public IInputHandler {
mutable nuis::utility::TreeFile fInputTree;
mutable nuis::event::FullEvent fReaderEvent;
mutable std::vector fWeightCache;
mutable genie::NtpMCEventRecord *fGenieNtpl;
bool fKeepIntermediates;
bool fKeepNuclearParticles;
+ double fFileWeight;
+
public:
NEW_NUIS_EXCEPT(weight_cache_miss);
GENIEInputHandler();
GENIEInputHandler(GENIEInputHandler const &) = delete;
GENIEInputHandler(GENIEInputHandler &&);
void Initialize(fhicl::ParameterSet const &);
nuis::event::MinimalEvent const &GetMinimalEvent(ev_index_t idx) const;
nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const;
double GetEventWeight(ev_index_t idx) const;
size_t GetNEvents() const;
double GetXSecScaleFactor(std::pair const &EnuRange) const;
nuis::GeneratorManager::Generator_id_t GetGeneratorId() const;
};
diff --git a/src/generator/input/NEUTInputHandler.cxx b/src/generator/input/NEUTInputHandler.cxx
index ea62b19..9fc688d 100644
--- a/src/generator/input/NEUTInputHandler.cxx
+++ b/src/generator/input/NEUTInputHandler.cxx
@@ -1,257 +1,260 @@
#include "generator/input/NEUTInputHandler.hxx"
#include "persistency/ROOTOutput.hxx"
#include "utility/HistogramUtility.hxx"
#include "utility/InteractionChannelUtility.hxx"
#include "utility/PDGCodeUtility.hxx"
#include "generator/utility/NEUTUtility.hxx"
#include "neutpart.h"
#include "neutvect.h"
#include "fhiclcpp/ParameterSet.h"
using namespace nuis;
using namespace nuis::event;
using namespace nuis::utility;
using namespace nuis::neuttools;
NEUTInputHandler::NEUTInputHandler() : fInputTreeFile() {}
NEUTInputHandler::NEUTInputHandler(NEUTInputHandler &&other)
: fInputTreeFile(std::move(other.fInputTreeFile)),
fReaderEvent(std::move(other.fReaderEvent)) {}
void NEUTInputHandler::Initialize(fhicl::ParameterSet const &ps) {
fInputTreeFile = CheckGetTTree(ps.get("file"), "neuttree");
fNeutVect = nullptr;
fInputTreeFile.tree->SetBranchAddress("vectorbranch", &fNeutVect);
fKeepIntermediates = ps.get("keep_intermediates", false);
std::string flux_name = ps.get("flux_hist_name", "flux_numu");
std::string evtrt_name = ps.get("evtrt_hist_name", "evtrt_numu");
fXSecRescaleFactor = ps.get("xsec_rescale_factor",1);
std::string override_flux_file =
ps.get("override_flux_file", "");
if (override_flux_file.size()) {
fFlux = GetHistogramFromROOTFile(override_flux_file, flux_name);
RebuildEventRate(!ps.get("flux_hist_in_MeV", false));
} else {
fFlux =
GetHistogramFromROOTFile(fInputTreeFile.file, flux_name, false);
}
if (fFlux) {
if (!fEvtrt) {
fEvtrt = GetHistogramFromROOTFile(fInputTreeFile.file, evtrt_name);
}
fFileWeight = fEvtrt->Integral() * 1.0E-38 /
(fFlux->Integral() * double(GetNEvents()));
std::cout << "[INFO]: Average NEUT XSecWeight = " << fFileWeight << ""
<< std::endl;
} else {
std::cout
<< "[INFO]: Input NEUT file doesn't have flux information, assuming "
"mono energetic and calculating average cross-section..."
<< std::endl;
fFileWeight = GetMonoEXSecWeight() / double(GetNEvents());
std::cout << "[INFO]: Done (Average NEUT XSecWeight = " << fFileWeight
<< ")!" << std::endl;
}
fFileWeight *= fXSecRescaleFactor;
}
NEW_NUIS_EXCEPT(non_mono_energetic_input_file);
double NEUTInputHandler::GetMonoEXSecWeight() {
double xsec = 0;
double count = 0;
double E_first = std::numeric_limits::max();
size_t NEvents = GetNEvents();
size_t ShoutEvery = NEvents / 100;
std::cout << "[INFO]: Read " << 0 << "/" << NEvents << " NEUT events."
<< std::flush;
for (size_t nev_it = 0; nev_it < NEvents; ++nev_it) {
if (ShoutEvery && !(nev_it % ShoutEvery)) {
std::cout << "\r[INFO]: Read " << nev_it << "/" << NEvents
<< " NEUT events." << std::flush;
}
fInputTreeFile.tree->GetEntry(nev_it);
xsec += fNeutVect->Totcrs;
+ //Assume monoE.
+ return xsec * 1E-38;
+
count++;
if (E_first == std::numeric_limits::max()) {
NeutPart *part = fNeutVect->PartInfo(0);
E_first = part->fP.E();
} else {
NeutPart *part = fNeutVect->PartInfo(0);
if (E_first != part->fP.E()) {
throw non_mono_energetic_input_file()
<< "[ERROR]: When calculating NEUT xsec weight for a mono "
"energetic beam, found different energies: "
<< E_first << " and " << part->fP.E()
<< ". Cannot continue with this input file.";
}
}
}
std::cout << "\r[INFO]: Read " << NEvents << "/" << NEvents << " NEUT events."
<< std::endl;
return (xsec / count) * 1E-38;
}
void NEUTInputHandler::RebuildEventRate(bool FluxInGeV) {
auto XSec = Clone(fFlux, true, "xsec");
auto Entry = Clone(fFlux, true, "entry");
std::cout << "[INFO]: Rebuilding total cross-section prediction..."
<< std::endl;
size_t NEvents = GetNEvents();
size_t ShoutEvery = NEvents / 100;
std::cout << "[INFO]: Read " << 0 << "/" << NEvents << " NEUT events."
<< std::flush;
for (size_t nev_it = 0; nev_it < NEvents; ++nev_it) {
if (ShoutEvery && !(nev_it % ShoutEvery)) {
std::cout << "\r[INFO]: Read " << nev_it << "/" << NEvents
<< " NEUT events." << std::flush;
}
fInputTreeFile.tree->GetEntry(nev_it);
NeutPart *part = fNeutVect->PartInfo(0);
double E = part->fP.E() * (FluxInGeV ? 1E-3 : 1);
XSec->Fill(E, fNeutVect->Totcrs);
Entry->Fill(E);
}
std::cout << std::endl << "[INFO]: Done!" << std::endl;
// Binned total xsec
XSec->Divide(Entry.get());
fEvtrt = Clone(XSec, false, "evtrt");
// Event rate prediction
fEvtrt->Multiply(fFlux.get());
}
double
NEUTInputHandler::GetXSecScaleFactor(std::pair const &) const {
throw input_handler_feature_unimplemented()
<< "[ERROR]: Flux cuts not yet implemented for NEUT input handler.";
}
NeutVect *NEUTInputHandler::GetNeutEvent(ev_index_t idx) const {
if (idx >= GetNEvents()) {
throw IInputHandler::invalid_entry()
<< "[ERROR]: Attempted to get entry " << idx
<< " from an InputHandler with only " << GetNEvents();
}
fInputTreeFile.tree->GetEntry(idx);
return fNeutVect;
}
MinimalEvent const &NEUTInputHandler::GetMinimalEvent(ev_index_t idx) const {
if (idx >= GetNEvents()) {
throw IInputHandler::invalid_entry()
<< "[ERROR]: Attempted to get entry " << idx
<< " from an InputHandler with only " << GetNEvents();
}
fInputTreeFile.tree->GetEntry(idx);
fReaderEvent.mode = IntToChannel(fNeutVect->Mode);
size_t NPart = fNeutVect->Npart();
for (size_t part_it = 0; part_it < NPart; part_it++) {
NeutPart const &part = (*fNeutVect->PartInfo(part_it));
if ((part.fIsAlive == false) && (part.fStatus == -1) &&
IsNeutralLepton(part.fPID)) {
fReaderEvent.probe_E = part.fP.T();
fReaderEvent.probe_pdg = part.fPID;
break;
}
}
fReaderEvent.XSecWeight = fFileWeight;
if (fWeightCache.size() <= idx) {
fWeightCache.push_back(fReaderEvent.XSecWeight);
}
fReaderEvent.fGenEvent = fNeutVect;
return fReaderEvent;
}
FullEvent const &NEUTInputHandler::GetFullEvent(ev_index_t idx) const {
(void)GetMinimalEvent(idx);
fReaderEvent.ClearParticleStack();
if (fNeutVect->Ibound) {
fReaderEvent.target_pdg =
MakeNuclearPDG(fNeutVect->TargetA, fNeutVect->TargetZ);
} else {
fReaderEvent.target_pdg = MakeNuclearPDG(1, 1);
}
size_t NPart = fNeutVect->Npart();
size_t NPrimary = fNeutVect->Nprimary();
for (size_t part_it = 0; part_it < NPart; part_it++) {
NeutPart const &part = (*fNeutVect->PartInfo(part_it));
Particle nuis_part;
nuis_part.pdg = part.fPID;
nuis_part.P4 = part.fP;
Particle::Status_t state = GetNeutParticleStatus(part, fReaderEvent.mode);
if (!fKeepIntermediates && (state == Particle::Status_t::kIntermediate)) {
continue;
}
size_t state_int = static_cast(state);
// Add status == 0 particles as pre-FSI particles until we find an
// intermediate state particle
if ((part_it < NPrimary) &&
(state != Particle::Status_t::kPrimaryInitialState)) {
fReaderEvent
.ParticleStack[static_cast(
Particle::Status_t::kPrimaryFinalState)]
.particles.push_back(nuis_part);
}
fReaderEvent.ParticleStack[state_int].particles.push_back(nuis_part);
}
return fReaderEvent;
}
double NEUTInputHandler::GetEventWeight(ev_index_t idx) const {
if (idx > fWeightCache.size()) {
throw weight_cache_miss()
<< "[ERROR]: Failed to get cached weight for event index: " << idx;
}
return fWeightCache[idx];
}
size_t NEUTInputHandler::GetNEvents() const {
return fInputTreeFile.tree->GetEntries();
}
GeneratorManager::Generator_id_t NEUTInputHandler::GetGeneratorId() const {
return GeneratorManager::Get().EnsureGeneratorRegistered("NEUT");
}
DECLARE_PLUGIN(IInputHandler, NEUTInputHandler);
diff --git a/src/generator/utility/CMakeLists.txt b/src/generator/utility/CMakeLists.txt
index 00e5686..bca56b8 100644
--- a/src/generator/utility/CMakeLists.txt
+++ b/src/generator/utility/CMakeLists.txt
@@ -1,57 +1,63 @@
if(USE_NUWRO)
LIST(APPEND GENERATOR_UTILS_IMPL NuWroUtility.cxx)
LIST(APPEND GENERATOR_UTILS_HDR NuWroUtility.hxx)
endif(USE_NUWRO)
if(USE_NEUT)
LIST(APPEND GENERATOR_UTILS_IMPL NEUTUtility.cxx)
LIST(APPEND GENERATOR_UTILS_HDR NEUTUtility.hxx)
endif(USE_NEUT)
if(USE_GENIE)
- LIST(APPEND GENERATOR_UTILS_IMPL GENIEUtility.cxx)
- LIST(APPEND GENERATOR_UTILS_HDR GENIEUtility.hxx)
+ LIST(APPEND GENERATOR_UTILS_IMPL
+ GENIEUtility.cxx
+ GENIESplineReader.cxx
+ )
+ LIST(APPEND GENERATOR_UTILS_HDR
+ GENIEUtility.hxx
+ GENIESplineReader.hxx
+ )
endif(USE_GENIE)
if(GENERATOR_UTILS_IMPL)
add_library(nuis_generator_utility SHARED ${GENERATOR_UTILS_IMPL})
if(USE_GENIE)
target_compile_options(nuis_generator_utility PUBLIC ${GENIE_CXX_FLAGS})
include_directories(${GENIE_INCLUDE_DIRS})
target_link_libraries(nuis_generator_utility ${GENIE_LINK_DIRS})
target_link_libraries(nuis_generator_utility -Wl,--no-as-needed ${GENIE_LIBS})
endif()
if(USE_NEUT)
target_compile_options(nuis_generator_utility PUBLIC ${NEUT_CXX_FLAGS})
include_directories(${NEUT_INCLUDE_DIRS})
target_link_libraries(nuis_generator_utility ${NEUT_LINK_DIRS})
target_link_libraries(nuis_generator_utility -Wl,--no-as-needed ${NEUT_LIBS})
target_link_libraries(nuis_generator_utility ${NEUT_IMPORTED_TARGETS})
endif()
if(USE_NUWRO)
target_compile_options(nuis_generator_utility PUBLIC ${NUWRO_CXX_FLAGS})
include_directories(${NUWRO_INCLUDE_DIRS})
target_link_libraries(nuis_generator_utility -Wl,--no-as-needed ${NUWRO_IMPORTED_TARGETS})
endif()
if(${NEED_PYTHIA6})
target_link_libraries(nuis_generator_utility -L${PYTHIA6})
endif()
target_link_libraries(nuis_generator_utility -L${ROOT_LIBDIR} ${ROOT_LIBS})
install(TARGETS nuis_generator_utility DESTINATION lib)
endif()
if(GENERATOR_UTILS_HDR)
install(FILES ${GENERATOR_UTILS_HDR} DESTINATION include/generator/utility)
endif()
LIST(APPEND GENERATOR_DEPENDENT_TARGETS nuis_generator_utility)
SET(GENERATOR_DEPENDENT_TARGETS ${GENERATOR_DEPENDENT_TARGETS} PARENT_SCOPE)
diff --git a/src/generator/utility/GENIESplineReader.cxx b/src/generator/utility/GENIESplineReader.cxx
new file mode 100644
index 0000000..9f9657c
--- /dev/null
+++ b/src/generator/utility/GENIESplineReader.cxx
@@ -0,0 +1,483 @@
+#include "GENIESplineReader.hxx"
+
+#include "exception/exception.hxx"
+
+#include "utility/StringUtility.hxx"
+
+#include "string_parsers/from_string.hxx"
+#include "string_parsers/utility.hxx"
+
+#include
+
+#include
+#include
+
+// #define SPLINE_NAME_MATCHER_VERBOSE
+
+namespace nuis {
+namespace genietools {
+
+SSAX::SSAX()
+ : EventOnAllTags(true), LeavingAlready(false), InterestingTags(),
+ buf_size(1), buf_usage(0), buf(NULL), depth(0) {
+ buf = new char[buf_size];
+};
+
+SSAX::~SSAX() { delete buf; }
+
+void SSAX::OnStartDocument() {}
+void SSAX::OnEndDocument() {}
+void SSAX::OnStartElement(
+ const char *name,
+ std::vector> &attributes) {
+ std::cout << "[SSAX]: OnStartElement(\"" << name << "\", " << std::endl;
+ for (size_t a_it = 0; a_it < attributes.size(); ++a_it) {
+ std::cout << "\t\"" << attributes[a_it].first << "\" = \""
+ << attributes[a_it].second << "\"" << std::endl;
+ }
+ std::cout << "\t)" << std::endl;
+}
+void SSAX::OnCharacters(const char *characters) {
+ std::cout << "[SSAX]: OnCharacters(\"" << characters << "\")" << std::endl;
+}
+void SSAX::OnEndElement(const char *name) {
+ std::cout << "[SSAX]: OnEndElement(\"" << name << "\")" << std::endl;
+}
+void SSAX::OnComment(const char *) {}
+
+void SSAX::NotifyInterestingTags(std::vector it) {
+ InterestingTags = it;
+ EventOnAllTags = !InterestingTags.size();
+}
+void SSAX::NotifyEarlyExit() { LeavingAlready = true; }
+
+void SSAX::ExpandBuf() {
+ size_t new_buf_size = (buf_size <= 1E6) ? buf_size * 10 : (buf_size + 1E6);
+
+ std::cout << "[SSAX]: Expanding internal buffer to: " << new_buf_size
+ << " bytes." << std::endl;
+
+ char *new_buf = new char[new_buf_size];
+ memcpy(new_buf, buf, buf_size);
+ delete[] buf;
+
+ buf = new_buf;
+ buf_size = new_buf_size;
+}
+
+void SSAX::AddToBuf(char character) {
+ if (buf_usage == buf_size) {
+ ExpandBuf();
+ }
+ buf[buf_usage++] = character;
+}
+
+void SSAX::NullTermBuf() { AddToBuf('\0'); }
+
+bool SSAX::IsNameStartChar(char c) {
+ return (isalpha(c) || (c == ':') || (c == '_'));
+}
+
+bool SSAX::IsNameChar(char c) {
+ return (IsNameStartChar(c) || isdigit(c) || (c == '-') || (c == '.'));
+}
+
+void SSAX::IgnoreAllToRightBracket(std::ifstream &ifs) {
+ char prev, next;
+ while (EOF != (next = ifs.get())) {
+ if (next == '>') {
+ if (prev == '/') { // Empty element
+ depth--;
+ }
+
+ return;
+ }
+ prev = next;
+ }
+ std::cerr << "[SSAX - ERROR]: encountered EOF before ending '>' at depth: "
+ << depth << std::endl;
+ throw;
+}
+
+char SSAX::IgnoreSpace(std::ifstream &ifs) {
+ char next;
+ while (EOF != (next = ifs.get())) {
+ if (isspace(next)) {
+ continue;
+ } else {
+ return next;
+ }
+ }
+ std::cerr << "[SSAX - ERROR]: encountered EOF early at depth: " << depth
+ << std::endl;
+ throw;
+}
+
+char SSAX::BufNameToBoundary(std::ifstream &ifs) {
+ char next;
+ while (EOF != (next = ifs.get())) {
+ if (isspace(next) || (next == '>') || (next == '/') || (next == '=')) {
+ return next;
+ }
+ if (IsNameChar(next)) {
+ AddToBuf(next);
+ } else {
+ std::cerr << "[SSAX - ERROR]: encountered bad character " << next
+ << std::endl;
+ throw;
+ }
+ }
+ std::cerr << "[SSAX - ERROR]: encountered EOF early at depth: " << depth
+ << std::endl;
+ throw;
+}
+
+char SSAX::BufAttrValueToMatchingQuote(std::ifstream &ifs, char quot) {
+ char next;
+ while (EOF != (next = ifs.get())) {
+ if (next == quot) {
+ return next;
+ }
+ AddToBuf(next);
+ }
+ std::cerr << "[SSAX - ERROR]: encountered EOF early at depth: " << depth
+ << std::endl;
+ throw;
+}
+
+std::pair
+SSAX::HandleAttribute(std::ifstream &ifs) {
+ return std::pair(nullptr, nullptr);
+}
+
+void SSAX::HandleOpenTag(std::ifstream &ifs) {
+ char next = ifs.get();
+ size_t AttrNameStart = 0;
+ size_t AttrValueStart = 0;
+ attrs.clear();
+ bool IsInteresting = true;
+
+ if (next == EOF) {
+ std::cerr << "[SSAX - ERROR]: encountered EOF directly after element "
+ "opening at depth: "
+ << depth << std::endl;
+ throw;
+ }
+
+ if (!IsNameStartChar(next)) {
+ std::cerr << "[SSAX - ERROR]: Bad element name first character: " << next
+ << std::endl;
+ throw;
+ }
+ AddToBuf(next);
+
+ next = BufNameToBoundary(ifs);
+ NullTermBuf();
+ AttrNameStart = buf_usage;
+
+ if (!EventOnAllTags) {
+ IsInteresting = false;
+ for (size_t it_it = 0; it_it < InterestingTags.size(); ++it_it) {
+ if (InterestingTags[it_it] == buf) {
+ IsInteresting = true;
+ break;
+ }
+ }
+ }
+
+ if ((next == '>')) { // Handles elements like
+ if (IsInteresting) {
+ OnStartElement(buf, attrs);
+ }
+ depth++;
+ return;
+ }
+ if ((next == '/')) { // Handles elements like
+ next = ifs.get();
+ if (EOF == next) {
+ std::cerr << "[SSAX - ERROR]: encountered EOF early at depth: " << depth
+ << std::endl;
+ throw;
+ }
+
+ if ((next == '>')) {
+ if (IsInteresting) {
+ OnStartElement(buf, attrs);
+ }
+ return;
+ } else {
+ std::cerr << "[SSAX - ERROR]: encountered unexpected character " << next
+ << " after /. Expected >. Exiting early at depth: " << depth
+ << std::endl;
+ throw;
+ }
+ }
+
+ while (EOF != (next = ifs.get())) {
+ if (isspace(next)) {
+ next = IgnoreSpace(ifs);
+ }
+
+ if (IsNameStartChar(next)) {
+ AddToBuf(next);
+ next = BufNameToBoundary(ifs);
+ if (next != '=') {
+ std::cerr << "[SSAX - ERROR]: encountered unexpected character \""
+ << next
+ << "\". Expected \"=\". Exiting early at depth: " << depth
+ << std::endl;
+ throw;
+ }
+ NullTermBuf();
+ AttrValueStart = buf_usage;
+
+ next = ifs.get();
+ if ((next != '\"') && (next != '\'')) {
+ std::cerr << "[SSAX - ERROR]: encountered unexpected character \""
+ << next << "\" (" << int(next) << "). Expected \"\'\"("
+ << int('\'') << ") or \"\"\"(" << int('\"')
+ << "). Exiting early at depth: " << depth << std::endl;
+ throw;
+ }
+ next = BufAttrValueToMatchingQuote(ifs, next);
+ NullTermBuf();
+ attrs.push_back(
+ std::make_pair(&buf[AttrNameStart], &buf[AttrValueStart]));
+ AttrNameStart = buf_usage;
+ continue;
+ }
+
+ if ((next == '>')) { // Handles elements like
+ if (IsInteresting) {
+ OnStartElement(buf, attrs);
+ }
+ depth++;
+ return;
+ }
+ if ((next == '/')) { // Handles elements like
+ next = ifs.get();
+ if (EOF == next) {
+ std::cerr << "[SSAX - ERROR]: encountered EOF early at depth: " << depth
+ << std::endl;
+ throw;
+ }
+
+ if ((next == '>')) {
+ if (IsInteresting) {
+ OnStartElement(buf, attrs);
+ }
+ // No depth ++;
+ return;
+ } else {
+ std::cerr << "[SSAX - ERROR]: encountered unexpected character " << next
+ << " after /. Expected >. Exiting early at depth: " << depth
+ << std::endl;
+ throw;
+ }
+ }
+
+ std::cerr << "[SSAX - ERROR]: encountered unexpected character " << next
+ << " after /. Expected >. Exiting early at depth: " << depth
+ << std::endl;
+ throw;
+ }
+ std::cerr << "[SSAX - ERROR]: encountered EOF early at depth: " << depth
+ << std::endl;
+ throw;
+}
+
+void SSAX::HandleCloseTag(std::ifstream &ifs) {
+ ifs.ignore();
+ char next = ifs.get();
+ bool IsInteresting = true;
+
+ if (!IsNameStartChar(next)) {
+ std::cerr << "[SSAX - ERROR]: Bad element name first character: " << next
+ << std::endl;
+ throw;
+ }
+ AddToBuf(next);
+
+ next = BufNameToBoundary(ifs);
+ NullTermBuf();
+
+ if (!EventOnAllTags) {
+ IsInteresting = false;
+ for (size_t it_it = 0; it_it < InterestingTags.size(); ++it_it) {
+ if (InterestingTags[it_it] == buf) {
+ IsInteresting = true;
+ break;
+ }
+ }
+ }
+
+ if (IsInteresting) {
+ OnEndElement(buf);
+ }
+}
+
+void SSAX::HandleCommentTag(std::ifstream &ifs) {
+ IgnoreAllToRightBracket(ifs);
+}
+
+void SSAX::FlushBufferToCharacterHandler() {
+ NullTermBuf();
+ OnCharacters(buf);
+}
+void SSAX::ClearBuf() { buf_usage = 0; }
+
+void SSAX::ParseFile(const char *filename) {
+ std::ifstream ifs(filename);
+
+ if (!ifs.good()) {
+ std::cerr << "[SSAX]: Failed to open " << filename << " for reading."
+ << std::endl;
+ throw;
+ }
+ ifs.clear();
+ ifs.seekg(0);
+
+ OnStartDocument();
+
+ char next;
+ while (EOF != (next = ifs.get())) {
+ if (LeavingAlready) {
+ return;
+ }
+ if (next == '<') {
+ FlushBufferToCharacterHandler();
+ ClearBuf();
+ switch (ifs.peek()) {
+ case '!':
+ case '?': {
+ HandleCommentTag(ifs);
+ ClearBuf();
+ break;
+ }
+ case '/': {
+ HandleCloseTag(ifs);
+ ClearBuf();
+ break;
+ }
+ default: {
+ HandleOpenTag(ifs);
+ ClearBuf();
+ }
+ }
+ continue;
+ }
+ AddToBuf(next);
+ }
+ OnEndDocument();
+}
+
+GENIESplineGetter::GENIESplineGetter(std::vector sps) : SSAX() {
+ SplinePatternStrings = sps;
+ for (auto const &sp : SplinePatternStrings) {
+ SplinePatterns.emplace_back(nuis::utility::DeGlobPattern(sp));
+ Splines.push_back({});
+ }
+
+ ReadingDocument = false;
+ NSplines = 0;
+
+ InSplineElement = false;
+ InEElement = false;
+ InxsecElement = false;
+}
+
+void GENIESplineGetter::OnStartDocument() {
+ ReadingDocument = true;
+ std::cout << "[GENIESplineGetter]: Begin document..." << std::endl;
+}
+void GENIESplineGetter::OnEndDocument() {
+ ReadingDocument = false;
+ std::cout << "[GENIESplineGetter]: Finished document, read " << NSplines
+ << " splines." << std::endl;
+}
+
+void GENIESplineGetter::OnStartElement(
+ const char *name,
+ std::vector> &attributes) {
+ std::string name_str = std::string(name);
+
+ if (InSplineElement) {
+ if (name_str == "E") {
+ InEElement = true;
+ }
+ if (name_str == "xsec") {
+ InxsecElement = true;
+ }
+ }
+
+ if (name_str == "spline") {
+ for (size_t attr_it = 0; attr_it < attributes.size(); ++attr_it) {
+ if (std::string(attributes[attr_it].first) == "name") {
+ std::string attr_val = std::string(attributes[attr_it].second);
+#ifdef SPLINE_NAME_MATCHER_VERBOSE
+ std::cout << "[DEBUG]: Found spline named: " << attr_val << std::endl;
+#endif
+ for (sp_it = 0; sp_it < SplinePatterns.size(); ++sp_it) {
+#ifdef SPLINE_NAME_MATCHER_VERBOSE
+ std::cout << "[DEBUG]: Checking against "
+ << SplinePatternStrings[sp_it] << std::endl;
+#endif
+ if (std::regex_match(attr_val, SplinePatterns[sp_it])) { // match
+#ifdef SPLINE_NAME_MATCHER_VERBOSE
+ std::cout << "[DEBUG]: Matched! " << std::endl;
+#endif
+ InSplineElement = true;
+ Splines[sp_it].push_back({{}, {}, attr_val});
+ break;
+ }
+ }
+ return;
+ }
+ }
+ }
+}
+
+void GENIESplineGetter::OnCharacters(const char *characters) {
+ if (InEElement) {
+ std::string ch(characters);
+ fhicl::string_parsers::trim(ch);
+ Splines[sp_it].back().EValues.push_back(
+ fhicl::string_parsers::str2T(ch));
+ InEElement = false;
+ }
+ if (InxsecElement) {
+ std::string ch(characters);
+ fhicl::string_parsers::trim(ch);
+ Splines[sp_it].back().XSecValues.push_back(
+ fhicl::string_parsers::str2T(ch));
+ InxsecElement = false;
+ }
+}
+
+void GENIESplineGetter::OnEndElement(const char *name) {
+ if (InSplineElement && (std::string(name) == "spline")) {
+ std::cout << "[GENIESplineGetter]: Finished reading spline element: "
+ << Splines[sp_it].back().SplineName << std::endl;
+ InSplineElement = false;
+ NSplines++;
+ }
+}
+
+void GENIESplineGetter::OnComment(const char *) {}
+
+std::vector> GENIESplineGetter::GetTGraphs() {
+ std::vector> graphs;
+ for (sp_it = 0; sp_it < SplinePatterns.size(); ++sp_it) {
+ graphs.push_back({});
+ for (size_t spl_it = 0; spl_it < Splines[sp_it].size(); ++spl_it) {
+ graphs.back().emplace_back(Splines[sp_it][spl_it].EValues.size(),
+ Splines[sp_it][spl_it].EValues.data(),
+ Splines[sp_it][spl_it].XSecValues.data());
+ graphs.back().back().SetName(Splines[sp_it][spl_it].SplineName.c_str());
+ }
+ }
+ return graphs;
+}
+
+} // namespace genietools
+} // namespace nuis
diff --git a/src/generator/utility/GENIESplineReader.hxx b/src/generator/utility/GENIESplineReader.hxx
new file mode 100644
index 0000000..a3d798e
--- /dev/null
+++ b/src/generator/utility/GENIESplineReader.hxx
@@ -0,0 +1,110 @@
+#pragma once
+
+#include "TGraph.h"
+
+#include
+#include
+#include
+#include
+#include
+
+namespace nuis {
+namespace genietools {
+
+class SSAX {
+ bool EventOnAllTags;
+ bool LeavingAlready;
+ std::vector InterestingTags;
+
+ size_t buf_size;
+ size_t buf_usage;
+ char *buf;
+ std::vector> attrs;
+
+ size_t depth;
+
+public:
+ SSAX();
+
+ virtual ~SSAX();
+
+ virtual void OnStartDocument();
+ virtual void OnEndDocument();
+ virtual void OnStartElement(
+ const char *name,
+ std::vector> &attributes);
+ virtual void OnCharacters(const char *characters);
+ virtual void OnEndElement(const char *name);
+ virtual void OnComment(const char *);
+
+ void NotifyInterestingTags(std::vector it);
+ void NotifyEarlyExit();
+
+private:
+ void ExpandBuf();
+
+ void AddToBuf(char character);
+ void NullTermBuf();
+
+ bool IsNameStartChar(char c);
+ bool IsNameChar(char c);
+
+ void IgnoreAllToRightBracket(std::ifstream &ifs);
+
+ char IgnoreSpace(std::ifstream &ifs);
+
+ char BufNameToBoundary(std::ifstream &ifs);
+
+ char BufAttrValueToMatchingQuote(std::ifstream &ifs, char quot);
+
+ std::pair HandleAttribute(std::ifstream &ifs);
+
+ void HandleOpenTag(std::ifstream &ifs);
+ void HandleCloseTag(std::ifstream &ifs);
+ void HandleCommentTag(std::ifstream &ifs);
+
+ void FlushBufferToCharacterHandler();
+
+ void ClearBuf();
+
+public:
+ void ParseFile(const char *filename);
+};
+
+struct spline {
+ std::vector EValues;
+ std::vector XSecValues;
+ std::string SplineName;
+};
+
+class GENIESplineGetter : public SSAX {
+public:
+ std::vector SplinePatternStrings;
+ std::vector SplinePatterns;
+ std::vector> Splines;
+
+ bool InEElement;
+ bool InxsecElement;
+ bool InSplineElement;
+
+ bool ReadingDocument;
+
+ size_t sp_it;
+ size_t NSplines;
+
+ GENIESplineGetter(std::vector sn);
+
+ void OnStartDocument();
+ void OnEndDocument();
+ void OnStartElement(
+ const char *name,
+ std::vector> &attributes);
+ void OnCharacters(const char *characters);
+ void OnEndElement(const char *name);
+ void OnComment(const char *);
+
+ std::vector> GetTGraphs();
+};
+
+} // namespace genietools
+} // namespace nuis
diff --git a/src/generator/utility/GENIEUtility.cxx b/src/generator/utility/GENIEUtility.cxx
index 2d93aea..b3be697 100644
--- a/src/generator/utility/GENIEUtility.cxx
+++ b/src/generator/utility/GENIEUtility.cxx
@@ -1,298 +1,465 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#include "generator/utility/GENIEUtility.hxx"
+#include "generator/utility/GENIESplineReader.hxx"
+#include "utility/HistogramUtility.hxx"
+#include "utility/StringUtility.hxx"
+#include "utility/InteractionChannelUtility.hxx"
#include "utility/PDGCodeUtility.hxx"
#ifdef GENIE_V3_INTERFACE
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/GHEP/GHepRecord.h"
#include "Framework/GHEP/GHepUtils.h"
#include "Framework/ParticleData/PDGCodes.h"
#else
#include "GHEP/GHepParticle.h"
#include "GHEP/GHepRecord.h"
#include "GHEP/GHepUtils.h"
#include "PDG/PDGCodes.h"
#endif
+#include "fhiclcpp/ParameterSet.h"
+
+#include "string_parsers/from_string.hxx"
+
#include "TGraph.h"
namespace nuis {
using namespace event;
namespace genietools {
static std::map SplineCache;
TGraph dum;
TGraph const &GetGENIESpline(std::string const &SplineFile,
std::string const &SplineIdentifier) {
if (SplineCache.find(SplineFile + SplineIdentifier) != SplineCache.end()) {
return SplineCache.find(SplineFile + SplineIdentifier)->second;
}
return dum;
}
struct NFSParticleCount {
size_t NProton;
size_t NNeutron;
size_t NPip;
size_t NPi0;
size_t NPim;
size_t NOther;
};
NFSParticleCount CountPreFSIParticles(genie::GHepRecord const &ev) {
// This code in this method is adapted from the GENIE source code found in
// GHep/GHepUtils.cxx This method therefore carries the GENIE copyright
// licence as copied below:
//
/// Copyright (c) 2003-2017, GENIE Neutrino MC Generator Collaboration
/// For the full text of the license visit http://copyright.genie-mc.org
/// or see $GENIE/LICENSE
//
genie::Target const &tgt = ev.Summary()->InitState().Tgt();
if (!tgt.HitNucIsSet()) {
throw invalid_GENIE_event()
<< "[ERROR]: Failed to get hit nucleon kinematics as it was not "
"included in this GHep event. This is a fatal error.";
}
genie::GHepParticle *FSLep = ev.FinalStatePrimaryLepton();
genie::GHepParticle *ISLep = ev.Probe();
if (!FSLep || !ISLep) {
throw invalid_GENIE_event()
<< "[ERROR]: Failed to find IS and FS lepton in event: "
<< ev.Summary()->AsString();
}
size_t NPi0 = 0, NPip = 0, NPim = 0, NProton = 0, NNeutron = 0, NOther = 0;
bool nuclear_target = tgt.IsNucleus();
TIter event_iter(&ev);
genie::GHepParticle *p = 0;
while ((p = dynamic_cast(event_iter.Next()))) {
genie::GHepStatus_t ghep_ist = (genie::GHepStatus_t)p->Status();
int ghep_pdgc = p->Pdg();
int ghep_fm = p->FirstMother();
int ghep_fmpdgc = (ghep_fm == -1) ? 0 : ev.Particle(ghep_fm)->Pdg();
// For nuclear targets use hadrons marked as 'hadron in the nucleus'
// which are the ones passed in the intranuclear rescattering
// For free nucleon targets use particles marked as 'final state'
// but make an exception for decayed pi0's,eta's (count them and not their
// daughters)
bool decayed =
(ghep_ist == genie::kIStDecayedState &&
(ghep_pdgc == genie::kPdgPi0 || ghep_pdgc == genie::kPdgEta));
bool parent_included =
(ghep_fmpdgc == genie::kPdgPi0 || ghep_fmpdgc == genie::kPdgEta);
bool count_it =
(nuclear_target && ghep_ist == genie::kIStHadronInTheNucleus) ||
(!nuclear_target && decayed) ||
(!nuclear_target && ghep_ist == genie::kIStStableFinalState &&
!parent_included);
if (!count_it) {
continue;
}
if (ghep_pdgc == genie::kPdgPiP) {
NPip++;
} else if (ghep_pdgc == genie::kPdgPiM) {
NPim++;
} else if (ghep_pdgc == genie::kPdgPi0) {
NPi0++;
} else if (ghep_pdgc == genie::kPdgProton) {
NProton++;
} else if (ghep_pdgc == genie::kPdgNeutron) {
NNeutron++;
} else if (!utility::IsNeutralLepton(
ghep_pdgc, utility::pdgcodes::kMatterAntimatter) &&
!utility::IsChargedLepton(
ghep_pdgc, utility::pdgcodes::kMatterAntimatter)) {
NOther++;
}
}
return NFSParticleCount{NProton, NNeutron, NPip, NPi0, NPim, NOther};
}
Channel_t GetEventChannel(genie::GHepRecord const &gev) {
// Electron Scattering
if (gev.Summary()->ProcInfo().IsEM()) {
if (gev.Summary()->InitState().ProbePdg() == utility::pdgcodes::kElectron) {
if (gev.Summary()->ProcInfo().IsQuasiElastic()) {
NFSParticleCount fsparts = CountPreFSIParticles(gev);
if (fsparts.NProton) {
return Channel_t::kNCELP;
} else {
return Channel_t::kNCELN;
}
} else if (gev.Summary()->ProcInfo().IsMEC()) {
return Channel_t::kNC2p2h;
} else if (gev.Summary()->ProcInfo().IsResonant()) {
NFSParticleCount fsparts = CountPreFSIParticles(gev);
if (fsparts.NOther ||
((fsparts.NPip + fsparts.NPi0 + fsparts.NPim) > 1)) {
return Channel_t::kNCTransitionMPi;
} else if (fsparts.NPip) {
return Channel_t::kNCSPP_NPip;
} else if (fsparts.NPi0) {
return fsparts.NProton ? Channel_t::kNCSPP_PPi0
: Channel_t::kNCSPP_NPi0;
} else if (fsparts.NPim) {
return Channel_t::kNCSPP_PPim;
}
return Channel_t::kNCTransitionMPi;
} else if (gev.Summary()->ProcInfo().IsDeepInelastic()) {
return Channel_t::kNCDIS;
} else {
std::cout << "Unknown GENIE Electron Scattering Mode!" << std::endl
<< "ScatteringTypeId = "
<< gev.Summary()->ProcInfo().ScatteringTypeId() << " "
<< "InteractionTypeId = "
<< gev.Summary()->ProcInfo().InteractionTypeId() << std::endl
<< genie::ScatteringType::AsString(
gev.Summary()->ProcInfo().ScatteringTypeId())
<< " "
<< genie::InteractionType::AsString(
gev.Summary()->ProcInfo().InteractionTypeId())
<< " " << gev.Summary()->ProcInfo().IsMEC() << std::endl;
return Channel_t::kUndefined;
}
}
// Weak CC
} else if (gev.Summary()->ProcInfo().IsWeakCC()) {
// CC MEC
if (gev.Summary()->ProcInfo().IsMEC()) {
if (utility::IsNeutralLepton(gev.Summary()->InitState().ProbePdg(),
utility::pdgcodes::kMatter)) {
return Channel_t::kCC2p2h;
} else if (utility::IsNeutralLepton(gev.Summary()->InitState().ProbePdg(),
utility::pdgcodes::kAntimatter)) {
return Channel_t::kCC2p2h_nub;
}
// CC OTHER
} else {
return FromNEUTCode(genie::utils::ghep::NeutReactionCode(&gev));
}
// Weak NC
} else if (gev.Summary()->ProcInfo().IsWeakNC()) {
// NC MEC
if (gev.Summary()->ProcInfo().IsMEC()) {
if (utility::IsNeutralLepton(gev.Summary()->InitState().ProbePdg(),
utility::pdgcodes::kMatter)) {
return Channel_t::kNC2p2h;
} else if (utility::IsNeutralLepton(gev.Summary()->InitState().ProbePdg(),
utility::pdgcodes::kAntimatter)) {
return Channel_t::kNC2p2h_nub;
}
// NC OTHER
} else {
return FromNEUTCode(genie::utils::ghep::NeutReactionCode(&gev));
}
}
return Channel_t::kUndefined;
}
Particle::Status_t GetParticleStatus(genie::GHepParticle const &p,
- Channel_t chan) {
+ Channel_t chan) {
/*
kIStUndefined = -1,
kIStInitialState = 0, / generator-level initial state /
kIStStableFinalState = 1, / generator-level final state:
particles to be tracked by detector-level MC /
kIStIntermediateState = 2,
kIStDecayedState = 3,
kIStCorrelatedNucleon = 10,
kIStNucleonTarget = 11,
kIStDISPreFragmHadronicState = 12,
kIStPreDecayResonantState = 13,
kIStHadronInTheNucleus = 14, / hadrons inside the nucleus: marked
for hadron transport modules to act on /
kIStFinalStateNuclearRemnant = 15, / low energy nuclear fragments
entering the record collectively as a 'hadronic blob' pseudo-particle /
kIStNucleonClusterTarget = 16, // for composite nucleons before
phase space decay
*/
Particle::Status_t state = Particle::Status_t::kUnknown;
switch (p.Status()) {
case genie::kIStNucleonTarget:
case genie::kIStInitialState:
case genie::kIStCorrelatedNucleon:
case genie::kIStNucleonClusterTarget: {
state = Particle::Status_t::kPrimaryInitialState;
break;
}
case genie::kIStStableFinalState: {
state = Particle::Status_t::kNuclearLeaving;
break;
}
case genie::kIStHadronInTheNucleus: {
- state = Is2p2h(chan) ? Particle::Status_t::kPrimaryInitialState
- : Particle::Status_t::kIntermediate;
+ state = utility::Is2p2h(chan) ? Particle::Status_t::kPrimaryInitialState
+ : Particle::Status_t::kIntermediate;
break;
}
case genie::kIStPreDecayResonantState:
case genie::kIStDISPreFragmHadronicState:
case genie::kIStIntermediateState: {
state = Particle::Status_t::kIntermediate;
break;
}
case genie::kIStFinalStateNuclearRemnant:
case genie::kIStUndefined:
case genie::kIStDecayedState:
default: { state = Particle::Status_t::kUnknown; }
}
if (utility::IsNuclearPDG(p.Pdg())) {
if (state == Particle::Status_t::kPrimaryInitialState) {
state = Particle::Status_t::kPrimaryInitialState;
} else if (state == Particle::Status_t::kNuclearLeaving) {
state = Particle::Status_t::kPrimaryFinalState;
}
}
return state;
}
+struct TargetSplineBlob {
+ double MolecularWeight;
+ size_t NNucleons;
+ std::string search_term;
+ TGraph TotSpline;
+};
+
+NEW_NUIS_EXCEPT(Invalid_target_specifier);
+
+double GetFileWeight(fhicl::ParameterSet const &xsecinfo) {
+ double mec_scale = xsecinfo.get("mec_scale", 1);
+ double EMin = xsecinfo.get("EMin", 0);
+ double EMax = xsecinfo.get("EMax", 10);
+ size_t NKnots = xsecinfo.get("NKnots", 100);
+
+ std::vector TargetSplines;
+ int probe = xsecinfo.get("nu_pdg");
+
+ for (std::string const &target :
+ xsecinfo.get>("target")) {
+ std::vector splits = nuis::utility::split(target, ";");
+ if (splits.size() != 1 && splits.size() != 2) {
+ throw Invalid_target_specifier()
+ << "Expected to find 100ZZZAAA0;, e.g. "
+ "1000060120;1, "
+ "but found: \""
+ << target << "\"";
+ }
+ std::stringstream ss;
+ ss << ".*nu:" << probe << ";tgt:" << splits[0] << ".*";
+
+ size_t nuc_pdg = fhicl::string_parsers::str2T(splits[0]);
+
+ double molwght = splits.size() == 2
+ ? fhicl::string_parsers::str2T(splits[1])
+ : 1;
+
+ TargetSplines.push_back(
+ {molwght, (nuc_pdg % 1000) / 10, ss.str(), TGraph()});
+
+ std::cout << "[INFO]:\tGetting splines that match: \"" << ss.str()
+ << "\" with A = " << TargetSplines.back().NNucleons
+ << " and multiplied by "
+ << molwght / double(TargetSplines.back().NNucleons) << std::endl;
+ }
+
+ size_t NTotNucleons = 0;
+ std::vector spline_patterns;
+ for (auto const &spb : TargetSplines) {
+ spline_patterns.push_back(spb.search_term);
+ NTotNucleons += spb.MolecularWeight * spb.NNucleons;
+ }
+ double WeightToPerNucleon = 1.0 / double(NTotNucleons);
+
+ // Read the spline file
+ std::cout << "[INFO]: Reading input file: "
+ << xsecinfo.get("spline_file")
+ << ", this may take a while..." << std::endl;
+ GENIESplineGetter spg(spline_patterns);
+ spg.ParseFile(xsecinfo.get("spline_file").c_str());
+ std::vector> all_splines = spg.GetTGraphs();
+
+ double step = (EMax - EMin) / double(NKnots);
+
+ // Sum each targets-worth of interaction splines together, weighting for
+ // molecular fraction and possibly fixing the MEC spline oddity
+ for (size_t ts_it = 0; ts_it < TargetSplines.size(); ++ts_it) {
+
+ TargetSplines[ts_it].TotSpline = TGraph(NKnots);
+
+ for (Int_t p_it = 0; p_it < NKnots; ++p_it) {
+ TargetSplines[ts_it].TotSpline.SetPoint(p_it, EMin + p_it * step, 0);
+ }
+
+ double target_weight = (TargetSplines[ts_it].MolecularWeight /
+ double(TargetSplines[ts_it].NNucleons));
+ for (size_t c_it = 0; c_it < all_splines[ts_it].size(); ++c_it) {
+ double weight = target_weight;
+ if (std::string(all_splines[ts_it][c_it].GetName()).find("MEC") !=
+ std::string::npos) {
+ if (mec_scale != 1) {
+ std::cout << "[INFO]: Weighting MEC splines by " << mec_scale
+ << std::endl;
+ weight *= mec_scale;
+ }
+ }
+ for (Int_t p_it = 0; p_it < NKnots; ++p_it) {
+ double E, XSec;
+ TargetSplines[ts_it].TotSpline.GetPoint(p_it, E, XSec);
+
+ // Copied from GENIE/Convents/Units.h
+ static const double gigaelectronvolt = 1.;
+ static const double GeV = gigaelectronvolt;
+ static const double meter = 5.07e+15 / GeV;
+ static const double centimeter = 0.01 * meter;
+ static const double centimeter2 = centimeter * centimeter;
+
+ XSec += all_splines[ts_it][c_it].Eval(E) * weight / centimeter2;
+
+ TargetSplines[ts_it].TotSpline.SetPoint(p_it, E, XSec);
+ }
+ }
+ }
+
+ // Sum all the correctly weighted per-target splines
+ TGraph MasterSpline(NKnots);
+
+ for (Int_t p_it = 0; p_it < NKnots; ++p_it) {
+ MasterSpline.SetPoint(p_it, EMin + p_it * step, 0);
+ }
+ for (size_t c_it = 0; c_it < TargetSplines.size(); ++c_it) {
+ for (Int_t p_it = 0; p_it < NKnots; ++p_it) {
+ double E, XSec;
+ MasterSpline.GetPoint(p_it, E, XSec);
+ XSec += TargetSplines[c_it].TotSpline.Eval(E) * WeightToPerNucleon;
+
+ MasterSpline.SetPoint(p_it, E, XSec);
+ }
+ }
+
+ fhicl::ParameterSet fluxps = xsecinfo.get("flux");
+ // If MonoE
+ if (fluxps.has_key("energy_GeV")) {
+ double monoE = fluxps.get("energy_GeV");
+
+ return MasterSpline.Eval(monoE);
+ } else {
+ // Use the master spline to get flux*xsec (which the events are thrown
+ // according to) and return the flux-averaged total xsec which is used for
+ // weighting
+ std::unique_ptr Flux =
+ nuis::utility::GetHistogram(fluxps.get("histogram"));
+ bool per_width = fluxps.get("is_divided_by_bin_width", true);
+ double FluxIntegral = Flux->Integral(per_width ? "width" : "");
+
+ for (int bi_it = 0; bi_it < Flux->GetXaxis()->GetNbins(); ++bi_it) {
+ double bc = Flux->GetBinContent(bi_it + 1);
+ double low_e = Flux->GetXaxis()->GetBinLowEdge(bi_it + 1);
+ double up_e = Flux->GetXaxis()->GetBinUpEdge(bi_it + 1);
+
+ size_t NIntSteps = 100;
+ double step = (up_e - low_e) / double(NIntSteps);
+
+ double avg_xsec = 0;
+ for (size_t i = 0; i < NIntSteps; ++i) {
+ double e = low_e + (double(i) + 0.5) * step;
+ avg_xsec += MasterSpline.Eval(e);
+ }
+ avg_xsec /= double(NIntSteps);
+
+ Flux->SetBinContent(bi_it + 1, bc * avg_xsec);
+ }
+
+ double EvRateIntegral = Flux->Integral(per_width ? "width" : "");
+
+ return EvRateIntegral / FluxIntegral;
+ }
+}
+
} // namespace genietools
} // namespace nuis
diff --git a/src/generator/utility/GENIEUtility.hxx b/src/generator/utility/GENIEUtility.hxx
index b95f450..3bc447a 100644
--- a/src/generator/utility/GENIEUtility.hxx
+++ b/src/generator/utility/GENIEUtility.hxx
@@ -1,47 +1,53 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#pragma once
#include "event/Particle.hxx"
#include "event/types.hxx"
class TGraph;
namespace genie {
class GHepRecord;
class GHepParticle;
} // namespace genie
+namespace fhicl {
+ class ParameterSet;
+}
+
namespace nuis {
namespace genietools {
NEW_NUIS_EXCEPT(invalid_GENIE_event);
TGraph const &GetGENIESpline(std::string const &SplineFile,
std::string const &SplineIdentifier);
event::Channel_t GetEventChannel(genie::GHepRecord const &);
event::Particle::Status_t GetParticleStatus(genie::GHepParticle const &p,
- event::Channel_t chan);
+ event::Channel_t chan);
+
+double GetFileWeight(fhicl::ParameterSet const &xsecinfo);
} // namespace genietools
} // namespace nuis
diff --git a/src/generator/utility/NEUTUtility.cxx b/src/generator/utility/NEUTUtility.cxx
index 5ff5b77..0b0e04a 100644
--- a/src/generator/utility/NEUTUtility.cxx
+++ b/src/generator/utility/NEUTUtility.cxx
@@ -1,95 +1,95 @@
#include "generator/utility/NEUTUtility.hxx"
#include "generator/input/NEUTInputHandler.hxx"
#include "exception/exception.hxx"
#include "utility/InteractionChannelUtility.hxx"
#include "utility/PDGCodeUtility.hxx"
#include "neutpart.h"
#include "neutvect.h"
using namespace nuis::event;
using namespace nuis::utility;
namespace nuis {
namespace neuttools {
NEW_NUIS_EXCEPT(unexpected_NEUT_particle_state);
Particle::Status_t GetNeutParticleStatus(NeutPart const &part, Channel_t mode) {
#ifdef DEBUG_NEUT_UTILITY
std::cout << "[DEBUG]: Mode: " << mode << ", Part: { Status: " << part.fStatus
<< ", IsAlive: " << part.fIsAlive << ", PDG: " << part.fPID << " }."
<< std::endl;
#endif
// Remove Pauli blocked events, probably just single pion events
if (part.fStatus == 5) {
return Particle::Status_t::kBlocked;
// fStatus == -1 means initial state
} else if (part.fIsAlive == false && part.fStatus == -1) {
return Particle::Status_t::kPrimaryInitialState;
// NEUT has a bit of a strange convention for fIsAlive and fStatus
// combinations
// for NC and neutrino particle isAlive true/false and status 2 means
// final state particle
// for other particles in NC status 2 means it's an FSI particle
// for CC it means it was an FSI particle
} else if (part.fStatus == 2) {
// NC case is a little strange... The outgoing neutrino might be alive or
// not alive. Remaining particles with status 2 are FSI particles that
// reinteracted
if (IsNC(mode) && IsNeutralLepton(part.fPID)) {
return Particle::Status_t::kNuclearLeaving;
// The usual CC case
} else if (part.fIsAlive == true) {
// return Particle::Status_t::kIntermediate;
throw unexpected_NEUT_particle_state()
<< "[ERROR] Found unexpected NEUT particle in neutvect stack: Mode: "
<< mode << " (IsNC: " << IsNC(mode)
<< "), Part: { Status: " << part.fStatus
<< ", IsAlive: " << part.fIsAlive << ", PDG: " << part.fPID
<< ", IsNeutralLepton: " << IsNeutralLepton(part.fPID) << " }.";
}
} else if ((part.fIsAlive == true) && (part.fStatus == 2) &&
IsNeutralLepton(part.fPID)) {
return Particle::Status_t::kNuclearLeaving;
} else if ((part.fIsAlive == true) && (part.fStatus == 0)) {
return Particle::Status_t::kNuclearLeaving;
} else if (!part.fIsAlive && ((part.fStatus == 1) || (part.fStatus == 3) ||
(part.fStatus == 4) || (part.fStatus == 7) ||
(part.fStatus == 8))) {
return Particle::Status_t::kIntermediate;
// There's one hyper weird case where fStatus = -3. This apparently
// corresponds to a nucleon being ejected via pion FSI when there is "data
// available"
} else if (!part.fIsAlive && (part.fStatus == -3)) {
return Particle::Status_t::kUnknown;
// NC neutrino outgoing
} else if (!part.fIsAlive && part.fStatus == 0 &&
IsNeutralLepton(part.fPID)) {
return Particle::Status_t::kNuclearLeaving;
// Warn if we still find alive particles without classifying them
} else if (part.fIsAlive == true) {
- throw unexpected_NEUT_particle_state()
- << "[ERROR]: Undefined NEUT state "
- << " Alive: " << part.fIsAlive << " Status: " << part.fStatus
- << " PDG: " << part.fPID;
+ std::cout << "[WARN]: Undefined NEUT state "
+ << " Alive: " << part.fIsAlive << " Status: " << part.fStatus
+ << " PDG: " << part.fPID << std::endl;
+ return Particle::Status_t::kUnknown;
}
// Warn if we find dead particles that we haven't classified
- throw unexpected_NEUT_particle_state()
- << "[ERROR]: Undefined NEUT state "
- << " Alive: " << part.fIsAlive << " Status: " << part.fStatus
- << " PDG: " << part.fPID;
+ std::cout << "[WARN]: Undefined NEUT state "
+ << " Alive: " << part.fIsAlive << " Status: " << part.fStatus
+ << " PDG: " << part.fPID << std::endl;
+ return Particle::Status_t::kUnknown;
}
} // namespace neuttools
} // namespace nuis
diff --git a/src/persistency/ROOTOutput.cxx b/src/persistency/ROOTOutput.cxx
index 2de36b1..34ac610 100644
--- a/src/persistency/ROOTOutput.cxx
+++ b/src/persistency/ROOTOutput.cxx
@@ -1,81 +1,88 @@
#include "persistency/ROOTOutput.hxx"
#include "utility/ROOTUtility.hxx"
#include "config/GlobalConfiguration.hxx"
#include "fhiclcpp/ParameterSet.h"
#include "TFile.h"
namespace nuis {
namespace persistency {
struct NamedTFile {
std::string name;
std::unique_ptr file;
NamedTFile() : name(""), file(nullptr) {}
NamedTFile(NamedTFile &&other)
: name(std::move(other.name)), file(std::move(other.file)) {}
~NamedTFile() {
if (file) {
file->Write();
file->Close();
}
}
};
static std::vector Files;
-std::unique_ptr &GetOutputFile(std::string const &name) {
- for (NamedTFile &file : Files) {
- if (file.name == name) {
- return file.file;
- }
- }
-
- fhicl::ParameterSet const &persistency =
- config::GetDocument().get("persistency");
- std::string file_name = persistency.get(name + ".output_file");
- std::string open_opts =
- persistency.get(name + ".open_mode", "CREATE");
-
+void NewStream(std::string const &name,
+ std::string file_name,
+ std::string opts) {
NamedTFile ntf;
ntf.name = name;
- ntf.file = std::make_unique(file_name.c_str(), open_opts.c_str());
+ ntf.file = std::make_unique(file_name.c_str(), opts.c_str());
if (!ntf.file || !ntf.file->IsOpen()) {
- if ((name == "default") && (open_opts == "CREATE")) {
+ if ((name == "default") && (opts == "CREATE")) {
std::cout
<< "[WARN]: It appears the default file cannot be opened because it "
"exists already, to stop you wasting processing time, the default "
"output stream will be written to nuis.default.tmp.root in the "
"current directory in RECREATE mode."
<< std::endl;
file_name = "nuis.default.tmp.root";
- open_opts = "RECREATE";
- ntf.file = std::make_unique(file_name.c_str(), open_opts.c_str());
+ opts = "RECREATE";
+ ntf.file = std::make_unique(file_name.c_str(), opts.c_str());
if (!ntf.file || !ntf.file->IsOpen()) {
throw utility::failed_to_open_TFile()
<< "[ERROR]: Failed to open output file: " << std::quoted(file_name)
- << " in write mode with opts = " << std::quoted(open_opts);
+ << " in write mode with opts = " << std::quoted(opts);
}
} else {
throw utility::failed_to_open_TFile()
<< "[ERROR]: Failed to open output file: " << std::quoted(file_name)
- << " in write mode with opts = " << std::quoted(open_opts);
+ << " in write mode with opts = " << std::quoted(opts);
}
}
Files.push_back(std::move(ntf));
- return Files.back().file;
+}
+
+std::unique_ptr &GetOutputFile(std::string const &name) {
+ for (NamedTFile &file : Files) {
+ if (file.name == name) {
+ return file.file;
+ }
+ }
+
+ fhicl::ParameterSet const &persistency =
+ config::GetDocument().get("persistency");
+ std::string file_name = persistency.get(name + ".output_file");
+ std::string opts =
+ persistency.get(name + ".open_mode", "CREATE");
+
+ NewStream(name, file_name, opts);
+
+ return GetOutputFile(name);
}
void CloseOpenTFiles() {
for (NamedTFile &f : Files) {
std::cout << "[INFO]: Closing open TFile: " << f.name << " "
<< f.file->GetName() << std::endl;
}
Files.clear();
}
} // namespace persistency
} // namespace nuis
diff --git a/src/persistency/ROOTOutput.hxx b/src/persistency/ROOTOutput.hxx
index 1aa5bb6..f88da22 100644
--- a/src/persistency/ROOTOutput.hxx
+++ b/src/persistency/ROOTOutput.hxx
@@ -1,76 +1,80 @@
#ifndef PERSITENCY_ROOTOUTPUT_HXX_SEEN
#define PERSITENCY_ROOTOUTPUT_HXX_SEEN
#include "exception/exception.hxx"
#include "TFile.h"
#include
#include
#include
namespace nuis {
namespace persistency {
NEW_NUIS_EXCEPT(WriteToOutputFile_nullptr);
/// Will get/open a TFile that is described in the global config
///
/// The named streams will be used to configure the file name and open mode from
/// the global config element persistency.: {file: output.root opts:
/// CREATE}
std::unique_ptr &GetOutputFile(std::string const &name = "default");
+void NewStream(std::string const &name = "default",
+ std::string file_name = "default.nuis.root",
+ std::string opts = "CREATE");
+
template
inline void WriteToOutputFile(T *object, std::string const &object_name,
std::string dir_name = "",
std::string const &file_name = "default") {
if (!object) {
throw WriteToOutputFile_nullptr();
}
TDirectory *ogdir = gDirectory;
std::unique_ptr &f = GetOutputFile(file_name);
TDirectory *d = f.get();
while (dir_name.length()) {
size_t next_slash = dir_name.find_first_of('/');
std::string next_dir = dir_name.substr(0, next_slash);
if (next_slash != std::string::npos) {
dir_name = dir_name.substr(next_slash + 1);
} else {
dir_name = "";
}
TDirectory *nd = d->GetDirectory(next_dir.c_str());
if (!nd) {
nd = d->mkdir(next_dir.c_str());
}
nd->cd();
d = nd;
}
d->cd();
object->Write(object_name.c_str(), TObject::kOverwrite);
if (ogdir) {
ogdir->cd();
}
}
template
inline void WriteToOutputFile(std::unique_ptr &object,
std::string const &object_name,
std::string dir_name = "",
std::string const &file_name = "default") {
return WriteToOutputFile(object.get(), object_name, dir_name, file_name);
}
void CloseOpenTFiles();
} // namespace persistency
} // namespace nuis
#endif
diff --git a/src/plugins/Instantiate.hxx b/src/plugins/Instantiate.hxx
index 818d778..f4ce598 100644
--- a/src/plugins/Instantiate.hxx
+++ b/src/plugins/Instantiate.hxx
@@ -1,252 +1,253 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#pragma once
#include "plugins/NamedSO.hxx"
#include "plugins/traits.hxx"
#include "config/GlobalConfiguration.hxx"
#include "utility/FileSystemUtility.hxx"
+#include "utility/StringUtility.hxx"
#include "exception/exception.hxx"
#include "fhiclcpp/ParameterSet.h"
#include "string_parsers/to_string.hxx"
// linux
#include
#include
#include
#include
#include
#include
#include
// #define DEBUG_INSTANTIATE
namespace nuis {
namespace plugins {
extern std::vector LoadedSharedObjects;
NEW_NUIS_EXCEPT(failed_to_find_instantiator);
NEW_NUIS_EXCEPT(malformed_plugin_interface);
NEW_NUIS_EXCEPT(failed_to_load_so);
typedef void *(*inst_fcn)();
typedef void (*dltr_fcn)(void *);
template struct PluginInstantiator {
std::string FQ_so_path;
std::string Base_classname;
std::string Classname;
void *dllib;
inst_fcn Instantiator;
dltr_fcn Deleter;
PluginInstantiator()
: FQ_so_path(""), Base_classname(""), Classname(""), dllib(nullptr),
Instantiator(nullptr), Deleter(nullptr) {}
PluginInstantiator(PluginInstantiator const &) = delete;
PluginInstantiator(PluginInstantiator &&other) {
FQ_so_path = std::move(other.FQ_so_path);
Base_classname = std::move(other.Base_classname);
Classname = std::move(other.Classname);
dllib = other.dllib;
Instantiator = other.Instantiator;
Deleter = other.Deleter;
other.FQ_so_path = "";
other.Base_classname = "";
other.Classname = "";
other.dllib = nullptr;
other.Instantiator = nullptr;
other.Deleter = nullptr;
}
typename plugin_traits::unique_ptr_t Instantiate() {
T *inst = reinterpret_cast((*Instantiator)());
dltr_fcn dltr = Deleter;
std::string cln = Classname;
std::function deleter = [=](T *inst) {
#ifdef DEBUG_INSTANTIATE
std::cout << "[INFO]: Deleting instance of " << cln << " with "
<< (void *)dltr << std::endl;
#endif
(*dltr)(inst);
};
return typename plugin_traits::unique_ptr_t(inst, deleter);
}
};
NamedSO &GetSharedObject(std::string const &FQPath) {
for (NamedSO &so : LoadedSharedObjects) {
if (so.name == FQPath) {
return so;
}
}
NamedSO so;
so.name = FQPath;
so.dllib = dlopen(FQPath.c_str(), RTLD_NOW | RTLD_GLOBAL);
char const *dlerr_cstr = dlerror();
std::string dlerr;
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr.length()) {
throw failed_to_load_so()
<< "[INFO]: Failed to load shared object: " << FQPath
<< " with dlerror: " << dlerr;
} else {
#ifdef DEBUG_INSTANTIATE
std::cout << "[INFO]: Loaded shared object " << FQPath << std::endl;
#endif
}
LoadedSharedObjects.push_back(std::move(so));
return LoadedSharedObjects.back();
}
template
typename plugin_traits::unique_ptr_t
Instantiate(std::string const &classname) {
static std::vector> LoadedPlugins;
fhicl::ParameterSet const &plugins =
config::GetDocument().get("plugins");
fhicl::ParameterSet const &search_paths =
plugins.get("search_paths");
std::vector plugin_search_dirs;
// Look for plugin search paths in sequence elements of the
// plugins.search_paths table
for (std::string const &key : search_paths.get_names()) {
if (!search_paths.is_key_to_sequence(key)) {
continue;
}
for (std::string const &path :
search_paths.get>(key)) {
plugin_search_dirs.push_back(path);
}
}
for (std::string path : plugin_search_dirs) {
path = utility::EnsureTrailingSlash(path);
for (std::string const &so_name :
utility::GetMatchingFiles(path, ".*\\.so")) {
for (PluginInstantiator &plugin : LoadedPlugins) {
if (plugin.FQ_so_path == (path + so_name) &&
(plugin.Base_classname == plugin_traits::interface_name()) &&
(plugin.Classname == classname)) {
#ifdef DEBUG_INSTANTIATE
std::cout << "[INFO]: Using already loaded PluginInstantiator"
<< std::endl;
#endif
return plugin.Instantiate();
}
}
PluginInstantiator plugin;
plugin.FQ_so_path = path + so_name;
plugin.Base_classname = plugin_traits::interface_name();
plugin.Classname = classname;
plugin.dllib = GetSharedObject(plugin.FQ_so_path).dllib;
char const *dlerr_cstr = nullptr;
std::string dlerr("");
plugin.Instantiator = reinterpret_cast(dlsym(
plugin.dllib,
plugin_traits::instantiator_function_name(classname).c_str()));
dlerr_cstr = dlerror();
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr_cstr) {
#ifdef DEBUG_INSTANTIATE
std::cout << "[INFO]: Failed to load appropriate instantiator method: "
<< plugin_traits::instantiator_function_name(classname)
<< " from shared object " << plugin.FQ_so_path << std::endl;
#endif
continue;
} else {
#ifdef DEBUG_INSTANTIATE
std::cout << "[INFO]: Loaded instantiator method: "
<< plugin_traits::instantiator_function_name(classname)
<< " from shared object " << plugin.FQ_so_path << std::endl;
#endif
}
plugin.Deleter = reinterpret_cast(
dlsym(plugin.dllib,
plugin_traits::deleter_function_name(classname).c_str()));
dlerr_cstr = dlerror();
if (dlerr_cstr) {
dlerr = dlerr_cstr;
}
if (dlerr_cstr) {
throw malformed_plugin_interface()
<< "[ERROR]: Failed to load appropriate deleter method: "
<< plugin_traits::deleter_function_name(classname)
<< " from shared object " << plugin.FQ_so_path
<< " with error: " << std::quoted(dlerr);
} else {
#ifdef DEBUG_INSTANTIATE
std::cout << "[INFO]: Loaded deleter method: "
<< plugin_traits::deleter_function_name(classname)
<< " from shared object " << plugin.FQ_so_path << std::endl;
#endif
}
#ifdef DEBUG_INSTANTIATE
std::cout << "[INFO]: Checking if shared object "
<< std::quoted(plugin.FQ_so_path)
<< " knows how to instantiate class " << std::quoted(classname)
<< " via interface "
<< std::quoted(plugin_traits::interface_name()) << std::endl;
#endif
LoadedPlugins.push_back(std::move(plugin));
return LoadedPlugins.back().Instantiate();
}
}
throw failed_to_find_instantiator()
<< "[ERROR]: Failed to find instantiator for classname: "
<< std::quoted(classname) << " using interface "
<< std::quoted(plugin_traits::interface_name())
<< " from configured search paths: "
<< fhicl::string_parsers::T2Str>(
plugin_search_dirs);
}
} // namespace plugins
} // namespace nuis
diff --git a/src/samples/CMakeLists.txt b/src/samples/CMakeLists.txt
index 6313da0..c1a9654 100644
--- a/src/samples/CMakeLists.txt
+++ b/src/samples/CMakeLists.txt
@@ -1,24 +1,25 @@
set(samples_header_files
IEventProcessor.hxx
IDataComparison.hxx
SimpleDataComparison.hxx
+ MultiDataComparison.hxx
SimpleMCStudy.hxx)
install(FILES ${samples_header_files} DESTINATION include/samples)
add_subdirectory(MCTools)
add_subdirectory(nuA)
cmessage(DEBUG "INuADataComparisons: ${INuADataComparisons}")
SET(INuADataComparisons_List)
if(NOT IDataComparisons STREQUAL "")
string(REPLACE ";" ", " INuADataComparisons_List "${INuADataComparisons}")
endif()
cmessage(DEBUG "INuADataComparisons_List: ${INuADataComparisons_List}")
SET(INuADataComparisons_List ${INuADataComparisons_List} PARENT_SCOPE)
SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE)
diff --git a/src/samples/IEventProcessor.hxx b/src/samples/IEventProcessor.hxx
index 90eff1a..a3693ce 100644
--- a/src/samples/IEventProcessor.hxx
+++ b/src/samples/IEventProcessor.hxx
@@ -1,119 +1,130 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef SAMPLES_IEventProcessor_HXX_SEEN
#define SAMPLES_IEventProcessor_HXX_SEEN
#include "plugins/traits.hxx"
#include "exception/exception.hxx"
#include "config/GlobalConfiguration.hxx"
#include "fhiclcpp/ParameterSet.h"
#include "TH1.h"
#include
#include
#include
namespace nuis {
namespace event {
class FullEvent;
class MinimalEvent;
} // namespace event
} // namespace nuis
#define IEventProcessor_DEBUG(v) \
if (verb > 2) { \
std::cout << "[DEBUG]: " << v << std::endl; \
}
#define IEventProcessor_INFO(v) \
if (verb > 1) { \
std::cout << "[INFO]: " << v << std::endl; \
}
#define IEventProcessor_WARN(v) \
if (verb > 0) { \
std::cout << "[WARN] " << __FILENAME__ << ":" << __LINE__ << " : " << v \
<< std::endl; \
}
class IEventProcessor {
protected:
NEW_NUIS_EXCEPT(unknown_IEventProcessor_verbosity);
enum sample_verbosity { kSilent = 0, kWarn = 1, kReticent = 2, kVerbose = 3 };
sample_verbosity verb;
void SetSampleVerbosity(std::string v) {
v = nuis::config::GetDocument().get(
"global.sample.verbosity_override", v);
if (v == "Silent") {
verb = kSilent;
} else if (v == "Warn") {
verb = kWarn;
} else if (v == "Reticent") {
verb = kReticent;
} else if (v == "Verbose") {
verb = kVerbose;
} else {
throw unknown_IEventProcessor_verbosity()
<< "[ERROR]: Failed to parse IEventProcessor verbosity setting from: "
<< std::quoted(v);
}
}
public:
NEW_NUIS_EXCEPT(uninitialized_IEventProcessor);
NEW_NUIS_EXCEPT(unimplemented_IEventProcessor_optional_method);
IEventProcessor() {
SetSampleVerbosity(nuis::config::GetDocument().get(
"global.sample.verbosity_default", "Silent"));
TH1::SetDefaultSumw2();
}
virtual void Initialize(fhicl::ParameterSet const &) = 0;
// Interface for processing a single, external event
//
// IEventProcessors are not required to implement processing events from
// 'outside'.
- virtual void ProcessEvent(nuis::event::FullEvent const &) {
+ virtual void ProcessEvent(nuis::event::FullEvent const &, double weight = 1) {
throw unimplemented_IEventProcessor_optional_method()
<< "[ERROR]: IEventProcessor " << Name()
<< " does not implement ProcessEvent.";
}
+ // Interface for processing a pre-selected single, external event
+ //
+ // IEventProcessors are not required to implement processing events from
+ // 'outside'.
+ virtual void ProcessSignalEvent(nuis::event::FullEvent const &,
+ double weight = 1) {
+ throw unimplemented_IEventProcessor_optional_method()
+ << "[ERROR]: IEventProcessor " << Name()
+ << " does not implement ProcessSignalEvent.";
+ }
+
virtual void
ProcessSample(size_t nmax = std::numeric_limits::max()) = 0;
virtual void Write() = 0;
virtual std::string Name() = 0;
virtual ~IEventProcessor() {}
};
DECLARE_PLUGIN_INTERFACE(IEventProcessor);
#endif
diff --git a/src/samples/MCTools/CMakeLists.txt b/src/samples/MCTools/CMakeLists.txt
index 6d1455b..b325ae9 100644
--- a/src/samples/MCTools/CMakeLists.txt
+++ b/src/samples/MCTools/CMakeLists.txt
@@ -1,16 +1,16 @@
if(USE_NUWRO)
LIST(APPEND MC_TOOL_IMPL NuisToNuWro.cxx)
endif(USE_NUWRO)
-LIST(APPEND MC_TOOL_IMPL VerboseEventSummary.cxx NuisToCAFTruth)
+LIST(APPEND MC_TOOL_IMPL VerboseEventSummary.cxx NuisToCAFTruth.cxx EventSummary_ECTJune2019.cxx)
add_library(MCTools SHARED ${MC_TOOL_IMPL})
target_link_libraries(MCTools ${SAMPLES_MCTOOLS_LIBS})
if(USE_NUWRO)
target_compile_options(MCTools PUBLIC ${NUWRO_CXX_FLAGS})
include_directories(${NUWRO_INCLUDE_DIRS})
target_link_libraries(MCTools -Wl,--no-as-needed ${NUWRO_IMPORTED_TARGETS})
endif()
install(TARGETS MCTools DESTINATION plugins)
diff --git a/src/samples/MCTools/EventSummary_ECTJune2019.cxx b/src/samples/MCTools/EventSummary_ECTJune2019.cxx
new file mode 100644
index 0000000..d029fc5
--- /dev/null
+++ b/src/samples/MCTools/EventSummary_ECTJune2019.cxx
@@ -0,0 +1,311 @@
+// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
+
+/*******************************************************************************
+ * This file is part of NUISANCE.
+ *
+ * NUISANCE is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * NUISANCE is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with NUISANCE. If not, see .
+ *******************************************************************************/
+
+//********************************************************************
+
+#include "samples/SimpleMCStudy.hxx"
+
+#include "utility/EventTopologyUtility.hxx"
+#include "utility/FullEventUtility.hxx"
+#include "utility/InteractionChannelUtility.hxx"
+#include "utility/KinematicUtility.hxx"
+#include "utility/PDGCodeUtility.hxx"
+#include "utility/ROOTUtility.hxx"
+
+#include "utility/experimental/MINERvAUtility.hxx"
+#include "utility/experimental/T2KUtility.hxx"
+
+#include "persistency/ROOTOutput.hxx"
+
+using namespace nuis::event;
+using namespace nuis::utility;
+
+class EventSummary_ECTJune2019 : public SimpleMCStudy {
+
+ TreeFile EvSumTree;
+
+ struct EvSum {
+ double e_nu;
+ int pdg_nu;
+
+ double e_fslep;
+ double costheta_fslep;
+ int pdg_fslep;
+
+ double e_hmfsproton;
+ double costheta_hmfsproton;
+
+ double e_hmfspi;
+ double costheta_hmfspi;
+ int pdg_hmfspi;
+
+ double dpt;
+ double dat;
+ double dphit;
+ double dpt_mnv;
+ double dat_mnv;
+ double dphit_mnv;
+ double pn_mnv;
+ double q0;
+ double q3;
+ double dpn;
+ double EAvailProxy;
+ double Q2;
+ double W_rec;
+
+ int mode;
+
+ int NFSPip;
+ int NFSPim;
+ int NFSPi0;
+ int NFSProton;
+ int NFSNeutron;
+ int NFSGamma;
+ int NFSKaon;
+ int NFSOther;
+
+ bool isCC;
+
+ bool isQE;
+ bool isQEL;
+ bool is2p2h;
+ bool isNucleonSPP;
+ bool isCoh;
+ bool isMultiPi;
+ bool isDIS;
+ bool isOtherChannel;
+
+ bool is0Pi;
+ bool is1Pi;
+ bool isNPi;
+ bool isOtherTopo;
+
+ bool isT2K0Pi;
+ bool isT2K1Pip_CH;
+ bool isT2K_STV;
+
+ bool isMINERvA0Pi;
+ bool isMINERvA1CPi;
+ bool isMINERvA1Pi0;
+ bool isMINERvALowRecoil;
+ bool isMINERvASTV;
+
+ double xsweight;
+ };
+
+ EvSum es;
+
+public:
+ EventSummary_ECTJune2019() { ReadGlobalConfigDefaults(); }
+
+ void SetBranchAddresses() {
+ EvSumTree->Branch("e_nu", &es.e_nu, "e_nu/D");
+ EvSumTree->Branch("pdg_nu", &es.pdg_nu, "pdg_nu/I");
+
+ EvSumTree->Branch("e_fslep", &es.e_fslep, "e_fslep/D");
+ EvSumTree->Branch("costheta_fslep", &es.costheta_fslep, "costheta_fslep/D");
+ EvSumTree->Branch("pdg_fslep", &es.pdg_fslep, "pdg_fslep/I");
+
+ EvSumTree->Branch("e_hmfsproton", &es.e_hmfsproton, "e_hmfsproton/D");
+ EvSumTree->Branch("costheta_hmfsproton", &es.costheta_hmfsproton,
+ "costheta_hmfsproton/D");
+
+ EvSumTree->Branch("e_hmfspi", &es.e_hmfspi, "e_hmfspi/D");
+ EvSumTree->Branch("costheta_hmfspi", &es.costheta_hmfspi,
+ "costheta_hmfspi/D");
+ EvSumTree->Branch("pdg_hmfspi", &es.pdg_hmfspi, "pdg_hmfspi/I");
+
+ EvSumTree->Branch("dpt", &es.dpt, "dpt/D");
+ EvSumTree->Branch("dat", &es.dat, "dat/D");
+ EvSumTree->Branch("dphit", &es.dphit, "dphit/D");
+ EvSumTree->Branch("dpt_mnv", &es.dpt_mnv, "dpt_mnv/D");
+ EvSumTree->Branch("dat_mnv", &es.dat_mnv, "dat_mnv/D");
+ EvSumTree->Branch("dphit_mnv", &es.dphit_mnv, "dphit_mnv/D");
+ EvSumTree->Branch("pn_mnv", &es.pn_mnv, "pn_mnv/D");
+ EvSumTree->Branch("q0", &es.q0, "q0/D");
+ EvSumTree->Branch("q3", &es.q3, "q3/D");
+ EvSumTree->Branch("EAvailProxy", &es.EAvailProxy, "EAvailProxy/D");
+ EvSumTree->Branch("Q2", &es.Q2, "Q2/D");
+ EvSumTree->Branch("W_rec", &es.W_rec, "W_rec/D");
+
+ EvSumTree->Branch("mode", &es.mode, "mode/I");
+
+ EvSumTree->Branch("NFSPip", &es.NFSPip, "NFSPip/I");
+ EvSumTree->Branch("NFSPim", &es.NFSPim, "NFSPim/I");
+ EvSumTree->Branch("NFSPi0", &es.NFSPi0, "NFSPi0/I");
+ EvSumTree->Branch("NFSProton", &es.NFSProton, "NFSProton/I");
+ EvSumTree->Branch("NFSNeutron", &es.NFSNeutron, "NFSNeutron/I");
+ EvSumTree->Branch("NFSGamma", &es.NFSGamma, "NFSGamma/I");
+ EvSumTree->Branch("NFSKaon", &es.NFSKaon, "NFSKaon/I");
+ EvSumTree->Branch("NFSOther", &es.NFSOther, "NFSOther/I");
+
+ EvSumTree->Branch("isCC", &es.isCC, "isCC/O");
+
+ EvSumTree->Branch("isQE", &es.isQE, "isQE/O");
+ EvSumTree->Branch("isQEL", &es.isQEL, "isQEL/O");
+ EvSumTree->Branch("is2p2h", &es.is2p2h, "is2p2h/O");
+ EvSumTree->Branch("isNucleonSPP", &es.isNucleonSPP, "isNucleonSPP/O");
+ EvSumTree->Branch("isCoh", &es.isCoh, "isCoh/O");
+ EvSumTree->Branch("isMultiPi", &es.isMultiPi, "isMultiPi/O");
+ EvSumTree->Branch("isDIS", &es.isDIS, "isDIS/O");
+ EvSumTree->Branch("isOtherChannel", &es.isOtherChannel, "isOtherChannel/O");
+
+ EvSumTree->Branch("is0Pi", &es.is0Pi, "is0Pi/O");
+ EvSumTree->Branch("is1Pi", &es.is1Pi, "is1Pi/O");
+ EvSumTree->Branch("isNPi", &es.isNPi, "isNPi/O");
+ EvSumTree->Branch("isOtherTopo", &es.isOtherTopo, "isOtherTopo/O");
+
+ EvSumTree->Branch("isT2K0Pi", &es.isT2K0Pi, "isT2K0Pi/O");
+ EvSumTree->Branch("isT2K1Pip_CH", &es.isT2K1Pip_CH, "isT2K1Pip_CH/O");
+ EvSumTree->Branch("isT2K_STV", &es.isT2K_STV, "isT2K_STV/O");
+
+ EvSumTree->Branch("isMINERvA0Pi", &es.isMINERvA0Pi, "isMINERvA0Pi/O");
+ EvSumTree->Branch("isMINERvA1CPi", &es.isMINERvA1CPi, "isMINERvA1CPi/O");
+ EvSumTree->Branch("isMINERvA1Pi0", &es.isMINERvA1Pi0, "isMINERvA1Pi0/O");
+ EvSumTree->Branch("isMINERvALowRecoil", &es.isMINERvALowRecoil,
+ "isMINERvALowRecoil/O");
+ EvSumTree->Branch("isMINERvASTV", &es.isMINERvASTV, "isMINERvASTV/O");
+
+ EvSumTree->Branch("xsweight", &es.xsweight, "xsweight/D");
+ }
+
+ void Initialize(fhicl::ParameterSet const &instance_sample_configuration) {
+
+ // Perform any per-sample configuration in the base class
+ SimpleMCStudy::Initialize(instance_sample_configuration);
+
+ EvSumTree = AddNewTTreeToFile(nuis::persistency::GetOutputFile().get(),
+ "ECTEvSumTree", write_directory);
+ SetBranchAddresses();
+
+ //! Here you do whatever you want to do on a per-event basis.
+ //! This method will be run for every event given by the input handler, the
+ //! weight allows differential xsecs to be plotted easily and also contains
+ //! any requested reweightable parameter variations.
+ ProcessEventFunction = [&](nuis::event::FullEvent const &ev,
+ double weight) -> void {
+ Particle ISNu = GetHMISNeutralLepton(ev);
+ if (!ISNu) {
+ return;
+ }
+
+ Particle FSLep = GetHMFSLepton(ev);
+ if (!FSLep) {
+ return;
+ }
+
+ es.e_nu = ISNu.E();
+ es.pdg_nu = ISNu.pdg;
+
+ es.e_fslep = FSLep.E();
+ es.costheta_fslep = FSLep.CosTheta();
+ es.pdg_fslep = FSLep.pdg;
+
+ Particle FSPi = GetHMFSPion(ev);
+ if (!FSPi) {
+ es.e_hmfspi = 0;
+ es.costheta_hmfspi = 0;
+ es.pdg_hmfspi = 0;
+ } else {
+ es.e_hmfspi = FSPi.E();
+ es.costheta_hmfspi = FSPi.CosTheta();
+ es.pdg_hmfspi = FSPi.pdg;
+ }
+
+ Particle FSProton = GetHMFSProton(ev);
+ if (!FSProton) {
+ es.e_hmfsproton = 0;
+ es.costheta_hmfsproton = 0;
+ } else {
+ es.e_hmfsproton = FSProton.E();
+ es.costheta_hmfsproton = FSProton.CosTheta();
+ }
+
+ es.dpt = GetDeltaPT_CC0PiN(ev, ISNu.pdg).Mag();
+ es.dat = GetDeltaPhiT_CC0PiN(ev, ISNu.pdg);
+ es.dphit = GetDeltaAlphaT_CC0PiN(ev, ISNu.pdg);
+ es.dpt_mnv = mnv::GetDeltaPT_CC0PiN_mnv(ev).Mag();
+ es.dat_mnv = mnv::GetDeltaAlphaT_CC0PiN_mnv(ev);
+ es.dphit_mnv = mnv::GetDeltaPhiT_CC0PiN_mnv(ev);
+ es.pn_mnv = mnv::GetNeutronMomentumReco_CC0PiN_mnv(ev);
+ es.EAvailProxy = GetEAvailProxy(ev);
+ TLorentzVector EMTransfer_lep = GetEnergyMomentumTransfer(ev);
+ es.Q2 = -EMTransfer_lep.Mag2();
+ es.q0 = EMTransfer_lep.E();
+ es.q3 = EMTransfer_lep.Vect().Mag();
+ es.W_rec = GetNeutrinoWRec(ev);
+
+ es.mode = ChannelToInt(ev.mode);
+
+ es.NFSPip = GetNParticles(ev, {pdgcodes::kPiPlus});
+ es.NFSPim = GetNParticles(ev, {pdgcodes::kPiMinus});
+ es.NFSPi0 = GetNParticles(ev, {pdgcodes::kPi0});
+ es.NFSProton = GetNParticles(ev, {pdgcodes::kProton});
+ es.NFSNeutron = GetNParticles(ev, {pdgcodes::kNeutron});
+ es.NFSGamma = GetNParticles(ev, {pdgcodes::kGamma});
+ es.NFSKaon = GetNParticles(ev, pdgcodes::Kaons);
+ es.NFSOther = GetNFSOthers(ev) - (es.NFSGamma + es.NFSKaon);
+
+ es.isCC = IsCC(ev.mode);
+
+ es.isQE = IsQE(ev.mode);
+ es.isQEL = IsQEL(ev.mode);
+ es.is2p2h = Is2p2h(ev.mode);
+ es.isNucleonSPP = IsNucleonSPP(ev.mode);
+ es.isCoh = IsCoh(ev.mode);
+ es.isMultiPi = IsMultiPi(ev.mode);
+ es.isDIS = IsDIS(ev.mode);
+ es.isOtherChannel =
+ !(es.isQE || es.isQEL || es.is2p2h || es.isNucleonSPP || es.isCoh ||
+ es.isMultiPi || es.isDIS);
+
+ size_t NPi = GetNFSPions(ev);
+ size_t NOthers = GetNFSOthers(ev);
+ es.is0Pi = (NPi == 0) && (!NOthers);
+ es.is1Pi = (NPi == 1) && (!NOthers);
+ es.isNPi = (NPi > 1) && (!NOthers);
+ es.isOtherTopo = NOthers;
+
+ es.isT2K0Pi =
+ t2k::IsCC0Pi_NumProtons(ev).first; // Don't care about n protons
+ es.isT2K1Pip_CH = t2k::IsCC1Pip_CH_RecPi(ev);
+ es.isT2K_STV = t2k::IsCC0Pi_STV(ev);
+
+ es.isMINERvA0Pi =
+ mnv::IsCC0Pi_NumProtons(ev).first; // Don't care about n protons
+ es.isMINERvA1CPi = mnv::IsCC1CPi_2017(ev);
+ es.isMINERvA1Pi0 = mnv::IsCC1Pi0_2016(ev);
+ es.isMINERvALowRecoil = mnv::IsCCIncLowRecoil(ev);
+ es.isMINERvASTV = mnv::IsCC0PiNp_STV(ev);
+
+ es.xsweight = weight;
+
+ EvSumTree->Fill();
+ };
+ }
+
+ std::string Name() { return "EventSummary_ECTJune2019"; }
+
+ //! Here you can write any custom histograms to TTrees that your sample has
+ //! been handling.
+ void Write() {}
+};
+
+//! These declarations allow your class to be loaded dynamically by NUISANCE
+DECLARE_PLUGIN(IEventProcessor, EventSummary_ECTJune2019);
diff --git a/src/samples/MultiDataComparison.hxx b/src/samples/MultiDataComparison.hxx
new file mode 100644
index 0000000..cb9472c
--- /dev/null
+++ b/src/samples/MultiDataComparison.hxx
@@ -0,0 +1,231 @@
+// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
+
+/*******************************************************************************
+ * This file is part of NUISANCE.
+ *
+ * NUISANCE is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * NUISANCE is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with NUISANCE. If not, see .
+ *******************************************************************************/
+
+#pragma once
+
+#include "samples/IDataComparison.hxx"
+
+#include "event/FullEvent.hxx"
+
+#include "input/InputManager.hxx"
+
+#include "persistency/ROOTOutput.hxx"
+
+#include "utility/FileSystemUtility.hxx"
+#include "utility/HistogramUtility.hxx"
+#include "utility/KinematicUtility.hxx"
+#include "utility/ROOTUtility.hxx"
+#include "utility/StatsUtility.hxx"
+
+#include
+#include
+#include
+#include
+#include
+
+///Wraps multiple projections of the same signal selection
+class MultiDataComparison : public IDataComparison {
+
+ NEW_NUIS_EXCEPT(invalid_MultiDataComparison_initialization);
+ NEW_NUIS_EXCEPT(MultiDataComparison_already_finalized);
+
+protected:
+ std::string fName;
+
+ nuis::input::InputManager::Input_id_t fIH_id;
+ std::string write_directory;
+ size_t NMaxSample_override;
+
+ std::string fJournalReference;
+ std::string fDOI;
+ std::string fYear;
+ std::string fTargetMaterial;
+ std::string fFluxDescription;
+ std::string fSignalDescription;
+
+ std::function IsSigFunc;
+
+ std::vector>>
+ Comparisons;
+
+ nuis::utility::KinematicRange energy_cut;
+
+public:
+ MultiDataComparison(std::string name) {
+ fName = std::move(name);
+ fIH_id = std::numeric_limits::max();
+ write_directory = "";
+ NMaxSample_override = std::numeric_limits::max();
+
+ fJournalReference = "";
+ fDOI = "";
+ fYear = "";
+ fTargetMaterial = "";
+ fFluxDescription = "";
+ fSignalDescription = "";
+
+ energy_cut = nuis::utility::KinematicRange{
+ std::numeric_limits::max(), std::numeric_limits::max()};
+ }
+
+ std::string Name() { return fName; }
+
+ fhicl::ParameterSet fGlobalConfig;
+ fhicl::ParameterSet fInstanceConfig;
+
+ void ReadGlobalConfigDefaults() {
+ fGlobalConfig = nuis::config::GetDocument().get(
+ std::string("global.sample_configuration.") + Name(),
+ fhicl::ParameterSet());
+
+ if (!fJournalReference.length()) {
+ fJournalReference = fGlobalConfig.get("journal_reference",
+ "Not yet published");
+ }
+ if (!fDOI.length()) {
+ fDOI = fGlobalConfig.get("DOI", "Unknown DOI");
+ }
+ if (!fYear.length()) {
+ fYear = fGlobalConfig.get("year", "Unknown year");
+ }
+
+ if (!fTargetMaterial.length()) {
+ fTargetMaterial = fGlobalConfig.get(
+ "target_material", "Unknown target material");
+ }
+ if (!fFluxDescription.length()) {
+ fFluxDescription =
+ fGlobalConfig.get("flux_description", "Unknown flux");
+ }
+ if (!fSignalDescription.length()) {
+ fSignalDescription = fGlobalConfig.get("signal_description",
+ "Unknown Signal");
+ }
+
+ if ((energy_cut.first == std::numeric_limits::max()) &&
+ (fGlobalConfig.has_key("enu_range"))) {
+ energy_cut = fGlobalConfig.get>("enu_range");
+ }
+ }
+
+ virtual std::string GetJournalReference() { return fJournalReference; }
+ virtual std::string GetDOI() { return fDOI; }
+ virtual std::string GetYear() { return fYear; }
+ virtual std::string GetTargetMaterial() { return fTargetMaterial; }
+ virtual std::string GetFluxDescription() { return fFluxDescription; }
+ virtual std::string GetSignalDescription() { return fSignalDescription; }
+
+ void Initialize(fhicl::ParameterSet const &instance_sample_configuration) {
+
+ fInstanceConfig = instance_sample_configuration;
+
+ SetSampleVerbosity(fInstanceConfig.get(
+ "verbosity", nuis::config::GetDocument().get(
+ "global.sample.verbosity_default", "Reticent")));
+
+ ReadGlobalConfigDefaults();
+
+ if (fInstanceConfig.has_key("verbosity")) {
+ SetSampleVerbosity(fInstanceConfig.get("verbosity"));
+ }
+
+ NMaxSample_override =
+ fInstanceConfig.get("nmax", std::numeric_limits::max());
+
+ write_directory =
+ fInstanceConfig.get("write_directory", Name());
+ }
+
+ void ProcessSample(size_t nmax = std::numeric_limits::max()) {
+ if (fIH_id ==
+ std::numeric_limits::max()) {
+ fIH_id =
+ nuis::input::InputManager::Get().EnsureInputLoaded(fInstanceConfig);
+ }
+ IInputHandler const &IH =
+ nuis::input::InputManager::Get().GetInputHandler(fIH_id);
+
+ nmax = std::min(NMaxSample_override, nmax);
+
+ size_t NEvsToProcess = std::min(nmax, IH.GetNEvents());
+
+ double nmax_scaling = double(IH.GetNEvents()) / double(NEvsToProcess);
+
+ size_t NToShout = NEvsToProcess / 10;
+ IEventProcessor_INFO("Sample " << std::quoted(Name()) << " processing "
+ << NEvsToProcess << " events.");
+
+ IInputHandler::ev_index_t ev_idx = 0;
+
+ size_t NSig = 0;
+ while (ev_idx < NEvsToProcess) {
+ bool is_sig = false;
+
+ nuis::event::FullEvent const &fev = IH.GetFullEvent(ev_idx);
+
+ if (NToShout && !(ev_idx % NToShout)) {
+ IEventProcessor_INFO("\tProcessed " << ev_idx << "/" << NEvsToProcess
+ << " events.");
+ }
+
+ if (IsSigFunc(fev)) {
+ NSig++;
+ for (auto &comp : Comparisons) {
+ comp.second->ProcessSignalEvent(fev, IH.GetEventWeight(ev_idx) *
+ nmax_scaling);
+ }
+ }
+
+ ev_idx++;
+ }
+ IEventProcessor_INFO("\t" << NSig << "/" << NEvsToProcess
+ << " events passed selection.");
+ }
+ void Write() {}
+
+ double GetGOF() {
+ double totGOF = 0;
+ for (auto const &comp : Comparisons) {
+ double sGOF = comp.second->GetGOF();
+ if (sGOF != std::numeric_limits::max()) {
+ totGOF += sGOF;
+ } else {
+ std::cout << "[WARN]: Projection \"" << comp.first << "\" of sample "
+ << Name() << " produced bad GOF." << std::endl;
+ }
+ }
+ return totGOF;
+ }
+ double GetNDOGuess() {
+ size_t totNDOGuess = 0;
+ for (auto const &comp : Comparisons) {
+ totNDOGuess += comp.second->GetNDOGuess();
+ }
+ return totNDOGuess;
+ }
+
+ fhicl::ParameterSet GetExampleConfiguration() {
+ fhicl::ParameterSet exps;
+ exps.put("name", Name());
+ exps.put("input_type", "NuWro/NEUT/GENIE");
+ exps.put("file", "input_events.root");
+ exps.put("write_directory", Name());
+ return exps;
+ }
+};
diff --git a/src/samples/SimpleDataComparison.hxx b/src/samples/SimpleDataComparison.hxx
index b8f91e7..fbfb548 100644
--- a/src/samples/SimpleDataComparison.hxx
+++ b/src/samples/SimpleDataComparison.hxx
@@ -1,451 +1,493 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#pragma once
#include "samples/IDataComparison.hxx"
#include "event/FullEvent.hxx"
#include "input/InputManager.hxx"
#include "persistency/ROOTOutput.hxx"
#include "utility/FileSystemUtility.hxx"
#include "utility/HistogramUtility.hxx"
#include "utility/KinematicUtility.hxx"
#include "utility/ROOTUtility.hxx"
#include "utility/StatsUtility.hxx"
#include
#include
#include
#include
#include
template ::type>
class SimpleDataComparison : public IDataComparison {
NEW_NUIS_EXCEPT(invalid_SimpleDataComparison_initialization);
NEW_NUIS_EXCEPT(SimpleDataComparison_already_finalized);
+public:
+ static size_t const NDim = nd;
+ std::string write_directory;
+
protected:
using HistType = HT;
using TH_Help = typename nuis::utility::TH_Helper;
+ nuis::input::InputManager::Input_id_t fIH_id;
- static size_t const NDim = nd;
+ std::string fName;
- nuis::input::InputManager::Input_id_t fIH_id;
- std::string write_directory;
size_t NMaxSample_override;
int fIsShapeOnly;
int fIsFluxUnfolded;
std::vector fSignalCache;
std::vector> fProjectionCache;
bool fUseCache;
+ bool fModeHists;
std::string fDataInputDescriptor;
std::unique_ptr fData;
std::string fMaskInputDescriptor;
std::unique_ptr fMask;
std::string fCovarianceInputDescriptor;
std::unique_ptr fCovariance;
std::unique_ptr fPrediction;
+ std::map> fPrediction_modes;
std::unique_ptr fPrediction_xsec;
std::unique_ptr fPrediction_shape;
std::unique_ptr fPrediction_comparison;
bool fComparisonFinalized;
std::string fJournalReference;
std::string fDOI;
std::string fYear;
std::string fTargetMaterial;
std::string fFluxDescription;
std::string fSignalDescription;
- nuis::utility::ENuRange energy_cut;
+ nuis::utility::KinematicRange energy_cut;
std::function IsSigFunc;
- std::function(nuis::event::FullEvent const &)>
- CompProjFunc;
// If assigned by subclass will be called on for all events, bool signifies
// whether the event was selected.
std::function
ProcessExtraFunc;
public:
- SimpleDataComparison() {
+ std::function(nuis::event::FullEvent const &)>
+ CompProjFunc;
+
+ SimpleDataComparison(std::string name) {
+ fName = std::move(name);
fIH_id = std::numeric_limits::max();
fUseCache = false;
write_directory = "";
NMaxSample_override = std::numeric_limits::max();
fDataInputDescriptor = "";
fData = nullptr;
fMaskInputDescriptor = "";
fMask = nullptr;
fCovarianceInputDescriptor = "";
fCovariance = nullptr;
fPrediction = nullptr;
fPrediction_xsec = nullptr;
fPrediction_shape = nullptr;
fPrediction_comparison = nullptr;
fComparisonFinalized = false;
IsSigFunc = [](nuis::event::FullEvent const &) -> bool { return true; };
CompProjFunc =
[](nuis::event::FullEvent const &) -> std::array {
std::array arr;
for (NumericT &el : arr) {
el = 0;
}
return arr;
};
ProcessExtraFunc =
std::function();
fJournalReference = "";
fDOI = "";
fYear = "";
fTargetMaterial = "";
fFluxDescription = "";
fSignalDescription = "";
fIsShapeOnly = -1;
fIsFluxUnfolded = -1;
- energy_cut = nuis::utility::ENuRange{std::numeric_limits::max(),
- std::numeric_limits::max()};
+ energy_cut = nuis::utility::KinematicRange{
+ std::numeric_limits::max(), std::numeric_limits::max()};
}
+ std::string Name() { return fName; }
+
void SetUseCache(bool uc = true) {
fUseCache = uc;
if (!uc) {
fSignalCache.resize(0);
fProjectionCache.resize(0);
}
}
fhicl::ParameterSet fGlobalConfig;
fhicl::ParameterSet fInstanceConfig;
void ReadGlobalConfigDefaults() {
fGlobalConfig = nuis::config::GetDocument().get(
std::string("global.sample_configuration.") + Name(),
fhicl::ParameterSet());
if (!fJournalReference.length()) {
fJournalReference = fGlobalConfig.get("journal_reference",
"Not yet published");
}
if (!fDOI.length()) {
fDOI = fGlobalConfig.get("DOI", "Unknown DOI");
}
if (!fYear.length()) {
fYear = fGlobalConfig.get("year", "Unknown year");
}
if (!fTargetMaterial.length()) {
fTargetMaterial = fGlobalConfig.get(
"target_material", "Unknown target material");
}
if (!fFluxDescription.length()) {
fFluxDescription =
fGlobalConfig.get("flux_description", "Unknown flux");
}
if (!fSignalDescription.length()) {
fSignalDescription = fGlobalConfig.get("signal_description",
"Unknown Signal");
}
if (fIsShapeOnly == -1) {
fIsShapeOnly = fGlobalConfig.get("shape_only", false);
}
if (fIsFluxUnfolded == -1) {
fIsFluxUnfolded = fGlobalConfig.get("flux_unfolded", false);
}
if ((energy_cut.first == std::numeric_limits::max()) &&
(fGlobalConfig.has_key("enu_range"))) {
energy_cut = fGlobalConfig.get>("enu_range");
}
}
virtual std::string GetJournalReference() { return fJournalReference; }
virtual std::string GetDOI() { return fDOI; }
virtual std::string GetYear() { return fYear; }
virtual std::string GetTargetMaterial() { return fTargetMaterial; }
virtual std::string GetFluxDescription() { return fFluxDescription; }
virtual std::string GetSignalDescription() { return fSignalDescription; }
void SetShapeOnly(bool iso) { fIsShapeOnly = iso; }
void SetFluxUnfolded(bool ifo) { fIsFluxUnfolded = ifo; }
void SetData(std::string const &data_descriptor) {
fDataInputDescriptor = data_descriptor;
}
void SetMask(std::string const &mask_descriptor) {
fMaskInputDescriptor = mask_descriptor;
}
void SetCovariance(std::string const &cov_descriptor) {
fCovarianceInputDescriptor = cov_descriptor;
}
virtual void FillProjection(std::array const &proj,
NumericT event_weight) {
TH_Help::Fill(fPrediction, proj, event_weight);
}
+ virtual void FillModeProjection(std::array const &proj,
+ NumericT event_weight,
+ nuis::event::Channel_t mode) {
+ if (!fPrediction_modes.count(mode)) {
+ std::stringstream ss;
+ ss << "Prediction_" << mode;
+ fPrediction_modes[mode] =
+ nuis::utility::Clone(fPrediction, true, ss.str());
+ }
+ TH_Help::Fill(fPrediction_modes[mode], proj, event_weight);
+ }
+
+ virtual void ProcessSignalEvent(nuis::event::FullEvent const &fev,
+ double weight = 1) {
+ auto const &proj = CompProjFunc(fev);
+ FillProjection(proj, weight);
+ if (fModeHists) {
+ FillModeProjection(proj, weight, fev.mode);
+ }
+ }
virtual void FinalizeComparison() {
if (fComparisonFinalized) {
throw SimpleDataComparison_already_finalized()
<< "[ERROR]: Attempted to re-finalize a comparison for "
"SimpleDataComparison: "
<< std::quoted(Name());
}
fPrediction_xsec =
nuis::utility::Clone(fPrediction, false, "Prediction_xsec");
- IInputHandler const &IH =
- nuis::input::InputManager::Get().GetInputHandler(fIH_id);
-
TH_Help::Scale(fPrediction_xsec, 1.0, "width");
-
+ if (fModeHists) {
+ for (auto &mh : fPrediction_modes) {
+ TH_Help::Scale(mh.second, 1.0, "width");
+ }
+ }
// If we have a flux cut
if (energy_cut.first != std::numeric_limits::max()) {
+ IInputHandler const &IH =
+ nuis::input::InputManager::Get().GetInputHandler(fIH_id);
TH_Help::Scale(fPrediction_xsec, IH.GetXSecScaleFactor(energy_cut));
}
fPrediction_shape =
nuis::utility::Clone(fPrediction_xsec, false, "Prediction_shape");
if (fData) {
TH_Help::Scale(fPrediction_shape,
TH_Help::Integral(fData, "width") /
TH_Help::Integral(fPrediction_shape, "width"));
} else {
IEventProcessor_WARN(
"When Finalizing comparison, no Data histogram available.");
}
if (fIsFluxUnfolded) {
// fPrediction_comparison
} else if (fIsShapeOnly) {
fPrediction_comparison = nuis::utility::Clone(fPrediction_shape, false,
"Prediction_comparison");
} else {
fPrediction_comparison = nuis::utility::Clone(fPrediction_xsec, false,
"Prediction_comparison");
}
fComparisonFinalized = true;
}
void Initialize(fhicl::ParameterSet const &instance_sample_configuration) {
fInstanceConfig = instance_sample_configuration;
SetSampleVerbosity(fInstanceConfig.get(
"verbosity", nuis::config::GetDocument().get(
"global.sample.verbosity_default", "Reticent")));
ReadGlobalConfigDefaults();
if (fInstanceConfig.has_key("fake_data")) {
fData = nuis::utility::GetHistogram(
fInstanceConfig.get("fake_data_histogram"));
} else if (!fGlobalConfig.get("has_data", true) ||
!fInstanceConfig.get("has_data", true)) {
// Explicitly not expecting data
} else {
if (!fDataInputDescriptor.length()) {
if (!fGlobalConfig.has_key("data_descriptor")) {
throw invalid_SimpleDataComparison_initialization()
<< "[ERROR]: SimpleDataComparison::Initialize for "
"IDataComparison: "
<< std::quoted(Name())
<< " failed as no input data was set by a call to "
"SimpleDataComparison::SetData and no data_descriptor for "
"this SimpleDataComparison could be found in the global "
"configuration.";
}
fDataInputDescriptor =
fGlobalConfig.get("data_descriptor");
}
fData = nuis::utility::GetHistogram