Skip to content

Commit 302ecee

Browse files
pstahlhupstahlhu
committed
Store fraction for double Gaussian and add possibility to fix it
1 parent 2fcfa1a commit 302ecee

File tree

3 files changed

+49
-20
lines changed

3 files changed

+49
-20
lines changed

‎PWGHF/D2H/Macros/HFInvMassFitter.cxx‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -106,6 +106,7 @@ HFInvMassFitter::HFInvMassFitter() : TNamed(),
106106
mRooMeanSgn(nullptr),
107107
mRooSigmaSgn(nullptr),
108108
mRooSecSigmaSgn(nullptr),
109+
mRooFracDoubleGaus(nullptr),
109110
mSgnPdf(nullptr),
110111
mBkgPdf(nullptr),
111112
mReflPdf(nullptr),
@@ -179,6 +180,7 @@ HFInvMassFitter::HFInvMassFitter(const TH1* histoToFit, Double_t minValue, Doubl
179180
mRooMeanSgn(nullptr),
180181
mRooSigmaSgn(nullptr),
181182
mRooSecSigmaSgn(nullptr),
183+
mRooFracDoubleGaus(nullptr),
182184
mSgnPdf(nullptr),
183185
mBkgPdf(nullptr),
184186
mReflPdf(nullptr),
@@ -213,6 +215,7 @@ HFInvMassFitter::~HFInvMassFitter()
213215
delete mRooMeanSgn;
214216
delete mRooSigmaSgn;
215217
delete mRooSecSigmaSgn;
218+
delete mRooFracDoubleGaus;
216219
delete mSgnPdf;
217220
delete mBkgPdf;
218221
delete mReflPdf;
@@ -812,12 +815,14 @@ RooAbsPdf* HFInvMassFitter::createSignalFitFunction(RooWorkspace* workspace)
812815
mRooSigmaSgn = workspace->var("sigma");
813816
mRooSecSigmaSgn = workspace->var("sigmaDoubleGaus");
814817
mRooMeanSgn = workspace->var("mean");
818+
mRooFracDoubleGaus = workspace->var("fracDoubleGaus");
815819
} break;
816820
case 2: {
817821
sgnPdf = workspace->pdf("sgnFuncGausRatio");
818822
mRooSigmaSgn = workspace->var("sigma");
819823
mRooSecSigmaSgn = workspace->var("sigmaDoubleGausRatio");
820824
mRooMeanSgn = workspace->var("mean");
825+
mRooFracDoubleGaus = workspace->var("fracDoubleGausRatio");
821826
} break;
822827
case 3: {
823828
sgnPdf = workspace->pdf("sgnFuncDoublePeak");

‎PWGHF/D2H/Macros/HFInvMassFitter.h‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -213,6 +213,8 @@ class HFInvMassFitter : public TNamed
213213
Double_t getSigmaUncertainty() const { return mRooSigmaSgn->getError(); }
214214
Double_t getSecSigma() const { return mRooSecSigmaSgn->getVal(); }
215215
Double_t getSecSigmaUncertainty() const { return mRooSecSigmaSgn->getError(); }
216+
Double_t getFracDoubleGaus() const { return mRooFracDoubleGaus->getVal(); }
217+
Double_t getFracDoubleGausUncertainty() const { return mRooFracDoubleGaus->getError(); }
216218
Double_t getReflOverSig() const
217219
{
218220
if (mReflPdf) {
@@ -287,6 +289,7 @@ class HFInvMassFitter : public TNamed
287289
RooRealVar* mRooMeanSgn; /// mean for gaussian of signal
288290
RooRealVar* mRooSigmaSgn; /// sigma for gaussian of signal
289291
RooRealVar* mRooSecSigmaSgn; /// second sigma for composite gaussian of signal
292+
RooRealVar* mRooFracDoubleGaus; /// fraction of second gaussian for composite gaussian of signal
290293
RooAbsPdf* mSgnPdf; /// signal fit function
291294
RooAbsPdf* mBkgPdf; /// background fit function
292295
RooAbsPdf* mReflPdf; /// reflection fit function

‎PWGHF/D2H/Macros/runMassFitter.C‎

Lines changed: 41 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -98,9 +98,10 @@ int runMassFitter(const TString& configFileName)
9898
std::vector<double> sliceVarMax;
9999
std::vector<double> massMin;
100100
std::vector<double> massMax;
101+
std::vector<double> fixMeanManual;
101102
std::vector<double> fixSigmaManual;
102103
std::vector<double> fixSecondSigmaManual;
103-
std::vector<double> fixMeanManual;
104+
std::vector<double> fixFracDoubleGausManual;
104105
std::vector<int> nRebin;
105106
std::vector<int> bkgFuncConfig;
106107
std::vector<int> sgnFuncConfig;
@@ -123,12 +124,15 @@ int runMassFitter(const TString& configFileName)
123124
const Value& fdSecPeakHistoNameValue = config["FDSecPeakHistoName"];
124125
parseStringArray(fdSecPeakHistoNameValue, fdSecPeakHistoName);
125126

126-
const bool fixSigma = config["FixSigma"].GetBool();
127-
const std::string sigmaFile = config["SigmaFile"].GetString();
128-
129127
const bool fixMean = config["FixMean"].GetBool();
130128
const std::string meanFile = config["MeanFile"].GetString();
131129

130+
const Value& fixMeanManualValue = config["FixMeanManual"];
131+
readArray(fixMeanManualValue, fixMeanManual);
132+
133+
const bool fixSigma = config["FixSigma"].GetBool();
134+
const std::string sigmaFile = config["SigmaFile"].GetString();
135+
132136
const Value& fixSigmaManualValue = config["FixSigmaManual"];
133137
readArray(fixSigmaManualValue, fixSigmaManual);
134138

@@ -138,8 +142,11 @@ int runMassFitter(const TString& configFileName)
138142
const Value& fixSecondSigmaManualValue = config["FixSecondSigmaManual"];
139143
readArray(fixSecondSigmaManualValue, fixSecondSigmaManual);
140144

141-
const Value& fixMeanManualValue = config["FixMeanManual"];
142-
readArray(fixMeanManualValue, fixMeanManual);
145+
const bool fixFracDoubleGaus = config["FixFracDoubleGaus"].GetBool();
146+
const std::string fracDoubleGausFile = config["FracDoubleGausFile"].GetString();
147+
148+
const Value& fixFracDoubleGausManualValue = config["FixFracDoubleGausManual"];
149+
readArray(fixFracDoubleGausManualValue, fixFracDoubleGausManual);
143150

144151
sliceVarName = config["SliceVarName"].GetString();
145152
sliceVarUnit = config["SliceVarUnit"].GetString();
@@ -259,14 +266,17 @@ int runMassFitter(const TString& configFileName)
259266
nSliceVarBins, sliceVarLimits.data());
260267
auto hRawYieldsSignalCounted = new TH1D("hRawYieldsSignalCounted", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield via bin count",
261268
nSliceVarBins, sliceVarLimits.data());
269+
auto hRawYieldsMean = new TH1D(
270+
"hRawYieldsMean", ";" + sliceVarName + "(" + sliceVarUnit + ");mean (GeV/#it{c}^{2})",
271+
nSliceVarBins, sliceVarLimits.data());
262272
auto hRawYieldsSigma = new TH1D(
263273
"hRawYieldsSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})",
264274
nSliceVarBins, sliceVarLimits.data());
265275
auto hRawYieldsSecSigma = new TH1D(
266276
"hRawYieldsSecSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})",
267277
nSliceVarBins, sliceVarLimits.data());
268-
auto hRawYieldsMean = new TH1D(
269-
"hRawYieldsMean", ";" + sliceVarName + "(" + sliceVarUnit + ");mean (GeV/#it{c}^{2})",
278+
auto hRawYieldsFracDoubleGaus = new TH1D(
279+
"hRawYieldsFracDoubleGaus", ";" + sliceVarName + "(" + sliceVarUnit + ");fraction of double gaussian",
270280
nSliceVarBins, sliceVarLimits.data());
271281
auto hRawYieldsSignificance = new TH1D(
272282
"hRawYieldsSignificance",
@@ -300,9 +310,10 @@ int runMassFitter(const TString& configFileName)
300310

301311
setHistoStyle(hRawYieldsSignal);
302312
setHistoStyle(hRawYieldsSignalCounted);
313+
setHistoStyle(hRawYieldsMean);
303314
setHistoStyle(hRawYieldsSigma);
304315
setHistoStyle(hRawYieldsSecSigma);
305-
setHistoStyle(hRawYieldsMean);
316+
setHistoStyle(hRawYieldsFracDoubleGaus);
306317
setHistoStyle(hRawYieldsSignificance);
307318
setHistoStyle(hRawYieldsSgnOverBkg);
308319
setHistoStyle(hRawYieldsBkg);
@@ -332,7 +343,8 @@ int runMassFitter(const TString& configFileName)
332343

333344
TH1* hSigmaToFix = getHistToFix(fixSigma, fixSigmaManual, sigmaFile, "Sigma");
334345
TH1* hMeanToFix = getHistToFix(fixMean, fixMeanManual, meanFile, "Mean");
335-
TH1* hSecondSigmaToFix = getHistToFix(fixSecondSigma, fixSecondSigmaManual, secondSigmaFile, "Sigma");
346+
TH1* hSecondSigmaToFix = getHistToFix(fixSecondSigma, fixSecondSigmaManual, secondSigmaFile, "SecSigma");
347+
TH1* hFracDoubleGausToFix = getHistToFix(fixFracDoubleGaus, fixFracDoubleGausManual, fracDoubleGausFile, "FracDoubleGaus");
336348

337349
// fit histograms
338350

@@ -418,26 +430,29 @@ int runMassFitter(const TString& configFileName)
418430
const Double_t rawYieldErr = massFitter->getRawYieldError();
419431
const Double_t rawYieldCounted = massFitter->getRawYieldCounted();
420432
const Double_t rawYieldCountedErr = massFitter->getRawYieldCountedError();
421-
433+
const Double_t mean = massFitter->getMean();
434+
const Double_t meanErr = massFitter->getMeanUncertainty();
422435
const Double_t sigma = massFitter->getSigma();
423436
const Double_t sigmaErr = massFitter->getSigmaUncertainty();
424437
const Double_t secSigma = massFitter->getSecSigma();
425438
const Double_t secSigmaErr = massFitter->getSecSigmaUncertainty();
426-
const Double_t mean = massFitter->getMean();
427-
const Double_t meanErr = massFitter->getMeanUncertainty();
439+
const Double_t fracDoubleGaus = massFitter->getFracDoubleGaus();
440+
const Double_t fracDoubleGausErr = massFitter->getFracDoubleGausUncertainty();
428441
const Double_t reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg();
429442
const Double_t reducedChiSquareTotal = massFitter->getChiSquareOverNDFTotal();
430443

431444
hRawYieldsSignal->SetBinContent(iSliceVar + 1, rawYield);
432445
hRawYieldsSignal->SetBinError(iSliceVar + 1, rawYieldErr);
433446
hRawYieldsSignalCounted->SetBinContent(iSliceVar + 1, rawYieldCounted);
434447
hRawYieldsSignalCounted->SetBinError(iSliceVar + 1, rawYieldCountedErr);
448+
hRawYieldsMean->SetBinContent(iSliceVar + 1, mean);
449+
hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr);
435450
hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma);
436451
hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr);
437452
hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma);
438453
hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr);
439-
hRawYieldsMean->SetBinContent(iSliceVar + 1, mean);
440-
hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr);
454+
hRawYieldsFracDoubleGaus->SetBinContent(iSliceVar + 1, fracDoubleGaus);
455+
hRawYieldsFracDoubleGaus->SetBinError(iSliceVar + 1, fracDoubleGausErr);
441456
hRawYieldsChiSquareBkg->SetBinContent(iSliceVar + 1, reducedChiSquareBkg);
442457
hRawYieldsChiSquareBkg->SetBinError(iSliceVar + 1, 0.);
443458
hRawYieldsChiSquareTotal->SetBinContent(iSliceVar + 1, reducedChiSquareTotal);
@@ -474,6 +489,7 @@ int runMassFitter(const TString& configFileName)
474489
setFixedValue(fixMean, fixMeanManual, hMeanToFix, std::bind(&HFInvMassFitter::setFixGaussianMean, massFitter, std::placeholders::_1), "MEAN");
475490
setFixedValue(fixSigma, fixSigmaManual, hSigmaToFix, std::bind(&HFInvMassFitter::setFixGaussianSigma, massFitter, std::placeholders::_1), "SIGMA");
476491
setFixedValue(fixSecondSigma, fixSecondSigmaManual, hSecondSigmaToFix, std::bind(&HFInvMassFitter::setFixSecondGaussianSigma, massFitter, std::placeholders::_1), "SECOND SIGMA");
492+
setFixedValue(fixFracDoubleGaus, fixFracDoubleGausManual, nullptr, std::bind(&HFInvMassFitter::setFixFrac2Gaus, massFitter, std::placeholders::_1), "FRAC DOUBLE GAUS");
477493

