diff --git a/src/generator/input/CMakeLists.txt b/src/generator/input/CMakeLists.txt index fc9daea..2c19d81 100644 --- a/src/generator/input/CMakeLists.txt +++ b/src/generator/input/CMakeLists.txt @@ -1,46 +1,51 @@ if(USE_NUWRO) LIST(APPEND INPUT_HANDLERS_IMPL NuWroInputHandler.cxx) + LIST(APPEND INPUT_HANDLERS_HDRS NuWroInputHandler.hxx) endif(USE_NUWRO) if(USE_NEUT) LIST(APPEND INPUT_HANDLERS_IMPL NEUTInputHandler.cxx) + LIST(APPEND INPUT_HANDLERS_HDRS NEUTInputHandler.hxx) endif(USE_NEUT) if(USE_GENIE) LIST(APPEND INPUT_HANDLERS_IMPL GENIEInputHandler.cxx) + LIST(APPEND INPUT_HANDLERS_HDRS GENIEInputHandler.hxx) endif(USE_GENIE) if(INPUT_HANDLERS_IMPL) add_library(InputHandlers SHARED ${INPUT_HANDLERS_IMPL}) target_link_libraries(InputHandlers ${INPUT_HANDLERS_LINK_LIBS}) if(USE_GENIE) target_compile_options(InputHandlers PUBLIC ${GENIE_CXX_FLAGS}) include_directories(${GENIE_INCLUDE_DIRS}) target_link_libraries(InputHandlers ${GENIE_LINK_DIRS}) target_link_libraries(InputHandlers -Wl,--no-as-needed ${GENIE_LIBS}) endif() if(USE_NEUT) target_compile_options(InputHandlers PUBLIC ${NEUT_CXX_FLAGS}) include_directories(${NEUT_INCLUDE_DIRS}) target_link_libraries(InputHandlers ${NEUT_LINK_DIRS}) target_link_libraries(InputHandlers -Wl,--no-as-needed ${NEUT_LIBS}) target_link_libraries(InputHandlers ${NEUT_IMPORTED_TARGETS}) endif() if(USE_NUWRO) target_compile_options(InputHandlers PUBLIC ${NUWRO_CXX_FLAGS}) include_directories(${NUWRO_INCLUDE_DIRS}) target_link_libraries(InputHandlers -Wl,--no-as-needed ${NUWRO_IMPORTED_TARGETS}) endif() if(${NEED_PYTHIA6}) target_link_libraries(InputHandlers -L${PYTHIA6}) endif() target_link_libraries(InputHandlers -L${ROOT_LIBDIR} ${ROOT_LIBS}) install(TARGETS InputHandlers DESTINATION plugins) + + install(FILES ${INPUT_HANDLERS_HDRS} DESTINATION include/generator/input) endif() diff --git a/src/generator/input/GENIEInputHandler.cxx b/src/generator/input/GENIEInputHandler.cxx index 68a169c..baa784a 100644 --- a/src/generator/input/GENIEInputHandler.cxx +++ b/src/generator/input/GENIEInputHandler.cxx @@ -1,149 +1,149 @@ #include "generator/input/GENIEInputHandler.hxx" #include "generator/utility/GENIEUtility.hxx" #include "utility/PDGCodeUtility.hxx" #ifdef GENIE_V3_INTERFACE #include "Framework/EventGen/EventRecord.h" #include "Framework/GHEP/GHepParticle.h" #include "Framework/GHEP/GHepRecord.h" #else #include "EVGCore/EventRecord.h" #include "GHEP/GHepParticle.h" #include "GHEP/GHepRecord.h" #endif #include "fhiclcpp/ParameterSet.h" using namespace nuis; using namespace nuis::event; using namespace nuis::utility; using namespace nuis::genietools; GENIEInputHandler::GENIEInputHandler() : fInputTree(), fGenieNtpl(nullptr) {} GENIEInputHandler::GENIEInputHandler(GENIEInputHandler &&other) : fInputTree(std::move(other.fInputTree)), fReaderEvent(std::move(other.fReaderEvent)), fGenieNtpl(nullptr) {} void GENIEInputHandler::Initialize(fhicl::ParameterSet const &ps) { fInputTree = CheckGetTTree(ps.get("file"), "gtree"); fInputTree.tree->SetBranchAddress("gmcrec", &fGenieNtpl); fKeepIntermediates = ps.get("keep_intermediates", false); fKeepNuclearParticles = ps.get("keep_nuclear_particles", false); } MinimalEvent const &GENIEInputHandler::GetMinimalEvent(ev_index_t idx) const { if (idx >= GetNEvents()) { throw IInputHandler::invalid_entry() << "[ERROR]: Attempted to get entry " << idx << " from an InputHandler with only " << GetNEvents(); } fInputTree.tree->GetEntry(idx); genie::GHepRecord *GHep = static_cast(fGenieNtpl->event); if (!GHep) { throw invalid_GENIE_event() << "[ERROR]: GENIE event " << idx << " failed to contain a GHepRecord"; } fReaderEvent.mode = GetEventChannel(*GHep); TObjArrayIter iter(GHep); genie::GHepParticle *p; while ((p = (dynamic_cast((iter).Next())))) { if (!p) { continue; } Particle::Status_t state = GetParticleStatus(*p, fReaderEvent.mode); if (state != Particle::Status_t::kPrimaryInitialState) { continue; } if (!IsNeutralLepton(p->Pdg(), pdgcodes::kMatterAntimatter) && !IsChargedLepton(p->Pdg(), pdgcodes::kMatterAntimatter)) { continue; } fReaderEvent.probe_E = p->E() * 1.E3; fReaderEvent.probe_pdg = p->Pdg(); break; } fReaderEvent.XSecWeight = 1; if (fWeightCache.size() <= idx) { fWeightCache.push_back(fReaderEvent.XSecWeight); } return fReaderEvent; } FullEvent const &GENIEInputHandler::GetFullEvent(ev_index_t idx) const { (void)GetMinimalEvent(idx); fReaderEvent.ClearParticleStack(); genie::GHepRecord *GHep = static_cast(fGenieNtpl->event); unsigned int npart = GHep->GetEntries(); // Fill Particle Stack genie::GHepParticle *p = 0; TObjArrayIter iter(GHep); // Loop over all particles while ((p = (dynamic_cast((iter).Next())))) { if (!p) { continue; } // Get Status Particle::Status_t state = GetParticleStatus(*p, fReaderEvent.mode); if (!fKeepIntermediates) { continue; } if (!fKeepNuclearParticles && IsNuclearPDG(p->Pdg())) { continue; } Particle nuis_part; nuis_part.pdg = p->Pdg(); nuis_part.P4 = TLorentzVector(p->Px(), p->Py(), p->Pz(), p->E()) * 1E3; fReaderEvent.ParticleStack[static_cast(state)].particles.push_back( nuis_part); } return fReaderEvent; } double GENIEInputHandler::GetEventWeight(ev_index_t idx) const { if (idx > fWeightCache.size()) { throw weight_cache_miss() << "[ERROR]: Failed to get cached weight for event index: " << idx; } return fWeightCache[idx]; } double GENIEInputHandler::GetXSecScaleFactor( std::pair const &EnuRange) const { throw input_handler_feature_unimplemented() << "[ERROR]: Flux cuts not yet implemented for GENIE input handler."; } size_t GENIEInputHandler::GetNEvents() const { return fInputTree.tree->GetEntries(); } -GeneratorManager::Generator_id_t GENIEInputHandler::GetGeneratorId() { +GeneratorManager::Generator_id_t GENIEInputHandler::GetGeneratorId() const { return GeneratorManager::Get().EnsureGeneratorRegistered("GENIE"); } DECLARE_PLUGIN(IInputHandler, GENIEInputHandler); diff --git a/src/generator/input/GENIEInputHandler.hxx b/src/generator/input/GENIEInputHandler.hxx index f1086f1..7c87f17 100644 --- a/src/generator/input/GENIEInputHandler.hxx +++ b/src/generator/input/GENIEInputHandler.hxx @@ -1,67 +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 "event/FullEvent.hxx" #include "input/IInputHandler.hxx" #include "exception/exception.hxx" #include "utility/ROOTUtility.hxx" #ifdef GENIE_V3_INTERFACE #include "Framework/Ntuple/NtpMCEventRecord.h" #else #include "Ntuple/NtpMCEventRecord.h" #endif #include namespace fhicl { class ParameterSet; } class GENIEInputHandler : public IInputHandler { mutable nuis::utility::TreeFile fInputTree; mutable nuis::event::FullEvent fReaderEvent; mutable std::vector fWeightCache; mutable genie::NtpMCEventRecord *fGenieNtpl; bool fKeepIntermediates; bool fKeepNuclearParticles; public: NEW_NUIS_EXCEPT(weight_cache_miss); GENIEInputHandler(); GENIEInputHandler(GENIEInputHandler const &) = delete; GENIEInputHandler(GENIEInputHandler &&); void Initialize(fhicl::ParameterSet 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(); + nuis::GeneratorManager::Generator_id_t GetGeneratorId() const; }; diff --git a/src/generator/input/NEUTInputHandler.cxx b/src/generator/input/NEUTInputHandler.cxx index 8968f3a..1d9ea8e 100644 --- a/src/generator/input/NEUTInputHandler.cxx +++ b/src/generator/input/NEUTInputHandler.cxx @@ -1,252 +1,254 @@ #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"); 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; } } 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() { +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 456b8cc..5418829 100644 --- a/src/generator/input/NEUTInputHandler.hxx +++ b/src/generator/input/NEUTInputHandler.hxx @@ -1,70 +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(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(); + nuis::GeneratorManager::Generator_id_t GetGeneratorId() const; }; diff --git a/src/generator/input/NuWroInputHandler.cxx b/src/generator/input/NuWroInputHandler.cxx index 0b50cb7..15e3814 100644 --- a/src/generator/input/NuWroInputHandler.cxx +++ b/src/generator/input/NuWroInputHandler.cxx @@ -1,132 +1,134 @@ #include "generator/input/NuWroInputHandler.hxx" #include "generator/utility/NuWroUtility.hxx" #include "fhiclcpp/ParameterSet.h" #include "particle.h" using NuWroParticle = ::particle; using namespace nuis; using namespace nuis::event; using namespace nuis::utility; using namespace nuis::nuwrotools; NuWroInputHandler::NuWroInputHandler() : fInputTree(), fTreeEvent(nullptr) {} NuWroInputHandler::NuWroInputHandler(NuWroInputHandler &&other) : fInputTree(std::move(other.fInputTree)), fReaderEvent(std::move(other.fReaderEvent)), fTreeEvent(other.fTreeEvent) {} void NuWroInputHandler::Initialize(fhicl::ParameterSet const &ps) { fInputTree = CheckGetTTree(ps.get("file"), "treeout"); fTreeEvent = nullptr; fInputTree.tree->SetBranchAddress("e", &fTreeEvent); fKeepIntermediates = ps.get("keep_intermediates", false); if(GetNEvents()){ fInputTree.tree->GetEntry(0); std::cout << "[INFO]: Average NuWro XSec weight: " << (fTreeEvent->weight / double(GetNEvents())) << std::endl; } } MinimalEvent const &NuWroInputHandler::GetMinimalEvent(ev_index_t idx) const { if (idx >= GetNEvents()) { throw IInputHandler::invalid_entry() << "[ERROR]: Attempted to get entry " << idx << " from an InputHandler with only " << GetNEvents(); } fInputTree.tree->GetEntry(idx); fReaderEvent.mode = NuWroEventChannel(*fTreeEvent); fReaderEvent.probe_E = fTreeEvent->in[0].E(); fReaderEvent.probe_pdg = fTreeEvent->in[0].pdg; fReaderEvent.XSecWeight = fTreeEvent->weight / double(GetNEvents()); if (fWeightCache.size() <= idx) { fWeightCache.push_back(fReaderEvent.XSecWeight); } return fReaderEvent; } FullEvent const &NuWroInputHandler::GetFullEvent(ev_index_t idx) const { (void)GetMinimalEvent(idx); fReaderEvent.ClearParticleStack(); for (size_t p_it = 0; p_it < fTreeEvent->in.size(); ++p_it) { NuWroParticle &part = fTreeEvent->in[p_it]; Particle nuis_part; nuis_part.pdg = part.pdg; nuis_part.P4 = TLorentzVector(part[1], part[2], part[3], part[0]); fReaderEvent .ParticleStack[static_cast( Particle::Status_t::kPrimaryInitialState)] .particles.push_back(nuis_part); } for (size_t p_it = 0; p_it < fKeepIntermediates && fTreeEvent->out.size(); ++p_it) { NuWroParticle &part = fTreeEvent->out[p_it]; Particle nuis_part; nuis_part.pdg = part.pdg; nuis_part.P4 = TLorentzVector(part[1], part[2], part[3], part[0]); fReaderEvent .ParticleStack[static_cast( Particle::Status_t::kPrimaryFinalState)] .particles.push_back(nuis_part); } for (size_t p_it = 0; (p_it < fTreeEvent->post.size()); ++p_it) { NuWroParticle &part = fTreeEvent->post[p_it]; Particle nuis_part; nuis_part.pdg = part.pdg; nuis_part.P4 = TLorentzVector(part[1], part[2], part[3], part[0]); fReaderEvent .ParticleStack[static_cast(Particle::Status_t::kNuclearLeaving)] .particles.push_back(nuis_part); } + fReaderEvent.fGenEvent = fTreeEvent; + return fReaderEvent; } NuWroEvent const &NuWroInputHandler::GetNuWroEvent(ev_index_t idx) const { (void)GetMinimalEvent(idx); return *fTreeEvent; } double NuWroInputHandler::GetEventWeight(ev_index_t idx) const { if (idx > fWeightCache.size()) { throw weight_cache_miss() << "[ERROR]: Failed to get cached weight for event index: " << idx; } return fWeightCache[idx]; } double NuWroInputHandler::GetXSecScaleFactor( std::pair const &EnuRange) const { throw input_handler_feature_unimplemented() << "[ERROR]: Flux cuts not yet implemented for NuWro input handler."; } size_t NuWroInputHandler::GetNEvents() const { return fInputTree.tree->GetEntries(); } -GeneratorManager::Generator_id_t NuWroInputHandler::GetGeneratorId(){ +GeneratorManager::Generator_id_t NuWroInputHandler::GetGeneratorId() const { return GeneratorManager::Get().EnsureGeneratorRegistered("NuWro"); } DECLARE_PLUGIN(IInputHandler, NuWroInputHandler); diff --git a/src/generator/input/NuWroInputHandler.hxx b/src/generator/input/NuWroInputHandler.hxx index d985f2d..ec36b0c 100644 --- a/src/generator/input/NuWroInputHandler.hxx +++ b/src/generator/input/NuWroInputHandler.hxx @@ -1,66 +1,66 @@ // 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 "exception/exception.hxx" #include "utility/ROOTUtility.hxx" #include "event1.h" using NuWroEvent = ::event; #include namespace fhicl { class ParameterSet; } class NuWroInputHandler : public IInputHandler { mutable nuis::utility::TreeFile fInputTree; mutable nuis::event::FullEvent fReaderEvent; mutable std::vector fWeightCache; mutable NuWroEvent *fTreeEvent; bool fKeepIntermediates; public: NEW_NUIS_EXCEPT(weight_cache_miss); NuWroInputHandler(); NuWroInputHandler(NuWroInputHandler const &) = delete; NuWroInputHandler(NuWroInputHandler &&); void Initialize(fhicl::ParameterSet const &); nuis::event::MinimalEvent const &GetMinimalEvent(ev_index_t idx) const; nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const; NuWroEvent const &GetNuWroEvent(ev_index_t idx) const; double GetEventWeight(ev_index_t idx) const; size_t GetNEvents() const; double GetXSecScaleFactor(std::pair const &EnuRange) const; - nuis::GeneratorManager::Generator_id_t GetGeneratorId(); + nuis::GeneratorManager::Generator_id_t GetGeneratorId() const; }; diff --git a/src/generator/utility/NEUTUtility.cxx b/src/generator/utility/NEUTUtility.cxx index b7f6310..5ff5b77 100644 --- a/src/generator/utility/NEUTUtility.cxx +++ b/src/generator/utility/NEUTUtility.cxx @@ -1,93 +1,95 @@ #include "generator/utility/NEUTUtility.hxx" #include "generator/input/NEUTInputHandler.hxx" #include "exception/exception.hxx" #include "utility/InteractionChannelUtility.hxx" #include "utility/PDGCodeUtility.hxx" #include "neutpart.h" #include "neutvect.h" using namespace nuis::event; using namespace nuis::utility; namespace nuis { namespace neuttools { NEW_NUIS_EXCEPT(unexpected_NEUT_particle_state); Particle::Status_t GetNeutParticleStatus(NeutPart const &part, Channel_t mode) { #ifdef DEBUG_NEUT_UTILITY std::cout << "[DEBUG]: Mode: " << mode << ", Part: { Status: " << part.fStatus << ", IsAlive: " << part.fIsAlive << ", PDG: " << part.fPID << " }." << std::endl; #endif // Remove Pauli blocked events, probably just single pion events if (part.fStatus == 5) { return Particle::Status_t::kBlocked; // fStatus == -1 means initial state } else if (part.fIsAlive == false && part.fStatus == -1) { return Particle::Status_t::kPrimaryInitialState; // NEUT has a bit of a strange convention for fIsAlive and fStatus // combinations // for NC and neutrino particle isAlive true/false and status 2 means // final state particle // for other particles in NC status 2 means it's an FSI particle // for CC it means it was an FSI particle } else if (part.fStatus == 2) { // NC case is a little strange... The outgoing neutrino might be alive or // not alive. Remaining particles with status 2 are FSI particles that // reinteracted if (IsNC(mode) && IsNeutralLepton(part.fPID)) { return Particle::Status_t::kNuclearLeaving; // The usual CC case } else if (part.fIsAlive == true) { // return Particle::Status_t::kIntermediate; throw unexpected_NEUT_particle_state() << "[ERROR] Found unexpected NEUT particle in neutvect stack: Mode: " - << mode << ", Part: { Status: " << part.fStatus - << ", IsAlive: " << part.fIsAlive << ", PDG: " << part.fPID << " }."; + << mode << " (IsNC: " << IsNC(mode) + << "), Part: { Status: " << part.fStatus + << ", IsAlive: " << part.fIsAlive << ", PDG: " << part.fPID + << ", IsNeutralLepton: " << IsNeutralLepton(part.fPID) << " }."; } } else if ((part.fIsAlive == true) && (part.fStatus == 2) && IsNeutralLepton(part.fPID)) { return Particle::Status_t::kNuclearLeaving; } else if ((part.fIsAlive == true) && (part.fStatus == 0)) { return Particle::Status_t::kNuclearLeaving; } else if (!part.fIsAlive && ((part.fStatus == 1) || (part.fStatus == 3) || (part.fStatus == 4) || (part.fStatus == 7) || (part.fStatus == 8))) { return Particle::Status_t::kIntermediate; // There's one hyper weird case where fStatus = -3. This apparently // corresponds to a nucleon being ejected via pion FSI when there is "data // available" } else if (!part.fIsAlive && (part.fStatus == -3)) { return Particle::Status_t::kUnknown; // NC neutrino outgoing } else if (!part.fIsAlive && part.fStatus == 0 && IsNeutralLepton(part.fPID)) { return Particle::Status_t::kNuclearLeaving; // Warn if we still find alive particles without classifying them } else if (part.fIsAlive == true) { throw unexpected_NEUT_particle_state() << "[ERROR]: Undefined NEUT state " << " Alive: " << part.fIsAlive << " Status: " << part.fStatus << " PDG: " << part.fPID; } // Warn if we find dead particles that we haven't classified throw unexpected_NEUT_particle_state() << "[ERROR]: Undefined NEUT state " << " Alive: " << part.fIsAlive << " Status: " << part.fStatus << " PDG: " << part.fPID; } } // namespace neuttools } // namespace nuis diff --git a/src/input/IInputHandler.hxx b/src/input/IInputHandler.hxx index cb0babd..fcb467e 100644 --- a/src/input/IInputHandler.hxx +++ b/src/input/IInputHandler.hxx @@ -1,122 +1,122 @@ // 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 "plugins/traits.hxx" #include "exception/exception.hxx" #include "generator/GeneratorManager.hxx" #include #include namespace fhicl { class ParameterSet; } namespace nuis { namespace event { class MinimalEvent; class FullEvent; } // namespace event } // namespace nuis class IInputHandler { public: struct FullEvent_const_iterator : public std::iterator< std::input_iterator_tag, nuis::event::FullEvent const, size_t, nuis::event::FullEvent const *, nuis::event::FullEvent const &> { FullEvent_const_iterator(size_t _idx, IInputHandler const *_ih) { idx = _idx; ih = _ih; } FullEvent_const_iterator(FullEvent_const_iterator const &other) { idx = other.idx; ih = other.ih; } FullEvent_const_iterator operator=(FullEvent_const_iterator const &other) { idx = other.idx; ih = other.ih; return (*this); } bool operator==(FullEvent_const_iterator const &other) { return (idx == other.idx); } bool operator!=(FullEvent_const_iterator const &other) { return !(*this == other); } nuis::event::FullEvent const &operator*() { return ih->GetFullEvent(idx); } nuis::event::FullEvent const *operator->() { return &ih->GetFullEvent(idx); } FullEvent_const_iterator operator++() { idx++; return *this; } FullEvent_const_iterator operator++(int) { FullEvent_const_iterator tmp(*this); idx++; return tmp; } private: size_t idx; IInputHandler const *ih; }; NEW_NUIS_EXCEPT(invalid_input_file); NEW_NUIS_EXCEPT(invalid_entry); NEW_NUIS_EXCEPT(input_handler_feature_unimplemented); typedef size_t ev_index_t; virtual void Initialize(fhicl::ParameterSet const &) = 0; virtual nuis::event::MinimalEvent const & GetMinimalEvent(ev_index_t idx) const = 0; virtual nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const = 0; virtual void RecalculateEventWeights(){}; virtual double GetEventWeight(ev_index_t) const { return 1; }; /// Allows samples that implement flux cuts to request an appropriate virtual double GetXSecScaleFactor( std::pair const &EnuRange = std::pair{ std::numeric_limits::max(), std::numeric_limits::max()}) const = 0; virtual size_t GetNEvents() const = 0; FullEvent_const_iterator begin() const { return FullEvent_const_iterator(0, this); } FullEvent_const_iterator end() const { return FullEvent_const_iterator(GetNEvents(), this); } virtual ~IInputHandler() {} - virtual nuis::GeneratorManager::Generator_id_t GetGeneratorId() = 0; + virtual nuis::GeneratorManager::Generator_id_t GetGeneratorId() const = 0; }; DECLARE_PLUGIN_INTERFACE(IInputHandler); diff --git a/src/utility/HistogramUtility.cxx b/src/utility/HistogramUtility.cxx index 9011cef..4246fe1 100644 --- a/src/utility/HistogramUtility.cxx +++ b/src/utility/HistogramUtility.cxx @@ -1,102 +1,152 @@ #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()) { + continue; + } + + 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); + } 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); + } + } + } +} + } // namespace utility } // namespace nuis diff --git a/src/utility/HistogramUtility.hxx b/src/utility/HistogramUtility.hxx index 6010309..c84d620 100644 --- a/src/utility/HistogramUtility.hxx +++ b/src/utility/HistogramUtility.hxx @@ -1,461 +1,467 @@ // 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 "utility/ROOTUtility.hxx" #include "exception/exception.hxx" #include "string_parsers/from_string.hxx" #include "fhiclcpp/ParameterSet.h" #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); bool IsFlowBin(TAxis const &ax, Int_t 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 void Clear(typename std::enable_if::NDim != 0, HT>::type &h) { for (Int_t bin_it = 0; bin_it < TH_Helper::NbinsIncludeFlow(h); ++bin_it) { h.SetBinContent(bin_it, 0); h.SetBinError(bin_it, 0); } } template void Clear( typename std::enable_if::value, HT>::type &h) { h.ClearBinContents(); } template inline std::unique_ptr Clone(HT const &source, bool clear = false, std::string const &clone_name = "") { std::unique_ptr target(dynamic_cast( source.Clone(clone_name.size() ? clone_name.c_str() : ""))); if (!target) { throw failed_to_clone() << "[ERROR]: Failed to clone a " << TH_Helper::name() << "."; } target->SetDirectory(nullptr); if (clear) { Clear(*target); } return target; } template inline std::unique_ptr Clone(std::unique_ptr const &source, bool clear = false, std::string const &clone_name = "") { return Clone(*source, clear, clone_name); } template inline std::unique_ptr GetHistogramFromROOTFile(TFile_ptr &f, std::string const &hname, bool ThrowIfMissing = true) { f->Get(hname.c_str()); HT *h = dynamic_cast(f->Get(hname.c_str())); if (!h) { if (ThrowIfMissing) { throw invalid_histogram_name() << "[ERROR]: Failed to get " << TH_Helper::name() << " named " << std::quoted(hname) << " from input file " << std::quoted(f->GetName()); } else { return nullptr; } } std::unique_ptr clone = Clone(*h); return clone; } template inline std::unique_ptr GetHistogramFromROOTFile(std::string const &fname, std::string const &hname, bool ThrowIfMissing = true) { TFile_ptr temp = CheckOpenTFile(fname, "READ"); return GetHistogramFromROOTFile(temp, hname, ThrowIfMissing); } template std::unique_ptr GetHistogram(std::string const &input_descriptor) { std::vector split_descriptor = fhicl::string_parsers::ParseToVect(input_descriptor, ";", true, true); if (split_descriptor.size() == 1) { // assume text throw unimplemented_GetHistogram_method(); } else if (split_descriptor.size() == 2) { return GetHistogramFromROOTFile(split_descriptor[0], split_descriptor[1]); } else { throw invalid_histogram_descriptor() << "[ERROR]: Failed to parse histogram descriptor: " << std::quoted(input_descriptor) << " as an input histogram (Text/ROOT)."; } } struct PolyBinSpecifier { double X, Y; bool UseXAxis; }; constexpr PolyBinSpecifier XPolyBinSpec(double X, double Y) { return PolyBinSpecifier{X + std::numeric_limits::epsilon() * 1E2, Y + std::numeric_limits::epsilon() * 1E2, true}; } constexpr PolyBinSpecifier YPolyBinSpec(double X, double Y) { return PolyBinSpecifier{X + std::numeric_limits::epsilon() * 1E2, Y + std::numeric_limits::epsilon() * 1E2, false}; } std::vector> GetTH2PolySlices( std::unique_ptr &hinp, std::vector> const &BinsSpecifiers); template inline typename std::enable_if::NDim == 1, std::unique_ptr>::type BuildHistFromFHiCL(fhicl::ParameterSet const &ps) { std::array xaxis = ps.get>("xaxis_descriptor"); std::unique_ptr rtn = std::make_unique( ps.get("name", "").c_str(), ps.get("title", "").c_str(), xaxis[0], xaxis[1], xaxis[2]); rtn->SetDirectory(nullptr); return rtn; } template inline typename std::enable_if::NDim == 2, std::unique_ptr>::type BuildHistFromFHiCL(fhicl::ParameterSet const &ps) { std::array xaxis = ps.get>("xaxis_descriptor"); std::array yaxis = ps.get>("yaxis_descriptor"); std::unique_ptr rtn = std::make_unique(ps.get("name", "").c_str(), ps.get("title", "").c_str(), xaxis[0], xaxis[1], xaxis[2], yaxis[0], yaxis[1], yaxis[2]); rtn->SetDirectory(nullptr); return rtn; } + +static bool const kYSlice = true; +static bool const kXSlice = false; +void SliceNorm(std::unique_ptr &hist, bool AlongY = kYSlice, + char const *opt = ""); + } // namespace utility } // namespace nuis diff --git a/src/utility/InteractionChannelUtility.hxx b/src/utility/InteractionChannelUtility.hxx index 1bced13..00b8a8c 100644 --- a/src/utility/InteractionChannelUtility.hxx +++ b/src/utility/InteractionChannelUtility.hxx @@ -1,74 +1,74 @@ #ifndef UTILITY_CHANNELUTILITY_HXX_SEEN #define UTILITY_CHANNELUTILITY_HXX_SEEN #include "event/types.hxx" #include "exception/exception.hxx" #include namespace nuis { namespace utility { NEW_NUIS_EXCEPT(invalid_channel_conversion); #define X(A, B) \ case B: { \ return nuis::event::Channel_t::A; \ } inline event::Channel_t IntToChannel(int mode) { switch (mode) { NUIS_INTERACTION_CHANNEL_LIST default: { throw invalid_channel_conversion() << "[ERROR]: Failed to parse mode integer " << mode << " as a nuis::event::Channel_t."; } } } #undef X #define X(A, B) \ case event::Channel_t::A: { \ return B; \ } inline int ChannelToInt(event::Channel_t mode) { switch (mode) { NUIS_INTERACTION_CHANNEL_LIST default: { throw invalid_channel_conversion() << "[ERROR]: Attempting to convert " "undefined nuis::event::Channel_t to an " "integer."; } } } #undef X -inline bool IsNC(event::Channel_t mode) { return abs(ChannelToInt(mode) > 30); } +inline bool IsNC(event::Channel_t mode) { return abs(ChannelToInt(mode)) > 30; } inline bool IsCC(event::Channel_t mode) { return !IsNC(mode); } inline bool IsNu(event::Channel_t mode) { return ChannelToInt(mode) > 0; } inline bool IsNub(event::Channel_t mode) { return !IsNu(mode); } inline bool IsCoh(event::Channel_t mode) { return ((mode == event::Channel_t::kCCCohPi) || (mode == event::Channel_t::kNCCohPi) || (mode == event::Channel_t::kCCCohPi_nub) || (mode == event::Channel_t::kNCCohPi_nub)); } inline bool IsQE(event::Channel_t mode) { return ((mode == event::Channel_t::kCCQE) || (mode == event::Channel_t::kNCELP) || (mode == event::Channel_t::kNCELN) || (mode == event::Channel_t::kCCQE_nub) || (mode == event::Channel_t::kNCELP_nub) || (mode == event::Channel_t::kNCELN_nub)); } } // namespace utility } // namespace nuis #endif