Skip to content

[PWGHF] Improvements to HFInvMassFitter #13002

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
zhangbiao-phy merged 12 commits into from
Nov 17, 2025
+327 −157
Merged

[PWGHF] Improvements to HFInvMassFitter #13002

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
19475fa
HFInvMassFitter: Fix canvas division, add ratio histograms
pstahlhu Aug 26, 2025
5488306
Improve handling of second sigma for double Gaussian fits
pstahlhu Sep 14, 2025
b464a70
Store fraction for double Gaussian and add possibility to fix it
pstahlhu Sep 16, 2025
92823fc
Store signal parameters to hist only when present
pstahlhu Sep 16, 2025
56db18f
Improve plotting style of the mass canvases
pstahlhu Sep 19, 2025
0db0b7f
Implement Oleksii's comments
pstahlhu Sep 19, 2025
4a19c46
adapt config_massfitter for FixFracDoubleGaus
lubynets Sep 18, 2025
be724be
Fix includes in runMassFitter.C
pstahlhu Sep 19, 2025
02f1bba
Restore S_count metric to mass fit histograms
pstahlhu Nov 7, 2025
e0cd945
Comply with Vit's clang-tidy fixes
pstahlhu Nov 13, 2025
f61f7be
Unify histogram order and fix if condition for plot labels
pstahlhu Nov 13, 2025
6f24caa
Fix lambda to fix fraction of double Gaussian
pstahlhu Nov 14, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
198 changes: 142 additions & 56 deletions PWGHF/D2H/Macros/HFInvMassFitter.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.

Check failure on line 1 in PWGHF/D2H/Macros/HFInvMassFitter.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[name/workflow-file] 

Name of a workflow file must match the name of the main struct in it (without the PWG prefix). (Class implementation files should be in "Core" directories.)
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
Expand All @@ -17,6 +17,7 @@
/// \author Xinye Peng <[email protected]>
/// \author Biao Zhang <[email protected]>
/// \author Oleksii Lubynets <[email protected]>
/// \author Phil Stahlhut <[email protected]>

#include "HFInvMassFitter.h"

Expand Down Expand Up @@ -50,6 +51,7 @@
#include <cstring>
#include <stdexcept>
#include <string>
#include <vector>

using namespace RooFit;

Expand All @@ -74,7 +76,7 @@
mNSigmaForSidebands(4.),
mNSigmaForSgn(3.),
mSigmaSgnErr(0.),
mSigmaSgnDoubleGaus(0.012),
mSigmaSgnDoubleGaus(0.025),
mFixedMean(kFALSE),
mBoundMean(kFALSE),
mBoundReflMean(kFALSE),
Expand Down Expand Up @@ -103,6 +105,8 @@
mFixReflOverSgn(kFALSE),
mRooMeanSgn(nullptr),
mRooSigmaSgn(nullptr),
mRooSecSigmaSgn(nullptr),
mRooFracDoubleGaus(nullptr),
mSgnPdf(nullptr),
mBkgPdf(nullptr),
mReflPdf(nullptr),
Expand Down Expand Up @@ -145,7 +149,7 @@
mNSigmaForSidebands(3.),
mNSigmaForSgn(3.),
mSigmaSgnErr(0.),
mSigmaSgnDoubleGaus(0.012),
mSigmaSgnDoubleGaus(0.025),
mFixedMean(kFALSE),
mBoundMean(kFALSE),
mBoundReflMean(kFALSE),
Expand Down Expand Up @@ -174,6 +178,8 @@
mFixReflOverSgn(kFALSE),
mRooMeanSgn(nullptr),
mRooSigmaSgn(nullptr),
mRooSecSigmaSgn(nullptr),
mRooFracDoubleGaus(nullptr),
mSgnPdf(nullptr),
mBkgPdf(nullptr),
mReflPdf(nullptr),
Expand Down Expand Up @@ -207,6 +213,8 @@
delete mHistoTemplateRefl;
delete mRooMeanSgn;
delete mRooSigmaSgn;
delete mRooSecSigmaSgn;
delete mRooFracDoubleGaus;
delete mSgnPdf;
delete mBkgPdf;
delete mReflPdf;
Expand Down Expand Up @@ -356,13 +364,14 @@
mTotalPdf->plotOn(mInvMassFrame, Name("Tot_c"), LineColor(kBlue));
mSgnPdf->plotOn(mInvMassFrame, Normalization(1.0, RooAbsReal::RelativeExpected), DrawOption("F"), FillColor(TColor::GetColorTransparent(kBlue, 0.2)), VLines());
mChiSquareOverNdfTotal = mInvMassFrame->chiSquare("Tot_c", "data_c"); // calculate reduced chi2 / DNF

