15 #include "TDirectory.h" 16 #include "TObjString.h" 17 #include "Utilities/func/MathUtil.h" 30 const Cut& firstSample,
37 firstSample, isSignal, anaCut, shift, weight)
46 const Cut& firstSample,
56 fFirstSig (loaderMC, axis, firstSample && isSignal && anaCut,
58 fSecondSig(loaderMC, axis, !firstSample && isSignal && anaCut,
60 fFirstBkg (loaderMC, axis, firstSample && !isSignal && anaCut,
62 fSecondBkg(loaderMC, axis, !firstSample && !isSignal && anaCut,
64 fTrueSig (loaderMC, axis, isSignal && anaCut,
66 fTrueBkg (loaderMC, axis, !isSignal && anaCut,
68 fFirst (loaderData, axis, firstSample && anaCut,
70 fSecond (loaderData, axis, !firstSample && anaCut,
79 std::cout <<
"Error: TwoSampleDecomp not equipped to handle multi-dimensional axis." <<
std::endl;
108 TDirectory*
tmp = gDirectory;
110 dir = dir->mkdir(name.c_str());
113 TObjString(
"TwoSampleDecomp").Write(
"type");
135 dir = dir->GetDirectory(name.c_str());
164 &&
"Call Go() on MC SpectrumLoader before TwoSampleDecomp::Decomp");
166 &&
"Call Go() on Data SpectrumLoader before TwoSampleDecomp::Decomp");
174 int nBins = hFirst->GetNbinsX();
177 assert(hFirst->GetNbinsX() == nBins &&
178 "Inconsistent binning across spectra.");
179 assert(hSecond->GetNbinsX() == nBins &&
180 "Inconsistent binning across spectra.");
182 "MC and data have different numbers of bins.");
184 "MC and data have different numbers of bins.");
186 std::stringstream
title;
187 title <<
std::string(
";") << hFirst->GetXaxis()->GetTitle()
191 float axisMin = hFirst->GetBinLowEdge(1);
192 float axisMax = hFirst->GetBinLowEdge(nBins+1);
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);
212 double nFirst = hFirst->GetBinContent(
i);
213 double nSecond = hSecond->GetBinContent(
i);
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);
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));
242 (nFirst*sBkg)/(rSig - rBkg)) +
243 util::sqr((-nSecond + nFirst*rSig)/(rSig - rBkg) +
244 ((-nSecond + nFirst*rSig) *sBkg)/d2)*
util::sqr(rBkgErr));
248 hSig->SetBinContent(
i, nSig);
249 hBkg->SetBinContent(
i, nBkg);
250 hSig->SetBinError(
i, errSig);
251 hBkg->SetBinError(
i, errBkg);
254 fSig =
Spectrum(Eigen::ArrayXd(Eigen::Map<Eigen::ArrayXd>(hSig->GetArray(), hSig->GetNbinsX()+2)),
257 fBkg =
Spectrum(Eigen::ArrayXd(Eigen::Map<Eigen::ArrayXd>(hBkg->GetArray(), hBkg->GetNbinsX()+2)),
282 std::stringstream
title;
283 title <<
std::string(
";") << hFirstSig->GetXaxis()->GetTitle()
287 int nBins = hFirstSig->GetNbinsX();
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.");
303 bins.Min(), bins.Max());
305 bins.Min(), bins.Max());
312 if(hFirstSig->GetBinContent(
i) != 0 and hSecondSig->GetBinContent(
i) != 0)
315 RatioError ratioError(hSecondSig->GetBinContent(
i),
316 hFirstSig->GetBinContent(
i));
318 fSigRatio->SetBinContent(
i, ratioError.fRatio);
319 fSigRatio->SetBinError (
i, ratioError.fError);
322 if(hFirstBkg->GetBinContent(
i) != 0 and hSecondBkg->GetBinContent(
i) != 0)
324 RatioError ratioError(hSecondBkg->GetBinContent(
i),
325 hFirstBkg->GetBinContent(
i));
326 fBkgRatio->SetBinContent(
i, ratioError.fRatio);
327 fBkgRatio->SetBinError (
i, ratioError.fError);
339 TH1D* hData = data.
ToTH1(data.
POT());
345 hSig->SetLineColor(
kGreen+2);
346 hBkg->SetLineColor(
kBlue-4);
347 hSig->SetMarkerColor(
kGreen+2);
348 hBkg->SetMarkerColor(
kBlue-4);
357 hTrueSig->SetLineColor(
kGreen+2);
358 hTrueBkg->SetLineColor(
kBlue-4);
360 hTrueSig->Draw(
"hist same");
361 hTrueBkg->Draw(
"hist same");
365 hTrueTotal->SetLineColor(
kRed);
366 hTrueTotal->Draw(
"hist same");
382 std::vector<TH1D*>
hists = {hFirstSig, hSecondSig, hFirstBkg, hSecondBkg};
385 for(
const auto&
h:hists){
386 float m =
h->GetBinContent(
h->GetMaximumBin());
390 for(
const auto&
h:hists)
h->SetMaximum(1.2*max);
392 hFirstSig->SetLineColor(
kRed);
393 hSecondSig->SetLineColor(
kBlue-4);
394 hFirstBkg->SetLineColor(
kRed);
395 hSecondBkg->SetLineColor(
kBlue-4);
397 hFirstBkg->SetLineStyle(2);
398 hSecondBkg->SetLineStyle(2);
400 hFirstSig->Draw(
"hist");
401 hSecondSig->Draw(
"hist same");
402 hFirstBkg->Draw(
"hist same");
403 hSecondBkg->Draw(
"hist same");
418 fRatio(numerator/denominator),
420 fError(fRatio *
sqrt(1/ numerator + 1/denominator))
Spectrum fSig
Decomposed spectrum of signal events.
Spectrum fBkg
Decomposed spectrum of signal events.
const std::vector< Binning > & GetBinnings() const
_HistAxis< Var > HistAxis
Represent the binning of a Spectrum's x-axis.
Cuts and Vars for the 2020 FD DiF Study.
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.
Simple record of shifts applied to systematic parameters.
Spectrum fTrueSig
True signal spectrum.
TH1D * fBkgRatio
Ratio of second/first spectra for bkg.
Spectrum fFirstSig
Spectrum for first sample, signal.
virtual Spectrum NCAntiComponent() const override
bool fDecomposed
Indicates whether decomposition was done.
T sqr(T x)
More efficient square function than pow(x,2)
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.
const XML_Char const XML_Char * data
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)
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
Spectrum fTrueBkg
True backgroundh spectrum.
std::vector< float > Spectrum
Base class for the various types of spectrum loader.
const std::vector< Binning > & GetBinnings() const
unsigned int NDimensions() const
const SpectrumLoaderBase & fLoaderMC
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 DrawTwoSamplesWithRatios()
void SaveTo(TDirectory *dir, const std::string &name) const override
Spectrum fFirstBkg
Spectrum for first sample, background.
const std::vector< std::string > & GetLabels() const
Spectrum fSecondBkg
Spectrum for second sample, background.
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
virtual bool Gone() const
Indicate whether or not Go has been called.
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
static NullLoader kNullLoader
Dummy loader that doesn't load any files.