SystematicsWeightTest_module.cc
Go to the documentation of this file.
1 /*
2  * \file SystematicsWeightTest_module.cc
3  * \brief Module for testing that the interpolation used
4  * in SystematicsWeightFn type fnex::Weighters
5  * works correctly.
6  *
7  * Created on: \date May 23, 2016
8  * Original author: \author J. Wolcott <jwolcott@fnal.gov>
9  *
10  */
11 
12 #include "TFile.h"
13 
20 
21 #include "FNEX/core/Weighter.h"
22 #include "FNEX/core/Spectrum.h"
23 #include "FNEX/core/Weighters.h"
25 
26 namespace fnex
27 {
28  /// Stupid test weight that produces a linear shift for every event, where 1 sigma = 50%.
29  // class StupidWeight_Linear : public fnex::SystematicsWeightFn
30  class StupidWeight_Linear : public fnex::InterpSystWgtFn
31  {
32  private:
33  //std::map<double, double> CalcWeights(const fnex::Event&) const override { return { {-1.0, 0.5}, {1.0, 1.5} }; };
34  std::map<float, float> CalcWeights(const fnex::Event&) const override { return { {-1.0, 0.5}, {1.0, 1.5} }; };
35  };
36 
37 
39  {
40  public:
41  //======= Class constructor and destructor
42  explicit SystematicsWeightTest(fhicl::ParameterSet const& pset);
43  virtual ~SystematicsWeightTest() = default;
44 
45  //======= ART methods
46  void beginJob() override;
47  void beginRun(const art::Run& run) override {};
48  void beginSubRun(const art::SubRun& sr) override {};
49  void analyze(const art::Event& evt) override {};
50  void reconfigure(fhicl::ParameterSet const& p) override;
51  void endSubRun(const art::SubRun& sr) override {};
52  void endRun(const art::Run& run) override {};
53  void endJob() override;
54 
55  private:
56  void CreateSpectra();
57  void InitializeWeighters(fhicl::ParameterSet const& pset);
58  void MakePlots();
59  void SavePlots();
60 
61 
62  std::unique_ptr<fnex::EventListManipulator> fManipulator; ///< helper object to work with event lists from TTrees
63  fnex::EventListMap fEventLists; ///< event lists
64  std::unordered_map<std::string, std::shared_ptr<fnex::Weighter>> fWeighters;
65  std::unordered_map<fnex::MetaData, fnex::Spectrum, fnex::MetaDataHasher> fSpectra;
66 
67  // fhicl properties
68  std::vector<std::string> fVarNames;
69  };
70 
71  // ====================
72 
74  : art::EDAnalyzer(pset)
75  {
76  // do the necessary setup
77  this->reconfigure(pset);
78  this->InitializeWeighters(pset);
79  }
80 
81  // ====================
83  {
84  this->fEventLists = this->fManipulator->Deserialize(); // copy the event lists from the deserializer.
85 
86  mf::LogDebug("SystematicsWeightTest") << "Got " << this->fEventLists.size() << " event lists from the deserializer.";
87 
88  this->CreateSpectra();
89  }
90 
91  // ====================
93  {
94  this->SavePlots();
95  }
96 
97  // ====================
99  {
100  this->fManipulator = std::make_unique<fnex::EventListManipulator>(pset);
101 
102  fVarNames = pset.get< std::vector<std::string> >("SpectrumVars");
103  }
104 
105  // ====================
107  {
108  for (auto & evListPair : this->fEventLists)
109  {
110  this->fSpectra.emplace(std::piecewise_construct,
111  std::make_tuple(evListPair.first),
112  std::make_tuple(evListPair.second)
113  );
114  }
115  } // SystematicsWeightTest::CreateSpectra()
116 
117  // ====================
118 
120  {
121  StupidWeight_Linear * s = new StupidWeight_Linear; // this is leaked. in a real module you'd want to keep track of this guy
122  this->fWeighters.emplace("StupidLinearWgt", std::make_shared<fnex::Weighter>("StupidLinearWgt", s));
123  std::cout << " weight function address:" << this->fWeighters.at("StupidLinearWgt")->WeightFunctionObj() << std::endl;
124  }
125 
126  // ====================
127 
129  {
130  std::string varName = fnex::var::Nu_RecoE::Name();
131 
132  std::unique_ptr<TFile> f(TFile::Open("/nova/ana/users/jwolcott/scratch/weight_test_hists.root", "recreate"));
133  mf::LogDebug("SystematicsWeightTest") << "Got spectra for metadata:";
134  for (auto & specPair : this->fSpectra)
135  {
136  if (!specPair.first.isMC)
137  continue;
138 
139  mf::LogDebug("SystematicsWeightTest") << " " << specPair.first.ToString();
140 
141  const TH1F * hSpec = dynamic_cast<const TH1F*>(specPair.second.Histogram(varName));
142  TH1F * hnew = new TH1F(*hSpec);
143  hnew->SetName( (varName + "_" + specPair.first.ToString()).c_str() );
144  hnew->Write();
145 
146  for (auto & wgtPair : this->fWeighters)
147  {
148  const auto & wgtName = wgtPair.first;
149  auto & wgtr = wgtPair.second;
150  for (const auto sigma : {-0.5, 0.0, 1.0, 1.5})
151  {
152  //auto wgtFnPtr = dynamic_cast<fnex::SystematicsWeightFn*>(wgtr->WeightFunctionObj());
153  auto wgtFnPtr = dynamic_cast<fnex::InterpSystWgtFn*>(wgtr->WeightFunctionObj());
154  if (wgtFnPtr)
155  wgtFnPtr->SetSigma(sigma);
156  else
157  {
158  std::cout << " wgt function pointer address: " << wgtr->WeightFunctionObj() << std::endl;
159  LOG_ERROR("SystematicsWeightTest") << "Weight function pointer is of type: " << typeid(wgtr->WeightFunctionObj()).name();
160  throw cet::exception("SystematicsWeightTest::SavePlots(): couldn't set sigma for weight function pointer!");
161  }
162 
163  std::cout << " calculating weights for sigma = " << sigma << std::endl;
164  std::cout << " (weighter state = " << wgtr->State() << ")" << std::endl;
165  const TH1F * hSpecWgtd = dynamic_cast<const TH1F*>(specPair.second.Histogram(varName, {{}}, {{wgtr, wgtr->State()}}));
166  TH1F * hnewWgtd = new TH1F(*hSpecWgtd);
167  hnewWgtd->SetName( Form("%s_%s_%s_sigma=%.2f", varName.c_str(), specPair.first.ToString().c_str(), wgtName.c_str(), sigma) );
168  hnewWgtd->Write();
169  }
170  }
171  }
172  }
173 
175 } // namespace fnex
const XML_Char * name
Definition: expat.h:151
void MakePlots(std::string fname)
void beginRun(const art::Run &run) override
void endRun(const art::Run &run) override
const char * p
Definition: xmltok.h:285
void analyze(const art::Event &evt) override
void reconfigure(fhicl::ParameterSet const &p) override
std::unordered_map< fnex::MetaData, fnex::Spectrum, fnex::MetaDataHasher > fSpectra
void InitializeWeighters(fhicl::ParameterSet const &pset)
void beginSubRun(const art::SubRun &sr) override
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::map< fnex::MetaData, fnex::EventList > EventListMap
Definition: Event.h:186
std::vector< std::string > fVarNames
Create a list of fnex::Events to be used in fits.
Stupid test weight that produces a linear shift for every event, where 1 sigma = 50%.
DEFINE_ART_MODULE(SpectrumTest)
Definition: Run.h:31
const XML_Char * s
Definition: expat.h:262
std::map< float, float > CalcWeights(const fnex::Event &) const override
std::unordered_map< std::string, std::shared_ptr< fnex::Weighter > > fWeighters
void beginJob()
T get(std::string const &key) const
Definition: ParameterSet.h:231
void endSubRun(const art::SubRun &sr) override
int evt
caf::StandardRecord * sr
double sigma(TH1F *hist, double percentile)
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
fnex::EventListMap fEventLists
event lists
Service to store calibration data products (CDP) in the SQLite3 metadatabase of a file...
Definition: FillParentInfo.h:8
SystematicsWeightTest(fhicl::ParameterSet const &pset)
#define LOG_ERROR(stream)
Definition: Messenger.h:129
std::unique_ptr< fnex::EventListManipulator > fManipulator
helper object to work with event lists from TTrees
enum BeamMode string