TwoSampleDecomp.cxx
Go to the documentation of this file.
2 
3 #include <iostream>
4 #include <cmath>
5 #include <cassert>
6 #include <exception>
7 #include <iostream>
8 #include <string>
9 #include <sstream>
10 
11 #include "CAFAna/Core/Spectrum.h"
12 
13 #include "TH1D.h"
14 #include "TDirectory.h"
15 #include "TObjString.h"
16 #include "Utilities/func/MathUtil.h"
17 #include <iostream>
18 
19 namespace ana
20 {
21 
23  SpectrumLoaderBase& loaderMC, ///< Spectrum loader for MC
24  SpectrumLoaderBase& loaderData,///< Spectrum loader for data
25  std::string label, ///< Name of the thing
26  const Binning& bins, ///< Desired binning for spectra
27  const Var& var, ///< Variable to be plotted in spectra
28  const Cut& firstSample, ///< Cut to define samples,
29  ///< true: first, false: second
30  const Cut& isSignal, ///< Signal or not? true:sig, false: bkg
31  const Cut& anaCut, ///< Cut to define selection (e.g. numu CC)
32  const SystShifts& shift,
33  const Var& weight ///< Weight for events in spectrum
34  ):TwoSampleDecomp(loaderMC, loaderData, HistAxis(label, bins, var),
35  firstSample, isSignal, anaCut, shift, weight)
36  {}
37 
38 
39 
41  SpectrumLoaderBase& loaderMC, ///< Spectrum loader for MC
42  SpectrumLoaderBase& loaderData, ///< Spectrum loader for data
43  const HistAxis& axis, ///< Variable to be plotted in spectra
44  const Cut& firstSample, ///< Cut to define samples,
45  ///< true: first, false: second
46  const Cut& isSignal, ///< Signal or not? true:sig, false: bkg
47  const Cut& anaCut, ///< Cut to define selection (e.g. numu CC)
48  const SystShifts& shift,
49  const Var& weight ///< Weight for events in spectrum
50  ):
51  fLoaderMC (loaderMC),
52  fLoaderData(loaderData),
53  fAxis(&axis),
54  fFirstSig (loaderMC, axis, firstSample && isSignal && anaCut,
55  shift, weight),
56  fSecondSig(loaderMC, axis, !firstSample && isSignal && anaCut,
57  shift, weight),
58  fFirstBkg (loaderMC, axis, firstSample && !isSignal && anaCut,
59  shift, weight),
60  fSecondBkg(loaderMC, axis, !firstSample && !isSignal && anaCut,
61  shift, weight),
62  fTrueSig (loaderMC, axis, isSignal && anaCut,
63  shift, weight),
64  fTrueBkg (loaderMC, axis, !isSignal && anaCut,
65  shift, weight),
66  fFirst (loaderData, axis, firstSample && anaCut,
67  shift, weight),
68  fSecond (loaderData, axis, !firstSample && anaCut,
69  shift, weight),
70  fSig (Spectrum::Uninitialized()),
71  fBkg (Spectrum::Uninitialized()),
72  fSigRatio (0),
73  fBkgRatio (0),
74  fDecomposed(false)
75  {
76  if(axis.NDimensions() != 1){
77  std::cout << "Error: TwoSampleDecomp not equipped to handle multi-dimensional axis." << std::endl;
78  abort();
79  }
80  }
81 
82 
83 
87  fFirstSig (Spectrum::Uninitialized()),
88  fSecondSig(Spectrum::Uninitialized()),
89  fFirstBkg (Spectrum::Uninitialized()),
90  fSecondBkg(Spectrum::Uninitialized()),
91  fTrueSig (Spectrum::Uninitialized()),
92  fTrueBkg (Spectrum::Uninitialized()),
93  fFirst (Spectrum::Uninitialized()),
94  fSecond (Spectrum::Uninitialized()),
95  fSig (Spectrum::Uninitialized()),
96  fBkg (Spectrum::Uninitialized()),
97  fSigRatio (0),
98  fBkgRatio (0),
99  fDecomposed(false)
100  { }
101 
102 
103 
104  void TwoSampleDecomp::SaveTo(TDirectory* dir, const std::string& name) const
105  {
106  TDirectory* tmp = gDirectory;
107 
108  dir = dir->mkdir(name.c_str()); // switch to subdir
109  dir->cd();
110 
111  TObjString("TwoSampleDecomp").Write("type");
112 
113  fFirstSig.SaveTo(dir, "FirstSig");
114  fSecondSig.SaveTo(dir, "SecondSig");
115  fFirstBkg.SaveTo(dir, "FirstBkg");
116  fSecondBkg.SaveTo(dir, "SecondBkg");
117  fTrueSig.SaveTo(dir, "TrueSig");
118  fTrueBkg.SaveTo(dir, "TrueBkg");
119  fFirst.SaveTo(dir, "First");
120  fSecond.SaveTo(dir, "Second");
121  fSig.SaveTo(dir, "Sig");
122  fBkg.SaveTo(dir, "Bkg");
123 
124  dir->Write();
125  delete dir;
126 
127  tmp->cd();
128  }
129 
130  //----------------------------------------------------------------------
131  std::unique_ptr<TwoSampleDecomp> TwoSampleDecomp::LoadFrom(TDirectory* dir, const std::string& name)
132  {
133  dir = dir->GetDirectory(name.c_str()); // switch to subdir
134  assert(dir);
135 
136  std::unique_ptr<TwoSampleDecomp> ret(new TwoSampleDecomp);
137 
138  ret->fFirstSig = *Spectrum::LoadFrom(dir, "FirstSig");
139  ret->fSecondSig = *Spectrum::LoadFrom(dir, "SecondSig");
140  ret->fFirstBkg = *Spectrum::LoadFrom(dir, "FirstBkg");
141  ret->fSecondBkg = *Spectrum::LoadFrom(dir, "SecondBkg");
142  ret->fTrueSig = *Spectrum::LoadFrom(dir, "TrueSig");
143  ret->fTrueBkg = *Spectrum::LoadFrom(dir, "TrueBkg");
144  ret->fFirst = *Spectrum::LoadFrom(dir, "First");
145  ret->fSecond = *Spectrum::LoadFrom(dir, "Second");
146  ret->fSig = *Spectrum::LoadFrom(dir, "Sig");
147  ret->fBkg = *Spectrum::LoadFrom(dir, "Bkg");
148 
149  delete dir;
150 
151  return ret;
152  }
153 
154  Spectrum TwoSampleDecomp::NCComponent() const {std::cout << "TwoSampleDecomp::NCComponent() not implemented" << std::endl; abort();}
155 
156  Spectrum TwoSampleDecomp::NCAntiComponent() const {std::cout << "TwoSampleDecomp::NCAntiComponent() not implemented" << std::endl; abort();}
157 
159  if(fDecomposed) return;
160 
162  && "Call Go() on MC SpectrumLoader before TwoSampleDecomp::Decomp");
164  && "Call Go() on Data SpectrumLoader before TwoSampleDecomp::Decomp");
165 
166  RatioCalc();
167 
168  TH1D* hFirst = fFirst.ToTH1(fFirst.POT());
169  TH1D* hSecond = fSecond.ToTH1(fSecond.POT());
170 
171  // Get the number of bins from the axis
172  int nBins = hFirst->GetNbinsX();
173 
174  // Make sure the histogram binning is sane
175  assert(hFirst->GetNbinsX() == nBins &&
176  "Inconsistent binning across spectra.");
177  assert(hSecond->GetNbinsX() == nBins &&
178  "Inconsistent binning across spectra.");
179  assert(fSigRatio->GetNbinsX() == nBins &&
180  "MC and data have different numbers of bins.");
181  assert(fBkgRatio->GetNbinsX() == nBins &&
182  "MC and data have different numbers of bins.");
183 
184  std::stringstream title;
185  title << std::string(";") << hFirst->GetXaxis()->GetTitle()
186  << std::string(";Events");
187 
188 
189  float axisMin = hFirst->GetBinLowEdge(1);
190  float axisMax = hFirst->GetBinLowEdge(nBins+1);
191 
192  TH1D* hSig = new TH1D("Sig", title.str().c_str(), nBins, axisMin, axisMax);
193  TH1D* hBkg = new TH1D("Bkg", title.str().c_str(), nBins, axisMin, axisMax);
194 
195  for(int i = 1; i <= nBins; ++i)
196  {
197  //These guys are what we're after, signal and background. Initialize them.
198  double nSig = 0;
199  double nBkg = 0;
200  double errSig = 0;
201  double errBkg = 0;
202 
203 
204  // Grab all of the things out of the histograms and store them so that
205  // the math looks a bit nicer.
206  double rSig = fSigRatio->GetBinContent(i);
207  double rSigErr = fSigRatio->GetBinError(i);
208  double rBkg = fBkgRatio->GetBinContent(i);
209  double rBkgErr = fBkgRatio->GetBinError(i);
210  double nFirst = hFirst->GetBinContent(i);
211  double nSecond = hSecond->GetBinContent(i);
212 
213  // d2 is the square of the determinant of the matrix we're inverting.
214  // It shows up in the math a lot, so we'll store it for neatness.
215  double d2 = util::sqr(rSig - rBkg);
216 
217 
218  //We're about to divide by d2 until the cows come home, so make sure it's
219  //not zero. When it is zero, it's probably because those bins were empty.
220  if (d2 != 0)
221  {
222 
223  //Do the math in a way that stores the common terms in the error formula
224  double sSig = (1+rSig);
225  double sBkg = (1+rBkg);
226  nSig = sSig * (nSecond - rBkg*nFirst) / (rSig - rBkg);
227  nBkg = sBkg * (nFirst*rSig - nSecond) / (rSig - rBkg);
228 
229  // Here's the error formula
230  errSig = sqrt((nSecond * util::sqr(sSig))/d2 +
231  (nFirst*util::sqr(sSig) *util::sqr(rBkg))/d2 +
232  util::sqr(rSigErr) * (-((sSig*(nSecond - nFirst*rBkg))/d2) +
233  util::sqr(nSecond - nFirst*rBkg)/(rSig - rBkg)) +
234  util::sqr(-((nFirst*(sSig))/(rSig - rBkg)) +
235  (sSig*(nSecond - nFirst*rBkg))/d2) * util::sqr(rBkgErr));
236  errBkg = sqrt((nSecond * util::sqr(sBkg))/d2 +
237  (nFirst*util::sqr(rSig) * util::sqr(sBkg))/d2 +
238  util::sqr(rSigErr)*util::sqr(-(((-nSecond + nFirst*rSig)*sBkg)/
239  (rSig - rBkg)*2) +
240  (nFirst*sBkg)/(rSig - rBkg)) +
241  util::sqr((-nSecond + nFirst*rSig)/(rSig - rBkg) +
242  ((-nSecond + nFirst*rSig) *sBkg)/d2)*util::sqr(rBkgErr));
243  }
244 
245  // Set the values in the resultant histograms, signal and background.
246  hSig->SetBinContent(i, nSig);
247  hBkg->SetBinContent(i, nBkg);
248  hSig->SetBinError(i, errSig);
249  hBkg->SetBinError(i, errBkg);
250  }
251 
252  fSig = Spectrum(Eigen::ArrayXd(Eigen::Map<Eigen::ArrayXd>(hSig->GetArray(), hSig->GetNbinsX()+2)),
254  fFirst.POT(), fFirst.Livetime());
255  fBkg = Spectrum(Eigen::ArrayXd(Eigen::Map<Eigen::ArrayXd>(hBkg->GetArray(), hBkg->GetNbinsX()+2)),
257  fFirst.POT(), fFirst.Livetime());
258 
259  fDecomposed = true;
260  }
261 
262 
263 
265  {
266 
267 
268  // Check to make sure we haven't already computed ratios
269  // histogram pointers should be null
270  assert(!fSigRatio && "Ratio histograms already constructed.");
271  assert(!fBkgRatio && "Ratio histograms already constructed.");
272 
273  // Pull histograms out of the required Spectrum object
274  TH1D* hFirstSig = fFirstSig .ToTH1(fFirstSig.POT());
275  TH1D* hSecondSig = fSecondSig.ToTH1(fSecondSig.POT());
276  TH1D* hFirstBkg = fFirstBkg .ToTH1(fFirstBkg.POT());
277  TH1D* hSecondBkg = fSecondBkg.ToTH1(fSecondBkg.POT());
278 
279  // Make the axis labels for the ratio histograms out of the
280  std::stringstream title;
281  title << std::string(";") << hFirstSig->GetXaxis()->GetTitle()
282  << std::string(";Ratio");
283 
284  // Get the number of bins from the axis
285  int nBins = hFirstSig->GetNbinsX();
286 
287  // Make sure the histogram binning is sane
288  assert(hFirstSig->GetNbinsX() == nBins &&
289  "Inconsistent binning across spectra.");
290  assert(hSecondSig->GetNbinsX() == nBins &&
291  "Inconsistent binning across spectra.");
292  assert(hFirstBkg->GetNbinsX() == nBins &&
293  "Inconsistent binning across spectra.");
294  assert(hSecondBkg->GetNbinsX() == nBins &&
295  "Inconsistent binning across spectra.");
296 
297  const Binning bins = fAxis->GetBinnings()[0];
298 
299  // Construct the ratio histogram objects
300  fSigRatio = new TH1D("SigRatio", title.str().c_str(), nBins,
301  bins.Min(), bins.Max());
302  fBkgRatio = new TH1D("BkgRatio", title.str().c_str(), nBins,
303  bins.Min(), bins.Max());
304 
305 
306  // Loop over the bins in the histograms to compute the ratios bin-by-bin
307  for (int i = 1; i <= nBins; ++i)
308  {
309  // Ratios for signal
310  if(hFirstSig->GetBinContent(i) != 0 and hSecondSig->GetBinContent(i) != 0)
311  {
312  // if non zero, make an object that does the math
313  RatioError ratioError(hSecondSig->GetBinContent(i),
314  hFirstSig->GetBinContent(i));
315  // Set the bins
316  fSigRatio->SetBinContent(i, ratioError.fRatio);
317  fSigRatio->SetBinError (i, ratioError.fError);
318  }
319  // Repeat for background
320  if(hFirstBkg->GetBinContent(i) != 0 and hSecondBkg->GetBinContent(i) != 0)
321  {
322  RatioError ratioError(hSecondBkg->GetBinContent(i),
323  hFirstBkg->GetBinContent(i));
324  fBkgRatio->SetBinContent(i, ratioError.fRatio);
325  fBkgRatio->SetBinError (i, ratioError.fError);
326  }
327  }
328  }
329 
330 
332  {
333  Decomp();
334 
336 
337  TH1D* hData = data.ToTH1(data.POT());
338  hData->Draw();
339 
340  TH1D* hSig = fSig.ToTH1(data.POT());
341  TH1D* hBkg = fBkg.ToTH1(data.POT());
342 
343  hSig->SetLineColor(kGreen+2);
344  hBkg->SetLineColor(kBlue-4);
345  hSig->SetMarkerColor(kGreen+2);
346  hBkg->SetMarkerColor(kBlue-4);
347 
348  hSig->Draw("same");
349  hBkg->Draw("same");
350 
351 
352  TH1D* hTrueSig = fTrueSig.ToTH1(data.POT());
353  TH1D* hTrueBkg = fTrueBkg.ToTH1(data.POT());
354 
355  hTrueSig->SetLineColor(kGreen+2);
356  hTrueBkg->SetLineColor(kBlue-4);
357 
358  hTrueSig->Draw("hist same");
359  hTrueBkg->Draw("hist same");
360 
361  TH1D* hTrueTotal = (fTrueSig + fTrueBkg).ToTH1(data.POT());
362 
363  hTrueTotal->SetLineColor(kRed);
364  hTrueTotal->Draw("hist same");
365 
366 
367  }
368 
369 
371  {
372  Decomp();
373 
374  TH1D* hFirstSig = fFirstSig.ToTH1(fFirst.POT());
375  TH1D* hSecondSig = fSecondSig.ToTH1(fFirst.POT());
376  TH1D* hFirstBkg = fFirstBkg.ToTH1(fFirst.POT());
377  TH1D* hSecondBkg = fSecondBkg.ToTH1(fFirst.POT());
378 
379 
380  std::vector<TH1D*> hists = {hFirstSig, hSecondSig, hFirstBkg, hSecondBkg};
381 
382  float max = 0;
383  for(const auto& h:hists){
384  float m = h->GetBinContent(h->GetMaximumBin());
385  if(m > max) max = m;
386  }
387 
388  for(const auto& h:hists) h->SetMaximum(1.2*max);
389 
390  hFirstSig->SetLineColor(kRed);
391  hSecondSig->SetLineColor(kBlue-4);
392  hFirstBkg->SetLineColor(kRed);
393  hSecondBkg->SetLineColor(kBlue-4);
394 
395  hFirstBkg->SetLineStyle(2);
396  hSecondBkg->SetLineStyle(2);
397 
398  hFirstSig->Draw("hist");
399  hSecondSig->Draw("hist same");
400  hFirstBkg->Draw("hist same");
401  hSecondBkg->Draw("hist same");
402 
403 
404  }
405 
407  {
408  fSigRatio->Draw("");
409  fBkgRatio->SetLineStyle(2);
410  fBkgRatio->Draw("same");
411  }
412 
413 
414 
415  RatioError::RatioError(double numerator, double denominator):
416  fRatio(numerator/denominator), // Ratio is the ratio
417  // Error for gaussian counting errors
418  fError(fRatio * sqrt(1/ numerator + 1/denominator))
419  { }
420 
421 }
422 
423 
424 
Spectrum fSig
Decomposed spectrum of signal events.
Spectrum fBkg
Decomposed spectrum of signal events.
const std::vector< Binning > & GetBinnings() const
Definition: Spectrum.h:268
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
virtual Spectrum NCComponent() const override
unsigned int NDimensions() const
Definition: HistAxis.h:87
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:209
const Var weight
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Spectrum fTrueSig
True signal spectrum.
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH1D * fBkgRatio
Ratio of second/first spectra for bkg.
Spectrum fFirstSig
Spectrum for first sample, signal.
Float_t tmp
Definition: plot.C:36
virtual Spectrum NCAntiComponent() const override
TString hists[nhists]
Definition: bdt_com.C:3
bool fDecomposed
Indicates whether decomposition was done.
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Class for data driven decomposition of signal and background using two different samples with slightl...
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
const char * label
const XML_Char const XML_Char * data
Definition: expat.h:268
Given a numerator and denominator, this class calculates the ratio and associated error on that ratio...
Spectrum fSecondSig
Spectrum for second sample, signal.
TH1D * fSigRatio
Ratio of second/first spectra for signal.
const SpectrumLoaderBase & fLoaderData
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:607
void Decomp() const
Preform the decomposition using the ratios in the MC sample.
Spectrum fFirst
Spectrum of events that pass the firstSample cut.
const HistAxis * fAxis
Axis used for spectra.
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:578
Spectrum fTrueBkg
True backgroundh spectrum.
std::vector< float > Spectrum
Definition: Constants.h:527
double POT() const
Definition: Spectrum.h:231
OStream cout
Definition: OStream.cxx:6
Base class for the various types of spectrum loader.
const Binning bins
Definition: NumuCC_CPiBin.h:8
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const SpectrumLoaderBase & fLoaderMC
TDirectory * dir
Definition: macro.C:5
static std::unique_ptr< TwoSampleDecomp > LoadFrom(TDirectory *dir, const std::string &name)
RatioError(double numerator, double denominator)
GenericHistAxis< Var > HistAxis
Definition: HistAxis.h:111
assert(nhit_max >=nhit_nbins)
void RatioCalc() const
Calculate the ratios for the signal and background, second/first.
Spectrum fSecond
Spectrum of events that fail the firstSample cut.
void SaveTo(TDirectory *dir, const std::string &name) const override
Spectrum fFirstBkg
Spectrum for first sample, background.
const std::vector< std::string > & GetLabels() const
Definition: Spectrum.h:267
Spectrum fSecondBkg
Spectrum for second sample, background.
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
virtual bool Gone() const
Indicate whether or not Go has been called.
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:234
static NullLoader kNullLoader
Dummy loader that doesn&#39;t load any files.
const std::vector< Binning > & GetBinnings() const
Definition: HistAxis.h:91