From 19475fab4325133425af90c197dd7e710ca0395a Mon Sep 17 00:00:00 2001 From: pstahlhu Date: Tue, 26 Aug 2025 11:11:25 +0200 Subject: [PATCH 01/12] HFInvMassFitter: Fix canvas division, add ratio histograms --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 69 ++++++++++++++++++++++++++++ PWGHF/D2H/Macros/HFInvMassFitter.h | 3 ++ PWGHF/D2H/Macros/runMassFitter.C | 38 ++++++++++----- 3 files changed, 99 insertions(+), 11 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 957ec2bc523..ae5bd7534d7 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -356,6 +356,7 @@ void HFInvMassFitter::doFit() mTotalPdf->plotOn(mInvMassFrame, Name("Tot_c"), LineColor(kBlue)); mSgnPdf->plotOn(mInvMassFrame, Normalization(1.0, RooAbsReal::RelativeExpected), DrawOption("F"), FillColor(TColor::GetColorTransparent(kBlue, 0.2)), VLines()); mChiSquareOverNdfTotal = mInvMassFrame->chiSquare("Tot_c", "data_c"); // calculate reduced chi2 / DNF + // plot residual distribution mResidualFrame = mass->frame(Title("Residual Distribution")); RooHist* residualHistogram = mInvMassFrame->residHist("data_c", "Bkg_c"); @@ -372,6 +373,7 @@ void HFInvMassFitter::doFit() calculateSignal(mRawYield, mRawYieldErr); countSignal(mRawYieldCounted, mRawYieldCountedErr); calculateSignificance(mSignificance, mSignificanceErr); + calculateFitToDataRatio(); } } @@ -637,6 +639,37 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad) highlightPeakRegion(mResidualFrame); } +// draw ratio on canvas +void HFInvMassFitter::drawRatio(TVirtualPad* pad) +{ + pad->cd(); + mRatioFrame->GetYaxis()->SetTitle(""); + TPaveText* textInfo = new TPaveText(0.12, 0.65, 0.47, .89, "NDC"); + textInfo->SetBorderSize(0); + textInfo->SetFillStyle(0); + 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())); + if (mTypeOfSgnPdf == DoubleGaus) { + auto const& baseSigmaSgn = mWorkspace->var("sigma"); + textInfo->AddText(Form("sigma = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError())); + textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + } else { + textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + } + mRatioFrame->addObject(textInfo); + double xMin = mRatioFrame->GetXaxis()->GetXmin(); + double xMax = mRatioFrame->GetXaxis()->GetXmax(); + TLine* line = new TLine(xMin, 1.0, xMax, 1.0); + line->SetLineColor(kGray); + line->SetLineStyle(2); + line->SetLineWidth(2); + mRatioFrame->addObject(line); + mRatioFrame->Draw(); + highlightPeakRegion(mRatioFrame); +} + // draw peak region with vertical lines void HFInvMassFitter::highlightPeakRegion(const RooPlot* plot, Color_t color, Width_t width, Style_t style) const { @@ -818,6 +851,42 @@ void HFInvMassFitter::plotRefl(RooAbsPdf* pdf) pdf->plotOn(mInvMassFrame, Components(namesOfReflPdf.at(mTypeOfReflPdf).c_str()), Name("Refl_c"), LineColor(kGreen)); } +// Calculate fit to data ratio +void HFInvMassFitter::calculateFitToDataRatio() +{ + if (!mInvMassFrame) + return; + // Create a new RooPlot for the ratio + mRatioFrame = mWorkspace->var("mass")->frame(Title("Fit/Data Ratio")); + + // Get the data and fit curves from the frame + auto* dataHist = dynamic_cast(mInvMassFrame->findObject("data_c")); + auto* fitCurve = dynamic_cast(mInvMassFrame->findObject("Tot_c")); // or the relevant fit curve + + if (!dataHist || !fitCurve) + return; + + RooHist* ratioHist = new RooHist(); + + for (int i = 0; i < dataHist->GetN(); ++i) { + double x, dataY, dataErr; + dataHist->GetPoint(i, x, dataY); + dataErr = dataHist->GetErrorY(i); + + double fitY = fitCurve->Eval(x); + + double ratio = dataY != 0 ? fitY / dataY : 0; + double err = dataY != 0 ? ratio * dataErr / dataY : 0; + + ratioHist->SetPoint(i, x, ratio); + ratioHist->SetPointError(i, 0, 0, err, err); + } + + mRatioFrame->SetMinimum(0.0); + mRatioFrame->SetMaximum(2.0); + mRatioFrame->addPlotable(ratioHist, "P"); +} + // Fix reflection pdf void HFInvMassFitter::setReflFuncFixed() { diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 1a245fe0db5..726b962db64 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -224,8 +224,10 @@ class HFInvMassFitter : public TNamed void calculateBackground(Double_t& bkg, Double_t& bkgErr) const; void calculateSignificance(Double_t& significance, Double_t& significanceErr) const; void checkForSignal(Double_t& estimatedSignal); + void calculateFitToDataRatio(); void drawFit(TVirtualPad* c, Int_t writeFitInfo = 2); void drawResidual(TVirtualPad* c); + void drawRatio(TVirtualPad* c); void drawReflection(TVirtualPad* c); private: @@ -293,6 +295,7 @@ class HFInvMassFitter : public TNamed RooPlot* mReflFrame; /// reflection frame RooPlot* mReflOnlyFrame; /// reflection frame plot on reflection only RooPlot* mResidualFrame; /// residual frame + RooPlot* mRatioFrame; /// fit/data ratio frame RooPlot* mResidualFrameForCalculation; RooWorkspace* mWorkspace; /// workspace Double_t mIntegralHisto; /// integral of histogram to fit diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 0b926b321f0..731e568910b 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -353,6 +353,7 @@ int runMassFitter(const TString& configFileName) const Int_t nCanvases = std::ceil(static_cast(nSliceVarBins) / 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; @@ -363,9 +364,16 @@ int runMassFitter(const TString& configFileName) canvasResiduals[iCanvas] = new TCanvas(Form("canvasResiduals%d", iCanvas), Form("canvasResiduals%d", iCanvas), canvasSize[0], canvasSize[1]); divideCanvas(canvasResiduals[iCanvas], nPads); - canvasRefl[iCanvas] = new TCanvas(Form("canvasRefl%d", iCanvas), Form("canvasRefl%d", iCanvas), - canvasSize[0], canvasSize[1]); - divideCanvas(canvasRefl[iCanvas], nPads); + + canvasRatio[iCanvas] = new TCanvas(Form("canvasRatio%d", iCanvas), Form("canvasRatio%d", iCanvas), + canvasSize[0], canvasSize[1]); + divideCanvas(canvasRatio[iCanvas], nPads); + + if (enableRefl) { + canvasRefl[iCanvas] = new TCanvas(Form("canvasRefl%d", iCanvas), Form("canvasRefl%d", iCanvas), + canvasSize[0], canvasSize[1]); + divideCanvas(canvasRefl[iCanvas], nPads); + } } for (unsigned int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) { @@ -532,6 +540,15 @@ int runMassFitter(const TString& configFileName) massFitter->drawResidual(gPad); canvasResiduals[iCanvas]->Modified(); canvasResiduals[iCanvas]->Update(); + + if (nSliceVarBins > 1) { + canvasRatio[iCanvas]->cd(iSliceVar - nCanvasesMax * iCanvas + 1); + } else { + canvasRatio[iCanvas]->cd(); + } + massFitter->drawRatio(gPad); + canvasRatio[iCanvas]->Modified(); + canvasRatio[iCanvas]->Update(); } hFitConfig->SetBinContent(1, iSliceVar + 1, massMin[iSliceVar]); @@ -554,7 +571,10 @@ int runMassFitter(const TString& configFileName) canvasMass[iCanvas]->Write(); if (!isMc) { canvasResiduals[iCanvas]->Write(); - canvasRefl[iCanvas]->Write(); + canvasRatio[iCanvas]->Write(); + if (enableRefl) { + canvasRefl[iCanvas]->Write(); + } } } @@ -613,13 +633,9 @@ void setHistoStyle(TH1* histo, Color_t color, Size_t markerSize) void divideCanvas(TCanvas* canvas, int nSliceVarBins) { - const int rectangularSideMin = std::floor(std::sqrt(nSliceVarBins)); - constexpr int RectangularSidesDiffMax = 2; - for (int rectangularSidesDiff = 0; rectangularSidesDiff < RectangularSidesDiffMax; ++rectangularSidesDiff) { - if (rectangularSideMin * (rectangularSideMin + rectangularSidesDiff) >= nSliceVarBins) { - canvas->Divide(rectangularSideMin + rectangularSidesDiff, rectangularSideMin); - } - } + int nCols = std::ceil(std::sqrt(nSliceVarBins)); + int nRows = std::ceil(double(nSliceVarBins) / nCols); + canvas->Divide(nCols, nRows); } int main(int argc, const char* argv[]) From 54883060ade14003f5fe3192bbbb71777eaff44d Mon Sep 17 00:00:00 2001 From: pstahlhu Date: Sun, 14 Sep 2025 19:41:07 +0200 Subject: [PATCH 02/12] Improve handling of second sigma for double Gaussian fits --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 60 ++++++++++++++++------------ PWGHF/D2H/Macros/HFInvMassFitter.h | 4 ++ PWGHF/D2H/Macros/runMassFitter.C | 49 +++++++++++++++-------- 3 files changed, 70 insertions(+), 43 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index ae5bd7534d7..b7c132eaf74 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -74,7 +74,7 @@ HFInvMassFitter::HFInvMassFitter() : mHistoInvMass(nullptr), mNSigmaForSidebands(4.), mNSigmaForSgn(3.), mSigmaSgnErr(0.), - mSigmaSgnDoubleGaus(0.012), + mSigmaSgnDoubleGaus(0.025), mFixedMean(kFALSE), mBoundMean(kFALSE), mBoundReflMean(kFALSE), @@ -103,6 +103,7 @@ HFInvMassFitter::HFInvMassFitter() : mHistoInvMass(nullptr), mFixReflOverSgn(kFALSE), mRooMeanSgn(nullptr), mRooSigmaSgn(nullptr), + mRooSecSigmaSgn(nullptr), mSgnPdf(nullptr), mBkgPdf(nullptr), mReflPdf(nullptr), @@ -145,7 +146,7 @@ HFInvMassFitter::HFInvMassFitter(const TH1* histoToFit, Double_t minValue, Doubl mNSigmaForSidebands(3.), mNSigmaForSgn(3.), mSigmaSgnErr(0.), - mSigmaSgnDoubleGaus(0.012), + mSigmaSgnDoubleGaus(0.025), mFixedMean(kFALSE), mBoundMean(kFALSE), mBoundReflMean(kFALSE), @@ -174,6 +175,7 @@ HFInvMassFitter::HFInvMassFitter(const TH1* histoToFit, Double_t minValue, Doubl mFixReflOverSgn(kFALSE), mRooMeanSgn(nullptr), mRooSigmaSgn(nullptr), + mRooSecSigmaSgn(nullptr), mSgnPdf(nullptr), mBkgPdf(nullptr), mReflPdf(nullptr), @@ -207,6 +209,7 @@ HFInvMassFitter::~HFInvMassFitter() delete mHistoTemplateRefl; delete mRooMeanSgn; delete mRooSigmaSgn; + delete mRooSecSigmaSgn; delete mSgnPdf; delete mBkgPdf; delete mReflPdf; @@ -363,7 +366,7 @@ void HFInvMassFitter::doFit() mResidualFrame->addPlotable(residualHistogram, "P"); mSgnPdf->plotOn(mResidualFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineColor(kBlue)); } - mass->setRange("bkgForSignificance", mRooMeanSgn->getVal() - mNSigmaForSgn * mRooSigmaSgn->getVal(), mRooMeanSgn->getVal() + mNSigmaForSgn * mRooSigmaSgn->getVal()); + mass->setRange("bkgForSignificance", mRooMeanSgn->getVal() - mNSigmaForSgn * mRooSecSigmaSgn->getVal(), mRooMeanSgn->getVal() + mNSigmaForSgn * mRooSecSigmaSgn->getVal()); bkgIntegral = mBkgPdf->createIntegral(*mass, NormSet(*mass), Range("bkgForSignificance")); mIntegralBkg = bkgIntegral->getValV(); calculateBackground(mBkgYield, mBkgYieldErr); @@ -373,6 +376,8 @@ void HFInvMassFitter::doFit() calculateSignal(mRawYield, mRawYieldErr); countSignal(mRawYieldCounted, mRawYieldCountedErr); calculateSignificance(mSignificance, mSignificanceErr); + // Fit to data ratio + mRatioFrame = mass->frame(Title("Fit/Data Ratio")); calculateFitToDataRatio(); } } @@ -442,10 +447,10 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const workspace.import(*sgnFuncGaus); delete sgnFuncGaus; // signal double Gaussian - RooRealVar sigmaDoubleGaus("sigmaDoubleGaus", "sigma2Gaus", mSigmaSgn, mSigmaSgn - 0.01, mSigmaSgn + 0.01); + RooRealVar sigmaDoubleGaus("sigmaDoubleGaus", "sigma2Gaus", mSigmaSgnDoubleGaus, mSigmaSgnDoubleGaus - 0.003, mSigmaSgnDoubleGaus + 0.003); if (mBoundSigma) { - sigmaDoubleGaus.setMax(mSigmaSgn * (1 + mParamSgn)); - sigmaDoubleGaus.setMin(mSigmaSgn * (1 - mParamSgn)); + sigmaDoubleGaus.setMax(mSigmaSgnDoubleGaus * (1 + mParamSgn)); + sigmaDoubleGaus.setMin(mSigmaSgnDoubleGaus * (1 - mParamSgn)); } if (mFixedSigma) { sigma.setVal(mSigmaSgn); @@ -554,6 +559,7 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const workspace.import(*reflFuncPoly6); delete reflFuncPoly6; } + // draw fit output void HFInvMassFitter::drawFit(TVirtualPad* pad, Int_t writeFitInfo) { @@ -588,13 +594,15 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, Int_t writeFitInfo) textInfoRight->AddText(Form("mean(free) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError())); } if (mTypeOfSgnPdf == DoubleGaus) { - auto const& baseSigmaSgn = mWorkspace->var("sigma"); + if (mFixedSigma) { + textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + } else { + textInfoRight->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + } if (mFixedSigmaDoubleGaus) { - textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError())); - textInfoRight->AddText(Form("sigma 2(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textInfoRight->AddText(Form("sigma 2(fixed) = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); } else { - textInfoRight->AddText(Form("sigma(free) = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError())); - textInfoRight->AddText(Form("sigma 2(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textInfoRight->AddText(Form("sigma 2(free) = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); } } else if (mFixedSigma) { textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); @@ -628,9 +636,8 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad) textInfo->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr)); textInfo->AddText(Form("mean = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError())); if (mTypeOfSgnPdf == DoubleGaus) { - auto const& baseSigmaSgn = mWorkspace->var("sigma"); - textInfo->AddText(Form("sigma = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError())); - textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); } else { textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); } @@ -652,9 +659,8 @@ void HFInvMassFitter::drawRatio(TVirtualPad* pad) textInfo->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr)); textInfo->AddText(Form("mean = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError())); if (mTypeOfSgnPdf == DoubleGaus) { - auto const& baseSigmaSgn = mWorkspace->var("sigma"); - textInfo->AddText(Form("sigma = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError())); - textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); } else { textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); } @@ -679,7 +685,7 @@ void HFInvMassFitter::highlightPeakRegion(const RooPlot* plot, Color_t color, Wi double const yMin = plot->GetMinimum(); double const yMax = plot->GetMaximum(); const Double_t mean = mRooMeanSgn->getVal(); - const Double_t sigma = mRooSigmaSgn->getVal(); + const Double_t sigma = mRooSecSigmaSgn->getVal(); const Double_t minForSgn = mean - mNSigmaForSidebands * sigma; const Double_t maxForSgn = mean + mNSigmaForSidebands * sigma; auto* leftLine = new TLine(minForSgn, yMin, minForSgn, yMax); @@ -704,7 +710,7 @@ void HFInvMassFitter::drawReflection(TVirtualPad* pad) void HFInvMassFitter::countSignal(Double_t& signal, Double_t& signalErr) const { const Double_t mean = mRooMeanSgn->getVal(); - const Double_t sigma = mRooSigmaSgn->getVal(); + const Double_t sigma = mRooSecSigmaSgn->getVal(); const Double_t minForSgn = mean - mNSigmaForSidebands * sigma; const Double_t maxForSgn = mean + mNSigmaForSidebands * sigma; const Int_t binForMinSgn = mHistoInvMass->FindBin(minForSgn); @@ -796,21 +802,25 @@ RooAbsPdf* HFInvMassFitter::createSignalFitFunction(RooWorkspace* workspace) case 0: { sgnPdf = workspace->pdf("sgnFuncGaus"); mRooSigmaSgn = workspace->var("sigma"); + mRooSecSigmaSgn = workspace->var("sigma"); mRooMeanSgn = workspace->var("mean"); } break; case 1: { sgnPdf = workspace->pdf("sgnFuncDoubleGaus"); - mRooSigmaSgn = workspace->var("sigmaDoubleGaus"); + mRooSigmaSgn = workspace->var("sigma"); + mRooSecSigmaSgn = workspace->var("sigmaDoubleGaus"); mRooMeanSgn = workspace->var("mean"); } break; case 2: { sgnPdf = workspace->pdf("sgnFuncGausRatio"); - mRooSigmaSgn = workspace->var("sigmaDoubleGausRatio"); + mRooSigmaSgn = workspace->var("sigma"); + mRooSecSigmaSgn = workspace->var("sigmaDoubleGausRatio"); mRooMeanSgn = workspace->var("mean"); } break; case 3: { sgnPdf = workspace->pdf("sgnFuncDoublePeak"); - mRooSigmaSgn = workspace->var("sigmaSec"); + mRooSigmaSgn = workspace->var("sigma"); + mRooSecSigmaSgn = workspace->var("sigmaSec"); mRooMeanSgn = workspace->var("meanSec"); } break; default: @@ -856,8 +866,6 @@ void HFInvMassFitter::calculateFitToDataRatio() { if (!mInvMassFrame) return; - // Create a new RooPlot for the ratio - mRatioFrame = mWorkspace->var("mass")->frame(Title("Fit/Data Ratio")); // Get the data and fit curves from the frame auto* dataHist = dynamic_cast(mInvMassFrame->findObject("data_c")); @@ -882,8 +890,8 @@ void HFInvMassFitter::calculateFitToDataRatio() ratioHist->SetPointError(i, 0, 0, err, err); } - mRatioFrame->SetMinimum(0.0); - mRatioFrame->SetMaximum(2.0); + mRatioFrame->SetMinimum(0.5); + mRatioFrame->SetMaximum(1.5); mRatioFrame->addPlotable(ratioHist, "P"); } diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 726b962db64..112e1c6a0db 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -212,7 +212,10 @@ class HFInvMassFitter : public TNamed [[nodiscard]] Double_t getMeanUncertainty() const { return mRooMeanSgn->getError(); } [[nodiscard]] Double_t getSigma() const { return mRooSigmaSgn->getVal(); } [[nodiscard]] Double_t getSigmaUncertainty() const { return mRooSigmaSgn->getError(); } + [[nodiscard]] Double_t getSecSigma() const { return mRooSecSigmaSgn->getVal(); } + [[nodiscard]] Double_t getSecSigmaUncertainty() const { return mRooSecSigmaSgn->getError(); } [[nodiscard]] Double_t getReflOverSig() const + { if (mReflPdf != nullptr) { return mReflOverSgn; @@ -284,6 +287,7 @@ class HFInvMassFitter : public TNamed Bool_t mFixReflOverSgn; /// switch for fix refl/signal RooRealVar* mRooMeanSgn; /// mean for gaussian of signal RooRealVar* mRooSigmaSgn; /// sigma for gaussian of signal + RooRealVar* mRooSecSigmaSgn; /// second sigma for composite gaussian of signal RooAbsPdf* mSgnPdf; /// signal fit function RooAbsPdf* mBkgPdf; /// background fit function RooAbsPdf* mReflPdf; /// reflection fit function diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 731e568910b..825be0ef7e7 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -269,7 +269,10 @@ int runMassFitter(const TString& configFileName) auto* hRawYieldsSigma = new TH1D( "hRawYieldsSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsMean = new TH1D( + auto hRawYieldsSecSigma = new TH1D( + "hRawYieldsSecSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", + nSliceVarBins, sliceVarLimits.data()); + auto hRawYieldsMean = new TH1D( "hRawYieldsMean", ";" + sliceVarName + "(" + sliceVarUnit + ");mean (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data()); auto* hRawYieldsSignificance = new TH1D( @@ -305,6 +308,7 @@ int runMassFitter(const TString& configFileName) setHistoStyle(hRawYieldsSignal); setHistoStyle(hRawYieldsSignalCounted); setHistoStyle(hRawYieldsSigma); + setHistoStyle(hRawYieldsSecSigma); setHistoStyle(hRawYieldsMean); setHistoStyle(hRawYieldsSignificance); setHistoStyle(hRawYieldsSgnOverBkg); @@ -419,6 +423,8 @@ int runMassFitter(const TString& configFileName) const Double_t sigma = massFitter->getSigma(); const Double_t sigmaErr = massFitter->getSigmaUncertainty(); + const Double_t secSigma = massFitter->getSecSigma(); + const Double_t secSigmaErr = massFitter->getSecSigmaUncertainty(); const Double_t mean = massFitter->getMean(); const Double_t meanErr = massFitter->getMeanUncertainty(); const Double_t reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg(); @@ -430,6 +436,8 @@ int runMassFitter(const TString& configFileName) hRawYieldsSignalCounted->SetBinError(iSliceVar + 1, rawYieldCountedErr); hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma); hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr); + hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma); + hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr); hRawYieldsMean->SetBinContent(iSliceVar + 1, mean); hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr); hRawYieldsChiSquareBkg->SetBinContent(iSliceVar + 1, reducedChiSquareBkg); @@ -484,6 +492,8 @@ int runMassFitter(const TString& configFileName) const double rawYieldCountedErr = massFitter->getRawYieldCountedError(); const double sigma = massFitter->getSigma(); const double sigmaErr = massFitter->getSigmaUncertainty(); + const double secSigma = massFitter->getSecSigma(); + const double secSigmaErr = massFitter->getSecSigmaUncertainty(); const double mean = massFitter->getMean(); const double meanErr = massFitter->getMeanUncertainty(); const double reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg(); @@ -499,6 +509,8 @@ int runMassFitter(const TString& configFileName) hRawYieldsSignalCounted->SetBinError(iSliceVar + 1, rawYieldCountedErr); hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma); hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr); + hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma); + hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr); hRawYieldsMean->SetBinContent(iSliceVar + 1, mean); hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr); hRawYieldsSignificance->SetBinContent(iSliceVar + 1, significance); @@ -532,23 +544,25 @@ int runMassFitter(const TString& configFileName) canvasMass[iCanvas]->Modified(); canvasMass[iCanvas]->Update(); - if (nSliceVarBins > 1) { - canvasResiduals[iCanvas]->cd(iSliceVar - nCanvasesMax * iCanvas + 1); - } else { - canvasResiduals[iCanvas]->cd(); - } - massFitter->drawResidual(gPad); - canvasResiduals[iCanvas]->Modified(); - canvasResiduals[iCanvas]->Update(); + if (bkgFunc[iSliceVar] != HFInvMassFitter::NoBkg) { + if (nSliceVarBins > 1) { + canvasResiduals[iCanvas]->cd(iSliceVar - nCanvasesMax * iCanvas + 1); + } else { + canvasResiduals[iCanvas]->cd(); + } + massFitter->drawResidual(gPad); + canvasResiduals[iCanvas]->Modified(); + canvasResiduals[iCanvas]->Update(); - if (nSliceVarBins > 1) { - canvasRatio[iCanvas]->cd(iSliceVar - nCanvasesMax * iCanvas + 1); - } else { - canvasRatio[iCanvas]->cd(); + if (nSliceVarBins > 1) { + canvasRatio[iCanvas]->cd(iSliceVar - nCanvasesMax * iCanvas + 1); + } else { + canvasRatio[iCanvas]->cd(); + } + massFitter->drawRatio(gPad); + canvasRatio[iCanvas]->Modified(); + canvasRatio[iCanvas]->Update(); } - massFitter->drawRatio(gPad); - canvasRatio[iCanvas]->Modified(); - canvasRatio[iCanvas]->Update(); } hFitConfig->SetBinContent(1, iSliceVar + 1, massMin[iSliceVar]); @@ -584,6 +598,7 @@ int runMassFitter(const TString& configFileName) hRawYieldsSignal->Write(); hRawYieldsSignalCounted->Write(); hRawYieldsSigma->Write(); + hRawYieldsSecSigma->Write(); hRawYieldsMean->Write(); hRawYieldsSignificance->Write(); hRawYieldsSgnOverBkg->Write(); @@ -634,7 +649,7 @@ void setHistoStyle(TH1* histo, Color_t color, Size_t markerSize) void divideCanvas(TCanvas* canvas, int nSliceVarBins) { int nCols = std::ceil(std::sqrt(nSliceVarBins)); - int nRows = std::ceil(double(nSliceVarBins) / nCols); + int nRows = std::ceil(static_cast(nSliceVarBins) / nCols); canvas->Divide(nCols, nRows); } From b464a701dd2c9735174877cd2ce7507f66f4c52d Mon Sep 17 00:00:00 2001 From: pstahlhu Date: Tue, 16 Sep 2025 12:23:41 +0200 Subject: [PATCH 03/12] Store fraction for double Gaussian and add possibility to fix it --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 5 ++ PWGHF/D2H/Macros/HFInvMassFitter.h | 5 +- PWGHF/D2H/Macros/runMassFitter.C | 71 ++++++++++++++++++---------- 3 files changed, 55 insertions(+), 26 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index b7c132eaf74..23beac876be 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -104,6 +104,7 @@ HFInvMassFitter::HFInvMassFitter() : mHistoInvMass(nullptr), mRooMeanSgn(nullptr), mRooSigmaSgn(nullptr), mRooSecSigmaSgn(nullptr), + mRooFracDoubleGaus(nullptr), mSgnPdf(nullptr), mBkgPdf(nullptr), mReflPdf(nullptr), @@ -176,6 +177,7 @@ HFInvMassFitter::HFInvMassFitter(const TH1* histoToFit, Double_t minValue, Doubl mRooMeanSgn(nullptr), mRooSigmaSgn(nullptr), mRooSecSigmaSgn(nullptr), + mRooFracDoubleGaus(nullptr), mSgnPdf(nullptr), mBkgPdf(nullptr), mReflPdf(nullptr), @@ -210,6 +212,7 @@ HFInvMassFitter::~HFInvMassFitter() delete mRooMeanSgn; delete mRooSigmaSgn; delete mRooSecSigmaSgn; + delete mRooFracDoubleGaus; delete mSgnPdf; delete mBkgPdf; delete mReflPdf; @@ -810,12 +813,14 @@ RooAbsPdf* HFInvMassFitter::createSignalFitFunction(RooWorkspace* workspace) mRooSigmaSgn = workspace->var("sigma"); mRooSecSigmaSgn = workspace->var("sigmaDoubleGaus"); mRooMeanSgn = workspace->var("mean"); + mRooFracDoubleGaus = workspace->var("fracDoubleGaus"); } break; case 2: { sgnPdf = workspace->pdf("sgnFuncGausRatio"); mRooSigmaSgn = workspace->var("sigma"); mRooSecSigmaSgn = workspace->var("sigmaDoubleGausRatio"); mRooMeanSgn = workspace->var("mean"); + mRooFracDoubleGaus = workspace->var("fracDoubleGausRatio"); } break; case 3: { sgnPdf = workspace->pdf("sgnFuncDoublePeak"); diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 112e1c6a0db..2141f855087 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -214,8 +214,10 @@ class HFInvMassFitter : public TNamed [[nodiscard]] Double_t getSigmaUncertainty() const { return mRooSigmaSgn->getError(); } [[nodiscard]] Double_t getSecSigma() const { return mRooSecSigmaSgn->getVal(); } [[nodiscard]] Double_t getSecSigmaUncertainty() const { return mRooSecSigmaSgn->getError(); } + [[nodiscard]] Double_t getFracDoubleGaus() const { return mRooFracDoubleGaus->getVal(); } + [[nodiscard]] Double_t getFracDoubleGausUncertainty() const { return mRooFracDoubleGaus->getError(); } [[nodiscard]] Double_t getReflOverSig() const - + { if (mReflPdf != nullptr) { return mReflOverSgn; @@ -288,6 +290,7 @@ class HFInvMassFitter : public TNamed RooRealVar* mRooMeanSgn; /// mean for gaussian of signal RooRealVar* mRooSigmaSgn; /// sigma for gaussian of signal RooRealVar* mRooSecSigmaSgn; /// second sigma for composite gaussian of signal + RooRealVar* mRooFracDoubleGaus; /// fraction of second gaussian for composite gaussian of signal RooAbsPdf* mSgnPdf; /// signal fit function RooAbsPdf* mBkgPdf; /// background fit function RooAbsPdf* mReflPdf; /// reflection fit function diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 825be0ef7e7..790bc2fcb2e 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -105,9 +105,10 @@ int runMassFitter(const TString& configFileName) std::vector sliceVarMax; std::vector massMin; std::vector massMax; + std::vector fixMeanManual; std::vector fixSigmaManual; std::vector fixSecondSigmaManual; - std::vector fixMeanManual; + std::vector fixFracDoubleGausManual; std::vector nRebin; std::vector bkgFuncConfig; std::vector sgnFuncConfig; @@ -130,12 +131,15 @@ int runMassFitter(const TString& configFileName) const Value& fdSecPeakHistoNameValue = config["FDSecPeakHistoName"]; parseStringArray(fdSecPeakHistoNameValue, fdSecPeakHistoName); - const bool fixSigma = config["FixSigma"].GetBool(); - const std::string sigmaFile = config["SigmaFile"].GetString(); - const bool fixMean = config["FixMean"].GetBool(); const std::string meanFile = config["MeanFile"].GetString(); + const Value& fixMeanManualValue = config["FixMeanManual"]; + readArray(fixMeanManualValue, fixMeanManual); + + const bool fixSigma = config["FixSigma"].GetBool(); + const std::string sigmaFile = config["SigmaFile"].GetString(); + const Value& fixSigmaManualValue = config["FixSigmaManual"]; readArray(fixSigmaManualValue, fixSigmaManual); @@ -145,8 +149,11 @@ int runMassFitter(const TString& configFileName) const Value& fixSecondSigmaManualValue = config["FixSecondSigmaManual"]; readArray(fixSecondSigmaManualValue, fixSecondSigmaManual); - const Value& fixMeanManualValue = config["FixMeanManual"]; - readArray(fixMeanManualValue, fixMeanManual); + const bool fixFracDoubleGaus = config["FixFracDoubleGaus"].GetBool(); + const std::string fracDoubleGausFile = config["FracDoubleGausFile"].GetString(); + + const Value& fixFracDoubleGausManualValue = config["FixFracDoubleGausManual"]; + readArray(fixFracDoubleGausManualValue, fixFracDoubleGausManual); sliceVarName = config["SliceVarName"].GetString(); sliceVarUnit = config["SliceVarUnit"].GetString(); @@ -262,18 +269,21 @@ int runMassFitter(const TString& configFileName) inputFile->Close(); // 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* hRawYieldsSigma = new TH1D( + 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 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 hRawYieldsMean = new TH1D( - "hRawYieldsMean", ";" + sliceVarName + "(" + sliceVarUnit + ");mean (GeV/#it{c}^{2})", + auto hRawYieldsFracDoubleGaus = new TH1D( + "hRawYieldsFracDoubleGaus", ";" + sliceVarName + "(" + sliceVarUnit + ");fraction of double gaussian", nSliceVarBins, sliceVarLimits.data()); auto* hRawYieldsSignificance = new TH1D( "hRawYieldsSignificance", @@ -307,9 +317,10 @@ int runMassFitter(const TString& configFileName) setHistoStyle(hRawYieldsSignal); setHistoStyle(hRawYieldsSignalCounted); + setHistoStyle(hRawYieldsMean); setHistoStyle(hRawYieldsSigma); setHistoStyle(hRawYieldsSecSigma); - setHistoStyle(hRawYieldsMean); + setHistoStyle(hRawYieldsFracDoubleGaus); setHistoStyle(hRawYieldsSignificance); setHistoStyle(hRawYieldsSgnOverBkg); setHistoStyle(hRawYieldsBkg); @@ -339,7 +350,8 @@ int runMassFitter(const TString& configFileName) TH1* hSigmaToFix = getHistToFix(fixSigma, fixSigmaManual, sigmaFile, "Sigma"); TH1* hMeanToFix = getHistToFix(fixMean, fixMeanManual, meanFile, "Mean"); - TH1* hSecondSigmaToFix = getHistToFix(fixSecondSigma, fixSecondSigmaManual, secondSigmaFile, "Sigma"); + TH1* hSecondSigmaToFix = getHistToFix(fixSecondSigma, fixSecondSigmaManual, secondSigmaFile, "SecSigma"); + TH1* hFracDoubleGausToFix = getHistToFix(fixFracDoubleGaus, fixFracDoubleGausManual, fracDoubleGausFile, "FracDoubleGaus"); // fit histograms @@ -420,13 +432,14 @@ int runMassFitter(const TString& configFileName) const Double_t rawYieldErr = massFitter->getRawYieldError(); const Double_t rawYieldCounted = massFitter->getRawYieldCounted(); const Double_t rawYieldCountedErr = massFitter->getRawYieldCountedError(); - + const Double_t mean = massFitter->getMean(); + const Double_t meanErr = massFitter->getMeanUncertainty(); const Double_t sigma = massFitter->getSigma(); const Double_t sigmaErr = massFitter->getSigmaUncertainty(); const Double_t secSigma = massFitter->getSecSigma(); const Double_t secSigmaErr = massFitter->getSecSigmaUncertainty(); - const Double_t mean = massFitter->getMean(); - const Double_t meanErr = massFitter->getMeanUncertainty(); + const Double_t fracDoubleGaus = massFitter->getFracDoubleGaus(); + const Double_t fracDoubleGausErr = massFitter->getFracDoubleGausUncertainty(); const Double_t reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg(); const Double_t reducedChiSquareTotal = massFitter->getChiSquareOverNDFTotal(); @@ -434,12 +447,14 @@ int runMassFitter(const TString& configFileName) hRawYieldsSignal->SetBinError(iSliceVar + 1, rawYieldErr); hRawYieldsSignalCounted->SetBinContent(iSliceVar + 1, rawYieldCounted); hRawYieldsSignalCounted->SetBinError(iSliceVar + 1, rawYieldCountedErr); + hRawYieldsMean->SetBinContent(iSliceVar + 1, mean); + hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr); hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma); hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr); hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma); hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr); - hRawYieldsMean->SetBinContent(iSliceVar + 1, mean); - hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr); + hRawYieldsFracDoubleGaus->SetBinContent(iSliceVar + 1, fracDoubleGaus); + hRawYieldsFracDoubleGaus->SetBinError(iSliceVar + 1, fracDoubleGausErr); hRawYieldsChiSquareBkg->SetBinContent(iSliceVar + 1, reducedChiSquareBkg); hRawYieldsChiSquareBkg->SetBinError(iSliceVar + 1, 0.); hRawYieldsChiSquareTotal->SetBinContent(iSliceVar + 1, reducedChiSquareTotal); @@ -476,6 +491,7 @@ int runMassFitter(const TString& configFileName) setFixedValue(fixMean, fixMeanManual, hMeanToFix, std::bind(&HFInvMassFitter::setFixGaussianMean, massFitter, std::placeholders::_1), "MEAN"); setFixedValue(fixSigma, fixSigmaManual, hSigmaToFix, std::bind(&HFInvMassFitter::setFixGaussianSigma, massFitter, std::placeholders::_1), "SIGMA"); setFixedValue(fixSecondSigma, fixSecondSigmaManual, hSecondSigmaToFix, std::bind(&HFInvMassFitter::setFixSecondGaussianSigma, massFitter, std::placeholders::_1), "SECOND SIGMA"); + setFixedValue(fixFracDoubleGaus, fixFracDoubleGausManual, nullptr, std::bind(&HFInvMassFitter::setFixFrac2Gaus, massFitter, std::placeholders::_1), "FRAC DOUBLE GAUS"); if (enableRefl) { reflOverSgn = hMassForSgn[iSliceVar]->Integral(hMassForSgn[iSliceVar]->FindBin(massMin[iSliceVar] * 1.0001), hMassForSgn[iSliceVar]->FindBin(massMax[iSliceVar] * 0.999)); @@ -490,12 +506,14 @@ int runMassFitter(const TString& configFileName) const double rawYieldErr = massFitter->getRawYieldError(); const double rawYieldCounted = massFitter->getRawYieldCounted(); const double rawYieldCountedErr = massFitter->getRawYieldCountedError(); + const double mean = massFitter->getMean(); + const double meanErr = massFitter->getMeanUncertainty(); const double sigma = massFitter->getSigma(); const double sigmaErr = massFitter->getSigmaUncertainty(); const double secSigma = massFitter->getSecSigma(); const double secSigmaErr = massFitter->getSecSigmaUncertainty(); - const double mean = massFitter->getMean(); - const double meanErr = massFitter->getMeanUncertainty(); + const double fracDoubleGaus = massFitter->getFracDoubleGaus(); + const double fracDoubleGausErr = massFitter->getFracDoubleGausUncertainty(); const double reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg(); const double reducedChiSquareTotal = massFitter->getChiSquareOverNDFTotal(); const double significance = massFitter->getSignificance(); @@ -507,12 +525,14 @@ int runMassFitter(const TString& configFileName) hRawYieldsSignal->SetBinError(iSliceVar + 1, rawYieldErr); hRawYieldsSignalCounted->SetBinContent(iSliceVar + 1, rawYieldCounted); hRawYieldsSignalCounted->SetBinError(iSliceVar + 1, rawYieldCountedErr); + hRawYieldsMean->SetBinContent(iSliceVar + 1, mean); + hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr); hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma); hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr); hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma); hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr); - hRawYieldsMean->SetBinContent(iSliceVar + 1, mean); - hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr); + hRawYieldsFracDoubleGaus->SetBinContent(iSliceVar + 1, fracDoubleGaus); + hRawYieldsFracDoubleGaus->SetBinError(iSliceVar + 1, fracDoubleGausErr); hRawYieldsSignificance->SetBinContent(iSliceVar + 1, significance); hRawYieldsSignificance->SetBinError(iSliceVar + 1, significanceErr); hRawYieldsSgnOverBkg->SetBinContent(iSliceVar + 1, rawYield / bkg); @@ -597,9 +617,10 @@ int runMassFitter(const TString& configFileName) } hRawYieldsSignal->Write(); hRawYieldsSignalCounted->Write(); + hRawYieldsMean->Write(); hRawYieldsSigma->Write(); hRawYieldsSecSigma->Write(); - hRawYieldsMean->Write(); + hRawYieldsFracDoubleGaus->Write(); hRawYieldsSignificance->Write(); hRawYieldsSgnOverBkg->Write(); hRawYieldsBkg->Write(); From 92823fceda2d7c20d80fa6175bd7a7394adc5a40 Mon Sep 17 00:00:00 2001 From: pstahlhu Date: Tue, 16 Sep 2025 18:15:03 +0200 Subject: [PATCH 04/12] Store signal parameters to hist only when present --- PWGHF/D2H/Macros/runMassFitter.C | 93 ++++++++++++++++++-------------- 1 file changed, 52 insertions(+), 41 deletions(-) diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 790bc2fcb2e..6fbd8d2604e 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -432,33 +432,38 @@ int runMassFitter(const TString& configFileName) const Double_t rawYieldErr = massFitter->getRawYieldError(); const Double_t rawYieldCounted = massFitter->getRawYieldCounted(); const Double_t rawYieldCountedErr = massFitter->getRawYieldCountedError(); + const Double_t reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg(); + const Double_t reducedChiSquareTotal = massFitter->getChiSquareOverNDFTotal(); const Double_t mean = massFitter->getMean(); const Double_t meanErr = massFitter->getMeanUncertainty(); const Double_t sigma = massFitter->getSigma(); const Double_t sigmaErr = massFitter->getSigmaUncertainty(); - const Double_t secSigma = massFitter->getSecSigma(); - const Double_t secSigmaErr = massFitter->getSecSigmaUncertainty(); - const Double_t fracDoubleGaus = massFitter->getFracDoubleGaus(); - const Double_t fracDoubleGausErr = massFitter->getFracDoubleGausUncertainty(); - const Double_t reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg(); - const Double_t reducedChiSquareTotal = massFitter->getChiSquareOverNDFTotal(); hRawYieldsSignal->SetBinContent(iSliceVar + 1, rawYield); hRawYieldsSignal->SetBinError(iSliceVar + 1, rawYieldErr); hRawYieldsSignalCounted->SetBinContent(iSliceVar + 1, rawYieldCounted); hRawYieldsSignalCounted->SetBinError(iSliceVar + 1, rawYieldCountedErr); - hRawYieldsMean->SetBinContent(iSliceVar + 1, mean); - hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr); - hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma); - hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr); - hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma); - hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr); - hRawYieldsFracDoubleGaus->SetBinContent(iSliceVar + 1, fracDoubleGaus); - hRawYieldsFracDoubleGaus->SetBinError(iSliceVar + 1, fracDoubleGausErr); hRawYieldsChiSquareBkg->SetBinContent(iSliceVar + 1, reducedChiSquareBkg); hRawYieldsChiSquareBkg->SetBinError(iSliceVar + 1, 0.); hRawYieldsChiSquareTotal->SetBinContent(iSliceVar + 1, reducedChiSquareTotal); hRawYieldsChiSquareTotal->SetBinError(iSliceVar + 1, 0.); + hRawYieldsMean->SetBinContent(iSliceVar + 1, mean); + hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr); + hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma); + hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr); + + if (sgnFunc[iSliceVar] != HFInvMassFitter::SingleGaus) { + const Double_t secSigma = massFitter->getSecSigma(); + const Double_t secSigmaErr = massFitter->getSecSigmaUncertainty(); + hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma); + hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr); + } + if (sgnFunc[iSliceVar] == HFInvMassFitter::DoubleGaus || sgnFunc[iSliceVar] == HFInvMassFitter::DoubleGausSigmaRatioPar) { + const Double_t fracDoubleGaus = massFitter->getFracDoubleGaus(); + const Double_t fracDoubleGausErr = massFitter->getFracDoubleGausUncertainty(); + hRawYieldsFracDoubleGaus->SetBinContent(iSliceVar + 1, fracDoubleGaus); + hRawYieldsFracDoubleGaus->SetBinError(iSliceVar + 1, fracDoubleGausErr); + } } else { HFInvMassFitter* massFitter; massFitter = new HFInvMassFitter(hMassForFit[iSliceVar], massMin[iSliceVar], massMax[iSliceVar], @@ -506,43 +511,49 @@ int runMassFitter(const TString& configFileName) const double rawYieldErr = massFitter->getRawYieldError(); const double rawYieldCounted = massFitter->getRawYieldCounted(); const double rawYieldCountedErr = massFitter->getRawYieldCountedError(); + const double bkg = massFitter->getBkgYield(); + const double bkgErr = massFitter->getBkgYieldError(); + const double significance = massFitter->getSignificance(); + const double significanceErr = massFitter->getSignificanceError(); + const double reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg(); + const double reducedChiSquareTotal = massFitter->getChiSquareOverNDFTotal(); const double mean = massFitter->getMean(); const double meanErr = massFitter->getMeanUncertainty(); const double sigma = massFitter->getSigma(); const double sigmaErr = massFitter->getSigmaUncertainty(); - const double secSigma = massFitter->getSecSigma(); - const double secSigmaErr = massFitter->getSecSigmaUncertainty(); - const double fracDoubleGaus = massFitter->getFracDoubleGaus(); - const double fracDoubleGausErr = massFitter->getFracDoubleGausUncertainty(); - const double reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg(); - const double reducedChiSquareTotal = massFitter->getChiSquareOverNDFTotal(); - const double significance = massFitter->getSignificance(); - const double significanceErr = massFitter->getSignificanceError(); - const double bkg = massFitter->getBkgYield(); - const double bkgErr = massFitter->getBkgYieldError(); hRawYieldsSignal->SetBinContent(iSliceVar + 1, rawYield); hRawYieldsSignal->SetBinError(iSliceVar + 1, rawYieldErr); hRawYieldsSignalCounted->SetBinContent(iSliceVar + 1, rawYieldCounted); hRawYieldsSignalCounted->SetBinError(iSliceVar + 1, rawYieldCountedErr); - hRawYieldsMean->SetBinContent(iSliceVar + 1, mean); - hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr); - hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma); - hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr); - hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma); - hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr); - hRawYieldsFracDoubleGaus->SetBinContent(iSliceVar + 1, fracDoubleGaus); - hRawYieldsFracDoubleGaus->SetBinError(iSliceVar + 1, fracDoubleGausErr); - hRawYieldsSignificance->SetBinContent(iSliceVar + 1, significance); - hRawYieldsSignificance->SetBinError(iSliceVar + 1, significanceErr); - hRawYieldsSgnOverBkg->SetBinContent(iSliceVar + 1, rawYield / bkg); - hRawYieldsSgnOverBkg->SetBinError(iSliceVar + 1, rawYield / bkg * std::sqrt(rawYieldErr / rawYield * rawYieldErr / rawYield + bkgErr / bkg * bkgErr / bkg)); hRawYieldsBkg->SetBinContent(iSliceVar + 1, bkg); hRawYieldsBkg->SetBinError(iSliceVar + 1, bkgErr); + hRawYieldsSgnOverBkg->SetBinContent(iSliceVar + 1, rawYield / bkg); + hRawYieldsSgnOverBkg->SetBinError(iSliceVar + 1, rawYield / bkg * std::sqrt(rawYieldErr / rawYield * rawYieldErr / rawYield + bkgErr / bkg * bkgErr / bkg)); + hRawYieldsSignificance->SetBinContent(iSliceVar + 1, significance); + hRawYieldsSignificance->SetBinError(iSliceVar + 1, significanceErr); hRawYieldsChiSquareBkg->SetBinContent(iSliceVar + 1, reducedChiSquareBkg); hRawYieldsChiSquareBkg->SetBinError(iSliceVar + 1, 1.e-20); hRawYieldsChiSquareTotal->SetBinContent(iSliceVar + 1, reducedChiSquareTotal); hRawYieldsChiSquareTotal->SetBinError(iSliceVar + 1, 1.e-20); + hRawYieldsMean->SetBinContent(iSliceVar + 1, mean); + hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr); + hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma); + hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr); + + if (sgnFunc[iSliceVar] != HFInvMassFitter::SingleGaus) { + const double secSigma = massFitter->getSecSigma(); + const double secSigmaErr = massFitter->getSecSigmaUncertainty(); + hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma); + hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr); + } + if (sgnFunc[iSliceVar] == HFInvMassFitter::DoubleGaus || sgnFunc[iSliceVar] == HFInvMassFitter::DoubleGausSigmaRatioPar) { + const double fracDoubleGaus = massFitter->getFracDoubleGaus(); + const double fracDoubleGausErr = massFitter->getFracDoubleGausUncertainty(); + hRawYieldsFracDoubleGaus->SetBinContent(iSliceVar + 1, fracDoubleGaus); + hRawYieldsFracDoubleGaus->SetBinError(iSliceVar + 1, fracDoubleGausErr); + } + if (enableRefl) { hReflectionOverSignal->SetBinContent(iSliceVar + 1, reflOverSgn); if (nSliceVarBins > 1) { @@ -617,15 +628,15 @@ int runMassFitter(const TString& configFileName) } hRawYieldsSignal->Write(); hRawYieldsSignalCounted->Write(); + hRawYieldsBkg->Write(); + hRawYieldsSgnOverBkg->Write(); + hRawYieldsSignificance->Write(); + hRawYieldsChiSquareBkg->Write(); + hRawYieldsChiSquareTotal->Write(); hRawYieldsMean->Write(); hRawYieldsSigma->Write(); hRawYieldsSecSigma->Write(); hRawYieldsFracDoubleGaus->Write(); - hRawYieldsSignificance->Write(); - hRawYieldsSgnOverBkg->Write(); - hRawYieldsBkg->Write(); - hRawYieldsChiSquareBkg->Write(); - hRawYieldsChiSquareTotal->Write(); if (enableRefl) { hReflectionOverSignal->Write(); } From 56db18f6ebd60d7c83af47ac6b594c9b475507f1 Mon Sep 17 00:00:00 2001 From: pstahlhu Date: Fri, 19 Sep 2025 16:47:49 +0200 Subject: [PATCH 05/12] Improve plotting style of the mass canvases --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 98 ++++++++++++++++------------ PWGHF/D2H/Macros/HFInvMassFitter.h | 2 +- PWGHF/D2H/Macros/runMassFitter.C | 36 +++++----- 3 files changed, 78 insertions(+), 58 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 23beac876be..98d5979aeac 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -50,6 +50,7 @@ #include #include #include +#include using namespace RooFit; @@ -564,65 +565,80 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const } // draw fit output -void HFInvMassFitter::drawFit(TVirtualPad* pad, Int_t writeFitInfo) +void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& plotLabels, Bool_t writeParInfo) { gStyle->SetOptStat(0); gStyle->SetCanvasColor(0); gStyle->SetFrameFillColor(0); pad->cd(); - if (writeFitInfo > 0) { - auto* textInfoLeft = new TPaveText(0.12, 0.65, 0.47, 0.89, "NDC"); - auto* textInfoRight = new TPaveText(0.6, 0.7, 1., .87, "NDC"); - textInfoLeft->SetBorderSize(0); - textInfoLeft->SetFillStyle(0); - textInfoRight->SetBorderSize(0); - textInfoRight->SetFillStyle(0); - textInfoRight->SetTextColor(kBlue); - textInfoLeft->AddText(Form("S = %.0f #pm %.0f ", mRawYield, mRawYieldErr)); - textInfoLeft->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr)); - if (mTypeOfBkgPdf != NoBkg) { - textInfoLeft->AddText(Form("B (%d#sigma) = %.0f #pm %.0f", mNSigmaForSidebands, mBkgYield, mBkgYieldErr)); - textInfoLeft->AddText(Form("S/B (%d#sigma) = %.4g ", mNSigmaForSidebands, mRawYield / mBkgYield)); - } - if (mReflPdf != nullptr) { - textInfoLeft->AddText(Form("Refl/Sig = %.3f #pm %.3f ", mReflOverSgn, 0.0)); - } - if (mTypeOfBkgPdf != NoBkg) { - textInfoLeft->AddText(Form("Signif (%d#sigma) = %.1f #pm %.1f ", mNSigmaForSidebands, mSignificance, mSignificanceErr)); - textInfoLeft->AddText(Form("#chi^{2} / ndf = %.3f", mChiSquareOverNdfTotal)); - } + // Fit metrics + TPaveText* textFitMetrics = new TPaveText(0.65, 0.7, 0.9, 0.88, "NDC"); + textFitMetrics->SetBorderSize(0); + textFitMetrics->SetFillStyle(0); + textFitMetrics->SetTextSize(0.04); + textFitMetrics->SetTextAlign(33); + textFitMetrics->AddText(Form("S = %.0f #pm %.0f ", mRawYield, mRawYieldErr)); + if (mTypeOfBkgPdf != NoBkg) { + textFitMetrics->AddText(Form("B (%d#sigma) = %.0f #pm %.0f", mNSigmaForSidebands, mBkgYield, mBkgYieldErr)); + textFitMetrics->AddText(Form("S/B (%d#sigma) = %.4g ", mNSigmaForSidebands, mRawYield / mBkgYield)); + } + if (mReflPdf) { + textFitMetrics->AddText(Form("Refl/Sig = %.3f #pm %.3f ", mReflOverSgn, 0.0)); + } + if (mTypeOfBkgPdf != NoBkg) { + textFitMetrics->AddText(Form("Significance (%d#sigma) = %.1f #pm %.1f ", mNSigmaForSidebands, mSignificance, mSignificanceErr)); + textFitMetrics->AddText(Form("#chi^{2} / ndf = %.3f", mChiSquareOverNdfTotal)); + } + mInvMassFrame->addObject(textFitMetrics); + // Analysis information + TPaveText* textAnalysisInfo = new TPaveText(0.18, 0.78, 0.35, 0.88, "NDC"); + textAnalysisInfo->SetBorderSize(0); + textAnalysisInfo->SetFillStyle(0); + textAnalysisInfo->SetTextSize(0.05); + textAnalysisInfo->SetTextAlign(13); + for (const auto& label : plotLabels) { + textAnalysisInfo->AddText(label.c_str()); + } + mInvMassFrame->addObject(textAnalysisInfo); + if (writeParInfo) { + // right text box + TPaveText* textSignalPar = new TPaveText(0.18, 0.65, 0.4, 0.75, "NDC"); + textSignalPar->SetBorderSize(0); + textSignalPar->SetFillStyle(0); + textSignalPar->SetTextColor(kBlue); + textSignalPar->SetTextAlign(13); if (mFixedMean) { - textInfoRight->AddText(Form("mean(fixed) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError())); + textSignalPar->AddText(Form("mean(fixed) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError())); } else { - textInfoRight->AddText(Form("mean(free) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError())); + textSignalPar->AddText(Form("mean(free) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError())); } if (mTypeOfSgnPdf == DoubleGaus) { if (mFixedSigma) { - textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textSignalPar->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); } else { - textInfoRight->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textSignalPar->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); } if (mFixedSigmaDoubleGaus) { - textInfoRight->AddText(Form("sigma 2(fixed) = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); + textSignalPar->AddText(Form("sigma 2(fixed) = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); } else { - textInfoRight->AddText(Form("sigma 2(free) = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); + textSignalPar->AddText(Form("sigma 2(free) = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); } } else if (mFixedSigma) { - textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textSignalPar->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); } else { - textInfoRight->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); - } - mInvMassFrame->addObject(textInfoLeft); - mInvMassFrame->addObject(textInfoRight); - mInvMassFrame->GetYaxis()->SetTitleOffset(1.8); - gPad->SetLeftMargin(0.15); - mInvMassFrame->GetYaxis()->SetTitle(Form("%s", mHistoInvMass->GetYaxis()->GetTitle())); - mInvMassFrame->GetXaxis()->SetTitle(Form("%s", mHistoInvMass->GetXaxis()->GetTitle())); - mInvMassFrame->Draw(); - highlightPeakRegion(mInvMassFrame); - if (mHistoTemplateRefl != nullptr) { - mReflFrame->Draw("same"); + textSignalPar->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); } + mInvMassFrame->addObject(textSignalPar); + } + mInvMassFrame->GetXaxis()->SetTitleOffset(1.2); + mInvMassFrame->GetYaxis()->SetTitleOffset(1.8); + gPad->SetLeftMargin(0.15); + mInvMassFrame->GetYaxis()->SetTitle(Form("%s", mHistoInvMass->GetYaxis()->GetTitle())); + mInvMassFrame->GetXaxis()->SetTitle(Form("%s", mHistoInvMass->GetXaxis()->GetTitle())); + mInvMassFrame->Draw(); + highlightPeakRegion(mInvMassFrame); + if (mHistoTemplateRefl) { + mReflFrame->Draw("same"); } } diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 2141f855087..3ca571e1573 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -230,7 +230,7 @@ class HFInvMassFitter : public TNamed void calculateSignificance(Double_t& significance, Double_t& significanceErr) const; void checkForSignal(Double_t& estimatedSignal); void calculateFitToDataRatio(); - void drawFit(TVirtualPad* c, Int_t writeFitInfo = 2); + void drawFit(TVirtualPad* c, const std::vector& plotLabels, Bool_t writeParInfo = true); void drawResidual(TVirtualPad* c); void drawRatio(TVirtualPad* c); void drawReflection(TVirtualPad* c); diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 6fbd8d2604e..22197139dd1 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -87,11 +87,13 @@ int runMassFitter(const TString& configFileName) config.ParseStream(is); fclose(configFile); - Bool_t const isMc = config["IsMC"].GetBool(); - TString const inputFileName = config["InFileName"].GetString(); - TString const reflFileName = config["ReflFileName"].GetString(); + Bool_t isMc = config["IsMC"].GetBool(); + Bool_t writeSignalPar = config["WriteSignalPar"].GetBool(); + TString inputFileName = config["InFileName"].GetString(); + TString reflFileName = config["ReflFileName"].GetString(); TString outputFileName = config["OutFileName"].GetString(); - TString const particleName = config["Particle"].GetString(); + TString particleName = config["Particle"].GetString(); + TString collisionSystem = config["CollisionSystem"].GetString(); std::vector inputHistoName; std::vector promptHistoName; @@ -207,19 +209,21 @@ int runMassFitter(const TString& configFileName) sgnFunc[iSliceVar] = sgnFuncConfig[iSliceVar]; } - std::map> particles{ - {"Dplus", {"K#pi#pi", "D+"}}, - {"D0", {"K#pi", "D0"}}, - {"Ds", {"KK#pi", "D_s+"}}, - {"LcToPKPi", {"pK#pi", "Lambda_c+"}}, - {"LcToPK0s", {"pK^{0}_{s}", "Lambda_c+"}}, - {"Dstar", {"D^{0}pi^{+}", "D*+"}}, - {"XicToXiPiPi", {"#Xi#pi#pi", "Xi_c+"}}}; + std::map> particles{ + {"Dplus", {"K#pi#pi", "D+", "D^{+} #rightarrow K^{-}#pi^{+}#pi^{+} + c.c."}}, + {"D0", {"K#pi", "D0", "D^{0} #rightarrow K^{-}#pi^{+} + c.c."}}, + {"Ds", {"KK#pi", "D_s+", "D_{s}^{+} #rightarrow K^{-}K^{+}#pi^{+} + c.c."}}, + {"LcToPKPi", {"pK#pi", "Lambda_c+", "#Lambda_{c}^{+} #rightarrow pK^{-}#pi^{+} + c.c."}}, + {"LcToPK0s", {"pK^{0}_{s}", "Lambda_c+", "#Lambda_{c}^{+} #rightarrow pK^{0}_{s} + c.c."}}, + {"Dstar", {"D^{0}pi^{+}", "D*+", "D^{*+} #rightarrow D^{0}#pi^{+} + c.c."}}, + {"XicToXiPiPi", {"#Xi#pi#pi", "Xi_c+", "#Xi_{c}^{+} #rightarrow #Xi^{-}#pi^{+}#pi^{+} + c.c."}}}; if (particles.find(particleName.Data()) == particles.end()) { throw std::runtime_error("ERROR: only Dplus, D0, Ds, LcToPKPi, LcToPK0s, Dstar and XicToXiPiPi particles supported! Exit"); } - const TString massAxisTitle = "#it{M}(" + particles[particleName.Data()].first + ") (GeV/#it{c}^{2})"; - const double massPDG = TDatabasePDG::Instance()->GetParticle(particles[particleName.Data()].second.c_str())->Mass(); + const auto& particleTuple = particles[particleName.Data()]; + const TString massAxisTitle = "#it{M}(" + std::get<0>(particleTuple) + ") (GeV/#it{c}^{2})"; + const double massPDG = TDatabasePDG::Instance()->GetParticle(std::get<1>(particleTuple).c_str())->Mass(); + const std::vector plotLabels = {std::get<2>(particleTuple), collisionSystem.Data()}; // load inv-mass histograms auto* inputFile = TFile::Open(inputFileName.Data()); @@ -426,7 +430,7 @@ int runMassFitter(const TString& configFileName) canvasMass[iCanvas]->cd(); } - massFitter->drawFit(gPad); + massFitter->drawFit(gPad, plotLabels, writeSignalPar); const Double_t rawYield = massFitter->getRawYield(); const Double_t rawYieldErr = massFitter->getRawYieldError(); @@ -571,7 +575,7 @@ int runMassFitter(const TString& configFileName) } else { canvasMass[iCanvas]->cd(); } - massFitter->drawFit(gPad); + massFitter->drawFit(gPad, plotLabels, writeSignalPar); canvasMass[iCanvas]->Modified(); canvasMass[iCanvas]->Update(); From 0db0b7fc9fb35605472fdd4ebdc25be0837089d2 Mon Sep 17 00:00:00 2001 From: pstahlhu Date: Fri, 19 Sep 2025 17:27:28 +0200 Subject: [PATCH 06/12] Implement Oleksii's comments --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 24 ++++++------------------ PWGHF/D2H/Macros/HFInvMassFitter.h | 2 +- PWGHF/D2H/Macros/config_massfitter.json | 4 ++++ PWGHF/D2H/Macros/runMassFitter.C | 11 +++++++++++ 4 files changed, 22 insertions(+), 19 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 98d5979aeac..88b191c409b 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -381,7 +381,7 @@ void HFInvMassFitter::doFit() countSignal(mRawYieldCounted, mRawYieldCountedErr); calculateSignificance(mSignificance, mSignificanceErr); // Fit to data ratio - mRatioFrame = mass->frame(Title("Fit/Data Ratio")); + mRatioFrame = mass->frame(Title(Form("%s", mHistoInvMass->GetTitle()))); calculateFitToDataRatio(); } } @@ -669,21 +669,9 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad) void HFInvMassFitter::drawRatio(TVirtualPad* pad) { pad->cd(); - mRatioFrame->GetYaxis()->SetTitle(""); - TPaveText* textInfo = new TPaveText(0.12, 0.65, 0.47, .89, "NDC"); - textInfo->SetBorderSize(0); - textInfo->SetFillStyle(0); - 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())); - if (mTypeOfSgnPdf == DoubleGaus) { - textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); - textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); - } else { - textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); - } - mRatioFrame->addObject(textInfo); + mRatioFrame->GetXaxis()->SetTitleOffset(1.2); + mRatioFrame->GetYaxis()->SetTitleOffset(1.5); + mRatioFrame->GetYaxis()->SetTitle("Fit / Data"); double xMin = mRatioFrame->GetXaxis()->GetXmin(); double xMax = mRatioFrame->GetXaxis()->GetXmax(); TLine* line = new TLine(xMin, 1.0, xMax, 1.0); @@ -883,7 +871,7 @@ void HFInvMassFitter::plotRefl(RooAbsPdf* pdf) } // Calculate fit to data ratio -void HFInvMassFitter::calculateFitToDataRatio() +void HFInvMassFitter::calculateFitToDataRatio() const { if (!mInvMassFrame) return; @@ -911,9 +899,9 @@ void HFInvMassFitter::calculateFitToDataRatio() ratioHist->SetPointError(i, 0, 0, err, err); } + mRatioFrame->addPlotable(ratioHist, "P"); mRatioFrame->SetMinimum(0.5); mRatioFrame->SetMaximum(1.5); - mRatioFrame->addPlotable(ratioHist, "P"); } // Fix reflection pdf diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 3ca571e1573..6fa614e961e 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -229,7 +229,7 @@ class HFInvMassFitter : public TNamed void calculateBackground(Double_t& bkg, Double_t& bkgErr) const; void calculateSignificance(Double_t& significance, Double_t& significanceErr) const; void checkForSignal(Double_t& estimatedSignal); - void calculateFitToDataRatio(); + void calculateFitToDataRatio() const; void drawFit(TVirtualPad* c, const std::vector& plotLabels, Bool_t writeParInfo = true); void drawResidual(TVirtualPad* c); void drawRatio(TVirtualPad* c); diff --git a/PWGHF/D2H/Macros/config_massfitter.json b/PWGHF/D2H/Macros/config_massfitter.json index 43e67b76f5d..6ddb0626623 100644 --- a/PWGHF/D2H/Macros/config_massfitter.json +++ b/PWGHF/D2H/Macros/config_massfitter.json @@ -1,5 +1,7 @@ { "IsMC": false, + "WriteSignalPar": true, + "_WriteSignalPar": "write signal parameter values on mass canvas", "InFileName": "hMasses.root", "_InFileName": "download example file here https://cernbox.cern.ch/s/uoWicuRAVGArDsV", "ReflFileName": "", @@ -35,6 +37,8 @@ "Dstar", "XicToXiPiPi" ], + "CollisionSystem": "pp #sqrt{#it{s}} = 13.6 TeV", + "_CollisionSystems": "Label for mass canvases", "EnableRefl": false, "FixSigma": false, "SigmaFile": "", diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 22197139dd1..a6f170d568a 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -651,6 +651,8 @@ int runMassFitter(const TString& configFileName) outputFileName.ReplaceAll(".root", ".pdf"); TString outputFileNameResidual = outputFileName; outputFileNameResidual.ReplaceAll(".pdf", "_Residuals.pdf"); + TString outputFileRatio = outputFileName; + outputFileRatio.ReplaceAll(".pdf", "_Ratio.pdf"); for (int iCanvas = 0; iCanvas < nCanvases; iCanvas++) { if (iCanvas == 0 && nCanvases > 1) { canvasMass[iCanvas]->SaveAs(Form("%s[", outputFileName.Data())); @@ -660,6 +662,7 @@ int runMassFitter(const TString& configFileName) canvasMass[iCanvas]->SaveAs(Form("%s]", outputFileName.Data())); } if (!isMc) { + // residuals if (iCanvas == 0 && nCanvases > 1) { canvasResiduals[iCanvas]->SaveAs(Form("%s[", outputFileNameResidual.Data())); } @@ -667,6 +670,14 @@ int runMassFitter(const TString& configFileName) if (iCanvas == nCanvases - 1 && nCanvases > 1) { canvasResiduals[iCanvas]->SaveAs(Form("%s]", outputFileNameResidual.Data())); } + // ratio + if (iCanvas == 0 && nCanvases > 1) { + canvasRatio[iCanvas]->SaveAs(Form("%s[", outputFileRatio.Data())); + } + canvasRatio[iCanvas]->SaveAs(outputFileRatio.Data()); + if (iCanvas == nCanvases - 1 && nCanvases > 1) { + canvasRatio[iCanvas]->SaveAs(Form("%s]", outputFileRatio.Data())); + } } } return 0; From 4a19c4611202627b279b034d20100497ab8b3d64 Mon Sep 17 00:00:00 2001 From: Oleksii Lubynets Date: Thu, 18 Sep 2025 23:52:09 +0200 Subject: [PATCH 07/12] adapt config_massfitter for FixFracDoubleGaus --- PWGHF/D2H/Macros/config_massfitter.json | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/PWGHF/D2H/Macros/config_massfitter.json b/PWGHF/D2H/Macros/config_massfitter.json index 6ddb0626623..393fb8aa0f3 100644 --- a/PWGHF/D2H/Macros/config_massfitter.json +++ b/PWGHF/D2H/Macros/config_massfitter.json @@ -61,7 +61,7 @@ 0, 0 ], - "_FixMeanManual": "fix mean mannually", + "_FixMeanManual": "fix mean manually", "FixSecondSigma": false, "SecondSigmaFile": "", "_SecondSigmaFile": "fix second sigma for double gauss from file", @@ -73,6 +73,17 @@ 0 ], "_FixSecondSigmaManual": "fix second sigma for double gauss manually", + "FixFracDoubleGaus": false, + "FracDoubleGausFile": "", + "_FracDoubleGausFile": "fix fraction of sigmas for double gauss from file", + "FixFracDoubleGausManual": [ + 0, + 0, + 0, + 0, + 0 + ], + "_FixFracDoubleGausManual": "fix fraction of sigmas for double gauss manually", "SliceVarName": "T", "SliceVarUnit": "ps", "_SliceVarName, _SliceVarUnit": "e.g. pT, GeV/c or something else depending on user's needs", From be724be154d9fe9706d49462994a58a2e5cfce2e Mon Sep 17 00:00:00 2001 From: pstahlhu Date: Fri, 19 Sep 2025 17:56:54 +0200 Subject: [PATCH 08/12] Fix includes in runMassFitter.C --- PWGHF/D2H/Macros/runMassFitter.C | 1 + 1 file changed, 1 insertion(+) diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index a6f170d568a..54f6062435e 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -42,6 +42,7 @@ #include #include #include // std::string +#include #include #include // std::vector From 02f1bbadecae64831c9fc365e18dbb690e0da761 Mon Sep 17 00:00:00 2001 From: pstahlhu Date: Fri, 7 Nov 2025 13:43:03 +0100 Subject: [PATCH 09/12] Restore S_count metric to mass fit histograms --- 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 88b191c409b..b85a8f738d2 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -578,6 +578,7 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& textFitMetrics->SetTextSize(0.04); textFitMetrics->SetTextAlign(33); textFitMetrics->AddText(Form("S = %.0f #pm %.0f ", mRawYield, mRawYieldErr)); + textFitMetrics->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr)); if (mTypeOfBkgPdf != NoBkg) { textFitMetrics->AddText(Form("B (%d#sigma) = %.0f #pm %.0f", mNSigmaForSidebands, mBkgYield, mBkgYieldErr)); textFitMetrics->AddText(Form("S/B (%d#sigma) = %.4g ", mNSigmaForSidebands, mRawYield / mBkgYield)); From e0cd9453f4153daa2810ecdbffdf9e1c9163127c Mon Sep 17 00:00:00 2001 From: pstahlhu Date: Thu, 13 Nov 2025 21:00:24 +0100 Subject: [PATCH 10/12] Comply with Vit's clang-tidy fixes --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 17 ++++--- PWGHF/D2H/Macros/HFInvMassFitter.h | 1 + PWGHF/D2H/Macros/runMassFitter.C | 71 ++++++++++++---------------- 3 files changed, 40 insertions(+), 49 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index b85a8f738d2..0d1dc8b2405 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -17,6 +17,7 @@ /// \author Xinye Peng /// \author Biao Zhang /// \author Oleksii Lubynets +/// \author Phil Stahlhut #include "HFInvMassFitter.h" @@ -572,7 +573,7 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& gStyle->SetFrameFillColor(0); pad->cd(); // Fit metrics - TPaveText* textFitMetrics = new TPaveText(0.65, 0.7, 0.9, 0.88, "NDC"); + auto* textFitMetrics = new TPaveText(0.65, 0.7, 0.9, 0.88, "NDC"); textFitMetrics->SetBorderSize(0); textFitMetrics->SetFillStyle(0); textFitMetrics->SetTextSize(0.04); @@ -582,17 +583,15 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& if (mTypeOfBkgPdf != NoBkg) { textFitMetrics->AddText(Form("B (%d#sigma) = %.0f #pm %.0f", mNSigmaForSidebands, mBkgYield, mBkgYieldErr)); textFitMetrics->AddText(Form("S/B (%d#sigma) = %.4g ", mNSigmaForSidebands, mRawYield / mBkgYield)); - } - if (mReflPdf) { - textFitMetrics->AddText(Form("Refl/Sig = %.3f #pm %.3f ", mReflOverSgn, 0.0)); - } - if (mTypeOfBkgPdf != NoBkg) { textFitMetrics->AddText(Form("Significance (%d#sigma) = %.1f #pm %.1f ", mNSigmaForSidebands, mSignificance, mSignificanceErr)); + if (mReflPdf != nullptr) { + textFitMetrics->AddText(Form("Refl/Sig = %.3f #pm %.3f ", mReflOverSgn, 0.0)); + } textFitMetrics->AddText(Form("#chi^{2} / ndf = %.3f", mChiSquareOverNdfTotal)); } mInvMassFrame->addObject(textFitMetrics); // Analysis information - TPaveText* textAnalysisInfo = new TPaveText(0.18, 0.78, 0.35, 0.88, "NDC"); + auto* textAnalysisInfo = new TPaveText(0.18, 0.78, 0.35, 0.88, "NDC"); textAnalysisInfo->SetBorderSize(0); textAnalysisInfo->SetFillStyle(0); textAnalysisInfo->SetTextSize(0.05); @@ -603,7 +602,7 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& mInvMassFrame->addObject(textAnalysisInfo); if (writeParInfo) { // right text box - TPaveText* textSignalPar = new TPaveText(0.18, 0.65, 0.4, 0.75, "NDC"); + auto* textSignalPar = new TPaveText(0.18, 0.65, 0.4, 0.75, "NDC"); textSignalPar->SetBorderSize(0); textSignalPar->SetFillStyle(0); textSignalPar->SetTextColor(kBlue); @@ -675,7 +674,7 @@ void HFInvMassFitter::drawRatio(TVirtualPad* pad) mRatioFrame->GetYaxis()->SetTitle("Fit / Data"); double xMin = mRatioFrame->GetXaxis()->GetXmin(); double xMax = mRatioFrame->GetXaxis()->GetXmax(); - TLine* line = new TLine(xMin, 1.0, xMax, 1.0); + auto* line = new TLine(xMin, 1.0, xMax, 1.0); line->SetLineColor(kGray); line->SetLineStyle(2); line->SetLineWidth(2); diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 6fa614e961e..daa00ec0b09 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -17,6 +17,7 @@ /// \author Xinye Peng /// \author Biao Zhang /// \author Oleksii Lubynets +/// \author Phil Stahlhut #ifndef PWGHF_D2H_MACROS_HFINVMASSFITTER_H_ #define PWGHF_D2H_MACROS_HFINVMASSFITTER_H_ diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 54f6062435e..d02fe292ac4 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -17,6 +17,7 @@ /// \author Xinye Peng /// \author Biao Zhang /// \author Oleksii Lubynets +/// \author Phil Stahlhut #if !defined(__CINT__) || defined(__CLING__) @@ -88,13 +89,13 @@ int runMassFitter(const TString& configFileName) config.ParseStream(is); fclose(configFile); - Bool_t isMc = config["IsMC"].GetBool(); - Bool_t writeSignalPar = config["WriteSignalPar"].GetBool(); - TString inputFileName = config["InFileName"].GetString(); - TString reflFileName = config["ReflFileName"].GetString(); + Bool_t const isMc = config["IsMC"].GetBool(); + Bool_t const writeSignalPar = config["WriteSignalPar"].GetBool(); + TString const inputFileName = config["InFileName"].GetString(); + TString const reflFileName = config["ReflFileName"].GetString(); TString outputFileName = config["OutFileName"].GetString(); - TString particleName = config["Particle"].GetString(); - TString collisionSystem = config["CollisionSystem"].GetString(); + TString const particleName = config["Particle"].GetString(); + TString const collisionSystem = config["CollisionSystem"].GetString(); std::vector inputHistoName; std::vector promptHistoName; @@ -274,40 +275,30 @@ int runMassFitter(const TString& configFileName) inputFile->Close(); // 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 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* hRawYieldsSignificance = new TH1D( - "hRawYieldsSignificance", - ";" + sliceVarName + "(" + sliceVarUnit + ");significance (3#sigma)", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsSgnOverBkg = - new TH1D("hRawYieldsSgnOverBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");S/B (3#sigma)", - nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsBkg = - new TH1D("hRawYieldsBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");Background (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* 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* 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* hRawYieldsSignificance = new TH1D("hRawYieldsSignificance", + ";" + sliceVarName + "(" + sliceVarUnit + ");significance (3#sigma)", nSliceVarBins, sliceVarLimits.data()); + auto* hRawYieldsSgnOverBkg = new TH1D("hRawYieldsSgnOverBkg", + ";" + sliceVarName + "(" + sliceVarUnit + ");S/B (3#sigma)", nSliceVarBins, sliceVarLimits.data()); + auto* hRawYieldsBkg = new TH1D("hRawYieldsBkg", + ";" + sliceVarName + "(" + sliceVarUnit + ");Background (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()); const Int_t nConfigsToSave = 6; auto* hFitConfig = new TH2F("hfitConfig", "Fit Configurations", nConfigsToSave, 0, 6, nSliceVarBins, sliceVarLimits.data()); From f61f7be322715762652365f119d3770e65ae37a9 Mon Sep 17 00:00:00 2001 From: pstahlhu Date: Thu, 13 Nov 2025 21:31:52 +0100 Subject: [PATCH 11/12] Unify histogram order and fix if condition for plot labels --- PWGHF/D2H/Macros/HFInvMassFitter.cxx | 6 ++-- PWGHF/D2H/Macros/runMassFitter.C | 52 +++++++++++----------------- 2 files changed, 23 insertions(+), 35 deletions(-) diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 0d1dc8b2405..843712e2245 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -584,11 +584,11 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& textFitMetrics->AddText(Form("B (%d#sigma) = %.0f #pm %.0f", mNSigmaForSidebands, mBkgYield, mBkgYieldErr)); textFitMetrics->AddText(Form("S/B (%d#sigma) = %.4g ", mNSigmaForSidebands, mRawYield / mBkgYield)); textFitMetrics->AddText(Form("Significance (%d#sigma) = %.1f #pm %.1f ", mNSigmaForSidebands, mSignificance, mSignificanceErr)); - if (mReflPdf != nullptr) { - textFitMetrics->AddText(Form("Refl/Sig = %.3f #pm %.3f ", mReflOverSgn, 0.0)); - } textFitMetrics->AddText(Form("#chi^{2} / ndf = %.3f", mChiSquareOverNdfTotal)); } + if (mReflPdf != nullptr) { + textFitMetrics->AddText(Form("Refl/Sig = %.3f #pm %.3f ", mReflOverSgn, 0.0)); + } mInvMassFrame->addObject(textFitMetrics); // Analysis information auto* textAnalysisInfo = new TPaveText(0.18, 0.78, 0.35, 0.88, "NDC"); diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index d02fe292ac4..6c5b5a6ecf8 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -91,11 +91,11 @@ int runMassFitter(const TString& configFileName) Bool_t const isMc = config["IsMC"].GetBool(); Bool_t const writeSignalPar = config["WriteSignalPar"].GetBool(); + TString const particleName = config["Particle"].GetString(); + TString const collisionSystem = config["CollisionSystem"].GetString(); TString const inputFileName = config["InFileName"].GetString(); TString const reflFileName = config["ReflFileName"].GetString(); TString outputFileName = config["OutFileName"].GetString(); - TString const particleName = config["Particle"].GetString(); - TString const collisionSystem = config["CollisionSystem"].GetString(); std::vector inputHistoName; std::vector promptHistoName; @@ -275,30 +275,18 @@ int runMassFitter(const TString& configFileName) inputFile->Close(); // 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* 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* hRawYieldsSignificance = new TH1D("hRawYieldsSignificance", - ";" + sliceVarName + "(" + sliceVarUnit + ");significance (3#sigma)", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsSgnOverBkg = new TH1D("hRawYieldsSgnOverBkg", - ";" + sliceVarName + "(" + sliceVarUnit + ");S/B (3#sigma)", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsBkg = new TH1D("hRawYieldsBkg", - ";" + sliceVarName + "(" + sliceVarUnit + ");Background (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* 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()); const Int_t nConfigsToSave = 6; auto* hFitConfig = new TH2F("hfitConfig", "Fit Configurations", nConfigsToSave, 0, 6, nSliceVarBins, sliceVarLimits.data()); @@ -313,16 +301,16 @@ int runMassFitter(const TString& configFileName) setHistoStyle(hRawYieldsSignal); setHistoStyle(hRawYieldsSignalCounted); - setHistoStyle(hRawYieldsMean); - setHistoStyle(hRawYieldsSigma); - setHistoStyle(hRawYieldsSecSigma); - setHistoStyle(hRawYieldsFracDoubleGaus); - setHistoStyle(hRawYieldsSignificance); - setHistoStyle(hRawYieldsSgnOverBkg); setHistoStyle(hRawYieldsBkg); + setHistoStyle(hRawYieldsSgnOverBkg); + setHistoStyle(hRawYieldsSignificance); setHistoStyle(hRawYieldsChiSquareBkg); setHistoStyle(hRawYieldsChiSquareTotal); setHistoStyle(hReflectionOverSignal, kRed + 1); + setHistoStyle(hRawYieldsMean); + setHistoStyle(hRawYieldsSigma); + setHistoStyle(hRawYieldsSecSigma); + setHistoStyle(hRawYieldsFracDoubleGaus); auto getHistToFix = [&nSliceVarBins](bool const& isFix, std::vector const& fixManual, std::string const& fixFileName, std::string const& var) -> TH1* { TH1* histToFix = nullptr; From 6f24caae01102da97210eb248a9ba98172943f54 Mon Sep 17 00:00:00 2001 From: pstahlhu Date: Fri, 14 Nov 2025 20:46:58 +0100 Subject: [PATCH 12/12] Fix lambda to fix fraction of double Gaussian --- PWGHF/D2H/Macros/runMassFitter.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 6c5b5a6ecf8..e4baa2a856a 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -480,7 +480,7 @@ int runMassFitter(const TString& configFileName) setFixedValue(fixMean, fixMeanManual, hMeanToFix, std::bind(&HFInvMassFitter::setFixGaussianMean, massFitter, std::placeholders::_1), "MEAN"); setFixedValue(fixSigma, fixSigmaManual, hSigmaToFix, std::bind(&HFInvMassFitter::setFixGaussianSigma, massFitter, std::placeholders::_1), "SIGMA"); setFixedValue(fixSecondSigma, fixSecondSigmaManual, hSecondSigmaToFix, std::bind(&HFInvMassFitter::setFixSecondGaussianSigma, massFitter, std::placeholders::_1), "SECOND SIGMA"); - setFixedValue(fixFracDoubleGaus, fixFracDoubleGausManual, nullptr, std::bind(&HFInvMassFitter::setFixFrac2Gaus, massFitter, std::placeholders::_1), "FRAC DOUBLE GAUS"); + setFixedValue(fixFracDoubleGaus, fixFracDoubleGausManual, hFracDoubleGausToFix, std::bind(&HFInvMassFitter::setFixFrac2Gaus, massFitter, std::placeholders::_1), "FRAC DOUBLE GAUS"); if (enableRefl) { reflOverSgn = hMassForSgn[iSliceVar]->Integral(hMassForSgn[iSliceVar]->FindBin(massMin[iSliceVar] * 1.0001), hMassForSgn[iSliceVar]->FindBin(massMax[iSliceVar] * 0.999));