diff --git a/data/nuA/Nuclear/T2K/CC0Pi/datResults.root b/data/nuA/Nuclear/T2K/CC0Pi/datResults.root
new file mode 100644
index 0000000..694e5d6
Binary files /dev/null and b/data/nuA/Nuclear/T2K/CC0Pi/datResults.root differ
diff --git a/data/nuA/Nuclear/T2K/CC0Pi/dphitResults.root b/data/nuA/Nuclear/T2K/CC0Pi/dphitResults.root
new file mode 100644
index 0000000..d17428b
Binary files /dev/null and b/data/nuA/Nuclear/T2K/CC0Pi/dphitResults.root differ
diff --git a/data/nuA/Nuclear/T2K/CC0Pi/dptResults.root b/data/nuA/Nuclear/T2K/CC0Pi/dptResults.root
new file mode 100644
index 0000000..8d51874
Binary files /dev/null and b/data/nuA/Nuclear/T2K/CC0Pi/dptResults.root differ
diff --git a/scripts/to_configure/SimpleMCStudy.cxx.in b/scripts/to_configure/SimpleMCStudy.cxx.in
index e52906a..d8bacdd 100644
--- a/scripts/to_configure/SimpleMCStudy.cxx.in
+++ b/scripts/to_configure/SimpleMCStudy.cxx.in
@@ -1,117 +1,121 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
//********************************************************************
#include "samples/SimpleMCStudy.hxx"
#include "utility/EventTopologyUtility.hxx"
#include "utility/FullEventUtility.hxx"
#include "utility/KinematicUtility.hxx"
#include "utility/PDGCodeUtility.hxx"
using namespace nuis::event;
using namespace nuis::utility;
class __SAMPLE_NAME__ : public SimpleMCStudy {
//! Add any plots that you want to fill as data members of the
//! comparisons class.
// std::unique_ptr myPlot;
public:
__SAMPLE_NAME__() { ReadGlobalConfigDefaults(); }
void Initialize(fhicl::ParameterSet const &instance_sample_configuration) {
- //! Set the verbosity of the sample logging macros.
+ //! This is also done by SimpleMCStudy::Initialize so can be removed here,
+ //! but left to demonstrate the simpel logging features. Set the verbosity
+ //! of the sample logging macros.
if (instance_sample_configuration.has_key("verbosity")) {
SetSampleVerbosity(
instance_sample_configuration.get("verbosity"));
//! Use these macros within member functions to log to the terminal
//! See src/samples/IEventProcessor.hxx for more details.
IEventProcessor_INFO(
"Verbosity set: "
<< instance_sample_configuration.get("verbosity"));
}
- // Get the global configuration for this sample, if it exists.
+ //! If your sample doesn't need any global configuration, this can be
+ //! removed.
+ //! Get the global configuration for this sample, if it exists.
fhicl::ParameterSet const &global_sample_configuration =
nuis::config::GetDocument().get(
std::string("global.sample_configuration.") + Name(),
fhicl::ParameterSet());
//! Instantiate any plots that you want to fill
// myPlot = std::unique_ptr(
// new TH1D("myPlot",
// "title;xlabel;ylabel",
// nbinsx,xlow,xhigh));
// Perform any per-sample configuration in the base class
SimpleMCStudy::Initialize(instance_sample_configuration);
//! Here you do whatever you want to do on a per-event basis.
//! This method will be run for every event given by the input handler, the
//! weight allows differential xsecs to be plotted easily and also contains
//! any requested reweightable parameter variations.
ProcessEventFunction = [&](nuis::event::FullEvent const &ev,
double weight) -> void {
//! Select events
//! See src/utility/EventTopologyUtility.hxx for more pre-defined
//! topological signals.
// if (!IsCC0Pi(fev)) {
// return false;
// }
//! Build composite variables
//! See src/event/FullEvent.hxx for the full event class definition.
//! See src/utility/FullEventUtility.hxx for more helper methods for
//! interacting with the event class.
// Particle FSMuPlus = GetHMFSParticle(fev, {pdgcodes::kMuPlus});
//! Fill your histograms!
// myPlot->Fill(FSMuPlus.CosTheta(), weight);
};
}
std::string Name() { return "__SAMPLE_NAME__"; }
//! Here you can write any custom histograms to TTrees that your sample has
//! been handling.
void Write() {
//! The nuis::persistency helper methods will help write your objects to the
//! correct place and in a TFile that will be closed before the program
//! exits.
//! They will write a copy of the TObject that is passed, so you do not need
//! to release/worry about double freeing any histograms passed.
//! write_directory is set up by the baseclass and is configurable in the
//! instance configuration.
// nuis::persistency::WriteToOutputFile(Extra_plot,
// Extra_plot->GetName(),
// write_directory);
//! It probably isn't a good idea to write a TTree with the above method,
//! instead, try using nuis::persistency::GetOutputFile() to get the TFile
//! and then manually writing/giving ownership of your TTree to the output
//! file.
}
};
//! These declarations allow your class to be loaded dynamically by NUISANCE
DECLARE_PLUGIN(IEventProcessor, __SAMPLE_NAME__);
diff --git a/src/app/nuissamples.cxx b/src/app/nuissamples.cxx
index d3beeb5..809f275 100644
--- a/src/app/nuissamples.cxx
+++ b/src/app/nuissamples.cxx
@@ -1,115 +1,155 @@
#include "config/GlobalConfiguration.hxx"
#include "input/IInputHandler.hxx"
#include "plugins/Instantiate.hxx"
#include "samples/IDataComparison.hxx"
#include "utility/StringUtility.hxx"
#include
#include
#include
#include
-std::string search_term = ".*";
-bool fStrictRegex = false;
+std::string name_search_term = ".*";
+bool strict_name_regex = false;
+
+std::string target_search_term = ".*";
+bool strict_target_regex = false;
+
+std::string year_search_term = ".*";
+
std::string config_out_filename = "";
+bool NameOnly = false;
+
void SayUsage(char const *argv[]) {
- std::cout << "[USAGE]: " << argv[0]
- << "\n"
- "\t-s : Used to filter known IDataComparisons. "
- "\n\t\tWild cards are added to either side of the search term.\n"
- "\t-S : Used to filter known IDataComparisons. "
- "\n\t\tThe exact passed term is used.\n"
- "\t-o : Dump example sample configuration file\n"
- "\t\tfor matching samples.\n"
- << std::endl;
+ std::stringstream ss("");
+ ss << "[USAGE]: " << argv[0]
+ << "\n"
+ "\t-n,-N,-t,-T,-y : Filters known IDataComparisons by "
+ "the search term (-n: Name, -t: Target, -y: Year). Capitalized "
+ "versions filter on exact match."
+ "\t-o : Dump example sample configuration "
+ "file for matching samples.\n"
+ "\t--name-only : Only write out matching sample "
+ "names. (Still applies all search terms)\n";
+
+ std::cout << nuis::utility::indent_apply_width(ss.str()) << std::endl;
}
void handleOpts(int argc, char const *argv[]) {
int opt = 1;
while (opt < argc) {
if ((std::string(argv[opt]) == "-?") ||
(std::string(argv[opt]) == "--help")) {
SayUsage(argv);
exit(0);
} else if (std::string(argv[opt]) == "-s") {
- search_term = argv[++opt];
+ name_search_term = argv[++opt];
} else if (std::string(argv[opt]) == "-S") {
- fStrictRegex = true;
- search_term = argv[++opt];
+ strict_name_regex = true;
+ name_search_term = argv[++opt];
+ } else if (std::string(argv[opt]) == "-t") {
+ target_search_term = argv[++opt];
+ } else if (std::string(argv[opt]) == "-T") {
+ strict_target_regex = true;
+ target_search_term = argv[++opt];
+ } else if (std::string(argv[opt]) == "-y") {
+ year_search_term = argv[++opt];
} else if (std::string(argv[opt]) == "-o") {
config_out_filename = argv[++opt];
+ } else if (std::string(argv[opt]) == "--name-only") {
+ NameOnly = true;
} else {
std::cout << "[ERROR]: Unknown option: " << argv[opt] << std::endl;
SayUsage(argv);
exit(1);
}
opt++;
}
}
NEW_NUIS_EXCEPT(invalid_output_file);
int main(int argc, char const *argv[]) {
nuis::config::EnsureConfigurationRead("nuis.global.config.fcl");
nuis::config::EnsureConfigurationRead("nuis.datacomparisons.fcl");
handleOpts(argc, argv);
- std::regex rpattern(fStrictRegex ? search_term
- : std::string(".*") + search_term + ".*");
+ std::regex rpattern_name(strict_name_regex
+ ? name_search_term
+ : std::string(".*") + name_search_term + ".*");
+ std::regex rpattern_target(
+ strict_target_regex ? target_search_term
+ : std::string(".*") + target_search_term + ".*");
+
+ std::regex rpattern_year(std::string(".*") + year_search_term + ".*");
std::vector example_sample_configs;
for (std::string const &comparison_set_key :
nuis::config::GetDocument()
.get("data_comparisons")
.get_names()) {
for (std::string const &sample_name :
nuis::config::GetDocument().get>(
std::string("data_comparisons.") + comparison_set_key)) {
- if (!std::regex_match(sample_name, rpattern)) {
+ if (!std::regex_match(sample_name, rpattern_name)) {
continue;
}
nuis::plugins::plugin_traits::unique_ptr_t sample =
nuis::plugins::Instantiate(sample_name);
+ if (!std::regex_match(sample->GetTargetMaterial(), rpattern_target)) {
+ continue;
+ }
+
+ if (!std::regex_match(sample->GetYear(), rpattern_year)) {
+ continue;
+ }
+
std::cout << sample->Name() << std::endl;
- std::cout << "\tJournal: " << sample->GetJournalReference() << std::endl;
- std::cout << "\tTarget: " << sample->GetTargetMaterial() << std::endl;
- std::cout << "\tFlux: " << sample->GetFluxDescription() << std::endl;
- std::cout << "\tSignal: " << sample->GetSignalDescription() << std::endl;
- std::cout << "\tDocs: \n"
- << nuis::utility::indent_apply_width(sample->GetDocumentation(),
- 10)
- << std::endl;
- std::cout << "\tExample_Config: {\n"
- << nuis::utility::indent_apply_width(
- sample->GetExampleConfiguration().to_indented_string(),
- 12)
- << "\n\t}\n"
- << std::endl;
- ;
+ if (!NameOnly) {
+ std::cout << "\tJournal: " << sample->GetJournalReference()
+ << std::endl;
+ std::cout << "\tDOI: " << sample->GetDOI() << std::endl;
+ std::cout << "\tYear: " << sample->GetYear() << std::endl;
+ std::cout << "\tTarget: " << sample->GetTargetMaterial() << std::endl;
+ std::cout << "\tFlux: " << sample->GetFluxDescription() << std::endl;
+ std::cout << "\tSignal: " << sample->GetSignalDescription()
+ << std::endl;
+ std::cout << "\tDocs: \n"
+ << nuis::utility::indent_apply_width(
+ sample->GetDocumentation(), 10)
+ << std::endl;
+ std::cout << "\tExample_Config: {\n"
+ << nuis::utility::indent_apply_width(
+ sample->GetExampleConfiguration().to_indented_string(),
+ 12)
+ << "\n\t}\n"
+ << std::endl;
+ }
+
example_sample_configs.push_back(sample->GetExampleConfiguration());
}
}
if (config_out_filename.length()) {
std::ofstream out_file(config_out_filename);
if (!out_file) {
throw invalid_output_file() << "[ERROR]: Failed to open output file: "
<< std::quoted(config_out_filename);
}
fhicl::ParameterSet example_config;
example_config.put("samples", example_sample_configs);
out_file << example_config.to_indented_string();
}
}
diff --git a/src/event/types.hxx b/src/event/types.hxx
index f40e5f0..8850d33 100644
--- a/src/event/types.hxx
+++ b/src/event/types.hxx
@@ -1,125 +1,129 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef EVENT_TYPES_HXX_SEEN
#define EVENT_TYPES_HXX_SEEN
#include "exception/exception.hxx"
namespace nuis {
namespace event {
// Modelled off the NEUT interaction codes.
#define NUIS_INTERACTION_CHANNEL_LIST \
X(kCCQE, 1) \
X(kCC2p2h, 2) \
X(kCCSPP_PPip, 11) \
X(kCCSPP_PPi0, 12) \
X(kCCSPP_NPip, 13) \
+ X(kCCDiffPP, 15) \
X(kCCCohPi, 16) \
X(kCCResGamma, 17) \
X(kCCTransitionMPi, 21) \
X(kCCResEta0, 22) \
X(kCCResK, 23) \
X(kCCDIS, 26) \
\
X(kNCSPP_NPi0, 31) \
X(kNCSPP_PPi0, 32) \
X(kNCSPP_PPim, 33) \
X(kNCSPP_NPip, 34) \
+ X(kNCDiffPP, 35) \
X(kNCCohPi, 36) \
X(kNCResNGamma, 38) \
X(kNCResPGamma, 39) \
X(kNCTransitionMPi, 41) \
X(kNCResNEta0, 42) \
X(kNCResPEta0, 43) \
X(kNCResK0, 44) \
X(kNCResKp, 45) \
X(kNCDIS, 46) \
X(kNCELP, 51) \
X(kNCELN, 52) \
X(kNC2p2h, 53) \
\
X(kCCQE_nub, -1) \
X(kCC2p2h_nub, -2) \
X(kCCSPP_NPim_nub, -11) \
X(kCCSPP_NPi0_nub, -12) \
X(kCCSPP_PPim_nub, -13) \
+ X(kCCDiffPP_nub, -15) \
X(kCCCohPi_nub, -16) \
X(kCCResGamma_nub, -17) \
X(kCCTransitionMPi_nub, -21) \
X(kCCResEta0_nub, -22) \
X(kCCResK_nub, -23) \
X(kCCDIS_nub, -26) \
\
X(kNCSPP_NPi0_nub, -31) \
X(kNCSPP_PPi0_nub, -32) \
X(kNCSPP_PPim_nub, -33) \
X(kNCSPP_NPip_nub, -34) \
+ X(kNCDiffPP_nub, -35) \
X(kNCCohPi_nub, -36) \
X(kNCResNGamma_nub, -38) \
X(kNCResPGamma_nub, -39) \
X(kNCTransitionMPi_nub, -41) \
X(kNCResNEta0_nub, -42) \
X(kNCResPEta0_nub, -43) \
X(kNCResK0_nub, -44) \
X(kNCResKp_nub, -45) \
X(kNCDIS_nub, -46) \
X(kNCELP_nub, -51) \
X(kNCELN_nub, -52) \
X(kNC2p2h_nub, -53) \
\
X(kUndefined, 0)
#define X(A, B) A = B,
enum class Channel_t { NUIS_INTERACTION_CHANNEL_LIST };
#undef X
#define X(A, B) \
case B: { \
return Channel_t::A; \
}
inline Channel_t FromNEUTCode(int nc) {
switch (nc) {
NUIS_INTERACTION_CHANNEL_LIST
default: { return Channel_t::kUndefined; }
}
}
#undef X
using PDG_t = long;
inline bool Is2p2h(Channel_t chan) {
return ((chan == Channel_t::kCC2p2h) || (chan == Channel_t::kNC2p2h) ||
(chan == Channel_t::kCC2p2h_nub) || (chan == Channel_t::kNC2p2h_nub));
}
} // namespace event
} // namespace nuis
#define X(A, B) \
case nuis::event::Channel_t::A: { \
return os << #A; \
}
inline std::ostream &operator<<(std::ostream &os, nuis::event::Channel_t te) {
switch (te) { NUIS_INTERACTION_CHANNEL_LIST }
return os;
}
#undef X
#endif
diff --git a/src/generator/input/NEUTInputHandler.cxx b/src/generator/input/NEUTInputHandler.cxx
index dbd85bd..52cadce 100644
--- a/src/generator/input/NEUTInputHandler.cxx
+++ b/src/generator/input/NEUTInputHandler.cxx
@@ -1,198 +1,250 @@
#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", false);
+ std::string flux_name = ps.get("flux_hist_name", "flux_numu");
+ std::string evtrt_name = ps.get("evtrt_hist_name", "evtrt_numu");
+
+ std::string override_flux_file =
+ ps.get("override_flux_file", "");
+
+ if (override_flux_file.size()) {
+ fFlux = GetHistogramFromROOTFile(override_flux_file, flux_name);
+ RebuildEventRate();
+ } else {
+ fFlux =
+ GetHistogramFromROOTFile(fInputTreeFile.file, flux_name, false);
+ }
+
if (fFlux) {
- fEvtrt = GetHistogramFromROOTFile(fInputTreeFile.file, "evtrt_numu");
+ if (!fEvtrt) {
+ fEvtrt = GetHistogramFromROOTFile(fInputTreeFile.file, evtrt_name);
+ }
- fFileWeight =
- fEvtrt->Integral() / (fFlux->Integral() * double(GetNEvents()));
+ 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;
+ 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;
+ return (xsec / count) * 1E-38;
+}
+
+void NEUTInputHandler::RebuildEventRate() {
+ 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();
+
+ XSec->Fill(E, fNeutVect->Totcrs);
+ Entry->Fill(E);
+ }
+ std::cout << std::endl << "[INFO]: Done!" << std::endl;
+
+ // Binned total xsec
+ XSec->Divide(Entry.get());
+
+ fEvtrt = Clone(XSec, false, "evtrt");
+ // Event rate prediction
+ fEvtrt->Multiply(fFlux.get());
}
double
NEUTInputHandler::GetXSecScaleFactor(std::pair const &) const {
throw input_handler_feature_unimplemented()
<< "[ERROR]: Flux cuts not yet implemented for NEUT input handler.";
}
NeutVect *NEUTInputHandler::GetNeutEvent(ev_index_t idx) const {
if (idx >= GetNEvents()) {
throw IInputHandler::invalid_entry()
<< "[ERROR]: Attempted to get entry " << idx
<< " from an InputHandler with only " << GetNEvents();
}
fInputTreeFile.tree->GetEntry(idx);
return fNeutVect;
}
MinimalEvent const &NEUTInputHandler::GetMinimalEvent(ev_index_t idx) const {
if (idx >= GetNEvents()) {
throw IInputHandler::invalid_entry()
<< "[ERROR]: Attempted to get entry " << idx
<< " from an InputHandler with only " << GetNEvents();
}
fInputTreeFile.tree->GetEntry(idx);
fReaderEvent.mode = IntToChannel(fNeutVect->Mode);
size_t NPart = fNeutVect->Npart();
for (size_t part_it = 0; part_it < NPart; part_it++) {
NeutPart const &part = (*fNeutVect->PartInfo(part_it));
if ((part.fIsAlive == false) && (part.fStatus == -1) &&
IsNeutralLepton(part.fPID)) {
fReaderEvent.probe_E = part.fP.T();
fReaderEvent.probe_pdg = part.fPID;
break;
}
}
fReaderEvent.XSecWeight = fFileWeight;
if (fWeightCache.size() <= idx) {
fWeightCache.push_back(fReaderEvent.XSecWeight);
}
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() {
return GeneratorManager::Get().EnsureGeneratorRegistered("NEUT");
}
DECLARE_PLUGIN(IInputHandler, NEUTInputHandler);
diff --git a/src/generator/input/NEUTInputHandler.hxx b/src/generator/input/NEUTInputHandler.hxx
index 7fc1c3f..36050cb 100644
--- a/src/generator/input/NEUTInputHandler.hxx
+++ b/src/generator/input/NEUTInputHandler.hxx
@@ -1,69 +1,70 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#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();
+ void RebuildEventRate();
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/persistency/ROOTOutput.cxx b/src/persistency/ROOTOutput.cxx
index 8103ee8..619a820 100644
--- a/src/persistency/ROOTOutput.cxx
+++ b/src/persistency/ROOTOutput.cxx
@@ -1,64 +1,81 @@
#include "persistency/ROOTOutput.hxx"
#include "utility/ROOTUtility.hxx"
#include "config/GlobalConfiguration.hxx"
#include "fhiclcpp/ParameterSet.h"
#include "TFile.h"
namespace nuis {
namespace persistency {
struct NamedTFile {
std::string name;
std::unique_ptr file;
NamedTFile() : name(""), file(nullptr) {}
NamedTFile(NamedTFile &&other)
: name(std::move(other.name)), file(std::move(other.file)) {}
~NamedTFile() {
if (file) {
file->Write();
file->Close();
}
}
};
static std::vector Files;
std::unique_ptr &GetOutputFile(std::string const &name) {
for (NamedTFile &file : Files) {
if (file.name == name) {
return file.file;
}
}
fhicl::ParameterSet const &persistency =
config::GetDocument().get("persistency");
std::string file_name = persistency.get(name + ".output_file");
std::string open_opts =
persistency.get(name + ".open_mode", "CREATE");
NamedTFile ntf;
ntf.name = name;
ntf.file = std::make_unique(file_name.c_str(), open_opts.c_str());
if (!ntf.file || !ntf.file->IsOpen()) {
- throw utility::failed_to_open_TFile()
- << "[ERROR]: Failed to open output file: " << std::quoted(file_name)
- << " in write mode with opts = " << std::quoted(open_opts);
+ if ((name == "default") && (open_opts == "CREATE")) {
+ std::cout
+ << "[WARN]: It appears the default file cannot be opened because it "
+ "exists already, to stop you wasting processing time, the default "
+ "output stream will be written to nuis.default.tmp.root in the "
+ "current director in RECREATE mode."
+ << std::endl;
+ file_name = "nuis.default.tmp.root";
+ open_opts = "RECREATE";
+ ntf.file = std::make_unique(file_name.c_str(), open_opts.c_str());
+ if (!ntf.file || !ntf.file->IsOpen()) {
+ throw utility::failed_to_open_TFile()
+ << "[ERROR]: Failed to open output file: " << std::quoted(file_name)
+ << " in write mode with opts = " << std::quoted(open_opts);
+ }
+ } else {
+ throw utility::failed_to_open_TFile()
+ << "[ERROR]: Failed to open output file: " << std::quoted(file_name)
+ << " in write mode with opts = " << std::quoted(open_opts);
+ }
}
Files.push_back(std::move(ntf));
return Files.back().file;
}
void CloseOpenTFiles() {
for (NamedTFile &f : Files) {
std::cout << "[INFO]: Closing open TFile: " << f.name << " "
<< f.file->GetName() << std::endl;
}
Files.clear();
}
} // namespace persistency
} // namespace nuis
diff --git a/src/samples/IDataComparison.hxx b/src/samples/IDataComparison.hxx
index 1fdc67d..b7cd809 100644
--- a/src/samples/IDataComparison.hxx
+++ b/src/samples/IDataComparison.hxx
@@ -1,73 +1,84 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#ifndef SAMPLES_IDATACOMPARISON_HXX_SEEN
#define SAMPLES_IDATACOMPARISON_HXX_SEEN
#include "samples/IEventProcessor.hxx"
#include "fhiclcpp/ParameterSet.h"
#include
#include
class IDataComparison : public IEventProcessor {
public:
virtual double GetGOF() = 0;
virtual double GetNDOGuess() = 0;
virtual std::string GetJournalReference() {
std::stringstream ss("");
ss << "Unknown Journal Ref. for IDataComparison: " << std::quoted(Name());
return ss.str();
}
-
+ virtual std::string GetDOI() {
+ std::stringstream ss("");
+ ss << "Unknown DOI for IDataComparison: " << std::quoted(Name());
+ return ss.str();
+ }
+ virtual std::string GetYear() {
+ std::stringstream ss("");
+ ss << "Unknown Year for IDataComparison: " << std::quoted(Name());
+ return ss.str();
+ }
virtual std::string GetTargetMaterial() {
std::stringstream ss("");
- ss << "Unknown Target material for IDataComparison: " << std::quoted(Name());
+ ss << "Unknown Target material for IDataComparison: "
+ << std::quoted(Name());
return ss.str();
}
-
virtual std::string GetFluxDescription() {
std::stringstream ss("");
- ss << "Unknown Flux description for IDataComparison: " << std::quoted(Name());
+ ss << "Unknown Flux description for IDataComparison: "
+ << std::quoted(Name());
return ss.str();
}
virtual std::string GetSignalDescription() {
std::stringstream ss("");
- ss << "Unknown Signal description for IDataComparison: " << std::quoted(
- Name());
+ ss << "Unknown Signal description for IDataComparison: "
+ << std::quoted(Name());
return ss.str();
}
virtual std::string GetDocumentation() {
std::stringstream ss("");
- ss << "No documentation provided for IDataComparison: " << std::quoted(Name());
+ ss << "No documentation provided for IDataComparison: "
+ << std::quoted(Name());
return ss.str();
}
virtual fhicl::ParameterSet GetExampleConfiguration() {
return fhicl::ParameterSet();
}
};
DECLARE_PLUGIN_INTERFACE(IDataComparison);
#endif
diff --git a/src/samples/SimpleDataComparison.hxx b/src/samples/SimpleDataComparison.hxx
index a63e4df..b3a5e17 100644
--- a/src/samples/SimpleDataComparison.hxx
+++ b/src/samples/SimpleDataComparison.hxx
@@ -1,426 +1,435 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
#pragma once
#include "samples/IDataComparison.hxx"
#include "event/FullEvent.hxx"
#include "input/InputManager.hxx"
#include "persistency/ROOTOutput.hxx"
#include "utility/FileSystemUtility.hxx"
#include "utility/HistogramUtility.hxx"
#include "utility/KinematicUtility.hxx"
#include "utility/ROOTUtility.hxx"
#include "utility/StatsUtility.hxx"
#include
#include
#include
#include
#include
template ::type>
class SimpleDataComparison : public IDataComparison {
NEW_NUIS_EXCEPT(invalid_SimpleDataComparison_initialization);
NEW_NUIS_EXCEPT(SimpleDataComparison_already_finalized);
protected:
using HistType = HT;
using TH_Help = typename nuis::utility::TH_Helper;
static size_t const NDim = nd;
nuis::input::InputManager::Input_id_t fIH_id;
std::string write_directory;
size_t NMaxSample_override;
int fIsShapeOnly;
int fIsFluxUnfolded;
std::vector fSignalCache;
std::vector> fProjectionCache;
std::string fDataInputDescriptor;
std::unique_ptr fData;
std::string fMaskInputDescriptor;
std::unique_ptr fMask;
std::string fCovarianceInputDescriptor;
std::unique_ptr fCovariance;
std::unique_ptr fPrediction;
std::unique_ptr fPrediction_xsec;
std::unique_ptr fPrediction_shape;
std::unique_ptr fPrediction_comparison;
bool fComparisonFinalized;
std::string fJournalReference;
+ std::string fDOI;
+ std::string fYear;
std::string fTargetMaterial;
std::string fFluxDescription;
std::string fSignalDescription;
nuis::utility::ENuRange energy_cut;
std::function IsSigFunc;
std::function(nuis::event::FullEvent const &)>
CompProjFunc;
// If assigned by subclass will be called on for all events, bool signifies
// whether the event was selected.
std::function
ProcessExtraFunc;
public:
SimpleDataComparison() {
fIH_id = std::numeric_limits::max();
write_directory = "";
NMaxSample_override = std::numeric_limits::max();
fDataInputDescriptor = "";
fData = nullptr;
fMaskInputDescriptor = "";
fMask = nullptr;
fCovarianceInputDescriptor = "";
fCovariance = nullptr;
fPrediction = nullptr;
fPrediction_xsec = nullptr;
fPrediction_shape = nullptr;
fPrediction_comparison = nullptr;
fComparisonFinalized = false;
IsSigFunc = [](nuis::event::FullEvent const &) -> bool { return true; };
CompProjFunc =
[](nuis::event::FullEvent const &) -> std::array {
std::array arr;
for (NumericT &el : arr) {
el = 0;
}
return arr;
};
ProcessExtraFunc =
std::function();
fJournalReference = "";
+ fDOI = "";
+ fYear = "";
fTargetMaterial = "";
fFluxDescription = "";
fSignalDescription = "";
fIsShapeOnly = -1;
fIsFluxUnfolded = -1;
energy_cut = nuis::utility::ENuRange{std::numeric_limits::max(),
std::numeric_limits::max()};
}
fhicl::ParameterSet fGlobalConfig;
fhicl::ParameterSet fInstanceConfig;
void ReadGlobalConfigDefaults() {
fGlobalConfig = nuis::config::GetDocument().get(
std::string("global.sample_configuration.") + Name(),
fhicl::ParameterSet());
if (!fJournalReference.length()) {
fJournalReference = fGlobalConfig.get("journal_reference",
- fJournalReference);
+ "Not yet published");
}
+ if (!fDOI.length()) {
+ fDOI = fGlobalConfig.get("DOI", "Unknown DOI");
+ }
+ if (!fYear.length()) {
+ fYear = fGlobalConfig.get("year", "Unknown year");
+ }
+
if (!fTargetMaterial.length()) {
- fTargetMaterial =
- fGlobalConfig.get("target_material", fTargetMaterial);
+ fTargetMaterial = fGlobalConfig.get(
+ "target_material", "Unknown target material");
}
if (!fFluxDescription.length()) {
fFluxDescription =
- fGlobalConfig.get("flux_description", fFluxDescription);
+ fGlobalConfig.get("flux_description", "Unknown flux");
}
if (!fSignalDescription.length()) {
fSignalDescription = fGlobalConfig.get("signal_description",
- fSignalDescription);
+ "Unknown Signal");
}
if (fIsShapeOnly == -1) {
fIsShapeOnly = fGlobalConfig.get("shape_only", false);
}
if (fIsFluxUnfolded == -1) {
fIsFluxUnfolded = fGlobalConfig.get("flux_unfolded", false);
}
if ((energy_cut.first == std::numeric_limits::max()) &&
(fGlobalConfig.has_key("enu_range"))) {
energy_cut = fGlobalConfig.get>("enu_range");
}
}
virtual std::string GetJournalReference() { return fJournalReference; }
+ virtual std::string GetDOI() { return fDOI; }
+ virtual std::string GetYear() { return fYear; }
virtual std::string GetTargetMaterial() { return fTargetMaterial; }
virtual std::string GetFluxDescription() { return fFluxDescription; }
virtual std::string GetSignalDescription() { return fSignalDescription; }
void SetShapeOnly(bool iso) { fIsShapeOnly = iso; }
void SetFluxUnfolded(bool ifo) { fIsFluxUnfolded = ifo; }
void SetData(std::string const &data_descriptor) {
fDataInputDescriptor = data_descriptor;
}
void SetMask(std::string const &mask_descriptor) {
fMaskInputDescriptor = mask_descriptor;
}
void SetCovariance(std::string const &cov_descriptor) {
fCovarianceInputDescriptor = cov_descriptor;
}
virtual void FillProjection(std::array const &proj,
NumericT event_weight) {
TH_Help::Fill(fPrediction, proj, event_weight);
}
virtual void FinalizeComparison() {
if (fComparisonFinalized) {
throw SimpleDataComparison_already_finalized()
<< "[ERROR]: Attempted to re-finalize a comparison for "
"SimpleDataComparison: "
<< std::quoted(Name());
}
fPrediction_xsec =
nuis::utility::Clone(fPrediction, false, "Prediction_xsec");
IInputHandler const &IH =
nuis::input::InputManager::Get().GetInputHandler(fIH_id);
TH_Help::Scale(fPrediction_xsec, 1.0, "width");
// If we have a flux cut
if (energy_cut.first != std::numeric_limits::max()) {
TH_Help::Scale(fPrediction_xsec, IH.GetXSecScaleFactor(energy_cut));
}
fPrediction_shape =
nuis::utility::Clone(fPrediction_xsec, false, "Prediction_shape");
if (fData) {
TH_Help::Scale(fPrediction_shape,
TH_Help::Integral(fData, "width") /
TH_Help::Integral(fPrediction_shape, "width"));
} else {
IEventProcessor_WARN(
"When Finalizing comparison, no Data histogram available.");
}
if (fIsFluxUnfolded) {
// fPrediction_comparison
} else if (fIsShapeOnly) {
fPrediction_comparison = nuis::utility::Clone(fPrediction_shape, false,
"Prediction_comparison");
} else {
fPrediction_comparison = nuis::utility::Clone(fPrediction_xsec, false,
"Prediction_comparison");
}
fComparisonFinalized = true;
}
void Initialize(fhicl::ParameterSet const &instance_sample_configuration) {
fInstanceConfig = instance_sample_configuration;
- 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();
if (fInstanceConfig.has_key("fake_data")) {
fData = nuis::utility::GetHistogram(
fInstanceConfig.get("fake_data_histogram"));
} else if (!fGlobalConfig.get("has_data", true) ||
!fInstanceConfig.get("has_data", true)) {
// Explicitly not expecting data
} else {
if (!fDataInputDescriptor.length()) {
if (!fGlobalConfig.has_key("data_descriptor")) {
throw invalid_SimpleDataComparison_initialization()
<< "[ERROR]: SimpleDataComparison::Initialize for "
"IDataComparison: "
<< std::quoted(Name())
<< " failed as no input data was set by a call to "
"SimpleDataComparison::SetData and no data_descriptor for "
"this SimpleDataComparison could be found in the global "
"configuration.";
}
fDataInputDescriptor =
fGlobalConfig.get("data_descriptor");
}
fData = nuis::utility::GetHistogram(fDataInputDescriptor);
}
if (!fPrediction) {
if (!fData) {
throw invalid_SimpleDataComparison_initialization()
<< "[ERROR]: SimpleDataComparison::Initialize for "
"IDataComparison: "
<< std::quoted(Name())
<< " failed. As `has_data: false` was set in the configuration "
"(global: "
<< (!fGlobalConfig.get("has_data", true))
<< ", instance: " << !fInstanceConfig.get("has_data", true)
<< "), the instance constructor must supply the fPrediction "
"binning, and it wasn't.";
}
fPrediction = nuis::utility::Clone(fData, true, "Prediction");
}
if (fCovarianceInputDescriptor.length()) {
fCovariance =
nuis::utility::GetHistogram(fCovarianceInputDescriptor);
} else if (fGlobalConfig.has_key("covariance_descriptor")) {
fCovariance = nuis::utility::GetHistogram(
fGlobalConfig.get("covariance_descriptor"));
}
if (fMaskInputDescriptor.length()) {
fMask = nuis::utility::GetHistogram(fMaskInputDescriptor);
} else if (fGlobalConfig.has_key("mask_descriptor")) {
fMask = nuis::utility::GetHistogram(
fGlobalConfig.get("mask_descriptor"));
}
if (fInstanceConfig.has_key("verbosity")) {
SetSampleVerbosity(fInstanceConfig.get("verbosity"));
}
NMaxSample_override =
fInstanceConfig.get("nmax", std::numeric_limits::max());
write_directory =
fInstanceConfig.get("write_directory", Name());
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;
bool DetermineSignalEvents = !fSignalCache.size();
nuis::utility::Clear(*fPrediction);
fComparisonFinalized = false;
size_t cache_ctr = 0;
while (ev_idx < NEvsToProcess) {
if (DetermineSignalEvents) {
nuis::event::FullEvent const &fev = IH.GetFullEvent(ev_idx);
bool is_sig = IsSigFunc(fev);
fSignalCache.push_back(is_sig);
if (is_sig) {
fProjectionCache.push_back(CompProjFunc(fev));
}
if (ProcessExtraFunc) {
ProcessExtraFunc(fev, is_sig,
IH.GetEventWeight(ev_idx) * nmax_scaling);
}
}
if (NToShout && !(ev_idx % NToShout)) {
IEventProcessor_INFO("\tProcessed " << ev_idx << "/" << NEvsToProcess
<< " events.");
}
if (fSignalCache[ev_idx]) {
FillProjection(fProjectionCache[cache_ctr++],
IH.GetEventWeight(ev_idx) * nmax_scaling);
}
ev_idx++;
}
IEventProcessor_INFO("\t" << fProjectionCache.size() << "/" << NEvsToProcess
<< " events passed selection.");
}
void Write() {
if (!fComparisonFinalized) {
FinalizeComparison();
}
nuis::persistency::WriteToOutputFile(
fPrediction_comparison, "Prediction", write_directory);
nuis::persistency::WriteToOutputFile(
fPrediction_xsec, "Prediction_xsec", write_directory);
if (fData) {
nuis::persistency::WriteToOutputFile(fData, "Data",
write_directory);
nuis::persistency::WriteToOutputFile(
fPrediction_shape, "Prediction_shape", write_directory);
}
}
double GetGOF() {
if (!fComparisonFinalized) {
FinalizeComparison();
}
if (fData && fPrediction_comparison) {
return nuis::utility::GetChi2(fData, fPrediction_comparison);
} else
return std::numeric_limits::max();
}
double GetNDOGuess() {
if (fData) {
return TH_Help::Nbins(fData);
}
return 0;
}
fhicl::ParameterSet GetExampleConfiguration() {
fhicl::ParameterSet exps;
exps.put("name", Name());
- exps.put("input_type", "Generator");
- exps.put("file", "Events.root");
+ exps.put("input_type", "NuWro/NEUT/GENIE");
+ exps.put("file", "input_events.root");
exps.put("write_directory", Name());
- exps.put("fake_data_histogram", "/path/to/file.root;h_name");
-
return exps;
}
};
typedef SimpleDataComparison<1, double, TH1D> SimpleDataComparison_1D;
typedef SimpleDataComparison<2, double, TH2D> SimpleDataComparison_2D;
typedef SimpleDataComparison<1, float, TH1F> SimpleDataComparison_1F;
typedef SimpleDataComparison<2, float, TH2F> SimpleDataComparison_2F;
typedef SimpleDataComparison<2, double, TH2Poly> SimpleDataComparison_2DPoly;
diff --git a/src/samples/nuA/BubbleChamber/ANL/ANL_CCQE_Evt_1DQ2_nu.cxx b/src/samples/nuA/BubbleChamber/ANL/ANL_CCQE_Evt_1DQ2_nu.cxx
index 5562945..53e7106 100644
--- a/src/samples/nuA/BubbleChamber/ANL/ANL_CCQE_Evt_1DQ2_nu.cxx
+++ b/src/samples/nuA/BubbleChamber/ANL/ANL_CCQE_Evt_1DQ2_nu.cxx
@@ -1,195 +1,190 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
//********************************************************************
#include "samples/SimpleDataComparison.hxx"
#include "utility/FullEventUtility.hxx"
#include "utility/KinematicUtility.hxx"
#include "utility/PDGCodeUtility.hxx"
using namespace nuis::event;
using namespace nuis::utility;
class ANL_CCQE_Evt_1DQ2_nu : public SimpleDataComparison_1D {
public:
NEW_NUIS_EXCEPT(invalid_publication_specifier);
enum Publication { kPRL31, kPRD16, kPRD26 };
Publication Pub;
std::string Pub_str;
bool UseD2Corr;
std::unique_ptr fD2CorrHist;
std::unique_ptr fPrediction_Uncorr;
ANL_CCQE_Evt_1DQ2_nu()
: Pub(kPRD26), Pub_str(""), UseD2Corr(false), fD2CorrHist(nullptr) {
ReadGlobalConfigDefaults();
}
std::string GetDocumentation() {
return "Can specify \"publication: \", where is one of [ PRL31, "
"PRD16, PRD26 ] to clarify a publication for comparison. Defaults "
"to PRD26.\n"
"Can enable deuterium Q2 correction by specifying "
"\"use_D2_correction: true\"";
}
fhicl::ParameterSet GetExampleConfiguration() {
fhicl::ParameterSet exps =
SimpleDataComparison_1D::GetExampleConfiguration();
exps.put("publication", "PRD26");
exps.put("use_D2_correction", false);
return exps;
}
void Initialize(fhicl::ParameterSet const &instance_sample_configuration) {
- if (instance_sample_configuration.has_key("verbosity")) {
- SetSampleVerbosity(
- instance_sample_configuration.get("verbosity"));
- }
-
std::string publication =
instance_sample_configuration.get("publication", "PRD26");
if (publication == "PRL31") {
Pub = kPRL31;
} else if (publication == "PRD16") {
Pub = kPRD16;
} else if (publication == "PRD26") {
Pub = kPRD26;
} else {
throw invalid_publication_specifier()
<< "[ERROR]: Found unexpected publication specifier "
<< std::quoted(publication)
<< ". Expected one of [ PRL31, PRD16, PRD26 ]";
}
switch (Pub) {
case kPRL31: {
Pub_str = "PRL31_844";
energy_cut = std::pair{0, 3E3};
IEventProcessor_INFO(
"Sample " << Name() << " specialized for publication: " << Pub_str);
break;
}
case kPRD16: {
Pub_str = "PRD16_3103";
energy_cut = std::pair{0, 6E3};
IEventProcessor_INFO(
"Sample " << Name() << " specialized for publication: " << Pub_str);
break;
}
case kPRD26: {
Pub_str = "PRD26_537";
energy_cut = std::pair{0, 6E3};
IEventProcessor_INFO(
"Sample " << Name() << " specialized for publication: " << Pub_str);
break;
}
}
fhicl::ParameterSet const &global_sample_configuration =
nuis::config::GetDocument().get(
std::string("global.sample_configuration.") + Name(),
fhicl::ParameterSet());
SetData(GetDataDir() + "nuA/BubbleChamber/ANL/CCQE/ANL_CCQE_Data_" +
Pub_str + ".root;ANL_1DQ2_Data");
SimpleDataComparison_1D::Initialize(instance_sample_configuration);
UseD2Corr = instance_sample_configuration.get(
"use_D2_correction",
global_sample_configuration.get("use_D2_correction", false));
if (UseD2Corr) {
fD2CorrHist = nuis::utility::GetHistogram(
GetDataDir() + "nuA/BubbleChamber/ANL/CCQE/"
"ANL_CCQE_Data_PRL31_844.root;ANL_1DQ2_Correction");
fPrediction_Uncorr = Clone(fPrediction, true);
}
// Signal selection function
IsSigFunc = [&](FullEvent const &fev) -> bool {
if (fev.mode != Channel_t::kCCQE) {
return false;
}
Particle ISNumu = GetHMISNeutralLepton(fev);
if (!ISNumu) {
return false;
}
if (ISNumu.pdg != pdgcodes::kNuMu) {
return false;
}
if (!energy_cut.IsInRange(ISNumu.P4.E())) {
return false;
}
double Q2 = GetNeutrinoQ2QERec(fev, 0);
if (Q2 <= 0) {
return false;
}
return true;
};
// 1D Projection function
CompProjFunc = [](FullEvent const &fev) -> std::array {
return {GetNeutrinoQ2QERec(fev, 0)};
};
}
// Used to apply D2 correction if requested
virtual void FillProjection(std::array const &proj,
double event_weight) {
if (UseD2Corr) {
TH_Help::Fill(fPrediction_Uncorr, proj, event_weight);
event_weight *= fD2CorrHist->Interpolate(proj[0]);
}
TH_Help::Fill(fPrediction, proj, event_weight);
}
void FinalizeComparison() {
SimpleDataComparison_1D::FinalizeComparison();
if (UseD2Corr) {
fPrediction_Uncorr->Scale(1.0, "width");
}
}
void Write() {
SimpleDataComparison_1D::Write();
if (UseD2Corr) {
nuis::persistency::WriteToOutputFile(
fPrediction_Uncorr, "Prediction_Uncorr", write_directory);
}
}
std::string Name() { return "ANL_CCQE_Evt_1DQ2_nu"; }
};
DECLARE_PLUGIN(IDataComparison, ANL_CCQE_Evt_1DQ2_nu);
DECLARE_PLUGIN(IEventProcessor, ANL_CCQE_Evt_1DQ2_nu);
diff --git a/src/samples/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt b/src/samples/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt
index d62dda3..39a9f8b 100644
--- a/src/samples/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt
+++ b/src/samples/nuA/Nuclear/T2K/CC0Pi/CMakeLists.txt
@@ -1,15 +1,18 @@
-SET(SAMPLES T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar)
+SET(SAMPLES
+ T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar
+ T2K_CC0PiNProt_CH_xsec_STV_nu)
foreach(S ${SAMPLES})
LIST(APPEND INuADataComparisons ${S})
LIST(APPEND SAMPLES_src ${S}.cxx)
endforeach()
add_library(T2KCC0PiDataComparisons SHARED ${SAMPLES_src})
-target_link_libraries(T2KCC0PiDataComparisons nuis_event nuis_input nuis_config nuis_persistency)
+target_link_libraries(T2KCC0PiDataComparisons
+ nuis_event nuis_input nuis_config nuis_persistency)
install(TARGETS T2KCC0PiDataComparisons DESTINATION plugins)
install(FILES T2K.CC0Pi.sample.config.fcl DESTINATION fcl)
SET(INuADataComparisons ${INuADataComparisons} PARENT_SCOPE)
diff --git a/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K.CC0Pi.sample.config.fcl b/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K.CC0Pi.sample.config.fcl
index fb2a079..6257eab 100644
--- a/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K.CC0Pi.sample.config.fcl
+++ b/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K.CC0Pi.sample.config.fcl
@@ -1,9 +1,20 @@
global.sample_configuration.T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar: {
journal_reference: "Not yet published"
target_material: "H20"
- flux_description: "ND280 FHC anti-muon neutrino"
+ flux_description: "ND280 RHC anti-muon neutrino"
signal_description: "CC0Pi"
shape_only: false
flux_unfolded: false
#has_data: false
}
+
+global.sample_configuration.T2K_CC0PiNProt_CH_xsec_STV_nu: {
+ journal_reference: "Phys. Rev. D 98, 032003"
+ doi: "https://doi.org/10.1103/PhysRevD.98.032003"
+ target_material: "CH"
+ flux_description: "ND280 FHC muon neutrino"
+ signal_description: "CC0Pi + N protons"
+ shape_only: false
+ flux_unfolded: false
+ #has_data: false
+}
diff --git a/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0PiNProt_CH_xsec_STV_nu.cxx b/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0PiNProt_CH_xsec_STV_nu.cxx
new file mode 100644
index 0000000..8f7ebd6
--- /dev/null
+++ b/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0PiNProt_CH_xsec_STV_nu.cxx
@@ -0,0 +1,190 @@
+// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
+
+/*******************************************************************************
+ * This file is part of NUISANCE.
+ *
+ * NUISANCE is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * NUISANCE is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with NUISANCE. If not, see .
+ *******************************************************************************/
+
+//********************************************************************
+
+#include "samples/SimpleDataComparison.hxx"
+
+#include "utility/EventTopologyUtility.hxx"
+#include "utility/FullEventUtility.hxx"
+#include "utility/KinematicUtility.hxx"
+#include "utility/PDGCodeUtility.hxx"
+
+using namespace nuis::event;
+using namespace nuis::utility;
+
+class T2K_CC0PiNProt_CH_xsec_STV_nu : public SimpleDataComparison_1D {
+
+ enum STVProjection { kDeltaPT, kDeltaPhiT, kDeltaAlphaT };
+
+public:
+ T2K_CC0PiNProt_CH_xsec_STV_nu() { ReadGlobalConfigDefaults(); }
+
+ std::string GetDocumentation() {
+ return "Can specify \"projection: \", where is one of [ "
+ "DeltaPT, DeltaPhiT, DeltaAlphaT ] to clarify a projection for "
+ "comparison. Defaults to DeltaPT.\n";
+ }
+
+ fhicl::ParameterSet GetExampleConfiguration() {
+ fhicl::ParameterSet exps =
+ SimpleDataComparison_1D::GetExampleConfiguration();
+
+ exps.put("projection", "DeltaPT");
+
+ return exps;
+ }
+
+ NEW_NUIS_EXCEPT(invalid_projection_specifier);
+
+ STVProjection GetProjection(std::string const &pstring) {
+ if ((pstring == "DeltaPT") || (pstring == "deltapt") ||
+ (pstring == "dpt")) {
+ return kDeltaPT;
+ } else if ((pstring == "DeltaPhiT") || (pstring == "deltaphit") ||
+ (pstring == "dphit")) {
+ return kDeltaPhiT;
+ } else if ((pstring == "DeltaAlphaT") || (pstring == "deltaalphat") ||
+ (pstring == "dat")) {
+ return kDeltaAlphaT;
+ }
+ throw invalid_projection_specifier()
+ << "Unknown projection " << std::quoted(pstring)
+ << " for comparison: " << std::quoted(Name());
+ }
+
+ void Initialize(fhicl::ParameterSet const &instance_sample_configuration) {
+
+ //! Set the verbosity of the sample logging macros.
+ if (instance_sample_configuration.has_key("verbosity")) {
+ SetSampleVerbosity(
+ instance_sample_configuration.get("verbosity"));
+ }
+
+ STVProjection proj =
+ GetProjection(instance_sample_configuration.get(
+ "projection", "DeltaPT"));
+
+ //! This will automatically set the data histogram to be loaded.
+ SetData(instance_sample_configuration.get("data_specifier"));
+
+ std::string projstr = "";
+ switch (proj) {
+ case kDeltaPT: {
+ projstr = "dpt";
+ break;
+ }
+ case kDeltaPhiT: {
+ projstr = "dphit";
+ break;
+ }
+ case kDeltaAlphaT: {
+ projstr = "dat";
+ break;
+ }
+ }
+
+ SetData(GetDataDir() + "nuA/Nuclear/T2K/CC0Pi/" + projstr +
+ "Results.root;Result");
+
+ // Perform any per-sample configuration in the base class
+ SimpleDataComparison_1D::Initialize(instance_sample_configuration);
+
+ //! Define your event signal here.
+ IsSigFunc = [](FullEvent const &fev) -> bool {
+ //! See src/utility/EventTopologyUtility.hxx for more pre-defined
+ //! topological signals.
+ if (!IsCC0Pi(fev)) {
+ return false;
+ }
+
+ //! See src/event/FullEvent.hxx for the full event class definition.
+ //! See src/utility/FullEventUtility.hxx for more helper methods for
+ //! interacting with the event class.
+
+ //! Get the initial state muon neutrino
+ Particle ISNumu = GetHMISParticle(fev, {pdgcodes::kNuMu});
+ //! An nuis::event::Particle that return true for !part is invalid and
+ //! does not exist on the particle stack. i.e. here, the selection fails
+ //! if the event didn't have an initial state numubar.
+ if (!ISNumu) {
+ return false;
+ }
+
+ //! Get the final state muon
+ Particle FSMu = GetHMFSParticle(fev, {pdgcodes::kMu});
+ if (!FSMu) {
+ return false;
+ }
+
+ //! Get the highest momentum final state proton
+ Particle FSProton = GetHMFSParticle(fev, {pdgcodes::kProton});
+ if (!FSProton) {
+ return false;
+ }
+
+ //! Cut on kinematic properties of the true final state.
+
+ // Muon phase space
+ // Pmu > 250 MeV, cos(theta_mu) > -0.6 (Sweet phase space!)
+ if ((FSMu.P() < 250) || (cos(ISNumu.P3().Angle(FSMu.P3())) < -0.6)) {
+ return false;
+ }
+
+ // Proton phase space
+ // Pprot > 450 MeV, cos(theta_proton) > 0.4
+ if ((FSProton.P() < 450) || (FSProton.P() > 1E3) ||
+ (cos(ISNumu.P3().Angle(FSProton.P3())) < 0.4)) {
+ return false;
+ }
+
+ //! Select the event!
+ return true;
+ };
+
+ //! 1D Projection function
+ //! This function takes selected events and returns an array of size, the
+ //! dimensionality of the comparisons (SimpleDataComparison_1D::NDim)
+ CompProjFunc = [=](FullEvent const &fev)
+ -> std::array {
+ switch (proj) {
+ case kDeltaPT: {
+ return {GetDeltaPT_CC0PiN(fev).Mag() * 1E-3};
+ }
+ case kDeltaPhiT: {
+ return {GetDeltaPhiT_CC0PiN(fev)};
+ }
+ case kDeltaAlphaT: {
+ return {GetDeltaAlphaT_CC0PiN(fev)};
+ }
+ }
+ return {-std::numeric_limits::max()};
+ };
+ }
+
+ std::string Name() { return "T2K_CC0PiNProt_CH_xsec_STV_nu"; }
+
+ //! Here you can write any custom histograms to TTrees that your sample has
+ //! been handling.
+ void Write() { SimpleDataComparison_1D::Write(); }
+};
+
+//! These declarations allow your class to be loaded dynamically by NUISANCE
+DECLARE_PLUGIN(IDataComparison, T2K_CC0PiNProt_CH_xsec_STV_nu);
+DECLARE_PLUGIN(IEventProcessor, T2K_CC0PiNProt_CH_xsec_STV_nu);
diff --git a/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar.cxx b/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar.cxx
index 7576869..a225472 100644
--- a/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar.cxx
+++ b/src/samples/nuA/Nuclear/T2K/CC0Pi/T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar.cxx
@@ -1,194 +1,182 @@
// Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret
/*******************************************************************************
* This file is part of NUISANCE.
*
* NUISANCE is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* NUISANCE is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NUISANCE. If not, see .
*******************************************************************************/
//********************************************************************
#include "samples/SimpleDataComparison.hxx"
#include "utility/EventTopologyUtility.hxx"
#include "utility/FullEventUtility.hxx"
#include "utility/KinematicUtility.hxx"
#include "utility/PDGCodeUtility.hxx"
using namespace nuis::event;
using namespace nuis::utility;
-using SimpleDataComparison_2DPoly = SimpleDataComparison<2, double, TH2Poly>;
-
class T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar
: public SimpleDataComparison_2DPoly {
std::unique_ptr fPrediction_pmu;
std::unique_ptr fPrediction_ctheta;
std::unique_ptr fPrediction_fine;
std::unique_ptr fPrediction_unscale;
public:
T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar() { ReadGlobalConfigDefaults(); }
std::string GetDocumentation() { return ""; }
fhicl::ParameterSet GetExampleConfiguration() {
fhicl::ParameterSet exps =
SimpleDataComparison_2DPoly::GetExampleConfiguration();
return exps;
}
void Initialize(fhicl::ParameterSet const &instance_sample_configuration) {
- if (instance_sample_configuration.has_key("verbosity")) {
- SetSampleVerbosity(
- instance_sample_configuration.get("verbosity"));
- }
-
- fhicl::ParameterSet const &global_sample_configuration =
- nuis::config::GetDocument().get(
- std::string("global.sample_configuration.") + Name(),
- fhicl::ParameterSet());
-
SetData(
GetDataDir() +
"nuA/Nuclear/T2K/CC0Pi/H2O_xsec_2Dpmuthetamu_numubar.root;datapoly");
fPrediction_fine = std::unique_ptr(
new TH2D("fPrediction_fine",
";#it{p}_{#mu} "
"(MeV/#it{c});cos(#theta_{#mu});d^{2}#sigma/"
"d#it{p}_{#mu}dcos(#theta_{#mu}) (cm^{2} MeV^{-1} A^{-1})",
100, 0, 3000, 50, 0.84, 1));
fPrediction_pmu =
std::unique_ptr(new TH1D("fPrediction_pmu",
";#it{p}_{#mu} "
"(MeV/#it{c});d#sigma/"
"d#it{p}_{#mu} (cm^{2} MeV^{-1} A^{-1})",
100, 0, 3000));
fPrediction_ctheta =
std::unique_ptr(new TH1D("fPrediction_ctheta",
";cos(#theta_{#mu});d#sigma/"
"ddcos(#theta_{#mu}) (cm^{2} A^{-1})",
50, 0.84, 1));
SimpleDataComparison_2DPoly::Initialize(instance_sample_configuration);
fPrediction_unscale =
nuis::utility::Clone(fPrediction, false, "Prediction_unscale");
// Signal selection function
IsSigFunc = [](FullEvent const &fev) -> bool {
if (!IsCC0Pi(fev)) {
return false;
}
Particle ISNumuBar = GetHMISParticle(fev, {pdgcodes::kNuMuBar});
if (!ISNumuBar) {
return false;
}
Particle FSMuPlus = GetHMFSParticle(fev, {pdgcodes::kMuPlus});
if (!FSMuPlus) {
return false;
}
if (FSMuPlus.CosTheta() < 0.84) {
return false;
}
return true;
};
// 1D Projection function
CompProjFunc = [](FullEvent const &fev) -> std::array {
Particle FSMuPlus = GetHMFSParticle(fev, {pdgcodes::kMuPlus});
return {FSMuPlus.P(), FSMuPlus.CosTheta()};
};
ProcessExtraFunc = [&](FullEvent const &fev, bool isSel, double weight) {
if (isSel) {
Particle FSMuPlus = GetHMFSParticle(fev, {pdgcodes::kMuPlus});
TH_Helper::Fill(fPrediction_fine,
{FSMuPlus.P(), FSMuPlus.CosTheta()}, weight);
TH_Help::Fill(fPrediction_unscale, {FSMuPlus.P(), FSMuPlus.CosTheta()});
}
Particle ISNumuBar = GetHMISParticle(fev, {pdgcodes::kNuMuBar});
if (!ISNumuBar) {
return;
}
Particle FSMuPlus = GetHMFSParticle(fev, {pdgcodes::kMuPlus});
if (!FSMuPlus) {
return;
}
TH_Helper::Fill(fPrediction_pmu, {FSMuPlus.P()}, weight);
TH_Helper::Fill(fPrediction_ctheta, {FSMuPlus.CosTheta()}, weight);
};
}
std::string Name() { return "T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar"; }
void Write() {
SimpleDataComparison_2DPoly::Write();
nuis::persistency::WriteToOutputFile(
fPrediction_fine, "Prediction_fine", write_directory);
nuis::persistency::WriteToOutputFile(
fPrediction_pmu, "fPrediction_pmu", write_directory);
nuis::persistency::WriteToOutputFile(
fPrediction_ctheta, "fPrediction_ctheta", write_directory);
nuis::persistency::WriteToOutputFile(
fPrediction_unscale, "Prediction_unscale", write_directory);
static std::vector> const
binning = {
{// Slice_0
YPolyBinSpec(450, 0.85), YPolyBinSpec(450, 0.95)},
{// Slice_1
YPolyBinSpec(550, 0.86), YPolyBinSpec(550, 0.93),
YPolyBinSpec(550, 0.97)},
{// Slice_2
YPolyBinSpec(700, 0.89), YPolyBinSpec(700, 0.94),
YPolyBinSpec(700, 0.98)},
{// Slice_3
YPolyBinSpec(900, 0.91), YPolyBinSpec(900, 0.95),
YPolyBinSpec(900, 0.98)},
{// Slice_4
YPolyBinSpec(1200, 0.92), YPolyBinSpec(1200, 0.96),
YPolyBinSpec(1200, 0.98)},
{// Slice_4
YPolyBinSpec(1600, 0.93), YPolyBinSpec(1600, 0.97),
YPolyBinSpec(1600, 0.99)},
{// Slice_5
YPolyBinSpec(2500, 0.96), YPolyBinSpec(2500, 0.99)},
};
for (auto &slice : GetTH2PolySlices(fData, binning)) {
nuis::persistency::WriteToOutputFile(slice, slice->GetName(),
write_directory);
}
for (auto &slice : GetTH2PolySlices(fPrediction_xsec, binning)) {
nuis::persistency::WriteToOutputFile(slice, slice->GetName(),
write_directory);
}
}
};
DECLARE_PLUGIN(IDataComparison, T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar);
DECLARE_PLUGIN(IEventProcessor, T2K_CC0Pi_H2O_xsec_2Dpcthetamu_numubar);
diff --git a/src/utility/CMakeLists.txt b/src/utility/CMakeLists.txt
index 348e8e9..5f8c5d6 100644
--- a/src/utility/CMakeLists.txt
+++ b/src/utility/CMakeLists.txt
@@ -1,27 +1,28 @@
set(utility_implementation_files
FileSystemUtility.cxx
FullEventUtility.cxx
KinematicUtility.cxx
PDGCodeUtility.cxx
StringUtility.cxx
- EventTopologyUtility.cxx)
+ EventTopologyUtility.cxx
+ HistogramUtility.cxx)
set(utility_header_files
FileSystemUtility.hxx
FullEventUtility.hxx
InteractionChannelUtility.hxx
KinematicUtility.hxx
PDGCodeUtility.hxx
ROOTUtility.hxx
HistogramUtility.hxx
StringUtility.hxx
TerminalUtility.hxx
StatsUtility.hxx
EventTopologyUtility.hxx)
add_library(nuis_utility SHARED ${utility_implementation_files})
target_link_libraries(nuis_utility)
install(TARGETS nuis_utility DESTINATION lib)
install(FILES ${utility_header_files} DESTINATION include/utility)
diff --git a/src/utility/HistogramUtility.cxx b/src/utility/HistogramUtility.cxx
new file mode 100644
index 0000000..9011cef
--- /dev/null
+++ b/src/utility/HistogramUtility.cxx
@@ -0,0 +1,102 @@
+#include "utility/HistogramUtility.hxx"
+
+namespace nuis {
+namespace utility {
+bool IsFlowBin(TAxis const &ax, Int_t bin_it) {
+ return ((bin_it <= 0) || (bin_it >= (ax.GetNbins() + 1)));
+}
+
+bool IsInHistogramRange(TAxis const &ax, double v) {
+ Int_t bin_it = ax.FindFixBin(v);
+ return !IsFlowBin(ax, bin_it);
+}
+
+NEW_NUIS_EXCEPT(invalid_TH2Poly);
+NEW_NUIS_EXCEPT(invalid_PolyBinSpecifierList);
+
+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
diff --git a/src/utility/HistogramUtility.hxx b/src/utility/HistogramUtility.hxx
index 851d84c..40a8273 100644
--- a/src/utility/HistogramUtility.hxx
+++ b/src/utility/HistogramUtility.hxx
@@ -1,522 +1,430 @@
// 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
+#pragma once
#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)));
-}
+bool IsFlowBin(TAxis const &ax, Int_t bin_it);
-inline bool IsInHistogramRange(TAxis const &ax, double v) {
- Int_t bin_it = ax.FindFixBin(v);
- return !IsFlowBin(ax, bin_it);
-}
+bool IsInHistogramRange(TAxis const &ax, double v);
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) {
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) {
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) {
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) {
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) {
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) {
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