From aeba3d41cdc2516337beaa3b0466dc939c80aa8b Mon Sep 17 00:00:00 2001 From: Marius Pfeiffer Date: Thu, 29 Feb 2024 15:54:28 +0100 Subject: [PATCH] plot styling --- basic_analysis.h | 65 ++++++++++++++++++++++++----------- new_analysis_b02hphmmumu.cpp | 28 +++++++++------ new_analysis_b02kppimmumu.cpp | 57 +++++++++++++++++++++++------- new_analysis_bu2hpmumu.cpp | 24 ++++++------- new_analysis_bu2kpmumu.cpp | 11 +++--- 5 files changed, 123 insertions(+), 62 deletions(-) diff --git a/basic_analysis.h b/basic_analysis.h index 142c716..df8d5a4 100644 --- a/basic_analysis.h +++ b/basic_analysis.h @@ -151,6 +151,23 @@ void DrawInDefaultCanvas(RooPlot *rooPlot, const char *folder, double margin_lef c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data()); } +void PrintStyles(RooPlot* plt) { + auto Xaxis = plt->GetXaxis(); + auto Yaxis = plt->GetYaxis(); + std::cout << "## Styles for plot " << plt->GetName() << std::endl; + std::cout << "## X axis\n" + << "Label Size: " << Xaxis->GetLabelSize() << ",\n" + << "Label Offset: " << Xaxis->GetLabelOffset() << ",\n" + << "Title Size: " << Xaxis->GetTitleSize() << ",\n" + << "Title Offset: " << Xaxis->GetTitleOffset() << "." << std::endl; + + std::cout << "## Y axis\n" + << "Label Size: " << Yaxis->GetLabelSize() << ",\n" + << "Label Offset: " << Yaxis->GetLabelOffset() << ",\n" + << "Title Size: " << Yaxis->GetTitleSize() << ",\n" + << "Title Offset: " << Yaxis->GetTitleOffset() << "." << std::endl; +} + void DrawInDefaultCanvas(RooFitSummary fitSummary, const char *folder) { std::filesystem::create_directory(TString::Format("output_files/analysis/%s", folder).Data()); @@ -161,29 +178,26 @@ void DrawInDefaultCanvas(RooFitSummary fitSummary, const char *folder) c->GetPad(0)->GetPadPar(xlow, ylow, xup, yup); c->Divide(1, 2); - Double_t upPad_ylow = ylow + 0.30 * (yup - ylow); - Double_t dwPad_yup = ylow + 0.30 * (yup - ylow); + Double_t upPad_ylow = ylow + 0.25 * (yup - ylow); + Double_t dwPad_yup = ylow + 0.25 * (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); + dwPad->SetTopMargin(0); 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; + // 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); @@ -201,13 +215,15 @@ void DrawInDefaultCanvas(RooFitSummary fitSummary, const char *folder) auto fit_Xaxis = fitSummary.fit_histogram->GetXaxis(); auto fit_Yaxis = fitSummary.fit_histogram->GetYaxis(); + fit_Xaxis->SetLabelOffset(0.02); fit_Xaxis->SetLabelSize(fit_label_size); - fit_Xaxis->SetTitleOffset(fit_title_offset_x); + fit_Xaxis->SetTitleOffset(1.4); fit_Xaxis->SetTitleSize(fit_title_size); + fit_Yaxis->SetLabelOffset(0.01); fit_Yaxis->SetLabelSize(fit_label_size); fit_Yaxis->SetTitleSize(fit_title_size); - fit_Yaxis->SetTitleOffset(fit_title_offset_y); + fit_Yaxis->SetTitleOffset(0); fitSummary.fit_histogram->Draw(); @@ -232,12 +248,10 @@ void DrawInDefaultCanvas(RooFitSummary fitSummary, const char *folder) c->cd(2); - RooHist *pull_hist = (RooHist *)fitSummary.pull_histogram->getObject(0); + Double_t pull_min = fitSummary.pull_histogram->GetMinimum(); + Double_t pull_max = fitSummary.pull_histogram->GetMaximum(); - 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; + // 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)) @@ -251,26 +265,35 @@ void DrawInDefaultCanvas(RooFitSummary fitSummary, const char *folder) auto pull_Xaxis = fitSummary.pull_histogram->GetXaxis(); auto pull_Yaxis = fitSummary.pull_histogram->GetYaxis(); + pull_Xaxis->SetLabelOffset(0.02); pull_Xaxis->SetLabelSize(pull_label_size); - pull_Xaxis->SetTitleOffset(pull_title_offset_x); pull_Xaxis->SetTitleSize(pull_title_size); pull_Xaxis->SetTitle(""); + pull_Yaxis->SetLabelOffset(0.01); pull_Yaxis->SetLabelSize(pull_label_size); pull_Yaxis->SetTitleSize(pull_title_size); - pull_Yaxis->SetTitleOffset(pull_title_offset_y); + pull_Yaxis->SetTitleOffset(0.45); pull_Yaxis->SetTitle("Pull Distribution"); fitSummary.pull_histogram->SetTitle(""); + fitSummary.pull_histogram->SetMaximum(bound); + fitSummary.pull_histogram->SetMinimum(-bound); + fitSummary.pull_histogram->Draw(); auto line_mid = new TLine(B_MASS_VAR_MIN, 0, B_MASS_VAR_MAX, 0); + line_mid->SetLineColor(kOrange+4); line_mid->Draw(); - auto line_up = new TLine(B_MASS_VAR_MIN, bound, B_MASS_VAR_MAX, bound); + auto line_up = new TLine(B_MASS_VAR_MIN, bound / 2, B_MASS_VAR_MAX, bound / 2); + line_up->SetLineStyle(kDashed); + line_up->SetLineColor(kOrange+4); line_up->Draw(); - auto line_low = new TLine(B_MASS_VAR_MIN, -bound, B_MASS_VAR_MAX, -bound); + auto line_low = new TLine(B_MASS_VAR_MIN, -bound / 2, B_MASS_VAR_MAX, -bound / 2); + line_low->SetLineStyle(kDashed); + line_low->SetLineColor(kOrange+4); line_low->Draw(); c->Draw(); @@ -400,10 +423,10 @@ RooFitSummary CreateRooFitHistogramAndFitCB(TH1D *hist, bool hasExpBkg, bool use 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(); + // 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; + // 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 1aefe2e..74811b8 100644 --- a/new_analysis_b02hphmmumu.cpp +++ b/new_analysis_b02hphmmumu.cpp @@ -47,8 +47,8 @@ int new_analysis_b02hphmmumu() const char *analysis_name = "B0ToHpHmMuMu"; const char *data_tree_name = "B0ToHpHmMuMu"; const char *sim_tree_name = "B0ToHpHmMuMu_noPID"; - const char *end_state_mass_literal = "m(#pi^{+}_{(#rightarrow K^{+})}#pi^{-}_{(#rightarrow K^{-})}#mu^{+}#mu^{-})"; - const bool retrain_bdt = true; + const char *end_state_mass_literal = "m(#pi^{+}#pi^{-}_{(#rightarrow K^{-})}#mu^{+}#mu^{-} & #pi^{+}_{(#rightarrow K^{+})}#pi^{-}#mu^{+}#mu^{-})"; + const bool retrain_bdt = false; TChain *data_chain = new TChain(TString::Format("%s/DecayTree", data_tree_name)); // data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0_tuple_90000000_v0r0p6288631.root"); @@ -88,6 +88,7 @@ int new_analysis_b02hphmmumu() h1_B_Mass_jpsi->GetXaxis()->SetTitle(end_state_mass_literal); h1_B_Mass_psi2s->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_sim->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); @@ -119,9 +120,10 @@ int new_analysis_b02hphmmumu() ConnectVarsToData(vars, data_chain, sim_chain, sig_tree, bkg_tree); unsigned int data_entries = data_chain->GetEntries(); - unsigned int sim_entries = 100000; // sim_chain->GetEntries(); + unsigned int sim_entries = sim_chain->GetEntries(); unsigned int bkg_events = 0; + unsigned int sig_events = 0; if (retrain_bdt) { @@ -159,10 +161,12 @@ int new_analysis_b02hphmmumu() PrintProgress(TString::Format("%s BKG Collection", analysis_name), data_entries, 10000, i); } + std::cout << "# Added " << bkg_events << " background events." << std::endl; } else { std::cout << "# Starting without BDT retrain." << std::endl; + bkg_events = data_entries; } for (unsigned int i = 0; i < sim_entries; i++) @@ -188,15 +192,17 @@ int new_analysis_b02hphmmumu() { Double_t reconstructed_B_Mass = (reconstructed_Kstar + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M(); - if (std::all_of(vars.begin(), vars.end(), [](TV *v) - { return v->IsMCFinite(); })) + if (sig_events < bkg_events) { - if (retrain_bdt) + if (retrain_bdt && std::all_of(vars.begin(), vars.end(), [](TV *v) + { return v->IsMCFinite(); })) { sig_tree->Fill(); + sig_events++; } - h1_B_Mass_sim->Fill(reconstructed_B_Mass); } + + h1_B_Mass_sim->Fill(reconstructed_B_Mass); } } @@ -214,7 +220,7 @@ int new_analysis_b02hphmmumu() Float_t *train_vars = new Float_t[vars.size()]; auto reader = SetupReader(vars, train_vars, analysis_name); - const double mva_cut_value = -0.12; + const double mva_cut_value = -0.0508; for (unsigned int i = 0; i < data_entries; i++) { @@ -265,8 +271,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) < 50.) + { if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) { h1_B_Mass_jpsi->Fill(reconstructed_B_Mass); @@ -275,7 +281,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 c933f19..f2b0c0c 100644 --- a/new_analysis_b02kppimmumu.cpp +++ b/new_analysis_b02kppimmumu.cpp @@ -53,10 +53,14 @@ int new_analysis_b02kppimmumu() TChain *data_chain = new TChain(TString::Format("%s/DecayTree", data_tree_name)); data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_RD.root"); + Int_t Hp_Q, Hm_Q; + FourVect *l14v_data = FourVect::Init(data_chain, "muplus"); FourVect *l24v_data = FourVect::Init(data_chain, "muminus"); FourVect *hp4v_data = FourVect::Init(data_chain, "Kplus"); FourVect *hm4v_data = FourVect::Init(data_chain, "piminus"); + data_chain->SetBranchAddress("Kplus_Q", &Hp_Q); + data_chain->SetBranchAddress("piminus_Q", &Hm_Q); TChain *sim_chain = new TChain(TString::Format("%s/DecayTree", sim_tree_name)); sim_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/rd_btoxll_simulation_turbo_v0r0p6657752_B0ToKpPimMuMu_11144002_magdown.root"); @@ -82,6 +86,7 @@ int new_analysis_b02kppimmumu() h1_B_Mass_jpsi->GetXaxis()->SetTitle(end_state_mass_literal); h1_B_Mass_psi2s->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_sim->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); @@ -106,6 +111,11 @@ int new_analysis_b02kppimmumu() std::cout << "# Starting evaluation of data." << std::endl; + int kplus = 0; + int kminus = 0; + int piplus = 0; + int piminus = 0; + for (unsigned int i = 0; i < data_entries; i++) { data_chain->GetEntry(i); @@ -120,15 +130,36 @@ int new_analysis_b02kppimmumu() FillHlt1DecisionHistos(hlt1_decision_histos, reconstructed_B_Mass); } - if (TMath::Abs(reconstructed_Kstar.M() - KSTAR_MASS) < 50.) + if (B_MASS_VAR_MIN <= reconstructed_B_Mass && reconstructed_B_Mass <= B_MASS_VAR_MAX) { - if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) + if (TMath::Abs(reconstructed_Kstar.M() - KSTAR_MASS) < 50.) { - 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); + if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) + { + if (Hp_Q == 1) + { + kplus++; + } + else if (Hp_Q == -1) + { + kminus++; + } + + if (Hm_Q == 1) + { + piplus++; + } + else if (Hm_Q == -1) + { + piminus++; + } + + 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); + } } } @@ -141,10 +172,10 @@ int new_analysis_b02kppimmumu() 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); + // 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_sim = CreateRooFitHistogramAndFitCB(h1_B_Mass_sim, false, false, ShapeParamters{}); auto roofit_hist_jpsi_fitsum = CreateRooFitHistogramAndFitCB(h1_B_Mass_jpsi, true, true, roofit_hist_sim.shape_parameters); @@ -154,7 +185,9 @@ int new_analysis_b02kppimmumu() 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); + + std::cout << "hist entries: " << h1_B_Mass_jpsi->GetEntries() << ", kplus: " << kplus << ", kminus: " << kminus << ", piplus: " << piplus << ", piminus: " << piminus << std::endl; return 0; } \ No newline at end of file diff --git a/new_analysis_bu2hpmumu.cpp b/new_analysis_bu2hpmumu.cpp index 7fbf69b..a4226e4 100644 --- a/new_analysis_bu2hpmumu.cpp +++ b/new_analysis_bu2hpmumu.cpp @@ -48,7 +48,7 @@ int new_analysis_bu2hpmumu() const char *data_tree_name = "BuToHpMuMu"; const char *sim_tree_name = "BuToHpMuMu_noPID"; const char *end_state_mass_literal = "m(#pi^{+}_{(#rightarrow K^{+})}#mu^{+}#mu^{-})"; - const bool retrain_bdt = false; + const bool retrain_bdt = true; TChain *data_chain = new TChain(TString::Format("%s/DecayTree", data_tree_name)); // data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0_tuple_90000000_v0r0p6288631.root"); @@ -81,6 +81,7 @@ int new_analysis_bu2hpmumu() h1_B_Mass_jpsi->GetXaxis()->SetTitle(end_state_mass_literal); h1_B_Mass_psi2s->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_sim->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); @@ -108,8 +109,8 @@ int new_analysis_bu2hpmumu() ConnectVarsToData(vars, data_chain, sim_chain, sig_tree, bkg_tree); - unsigned int data_entries = 100000;// data_chain->GetEntries(); - unsigned int sim_entries = 100000; //sim_chain->GetEntries(); + unsigned int data_entries = data_chain->GetEntries(); + unsigned int sim_entries = sim_chain->GetEntries(); std::cout << "# Got " << data_entries << " data and " << sim_entries << " simulated events." << std::endl; @@ -150,22 +151,19 @@ 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 (retrain_bdt && std::all_of(vars.begin(), vars.end(), [](TV *v) - { return v->IsMCFinite(); })) + if (sig_events < bkg_events) { - - sig_tree->Fill(); + if (retrain_bdt && std::all_of(vars.begin(), vars.end(), [](TV *v) + { return v->IsMCFinite(); })) + { + sig_tree->Fill(); + sig_events++; + } } - 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); } diff --git a/new_analysis_bu2kpmumu.cpp b/new_analysis_bu2kpmumu.cpp index 0adf66a..e68443c 100644 --- a/new_analysis_bu2kpmumu.cpp +++ b/new_analysis_bu2kpmumu.cpp @@ -79,6 +79,7 @@ int new_analysis_bu2kpmumu() h1_B_Mass_jpsi->GetXaxis()->SetTitle(end_state_mass_literal); h1_B_Mass_psi2s->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_sim->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); @@ -134,10 +135,10 @@ int new_analysis_bu2kpmumu() 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); + // 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_sim = CreateRooFitHistogramAndFitCB(h1_B_Mass_sim, false, false, ShapeParamters{}); auto roofit_hist_jpsi_fitsum = CreateRooFitHistogramAndFitCB(h1_B_Mass_jpsi, true, true, roofit_hist_sim.shape_parameters); @@ -147,7 +148,7 @@ int new_analysis_bu2kpmumu() 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); return 0; } \ No newline at end of file