OscCovMxExperiment.cxx
Go to the documentation of this file.
2 
5 
6 #include "OscLib/IOscCalc.h"
7 
8 #include "TDirectory.h"
9 #include "TObjString.h"
10 #include "TH1.h"
11 #include "TMatrixDEigen.h"
12 #include "TDecompChol.h"
13 
14 #include <iomanip>
15 #include <iostream>
16 #include <chrono>
17 
18 #include <Eigen/Dense>
19 
20 using Eigen::MatrixXd;
21 using Eigen::VectorXd;
22 
23 namespace ana
24 {
25  //----------------------------------------------------------------------
26  OscCovMxExperiment::OscCovMxExperiment(std::vector<covmx::Sample> samples,
28  : fSamples(samples), fCovMx(covmx), fIncludeStatError(true),
29  fCombinedNeymanPearson(false), fShapeOnly(false), fVerbose(false),
30  fInversionTimeHist(nullptr)
31  {
34  fNBins = 0;
35  fNBinsFull = 0;
36  for (covmx::Sample& s: fSamples) {
37  fBins.push_back(s.GetBinning());
38  fNBins += s.GetBinning().NBins();
39  fNBinsFull += s.GetBinning().NBins() * GetComponents(s).size();
40  // fData.push_back(s.GetData());
41  // fCosmic.push_back(s.HasCosmic()? s.GetCosmic() : nullptr);
42  }
43  delete calc;
44 
45  fUnoscillated.resize(fSamples.size(), nullptr);
46  }
47 
48  //----------------------------------------------------------------------
50  {}
51 
52  //----------------------------------------------------------------------
54  const SystShifts& syst) const
55  {
56  DontAddDirectory guard;
57 
58  // Get MC and data spectra
59  std::vector<double> pot(fSamples.size());
60  std::vector<TH1D*> hData(fSamples.size());
61  std::vector<TH1D*> hPred(fSamples.size());
62 
63  for (size_t i = 0; i < fSamples.size(); ++i) {
64 
65  pot[i] = fSamples[i].GetData()->POT();
66  hData[i] = fSamples[i].GetData()->ToTH1(pot[i]);
67 
68  SystShifts shift = fSamples[i].GetSystShifts(syst);
69 
70  if (fUnoscillated[i]) hPred[i] = (TH1D*)fUnoscillated[i]->Clone();
71  else {
72  const Spectrum pred = fSamples[i].GetPrediction()->PredictSyst(osc, shift);
73  hPred[i] = pred.ToTH1(pot[i]);
74  }
75 
76  if (fSamples[i].HasCosmic()) {
77  hPred[i]->Add(fSamples[i].GetCosmic()->ToTH1(fSamples[i].GetData()->Livetime(), kLivetime));
78  }
79 
80  if (fShapeOnly) {
81  hPred[i]->Scale(hData[i]->Integral()/hPred[i]->Integral());
82  }
83  }
84 
85  // Number of bins
86  if (fCovMx && fNBins != (size_t)fCovMx->GetBinning().NBins()) {
87  throw std::runtime_error(
88  "Number of bins in predictions does not match covariance matrix!");
89  }
90 
91  // Define difference vector and covariance matrix
92  VectorXd difference(fNBins);
93  MatrixXd covMx(fNBins, fNBins);
94  if (not fCovMx) covMx.setZero();
95  else covMx = fCovMx->Predict(fSamples, osc);//, fSystConcat, syst);
96 
97  size_t i = 0; // loop over every bin
98  for (size_t iPred = 0; iPred < fSamples.size(); ++iPred) {
99  for (size_t iBin = 1; iBin <= (size_t)fBins[iPred].NBins(); ++iBin) {
100  difference(i) = hData[iPred]->GetBinContent(iBin)
101  - hPred[iPred]->GetBinContent(iBin);
102  if (fIncludeStatError) { // add in statistical errors
104  covMx(i,i) += 3 / ((1./hData[iPred]->GetBinContent(iBin))+(2./hPred[iPred]->GetBinContent(iBin)));
105  } else {
106  covMx(i,i) += hPred[iPred]->GetBinContent(iBin);
107  }
108  } // if including statistical errors
109  ++i;
110  } // for bin i
111  } // for sample i
112 
113  // Decompose and invert matrix
114 
115  MatrixXd covMxInv = covMx.inverse();
116  double chi2 = (difference*covMxInv*difference).value();
117 
118  // Debug output
119  if (fVerbose) {
120  for (size_t i = 0; i < fNBins; ++i) {
121  for (size_t j = 0; j < fNBins; ++j) {
122  std::cout << "Chi^2 contribution: " << difference(i)
123  << " x " << covMxInv(i,j) << " x " << difference(j)
124  << " = " << difference(i) * covMxInv(i,j) * difference(j)
125  << std::endl;
126  }
127  }
128  }
129 
130  // Clean up and return
131  for (auto h : hPred) delete h;
132  for (auto h : hData) delete h;
133 
134  return chi2;
135 
136  } // function OscCovMxExperiment::ChiSq
137 
138  //----------------------------------------------------------------------
139  double OscCovMxExperiment::EvaluateChiSq(std::vector<double> vBetaMu,
140  std::vector<double> vPred, std::vector<double> vData) const {
141 
142  // Set up inputs
143  MatrixXd mx = fCovMx->Predict(vBetaMu);
144  VectorXd difference(fNBins);
145 
146  // Handle differemce vector and statistical uncertainties
147  for (size_t i = 0; i < vData.size(); ++i) {
148  difference(i) = vPred[i] - vData[i];
149  if (fIncludeStatError) { // add in statistical errors
151  mx(i,i) += 3 / ((1./vData[i])+(2./vPred[i]));
152  } else {
153  mx(i,i) += vPred[i];
154  }
155  } // if including statistical error
156  } // for bin
157 
158  return difference.transpose() * mx.inverse() * difference;
159 
160  } // function OscCovMxExperiment::EvaluateChiSq
161 
162  //----------------------------------------------------------------------
164  {
165  if( !fInversionTimeHist ){
166  std::cerr << "Inversion Time histogram not configured to fill, "
167  << "use OscCovMxExperiment::SetInversionTimeHist() before running the fit"
168  << std::endl;
169  abort();
170  }
171  return fInversionTimeHist;
172  }
173  //----------------------------------------------------------------------
174  void OscCovMxExperiment::SetInversionTimeHist(bool timeit, int tmax, int tmin)
175  {
176  if( !timeit ){
177  fInversionTimeHist = nullptr;
178  } else {
179  fInversionTimeHist = new TH1I("inversion_time",
180  ";Single Inversion Time (ms);Number of Inversions",
181  1000, tmin, tmax);
182  }
183  }
184 
185 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< std::tuple< ana::Flavors::Flavors_t, ana::Current::Current_t, ana::Sign::Sign_t > > GetComponents()
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
double Integral(const Spectrum &s, const double pot, int cvnregion)
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
std::vector< TH1D * > fUnoscillated
double EvaluateChiSq(std::vector< double > vBetaMu, std::vector< double > vPred, std::vector< double > vData) const
OStream cerr
Definition: OStream.cxx:7
Version of OscCalcSterile that always returns probability of 1.
osc::OscCalcDumb calc
double chi2()
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void SetInversionTimeHist(bool timeit=true, int tmax=1000, int tmin=0)
const XML_Char * s
Definition: expat.h:262
const XML_Char int const XML_Char * value
Definition: expat.h:331
#define pot
const double j
Definition: BetheBloch.cxx:29
Oscillation probability calculators.
Definition: Calcs.h:5
OStream cout
Definition: OStream.cxx:6
covmx::CovarianceMatrix * fCovMx
std::vector< Binning > fBins
int NBins() const
Definition: Binning.h:29
std::vector< covmx::Sample > fSamples
OscCovMxExperiment(std::vector< covmx::Sample > samples, covmx::CovarianceMatrix *covmx=nullptr)
virtual double ChiSq(osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const
Eigen::MatrixXd Predict(std::vector< Sample > samples, osc::IOscCalcAdjustable *calc, const SystShifts &systs=SystShifts::Nominal())
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
A class for generating a covariance matrices as a function of oscillation parameters and systematics ...