Skip to content

Commit cc8fbbd

Browse files
lubynetslubynets
committed
add dscb i/o in runMassFitter
1 parent be6c3e3 commit cc8fbbd

1 file changed

Lines changed: 123 additions & 1 deletion

File tree

‎PWGHF/D2H/Macros/runMassFitter.C‎

Lines changed: 123 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,6 @@
4343
#include <map>
4444
#include <stdexcept>
4545
#include <string>
46-
#include <tuple>
4746
#include <type_traits>
4847
#include <vector>
4948

@@ -67,6 +66,8 @@ T readJsonField(const Document& config, const std::string& fieldName);
6766
template <typename T>
6867
void readJsonVector(std::vector<T>& vec, const Document& config, const std::string& fieldName, bool isRequired = false);
6968

69+
void readJsonVectorFromHisto(std::vector<double>& vec, const Document& config, const std::string& fileNameFieldName, const std::string& histoNameFieldName);
70+
7071
void divideCanvas(TCanvas* c, int nSliceVarBins);
7172

7273
void setHistoStyle(TH1* histo, Color_t color = kBlack, Size_t markerSize = 1);
@@ -110,6 +111,18 @@ void runMassFitter(const std::string& configFileName)
110111
std::vector<int> nRebin;
111112
std::vector<int> bkgFunc;
112113
std::vector<int> sgnFunc;
114+
std::vector<double> dscbAlphaLInitial;
115+
std::vector<double> dscbAlphaLLower;
116+
std::vector<double> dscbAlphaLUpper;
117+
std::vector<double> dscbAlphaRInitial;
118+
std::vector<double> dscbAlphaRLower;
119+
std::vector<double> dscbAlphaRUpper;
120+
std::vector<double> dscbNLInitial;
121+
std::vector<double> dscbNLLower;
122+
std::vector<double> dscbNLUpper;
123+
std::vector<double> dscbNRInitial;
124+
std::vector<double> dscbNRLower;
125+
std::vector<double> dscbNRUpper;
113126

114127
readJsonVector(inputHistoName, config, "InputHistoName", true);
115128
readJsonVector(promptHistoName, config, "PromptHistoName");
@@ -134,6 +147,8 @@ void runMassFitter(const std::string& configFileName)
134147
const std::string fracDoubleGausFile = readJsonField<std::string>(config, "FracDoubleGausFile", "");
135148
readJsonVector(fixFracDoubleGausManual, config, "FixFracDoubleGausManual");
136149

150+
const bool fixDscbTailParams = readJsonField<bool>(config, "FixDscbTailParams", false);
151+
137152
const TString sliceVarName = readJsonField<std::string>(config, "SliceVarName");
138153
const TString sliceVarUnit = readJsonField<std::string>(config, "SliceVarUnit");
139154

@@ -156,6 +171,31 @@ void runMassFitter(const std::string& configFileName)
156171
const bool highlightPeakRegion = readJsonField<bool>(config, "HighlightPeakRegion", true);
157172
const int randomSeed = readJsonField<int>(config, "RandomSeed", -1);
158173

174+
readJsonVector(dscbAlphaLInitial, config, "DscbAlphaLInitial");
175+
readJsonVector(dscbAlphaLLower, config, "DscbAlphaLLower");
176+
readJsonVector(dscbAlphaLUpper, config, "DscbAlphaLUpper");
177+
readJsonVector(dscbAlphaRInitial, config, "DscbAlphaRInitial");
178+
readJsonVector(dscbAlphaRLower, config, "DscbAlphaRLower");
179+
readJsonVector(dscbAlphaRUpper, config, "DscbAlphaRUpper");
180+
readJsonVector(dscbNLInitial, config, "DscbNLInitial");
181+
readJsonVector(dscbNLLower, config, "DscbNLLower");
182+
readJsonVector(dscbNLUpper, config, "DscbNLUpper");
183+
readJsonVector(dscbNRInitial, config, "DscbNRInitial");
184+
readJsonVector(dscbNRLower, config, "DscbNRLower");
185+
readJsonVector(dscbNRUpper, config, "DscbNRUpper");
186+
readJsonVectorFromHisto(dscbAlphaLInitial, config, "DscbParametersFile", "DscbAlphaLInitialHisto");
187+
readJsonVectorFromHisto(dscbAlphaLLower, config, "DscbParametersFile", "DscbAlphaLLowerHisto");
188+
readJsonVectorFromHisto(dscbAlphaLUpper, config, "DscbParametersFile", "DscbAlphaLUpperHisto");
189+
readJsonVectorFromHisto(dscbAlphaRInitial, config, "DscbParametersFile", "DscbAlphaRInitialHisto");
190+
readJsonVectorFromHisto(dscbAlphaRLower, config, "DscbParametersFile", "DscbAlphaRLowerHisto");
191+
readJsonVectorFromHisto(dscbAlphaRUpper, config, "DscbParametersFile", "DscbAlphaRUpperHisto");
192+
readJsonVectorFromHisto(dscbNLInitial, config, "DscbParametersFile", "DscbNLInitialHisto");
193+
readJsonVectorFromHisto(dscbNLLower, config, "DscbParametersFile", "DscbNLLowerHisto");
194+
readJsonVectorFromHisto(dscbNLUpper, config, "DscbParametersFile", "DscbNLUpperHisto");
195+
readJsonVectorFromHisto(dscbNRInitial, config, "DscbParametersFile", "DscbNRInitialHisto");
196+
readJsonVectorFromHisto(dscbNRLower, config, "DscbParametersFile", "DscbNRLowerHisto");
197+
readJsonVectorFromHisto(dscbNRUpper, config, "DscbParametersFile", "DscbNRUpperHisto");
198+
159199
const int nSliceVarBins = static_cast<int>(sliceVarMin.size());
160200
std::vector<double> sliceVarLimits(nSliceVarBins + 1);
161201

