From 1a2e5569c1cb83ce43a365c19098eb36e65259c4 Mon Sep 17 00:00:00 2001 From: Marius Pfeiffer Date: Fri, 26 Jan 2024 14:16:59 +0100 Subject: [PATCH] macros for data and simulation evaluation, start of do_all --- analysis_fullstream_b02hphmmumu.cpp | 71 +++++-- analysis_fullstream_bu2hpmumu.cpp | 71 ++++++- analysis_turbo.cpp | 67 +++++- do_all.cpp | 233 +++++++++++++++++++++ simulation.cpp | 135 ++++++++++--- simulation_fullstream_B02KpPimMuMu.cpp | 269 +++++++++++++++++++++++++ simulation_fullstream_Bu2KpMuMu.cpp | 243 ++++++++++++++++++++++ status_report_plots.cpp | 3 + 8 files changed, 1034 insertions(+), 58 deletions(-) create mode 100644 do_all.cpp create mode 100644 simulation_fullstream_B02KpPimMuMu.cpp create mode 100644 simulation_fullstream_Bu2KpMuMu.cpp diff --git a/analysis_fullstream_b02hphmmumu.cpp b/analysis_fullstream_b02hphmmumu.cpp index a10174f..9ca0418 100644 --- a/analysis_fullstream_b02hphmmumu.cpp +++ b/analysis_fullstream_b02hphmmumu.cpp @@ -40,9 +40,9 @@ 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_Pim2Km"; -const char* MASS_LITERAL = "m(#pi^{+}#pi^{-}_{(#rightarrow K^{-})}#mu^{+}#mu^{-})"; +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; @@ -173,14 +173,14 @@ int analysis_fullstream_b02hphmmumu() 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); + 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); + // 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); @@ -198,9 +198,9 @@ int analysis_fullstream_b02hphmmumu() // (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) & + (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) + (Hp_BPVIPCHI2 > 6) & (Hp_PT > 250) & (Hp_P > 2000) & (Hp_PID_K > 0) ) & // (Hlt2RD_B0ToKpPimMuMuDecision) & ((Hlt2_InclDetDiMuon_4BodyDecision) | (Hlt2_InclDetDiMuon_3BodyDecision) | (Hlt2_InclDetDiMuonDecision))) @@ -247,18 +247,63 @@ void CreateRooFitAndDraw(TH1D *hist, int fitting_entries) roohist_B0_M.plotOn(roo_frame_mass, RooFit::Binning(nBins), RooFit::Name("B Mass Distribution")); roo_frame_mass->GetXaxis()->SetTitle(MASS_LITERAL); + // Crystal Ball for Signal + RooRealVar var_mass_x0("var_mass_x0", "#mu", 5278., 5170., 5500.); + RooRealVar var_mass_sigmaLR("var_mass_sigmaLR", "#sigma_{LR}", 16., 5., 40.); + + // Same Variables for Left and Right Tail + // RooRealVar var_mass_alphaL("var_mass_alphaL", "#alpha_{L}", 2., 0., 4.); + // RooRealVar var_mass_nL("var_mass_nL", "n_{L}", 5., 0., 15.); + + // RooRealVar var_mass_alphaR("var_mass_alphaR", "#alpha_{R}", 2., 0., 4.); + // RooRealVar var_mass_nR("var_mass_nR", "n_{R}", 5., 0., 15.); + + auto var_mass_alphaL = RooRealConstant::value(1.915); + auto var_mass_nL = RooRealConstant::value(0.573); + auto var_mass_alphaR = RooRealConstant::value(2.257); + auto var_mass_nR = RooRealConstant::value(0.418); + + RooCrystalBall sig_cb("sig_cb", "Signal Crystal Ball", roo_var_mass, var_mass_x0, var_mass_sigmaLR, var_mass_alphaL, var_mass_nL, var_mass_alphaR, var_mass_nR); + + // Exponential for Background + RooRealVar var_mass_bkg_c("var_mass_bkg_c", "#lambda", -0.0014, -0.004, -0.000); + RooExponential bkg_exp("bkg_exp", "Exp Background", roo_var_mass, var_mass_bkg_c); + + RooRealVar var_mass_nsig("nsig", "Mass N Signal", 0., hist->GetEntries()); + RooRealVar var_mass_nbkg("nbkg", "Mass N Background", 0., hist->GetEntries()); + + // TString pdf_name = TString::Format("%s_sigplusbkg", var_id.c_str()); + RooAddPdf sigplusbkg("sigplusbkg", "Sig and Bkg PDF", RooArgList(sig_cb, bkg_exp), RooArgList(var_mass_nsig, var_mass_nbkg)); + + RooFitResult *fitres = sigplusbkg.fitTo(roohist_B0_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range("fitting_range")); + + sigplusbkg.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144)); + sigplusbkg.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range("fitting_range"), RooFit::Name("sigplusbkg")); + sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(bkg_exp)), RooFit::LineColor(kBlue-7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("bkg_exp")); + sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(sig_cb)), RooFit::LineColor(kRed-7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("sig_cb")); + 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); + TLegend *leg1 = new TLegend(0.50, 0.58, 0.87, 0.89); leg1->SetFillColor(kWhite); leg1->SetLineColor(kBlack); + leg1->AddEntry(roo_frame_mass->findObject("sigplusbkg"), "Signal + Background", "LP"); + leg1->AddEntry(roo_frame_mass->findObject("sig_cb"), "Signal", "LP"); + leg1->AddEntry(roo_frame_mass->findObject("bkg_exp"), "Background", "LP"); + leg1->AddEntry((TObject *)0, "", ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_x0.getTitle().Data(), var_mass_x0.getVal(), var_mass_x0.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_sigmaLR.getTitle().Data(), var_mass_sigmaLR.getVal(), var_mass_sigmaLR.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.6f #pm %.6f", var_mass_bkg_c.getTitle().Data(), var_mass_bkg_c.getVal(), var_mass_bkg_c.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "S", var_mass_nsig.getVal(), var_mass_nsig.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "B", var_mass_nbkg.getVal(), var_mass_nbkg.getError()).Data(), ""); + 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()); + c->SaveAs(TString::Format("images/data/roofit_hist_%s_fitted_bmass.pdf", FILE_NAME).Data()); } \ No newline at end of file diff --git a/analysis_fullstream_bu2hpmumu.cpp b/analysis_fullstream_bu2hpmumu.cpp index 489f4c0..e85e4ec 100644 --- a/analysis_fullstream_bu2hpmumu.cpp +++ b/analysis_fullstream_bu2hpmumu.cpp @@ -43,9 +43,9 @@ 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^{-})"; +const char *TITLE = "SpruceRD_BuToHpMuMu (#pi^{+} #rightarrow K^{+}, Small Cuts 2)"; +const char *FILE_NAME = "SpruceRD_BuToHpMuMu_Pip2Kp_small2"; +const char *MASS_LITERAL = "m(#pi^{+}_{(#rightarrow K^{+})}#mu^{+}#mu^{-})"; int analysis_fullstream_bu2hpmumu() { @@ -96,7 +96,7 @@ int analysis_fullstream_bu2hpmumu() 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); @@ -145,15 +145,22 @@ int analysis_fullstream_bu2hpmumu() 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))) + + if ((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) { + if (MASS_HIST_FIT_MIN <= reconstructed_B_Mass && reconstructed_B_Mass <= MASS_HIST_FIT_MAX) + { fitting_entries++; } } + // 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))))) + // { + + // } + if (i % 10000 == 0) { std::cout << "[" @@ -183,18 +190,64 @@ void CreateRooFitAndDraw(TH1D *hist, int fitting_entries) roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(nBins), RooFit::Name("B Mass Distribution")); roo_frame_mass->GetXaxis()->SetTitle(MASS_LITERAL); + + // Crystal Ball for Signal + RooRealVar var_mass_x0("var_mass_x0", "#mu", 5278., 5170., 5500.); + RooRealVar var_mass_sigmaLR("var_mass_sigmaLR", "#sigma_{LR}", 16., 5., 40.); + + // Same Variables for Left and Right Tail + // RooRealVar var_mass_alphaL("var_mass_alphaL", "#alpha_{L}", 2., 0., 4.); + // RooRealVar var_mass_nL("var_mass_nL", "n_{L}", 5., 0., 15.); + + // RooRealVar var_mass_alphaR("var_mass_alphaR", "#alpha_{R}", 2., 0., 4.); + // RooRealVar var_mass_nR("var_mass_nR", "n_{R}", 5., 0., 15.); + + auto var_mass_alphaL = RooRealConstant::value(1.309); + auto var_mass_nL = RooRealConstant::value(0.72); + auto var_mass_alphaR = RooRealConstant::value(1.226); + auto var_mass_nR = RooRealConstant::value(0.992); + + RooCrystalBall sig_cb("sig_cb", "Signal Crystal Ball", roo_var_mass, var_mass_x0, var_mass_sigmaLR, var_mass_alphaL, var_mass_nL, var_mass_alphaR, var_mass_nR); + + // Exponential for Background + RooRealVar var_mass_bkg_c("var_mass_bkg_c", "#lambda", -0.0014, -0.004, -0.000); + RooExponential bkg_exp("bkg_exp", "Exp Background", roo_var_mass, var_mass_bkg_c); + + RooRealVar var_mass_nsig("nsig", "Mass N Signal", 0., hist->GetEntries()); + RooRealVar var_mass_nbkg("nbkg", "Mass N Background", 0., hist->GetEntries()); + + // TString pdf_name = TString::Format("%s_sigplusbkg", var_id.c_str()); + RooAddPdf sigplusbkg("sigplusbkg", "Sig and Bkg PDF", RooArgList(sig_cb, bkg_exp), RooArgList(var_mass_nsig, var_mass_nbkg)); + + RooFitResult *fitres = sigplusbkg.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range("fitting_range")); + + sigplusbkg.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144)); + sigplusbkg.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range("fitting_range"), RooFit::Name("sigplusbkg")); + sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(bkg_exp)), RooFit::LineColor(kBlue-7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("bkg_exp")); + sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(sig_cb)), RooFit::LineColor(kRed-7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("sig_cb")); + 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); + TLegend *leg1 = new TLegend(0.50, 0.58, 0.87, 0.89); leg1->SetFillColor(kWhite); leg1->SetLineColor(kBlack); - leg1->AddEntry((TObject*)0, TString::Format("Entries: %d", fitting_entries), ""); + leg1->AddEntry(roo_frame_mass->findObject("sigplusbkg"), "Signal + Background", "LP"); + leg1->AddEntry(roo_frame_mass->findObject("sig_cb"), "Signal", "LP"); + leg1->AddEntry(roo_frame_mass->findObject("bkg_exp"), "Background", "LP"); + leg1->AddEntry((TObject *)0, "", ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_x0.getTitle().Data(), var_mass_x0.getVal(), var_mass_x0.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_sigmaLR.getTitle().Data(), var_mass_sigmaLR.getVal(), var_mass_sigmaLR.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.6f #pm %.6f", var_mass_bkg_c.getTitle().Data(), var_mass_bkg_c.getVal(), var_mass_bkg_c.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "S", var_mass_nsig.getVal(), var_mass_nsig.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "B", var_mass_nbkg.getVal(), var_mass_nbkg.getError()).Data(), ""); + + 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()); + c->SaveAs(TString::Format("images/data/roofit_hist_%s_bmass.pdf", FILE_NAME).Data()); } \ No newline at end of file diff --git a/analysis_turbo.cpp b/analysis_turbo.cpp index 46f8638..482186d 100644 --- a/analysis_turbo.cpp +++ b/analysis_turbo.cpp @@ -34,8 +34,8 @@ const int nBins = 70; -const Double_t MASS_HIST_FIT_MIN = 5100.; -const Double_t MASS_HIST_FIT_MAX = 6000.; +const Double_t MASS_HIST_MIN = 5100.; +const Double_t MASS_HIST_MAX = 6000.; const char *B0_DECAY = "Hlt2RD_B0ToKpPimMuMu"; const char *Bu_DECAY = "Hlt2RD_BuToKpMuMu"; @@ -72,7 +72,7 @@ int analysis_turbo() data_chain->SetBranchAddress("Hlt1TwoTrackMVADecision", &Hlt1TwoTrackMVADecision); } - TH1D *h1_B_M = new TH1D("h1_B_M", TITLE, nBins, MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX); + TH1D *h1_B_M = new TH1D("h1_B_M", TITLE, nBins, MASS_HIST_MIN, MASS_HIST_MAX); h1_B_M->GetXaxis()->SetTitle(X_AXIS); int fitting_entries = 0; @@ -86,7 +86,7 @@ int analysis_turbo() 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); - if (MASS_HIST_FIT_MIN <= B0_M && B0_M <= MASS_HIST_FIT_MAX) + if (MASS_HIST_MIN <= B0_M && B0_M <= MASS_HIST_MAX) { fitting_entries++; } @@ -97,7 +97,7 @@ int analysis_turbo() if ((TMath::Abs(Jpsi_M - 3096.9) < 100.) && (B_M > 4500.) && (B_M < 6000.) && ((Hlt1TrackMVADecision) || (Hlt1TwoTrackMVADecision))) { h1_B_M->Fill(B_M); - if (MASS_HIST_FIT_MIN <= B_M && B_M <= MASS_HIST_FIT_MAX) + if (MASS_HIST_MIN <= B_M && B_M <= MASS_HIST_MAX) { fitting_entries++; } @@ -112,6 +112,10 @@ int analysis_turbo() } } + std::cout << "[" + << TITLE + << "] Processed event: " << entries << " (" << TString::Format("%.2f", 100.) << "%)" << std::endl; + TCanvas *c1 = new TCanvas("c1", "c1", 0, 0, 800, 600); h1_B_M->Draw(); c1->Draw(); @@ -124,8 +128,8 @@ int analysis_turbo() 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); + RooRealVar roo_var_mass("var_mass", "B Mass Variable", MASS_HIST_MIN, MASS_HIST_MAX); + roo_var_mass.setRange("fitting_range", MASS_HIST_MIN, MASS_HIST_MAX); RooDataHist roohist_B_M("roohist_B_M", "B Mass Histogram", roo_var_mass, RooFit::Import(*hist)); @@ -133,18 +137,63 @@ void CreateRooFitAndDraw(TH1D *hist, int fitting_entries) roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(nBins), RooFit::Name("B Mass Distribution")); roo_frame_mass->GetXaxis()->SetTitle(X_AXIS); + // Crystal Ball for Signal + RooRealVar var_mass_x0("var_mass_x0", "#mu", 5278., 5170., 5500.); + RooRealVar var_mass_sigmaLR("var_mass_sigmaLR", "#sigma_{LR}", 16., 5., 40.); + + // Same Variables for Left and Right Tail + // RooRealVar var_mass_alphaL("var_mass_alphaL", "#alpha_{L}", 2., 0., 4.); + // RooRealVar var_mass_nL("var_mass_nL", "n_{L}", 5., 0., 15.); + + // RooRealVar var_mass_alphaR("var_mass_alphaR", "#alpha_{R}", 2., 0., 4.); + // RooRealVar var_mass_nR("var_mass_nR", "n_{R}", 5., 0., 15.); + + auto var_mass_alphaL = RooRealConstant::value(1.948); + auto var_mass_nL = RooRealConstant::value(0.618); + auto var_mass_alphaR = RooRealConstant::value(2.320); + auto var_mass_nR = RooRealConstant::value(0.473); + + RooCrystalBall sig_cb("sig_cb", "Signal Crystal Ball", roo_var_mass, var_mass_x0, var_mass_sigmaLR, var_mass_alphaL, var_mass_nL, var_mass_alphaR, var_mass_nR); + + // Exponential for Background + RooRealVar var_mass_bkg_c("var_mass_bkg_c", "#lambda", -0.0014, -0.004, -0.000); + RooExponential bkg_exp("bkg_exp", "Exp Background", roo_var_mass, var_mass_bkg_c); + + RooRealVar var_mass_nsig("nsig", "Mass N Signal", 0., hist->GetEntries()); + RooRealVar var_mass_nbkg("nbkg", "Mass N Background", 0., hist->GetEntries()); + + // TString pdf_name = TString::Format("%s_sigplusbkg", var_id.c_str()); + RooAddPdf sigplusbkg("sigplusbkg", "Sig and Bkg PDF", RooArgList(sig_cb, bkg_exp), RooArgList(var_mass_nsig, var_mass_nbkg)); + + RooFitResult *fitres = sigplusbkg.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range("fitting_range")); + + sigplusbkg.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144)); + sigplusbkg.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range("fitting_range"), RooFit::Name("sigplusbkg")); + sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(bkg_exp)), RooFit::LineColor(kBlue-7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("bkg_exp")); + sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(sig_cb)), RooFit::LineColor(kRed-7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("sig_cb")); + 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); + TLegend *leg1 = new TLegend(0.50, 0.58, 0.87, 0.89); leg1->SetFillColor(kWhite); leg1->SetLineColor(kBlack); + leg1->AddEntry(roo_frame_mass->findObject("sigplusbkg"), "Signal + Background", "LP"); + leg1->AddEntry(roo_frame_mass->findObject("sig_cb"), "Signal", "LP"); + leg1->AddEntry(roo_frame_mass->findObject("bkg_exp"), "Background", "LP"); + leg1->AddEntry((TObject *)0, "", ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_x0.getTitle().Data(), var_mass_x0.getVal(), var_mass_x0.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_sigmaLR.getTitle().Data(), var_mass_sigmaLR.getVal(), var_mass_sigmaLR.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.6f #pm %.6f", var_mass_bkg_c.getTitle().Data(), var_mass_bkg_c.getVal(), var_mass_bkg_c.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "S", var_mass_nsig.getVal(), var_mass_nsig.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "B", var_mass_nbkg.getVal(), var_mass_nbkg.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("Entries: %d", fitting_entries), ""); leg1->Draw(); c->Draw(); - c->SaveAs(TString::Format("images/roofit_hist_%s_bmass.pdf", TITLE).Data()); + c->SaveAs(TString::Format("images/data/roofit_hist_%s_fitted_bmass.pdf", TITLE).Data()); } \ No newline at end of file diff --git a/do_all.cpp b/do_all.cpp new file mode 100644 index 0000000..d252209 --- /dev/null +++ b/do_all.cpp @@ -0,0 +1,233 @@ +#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 "TString.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 + +template +class Cut { + public: + virtual bool Check(T y) = 0; + virtual std::string ToString() = 0; +}; + +template +class LT : public Cut { + private: + T value; + + public: + LT(T val) : value{val} { + + } + + bool Check(T y) { + return y < value; + } + + std::string ToString() { + return TString::Format("(x < %s)", std::to_string(value).c_str()).Data(); + } +}; + +template +class GT : public Cut { + private: + T value; + + public: + GT(T val) : value{val} {} + + bool Check(T y) { + return y > value; + } + + std::string ToString() { + return TString::Format("(x > %s)", std::to_string(value).c_str()).Data(); + } +}; + +template +class EQ : public Cut { + private: + T value; + + public: + EQ(T val) : value{val} {} + + bool Check(T y) { + return y == value; + } + + std::string ToString() { + return TString::Format("(x == %s)", std::to_string(value).c_str()).Data(); + } +}; + +template +class BT : public Cut { + private: + T low; + T up; + + public: + BT(T lower, T upper) : low{lower}, up{upper} {} + + bool Check(T y) { + return low < y && y < up; + } + + std::string ToString() { + return TString::Format("(x in [%s, %s])", std::to_string(low).c_str(), std::to_string(up).c_str()).Data(); + } +}; + +template +class Variable { + private: + const char* name; + std::vector*> cuts; + + T value; + + public: + + Variable() { + } + + Variable(const char* name, std::vector*> cuts) { + this->name = name; + this->cuts = cuts; + } + + T* GetValuePtr() { + return &value; + } + + T GetValue() const { + return value; + } +}; + +class VariableCollection { + private: + std::map> doubles; + std::map> floats; + std::map> ints; + std::map> bools; + + public: + VariableCollection() { + doubles = std::map>{}; + floats = std::map>{}; + ints = std::map>{}; + bools = std::map>{}; + } + + void Connect(TChain *chain) { + for (auto& [key, var] : doubles) { + chain->SetBranchAddress(key.c_str(), var.GetValuePtr()); + } + + for (auto& [key, var] : floats) { + chain->SetBranchAddress(key.c_str(), var.GetValuePtr()); + } + + for (auto& [key, var] : ints) { + chain->SetBranchAddress(key.c_str(), var.GetValuePtr()); + } + + for (auto& [key, var] : bools) { + chain->SetBranchAddress(key.c_str(), var.GetValuePtr()); + } + } + + void AddDouble(const char* name, const std::vector*>& cuts) { + doubles[name] = Variable(name, cuts); + } + + void AddFloat(const char* name, const std::vector*>& cuts) { + floats[name] = Variable(name, cuts); + } + + void AddInt(const char* name, const std::vector*>& cuts) { + ints[name] = Variable(name, cuts); + } + + void AddBool(const char* name, const std::vector*>& cuts) { + bools[name] = Variable(name, cuts); + } + + void PrintValues() { + for (const auto& [key, var] : doubles) { + std::cout << key << ": " << var.GetValue() << std::endl; + } + + for (const auto& [key, var] : floats) { + std::cout << key << ": " << var.GetValue() << std::endl; + } + + for (const auto& [key, var] : ints) { + std::cout << key << ": " << var.GetValue() << std::endl; + } + + for (const auto& [key, var] : bools) { + std::cout << key << ": " << var.GetValue() << std::endl; + } + } +}; + +int do_all() { + + VariableCollection col{}; + col.AddDouble("B_M", {}); + col.AddDouble("Jpsi_M", {}); + col.AddFloat("Kplus_PT", {}); + col.AddBool("muminus_ISMUON", {}); + col.AddInt("muplus_Q", {}); + + TChain *data_chain = new TChain(TString::Format("%s/DecayTree", "Hlt2RD_BuToKpMuMu").Data()); + data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_RD.root"); + + col.Connect(data_chain); + + data_chain->GetEntry(1); + + col.PrintValues(); + + data_chain->GetEntry(2); + + col.PrintValues(); + + return 0; +} \ No newline at end of file diff --git a/simulation.cpp b/simulation.cpp index 41dae4c..2166ddb 100644 --- a/simulation.cpp +++ b/simulation.cpp @@ -33,30 +33,32 @@ #include #include -const int nBins = 170;//70; +const int nBins = 70; -const Double_t MASS_HIST_MIN = 2000.; +const Double_t MASS_HIST_MIN = 5100.; const Double_t MASS_HIST_MAX = 6000.; -const Double_t MASS_HIST_FIT_MIN = 2000.; -const Double_t MASS_HIST_FIT_MAX = 6000.; const Double_t K_MASS = 493.677; +const char *NAME = "B0ToKpPimMuMu_Trubo"; +const char *X_AXIS = "m(K^{+}#pi^{-}#mu^{+}#mu^{-})"; +const bool USE_HYP_REPLACE = false; +const char *B_NAME = "B0"; + +TString quickformat_name(const char *format) +{ + return TString::Format(format, NAME); +}; + +void CreateRooFitAndDraw(TH1D *hist); + int simulation() { - const char *name = "BuToHpMuMu_Fullstream"; - const bool use_hyp_replace = true; - const char *b_name = "B"; + TChain *data_chain = new TChain("B0ToKpPimMuMu_noPID/DecayTree"); + // data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/rd_btoxll_simulation_fullstream_v0r0p6671378_BuToKpMuMu_12143001_magdown.root"); + // data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/rd_btoxll_simulation_turbo_v0r0p6657752_BuToKpMuMu_12143001_magdown.root"); + data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/rd_btoxll_simulation_turbo_v0r0p6657752_B0ToKpPimMuMu_11144002_magdown.root"); - auto quickformat_name = [name](const char *format) - { - return TString::Format(format, name); - }; - - TChain *data_chain = new TChain("BuToHpMuMu_noPID/DecayTree"); - data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/rd_btoxll_simulation_fullstream_v0r0p6671378_BuToKpMuMu_12143001_magdown.root"); - - // TRUE ID: Muon sind Muons, Hadronen sind K, pi /* rd_btoxll_simulation_fullstream_v0r0p6671378_B0ToKpPimMuMu_11144002_magdown.root rd_btoxll_simulation_fullstream_v0r0p6671378_BuToKpMuMu_12143001_magdown.root @@ -65,20 +67,20 @@ int simulation() */ Int_t B_BKGCAT; - data_chain->SetBranchAddress(TString::Format("%s_BKGCAT", b_name).Data(), &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; Double_t B_M; - Int_t L1_TRUEID, L2_TRUEID, Hp_TRUEID; + Int_t L1_TRUEID, L2_TRUEID, Hp_TRUEID, Hm_TRUEID; data_chain->SetBranchAddress("L1_TRUEID", &L1_TRUEID); data_chain->SetBranchAddress("L2_TRUEID", &L2_TRUEID); - data_chain->SetBranchAddress("Hp_TRUEID", &Hp_TRUEID); + - if (use_hyp_replace) + if (USE_HYP_REPLACE) { data_chain->SetBranchAddress("L1_PX", &L1_PX); data_chain->SetBranchAddress("L1_PY", &L1_PY); @@ -92,10 +94,14 @@ int simulation() data_chain->SetBranchAddress("Hp_PY", &Hp_PY); data_chain->SetBranchAddress("Hp_PZ", &Hp_PZ); data_chain->SetBranchAddress("Hp_ENERGY", &Hp_ENERGY); + + data_chain->SetBranchAddress("Hp_TRUEID", &Hp_TRUEID); } else { - data_chain->SetBranchAddress(TString::Format("%s_M", b_name).Data(), &B_M); + data_chain->SetBranchAddress(TString::Format("%s_M", B_NAME).Data(), &B_M); + data_chain->SetBranchAddress("K_TRUEID", &Hp_TRUEID); + // data_chain->SetBranchAddress("Pi_TRUEID", &Hm_TRUEID); } TH1D *h1_B_M = new TH1D("h1_B_M", "B Mass", nBins, MASS_HIST_MIN, MASS_HIST_MAX); @@ -103,6 +109,9 @@ int simulation() 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); + h1_B_M->GetXaxis()->SetTitle(X_AXIS); + h1_B_M_trueid->GetXaxis()->SetTitle(X_AXIS); + unsigned int entries = data_chain->GetEntries(); for (size_t i = 0; i < entries; i++) { @@ -111,7 +120,7 @@ int simulation() Double_t reconstructed_B_Mass = 0; - if (use_hyp_replace) + if (USE_HYP_REPLACE) { TVector3 K_momentum(Hp_PX, Hp_PY, Hp_PZ); @@ -126,10 +135,10 @@ int simulation() reconstructed_B_Mass = B_M; } - if (TMath::Abs(L1_TRUEID) == 13 && L2_TRUEID == -L1_TRUEID && TMath::Abs(Hp_TRUEID) == 321) + // std::cout << "L1: " << L1_TRUEID << ", L2: " << L2_TRUEID << ", Hp: " << Hp_TRUEID << std::endl; + if (TMath::Abs(L1_TRUEID) == 13 && L2_TRUEID == -L1_TRUEID && TMath::Abs(Hp_TRUEID) == 321/* && TMath::Abs(Hm_TRUEID) == 211*/) { - std::cout << "L1: " << L1_TRUEID << ", L2: " << L2_TRUEID << ", Hp: " << Hp_TRUEID << std::endl; - std::cout << "Mass: " << reconstructed_B_Mass << std::endl; + // std::cout << "Mass: " << reconstructed_B_Mass << std::endl; h1_B_M_trueid->Fill(reconstructed_B_Mass); } @@ -140,11 +149,15 @@ int simulation() if (i % 10000 == 0) { std::cout << "[" - << name + << NAME << "] Processed event: " << i << " (" << TString::Format("%.2f", ((double)i / (double)entries) * 100.) << "%)" << std::endl; } } + std::cout << "[" + << NAME + << "] Processed event: " << entries << " (" << TString::Format("%.2f", 100.) << "%)" << std::endl; + h1_B_M->SetLineColor(kBlue); h1_B_M_trueid->SetLineColor(kMagenta + 2); h1_B_M_trueid->SetFillColor(kMagenta + 2); @@ -154,7 +167,7 @@ int simulation() h1_B_BKGCAT->SetFillColor(kBlue); h1_B_BKGCAT->SetFillStyle(3351); - std::filesystem::create_directory(TString::Format("images/sim/%s", name).Data()); + 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); @@ -178,5 +191,73 @@ int simulation() c3->Draw(); c3->SaveAs(quickformat_name("images/sim/%s/h2_B_M_vs_B_BKGCAT.pdf")); + CreateRooFitAndDraw(h1_B_M_trueid); + return 0; +} + +void CreateRooFitAndDraw(TH1D *hist) +{ + // std::cout << "### Entries: " << hist->GetEntries() << std::endl; + RooRealVar roo_var_mass("var_mass", "B Mass Variable", MASS_HIST_MIN, MASS_HIST_MAX); + roo_var_mass.setRange("fitting_range", MASS_HIST_MIN, MASS_HIST_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(quickformat_name("Simulation of %s"))); + roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(nBins), RooFit::Name("B Mass Distribution")); + roo_frame_mass->GetXaxis()->SetTitle(X_AXIS); + + // Crystal Ball for Signal + RooRealVar var_mass_x0("var_mass_x0", "#mu", 5278., 5170., 5500.); + RooRealVar var_mass_sigmaLR("var_mass_sigmaLR", "#sigma_{LR}", 16., 5., 25.); + + // Same Variables for Left and Right Tail + RooRealVar var_mass_alphaL("var_mass_alphaL", "#alpha_{L}", 2., 0., 4.); + RooRealVar var_mass_nL("var_mass_nL", "n_{L}", 5., 0., 15.); + + RooRealVar var_mass_alphaR("var_mass_alphaR", "#alpha_{R}", 2., 0., 4.); + RooRealVar var_mass_nR("var_mass_nR", "n_{R}", 5., 0., 15.); + + RooCrystalBall sig_cb("sig_cb", "Signal Crystal Ball", roo_var_mass, var_mass_x0, var_mass_sigmaLR, var_mass_alphaL, var_mass_nL, var_mass_alphaR, var_mass_nR); + + // Exponential for Background + // RooRealVar var_mass_bkg_c("var_mass_bkg_c", "Background C", -0.00231049, -0.003, -0.002); + // RooExponential bkg_exp("bkg_exp", "Exp Background", var_mass, var_mass_bkg_c); + + // RooRealVar var_mass_nsig(TString::Format("var_%s_nsig", var_id.c_str()), "Mass N Signal", 0., hist->GetEntries()); + // RooRealVar var_mass_nbkg(TString::Format("var_%s_nbkg", var_id.c_str()), "Mass N Background", 0., hist->GetEntries()); + + // TString pdf_name = TString::Format("%s_sigplusbkg", var_id.c_str()); + // RooAddPdf sigplusbkg(pdf_name, "Sig and Bkg PDF", RooArgList(sig_cb, bkg_exp), RooArgList(var_mass_nsig, var_mass_nbkg)); + + RooFitResult *fitres = sig_cb.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range("fitting_range")); + + sig_cb.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144)); + sig_cb.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range("fitting_range"), RooFit::Name("sig_cb")); + // sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(bkg_exp)), RooFit::LineColor(kBlue), RooFit::LineStyle(kDashed), RooFit::Range(5170., 5500.)); + + TCanvas *c = new TCanvas("roofit_c", "roofit_c", 0, 0, 800, 600); + + roo_frame_mass->Draw(); + + TLegend *leg1 = new TLegend(0.58, 0.50, 0.87, 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("sig_cb"), "Signal", "LP"); + //leg1->AddEntry(roo_frame_mass->findObject(name_fit_func_bkg.c_str()), "Background", "LP"); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_x0.getTitle().Data(), var_mass_x0.getVal(), var_mass_x0.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_sigmaLR.getTitle().Data(), var_mass_sigmaLR.getVal(), var_mass_sigmaLR.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 = %.3f #pm %.3f", var_mass_alphaL.getTitle().Data(), var_mass_alphaL.getVal(), var_mass_alphaL.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_nL.getTitle().Data(), var_mass_nL.getVal(), var_mass_nL.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_alphaR.getTitle().Data(), var_mass_alphaR.getVal(), var_mass_alphaR.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_nR.getTitle().Data(), var_mass_nR.getVal(), var_mass_nR.getError()).Data(), ""); + + leg1->Draw(); + + c->Draw(); + c->SaveAs(quickformat_name("images/sim/%s/roof_B_M.pdf")); } \ No newline at end of file diff --git a/simulation_fullstream_B02KpPimMuMu.cpp b/simulation_fullstream_B02KpPimMuMu.cpp new file mode 100644 index 0000000..8688b38 --- /dev/null +++ b/simulation_fullstream_B02KpPimMuMu.cpp @@ -0,0 +1,269 @@ +#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 K_MASS = 493.677; + +const int PID_KAON = 321; +const int PID_PION = 211; +const int PID_MUON = 13; + +const char *NAME = "B0ToHpHmMuMu_Fullstream"; +const char *X_AXIS = "m(#pi^{+}_{#rightarrow K^{+}}#pi^{-}#mu^{+}#mu^{-})"; +const bool USE_HYP_REPLACE = false; +const char *B_NAME = "B0"; + +TString quickformat_name(const char *format) +{ + return TString::Format(format, NAME); +}; + +void CreateRooFitAndDraw(TH1D *hist); + +int simulation_fullstream_B02KpPimMuMu() +{ + TChain *data_chain = new TChain("B0ToHpHmMuMu_noPID/DecayTree"); + data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/rd_btoxll_simulation_fullstream_v0r0p6671378_B0ToKpPimMuMu_11144002_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; + + Int_t L1_TRUEID, L2_TRUEID, Hp_TRUEID, Hm_TRUEID; + + data_chain->SetBranchAddress("L1_TRUEID", &L1_TRUEID); + data_chain->SetBranchAddress("L2_TRUEID", &L2_TRUEID); + + 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); + + data_chain->SetBranchAddress("Hp_TRUEID", &Hp_TRUEID); + data_chain->SetBranchAddress("Hm_TRUEID", &Hm_TRUEID); + + TH1D *h1_B_M = new TH1D("h1_B_M", "B Mass", nBins, MASS_HIST_MIN, MASS_HIST_MAX); + TH1D *h1_B_M_trueid = new TH1D("h1_B_M_trueid", "B Mass, TrueID matched", 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); + + h1_B_M->GetXaxis()->SetTitle(X_AXIS); + h1_B_M_trueid->GetXaxis()->SetTitle(X_AXIS); + + 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 (TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID) + { + if (TMath::Abs(Hp_TRUEID) == PID_KAON && TMath::Abs(Hm_TRUEID) == PID_PION) + { + 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 Pi_4v(Hm_PX, Hm_PY, Hm_PZ, Hm_ENERGY); + reconstructed_B_Mass = (K_4v + l1_4v + l2_4v + Pi_4v).M(); + } + else if (TMath::Abs(Hp_TRUEID) == PID_PION && TMath::Abs(Hm_TRUEID) == PID_KAON) + { + TVector3 K_momentum(Hm_PX, Hm_PY, Hm_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 Pi_4v(Hp_PX, Hp_PY, Hp_PZ, Hp_ENERGY); + reconstructed_B_Mass = (K_4v + l1_4v + l2_4v + Pi_4v).M(); + } + } + + if (reconstructed_B_Mass != 0) + { + h1_B_M_trueid->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; + } + } + + std::cout << "[" + << NAME + << "] Processed event: " << entries << " (" << TString::Format("%.2f", 100.) << "%)" << std::endl; + + h1_B_M->SetLineColor(kBlue); + h1_B_M_trueid->SetLineColor(kMagenta + 2); + h1_B_M_trueid->SetFillColor(kMagenta + 2); + h1_B_M_trueid->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_trueid->SetStats(0); + h1_B_M->Draw(); + h1_B_M_trueid->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")); + + CreateRooFitAndDraw(h1_B_M_trueid); + + return 0; +} + +void CreateRooFitAndDraw(TH1D *hist) +{ + // std::cout << "### Entries: " << hist->GetEntries() << std::endl; + RooRealVar roo_var_mass("var_mass", "B Mass Variable", MASS_HIST_MIN, MASS_HIST_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(quickformat_name("Simulation of %s"))); + roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(nBins), RooFit::Name("B Mass Distribution")); + roo_frame_mass->GetXaxis()->SetTitle(X_AXIS); + + // Crystal Ball for Signal + RooRealVar var_mass_x0("var_mass_x0", "#mu", 5278., 5170., 5500.); + RooRealVar var_mass_sigmaLR("var_mass_sigmaLR", "#sigma_{LR}", 16., 5., 25.); + + RooRealVar var_mass_alphaL("var_mass_alphaL", "#alpha_{L}", 2., 0., 4.); + RooRealVar var_mass_nL("var_mass_nL", "n_{L}", 5., 0., 15.); + + RooRealVar var_mass_alphaR("var_mass_alphaR", "#alpha_{R}", 2., 0., 4.); + RooRealVar var_mass_nR("var_mass_nR", "n_{R}", 5., 0., 15.); + + RooCrystalBall sig_cb("sig_cb", "Signal Crystal Ball", roo_var_mass, var_mass_x0, var_mass_sigmaLR, var_mass_alphaL, var_mass_nL, var_mass_alphaR, var_mass_nR); + + // Exponential for Background + // RooRealVar var_mass_bkg_c("var_mass_bkg_c", "Background C", -0.00231049, -0.003, -0.002); + // RooExponential bkg_exp("bkg_exp", "Exp Background", var_mass, var_mass_bkg_c); + + // RooRealVar var_mass_nsig(TString::Format("var_%s_nsig", var_id.c_str()), "Mass N Signal", 0., hist->GetEntries()); + // RooRealVar var_mass_nbkg(TString::Format("var_%s_nbkg", var_id.c_str()), "Mass N Background", 0., hist->GetEntries()); + + // TString pdf_name = TString::Format("%s_sigplusbkg", var_id.c_str()); + // RooAddPdf sigplusbkg(pdf_name, "Sig and Bkg PDF", RooArgList(sig_cb, bkg_exp), RooArgList(var_mass_nsig, var_mass_nbkg)); + + RooFitResult *fitres = sig_cb.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range("fitting_range")); + + sig_cb.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144)); + sig_cb.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range("fitting_range"), RooFit::Name("sig_cb")); + // sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(bkg_exp)), RooFit::LineColor(kBlue), RooFit::LineStyle(kDashed), RooFit::Range(5170., 5500.)); + + TCanvas *c = new TCanvas("roofit_c", "roofit_c", 0, 0, 800, 600); + + roo_frame_mass->Draw(); + + TLegend *leg1 = new TLegend(0.58, 0.50, 0.87, 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("sig_cb"), "Signal", "LP"); + // leg1->AddEntry(roo_frame_mass->findObject(name_fit_func_bkg.c_str()), "Background", "LP"); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_x0.getTitle().Data(), var_mass_x0.getVal(), var_mass_x0.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_sigmaLR.getTitle().Data(), var_mass_sigmaLR.getVal(), var_mass_sigmaLR.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 = %.3f #pm %.3f", var_mass_alphaL.getTitle().Data(), var_mass_alphaL.getVal(), var_mass_alphaL.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_nL.getTitle().Data(), var_mass_nL.getVal(), var_mass_nL.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_alphaR.getTitle().Data(), var_mass_alphaR.getVal(), var_mass_alphaR.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_nR.getTitle().Data(), var_mass_nR.getVal(), var_mass_nR.getError()).Data(), ""); + + leg1->Draw(); + + c->Draw(); + c->SaveAs(quickformat_name("images/sim/%s/roof_B_M.pdf")); +} \ No newline at end of file diff --git a/simulation_fullstream_Bu2KpMuMu.cpp b/simulation_fullstream_Bu2KpMuMu.cpp new file mode 100644 index 0000000..267eb6b --- /dev/null +++ b/simulation_fullstream_Bu2KpMuMu.cpp @@ -0,0 +1,243 @@ +#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 K_MASS = 493.677; + +const int PID_KAON = 321; +const int PID_PION = 211; +const int PID_MUON = 13; + +const char *NAME = "BuToHpMuMu_Fullstream_NoKTrueID"; +const char *X_AXIS = "m(#pi^{+}_{#rightarrow K^{+}}#mu^{+}#mu^{-})"; +const char *B_NAME = "B"; + +TString quickformat_name(const char *format) +{ + return TString::Format(format, NAME); +}; + +void CreateRooFitAndDraw(TH1D *hist); + +int simulation_fullstream_Bu2KpMuMu() +{ + TChain *data_chain = new TChain("BuToHpMuMu_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; + Double_t B_M; + + Int_t L1_TRUEID, L2_TRUEID, Hp_TRUEID; + + data_chain->SetBranchAddress("L1_TRUEID", &L1_TRUEID); + data_chain->SetBranchAddress("L2_TRUEID", &L2_TRUEID); + + 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("Hp_TRUEID", &Hp_TRUEID); + + TH1D *h1_B_M = new TH1D("h1_B_M", "B Mass", nBins, MASS_HIST_MIN, MASS_HIST_MAX); + TH1D *h1_B_M_trueid = new TH1D("h1_B_M_trueid", "B Mass, TrueID matched", 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); + + h1_B_M->GetXaxis()->SetTitle(X_AXIS); + h1_B_M_trueid->GetXaxis()->SetTitle(X_AXIS); + + 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 (TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID/* && TMath::Abs(Hp_TRUEID) == PID_KAON*/) + { + 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); + reconstructed_B_Mass = (K_4v + l1_4v + l2_4v).M(); + + h1_B_M_trueid->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; + } + } + + std::cout << "[" + << NAME + << "] Processed event: " << entries << " (" << TString::Format("%.2f", 100.) << "%)" << std::endl; + + h1_B_M->SetLineColor(kBlue); + h1_B_M_trueid->SetLineColor(kMagenta + 2); + h1_B_M_trueid->SetFillColor(kMagenta + 2); + h1_B_M_trueid->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_trueid->SetStats(0); + h1_B_M->Draw(); + h1_B_M_trueid->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")); + + CreateRooFitAndDraw(h1_B_M_trueid); + + return 0; +} + +void CreateRooFitAndDraw(TH1D *hist) +{ + // std::cout << "### Entries: " << hist->GetEntries() << std::endl; + RooRealVar roo_var_mass("var_mass", "B Mass Variable", MASS_HIST_MIN, MASS_HIST_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(quickformat_name("Simulation of %s"))); + roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(nBins), RooFit::Name("B Mass Distribution")); + roo_frame_mass->GetXaxis()->SetTitle(X_AXIS); + + // Crystal Ball for Signal + RooRealVar var_mass_x0("var_mass_x0", "#mu", 5278., 5170., 5500.); + RooRealVar var_mass_sigmaLR("var_mass_sigmaLR", "#sigma_{LR}", 16., 5., 25.); + + RooRealVar var_mass_alphaL("var_mass_alphaL", "#alpha_{L}", 2., 0., 4.); + RooRealVar var_mass_nL("var_mass_nL", "n_{L}", 5., 0., 15.); + + RooRealVar var_mass_alphaR("var_mass_alphaR", "#alpha_{R}", 2., 0., 4.); + RooRealVar var_mass_nR("var_mass_nR", "n_{R}", 5., 0., 15.); + + RooCrystalBall sig_cb("sig_cb", "Signal Crystal Ball", roo_var_mass, var_mass_x0, var_mass_sigmaLR, var_mass_alphaL, var_mass_nL, var_mass_alphaR, var_mass_nR); + + // Exponential for Background + // RooRealVar var_mass_bkg_c("var_mass_bkg_c", "Background C", -0.00231049, -0.003, -0.002); + // RooExponential bkg_exp("bkg_exp", "Exp Background", var_mass, var_mass_bkg_c); + + // RooRealVar var_mass_nsig(TString::Format("var_%s_nsig", var_id.c_str()), "Mass N Signal", 0., hist->GetEntries()); + // RooRealVar var_mass_nbkg(TString::Format("var_%s_nbkg", var_id.c_str()), "Mass N Background", 0., hist->GetEntries()); + + // TString pdf_name = TString::Format("%s_sigplusbkg", var_id.c_str()); + // RooAddPdf sigplusbkg(pdf_name, "Sig and Bkg PDF", RooArgList(sig_cb, bkg_exp), RooArgList(var_mass_nsig, var_mass_nbkg)); + + RooFitResult *fitres = sig_cb.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range("fitting_range")); + + sig_cb.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144)); + sig_cb.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range("fitting_range"), RooFit::Name("sig_cb")); + // sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(bkg_exp)), RooFit::LineColor(kBlue), RooFit::LineStyle(kDashed), RooFit::Range(5170., 5500.)); + + TCanvas *c = new TCanvas("roofit_c", "roofit_c", 0, 0, 800, 600); + + roo_frame_mass->Draw(); + + TLegend *leg1 = new TLegend(0.58, 0.50, 0.87, 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("sig_cb"), "Signal", "LP"); + // leg1->AddEntry(roo_frame_mass->findObject(name_fit_func_bkg.c_str()), "Background", "LP"); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_x0.getTitle().Data(), var_mass_x0.getVal(), var_mass_x0.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_sigmaLR.getTitle().Data(), var_mass_sigmaLR.getVal(), var_mass_sigmaLR.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 = %.3f #pm %.3f", var_mass_alphaL.getTitle().Data(), var_mass_alphaL.getVal(), var_mass_alphaL.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_nL.getTitle().Data(), var_mass_nL.getVal(), var_mass_nL.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_alphaR.getTitle().Data(), var_mass_alphaR.getVal(), var_mass_alphaR.getError()).Data(), ""); + leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_nR.getTitle().Data(), var_mass_nR.getVal(), var_mass_nR.getError()).Data(), ""); + + leg1->Draw(); + + c->Draw(); + c->SaveAs(quickformat_name("images/sim/%s/roof_B_M.pdf")); +} \ No newline at end of file diff --git a/status_report_plots.cpp b/status_report_plots.cpp index 915a091..791698f 100644 --- a/status_report_plots.cpp +++ b/status_report_plots.cpp @@ -3,6 +3,7 @@ #include "THStack.h" #include "TGraph.h" #include "TTree.h" +#include "TF1.h" #include "TChain.h" #include "TFile.h" #include "TCanvas.h" @@ -376,6 +377,8 @@ void CreateRooFitAndSavePDF(TH1D *hist, AnalysisOutput ana, const char *name) RooRealVar roo_bkg_exp_c(ana.suf("bkg_exp_c").c_str(), "#lambda_{bkg}", -0.001145, -0.00199, -0.00100); RooExponential roo_bkg_exp(ana.suf("bkg_exp").c_str(), "B Mass Background Exp", roo_var_mass, roo_bkg_exp_c); + roo_bkg_exp.asTF() + RooRealVar roo_var_mass_sig_yield(ana.suf("sig_yield").c_str(), "N_{Sig}", 0., hist->GetEntries()); RooRealVar roo_var_mass_bkg_yield(ana.suf("bkg_yield").c_str(), "N_{Bkg}", 0., hist->GetEntries());