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 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)."; } } -NEW_NUIS_EXCEPT(invalid_TH2Poly); -NEW_NUIS_EXCEPT(invalid_PolyBinSpecifierList); - -// List of bins to put in a row, X/Y struct PolyBinSpecifier { double X, Y; bool UseXAxis; }; constexpr PolyBinSpecifier XPolyBinSpec(double X, double Y) { return PolyBinSpecifier{X + std::numeric_limits::epsilon() * 1E2, Y + std::numeric_limits::epsilon() * 1E2, true}; } constexpr PolyBinSpecifier YPolyBinSpec(double X, double Y) { return PolyBinSpecifier{X + std::numeric_limits::epsilon() * 1E2, Y + std::numeric_limits::epsilon() * 1E2, false}; } std::vector> GetTH2PolySlices( std::unique_ptr &hinp, - std::vector> const &BinsSpecifiers) { - - std::vector> slices; - - size_t sl_it = 0; - for (auto &slice_spec : BinsSpecifiers) { - std::vector Binning; - std::vector BinContent; - std::vector BinError; - bool UseXAxis = false; - size_t bin_ctr = 0; - for (auto poly_bin_spec : slice_spec) { - Int_t bin_it = hinp->FindBin(poly_bin_spec.X, poly_bin_spec.Y); - if (bin_it < 1) { - std::cout << "[WARN]: When searching for matching bin: { X: " - << poly_bin_spec.X << ", Y: " << poly_bin_spec.Y - << "} got flow bin: " << bin_it << std::endl; - continue; - } - TH2PolyBin *poly_bin = - dynamic_cast(hinp->GetBins()->At(bin_it - 1)); - - if (!bin_ctr) { - UseXAxis = poly_bin_spec.UseXAxis; - } else if (UseXAxis != poly_bin_spec.UseXAxis) { - throw invalid_PolyBinSpecifierList() - << "[ERROR]: For slice: " << sl_it - << " of TH2Poly: " << std::quoted(hinp->GetName()) - << " bin specifier: " << bin_ctr << " was set to use the " - << (poly_bin_spec.UseXAxis ? "X" : "Y") - << " axis as the dependent axis of the slice, but previous bins " - "were set to use the " - << (poly_bin_spec.UseXAxis ? "X" : "Y") << " axis."; - } - - if ((poly_bin_spec.X < poly_bin->GetXMin()) || - (poly_bin_spec.X >= poly_bin->GetXMax()) || - (poly_bin_spec.Y < poly_bin->GetYMin()) || - (poly_bin_spec.Y >= poly_bin->GetYMax())) { - std::cout - << "[WARN]: Found bin doesn't seem to contain expected point: { X: " - << poly_bin_spec.X << ", Y: " << poly_bin_spec.Y - << "}, got bin_it = " << bin_it - << " which had x_low: " << poly_bin->GetXMin() - << ", and x_up: " << poly_bin->GetXMax() - << ", y_low: " << poly_bin->GetYMin() - << ", y_up: " << poly_bin->GetYMax() << std::endl; - } - - double low = UseXAxis ? poly_bin->GetXMin() : poly_bin->GetYMin(); - double up = UseXAxis ? poly_bin->GetXMax() : poly_bin->GetYMax(); - - if (!Binning.size()) { // Add low edge - Binning.push_back(low); - } else if (std::abs(Binning.back() - low) > - (std::numeric_limits::epsilon() * 1E2)) { - BinContent.push_back(0); - BinError.push_back(0); - Binning.push_back(low); - } - - BinContent.push_back(hinp->GetBinContent(bin_it)); - BinError.push_back(hinp->GetBinError(bin_it)); - Binning.push_back(up); - } - slices.emplace_back(new TH1D( - (std::string(hinp->GetName()) + "_slice" + std::to_string(sl_it++)) - .c_str(), - (std::string(";") + - (UseXAxis ? hinp->GetXaxis() : hinp->GetYaxis())->GetTitle() + ";" + - hinp->GetZaxis()->GetTitle()) - .c_str(), - Binning.size() - 1, Binning.data())); - - for (size_t bin_it = 0; bin_it < BinContent.size(); ++bin_it) { - slices.back()->SetBinContent(bin_it + 1, BinContent[bin_it]); - slices.back()->SetBinError(bin_it + 1, BinError[bin_it]); - } - } - return slices; -} - + std::vector> const &BinsSpecifiers); } // namespace utility } // namespace nuis -#endif diff --git a/src/utility/KinematicUtility.cxx b/src/utility/KinematicUtility.cxx index 8fda873..cfdbd69 100644 --- a/src/utility/KinematicUtility.cxx +++ b/src/utility/KinematicUtility.cxx @@ -1,106 +1,228 @@ #include "utility/KinematicUtility.hxx" #include "event/FullEvent.hxx" #include "event/Particle.hxx" #include "utility/FullEventUtility.hxx" #include "utility/PDGCodeUtility.hxx" using namespace nuis::event; namespace nuis { namespace utility { double GetNeutrinoEQERec(FullEvent const &fev, double SeparationEnergy_MeV) { Particle const &neutrino = GetHMISNeutralLepton(fev); if (!neutrino) { throw Particle::invalid_particle() << "[ERROR]: In GetNeutrinoEQERec, expected to be able to get IS " "neutral lepton, but found none:" << "\n" << fev.to_string(); } Particle const &charged_lepton = GetHMFSChargedLepton(fev); if (!charged_lepton) { throw Particle::invalid_particle() << "[ERROR]: In GetNeutrinoEQERec, expected to be able to get FS " "charged lepton, but found none:" << "\n" << fev.to_string(); } double mass_ISN = (IsMatter(neutrino.pdg) ? (pdgmasses::kNeutronMass_MeV - SeparationEnergy_MeV) : (pdgmasses::kProtonMass_MeV - SeparationEnergy_MeV)) / 1000.0; double mass_FSN = (IsMatter(neutrino.pdg) ? pdgmasses::kProtonMass_MeV : pdgmasses::kNeutronMass_MeV) / 1000.0; double el_GeV = charged_lepton.E() / 1000.0; double pl_GeV = charged_lepton.P() / 1000.0; // momentum of lepton double ml_GeV = charged_lepton.M() / 1000.0; // lepton mass double Theta_nu_mu = neutrino.P3().Angle(charged_lepton.P3()); return (2.0 * mass_ISN * el_GeV - ml_GeV * ml_GeV + mass_FSN * mass_FSN - mass_ISN * mass_ISN) / (2.0 * (mass_ISN - el_GeV + pl_GeV * cos(Theta_nu_mu))); } double GetNeutrinoQ2QERec(FullEvent const &fev, double SeparationEnergy_MeV) { Particle const &neutrino = GetHMISNeutralLepton(fev); if (!neutrino) { throw Particle::invalid_particle() << "[ERROR]: In GetNeutrinoQ2QERec, expected to be able to get IS " "neutral lepton, but found none:" << "\n" << fev.to_string(); } Particle const &charged_lepton = GetHMFSChargedLepton(fev); if (!charged_lepton) { throw Particle::invalid_particle() << "[ERROR]: In GetNeutrinoQ2QERec, expected to be able to get FS " "charged lepton, but found none:" << "\n" << fev.to_string(); } double el_GeV = charged_lepton.E() / 1000.0; double pl_GeV = charged_lepton.P() / 1000.0; // momentum of lepton double ml_GeV = charged_lepton.M() / 1000.0; // lepton mass double Theta_nu_mu = neutrino.P3().Angle(charged_lepton.P3()); return -ml_GeV * ml_GeV + 2.0 * GetNeutrinoEQERec(fev, SeparationEnergy_MeV) * (el_GeV - pl_GeV * cos(Theta_nu_mu)); } +TVector3 GetVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal) { + TVector3 pnUnit = planarNormal.Unit(); + double inpProjectPN = inp.Dot(pnUnit); + + TVector3 InPlanarInput = inp - (inpProjectPN * pnUnit); + return InPlanarInput; +} + +TVector3 GetUnitVectorInTPlane(const TVector3 &inp, + const TVector3 &planarNormal) { + return GetVectorInTPlane(inp, planarNormal).Unit(); +} + +double GetDeltaPhiT(TVector3 const &V_lepton, TVector3 const &V_other, + TVector3 const &Normal, bool PiMinus) { + TVector3 V_lepton_T = GetUnitVectorInTPlane(V_lepton, Normal); + + TVector3 V_other_T = GetUnitVectorInTPlane(V_other, Normal); + + return PiMinus ? acos(V_lepton_T.Dot(V_other_T)) + : (M_PI - acos(V_lepton_T.Dot(V_other_T))); +} + +TVector3 GetDeltaPT(TVector3 const &V_lepton, TVector3 const &V_other, + TVector3 const &Normal) { + TVector3 V_lepton_T = GetVectorInTPlane(V_lepton, Normal); + + TVector3 V_other_T = GetVectorInTPlane(V_other, Normal); + + return V_lepton_T + V_other_T; +} + +double GetDeltaAlphaT(TVector3 const &V_lepton, TVector3 const &V_other, + TVector3 const &Normal, bool PiMinus) { + TVector3 DeltaPT = GetDeltaPT(V_lepton, V_other, Normal); + + return GetDeltaPhiT(V_lepton, DeltaPT, Normal, PiMinus); +} + +TVector3 GetDeltaPT_CC0PiN(event::FullEvent const &fev, event::PDG_t islep_pdg, + event::PDG_t fslep_pdg, event::PDG_t fsnuc_pdg) { + if (fslep_pdg == nuis::utility::pdgcodes::kDefault) { // Assume neutrino CC + fslep_pdg = islep_pdg > 0 ? islep_pdg - 1 : islep_pdg + 1; + } + if (fsnuc_pdg == nuis::utility::pdgcodes::kDefault) { + fsnuc_pdg = islep_pdg > 0 ? nuis::utility::pdgcodes::kProton + : nuis::utility::pdgcodes::kNeutron; + } + + Particle ISSLep = GetHMISParticle(fev, {islep_pdg}); + if (!ISSLep) { + return {-std::numeric_limits::max(), + -std::numeric_limits::max(), + -std::numeric_limits::max()}; + } + Particle FSLep = GetHMFSParticle(fev, {fslep_pdg}); + if (!FSLep) { + return {-std::numeric_limits::max(), + -std::numeric_limits::max(), + -std::numeric_limits::max()}; + } + Particle FSNuc = GetHMFSParticle(fev, {fsnuc_pdg}); + if (!FSNuc) { + return {-std::numeric_limits::max(), + -std::numeric_limits::max(), + -std::numeric_limits::max()}; + } + + return GetDeltaPT(FSLep.P3(), FSNuc.P3(), ISSLep.P3()); +} +double GetDeltaPhiT_CC0PiN(event::FullEvent const &fev, event::PDG_t islep_pdg, + event::PDG_t fslep_pdg, event::PDG_t fsnuc_pdg) { + if (fslep_pdg == nuis::utility::pdgcodes::kDefault) { // Assume neutrino CC + fslep_pdg = islep_pdg > 0 ? islep_pdg - 1 : islep_pdg + 1; + } + if (fsnuc_pdg == nuis::utility::pdgcodes::kDefault) { + fsnuc_pdg = islep_pdg > 0 ? nuis::utility::pdgcodes::kProton + : nuis::utility::pdgcodes::kNeutron; + } + + Particle ISSLep = GetHMISParticle(fev, {islep_pdg}); + if (!ISSLep) { + return -std::numeric_limits::max(); + } + Particle FSLep = GetHMFSParticle(fev, {fslep_pdg}); + if (!FSLep) { + return -std::numeric_limits::max(); + } + Particle FSNuc = GetHMFSParticle(fev, {fsnuc_pdg}); + if (!FSNuc) { + return -std::numeric_limits::max(); + } + + return GetDeltaPhiT(FSLep.P3(), FSNuc.P3(), ISSLep.P3()); +} +double GetDeltaAlphaT_CC0PiN(event::FullEvent const &fev, + event::PDG_t islep_pdg, event::PDG_t fslep_pdg, + event::PDG_t fsnuc_pdg) { + if (fslep_pdg == nuis::utility::pdgcodes::kDefault) { // Assume neutrino CC + fslep_pdg = islep_pdg > 0 ? islep_pdg - 1 : islep_pdg + 1; + } + if (fsnuc_pdg == nuis::utility::pdgcodes::kDefault) { + fsnuc_pdg = islep_pdg > 0 ? nuis::utility::pdgcodes::kProton + : nuis::utility::pdgcodes::kNeutron; + } + + Particle ISSLep = GetHMISParticle(fev, {islep_pdg}); + if (!ISSLep) { + return -std::numeric_limits::max(); + } + Particle FSLep = GetHMFSParticle(fev, {fslep_pdg}); + if (!FSLep) { + return -std::numeric_limits::max(); + } + Particle FSNuc = GetHMFSParticle(fev, {fsnuc_pdg}); + if (!FSNuc) { + return -std::numeric_limits::max(); + } + + return GetDeltaAlphaT(FSLep.P3(), FSNuc.P3(), ISSLep.P3()); +} + TLorentzVector GetEnergyMomentumTransfer(event::FullEvent const &fev) { Particle const &neutrino = GetHMISNeutralLepton(fev); if (!neutrino) { throw Particle::invalid_particle() << "[ERROR]: In GetEnergyMomentumTransfer, expected to be able to get " "IS neutral lepton, but found none: \n" << fev.to_string(); } Particle const &charged_lepton = GetHMFSChargedLepton(fev); if (!charged_lepton) { throw Particle::invalid_particle() << "[ERROR]: In GetEnergyMomentumTransfer, expected to be able to get " "FS charged lepton, but found none: \n" << fev.to_string(); } - return (neutrino.P4 - charged_lepton.P4); + return (neutrino.P4 - charged_lepton.P4); } } // namespace utility } // namespace nuis diff --git a/src/utility/KinematicUtility.hxx b/src/utility/KinematicUtility.hxx index f44547a..52c0536 100644 --- a/src/utility/KinematicUtility.hxx +++ b/src/utility/KinematicUtility.hxx @@ -1,57 +1,86 @@ // 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_KINEMATICUTILITY_HXX_SEEN #define UTILITY_KINEMATICUTILITY_HXX_SEEN +#include "utility/PDGCodeUtility.hxx" + #include "TLorentzVector.h" +#include "TVector3.h" #include #include namespace nuis { namespace event { class FullEvent; } } // namespace nuis namespace nuis { namespace utility { double GetNeutrinoEQERec(event::FullEvent const &fev, double SeparationEnergy_MeV); double GetNeutrinoQ2QERec(event::FullEvent const &fev, double SeparationEnergy_MeV); +TVector3 GetVectorInTPlane(const TVector3 &inp, const TVector3 &planarNormal); +TVector3 GetUnitVectorInTPlane(const TVector3 &inp, + const TVector3 &planarNormal); +double GetDeltaPhiT(TVector3 const &V_lepton, TVector3 const &V_other, + TVector3 const &Normal, bool PiMinus = false); +TVector3 GetDeltaPT(TVector3 const &V_lepton, TVector3 const &V_other, + TVector3 const &Normal); +double GetDeltaAlphaT(TVector3 const &V_lepton, TVector3 const &V_other, + TVector3 const &Normal, bool PiMinus = false); + +TVector3 +GetDeltaPT_CC0PiN(event::FullEvent const &fev, + event::PDG_t islep_pdg = nuis::utility::pdgcodes::kNuMu, + event::PDG_t fslep_pdg = nuis::utility::pdgcodes::kDefault, + event::PDG_t fsnuc_pdg = nuis::utility::pdgcodes::kDefault); +double +GetDeltaPhiT_CC0PiN(event::FullEvent const &fev, + event::PDG_t islep_pdg = nuis::utility::pdgcodes::kNuMu, + event::PDG_t fslep_pdg = nuis::utility::pdgcodes::kDefault, + event::PDG_t fsnuc_pdg = nuis::utility::pdgcodes::kDefault); +double GetDeltaAlphaT_CC0PiN( + event::FullEvent const &fev, + event::PDG_t islep_pdg = nuis::utility::pdgcodes::kNuMu, + event::PDG_t fslep_pdg = nuis::utility::pdgcodes::kDefault, + event::PDG_t fsnuc_pdg = nuis::utility::pdgcodes::kDefault); + TLorentzVector GetEnergyMomentumTransfer(event::FullEvent const &fev); struct ENuRange : public std::pair { using std::pair::pair; using std::pair::operator=; ENuRange() : std::pair(0, std::numeric_limits::max()) {} bool IsInRange(double enu) { return (enu > second) || (enu < first); } }; } // namespace utility } // namespace nuis #endif diff --git a/src/utility/PDGCodeUtility.hxx b/src/utility/PDGCodeUtility.hxx index e76aa7f..3cf290c 100644 --- a/src/utility/PDGCodeUtility.hxx +++ b/src/utility/PDGCodeUtility.hxx @@ -1,131 +1,133 @@ // 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/types.hxx" #include "exception/exception.hxx" #include namespace nuis { namespace utility { namespace pdgcodes { NEW_NUIS_EXCEPT(unhandled_pdg_code); enum MatterType { kMatter = 1, kMatterAntimatter = 0, kAntimatter = -1 }; +static event::PDG_t const kDefault = 0; + static event::PDG_t const kNuMu = 14; static event::PDG_t const kNuMuBar = -14; static event::PDG_t const kMu = 13; static event::PDG_t const kMuPlus = -13; static event::PDG_t const kElectron = 11; static event::PDG_t const kPositron = -11; static event::PDG_t const kNue = 12; static event::PDG_t const kNueBar = 12; static event::PDG_t const kPiPlus = 211; static event::PDG_t const kPiMinus = -211; static event::PDG_t const kPi0 = 111; static event::PDG_t const kProton = 2212; static event::PDG_t const kNeutron = 2112; static std::vector const ChargedLeptons{kElectron, kMu, 15, kPositron, kMuPlus, -15}; static std::vector const ChargedLeptons_matter{kElectron, kMu, 15}; static std::vector const ChargedLeptons_antimatter{kPositron, kMuPlus, -15}; static std::vector const NeutralLeptons{kNue, kNuMu, 16, kNueBar, kNuMuBar, -16}; static std::vector const NeutralLeptons_matter{kNue, kNuMu, 16}; static std::vector const NeutralLeptons_antimatter{kNueBar, kNuMuBar, -16}; static std::vector const ChargedPions{kPiPlus, kPiMinus}; static std::vector const NeutralPions{kPi0}; static std::vector const Pions{kPiPlus, kPiMinus, kPi0}; static std::vector const Protons{kProton, -kProton}; static std::vector const Proton_matter{kProton}; static std::vector const Proton_antimatter{-kProton}; static std::vector const Neutron{kNeutron}; static std::vector const Nucleons{kProton, kNeutron, -kProton}; static std::vector const Nucleons_matter{kProton, kNeutron}; static std::vector const Nucleons_antimatter{-kProton, kNeutron}; static std::vector const CommonParticles{ kElectron, kMu, 15, kPositron, kMuPlus, -15, kNue, kNuMu, 16, kNueBar, kNuMuBar, -16, kPiPlus, kPiMinus, kPi0, kProton, kNeutron}; } // namespace pdgcodes namespace pdgmasses { static double const kNuMuMass_MeV = 0; static double const kNuMuBarMass_MeV = 0; static double const kMuMass_MeV = 105.65; static double const kMuPlusMass_MeV = 105.65; static double const kNueMass_MeV = 0; static double const kNueBarMass_MeV = 0; static double const kPiPlusMass_MeV = 139.57; static double const kPiMinusMass_MeV = 139.57; static double const kPi0Mass_MeV = 134.97; static double const kNeutronMass_MeV = 939.56; static double const kProtonMass_MeV = 938.27; } // namespace pdgmasses double GetPDGMass(event::PDG_t); bool IsInPDGList(event::PDG_t pdg, std::vector const &MatterList, std::vector const &AntiMatterList, pdgcodes::MatterType type = pdgcodes::kMatterAntimatter); bool IsNeutralLepton(event::PDG_t pdg, pdgcodes::MatterType type = pdgcodes::kMatterAntimatter); bool IsChargedLepton(event::PDG_t pdg, pdgcodes::MatterType type = pdgcodes::kMatterAntimatter); bool IsProton(event::PDG_t pdg, pdgcodes::MatterType type = pdgcodes::kMatterAntimatter); bool IsNeutron(event::PDG_t pdg); bool IsChargedPion(event::PDG_t pdg); bool IsNeutralPion(event::PDG_t pdg); bool IsPion(event::PDG_t pdg); bool IsOther(event::PDG_t pdg); bool IsMatter(event::PDG_t pdg); bool IsAntiMatter(event::PDG_t pdg); bool IsNuclearPDG(event::PDG_t pdg); event::PDG_t MakeNuclearPDG(size_t A, size_t Z); size_t GetA(event::PDG_t); size_t GetZ(event::PDG_t); } // namespace utility } // namespace nuis