// plot residual distribution
mResidualFrame = mass->frame(Title("Residual Distribution"));
RooHist* residualHistogram = mInvMassFrame->residHist("data_c", "Bkg_c");
mResidualFrame->addPlotable(residualHistogram, "P");
mSgnPdf->plotOn(mResidualFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineColor(kBlue));
}
mass->setRange("bkgForSignificance", mRooMeanSgn->getVal() - mNSigmaForSgn * mRooSigmaSgn->getVal(), mRooMeanSgn->getVal() + mNSigmaForSgn * mRooSigmaSgn->getVal());
mass->setRange("bkgForSignificance", mRooMeanSgn->getVal() - mNSigmaForSgn * mRooSecSigmaSgn->getVal(), mRooMeanSgn->getVal() + mNSigmaForSgn * mRooSecSigmaSgn->getVal());
bkgIntegral = mBkgPdf->createIntegral(*mass, NormSet(*mass), Range("bkgForSignificance"));
mIntegralBkg = bkgIntegral->getValV();
calculateBackground(mBkgYield, mBkgYieldErr);
Expand All @@ -372,6 +381,9 @@
calculateSignal(mRawYield, mRawYieldErr);
countSignal(mRawYieldCounted, mRawYieldCountedErr);
calculateSignificance(mSignificance, mSignificanceErr);
// Fit to data ratio
mRatioFrame = mass->frame(Title(Form("%s", mHistoInvMass->GetTitle())));
calculateFitToDataRatio();
}
}

Expand Down Expand Up @@ -440,10 +452,10 @@
workspace.import(*sgnFuncGaus);
delete sgnFuncGaus;
// signal double Gaussian
RooRealVar sigmaDoubleGaus("sigmaDoubleGaus", "sigma2Gaus", mSigmaSgn, mSigmaSgn - 0.01, mSigmaSgn + 0.01);
RooRealVar sigmaDoubleGaus("sigmaDoubleGaus", "sigma2Gaus", mSigmaSgnDoubleGaus, mSigmaSgnDoubleGaus - 0.003, mSigmaSgnDoubleGaus + 0.003);
if (mBoundSigma) {
sigmaDoubleGaus.setMax(mSigmaSgn * (1 + mParamSgn));
sigmaDoubleGaus.setMin(mSigmaSgn * (1 - mParamSgn));
sigmaDoubleGaus.setMax(mSigmaSgnDoubleGaus * (1 + mParamSgn));
sigmaDoubleGaus.setMin(mSigmaSgnDoubleGaus * (1 - mParamSgn));
}
if (mFixedSigma) {
sigma.setVal(mSigmaSgn);
Expand Down Expand Up @@ -552,64 +564,81 @@
workspace.import(*reflFuncPoly6);
delete reflFuncPoly6;
}

