Skip to content

Commit 095d53a

Browse files
committed
address O2 linter issues; re-implement divideCanvas()
1 parent fb90b1a commit 095d53a

2 files changed

Lines changed: 61 additions & 97 deletions

File tree

PWGHF/D2H/Macros/HFInvMassFitter.cxx

Lines changed: 49 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@
1919

2020
#include "HFInvMassFitter.h"
2121

22+
#include <TDatabasePDG.h>
2223
#include <TLine.h>
2324

2425
#include <RooRealVar.h>
@@ -45,7 +46,7 @@ HFInvMassFitter::HFInvMassFitter() : TNamed(),
4546
mMinMass(0),
4647
mMaxMass(5),
4748
mTypeOfBkgPdf(Expo),
48-
mMassParticle(1.864),
49+
mMassParticle(TDatabasePDG::Instance()->GetParticle("D0")->Mass()),
4950
mTypeOfSgnPdf(SingleGaus),
5051
mTypeOfReflPdf(1),
5152
mMass(1.865),
@@ -112,7 +113,7 @@ HFInvMassFitter::HFInvMassFitter(const TH1* histoToFit, Double_t minValue, Doubl
112113
mMinMass(minValue),
113114
mMaxMass(maxValue),
114115
mTypeOfBkgPdf(fitTypeBkg),
115-
mMassParticle(1.864),
116+
mMassParticle(TDatabasePDG::Instance()->GetParticle("D0")->Mass()),
116117
mTypeOfSgnPdf(fitTypeSgn),
117118
mTypeOfReflPdf(1),
118119
mMass(1.865),
@@ -205,10 +206,10 @@ void HFInvMassFitter::doFit()
205206
RooRealVar* mass = mWorkspace->var("mass");
206207
RooDataHist dataHistogram("dataHistogram", "data", *mass, Import(*mHistoInvMass));
207208

208-
if (mTypeOfBkgPdf == 6) { // MC
209+
if (mTypeOfBkgPdf == NoBkg) { // MC
209210
mass->setRange("signal", mMass - 3. * mSigmaSgn, mMass + 3. * mSigmaSgn);
210211
} else {
211-
if (mTypeOfSgnPdf == 3) { // Second Peak fit range
212+
if (mTypeOfSgnPdf == GausSec) { // Second Peak fit range
212213
mass->setRange("SBL", mMinMass, mMass - mNSigmaForSidebands * mSigmaSgn);
213214
mass->setRange("SBR", mMass + mNSigmaForSidebands * mSigmaSgn, mSecMass - mNSigmaForSidebands * mSecSigma);
214215
mass->setRange("SEC", mSecMass + mNSigmaForSidebands * mSecSigma, mMaxMass);
@@ -230,7 +231,7 @@ void HFInvMassFitter::doFit()
230231
RooAbsPdf* sgnPdf = createSignalFitFunction(mWorkspace); // Create signal pdf
231232

232233
// fit MC or Data
233-
if (mTypeOfBkgPdf == 6) { // MC
234+
if (mTypeOfBkgPdf == NoBkg) { // MC
234235
mRooNSgn = new RooRealVar("mRooNSig", "number of signal", 0.3 * mIntegralHisto, 0., 1.2 * mIntegralHisto); // signal yield
235236
mTotalPdf = new RooAddPdf("mMCFunc", "MC fit function", RooArgList(*sgnPdf), RooArgList(*mRooNSgn)); // create total pdf
236237
if (!strcmp(mFitOption.Data(), "Chi2")) {
@@ -244,7 +245,7 @@ void HFInvMassFitter::doFit()
244245
mTotalPdf->plotOn(mInvMassFrame, Name("Tot_c")); // plot total function
245246
} else { // data
246247
mBkgPdf = new RooAddPdf("mBkgPdf", "background fit function", RooArgList(*bkgPdf), RooArgList(*mRooNBkg));
247-
if (mTypeOfSgnPdf == 3) { // two peak fit
248+
if (mTypeOfSgnPdf == GausSec) { // two peak fit
248249
if (!strcmp(mFitOption.Data(), "Chi2")) {
249250
mBkgPdf->chi2FitTo(dataHistogram, Range("SBL,SBR,SEC"), Save());
250251
} else {
@@ -355,31 +356,31 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const
355356
workspace.import(*bkgFuncExpo);
356357
delete bkgFuncExpo;
357358
// bkg poly1
358-
RooRealVar PolyParam0("PolyParam0", "Parameter of Poly function", 0.5, -5., 5.);
359-
RooRealVar PolyParam1("PolyParam1", "Parameter of Poly function", 0.2, -5., 5.);
360-
RooAbsPdf* bkgFuncPoly1 = new RooPolynomial("bkgFuncPoly1", "background fit function", mass, RooArgSet(PolyParam0, PolyParam1));
359+
RooRealVar polyParam0("polyParam0", "Parameter of Poly function", 0.5, -5., 5.);
360+
RooRealVar polyParam1("polyParam1", "Parameter of Poly function", 0.2, -5., 5.);
361+
RooAbsPdf* bkgFuncPoly1 = new RooPolynomial("bkgFuncPoly1", "background fit function", mass, RooArgSet(polyParam0, polyParam1));
361362
workspace.import(*bkgFuncPoly1);
362363
delete bkgFuncPoly1;
363364
// bkg poly2
364-
RooRealVar PolyParam2("PolyParam2", "Parameter of Poly function", 0.2, -5., 5.);
365-
RooAbsPdf* bkgFuncPoly2 = new RooPolynomial("bkgFuncPoly2", "background fit function", mass, RooArgSet(PolyParam0, PolyParam1, PolyParam2));
365+
RooRealVar polyParam2("polyParam2", "Parameter of Poly function", 0.2, -5., 5.);
366+
RooAbsPdf* bkgFuncPoly2 = new RooPolynomial("bkgFuncPoly2", "background fit function", mass, RooArgSet(polyParam0, polyParam1, polyParam2));
366367
workspace.import(*bkgFuncPoly2);
367368
delete bkgFuncPoly2;
368369
// bkg poly3
369-
RooRealVar PolyParam3("PolyParam3", "Parameter of Poly function", 0.2, -1., 1.);
370-
RooAbsPdf* bkgFuncPoly3 = new RooPolynomial("bkgFuncPoly3", "background pdf", mass, RooArgSet(PolyParam0, PolyParam1, PolyParam2, PolyParam3));
370+
RooRealVar polyParam3("polyParam3", "Parameter of Poly function", 0.2, -1., 1.);
371+
RooAbsPdf* bkgFuncPoly3 = new RooPolynomial("bkgFuncPoly3", "background pdf", mass, RooArgSet(polyParam0, polyParam1, polyParam2, polyParam3));
371372
workspace.import(*bkgFuncPoly3);
372373
delete bkgFuncPoly3;
373374
// bkg power law
374-
RooRealVar PowParam1("PowParam1", "Parameter of Pow function", 0.13957);
375+
RooRealVar PowParam1("PowParam1", "Parameter of Pow function", TDatabasePDG::Instance()->GetParticle("pi+")->Mass());
375376
RooRealVar PowParam2("PowParam2", "Parameter of Pow function", 1., -10, 10);
376377
RooAbsPdf* bkgFuncPow = new RooGenericPdf("bkgFuncPow", "bkgFuncPow", "(mass-PowParam1)^PowParam2", RooArgSet(mass, PowParam1, PowParam2));
377378
workspace.import(*bkgFuncPow);
378379
delete bkgFuncPow;
379380
// pow * exp
380381
RooRealVar PowExpoParam1("PowExpoParam1", "Parameter of PowExpo function", 1 / 2);
381382
RooRealVar PowExpoParam2("PowExpoParam2", "Parameter of PowExpo function", 1, -10, 10);
382-
RooRealVar massPi("massPi", "mass of pion", 0.13957);
383+
RooRealVar massPi("massPi", "mass of pion", TDatabasePDG::Instance()->GetParticle("pi+")->Mass());
383384
RooFormulaVar PowExpoParam3("PowExpoParam3", "PowExpoParam1 + 1", RooArgList(PowExpoParam1));
384385
RooFormulaVar PowExpoParam4("PowExpoParam4", "1./PowExpoParam2", RooArgList(PowExpoParam2));
385386
RooAbsPdf* bkgFuncPowExpo = new RooGamma("bkgFuncPowExpo", "background pdf", mass, PowExpoParam3, PowExpoParam4, massPi);
@@ -507,18 +508,18 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const
507508
workspace.import(*reflFuncDoubleGaus);
508509
delete reflFuncDoubleGaus;
509510
// reflection poly3
510-
RooRealVar PolyReflParam0("PolyReflParam0", "PolyReflParam0", 0.5, -1., 1.);
511-
RooRealVar PolyReflParam1("PolyReflParam1", "PolyReflParam1", 0.2, -1., 1.);
512-
RooRealVar PolyReflParam2("PolyReflParam2", "PolyReflParam2", 0.2, -1., 1.);
513-
RooRealVar PolyReflParam3("PolyReflParam3", "PolyReflParam3", 0.2, -1., 1.);
514-
RooAbsPdf* reflFuncPoly3 = new RooPolynomial("reflFuncPoly3", "reflection PDF", mass, RooArgSet(PolyReflParam0, PolyReflParam1, PolyReflParam2, PolyReflParam3));
511+
RooRealVar polyReflParam0("polyReflParam0", "polyReflParam0", 0.5, -1., 1.);
512+
RooRealVar polyReflParam1("polyReflParam1", "polyReflParam1", 0.2, -1., 1.);
513+
RooRealVar polyReflParam2("polyReflParam2", "polyReflParam2", 0.2, -1., 1.);
514+
RooRealVar polyReflParam3("polyReflParam3", "polyReflParam3", 0.2, -1., 1.);
515+
RooAbsPdf* reflFuncPoly3 = new RooPolynomial("reflFuncPoly3", "reflection PDF", mass, RooArgSet(polyReflParam0, polyReflParam1, polyReflParam2, polyReflParam3));
515516
workspace.import(*reflFuncPoly3);
516517
delete reflFuncPoly3;
517518
// reflection poly6
518-
RooRealVar PolyReflParam4("PolyReflParam4", "PolyReflParam4", 0.2, -1., 1.);
519-
RooRealVar PolyReflParam5("PolyReflParam5", "PolyReflParam5", 0.2, -1., 1.);
520-
RooRealVar PolyReflParam6("PolyReflParam6", "PolyReflParam6", 0.2, -1., 1.);
521-
RooAbsPdf* reflFuncPoly6 = new RooPolynomial("reflFuncPoly6", "reflection pdf", mass, RooArgSet(PolyReflParam0, PolyReflParam1, PolyReflParam2, PolyReflParam3, PolyReflParam4, PolyReflParam5, PolyReflParam6));
519+
RooRealVar polyReflParam4("polyReflParam4", "polyReflParam4", 0.2, -1., 1.);
520+
RooRealVar polyReflParam5("polyReflParam5", "polyReflParam5", 0.2, -1., 1.);
521+
RooRealVar polyReflParam6("polyReflParam6", "polyReflParam6", 0.2, -1., 1.);
522+
RooAbsPdf* reflFuncPoly6 = new RooPolynomial("reflFuncPoly6", "reflection pdf", mass, RooArgSet(polyReflParam0, polyReflParam1, polyReflParam2, polyReflParam3, polyReflParam4, polyReflParam5, polyReflParam6));
522523
workspace.import(*reflFuncPoly6);
523524
delete reflFuncPoly6;
524525
}
@@ -539,14 +540,14 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, Int_t writeFitInfo)
539540
textInfoRight->SetTextColor(kBlue);
540541
textInfoLeft->AddText(Form("S = %.0f #pm %.0f ", mRawYield, mRawYieldErr));
541542
textInfoLeft->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr));
542-
if (mTypeOfBkgPdf != 6) {
543+
if (mTypeOfBkgPdf != NoBkg) {
543544
textInfoLeft->AddText(Form("B (%d#sigma) = %.0f #pm %.0f", mNSigmaForSidebands, mBkgYield, mBkgYieldErr));
544545
textInfoLeft->AddText(Form("S/B (%d#sigma) = %.4g ", mNSigmaForSidebands, mRawYield / mBkgYield));
545546
}
546547
if (mReflPdf) {
547548
textInfoLeft->AddText(Form("Refl/Sig = %.3f #pm %.3f ", mReflOverSgn, 0.0));
548549
}
549-
if (mTypeOfBkgPdf != 6) {
550+
if (mTypeOfBkgPdf != NoBkg) {
550551
textInfoLeft->AddText(Form("Signif (%d#sigma) = %.1f #pm %.1f ", mNSigmaForSidebands, mSignificance, mSignificanceErr));
551552
textInfoLeft->AddText(Form("#chi^{2} / ndf = %.3f", mChiSquareOverNdf));
552553
}
@@ -856,30 +857,30 @@ void HFInvMassFitter::setReflFuncFixed()
856857
fracRefl->setConstant(kTRUE);
857858
} break;
858859
case 2: {
859-
RooRealVar* PolyReflParam0 = mWorkspace->var("PolyReflParam0");
860-
RooRealVar* PolyReflParam1 = mWorkspace->var("PolyReflParam1");
861-
RooRealVar* PolyReflParam2 = mWorkspace->var("PolyReflParam2");
862-
RooRealVar* PolyReflParam3 = mWorkspace->var("PolyReflParam3");
863-
PolyReflParam0->setConstant(kTRUE);
864-
PolyReflParam1->setConstant(kTRUE);
865-
PolyReflParam2->setConstant(kTRUE);
866-
PolyReflParam3->setConstant(kTRUE);
860+
RooRealVar* polyReflParam0 = mWorkspace->var("polyReflParam0");
861+
RooRealVar* polyReflParam1 = mWorkspace->var("polyReflParam1");
862+
RooRealVar* polyReflParam2 = mWorkspace->var("polyReflParam2");
863+
RooRealVar* polyReflParam3 = mWorkspace->var("polyReflParam3");
864+
polyReflParam0->setConstant(kTRUE);
865+
polyReflParam1->setConstant(kTRUE);
866+
polyReflParam2->setConstant(kTRUE);
867+
polyReflParam3->setConstant(kTRUE);
867868
} break;
868869
case 3: {
869-
RooRealVar* PolyReflParam0 = mWorkspace->var("PolyReflParam0");
870-
RooRealVar* PolyReflParam1 = mWorkspace->var("PolyReflParam1");
871-
RooRealVar* PolyReflParam2 = mWorkspace->var("PolyReflParam2");
872-
RooRealVar* PolyReflParam3 = mWorkspace->var("PolyReflParam3");
873-
RooRealVar* PolyReflParam4 = mWorkspace->var("PolyReflParam4");
874-
RooRealVar* PolyReflParam5 = mWorkspace->var("PolyReflParam5");
875-
RooRealVar* PolyReflParam6 = mWorkspace->var("PolyReflParam6");
876-
PolyReflParam0->setConstant(kTRUE);
877-
PolyReflParam1->setConstant(kTRUE);
878-
PolyReflParam2->setConstant(kTRUE);
879-
PolyReflParam3->setConstant(kTRUE);
880-
PolyReflParam4->setConstant(kTRUE);
881-
PolyReflParam5->setConstant(kTRUE);
882-
PolyReflParam6->setConstant(kTRUE);
870+
RooRealVar* polyReflParam0 = mWorkspace->var("polyReflParam0");
871+
RooRealVar* polyReflParam1 = mWorkspace->var("polyReflParam1");
872+
RooRealVar* polyReflParam2 = mWorkspace->var("polyReflParam2");
873+
RooRealVar* polyReflParam3 = mWorkspace->var("polyReflParam3");
874+
RooRealVar* polyReflParam4 = mWorkspace->var("polyReflParam4");
875+
RooRealVar* polyReflParam5 = mWorkspace->var("polyReflParam5");
876+
RooRealVar* polyReflParam6 = mWorkspace->var("polyReflParam6");
877+
polyReflParam0->setConstant(kTRUE);
878+
polyReflParam1->setConstant(kTRUE);
879+
polyReflParam2->setConstant(kTRUE);
880+
polyReflParam3->setConstant(kTRUE);
881+
polyReflParam4->setConstant(kTRUE);
882+
polyReflParam5->setConstant(kTRUE);
883+
polyReflParam6->setConstant(kTRUE);
883884
} break;
884885
default:
885886
break;

PWGHF/D2H/Macros/runMassFitter.C

Lines changed: 12 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -168,39 +168,21 @@ int runMassFitter(const TString& configFileName)
168168
sliceVarLimits[iSliceVar] = sliceVarMin[iSliceVar];
169169
sliceVarLimits[iSliceVar + 1] = sliceVarMax[iSliceVar];
170170

171-
if (bkgFuncConfig[iSliceVar] == 0) {
172-
bkgFunc[iSliceVar] = HFInvMassFitter::Expo;
173-
} else if (bkgFuncConfig[iSliceVar] == 1) {
174-
bkgFunc[iSliceVar] = HFInvMassFitter::Poly1;
175-
} else if (bkgFuncConfig[iSliceVar] == 2) {
176-
bkgFunc[iSliceVar] = HFInvMassFitter::Poly2;
177-
} else if (bkgFuncConfig[iSliceVar] == 3) {
178-
bkgFunc[iSliceVar] = HFInvMassFitter::Pow;
179-
} else if (bkgFuncConfig[iSliceVar] == 4) {
180-
bkgFunc[iSliceVar] = HFInvMassFitter::PowExpo;
181-
} else if (bkgFuncConfig[iSliceVar] == 5) {
182-
bkgFunc[iSliceVar] = HFInvMassFitter::Poly3;
183-
} else if (bkgFuncConfig[iSliceVar] == 6) {
184-
bkgFunc[iSliceVar] = HFInvMassFitter::NoBkg;
185-
} else {
171+
if (bkgFuncConfig[iSliceVar] < HFInvMassFitter::Expo || bkgFuncConfig[iSliceVar] > HFInvMassFitter::NoBkg) {
186172
std::cerr << "ERROR: only Expo, Poly1, Poly2, Pow and PowEx background "
187173
"functions supported! Exit"
188174
<< std::endl;
189175
return -1;
190176
}
177+
bkgFunc[iSliceVar] = bkgFuncConfig[iSliceVar];
191178

192-
if (sgnFuncConfig[iSliceVar] == 0) {
193-
sgnFunc[iSliceVar] = HFInvMassFitter::SingleGaus;
194-
} else if (sgnFuncConfig[iSliceVar] == 1) {
195-
sgnFunc[iSliceVar] = HFInvMassFitter::DoubleGaus;
196-
} else if (sgnFuncConfig[iSliceVar] == 2) {
197-
sgnFunc[iSliceVar] = HFInvMassFitter::DoubleGausSigmaRatioPar;
198-
} else {
179+
if (sgnFuncConfig[iSliceVar] < HFInvMassFitter::SingleGaus || sgnFuncConfig[iSliceVar] > HFInvMassFitter::DoubleGausSigmaRatioPar) {
199180
std::cerr << "ERROR: only SingleGaus, DoubleGaus and DoubleGausSigmaRatioPar signal "
200181
"functions supported! Exit"
201182
<< std::endl;
202183
return -1;
203184
}
185+
sgnFunc[iSliceVar] = sgnFuncConfig[iSliceVar];
204186
}
205187

206188
TString massAxisTitle = "";
@@ -457,7 +439,8 @@ int runMassFitter(const TString& configFileName)
457439

458440
Double_t reflOverSgn = 0;
459441
double markerSize = 1.;
460-
if (nSliceVarBins > 15) {
442+
const int NSliceVarBinsLarge = 15;
443+
if (nSliceVarBins > NSliceVarBinsLarge) {
461444
markerSize = 0.5;
462445
}
463446

@@ -694,31 +677,11 @@ void setHistoStyle(TH1* histo, Color_t color, Size_t markerSize)
694677

695678
void divideCanvas(TCanvas* canvas, int nSliceVarBins)
696679
{
697-
if (nSliceVarBins < 2) {
698-
canvas->cd();
699-
} else if (nSliceVarBins == 2 || nSliceVarBins == 3) {
700-
canvas->Divide(nSliceVarBins, 1);
701-
} else if (nSliceVarBins == 4 || nSliceVarBins == 6 || nSliceVarBins == 8) {
702-
canvas->Divide(nSliceVarBins / 2, 2);
703-
} else if (nSliceVarBins == 5 || nSliceVarBins == 7) {
704-
canvas->Divide((nSliceVarBins + 1) / 2, 2);
705-
} else if (nSliceVarBins == 9 || nSliceVarBins == 12 || nSliceVarBins == 15) {
706-
canvas->Divide(nSliceVarBins / 3, 3);
707-
} else if (nSliceVarBins == 10 || nSliceVarBins == 11) {
708-
canvas->Divide(4, 3);
709-
} else if (nSliceVarBins == 13 || nSliceVarBins == 14) {
710-
canvas->Divide(5, 3);
711-
} else if (nSliceVarBins > 15 && nSliceVarBins <= 20 && nSliceVarBins % 4 == 0) {
712-
canvas->Divide(nSliceVarBins / 4, 4);
713-
} else if (nSliceVarBins > 15 && nSliceVarBins <= 20 && nSliceVarBins % 4 != 0) {
714-
canvas->Divide(5, 4);
715-
} else if (nSliceVarBins == 21) {
716-
canvas->Divide(7, 3);
717-
} else if (nSliceVarBins > 21 && nSliceVarBins <= 25) {
718-
canvas->Divide(5, 5);
719-
} else if (nSliceVarBins > 25 && nSliceVarBins % 2 == 0) {
720-
canvas->Divide(nSliceVarBins / 2, 2);
721-
} else {
722-
canvas->Divide((nSliceVarBins + 1) / 2, 2);
680+
const int rectangularSideMin = std::floor(std::sqrt(nSliceVarBins));
681+
constexpr int rectangularSidesDiffMax = 2;
682+
for (int rectangularSidesDiff = 0; rectangularSidesDiff < rectangularSidesDiffMax; ++rectangularSidesDiff) {
683+
if (rectangularSideMin * (rectangularSideMin + rectangularSidesDiff) >= nSliceVarBins) {
684+
canvas->Divide(rectangularSideMin + rectangularSidesDiff, rectangularSideMin);
685+
}
723686
}
724687
}

0 commit comments

Comments
 (0)