From 4e167f2ea181677b4b288fbd5aa273dc40e5c2d8 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Fri, 12 Jun 2026 14:28:34 +0200 Subject: [PATCH 01/34] remove mIntegralHisto from class members --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 7 +++---- PWGHF/D2H/Macros/HFInvMassFitter.h | 1 - 2 files changed, 3 insertions(+), 5 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 304f7b50d7a..b647df811d1 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -147,7 +147,6 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mResidualHist(nullptr), mRatioFrame(nullptr), mWorkspace(nullptr), - mIntegralHisto(0), mIntegralBkg(0), mIntegralSgn(0), mHistoTemplateRefl(nullptr), @@ -192,7 +191,7 @@ HFInvMassFitter::~HFInvMassFitter() void HFInvMassFitter::doFit() { - mIntegralHisto = mHistoInvMass->Integral(mHistoInvMass->FindBin(mMinMass), mHistoInvMass->FindBin(mMaxMass)); + const double integralHisto = mHistoInvMass->Integral(mHistoInvMass->FindBin(mMinMass), mHistoInvMass->FindBin(mMaxMass)); mWorkspace = new RooWorkspace("mWorkspace"); fillWorkspace(*mWorkspace); RooRealVar* mass = mWorkspace->var("mass"); @@ -222,7 +221,7 @@ void HFInvMassFitter::doFit() // fit MC or Data if (mTypeOfBkgPdf == NoBkg) { // MC - const ParameterRanges rooNSgnParamRanges{0., 1.2 * mIntegralHisto, 0.3 * mIntegralHisto}; + const ParameterRanges rooNSgnParamRanges{0., 1.2 * integralHisto, 0.3 * integralHisto}; mRooNSgn = new RooRealVar("mRooNSig", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // signal yield mTotalPdf = new RooAddPdf("mMCFunc", "MC fit function", RooArgList(*sgnPdf), RooArgList(*mRooNSgn)); // create total pdf if (strcmp(mFitOption.c_str(), "Chi2") == 0) { @@ -239,7 +238,7 @@ void HFInvMassFitter::doFit() mRatioFrame = mass->frame(Title(Form("%s", mHistoInvMass->GetTitle()))); calculateFitToDataRatio(); } else { // data - const ParameterRanges rooNBkgParamRanges{0., 1.2 * mIntegralHisto, 0.3 * mIntegralHisto}; + const ParameterRanges rooNBkgParamRanges{0., 1.2 * integralHisto, 0.3 * integralHisto}; mRooNBkg = new RooRealVar("mRooNBkg", "number of background", randomizeInitialParameter(rooNBkgParamRanges), rooNBkgParamRanges.lower, rooNBkgParamRanges.upper); // background yield mBkgPdf = new RooAddPdf("mBkgPdf", "background fit function", RooArgList(*bkgPdf), RooArgList(*mRooNBkg)); if (mTypeOfSgnPdf == GausSec) { // two peak fit diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 6841cf15a55..54bbc83f018 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -261,7 +261,6 @@ class HFInvMassFitter : public TNamed RooHist* mResidualHist; /// residual histogram RooPlot* mRatioFrame; /// fit/data ratio frame RooWorkspace* mWorkspace; /// workspace - double mIntegralHisto; /// integral of histogram to fit double mIntegralBkg; /// integral of background fit function double mIntegralSgn; /// integral of signal fit function TH1* mHistoTemplateRefl; /// reflection histogram From b59dc9e6b512ba6eaad8a7e4d46f58317893d165 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Fri, 19 Jun 2026 08:22:01 +0200 Subject: [PATCH 02/34] checkVectorSize: do not default name parameter --- PWGHF/D2H/Macros/runMassFitter.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index b208e2dbf46..628f2a7d9de 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -218,7 +218,7 @@ void runMassFitter(const std::string& configFileName) std::vector sliceVarLimits(nHistograms + 1); - auto checkVectorSize = [&](const auto& vec, const std::string& name = "", const bool isEmptyOk = false) { + auto checkVectorSize = [&](const auto& vec, const std::string& name, const bool isEmptyOk = false) { if (vec.size() != static_cast(nHistograms)) { if (isEmptyOk && vec.empty()) { return; From ec5f91fb6301400e98a44cdbc0139fc818be9054 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Sat, 20 Jun 2026 21:17:03 +0200 Subject: [PATCH 03/34] DRY sideband ranges --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index b647df811d1..fe192206ea6 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -241,17 +241,19 @@ void HFInvMassFitter::doFit() const ParameterRanges rooNBkgParamRanges{0., 1.2 * integralHisto, 0.3 * integralHisto}; mRooNBkg = new RooRealVar("mRooNBkg", "number of background", randomizeInitialParameter(rooNBkgParamRanges), rooNBkgParamRanges.lower, rooNBkgParamRanges.upper); // background yield mBkgPdf = new RooAddPdf("mBkgPdf", "background fit function", RooArgList(*bkgPdf), RooArgList(*mRooNBkg)); + std::string sbRanges{"SBL,SBR"}; if (mTypeOfSgnPdf == GausSec) { // two peak fit + sbRanges.append(",SEC"); if (strcmp(mFitOption.c_str(), "Chi2") == 0) { - mBkgPdf->chi2FitTo(dataHistogram, Range("SBL,SBR,SEC"), Save()); + mBkgPdf->chi2FitTo(dataHistogram, Range(sbRanges.c_str()), Save()); } else { - mBkgPdf->fitTo(dataHistogram, Range("SBL,SBR,SEC"), Save()); + mBkgPdf->fitTo(dataHistogram, Range(sbRanges.c_str()), Save()); } } else { // single peak fit if (strcmp(mFitOption.c_str(), "Chi2") == 0) { - mBkgPdf->chi2FitTo(dataHistogram, Range("SBL,SBR"), Save()); + mBkgPdf->chi2FitTo(dataHistogram, Range(sbRanges.c_str()), Save()); } else { - mBkgPdf->fitTo(dataHistogram, Range("SBL,SBR"), Save()); + mBkgPdf->fitTo(dataHistogram, Range(sbRanges.c_str()), Save()); } } // define the frame to evaluate background sidebands chi2 (bg pdf needs to be plotted within sideband ranges) From 5e8e9ecf00eef97522fabfd7f558fc50a58bc619 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Sat, 20 Jun 2026 21:19:55 +0200 Subject: [PATCH 04/34] bugfix Ranges() when calculate BG chi2 sidebands --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index fe192206ea6..396642c7604 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -259,7 +259,7 @@ void HFInvMassFitter::doFit() // define the frame to evaluate background sidebands chi2 (bg pdf needs to be plotted within sideband ranges) RooPlot* frameTemporary = mass->frame(Title(Form("%s_temp", mHistoInvMass->GetTitle()))); dataHistogram.plotOn(frameTemporary, Name("data_for_bkgchi2")); - mBkgPdf->plotOn(frameTemporary, Range("SBL", true), Name("Bkg_sidebands")); + mBkgPdf->plotOn(frameTemporary, Range(sbRanges.c_str()), Name("Bkg_sidebands")); mChiSquareOverNdfBkg = frameTemporary->chiSquare("Bkg_sidebands", "data_for_bkgchi2"); // calculate reduced chi2 / NDF of background sidebands (pre-fit) delete frameTemporary; if (mDrawBgPrefit) { From 5aee216e9332599c082a4969db09906913033b65 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Sun, 21 Jun 2026 20:01:40 +0200 Subject: [PATCH 05/34] fix estimate and ranging of BG and SIG integral --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 58 ++++++++++++++++++++-------- PWGHF/D2H/Macros/HFInvMassFitter.h | 1 + 2 files changed, 42 insertions(+), 17 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 396642c7604..043a78d9c81 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -191,7 +191,6 @@ HFInvMassFitter::~HFInvMassFitter() void HFInvMassFitter::doFit() { - const double integralHisto = mHistoInvMass->Integral(mHistoInvMass->FindBin(mMinMass), mHistoInvMass->FindBin(mMaxMass)); mWorkspace = new RooWorkspace("mWorkspace"); fillWorkspace(*mWorkspace); RooRealVar* mass = mWorkspace->var("mass"); @@ -221,7 +220,8 @@ void HFInvMassFitter::doFit() // fit MC or Data if (mTypeOfBkgPdf == NoBkg) { // MC - const ParameterRanges rooNSgnParamRanges{0., 1.2 * integralHisto, 0.3 * integralHisto}; + const double integralHisto = integrateHistoInvMassOverWorkspaceRanges({"full"}); + const ParameterRanges rooNSgnParamRanges{0.5 * integralHisto, 1.5 * integralHisto, integralHisto}; mRooNSgn = new RooRealVar("mRooNSig", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // signal yield mTotalPdf = new RooAddPdf("mMCFunc", "MC fit function", RooArgList(*sgnPdf), RooArgList(*mRooNSgn)); // create total pdf if (strcmp(mFitOption.c_str(), "Chi2") == 0) { @@ -238,24 +238,23 @@ void HFInvMassFitter::doFit() mRatioFrame = mass->frame(Title(Form("%s", mHistoInvMass->GetTitle()))); calculateFitToDataRatio(); } else { // data - const ParameterRanges rooNBkgParamRanges{0., 1.2 * integralHisto, 0.3 * integralHisto}; + const double integralSidebands = integrateHistoInvMassOverWorkspaceRanges({"SBL", "SBR"}); + const ParameterRanges rooNBkgParamRanges{0.5 * integralSidebands, 1.5 * integralSidebands, integralSidebands}; + std::cout << "integralSidebands = " << integralSidebands << "\n"; + std::cout << "rooNBkgParamRanges:\n"; mRooNBkg = new RooRealVar("mRooNBkg", "number of background", randomizeInitialParameter(rooNBkgParamRanges), rooNBkgParamRanges.lower, rooNBkgParamRanges.upper); // background yield mBkgPdf = new RooAddPdf("mBkgPdf", "background fit function", RooArgList(*bkgPdf), RooArgList(*mRooNBkg)); std::string sbRanges{"SBL,SBR"}; if (mTypeOfSgnPdf == GausSec) { // two peak fit sbRanges.append(",SEC"); - if (strcmp(mFitOption.c_str(), "Chi2") == 0) { - mBkgPdf->chi2FitTo(dataHistogram, Range(sbRanges.c_str()), Save()); - } else { - mBkgPdf->fitTo(dataHistogram, Range(sbRanges.c_str()), Save()); - } - } else { // single peak fit - if (strcmp(mFitOption.c_str(), "Chi2") == 0) { - mBkgPdf->chi2FitTo(dataHistogram, Range(sbRanges.c_str()), Save()); - } else { - mBkgPdf->fitTo(dataHistogram, Range(sbRanges.c_str()), Save()); - } } + if (strcmp(mFitOption.c_str(), "Chi2") == 0) { + mBkgPdf->chi2FitTo(dataHistogram, Range(sbRanges.c_str()), Save()); + } else { + mBkgPdf->fitTo(dataHistogram, Range(sbRanges.c_str()), Save()); + } + + std::cout << "mRooNBkg->getVal() = " << mRooNBkg->getVal() << "\n"; // define the frame to evaluate background sidebands chi2 (bg pdf needs to be plotted within sideband ranges) RooPlot* frameTemporary = mass->frame(Title(Form("%s_temp", mHistoInvMass->GetTitle()))); dataHistogram.plotOn(frameTemporary, Name("data_for_bkgchi2")); @@ -275,7 +274,8 @@ void HFInvMassFitter::doFit() checkForSignal(estimatedSignal); // SIG's absolute integral in "bkg" range calculateBackground(mBkgYield, mBkgYieldErr); // BG's absolute integral in "bkg" range - const ParameterRanges rooNSgnParamRanges{0., 1.2 * estimatedSignal, 0.3 * estimatedSignal}; + std::cout << "estimatedSignal = " << estimatedSignal << "\n"; + const ParameterRanges rooNSgnParamRanges{0.5 * estimatedSignal, 1.5 * estimatedSignal, estimatedSignal}; mRooNSgn = new RooRealVar("mNSgn", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // estimated signal yield if (mFixedRawYield > 0) { mRooNSgn->setVal(mFixedRawYield); // fixed signal yield @@ -328,6 +328,8 @@ void HFInvMassFitter::doFit() } else { mTotalPdf->fitTo(dataHistogram); } + std::cout << "mRooNBkg->getVal() = " << mRooNBkg->getVal() << "\n"; + std::cout << "mRooNSgn->getVal() = " << mRooNSgn->getVal() << "\n"; plotBkg(mTotalPdf); 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()); @@ -744,13 +746,15 @@ void HFInvMassFitter::calculateSignal(double& signal, double& errSignal) const // calculate background yield void HFInvMassFitter::calculateBackground(double& bkg, double& errBkg) const { + std::cout << "calculateBackground()\n"; if (mTypeOfBkgPdf == NoBkg) { bkg = 0.; errBkg = 0.; return; } - bkg = mRooNBkg->getVal() * mIntegralBkg; - errBkg = mRooNBkg->getError() * mIntegralBkg; + const double bgCoefficient = mTypeOfSgnPdf == DoubleSidedCrystalBall ? 1. : mIntegralBkg; + bkg = mRooNBkg->getVal() * bgCoefficient; + errBkg = mRooNBkg->getError() * bgCoefficient; } // calculate significance @@ -770,6 +774,7 @@ void HFInvMassFitter::calculateSignificance(double& significance, double& errSig // estimate Signal void HFInvMassFitter::checkForSignal(double& estimatedSignal) { + std::cout << "checkForSignal()\n"; auto const [minForSgn, maxForSgn] = getRangesOfSignal(); int const binForMinSgn = mHistoInvMass->FindBin(minForSgn); int const binForMaxSgn = mHistoInvMass->FindBin(maxForSgn); @@ -780,6 +785,7 @@ void HFInvMassFitter::checkForSignal(double& estimatedSignal) } double bkg{}, errBkg{}; calculateBackground(bkg, errBkg); + std::cout << "sum = " << sum << ", bkg = " << bkg << "\n"; estimatedSignal = sum - bkg; } @@ -1148,5 +1154,23 @@ double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& paramet } } while (result < parameterRanges.lower || result > parameterRanges.upper); + std::cout << "randomizeInitialParameter():\nfrom " << parameterRanges.lower << "\nto " << parameterRanges.upper << "\ninitial " << parameterRanges.initial << "\nsigma " << sigma << "\n"; + std::cout << "randomized to " << result << "\n"; + return result; } + +double HFInvMassFitter::integrateHistoInvMassOverWorkspaceRanges(const std::vector& ranges) const +{ + double sumEntries{0.}; + double sumLengths{0.}; + for (const auto& range : ranges) { + const auto [lo, hi] = mWorkspace->var("mass")->getRange(range.c_str()); + sumEntries += mHistoInvMass->Integral(mHistoInvMass->FindBin(lo), mHistoInvMass->FindBin(hi)); + sumLengths += (hi - lo); + } + const auto [fullLo, fullHi] = mWorkspace->var("mass")->getRange("full"); + const double fullLength = fullHi - fullLo; + + return sumEntries / sumLengths * fullLength; +} diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 54bbc83f018..d23eec33d6c 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -179,6 +179,7 @@ class HFInvMassFitter : public TNamed void highlightPeakRegion(const RooPlot* plot, Color_t color = kGray + 1, Width_t width = 1, Style_t style = 2) const; [[nodiscard]] double randomizeInitialParameter(const ParameterRanges& parameterRanges) const; [[nodiscard]] std::pair getRangesOfSignal() const; + [[nodiscard]] double integrateHistoInvMassOverWorkspaceRanges(const std::vector& ranges) const; TH1* mHistoInvMass; // histogram to fit std::string mFitOption; From c4053d5686c9230d7fff0d8931ccd0edee020ec4 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Tue, 23 Jun 2026 00:21:05 +0200 Subject: [PATCH 06/34] do not crash when one of slice fits is failed --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 18 ++++++++++++++++++ PWGHF/D2H/Macros/runMassFitter.C | 6 +++++- 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 043a78d9c81..2c4496bd637 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -575,6 +575,7 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const // draw fit output void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& plotLabels, bool writeParInfo) { + std::cout << "drawFit() start\n"; gStyle->SetOptStat(0); gStyle->SetCanvasColor(0); gStyle->SetFrameFillColor(0); @@ -634,11 +635,17 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& if (mHistoTemplateRefl) { mReflFrame->Draw("same"); } + std::cout << "drawFit() finish\n"; } // draw residual distribution on canvas void HFInvMassFitter::drawResidual(TVirtualPad* pad) { + std::cout << "drawResidual() start\n"; + if (mResidualFrame == nullptr) { + printf("Warning HFInvMassFitter::drawResidual(): mResidualFrame == nullptr and will not be drawn\n"); + return; + } pad->cd(); mResidualFrame->GetYaxis()->SetTitle(""); auto* textInfo = new TPaveText(0.12, 0.65, 0.47, .89, "NDC"); @@ -655,11 +662,17 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad) mResidualFrame->addObject(textInfo); mResidualFrame->Draw(); highlightPeakRegion(mResidualFrame); + std::cout << "drawResidual() finish\n"; } // draw ratio on canvas void HFInvMassFitter::drawRatio(TVirtualPad* pad) { + std::cout << "drawRatio() start\n"; + if (mRatioFrame == nullptr) { + printf("Warning HFInvMassFitter::drawRatio(): mRatioFrame == nullptr and will not be drawn\n"); + return; + } pad->cd(); mRatioFrame->GetXaxis()->SetTitleOffset(1.2); mRatioFrame->GetYaxis()->SetTitleOffset(1.5); @@ -673,6 +686,7 @@ void HFInvMassFitter::drawRatio(TVirtualPad* pad) mRatioFrame->addObject(line); mRatioFrame->Draw(); highlightPeakRegion(mRatioFrame); + std::cout << "drawRatio() finish\n"; } // draw peak region with vertical lines @@ -1162,13 +1176,17 @@ double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& paramet double HFInvMassFitter::integrateHistoInvMassOverWorkspaceRanges(const std::vector& ranges) const { + std::cout << "integrateHistoInvMassOverWorkspaceRanges()\n"; double sumEntries{0.}; double sumLengths{0.}; for (const auto& range : ranges) { const auto [lo, hi] = mWorkspace->var("mass")->getRange(range.c_str()); sumEntries += mHistoInvMass->Integral(mHistoInvMass->FindBin(lo), mHistoInvMass->FindBin(hi)); sumLengths += (hi - lo); + std::cout << "lo = " << lo << ", hi = " << hi << "\n"; + std::cout << "sumEntries = " << sumEntries << "\n"; } + mHistoInvMass->SaveAs("mHistoInvMass.root"); const auto [fullLo, fullHi] = mWorkspace->var("mass")->getRange("full"); const double fullLength = fullHi - fullLo; diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 628f2a7d9de..746bb9114f3 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -542,7 +542,11 @@ void runMassFitter(const std::string& configFileName) setDscbParameter(dscbNRLower, &HFInvMassFitter::setDscbNRLowLimit); setDscbParameter(dscbNRUpper, &HFInvMassFitter::setDscbNRUpLimit); - massFitter->doFit(); + try { + massFitter->doFit(); + } catch (...) { + std::cout << "exception caught while doing fit\n"; + } auto drawOnCanvas = [&](std::vector& canvas, std::function drawer) { if (nHistograms > 1) { From cfc54110612b712078169b81fdce49db42d69036 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 24 Jun 2026 02:49:17 +0200 Subject: [PATCH 07/34] prefit BG sidebands with chi2 always --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 48 +++++++++++++++++++++------- PWGHF/D2H/Macros/HFInvMassFitter.h | 1 + 2 files changed, 38 insertions(+), 11 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 2c4496bd637..4ab13efc817 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -50,6 +50,7 @@ #include #include +#include #include #include #include @@ -199,13 +200,12 @@ void HFInvMassFitter::doFit() if (mTypeOfBkgPdf == NoBkg) { // MC mass->setRange("signal", mMass - 3. * mSigmaSgn, mMass + 3. * mSigmaSgn); } else { + mass->setRange("SBL", mMinMass, mMass - mNSigmaForSidebands * mSigmaSgn); if (mTypeOfSgnPdf == GausSec) { // Second Peak fit range - mass->setRange("SBL", mMinMass, mMass - mNSigmaForSidebands * mSigmaSgn); mass->setRange("SBR", mMass + mNSigmaForSidebands * mSigmaSgn, mSecMass - mNSigmaForSidebands * mSecSigma); mass->setRange("SEC", mSecMass + mNSigmaForSidebands * mSecSigma, mMaxMass); mass->setRange("signal", mSecMass - mNSigmaForSgn * mSecSigma, mSecMass + mNSigmaForSgn * mSecSigma); } else { // Single Peak fit range - mass->setRange("SBL", mMinMass, mMass - mNSigmaForSidebands * mSigmaSgn); mass->setRange("SBR", mMass + mNSigmaForSidebands * mSigmaSgn, mMaxMass); mass->setRange("signal", mMass - mNSigmaForSgn * mSigmaSgn, mMass + mNSigmaForSgn * mSigmaSgn); } @@ -215,12 +215,17 @@ void HFInvMassFitter::doFit() mInvMassFrame = mass->frame(Title(Form("%s", mHistoInvMass->GetTitle()))); // define the frame to plot dataHistogram.plotOn(mInvMassFrame, Name("data_c")); // plot data histogram on the frame + TH1* histoInvMassSB = dynamic_cast(mHistoInvMass->Clone()); + cutRangesFromHisto(histoInvMassSB, {"SBL", "SBR"}); + histoInvMassSB->SaveAs("histoInvMassSB.root"); + RooDataHist sbHistogram("sbHistogram", "sb", *mass, Import(*histoInvMassSB)); + RooAbsPdf* bkgPdf = createBackgroundFitFunction(mWorkspace); // Create background pdf RooAbsPdf* sgnPdf = createSignalFitFunction(mWorkspace); // Create signal pdf + const double integralHisto = integrateHistoInvMassOverWorkspaceRanges({"full"}); // fit MC or Data if (mTypeOfBkgPdf == NoBkg) { // MC - const double integralHisto = integrateHistoInvMassOverWorkspaceRanges({"full"}); const ParameterRanges rooNSgnParamRanges{0.5 * integralHisto, 1.5 * integralHisto, integralHisto}; mRooNSgn = new RooRealVar("mRooNSig", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // signal yield mTotalPdf = new RooAddPdf("mMCFunc", "MC fit function", RooArgList(*sgnPdf), RooArgList(*mRooNSgn)); // create total pdf @@ -248,13 +253,9 @@ void HFInvMassFitter::doFit() if (mTypeOfSgnPdf == GausSec) { // two peak fit sbRanges.append(",SEC"); } - if (strcmp(mFitOption.c_str(), "Chi2") == 0) { - mBkgPdf->chi2FitTo(dataHistogram, Range(sbRanges.c_str()), Save()); - } else { - mBkgPdf->fitTo(dataHistogram, Range(sbRanges.c_str()), Save()); - } - + mBkgPdf->chi2FitTo(sbHistogram, DataError(RooAbsData::SumW2), Save()); std::cout << "mRooNBkg->getVal() = " << mRooNBkg->getVal() << "\n"; + // define the frame to evaluate background sidebands chi2 (bg pdf needs to be plotted within sideband ranges) RooPlot* frameTemporary = mass->frame(Title(Form("%s_temp", mHistoInvMass->GetTitle()))); dataHistogram.plotOn(frameTemporary, Name("data_for_bkgchi2")); @@ -263,7 +264,7 @@ void HFInvMassFitter::doFit() delete frameTemporary; if (mDrawBgPrefit) { RooAbsPdf* bkgPdfPrefit = dynamic_cast(mBkgPdf->Clone()); - bkgPdfPrefit->plotOn(mInvMassFrame, Range("full"), Name("Bkg_c_prefit"), LineColor(kGray)); + bkgPdfPrefit->plotOn(mInvMassFrame, Range("full"), Normalization(mRooNBkg->getVal(), RooAbsReal::NumEvent), Name("Bkg_c_prefit"), LineColor(kGray)); delete bkgPdfPrefit; } @@ -275,7 +276,7 @@ void HFInvMassFitter::doFit() calculateBackground(mBkgYield, mBkgYieldErr); // BG's absolute integral in "bkg" range std::cout << "estimatedSignal = " << estimatedSignal << "\n"; - const ParameterRanges rooNSgnParamRanges{0.5 * estimatedSignal, 1.5 * estimatedSignal, estimatedSignal}; + const ParameterRanges rooNSgnParamRanges{0.1 * estimatedSignal, 20 * estimatedSignal, estimatedSignal}; mRooNSgn = new RooRealVar("mNSgn", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // estimated signal yield if (mFixedRawYield > 0) { mRooNSgn->setVal(mFixedRawYield); // fixed signal yield @@ -1192,3 +1193,28 @@ double HFInvMassFitter::integrateHistoInvMassOverWorkspaceRanges(const std::vect return sumEntries / sumLengths * fullLength; } + +void HFInvMassFitter::cutRangesFromHisto(TH1* histo, const std::vector& ranges) const +{ + auto* mass = mWorkspace->var("mass"); + + for (int iBin = 1, nBins = histo->GetNbinsX(); iBin <= nBins; ++iBin) { + const double binLow = histo->GetBinLowEdge(iBin); + const double binHigh = binLow + histo->GetBinWidth(iBin); + + bool overlapsAnyRange = false; + for (const auto& range : ranges) { + const double rangeMin = mass->getMin(range.c_str()); + const double rangeMax = mass->getMax(range.c_str()); + if (std::max(binLow, rangeMin) <= std::min(binHigh, rangeMax)) { + overlapsAnyRange = true; + break; + } + } + + if (!overlapsAnyRange) { + histo->SetBinContent(iBin, 0.); + histo->SetBinError(iBin, 1.e9); + } + } +} diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index d23eec33d6c..e7bdff26e97 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -180,6 +180,7 @@ class HFInvMassFitter : public TNamed [[nodiscard]] double randomizeInitialParameter(const ParameterRanges& parameterRanges) const; [[nodiscard]] std::pair getRangesOfSignal() const; [[nodiscard]] double integrateHistoInvMassOverWorkspaceRanges(const std::vector& ranges) const; + void cutRangesFromHisto(TH1* histo, const std::vector& ranges) const; TH1* mHistoInvMass; // histogram to fit std::string mFitOption; From 4720e6ba7d1f54317e6608da9b458ae85b9458bb Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 24 Jun 2026 05:25:39 +0200 Subject: [PATCH 08/34] revert saving debug histograms --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 2 -- 1 file changed, 2 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 4ab13efc817..3e1b0f8b9d4 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -217,7 +217,6 @@ void HFInvMassFitter::doFit() TH1* histoInvMassSB = dynamic_cast(mHistoInvMass->Clone()); cutRangesFromHisto(histoInvMassSB, {"SBL", "SBR"}); - histoInvMassSB->SaveAs("histoInvMassSB.root"); RooDataHist sbHistogram("sbHistogram", "sb", *mass, Import(*histoInvMassSB)); RooAbsPdf* bkgPdf = createBackgroundFitFunction(mWorkspace); // Create background pdf @@ -1187,7 +1186,6 @@ double HFInvMassFitter::integrateHistoInvMassOverWorkspaceRanges(const std::vect std::cout << "lo = " << lo << ", hi = " << hi << "\n"; std::cout << "sumEntries = " << sumEntries << "\n"; } - mHistoInvMass->SaveAs("mHistoInvMass.root"); const auto [fullLo, fullHi] = mWorkspace->var("mass")->getRange("full"); const double fullLength = fullHi - fullLo; From f487afd878ba88b7a122af06f006f6eab4bb147c Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Thu, 2 Jul 2026 14:02:02 +0200 Subject: [PATCH 09/34] DRY canvas init, add ticks on top and right axes --- PWGHF/D2H/Macros/runMassFitter.C | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 746bb9114f3..ca2741c07e0 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -442,24 +442,24 @@ void runMassFitter(const std::string& configFileName) int constexpr NCanvasesMax = 20; // do not put more than 20 bins per canvas to make them visible const int nCanvases = std::ceil(static_cast(nHistograms) / NCanvasesMax); + const int nPads = (nCanvases == 1) ? nHistograms : NCanvasesMax; std::vector canvasMass(nCanvases); std::vector canvasResiduals(nCanvases); std::vector canvasRatio(nCanvases); std::vector canvasRefl(nCanvases); - for (int iCanvas = 0; iCanvas < nCanvases; iCanvas++) { - const int nPads = (nCanvases == 1) ? nHistograms : NCanvasesMax; - canvasMass[iCanvas] = new TCanvas(Form("canvasMass%d", iCanvas), Form("canvasMass%d", iCanvas), canvasSize[0], canvasSize[1]); - divideCanvas(canvasMass[iCanvas], nPads); - - canvasResiduals[iCanvas] = new TCanvas(Form("canvasResiduals%d", iCanvas), Form("canvasResiduals%d", iCanvas), canvasSize[0], canvasSize[1]); - divideCanvas(canvasResiduals[iCanvas], nPads); - - canvasRatio[iCanvas] = new TCanvas(Form("canvasRatio%d", iCanvas), Form("canvasRatio%d", iCanvas), canvasSize[0], canvasSize[1]); - divideCanvas(canvasRatio[iCanvas], nPads); - if (enableRefl) { - canvasRefl[iCanvas] = new TCanvas(Form("canvasRefl%d", iCanvas), Form("canvasRefl%d", iCanvas), canvasSize[0], canvasSize[1]); - divideCanvas(canvasRefl[iCanvas], nPads); + std::vector*> canvasTypes{&canvasMass, &canvasResiduals, &canvasRatio}; + std::vector canvasTypeNames{"canvasMass", "canvasResiduals", "canvasRatio"}; + if (enableRefl) { + canvasTypes.push_back(&canvasRefl); + canvasTypeNames.push_back("canvasRefl"); + } + for (int iCanvasType = 0, nCanvasTypeNames = canvasTypes.size(); iCanvasType < nCanvasTypeNames; ++iCanvasType) { + for (int iCanvas = 0; iCanvas < nCanvases; iCanvas++) { + auto& canvas = (*canvasTypes.at(iCanvasType))[iCanvas]; + canvas = new TCanvas(Form((canvasTypeNames.at(iCanvasType) + "%d").c_str(), iCanvas), Form((canvasTypeNames.at(iCanvasType) + "%d").c_str(), iCanvas), canvasSize[0], canvasSize[1]); + canvas->SetTicks(1, 1); + divideCanvas(canvas, nPads); } } From 6959eba044aa66a2d9fe36a88a718e41b177742e Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Thu, 2 Jul 2026 14:18:35 +0200 Subject: [PATCH 10/34] clean previous commit f487afd --- PWGHF/D2H/Macros/runMassFitter.C | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index ca2741c07e0..d5c7c5de77a 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -441,7 +441,7 @@ void runMassFitter(const std::string& configFileName) } int constexpr NCanvasesMax = 20; // do not put more than 20 bins per canvas to make them visible - const int nCanvases = std::ceil(static_cast(nHistograms) / NCanvasesMax); + const int nCanvases = (nHistograms + NCanvasesMax - 1) / NCanvasesMax; const int nPads = (nCanvases == 1) ? nHistograms : NCanvasesMax; std::vector canvasMass(nCanvases); std::vector canvasResiduals(nCanvases); @@ -449,15 +449,17 @@ void runMassFitter(const std::string& configFileName) std::vector canvasRefl(nCanvases); std::vector*> canvasTypes{&canvasMass, &canvasResiduals, &canvasRatio}; - std::vector canvasTypeNames{"canvasMass", "canvasResiduals", "canvasRatio"}; + std::vector canvasTypeNames{"canvasMass", "canvasResiduals", "canvasRatio"}; if (enableRefl) { canvasTypes.push_back(&canvasRefl); canvasTypeNames.push_back("canvasRefl"); } - for (int iCanvasType = 0, nCanvasTypeNames = canvasTypes.size(); iCanvasType < nCanvasTypeNames; ++iCanvasType) { + for (int iCanvasType = 0, nCanvasTypeNames = static_cast(canvasTypes.size()); iCanvasType < nCanvasTypeNames; ++iCanvasType) { + const auto canvasTypeName = canvasTypeNames[iCanvasType]; for (int iCanvas = 0; iCanvas < nCanvases; iCanvas++) { - auto& canvas = (*canvasTypes.at(iCanvasType))[iCanvas]; - canvas = new TCanvas(Form((canvasTypeNames.at(iCanvasType) + "%d").c_str(), iCanvas), Form((canvasTypeNames.at(iCanvasType) + "%d").c_str(), iCanvas), canvasSize[0], canvasSize[1]); + auto& canvas = (*canvasTypes[iCanvasType])[iCanvas]; + const auto canvasName = Form("%s%d", canvasTypeName, iCanvas); + canvas = new TCanvas(canvasName, canvasName, canvasSize[0], canvasSize[1]); canvas->SetTicks(1, 1); divideCanvas(canvas, nPads); } From 63efdcb12f1164674055d2b5889d75d643bb1ac6 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Thu, 2 Jul 2026 14:38:45 +0200 Subject: [PATCH 11/34] simplify checkForSignal() --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 18 +++++------------- 1 file changed, 5 insertions(+), 13 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 3e1b0f8b9d4..82ca08cd2d8 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -275,7 +275,7 @@ void HFInvMassFitter::doFit() calculateBackground(mBkgYield, mBkgYieldErr); // BG's absolute integral in "bkg" range std::cout << "estimatedSignal = " << estimatedSignal << "\n"; - const ParameterRanges rooNSgnParamRanges{0.1 * estimatedSignal, 20 * estimatedSignal, estimatedSignal}; + const ParameterRanges rooNSgnParamRanges{0.001 * estimatedSignal, 1000 * estimatedSignal, estimatedSignal}; mRooNSgn = new RooRealVar("mNSgn", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // estimated signal yield if (mFixedRawYield > 0) { mRooNSgn->setVal(mFixedRawYield); // fixed signal yield @@ -789,18 +789,10 @@ void HFInvMassFitter::calculateSignificance(double& significance, double& errSig void HFInvMassFitter::checkForSignal(double& estimatedSignal) { std::cout << "checkForSignal()\n"; - auto const [minForSgn, maxForSgn] = getRangesOfSignal(); - int const binForMinSgn = mHistoInvMass->FindBin(minForSgn); - int const binForMaxSgn = mHistoInvMass->FindBin(maxForSgn); - - double sum = 0; - for (int i = binForMinSgn; i <= binForMaxSgn; i++) { - sum += mHistoInvMass->GetBinContent(i); - } - double bkg{}, errBkg{}; - calculateBackground(bkg, errBkg); - std::cout << "sum = " << sum << ", bkg = " << bkg << "\n"; - estimatedSignal = sum - bkg; + const double integralHisto = integrateHistoInvMassOverWorkspaceRanges({"full"}); + const double bkg = mRooNBkg->getVal(); + std::cout << "integralHisto = " << integralHisto << ", bkg = " << bkg << "\n"; + estimatedSignal = integralHisto - bkg; } // Estimate ranges where signal is located to be used in countSignal() and checkForSignal() From 144bb0e3049a714b6efbb606337faeafde651d9b Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Thu, 2 Jul 2026 15:26:06 +0200 Subject: [PATCH 12/34] integrateHistoInvMassOverWorkspaceRanges(): more accurate account for edge bins --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 82ca08cd2d8..d51e4a4423e 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -275,7 +275,7 @@ void HFInvMassFitter::doFit() calculateBackground(mBkgYield, mBkgYieldErr); // BG's absolute integral in "bkg" range std::cout << "estimatedSignal = " << estimatedSignal << "\n"; - const ParameterRanges rooNSgnParamRanges{0.001 * estimatedSignal, 1000 * estimatedSignal, estimatedSignal}; + const ParameterRanges rooNSgnParamRanges{0.1 * estimatedSignal, 10 * estimatedSignal, estimatedSignal}; mRooNSgn = new RooRealVar("mNSgn", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // estimated signal yield if (mFixedRawYield > 0) { mRooNSgn->setVal(mFixedRawYield); // fixed signal yield @@ -1173,7 +1173,11 @@ double HFInvMassFitter::integrateHistoInvMassOverWorkspaceRanges(const std::vect double sumLengths{0.}; for (const auto& range : ranges) { const auto [lo, hi] = mWorkspace->var("mass")->getRange(range.c_str()); - sumEntries += mHistoInvMass->Integral(mHistoInvMass->FindBin(lo), mHistoInvMass->FindBin(hi)); + const auto binLo = mHistoInvMass->FindBin(lo); + const auto binHi = mHistoInvMass->FindBin(hi); + sumEntries += mHistoInvMass->Integral(binLo + 1, binHi - 1); + sumEntries += mHistoInvMass->GetBinContent(binLo) * (mHistoInvMass->GetBinLowEdge(binLo + 1) - lo) / mHistoInvMass->GetBinWidth(binLo); + sumEntries += mHistoInvMass->GetBinContent(binHi) * (hi - mHistoInvMass->GetBinLowEdge(binHi)) / mHistoInvMass->GetBinWidth(binHi); sumLengths += (hi - lo); std::cout << "lo = " << lo << ", hi = " << hi << "\n"; std::cout << "sumEntries = " << sumEntries << "\n"; From 977fda38d35bda6908633f96edb32291d716ae03 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Thu, 2 Jul 2026 15:27:52 +0200 Subject: [PATCH 13/34] check massMin, massMax for being within hMass ranges --- PWGHF/D2H/Macros/runMassFitter.C | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index d5c7c5de77a..b1e94c56229 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -481,6 +481,17 @@ void runMassFitter(const std::string& configFileName) hMassSgn[iSliceVar]->Rebin(nRebin[iSliceVar]); } + const auto hMassLo = hMass[iSliceVar]->GetXaxis()->GetXmin(); + if (massMin[iSliceVar] < hMassLo) { + printf("Warning! massMin[%d] is less than hMass[%d] left edge (%f vs %f) and will be assigned the value of the latter\n", iSliceVar, iSliceVar, massMin[iSliceVar], hMassLo); + massMin[iSliceVar] = hMassLo; + } + const auto hMassUp = hMass[iSliceVar]->GetXaxis()->GetXmax(); + if (massMax[iSliceVar] > hMassUp) { + printf("Warning! massMax[%d] is greater than hMass[%d] right edge (%f vs %f) and will be assigned the value of the latter\n", iSliceVar, iSliceVar, massMax[iSliceVar], hMassUp); + massMax[iSliceVar] = hMassUp; + } + double reflOverSgn = 0; HFInvMassFitter* massFitter = new HFInvMassFitter(hMass[iSliceVar], massMin[iSliceVar], massMax[iSliceVar], bkgFunc[iSliceVar], sgnFunc[iSliceVar], randomSeed); From 991a667270f6aa85d2766ddcacbe2f952a2d7c69 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 8 Jul 2026 12:15:27 +0200 Subject: [PATCH 14/34] DRY fit w/ and w/o reflection --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 56 ++++++++++++---------------- 1 file changed, 23 insertions(+), 33 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index d51e4a4423e..2d20b4353e9 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -283,8 +283,9 @@ void HFInvMassFitter::doFit() } mSgnPdf = new RooAddPdf("mSgnPdf", "signal fit function", RooArgList(*sgnPdf), RooArgList(*mRooNSgn)); // create reflection template and fit to reflection + RooAbsPdf* reflPdf{nullptr}; if (mHistoTemplateRefl != nullptr) { - RooAbsPdf* reflPdf = createReflectionFitFunction(mWorkspace); // create reflection pdf + reflPdf = createReflectionFitFunction(mWorkspace); // create reflection pdf RooDataHist reflHistogram("reflHistogram", "refl for fit", *mass, Import(*mHistoTemplateRefl)); mReflFrame = mass->frame(); mReflOnlyFrame = mass->frame(Title(Form("%s", mHistoTemplateRefl->GetTitle()))); @@ -303,44 +304,33 @@ void HFInvMassFitter::doFit() mRooNRefl->setConstant(true); setReflFuncFixed(); // fix reflection pdf parameter mTotalPdf = new RooAddPdf("mTotalPdf", "background + signal + reflection fit function", RooArgList(*bkgPdf, *sgnPdf, *reflPdf), RooArgList(*mRooNBkg, *mRooNSgn, *mRooNRefl)); - if (strcmp(mFitOption.c_str(), "Chi2") == 0) { - mTotalPdf->chi2FitTo(dataHistogram); - } else { - mTotalPdf->fitTo(dataHistogram); - } - mTotalPdf->plotOn(mInvMassFrame, Name("Tot_c")); + } else { + mTotalPdf = new RooAddPdf("mTotalPdf", "background + signal pdf", RooArgList(*bkgPdf, *sgnPdf), RooArgList(*mRooNBkg, *mRooNSgn)); + } + if (strcmp(mFitOption.c_str(), "Chi2") == 0) { + mTotalPdf->chi2FitTo(dataHistogram); + } else { + mTotalPdf->fitTo(dataHistogram); + } + std::cout << "mRooNBkg->getVal() = " << mRooNBkg->getVal() << "\n"; + std::cout << "mRooNSgn->getVal() = " << mRooNSgn->getVal() << "\n"; + plotBkg(mTotalPdf); + mTotalPdf->plotOn(mInvMassFrame, Name("Tot_c"), LineColor(kBlue)); + if (mHistoTemplateRefl != nullptr) { mReflPdf = new RooAddPdf("mReflPdf", "reflection fit function", RooArgList(*reflPdf), RooArgList(*mRooNRefl)); RooAddPdf const reflBkgPdf("reflBkgPdf", "reflBkgPdf", RooArgList(*bkgPdf, *reflPdf), RooArgList(*mRooNBkg, *mRooNRefl)); reflBkgPdf.plotOn(mInvMassFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineStyle(7), LineColor(kRed + 1), Name("ReflBkg_c")); - plotBkg(mTotalPdf); // plot bkg pdf in total pdf - plotRefl(mTotalPdf); // plot reflection in total pdf - mChiSquareOverNdfTotal = mInvMassFrame->chiSquare("Tot_c", "data_c"); // calculate reduced chi2 / NDF - - // plot residual distribution - mResidualHist = mInvMassFrame->residHist("data_c", "ReflBkg_c"); - mResidualFrame = mass->frame(Title("Residual Distribution")); - mResidualFrame->addPlotable(mResidualHist, "p"); - mSgnPdf->plotOn(mResidualFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineColor(kBlue)); + plotRefl(mTotalPdf); // plot reflection in total pdf } else { - mTotalPdf = new RooAddPdf("mTotalPdf", "background + signal pdf", RooArgList(*bkgPdf, *sgnPdf), RooArgList(*mRooNBkg, *mRooNSgn)); - if (strcmp(mFitOption.c_str(), "Chi2") == 0) { - mTotalPdf->chi2FitTo(dataHistogram); - } else { - mTotalPdf->fitTo(dataHistogram); - } - std::cout << "mRooNBkg->getVal() = " << mRooNBkg->getVal() << "\n"; - std::cout << "mRooNSgn->getVal() = " << mRooNSgn->getVal() << "\n"; - plotBkg(mTotalPdf); - 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")); - mResidualHist = mInvMassFrame->residHist("data_c", "Bkg_c"); - mResidualFrame->addPlotable(mResidualHist, "P"); - mSgnPdf->plotOn(mResidualFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineColor(kBlue)); } + mChiSquareOverNdfTotal = mInvMassFrame->chiSquare("Tot_c", "data_c"); // calculate reduced chi2 / NDF + // plot residual distribution + mResidualFrame = mass->frame(Title("Residual Distribution")); + mResidualHist = mInvMassFrame->residHist("data_c", mHistoTemplateRefl ? "ReflBkg_c" : "Bkg_c"); + mResidualFrame->addPlotable(mResidualHist, "P"); + mSgnPdf->plotOn(mResidualFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineColor(kBlue)); + mass->setRange("bkgForSignificance", mRooMeanSgn->getVal() - mNSigmaForSgn * mRooSecSigmaSgn->getVal(), mRooMeanSgn->getVal() + mNSigmaForSgn * mRooSecSigmaSgn->getVal()); bkgIntegral = mBkgPdf->createIntegral(*mass, NormSet(*mass), Range("bkgForSignificance")); mIntegralBkg = bkgIntegral->getValV(); From 15521a89c41f753761e1bc79fc622e4b54e0074a Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 8 Jul 2026 12:21:28 +0200 Subject: [PATCH 15/34] make residual's frame title more informative --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 2d20b4353e9..661f127b88d 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -326,7 +326,7 @@ void HFInvMassFitter::doFit() } mChiSquareOverNdfTotal = mInvMassFrame->chiSquare("Tot_c", "data_c"); // calculate reduced chi2 / NDF // plot residual distribution - mResidualFrame = mass->frame(Title("Residual Distribution")); + mResidualFrame = mass->frame(Title(Form("%s", mHistoInvMass->GetTitle()))); mResidualHist = mInvMassFrame->residHist("data_c", mHistoTemplateRefl ? "ReflBkg_c" : "Bkg_c"); mResidualFrame->addPlotable(mResidualHist, "P"); mSgnPdf->plotOn(mResidualFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineColor(kBlue)); From 13863b9ab4d3a2e05f329508a9e25c4a035e6906 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Thu, 9 Jul 2026 16:27:51 +0200 Subject: [PATCH 16/34] residual plot: add Y-axis title --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 661f127b88d..4c06c8f94a8 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -618,6 +618,7 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& mInvMassFrame->GetXaxis()->SetTitleOffset(1.2); mInvMassFrame->GetYaxis()->SetTitleOffset(1.8); gPad->SetLeftMargin(0.15); + gPad->SetRightMargin(0.02); mInvMassFrame->GetYaxis()->SetTitle(Form("%s", mHistoInvMass->GetYaxis()->GetTitle())); mInvMassFrame->GetXaxis()->SetTitle(Form("%s", mHistoInvMass->GetXaxis()->GetTitle())); mInvMassFrame->Draw(); @@ -637,7 +638,7 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad) return; } pad->cd(); - mResidualFrame->GetYaxis()->SetTitle(""); + mResidualFrame->GetYaxis()->SetTitle(mHistoInvMass->GetYaxis()->GetTitle()); auto* textInfo = new TPaveText(0.12, 0.65, 0.47, .89, "NDC"); textInfo->SetBorderSize(0); textInfo->SetFillStyle(0); @@ -650,6 +651,9 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad) textInfo->AddText(Form("#sigma_{2} = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); } mResidualFrame->addObject(textInfo); + gPad->SetLeftMargin(0.15); + gPad->SetRightMargin(0.02); + mResidualFrame->GetYaxis()->SetTitleOffset(1.8); mResidualFrame->Draw(); highlightPeakRegion(mResidualFrame); std::cout << "drawResidual() finish\n"; @@ -674,6 +678,8 @@ void HFInvMassFitter::drawRatio(TVirtualPad* pad) line->SetLineStyle(2); line->SetLineWidth(2); mRatioFrame->addObject(line); + gPad->SetLeftMargin(0.15); + gPad->SetRightMargin(0.02); mRatioFrame->Draw(); highlightPeakRegion(mRatioFrame); std::cout << "drawRatio() finish\n"; From 2f3b816d210b36a838a91c80c76a9f2455229e54 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Tue, 21 Jul 2026 15:25:39 +0200 Subject: [PATCH 17/34] add hFitResult output histogram --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 20 +++++++++++++++++--- PWGHF/D2H/Macros/HFInvMassFitter.h | 4 ++++ PWGHF/D2H/Macros/runMassFitter.C | 24 ++++++++++++++++++++---- 3 files changed, 41 insertions(+), 7 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 4c06c8f94a8..8f2e17d8b6f 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -154,7 +154,9 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mDrawBgPrefit(false), mHighlightPeakRegion(false), mRandomSeed(randomSeed), - mRandomGen(nullptr) + mRandomGen(nullptr), + mFitStatus(-999), + mCovQual(-999) { // standard constructor mHistoInvMass = histoToFit; @@ -252,7 +254,9 @@ void HFInvMassFitter::doFit() if (mTypeOfSgnPdf == GausSec) { // two peak fit sbRanges.append(",SEC"); } + std::cout << "Start prefit of BG sidebands\n"; mBkgPdf->chi2FitTo(sbHistogram, DataError(RooAbsData::SumW2), Save()); + std::cout << "Finish prefit of BG sidebands\n"; std::cout << "mRooNBkg->getVal() = " << mRooNBkg->getVal() << "\n"; // define the frame to evaluate background sidebands chi2 (bg pdf needs to be plotted within sideband ranges) @@ -307,11 +311,21 @@ void HFInvMassFitter::doFit() } else { mTotalPdf = new RooAddPdf("mTotalPdf", "background + signal pdf", RooArgList(*bkgPdf, *sgnPdf), RooArgList(*mRooNBkg, *mRooNSgn)); } + std::cout << "Start total fit\n"; + RooFitResult* fitResult{nullptr}; if (strcmp(mFitOption.c_str(), "Chi2") == 0) { - mTotalPdf->chi2FitTo(dataHistogram); + fitResult = mTotalPdf->chi2FitTo(dataHistogram, Save()); } else { - mTotalPdf->fitTo(dataHistogram); + fitResult = mTotalPdf->fitTo(dataHistogram, Save()); } + std::cout << "Finish total fit\n"; + std::cout << "Status = " << fitResult->status() << "\n"; + std::cout << "CovQual = " << fitResult->covQual() << "\n"; + std::cout << "EDM = " << fitResult->edm() << "\n"; + std::cout << "minNLL = " << fitResult->minNll() << "\n"; + fitResult->Print("v"); + mFitStatus = fitResult->status(); + mCovQual = fitResult->covQual(); std::cout << "mRooNBkg->getVal() = " << mRooNBkg->getVal() << "\n"; std::cout << "mRooNSgn->getVal() = " << mRooNSgn->getVal() << "\n"; plotBkg(mTotalPdf); diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index e7bdff26e97..434bd5d87fc 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -161,6 +161,8 @@ class HFInvMassFitter : public TNamed [[nodiscard]] double getFracDoubleGaus() const { return mRooFracDoubleGaus->getVal(); } [[nodiscard]] double getFracDoubleGausUncertainty() const { return mRooFracDoubleGaus->getError(); } [[nodiscard]] double getReflOverSig() const { return mReflPdf != nullptr ? mReflOverSgn : 0.; } + [[nodiscard]] int getFitStatus() const { return mFitStatus; } + [[nodiscard]] int getCovQual() const { return mCovQual; } void calculateSignal(double& signal, double& signalErr) const; void countSignal(double& signal, double& signalErr) const; void calculateBackground(double& bkg, double& bkgErr) const; @@ -270,6 +272,8 @@ class HFInvMassFitter : public TNamed bool mHighlightPeakRegion; /// draw vertical lines showing the peak region (usually +- 3 sigma) int mRandomSeed; /// seed for random engine for fit's initial parameters randomization TRandom3* mRandomGen; /// engine for fit's initial parameters randomization + int mFitStatus; /// fit result status, see https://root-forum.cern.ch/t/meaning-of-values-returned-by-roofitresult-status/16355/2 + int mCovQual; /// fit result covariance matrix quality, see https://root.cern.ch/doc/v620/Minuit2Minimizer_8cxx_source.html#l01121 ClassDefOverride(HFInvMassFitter, 1); }; diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index b1e94c56229..339aebfb8cd 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -385,15 +385,27 @@ void runMassFitter(const std::string& configFileName) ConfigRandomSeed, NConfigsToSave }; + enum { + FitResultStatus = 1, + FitResultCovQual, + NFitResultsToSave + }; auto* hFitConfig = new TH2F("hfitConfig", "Fit Configurations", NConfigsToSave - 1, 0, NConfigsToSave - 1, nHistograms, sliceVarLimits.data()); const char* hFitConfigXLabel[NConfigsToSave - 1] = {"mass min", "mass max", "rebin num", "fix sigma", "bkg func", "sgn func", "rnd seed"}; - hFitConfig->SetStats(false); + auto* hFitResult = new TH2F("hfitResult", "Fit Result", NFitResultsToSave - 1, 0, NFitResultsToSave - 1, nHistograms, sliceVarLimits.data()); + const char* hFitResultXLabel[NConfigsToSave - 1] = {"status", "cov qual"}; for (int i = 0; i < NConfigsToSave - 1; i++) { hFitConfig->GetXaxis()->SetBinLabel(i + 1, hFitConfigXLabel[i]); } - hFitConfig->LabelsDeflate("X"); - hFitConfig->LabelsDeflate("Y"); - hFitConfig->LabelsOption("v"); + for (int i = 0; i < NFitResultsToSave - 1; i++) { + hFitResult->GetXaxis()->SetBinLabel(i + 1, hFitResultXLabel[i]); + } + for (const auto& h : {hFitConfig, hFitResult}) { + h->SetStats(false); + h->LabelsDeflate("X"); + h->LabelsDeflate("Y"); + h->LabelsOption("v"); + } setHistoStyle(hRawYieldsSignal); setHistoStyle(hRawYieldsSignalCounted); @@ -655,6 +667,9 @@ void runMassFitter(const std::string& configFileName) hFitConfig->SetBinContent(ConfigBkgFunc, iSliceVar + 1, bkgFunc[iSliceVar]); hFitConfig->SetBinContent(ConfigSgnFunc, iSliceVar + 1, sgnFunc[iSliceVar]); hFitConfig->SetBinContent(ConfigRandomSeed, iSliceVar + 1, randomSeed); + + hFitResult->SetBinContent(FitResultStatus, iSliceVar + 1, massFitter->getFitStatus()); + hFitResult->SetBinContent(FitResultCovQual, iSliceVar + 1, massFitter->getCovQual()); } // save output histograms @@ -694,6 +709,7 @@ void runMassFitter(const std::string& configFileName) hRawYieldsDscbNR->Write(); } hFitConfig->Write(); + hFitResult->Write(); outputFile.Close(); From 8379accb96c7490be1363394e9dd4056d44e734c Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Tue, 21 Jul 2026 15:40:53 +0200 Subject: [PATCH 18/34] cosmetics of output histograms --- PWGHF/D2H/Macros/runMassFitter.C | 37 ++++++++++++++++---------------- 1 file changed, 19 insertions(+), 18 deletions(-) diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 339aebfb8cd..4ffe4cf3851 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -358,22 +358,22 @@ void runMassFitter(const std::string& configFileName) } // define output histos - auto* hRawYieldsSignal = new TH1D("hRawYieldsSignal", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield", nHistograms, sliceVarLimits.data()); - auto* hRawYieldsSignalCounted = new TH1D("hRawYieldsSignalCounted", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield via bin count", nHistograms, sliceVarLimits.data()); - auto* hRawYieldsBkg = new TH1D("hRawYieldsBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");Background (3#sigma)", nHistograms, sliceVarLimits.data()); - auto* hRawYieldsSgnOverBkg = new TH1D("hRawYieldsSgnOverBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");S/B (3#sigma)", nHistograms, sliceVarLimits.data()); - auto* hRawYieldsSignificance = new TH1D("hRawYieldsSignificance", ";" + sliceVarName + "(" + sliceVarUnit + ");significance (3#sigma)", nHistograms, sliceVarLimits.data()); - auto* hRawYieldsChiSquareBkg = new TH1D("hRawYieldsChiSquareBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");#chi^{2}/#it{ndf}", nHistograms, sliceVarLimits.data()); - auto* hRawYieldsChiSquareTotal = new TH1D("hRawYieldsChiSquareTotal", ";" + sliceVarName + "(" + sliceVarUnit + ");#chi^{2}/#it{ndf}", nHistograms, sliceVarLimits.data()); - auto* hReflectionOverSignal = new TH1D("hReflectionOverSignal", ";" + sliceVarName + "(" + sliceVarUnit + ");Refl/Signal", nHistograms, sliceVarLimits.data()); - auto* hRawYieldsMean = new TH1D("hRawYieldsMean", ";" + sliceVarName + "(" + sliceVarUnit + ");mean (GeV/#it{c}^{2})", nHistograms, sliceVarLimits.data()); - auto* hRawYieldsSigma = new TH1D("hRawYieldsSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nHistograms, sliceVarLimits.data()); - auto* hRawYieldsSecSigma = new TH1D("hRawYieldsSecSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nHistograms, sliceVarLimits.data()); - auto* hRawYieldsFracDoubleGaus = new TH1D("hRawYieldsFracDoubleGaus", ";" + sliceVarName + "(" + sliceVarUnit + ");fraction of double gaussian", nHistograms, sliceVarLimits.data()); - auto* hRawYieldsDscbAlphaL = new TH1D("hRawYieldsDscbAlphaL", ";" + sliceVarName + "(" + sliceVarUnit + ");#alpha_{L}", nHistograms, sliceVarLimits.data()); - auto* hRawYieldsDscbAlphaR = new TH1D("hRawYieldsDscbAlphaR", ";" + sliceVarName + "(" + sliceVarUnit + ");#alpha_{R}", nHistograms, sliceVarLimits.data()); - auto* hRawYieldsDscbNL = new TH1D("hRawYieldsDscbNL", ";" + sliceVarName + "(" + sliceVarUnit + ");n_{L}", nHistograms, sliceVarLimits.data()); - auto* hRawYieldsDscbNR = new TH1D("hRawYieldsDscbNR", ";" + sliceVarName + "(" + sliceVarUnit + ");n_{R}", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsSignal = new TH1D("hRawYieldsSignal", ";" + sliceVarName + " (" + sliceVarUnit + ");raw yield", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsSignalCounted = new TH1D("hRawYieldsSignalCounted", ";" + sliceVarName + " (" + sliceVarUnit + ");raw yield via bin count", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsBkg = new TH1D("hRawYieldsBkg", ";" + sliceVarName + " (" + sliceVarUnit + ");Background (3#sigma)", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsSgnOverBkg = new TH1D("hRawYieldsSgnOverBkg", ";" + sliceVarName + " (" + sliceVarUnit + ");S/B (3#sigma)", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsSignificance = new TH1D("hRawYieldsSignificance", ";" + sliceVarName + " (" + sliceVarUnit + ");significance (3#sigma)", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsChiSquareBkg = new TH1D("hRawYieldsChiSquareBkg", ";" + sliceVarName + " (" + sliceVarUnit + ");#chi^{2}/#it{ndf}", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsChiSquareTotal = new TH1D("hRawYieldsChiSquareTotal", ";" + sliceVarName + " (" + sliceVarUnit + ");#chi^{2}/#it{ndf}", nHistograms, sliceVarLimits.data()); + auto* hReflectionOverSignal = new TH1D("hReflectionOverSignal", ";" + sliceVarName + " (" + sliceVarUnit + ");Refl/Signal", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsMean = new TH1D("hRawYieldsMean", ";" + sliceVarName + " (" + sliceVarUnit + ");mean (GeV/#it{c}^{2})", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsSigma = new TH1D("hRawYieldsSigma", ";" + sliceVarName + " (" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsSecSigma = new TH1D("hRawYieldsSecSigma", ";" + sliceVarName + " (" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsFracDoubleGaus = new TH1D("hRawYieldsFracDoubleGaus", ";" + sliceVarName + " (" + sliceVarUnit + ");fraction of double gaussian", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsDscbAlphaL = new TH1D("hRawYieldsDscbAlphaL", ";" + sliceVarName + " (" + sliceVarUnit + ");#alpha_{L}", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsDscbAlphaR = new TH1D("hRawYieldsDscbAlphaR", ";" + sliceVarName + " (" + sliceVarUnit + ");#alpha_{R}", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsDscbNL = new TH1D("hRawYieldsDscbNL", ";" + sliceVarName + " (" + sliceVarUnit + ");n_{L}", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsDscbNR = new TH1D("hRawYieldsDscbNR", ";" + sliceVarName + " (" + sliceVarUnit + ");n_{R}", nHistograms, sliceVarLimits.data()); enum { ConfigMassMin = 1, @@ -390,9 +390,9 @@ void runMassFitter(const std::string& configFileName) FitResultCovQual, NFitResultsToSave }; - auto* hFitConfig = new TH2F("hfitConfig", "Fit Configurations", NConfigsToSave - 1, 0, NConfigsToSave - 1, nHistograms, sliceVarLimits.data()); + auto* hFitConfig = new TH2F("hFitConfig", "Fit Configurations", NConfigsToSave - 1, 0, NConfigsToSave - 1, nHistograms, sliceVarLimits.data()); const char* hFitConfigXLabel[NConfigsToSave - 1] = {"mass min", "mass max", "rebin num", "fix sigma", "bkg func", "sgn func", "rnd seed"}; - auto* hFitResult = new TH2F("hfitResult", "Fit Result", NFitResultsToSave - 1, 0, NFitResultsToSave - 1, nHistograms, sliceVarLimits.data()); + auto* hFitResult = new TH2F("hFitResult", "Fit Result", NFitResultsToSave - 1, 0, NFitResultsToSave - 1, nHistograms, sliceVarLimits.data()); const char* hFitResultXLabel[NConfigsToSave - 1] = {"status", "cov qual"}; for (int i = 0; i < NConfigsToSave - 1; i++) { hFitConfig->GetXaxis()->SetBinLabel(i + 1, hFitConfigXLabel[i]); @@ -405,6 +405,7 @@ void runMassFitter(const std::string& configFileName) h->LabelsDeflate("X"); h->LabelsDeflate("Y"); h->LabelsOption("v"); + h->GetYaxis()->SetTitle(sliceVarName + " (" + sliceVarUnit + ")"); } setHistoStyle(hRawYieldsSignal); From a3e1784e1564ede3a089d2ed1456b8c8e29049ab Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Tue, 28 Jul 2026 17:34:25 +0200 Subject: [PATCH 19/34] add EDM and minNLL to hFitResult --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 6 +++++- PWGHF/D2H/Macros/HFInvMassFitter.h | 4 ++++ PWGHF/D2H/Macros/runMassFitter.C | 6 +++++- 3 files changed, 14 insertions(+), 2 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 8f2e17d8b6f..307d6f63646 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -156,7 +156,9 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mRandomSeed(randomSeed), mRandomGen(nullptr), mFitStatus(-999), - mCovQual(-999) + mCovQual(-999), + mEdm(-999.), + mMinNll(-999.) { // standard constructor mHistoInvMass = histoToFit; @@ -326,6 +328,8 @@ void HFInvMassFitter::doFit() fitResult->Print("v"); mFitStatus = fitResult->status(); mCovQual = fitResult->covQual(); + mEdm = fitResult->edm(); + mMinNll = fitResult->minNll(); std::cout << "mRooNBkg->getVal() = " << mRooNBkg->getVal() << "\n"; std::cout << "mRooNSgn->getVal() = " << mRooNSgn->getVal() << "\n"; plotBkg(mTotalPdf); diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 434bd5d87fc..43c88b11f9f 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -163,6 +163,8 @@ class HFInvMassFitter : public TNamed [[nodiscard]] double getReflOverSig() const { return mReflPdf != nullptr ? mReflOverSgn : 0.; } [[nodiscard]] int getFitStatus() const { return mFitStatus; } [[nodiscard]] int getCovQual() const { return mCovQual; } + [[nodiscard]] double getEDM() const { return mEdm; } + [[nodiscard]] double getMinNll() const { return mMinNll; } void calculateSignal(double& signal, double& signalErr) const; void countSignal(double& signal, double& signalErr) const; void calculateBackground(double& bkg, double& bkgErr) const; @@ -274,6 +276,8 @@ class HFInvMassFitter : public TNamed TRandom3* mRandomGen; /// engine for fit's initial parameters randomization int mFitStatus; /// fit result status, see https://root-forum.cern.ch/t/meaning-of-values-returned-by-roofitresult-status/16355/2 int mCovQual; /// fit result covariance matrix quality, see https://root.cern.ch/doc/v620/Minuit2Minimizer_8cxx_source.html#l01121 + double mEdm; /// fit quality metrics: Estimated Distance to Minimum + double mMinNll; /// fit quality metrics: minimum negative log-likelihood (NLL) value achieved at the best-fit parameter values ClassDefOverride(HFInvMassFitter, 1); }; diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 4ffe4cf3851..432cd63eba7 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -388,12 +388,14 @@ void runMassFitter(const std::string& configFileName) enum { FitResultStatus = 1, FitResultCovQual, + FitResultEdm, + FitResultMinNll, NFitResultsToSave }; auto* hFitConfig = new TH2F("hFitConfig", "Fit Configurations", NConfigsToSave - 1, 0, NConfigsToSave - 1, nHistograms, sliceVarLimits.data()); const char* hFitConfigXLabel[NConfigsToSave - 1] = {"mass min", "mass max", "rebin num", "fix sigma", "bkg func", "sgn func", "rnd seed"}; auto* hFitResult = new TH2F("hFitResult", "Fit Result", NFitResultsToSave - 1, 0, NFitResultsToSave - 1, nHistograms, sliceVarLimits.data()); - const char* hFitResultXLabel[NConfigsToSave - 1] = {"status", "cov qual"}; + const char* hFitResultXLabel[NConfigsToSave - 1] = {"status", "cov qual", "edm", "minNLL"}; for (int i = 0; i < NConfigsToSave - 1; i++) { hFitConfig->GetXaxis()->SetBinLabel(i + 1, hFitConfigXLabel[i]); } @@ -671,6 +673,8 @@ void runMassFitter(const std::string& configFileName) hFitResult->SetBinContent(FitResultStatus, iSliceVar + 1, massFitter->getFitStatus()); hFitResult->SetBinContent(FitResultCovQual, iSliceVar + 1, massFitter->getCovQual()); + hFitResult->SetBinContent(FitResultEdm, iSliceVar + 1, massFitter->getEDM()); + hFitResult->SetBinContent(FitResultMinNll, iSliceVar + 1, massFitter->getMinNll()); } // save output histograms From 8d283094acb6689902ba556660ec12c2f77e5eea Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 29 Jul 2026 12:36:08 +0200 Subject: [PATCH 20/34] add mNSgn global correlation coefficient to hFitResult --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 5 ++++- PWGHF/D2H/Macros/HFInvMassFitter.h | 2 ++ PWGHF/D2H/Macros/runMassFitter.C | 4 +++- 3 files changed, 9 insertions(+), 2 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 307d6f63646..728e7f69ae8 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -158,7 +158,8 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mFitStatus(-999), mCovQual(-999), mEdm(-999.), - mMinNll(-999.) + mMinNll(-999.), + mSgnGlobalCorrelCoeff(-999.) { // standard constructor mHistoInvMass = histoToFit; @@ -325,11 +326,13 @@ void HFInvMassFitter::doFit() std::cout << "CovQual = " << fitResult->covQual() << "\n"; std::cout << "EDM = " << fitResult->edm() << "\n"; std::cout << "minNLL = " << fitResult->minNll() << "\n"; + std::cout << "Global correlation of mNSgn = " << fitResult->globalCorr("mNSgn") << "\n"; fitResult->Print("v"); mFitStatus = fitResult->status(); mCovQual = fitResult->covQual(); mEdm = fitResult->edm(); mMinNll = fitResult->minNll(); + mSgnGlobalCorrelCoeff = fitResult->globalCorr("mNSgn"); std::cout << "mRooNBkg->getVal() = " << mRooNBkg->getVal() << "\n"; std::cout << "mRooNSgn->getVal() = " << mRooNSgn->getVal() << "\n"; plotBkg(mTotalPdf); diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 43c88b11f9f..42fa5310074 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -165,6 +165,7 @@ class HFInvMassFitter : public TNamed [[nodiscard]] int getCovQual() const { return mCovQual; } [[nodiscard]] double getEDM() const { return mEdm; } [[nodiscard]] double getMinNll() const { return mMinNll; } + [[nodiscard]] double getSgnGlobalCorrelCoeff() const { return mSgnGlobalCorrelCoeff; } void calculateSignal(double& signal, double& signalErr) const; void countSignal(double& signal, double& signalErr) const; void calculateBackground(double& bkg, double& bkgErr) const; @@ -278,6 +279,7 @@ class HFInvMassFitter : public TNamed int mCovQual; /// fit result covariance matrix quality, see https://root.cern.ch/doc/v620/Minuit2Minimizer_8cxx_source.html#l01121 double mEdm; /// fit quality metrics: Estimated Distance to Minimum double mMinNll; /// fit quality metrics: minimum negative log-likelihood (NLL) value achieved at the best-fit parameter values + double mSgnGlobalCorrelCoeff; /// global correlation coefficient of mRooNSgn with other fit parameters ClassDefOverride(HFInvMassFitter, 1); }; diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 432cd63eba7..4feb15c1a0d 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -390,12 +390,13 @@ void runMassFitter(const std::string& configFileName) FitResultCovQual, FitResultEdm, FitResultMinNll, + FitResultNSgnGlobalCorrelCoeff, NFitResultsToSave }; auto* hFitConfig = new TH2F("hFitConfig", "Fit Configurations", NConfigsToSave - 1, 0, NConfigsToSave - 1, nHistograms, sliceVarLimits.data()); const char* hFitConfigXLabel[NConfigsToSave - 1] = {"mass min", "mass max", "rebin num", "fix sigma", "bkg func", "sgn func", "rnd seed"}; auto* hFitResult = new TH2F("hFitResult", "Fit Result", NFitResultsToSave - 1, 0, NFitResultsToSave - 1, nHistograms, sliceVarLimits.data()); - const char* hFitResultXLabel[NConfigsToSave - 1] = {"status", "cov qual", "edm", "minNLL"}; + const char* hFitResultXLabel[NConfigsToSave - 1] = {"status", "cov qual", "edm", "minNLL", "N sig GCC"}; for (int i = 0; i < NConfigsToSave - 1; i++) { hFitConfig->GetXaxis()->SetBinLabel(i + 1, hFitConfigXLabel[i]); } @@ -675,6 +676,7 @@ void runMassFitter(const std::string& configFileName) hFitResult->SetBinContent(FitResultCovQual, iSliceVar + 1, massFitter->getCovQual()); hFitResult->SetBinContent(FitResultEdm, iSliceVar + 1, massFitter->getEDM()); hFitResult->SetBinContent(FitResultMinNll, iSliceVar + 1, massFitter->getMinNll()); + hFitResult->SetBinContent(FitResultNSgnGlobalCorrelCoeff, iSliceVar + 1, massFitter->getSgnGlobalCorrelCoeff()); } // save output histograms From 8fc7222325372d2ab96677c6d2e7403f3344d9e2 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 29 Jul 2026 18:15:35 +0200 Subject: [PATCH 21/34] add pt-wise hSgnCorr histograms --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 17 ++- PWGHF/D2H/Macros/HFInvMassFitter.h | 188 ++++++++++++++------------- PWGHF/D2H/Macros/runMassFitter.C | 16 ++- 3 files changed, 124 insertions(+), 97 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 728e7f69ae8..504e6322e91 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -159,7 +159,8 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mCovQual(-999), mEdm(-999.), mMinNll(-999.), - mSgnGlobalCorrelCoeff(-999.) + mSgnCorrelCoeffValues({}), + mSgnCorrelCoeffNames({}) { // standard constructor mHistoInvMass = histoToFit; @@ -321,6 +322,7 @@ void HFInvMassFitter::doFit() } else { fitResult = mTotalPdf->fitTo(dataHistogram, Save()); } + std::cout << "Finish total fit\n"; std::cout << "Status = " << fitResult->status() << "\n"; std::cout << "CovQual = " << fitResult->covQual() << "\n"; @@ -332,7 +334,18 @@ void HFInvMassFitter::doFit() mCovQual = fitResult->covQual(); mEdm = fitResult->edm(); mMinNll = fitResult->minNll(); - mSgnGlobalCorrelCoeff = fitResult->globalCorr("mNSgn"); + + mSgnCorrelCoeffValues.push_back(fitResult->globalCorr("mNSgn")); + mSgnCorrelCoeffNames.push_back("global"); + mSgnCorrelCoeffValues.push_back(fitResult->correlation("mNSgn", "mRooNBkg")); + mSgnCorrelCoeffNames.push_back("mRooNBkg"); + const auto& bkgPars = bkgPdf->getParameters(dataHistogram)->selectByAttrib("Constant", false); + for (const auto& bkgPar : *bkgPars) { + const std::string& bkgParName = bkgPar->GetName(); + mSgnCorrelCoeffValues.push_back(fitResult->correlation("mNSgn", bkgParName.c_str())); + mSgnCorrelCoeffNames.push_back(bkgParName); + } + std::cout << "mRooNBkg->getVal() = " << mRooNBkg->getVal() << "\n"; std::cout << "mRooNSgn->getVal() = " << mRooNSgn->getVal() << "\n"; plotBkg(mTotalPdf); diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 42fa5310074..82f8ddbdde6 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -22,6 +22,8 @@ #ifndef PWGHF_D2H_MACROS_HFINVMASSFITTER_H_ #define PWGHF_D2H_MACROS_HFINVMASSFITTER_H_ +#include +#include #include #include #include @@ -165,7 +167,8 @@ class HFInvMassFitter : public TNamed [[nodiscard]] int getCovQual() const { return mCovQual; } [[nodiscard]] double getEDM() const { return mEdm; } [[nodiscard]] double getMinNll() const { return mMinNll; } - [[nodiscard]] double getSgnGlobalCorrelCoeff() const { return mSgnGlobalCorrelCoeff; } + [[nodiscard]] const std::vector& getSgnCorrelCoeffValues() const { return mSgnCorrelCoeffValues; } + [[nodiscard]] const std::vector& getSgnCorrelCoeffNames() const { return mSgnCorrelCoeffNames; } void calculateSignal(double& signal, double& signalErr) const; void countSignal(double& signal, double& signalErr) const; void calculateBackground(double& bkg, double& bkgErr) const; @@ -189,97 +192,98 @@ class HFInvMassFitter : public TNamed TH1* mHistoInvMass; // histogram to fit std::string mFitOption; - double mMinMass; // lower mass limit - double mMaxMass; // upper mass limit - int mTypeOfBkgPdf; // background fit function - int mTypeOfSgnPdf; // signal fit function - int mTypeOfReflPdf; // reflection fit function - double mMassParticle; // pdg value of particle mass - double mMass; /// signal gaussian mean value - double mMassLowLimit; /// lower limit of the allowed mass range - double mMassUpLimit; /// upper limit of the allowed mass range - double mMassReflLowLimit; /// lower limit of the allowed mass range for reflection - double mMassReflUpLimit; /// upper limit of the allowed mass range for reflection - double mSecMass; /// Second peak mean value - double mSigmaSgn; /// signal gaussian sigma - double mSecSigma; /// Second peak gaussian sigma - double mNSigmaForSidebands; /// number of sigmas to veto the signal peak - double mNSigmaForSgn; /// number of sigmas to veto the signal peak - double mSigmaSgnErr; /// uncertainty on signal gaussian sigma - double mSigmaSgnDoubleGaus; /// signal 2gaussian sigma - bool mFixedMean; /// switch for fix mean of gaussian - bool mBoundMean; /// switch for bound mean of guassian - bool mBoundReflMean; /// switch for bound mean of guassian for reflection - bool mFixedSigma; /// fix sigma or not - bool mFixedSigmaDoubleGaus; /// fix sigma of 2gaussian or not - bool mBoundSigma; /// set bound sigma or not - bool mFixedDscbTailParams; /// switch for fix double sided Crystal Ball tail parameters - double mSigmaValue; /// value of sigma - double mParamSgn; /// +/- range variation of bound Sigma of gaussian in % - double mFracDoubleGaus; /// initialization for fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar - double mFixedRawYield; /// initialization for raw yield - bool mFixedFracDoubleGaus; /// switch for fixed fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar - double mRatioDoubleGausSigma; /// initialization for ratio between two gaussian sigmas in case of k2GausSigmaRatioPar - bool mFixedRatioDoubleGausSigma; /// switch for fixed ratio between two gaussian sigmas in case of k2GausSigmaRatioPar - double mReflOverSgn; /// reflection/signal - bool mEnableReflections; /// flag use/not use reflections - double mRawYield; /// signal gaussian integral - double mRawYieldErr; /// err on signal gaussian integral - double mRawYieldCounted; /// signal gaussian integral evaluated via bin counting - double mRawYieldCountedErr; /// err on signal gaussian integral evaluated via bin counting - double mBkgYield; /// background - double mBkgYieldErr; /// err on background - double mSignificance; /// significance - double mSignificanceErr; /// err on significance - double mChiSquareOverNdfTotal; /// chi2/ndf of the total fit - double mChiSquareOverNdfBkg; /// chi2/ndf of the background (sidebands) pre-fit - bool mFixReflOverSgn; /// switch for fix refl/signal - double mDscbAlphaLInitialValue; /// double sided Crystal Ball alpha left initial value - double mDscbAlphaLLowLimit; /// double sided Crystal Ball alpha left lower limit - double mDscbAlphaLUpLimit; /// double sided Crystal Ball alpha left upper limit - double mDscbAlphaRInitialValue; /// double sided Crystal Ball alpha right initial value - double mDscbAlphaRLowLimit; /// double sided Crystal Ball alpha right lower limit - double mDscbAlphaRUpLimit; /// double sided Crystal Ball alpha right upper limit - double mDscbNLInitialValue; /// double sided Crystal Ball n left initial value - double mDscbNLLowLimit; /// double sided Crystal Ball n left lower limit - double mDscbNLUpLimit; /// double sided Crystal Ball n left upper limit - double mDscbNRInitialValue; /// double sided Crystal Ball n right initial value - double mDscbNRLowLimit; /// double sided Crystal Ball n right lower limit - double mDscbNRUpLimit; /// double sided Crystal Ball n right upper limit - RooRealVar* mRooMeanSgn; /// mean for gaussian of signal - RooRealVar* mRooSigmaSgn; /// sigma for gaussian of signal - RooRealVar* mRooSecSigmaSgn; /// second sigma for composite gaussian of signal - RooRealVar* mRooFracDoubleGaus; /// fraction of second gaussian for composite gaussian of signal - RooAbsPdf* mSgnPdf; /// signal fit function - RooAbsPdf* mBkgPdf; /// background fit function - RooAbsPdf* mReflPdf; /// reflection fit function - RooRealVar* mRooNSgn; /// total Signal fit function integral - RooRealVar* mRooNBkg; /// total background fit function integral - RooRealVar* mRooNRefl; /// total reflection fit function integral - RooRealVar* mRooDscbAlphaL; /// double sided Crystal Ball alpha left - RooRealVar* mRooDscbAlphaR; /// double sided Crystal Ball alpha right - RooRealVar* mRooDscbNL; /// double sided Crystal Ball n left - RooRealVar* mRooDscbNR; /// double sided Crystal Ball n right - RooAbsPdf* mTotalPdf; /// total fit function - RooPlot* mInvMassFrame; /// frame of mass - RooPlot* mReflFrame; /// reflection frame - RooPlot* mReflOnlyFrame; /// reflection frame plot on reflection only - RooPlot* mResidualFrame; /// residual frame - RooHist* mResidualHist; /// residual histogram - RooPlot* mRatioFrame; /// fit/data ratio frame - RooWorkspace* mWorkspace; /// workspace - double mIntegralBkg; /// integral of background fit function - double mIntegralSgn; /// integral of signal fit function - TH1* mHistoTemplateRefl; /// reflection histogram - bool mDrawBgPrefit; /// draw background after fitting the sidebands - bool mHighlightPeakRegion; /// draw vertical lines showing the peak region (usually +- 3 sigma) - int mRandomSeed; /// seed for random engine for fit's initial parameters randomization - TRandom3* mRandomGen; /// engine for fit's initial parameters randomization - int mFitStatus; /// fit result status, see https://root-forum.cern.ch/t/meaning-of-values-returned-by-roofitresult-status/16355/2 - int mCovQual; /// fit result covariance matrix quality, see https://root.cern.ch/doc/v620/Minuit2Minimizer_8cxx_source.html#l01121 - double mEdm; /// fit quality metrics: Estimated Distance to Minimum - double mMinNll; /// fit quality metrics: minimum negative log-likelihood (NLL) value achieved at the best-fit parameter values - double mSgnGlobalCorrelCoeff; /// global correlation coefficient of mRooNSgn with other fit parameters + double mMinMass; // lower mass limit + double mMaxMass; // upper mass limit + int mTypeOfBkgPdf; // background fit function + int mTypeOfSgnPdf; // signal fit function + int mTypeOfReflPdf; // reflection fit function + double mMassParticle; // pdg value of particle mass + double mMass; /// signal gaussian mean value + double mMassLowLimit; /// lower limit of the allowed mass range + double mMassUpLimit; /// upper limit of the allowed mass range + double mMassReflLowLimit; /// lower limit of the allowed mass range for reflection + double mMassReflUpLimit; /// upper limit of the allowed mass range for reflection + double mSecMass; /// Second peak mean value + double mSigmaSgn; /// signal gaussian sigma + double mSecSigma; /// Second peak gaussian sigma + double mNSigmaForSidebands; /// number of sigmas to veto the signal peak + double mNSigmaForSgn; /// number of sigmas to veto the signal peak + double mSigmaSgnErr; /// uncertainty on signal gaussian sigma + double mSigmaSgnDoubleGaus; /// signal 2gaussian sigma + bool mFixedMean; /// switch for fix mean of gaussian + bool mBoundMean; /// switch for bound mean of guassian + bool mBoundReflMean; /// switch for bound mean of guassian for reflection + bool mFixedSigma; /// fix sigma or not + bool mFixedSigmaDoubleGaus; /// fix sigma of 2gaussian or not + bool mBoundSigma; /// set bound sigma or not + bool mFixedDscbTailParams; /// switch for fix double sided Crystal Ball tail parameters + double mSigmaValue; /// value of sigma + double mParamSgn; /// +/- range variation of bound Sigma of gaussian in % + double mFracDoubleGaus; /// initialization for fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar + double mFixedRawYield; /// initialization for raw yield + bool mFixedFracDoubleGaus; /// switch for fixed fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar + double mRatioDoubleGausSigma; /// initialization for ratio between two gaussian sigmas in case of k2GausSigmaRatioPar + bool mFixedRatioDoubleGausSigma; /// switch for fixed ratio between two gaussian sigmas in case of k2GausSigmaRatioPar + double mReflOverSgn; /// reflection/signal + bool mEnableReflections; /// flag use/not use reflections + double mRawYield; /// signal gaussian integral + double mRawYieldErr; /// err on signal gaussian integral + double mRawYieldCounted; /// signal gaussian integral evaluated via bin counting + double mRawYieldCountedErr; /// err on signal gaussian integral evaluated via bin counting + double mBkgYield; /// background + double mBkgYieldErr; /// err on background + double mSignificance; /// significance + double mSignificanceErr; /// err on significance + double mChiSquareOverNdfTotal; /// chi2/ndf of the total fit + double mChiSquareOverNdfBkg; /// chi2/ndf of the background (sidebands) pre-fit + bool mFixReflOverSgn; /// switch for fix refl/signal + double mDscbAlphaLInitialValue; /// double sided Crystal Ball alpha left initial value + double mDscbAlphaLLowLimit; /// double sided Crystal Ball alpha left lower limit + double mDscbAlphaLUpLimit; /// double sided Crystal Ball alpha left upper limit + double mDscbAlphaRInitialValue; /// double sided Crystal Ball alpha right initial value + double mDscbAlphaRLowLimit; /// double sided Crystal Ball alpha right lower limit + double mDscbAlphaRUpLimit; /// double sided Crystal Ball alpha right upper limit + double mDscbNLInitialValue; /// double sided Crystal Ball n left initial value + double mDscbNLLowLimit; /// double sided Crystal Ball n left lower limit + double mDscbNLUpLimit; /// double sided Crystal Ball n left upper limit + double mDscbNRInitialValue; /// double sided Crystal Ball n right initial value + double mDscbNRLowLimit; /// double sided Crystal Ball n right lower limit + double mDscbNRUpLimit; /// double sided Crystal Ball n right upper limit + RooRealVar* mRooMeanSgn; /// mean for gaussian of signal + RooRealVar* mRooSigmaSgn; /// sigma for gaussian of signal + RooRealVar* mRooSecSigmaSgn; /// second sigma for composite gaussian of signal + RooRealVar* mRooFracDoubleGaus; /// fraction of second gaussian for composite gaussian of signal + RooAbsPdf* mSgnPdf; /// signal fit function + RooAbsPdf* mBkgPdf; /// background fit function + RooAbsPdf* mReflPdf; /// reflection fit function + RooRealVar* mRooNSgn; /// total Signal fit function integral + RooRealVar* mRooNBkg; /// total background fit function integral + RooRealVar* mRooNRefl; /// total reflection fit function integral + RooRealVar* mRooDscbAlphaL; /// double sided Crystal Ball alpha left + RooRealVar* mRooDscbAlphaR; /// double sided Crystal Ball alpha right + RooRealVar* mRooDscbNL; /// double sided Crystal Ball n left + RooRealVar* mRooDscbNR; /// double sided Crystal Ball n right + RooAbsPdf* mTotalPdf; /// total fit function + RooPlot* mInvMassFrame; /// frame of mass + RooPlot* mReflFrame; /// reflection frame + RooPlot* mReflOnlyFrame; /// reflection frame plot on reflection only + RooPlot* mResidualFrame; /// residual frame + RooHist* mResidualHist; /// residual histogram + RooPlot* mRatioFrame; /// fit/data ratio frame + RooWorkspace* mWorkspace; /// workspace + double mIntegralBkg; /// integral of background fit function + double mIntegralSgn; /// integral of signal fit function + TH1* mHistoTemplateRefl; /// reflection histogram + bool mDrawBgPrefit; /// draw background after fitting the sidebands + bool mHighlightPeakRegion; /// draw vertical lines showing the peak region (usually +- 3 sigma) + int mRandomSeed; /// seed for random engine for fit's initial parameters randomization + TRandom3* mRandomGen; /// engine for fit's initial parameters randomization + int mFitStatus; /// fit result status, see https://root-forum.cern.ch/t/meaning-of-values-returned-by-roofitresult-status/16355/2 + int mCovQual; /// fit result covariance matrix quality, see https://root.cern.ch/doc/v620/Minuit2Minimizer_8cxx_source.html#l01121 + double mEdm; /// fit quality metrics: Estimated Distance to Minimum + double mMinNll; /// fit quality metrics: minimum negative log-likelihood (NLL) value achieved at the best-fit parameter values + std::vector mSgnCorrelCoeffValues; /// correlation coefficients of mRooNSgn (global and with BG fit parameters) values + std::vector mSgnCorrelCoeffNames; /// correlation coefficients of mRooNSgn (global and with BG fit parameters) names ClassDefOverride(HFInvMassFitter, 1); }; diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 4feb15c1a0d..71c04b1d8ff 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -321,6 +321,7 @@ void runMassFitter(const std::string& configFileName) std::vector hMassSgn(nHistograms); std::vector hMassRefl(nHistograms); std::vector hMass(nHistograms); + std::vector hSgnCorr(nHistograms); for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { if (!isMc) { @@ -390,13 +391,12 @@ void runMassFitter(const std::string& configFileName) FitResultCovQual, FitResultEdm, FitResultMinNll, - FitResultNSgnGlobalCorrelCoeff, NFitResultsToSave }; auto* hFitConfig = new TH2F("hFitConfig", "Fit Configurations", NConfigsToSave - 1, 0, NConfigsToSave - 1, nHistograms, sliceVarLimits.data()); const char* hFitConfigXLabel[NConfigsToSave - 1] = {"mass min", "mass max", "rebin num", "fix sigma", "bkg func", "sgn func", "rnd seed"}; auto* hFitResult = new TH2F("hFitResult", "Fit Result", NFitResultsToSave - 1, 0, NFitResultsToSave - 1, nHistograms, sliceVarLimits.data()); - const char* hFitResultXLabel[NConfigsToSave - 1] = {"status", "cov qual", "edm", "minNLL", "N sig GCC"}; + const char* hFitResultXLabel[NFitResultsToSave - 1] = {"status", "cov qual", "edm", "minNLL"}; for (int i = 0; i < NConfigsToSave - 1; i++) { hFitConfig->GetXaxis()->SetBinLabel(i + 1, hFitConfigXLabel[i]); } @@ -676,7 +676,16 @@ void runMassFitter(const std::string& configFileName) hFitResult->SetBinContent(FitResultCovQual, iSliceVar + 1, massFitter->getCovQual()); hFitResult->SetBinContent(FitResultEdm, iSliceVar + 1, massFitter->getEDM()); hFitResult->SetBinContent(FitResultMinNll, iSliceVar + 1, massFitter->getMinNll()); - hFitResult->SetBinContent(FitResultNSgnGlobalCorrelCoeff, iSliceVar + 1, massFitter->getSgnGlobalCorrelCoeff()); + + const std::string hSgnCorrName = "hSgnCorr" + std::to_string(iSliceVar + 1); + const auto& sgnCorrValues = massFitter->getSgnCorrelCoeffValues(); + const auto& sgnCorrNames = massFitter->getSgnCorrelCoeffNames(); + const int nSgnCorrValues = sgnCorrValues.size(); + hSgnCorr[iSliceVar] = new TH1D(hSgnCorrName.c_str(), hSgnCorrName.c_str(), nSgnCorrValues, 0, nSgnCorrValues); + for (int iSgnCorrValue = 0; iSgnCorrValue < nSgnCorrValues; ++iSgnCorrValue) { + hSgnCorr[iSliceVar]->SetBinContent(iSgnCorrValue + 1, sgnCorrValues.at(iSgnCorrValue)); + hSgnCorr[iSliceVar]->GetXaxis()->SetBinLabel(iSgnCorrValue + 1, sgnCorrNames.at(iSgnCorrValue).c_str()); + } } // save output histograms @@ -694,6 +703,7 @@ void runMassFitter(const std::string& configFileName) for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { hMass[iSliceVar]->Write(); + hSgnCorr[iSliceVar]->Write(); } hRawYieldsSignal->Write(); hRawYieldsSignalCounted->Write(); From 899e36b873d1d2758ceec3854e0289015e9cf97d Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Thu, 30 Jul 2026 10:48:14 +0200 Subject: [PATCH 22/34] save pt-wise histograms in dedicated TDirectory-s --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 1 - PWGHF/D2H/Macros/runMassFitter.C | 9 ++++++++- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 504e6322e91..556074749d4 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -164,7 +164,6 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, { // standard constructor mHistoInvMass = histoToFit; - mHistoInvMass->SetName("mHistoInvMass"); mHistoInvMass->SetDirectory(nullptr); if (mRandomSeed >= 0) { mRandomGen = new TRandom3(); diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 71c04b1d8ff..43949c22ba3 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -490,7 +490,7 @@ void runMassFitter(const std::string& configFileName) hMass[iSliceVar]->SetTitle(Form("%s;%s;Counts per %0.1f MeV/#it{c}^{2}", ptTitle.Data(), massAxisTitle.c_str(), hMass[iSliceVar]->GetBinWidth(1) * 1000)); - hMass[iSliceVar]->SetName(Form("MassForFit%d", iSliceVar)); + hMass[iSliceVar]->SetName(Form("hMassForFit%d", iSliceVar + 1)); if (enableRefl) { hMassRefl[iSliceVar]->Rebin(nRebin[iSliceVar]); @@ -702,9 +702,16 @@ void runMassFitter(const std::string& configFileName) } for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { + if (iSliceVar == 0) { + outputFile.mkdir("MassHistograms"); + outputFile.mkdir("SgnCorrHistograms"); + } + outputFile.cd("MassHistograms"); hMass[iSliceVar]->Write(); + outputFile.cd("SgnCorrHistograms"); hSgnCorr[iSliceVar]->Write(); } + outputFile.cd(); hRawYieldsSignal->Write(); hRawYieldsSignalCounted->Write(); hRawYieldsBkg->Write(); From b7f842d0d43fe82904824e75fbbb948f847f4160 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Thu, 30 Jul 2026 19:08:06 +0200 Subject: [PATCH 23/34] print more debug info (cov, corr matrix; status history) --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 556074749d4..473e6d9f17a 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -323,12 +323,21 @@ void HFInvMassFitter::doFit() } std::cout << "Finish total fit\n"; + + std::cout << "History:\t"; + for (unsigned i = 0; i < fitResult->numStatusHistory(); ++i) { + std::cout << fitResult->statusLabelHistory(i) << " : " << fitResult->statusCodeHistory(i) << "\t"; + } + std::cout << "\n"; + std::cout << "Status = " << fitResult->status() << "\n"; std::cout << "CovQual = " << fitResult->covQual() << "\n"; std::cout << "EDM = " << fitResult->edm() << "\n"; std::cout << "minNLL = " << fitResult->minNll() << "\n"; std::cout << "Global correlation of mNSgn = " << fitResult->globalCorr("mNSgn") << "\n"; fitResult->Print("v"); + fitResult->covarianceMatrix().Print(); + fitResult->correlationMatrix().Print(); mFitStatus = fitResult->status(); mCovQual = fitResult->covQual(); mEdm = fitResult->edm(); @@ -347,6 +356,12 @@ void HFInvMassFitter::doFit() std::cout << "mRooNBkg->getVal() = " << mRooNBkg->getVal() << "\n"; std::cout << "mRooNSgn->getVal() = " << mRooNSgn->getVal() << "\n"; + + std::cout << "Value = " << mRooNSgn->getVal() << "\n"; + std::cout << "HESSE error = " << mRooNSgn->getError() << "\n"; + std::cout << "MINOS lower = " << mRooNSgn->getAsymErrorLo() << "\n"; + std::cout << "MINOS upper = " << mRooNSgn->getAsymErrorHi() << "\n"; + plotBkg(mTotalPdf); mTotalPdf->plotOn(mInvMassFrame, Name("Tot_c"), LineColor(kBlue)); if (mHistoTemplateRefl != nullptr) { From d5b5bc4bee55dd67a8dcdf46aa0fd5a328cd9636 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Fri, 31 Jul 2026 13:47:01 +0200 Subject: [PATCH 24/34] save cov and corr matrix into dedicated TH2 --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 31 +++++++++++++++++++++++++++- PWGHF/D2H/Macros/HFInvMassFitter.h | 5 +++++ PWGHF/D2H/Macros/runMassFitter.C | 4 ++++ 3 files changed, 39 insertions(+), 1 deletion(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 473e6d9f17a..744aec9f759 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -40,6 +40,7 @@ #include #include #include +#include #include #include #include @@ -160,7 +161,8 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mEdm(-999.), mMinNll(-999.), mSgnCorrelCoeffValues({}), - mSgnCorrelCoeffNames({}) + mSgnCorrelCoeffNames({}), + mCovCorrMatrix(nullptr) { // standard constructor mHistoInvMass = histoToFit; @@ -343,6 +345,8 @@ void HFInvMassFitter::doFit() mEdm = fitResult->edm(); mMinNll = fitResult->minNll(); + mCovCorrMatrix = fillCovCorrMatrix(fitResult); + mSgnCorrelCoeffValues.push_back(fitResult->globalCorr("mNSgn")); mSgnCorrelCoeffNames.push_back("global"); mSgnCorrelCoeffValues.push_back(fitResult->correlation("mNSgn", "mRooNBkg")); @@ -1256,3 +1260,28 @@ void HFInvMassFitter::cutRangesFromHisto(TH1* histo, const std::vectorfloatParsFinal(); + const int nPars = pars.size(); + + TH2* hMatrix = new TH2D("covCorrMatrix", "covariance (upper left + diagonal) and correlation (lower right) matrix", nPars, 0, nPars, nPars, 0, nPars); + for (int iPar = 0; iPar < nPars; ++iPar) { + hMatrix->GetXaxis()->SetBinLabel(iPar + 1, pars[iPar].GetName()); + hMatrix->GetYaxis()->SetBinLabel(iPar + 1, pars[iPar].GetName()); + } + + const TMatrixDSym& covMatrix = fitResult->covarianceMatrix(); + const TMatrixDSym& corrMatrix = fitResult->correlationMatrix(); + + for (int iPar = 0; iPar < nPars; ++iPar) { + hMatrix->SetBinContent(iPar + 1, iPar + 1, covMatrix(iPar, iPar)); + for (int jPar = iPar + 1; jPar < nPars; ++jPar) { + hMatrix->SetBinContent(iPar + 1, jPar + 1, covMatrix(iPar, jPar)); + hMatrix->SetBinContent(jPar + 1, iPar + 1, corrMatrix(iPar, jPar)); + } + } + + return hMatrix; +} diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 82f8ddbdde6..abbfddee4f7 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -24,11 +24,13 @@ #include #include +#include #include #include #include #include #include +#include #include #include #include @@ -169,6 +171,7 @@ class HFInvMassFitter : public TNamed [[nodiscard]] double getMinNll() const { return mMinNll; } [[nodiscard]] const std::vector& getSgnCorrelCoeffValues() const { return mSgnCorrelCoeffValues; } [[nodiscard]] const std::vector& getSgnCorrelCoeffNames() const { return mSgnCorrelCoeffNames; } + [[nodiscard]] TH2* getCovCorrMatrix() const { return mCovCorrMatrix; } void calculateSignal(double& signal, double& signalErr) const; void countSignal(double& signal, double& signalErr) const; void calculateBackground(double& bkg, double& bkgErr) const; @@ -189,6 +192,7 @@ class HFInvMassFitter : public TNamed [[nodiscard]] std::pair getRangesOfSignal() const; [[nodiscard]] double integrateHistoInvMassOverWorkspaceRanges(const std::vector& ranges) const; void cutRangesFromHisto(TH1* histo, const std::vector& ranges) const; + static TH2* fillCovCorrMatrix(const RooFitResult* fitResult); TH1* mHistoInvMass; // histogram to fit std::string mFitOption; @@ -284,6 +288,7 @@ class HFInvMassFitter : public TNamed double mMinNll; /// fit quality metrics: minimum negative log-likelihood (NLL) value achieved at the best-fit parameter values std::vector mSgnCorrelCoeffValues; /// correlation coefficients of mRooNSgn (global and with BG fit parameters) values std::vector mSgnCorrelCoeffNames; /// correlation coefficients of mRooNSgn (global and with BG fit parameters) names + TH2* mCovCorrMatrix; /// covariance (upper left + diagonal) and correlation (lower right) matrix of free fit parameters ClassDefOverride(HFInvMassFitter, 1); }; diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 43949c22ba3..a70f5caa3fc 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -322,6 +322,7 @@ void runMassFitter(const std::string& configFileName) std::vector hMassRefl(nHistograms); std::vector hMass(nHistograms); std::vector hSgnCorr(nHistograms); + std::vector hCovCorr(nHistograms); for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { if (!isMc) { @@ -686,6 +687,8 @@ void runMassFitter(const std::string& configFileName) hSgnCorr[iSliceVar]->SetBinContent(iSgnCorrValue + 1, sgnCorrValues.at(iSgnCorrValue)); hSgnCorr[iSliceVar]->GetXaxis()->SetBinLabel(iSgnCorrValue + 1, sgnCorrNames.at(iSgnCorrValue).c_str()); } + + hCovCorr[iSliceVar] = massFitter->getCovCorrMatrix(); } // save output histograms @@ -710,6 +713,7 @@ void runMassFitter(const std::string& configFileName) hMass[iSliceVar]->Write(); outputFile.cd("SgnCorrHistograms"); hSgnCorr[iSliceVar]->Write(); + hCovCorr[iSliceVar]->Write(Form("hCovCorrMatrix%d", iSliceVar + 1)); } outputFile.cd(); hRawYieldsSignal->Write(); From 07fbabf508d6e37bd7e11ff2d4163fac0466c8bf Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Fri, 31 Jul 2026 13:53:01 +0200 Subject: [PATCH 25/34] rm sgn corr 1D histograms --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 13 -- PWGHF/D2H/Macros/HFInvMassFitter.h | 186 +++++++++++++-------------- PWGHF/D2H/Macros/runMassFitter.C | 14 +- 3 files changed, 92 insertions(+), 121 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 744aec9f759..b92919a8714 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -160,8 +160,6 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mCovQual(-999), mEdm(-999.), mMinNll(-999.), - mSgnCorrelCoeffValues({}), - mSgnCorrelCoeffNames({}), mCovCorrMatrix(nullptr) { // standard constructor @@ -347,17 +345,6 @@ void HFInvMassFitter::doFit() mCovCorrMatrix = fillCovCorrMatrix(fitResult); - mSgnCorrelCoeffValues.push_back(fitResult->globalCorr("mNSgn")); - mSgnCorrelCoeffNames.push_back("global"); - mSgnCorrelCoeffValues.push_back(fitResult->correlation("mNSgn", "mRooNBkg")); - mSgnCorrelCoeffNames.push_back("mRooNBkg"); - const auto& bkgPars = bkgPdf->getParameters(dataHistogram)->selectByAttrib("Constant", false); - for (const auto& bkgPar : *bkgPars) { - const std::string& bkgParName = bkgPar->GetName(); - mSgnCorrelCoeffValues.push_back(fitResult->correlation("mNSgn", bkgParName.c_str())); - mSgnCorrelCoeffNames.push_back(bkgParName); - } - std::cout << "mRooNBkg->getVal() = " << mRooNBkg->getVal() << "\n"; std::cout << "mRooNSgn->getVal() = " << mRooNSgn->getVal() << "\n"; diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index abbfddee4f7..6fe6538f6cc 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -169,8 +169,6 @@ class HFInvMassFitter : public TNamed [[nodiscard]] int getCovQual() const { return mCovQual; } [[nodiscard]] double getEDM() const { return mEdm; } [[nodiscard]] double getMinNll() const { return mMinNll; } - [[nodiscard]] const std::vector& getSgnCorrelCoeffValues() const { return mSgnCorrelCoeffValues; } - [[nodiscard]] const std::vector& getSgnCorrelCoeffNames() const { return mSgnCorrelCoeffNames; } [[nodiscard]] TH2* getCovCorrMatrix() const { return mCovCorrMatrix; } void calculateSignal(double& signal, double& signalErr) const; void countSignal(double& signal, double& signalErr) const; @@ -196,99 +194,97 @@ class HFInvMassFitter : public TNamed TH1* mHistoInvMass; // histogram to fit std::string mFitOption; - double mMinMass; // lower mass limit - double mMaxMass; // upper mass limit - int mTypeOfBkgPdf; // background fit function - int mTypeOfSgnPdf; // signal fit function - int mTypeOfReflPdf; // reflection fit function - double mMassParticle; // pdg value of particle mass - double mMass; /// signal gaussian mean value - double mMassLowLimit; /// lower limit of the allowed mass range - double mMassUpLimit; /// upper limit of the allowed mass range - double mMassReflLowLimit; /// lower limit of the allowed mass range for reflection - double mMassReflUpLimit; /// upper limit of the allowed mass range for reflection - double mSecMass; /// Second peak mean value - double mSigmaSgn; /// signal gaussian sigma - double mSecSigma; /// Second peak gaussian sigma - double mNSigmaForSidebands; /// number of sigmas to veto the signal peak - double mNSigmaForSgn; /// number of sigmas to veto the signal peak - double mSigmaSgnErr; /// uncertainty on signal gaussian sigma - double mSigmaSgnDoubleGaus; /// signal 2gaussian sigma - bool mFixedMean; /// switch for fix mean of gaussian - bool mBoundMean; /// switch for bound mean of guassian - bool mBoundReflMean; /// switch for bound mean of guassian for reflection - bool mFixedSigma; /// fix sigma or not - bool mFixedSigmaDoubleGaus; /// fix sigma of 2gaussian or not - bool mBoundSigma; /// set bound sigma or not - bool mFixedDscbTailParams; /// switch for fix double sided Crystal Ball tail parameters - double mSigmaValue; /// value of sigma - double mParamSgn; /// +/- range variation of bound Sigma of gaussian in % - double mFracDoubleGaus; /// initialization for fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar - double mFixedRawYield; /// initialization for raw yield - bool mFixedFracDoubleGaus; /// switch for fixed fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar - double mRatioDoubleGausSigma; /// initialization for ratio between two gaussian sigmas in case of k2GausSigmaRatioPar - bool mFixedRatioDoubleGausSigma; /// switch for fixed ratio between two gaussian sigmas in case of k2GausSigmaRatioPar - double mReflOverSgn; /// reflection/signal - bool mEnableReflections; /// flag use/not use reflections - double mRawYield; /// signal gaussian integral - double mRawYieldErr; /// err on signal gaussian integral - double mRawYieldCounted; /// signal gaussian integral evaluated via bin counting - double mRawYieldCountedErr; /// err on signal gaussian integral evaluated via bin counting - double mBkgYield; /// background - double mBkgYieldErr; /// err on background - double mSignificance; /// significance - double mSignificanceErr; /// err on significance - double mChiSquareOverNdfTotal; /// chi2/ndf of the total fit - double mChiSquareOverNdfBkg; /// chi2/ndf of the background (sidebands) pre-fit - bool mFixReflOverSgn; /// switch for fix refl/signal - double mDscbAlphaLInitialValue; /// double sided Crystal Ball alpha left initial value - double mDscbAlphaLLowLimit; /// double sided Crystal Ball alpha left lower limit - double mDscbAlphaLUpLimit; /// double sided Crystal Ball alpha left upper limit - double mDscbAlphaRInitialValue; /// double sided Crystal Ball alpha right initial value - double mDscbAlphaRLowLimit; /// double sided Crystal Ball alpha right lower limit - double mDscbAlphaRUpLimit; /// double sided Crystal Ball alpha right upper limit - double mDscbNLInitialValue; /// double sided Crystal Ball n left initial value - double mDscbNLLowLimit; /// double sided Crystal Ball n left lower limit - double mDscbNLUpLimit; /// double sided Crystal Ball n left upper limit - double mDscbNRInitialValue; /// double sided Crystal Ball n right initial value - double mDscbNRLowLimit; /// double sided Crystal Ball n right lower limit - double mDscbNRUpLimit; /// double sided Crystal Ball n right upper limit - RooRealVar* mRooMeanSgn; /// mean for gaussian of signal - RooRealVar* mRooSigmaSgn; /// sigma for gaussian of signal - RooRealVar* mRooSecSigmaSgn; /// second sigma for composite gaussian of signal - RooRealVar* mRooFracDoubleGaus; /// fraction of second gaussian for composite gaussian of signal - RooAbsPdf* mSgnPdf; /// signal fit function - RooAbsPdf* mBkgPdf; /// background fit function - RooAbsPdf* mReflPdf; /// reflection fit function - RooRealVar* mRooNSgn; /// total Signal fit function integral - RooRealVar* mRooNBkg; /// total background fit function integral - RooRealVar* mRooNRefl; /// total reflection fit function integral - RooRealVar* mRooDscbAlphaL; /// double sided Crystal Ball alpha left - RooRealVar* mRooDscbAlphaR; /// double sided Crystal Ball alpha right - RooRealVar* mRooDscbNL; /// double sided Crystal Ball n left - RooRealVar* mRooDscbNR; /// double sided Crystal Ball n right - RooAbsPdf* mTotalPdf; /// total fit function - RooPlot* mInvMassFrame; /// frame of mass - RooPlot* mReflFrame; /// reflection frame - RooPlot* mReflOnlyFrame; /// reflection frame plot on reflection only - RooPlot* mResidualFrame; /// residual frame - RooHist* mResidualHist; /// residual histogram - RooPlot* mRatioFrame; /// fit/data ratio frame - RooWorkspace* mWorkspace; /// workspace - double mIntegralBkg; /// integral of background fit function - double mIntegralSgn; /// integral of signal fit function - TH1* mHistoTemplateRefl; /// reflection histogram - bool mDrawBgPrefit; /// draw background after fitting the sidebands - bool mHighlightPeakRegion; /// draw vertical lines showing the peak region (usually +- 3 sigma) - int mRandomSeed; /// seed for random engine for fit's initial parameters randomization - TRandom3* mRandomGen; /// engine for fit's initial parameters randomization - int mFitStatus; /// fit result status, see https://root-forum.cern.ch/t/meaning-of-values-returned-by-roofitresult-status/16355/2 - int mCovQual; /// fit result covariance matrix quality, see https://root.cern.ch/doc/v620/Minuit2Minimizer_8cxx_source.html#l01121 - double mEdm; /// fit quality metrics: Estimated Distance to Minimum - double mMinNll; /// fit quality metrics: minimum negative log-likelihood (NLL) value achieved at the best-fit parameter values - std::vector mSgnCorrelCoeffValues; /// correlation coefficients of mRooNSgn (global and with BG fit parameters) values - std::vector mSgnCorrelCoeffNames; /// correlation coefficients of mRooNSgn (global and with BG fit parameters) names - TH2* mCovCorrMatrix; /// covariance (upper left + diagonal) and correlation (lower right) matrix of free fit parameters + double mMinMass; // lower mass limit + double mMaxMass; // upper mass limit + int mTypeOfBkgPdf; // background fit function + int mTypeOfSgnPdf; // signal fit function + int mTypeOfReflPdf; // reflection fit function + double mMassParticle; // pdg value of particle mass + double mMass; /// signal gaussian mean value + double mMassLowLimit; /// lower limit of the allowed mass range + double mMassUpLimit; /// upper limit of the allowed mass range + double mMassReflLowLimit; /// lower limit of the allowed mass range for reflection + double mMassReflUpLimit; /// upper limit of the allowed mass range for reflection + double mSecMass; /// Second peak mean value + double mSigmaSgn; /// signal gaussian sigma + double mSecSigma; /// Second peak gaussian sigma + double mNSigmaForSidebands; /// number of sigmas to veto the signal peak + double mNSigmaForSgn; /// number of sigmas to veto the signal peak + double mSigmaSgnErr; /// uncertainty on signal gaussian sigma + double mSigmaSgnDoubleGaus; /// signal 2gaussian sigma + bool mFixedMean; /// switch for fix mean of gaussian + bool mBoundMean; /// switch for bound mean of guassian + bool mBoundReflMean; /// switch for bound mean of guassian for reflection + bool mFixedSigma; /// fix sigma or not + bool mFixedSigmaDoubleGaus; /// fix sigma of 2gaussian or not + bool mBoundSigma; /// set bound sigma or not + bool mFixedDscbTailParams; /// switch for fix double sided Crystal Ball tail parameters + double mSigmaValue; /// value of sigma + double mParamSgn; /// +/- range variation of bound Sigma of gaussian in % + double mFracDoubleGaus; /// initialization for fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar + double mFixedRawYield; /// initialization for raw yield + bool mFixedFracDoubleGaus; /// switch for fixed fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar + double mRatioDoubleGausSigma; /// initialization for ratio between two gaussian sigmas in case of k2GausSigmaRatioPar + bool mFixedRatioDoubleGausSigma; /// switch for fixed ratio between two gaussian sigmas in case of k2GausSigmaRatioPar + double mReflOverSgn; /// reflection/signal + bool mEnableReflections; /// flag use/not use reflections + double mRawYield; /// signal gaussian integral + double mRawYieldErr; /// err on signal gaussian integral + double mRawYieldCounted; /// signal gaussian integral evaluated via bin counting + double mRawYieldCountedErr; /// err on signal gaussian integral evaluated via bin counting + double mBkgYield; /// background + double mBkgYieldErr; /// err on background + double mSignificance; /// significance + double mSignificanceErr; /// err on significance + double mChiSquareOverNdfTotal; /// chi2/ndf of the total fit + double mChiSquareOverNdfBkg; /// chi2/ndf of the background (sidebands) pre-fit + bool mFixReflOverSgn; /// switch for fix refl/signal + double mDscbAlphaLInitialValue; /// double sided Crystal Ball alpha left initial value + double mDscbAlphaLLowLimit; /// double sided Crystal Ball alpha left lower limit + double mDscbAlphaLUpLimit; /// double sided Crystal Ball alpha left upper limit + double mDscbAlphaRInitialValue; /// double sided Crystal Ball alpha right initial value + double mDscbAlphaRLowLimit; /// double sided Crystal Ball alpha right lower limit + double mDscbAlphaRUpLimit; /// double sided Crystal Ball alpha right upper limit + double mDscbNLInitialValue; /// double sided Crystal Ball n left initial value + double mDscbNLLowLimit; /// double sided Crystal Ball n left lower limit + double mDscbNLUpLimit; /// double sided Crystal Ball n left upper limit + double mDscbNRInitialValue; /// double sided Crystal Ball n right initial value + double mDscbNRLowLimit; /// double sided Crystal Ball n right lower limit + double mDscbNRUpLimit; /// double sided Crystal Ball n right upper limit + RooRealVar* mRooMeanSgn; /// mean for gaussian of signal + RooRealVar* mRooSigmaSgn; /// sigma for gaussian of signal + RooRealVar* mRooSecSigmaSgn; /// second sigma for composite gaussian of signal + RooRealVar* mRooFracDoubleGaus; /// fraction of second gaussian for composite gaussian of signal + RooAbsPdf* mSgnPdf; /// signal fit function + RooAbsPdf* mBkgPdf; /// background fit function + RooAbsPdf* mReflPdf; /// reflection fit function + RooRealVar* mRooNSgn; /// total Signal fit function integral + RooRealVar* mRooNBkg; /// total background fit function integral + RooRealVar* mRooNRefl; /// total reflection fit function integral + RooRealVar* mRooDscbAlphaL; /// double sided Crystal Ball alpha left + RooRealVar* mRooDscbAlphaR; /// double sided Crystal Ball alpha right + RooRealVar* mRooDscbNL; /// double sided Crystal Ball n left + RooRealVar* mRooDscbNR; /// double sided Crystal Ball n right + RooAbsPdf* mTotalPdf; /// total fit function + RooPlot* mInvMassFrame; /// frame of mass + RooPlot* mReflFrame; /// reflection frame + RooPlot* mReflOnlyFrame; /// reflection frame plot on reflection only + RooPlot* mResidualFrame; /// residual frame + RooHist* mResidualHist; /// residual histogram + RooPlot* mRatioFrame; /// fit/data ratio frame + RooWorkspace* mWorkspace; /// workspace + double mIntegralBkg; /// integral of background fit function + double mIntegralSgn; /// integral of signal fit function + TH1* mHistoTemplateRefl; /// reflection histogram + bool mDrawBgPrefit; /// draw background after fitting the sidebands + bool mHighlightPeakRegion; /// draw vertical lines showing the peak region (usually +- 3 sigma) + int mRandomSeed; /// seed for random engine for fit's initial parameters randomization + TRandom3* mRandomGen; /// engine for fit's initial parameters randomization + int mFitStatus; /// fit result status, see https://root-forum.cern.ch/t/meaning-of-values-returned-by-roofitresult-status/16355/2 + int mCovQual; /// fit result covariance matrix quality, see https://root.cern.ch/doc/v620/Minuit2Minimizer_8cxx_source.html#l01121 + double mEdm; /// fit quality metrics: Estimated Distance to Minimum + double mMinNll; /// fit quality metrics: minimum negative log-likelihood (NLL) value achieved at the best-fit parameter values + TH2* mCovCorrMatrix; /// covariance (upper left + diagonal) and correlation (lower right) matrix of free fit parameters ClassDefOverride(HFInvMassFitter, 1); }; diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index a70f5caa3fc..135fff1c72b 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -321,7 +321,6 @@ void runMassFitter(const std::string& configFileName) std::vector hMassSgn(nHistograms); std::vector hMassRefl(nHistograms); std::vector hMass(nHistograms); - std::vector hSgnCorr(nHistograms); std::vector hCovCorr(nHistograms); for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { @@ -678,16 +677,6 @@ void runMassFitter(const std::string& configFileName) hFitResult->SetBinContent(FitResultEdm, iSliceVar + 1, massFitter->getEDM()); hFitResult->SetBinContent(FitResultMinNll, iSliceVar + 1, massFitter->getMinNll()); - const std::string hSgnCorrName = "hSgnCorr" + std::to_string(iSliceVar + 1); - const auto& sgnCorrValues = massFitter->getSgnCorrelCoeffValues(); - const auto& sgnCorrNames = massFitter->getSgnCorrelCoeffNames(); - const int nSgnCorrValues = sgnCorrValues.size(); - hSgnCorr[iSliceVar] = new TH1D(hSgnCorrName.c_str(), hSgnCorrName.c_str(), nSgnCorrValues, 0, nSgnCorrValues); - for (int iSgnCorrValue = 0; iSgnCorrValue < nSgnCorrValues; ++iSgnCorrValue) { - hSgnCorr[iSliceVar]->SetBinContent(iSgnCorrValue + 1, sgnCorrValues.at(iSgnCorrValue)); - hSgnCorr[iSliceVar]->GetXaxis()->SetBinLabel(iSgnCorrValue + 1, sgnCorrNames.at(iSgnCorrValue).c_str()); - } - hCovCorr[iSliceVar] = massFitter->getCovCorrMatrix(); } @@ -711,8 +700,7 @@ void runMassFitter(const std::string& configFileName) } outputFile.cd("MassHistograms"); hMass[iSliceVar]->Write(); - outputFile.cd("SgnCorrHistograms"); - hSgnCorr[iSliceVar]->Write(); + outputFile.cd("CovCorrMatrices"); hCovCorr[iSliceVar]->Write(Form("hCovCorrMatrix%d", iSliceVar + 1)); } outputFile.cd(); From 59b028df462b301994c6d9c705269ca13d63e838 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Fri, 31 Jul 2026 14:05:47 +0200 Subject: [PATCH 26/34] add mNSgn global correlation coefficient to hFitResult (x2); bugfix directory name --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 2 ++ PWGHF/D2H/Macros/HFInvMassFitter.h | 2 ++ PWGHF/D2H/Macros/runMassFitter.C | 6 ++++-- 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index b92919a8714..775ecd7e01e 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -160,6 +160,7 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mCovQual(-999), mEdm(-999.), mMinNll(-999.), + mSgnGlobalCorrelCoeff(-999.), mCovCorrMatrix(nullptr) { // standard constructor @@ -342,6 +343,7 @@ void HFInvMassFitter::doFit() mCovQual = fitResult->covQual(); mEdm = fitResult->edm(); mMinNll = fitResult->minNll(); + mSgnGlobalCorrelCoeff = fitResult->globalCorr("mNSgn"); mCovCorrMatrix = fillCovCorrMatrix(fitResult); diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 6fe6538f6cc..6958c68ba36 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -169,6 +169,7 @@ class HFInvMassFitter : public TNamed [[nodiscard]] int getCovQual() const { return mCovQual; } [[nodiscard]] double getEDM() const { return mEdm; } [[nodiscard]] double getMinNll() const { return mMinNll; } + [[nodiscard]] double getSgnGlobalCorrelCoeff() const { return mSgnGlobalCorrelCoeff; } [[nodiscard]] TH2* getCovCorrMatrix() const { return mCovCorrMatrix; } void calculateSignal(double& signal, double& signalErr) const; void countSignal(double& signal, double& signalErr) const; @@ -284,6 +285,7 @@ class HFInvMassFitter : public TNamed int mCovQual; /// fit result covariance matrix quality, see https://root.cern.ch/doc/v620/Minuit2Minimizer_8cxx_source.html#l01121 double mEdm; /// fit quality metrics: Estimated Distance to Minimum double mMinNll; /// fit quality metrics: minimum negative log-likelihood (NLL) value achieved at the best-fit parameter values + double mSgnGlobalCorrelCoeff; /// global correlation coefficient of mRooNSgn with other fit parameters TH2* mCovCorrMatrix; /// covariance (upper left + diagonal) and correlation (lower right) matrix of free fit parameters ClassDefOverride(HFInvMassFitter, 1); diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 135fff1c72b..f231064a6e4 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -391,12 +391,13 @@ void runMassFitter(const std::string& configFileName) FitResultCovQual, FitResultEdm, FitResultMinNll, + FitResultNSgnGCC, NFitResultsToSave }; auto* hFitConfig = new TH2F("hFitConfig", "Fit Configurations", NConfigsToSave - 1, 0, NConfigsToSave - 1, nHistograms, sliceVarLimits.data()); const char* hFitConfigXLabel[NConfigsToSave - 1] = {"mass min", "mass max", "rebin num", "fix sigma", "bkg func", "sgn func", "rnd seed"}; auto* hFitResult = new TH2F("hFitResult", "Fit Result", NFitResultsToSave - 1, 0, NFitResultsToSave - 1, nHistograms, sliceVarLimits.data()); - const char* hFitResultXLabel[NFitResultsToSave - 1] = {"status", "cov qual", "edm", "minNLL"}; + const char* hFitResultXLabel[NFitResultsToSave - 1] = {"status", "cov qual", "edm", "minNLL", "N Sig GCC"}; for (int i = 0; i < NConfigsToSave - 1; i++) { hFitConfig->GetXaxis()->SetBinLabel(i + 1, hFitConfigXLabel[i]); } @@ -676,6 +677,7 @@ void runMassFitter(const std::string& configFileName) hFitResult->SetBinContent(FitResultCovQual, iSliceVar + 1, massFitter->getCovQual()); hFitResult->SetBinContent(FitResultEdm, iSliceVar + 1, massFitter->getEDM()); hFitResult->SetBinContent(FitResultMinNll, iSliceVar + 1, massFitter->getMinNll()); + hFitResult->SetBinContent(FitResultNSgnGCC, iSliceVar + 1, massFitter->getSgnGlobalCorrelCoeff()); hCovCorr[iSliceVar] = massFitter->getCovCorrMatrix(); } @@ -696,7 +698,7 @@ void runMassFitter(const std::string& configFileName) for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { if (iSliceVar == 0) { outputFile.mkdir("MassHistograms"); - outputFile.mkdir("SgnCorrHistograms"); + outputFile.mkdir("CovCorrMatrices"); } outputFile.cd("MassHistograms"); hMass[iSliceVar]->Write(); From 425a34e67eb03cbcecb03a7d40fa96cc9f77f6b4 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Tue, 4 Aug 2026 23:24:24 +0200 Subject: [PATCH 27/34] handle exception properly; do not write hCovCorr in case of esception --- PWGHF/D2H/Macros/runMassFitter.C | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index f231064a6e4..774c61aff96 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -574,8 +574,8 @@ void runMassFitter(const std::string& configFileName) try { massFitter->doFit(); - } catch (...) { - std::cout << "exception caught while doing fit\n"; + } catch (const std::exception& e) { + printf("Warinig! Exception \"%s\" caught while doing fit of the histogram no. %d. The fitting process will be continued without it.\n", e.what(), iSliceVar); } auto drawOnCanvas = [&](std::vector& canvas, std::function drawer) { @@ -703,7 +703,9 @@ void runMassFitter(const std::string& configFileName) outputFile.cd("MassHistograms"); hMass[iSliceVar]->Write(); outputFile.cd("CovCorrMatrices"); - hCovCorr[iSliceVar]->Write(Form("hCovCorrMatrix%d", iSliceVar + 1)); + if (hCovCorr[iSliceVar] != nullptr) { + hCovCorr[iSliceVar]->Write(Form("hCovCorrMatrix%d", iSliceVar + 1)); + } } outputFile.cd(); hRawYieldsSignal->Write(); From 18417db01ece0f1ecc171085e14ee22800e53bbf Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Tue, 4 Aug 2026 23:28:20 +0200 Subject: [PATCH 28/34] clean-up debug std::cout's --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 48 +--------------------------- 1 file changed, 1 insertion(+), 47 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 775ecd7e01e..7a3104a7249 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -250,18 +250,13 @@ void HFInvMassFitter::doFit() } else { // data const double integralSidebands = integrateHistoInvMassOverWorkspaceRanges({"SBL", "SBR"}); const ParameterRanges rooNBkgParamRanges{0.5 * integralSidebands, 1.5 * integralSidebands, integralSidebands}; - std::cout << "integralSidebands = " << integralSidebands << "\n"; - std::cout << "rooNBkgParamRanges:\n"; mRooNBkg = new RooRealVar("mRooNBkg", "number of background", randomizeInitialParameter(rooNBkgParamRanges), rooNBkgParamRanges.lower, rooNBkgParamRanges.upper); // background yield mBkgPdf = new RooAddPdf("mBkgPdf", "background fit function", RooArgList(*bkgPdf), RooArgList(*mRooNBkg)); std::string sbRanges{"SBL,SBR"}; if (mTypeOfSgnPdf == GausSec) { // two peak fit sbRanges.append(",SEC"); } - std::cout << "Start prefit of BG sidebands\n"; mBkgPdf->chi2FitTo(sbHistogram, DataError(RooAbsData::SumW2), Save()); - std::cout << "Finish prefit of BG sidebands\n"; - std::cout << "mRooNBkg->getVal() = " << mRooNBkg->getVal() << "\n"; // define the frame to evaluate background sidebands chi2 (bg pdf needs to be plotted within sideband ranges) RooPlot* frameTemporary = mass->frame(Title(Form("%s_temp", mHistoInvMass->GetTitle()))); @@ -282,7 +277,6 @@ void HFInvMassFitter::doFit() checkForSignal(estimatedSignal); // SIG's absolute integral in "bkg" range calculateBackground(mBkgYield, mBkgYieldErr); // BG's absolute integral in "bkg" range - std::cout << "estimatedSignal = " << estimatedSignal << "\n"; const ParameterRanges rooNSgnParamRanges{0.1 * estimatedSignal, 10 * estimatedSignal, estimatedSignal}; mRooNSgn = new RooRealVar("mNSgn", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // estimated signal yield if (mFixedRawYield > 0) { @@ -315,7 +309,7 @@ void HFInvMassFitter::doFit() } else { mTotalPdf = new RooAddPdf("mTotalPdf", "background + signal pdf", RooArgList(*bkgPdf, *sgnPdf), RooArgList(*mRooNBkg, *mRooNSgn)); } - std::cout << "Start total fit\n"; + RooFitResult* fitResult{nullptr}; if (strcmp(mFitOption.c_str(), "Chi2") == 0) { fitResult = mTotalPdf->chi2FitTo(dataHistogram, Save()); @@ -323,38 +317,13 @@ void HFInvMassFitter::doFit() fitResult = mTotalPdf->fitTo(dataHistogram, Save()); } - std::cout << "Finish total fit\n"; - - std::cout << "History:\t"; - for (unsigned i = 0; i < fitResult->numStatusHistory(); ++i) { - std::cout << fitResult->statusLabelHistory(i) << " : " << fitResult->statusCodeHistory(i) << "\t"; - } - std::cout << "\n"; - - std::cout << "Status = " << fitResult->status() << "\n"; - std::cout << "CovQual = " << fitResult->covQual() << "\n"; - std::cout << "EDM = " << fitResult->edm() << "\n"; - std::cout << "minNLL = " << fitResult->minNll() << "\n"; - std::cout << "Global correlation of mNSgn = " << fitResult->globalCorr("mNSgn") << "\n"; - fitResult->Print("v"); - fitResult->covarianceMatrix().Print(); - fitResult->correlationMatrix().Print(); mFitStatus = fitResult->status(); mCovQual = fitResult->covQual(); mEdm = fitResult->edm(); mMinNll = fitResult->minNll(); mSgnGlobalCorrelCoeff = fitResult->globalCorr("mNSgn"); - mCovCorrMatrix = fillCovCorrMatrix(fitResult); - std::cout << "mRooNBkg->getVal() = " << mRooNBkg->getVal() << "\n"; - std::cout << "mRooNSgn->getVal() = " << mRooNSgn->getVal() << "\n"; - - std::cout << "Value = " << mRooNSgn->getVal() << "\n"; - std::cout << "HESSE error = " << mRooNSgn->getError() << "\n"; - std::cout << "MINOS lower = " << mRooNSgn->getAsymErrorLo() << "\n"; - std::cout << "MINOS upper = " << mRooNSgn->getAsymErrorHi() << "\n"; - plotBkg(mTotalPdf); mTotalPdf->plotOn(mInvMassFrame, Name("Tot_c"), LineColor(kBlue)); if (mHistoTemplateRefl != nullptr) { @@ -606,7 +575,6 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const // draw fit output void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& plotLabels, bool writeParInfo) { - std::cout << "drawFit() start\n"; gStyle->SetOptStat(0); gStyle->SetCanvasColor(0); gStyle->SetFrameFillColor(0); @@ -667,13 +635,11 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& if (mHistoTemplateRefl) { mReflFrame->Draw("same"); } - std::cout << "drawFit() finish\n"; } // draw residual distribution on canvas void HFInvMassFitter::drawResidual(TVirtualPad* pad) { - std::cout << "drawResidual() start\n"; if (mResidualFrame == nullptr) { printf("Warning HFInvMassFitter::drawResidual(): mResidualFrame == nullptr and will not be drawn\n"); return; @@ -697,13 +663,11 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad) mResidualFrame->GetYaxis()->SetTitleOffset(1.8); mResidualFrame->Draw(); highlightPeakRegion(mResidualFrame); - std::cout << "drawResidual() finish\n"; } // draw ratio on canvas void HFInvMassFitter::drawRatio(TVirtualPad* pad) { - std::cout << "drawRatio() start\n"; if (mRatioFrame == nullptr) { printf("Warning HFInvMassFitter::drawRatio(): mRatioFrame == nullptr and will not be drawn\n"); return; @@ -723,7 +687,6 @@ void HFInvMassFitter::drawRatio(TVirtualPad* pad) gPad->SetRightMargin(0.02); mRatioFrame->Draw(); highlightPeakRegion(mRatioFrame); - std::cout << "drawRatio() finish\n"; } // draw peak region with vertical lines @@ -797,7 +760,6 @@ void HFInvMassFitter::calculateSignal(double& signal, double& errSignal) const // calculate background yield void HFInvMassFitter::calculateBackground(double& bkg, double& errBkg) const { - std::cout << "calculateBackground()\n"; if (mTypeOfBkgPdf == NoBkg) { bkg = 0.; errBkg = 0.; @@ -825,10 +787,8 @@ void HFInvMassFitter::calculateSignificance(double& significance, double& errSig // estimate Signal void HFInvMassFitter::checkForSignal(double& estimatedSignal) { - std::cout << "checkForSignal()\n"; const double integralHisto = integrateHistoInvMassOverWorkspaceRanges({"full"}); const double bkg = mRooNBkg->getVal(); - std::cout << "integralHisto = " << integralHisto << ", bkg = " << bkg << "\n"; estimatedSignal = integralHisto - bkg; } @@ -1197,15 +1157,11 @@ double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& paramet } } while (result < parameterRanges.lower || result > parameterRanges.upper); - std::cout << "randomizeInitialParameter():\nfrom " << parameterRanges.lower << "\nto " << parameterRanges.upper << "\ninitial " << parameterRanges.initial << "\nsigma " << sigma << "\n"; - std::cout << "randomized to " << result << "\n"; - return result; } double HFInvMassFitter::integrateHistoInvMassOverWorkspaceRanges(const std::vector& ranges) const { - std::cout << "integrateHistoInvMassOverWorkspaceRanges()\n"; double sumEntries{0.}; double sumLengths{0.}; for (const auto& range : ranges) { @@ -1216,8 +1172,6 @@ double HFInvMassFitter::integrateHistoInvMassOverWorkspaceRanges(const std::vect sumEntries += mHistoInvMass->GetBinContent(binLo) * (mHistoInvMass->GetBinLowEdge(binLo + 1) - lo) / mHistoInvMass->GetBinWidth(binLo); sumEntries += mHistoInvMass->GetBinContent(binHi) * (hi - mHistoInvMass->GetBinLowEdge(binHi)) / mHistoInvMass->GetBinWidth(binHi); sumLengths += (hi - lo); - std::cout << "lo = " << lo << ", hi = " << hi << "\n"; - std::cout << "sumEntries = " << sumEntries << "\n"; } const auto [fullLo, fullHi] = mWorkspace->var("mass")->getRange("full"); const double fullLength = fullHi - fullLo; From 6e2e56b1468805ab70024dce33cd1f48cb24eb0b Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 5 Aug 2026 00:01:27 +0200 Subject: [PATCH 29/34] address (part of) clang-tidy errors --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 15 +++++++-------- PWGHF/D2H/Macros/HFInvMassFitter.h | 26 +++++++++++++------------- 2 files changed, 20 insertions(+), 21 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 7a3104a7249..8c1bb01d5c6 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -68,7 +68,7 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, double maxValue, int fitTypeBkg, int fitTypeSgn, - int randomSeed) : mHistoInvMass(nullptr), + int randomSeed) : mHistoInvMass(histoToFit), mFitOption("L,E"), mMinMass(minValue), mMaxMass(maxValue), @@ -164,7 +164,6 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mCovCorrMatrix(nullptr) { // standard constructor - mHistoInvMass = histoToFit; mHistoInvMass->SetDirectory(nullptr); if (mRandomSeed >= 0) { mRandomGen = new TRandom3(); @@ -914,7 +913,7 @@ void HFInvMassFitter::calculateFitToDataRatio() const return; } - RooHist* ratioHist = new RooHist(); + auto* ratioHist = new RooHist(); for (int i = 0; i < dataHist->GetN(); ++i) { double x{}, dataY{}, dataErr{}; @@ -1145,15 +1144,15 @@ double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& paramet } const auto sigma = parameterRanges.sigma < 0 ? (parameterRanges.upper - parameterRanges.lower) / DefaultSigmaFraction : parameterRanges.sigma; - double result; + double result{}; int nIter{0}; do { result = mRandomGen->Gaus(parameterRanges.initial, sigma); ++nIter; if (nIter > MaximalNumberOfIterations) { - char errorMessage[200]; - std::snprintf(errorMessage, sizeof(errorMessage), "randomizeInitialParameter() - long while loop with lower = %f upper = %f initial = %f sigma = %f\n", parameterRanges.lower, parameterRanges.upper, parameterRanges.initial, sigma); - throw std::runtime_error(errorMessage); + std::array errorMessage{}; + std::snprintf(errorMessage.data(), errorMessage.size(), "randomizeInitialParameter() - long while loop with lower = %f upper = %f initial = %f sigma = %f\n", parameterRanges.lower, parameterRanges.upper, parameterRanges.initial, sigma); + throw std::runtime_error(errorMessage.data()); } } while (result < parameterRanges.lower || result > parameterRanges.upper); @@ -1207,7 +1206,7 @@ void HFInvMassFitter::cutRangesFromHisto(TH1* histo, const std::vectorfloatParsFinal(); - const int nPars = pars.size(); + const int nPars = static_cast(pars.size()); TH2* hMatrix = new TH2D("covCorrMatrix", "covariance (upper left + diagonal) and correlation (lower right) matrix", nPars, 0, nPars, nPars, 0, nPars); for (int iPar = 0; iPar < nPars; ++iPar) { diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 6958c68ba36..0a3afb0b983 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -56,28 +56,28 @@ class HFInvMassFitter : public TNamed public: enum TypeOfBkgPdf { Expo = 0, - Poly1 = 1, - Poly2 = 2, - Pow = 3, - PowExpo = 4, - Poly3 = 5, - NoBkg = 6, + Poly1, // 1 + Poly2, // 2 + Pow, // 3 + PowExpo, // 4 + Poly3, // 5 + NoBkg, // 6 NTypesOfBkgPdf }; std::array namesOfBkgPdf{"bkgFuncExpo", "bkgFuncPoly1", "bkgFuncPoly2", "bkgFuncPow", "bkgFuncPowExpo", "bkgFuncPoly3"}; enum TypeOfSgnPdf { SingleGaus = 0, - DoubleGaus = 1, - DoubleGausSigmaRatioPar = 2, - GausSec = 3, - DoubleSidedCrystalBall = 4, + DoubleGaus, // 1 + DoubleGausSigmaRatioPar, // 2 + GausSec, // 3 + DoubleSidedCrystalBall, // 4 NTypesOfSgnPdf }; enum TypeOfReflPdf { SingleGausRefl = 0, - DoubleGausRefl = 1, - Poly3Refl = 2, - Poly6Refl = 3, + DoubleGausRefl, // 1 + Poly3Refl, // 2 + Poly6Refl, // 3 NTypesOfReflPdf }; std::array namesOfReflPdf{"reflFuncGaus", "reflFuncDoubleGaus", "reflFuncPoly3", "reflFuncPoly6"}; From f8b5c34d8562b7e3a51c9ede259b7af16fd1d87d Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 5 Aug 2026 00:04:02 +0200 Subject: [PATCH 30/34] do not delete not owned mHistoInvMass --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 8c1bb01d5c6..cd22c8aff20 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -175,7 +175,6 @@ HFInvMassFitter::~HFInvMassFitter() { /// destructor - delete mHistoInvMass; delete mHistoTemplateRefl; delete mRooMeanSgn; delete mRooSigmaSgn; From 5c0fa661eaf456487601ec19d4367e73144411d9 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 5 Aug 2026 12:07:28 +0200 Subject: [PATCH 31/34] fix clang-tidy errors --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 11 ++-- PWGHF/D2H/Macros/HFInvMassFitter.h | 26 ++++---- PWGHF/D2H/Macros/runMassFitter.C | 93 ++++++++++++++-------------- 3 files changed, 63 insertions(+), 67 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index cd22c8aff20..3bcf783ba42 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -718,7 +718,7 @@ void HFInvMassFitter::drawReflection(TVirtualPad* pad) } // calculate signal yield via bin counting -void HFInvMassFitter::countSignal(double& signal, double& signalErr) const +void HFInvMassFitter::countSignal(double& signal, double& errSignal) const { const auto& histoForCounting = mTypeOfBkgPdf == NoBkg ? mInvMassFrame->getHist("data_c") : mResidualHist; @@ -745,7 +745,7 @@ void HFInvMassFitter::countSignal(double& signal, double& signalErr) const sumErrorsSquare += square(histoForCounting->GetErrorY(binForMaxSgn - 1) * binForMaxSgnFraction); signal = sumValues; - signalErr = std::sqrt(sumErrorsSquare); + errSignal = std::sqrt(sumErrorsSquare); } // calculate signal yield @@ -796,11 +796,10 @@ std::pair HFInvMassFitter::getRangesOfSignal() const { if (mTypeOfSgnPdf == DoubleSidedCrystalBall) { return std::make_pair(mMinMass, mMaxMass); - } else { - const double mean = mRooMeanSgn->getVal(); - const double sigma = mRooSecSigmaSgn->getVal(); - return std::make_pair(mean - mNSigmaForSgn * sigma, mean + mNSigmaForSgn * sigma); } + const double mean = mRooMeanSgn->getVal(); + const double sigma = mRooSecSigmaSgn->getVal(); + return std::make_pair(mean - mNSigmaForSgn * sigma, mean + mNSigmaForSgn * sigma); } // Create Background Fit Function diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 0a3afb0b983..be5c95be829 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -88,9 +88,9 @@ class HFInvMassFitter : public TNamed void setUseLikelihoodFit() { mFitOption = "L,E"; } void setUseChi2Fit() { mFitOption = "Chi2"; } void setFitOption(const std::string& opt) { mFitOption = opt; } - RooAbsPdf* createBackgroundFitFunction(RooWorkspace* w1) const; - RooAbsPdf* createSignalFitFunction(RooWorkspace* w1); - RooAbsPdf* createReflectionFitFunction(RooWorkspace* w1) const; + RooAbsPdf* createBackgroundFitFunction(RooWorkspace* workspace) const; + RooAbsPdf* createSignalFitFunction(RooWorkspace* workspace); + RooAbsPdf* createReflectionFitFunction(RooWorkspace* workspace) const; void setFitRange(double minValue, double maxValue); void setFitFunctions(int fitTypeBkg, int fitTypeSgn); @@ -129,8 +129,8 @@ class HFInvMassFitter : public TNamed void setDscbNRInitialValue(double value) { mDscbNRInitialValue = value; } void setDscbNRLowLimit(double value) { mDscbNRLowLimit = value; } void setDscbNRUpLimit(double value) { mDscbNRUpLimit = value; } - void plotBkg(RooAbsPdf* mFunc, Color_t color = kRed); - void plotRefl(RooAbsPdf* mFunc); + void plotBkg(RooAbsPdf* pdf, Color_t color = kRed); + void plotRefl(RooAbsPdf* pdf); void setReflFuncFixed(); void doFit(); void setInitialReflOverSgn(double reflOverSgn) { mReflOverSgn = reflOverSgn; } @@ -171,16 +171,16 @@ class HFInvMassFitter : public TNamed [[nodiscard]] double getMinNll() const { return mMinNll; } [[nodiscard]] double getSgnGlobalCorrelCoeff() const { return mSgnGlobalCorrelCoeff; } [[nodiscard]] TH2* getCovCorrMatrix() const { return mCovCorrMatrix; } - void calculateSignal(double& signal, double& signalErr) const; - void countSignal(double& signal, double& signalErr) const; - void calculateBackground(double& bkg, double& bkgErr) const; - void calculateSignificance(double& significance, double& significanceErr) const; + void calculateSignal(double& signal, double& errSignal) const; + void countSignal(double& signal, double& errSignal) const; + void calculateBackground(double& bkg, double& errBkg) const; + void calculateSignificance(double& significance, double& errSignificance) const; void checkForSignal(double& estimatedSignal); void calculateFitToDataRatio() const; - void drawFit(TVirtualPad* c, const std::vector& plotLabels, bool writeParInfo = true); - void drawResidual(TVirtualPad* c); - void drawRatio(TVirtualPad* c); - void drawReflection(TVirtualPad* c); + void drawFit(TVirtualPad* pad, const std::vector& plotLabels, bool writeParInfo = true); + void drawResidual(TVirtualPad* pad); + void drawRatio(TVirtualPad* pad); + void drawReflection(TVirtualPad* pad); private: HFInvMassFitter(const HFInvMassFitter& source); diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 774c61aff96..eec8983d9c8 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -52,8 +52,6 @@ using namespace rapidjson; -constexpr int UndefValueInt{-999}; - void runMassFitter(const std::string& configFileName = "config_massfitter.json"); TFile* openFileWithNullptrCheck(const std::string& fileName, const std::string& option = "read"); @@ -88,8 +86,8 @@ void runMassFitter(const std::string& configFileName) } Document config; - char readBuffer[65536]; - FileReadStream is(configFile, readBuffer, sizeof(readBuffer)); + std::array readBuffer{}; + FileReadStream is(configFile, readBuffer.data(), readBuffer.size()); config.ParseStream(is); fclose(configFile); @@ -97,12 +95,12 @@ void runMassFitter(const std::string& configFileName) throw std::runtime_error("ERROR: Parsing the configuration json file failed. Check the config for correct formatting"); } - bool const isMc = readJsonField(config, "IsMC"); - bool const writeSignalPar = readJsonField(config, "WriteSignalPar", true); - std::string const particleName = readJsonField(config, "Particle"); - std::string const collisionSystem = readJsonField(config, "CollisionSystem", ""); - std::string const inputFileName = readJsonField(config, "InFileName"); - std::string const reflFileName = readJsonField(config, "ReflFileName", ""); + auto const isMc = readJsonField(config, "IsMC"); + auto const writeSignalPar = readJsonField(config, "WriteSignalPar", true); + auto const particleName = readJsonField(config, "Particle"); + auto const collisionSystem = readJsonField(config, "CollisionSystem", ""); + auto const inputFileName = readJsonField(config, "InFileName"); + auto const reflFileName = readJsonField(config, "ReflFileName", ""); TString outputFileName = readJsonField(config, "OutFileName", "mInvFit.root"); std::vector inputHistoName; @@ -149,23 +147,23 @@ void runMassFitter(const std::string& configFileName) readJsonVector(fdSecPeakHistoName, config, "FDSecPeakHistoName"); readJsonVector(signalSecPeakHistoName, config, "SignalSecPeakHistoName"); - const bool fixMean = readJsonField(config, "FixMean", false); - const std::string meanFile = readJsonField(config, "MeanFile", ""); + const auto fixMean = readJsonField(config, "FixMean", false); + const auto meanFile = readJsonField(config, "MeanFile", ""); readJsonVectorFlexible(fixMeanManual, config, nHistograms, "FixMeanManual"); - const bool fixSigma = readJsonField(config, "FixSigma", false); - const std::string sigmaFile = readJsonField(config, "SigmaFile", ""); + const auto fixSigma = readJsonField(config, "FixSigma", false); + const auto sigmaFile = readJsonField(config, "SigmaFile", ""); readJsonVectorFlexible(fixSigmaManual, config, nHistograms, "FixSigmaManual"); - const bool fixSecondSigma = readJsonField(config, "FixSecondSigma", false); - const std::string secondSigmaFile = readJsonField(config, "SecondSigmaFile", ""); + const auto fixSecondSigma = readJsonField(config, "FixSecondSigma", false); + const auto secondSigmaFile = readJsonField(config, "SecondSigmaFile", ""); readJsonVectorFlexible(fixSecondSigmaManual, config, nHistograms, "FixSecondSigmaManual"); - const bool fixFracDoubleGaus = readJsonField(config, "FixFracDoubleGaus", false); - const std::string fracDoubleGausFile = readJsonField(config, "FracDoubleGausFile", ""); + const auto fixFracDoubleGaus = readJsonField(config, "FixFracDoubleGaus", false); + const auto fracDoubleGausFile = readJsonField(config, "FracDoubleGausFile", ""); readJsonVectorFlexible(fixFracDoubleGausManual, config, nHistograms, "FixFracDoubleGausManual"); - const bool fixDscbTailParams = readJsonField(config, "FixDscbTailParams", false); + const auto fixDscbTailParams = readJsonField(config, "FixDscbTailParams", false); const TString sliceVarName = readJsonField(config, "SliceVarName"); const TString sliceVarUnit = readJsonField(config, "SliceVarUnit"); @@ -178,18 +176,18 @@ void runMassFitter(const std::string& configFileName) readJsonVectorFlexible(nRebin, config, nHistograms, "Rebin", true); - bool const includeSecPeak = readJsonField(config, "InclSecPeak", false); - bool const useLikelihood = readJsonField(config, "UseLikelihood"); + auto const includeSecPeak = readJsonField(config, "InclSecPeak", false); + auto const useLikelihood = readJsonField(config, "UseLikelihood"); readJsonVectorFlexible(bkgFunc, config, nHistograms, "BkgFunc", true); readJsonVectorFlexible(sgnFunc, config, nHistograms, "SgnFunc", true); - const bool enableRefl = readJsonField(config, "EnableRefl", false); - const bool drawBgPrefit = readJsonField(config, "DrawBgPrefit", true); - const bool highlightPeakRegion = readJsonField(config, "HighlightPeakRegion", true); - const int randomSeed = readJsonField(config, "RandomSeed", -1); - const double nSigmaForSideband = readJsonField(config, "NSigmaForSideband", 3.); - const double nSigmaForSignal = readJsonField(config, "NSigmaForSignal", 3.); + const auto enableRefl = readJsonField(config, "EnableRefl", false); + const auto drawBgPrefit = readJsonField(config, "DrawBgPrefit", true); + const auto highlightPeakRegion = readJsonField(config, "HighlightPeakRegion", true); + const auto randomSeed = readJsonField(config, "RandomSeed", -1); + const auto nSigmaForSideband = readJsonField(config, "NSigmaForSideband", 3.); + const auto nSigmaForSignal = readJsonField(config, "NSigmaForSignal", 3.); readJsonVector(dscbAlphaLInitial, config, "DscbAlphaLInitial"); readJsonVector(dscbAlphaLLower, config, "DscbAlphaLLower"); @@ -305,13 +303,13 @@ void runMassFitter(const std::string& configFileName) {"LcToPK0s", {"pK^{0}_{s}", "Lambda_c+", "#Lambda_{c}^{+}", "pK^{0}_{s}"}}, {"Dstar", {"D^{0}pi^{+}", "D*+", "D^{*+}", "D^{0}#pi^{+}"}}, {"XicToXiPiPi", {"#Xi#pi#pi", "Xi_c+", "#Xi_{c}^{+}", "#Xi^{-}#pi^{+}#pi^{+}"}}}; - if (particles.find(particleName.c_str()) == particles.end()) { + if (particles.find(particleName) == particles.end()) { throw std::runtime_error("ERROR: only Dplus, D0, Ds, LcToPKPi, LcToPK0s, Dstar and XicToXiPiPi particles supported! Exit"); } const auto& particle = particles[particleName.c_str()]; const std::string massAxisTitle = "#it{M}(" + particle.decayProducts + ") (GeV/#it{c}^{2})"; const double massPDG = TDatabasePDG::Instance()->GetParticle(particle.pdgName.c_str())->Mass(); - const std::vector plotLabels = {(particle.decayFormulaLhs + " #rightarrow " + particle.decayFormulaRhs + " + c.c.").c_str(), collisionSystem}; + const std::vector plotLabels = {particle.decayFormulaLhs + " #rightarrow " + particle.decayFormulaRhs + " + c.c.", collisionSystem}; // load inv-mass histograms auto* inputFile = openFileWithNullptrCheck(inputFileName); @@ -395,14 +393,14 @@ void runMassFitter(const std::string& configFileName) NFitResultsToSave }; auto* hFitConfig = new TH2F("hFitConfig", "Fit Configurations", NConfigsToSave - 1, 0, NConfigsToSave - 1, nHistograms, sliceVarLimits.data()); - const char* hFitConfigXLabel[NConfigsToSave - 1] = {"mass min", "mass max", "rebin num", "fix sigma", "bkg func", "sgn func", "rnd seed"}; + constexpr std::array HFitConfigXLabel = {"mass min", "mass max", "rebin num", "fix sigma", "bkg func", "sgn func", "rnd seed"}; auto* hFitResult = new TH2F("hFitResult", "Fit Result", NFitResultsToSave - 1, 0, NFitResultsToSave - 1, nHistograms, sliceVarLimits.data()); - const char* hFitResultXLabel[NFitResultsToSave - 1] = {"status", "cov qual", "edm", "minNLL", "N Sig GCC"}; + constexpr std::array HFitResultXLabel = {"status", "cov qual", "edm", "minNLL", "N Sig GCC"}; for (int i = 0; i < NConfigsToSave - 1; i++) { - hFitConfig->GetXaxis()->SetBinLabel(i + 1, hFitConfigXLabel[i]); + hFitConfig->GetXaxis()->SetBinLabel(i + 1, HFitConfigXLabel[i]); } for (int i = 0; i < NFitResultsToSave - 1; i++) { - hFitResult->GetXaxis()->SetBinLabel(i + 1, hFitResultXLabel[i]); + hFitResult->GetXaxis()->SetBinLabel(i + 1, HFitResultXLabel[i]); } for (const auto& h : {hFitConfig, hFitResult}) { h->SetStats(false); @@ -451,7 +449,7 @@ void runMassFitter(const std::string& configFileName) TH1* hSecondSigmaToFix = getHistToFix(fixSecondSigma, fixSecondSigmaManual, secondSigmaFile, "SecSigma"); TH1* hFracDoubleGausToFix = getHistToFix(fixFracDoubleGaus, fixFracDoubleGausManual, fracDoubleGausFile, "FracDoubleGaus"); - int canvasSize[2] = {1920, 1080}; + std::array canvasSize = {1920, 1080}; if (nHistograms == 1) { canvasSize[0] = 500; canvasSize[1] = 500; @@ -511,7 +509,7 @@ void runMassFitter(const std::string& configFileName) double reflOverSgn = 0; - HFInvMassFitter* massFitter = new HFInvMassFitter(hMass[iSliceVar], massMin[iSliceVar], massMax[iSliceVar], bkgFunc[iSliceVar], sgnFunc[iSliceVar], randomSeed); + auto* massFitter = new HFInvMassFitter(hMass[iSliceVar], massMin[iSliceVar], massMax[iSliceVar], bkgFunc[iSliceVar], sgnFunc[iSliceVar], randomSeed); massFitter->setDrawBgPrefit(drawBgPrefit); massFitter->setNumberOfSigmaForSidebands(nSigmaForSideband); massFitter->setNumberOfSigmaForSignal(nSigmaForSignal); @@ -525,27 +523,27 @@ void runMassFitter(const std::string& configFileName) massFitter->setUseChi2Fit(); } - auto setFixedValue = [&iSliceVar](bool const& isFix, std::vector const& fixManual, const TH1* histToFix, std::function setFunc, std::string const& var) -> void { + auto setFixedValue = [&iSliceVar, massFitter](bool const& isFix, std::vector const& fixManual, const TH1* histToFix, void (HFInvMassFitter::*setter)(double), std::string const& var) -> void { if (isFix) { if (fixManual.empty() && histToFix == nullptr) { throw std::runtime_error("Histogram to fix " + var + " is null while isFix==true and fixManual is empty"); } const auto valueToFix = fixManual.empty() ? histToFix->GetBinContent(iSliceVar + 1) : fixManual[iSliceVar]; - setFunc(valueToFix); + (massFitter->*setter)(valueToFix); printf("*****************************\n"); printf("FIXED %s: %f\n", var.data(), valueToFix); printf("*****************************\n"); } }; - setFixedValue(fixMean, fixMeanManual, hMeanToFix, std::bind(&HFInvMassFitter::setFixGaussianMean, massFitter, std::placeholders::_1), "MEAN"); - setFixedValue(fixSigma, fixSigmaManual, hSigmaToFix, std::bind(&HFInvMassFitter::setFixGaussianSigma, massFitter, std::placeholders::_1), "SIGMA"); - setFixedValue(fixSecondSigma, fixSecondSigmaManual, hSecondSigmaToFix, std::bind(&HFInvMassFitter::setFixSecondGaussianSigma, massFitter, std::placeholders::_1), "SECOND SIGMA"); - setFixedValue(fixFracDoubleGaus, fixFracDoubleGausManual, hFracDoubleGausToFix, std::bind(&HFInvMassFitter::setFixFrac2Gaus, massFitter, std::placeholders::_1), "FRAC DOUBLE GAUS"); - setFixedValue(fixDscbTailParams, dscbAlphaLInitial, nullptr, std::bind(&HFInvMassFitter::setFixDscbAlphaL, massFitter, std::placeholders::_1), "DSCB ALPHA LEFT"); - setFixedValue(fixDscbTailParams, dscbAlphaRInitial, nullptr, std::bind(&HFInvMassFitter::setFixDscbAlphaR, massFitter, std::placeholders::_1), "DSCB ALPHA RIGHT"); - setFixedValue(fixDscbTailParams, dscbNLInitial, nullptr, std::bind(&HFInvMassFitter::setFixDscbNL, massFitter, std::placeholders::_1), "DSCB N LEFT"); - setFixedValue(fixDscbTailParams, dscbNRInitial, nullptr, std::bind(&HFInvMassFitter::setFixDscbNR, massFitter, std::placeholders::_1), "DSCB N RIGHT"); + setFixedValue(fixMean, fixMeanManual, hMeanToFix, &HFInvMassFitter::setFixGaussianMean, "MEAN"); + setFixedValue(fixSigma, fixSigmaManual, hSigmaToFix, &HFInvMassFitter::setFixGaussianSigma, "SIGMA"); + setFixedValue(fixSecondSigma, fixSecondSigmaManual, hSecondSigmaToFix, &HFInvMassFitter::setFixSecondGaussianSigma, "SECOND SIGMA"); + setFixedValue(fixFracDoubleGaus, fixFracDoubleGausManual, hFracDoubleGausToFix, &HFInvMassFitter::setFixFrac2Gaus, "FRAC DOUBLE GAUS"); + setFixedValue(fixDscbTailParams, dscbAlphaLInitial, nullptr, &HFInvMassFitter::setFixDscbAlphaL, "DSCB ALPHA LEFT"); + setFixedValue(fixDscbTailParams, dscbAlphaRInitial, nullptr, &HFInvMassFitter::setFixDscbAlphaR, "DSCB ALPHA RIGHT"); + setFixedValue(fixDscbTailParams, dscbNLInitial, nullptr, &HFInvMassFitter::setFixDscbNL, "DSCB N LEFT"); + setFixedValue(fixDscbTailParams, dscbNRInitial, nullptr, &HFInvMassFitter::setFixDscbNR, "DSCB N RIGHT"); if (!isMc && enableRefl) { reflOverSgn = hMassSgn[iSliceVar]->Integral(hMassSgn[iSliceVar]->FindBin(massMin[iSliceVar] * 1.0001), hMassSgn[iSliceVar]->FindBin(massMax[iSliceVar] * 0.999)); @@ -578,7 +576,7 @@ void runMassFitter(const std::string& configFileName) printf("Warinig! Exception \"%s\" caught while doing fit of the histogram no. %d. The fitting process will be continued without it.\n", e.what(), iSliceVar); } - auto drawOnCanvas = [&](std::vector& canvas, std::function drawer) { + auto drawOnCanvas = [&](std::vector& canvas, const std::function& drawer) { if (nHistograms > 1) { canvas[iCanvas]->cd(iSliceVar - NCanvasesMax * iCanvas + 1); } else { @@ -853,9 +851,8 @@ void readJsonVectorFlexible(std::vector& vec, const Document& config, int nHi if (!config.HasMember(fieldName.c_str())) { if (isRequired) { throw std::runtime_error("readJsonVectorFlexible(): missing required field " + fieldName); - } else { - return; } + return; } if (config[fieldName.c_str()].IsArray()) { readJsonVector(vec, config, fieldName); From d181de26fb4f160757735fdac5921df84b075e77 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 5 Aug 2026 13:46:25 +0200 Subject: [PATCH 32/34] handle exception caught by main() --- PWGHF/D2H/Macros/runMassFitter.C | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index eec8983d9c8..76c7b8a045f 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -889,7 +889,11 @@ int main(int argc, const char* argv[]) const std::string configFileName = argv[1]; - runMassFitter(configFileName); - - return 0; + try { + runMassFitter(configFileName); + return 0; + } catch (const std::exception& e) { + printf("Error: Exception \"%s\" caught during runMassFitter() call. Exit.\n", e.what()); + return 1; + } } From 0336f8632915183bb3890ffda886d1aabf9f9ed2 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 5 Aug 2026 15:42:39 +0200 Subject: [PATCH 33/34] extend hFitResult business on MC case --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 3bcf783ba42..b6987efe7bd 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -228,14 +228,15 @@ void HFInvMassFitter::doFit() const double integralHisto = integrateHistoInvMassOverWorkspaceRanges({"full"}); // fit MC or Data + RooFitResult* fitResult{nullptr}; if (mTypeOfBkgPdf == NoBkg) { // MC const ParameterRanges rooNSgnParamRanges{0.5 * integralHisto, 1.5 * integralHisto, integralHisto}; mRooNSgn = new RooRealVar("mRooNSig", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // signal yield mTotalPdf = new RooAddPdf("mMCFunc", "MC fit function", RooArgList(*sgnPdf), RooArgList(*mRooNSgn)); // create total pdf if (strcmp(mFitOption.c_str(), "Chi2") == 0) { - mTotalPdf->chi2FitTo(dataHistogram, Range("full")); + fitResult = mTotalPdf->chi2FitTo(dataHistogram, Range("full"), Save()); } else { - mTotalPdf->fitTo(dataHistogram, Range("full")); + fitResult = mTotalPdf->fitTo(dataHistogram, Range("full"), Save()); } RooAbsReal* signalIntegralMc = mTotalPdf->createIntegral(*mass, NormSet(*mass), Range("signal")); // sig yield from fit mIntegralSgn = signalIntegralMc->getValV(); @@ -276,7 +277,7 @@ void HFInvMassFitter::doFit() calculateBackground(mBkgYield, mBkgYieldErr); // BG's absolute integral in "bkg" range const ParameterRanges rooNSgnParamRanges{0.1 * estimatedSignal, 10 * estimatedSignal, estimatedSignal}; - mRooNSgn = new RooRealVar("mNSgn", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // estimated signal yield + mRooNSgn = new RooRealVar("mRooNSig", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // estimated signal yield if (mFixedRawYield > 0) { mRooNSgn->setVal(mFixedRawYield); // fixed signal yield mRooNSgn->setConstant(true); @@ -308,20 +309,12 @@ void HFInvMassFitter::doFit() mTotalPdf = new RooAddPdf("mTotalPdf", "background + signal pdf", RooArgList(*bkgPdf, *sgnPdf), RooArgList(*mRooNBkg, *mRooNSgn)); } - RooFitResult* fitResult{nullptr}; if (strcmp(mFitOption.c_str(), "Chi2") == 0) { fitResult = mTotalPdf->chi2FitTo(dataHistogram, Save()); } else { fitResult = mTotalPdf->fitTo(dataHistogram, Save()); } - mFitStatus = fitResult->status(); - mCovQual = fitResult->covQual(); - mEdm = fitResult->edm(); - mMinNll = fitResult->minNll(); - mSgnGlobalCorrelCoeff = fitResult->globalCorr("mNSgn"); - mCovCorrMatrix = fillCovCorrMatrix(fitResult); - plotBkg(mTotalPdf); mTotalPdf->plotOn(mInvMassFrame, Name("Tot_c"), LineColor(kBlue)); if (mHistoTemplateRefl != nullptr) { @@ -353,6 +346,12 @@ void HFInvMassFitter::doFit() mRatioFrame = mass->frame(Title(Form("%s", mHistoInvMass->GetTitle()))); calculateFitToDataRatio(); } + mFitStatus = fitResult->status(); + mCovQual = fitResult->covQual(); + mEdm = fitResult->edm(); + mMinNll = fitResult->minNll(); + mSgnGlobalCorrelCoeff = fitResult->globalCorr("mRooNSig"); + mCovCorrMatrix = fillCovCorrMatrix(fitResult); } void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const From c3737d3086e9e8270cb794224135f5c1c6888d9e Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Thu, 6 Aug 2026 01:22:47 +0200 Subject: [PATCH 34/34] fix code-check errors --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 3 ++- PWGHF/D2H/Macros/runMassFitter.C | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index b6987efe7bd..5f80f2cd444 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -42,6 +42,7 @@ #include #include #include +#include #include #include #include @@ -264,7 +265,7 @@ void HFInvMassFitter::doFit() mChiSquareOverNdfBkg = frameTemporary->chiSquare("Bkg_sidebands", "data_for_bkgchi2"); // calculate reduced chi2 / NDF of background sidebands (pre-fit) delete frameTemporary; if (mDrawBgPrefit) { - RooAbsPdf* bkgPdfPrefit = dynamic_cast(mBkgPdf->Clone()); + auto* bkgPdfPrefit = dynamic_cast(mBkgPdf->Clone()); bkgPdfPrefit->plotOn(mInvMassFrame, Range("full"), Normalization(mRooNBkg->getVal(), RooAbsReal::NumEvent), Name("Bkg_c_prefit"), LineColor(kGray)); delete bkgPdfPrefit; } diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 76c7b8a045f..6ed3ccee7c1 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -41,6 +41,7 @@ #include #include #include +#include #include #include #include