// draw fit output
void HFInvMassFitter::drawFit(TVirtualPad* pad, Int_t writeFitInfo)
void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector<std::string>& plotLabels, Bool_t writeParInfo)
{
gStyle->SetOptStat(0);
gStyle->SetCanvasColor(0);
gStyle->SetFrameFillColor(0);
pad->cd();
if (writeFitInfo > 0) {
auto* textInfoLeft = new TPaveText(0.12, 0.65, 0.47, 0.89, "NDC");
auto* textInfoRight = new TPaveText(0.6, 0.7, 1., .87, "NDC");
textInfoLeft->SetBorderSize(0);
textInfoLeft->SetFillStyle(0);
textInfoRight->SetBorderSize(0);
textInfoRight->SetFillStyle(0);
textInfoRight->SetTextColor(kBlue);
textInfoLeft->AddText(Form("S = %.0f #pm %.0f ", mRawYield, mRawYieldErr));
textInfoLeft->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr));
if (mTypeOfBkgPdf != NoBkg) {
textInfoLeft->AddText(Form("B (%d#sigma) = %.0f #pm %.0f", mNSigmaForSidebands, mBkgYield, mBkgYieldErr));
textInfoLeft->AddText(Form("S/B (%d#sigma) = %.4g ", mNSigmaForSidebands, mRawYield / mBkgYield));
}
if (mReflPdf != nullptr) {
textInfoLeft->AddText(Form("Refl/Sig = %.3f #pm %.3f ", mReflOverSgn, 0.0));
}
if (mTypeOfBkgPdf != NoBkg) {
textInfoLeft->AddText(Form("Signif (%d#sigma) = %.1f #pm %.1f ", mNSigmaForSidebands, mSignificance, mSignificanceErr));
textInfoLeft->AddText(Form("#chi^{2} / ndf = %.3f", mChiSquareOverNdfTotal));
}
// Fit metrics
auto* textFitMetrics = new TPaveText(0.65, 0.7, 0.9, 0.88, "NDC");
textFitMetrics->SetBorderSize(0);
textFitMetrics->SetFillStyle(0);
textFitMetrics->SetTextSize(0.04);
textFitMetrics->SetTextAlign(33);
textFitMetrics->AddText(Form("S = %.0f #pm %.0f ", mRawYield, mRawYieldErr));
textFitMetrics->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr));
if (mTypeOfBkgPdf != NoBkg) {
textFitMetrics->AddText(Form("B (%d#sigma) = %.0f #pm %.0f", mNSigmaForSidebands, mBkgYield, mBkgYieldErr));
textFitMetrics->AddText(Form("S/B (%d#sigma) = %.4g ", mNSigmaForSidebands, mRawYield / mBkgYield));
textFitMetrics->AddText(Form("Significance (%d#sigma) = %.1f #pm %.1f ", mNSigmaForSidebands, mSignificance, mSignificanceErr));
textFitMetrics->AddText(Form("#chi^{2} / ndf = %.3f", mChiSquareOverNdfTotal));
}
if (mReflPdf != nullptr) {
textFitMetrics->AddText(Form("Refl/Sig = %.3f #pm %.3f ", mReflOverSgn, 0.0));
}
mInvMassFrame->addObject(textFitMetrics);
// Analysis information
auto* textAnalysisInfo = new TPaveText(0.18, 0.78, 0.35, 0.88, "NDC");
textAnalysisInfo->SetBorderSize(0);
textAnalysisInfo->SetFillStyle(0);
textAnalysisInfo->SetTextSize(0.05);
textAnalysisInfo->SetTextAlign(13);
for (const auto& label : plotLabels) {
textAnalysisInfo->AddText(label.c_str());
}
mInvMassFrame->addObject(textAnalysisInfo);
if (writeParInfo) {
// right text box
auto* textSignalPar = new TPaveText(0.18, 0.65, 0.4, 0.75, "NDC");
textSignalPar->SetBorderSize(0);
textSignalPar->SetFillStyle(0);
textSignalPar->SetTextColor(kBlue);
textSignalPar->SetTextAlign(13);
if (mFixedMean) {
textInfoRight->AddText(Form("mean(fixed) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
textSignalPar->AddText(Form("mean(fixed) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
} else {
textInfoRight->AddText(Form("mean(free) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
textSignalPar->AddText(Form("mean(free) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
}
if (mTypeOfSgnPdf == DoubleGaus) {
auto const& baseSigmaSgn = mWorkspace->var("sigma");
if (mFixedSigma) {
textSignalPar->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
} else {
textSignalPar->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
}
if (mFixedSigmaDoubleGaus) {
textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
textInfoRight->AddText(Form("sigma 2(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
textSignalPar->AddText(Form("sigma 2(fixed) = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError()));
} else {
textInfoRight->AddText(Form("sigma(free) = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
textInfoRight->AddText(Form("sigma 2(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
textSignalPar->AddText(Form("sigma 2(free) = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError()));
}
} else if (mFixedSigma) {
textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
textSignalPar->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
} else {
textInfoRight->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
}
mInvMassFrame->addObject(textInfoLeft);
mInvMassFrame->addObject(textInfoRight);
mInvMassFrame->GetYaxis()->SetTitleOffset(1.8);
gPad->SetLeftMargin(0.15);
mInvMassFrame->GetYaxis()->SetTitle(Form("%s", mHistoInvMass->GetYaxis()->GetTitle()));
mInvMassFrame->GetXaxis()->SetTitle(Form("%s", mHistoInvMass->GetXaxis()->GetTitle()));
mInvMassFrame->Draw();
highlightPeakRegion(mInvMassFrame);
if (mHistoTemplateRefl != nullptr) {
mReflFrame->Draw("same");
textSignalPar->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
}
mInvMassFrame->addObject(textSignalPar);
}
mInvMassFrame->GetXaxis()->SetTitleOffset(1.2);
mInvMassFrame->GetYaxis()->SetTitleOffset(1.8);
gPad->SetLeftMargin(0.15);
mInvMassFrame->GetYaxis()->SetTitle(Form("%s", mHistoInvMass->GetYaxis()->GetTitle()));
mInvMassFrame->GetXaxis()->SetTitle(Form("%s", mHistoInvMass->GetXaxis()->GetTitle()));
mInvMassFrame->Draw();
highlightPeakRegion(mInvMassFrame);
if (mHistoTemplateRefl) {
mReflFrame->Draw("same");
}
}

Expand All @@ -626,9 +655,8 @@
textInfo->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr));
textInfo->AddText(Form("mean = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
if (mTypeOfSgnPdf == DoubleGaus) {
auto const& baseSigmaSgn = mWorkspace->var("sigma");
textInfo->AddText(Form("sigma = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError()));
} else {
textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
}
Expand All @@ -637,6 +665,24 @@
highlightPeakRegion(mResidualFrame);
}

// draw ratio on canvas
void HFInvMassFitter::drawRatio(TVirtualPad* pad)
{
pad->cd();
mRatioFrame->GetXaxis()->SetTitleOffset(1.2);
mRatioFrame->GetYaxis()->SetTitleOffset(1.5);
mRatioFrame->GetYaxis()->SetTitle("Fit / Data");
double xMin = mRatioFrame->GetXaxis()->GetXmin();
double xMax = mRatioFrame->GetXaxis()->GetXmax();
auto* line = new TLine(xMin, 1.0, xMax, 1.0);
line->SetLineColor(kGray);
line->SetLineStyle(2);
line->SetLineWidth(2);
mRatioFrame->addObject(line);
mRatioFrame->Draw();
highlightPeakRegion(mRatioFrame);
}

// draw peak region with vertical lines
void HFInvMassFitter::highlightPeakRegion(const RooPlot* plot, Color_t color, Width_t width, Style_t style) const
{
Expand All @@ -646,7 +692,7 @@
double const yMin = plot->GetMinimum();
double const yMax = plot->GetMaximum();
const Double_t mean = mRooMeanSgn->getVal();
const Double_t sigma = mRooSigmaSgn->getVal();
const Double_t sigma = mRooSecSigmaSgn->getVal();
const Double_t minForSgn = mean - mNSigmaForSidebands * sigma;
const Double_t maxForSgn = mean + mNSigmaForSidebands * sigma;
auto* leftLine = new TLine(minForSgn, yMin, minForSgn, yMax);
Expand All @@ -671,7 +717,7 @@
void HFInvMassFitter::countSignal(Double_t& signal, Double_t& signalErr) const
{
const Double_t mean = mRooMeanSgn->getVal();
const Double_t sigma = mRooSigmaSgn->getVal();
const Double_t sigma = mRooSecSigmaSgn->getVal();
const Double_t minForSgn = mean - mNSigmaForSidebands * sigma;
const Double_t maxForSgn = mean + mNSigmaForSidebands * sigma;
const Int_t binForMinSgn = mHistoInvMass->FindBin(minForSgn);
Expand Down Expand Up @@ -763,21 +809,27 @@
case 0: {
sgnPdf = workspace->pdf("sgnFuncGaus");
mRooSigmaSgn = workspace->var("sigma");
mRooSecSigmaSgn = workspace->var("sigma");
mRooMeanSgn = workspace->var("mean");
} break;
case 1: {
sgnPdf = workspace->pdf("sgnFuncDoubleGaus");
mRooSigmaSgn = workspace->var("sigmaDoubleGaus");
mRooSigmaSgn = workspace->var("sigma");
mRooSecSigmaSgn = workspace->var("sigmaDoubleGaus");
mRooMeanSgn = workspace->var("mean");
mRooFracDoubleGaus = workspace->var("fracDoubleGaus");
} break;
case 2: {
sgnPdf = workspace->pdf("sgnFuncGausRatio");
mRooSigmaSgn = workspace->var("sigmaDoubleGausRatio");
mRooSigmaSgn = workspace->var("sigma");
mRooSecSigmaSgn = workspace->var("sigmaDoubleGausRatio");
mRooMeanSgn = workspace->var("mean");
mRooFracDoubleGaus = workspace->var("fracDoubleGausRatio");
} break;
case 3: {
sgnPdf = workspace->pdf("sgnFuncDoublePeak");
mRooSigmaSgn = workspace->var("sigmaSec");
mRooSigmaSgn = workspace->var("sigma");
mRooSecSigmaSgn = workspace->var("sigmaSec");
mRooMeanSgn = workspace->var("meanSec");
} break;
default:
Expand Down Expand Up @@ -818,6 +870,40 @@
pdf->plotOn(mInvMassFrame, Components(namesOfReflPdf.at(mTypeOfReflPdf).c_str()), Name("Refl_c"), LineColor(kGreen));
}

// Calculate fit to data ratio
void HFInvMassFitter::calculateFitToDataRatio() const
{
if (!mInvMassFrame)
return;

// Get the data and fit curves from the frame
auto* dataHist = dynamic_cast<RooHist*>(mInvMassFrame->findObject("data_c"));
auto* fitCurve = dynamic_cast<RooCurve*>(mInvMassFrame->findObject("Tot_c")); // or the relevant fit curve

if (!dataHist || !fitCurve)
return;

RooHist* ratioHist = new RooHist();

for (int i = 0; i < dataHist->GetN(); ++i) {
double x, dataY, dataErr;
dataHist->GetPoint(i, x, dataY);
dataErr = dataHist->GetErrorY(i);

double fitY = fitCurve->Eval(x);

double ratio = dataY != 0 ? fitY / dataY : 0;
double err = dataY != 0 ? ratio * dataErr / dataY : 0;

ratioHist->SetPoint(i, x, ratio);
ratioHist->SetPointError(i, 0, 0, err, err);
}

mRatioFrame->addPlotable(ratioHist, "P");
mRatioFrame->SetMinimum(0.5);
mRatioFrame->SetMaximum(1.5);
}

// Fix reflection pdf
void HFInvMassFitter::setReflFuncFixed()
{
Expand Down
Loading
Loading