add a bunch of new analysis macros

.gitignore

@@ -2,3 +2,4 @@
 dataloader/**
 status_report/**
 branchnames.txt
+images/**

analysis_fullstream_b02hphmmumu.cpp

@@ -0,0 +1,264 @@
+#include "TH1D.h"
+#include "TH2D.h"
+#include "THStack.h"
+#include "TGraph.h"
+#include "TTree.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 "TLegend.h"
+#include "TLatex.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 = 5100.;
+const Double_t MASS_HIST_MAX = 5700.;
+const Double_t MASS_HIST_FIT_MIN = 5100.;
+const Double_t MASS_HIST_FIT_MAX = 5700.;
+
+const char* TITLE = "SpruceRD_B0ToHpHmMuMu (#pi^{+} #rightarrow K^{+})";
+const char* FILE_NAME = "SpruceRD_B0ToHpHmMuMu_Pip2Kp";
+const char* MASS_LITERAL = "m(#pi^{+}_{(#rightarrow K^{+})}#pi^{-}#mu^{+}#mu^{-})";
+
+const Double_t K_MASS = 493.677;
+
+void CreateRooFitAndDraw(TH1D *hist, int fitting_entries);
+
+int analysis_fullstream_b02hphmmumu()
+{
+    TChain *data_chain = new TChain("B0ToHpHmMuMu/DecayTree");
+    data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0r1_tuple_90000000_2023_v0r0p6288631.root");
+    data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0_tuple_90000000_v0r0p6288631.root");
+
+    Float_t B0_BPVFDCHI2,
+        B0_BPVIPCHI2,
+        L1_BPVIPCHI2,
+        L2_BPVIPCHI2,
+        L1_PT,
+        L2_PT,
+        Jpsi_BPVFDCHI2,
+        Hp_PT,
+        Hp_BPVIPCHI2,
+        Hp_P,
+        Hm_PT,
+        Hm_BPVIPCHI2,
+        Hm_P,
+        Res_PT,
+        Res_BPVIPCHI2,
+        Res_P;
+
+    Double_t L1_PID_MU,
+        L2_PID_MU,
+        B0_CHI2VXNDOF,
+        Jpsi_MAXDOCACHI2,
+        Jpsi_CHI2DOF,
+        Res_CHI2DOF,
+        Hp_PID_K,
+        Hm_PID_K,
+        Jpsi_M,
+        B0_M,
+        Res_MAXDOCACHI2,
+        Res_M;
+
+    Bool_t L1_ISMUON,
+        L2_ISMUON,
+        Hlt2RD_B0ToKpPimMuMuDecision,
+        Hlt2_InclDetDiMuon_4BodyDecision,
+        Hlt2_InclDetDiMuon_3BodyDecision,
+        Hlt2_InclDetDiMuonDecision,
+        Hlt1TrackMVADecision,
+        Hlt1TwoTrackMVADecision;
+
+    data_chain->SetBranchAddress("B0_M", &B0_M);
+
+    data_chain->SetBranchAddress("B0_BPVFDCHI2", &B0_BPVFDCHI2);
+
+    data_chain->SetBranchAddress("B0_BPVIPCHI2", &B0_BPVIPCHI2);
+
+    data_chain->SetBranchAddress("L1_BPVIPCHI2", &L1_BPVIPCHI2);
+    data_chain->SetBranchAddress("L2_BPVIPCHI2", &L2_BPVIPCHI2);
+    data_chain->SetBranchAddress("L1_PID_MU", &L1_PID_MU);
+    data_chain->SetBranchAddress("L2_PID_MU", &L2_PID_MU);
+    data_chain->SetBranchAddress("L1_ISMUON", &L1_ISMUON);
+    data_chain->SetBranchAddress("L2_ISMUON", &L2_ISMUON);
+    data_chain->SetBranchAddress("L1_PT", &L1_PT);
+    data_chain->SetBranchAddress("L2_PT", &L2_PT);
+
+    data_chain->SetBranchAddress("B0_CHI2VXNDOF", &B0_CHI2VXNDOF);
+    data_chain->SetBranchAddress("Jpsi_MAXDOCACHI2", &Jpsi_MAXDOCACHI2);
+    data_chain->SetBranchAddress("Jpsi_CHI2DOF", &Jpsi_CHI2DOF);
+    data_chain->SetBranchAddress("Res_CHI2DOF", &Res_CHI2DOF);
+
+    data_chain->SetBranchAddress("Hp_PT", &Hp_PT);
+    data_chain->SetBranchAddress("Hp_BPVIPCHI2", &Hp_BPVIPCHI2);
+    data_chain->SetBranchAddress("Hp_P", &Hp_P);
+    data_chain->SetBranchAddress("Hp_PID_K", &Hp_PID_K);
+    data_chain->SetBranchAddress("Hm_PT", &Hm_PT);
+    data_chain->SetBranchAddress("Hm_BPVIPCHI2", &Hm_BPVIPCHI2);
+    data_chain->SetBranchAddress("Hm_P", &Hm_P);
+    data_chain->SetBranchAddress("Hm_PID_K", &Hm_PID_K);
+    data_chain->SetBranchAddress("Res_PT", &Res_PT);
+    data_chain->SetBranchAddress("Res_BPVIPCHI2", &Res_BPVIPCHI2);
+    data_chain->SetBranchAddress("Res_P", &Res_P);
+    data_chain->SetBranchAddress("Res_MAXDOCACHI2", &Res_MAXDOCACHI2);
+    data_chain->SetBranchAddress("Res_M", &Res_M);
+
+
+    data_chain->SetBranchAddress("Jpsi_BPVFDCHI2", &Jpsi_BPVFDCHI2);
+    data_chain->SetBranchAddress("Jpsi_M", &Jpsi_M);
+
+    data_chain->SetBranchAddress("Hlt2RD_B0ToKpPimMuMuDecision", &Hlt2RD_B0ToKpPimMuMuDecision);
+    data_chain->SetBranchAddress("Hlt2_InclDetDiMuon_4BodyDecision", &Hlt2_InclDetDiMuon_4BodyDecision);
+    data_chain->SetBranchAddress("Hlt2_InclDetDiMuon_3BodyDecision", &Hlt2_InclDetDiMuon_3BodyDecision);
+    data_chain->SetBranchAddress("Hlt2_InclDetDiMuonDecision", &Hlt2_InclDetDiMuonDecision);
+    data_chain->SetBranchAddress("Hlt1TrackMVADecision", &Hlt1TrackMVADecision);
+    data_chain->SetBranchAddress("Hlt1TwoTrackMVADecision", &Hlt1TwoTrackMVADecision);
+
+    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,
+        Hm_PX, Hm_PY, Hm_PZ, Hm_ENERGY;
+    
+    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);
+    data_chain->SetBranchAddress("Hm_PX", &Hm_PX);
+    data_chain->SetBranchAddress("Hm_PY", &Hm_PY);
+    data_chain->SetBranchAddress("Hm_PZ", &Hm_PZ);
+    data_chain->SetBranchAddress("Hm_ENERGY", &Hm_ENERGY);
+
+    TH1D *h1_B0_M = new TH1D("h1_B0_M", TITLE, nBins, MASS_HIST_MIN, MASS_HIST_MAX);
+
+    unsigned int entries = data_chain->GetEntries();
+    unsigned int fitting_entries = 0;
+    for (unsigned int i = 0; i < entries; i++)
+    {
+        data_chain->GetEntry(i);
+
+        TLorentzVector l1_4v(L1_PX, L1_PY, L1_PZ, L1_ENERGY);
+        TLorentzVector l2_4v(L2_PX, L2_PY, L2_PZ, L2_ENERGY);
+
+        // Pi+ -> K+
+        TVector3 K_momentum(Hp_PX, Hp_PY, Hp_PZ);
+        double K_energy = TMath::Sqrt(TMath::Sq(K_MASS) + K_momentum.Mag2());
+        TLorentzVector Pi_4v(Hm_PX, Hm_PY, Hm_PZ, Hm_ENERGY);
+
+        // Pi- -> K-
+        // TVector3 K_momentum(Hm_PX, Hm_PY, Hm_PZ);
+        // double K_energy = TMath::Sqrt(TMath::Sq(K_MASS) + K_momentum.Mag2());
+        // TLorentzVector Pi_4v(Hp_PX, Hp_PY, Hp_PZ, Hp_ENERGY);
+
+        TLorentzVector K_4v(K_momentum, K_energy);
+        
+        Double_t reconstructed_Res_Mass = (K_4v + Pi_4v).M();
+        Double_t reconstructed_B0_Mass = (K_4v + Pi_4v + l1_4v + l2_4v).M();
+
+        if ((
+                (B0_BPVFDCHI2 > 64) & (B0_CHI2VXNDOF < 9) & (B0_BPVIPCHI2 < 25) &
+                //
+                (Jpsi_MAXDOCACHI2 < 36) &
+                //
+                (L1_BPVIPCHI2 > 9) & (L2_BPVIPCHI2 > 9) & (L1_PID_MU > -3) & (L2_PID_MU > -3) & (L1_ISMUON == 1) & (L2_ISMUON == 1) & (L1_PT > 350) & (L2_PT > 350) &
+                //
+                (Jpsi_CHI2DOF < 9) & (Jpsi_BPVFDCHI2 > 16) & (Jpsi_M < 5500) &
+                //
+                (Res_PT > 400) & (Res_MAXDOCACHI2 < 36) & (reconstructed_Res_Mass < 2600) & (Res_CHI2DOF < 25) & (Res_BPVIPCHI2 > 4) &
+                //
+                (Hm_BPVIPCHI2 > 6) & (Hm_PT > 250) & (Hm_P > 2000) & (Hm_PID_K < 0) &
+                //
+                (Hp_BPVIPCHI2 > 6) & (Hp_PT > 250) & (Hp_P > 2000) & (Hp_PID_K > 0)
+            ) &
+            //
+            (Hlt2RD_B0ToKpPimMuMuDecision) & ((Hlt2_InclDetDiMuon_4BodyDecision) | (Hlt2_InclDetDiMuon_3BodyDecision) | (Hlt2_InclDetDiMuonDecision)))
+        {
+            if ((TMath::Abs(Jpsi_M - 3096.9) < 100) & (Hlt1TrackMVADecision) & (Hlt1TwoTrackMVADecision))
+            {
+                h1_B0_M->Fill(reconstructed_B0_Mass);
+
+                if (MASS_HIST_FIT_MIN <= reconstructed_B0_Mass && reconstructed_B0_Mass <= MASS_HIST_FIT_MAX)
+                {
+                    fitting_entries++;
+                }
+            }
+        }
+
+        if (i % 10000 == 0)
+        {
+            std::cout << "["
+                      << FILE_NAME
+                      << "] Processed event: " << i << " (" << TString::Format("%.2f", ((double)i / (double)entries) * 100.) << "%)" << std::endl;
+        }
+    }
+
+    h1_B0_M->GetXaxis()->SetTitle(MASS_LITERAL);
+
+    TCanvas *c1 = new TCanvas("c1", "c1", 0, 0, 800, 600);
+    h1_B0_M->Draw();
+    c1->Draw();
+    c1->SaveAs(TString::Format("images/root_hist_%s_bmass.pdf", FILE_NAME).Data());
+
+    CreateRooFitAndDraw(h1_B0_M, fitting_entries);
+
+    return 0;
+}
+
+void CreateRooFitAndDraw(TH1D *hist, int fitting_entries)
+{
+    RooRealVar roo_var_mass("var_mass", "B0 Mass Variable", MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX);
+    roo_var_mass.setRange("fitting_range", MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX);
+
+    RooDataHist roohist_B0_M("roohist_B0_M", "B0 Mass Histogram", roo_var_mass, RooFit::Import(*hist));
+
+    RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title(TITLE));
+    roohist_B0_M.plotOn(roo_frame_mass, RooFit::Binning(nBins), RooFit::Name("B Mass Distribution"));
+    roo_frame_mass->GetXaxis()->SetTitle(MASS_LITERAL);
+
+    TCanvas *c = new TCanvas("roofit_c", "roofit_c", 0, 0, 800, 600);
+
+    roo_frame_mass->Draw();
+
+    TLegend *leg1 = new TLegend(0.65, 0.7, 0.87, 0.8);
+    leg1->SetFillColor(kWhite);
+    leg1->SetLineColor(kBlack);
+
+    leg1->AddEntry((TObject *)0, TString::Format("Entries: %d", fitting_entries), "");
+
+    leg1->Draw();
+
+    c->Draw();
+    c->SaveAs(TString::Format("images/roofit_hist_%s_bmass.pdf", FILE_NAME).Data());
+}

analysis_fullstream_bu2hpmumu.cpp

@@ -0,0 +1,200 @@
+#include "TH1D.h"
+#include "TH2D.h"
+#include "THStack.h"
+#include "TGraph.h"
+#include "TTree.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 "TLegend.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 = 5100.;
+const Double_t MASS_HIST_MAX = 6000.;
+const Double_t MASS_HIST_FIT_MIN = 5100.;
+const Double_t MASS_HIST_FIT_MAX = 6000.;
+
+const Double_t K_MASS = 493.677;
+
+void CreateRooFitAndDraw(TH1D *hist, int fitting_entries);
+
+const char* TITLE = "SpruceRD_BuToHpMuMu (#pi^{+} #rightarrow K^{+})";
+const char* FILE_NAME = "SpruceRD_BuToHpMuMu_Pip2Kp";
+const char* MASS_LITERAL = "m(#pi^{+}_{(#rightarrow K^{+})}#mu^{+}#mu^{-})";
+
+int analysis_fullstream_bu2hpmumu()
+{
+    TChain *data_chain = new TChain("BuToHpMuMu/DecayTree");
+    data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0r1_tuple_90000000_2023_v0r0p6288631.root");
+    data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0_tuple_90000000_v0r0p6288631.root");
+
+    Float_t B_BPVFDCHI2,
+        B_BPVIPCHI2,
+        L1_BPVIPCHI2,
+        L2_BPVIPCHI2,
+        L1_PT,
+        L2_PT,
+        Jpsi_BPVFDCHI2,
+        Hp_PT,
+        Hp_BPVIPCHI2,
+        Hp_P;
+
+    Double_t L1_PID_MU,
+        L2_PID_MU,
+        B_CHI2VXNDOF,
+        Jpsi_MAXDOCACHI2,
+        Jpsi_CHI2DOF,
+        Hp_PID_K,
+        Jpsi_M,
+        B_M;
+
+    Bool_t L1_ISMUON,
+        L2_ISMUON,
+        Hlt2RD_BuToKpMuMuDecision,
+        Hlt2_InclDetDiMuon_4BodyDecision,
+        Hlt2_InclDetDiMuon_3BodyDecision,
+        Hlt2_InclDetDiMuonDecision,
+        Hlt1TrackMVADecision,
+        Hlt1TwoTrackMVADecision;
+
+    data_chain->SetBranchAddress("B_M", &B_M);
+
+    data_chain->SetBranchAddress("B_BPVFDCHI2", &B_BPVFDCHI2);
+
+    data_chain->SetBranchAddress("B_BPVIPCHI2", &B_BPVIPCHI2);
+
+    data_chain->SetBranchAddress("L1_BPVIPCHI2", &L1_BPVIPCHI2);
+    data_chain->SetBranchAddress("L2_BPVIPCHI2", &L2_BPVIPCHI2);
+    data_chain->SetBranchAddress("L1_PID_MU", &L1_PID_MU);
+    data_chain->SetBranchAddress("L2_PID_MU", &L2_PID_MU);
+    data_chain->SetBranchAddress("L1_ISMUON", &L1_ISMUON);
+    data_chain->SetBranchAddress("L2_ISMUON", &L2_ISMUON);
+    data_chain->SetBranchAddress("L1_PT", &L1_PT);
+    data_chain->SetBranchAddress("L2_PT", &L2_PT);
+    
+    data_chain->SetBranchAddress("B_CHI2VXNDOF", &B_CHI2VXNDOF);
+    data_chain->SetBranchAddress("Jpsi_MAXDOCACHI2", &Jpsi_MAXDOCACHI2);
+    data_chain->SetBranchAddress("Jpsi_CHI2DOF", &Jpsi_CHI2DOF);
+
+    data_chain->SetBranchAddress("Hp_PT", &Hp_PT);
+    data_chain->SetBranchAddress("Hp_BPVIPCHI2", &Hp_BPVIPCHI2);
+    data_chain->SetBranchAddress("Hp_P", &Hp_P);
+    data_chain->SetBranchAddress("Hp_PID_K", &Hp_PID_K);
+
+    data_chain->SetBranchAddress("Jpsi_BPVFDCHI2", &Jpsi_BPVFDCHI2);
+    data_chain->SetBranchAddress("Jpsi_M", &Jpsi_M);
+
+    data_chain->SetBranchAddress("Hlt2RD_BuToKpMuMuDecision", &Hlt2RD_BuToKpMuMuDecision);
+    data_chain->SetBranchAddress("Hlt2_InclDetDiMuon_4BodyDecision", &Hlt2_InclDetDiMuon_4BodyDecision);
+    data_chain->SetBranchAddress("Hlt2_InclDetDiMuon_3BodyDecision", &Hlt2_InclDetDiMuon_3BodyDecision);
+    data_chain->SetBranchAddress("Hlt2_InclDetDiMuonDecision", &Hlt2_InclDetDiMuonDecision);
+    data_chain->SetBranchAddress("Hlt1TrackMVADecision", &Hlt1TrackMVADecision);
+    data_chain->SetBranchAddress("Hlt1TwoTrackMVADecision", &Hlt1TwoTrackMVADecision);
+
+    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;
+    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);
+
+    TH1D *h1_B_M = new TH1D("h1_B_M", "B Mass", nBins, MASS_HIST_MIN, MASS_HIST_MAX);
+
+    unsigned int entries = data_chain->GetEntries();
+    unsigned int fitting_entries = 0;
+    for (unsigned int i = 0; i < entries; i++)
+    {
+        data_chain->GetEntry(i);
+
+        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);
+        Double_t reconstructed_B_Mass = (K_4v + l1_4v + l2_4v).M();
+
+        if ((((B_BPVFDCHI2 > 36) & (B_CHI2VXNDOF<16) & ( B_BPVIPCHI2<25) & (Jpsi_MAXDOCACHI2 < 36 )&(L1_BPVIPCHI2 > 9) & (L2_BPVIPCHI2 > 9) & (L1_PID_MU > -3) & (L2_PID_MU > -3)& (L1_ISMUON==1) & (L2_ISMUON==1) & (L1_PT>350) & (L2_PT>350) & (Jpsi_CHI2DOF<9 ) & (Jpsi_BPVFDCHI2>16) & (Jpsi_M<5500) & (Hp_PT>400) &( Hp_BPVIPCHI2>6) & (Hp_PT>400) & (Hp_BPVIPCHI2>6) & (Hp_P>2000) & (Hp_PID_K > -4)) & (Hlt2RD_BuToKpMuMuDecision==1)  &(((Hlt2_InclDetDiMuon_4BodyDecision==1) | (Hlt2_InclDetDiMuon_3BodyDecision==1) | (Hlt2_InclDetDiMuonDecision==1))))&(TMath::Abs(Jpsi_M - 3096.9) < 100) & ((Hlt1TrackMVADecision==1) | (Hlt1TwoTrackMVADecision==1))) 
+        {
+            h1_B_M->Fill(reconstructed_B_Mass);
+
+            if (MASS_HIST_FIT_MIN <= reconstructed_B_Mass && reconstructed_B_Mass <= MASS_HIST_FIT_MAX) {
+                fitting_entries++;
+            }
+        }
+
+        if (i % 10000 == 0)
+        {
+            std::cout << "["
+                      << "SpruceRD_BuToHpMuMu"
+                      << "] Processed event: " << i << " (" << TString::Format("%.2f", ((double)i / (double)entries) * 100.) << "%)" << std::endl;
+        }
+    }
+
+    TCanvas *c1 = new TCanvas("c1", "c1", 0, 0, 800, 600);
+    h1_B_M->Draw();
+    c1->Draw();
+    c1->SaveAs(TString::Format("images/root_hist_%s_bmass.png", FILE_NAME).Data());
+
+    CreateRooFitAndDraw(h1_B_M, fitting_entries);
+
+    return 0;
+}
+
+void CreateRooFitAndDraw(TH1D *hist, int fitting_entries)
+{
+    RooRealVar roo_var_mass("var_mass", "B Mass Variable", MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX);
+    roo_var_mass.setRange("fitting_range", MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX);
+
+    RooDataHist roohist_B_M("roohist_B_M", "B Mass Histogram", roo_var_mass, RooFit::Import(*hist));
+
+    RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title(TITLE));
+    roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(nBins), RooFit::Name("B Mass Distribution"));
+    roo_frame_mass->GetXaxis()->SetTitle(MASS_LITERAL);
+
+    TCanvas *c = new TCanvas("roofit_c", "roofit_c", 0, 0, 800, 600);
+
+    roo_frame_mass->Draw();
+
+    TLegend *leg1 = new TLegend(0.65, 0.7, 0.87, 0.8);
+    leg1->SetFillColor(kWhite);
+    leg1->SetLineColor(kBlack);
+
+    leg1->AddEntry((TObject*)0, TString::Format("Entries: %d", fitting_entries), "");
+
+    leg1->Draw();
+
+    c->Draw();
+    c->SaveAs(TString::Format("images/roofit_hist_%s_bmass.png", FILE_NAME).Data());
+}

analysis_fullstream_michele.cpp

@@ -0,0 +1,193 @@
+#include "TH1D.h"
+#include "TH2D.h"
+#include "THStack.h"
+#include "TGraph.h"
+#include "TTree.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 "TLegend.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 = 4500.;
+const Double_t MASS_HIST_MAX = 7000.;
+const Double_t MASS_HIST_FIT_MIN = 4500.;
+const Double_t MASS_HIST_FIT_MAX = 7000.;
+
+void CreateRooFitAndDraw(TH1D *hist, const char *title);
+
+int analysis_fullstream()
+{
+
+    std::string title {"SpruceRD_BuToHpMuMu"};
+
+    TChain *data_chain = new TChain("BuToHpMuMu/DecayTree");
+    // data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_Sprucing23r1_90000000_RD.root");
+    data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0r1_tuple_90000000_2023_v0r0p6288631.root");
+    // data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0_tuple_90000000_v0r0p6288631.root");
+
+    Float_t B_BPVFDCHI2,
+        B_BPVIPCHI2,
+        L1_BPVIPCHI2,
+        L2_BPVIPCHI2,
+        L1_PT,
+        L2_PT,
+        Jpsi_BPVFDCHI2,
+        Hp_PT,
+        Hp_BPVIPCHI2,
+        Hp_P;
+
+    Double_t L1_PID_MU,
+        L2_PID_MU,
+        B_CHI2VXNDOF,
+        Jpsi_MAXDOCACHI2,
+        Jpsi_CHI2DOF,
+        Hp_PID_K,
+        Jpsi_M,
+        B_M;
+
+    Bool_t L1_ISMUON,
+        L2_ISMUON,
+        Hlt2RD_BuToKpMuMuDecision,
+        Hlt2_InclDetDiMuon_4BodyDecision,
+        Hlt2_InclDetDiMuon_3BodyDecision,
+        Hlt2_InclDetDiMuonDecision,
+        Hlt1TrackMVADecision,
+        Hlt1TwoTrackMVADecision;
+
+    data_chain->SetBranchAddress("B_M", &B_M);
+
+    data_chain->SetBranchAddress("B_BPVFDCHI2", &B_BPVFDCHI2);
+
+    data_chain->SetBranchAddress("B_BPVIPCHI2", &B_BPVIPCHI2);
+
+    // data_chain->SetBranchAddress("L1_BPVIPCHI2", &L1_BPVIPCHI2);
+    // data_chain->SetBranchAddress("L2_BPVIPCHI2", &L2_BPVIPCHI2);
+    // data_chain->SetBranchAddress("L1_PID_MU", &L1_PID_MU);
+    // data_chain->SetBranchAddress("L2_PID_MU", &L2_PID_MU);
+    // data_chain->SetBranchAddress("L1_ISMUON", &L1_ISMUON);
+    // data_chain->SetBranchAddress("L2_ISMUON", &L2_ISMUON);
+    // data_chain->SetBranchAddress("L1_PT", &L1_PT);
+    // data_chain->SetBranchAddress("L2_PT", &L2_PT);
+    // data_chain->SetBranchAddress("Hp_PT", &Hp_PT);
+    // data_chain->SetBranchAddress("Hp_BPVIPCHI2", &Hp_BPVIPCHI2);
+    // data_chain->SetBranchAddress("Hp_P", &Hp_P);
+    // data_chain->SetBranchAddress("Hp_PID_K", &Hp_PID_K);
+
+    data_chain->SetBranchAddress("L1_BPVIPCHI2", &L1_BPVIPCHI2);
+    data_chain->SetBranchAddress("L2_BPVIPCHI2", &L2_BPVIPCHI2);
+    data_chain->SetBranchAddress("L1_PID_MU", &L1_PID_MU);
+    data_chain->SetBranchAddress("L2_PID_MU", &L2_PID_MU);
+    data_chain->SetBranchAddress("L1_ISMUON", &L1_ISMUON);
+    data_chain->SetBranchAddress("L2_ISMUON", &L2_ISMUON);
+    data_chain->SetBranchAddress("L1_PT", &L1_PT);
+    data_chain->SetBranchAddress("L2_PT", &L2_PT);
+    
+    data_chain->SetBranchAddress("B_CHI2VXNDOF", &B_CHI2VXNDOF);
+    data_chain->SetBranchAddress("Jpsi_MAXDOCACHI2", &Jpsi_MAXDOCACHI2);
+    data_chain->SetBranchAddress("Jpsi_CHI2DOF", &Jpsi_CHI2DOF);
+
+    data_chain->SetBranchAddress("Hp_PT", &Hp_PT);
+    data_chain->SetBranchAddress("Hp_BPVIPCHI2", &Hp_BPVIPCHI2);
+    data_chain->SetBranchAddress("Hp_P", &Hp_P);
+    data_chain->SetBranchAddress("Hp_PID_K", &Hp_PID_K);
+
+    data_chain->SetBranchAddress("Jpsi_BPVFDCHI2", &Jpsi_BPVFDCHI2);
+    data_chain->SetBranchAddress("Jpsi_M", &Jpsi_M);
+
+    data_chain->SetBranchAddress("Hlt2RD_BuToKpMuMuDecision", &Hlt2RD_BuToKpMuMuDecision);
+    data_chain->SetBranchAddress("Hlt2_InclDetDiMuon_4BodyDecision", &Hlt2_InclDetDiMuon_4BodyDecision);
+    data_chain->SetBranchAddress("Hlt2_InclDetDiMuon_3BodyDecision", &Hlt2_InclDetDiMuon_3BodyDecision);
+    data_chain->SetBranchAddress("Hlt2_InclDetDiMuonDecision", &Hlt2_InclDetDiMuonDecision);
+    data_chain->SetBranchAddress("Hlt1TrackMVADecision", &Hlt1TrackMVADecision);
+    data_chain->SetBranchAddress("Hlt1TwoTrackMVADecision", &Hlt1TwoTrackMVADecision);
+
+    TH1D *h1_B_M = new TH1D("h1_B_M", "B Mass", nBins, MASS_HIST_MIN, MASS_HIST_MAX);
+
+    unsigned int entries = data_chain->GetEntries();
+    for (unsigned int i = 0; i < entries; i++)
+    {
+        data_chain->GetEntry(i);
+
+        // if (((B_BPVFDCHI2 > 36) && (B_CHI2VXNDOF < 16) && (B_BPVIPCHI2 < 25) && (Jpsi_MAXDOCACHI2 < 36) && (L1_BPVIPCHI2 > 9) && (L2_BPVIPCHI2 > 9) && (L1_PID_MU > -3) && (L2_PID_MU > -3) && (L1_ISMUON) && (L2_ISMUON) && (L1_PT > 350) && (L2_PT > 350) && (Jpsi_CHI2DOF < 9) && (Jpsi_BPVFDCHI2 > 16) && (Jpsi_M < 5500) && (Hp_PT > 400) && (Hp_BPVIPCHI2 > 6) && (Hp_PT > 400) && (Hp_BPVIPCHI2 > 6) && (Hp_P > 2000) && (Hp_PID_K > -4)) && (Hlt2RD_BuToKpMuMuDecision) && (((Hlt2_InclDetDiMuon_4BodyDecision) || (Hlt2_InclDetDiMuon_3BodyDecision) || (Hlt2_InclDetDiMuonDecision))))
+        // {
+        //     if (TMath::Abs(Jpsi_M - 3096.9) < 100 && ((Hlt1TrackMVADecision) || (Hlt1TwoTrackMVADecision)))
+        //     {
+        //         h1_B_M->Fill(B_M);
+        //     }
+        // }
+
+        if ((((B_BPVFDCHI2 > 36) & (B_CHI2VXNDOF<16) & ( B_BPVIPCHI2<25) & (Jpsi_MAXDOCACHI2 < 36 )&(L1_BPVIPCHI2 > 9) & (L2_BPVIPCHI2 > 9) & (L1_PID_MU > -3) & (L2_PID_MU > -3)& (L1_ISMUON==1) & (L2_ISMUON==1) & (L1_PT>350) & (L2_PT>350) & (Jpsi_CHI2DOF<9 ) & (Jpsi_BPVFDCHI2>16) & (Jpsi_M<5500) & (Hp_PT>400) &( Hp_BPVIPCHI2>6) & (Hp_PT>400) & (Hp_BPVIPCHI2>6) & (Hp_P>2000) & (Hp_PID_K > -4)) & (Hlt2RD_BuToKpMuMuDecision==1)  &(((Hlt2_InclDetDiMuon_4BodyDecision==1) | (Hlt2_InclDetDiMuon_3BodyDecision==1) | (Hlt2_InclDetDiMuonDecision==1))))&(abs(Jpsi_M - 3096.9) < 100) & ((Hlt1TrackMVADecision==1) | (Hlt1TwoTrackMVADecision==1))) 
+        {
+            h1_B_M->Fill(B_M);
+        }
+
+        if (i % 10000 == 0)
+        {
+            std::cout << "["
+                      << "SpruceRD_BuToHpMuMu"
+                      << "] Processed event: " << i << " (" << TString::Format("%.2f", ((double)i / (double)entries) * 100.) << "%)" << std::endl;
+        }
+    }
+
+    TCanvas *c1 = new TCanvas("c1", "c1", 0, 0, 800, 600);
+    h1_B_M->Draw();
+    c1->Draw();
+    c1->SaveAs(TString::Format("images/root_hist_%s_bmass.png", title.c_str()).Data());
+
+    CreateRooFitAndDraw(h1_B_M, title.c_str());
+
+    return 0;
+}
+
+void CreateRooFitAndDraw(TH1D *hist, const char *title)
+{
+    RooRealVar roo_var_mass("var_mass", "B Mass Variable", MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX);
+    roo_var_mass.setRange("fitting_range", MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX);
+
+    RooDataHist roohist_B_M("roohist_B_M", "B Mass Histogram", roo_var_mass, RooFit::Import(*hist));
+
+    RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title("SpruceRD_BuToHpMuMu"));
+    roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(nBins), RooFit::Name("B Mass Distribution"));
+    roo_frame_mass->GetXaxis()->SetTitle("m(#pi^{+}#mu^{+}#mu^{-})");
+
+    TCanvas *c = new TCanvas("roofit_c", "roofit_c", 0, 0, 800, 600);
+
+    roo_frame_mass->Draw();
+
+    TLegend *leg1 = new TLegend(0.65, 0.7, 0.87, 0.8);
+    leg1->SetFillColor(kWhite);
+    leg1->SetLineColor(kBlack);
+
+    leg1->AddEntry((TObject*)0, TString::Format("Entries: %d", (int)hist->GetEntries()), "");
+
+    leg1->Draw();
+
+    c->Draw();
+    c->SaveAs(TString::Format("images/roofit_hist_%s_bmass.png", title).Data());
+}

analysis_turbo.cpp

@@ -0,0 +1,141 @@
+#include "TH1D.h"
+#include "TH2D.h"
+#include "THStack.h"
+#include "TGraph.h"
+#include "TTree.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 = 5100.;
+const Double_t MASS_HIST_MAX = 6000.;
+const Double_t MASS_HIST_FIT_MIN = 5100.;
+const Double_t MASS_HIST_FIT_MAX = 6000.;
+
+void CreateRooFitAndDraw(TH1D *hist, const char *title, const char *x_axis);
+
+int analysis_turbo()
+{
+    // std::string title{"Hlt2RD_B0ToKpPimMuMu"};
+    // std::string x_axis{"m(#pi^{+}_{(#rightarrow K^{+})} #pi^{-} #mu^{+}#mu^{-})"};
+
+    std::string title{"Hlt2RD_BuToKpMuMu"};
+    std::string x_axis{"m(#pi^{+}_{(#rightarrow K^{+})}#mu^{+}#mu^{-})"};
+
+    TChain *data_chain = new TChain(TString::Format("%s/DecayTree", title.c_str()).Data());
+    data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_RD.root");
+
+    Double_t B_M, Jpsi_M, B0_M, Kst_M;
+
+    Bool_t Hlt1TrackMVADecision, Hlt1TwoTrackMVADecision;
+
+    if (title == "Hlt2RD_B0ToKpPimMuMu")
+    {
+        data_chain->SetBranchAddress("B0_M", &B0_M);
+        data_chain->SetBranchAddress("Jpsi_M", &Jpsi_M);
+        data_chain->SetBranchAddress("Kst0_M", &Kst_M);
+        data_chain->SetBranchAddress("Hlt1TrackMVADecision", &Hlt1TrackMVADecision);
+        data_chain->SetBranchAddress("Hlt1TwoTrackMVADecision", &Hlt1TwoTrackMVADecision);
+    }
+    else
+    {
+        data_chain->SetBranchAddress("B_M", &B_M);
+        data_chain->SetBranchAddress("Jpsi_M", &Jpsi_M);
+        data_chain->SetBranchAddress("Hlt1TrackMVADecision", &Hlt1TrackMVADecision);
+        data_chain->SetBranchAddress("Hlt1TwoTrackMVADecision", &Hlt1TwoTrackMVADecision);
+    }
+
+    TH1D *h1_B_M = new TH1D("h1_B_M", title.c_str(), nBins, MASS_HIST_MIN, MASS_HIST_MAX);
+    h1_B_M->GetXaxis()->SetTitle(x_axis.c_str());
+
+    unsigned int entries = data_chain->GetEntries();
+    for (unsigned int i = 0; i < entries; i++)
+    {
+        data_chain->GetEntry(i);
+
+        if (title == "Hlt2RD_B0ToKpPimMuMu")
+        {
+            if ((TMath::Abs(Jpsi_M - 3096.9) < 100) && (B0_M > 4500) && (B0_M < 6000) && (TMath::Abs(Kst_M - 895.55) < 50) && ((Hlt1TrackMVADecision) || (Hlt1TwoTrackMVADecision)))
+            {
+                h1_B_M->Fill(B0_M);
+            }
+        }
+        else
+        {
+            if ((TMath::Abs(Jpsi_M - 3096.9) < 100.) && (B_M > 4500.) && (B_M < 6000.) && ((Hlt1TrackMVADecision) || (Hlt1TwoTrackMVADecision)))
+            {
+                h1_B_M->Fill(B_M);
+            }
+        }
+
+        if (i % 10000 == 0)
+        {
+            std::cout << "["
+                      << title
+                      << "] Processed event: " << i << " (" << TString::Format("%.2f", ((double)i / (double)entries) * 100.) << "%)" << std::endl;
+        }
+    }
+
+    TCanvas *c1 = new TCanvas("c1", "c1", 0, 0, 800, 600);
+    h1_B_M->Draw();
+    c1->Draw();
+    c1->SaveAs(TString::Format("images/root_hist_%s_bmass.pdf", title.c_str()).Data());
+
+    CreateRooFitAndDraw(h1_B_M, title.c_str(), x_axis.c_str());
+
+    return 0;
+}
+
+void CreateRooFitAndDraw(TH1D *hist, const char *title, const char *x_axis)
+{
+    RooRealVar roo_var_mass("var_mass", "B Mass Variable", MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX);
+    roo_var_mass.setRange("fitting_range", MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX);
+
+    RooDataHist roohist_B_M("roohist_B_M", "B Mass Histogram", roo_var_mass, RooFit::Import(*hist));
+
+    RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title(title));
+    roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(nBins), RooFit::Name("B Mass Distribution"));
+    roo_frame_mass->GetXaxis()->SetTitle(x_axis);
+
+    TCanvas *c = new TCanvas("roofit_c", "roofit_c", 0, 0, 800, 600);
+
+    roo_frame_mass->Draw();
+
+    TLegend *leg1 = new TLegend(0.65, 0.7, 0.87, 0.8);
+    leg1->SetFillColor(kWhite);
+    leg1->SetLineColor(kBlack);
+
+    leg1->AddEntry((TObject*)0, TString::Format("Entries: %d", (int)hist->GetEntries()), "");
+
+    leg1->Draw();
+
+    c->Draw();
+    c->SaveAs(TString::Format("images/roofit_hist_%s_bmass.pdf", title).Data());
+}

simulation.cpp

@@ -0,0 +1,184 @@
+#include "TH1D.h"
+#include "TH2D.h"
+#include "THStack.h"
+#include "TGraph.h"
+#include "TTree.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>
+#include <filesystem>
+
+const int nBins = 70;
+
+const Double_t MASS_HIST_MIN = 5100.;
+const Double_t MASS_HIST_MAX = 6000.;
+const Double_t MASS_HIST_FIT_MIN = 5100.;
+const Double_t MASS_HIST_FIT_MAX = 6000.;
+
+const Double_t K_MASS = 493.677;
+
+int simulation()
+{
+    const char *name = "BuToKpMuMu_Fullstream";
+    const bool use_hyp_replace = true;
+    const char *b_name = "B";
+    const bool has_H_minus = false;
+
+    auto quickformat_name = [name](const char * format) {
+        return TString::Format(format, name);
+    };
+
+    TChain *data_chain = new TChain("B0ToHpHmMuMu_noPID/DecayTree");
+    data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/rd_btoxll_simulation_fullstream_v0r0p6671378_BuToKpMuMu_12143001_magdown.root");
+
+    /*
+    rd_btoxll_simulation_fullstream_v0r0p6671378_B0ToKpPimMuMu_11144002_magdown.root
+    rd_btoxll_simulation_fullstream_v0r0p6671378_BuToKpMuMu_12143001_magdown.root
+    rd_btoxll_simulation_turbo_v0r0p6657752_B0ToKpPimMuMu_11144002_magdown.root
+    rd_btoxll_simulation_turbo_v0r0p6657752_BuToKpMuMu_12143001_magdown.root
+    */
+
+    Int_t B_BKGCAT;
+    data_chain->SetBranchAddress(TString::Format("%s_BKGCAT", b_name).Data(), &B_BKGCAT);
+
+    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,
+            Hm_PX, Hm_PY, Hm_PZ, Hm_ENERGY;
+    Double_t B_M;
+
+    if (use_hyp_replace)
+    {
+        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);
+        if (has_H_minus) {
+            data_chain->SetBranchAddress("Hm_PX", &Hp_PX);
+            data_chain->SetBranchAddress("Hm_PY", &Hp_PY);
+            data_chain->SetBranchAddress("Hm_PZ", &Hp_PZ);
+            data_chain->SetBranchAddress("Hm_ENERGY", &Hp_ENERGY);
+        }
+    }
+    else
+    {
+        data_chain->SetBranchAddress(TString::Format("%s_M", b_name).Data(), &B_M);
+    }
+
+    TH1D *h1_B_M = new TH1D("h1_B_M", "B Mass", nBins, MASS_HIST_MIN, MASS_HIST_MAX);
+    TH1D *h1_B_M_bkgcat = new TH1D("h1_B_M_bkgcat", "B Mass with BKGCAT <= 10", nBins, MASS_HIST_MIN, MASS_HIST_MAX);
+    TH1I *h1_B_BKGCAT = new TH1I("h1_B_BKGCAT", "B Background Category", nBins, 0, 120);
+    TH2D *h2_B_M_vs_B_BKGCAT = new TH2D("h2_B_M_vs_B_BKGCAT", "B Mass vs Background Category", nBins, MASS_HIST_MIN, MASS_HIST_MAX, nBins, 0, 120);
+
+    unsigned int entries = data_chain->GetEntries();
+    for (size_t i = 0; i < entries; i++)
+    {
+        data_chain->GetEntry(i);
+        // std::cout << B_BKGCAT << std::endl;
+
+        Double_t reconstructed_B_Mass = 0;
+
+        if (use_hyp_replace)
+        {
+            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);
+            TLorentzVector sum = K_4v + l1_4v + l2_4v;
+            if (has_H_minus) {
+                TLorentzVector Hm_4v(Hm_PX, Hm_PY, Hm_PZ, Hm_ENERGY);
+                sum += Hm_4v;
+            }
+            reconstructed_B_Mass = sum.M();
+        }
+        else
+        {
+            reconstructed_B_Mass = B_M;
+        }
+
+        if (B_BKGCAT <= 10)
+        {
+            h1_B_M_bkgcat->Fill(reconstructed_B_Mass);
+        }
+
+        h1_B_M->Fill(reconstructed_B_Mass);
+        h1_B_BKGCAT->Fill(B_BKGCAT);
+        h2_B_M_vs_B_BKGCAT->Fill(reconstructed_B_Mass, B_BKGCAT);
+
+        if (i % 10000 == 0)
+        {
+            std::cout << "["
+                      << name
+                      << "] Processed event: " << i << " (" << TString::Format("%.2f", ((double)i / (double)entries) * 100.) << "%)" << std::endl;
+        }
+    }
+
+    h1_B_M->SetLineColor(kBlue);
+    h1_B_M_bkgcat->SetLineColor(kMagenta + 2);
+    h1_B_M_bkgcat->SetFillColor(kMagenta + 2);
+    h1_B_M_bkgcat->SetFillStyle(3244);
+
+    h1_B_BKGCAT->SetLineColor(kBlue);
+    h1_B_BKGCAT->SetFillColor(kBlue);
+    h1_B_BKGCAT->SetFillStyle(3351);
+
+    std::filesystem::create_directory(TString::Format("images/sim/%s", name).Data());
+
+    TCanvas *c1 = new TCanvas("c1", "c1", 0, 0, 800, 600);
+    h1_B_M->SetStats(0);
+    h1_B_M_bkgcat->SetStats(0);
+    h1_B_M->Draw();
+    h1_B_M_bkgcat->Draw("SAME");
+    c1->BuildLegend(0.58, 0.65, 0.85, 0.87);
+    c1->Draw();
+    c1->SaveAs(quickformat_name("images/sim/%s/h1_B_M.pdf"));
+
+    TCanvas *c2 = new TCanvas("c2", "c2", 0, 0, 800, 600);
+    h1_B_BKGCAT->SetStats(0);
+    h1_B_BKGCAT->Draw();
+    c2->BuildLegend(0.58, 0.77, 0.85, 0.87);
+    c2->Draw();
+    c2->SaveAs(quickformat_name("images/sim/%s/h1_B_BKGCAT.pdf"));
+
+    TCanvas *c3 = new TCanvas("c3", "c3", 0, 0, 800, 600);
+    h2_B_M_vs_B_BKGCAT->SetStats(0);
+    h2_B_M_vs_B_BKGCAT->Draw("COLZ");
+    c3->Draw();
+    c3->SaveAs(quickformat_name("images/sim/%s/h2_B_M_vs_B_BKGCAT.pdf"));
+
+    return 0;
+}

