Skip to content

Commit 5488306

Browse files
pstahlhupstahlhu
committed
Improve handling of second sigma for double Gaussian fits
1 parent 19475fa commit 5488306

File tree

3 files changed

+70
-43
lines changed

3 files changed

+70
-43
lines changed

‎PWGHF/D2H/Macros/HFInvMassFitter.cxx‎

Lines changed: 34 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -74,7 +74,7 @@ HFInvMassFitter::HFInvMassFitter() : mHistoInvMass(nullptr),
7474
mNSigmaForSidebands(4.),
7575
mNSigmaForSgn(3.),
7676
mSigmaSgnErr(0.),
77-
mSigmaSgnDoubleGaus(0.012),
77+
mSigmaSgnDoubleGaus(0.025),
7878
mFixedMean(kFALSE),
7979
mBoundMean(kFALSE),
8080
mBoundReflMean(kFALSE),
@@ -103,6 +103,7 @@ HFInvMassFitter::HFInvMassFitter() : mHistoInvMass(nullptr),
103103
mFixReflOverSgn(kFALSE),
104104
mRooMeanSgn(nullptr),
105105
mRooSigmaSgn(nullptr),
106+
mRooSecSigmaSgn(nullptr),
106107
mSgnPdf(nullptr),
107108
mBkgPdf(nullptr),
108109
mReflPdf(nullptr),
@@ -145,7 +146,7 @@ HFInvMassFitter::HFInvMassFitter(const TH1* histoToFit, Double_t minValue, Doubl
145146
mNSigmaForSidebands(3.),
146147
mNSigmaForSgn(3.),
147148
mSigmaSgnErr(0.),
148-
mSigmaSgnDoubleGaus(0.012),
149+
mSigmaSgnDoubleGaus(0.025),
149150
mFixedMean(kFALSE),
150151
mBoundMean(kFALSE),
151152
mBoundReflMean(kFALSE),
@@ -174,6 +175,7 @@ HFInvMassFitter::HFInvMassFitter(const TH1* histoToFit, Double_t minValue, Doubl
174175
mFixReflOverSgn(kFALSE),
175176
mRooMeanSgn(nullptr),
176177
mRooSigmaSgn(nullptr),
178+
mRooSecSigmaSgn(nullptr),
177179
mSgnPdf(nullptr),
178180
mBkgPdf(nullptr),
179181
mReflPdf(nullptr),
@@ -207,6 +209,7 @@ HFInvMassFitter::~HFInvMassFitter()
207209
delete mHistoTemplateRefl;
208210
delete mRooMeanSgn;
209211
delete mRooSigmaSgn;
212+
delete mRooSecSigmaSgn;
210213
delete mSgnPdf;
211214
delete mBkgPdf;
212215
delete mReflPdf;
@@ -363,7 +366,7 @@ void HFInvMassFitter::doFit()
363366
mResidualFrame->addPlotable(residualHistogram, "P");
364367
mSgnPdf->plotOn(mResidualFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineColor(kBlue));
365368
}
366-
mass->setRange("bkgForSignificance", mRooMeanSgn->getVal() - mNSigmaForSgn * mRooSigmaSgn->getVal(), mRooMeanSgn->getVal() + mNSigmaForSgn * mRooSigmaSgn->getVal());
369+
mass->setRange("bkgForSignificance", mRooMeanSgn->getVal() - mNSigmaForSgn * mRooSecSigmaSgn->getVal(), mRooMeanSgn->getVal() + mNSigmaForSgn * mRooSecSigmaSgn->getVal());
367370
bkgIntegral = mBkgPdf->createIntegral(*mass, NormSet(*mass), Range("bkgForSignificance"));
368371
mIntegralBkg = bkgIntegral->getValV();
369372
calculateBackground(mBkgYield, mBkgYieldErr);
@@ -373,6 +376,8 @@ void HFInvMassFitter::doFit()
373376
calculateSignal(mRawYield, mRawYieldErr);
374377
countSignal(mRawYieldCounted, mRawYieldCountedErr);
375378
calculateSignificance(mSignificance, mSignificanceErr);
379+
// Fit to data ratio
380+
mRatioFrame = mass->frame(Title("Fit/Data Ratio"));
376381
calculateFitToDataRatio();
377382
}
378383
}
@@ -442,10 +447,10 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const
442447
workspace.import(*sgnFuncGaus);
443448
delete sgnFuncGaus;
444449
// signal double Gaussian
445-
RooRealVar sigmaDoubleGaus("sigmaDoubleGaus", "sigma2Gaus", mSigmaSgn, mSigmaSgn - 0.01, mSigmaSgn + 0.01);
450+
RooRealVar sigmaDoubleGaus("sigmaDoubleGaus", "sigma2Gaus", mSigmaSgnDoubleGaus, mSigmaSgnDoubleGaus - 0.003, mSigmaSgnDoubleGaus + 0.003);
446451
if (mBoundSigma) {
447-
sigmaDoubleGaus.setMax(mSigmaSgn * (1 + mParamSgn));
448-
sigmaDoubleGaus.setMin(mSigmaSgn * (1 - mParamSgn));
452+
sigmaDoubleGaus.setMax(mSigmaSgnDoubleGaus * (1 + mParamSgn));
453+
sigmaDoubleGaus.setMin(mSigmaSgnDoubleGaus * (1 - mParamSgn));
449454
}
450455
if (mFixedSigma) {
451456
sigma.setVal(mSigmaSgn);
@@ -554,6 +559,7 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const
554559
workspace.import(*reflFuncPoly6);
555560
delete reflFuncPoly6;
556561
}
562+
557563
// draw fit output
558564
void HFInvMassFitter::drawFit(TVirtualPad* pad, Int_t writeFitInfo)
559565
{
@@ -588,13 +594,15 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, Int_t writeFitInfo)
588594
textInfoRight->AddText(Form("mean(free) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
589595
}
590596
if (mTypeOfSgnPdf == DoubleGaus) {
591-
auto const& baseSigmaSgn = mWorkspace->var("sigma");
597+
if (mFixedSigma) {
598+
textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
599+
} else {
600+
textInfoRight->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
601+
}
592602
if (mFixedSigmaDoubleGaus) {
593-
textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
594-
textInfoRight->AddText(Form("sigma 2(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
603+
textInfoRight->AddText(Form("sigma 2(fixed) = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError()));
595604
} else {
596-
textInfoRight->AddText(Form("sigma(free) = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
597-
textInfoRight->AddText(Form("sigma 2(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
605+
textInfoRight->AddText(Form("sigma 2(free) = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError()));
598606
}
599607
} else if (mFixedSigma) {
600608
textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
@@ -628,9 +636,8 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad)
628636
textInfo->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr));
629637
textInfo->AddText(Form("mean = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
630638
if (mTypeOfSgnPdf == DoubleGaus) {
631-
auto const& baseSigmaSgn = mWorkspace->var("sigma");
632-
textInfo->AddText(Form("sigma = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
633-
textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
639+
textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
640+
textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError()));
634641
} else {
635642
textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
636643
}
@@ -652,9 +659,8 @@ void HFInvMassFitter::drawRatio(TVirtualPad* pad)
652659
textInfo->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr));
653660
textInfo->AddText(Form("mean = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
654661
if (mTypeOfSgnPdf == DoubleGaus) {
655-
auto const& baseSigmaSgn = mWorkspace->var("sigma");
656-
textInfo->AddText(Form("sigma = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
657-
textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
662+
textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
663+
textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError()));
658664
} else {
659665
textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
660666
}
@@ -679,7 +685,7 @@ void HFInvMassFitter::highlightPeakRegion(const RooPlot* plot, Color_t color, Wi
679685
double const yMin = plot->GetMinimum();
680686
double const yMax = plot->GetMaximum();
681687
const Double_t mean = mRooMeanSgn->getVal();
682-
const Double_t sigma = mRooSigmaSgn->getVal();
688+
const Double_t sigma = mRooSecSigmaSgn->getVal();
683689
const Double_t minForSgn = mean - mNSigmaForSidebands * sigma;
684690
const Double_t maxForSgn = mean + mNSigmaForSidebands * sigma;
685691
auto* leftLine = new TLine(minForSgn, yMin, minForSgn, yMax);
@@ -704,7 +710,7 @@ void HFInvMassFitter::drawReflection(TVirtualPad* pad)
704710
void HFInvMassFitter::countSignal(Double_t& signal, Double_t& signalErr) const
705711
{
706712
const Double_t mean = mRooMeanSgn->getVal();
707-
const Double_t sigma = mRooSigmaSgn->getVal();
713+
const Double_t sigma = mRooSecSigmaSgn->getVal();
708714
const Double_t minForSgn = mean - mNSigmaForSidebands * sigma;
709715
const Double_t maxForSgn = mean + mNSigmaForSidebands * sigma;
710716
const Int_t binForMinSgn = mHistoInvMass->FindBin(minForSgn);
@@ -796,21 +802,25 @@ RooAbsPdf* HFInvMassFitter::createSignalFitFunction(RooWorkspace* workspace)
796802
case 0: {
797803
sgnPdf = workspace->pdf("sgnFuncGaus");
798804
mRooSigmaSgn = workspace->var("sigma");
805+
mRooSecSigmaSgn = workspace->var("sigma");
799806
mRooMeanSgn = workspace->var("mean");
800807
} break;
801808
case 1: {
802809
sgnPdf = workspace->pdf("sgnFuncDoubleGaus");
803-
mRooSigmaSgn = workspace->var("sigmaDoubleGaus");
810+
mRooSigmaSgn = workspace->var("sigma");
811+
mRooSecSigmaSgn = workspace->var("sigmaDoubleGaus");
804812
mRooMeanSgn = workspace->var("mean");
805813
} break;
806814
case 2: {
807815
sgnPdf = workspace->pdf("sgnFuncGausRatio");
808-
mRooSigmaSgn = workspace->var("sigmaDoubleGausRatio");
816+
mRooSigmaSgn = workspace->var("sigma");
817+
mRooSecSigmaSgn = workspace->var("sigmaDoubleGausRatio");
809818
mRooMeanSgn = workspace->var("mean");
810819
} break;
811820
case 3: {
812821
sgnPdf = workspace->pdf("sgnFuncDoublePeak");
813-
mRooSigmaSgn = workspace->var("sigmaSec");
822+
mRooSigmaSgn = workspace->var("sigma");
823+
mRooSecSigmaSgn = workspace->var("sigmaSec");
814824
mRooMeanSgn = workspace->var("meanSec");
815825
} break;
816826
default:
@@ -856,8 +866,6 @@ void HFInvMassFitter::calculateFitToDataRatio()
856866
{
857867
if (!mInvMassFrame)
858868
return;
859-
// Create a new RooPlot for the ratio
860-
mRatioFrame = mWorkspace->var("mass")->frame(Title("Fit/Data Ratio"));
861869

862870
// Get the data and fit curves from the frame
863871
auto* dataHist = dynamic_cast<RooHist*>(mInvMassFrame->findObject("data_c"));
@@ -882,8 +890,8 @@ void HFInvMassFitter::calculateFitToDataRatio()
882890
ratioHist->SetPointError(i, 0, 0, err, err);
883891
}
884892

885-
mRatioFrame->SetMinimum(0.0);
886-
mRatioFrame->SetMaximum(2.0);
893+
mRatioFrame->SetMinimum(0.5);
894+
mRatioFrame->SetMaximum(1.5);
887895
mRatioFrame->addPlotable(ratioHist, "P");
888896
}
889897

‎PWGHF/D2H/Macros/HFInvMassFitter.h‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -212,7 +212,10 @@ class HFInvMassFitter : public TNamed
212212
[[nodiscard]] Double_t getMeanUncertainty() const { return mRooMeanSgn->getError(); }
213213
[[nodiscard]] Double_t getSigma() const { return mRooSigmaSgn->getVal(); }
214214
[[nodiscard]] Double_t getSigmaUncertainty() const { return mRooSigmaSgn->getError(); }
215+
[[nodiscard]] Double_t getSecSigma() const { return mRooSecSigmaSgn->getVal(); }
216+
[[nodiscard]] Double_t getSecSigmaUncertainty() const { return mRooSecSigmaSgn->getError(); }
215217
[[nodiscard]] Double_t getReflOverSig() const
218+
216219
{
217220
if (mReflPdf != nullptr) {
218221
return mReflOverSgn;
@@ -284,6 +287,7 @@ class HFInvMassFitter : public TNamed
284287
Bool_t mFixReflOverSgn; /// switch for fix refl/signal
285288
RooRealVar* mRooMeanSgn; /// mean for gaussian of signal
286289
RooRealVar* mRooSigmaSgn; /// sigma for gaussian of signal
290+
RooRealVar* mRooSecSigmaSgn; /// second sigma for composite gaussian of signal
287291
RooAbsPdf* mSgnPdf; /// signal fit function
288292
RooAbsPdf* mBkgPdf; /// background fit function
289293
RooAbsPdf* mReflPdf; /// reflection fit function

‎PWGHF/D2H/Macros/runMassFitter.C‎

Lines changed: 32 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -269,7 +269,10 @@ int runMassFitter(const TString& configFileName)
269269
auto* hRawYieldsSigma = new TH1D(
270270
"hRawYieldsSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})",
271271
nSliceVarBins, sliceVarLimits.data());
272-
auto* hRawYieldsMean = new TH1D(
272+
auto hRawYieldsSecSigma = new TH1D(
273+
"hRawYieldsSecSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})",
274+
nSliceVarBins, sliceVarLimits.data());
275+
auto hRawYieldsMean = new TH1D(
273276
"hRawYieldsMean", ";" + sliceVarName + "(" + sliceVarUnit + ");mean (GeV/#it{c}^{2})",
274277
nSliceVarBins, sliceVarLimits.data());
275278
auto* hRawYieldsSignificance = new TH1D(
@@ -305,6 +308,7 @@ int runMassFitter(const TString& configFileName)
305308
setHistoStyle(hRawYieldsSignal);
306309
setHistoStyle(hRawYieldsSignalCounted);
307310
setHistoStyle(hRawYieldsSigma);
311+
setHistoStyle(hRawYieldsSecSigma);
308312
setHistoStyle(hRawYieldsMean);
309313
setHistoStyle(hRawYieldsSignificance);
310314
setHistoStyle(hRawYieldsSgnOverBkg);
@@ -419,6 +423,8 @@ int runMassFitter(const TString& configFileName)
419423

420424
const Double_t sigma = massFitter->getSigma();
421425
const Double_t sigmaErr = massFitter->getSigmaUncertainty();
426+
const Double_t secSigma = massFitter->getSecSigma();
427+
const Double_t secSigmaErr = massFitter->getSecSigmaUncertainty();
422428
const Double_t mean = massFitter->getMean();
423429
const Double_t meanErr = massFitter->getMeanUncertainty();
424430
const Double_t reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg();
@@ -430,6 +436,8 @@ int runMassFitter(const TString& configFileName)
430436
hRawYieldsSignalCounted->SetBinError(iSliceVar + 1, rawYieldCountedErr);
431437
hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma);
432438
hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr);
439+
hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma);
440+
hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr);
433441
hRawYieldsMean->SetBinContent(iSliceVar + 1, mean);
434442
hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr);
435443
hRawYieldsChiSquareBkg->SetBinContent(iSliceVar + 1, reducedChiSquareBkg);
@@ -484,6 +492,8 @@ int runMassFitter(const TString& configFileName)
484492
const double rawYieldCountedErr = massFitter->getRawYieldCountedError();
485493
const double sigma = massFitter->getSigma();
486494
const double sigmaErr = massFitter->getSigmaUncertainty();
495+
const double secSigma = massFitter->getSecSigma();
496+
const double secSigmaErr = massFitter->getSecSigmaUncertainty();
487497
const double mean = massFitter->getMean();
488498
const double meanErr = massFitter->getMeanUncertainty();
489499
const double reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg();
@@ -499,6 +509,8 @@ int runMassFitter(const TString& configFileName)
499509
hRawYieldsSignalCounted->SetBinError(iSliceVar + 1, rawYieldCountedErr);
500510
hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma);
501511
hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr);
512+
hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma);
513+
hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr);
502514
hRawYieldsMean->SetBinContent(iSliceVar + 1, mean);
503515
hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr);
504516
hRawYieldsSignificance->SetBinContent(iSliceVar + 1, significance);
@@ -532,23 +544,25 @@ int runMassFitter(const TString& configFileName)
532544
canvasMass[iCanvas]->Modified();
533545
canvasMass[iCanvas]->Update();
534546

535-
if (nSliceVarBins > 1) {
536-
canvasResiduals[iCanvas]->cd(iSliceVar - nCanvasesMax * iCanvas + 1);
537-
} else {
538-
canvasResiduals[iCanvas]->cd();
539-
}
540-
massFitter->drawResidual(gPad);
541-
canvasResiduals[iCanvas]->Modified();
542-
canvasResiduals[iCanvas]->Update();
547+
if (bkgFunc[iSliceVar] != HFInvMassFitter::NoBkg) {
548+
if (nSliceVarBins > 1) {
549+
canvasResiduals[iCanvas]->cd(iSliceVar - nCanvasesMax * iCanvas + 1);
550+
} else {
551+
canvasResiduals[iCanvas]->cd();
552+
}
553+
massFitter->drawResidual(gPad);
554+
canvasResiduals[iCanvas]->Modified();
555+
canvasResiduals[iCanvas]->Update();
543556

544-
if (nSliceVarBins > 1) {
545-
canvasRatio[iCanvas]->cd(iSliceVar - nCanvasesMax * iCanvas + 1);
546-
} else {
547-
canvasRatio[iCanvas]->cd();
557+
if (nSliceVarBins > 1) {
558+
canvasRatio[iCanvas]->cd(iSliceVar - nCanvasesMax * iCanvas + 1);
559+
} else {
560+
canvasRatio[iCanvas]->cd();
561+
}
562+
massFitter->drawRatio(gPad);
563+
canvasRatio[iCanvas]->Modified();
564+
canvasRatio[iCanvas]->Update();
548565
}
549-
massFitter->drawRatio(gPad);
550-
canvasRatio[iCanvas]->Modified();
551-
canvasRatio[iCanvas]->Update();
552566
}
553567

554568
hFitConfig->SetBinContent(1, iSliceVar + 1, massMin[iSliceVar]);
@@ -584,6 +598,7 @@ int runMassFitter(const TString& configFileName)
584598
hRawYieldsSignal->Write();
585599
hRawYieldsSignalCounted->Write();
586600
hRawYieldsSigma->Write();
601+
hRawYieldsSecSigma->Write();
587602
hRawYieldsMean->Write();
588603
hRawYieldsSignificance->Write();
589604
hRawYieldsSgnOverBkg->Write();
@@ -634,7 +649,7 @@ void setHistoStyle(TH1* histo, Color_t color, Size_t markerSize)
634649
void divideCanvas(TCanvas* canvas, int nSliceVarBins)
635650
{
636651
int nCols = std::ceil(std::sqrt(nSliceVarBins));
637-
int nRows = std::ceil(double(nSliceVarBins) / nCols);
652+
int nRows = std::ceil(static_cast<double>(nSliceVarBins) / nCols);
638653
canvas->Divide(nCols, nRows);
639654
}
640655

0 commit comments

Comments
 (0)