CovarianceMatrix.cxx
Go to the documentation of this file.
1 
4 #include "CAFAna/Core/Progress.h"
8 
9 #include "TH1D.h"
10 #include "TRandom3.h"
11 #include "TDirectory.h"
12 #include "TString.h"
13 #include "TMatrixD.h"
14 
15 #include <iostream>
16 #include <algorithm>
17 #include <vector>
18 #include <tuple>
19 #include <sstream>
20 
21 using std::back_inserter;
22 using std::begin;
23 using std::cout;
24 using std::end;
25 using std::endl;
26 using std::get;
27 using std::make_pair;
28 using std::make_tuple;
29 using std::make_unique;
30 using std::ostringstream;
31 using std::pair;
32 using std::runtime_error;
33 using std::transform;
34 using std::unique_ptr;
35 using std::vector;
36 
37 using Eigen::MatrixXd;
38 using Eigen::VectorXd;
39 using Eigen::Map;
40 
41 namespace ana
42 {
43  namespace covmx
44  {
45  //-------------------------------------------------------------------------
46  vector<Component> GetComponents(Sample sample)
47  {
48  vector<Component> ret;
49 
50  // Define vectors of flavour and sign
51  vector<Flavors::Flavors_t> flavours = { Flavors::kNuEToNuE,
53  if (sample.detector == kFarDet) {
54  vector<Flavors::Flavors_t> fdFlavours = { Flavors::kNuEToNuMu,
56  flavours.insert(end(flavours), begin(fdFlavours), end(fdFlavours));
57  }
58  vector<Sign::Sign_t> signs = { Sign::kNu, Sign::kAntiNu };
59 
60  // Add each combination of flavour and sign
61  for (auto & flavour: flavours) {
62  for (auto & sign: signs) {
63  ret.push_back(make_tuple(flavour, Current::kCC, sign));
64  } // for flavour
65  } // for sign
66 
67  // Add the neutral currents
68  ret.push_back(make_tuple(Flavors::kAll, Current::kNC, Sign::kBoth));
69  return ret;
70 
71  } // function GetComponents
72 
73  //-------------------------------------------------------------------------
74  CovarianceMatrix::CovarianceMatrix(vector<Sample> samples,
75  vector<const ISyst*> systs, int nUniverses)
76  : fNUniverses(nUniverses)
77  {
78  BuildFullCovMx(samples, systs);
79 
80  } // function CovarianceMatrix::CovarianceMatrix
81 
82  //-------------------------------------------------------------------------
83  CovarianceMatrix::CovarianceMatrix(MatrixXd fullCovMx,
84  MatrixXd covMxAbsolute, vector<pair<Binning, unsigned int>> bins)
85  : fBins(bins), fFullCovMx(fullCovMx), fCovMxAbsolute(covMxAbsolute),
86  fNUniverses(1e5)
87  {} // CovarianceMatrix constructor
88 
89  //-------------------------------------------------------------------------
90  CovarianceMatrix::CovarianceMatrix(MatrixXd fullCovMx,
91  vector<Sample> samples)
92  : fFullCovMx(fullCovMx), fNUniverses(1e5)
93  {
94  SetBinning(samples);
95  unsigned int nBins = GetBinning().NBins();
96  fCovMxAbsolute.resize(nBins, nBins);
97  fCovMxAbsolute.setZero();
98 
99  } // CovarianceMatrix constructor
100 
101  //-------------------------------------------------------------------------
102  /// \brief Constructor that makes an empty matrix from given bins
103  CovarianceMatrix::CovarianceMatrix(vector<Sample> samples)
104  : fNUniverses(1e5)
105  {
106  SetBinning(samples);
107  // Create empty covariance matrices
108  unsigned int nBinsFull = GetFullBinning().NBins();
109  fFullCovMx.resize(nBinsFull, nBinsFull);
110  fFullCovMx.setZero();
111  unsigned int nBins = GetBinning().NBins();
112  fCovMxAbsolute.resize(nBins, nBins);
113  fCovMxAbsolute.setZero();
114 
115  } // CovarianceMatrix constructor
116 
117  //----------------------------------------------------------------------
118  void CovarianceMatrix::BuildFullCovMx(vector<Sample> samples,
119  vector<const ISyst*> systs)
120  {
122 
123  // Set matrix binning
124  if (!fBins.empty()) {
125  throw runtime_error("Cannot reset matrix binning!");
126  }
127  for (Sample s: samples) {
128  Binning bins = s.GetPrediction()->Predict(calc).GetBinnings()[0];
129  unsigned int nComps = GetComponents(s).size();
130  fBins.push_back(make_pair(bins, nComps));
131  }
132 
133  // Instantiate empty matrices
134  vector<Binning> bins = GetBinnings();
135  unsigned int nBinsFull = GetFullBinning().NBins();
136  fFullCovMx.resize(nBinsFull, nBinsFull);
137  fFullCovMx.setZero();
138  unsigned int nBins = GetBinning().NBins();
139  fCovMxAbsolute.resize(nBins, nBins);
140  fCovMxAbsolute.setZero();
141 
142  vector<double> nominalPrediction(nBinsFull);
143  vector<double> fractionalShift(nBinsFull);
144 
145  size_t i = 0;
146  // Loop over samples
147  for (size_t iPred = 0; iPred < samples.size(); ++iPred) {
148  // Loop over components
149  for (auto& component: GetComponents(samples[iPred])) {
150  // Get spectrum and use it to set nominal prediction
151  Spectrum spec = samples[iPred].GetPrediction()->PredictComponent(calc,
152  get<0>(component), get<1>(component), get<2>(component));
153  TH1D* hist = spec.ToTH1(samples[iPred].GetPOT());
154  for (size_t iBin = 1; iBin <= (size_t)bins[iPred].NBins(); ++iBin) {
155  nominalPrediction[i] = hist->GetBinContent(iBin);
156  ++i;
157  }
158  delete hist;
159  }
160  }
161 
162  TRandom3 rnd(0);
163 
164  Progress progress("Constructing the full covariance matrix");
165  SystShifts shifts;
166 
167  // Randomly throw many universes with systematic shifts
168  for (Long64_t iThrow = 0; iThrow < fNUniverses; ++iThrow) {
169  progress.SetProgress(double(iThrow)/fNUniverses);
170 
171  // Set random Gaussian shifts for all systematics
172  for (const ISyst* syst : systs) {
173  double shift = rnd.Gaus();
174  const KeySyst* tmp = dynamic_cast<const KeySyst*>(syst); // Cast to key syst for one-sidedness
175  if (tmp && tmp->IsOneSided()) shift = fabs(shift);
176  shifts.SetShift(syst, shift);
177  }
178 
179  // Loop over components
180  i = 0;
181  for (size_t iPred = 0; iPred < samples.size(); ++iPred) {
182  SystShifts sampleShifts = samples[iPred].GetSystShifts(shifts);
183  auto comps = covmx::GetComponents(samples[iPred]);
184  for (auto& component : comps) {
185 
186  Spectrum shiftedSpec = samples[iPred].GetPrediction()->PredictComponentSyst(calc,
187  sampleShifts, get<0>(component), get<1>(component), get<2>(component));
188  auto shiftedHist = shiftedSpec.ToTH1(samples[iPred].GetPOT());
189 
190  // Get fractional shift in each energy bin
191  for (size_t iBin = 1; iBin <= (size_t)bins[iPred].NBins(); ++iBin) {
192 
193  // Get systematic-shifted prediction
194  double shiftedPrediction = shiftedHist->GetBinContent(iBin);
195  if (std::isnan(shiftedPrediction))
196  throw runtime_error("NaN value found in shifted prediction.");
197 
198  // Set fractional shift for this bin
199  if (nominalPrediction[i] > 0) fractionalShift[i] = ((shiftedPrediction+0.01)/(nominalPrediction[i]+0.01)) - 1;
200  else fractionalShift[i] = 0;
201  ++i;
202  }
203  delete shiftedHist;
204  } // for prediction
205  } // for component
206 
207  // Add fractional shifts to covariance matrix
208  for (size_t iBin = 0; iBin < nBinsFull; ++iBin) {
209  for (size_t jBin = 0; jBin < nBinsFull; ++jBin) {
210  fFullCovMx(iBin, jBin) += fractionalShift[iBin]*fractionalShift[jBin]*(1.0/(fNUniverses));
211  }
212  }
213  }
214 
215  progress.Done();
216 
217  delete calc;
218 
219  } // function CovarianceMatrix::BuildFullCovMx
220 
221  //-------------------------------------------------------------------------
222  MatrixXd CovarianceMatrix::Predict(vector<Sample> samples,
224  const SystShifts& shifts)
225  {
226  DontAddDirectory guard;
227 
228  fCovMxAbsolute.setZero();
229 
230  // Prepare oscillated prediction vector
231  size_t nBinsFull = GetFullBinning().NBins();
232  vector<double> oscPrediction(nBinsFull);
233 
234  size_t i = 0;
235  for (size_t iPred = 0; iPred < samples.size(); ++iPred) {
236  auto comps = covmx::GetComponents(samples[iPred]);
237  for (auto& component : comps) {
238  SystShifts sampleShifts = samples[iPred].GetSystShifts(shifts);
239  Spectrum spec = samples[iPred].GetPrediction()->PredictComponentSyst(calc,
240  sampleShifts, get<0>(component), get<1>(component), get<2>(component));
241  TH1D* hist = spec.ToTH1(samples[iPred].GetPOT());
242  for (int iBin = 1; iBin <= spec.GetBinnings()[0].NBins(); ++iBin) {
243  oscPrediction[i] = hist->GetBinContent(iBin);
244  ++i;
245  }
246  delete hist;
247  } // for component
248  } // for sample
249 
250  // Loop over prediction, component & bin in i
251  size_t iOffset(0), iFullOffset(0);
252  vector<unsigned int> nComps = GetNComponents();
253  for (size_t iPred = 0; iPred < samples.size(); ++iPred) {
254  size_t iBins = fBins[iPred].first.NBins();
255  for (size_t iComp = 0; iComp < nComps[iPred]; ++iComp) {
256  // Loop over prediction, component & bin in j
257  for (size_t iBin = 0; iBin < iBins; ++iBin) {
258  size_t jOffset(0), jFullOffset(0);
259  for (size_t jPred = 0; jPred < samples.size(); ++jPred) {
260  size_t jBins = fBins[jPred].first.NBins();
261  for (size_t jComp = 0; jComp < nComps[jPred]; ++jComp) {
262  for (size_t jBin = 0; jBin < jBins; ++jBin) {
263  size_t i(iOffset+iBin), iFull(iFullOffset+iBin),
264  j(jOffset+jBin), jFull(jFullOffset+jBin);
265  fCovMxAbsolute(i, j) += fFullCovMx(iFull, jFull) * oscPrediction[iFull] * oscPrediction[jFull];
266  } // for jBin
267  jFullOffset += jBins;
268  } // for jComp
269  jOffset += jBins;
270  } // for jPred
271  } // for iBin
272  iFullOffset += iBins;
273  } // for iComp
274  iOffset += iBins;
275  }// for iPred
276 
277  return fCovMxAbsolute;
278 
279  } // function CovarianceMatrix::Predict
280 
281  //----------------------------------------------------------------------
282  MatrixXd CovarianceMatrix::Predict(vector<double> vBetaMu) {
283  DontAddDirectory guard;
284  MatrixXd ret(vBetaMu.size(), vBetaMu.size());
285  ret.setZero();
286  size_t iOffset(0), iFullOffset(0);
287  vector<unsigned int> nComps = GetNComponents();
288  for (size_t iPred = 0; iPred < fBins.size(); ++iPred) {
289  size_t iBins = fBins[iPred].first.NBins();
290  for (size_t iComp = 0; iComp < nComps[iPred]; ++iComp) {
291  // Loop over prediction, component & bin in j
292  for (size_t iBin = 0; iBin < iBins; ++iBin) {
293  size_t jOffset(0), jFullOffset(0);
294  for (size_t jPred = 0; jPred < fBins.size(); ++jPred) {
295  size_t jBins = fBins[jPred].first.NBins();
296  for (size_t jComp = 0; jComp < nComps[jPred]; ++jComp) {
297  for (size_t jBin = 0; jBin < jBins; ++jBin) {
298  size_t i(iOffset+iBin), iFull(iFullOffset+iBin),
299  j(jOffset+jBin), jFull(jFullOffset+jBin);
300  ret(i, j) += fFullCovMx(iFull, jFull) * vBetaMu[iFull] * vBetaMu[jFull];
301  } // for jBin
302  jFullOffset += jBins;
303  } // for jComp
304  jOffset += jBins;
305  } // for jPred
306  } // for iBin
307  iFullOffset += iBins;
308  } // for iComp
309  iOffset += iBins;
310  }// for iPred
311  return ret;
312  } // function CovarianceMatrix::Predict
313 
314  //----------------------------------------------------------------------
316  {
317  unsigned int nBins = GetBinning().NBins();
318  unsigned int nBinsComp = gen->GetBinning().NBins();
319  unsigned int nBinsFull = GetFullBinning().NBins();
320  unsigned int nBinsFullComp = gen->GetFullBinning().NBins();
321  // Sanity check
322  if (nBins != nBinsComp || nBinsFull != nBinsFullComp) {
323  throw runtime_error("Cannot add matrices; bins do not match.");
324  }
325 
326  fFullCovMx += gen->GetFullCovMx();
328 
329  } // function CovarianceMatrix::AddMatrix
330 
331  //----------------------------------------------------------------------
332  MatrixXd CovarianceMatrix::ForcePosDef(MatrixXd mx, double epsilon)
333  {
334  MatrixXd ret(mx);
335  // Add on epsilon
336  for (size_t i = 0; i < (size_t)ret.rows(); ++i) {
337  ret(i,i) += epsilon;
338  }
339  return ret;
340 
341  } // function CovarianceMatrix::ForcePosDef
342 
343  //---------------------------------------------------------------------------
345  {
346  vector<Binning> ret;
347  transform(begin(fBins), end(fBins), back_inserter(ret),
348  [](auto const& pair){ return pair.first; });
349  return ret;
350 
351  } // function CovarianceMatrix::GetBinnings
352 
353  //---------------------------------------------------------------------------
355  {
356  // Concatenate together all the individual binning schemes
357  vector<double> mergedBins = { 0. };
358  vector<Binning> bins = GetBinnings();
359  for (Binning b : bins) {
360  for (size_t i = 1; i < b.Edges().size(); ++i) {
361  double width = b.Edges()[i] - b.Edges()[i-1];
362  mergedBins.push_back(mergedBins.back() + width);
363  } // for bin
364  } // for sample
365  return Binning::Custom(mergedBins);
366 
367  } // function CovarianceMatrix::GetBinning
368 
369  //---------------------------------------------------------------------------
371  {
372 
373  // Concatenate together all the binning schemes component-wise
374  vector<double> mergedBins = { 0. };
375  auto bins = fBins;
376  for (auto const& b : bins) {
377  vector<double> edges = b.first.Edges();
378  for (size_t iComp = 0; iComp < b.second; ++iComp) {
379  for (size_t iBin = 1; iBin < edges.size(); ++iBin) {
380  double width = edges[iBin] - edges[iBin-1];
381  mergedBins.push_back(mergedBins.back() + width);
382  } // for iBin
383  } // for iComp
384  } // for sample
385  return Binning::Custom(mergedBins);
386 
387  } // function CovarianceMatrix::GetFullBinning
388 
389  //---------------------------------------------------------------------------
390  vector<unsigned int> CovarianceMatrix::GetNBins()
391  {
392  vector<unsigned int> ret;
393  transform(begin(fBins), end(fBins), back_inserter(ret),
394  [](auto const& pair){ return pair.first.NBins(); });
395  return ret;
396 
397  } // function CovarianceMatrix::GetNBins
398 
399  //---------------------------------------------------------------------------
400  vector<unsigned int> CovarianceMatrix::GetNComponents()
401  {
402  vector<unsigned int> ret;
403  transform(begin(fBins), end(fBins), back_inserter(ret),
404  [](auto const& pair){ return pair.second; });
405  return ret;
406 
407  } // function CovarianceMatrix::GetNComponents
408 
409  //---------------------------------------------------------------------------
410  MatrixXd CovarianceMatrix::GetCovMxRelative(vector<Sample>& samples,
412  {
413  Predict(samples, calc);
414  MatrixXd ret(fCovMxAbsolute);
415  int c1 = 0;
416  for (Sample& s1 : samples) {
417  TH1D* h1 = s1.GetPrediction()->Predict(calc).ToTH1(s1.GetPOT());
418  for (int b1 = 1; b1 <= h1->GetNbinsX(); ++b1) {
419  double iVal = h1->GetBinContent(b1);
420  int c2 = 0;
421  for (Sample& s2 : samples) {
422  TH1D* h2 = s2.GetPrediction()->Predict(calc).ToTH1(s2.GetPOT());
423  for (int b2 = 1; b2 <= h2->GetNbinsX(); ++b2) {
424  double jVal = h2->GetBinContent(b2);
425  ret(c1, c2) /= iVal * jVal;
426  ++c2;
427  } // for bin b2
428  delete h2;
429  } // for sample s2
430  ++c1;
431  } // for bin b1
432  delete h1;
433  } // for sample s1
434 
435  return ret;
436 
437  } // function CovarianceMatrix::GetCovMxRelative
438 
439  //----------------------------------------------------------------------
441  {
442  // Populate the histogram with the elements of the full correlation
443  // matrix
444  MatrixXd corrMx = fFullCovMx;
445  corrMx.setZero();
446  for (int i = 0; i < corrMx.rows(); ++i) {
447  for (int j = 0; j < corrMx.cols(); ++j) {
448  double cov = fFullCovMx(i, j);
449  double sigma_1 = fFullCovMx(i, i);
450  double sigma_2 = fFullCovMx(j, j);
451  if (sigma_1 > 0 && sigma_2 > 0) {
452  corrMx(i, j) = cov/(sqrt(sigma_1) * sqrt(sigma_2));
453  } else {
454  corrMx(i, j) = 0;
455  }
456  }
457  }
458 
459  return corrMx;
460 
461  } // function CovarianceMatrix::GetFullCorrMx
462 
463  //----------------------------------------------------------------------
465  {
466  DontAddDirectory guard;
467  // Get full binning and create histogram
469  size_t nBins = bins.NBins();
470  vector<double> edges = bins.Edges();
471  TH2D* hist = new TH2D(UniqueName().c_str(),
472  ";Reco E bin (components);Reco E bin (components)",
473  nBins, &edges[0], nBins, &edges[0]);
474  // Populate histogram from matrix
475  for (size_t i = 0; i < nBins; ++i) {
476  for (size_t j = 0; j < nBins; ++j) {
477  hist->SetBinContent(i+1, j+1, fFullCovMx(i,j));
478  }
479  }
480  return hist;
481 
482  } // function CovarianceMatrix::GetFullCovMxTH2
483 
484  //----------------------------------------------------------------------
486  {
487  DontAddDirectory guard;
488  // Get full binning and create histogram
489  Binning bins = GetBinning();
490  size_t nBins = bins.NBins();
491  vector<double> edges = bins.Edges();
492  TH2D* hist = new TH2D(UniqueName().c_str(),
493  ";Reco E bin;Reco E bin",
494  nBins, &edges[0], nBins, &edges[0]);
495  // Populate histogram from matrix
496  for (size_t i = 0; i < nBins; ++i) {
497  for (size_t j = 0; j < nBins; ++j) {
498  hist->SetBinContent(i+1, j+1, fCovMxAbsolute(i,j));
499  }
500  }
501  return hist;
502 
503  } // function CovarianceMatrix::GetOscCovMxAbsTH2
504 
505  //----------------------------------------------------------------------
506  TH2D* CovarianceMatrix::GetCovMxRelativeTH2(vector<Sample>& samples,
508  {
509  DontAddDirectory guard;
510  // Get full binning and create histogram
511  Binning bins = GetBinning();
512  size_t nBins = bins.NBins();
513  vector<double> edges = bins.Edges();
514  TH2D* hist = new TH2D(UniqueName().c_str(),
515  ";Reco E bin;Reco E bin",
516  nBins, &edges[0], nBins, &edges[0]);
517  // Populate histogram from matrix
518  MatrixXd mx = GetCovMxRelative(samples, calc);
519  for (size_t i = 0; i < nBins; ++i) {
520  for (size_t j = 0; j < nBins; ++j) {
521  hist->SetBinContent(i+1, j+1, mx(i,j));
522  }
523  }
524  return hist;
525 
526  } // function CovarianceMatrix::GetOscCovMxRelTH2
527 
528  //----------------------------------------------------------------------
530  {
531  DontAddDirectory guard;
532 
533  Binning mxBins = GetBinning();
534  TH2D* hist = new TH2D(UniqueName().c_str(),
535  ";Reco E bin;Reco E bin",
536  mxBins.NBins(), &mxBins.Edges()[0],
537  mxBins.NBins(), &mxBins.Edges()[0]);
538 
539  for (int i = 0; i < mxBins.NBins(); ++i) {
540  for (int j = 0; j < mxBins.NBins(); ++j) {
541  double cov = fCovMxAbsolute(i,j);
542  double sigma_1 = fCovMxAbsolute(i,i);
543  double sigma_2 = fCovMxAbsolute(j,j);
544  if(sigma_1 > 0 && sigma_2 > 0)
545  hist->SetBinContent(i+1, j+1, cov/(sqrt(sigma_1)*sqrt(sigma_2)));
546  else hist->SetBinContent(i+1, j+1, 0);
547  }
548  }
549  return hist;
550 
551  } // function CovarianceMatrix::GetCorrMxTH2
552 
553  //----------------------------------------------------------------------
555  {
556  DontAddDirectory guard;
557 
558  // Get the binning of the full matrix, and then instantiate an empty
559  // 2D histogram
560  Binning mxBins = GetFullBinning();
561  TH2D* hist = new TH2D(UniqueName().c_str(),
562  ";Reco E bin (components);Reco E bin (components)",
563  mxBins.NBins(), &mxBins.Edges()[0],
564  mxBins.NBins(), &mxBins.Edges()[0]);
565 
566  // Populate the histogram with the elements of the full correlation
567  // matrix
568  MatrixXd corrMx = GetFullCorrMx();
569  for (int i = 0; i < mxBins.NBins(); ++i) {
570  for (int j = 0; j < mxBins.NBins(); ++j) {
571  hist->SetBinContent(i+1, j+1, corrMx(i, j));
572  }
573  }
574 
575  return hist;
576 
577  } // function CovarianceMatrix::GetFullCorrMxTH2
578 
579  //-------------------------------------------------------------------------
580  void CovarianceMatrix::SetBinning(vector<Sample> samples)
581  {
582  // Binning only gets set once, don't let it be reset
583  if (!fBins.empty()) {
584  throw runtime_error("Cannot reset matrix binning!");
585  }
586 
587  // Loop over samples and fill binning scheme & number of beam
588  // components for each one
589  for (Sample& s : samples) {
590  fBins.push_back(make_pair(s.GetBinning(),
591  GetComponents(s).size()));
592  }
593 
594  } // function CovarianceMatrix::SetBinning
595 
596  //----------------------------------------------------------------------
597  void CovarianceMatrix::SaveTo(TDirectory* dir) const
598  {
599  TDirectory* tmp = gDirectory;
600 
601  dir->cd();
602 
603  TObjString("CovarianceMatrix").Write("type");
604 
605  TMatrixD fullCovMx(fFullCovMx.rows(), fFullCovMx.cols(),
606  fFullCovMx.data());
607  fullCovMx.Write("fullcovmx");
608  TMatrixD covMxAbsolute(fCovMxAbsolute.rows(), fCovMxAbsolute.cols());
609  for (int i = 0; i < fCovMxAbsolute.rows(); ++i) {
610  for (int j = 0; j < fCovMxAbsolute.cols(); ++j) {
611  covMxAbsolute(i, j) = fCovMxAbsolute(i, j);
612  }
613  }
614  covMxAbsolute.Write("covmxabsolute");
615 
616  // Save binning information
617  vector<vector<double>> edges(fBins.size());
618  vector<unsigned int> nComps(fBins.size(), 0);
619  for (size_t it = 0; it < fBins.size(); ++it) {
620  edges[it] = fBins[it].first.Edges();
621  nComps[it] = fBins[it].second;
622  }
623  dir->WriteObject(&edges, "binedges");
624  dir->WriteObject(&nComps, "ncomps");
625 
626  tmp->cd();
627 
628  } // function CovarianceMatrix::SaveTo
629 
630  //----------------------------------------------------------------------
631  unique_ptr<CovarianceMatrix> CovarianceMatrix::LoadFrom(TDirectory* dir)
632  {
633  TObjString* tag = (TObjString*)dir->Get("type");
634  assert(tag);
635  assert(tag->GetString() == "CovarianceMatrix");
636 
637  // Load matrices
638  TMatrixD* tmp = (TMatrixD*)dir->Get("fullcovmx");
639  assert(tmp);
640  MatrixXd fullCovMx = Map<MatrixXd>(tmp->GetMatrixArray(),
641  tmp->GetNrows(), tmp->GetNcols());
642  delete tmp;
643  tmp = (TMatrixD*)dir->Get("covmxabsolute");
644  assert(tmp);
645  MatrixXd covMxAbsolute = Map<MatrixXd>(tmp->GetMatrixArray(),
646  tmp->GetNrows(), tmp->GetNcols());
647  delete tmp;
648 
649  // Load binning information
650  vector<vector<double>>* edges;
651  vector<unsigned int>* nComps;
652  vector<pair<Binning, unsigned int>> bins;
653  dir->GetObject("binedges", edges);
654  dir->GetObject("ncomps", nComps);
655  for (size_t it = 0; it < edges->size(); ++it) {
656  bins.push_back(make_pair(Binning::Custom(edges->at(it)),
657  nComps->at(it)));
658  }
659 
660  unique_ptr<CovarianceMatrix> ret(make_unique<CovarianceMatrix>
661  (fullCovMx, covMxAbsolute, bins));
662 
663  return ret;
664 
665  } // function CovarianceMatrix::LoadFrom
666 
667  } // namespace covmx
668 } // namespace ana
669 
TH2D * GetCovMxRelativeTH2(std::vector< Sample > &samples, osc::IOscCalcAdjustable *calc)
const std::vector< Binning > & GetBinnings() const
Definition: Spectrum.h:256
int nBins
Definition: plotROC.py:16
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
vector< Component > GetComponents(Sample sample)
set< int >::iterator it
Antineutrinos-only.
Definition: IPrediction.h:50
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
(&#39; appearance&#39;)
Definition: IPrediction.h:18
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
Eigen::MatrixXd GetFullCovMx()
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
(&#39;beam &#39;)
Definition: IPrediction.h:15
T sqrt(T number)
Definition: d0nt_math.hpp:156
Eigen::MatrixXd fCovMxAbsolute
Oscillated absolute covariance matrix.
Eigen::MatrixXd GetCovMxAbsolute()
void AddMatrix(CovarianceMatrix *gen)
bool progress
Insert implicit nodes #####.
Version of OscCalcSterile that always returns probability of 1.
Float_t tmp
Definition: plot.C:36
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
osc::OscCalcDumb calc
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
c2
Definition: demo5.py:33
const XML_Char * s
Definition: expat.h:262
Charged-current interactions.
Definition: IPrediction.h:39
bool IsOneSided() const
Definition: CovMxSysts.h:41
Eigen::MatrixXd fFullCovMx
Full covariance matrix.
std::vector< Binning > GetBinnings()
Eigen::MatrixXd GetCovMxRelative(std::vector< Sample > &samples, osc::IOscCalcAdjustable *calc)
virtual void SaveTo(TDirectory *dir) const
Detector detector
Definition: Sample.h:100
static Eigen::MatrixXd ForcePosDef(Eigen::MatrixXd mx, double epsilon=0.1)
const double j
Definition: BetheBloch.cxx:29
std::vector< std::string > comps
static Binning Custom(const std::vector< double > &edges)
Definition: Binning.cxx:161
Neutrinos-only.
Definition: IPrediction.h:49
const std::vector< double > & Edges() const
Definition: Binning.h:34
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
TH1F * h2
Definition: plot.C:45
std::vector< unsigned int > GetNComponents()
const Binning bins
Definition: NumuCC_CPiBin.h:8
TH1F * h1
std::vector< std::pair< Binning, unsigned int > > fBins
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
TDirectory * dir
Definition: macro.C:5
static std::unique_ptr< CovarianceMatrix > LoadFrom(TDirectory *dir)
int NBins() const
Definition: Binning.h:29
std::vector< unsigned int > GetNBins()
double epsilon
const hit & b
Definition: hits.cxx:21
Neutral-current interactions.
Definition: IPrediction.h:40
void SetBinning(std::vector< Sample > samples)
assert(nhit_max >=nhit_nbins)
A simple ascii-art progress bar.
Definition: Progress.h:9
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
c1
Definition: demo5.py:24
double GetPOT(bool isFHC=true)
void BuildFullCovMx(std::vector< Sample > samples, std::vector< const ISyst * > systs)
Eigen::MatrixXd Predict(std::vector< Sample > samples, osc::IOscCalcAdjustable *calc, const SystShifts &systs=SystShifts::Nominal())
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
void Done()
Call this when action is completed.
Definition: Progress.cxx:90
unsigned int fNUniverses
Number of universes to simulate.
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
def sign(x)
Definition: canMan.py:197
A class for generating a covariance matrices as a function of oscillation parameters and systematics ...