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