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"
28 
29 #include "TH1.h"
30 #include "TH2.h"
31 #include "TRandom3.h"
32 
33 namespace cmf {
34 
36  public:
37  explicit CovarianceMatrixMaker(fhicl::ParameterSet const& pset);
38  ~CovarianceMatrixMaker() = default;
39 
40  // Plugins should not be copied or assigned.
45 
46  void readResults (art::Results const& r) override;
47  void writeResults(art::Results & r) override;
48  void clear() override;
49  void beginJob() override;
50 
51  void reconfigure (fhicl::ParameterSet const&pset);
52 
53  private:
54 
56  void FillSpectrum (cmf::EventListMap const & eventLists,
57  TH1D* spectrumBin,
58  std::map<long, TH1D*> & spectraEnergy);
59 
61  int fIterations; ///< number of iterations for picking sigmas
62  TRandom3 fRand; ///< RNG for picking the variation in parameters
63  fhicl::ParameterSet fManipulatorParameters; ///< the fake point parameters
64  TH2D* fSystCovarianceMatrix; ///< Syst Covariance matrix
65  TH1D* fNominalSpectrum; ///< Nominal spectrum for each detector/beam/selection
66  std::map<long, TH1D*> fNominalEnergySpectra; ///< Nominal energy spectrum for each detector/file/beam/selection
67  TH1D* fShiftedSpectrum; ///< Shifted spectrum for each detector/beam/selection
68  TH2D* fShiftedSpectrumByUniverse; ///< Shifted spectrum as a function of universe
69  std::map<long, TH1D*> fShiftedEnergySpectra; ///< Shifted energy spectrum for each detector/file/beam/selection
70  std::map<std::string, TH1D*> fShifts; ///< Shifts selected for each parameter adjusted
71  TH2D* fShiftWeights; ///< Weights from 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  , fSystCovarianceMatrix(nullptr)
80  , fNominalSpectrum(nullptr)
81  , fShiftedSpectrum(nullptr)
83  , fShiftWeights(nullptr)
84  {
85  this->reconfigure(pset);
86  }
87 
88  //......................................................................
90  {
91  // create the manipulator
92  fManipulatorParameters = pset.get< fhicl::ParameterSet >("EventListManipulator");
93  fIterations = pset.get< int >("Iterations", 1000);
94  fInputPointParameters = pset.get< fhicl::ParameterSet >("InitialGuessInputPoint");
95  auto randSeed = pset.get< unsigned long int >("RandomSeed", 0);
96 
97  if(randSeed > 0) fRand.SetSeed(randSeed);
98 
99  auto systematicParameters = fInputPointParameters.get< fhicl::ParameterSet >("SystematicPullInitial");
100  auto parNames = systematicParameters.get_names();
101 
102  fhicl::ParameterSet systPars;
103  for(const auto &name : parNames){
104  LOG_VERBATIM("cmf::CovarianceMatrixMaker")
105  << "Running with systematic parameter: " << name;
106  systPars = systematicParameters.get<fhicl::ParameterSet>(name);
107 
108  // don't bother with any systematics parameters that have been set to fixed
109  if(systPars.get<bool>("Fixed")) continue;
110 
111  fSystParameterNames.push_back(systPars.get<std::string>("Parameter"));
112  }
113  LOG_VERBATIM("cmf::CovarianceMatrixMaker")
114  << "Running with " << fSystParameterNames.size() << " systematic parameters";
115  }
116 
117  //......................................................................
118  // Method to clear out the collections of data products after the
119  // writeResults method is called.
121  {
122  }
123 
124 
125  //......................................................................
127  {
128  // Create histograms for the covariance matrix, the nominal spectrum and
129  // the shifted values
130  // Each axis accounts for both detectors, both beam types and all desired selections
131  // Don't break out the interaction types or file types because they aren't known in
132  // the data
134 
135  fSystCovarianceMatrix = fTFS->make<TH2D>("CovarianceMatrix",
136  ";Energy Bin;Energy Bin",
137  numBins,
138  1,
139  numBins + 1,
140  numBins,
141  1,
142  numBins + 1);
143 
144  fNominalSpectrum = fTFS->make<TH1D>("NominalBin",
145  ";Energy Bin;Events",
146  numBins,
147  1,
148  numBins + 1);
149 
150  fShiftedSpectrum = fTFS->make<TH1D>("ShiftedBin",
151  ";Energy Bin;Events",
152  numBins,
153  1,
154  numBins + 1);
155 
156  fShiftedSpectrumByUniverse = fTFS->make<TH2D>("ShiftedBinByUniverse",
157  ";Energy Bin;Universe",
158  numBins,
159  1,
160  numBins + 1,
161  fIterations,
162  0,
163  fIterations);
164 
165  std::stringstream shifttitles;
166  shifttitles << "Seed: "
167  << fRand.GetSeed()
168  << ";Shifts;Iterations";
169  for(auto const& name : fSystParameterNames){
170  fShifts[name] = fTFS->make<TH1D>((name+"Shifts").c_str(),
171  shifttitles.str().c_str(),
172  500,
173  -5.,
174  5);
175 
176  }
177 
178  // histogram of the weights as a function of energy bin
179  std::stringstream shiftwgttitles;
180  shiftwgttitles << "Seed: "
181  << fRand.GetSeed()
182  << ";Energy Bin;Total Shift Weights";
183  fShiftWeights = fTFS->make<TH2D>("TotalShiftWeights",
184  shiftwgttitles.str().c_str(),
185  numBins,
186  1,
187  numBins + 1,
188  500,
189  -1.,
190  1);
191 
192 
193  long key = 0;
196  std::set<cmf::BeamType_t> beamTypes({cmf::kFHC,
197  cmf::kRHC});
198  std::set<novadaq::cnv::DetId> detectors({novadaq::cnv::kNEARDET,
200  std::set<cmf::FileType_t> fileTypes({cmf::kBeam,
201  cmf::kSwap,
206  std::set<cmf::SelectionType_t> selectionTypes({kNuESelectionLowPID,
213  kNCSelection});
214 
216  std::vector<double> bins;
217 
218  // loop over beams
219  for(auto const& bt : beamTypes){
220  md.epoch = (bt == cmf::kFHC) ? 1 : 4;
221 
222  // have to have the same bins in the ND and FD
223  for(auto const& det : detectors){
224  md.detector = det;
225 
226  title = (det == novadaq::cnv::kNEARDET) ? ";Energy (GeV);10^{3} Events/10^{20} POT" : ";Energy (GeV);Events/10^{20} POT";
227 
228  // loop over the selection types in each detector
229  for(auto const& sel : selectionTypes){
230 
231  // no peripheral nue selection in the ND
232  if(det == novadaq::cnv::kNEARDET &&
233  sel == cmf::kNuESelectionPeripheral) continue;
234 
235  md.selectionType = sel;
236 
237  auto const& highEdges = cmf::CovarianceBinUtility::Instance()->SelectionHighEdges(md);
238  bins.clear();
239  if(md.IsNuESelected() &&
240  md.selectionType != cmf::kNuESelectionPeripheral &&
241  md.detector == novadaq::cnv::kFARDET)
242  bins.insert(bins.begin(), 1.);
243  else
244  bins.insert(bins.begin(), 0.);
245 
246  for(auto const& itr : highEdges ) bins.push_back(itr.first);
247 
248  // loop over file types - the ND only cares about
249  // the nonswap files, FD has nonswap, fluxswap, tauswap,
250  // rocknonswap, rockfluxswap and cosmic background
251  for(auto fitr : fileTypes){
252 
253  // in the ND, we only care about the nonswap file type
254  if(det == novadaq::cnv::kNEARDET && fitr != cmf::kBeam) continue;
255 
256  key = md.DetectorBeamSelectionKey();
257  key += cmf::FileTypeKey(fitr);
258 
259  name = (cmf::cBeamType_Strings[bt] +
261  cmf::cFileTypeStrings[fitr] +
263 
264  fNominalEnergySpectra[key] = fTFS->make<TH1D>(name.c_str(),
265  title.c_str(),
266  bins.size() - 1,
267  &bins[0]);
268 
269  } // end loop over file types
270  } // end loop over selections
271  } // end loop over detectors
272  } // end loop over beam types
273 
274  LOG_DEBUG("CovarianceMatrixMaker")
275  << "histograms made";
276 
277  }
278 
279  //......................................................................
281  {
282  }
283 
284  // Function to fill the spectra for each meta data type
285  //......................................................................
287  TH1D* spectrumBin,
288  std::map<long, TH1D*> & spectraEnergy)
289  {
290  // first reset the spectra
291  spectrumBin->Reset();
292  for(auto hitr : spectraEnergy)
293  hitr.second->Reset();
294 
295  // declare some variables we will be using a lot in the for loops below
296  double weight;
297  double energy;
298  double normalization;
299  long key = 0;
300  int bin = 0;
301  int offset = 0;
302  double potForNorm = 1.e20; // spill summary POT values are in units of 10^12 POT,
303  // so this normalization would be to 10^20 POT
304 
305  // loop over each list and event in the list
306  for(auto const& itr : eventLists){
307  auto const& md = itr.first;
308 
309  normalization = potForNorm/itr.second.ListSpillSummary().goodPOT;
310 
311  // normalize the numu selected ND down to a reasonable number of events
312  // to have about the same peak as the highest fd bins. In the end we
313  // are just looking to see relative changes in each bin anyway.
314  if(md.detector == novadaq::cnv::kNEARDET){
315  if (md.IsNuMuSelected()) normalization *= 1.e-3;
316  else if(md.IsNuESelected() ) normalization *= 1.e-2;
317  else if(md.IsNCSelected() ) normalization *= 1.e-3;
318  }
319 
320  key = md.DetectorBeamSelectionKey();
321 
322  // figure out the offsets
324 
325  LOG_DEBUG("CovarianceMatrixMaker")
326  << md.ToString()
327  << " key: "
328  << key
329  << " bin offset: "
330  << offset;
331 
332  // add the file type to the key for the energy spectra
333  key += cmf::FileTypeKey(md.fileType);
334 
335  for(auto const& evt : itr.second){
336 
337  energy = evt->DataVals().val_at(cmf::kNu_RecoE, md);
339 
340  // if bin is negative then the energy is out of range for
341  // the selection type and we'll ignore it.
342  if(bin < 0) continue;
343 
345 
346  // now fill the spectrum
347  spectrumBin->Fill(bin + offset, weight * normalization);
348 
349  if(spectraEnergy.count(key) > 0)
350  spectraEnergy[key]->Fill(energy, weight * normalization);
351 
352  // fill the histogram of weights as a function of energy bin
353  fShiftWeights->Fill(bin + offset, weight * normalization);
354  } // end loop over events
355  } // end loop over lists of events
356  }
357 
358  // This function will pick a random value for the specified parameter
359  // from a guassian distribution. The random value will be within +/- 1 sigma
360  // of the nominal value
361  //......................................................................
363  {
364 
365  // we will pick a random value from a unit gaussian, ie gaussian with
366  // sigma = 1. That will be the fractional change on the sigma we want
367  // to apply to the new input point
368  double sigFrac = 0.;
369 
370  LOG_DEBUG("CovarianceMatrixMaker")
371  << "sigfrac = "
372  << sigFrac
373  << " original point "
374  << point;
375 
376  auto systematicPulls = point.SystematicPulls();
377 
378  for(auto const& name : fSystParameterNames) {
381 
382  // each systematic parameter should have a different random number
383  // change in sigma so that we don't artificially make them change
384  // in a correlated way if doing more than one parameter in a job
385  sigFrac = fRand.Gaus();
386 
387  // The calibration and light level systematic uncertainties are
388  // really only defined for a 1 sigma shift because their excursions
389  // are calculated from special data sets with the shift applied
390  // Also pull them from a uniform distribution as they are definitely
391  // not Gaussian.
392  if(name.find("lightmodel") != std::string::npos ||
393  name.find("ckv-proton" ) != std::string::npos ||
394  name.find("Calibration") != std::string::npos ){
395  sigFrac = fRand.Rndm();
396 
397  LOG_DEBUG("CovarianceMatrixMaker")
398  << "pulling from uniform distribution "
399  << sigFrac;
400  }
401 
402  fShifts[name]->Fill(sigFrac);
403 
404  LOG_DEBUG("CovarianceMatrixMaker")
405  << " parameter: "
406  << name
407  << " sigFrac "
408  << sigFrac;
409 
410  if(systematicPulls.count(tempPairFar) != 0){
411  systematicPulls[tempPairFar] = sigFrac;
412 
413  LOG_DEBUG("CovarianceMatrixMaker")
414  << "FD name: "
415  << name
416  << " "
417  << systematicPulls[tempPairFar];
418  }
419  if(systematicPulls.count(tempPairNear) != 0){
420  systematicPulls[tempPairNear] = sigFrac;
421 
422  LOG_DEBUG("CovarianceMatrixMaker")
423  << "ND name: "
424  << name
425  << " "
426  << systematicPulls[tempPairFar];
427  }
428  } // end loop over systematic parameters to change
429 
430  // Make the new input point
431  //
432  point = cmf::InputPoint(point.ParameterSpaceLocation(),
433  systematicPulls,
434  point.FixedParameters(),
435  point.FixedSystematics() );
436 
437  LOG_DEBUG("CovarianceMatrixMaker")
438  << " new point "
439  << point;
440 
442  }
443 
444  //......................................................................
446  {
447  // declare some variables we will be using a lot in the for loops below
448  int nbins = fNominalSpectrum->GetNbinsX();
449  double diffI;
450  double diffJ;
451  double iCenter;
452 
454 
455  // get the lists of events from the TTree and put them into a map
456  // of event lists
458  point);
459 
460 
461  // we only want the MC lists
462  cmf::EventListMap eventLists;
463  manipulator.Deserialize(eventLists, cmf::kMC);
464 
465  // The shifter and weighter was initialized in the EventListManipulator::Deserialize() method
466  // fill the nominal histograms first so we don't do this for every time we
467  // change the systematic parameters
468  this->FillSpectrum(eventLists,
471 
472  LOG_DEBUG("CovarianceMatrixMaker")
473  << "nominal spectrum filled";
474 
475  // loop over all the events and fill the matrices for many different
476  // values of the systematic parameters
477  for(int k = 0; k < fIterations; ++k){
478 
479  LOG_VERBATIM("CovarianceMatrixMaker")
480  << "iteration: "
481  << k;
482 
483  this->SetSystematicParameter(point);
484 
485  // this is a new shift for the parameter, so reset and refill the spectrum
487 
488  // now to fill the covariance matrices
489  // For ROOT histograms bin 0 is the underflow bin, so start at bin 1
490  // if either diffI or diffJ are 0, then we don't need to fill
491  // the bin with a 0 (ie the product) and can continue on
492  for(int i = 1; i < nbins + 1; ++i){
493  iCenter = fNominalSpectrum->GetBinCenter(i);
494  diffI = (fShiftedSpectrum->GetBinContent(i) - fNominalSpectrum->GetBinContent(i));
495  if(diffI != 0.) diffI /= fNominalSpectrum->GetBinContent(i);
496  else continue;
497 
498  for(int j = 1; j < nbins + 1; ++j){
499  diffJ = (fShiftedSpectrum->GetBinContent(j) - fNominalSpectrum->GetBinContent(j));
500  if (diffJ != 0.) diffJ /= fNominalSpectrum->GetBinContent(j);
501  else continue;
502  /*
503  if(i == 1)
504  LOG_VERBATIM("CovarianceMatrixMaker")
505  << "i: "
506  << i
507  << " shifted i: "
508  << fShiftedSpectrum->GetBinContent(i)
509  << " nominal i: "
510  << fNominalSpectrum->GetBinContent(i)
511  << " j: "
512  << j
513  << " shifted j: "
514  << fShiftedSpectrum->GetBinContent(j)
515  << " nominal j: "
516  << fNominalSpectrum->GetBinContent(j)
517  << " diffI: "
518  << diffI
519  << " diffJ: "
520  << diffJ
521  << " product: "
522  << diffI * diffJ;
523  */
524  fSystCovarianceMatrix->Fill(iCenter, fNominalSpectrum->GetBinCenter(j), diffI * diffJ);
525 
526  } // end loop over j bins of spectrum
527  } // end loop over i bins of spectrum
528 
529  // loop over each bin and fill fShiftedSpectrumByUniverse
530  for (int i_bin = 1; i_bin < fShiftedSpectrum->GetNbinsX(); ++i_bin){
531  fShiftedSpectrumByUniverse->SetBinContent(i_bin, k + 1, fShiftedSpectrum->GetBinContent(i_bin));
532  }
533  } // end loop over the iterations on the systematic parameter
534 
535  // if we're not running any systematic uncertainties,
536  // then assume we want stat only
537  // n.b. we're usually running the xsec weights
538  // so take that into account
539  if (fSystParameterNames.empty() || (fSystParameterNames.size() == 1 && fSystParameterNames[0] == "xseccvwgt")){
540  for (int i_bin = 0; i_bin < fNominalSpectrum->GetNbinsX(); ++i_bin){
541 
542  float binContent = fNominalSpectrum->GetBinContent(i_bin);
543 
544  if (binContent == 0){
545  fSystCovarianceMatrix->SetBinContent(i_bin, i_bin, std::numeric_limits<float>::min());
546  }
547  else {
548  fSystCovarianceMatrix->SetBinContent(
549  i_bin,i_bin,
550  binContent/(binContent * binContent));
551  }
552  }
553  }
554 
555  // protect against very low and very high covariances
556  // these usually happen because nentries in these bins for nominal
557  // and shifted spectra are identially 0
558  for (int i = 1; i < nbins + 1; ++i){
559  for (int j = 1; j < fSystCovarianceMatrix->GetNbinsY()+1; ++j){
560  double binContent = fSystCovarianceMatrix->GetBinContent(i, j);
561  if (std::abs(binContent) < std::pow(10,-14) ||
562  std::abs(binContent) > std::pow(10,14) ||
563  std::isnan(binContent) ){
564  // LOG_VERBATIM("CovarianceMatrixMaker")
565  // << "bin " << i
566  // << ", covariance: " << binContent
567  // << "\n -- nominal i: " << fNominalSpectrum->GetBinContent(i)
568  // << "\n -- shifted i: " << fShiftedSpectrum->GetBinContent(i)
569  // << "\n -- nominal minus shifted i: " << fNominalSpectrum->GetBinContent(i) - fShiftedSpectrum->GetBinContent(i)
570  // << "\n -- nominal j: " << fNominalSpectrum->GetBinContent(j)
571  // << "\n -- shifted j: " << fShiftedSpectrum->GetBinContent(j)
572  // << "\n -- nominal minus shifted j: " << fNominalSpectrum->GetBinContent(j) - fShiftedSpectrum->GetBinContent(j)
573  // << "\n -- diffI*diffJ" << (fNominalSpectrum->GetBinContent(i)-fShiftedSpectrum->GetBinContent(i))*(fNominalSpectrum->GetBinContent(j)-fShiftedSpectrum->GetBinContent(j));
574 
575  fSystCovarianceMatrix->SetBinContent(i, j, 0);
576  fSystCovarianceMatrix->SetBinContent(j, i, 0);
577  } // end if the bin content is acceptable
578  } // end loop over j bins
579  } // end loop over i bins
580  }
581 
583 
584 } // end cmf namespace
static std::string GetName(int id)
CovarianceMatrixMaker & operator=(CovarianceMatrixMaker const &)=delete
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
const XML_Char * name
Definition: expat.h:151
int fIterations
number of iterations for picking sigmas
void writeResults(art::Results &r) override
const std::string cSelectionType_Strings[12]
Definition: Constants.h:115
static const unsigned char kNu_RecoE
Definition: VarVals.h:35
std::map< long, TH1D * > fNominalEnergySpectra
Nominal energy spectrum for each detector/file/beam/selection.
void SetCurrentVals(cmf::InputPoint const &curPoint)
int KeyToOffset(long const &key)
const Var weight
ParameterMap const SystematicPulls(novadaq::cnv::DetId const &det) const
Definition: InputPoint.cxx:352
constexpr T pow(T x)
Definition: pow.h:75
void SetSystematicParameter(cmf::InputPoint &point)
TH1D * fNominalSpectrum
Nominal spectrum for each detector/beam/selection.
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
TH2D * fShiftWeights
Weights from shifts selected for each parameter adjusted.
void reconfigure(fhicl::ParameterSet const &pset)
void FillSpectrum(cmf::EventListMap const &eventLists, TH1D *spectrumBin, std::map< long, TH1D * > &spectraEnergy)
float abs(float number)
Definition: d0nt_math.hpp:39
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
void Deserialize(cmf::EventListMap &eventLists, cmf::DataMC_t dataMC=cmf::kBoth)
TH1D * fShiftedSpectrum
Shifted spectrum for each detector/beam/selection.
const int nbins
Definition: cellShifts.C:15
const std::string cFileTypeStrings[10]
Definition: Constants.h:380
TH2D * fShiftedSpectrumByUniverse
Shifted spectrum as a function of universe.
Far Detector at Ash River, MN.
std::pair< std::string, novadaq::cnv::DetId > ParameterDet
Definition: Constants.h:38
fhicl::ParameterSet fManipulatorParameters
the fake point 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
double energy
Definition: plottest35.C:25
fhicl::ParameterSet fInputPointParameters
input parameters
Near Detector in the NuMI cavern.
const double j
Definition: BetheBloch.cxx:29
const Binning bins
static ShifterAndWeighter * Instance()
std::vector< std::string > get_names() const
float bin[41]
Definition: plottest35.C:14
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
void readResults(art::Results const &r) override
Module to combine a set of results into a single file.
Definition: Event.cxx:24
CovarianceMatrixMaker(fhicl::ParameterSet const &pset)
std::map< std::string, TH1D * > fShifts
Shifts selected for each parameter adjusted.
std::map< cmf::MetaData, cmf::EventList > EventListMap
Definition: Event.h:148
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
art::ServiceHandle< art::TFileService > fTFS
TFileService.
std::map< double, int > const & SelectionHighEdges(cmf::MetaData const &md)
static long FileTypeKey(cmf::FileType_t const &ft)
Definition: StaticFuncs.h:18
TRandom3 r(0)
std::map< long, TH1D * > fShiftedEnergySpectra
Shifted energy spectrum for each detector/file/beam/selection.
std::vector< std::string > fSystParameterNames
names of parameters we are floating from the input point
int EnergyToHistogramBin(double const &energy, cmf::MetaData const &md)
#define LOG_VERBATIM(category)
Float_t e
Definition: plot.C:35
const std::string cBeamType_Strings[4]
Definition: Constants.h:75
double Weight(cmf::EventPtr const &curEvent, cmf::MetaData const &md, cmf::WeightType_t wgtType=kAllWgt)
TRandom3 fRand
RNG for picking the variation in parameters.
static CovarianceBinUtility * Instance()