Skip to content

Commit 2fcfa1a

Browse files
pstahlhupstahlhu
committed
Improve handling of second sigma for double Gaussian fits
1 parent d252d43 commit 2fcfa1a

File tree

3 files changed

+68
-42
lines changed

3 files changed

+68
-42
lines changed

‎PWGHF/D2H/Macros/HFInvMassFitter.cxx‎

Lines changed: 34 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -76,7 +76,7 @@ HFInvMassFitter::HFInvMassFitter() : TNamed(),
7676
mNSigmaForSidebands(4.),
7777
mNSigmaForSgn(3.),
7878
mSigmaSgnErr(0.),
79-
mSigmaSgnDoubleGaus(0.012),
79+
mSigmaSgnDoubleGaus(0.025),
8080
mFixedMean(kFALSE),
8181
mBoundMean(kFALSE),
8282
mBoundReflMean(kFALSE),
@@ -105,6 +105,7 @@ HFInvMassFitter::HFInvMassFitter() : TNamed(),
105105
mFixReflOverSgn(kFALSE),
106106
mRooMeanSgn(nullptr),
107107
mRooSigmaSgn(nullptr),
108+
mRooSecSigmaSgn(nullptr),
108109
mSgnPdf(nullptr),
109110
mBkgPdf(nullptr),
110111
mReflPdf(nullptr),
@@ -148,7 +149,7 @@ HFInvMassFitter::HFInvMassFitter(const TH1* histoToFit, Double_t minValue, Doubl
148149
mNSigmaForSidebands(3.),
149150
mNSigmaForSgn(3.),
150151
mSigmaSgnErr(0.),
151-
mSigmaSgnDoubleGaus(0.012),
152+
mSigmaSgnDoubleGaus(0.025),
152153
mFixedMean(kFALSE),
153154
mBoundMean(kFALSE),
154155
mBoundReflMean(kFALSE),
@@ -177,6 +178,7 @@ HFInvMassFitter::HFInvMassFitter(const TH1* histoToFit, Double_t minValue, Doubl
177178
mFixReflOverSgn(kFALSE),
178179
mRooMeanSgn(nullptr),
179180
mRooSigmaSgn(nullptr),
181+
mRooSecSigmaSgn(nullptr),
180182
mSgnPdf(nullptr),
181183
mBkgPdf(nullptr),
182184
mReflPdf(nullptr),
@@ -210,6 +212,7 @@ HFInvMassFitter::~HFInvMassFitter()
210212
delete mHistoTemplateRefl;
211213
delete mRooMeanSgn;
212214
delete mRooSigmaSgn;
215+
delete mRooSecSigmaSgn;
213216
delete mSgnPdf;
214217
delete mBkgPdf;
215218
delete mReflPdf;
@@ -366,7 +369,7 @@ void HFInvMassFitter::doFit()
366369
mResidualFrame->addPlotable(residualHistogram, "P");
367370
mSgnPdf->plotOn(mResidualFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineColor(kBlue));
368371
}
369-
mass->setRange("bkgForSignificance", mRooMeanSgn->getVal() - mNSigmaForSgn * mRooSigmaSgn->getVal(), mRooMeanSgn->getVal() + mNSigmaForSgn * mRooSigmaSgn->getVal());
372+
mass->setRange("bkgForSignificance", mRooMeanSgn->getVal() - mNSigmaForSgn * mRooSecSigmaSgn->getVal(), mRooMeanSgn->getVal() + mNSigmaForSgn * mRooSecSigmaSgn->getVal());
370373
bkgIntegral = mBkgPdf->createIntegral(*mass, NormSet(*mass), Range("bkgForSignificance"));
371374
mIntegralBkg = bkgIntegral->getValV();
372375
calculateBackground(mBkgYield, mBkgYieldErr);
@@ -376,6 +379,8 @@ void HFInvMassFitter::doFit()
376379
calculateSignal(mRawYield, mRawYieldErr);
377380
countSignal(mRawYieldCounted, mRawYieldCountedErr);
378381
calculateSignificance(mSignificance, mSignificanceErr);
382+
// Fit to data ratio
383+
mRatioFrame = mass->frame(Title("Fit/Data Ratio"));
379384
calculateFitToDataRatio();
380385
}
381386
}
@@ -445,10 +450,10 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const
445450
workspace.import(*sgnFuncGaus);
446451
delete sgnFuncGaus;
447452
// signal double Gaussian
448-
RooRealVar sigmaDoubleGaus("sigmaDoubleGaus", "sigma2Gaus", mSigmaSgn, mSigmaSgn - 0.01, mSigmaSgn + 0.01);
453+
RooRealVar sigmaDoubleGaus("sigmaDoubleGaus", "sigma2Gaus", mSigmaSgnDoubleGaus, mSigmaSgnDoubleGaus - 0.003, mSigmaSgnDoubleGaus + 0.003);
449454
if (mBoundSigma) {
450-
sigmaDoubleGaus.setMax(mSigmaSgn * (1 + mParamSgn));
451-
sigmaDoubleGaus.setMin(mSigmaSgn * (1 - mParamSgn));
455+
sigmaDoubleGaus.setMax(mSigmaSgnDoubleGaus * (1 + mParamSgn));
456+
sigmaDoubleGaus.setMin(mSigmaSgnDoubleGaus * (1 - mParamSgn));
452457
}
453458
if (mFixedSigma) {
454459
sigma.setVal(mSigmaSgn);
@@ -557,6 +562,7 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const
557562
workspace.import(*reflFuncPoly6);
558563
delete reflFuncPoly6;
559564
}
565+
560566
// draw fit output
561567
void HFInvMassFitter::drawFit(TVirtualPad* pad, Int_t writeFitInfo)
562568
{
@@ -591,13 +597,15 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, Int_t writeFitInfo)
591597
textInfoRight->AddText(Form("mean(free) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
592598
}
593599
if (mTypeOfSgnPdf == DoubleGaus) {
594-
auto const& baseSigmaSgn = mWorkspace->var("sigma");
600+
if (mFixedSigma) {
601+
textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
602+
} else {
603+
textInfoRight->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
604+
}
595605
if (mFixedSigmaDoubleGaus) {
596-
textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
597-
textInfoRight->AddText(Form("sigma 2(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
606+
textInfoRight->AddText(Form("sigma 2(fixed) = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError()));
598607
} else {
599-
textInfoRight->AddText(Form("sigma(free) = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
600-
textInfoRight->AddText(Form("sigma 2(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
608+
textInfoRight->AddText(Form("sigma 2(free) = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError()));
601609
}
602610
} else if (mFixedSigma) {
603611
textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
@@ -631,9 +639,8 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad)
631639
textInfo->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr));
632640
textInfo->AddText(Form("mean = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
633641
if (mTypeOfSgnPdf == DoubleGaus) {
634-
auto const& baseSigmaSgn = mWorkspace->var("sigma");
635-
textInfo->AddText(Form("sigma = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
636-
textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
642+
textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
643+
textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError()));
637644
} else {
638645
textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
639646
}
@@ -655,9 +662,8 @@ void HFInvMassFitter::drawRatio(TVirtualPad* pad)
655662
textInfo->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr));
656663
textInfo->AddText(Form("mean = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
657664
if (mTypeOfSgnPdf == DoubleGaus) {
658-
auto const& baseSigmaSgn = mWorkspace->var("sigma");
659-
textInfo->AddText(Form("sigma = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
660-
textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
665+
textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
666+
textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError()));
661667
} else {
662668
textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
663669
}
@@ -681,7 +687,7 @@ void HFInvMassFitter::highlightPeakRegion(const RooPlot* plot, Color_t color, Wi
681687
double yMin = plot->GetMinimum();
682688
double yMax = plot->GetMaximum();
683689
const Double_t mean = mRooMeanSgn->getVal();
684-
const Double_t sigma = mRooSigmaSgn->getVal();
690+
const Double_t sigma = mRooSecSigmaSgn->getVal();
685691
const Double_t minForSgn = mean - mNSigmaForSidebands * sigma;
686692
const Double_t maxForSgn = mean + mNSigmaForSidebands * sigma;
687693
TLine* leftLine = new TLine(minForSgn, yMin, minForSgn, yMax);
@@ -706,7 +712,7 @@ void HFInvMassFitter::drawReflection(TVirtualPad* pad)
706712
void HFInvMassFitter::countSignal(Double_t& signal, Double_t& signalErr) const
707713
{
708714
const Double_t mean = mRooMeanSgn->getVal();
709-
const Double_t sigma = mRooSigmaSgn->getVal();
715+
const Double_t sigma = mRooSecSigmaSgn->getVal();
710716
const Double_t minForSgn = mean - mNSigmaForSidebands * sigma;
711717
const Double_t maxForSgn = mean + mNSigmaForSidebands * sigma;
712718
const Int_t binForMinSgn = mHistoInvMass->FindBin(minForSgn);
@@ -798,21 +804,25 @@ RooAbsPdf* HFInvMassFitter::createSignalFitFunction(RooWorkspace* workspace)
798804
case 0: {
799805
sgnPdf = workspace->pdf("sgnFuncGaus");
800806
mRooSigmaSgn = workspace->var("sigma");
807+
mRooSecSigmaSgn = workspace->var("sigma");
801808
mRooMeanSgn = workspace->var("mean");
802809
} break;
803810
case 1: {
804811
sgnPdf = workspace->pdf("sgnFuncDoubleGaus");
805-
mRooSigmaSgn = workspace->var("sigmaDoubleGaus");
812+
mRooSigmaSgn = workspace->var("sigma");
813+
mRooSecSigmaSgn = workspace->var("sigmaDoubleGaus");
806814
mRooMeanSgn = workspace->var("mean");
807815
} break;
808816
case 2: {
809817
sgnPdf = workspace->pdf("sgnFuncGausRatio");
810-
mRooSigmaSgn = workspace->var("sigmaDoubleGausRatio");
818+
mRooSigmaSgn = workspace->var("sigma");
819+
mRooSecSigmaSgn = workspace->var("sigmaDoubleGausRatio");
811820
mRooMeanSgn = workspace->var("mean");
812821
} break;
813822
case 3: {
814823
sgnPdf = workspace->pdf("sgnFuncDoublePeak");
815-
mRooSigmaSgn = workspace->var("sigmaSec");
824+
mRooSigmaSgn = workspace->var("sigma");
825+
mRooSecSigmaSgn = workspace->var("sigmaSec");
816826
mRooMeanSgn = workspace->var("meanSec");
817827
} break;
818828
default:
@@ -858,8 +868,6 @@ void HFInvMassFitter::calculateFitToDataRatio()
858868
{
859869
if (!mInvMassFrame)
860870
return;
861-
// Create a new RooPlot for the ratio
862-
mRatioFrame = mWorkspace->var("mass")->frame(Title("Fit/Data Ratio"));
863871

864872
// Get the data and fit curves from the frame
865873
auto* dataHist = dynamic_cast<RooHist*>(mInvMassFrame->findObject("data_c"));
@@ -884,8 +892,8 @@ void HFInvMassFitter::calculateFitToDataRatio()
884892
ratioHist->SetPointError(i, 0, 0, err, err);
885893
}
886894

887-
mRatioFrame->SetMinimum(0.0);
888-
mRatioFrame->SetMaximum(2.0);
895+
mRatioFrame->SetMinimum(0.5);
896+
mRatioFrame->SetMaximum(1.5);
889897
mRatioFrame->addPlotable(ratioHist, "P");
890898
}
891899

‎PWGHF/D2H/Macros/HFInvMassFitter.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -211,6 +211,8 @@ class HFInvMassFitter : public TNamed
211211
Double_t getMeanUncertainty() const { return mRooMeanSgn->getError(); }
212212
Double_t getSigma() const { return mRooSigmaSgn->getVal(); }
213213
Double_t getSigmaUncertainty() const { return mRooSigmaSgn->getError(); }
214+
Double_t getSecSigma() const { return mRooSecSigmaSgn->getVal(); }
215+
Double_t getSecSigmaUncertainty() const { return mRooSecSigmaSgn->getError(); }
214216
Double_t getReflOverSig() const
215217
{
216218
if (mReflPdf) {
@@ -284,6 +286,7 @@ class HFInvMassFitter : public TNamed
284286
Bool_t mFixReflOverSgn; /// switch for fix refl/signal
285287
RooRealVar* mRooMeanSgn; /// mean for gaussian of signal
286288
RooRealVar* mRooSigmaSgn; /// sigma for gaussian of signal
289+
RooRealVar* mRooSecSigmaSgn; /// second sigma for composite gaussian of signal
287290
RooAbsPdf* mSgnPdf; /// signal fit function
288291
RooAbsPdf* mBkgPdf; /// background fit function
289292
RooAbsPdf* mReflPdf; /// reflection fit function

‎PWGHF/D2H/Macros/runMassFitter.C‎

Lines changed: 31 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -262,6 +262,9 @@ int runMassFitter(const TString& configFileName)
262262
auto hRawYieldsSigma = new TH1D(
263263
"hRawYieldsSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})",
264264
nSliceVarBins, sliceVarLimits.data());
265+
auto hRawYieldsSecSigma = new TH1D(
266+
"hRawYieldsSecSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})",
267+
nSliceVarBins, sliceVarLimits.data());
265268
auto hRawYieldsMean = new TH1D(
266269
"hRawYieldsMean", ";" + sliceVarName + "(" + sliceVarUnit + ");mean (GeV/#it{c}^{2})",
267270
nSliceVarBins, sliceVarLimits.data());
@@ -298,6 +301,7 @@ int runMassFitter(const TString& configFileName)
298301
setHistoStyle(hRawYieldsSignal);
299302
setHistoStyle(hRawYieldsSignalCounted);
300303
setHistoStyle(hRawYieldsSigma);
304+
setHistoStyle(hRawYieldsSecSigma);
301305
setHistoStyle(hRawYieldsMean);
302306
setHistoStyle(hRawYieldsSignificance);
303307
setHistoStyle(hRawYieldsSgnOverBkg);
@@ -417,6 +421,8 @@ int runMassFitter(const TString& configFileName)
417421

418422
const Double_t sigma = massFitter->getSigma();
419423
const Double_t sigmaErr = massFitter->getSigmaUncertainty();
424+
const Double_t secSigma = massFitter->getSecSigma();
425+
const Double_t secSigmaErr = massFitter->getSecSigmaUncertainty();
420426
const Double_t mean = massFitter->getMean();
421427
const Double_t meanErr = massFitter->getMeanUncertainty();
422428
const Double_t reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg();
@@ -428,6 +434,8 @@ int runMassFitter(const TString& configFileName)
428434
hRawYieldsSignalCounted->SetBinError(iSliceVar + 1, rawYieldCountedErr);
429435
hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma);
430436
hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr);
437+
hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma);
438+
hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr);
431439
hRawYieldsMean->SetBinContent(iSliceVar + 1, mean);
432440
hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr);
433441
hRawYieldsChiSquareBkg->SetBinContent(iSliceVar + 1, reducedChiSquareBkg);
@@ -482,6 +490,8 @@ int runMassFitter(const TString& configFileName)
482490
const double rawYieldCountedErr = massFitter->getRawYieldCountedError();
483491
const double sigma = massFitter->getSigma();
484492
const double sigmaErr = massFitter->getSigmaUncertainty();
493+
const double secSigma = massFitter->getSecSigma();
494+
const double secSigmaErr = massFitter->getSecSigmaUncertainty();
485495
const double mean = massFitter->getMean();
486496
const double meanErr = massFitter->getMeanUncertainty();
487497
const double reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg();
@@ -497,6 +507,8 @@ int runMassFitter(const TString& configFileName)
497507
hRawYieldsSignalCounted->SetBinError(iSliceVar + 1, rawYieldCountedErr);
498508
hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma);
499509
hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr);
510+
hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma);
511+
hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr);
500512
hRawYieldsMean->SetBinContent(iSliceVar + 1, mean);
501513
hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr);
502514
hRawYieldsSignificance->SetBinContent(iSliceVar + 1, significance);
@@ -533,23 +545,25 @@ int runMassFitter(const TString& configFileName)
533545
canvasMass[iCanvas]->Modified();
534546
canvasMass[iCanvas]->Update();
535547

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

545-
if (nSliceVarBins > 1) {
546-
canvasRatio[iCanvas]->cd(iSliceVar - nCanvasesMax * iCanvas + 1);
547-
} else {
548-
canvasRatio[iCanvas]->cd();
558+
if (nSliceVarBins > 1) {
559+
canvasRatio[iCanvas]->cd(iSliceVar - nCanvasesMax * iCanvas + 1);
560+
} else {
561+
canvasRatio[iCanvas]->cd();
562+
}
563+
massFitter->drawRatio(gPad);
564+
canvasRatio[iCanvas]->Modified();
565+
canvasRatio[iCanvas]->Update();
549566
}
550-
massFitter->drawRatio(gPad);
551-
canvasRatio[iCanvas]->Modified();
552-
canvasRatio[iCanvas]->Update();
553567
}
554568

555569
hFitConfig->SetBinContent(1, iSliceVar + 1, massMin[iSliceVar]);
@@ -585,6 +599,7 @@ int runMassFitter(const TString& configFileName)
585599
hRawYieldsSignal->Write();
586600
hRawYieldsSignalCounted->Write();
587601
hRawYieldsSigma->Write();
602+
hRawYieldsSecSigma->Write();
588603
hRawYieldsMean->Write();
589604
hRawYieldsSignificance->Write();
590605
hRawYieldsSgnOverBkg->Write();
@@ -635,7 +650,7 @@ void setHistoStyle(TH1* histo, Color_t color, Size_t markerSize)
635650
void divideCanvas(TCanvas* canvas, int nSliceVarBins)
636651
{
637652
int nCols = std::ceil(std::sqrt(nSliceVarBins));
638-
int nRows = std::ceil(double(nSliceVarBins) / nCols);
653+
int nRows = std::ceil(static_cast<double>(nSliceVarBins) / nCols);
639654
canvas->Divide(nCols, nRows);
640655
}
641656

0 commit comments

Comments
 (0)