status_report_plots.cpp

@@ -24,78 +24,93 @@
 #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 AnalysisOutput {
+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;
-    std::string paper_yield;
+    int k_charge;
 
-    const std::string suf(std::string input) {
+    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) {
+    const std::string title_suf(std::string input)
+    {
         return TString::Format("%s [%s]", input.c_str(), title.c_str()).Data();
     }
 };
 
-std::vector<RooPlot*> CreateRooFit(TH1D *hist, AnalysisOutput ana);
+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) {
+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) {
+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 = "") {
+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, 800);
+    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();
 }
 
-void savePDF(TObject* o1, TObject* o2, const char* decay, const char* name) {
-    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, 800);
-    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();
-
-    o1->Draw();
-
-    p2->cd();
-
-    o2->Draw();
-    c->SaveAs(TString::Format("%s/%s.jpg", SAVE_PATH.c_str(), cname.Data()));
-}
-
 RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING);
 
 int status_report_plots()
@@ -103,15 +118,16 @@ 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/BuToKpMuMu_Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_Sprucing23r1_90000000_RD.root", "SpruceRD_BuToKpMuMu/DecayTree", "Bplus", ""},
-        // AnalysisOutput{"BuToKpMuMu Inclusive, AnaProd NTuple", "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/BuToKpMuMu_Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_RD.root", "Hlt2RD_BuToKpMuMu_2023/DecayTree", "Bplus", "803.769 +/- 34" },
-        // 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", "B0ToKpPimMuMu_incl", "/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/B0ToKpPimMuMu_Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_Sprucing23r1_90000000_RD.root", "SpruceRD_B0ToKpPimMuMu/DecayTree", "B0", "" },
-        // AnalysisOutput{"B0ToKpPimMuMu_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/B0ToKpPimMuMu_Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_RD.root", "Hlt2RD_B0ToKpPimMuMu_2023/DecayTree", "B0", "" },
-        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{"BuToJpsiKplus_det", "/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/BuToJpsiKplus_Detached_Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_B2CC.root", "Hlt2B2CC_BuToJpsiKplus_JpsiToMuMu_Detached/DecayTree", "Bplus", "unknown" },
+        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++)
@@ -122,34 +138,96 @@ int status_report_plots()
         data_chain->Add(ana.root_file.c_str());
 
         Double_t B_M, // 4400 -> 8200
