diff --git a/src/app/nuiscomp.cxx b/src/app/nuiscomp.cxx index cfb36c4..da163f3 100644 --- a/src/app/nuiscomp.cxx +++ b/src/app/nuiscomp.cxx @@ -1,51 +1,53 @@ #include "config/GlobalConfiguration.hxx" #include "input/IInputHandler.hxx" #include "event/MinimalEvent.hxx" -#include "samples/ISample.hxx" +#include "samples/IDataComparison.hxx" #include "plugins/Instantiate.hxx" #include "exception/exception.hxx" #include "fhiclcpp/make_ParameterSet.h" #include NEW_NUIS_EXCEPT(invalid_cli_arguments); int main(int argc, char const *argv[]) { nuis::config::EnsureConfigurationRead("nuis.global.config.fcl"); nuis::config::EnsureConfigurationRead("nuis.datacomparisons.fcl"); if (argc != 2) { throw invalid_cli_arguments() << "[ERROR]: Expected to be passed a single FHiCL file name or " "absolute or relative path. N.B. Files in the local directory must " "be fully qualified like \"$ " << argv[0] << " ./myconf.fcl\"."; } nuis::config::EnsureConfigurationRead(argv[1]); size_t NMax = nuis::config::GetDocument().get( "nmax", std::numeric_limits::max()); for (fhicl::ParameterSet const &samp_config : nuis::config::GetDocument().get>( "samples")) { std::cout << "[INFO]: Reading sample: " << samp_config.get("name") << std::endl; - nuis::plugins::plugin_traits::unique_ptr_t sample = - nuis::plugins::Instantiate( + nuis::plugins::plugin_traits::unique_ptr_t sample = + nuis::plugins::Instantiate( samp_config.get("name")); sample->Initialize(samp_config); sample->ProcessSample(NMax); + std::cout << "[INFO]:\t Sample GOF = " << sample->GetGOF() << " / " + << sample->GetNDOGuess() << std::endl; sample->Write(); } } diff --git a/src/app/nuissamples.cxx b/src/app/nuissamples.cxx index d61d1e3..3fdf60d 100644 --- a/src/app/nuissamples.cxx +++ b/src/app/nuissamples.cxx @@ -1,36 +1,111 @@ #include "config/GlobalConfiguration.hxx" #include "input/IInputHandler.hxx" #include "plugins/Instantiate.hxx" #include "samples/IDataComparison.hxx" +#include "utility/StringUtility.hxx" + +#include #include +#include #include -int main() { +std::string search_term = ".*"; +bool fStrictRegex = false; +std::string config_out_filename = ""; + +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; +} + +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]; + } else if (std::string(argv[opt]) == "-S") { + fStrictRegex = true; + search_term = argv[++opt]; + } else if (std::string(argv[opt]) == "-o") { + config_out_filename = argv[++opt]; + } 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::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)) { + continue; + } + nuis::plugins::plugin_traits::unique_ptr_t sample = nuis::plugins::Instantiate(sample_name); - + 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: " << sample->GetDocumentation() << std::endl; + std::cout << "\tDocs: \n" + << nuis::utility::indent_apply_width(sample->GetDocumentation(), + 10) + << std::endl; + std::cout << "Config: " + << sample->GetExampleConfiguration().to_indented_string() + << 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/config/GlobalConfiguration.cxx b/src/config/GlobalConfiguration.cxx index c178fcd..ec35a15 100644 --- a/src/config/GlobalConfiguration.cxx +++ b/src/config/GlobalConfiguration.cxx @@ -1,40 +1,41 @@ #include "config/GlobalConfiguration.hxx" #include "fhiclcpp/ParameterSet.h" #include "fhiclcpp/make_ParameterSet.h" namespace nuis { namespace config { struct NamedDoc { std::string name; fhicl::ParameterSet document; std::vector files_read; }; std::vector Configurations; NamedDoc &GetDocument_modifyable(std::string const &name) { for (NamedDoc &doc : Configurations) { if (doc.name == name) { return doc; } } Configurations.push_back(NamedDoc{name, fhicl::ParameterSet()}); return GetDocument_modifyable(name); } fhicl::ParameterSet const &GetDocument(std::string const &name) { return GetDocument_modifyable(name).document; } bool EnsureConfigurationRead(std::string const &fhicl_file, std::string const &name) { NamedDoc &doc = GetDocument_modifyable(name); if (std::count(doc.files_read.begin(), doc.files_read.end(), fhicl_file)) { return false; } doc.document.splice(fhicl::make_ParameterSet(fhicl_file)); return true; } + } // namespace config } // namespace nuis diff --git a/src/config/GlobalConfiguration.hxx b/src/config/GlobalConfiguration.hxx index c5a75dd..ffab64c 100644 --- a/src/config/GlobalConfiguration.hxx +++ b/src/config/GlobalConfiguration.hxx @@ -1,37 +1,38 @@ // 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 CONFIG_GLOBALCONFIGURATION_HXX_SEEN #define CONFIG_GLOBALCONFIGURATION_HXX_SEEN #include namespace fhicl { class ParameterSet; } namespace nuis { namespace config { fhicl::ParameterSet const &GetDocument(std::string const &name = "global"); bool EnsureConfigurationRead(std::string const &fhicl_file, std::string const &name = "global"); + } // namespace config } // namespace nuis #endif diff --git a/src/input/IInputHandler.hxx b/src/input/IInputHandler.hxx index f0a2890..e558904 100644 --- a/src/input/IInputHandler.hxx +++ b/src/input/IInputHandler.hxx @@ -1,113 +1,113 @@ // 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 CORE_IINPUTHANDLER_HXX_SEEN -#define CORE_IINPUTHANDLER_HXX_SEEN +#ifndef INPUT_IINPUTHANDLER_HXX_SEEN +#define INPUT_IINPUTHANDLER_HXX_SEEN #include "plugins/traits.hxx" #include "exception/exception.hxx" #include namespace fhicl { class ParameterSet; } namespace nuis { namespace event { class MinimalEvent; class FullEvent; } // namespace event } // namespace nuis class IInputHandler { public: struct FullEvent_const_iterator : public std::iterator< std::input_iterator_tag, nuis::event::FullEvent const, size_t, nuis::event::FullEvent const *, nuis::event::FullEvent const &> { FullEvent_const_iterator(size_t _idx, IInputHandler const *_ih) { idx = _idx; ih = _ih; } FullEvent_const_iterator(FullEvent_const_iterator const &other) { idx = other.idx; ih = other.ih; } FullEvent_const_iterator operator=(FullEvent_const_iterator const &other) { idx = other.idx; ih = other.ih; return (*this); } bool operator==(FullEvent_const_iterator const &other) { return (idx == other.idx); } bool operator!=(FullEvent_const_iterator const &other) { return !(*this == other); } nuis::event::FullEvent const &operator*() { return ih->GetFullEvent(idx); } nuis::event::FullEvent const *operator->() { return &ih->GetFullEvent(idx); } FullEvent_const_iterator operator++() { idx++; return *this; } FullEvent_const_iterator operator++(int) { FullEvent_const_iterator tmp(*this); idx++; return tmp; } private: size_t idx; IInputHandler const *ih; }; NEW_NUIS_EXCEPT(invalid_input_file); NEW_NUIS_EXCEPT(invalid_entry); typedef size_t ev_index_t; virtual void Initialize(fhicl::ParameterSet const &) = 0; virtual nuis::event::MinimalEvent const & GetMinimalEvent(ev_index_t idx) const = 0; virtual nuis::event::FullEvent const &GetFullEvent(ev_index_t idx) const = 0; virtual void RecalculateEventWeights(){}; virtual double GetEventWeight(ev_index_t idx) const { return 1; }; virtual size_t GetNEvents() const = 0; FullEvent_const_iterator begin() const { return FullEvent_const_iterator(0, this); } FullEvent_const_iterator end() const { return FullEvent_const_iterator(GetNEvents(), this); } virtual ~IInputHandler() {} }; DECLARE_PLUGIN_INTERFACE(IInputHandler); #endif diff --git a/src/input/InputManager.hxx b/src/input/InputManager.hxx index 0a119f5..80a421d 100644 --- a/src/input/InputManager.hxx +++ b/src/input/InputManager.hxx @@ -1,64 +1,64 @@ // Copyright 2018 L. Pickering, P Stowell, R. Terri, C. Wilkinson, C. Wret /******************************************************************************* * This file is part of NUISANCE. * * NUISANCE is free software: you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * NUISANCE is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with NUISANCE. If not, see . *******************************************************************************/ -#ifndef CORE_INPUTMANAGER_HXX_SEEN -#define CORE_INPUTMANAGER_HXX_SEEN +#ifndef INPUT_INPUTMANAGER_HXX_SEEN +#define INPUT_INPUTMANAGER_HXX_SEEN #include "input/IInputHandler.hxx" #include "plugins/traits.hxx" #include "exception/exception.hxx" #include #include namespace fhicl { class ParameterSet; } namespace nuis { namespace input { class InputManager { struct NamedInputHandler { NamedInputHandler(std::string const &, plugins::plugin_traits::unique_ptr_t &&); std::string name; plugins::plugin_traits::unique_ptr_t handler; }; std::vector Inputs; InputManager(); static InputManager *_global_inst; public: static InputManager &Get(); NEW_NUIS_EXCEPT(unknown_input); typedef size_t Input_id_t; Input_id_t EnsureInputLoaded(fhicl::ParameterSet const &); Input_id_t GetInputId(std::string const &) const; IInputHandler const &GetInputHandler(Input_id_t) const; }; } // namespace input } // namespace nuis #endif diff --git a/src/samples/SimpleDataComparison.hxx b/src/samples/SimpleDataComparison.hxx index 15559bd..bdce072 100644 --- a/src/samples/SimpleDataComparison.hxx +++ b/src/samples/SimpleDataComparison.hxx @@ -1,346 +1,375 @@ // 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_SIMPLEDATACOMPARISON_HXX_SEEN #define SAMPLES_SIMPLEDATACOMPARISON_HXX_SEEN #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/ROOTUtility.hxx" +#include "utility/StatsUtility.hxx" #include #include #include #include #include template struct TH_dim_helper {}; template <> struct TH_dim_helper<1> { typedef TH1D type; }; template <> struct TH_dim_helper<2> { typedef TH2D type; }; template class SimpleDataComparison : public IDataComparison { NEW_NUIS_EXCEPT(invalid_SimpleDataComparison_initialization); + NEW_NUIS_EXCEPT(SimpleDataComparison_already_finalized); protected: 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::type> fData; + std::string fMaskInputDescriptor; + std::unique_ptr::type> fMask; std::string fCovarianceInputDescriptor; std::unique_ptr fCovariance; std::unique_ptr::type> fPrediction; std::unique_ptr::type> fPrediction_xsec; std::unique_ptr::type> fPrediction_shape; std::unique_ptr::type> fPrediction_comparison; bool fComparisonFinalized; std::string fJournalReference; std::string fTargetMaterial; std::string fFluxDescription; std::string fSignalDescription; std::pair EnuRange; std::function IsSigFunc; std::function(nuis::event::FullEvent const &)> CompProjFunc; 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 (double &el : arr) { el = 0; } return arr; }; fJournalReference = ""; fTargetMaterial = ""; fFluxDescription = ""; fSignalDescription = ""; fIsShapeOnly = -1; fIsFluxUnfolded = -1; EnuRange = std::pair{std::numeric_limits::max(), std::numeric_limits::max()}; } void ReadGlobalConfigDefaults() { fhicl::ParameterSet const &global_sample_configuration = nuis::config::GetDocument().get( std::string("global.sample_configuration.") + Name(), fhicl::ParameterSet()); if (!fJournalReference.length()) { fJournalReference = global_sample_configuration.get( "journal_reference", fJournalReference); } if (!fTargetMaterial.length()) { fTargetMaterial = global_sample_configuration.get( "target_material", fTargetMaterial); } if (!fFluxDescription.length()) { fFluxDescription = global_sample_configuration.get( "flux_description", fFluxDescription); } if (!fSignalDescription.length()) { fSignalDescription = global_sample_configuration.get( "signal_description", fSignalDescription); } if (fIsShapeOnly == -1) { fIsShapeOnly = global_sample_configuration.get("shape_only", false); } if (fIsFluxUnfolded == -1) { fIsFluxUnfolded = global_sample_configuration.get("flux_unfolded", false); } if ((EnuRange.first == std::numeric_limits::max()) && (global_sample_configuration.has_key("enu_range"))) { EnuRange = global_sample_configuration.get>( "enu_range"); } } virtual std::string GetJournalReference() { return fJournalReference; } 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, double event_weight) { nuis::utility::Fill(fPrediction.get(), proj, event_weight); } virtual void FinalizeComparison() { + if (fComparisonFinalized) { + throw SimpleDataComparison_already_finalized() + << "[ERROR]: Attempted to re-finalize a comparison for " + "IDataComparison: " + << std::quoted(Name()); + } fPrediction_xsec = nuis::utility::Clone(fPrediction); fPrediction_xsec->Scale(1.0, "width"); fPrediction_shape = nuis::utility::Clone(fPrediction_xsec); if (fData) { fPrediction_shape->Scale(fData->Integral() / fPrediction_shape->Integral()); } else { ISAMPLE_WARN("When Finalizing comparison, no Data histogram available."); } if (fIsFluxUnfolded) { // fPrediction_comparison } else if (fIsShapeOnly) { fPrediction_comparison = nuis::utility::Clone(fPrediction_shape); } else { fPrediction_comparison = nuis::utility::Clone(fPrediction_xsec); } fComparisonFinalized = true; } void Initialize(fhicl::ParameterSet const &ps) { if (ps.has_key("verbosity")) { SetSampleVerbosity(ps.get("verbosity")); } else { SetSampleVerbosity("Reticent"); } ReadGlobalConfigDefaults(); fhicl::ParameterSet const &global_sample_configuration = nuis::config::GetDocument().get( std::string("global.sample_configuration.") + Name(), fhicl::ParameterSet()); if (ps.has_key("fake_data")) { fData = nuis::utility::GetHistogram::type>( ps.get("fake_data_histogram")); } else { if (!fDataInputDescriptor.length()) { if (!global_sample_configuration.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 = global_sample_configuration.get("data_descriptor"); } fData = nuis::utility::GetHistogram::type>( fDataInputDescriptor); } fPrediction = nuis::utility::Clone(fData, true); if (fCovarianceInputDescriptor.length()) { fCovariance = nuis::utility::GetHistogram(fCovarianceInputDescriptor); } else if (global_sample_configuration.has_key("covariance_descriptor")) { fCovariance = nuis::utility::GetHistogram( global_sample_configuration.get( "covariance_descriptor")); } + if (fMaskInputDescriptor.length()) { + fMask = nuis::utility::GetHistogram::type>( + fMaskInputDescriptor); + } else if (global_sample_configuration.has_key("mask_descriptor")) { + fMask = nuis::utility::GetHistogram::type>( + global_sample_configuration.get("mask_descriptor")); + } + if (ps.has_key("verbosity")) { SetSampleVerbosity(ps.get("verbosity")); } NMaxSample_override = ps.get("nmax", std::numeric_limits::max()); write_directory = ps.get("write_directory", ""); fIH_id = nuis::input::InputManager::Get().EnsureInputLoaded(ps); } void ProcessSample(size_t nmax = std::numeric_limits::max()) { if (fIH_id == std::numeric_limits::max()) { throw uninitialized_ISample(); } 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; ISAMPLE_INFO("Sample " << std::quoted(Name()) << " processing " << NEvsToProcess << " events."); IInputHandler::ev_index_t ev_idx = 0; size_t NSigEvents = 0; bool DetermineSignalEvents = !fSignalCache.size(); nuis::utility::Clear(fPrediction.get()); fComparisonFinalized = false; 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 (NToShout && !(ev_idx % NToShout)) { ISAMPLE_INFO("\tProcessed " << ev_idx << "/" << NEvsToProcess << " events."); } if (fSignalCache[ev_idx]) { FillProjection(fProjectionCache[ev_idx], IH.GetEventWeight(ev_idx) * nmax_scaling); } ev_idx++; } ISAMPLE_INFO("\t" << fProjectionCache.size() << "/" << NEvsToProcess << " events passed selection."); } void Write() { if (!fComparisonFinalized) { FinalizeComparison(); } nuis::persistency::WriteToOutputFile::type>( fPrediction_comparison.get(), "Prediction", write_directory); nuis::persistency::WriteToOutputFile::type>( fPrediction_xsec.get(), "Prediction_xsec", write_directory); if (fData) { nuis::persistency::WriteToOutputFile::type>( fData.get(), "Data", write_directory); nuis::persistency::WriteToOutputFile::type>( fPrediction_shape.get(), "Prediction_shape", write_directory); } } - double GetGOF() { return 1; } + double GetGOF() { + if (!fComparisonFinalized) { + FinalizeComparison(); + } + return nuis::utility::GetChi2(fData.get(), fPrediction_comparison.get()); + } double GetNDOGuess() { if (fData) { return nuis::utility::TH_traits< typename TH_dim_helper::type>::NbinsIncludeFlow(fData.get()); } return 0; } fhicl::ParameterSet GetExampleConfiguration() { fhicl::ParameterSet exps; exps.put("name", Name()); exps.put("input_type", "Generator"); exps.put("file", "Events.root"); exps.put("write_directory", Name()); exps.put("fake_data_histogram", "/path/to/file.root;h_name"); return exps; } }; typedef SimpleDataComparison<1> SimpleDataComparison_1D; typedef SimpleDataComparison<2> SimpleDataComparison_2D; #endif 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 6fa331f..fa99d28 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,158 +1,192 @@ // 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) { + 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 &ps) { if (ps.has_key("verbosity")) { SetSampleVerbosity(ps.get("verbosity")); } std::string publication = ps.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"; EnuRange = std::pair{0, 3E3}; ISAMPLE_INFO("Sample " << Name() << " specialized for publication: " << Pub_str); break; } case kPRD16: { Pub_str = "PRD16_3103"; EnuRange = std::pair{0, 6E3}; ISAMPLE_INFO("Sample " << Name() << " specialized for publication: " << Pub_str); break; } case kPRD26: { Pub_str = "PRD26_537"; EnuRange = std::pair{0, 6E3}; ISAMPLE_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()); - UseD2Corr = ps.get( - "use_D2_correction", - global_sample_configuration.get("use_D2_correction", false)); - SetData(GetDataDir() + "nuA/BubbleChamber/ANL/CCQE/ANL_CCQE_Data_" + Pub_str + ".root;ANL_1DQ2_Data"); SimpleDataComparison_1D::Initialize(ps); + UseD2Corr = ps.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 ((ISNumu.P4.E() < EnuRange.first) || (ISNumu.P4.E() > EnuRange.second)) { 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) { + nuis::utility::Fill(fPrediction_Uncorr.get(), proj, event_weight); + event_weight *= fD2CorrHist->Interpolate(proj[0]); + } + nuis::utility::Fill(fPrediction.get(), proj, event_weight); + } + + void FinalizeComparison() { + SimpleDataComparison_1D::FinalizeComparison(); + fPrediction_Uncorr->Scale(1.0, "width"); + } + + void Write() { + SimpleDataComparison_1D::Write(); + if (UseD2Corr) { + nuis::persistency::WriteToOutputFile( + fPrediction_Uncorr.get(), "Prediction_Uncorr", write_directory); + } + } + std::string Name() { return "ANL_CCQE_Evt_1DQ2_nu"; } }; DECLARE_PLUGIN(IDataComparison, ANL_CCQE_Evt_1DQ2_nu); DECLARE_PLUGIN(ISample, ANL_CCQE_Evt_1DQ2_nu); diff --git a/src/utility/CMakeLists.txt b/src/utility/CMakeLists.txt index 7b26067..abde72b 100644 --- a/src/utility/CMakeLists.txt +++ b/src/utility/CMakeLists.txt @@ -1,21 +1,25 @@ set(utility_implementation_files FileSystemUtility.cxx FullEventUtility.cxx KinematicUtility.cxx - PDGCodeUtility.cxx) + PDGCodeUtility.cxx + StringUtility.cxx) set(utility_header_files FileSystemUtility.hxx FullEventUtility.hxx InteractionChannelUtility.hxx KinematicUtility.hxx PDGCodeUtility.hxx ROOTUtility.hxx - HistogramUtility.hxx) + HistogramUtility.hxx + StringUtility.hxx + TerminalUtility.hxx + StatsUtility.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/StatsUtility.hxx b/src/utility/StatsUtility.hxx new file mode 100644 index 0000000..4e207c4 --- /dev/null +++ b/src/utility/StatsUtility.hxx @@ -0,0 +1,45 @@ +#ifndef UTILITY_STATSUTILITY_HXX_SEEN +#define UTILITY_STATSUTILITY_HXX_SEEN + +namespace nuis { +namespace utility { + +NEW_NUIS_EXCEPT(unimplemented_covariance_usage); +NEW_NUIS_EXCEPT(Mismatched_NBins); + +template +typename std::enable_if::NDims == 1, double>::type +GetChi2(HT const *a, HT const *b, TH2D const *Covariance = nullptr) { + if (Covariance) { + throw unimplemented_covariance_usage() + << "[ERROR]: Using a covariance in the Chi2 evaluation is not yet " + "implemented."; + } + + if (a->GetXaxis()->GetNbins() != b->GetXaxis()->GetNbins()) { + Mismatched_NBins() << "[ERROR]: Attempted to evaluate Chi2 between two " + "histograms with differing bin contents: NBins(a) = " + << a->GetXaxis()->GetNbins() + << ", NBins(b) = " << b->GetXaxis()->GetNbins(); + } + + double chi2 = 0; + for (Int_t bin_it_i = 0; bin_it_i < a->GetXaxis()->GetNbins(); ++bin_it_i) { + for (Int_t bin_it_j = 0; bin_it_j < a->GetXaxis()->GetNbins(); ++bin_it_j) { + double err = + 1.0 / (b->GetBinError(bin_it_i + 1) * b->GetBinError(bin_it_j + 1)); + + double contrib = + (a->GetBinContent(bin_it_i + 1) - b->GetBinContent(bin_it_i + 1)) * + (err) * + (a->GetBinContent(bin_it_j + 1) - b->GetBinContent(bin_it_j + 1)); + + chi2 += contrib; + } + } + return chi2; +} +} // namespace utility +} // namespace nuis + +#endif diff --git a/src/utility/StringUtility.cxx b/src/utility/StringUtility.cxx new file mode 100644 index 0000000..c0ceec6 --- /dev/null +++ b/src/utility/StringUtility.cxx @@ -0,0 +1,71 @@ +#include "utility/StringUtility.hxx" + +#include "string_parsers/utility.hxx" + +#include +#include +#include + +namespace nuis { +namespace utility { + +// #define DEBUG_INDENT_APPLY_WIDTH + +std::string indent_apply_width(std::string input, size_t indent, size_t width) { + std::stringstream output_stream; + std::string indent_string(indent, ' '); + + if (width) { + width--; + } + + if (width <= indent) { + return input; + } + +#ifdef DEBUG_INDENT_APPLY_WIDTH + std::cout << "[IDEBUG]: Indenting and forcing width on: " + << std::quoted(input) << ", width = " << width + << ", indent = " << indent << std::endl; +#endif + + while (input.length()) { +#ifdef DEBUG_INDENT_APPLY_WIDTH + std::cout << "[DEBUG]\tRemaining string: " << std::quoted(input) + << ", width = " << width << ", indent = " << indent << std::endl; +#endif + if (input.length() < width - indent) { +#ifdef DEBUG_INDENT_APPLY_WIDTH + std::cout << "[DEBUG]\tRemaining string is less than allowed width " + << std::quoted(input) << std::endl; +#endif + output_stream << indent_string << input << std::endl; + break; + } + size_t last_ws = input.find_last_of(" ", width - indent); + size_t first_nl = input.find_first_of("\n"); +#ifdef DEBUG_INDENT_APPLY_WIDTH + std::cout << "[DEBUG]: last_ws = " << last_ws << ", first_nl = " << first_nl + << std::endl; +#endif + size_t next_split = std::min(last_ws, first_nl); + +#ifdef DEBUG_INDENT_APPLY_WIDTH + std::cout << "[DEBUG]: Found split at " << next_split << std::endl; +#endif + if (next_split == std::string::npos) { + next_split = std::min(width - indent, input.length()); + } + std::string next_line = input.substr(0, next_split); + fhicl::string_parsers::trim(next_line); +#ifdef DEBUG_INDENT_APPLY_WIDTH + std::cout << "[DEBUG]: Streaming " << std::quoted(next_line) << std::endl; +#endif + output_stream << indent_string << next_line << std::endl; + input = input.substr(next_split + 1); + } + + return output_stream.str(); +} +} // namespace utility +} // namespace nuis diff --git a/src/config/GlobalConfiguration.hxx b/src/utility/StringUtility.hxx similarity index 73% copy from src/config/GlobalConfiguration.hxx copy to src/utility/StringUtility.hxx index c5a75dd..541393e 100644 --- a/src/config/GlobalConfiguration.hxx +++ b/src/utility/StringUtility.hxx @@ -1,37 +1,34 @@ // 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 CONFIG_GLOBALCONFIGURATION_HXX_SEEN -#define CONFIG_GLOBALCONFIGURATION_HXX_SEEN +#ifndef UTILITY_STRINGUTILITY_HXX_SEEN +#define UTILITY_STRINGUTILITY_HXX_SEEN -#include +#include "utility/TerminalUtility.hxx" -namespace fhicl { -class ParameterSet; -} +#include namespace nuis { -namespace config { -fhicl::ParameterSet const &GetDocument(std::string const &name = "global"); -bool EnsureConfigurationRead(std::string const &fhicl_file, - std::string const &name = "global"); -} // namespace config +namespace utility { +std::string indent_apply_width(std::string, + size_t indent = 0, size_t width = GetWindowWidth()); +} } // namespace nuis #endif diff --git a/src/utility/TerminalUtility.hxx b/src/utility/TerminalUtility.hxx new file mode 100644 index 0000000..3261223 --- /dev/null +++ b/src/utility/TerminalUtility.hxx @@ -0,0 +1,18 @@ +#ifndef UTILITY_TERMINALUTILITY_HXX_SEEN +#define UTILITY_TERMINALUTILITY_HXX_SEEN + +#include + +#include + +namespace nuis { +namespace utility { +inline size_t GetWindowWidth() { + struct winsize w; + ioctl(STDOUT_FILENO, TIOCGWINSZ, &w); + return w.ws_col; +} +} // namespace utility +} // namespace nuis + +#endif diff --git a/src/variation/CMakeLists.txt b/src/variation/CMakeLists.txt new file mode 100644 index 0000000..e69de29 diff --git a/src/input/InputManager.hxx b/src/variation/IVariationProvider.hxx similarity index 58% copy from src/input/InputManager.hxx copy to src/variation/IVariationProvider.hxx index 0a119f5..5f66584 100644 --- a/src/input/InputManager.hxx +++ b/src/variation/IVariationProvider.hxx @@ -1,64 +1,60 @@ // 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 CORE_INPUTMANAGER_HXX_SEEN -#define CORE_INPUTMANAGER_HXX_SEEN - -#include "input/IInputHandler.hxx" +#ifndef VARIATION_IVARIATIONPROVIDER_HXX_SEEN +#define VARIATION_IVARIATIONPROVIDER_HXX_SEEN #include "plugins/traits.hxx" #include "exception/exception.hxx" -#include -#include - namespace fhicl { class ParameterSet; } namespace nuis { -namespace input { -class InputManager { - struct NamedInputHandler { - NamedInputHandler(std::string const &, - plugins::plugin_traits::unique_ptr_t &&); - std::string name; - plugins::plugin_traits::unique_ptr_t handler; - }; - std::vector Inputs; - - InputManager(); +namespace event { +class MinimalEvent; +} // namespace event +} // namespace nuis - static InputManager *_global_inst; +class IVariationProvider { public: + typedef int paramId_t; + static paramId_t const kParamUnhandled = + std::numeric_limits::max(); - static InputManager &Get(); + virtual void Initialize(fhicl::ParameterSet const &) = 0; - NEW_NUIS_EXCEPT(unknown_input); - typedef size_t Input_id_t; + paramId_t GetParameterId(std::string const &) = 0; - Input_id_t EnsureInputLoaded(fhicl::ParameterSet const &); - Input_id_t GetInputId(std::string const &) const; - IInputHandler const &GetInputHandler(Input_id_t) const; + bool HandlesParameter(std::string const ¶m_name) { + return (GetParameterId(param_name) != kParamUnhandled); + } + + void SetParameterValue() = 0; + bool ParametersVaried() = 0; + void Reconfigure() = 0; + + virtual ~IVariationProvider() {} }; -} // namespace input -} // namespace nuis + +DECLARE_PLUGIN_INTERFACE(IVariationProvider); #endif diff --git a/src/config/GlobalConfiguration.hxx b/src/variation/IWeightProvider.hxx similarity index 69% copy from src/config/GlobalConfiguration.hxx copy to src/variation/IWeightProvider.hxx index c5a75dd..dfaf751 100644 --- a/src/config/GlobalConfiguration.hxx +++ b/src/variation/IWeightProvider.hxx @@ -1,37 +1,34 @@ // 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 CONFIG_GLOBALCONFIGURATION_HXX_SEEN -#define CONFIG_GLOBALCONFIGURATION_HXX_SEEN +#ifndef VARIATION_IWEIGHTPROVIDER_HXX_SEEN +#define VARIATION_IWEIGHTPROVIDER_HXX_SEEN -#include +#include "variation/IVariationProvider.hxx" -namespace fhicl { -class ParameterSet; -} +class IWeightProvider { +public: + double GetEventWeight(nuis::event::MinimalEvent) = 0; -namespace nuis { -namespace config { -fhicl::ParameterSet const &GetDocument(std::string const &name = "global"); -bool EnsureConfigurationRead(std::string const &fhicl_file, - std::string const &name = "global"); -} // namespace config -} // namespace nuis + virtual ~IWeightProvider() {} +}; + +DECLARE_PLUGIN_INTERFACE(IWeightProvider); #endif diff --git a/src/variation/VariationManager.cxx b/src/variation/VariationManager.cxx new file mode 100644 index 0000000..e69de29 diff --git a/src/input/InputManager.hxx b/src/variation/VariationManager.hxx similarity index 55% copy from src/input/InputManager.hxx copy to src/variation/VariationManager.hxx index 0a119f5..773d1c6 100644 --- a/src/input/InputManager.hxx +++ b/src/variation/VariationManager.hxx @@ -1,64 +1,76 @@ // 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 CORE_INPUTMANAGER_HXX_SEEN -#define CORE_INPUTMANAGER_HXX_SEEN +#ifndef VARIATION_VARIATIONMANAGER_HXX_SEEN +#define VARIATION_VARIATIONMANAGER_HXX_SEEN -#include "input/IInputHandler.hxx" +#include "variation/IWeightProvider.hxx" #include "plugins/traits.hxx" #include "exception/exception.hxx" #include #include namespace fhicl { class ParameterSet; } namespace nuis { namespace input { -class InputManager { - struct NamedInputHandler { - NamedInputHandler(std::string const &, - plugins::plugin_traits::unique_ptr_t &&); +class VariationManager { + struct NamedVariationProvider { + NamedVariationProvider( + std::string const &, + plugins::plugin_traits::unique_ptr_t &&); std::string name; - plugins::plugin_traits::unique_ptr_t handler; + plugins::plugin_traits::unique_ptr_t handler; }; - std::vector Inputs; + std::vector VarProvs; - InputManager(); + VariationManager(); + + static VariationManager *_global_inst; - static InputManager *_global_inst; public: + typedef size_t VarProv_id_t; + typedef size_t paramId_t; + + struct VariationProviderParameter { + std::string name; + VarProv_id_t providerId; + IVariationProvider::paramId_t providerParameterId; + }; + + std::vector VarParams; + - static InputManager &Get(); + static VariationManager &Get(); - NEW_NUIS_EXCEPT(unknown_input); - typedef size_t Input_id_t; + paramId_t EnsureParameterHandled(fhicl::ParameterSet const &); + void SetParameterValue(paramId_t, double); + void GetParameterPull(paramId_t); - Input_id_t EnsureInputLoaded(fhicl::ParameterSet const &); - Input_id_t GetInputId(std::string const &) const; - IInputHandler const &GetInputHandler(Input_id_t) const; + double GetEventWeight(); }; } // namespace input } // namespace nuis #endif