diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index e56c9b87a28..304f7b50d7a 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -19,9 +19,12 @@ /// \author Oleksii Lubynets /// \author Phil Stahlhut +// o2-linter: disable=name/workflow-file (not an O2Physics workflow) + #include "HFInvMassFitter.h" #include +#include #include #include #include @@ -39,6 +42,7 @@ #include #include #include +#include #include #include #include @@ -52,6 +56,7 @@ #include #include #include +#include #include using namespace RooFit; @@ -60,7 +65,8 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, double minValue, double maxValue, int fitTypeBkg, - int fitTypeSgn) : mHistoInvMass(nullptr), + int fitTypeSgn, + int randomSeed) : mHistoInvMass(nullptr), mFitOption("L,E"), mMinMass(minValue), mMaxMass(maxValue), @@ -86,6 +92,7 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mFixedSigma(false), mFixedSigmaDoubleGaus(false), mBoundSigma(false), + mFixedDscbTailParams(false), mSigmaValue(0.012), mParamSgn(0.1), mFracDoubleGaus(0.2), @@ -106,6 +113,18 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mChiSquareOverNdfTotal(0), mChiSquareOverNdfBkg(0), mFixReflOverSgn(false), + mDscbAlphaLInitialValue(1.5), + mDscbAlphaLLowLimit(0.1), + mDscbAlphaLUpLimit(5.0), + mDscbAlphaRInitialValue(1.5), + mDscbAlphaRLowLimit(0.1), + mDscbAlphaRUpLimit(5.0), + mDscbNLInitialValue(2.0), + mDscbNLLowLimit(0.5), + mDscbNLUpLimit(50.), + mDscbNRInitialValue(2.0), + mDscbNRLowLimit(0.5), + mDscbNRUpLimit(50.), mRooMeanSgn(nullptr), mRooSigmaSgn(nullptr), mRooSecSigmaSgn(nullptr), @@ -116,25 +135,35 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mRooNSgn(nullptr), mRooNBkg(nullptr), mRooNRefl(nullptr), + mRooDscbAlphaL(nullptr), + mRooDscbAlphaR(nullptr), + mRooDscbNL(nullptr), + mRooDscbNR(nullptr), mTotalPdf(nullptr), mInvMassFrame(nullptr), mReflFrame(nullptr), mReflOnlyFrame(nullptr), mResidualFrame(nullptr), + mResidualHist(nullptr), mRatioFrame(nullptr), - mResidualFrameForCalculation(nullptr), mWorkspace(nullptr), mIntegralHisto(0), mIntegralBkg(0), mIntegralSgn(0), mHistoTemplateRefl(nullptr), mDrawBgPrefit(false), - mHighlightPeakRegion(false) + mHighlightPeakRegion(false), + mRandomSeed(randomSeed), + mRandomGen(nullptr) { // standard constructor mHistoInvMass = histoToFit; mHistoInvMass->SetName("mHistoInvMass"); mHistoInvMass->SetDirectory(nullptr); + if (mRandomSeed >= 0) { + mRandomGen = new TRandom3(); + mRandomGen->SetSeed(mRandomSeed); + } } HFInvMassFitter::~HFInvMassFitter() @@ -176,40 +205,42 @@ void HFInvMassFitter::doFit() mass->setRange("SBL", mMinMass, mMass - mNSigmaForSidebands * mSigmaSgn); mass->setRange("SBR", mMass + mNSigmaForSidebands * mSigmaSgn, mSecMass - mNSigmaForSidebands * mSecSigma); mass->setRange("SEC", mSecMass + mNSigmaForSidebands * mSecSigma, mMaxMass); - mass->setRange("signal", mSecMass - mNSigmaForSidebands * mSecSigma, mSecMass + mNSigmaForSidebands * mSecSigma); + mass->setRange("signal", mSecMass - mNSigmaForSgn * mSecSigma, mSecMass + mNSigmaForSgn * mSecSigma); } else { // Single Peak fit range mass->setRange("SBL", mMinMass, mMass - mNSigmaForSidebands * mSigmaSgn); mass->setRange("SBR", mMass + mNSigmaForSidebands * mSigmaSgn, mMaxMass); mass->setRange("signal", mMass - mNSigmaForSgn * mSigmaSgn, mMass + mNSigmaForSgn * mSigmaSgn); } } - mass->setRange("bkg", mMass - 4 * mSigmaSgn, mMass + 4 * mSigmaSgn); + mass->setRange("bkg", mMass - mNSigmaForSgn * mSigmaSgn, mMass + mNSigmaForSgn * mSigmaSgn); mass->setRange("full", mMinMass, mMaxMass); mInvMassFrame = mass->frame(Title(Form("%s", mHistoInvMass->GetTitle()))); // define the frame to plot dataHistogram.plotOn(mInvMassFrame, Name("data_c")); // plot data histogram on the frame - // define number of background and background fit function - mRooNBkg = new RooRealVar("mRooNBkg", "number of background", 0.3 * mIntegralHisto, 0., 1.2 * mIntegralHisto); // background yield - RooAbsPdf* bkgPdf = createBackgroundFitFunction(mWorkspace); // Create background pdf - RooAbsPdf* sgnPdf = createSignalFitFunction(mWorkspace); // Create signal pdf + RooAbsPdf* bkgPdf = createBackgroundFitFunction(mWorkspace); // Create background pdf + RooAbsPdf* sgnPdf = createSignalFitFunction(mWorkspace); // Create signal pdf // fit MC or Data - if (mTypeOfBkgPdf == NoBkg) { // MC - mRooNSgn = new RooRealVar("mRooNSig", "number of signal", 0.3 * mIntegralHisto, 0., 1.2 * mIntegralHisto); // signal yield - mTotalPdf = new RooAddPdf("mMCFunc", "MC fit function", RooArgList(*sgnPdf), RooArgList(*mRooNSgn)); // create total pdf + if (mTypeOfBkgPdf == NoBkg) { // MC + const ParameterRanges rooNSgnParamRanges{0., 1.2 * mIntegralHisto, 0.3 * mIntegralHisto}; + mRooNSgn = new RooRealVar("mRooNSig", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // signal yield + mTotalPdf = new RooAddPdf("mMCFunc", "MC fit function", RooArgList(*sgnPdf), RooArgList(*mRooNSgn)); // create total pdf if (strcmp(mFitOption.c_str(), "Chi2") == 0) { - mTotalPdf->chi2FitTo(dataHistogram, Range("signal")); + mTotalPdf->chi2FitTo(dataHistogram, Range("full")); } else { - mTotalPdf->fitTo(dataHistogram, Range("signal")); + mTotalPdf->fitTo(dataHistogram, Range("full")); } RooAbsReal* signalIntegralMc = mTotalPdf->createIntegral(*mass, NormSet(*mass), Range("signal")); // sig yield from fit mIntegralSgn = signalIntegralMc->getValV(); - calculateSignal(mRawYield, mRawYieldErr); // calculate signal and signal error + calculateSignal(mRawYield, mRawYieldErr); // calculate signal and signal error + countSignal(mRawYieldCounted, mRawYieldCountedErr); mTotalPdf->plotOn(mInvMassFrame, Name("Tot_c")); // plot total function // Fit to data ratio mRatioFrame = mass->frame(Title(Form("%s", mHistoInvMass->GetTitle()))); calculateFitToDataRatio(); } else { // data + const ParameterRanges rooNBkgParamRanges{0., 1.2 * mIntegralHisto, 0.3 * mIntegralHisto}; + mRooNBkg = new RooRealVar("mRooNBkg", "number of background", randomizeInitialParameter(rooNBkgParamRanges), rooNBkgParamRanges.lower, rooNBkgParamRanges.upper); // background yield mBkgPdf = new RooAddPdf("mBkgPdf", "background fit function", RooArgList(*bkgPdf), RooArgList(*mRooNBkg)); if (mTypeOfSgnPdf == GausSec) { // two peak fit if (strcmp(mFitOption.c_str(), "Chi2") == 0) { @@ -243,7 +274,8 @@ void HFInvMassFitter::doFit() checkForSignal(estimatedSignal); // SIG's absolute integral in "bkg" range calculateBackground(mBkgYield, mBkgYieldErr); // BG's absolute integral in "bkg" range - mRooNSgn = new RooRealVar("mNSgn", "number of signal", 0.3 * estimatedSignal, 0., 1.2 * estimatedSignal); // estimated signal yield + const ParameterRanges rooNSgnParamRanges{0., 1.2 * estimatedSignal, 0.3 * estimatedSignal}; + mRooNSgn = new RooRealVar("mNSgn", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // estimated signal yield if (mFixedRawYield > 0) { mRooNSgn->setVal(mFixedRawYield); // fixed signal yield mRooNSgn->setConstant(true); @@ -256,7 +288,8 @@ void HFInvMassFitter::doFit() mReflFrame = mass->frame(); mReflOnlyFrame = mass->frame(Title(Form("%s", mHistoTemplateRefl->GetTitle()))); reflHistogram.plotOn(mReflOnlyFrame); - mRooNRefl = new RooRealVar("mNRefl", "number of reflection", 0.5 * mHistoTemplateRefl->Integral(), 0, mHistoTemplateRefl->Integral()); + const ParameterRanges rooNReflParamRanges{0., mHistoTemplateRefl->Integral(), 0.5 * mHistoTemplateRefl->Integral()}; + mRooNRefl = new RooRealVar("mNRefl", "number of reflection", randomizeInitialParameter(rooNReflParamRanges), rooNReflParamRanges.lower, rooNReflParamRanges.upper); RooAddPdf reflFuncTemp("reflFuncTemp", "template reflection fit function", RooArgList(*reflPdf), RooArgList(*mRooNRefl)); if (strcmp(mFitOption.c_str(), "Chi2") == 0) { reflFuncTemp.chi2FitTo(reflHistogram); @@ -283,9 +316,9 @@ void HFInvMassFitter::doFit() mChiSquareOverNdfTotal = mInvMassFrame->chiSquare("Tot_c", "data_c"); // calculate reduced chi2 / NDF // plot residual distribution - RooHist* residualHistogram = mInvMassFrame->residHist("data_c", "ReflBkg_c"); + mResidualHist = mInvMassFrame->residHist("data_c", "ReflBkg_c"); mResidualFrame = mass->frame(Title("Residual Distribution")); - mResidualFrame->addPlotable(residualHistogram, "p"); + mResidualFrame->addPlotable(mResidualHist, "p"); mSgnPdf->plotOn(mResidualFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineColor(kBlue)); } else { mTotalPdf = new RooAddPdf("mTotalPdf", "background + signal pdf", RooArgList(*bkgPdf, *sgnPdf), RooArgList(*mRooNBkg, *mRooNSgn)); @@ -301,8 +334,8 @@ void HFInvMassFitter::doFit() // plot residual distribution mResidualFrame = mass->frame(Title("Residual Distribution")); - RooHist* residualHistogram = mInvMassFrame->residHist("data_c", "Bkg_c"); - mResidualFrame->addPlotable(residualHistogram, "P"); + mResidualHist = mInvMassFrame->residHist("data_c", "Bkg_c"); + mResidualFrame->addPlotable(mResidualHist, "P"); mSgnPdf->plotOn(mResidualFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineColor(kBlue)); } mass->setRange("bkgForSignificance", mRooMeanSgn->getVal() - mNSigmaForSgn * mRooSecSigmaSgn->getVal(), mRooMeanSgn->getVal() + mNSigmaForSgn * mRooSecSigmaSgn->getVal()); @@ -326,35 +359,42 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const // Declare observable variable RooRealVar mass("mass", "mass", mMinMass, mMaxMass, "GeV/c^{2}"); // bkg expo - RooRealVar tau("tau", "tau", -1, -5., 5.); + const ParameterRanges tauParamRanges{-5., 5., -1., 0.1}; + RooRealVar tau("tau", "tau", randomizeInitialParameter(tauParamRanges), tauParamRanges.lower, tauParamRanges.upper); RooAbsPdf* bkgFuncExpo = new RooExponential("bkgFuncExpo", "background fit function", mass, tau); workspace.import(*bkgFuncExpo); delete bkgFuncExpo; // bkg poly1 - RooRealVar const polyParam0("polyParam0", "Parameter of Poly function", 0.5, -5., 5.); - RooRealVar const polyParam1("polyParam1", "Parameter of Poly function", 0.2, -5., 5.); + const ParameterRanges polyParam0ParamRanges{-5., 5., 0.5, 0.1}; + RooRealVar const polyParam0("polyParam0", "Parameter of Poly function", randomizeInitialParameter(polyParam0ParamRanges), polyParam0ParamRanges.lower, polyParam0ParamRanges.upper); + const ParameterRanges polyParam1ParamRanges{-5., 5., 0.2, 0.05}; + RooRealVar const polyParam1("polyParam1", "Parameter of Poly function", randomizeInitialParameter(polyParam1ParamRanges), polyParam1ParamRanges.lower, polyParam1ParamRanges.upper); RooAbsPdf* bkgFuncPoly1 = new RooPolynomial("bkgFuncPoly1", "background fit function", mass, RooArgSet(polyParam0, polyParam1)); workspace.import(*bkgFuncPoly1); delete bkgFuncPoly1; // bkg poly2 - RooRealVar const polyParam2("polyParam2", "Parameter of Poly function", 0.2, -5., 5.); + const ParameterRanges polyParam2ParamRanges{-5., 5., 0.2, 0.05}; + RooRealVar const polyParam2("polyParam2", "Parameter of Poly function", randomizeInitialParameter(polyParam2ParamRanges), polyParam2ParamRanges.lower, polyParam2ParamRanges.upper); RooAbsPdf* bkgFuncPoly2 = new RooPolynomial("bkgFuncPoly2", "background fit function", mass, RooArgSet(polyParam0, polyParam1, polyParam2)); workspace.import(*bkgFuncPoly2); delete bkgFuncPoly2; // bkg poly3 - RooRealVar const polyParam3("polyParam3", "Parameter of Poly function", 0.2, -1., 1.); + const ParameterRanges polyParam3ParamRanges{-1., 1., 0.2, 0.05}; + RooRealVar const polyParam3("polyParam3", "Parameter of Poly function", randomizeInitialParameter(polyParam3ParamRanges), polyParam3ParamRanges.lower, polyParam3ParamRanges.upper); RooAbsPdf* bkgFuncPoly3 = new RooPolynomial("bkgFuncPoly3", "background pdf", mass, RooArgSet(polyParam0, polyParam1, polyParam2, polyParam3)); workspace.import(*bkgFuncPoly3); delete bkgFuncPoly3; // bkg power law RooRealVar const powParam1("powParam1", "Parameter of Pow function", TDatabasePDG::Instance()->GetParticle("pi+")->Mass()); - RooRealVar const powParam2("powParam2", "Parameter of Pow function", 1., -10, 10); + const ParameterRanges powParam2ParamRanges{-10., 10., 1., 0.2}; + RooRealVar const powParam2("powParam2", "Parameter of Pow function", randomizeInitialParameter(powParam2ParamRanges), powParam2ParamRanges.lower, powParam2ParamRanges.upper); RooAbsPdf* bkgFuncPow = new RooGenericPdf("bkgFuncPow", "bkgFuncPow", "(mass-powParam1)^powParam2", RooArgSet(mass, powParam1, powParam2)); workspace.import(*bkgFuncPow); delete bkgFuncPow; // pow * exp RooRealVar const powExpoParam1("powExpoParam1", "Parameter of PowExpo function", 1. / 2.); - RooRealVar const powExpoParam2("powExpoParam2", "Parameter of PowExpo function", 1, -10, 10); + const ParameterRanges powExpoParam2ParamRanges{-10., 10., 1., 0.2}; + RooRealVar const powExpoParam2("powExpoParam2", "Parameter of PowExpo function", randomizeInitialParameter(powExpoParam2ParamRanges), powExpoParam2ParamRanges.lower, powExpoParam2ParamRanges.upper); RooRealVar massPi("massPi", "mass of pion", TDatabasePDG::Instance()->GetParticle("pi+")->Mass()); RooFormulaVar powExpoParam3("powExpoParam3", "powExpoParam1 + 1", RooArgList(powExpoParam1)); RooFormulaVar powExpoParam4("powExpoParam4", "1./powExpoParam2", RooArgList(powExpoParam2)); @@ -482,18 +522,48 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const RooAbsPdf* reflFuncDoubleGaus = new RooAddPdf("reflFuncDoubleGaus", "reflection pdf", RooArgList(gausRefl1, gausRefl2), fracRefl); workspace.import(*reflFuncDoubleGaus); delete reflFuncDoubleGaus; + // signal DSCB pdf + const ParameterRanges dscbAlphaLParamRanges{mDscbAlphaLLowLimit, mDscbAlphaLUpLimit, mDscbAlphaLInitialValue}; + const ParameterRanges dscbNLParamRanges{mDscbNLLowLimit, mDscbNLUpLimit, mDscbNLInitialValue}; + const ParameterRanges dscbAlphaRParamRanges{mDscbAlphaRLowLimit, mDscbAlphaRUpLimit, mDscbAlphaRInitialValue}; + const ParameterRanges dscbNRParamRanges{mDscbNRLowLimit, mDscbNRUpLimit, mDscbNRInitialValue}; + RooRealVar alphaL("alphaL", "left tail alpha", randomizeInitialParameter(dscbAlphaLParamRanges), dscbAlphaLParamRanges.lower, dscbAlphaLParamRanges.upper); + RooRealVar nL("nL", "left tail n", randomizeInitialParameter(dscbNLParamRanges), dscbNLParamRanges.lower, dscbNLParamRanges.upper); + RooRealVar alphaR("alphaR", "right tail alpha", randomizeInitialParameter(dscbAlphaRParamRanges), dscbAlphaRParamRanges.lower, dscbAlphaRParamRanges.upper); + RooRealVar nR("nR", "right tail n", randomizeInitialParameter(dscbNRParamRanges), dscbNRParamRanges.lower, dscbNRParamRanges.upper); + if (mFixedDscbTailParams) { + printf("Fixing DSCB tail parameters to initial values.\n"); + alphaL.setVal(mDscbAlphaLInitialValue); + alphaL.setConstant(true); + alphaR.setVal(mDscbAlphaRInitialValue); + alphaR.setConstant(true); + nL.setVal(mDscbNLInitialValue); + nL.setConstant(true); + nR.setVal(mDscbNRInitialValue); + nR.setConstant(true); + } + RooAbsPdf* sgnFuncDSCB = new RooCrystalBall("sgnFuncDSCB", "double sided crystal ball signal", mass, mean, sigma, alphaL, nL, alphaR, nR); + workspace.import(*sgnFuncDSCB); + delete sgnFuncDSCB; // reflection poly3 - RooRealVar const polyReflParam0("polyReflParam0", "polyReflParam0", 0.5, -1., 1.); - RooRealVar const polyReflParam1("polyReflParam1", "polyReflParam1", 0.2, -1., 1.); - RooRealVar const polyReflParam2("polyReflParam2", "polyReflParam2", 0.2, -1., 1.); - RooRealVar const polyReflParam3("polyReflParam3", "polyReflParam3", 0.2, -1., 1.); + const ParameterRanges polyReflParam0ParamRanges{-1., 1., 0.5, 0.1}; + RooRealVar const polyReflParam0("polyReflParam0", "polyReflParam0", randomizeInitialParameter(polyReflParam0ParamRanges), polyReflParam0ParamRanges.lower, polyReflParam0ParamRanges.upper); + const ParameterRanges polyReflParam1ParamRanges{-1., 1., 0.2, 0.05}; + RooRealVar const polyReflParam1("polyReflParam1", "polyReflParam1", randomizeInitialParameter(polyReflParam1ParamRanges), polyReflParam1ParamRanges.lower, polyReflParam1ParamRanges.upper); + const ParameterRanges polyReflParam2ParamRanges{-1., 1., 0.2, 0.05}; + RooRealVar const polyReflParam2("polyReflParam2", "polyReflParam2", randomizeInitialParameter(polyReflParam2ParamRanges), polyReflParam2ParamRanges.lower, polyReflParam2ParamRanges.upper); + const ParameterRanges polyReflParam3ParamRanges{-1., 1., 0.2, 0.05}; + RooRealVar const polyReflParam3("polyReflParam3", "polyReflParam3", randomizeInitialParameter(polyReflParam3ParamRanges), polyReflParam3ParamRanges.lower, polyReflParam3ParamRanges.upper); RooAbsPdf* reflFuncPoly3 = new RooPolynomial("reflFuncPoly3", "reflection PDF", mass, RooArgSet(polyReflParam0, polyReflParam1, polyReflParam2, polyReflParam3)); workspace.import(*reflFuncPoly3); delete reflFuncPoly3; // reflection poly6 - RooRealVar const polyReflParam4("polyReflParam4", "polyReflParam4", 0.2, -1., 1.); - RooRealVar const polyReflParam5("polyReflParam5", "polyReflParam5", 0.2, -1., 1.); - RooRealVar const polyReflParam6("polyReflParam6", "polyReflParam6", 0.2, -1., 1.); + const ParameterRanges polyReflParam4ParamRanges{-1., 1., 0.2, 0.05}; + RooRealVar const polyReflParam4("polyReflParam4", "polyReflParam4", randomizeInitialParameter(polyReflParam4ParamRanges), polyReflParam4ParamRanges.lower, polyReflParam4ParamRanges.upper); + const ParameterRanges polyReflParam5ParamRanges{-1., 1., 0.2, 0.05}; + RooRealVar const polyReflParam5("polyReflParam5", "polyReflParam5", randomizeInitialParameter(polyReflParam5ParamRanges), polyReflParam5ParamRanges.lower, polyReflParam5ParamRanges.upper); + const ParameterRanges polyReflParam6ParamRanges{-1., 1., 0.2, 0.05}; + RooRealVar const polyReflParam6("polyReflParam6", "polyReflParam6", randomizeInitialParameter(polyReflParam6ParamRanges), polyReflParam6ParamRanges.lower, polyReflParam6ParamRanges.upper); RooAbsPdf* reflFuncPoly6 = new RooPolynomial("reflFuncPoly6", "reflection pdf", mass, RooArgSet(polyReflParam0, polyReflParam1, polyReflParam2, polyReflParam3, polyReflParam4, polyReflParam5, polyReflParam6)); workspace.import(*reflFuncPoly6); delete reflFuncPoly6; @@ -515,9 +585,9 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& textFitMetrics->AddText(Form("S = %.0f #pm %.0f ", mRawYield, mRawYieldErr)); textFitMetrics->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr)); if (mTypeOfBkgPdf != NoBkg) { - textFitMetrics->AddText(Form("B (%.1f#sigma) = %.0f #pm %.0f", mNSigmaForSidebands, mBkgYield, mBkgYieldErr)); - textFitMetrics->AddText(Form("S/B (%.1f#sigma) = %.4g ", mNSigmaForSidebands, mRawYield / mBkgYield)); - textFitMetrics->AddText(Form("Significance (%.1f#sigma) = %.1f #pm %.1f ", mNSigmaForSidebands, mSignificance, mSignificanceErr)); + textFitMetrics->AddText(Form("B (%.1f#sigma) = %.0f #pm %.0f", mNSigmaForSgn, mBkgYield, mBkgYieldErr)); + textFitMetrics->AddText(Form("S/B (%.1f#sigma) = %.4g ", mNSigmaForSgn, mRawYield / mBkgYield)); + textFitMetrics->AddText(Form("Significance (%.1f#sigma) = %.1f #pm %.1f ", mNSigmaForSgn, mSignificance, mSignificanceErr)); textFitMetrics->AddText(Form("#chi^{2} / ndf = %.3f", mChiSquareOverNdfTotal)); } if (mReflPdf != nullptr) { @@ -544,10 +614,10 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& const std::string fixMeanStatus = mFixedMean ? "fixed" : "free"; const std::string fixSigmaStatus = mFixedSigma ? "fixed" : "free"; const std::string fixSigmaDoubleGausStatus = mFixedSigmaDoubleGaus ? "fixed" : "free"; - textSignalPar->AddText(Form("mean(%s) = %.3f #pm %.3f", fixMeanStatus.c_str(), mRooMeanSgn->getVal(), mRooMeanSgn->getError())); - textSignalPar->AddText(Form("sigma(%s) = %.3f #pm %.3f", fixSigmaStatus.c_str(), mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textSignalPar->AddText(Form("#mu(%s) = %.3f #pm %.3f", fixMeanStatus.c_str(), mRooMeanSgn->getVal(), mRooMeanSgn->getError())); + textSignalPar->AddText(Form("#sigma(%s) = %.3f #pm %.3f", fixSigmaStatus.c_str(), mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); if (mTypeOfSgnPdf == DoubleGaus) { - textSignalPar->AddText(Form("sigma 2(%s) = %.3f #pm %.3f", fixSigmaDoubleGausStatus.c_str(), mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); + textSignalPar->AddText(Form("#sigma_{2}(%s) = %.3f #pm %.3f", fixSigmaDoubleGausStatus.c_str(), mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); } mInvMassFrame->addObject(textSignalPar); } @@ -574,10 +644,10 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad) textInfo->SetTextColor(kBlue); textInfo->AddText(Form("S = %.0f #pm %.0f ", mRawYield, mRawYieldErr)); textInfo->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr)); - textInfo->AddText(Form("mean = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError())); - textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textInfo->AddText(Form("#mu = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError())); + textInfo->AddText(Form("#sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); if (mTypeOfSgnPdf == DoubleGaus) { - textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); + textInfo->AddText(Form("#sigma_{2} = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); } mResidualFrame->addObject(textInfo); mResidualFrame->Draw(); @@ -635,29 +705,32 @@ void HFInvMassFitter::drawReflection(TVirtualPad* pad) // calculate signal yield via bin counting void HFInvMassFitter::countSignal(double& signal, double& signalErr) const { - const double mean = mRooMeanSgn->getVal(); - const double sigma = mRooSecSigmaSgn->getVal(); - const double minForSgn = mean - mNSigmaForSidebands * sigma; - const double maxForSgn = mean + mNSigmaForSidebands * sigma; - const int binForMinSgn = mHistoInvMass->FindBin(minForSgn); - const int binForMaxSgn = mHistoInvMass->FindBin(maxForSgn); - const double binForMinSgnUpperEdge = mHistoInvMass->GetBinLowEdge(binForMinSgn + 1); - const double binForMaxSgnLowerEdge = mHistoInvMass->GetBinLowEdge(binForMaxSgn); - const double binForMinSgnFraction = (binForMinSgnUpperEdge - minForSgn) / mHistoInvMass->GetBinWidth(binForMinSgn); - const double binForMaxSgnFraction = (maxForSgn - binForMaxSgnLowerEdge) / mHistoInvMass->GetBinWidth(binForMaxSgn); - - double sum = 0; - sum += mHistoInvMass->GetBinContent(binForMinSgn) * binForMinSgnFraction; + const auto& histoForCounting = mTypeOfBkgPdf == NoBkg ? mInvMassFrame->getHist("data_c") : mResidualHist; + + const double binWidth = histoForCounting->GetX()[1] - histoForCounting->GetX()[0]; + const double firstBinLowEdge = histoForCounting->GetX()[0] - binWidth / 2; + const auto [minForSgn, maxForSgn] = getRangesOfSignal(); + const int binForMinSgn = static_cast((minForSgn - firstBinLowEdge) / binWidth) + 1; + const int binForMaxSgn = static_cast((maxForSgn - firstBinLowEdge) / binWidth) + 1; + const double binForMinSgnUpperEdge = firstBinLowEdge + (binForMinSgn - 1) * binWidth + binWidth; + const double binForMaxSgnLowerEdge = firstBinLowEdge + (binForMaxSgn - 1) * binWidth; + const double binForMinSgnFraction = (binForMinSgnUpperEdge - minForSgn) / binWidth; + const double binForMaxSgnFraction = (maxForSgn - binForMaxSgnLowerEdge) / binWidth; + + auto square = [](double value) { return value * value; }; + + double sumValues{}, sumErrorsSquare{}; + sumValues += histoForCounting->GetY()[binForMinSgn - 1] * binForMinSgnFraction; + sumErrorsSquare += square(histoForCounting->GetErrorY(binForMinSgn - 1) * binForMinSgnFraction); for (int iBin = binForMinSgn + 1; iBin <= binForMaxSgn - 1; iBin++) { - sum += mHistoInvMass->GetBinContent(iBin); + sumValues += histoForCounting->GetY()[iBin - 1]; + sumErrorsSquare += square(histoForCounting->GetErrorY(iBin - 1)); } - sum += mHistoInvMass->GetBinContent(binForMaxSgn) * binForMaxSgnFraction; + sumValues += histoForCounting->GetY()[binForMaxSgn - 1] * binForMaxSgnFraction; + sumErrorsSquare += square(histoForCounting->GetErrorY(binForMaxSgn - 1) * binForMaxSgnFraction); - double bkg{}, errBkg{}; - calculateBackground(bkg, errBkg); - - signal = sum - bkg; - signalErr = std::sqrt(sum + errBkg * errBkg); // sum error squared is equal to sum + signal = sumValues; + signalErr = std::sqrt(sumErrorsSquare); } // calculate signal yield @@ -670,6 +743,11 @@ void HFInvMassFitter::calculateSignal(double& signal, double& errSignal) const // calculate background yield void HFInvMassFitter::calculateBackground(double& bkg, double& errBkg) const { + if (mTypeOfBkgPdf == NoBkg) { + bkg = 0.; + errBkg = 0.; + return; + } bkg = mRooNBkg->getVal() * mIntegralBkg; errBkg = mRooNBkg->getError() * mIntegralBkg; } @@ -691,8 +769,7 @@ void HFInvMassFitter::calculateSignificance(double& significance, double& errSig // estimate Signal void HFInvMassFitter::checkForSignal(double& estimatedSignal) { - double const minForSgn = mMass - 4 * mSigmaSgn; - double const maxForSgn = mMass + 4 * mSigmaSgn; + auto const [minForSgn, maxForSgn] = getRangesOfSignal(); int const binForMinSgn = mHistoInvMass->FindBin(minForSgn); int const binForMaxSgn = mHistoInvMass->FindBin(maxForSgn); @@ -705,6 +782,19 @@ void HFInvMassFitter::checkForSignal(double& estimatedSignal) estimatedSignal = sum - bkg; } +// Estimate ranges where signal is located to be used in countSignal() and checkForSignal() +// It is mu +- n*sigma for Gaussian, and the entire mInv histogram for DoubleSidedCrystalBall (due to heavy tails) +std::pair HFInvMassFitter::getRangesOfSignal() const +{ + if (mTypeOfSgnPdf == DoubleSidedCrystalBall) { + return std::make_pair(mMinMass, mMaxMass); + } else { + const double mean = mRooMeanSgn->getVal(); + const double sigma = mRooSecSigmaSgn->getVal(); + return std::make_pair(mean - mNSigmaForSgn * sigma, mean + mNSigmaForSgn * sigma); + } +} + // Create Background Fit Function RooAbsPdf* HFInvMassFitter::createBackgroundFitFunction(RooWorkspace* workspace) const { @@ -751,6 +841,16 @@ RooAbsPdf* HFInvMassFitter::createSignalFitFunction(RooWorkspace* workspace) mRooSecSigmaSgn = workspace->var("sigmaSec"); mRooMeanSgn = workspace->var("meanSec"); } break; + case DoubleSidedCrystalBall: { + sgnPdf = workspace->pdf("sgnFuncDSCB"); + mRooSigmaSgn = workspace->var("sigma"); + mRooSecSigmaSgn = workspace->var("sigma"); + mRooMeanSgn = workspace->var("mean"); + mRooDscbAlphaL = workspace->var("alphaL"); + mRooDscbNL = workspace->var("nL"); + mRooDscbAlphaR = workspace->var("alphaR"); + mRooDscbNR = workspace->var("nR"); + } break; default: break; } @@ -991,3 +1091,61 @@ void HFInvMassFitter::setTemplateReflections(TH1* histoRefl) mHistoTemplateRefl = histoRefl; mHistoTemplateRefl->SetName("mHistoTemplateRefl"); } + +void HFInvMassFitter::setFixDscbAlphaL(double alphaL) +{ + if (mTypeOfSgnPdf != DoubleSidedCrystalBall) { + printf("Fit type should be DSCB!\n"); + } + mFixedDscbTailParams = kTRUE; + mDscbAlphaLInitialValue = alphaL; +} +void HFInvMassFitter::setFixDscbAlphaR(double alphaR) +{ + if (mTypeOfSgnPdf != DoubleSidedCrystalBall) { + printf("Fit type should be DSCB!\n"); + } + mFixedDscbTailParams = kTRUE; + mDscbAlphaRInitialValue = alphaR; +} +void HFInvMassFitter::setFixDscbNL(double nL) +{ + if (mTypeOfSgnPdf != DoubleSidedCrystalBall) { + printf("Fit type should be DSCB!\n"); + } + mFixedDscbTailParams = kTRUE; + mDscbNLInitialValue = nL; +} +void HFInvMassFitter::setFixDscbNR(double nR) +{ + if (mTypeOfSgnPdf != DoubleSidedCrystalBall) { + printf("Fit type should be DSCB!\n"); + } + mFixedDscbTailParams = kTRUE; + mDscbNRInitialValue = nR; +} + +double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& parameterRanges) const +{ + constexpr double DefaultSigmaFraction{10.}; + constexpr int MaximalNumberOfIterations{20}; + + if (mRandomSeed < 0) { + return parameterRanges.initial; + } + const auto sigma = parameterRanges.sigma < 0 ? (parameterRanges.upper - parameterRanges.lower) / DefaultSigmaFraction : parameterRanges.sigma; + + double result; + int nIter{0}; + do { + result = mRandomGen->Gaus(parameterRanges.initial, sigma); + ++nIter; + if (nIter > MaximalNumberOfIterations) { + char errorMessage[200]; + std::snprintf(errorMessage, sizeof(errorMessage), "randomizeInitialParameter() - long while loop with lower = %f upper = %f initial = %f sigma = %f\n", parameterRanges.lower, parameterRanges.upper, parameterRanges.initial, sigma); + throw std::runtime_error(errorMessage); + } + } while (result < parameterRanges.lower || result > parameterRanges.upper); + + return result; +} diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 28ee858bba6..6841cf15a55 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -28,6 +28,7 @@ #include #include #include +#include #include #include @@ -35,10 +36,19 @@ #include #include +#include #include class HFInvMassFitter : public TNamed { + private: + struct ParameterRanges { + double lower{}; + double upper{}; + double initial{}; + double sigma{-999.}; + }; + public: enum TypeOfBkgPdf { Expo = 0, @@ -56,6 +66,7 @@ class HFInvMassFitter : public TNamed DoubleGaus = 1, DoubleGausSigmaRatioPar = 2, GausSec = 3, + DoubleSidedCrystalBall = 4, NTypesOfSgnPdf }; enum TypeOfReflPdf { @@ -67,7 +78,7 @@ class HFInvMassFitter : public TNamed }; std::array namesOfReflPdf{"reflFuncGaus", "reflFuncDoubleGaus", "reflFuncPoly3", "reflFuncPoly6"}; HFInvMassFitter() = delete; - HFInvMassFitter(TH1* histoToFit, double minValue, double maxValue, int fitTypeBkg = Expo, int fitTypeSgn = SingleGaus); + HFInvMassFitter(TH1* histoToFit, double minValue, double maxValue, int fitTypeBkg = Expo, int fitTypeSgn = SingleGaus, int randomSeed = -1); ~HFInvMassFitter() override; void setHistogramForFit(TH1* histoToFit); void setUseLikelihoodFit() { mFitOption = "L,E"; } @@ -97,6 +108,23 @@ class HFInvMassFitter : public TNamed void setFixRatioToGausSigma(double sigmaFrac); void setFixSignalYield(double yield) { mFixedRawYield = yield; } void setNumberOfSigmaForSidebands(double numberOfSigma) { mNSigmaForSidebands = numberOfSigma; } + void setNumberOfSigmaForSignal(double numberOfSigma) { mNSigmaForSgn = numberOfSigma; } + void setFixDscbAlphaL(double alphaL); + void setFixDscbAlphaR(double alphaR); + void setFixDscbNL(double nL); + void setFixDscbNR(double nR); + void setDscbAlphaLInitialValue(double value) { mDscbAlphaLInitialValue = value; } + void setDscbAlphaLLowLimit(double value) { mDscbAlphaLLowLimit = value; } + void setDscbAlphaLUpLimit(double value) { mDscbAlphaLUpLimit = value; } + void setDscbAlphaRInitialValue(double value) { mDscbAlphaRInitialValue = value; } + void setDscbAlphaRLowLimit(double value) { mDscbAlphaRLowLimit = value; } + void setDscbAlphaRUpLimit(double value) { mDscbAlphaRUpLimit = value; } + void setDscbNLInitialValue(double value) { mDscbNLInitialValue = value; } + void setDscbNLLowLimit(double value) { mDscbNLLowLimit = value; } + void setDscbNLUpLimit(double value) { mDscbNLUpLimit = value; } + void setDscbNRInitialValue(double value) { mDscbNRInitialValue = value; } + void setDscbNRLowLimit(double value) { mDscbNRLowLimit = value; } + void setDscbNRUpLimit(double value) { mDscbNRUpLimit = value; } void plotBkg(RooAbsPdf* mFunc, Color_t color = kRed); void plotRefl(RooAbsPdf* mFunc); void setReflFuncFixed(); @@ -120,6 +148,14 @@ class HFInvMassFitter : public TNamed [[nodiscard]] double getMeanUncertainty() const { return mRooMeanSgn->getError(); } [[nodiscard]] double getSigma() const { return mRooSigmaSgn->getVal(); } [[nodiscard]] double getSigmaUncertainty() const { return mRooSigmaSgn->getError(); } + [[nodiscard]] double getDscbAlphaL() const { return mRooDscbAlphaL ? mRooDscbAlphaL->getVal() : 0.; } + [[nodiscard]] double getDscbAlphaR() const { return mRooDscbAlphaR ? mRooDscbAlphaR->getVal() : 0.; } + [[nodiscard]] double getDscbNL() const { return mRooDscbNL ? mRooDscbNL->getVal() : 0.; } + [[nodiscard]] double getDscbNR() const { return mRooDscbNR ? mRooDscbNR->getVal() : 0.; } + [[nodiscard]] double getDscbAlphaLUncertainty() const { return mRooDscbAlphaL ? mRooDscbAlphaL->getError() : 0.; } + [[nodiscard]] double getDscbAlphaRUncertainty() const { return mRooDscbAlphaR ? mRooDscbAlphaR->getError() : 0.; } + [[nodiscard]] double getDscbNLUncertainty() const { return mRooDscbNL ? mRooDscbNL->getError() : 0.; } + [[nodiscard]] double getDscbNRUncertainty() const { return mRooDscbNR ? mRooDscbNR->getError() : 0.; } [[nodiscard]] double getSecSigma() const { return mRooSecSigmaSgn->getVal(); } [[nodiscard]] double getSecSigmaUncertainty() const { return mRooSecSigmaSgn->getError(); } [[nodiscard]] double getFracDoubleGaus() const { return mRooFracDoubleGaus->getVal(); } @@ -141,6 +177,8 @@ class HFInvMassFitter : public TNamed HFInvMassFitter& operator=(const HFInvMassFitter& source); void fillWorkspace(RooWorkspace& w) const; void highlightPeakRegion(const RooPlot* plot, Color_t color = kGray + 1, Width_t width = 1, Style_t style = 2) const; + [[nodiscard]] double randomizeInitialParameter(const ParameterRanges& parameterRanges) const; + [[nodiscard]] std::pair getRangesOfSignal() const; TH1* mHistoInvMass; // histogram to fit std::string mFitOption; @@ -168,6 +206,7 @@ class HFInvMassFitter : public TNamed bool mFixedSigma; /// fix sigma or not bool mFixedSigmaDoubleGaus; /// fix sigma of 2gaussian or not bool mBoundSigma; /// set bound sigma or not + bool mFixedDscbTailParams; /// switch for fix double sided Crystal Ball tail parameters double mSigmaValue; /// value of sigma double mParamSgn; /// +/- range variation of bound Sigma of gaussian in % double mFracDoubleGaus; /// initialization for fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar @@ -188,6 +227,18 @@ class HFInvMassFitter : public TNamed double mChiSquareOverNdfTotal; /// chi2/ndf of the total fit double mChiSquareOverNdfBkg; /// chi2/ndf of the background (sidebands) pre-fit bool mFixReflOverSgn; /// switch for fix refl/signal + double mDscbAlphaLInitialValue; /// double sided Crystal Ball alpha left initial value + double mDscbAlphaLLowLimit; /// double sided Crystal Ball alpha left lower limit + double mDscbAlphaLUpLimit; /// double sided Crystal Ball alpha left upper limit + double mDscbAlphaRInitialValue; /// double sided Crystal Ball alpha right initial value + double mDscbAlphaRLowLimit; /// double sided Crystal Ball alpha right lower limit + double mDscbAlphaRUpLimit; /// double sided Crystal Ball alpha right upper limit + double mDscbNLInitialValue; /// double sided Crystal Ball n left initial value + double mDscbNLLowLimit; /// double sided Crystal Ball n left lower limit + double mDscbNLUpLimit; /// double sided Crystal Ball n left upper limit + double mDscbNRInitialValue; /// double sided Crystal Ball n right initial value + double mDscbNRLowLimit; /// double sided Crystal Ball n right lower limit + double mDscbNRUpLimit; /// double sided Crystal Ball n right upper limit RooRealVar* mRooMeanSgn; /// mean for gaussian of signal RooRealVar* mRooSigmaSgn; /// sigma for gaussian of signal RooRealVar* mRooSecSigmaSgn; /// second sigma for composite gaussian of signal @@ -198,20 +249,26 @@ class HFInvMassFitter : public TNamed RooRealVar* mRooNSgn; /// total Signal fit function integral RooRealVar* mRooNBkg; /// total background fit function integral RooRealVar* mRooNRefl; /// total reflection fit function integral + RooRealVar* mRooDscbAlphaL; /// double sided Crystal Ball alpha left + RooRealVar* mRooDscbAlphaR; /// double sided Crystal Ball alpha right + RooRealVar* mRooDscbNL; /// double sided Crystal Ball n left + RooRealVar* mRooDscbNR; /// double sided Crystal Ball n right RooAbsPdf* mTotalPdf; /// total fit function RooPlot* mInvMassFrame; /// frame of mass RooPlot* mReflFrame; /// reflection frame RooPlot* mReflOnlyFrame; /// reflection frame plot on reflection only RooPlot* mResidualFrame; /// residual frame + RooHist* mResidualHist; /// residual histogram RooPlot* mRatioFrame; /// fit/data ratio frame - RooPlot* mResidualFrameForCalculation; - RooWorkspace* mWorkspace; /// workspace - double mIntegralHisto; /// integral of histogram to fit - double mIntegralBkg; /// integral of background fit function - double mIntegralSgn; /// integral of signal fit function - TH1* mHistoTemplateRefl; /// reflection histogram - bool mDrawBgPrefit; /// draw background after fitting the sidebands - bool mHighlightPeakRegion; /// draw vertical lines showing the peak region (usually +- 3 sigma) + RooWorkspace* mWorkspace; /// workspace + double mIntegralHisto; /// integral of histogram to fit + double mIntegralBkg; /// integral of background fit function + double mIntegralSgn; /// integral of signal fit function + TH1* mHistoTemplateRefl; /// reflection histogram + bool mDrawBgPrefit; /// draw background after fitting the sidebands + bool mHighlightPeakRegion; /// draw vertical lines showing the peak region (usually +- 3 sigma) + int mRandomSeed; /// seed for random engine for fit's initial parameters randomization + TRandom3* mRandomGen; /// engine for fit's initial parameters randomization ClassDefOverride(HFInvMassFitter, 1); }; diff --git a/PWGHF/D2H/Macros/config_massfitter.json b/PWGHF/D2H/Macros/config_massfitter.json index 25dc9d50bcc..0baf94cabad 100644 --- a/PWGHF/D2H/Macros/config_massfitter.json +++ b/PWGHF/D2H/Macros/config_massfitter.json @@ -101,21 +101,24 @@ 0.9, 1.6 ], - "MassMin": [ + "MassMin": 2.12, + "_MassMin": [ 2.12, 2.12, 2.12, 2.12, 2.12 ], - "MassMax": [ + "MassMax": 2.42, + "_MassMax": [ 2.42, 2.42, 2.42, 2.42, 2.42 ], - "Rebin": [ + "Rebin": 4, + "_Rebin": [ 4, 4, 4, @@ -128,7 +131,8 @@ "true: likelihood fit", "false: chi2 fit" ], - "BkgFunc": [ + "BkgFunc": 2, + "_BkgFunc": [ 2, 2, 2, @@ -144,7 +148,8 @@ "5 for Poly3", "6 for NoBkg" ], - "SgnFunc": [ + "SgnFunc": 0, + "_SgnFunc": [ 0, 0, 0, @@ -155,8 +160,12 @@ "0 for SingleGaus", "1 for DoubleGaus", "2 for DoubleGausSigmaRatioPar", - "3 for GausSec" + "3 for GausSec", + "4 for DoubleSidedCrystalBall" ], + "NSigmaForSideband": 3, + "NSigmaForSignal": 3, "DrawBgPrefit": true, - "HighlightPeakRegion": true + "HighlightPeakRegion": true, + "RandomSeed": -1 } diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 5937d12e0ba..3b05aa6f9ed 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -37,6 +37,7 @@ #include #include +#include #include #include #include @@ -50,6 +51,8 @@ using namespace rapidjson; +constexpr int UndefValueInt{-999}; + void runMassFitter(const std::string& configFileName = "config_massfitter.json"); TFile* openFileWithNullptrCheck(const std::string& fileName, const std::string& option = "read"); @@ -66,7 +69,12 @@ T readJsonField(const Document& config, const std::string& fieldName); template void readJsonVector(std::vector& vec, const Document& config, const std::string& fieldName, bool isRequired = false); -void divideCanvas(TCanvas* c, int nSliceVarBins); +template +void readJsonVectorFlexible(std::vector& vec, const Document& config, int nHistograms, const std::string& fieldName, bool isRequired = false); + +void readJsonVectorFromHisto(std::vector& vec, const Document& config, const std::string& fileNameFieldName, const std::string& histoNameFieldName); + +void divideCanvas(TCanvas* c, int nHistograms); void setHistoStyle(TH1* histo, Color_t color = kBlack, Size_t markerSize = 1); @@ -84,6 +92,10 @@ void runMassFitter(const std::string& configFileName) config.ParseStream(is); fclose(configFile); + if (config.HasParseError()) { + throw std::runtime_error("ERROR: Parsing the configuration json file failed. Check the config for correct formatting"); + } + bool const isMc = readJsonField(config, "IsMC"); bool const writeSignalPar = readJsonField(config, "WriteSignalPar", true); std::string const particleName = readJsonField(config, "Particle"); @@ -95,9 +107,11 @@ void runMassFitter(const std::string& configFileName) std::vector inputHistoName; std::vector promptHistoName; std::vector fdHistoName; + std::vector signalHistoName; std::vector reflHistoName; std::vector promptSecPeakHistoName; std::vector fdSecPeakHistoName; + std::vector signalSecPeakHistoName; std::vector sliceVarMin; std::vector sliceVarMax; std::vector massMin; @@ -109,29 +123,48 @@ void runMassFitter(const std::string& configFileName) std::vector nRebin; std::vector bkgFunc; std::vector sgnFunc; - - readJsonVector(inputHistoName, config, "InputHistoName", true); + std::vector dscbAlphaLInitial; + std::vector dscbAlphaLLower; + std::vector dscbAlphaLUpper; + std::vector dscbAlphaRInitial; + std::vector dscbAlphaRLower; + std::vector dscbAlphaRUpper; + std::vector dscbNLInitial; + std::vector dscbNLLower; + std::vector dscbNLUpper; + std::vector dscbNRInitial; + std::vector dscbNRLower; + std::vector dscbNRUpper; + + readJsonVector(inputHistoName, config, "InputHistoName"); readJsonVector(promptHistoName, config, "PromptHistoName"); readJsonVector(fdHistoName, config, "FDHistoName"); + readJsonVector(signalHistoName, config, "SignalHistoName"); + const std::array possibleInputHistogramSizes{inputHistoName.size(), promptHistoName.size(), fdHistoName.size(), signalHistoName.size()}; + const int nHistograms = static_cast(*std::max_element(possibleInputHistogramSizes.begin(), possibleInputHistogramSizes.end())); + readJsonVector(reflHistoName, config, "ReflHistoName"); readJsonVector(promptSecPeakHistoName, config, "PromptSecPeakHistoName"); readJsonVector(fdSecPeakHistoName, config, "FDSecPeakHistoName"); + readJsonVector(signalSecPeakHistoName, config, "SignalSecPeakHistoName"); const bool fixMean = readJsonField(config, "FixMean", false); const std::string meanFile = readJsonField(config, "MeanFile", ""); - readJsonVector(fixMeanManual, config, "FixMeanManual"); + readJsonVectorFlexible(fixMeanManual, config, nHistograms, "FixMeanManual"); const bool fixSigma = readJsonField(config, "FixSigma", false); const std::string sigmaFile = readJsonField(config, "SigmaFile", ""); - readJsonVector(fixSigmaManual, config, "FixSigmaManual"); + readJsonVectorFlexible(fixSigmaManual, config, nHistograms, "FixSigmaManual"); const bool fixSecondSigma = readJsonField(config, "FixSecondSigma", false); const std::string secondSigmaFile = readJsonField(config, "SecondSigmaFile", ""); - readJsonVector(fixSecondSigmaManual, config, "FixSecondSigmaManual"); + readJsonVectorFlexible(fixSecondSigmaManual, config, nHistograms, "FixSecondSigmaManual"); const bool fixFracDoubleGaus = readJsonField(config, "FixFracDoubleGaus", false); const std::string fracDoubleGausFile = readJsonField(config, "FracDoubleGausFile", ""); - readJsonVector(fixFracDoubleGausManual, config, "FixFracDoubleGausManual"); + readJsonVectorFlexible(fixFracDoubleGausManual, config, nHistograms, "FixFracDoubleGausManual"); + + const bool fixDscbTailParams = readJsonField(config, "FixDscbTailParams", false); const TString sliceVarName = readJsonField(config, "SliceVarName"); const TString sliceVarUnit = readJsonField(config, "SliceVarUnit"); @@ -139,26 +172,53 @@ void runMassFitter(const std::string& configFileName) readJsonVector(sliceVarMin, config, "SliceVarMin", true); readJsonVector(sliceVarMax, config, "SliceVarMax", true); - readJsonVector(massMin, config, "MassMin", true); - readJsonVector(massMax, config, "MassMax", true); + readJsonVectorFlexible(massMin, config, nHistograms, "MassMin", true); + readJsonVectorFlexible(massMax, config, nHistograms, "MassMax", true); - readJsonVector(nRebin, config, "Rebin", true); + readJsonVectorFlexible(nRebin, config, nHistograms, "Rebin", true); bool const includeSecPeak = readJsonField(config, "InclSecPeak", false); bool const useLikelihood = readJsonField(config, "UseLikelihood"); - readJsonVector(bkgFunc, config, "BkgFunc", true); - readJsonVector(sgnFunc, config, "SgnFunc", true); + readJsonVectorFlexible(bkgFunc, config, nHistograms, "BkgFunc", true); + readJsonVectorFlexible(sgnFunc, config, nHistograms, "SgnFunc", true); const bool enableRefl = readJsonField(config, "EnableRefl", false); const bool drawBgPrefit = readJsonField(config, "DrawBgPrefit", true); const bool highlightPeakRegion = readJsonField(config, "HighlightPeakRegion", true); - - const int nSliceVarBins = static_cast(sliceVarMin.size()); - std::vector sliceVarLimits(nSliceVarBins + 1); + const int randomSeed = readJsonField(config, "RandomSeed", -1); + const double nSigmaForSideband = readJsonField(config, "NSigmaForSideband", 3.); + const double nSigmaForSignal = readJsonField(config, "NSigmaForSignal", 3.); + + readJsonVector(dscbAlphaLInitial, config, "DscbAlphaLInitial"); + readJsonVector(dscbAlphaLLower, config, "DscbAlphaLLower"); + readJsonVector(dscbAlphaLUpper, config, "DscbAlphaLUpper"); + readJsonVector(dscbAlphaRInitial, config, "DscbAlphaRInitial"); + readJsonVector(dscbAlphaRLower, config, "DscbAlphaRLower"); + readJsonVector(dscbAlphaRUpper, config, "DscbAlphaRUpper"); + readJsonVector(dscbNLInitial, config, "DscbNLInitial"); + readJsonVector(dscbNLLower, config, "DscbNLLower"); + readJsonVector(dscbNLUpper, config, "DscbNLUpper"); + readJsonVector(dscbNRInitial, config, "DscbNRInitial"); + readJsonVector(dscbNRLower, config, "DscbNRLower"); + readJsonVector(dscbNRUpper, config, "DscbNRUpper"); + readJsonVectorFromHisto(dscbAlphaLInitial, config, "DscbParametersFile", "DscbAlphaLInitialHisto"); + readJsonVectorFromHisto(dscbAlphaLLower, config, "DscbParametersFile", "DscbAlphaLLowerHisto"); + readJsonVectorFromHisto(dscbAlphaLUpper, config, "DscbParametersFile", "DscbAlphaLUpperHisto"); + readJsonVectorFromHisto(dscbAlphaRInitial, config, "DscbParametersFile", "DscbAlphaRInitialHisto"); + readJsonVectorFromHisto(dscbAlphaRLower, config, "DscbParametersFile", "DscbAlphaRLowerHisto"); + readJsonVectorFromHisto(dscbAlphaRUpper, config, "DscbParametersFile", "DscbAlphaRUpperHisto"); + readJsonVectorFromHisto(dscbNLInitial, config, "DscbParametersFile", "DscbNLInitialHisto"); + readJsonVectorFromHisto(dscbNLLower, config, "DscbParametersFile", "DscbNLLowerHisto"); + readJsonVectorFromHisto(dscbNLUpper, config, "DscbParametersFile", "DscbNLUpperHisto"); + readJsonVectorFromHisto(dscbNRInitial, config, "DscbParametersFile", "DscbNRInitialHisto"); + readJsonVectorFromHisto(dscbNRLower, config, "DscbParametersFile", "DscbNRLowerHisto"); + readJsonVectorFromHisto(dscbNRUpper, config, "DscbParametersFile", "DscbNRUpperHisto"); + + std::vector sliceVarLimits(nHistograms + 1); auto checkVectorSize = [&](const auto& vec, const std::string& name = "", const bool isEmptyOk = false) { - if (vec.size() != static_cast(nSliceVarBins)) { + if (vec.size() != static_cast(nHistograms)) { if (isEmptyOk && vec.empty()) { return; } @@ -169,9 +229,11 @@ void runMassFitter(const std::string& configFileName) checkVectorSize(inputHistoName, "inputHistoName", true); checkVectorSize(promptHistoName, "promptHistoName", true); checkVectorSize(fdHistoName, "fdHistoName", true); + checkVectorSize(signalHistoName, "signalHistoName", true); checkVectorSize(reflHistoName, "reflHistoName", true); checkVectorSize(promptSecPeakHistoName, "promptSecPeakHistoName", true); checkVectorSize(fdSecPeakHistoName, "fdSecPeakHistoName", true); + checkVectorSize(signalSecPeakHistoName, "signalSecPeakHistoName", true); checkVectorSize(sliceVarMin, "sliceVarMin"); checkVectorSize(sliceVarMax, "sliceVarMax"); checkVectorSize(massMin, "massMin"); @@ -183,8 +245,36 @@ void runMassFitter(const std::string& configFileName) checkVectorSize(nRebin, "nRebin"); checkVectorSize(bkgFunc, "bkgFunc"); checkVectorSize(sgnFunc, "sgnFunc"); + checkVectorSize(dscbAlphaLInitial, "dscbAlphaLInitial", true); + checkVectorSize(dscbAlphaLLower, "dscbAlphaLLower", true); + checkVectorSize(dscbAlphaLUpper, "dscbAlphaLUpper", true); + checkVectorSize(dscbAlphaRInitial, "dscbAlphaRInitial", true); + checkVectorSize(dscbAlphaRLower, "dscbAlphaRLower", true); + checkVectorSize(dscbAlphaRUpper, "dscbAlphaRUpper", true); + checkVectorSize(dscbNLInitial, "dscbNLInitial", true); + checkVectorSize(dscbNLLower, "dscbNLLower", true); + checkVectorSize(dscbNLUpper, "dscbNLUpper", true); + checkVectorSize(dscbNRInitial, "dscbNRInitial", true); + checkVectorSize(dscbNRLower, "dscbNRLower", true); + checkVectorSize(dscbNRUpper, "dscbNRUpper", true); + + auto checkVectorSizeMcHistograms = [](const auto& vecSignal, const auto& vecPrompt, const auto& vecFd) { + const auto signalSize = vecSignal.size(); + const auto promptSize = vecPrompt.size(); + const auto fdSize = vecFd.size(); + if (!((signalSize > 0 && promptSize == 0 && fdSize == 0) || (signalSize == 0 && promptSize > 0 && fdSize > 0))) { + throw std::runtime_error("ERROR: either signal histogram must be provided or both prompt and fd, but not all three. Exit"); + } + }; + + if ((!isMc && enableRefl) || isMc) { + checkVectorSizeMcHistograms(signalHistoName, promptHistoName, fdHistoName); + } + if (isMc && includeSecPeak) { + checkVectorSizeMcHistograms(signalSecPeakHistoName, promptSecPeakHistoName, fdSecPeakHistoName); + } - for (int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) { + for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { sliceVarLimits[iSliceVar] = sliceVarMin[iSliceVar]; if (bkgFunc[iSliceVar] < 0 || bkgFunc[iSliceVar] >= HFInvMassFitter::NTypesOfBkgPdf) { @@ -197,7 +287,7 @@ void runMassFitter(const std::string& configFileName) throw std::runtime_error("ERROR: only SingleGaus, DoubleGaus and DoubleGausSigmaRatioPar signal functions supported! Exit"); } } - sliceVarLimits[nSliceVarBins] = sliceVarMax[nSliceVarBins - 1]; + sliceVarLimits[nHistograms] = sliceVarMax[nHistograms - 1]; struct DecayInfo { std::string decayProducts; @@ -227,24 +317,36 @@ void runMassFitter(const std::string& configFileName) TFile* inputFileRefl = enableRefl ? openFileWithNullptrCheck(reflFileName) : nullptr; - std::vector hMassSgn(nSliceVarBins); - std::vector hMassRefl(nSliceVarBins); - std::vector hMass(nSliceVarBins); + std::vector hMassSgn(nHistograms); + std::vector hMassRefl(nHistograms); + std::vector hMass(nHistograms); - for (int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) { + for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { if (!isMc) { hMass[iSliceVar] = getObjectWithNullPtrCheck(inputFile, inputHistoName[iSliceVar]); if (enableRefl) { hMassRefl[iSliceVar] = getObjectWithNullPtrCheck(inputFileRefl, reflHistoName[iSliceVar]); - hMassSgn[iSliceVar] = getObjectWithNullPtrCheck(inputFileRefl, fdHistoName[iSliceVar]); - hMassSgn[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFileRefl, promptHistoName[iSliceVar])); + if (!signalHistoName.empty()) { + hMassSgn[iSliceVar] = getObjectWithNullPtrCheck(inputFileRefl, signalHistoName[iSliceVar]); + } else { + hMassSgn[iSliceVar] = getObjectWithNullPtrCheck(inputFileRefl, fdHistoName[iSliceVar]); + hMassSgn[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFileRefl, promptHistoName[iSliceVar])); + } } } else { - hMass[iSliceVar] = getObjectWithNullPtrCheck(inputFile, promptHistoName[iSliceVar]); - hMass[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFile, fdHistoName[iSliceVar])); + if (!signalHistoName.empty()) { + hMass[iSliceVar] = getObjectWithNullPtrCheck(inputFile, signalHistoName[iSliceVar]); + } else { + hMass[iSliceVar] = getObjectWithNullPtrCheck(inputFile, promptHistoName[iSliceVar]); + hMass[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFile, fdHistoName[iSliceVar])); + } if (includeSecPeak) { - hMass[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFile, promptSecPeakHistoName[iSliceVar])); - hMass[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFile, fdSecPeakHistoName[iSliceVar])); + if (!signalHistoName.empty()) { + hMass[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFile, signalSecPeakHistoName[iSliceVar])); + } else { + hMass[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFile, promptSecPeakHistoName[iSliceVar])); + hMass[iSliceVar]->Add(getObjectWithNullPtrCheck(inputFile, fdSecPeakHistoName[iSliceVar])); + } } } hMass[iSliceVar]->SetDirectory(nullptr); @@ -255,18 +357,22 @@ void runMassFitter(const std::string& configFileName) } // define output histos - auto* hRawYieldsSignal = new TH1D("hRawYieldsSignal", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsSignalCounted = new TH1D("hRawYieldsSignalCounted", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield via bin count", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsBkg = new TH1D("hRawYieldsBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");Background (3#sigma)", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsSgnOverBkg = new TH1D("hRawYieldsSgnOverBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");S/B (3#sigma)", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsSignificance = new TH1D("hRawYieldsSignificance", ";" + sliceVarName + "(" + sliceVarUnit + ");significance (3#sigma)", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsChiSquareBkg = new TH1D("hRawYieldsChiSquareBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");#chi^{2}/#it{ndf}", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsChiSquareTotal = new TH1D("hRawYieldsChiSquareTotal", ";" + sliceVarName + "(" + sliceVarUnit + ");#chi^{2}/#it{ndf}", nSliceVarBins, sliceVarLimits.data()); - auto* hReflectionOverSignal = new TH1D("hReflectionOverSignal", ";" + sliceVarName + "(" + sliceVarUnit + ");Refl/Signal", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsMean = new TH1D("hRawYieldsMean", ";" + sliceVarName + "(" + sliceVarUnit + ");mean (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsSigma = new TH1D("hRawYieldsSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsSecSigma = new TH1D("hRawYieldsSecSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsFracDoubleGaus = new TH1D("hRawYieldsFracDoubleGaus", ";" + sliceVarName + "(" + sliceVarUnit + ");fraction of double gaussian", nSliceVarBins, sliceVarLimits.data()); + auto* hRawYieldsSignal = new TH1D("hRawYieldsSignal", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsSignalCounted = new TH1D("hRawYieldsSignalCounted", ";" + sliceVarName + "(" + sliceVarUnit + ");raw yield via bin count", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsBkg = new TH1D("hRawYieldsBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");Background (3#sigma)", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsSgnOverBkg = new TH1D("hRawYieldsSgnOverBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");S/B (3#sigma)", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsSignificance = new TH1D("hRawYieldsSignificance", ";" + sliceVarName + "(" + sliceVarUnit + ");significance (3#sigma)", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsChiSquareBkg = new TH1D("hRawYieldsChiSquareBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");#chi^{2}/#it{ndf}", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsChiSquareTotal = new TH1D("hRawYieldsChiSquareTotal", ";" + sliceVarName + "(" + sliceVarUnit + ");#chi^{2}/#it{ndf}", nHistograms, sliceVarLimits.data()); + auto* hReflectionOverSignal = new TH1D("hReflectionOverSignal", ";" + sliceVarName + "(" + sliceVarUnit + ");Refl/Signal", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsMean = new TH1D("hRawYieldsMean", ";" + sliceVarName + "(" + sliceVarUnit + ");mean (GeV/#it{c}^{2})", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsSigma = new TH1D("hRawYieldsSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsSecSigma = new TH1D("hRawYieldsSecSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsFracDoubleGaus = new TH1D("hRawYieldsFracDoubleGaus", ";" + sliceVarName + "(" + sliceVarUnit + ");fraction of double gaussian", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsDscbAlphaL = new TH1D("hRawYieldsDscbAlphaL", ";" + sliceVarName + "(" + sliceVarUnit + ");#alpha_{L}", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsDscbAlphaR = new TH1D("hRawYieldsDscbAlphaR", ";" + sliceVarName + "(" + sliceVarUnit + ");#alpha_{R}", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsDscbNL = new TH1D("hRawYieldsDscbNL", ";" + sliceVarName + "(" + sliceVarUnit + ");n_{L}", nHistograms, sliceVarLimits.data()); + auto* hRawYieldsDscbNR = new TH1D("hRawYieldsDscbNR", ";" + sliceVarName + "(" + sliceVarUnit + ");n_{R}", nHistograms, sliceVarLimits.data()); enum { ConfigMassMin = 1, @@ -275,10 +381,11 @@ void runMassFitter(const std::string& configFileName) ConfigFixSigma, ConfigBkgFunc, ConfigSgnFunc, + ConfigRandomSeed, NConfigsToSave }; - auto* hFitConfig = new TH2F("hfitConfig", "Fit Configurations", NConfigsToSave - 1, 0, NConfigsToSave - 1, nSliceVarBins, sliceVarLimits.data()); - const char* hFitConfigXLabel[NConfigsToSave - 1] = {"mass min", "mass max", "rebin num", "fix sigma", "bkg func", "sgn func"}; + auto* hFitConfig = new TH2F("hfitConfig", "Fit Configurations", NConfigsToSave - 1, 0, NConfigsToSave - 1, nHistograms, sliceVarLimits.data()); + const char* hFitConfigXLabel[NConfigsToSave - 1] = {"mass min", "mass max", "rebin num", "fix sigma", "bkg func", "sgn func", "rnd seed"}; hFitConfig->SetStats(false); for (int i = 0; i < NConfigsToSave - 1; i++) { hFitConfig->GetXaxis()->SetBinLabel(i + 1, hFitConfigXLabel[i]); @@ -299,8 +406,12 @@ void runMassFitter(const std::string& configFileName) setHistoStyle(hRawYieldsSigma); setHistoStyle(hRawYieldsSecSigma); setHistoStyle(hRawYieldsFracDoubleGaus); + setHistoStyle(hRawYieldsDscbAlphaL); + setHistoStyle(hRawYieldsDscbAlphaR); + setHistoStyle(hRawYieldsDscbNL); + setHistoStyle(hRawYieldsDscbNR); - auto getHistToFix = [&nSliceVarBins](bool const& isFix, std::vector const& fixManual, std::string const& fixFileName, std::string const& var) -> TH1* { + auto getHistToFix = [&nHistograms](bool const& isFix, std::vector const& fixManual, std::string const& fixFileName, std::string const& var) -> TH1* { TH1* histToFix = nullptr; if (isFix) { if (fixManual.empty()) { @@ -308,7 +419,7 @@ void runMassFitter(const std::string& configFileName) const std::string histName = "hRawYields" + var; histToFix = getObjectWithNullPtrCheck(fixInputFile, histName); histToFix->SetDirectory(nullptr); - if (histToFix->GetNbinsX() != nSliceVarBins) { + if (histToFix->GetNbinsX() != nHistograms) { throw std::runtime_error("Different number of bins for this analysis and histo for fixed " + var); } fixInputFile->Close(); @@ -323,19 +434,19 @@ void runMassFitter(const std::string& configFileName) TH1* hFracDoubleGausToFix = getHistToFix(fixFracDoubleGaus, fixFracDoubleGausManual, fracDoubleGausFile, "FracDoubleGaus"); int canvasSize[2] = {1920, 1080}; - if (nSliceVarBins == 1) { + if (nHistograms == 1) { canvasSize[0] = 500; canvasSize[1] = 500; } int constexpr NCanvasesMax = 20; // do not put more than 20 bins per canvas to make them visible - const int nCanvases = std::ceil(static_cast(nSliceVarBins) / NCanvasesMax); + const int nCanvases = std::ceil(static_cast(nHistograms) / NCanvasesMax); std::vector canvasMass(nCanvases); std::vector canvasResiduals(nCanvases); std::vector canvasRatio(nCanvases); std::vector canvasRefl(nCanvases); for (int iCanvas = 0; iCanvas < nCanvases; iCanvas++) { - const int nPads = (nCanvases == 1) ? nSliceVarBins : NCanvasesMax; + const int nPads = (nCanvases == 1) ? nHistograms : NCanvasesMax; canvasMass[iCanvas] = new TCanvas(Form("canvasMass%d", iCanvas), Form("canvasMass%d", iCanvas), canvasSize[0], canvasSize[1]); divideCanvas(canvasMass[iCanvas], nPads); @@ -351,7 +462,7 @@ void runMassFitter(const std::string& configFileName) } } - for (int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) { + for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { const int iCanvas = std::floor(static_cast(iSliceVar) / NCanvasesMax); hMass[iSliceVar]->Rebin(nRebin[iSliceVar]); @@ -369,8 +480,10 @@ void runMassFitter(const std::string& configFileName) double reflOverSgn = 0; - HFInvMassFitter* massFitter = new HFInvMassFitter(hMass[iSliceVar], massMin[iSliceVar], massMax[iSliceVar], bkgFunc[iSliceVar], sgnFunc[iSliceVar]); + HFInvMassFitter* massFitter = new HFInvMassFitter(hMass[iSliceVar], massMin[iSliceVar], massMax[iSliceVar], bkgFunc[iSliceVar], sgnFunc[iSliceVar], randomSeed); massFitter->setDrawBgPrefit(drawBgPrefit); + massFitter->setNumberOfSigmaForSidebands(nSigmaForSideband); + massFitter->setNumberOfSigmaForSignal(nSigmaForSignal); massFitter->setHighlightPeakRegion(highlightPeakRegion); massFitter->setInitialGaussianMean(massPDG); massFitter->setParticlePdgMass(massPDG); @@ -398,6 +511,10 @@ void runMassFitter(const std::string& configFileName) setFixedValue(fixSigma, fixSigmaManual, hSigmaToFix, std::bind(&HFInvMassFitter::setFixGaussianSigma, massFitter, std::placeholders::_1), "SIGMA"); setFixedValue(fixSecondSigma, fixSecondSigmaManual, hSecondSigmaToFix, std::bind(&HFInvMassFitter::setFixSecondGaussianSigma, massFitter, std::placeholders::_1), "SECOND SIGMA"); setFixedValue(fixFracDoubleGaus, fixFracDoubleGausManual, hFracDoubleGausToFix, std::bind(&HFInvMassFitter::setFixFrac2Gaus, massFitter, std::placeholders::_1), "FRAC DOUBLE GAUS"); + setFixedValue(fixDscbTailParams, dscbAlphaLInitial, nullptr, std::bind(&HFInvMassFitter::setFixDscbAlphaL, massFitter, std::placeholders::_1), "DSCB ALPHA LEFT"); + setFixedValue(fixDscbTailParams, dscbAlphaRInitial, nullptr, std::bind(&HFInvMassFitter::setFixDscbAlphaR, massFitter, std::placeholders::_1), "DSCB ALPHA RIGHT"); + setFixedValue(fixDscbTailParams, dscbNLInitial, nullptr, std::bind(&HFInvMassFitter::setFixDscbNL, massFitter, std::placeholders::_1), "DSCB N LEFT"); + setFixedValue(fixDscbTailParams, dscbNRInitial, nullptr, std::bind(&HFInvMassFitter::setFixDscbNR, massFitter, std::placeholders::_1), "DSCB N RIGHT"); if (!isMc && enableRefl) { reflOverSgn = hMassSgn[iSliceVar]->Integral(hMassSgn[iSliceVar]->FindBin(massMin[iSliceVar] * 1.0001), hMassSgn[iSliceVar]->FindBin(massMax[iSliceVar] * 0.999)); @@ -406,10 +523,28 @@ void runMassFitter(const std::string& configFileName) massFitter->setTemplateReflections(hMassRefl[iSliceVar]); } + auto setDscbParameter = [&](const std::vector& vec, void (HFInvMassFitter::*setter)(double)) { + if (static_cast(vec.size()) == nHistograms) { + (massFitter->*setter)(vec[iSliceVar]); + } + }; + setDscbParameter(dscbAlphaLInitial, &HFInvMassFitter::setDscbAlphaLInitialValue); + setDscbParameter(dscbAlphaLLower, &HFInvMassFitter::setDscbAlphaLLowLimit); + setDscbParameter(dscbAlphaLUpper, &HFInvMassFitter::setDscbAlphaLUpLimit); + setDscbParameter(dscbAlphaRInitial, &HFInvMassFitter::setDscbAlphaRInitialValue); + setDscbParameter(dscbAlphaRLower, &HFInvMassFitter::setDscbAlphaRLowLimit); + setDscbParameter(dscbAlphaRUpper, &HFInvMassFitter::setDscbAlphaRUpLimit); + setDscbParameter(dscbNLInitial, &HFInvMassFitter::setDscbNLInitialValue); + setDscbParameter(dscbNLLower, &HFInvMassFitter::setDscbNLLowLimit); + setDscbParameter(dscbNLUpper, &HFInvMassFitter::setDscbNLUpLimit); + setDscbParameter(dscbNRInitial, &HFInvMassFitter::setDscbNRInitialValue); + setDscbParameter(dscbNRLower, &HFInvMassFitter::setDscbNRLowLimit); + setDscbParameter(dscbNRUpper, &HFInvMassFitter::setDscbNRUpLimit); + massFitter->doFit(); auto drawOnCanvas = [&](std::vector& canvas, std::function drawer) { - if (nSliceVarBins > 1) { + if (nHistograms > 1) { canvas[iCanvas]->cd(iSliceVar - NCanvasesMax * iCanvas + 1); } else { canvas[iCanvas]->cd(); @@ -442,6 +577,14 @@ void runMassFitter(const std::string& configFileName) const double meanErr = massFitter->getMeanUncertainty(); const double sigma = massFitter->getSigma(); const double sigmaErr = massFitter->getSigmaUncertainty(); + const double dscbAlphaL = massFitter->getDscbAlphaL(); + const double dscbAlphaR = massFitter->getDscbAlphaR(); + const double dscbNL = massFitter->getDscbNL(); + const double dscbNR = massFitter->getDscbNR(); + const double dscbAlphaLErr = massFitter->getDscbAlphaLUncertainty(); + const double dscbAlphaRErr = massFitter->getDscbAlphaRUncertainty(); + const double dscbNErrL = massFitter->getDscbNLUncertainty(); + const double dscbNErrR = massFitter->getDscbNRUncertainty(); hRawYieldsSignal->SetBinContent(iSliceVar + 1, rawYield); hRawYieldsSignal->SetBinError(iSliceVar + 1, rawYieldErr); @@ -462,6 +605,14 @@ void runMassFitter(const std::string& configFileName) hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma); hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr); hReflectionOverSignal->SetBinContent(iSliceVar + 1, reflOverSgn); + hRawYieldsDscbAlphaL->SetBinContent(iSliceVar + 1, dscbAlphaL); + hRawYieldsDscbAlphaL->SetBinError(iSliceVar + 1, dscbAlphaLErr); + hRawYieldsDscbAlphaR->SetBinContent(iSliceVar + 1, dscbAlphaR); + hRawYieldsDscbAlphaR->SetBinError(iSliceVar + 1, dscbAlphaRErr); + hRawYieldsDscbNL->SetBinContent(iSliceVar + 1, dscbNL); + hRawYieldsDscbNL->SetBinError(iSliceVar + 1, dscbNErrL); + hRawYieldsDscbNR->SetBinContent(iSliceVar + 1, dscbNR); + hRawYieldsDscbNR->SetBinError(iSliceVar + 1, dscbNErrR); if (sgnFunc[iSliceVar] != HFInvMassFitter::SingleGaus) { // TODO foresee DSCB and Voigt cases const double secSigma = massFitter->getSecSigma(); @@ -485,6 +636,7 @@ void runMassFitter(const std::string& configFileName) } hFitConfig->SetBinContent(ConfigBkgFunc, iSliceVar + 1, bkgFunc[iSliceVar]); hFitConfig->SetBinContent(ConfigSgnFunc, iSliceVar + 1, sgnFunc[iSliceVar]); + hFitConfig->SetBinContent(ConfigRandomSeed, iSliceVar + 1, randomSeed); } // save output histograms @@ -500,7 +652,7 @@ void runMassFitter(const std::string& configFileName) } } - for (int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) { + for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { hMass[iSliceVar]->Write(); } hRawYieldsSignal->Write(); @@ -517,6 +669,12 @@ void runMassFitter(const std::string& configFileName) if (enableRefl) { hReflectionOverSignal->Write(); } + if (std::find(sgnFunc.begin(), sgnFunc.end(), HFInvMassFitter::DoubleSidedCrystalBall) != sgnFunc.end()) { + hRawYieldsDscbAlphaL->Write(); + hRawYieldsDscbAlphaR->Write(); + hRawYieldsDscbNL->Write(); + hRawYieldsDscbNR->Write(); + } hFitConfig->Write(); outputFile.Close(); @@ -548,10 +706,10 @@ void setHistoStyle(TH1* histo, Color_t color, Size_t markerSize) histo->SetLineColor(color); } -void divideCanvas(TCanvas* canvas, int nSliceVarBins) +void divideCanvas(TCanvas* canvas, int nHistograms) { - int nCols = std::ceil(std::sqrt(nSliceVarBins)); - int nRows = std::ceil(static_cast(nSliceVarBins) / nCols); + int nCols = std::ceil(std::sqrt(nHistograms)); + int nRows = std::ceil(static_cast(nHistograms) / nCols); canvas->Divide(nCols, nRows); } @@ -629,6 +787,48 @@ void readJsonVector(std::vector& vec, const Document& config, const std::stri } } +template +void readJsonVectorFlexible(std::vector& vec, const Document& config, int nHistograms, const std::string& fieldName, bool isRequired) +{ + if constexpr (!(std::is_same_v, int> || std::is_same_v, double>)) { + static_assert(AlwaysFalse, "readJsonVectorFlexible(): unsupported type!"); + } + if (!vec.empty()) { + throw std::runtime_error("readJsonVectorFlexible(): vector is not empty!"); + } + if (!config.HasMember(fieldName.c_str())) { + if (isRequired) { + throw std::runtime_error("readJsonVectorFlexible(): missing required field " + fieldName); + } else { + return; + } + } + if (config[fieldName.c_str()].IsArray()) { + readJsonVector(vec, config, fieldName); + } else { + const T value = readJsonField(config, fieldName); + vec.assign(nHistograms, value); + } +} + +void readJsonVectorFromHisto(std::vector& vec, const Document& config, const std::string& fileNameFieldName, const std::string& histoNameFieldName) +{ + if (!vec.empty()) { + throw std::runtime_error("readJsonVectorFromHisto(): vector is not empty!"); + } + const auto fileName = readJsonField(config, fileNameFieldName, ""); + const auto histoName = readJsonField(config, histoNameFieldName, ""); + if (fileName.empty() || histoName.empty()) { + return; + } + TFile* inputFile = openFileWithNullptrCheck(fileName); + TH1* histo = getObjectWithNullPtrCheck(inputFile, histoName); + for (int iBin = 1; iBin <= histo->GetNbinsX(); iBin++) { + vec.push_back(histo->GetBinContent(iBin)); + } + inputFile->Close(); +} + int main(int argc, const char* argv[]) { if (argc == 1) {