Skip to content

Commit b464a70

Browse files
pstahlhupstahlhu
committed
Store fraction for double Gaussian and add possibility to fix it
1 parent 5488306 commit b464a70

File tree

3 files changed

+55
-26
lines changed

3 files changed

+55
-26
lines changed

‎PWGHF/D2H/Macros/HFInvMassFitter.cxx‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -104,6 +104,7 @@ HFInvMassFitter::HFInvMassFitter() : mHistoInvMass(nullptr),
104104
mRooMeanSgn(nullptr),
105105
mRooSigmaSgn(nullptr),
106106
mRooSecSigmaSgn(nullptr),
107+
mRooFracDoubleGaus(nullptr),
107108
mSgnPdf(nullptr),
108109
mBkgPdf(nullptr),
109110
mReflPdf(nullptr),
@@ -176,6 +177,7 @@ HFInvMassFitter::HFInvMassFitter(const TH1* histoToFit, Double_t minValue, Doubl
176177
mRooMeanSgn(nullptr),
177178
mRooSigmaSgn(nullptr),
178179
mRooSecSigmaSgn(nullptr),
180+
mRooFracDoubleGaus(nullptr),
179181
mSgnPdf(nullptr),
180182
mBkgPdf(nullptr),
181183
mReflPdf(nullptr),
@@ -210,6 +212,7 @@ HFInvMassFitter::~HFInvMassFitter()
210212
delete mRooMeanSgn;
211213
delete mRooSigmaSgn;
212214
delete mRooSecSigmaSgn;
215+
delete mRooFracDoubleGaus;
213216
delete mSgnPdf;
214217
delete mBkgPdf;
215218
delete mReflPdf;
@@ -810,12 +813,14 @@ RooAbsPdf* HFInvMassFitter::createSignalFitFunction(RooWorkspace* workspace)
810813
mRooSigmaSgn = workspace->var("sigma");
811814
mRooSecSigmaSgn = workspace->var("sigmaDoubleGaus");
812815
mRooMeanSgn = workspace->var("mean");
816+
mRooFracDoubleGaus = workspace->var("fracDoubleGaus");
813817
} break;
814818
case 2: {
815819
sgnPdf = workspace->pdf("sgnFuncGausRatio");
816820
mRooSigmaSgn = workspace->var("sigma");
817821
mRooSecSigmaSgn = workspace->var("sigmaDoubleGausRatio");
818822
mRooMeanSgn = workspace->var("mean");
823+
mRooFracDoubleGaus = workspace->var("fracDoubleGausRatio");
819824
} break;
820825
case 3: {
821826
sgnPdf = workspace->pdf("sgnFuncDoublePeak");

‎PWGHF/D2H/Macros/HFInvMassFitter.h‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -214,8 +214,10 @@ class HFInvMassFitter : public TNamed
214214
[[nodiscard]] Double_t getSigmaUncertainty() const { return mRooSigmaSgn->getError(); }
215215
[[nodiscard]] Double_t getSecSigma() const { return mRooSecSigmaSgn->getVal(); }
216216
[[nodiscard]] Double_t getSecSigmaUncertainty() const { return mRooSecSigmaSgn->getError(); }
217+
[[nodiscard]] Double_t getFracDoubleGaus() const { return mRooFracDoubleGaus->getVal(); }
218+
[[nodiscard]] Double_t getFracDoubleGausUncertainty() const { return mRooFracDoubleGaus->getError(); }
217219
[[nodiscard]] Double_t getReflOverSig() const
218-
220+
219221
{
220222
if (mReflPdf != nullptr) {
221223
return mReflOverSgn;
@@ -288,6 +290,7 @@ class HFInvMassFitter : public TNamed
288290
RooRealVar* mRooMeanSgn; /// mean for gaussian of signal
289291
RooRealVar* mRooSigmaSgn; /// sigma for gaussian of signal
290292
RooRealVar* mRooSecSigmaSgn; /// second sigma for composite gaussian of signal
293+
RooRealVar* mRooFracDoubleGaus; /// fraction of second gaussian for composite gaussian of signal
291294
RooAbsPdf* mSgnPdf; /// signal fit function
292295
RooAbsPdf* mBkgPdf; /// background fit function
293296
RooAbsPdf* mReflPdf; /// reflection fit function

‎PWGHF/D2H/Macros/runMassFitter.C‎

Lines changed: 46 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -105,9 +105,10 @@ int runMassFitter(const TString& configFileName)
105105
std::vector<double> sliceVarMax;
106106
std::vector<double> massMin;
107107
std::vector<double> massMax;
108+
std::vector<double> fixMeanManual;
108109
std::vector<double> fixSigmaManual;
109110
std::vector<double> fixSecondSigmaManual;
110-
std::vector<double> fixMeanManual;
111+
std::vector<double> fixFracDoubleGausManual;
111112
std::vector<int> nRebin;
112113
std::vector<int> bkgFuncConfig;
113114
std::vector<int> sgnFuncConfig;
@@ -130,12 +131,15 @@ int runMassFitter(const TString& configFileName)
130131
const Value& fdSecPeakHistoNameValue = config["FDSecPeakHistoName"];
131132
parseStringArray(fdSecPeakHistoNameValue, fdSecPeakHistoName);
132133

133-
const bool fixSigma = config["FixSigma"].GetBool();
134-
const std::string sigmaFile = config["SigmaFile"].GetString();
135-
136134
const bool fixMean = config["FixMean"].GetBool();
137135
const std::string meanFile = config["MeanFile"].GetString();
138136

137+
const Value& fixMeanManualValue = config["FixMeanManual"];
138+
readArray(fixMeanManualValue, fixMeanManual);
139+
140+
const bool fixSigma = config["FixSigma"].GetBool();
141+
const std::string sigmaFile = config["SigmaFile"].GetString();
142+
139143
const Value& fixSigmaManualValue = config["FixSigmaManual"];
140144
readArray(fixSigmaManualValue, fixSigmaManual);
141145

@@ -145,8 +149,11 @@ int runMassFitter(const TString& configFileName)
145149
const Value& fixSecondSigmaManualValue = config["FixSecondSigmaManual"];
146150
readArray(fixSecondSigmaManualValue, fixSecondSigmaManual);
147151

148-
const Value& fixMeanManualValue = config["FixMeanManual"];
149-
readArray(fixMeanManualValue, fixMeanManual);
152+
const bool fixFracDoubleGaus = config["FixFracDoubleGaus"].GetBool();
153+
const std::string fracDoubleGausFile = config["FracDoubleGausFile"].GetString();
154+
155+
const Value& fixFracDoubleGausManualValue = config["FixFracDoubleGausManual"];
156+
readArray(fixFracDoubleGausManualValue, fixFracDoubleGausManual);
150157

151158
sliceVarName = config["SliceVarName"].GetString();
152159
sliceVarUnit = config["SliceVarUnit"].GetString();
@@ -262,18 +269,21 @@ int runMassFitter(const TString& configFileName)
262269
inputFile->Close();
263270

264271
// define output histos
265-
auto* hRawYieldsSignal = new TH1D("hRawYieldsSignal", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield",
266-
nSliceVarBins, sliceVarLimits.data());
267-
auto* hRawYieldsSignalCounted = new TH1D("hRawYieldsSignalCounted", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield via bin count",
268-
nSliceVarBins, sliceVarLimits.data());
269-
auto* hRawYieldsSigma = new TH1D(
272+
auto hRawYieldsSignal = new TH1D("hRawYieldsSignal", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield",
273+
nSliceVarBins, sliceVarLimits.data());
274+
auto hRawYieldsSignalCounted = new TH1D("hRawYieldsSignalCounted", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield via bin count",
275+
nSliceVarBins, sliceVarLimits.data());
276+
auto hRawYieldsMean = new TH1D(
277+
"hRawYieldsMean", ";" + sliceVarName + "(" + sliceVarUnit + ");mean (GeV/#it{c}^{2})",
278+
nSliceVarBins, sliceVarLimits.data());
279+
auto hRawYieldsSigma = new TH1D(
270280
"hRawYieldsSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})",
271281
nSliceVarBins, sliceVarLimits.data());
272282
auto hRawYieldsSecSigma = new TH1D(
273283
"hRawYieldsSecSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})",
274284
nSliceVarBins, sliceVarLimits.data());
275-
auto hRawYieldsMean = new TH1D(
276-
"hRawYieldsMean", ";" + sliceVarName + "(" + sliceVarUnit + ");mean (GeV/#it{c}^{2})",
285+
auto hRawYieldsFracDoubleGaus = new TH1D(
286+
"hRawYieldsFracDoubleGaus", ";" + sliceVarName + "(" + sliceVarUnit + ");fraction of double gaussian",
277287
nSliceVarBins, sliceVarLimits.data());
278288
auto* hRawYieldsSignificance = new TH1D(
279289
"hRawYieldsSignificance",
@@ -307,9 +317,10 @@ int runMassFitter(const TString& configFileName)
307317

308318
setHistoStyle(hRawYieldsSignal);
309319
setHistoStyle(hRawYieldsSignalCounted);
320+
setHistoStyle(hRawYieldsMean);
310321
setHistoStyle(hRawYieldsSigma);
311322
setHistoStyle(hRawYieldsSecSigma);
312-
setHistoStyle(hRawYieldsMean);
323+
setHistoStyle(hRawYieldsFracDoubleGaus);
313324
setHistoStyle(hRawYieldsSignificance);
314325
setHistoStyle(hRawYieldsSgnOverBkg);
315326
setHistoStyle(hRawYieldsBkg);
@@ -339,7 +350,8 @@ int runMassFitter(const TString& configFileName)
339350

340351
TH1* hSigmaToFix = getHistToFix(fixSigma, fixSigmaManual, sigmaFile, "Sigma");
341352
TH1* hMeanToFix = getHistToFix(fixMean, fixMeanManual, meanFile, "Mean");
342-
TH1* hSecondSigmaToFix = getHistToFix(fixSecondSigma, fixSecondSigmaManual, secondSigmaFile, "Sigma");
353+
TH1* hSecondSigmaToFix = getHistToFix(fixSecondSigma, fixSecondSigmaManual, secondSigmaFile, "SecSigma");
354+
TH1* hFracDoubleGausToFix = getHistToFix(fixFracDoubleGaus, fixFracDoubleGausManual, fracDoubleGausFile, "FracDoubleGaus");
343355

344356
// fit histograms
345357

@@ -420,26 +432,29 @@ int runMassFitter(const TString& configFileName)
420432
const Double_t rawYieldErr = massFitter->getRawYieldError();
421433
const Double_t rawYieldCounted = massFitter->getRawYieldCounted();
422434
const Double_t rawYieldCountedErr = massFitter->getRawYieldCountedError();
423-
435+
const Double_t mean = massFitter->getMean();
436+
const Double_t meanErr = massFitter->getMeanUncertainty();
424437
const Double_t sigma = massFitter->getSigma();
425438
const Double_t sigmaErr = massFitter->getSigmaUncertainty();
426439
const Double_t secSigma = massFitter->getSecSigma();
427440
const Double_t secSigmaErr = massFitter->getSecSigmaUncertainty();
428-
const Double_t mean = massFitter->getMean();
429-
const Double_t meanErr = massFitter->getMeanUncertainty();
441+
const Double_t fracDoubleGaus = massFitter->getFracDoubleGaus();
442+
const Double_t fracDoubleGausErr = massFitter->getFracDoubleGausUncertainty();
430443
const Double_t reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg();
431444
const Double_t reducedChiSquareTotal = massFitter->getChiSquareOverNDFTotal();
432445

433446
hRawYieldsSignal->SetBinContent(iSliceVar + 1, rawYield);
434447
hRawYieldsSignal->SetBinError(iSliceVar + 1, rawYieldErr);
435448
hRawYieldsSignalCounted->SetBinContent(iSliceVar + 1, rawYieldCounted);
436449
hRawYieldsSignalCounted->SetBinError(iSliceVar + 1, rawYieldCountedErr);
450+
hRawYieldsMean->SetBinContent(iSliceVar + 1, mean);
451+
hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr);
437452
hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma);
438453
hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr);
439454
hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma);
440455
hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr);
441-
hRawYieldsMean->SetBinContent(iSliceVar + 1, mean);
442-
hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr);
456+
hRawYieldsFracDoubleGaus->SetBinContent(iSliceVar + 1, fracDoubleGaus);
457+
hRawYieldsFracDoubleGaus->SetBinError(iSliceVar + 1, fracDoubleGausErr);
443458
hRawYieldsChiSquareBkg->SetBinContent(iSliceVar + 1, reducedChiSquareBkg);
444459
hRawYieldsChiSquareBkg->SetBinError(iSliceVar + 1, 0.);
445460
hRawYieldsChiSquareTotal->SetBinContent(iSliceVar + 1, reducedChiSquareTotal);
@@ -476,6 +491,7 @@ int runMassFitter(const TString& configFileName)
476491
setFixedValue(fixMean, fixMeanManual, hMeanToFix, std::bind(&HFInvMassFitter::setFixGaussianMean, massFitter, std::placeholders::_1), "MEAN");
477492
setFixedValue(fixSigma, fixSigmaManual, hSigmaToFix, std::bind(&HFInvMassFitter::setFixGaussianSigma, massFitter, std::placeholders::_1), "SIGMA");
478493
setFixedValue(fixSecondSigma, fixSecondSigmaManual, hSecondSigmaToFix, std::bind(&HFInvMassFitter::setFixSecondGaussianSigma, massFitter, std::placeholders::_1), "SECOND SIGMA");
494+
setFixedValue(fixFracDoubleGaus, fixFracDoubleGausManual, nullptr, std::bind(&HFInvMassFitter::setFixFrac2Gaus, massFitter, std::placeholders::_1), "FRAC DOUBLE GAUS");
479495

