Skip to content

Commit 68d5234

Browse files
pstahlhulubynets
pstahlhu
and
lubynets
authored
[PWGHF] Improvements to HFInvMassFitter (#13002)
Co-authored-by: Oleksii Lubynets <O.Lubynets@gsi.de>
1 parent 3bf89aa commit 68d5234

File tree

4 files changed

+327
-157
lines changed

4 files changed

+327
-157
lines changed

‎PWGHF/D2H/Macros/HFInvMassFitter.cxx‎

Lines changed: 142 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717
/// \author Xinye Peng <xinye.peng@cern.ch>
1818
/// \author Biao Zhang <biao.zhang@cern.ch>
1919
/// \author Oleksii Lubynets <oleksii.lubynets@cern.ch>
20+
/// \author Phil Stahlhut <phil.lennart.stahlhut@cern.ch>
2021

2122
#include "HFInvMassFitter.h"
2223

@@ -50,6 +51,7 @@
5051
#include <cstring>
5152
#include <stdexcept>
5253
#include <string>
54+
#include <vector>
5355

5456
using namespace RooFit;
5557

@@ -74,7 +76,7 @@ HFInvMassFitter::HFInvMassFitter() : mHistoInvMass(nullptr),
7476
mNSigmaForSidebands(4.),
7577
mNSigmaForSgn(3.),
7678
mSigmaSgnErr(0.),
77-
mSigmaSgnDoubleGaus(0.012),
79+
mSigmaSgnDoubleGaus(0.025),
7880
mFixedMean(kFALSE),
7981
mBoundMean(kFALSE),
8082
mBoundReflMean(kFALSE),
@@ -103,6 +105,8 @@ HFInvMassFitter::HFInvMassFitter() : mHistoInvMass(nullptr),
103105
mFixReflOverSgn(kFALSE),
104106
mRooMeanSgn(nullptr),
105107
mRooSigmaSgn(nullptr),
108+
mRooSecSigmaSgn(nullptr),
109+
mRooFracDoubleGaus(nullptr),
106110
mSgnPdf(nullptr),
107111
mBkgPdf(nullptr),
108112
mReflPdf(nullptr),
@@ -145,7 +149,7 @@ HFInvMassFitter::HFInvMassFitter(const TH1* histoToFit, Double_t minValue, Doubl
145149
mNSigmaForSidebands(3.),
146150
mNSigmaForSgn(3.),
147151
mSigmaSgnErr(0.),
148-
mSigmaSgnDoubleGaus(0.012),
152+
mSigmaSgnDoubleGaus(0.025),
149153
mFixedMean(kFALSE),
150154
mBoundMean(kFALSE),
151155
mBoundReflMean(kFALSE),
@@ -174,6 +178,8 @@ HFInvMassFitter::HFInvMassFitter(const TH1* histoToFit, Double_t minValue, Doubl
174178
mFixReflOverSgn(kFALSE),
175179
mRooMeanSgn(nullptr),
176180
mRooSigmaSgn(nullptr),
181+
mRooSecSigmaSgn(nullptr),
182+
mRooFracDoubleGaus(nullptr),
177183
mSgnPdf(nullptr),
178184
mBkgPdf(nullptr),
179185
mReflPdf(nullptr),
@@ -207,6 +213,8 @@ HFInvMassFitter::~HFInvMassFitter()
207213
delete mHistoTemplateRefl;
208214
delete mRooMeanSgn;
209215
delete mRooSigmaSgn;
216+
delete mRooSecSigmaSgn;
217+
delete mRooFracDoubleGaus;
210218
delete mSgnPdf;
211219
delete mBkgPdf;
212220
delete mReflPdf;
@@ -356,13 +364,14 @@ void HFInvMassFitter::doFit()
356364
mTotalPdf->plotOn(mInvMassFrame, Name("Tot_c"), LineColor(kBlue));
357365
mSgnPdf->plotOn(mInvMassFrame, Normalization(1.0, RooAbsReal::RelativeExpected), DrawOption("F"), FillColor(TColor::GetColorTransparent(kBlue, 0.2)), VLines());
358366
mChiSquareOverNdfTotal = mInvMassFrame->chiSquare("Tot_c", "data_c"); // calculate reduced chi2 / DNF
367+
359368
// plot residual distribution
360369
mResidualFrame = mass->frame(Title("Residual Distribution"));
361370
RooHist* residualHistogram = mInvMassFrame->residHist("data_c", "Bkg_c");
362371
mResidualFrame->addPlotable(residualHistogram, "P");
363372
mSgnPdf->plotOn(mResidualFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineColor(kBlue));
364373
}
365-
mass->setRange("bkgForSignificance", mRooMeanSgn->getVal() - mNSigmaForSgn * mRooSigmaSgn->getVal(), mRooMeanSgn->getVal() + mNSigmaForSgn * mRooSigmaSgn->getVal());
374+
mass->setRange("bkgForSignificance", mRooMeanSgn->getVal() - mNSigmaForSgn * mRooSecSigmaSgn->getVal(), mRooMeanSgn->getVal() + mNSigmaForSgn * mRooSecSigmaSgn->getVal());
366375
bkgIntegral = mBkgPdf->createIntegral(*mass, NormSet(*mass), Range("bkgForSignificance"));
367376
mIntegralBkg = bkgIntegral->getValV();
368377
calculateBackground(mBkgYield, mBkgYieldErr);
@@ -372,6 +381,9 @@ void HFInvMassFitter::doFit()
372381
calculateSignal(mRawYield, mRawYieldErr);
373382
countSignal(mRawYieldCounted, mRawYieldCountedErr);
374383
calculateSignificance(mSignificance, mSignificanceErr);
384+
// Fit to data ratio
385+
mRatioFrame = mass->frame(Title(Form("%s", mHistoInvMass->GetTitle())));
386+
calculateFitToDataRatio();
375387
}
376388
}
377389