-                Jpsi_M; // 200 -> 6600
+            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;
 
-        Float_t muplus_PX,    // 0 -> 30 * 10^3
-                muplus_PY, // 0 -> 7000
-                muplus_PZ,  // 0 -> 3
-                muplus_ENERGY,         // 0 -> 26 * 10^3
-                muminus_PX, // 0 -> 7
-                muminus_PY,        // 0 -> 11 * 10^3
-                muminus_PZ, // 0 -> 7
-                muminus_ENERGY;        // 0 -> 11 * 10^3
+        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);
+        }
 
-        // 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);
+        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.);
-        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_MIN, MASS_HIST_MAX, 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();
 
@@ -159,78 +237,245 @@ int status_report_plots()
             {
                 std::cout << "[" << ana.name << "] Processing event: " << i << " (" << TString::Format("%.2f", ((double)i / (double)entries) * 100.) << "%)" << std::endl;
             }
+            data_chain->GetEntry(i);
 
-            // TLorentzVector v4_muplus(muplus_PX, muplus_PY, muplus_PZ, muplus_ENERGY);
-            // TLorentzVector v4_muminus(muminus_PX, muminus_PY, muminus_PZ, muminus_ENERGY);
-            // Double_t mumu_inv_mass = (v4_muplus + v4_muminus).M();
+            if (ana.k_charge != 0 && Kplus_Q != ana.k_charge)
+            {
+                continue;
+            }
 
