CMFRandomUniverses_plugin.cc
Go to the documentation of this file.
1 //
2 // Created by Brian Rebel on 10/29/19.
3 //
4 
5 // Module to create fake data for multiple random universes
6 //
7 
8 #include <memory>
9 #include <string>
10 #include <vector>
11 #include <map>
12 #include <set>
13 
14 // Framework includes
21 #include "fhiclcpp/ParameterSet.h"
23 
24 // NOvA includes
36 
37 #include "TFile.h"
38 #include "TH1.h"
39 #include "TH2.h"
40 #include "TRandom3.h"
41 
42 namespace cmf {
43 
45  public:
46  explicit RandomUniverses(fhicl::ParameterSet const& pset);
47  virtual ~RandomUniverses();
48 
49  // Plugins should not be copied or assigned.
50  RandomUniverses(RandomUniverses const &) = delete;
51  RandomUniverses(RandomUniverses &&) = delete;
52  RandomUniverses & operator = (RandomUniverses const &) = delete;
54 
55  void readResults (art::Results const& r) override;
56  void writeResults(art::Results & r) override;
57  void clear() override;
58  void beginJob() override;
59  void reconfigure(fhicl::ParameterSet const& p);
60 
61  private:
62 
63  std::map<cmf::DetType_t, fhicl::ParameterSet> fManipulatorPars; ///< parameters to make object for deserializing events
64  fhicl::ParameterSet fSAWPars; ///< parameters for the shifter and weighter
65  bool fIsFCJob; ///< are we making Feldman Cousins points or not?
66  unsigned int fStartUniverse; ///< which universe to start?
67  TH2D* fPoissonDraws; ///< histogram recording the Poisson draws for each logical bin
68  std::vector<cmf::Location> fUniverseLocs; ///< vector of input points for the random universes
69  std::vector<cmf::Location> fAsimovLocs; ///< the Asimov input points, ie no systematic shifts
70  };
71 
72  //......................................................................
74  : fSAWPars (pset.get<fhicl::ParameterSet>("ShifterAndWeighterParameters"))
75  , fPoissonDraws(nullptr)
76  {
77  this->reconfigure(pset);
78 
79  LOG_DEBUG("RandomUniverses")
80  << "ShifterAndWeighter parameters: "
81  << fSAWPars.to_string();
82 
83  // produce a vector of FakeUniverses
84  produces< std::vector<cmf::FakeUniverse> >();
85  }
86 
87  //......................................................................
89  = default;
90 
91  //......................................................................
93  {
94  cmf::SelectionUtility::Instance() ->Initialize(p.get< std::vector<fhicl::ParameterSet> >("SelectionsToUse", std::vector<fhicl::ParameterSet>() ));
97 
98  fManipulatorPars.emplace(cmf::kFARDET, p.get<fhicl::ParameterSet>("FarDetEventListManipulator" ));
99  fManipulatorPars.emplace(cmf::kNEARDET, p.get<fhicl::ParameterSet>("NearDetEventListManipulator"));
100  fIsFCJob = p.get<bool >("IsFeldmanCousinsJob", false);
101  fStartUniverse = p.get<unsigned int>("StartUniverse");
102 
103  // set the seed for the random number to be what is configured, plus the
104  // start universe value in case lots are jobs are submitted with the same seed
105  // by accident
107  }
108 
109  //......................................................................
111  {
113 
114  LOG_DEBUG("RandomUniverse")
115  << "nominal location is\n"
116  << nominalLoc;
117 
118  // get a list of the oscillation parameter names
119  std::set<std::string> oscParNames;
120  for(auto const& itr : nominalLoc.FDLoc){
121  if(itr.second.IsOscPar())
122  oscParNames.emplace(itr.first);
123  }
124 
125  if(fIsFCJob){
126  // use the RandomUniverseUtility to get our random values for the hidden
127  // parameters in oscillation space first then use those oscillation
128  // parameters to create systematically varied points as well
130  fAsimovLocs);
131 
132  // now loop over the Asimov points and for each one create a template
133  // to be used to get the systematic uncertainty parameter variations
134  fUniverseLocs.resize(fAsimovLocs.size(), nominalLoc);
135 
136  for(size_t u = 0; u < fAsimovLocs.size(); ++u){
137  // grab the oscillation parameters from this Asimov point
138  // the oscillation parameter values are the same in both detectors
139  for(auto const& itr : oscParNames){
140  fUniverseLocs[u].NDLoc.find(itr)->second.SetValue(fAsimovLocs[u].NDLoc.find(itr)->second.Value());
141 
142  fUniverseLocs[u].FDLoc.find(itr)->second.SetValue(fAsimovLocs[u].FDLoc.find(itr)->second.Value());
143  }
144 
146  fUniverseLocs[u],
147  u);
148  }
149  }
150  else{
151  // In this mode we are only going to change the systematic parameters but not the
152  // oscillation parameters so we only need 1 Asimov point
153  fAsimovLocs.emplace_back(nominalLoc);
154 
155  // get the set of locations for our universes
157  fUniverseLocs);
158  }
159 
160  for(size_t u = 0; u < fUniverseLocs.size(); ++u)
161  LOG_VERBATIM("RandomUniverses")
162  << " new location "
163  << u
164  << "\n"
165  << fUniverseLocs[u];
166 
167  // now make some histograms
169 
170  fPoissonDraws = tfs->make<TH2D>("PoissonDraws",
171  ";Energy Bin;Poisson Draw",
173  0,
175  fUniverseLocs.size(),
176  0.,
177  fUniverseLocs.size());
178 
179  }
180 
181  //......................................................................
182  // Method to clear out the collections of data products after the
183  // writeResults method is called.
185  {
186  }
187 
188  //......................................................................
190  {
191  }
192 
193 
194  //......................................................................
196  {
197  cmf::EventListColl mcLists;
198 
199  std::vector<cmf::Spectrum> spectrum (fUniverseLocs.size());
200  std::vector<cmf::Spectrum> poissonSpectrum(fUniverseLocs.size());
201  std::vector<cmf::Spectrum> fakeCount (fUniverseLocs.size());
202  std::vector<cmf::Spectrum> count (fUniverseLocs.size());
203  std::vector<cmf::Spectrum> asimov (fAsimovLocs.size() );
204 
205  std::set<cmf::DetType_t> dets;
206  for(auto const& mpItr : fManipulatorPars){
207 
208  mcLists.clear();
209 
210  dets.clear();
211  dets.insert(mpItr.first);
212 
213  LOG_VERBATIM("RandomUniverses")
214  << "Fill spectrum from "
215  << cmf::cDetType_Strings[mpItr.first]
216  << " starting universe "
217  << fStartUniverse
218  << " "
219  << fUniverseLocs.size();
220 
221  // get the manipulator to deserialize the data and MC
222  // we'll keep the MC lists around, but will drop the data lists as
223  // soon as we leave this method
224  cmf::EventListManipulator elm(mpItr.second);
225 
226  elm.Deserialize(mcLists, cmf::kMC, dets);
227  auto const& exposureMap = elm.ExposureMap();
228 
229  // now make the spectra for all the universes
230  unsigned int uni = 0;
231  for(unsigned int u = fStartUniverse; u < fStartUniverse + fUniverseLocs.size(); ++u){
232 
233  uni = u - fStartUniverse;
234 
235  if(fIsFCJob){
236  LOG_DEBUG("RandomUniverses")
237  << "Filling FC Asimov Spectra";
238 
240  fSAWPars);
241 
242  // fill the asimov spectrum
243  cmf::FillSpectrum(mcLists,
244  exposureMap,
245  asimov[uni],
246  count[uni]);
247  }
248  else if(!fIsFCJob && u == 0){
249  LOG_DEBUG("RandomUniverses")
250  << fIsFCJob
251  << " "
252  << u
253  << " fill Asimov spectrum for "
254  << fAsimovLocs.size()
255  << " locations";
256 
258  fSAWPars);
259  // fill the asimov spectrum
260  cmf::FillSpectrum(mcLists,
261  exposureMap,
262  asimov[u],
263  count[u]);
264 
265  LOG_DEBUG("RandomUniverses")
266  << "Asimov spectrum filled "
267  << asimov.size()
268  << " "
269  << count.size();
270 
271  }
272 
273  LOG_DEBUG("CMFRandomUniverses")
274  << "Finished filling asimov spectrum";
275 
277  fSAWPars);
278 
279  LOG_DEBUG("CMFRandomUniverses")
280  << "Filling universe spectrum";
281 
282  cmf::FillSpectrum(mcLists,
283  exposureMap,
284  spectrum[uni],
285  fakeCount[uni]);
286 
287  cmf::RandomUniverseUtility::Instance()->FillPoissonSpectrum(spectrum[uni], poissonSpectrum[uni]);
288 
289  for(size_t b = 0; b < poissonSpectrum[uni].size(); ++b){
290  // the ND bins fill twice without the following check because
291  // the ND portion of the spectrum is filled before the FD portion
292  if(mcLists.begin()->first.detector == cmf::CovarianceBinUtility::Instance()->BinToDetector(b))
293  fPoissonDraws->Fill(b, uni, poissonSpectrum[uni][b]);
294  }
295 
296  // for(size_t b = 0; b < spectrum[uni].size(); ++b){
297  // LOG_VERBATIM("RandomUniverses")
298  // << "universe "
299  // << u
300  // << " spectrum bin "
301  // << b
302  // << " fluctuated spectrum: "
303  // << poissonSpectrum[uni][b]
304  // << " unfluctuated/asimov: "
305  // << spectrum[uni][b]
306  // << " / "
307  // << asimov[0][b];
308  // }
309 
310  LOG_VERBATIM("RandomUniverses")
311  << " Finished universe "
312  << uni
313  << " for "
314  << cmf::cDetType_Strings[mpItr.first];
315 
316  } // end loop to make fake universe spectra
317  } // end loop over the manipulators
318 
319  // actually make the fake universes
320  std::unique_ptr<std::vector<cmf::FakeUniverse>> fakeUniverses = std::make_unique<std::vector<cmf::FakeUniverse>>();
321 
322  for(size_t u = 0; u < fUniverseLocs.size(); ++u){
323  cmf::Spectrum asimovSpec(asimov[0]);
324  if(fIsFCJob) asimovSpec.swap(asimov[u]);
325  fakeUniverses->emplace_back(fUniverseLocs[u],
326  asimovSpec,
327  poissonSpectrum[u],
328  fakeCount[u]);
329 
330  LOG_DEBUG("RandomUniverses")
331  << "added fake universe "
332  << u;
333 
334  // make histograms for the fake spectra and the total spectrum to compare
335  // them against each other.
337  u + fStartUniverse);
338 
339  LOG_DEBUG("RandomUniverses")
340  << fakeUniverses->back();
341  }
342 
343  r.put(std::move(fakeUniverses));
344  } // end writeResults
345 
347 
348 } // end cmf namespace
349 
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
void writeResults(art::Results &r) override
void Initialize(fhicl::ParameterSet const &pset)
void Initialize(fhicl::ParameterSet const &pset)
std::vector< cmf::Location > fAsimovLocs
the Asimov input points, ie no systematic shifts
const char * p
Definition: xmltok.h:285
void FillPoissonSpectrum(std::vector< float > &spectrum, std::vector< float > &poissonSpectrum)
static SelectionUtility * Instance()
ParameterSpaceLoc FDLoc
Definition: Parameter.h:31
art::ProductID put(std::unique_ptr< PROD > &&product)
Definition: Results.h:78
virtual ~RandomUniverses()
std::vector< cmf::Location > fUniverseLocs
vector of input points for the random universes
TH2D * fPoissonDraws
histogram recording the Poisson draws for each logical bin
fhicl::ParameterSet fSAWPars
parameters for the shifter and weighter
static ParameterUtility * Instance()
void CreateLocationWithVariedSystematics(cmf::Location const &nominalLoc, cmf::Location &vecLoc, size_t universeNum)
unsigned int fStartUniverse
which universe to start?
RandomUniverses(fhicl::ParameterSet const &pset)
void FillSpectrum(cmf::EventListColl const &eventLists, cmf::ExposureMap const &exposureMap, std::vector< float > &spectrum, std::vector< float > &count, cmf::InteractionType_t intType, bool applyOscWeight, bool applyExposureNorm)
RandomUniverses & operator=(RandomUniverses const &)=delete
void Deserialize(cmf::EventListColl &eventLists, cmf::DataMC_t dataMC=cmf::kBoth, std::set< cmf::DetType_t > const &detectors=std::set< cmf::DetType_t >({cmf::kNEARDET, cmf::kFARDET}))
void reconfigure(fhicl::ParameterSet const &p)
#define DEFINE_ART_RESULTS_PLUGIN(klass)
T get(std::string const &key) const
Definition: ParameterSet.h:231
cmf::ExposureMap const & ExposureMap() const
void Initialize(fhicl::ParameterSet const &pset)
static ShifterAndWeighter * Instance()
std::vector< float > Spectrum
Definition: Constants.h:527
void MakeUniverseHistograms(cmf::FakeUniverse &fakeU, unsigned int uNumber)
size_t TotalBins(bool allSels=false)
Module to combine a set of results into a single file currently only does one data product type at a ...
Definition: Event.cxx:24
bool fIsFCJob
are we making Feldman Cousins points or not?
void CreateLocationsWithVariedSystematics(cmf::Location const &nominalLoc, std::vector< cmf::Location > &vecLocs)
void InitShiftsAndWeightsToUse(cmf::Location const &loc, fhicl::ParameterSet const &parset)
static RandomUniverseUtility * Instance()
T * make(ARGS...args) const
void Initialize(std::vector< fhicl::ParameterSet > const &pset)
cmf::DetType_t BinToDetector(int const &bin, bool allSels=false)
const hit & b
Definition: hits.cxx:21
cmf::Location ParametersAsLocation()
void readResults(art::Results const &r) override
TRandom3 r(0)
std::map< cmf::DetType_t, fhicl::ParameterSet > fManipulatorPars
parameters to make object for deserializing events
#define LOG_VERBATIM(category)
std::map< cmf::MetaData, cmf::EventList > EventListColl
Definition: Event.h:147
const std::string cDetType_Strings[5]
Definition: Constants.h:509
cons_index_list< index_uni, nil_index_list > uni
static CovarianceBinUtility * Instance()
std::string to_string() const
Definition: ParameterSet.h:137
void CreateLocationsWithVariedOscParameters(cmf::Location const &nominalLoc, std::vector< cmf::Location > &vecLoc)