Public Member Functions | Static Public Member Functions | Protected Member Functions | Protected Attributes | Friends | List of all members
ana::covmx::CovarianceMatrix Class Reference

A class for generating a covariance matrices as a function of oscillation parameters and systematics as nuisance parameters. More...

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-04-11/CAFAna/Prediction/CovarianceMatrix.h"

Public Member Functions

virtual ~CovarianceMatrix ()
 
 CovarianceMatrix (std::vector< Sample > samples, std::vector< const ISyst * > systs, int nUniverses=1e5)
 
 CovarianceMatrix (Eigen::MatrixXd fullCovMx, Eigen::MatrixXd covMxAbsolute, std::vector< std::pair< Binning, unsigned int >> bins)
 
 CovarianceMatrix (Eigen::MatrixXd fullCovMx, std::vector< Sample > samples)
 
 CovarianceMatrix (std::vector< Sample > samples)
 Constructor that makes an empty matrix from given bins. More...
 
Eigen::MatrixXd Predict (std::vector< Sample > samples, osc::IOscCalcAdjustable *calc, const SystShifts &systs=SystShifts::Nominal())
 
Eigen::MatrixXd Predict (std::vector< double > vBetaMu)
 
void AddMatrix (CovarianceMatrix *gen)
 
std::vector< SpectrumGetUniverse (osc::OscCalcSterile *calc, std::vector< covmx::Sample > samples)
 
Eigen::MatrixXd GetFullCovMx ()
 
Eigen::MatrixXd GetCovMxAbsolute ()
 
Eigen::MatrixXd GetCovMxRelative (std::vector< Sample > &samples, osc::IOscCalcAdjustable *calc)
 
Eigen::MatrixXd GetFullCorrMx ()
 
void SetName (std::string name)
 
std::string GetName ()
 
std::vector< covmx::SampleGetSamples ()
 
std::vector< BinningGetBinnings ()
 
Binning GetBinning ()
 
Binning GetFullBinning ()
 
std::vector< unsigned intGetNBins ()
 
std::vector< unsigned intGetNComponents ()
 
void SetNUniverses (unsigned int i)
 
TH2D * GetFullCovMxTH2 ()
 
TH2D * GetCovMxAbsoluteTH2 ()
 
TH2D * GetCovMxRelativeTH2 (std::vector< Sample > &samples, osc::IOscCalcAdjustable *calc)
 
TH2D * GetCorrMxTH2 ()
 
TH2D * GetFullCorrMxTH2 ()
 
virtual void SaveTo (TDirectory *dir, const std::string &name) const
 
 CovarianceMatrix ()=delete
 

Static Public Member Functions

static Eigen::MatrixXd ForcePosDef (Eigen::MatrixXd mx, double epsilon=0.1)
 
static TMatrixD ForcePosDefROOT (Eigen::MatrixXd mx, double epsilon=0.1)
 
static std::unique_ptr< CovarianceMatrixLoadFrom (TDirectory *dir, const std::string &name)
 

Protected Member Functions

 CovarianceMatrix (Eigen::MatrixXd fullCovMx, Eigen::MatrixXd covMxAbsolute, std::vector< std::vector< double >> edges, std::vector< unsigned int > comps)
 
void BuildFullCovMx (std::vector< Sample > samples, std::vector< const ISyst * > systs)
 
void SetBinning (std::vector< Sample > samples)
 

Protected Attributes

std::string fName
 String name identifier. More...
 
std::vector< std::pair< Binning, unsigned int > > fBins
 
Eigen::MatrixXd fFullCovMx
 Full covariance matrix. More...
 
Eigen::MatrixXd fCovMxAbsolute
 Oscillated absolute covariance matrix. More...
 
unsigned int fNUniverses
 Number of universes to simulate. More...
 

Friends

class CovMxManager
 

Detailed Description

A class for generating a covariance matrices as a function of oscillation parameters and systematics as nuisance parameters.

Definition at line 29 of file CovarianceMatrix.h.

Constructor & Destructor Documentation

