Skip to content

Commit a95b7a4

Browse files
lubynetslubynets
committed
randomize initial parameters
1 parent 38b0367 commit a95b7a4

2 files changed

Lines changed: 43 additions & 25 deletions

File tree

‎PWGHF/D2H/Macros/HFInvMassFitter.cxx‎

Lines changed: 42 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -196,14 +196,16 @@ void HFInvMassFitter::doFit()
196196
dataHistogram.plotOn(mInvMassFrame, Name("data_c")); // plot data histogram on the frame
197197

198198
// define number of background and background fit function
199-
mRooNBkg = new RooRealVar("mRooNBkg", "number of background", 0.3 * mIntegralHisto, 0., 1.2 * mIntegralHisto); // background yield
200-
RooAbsPdf* bkgPdf = createBackgroundFitFunction(mWorkspace); // Create background pdf
201-
RooAbsPdf* sgnPdf = createSignalFitFunction(mWorkspace); // Create signal pdf
199+
const ParameterRanges rooNBkgParamRanges{0., 1.2 * mIntegralHisto, 0.3 * mIntegralHisto};
200+
mRooNBkg = new RooRealVar("mRooNBkg", "number of background", randomizeInitialParameter(rooNBkgParamRanges), rooNBkgParamRanges.lower, rooNBkgParamRanges.upper); // background yield
201+
RooAbsPdf* bkgPdf = createBackgroundFitFunction(mWorkspace); // Create background pdf
202+
RooAbsPdf* sgnPdf = createSignalFitFunction(mWorkspace); // Create signal pdf
202203

