From 512ef65d1c40593b5850c8cb5e84481f8c855af7 Mon Sep 17 00:00:00 2001 From: Marius Pfeiffer Date: Wed, 28 Feb 2024 15:50:50 +0100 Subject: [PATCH] cut changes, styling --- basic_analysis.h | 112 ++++++++++-- new_analysis_b02hphmmumu.cpp | 6 +- new_analysis_b02kppimmumu.cpp | 2 +- new_analysis_bu2hpmumu.cpp | 325 +++------------------------------- 4 files changed, 128 insertions(+), 317 deletions(-) diff --git a/basic_analysis.h b/basic_analysis.h index 5aacbbd..142c716 100644 --- a/basic_analysis.h +++ b/basic_analysis.h @@ -11,6 +11,10 @@ #include "TTree.h" +#include "RooHist.h" + +#include "constants.h" + class FourVect { private: @@ -153,16 +157,57 @@ void DrawInDefaultCanvas(RooFitSummary fitSummary, const char *folder) TString name = TString::Format("%s_canvas", fitSummary.fit_histogram->GetName()); TCanvas *c = new TCanvas(name, fitSummary.fit_histogram->GetTitle(), 0, 0, 800, 600); - TPad *pull_pad = new TPad(TString::Format("%s_canv_pull", fitSummary.fit_histogram->GetName()), "Pull Pad", 0., 0., 1., 0.3); - pull_pad->Draw(); - pull_pad->SetTopMargin(0.001); - pull_pad->SetBottomMargin(0.3); - pull_pad->SetGrid(); + Double_t xlow, ylow, xup, yup; + c->GetPad(0)->GetPadPar(xlow, ylow, xup, yup); + c->Divide(1, 2); - TPad *fit_pad = new TPad(TString::Format("%s_canv_fit", fitSummary.fit_histogram->GetName()), "Fit Pad", 0., 0.3, 1., 1.); - fit_pad->Draw(); - fit_pad->SetBottomMargin(0.001); - fit_pad->cd(); + Double_t upPad_ylow = ylow + 0.30 * (yup - ylow); + Double_t dwPad_yup = ylow + 0.30 * (yup - ylow); + + TVirtualPad *upPad = c->GetPad(1); + upPad->SetPad(xlow, upPad_ylow, xup, yup); + + TVirtualPad *dwPad = c->GetPad(2); + dwPad->SetPad(xlow, ylow, xup, dwPad_yup); + + Double_t upPad_area = (yup - upPad_ylow) * (xup - xlow); + Double_t dwPad_area = (dwPad_yup - ylow) * (xup - xlow); + + std::cout << "### UP AREA: " << upPad_area << " LOW AREA: " << dwPad_area << std::endl; + + const Double_t pull_title_size = 0.09; + const Double_t pull_label_size = 0.09; + const Double_t pull_title_offset_x = 3.0; + const Double_t pull_title_offset_y = 2.4; + + Double_t fit_title_size = (pull_title_size * dwPad_area) / upPad_area; + Double_t fit_label_size = (pull_label_size * dwPad_area) / upPad_area; + Double_t fit_title_offset_x = (pull_title_offset_x * dwPad_area) / upPad_area; + Double_t fit_title_offset_y = (pull_title_offset_y * dwPad_area) / upPad_area; + + // canvas->Update(); + c->cd(1); + + // TPad *pull_pad = new TPad(TString::Format("%s_canv_pull", fitSummary.fit_histogram->GetName()), "Pull Pad", 0., 0., 1., 0.3); + // pull_pad->Draw(); + // pull_pad->SetTopMargin(0.001); + // pull_pad->SetBottomMargin(0.3); + // pull_pad->SetGrid(); + + // TPad *fit_pad = new TPad(TString::Format("%s_canv_fit", fitSummary.fit_histogram->GetName()), "Fit Pad", 0., 0.3, 1., 1.); + // fit_pad->Draw(); + // fit_pad->SetBottomMargin(0.001); + // fit_pad->cd(); + + auto fit_Xaxis = fitSummary.fit_histogram->GetXaxis(); + auto fit_Yaxis = fitSummary.fit_histogram->GetYaxis(); + fit_Xaxis->SetLabelSize(fit_label_size); + fit_Xaxis->SetTitleOffset(fit_title_offset_x); + fit_Xaxis->SetTitleSize(fit_title_size); + + fit_Yaxis->SetLabelSize(fit_label_size); + fit_Yaxis->SetTitleSize(fit_title_size); + fit_Yaxis->SetTitleOffset(fit_title_offset_y); fitSummary.fit_histogram->Draw(); @@ -185,10 +230,49 @@ void DrawInDefaultCanvas(RooFitSummary fitSummary, const char *folder) leg1->Draw(); - pull_pad->cd(); + c->cd(2); + + RooHist *pull_hist = (RooHist *)fitSummary.pull_histogram->getObject(0); + + Double_t pull_min = pull_hist->GetMinimum(); + Double_t pull_max = pull_hist->GetMaximum(); + + std::cout << "### (" << fitSummary.pull_histogram->GetName() << ") PULL MIN: " << pull_min << " PULL MAX: " << pull_max << std::endl; + + double bound = 0; + if (TMath::Abs(pull_min) > TMath::Abs(pull_max)) + { + bound = TMath::Abs(pull_min); + } + else + { + bound = TMath::Abs(pull_max); + } + + auto pull_Xaxis = fitSummary.pull_histogram->GetXaxis(); + auto pull_Yaxis = fitSummary.pull_histogram->GetYaxis(); + pull_Xaxis->SetLabelSize(pull_label_size); + pull_Xaxis->SetTitleOffset(pull_title_offset_x); + pull_Xaxis->SetTitleSize(pull_title_size); + pull_Xaxis->SetTitle(""); + + pull_Yaxis->SetLabelSize(pull_label_size); + pull_Yaxis->SetTitleSize(pull_title_size); + pull_Yaxis->SetTitleOffset(pull_title_offset_y); + pull_Yaxis->SetTitle("Pull Distribution"); + fitSummary.pull_histogram->SetTitle(""); fitSummary.pull_histogram->Draw(); + auto line_mid = new TLine(B_MASS_VAR_MIN, 0, B_MASS_VAR_MAX, 0); + line_mid->Draw(); + + auto line_up = new TLine(B_MASS_VAR_MIN, bound, B_MASS_VAR_MAX, bound); + line_up->Draw(); + + auto line_low = new TLine(B_MASS_VAR_MIN, -bound, B_MASS_VAR_MAX, -bound); + line_low->Draw(); + c->Draw(); c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data()); } @@ -313,7 +397,13 @@ RooFitSummary CreateRooFitHistogramAndFitCB(TH1D *hist, bool hasExpBkg, bool use fitted_params.push_back(FittedParam(var_mass_nR, 2)); RooPlot *roo_frame_pull = roo_var_mass.frame(RooFit::Title("Pull Distribution"), RooFit::Name(TString::Format("%s_pull_rplt", hist->GetName()))); - roo_frame_pull->addPlotable(roo_frame_mass->pullHist(hist_name, pull_compare_name), "P"); + RooHist *pull_hist = roo_frame_mass->pullHist(hist_name, pull_compare_name); + roo_frame_pull->addPlotable(pull_hist, "P"); + + Double_t pull_min = pull_hist->GetMinimum(); + Double_t pull_max = pull_hist->GetMaximum(); + + std::cout << "##### (" << pull_hist->GetName() << ") PULL MIN: " << pull_min << " PULL MAX: " << pull_max << std::endl; return RooFitSummary{ roo_frame_mass, diff --git a/new_analysis_b02hphmmumu.cpp b/new_analysis_b02hphmmumu.cpp index 5bfaa3e..1aefe2e 100644 --- a/new_analysis_b02hphmmumu.cpp +++ b/new_analysis_b02hphmmumu.cpp @@ -265,8 +265,8 @@ int new_analysis_b02hphmmumu() if (mva_response > mva_cut_value) { - if (TMath::Abs(reconstructed_Kstar.M() - KSTAR_MASS) < 100.) - { + // if (TMath::Abs(reconstructed_Kstar.M() - KSTAR_MASS) < 100.) + // { if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) { h1_B_Mass_jpsi->Fill(reconstructed_B_Mass); @@ -275,7 +275,7 @@ int new_analysis_b02hphmmumu() { h1_B_Mass_psi2s->Fill(reconstructed_B_Mass); } - } + // } } } diff --git a/new_analysis_b02kppimmumu.cpp b/new_analysis_b02kppimmumu.cpp index 3921230..c933f19 100644 --- a/new_analysis_b02kppimmumu.cpp +++ b/new_analysis_b02kppimmumu.cpp @@ -120,7 +120,7 @@ int new_analysis_b02kppimmumu() FillHlt1DecisionHistos(hlt1_decision_histos, reconstructed_B_Mass); } - if (TMath::Abs(reconstructed_Kstar.M() - KSTAR_MASS) < 100.) + if (TMath::Abs(reconstructed_Kstar.M() - KSTAR_MASS) < 50.) { if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) { diff --git a/new_analysis_bu2hpmumu.cpp b/new_analysis_bu2hpmumu.cpp index 0fd525c..7fbf69b 100644 --- a/new_analysis_bu2hpmumu.cpp +++ b/new_analysis_bu2hpmumu.cpp @@ -108,10 +108,13 @@ int new_analysis_bu2hpmumu() ConnectVarsToData(vars, data_chain, sim_chain, sig_tree, bkg_tree); - unsigned int data_entries = data_chain->GetEntries(); - unsigned int sim_entries = sim_chain->GetEntries(); + unsigned int data_entries = 100000;// data_chain->GetEntries(); + unsigned int sim_entries = 100000; //sim_chain->GetEntries(); + + std::cout << "# Got " << data_entries << " data and " << sim_entries << " simulated events." << std::endl; unsigned int bkg_events = 0; + unsigned int sig_events = 0; if (retrain_bdt) { @@ -133,9 +136,12 @@ int new_analysis_bu2hpmumu() PrintProgress(TString::Format("%s BKG Collection", analysis_name), data_entries, 10000, i); } + std::cout << "# Added " << bkg_events << " background events." << std::endl; } - else { + else + { std::cout << "# Starting without BDT retrain." << std::endl; + bkg_events = data_entries; } for (unsigned int i = 0; i < sim_entries; i++) @@ -144,15 +150,20 @@ int new_analysis_bu2hpmumu() Double_t reconstructed_B_Mass = (hp4v_sim->LorentzVector(K_MASS) + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M(); if (B_BKGCAT == 30 && TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID && TMath::Abs(Hp_TRUEID) == PID_KAON) { - if (std::all_of(vars.begin(), vars.end(), [](TV *v) + if (retrain_bdt && std::all_of(vars.begin(), vars.end(), [](TV *v) { return v->IsMCFinite(); })) { - if (retrain_bdt) - { - sig_tree->Fill(); - } - h1_B_Mass_sim->Fill(reconstructed_B_Mass); + + sig_tree->Fill(); } + + sig_events++; + h1_B_Mass_sim->Fill(reconstructed_B_Mass); + } + + if (sig_events >= bkg_events) + { + break; } PrintProgress(TString::Format("%s SIG Collection", analysis_name), sim_entries, 10000, i); @@ -224,8 +235,8 @@ int new_analysis_bu2hpmumu() std::cout << line << ": " << hits << std::endl; } - DrawInDefaultCanvas(h2_Hlt1_flags_B_Mass, analysis_name, 0.16, "COLZ"); - DrawInDefaultCanvas(h2_Hlt1_flags_excl_B_Mass, analysis_name, 0.16, "COLZ"); + // DrawInDefaultCanvas(h2_Hlt1_flags_B_Mass, analysis_name, 0.16, "COLZ"); + // DrawInDefaultCanvas(h2_Hlt1_flags_excl_B_Mass, analysis_name, 0.16, "COLZ"); DrawInDefaultCanvas(h1_B_Mass_jpsi, analysis_name, 0.1); DrawInDefaultCanvas(h1_B_Mass_psi2s, analysis_name, 0.1); @@ -237,7 +248,7 @@ int new_analysis_bu2hpmumu() DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name); DrawInDefaultCanvas(roofit_hist_sim, analysis_name); - DrawHlt1DecisionHistos(analysis_name, hlt1_decision_histos); + // DrawHlt1DecisionHistos(analysis_name, hlt1_decision_histos); DrawBDTProbs(h1_bdt_probs, mva_cut_value, analysis_name); @@ -332,293 +343,3 @@ int analyze_bu2hpmumu_simulation() return 0; } - -// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% B0 To Hp Hm Mu Mu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -int analyze_b02hphmmumu() -{ - const char *analysis_name = "B0ToHpHmMuMu"; - const char *end_state_mass_literal = "m(#pi^{+}_{(#rightarrow K^{+})} #pi^{-} #mu^{+}#mu^{-})"; - - TChain *data_chain = new TChain(TString::Format("%s/DecayTree", analysis_name)); - data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0_tuple_90000000_v0r0p6288631.root"); - data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0r1_tuple_90000000_2023_v0r0p6288631.root"); - - Double_t Hp_PID_K, Hm_PID_K, Hp_PID_PI, Hm_PID_PI; - - data_chain->SetBranchAddress("Hp_PID_K", &Hp_PID_K); - data_chain->SetBranchAddress("Hm_PID_K", &Hm_PID_K); - data_chain->SetBranchAddress("Hp_PID_PI", &Hp_PID_PI); - data_chain->SetBranchAddress("Hm_PID_PI", &Hm_PID_PI); - - FourVect *l14v = FourVect::Init(data_chain, "L1"); - FourVect *l24v = FourVect::Init(data_chain, "L2"); - FourVect *hp4v = FourVect::Init(data_chain, "Hp"); - FourVect *hm4v = FourVect::Init(data_chain, "Hm"); - - TH1D *h1_B_Mass_jpsi = new TH1D("h1_B_Mass_jpsi", TString::Format("B Mass, J/#psi Mode (%s)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); - TH1D *h1_B_Mass_psi2s = new TH1D("h1_B_Mass_psi2s", TString::Format("B Mass, #psi(2S) Mode (%s)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); - TH2D *h2_Hlt1_flags_B_Mass = new TH2D("h2_Hlt1_flags_B_Mass", "Hlt1 Decision vs B Mass", 50, 5100, 5400, 13, 1., 14.); - TH2D *h2_Hlt1_flags_excl_B_Mass = new TH2D("h2_Hlt1_flags_excl_B_Mass", "Excl Hlt1 Decision vs B Mass", 50, 5100, 5400, 13, 1., 14.); - - h1_B_Mass_jpsi->GetXaxis()->SetTitle(end_state_mass_literal); - h1_B_Mass_psi2s->GetXaxis()->SetTitle(end_state_mass_literal); - h2_Hlt1_flags_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); - h2_Hlt1_flags_excl_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); - - ConnectHlt1Decisions(data_chain, h2_Hlt1_flags_B_Mass, h2_Hlt1_flags_excl_B_Mass); - - auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name); - - std::map exclusive_hits{}; - - unsigned int entries = data_chain->GetEntries(); - for (unsigned int i = 0; i < entries; i++) - { - data_chain->GetEntry(i); - - TLorentzVector dimuon = l14v->LorentzVector() + l24v->LorentzVector(); - - TLorentzVector lv_hp{}; - TLorentzVector lv_hm{}; - - if (Hp_PID_K > 0 && Hm_PID_K < 0) - { - continue; - lv_hp = hp4v->LorentzVector(K_MASS); - lv_hm = hm4v->LorentzVector(); - } - else if (Hp_PID_K < 0 && Hm_PID_K > 0) - { - lv_hp = hp4v->LorentzVector(); - lv_hm = hm4v->LorentzVector(K_MASS); - } - else - { - continue; - } - - Double_t reconstructed_Kstar_Mass = (lv_hp + lv_hm).M(); - Double_t reconstructed_B_Mass = (lv_hp + lv_hm + dimuon).M(); - - if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) - { - CheckHlt1Decisioins(h2_Hlt1_flags_B_Mass, h2_Hlt1_flags_excl_B_Mass, exclusive_hits, reconstructed_B_Mass); - FillHlt1DecisionHistos(hlt1_decision_histos, reconstructed_B_Mass); - } - - if (CutHlt1DecisionsOr({"TwoTrackMuonMVA", "DiMuonHighMass"})) - { - if (TMath::Abs(reconstructed_Kstar_Mass - KSTAR_MASS) < 50.) - { - if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) - { - h1_B_Mass_jpsi->Fill(reconstructed_B_Mass); - } - else if (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.) - { - h1_B_Mass_psi2s->Fill(reconstructed_B_Mass); - } - } - } - - PrintProgress(analysis_name, entries, 10000, i); - } - - std::cout << "# Exclusive Hits" << std::endl; - for (const auto &[line, hits] : exclusive_hits) - { - std::cout << line << ": " << hits << std::endl; - } - - DrawInDefaultCanvas(h2_Hlt1_flags_B_Mass, analysis_name, 0.16, "COLZ"); - DrawInDefaultCanvas(h2_Hlt1_flags_excl_B_Mass, analysis_name, 0.16, "COLZ"); - DrawInDefaultCanvas(h1_B_Mass_jpsi, analysis_name, 0.1); - DrawInDefaultCanvas(h1_B_Mass_psi2s, analysis_name, 0.1); - - auto roofit_hist_jpsi = CreateRooFitHistogram(h1_B_Mass_jpsi); - auto roofit_hist_psi2s = CreateRooFitHistogram(h1_B_Mass_psi2s); - DrawInDefaultCanvas(roofit_hist_jpsi, analysis_name, 0.1); - DrawInDefaultCanvas(roofit_hist_psi2s, analysis_name, 0.1); - - DrawHlt1DecisionHistos(analysis_name, hlt1_decision_histos); - - return 0; -} - -// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Hlt2RD_BuToKpMuMu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -int analyze_Hlt2RD_BuToKpMuMu() -{ - const char *analysis_name = "Hlt2RD_BuToKpMuMu"; - const char *end_state_mass_literal = "m(K^{+} #mu^{+}#mu^{-})"; - - TChain *data_chain = new TChain(TString::Format("%s/DecayTree", analysis_name)); - data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_RD.root"); - - FourVect *l14v = FourVect::Init(data_chain, "muplus"); - FourVect *l24v = FourVect::Init(data_chain, "muminus"); - FourVect *hp4v = FourVect::Init(data_chain, "Kplus"); - - TH1D *h1_B_Mass_jpsi = new TH1D("h1_B_Mass_jpsi", TString::Format("B Mass, J/#psi Mode (%s)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); - TH1D *h1_B_Mass_psi2s = new TH1D("h1_B_Mass_psi2s", TString::Format("B Mass, #psi(2S) Mode (%s)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); - TH2D *h2_Hlt1_flags_B_Mass = new TH2D("h2_Hlt1_flags_B_Mass", "Hlt1 Decision vs B Mass", 50, 5100, 5400, 13, 1., 14.); - TH2D *h2_Hlt1_flags_excl_B_Mass = new TH2D("h2_Hlt1_flags_excl_B_Mass", "Excl Hlt1 Decision vs B Mass", 50, 5100, 5400, 13, 1., 14.); - - h1_B_Mass_jpsi->GetXaxis()->SetTitle(end_state_mass_literal); - h1_B_Mass_psi2s->GetXaxis()->SetTitle(end_state_mass_literal); - h2_Hlt1_flags_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); - h2_Hlt1_flags_excl_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); - - ConnectHlt1Decisions(data_chain, h2_Hlt1_flags_B_Mass, h2_Hlt1_flags_excl_B_Mass); - - auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name); - - std::map exclusive_hits{}; - - unsigned int entries = data_chain->GetEntries(); - for (unsigned int i = 0; i < entries; i++) - { - data_chain->GetEntry(i); - - TLorentzVector dimuon = l14v->LorentzVector() + l24v->LorentzVector(); - - TLorentzVector lv_hp = hp4v->LorentzVector(); - - Double_t reconstructed_B_Mass = (lv_hp + dimuon).M(); - - if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) - { - CheckHlt1Decisioins(h2_Hlt1_flags_B_Mass, h2_Hlt1_flags_excl_B_Mass, exclusive_hits, reconstructed_B_Mass); - FillHlt1DecisionHistos(hlt1_decision_histos, reconstructed_B_Mass); - } - - if (CutHlt1DecisionsOr({"TwoTrackMuonMVA", "DiMuonHighMass"})) - { - if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) - { - h1_B_Mass_jpsi->Fill(reconstructed_B_Mass); - } - else if (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.) - { - h1_B_Mass_psi2s->Fill(reconstructed_B_Mass); - } - } - - PrintProgress("Hlt2RD_BuToKpMuMu", entries, 10000, i); - } - - std::cout << "# Exclusive Hits" << std::endl; - for (const auto &[line, hits] : exclusive_hits) - { - std::cout << line << ": " << hits << std::endl; - } - - DrawInDefaultCanvas(h2_Hlt1_flags_B_Mass, analysis_name, 0.16, "COLZ"); - DrawInDefaultCanvas(h2_Hlt1_flags_excl_B_Mass, analysis_name, 0.16, "COLZ"); - DrawInDefaultCanvas(h1_B_Mass_jpsi, analysis_name, 0.1); - DrawInDefaultCanvas(h1_B_Mass_psi2s, analysis_name, 0.1); - - auto roofit_hist_jpsi = CreateRooFitHistogram(h1_B_Mass_jpsi); - auto roofit_hist_psi2s = CreateRooFitHistogram(h1_B_Mass_psi2s); - DrawInDefaultCanvas(roofit_hist_jpsi, analysis_name, 0.1); - DrawInDefaultCanvas(roofit_hist_psi2s, analysis_name, 0.1); - - DrawHlt1DecisionHistos(analysis_name, hlt1_decision_histos); - - return 0; -} - -// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Hlt2RD_B0ToKpPimMuMu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -int analyze_Hlt2RD_B0ToKpPimMuMu() -{ - const char *analysis_name = "Hlt2RD_B0ToKpPimMuMu"; - const char *end_state_mass_literal = "m(K^{+} #pi^{-} #mu^{+}#mu^{-})"; - - TChain *data_chain = new TChain(TString::Format("%s/DecayTree", analysis_name)); - data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_RD.root"); - - FourVect *l14v = FourVect::Init(data_chain, "muplus"); - FourVect *l24v = FourVect::Init(data_chain, "muminus"); - FourVect *hp4v = FourVect::Init(data_chain, "Kplus"); - FourVect *hm4v = FourVect::Init(data_chain, "piminus"); - - TH1D *h1_B_Mass_jpsi = new TH1D("h1_B_Mass_jpsi", TString::Format("B Mass, J/#psi Mode (%s)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); - TH1D *h1_B_Mass_psi2s = new TH1D("h1_B_Mass_psi2s", TString::Format("B Mass, #psi(2S) Mode (%s)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); - TH2D *h2_Hlt1_flags_B_Mass = new TH2D("h2_Hlt1_flags_B_Mass", "Hlt1 Decision vs B Mass", 50, 5100, 5400, 13, 1., 14.); - TH2D *h2_Hlt1_flags_excl_B_Mass = new TH2D("h2_Hlt1_flags_excl_B_Mass", "Excl Hlt1 Decision vs B Mass", 50, 5100, 5400, 13, 1., 14.); - - h1_B_Mass_jpsi->GetXaxis()->SetTitle(end_state_mass_literal); - h1_B_Mass_psi2s->GetXaxis()->SetTitle(end_state_mass_literal); - h2_Hlt1_flags_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); - h2_Hlt1_flags_excl_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); - - ConnectHlt1Decisions(data_chain, h2_Hlt1_flags_B_Mass, h2_Hlt1_flags_excl_B_Mass); - - auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name); - - std::map exclusive_hits{}; - - unsigned int entries = data_chain->GetEntries(); - for (unsigned int i = 0; i < entries; i++) - { - data_chain->GetEntry(i); - - TLorentzVector dimuon = l14v->LorentzVector() + l24v->LorentzVector(); - - TLorentzVector lv_hp = hp4v->LorentzVector(); - TLorentzVector lv_hm = hm4v->LorentzVector(); - - Double_t reconstructed_Kstar_Mass = (lv_hp + lv_hm).M(); - Double_t reconstructed_B_Mass = (lv_hp + lv_hm + dimuon).M(); - - if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) - { - CheckHlt1Decisioins(h2_Hlt1_flags_B_Mass, h2_Hlt1_flags_excl_B_Mass, exclusive_hits, reconstructed_B_Mass); - FillHlt1DecisionHistos(hlt1_decision_histos, reconstructed_B_Mass); - } - - if (CutHlt1DecisionsOr({"TwoTrackMuonMVA", "DiMuonHighMass"})) - { - if (TMath::Abs(reconstructed_Kstar_Mass - KSTAR_MASS) < 50.) - { - if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) - { - h1_B_Mass_jpsi->Fill(reconstructed_B_Mass); - } - else if (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.) - { - h1_B_Mass_psi2s->Fill(reconstructed_B_Mass); - } - } - } - - PrintProgress(analysis_name, entries, 10000, i); - } - - std::cout << "# Exclusive Hits" << std::endl; - for (const auto &[line, hits] : exclusive_hits) - { - std::cout << line << ": " << hits << std::endl; - } - - DrawInDefaultCanvas(h2_Hlt1_flags_B_Mass, analysis_name, 0.16, "COLZ"); - DrawInDefaultCanvas(h2_Hlt1_flags_excl_B_Mass, analysis_name, 0.16, "COLZ"); - DrawInDefaultCanvas(h1_B_Mass_jpsi, analysis_name, 0.1); - DrawInDefaultCanvas(h1_B_Mass_psi2s, analysis_name, 0.1); - - auto roofit_hist_jpsi = CreateRooFitHistogram(h1_B_Mass_jpsi); - auto roofit_hist_psi2s = CreateRooFitHistogram(h1_B_Mass_psi2s); - DrawInDefaultCanvas(roofit_hist_jpsi, analysis_name, 0.1); - DrawInDefaultCanvas(roofit_hist_psi2s, analysis_name, 0.1); - - DrawHlt1DecisionHistos(analysis_name, hlt1_decision_histos); - - return 0; -} - -// int new_analysis() -// { -// return analyze_bu2hpmumu_data(); -// } \ No newline at end of file