diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2b088f9..6b75c51 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,78 +1,79 @@
# 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 <http://www.gnu.org/licenses/>.
################################################################################
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)
################################## COMPILER ####################################
+include(${CMAKE_SOURCE_DIR}/cmake/gperfSetup.cmake)
include(${CMAKE_SOURCE_DIR}/cmake/c++CompilerSetup.cmake)
################################################################################
add_subdirectory(src)
add_subdirectory(scripts)
add_subdirectory(config)
add_subdirectory(data)
diff --git a/cmake/gperfSetup.cmake b/cmake/gperfSetup.cmake
new file mode 100644
index 0000000..3d54519
--- /dev/null
+++ b/cmake/gperfSetup.cmake
@@ -0,0 +1,60 @@
+# 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 <http://www.gnu.org/licenses/>.
+################################################################################
+
+if(DEFINED USE_GPERFTOOLS AND USE_GPERFTOOLS)
+
+ ExternalProject_Add(libunwind
+ PREFIX "${CMAKE_BINARY_DIR}/Ext"
+ GIT_REPOSITORY "https://github.com/libunwind/libunwind.git"
+ GIT_TAG v1.3.1
+ CONFIGURE_COMMAND ${CMAKE_BINARY_DIR}/Ext/src/libunwind/autogen.sh --prefix=${CMAKE_INSTALL_PREFIX}
+ UPDATE_DISCONNECTED 1
+ UPDATE_COMMAND ""
+ BUILD_COMMAND make
+ INSTALL_COMMAND make install
+ )
+
+ ExternalProject_Add(gperftools
+ PREFIX "${CMAKE_BINARY_DIR}/Ext"
+ GIT_REPOSITORY "https://github.com/gperftools/gperftools.git"
+ GIT_TAG "gperftools-2.7"
+ CONFIGURE_COMMAND ./autogen.sh && ./configure --prefix=${CMAKE_INSTALL_PREFIX} CPPFLAGS=-I${CMAKE_INSTALL_PREFIX}/include LDFLAGS=-L${CMAKE_INSTALL_PREFIX}/lib
+ BUILD_IN_SOURCE 1
+ UPDATE_DISCONNECTED 1
+ UPDATE_COMMAND ""
+ BUILD_COMMAND make
+ INSTALL_COMMAND make install
+ )
+
+ add_dependencies(gperftools libunwind)
+
+ LIST(APPEND EXTRA_CXX_FLAGS
+ -fno-builtin-malloc
+ -fno-builtin-calloc
+ -fno-builtin-realloc
+ -fno-builtin-free)
+
+
+ LIST(APPEND EXTRA_LINK_DIRS ${CMAKE_INSTALL_PREFIX}/lib)
+
+ ##Want to prepend them
+ LIST(APPEND EXTRA_LIBS tcmalloc_and_profiler)
+
+ cmessage(STATUS "Using google performance libraries")
+endif()
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 63767ba..68b52bd 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -1,18 +1,19 @@
add_subdirectory(config)
add_subdirectory(event)
add_subdirectory(exception)
add_subdirectory(generator)
add_subdirectory(input)
add_subdirectory(parameters)
add_subdirectory(persistency)
add_subdirectory(plugins)
add_subdirectory(samples)
add_subdirectory(utility)
add_subdirectory(variation)
SET(INuADataComparisons_List ${INuADataComparisons_List} PARENT_SCOPE)
SET(INuADataComparisons_FHiCL ${INuADataComparisons_FHiCL} PARENT_SCOPE)
SET(GENERATOR_DEPENDENT_TARGETS ${GENERATOR_DEPENDENT_TARGETS} PARENT_SCOPE)
+add_subdirectory(tests)
add_subdirectory(app)
diff --git a/src/app/CMakeLists.txt b/src/app/CMakeLists.txt
index c84c5fe..34feffb 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_experimental 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_generator_utility 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/generator/input/GENIEInputHandler.cxx b/src/generator/input/GENIEInputHandler.cxx
index cc5c146..1bfb67d 100644
--- a/src/generator/input/GENIEInputHandler.cxx
+++ b/src/generator/input/GENIEInputHandler.cxx
@@ -1,165 +1,215 @@
#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), fFileWeight(1) {}
GENIEInputHandler::GENIEInputHandler(GENIEInputHandler &&other)
: fInputTree(std::move(other.fInputTree)),
fReaderEvent(std::move(other.fReaderEvent)), fGenieNtpl(nullptr),
fFileWeight(1) {}
+std::set<std::string> GENIEInputHandler::GetSplineList() {
+ std::set<std::string> splinenames;
+ size_t NEvents = 1000; // GetNEvents();
+ size_t ShoutEvery = NEvents / 100;
+ 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<genie::GHepRecord *>(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<std::string>("file"), "gtree");
fInputTree.tree->SetBranchAddress("gmcrec", &fGenieNtpl);
fKeepIntermediates = ps.get<bool>("keep_intermediates", false);
fKeepNuclearParticles = ps.get<bool>("keep_nuclear_particles", false);
fhicl::ParameterSet XSecInfo = ps.get<fhicl::ParameterSet>("xsec_info", {});
if (XSecInfo.has_key("weight")) {
fFileWeight = XSecInfo.get<double>("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("spline_file") ||
+ XSecInfo.has_key("filter_splines_by_event_tree")) {
std::cout
<< "[INFO]: Attempting to build GENIE file weight with input splines."
<< std::endl;
- fFileWeight = genietools::GetFileWeight(XSecInfo) / double(GetNEvents());
+
+ if (XSecInfo.get<bool>("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;
}
}
-MinimalEvent const &GENIEInputHandler::GetMinimalEvent(ev_index_t idx) const {
+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<genie::GHepRecord *>(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<genie::GHepParticle *>((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<genie::GHepRecord *>(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<genie::GHepParticle *>((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<size_t>(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<double, double> 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 312b6b7..a1f3456 100644
--- a/src/generator/input/GENIEInputHandler.hxx
+++ b/src/generator/input/GENIEInputHandler.hxx
@@ -1,69 +1,72 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#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
namespace fhicl {
class ParameterSet;
}
#include <memory>
+#include <set>
class GENIEInputHandler : public IInputHandler {
mutable nuis::utility::TreeFile fInputTree;
mutable nuis::event::FullEvent fReaderEvent;
mutable std::vector<double> 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 &);
+ void GetEntry(ev_index_t) 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;
+ std::set<std::string> GetSplineList();
double GetXSecScaleFactor(std::pair<double, double> 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 9fc688d..476b383 100644
--- a/src/generator/input/NEUTInputHandler.cxx
+++ b/src/generator/input/NEUTInputHandler.cxx
@@ -1,260 +1,264 @@
#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<std::string>("file"), "neuttree");
fNeutVect = nullptr;
fInputTreeFile.tree->SetBranchAddress("vectorbranch", &fNeutVect);
fKeepIntermediates = ps.get<bool>("keep_intermediates", false);
std::string flux_name = ps.get<std::string>("flux_hist_name", "flux_numu");
std::string evtrt_name = ps.get<std::string>("evtrt_hist_name", "evtrt_numu");
fXSecRescaleFactor = ps.get<double>("xsec_rescale_factor",1);
std::string override_flux_file =
ps.get<std::string>("override_flux_file", "");
if (override_flux_file.size()) {
fFlux = GetHistogramFromROOTFile<TH1>(override_flux_file, flux_name);
RebuildEventRate(!ps.get<bool>("flux_hist_in_MeV", false));
} else {
fFlux =
GetHistogramFromROOTFile<TH1>(fInputTreeFile.file, flux_name, false);
}
if (fFlux) {
if (!fEvtrt) {
fEvtrt = GetHistogramFromROOTFile<TH1>(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<double>::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<double>::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<double, double> 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 {
+void NEUTInputHandler::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();
}
fInputTreeFile.tree->GetEntry(idx);
+}
+
+MinimalEvent const &NEUTInputHandler::GetMinimalEvent(ev_index_t idx) const {
+ 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<size_t>(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<size_t>(
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/input/NEUTInputHandler.hxx b/src/generator/input/NEUTInputHandler.hxx
index a2cc915..574e2fa 100644
--- a/src/generator/input/NEUTInputHandler.hxx
+++ b/src/generator/input/NEUTInputHandler.hxx
@@ -1,71 +1,72 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#pragma once
#include "event/FullEvent.hxx"
#include "input/IInputHandler.hxx"
#include "utility/ROOTUtility.hxx"
#include <memory>
class NeutVect;
class TH1;
namespace fhicl {
class ParameterSet;
}
class NEUTInputHandler : public IInputHandler {
mutable nuis::utility::TreeFile fInputTreeFile;
mutable nuis::event::FullEvent fReaderEvent;
mutable NeutVect *fNeutVect;
mutable std::vector<double> fWeightCache;
double fXSecRescaleFactor;
std::unique_ptr<TH1> fFlux;
std::unique_ptr<TH1> fEvtrt;
double fFileWeight;
bool fKeepIntermediates;
double GetMonoEXSecWeight();
void RebuildEventRate(bool FluxInGeV=true);
public:
NEW_NUIS_EXCEPT(weight_cache_miss);
NEUTInputHandler();
NEUTInputHandler(NEUTInputHandler const &) = delete;
NEUTInputHandler(NEUTInputHandler &&);
void Initialize(fhicl::ParameterSet const &);
NeutVect *GetNeutEvent(ev_index_t) const;
+ void GetEntry(ev_index_t) 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<double, double> const &EnuRange) const;
nuis::GeneratorManager::Generator_id_t GetGeneratorId() const;
};
diff --git a/src/generator/input/NuWroInputHandler.cxx b/src/generator/input/NuWroInputHandler.cxx
index 15e3814..ae0d48c 100644
--- a/src/generator/input/NuWroInputHandler.cxx
+++ b/src/generator/input/NuWroInputHandler.cxx
@@ -1,134 +1,139 @@
#include "generator/input/NuWroInputHandler.hxx"
#include "generator/utility/NuWroUtility.hxx"
#include "fhiclcpp/ParameterSet.h"
#include "particle.h"
using NuWroParticle = ::particle;
using namespace nuis;
using namespace nuis::event;
using namespace nuis::utility;
using namespace nuis::nuwrotools;
NuWroInputHandler::NuWroInputHandler() : fInputTree(), fTreeEvent(nullptr) {}
NuWroInputHandler::NuWroInputHandler(NuWroInputHandler &&other)
: fInputTree(std::move(other.fInputTree)),
fReaderEvent(std::move(other.fReaderEvent)),
fTreeEvent(other.fTreeEvent) {}
void NuWroInputHandler::Initialize(fhicl::ParameterSet const &ps) {
fInputTree = CheckGetTTree(ps.get<std::string>("file"), "treeout");
fTreeEvent = nullptr;
fInputTree.tree->SetBranchAddress("e", &fTreeEvent);
fKeepIntermediates = ps.get<bool>("keep_intermediates", false);
- if(GetNEvents()){
+ if (GetNEvents()) {
fInputTree.tree->GetEntry(0);
- std::cout << "[INFO]: Average NuWro XSec weight: " << (fTreeEvent->weight / double(GetNEvents())) << std::endl;
+ std::cout << "[INFO]: Average NuWro XSec weight: "
+ << (fTreeEvent->weight / double(GetNEvents())) << std::endl;
}
}
-MinimalEvent const &NuWroInputHandler::GetMinimalEvent(ev_index_t idx) const {
+
+void NuWroInputHandler::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();
}
fInputTree.tree->GetEntry(idx);
+}
+
+MinimalEvent const &NuWroInputHandler::GetMinimalEvent(ev_index_t idx) const {
fReaderEvent.mode = NuWroEventChannel(*fTreeEvent);
fReaderEvent.probe_E = fTreeEvent->in[0].E();
fReaderEvent.probe_pdg = fTreeEvent->in[0].pdg;
fReaderEvent.XSecWeight = fTreeEvent->weight / double(GetNEvents());
if (fWeightCache.size() <= idx) {
fWeightCache.push_back(fReaderEvent.XSecWeight);
}
return fReaderEvent;
}
FullEvent const &NuWroInputHandler::GetFullEvent(ev_index_t idx) const {
(void)GetMinimalEvent(idx);
fReaderEvent.ClearParticleStack();
for (size_t p_it = 0; p_it < fTreeEvent->in.size(); ++p_it) {
NuWroParticle &part = fTreeEvent->in[p_it];
Particle nuis_part;
nuis_part.pdg = part.pdg;
nuis_part.P4 = TLorentzVector(part[1], part[2], part[3], part[0]);
fReaderEvent
.ParticleStack[static_cast<size_t>(
Particle::Status_t::kPrimaryInitialState)]
.particles.push_back(nuis_part);
}
for (size_t p_it = 0; p_it < fKeepIntermediates && fTreeEvent->out.size();
++p_it) {
NuWroParticle &part = fTreeEvent->out[p_it];
Particle nuis_part;
nuis_part.pdg = part.pdg;
nuis_part.P4 = TLorentzVector(part[1], part[2], part[3], part[0]);
fReaderEvent
.ParticleStack[static_cast<size_t>(
Particle::Status_t::kPrimaryFinalState)]
.particles.push_back(nuis_part);
}
for (size_t p_it = 0; (p_it < fTreeEvent->post.size()); ++p_it) {
NuWroParticle &part = fTreeEvent->post[p_it];
Particle nuis_part;
nuis_part.pdg = part.pdg;
nuis_part.P4 = TLorentzVector(part[1], part[2], part[3], part[0]);
fReaderEvent
.ParticleStack[static_cast<size_t>(Particle::Status_t::kNuclearLeaving)]
.particles.push_back(nuis_part);
}
fReaderEvent.fGenEvent = fTreeEvent;
return fReaderEvent;
}
NuWroEvent const &NuWroInputHandler::GetNuWroEvent(ev_index_t idx) const {
(void)GetMinimalEvent(idx);
return *fTreeEvent;
}
double NuWroInputHandler::GetEventWeight(ev_index_t idx) const {
if (idx > fWeightCache.size()) {
throw weight_cache_miss()
<< "[ERROR]: Failed to get cached weight for event index: " << idx;
}
return fWeightCache[idx];
}
double NuWroInputHandler::GetXSecScaleFactor(
std::pair<double, double> const &EnuRange) const {
throw input_handler_feature_unimplemented()
<< "[ERROR]: Flux cuts not yet implemented for NuWro input handler.";
}
size_t NuWroInputHandler::GetNEvents() const {
return fInputTree.tree->GetEntries();
}
GeneratorManager::Generator_id_t NuWroInputHandler::GetGeneratorId() const {
return GeneratorManager::Get().EnsureGeneratorRegistered("NuWro");
}
DECLARE_PLUGIN(IInputHandler, NuWroInputHandler);
diff --git a/src/generator/input/NuWroInputHandler.hxx b/src/generator/input/NuWroInputHandler.hxx
index ec36b0c..143a39c 100644
--- a/src/generator/input/NuWroInputHandler.hxx
+++ b/src/generator/input/NuWroInputHandler.hxx
@@ -1,66 +1,67 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#pragma once
#include "event/FullEvent.hxx"
#include "input/IInputHandler.hxx"
#include "exception/exception.hxx"
#include "utility/ROOTUtility.hxx"
#include "event1.h"
using NuWroEvent = ::event;
#include <memory>
namespace fhicl {
class ParameterSet;
}
class NuWroInputHandler : public IInputHandler {
mutable nuis::utility::TreeFile fInputTree;
mutable nuis::event::FullEvent fReaderEvent;
mutable std::vector<double> fWeightCache;
mutable NuWroEvent *fTreeEvent;
bool fKeepIntermediates;
public:
NEW_NUIS_EXCEPT(weight_cache_miss);
NuWroInputHandler();
NuWroInputHandler(NuWroInputHandler const &) = delete;
NuWroInputHandler(NuWroInputHandler &&);
void Initialize(fhicl::ParameterSet const &);
+ void GetEntry(ev_index_t) const;
nuis::event::MinimalEvent const &GetMinimalEvent(ev_index_t idx) const;
nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const;
NuWroEvent const &GetNuWroEvent(ev_index_t idx) const;
double GetEventWeight(ev_index_t idx) const;
size_t GetNEvents() const;
double GetXSecScaleFactor(std::pair<double, double> const &EnuRange) const;
nuis::GeneratorManager::Generator_id_t GetGeneratorId() const;
};
diff --git a/src/generator/utility/GENIEUtility.cxx b/src/generator/utility/GENIEUtility.cxx
index b3be697..6011309 100644
--- a/src/generator/utility/GENIEUtility.cxx
+++ b/src/generator/utility/GENIEUtility.cxx
@@ -1,465 +1,479 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#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"
+#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 {
-static std::map<std::string, TGraph> 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<genie::GHepParticle *>(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) {
+double GetFileWeight(fhicl::ParameterSet const &xsecinfo,
+ std::set<std::string> const &spline_list) {
double mec_scale = xsecinfo.get<double>("mec_scale", 1);
double EMin = xsecinfo.get<double>("EMin", 0);
double EMax = xsecinfo.get<double>("EMax", 10);
size_t NKnots = xsecinfo.get<size_t>("NKnots", 100);
std::vector<TargetSplineBlob> TargetSplines;
int probe = xsecinfo.get<int>("nu_pdg");
for (std::string const &target :
xsecinfo.get<std::vector<std::string>>("target")) {
std::vector<std::string> splits = nuis::utility::split(target, ";");
if (splits.size() != 1 && splits.size() != 2) {
throw Invalid_target_specifier()
<< "Expected to find 100ZZZAAA0;<MolecularWeight>, e.g. "
"1000060120;1, "
"but found: \""
<< target << "\"";
}
std::stringstream ss;
ss << ".*nu:" << probe << ";tgt:" << splits[0] << ".*";
size_t nuc_pdg = fhicl::string_parsers::str2T<size_t>(splits[0]);
double molwght = splits.size() == 2
? fhicl::string_parsers::str2T<double>(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<std::string> 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<std::string>("spline_file")
<< ", this may take a while..." << std::endl;
GENIESplineGetter spg(spline_patterns);
spg.ParseFile(xsecinfo.get<std::string>("spline_file").c_str());
std::vector<std::vector<TGraph>> 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) {
+ 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;
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) {
+ 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
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) {
+ 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 (Int_t p_it = 0; p_it < NKnots; ++p_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;
MasterSpline.SetPoint(p_it, E, XSec);
}
}
+ nuis::utility::TFile_ptr outfile(nullptr, [](TFile *) {});
+ if (xsecinfo.has_key("spline_output_file")) {
+ outfile = nuis::utility::CheckOpenTFile(
+ xsecinfo.get<std::string>("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");
+ }
+
fhicl::ParameterSet fluxps = xsecinfo.get<fhicl::ParameterSet>("flux");
// If MonoE
if (fluxps.has_key("energy_GeV")) {
double monoE = fluxps.get<double>("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<TH1> Flux =
nuis::utility::GetHistogram<TH1>(fluxps.get<std::string>("histogram"));
+
+ if (outfile) {
+ std::unique_ptr<TH1> Fluxcp = nuis::utility::Clone<TH1>(Flux);
+ nuis::utility::WriteToTFile(outfile, Fluxcp.get(), "Flux");
+ }
+
bool per_width = fluxps.get<bool>("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);
}
+ if (outfile) {
+ std::unique_ptr<TH1> Fluxcp = nuis::utility::Clone<TH1>(Flux);
+ nuis::utility::WriteToTFile(outfile, Fluxcp.get(), "EvRate");
+ }
+
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 3bc447a..bb25aef 100644
--- a/src/generator/utility/GENIEUtility.hxx
+++ b/src/generator/utility/GENIEUtility.hxx
@@ -1,53 +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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#pragma once
#include "event/Particle.hxx"
#include "event/types.hxx"
+#include <set>
+
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);
-double GetFileWeight(fhicl::ParameterSet const &xsecinfo);
+double GetFileWeight(fhicl::ParameterSet const &xsecinfo,
+ std::set<std::string> const &spline_list = {});
} // namespace genietools
} // namespace nuis
diff --git a/src/input/IInputHandler.hxx b/src/input/IInputHandler.hxx
index fcb467e..205ea65 100644
--- a/src/input/IInputHandler.hxx
+++ b/src/input/IInputHandler.hxx
@@ -1,122 +1,123 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#pragma once
#include "plugins/traits.hxx"
#include "exception/exception.hxx"
#include "generator/GeneratorManager.hxx"
#include <iterator>
#include <limits>
namespace fhicl {
class ParameterSet;
}
namespace nuis {
namespace event {
class MinimalEvent;
class FullEvent;
} // namespace event
} // namespace nuis
class IInputHandler {
public:
struct FullEvent_const_iterator
: public std::iterator<
std::input_iterator_tag, nuis::event::FullEvent const, size_t,
nuis::event::FullEvent const *, nuis::event::FullEvent const &> {
FullEvent_const_iterator(size_t _idx, IInputHandler const *_ih) {
idx = _idx;
ih = _ih;
}
FullEvent_const_iterator(FullEvent_const_iterator const &other) {
idx = other.idx;
ih = other.ih;
}
FullEvent_const_iterator operator=(FullEvent_const_iterator const &other) {
idx = other.idx;
ih = other.ih;
return (*this);
}
bool operator==(FullEvent_const_iterator const &other) {
return (idx == other.idx);
}
bool operator!=(FullEvent_const_iterator const &other) {
return !(*this == other);
}
nuis::event::FullEvent const &operator*() { return ih->GetFullEvent(idx); }
nuis::event::FullEvent const *operator->() {
return &ih->GetFullEvent(idx);
}
FullEvent_const_iterator operator++() {
idx++;
return *this;
}
FullEvent_const_iterator operator++(int) {
FullEvent_const_iterator tmp(*this);
idx++;
return tmp;
}
private:
size_t idx;
IInputHandler const *ih;
};
NEW_NUIS_EXCEPT(invalid_input_file);
NEW_NUIS_EXCEPT(invalid_entry);
NEW_NUIS_EXCEPT(input_handler_feature_unimplemented);
typedef size_t ev_index_t;
virtual void Initialize(fhicl::ParameterSet const &) = 0;
+ virtual void GetEntry(ev_index_t) const = 0;
virtual nuis::event::MinimalEvent const &
GetMinimalEvent(ev_index_t idx) const = 0;
virtual nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const = 0;
virtual void RecalculateEventWeights(){};
virtual double GetEventWeight(ev_index_t) const { return 1; };
/// Allows samples that implement flux cuts to request an appropriate
virtual double GetXSecScaleFactor(
std::pair<double, double> const &EnuRange = std::pair<double, double>{
std::numeric_limits<double>::max(),
std::numeric_limits<double>::max()}) const = 0;
virtual size_t GetNEvents() const = 0;
FullEvent_const_iterator begin() const {
return FullEvent_const_iterator(0, this);
}
FullEvent_const_iterator end() const {
return FullEvent_const_iterator(GetNEvents(), this);
}
virtual ~IInputHandler() {}
virtual nuis::GeneratorManager::Generator_id_t GetGeneratorId() const = 0;
};
DECLARE_PLUGIN_INTERFACE(IInputHandler);
diff --git a/src/tests/CMakeLists.txt b/src/tests/CMakeLists.txt
new file mode 100644
index 0000000..c0f8bcb
--- /dev/null
+++ b/src/tests/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(memcheck)
diff --git a/src/tests/memcheck/CMakeLists.txt b/src/tests/memcheck/CMakeLists.txt
new file mode 100644
index 0000000..c3fb4e8
--- /dev/null
+++ b/src/tests/memcheck/CMakeLists.txt
@@ -0,0 +1,30 @@
+SET(TESTS InputHandlerCheck)
+
+foreach(a ${TESTS})
+ add_executable(${a} ${a}.cxx)
+ 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 InputHandlers)
+ 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()
+
+ if(USE_GENIE)
+ target_compile_options(${a} PUBLIC ${GENIE_CXX_FLAGS})
+ include_directories(${GENIE_INCLUDE_DIRS})
+ endif()
+
+ if(USE_NEUT)
+ target_compile_options(${a} PUBLIC ${NEUT_CXX_FLAGS})
+ include_directories(${NEUT_INCLUDE_DIRS})
+ endif()
+
+ if(USE_NUWRO)
+ target_compile_options(${a} PUBLIC ${NUWRO_CXX_FLAGS})
+ include_directories(${NUWRO_INCLUDE_DIRS})
+ endif()
+
+ install(TARGETS ${a} DESTINATION test)
+endforeach()
diff --git a/src/tests/memcheck/InputHandlerCheck.cxx b/src/tests/memcheck/InputHandlerCheck.cxx
new file mode 100644
index 0000000..9a795ad
--- /dev/null
+++ b/src/tests/memcheck/InputHandlerCheck.cxx
@@ -0,0 +1,156 @@
+#include "config/GlobalConfiguration.hxx"
+
+#ifdef USE_GENIE
+#include "generator/input/GENIEInputHandler.hxx"
+#endif
+#ifdef USE_NEUT
+#include "generator/input/NEUTInputHandler.hxx"
+#endif
+#include "input/InputManager.hxx"
+
+#include "event/MinimalEvent.hxx"
+
+#include "exception/exception.hxx"
+
+#include "utility/StringUtility.hxx"
+
+#include "fhiclcpp/make_ParameterSet.h"
+#include "string_parsers/from_string.hxx"
+
+#include <string>
+
+NEW_NUIS_EXCEPT(invalid_cli_arguments);
+
+void SayUsage(char const *argv[]) {
+ std::cout
+ << "[USAGE]: " << argv[0]
+ << "\n"
+ "\t-c <config.fcl> : FHiCL file containing study "
+ "configuration. \n"
+ "\t-s <sample name> : FHiCL key of a single sample "
+ "to run from the -c argument. \n"
+ "\t-n <NTries> : Number of events toload. \n"
+#ifdef USE_GENIE
+ "\t--genie : Just use GENIEInputHandler directly "
+#endif
+#ifdef USE_NEUT
+ "\t--neut : Just use NEUTInputHandler directly "
+#endif
+ "\t-G : Just use GetEntry directly "
+ "rather than the GetMinimalEvent interface."
+ << std::endl;
+}
+
+std::string fhicl_file = "";
+std::string named_sample = "";
+size_t NTries = 0;
+bool GetEntry = false;
+
+enum EGen {
+ kManager,
+#ifdef USE_NEUT
+ kNEUT,
+#endif
+#ifdef USE_GENIE
+ kGENIE,
+#endif
+};
+EGen WhichGen = kManager;
+
+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]) == "-n") {
+ NTries = fhicl::string_parsers::str2T<size_t>(argv[++opt]);
+ } else if (std::string(argv[opt]) == "-G") {
+ GetEntry = true;
+ }
+#ifdef USE_GENIE
+ else if (std::string(argv[opt]) == "--genie") {
+ WhichGen = kGENIE;
+ }
+#endif
+#ifdef USE_NEUT
+ else if (std::string(argv[opt]) == "--neut") {
+ WhichGen = kNEUT;
+ }
+#endif
+ 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);
+
+ fhicl::ParameterSet inps =
+ nuis::config::GetDocument().get<fhicl::ParameterSet>(named_sample);
+ IInputHandler const *IH;
+
+ switch (WhichGen) {
+ case kManager: {
+ auto id = nuis::input::InputManager::Get().EnsureInputLoaded(inps);
+ IH = &nuis::input::InputManager::Get().GetInputHandler(id);
+ break;
+ }
+#ifdef USE_NEUT
+ case kNEUT: {
+ IInputHandler *IIH = new NEUTInputHandler();
+ IIH->Initialize(inps);
+ IH = IIH;
+ break;
+ }
+#endif
+#ifdef USE_GENIE
+ case kGENIE: {
+ IInputHandler *IIH = new GENIEInputHandler();
+ IIH->Initialize(inps);
+ IH = IIH;
+ break;
+ }
+#endif
+ default: {
+ std::cout << "[ERROR]: Invalid input handler type." << std::endl;
+ return 1;
+ }
+ }
+
+ size_t ShoutEvery = NTries / 100;
+
+ std::cout << "[INFO]: Read " << 0 << "/" << NTries << " GENIE events."
+ << std::flush;
+ for (size_t t_it = 0; t_it < NTries; ++t_it) {
+ if (ShoutEvery && !(t_it % ShoutEvery)) {
+ std::cout << "\r[INFO]: Read " << t_it << "/" << NTries
+ << " GENIE events." << std::flush;
+ }
+ if (GetEntry) {
+ IH->GetEntry(0);
+ } else {
+ IH->GetMinimalEvent(0);
+ }
+ }
+ std::cout << std::endl;
+}
diff --git a/src/utility/ROOTUtility.hxx b/src/utility/ROOTUtility.hxx
index 75b0d47..2099789 100644
--- a/src/utility/ROOTUtility.hxx
+++ b/src/utility/ROOTUtility.hxx
@@ -1,164 +1,181 @@
// 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 <http://www.gnu.org/licenses/>.
*******************************************************************************/
#ifndef UTILITY_ROOTUTILITY_HXX_SEEN
#define UTILITY_ROOTUTILITY_HXX_SEEN
#include "exception/exception.hxx"
#include "TFile.h"
#include "TTree.h"
#include <iomanip>
#include <iostream>
#include <memory>
#include <string>
namespace nuis {
namespace utility {
NEW_NUIS_EXCEPT(failed_to_open_TFile);
NEW_NUIS_EXCEPT(failed_to_get_TTree);
+NEW_NUIS_EXCEPT(WriteToTFile_nullptr);
inline void CloseTFile(TFile *f = nullptr) {
if (f) {
std::cout << "[INFO]: Shutting TFile: " << f->GetName() << ", " << f
<< std::endl;
f->Close();
}
delete f;
}
using TFile_ptr = std::unique_ptr<TFile, decltype(&CloseTFile)>;
inline TFile_ptr make_unique_TFile(TFile *f) {
return TFile_ptr(f, &CloseTFile);
}
inline TFile_ptr CheckOpenTFile(std::string const &fname,
char const *opts = "") {
TDirectory *ogDir = gDirectory;
TFile *inpF = new TFile(fname.c_str(), opts);
if (!inpF || !inpF->IsOpen()) {
throw failed_to_open_TFile()
<< "[ERROR]: Couldn't open input file: " << std::quoted(fname);
}
if (ogDir) {
ogDir->cd();
}
return make_unique_TFile(inpF);
}
struct TreeFile {
TFile_ptr file;
TTree *tree;
bool file_owned;
TreeFile()
: file(make_unique_TFile(nullptr)), tree(nullptr), file_owned(false) {}
TreeFile(TreeFile const &other) = delete;
TreeFile &operator=(TreeFile &&other) {
file = std::move(other.file);
tree = other.tree;
file_owned = other.file_owned;
other.file = nullptr;
other.tree = nullptr;
other.file_owned = false;
return *this;
}
TreeFile(TreeFile &&other)
: file(std::move(other.file)), tree(other.tree),
file_owned(other.file_owned) {
other.file = nullptr;
other.tree = nullptr;
other.file_owned = false;
}
TTree *operator->() { return tree; }
~TreeFile() {
if (!file_owned) {
file.release();
}
}
};
inline TreeFile CheckGetTTree(TFile *file, std::string const &tname) {
TreeFile tf;
tf.file = make_unique_TFile(file);
tf.tree = dynamic_cast<TTree *>(tf.file->Get(tname.c_str()));
tf.file_owned = false;
if (!tf.tree) {
throw failed_to_get_TTree()
<< "[ERROR]: Failed to get TTree named: " << std::quoted(tname)
<< " from TFile: " << std::quoted(file->GetName());
}
std::cout << "[INFO]: Opened TFile: " << file->GetName() << ", " << file
<< std::endl;
return tf;
}
inline TreeFile CheckGetTTree(std::string const &fname,
std::string const &tname, char const *opts = "") {
TreeFile tf = CheckGetTTree(CheckOpenTFile(fname, opts).release(), tname);
tf.file_owned = true;
return tf;
}
inline TreeFile MakeNewTTree(std::string const &fname, std::string const &tname,
char const *opts = "") {
TreeFile tf;
tf.file = CheckOpenTFile(fname, opts);
tf.tree = new TTree(tname.c_str(), "");
tf.tree->SetDirectory(tf.file.get());
tf.file_owned = true;
return tf;
}
inline TreeFile AddNewTTreeToFile(TFile *f, std::string const &tname,
std::string const &dirname = "") {
TreeFile tf;
// This is pretty dumb...
tf.file = make_unique_TFile(f);
tf.tree = new TTree(tname.c_str(), "");
if (dirname.size()) {
TDirectory *dir = tf.file->GetDirectory(dirname.c_str());
if (!dir) {
dir = tf.file->mkdir(dirname.c_str());
}
tf.tree->SetDirectory(dir);
} else {
tf.tree->SetDirectory(tf.file.get());
}
tf.file_owned = false;
return tf;
}
inline std::string SanitizeROOTObjectName(std::string name) { return name; }
+inline void WriteToTFile(TFile_ptr &tf, TObject *object,
+ std::string const &object_name) {
+ if (!object) {
+ throw WriteToTFile_nullptr();
+ }
+
+ TDirectory *ogdir = gDirectory;
+
+ tf->cd();
+ object->Write(object_name.c_str(), TObject::kOverwrite);
+
+ if (ogdir) {
+ ogdir->cd();
+ }
+}
+
} // namespace utility
} // namespace nuis
#endif