CovarianceMatrixMaker_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 <sstream>
9 
10 // Framework includes
17 #include "fhiclcpp/ParameterSet.h"
19 
20 // NOvA includes
21 #include "NovaDAQConventions/DAQConventions.h"
30 
31 #include "TH1.h"
32 #include "TH2.h"
33 #include "TRandom3.h"
34 
35 namespace cmf {
36 
38  public:
39  explicit CovarianceMatrixMaker(fhicl::ParameterSet const& pset);
40  ~CovarianceMatrixMaker() = default;
41 
42  // Plugins should not be copied or assigned.
47 
48  void readResults (art::Results const& r) override;
49  void writeResults(art::Results & r) override;
50  void clear() override;
51  void beginJob() override;
52 
53  void reconfigure (fhicl::ParameterSet const&pset);
54 
55  private:
56 
58  void FillSpectrum (cmf::EventListColl const& eventLists,
59  cmf::ExposureMap const& exposureMap,
60  cmf::Location const& loc,
61  TH1D* spectrumBin,
62  std::map<long, TH1D*> & spectraEnergy);
63 
65  int fIterations; ///< number of iterations for picking sigmas
66  TRandom3 fRand; ///< RNG for picking the variation in parameters
67  fhicl::ParameterSet fManipulatorParameters; ///< the fake point parameters
68  TH2D* fSystCovarianceMatrix; ///< Syst Covariance matrix
69  TH1D* fNominalSpectrum; ///< Nominal spectrum for each detector/beam/selection
70  std::map<long, TH1D*> fNominalEnergySpectra; ///< Nominal energy spectrum for each detector/file/beam/selection
71  TH1D* fPlus1SigmaSpectrum; ///< Plus 1 Sigma shifted spectrum for each detector/beam/selection
72  std::map<long, TH1D*> fPlus1SigmaEnergySpectra; ///< Plus 1 Sigma shifted energy spectrum for each detector/file/beam/selection
73  TH1D* fMinus1SigmaSpectrum; ///< Minus 1 Sigma shifted spectrum for each detector/beam/selection
74  std::map<long, TH1D*> fMinus1SigmaEnergySpectra; ///< Minus 1 Sigma shifted energy spectrum for each detector/file/beam/selection
75  TH1D* fShiftedSpectrum; ///< Shifted spectrum for each detector/beam/selection
76  TH2D* fShiftedSpectrumByUniverse; ///< Shifted spectrum as a function of universe
77  std::map<long, TH1D*> fShiftedEnergySpectra; ///< Shifted energy spectrum for each detector/file/beam/selection
78  std::map<std::string, TH1D*> fShifts; ///< Shifts selected for each parameter adjusted
79  TH2D* fShiftWeights; ///< Weights from shifts selected for each parameter adjusted
80  fhicl::ParameterSet fSAWParameters; ///< shifter and weighter parameters
81  };
82 
83  //......................................................................
85  : fRand(TRandom3(0))
86  , fSystCovarianceMatrix(nullptr)
87  , fNominalSpectrum(nullptr)
88  , fPlus1SigmaSpectrum(nullptr)
89  , fMinus1SigmaSpectrum(nullptr)
90  , fShiftedSpectrum(nullptr)
92  , fShiftWeights(nullptr)
93  {
94  this->reconfigure(pset);
95  }
96 
97  //......................................................................
99  {
102  cmf::ParameterUtility::Instance() ->Initialize(pset.get< fhicl::ParameterSet >("ParametersToUse" ));
103 
104  fManipulatorParameters = pset.get< fhicl::ParameterSet >("EventListManipulator");
105  fIterations = pset.get< int >("Iterations", 1000);
106  fSAWParameters = pset.get< fhicl::ParameterSet >("ShifterAndWeighterParameters");
107  auto randSeed = pset.get< unsigned long int >("RandomSeed", 0);
108 
110  fSAWParameters);
111 
112  if(randSeed > 0) fRand.SetSeed(randSeed);
113  }
114 
115  //......................................................................
116  // Method to clear out the collections of data products after the
117  // writeResults method is called.
119  {
120  }
121 
122 
123  //......................................................................
125  {
126  // Create histograms for the covariance matrix, the nominal spectrum and
127  // the shifted values
128  // Each axis accounts for both detectors, both beam types and all desired selections
129  // Don't break out the interaction types or file types because they aren't known in
130  // the data
132 
133  fSystCovarianceMatrix = fTFS->make<TH2D>("CovarianceMatrix",
134  ";Energy Bin;Energy Bin",
135  numBins,
136  1,
137  numBins + 1,
138  numBins,
139  1,
140  numBins + 1);
141 
142  fNominalSpectrum = fTFS->make<TH1D>("NominalBin",
143  ";Energy Bin;Events",
144  numBins,
145  1,
146  numBins + 1);
147 
148  fShiftedSpectrum = fTFS->make<TH1D>("ShiftedBin",
149  ";Energy Bin;Events",
150  numBins,
151  1,
152  numBins + 1);
153 
154  fPlus1SigmaSpectrum = fTFS->make<TH1D>("Plus1SigmaSpectrum",
155  ";Energy Bin;Events",
156  numBins,
157  1,
158  numBins + 1);
159 
160  fMinus1SigmaSpectrum = fTFS->make<TH1D>("Minus1SigmaSpectrum",
161  ";Energy Bin;Events",
162  numBins,
163  1,
164  numBins + 1);
165 
166  fShiftedSpectrumByUniverse = fTFS->make<TH2D>("ShiftedBinByUniverse",
167  ";Energy Bin;Universe",
168  numBins,
169  1,
170  numBins + 1,
171  fIterations,
172  0,
173  fIterations);
174 
175  std::stringstream shifttitles;
176  shifttitles << "Seed: "
177  << fRand.GetSeed()
178  << ";Shifts;Iterations";
179  for(auto const& itr : cmf::ParameterUtility::Instance()->SysParNames()){
180  fShifts[itr] = fTFS->make<TH1D>((itr+"Shifts").c_str(),
181  shifttitles.str().c_str(),
182  500,
183  -5.,
184  5);
185 
186  }
187 
188  // histogram of the weights as a function of energy bin
189  std::stringstream shiftwgttitles;
190  shiftwgttitles << "Seed: "
191  << fRand.GetSeed()
192  << ";Energy Bin;Total Shift Weights";
193  fShiftWeights = fTFS->make<TH2D>("TotalShiftWeights",
194  shiftwgttitles.str().c_str(),
195  numBins,
196  1,
197  numBins + 1,
198  500,
199  -1.,
200  1);
201 
202 
203  long key = 0;
206 
208  std::vector<double> bins;
209 
210  // loop over the selections to use
211  auto const& dbs = cmf::SelectionUtility::Instance()->SelectionsToUse();
212 
213  for(auto const& itr : dbs){
214  md.period = (itr.BeamType() == cmf::kFHC) ? 1 : 4;
215  md.detector = itr.Detector();
216 
217  title = (md.detector == cmf::kNEARDET) ? ";Energy (GeV);10^{3} Events/10^{20} POT" : ";Energy (GeV);Events/10^{20} POT";
218 
219  // loop over the selection types in each detector
220  for(auto const& sel : itr.Selections()){
221 
222  // no peripheral nue selection in the ND
223  if(md.detector == cmf::kNEARDET &&
224  sel == cmf::kNuESelectionPeripheral) continue;
225 
226  md.selectionType = sel;
227 
228  auto const& highEdges = cmf::CovarianceBinUtility::Instance()->SelectionHighEdges(md);
229  bins.clear();
230  bins.insert(bins.begin(), cmf::CovarianceBinUtility::Instance()->SelectionLowEdge(md));
231 
232  for(auto const& itr : highEdges ) bins.push_back(itr.first);
233 
234  key = md.DetectorBeamSelectionKey();
235 
236  name = (cmf::cBeamType_Strings[itr.BeamType()] +
237  cmf::cDetType_Strings[itr.Detector()] +
239 
240 // LOG_VERBATIM("tpm")
241 // << "key: " << cmf::KeyToString(key) << " bins: ";
242 // for (auto& b : bins)
243 // LOG_VERBATIM("tmp") << b;
244 
245  fNominalEnergySpectra[key] = fTFS->make<TH1D>(name.c_str(),
246  title.c_str(),
247  bins.size() - 1,
248  &bins[0]);
249 
250  } // end loop over selections
251  } // end loop over dbs
252 
253  LOG_DEBUG("CovarianceMatrixMaker")
254  << "histograms made";
255 
256  }
257 
258  //......................................................................
260  {
261  }
262 
263  // Function to fill the spectra for each meta data type
264  //......................................................................
266  cmf::ExposureMap const& exposureMap,
267  cmf::Location const& loc,
268  TH1D* spectrumBin,
269  std::map<long, TH1D*> & spectraEnergy)
270  {
271  // first reset the spectra
272  spectrumBin->Reset();
273  for(auto hitr : spectraEnergy) hitr.second->Reset();
274 
275  // now call the SpectrumTools.
276  std::vector<float> spectrum;
277  std::vector<float> count;
278 
280  cmf::FillSpectrum(eventLists,
281  exposureMap,
282  spectrum,
283  count);
284 
285  // now fill the spectrum histogram
286  long key;
287  for(size_t b = 0; b < spectrum.size(); ++b){
288  spectrumBin->Fill(b, spectrum[b]);
289 
291  if(spectraEnergy.count(key) > 0)
292  spectraEnergy[key]->Fill(cmf::CovarianceBinUtility::Instance()->BinToEnergy(b),
293  spectrum[b]);
294 
295  // fill the histogram of weights as a function of energy bin
296  if(count[b] > 0) fShiftWeights->Fill(b, spectrum[b]/count[b]);
297  } // end loop over spectrum bins
298  }
299 
300  // This function will pick a random value for the specified parameter
301  // from a guassian distribution. The random value will be within +/- 1 sigma
302  // of the nominal value
303  //......................................................................
305  {
306 
307  // we will pick a random value from a unit gaussian, ie gaussian with
308  // sigma = 1. That will be the fractional change on the sigma we want
309  // to apply to the new input point
310  double sigFrac = 0.;
311 
312  LOG_DEBUG("CovarianceMatrixMaker")
313  << "sigfrac = "
314  << sigFrac;
315 
316  for(auto const& name : cmf::ParameterUtility::Instance()->SysParNames()){
317 
318  // each systematic parameter should have a different random number
319  // change in sigma so that we don't artificially make them change
320  // in a correlated way if doing more than one parameter in a job
321  sigFrac = fRand.Gaus();
322 
323  // The calibration and light level systematic uncertainties are
324  // really only defined for a 1 sigma shift because their excursions
325  // are calculated from special data sets with the shift applied
326  // Also pull them from a uniform distribution as they are definitely
327  // not Gaussian.
328  if(name.find("lightmodel") != std::string::npos ||
329  name.find("ckv-proton" ) != std::string::npos ||
330  name.find("Calibration") != std::string::npos ){
331  sigFrac = fRand.Rndm();
332 
333  LOG_DEBUG("CovarianceMatrixMaker")
334  << "pulling from uniform distribution "
335  << sigFrac;
336  }
337 
338  fShifts[name]->Fill(sigFrac);
339 
340  LOG_DEBUG("CovarianceMatrixMaker")
341  << " parameter: "
342  << name
343  << " sigFrac "
344  << sigFrac;
345 
346  loc.SetParameterValue(name, sigFrac, cmf::kFARDET);
347  LOG_DEBUG("CovarianceMatrixMaker")
348  << "FD name: "
349  << name
350  << " "
351  << loc.FDLocation().find(name)->second.Value();
352 
353  loc.SetParameterValue(name, sigFrac, cmf::kNEARDET);
354  LOG_DEBUG("CovarianceMatrixMaker")
355  << "ND name: "
356  << name
357  << " "
358  << loc.NDLocation().find(name)->second.Value();
359 
360  } // end loop over systematic parameters to change
361  }
362 
363  //......................................................................
365  {
366  // declare some variables we will be using a lot in the for loops below
367  int nbins = fNominalSpectrum->GetNbinsX();
368  double diffI;
369  double diffJ;
370  double iCenter;
371 
372  // containers for different oscillation and systematic parameter locations
374 
375  // get the lists of events from the TTree and put them into a map
376  // of event lists
378 
379  // we only want the MC lists
380  cmf::EventListColl eventLists;
381  manipulator.Deserialize(eventLists, cmf::kMC);
382 
383  cmf::ExposureMap exposureMap = manipulator.ExposureMap();
384 
385  // fill the nominal histograms first so we don't do this for every time we
386  // change the systematic parameters
387  this->FillSpectrum(eventLists,
388  exposureMap,
389  loc,
392 
393  LOG_DEBUG("CovarianceMatrixMaker")
394  << "nominal spectrum filled";
395 
396  // now fake input points for the +/- 1 sigma shifts
397 
398  auto parNames = cmf::ParameterUtility::Instance()->SysParNames();
399 
400  // plus 1 sigma
401  for(auto const& itr : parNames){
402  loc.SetParameterValue(itr, 1, cmf::kNEARDET);
403  loc.SetParameterValue(itr, 1, cmf::kFARDET);
404  }
405 
406  this->FillSpectrum(eventLists,
407  exposureMap,
408  loc,
411 
412  // minus 1 sigma
413  for(auto const& itr : parNames){
414  loc.SetParameterValue(itr, -1, cmf::kNEARDET);
415  loc.SetParameterValue(itr, -1, cmf::kFARDET);
416  }
417 
418  this->FillSpectrum(eventLists,
419  exposureMap,
420  loc,
423 
424  // loop over all the events and fill the matrices for many different
425  // values of the systematic parameters
426  for(int k = 0; k < fIterations; ++k){
427 
428  LOG_VERBATIM("CovarianceMatrixMaker")
429  << "iteration: "
430  << k;
431 
433 
435 
436  // this is a new shift for the parameter, so reset and refill the spectrum
437  this->FillSpectrum(eventLists,
438  exposureMap,
439  loc,
442 
443  // now to fill the covariance matrices
444  // For ROOT histograms bin 0 is the underflow bin, so start at bin 1
445  // if either diffI or diffJ are 0, then we don't need to fill
446  // the bin with a 0 (ie the product) and can continue on
447  for(int i = 1; i < nbins + 1; ++i){
448  iCenter = fNominalSpectrum->GetBinCenter(i);
449  diffI = (fShiftedSpectrum->GetBinContent(i) - fNominalSpectrum->GetBinContent(i));
450  if(diffI != 0.) diffI /= fNominalSpectrum->GetBinContent(i);
451  else continue;
452 
453  for(int j = 1; j < nbins + 1; ++j){
454  diffJ = (fShiftedSpectrum->GetBinContent(j) - fNominalSpectrum->GetBinContent(j));
455  if (diffJ != 0.) diffJ /= fNominalSpectrum->GetBinContent(j);
456  else continue;
457  /*
458  if(i == 1)
459  LOG_VERBATIM("CovarianceMatrixMaker")
460  << "i: "
461  << i
462  << " shifted i: "
463  << fShiftedSpectrum->GetBinContent(i)
464  << " nominal i: "
465  << fNominalSpectrum->GetBinContent(i)
466  << " j: "
467  << j
468  << " shifted j: "
469  << fShiftedSpectrum->GetBinContent(j)
470  << " nominal j: "
471  << fNominalSpectrum->GetBinContent(j)
472  << " diffI: "
473  << diffI
474  << " diffJ: "
475  << diffJ
476  << " product: "
477  << diffI * diffJ;
478  */
479  fSystCovarianceMatrix->Fill(iCenter, fNominalSpectrum->GetBinCenter(j), diffI * diffJ);
480 
481  } // end loop over j bins of spectrum
482  } // end loop over i bins of spectrum
483 
484  // loop over each bin and fill fShiftedSpectrumByUniverse
485  for (int i_bin = 1; i_bin < fShiftedSpectrum->GetNbinsX(); ++i_bin){
486  fShiftedSpectrumByUniverse->SetBinContent(i_bin, k + 1, fShiftedSpectrum->GetBinContent(i_bin));
487  }
488  } // end loop over the iterations on the systematic parameter
489 
490  // if we're not running any systematic uncertainties,
491  // then assume we want stat only
492  // n.b. we're usually running the xsec weights
493  // so take that into account
494  float binContent;
495  if (parNames.empty() || (parNames.size() == 1 && *(parNames.begin()) == "XSecCVWgt")){
496  for (int i_bin = 0; i_bin < fNominalSpectrum->GetNbinsX(); ++i_bin){
497 
498  binContent = fNominalSpectrum->GetBinContent(i_bin);
499 
500  if(binContent == 0){
501  fSystCovarianceMatrix->SetBinContent(i_bin,
502  i_bin,
504  }
505  else {
506  fSystCovarianceMatrix->SetBinContent(i_bin,
507  i_bin,
508  1./binContent);
509  }
510  }
511  }
512 
513  // protect against very low and very high covariances
514  // these usually happen because nentries in these bins for nominal
515  // and shifted spectra are identially 0
516  for (int i = 1; i < nbins + 1; ++i){
517  for (int j = 1; j < fSystCovarianceMatrix->GetNbinsY()+1; ++j){
518  binContent = fSystCovarianceMatrix->GetBinContent(i, j);
519  if (std::abs(binContent) < std::pow(10,-14) ||
520  std::abs(binContent) > std::pow(10,14) ||
521  std::isnan(binContent) ){
522  // LOG_VERBATIM("CovarianceMatrixMaker")
523  // << "bin " << i
524  // << ", covariance: " << binContent
525  // << "\n -- nominal i: " << fNominalSpectrum->GetBinContent(i)
526  // << "\n -- shifted i: " << fShiftedSpectrum->GetBinContent(i)
527  // << "\n -- nominal minus shifted i: " << fNominalSpectrum->GetBinContent(i) - fShiftedSpectrum->GetBinContent(i)
528  // << "\n -- nominal j: " << fNominalSpectrum->GetBinContent(j)
529  // << "\n -- shifted j: " << fShiftedSpectrum->GetBinContent(j)
530  // << "\n -- nominal minus shifted j: " << fNominalSpectrum->GetBinContent(j) - fShiftedSpectrum->GetBinContent(j)
531  // << "\n -- diffI*diffJ" << (fNominalSpectrum->GetBinContent(i)-fShiftedSpectrum->GetBinContent(i))*(fNominalSpectrum->GetBinContent(j)-fShiftedSpectrum->GetBinContent(j));
532 
533  fSystCovarianceMatrix->SetBinContent(i, j, 0);
534  fSystCovarianceMatrix->SetBinContent(j, i, 0);
535  } // end if the bin content is acceptable
536  } // end loop over j bins
537  } // end loop over i bins
538  }
539 
541 
542 } // end cmf namespace
CovarianceMatrixMaker & operator=(CovarianceMatrixMaker const &)=delete
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
const XML_Char * name
Definition: expat.h:151
std::map< cmf::MetaData, cmf::SpillSummary > ExposureMap
Definition: Structs.h:215
int fIterations
number of iterations for picking sigmas
void writeResults(art::Results &r) override
const std::string cSelectionType_Strings[12]
Definition: Constants.h:93
void Initialize(fhicl::ParameterSet const &pset)
cmf::DetType_t detector
Definition: Structs.h:114
std::map< long, TH1D * > fNominalEnergySpectra
Nominal energy spectrum for each detector/file/beam/selection.
double const & SelectionLowEdge(cmf::MetaData const &md)
constexpr T pow(T x)
Definition: pow.h:75
long DetectorBeamSelectionKey() const
Definition: Structs.cxx:194
static SelectionUtility * Instance()
TH1D * fNominalSpectrum
Nominal spectrum for each detector/beam/selection.
long BinToKey(int const &bin, bool allSels=false)
TH2D * fShiftWeights
Weights from shifts selected for each parameter adjusted.
void reconfigure(fhicl::ParameterSet const &pset)
float abs(float number)
Definition: d0nt_math.hpp:39
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
std::vector< std::string > SysParNames() const
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}))
double BinToEnergy(int const &bin, bool allSels=false)
static ParameterUtility * Instance()
TH1D * fShiftedSpectrum
Shifted spectrum for each detector/beam/selection.
const int nbins
Definition: cellShifts.C:15
TH2D * fShiftedSpectrumByUniverse
Shifted spectrum as a function of universe.
std::map< long, TH1D * > fPlus1SigmaEnergySpectra
Plus 1 Sigma shifted energy spectrum for each detector/file/beam/selection.
cmf::SelectionType_t selectionType
Definition: Structs.h:116
fhicl::ParameterSet fManipulatorParameters
the fake point parameters
cmf::ParameterSpaceLoc const & NDLocation() const
Definition: Parameter.h:160
fhicl::ParameterSet fSAWParameters
shifter and weighter parameters
TH2D * fSystCovarianceMatrix
Syst Covariance matrix.
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
#define DEFINE_ART_RESULTS_PLUGIN(klass)
T get(std::string const &key) const
Definition: ParameterSet.h:231
cmf::ParameterSpaceLoc const & FDLocation() const
Definition: Parameter.h:161
cmf::ExposureMap const & ExposureMap() const
void Initialize(fhicl::ParameterSet const &pset)
const double j
Definition: BetheBloch.cxx:29
void FillSpectrum(cmf::EventListColl const &eventLists, cmf::ExposureMap const &exposureMap, std::vector< float > &spectrum, std::vector< float > &count, cmf::InteractionType_t intType, bool applyExposureNorm)
const Binning bins
static ShifterAndWeighter * Instance()
std::map< long, TH1D * > fMinus1SigmaEnergySpectra
Minus 1 Sigma shifted energy spectrum for each detector/file/beam/selection.
TH1D * fMinus1SigmaSpectrum
Minus 1 Sigma shifted spectrum for each detector/beam/selection.
TH1D * fPlus1SigmaSpectrum
Plus 1 Sigma shifted spectrum for each detector/beam/selection.
size_t TotalBins(bool allSels=false)
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
void FillSpectrum(cmf::EventListColl const &eventLists, cmf::ExposureMap const &exposureMap, cmf::Location const &loc, TH1D *spectrumBin, std::map< long, TH1D * > &spectraEnergy)
void readResults(art::Results const &r) override
Module to combine a set of results into a single file currently only does one data product type at a ...
Definition: Event.cxx:24
CovarianceMatrixMaker(fhicl::ParameterSet const &pset)
std::map< std::string, TH1D * > fShifts
Shifts selected for each parameter adjusted.
std::vector< cmf::EventList > EventListColl
Definition: Event.h:147
void InitShiftsAndWeightsToUse(cmf::Location const &loc, fhicl::ParameterSet const &parset)
T * make(ARGS...args) const
art::ServiceHandle< art::TFileService > fTFS
TFileService.
std::map< double, int > const & SelectionHighEdges(cmf::MetaData const &md)
const hit & b
Definition: hits.cxx:21
cmf::Location ParametersAsLocation()
TRandom3 r(0)
std::map< long, TH1D * > fShiftedEnergySpectra
Shifted energy spectrum for each detector/file/beam/selection.
const DetBeamSelSet & SelectionsToUse()
long period
Definition: Structs.h:118
#define LOG_VERBATIM(category)
void Initialize(fhicl::ParameterSet const &pset)
const std::string cDetType_Strings[5]
Definition: Constants.h:29
const std::string cBeamType_Strings[4]
Definition: Constants.h:51
void SetCurrentVals(cmf::Location const &loc)
void SetSystematicParameter(cmf::Location &loc)
TRandom3 fRand
RNG for picking the variation in parameters.
void SetParameterValue(std::string const &parName, double const &value, cmf::DetType_t const &det)
Definition: Parameter.cxx:236
static CovarianceBinUtility * Instance()
enum BeamMode string