diff --git a/cmake/NEUTSetup.cmake b/cmake/NEUTSetup.cmake
index 098c889..acf03aa 100644
--- a/cmake/NEUTSetup.cmake
+++ b/cmake/NEUTSetup.cmake
@@ -1,214 +1,217 @@
# Copyright 2016 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
-find_program(neut-config NEUTCONFIGFOUND)
+find_program(NEUTCONFIGFOUND NAMES neut-config)
LIST(APPEND EXTRA_CXX_FLAGS -DNEED_FILL_NEUT_COMMONS)
SET(HAVENEUTCONFIG FALSE)
# We are dealing with shiny NEUT
if(NOT "${NEUTCONFIGFOUND} " STREQUAL " ")
SET(HAVENEUTCONFIG TRUE)
+ cmessage(STATUS "Found neut-config, using it to determine configuration.")
+else()
+ cmessage(STATUS "Failed to find neut-config, assuming older NEUT build.")
endif()
if(HAVENEUTCONFIG)
execute_process (COMMAND neut-config
--version OUTPUT_VARIABLE NEUT_VERSION
OUTPUT_STRIP_TRAILING_WHITESPACE)
execute_process (COMMAND neut-config
--incdir OUTPUT_VARIABLE NEUT_INCLUDE_DIRS
OUTPUT_STRIP_TRAILING_WHITESPACE)
execute_process (COMMAND neut-config
--libdir OUTPUT_VARIABLE NEUT_LINK_DIRS
OUTPUT_STRIP_TRAILING_WHITESPACE)
GETLIBDIRS(neut-config --cernflags CERN_LIB_DIR)
LIST(APPEND NEUT_LINK_DIRS ${CERN_LIB_DIR})
GETLIBS(neut-config --cernflags CERN_LIBS)
if(USE_GENERATOR_REWEIGHT)
execute_process (COMMAND neut-config
--rwlibflags OUTPUT_VARIABLE NEUT_RWLIBS
OUTPUT_STRIP_TRAILING_WHITESPACE)
GETLIBS(neut-config --rwlibflags NEUT_RWLIBS)
LIST(APPEND NEUT_LIBS ${NEUT_RWLIBS})
else()
GETLIBS(neut-config --libflags NEUT_GENLIBS)
GETLIBS(neut-config --iolibflags NEUT_LIBS)
LIST(APPEND NEUT_LIBS ${NEUT_IOLIBS})
LIST(APPEND NEUT_LIBS ${NEUT_GENLIBS})
endif()
LIST(APPEND NEUT_LIBS ${CERN_LIBS};gfortran)
PrefixList(NEUT_LINK_DIRS "-L" ${NEUT_LINK_DIRS})
LIST(APPEND NEUT_CXX_FLAGS -DUSE_NEUT -DNEUT_VERSION=${NEUT_VERSION})
cmessage(STATUS "NEUT")
cmessage(STATUS " Version : ${NEUT_VERSION}")
cmessage(STATUS " Flags : ${NEUT_CXX_FLAGS}")
cmessage(STATUS " Includes : ${NEUT_INCLUDE_DIRS}")
cmessage(STATUS " Link Dirs : ${NEUT_LINK_DIRS}")
cmessage(STATUS " Libs : ${NEUT_LIBS}")
else() # Everything better be set up already
if(NEUT_ROOT STREQUAL "")
cmessage(FATAL_ERROR "Variable NEUT_ROOT is not defined. Please export environment variable NEUT_ROOT or configure with -DNEUT_ROOT=/path/to/NEUT. This must be set to point to a prebuilt NEUT instance.")
endif()
if(CERN STREQUAL "")
cmessage(FATAL_ERROR "Variable CERN is not defined. Please export environment variable CERN or configure with -DCERN=/path/to/CERNLIB. This must be set to point to a prebuilt CERNLIB instance.")
endif()
if(CERN_LEVEL STREQUAL "")
cmessage(FATAL_ERROR "Variable CERN_LEVEL is not defined. Please export environment variable CERN_LEVEL or configure with -DCERN_LEVEL=XXXX (likely to be 2005).")
endif()
if(${NEUT_VERSION} VERSION_LESS 5.4.0)
set(NEUT_LIB_DIR ${NEUT_ROOT}/lib/Linux_pc)
else()
set(NEUT_LIB_DIR ${NEUT_ROOT}/lib)
endif()
set(NEUT_CLASS ${NEUT_ROOT}/src/neutclass)
LIST(APPEND NEUT_INCLUDE_DIRS
${NEUT_ROOT}/include
${NEUT_ROOT}/src/neutclass)
LIST(APPEND NEUT_LINK_DIRS
${NEUT_LIB_DIR}
${CERN}/${CERN_LEVEL}/lib)
if(${NEUT_VERSION} VERSION_GREATER 5.4.1.999)
LIST(APPEND NEUT_LIBS
neutcore_5.4.2
nuccorspl_5.4.2 #typo in NEUT, may hopefully disappear
nuceff_5.4.2
partnuck_5.4.2
skmcsvc_5.4.2
tauola_5.4.2
HT2p2h_5.4.0
N1p1h_5.4.0)
LIST(APPEND NEUT_CXX_FLAGS -DNEUT_COMMON_QEAV)
elseif(${NEUT_VERSION} VERSION_EQUAL 5.4.1)
LIST(APPEND NEUT_LIBS
neutcore_5.4.1
nuccorspl_5.4.1 #typo in NEUT, may hopefully disappear
nuceff_5.4.1
partnuck_5.4.1
skmcsvc_5.4.1
tauola_5.4.1
HT2p2h_5.4.1
N1p1h_5.4.1)
elseif(${NEUT_VERSION} VERSION_GREATER 5.3.999)
LIST(APPEND NEUT_LIBS
neutcore_5.4.0
nuccorspl_5.4.0 #typo in NEUT, may hopefully disappear
nuceff_5.4.0
partnuck_5.4.0
skmcsvc_5.4.0
tauola_5.4.0
HT2p2h_5.4.0
N1p1h_5.4.0)
else()
LIST(APPEND NEUT_LIBS
neutcore
nuccorrspl
nuceff
partnuck
skmcsvc
tauola)
endif()
LIST(APPEND NEUT_LIBS
jetset74
pdflib804
mathlib
packlib
pawlib
gfortran)
if(USE_GENERATOR_REWEIGHT)
LIST(APPEND NEUT_INCLUDE_DIRS
${NEUT_ROOT}/src/reweight)
LIST(APPEND NEUT_LINK_DIRS
${NEUT_ROOT}/src/reweight)
LIST(REVERSE NEUT_LIBS)
LIST(APPEND NEUT_LIBS NReWeight)
LIST(REVERSE NEUT_LIBS)
endif()
set(NEUT_ROOT_SHAREDOBJS)
LIST(APPEND NEUT_ROOT_SHAREDOBJS
${NEUT_CLASS}/neutctrl.so
${NEUT_CLASS}/neutfsivert.so)
# Check for new versions of NEUT with NUCLEON FSI
if(EXISTS "${NEUT_CLASS}/neutnucfsistep.so")
set(NEUT_NUCFSI 1)
LIST(APPEND NEUT_CXX_FLAGS -DNEUT_NUCFSI_ENABLED)
LIST(APPEND NEUT_ROOT_SHAREDOBJS
${NEUT_CLASS}/neutnucfsistep.so
${NEUT_CLASS}/neutnucfsivert.so
)
endif()
if(${NEUT_VERSION} VERSION_LESS 5.4.0)
LIST(APPEND NEUT_ROOT_SHAREDOBJS
${NEUT_CLASS}/neutrootTreeSingleton.so)
endif()
LIST(APPEND NEUT_ROOT_SHAREDOBJS
${NEUT_CLASS}/neutvtx.so
${NEUT_CLASS}/neutfsipart.so
${NEUT_CLASS}/neutpart.so
${NEUT_CLASS}/neutvect.so
)
foreach(so ${NEUT_ROOT_SHAREDOBJS})
get_filename_component(SONAME ${so} NAME_WE)
add_library(${SONAME} SHARED IMPORTED)
set_property(TARGET ${SONAME} PROPERTY IMPORTED_LOCATION ${so})
LIST(APPEND NEUT_IMPORTED_TARGETS ${SONAME})
endforeach()
PrefixList(NEUT_LINK_DIRS "-L" ${NEUT_LINK_DIRS})
LIST(APPEND NEUT_CXX_FLAGS -DUSE_NEUT -DNEUT_VERSION=${NEUT_VERSION})
cmessage(STATUS "NEUT")
cmessage(STATUS " Version : ${NEUT_VERSION}")
cmessage(STATUS " Flags : ${NEUT_CXX_FLAGS}")
cmessage(STATUS " Includes : ${NEUT_INCLUDE_DIRS}")
cmessage(STATUS " Link Dirs : ${NEUT_LINK_DIRS}")
cmessage(STATUS " Libs : ${NEUT_LIBS}")
cmessage(STATUS " SOs : ${NEUT_IMPORTED_TARGETS}")
endif()
SET(USE_NEUT TRUE CACHE BOOL "Whether to enable NEUT support. Requires external libraries. " FORCE)
diff --git a/cmake/c++CompilerSetup.cmake b/cmake/c++CompilerSetup.cmake
index 5083311..20ffe38 100644
--- a/cmake/c++CompilerSetup.cmake
+++ b/cmake/c++CompilerSetup.cmake
@@ -1,57 +1,58 @@
# Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
################################################################################
# This file is part of NUISANCE.
#
# NUISANCE is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# NUISANCE is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with NUISANCE. If not, see .
################################################################################
set(CXX_WARNINGS -Wall -Wextra)
+set(CXX_IGNORE_WARNINGS -Wno-delete-non-virtual-dtor -Wno-unused -Wno-misleading-indentation)
-LIST(APPEND EXTRA_CXX_FLAGS ${CXX_WARNINGS} -Werror -Wno-delete-non-virtual-dtor -Wno-unused "-D__FILENAME__=\"$(subst ${CMAKE_SOURCE_DIR}/,,$(abspath $<))\"")
+LIST(APPEND EXTRA_CXX_FLAGS ${CXX_WARNINGS} -Werror ${CXX_IGNORE_WARNINGS} "-D__FILENAME__=\"$(subst ${CMAKE_SOURCE_DIR}/,,$(abspath $<))\"")
BuildFlagString(NUISANCE_LINK_DIRS "-L" ${EXTRA_LINK_DIRS})
LIST(APPEND EXTRA_LIBS dl)
BuildLibraryFlagString(STR_EXTRA_LIBS ${EXTRA_LIBS})
BuildFlagString(STR_EXTRA_SHAREDOBJS " " ${EXTRA_SHAREDOBJS})
#This ends up holding all of the libraries and search paths for extenal dependencies
CatStringsIfNotEmpty(NUISANCE_DEPEND_LIBS
${STR_EXTRA_SHAREDOBJS}
${STR_EXTRA_LIBS})
BuildFlagString(STR_EXTRA_LINK_FLAGS " " ${EXTRA_LINK_FLAGS})
CatStringsIfNotEmpty(CMAKE_LINK_FLAGS
${CMAKE_LINK_FLAGS}
${STR_EXTRA_LINK_FLAGS})
get_directory_property(NUISANCE_INCLUDE_DIRS INCLUDE_DIRECTORIES)
BuildFlagString(NUISANCE_CXX_FLAGS " " ${EXTRA_CXX_FLAGS})
CatStringsIfNotEmpty(NUISANCE_CXX_FLAGS ${CMAKE_CXX_FLAGS} ${NUISANCE_CXX_FLAGS} )
if (VERBOSE)
cmessage (STATUS "C++ Compiler : ${CXX_COMPILER_NAME}")
cmessage (STATUS " Flags : ${NUISANCE_CXX_FLAGS}")
cmessage (STATUS " Release Flags : ${CMAKE_CXX_FLAGS_RELEASE}")
cmessage (STATUS " Debug Flags : ${CMAKE_CXX_FLAGS_DEBUG}")
cmessage (STATUS " Include Dirs : ${NUISANCE_INCLUDE_DIRS}")
cmessage (STATUS " Linker Flags : ${CMAKE_LINK_FLAGS}")
cmessage (STATUS " Link Dirs : ${NUISANCE_LINK_DIRS}")
cmessage (STATUS " Lib Flags : ${NUISANCE_DEPEND_LIBS}")
endif()
add_compile_options(${EXTRA_CXX_FLAGS})
diff --git a/src/app/nuisevsum.cxx b/src/app/nuisevsum.cxx
index 5461e22..c50f41b 100644
--- a/src/app/nuisevsum.cxx
+++ b/src/app/nuisevsum.cxx
@@ -1,79 +1,88 @@
#include "config/GlobalConfiguration.hxx"
#include "input/IInputHandler.hxx"
#include "event/MinimalEvent.hxx"
#include "samples/IEventProcessor.hxx"
#include "plugins/Instantiate.hxx"
#include "exception/exception.hxx"
#include "fhiclcpp/make_ParameterSet.h"
#include "string_parsers/from_string.hxx"
#include
NEW_NUIS_EXCEPT(invalid_cli_arguments);
size_t NMax = std::numeric_limits::max();
std::string input_file;
std::string input_type;
+bool show_intermediates = false;
void SayUsage(char const *argv[]) {
std::cout << "[USAGE]: " << argv[0]
<< "\n"
"\t-i : Input file passed to named "
"IInputHandler instance \n"
"\t-H : Name of IInputHandler subclass "
"capable of reading NUISANCE events from the argument of -i.\n"
"\t-n : Maximum number of events to "
"read. Will read entire input file by default.\n"
+ "\t--show-intermediates : Show pre-FSI particles where"
+ " relevant.\n"
<< std::endl;
}
void handleOpts(int argc, char const *argv[]) {
int opt = 1;
while (opt < argc) {
if ((std::string(argv[opt]) == "-?") ||
(std::string(argv[opt]) == "--help")) {
SayUsage(argv);
exit(0);
} else if (std::string(argv[opt]) == "-i") {
input_file = argv[++opt];
} else if (std::string(argv[opt]) == "-H") {
input_type = argv[++opt];
} else if (std::string(argv[opt]) == "-n") {
NMax = fhicl::string_parsers::str2T(argv[++opt]);
+ } else if (std::string(argv[opt]) == "--show-intermediates") {
+ show_intermediates = true;
} else {
std::cout << "[ERROR]: Unknown option: " << argv[opt] << std::endl;
SayUsage(argv);
exit(1);
}
opt++;
}
}
int main(int argc, char const *argv[]) {
nuis::config::EnsureConfigurationRead("nuis.global.config.fcl");
handleOpts(argc, argv);
if (!input_type.length() || !input_file.length()) {
SayUsage(argv);
throw invalid_cli_arguments()
<< "[ERROR]: Require both -i and -H cli options to be passed.";
}
fhicl::ParameterSet sample_config;
sample_config.put("input_type", input_type);
sample_config.put("file", input_file);
+ if (show_intermediates) {
+ sample_config.put("keep_intermediates", true);
+ }
- nuis::plugins::plugin_traits::unique_ptr_t VerboseEventSummary =
- nuis::plugins::Instantiate("VerboseEventSummary");
+ nuis::plugins::plugin_traits::unique_ptr_t
+ VerboseEventSummary =
+ nuis::plugins::Instantiate("VerboseEventSummary");
VerboseEventSummary->Initialize(sample_config);
VerboseEventSummary->ProcessSample(NMax);
}
diff --git a/src/generator/input/GENIEInputHandler.cxx b/src/generator/input/GENIEInputHandler.cxx
index 1bfb67d..f820fd9 100644
--- a/src/generator/input/GENIEInputHandler.cxx
+++ b/src/generator/input/GENIEInputHandler.cxx
@@ -1,215 +1,222 @@
#include "generator/input/GENIEInputHandler.hxx"
#include "generator/utility/GENIEUtility.hxx"
#include "utility/PDGCodeUtility.hxx"
#ifdef GENIE_V3_INTERFACE
+#include "Framework/Conventions/Units.h"
#include "Framework/EventGen/EventRecord.h"
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/GHEP/GHepRecord.h"
#else
+#include "Conventions/Units.h"
#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), fFileWeight(1) {}
GENIEInputHandler::GENIEInputHandler(GENIEInputHandler &&other)
: fInputTree(std::move(other.fInputTree)),
fReaderEvent(std::move(other.fReaderEvent)), fGenieNtpl(nullptr),
fFileWeight(1) {}
std::set GENIEInputHandler::GetSplineList() {
std::set splinenames;
- size_t NEvents = 1000; // GetNEvents();
+ size_t NEvents = std::min(size_t(10000), GetNEvents());
size_t ShoutEvery = NEvents / 100;
+ std::cout << "[INFO]: Determining list of required GENIE splines."
+ << std::endl;
std::cout << "[INFO]: Read " << 0 << "/" << NEvents << " GENIE events."
<< std::flush;
for (size_t ev_it = 0; ev_it < NEvents; ++ev_it) {
if (fGenieNtpl) {
fGenieNtpl->Clear();
}
fInputTree.tree->GetEntry(ev_it);
if (ShoutEvery && !(ev_it % ShoutEvery)) {
std::cout << "\r[INFO]: Read " << ev_it << "/" << NEvents
<< " GENIE events." << std::flush;
}
genie::GHepRecord *GHep =
static_cast(fGenieNtpl->event);
if (!GHep) {
throw invalid_GENIE_event() << "[ERROR]: GENIE event " << ev_it
<< " failed to contain a GHepRecord";
}
splinenames.insert(GHep->Summary()->AsString());
}
std::cout << "[INFO]: Having read GENIE events, found " << splinenames.size()
<< " splines: " << std::endl;
for (auto const &sn : splinenames) {
std::cout << "\t" << sn << std::endl;
}
return splinenames;
}
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") ||
- XSecInfo.has_key("filter_splines_by_event_tree")) {
+ } else if (XSecInfo.has_key("spline_file")) {
std::cout
<< "[INFO]: Attempting to build GENIE file weight with input splines."
<< std::endl;
if (XSecInfo.get("filter_splines_by_event_tree", false)) {
std::cout << "[INFO]: Ensuring we only read relevant splines... this "
"will take a while."
<< std::endl;
fFileWeight = genietools::GetFileWeight(XSecInfo, GetSplineList()) /
double(GetNEvents());
} else {
fFileWeight = genietools::GetFileWeight(XSecInfo) / double(GetNEvents());
}
std::cout << "[INFO]: Average GENIE XSecWeight = " << fFileWeight << ""
<< std::endl;
+ } else if (XSecInfo.has_key("flux")) {
+ fFileWeight = genietools::GetFileWeightFromReconstructedSplines(
+ XSecInfo, fGenieNtpl, fInputTree.tree);
+ } else {
+ fFileWeight = 1;
}
}
void GENIEInputHandler::GetEntry(ev_index_t idx) const {
if (idx >= GetNEvents()) {
throw IInputHandler::invalid_entry()
<< "[ERROR]: Attempted to get entry " << idx
<< " from an InputHandler with only " << GetNEvents();
}
if (fGenieNtpl) {
fGenieNtpl->Clear();
}
fInputTree.tree->GetEntry(idx);
}
MinimalEvent const &GENIEInputHandler::GetMinimalEvent(ev_index_t idx) const {
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 = 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 && (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/utility/GENIEUtility.cxx b/src/generator/utility/GENIEUtility.cxx
index 6011309..741d044 100644
--- a/src/generator/utility/GENIEUtility.cxx
+++ b/src/generator/utility/GENIEUtility.cxx
@@ -1,479 +1,588 @@
// 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/InteractionChannelUtility.hxx"
#include "utility/PDGCodeUtility.hxx"
#include "utility/ROOTUtility.hxx"
#include "utility/StringUtility.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 {
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) {
/*
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 = 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,
std::set const &spline_list) {
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;
+ std::set FoundSplines;
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;
+ << " and multiplied by " << molwght << 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 (size_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;
+ double weight = TargetSplines[ts_it].MolecularWeight;
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;
}
}
+ FoundSplines.insert(all_splines[ts_it][c_it].GetName());
for (size_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
+ // Copied from GENIE/Conventions/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);
}
}
}
+ for (auto const &sn : spline_list) {
+ std::cout << "Want: " << sn << std::endl;
+ if (FoundSplines.count(sn)) {
+ std::cout << "\tFound!" << std::endl;
+ } else {
+ std::cout << "Didn't find." << std::endl;
+ }
+ }
+
// Sum all the correctly weighted per-target splines
TGraph MasterSpline(NKnots);
+ TGraph MasterSpline_pernuc(NKnots);
for (size_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 (size_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;
-
+ XSec += TargetSplines[c_it].TotSpline.Eval(E);
MasterSpline.SetPoint(p_it, E, XSec);
}
}
+ for (size_t p_it = 0; p_it < NKnots; ++p_it) {
+ double E, XSec;
+ MasterSpline.GetPoint(p_it, E, XSec);
+ MasterSpline_pernuc.SetPoint(p_it, E, XSec * WeightToPerNucleon);
+ }
+
nuis::utility::TFile_ptr outfile(nullptr, [](TFile *) {});
if (xsecinfo.has_key("spline_output_file")) {
outfile = nuis::utility::CheckOpenTFile(
xsecinfo.get("spline_output_file"), "RECREATE");
for (auto &sp : TargetSplines) {
std::string spname = nuis::utility::SanitizeROOTObjectName(
sp.search_term.substr(2, sp.search_term.length() - 4));
std::cout << "[INFO]: Dumping spline: " << spname << std::endl;
nuis::utility::WriteToTFile(outfile, &sp.TotSpline, spname.c_str());
}
nuis::utility::WriteToTFile(outfile, &MasterSpline, "TotalXSec");
+ nuis::utility::WriteToTFile(outfile, &MasterSpline_pernuc,
+ "TotalXSec_pernuc");
}
fhicl::ParameterSet fluxps = xsecinfo.get("flux");
// If MonoE
if (fluxps.has_key("energy_GeV")) {
double monoE = fluxps.get("energy_GeV");
- return MasterSpline.Eval(monoE);
+ return MasterSpline_pernuc.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"));
if (outfile) {
std::unique_ptr Fluxcp = nuis::utility::Clone(Flux);
nuis::utility::WriteToTFile(outfile, Fluxcp.get(), "Flux");
}
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);
+ if ((e < EMin) || (e > EMax)) {
+ continue;
+ }
+ avg_xsec += MasterSpline_pernuc.Eval(e);
}
avg_xsec /= double(NIntSteps);
Flux->SetBinContent(bi_it + 1, bc * avg_xsec);
}
if (outfile) {
std::unique_ptr Fluxcp = nuis::utility::Clone(Flux);
nuis::utility::WriteToTFile(outfile, Fluxcp.get(), "EvRate");
}
double EvRateIntegral = Flux->Integral(per_width ? "width" : "");
return EvRateIntegral / FluxIntegral;
}
}
+struct TargetXSecHistBlob {
+ double MolecularWeight;
+ size_t NNucleons;
+ std::string search_term;
+ TH1 *TotHist;
+};
+
+std::set GENIEInputHandler::GetFileWeightFromReconstructedSplines(
+ fhicl::ParameterSet const &xsecinfo, TTree *GENIETree,
+ genie::NtpMCEventRecord *GenieNtpl) {
+
+
+
+ double EMin = xsecinfo.get("EMin", 0);
+ double EMax = xsecinfo.get("EMax", 10);
+ std::vector TargetXSecHists;
+ 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;
+
+ TargetXSecHists.push_back(
+ {molwght, (nuc_pdg % 1000) / 10, ss.str(), TGraph()});
+
+ std::cout << "[INFO]:\tGetting splines that match: \"" << ss.str()
+ << "\" with A = " << TargetXSecHists.back().NNucleons
+ << " and multiplied by " << molwght << std::endl;
+ }
+
+ size_t NTotNucleons = 0;
+ std::vector spline_patterns;
+ for (auto const &spb : TargetXSecHists) {
+ spline_patterns.push_back(spb.search_term);
+ NTotNucleons += spb.MolecularWeight * spb.NNucleons;
+ }
+ double WeightToPerNucleon = 1.0 / double(NTotNucleons);
+
+ std::map splines;
+ size_t NEvents = GetNEvents();
+ size_t ShoutEvery = NEvents / 100;
+ std::cout << "[INFO]: Rebuilding GENIE splines." << std::endl;
+ std::cout << "[INFO]: Read " << 0 << "/" << NEvents << " GENIE events."
+ << std::flush;
+ for (size_t ev_it = 0; ev_it < NEvents; ++ev_it) {
+ if (fGenieNtpl) {
+ fGenieNtpl->Clear();
+ }
+ fInputTree.tree->GetEntry(ev_it);
+
+ if (ShoutEvery && !(ev_it % ShoutEvery)) {
+ std::cout << "\r[INFO]: Read " << ev_it << "/" << NEvents
+ << " GENIE events." << std::flush;
+ }
+
+ genie::GHepRecord *GHep =
+ static_cast(fGenieNtpl->event);
+ if (!GHep) {
+ throw invalid_GENIE_event() << "[ERROR]: GENIE event " << ev_it
+ << " failed to contain a GHepRecord";
+ }
+ std::string const &spline_name = GHep->Summary()->AsString();
+ if (!splinenames.count(spline_name)) {
+ splinenames.emplace(spline_name, TGraph());
+ }
+ double E_GeV = GHep->Probe()->E();
+ double XSec = (genie_record.XSec() / (1E-38 * genie::units::cm2));
+ splinenames[spline_name].SetPoint(E_GeV, XSec);
+ }
+ std::cout << "[INFO]: Having read GENIE events, found " << splinenames.size()
+ << " splines: " << std::endl;
+ for (auto const &sn : splinenames) {
+ std::cout << "\t" << sn << std::endl;
+ }
+ return splinenames;
+}
+
} // namespace genietools
} // namespace nuis
diff --git a/src/generator/utility/NEUTUtility.cxx b/src/generator/utility/NEUTUtility.cxx
index 0b0e04a..d1cf057 100644
--- a/src/generator/utility/NEUTUtility.cxx
+++ b/src/generator/utility/NEUTUtility.cxx
@@ -1,95 +1,97 @@
#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);
+// #define DEBUG_NEUT_UTILITY
+
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) {
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
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/samples/nuA/Nuclear/MINERvA/CC0Pi/MINERvA_CC0Pi_CH_xsec_ptpl_nu.cxx b/src/samples/nuA/Nuclear/MINERvA/CC0Pi/MINERvA_CC0Pi_CH_xsec_ptpl_nu.cxx
new file mode 100644
index 0000000..e69de29
diff --git a/src/samples/nuA/Nuclear/MINERvA/CCInclusive/CMakeLists.txt b/src/samples/nuA/Nuclear/MINERvA/CCInclusive/CMakeLists.txt
new file mode 100644
index 0000000..01227b9
--- /dev/null
+++ b/src/samples/nuA/Nuclear/MINERvA/CCInclusive/CMakeLists.txt
@@ -0,0 +1,17 @@
+SET(SAMPLES
+ MINERvA_CCInclusive_CH_xsec_EAvailq3_nu)
+
+foreach(S ${SAMPLES})
+ LIST(APPEND INuADataComparisons ${S})
+ LIST(APPEND SAMPLES_src ${S}.cxx)
+endforeach()
+
+add_library(MINERvADataComparisons SHARED ${SAMPLES_src})
+target_link_libraries(MINERvADataComparisons
+ nuis_event nuis_input nuis_config nuis_persistency)
+
+install(TARGETS MINERvADataComparisons DESTINATION plugins)
+
+install(FILES MINERvA.CCInclusive.sample.config.fcl DESTINATION fcl)
+
+SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE)
diff --git a/src/samples/nuA/Nuclear/MINERvA/CCInclusive/MINERvA.CCInclusive.sample.config.fcl b/src/samples/nuA/Nuclear/MINERvA/CCInclusive/MINERvA.CCInclusive.sample.config.fcl
new file mode 100644
index 0000000..e69de29
diff --git a/src/samples/nuA/Nuclear/MINERvA/CCInclusive/MINERvA_CCInclusive_CH_xsec_EAvailq3_nu.cxx b/src/samples/nuA/Nuclear/MINERvA/CCInclusive/MINERvA_CCInclusive_CH_xsec_EAvailq3_nu.cxx
new file mode 100644
index 0000000..e69de29
diff --git a/src/samples/nuA/Nuclear/MINERvA/MINERvA.sample.config.fcl b/src/samples/nuA/Nuclear/MINERvA/MINERvA.sample.config.fcl
index 09200c4..afeeb92 100644
--- a/src/samples/nuA/Nuclear/MINERvA/MINERvA.sample.config.fcl
+++ b/src/samples/nuA/Nuclear/MINERvA/MINERvA.sample.config.fcl
@@ -1 +1,2 @@
#include "MINERvA.CC0Pi.sample.config.fcl"
+# include "MINERvA.CCInclusive.sample.config.fcl"
diff --git a/src/samples/nuA/Nuclear/MicroBooNE/CCInclusive/CMakeLists.txt b/src/samples/nuA/Nuclear/MicroBooNE/CCInclusive/CMakeLists.txt
new file mode 100644
index 0000000..dd82b30
--- /dev/null
+++ b/src/samples/nuA/Nuclear/MicroBooNE/CCInclusive/CMakeLists.txt
@@ -0,0 +1,17 @@
+SET(SAMPLES
+MicroBooNE_CCInclusive_Ar40_xsec_reco_pmu_thetamu_nu)
+
+foreach(S ${SAMPLES})
+ LIST(APPEND INuADataComparisons ${S})
+ LIST(APPEND SAMPLES_src ${S}.cxx)
+endforeach()
+
+add_library(MicroBooNEDataComparisons SHARED ${SAMPLES_src})
+target_link_libraries(MicroBooNEDataComparisons
+ nuis_event nuis_input nuis_config nuis_persistency)
+
+install(TARGETS MicroBooNEDataComparisons DESTINATION plugins)
+
+install(FILES MicroBooNE.CCInclusive.sample.config.fcl DESTINATION fcl)
+
+SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE)
diff --git a/src/samples/nuA/Nuclear/MicroBooNE/CCInclusive/MicroBooNE.CCInclusive.sample.config.fcl b/src/samples/nuA/Nuclear/MicroBooNE/CCInclusive/MicroBooNE.CCInclusive.sample.config.fcl
new file mode 100644
index 0000000..e69de29
diff --git a/src/samples/nuA/Nuclear/MicroBooNE/CCInclusive/MicroBooNE_CCInclusive_Ar40_xsec_reco_pmu_thetamu_nu.cxx b/src/samples/nuA/Nuclear/MicroBooNE/CCInclusive/MicroBooNE_CCInclusive_Ar40_xsec_reco_pmu_thetamu_nu.cxx
new file mode 100644
index 0000000..e69de29
diff --git a/src/samples/nuA/Nuclear/MicroBooNE/CMakeLists.txt b/src/samples/nuA/Nuclear/MicroBooNE/CMakeLists.txt
new file mode 100644
index 0000000..01ffc5a
--- /dev/null
+++ b/src/samples/nuA/Nuclear/MicroBooNE/CMakeLists.txt
@@ -0,0 +1,7 @@
+add_subdirectory(CCInclusive)
+
+LIST(APPEND INuADataComparisons_FHiCL MicroBooNE.sample.config.fcl)
+install(FILES MicroBooNE.sample.config.fcl DESTINATION fcl)
+
+SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE)
+SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE)
diff --git a/src/samples/nuA/Nuclear/MicroBooNE/MicroBooNE.sample.config.fcl b/src/samples/nuA/Nuclear/MicroBooNE/MicroBooNE.sample.config.fcl
new file mode 100644
index 0000000..675ccaf
--- /dev/null
+++ b/src/samples/nuA/Nuclear/MicroBooNE/MicroBooNE.sample.config.fcl
@@ -0,0 +1 @@
+#include "MiniBooNE.CC0Pi.sample.config.fcl"
diff --git a/src/samples/nuA/Nuclear/MiniBooNE/CC0Pi/CMakeLists.txt b/src/samples/nuA/Nuclear/MiniBooNE/CC0Pi/CMakeLists.txt
new file mode 100644
index 0000000..d2e924d
--- /dev/null
+++ b/src/samples/nuA/Nuclear/MiniBooNE/CC0Pi/CMakeLists.txt
@@ -0,0 +1,17 @@
+SET(SAMPLES
+ MiniBooNE_CCQE_CH2_xsec_nu)
+
+foreach(S ${SAMPLES})
+ LIST(APPEND INuADataComparisons ${S})
+ LIST(APPEND SAMPLES_src ${S}.cxx)
+endforeach()
+
+add_library(MiniBooNEDataComparisons SHARED ${SAMPLES_src})
+target_link_libraries(MiniBooNEDataComparisons
+ nuis_event nuis_input nuis_config nuis_persistency)
+
+install(TARGETS MiniBooNEDataComparisons DESTINATION plugins)
+
+install(FILES MiniBooNE.CC0Pi.sample.config.fcl DESTINATION fcl)
+
+SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE)
diff --git a/src/samples/nuA/Nuclear/MiniBooNE/CC0Pi/MiniBooNE.CC0Pi.sample.config.fcl b/src/samples/nuA/Nuclear/MiniBooNE/CC0Pi/MiniBooNE.CC0Pi.sample.config.fcl
new file mode 100644
index 0000000..e69de29
diff --git a/src/samples/nuA/Nuclear/MiniBooNE/CC0Pi/MiniBooNE_CCQE_CH2_xsec_nu.cxx b/src/samples/nuA/Nuclear/MiniBooNE/CC0Pi/MiniBooNE_CCQE_CH2_xsec_nu.cxx
new file mode 100644
index 0000000..e69de29
diff --git a/src/samples/nuA/Nuclear/MiniBooNE/CMakeLists.txt b/src/samples/nuA/Nuclear/MiniBooNE/CMakeLists.txt
new file mode 100644
index 0000000..f226c0e
--- /dev/null
+++ b/src/samples/nuA/Nuclear/MiniBooNE/CMakeLists.txt
@@ -0,0 +1,7 @@
+add_subdirectory(CC0Pi)
+
+LIST(APPEND INuADataComparisons_FHiCL MiniBooNE.sample.config.fcl)
+install(FILES MiniBooNE.sample.config.fcl DESTINATION fcl)
+
+SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE)
+SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE)
diff --git a/src/samples/nuA/Nuclear/MiniBooNE/MiniBooNE.sample.config.fcl b/src/samples/nuA/Nuclear/MiniBooNE/MiniBooNE.sample.config.fcl
new file mode 100644
index 0000000..675ccaf
--- /dev/null
+++ b/src/samples/nuA/Nuclear/MiniBooNE/MiniBooNE.sample.config.fcl
@@ -0,0 +1 @@
+#include "MiniBooNE.CC0Pi.sample.config.fcl"
diff --git a/src/variation/engines/T2KReWeight.cxx b/src/variation/engines/T2KReWeight.cxx
index 71eda4e..e9d6d77 100644
--- a/src/variation/engines/T2KReWeight.cxx
+++ b/src/variation/engines/T2KReWeight.cxx
@@ -1,102 +1,103 @@
#include "variation/engines/T2KReWeight.hxx"
#include "fhiclcpp/ParameterSet.h"
#include "event/MinimalEvent.hxx"
// T2K Engine includes
#include "T2KNeutReWeight.h"
#include "T2KNIWGReWeight.h"
#include "T2KReWeight.h"
#include
using namespace nuis;
using namespace nuis::params;
using namespace nuis::event;
void T2KReWeightEngine::Initialize(fhicl::ParameterSet const &ps) {
fT2KRW = std::make_unique();
TDirectory *dir = gDirectory;
fT2KRW->AdoptWghtEngine("NEUT", new t2krew::T2KNeutReWeight);
fT2KRW->AdoptWghtEngine("NIWG", new t2krew::T2KNIWGReWeight);
dir->cd();
for (fhicl::ParameterSet param_ps :
ps.get>("parameters")) {
param_ps.put("type", GetName());
std::string const ¶m_name = param_ps.get("name");
T2KSystParam tsp;
tsp.t2ksyst = t2krew::T2KSyst::FromString(param_name);
if (tsp.t2ksyst == t2krew::kSystNull) {
throw invalid_T2K_syst_name()
<< "[ERROR]: T2KReWeight failed to parse " << std::quoted(param_name)
<< " as a T2K dial.";
}
fT2KRW->Systematics().Include(tsp.t2ksyst);
tsp.pid = ParameterManager::Get().EnsureParameterRegistered(param_ps);
fT2KSysts.push_back(tsp);
}
Reconfigure();
}
void T2KReWeightEngine::Reconfigure() {
for (T2KSystParam const &tsp : fT2KSysts) {
double val = ParameterManager::Get().GetParameterValue(tsp.pid);
fT2KRW->Systematics().SetTwkDial(tsp.t2ksyst, val);
+ std::cout << "Reconf: " << t2krew::T2KSyst::AsString(tsp.t2ksyst) << " -> " << val << std::endl;
}
fT2KRW->Reconfigure();
}
double T2KReWeightEngine::GetEventWeight(nuis::event::MinimalEvent const &ev) {
NeutVect const *nv = static_cast(ev.fGenEvent);
if (!nv) {
return 1.0;
}
NeutVect *nv_nc = const_cast(nv);
return fT2KRW->CalcWeight(nv_nc);
}
std::string T2KReWeightEngine::GetName() { return "T2KReWeightEngine"; }
std::string T2KReWeightEngine::GetDocumentation() { return ""; }
fhicl::ParameterSet T2KReWeightEngine::GetExampleConfiguration() {
fhicl::ParameterSet ps;
ps.put("name", GetName());
fhicl::ParameterSet dial_maqe;
dial_maqe.put("name", "NXSec_MaCCQE");
dial_maqe.put("start", 0);
dial_maqe.put("min", -3);
dial_maqe.put("max", 3);
dial_maqe.put("step", 0.1);
fhicl::ParameterSet dial_mares;
dial_mares.put("name", "NXSec_MARES");
dial_mares.put("start", 0);
dial_mares.put("min", -3);
dial_mares.put("max", 3);
dial_mares.put("step", 0.1);
ps.put>(
"parameters", std::vector{{dial_maqe, dial_mares}});
return ps;
}
GeneratorManager::Generator_id_t T2KReWeightEngine::GetGeneratorId() {
return GeneratorManager::Get().EnsureGeneratorRegistered("NEUT");
}
DECLARE_PLUGIN(IWeightProvider, T2KReWeightEngine);