203204
// fit MC or Data
204-
if (mTypeOfBkgPdf == NoBkg) { // MC
205-
mRooNSgn = new RooRealVar("mRooNSig", "number of signal", 0.3 * mIntegralHisto, 0., 1.2 * mIntegralHisto); // signal yield
206-
mTotalPdf = new RooAddPdf("mMCFunc", "MC fit function", RooArgList(*sgnPdf), RooArgList(*mRooNSgn)); // create total pdf
205+
if (mTypeOfBkgPdf == NoBkg) { // MC
206+
const ParameterRanges rooNSgnParamRanges{0., 1.2 * mIntegralHisto, 0.3 * mIntegralHisto};
207+
mRooNSgn = new RooRealVar("mRooNSig", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // signal yield
208+
mTotalPdf = new RooAddPdf("mMCFunc", "MC fit function", RooArgList(*sgnPdf), RooArgList(*mRooNSgn)); // create total pdf
207209
if (strcmp(mFitOption.c_str(), "Chi2") == 0) {
208210
mTotalPdf->chi2FitTo(dataHistogram, Range("signal"));
209211
} else {
@@ -250,7 +252,8 @@ void HFInvMassFitter::doFit()
250252
checkForSignal(estimatedSignal); // SIG's absolute integral in "bkg" range
251253
calculateBackground(mBkgYield, mBkgYieldErr); // BG's absolute integral in "bkg" range
252254

253-
mRooNSgn = new RooRealVar("mNSgn", "number of signal", 0.3 * estimatedSignal, 0., 1.2 * estimatedSignal); // estimated signal yield
255+
const ParameterRanges rooNSgnParamRanges{0., 1.2 * estimatedSignal, 0.3 * estimatedSignal};
256+
mRooNSgn = new RooRealVar("mNSgn", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // estimated signal yield
254257
if (mFixedRawYield > 0) {
255258
mRooNSgn->setVal(mFixedRawYield); // fixed signal yield
256259
mRooNSgn->setConstant(true);
@@ -263,7 +266,8 @@ void HFInvMassFitter::doFit()
263266
mReflFrame = mass->frame();
264267
mReflOnlyFrame = mass->frame(Title(Form("%s", mHistoTemplateRefl->GetTitle())));
265268
reflHistogram.plotOn(mReflOnlyFrame);
266-
mRooNRefl = new RooRealVar("mNRefl", "number of reflection", 0.5 * mHistoTemplateRefl->Integral(), 0, mHistoTemplateRefl->Integral());
269+
const ParameterRanges rooNReflParamRanges{0., mHistoTemplateRefl->Integral(), 0.5 * mHistoTemplateRefl->Integral()};
270+
mRooNRefl = new RooRealVar("mNRefl", "number of reflection", randomizeInitialParameter(rooNReflParamRanges), rooNReflParamRanges.lower, rooNReflParamRanges.upper);
267271
RooAddPdf reflFuncTemp("reflFuncTemp", "template reflection fit function", RooArgList(*reflPdf), RooArgList(*mRooNRefl));
268272
if (strcmp(mFitOption.c_str(), "Chi2") == 0) {
269273
reflFuncTemp.chi2FitTo(reflHistogram);
@@ -333,35 +337,42 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const
333337
// Declare observable variable
334338
RooRealVar mass("mass", "mass", mMinMass, mMaxMass, "GeV/c^{2}");
335339
// bkg expo
336-
RooRealVar tau("tau", "tau", -1, -5., 5.);
340+
const ParameterRanges tauParamRanges{-5., 5., -1., 0.1};
341+
RooRealVar tau("tau", "tau", randomizeInitialParameter(tauParamRanges), tauParamRanges.lower, tauParamRanges.upper);
337342
RooAbsPdf* bkgFuncExpo = new RooExponential("bkgFuncExpo", "background fit function", mass, tau);
338343
workspace.import(*bkgFuncExpo);
339344
delete bkgFuncExpo;
340345
// bkg poly1
341-
RooRealVar const polyParam0("polyParam0", "Parameter of Poly function", 0.5, -5., 5.);
342-
RooRealVar const polyParam1("polyParam1", "Parameter of Poly function", 0.2, -5., 5.);
346+
const ParameterRanges polyParam0ParamRanges{-5., 5., 0.5, 0.1};
347+
RooRealVar const polyParam0("polyParam0", "Parameter of Poly function", randomizeInitialParameter(polyParam0ParamRanges), polyParam0ParamRanges.lower, polyParam0ParamRanges.upper);
348+
const ParameterRanges polyParam1ParamRanges{-5., 5., 0.2, 0.05};
349+
RooRealVar const polyParam1("polyParam1", "Parameter of Poly function", randomizeInitialParameter(polyParam1ParamRanges), polyParam1ParamRanges.lower, polyParam1ParamRanges.upper);
343350
RooAbsPdf* bkgFuncPoly1 = new RooPolynomial("bkgFuncPoly1", "background fit function", mass, RooArgSet(polyParam0, polyParam1));
344351
workspace.import(*bkgFuncPoly1);
345352
delete bkgFuncPoly1;
346353
// bkg poly2
347-
RooRealVar const polyParam2("polyParam2", "Parameter of Poly function", 0.2, -5., 5.);
354+
const ParameterRanges polyParam2ParamRanges{-5., 5., 0.2, 0.05};
355+
RooRealVar const polyParam2("polyParam2", "Parameter of Poly function", randomizeInitialParameter(polyParam2ParamRanges), polyParam2ParamRanges.lower, polyParam2ParamRanges.upper);
348356
RooAbsPdf* bkgFuncPoly2 = new RooPolynomial("bkgFuncPoly2", "background fit function", mass, RooArgSet(polyParam0, polyParam1, polyParam2));
349357
workspace.import(*bkgFuncPoly2);
350358
delete bkgFuncPoly2;
351359
// bkg poly3
352-
RooRealVar const polyParam3("polyParam3", "Parameter of Poly function", 0.2, -1., 1.);
360+
const ParameterRanges polyParam3ParamRanges{-1., 1., 0.2, 0.05};
361+
RooRealVar const polyParam3("polyParam3", "Parameter of Poly function", randomizeInitialParameter(polyParam3ParamRanges), polyParam3ParamRanges.lower, polyParam3ParamRanges.upper);
353362
RooAbsPdf* bkgFuncPoly3 = new RooPolynomial("bkgFuncPoly3", "background pdf", mass, RooArgSet(polyParam0, polyParam1, polyParam2, polyParam3));
354363
workspace.import(*bkgFuncPoly3);
355364
delete bkgFuncPoly3;
356365
// bkg power law
357366
RooRealVar const powParam1("powParam1", "Parameter of Pow function", TDatabasePDG::Instance()->GetParticle("pi+")->Mass());
358-
RooRealVar const powParam2("powParam2", "Parameter of Pow function", 1., -10, 10);
367+
const ParameterRanges powParam2ParamRanges{-10., 10., 1., 0.2};
368+
RooRealVar const powParam2("powParam2", "Parameter of Pow function", randomizeInitialParameter(powParam2ParamRanges), powParam2ParamRanges.lower, powParam2ParamRanges.upper);
359369
RooAbsPdf* bkgFuncPow = new RooGenericPdf("bkgFuncPow", "bkgFuncPow", "(mass-powParam1)^powParam2", RooArgSet(mass, powParam1, powParam2));
360370
workspace.import(*bkgFuncPow);
361371
delete bkgFuncPow;
362372
// pow * exp
363373
RooRealVar const powExpoParam1("powExpoParam1", "Parameter of PowExpo function", 1. / 2.);
364-
RooRealVar const powExpoParam2("powExpoParam2", "Parameter of PowExpo function", 1, -10, 10);
374+
const ParameterRanges powExpoParam2ParamRanges{-10., 10., 1., 0.2};
375+
RooRealVar const powExpoParam2("powExpoParam2", "Parameter of PowExpo function", randomizeInitialParameter(powExpoParam2ParamRanges), powExpoParam2ParamRanges.lower, powExpoParam2ParamRanges.upper);
365376
RooRealVar massPi("massPi", "mass of pion", TDatabasePDG::Instance()->GetParticle("pi+")->Mass());
366377
RooFormulaVar powExpoParam3("powExpoParam3", "powExpoParam1 + 1", RooArgList(powExpoParam1));
367378
RooFormulaVar powExpoParam4("powExpoParam4", "1./powExpoParam2", RooArgList(powExpoParam2));
@@ -490,17 +501,24 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const
490501
workspace.import(*reflFuncDoubleGaus);
491502
delete reflFuncDoubleGaus;
492503
// reflection poly3
493-
RooRealVar const polyReflParam0("polyReflParam0", "polyReflParam0", 0.5, -1., 1.);
494-
RooRealVar const polyReflParam1("polyReflParam1", "polyReflParam1", 0.2, -1., 1.);
495-
RooRealVar const polyReflParam2("polyReflParam2", "polyReflParam2", 0.2, -1., 1.);
496-
RooRealVar const polyReflParam3("polyReflParam3", "polyReflParam3", 0.2, -1., 1.);
504+
const ParameterRanges polyReflParam0ParamRanges{-1., 1., 0.5, 0.1};
505+
RooRealVar const polyReflParam0("polyReflParam0", "polyReflParam0", randomizeInitialParameter(polyReflParam0ParamRanges), polyReflParam0ParamRanges.lower, polyReflParam0ParamRanges.upper);
506+
const ParameterRanges polyReflParam1ParamRanges{-1., 1., 0.2, 0.05};
507+
RooRealVar const polyReflParam1("polyReflParam1", "polyReflParam1", randomizeInitialParameter(polyReflParam1ParamRanges), polyReflParam1ParamRanges.lower, polyReflParam1ParamRanges.upper);
508+
const ParameterRanges polyReflParam2ParamRanges{-1., 1., 0.2, 0.05};
509+
RooRealVar const polyReflParam2("polyReflParam2", "polyReflParam2", randomizeInitialParameter(polyReflParam2ParamRanges), polyReflParam2ParamRanges.lower, polyReflParam2ParamRanges.upper);
510+
const ParameterRanges polyReflParam3ParamRanges{-1., 1., 0.2, 0.05};
511+
RooRealVar const polyReflParam3("polyReflParam3", "polyReflParam3", randomizeInitialParameter(polyReflParam3ParamRanges), polyReflParam3ParamRanges.lower, polyReflParam3ParamRanges.upper);
497512
RooAbsPdf* reflFuncPoly3 = new RooPolynomial("reflFuncPoly3", "reflection PDF", mass, RooArgSet(polyReflParam0, polyReflParam1, polyReflParam2, polyReflParam3));
498513
workspace.import(*reflFuncPoly3);
499514
delete reflFuncPoly3;
500515
// reflection poly6
501-
RooRealVar const polyReflParam4("polyReflParam4", "polyReflParam4", 0.2, -1., 1.);
502-
RooRealVar const polyReflParam5("polyReflParam5", "polyReflParam5", 0.2, -1., 1.);
503-
RooRealVar const polyReflParam6("polyReflParam6", "polyReflParam6", 0.2, -1., 1.);
516+
const ParameterRanges polyReflParam4ParamRanges{-1., 1., 0.2, 0.05};
517+
RooRealVar const polyReflParam4("polyReflParam4", "polyReflParam4", randomizeInitialParameter(polyReflParam4ParamRanges), polyReflParam4ParamRanges.lower, polyReflParam4ParamRanges.upper);
518+
const ParameterRanges polyReflParam5ParamRanges{-1., 1., 0.2, 0.05};
519+
RooRealVar const polyReflParam5("polyReflParam5", "polyReflParam5", randomizeInitialParameter(polyReflParam5ParamRanges), polyReflParam5ParamRanges.lower, polyReflParam5ParamRanges.upper);
520+
const ParameterRanges polyReflParam6ParamRanges{-1., 1., 0.2, 0.05};
521+
RooRealVar const polyReflParam6("polyReflParam6", "polyReflParam6", randomizeInitialParameter(polyReflParam6ParamRanges), polyReflParam6ParamRanges.lower, polyReflParam6ParamRanges.upper);
504522
RooAbsPdf* reflFuncPoly6 = new RooPolynomial("reflFuncPoly6", "reflection pdf", mass, RooArgSet(polyReflParam0, polyReflParam1, polyReflParam2, polyReflParam3, polyReflParam4, polyReflParam5, polyReflParam6));
505523
workspace.import(*reflFuncPoly6);
506524
delete reflFuncPoly6;
@@ -999,7 +1017,7 @@ void HFInvMassFitter::setTemplateReflections(TH1* histoRefl)
9991017
mHistoTemplateRefl->SetName("mHistoTemplateRefl");
10001018
}
10011019

1002-
double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& parameterRanges)
1020+
double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& parameterRanges) const
10031021
{
10041022
constexpr double DefaultSigmaFraction{10.};
10051023
constexpr int MaximalNumberOfIterations{20};
@@ -1015,7 +1033,7 @@ double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& paramet
10151033
result = mRandomGen->Gaus(parameterRanges.initial, sigma);
10161034
++nIter;
10171035
if (nIter > MaximalNumberOfIterations) {
1018-
printf("randomizeInitialFitParameter() - long while loop with lower = %f upper = %f initial = %f sigma = %f\n", parameterRanges.lower, parameterRanges.upper, parameterRanges.initial, sigma);
1036+
printf("randomizeInitialParameter() - long while loop with lower = %f upper = %f initial = %f sigma = %f\n", parameterRanges.lower, parameterRanges.upper, parameterRanges.initial, sigma);
10191037
throw;
10201038
}
10211039
} while (result < parameterRanges.lower || result > parameterRanges.upper);

‎PWGHF/D2H/Macros/HFInvMassFitter.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -151,7 +151,7 @@ class HFInvMassFitter : public TNamed
151151
HFInvMassFitter& operator=(const HFInvMassFitter& source);
152152
void fillWorkspace(RooWorkspace& w) const;
153153
void highlightPeakRegion(const RooPlot* plot, Color_t color = kGray + 1, Width_t width = 1, Style_t style = 2) const;
154-
double randomizeInitialParameter(const ParameterRanges& parameterRanges);
154+
double randomizeInitialParameter(const ParameterRanges& parameterRanges) const;
155155

156156
TH1* mHistoInvMass; // histogram to fit
157157
std::string mFitOption;

0 commit comments

Comments
 (0)