virtual ana::covmx::CovarianceMatrix::~CovarianceMatrix ( )
inlinevirtual
ana::covmx::CovarianceMatrix::CovarianceMatrix ( std::vector< Sample samples,
std::vector< const ISyst * >  systs,
int  nUniverses = 1e5 
)

Definition at line 75 of file CovarianceMatrix.cxx.

References ana::bins, BuildFullCovMx(), CovarianceMatrix(), fBins, fCovMxAbsolute, fFullCovMx, fNUniverses, GetBinning(), plotROC::nBins, ana::Binning::NBins(), and SetBinning().

77  : fNUniverses(nUniverses)
78  {
79  BuildFullCovMx(samples, systs);
80 
81  } // function CovarianceMatrix::CovarianceMatrix
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void BuildFullCovMx(std::vector< Sample > samples, std::vector< const ISyst * > systs)
unsigned int fNUniverses
Number of universes to simulate.
ana::covmx::CovarianceMatrix::CovarianceMatrix ( Eigen::MatrixXd  fullCovMx,
Eigen::MatrixXd  covMxAbsolute,
std::vector< std::pair< Binning, unsigned int >>  bins 
)
ana::covmx::CovarianceMatrix::CovarianceMatrix ( Eigen::MatrixXd  fullCovMx,
std::vector< Sample samples 
)
ana::covmx::CovarianceMatrix::CovarianceMatrix ( std::vector< Sample samples)

Constructor that makes an empty matrix from given bins.

Definition at line 104 of file CovarianceMatrix.cxx.

References fCovMxAbsolute, fFullCovMx, GetBinning(), GetFullBinning(), plotROC::nBins, ana::Binning::NBins(), and SetBinning().

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
int nBins
Definition: plotROC.py:16
Eigen::MatrixXd fCovMxAbsolute
Oscillated absolute covariance matrix.
Eigen::MatrixXd fFullCovMx
Full covariance matrix.
int NBins() const
Definition: Binning.h:29
void SetBinning(std::vector< Sample > samples)
unsigned int fNUniverses
Number of universes to simulate.
ana::covmx::CovarianceMatrix::CovarianceMatrix ( )
delete
ana::covmx::CovarianceMatrix::CovarianceMatrix ( Eigen::MatrixXd  fullCovMx,
Eigen::MatrixXd  covMxAbsolute,
std::vector< std::vector< double >>  edges,
std::vector< unsigned int comps 
)
protected

Member Function Documentation

void ana::covmx::CovarianceMatrix::AddMatrix ( CovarianceMatrix gen)

Definition at line 316 of file CovarianceMatrix.cxx.

References fCovMxAbsolute, fFullCovMx, GetBinning(), GetCovMxAbsolute(), GetFullBinning(), GetFullCovMx(), plotROC::nBins, and ana::Binning::NBins().

Referenced by MakePlots(), and ~CovarianceMatrix().

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();
328  fCovMxAbsolute += gen->GetCovMxAbsolute();
329 
330  } // function CovarianceMatrix::AddMatrix
int nBins
Definition: plotROC.py:16
Eigen::MatrixXd fCovMxAbsolute
Oscillated absolute covariance matrix.
Eigen::MatrixXd fFullCovMx
Full covariance matrix.
int NBins() const
Definition: Binning.h:29
void ana::covmx::CovarianceMatrix::BuildFullCovMx ( std::vector< Sample samples,
std::vector< const ISyst * >  systs 
)
protected

Definition at line 119 of file CovarianceMatrix.cxx.

References ana::bins, calc, comps, ana::Progress::Done(), stan::math::fabs(), fBins, fCovMxAbsolute, fFullCovMx, fNUniverses, GetBinning(), GetBinnings(), ana::covmx::GetComponents(), GetFullBinning(), ana::GetPOT(), analysePickle::hist, MECModelEnuComparisons::i, std::isnan(), ana::KeySyst::IsOneSided(), make_pair(), plotROC::nBins, ana::Binning::NBins(), file_size_ana::progress, generate_hists::rnd, ana::Progress::SetProgress(), ana::SystShifts::SetShift(), tmp, and ana::Spectrum::ToTH1().

