diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 304f7b50d7a..5f80f2cd444 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -40,7 +40,9 @@ #include #include #include +#include #include +#include #include #include #include @@ -50,6 +52,7 @@ #include #include +#include #include #include #include @@ -66,7 +69,7 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, double maxValue, int fitTypeBkg, int fitTypeSgn, - int randomSeed) : mHistoInvMass(nullptr), + int randomSeed) : mHistoInvMass(histoToFit), mFitOption("L,E"), mMinMass(minValue), mMaxMass(maxValue), @@ -147,18 +150,21 @@ HFInvMassFitter::HFInvMassFitter(TH1* histoToFit, mResidualHist(nullptr), mRatioFrame(nullptr), mWorkspace(nullptr), - mIntegralHisto(0), mIntegralBkg(0), mIntegralSgn(0), mHistoTemplateRefl(nullptr), mDrawBgPrefit(false), mHighlightPeakRegion(false), mRandomSeed(randomSeed), - mRandomGen(nullptr) + mRandomGen(nullptr), + mFitStatus(-999), + mCovQual(-999), + mEdm(-999.), + mMinNll(-999.), + mSgnGlobalCorrelCoeff(-999.), + mCovCorrMatrix(nullptr) { // standard constructor - mHistoInvMass = histoToFit; - mHistoInvMass->SetName("mHistoInvMass"); mHistoInvMass->SetDirectory(nullptr); if (mRandomSeed >= 0) { mRandomGen = new TRandom3(); @@ -170,7 +176,6 @@ HFInvMassFitter::~HFInvMassFitter() { /// destructor - delete mHistoInvMass; delete mHistoTemplateRefl; delete mRooMeanSgn; delete mRooSigmaSgn; @@ -192,7 +197,6 @@ HFInvMassFitter::~HFInvMassFitter() void HFInvMassFitter::doFit() { - mIntegralHisto = mHistoInvMass->Integral(mHistoInvMass->FindBin(mMinMass), mHistoInvMass->FindBin(mMaxMass)); mWorkspace = new RooWorkspace("mWorkspace"); fillWorkspace(*mWorkspace); RooRealVar* mass = mWorkspace->var("mass"); @@ -201,13 +205,12 @@ void HFInvMassFitter::doFit() if (mTypeOfBkgPdf == NoBkg) { // MC mass->setRange("signal", mMass - 3. * mSigmaSgn, mMass + 3. * mSigmaSgn); } else { + mass->setRange("SBL", mMinMass, mMass - mNSigmaForSidebands * mSigmaSgn); if (mTypeOfSgnPdf == GausSec) { // Second Peak fit range - 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 - 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); } @@ -217,18 +220,24 @@ void HFInvMassFitter::doFit() 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 + TH1* histoInvMassSB = dynamic_cast(mHistoInvMass->Clone()); + cutRangesFromHisto(histoInvMassSB, {"SBL", "SBR"}); + RooDataHist sbHistogram("sbHistogram", "sb", *mass, Import(*histoInvMassSB)); + RooAbsPdf* bkgPdf = createBackgroundFitFunction(mWorkspace); // Create background pdf RooAbsPdf* sgnPdf = createSignalFitFunction(mWorkspace); // Create signal pdf + const double integralHisto = integrateHistoInvMassOverWorkspaceRanges({"full"}); // fit MC or Data + RooFitResult* fitResult{nullptr}; if (mTypeOfBkgPdf == NoBkg) { // MC - const ParameterRanges rooNSgnParamRanges{0., 1.2 * mIntegralHisto, 0.3 * mIntegralHisto}; + const ParameterRanges rooNSgnParamRanges{0.5 * integralHisto, 1.5 * integralHisto, integralHisto}; 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("full")); + fitResult = mTotalPdf->chi2FitTo(dataHistogram, Range("full"), Save()); } else { - mTotalPdf->fitTo(dataHistogram, Range("full")); + fitResult = mTotalPdf->fitTo(dataHistogram, Range("full"), Save()); } RooAbsReal* signalIntegralMc = mTotalPdf->createIntegral(*mass, NormSet(*mass), Range("signal")); // sig yield from fit mIntegralSgn = signalIntegralMc->getValV(); @@ -239,31 +248,25 @@ void HFInvMassFitter::doFit() mRatioFrame = mass->frame(Title(Form("%s", mHistoInvMass->GetTitle()))); calculateFitToDataRatio(); } else { // data - const ParameterRanges rooNBkgParamRanges{0., 1.2 * mIntegralHisto, 0.3 * mIntegralHisto}; + const double integralSidebands = integrateHistoInvMassOverWorkspaceRanges({"SBL", "SBR"}); + const ParameterRanges rooNBkgParamRanges{0.5 * integralSidebands, 1.5 * integralSidebands, integralSidebands}; 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)); + std::string sbRanges{"SBL,SBR"}; if (mTypeOfSgnPdf == GausSec) { // two peak fit - if (strcmp(mFitOption.c_str(), "Chi2") == 0) { - mBkgPdf->chi2FitTo(dataHistogram, Range("SBL,SBR,SEC"), Save()); - } else { - mBkgPdf->fitTo(dataHistogram, Range("SBL,SBR,SEC"), Save()); - } - } else { // single peak fit - if (strcmp(mFitOption.c_str(), "Chi2") == 0) { - mBkgPdf->chi2FitTo(dataHistogram, Range("SBL,SBR"), Save()); - } else { - mBkgPdf->fitTo(dataHistogram, Range("SBL,SBR"), Save()); - } + sbRanges.append(",SEC"); } + mBkgPdf->chi2FitTo(sbHistogram, DataError(RooAbsData::SumW2), Save()); + // define the frame to evaluate background sidebands chi2 (bg pdf needs to be plotted within sideband ranges) RooPlot* frameTemporary = mass->frame(Title(Form("%s_temp", mHistoInvMass->GetTitle()))); dataHistogram.plotOn(frameTemporary, Name("data_for_bkgchi2")); - mBkgPdf->plotOn(frameTemporary, Range("SBL", true), Name("Bkg_sidebands")); + mBkgPdf->plotOn(frameTemporary, Range(sbRanges.c_str()), Name("Bkg_sidebands")); mChiSquareOverNdfBkg = frameTemporary->chiSquare("Bkg_sidebands", "data_for_bkgchi2"); // calculate reduced chi2 / NDF of background sidebands (pre-fit) delete frameTemporary; if (mDrawBgPrefit) { - RooAbsPdf* bkgPdfPrefit = dynamic_cast(mBkgPdf->Clone()); - bkgPdfPrefit->plotOn(mInvMassFrame, Range("full"), Name("Bkg_c_prefit"), LineColor(kGray)); + auto* bkgPdfPrefit = dynamic_cast(mBkgPdf->Clone()); + bkgPdfPrefit->plotOn(mInvMassFrame, Range("full"), Normalization(mRooNBkg->getVal(), RooAbsReal::NumEvent), Name("Bkg_c_prefit"), LineColor(kGray)); delete bkgPdfPrefit; } @@ -274,16 +277,17 @@ void HFInvMassFitter::doFit() checkForSignal(estimatedSignal); // SIG's absolute integral in "bkg" range calculateBackground(mBkgYield, mBkgYieldErr); // BG's absolute integral in "bkg" range - 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 + const ParameterRanges rooNSgnParamRanges{0.1 * estimatedSignal, 10 * estimatedSignal, estimatedSignal}; + mRooNSgn = new RooRealVar("mRooNSig", "number of signal", randomizeInitialParameter(rooNSgnParamRanges), rooNSgnParamRanges.lower, rooNSgnParamRanges.upper); // estimated signal yield if (mFixedRawYield > 0) { mRooNSgn->setVal(mFixedRawYield); // fixed signal yield mRooNSgn->setConstant(true); } mSgnPdf = new RooAddPdf("mSgnPdf", "signal fit function", RooArgList(*sgnPdf), RooArgList(*mRooNSgn)); // create reflection template and fit to reflection + RooAbsPdf* reflPdf{nullptr}; if (mHistoTemplateRefl != nullptr) { - RooAbsPdf* reflPdf = createReflectionFitFunction(mWorkspace); // create reflection pdf + reflPdf = createReflectionFitFunction(mWorkspace); // create reflection pdf RooDataHist reflHistogram("reflHistogram", "refl for fit", *mass, Import(*mHistoTemplateRefl)); mReflFrame = mass->frame(); mReflOnlyFrame = mass->frame(Title(Form("%s", mHistoTemplateRefl->GetTitle()))); @@ -302,42 +306,33 @@ void HFInvMassFitter::doFit() mRooNRefl->setConstant(true); setReflFuncFixed(); // fix reflection pdf parameter mTotalPdf = new RooAddPdf("mTotalPdf", "background + signal + reflection fit function", RooArgList(*bkgPdf, *sgnPdf, *reflPdf), RooArgList(*mRooNBkg, *mRooNSgn, *mRooNRefl)); - if (strcmp(mFitOption.c_str(), "Chi2") == 0) { - mTotalPdf->chi2FitTo(dataHistogram); - } else { - mTotalPdf->fitTo(dataHistogram); - } - mTotalPdf->plotOn(mInvMassFrame, Name("Tot_c")); + } else { + mTotalPdf = new RooAddPdf("mTotalPdf", "background + signal pdf", RooArgList(*bkgPdf, *sgnPdf), RooArgList(*mRooNBkg, *mRooNSgn)); + } + + if (strcmp(mFitOption.c_str(), "Chi2") == 0) { + fitResult = mTotalPdf->chi2FitTo(dataHistogram, Save()); + } else { + fitResult = mTotalPdf->fitTo(dataHistogram, Save()); + } + + plotBkg(mTotalPdf); + mTotalPdf->plotOn(mInvMassFrame, Name("Tot_c"), LineColor(kBlue)); + if (mHistoTemplateRefl != nullptr) { mReflPdf = new RooAddPdf("mReflPdf", "reflection fit function", RooArgList(*reflPdf), RooArgList(*mRooNRefl)); RooAddPdf const reflBkgPdf("reflBkgPdf", "reflBkgPdf", RooArgList(*bkgPdf, *reflPdf), RooArgList(*mRooNBkg, *mRooNRefl)); reflBkgPdf.plotOn(mInvMassFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineStyle(7), LineColor(kRed + 1), Name("ReflBkg_c")); - plotBkg(mTotalPdf); // plot bkg pdf in total pdf - plotRefl(mTotalPdf); // plot reflection in total pdf - mChiSquareOverNdfTotal = mInvMassFrame->chiSquare("Tot_c", "data_c"); // calculate reduced chi2 / NDF - - // plot residual distribution - mResidualHist = mInvMassFrame->residHist("data_c", "ReflBkg_c"); - mResidualFrame = mass->frame(Title("Residual Distribution")); - mResidualFrame->addPlotable(mResidualHist, "p"); - mSgnPdf->plotOn(mResidualFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineColor(kBlue)); + plotRefl(mTotalPdf); // plot reflection in total pdf } else { - mTotalPdf = new RooAddPdf("mTotalPdf", "background + signal pdf", RooArgList(*bkgPdf, *sgnPdf), RooArgList(*mRooNBkg, *mRooNSgn)); - if (strcmp(mFitOption.c_str(), "Chi2") == 0) { - mTotalPdf->chi2FitTo(dataHistogram); - } else { - mTotalPdf->fitTo(dataHistogram); - } - plotBkg(mTotalPdf); - mTotalPdf->plotOn(mInvMassFrame, Name("Tot_c"), LineColor(kBlue)); mSgnPdf->plotOn(mInvMassFrame, Normalization(1.0, RooAbsReal::RelativeExpected), DrawOption("F"), FillColor(TColor::GetColorTransparent(kBlue, 0.2)), VLines()); - mChiSquareOverNdfTotal = mInvMassFrame->chiSquare("Tot_c", "data_c"); // calculate reduced chi2 / DNF - - // plot residual distribution - mResidualFrame = mass->frame(Title("Residual Distribution")); - mResidualHist = mInvMassFrame->residHist("data_c", "Bkg_c"); - mResidualFrame->addPlotable(mResidualHist, "P"); - mSgnPdf->plotOn(mResidualFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineColor(kBlue)); } + mChiSquareOverNdfTotal = mInvMassFrame->chiSquare("Tot_c", "data_c"); // calculate reduced chi2 / NDF + // plot residual distribution + mResidualFrame = mass->frame(Title(Form("%s", mHistoInvMass->GetTitle()))); + mResidualHist = mInvMassFrame->residHist("data_c", mHistoTemplateRefl ? "ReflBkg_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()); bkgIntegral = mBkgPdf->createIntegral(*mass, NormSet(*mass), Range("bkgForSignificance")); mIntegralBkg = bkgIntegral->getValV(); @@ -352,6 +347,12 @@ void HFInvMassFitter::doFit() mRatioFrame = mass->frame(Title(Form("%s", mHistoInvMass->GetTitle()))); calculateFitToDataRatio(); } + mFitStatus = fitResult->status(); + mCovQual = fitResult->covQual(); + mEdm = fitResult->edm(); + mMinNll = fitResult->minNll(); + mSgnGlobalCorrelCoeff = fitResult->globalCorr("mRooNSig"); + mCovCorrMatrix = fillCovCorrMatrix(fitResult); } void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const @@ -624,6 +625,7 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& mInvMassFrame->GetXaxis()->SetTitleOffset(1.2); mInvMassFrame->GetYaxis()->SetTitleOffset(1.8); gPad->SetLeftMargin(0.15); + gPad->SetRightMargin(0.02); mInvMassFrame->GetYaxis()->SetTitle(Form("%s", mHistoInvMass->GetYaxis()->GetTitle())); mInvMassFrame->GetXaxis()->SetTitle(Form("%s", mHistoInvMass->GetXaxis()->GetTitle())); mInvMassFrame->Draw(); @@ -636,8 +638,12 @@ void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& // draw residual distribution on canvas void HFInvMassFitter::drawResidual(TVirtualPad* pad) { + if (mResidualFrame == nullptr) { + printf("Warning HFInvMassFitter::drawResidual(): mResidualFrame == nullptr and will not be drawn\n"); + return; + } pad->cd(); - mResidualFrame->GetYaxis()->SetTitle(""); + mResidualFrame->GetYaxis()->SetTitle(mHistoInvMass->GetYaxis()->GetTitle()); auto* textInfo = new TPaveText(0.12, 0.65, 0.47, .89, "NDC"); textInfo->SetBorderSize(0); textInfo->SetFillStyle(0); @@ -650,6 +656,9 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad) textInfo->AddText(Form("#sigma_{2} = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); } mResidualFrame->addObject(textInfo); + gPad->SetLeftMargin(0.15); + gPad->SetRightMargin(0.02); + mResidualFrame->GetYaxis()->SetTitleOffset(1.8); mResidualFrame->Draw(); highlightPeakRegion(mResidualFrame); } @@ -657,6 +666,10 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad) // draw ratio on canvas void HFInvMassFitter::drawRatio(TVirtualPad* pad) { + if (mRatioFrame == nullptr) { + printf("Warning HFInvMassFitter::drawRatio(): mRatioFrame == nullptr and will not be drawn\n"); + return; + } pad->cd(); mRatioFrame->GetXaxis()->SetTitleOffset(1.2); mRatioFrame->GetYaxis()->SetTitleOffset(1.5); @@ -668,6 +681,8 @@ void HFInvMassFitter::drawRatio(TVirtualPad* pad) line->SetLineStyle(2); line->SetLineWidth(2); mRatioFrame->addObject(line); + gPad->SetLeftMargin(0.15); + gPad->SetRightMargin(0.02); mRatioFrame->Draw(); highlightPeakRegion(mRatioFrame); } @@ -703,7 +718,7 @@ void HFInvMassFitter::drawReflection(TVirtualPad* pad) } // calculate signal yield via bin counting -void HFInvMassFitter::countSignal(double& signal, double& signalErr) const +void HFInvMassFitter::countSignal(double& signal, double& errSignal) const { const auto& histoForCounting = mTypeOfBkgPdf == NoBkg ? mInvMassFrame->getHist("data_c") : mResidualHist; @@ -730,7 +745,7 @@ void HFInvMassFitter::countSignal(double& signal, double& signalErr) const sumErrorsSquare += square(histoForCounting->GetErrorY(binForMaxSgn - 1) * binForMaxSgnFraction); signal = sumValues; - signalErr = std::sqrt(sumErrorsSquare); + errSignal = std::sqrt(sumErrorsSquare); } // calculate signal yield @@ -748,8 +763,9 @@ void HFInvMassFitter::calculateBackground(double& bkg, double& errBkg) const errBkg = 0.; return; } - bkg = mRooNBkg->getVal() * mIntegralBkg; - errBkg = mRooNBkg->getError() * mIntegralBkg; + const double bgCoefficient = mTypeOfSgnPdf == DoubleSidedCrystalBall ? 1. : mIntegralBkg; + bkg = mRooNBkg->getVal() * bgCoefficient; + errBkg = mRooNBkg->getError() * bgCoefficient; } // calculate significance @@ -769,17 +785,9 @@ void HFInvMassFitter::calculateSignificance(double& significance, double& errSig // estimate Signal void HFInvMassFitter::checkForSignal(double& estimatedSignal) { - auto const [minForSgn, maxForSgn] = getRangesOfSignal(); - int const binForMinSgn = mHistoInvMass->FindBin(minForSgn); - int const binForMaxSgn = mHistoInvMass->FindBin(maxForSgn); - - double sum = 0; - for (int i = binForMinSgn; i <= binForMaxSgn; i++) { - sum += mHistoInvMass->GetBinContent(i); - } - double bkg{}, errBkg{}; - calculateBackground(bkg, errBkg); - estimatedSignal = sum - bkg; + const double integralHisto = integrateHistoInvMassOverWorkspaceRanges({"full"}); + const double bkg = mRooNBkg->getVal(); + estimatedSignal = integralHisto - bkg; } // Estimate ranges where signal is located to be used in countSignal() and checkForSignal() @@ -788,11 +796,10 @@ 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); } + const double mean = mRooMeanSgn->getVal(); + const double sigma = mRooSecSigmaSgn->getVal(); + return std::make_pair(mean - mNSigmaForSgn * sigma, mean + mNSigmaForSgn * sigma); } // Create Background Fit Function @@ -904,7 +911,7 @@ void HFInvMassFitter::calculateFitToDataRatio() const return; } - RooHist* ratioHist = new RooHist(); + auto* ratioHist = new RooHist(); for (int i = 0; i < dataHist->GetN(); ++i) { double x{}, dataY{}, dataErr{}; @@ -1135,17 +1142,86 @@ double HFInvMassFitter::randomizeInitialParameter(const ParameterRanges& paramet } const auto sigma = parameterRanges.sigma < 0 ? (parameterRanges.upper - parameterRanges.lower) / DefaultSigmaFraction : parameterRanges.sigma; - double result; + 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); + std::array errorMessage{}; + std::snprintf(errorMessage.data(), errorMessage.size(), "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.data()); } } while (result < parameterRanges.lower || result > parameterRanges.upper); return result; } + +double HFInvMassFitter::integrateHistoInvMassOverWorkspaceRanges(const std::vector& ranges) const +{ + double sumEntries{0.}; + double sumLengths{0.}; + for (const auto& range : ranges) { + const auto [lo, hi] = mWorkspace->var("mass")->getRange(range.c_str()); + const auto binLo = mHistoInvMass->FindBin(lo); + const auto binHi = mHistoInvMass->FindBin(hi); + sumEntries += mHistoInvMass->Integral(binLo + 1, binHi - 1); + sumEntries += mHistoInvMass->GetBinContent(binLo) * (mHistoInvMass->GetBinLowEdge(binLo + 1) - lo) / mHistoInvMass->GetBinWidth(binLo); + sumEntries += mHistoInvMass->GetBinContent(binHi) * (hi - mHistoInvMass->GetBinLowEdge(binHi)) / mHistoInvMass->GetBinWidth(binHi); + sumLengths += (hi - lo); + } + const auto [fullLo, fullHi] = mWorkspace->var("mass")->getRange("full"); + const double fullLength = fullHi - fullLo; + + return sumEntries / sumLengths * fullLength; +} + +void HFInvMassFitter::cutRangesFromHisto(TH1* histo, const std::vector& ranges) const +{ + auto* mass = mWorkspace->var("mass"); + + for (int iBin = 1, nBins = histo->GetNbinsX(); iBin <= nBins; ++iBin) { + const double binLow = histo->GetBinLowEdge(iBin); + const double binHigh = binLow + histo->GetBinWidth(iBin); + + bool overlapsAnyRange = false; + for (const auto& range : ranges) { + const double rangeMin = mass->getMin(range.c_str()); + const double rangeMax = mass->getMax(range.c_str()); + if (std::max(binLow, rangeMin) <= std::min(binHigh, rangeMax)) { + overlapsAnyRange = true; + break; + } + } + + if (!overlapsAnyRange) { + histo->SetBinContent(iBin, 0.); + histo->SetBinError(iBin, 1.e9); + } + } +} + +TH2* HFInvMassFitter::fillCovCorrMatrix(const RooFitResult* fitResult) +{ + const RooArgList& pars = fitResult->floatParsFinal(); + const int nPars = static_cast(pars.size()); + + TH2* hMatrix = new TH2D("covCorrMatrix", "covariance (upper left + diagonal) and correlation (lower right) matrix", nPars, 0, nPars, nPars, 0, nPars); + for (int iPar = 0; iPar < nPars; ++iPar) { + hMatrix->GetXaxis()->SetBinLabel(iPar + 1, pars[iPar].GetName()); + hMatrix->GetYaxis()->SetBinLabel(iPar + 1, pars[iPar].GetName()); + } + + const TMatrixDSym& covMatrix = fitResult->covarianceMatrix(); + const TMatrixDSym& corrMatrix = fitResult->correlationMatrix(); + + for (int iPar = 0; iPar < nPars; ++iPar) { + hMatrix->SetBinContent(iPar + 1, iPar + 1, covMatrix(iPar, iPar)); + for (int jPar = iPar + 1; jPar < nPars; ++jPar) { + hMatrix->SetBinContent(iPar + 1, jPar + 1, covMatrix(iPar, jPar)); + hMatrix->SetBinContent(jPar + 1, iPar + 1, corrMatrix(iPar, jPar)); + } + } + + return hMatrix; +} diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 6841cf15a55..be5c95be829 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -22,11 +22,15 @@ #ifndef PWGHF_D2H_MACROS_HFINVMASSFITTER_H_ #define PWGHF_D2H_MACROS_HFINVMASSFITTER_H_ +#include +#include +#include #include #include #include #include #include +#include #include #include #include @@ -52,28 +56,28 @@ class HFInvMassFitter : public TNamed public: enum TypeOfBkgPdf { Expo = 0, - Poly1 = 1, - Poly2 = 2, - Pow = 3, - PowExpo = 4, - Poly3 = 5, - NoBkg = 6, + Poly1, // 1 + Poly2, // 2 + Pow, // 3 + PowExpo, // 4 + Poly3, // 5 + NoBkg, // 6 NTypesOfBkgPdf }; std::array namesOfBkgPdf{"bkgFuncExpo", "bkgFuncPoly1", "bkgFuncPoly2", "bkgFuncPow", "bkgFuncPowExpo", "bkgFuncPoly3"}; enum TypeOfSgnPdf { SingleGaus = 0, - DoubleGaus = 1, - DoubleGausSigmaRatioPar = 2, - GausSec = 3, - DoubleSidedCrystalBall = 4, + DoubleGaus, // 1 + DoubleGausSigmaRatioPar, // 2 + GausSec, // 3 + DoubleSidedCrystalBall, // 4 NTypesOfSgnPdf }; enum TypeOfReflPdf { SingleGausRefl = 0, - DoubleGausRefl = 1, - Poly3Refl = 2, - Poly6Refl = 3, + DoubleGausRefl, // 1 + Poly3Refl, // 2 + Poly6Refl, // 3 NTypesOfReflPdf }; std::array namesOfReflPdf{"reflFuncGaus", "reflFuncDoubleGaus", "reflFuncPoly3", "reflFuncPoly6"}; @@ -84,9 +88,9 @@ class HFInvMassFitter : public TNamed void setUseLikelihoodFit() { mFitOption = "L,E"; } void setUseChi2Fit() { mFitOption = "Chi2"; } void setFitOption(const std::string& opt) { mFitOption = opt; } - RooAbsPdf* createBackgroundFitFunction(RooWorkspace* w1) const; - RooAbsPdf* createSignalFitFunction(RooWorkspace* w1); - RooAbsPdf* createReflectionFitFunction(RooWorkspace* w1) const; + RooAbsPdf* createBackgroundFitFunction(RooWorkspace* workspace) const; + RooAbsPdf* createSignalFitFunction(RooWorkspace* workspace); + RooAbsPdf* createReflectionFitFunction(RooWorkspace* workspace) const; void setFitRange(double minValue, double maxValue); void setFitFunctions(int fitTypeBkg, int fitTypeSgn); @@ -125,8 +129,8 @@ class HFInvMassFitter : public TNamed 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 plotBkg(RooAbsPdf* pdf, Color_t color = kRed); + void plotRefl(RooAbsPdf* pdf); void setReflFuncFixed(); void doFit(); void setInitialReflOverSgn(double reflOverSgn) { mReflOverSgn = reflOverSgn; } @@ -161,16 +165,22 @@ class HFInvMassFitter : public TNamed [[nodiscard]] double getFracDoubleGaus() const { return mRooFracDoubleGaus->getVal(); } [[nodiscard]] double getFracDoubleGausUncertainty() const { return mRooFracDoubleGaus->getError(); } [[nodiscard]] double getReflOverSig() const { return mReflPdf != nullptr ? mReflOverSgn : 0.; } - void calculateSignal(double& signal, double& signalErr) const; - void countSignal(double& signal, double& signalErr) const; - void calculateBackground(double& bkg, double& bkgErr) const; - void calculateSignificance(double& significance, double& significanceErr) const; + [[nodiscard]] int getFitStatus() const { return mFitStatus; } + [[nodiscard]] int getCovQual() const { return mCovQual; } + [[nodiscard]] double getEDM() const { return mEdm; } + [[nodiscard]] double getMinNll() const { return mMinNll; } + [[nodiscard]] double getSgnGlobalCorrelCoeff() const { return mSgnGlobalCorrelCoeff; } + [[nodiscard]] TH2* getCovCorrMatrix() const { return mCovCorrMatrix; } + void calculateSignal(double& signal, double& errSignal) const; + void countSignal(double& signal, double& errSignal) const; + void calculateBackground(double& bkg, double& errBkg) const; + void calculateSignificance(double& significance, double& errSignificance) const; void checkForSignal(double& estimatedSignal); void calculateFitToDataRatio() const; - void drawFit(TVirtualPad* c, const std::vector& plotLabels, bool writeParInfo = true); - void drawResidual(TVirtualPad* c); - void drawRatio(TVirtualPad* c); - void drawReflection(TVirtualPad* c); + void drawFit(TVirtualPad* pad, const std::vector& plotLabels, bool writeParInfo = true); + void drawResidual(TVirtualPad* pad); + void drawRatio(TVirtualPad* pad); + void drawReflection(TVirtualPad* pad); private: HFInvMassFitter(const HFInvMassFitter& source); @@ -179,6 +189,9 @@ class HFInvMassFitter : public TNamed 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; + [[nodiscard]] double integrateHistoInvMassOverWorkspaceRanges(const std::vector& ranges) const; + void cutRangesFromHisto(TH1* histo, const std::vector& ranges) const; + static TH2* fillCovCorrMatrix(const RooFitResult* fitResult); TH1* mHistoInvMass; // histogram to fit std::string mFitOption; @@ -261,7 +274,6 @@ class HFInvMassFitter : public TNamed RooHist* mResidualHist; /// residual histogram RooPlot* mRatioFrame; /// fit/data ratio frame 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 @@ -269,6 +281,12 @@ class HFInvMassFitter : public TNamed 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 + int mFitStatus; /// fit result status, see https://root-forum.cern.ch/t/meaning-of-values-returned-by-roofitresult-status/16355/2 + int mCovQual; /// fit result covariance matrix quality, see https://root.cern.ch/doc/v620/Minuit2Minimizer_8cxx_source.html#l01121 + double mEdm; /// fit quality metrics: Estimated Distance to Minimum + double mMinNll; /// fit quality metrics: minimum negative log-likelihood (NLL) value achieved at the best-fit parameter values + double mSgnGlobalCorrelCoeff; /// global correlation coefficient of mRooNSgn with other fit parameters + TH2* mCovCorrMatrix; /// covariance (upper left + diagonal) and correlation (lower right) matrix of free fit parameters ClassDefOverride(HFInvMassFitter, 1); }; diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index b208e2dbf46..6ed3ccee7c1 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -41,6 +41,7 @@ #include #include #include +#include #include #include #include @@ -52,8 +53,6 @@ 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"); @@ -88,8 +87,8 @@ void runMassFitter(const std::string& configFileName) } Document config; - char readBuffer[65536]; - FileReadStream is(configFile, readBuffer, sizeof(readBuffer)); + std::array readBuffer{}; + FileReadStream is(configFile, readBuffer.data(), readBuffer.size()); config.ParseStream(is); fclose(configFile); @@ -97,12 +96,12 @@ void runMassFitter(const std::string& configFileName) 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"); - std::string const collisionSystem = readJsonField(config, "CollisionSystem", ""); - std::string const inputFileName = readJsonField(config, "InFileName"); - std::string const reflFileName = readJsonField(config, "ReflFileName", ""); + auto const isMc = readJsonField(config, "IsMC"); + auto const writeSignalPar = readJsonField(config, "WriteSignalPar", true); + auto const particleName = readJsonField(config, "Particle"); + auto const collisionSystem = readJsonField(config, "CollisionSystem", ""); + auto const inputFileName = readJsonField(config, "InFileName"); + auto const reflFileName = readJsonField(config, "ReflFileName", ""); TString outputFileName = readJsonField(config, "OutFileName", "mInvFit.root"); std::vector inputHistoName; @@ -149,23 +148,23 @@ void runMassFitter(const std::string& configFileName) readJsonVector(fdSecPeakHistoName, config, "FDSecPeakHistoName"); readJsonVector(signalSecPeakHistoName, config, "SignalSecPeakHistoName"); - const bool fixMean = readJsonField(config, "FixMean", false); - const std::string meanFile = readJsonField(config, "MeanFile", ""); + const auto fixMean = readJsonField(config, "FixMean", false); + const auto meanFile = readJsonField(config, "MeanFile", ""); readJsonVectorFlexible(fixMeanManual, config, nHistograms, "FixMeanManual"); - const bool fixSigma = readJsonField(config, "FixSigma", false); - const std::string sigmaFile = readJsonField(config, "SigmaFile", ""); + const auto fixSigma = readJsonField(config, "FixSigma", false); + const auto sigmaFile = readJsonField(config, "SigmaFile", ""); readJsonVectorFlexible(fixSigmaManual, config, nHistograms, "FixSigmaManual"); - const bool fixSecondSigma = readJsonField(config, "FixSecondSigma", false); - const std::string secondSigmaFile = readJsonField(config, "SecondSigmaFile", ""); + const auto fixSecondSigma = readJsonField(config, "FixSecondSigma", false); + const auto secondSigmaFile = readJsonField(config, "SecondSigmaFile", ""); readJsonVectorFlexible(fixSecondSigmaManual, config, nHistograms, "FixSecondSigmaManual"); - const bool fixFracDoubleGaus = readJsonField(config, "FixFracDoubleGaus", false); - const std::string fracDoubleGausFile = readJsonField(config, "FracDoubleGausFile", ""); + const auto fixFracDoubleGaus = readJsonField(config, "FixFracDoubleGaus", false); + const auto fracDoubleGausFile = readJsonField(config, "FracDoubleGausFile", ""); readJsonVectorFlexible(fixFracDoubleGausManual, config, nHistograms, "FixFracDoubleGausManual"); - const bool fixDscbTailParams = readJsonField(config, "FixDscbTailParams", false); + const auto fixDscbTailParams = readJsonField(config, "FixDscbTailParams", false); const TString sliceVarName = readJsonField(config, "SliceVarName"); const TString sliceVarUnit = readJsonField(config, "SliceVarUnit"); @@ -178,18 +177,18 @@ void runMassFitter(const std::string& configFileName) readJsonVectorFlexible(nRebin, config, nHistograms, "Rebin", true); - bool const includeSecPeak = readJsonField(config, "InclSecPeak", false); - bool const useLikelihood = readJsonField(config, "UseLikelihood"); + auto const includeSecPeak = readJsonField(config, "InclSecPeak", false); + auto const useLikelihood = readJsonField(config, "UseLikelihood"); 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 randomSeed = readJsonField(config, "RandomSeed", -1); - const double nSigmaForSideband = readJsonField(config, "NSigmaForSideband", 3.); - const double nSigmaForSignal = readJsonField(config, "NSigmaForSignal", 3.); + const auto enableRefl = readJsonField(config, "EnableRefl", false); + const auto drawBgPrefit = readJsonField(config, "DrawBgPrefit", true); + const auto highlightPeakRegion = readJsonField(config, "HighlightPeakRegion", true); + const auto randomSeed = readJsonField(config, "RandomSeed", -1); + const auto nSigmaForSideband = readJsonField(config, "NSigmaForSideband", 3.); + const auto nSigmaForSignal = readJsonField(config, "NSigmaForSignal", 3.); readJsonVector(dscbAlphaLInitial, config, "DscbAlphaLInitial"); readJsonVector(dscbAlphaLLower, config, "DscbAlphaLLower"); @@ -218,7 +217,7 @@ void runMassFitter(const std::string& configFileName) std::vector sliceVarLimits(nHistograms + 1); - auto checkVectorSize = [&](const auto& vec, const std::string& name = "", const bool isEmptyOk = false) { + auto checkVectorSize = [&](const auto& vec, const std::string& name, const bool isEmptyOk = false) { if (vec.size() != static_cast(nHistograms)) { if (isEmptyOk && vec.empty()) { return; @@ -305,13 +304,13 @@ void runMassFitter(const std::string& configFileName) {"LcToPK0s", {"pK^{0}_{s}", "Lambda_c+", "#Lambda_{c}^{+}", "pK^{0}_{s}"}}, {"Dstar", {"D^{0}pi^{+}", "D*+", "D^{*+}", "D^{0}#pi^{+}"}}, {"XicToXiPiPi", {"#Xi#pi#pi", "Xi_c+", "#Xi_{c}^{+}", "#Xi^{-}#pi^{+}#pi^{+}"}}}; - if (particles.find(particleName.c_str()) == particles.end()) { + if (particles.find(particleName) == particles.end()) { throw std::runtime_error("ERROR: only Dplus, D0, Ds, LcToPKPi, LcToPK0s, Dstar and XicToXiPiPi particles supported! Exit"); } const auto& particle = particles[particleName.c_str()]; const std::string massAxisTitle = "#it{M}(" + particle.decayProducts + ") (GeV/#it{c}^{2})"; const double massPDG = TDatabasePDG::Instance()->GetParticle(particle.pdgName.c_str())->Mass(); - const std::vector plotLabels = {(particle.decayFormulaLhs + " #rightarrow " + particle.decayFormulaRhs + " + c.c.").c_str(), collisionSystem}; + const std::vector plotLabels = {particle.decayFormulaLhs + " #rightarrow " + particle.decayFormulaRhs + " + c.c.", collisionSystem}; // load inv-mass histograms auto* inputFile = openFileWithNullptrCheck(inputFileName); @@ -321,6 +320,7 @@ void runMassFitter(const std::string& configFileName) std::vector hMassSgn(nHistograms); std::vector hMassRefl(nHistograms); std::vector hMass(nHistograms); + std::vector hCovCorr(nHistograms); for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { if (!isMc) { @@ -358,22 +358,22 @@ void runMassFitter(const std::string& configFileName) } // define output histos - 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()); + 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, @@ -385,15 +385,31 @@ void runMassFitter(const std::string& configFileName) ConfigRandomSeed, NConfigsToSave }; - 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); + enum { + FitResultStatus = 1, + FitResultCovQual, + FitResultEdm, + FitResultMinNll, + FitResultNSgnGCC, + NFitResultsToSave + }; + auto* hFitConfig = new TH2F("hFitConfig", "Fit Configurations", NConfigsToSave - 1, 0, NConfigsToSave - 1, nHistograms, sliceVarLimits.data()); + constexpr std::array HFitConfigXLabel = {"mass min", "mass max", "rebin num", "fix sigma", "bkg func", "sgn func", "rnd seed"}; + auto* hFitResult = new TH2F("hFitResult", "Fit Result", NFitResultsToSave - 1, 0, NFitResultsToSave - 1, nHistograms, sliceVarLimits.data()); + constexpr std::array HFitResultXLabel = {"status", "cov qual", "edm", "minNLL", "N Sig GCC"}; for (int i = 0; i < NConfigsToSave - 1; i++) { - hFitConfig->GetXaxis()->SetBinLabel(i + 1, hFitConfigXLabel[i]); + hFitConfig->GetXaxis()->SetBinLabel(i + 1, HFitConfigXLabel[i]); + } + for (int i = 0; i < NFitResultsToSave - 1; i++) { + hFitResult->GetXaxis()->SetBinLabel(i + 1, HFitResultXLabel[i]); + } + for (const auto& h : {hFitConfig, hFitResult}) { + h->SetStats(false); + h->LabelsDeflate("X"); + h->LabelsDeflate("Y"); + h->LabelsOption("v"); + h->GetYaxis()->SetTitle(sliceVarName + " (" + sliceVarUnit + ")"); } - hFitConfig->LabelsDeflate("X"); - hFitConfig->LabelsDeflate("Y"); - hFitConfig->LabelsOption("v"); setHistoStyle(hRawYieldsSignal); setHistoStyle(hRawYieldsSignalCounted); @@ -434,32 +450,34 @@ void runMassFitter(const std::string& configFileName) TH1* hSecondSigmaToFix = getHistToFix(fixSecondSigma, fixSecondSigmaManual, secondSigmaFile, "SecSigma"); TH1* hFracDoubleGausToFix = getHistToFix(fixFracDoubleGaus, fixFracDoubleGausManual, fracDoubleGausFile, "FracDoubleGaus"); - int canvasSize[2] = {1920, 1080}; + std::array canvasSize = {1920, 1080}; 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(nHistograms) / NCanvasesMax); + const int nCanvases = (nHistograms + NCanvasesMax - 1) / NCanvasesMax; + const int nPads = (nCanvases == 1) ? 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) ? nHistograms : NCanvasesMax; - canvasMass[iCanvas] = new TCanvas(Form("canvasMass%d", iCanvas), Form("canvasMass%d", iCanvas), canvasSize[0], canvasSize[1]); - divideCanvas(canvasMass[iCanvas], nPads); - canvasResiduals[iCanvas] = new TCanvas(Form("canvasResiduals%d", iCanvas), Form("canvasResiduals%d", iCanvas), canvasSize[0], canvasSize[1]); - divideCanvas(canvasResiduals[iCanvas], nPads); - - canvasRatio[iCanvas] = new TCanvas(Form("canvasRatio%d", iCanvas), Form("canvasRatio%d", iCanvas), canvasSize[0], canvasSize[1]); - divideCanvas(canvasRatio[iCanvas], nPads); - - if (enableRefl) { - canvasRefl[iCanvas] = new TCanvas(Form("canvasRefl%d", iCanvas), Form("canvasRefl%d", iCanvas), canvasSize[0], canvasSize[1]); - divideCanvas(canvasRefl[iCanvas], nPads); + std::vector*> canvasTypes{&canvasMass, &canvasResiduals, &canvasRatio}; + std::vector canvasTypeNames{"canvasMass", "canvasResiduals", "canvasRatio"}; + if (enableRefl) { + canvasTypes.push_back(&canvasRefl); + canvasTypeNames.push_back("canvasRefl"); + } + for (int iCanvasType = 0, nCanvasTypeNames = static_cast(canvasTypes.size()); iCanvasType < nCanvasTypeNames; ++iCanvasType) { + const auto canvasTypeName = canvasTypeNames[iCanvasType]; + for (int iCanvas = 0; iCanvas < nCanvases; iCanvas++) { + auto& canvas = (*canvasTypes[iCanvasType])[iCanvas]; + const auto canvasName = Form("%s%d", canvasTypeName, iCanvas); + canvas = new TCanvas(canvasName, canvasName, canvasSize[0], canvasSize[1]); + canvas->SetTicks(1, 1); + divideCanvas(canvas, nPads); } } @@ -472,16 +490,27 @@ void runMassFitter(const std::string& configFileName) hMass[iSliceVar]->SetTitle(Form("%s;%s;Counts per %0.1f MeV/#it{c}^{2}", ptTitle.Data(), massAxisTitle.c_str(), hMass[iSliceVar]->GetBinWidth(1) * 1000)); - hMass[iSliceVar]->SetName(Form("MassForFit%d", iSliceVar)); + hMass[iSliceVar]->SetName(Form("hMassForFit%d", iSliceVar + 1)); if (enableRefl) { hMassRefl[iSliceVar]->Rebin(nRebin[iSliceVar]); hMassSgn[iSliceVar]->Rebin(nRebin[iSliceVar]); } + const auto hMassLo = hMass[iSliceVar]->GetXaxis()->GetXmin(); + if (massMin[iSliceVar] < hMassLo) { + printf("Warning! massMin[%d] is less than hMass[%d] left edge (%f vs %f) and will be assigned the value of the latter\n", iSliceVar, iSliceVar, massMin[iSliceVar], hMassLo); + massMin[iSliceVar] = hMassLo; + } + const auto hMassUp = hMass[iSliceVar]->GetXaxis()->GetXmax(); + if (massMax[iSliceVar] > hMassUp) { + printf("Warning! massMax[%d] is greater than hMass[%d] right edge (%f vs %f) and will be assigned the value of the latter\n", iSliceVar, iSliceVar, massMax[iSliceVar], hMassUp); + massMax[iSliceVar] = hMassUp; + } + double reflOverSgn = 0; - HFInvMassFitter* massFitter = new HFInvMassFitter(hMass[iSliceVar], massMin[iSliceVar], massMax[iSliceVar], bkgFunc[iSliceVar], sgnFunc[iSliceVar], randomSeed); + auto* massFitter = new HFInvMassFitter(hMass[iSliceVar], massMin[iSliceVar], massMax[iSliceVar], bkgFunc[iSliceVar], sgnFunc[iSliceVar], randomSeed); massFitter->setDrawBgPrefit(drawBgPrefit); massFitter->setNumberOfSigmaForSidebands(nSigmaForSideband); massFitter->setNumberOfSigmaForSignal(nSigmaForSignal); @@ -495,27 +524,27 @@ void runMassFitter(const std::string& configFileName) massFitter->setUseChi2Fit(); } - auto setFixedValue = [&iSliceVar](bool const& isFix, std::vector const& fixManual, const TH1* histToFix, std::function setFunc, std::string const& var) -> void { + auto setFixedValue = [&iSliceVar, massFitter](bool const& isFix, std::vector const& fixManual, const TH1* histToFix, void (HFInvMassFitter::*setter)(double), std::string const& var) -> void { if (isFix) { if (fixManual.empty() && histToFix == nullptr) { throw std::runtime_error("Histogram to fix " + var + " is null while isFix==true and fixManual is empty"); } const auto valueToFix = fixManual.empty() ? histToFix->GetBinContent(iSliceVar + 1) : fixManual[iSliceVar]; - setFunc(valueToFix); + (massFitter->*setter)(valueToFix); printf("*****************************\n"); printf("FIXED %s: %f\n", var.data(), valueToFix); printf("*****************************\n"); } }; - setFixedValue(fixMean, fixMeanManual, hMeanToFix, std::bind(&HFInvMassFitter::setFixGaussianMean, massFitter, std::placeholders::_1), "MEAN"); - 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"); + setFixedValue(fixMean, fixMeanManual, hMeanToFix, &HFInvMassFitter::setFixGaussianMean, "MEAN"); + setFixedValue(fixSigma, fixSigmaManual, hSigmaToFix, &HFInvMassFitter::setFixGaussianSigma, "SIGMA"); + setFixedValue(fixSecondSigma, fixSecondSigmaManual, hSecondSigmaToFix, &HFInvMassFitter::setFixSecondGaussianSigma, "SECOND SIGMA"); + setFixedValue(fixFracDoubleGaus, fixFracDoubleGausManual, hFracDoubleGausToFix, &HFInvMassFitter::setFixFrac2Gaus, "FRAC DOUBLE GAUS"); + setFixedValue(fixDscbTailParams, dscbAlphaLInitial, nullptr, &HFInvMassFitter::setFixDscbAlphaL, "DSCB ALPHA LEFT"); + setFixedValue(fixDscbTailParams, dscbAlphaRInitial, nullptr, &HFInvMassFitter::setFixDscbAlphaR, "DSCB ALPHA RIGHT"); + setFixedValue(fixDscbTailParams, dscbNLInitial, nullptr, &HFInvMassFitter::setFixDscbNL, "DSCB N LEFT"); + setFixedValue(fixDscbTailParams, dscbNRInitial, nullptr, &HFInvMassFitter::setFixDscbNR, "DSCB N RIGHT"); if (!isMc && enableRefl) { reflOverSgn = hMassSgn[iSliceVar]->Integral(hMassSgn[iSliceVar]->FindBin(massMin[iSliceVar] * 1.0001), hMassSgn[iSliceVar]->FindBin(massMax[iSliceVar] * 0.999)); @@ -542,9 +571,13 @@ void runMassFitter(const std::string& configFileName) setDscbParameter(dscbNRLower, &HFInvMassFitter::setDscbNRLowLimit); setDscbParameter(dscbNRUpper, &HFInvMassFitter::setDscbNRUpLimit); - massFitter->doFit(); + try { + massFitter->doFit(); + } catch (const std::exception& e) { + printf("Warinig! Exception \"%s\" caught while doing fit of the histogram no. %d. The fitting process will be continued without it.\n", e.what(), iSliceVar); + } - auto drawOnCanvas = [&](std::vector& canvas, std::function drawer) { + auto drawOnCanvas = [&](std::vector& canvas, const std::function& drawer) { if (nHistograms > 1) { canvas[iCanvas]->cd(iSliceVar - NCanvasesMax * iCanvas + 1); } else { @@ -638,6 +671,14 @@ 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); + + hFitResult->SetBinContent(FitResultStatus, iSliceVar + 1, massFitter->getFitStatus()); + hFitResult->SetBinContent(FitResultCovQual, iSliceVar + 1, massFitter->getCovQual()); + hFitResult->SetBinContent(FitResultEdm, iSliceVar + 1, massFitter->getEDM()); + hFitResult->SetBinContent(FitResultMinNll, iSliceVar + 1, massFitter->getMinNll()); + hFitResult->SetBinContent(FitResultNSgnGCC, iSliceVar + 1, massFitter->getSgnGlobalCorrelCoeff()); + + hCovCorr[iSliceVar] = massFitter->getCovCorrMatrix(); } // save output histograms @@ -654,8 +695,18 @@ void runMassFitter(const std::string& configFileName) } for (int iSliceVar = 0; iSliceVar < nHistograms; iSliceVar++) { + if (iSliceVar == 0) { + outputFile.mkdir("MassHistograms"); + outputFile.mkdir("CovCorrMatrices"); + } + outputFile.cd("MassHistograms"); hMass[iSliceVar]->Write(); + outputFile.cd("CovCorrMatrices"); + if (hCovCorr[iSliceVar] != nullptr) { + hCovCorr[iSliceVar]->Write(Form("hCovCorrMatrix%d", iSliceVar + 1)); + } } + outputFile.cd(); hRawYieldsSignal->Write(); hRawYieldsSignalCounted->Write(); hRawYieldsBkg->Write(); @@ -677,6 +728,7 @@ void runMassFitter(const std::string& configFileName) hRawYieldsDscbNR->Write(); } hFitConfig->Write(); + hFitResult->Write(); outputFile.Close(); @@ -800,9 +852,8 @@ void readJsonVectorFlexible(std::vector& vec, const Document& config, int nHi if (!config.HasMember(fieldName.c_str())) { if (isRequired) { throw std::runtime_error("readJsonVectorFlexible(): missing required field " + fieldName); - } else { - return; } + return; } if (config[fieldName.c_str()].IsArray()) { readJsonVector(vec, config, fieldName); @@ -839,7 +890,11 @@ int main(int argc, const char* argv[]) const std::string configFileName = argv[1]; - runMassFitter(configFileName); - - return 0; + try { + runMassFitter(configFileName); + return 0; + } catch (const std::exception& e) { + printf("Error: Exception \"%s\" caught during runMassFitter() call. Exit.\n", e.what()); + return 1; + } }