#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};
// }