SpectrumTest_module.cc
Go to the documentation of this file.
1 /*
2  * \file SpectrumTester_module.cc
3  * \brief Module for testing that the Spectrum class works as designed.
4  *
5  * Created on: \date Mar 21, 2016
6  * Original author: \author J. Wolcott <jwolcott@fnal.gov>
7  *
8  */
9 
10 #include <memory>
11 #include <unordered_map>
12 
13 #include "TFile.h"
14 #include "TH1F.h"
15 
22 
23 #include "FNEX/core/Spectrum.h"
25 #include "FNEX/core/Cuts.h"
26 #include "FNEX/core/Shifters.h"
27 #include "FNEX/core/Weighters.h"
28 #include "FNEX/vars/CommonVars.h"
29 
30 namespace fnex
31 {
32 
34  {
35  public:
36  //======= Class constructor and destructor
37  explicit SpectrumTest(fhicl::ParameterSet const& pset);
38  virtual ~SpectrumTest() = default;
39 
40  //======= ART methods
41  void beginJob() override;
42  void beginRun(const art::Run& run) override {};
43  void beginSubRun(const art::SubRun& sr) override {};
44  void analyze(const art::Event& evt) override {};
45  void reconfigure(fhicl::ParameterSet const& p) override;
46  void endSubRun(const art::SubRun& sr) override {};
47  void endRun(const art::Run& run) override {};
48  void endJob() override;
49 
50  private:
51  void CreateSpectra();
52  void InitializeWeighters(fhicl::ParameterSet const& pset);
53  void InitializeShifters(const fhicl::ParameterSet& pset);
54  void InspectEventLists() const;
55  void SavePlots() const;
56 
57  std::unique_ptr<fnex::EventListManipulator> fManipulator; ///< helper object to work with event lists from TTrees
58  fnex::EventListMap fEventLists; ///< event lists
59  std::unordered_map<std::string, std::shared_ptr<fnex::Weighter>> fWeighters;
60  std::unordered_map<std::string, std::shared_ptr<fnex::Shifter>> fShifters;
61  std::unordered_map<fnex::MetaData, fnex::Spectrum, fnex::MetaDataHasher> fSpectra;
62 
63  // fhicl properties
64  std::vector<std::string> fVarNames;
65  };
66 
67  // ====================
68 
70  : art::EDAnalyzer(pset)
71  {
72  // do the necessary setup
73  this->reconfigure(pset);
74  this->InitializeWeighters(pset);
75  this->InitializeShifters(pset);
76  }
77 
78  // ====================
79 
81  {
82  this->fManipulator = std::make_unique<fnex::EventListManipulator>(pset);
83 
84  fVarNames = pset.get< std::vector<std::string> >("SpectrumVars");
85  }
86 
87  // ====================
88 
90  {
91  this->fEventLists = this->fManipulator->Deserialize(); // copy the event lists from the deserializer.
92 
93  mf::LogDebug("SpectrumTest") << "Got " << this->fEventLists.size() << " event lists from the deserializer.";
94 
95  this->CreateSpectra();
96  }
97 
98  // ====================
99 
101  {
102  this->SavePlots();
103  this->InspectEventLists();
104  } // SpectrumTest::endJob()
105 
106  // ====================
107 
109  {
110  for (auto & evListPair : this->fEventLists)
111  {
112  this->fSpectra.emplace(std::piecewise_construct,
113  std::make_tuple(evListPair.first),
114  std::make_tuple(evListPair.second)
115  );
116  }
117  } // SpectrumTest::CreateSpectra()
118 
119  // ====================
120 
122  {
123  auto oscParams_config = pset.get<fhicl::ParameterSet>("OscillationParams");
124  std::unordered_map<std::string, double> oscParams;
125  // unfortunately you can't get a map back from a fhicl::ParameterSet, so we're forced to make do with copying
126  for (const auto & pName : oscParams_config.get_all_keys())
127  {
128 // std::cout << " loading oscillation parameter with name: " << pName << std::endl;
129  oscParams[pName] = oscParams_config.get<double>(pName);
130  }
131  this->fWeighters.emplace("OscillationWgt", std::make_shared<fnex::Weighter>("OscillationWgt", new fnex::OscillationWeightFn(oscParams)));
132  }
133 
135  {
136  this->fShifters["EnergyShift"] = fnex::kShifterRegistry.at("RecoHadEShifter");
137  }
138 
139  // ====================
140 
142  {
144  const fnex::Spectrum & spec = this->fSpectra.at(md);
145 
146 // std::cout << "Stored histograms for spectrum " << md.ToString() << ":" << std::endl;
147 // for (const auto & hsetPair : spec.Histograms())
148 // {
149 // std::cout << " Variable '" << hsetPair.first << "':" << std::endl;
150 //
151 // for (const auto & histObj : hsetPair.second)
152 // std::cout << " hash = " << histObj.histhash.hash.to_string() << ", name = " << histObj.hist->GetName() << ", ptr = " << histObj.hist.get() << std::endl;
153 // } // for (hsetPair)
154 
155 
156  const auto & ev = spec.Events().EventNum(35);
157  std::cout << "Event tree for event #35:" << std::endl;
158  std::cout << " Weights:" << std::endl;
159  for (const auto & wgtPair : ev->weights)
160  std::cout << " name = " << wgtPair.first->Name() << "; wgt = " << wgtPair.second.first << std::endl;
161 
162  std::cout << "\n CV vals:" << std::endl;
163  for (const auto & varPair : ev->cvVals)
164  std::cout << " name = " << varPair.first << ", ptr = " << varPair.second.get() << ", val = " << std::dynamic_pointer_cast<fnex::SimpleVar>(varPair.second)->val() << std::endl;
165 
166  std::cout << "\n shifted vals:" << std::endl;
167  fnex::ShiftCollectionHasher hasher;
168  for (const auto & shiftPair : ev->shiftedVals)
169  {
170  std::cout << " shift hash: " << hasher(shiftPair.first) << std::endl;
171  for (const auto & varPair : shiftPair.second)
172  std::cout << " name = " << varPair.first << ", ptr = " << varPair.second.get() << ", val = " << std::dynamic_pointer_cast<fnex::SimpleVar>(varPair.second)->val() << std::endl;
173  }
174  } // SpectrumTest::InspectEventLists()
175 
176  // ====================
177 
179  {
180 // std::cout << "Weight parameter names:\n";
181 // for (const auto & paramPair : this->fWeighters.at("OscillationWgt").GetDefaultConfig())
182 // std::cout << " " << paramPair.first << " = " << paramPair.second << "\n";
183 // std::cout << std::endl;
184 
185  std::string varName = fnex::var::Nu_RecoE::Name();
186 
187  std::unique_ptr<TFile> f(TFile::Open("/nova/ana/users/jwolcott/scratch/spectrum_test_hists.root", "recreate"));
188  mf::LogDebug("SpectrumTest") << "Got spectra for metadata:";
189  for (auto & specPair : this->fSpectra)
190  {
191  mf::LogDebug("SpectrumTest") << " " << specPair.first.ToString();
192 
193  const TH1F * hSpec = dynamic_cast<const TH1F*>(specPair.second.Histogram(varName));
194 
195  TH1F * hnew = new TH1F(*hSpec);
196  hnew->SetName( (varName + "_" + specPair.first.ToString()).c_str() );
197  hnew->Write();
198 
199  if (!specPair.first.isMC)
200  continue;
201 
202  const std::string wgtName = "OscillationWgt";
203 // std::string wgtName = "TestWgt";
204 
205  const TH1F * hSpecWgtd = dynamic_cast<const TH1F*>(specPair.second.Histogram(varName, {{}}, {{this->fWeighters.at(wgtName), this->fWeighters.at(wgtName)->State()}}));
206  TH1F * hnewWgtd = new TH1F(*hSpecWgtd);
207  hnewWgtd->SetName( (varName + "_" + specPair.first.ToString() + "_" + this->fWeighters.at(wgtName)->Name()).c_str() );
208  hnewWgtd->Write();
209 
210  const std::string shiftName = "EnergyShift";
211  int sigma, sigmaMax;
212  if (specPair.first.ToString() == "MCNearNonSwapNuMuSel")
213  sigma = -5, sigmaMax = 6;
214  else
215  sigma = 0; sigmaMax = 1;
216  for (; sigma < sigmaMax; sigma++)
217  {
218  const TH1F * hSpecShifted = dynamic_cast<const TH1F*>(specPair.second.Histogram(varName, fnex::ShiftCollection{{this->fShifters.at(shiftName), sigma}}));
219  TH1F * hnewShifted = new TH1F(*hSpecShifted);
220  hnewShifted->SetName( Form("%s_%s_%s_shift=%dsigma", varName.c_str(), specPair.first.ToString().c_str(), this->fShifters.at(shiftName)->Name().c_str(), sigma) );
221  hnewShifted->Write();
222  }
223 
224  const TH1F * hSpecNC = dynamic_cast<const TH1F*>(specPair.second.Histogram(varName, {{}}, {{}}, fnex::CutCollection{std::make_shared<fnex::Cut>(fnex::kIsNC)}));
225  // might be NULL if no events passed...
226  if (hSpecNC)
227  {
228  TH1F * hnewNC = new TH1F(*hSpecNC);
229  hnewNC->SetName( Form("%s_%s_cut=NC", varName.c_str(), specPair.first.ToString().c_str()) );
230  hnewNC->Write();
231  }
232  } // for (specPair)
233  f->Close();
234 
235  } // SpectrumTest::SavePlots()
236 
238 
239 } // namespace fnex
240 
void endSubRun(const art::SubRun &sr) override
void analyze(const art::Event &evt) override
const char * p
Definition: xmltok.h:285
void InitializeShifters(const fhicl::ParameterSet &pset)
std::map< fnex::MetaData, fnex::EventList > EventListMap
Definition: Event.h:186
SpectrumTest(fhicl::ParameterSet const &pset)
Create a list of fnex::Events to be used in fits.
DEFINE_ART_MODULE(SpectrumTest)
virtual ~SpectrumTest()=default
Definition: Run.h:31
std::vector< std::string > fVarNames
const Cut kIsNC([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return!sr->mc.nu[0].iscc;})
Is this a Neutral Current event?
Definition: TruthCuts.h:8
The Spectrum is the means for interaction with Event collections.
Definition: Spectrum.h:29
void beginSubRun(const art::SubRun &sr) override
void InspectEventLists() const
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
Near Detector in the NuMI cavern.
caf::StandardRecord * sr
void endRun(const art::Run &run) override
void reconfigure(fhicl::ParameterSet const &p) override
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
double sigma(TH1F *hist, double percentile)
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
std::unordered_map< fnex::MetaData, fnex::Spectrum, fnex::MetaDataHasher > fSpectra
std::unique_ptr< fnex::EventListManipulator > fManipulator
helper object to work with event lists from TTrees
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
void InitializeWeighters(fhicl::ParameterSet const &pset)
Service to store calibration data products (CDP) in the SQLite3 metadatabase of a file...
Definition: FillParentInfo.h:8
fnex::EventList const & Events() const
Get the EventList that was supplied to this Spectrum.
Definition: Spectrum.h:57
std::unordered_map< std::string, std::shared_ptr< fnex::Shifter > > fShifters
fnex::EventListMap fEventLists
event lists
std::unordered_map< std::string, std::shared_ptr< fnex::Weighter > > fWeighters
void beginJob() override
EventPtr EventNum(std::size_t idx)
Definition: Event.h:113
void beginRun(const art::Run &run) override
enum BeamMode string