diff --git a/scripts/to_configure/nuisbuild.in b/scripts/to_configure/nuisbuild.in
index 316b790..1daaeb1 100644
--- a/scripts/to_configure/nuisbuild.in
+++ b/scripts/to_configure/nuisbuild.in
@@ -1,101 +1,109 @@
# 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 .
################################################################################
#!/bin/bash
EXTRA_COMPILE_FLAGS=""
EXTRA_LINK_FLAGS=""
INPUT_FILE_NAME=""
OUTPUT_LIB_NAME=""
+DEBUG_FLAGS=""
SCRIPTNAME=${0}
while [[ ${#} -gt 0 ]]; do
key="$1"
case $key in
-i|--implementation-file)
if [[ ${#} -lt 2 ]]; then
echo "[ERROR]: ${1} expected a value."
exit 1
fi
INPUT_FILE_NAME="$2"
echo "[OPT]: Reading event processor implementation from: ${INPUT_FILE_NAME}"
shift # past argument
;;
--compile-flags)
if [[ ${#} -lt 2 ]]; then
echo "[ERROR]: ${1} expected a value."
exit 1
fi
EXTRA_COMPILE_FLAGS="$2"
echo "[OPT]: Including extra compile flags: ${EXTRA_COMPILE_FLAGS}"
shift # past argument
;;
+ -g|--debug)
+
+ DEBUG_FLAGS=" -g -O0 "
+ echo "[OPT]: Building with debugging symbols."
+ ;;
+
--link-flags)
if [[ ${#} -lt 2 ]]; then
echo "[ERROR]: ${1} expected a value."
exit 1
fi
EXTRA_LINK_FLAGS="$2"
echo "[OPT]: Including extra linker flags: ${EXTRA_LINK_FLAGS}"
shift # past argument
;;
-?|--help)
echo "[RUNLIKE] ${SCRIPTNAME}"
echo -e "\t-i|--implementation-file : Path to IEventProcessor subclass to build into dynamically loadable shared object."
echo -e "\t--compile-flags : Extra compile flags to pass to the compiler."
echo -e "\t--link-flags : Extra linker flags to pass to the linker."
+ echo -e "\t-g|--debug : Add default debugging compiler flags.."
echo -e "\t-?|--help : Print this message."
exit 0
;;
*)
# unknown option
echo "Unknown option $1"
exit 1
;;
esac
shift
done
OUTPUT_LIB_NAME=${INPUT_FILE_NAME%.*}.so
#Removes the NUISANCE filename descriptor that uses a gmake trick to insert the file name into the compiler options
NUISANCE_CXX_FLAGS_GMAKE_SCRUB=$(echo '@NUISANCE_CXX_FLAGS@' | sed s/-D__FILENAME__=\".*\"//g)
-echo "@CMAKE_CXX_COMPILER@ ${INPUT_FILE_NAME} -shared -fPIC -o ${OUTPUT_LIB_NAME} ${NUISANCE_CXX_FLAGS_GMAKE_SCRUB} ${EXTRA_COMPILE_FLAGS} -D__FILENAME__=\"${INPUT_FILE_NAME##*/}\" -I. -I@CMAKE_INSTALL_PREFIX@/include @CMAKE_LINK_FLAGS@ ${EXTRA_LINK_FLAGS} -L@CMAKE_INSTALL_PREFIX@/lib @NUISANCE_LINK_DIRS@ @NUISANCE_DEPEND_LIBS@"
+echo "@CMAKE_CXX_COMPILER@ ${INPUT_FILE_NAME} -shared -fPIC -o ${OUTPUT_LIB_NAME} ${NUISANCE_CXX_FLAGS_GMAKE_SCRUB} ${EXTRA_COMPILE_FLAGS} -D__FILENAME__=\"${INPUT_FILE_NAME##*/}\" -I. -I@CMAKE_INSTALL_PREFIX@/include ${DEBUG_FLAGS} @CMAKE_LINK_FLAGS@ ${EXTRA_LINK_FLAGS} -L@CMAKE_INSTALL_PREFIX@/lib @NUISANCE_LINK_DIRS@ @NUISANCE_DEPEND_LIBS@"
-if ! @CMAKE_CXX_COMPILER@ ${INPUT_FILE_NAME} -shared -fPIC -o ${OUTPUT_LIB_NAME} ${NUISANCE_CXX_FLAGS_GMAKE_SCRUB} ${EXTRA_COMPILE_FLAGS} -D__FILENAME__=\"${INPUT_FILE_NAME##*/}\" -I. -I@CMAKE_INSTALL_PREFIX@/include @CMAKE_LINK_FLAGS@ ${EXTRA_LINK_FLAGS} -L@CMAKE_INSTALL_PREFIX@/lib @NUISANCE_LINK_DIRS@ @NUISANCE_DEPEND_LIBS@; then
+if ! @CMAKE_CXX_COMPILER@ ${INPUT_FILE_NAME} -shared -fPIC -o ${OUTPUT_LIB_NAME} ${NUISANCE_CXX_FLAGS_GMAKE_SCRUB} ${EXTRA_COMPILE_FLAGS} -D__FILENAME__=\"${INPUT_FILE_NAME##*/}\" -I. -I@CMAKE_INSTALL_PREFIX@/include ${DEBUG_FLAGS} @CMAKE_LINK_FLAGS@ ${EXTRA_LINK_FLAGS} -L@CMAKE_INSTALL_PREFIX@/lib @NUISANCE_LINK_DIRS@ @NUISANCE_DEPEND_LIBS@; then
echo "[ERROR]: Failed to compile ${INPUT_FILE_NAME} -> ${OUTPUT_LIB_NAME}"
else
echo "[INFO]: Successfully build: ${OUTPUT_LIB_NAME}."
fi
diff --git a/src/generator/input/NEUTInputHandler.cxx b/src/generator/input/NEUTInputHandler.cxx
index dda5a6d..dbd85bd 100644
--- a/src/generator/input/NEUTInputHandler.cxx
+++ b/src/generator/input/NEUTInputHandler.cxx
@@ -1,129 +1,198 @@
#include "generator/input/NEUTInputHandler.hxx"
#include "utility/HistogramUtility.hxx"
#include "utility/InteractionChannelUtility.hxx"
#include "utility/PDGCodeUtility.hxx"
#include "generator/utility/NEUTUtility.hxx"
#include "neutpart.h"
#include "neutvect.h"
#include "fhiclcpp/ParameterSet.h"
using namespace nuis;
using namespace nuis::event;
using namespace nuis::utility;
using namespace nuis::neuttools;
NEUTInputHandler::NEUTInputHandler() : fInputTreeFile() {}
NEUTInputHandler::NEUTInputHandler(NEUTInputHandler &&other)
: fInputTreeFile(std::move(other.fInputTreeFile)),
fReaderEvent(std::move(other.fReaderEvent)) {}
void NEUTInputHandler::Initialize(fhicl::ParameterSet const &ps) {
fInputTreeFile = CheckGetTTree(ps.get("file"), "neuttree");
fNeutVect = nullptr;
fInputTreeFile.tree->SetBranchAddress("vectorbranch", &fNeutVect);
fKeepIntermediates = ps.get("keep_intermediates", false);
- fFlux = GetHistogramFromROOTFile(fInputTreeFile.file, "flux_numu");
- fEvtrt = GetHistogramFromROOTFile(fInputTreeFile.file, "evtrt_numu");
+ fFlux =
+ GetHistogramFromROOTFile(fInputTreeFile.file, "flux_numu", false);
+ if (fFlux) {
+ fEvtrt = GetHistogramFromROOTFile(fInputTreeFile.file, "evtrt_numu");
- fFileWeight = fEvtrt->Integral() / (fFlux->Integral() * double(GetNEvents()));
+ fFileWeight =
+ fEvtrt->Integral() / (fFlux->Integral() * double(GetNEvents()));
+ } 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;
+ }
+}
+
+NEW_NUIS_EXCEPT(non_mono_energetic_input_file);
+double NEUTInputHandler::GetMonoEXSecWeight() {
+ double xsec = 0;
+ double count = 0;
+ double E_first = std::numeric_limits::max();
+ size_t NEvents = GetNEvents();
+ size_t ShoutEvery = NEvents / 100;
+ std::cout << "[INFO]: Read " << 0 << "/" << NEvents << " NEUT events."
+ << std::flush;
+
+ for (size_t nev_it = 0; nev_it < NEvents; ++nev_it) {
+ if (ShoutEvery && !(nev_it % ShoutEvery)) {
+ std::cout << "\r[INFO]: Read " << nev_it << "/" << NEvents
+ << " NEUT events." << std::flush;
+ }
+
+ fInputTreeFile.tree->GetEntry(nev_it);
+ xsec += fNeutVect->Totcrs;
+ count++;
+ if (E_first == std::numeric_limits::max()) {
+ NeutPart *part = fNeutVect->PartInfo(0);
+ E_first = part->fP.E();
+ } else {
+ NeutPart *part = fNeutVect->PartInfo(0);
+ if (E_first != part->fP.E()) {
+ throw non_mono_energetic_input_file()
+ << "[ERROR]: When calculating NEUT xsec weight for a mono "
+ "energetic beam, found different energies: "
+ << E_first << " and " << part->fP.E()
+ << ". Cannot continue with this input file.";
+ }
+ }
+ }
+ std::cout << "\r[INFO]: Read " << NEvents << "/" << NEvents << " NEUT events."
+ << std::endl;
+
+ return (xsec / count)*1E-38;
}
double
NEUTInputHandler::GetXSecScaleFactor(std::pair const &) const {
throw input_handler_feature_unimplemented()
<< "[ERROR]: Flux cuts not yet implemented for NEUT input handler.";
}
+NeutVect *NEUTInputHandler::GetNeutEvent(ev_index_t idx) const {
+ if (idx >= GetNEvents()) {
+ throw IInputHandler::invalid_entry()
+ << "[ERROR]: Attempted to get entry " << idx
+ << " from an InputHandler with only " << GetNEvents();
+ }
+ fInputTreeFile.tree->GetEntry(idx);
+
+ return fNeutVect;
+}
+
MinimalEvent const &NEUTInputHandler::GetMinimalEvent(ev_index_t idx) const {
if (idx >= GetNEvents()) {
throw IInputHandler::invalid_entry()
<< "[ERROR]: Attempted to get entry " << idx
<< " from an InputHandler with only " << GetNEvents();
}
fInputTreeFile.tree->GetEntry(idx);
fReaderEvent.mode = IntToChannel(fNeutVect->Mode);
size_t NPart = fNeutVect->Npart();
for (size_t part_it = 0; part_it < NPart; part_it++) {
NeutPart const &part = (*fNeutVect->PartInfo(part_it));
if ((part.fIsAlive == false) && (part.fStatus == -1) &&
IsNeutralLepton(part.fPID)) {
fReaderEvent.probe_E = part.fP.T();
fReaderEvent.probe_pdg = part.fPID;
break;
}
}
fReaderEvent.XSecWeight = fFileWeight;
if (fWeightCache.size() <= idx) {
fWeightCache.push_back(fReaderEvent.XSecWeight);
}
return fReaderEvent;
}
FullEvent const &NEUTInputHandler::GetFullEvent(ev_index_t idx) const {
(void)GetMinimalEvent(idx);
fReaderEvent.ClearParticleStack();
if (fNeutVect->Ibound) {
fReaderEvent.target_pdg =
MakeNuclearPDG(fNeutVect->TargetA, fNeutVect->TargetZ);
} else {
fReaderEvent.target_pdg = MakeNuclearPDG(1, 1);
}
size_t NPart = fNeutVect->Npart();
size_t NPrimary = fNeutVect->Nprimary();
for (size_t part_it = 0; part_it < NPart; part_it++) {
NeutPart const &part = (*fNeutVect->PartInfo(part_it));
Particle nuis_part;
nuis_part.pdg = part.fPID;
nuis_part.P4 = part.fP;
Particle::Status_t state = GetNeutParticleStatus(part, fReaderEvent.mode);
if (!fKeepIntermediates && (state == Particle::Status_t::kIntermediate)) {
continue;
}
size_t state_int = static_cast(state);
// Add status == 0 particles as pre-FSI particles until we find an
// intermediate state particle
if ((part_it < NPrimary) &&
(state != Particle::Status_t::kPrimaryInitialState)) {
fReaderEvent
.ParticleStack[static_cast(
Particle::Status_t::kPrimaryFinalState)]
.particles.push_back(nuis_part);
}
fReaderEvent.ParticleStack[state_int].particles.push_back(nuis_part);
}
return fReaderEvent;
}
+double NEUTInputHandler::GetEventWeight(ev_index_t idx) const {
+ if (idx > fWeightCache.size()) {
+ throw weight_cache_miss()
+ << "[ERROR]: Failed to get cached weight for event index: " << idx;
+ }
+ return fWeightCache[idx];
+}
+
size_t NEUTInputHandler::GetNEvents() const {
return fInputTreeFile.tree->GetEntries();
}
-GeneratorManager::Generator_id_t NEUTInputHandler::GetGeneratorId(){
+GeneratorManager::Generator_id_t NEUTInputHandler::GetGeneratorId() {
return GeneratorManager::Get().EnsureGeneratorRegistered("NEUT");
}
DECLARE_PLUGIN(IInputHandler, NEUTInputHandler);
diff --git a/src/generator/input/NEUTInputHandler.hxx b/src/generator/input/NEUTInputHandler.hxx
index 55cbf8c..7fc1c3f 100644
--- a/src/generator/input/NEUTInputHandler.hxx
+++ b/src/generator/input/NEUTInputHandler.hxx
@@ -1,61 +1,69 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#pragma once
#include "event/FullEvent.hxx"
#include "input/IInputHandler.hxx"
#include "utility/ROOTUtility.hxx"
#include
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 fWeightCache;
std::unique_ptr fFlux;
std::unique_ptr fEvtrt;
double fFileWeight;
bool fKeepIntermediates;
+ double GetMonoEXSecWeight();
+
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;
nuis::event::MinimalEvent const &GetMinimalEvent(ev_index_t idx) const;
nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const;
+
+ double GetEventWeight(ev_index_t idx) const;
size_t GetNEvents() const;
+
double GetXSecScaleFactor(std::pair const &EnuRange) const;
nuis::GeneratorManager::Generator_id_t GetGeneratorId();
};
diff --git a/src/generator/input/NuWroInputHandler.cxx b/src/generator/input/NuWroInputHandler.cxx
index 7769c14..0b50cb7 100644
--- a/src/generator/input/NuWroInputHandler.cxx
+++ b/src/generator/input/NuWroInputHandler.cxx
@@ -1,127 +1,132 @@
#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("file"), "treeout");
fTreeEvent = nullptr;
fInputTree.tree->SetBranchAddress("e", &fTreeEvent);
fKeepIntermediates = ps.get("keep_intermediates", false);
+
+ if(GetNEvents()){
+ fInputTree.tree->GetEntry(0);
+ std::cout << "[INFO]: Average NuWro XSec weight: " << (fTreeEvent->weight / double(GetNEvents())) << std::endl;
+ }
}
MinimalEvent const &NuWroInputHandler::GetMinimalEvent(ev_index_t idx) const {
if (idx >= GetNEvents()) {
throw IInputHandler::invalid_entry()
<< "[ERROR]: Attempted to get entry " << idx
<< " from an InputHandler with only " << GetNEvents();
}
fInputTree.tree->GetEntry(idx);
fReaderEvent.mode = NuWroEventChannel(*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(
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(
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(Particle::Status_t::kNuclearLeaving)]
.particles.push_back(nuis_part);
}
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 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(){
return GeneratorManager::Get().EnsureGeneratorRegistered("NuWro");
}
DECLARE_PLUGIN(IInputHandler, NuWroInputHandler);
diff --git a/src/generator/utility/NEUTUtility.cxx b/src/generator/utility/NEUTUtility.cxx
index 56cd4a9..b7f6310 100644
--- a/src/generator/utility/NEUTUtility.cxx
+++ b/src/generator/utility/NEUTUtility.cxx
@@ -1,90 +1,93 @@
#include "generator/utility/NEUTUtility.hxx"
+#include "generator/input/NEUTInputHandler.hxx"
+
#include "exception/exception.hxx"
#include "utility/InteractionChannelUtility.hxx"
#include "utility/PDGCodeUtility.hxx"
#include "neutpart.h"
+#include "neutvect.h"
using namespace nuis::event;
using namespace nuis::utility;
namespace nuis {
namespace neuttools {
NEW_NUIS_EXCEPT(unexpected_NEUT_particle_state);
Particle::Status_t GetNeutParticleStatus(NeutPart const &part, Channel_t mode) {
#ifdef DEBUG_NEUT_UTILITY
std::cout << "[DEBUG]: Mode: " << mode << ", Part: { Status: " << part.fStatus
<< ", IsAlive: " << part.fIsAlive << ", PDG: " << part.fPID << " }."
<< std::endl;
#endif
// Remove Pauli blocked events, probably just single pion events
if (part.fStatus == 5) {
return Particle::Status_t::kBlocked;
// fStatus == -1 means initial state
} else if (part.fIsAlive == false && part.fStatus == -1) {
return Particle::Status_t::kPrimaryInitialState;
// NEUT has a bit of a strange convention for fIsAlive and fStatus
// combinations
// for NC and neutrino particle isAlive true/false and status 2 means
// final state particle
// for other particles in NC status 2 means it's an FSI particle
// for CC it means it was an FSI particle
} else if (part.fStatus == 2) {
// NC case is a little strange... The outgoing neutrino might be alive or
// not alive. Remaining particles with status 2 are FSI particles that
// reinteracted
if (IsNC(mode) && IsNeutralLepton(part.fPID)) {
return Particle::Status_t::kNuclearLeaving;
// The usual CC case
} else if (part.fIsAlive == true) {
// return Particle::Status_t::kIntermediate;
throw unexpected_NEUT_particle_state()
<< "[ERROR] Found unexpected NEUT particle in neutvect stack: Mode: "
<< mode << ", Part: { Status: " << part.fStatus
<< ", IsAlive: " << part.fIsAlive << ", PDG: " << part.fPID << " }.";
}
} else if ((part.fIsAlive == true) && (part.fStatus == 2) &&
IsNeutralLepton(part.fPID)) {
return Particle::Status_t::kNuclearLeaving;
} else if ((part.fIsAlive == true) && (part.fStatus == 0)) {
return Particle::Status_t::kNuclearLeaving;
} else if (!part.fIsAlive && ((part.fStatus == 1) || (part.fStatus == 3) ||
(part.fStatus == 4) || (part.fStatus == 7) ||
(part.fStatus == 8))) {
return Particle::Status_t::kIntermediate;
// There's one hyper weird case where fStatus = -3. This apparently
// corresponds to a nucleon being ejected via pion FSI when there is "data
// available"
} else if (!part.fIsAlive && (part.fStatus == -3)) {
return Particle::Status_t::kUnknown;
// NC neutrino outgoing
} else if (!part.fIsAlive && part.fStatus == 0 &&
IsNeutralLepton(part.fPID)) {
return Particle::Status_t::kNuclearLeaving;
// Warn if we still find alive particles without classifying them
} else if (part.fIsAlive == true) {
throw unexpected_NEUT_particle_state()
<< "[ERROR]: Undefined NEUT state "
<< " Alive: " << part.fIsAlive << " Status: " << part.fStatus
<< " PDG: " << part.fPID;
}
// Warn if we find dead particles that we haven't classified
throw unexpected_NEUT_particle_state()
<< "[ERROR]: Undefined NEUT state "
<< " Alive: " << part.fIsAlive << " Status: " << part.fStatus
<< " PDG: " << part.fPID;
}
} // namespace neuttools
} // namespace nuis
diff --git a/src/generator/utility/NEUTUtility.hxx b/src/generator/utility/NEUTUtility.hxx
index c191e52..fdcf697 100644
--- a/src/generator/utility/NEUTUtility.hxx
+++ b/src/generator/utility/NEUTUtility.hxx
@@ -1,15 +1,17 @@
-#ifndef GENERATOR_UTILITY_NEUTUTILITY_HXX_SEEN
-#define GENERATOR_UTILITY_NEUTUTILITY_HXX_SEEN
+#pragma once
#include "event/Particle.hxx"
#include "event/types.hxx"
+#include "exception/exception.hxx"
+
+#include
+
class NeutPart;
namespace nuis {
namespace neuttools {
nuis::event::Particle::Status_t GetNeutParticleStatus(NeutPart const &,
nuis::event::Channel_t);
}
} // namespace nuis
-#endif
diff --git a/src/samples/SimpleMCStudy.hxx b/src/samples/SimpleMCStudy.hxx
index c8a55d5..b6da7ae 100644
--- a/src/samples/SimpleMCStudy.hxx
+++ b/src/samples/SimpleMCStudy.hxx
@@ -1,108 +1,111 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#pragma once
#include "samples/IEventProcessor.hxx"
#include "event/FullEvent.hxx"
#include "input/InputManager.hxx"
#include
class SimpleMCStudy : public IEventProcessor {
protected:
nuis::input::InputManager::Input_id_t fIH_id;
std::string write_directory;
size_t NMaxSample_override;
std::function
ProcessEventFunction;
public:
SimpleMCStudy() {
fIH_id = std::numeric_limits::max();
write_directory = "";
NMaxSample_override = std::numeric_limits::max();
- ProcessEventFunction = std::function();
+ ProcessEventFunction =
+ std::function();
}
fhicl::ParameterSet fGlobalConfig;
fhicl::ParameterSet fInstanceConfig;
void ReadGlobalConfigDefaults() {
fGlobalConfig = nuis::config::GetDocument().get(
std::string("global.sample_configuration.") + Name(),
fhicl::ParameterSet());
}
void Initialize(fhicl::ParameterSet const &instance_sample_configuration) {
fInstanceConfig = instance_sample_configuration;
- if (fInstanceConfig.has_key("verbosity")) {
- SetSampleVerbosity(fInstanceConfig.get("verbosity"));
- } else {
- SetSampleVerbosity("Reticent");
- }
+ SetSampleVerbosity(fInstanceConfig.get(
+ "verbosity", nuis::config::GetDocument().get(
+ "global.sample.verbosity_default", "Reticent")));
ReadGlobalConfigDefaults();
NMaxSample_override =
fInstanceConfig.get("nmax", std::numeric_limits::max());
write_directory =
fInstanceConfig.get("write_directory", Name());
fIH_id =
nuis::input::InputManager::Get().EnsureInputLoaded(fInstanceConfig);
}
void ProcessSample(size_t nmax = std::numeric_limits::max()) {
if (fIH_id ==
std::numeric_limits::max()) {
throw uninitialized_IEventProcessor();
}
IInputHandler const &IH =
nuis::input::InputManager::Get().GetInputHandler(fIH_id);
nmax = std::min(NMaxSample_override, nmax);
size_t NEvsToProcess = std::min(nmax, IH.GetNEvents());
double nmax_scaling = double(IH.GetNEvents()) / double(NEvsToProcess);
size_t NToShout = NEvsToProcess / 10;
IEventProcessor_INFO("Sample " << std::quoted(Name()) << " processing "
<< NEvsToProcess << " events.");
IInputHandler::ev_index_t ev_idx = 0;
while (ev_idx < NEvsToProcess) {
+ if (NToShout && !(ev_idx % NToShout)) {
+ IEventProcessor_INFO("\t\t Done " << ev_idx << "/" << NEvsToProcess
+ << " events.");
+ }
nuis::event::FullEvent const &fev = IH.GetFullEvent(ev_idx);
ProcessEventFunction(fev, IH.GetEventWeight(ev_idx) * nmax_scaling);
ev_idx++;
}
IEventProcessor_INFO("Sample " << std::quoted(Name()) << " processed "
<< NEvsToProcess << " events.");
}
};
diff --git a/src/utility/HistogramUtility.hxx b/src/utility/HistogramUtility.hxx
index b3f8ef4..851d84c 100644
--- a/src/utility/HistogramUtility.hxx
+++ b/src/utility/HistogramUtility.hxx
@@ -1,519 +1,522 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef UTILITY_HISTOGRAMUTILITY_HXX_SEEN
#define UTILITY_HISTOGRAMUTILITY_HXX_SEEN
#include "utility/ROOTUtility.hxx"
#include "exception/exception.hxx"
#include "string_parsers/from_string.hxx"
#include "TAxis.h"
#include "TH1D.h"
#include "TH1F.h"
#include "TH2D.h"
#include "TH2F.h"
#include "TH2Poly.h"
#include
#include
#include
#include
namespace nuis {
namespace utility {
NEW_NUIS_EXCEPT(unimplemented_GetHistogram_method);
NEW_NUIS_EXCEPT(invalid_histogram_descriptor);
NEW_NUIS_EXCEPT(invalid_histogram_name);
NEW_NUIS_EXCEPT(failed_to_clone);
inline bool IsFlowBin(TAxis const &ax, Int_t bin_it) {
return ((bin_it <= 0) || (bin_it >= (ax.GetNbins() + 1)));
}
inline bool IsInHistogramRange(TAxis const &ax, double v) {
Int_t bin_it = ax.FindFixBin(v);
return !IsFlowBin(ax, bin_it);
}
template struct HType_traits {};
template <> struct HType_traits {
using type = TH1;
static size_t const NDim = 1;
using NumericT = double;
static std::string name() { return "TH1"; }
};
template <> struct HType_traits {
using type = TH1D;
static size_t const NDim = 1;
using NumericT = double;
static std::string name() { return "TH1D"; }
};
template <> struct HType_traits {
using type = TH1F;
static size_t const NDim = 1;
using NumericT = float;
static std::string name() { return "TH1F"; }
};
template <> struct HType_traits {
using type = TH2;
static size_t const NDim = 2;
using NumericT = double;
static std::string name() { return "TH2"; }
};
template <> struct HType_traits {
using type = TH2D;
static size_t const NDim = 2;
using NumericT = double;
static std::string name() { return "TH2D"; }
};
template <> struct HType_traits {
using type = TH2F;
static size_t const NDim = 2;
using NumericT = float;
static std::string name() { return "TH2F"; }
};
template <> struct HType_traits {
using type = TH2Poly;
static size_t const NDim = 0;
using NumericT = double;
static std::string name() { return "TH2Poly"; }
};
template struct HType_Helper {};
template <> struct HType_Helper<1, void> {
using type = TH1;
static size_t const NDim = HType_traits::NDim;
using NumericT = HType_traits::NumericT;
static std::string name() { return HType_traits::name(); }
};
template <> struct HType_Helper<1, double> {
using type = TH1D;
static size_t const NDim = HType_traits::NDim;
using NumericT = HType_traits::NumericT;
static std::string name() { return HType_traits::name(); }
};
template <> struct HType_Helper<1, float> {
using type = TH1F;
static size_t const NDim = HType_traits::NDim;
using NumericT = HType_traits::NumericT;
static std::string name() { return HType_traits::name(); }
};
template <> struct HType_Helper<2, void> {
using type = TH2;
static size_t const NDim = HType_traits::NDim;
using NumericT = HType_traits::NumericT;
static std::string name() { return HType_traits::name(); }
};
template <> struct HType_Helper<2, double> {
using type = TH2D;
static size_t const NDim = HType_traits::NDim;
using NumericT = HType_traits::NumericT;
static std::string name() { return HType_traits::name(); }
};
template <> struct HType_Helper<2, float> {
using type = TH2F;
static size_t const NDim = HType_traits::NDim;
using NumericT = HType_traits::NumericT;
static std::string name() { return HType_traits::name(); }
};
template struct TH_Helper {};
template
struct TH_Helper::NDim == 1>::type> {
static size_t const NDim = HType_traits::NDim;
using NumericT = typename HType_traits::NumericT;
static std::string name() { return HType_traits::name(); }
static Int_t NbinsIncludeFlow(HT const &h) {
return h.GetXaxis()->GetNbins() + 2;
}
static Int_t Nbins(HT const &h) { return h.GetXaxis()->GetNbins(); }
static bool IsFlowBin(HT const &h, Int_t bin_it) {
return nuis::utility::IsFlowBin(*h.GetXaxis(), bin_it);
}
- static void Fill(HT &h, std::array const &v,
- double w = 1) {
+ static void Fill(HT &h, std::array const &v, double w = 1) {
h.Fill(v[0], w);
}
static void Scale(HT &h, NumericT SF, char const *opt = "") {
h.Scale(SF, opt);
}
static double Integral(HT const &h, char const *opt = "") {
return h.Integral(opt);
}
static Int_t NbinsIncludeFlow(std::unique_ptr const &h) {
return NbinsIncludeFlow(*h);
}
static Int_t Nbins(std::unique_ptr const &h) { return Nbins(*h); }
static bool IsFlowBin(std::unique_ptr const &h, Int_t bin_it) {
return IsFlowBin(*h, bin_it);
}
- static void Fill(std::unique_ptr &h,
- std::array const &v, double w = 1) {
+ static void Fill(std::unique_ptr &h, std::array const &v,
+ double w = 1) {
Fill(*h, v, w);
}
static void Scale(std::unique_ptr &h, NumericT SF, char const *opt = "") {
Scale(*h, SF, opt);
}
static double Integral(std::unique_ptr const &h, char const *opt = "") {
return Integral(*h, opt);
}
};
template
struct TH_Helper::NDim == 2>::type> {
static size_t const NDim = HType_traits::NDim;
using NumericT = typename HType_traits::NumericT;
static std::string name() { return HType_traits::name(); }
// TH2 ***************************************************************
static Int_t NbinsIncludeFlow(HT const &h) {
return (h.GetXaxis()->GetNbins() + 2) * (h.GetYaxis()->GetNbins() + 2);
}
static Int_t Nbins(HT const &h) {
return (h.GetXaxis()->GetNbins()) * (h.GetYaxis()->GetNbins());
}
static bool IsFlowBin(HT const &h, Int_t xbin_it, Int_t ybin_it) {
return nuis::utility::IsFlowBin(*h.GetXaxis(), xbin_it) ||
nuis::utility::IsFlowBin(*h.GetYaxis(), ybin_it);
}
- static void Fill(HT &h, std::array const &v,
- double w = 1) {
+ static void Fill(HT &h, std::array const &v, double w = 1) {
h.Fill(v[0], v[1], w);
}
static void Scale(HT &h, NumericT SF, char const *opt = "") {
h.Scale(SF, opt);
}
static double Integral(HT const &h, char const *opt = "") {
return h.Integral(opt);
}
static Int_t NbinsIncludeFlow(std::unique_ptr const &h) {
return NbinsIncludeFlow(*h);
}
static Int_t Nbins(std::unique_ptr const &h) { return Nbins(*h); }
static bool IsFlowBin(std::unique_ptr const &h, Int_t xbin_it,
Int_t ybin_it) {
return IsFlowBin(*h, xbin_it, ybin_it);
}
- static void Fill(std::unique_ptr &h,
- std::array const &v, double w = 1) {
+ static void Fill(std::unique_ptr &h, std::array const &v,
+ double w = 1) {
Fill(*h, v, w);
}
static void Scale(std::unique_ptr &h, NumericT SF, char const *opt = "") {
Scale(*h, SF, opt);
}
static double Integral(std::unique_ptr const &h, char const *opt = "") {
return Integral(*h, opt);
}
};
template
struct TH_Helper<
HT, typename std::enable_if::value>::type> {
static size_t const NDim = 2;
using NumericT = typename HType_traits::NumericT;
static std::string name() { return HType_traits::name(); }
// TH2Poly ***************************************************************
static Int_t NbinsIncludeFlow(HT const &h) { return h.GetNumberOfBins() + 9; }
static Int_t Nbins(HT const &h) { return h.GetNumberOfBins(); }
static bool IsFlowBin(HT const &h, Int_t bin_it) { return (bin_it < 0); }
- static void Fill(HT &h, std::array const &v,
- double w = 1) {
+ static void Fill(HT &h, std::array const &v, double w = 1) {
h.Fill(v[0], v[1], w);
}
static void Scale(HT &h, NumericT SF, char const *opt = "") {
bool width = (std::string(opt).find("width") != std::string::npos);
size_t nbins = Nbins(h);
for (size_t bin_it = 0; bin_it < nbins; ++bin_it) {
double bin_area = 1;
if (width) {
TH2PolyBin *poly_bin =
dynamic_cast(h.GetBins()->At(bin_it));
bin_area = poly_bin->GetArea();
}
h.SetBinContent(bin_it + 1,
h.GetBinContent(bin_it + 1) * (SF / bin_area));
h.SetBinError(bin_it + 1, h.GetBinError(bin_it + 1) * (SF / bin_area));
}
}
static double Integral(HT &h, char const *opt = "") {
bool width = (std::string(opt).find("width") != std::string::npos);
size_t nbins = Nbins(h);
double integral = 0;
for (size_t bin_it = 0; bin_it < nbins; ++bin_it) {
double bin_area = 1;
if (width) {
TH2PolyBin *poly_bin =
dynamic_cast(h.GetBins()->At(bin_it));
bin_area = poly_bin->GetArea();
}
integral += h.GetBinContent(bin_it + 1) * bin_area;
}
return integral;
}
static Int_t NbinsIncludeFlow(std::unique_ptr const &h) {
return NbinsIncludeFlow(*h);
}
static Int_t Nbins(std::unique_ptr const &h) { return Nbins(*h); }
static bool IsFlowBin(std::unique_ptr const &h, Int_t bin_it) {
return IsFlowBin(*h, bin_it);
}
- static void Fill(std::unique_ptr &h,
- std::array const &v, double w = 1) {
+ static void Fill(std::unique_ptr &h, std::array const &v,
+ double w = 1) {
Fill(*h, v, w);
}
static void Scale(std::unique_ptr &h, NumericT SF, char const *opt = "") {
Scale(*h, SF, opt);
}
static double Integral(std::unique_ptr &h, char const *opt = "") {
return Integral(*h, opt);
}
};
template
void Clear(typename std::enable_if::NDim != 0, HT>::type &h) {
for (Int_t bin_it = 0; bin_it < TH_Helper::NbinsIncludeFlow(h);
++bin_it) {
h.SetBinContent(bin_it, 0);
h.SetBinError(bin_it, 0);
}
}
template
void Clear(
typename std::enable_if::value, HT>::type &h) {
h.ClearBinContents();
}
template
inline std::unique_ptr Clone(HT const &source, bool clear = false,
std::string const &clone_name = "") {
std::unique_ptr target(dynamic_cast(
source.Clone(clone_name.size() ? clone_name.c_str() : "")));
if (!target) {
throw failed_to_clone()
<< "[ERROR]: Failed to clone a " << TH_Helper::name() << ".";
}
target->SetDirectory(nullptr);
if (clear) {
Clear(*target);
}
return target;
}
template
inline std::unique_ptr Clone(std::unique_ptr const &source,
bool clear = false,
std::string const &clone_name = "") {
return Clone(*source, clear, clone_name);
}
template
-inline std::unique_ptr GetHistogramFromROOTFile(TFile_ptr &f,
- std::string const &hname) {
+inline std::unique_ptr
+GetHistogramFromROOTFile(TFile_ptr &f, std::string const &hname,
+ bool ThrowIfMissing = true) {
f->Get(hname.c_str());
HT *h = dynamic_cast(f->Get(hname.c_str()));
if (!h) {
- throw invalid_histogram_name()
- << "[ERROR]: Failed to get " << TH_Helper::name() << " named "
- << std::quoted(hname) << " from input file "
- << std::quoted(f->GetName());
+ if (ThrowIfMissing) {
+ throw invalid_histogram_name()
+ << "[ERROR]: Failed to get " << TH_Helper::name() << " named "
+ << std::quoted(hname) << " from input file "
+ << std::quoted(f->GetName());
+ } else {
+ return nullptr;
+ }
}
std::unique_ptr clone = Clone(*h);
return clone;
}
template
-inline std::unique_ptr GetHistogramFromROOTFile(std::string const &fname,
- std::string const &hname) {
+inline std::unique_ptr
+GetHistogramFromROOTFile(std::string const &fname, std::string const &hname,
+ bool ThrowIfMissing = true) {
TFile_ptr temp = CheckOpenTFile(fname, "READ");
- return GetHistogramFromROOTFile(temp, hname);
+ return GetHistogramFromROOTFile(temp, hname, ThrowIfMissing);
}
template
std::unique_ptr GetHistogram(std::string const &input_descriptor) {
std::vector split_descriptor =
fhicl::string_parsers::ParseToVect(input_descriptor, ";",
true, true);
if (split_descriptor.size() == 1) { // assume text
throw unimplemented_GetHistogram_method();
} else if (split_descriptor.size() == 2) {
return GetHistogramFromROOTFile(split_descriptor[0],
split_descriptor[1]);
} else {
throw invalid_histogram_descriptor()
<< "[ERROR]: Failed to parse histogram descriptor: "
<< std::quoted(input_descriptor)
<< " as an input histogram (Text/ROOT).";
}
}
NEW_NUIS_EXCEPT(invalid_TH2Poly);
NEW_NUIS_EXCEPT(invalid_PolyBinSpecifierList);
// List of bins to put in a row, X/Y
struct PolyBinSpecifier {
double X, Y;
bool UseXAxis;
};
constexpr PolyBinSpecifier XPolyBinSpec(double X, double Y) {
return PolyBinSpecifier{X + std::numeric_limits::epsilon() * 1E2,
Y + std::numeric_limits::epsilon() * 1E2,
true};
}
constexpr PolyBinSpecifier YPolyBinSpec(double X, double Y) {
return PolyBinSpecifier{X + std::numeric_limits::epsilon() * 1E2,
Y + std::numeric_limits::epsilon() * 1E2,
false};
}
std::vector> GetTH2PolySlices(
std::unique_ptr &hinp,
std::vector> const &BinsSpecifiers) {
std::vector> slices;
size_t sl_it = 0;
for (auto &slice_spec : BinsSpecifiers) {
std::vector Binning;
std::vector BinContent;
std::vector BinError;
bool UseXAxis = false;
size_t bin_ctr = 0;
for (auto poly_bin_spec : slice_spec) {
Int_t bin_it = hinp->FindBin(poly_bin_spec.X, poly_bin_spec.Y);
if (bin_it < 1) {
std::cout << "[WARN]: When searching for matching bin: { X: "
<< poly_bin_spec.X << ", Y: " << poly_bin_spec.Y
<< "} got flow bin: " << bin_it << std::endl;
continue;
}
TH2PolyBin *poly_bin =
dynamic_cast(hinp->GetBins()->At(bin_it - 1));
if (!bin_ctr) {
UseXAxis = poly_bin_spec.UseXAxis;
} else if (UseXAxis != poly_bin_spec.UseXAxis) {
throw invalid_PolyBinSpecifierList()
<< "[ERROR]: For slice: " << sl_it
<< " of TH2Poly: " << std::quoted(hinp->GetName())
<< " bin specifier: " << bin_ctr << " was set to use the "
<< (poly_bin_spec.UseXAxis ? "X" : "Y")
<< " axis as the dependent axis of the slice, but previous bins "
"were set to use the "
<< (poly_bin_spec.UseXAxis ? "X" : "Y") << " axis.";
}
if ((poly_bin_spec.X < poly_bin->GetXMin()) ||
(poly_bin_spec.X >= poly_bin->GetXMax()) ||
(poly_bin_spec.Y < poly_bin->GetYMin()) ||
(poly_bin_spec.Y >= poly_bin->GetYMax())) {
std::cout
<< "[WARN]: Found bin doesn't seem to contain expected point: { X: "
<< poly_bin_spec.X << ", Y: " << poly_bin_spec.Y
<< "}, got bin_it = " << bin_it
<< " which had x_low: " << poly_bin->GetXMin()
<< ", and x_up: " << poly_bin->GetXMax()
<< ", y_low: " << poly_bin->GetYMin()
<< ", y_up: " << poly_bin->GetYMax() << std::endl;
}
double low = UseXAxis ? poly_bin->GetXMin() : poly_bin->GetYMin();
double up = UseXAxis ? poly_bin->GetXMax() : poly_bin->GetYMax();
if (!Binning.size()) { // Add low edge
Binning.push_back(low);
} else if (std::abs(Binning.back() - low) >
(std::numeric_limits::epsilon() * 1E2)) {
BinContent.push_back(0);
BinError.push_back(0);
Binning.push_back(low);
}
BinContent.push_back(hinp->GetBinContent(bin_it));
BinError.push_back(hinp->GetBinError(bin_it));
Binning.push_back(up);
}
slices.emplace_back(new TH1D(
(std::string(hinp->GetName()) + "_slice" + std::to_string(sl_it++))
.c_str(),
(std::string(";") +
(UseXAxis ? hinp->GetXaxis() : hinp->GetYaxis())->GetTitle() + ";" +
hinp->GetZaxis()->GetTitle())
.c_str(),
Binning.size() - 1, Binning.data()));
for (size_t bin_it = 0; bin_it < BinContent.size(); ++bin_it) {
slices.back()->SetBinContent(bin_it + 1, BinContent[bin_it]);
slices.back()->SetBinError(bin_it + 1, BinError[bin_it]);
}
}
return slices;
}
} // namespace utility
} // namespace nuis
#endif