Referenced by CovarianceMatrix(), and SetNUniverses().

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
int nBins
Definition: plotROC.py:16
vector< Component > GetComponents(Sample sample)
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Eigen::MatrixXd fCovMxAbsolute
Oscillated absolute covariance matrix.
std::vector< double > Spectrum
Definition: Constants.h:743
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
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
const XML_Char * s
Definition: expat.h:262
Eigen::MatrixXd fFullCovMx
Full covariance matrix.
std::vector< Binning > GetBinnings()
std::vector< std::string > comps
const Binning bins
Definition: NumuCC_CPiBin.h:8
std::vector< std::pair< Binning, unsigned int > > fBins
int NBins() const
Definition: Binning.h:29
double GetPOT(bool isFHC=true)
unsigned int fNUniverses
Number of universes to simulate.
MatrixXd ana::covmx::CovarianceMatrix::ForcePosDef ( Eigen::MatrixXd  mx,
double  epsilon = 0.1 
)
static

Definition at line 333 of file CovarianceMatrix.cxx.

References epsilon, and MECModelEnuComparisons::i.

Referenced by ~CovarianceMatrix().

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
double epsilon
static TMatrixD ana::covmx::CovarianceMatrix::ForcePosDefROOT ( Eigen::MatrixXd  mx,
double  epsilon = 0.1 
)
static

Referenced by ~CovarianceMatrix().

Binning ana::covmx::CovarianceMatrix::GetBinning ( )

Definition at line 358 of file CovarianceMatrix.cxx.

References b, ana::bins, ana::Binning::Custom(), GetBinnings(), and MECModelEnuComparisons::i.

Referenced by AddMatrix(), BuildFullCovMx(), ana::LikelihoodCovMxExperiment::ChiSq(), ana::OscCovMxExperiment::ChiSq(), CovarianceMatrix(), GetCorrMxTH2(), GetCovMxAbsoluteTH2(), GetCovMxRelativeTH2(), GetName(), and ana::LikelihoodCovMxExperiment::SaveHists().

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
std::vector< Binning > GetBinnings()
static Binning Custom(const std::vector< double > &edges)
Definition: Binning.cxx:145
const Binning bins
Definition: NumuCC_CPiBin.h:8
const hit & b
Definition: hits.cxx:21
vector< Binning > ana::covmx::CovarianceMatrix::GetBinnings ( )

Definition at line 348 of file CovarianceMatrix.cxx.

References febshutoff_auto::end, fBins, runNovaSAM::ret, and PandAna.Demos.pi0_spectra::transform.

Referenced by BuildFullCovMx(), GetBinning(), and GetName().

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
std::vector< std::pair< Binning, unsigned int > > fBins
TH2D * ana::covmx::CovarianceMatrix::GetCorrMxTH2 ( )

Definition at line 533 of file CovarianceMatrix.cxx.

References ana::Binning::Edges(), fCovMxAbsolute, GetBinning(), analysePickle::hist, MECModelEnuComparisons::i, calib::j, ana::Binning::NBins(), std::sqrt(), and ana::UniqueName().

Referenced by SetNUniverses().

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
T sqrt(T number)
Definition: d0nt_math.hpp:156
Eigen::MatrixXd fCovMxAbsolute
Oscillated absolute covariance matrix.
const double j
Definition: BetheBloch.cxx:29
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
Eigen::MatrixXd ana::covmx::CovarianceMatrix::GetCovMxAbsolute ( )
inline

Definition at line 62 of file CovarianceMatrix.h.

References fCovMxAbsolute, GetCovMxRelative(), and GetFullCorrMx().

Referenced by AddMatrix().

62 { return fCovMxAbsolute; }
Eigen::MatrixXd fCovMxAbsolute
Oscillated absolute covariance matrix.
TH2D * ana::covmx::CovarianceMatrix::GetCovMxAbsoluteTH2 ( )

Definition at line 489 of file CovarianceMatrix.cxx.

References ana::bins, ana::Binning::Edges(), fCovMxAbsolute, GetBinning(), analysePickle::hist, MECModelEnuComparisons::i, calib::j, plotROC::nBins, ana::Binning::NBins(), and ana::UniqueName().