480496
if (enableRefl) {
481497
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)
490506
const double rawYieldErr = massFitter->getRawYieldError();
491507
const double rawYieldCounted = massFitter->getRawYieldCounted();
492508
const double rawYieldCountedErr = massFitter->getRawYieldCountedError();
509+
const double mean = massFitter->getMean();
510+
const double meanErr = massFitter->getMeanUncertainty();
493511
const double sigma = massFitter->getSigma();
494512
const double sigmaErr = massFitter->getSigmaUncertainty();
495513
const double secSigma = massFitter->getSecSigma();
496514
const double secSigmaErr = massFitter->getSecSigmaUncertainty();
497-
const double mean = massFitter->getMean();
498-
const double meanErr = massFitter->getMeanUncertainty();
515+
const double fracDoubleGaus = massFitter->getFracDoubleGaus();
516+
const double fracDoubleGausErr = massFitter->getFracDoubleGausUncertainty();
499517
const double reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg();
500518
const double reducedChiSquareTotal = massFitter->getChiSquareOverNDFTotal();
501519
const double significance = massFitter->getSignificance();
@@ -507,12 +525,14 @@ int runMassFitter(const TString& configFileName)
507525
hRawYieldsSignal->SetBinError(iSliceVar + 1, rawYieldErr);
508526
hRawYieldsSignalCounted->SetBinContent(iSliceVar + 1, rawYieldCounted);
509527
hRawYieldsSignalCounted->SetBinError(iSliceVar + 1, rawYieldCountedErr);
528+
hRawYieldsMean->SetBinContent(iSliceVar + 1, mean);
529+
hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr);
510530
hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma);
511531
hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr);
512532
hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma);
513533
hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr);
514-
hRawYieldsMean->SetBinContent(iSliceVar + 1, mean);
515-
hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr);
534+
hRawYieldsFracDoubleGaus->SetBinContent(iSliceVar + 1, fracDoubleGaus);
535+
hRawYieldsFracDoubleGaus->SetBinError(iSliceVar + 1, fracDoubleGausErr);
516536
hRawYieldsSignificance->SetBinContent(iSliceVar + 1, significance);
517537
hRawYieldsSignificance->SetBinError(iSliceVar + 1, significanceErr);
518538
hRawYieldsSgnOverBkg->SetBinContent(iSliceVar + 1, rawYield / bkg);
@@ -597,9 +617,10 @@ int runMassFitter(const TString& configFileName)
597617
}
598618
hRawYieldsSignal->Write();
599619
hRawYieldsSignalCounted->Write();
620+
hRawYieldsMean->Write();
600621
hRawYieldsSigma->Write();
601622
hRawYieldsSecSigma->Write();
602-
hRawYieldsMean->Write();
623+
hRawYieldsFracDoubleGaus->Write();
603624
hRawYieldsSignificance->Write();
604625
hRawYieldsSgnOverBkg->Write();
605626
hRawYieldsBkg->Write();

0 commit comments

Comments
 (0)