Skip to content

Commit cce147c

Browse files
committed
read vector from json: allow for a single value if all are equal
1 parent 053cdd8 commit cce147c

2 files changed

Lines changed: 88 additions & 53 deletions

File tree

PWGHF/D2H/Macros/config_massfitter.json

Lines changed: 10 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -101,21 +101,24 @@
101101
0.9,
102102
1.6
103103
],
104-
"MassMin": [
104+
"MassMin": 2.12,
105+
"_MassMin": [
105106
2.12,
106107
2.12,
107108
2.12,
108109
2.12,
109110
2.12
110111
],
111-
"MassMax": [
112+
"MassMax": 2.42,
113+
"_MassMax": [
112114
2.42,
113115
2.42,
114116
2.42,
115117
2.42,
116118
2.42
117119
],
118-
"Rebin": [
120+
"Rebin": 4,
121+
"_Rebin": [
119122
4,
120123
4,
121124
4,
@@ -128,7 +131,8 @@
128131
"true: likelihood fit",
129132
"false: chi2 fit"
130133
],
131-
"BkgFunc": [
134+
"BkgFunc": 2,
135+
"_BkgFunc": [
132136
2,
133137
2,
134138
2,
@@ -144,7 +148,8 @@
144148
"5 for Poly3",
145149
"6 for NoBkg"
146150
],
147-
"SgnFunc": [
151+
"SgnFunc": 0,
152+
"_SgnFunc": [
148153
0,
149154
0,
150155
0,

PWGHF/D2H/Macros/runMassFitter.C

Lines changed: 78 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,8 @@
5050

5151
using namespace rapidjson;
5252

53+
constexpr int UndefValueInt{-999};
54+
5355
void runMassFitter(const std::string& configFileName = "config_massfitter.json");
5456

5557
TFile* openFileWithNullptrCheck(const std::string& fileName, const std::string& option = "read");
@@ -66,9 +68,12 @@ T readJsonField(const Document& config, const std::string& fieldName);
6668
template <typename T>
6769
void readJsonVector(std::vector<T>& vec, const Document& config, const std::string& fieldName, bool isRequired = false);
6870

71+
template <typename T>
72+
void readJsonVectorFlexible(std::vector<T>& vec, const Document& config, int nHistograms, const std::string& fieldName, bool isRequired = false);
73+
6974
void readJsonVectorFromHisto(std::vector<double>& vec, const Document& config, const std::string& fileNameFieldName, const std::string& histoNameFieldName);
7075

71-
void divideCanvas(TCanvas* c, int nSliceVarBins);
76+
void divideCanvas(TCanvas* c, int nHistograms);
7277

7378
void setHistoStyle(TH1* histo, Color_t color = kBlack, Size_t markerSize = 1);
7479

@@ -127,6 +132,8 @@ void runMassFitter(const std::string& configFileName)
127132
std::vector<double> dscbNRUpper;
128133

129134
readJsonVector(inputHistoName, config, "InputHistoName", true);
135+
const int nHistograms = static_cast<int>(inputHistoName.size());
136+
130137
readJsonVector(promptHistoName, config, "PromptHistoName");
131138
readJsonVector(fdHistoName, config, "FDHistoName");
132139
readJsonVector(signalHistoName, config, "SignalHistoName");
@@ -137,19 +144,19 @@ void runMassFitter(const std::string& configFileName)
137144

138145
const bool fixMean = readJsonField<bool>(config, "FixMean", false);
139146
const std::string meanFile = readJsonField<std::string>(config, "MeanFile", "");
140-
readJsonVector(fixMeanManual, config, "FixMeanManual");
147+
readJsonVectorFlexible(fixMeanManual, config, nHistograms, "FixMeanManual");
141148

142149
const bool fixSigma = readJsonField<bool>(config, "FixSigma", false);
143150
const std::string sigmaFile = readJsonField<std::string>(config, "SigmaFile", "");
144-
readJsonVector(fixSigmaManual, config, "FixSigmaManual");
151+
readJsonVectorFlexible(fixSigmaManual, config, nHistograms, "FixSigmaManual");
145152

146153
const bool fixSecondSigma = readJsonField<bool>(config, "FixSecondSigma", false);
147154
const std::string secondSigmaFile = readJsonField<std::string>(config, "SecondSigmaFile", "");
148-
readJsonVector(fixSecondSigmaManual, config, "FixSecondSigmaManual");
155+
readJsonVectorFlexible(fixSecondSigmaManual, config, nHistograms, "FixSecondSigmaManual");
149156

150157
const bool fixFracDoubleGaus = readJsonField<bool>(config, "FixFracDoubleGaus", false);
151158
const std::string fracDoubleGausFile = readJsonField<std::string>(config, "FracDoubleGausFile", "");
152-
readJsonVector(fixFracDoubleGausManual, config, "FixFracDoubleGausManual");
159+
readJsonVectorFlexible(fixFracDoubleGausManual, config, nHistograms, "FixFracDoubleGausManual");
153160

154161
const bool fixDscbTailParams = readJsonField<bool>(config, "FixDscbTailParams", false);
155162

@@ -159,16 +166,16 @@ void runMassFitter(const std::string& configFileName)
159166
readJsonVector(sliceVarMin, config, "SliceVarMin", true);
160167
readJsonVector(sliceVarMax, config, "SliceVarMax", true);
161168

162-
readJsonVector(massMin, config, "MassMin", true);
163-
readJsonVector(massMax, config, "MassMax", true);
169+
readJsonVectorFlexible(massMin, config, nHistograms, "MassMin", true);
170+
readJsonVectorFlexible(massMax, config, nHistograms, "MassMax", true);
164171

165-
readJsonVector(nRebin, config, "Rebin", true);
172+
readJsonVectorFlexible(nRebin, config, nHistograms, "Rebin", true);
166173

167174
bool const includeSecPeak = readJsonField<bool>(config, "InclSecPeak", false);
168175
bool const useLikelihood = readJsonField<bool>(config, "UseLikelihood");
169176

170-
readJsonVector(bkgFunc, config, "BkgFunc", true);
171-
readJsonVector(sgnFunc, config, "SgnFunc", true);
177+
readJsonVectorFlexible(bkgFunc, config, nHistograms, "BkgFunc", true);
178+
readJsonVectorFlexible(sgnFunc, config, nHistograms, "SgnFunc", true);
172179

173180
const bool enableRefl = readJsonField<bool>(config, "EnableRefl", false);
174181
const bool drawBgPrefit = readJsonField<bool>(config, "DrawBgPrefit", true);
@@ -200,11 +207,10 @@ void runMassFitter(const std::string& configFileName)
200207
readJsonVectorFromHisto(dscbNRLower, config, "DscbParametersFile", "DscbNRLowerHisto");
201208
readJsonVectorFromHisto(dscbNRUpper, config, "DscbParametersFile", "DscbNRUpperHisto");
202209

203-
const int nSliceVarBins = static_cast<int>(sliceVarMin.size());
204-
std::vector<double> sliceVarLimits(nSliceVarBins + 1);
210+
std::vector<double> sliceVarLimits(nHistograms + 1);
205211

206212
auto checkVectorSize = [&](const auto& vec, const std::string& name = "", const bool isEmptyOk = false) {
207-
if (vec.size() != static_cast<size_t>(nSliceVarBins)) {
213+
if (vec.size() != static_cast<size_t>(nHistograms)) {
208214
if (isEmptyOk && vec.empty()) {
209215
return;
210216
}
@@ -253,7 +259,7 @@ void runMassFitter(const std::string& configFileName)
253259
}
254260
};
255261

256-
for (int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) {
262+
for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) {
257263
sliceVarLimits[iSliceVar] = sliceVarMin[iSliceVar];
258264

259265
if (bkgFunc[iSliceVar] < 0 || bkgFunc[iSliceVar] >= HFInvMassFitter::NTypesOfBkgPdf) {
@@ -266,7 +272,7 @@ void runMassFitter(const std::string& configFileName)
266272
throw std::runtime_error("ERROR: only SingleGaus, DoubleGaus and DoubleGausSigmaRatioPar signal functions supported! Exit");
267273
}
268274
}
269-
sliceVarLimits[nSliceVarBins] = sliceVarMax[nSliceVarBins - 1];
275+
sliceVarLimits[nHistograms] = sliceVarMax[nHistograms - 1];
270276

271277
struct DecayInfo {
272278
std::string decayProducts;
@@ -296,11 +302,11 @@ void runMassFitter(const std::string& configFileName)
296302

297303
TFile* inputFileRefl = enableRefl ? openFileWithNullptrCheck(reflFileName) : nullptr;
298304

299-
std::vector<TH1*> hMassSgn(nSliceVarBins);
300-
std::vector<TH1*> hMassRefl(nSliceVarBins);
301-
std::vector<TH1*> hMass(nSliceVarBins);
305+
std::vector<TH1*> hMassSgn(nHistograms);
306+
std::vector<TH1*> hMassRefl(nHistograms);
307+
std::vector<TH1*> hMass(nHistograms);
302308

303-
for (int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) {
309+
for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) {
304310
if (!isMc) {
305311
hMass[iSliceVar] = getObjectWithNullPtrCheck<TH1>(inputFile, inputHistoName[iSliceVar]);
306312
if (enableRefl) {
@@ -340,22 +346,22 @@ void runMassFitter(const std::string& configFileName)
340346
}
341347

342348
// define output histos
343-
auto* hRawYieldsSignal = new TH1D("hRawYieldsSignal", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield", nSliceVarBins, sliceVarLimits.data());
344-
auto* hRawYieldsSignalCounted = new TH1D("hRawYieldsSignalCounted", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield via bin count", nSliceVarBins, sliceVarLimits.data());
345-
auto* hRawYieldsBkg = new TH1D("hRawYieldsBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");Background (3#sigma)", nSliceVarBins, sliceVarLimits.data());
346-
auto* hRawYieldsSgnOverBkg = new TH1D("hRawYieldsSgnOverBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");S/B (3#sigma)", nSliceVarBins, sliceVarLimits.data());
347-
auto* hRawYieldsSignificance = new TH1D("hRawYieldsSignificance", ";" + sliceVarName + "(" + sliceVarUnit + ");significance (3#sigma)", nSliceVarBins, sliceVarLimits.data());
348-
auto* hRawYieldsChiSquareBkg = new TH1D("hRawYieldsChiSquareBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");#chi^{2}/#it{ndf}", nSliceVarBins, sliceVarLimits.data());
349-
auto* hRawYieldsChiSquareTotal = new TH1D("hRawYieldsChiSquareTotal", ";" + sliceVarName + "(" + sliceVarUnit + ");#chi^{2}/#it{ndf}", nSliceVarBins, sliceVarLimits.data());
350-
auto* hReflectionOverSignal = new TH1D("hReflectionOverSignal", ";" + sliceVarName + "(" + sliceVarUnit + ");Refl/Signal", nSliceVarBins, sliceVarLimits.data());
351-
auto* hRawYieldsMean = new TH1D("hRawYieldsMean", ";" + sliceVarName + "(" + sliceVarUnit + ");mean (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data());
352-
auto* hRawYieldsSigma = new TH1D("hRawYieldsSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data());
353-
auto* hRawYieldsSecSigma = new TH1D("hRawYieldsSecSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data());
354-
auto* hRawYieldsFracDoubleGaus = new TH1D("hRawYieldsFracDoubleGaus", ";" + sliceVarName + "(" + sliceVarUnit + ");fraction of double gaussian", nSliceVarBins, sliceVarLimits.data());
355-
auto* hRawYieldsDscbAlphaL = new TH1D("hRawYieldsDscbAlphaL", ";" + sliceVarName + "(" + sliceVarUnit + ");#alpha_{L}", nSliceVarBins, sliceVarLimits.data());
356-
auto* hRawYieldsDscbAlphaR = new TH1D("hRawYieldsDscbAlphaR", ";" + sliceVarName + "(" + sliceVarUnit + ");#alpha_{R}", nSliceVarBins, sliceVarLimits.data());
357-
auto* hRawYieldsDscbNL = new TH1D("hRawYieldsDscbNL", ";" + sliceVarName + "(" + sliceVarUnit + ");n_{L}", nSliceVarBins, sliceVarLimits.data());
358-
auto* hRawYieldsDscbNR = new TH1D("hRawYieldsDscbNR", ";" + sliceVarName + "(" + sliceVarUnit + ");n_{R}", nSliceVarBins, sliceVarLimits.data());
349+
auto* hRawYieldsSignal = new TH1D("hRawYieldsSignal", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield", nHistograms, sliceVarLimits.data());
350+
auto* hRawYieldsSignalCounted = new TH1D("hRawYieldsSignalCounted", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield via bin count", nHistograms, sliceVarLimits.data());
351+
auto* hRawYieldsBkg = new TH1D("hRawYieldsBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");Background (3#sigma)", nHistograms, sliceVarLimits.data());
352+
auto* hRawYieldsSgnOverBkg = new TH1D("hRawYieldsSgnOverBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");S/B (3#sigma)", nHistograms, sliceVarLimits.data());
353+
auto* hRawYieldsSignificance = new TH1D("hRawYieldsSignificance", ";" + sliceVarName + "(" + sliceVarUnit + ");significance (3#sigma)", nHistograms, sliceVarLimits.data());
354+
auto* hRawYieldsChiSquareBkg = new TH1D("hRawYieldsChiSquareBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");#chi^{2}/#it{ndf}", nHistograms, sliceVarLimits.data());
355+
auto* hRawYieldsChiSquareTotal = new TH1D("hRawYieldsChiSquareTotal", ";" + sliceVarName + "(" + sliceVarUnit + ");#chi^{2}/#it{ndf}", nHistograms, sliceVarLimits.data());
356+
auto* hReflectionOverSignal = new TH1D("hReflectionOverSignal", ";" + sliceVarName + "(" + sliceVarUnit + ");Refl/Signal", nHistograms, sliceVarLimits.data());
357+
auto* hRawYieldsMean = new TH1D("hRawYieldsMean", ";" + sliceVarName + "(" + sliceVarUnit + ");mean (GeV/#it{c}^{2})", nHistograms, sliceVarLimits.data());
358+
auto* hRawYieldsSigma = new TH1D("hRawYieldsSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nHistograms, sliceVarLimits.data());
359+
auto* hRawYieldsSecSigma = new TH1D("hRawYieldsSecSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nHistograms, sliceVarLimits.data());
360+
auto* hRawYieldsFracDoubleGaus = new TH1D("hRawYieldsFracDoubleGaus", ";" + sliceVarName + "(" + sliceVarUnit + ");fraction of double gaussian", nHistograms, sliceVarLimits.data());
361+
auto* hRawYieldsDscbAlphaL = new TH1D("hRawYieldsDscbAlphaL", ";" + sliceVarName + "(" + sliceVarUnit + ");#alpha_{L}", nHistograms, sliceVarLimits.data());
362+
auto* hRawYieldsDscbAlphaR = new TH1D("hRawYieldsDscbAlphaR", ";" + sliceVarName + "(" + sliceVarUnit + ");#alpha_{R}", nHistograms, sliceVarLimits.data());
363+
auto* hRawYieldsDscbNL = new TH1D("hRawYieldsDscbNL", ";" + sliceVarName + "(" + sliceVarUnit + ");n_{L}", nHistograms, sliceVarLimits.data());
364+
auto* hRawYieldsDscbNR = new TH1D("hRawYieldsDscbNR", ";" + sliceVarName + "(" + sliceVarUnit + ");n_{R}", nHistograms, sliceVarLimits.data());
359365

360366
enum {
361367
ConfigMassMin = 1,
@@ -366,7 +372,7 @@ void runMassFitter(const std::string& configFileName)
366372
ConfigSgnFunc,
367373
NConfigsToSave
368374
};
369-
auto* hFitConfig = new TH2F("hfitConfig", "Fit Configurations", NConfigsToSave - 1, 0, NConfigsToSave - 1, nSliceVarBins, sliceVarLimits.data());
375+
auto* hFitConfig = new TH2F("hfitConfig", "Fit Configurations", NConfigsToSave - 1, 0, NConfigsToSave - 1, nHistograms, sliceVarLimits.data());
370376
const char* hFitConfigXLabel[NConfigsToSave - 1] = {"mass min", "mass max", "rebin num", "fix sigma", "bkg func", "sgn func"};
371377
hFitConfig->SetStats(false);
372378
for (int i = 0; i < NConfigsToSave - 1; i++) {
@@ -393,15 +399,15 @@ void runMassFitter(const std::string& configFileName)
393399
setHistoStyle(hRawYieldsDscbNL);
394400
setHistoStyle(hRawYieldsDscbNR);
395401

396-
auto getHistToFix = [&nSliceVarBins](bool const& isFix, std::vector<double> const& fixManual, std::string const& fixFileName, std::string const& var) -> TH1* {
402+
auto getHistToFix = [&nHistograms](bool const& isFix, std::vector<double> const& fixManual, std::string const& fixFileName, std::string const& var) -> TH1* {
397403
TH1* histToFix = nullptr;
398404
if (isFix) {
399405
if (fixManual.empty()) {
400406
auto* fixInputFile = openFileWithNullptrCheck(fixFileName);
401407
const std::string histName = "hRawYields" + var;
402408
histToFix = getObjectWithNullPtrCheck<TH1>(fixInputFile, histName);
403409
histToFix->SetDirectory(nullptr);
404-
if (histToFix->GetNbinsX() != nSliceVarBins) {
410+
if (histToFix->GetNbinsX() != nHistograms) {
405411
throw std::runtime_error("Different number of bins for this analysis and histo for fixed " + var);
406412
}
407413
fixInputFile->Close();
@@ -416,19 +422,19 @@ void runMassFitter(const std::string& configFileName)
416422
TH1* hFracDoubleGausToFix = getHistToFix(fixFracDoubleGaus, fixFracDoubleGausManual, fracDoubleGausFile, "FracDoubleGaus");
417423

418424
int canvasSize[2] = {1920, 1080};
419-
if (nSliceVarBins == 1) {
425+
if (nHistograms == 1) {
420426
canvasSize[0] = 500;
421427
canvasSize[1] = 500;
422428
}
423429

424430
int constexpr NCanvasesMax = 20; // do not put more than 20 bins per canvas to make them visible
425-
const int nCanvases = std::ceil(static_cast<float>(nSliceVarBins) / NCanvasesMax);
431+
const int nCanvases = std::ceil(static_cast<float>(nHistograms) / NCanvasesMax);
426432
std::vector<TCanvas*> canvasMass(nCanvases);
427433
std::vector<TCanvas*> canvasResiduals(nCanvases);
428434
std::vector<TCanvas*> canvasRatio(nCanvases);
429435
std::vector<TCanvas*> canvasRefl(nCanvases);
430436
for (int iCanvas = 0; iCanvas < nCanvases; iCanvas++) {
431-
const int nPads = (nCanvases == 1) ? nSliceVarBins : NCanvasesMax;
437+
const int nPads = (nCanvases == 1) ? nHistograms : NCanvasesMax;
432438
canvasMass[iCanvas] = new TCanvas(Form("canvasMass%d", iCanvas), Form("canvasMass%d", iCanvas), canvasSize[0], canvasSize[1]);
433439
divideCanvas(canvasMass[iCanvas], nPads);
434440

@@ -444,7 +450,7 @@ void runMassFitter(const std::string& configFileName)
444450
}
445451
}
446452

447-
for (int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) {
453+
for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) {
448454
const int iCanvas = std::floor(static_cast<float>(iSliceVar) / NCanvasesMax);
449455

450456
hMass[iSliceVar]->Rebin(nRebin[iSliceVar]);
@@ -504,7 +510,7 @@ void runMassFitter(const std::string& configFileName)
504510
}
505511

506512
auto setDscbParameter = [&](const std::vector<double>& vec, void (HFInvMassFitter::*setter)(double)) {
507-
if (static_cast<int>(vec.size()) == nSliceVarBins) {
513+
if (static_cast<int>(vec.size()) == nHistograms) {
508514
(massFitter->*setter)(vec[iSliceVar]);
509515
}
510516
};
@@ -524,7 +530,7 @@ void runMassFitter(const std::string& configFileName)
524530
massFitter->doFit();
525531

526532
auto drawOnCanvas = [&](std::vector<TCanvas*>& canvas, std::function<void()> drawer) {
527-
if (nSliceVarBins > 1) {
533+
if (nHistograms > 1) {
528534
canvas[iCanvas]->cd(iSliceVar - NCanvasesMax * iCanvas + 1);
529535
} else {
530536
canvas[iCanvas]->cd();
@@ -631,7 +637,7 @@ void runMassFitter(const std::string& configFileName)
631637
}
632638
}
633639

634-
for (int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) {
640+
for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) {
635641
hMass[iSliceVar]->Write();
636642
}
637643
hRawYieldsSignal->Write();
@@ -685,10 +691,10 @@ void setHistoStyle(TH1* histo, Color_t color, Size_t markerSize)
685691
histo->SetLineColor(color);
686692
}
687693

688-
void divideCanvas(TCanvas* canvas, int nSliceVarBins)
694+
void divideCanvas(TCanvas* canvas, int nHistograms)
689695
{
690-
int nCols = std::ceil(std::sqrt(nSliceVarBins));
691-
int nRows = std::ceil(static_cast<double>(nSliceVarBins) / nCols);
696+
int nCols = std::ceil(std::sqrt(nHistograms));
697+
int nRows = std::ceil(static_cast<double>(nHistograms) / nCols);
692698
canvas->Divide(nCols, nRows);
693699
}
694700

@@ -766,6 +772,30 @@ void readJsonVector(std::vector<T>& vec, const Document& config, const std::stri
766772
}
767773
}
768774

775+
template <typename T>
776+
void readJsonVectorFlexible(std::vector<T>& vec, const Document& config, int nHistograms, const std::string& fieldName, bool isRequired)
777+
{
778+
if constexpr (!(std::is_same_v<std::decay_t<T>, int> || std::is_same_v<std::decay_t<T>, double>)) {
779+
static_assert(AlwaysFalse<T>, "readJsonVectorFlexible(): unsupported type!");
780+
}
781+
if (!vec.empty()) {
782+
throw std::runtime_error("readJsonVectorFlexible(): vector is not empty!");
783+
}
784+
if (!config.HasMember(fieldName.c_str())) {
785+
if (isRequired) {
786+
throw std::runtime_error("readJsonVectorFlexible(): missing required field " + fieldName);
787+
} else {
788+
return;
789+
}
790+
}
791+
if (config[fieldName.c_str()].IsArray()) {
792+
readJsonVector(vec, config, fieldName);
793+
} else {
794+
const T value = readJsonField<T>(config, fieldName);
795+
vec.assign(nHistograms, value);
796+
}
797+
}
798+
769799
void readJsonVectorFromHisto(std::vector<double>& vec, const Document& config, const std::string& fileNameFieldName, const std::string& histoNameFieldName)
770800
{
771801
if (!vec.empty()) {

0 commit comments

Comments
 (0)