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