-            data_chain->GetEntry(i);
-            h1_B_M->Fill(B_M);
+            // 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);
-            h2_B_M_Jpsi_M->Fill(B_M, Jpsi_M);
-            // std::cout << mumu_inv_mass << "-" << Jpsi_M << "=" << TMath::Abs(mumu_inv_mass-Jpsi_M) << std::endl;
+            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(Jpsi_M - J_PSI_MASS) < 100.)
+                if (TMath::Abs(kpi_mass - K_STAR_0_MASS) < 100.)
                 {
-                h1_B_M_JPsi_cut->Fill(B_M);
+                    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);
 
-        auto fitRes = CreateRooFit(h1_B_M_JPsi_cut, ana);
+        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(fitRes[0], fitRes[1], ana.name.c_str(), "B_mass_JPsi_cut_fit");
+        // 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;
 }
 
-
-std::vector<RooPlot*> CreateRooFit(TH1D *hist, AnalysisOutput ana)
+void CreateRooFitAndSavePDF(TH1D *hist, AnalysisOutput ana, const char *name)
 {
-    RooRealVar roo_var_mass(ana.suf("roo_var_mass").c_str(), "B Mass Variable", MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX);
+    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("roohist_B_M");
+    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("roo_sig_bw_mean").c_str(), "Mass BW Mean", 5250., 5100., 5400.);
-    RooRealVar roo_sig_bw_with(ana.suf("roo_sig_bw_with").c_str(), "Mass BW Width", 20., 0., 50.);
+    // 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);
 