@@ -440,10 +452,10 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const
440452
workspace.import(*sgnFuncGaus);
441453
delete sgnFuncGaus;
442454
// signal double Gaussian
443-
RooRealVar sigmaDoubleGaus("sigmaDoubleGaus", "sigma2Gaus", mSigmaSgn, mSigmaSgn - 0.01, mSigmaSgn + 0.01);
455+
RooRealVar sigmaDoubleGaus("sigmaDoubleGaus", "sigma2Gaus", mSigmaSgnDoubleGaus, mSigmaSgnDoubleGaus - 0.003, mSigmaSgnDoubleGaus + 0.003);
444456
if (mBoundSigma) {
445-
sigmaDoubleGaus.setMax(mSigmaSgn * (1 + mParamSgn));
446-
sigmaDoubleGaus.setMin(mSigmaSgn * (1 - mParamSgn));
457+
sigmaDoubleGaus.setMax(mSigmaSgnDoubleGaus * (1 + mParamSgn));
458+
sigmaDoubleGaus.setMin(mSigmaSgnDoubleGaus * (1 - mParamSgn));
447459
}
448460
if (mFixedSigma) {
449461
sigma.setVal(mSigmaSgn);
@@ -552,64 +564,81 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const
552564
workspace.import(*reflFuncPoly6);
553565
delete reflFuncPoly6;
554566
}
567+
555568
// draw fit output
556-
void HFInvMassFitter::drawFit(TVirtualPad* pad, Int_t writeFitInfo)
569+
void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector<std::string>& plotLabels, Bool_t writeParInfo)
557570
{
558571
gStyle->SetOptStat(0);
559572
gStyle->SetCanvasColor(0);
560573
gStyle->SetFrameFillColor(0);
561574
pad->cd();
562-
if (writeFitInfo > 0) {
563-
auto* textInfoLeft = new TPaveText(0.12, 0.65, 0.47, 0.89, "NDC");
564-
auto* textInfoRight = new TPaveText(0.6, 0.7, 1., .87, "NDC");
565-
textInfoLeft->SetBorderSize(0);
566-
textInfoLeft->SetFillStyle(0);
567-
textInfoRight->SetBorderSize(0);
568-
textInfoRight->SetFillStyle(0);
569-
textInfoRight->SetTextColor(kBlue);
570-
textInfoLeft->AddText(Form("S = %.0f #pm %.0f ", mRawYield, mRawYieldErr));
571-
textInfoLeft->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr));
572-
if (mTypeOfBkgPdf != NoBkg) {
573-
textInfoLeft->AddText(Form("B (%d#sigma) = %.0f #pm %.0f", mNSigmaForSidebands, mBkgYield, mBkgYieldErr));
574-
textInfoLeft->AddText(Form("S/B (%d#sigma) = %.4g ", mNSigmaForSidebands, mRawYield / mBkgYield));
575-
}
576-
if (mReflPdf != nullptr) {
577-
textInfoLeft->AddText(Form("Refl/Sig = %.3f #pm %.3f ", mReflOverSgn, 0.0));
578-
}
579-
if (mTypeOfBkgPdf != NoBkg) {
580-
textInfoLeft->AddText(Form("Signif (%d#sigma) = %.1f #pm %.1f ", mNSigmaForSidebands, mSignificance, mSignificanceErr));
581-
textInfoLeft->AddText(Form("#chi^{2} / ndf = %.3f", mChiSquareOverNdfTotal));
582-
}
575+
// Fit metrics
576+
auto* textFitMetrics = new TPaveText(0.65, 0.7, 0.9, 0.88, "NDC");
577+
textFitMetrics->SetBorderSize(0);
578+
textFitMetrics->SetFillStyle(0);
579+
textFitMetrics->SetTextSize(0.04);
580+
textFitMetrics->SetTextAlign(33);
581+
textFitMetrics->AddText(Form("S = %.0f #pm %.0f ", mRawYield, mRawYieldErr));
582+
textFitMetrics->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr));
583+
if (mTypeOfBkgPdf != NoBkg) {
584+
textFitMetrics->AddText(Form("B (%d#sigma) = %.0f #pm %.0f", mNSigmaForSidebands, mBkgYield, mBkgYieldErr));
585+
textFitMetrics->AddText(Form("S/B (%d#sigma) = %.4g ", mNSigmaForSidebands, mRawYield / mBkgYield));
586+
textFitMetrics->AddText(Form("Significance (%d#sigma) = %.1f #pm %.1f ", mNSigmaForSidebands, mSignificance, mSignificanceErr));
587+
textFitMetrics->AddText(Form("#chi^{2} / ndf = %.3f", mChiSquareOverNdfTotal));
588+
}
589+
if (mReflPdf != nullptr) {
590+
textFitMetrics->AddText(Form("Refl/Sig = %.3f #pm %.3f ", mReflOverSgn, 0.0));
591+
}
592+
mInvMassFrame->addObject(textFitMetrics);
593+
// Analysis information
594+
auto* textAnalysisInfo = new TPaveText(0.18, 0.78, 0.35, 0.88, "NDC");
595+
textAnalysisInfo->SetBorderSize(0);
596+
textAnalysisInfo->SetFillStyle(0);
597+
textAnalysisInfo->SetTextSize(0.05);
598+
textAnalysisInfo->SetTextAlign(13);
599+
for (const auto& label : plotLabels) {
600+
textAnalysisInfo->AddText(label.c_str());
601+
}
602+
mInvMassFrame->addObject(textAnalysisInfo);
603+
if (writeParInfo) {
604+
// right text box
605+
auto* textSignalPar = new TPaveText(0.18, 0.65, 0.4, 0.75, "NDC");
606+
textSignalPar->SetBorderSize(0);
607+
textSignalPar->SetFillStyle(0);
608+
textSignalPar->SetTextColor(kBlue);
609+
textSignalPar->SetTextAlign(13);
583610
if (mFixedMean) {
584-
textInfoRight->AddText(Form("mean(fixed) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
611+
textSignalPar->AddText(Form("mean(fixed) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
585612
} else {
586-
textInfoRight->AddText(Form("mean(free) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
613+
textSignalPar->AddText(Form("mean(free) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
587614
}
588615
if (mTypeOfSgnPdf == DoubleGaus) {
589-
auto const& baseSigmaSgn = mWorkspace->var("sigma");
616+
if (mFixedSigma) {
617+
textSignalPar->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
618+
} else {
619+
textSignalPar->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
620+
}
590621
if (mFixedSigmaDoubleGaus) {
591-
textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
592-
textInfoRight->AddText(Form("sigma 2(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
622+
textSignalPar->AddText(Form("sigma 2(fixed) = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError()));
593623
} else {
594-
textInfoRight->AddText(Form("sigma(free) = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
595-
textInfoRight->AddText(Form("sigma 2(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
624+
textSignalPar->AddText(Form("sigma 2(free) = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError()));
596625
}
597626
} else if (mFixedSigma) {
598-
textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
627+
textSignalPar->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
599628
} else {
600-
textInfoRight->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
601-
}
602-
mInvMassFrame->addObject(textInfoLeft);
603-
mInvMassFrame->addObject(textInfoRight);
604-
mInvMassFrame->GetYaxis()->SetTitleOffset(1.8);
605-
gPad->SetLeftMargin(0.15);
606-
mInvMassFrame->GetYaxis()->SetTitle(Form("%s", mHistoInvMass->GetYaxis()->GetTitle()));
607-
mInvMassFrame->GetXaxis()->SetTitle(Form("%s", mHistoInvMass->GetXaxis()->GetTitle()));
608-
mInvMassFrame->Draw();
609-
highlightPeakRegion(mInvMassFrame);
610-
if (mHistoTemplateRefl != nullptr) {
611-
mReflFrame->Draw("same");
629+
textSignalPar->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
612630
}
631+
mInvMassFrame->addObject(textSignalPar);
632+
}
633+
mInvMassFrame->GetXaxis()->SetTitleOffset(1.2);
634+
mInvMassFrame->GetYaxis()->SetTitleOffset(1.8);
635+
gPad->SetLeftMargin(0.15);
636+
mInvMassFrame->GetYaxis()->SetTitle(Form("%s", mHistoInvMass->GetYaxis()->GetTitle()));
637+
mInvMassFrame->GetXaxis()->SetTitle(Form("%s", mHistoInvMass->GetXaxis()->GetTitle()));
638+
mInvMassFrame->Draw();
639+
highlightPeakRegion(mInvMassFrame);
640+
if (mHistoTemplateRefl) {
641+
mReflFrame->Draw("same");
613642
}
614643
}
615644

@@ -626,9 +655,8 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad)
626655
textInfo->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr));
627656
textInfo->AddText(Form("mean = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
628657
if (mTypeOfSgnPdf == DoubleGaus) {
629-
auto const& baseSigmaSgn = mWorkspace->var("sigma");
630-
textInfo->AddText(Form("sigma = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
631-
textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
658+
textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
659+
textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError()));
632660
} else {
633661
textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
634662
}
@@ -637,6 +665,24 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad)
637665
highlightPeakRegion(mResidualFrame);
638666
}
639667

668+
// draw ratio on canvas
669+
void HFInvMassFitter::drawRatio(TVirtualPad* pad)
670+
{
671+
pad->cd();
672+
mRatioFrame->GetXaxis()->SetTitleOffset(1.2);
673+
mRatioFrame->GetYaxis()->SetTitleOffset(1.5);
674+
mRatioFrame->GetYaxis()->SetTitle("Fit / Data");
675+
double xMin = mRatioFrame->GetXaxis()->GetXmin();
676+
double xMax = mRatioFrame->GetXaxis()->GetXmax();
677+
auto* line = new TLine(xMin, 1.0, xMax, 1.0);
678+
line->SetLineColor(kGray);
679+
line->SetLineStyle(2);
680+
line->SetLineWidth(2);
681+
mRatioFrame->addObject(line);
682+
mRatioFrame->Draw();
683+
highlightPeakRegion(mRatioFrame);
684+
}
685+
640686
// draw peak region with vertical lines
641687
void HFInvMassFitter::highlightPeakRegion(const RooPlot* plot, Color_t color, Width_t width, Style_t style) const
642688
{
@@ -646,7 +692,7 @@ void HFInvMassFitter::highlightPeakRegion(const RooPlot* plot, Color_t color, Wi
646692
double const yMin = plot->GetMinimum();
647693
double const yMax = plot->GetMaximum();
648694
const Double_t mean = mRooMeanSgn->getVal();
649-
const Double_t sigma = mRooSigmaSgn->getVal();
695+
const Double_t sigma = mRooSecSigmaSgn->getVal();
650696
const Double_t minForSgn = mean - mNSigmaForSidebands * sigma;
651697
const Double_t maxForSgn = mean + mNSigmaForSidebands * sigma;
652698
auto* leftLine = new TLine(minForSgn, yMin, minForSgn, yMax);
@@ -671,7 +717,7 @@ void HFInvMassFitter::drawReflection(TVirtualPad* pad)
671717
void HFInvMassFitter::countSignal(Double_t& signal, Double_t& signalErr) const
672718
{
673719
const Double_t mean = mRooMeanSgn->getVal();
674-
const Double_t sigma = mRooSigmaSgn->getVal();
720+
const Double_t sigma = mRooSecSigmaSgn->getVal();
675721
const Double_t minForSgn = mean - mNSigmaForSidebands * sigma;
676722
const Double_t maxForSgn = mean + mNSigmaForSidebands * sigma;
677723
const Int_t binForMinSgn = mHistoInvMass->FindBin(minForSgn);
@@ -763,21 +809,27 @@ RooAbsPdf* HFInvMassFitter::createSignalFitFunction(RooWorkspace* workspace)
763809
case 0: {
764810
sgnPdf = workspace->pdf("sgnFuncGaus");
765811
mRooSigmaSgn = workspace->var("sigma");
812+
mRooSecSigmaSgn = workspace->var("sigma");
766813
mRooMeanSgn = workspace->var("mean");
767814
} break;
768815
case 1: {
769816
sgnPdf = workspace->pdf("sgnFuncDoubleGaus");
770-
mRooSigmaSgn = workspace->var("sigmaDoubleGaus");
817+
mRooSigmaSgn = workspace->var("sigma");
818+
mRooSecSigmaSgn = workspace->var("sigmaDoubleGaus");
771819
mRooMeanSgn = workspace->var("mean");
820+
mRooFracDoubleGaus = workspace->var("fracDoubleGaus");
772821
} break;
773822
case 2: {
774823
sgnPdf = workspace->pdf("sgnFuncGausRatio");
775-
mRooSigmaSgn = workspace->var("sigmaDoubleGausRatio");
824+
mRooSigmaSgn = workspace->var("sigma");
825+
mRooSecSigmaSgn = workspace->var("sigmaDoubleGausRatio");
776826
mRooMeanSgn = workspace->var("mean");
827+
mRooFracDoubleGaus = workspace->var("fracDoubleGausRatio");
777828
} break;
778829
case 3: {
779830
sgnPdf = workspace->pdf("sgnFuncDoublePeak");
780-
mRooSigmaSgn = workspace->var("sigmaSec");
831+
mRooSigmaSgn = workspace->var("sigma");
832+
mRooSecSigmaSgn = workspace->var("sigmaSec");
781833
mRooMeanSgn = workspace->var("meanSec");
782834
} break;
783835
default:
@@ -818,6 +870,40 @@ void HFInvMassFitter::plotRefl(RooAbsPdf* pdf)
818870
pdf->plotOn(mInvMassFrame, Components(namesOfReflPdf.at(mTypeOfReflPdf).c_str()), Name("Refl_c"), LineColor(kGreen));
819871
}
820872

873+
// Calculate fit to data ratio
874+
void HFInvMassFitter::calculateFitToDataRatio() const
875+
{
876+
if (!mInvMassFrame)
877+
return;
878+
879+
// Get the data and fit curves from the frame
880+
auto* dataHist = dynamic_cast<RooHist*>(mInvMassFrame->findObject("data_c"));
881+
auto* fitCurve = dynamic_cast<RooCurve*>(mInvMassFrame->findObject("Tot_c")); // or the relevant fit curve
882+
883+
if (!dataHist || !fitCurve)
884+
return;
885+
886+
RooHist* ratioHist = new RooHist();
887+
888+
for (int i = 0; i < dataHist->GetN(); ++i) {
889+
double x, dataY, dataErr;
890+
dataHist->GetPoint(i, x, dataY);
891+
dataErr = dataHist->GetErrorY(i);
892+
893+
double fitY = fitCurve->Eval(x);
894+
895+
double ratio = dataY != 0 ? fitY / dataY : 0;
896+
double err = dataY != 0 ? ratio * dataErr / dataY : 0;
897+
898+
ratioHist->SetPoint(i, x, ratio);
899+
ratioHist->SetPointError(i, 0, 0, err, err);
900+
}
901+
902+
mRatioFrame->addPlotable(ratioHist, "P");
903+
mRatioFrame->SetMinimum(0.5);
904+
mRatioFrame->SetMaximum(1.5);
905+
}
906+
821907
// Fix reflection pdf
822908
void HFInvMassFitter::setReflFuncFixed()
823909
{

0 commit comments

Comments
 (0)