diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.cxx b/PWGHF/D2H/Macros/HFInvMassFitter.cxx index 957ec2bc523..843712e2245 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.cxx +++ b/PWGHF/D2H/Macros/HFInvMassFitter.cxx @@ -17,6 +17,7 @@ /// \author Xinye Peng /// \author Biao Zhang /// \author Oleksii Lubynets +/// \author Phil Stahlhut #include "HFInvMassFitter.h" @@ -50,6 +51,7 @@ #include #include #include +#include using namespace RooFit; @@ -74,7 +76,7 @@ HFInvMassFitter::HFInvMassFitter() : mHistoInvMass(nullptr), mNSigmaForSidebands(4.), mNSigmaForSgn(3.), mSigmaSgnErr(0.), - mSigmaSgnDoubleGaus(0.012), + mSigmaSgnDoubleGaus(0.025), mFixedMean(kFALSE), mBoundMean(kFALSE), mBoundReflMean(kFALSE), @@ -103,6 +105,8 @@ HFInvMassFitter::HFInvMassFitter() : mHistoInvMass(nullptr), mFixReflOverSgn(kFALSE), mRooMeanSgn(nullptr), mRooSigmaSgn(nullptr), + mRooSecSigmaSgn(nullptr), + mRooFracDoubleGaus(nullptr), mSgnPdf(nullptr), mBkgPdf(nullptr), mReflPdf(nullptr), @@ -145,7 +149,7 @@ HFInvMassFitter::HFInvMassFitter(const TH1* histoToFit, Double_t minValue, Doubl mNSigmaForSidebands(3.), mNSigmaForSgn(3.), mSigmaSgnErr(0.), - mSigmaSgnDoubleGaus(0.012), + mSigmaSgnDoubleGaus(0.025), mFixedMean(kFALSE), mBoundMean(kFALSE), mBoundReflMean(kFALSE), @@ -174,6 +178,8 @@ HFInvMassFitter::HFInvMassFitter(const TH1* histoToFit, Double_t minValue, Doubl mFixReflOverSgn(kFALSE), mRooMeanSgn(nullptr), mRooSigmaSgn(nullptr), + mRooSecSigmaSgn(nullptr), + mRooFracDoubleGaus(nullptr), mSgnPdf(nullptr), mBkgPdf(nullptr), mReflPdf(nullptr), @@ -207,6 +213,8 @@ HFInvMassFitter::~HFInvMassFitter() delete mHistoTemplateRefl; delete mRooMeanSgn; delete mRooSigmaSgn; + delete mRooSecSigmaSgn; + delete mRooFracDoubleGaus; delete mSgnPdf; delete mBkgPdf; delete mReflPdf; @@ -356,13 +364,14 @@ void HFInvMassFitter::doFit() 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")); RooHist* residualHistogram = mInvMassFrame->residHist("data_c", "Bkg_c"); mResidualFrame->addPlotable(residualHistogram, "P"); mSgnPdf->plotOn(mResidualFrame, Normalization(1.0, RooAbsReal::RelativeExpected), LineColor(kBlue)); } - mass->setRange("bkgForSignificance", mRooMeanSgn->getVal() - mNSigmaForSgn * mRooSigmaSgn->getVal(), mRooMeanSgn->getVal() + mNSigmaForSgn * mRooSigmaSgn->getVal()); + mass->setRange("bkgForSignificance", mRooMeanSgn->getVal() - mNSigmaForSgn * mRooSecSigmaSgn->getVal(), mRooMeanSgn->getVal() + mNSigmaForSgn * mRooSecSigmaSgn->getVal()); bkgIntegral = mBkgPdf->createIntegral(*mass, NormSet(*mass), Range("bkgForSignificance")); mIntegralBkg = bkgIntegral->getValV(); calculateBackground(mBkgYield, mBkgYieldErr); @@ -372,6 +381,9 @@ void HFInvMassFitter::doFit() calculateSignal(mRawYield, mRawYieldErr); countSignal(mRawYieldCounted, mRawYieldCountedErr); calculateSignificance(mSignificance, mSignificanceErr); + // Fit to data ratio + mRatioFrame = mass->frame(Title(Form("%s", mHistoInvMass->GetTitle()))); + calculateFitToDataRatio(); } } @@ -440,10 +452,10 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const workspace.import(*sgnFuncGaus); delete sgnFuncGaus; // signal double Gaussian - RooRealVar sigmaDoubleGaus("sigmaDoubleGaus", "sigma2Gaus", mSigmaSgn, mSigmaSgn - 0.01, mSigmaSgn + 0.01); + RooRealVar sigmaDoubleGaus("sigmaDoubleGaus", "sigma2Gaus", mSigmaSgnDoubleGaus, mSigmaSgnDoubleGaus - 0.003, mSigmaSgnDoubleGaus + 0.003); if (mBoundSigma) { - sigmaDoubleGaus.setMax(mSigmaSgn * (1 + mParamSgn)); - sigmaDoubleGaus.setMin(mSigmaSgn * (1 - mParamSgn)); + sigmaDoubleGaus.setMax(mSigmaSgnDoubleGaus * (1 + mParamSgn)); + sigmaDoubleGaus.setMin(mSigmaSgnDoubleGaus * (1 - mParamSgn)); } if (mFixedSigma) { sigma.setVal(mSigmaSgn); @@ -552,64 +564,81 @@ void HFInvMassFitter::fillWorkspace(RooWorkspace& workspace) const workspace.import(*reflFuncPoly6); delete reflFuncPoly6; } + // draw fit output -void HFInvMassFitter::drawFit(TVirtualPad* pad, Int_t writeFitInfo) +void HFInvMassFitter::drawFit(TVirtualPad* pad, const std::vector& plotLabels, Bool_t writeParInfo) { gStyle->SetOptStat(0); gStyle->SetCanvasColor(0); gStyle->SetFrameFillColor(0); pad->cd(); - if (writeFitInfo > 0) { - auto* textInfoLeft = new TPaveText(0.12, 0.65, 0.47, 0.89, "NDC"); - auto* textInfoRight = new TPaveText(0.6, 0.7, 1., .87, "NDC"); - textInfoLeft->SetBorderSize(0); - textInfoLeft->SetFillStyle(0); - textInfoRight->SetBorderSize(0); - textInfoRight->SetFillStyle(0); - textInfoRight->SetTextColor(kBlue); - textInfoLeft->AddText(Form("S = %.0f #pm %.0f ", mRawYield, mRawYieldErr)); - textInfoLeft->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr)); - if (mTypeOfBkgPdf != NoBkg) { - textInfoLeft->AddText(Form("B (%d#sigma) = %.0f #pm %.0f", mNSigmaForSidebands, mBkgYield, mBkgYieldErr)); - textInfoLeft->AddText(Form("S/B (%d#sigma) = %.4g ", mNSigmaForSidebands, mRawYield / mBkgYield)); - } - if (mReflPdf != nullptr) { - textInfoLeft->AddText(Form("Refl/Sig = %.3f #pm %.3f ", mReflOverSgn, 0.0)); - } - if (mTypeOfBkgPdf != NoBkg) { - textInfoLeft->AddText(Form("Signif (%d#sigma) = %.1f #pm %.1f ", mNSigmaForSidebands, mSignificance, mSignificanceErr)); - textInfoLeft->AddText(Form("#chi^{2} / ndf = %.3f", mChiSquareOverNdfTotal)); - } + // Fit metrics + auto* textFitMetrics = new TPaveText(0.65, 0.7, 0.9, 0.88, "NDC"); + textFitMetrics->SetBorderSize(0); + textFitMetrics->SetFillStyle(0); + textFitMetrics->SetTextSize(0.04); + textFitMetrics->SetTextAlign(33); + 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 (%d#sigma) = %.0f #pm %.0f", mNSigmaForSidebands, mBkgYield, mBkgYieldErr)); + textFitMetrics->AddText(Form("S/B (%d#sigma) = %.4g ", mNSigmaForSidebands, mRawYield / mBkgYield)); + textFitMetrics->AddText(Form("Significance (%d#sigma) = %.1f #pm %.1f ", mNSigmaForSidebands, mSignificance, mSignificanceErr)); + textFitMetrics->AddText(Form("#chi^{2} / ndf = %.3f", mChiSquareOverNdfTotal)); + } + if (mReflPdf != nullptr) { + textFitMetrics->AddText(Form("Refl/Sig = %.3f #pm %.3f ", mReflOverSgn, 0.0)); + } + mInvMassFrame->addObject(textFitMetrics); + // Analysis information + auto* textAnalysisInfo = new TPaveText(0.18, 0.78, 0.35, 0.88, "NDC"); + textAnalysisInfo->SetBorderSize(0); + textAnalysisInfo->SetFillStyle(0); + textAnalysisInfo->SetTextSize(0.05); + textAnalysisInfo->SetTextAlign(13); + for (const auto& label : plotLabels) { + textAnalysisInfo->AddText(label.c_str()); + } + mInvMassFrame->addObject(textAnalysisInfo); + if (writeParInfo) { + // right text box + auto* textSignalPar = new TPaveText(0.18, 0.65, 0.4, 0.75, "NDC"); + textSignalPar->SetBorderSize(0); + textSignalPar->SetFillStyle(0); + textSignalPar->SetTextColor(kBlue); + textSignalPar->SetTextAlign(13); if (mFixedMean) { - textInfoRight->AddText(Form("mean(fixed) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError())); + textSignalPar->AddText(Form("mean(fixed) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError())); } else { - textInfoRight->AddText(Form("mean(free) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError())); + textSignalPar->AddText(Form("mean(free) = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError())); } if (mTypeOfSgnPdf == DoubleGaus) { - auto const& baseSigmaSgn = mWorkspace->var("sigma"); + if (mFixedSigma) { + textSignalPar->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + } else { + textSignalPar->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + } if (mFixedSigmaDoubleGaus) { - textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError())); - textInfoRight->AddText(Form("sigma 2(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textSignalPar->AddText(Form("sigma 2(fixed) = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); } else { - textInfoRight->AddText(Form("sigma(free) = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError())); - textInfoRight->AddText(Form("sigma 2(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textSignalPar->AddText(Form("sigma 2(free) = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); } } else if (mFixedSigma) { - textInfoRight->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textSignalPar->AddText(Form("sigma(fixed) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); } else { - textInfoRight->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); - } - mInvMassFrame->addObject(textInfoLeft); - mInvMassFrame->addObject(textInfoRight); - mInvMassFrame->GetYaxis()->SetTitleOffset(1.8); - gPad->SetLeftMargin(0.15); - mInvMassFrame->GetYaxis()->SetTitle(Form("%s", mHistoInvMass->GetYaxis()->GetTitle())); - mInvMassFrame->GetXaxis()->SetTitle(Form("%s", mHistoInvMass->GetXaxis()->GetTitle())); - mInvMassFrame->Draw(); - highlightPeakRegion(mInvMassFrame); - if (mHistoTemplateRefl != nullptr) { - mReflFrame->Draw("same"); + textSignalPar->AddText(Form("sigma(free) = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); } + mInvMassFrame->addObject(textSignalPar); + } + mInvMassFrame->GetXaxis()->SetTitleOffset(1.2); + mInvMassFrame->GetYaxis()->SetTitleOffset(1.8); + gPad->SetLeftMargin(0.15); + mInvMassFrame->GetYaxis()->SetTitle(Form("%s", mHistoInvMass->GetYaxis()->GetTitle())); + mInvMassFrame->GetXaxis()->SetTitle(Form("%s", mHistoInvMass->GetXaxis()->GetTitle())); + mInvMassFrame->Draw(); + highlightPeakRegion(mInvMassFrame); + if (mHistoTemplateRefl) { + mReflFrame->Draw("same"); } } @@ -626,9 +655,8 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad) textInfo->AddText(Form("S_{count} = %.0f #pm %.0f ", mRawYieldCounted, mRawYieldCountedErr)); textInfo->AddText(Form("mean = %.3f #pm %.3f", mRooMeanSgn->getVal(), mRooMeanSgn->getError())); if (mTypeOfSgnPdf == DoubleGaus) { - auto const& baseSigmaSgn = mWorkspace->var("sigma"); - textInfo->AddText(Form("sigma = %.3f #pm %.3f", baseSigmaSgn->getVal(), baseSigmaSgn->getError())); - textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); + textInfo->AddText(Form("sigma 2 = %.3f #pm %.3f", mRooSecSigmaSgn->getVal(), mRooSecSigmaSgn->getError())); } else { textInfo->AddText(Form("sigma = %.3f #pm %.3f", mRooSigmaSgn->getVal(), mRooSigmaSgn->getError())); } @@ -637,6 +665,24 @@ void HFInvMassFitter::drawResidual(TVirtualPad* pad) highlightPeakRegion(mResidualFrame); } +// draw ratio on canvas +void HFInvMassFitter::drawRatio(TVirtualPad* pad) +{ + pad->cd(); + mRatioFrame->GetXaxis()->SetTitleOffset(1.2); + mRatioFrame->GetYaxis()->SetTitleOffset(1.5); + mRatioFrame->GetYaxis()->SetTitle("Fit / Data"); + double xMin = mRatioFrame->GetXaxis()->GetXmin(); + double xMax = mRatioFrame->GetXaxis()->GetXmax(); + auto* line = new TLine(xMin, 1.0, xMax, 1.0); + line->SetLineColor(kGray); + line->SetLineStyle(2); + line->SetLineWidth(2); + mRatioFrame->addObject(line); + mRatioFrame->Draw(); + highlightPeakRegion(mRatioFrame); +} + // draw peak region with vertical lines void HFInvMassFitter::highlightPeakRegion(const RooPlot* plot, Color_t color, Width_t width, Style_t style) const { @@ -646,7 +692,7 @@ void HFInvMassFitter::highlightPeakRegion(const RooPlot* plot, Color_t color, Wi double const yMin = plot->GetMinimum(); double const yMax = plot->GetMaximum(); const Double_t mean = mRooMeanSgn->getVal(); - const Double_t sigma = mRooSigmaSgn->getVal(); + const Double_t sigma = mRooSecSigmaSgn->getVal(); const Double_t minForSgn = mean - mNSigmaForSidebands * sigma; const Double_t maxForSgn = mean + mNSigmaForSidebands * sigma; auto* leftLine = new TLine(minForSgn, yMin, minForSgn, yMax); @@ -671,7 +717,7 @@ void HFInvMassFitter::drawReflection(TVirtualPad* pad) void HFInvMassFitter::countSignal(Double_t& signal, Double_t& signalErr) const { const Double_t mean = mRooMeanSgn->getVal(); - const Double_t sigma = mRooSigmaSgn->getVal(); + const Double_t sigma = mRooSecSigmaSgn->getVal(); const Double_t minForSgn = mean - mNSigmaForSidebands * sigma; const Double_t maxForSgn = mean + mNSigmaForSidebands * sigma; const Int_t binForMinSgn = mHistoInvMass->FindBin(minForSgn); @@ -763,21 +809,27 @@ RooAbsPdf* HFInvMassFitter::createSignalFitFunction(RooWorkspace* workspace) case 0: { sgnPdf = workspace->pdf("sgnFuncGaus"); mRooSigmaSgn = workspace->var("sigma"); + mRooSecSigmaSgn = workspace->var("sigma"); mRooMeanSgn = workspace->var("mean"); } break; case 1: { sgnPdf = workspace->pdf("sgnFuncDoubleGaus"); - mRooSigmaSgn = workspace->var("sigmaDoubleGaus"); + mRooSigmaSgn = workspace->var("sigma"); + mRooSecSigmaSgn = workspace->var("sigmaDoubleGaus"); mRooMeanSgn = workspace->var("mean"); + mRooFracDoubleGaus = workspace->var("fracDoubleGaus"); } break; case 2: { sgnPdf = workspace->pdf("sgnFuncGausRatio"); - mRooSigmaSgn = workspace->var("sigmaDoubleGausRatio"); + mRooSigmaSgn = workspace->var("sigma"); + mRooSecSigmaSgn = workspace->var("sigmaDoubleGausRatio"); mRooMeanSgn = workspace->var("mean"); + mRooFracDoubleGaus = workspace->var("fracDoubleGausRatio"); } break; case 3: { sgnPdf = workspace->pdf("sgnFuncDoublePeak"); - mRooSigmaSgn = workspace->var("sigmaSec"); + mRooSigmaSgn = workspace->var("sigma"); + mRooSecSigmaSgn = workspace->var("sigmaSec"); mRooMeanSgn = workspace->var("meanSec"); } break; default: @@ -818,6 +870,40 @@ void HFInvMassFitter::plotRefl(RooAbsPdf* pdf) pdf->plotOn(mInvMassFrame, Components(namesOfReflPdf.at(mTypeOfReflPdf).c_str()), Name("Refl_c"), LineColor(kGreen)); } +// Calculate fit to data ratio +void HFInvMassFitter::calculateFitToDataRatio() const +{ + if (!mInvMassFrame) + return; + + // Get the data and fit curves from the frame + auto* dataHist = dynamic_cast(mInvMassFrame->findObject("data_c")); + auto* fitCurve = dynamic_cast(mInvMassFrame->findObject("Tot_c")); // or the relevant fit curve + + if (!dataHist || !fitCurve) + return; + + RooHist* ratioHist = new RooHist(); + + for (int i = 0; i < dataHist->GetN(); ++i) { + double x, dataY, dataErr; + dataHist->GetPoint(i, x, dataY); + dataErr = dataHist->GetErrorY(i); + + double fitY = fitCurve->Eval(x); + + double ratio = dataY != 0 ? fitY / dataY : 0; + double err = dataY != 0 ? ratio * dataErr / dataY : 0; + + ratioHist->SetPoint(i, x, ratio); + ratioHist->SetPointError(i, 0, 0, err, err); + } + + mRatioFrame->addPlotable(ratioHist, "P"); + mRatioFrame->SetMinimum(0.5); + mRatioFrame->SetMaximum(1.5); +} + // Fix reflection pdf void HFInvMassFitter::setReflFuncFixed() { diff --git a/PWGHF/D2H/Macros/HFInvMassFitter.h b/PWGHF/D2H/Macros/HFInvMassFitter.h index 1a245fe0db5..daa00ec0b09 100644 --- a/PWGHF/D2H/Macros/HFInvMassFitter.h +++ b/PWGHF/D2H/Macros/HFInvMassFitter.h @@ -17,6 +17,7 @@ /// \author Xinye Peng /// \author Biao Zhang /// \author Oleksii Lubynets +/// \author Phil Stahlhut #ifndef PWGHF_D2H_MACROS_HFINVMASSFITTER_H_ #define PWGHF_D2H_MACROS_HFINVMASSFITTER_H_ @@ -212,7 +213,12 @@ class HFInvMassFitter : public TNamed [[nodiscard]] Double_t getMeanUncertainty() const { return mRooMeanSgn->getError(); } [[nodiscard]] Double_t getSigma() const { return mRooSigmaSgn->getVal(); } [[nodiscard]] Double_t getSigmaUncertainty() const { return mRooSigmaSgn->getError(); } + [[nodiscard]] Double_t getSecSigma() const { return mRooSecSigmaSgn->getVal(); } + [[nodiscard]] Double_t getSecSigmaUncertainty() const { return mRooSecSigmaSgn->getError(); } + [[nodiscard]] Double_t getFracDoubleGaus() const { return mRooFracDoubleGaus->getVal(); } + [[nodiscard]] Double_t getFracDoubleGausUncertainty() const { return mRooFracDoubleGaus->getError(); } [[nodiscard]] Double_t getReflOverSig() const + { if (mReflPdf != nullptr) { return mReflOverSgn; @@ -224,8 +230,10 @@ class HFInvMassFitter : public TNamed void calculateBackground(Double_t& bkg, Double_t& bkgErr) const; void calculateSignificance(Double_t& significance, Double_t& significanceErr) const; void checkForSignal(Double_t& estimatedSignal); - void drawFit(TVirtualPad* c, Int_t writeFitInfo = 2); + void calculateFitToDataRatio() const; + void drawFit(TVirtualPad* c, const std::vector& plotLabels, Bool_t writeParInfo = true); void drawResidual(TVirtualPad* c); + void drawRatio(TVirtualPad* c); void drawReflection(TVirtualPad* c); private: @@ -282,6 +290,8 @@ class HFInvMassFitter : public TNamed Bool_t mFixReflOverSgn; /// switch for fix refl/signal RooRealVar* mRooMeanSgn; /// mean for gaussian of signal RooRealVar* mRooSigmaSgn; /// sigma for gaussian of signal + RooRealVar* mRooSecSigmaSgn; /// second sigma for composite gaussian of signal + RooRealVar* mRooFracDoubleGaus; /// fraction of second gaussian for composite gaussian of signal RooAbsPdf* mSgnPdf; /// signal fit function RooAbsPdf* mBkgPdf; /// background fit function RooAbsPdf* mReflPdf; /// reflection fit function @@ -293,6 +303,7 @@ class HFInvMassFitter : public TNamed RooPlot* mReflFrame; /// reflection frame RooPlot* mReflOnlyFrame; /// reflection frame plot on reflection only RooPlot* mResidualFrame; /// residual frame + RooPlot* mRatioFrame; /// fit/data ratio frame RooPlot* mResidualFrameForCalculation; RooWorkspace* mWorkspace; /// workspace Double_t mIntegralHisto; /// integral of histogram to fit diff --git a/PWGHF/D2H/Macros/config_massfitter.json b/PWGHF/D2H/Macros/config_massfitter.json index 43e67b76f5d..393fb8aa0f3 100644 --- a/PWGHF/D2H/Macros/config_massfitter.json +++ b/PWGHF/D2H/Macros/config_massfitter.json @@ -1,5 +1,7 @@ { "IsMC": false, + "WriteSignalPar": true, + "_WriteSignalPar": "write signal parameter values on mass canvas", "InFileName": "hMasses.root", "_InFileName": "download example file here https://cernbox.cern.ch/s/uoWicuRAVGArDsV", "ReflFileName": "", @@ -35,6 +37,8 @@ "Dstar", "XicToXiPiPi" ], + "CollisionSystem": "pp #sqrt{#it{s}} = 13.6 TeV", + "_CollisionSystems": "Label for mass canvases", "EnableRefl": false, "FixSigma": false, "SigmaFile": "", @@ -57,7 +61,7 @@ 0, 0 ], - "_FixMeanManual": "fix mean mannually", + "_FixMeanManual": "fix mean manually", "FixSecondSigma": false, "SecondSigmaFile": "", "_SecondSigmaFile": "fix second sigma for double gauss from file", @@ -69,6 +73,17 @@ 0 ], "_FixSecondSigmaManual": "fix second sigma for double gauss manually", + "FixFracDoubleGaus": false, + "FracDoubleGausFile": "", + "_FracDoubleGausFile": "fix fraction of sigmas for double gauss from file", + "FixFracDoubleGausManual": [ + 0, + 0, + 0, + 0, + 0 + ], + "_FixFracDoubleGausManual": "fix fraction of sigmas for double gauss manually", "SliceVarName": "T", "SliceVarUnit": "ps", "_SliceVarName, _SliceVarUnit": "e.g. pT, GeV/c or something else depending on user's needs", diff --git a/PWGHF/D2H/Macros/runMassFitter.C b/PWGHF/D2H/Macros/runMassFitter.C index 0b926b321f0..e4baa2a856a 100644 --- a/PWGHF/D2H/Macros/runMassFitter.C +++ b/PWGHF/D2H/Macros/runMassFitter.C @@ -17,6 +17,7 @@ /// \author Xinye Peng /// \author Biao Zhang /// \author Oleksii Lubynets +/// \author Phil Stahlhut #if !defined(__CINT__) || defined(__CLING__) @@ -42,6 +43,7 @@ #include #include #include // std::string +#include #include #include // std::vector @@ -88,10 +90,12 @@ int runMassFitter(const TString& configFileName) fclose(configFile); Bool_t const isMc = config["IsMC"].GetBool(); + Bool_t const writeSignalPar = config["WriteSignalPar"].GetBool(); + TString const particleName = config["Particle"].GetString(); + TString const collisionSystem = config["CollisionSystem"].GetString(); TString const inputFileName = config["InFileName"].GetString(); TString const reflFileName = config["ReflFileName"].GetString(); TString outputFileName = config["OutFileName"].GetString(); - TString const particleName = config["Particle"].GetString(); std::vector inputHistoName; std::vector promptHistoName; @@ -105,9 +109,10 @@ int runMassFitter(const TString& configFileName) std::vector sliceVarMax; std::vector massMin; std::vector massMax; + std::vector fixMeanManual; std::vector fixSigmaManual; std::vector fixSecondSigmaManual; - std::vector fixMeanManual; + std::vector fixFracDoubleGausManual; std::vector nRebin; std::vector bkgFuncConfig; std::vector sgnFuncConfig; @@ -130,12 +135,15 @@ int runMassFitter(const TString& configFileName) const Value& fdSecPeakHistoNameValue = config["FDSecPeakHistoName"]; parseStringArray(fdSecPeakHistoNameValue, fdSecPeakHistoName); - const bool fixSigma = config["FixSigma"].GetBool(); - const std::string sigmaFile = config["SigmaFile"].GetString(); - const bool fixMean = config["FixMean"].GetBool(); const std::string meanFile = config["MeanFile"].GetString(); + const Value& fixMeanManualValue = config["FixMeanManual"]; + readArray(fixMeanManualValue, fixMeanManual); + + const bool fixSigma = config["FixSigma"].GetBool(); + const std::string sigmaFile = config["SigmaFile"].GetString(); + const Value& fixSigmaManualValue = config["FixSigmaManual"]; readArray(fixSigmaManualValue, fixSigmaManual); @@ -145,8 +153,11 @@ int runMassFitter(const TString& configFileName) const Value& fixSecondSigmaManualValue = config["FixSecondSigmaManual"]; readArray(fixSecondSigmaManualValue, fixSecondSigmaManual); - const Value& fixMeanManualValue = config["FixMeanManual"]; - readArray(fixMeanManualValue, fixMeanManual); + const bool fixFracDoubleGaus = config["FixFracDoubleGaus"].GetBool(); + const std::string fracDoubleGausFile = config["FracDoubleGausFile"].GetString(); + + const Value& fixFracDoubleGausManualValue = config["FixFracDoubleGausManual"]; + readArray(fixFracDoubleGausManualValue, fixFracDoubleGausManual); sliceVarName = config["SliceVarName"].GetString(); sliceVarUnit = config["SliceVarUnit"].GetString(); @@ -200,19 +211,21 @@ int runMassFitter(const TString& configFileName) sgnFunc[iSliceVar] = sgnFuncConfig[iSliceVar]; } - std::map> particles{ - {"Dplus", {"K#pi#pi", "D+"}}, - {"D0", {"K#pi", "D0"}}, - {"Ds", {"KK#pi", "D_s+"}}, - {"LcToPKPi", {"pK#pi", "Lambda_c+"}}, - {"LcToPK0s", {"pK^{0}_{s}", "Lambda_c+"}}, - {"Dstar", {"D^{0}pi^{+}", "D*+"}}, - {"XicToXiPiPi", {"#Xi#pi#pi", "Xi_c+"}}}; + std::map> particles{ + {"Dplus", {"K#pi#pi", "D+", "D^{+} #rightarrow K^{-}#pi^{+}#pi^{+} + c.c."}}, + {"D0", {"K#pi", "D0", "D^{0} #rightarrow K^{-}#pi^{+} + c.c."}}, + {"Ds", {"KK#pi", "D_s+", "D_{s}^{+} #rightarrow K^{-}K^{+}#pi^{+} + c.c."}}, + {"LcToPKPi", {"pK#pi", "Lambda_c+", "#Lambda_{c}^{+} #rightarrow pK^{-}#pi^{+} + c.c."}}, + {"LcToPK0s", {"pK^{0}_{s}", "Lambda_c+", "#Lambda_{c}^{+} #rightarrow pK^{0}_{s} + c.c."}}, + {"Dstar", {"D^{0}pi^{+}", "D*+", "D^{*+} #rightarrow D^{0}#pi^{+} + c.c."}}, + {"XicToXiPiPi", {"#Xi#pi#pi", "Xi_c+", "#Xi_{c}^{+} #rightarrow #Xi^{-}#pi^{+}#pi^{+} + c.c."}}}; if (particles.find(particleName.Data()) == particles.end()) { throw std::runtime_error("ERROR: only Dplus, D0, Ds, LcToPKPi, LcToPK0s, Dstar and XicToXiPiPi particles supported! Exit"); } - const TString massAxisTitle = "#it{M}(" + particles[particleName.Data()].first + ") (GeV/#it{c}^{2})"; - const double massPDG = TDatabasePDG::Instance()->GetParticle(particles[particleName.Data()].second.c_str())->Mass(); + const auto& particleTuple = particles[particleName.Data()]; + const TString massAxisTitle = "#it{M}(" + std::get<0>(particleTuple) + ") (GeV/#it{c}^{2})"; + const double massPDG = TDatabasePDG::Instance()->GetParticle(std::get<1>(particleTuple).c_str())->Mass(); + const std::vector plotLabels = {std::get<2>(particleTuple), collisionSystem.Data()}; // load inv-mass histograms auto* inputFile = TFile::Open(inputFileName.Data()); @@ -262,34 +275,18 @@ int runMassFitter(const TString& configFileName) inputFile->Close(); // 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* hRawYieldsSigma = new TH1D( - "hRawYieldsSigma", ";" + sliceVarName + "(" + sliceVarUnit + ");width (GeV/#it{c}^{2})", - nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsMean = new TH1D( - "hRawYieldsMean", ";" + sliceVarName + "(" + sliceVarUnit + ");mean (GeV/#it{c}^{2})", - nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsSignificance = new TH1D( - "hRawYieldsSignificance", - ";" + sliceVarName + "(" + sliceVarUnit + ");significance (3#sigma)", nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsSgnOverBkg = - new TH1D("hRawYieldsSgnOverBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");S/B (3#sigma)", - nSliceVarBins, sliceVarLimits.data()); - auto* hRawYieldsBkg = - new TH1D("hRawYieldsBkg", ";" + sliceVarName + "(" + sliceVarUnit + ");Background (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* 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()); const Int_t nConfigsToSave = 6; auto* hFitConfig = new TH2F("hfitConfig", "Fit Configurations", nConfigsToSave, 0, 6, nSliceVarBins, sliceVarLimits.data()); @@ -304,14 +301,16 @@ int runMassFitter(const TString& configFileName) setHistoStyle(hRawYieldsSignal); setHistoStyle(hRawYieldsSignalCounted); - setHistoStyle(hRawYieldsSigma); - setHistoStyle(hRawYieldsMean); - setHistoStyle(hRawYieldsSignificance); - setHistoStyle(hRawYieldsSgnOverBkg); setHistoStyle(hRawYieldsBkg); + setHistoStyle(hRawYieldsSgnOverBkg); + setHistoStyle(hRawYieldsSignificance); setHistoStyle(hRawYieldsChiSquareBkg); setHistoStyle(hRawYieldsChiSquareTotal); setHistoStyle(hReflectionOverSignal, kRed + 1); + setHistoStyle(hRawYieldsMean); + setHistoStyle(hRawYieldsSigma); + setHistoStyle(hRawYieldsSecSigma); + setHistoStyle(hRawYieldsFracDoubleGaus); auto getHistToFix = [&nSliceVarBins](bool const& isFix, std::vector const& fixManual, std::string const& fixFileName, std::string const& var) -> TH1* { TH1* histToFix = nullptr; @@ -335,7 +334,8 @@ int runMassFitter(const TString& configFileName) TH1* hSigmaToFix = getHistToFix(fixSigma, fixSigmaManual, sigmaFile, "Sigma"); TH1* hMeanToFix = getHistToFix(fixMean, fixMeanManual, meanFile, "Mean"); - TH1* hSecondSigmaToFix = getHistToFix(fixSecondSigma, fixSecondSigmaManual, secondSigmaFile, "Sigma"); + TH1* hSecondSigmaToFix = getHistToFix(fixSecondSigma, fixSecondSigmaManual, secondSigmaFile, "SecSigma"); + TH1* hFracDoubleGausToFix = getHistToFix(fixFracDoubleGaus, fixFracDoubleGausManual, fracDoubleGausFile, "FracDoubleGaus"); // fit histograms @@ -353,6 +353,7 @@ int runMassFitter(const TString& configFileName) const Int_t nCanvases = std::ceil(static_cast(nSliceVarBins) / 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; @@ -363,9 +364,16 @@ int runMassFitter(const TString& configFileName) canvasResiduals[iCanvas] = new TCanvas(Form("canvasResiduals%d", iCanvas), Form("canvasResiduals%d", iCanvas), canvasSize[0], canvasSize[1]); divideCanvas(canvasResiduals[iCanvas], nPads); - canvasRefl[iCanvas] = new TCanvas(Form("canvasRefl%d", iCanvas), Form("canvasRefl%d", iCanvas), - canvasSize[0], canvasSize[1]); - divideCanvas(canvasRefl[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); + } } for (unsigned int iSliceVar = 0; iSliceVar < nSliceVarBins; iSliceVar++) { @@ -402,32 +410,44 @@ int runMassFitter(const TString& configFileName) canvasMass[iCanvas]->cd(); } - massFitter->drawFit(gPad); + massFitter->drawFit(gPad, plotLabels, writeSignalPar); const Double_t rawYield = massFitter->getRawYield(); const Double_t rawYieldErr = massFitter->getRawYieldError(); const Double_t rawYieldCounted = massFitter->getRawYieldCounted(); const Double_t rawYieldCountedErr = massFitter->getRawYieldCountedError(); - - const Double_t sigma = massFitter->getSigma(); - const Double_t sigmaErr = massFitter->getSigmaUncertainty(); - const Double_t mean = massFitter->getMean(); - const Double_t meanErr = massFitter->getMeanUncertainty(); const Double_t reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg(); const Double_t reducedChiSquareTotal = massFitter->getChiSquareOverNDFTotal(); + const Double_t mean = massFitter->getMean(); + const Double_t meanErr = massFitter->getMeanUncertainty(); + const Double_t sigma = massFitter->getSigma(); + const Double_t sigmaErr = massFitter->getSigmaUncertainty(); hRawYieldsSignal->SetBinContent(iSliceVar + 1, rawYield); hRawYieldsSignal->SetBinError(iSliceVar + 1, rawYieldErr); hRawYieldsSignalCounted->SetBinContent(iSliceVar + 1, rawYieldCounted); hRawYieldsSignalCounted->SetBinError(iSliceVar + 1, rawYieldCountedErr); - hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma); - hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr); - hRawYieldsMean->SetBinContent(iSliceVar + 1, mean); - hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr); hRawYieldsChiSquareBkg->SetBinContent(iSliceVar + 1, reducedChiSquareBkg); hRawYieldsChiSquareBkg->SetBinError(iSliceVar + 1, 0.); hRawYieldsChiSquareTotal->SetBinContent(iSliceVar + 1, reducedChiSquareTotal); hRawYieldsChiSquareTotal->SetBinError(iSliceVar + 1, 0.); + hRawYieldsMean->SetBinContent(iSliceVar + 1, mean); + hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr); + hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma); + hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr); + + if (sgnFunc[iSliceVar] != HFInvMassFitter::SingleGaus) { + const Double_t secSigma = massFitter->getSecSigma(); + const Double_t secSigmaErr = massFitter->getSecSigmaUncertainty(); + hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma); + hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr); + } + if (sgnFunc[iSliceVar] == HFInvMassFitter::DoubleGaus || sgnFunc[iSliceVar] == HFInvMassFitter::DoubleGausSigmaRatioPar) { + const Double_t fracDoubleGaus = massFitter->getFracDoubleGaus(); + const Double_t fracDoubleGausErr = massFitter->getFracDoubleGausUncertainty(); + hRawYieldsFracDoubleGaus->SetBinContent(iSliceVar + 1, fracDoubleGaus); + hRawYieldsFracDoubleGaus->SetBinError(iSliceVar + 1, fracDoubleGausErr); + } } else { HFInvMassFitter* massFitter; massFitter = new HFInvMassFitter(hMassForFit[iSliceVar], massMin[iSliceVar], massMax[iSliceVar], @@ -460,6 +480,7 @@ int runMassFitter(const TString& configFileName) 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"); if (enableRefl) { reflOverSgn = hMassForSgn[iSliceVar]->Integral(hMassForSgn[iSliceVar]->FindBin(massMin[iSliceVar] * 1.0001), hMassForSgn[iSliceVar]->FindBin(massMax[iSliceVar] * 0.999)); @@ -474,35 +495,49 @@ int runMassFitter(const TString& configFileName) const double rawYieldErr = massFitter->getRawYieldError(); const double rawYieldCounted = massFitter->getRawYieldCounted(); const double rawYieldCountedErr = massFitter->getRawYieldCountedError(); - const double sigma = massFitter->getSigma(); - const double sigmaErr = massFitter->getSigmaUncertainty(); - const double mean = massFitter->getMean(); - const double meanErr = massFitter->getMeanUncertainty(); - const double reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg(); - const double reducedChiSquareTotal = massFitter->getChiSquareOverNDFTotal(); - const double significance = massFitter->getSignificance(); - const double significanceErr = massFitter->getSignificanceError(); const double bkg = massFitter->getBkgYield(); const double bkgErr = massFitter->getBkgYieldError(); + const double significance = massFitter->getSignificance(); + const double significanceErr = massFitter->getSignificanceError(); + const double reducedChiSquareBkg = massFitter->getChiSquareOverNDFBkg(); + const double reducedChiSquareTotal = massFitter->getChiSquareOverNDFTotal(); + const double mean = massFitter->getMean(); + const double meanErr = massFitter->getMeanUncertainty(); + const double sigma = massFitter->getSigma(); + const double sigmaErr = massFitter->getSigmaUncertainty(); hRawYieldsSignal->SetBinContent(iSliceVar + 1, rawYield); hRawYieldsSignal->SetBinError(iSliceVar + 1, rawYieldErr); hRawYieldsSignalCounted->SetBinContent(iSliceVar + 1, rawYieldCounted); hRawYieldsSignalCounted->SetBinError(iSliceVar + 1, rawYieldCountedErr); - hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma); - hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr); - hRawYieldsMean->SetBinContent(iSliceVar + 1, mean); - hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr); - hRawYieldsSignificance->SetBinContent(iSliceVar + 1, significance); - hRawYieldsSignificance->SetBinError(iSliceVar + 1, significanceErr); - hRawYieldsSgnOverBkg->SetBinContent(iSliceVar + 1, rawYield / bkg); - hRawYieldsSgnOverBkg->SetBinError(iSliceVar + 1, rawYield / bkg * std::sqrt(rawYieldErr / rawYield * rawYieldErr / rawYield + bkgErr / bkg * bkgErr / bkg)); hRawYieldsBkg->SetBinContent(iSliceVar + 1, bkg); hRawYieldsBkg->SetBinError(iSliceVar + 1, bkgErr); + hRawYieldsSgnOverBkg->SetBinContent(iSliceVar + 1, rawYield / bkg); + hRawYieldsSgnOverBkg->SetBinError(iSliceVar + 1, rawYield / bkg * std::sqrt(rawYieldErr / rawYield * rawYieldErr / rawYield + bkgErr / bkg * bkgErr / bkg)); + hRawYieldsSignificance->SetBinContent(iSliceVar + 1, significance); + hRawYieldsSignificance->SetBinError(iSliceVar + 1, significanceErr); hRawYieldsChiSquareBkg->SetBinContent(iSliceVar + 1, reducedChiSquareBkg); hRawYieldsChiSquareBkg->SetBinError(iSliceVar + 1, 1.e-20); hRawYieldsChiSquareTotal->SetBinContent(iSliceVar + 1, reducedChiSquareTotal); hRawYieldsChiSquareTotal->SetBinError(iSliceVar + 1, 1.e-20); + hRawYieldsMean->SetBinContent(iSliceVar + 1, mean); + hRawYieldsMean->SetBinError(iSliceVar + 1, meanErr); + hRawYieldsSigma->SetBinContent(iSliceVar + 1, sigma); + hRawYieldsSigma->SetBinError(iSliceVar + 1, sigmaErr); + + if (sgnFunc[iSliceVar] != HFInvMassFitter::SingleGaus) { + const double secSigma = massFitter->getSecSigma(); + const double secSigmaErr = massFitter->getSecSigmaUncertainty(); + hRawYieldsSecSigma->SetBinContent(iSliceVar + 1, secSigma); + hRawYieldsSecSigma->SetBinError(iSliceVar + 1, secSigmaErr); + } + if (sgnFunc[iSliceVar] == HFInvMassFitter::DoubleGaus || sgnFunc[iSliceVar] == HFInvMassFitter::DoubleGausSigmaRatioPar) { + const double fracDoubleGaus = massFitter->getFracDoubleGaus(); + const double fracDoubleGausErr = massFitter->getFracDoubleGausUncertainty(); + hRawYieldsFracDoubleGaus->SetBinContent(iSliceVar + 1, fracDoubleGaus); + hRawYieldsFracDoubleGaus->SetBinError(iSliceVar + 1, fracDoubleGausErr); + } + if (enableRefl) { hReflectionOverSignal->SetBinContent(iSliceVar + 1, reflOverSgn); if (nSliceVarBins > 1) { @@ -520,18 +555,29 @@ int runMassFitter(const TString& configFileName) } else { canvasMass[iCanvas]->cd(); } - massFitter->drawFit(gPad); + massFitter->drawFit(gPad, plotLabels, writeSignalPar); canvasMass[iCanvas]->Modified(); canvasMass[iCanvas]->Update(); - if (nSliceVarBins > 1) { - canvasResiduals[iCanvas]->cd(iSliceVar - nCanvasesMax * iCanvas + 1); - } else { - canvasResiduals[iCanvas]->cd(); + if (bkgFunc[iSliceVar] != HFInvMassFitter::NoBkg) { + if (nSliceVarBins > 1) { + canvasResiduals[iCanvas]->cd(iSliceVar - nCanvasesMax * iCanvas + 1); + } else { + canvasResiduals[iCanvas]->cd(); + } + massFitter->drawResidual(gPad); + canvasResiduals[iCanvas]->Modified(); + canvasResiduals[iCanvas]->Update(); + + if (nSliceVarBins > 1) { + canvasRatio[iCanvas]->cd(iSliceVar - nCanvasesMax * iCanvas + 1); + } else { + canvasRatio[iCanvas]->cd(); + } + massFitter->drawRatio(gPad); + canvasRatio[iCanvas]->Modified(); + canvasRatio[iCanvas]->Update(); } - massFitter->drawResidual(gPad); - canvasResiduals[iCanvas]->Modified(); - canvasResiduals[iCanvas]->Update(); } hFitConfig->SetBinContent(1, iSliceVar + 1, massMin[iSliceVar]); @@ -554,7 +600,10 @@ int runMassFitter(const TString& configFileName) canvasMass[iCanvas]->Write(); if (!isMc) { canvasResiduals[iCanvas]->Write(); - canvasRefl[iCanvas]->Write(); + canvasRatio[iCanvas]->Write(); + if (enableRefl) { + canvasRefl[iCanvas]->Write(); + } } } @@ -563,13 +612,15 @@ int runMassFitter(const TString& configFileName) } hRawYieldsSignal->Write(); hRawYieldsSignalCounted->Write(); - hRawYieldsSigma->Write(); - hRawYieldsMean->Write(); - hRawYieldsSignificance->Write(); - hRawYieldsSgnOverBkg->Write(); hRawYieldsBkg->Write(); + hRawYieldsSgnOverBkg->Write(); + hRawYieldsSignificance->Write(); hRawYieldsChiSquareBkg->Write(); hRawYieldsChiSquareTotal->Write(); + hRawYieldsMean->Write(); + hRawYieldsSigma->Write(); + hRawYieldsSecSigma->Write(); + hRawYieldsFracDoubleGaus->Write(); if (enableRefl) { hReflectionOverSignal->Write(); } @@ -580,6 +631,8 @@ int runMassFitter(const TString& configFileName) outputFileName.ReplaceAll(".root", ".pdf"); TString outputFileNameResidual = outputFileName; outputFileNameResidual.ReplaceAll(".pdf", "_Residuals.pdf"); + TString outputFileRatio = outputFileName; + outputFileRatio.ReplaceAll(".pdf", "_Ratio.pdf"); for (int iCanvas = 0; iCanvas < nCanvases; iCanvas++) { if (iCanvas == 0 && nCanvases > 1) { canvasMass[iCanvas]->SaveAs(Form("%s[", outputFileName.Data())); @@ -589,6 +642,7 @@ int runMassFitter(const TString& configFileName) canvasMass[iCanvas]->SaveAs(Form("%s]", outputFileName.Data())); } if (!isMc) { + // residuals if (iCanvas == 0 && nCanvases > 1) { canvasResiduals[iCanvas]->SaveAs(Form("%s[", outputFileNameResidual.Data())); } @@ -596,6 +650,14 @@ int runMassFitter(const TString& configFileName) if (iCanvas == nCanvases - 1 && nCanvases > 1) { canvasResiduals[iCanvas]->SaveAs(Form("%s]", outputFileNameResidual.Data())); } + // ratio + if (iCanvas == 0 && nCanvases > 1) { + canvasRatio[iCanvas]->SaveAs(Form("%s[", outputFileRatio.Data())); + } + canvasRatio[iCanvas]->SaveAs(outputFileRatio.Data()); + if (iCanvas == nCanvases - 1 && nCanvases > 1) { + canvasRatio[iCanvas]->SaveAs(Form("%s]", outputFileRatio.Data())); + } } } return 0; @@ -613,13 +675,9 @@ void setHistoStyle(TH1* histo, Color_t color, Size_t markerSize) void divideCanvas(TCanvas* canvas, int nSliceVarBins) { - const int rectangularSideMin = std::floor(std::sqrt(nSliceVarBins)); - constexpr int RectangularSidesDiffMax = 2; - for (int rectangularSidesDiff = 0; rectangularSidesDiff < RectangularSidesDiffMax; ++rectangularSidesDiff) { - if (rectangularSideMin * (rectangularSideMin + rectangularSidesDiff) >= nSliceVarBins) { - canvas->Divide(rectangularSideMin + rectangularSidesDiff, rectangularSideMin); - } - } + int nCols = std::ceil(std::sqrt(nSliceVarBins)); + int nRows = std::ceil(static_cast(nSliceVarBins) / nCols); + canvas->Divide(nCols, nRows); } int main(int argc, const char* argv[])