@@ -185,6 +225,18 @@ void runMassFitter(const std::string& configFileName)
185225
checkVectorSize(nRebin, "nRebin");
186226
checkVectorSize(bkgFunc, "bkgFunc");
187227
checkVectorSize(sgnFunc, "sgnFunc");
228+
checkVectorSize(dscbAlphaLInitial, "dscbAlphaLInitial", true);
229+
checkVectorSize(dscbAlphaLLower, "dscbAlphaLLower", true);
230+
checkVectorSize(dscbAlphaLUpper, "dscbAlphaLUpper", true);
231+
checkVectorSize(dscbAlphaRInitial, "dscbAlphaRInitial", true);
232+
checkVectorSize(dscbAlphaRLower, "dscbAlphaRLower", true);
233+
checkVectorSize(dscbAlphaRUpper, "dscbAlphaRUpper", true);
234+
checkVectorSize(dscbNLInitial, "dscbNLInitial", true);
235+
checkVectorSize(dscbNLLower, "dscbNLLower", true);
236+
checkVectorSize(dscbNLUpper, "dscbNLUpper", true);
237+
checkVectorSize(dscbNRInitial, "dscbNRInitial", true);
238+
checkVectorSize(dscbNRLower, "dscbNRLower", true);
239+
checkVectorSize(dscbNRUpper, "dscbNRUpper", true);
188240

189241
for (int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) {
190242
sliceVarLimits[iSliceVar] = sliceVarMin[iSliceVar];
@@ -269,6 +321,10 @@ void runMassFitter(const std::string& configFileName)
269321
auto* hRawYieldsSigma = new TH1D("hRawYieldsSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data());
270322
auto* hRawYieldsSecSigma = new TH1D("hRawYieldsSecSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data());
271323
auto* hRawYieldsFracDoubleGaus = new TH1D("hRawYieldsFracDoubleGaus", ";" + sliceVarName + "(" + sliceVarUnit + ");fraction of double gaussian", nSliceVarBins, sliceVarLimits.data());
324+
auto* hRawYieldsDscbAlphaL = new TH1D("hRawYieldsDscbAlphaL", ";" + sliceVarName + "(" + sliceVarUnit + ");#alpha_{L}", nSliceVarBins, sliceVarLimits.data());
325+
auto* hRawYieldsDscbAlphaR = new TH1D("hRawYieldsDscbAlphaR", ";" + sliceVarName + "(" + sliceVarUnit + ");#alpha_{R}", nSliceVarBins, sliceVarLimits.data());
326+
auto* hRawYieldsDscbNL = new TH1D("hRawYieldsDscbNL", ";" + sliceVarName + "(" + sliceVarUnit + ");n_{L}", nSliceVarBins, sliceVarLimits.data());
327+
auto* hRawYieldsDscbNR = new TH1D("hRawYieldsDscbNR", ";" + sliceVarName + "(" + sliceVarUnit + ");n_{R}", nSliceVarBins, sliceVarLimits.data());
272328

273329
enum {
274330
ConfigMassMin = 1,
@@ -301,6 +357,10 @@ void runMassFitter(const std::string& configFileName)
301357
setHistoStyle(hRawYieldsSigma);
302358
setHistoStyle(hRawYieldsSecSigma);
303359
setHistoStyle(hRawYieldsFracDoubleGaus);
360+
setHistoStyle(hRawYieldsDscbAlphaL);
361+
setHistoStyle(hRawYieldsDscbAlphaR);
362+
setHistoStyle(hRawYieldsDscbNL);
363+
setHistoStyle(hRawYieldsDscbNR);
304364

305365
auto getHistToFix = [&nSliceVarBins](bool const& isFix, std::vector<double> const& fixManual, std::string const& fixFileName, std::string const& var) -> TH1* {
306366
TH1* histToFix = nullptr;
@@ -401,6 +461,10 @@ void runMassFitter(const std::string& configFileName)
401461
setFixedValue(fixSigma, fixSigmaManual, hSigmaToFix, std::bind(&HFInvMassFitter::setFixGaussianSigma, massFitter, std::placeholders::_1), "SIGMA");
402462
setFixedValue(fixSecondSigma, fixSecondSigmaManual, hSecondSigmaToFix, std::bind(&HFInvMassFitter::setFixSecondGaussianSigma, massFitter, std::placeholders::_1), "SECOND SIGMA");
403463
setFixedValue(fixFracDoubleGaus, fixFracDoubleGausManual, hFracDoubleGausToFix, std::bind(&HFInvMassFitter::setFixFrac2Gaus, massFitter, std::placeholders::_1), "FRAC DOUBLE GAUS");
464+
setFixedValue(fixDscbTailParams, dscbAlphaLInitial, nullptr, std::bind(&HFInvMassFitter::setFixDscbAlphaL, massFitter, std::placeholders::_1), "DSCB ALPHA LEFT");
465+
setFixedValue(fixDscbTailParams, dscbAlphaRInitial, nullptr, std::bind(&HFInvMassFitter::setFixDscbAlphaR, massFitter, std::placeholders::_1), "DSCB ALPHA RIGHT");
466+
setFixedValue(fixDscbTailParams, dscbNLInitial, nullptr, std::bind(&HFInvMassFitter::setFixDscbNL, massFitter, std::placeholders::_1), "DSCB N LEFT");
467+
setFixedValue(fixDscbTailParams, dscbNRInitial, nullptr, std::bind(&HFInvMassFitter::setFixDscbNR, massFitter, std::placeholders::_1), "DSCB N RIGHT");
404468

405469
if (!isMc && enableRefl) {
406470
reflOverSgn = hMassSgn[iSliceVar]->Integral(hMassSgn[iSliceVar]->FindBin(massMin[iSliceVar] * 1.0001), hMassSgn[iSliceVar]->FindBin(massMax[iSliceVar] * 0.999));
@@ -409,6 +473,24 @@ void runMassFitter(const std::string& configFileName)
409473
massFitter->setTemplateReflections(hMassRefl[iSliceVar]);
410474
}
411475

476+
auto setDscbParameter = [&](const std::vector<double>& vec, void (HFInvMassFitter::*setter)(double)) {
477+
if (static_cast<int>(vec.size()) == nSliceVarBins) {
478+
(massFitter->*setter)(vec[iSliceVar]);
479+
}
480+
};
481+
setDscbParameter(dscbAlphaLInitial, &HFInvMassFitter::setDscbAlphaLInitialValue);
482+
setDscbParameter(dscbAlphaLLower, &HFInvMassFitter::setDscbAlphaLLowLimit);
483+
setDscbParameter(dscbAlphaLUpper, &HFInvMassFitter::setDscbAlphaLUpLimit);
484+
setDscbParameter(dscbAlphaRInitial, &HFInvMassFitter::setDscbAlphaRInitialValue);
485+
setDscbParameter(dscbAlphaRLower, &HFInvMassFitter::setDscbAlphaRLowLimit);
486+
setDscbParameter(dscbAlphaRUpper, &HFInvMassFitter::setDscbAlphaRUpLimit);
487+
setDscbParameter(dscbNLInitial, &HFInvMassFitter::setDscbNLInitialValue);
488+
setDscbParameter(dscbNLLower, &HFInvMassFitter::setDscbNLLowLimit);
489+
setDscbParameter(dscbNLUpper, &HFInvMassFitter::setDscbNLUpLimit);
490+
setDscbParameter(dscbNRInitial, &HFInvMassFitter::setDscbNRInitialValue);
491+
setDscbParameter(dscbNRLower, &HFInvMassFitter::setDscbNRLowLimit);
492+
setDscbParameter(dscbNRUpper, &HFInvMassFitter::setDscbNRUpLimit);
493+
412494
massFitter->doFit();
413495

414496
auto drawOnCanvas = [&](std::vector<TCanvas*>& canvas, std::function<void()> drawer) {
@@ -445,6 +527,14 @@ void runMassFitter(const std::string& configFileName)
445527
const double meanErr = massFitter->getMeanUncertainty();
446528
const double sigma = massFitter->getSigma();
447529
const double sigmaErr = massFitter->getSigmaUncertainty();
530+
const double dscbAlphaL = massFitter->getDscbAlphaL();
531+
const double dscbAlphaR = massFitter->getDscbAlphaR();
532+
const double dscbNL = massFitter->getDscbNL();
533+
const double dscbNR = massFitter->getDscbNR();
534+
const double dscbAlphaLErr = massFitter->getDscbAlphaLUncertainty();
535+
const double dscbAlphaRErr = massFitter->getDscbAlphaRUncertainty();
536+
const double dscbNErrL = massFitter->getDscbNLUncertainty();
537+
const double dscbNErrR = massFitter->getDscbNRUncertainty();
448538

449539
hRawYieldsSignal->SetBinContent(iSliceVar + 1, rawYield);
450540
hRawYieldsSignal->SetBinError(iSliceVar + 1, rawYieldErr);
@@ -465,6 +555,14 @@ void runMassFitter(const std::string& configFileName)
465555
hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma);
466556
hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr);
467557
hReflectionOverSignal->SetBinContent(iSliceVar + 1, reflOverSgn);
558+
hRawYieldsDscbAlphaL->SetBinContent(iSliceVar + 1, dscbAlphaL);
559+
hRawYieldsDscbAlphaL->SetBinError(iSliceVar + 1, dscbAlphaLErr);
560+
hRawYieldsDscbAlphaR->SetBinContent(iSliceVar + 1, dscbAlphaR);
561+
hRawYieldsDscbAlphaR->SetBinError(iSliceVar + 1, dscbAlphaRErr);
562+
hRawYieldsDscbNL->SetBinContent(iSliceVar + 1, dscbNL);
563+
hRawYieldsDscbNL->SetBinError(iSliceVar + 1, dscbNErrL);
564+
hRawYieldsDscbNR->SetBinContent(iSliceVar + 1, dscbNR);
565+
hRawYieldsDscbNR->SetBinError(iSliceVar + 1, dscbNErrR);
468566

469567
if (sgnFunc[iSliceVar] != HFInvMassFitter::SingleGaus) { // TODO foresee DSCB and Voigt cases
470568
const double secSigma = massFitter->getSecSigma();
@@ -520,6 +618,12 @@ void runMassFitter(const std::string& configFileName)
520618
if (enableRefl) {
521619
hReflectionOverSignal->Write();
522620
}
621+
if (std::find(sgnFunc.begin(), sgnFunc.end(), HFInvMassFitter::DoubleSidedCrystalBall) != sgnFunc.end()) {
622+
hRawYieldsDscbAlphaL->Write();
623+
hRawYieldsDscbAlphaR->Write();
624+
hRawYieldsDscbNL->Write();
625+
hRawYieldsDscbNR->Write();
626+
}
523627
hFitConfig->Write();
524628

525629
outputFile.Close();
@@ -632,6 +736,24 @@ void readJsonVector(std::vector<T>& vec, const Document& config, const std::stri
632736
}
633737
}
634738

739+
void readJsonVectorFromHisto(std::vector<double>& vec, const Document& config, const std::string& fileNameFieldName, const std::string& histoNameFieldName)
740+
{
741+
if (!vec.empty()) {
742+
throw std::runtime_error("readJsonVectorFromHisto(): vector is not empty!");
743+
}
744+
const auto fileName = readJsonField<std::string>(config, fileNameFieldName);
745+
const auto histoName = readJsonField<std::string>(config, histoNameFieldName);
746+
if (fileName.empty() || histoName.empty()) {
747+
return;
748+
}
749+
TFile* inputFile = openFileWithNullptrCheck(fileName);
750+
TH1* histo = getObjectWithNullPtrCheck<TH1>(inputFile, histoName);
751+
for (int iBin = 1; iBin <= histo->GetNbinsX(); iBin++) {
752+
vec.push_back(histo->GetBinContent(iBin));
753+
}
754+
inputFile->Close();
755+
}
756+
635757
int main(int argc, const char* argv[])
636758
{
637759
if (argc == 1) {

0 commit comments

Comments
 (0)