LightFlavorRatios/unbinnedMaxLikelihood_fitter/RDataframeToRoofit.C

0001 #ifndef RDATAFRAMETOROOFIT_C
0002 #define RDATAFRAMETOROOFIT_C
0003
0004 #include <algorithm>
0005 #include <cmath>
0006 #include <chrono>
0007 #include <iostream>
0008 #include <memory>
0009 #include <string>
0010 #include <vector>
0011
0012 #ifndef __CINT__
0013 #include <RooGlobalFunc.h>
0014 #endif
0015
0016 #include <RooAbsData.h>
0017 #include <RooAbsPdf.h>
0018 #include <RooAddPdf.h>
0019 #include <RooArgList.h>
0020 #include <RooArgSet.h>
0021 #include <RooBinning.h>
0022 #include <RooDataSet.h>
0023 #include <RooFitResult.h>
0024 #include <RooHist.h>
0025 #include <RooMsgService.h>
0026 #include <RooPlot.h>
0027 #include <RooRealVar.h>
0028 #include <RooAbsReal.h>
0029
0030 #include <RooAbsDataHelper.h>
0031
0032 #include <ROOT/RDataFrame.hxx>
0033 #include <ROOT/RDFHelpers.hxx>
0034
0035 #include <TCanvas.h>
0036 #include <TCut.h>
0037 #include <TFile.h>
0038 #include <TH1.h>
0039 #include <TLegend.h>
0040 #include <TLine.h>
0041 #include <TLatex.h>
0042 #include <TPad.h>
0043 #include <TSystem.h>
0044 #include <TTree.h>
0045
0046 #include "./fitutil.h"
0047 #include "./plotutil.h"
0048
0049 using namespace RooFit;
0050 using namespace std;
0051
0052 namespace
0053 {
0054 inline void mkdir_p(const std::string &path)
0055 {
0056     if (path.empty())
0057         return;
0058     gSystem->mkdir(path.c_str(), kTRUE);
0059 }
0060
0061 inline std::string dirname_of(const std::string &p)
0062 {
0063     auto pos = p.find_last_of('/');
0064     if (pos == std::string::npos)
0065         return {};
0066     return p.substr(0, pos);
0067 }
0068
0069 struct SeedGuess
0070 {
0071     double mean = 0.0;
0072     double sigma = 0.01;
0073 };
0074
0075 SeedGuess seed_from_dataset(RooAbsData &data, RooRealVar &mass, int nbins, double window_half_width)
0076 {
0077     SeedGuess out;
0078
0079     std::unique_ptr<TH1> h(data.createHistogram("h_seed", mass, RooFit::Binning(nbins, mass.getMin(), mass.getMax())));
0080     if (!h || h->GetEntries() <= 0)
0081     {
0082         out.mean = 0.5 * (mass.getMin() + mass.getMax());
0083         out.sigma = std::max(1e-4, 0.02 * (mass.getMax() - mass.getMin()));
0084         return out;
0085     }
0086
0087     const int imax = h->GetMaximumBin();
0088     const double peak = h->GetXaxis()->GetBinCenter(imax);
0089
0090     const double loX = std::max(mass.getMin(), peak - window_half_width);
0091     const double hiX = std::min(mass.getMax(), peak + window_half_width);
0092
0093     int lo = std::max(1, h->GetXaxis()->FindBin(loX));
0094     int hi = std::min(h->GetNbinsX(), h->GetXaxis()->FindBin(hiX));
0095
0096     double sw = 0.0, sx = 0.0, sx2 = 0.0;
0097     for (int i = lo; i <= hi; ++i)
0098     {
0099         const double x = h->GetXaxis()->GetBinCenter(i);
0100         const double w = h->GetBinContent(i);
0101         sw += w;
0102         sx += w * x;
0103         sx2 += w * x * x;
0104     }
0105
0106     if (sw > 0.0)
0107     {
0108         out.mean = sx / sw;
0109         const double var = std::max(0.0, sx2 / sw - out.mean * out.mean);
0110         out.sigma = std::sqrt(std::max(var, 1e-12));
0111     }
0112     else
0113     {
0114         out.mean = peak;
0115         out.sigma = std::max(1e-4, window_half_width / 3.0);
0116     }
0117
0118     const double fullW = mass.getMax() - mass.getMin();
0119     out.sigma = std::clamp(out.sigma, 1e-4, 0.25 * fullW);
0120     return out;
0121 }
0122
0123 inline bool fit_is_good(const RooFitResult &fr) { return (fr.status() == 0) && (fr.covQual() >= 2); }
0124
0125 std::unique_ptr<RooFitResult> fit_pdf(RooAbsPdf &pdf, RooAbsData &data, bool extended, int strategy, int printLevel)
0126 {
0127     RooCmdArg ext = extended ? Extended(true) : Extended(false);
0128
0129     RooFitResult *r = pdf.fitTo(data, Save(true), ext, Minimizer("Minuit2", "migrad"), Strategy(strategy), Offset(true), PrintLevel(printLevel));
0130     return std::unique_ptr<RooFitResult>(r);
0131 }
0132
0133 std::string define_sideband_ranges(RooRealVar &mass, double minMass, double maxMass, double mu, double sigma, double gapNSigma)
0134 {
0135     const double sig = std::max(1e-9, std::abs(sigma));
0136     const double gap = std::max(1e-9, gapNSigma * sig);
0137
0138     double lo_hi = std::clamp(mu - gap, minMass, maxMass);
0139     double hi_lo = std::clamp(mu + gap, minMass, maxMass);
0140
0141     bool has_lo = (lo_hi > minMass + 1e-9);
0142     bool has_hi = (maxMass > hi_lo + 1e-9);
0143
0144     std::string ranges;
0145     if (has_lo)
0146     {
0147         mass.setRange("sb_lo", minMass, lo_hi);
0148         ranges += "sb_lo";
0149     }
0150     if (has_hi)
0151     {
0152         mass.setRange("sb_hi", hi_lo, maxMass);
0153         if (!ranges.empty())
0154             ranges += ",";
0155         ranges += "sb_hi";
0156     }
0157     return ranges;
0158 }
0159
0160 void configure_signal_params(const fitparam_config &fit_conf)
0161 {
0162     if (fit_conf.sigmodel == "Gaussian")
0163     {
0164         FitParams::mean.setVal(fit_conf.mean);
0165         FitParams::mean.setRange(fit_conf.mean_low, fit_conf.mean_high);
0166         FitParams::sigma.setVal(fit_conf.sigma);
0167         FitParams::sigma.setRange(fit_conf.sigma_low, fit_conf.sigma_high);
0168     }
0169     else if (fit_conf.sigmodel == "Voigtian")
0170     {
0171         FitParams::mean.setVal(fit_conf.mean);
0172         FitParams::mean.setRange(fit_conf.mean_low, fit_conf.mean_high);
0173         FitParams::sigma.setVal(fit_conf.sigma);
0174         FitParams::sigma.setRange(fit_conf.sigma_low, fit_conf.sigma_high);
0175         FitParams::width.setVal(fit_conf.width);
0176         FitParams::width.setRange(fit_conf.width_low, fit_conf.width_high);
0177     }
0178     else if (fit_conf.sigmodel == "CrystalBall")
0179     {
0180         FitParams::mean.setVal(fit_conf.mean);
0181         FitParams::mean.setRange(fit_conf.mean_low, fit_conf.mean_high);
0182         FitParams::sigma.setVal(fit_conf.sigma);
0183         FitParams::sigma.setRange(fit_conf.sigma_low, fit_conf.sigma_high);
0184         FitParams::alpha1.setVal(fit_conf.alpha1);
0185         FitParams::alpha1.setRange(fit_conf.alpha1_low, fit_conf.alpha1_high);
0186         FitParams::n1.setVal(fit_conf.n1);
0187         FitParams::n1.setRange(fit_conf.n1_low, fit_conf.n1_high);
0188     }
0189     else if (fit_conf.sigmodel == "DoubleCrystalBall")
0190     {
0191         FitParams::mean.setVal(fit_conf.mean);
0192         FitParams::mean.setRange(fit_conf.mean_low, fit_conf.mean_high);
0193         FitParams::sigma.setVal(fit_conf.sigma);
0194         FitParams::sigma.setRange(fit_conf.sigma_low, fit_conf.sigma_high);
0195         FitParams::alpha1.setVal(fit_conf.alpha1);
0196         FitParams::alpha1.setRange(fit_conf.alpha1_low, fit_conf.alpha1_high);
0197         FitParams::n1.setVal(fit_conf.n1);
0198         FitParams::n1.setRange(fit_conf.n1_low, fit_conf.n1_high);
0199         FitParams::alpha2.setVal(fit_conf.alpha2);
0200         FitParams::alpha2.setRange(fit_conf.alpha2_low, fit_conf.alpha2_high);
0201         FitParams::n2.setVal(fit_conf.n2);
0202         FitParams::n2.setRange(fit_conf.n2_low, fit_conf.n2_high);
0203         FitParams::frac.setVal(fit_conf.frac);
0204         FitParams::frac.setRange(0.0, 1.0);
0205     }
0206     else
0207     {
0208         std::cerr << "[RDataframeToRoofit] ERROR: Unknown signal model type: " << fit_conf.sigmodel << std::endl;
0209     }
0210 }
0211
0212 void configure_background_params(const fitparam_config &fit_conf)
0213 {
0214     if (fit_conf.bkgmodel == "Argus")
0215     {
0216         FitParams::k.setVal(fit_conf.k);
0217         FitParams::k.setRange(fit_conf.k_low, fit_conf.k_high);
0218         FitParams::edp.setVal(fit_conf.edp);
0219         FitParams::edp.setRange(fit_conf.edp_low, fit_conf.edp_high);
0220         FitParams::ArgusShift = fit_conf.argusShift;
0221         FitParams::frac_argus.setVal(fit_conf.argusfrac);
0222         FitParams::frac_argus.setRange(0.0, 1.0);
0223         return;
0224     }
0225
0226     const int order = FitParams::ParsePolyOrder(fit_conf.bkgmodel);
0227     if (order < 1)
0228     {
0229         std::cerr << "[RDataframeToRoofit] ERROR: Unknown background model type: " << fit_conf.bkgmodel << std::endl;
0230         return;
0231     }
0232
0233     FitParams::p1.setVal(fit_conf.p1);
0234     FitParams::p1.setRange(fit_conf.p1_low, fit_conf.p1_high);
0235     if (order >= 2)
0236     {
0237         FitParams::p2.setVal(fit_conf.p2);
0238         FitParams::p2.setRange(fit_conf.p2_low, fit_conf.p2_high);
0239     }
0240     if (order >= 3)
0241     {
0242         FitParams::p3.setVal(fit_conf.p3);
0243         FitParams::p3.setRange(fit_conf.p3_low, fit_conf.p3_high);
0244     }
0245     if (order >= 4)
0246     {
0247         FitParams::p4.setVal(fit_conf.p4);
0248         FitParams::p4.setRange(fit_conf.p4_low, fit_conf.p4_high);
0249     }
0250     if (order >= 5)
0251     {
0252         FitParams::p5.setVal(fit_conf.p5);
0253         FitParams::p5.setRange(fit_conf.p5_low, fit_conf.p5_high);
0254     }
0255     if (order >= 6)
0256     {
0257         FitParams::p6.setVal(fit_conf.p6);
0258         FitParams::p6.setRange(fit_conf.p6_low, fit_conf.p6_high);
0259     }
0260     if (order >= 7)
0261     {
0262         FitParams::p7.setVal(fit_conf.p7);
0263         FitParams::p7.setRange(fit_conf.p7_low, fit_conf.p7_high);
0264     }
0265     if (order >= 8)
0266     {
0267         FitParams::p8.setVal(fit_conf.p8);
0268         FitParams::p8.setRange(fit_conf.p8_low, fit_conf.p8_high);
0269     }
0270     if (order >= 9)
0271     {
0272         FitParams::p9.setVal(fit_conf.p9);
0273         FitParams::p9.setRange(fit_conf.p9_low, fit_conf.p9_high);
0274     }
0275     if (order >= 10)
0276     {
0277         FitParams::p10.setVal(fit_conf.p10);
0278         FitParams::p10.setRange(fit_conf.p10_low, fit_conf.p10_high);
0279     }
0280 }
0281
0282 } // namespace
0283
0284 double RDataframeToRoofit(const bool doSnapshot, std::string snapshotName, std::string inputfilename, TCut selections, fitparam_config &fit_conf, const std::string plotdir = "./figure", const std::vector<std::string> &cut_legend_entries = {})
0285 {
0286     RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
0287
0288     mkdir_p(plotdir);
0289     if (doSnapshot)
0290         mkdir_p(dirname_of(snapshotName));
0291
0292     auto start = std::chrono::high_resolution_clock::now();
0293
0294     static bool mt_enabled = false;
0295     if (!mt_enabled)
0296     {
0297         ROOT::EnableImplicitMT();
0298         mt_enabled = true;
0299     }
0300
0301     ROOT::RDataFrame df("DecayTree", inputfilename.c_str());
0302     auto filtered_df = df.Filter(selections.GetTitle());
0303
0304     FitParams::branch = fit_conf.branch;
0305     FitParams::minMass = fit_conf.minMass;
0306     FitParams::maxMass = fit_conf.maxMass;
0307
0308     FitParams::mass.SetName(FitParams::branch.c_str());
0309     FitParams::mass.SetTitle("mass");
0310     FitParams::mass.setRange(FitParams::minMass, FitParams::maxMass);
0311
0312     auto dsAction = filtered_df.Book<float>(RooDataSetHelper("dataset", "dataset", RooArgSet(FitParams::mass)), {FitParams::branch.c_str()});
0313
0314     std::vector<ROOT::RDF::RResultHandle> handles;
0315     handles.emplace_back(dsAction);
0316
0317     if (doSnapshot)
0318     {
0319         auto snapAction = filtered_df.Snapshot("DecayTree", snapshotName.c_str());
0320         handles.emplace_back(snapAction);
0321     }
0322
0323     ROOT::RDF::RunGraphs(handles);
0324
0325     RooDataSet *dataset = dsAction.GetPtr();
0326     const int nEntries = dataset ? dataset->numEntries() : 0;
0327     std::cout << "Number of entries in the dataset: " << nEntries << std::endl;
0328
0329     if (!dataset || nEntries <= 0)
0330     {
0331         std::cerr << "[RDataframeToRoofit] WARNING: Empty dataset after selections.\n";
0332         return 0.0;
0333     }
0334
0335     FitParams::BeginFit();
0336
0337     configure_signal_params(fit_conf);
0338     configure_background_params(fit_conf);
0339
0340     RooAbsPdf *signal = signalModel(fit_conf.sigmodel, FitParams::mass);
0341     RooAbsPdf *background = backgroundModel(fit_conf.bkgmodel, FitParams::mass);
0342     if (!signal || !background)
0343     {
0344         std::cerr << "[RDataframeToRoofit] ERROR: Failed to build signal/background PDFs.\n";
0345         return 0.0;
0346     }
0347
0348     const double nTot = static_cast<double>(nEntries);
0349
0350     const double nSigStart = (fit_conf.nSig > 0.0) ? fit_conf.nSig : fit_conf.nSig_fracTol * nTot;
0351     FitParams::nSig.setVal(std::clamp(nSigStart, 0.0, nTot));
0352     if (fit_conf.nSig_low > 0.0)
0353     {
0354         FitParams::nSig.setRange(fit_conf.nSig_low, nTot);
0355     }
0356
0357     const double nBkgStart = (fit_conf.nBkg > 0.0) ? fit_conf.nBkg : fit_conf.nBkg_fracTol * nTot;
0358     FitParams::nBkg.setVal(std::clamp(nBkgStart, 0.0, nTot));
0359     if (fit_conf.nBkg_low > 0.0)
0360     {
0361         FitParams::nBkg.setRange(fit_conf.nBkg_low, nTot);
0362     }
0363
0364     RooAddPdf model("model", "model", RooArgList(*signal, *background), RooArgList(FitParams::nSig, FitParams::nBkg));
0365
0366     const double fullW = FitParams::maxMass - FitParams::minMass;
0367     SeedGuess seed0 = {FitParams::mean.getVal(), FitParams::sigma.getVal()};
0368     if (!(seed0.sigma > 0.0) || seed0.sigma > 0.3 * fullW)
0369     {
0370         const double win = std::min(0.06, 0.25 * fullW);
0371         seed0 = seed_from_dataset(*dataset, FitParams::mass, std::max(40, fit_conf.nBins), win);
0372         FitParams::mean.setVal(seed0.mean);
0373         FitParams::sigma.setVal(seed0.sigma);
0374     }
0375
0376     std::unique_ptr<RooFitResult> bkg_prefit;
0377     if (fit_conf.useStagedFit)
0378     {
0379         const std::string sbRanges = define_sideband_ranges(FitParams::mass, FitParams::minMass, FitParams::maxMass, seed0.mean, seed0.sigma, fit_conf.stagedGapNSigma);
0380
0381         if (!sbRanges.empty())
0382         {
0383             std::unique_ptr<RooAbsData> sbData(dataset->reduce(CutRange(sbRanges.c_str())));
0384             if (sbData && sbData->numEntries() > 0)
0385             {
0386                 bkg_prefit = fit_pdf(*background, *sbData, false, 1, -1);
0387                 std::cout << "[RDataframeToRoofit] Sideband bkg prefit: status=" << bkg_prefit->status() << ", covQual=" << bkg_prefit->covQual() << "\n";
0388             }
0389             else
0390             {
0391                 std::cout << "[RDataframeToRoofit] Sideband bkg prefit skipped (empty sideband data).\n";
0392             }
0393         }
0394         else
0395         {
0396             std::cout << "[RDataframeToRoofit] Sideband bkg prefit skipped (sidebands not usable).\n";
0397         }
0398     }
0399
0400     std::unique_ptr<RooFitResult> fitres = fit_pdf(model, *dataset, true, 1, -1);
0401
0402     if (!fit_is_good(*fitres))
0403     {
0404         std::cerr << "[RDataframeToRoofit] WARNING: Full fit not good (status=" << fitres->status() << ", covQual=" << fitres->covQual() << "). Reseed + retry.\n";
0405
0406         const double win = std::min(0.06, 0.25 * fullW);
0407         SeedGuess seed = seed_from_dataset(*dataset, FitParams::mass, std::max(40, fit_conf.nBins), win);
0408
0409         FitParams::mean.setVal(seed.mean);
0410         FitParams::mean.setRange(std::max(FitParams::minMass, seed.mean - 0.5 * win), std::min(FitParams::maxMass, seed.mean + 0.5 * win));
0411
0412         FitParams::sigma.setVal(seed.sigma);
0413         FitParams::sigma.setRange(fit_conf.sigma_low, fit_conf.sigma_high);
0414
0415         if (fit_conf.useStagedFit)
0416         {
0417             const std::string sbRanges = define_sideband_ranges(FitParams::mass, FitParams::minMass, FitParams::maxMass, seed.mean, seed.sigma, fit_conf.stagedGapNSigma);
0418
0419             if (!sbRanges.empty())
0420             {
0421                 std::unique_ptr<RooAbsData> sbData(dataset->reduce(CutRange(sbRanges.c_str())));
0422                 if (sbData && sbData->numEntries() > 0)
0423                     bkg_prefit = fit_pdf(*background, *sbData, false, 2, -1);
0424             }
0425         }
0426
0427         std::unique_ptr<RooFitResult> fitres2 = fit_pdf(model, *dataset, true, 2, -1);
0428
0429         const bool good1 = fit_is_good(*fitres);
0430         const bool good2 = fit_is_good(*fitres2);
0431
0432         if (good2 && !good1)
0433             fitres.swap(fitres2);
0434         else if (!good1 && !good2 && fitres2->minNll() < fitres->minNll())
0435             fitres.swap(fitres2);
0436     }
0437
0438     fitres->Print();
0439
0440     {
0441         TFile out(Form("%s/fitresults.root", plotdir.c_str()), "RECREATE");
0442         out.cd();
0443         fitres->Write("fitres");
0444         if (bkg_prefit)
0445             bkg_prefit->Write("bkg_prefit_sidebands");
0446     }
0447
0448     const double mu = FitParams::mean.getVal();
0449     const double sg = FitParams::sigma.getVal();
0450
0451     double lo = std::max(FitParams::minMass, mu - 3.0 * sg);
0452     double hi = std::min(FitParams::maxMass, mu + 3.0 * sg);
0453     if (!(hi > lo))
0454     {
0455         lo = FitParams::minMass;
0456         hi = FitParams::maxMass;
0457     }
0458
0459     FitParams::mass.setRange("integralrange", lo, hi);
0460
0461     std::unique_ptr<RooAbsReal> signalFrac(signal->createIntegral(FitParams::mass, NormSet(FitParams::mass), Range("integralrange")));
0462     std::unique_ptr<RooAbsReal> bkgFrac(background->createIntegral(FitParams::mass, NormSet(FitParams::mass), Range("integralrange")));
0463
0464     const double S = signalFrac ? signalFrac->getVal() * FitParams::nSig.getVal() : 0.0;
0465     const double B = bkgFrac ? bkgFrac->getVal() * FitParams::nBkg.getVal() : 0.0;
0466     const double significance = (S + B > 0.0) ? (S / std::sqrt(S + B)) : 0.0;
0467
0468     std::string xAxisTitle = fit_conf.decaystring + " candidate mass [GeV]";
0469
0470     RooBinning bins(FitParams::minMass, FitParams::maxMass);
0471     bins.addUniform(fit_conf.nBins, FitParams::minMass, FitParams::maxMass);
0472
0473     RooPlot *frame_display = FitParams::mass.frame(Title(""));
0474     RooPlot *frame_calc = FitParams::mass.frame(Title(""));
0475
0476     dataset->plotOn(frame_calc, Binning(bins), DataError(RooAbsData::SumW2), Name("data_calc"));
0477     model.plotOn(frame_calc, LineColor(kBlack), Name("total_fit_calc"));
0478     const double chi2ndf = frame_calc->chiSquare("total_fit_calc", "data_calc");
0479
0480     RooHist *pull = frame_calc->pullHist("data_calc", "total_fit_calc");
0481     RooPlot *frame_pull = FitParams::mass.frame(Title(""));
0482     frame_pull->addPlotable(pull, "PE1");
0483
0484     dataset->plotOn(frame_display, Binning(bins), XErrorSize(0), DataError(RooAbsData::SumW2), Name("data_binned"));
0485     model.plotOn(frame_display, Components(*background), LineColor(kGray), DrawOption("F"), FillColor(kGray), Name("bkg"));
0486     model.plotOn(frame_display, Components(RooArgSet(*signal, *background)), LineColor(kAzure + 8), DrawOption("F"), FillColor(kAzure + 8), MoveToBack(), Name("sig_plus_bkg"));
0487     model.plotOn(frame_display, LineColor(kBlack), Name("total_fit"));
0488     dataset->plotOn(frame_display, DrawOption("PE1"), Binning(bins), XErrorSize(0), DataError(RooAbsData::SumW2), Name("data"));
0489
0490     TCanvas *c = new TCanvas("massFitCanvas", "massFitCanvas", 800, 800);
0491
0492     TPad mainPad("mainPad", "mainPad", 0., 0.3, 1., 1.);
0493     mainPad.SetTopMargin(TopMargin);
0494     mainPad.SetBottomMargin(0);
0495     mainPad.Draw();
0496
0497     TPad pullPad("pullPad", "pullPad", 0., 0.0, 1., 0.3);
0498     pullPad.SetBottomMargin(0.5);
0499     pullPad.SetTopMargin(0);
0500     pullPad.Draw();
0501
0502     mainPad.cd();
0503     frame_display->SetMarkerStyle(kCircle);
0504     frame_display->SetMarkerSize(0.02);
0505     frame_display->SetLineWidth(1);
0506     frame_display->GetXaxis()->SetTitleSize(0);
0507     frame_display->GetXaxis()->SetLabelSize(0);
0508
0509     frame_display->GetYaxis()->SetTitleSize(AxisTitleSize * textscale_pad1);
0510     frame_display->GetYaxis()->SetLabelSize(AxisLabelSize * textscale_pad1);
0511     frame_display->GetYaxis()->SetTitleOffset(1.2);
0512     frame_display->GetYaxis()->SetTitleFont(42);
0513     frame_display->GetYaxis()->SetLabelFont(42);
0514
0515     std::unique_ptr<TH1> hdataset(dataset->createHistogram("hdataset", FitParams::mass, Binning(fit_conf.nBins, FitParams::minMass, FitParams::maxMass)));
0516     frame_display->GetYaxis()->SetRangeUser(0.1, hdataset ? hdataset->GetMaximum() * 1.8 : 1.0);
0517
0518     float binWidth = 1000.f * float(FitParams::maxMass - FitParams::minMass) / float(fit_conf.nBins);
0519     string yAxisTitle = "Candidates / (" + to_string_with_precision(binWidth, 1) + " MeV)";
0520     frame_display->GetYaxis()->SetTitle(yAxisTitle.c_str());
0521     frame_display->Draw();
0522     c->RedrawAxis();
0523
0524     TLatex *datestamp = new TLatex();
0525     datestamp->SetTextSize(0.06);
0526     datestamp->SetTextAlign(kHAlignRight + kVAlignBottom);
0527     datestamp->SetNDC();
0528     datestamp->DrawLatex(1 - mainPad.GetRightMargin(), 1 - mainPad.GetTopMargin() + 0.01, getTime().c_str());
0529
0530     TLegend *sphnxleg = new TLegend(mainPad.GetLeftMargin() + 0.03, 1 - mainPad.GetTopMargin() - 0.2, mainPad.GetLeftMargin() + 0.2, 1 - mainPad.GetTopMargin() - 0.05);
0531     sphnxleg->SetTextAlign(kHAlignLeft + kVAlignCenter);
0532     sphnxleg->SetTextSize(0.06);
0533     sphnxleg->SetFillStyle(0);
0534     sphnxleg->AddEntry("", Form("#it{#bf{sPHENIX}} %s", prelimtext.c_str()), "");
0535     sphnxleg->AddEntry("", "p+p #sqrt{s_{NN}}=200 GeV", "");
0536     sphnxleg->Draw();
0537
0538     TLegend *leg = new TLegend(1 - mainPad.GetRightMargin() - 0.33, 1 - mainPad.GetTopMargin() - 0.33, 1 - mainPad.GetRightMargin() - 0.1, 1 - mainPad.GetTopMargin() - 0.06);
0539     leg->AddEntry(frame_display->findObject("data"), "Data", "PE2");
0540     leg->AddEntry(frame_display->findObject("total_fit"), "Fit", "L");
0541     leg->AddEntry(frame_display->findObject("sig_plus_bkg"), fit_conf.decaystring.c_str(), "f");
0542     leg->AddEntry(frame_display->findObject("bkg"), "Comb. Bkg.", "f");
0543     leg->SetFillColor(0);
0544     leg->SetFillStyle(0);
0545     leg->SetBorderSize(0);
0546     leg->SetTextSize(0.055);
0547     leg->Draw();
0548
0549     TLatex *text = new TLatex();
0550     text->SetTextSize(0.045);
0551     text->SetTextAlign(kHAlignLeft + kVAlignCenter);
0552     text->SetNDC();
0553
0554     double y0 = 1 - mainPad.GetTopMargin() - 0.25;
0555     text->DrawLatex(mainPad.GetLeftMargin() + 0.07, y0, Form("S/#sqrt{S+B} = %.2f", significance));
0556     text->DrawLatex(mainPad.GetLeftMargin() + 0.07, y0 - 0.06, Form("#mu = %.0f #pm %.0f MeV", FitParams::mean.getVal() * 1000.0, FitParams::mean.getError() * 1000.0));
0557     text->DrawLatex(mainPad.GetLeftMargin() + 0.07, y0 - 0.115, Form("#sigma = %.2f #pm %.2f MeV", FitParams::sigma.getVal() * 1000.0, FitParams::sigma.getError() * 1000.0));
0558     text->DrawLatex(mainPad.GetLeftMargin() + 0.07, y0 - 0.170, Form("Yield = %.0f #pm %.0f", FitParams::nSig.getVal(), FitParams::nSig.getError()));
0559     // text->DrawLatex(mainPad.GetLeftMargin() + 0.07, y0 - 0.225, Form("#chi^{2}/ndf #approx %.2f", chi2ndf));
0560
0561     if (!cut_legend_entries.empty())
0562     {
0563         const double cutTop = 1 - mainPad.GetTopMargin() - 0.35;
0564         const double cutBottom = cutTop - 0.045*cut_legend_entries.size();
0565         TLegend *cutleg = new TLegend(1 - mainPad.GetRightMargin() - 0.33, cutBottom, 1 - mainPad.GetRightMargin() - 0.1, cutTop);
0566         cutleg->SetFillColor(0);
0567         cutleg->SetFillStyle(0);
0568         cutleg->SetBorderSize(0);
0569         cutleg->SetTextSize(0.04);
0570         cutleg->SetTextAlign(kHAlignLeft + kVAlignCenter);
0571         for (const auto &entry : cut_legend_entries)
0572             cutleg->AddEntry("", entry.c_str(), "");
0573         cutleg->Draw();
0574     }
0575
0576     // TLatex *qual = new TLatex();
0577     // qual->SetTextSize(0.04);
0578     // qual->SetTextAlign(kHAlignRight + kVAlignTop);
0579     // qual->SetNDC();
0580     // qual->DrawLatex(1 - mainPad.GetRightMargin(), 1 - mainPad.GetTopMargin() - 0.02, Form("status=%d, covQual=%d", fitres->status(), fitres->covQual()));
0581
0582     pullPad.cd();
0583     frame_pull->SetMarkerStyle(kCircle);
0584     frame_pull->SetMarkerSize(0.02);
0585     frame_pull->SetTitle("");
0586     frame_pull->GetXaxis()->SetTitle(xAxisTitle.c_str());
0587     frame_pull->GetXaxis()->SetTitleOffset(1.3);
0588     frame_pull->GetXaxis()->SetTitleFont(42);
0589     frame_pull->GetXaxis()->SetTitleSize(AxisTitleSize * textscale_pad2);
0590     frame_pull->GetXaxis()->SetLabelFont(42);
0591     frame_pull->GetXaxis()->SetLabelSize(AxisLabelSize * textscale_pad2);
0592
0593     frame_pull->GetYaxis()->SetTitle("Pull");
0594     frame_pull->GetYaxis()->SetTitleOffset(0.5);
0595     frame_pull->GetYaxis()->SetTitleFont(42);
0596     frame_pull->GetYaxis()->SetTitleSize(AxisTitleSize * textscale_pad2);
0597     frame_pull->GetYaxis()->SetLabelFont(42);
0598     frame_pull->GetYaxis()->SetLabelSize(AxisLabelSize * textscale_pad2);
0599     frame_pull->GetYaxis()->SetRangeUser(-6, 6);
0600     frame_pull->GetYaxis()->SetNdivisions(5);
0601
0602     frame_pull->Draw();
0603
0604     TLine *plusThreeLine = new TLine(FitParams::minMass, 3, FitParams::maxMass, 3);
0605     plusThreeLine->SetLineColor(1);
0606     plusThreeLine->SetLineStyle(2);
0607     plusThreeLine->SetLineWidth(2);
0608     plusThreeLine->Draw("same");
0609     TLine *zeroLine = new TLine(FitParams::minMass, 0, FitParams::maxMass, 0);
0610     zeroLine->SetLineColor(1);
0611     zeroLine->SetLineStyle(2);
0612     zeroLine->SetLineWidth(2);
0613     zeroLine->Draw("same");
0614     TLine *minusThreeLine = new TLine(FitParams::minMass, -3, FitParams::maxMass, -3);
0615     minusThreeLine->SetLineColor(1);
0616     minusThreeLine->SetLineStyle(2);
0617     minusThreeLine->SetLineWidth(2);
0618     minusThreeLine->Draw("same");
0619
0620     const std::vector<std::string> exts = {".C", ".pdf", ".png"};
0621     for (const auto &ext : exts)
0622         c->SaveAs(Form("%s/RDataframeToRoofit_%s%s", plotdir.c_str(), FitParams::branch.c_str(), ext.c_str()));
0623
0624     auto end = std::chrono::high_resolution_clock::now();
0625     std::chrono::duration<double> elapsed = end - start;
0626     std::cout << "Elapsed time: " << elapsed.count() / 60.0 << " minutes\n";
0627
0628     return significance;
0629 }
0630
0631 #endif // RDATAFRAMETOROOFIT_C