FNEXCovarianceMatrixMaker_plugin.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////
2 /// \brief Create a set of Covariance matrices to be used in fits
3 /// \author brebel@fnal.gov
4 ///////////////////////////////////////////////////////////////////////
5 #include <cmath>
6 #include <string>
7 #include <vector>
8 #include <unordered_map>
9 
10 // Framework includes
17 #include "fhiclcpp/ParameterSet.h"
19 
20 // NOvA includes
21 #include "NovaDAQConventions/DAQConventions.h"
22 #include "FNEX/core/Event.h"
23 #include "FNEX/core/VarVals.h"
29 
30 #include "TH1.h"
31 #include "TH2.h"
32 #include "TRandom3.h"
33 
34 namespace fnex {
35 
37  public:
38  explicit CovarianceMatrixMaker(fhicl::ParameterSet const& pset);
39  virtual ~CovarianceMatrixMaker();
40 
41  // Plugins should not be copied or assigned.
46 
47  void readResults (art::Results const& r) override;
48  void writeResults(art::Results & r) override;
49  void clear() override;
50  void beginJob() override;
51 
52  void reconfigure (fhicl::ParameterSet const&pset);
53 
54  private:
55 
57  void FillSpectrum (TH1D* spectrumBin,
58  std::map<int, TH1D*> & spectraEnergy);
59 
60  std::unique_ptr<fnex::CovarianceBinUtility> fBinUtil; ///< useful functions for mapping bins to energy and back
62  int fIterations; ///< number of iterations for picking sigmas
63  TRandom3 fRand; ///< RNG for picking the variation in parameters
64  fhicl::ParameterSet fManipulatorParameters; ///< the fake point parameters
65  fnex::EventListMap fEventLists; ///< the lists of events we are going to use, data and MC
66  TH2D* fCovarianceMatrix; ///< Covariance matrix
67  TH1D* fNominalSpectrum; ///< Nominal spectrum for each detector/beam/selection
68  std::map<int, TH1D*> fNominalEnergySpectra; ///< Nominal energy spectrum for each detector/file/beam/selection
69  TH1D* fShiftedSpectrum; ///< Shifted spectrum for each detector/beam/selection
70  std::map<int, TH1D*> fShiftedEnergySpectra; ///< Shifted energy spectrum for each detector/file/beam/selection
71  std::map<std::string, TH1D*> fShifts; ///< Shifts selected for each parameter adjusted
73  std::vector<std::string> fSystParameterNames; ///< names of parameters we are floating from the input point
74  };
75 
76  //......................................................................
78  : fRand(TRandom3(0))
79  {
80  this->reconfigure(pset);
81 
82  return;
83  }
84 
85  //......................................................................
87  {
88  }
89 
90  //......................................................................
92  {
93  // create the manipulator
94  fManipulatorParameters = pset.get< fhicl::ParameterSet >("EventListManipulator");
95  fIterations = pset.get< int >("Iterations", 1000);
96  fInputPointParameters = pset.get< fhicl::ParameterSet >("InitialGuessInputPoint");
97  auto randSeed = pset.get< unsigned long int >("RandomSeed", 0);
98 
99  if(randSeed > 0) fRand.SetSeed(randSeed);
100 
101  // get data samples to use from the manipulatorParameters
102  fBinUtil = std::make_unique<fnex::CovarianceBinUtility>();
103 
104  auto systematicParameters = fInputPointParameters.get< fhicl::ParameterSet >("SystematicPullInitial");
105  auto parNames = systematicParameters.get_names();
106 
107  fhicl::ParameterSet systPars;
108  for(const auto &name : parNames){
109  systPars = systematicParameters.get<fhicl::ParameterSet>(name);
110 
111  // don't bother with any systematics parameters that have been set to fixed
112  if(systPars.get<bool>("Fixed")) continue;
113 
114  fSystParameterNames.push_back(systPars.get<std::string>("Parameter"));
115  }
116 
117  return;
118  }
119 
120  //......................................................................
121  // Method to clear out the collections of data products after the
122  // writeResults method is called.
124  {
125  fEventLists.clear();
126 
127  return;
128  }
129 
130 
131  //......................................................................
133  {
134  // Create histograms for the covariance matrix, the nominal spectrum and
135  // the shifted values
136  // Each axis accounts for both detectors, both beam types and all desired selections
137  // Don't break out the interaction types or file types because they aren't known in
138  // the data
139  // The ND has 4 numu spectra and 2 nue
140  // The FD has 4 numu spectra and 2 nue and 1 peripheral
141  // For each beam type I need 8 * fBinUtil->NuMuHighEdges().size(),
142  // 4 * fBinUtil->NuEHighEdges.size() + 1
143  int numBins = 2 * (4 * fBinUtil->NuEHighEdges().size() +
144  8 * fBinUtil->NuMuHighEdges().size() +
145  1);
146  fCovarianceMatrix = fTFS->make<TH2D>("CovarianceMatrix",
147  ";Energy Bin;Energy Bin",
148  numBins,
149  1,
150  numBins + 1,
151  numBins,
152  1,
153  numBins + 1);
154 
155  fNominalSpectrum = fTFS->make<TH1D>("NominalBin",
156  ";Energy Bin;Events",
157  numBins,
158  1,
159  numBins + 1);
160 
161  fShiftedSpectrum = fTFS->make<TH1D>("ShiftedBin",
162  ";Energy Bin;Events",
163  numBins,
164  1,
165  numBins + 1);
166 
167 
168  std::stringstream titles;
169  titles << "Seed: "
170  << fRand.GetSeed()
171  << ";Shifts;Iterations";
172  for(auto name : fSystParameterNames){
173  fShifts[name] = fTFS->make<TH1D>((name+"Shifts").c_str(),
174  titles.str().c_str(),
175  500,
176  -5.,
177  5);
178  }
179 
180 
181  // make histograms for the actual energy spectrum too
182  // first have to determine the boundaries
183  int key = 0;
184  int bdf = 0;
185  double maxe = 0.;
188  std::set<fnex::FileType_t> fileTypes({fnex::kBeam,
189  fnex::kSwap,
192 
193  // loop over beams
194  for(int b = 0; b < 2; ++b){
195  // have to have the same bins in the ND and FD
196  for(int d = 0; d < 2; ++d){
197 
198  // loop over file types - the ND only cares about
199  // the nonswap files, FD has nonswap, fluxswap and tauswap
200  for(auto fitr : fileTypes){
201 
202  // in the ND, we only care about the nonswap file type
203  if(d == 0 && fitr != fnex::kBeam) continue;
204 
205  bdf = (b + 1) * 1000 + (d + 1) * 100 + (fitr + 1) * 10;
206 
208 
209  // MidPID selection is not used after 2017
210  // We only have peripheral samples in the FD
212  (e == fnex::kNuESelectionPeripheral && d == 0)) continue;
213 
214  key = bdf + e;
215 
216  maxe = (e == fnex::kNuESelectionPeripheral) ? 10. : fBinUtil->NuEHighEdges().rbegin()->first;
217  name = (fnex::cBeamType_Strings[b] +
219  fnex::kFileTypeStrings[fitr] +
221 
222  title = ";Energy (GeV);Events/10^{20} POT";
223 
224  fNominalEnergySpectra[key] = fTFS->make<TH1D>(name.c_str(), title.c_str(),
225  fBinUtil->NuEHighEdges().size() + 1,
226  0,
227  maxe);
228  } // end loop over nue selections
229 
230  // now the numu bins
231  for(size_t m = fnex::kNuMuSelectionQ1; m < fnex::kNuMuSelectionQ4 + 1; ++m){
232 
233  key = bdf + m;
234 
235  name = (fnex::cBeamType_Strings[b] +
237  fnex::kFileTypeStrings[fitr] +
239 
240  title = (d < 1) ? ";Energy (GeV);10^{3} Events/10^{20} POT" : ";Energy (GeV);Events/10^{20} POT";
241 
242  fNominalEnergySpectra[key] = fTFS->make<TH1D>(name.c_str(), title.c_str(),
243  fBinUtil->NuMuHighEdges().size() + 1,
244  0,
245  fBinUtil->NuMuHighEdges().rbegin()->first);
246  } // end loop over numu selections
247 
248  } // end loop over file types
249  } // end loop over detectors
250  } // end loop over beam types
251 
252  LOG_DEBUG("CovarianceMatrixMaker")
253  << "histograms made";
254 
255  // get the lists of events from the TTree and put them into a map
256  // of event lists
257  std::unique_ptr<fnex::EventListManipulator> manipulator = std::make_unique<fnex::EventListManipulator>(fManipulatorParameters,
259 
260  // we only want the MC lists
261  fEventLists = manipulator->DeserializeMC();
262 
263  return;
264  }
265 
266  //......................................................................
268  {
269  // The shifter and weighter was initialized in the EventListManipulator::Deserialize() method
270 
271  // fill the nominal histograms first so we don't do this for every time we
272  // change the systematic parameters
274 
275  LOG_DEBUG("CovarianceMatrixMaker")
276  << "nominal spectrum filled";
277 
278  return;
279  }
280 
281  // Function to fill the spectra for each meta data type
282  //......................................................................
283  void CovarianceMatrixMaker::FillSpectrum(TH1D* spectrumBin,
284  std::map<int, TH1D*> & spectraEnergy)
285  {
286  // first reset the spectra
287  spectrumBin ->Reset();
288  for(auto hitr : spectraEnergy)
289  hitr.second->Reset();
290 
291  // declare some variables we will be using a lot in the for loops below
292  double weight = 1.;
293  double energy = 0.;
294  int key = 0;
295  int bin = 0;
296  int offset = 0;
297  double normalization = 1.;
298  double potForNorm = 1.e8; // spill summary POT values are in units of 10^12 POT,
299  // so this normalization would be to 10^20 POT
300 
301  // loop over each list and event in the list
302  for(auto itr : fEventLists){
303  auto const md = itr.first;
304 
305  normalization = potForNorm/itr.second.ListSpillSummary().goodPOT;
306 
307  // normalize the numu selected ND down to a reasonable number of events
308  // to have about the same peak as the highest fd bins. In the end we
309  // are just looking to see relative changes in each bin anyway.
310  if(md.detector == novadaq::cnv::kNEARDET){
311  if(md.IsNuMuSelected()) normalization *= 1.e-3;
312  else if(md.IsNuESelected() ) normalization *= 1.e-2;
313  }
314 
315  key = fBinUtil->MetaDataToKey(md);
316 
317  // figure out the offsets
318  offset = fBinUtil->KeyToOffset(key);
319 
320  LOG_DEBUG("CovarianceMatrixMaker")
321  << md.ToString()
322  << " key: "
323  << fBinUtil->MetaDataToKey(md)
324  << " bin offset: "
325  << offset;
326 
327  // add the file type to the key for the energy spectra
328  key += (md.fileType + 1) * 10;
329 
330  for(auto const& evt : itr.second){
331 
332  energy = evt->dataVals->val_at(fnex::kNu_RecoE, md);
333  bin = fBinUtil->EnergyToBin(energy, md);
334 
335  // if bin is negative then the energy is out of range for
336  // the selection type and we'll ignore it.
337  if(bin < 0) continue;
338 
340 
341  // now fill the spectrum
342  spectrumBin->Fill(bin + offset, weight * normalization);
343 
344  if(spectraEnergy.count(key) > 0)
345  spectraEnergy[key]->Fill(energy, weight * normalization);
346 
347  } // end loop over events
348  } // end loop over lists of events
349 
350  return;
351  }
352 
353 
354  // This function will pick a random value for the specified parameter
355  // from a guassian distribution. The random value will be within +/- 1 sigma
356  // of the nominal value
357  //......................................................................
359  {
360 
361  // we will pick a random value from a unit gaussian, ie gaussian with
362  // sigma = 1. That will be the fractional change on the sigma we want
363  // to apply to the new input point
364  double sigFrac = 0.;
365 
366  LOG_DEBUG("CovarianceMatrixMaker")
367  << "sigfrac = "
368  << sigFrac
369  << " original point "
370  << point;
371 
372  ParameterDetMap systematicPulls = point.SystematicPulls();
373 
374  for(auto name : fSystParameterNames) {
377 
378  // each systematic parameter should have a different random number
379  // change in sigma so that we don't artificially make them change
380  // in a correlated way if doing more than one parameter in a job
381  sigFrac = fRand.Gaus();
382 
383  if((name.find("CalibShape") != std::string::npos ||
384  name.find("Cherenkov") != std::string::npos) &&
385  sigFrac < 0.) sigFrac *= -1.;
386 
387  fShifts[name]->Fill(sigFrac);
388 
389  LOG_DEBUG("CovarianceMatrixMaker")
390  << " parameter: "
391  << name
392  << " sigFrac "
393  << sigFrac;
394 
395  if(systematicPulls.count(tempPairFar) != 0){
396  systematicPulls[tempPairFar] = sigFrac;
397 
398  LOG_DEBUG("CovarianceMatrixMaker")
399  << "FD name: "
400  << name
401  << " "
402  << systematicPulls[tempPairFar];
403  }
404  if(systematicPulls.count(tempPairNear) != 0){
405  systematicPulls[tempPairNear] = sigFrac;
406 
407  LOG_DEBUG("CovarianceMatrixMaker")
408  << "ND name: "
409  << name
410  << " "
411  << systematicPulls[tempPairFar];
412  }
413  } // end loop over systematic parameters to change
414 
415  // Make the new input point
416  //
417  point = fnex::InputPoint(point.ParameterSpaceLocation(),
418  systematicPulls,
419  point.FixedParameters(),
420  point.FixedSystematics() );
421 
422  LOG_DEBUG("CovarianceMatrixMaker")
423  << " new point "
424  << point;
425 
427 
428  return;
429  }
430 
431  //......................................................................
433  {
434  // declare some variables we will be using a lot in the for loops below
435  int nbins = fNominalSpectrum->GetNbinsX();
436  double diffI = 0.;
437  double diffJ = 0.;
438  double iCenter = 0.;
439 
441 
442  // loop over all the events and fill the matrices for many different
443  // values of the systematic parameters
444  for(int k = 0; k < fIterations; ++k){
445 
446  LOG_VERBATIM("CovarianceMatrixMaker")
447  << "iteration: "
448  << k;
449 
450  this->SetSystematicParameter(point);
451 
452  // this is a new shift for the parameter, so reset and refill the spectrum
454 
455  // now to fill the covariance matrices
456  // For ROOT histograms bin 0 is the underflow bin, so start at bin 1
457  // if either diffI or diffJ are 0, then we don't need to fill
458  // the bin with a 0 (ie the product) and can continue on
459  for(int i = 1; i < nbins + 1; ++i){
460  iCenter = fNominalSpectrum->GetBinCenter(i);
461  diffI = (fShiftedSpectrum->GetBinContent(i) - fNominalSpectrum->GetBinContent(i));
462  if(diffI != 0.) diffI /= fNominalSpectrum->GetBinContent(i);
463  else continue;
464 
465  for(int j = 1; j < nbins + 1; ++j){
466  diffJ = (fShiftedSpectrum->GetBinContent(j) - fNominalSpectrum->GetBinContent(j));
467  if (diffJ != 0.) diffJ /= fNominalSpectrum->GetBinContent(j);
468  else continue;
469  /*
470  if(i == 1)
471  LOG_VERBATIM("CovarianceMatrixMaker")
472  << "i: "
473  << i
474  << " shifted i: "
475  << fShiftedSpectrum->GetBinContent(i)
476  << " nominal i: "
477  << fNominalSpectrum->GetBinContent(i)
478  << " j: "
479  << j
480  << " shifted j: "
481  << fShiftedSpectrum->GetBinContent(j)
482  << " nominal j: "
483  << fNominalSpectrum->GetBinContent(j)
484  << " diffI: "
485  << diffI
486  << " diffJ: "
487  << diffJ
488  << " product: "
489  << diffI * diffJ;
490  */
491  fCovarianceMatrix->Fill(iCenter, fNominalSpectrum->GetBinCenter(j), diffI * diffJ);
492  } // end loop over j bins of spectrum
493  } // end loop over i bins of spectrum
494 
495  } // end loop over the iterations on the systematic parameter
496 
497  // protect against very low and very high covariances
498  // these usually happen because nentries in these bins for nominal
499  // and shifted spectra are identially 0
500  for (int i = 1; i < nbins + 1; ++i){
501  for (int j = 1; j < fCovarianceMatrix->GetNbinsY()+1; j++){
502  double binContent = fCovarianceMatrix->GetBinContent(i,j);
503  if (std::abs(binContent) < std::pow(10,-14) || std::abs(binContent) > std::pow(10,14) || isnan(binContent)) {
504  // LOG_VERBATIM("CovarianceMatrixMaker")
505  // << "bin " << i
506  // << ", covariance: " << binContent
507  // << "\n -- nominal i: " << fNominalSpectrum->GetBinContent(i)
508  // << "\n -- shifted i: " << fShiftedSpectrum->GetBinContent(i)
509  // << "\n -- nominal minus shifted i: " << fNominalSpectrum->GetBinContent(i) - fShiftedSpectrum->GetBinContent(i)
510  // << "\n -- nominal j: " << fNominalSpectrum->GetBinContent(j)
511  // << "\n -- shifted j: " << fShiftedSpectrum->GetBinContent(j)
512  // << "\n -- nominal minus shifted j: " << fNominalSpectrum->GetBinContent(j) - fShiftedSpectrum->GetBinContent(j)
513  // << "\n -- diffI*diffJ" << (fNominalSpectrum->GetBinContent(i)-fShiftedSpectrum->GetBinContent(i))*(fNominalSpectrum->GetBinContent(j)-fShiftedSpectrum->GetBinContent(j));
514 
515  fCovarianceMatrix->SetBinContent(i,j,0);
516  fCovarianceMatrix->SetBinContent(j,i,0);
517  }
518  }
519  }
520 
521  return;
522  }
523 
524 
526 
527 } // end fnex namespace
void SetCurrentVals(fnex::InputPoint const &curPoint)
static std::string GetName(int id)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
std::map< int, TH1D * > fNominalEnergySpectra
Nominal energy spectrum for each detector/file/beam/selection.
const XML_Char * name
Definition: expat.h:151
static const unsigned char kNu_RecoE
Definition: VarVals.h:36
int fIterations
number of iterations for picking sigmas
const Var weight
fnex::EventListMap fEventLists
the lists of events we are going to use, data and MC
constexpr T pow(T x)
Definition: pow.h:75
std::map< std::string, TH1D * > fShifts
Shifts selected for each parameter adjusted.
TRandom3 fRand
RNG for picking the variation in parameters.
std::map< int, TH1D * > fShiftedEnergySpectra
Shifted energy spectrum for each detector/file/beam/selection.
std::map< fnex::MetaData, fnex::EventList > EventListMap
Definition: Event.h:186
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
std::vector< std::string > fSystParameterNames
names of parameters we are floating from the input point
Create a list of fnex::Events to be used in fits.
static ShifterAndWeighter * Instance()
float abs(float number)
Definition: d0nt_math.hpp:39
std::unique_ptr< fnex::CovarianceBinUtility > fBinUtil
useful functions for mapping bins to energy and back
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
const int nbins
Definition: cellShifts.C:15
Far Detector at Ash River, MN.
art::ServiceHandle< art::TFileService > fTFS
TFileService.
CovarianceMatrixMaker & operator=(CovarianceMatrixMaker const &)=delete
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
double Weight(fnex::EventPtr const &curEvent, fnex::MetaData const &md, fnex::WeightType_t const wgtType=kAllWgt)
#define DEFINE_ART_RESULTS_PLUGIN(klass)
T get(std::string const &key) const
Definition: ParameterSet.h:231
void SetSystematicParameter(fnex::InputPoint &point)
int evt
const std::string kFileTypeStrings[8]
Definition: Constants.h:386
void writeResults(art::Results &r) override
double energy
Definition: plottest35.C:25
Near Detector in the NuMI cavern.
Float_t d
Definition: plot.C:236
const double j
Definition: BetheBloch.cxx:29
void reconfigure(fhicl::ParameterSet const &pset)
std::vector< std::string > get_names() const
float bin[41]
Definition: plottest35.C:14
const std::string cBeamType_Strings[4]
Definition: Constants.h:76
std::unordered_map< ParameterDet, double, ParameterDetHasher > ParameterDetMap
Definition: Constants.h:63
CovarianceMatrixMaker(fhicl::ParameterSet const &pset)
TH1D * fShiftedSpectrum
Shifted spectrum for each detector/beam/selection.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
fhicl::ParameterSet fInputPointParameters
input parameters
void readResults(art::Results const &r) override
const hit & b
Definition: hits.cxx:21
std::pair< std::string, novadaq::cnv::DetId > ParameterDet
Definition: Constants.h:38
void FillSpectrum(TH1D *spectrumBin, std::map< int, TH1D * > &spectraEnergy)
TRandom3 r(0)
fhicl::ParameterSet fManipulatorParameters
the fake point parameters
#define LOG_VERBATIM(category)
const std::string cSelectionType_Strings[11]
Definition: Constants.h:101
Float_t e
Definition: plot.C:35
TH1D * fNominalSpectrum
Nominal spectrum for each detector/beam/selection.
ParameterMap const SystematicPulls(novadaq::cnv::DetId const &det) const
Definition: InputPoint.cxx:364