pfeiffer
/
inclusive_detached_dilepton
1
0
Fork 0
Code Issues Pull Requests 1 Projects Releases Wiki Activity
ROOT Analysis for the Inclusive Detachted Dilepton Trigger Lines
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
27 Commits
2 Branches
0 Tags
393 KiB
 
 
Branch: master
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from 'master'
${ noResults }
inclusive_detached_dilepton/old_analysis_macros/status_report_plots.cpp

484 lines
24 KiB
Raw Permalink Blame History

#include "TH1D.h"
#include "TH2D.h"
#include "THStack.h"
#include "TGraph.h"
#include "TTree.h"
#include "TF1.h"
#include "TChain.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TColor.h"
#include "TLorentzVector.h"
#include "TRandom3.h"
#include "TLorentzVector.h"
#include "RooDataHist.h"
#include "RooRealVar.h"
#include "RooPlot.h"
#include "RooGaussian.h"
#include "RooExponential.h"
#include "RooRealConstant.h"
#include "RooAddPdf.h"
#include "RooFitResult.h"
#include "RooProduct.h"
#include "RooCrystalBall.h"
#include "RooBreitWigner.h"
#include "RooArgSet.h"
#include "RooFFTConvPdf.h"
#include "RooNovosibirsk.h"
#include <string>
#include <iostream>
#include <cmath>
const int nBins = 70;
// const Double_t MASS_HIST_MIN = 5150.;
// const Double_t MASS_HIST_MAX = 5450.;
// const Double_t MASS_HIST_FIT_MIN = 5150.;
// const Double_t MASS_HIST_FIT_MAX = 5450.;
const Double_t MASS_HIST_MIN = 4000.;
const Double_t MASS_HIST_MAX = 8500.;
const Double_t MASS_HIST_FIT_MIN = 5100.;
const Double_t MASS_HIST_FIT_MAX = 6000.;
// PDG Values
const Double_t J_PSI_MASS = 3096.916;
const Double_t PSI_2S_MASS = 3686.09;
const Double_t K_MASS = 493.677;
const Double_t K_STAR_0_MASS = 891.67;
const std::string SAVE_PATH = "/work/pfeiffer/inclusive_detached_dilepton/status_report";
struct FitParams
{
Double_t lambda;
Double_t mean;
Double_t sigma;
Double_t sig_yield;
Double_t bkg_yield;
FitParams(Double_t lambda, Double_t mean, Double_t sigma, Double_t sig_yield, Double_t bkg_yield)
{
this->lambda = lambda;
this->mean = mean;
this->sigma = sigma;
this->sig_yield = sig_yield;
this->bkg_yield = bkg_yield;
}
};
struct AnalysisOutput
{
std::string title;
std::string name;
std::string root_file;
std::string root_file_tree;
std::string B_name;
int k_charge;
const std::string suf(std::string input)
{
return TString::Format("%s_%s", input.c_str(), name.c_str()).Data();
}
const std::string title_suf(std::string input)
{
return TString::Format("%s [%s]", input.c_str(), title.c_str()).Data();
}
};
void CreateRooFitAndSavePDF(TH1D *hist, AnalysisOutput ana, const char *name);
std::vector<RooPlot *> PlotWithParams(TH1D *hist, AnalysisOutput ana);
bool inRange(double value, double center, double low_intvl, double up_intvl)
{
return center - low_intvl < value && value < center + up_intvl;
}
bool inRange(double value, double center, double intvl)
{
return inRange(value, center, intvl, intvl);
}
void savePDF(TObject *o, const char *decay, const char *name, Option_t *opt = "")
{
auto cname = TString::Format("%s_%s", decay, name);
std::cout << " ----- " << cname.Data() << " - " << decay << " - " << name << std::endl;
auto c = new TCanvas(cname.Data(), cname.Data(), 0, 0, 800, 600);
o->Draw(opt);
c->SaveAs(TString::Format("%s/%s.jpg", SAVE_PATH.c_str(), cname.Data()));
// c->Draw();
}
RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING);
int status_report_plots()
{
// gROOT->ProcessLine(".x /work/pfeiffer/lhcbStyle.C");
std::vector<AnalysisOutput> anas{
AnalysisOutput{"SpruceRD_BuToHpMuMu", "BuToKpMuMu_incl", "/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_Sprucing23r1_90000000_RD.root", "SpruceRD_BuToHpMuMu/DecayTree", "B", 0},
// AnalysisOutput{"SpruceRD_BuToHpMuMu (ap)", "BuToKpMuMu_incl_ap", "/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0r1_tuple_90000000_2023_v0r0p6288631.root", "BuToHpMuMu/DecayTree", "B", "1395.27 +/- 58"},
// AnalysisOutput{"Hlt2RD_BuToKpMuMu_2023", "BuToKpMuMu_excl", "/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_RD.root", "Hlt2RD_BuToKpMuMu_2023/DecayTree", "B", 0 },
// // AnalysisOutput{"Hlt2RD_", "BuToKpMuMu_excl_ap23", "/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/turbopass_magdown_2023_v1_tuple_94000000_2023_v0r0p6201764.root", "BuToKpMuMu23/DecayTree", "B", "803.769 +/- 34"},
// AnalysisOutput{"SpruceRD_B0ToHpHmMuMu K+", "B0ToKpPimMuMu_Kp_incl", "/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_Sprucing23r1_90000000_RD.root", "SpruceRD_B0ToHpHmMuMu/DecayTree", "B0", 1 },
// AnalysisOutput{"SpruceRD_B0ToHpHmMuMu K-", "B0ToKpPimMuMu_Km_incl", "/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_Sprucing23r1_90000000_RD.root", "SpruceRD_B0ToHpHmMuMu/DecayTree", "B0", -1 },
// AnalysisOutput{"SpruceRD_B0ToHpHmMuMu (ap)", "B0ToKpPimMuMu_incl_ap", "/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0r1_tuple_90000000_2023_v0r0p6288631.root", "B0ToHpHmMuMu/DecayTree", "B0", ""},
// AnalysisOutput{"Hlt2RD_B0ToKpPimMuMu_2023", "B0ToKpPimMuMu_excl", "/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_RD.root", "Hlt2RD_B0ToKpPimMuMu_2023/DecayTree", "B0", 0 },
// AnalysisOutput{"Hlt2RD_B0ToKpPimMuMu_2023 AP", "B0ToKpPimMuMu_excl_ap23", "/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/turbopass_magdown_2023_v1_tuple_94000000_2023_v0r0p6201764.root", "B0ToKpPimMuMu23/DecayTree", "B0", ""},
// AnalysisOutput{"Hlt2B2CC_BuToJpsiKplus_JpsiToMuMu_Detached", "BuToJpsiKplus_det", "/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_B2CC.root", "Hlt2B2CC_BuToJpsiKplus_JpsiToMuMu_Detached/DecayTree", "B", 0 },
};
for (size_t a = 0; a < anas.size(); a++)
{
auto ana = anas[a];
TChain *data_chain = new TChain(ana.root_file_tree.c_str());
data_chain->Add(ana.root_file.c_str());
Double_t B_M, // 4400 -> 8200
Jpsi_M, // 200 -> 6600
muminus_PID_MU,
muplus_PID_MU,
muminus_PID_K,
muplus_PID_K,
Kplus_PID_K;
Float_t muminus_PT,
muplus_PT,
Kplus_PT,
Kplus_P;
Int_t Kplus_Q;
Bool_t muplus_ISMUON, muminus_ISMUON, Hlt1TrackMVADecision, Hlt1TwoTrackMVADecision;
data_chain->SetBranchAddress(TString::Format("%s_M", ana.B_name.c_str()).Data(), &B_M);
data_chain->SetBranchAddress("Jpsi_M", &Jpsi_M);
data_chain->SetBranchAddress("muplus_ISMUON", &muplus_ISMUON);
data_chain->SetBranchAddress("muminus_ISMUON", &muminus_ISMUON);
data_chain->SetBranchAddress("Hlt1TrackMVADecision", &Hlt1TrackMVADecision);
data_chain->SetBranchAddress("Hlt1TwoTrackMVADecision", &Hlt1TwoTrackMVADecision);
data_chain->SetBranchAddress("Kplus_Q", &Kplus_Q);
// manually sub mass hyp
Double_t L1_M, L2_M, Hp_M;
Float_t L1_PX, L1_PY, L1_PZ, L1_ENERGY, L2_PX, L2_PY, L2_PZ, L2_ENERGY, Hp_PX, Hp_PY, Hp_PZ, Hp_ENERGY;
if (ana.name == "BuToKpMuMu_incl_ap" || ana.name == "B0ToKpPimMuMu_incl_ap")
{
data_chain->SetBranchAddress("L1_M", &L1_M);
data_chain->SetBranchAddress("L2_M", &L2_M);
data_chain->SetBranchAddress("Hp_M", &Hp_M);
data_chain->SetBranchAddress("L1_PX", &L1_PX);
data_chain->SetBranchAddress("L1_PY", &L1_PY);
data_chain->SetBranchAddress("L1_PZ", &L1_PZ);
data_chain->SetBranchAddress("L1_ENERGY", &L1_ENERGY);
data_chain->SetBranchAddress("L2_PX", &L2_PX);
data_chain->SetBranchAddress("L2_PY", &L2_PY);
data_chain->SetBranchAddress("L2_PZ", &L2_PZ);
data_chain->SetBranchAddress("L2_ENERGY", &L2_ENERGY);
data_chain->SetBranchAddress("Hp_PX", &Hp_PX);
data_chain->SetBranchAddress("Hp_PY", &Hp_PY);
data_chain->SetBranchAddress("Hp_PZ", &Hp_PZ);
data_chain->SetBranchAddress("Hp_ENERGY", &Hp_ENERGY);
}
Float_t muplus_PX, muplus_PY, muplus_PZ, muplus_ENERGY, muminus_PX, muminus_PY, muminus_PZ, muminus_ENERGY;
data_chain->SetBranchAddress("muplus_PX", &muplus_PX);
data_chain->SetBranchAddress("muplus_PY", &muplus_PY);
data_chain->SetBranchAddress("muplus_PZ", &muplus_PZ);
data_chain->SetBranchAddress("muplus_ENERGY", &muplus_ENERGY);
data_chain->SetBranchAddress("muminus_PX", &muminus_PX);
data_chain->SetBranchAddress("muminus_PY", &muminus_PY);
data_chain->SetBranchAddress("muminus_PZ", &muminus_PZ);
data_chain->SetBranchAddress("muminus_ENERGY", &muminus_ENERGY);
Double_t Kst0_M, Kplus_M, piminus_M;
Float_t Kst0_PX, Kst0_PY, Kst0_PZ, Kst0_ENERGY, Kplus_PX, Kplus_PY, Kplus_PZ, Kplus_ENERGY, piminus_PX, piminus_PY, piminus_PZ, piminus_ENERGY;
if (ana.B_name == "B0")
{
data_chain->SetBranchAddress("Kst0_M", &Kst0_M);
data_chain->SetBranchAddress("Kst0_PX", &Kst0_PX);
data_chain->SetBranchAddress("Kst0_PY", &Kst0_PY);
data_chain->SetBranchAddress("Kst0_PZ", &Kst0_PZ);
data_chain->SetBranchAddress("Kst0_ENERGY", &Kst0_ENERGY);
data_chain->SetBranchAddress("Kplus_M", &Kplus_M);
data_chain->SetBranchAddress("Kplus_PX", &Kplus_PX);
data_chain->SetBranchAddress("Kplus_PY", &Kplus_PY);
data_chain->SetBranchAddress("Kplus_PZ", &Kplus_PZ);
data_chain->SetBranchAddress("Kplus_ENERGY", &Kplus_ENERGY);
data_chain->SetBranchAddress("piminus_M", &piminus_M);
data_chain->SetBranchAddress("piminus_PX", &piminus_PX);
data_chain->SetBranchAddress("piminus_PY", &piminus_PY);
data_chain->SetBranchAddress("piminus_PZ", &piminus_PZ);
data_chain->SetBranchAddress("piminus_ENERGY", &piminus_ENERGY);
}
TH1D *h1_B_M = new TH1D(ana.suf("h1_B_M").c_str(), ana.title_suf("B Mass").c_str(), nBins, MASS_HIST_MIN, MASS_HIST_MAX);
TH1D *h1_B_M_JPsi_cut = new TH1D(ana.suf("h1_B_M_JPsi_cut").c_str(), ana.title_suf("B Mass").c_str(), nBins, MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX);
TH1D *h1_B_M_Psi2s_cut = new TH1D(ana.suf("h1_B_M_Psi2s_cut").c_str(), ana.title_suf("B Mass").c_str(), nBins, MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX);
TH1D *h1_Jpsi_M = new TH1D(ana.suf("h1_Jpsi_M").c_str(), ana.title_suf("J/#psi Mass").c_str(), nBins, 200., 5500.);
TH1D *h1_Jpsi_M_align = new TH1D(ana.suf("h1_Jpsi_M_align").c_str(), ana.title_suf("J/#psi Mass Align").c_str(), nBins, J_PSI_MASS - 200., J_PSI_MASS + 200.);
TH2D *h2_B_M_Jpsi_M = new TH2D(ana.suf("h2_B_M_Jpsi_M").c_str(), ana.title_suf("B Mass vs. J/#psi Mass").c_str(), nBins, MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX, nBins, 200., 5500.);
TH1D *h1_Kst0_M = new TH1D(ana.suf("h1_Kst0_M").c_str(), ana.title_suf("K*0 Mass").c_str(), nBins, 600., 3000.);
TH1D *h1_Kpi_M = new TH1D(ana.suf("h1_Kpi_M").c_str(), ana.title_suf("K #pi Mass").c_str(), nBins, 600., 3000.);
h1_B_M->GetXaxis()->SetTitle(TString::Format("m(%s) / MeV", ana.B_name.c_str()).Data());
h1_B_M_JPsi_cut->GetXaxis()->SetTitle(TString::Format("m(%s) / MeV", ana.B_name.c_str()).Data());
h1_Jpsi_M->GetXaxis()->SetTitle("m(#mu#mu) / MeV");
unsigned int entries = data_chain->GetEntries();
for (unsigned int i = 0; i < entries; i++)
{
if (i % 10000 == 0)
{
std::cout << "[" << ana.name << "] Processing event: " << i << " (" << TString::Format("%.2f", ((double)i / (double)entries) * 100.) << "%)" << std::endl;
}
data_chain->GetEntry(i);
if (ana.k_charge != 0 && Kplus_Q != ana.k_charge)
{
continue;
}
// if (!(muplus_ISMUON && muminus_ISMUON && (Hlt1TrackMVADecision | Hlt1TwoTrackMVADecision)))
// {
// continue;
// }
Double_t used_B_Mass = 0;
// manually sub mass hyp
if (ana.name == "BuToKpMuMu_incl_ap")
{
TVector3 K_momentum(Hp_PX, Hp_PY, Hp_PZ);
double K_energy = TMath::Sqrt(TMath::Sq(K_MASS) + K_momentum.Mag2());
TLorentzVector K_4v(K_momentum, K_energy);
TLorentzVector l1_4v(L1_PX, L1_PY, L1_PZ, L1_ENERGY);
TLorentzVector l2_4v(L2_PX, L2_PY, L2_PZ, L2_ENERGY);
used_B_Mass = (K_4v + l1_4v + l2_4v).M();
}
else
{
used_B_Mass = B_M;
}
TLorentzVector muplus_4v(muplus_PX, muplus_PY, muplus_PZ, muplus_ENERGY);
TLorentzVector muminus_4v(muminus_PX, muminus_PY, muminus_PZ, muminus_ENERGY);
Double_t calc_q2 =
h1_B_M->Fill(used_B_Mass);
h1_Jpsi_M->Fill(Jpsi_M);
h1_Jpsi_M_align->Fill(Jpsi_M);
h2_B_M_Jpsi_M->Fill(used_B_Mass, Jpsi_M);
if (ana.B_name == "B0")
{
TLorentzVector kplus_4v(Kplus_PX, Kplus_PY, Kplus_PZ, Kplus_ENERGY);
TLorentzVector piminus_4v(piminus_PX, piminus_PY, piminus_PZ, piminus_ENERGY);
TLorentzVector kst0_4v(Kst0_PX, Kst0_PY, Kst0_PZ, Kst0_ENERGY);
Double_t kpi_mass = (kplus_4v + piminus_4v).M();
Double_t kst0_mass = kst0_4v.M();
if (TMath::Abs(kpi_mass - K_STAR_0_MASS) < 100.)
{
h1_Kst0_M->Fill(kst0_mass);
h1_Kpi_M->Fill(kpi_mass);
if (TMath::Abs(Jpsi_M - J_PSI_MASS) < 60.)
{
h1_B_M_JPsi_cut->Fill(used_B_Mass);
}
else if (TMath::Abs(Jpsi_M - PSI_2S_MASS) < 60.)
{
h1_B_M_Psi2s_cut->Fill(used_B_Mass);
}
}
}
else
{
if (TMath::Abs(Jpsi_M - J_PSI_MASS) < 60.)
{
h1_B_M_JPsi_cut->Fill(used_B_Mass);
}
else if (TMath::Abs(Jpsi_M - PSI_2S_MASS) < 60.)
{
h1_B_M_Psi2s_cut->Fill(used_B_Mass);
}
}
}
h1_B_M->GetYaxis()->SetTitle(TString::Format("Events").Data());
h1_B_M_JPsi_cut->GetYaxis()->SetTitle(TString::Format("Events").Data());
h1_Jpsi_M->GetYaxis()->SetTitle(TString::Format("Events").Data());
h1_B_M->SetMinimum(0);
h1_B_M_JPsi_cut->SetMinimum(0);
h1_Jpsi_M->SetMinimum(0);
h2_B_M_Jpsi_M->SetMinimum(0);
h1_B_M->SetStats(0);
h1_B_M_JPsi_cut->SetStats(0);
h1_Jpsi_M->SetStats(0);
h2_B_M_Jpsi_M->SetStats(0);
// auto fitRes = CreateRooFit(h1_B_M_JPsi_cut, ana);
// auto paramsPlot = PlotWithParams(h1_B_M_JPsi_cut, ana);
// savePDF(h1_B_M, ana.name.c_str(), "B_M_uncut");
// savePDF(h1_Jpsi_M, ana.name.c_str(), "JPsi_M_uncut");
// savePDF(h1_B_M_JPsi_cut, ana.name.c_str(), "B_M_JPsi_cut");
// savePDF(h1_Jpsi_M_align, ana.name.c_str(), "JPsi_M_uncut_al");
// savePDF(h2_B_M_Jpsi_M, ana.name.c_str(), "B_M_vs_Jpsi_M");
CreateRooFitAndSavePDF(h1_B_M_JPsi_cut, ana, "B_mass_JPsi_cut_fit");
CreateRooFitAndSavePDF(h1_B_M_Psi2s_cut, ana, "B_mass_Psi2s_cut_fit");
if (ana.B_name == "B0")
{
savePDF(h1_Kst0_M, ana.name.c_str(), "Kst0_M_uncut");
savePDF(h1_Kpi_M, ana.name.c_str(), "Kpi_M_uncut");
}
}
return 0;
}
void CreateRooFitAndSavePDF(TH1D *hist, AnalysisOutput ana, const char *name)
{
RooRealVar roo_var_mass(ana.suf("var_mass").c_str(), TString::Format("m(%s)", ana.B_name.c_str()).Data(), MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX);
roo_var_mass.setRange("fitting_range", MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX);
std::string hist_name = ana.suf("hist_B_M");
RooDataHist roohist_B_M(hist_name.c_str(), "B Mass Histogram", roo_var_mass, RooFit::Import(*hist));
// RooRealVar roo_sig_bw_mean(ana.suf("sig_mean").c_str(), "Mass BW Mean", 5250., 5100., 5400.);
// RooRealVar roo_sig_bw_with(ana.suf("sig_width").c_str(), "Mass BW Width", 20., 0., 50.);
// RooBreitWigner roo_sig_bw(ana.suf("roo_sig_bw").c_str(), "B Signal Breit Wigner", roo_var_mass, roo_sig_bw_mean, roo_sig_bw_with);
RooRealVar roo_sig_gauss_mean(ana.suf("sig_mean").c_str(), "#mu", 5250., 5100., 5400.);
RooRealVar roo_sig_gauss_sigma(ana.suf("sig_sigma").c_str(), "#sigma", 20., 0., 50.);
RooGaussian roo_sig_gauss(ana.suf("sig_gauss").c_str(), "B Signal Gaussian", roo_var_mass, roo_sig_gauss_mean, roo_sig_gauss_sigma);
// RooRealVar roo_sig_tail(ana.suf("sig_tail").c_str(), "#lambda_{sig}", -0.5, -1., 0.);
// RooNovosibirsk roo_sig_nov(ana.suf("sig_nov").c_str(), "B Signal Nov", roo_var_mass, roo_sig_gauss_mean, roo_sig_gauss_sigma, roo_sig_tail);
// RooRealVar roo_sig_add_gau_exp_frac(ana.suf("sig_add_gau_exp_frac").c_str(), "sig exp gau frac", 0.5, 0., 1.);
// RooAddPdf roo_sig_add_gau_exp(ana.suf("sig_add_gau_exp").c_str(), "B Mass Signal Gaus + Exp", roo_sig_gauss, roo_sig_exp, roo_sig_add_gau_exp_frac);
// RooFFTConvPdf roo_sig_conv_gau_exp(ana.suf("sig_conv_gau_exp").c_str(), "Exp (x) Gauss", roo_var_mass, roo_sig_gauss, roo_sig_exp);
RooRealVar roo_bkg_exp_c(ana.suf("bkg_exp_c").c_str(), "#lambda_{bkg}", -0.001145, -0.00199, -0.00100);
RooExponential roo_bkg_exp(ana.suf("bkg_exp").c_str(), "B Mass Background Exp", roo_var_mass, roo_bkg_exp_c);
roo_bkg_exp.asTF()
RooRealVar roo_var_mass_sig_yield(ana.suf("sig_yield").c_str(), "N_{Sig}", 0., hist->GetEntries());
RooRealVar roo_var_mass_bkg_yield(ana.suf("bkg_yield").c_str(), "N_{Bkg}", 0., hist->GetEntries());
std::string pdf_name = ana.suf("pdf_sig_plus_bkg");
RooAddPdf roo_pdf_sig_plus_bkg(pdf_name.c_str(), "Sig + Bkg PDF",
RooArgList(roo_sig_gauss, roo_bkg_exp),
RooArgList(roo_var_mass_sig_yield, roo_var_mass_bkg_yield));
RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title(ana.title_suf("B Mass Fit").c_str()));
roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(nBins), RooFit::Name(hist_name.c_str()));
RooFitResult *fitres = roo_pdf_sig_plus_bkg.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range("fitting_range"));
auto name_fit_func_sig = ana.suf("fit_fsig");
auto name_fit_func_bkg = ana.suf("fit_fbkg");
roo_pdf_sig_plus_bkg.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range("fitting_range"), RooFit::Name(pdf_name.c_str()));
roo_pdf_sig_plus_bkg.plotOn(roo_frame_mass, RooFit::Name(name_fit_func_bkg.c_str()), RooFit::Components(RooArgSet(roo_bkg_exp)), RooFit::LineColor(kBlue), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"));
roo_pdf_sig_plus_bkg.plotOn(roo_frame_mass, RooFit::Name(name_fit_func_sig.c_str()), RooFit::Components(RooArgSet(roo_sig_gauss)), RooFit::FillStyle(3244), RooFit::LineColor(kRed - 7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"));
// roo_pdf_sig_plus_bkg.paramOn(roo_frame_mass, RooFit::Layout(0.45, 0.99, 0.90));
// roo_frame_mass->getAttText()->SetTextSize(0.030);
RooPlot *roo_frame_pull = roo_var_mass.frame(RooFit::Title("."));
roo_frame_pull->addPlotable(roo_frame_mass->pullHist(hist_name.c_str(), pdf_name.c_str()), "P");
Float_t title_size = 0.083;
Float_t label_size = 0.073;
roo_frame_pull->GetXaxis()->SetTitle(TString::Format("m(%s)", ana.B_name.c_str()).Data());
roo_frame_pull->GetXaxis()->SetTitleSize(title_size);
roo_frame_pull->GetYaxis()->SetTitleSize(title_size);
roo_frame_pull->GetXaxis()->SetLabelSize(label_size);
roo_frame_pull->GetYaxis()->SetLabelSize(label_size);
auto cname = TString::Format("%s_%s", ana.name.c_str(), name);
auto c = new TCanvas(cname.Data(), cname.Data(), 0, 0, 800, 600);
auto p2 = new TPad(TString::Format("%s_p2", name), "Lower Pad", 0., 0., 1., 0.3);
p2->Draw();
p2->SetTopMargin(0.001);
p2->SetBottomMargin(0.3);
p2->SetGrid();
auto *p1 = new TPad(TString::Format("%s_p1", name), "Upper Pad", 0., 0.32, 1., 1.);
p1->Draw();
p1->SetBottomMargin(0.001);
p1->cd();
roo_frame_mass->Draw();
TLegend *leg1 = new TLegend(0.58, 0.50, 0.96, 0.87);
// leg1->SetFillColor(kWhite);
leg1->SetLineColor(kWhite);
leg1->AddEntry(roo_frame_mass->findObject(pdf_name.c_str()), "Signal + Background", "LP");
leg1->AddEntry(roo_frame_mass->findObject(name_fit_func_sig.c_str()), "Signal", "LP");
leg1->AddEntry(roo_frame_mass->findObject(name_fit_func_bkg.c_str()), "Background", "LP");
leg1->AddEntry((TObject *)0, TString::Format("%s = %.2f #pm %.2f", roo_sig_gauss_mean.getTitle().Data(), roo_sig_gauss_mean.getVal(), roo_sig_gauss_mean.getError()).Data(), "");
leg1->AddEntry((TObject *)0, TString::Format("%s = %.2f #pm %.2f", roo_sig_gauss_sigma.getTitle().Data(), roo_sig_gauss_sigma.getVal(), roo_sig_gauss_sigma.getError()).Data(), "");
// leg1->AddEntry((TObject *)0, TString::Format("%s = %.8f #pm %.8f", roo_sig_tail.getTitle().Data(), roo_sig_tail.getVal(), roo_sig_tail.getError()).Data(), "");
leg1->AddEntry((TObject *)0, TString::Format("%s = %.8f #pm %.8f", roo_bkg_exp_c.getTitle().Data(), roo_bkg_exp_c.getVal(), roo_bkg_exp_c.getError()).Data(), "");
leg1->AddEntry((TObject *)0, TString::Format("%s = %.2f #pm %.2f", roo_var_mass_sig_yield.getTitle().Data(), roo_var_mass_sig_yield.getVal(), roo_var_mass_sig_yield.getError()).Data(), "");
leg1->AddEntry((TObject *)0, TString::Format("%s = %.2f #pm %.2f", roo_var_mass_bkg_yield.getTitle().Data(), roo_var_mass_bkg_yield.getVal(), roo_var_mass_bkg_yield.getError()).Data(), "");
leg1->Draw();
p2->cd();
roo_frame_pull->Draw();
c->SaveAs(TString::Format("%s/%s.png", SAVE_PATH.c_str(), cname.Data()));
}
// std::vector<RooPlot *> PlotWithParams(TH1D *hist, AnalysisOutput ana)
// {
// RooRealVar roo_var_mass(ana.suf("var_mass").c_str(), TString::Format("m(%s)", ana.B_name).Data(), MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX);
// roo_var_mass.setRange("fitting_range", MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX);
// std::string hist_name = ana.suf("hist_B_M");
// RooDataHist roohist_B_M(hist_name.c_str(), "B Mass Histogram", roo_var_mass, RooFit::Import(*hist));
// // RooRealVar roo_sig_bw_mean(ana.suf("sig_mean").c_str(), "Mass BW Mean", 5250., 5100., 5400.);
// // RooRealVar roo_sig_bw_with(ana.suf("sig_width").c_str(), "Mass BW Width", 20., 0., 50.);
// // RooBreitWigner roo_sig_bw(ana.suf("roo_sig_bw").c_str(), "B Signal Breit Wigner", roo_var_mass, roo_sig_bw_mean, roo_sig_bw_with);
// RooGaussian roo_sig_gauss(ana.suf("sig_gauss").c_str(), "B Signal Breit Wigner", roo_var_mass, RooRealConstant::value(ana.fit_params.mean), RooRealConstant::value(ana.fit_params.sigma));
// RooExponential roo_bkg_exp(ana.suf("bkg_exp").c_str(), "B Mass Background Exp", roo_var_mass, RooRealConstant::value(ana.fit_params.lambda));
// std::string pdf_name = ana.suf("pdf_sig_plus_bkg");
// RooAddPdf roo_pdf_sig_plus_bkg(pdf_name.c_str(), "Sig + Bkg PDF",
// RooArgList(roo_sig_gauss, roo_bkg_exp),
// RooArgList(RooRealConstant::value(ana.fit_params.sig_yield), RooRealConstant::value(ana.fit_params.bkg_yield)));
// RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title(TString::Format("%s [%s]", ana.name.c_str(), ana.title.c_str()).Data()));
// roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(nBins), RooFit::Name(hist_name.c_str()));
// roo_pdf_sig_plus_bkg.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name(pdf_name.c_str()));
// roo_pdf_sig_plus_bkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(roo_bkg_exp)), RooFit::LineColor(kBlue), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"));
// roo_pdf_sig_plus_bkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(roo_sig_gauss)), RooFit::FillStyle(3244), RooFit::LineColor(kRed - 7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"));
// return std::vector<RooPlot *>{roo_frame_mass};
// }