Skip to content

Commit ba58e6f

Browse files
saganattsaganatt
authored
[PWGHF] D2H fitter: Add possibility to set fixed manual mean and second sigma for double gauss (#12060)
1 parent 31b13b3 commit ba58e6f

File tree

3 files changed

+90
-49
lines changed

3 files changed

+90
-49
lines changed

‎PWGHF/D2H/Macros/HFInvMassFitter.cxx‎

Lines changed: 17 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -588,7 +588,16 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, Int_t writeFitInfo)
588588
} else {
589589
textInfoRight->AddText(Form("mean(free) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
590590
}
591-
if (mFixedSigma) {
591+
if (mTypeOfSgnPdf == DoubleGaus) {
592+
auto const& baseSigmaSgn = mWorkspace->var("sigma");
593+
if (mFixedSigmaDoubleGaus) {
594+
textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
595+
textInfoRight->AddText(Form("sigma 2(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
596+
} else {
597+
textInfoRight->AddText(Form("sigma(free) = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
598+
textInfoRight->AddText(Form("sigma 2(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
599+
}
600+
} else if (mFixedSigma) {
592601
textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
593602
} else {
594603
textInfoRight->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
@@ -619,7 +628,13 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad)
619628
textInfo->AddText(Form("S = %.0f #pm %.0f ", mRawYield, mRawYieldErr));
620629
textInfo->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr));
621630
textInfo->AddText(Form("mean = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError()));
622-
textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
631+
if (mTypeOfSgnPdf == DoubleGaus) {
632+
auto const& baseSigmaSgn = mWorkspace->var("sigma");
633+
textInfo->AddText(Form("sigma = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError()));
634+
textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
635+
} else {
636+
textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError()));
637+
}
623638
mResidualFrame->addObject(textInfo);
624639
mResidualFrame->Draw();
625640
highlightPeakRegion(mResidualFrame);

‎PWGHF/D2H/Macros/config_massfitter.json‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,25 @@
5050
"FixMean": false,
5151
"MeanFile": "",
5252
"_MeanFile": "fix mean from file",
53+
"FixMeanManual": [
54+
0,
55+
0,
56+
0,
57+
0,
58+
0
59+
],
60+
"_FixMeanManual": "fix mean mannually",
61+
"FixSecondSigma": false,
62+
"SecondSigmaFile": "",
63+
"_SecondSigmaFile": "fix second sigma for double gauss from file",
64+
"FixSecondSigmaManual": [
65+
0,
66+
0,
67+
0,
68+
0,
69+
0
70+
],
71+
"_FixSecondSigmaManual": "fix second sigma for double gauss manually",
5372
"SliceVarName": "T",
5473
"SliceVarUnit": "ps",
5574
"_SliceVarName, _SliceVarUnit": "e.g. pT, GeV/c or something else depending on user's needs",

‎PWGHF/D2H/Macros/runMassFitter.C‎

Lines changed: 54 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -99,6 +99,8 @@ int runMassFitter(const TString& configFileName)
9999
std::vector<double> massMin;
100100
std::vector<double> massMax;
101101
std::vector<double> fixSigmaManual;
102+
std::vector<double> fixSecondSigmaManual;
103+
std::vector<double> fixMeanManual;
102104
std::vector<int> nRebin;
103105
std::vector<int> bkgFuncConfig;
104106
std::vector<int> sgnFuncConfig;
@@ -121,15 +123,24 @@ int runMassFitter(const TString& configFileName)
121123
const Value& fdSecPeakHistoNameValue = config["FDSecPeakHistoName"];
122124
parseStringArray(fdSecPeakHistoNameValue, fdSecPeakHistoName);
123125

124-
bool fixSigma = config["FixSigma"].GetBool();
125-
std::string sigmaFile = config["SigmaFile"].GetString();
126+
const bool fixSigma = config["FixSigma"].GetBool();
127+
const std::string sigmaFile = config["SigmaFile"].GetString();
126128

127-
bool fixMean = config["FixMean"].GetBool();
128-
std::string meanFile = config["MeanFile"].GetString();
129+
const bool fixMean = config["FixMean"].GetBool();
130+
const std::string meanFile = config["MeanFile"].GetString();
129131

130132
const Value& fixSigmaManualValue = config["FixSigmaManual"];
131133
readArray(fixSigmaManualValue, fixSigmaManual);
132134

135+
const bool fixSecondSigma = config["FixSecondSigma"].GetBool();
136+
const std::string secondSigmaFile = config["SecondSigmaFile"].GetString();
137+
138+
const Value& fixSecondSigmaManualValue = config["FixSecondSigmaManual"];
139+
readArray(fixSecondSigmaManualValue, fixSecondSigmaManual);
140+
141+
const Value& fixMeanManualValue = config["FixMeanManual"];
142+
readArray(fixMeanManualValue, fixMeanManual);
143+
133144
sliceVarName = config["SliceVarName"].GetString();
134145
sliceVarUnit = config["SliceVarUnit"].GetString();
135146

@@ -295,35 +306,29 @@ int runMassFitter(const TString& configFileName)
295306
setHistoStyle(hRawYieldsChiSquareTotal);
296307
setHistoStyle(hReflectionOverSignal, kRed + 1);
297308

298-
TH1* hSigmaToFix = nullptr;
299-
if (fixSigma) {
300-
if (fixSigmaManual.empty()) {
301-
auto inputFileSigma = TFile::Open(sigmaFile.data());
302-
if (!inputFileSigma) {
303-
return -2;
304-
}
305-
hSigmaToFix = inputFileSigma->Get<TH1>("hRawYieldsSigma");
306-
hSigmaToFix->SetDirectory(0);
307-
if (static_cast<unsigned int>(hSigmaToFix->GetNbinsX()) != nSliceVarBins) {
308-
printf("WARNING: Different number of bins for this analysis and histo for fix sigma!\n");
309+
auto getHistToFix = [&nSliceVarBins](bool const& isFix, std::vector<double> const& fixManual, std::string const& fixFileName, std::string const& var) -> TH1* {
310+
TH1* histToFix = nullptr;
311+
if (isFix) {
312+
if (fixManual.empty()) {
313+
auto fixInputFile = TFile::Open(fixFileName.data());
314+
if (!fixInputFile) {
315+
throw std::runtime_error("Cannot open file for fixed " + var);
316+
}
317+
const std::string histName = "hRawYields" + var;
318+
histToFix = fixInputFile->Get<TH1>(histName.data());
319+
histToFix->SetDirectory(nullptr);
320+
if (static_cast<unsigned int>(histToFix->GetNbinsX()) != nSliceVarBins) {
321+
throw std::runtime_error("Different number of bins for this analysis and histo for fixed " + var);
322+
}
323+
fixInputFile->Close();
309324
}
310-
inputFileSigma->Close();
311325
}
312-
}
326+
return histToFix;
327+
};
313328

314-
TH1* hMeanToFix = nullptr;
315-
if (fixMean) {
316-
auto inputFileMean = TFile::Open(meanFile.data());
317-
if (!inputFileMean) {
318-
return -3;
319-
}
320-
hMeanToFix = inputFileMean->Get<TH1>("hRawYieldsMean");
321-
hMeanToFix->SetDirectory(0);
322-
if (static_cast<unsigned int>(hMeanToFix->GetNbinsX()) != nSliceVarBins) {
323-
printf("WARNING: Different number of bins for this analysis and histo for fix mean\n");
324-
}
325-
inputFileMean->Close();
326-
}
329+
TH1* hSigmaToFix = getHistToFix(fixSigma, fixSigmaManual, sigmaFile, "Sigma");
330+
TH1* hMeanToFix = getHistToFix(fixMean, fixMeanManual, meanFile, "Mean");
331+
TH1* hSecondSigmaToFix = getHistToFix(fixSecondSigma, fixSecondSigmaManual, secondSigmaFile, "Sigma");
327332

328333
// fit histograms
329334

@@ -433,24 +438,26 @@ int runMassFitter(const TString& configFileName)
433438
if (useLikelihood) {
434439
massFitter->setUseLikelihoodFit();
435440
}
436-
if (fixMean) {
437-
massFitter->setFixGaussianMean(hMeanToFix->GetBinContent(iSliceVar + 1));
438-
}
439-
if (fixSigma) {
440-
if (fixSigmaManual.empty()) {
441-
massFitter->setFixGaussianSigma(hSigmaToFix->GetBinContent(iSliceVar + 1));
442-
printf("*****************************\n");
443-
printf("FIXED SIGMA: %f\n", hSigmaToFix->GetBinContent(iSliceVar + 1));
444-
printf("*****************************\n");
445-
} else if (!fixSigmaManual.empty()) {
446-
massFitter->setFixGaussianSigma(fixSigmaManual[iSliceVar]);
447-
printf("*****************************\n");
448-
printf("FIXED SIGMA: %f\n", fixSigmaManual[iSliceVar]);
449-
printf("*****************************\n");
450-
} else {
451-
printf("WARNING: impossible to fix sigma! Wrong fix sigma file or value!\n");
441+
442+
auto setFixedValue = [&massFitter, &iSliceVar](bool const& isFix, std::vector<double> const& fixManual, const TH1* histToFix, std::function<void(Double_t)> setFunc, std::string const& var) -> void {
443+
if (isFix) {
444+
if (fixManual.empty()) {
445+
setFunc(histToFix->GetBinContent(iSliceVar + 1));
446+
printf("*****************************\n");
447+
printf("FIXED %s: %f\n", var.data(), histToFix->GetBinContent(iSliceVar + 1));
448+
printf("*****************************\n");
449+
} else {
450+
setFunc(fixManual[iSliceVar]);
451+
printf("*****************************\n");
452+
printf("FIXED %s: %f\n", var.data(), fixManual[iSliceVar]);
453+
printf("*****************************\n");
454+
}
452455
}
453-
}
456+
};
457+
458+
setFixedValue(fixMean, fixMeanManual, hMeanToFix, std::bind(&HFInvMassFitter::setFixGaussianMean, massFitter, std::placeholders::_1), "MEAN");
459+
setFixedValue(fixSigma, fixSigmaManual, hSigmaToFix, std::bind(&HFInvMassFitter::setFixGaussianSigma, massFitter, std::placeholders::_1), "SIGMA");
460+
setFixedValue(fixSecondSigma, fixSecondSigmaManual, hSecondSigmaToFix, std::bind(&HFInvMassFitter::setFixSecondGaussianSigma, massFitter, std::placeholders::_1), "SECOND SIGMA");
454461

455462
if (enableRefl) {
456463
reflOverSgn = hMassForSgn[iSliceVar]->Integral(hMassForSgn[iSliceVar]->FindBin(massMin[iSliceVar] * 1.0001), hMassForSgn[iSliceVar]->FindBin(massMax[iSliceVar] * 0.999));

0 commit comments

Comments
 (0)