diff --git a/src/app/nuisstudy.cxx b/src/app/nuisstudy.cxx index 17e18b8..bd54382 100644 --- a/src/app/nuisstudy.cxx +++ b/src/app/nuisstudy.cxx @@ -1,56 +1,59 @@ #include "config/GlobalConfiguration.hxx" -#include "input/IInputHandler.hxx" +#include "input/InputManager.hxx" #include "event/MinimalEvent.hxx" #include "samples/IEventProcessor.hxx" #include "plugins/Instantiate.hxx" #include "exception/exception.hxx" #include "persistency/ROOTOutput.hxx" #include "fhiclcpp/make_ParameterSet.h" #include NEW_NUIS_EXCEPT(invalid_cli_arguments); int main(int argc, char const *argv[]) { nuis::config::EnsureConfigurationRead("nuis.global.config.fcl"); nuis::config::EnsureConfigurationRead("nuis.datacomparisons.fcl"); if (argc != 2 || (std::string(argv[1]) == "-?") || (std::string(argv[1]) == "--help")) { throw invalid_cli_arguments() << "[ERROR]: Expected to be passed a single FHiCL file name or " "absolute or relative path. N.B. Files in the local directory must " "be fully qualified like \"$ " << argv[0] << " ./myconf.fcl\"."; } nuis::config::EnsureConfigurationRead(argv[1]); size_t NMax = nuis::config::GetDocument().get( "nmax", std::numeric_limits::max()); for (fhicl::ParameterSet const &samp_config : nuis::config::GetDocument().get>( "samples")) { std::cout << "[INFO]: Reading sample: " << samp_config.get("name") << std::endl; nuis::plugins::plugin_traits::unique_ptr_t sample = nuis::plugins::Instantiate( samp_config.get("name")); sample->Initialize(samp_config); sample->ProcessSample(NMax); sample->Write(); + + //Ensures no re-use of samples but cleans up the memory. + nuis::input::InputManager::Get().Clear(); } nuis::persistency::CloseOpenTFiles(); } diff --git a/src/generator/input/NEUTInputHandler.cxx b/src/generator/input/NEUTInputHandler.cxx index 1d9ea8e..ea62b19 100644 --- a/src/generator/input/NEUTInputHandler.cxx +++ b/src/generator/input/NEUTInputHandler.cxx @@ -1,254 +1,257 @@ #include "generator/input/NEUTInputHandler.hxx" #include "persistency/ROOTOutput.hxx" #include "utility/HistogramUtility.hxx" #include "utility/InteractionChannelUtility.hxx" #include "utility/PDGCodeUtility.hxx" #include "generator/utility/NEUTUtility.hxx" #include "neutpart.h" #include "neutvect.h" #include "fhiclcpp/ParameterSet.h" using namespace nuis; using namespace nuis::event; using namespace nuis::utility; using namespace nuis::neuttools; NEUTInputHandler::NEUTInputHandler() : fInputTreeFile() {} NEUTInputHandler::NEUTInputHandler(NEUTInputHandler &&other) : fInputTreeFile(std::move(other.fInputTreeFile)), fReaderEvent(std::move(other.fReaderEvent)) {} void NEUTInputHandler::Initialize(fhicl::ParameterSet const &ps) { fInputTreeFile = CheckGetTTree(ps.get("file"), "neuttree"); fNeutVect = nullptr; fInputTreeFile.tree->SetBranchAddress("vectorbranch", &fNeutVect); fKeepIntermediates = ps.get("keep_intermediates", false); std::string flux_name = ps.get("flux_hist_name", "flux_numu"); std::string evtrt_name = ps.get("evtrt_hist_name", "evtrt_numu"); + fXSecRescaleFactor = ps.get("xsec_rescale_factor",1); std::string override_flux_file = ps.get("override_flux_file", ""); if (override_flux_file.size()) { fFlux = GetHistogramFromROOTFile(override_flux_file, flux_name); RebuildEventRate(!ps.get("flux_hist_in_MeV", false)); } else { fFlux = GetHistogramFromROOTFile(fInputTreeFile.file, flux_name, false); } if (fFlux) { if (!fEvtrt) { fEvtrt = GetHistogramFromROOTFile(fInputTreeFile.file, evtrt_name); } fFileWeight = fEvtrt->Integral() * 1.0E-38 / (fFlux->Integral() * double(GetNEvents())); std::cout << "[INFO]: Average NEUT XSecWeight = " << fFileWeight << "" << std::endl; } else { std::cout << "[INFO]: Input NEUT file doesn't have flux information, assuming " "mono energetic and calculating average cross-section..." << std::endl; fFileWeight = GetMonoEXSecWeight() / double(GetNEvents()); std::cout << "[INFO]: Done (Average NEUT XSecWeight = " << fFileWeight << ")!" << std::endl; } + + fFileWeight *= fXSecRescaleFactor; } NEW_NUIS_EXCEPT(non_mono_energetic_input_file); double NEUTInputHandler::GetMonoEXSecWeight() { double xsec = 0; double count = 0; double E_first = std::numeric_limits::max(); size_t NEvents = GetNEvents(); size_t ShoutEvery = NEvents / 100; std::cout << "[INFO]: Read " << 0 << "/" << NEvents << " NEUT events." << std::flush; for (size_t nev_it = 0; nev_it < NEvents; ++nev_it) { if (ShoutEvery && !(nev_it % ShoutEvery)) { std::cout << "\r[INFO]: Read " << nev_it << "/" << NEvents << " NEUT events." << std::flush; } fInputTreeFile.tree->GetEntry(nev_it); xsec += fNeutVect->Totcrs; count++; if (E_first == std::numeric_limits::max()) { NeutPart *part = fNeutVect->PartInfo(0); E_first = part->fP.E(); } else { NeutPart *part = fNeutVect->PartInfo(0); if (E_first != part->fP.E()) { throw non_mono_energetic_input_file() << "[ERROR]: When calculating NEUT xsec weight for a mono " "energetic beam, found different energies: " << E_first << " and " << part->fP.E() << ". Cannot continue with this input file."; } } } std::cout << "\r[INFO]: Read " << NEvents << "/" << NEvents << " NEUT events." << std::endl; return (xsec / count) * 1E-38; } void NEUTInputHandler::RebuildEventRate(bool FluxInGeV) { auto XSec = Clone(fFlux, true, "xsec"); auto Entry = Clone(fFlux, true, "entry"); std::cout << "[INFO]: Rebuilding total cross-section prediction..." << std::endl; size_t NEvents = GetNEvents(); size_t ShoutEvery = NEvents / 100; std::cout << "[INFO]: Read " << 0 << "/" << NEvents << " NEUT events." << std::flush; for (size_t nev_it = 0; nev_it < NEvents; ++nev_it) { if (ShoutEvery && !(nev_it % ShoutEvery)) { std::cout << "\r[INFO]: Read " << nev_it << "/" << NEvents << " NEUT events." << std::flush; } fInputTreeFile.tree->GetEntry(nev_it); NeutPart *part = fNeutVect->PartInfo(0); double E = part->fP.E() * (FluxInGeV ? 1E-3 : 1); XSec->Fill(E, fNeutVect->Totcrs); Entry->Fill(E); } std::cout << std::endl << "[INFO]: Done!" << std::endl; // Binned total xsec XSec->Divide(Entry.get()); fEvtrt = Clone(XSec, false, "evtrt"); // Event rate prediction fEvtrt->Multiply(fFlux.get()); } double NEUTInputHandler::GetXSecScaleFactor(std::pair const &) const { throw input_handler_feature_unimplemented() << "[ERROR]: Flux cuts not yet implemented for NEUT input handler."; } NeutVect *NEUTInputHandler::GetNeutEvent(ev_index_t idx) const { if (idx >= GetNEvents()) { throw IInputHandler::invalid_entry() << "[ERROR]: Attempted to get entry " << idx << " from an InputHandler with only " << GetNEvents(); } fInputTreeFile.tree->GetEntry(idx); return fNeutVect; } MinimalEvent const &NEUTInputHandler::GetMinimalEvent(ev_index_t idx) const { if (idx >= GetNEvents()) { throw IInputHandler::invalid_entry() << "[ERROR]: Attempted to get entry " << idx << " from an InputHandler with only " << GetNEvents(); } fInputTreeFile.tree->GetEntry(idx); fReaderEvent.mode = IntToChannel(fNeutVect->Mode); size_t NPart = fNeutVect->Npart(); for (size_t part_it = 0; part_it < NPart; part_it++) { NeutPart const &part = (*fNeutVect->PartInfo(part_it)); if ((part.fIsAlive == false) && (part.fStatus == -1) && IsNeutralLepton(part.fPID)) { fReaderEvent.probe_E = part.fP.T(); fReaderEvent.probe_pdg = part.fPID; break; } } fReaderEvent.XSecWeight = fFileWeight; if (fWeightCache.size() <= idx) { fWeightCache.push_back(fReaderEvent.XSecWeight); } fReaderEvent.fGenEvent = fNeutVect; return fReaderEvent; } FullEvent const &NEUTInputHandler::GetFullEvent(ev_index_t idx) const { (void)GetMinimalEvent(idx); fReaderEvent.ClearParticleStack(); if (fNeutVect->Ibound) { fReaderEvent.target_pdg = MakeNuclearPDG(fNeutVect->TargetA, fNeutVect->TargetZ); } else { fReaderEvent.target_pdg = MakeNuclearPDG(1, 1); } size_t NPart = fNeutVect->Npart(); size_t NPrimary = fNeutVect->Nprimary(); for (size_t part_it = 0; part_it < NPart; part_it++) { NeutPart const &part = (*fNeutVect->PartInfo(part_it)); Particle nuis_part; nuis_part.pdg = part.fPID; nuis_part.P4 = part.fP; Particle::Status_t state = GetNeutParticleStatus(part, fReaderEvent.mode); if (!fKeepIntermediates && (state == Particle::Status_t::kIntermediate)) { continue; } size_t state_int = static_cast(state); // Add status == 0 particles as pre-FSI particles until we find an // intermediate state particle if ((part_it < NPrimary) && (state != Particle::Status_t::kPrimaryInitialState)) { fReaderEvent .ParticleStack[static_cast( Particle::Status_t::kPrimaryFinalState)] .particles.push_back(nuis_part); } fReaderEvent.ParticleStack[state_int].particles.push_back(nuis_part); } return fReaderEvent; } double NEUTInputHandler::GetEventWeight(ev_index_t idx) const { if (idx > fWeightCache.size()) { throw weight_cache_miss() << "[ERROR]: Failed to get cached weight for event index: " << idx; } return fWeightCache[idx]; } size_t NEUTInputHandler::GetNEvents() const { return fInputTreeFile.tree->GetEntries(); } GeneratorManager::Generator_id_t NEUTInputHandler::GetGeneratorId() const { return GeneratorManager::Get().EnsureGeneratorRegistered("NEUT"); } DECLARE_PLUGIN(IInputHandler, NEUTInputHandler); diff --git a/src/generator/input/NEUTInputHandler.hxx b/src/generator/input/NEUTInputHandler.hxx index 5418829..a2cc915 100644 --- a/src/generator/input/NEUTInputHandler.hxx +++ b/src/generator/input/NEUTInputHandler.hxx @@ -1,70 +1,71 @@ // 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; + double fXSecRescaleFactor; std::unique_ptr fFlux; std::unique_ptr fEvtrt; double fFileWeight; bool fKeepIntermediates; double GetMonoEXSecWeight(); void RebuildEventRate(bool FluxInGeV=true); public: NEW_NUIS_EXCEPT(weight_cache_miss); NEUTInputHandler(); NEUTInputHandler(NEUTInputHandler const &) = delete; NEUTInputHandler(NEUTInputHandler &&); void Initialize(fhicl::ParameterSet const &); NeutVect *GetNeutEvent(ev_index_t) const; nuis::event::MinimalEvent const &GetMinimalEvent(ev_index_t idx) const; nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const; double GetEventWeight(ev_index_t idx) const; size_t GetNEvents() const; double GetXSecScaleFactor(std::pair const &EnuRange) const; nuis::GeneratorManager::Generator_id_t GetGeneratorId() const; }; diff --git a/src/input/InputManager.hxx b/src/input/InputManager.hxx index 59b75f5..29fa090 100644 --- a/src/input/InputManager.hxx +++ b/src/input/InputManager.hxx @@ -1,62 +1,64 @@ // 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 "input/IInputHandler.hxx" #include "plugins/traits.hxx" #include "exception/exception.hxx" #include #include namespace fhicl { class ParameterSet; } namespace nuis { namespace input { class InputManager { struct NamedInputHandler { NamedInputHandler(std::string const &, plugins::plugin_traits::unique_ptr_t &&); std::string name; plugins::plugin_traits::unique_ptr_t handler; }; std::vector Inputs; InputManager(); static InputManager *_global_inst; public: static InputManager &Get(); NEW_NUIS_EXCEPT(unknown_input); NEW_NUIS_EXCEPT(unsupported_input_type); using Input_id_t = size_t; Input_id_t EnsureInputLoaded(fhicl::ParameterSet const &); Input_id_t GetInputId(std::string const &) const; IInputHandler const &GetInputHandler(Input_id_t) const; + + void Clear() { Inputs.clear(); } }; } // namespace input } // namespace nuis diff --git a/src/plugins/traits.hxx b/src/plugins/traits.hxx index b797f67..9a884a8 100644 --- a/src/plugins/traits.hxx +++ b/src/plugins/traits.hxx @@ -1,73 +1,67 @@ // Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ #pragma once #include #include #include namespace nuis { namespace plugins { template struct plugin_traits {}; } // namespace plugins } // namespace nuis #define DECLARE_PLUGIN_INTERFACE(INTERFACE_CLASS_NAME) \ namespace nuis { \ namespace plugins { \ template <> struct plugin_traits { \ typedef std::unique_ptr> \ unique_ptr_t; \ static std::string interface_name() { return #INTERFACE_CLASS_NAME; } \ static std::string \ instantiator_function_name(std::string const &classname) { \ return std::string("nuis_plugins_") + \ plugin_traits::interface_name() + "_" + \ classname + "_instantiator"; \ } \ - static std::string \ - check_manifest_function_name(std::string const &classname) { \ - return std::string("nuis_plugins_") + \ - plugin_traits::interface_name() + "_" + \ - classname + "_check_manifest"; \ - } \ static std::string deleter_function_name(std::string const &classname) { \ return std::string("nuis_plugins_") + \ plugin_traits::interface_name() + "_" + \ classname + "_deleter"; \ } \ }; \ } \ } #define BUILDMETHODNAME(CN, PN, POST) nuis_plugins_##CN##_##PN##POST #define DECLARE_PLUGIN(INTERFACE_CLASS_NAME, PLUGIN_CLASS_NAME) \ extern "C" { \ void *BUILDMETHODNAME(INTERFACE_CLASS_NAME, PLUGIN_CLASS_NAME, \ _instantiator)() { \ return new PLUGIN_CLASS_NAME(); \ } \ void BUILDMETHODNAME(INTERFACE_CLASS_NAME, PLUGIN_CLASS_NAME, \ _deleter)(void *instance) { \ delete reinterpret_cast(instance); \ } \ } diff --git a/src/utility/HistogramUtility.cxx b/src/utility/HistogramUtility.cxx index 4246fe1..0528485 100644 --- a/src/utility/HistogramUtility.cxx +++ b/src/utility/HistogramUtility.cxx @@ -1,152 +1,154 @@ #include "utility/HistogramUtility.hxx" #include #include 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; } void SliceNorm(std::unique_ptr &hist, bool AlongY, char const *opt) { std::string opt_str = opt; std::transform(opt_str.begin(), opt_str.end(), opt_str.begin(), ::tolower); bool width = (opt_str.find("width") != std::string::npos); for (int slice_it = 0; slice_it < (AlongY ? hist->GetXaxis() : hist->GetYaxis())->GetNbins(); ++slice_it) { int xl = AlongY ? slice_it + 1 : 1; int xu = AlongY ? slice_it + 1 : hist->GetXaxis()->GetNbins(); int yl = AlongY ? 1 : slice_it + 1; int yu = AlongY ? hist->GetYaxis()->GetNbins() : slice_it + 1; double integral = hist->Integral(xl, xu, yl, yu, opt); - if (integral < std::numeric_limits::epsilon()) { + if (integral <= 0) { continue; } + double integ_check= 0; for (int bin_it = 0; bin_it < (AlongY ? hist->GetYaxis() : hist->GetXaxis())->GetNbins(); ++bin_it) { double s = 1.0 / integral; if (width) { if (AlongY) { s /= (hist->GetXaxis()->GetBinWidth(slice_it + 1) * hist->GetYaxis()->GetBinWidth(bin_it + 1)); } else { s /= (hist->GetYaxis()->GetBinWidth(slice_it + 1) * hist->GetXaxis()->GetBinWidth(bin_it + 1)); } } if (AlongY) { hist->SetBinContent(slice_it + 1, bin_it + 1, hist->GetBinContent(slice_it + 1, bin_it + 1) * s); - hist->SetBinContent(slice_it + 1, bin_it + 1, - hist->GetBinContent(slice_it + 1, bin_it + 1) * s); + hist->SetBinError(slice_it + 1, bin_it + 1, + hist->GetBinError(slice_it + 1, bin_it + 1) * s); } else { hist->SetBinContent(bin_it + 1, slice_it + 1, hist->GetBinContent(bin_it + 1, slice_it + 1) * s); - hist->SetBinContent(bin_it + 1, slice_it + 1, - hist->GetBinContent(bin_it + 1, slice_it + 1) * s); + hist->SetBinError(bin_it + 1, slice_it + 1, + hist->GetBinError(bin_it + 1, slice_it + 1) * s); } } + } } } // namespace utility } // namespace nuis