CovMxExperiment.cxx
Go to the documentation of this file.
1 /*
2  * NumuCovMx.cxx
3  *
4  * Created on: Jan 10, 2017
5  * Author: J. Wolcott <jwolcott@fnal.gov>
6  */
7 
9 
10 #include "TCanvas.h"
11 #include "TH1D.h"
12 #include "THStack.h"
13 #include "TMatrixD.h"
14 #include "TObjString.h"
15 
16 #include "CAFAna/Core/HistCache.h"
17 
19 
20 
21 namespace ana
22 {
24  {
25  // the covariance matrix replaces the systematic shift approach.
26  if (!syst.IsNominal())
27  throw std::runtime_error("ExperimentWithCovMx::ChiSq(): Systematic shifts cannot be specified when fitting against covariance matrices");
28 
29  const double pot = fData.POT();
30  TH1D* hdata = fData.ToTH1(pot);
31 
32  std::map<Component, TH1D*> predictions;
33  for (const auto & mxPair : this->fCovMxInv)
34  {
35  const auto & comp = mxPair.first;
36  if (comp.flavor == Flavors::kAll && comp.current == Current::kBoth && comp.sign == Sign::kBoth)
37  {
38  // if you specified one matrix as applying to the whole prediction,
39  // why did you give me more of them??
40  if (this->fCovMxInv.size() > 1)
41  throw std::runtime_error("CovMxExperiment::ChiSq(): one component corresponds to the whole prediction but you also provided others?!");
42 
43  predictions.emplace(comp, fMC->Predict(calc).ToTH1(pot));
44  }
45  else
46  predictions.emplace(mxPair.first, fMC->PredictComponent(calc, comp.flavor, comp.current, comp.sign).ToTH1(pot));
47  }
48 
49  // will want the total prediction in a moment,
50  // but since we assume the components add to the total
51  // in the following calculations, might as well ensure it here.
52  TH1D * totalPred = fMC->Predict(calc).ToTH1(pot);
53  TH1D diff(*totalPred);
54  for (const auto & predPair : predictions)
55  diff.Add(predPair.second, -1);
56  for (int bin = 0; bin <= totalPred->GetNbinsX() + 1; bin++)
57  {
58  if ( std::fabs(diff.GetBinContent(bin)) > 1e-4 )
59  throw std::runtime_error("CovMxExperiment::ChiSq(): components don't sum to total prediction!");
60  }
61 
62  // add the cosmics into the prediction.
63  if(fCosmic)
64  {
65  if(fCosmicScaleError != 0)
66  {
67  const double scale = 1 + syst.GetShift(&kCosmicBkgScaleSyst) * fCosmicScaleError;
68  totalPred->Add(fCosmic, scale);
69  }
70  else
71  totalPred->Add(fCosmic);
72 
73  }
74 
75  // first, compute the statistical uncertainty covariance matrix
76  // (which is diagonal). we'll need this repeatedly below.
77  // use the *total* prediction to get the sigma^2 here because
78  // that's the expected fluctuation size
79  TMatrixD stat_mx(totalPred->GetNbinsX(), totalPred->GetNbinsX());
80  for (int binNum = 1; binNum <= totalPred->GetNbinsX(); binNum++)
81  stat_mx[binNum-1][binNum-1] = totalPred->GetBinContent(binNum);
82 
83 
84  // We need to include the statistical covariance matrix
85  // to do the chi^2 calculation or it'll be way off.
86  // we need to recalculate the sum for *every* new configuration of
87  // the IPrediction, unfortunately--
88  // and since we have to then invert it, this could be very costly.
89  // StackExchange to the rescue: somebody there wanted to solve a similar problem
90  // (where we already have the inverse of a matrix A
91  // and want to compute (A+B)^{-1] with some other invertible B):
92  // http://math.stackexchange.com/a/17780
93  // method has apparently been discovered in various forms by multiple people,
94  // but this version is due to K. Miller, Math.Mag. 54, 67 (1981)
95  //
96 
97 
98  // Strictly speaking, it's not possible to divide the prediction up
99  // into pieces and then use the covariance matrices for those pieces
100  // to compute a chi^2 because it's not possible to know which was which
101  // in the data. However, we'll cheat a little and do the following:
102  // for each component, we'll subtract the rest of the prediction
103  // and compare background-subtracted data to that component.
104  // Because we have to include the statistical covariance matrix every time,
105  // this over-counts the statistical uncertainty by a factor of
106  // roughly (N-1), where N is the number of components.
107  // Probably this isn't exactly the right statistical thing to do,
108  // but for a first approximation, we'll take an average of the chi^2s
109  // weighted by how many events they contribute to the prediction.
110 
111  double chi2 = 0;
112  for (auto it_comp = this->fCovMxInv.begin(); it_comp != this->fCovMxInv.end(); ++it_comp)
113  {
114  const auto & comp = it_comp->first;
115  const auto & compCovMxInv = it_comp->second;
116  const auto & predComp = *(predictions.at(comp));
117  TH1D hdata_compSub = *hdata - ( *totalPred - predComp );
118 
119 
120  // now recursively update the matrix
121  // by calculating the inverse of (A+B_k),
122  // where B_k is one of the N rank-1 matrices
123  // (when N is the rank of B) that B can be decomposed into.
124  // our case is even easier because we know B is diagonal
125  // (it's the statistical uncertainty covariance matrix).
126  TMatrixD Bk(stat_mx.GetNrows(), stat_mx.GetNcols());
127  TMatrixD CkInv(*compCovMxInv); // here we're setting what Miller calls C_1 = A
128 
129  // the rank of the stat mx is just the number of rows:
130  // it's diagonal and the bins are orthogonal.
131  // note that we're counting from zero instead of 1
132  // for convenience in the TMatrixDs
133  for (int k = 0; k < stat_mx.GetNrows(); k++)
134  {
135  // reset the last iteration so we can re-use
136  if (k > 0)
137  Bk[k-1][k-1] = 0;
138 
139  // set current iteration B_k
140  Bk[k][k] = stat_mx[k][k];
141 
142  // can short-circuit the trace in Miller's algorithm
143  // since our Bk come from a diagonal matrix:
144  // they're single entries on the diagonal
145  double gk = 1. / (1 + CkInv[k][k]*Bk[k][k]);
146 
147  CkInv -= gk * CkInv * Bk * CkInv;
148  }
149 
150  // event-weighted chi^2
151  chi2 += predComp.Integral() * Chi2CovMx(predictions.at(comp), &hdata_compSub, &CkInv);
152 
153  }
154 
155  // divide out the total number of events to the get the weighted chi^2
156  chi2 /= totalPred->Integral();
157 
158 
159  // finally, clean up all the temp hists we got from Spectrum::ToTH1()
160  HistCache::Delete(hdata);
161  HistCache::Delete(totalPred);
162  for (auto & predPair : predictions)
163  HistCache::Delete(predPair.second);
164 
165 
166  return chi2;
167  }
168 
169 }
170 
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
bool IsNominal() const
Definition: SystShifts.h:43
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:553
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const CosmicBkgScaleSyst kCosmicBkgScaleSyst
double chi2()
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
Double_t scale
Definition: plot.C:25
General interface to any calculator that lets you set the parameters.
Interactions of both types.
Definition: IPrediction.h:42
static void Delete(TH1D *&h)
Definition: HistCache.cxx:215
#define pot
T GetShift(const ISyst *syst) const
virtual double ChiSq(osc::IOscCalculatorAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const override
float bin[41]
Definition: plottest35.C:14
double POT() const
Definition: Spectrum.h:263
predictions
std::map< Component, std::unique_ptr< TMatrixD > > fCovMxInv
string syst
Definition: plotSysts.py:176
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
All neutrinos, any flavor.
Definition: IPrediction.h:26
Float_t e
Definition: plot.C:35
double Chi2CovMx(const TVectorD *e, const TVectorD *o, const TMatrixD *covmxinv)
Definition: Utilities.cxx:176