Skip to content

Commit 570144e

Browse files
lubynetslubynets
committed
randomize initial parameters
1 parent 3ec2aed commit 570144e

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
@@ -195,14 +195,16 @@ void HFInvMassFitter::doFit()
195195
dataHistogram.plotOn(mInvMassFrame, Name("data_c")); // plot data histogram on the frame
196196

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

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

252-
mRooNSgn = new RooRealVar("mNSgn", "number of signal", 0.3 * estimatedSignal, 0., 1.2 * estimatedSignal); // estimated signal yield
254+
const ParameterRanges rooNSgnParamRanges{0., 1.2 * estimatedSignal, 0.3 * estimatedSignal};
255+
mRooNSgn = new RooRealVar("mNSgn", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // estimated signal yield
253256
if (mFixedRawYield > 0) {
254257
mRooNSgn->setVal(mFixedRawYield); // fixed signal yield
255258
mRooNSgn->setConstant(true);
@@ -262,7 +265,8 @@ void HFInvMassFitter::doFit()
262265
mReflFrame = mass->frame();
263266
mReflOnlyFrame = mass->frame(Title(Form("%s", mHistoTemplateRefl->GetTitle())));
264267
reflHistogram.plotOn(mReflOnlyFrame);
265-
mRooNRefl = new RooRealVar("mNRefl", "number of reflection", 0.5 * mHistoTemplateRefl->Integral(), 0, mHistoTemplateRefl->Integral());
268+
const ParameterRanges rooNReflParamRanges{0., mHistoTemplateRefl->Integral(), 0.5 * mHistoTemplateRefl->Integral()};
269+
mRooNRefl = new RooRealVar("mNRefl", "number of reflection", randomizeInitialParameter(rooNReflParamRanges), rooNReflParamRanges.lower, rooNReflParamRanges.upper);
266270
RooAddPdf reflFuncTemp("reflFuncTemp", "template reflection fit function", RooArgList(*reflPdf), RooArgList(*mRooNRefl));
267271
if (strcmp(mFitOption.c_str(), "Chi2") == 0) {
268272
reflFuncTemp.chi2FitTo(reflHistogram);
@@ -332,35 +336,42 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const
332336
// Declare observable variable
333337
RooRealVar mass("mass", "mass", mMinMass, mMaxMass, "GeV/c^{2}");
334338
// bkg expo
335-
RooRealVar tau("tau", "tau", -1, -5., 5.);
339+
const ParameterRanges tauParamRanges{-5., 5., -1., 0.1};
340+
RooRealVar tau("tau", "tau", randomizeInitialParameter(tauParamRanges), tauParamRanges.lower, tauParamRanges.upper);
336341
RooAbsPdf* bkgFuncExpo = new RooExponential("bkgFuncExpo", "background fit function", mass, tau);
337342
workspace.import(*bkgFuncExpo);
338343
delete bkgFuncExpo;
339344
// bkg poly1
340-
RooRealVar const polyParam0("polyParam0", "Parameter of Poly function", 0.5, -5., 5.);
341-
RooRealVar const polyParam1("polyParam1", "Parameter of Poly function", 0.2, -5., 5.);
345+
const ParameterRanges polyParam0ParamRanges{-5., 5., 0.5, 0.1};
346+
RooRealVar const polyParam0("polyParam0", "Parameter of Poly function", randomizeInitialParameter(polyParam0ParamRanges), polyParam0ParamRanges.lower, polyParam0ParamRanges.upper);
347+
const ParameterRanges polyParam1ParamRanges{-5., 5., 0.2, 0.05};
348+
RooRealVar const polyParam1("polyParam1", "Parameter of Poly function", randomizeInitialParameter(polyParam1ParamRanges), polyParam1ParamRanges.lower, polyParam1ParamRanges.upper);
342349
RooAbsPdf* bkgFuncPoly1 = new RooPolynomial("bkgFuncPoly1", "background fit function", mass, RooArgSet(polyParam0, polyParam1));
343350
workspace.import(*bkgFuncPoly1);
344351
delete bkgFuncPoly1;
345352
// bkg poly2
346-
RooRealVar const polyParam2("polyParam2", "Parameter of Poly function", 0.2, -5., 5.);
353+
const ParameterRanges polyParam2ParamRanges{-5., 5., 0.2, 0.05};
354+
RooRealVar const polyParam2("polyParam2", "Parameter of Poly function", randomizeInitialParameter(polyParam2ParamRanges), polyParam2ParamRanges.lower, polyParam2ParamRanges.upper);
347355
RooAbsPdf* bkgFuncPoly2 = new RooPolynomial("bkgFuncPoly2", "background fit function", mass, RooArgSet(polyParam0, polyParam1, polyParam2));
348356
workspace.import(*bkgFuncPoly2);
349357
delete bkgFuncPoly2;
350358
// bkg poly3
351-
RooRealVar const polyParam3("polyParam3", "Parameter of Poly function", 0.2, -1., 1.);
359+
const ParameterRanges polyParam3ParamRanges{-1., 1., 0.2, 0.05};
360+
RooRealVar const polyParam3("polyParam3", "Parameter of Poly function", randomizeInitialParameter(polyParam3ParamRanges), polyParam3ParamRanges.lower, polyParam3ParamRanges.upper);
352361
RooAbsPdf* bkgFuncPoly3 = new RooPolynomial("bkgFuncPoly3", "background pdf", mass, RooArgSet(polyParam0, polyParam1, polyParam2, polyParam3));
353362
workspace.import(*bkgFuncPoly3);
354363
delete bkgFuncPoly3;
355364
// bkg power law
356365
RooRealVar const powParam1("powParam1", "Parameter of Pow function", TDatabasePDG::Instance()->GetParticle("pi+")->Mass());
357-
RooRealVar const powParam2("powParam2", "Parameter of Pow function", 1., -10, 10);
366+
const ParameterRanges powParam2ParamRanges{-10., 10., 1., 0.2};
367+
RooRealVar const powParam2("powParam2", "Parameter of Pow function", randomizeInitialParameter(powParam2ParamRanges), powParam2ParamRanges.lower, powParam2ParamRanges.upper);
358368
RooAbsPdf* bkgFuncPow = new RooGenericPdf("bkgFuncPow", "bkgFuncPow", "(mass-powParam1)^powParam2", RooArgSet(mass, powParam1, powParam2));
359369
workspace.import(*bkgFuncPow);
360370
delete bkgFuncPow;
361371
// pow * exp
362372
RooRealVar const powExpoParam1("powExpoParam1", "Parameter of PowExpo function", 1. / 2.);
363-
RooRealVar const powExpoParam2("powExpoParam2", "Parameter of PowExpo function", 1, -10, 10);
373+
const ParameterRanges powExpoParam2ParamRanges{-10., 10., 1., 0.2};
374+
RooRealVar const powExpoParam2("powExpoParam2", "Parameter of PowExpo function", randomizeInitialParameter(powExpoParam2ParamRanges), powExpoParam2ParamRanges.lower, powExpoParam2ParamRanges.upper);
364375
RooRealVar massPi("massPi", "mass of pion", TDatabasePDG::Instance()->GetParticle("pi+")->Mass());
365376
RooFormulaVar powExpoParam3("powExpoParam3", "powExpoParam1 + 1", RooArgList(powExpoParam1));
366377
RooFormulaVar powExpoParam4("powExpoParam4", "1./powExpoParam2", RooArgList(powExpoParam2));
@@ -489,17 +500,24 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const
489500
workspace.import(*reflFuncDoubleGaus);
490501
delete reflFuncDoubleGaus;
491502
// reflection poly3
492-
RooRealVar const polyReflParam0("polyReflParam0", "polyReflParam0", 0.5, -1., 1.);
493-
RooRealVar const polyReflParam1("polyReflParam1", "polyReflParam1", 0.2, -1., 1.);
494-
RooRealVar const polyReflParam2("polyReflParam2", "polyReflParam2", 0.2, -1., 1.);
495-
RooRealVar const polyReflParam3("polyReflParam3", "polyReflParam3", 0.2, -1., 1.);
503+
const ParameterRanges polyReflParam0ParamRanges{-1., 1., 0.5, 0.1};
504+
RooRealVar const polyReflParam0("polyReflParam0", "polyReflParam0", randomizeInitialParameter(polyReflParam0ParamRanges), polyReflParam0ParamRanges.lower, polyReflParam0ParamRanges.upper);
505+
const ParameterRanges polyReflParam1ParamRanges{-1., 1., 0.2, 0.05};
506+
RooRealVar const polyReflParam1("polyReflParam1", "polyReflParam1", randomizeInitialParameter(polyReflParam1ParamRanges), polyReflParam1ParamRanges.lower, polyReflParam1ParamRanges.upper);
507+
const ParameterRanges polyReflParam2ParamRanges{-1., 1., 0.2, 0.05};
508+
RooRealVar const polyReflParam2("polyReflParam2", "polyReflParam2", randomizeInitialParameter(polyReflParam2ParamRanges), polyReflParam2ParamRanges.lower, polyReflParam2ParamRanges.upper);
509+
const ParameterRanges polyReflParam3ParamRanges{-1., 1., 0.2, 0.05};
510+
RooRealVar const polyReflParam3("polyReflParam3", "polyReflParam3", randomizeInitialParameter(polyReflParam3ParamRanges), polyReflParam3ParamRanges.lower, polyReflParam3ParamRanges.upper);
496511
RooAbsPdf* reflFuncPoly3 = new RooPolynomial("reflFuncPoly3", "reflection PDF", mass, RooArgSet(polyReflParam0, polyReflParam1, polyReflParam2, polyReflParam3));
497512
workspace.import(*reflFuncPoly3);
498513
delete reflFuncPoly3;
499514
// reflection poly6
500-
RooRealVar const polyReflParam4("polyReflParam4", "polyReflParam4", 0.2, -1., 1.);
501-
RooRealVar const polyReflParam5("polyReflParam5", "polyReflParam5", 0.2, -1., 1.);
502-
RooRealVar const polyReflParam6("polyReflParam6", "polyReflParam6", 0.2, -1., 1.);
515+
const ParameterRanges polyReflParam4ParamRanges{-1., 1., 0.2, 0.05};
516+
RooRealVar const polyReflParam4("polyReflParam4", "polyReflParam4", randomizeInitialParameter(polyReflParam4ParamRanges), polyReflParam4ParamRanges.lower, polyReflParam4ParamRanges.upper);
517+
const ParameterRanges polyReflParam5ParamRanges{-1., 1., 0.2, 0.05};
518+
RooRealVar const polyReflParam5("polyReflParam5", "polyReflParam5", randomizeInitialParameter(polyReflParam5ParamRanges), polyReflParam5ParamRanges.lower, polyReflParam5ParamRanges.upper);
519+
const ParameterRanges polyReflParam6ParamRanges{-1., 1., 0.2, 0.05};
520+
RooRealVar const polyReflParam6("polyReflParam6", "polyReflParam6", randomizeInitialParameter(polyReflParam6ParamRanges), polyReflParam6ParamRanges.lower, polyReflParam6ParamRanges.upper);
503521
RooAbsPdf* reflFuncPoly6 = new RooPolynomial("reflFuncPoly6", "reflection pdf", mass, RooArgSet(polyReflParam0, polyReflParam1, polyReflParam2, polyReflParam3, polyReflParam4, polyReflParam5, polyReflParam6));
504522
workspace.import(*reflFuncPoly6);
505523
delete reflFuncPoly6;
@@ -996,7 +1014,7 @@ void HFInvMassFitter::setTemplateReflections(TH1* histoRefl)
9961014
mHistoTemplateRefl->SetName("mHistoTemplateRefl");
9971015
}
9981016

999-
double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& parameterRanges)
1017+
double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& parameterRanges) const
10001018
{
10011019
constexpr double DefaultSigmaFraction{10.};
10021020
constexpr int MaximalNumberOfIterations{20};
@@ -1012,7 +1030,7 @@ double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& paramet
10121030
result = mRandomGen->Gaus(parameterRanges.initial, sigma);
10131031
++nIter;
10141032
if (nIter > MaximalNumberOfIterations) {
1015-
printf("randomizeInitialFitParameter() - long while loop with lower = %f upper = %f initial = %f sigma = %f\n", parameterRanges.lower, parameterRanges.upper, parameterRanges.initial, sigma);
1033+
printf("randomizeInitialParameter() - long while loop with lower = %f upper = %f initial = %f sigma = %f\n", parameterRanges.lower, parameterRanges.upper, parameterRanges.initial, sigma);
10161034
throw;
10171035
}
10181036
} 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)