diff --git a/.gitignore b/.gitignore index 6e9431b..d7f4219 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ *.root dataloader/** status_report/** -branchnames.txt \ No newline at end of file +branchnames.txt +images/** \ No newline at end of file diff --git a/analysis_fullstream_b02hphmmumu.cpp b/analysis_fullstream_b02hphmmumu.cpp new file mode 100644 index 0000000..fbbe8d7 --- /dev/null +++ b/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 +#include +#include + +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()); +} \ No newline at end of file diff --git a/analysis_fullstream_bu2hpmumu.cpp b/analysis_fullstream_bu2hpmumu.cpp new file mode 100644 index 0000000..ea1e5d5 --- /dev/null +++ b/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 +#include +#include + +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()); +} \ No newline at end of file diff --git a/analysis_fullstream_michele.cpp b/analysis_fullstream_michele.cpp new file mode 100644 index 0000000..47b425d --- /dev/null +++ b/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 +#include +#include + +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()); +} \ No newline at end of file diff --git a/analysis_turbo.cpp b/analysis_turbo.cpp new file mode 100644 index 0000000..2ea4f9e --- /dev/null +++ b/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 +#include +#include + +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()); +} \ No newline at end of file diff --git a/simulation.cpp b/simulation.cpp new file mode 100644 index 0000000..8ed2575 --- /dev/null +++ b/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 +#include +#include +#include + +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; +} \ No newline at end of file diff --git a/status_report_plots.cpp b/status_report_plots.cpp index c2054e4..915a091 100644 --- a/status_report_plots.cpp +++ b/status_report_plots.cpp @@ -24,94 +24,110 @@ #include "RooProduct.h" #include "RooCrystalBall.h" #include "RooBreitWigner.h" +#include "RooArgSet.h" +#include "RooFFTConvPdf.h" +#include "RooNovosibirsk.h" #include #include #include 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 CreateRooFit(TH1D *hist, AnalysisOutput ana); +void CreateRooFitAndSavePDF(TH1D *hist, AnalysisOutput ana, const char *name); +std::vector 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() { // gROOT->ProcessLine(".x /work/pfeiffer/lhcbStyle.C"); - std::vector 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" }, + std::vector 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++) @@ -121,35 +137,97 @@ int status_report_plots() 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 + 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; - 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); + + if (ana.k_charge != 0 && Kplus_Q != ana.k_charge) + { + continue; + } - // 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 (!(muplus_ISMUON && muminus_ISMUON && (Hlt1TrackMVADecision | Hlt1TwoTrackMVADecision))) + // { + // continue; + // } - data_chain->GetEntry(i); - h1_B_M->Fill(B_M); + 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 (TMath::Abs(Jpsi_M - J_PSI_MASS) < 100.) + if (ana.B_name == "B0") { - h1_B_M_JPsi_cut->Fill(B_M); + 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_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); - 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"); + // 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; } - -std::vector 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); - 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.); - 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); + 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_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_tail(ana.suf("sig_tail").c_str(), "#lambda_{sig}", -0.5, -1., 0.); - std::string pdf_name = ana.suf("roo_pdf_sig_plus_bkg"); + // 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); + + 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 { roo_frame_mass, roo_frame_pull }; -} \ No newline at end of file + 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 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{roo_frame_mass}; +// } \ No newline at end of file