-    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_tail(ana.suf("sig_tail").c_str(), "#lambda_{sig}", -0.5, -1., 0.);
 
-    RooRealVar roo_bkg_exp_c(ana.suf("roo_bkg_exp_c").c_str(), "Background C", -0.001145, -0.00199, -0.00100);
-    RooExponential roo_bkg_exp(ana.suf("roo_bkg_exp").c_str(), "B Mass Background Exp", roo_var_mass, roo_bkg_exp_c);
+    // 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_var_mass_sig_yield(ana.suf("roo_var_mass_sig_yield").c_str(), "B Mass Sig Yield", 0., hist->GetEntries());
-    RooRealVar roo_var_mass_bkg_yield(ana.suf("roo_var_mass_bkg_yield").c_str(), "B Mass Bckg Yield", 0., hist->GetEntries());
+    // 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);
 
-    std::string pdf_name = ana.suf("roo_pdf_sig_plus_bkg");
+    // 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);
+
+    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_bw, roo_bkg_exp),
+                                   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(TString::Format("%s [%s]", ana.name.c_str(), ana.title.c_str()).Data()));
+    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::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_bw)), RooFit::FillStyle(3244), RooFit::LineColor(kRed-7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"));
+    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.40, 0.99, 0.90));
-    roo_frame_mass->getAttText()->SetTextSize(0.024);
+    // 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("Pull Distribution"));
+    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");
 
-    return std::vector<RooPlot*> { roo_frame_mass, roo_frame_pull };
+    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};
+// }