Referenced by SetNUniverses(), and TestCovMxNew().

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
int nBins
Definition: plotROC.py:16
Eigen::MatrixXd fCovMxAbsolute
Oscillated absolute covariance matrix.
const double j
Definition: BetheBloch.cxx:29
const Binning bins
Definition: NumuCC_CPiBin.h:8
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
MatrixXd ana::covmx::CovarianceMatrix::GetCovMxRelative ( std::vector< Sample > &  samples,
osc::IOscCalcAdjustable calc 
)

Definition at line 414 of file CovarianceMatrix.cxx.

References demo5::c1, demo5::c2, fCovMxAbsolute, h1, h2, Predict(), and runNovaSAM::ret.

Referenced by GetCovMxAbsolute(), and GetCovMxRelativeTH2().

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
Eigen::MatrixXd fCovMxAbsolute
Oscillated absolute covariance matrix.
c2
Definition: demo5.py:33
TH1F * h2
Definition: plot.C:45
TH1F * h1
c1
Definition: demo5.py:24
Eigen::MatrixXd Predict(std::vector< Sample > samples, osc::IOscCalcAdjustable *calc, const SystShifts &systs=SystShifts::Nominal())
TH2D * ana::covmx::CovarianceMatrix::GetCovMxRelativeTH2 ( std::vector< Sample > &  samples,
osc::IOscCalcAdjustable calc 
)

Definition at line 510 of file CovarianceMatrix.cxx.

References ana::bins, ana::Binning::Edges(), GetBinning(), GetCovMxRelative(), analysePickle::hist, MECModelEnuComparisons::i, calib::j, plotROC::nBins, ana::Binning::NBins(), and ana::UniqueName().

Referenced by SetNUniverses(), and TestCovMxNew().

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
int nBins
Definition: plotROC.py:16
Eigen::MatrixXd GetCovMxRelative(std::vector< Sample > &samples, osc::IOscCalcAdjustable *calc)
const double j
Definition: BetheBloch.cxx:29
const Binning bins
Definition: NumuCC_CPiBin.h:8
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
Binning ana::covmx::CovarianceMatrix::GetFullBinning ( )

Definition at line 374 of file CovarianceMatrix.cxx.

References b, ana::bins, ana::Binning::Custom(), and fBins.

Referenced by AddMatrix(), BuildFullCovMx(), CovarianceMatrix(), GetFullCorrMxTH2(), GetFullCovMxTH2(), GetName(), Predict(), and ana::LikelihoodCovMxExperiment::SaveHists().

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
static Binning Custom(const std::vector< double > &edges)
Definition: Binning.cxx:145
const Binning bins
Definition: NumuCC_CPiBin.h:8
std::vector< std::pair< Binning, unsigned int > > fBins
const hit & b
Definition: hits.cxx:21
MatrixXd ana::covmx::CovarianceMatrix::GetFullCorrMx ( )

Definition at line 444 of file CovarianceMatrix.cxx.

References fFullCovMx, MECModelEnuComparisons::i, calib::j, and std::sqrt().

Referenced by GetCovMxAbsolute(), and GetFullCorrMxTH2().

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
T sqrt(T number)
Definition: d0nt_math.hpp:156
Eigen::MatrixXd fFullCovMx
Full covariance matrix.
const double j
Definition: BetheBloch.cxx:29
TH2D * ana::covmx::CovarianceMatrix::GetFullCorrMxTH2 ( )

Definition at line 558 of file CovarianceMatrix.cxx.

References ana::Binning::Edges(), GetFullBinning(), GetFullCorrMx(), analysePickle::hist, MECModelEnuComparisons::i, calib::j, ana::Binning::NBins(), and ana::UniqueName().

Referenced by SetNUniverses().

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
const double j
Definition: BetheBloch.cxx:29
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
Eigen::MatrixXd ana::covmx::CovarianceMatrix::GetFullCovMx ( )
inline

Definition at line 61 of file CovarianceMatrix.h.

References fFullCovMx.

Referenced by AddMatrix(), and ana::LikelihoodCovMxExperiment::LikelihoodCovMxExperiment().

