Skip to content

Commit 0feee38

Browse files
lubynetslubynets
committed
add dscb i/o in runMassFitter
1 parent 683f0f1 commit 0feee38

1 file changed

Lines changed: 123 additions & 0 deletions

File tree

‎PWGHF/D2H/Macros/runMassFitter.C‎

Lines changed: 123 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -66,6 +66,8 @@ T readJsonField(const Document& config, const std::string& fieldName);
6666
template <typename T>
6767
void readJsonVector(std::vector<T>& vec, const Document& config, const std::string& fieldName, bool isRequired = false);
6868

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

7173
void setHistoStyle(TH1* histo, Color_t color = kBlack, Size_t markerSize = 1);
@@ -109,6 +111,18 @@ void runMassFitter(const std::string& configFileName)
109111
std::vector<int> nRebin;
110112
std::vector<int> bkgFunc;
111113
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;
112126

113127
readJsonVector(inputHistoName, config, "InputHistoName", true);
114128
readJsonVector(promptHistoName, config, "PromptHistoName");
@@ -133,6 +147,8 @@ void runMassFitter(const std::string& configFileName)
133147
const std::string fracDoubleGausFile = readJsonField<std::string>(config, "FracDoubleGausFile", "");
134148
readJsonVector(fixFracDoubleGausManual, config, "FixFracDoubleGausManual");
135149

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

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

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+
158199
const int nSliceVarBins = static_cast<int>(sliceVarMin.size());
159200
std::vector<double> sliceVarLimits(nSliceVarBins + 1);
160201

@@ -184,6 +225,18 @@ void runMassFitter(const std::string& configFileName)
184225
checkVectorSize(nRebin, "nRebin");
185226
checkVectorSize(bkgFunc, "bkgFunc");
186227
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);
187240

188241
for (int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) {
189242
sliceVarLimits[iSliceVar] = sliceVarMin[iSliceVar];
@@ -268,6 +321,10 @@ void runMassFitter(const std::string& configFileName)
268321
auto* hRawYieldsSigma = new TH1D("hRawYieldsSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data());
269322
auto* hRawYieldsSecSigma = new TH1D("hRawYieldsSecSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data());
270323
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());
271328

272329
enum {
273330
ConfigMassMin = 1,
@@ -300,6 +357,10 @@ void runMassFitter(const std::string& configFileName)
300357
setHistoStyle(hRawYieldsSigma);
301358
setHistoStyle(hRawYieldsSecSigma);
302359
setHistoStyle(hRawYieldsFracDoubleGaus);
360+
setHistoStyle(hRawYieldsDscbAlphaL);
361+
setHistoStyle(hRawYieldsDscbAlphaR);
362+
setHistoStyle(hRawYieldsDscbNL);
363+
setHistoStyle(hRawYieldsDscbNR);
303364

304365
auto getHistToFix = [&nSliceVarBins](bool const& isFix, std::vector<double> const& fixManual, std::string const& fixFileName, std::string const& var) -> TH1* {
305366
TH1* histToFix = nullptr;
@@ -400,6 +461,10 @@ void runMassFitter(const std::string& configFileName)
400461
setFixedValue(fixSigma, fixSigmaManual, hSigmaToFix, std::bind(&HFInvMassFitter::setFixGaussianSigma, massFitter, std::placeholders::_1), "SIGMA");
401462
setFixedValue(fixSecondSigma, fixSecondSigmaManual, hSecondSigmaToFix, std::bind(&HFInvMassFitter::setFixSecondGaussianSigma, massFitter, std::placeholders::_1), "SECOND SIGMA");
402463
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");
403468

404469
if (!isMc && enableRefl) {
405470
reflOverSgn = hMassSgn[iSliceVar]->Integral(hMassSgn[iSliceVar]->FindBin(massMin[iSliceVar] * 1.0001), hMassSgn[iSliceVar]->FindBin(massMax[iSliceVar] * 0.999));
@@ -408,6 +473,24 @@ void runMassFitter(const std::string& configFileName)
408473
massFitter->setTemplateReflections(hMassRefl[iSliceVar]);
409474
}
410475

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+
411494
massFitter->doFit();
412495

413496
auto drawOnCanvas = [&](std::vector<TCanvas*>& canvas, std::function<void()> drawer) {
@@ -444,6 +527,14 @@ void runMassFitter(const std::string& configFileName)
444527
const double meanErr = massFitter->getMeanUncertainty();
445528
const double sigma = massFitter->getSigma();
446529
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();
447538

448539
hRawYieldsSignal->SetBinContent(iSliceVar + 1, rawYield);
449540
hRawYieldsSignal->SetBinError(iSliceVar + 1, rawYieldErr);
@@ -464,6 +555,14 @@ void runMassFitter(const std::string& configFileName)
464555
hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma);
465556
hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr);
466557
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);
467566

468567
if (sgnFunc[iSliceVar] != HFInvMassFitter::SingleGaus) { // TODO foresee DSCB and Voigt cases
469568
const double secSigma = massFitter->getSecSigma();
@@ -519,6 +618,12 @@ void runMassFitter(const std::string& configFileName)
519618
if (enableRefl) {
520619
hReflectionOverSignal->Write();
521620
}
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+
}
522627
hFitConfig->Write();
523628

524629
outputFile.Close();
@@ -631,6 +736,24 @@ void readJsonVector(std::vector<T>& vec, const Document& config, const std::stri
631736
}
632737
}
633738

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+
634757
int main(int argc, const char* argv[])
635758
{
636759
if (argc == 1) {

0 commit comments

Comments
 (0)