From 43fea073941f352176e9424d0813d4f88adacbac Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Sat, 18 Apr 2026 23:13:32 +0200 Subject: [PATCH 01/26] resolve conflict with master --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 33 +++++++++++++++++++++++++++- PWGHF/D2H/Macros/HFInvMassFitter.h | 12 ++++++++++ 2 files changed, 44 insertions(+), 1 deletion(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index e56c9b87a28..d2c2155dce5 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -39,6 +39,7 @@ #include #include #include +#include #include #include #include @@ -129,12 +130,18 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mIntegralSgn(0), mHistoTemplateRefl(nullptr), mDrawBgPrefit(false), - mHighlightPeakRegion(false) + mHighlightPeakRegion(false), + mRandomSeed(-1), + mRandomGen(nullptr) { // standard constructor mHistoInvMass = histoToFit; mHistoInvMass->SetName("mHistoInvMass"); mHistoInvMass->SetDirectory(nullptr); + if (mRandomSeed >= 0) { + mRandomGen = new TRandom3(); + mRandomGen->SetSeed(mRandomSeed); + } } HFInvMassFitter::~HFInvMassFitter() @@ -991,3 +998,27 @@ void HFInvMassFitter::setTemplateReflections(TH1* histoRefl) mHistoTemplateRefl = histoRefl; mHistoTemplateRefl->SetName("mHistoTemplateRefl"); } + +double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& parameterRanges) +{ + constexpr double DefaultSigmaFraction{10.}; + constexpr int MaximalNumberOfIterations{20}; + + if (mRandomSeed < 0) { + return parameterRanges.initial; + } + const auto sigma = parameterRanges.sigma < 0 ? (parameterRanges.upper - parameterRanges.lower) / DefaultSigmaFraction : parameterRanges.sigma; + + double result; + int nIter{0}; + do { + result = mRandomGen->Gaus(parameterRanges.initial, sigma); + ++nIter; + if (nIter > MaximalNumberOfIterations) { + printf("randomizeInitialFitParameter() - long while loop with lower = %f upper = %f initial = %f sigma = %f\n", parameterRanges.lower, parameterRanges.upper, parameterRanges.initial, sigma); + throw std::runtime_error(""); + } + } while (result < parameterRanges.lower || result > parameterRanges.upper); + + return result; +} diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 28ee858bba6..a3982817418 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -28,6 +28,7 @@ #include #include #include +#include #include #include @@ -39,6 +40,14 @@ class HFInvMassFitter : public TNamed { + private: + struct ParameterRanges { + double lower{}; + double upper{}; + double initial{}; + double sigma{-999.}; + }; + public: enum TypeOfBkgPdf { Expo = 0, @@ -141,6 +150,7 @@ class HFInvMassFitter : public TNamed HFInvMassFitter& operator=(const HFInvMassFitter& source); void fillWorkspace(RooWorkspace& w) const; void highlightPeakRegion(const RooPlot* plot, Color_t color = kGray + 1, Width_t width = 1, Style_t style = 2) const; + double randomizeInitialParameter(const ParameterRanges& parameterRanges); TH1* mHistoInvMass; // histogram to fit std::string mFitOption; @@ -212,6 +222,8 @@ class HFInvMassFitter : public TNamed 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 ClassDefOverride(HFInvMassFitter, 1); }; From 38b036743db99679e74df35e8cefb85e173fef71 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Fri, 16 Jan 2026 22:10:48 +0100 Subject: [PATCH 02/26] random seed: setter, read from json --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 2 +- PWGHF/D2H/Macros/HFInvMassFitter.h | 1 + PWGHF/D2H/Macros/config_massfitter.json | 3 ++- PWGHF/D2H/Macros/runMassFitter.C | 2 ++ 4 files changed, 6 insertions(+), 2 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index d2c2155dce5..041844fcb33 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -1016,7 +1016,7 @@ double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& paramet ++nIter; if (nIter > MaximalNumberOfIterations) { printf("randomizeInitialFitParameter() - long while loop with lower = %f upper = %f initial = %f sigma = %f\n", parameterRanges.lower, parameterRanges.upper, parameterRanges.initial, sigma); - throw std::runtime_error(""); + throw; } } while (result < parameterRanges.lower || result > parameterRanges.upper); diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index a3982817418..9fdaf397dfd 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -144,6 +144,7 @@ class HFInvMassFitter : public TNamed void drawResidual(TVirtualPad* c); void drawRatio(TVirtualPad* c); void drawReflection(TVirtualPad* c); + void setRandomSeed(int seed) { mRandomSeed = seed; } private: HFInvMassFitter(const HFInvMassFitter& source); diff --git a/PWGHF/D2H/Macros/config_massfitter.json b/PWGHF/D2H/Macros/config_massfitter.json index 25dc9d50bcc..08e3aa21a5a 100644 --- a/PWGHF/D2H/Macros/config_massfitter.json +++ b/PWGHF/D2H/Macros/config_massfitter.json @@ -158,5 +158,6 @@ "3 for GausSec" ], "DrawBgPrefit": true, - "HighlightPeakRegion": true + "HighlightPeakRegion": true, + "RandomSeed": -1 } diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 5937d12e0ba..d81b3c755a5 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -153,6 +153,7 @@ void runMassFitter(const std::string& configFileName) 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 int nSliceVarBins = static_cast(sliceVarMin.size()); std::vector sliceVarLimits(nSliceVarBins + 1); @@ -380,6 +381,7 @@ void runMassFitter(const std::string& configFileName) } else { massFitter->setUseChi2Fit(); } + massFitter->setRandomSeed(randomSeed); auto setFixedValue = [&iSliceVar](bool const& isFix, std::vector const& fixManual, const TH1* histToFix, std::function setFunc, std::string const& var) -> void { if (isFix) { From a95b7a4eea13531598604c7d1a1613163ade5c28 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Sun, 18 Jan 2026 22:20:41 +0100 Subject: [PATCH 03/26] randomize initial parameters --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 66 ++++++++++++++++++---------- PWGHF/D2H/Macros/HFInvMassFitter.h | 2 +- 2 files changed, 43 insertions(+), 25 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 041844fcb33..dfcd15be123 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -196,14 +196,16 @@ void HFInvMassFitter::doFit() dataHistogram.plotOn(mInvMassFrame, Name("data_c")); // plot data histogram on the frame // define number of background and background fit function - mRooNBkg = new RooRealVar("mRooNBkg", "number of background", 0.3 * mIntegralHisto, 0., 1.2 * mIntegralHisto); // background yield - RooAbsPdf* bkgPdf = createBackgroundFitFunction(mWorkspace); // Create background pdf - RooAbsPdf* sgnPdf = createSignalFitFunction(mWorkspace); // Create signal pdf + const ParameterRanges rooNBkgParamRanges{0., 1.2 * mIntegralHisto, 0.3 * mIntegralHisto}; + mRooNBkg = new RooRealVar("mRooNBkg", "number of background", randomizeInitialParameter(rooNBkgParamRanges), rooNBkgParamRanges.lower, rooNBkgParamRanges.upper); // background yield + RooAbsPdf* bkgPdf = createBackgroundFitFunction(mWorkspace); // Create background pdf + RooAbsPdf* sgnPdf = createSignalFitFunction(mWorkspace); // Create signal pdf // fit MC or Data - if (mTypeOfBkgPdf == NoBkg) { // MC - mRooNSgn = new RooRealVar("mRooNSig", "number of signal", 0.3 * mIntegralHisto, 0., 1.2 * mIntegralHisto); // signal yield - mTotalPdf = new RooAddPdf("mMCFunc", "MC fit function", RooArgList(*sgnPdf), RooArgList(*mRooNSgn)); // create total pdf + if (mTypeOfBkgPdf == NoBkg) { // MC + const ParameterRanges rooNSgnParamRanges{0., 1.2 * mIntegralHisto, 0.3 * mIntegralHisto}; + 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("signal")); } else { @@ -250,7 +252,8 @@ void HFInvMassFitter::doFit() checkForSignal(estimatedSignal); // SIG's absolute integral in "bkg" range calculateBackground(mBkgYield, mBkgYieldErr); // BG's absolute integral in "bkg" range - mRooNSgn = new RooRealVar("mNSgn", "number of signal", 0.3 * estimatedSignal, 0., 1.2 * estimatedSignal); // estimated signal yield + const ParameterRanges rooNSgnParamRanges{0., 1.2 * estimatedSignal, 0.3 * 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 mRooNSgn->setConstant(true); @@ -263,7 +266,8 @@ void HFInvMassFitter::doFit() mReflFrame = mass->frame(); mReflOnlyFrame = mass->frame(Title(Form("%s", mHistoTemplateRefl->GetTitle()))); reflHistogram.plotOn(mReflOnlyFrame); - mRooNRefl = new RooRealVar("mNRefl", "number of reflection", 0.5 * mHistoTemplateRefl->Integral(), 0, mHistoTemplateRefl->Integral()); + const ParameterRanges rooNReflParamRanges{0., mHistoTemplateRefl->Integral(), 0.5 * mHistoTemplateRefl->Integral()}; + mRooNRefl = new RooRealVar("mNRefl", "number of reflection", randomizeInitialParameter(rooNReflParamRanges), rooNReflParamRanges.lower, rooNReflParamRanges.upper); RooAddPdf reflFuncTemp("reflFuncTemp", "template reflection fit function", RooArgList(*reflPdf), RooArgList(*mRooNRefl)); if (strcmp(mFitOption.c_str(), "Chi2") == 0) { reflFuncTemp.chi2FitTo(reflHistogram); @@ -333,35 +337,42 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const // Declare observable variable RooRealVar mass("mass", "mass", mMinMass, mMaxMass, "GeV/c^{2}"); // bkg expo - RooRealVar tau("tau", "tau", -1, -5., 5.); + const ParameterRanges tauParamRanges{-5., 5., -1., 0.1}; + RooRealVar tau("tau", "tau", randomizeInitialParameter(tauParamRanges), tauParamRanges.lower, tauParamRanges.upper); RooAbsPdf* bkgFuncExpo = new RooExponential("bkgFuncExpo", "background fit function", mass, tau); workspace.import(*bkgFuncExpo); delete bkgFuncExpo; // bkg poly1 - RooRealVar const polyParam0("polyParam0", "Parameter of Poly function", 0.5, -5., 5.); - RooRealVar const polyParam1("polyParam1", "Parameter of Poly function", 0.2, -5., 5.); + const ParameterRanges polyParam0ParamRanges{-5., 5., 0.5, 0.1}; + RooRealVar const polyParam0("polyParam0", "Parameter of Poly function", randomizeInitialParameter(polyParam0ParamRanges), polyParam0ParamRanges.lower, polyParam0ParamRanges.upper); + const ParameterRanges polyParam1ParamRanges{-5., 5., 0.2, 0.05}; + RooRealVar const polyParam1("polyParam1", "Parameter of Poly function", randomizeInitialParameter(polyParam1ParamRanges), polyParam1ParamRanges.lower, polyParam1ParamRanges.upper); RooAbsPdf* bkgFuncPoly1 = new RooPolynomial("bkgFuncPoly1", "background fit function", mass, RooArgSet(polyParam0, polyParam1)); workspace.import(*bkgFuncPoly1); delete bkgFuncPoly1; // bkg poly2 - RooRealVar const polyParam2("polyParam2", "Parameter of Poly function", 0.2, -5., 5.); + const ParameterRanges polyParam2ParamRanges{-5., 5., 0.2, 0.05}; + RooRealVar const polyParam2("polyParam2", "Parameter of Poly function", randomizeInitialParameter(polyParam2ParamRanges), polyParam2ParamRanges.lower, polyParam2ParamRanges.upper); RooAbsPdf* bkgFuncPoly2 = new RooPolynomial("bkgFuncPoly2", "background fit function", mass, RooArgSet(polyParam0, polyParam1, polyParam2)); workspace.import(*bkgFuncPoly2); delete bkgFuncPoly2; // bkg poly3 - RooRealVar const polyParam3("polyParam3", "Parameter of Poly function", 0.2, -1., 1.); + const ParameterRanges polyParam3ParamRanges{-1., 1., 0.2, 0.05}; + RooRealVar const polyParam3("polyParam3", "Parameter of Poly function", randomizeInitialParameter(polyParam3ParamRanges), polyParam3ParamRanges.lower, polyParam3ParamRanges.upper); RooAbsPdf* bkgFuncPoly3 = new RooPolynomial("bkgFuncPoly3", "background pdf", mass, RooArgSet(polyParam0, polyParam1, polyParam2, polyParam3)); workspace.import(*bkgFuncPoly3); delete bkgFuncPoly3; // bkg power law RooRealVar const powParam1("powParam1", "Parameter of Pow function", TDatabasePDG::Instance()->GetParticle("pi+")->Mass()); - RooRealVar const powParam2("powParam2", "Parameter of Pow function", 1., -10, 10); + const ParameterRanges powParam2ParamRanges{-10., 10., 1., 0.2}; + RooRealVar const powParam2("powParam2", "Parameter of Pow function", randomizeInitialParameter(powParam2ParamRanges), powParam2ParamRanges.lower, powParam2ParamRanges.upper); RooAbsPdf* bkgFuncPow = new RooGenericPdf("bkgFuncPow", "bkgFuncPow", "(mass-powParam1)^powParam2", RooArgSet(mass, powParam1, powParam2)); workspace.import(*bkgFuncPow); delete bkgFuncPow; // pow * exp RooRealVar const powExpoParam1("powExpoParam1", "Parameter of PowExpo function", 1. / 2.); - RooRealVar const powExpoParam2("powExpoParam2", "Parameter of PowExpo function", 1, -10, 10); + const ParameterRanges powExpoParam2ParamRanges{-10., 10., 1., 0.2}; + RooRealVar const powExpoParam2("powExpoParam2", "Parameter of PowExpo function", randomizeInitialParameter(powExpoParam2ParamRanges), powExpoParam2ParamRanges.lower, powExpoParam2ParamRanges.upper); RooRealVar massPi("massPi", "mass of pion", TDatabasePDG::Instance()->GetParticle("pi+")->Mass()); RooFormulaVar powExpoParam3("powExpoParam3", "powExpoParam1 + 1", RooArgList(powExpoParam1)); RooFormulaVar powExpoParam4("powExpoParam4", "1./powExpoParam2", RooArgList(powExpoParam2)); @@ -490,17 +501,24 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const workspace.import(*reflFuncDoubleGaus); delete reflFuncDoubleGaus; // reflection poly3 - RooRealVar const polyReflParam0("polyReflParam0", "polyReflParam0", 0.5, -1., 1.); - RooRealVar const polyReflParam1("polyReflParam1", "polyReflParam1", 0.2, -1., 1.); - RooRealVar const polyReflParam2("polyReflParam2", "polyReflParam2", 0.2, -1., 1.); - RooRealVar const polyReflParam3("polyReflParam3", "polyReflParam3", 0.2, -1., 1.); + const ParameterRanges polyReflParam0ParamRanges{-1., 1., 0.5, 0.1}; + RooRealVar const polyReflParam0("polyReflParam0", "polyReflParam0", randomizeInitialParameter(polyReflParam0ParamRanges), polyReflParam0ParamRanges.lower, polyReflParam0ParamRanges.upper); + const ParameterRanges polyReflParam1ParamRanges{-1., 1., 0.2, 0.05}; + RooRealVar const polyReflParam1("polyReflParam1", "polyReflParam1", randomizeInitialParameter(polyReflParam1ParamRanges), polyReflParam1ParamRanges.lower, polyReflParam1ParamRanges.upper); + const ParameterRanges polyReflParam2ParamRanges{-1., 1., 0.2, 0.05}; + RooRealVar const polyReflParam2("polyReflParam2", "polyReflParam2", randomizeInitialParameter(polyReflParam2ParamRanges), polyReflParam2ParamRanges.lower, polyReflParam2ParamRanges.upper); + const ParameterRanges polyReflParam3ParamRanges{-1., 1., 0.2, 0.05}; + RooRealVar const polyReflParam3("polyReflParam3", "polyReflParam3", randomizeInitialParameter(polyReflParam3ParamRanges), polyReflParam3ParamRanges.lower, polyReflParam3ParamRanges.upper); RooAbsPdf* reflFuncPoly3 = new RooPolynomial("reflFuncPoly3", "reflection PDF", mass, RooArgSet(polyReflParam0, polyReflParam1, polyReflParam2, polyReflParam3)); workspace.import(*reflFuncPoly3); delete reflFuncPoly3; // reflection poly6 - RooRealVar const polyReflParam4("polyReflParam4", "polyReflParam4", 0.2, -1., 1.); - RooRealVar const polyReflParam5("polyReflParam5", "polyReflParam5", 0.2, -1., 1.); - RooRealVar const polyReflParam6("polyReflParam6", "polyReflParam6", 0.2, -1., 1.); + const ParameterRanges polyReflParam4ParamRanges{-1., 1., 0.2, 0.05}; + RooRealVar const polyReflParam4("polyReflParam4", "polyReflParam4", randomizeInitialParameter(polyReflParam4ParamRanges), polyReflParam4ParamRanges.lower, polyReflParam4ParamRanges.upper); + const ParameterRanges polyReflParam5ParamRanges{-1., 1., 0.2, 0.05}; + RooRealVar const polyReflParam5("polyReflParam5", "polyReflParam5", randomizeInitialParameter(polyReflParam5ParamRanges), polyReflParam5ParamRanges.lower, polyReflParam5ParamRanges.upper); + const ParameterRanges polyReflParam6ParamRanges{-1., 1., 0.2, 0.05}; + RooRealVar const polyReflParam6("polyReflParam6", "polyReflParam6", randomizeInitialParameter(polyReflParam6ParamRanges), polyReflParam6ParamRanges.lower, polyReflParam6ParamRanges.upper); RooAbsPdf* reflFuncPoly6 = new RooPolynomial("reflFuncPoly6", "reflection pdf", mass, RooArgSet(polyReflParam0, polyReflParam1, polyReflParam2, polyReflParam3, polyReflParam4, polyReflParam5, polyReflParam6)); workspace.import(*reflFuncPoly6); delete reflFuncPoly6; @@ -999,7 +1017,7 @@ void HFInvMassFitter::setTemplateReflections(TH1* histoRefl) mHistoTemplateRefl->SetName("mHistoTemplateRefl"); } -double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& parameterRanges) +double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& parameterRanges) const { constexpr double DefaultSigmaFraction{10.}; constexpr int MaximalNumberOfIterations{20}; @@ -1015,7 +1033,7 @@ double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& paramet result = mRandomGen->Gaus(parameterRanges.initial, sigma); ++nIter; if (nIter > MaximalNumberOfIterations) { - printf("randomizeInitialFitParameter() - long while loop with lower = %f upper = %f initial = %f sigma = %f\n", parameterRanges.lower, parameterRanges.upper, parameterRanges.initial, sigma); + printf("randomizeInitialParameter() - long while loop with lower = %f upper = %f initial = %f sigma = %f\n", parameterRanges.lower, parameterRanges.upper, parameterRanges.initial, sigma); throw; } } while (result < parameterRanges.lower || result > parameterRanges.upper); diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 9fdaf397dfd..8e7e23c6963 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -151,7 +151,7 @@ class HFInvMassFitter : public TNamed HFInvMassFitter& operator=(const HFInvMassFitter& source); void fillWorkspace(RooWorkspace& w) const; void highlightPeakRegion(const RooPlot* plot, Color_t color = kGray + 1, Width_t width = 1, Style_t style = 2) const; - double randomizeInitialParameter(const ParameterRanges& parameterRanges); + double randomizeInitialParameter(const ParameterRanges& parameterRanges) const; TH1* mHistoInvMass; // histogram to fit std::string mFitOption; From 6859e4a5921c51b5c419b87b6e0ea31f95f3e81e Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Mon, 19 Jan 2026 18:47:41 +0100 Subject: [PATCH 04/26] DSCB: add class members and member functions (set, get, fix) --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 50 ++++++++++++++++++++++++++++ PWGHF/D2H/Macros/HFInvMassFitter.h | 42 +++++++++++++++++++++++ 2 files changed, 92 insertions(+) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index dfcd15be123..af2dc90b007 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -87,6 +87,7 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mFixedSigma(false), mFixedSigmaDoubleGaus(false), mBoundSigma(false), + mFixedDscbTailParams(false), mSigmaValue(0.012), mParamSgn(0.1), mFracDoubleGaus(0.2), @@ -107,6 +108,18 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mChiSquareOverNdfTotal(0), mChiSquareOverNdfBkg(0), mFixReflOverSgn(false), + mDscbAlphaLInitialValue(1.5), + mDscbAlphaLLowLimit(0.1), + mDscbAlphaLUpLimit(5.0), + mDscbAlphaRInitialValue(1.5), + mDscbAlphaRLowLimit(0.1), + mDscbAlphaRUpLimit(5.0), + mDscbNLInitialValue(2.0), + mDscbNLLowLimit(0.5), + mDscbNLUpLimit(50.), + mDscbNRInitialValue(2.0), + mDscbNRLowLimit(0.5), + mDscbNRUpLimit(50.), mRooMeanSgn(nullptr), mRooSigmaSgn(nullptr), mRooSecSigmaSgn(nullptr), @@ -117,6 +130,10 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mRooNSgn(nullptr), mRooNBkg(nullptr), mRooNRefl(nullptr), + mRooDscbAlphaL(nullptr), + mRooDscbAlphaR(nullptr), + mRooDscbNL(nullptr), + mRooDscbNR(nullptr), mTotalPdf(nullptr), mInvMassFrame(nullptr), mReflFrame(nullptr), @@ -1017,6 +1034,39 @@ void HFInvMassFitter::setTemplateReflections(TH1* histoRefl) mHistoTemplateRefl->SetName("mHistoTemplateRefl"); } +void HFInvMassFitter::setFixDscbAlphaL(double alphaL) +{ + if (mTypeOfSgnPdf != DoubleSidedCrystalBall) { + printf("Fit type should be DSCB!\n"); + } + mFixedDscbTailParams = kTRUE; + mDscbAlphaLInitialValue = alphaL; +} +void HFInvMassFitter::setFixDscbAlphaR(double alphaR) +{ + if (mTypeOfSgnPdf != DoubleSidedCrystalBall) { + printf("Fit type should be DSCB!\n"); + } + mFixedDscbTailParams = kTRUE; + mDscbAlphaRInitialValue = alphaR; +} +void HFInvMassFitter::setFixDscbNL(double nL) +{ + if (mTypeOfSgnPdf != DoubleSidedCrystalBall) { + printf("Fit type should be DSCB!\n"); + } + mFixedDscbTailParams = kTRUE; + mDscbNLInitialValue = nL; +} +void HFInvMassFitter::setFixDscbNR(double nR) +{ + if (mTypeOfSgnPdf != DoubleSidedCrystalBall) { + printf("Fit type should be DSCB!\n"); + } + mFixedDscbTailParams = kTRUE; + mDscbNRInitialValue = nR; +} + double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& parameterRanges) const { constexpr double DefaultSigmaFraction{10.}; diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 8e7e23c6963..5e7133e08e2 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -65,6 +65,7 @@ class HFInvMassFitter : public TNamed DoubleGaus = 1, DoubleGausSigmaRatioPar = 2, GausSec = 3, + DoubleSidedCrystalBall = 4, NTypesOfSgnPdf }; enum TypeOfReflPdf { @@ -106,6 +107,22 @@ class HFInvMassFitter : public TNamed void setFixRatioToGausSigma(double sigmaFrac); void setFixSignalYield(double yield) { mFixedRawYield = yield; } void setNumberOfSigmaForSidebands(double numberOfSigma) { mNSigmaForSidebands = numberOfSigma; } + void setFixDscbAlphaL(double alphaL); + void setFixDscbAlphaR(double alphaR); + void setFixDscbNL(double nL); + void setFixDscbNR(double nR); + void setDscbAlphaLInitialValue(double value) { mDscbAlphaLInitialValue = value; } + void setDscbAlphaLLowLimit(double value) { mDscbAlphaLLowLimit = value; } + void setDscbAlphaLUpLimit(double value) { mDscbAlphaLUpLimit = value; } + void setDscbAlphaRInitialValue(double value) { mDscbAlphaRInitialValue = value; } + void setDscbAlphaRLowLimit(double value) { mDscbAlphaRLowLimit = value; } + void setDscbAlphaRUpLimit(double value) { mDscbAlphaRUpLimit = value; } + void setDscbNLInitialValue(double value) { mDscbNLInitialValue = value; } + void setDscbNLLowLimit(double value) { mDscbNLLowLimit = value; } + void setDscbNLUpLimit(double value) { mDscbNLUpLimit = value; } + 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 setReflFuncFixed(); @@ -129,6 +146,14 @@ class HFInvMassFitter : public TNamed [[nodiscard]] double getMeanUncertainty() const { return mRooMeanSgn->getError(); } [[nodiscard]] double getSigma() const { return mRooSigmaSgn->getVal(); } [[nodiscard]] double getSigmaUncertainty() const { return mRooSigmaSgn->getError(); } + [[nodiscard]] double getDscbAlphaL() const { return mRooDscbAlphaL ? mRooDscbAlphaL->getVal() : 0.; } + [[nodiscard]] double getDscbAlphaR() const { return mRooDscbAlphaR ? mRooDscbAlphaR->getVal() : 0.; } + [[nodiscard]] double getDscbNL() const { return mRooDscbNL ? mRooDscbNL->getVal() : 0.; } + [[nodiscard]] double getDscbNR() const { return mRooDscbNR ? mRooDscbNR->getVal() : 0.; } + [[nodiscard]] double getDscbAlphaLUncertainty() const { return mRooDscbAlphaL ? mRooDscbAlphaL->getError() : 0.; } + [[nodiscard]] double getDscbAlphaRUncertainty() const { return mRooDscbAlphaR ? mRooDscbAlphaR->getError() : 0.; } + [[nodiscard]] double getDscbNLUncertainty() const { return mRooDscbNL ? mRooDscbNL->getError() : 0.; } + [[nodiscard]] double getDscbNRUncertainty() const { return mRooDscbNR ? mRooDscbNR->getError() : 0.; } [[nodiscard]] double getSecSigma() const { return mRooSecSigmaSgn->getVal(); } [[nodiscard]] double getSecSigmaUncertainty() const { return mRooSecSigmaSgn->getError(); } [[nodiscard]] double getFracDoubleGaus() const { return mRooFracDoubleGaus->getVal(); } @@ -179,6 +204,7 @@ class HFInvMassFitter : public TNamed 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 @@ -199,6 +225,18 @@ class HFInvMassFitter : public TNamed 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 @@ -209,6 +247,10 @@ class HFInvMassFitter : public TNamed 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 From 683f0f18f6295ad72f25525fc496493d60d09f37 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Tue, 20 Jan 2026 22:50:46 +0100 Subject: [PATCH 05/26] signal DSCB pdf --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 24 ++++++++++++++++++++++++ PWGHF/D2H/Macros/config_massfitter.json | 3 ++- 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index af2dc90b007..5cd7151d3e4 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -22,6 +22,7 @@ #include "HFInvMassFitter.h" #include +#include #include #include #include @@ -517,6 +518,29 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const RooAbsPdf* reflFuncDoubleGaus = new RooAddPdf("reflFuncDoubleGaus", "reflection pdf", RooArgList(gausRefl1, gausRefl2), fracRefl); workspace.import(*reflFuncDoubleGaus); delete reflFuncDoubleGaus; + // signal DSCB pdf + const ParameterRanges dscbAlphaLParamRanges{mDscbAlphaLLowLimit, mDscbAlphaLUpLimit, mDscbAlphaLInitialValue}; + const ParameterRanges dscbNLParamRanges{mDscbNLLowLimit, mDscbNLUpLimit, mDscbNLInitialValue}; + const ParameterRanges dscbAlphaRParamRanges{mDscbAlphaRLowLimit, mDscbAlphaRUpLimit, mDscbAlphaRInitialValue}; + const ParameterRanges dscbNRParamRanges{mDscbNRLowLimit, mDscbNRUpLimit, mDscbNRInitialValue}; + RooRealVar alphaL("alphaL", "left tail alpha", randomizeInitialParameter(dscbAlphaLParamRanges), dscbAlphaLParamRanges.lower, dscbAlphaLParamRanges.upper); + RooRealVar nL("nL", "left tail n", randomizeInitialParameter(dscbNLParamRanges), dscbNLParamRanges.lower, dscbNLParamRanges.upper); + RooRealVar alphaR("alphaR", "right tail alpha", randomizeInitialParameter(dscbAlphaRParamRanges), dscbAlphaRParamRanges.lower, dscbAlphaRParamRanges.upper); + RooRealVar nR("nR", "right tail n", randomizeInitialParameter(dscbNRParamRanges), dscbNRParamRanges.lower, dscbNRParamRanges.upper); + if (mFixedDscbTailParams) { + printf("Fixing DSCB tail parameters to initial values.\n"); + alphaL.setVal(mDscbAlphaLInitialValue); + alphaL.setConstant(true); + alphaR.setVal(mDscbAlphaRInitialValue); + alphaR.setConstant(true); + nL.setVal(mDscbNLInitialValue); + nL.setConstant(true); + nR.setVal(mDscbNRInitialValue); + nR.setConstant(true); + } + RooAbsPdf* sgnFuncDSCB = new RooCrystalBall("sgnFuncDSCB", "double sided crystal ball signal", mass, mean, sigma, alphaL, nL, alphaR, nR); + workspace.import(*sgnFuncDSCB); + delete sgnFuncDSCB; // reflection poly3 const ParameterRanges polyReflParam0ParamRanges{-1., 1., 0.5, 0.1}; RooRealVar const polyReflParam0("polyReflParam0", "polyReflParam0", randomizeInitialParameter(polyReflParam0ParamRanges), polyReflParam0ParamRanges.lower, polyReflParam0ParamRanges.upper); diff --git a/PWGHF/D2H/Macros/config_massfitter.json b/PWGHF/D2H/Macros/config_massfitter.json index 08e3aa21a5a..7ee52cace9a 100644 --- a/PWGHF/D2H/Macros/config_massfitter.json +++ b/PWGHF/D2H/Macros/config_massfitter.json @@ -155,7 +155,8 @@ "0 for SingleGaus", "1 for DoubleGaus", "2 for DoubleGausSigmaRatioPar", - "3 for GausSec" + "3 for GausSec", + "4 for DoubleSidedCrystalBall" ], "DrawBgPrefit": true, "HighlightPeakRegion": true, From 0feee3898c0bffadb8a7e87344699c6880021837 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 21 Jan 2026 19:02:45 +0100 Subject: [PATCH 06/26] add dscb i/o in runMassFitter --- PWGHF/D2H/Macros/runMassFitter.C | 123 +++++++++++++++++++++++++++++++ 1 file changed, 123 insertions(+) diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index d81b3c755a5..e7ab0258307 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -66,6 +66,8 @@ T readJsonField(const Document& config, const std::string& fieldName); template void readJsonVector(std::vector& vec, const Document& config, const std::string& fieldName, bool isRequired = false); +void readJsonVectorFromHisto(std::vector& vec, const Document& config, const std::string& fileNameFieldName, const std::string& histoNameFieldName); + void divideCanvas(TCanvas* c, int nSliceVarBins); void setHistoStyle(TH1* histo, Color_t color = kBlack, Size_t markerSize = 1); @@ -109,6 +111,18 @@ void runMassFitter(const std::string& configFileName) std::vector nRebin; std::vector bkgFunc; std::vector sgnFunc; + std::vector dscbAlphaLInitial; + std::vector dscbAlphaLLower; + std::vector dscbAlphaLUpper; + std::vector dscbAlphaRInitial; + std::vector dscbAlphaRLower; + std::vector dscbAlphaRUpper; + std::vector dscbNLInitial; + std::vector dscbNLLower; + std::vector dscbNLUpper; + std::vector dscbNRInitial; + std::vector dscbNRLower; + std::vector dscbNRUpper; readJsonVector(inputHistoName, config, "InputHistoName", true); readJsonVector(promptHistoName, config, "PromptHistoName"); @@ -133,6 +147,8 @@ void runMassFitter(const std::string& configFileName) const std::string fracDoubleGausFile = readJsonField(config, "FracDoubleGausFile", ""); readJsonVector(fixFracDoubleGausManual, config, "FixFracDoubleGausManual"); + const bool fixDscbTailParams = readJsonField(config, "FixDscbTailParams", false); + const TString sliceVarName = readJsonField(config, "SliceVarName"); const TString sliceVarUnit = readJsonField(config, "SliceVarUnit"); @@ -155,6 +171,31 @@ void runMassFitter(const std::string& configFileName) const bool highlightPeakRegion = readJsonField(config, "HighlightPeakRegion", true); const int randomSeed = readJsonField(config, "RandomSeed", -1); + readJsonVector(dscbAlphaLInitial, config, "DscbAlphaLInitial"); + readJsonVector(dscbAlphaLLower, config, "DscbAlphaLLower"); + readJsonVector(dscbAlphaLUpper, config, "DscbAlphaLUpper"); + readJsonVector(dscbAlphaRInitial, config, "DscbAlphaRInitial"); + readJsonVector(dscbAlphaRLower, config, "DscbAlphaRLower"); + readJsonVector(dscbAlphaRUpper, config, "DscbAlphaRUpper"); + readJsonVector(dscbNLInitial, config, "DscbNLInitial"); + readJsonVector(dscbNLLower, config, "DscbNLLower"); + readJsonVector(dscbNLUpper, config, "DscbNLUpper"); + readJsonVector(dscbNRInitial, config, "DscbNRInitial"); + readJsonVector(dscbNRLower, config, "DscbNRLower"); + readJsonVector(dscbNRUpper, config, "DscbNRUpper"); + readJsonVectorFromHisto(dscbAlphaLInitial, config, "DscbParametersFile", "DscbAlphaLInitialHisto"); + readJsonVectorFromHisto(dscbAlphaLLower, config, "DscbParametersFile", "DscbAlphaLLowerHisto"); + readJsonVectorFromHisto(dscbAlphaLUpper, config, "DscbParametersFile", "DscbAlphaLUpperHisto"); + readJsonVectorFromHisto(dscbAlphaRInitial, config, "DscbParametersFile", "DscbAlphaRInitialHisto"); + readJsonVectorFromHisto(dscbAlphaRLower, config, "DscbParametersFile", "DscbAlphaRLowerHisto"); + readJsonVectorFromHisto(dscbAlphaRUpper, config, "DscbParametersFile", "DscbAlphaRUpperHisto"); + readJsonVectorFromHisto(dscbNLInitial, config, "DscbParametersFile", "DscbNLInitialHisto"); + readJsonVectorFromHisto(dscbNLLower, config, "DscbParametersFile", "DscbNLLowerHisto"); + readJsonVectorFromHisto(dscbNLUpper, config, "DscbParametersFile", "DscbNLUpperHisto"); + readJsonVectorFromHisto(dscbNRInitial, config, "DscbParametersFile", "DscbNRInitialHisto"); + readJsonVectorFromHisto(dscbNRLower, config, "DscbParametersFile", "DscbNRLowerHisto"); + readJsonVectorFromHisto(dscbNRUpper, config, "DscbParametersFile", "DscbNRUpperHisto"); + const int nSliceVarBins = static_cast(sliceVarMin.size()); std::vector sliceVarLimits(nSliceVarBins + 1); @@ -184,6 +225,18 @@ void runMassFitter(const std::string& configFileName) checkVectorSize(nRebin, "nRebin"); checkVectorSize(bkgFunc, "bkgFunc"); checkVectorSize(sgnFunc, "sgnFunc"); + checkVectorSize(dscbAlphaLInitial, "dscbAlphaLInitial", true); + checkVectorSize(dscbAlphaLLower, "dscbAlphaLLower", true); + checkVectorSize(dscbAlphaLUpper, "dscbAlphaLUpper", true); + checkVectorSize(dscbAlphaRInitial, "dscbAlphaRInitial", true); + checkVectorSize(dscbAlphaRLower, "dscbAlphaRLower", true); + checkVectorSize(dscbAlphaRUpper, "dscbAlphaRUpper", true); + checkVectorSize(dscbNLInitial, "dscbNLInitial", true); + checkVectorSize(dscbNLLower, "dscbNLLower", true); + checkVectorSize(dscbNLUpper, "dscbNLUpper", true); + checkVectorSize(dscbNRInitial, "dscbNRInitial", true); + checkVectorSize(dscbNRLower, "dscbNRLower", true); + checkVectorSize(dscbNRUpper, "dscbNRUpper", true); for (int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) { sliceVarLimits[iSliceVar] = sliceVarMin[iSliceVar]; @@ -268,6 +321,10 @@ void runMassFitter(const std::string& configFileName) auto* hRawYieldsSigma = new TH1D("hRawYieldsSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data()); auto* hRawYieldsSecSigma = new TH1D("hRawYieldsSecSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data()); auto* hRawYieldsFracDoubleGaus = new TH1D("hRawYieldsFracDoubleGaus", ";" + sliceVarName + "(" + sliceVarUnit + ");fraction of double gaussian", nSliceVarBins, sliceVarLimits.data()); + auto* hRawYieldsDscbAlphaL = new TH1D("hRawYieldsDscbAlphaL", ";" + sliceVarName + "(" + sliceVarUnit + ");#alpha_{L}", nSliceVarBins, sliceVarLimits.data()); + auto* hRawYieldsDscbAlphaR = new TH1D("hRawYieldsDscbAlphaR", ";" + sliceVarName + "(" + sliceVarUnit + ");#alpha_{R}", nSliceVarBins, sliceVarLimits.data()); + auto* hRawYieldsDscbNL = new TH1D("hRawYieldsDscbNL", ";" + sliceVarName + "(" + sliceVarUnit + ");n_{L}", nSliceVarBins, sliceVarLimits.data()); + auto* hRawYieldsDscbNR = new TH1D("hRawYieldsDscbNR", ";" + sliceVarName + "(" + sliceVarUnit + ");n_{R}", nSliceVarBins, sliceVarLimits.data()); enum { ConfigMassMin = 1, @@ -300,6 +357,10 @@ void runMassFitter(const std::string& configFileName) setHistoStyle(hRawYieldsSigma); setHistoStyle(hRawYieldsSecSigma); setHistoStyle(hRawYieldsFracDoubleGaus); + setHistoStyle(hRawYieldsDscbAlphaL); + setHistoStyle(hRawYieldsDscbAlphaR); + setHistoStyle(hRawYieldsDscbNL); + setHistoStyle(hRawYieldsDscbNR); auto getHistToFix = [&nSliceVarBins](bool const& isFix, std::vector const& fixManual, std::string const& fixFileName, std::string const& var) -> TH1* { TH1* histToFix = nullptr; @@ -400,6 +461,10 @@ void runMassFitter(const std::string& configFileName) 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"); if (!isMc && enableRefl) { reflOverSgn = hMassSgn[iSliceVar]->Integral(hMassSgn[iSliceVar]->FindBin(massMin[iSliceVar] * 1.0001), hMassSgn[iSliceVar]->FindBin(massMax[iSliceVar] * 0.999)); @@ -408,6 +473,24 @@ void runMassFitter(const std::string& configFileName) massFitter->setTemplateReflections(hMassRefl[iSliceVar]); } + auto setDscbParameter = [&](const std::vector& vec, void (HFInvMassFitter::*setter)(double)) { + if (static_cast(vec.size()) == nSliceVarBins) { + (massFitter->*setter)(vec[iSliceVar]); + } + }; + setDscbParameter(dscbAlphaLInitial, &HFInvMassFitter::setDscbAlphaLInitialValue); + setDscbParameter(dscbAlphaLLower, &HFInvMassFitter::setDscbAlphaLLowLimit); + setDscbParameter(dscbAlphaLUpper, &HFInvMassFitter::setDscbAlphaLUpLimit); + setDscbParameter(dscbAlphaRInitial, &HFInvMassFitter::setDscbAlphaRInitialValue); + setDscbParameter(dscbAlphaRLower, &HFInvMassFitter::setDscbAlphaRLowLimit); + setDscbParameter(dscbAlphaRUpper, &HFInvMassFitter::setDscbAlphaRUpLimit); + setDscbParameter(dscbNLInitial, &HFInvMassFitter::setDscbNLInitialValue); + setDscbParameter(dscbNLLower, &HFInvMassFitter::setDscbNLLowLimit); + setDscbParameter(dscbNLUpper, &HFInvMassFitter::setDscbNLUpLimit); + setDscbParameter(dscbNRInitial, &HFInvMassFitter::setDscbNRInitialValue); + setDscbParameter(dscbNRLower, &HFInvMassFitter::setDscbNRLowLimit); + setDscbParameter(dscbNRUpper, &HFInvMassFitter::setDscbNRUpLimit); + massFitter->doFit(); auto drawOnCanvas = [&](std::vector& canvas, std::function drawer) { @@ -444,6 +527,14 @@ void runMassFitter(const std::string& configFileName) const double meanErr = massFitter->getMeanUncertainty(); const double sigma = massFitter->getSigma(); const double sigmaErr = massFitter->getSigmaUncertainty(); + const double dscbAlphaL = massFitter->getDscbAlphaL(); + const double dscbAlphaR = massFitter->getDscbAlphaR(); + const double dscbNL = massFitter->getDscbNL(); + const double dscbNR = massFitter->getDscbNR(); + const double dscbAlphaLErr = massFitter->getDscbAlphaLUncertainty(); + const double dscbAlphaRErr = massFitter->getDscbAlphaRUncertainty(); + const double dscbNErrL = massFitter->getDscbNLUncertainty(); + const double dscbNErrR = massFitter->getDscbNRUncertainty(); hRawYieldsSignal->SetBinContent(iSliceVar + 1, rawYield); hRawYieldsSignal->SetBinError(iSliceVar + 1, rawYieldErr); @@ -464,6 +555,14 @@ void runMassFitter(const std::string& configFileName) hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma); hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr); hReflectionOverSignal->SetBinContent(iSliceVar + 1, reflOverSgn); + hRawYieldsDscbAlphaL->SetBinContent(iSliceVar + 1, dscbAlphaL); + hRawYieldsDscbAlphaL->SetBinError(iSliceVar + 1, dscbAlphaLErr); + hRawYieldsDscbAlphaR->SetBinContent(iSliceVar + 1, dscbAlphaR); + hRawYieldsDscbAlphaR->SetBinError(iSliceVar + 1, dscbAlphaRErr); + hRawYieldsDscbNL->SetBinContent(iSliceVar + 1, dscbNL); + hRawYieldsDscbNL->SetBinError(iSliceVar + 1, dscbNErrL); + hRawYieldsDscbNR->SetBinContent(iSliceVar + 1, dscbNR); + hRawYieldsDscbNR->SetBinError(iSliceVar + 1, dscbNErrR); if (sgnFunc[iSliceVar] != HFInvMassFitter::SingleGaus) { // TODO foresee DSCB and Voigt cases const double secSigma = massFitter->getSecSigma(); @@ -519,6 +618,12 @@ void runMassFitter(const std::string& configFileName) if (enableRefl) { hReflectionOverSignal->Write(); } + if (std::find(sgnFunc.begin(), sgnFunc.end(), HFInvMassFitter::DoubleSidedCrystalBall) != sgnFunc.end()) { + hRawYieldsDscbAlphaL->Write(); + hRawYieldsDscbAlphaR->Write(); + hRawYieldsDscbNL->Write(); + hRawYieldsDscbNR->Write(); + } hFitConfig->Write(); outputFile.Close(); @@ -631,6 +736,24 @@ void readJsonVector(std::vector& vec, const Document& config, const std::stri } } +void readJsonVectorFromHisto(std::vector& vec, const Document& config, const std::string& fileNameFieldName, const std::string& histoNameFieldName) +{ + if (!vec.empty()) { + throw std::runtime_error("readJsonVectorFromHisto(): vector is not empty!"); + } + const auto fileName = readJsonField(config, fileNameFieldName); + const auto histoName = readJsonField(config, histoNameFieldName); + if (fileName.empty() || histoName.empty()) { + return; + } + TFile* inputFile = openFileWithNullptrCheck(fileName); + TH1* histo = getObjectWithNullPtrCheck(inputFile, histoName); + for (int iBin = 1; iBin <= histo->GetNbinsX(); iBin++) { + vec.push_back(histo->GetBinContent(iBin)); + } + inputFile->Close(); +} + int main(int argc, const char* argv[]) { if (argc == 1) { From 5bf89d4f4843b2ab46dad6e82b1b4de0fab16d5b Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Sun, 1 Feb 2026 15:40:35 +0100 Subject: [PATCH 07/26] bugfix read vector from histo and init rndSeed --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 5 +++-- PWGHF/D2H/Macros/HFInvMassFitter.h | 3 +-- PWGHF/D2H/Macros/runMassFitter.C | 7 +++---- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 5cd7151d3e4..9333856288b 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -62,7 +62,8 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, double minValue, double maxValue, int fitTypeBkg, - int fitTypeSgn) : mHistoInvMass(nullptr), + int fitTypeSgn, + int randomSeed) : mHistoInvMass(nullptr), mFitOption("L,E"), mMinMass(minValue), mMaxMass(maxValue), @@ -149,7 +150,7 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mHistoTemplateRefl(nullptr), mDrawBgPrefit(false), mHighlightPeakRegion(false), - mRandomSeed(-1), + mRandomSeed(randomSeed), mRandomGen(nullptr) { // standard constructor diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 5e7133e08e2..6902bbc7d04 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -77,7 +77,7 @@ class HFInvMassFitter : public TNamed }; std::array namesOfReflPdf{"reflFuncGaus", "reflFuncDoubleGaus", "reflFuncPoly3", "reflFuncPoly6"}; HFInvMassFitter() = delete; - HFInvMassFitter(TH1* histoToFit, double minValue, double maxValue, int fitTypeBkg = Expo, int fitTypeSgn = SingleGaus); + HFInvMassFitter(TH1* histoToFit, double minValue, double maxValue, int fitTypeBkg = Expo, int fitTypeSgn = SingleGaus, int randomSeed = -1); ~HFInvMassFitter() override; void setHistogramForFit(TH1* histoToFit); void setUseLikelihoodFit() { mFitOption = "L,E"; } @@ -169,7 +169,6 @@ class HFInvMassFitter : public TNamed void drawResidual(TVirtualPad* c); void drawRatio(TVirtualPad* c); void drawReflection(TVirtualPad* c); - void setRandomSeed(int seed) { mRandomSeed = seed; } private: HFInvMassFitter(const HFInvMassFitter& source); diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index e7ab0258307..40d090f1aa0 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -431,7 +431,7 @@ void runMassFitter(const std::string& configFileName) double reflOverSgn = 0; - HFInvMassFitter* massFitter = new HFInvMassFitter(hMass[iSliceVar], massMin[iSliceVar], massMax[iSliceVar], bkgFunc[iSliceVar], sgnFunc[iSliceVar]); + HFInvMassFitter* massFitter = new HFInvMassFitter(hMass[iSliceVar], massMin[iSliceVar], massMax[iSliceVar], bkgFunc[iSliceVar], sgnFunc[iSliceVar], randomSeed); massFitter->setDrawBgPrefit(drawBgPrefit); massFitter->setHighlightPeakRegion(highlightPeakRegion); massFitter->setInitialGaussianMean(massPDG); @@ -442,7 +442,6 @@ void runMassFitter(const std::string& configFileName) } else { massFitter->setUseChi2Fit(); } - massFitter->setRandomSeed(randomSeed); auto setFixedValue = [&iSliceVar](bool const& isFix, std::vector const& fixManual, const TH1* histToFix, std::function setFunc, std::string const& var) -> void { if (isFix) { @@ -741,8 +740,8 @@ void readJsonVectorFromHisto(std::vector& vec, const Document& config, c if (!vec.empty()) { throw std::runtime_error("readJsonVectorFromHisto(): vector is not empty!"); } - const auto fileName = readJsonField(config, fileNameFieldName); - const auto histoName = readJsonField(config, histoNameFieldName); + const auto fileName = readJsonField(config, fileNameFieldName, ""); + const auto histoName = readJsonField(config, histoNameFieldName, ""); if (fileName.empty() || histoName.empty()) { return; } From b8cf48eddb40ab462d37958a62fa859bd02984bf Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Sun, 1 Feb 2026 16:02:03 +0100 Subject: [PATCH 08/26] add DSCB to fillWorkspace() --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 9333856288b..921e682c2ce 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -818,6 +818,15 @@ RooAbsPdf* HFInvMassFitter::createSignalFitFunction(RooWorkspace* workspace) mRooSecSigmaSgn = workspace->var("sigmaSec"); mRooMeanSgn = workspace->var("meanSec"); } break; + case DoubleSidedCrystalBall: { + sgnPdf = workspace->pdf("sgnFuncDSCB"); + mRooSigmaSgn = workspace->var("sigma"); + mRooMeanSgn = workspace->var("mean"); + mRooDscbAlphaL = workspace->var("alphaL"); + mRooDscbNL = workspace->var("nL"); + mRooDscbAlphaR = workspace->var("alphaR"); + mRooDscbNR = workspace->var("nR"); + } break; default: break; } From ab455e9b968d20349e86b633577da8335b78d83d Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Sun, 1 Feb 2026 16:33:53 +0100 Subject: [PATCH 09/26] DSCB: second sigma (concession) --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 921e682c2ce..15cabb9f6f8 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -821,6 +821,7 @@ RooAbsPdf* HFInvMassFitter::createSignalFitFunction(RooWorkspace* workspace) case DoubleSidedCrystalBall: { sgnPdf = workspace->pdf("sgnFuncDSCB"); mRooSigmaSgn = workspace->var("sigma"); + mRooSecSigmaSgn = workspace->var("sigma"); mRooMeanSgn = workspace->var("mean"); mRooDscbAlphaL = workspace->var("alphaL"); mRooDscbNL = workspace->var("nL"); From 07957b2cbea642697a15adc89181552111297f2c Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Sun, 1 Feb 2026 17:49:51 +0100 Subject: [PATCH 10/26] MC: allow single signal histogram instead of prompt+fd --- PWGHF/D2H/Macros/runMassFitter.C | 43 +++++++++++++++++++++++++++----- 1 file changed, 37 insertions(+), 6 deletions(-) diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 40d090f1aa0..741a80601a7 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -97,9 +97,11 @@ void runMassFitter(const std::string& configFileName) std::vector inputHistoName; std::vector promptHistoName; std::vector fdHistoName; + std::vector signalHistoName; std::vector reflHistoName; std::vector promptSecPeakHistoName; std::vector fdSecPeakHistoName; + std::vector signalSecPeakHistoName; std::vector sliceVarMin; std::vector sliceVarMax; std::vector massMin; @@ -127,9 +129,11 @@ void runMassFitter(const std::string& configFileName) readJsonVector(inputHistoName, config, "InputHistoName", true); readJsonVector(promptHistoName, config, "PromptHistoName"); readJsonVector(fdHistoName, config, "FDHistoName"); + readJsonVector(signalHistoName, config, "SignalHistoName"); readJsonVector(reflHistoName, config, "ReflHistoName"); readJsonVector(promptSecPeakHistoName, config, "PromptSecPeakHistoName"); readJsonVector(fdSecPeakHistoName, config, "FDSecPeakHistoName"); + readJsonVector(signalSecPeakHistoName, config, "SignalSecPeakHistoName"); const bool fixMean = readJsonField(config, "FixMean", false); const std::string meanFile = readJsonField(config, "MeanFile", ""); @@ -211,9 +215,11 @@ void runMassFitter(const std::string& configFileName) checkVectorSize(inputHistoName, "inputHistoName", true); checkVectorSize(promptHistoName, "promptHistoName", true); checkVectorSize(fdHistoName, "fdHistoName", true); + checkVectorSize(signalHistoName, "signalHistoName", true); checkVectorSize(reflHistoName, "reflHistoName", true); checkVectorSize(promptSecPeakHistoName, "promptSecPeakHistoName", true); checkVectorSize(fdSecPeakHistoName, "fdSecPeakHistoName", true); + checkVectorSize(signalSecPeakHistoName, "signalSecPeakHistoName", true); checkVectorSize(sliceVarMin, "sliceVarMin"); checkVectorSize(sliceVarMax, "sliceVarMax"); checkVectorSize(massMin, "massMin"); @@ -238,6 +244,15 @@ void runMassFitter(const std::string& configFileName) checkVectorSize(dscbNRLower, "dscbNRLower", true); checkVectorSize(dscbNRUpper, "dscbNRUpper", true); + auto checkVectorSizeMcHistograms = [](const auto& vecSignal, const auto& vecPrompt, const auto& vecFd) { + const auto signalSize = vecSignal.size(); + const auto promptSize = vecPrompt.size(); + const auto fdSize = vecFd.size(); + if (!((signalSize > 0 && promptSize == 0 && fdSize == 0) || (signalSize == 0 && promptSize > 0 && fdSize > 0))) { + throw std::runtime_error("ERROR: either signal histogram must be provided or both prompt and fd, but not all three. Exit"); + } + }; + for (int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) { sliceVarLimits[iSliceVar] = sliceVarMin[iSliceVar]; @@ -290,15 +305,31 @@ void runMassFitter(const std::string& configFileName) hMass[iSliceVar] = getObjectWithNullPtrCheck(inputFile, inputHistoName[iSliceVar]); if (enableRefl) { hMassRefl[iSliceVar] = getObjectWithNullPtrCheck(inputFileRefl, reflHistoName[iSliceVar]); - hMassSgn[iSliceVar] = getObjectWithNullPtrCheck(inputFileRefl, fdHistoName[iSliceVar]); - hMassSgn[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFileRefl, promptHistoName[iSliceVar])); + + checkVectorSizeMcHistograms(signalHistoName, promptHistoName, fdHistoName); + if (!signalHistoName.empty()) { + hMassSgn[iSliceVar] = getObjectWithNullPtrCheck(inputFileRefl, signalHistoName[iSliceVar]); + } else { + hMassSgn[iSliceVar] = getObjectWithNullPtrCheck(inputFileRefl, fdHistoName[iSliceVar]); + hMassSgn[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFileRefl, promptHistoName[iSliceVar])); + } } } else { - hMass[iSliceVar] = getObjectWithNullPtrCheck(inputFile, promptHistoName[iSliceVar]); - hMass[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFile, fdHistoName[iSliceVar])); + checkVectorSizeMcHistograms(signalHistoName, promptHistoName, fdHistoName); + if (!signalHistoName.empty()) { + hMass[iSliceVar] = getObjectWithNullPtrCheck(inputFile, signalHistoName[iSliceVar]); + } else { + hMass[iSliceVar] = getObjectWithNullPtrCheck(inputFile, promptHistoName[iSliceVar]); + hMass[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFile, fdHistoName[iSliceVar])); + } if (includeSecPeak) { - hMass[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFile, promptSecPeakHistoName[iSliceVar])); - hMass[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFile, fdSecPeakHistoName[iSliceVar])); + checkVectorSizeMcHistograms(signalSecPeakHistoName, promptSecPeakHistoName, fdSecPeakHistoName); + if (!signalHistoName.empty()) { + hMass[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFile, signalSecPeakHistoName[iSliceVar])); + } else { + hMass[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFile, promptSecPeakHistoName[iSliceVar])); + hMass[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFile, fdSecPeakHistoName[iSliceVar])); + } } } hMass[iSliceVar]->SetDirectory(nullptr); From b151616fb590fcfe1c21828037312600ea83f8ed Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Tue, 3 Feb 2026 19:30:08 +0100 Subject: [PATCH 11/26] read vector from json: allow for a single value if all are equal --- PWGHF/D2H/Macros/config_massfitter.json | 15 ++- PWGHF/D2H/Macros/runMassFitter.C | 126 +++++++++++++++--------- 2 files changed, 88 insertions(+), 53 deletions(-) diff --git a/PWGHF/D2H/Macros/config_massfitter.json b/PWGHF/D2H/Macros/config_massfitter.json index 7ee52cace9a..cb84f8aaafa 100644 --- a/PWGHF/D2H/Macros/config_massfitter.json +++ b/PWGHF/D2H/Macros/config_massfitter.json @@ -101,21 +101,24 @@ 0.9, 1.6 ], - "MassMin": [ + "MassMin": 2.12, + "_MassMin": [ 2.12, 2.12, 2.12, 2.12, 2.12 ], - "MassMax": [ + "MassMax": 2.42, + "_MassMax": [ 2.42, 2.42, 2.42, 2.42, 2.42 ], - "Rebin": [ + "Rebin": 4, + "_Rebin": [ 4, 4, 4, @@ -128,7 +131,8 @@ "true: likelihood fit", "false: chi2 fit" ], - "BkgFunc": [ + "BkgFunc": 2, + "_BkgFunc": [ 2, 2, 2, @@ -144,7 +148,8 @@ "5 for Poly3", "6 for NoBkg" ], - "SgnFunc": [ + "SgnFunc": 0, + "_SgnFunc": [ 0, 0, 0, diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 741a80601a7..0da6385ceca 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -50,6 +50,8 @@ 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"); @@ -66,9 +68,12 @@ T readJsonField(const Document& config, const std::string& fieldName); template void readJsonVector(std::vector& vec, const Document& config, const std::string& fieldName, bool isRequired = false); +template +void readJsonVectorFlexible(std::vector& vec, const Document& config, int nHistograms, const std::string& fieldName, bool isRequired = false); + void readJsonVectorFromHisto(std::vector& vec, const Document& config, const std::string& fileNameFieldName, const std::string& histoNameFieldName); -void divideCanvas(TCanvas* c, int nSliceVarBins); +void divideCanvas(TCanvas* c, int nHistograms); void setHistoStyle(TH1* histo, Color_t color = kBlack, Size_t markerSize = 1); @@ -127,6 +132,8 @@ void runMassFitter(const std::string& configFileName) std::vector dscbNRUpper; readJsonVector(inputHistoName, config, "InputHistoName", true); + const int nHistograms = static_cast(inputHistoName.size()); + readJsonVector(promptHistoName, config, "PromptHistoName"); readJsonVector(fdHistoName, config, "FDHistoName"); readJsonVector(signalHistoName, config, "SignalHistoName"); @@ -137,19 +144,19 @@ void runMassFitter(const std::string& configFileName) const bool fixMean = readJsonField(config, "FixMean", false); const std::string meanFile = readJsonField(config, "MeanFile", ""); - readJsonVector(fixMeanManual, config, "FixMeanManual"); + readJsonVectorFlexible(fixMeanManual, config, nHistograms, "FixMeanManual"); const bool fixSigma = readJsonField(config, "FixSigma", false); const std::string sigmaFile = readJsonField(config, "SigmaFile", ""); - readJsonVector(fixSigmaManual, config, "FixSigmaManual"); + readJsonVectorFlexible(fixSigmaManual, config, nHistograms, "FixSigmaManual"); const bool fixSecondSigma = readJsonField(config, "FixSecondSigma", false); const std::string secondSigmaFile = readJsonField(config, "SecondSigmaFile", ""); - readJsonVector(fixSecondSigmaManual, config, "FixSecondSigmaManual"); + readJsonVectorFlexible(fixSecondSigmaManual, config, nHistograms, "FixSecondSigmaManual"); const bool fixFracDoubleGaus = readJsonField(config, "FixFracDoubleGaus", false); const std::string fracDoubleGausFile = readJsonField(config, "FracDoubleGausFile", ""); - readJsonVector(fixFracDoubleGausManual, config, "FixFracDoubleGausManual"); + readJsonVectorFlexible(fixFracDoubleGausManual, config, nHistograms, "FixFracDoubleGausManual"); const bool fixDscbTailParams = readJsonField(config, "FixDscbTailParams", false); @@ -159,16 +166,16 @@ void runMassFitter(const std::string& configFileName) readJsonVector(sliceVarMin, config, "SliceVarMin", true); readJsonVector(sliceVarMax, config, "SliceVarMax", true); - readJsonVector(massMin, config, "MassMin", true); - readJsonVector(massMax, config, "MassMax", true); + readJsonVectorFlexible(massMin, config, nHistograms, "MassMin", true); + readJsonVectorFlexible(massMax, config, nHistograms, "MassMax", true); - readJsonVector(nRebin, config, "Rebin", true); + readJsonVectorFlexible(nRebin, config, nHistograms, "Rebin", true); bool const includeSecPeak = readJsonField(config, "InclSecPeak", false); bool const useLikelihood = readJsonField(config, "UseLikelihood"); - readJsonVector(bkgFunc, config, "BkgFunc", true); - readJsonVector(sgnFunc, config, "SgnFunc", true); + 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); @@ -200,11 +207,10 @@ void runMassFitter(const std::string& configFileName) readJsonVectorFromHisto(dscbNRLower, config, "DscbParametersFile", "DscbNRLowerHisto"); readJsonVectorFromHisto(dscbNRUpper, config, "DscbParametersFile", "DscbNRUpperHisto"); - const int nSliceVarBins = static_cast(sliceVarMin.size()); - std::vector sliceVarLimits(nSliceVarBins + 1); + std::vector sliceVarLimits(nHistograms + 1); auto checkVectorSize = [&](const auto& vec, const std::string& name = "", const bool isEmptyOk = false) { - if (vec.size() != static_cast(nSliceVarBins)) { + if (vec.size() != static_cast(nHistograms)) { if (isEmptyOk && vec.empty()) { return; } @@ -253,7 +259,7 @@ void runMassFitter(const std::string& configFileName) } }; - for (int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) { + for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { sliceVarLimits[iSliceVar] = sliceVarMin[iSliceVar]; if (bkgFunc[iSliceVar] < 0 || bkgFunc[iSliceVar] >= HFInvMassFitter::NTypesOfBkgPdf) { @@ -266,7 +272,7 @@ void runMassFitter(const std::string& configFileName) throw std::runtime_error("ERROR: only SingleGaus, DoubleGaus and DoubleGausSigmaRatioPar signal functions supported! Exit"); } } - sliceVarLimits[nSliceVarBins] = sliceVarMax[nSliceVarBins - 1]; + sliceVarLimits[nHistograms] = sliceVarMax[nHistograms - 1]; struct DecayInfo { std::string decayProducts; @@ -296,11 +302,11 @@ void runMassFitter(const std::string& configFileName) TFile* inputFileRefl = enableRefl ? openFileWithNullptrCheck(reflFileName) : nullptr; - std::vector hMassSgn(nSliceVarBins); - std::vector hMassRefl(nSliceVarBins); - std::vector hMass(nSliceVarBins); + std::vector hMassSgn(nHistograms); + std::vector hMassRefl(nHistograms); + std::vector hMass(nHistograms); - for (int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) { + for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { if (!isMc) { hMass[iSliceVar] = getObjectWithNullPtrCheck(inputFile, inputHistoName[iSliceVar]); if (enableRefl) { @@ -340,22 +346,22 @@ void runMassFitter(const std::string& configFileName) } // define output histos - auto* hRawYieldsSignal = new TH1D("hRawYieldsSignal", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsSignalCounted = new TH1D("hRawYieldsSignalCounted", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield via bin count", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsBkg = new TH1D("hRawYieldsBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");Background (3#sigma)", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsSgnOverBkg = new TH1D("hRawYieldsSgnOverBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");S/B (3#sigma)", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsSignificance = new TH1D("hRawYieldsSignificance", ";" + sliceVarName + "(" + sliceVarUnit + ");significance (3#sigma)", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsChiSquareBkg = new TH1D("hRawYieldsChiSquareBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");#chi^{2}/#it{ndf}", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsChiSquareTotal = new TH1D("hRawYieldsChiSquareTotal", ";" + sliceVarName + "(" + sliceVarUnit + ");#chi^{2}/#it{ndf}", nSliceVarBins, sliceVarLimits.data()); - auto* hReflectionOverSignal = new TH1D("hReflectionOverSignal", ";" + sliceVarName + "(" + sliceVarUnit + ");Refl/Signal", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsMean = new TH1D("hRawYieldsMean", ";" + sliceVarName + "(" + sliceVarUnit + ");mean (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsSigma = new TH1D("hRawYieldsSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsSecSigma = new TH1D("hRawYieldsSecSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsFracDoubleGaus = new TH1D("hRawYieldsFracDoubleGaus", ";" + sliceVarName + "(" + sliceVarUnit + ");fraction of double gaussian", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsDscbAlphaL = new TH1D("hRawYieldsDscbAlphaL", ";" + sliceVarName + "(" + sliceVarUnit + ");#alpha_{L}", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsDscbAlphaR = new TH1D("hRawYieldsDscbAlphaR", ";" + sliceVarName + "(" + sliceVarUnit + ");#alpha_{R}", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsDscbNL = new TH1D("hRawYieldsDscbNL", ";" + sliceVarName + "(" + sliceVarUnit + ");n_{L}", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsDscbNR = new TH1D("hRawYieldsDscbNR", ";" + sliceVarName + "(" + sliceVarUnit + ");n_{R}", nSliceVarBins, 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, @@ -366,7 +372,7 @@ void runMassFitter(const std::string& configFileName) ConfigSgnFunc, NConfigsToSave }; - auto* hFitConfig = new TH2F("hfitConfig", "Fit Configurations", NConfigsToSave - 1, 0, NConfigsToSave - 1, nSliceVarBins, 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"}; hFitConfig->SetStats(false); for (int i = 0; i < NConfigsToSave - 1; i++) { @@ -393,7 +399,7 @@ void runMassFitter(const std::string& configFileName) setHistoStyle(hRawYieldsDscbNL); setHistoStyle(hRawYieldsDscbNR); - auto getHistToFix = [&nSliceVarBins](bool const& isFix, std::vector const& fixManual, std::string const& fixFileName, std::string const& var) -> TH1* { + auto getHistToFix = [&nHistograms](bool const& isFix, std::vector const& fixManual, std::string const& fixFileName, std::string const& var) -> TH1* { TH1* histToFix = nullptr; if (isFix) { if (fixManual.empty()) { @@ -401,7 +407,7 @@ void runMassFitter(const std::string& configFileName) const std::string histName = "hRawYields" + var; histToFix = getObjectWithNullPtrCheck(fixInputFile, histName); histToFix->SetDirectory(nullptr); - if (histToFix->GetNbinsX() != nSliceVarBins) { + if (histToFix->GetNbinsX() != nHistograms) { throw std::runtime_error("Different number of bins for this analysis and histo for fixed " + var); } fixInputFile->Close(); @@ -416,19 +422,19 @@ void runMassFitter(const std::string& configFileName) TH1* hFracDoubleGausToFix = getHistToFix(fixFracDoubleGaus, fixFracDoubleGausManual, fracDoubleGausFile, "FracDoubleGaus"); int canvasSize[2] = {1920, 1080}; - if (nSliceVarBins == 1) { + if (nHistograms == 1) { canvasSize[0] = 500; canvasSize[1] = 500; } int constexpr NCanvasesMax = 20; // do not put more than 20 bins per canvas to make them visible - const int nCanvases = std::ceil(static_cast(nSliceVarBins) / NCanvasesMax); + const int nCanvases = std::ceil(static_cast(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) ? nSliceVarBins : NCanvasesMax; + 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); @@ -444,7 +450,7 @@ void runMassFitter(const std::string& configFileName) } } - for (int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) { + for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { const int iCanvas = std::floor(static_cast(iSliceVar) / NCanvasesMax); hMass[iSliceVar]->Rebin(nRebin[iSliceVar]); @@ -504,7 +510,7 @@ void runMassFitter(const std::string& configFileName) } auto setDscbParameter = [&](const std::vector& vec, void (HFInvMassFitter::*setter)(double)) { - if (static_cast(vec.size()) == nSliceVarBins) { + if (static_cast(vec.size()) == nHistograms) { (massFitter->*setter)(vec[iSliceVar]); } }; @@ -524,7 +530,7 @@ void runMassFitter(const std::string& configFileName) massFitter->doFit(); auto drawOnCanvas = [&](std::vector& canvas, std::function drawer) { - if (nSliceVarBins > 1) { + if (nHistograms > 1) { canvas[iCanvas]->cd(iSliceVar - NCanvasesMax * iCanvas + 1); } else { canvas[iCanvas]->cd(); @@ -631,7 +637,7 @@ void runMassFitter(const std::string& configFileName) } } - for (int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) { + for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { hMass[iSliceVar]->Write(); } hRawYieldsSignal->Write(); @@ -685,10 +691,10 @@ void setHistoStyle(TH1* histo, Color_t color, Size_t markerSize) histo->SetLineColor(color); } -void divideCanvas(TCanvas* canvas, int nSliceVarBins) +void divideCanvas(TCanvas* canvas, int nHistograms) { - int nCols = std::ceil(std::sqrt(nSliceVarBins)); - int nRows = std::ceil(static_cast(nSliceVarBins) / nCols); + int nCols = std::ceil(std::sqrt(nHistograms)); + int nRows = std::ceil(static_cast(nHistograms) / nCols); canvas->Divide(nCols, nRows); } @@ -766,6 +772,30 @@ void readJsonVector(std::vector& vec, const Document& config, const std::stri } } +template +void readJsonVectorFlexible(std::vector& vec, const Document& config, int nHistograms, const std::string& fieldName, bool isRequired) +{ + if constexpr (!(std::is_same_v, int> || std::is_same_v, double>)) { + static_assert(AlwaysFalse, "readJsonVectorFlexible(): unsupported type!"); + } + if (!vec.empty()) { + throw std::runtime_error("readJsonVectorFlexible(): vector is not empty!"); + } + if (!config.HasMember(fieldName.c_str())) { + if (isRequired) { + throw std::runtime_error("readJsonVectorFlexible(): missing required field " + fieldName); + } else { + return; + } + } + if (config[fieldName.c_str()].IsArray()) { + readJsonVector(vec, config, fieldName); + } else { + const T value = readJsonField(config, fieldName); + vec.assign(nHistograms, value); + } +} + void readJsonVectorFromHisto(std::vector& vec, const Document& config, const std::string& fileNameFieldName, const std::string& histoNameFieldName) { if (!vec.empty()) { From 68b3aeb48aea8fcc28b3f531b030495b9cdb7586 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Mon, 9 Feb 2026 11:11:47 +0100 Subject: [PATCH 12/26] fix 053cdd86: perform checkVectorSizeMcHistograms() only once --- PWGHF/D2H/Macros/runMassFitter.C | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 0da6385ceca..7ac7a468ea0 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -259,6 +259,13 @@ void runMassFitter(const std::string& configFileName) } }; + if ((!isMc && enableRefl) || isMc) { + checkVectorSizeMcHistograms(signalHistoName, promptHistoName, fdHistoName); + } + if (isMc && includeSecPeak) { + checkVectorSizeMcHistograms(signalSecPeakHistoName, promptSecPeakHistoName, fdSecPeakHistoName); + } + for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { sliceVarLimits[iSliceVar] = sliceVarMin[iSliceVar]; @@ -311,8 +318,6 @@ void runMassFitter(const std::string& configFileName) hMass[iSliceVar] = getObjectWithNullPtrCheck(inputFile, inputHistoName[iSliceVar]); if (enableRefl) { hMassRefl[iSliceVar] = getObjectWithNullPtrCheck(inputFileRefl, reflHistoName[iSliceVar]); - - checkVectorSizeMcHistograms(signalHistoName, promptHistoName, fdHistoName); if (!signalHistoName.empty()) { hMassSgn[iSliceVar] = getObjectWithNullPtrCheck(inputFileRefl, signalHistoName[iSliceVar]); } else { @@ -321,7 +326,6 @@ void runMassFitter(const std::string& configFileName) } } } else { - checkVectorSizeMcHistograms(signalHistoName, promptHistoName, fdHistoName); if (!signalHistoName.empty()) { hMass[iSliceVar] = getObjectWithNullPtrCheck(inputFile, signalHistoName[iSliceVar]); } else { @@ -329,7 +333,6 @@ void runMassFitter(const std::string& configFileName) hMass[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFile, fdHistoName[iSliceVar])); } if (includeSecPeak) { - checkVectorSizeMcHistograms(signalSecPeakHistoName, promptSecPeakHistoName, fdSecPeakHistoName); if (!signalHistoName.empty()) { hMass[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFile, signalSecPeakHistoName[iSliceVar])); } else { From 914032dc09cbac981b0fba2963ece1f5320e6973 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Mon, 9 Feb 2026 11:37:30 +0100 Subject: [PATCH 13/26] foresee empty InputHistoName when IsMc==true --- 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 7ac7a468ea0..801854782e9 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -37,6 +37,7 @@ #include #include +#include #include #include #include @@ -131,12 +132,13 @@ void runMassFitter(const std::string& configFileName) std::vector dscbNRLower; std::vector dscbNRUpper; - readJsonVector(inputHistoName, config, "InputHistoName", true); - const int nHistograms = static_cast(inputHistoName.size()); - + readJsonVector(inputHistoName, config, "InputHistoName"); readJsonVector(promptHistoName, config, "PromptHistoName"); readJsonVector(fdHistoName, config, "FDHistoName"); readJsonVector(signalHistoName, config, "SignalHistoName"); + const std::array possibleInputHistogramSizes{inputHistoName.size(), promptHistoName.size(), fdHistoName.size(), signalHistoName.size()}; + const int nHistograms = static_cast(*std::max_element(possibleInputHistogramSizes.begin(), possibleInputHistogramSizes.end())); + readJsonVector(reflHistoName, config, "ReflHistoName"); readJsonVector(promptSecPeakHistoName, config, "PromptSecPeakHistoName"); readJsonVector(fdSecPeakHistoName, config, "FDSecPeakHistoName"); From a1bc66ff7cceac5f769ba6a454a230ec40beb909 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Mon, 9 Feb 2026 17:44:51 +0100 Subject: [PATCH 14/26] operate with S&B ranges properly in case of DSCB --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 33 ++++++++++++++++++---------- PWGHF/D2H/Macros/HFInvMassFitter.h | 4 +++- 2 files changed, 25 insertions(+), 12 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 15cabb9f6f8..1e8d03a8f00 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -54,6 +54,7 @@ #include #include #include +#include #include using namespace RooFit; @@ -202,14 +203,14 @@ void HFInvMassFitter::doFit() 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 - mNSigmaForSidebands * mSecSigma, mSecMass + mNSigmaForSidebands * mSecSigma); + 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); } } - mass->setRange("bkg", mMass - 4 * mSigmaSgn, mMass + 4 * mSigmaSgn); + mass->setRange("bkg", mMass - mNSigmaForSgn * mSigmaSgn, mMass + mNSigmaForSgn * mSigmaSgn); mass->setRange("full", mMinMass, mMaxMass); 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 @@ -226,13 +227,14 @@ void HFInvMassFitter::doFit() 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("signal")); + mTotalPdf->chi2FitTo(dataHistogram, Range("full")); } else { - mTotalPdf->fitTo(dataHistogram, Range("signal")); + mTotalPdf->fitTo(dataHistogram, Range("full")); } RooAbsReal* signalIntegralMc = mTotalPdf->createIntegral(*mass, NormSet(*mass), Range("signal")); // sig yield from fit mIntegralSgn = signalIntegralMc->getValV(); - calculateSignal(mRawYield, mRawYieldErr); // calculate signal and signal error + calculateSignal(mRawYield, mRawYieldErr); // calculate signal and signal error + countSignal(mRawYieldCounted, mRawYieldCountedErr); mTotalPdf->plotOn(mInvMassFrame, Name("Tot_c")); // plot total function // Fit to data ratio mRatioFrame = mass->frame(Title(Form("%s", mHistoInvMass->GetTitle()))); @@ -702,10 +704,7 @@ void HFInvMassFitter::drawReflection(TVirtualPad* pad) // calculate signal yield via bin counting void HFInvMassFitter::countSignal(double& signal, double& signalErr) const { - const double mean = mRooMeanSgn->getVal(); - const double sigma = mRooSecSigmaSgn->getVal(); - const double minForSgn = mean - mNSigmaForSidebands * sigma; - const double maxForSgn = mean + mNSigmaForSidebands * sigma; + const auto [minForSgn, maxForSgn] = getRangesOfSignal(); const int binForMinSgn = mHistoInvMass->FindBin(minForSgn); const int binForMaxSgn = mHistoInvMass->FindBin(maxForSgn); const double binForMinSgnUpperEdge = mHistoInvMass->GetBinLowEdge(binForMinSgn + 1); @@ -758,8 +757,7 @@ void HFInvMassFitter::calculateSignificance(double& significance, double& errSig // estimate Signal void HFInvMassFitter::checkForSignal(double& estimatedSignal) { - double const minForSgn = mMass - 4 * mSigmaSgn; - double const maxForSgn = mMass + 4 * mSigmaSgn; + auto const [minForSgn, maxForSgn] = getRangesOfSignal(); int const binForMinSgn = mHistoInvMass->FindBin(minForSgn); int const binForMaxSgn = mHistoInvMass->FindBin(maxForSgn); @@ -772,6 +770,19 @@ void HFInvMassFitter::checkForSignal(double& estimatedSignal) estimatedSignal = sum - bkg; } +// Estimate ranges where signal is located to be used in countSignal() and checkForSignal() +// It is mu +- n*sigma for Gaussian, and the entire mInv histogram for DoubleSidedCrystalBall (due to heavy tails) +std::pair HFInvMassFitter::getRangesOfSignal() const +{ + if (mTypeOfSgnPdf == DoubleSidedCrystalBall) { + return std::make_pair(mMinMass, mMaxMass); + } else { + const double mean = mRooMeanSgn->getVal(); + const double sigma = mRooSigmaSgn->getVal(); + return std::make_pair(mean - mNSigmaForSgn * sigma, mean + mNSigmaForSgn * sigma); + } +} + // Create Background Fit Function RooAbsPdf* HFInvMassFitter::createBackgroundFitFunction(RooWorkspace* workspace) const { diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 6902bbc7d04..1107d3d2702 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -36,6 +36,7 @@ #include #include +#include #include class HFInvMassFitter : public TNamed @@ -175,7 +176,8 @@ class HFInvMassFitter : public TNamed HFInvMassFitter& operator=(const HFInvMassFitter& source); void fillWorkspace(RooWorkspace& w) const; void highlightPeakRegion(const RooPlot* plot, Color_t color = kGray + 1, Width_t width = 1, Style_t style = 2) const; - double randomizeInitialParameter(const ParameterRanges& parameterRanges) const; + [[nodiscard]] double randomizeInitialParameter(const ParameterRanges& parameterRanges) const; + [[nodiscard]] std::pair getRangesOfSignal() const; TH1* mHistoInvMass; // histogram to fit std::string mFitOption; From 416684602368e30e801f1f9abc7b1a3a79f68f72 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Mon, 9 Mar 2026 11:48:00 +0100 Subject: [PATCH 15/26] use the same sigma for sgn and bg in 2-gaus --- 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 1e8d03a8f00..96ccf00298e 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -778,7 +778,7 @@ std::pair HFInvMassFitter::getRangesOfSignal() const return std::make_pair(mMinMass, mMaxMass); } else { const double mean = mRooMeanSgn->getVal(); - const double sigma = mRooSigmaSgn->getVal(); + const double sigma = mRooSecSigmaSgn->getVal(); return std::make_pair(mean - mNSigmaForSgn * sigma, mean + mNSigmaForSgn * sigma); } } From 82ded42bc3adcd52dc1eb407adafb5b92d1ac3b8 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 11 Mar 2026 10:15:45 +0100 Subject: [PATCH 16/26] countSignal() via residual histo --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 45 +++++++++++++++------------- PWGHF/D2H/Macros/HFInvMassFitter.h | 1 + 2 files changed, 26 insertions(+), 20 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 96ccf00298e..61cc3a0fce2 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -142,6 +142,7 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mReflFrame(nullptr), mReflOnlyFrame(nullptr), mResidualFrame(nullptr), + mResidualHist(nullptr), mRatioFrame(nullptr), mResidualFrameForCalculation(nullptr), mWorkspace(nullptr), @@ -315,9 +316,9 @@ void HFInvMassFitter::doFit() mChiSquareOverNdfTotal = mInvMassFrame->chiSquare("Tot_c", "data_c"); // calculate reduced chi2 / NDF // plot residual distribution - RooHist* residualHistogram = mInvMassFrame->residHist("data_c", "ReflBkg_c"); + mResidualHist = mInvMassFrame->residHist("data_c", "ReflBkg_c"); mResidualFrame = mass->frame(Title("Residual Distribution")); - mResidualFrame->addPlotable(residualHistogram, "p"); + mResidualFrame->addPlotable(mResidualHist, "p"); mSgnPdf->plotOn(mResidualFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineColor(kBlue)); } else { mTotalPdf = new RooAddPdf("mTotalPdf", "background + signal pdf", RooArgList(*bkgPdf, *sgnPdf), RooArgList(*mRooNBkg, *mRooNSgn)); @@ -333,8 +334,8 @@ void HFInvMassFitter::doFit() // plot residual distribution mResidualFrame = mass->frame(Title("Residual Distribution")); - RooHist* residualHistogram = mInvMassFrame->residHist("data_c", "Bkg_c"); - mResidualFrame->addPlotable(residualHistogram, "P"); + mResidualHist = mInvMassFrame->residHist("data_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()); @@ -704,26 +705,30 @@ void HFInvMassFitter::drawReflection(TVirtualPad* pad) // calculate signal yield via bin counting void HFInvMassFitter::countSignal(double& signal, double& signalErr) const { + const double binWidth = mResidualHist->GetX()[1] - mResidualHist->GetX()[0]; + const double firstBinLowEdge = mResidualHist->GetX()[0] - binWidth / 2; const auto [minForSgn, maxForSgn] = getRangesOfSignal(); - const int binForMinSgn = mHistoInvMass->FindBin(minForSgn); - const int binForMaxSgn = mHistoInvMass->FindBin(maxForSgn); - const double binForMinSgnUpperEdge = mHistoInvMass->GetBinLowEdge(binForMinSgn + 1); - const double binForMaxSgnLowerEdge = mHistoInvMass->GetBinLowEdge(binForMaxSgn); - const double binForMinSgnFraction = (binForMinSgnUpperEdge - minForSgn) / mHistoInvMass->GetBinWidth(binForMinSgn); - const double binForMaxSgnFraction = (maxForSgn - binForMaxSgnLowerEdge) / mHistoInvMass->GetBinWidth(binForMaxSgn); - - double sum = 0; - sum += mHistoInvMass->GetBinContent(binForMinSgn) * binForMinSgnFraction; + const int binForMinSgn = static_cast((minForSgn - firstBinLowEdge) / binWidth) + 1; + const int binForMaxSgn = static_cast((maxForSgn - firstBinLowEdge) / binWidth) + 1; + const double binForMinSgnUpperEdge = firstBinLowEdge + (binForMinSgn - 1) * binWidth + binWidth; + const double binForMaxSgnLowerEdge = firstBinLowEdge + (binForMaxSgn - 1) * binWidth; + const double binForMinSgnFraction = (binForMinSgnUpperEdge - minForSgn) / binWidth; + const double binForMaxSgnFraction = (maxForSgn - binForMaxSgnLowerEdge) / binWidth; + + auto square = [](double value) { return value * value; }; + + double sumValues{}, sumErrorsSquare{}; + sumValues += mResidualHist->GetY()[binForMinSgn - 1] * binForMinSgnFraction; + sumErrorsSquare += square(mResidualHist->GetErrorY(binForMinSgn - 1) * binForMinSgnFraction); for (int iBin = binForMinSgn + 1; iBin <= binForMaxSgn - 1; iBin++) { - sum += mHistoInvMass->GetBinContent(iBin); + sumValues += mResidualHist->GetY()[iBin - 1]; + sumErrorsSquare += square(mResidualHist->GetErrorY(iBin - 1)); } - sum += mHistoInvMass->GetBinContent(binForMaxSgn) * binForMaxSgnFraction; - - double bkg{}, errBkg{}; - calculateBackground(bkg, errBkg); + sumValues += mResidualHist->GetY()[binForMaxSgn - 1] * binForMaxSgnFraction; + sumErrorsSquare += square(mResidualHist->GetErrorY(binForMaxSgn - 1) * binForMaxSgnFraction); - signal = sum - bkg; - signalErr = std::sqrt(sum + errBkg * errBkg); // sum error squared is equal to sum + signal = sumValues; + signalErr = std::sqrt(sumErrorsSquare); } // calculate signal yield diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 1107d3d2702..5c387fa48c8 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -257,6 +257,7 @@ class HFInvMassFitter : public TNamed 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 RooPlot* mResidualFrameForCalculation; RooWorkspace* mWorkspace; /// workspace From b6355eeb7fc9a3f407642d4f7598a4c9f509e08f Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 11 Mar 2026 10:18:46 +0100 Subject: [PATCH 17/26] rm unused class member --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 1 - PWGHF/D2H/Macros/HFInvMassFitter.h | 19 +++++++++---------- 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 61cc3a0fce2..3f30519a96d 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -144,7 +144,6 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mResidualFrame(nullptr), mResidualHist(nullptr), mRatioFrame(nullptr), - mResidualFrameForCalculation(nullptr), mWorkspace(nullptr), mIntegralHisto(0), mIntegralBkg(0), diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 5c387fa48c8..96763ff81c7 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -259,16 +259,15 @@ class HFInvMassFitter : public TNamed RooPlot* mResidualFrame; /// residual frame RooHist* mResidualHist; /// residual histogram RooPlot* mRatioFrame; /// fit/data ratio frame - RooPlot* mResidualFrameForCalculation; - 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 - 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 + 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 + 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 ClassDefOverride(HFInvMassFitter, 1); }; From 5d831e68f19e1df0355823c1f9ab5a79343c788d Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 11 Mar 2026 10:50:24 +0100 Subject: [PATCH 18/26] calc Bkg === 0 when NoBkg; do not allocate mRooNBkg when NoBkg --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 3f30519a96d..88f10dfda70 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -215,11 +215,8 @@ 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 - // define number of background and background fit function - const ParameterRanges rooNBkgParamRanges{0., 1.2 * mIntegralHisto, 0.3 * mIntegralHisto}; - mRooNBkg = new RooRealVar("mRooNBkg", "number of background", randomizeInitialParameter(rooNBkgParamRanges), rooNBkgParamRanges.lower, rooNBkgParamRanges.upper); // background yield - RooAbsPdf* bkgPdf = createBackgroundFitFunction(mWorkspace); // Create background pdf - RooAbsPdf* sgnPdf = createSignalFitFunction(mWorkspace); // Create signal pdf + RooAbsPdf* bkgPdf = createBackgroundFitFunction(mWorkspace); // Create background pdf + RooAbsPdf* sgnPdf = createSignalFitFunction(mWorkspace); // Create signal pdf // fit MC or Data if (mTypeOfBkgPdf == NoBkg) { // MC @@ -240,6 +237,8 @@ void HFInvMassFitter::doFit() mRatioFrame = mass->frame(Title(Form("%s", mHistoInvMass->GetTitle()))); calculateFitToDataRatio(); } else { // data + const ParameterRanges rooNBkgParamRanges{0., 1.2 * mIntegralHisto, 0.3 * mIntegralHisto}; + 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 if (strcmp(mFitOption.c_str(), "Chi2") == 0) { @@ -740,6 +739,11 @@ void HFInvMassFitter::calculateSignal(double& signal, double& errSignal) const // calculate background yield void HFInvMassFitter::calculateBackground(double& bkg, double& errBkg) const { + if (mTypeOfBkgPdf == NoBkg) { + bkg = 0.; + errBkg = 0.; + return; + } bkg = mRooNBkg->getVal() * mIntegralBkg; errBkg = mRooNBkg->getError() * mIntegralBkg; } From fcb2721c5c8310a48689ea194b78e644ae54b4ab Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 11 Mar 2026 11:51:26 +0100 Subject: [PATCH 19/26] bring to order N sigma [sgn, sidebands] --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 6 +++--- PWGHF/D2H/Macros/HFInvMassFitter.h | 1 + PWGHF/D2H/Macros/runMassFitter.C | 4 ++++ 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 88f10dfda70..7fda7335a56 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -583,9 +583,9 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& textFitMetrics->AddText(Form("S = %.0f #pm %.0f ", mRawYield, mRawYieldErr)); textFitMetrics->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr)); if (mTypeOfBkgPdf != NoBkg) { - textFitMetrics->AddText(Form("B (%.1f#sigma) = %.0f #pm %.0f", mNSigmaForSidebands, mBkgYield, mBkgYieldErr)); - textFitMetrics->AddText(Form("S/B (%.1f#sigma) = %.4g ", mNSigmaForSidebands, mRawYield / mBkgYield)); - textFitMetrics->AddText(Form("Significance (%.1f#sigma) = %.1f #pm %.1f ", mNSigmaForSidebands, mSignificance, mSignificanceErr)); + textFitMetrics->AddText(Form("B (%.1f#sigma) = %.0f #pm %.0f", mNSigmaForSgn, mBkgYield, mBkgYieldErr)); + textFitMetrics->AddText(Form("S/B (%.1f#sigma) = %.4g ", mNSigmaForSgn, mRawYield / mBkgYield)); + textFitMetrics->AddText(Form("Significance (%.1f#sigma) = %.1f #pm %.1f ", mNSigmaForSgn, mSignificance, mSignificanceErr)); textFitMetrics->AddText(Form("#chi^{2} / ndf = %.3f", mChiSquareOverNdfTotal)); } if (mReflPdf != nullptr) { diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 96763ff81c7..6841cf15a55 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -108,6 +108,7 @@ class HFInvMassFitter : public TNamed void setFixRatioToGausSigma(double sigmaFrac); void setFixSignalYield(double yield) { mFixedRawYield = yield; } void setNumberOfSigmaForSidebands(double numberOfSigma) { mNSigmaForSidebands = numberOfSigma; } + void setNumberOfSigmaForSignal(double numberOfSigma) { mNSigmaForSgn = numberOfSigma; } void setFixDscbAlphaL(double alphaL); void setFixDscbAlphaR(double alphaR); void setFixDscbNL(double nL); diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 801854782e9..ab8133b247d 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -183,6 +183,8 @@ void runMassFitter(const std::string& configFileName) 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.); readJsonVector(dscbAlphaLInitial, config, "DscbAlphaLInitial"); readJsonVector(dscbAlphaLLower, config, "DscbAlphaLLower"); @@ -475,6 +477,8 @@ void runMassFitter(const std::string& configFileName) HFInvMassFitter* massFitter = new HFInvMassFitter(hMass[iSliceVar], massMin[iSliceVar], massMax[iSliceVar], bkgFunc[iSliceVar], sgnFunc[iSliceVar], randomSeed); massFitter->setDrawBgPrefit(drawBgPrefit); + massFitter->setNumberOfSigmaForSidebands(nSigmaForSideband); + massFitter->setNumberOfSigmaForSignal(nSigmaForSignal); massFitter->setHighlightPeakRegion(highlightPeakRegion); massFitter->setInitialGaussianMean(massPDG); massFitter->setParticlePdgMass(massPDG); From 8163558ef7ff2ab3b51cb1bfef78645e65e709e1 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 11 Mar 2026 11:54:25 +0100 Subject: [PATCH 20/26] greek mu and sigma in plot text --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 7fda7335a56..875616f1018 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -612,10 +612,10 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& const std::string fixMeanStatus = mFixedMean ? "fixed" : "free"; const std::string fixSigmaStatus = mFixedSigma ? "fixed" : "free"; const std::string fixSigmaDoubleGausStatus = mFixedSigmaDoubleGaus ? "fixed" : "free"; - textSignalPar->AddText(Form("mean(%s) = %.3f #pm %.3f", fixMeanStatus.c_str(), mRooMeanSgn->getVal(), mRooMeanSgn->getError())); - textSignalPar->AddText(Form("sigma(%s) = %.3f #pm %.3f", fixSigmaStatus.c_str(), mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textSignalPar->AddText(Form("#mu(%s) = %.3f #pm %.3f", fixMeanStatus.c_str(), mRooMeanSgn->getVal(), mRooMeanSgn->getError())); + textSignalPar->AddText(Form("#sigma(%s) = %.3f #pm %.3f", fixSigmaStatus.c_str(), mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); if (mTypeOfSgnPdf == DoubleGaus) { - textSignalPar->AddText(Form("sigma 2(%s) = %.3f #pm %.3f", fixSigmaDoubleGausStatus.c_str(), mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); + textSignalPar->AddText(Form("#sigma_{2}(%s) = %.3f #pm %.3f", fixSigmaDoubleGausStatus.c_str(), mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); } mInvMassFrame->addObject(textSignalPar); } @@ -642,10 +642,10 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad) textInfo->SetTextColor(kBlue); textInfo->AddText(Form("S = %.0f #pm %.0f ", mRawYield, mRawYieldErr)); textInfo->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr)); - textInfo->AddText(Form("mean = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError())); - textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textInfo->AddText(Form("#mu = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError())); + textInfo->AddText(Form("#sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); if (mTypeOfSgnPdf == DoubleGaus) { - textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); + textInfo->AddText(Form("#sigma_{2} = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); } mResidualFrame->addObject(textInfo); mResidualFrame->Draw(); From d3d2474eab1b368f8b78a6529bb6966f695c538c Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 11 Mar 2026 12:02:36 +0100 Subject: [PATCH 21/26] add NSigmaFor[Sideband, Signal] configurables --- PWGHF/D2H/Macros/config_massfitter.json | 2 ++ 1 file changed, 2 insertions(+) diff --git a/PWGHF/D2H/Macros/config_massfitter.json b/PWGHF/D2H/Macros/config_massfitter.json index cb84f8aaafa..0baf94cabad 100644 --- a/PWGHF/D2H/Macros/config_massfitter.json +++ b/PWGHF/D2H/Macros/config_massfitter.json @@ -163,6 +163,8 @@ "3 for GausSec", "4 for DoubleSidedCrystalBall" ], + "NSigmaForSideband": 3, + "NSigmaForSignal": 3, "DrawBgPrefit": true, "HighlightPeakRegion": true, "RandomSeed": -1 From 5685b57fd81975434f7c7edfb79e269c4289012e Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 1 Apr 2026 13:06:15 +0200 Subject: [PATCH 22/26] add check for correct parsing the config --- PWGHF/D2H/Macros/runMassFitter.C | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index ab8133b247d..92d722a44d0 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -92,6 +92,10 @@ void runMassFitter(const std::string& configFileName) config.ParseStream(is); fclose(configFile); + if (config.HasParseError()) { + 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"); From 35d09694d6f70f36720eb258d30355acd2f8b6b6 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Wed, 1 Apr 2026 13:31:02 +0200 Subject: [PATCH 23/26] bugfix *nullptr when bin-count if NoBkg --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 875616f1018..a67c495093c 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -703,8 +703,10 @@ void HFInvMassFitter::drawReflection(TVirtualPad* pad) // calculate signal yield via bin counting void HFInvMassFitter::countSignal(double& signal, double& signalErr) const { - const double binWidth = mResidualHist->GetX()[1] - mResidualHist->GetX()[0]; - const double firstBinLowEdge = mResidualHist->GetX()[0] - binWidth / 2; + const auto& histoForCounting = mTypeOfBkgPdf == NoBkg ? mInvMassFrame->getHist("data_c") : mResidualHist; + + const double binWidth = histoForCounting->GetX()[1] - histoForCounting->GetX()[0]; + const double firstBinLowEdge = histoForCounting->GetX()[0] - binWidth / 2; const auto [minForSgn, maxForSgn] = getRangesOfSignal(); const int binForMinSgn = static_cast((minForSgn - firstBinLowEdge) / binWidth) + 1; const int binForMaxSgn = static_cast((maxForSgn - firstBinLowEdge) / binWidth) + 1; @@ -716,14 +718,14 @@ void HFInvMassFitter::countSignal(double& signal, double& signalErr) const auto square = [](double value) { return value * value; }; double sumValues{}, sumErrorsSquare{}; - sumValues += mResidualHist->GetY()[binForMinSgn - 1] * binForMinSgnFraction; - sumErrorsSquare += square(mResidualHist->GetErrorY(binForMinSgn - 1) * binForMinSgnFraction); + sumValues += histoForCounting->GetY()[binForMinSgn - 1] * binForMinSgnFraction; + sumErrorsSquare += square(histoForCounting->GetErrorY(binForMinSgn - 1) * binForMinSgnFraction); for (int iBin = binForMinSgn + 1; iBin <= binForMaxSgn - 1; iBin++) { - sumValues += mResidualHist->GetY()[iBin - 1]; - sumErrorsSquare += square(mResidualHist->GetErrorY(iBin - 1)); + sumValues += histoForCounting->GetY()[iBin - 1]; + sumErrorsSquare += square(histoForCounting->GetErrorY(iBin - 1)); } - sumValues += mResidualHist->GetY()[binForMaxSgn - 1] * binForMaxSgnFraction; - sumErrorsSquare += square(mResidualHist->GetErrorY(binForMaxSgn - 1) * binForMaxSgnFraction); + sumValues += histoForCounting->GetY()[binForMaxSgn - 1] * binForMaxSgnFraction; + sumErrorsSquare += square(histoForCounting->GetErrorY(binForMaxSgn - 1) * binForMaxSgnFraction); signal = sumValues; signalErr = std::sqrt(sumErrorsSquare); From dca01b77541ef5307786008021f4deb6e8d52855 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Tue, 7 Apr 2026 15:26:24 +0200 Subject: [PATCH 24/26] add randomSeed to hFitConfig --- PWGHF/D2H/Macros/runMassFitter.C | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 92d722a44d0..3b05aa6f9ed 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -381,10 +381,11 @@ void runMassFitter(const std::string& configFileName) ConfigFixSigma, ConfigBkgFunc, ConfigSgnFunc, + ConfigRandomSeed, NConfigsToSave }; 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"}; + const char* hFitConfigXLabel[NConfigsToSave - 1] = {"mass min", "mass max", "rebin num", "fix sigma", "bkg func", "sgn func", "rnd seed"}; hFitConfig->SetStats(false); for (int i = 0; i < NConfigsToSave - 1; i++) { hFitConfig->GetXaxis()->SetBinLabel(i + 1, hFitConfigXLabel[i]); @@ -635,6 +636,7 @@ 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); } // save output histograms From 284a373a29ff7b4fbf0661265d555290e4bb808b Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Sun, 19 Apr 2026 00:06:03 +0200 Subject: [PATCH 25/26] format error message when randomization hangs --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index a67c495093c..78a50a19e85 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -1139,8 +1139,9 @@ double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& paramet result = mRandomGen->Gaus(parameterRanges.initial, sigma); ++nIter; if (nIter > MaximalNumberOfIterations) { - printf("randomizeInitialParameter() - long while loop with lower = %f upper = %f initial = %f sigma = %f\n", parameterRanges.lower, parameterRanges.upper, parameterRanges.initial, sigma); - throw; + 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); } } while (result < parameterRanges.lower || result > parameterRanges.upper); From 4793f970ca62c9c15d527652fc78cc322dd9db14 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Sun, 19 Apr 2026 00:10:14 +0200 Subject: [PATCH 26/26] add o2-linter exception (not an O2Physics workflow) --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 78a50a19e85..304f7b50d7a 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -19,6 +19,8 @@ /// \author Oleksii Lubynets /// \author Phil Stahlhut +// o2-linter: disable=name/workflow-file (not an O2Physics workflow) + #include "HFInvMassFitter.h" #include