478494
if (enableRefl) {
479495
reflOverSgn = hMassForSgn[iSliceVar]->Integral(hMassForSgn[iSliceVar]->FindBin(massMin[iSliceVar] * 1.0001), hMassForSgn[iSliceVar]->FindBin(massMax[iSliceVar] * 0.999));
@@ -488,12 +504,14 @@ int runMassFitter(const TString& configFileName)
488504
const double rawYieldErr = massFitter->getRawYieldError();
489505
const double rawYieldCounted = massFitter->getRawYieldCounted();
490506
const double rawYieldCountedErr = massFitter->getRawYieldCountedError();
507+
const double mean = massFitter->getMean();
508+
const double meanErr = massFitter->getMeanUncertainty();
491509
const double sigma = massFitter->getSigma();
492510
const double sigmaErr = massFitter->getSigmaUncertainty();
493511
const double secSigma = massFitter->getSecSigma();
494512
const double secSigmaErr = massFitter->getSecSigmaUncertainty();
495-
const double mean = massFitter->getMean();
496-
const double meanErr = massFitter->getMeanUncertainty();
513+
const double fracDoubleGaus = massFitter->getFracDoubleGaus();
514+
const double fracDoubleGausErr = massFitter->getFracDoubleGausUncertainty();
497515
const double reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg();
498516
const double reducedChiSquareTotal = massFitter->getChiSquareOverNDFTotal();
499517
const double significance = massFitter->getSignificance();
@@ -505,12 +523,14 @@ int runMassFitter(const TString& configFileName)
505523
hRawYieldsSignal->SetBinError(iSliceVar + 1, rawYieldErr);
506524
hRawYieldsSignalCounted->SetBinContent(iSliceVar + 1, rawYieldCounted);
507525
hRawYieldsSignalCounted->SetBinError(iSliceVar + 1, rawYieldCountedErr);
526+
hRawYieldsMean->SetBinContent(iSliceVar + 1, mean);
527+
hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr);
508528
hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma);
509529
hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr);
510530
hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma);
511531
hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr);
512-
hRawYieldsMean->SetBinContent(iSliceVar + 1, mean);
513-
hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr);
532+
hRawYieldsFracDoubleGaus->SetBinContent(iSliceVar + 1, fracDoubleGaus);
533+
hRawYieldsFracDoubleGaus->SetBinError(iSliceVar + 1, fracDoubleGausErr);
514534
hRawYieldsSignificance->SetBinContent(iSliceVar + 1, significance);
515535
hRawYieldsSignificance->SetBinError(iSliceVar + 1, significanceErr);
516536
hRawYieldsSgnOverBkg->SetBinContent(iSliceVar + 1, rawYield / bkg);
@@ -598,9 +618,10 @@ int runMassFitter(const TString& configFileName)
598618
}
599619
hRawYieldsSignal->Write();
600620
hRawYieldsSignalCounted->Write();
621+
hRawYieldsMean->Write();
601622
hRawYieldsSigma->Write();
602623
hRawYieldsSecSigma->Write();
603-
hRawYieldsMean->Write();
624+
hRawYieldsFracDoubleGaus->Write();
604625
hRawYieldsSignificance->Write();
605626
hRawYieldsSgnOverBkg->Write();
606627
hRawYieldsBkg->Write();

0 commit comments

Comments
 (0)