61 { return fFullCovMx; }
Eigen::MatrixXd fFullCovMx
Full covariance matrix.
TH2D * ana::covmx::CovarianceMatrix::GetFullCovMxTH2 ( )

Definition at line 468 of file CovarianceMatrix.cxx.

References ana::bins, ana::Binning::Edges(), fFullCovMx, GetFullBinning(), analysePickle::hist, MECModelEnuComparisons::i, calib::j, plotROC::nBins, ana::Binning::NBins(), and ana::UniqueName().

Referenced by SetNUniverses(), and TestCovMxNew().

469  {
470  DontAddDirectory guard;
471  // Get full binning and create histogram
472  Binning bins = GetFullBinning();
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
int nBins
Definition: plotROC.py:16
Eigen::MatrixXd fFullCovMx
Full covariance matrix.
const double j
Definition: BetheBloch.cxx:29
const Binning bins
Definition: NumuCC_CPiBin.h:8
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
std::string ana::covmx::CovarianceMatrix::GetName ( )
inline

Definition at line 68 of file CovarianceMatrix.h.

References fName, GetBinning(), GetBinnings(), GetFullBinning(), GetNBins(), GetNComponents(), and GetSamples().

68 { return fName; }
std::string fName
String name identifier.
vector< unsigned int > ana::covmx::CovarianceMatrix::GetNBins ( )

Definition at line 394 of file CovarianceMatrix.cxx.

References febshutoff_auto::end, fBins, runNovaSAM::ret, and PandAna.Demos.pi0_spectra::transform.

Referenced by GetName().

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
std::vector< std::pair< Binning, unsigned int > > fBins
vector< unsigned int > ana::covmx::CovarianceMatrix::GetNComponents ( )

Definition at line 404 of file CovarianceMatrix.cxx.

References febshutoff_auto::end, fBins, runNovaSAM::ret, and PandAna.Demos.pi0_spectra::transform.

Referenced by GetName(), and Predict().

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
std::vector< std::pair< Binning, unsigned int > > fBins
std::vector<covmx::Sample> ana::covmx::CovarianceMatrix::GetSamples ( )

Referenced by GetName().

std::vector<Spectrum> ana::covmx::CovarianceMatrix::GetUniverse ( osc::OscCalcSterile calc,
std::vector< covmx::Sample samples 
)

Referenced by ~CovarianceMatrix().

unique_ptr< CovarianceMatrix > ana::covmx::CovarianceMatrix::LoadFrom ( TDirectory *  dir,
const std::string name 
)
static

Definition at line 639 of file CovarianceMatrix.cxx.

References ana::assert(), ana::bins, ana::Binning::Custom(), dir, it, make_pair(), runNovaSAM::ret, getGoodRuns4SAM::tag, and tmp.

Referenced by SetNUniverses().

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
const XML_Char * name
Definition: expat.h:151
set< int >::iterator it
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
static Binning Custom(const std::vector< double > &edges)
Definition: Binning.cxx:145
const Binning bins
Definition: NumuCC_CPiBin.h:8
TDirectory * dir
Definition: macro.C:5
assert(nhit_max >=nhit_nbins)
MatrixXd ana::covmx::CovarianceMatrix::Predict ( std::vector< Sample samples,
osc::IOscCalcAdjustable calc,
const SystShifts systs = SystShifts::Nominal() 
)

Definition at line 223 of file CovarianceMatrix.cxx.

References comps, fBins, fCovMxAbsolute, fFullCovMx, ana::Spectrum::GetBinnings(), ana::covmx::GetComponents(), GetFullBinning(), GetNComponents(), ana::GetPOT(), analysePickle::hist, MECModelEnuComparisons::i, calib::j, ana::Binning::NBins(), and ana::Spectrum::ToTH1().

Referenced by ana::OscCovMxExperiment::EvaluateChiSq(), GetCovMxRelative(), TestCovMxNew(), and ~CovarianceMatrix().

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
vector< Component > GetComponents(Sample sample)
Eigen::MatrixXd fCovMxAbsolute
Oscillated absolute covariance matrix.
std::vector< double > Spectrum
Definition: Constants.h:743
Eigen::MatrixXd fFullCovMx
Full covariance matrix.
const double j
Definition: BetheBloch.cxx:29
std::vector< std::string > comps
std::vector< unsigned int > GetNComponents()
std::vector< std::pair< Binning, unsigned int > > fBins
int NBins() const
Definition: Binning.h:29
double GetPOT(bool isFHC=true)
MatrixXd ana::covmx::CovarianceMatrix::Predict ( std::vector< double >  vBetaMu)

Definition at line 283 of file CovarianceMatrix.cxx.

References fBins, fFullCovMx, GetNComponents(), MECModelEnuComparisons::i, calib::j, and runNovaSAM::ret.

283  {
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
Eigen::MatrixXd fFullCovMx
Full covariance matrix.
const double j
Definition: BetheBloch.cxx:29
std::vector< unsigned int > GetNComponents()
std::vector< std::pair< Binning, unsigned int > > fBins
void ana::covmx::CovarianceMatrix::SaveTo ( TDirectory *  dir,
const std::string name 
) const
virtual

Definition at line 601 of file CovarianceMatrix.cxx.

References dir, fBins, fCovMxAbsolute, fFullCovMx, MECModelEnuComparisons::i, it, calib::j, and tmp.

Referenced by SetNUniverses(), and TestCovMxNew().

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
const XML_Char * name
Definition: expat.h:151
set< int >::iterator it
Eigen::MatrixXd fCovMxAbsolute
Oscillated absolute covariance matrix.
Float_t tmp
Definition: plot.C:36
Eigen::MatrixXd fFullCovMx
Full covariance matrix.
const double j
Definition: BetheBloch.cxx:29
std::vector< std::pair< Binning, unsigned int > > fBins
TDirectory * dir
Definition: macro.C:5
void ana::covmx::CovarianceMatrix::SetBinning ( std::vector< Sample samples)
protected

Definition at line 584 of file CovarianceMatrix.cxx.

References fBins, ana::covmx::GetComponents(), and make_pair().

Referenced by CovarianceMatrix(), and SetNUniverses().

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
vector< Component > GetComponents(Sample sample)
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
const XML_Char * s
Definition: expat.h:262
std::vector< std::pair< Binning, unsigned int > > fBins
void ana::covmx::CovarianceMatrix::SetName ( std::string  name)
inline

Definition at line 67 of file CovarianceMatrix.h.

References fName.

Referenced by MakePlots().

67 { fName = name; }
const XML_Char * name
Definition: expat.h:151
std::string fName
String name identifier.
void ana::covmx::CovarianceMatrix::SetNUniverses ( unsigned int  i)
inline

Friends And Related Function Documentation

friend class CovMxManager
friend

Definition at line 115 of file CovarianceMatrix.h.

Member Data Documentation

std::vector<std::pair<Binning, unsigned int> > ana::covmx::CovarianceMatrix::fBins
protected
Eigen::MatrixXd ana::covmx::CovarianceMatrix::fCovMxAbsolute
protected

Oscillated absolute covariance matrix.

Definition at line 109 of file CovarianceMatrix.h.

Referenced by AddMatrix(), BuildFullCovMx(), CovarianceMatrix(), GetCorrMxTH2(), GetCovMxAbsolute(), GetCovMxAbsoluteTH2(), GetCovMxRelative(), Predict(), and SaveTo().

Eigen::MatrixXd ana::covmx::CovarianceMatrix::fFullCovMx
protected

Full covariance matrix.

Definition at line 108 of file CovarianceMatrix.h.

Referenced by AddMatrix(), BuildFullCovMx(), CovarianceMatrix(), GetFullCorrMx(), GetFullCovMx(), GetFullCovMxTH2(), Predict(), and SaveTo().

std::string ana::covmx::CovarianceMatrix::fName
protected

String name identifier.

Definition at line 105 of file CovarianceMatrix.h.

Referenced by GetName(), and SetName().

unsigned int ana::covmx::CovarianceMatrix::fNUniverses
protected

Number of universes to simulate.

Definition at line 111 of file CovarianceMatrix.h.

Referenced by BuildFullCovMx(), CovarianceMatrix(), and SetNUniverses().


The documentation for this class was generated from the following files: