From 16c136b1098b4835eadb17bbd9d43f45942323c6 Mon Sep 17 00:00:00 2001 From: Marius Pfeiffer Date: Wed, 6 Mar 2024 16:17:59 +0100 Subject: [PATCH] trees, unbinned fit, stacked hists, b2cc line --- basic_analysis.h | 112 ++++++++++++++++----- new_analysis_b02hphmmumu.cpp | 63 ++++++++---- new_analysis_b02kppimmumu.cpp | 91 ++++++++++------- new_analysis_bu2hpmumu.cpp | 157 ++++++++++------------------- new_analysis_bu2jpsikplus.cpp | 179 ++++++++++++++++++++++++++++++++++ new_analysis_bu2kpmumu.cpp | 50 ++++++---- 6 files changed, 451 insertions(+), 201 deletions(-) create mode 100644 new_analysis_bu2jpsikplus.cpp diff --git a/basic_analysis.h b/basic_analysis.h index df8d5a4..6aa0dae 100644 --- a/basic_analysis.h +++ b/basic_analysis.h @@ -98,6 +98,16 @@ struct FittedParam this->decimals = decimals; } + FittedParam(std::string name, std::string title, double value, double err, int decimals) { + this->title = title; + this->name = name; + this->value = value; + this->error = err; + this->decimals = decimals; + this->constant = false; + this->at_limit = false; + } + std::string ToString() const { if (this->constant) @@ -128,7 +138,7 @@ struct RooFitSummary ShapeParamters shape_parameters; }; -void DrawInDefaultCanvas(TH1 *histogram, const char *folder, double margin_left = 0, Option_t *option = "") +void DrawInDefaultCanvas(TH1 *histogram, const char *folder, double margin_left = 0.1, Option_t *option = "") { std::filesystem::create_directory(TString::Format("output_files/analysis/%s", folder).Data()); TString name = TString::Format("%s_canvas", histogram->GetName()); @@ -140,6 +150,43 @@ void DrawInDefaultCanvas(TH1 *histogram, const char *folder, double margin_left c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data()); } +void DrawInDefaultCanvasStacked(std::vector histograms, std::vector colors, std::vector fill_style, const char *folder, double margin_left = 0.1, Option_t *option = "") +{ + std::filesystem::create_directory(TString::Format("output_files/analysis/%s", folder).Data()); + TString name = TString::Format("%s_stack_canvas", histograms[0]->GetName()); + TCanvas *c = new TCanvas(name, histograms[0]->GetName(), 0, 0, 800, 600); + c->SetLeftMargin(margin_left); + std::string drwopt_2 = std::string(option).empty() ? "SAME HIST" : TString::Format("%s SAME HIST", option).Data(); + for (size_t i = 0; i < histograms.size(); i++) + { + const Double_t scaling_factor = 1.; + auto hist_clone = (TH1*)histograms[i]->Clone(TString::Format("%s_clone", histograms[i]->GetName())); + hist_clone->Scale(scaling_factor / histograms[i]->GetMaximum()); + hist_clone->SetLineColor(colors[i]); + hist_clone->SetMaximum(scaling_factor + (scaling_factor * 0.05)); + hist_clone->SetMinimum(0.); + hist_clone->SetStats(0); + if (fill_style[i] != 0) { + hist_clone->SetFillStyle(fill_style[i]); + hist_clone->SetFillColor(colors[i]); + } + + if (i > 0) + { + hist_clone->Draw(drwopt_2.c_str()); + } + else + { + hist_clone->Draw(drwopt_2.c_str()); + } + } + + c->BuildLegend(0.55, 0.79, 0.87, 0.89); + + c->Draw(); + c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data()); +} + void DrawInDefaultCanvas(RooPlot *rooPlot, const char *folder, double margin_left = 0, Option_t *option = "") { std::filesystem::create_directory(TString::Format("output_files/analysis/%s", folder).Data()); @@ -151,7 +198,8 @@ 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) { +void PrintStyles(RooPlot *plt) +{ auto Xaxis = plt->GetXaxis(); auto Yaxis = plt->GetYaxis(); std::cout << "## Styles for plot " << plt->GetName() << std::endl; @@ -227,9 +275,9 @@ void DrawInDefaultCanvas(RooFitSummary fitSummary, const char *folder) fitSummary.fit_histogram->Draw(); - TLegend *leg1 = new TLegend(0.60, 0.55, 0.87, 0.89); + TLegend *leg1 = new TLegend(0.55, 0.45, 0.87, 0.89); leg1->SetFillColor(kWhite); - leg1->SetLineColor(kBlack); + leg1->SetLineColor(kWhite); for (const auto &[name, title] : fitSummary.pdf_names) { @@ -283,17 +331,17 @@ void DrawInDefaultCanvas(RooFitSummary fitSummary, const char *folder) 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->SetLineColor(kOrange + 4); line_mid->Draw(); 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->SetLineColor(kOrange + 4); line_up->Draw(); 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->SetLineColor(kOrange + 4); line_low->Draw(); c->Draw(); @@ -324,23 +372,27 @@ RooPlot *CreateRooFitHistogram(TH1D *hist) return roo_frame_mass; } -RooFitSummary CreateRooFitHistogramAndFitCB(TH1D *hist, bool hasExpBkg, bool useExtShape, ShapeParamters extShape) +RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString xAxis, bool hasExpBkg, bool useExtShape, ShapeParamters extShape, Double_t fit_low = 0, Double_t fit_up = 0) { - auto suffix_name = [name = hist->GetName()](const char *text) + auto suffix_name = [name = dataSet->GetName()](const char *text) { return TString::Format("%s_%s", text, name); }; - RooRealVar roo_var_mass(suffix_name("var_mass"), "B Mass Variable", B_MASS_VAR_MIN, B_MASS_VAR_MAX); + Double_t fitRangeUp = fit_low == 0 ? B_MASS_VAR_MIN : fit_low; + Double_t fitRangeLow = fit_up == 0 ? B_MASS_VAR_MAX : fit_up; + + RooRealVar roo_var_mass(var_name, "B Mass Variable", B_MASS_VAR_MIN, B_MASS_VAR_MAX); const char *fitting_range_name = "fitting_range"; - roo_var_mass.setRange(fitting_range_name, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + roo_var_mass.setRange(fitting_range_name, fitRangeUp, fitRangeLow); - TString hist_name = suffix_name("roohist_B_M"); - RooDataHist roohist_B_M(hist_name, "B Mass Histogram", roo_var_mass, RooFit::Import(*hist)); + TString dataset_name = suffix_name("roodataset_B_M"); + //RooDataHist roodataset_B_M(hist_name, "B Mass Histogram", roo_var_mass, RooFit::Import(*hist)); + RooDataSet roodataset_B_M(dataset_name, "B Mass Data Set", roo_var_mass, RooFit::Import(*dataSet)); - RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title(hist->GetTitle()), RooFit::Name(TString::Format("%s_rplt", hist->GetName()))); - roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(B_MASS_HIST_BINS), RooFit::Name(hist_name)); - roo_frame_mass->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle()); + RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title(dataSet->GetTitle()), RooFit::Name(TString::Format("%s_rplt", dataSet->GetName()))); + roodataset_B_M.plotOn(roo_frame_mass, RooFit::Name(dataset_name)); + roo_frame_mass->GetXaxis()->SetTitle(xAxis); // Crystal Ball for Signal RooRealVar var_mass_x0(suffix_name("var_mass_x0"), "#mu", 5278., 5170., 5500.); @@ -384,29 +436,43 @@ RooFitSummary CreateRooFitHistogramAndFitCB(TH1D *hist, bool hasExpBkg, bool use RooExponential bkg_exp(background_name, "Exp Background", roo_var_mass, var_mass_bkg_c); pdf_names[background_name.Data()] = bkg_exp.getTitle().Data(); - RooRealVar var_mass_nsig(suffix_name("nsig"), "N_{Sig}", 0., hist->GetEntries()); - RooRealVar var_mass_nbkg(suffix_name("nbkg"), "N_{Bkg}", 0., hist->GetEntries()); + RooRealVar var_mass_nsig(suffix_name("nsig"), "N_{Sig}", 0., dataSet->GetEntries()); + RooRealVar var_mass_nbkg(suffix_name("nbkg"), "N_{Bkg}", 0., dataSet->GetEntries()); TString fitted_pdf_name = suffix_name("sigplusbkg"); RooAddPdf fitted_pdf(fitted_pdf_name, "Sig + Bkg", RooArgList(sig_cb, bkg_exp), RooArgList(var_mass_nsig, var_mass_nbkg)); pdf_names[fitted_pdf_name.Data()] = fitted_pdf.getTitle().Data(); - RooFitResult *fitres = fitted_pdf.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range(fitting_range_name)); + RooFitResult *fitres = fitted_pdf.fitTo(roodataset_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range(fitting_range_name)); // sigplusbkg.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144)); fitted_pdf.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range(fitting_range_name), RooFit::Name(fitted_pdf_name)); fitted_pdf.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(bkg_exp)), RooFit::LineColor(kBlue - 7), RooFit::LineStyle(kDashed), RooFit::Range(fitting_range_name), RooFit::Name(background_name)); fitted_pdf.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(sig_cb)), RooFit::LineColor(kRed - 7), RooFit::LineStyle(kDashed), RooFit::Range(fitting_range_name), RooFit::Name(signal_name)); - fitted_params.push_back(FittedParam(var_mass_bkg_c, 5)); fitted_params.push_back(FittedParam(var_mass_nsig, 2)); fitted_params.push_back(FittedParam(var_mass_nbkg, 2)); + + + double sig_val = var_mass_nsig.getVal(); + double sig_err = var_mass_nsig.getError(); + double bkg_val = var_mass_nbkg.getVal(); + double bkg_err = var_mass_nbkg.getError(); + + double sig_over_bkg_val = sig_val / TMath::Sqrt(sig_val + bkg_val); + double err_prop_sig = (sig_val + 2 * bkg_val)/(2 * TMath::Power((sig_val + bkg_val), (3/2))); + double err_prop_bkg = -sig_val/(2 * TMath::Power((sig_val + bkg_val), (3/2))); + double sig_over_bkg_err = TMath::Sqrt(TMath::Sq(err_prop_sig * sig_err) + TMath::Sq(err_prop_bkg * bkg_err)); + + fitted_params.push_back(FittedParam("sig_over_bkg", "N_{Sig}/#sqrt{N_{Sig} + N_{Bkg}}", sig_over_bkg_val, sig_over_bkg_err, 2)); + + fitted_params.push_back(FittedParam(var_mass_bkg_c, 5)); pull_compare_name = fitted_pdf_name; } else { - RooFitResult *fitres = sig_cb.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range(fitting_range_name)); + RooFitResult *fitres = sig_cb.fitTo(roodataset_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range(fitting_range_name)); // sigplusbkg.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_name), RooFit::Name(signal_name)); pull_compare_name = signal_name; @@ -419,8 +485,8 @@ RooFitSummary CreateRooFitHistogramAndFitCB(TH1D *hist, bool hasExpBkg, bool use fitted_params.push_back(FittedParam(var_mass_alphaR, 2)); 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()))); - RooHist *pull_hist = roo_frame_mass->pullHist(hist_name, pull_compare_name); + RooPlot *roo_frame_pull = roo_var_mass.frame(RooFit::Title("Pull Distribution"), RooFit::Name(TString::Format("%s_pull_rplt", dataSet->GetName()))); + RooHist *pull_hist = roo_frame_mass->pullHist(dataset_name, pull_compare_name); roo_frame_pull->addPlotable(pull_hist, "P"); // Double_t pull_min = pull_hist->GetMinimum(); diff --git a/new_analysis_b02hphmmumu.cpp b/new_analysis_b02hphmmumu.cpp index 74811b8..d6a0d95 100644 --- a/new_analysis_b02hphmmumu.cpp +++ b/new_analysis_b02hphmmumu.cpp @@ -79,16 +79,30 @@ int new_analysis_b02hphmmumu() sim_chain->SetBranchAddress("Hm_TRUEID", &Hm_TRUEID); sim_chain->SetBranchAddress("B0_BKGCAT", &B_BKGCAT); - 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); - TH1D *h1_B_Mass_sim = new TH1D("h1_B_Mass_sim", TString::Format("B Mass, Simualted (%s_noPID)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + Double_t B_Mass_jpsi_var, B_Mass_psi2s_var, B_Mass_sim_var; + TString B_Mass_jpsi_var_name = "B_Mass_jpsi_var"; + TString B_Mass_psi2s_var_name = "B_Mass_psi2s_var"; + TString B_Mass_sim_var_name = "B_Mass_sim_var"; + + TTree *tree_B_Mass_jpsi = new TTree("tree_B_Mass_jpsi", TString::Format("B Mass, J/#psi Mode (%s)", analysis_name)); + TTree *tree_B_Mass_psi2s = new TTree("tree_B_Mass_psi2s", TString::Format("B Mass, #psi(2S) Mode (%s)", analysis_name)); + TTree *tree_B_Mass_sim = new TTree("tree_B_Mass_sim", TString::Format("B Mass, Simualted (%s_noPID)", analysis_name)); + + tree_B_Mass_jpsi->Branch(B_Mass_jpsi_var_name, &B_Mass_jpsi_var); + tree_B_Mass_psi2s->Branch(B_Mass_psi2s_var_name, &B_Mass_psi2s_var); + tree_B_Mass_sim->Branch(B_Mass_sim_var_name, &B_Mass_sim_var); + + TH1D *h1_B_Mass_unf = new TH1D("h1_B_Mass_unf", TString::Format("B Mass (%s), Unfiltered", data_tree_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TH1D *h1_B_Mass_bdtf = new TH1D("h1_B_Mass_bdtf", TString::Format("B Mass (%s), BDT Filter", data_tree_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TH1D *h1_B_Mass_sim_unf = new TH1D("h1_B_Mass_sim_unf", TString::Format("B Mass, Simualted (%s), Unfiltered", sim_tree_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.); TH1D *h1_bdt_probs = new TH1D("h1_bdt_probs", "BDT Probabilities", 100, -1, 1); - 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); + h1_B_Mass_unf->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_sim_unf->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_bdtf->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); @@ -173,6 +187,11 @@ int new_analysis_b02hphmmumu() { sim_chain->GetEntry(i); + Double_t reco_mass_pipkp = (hp4v_sim->LorentzVector(K_MASS) + hm4v_sim->LorentzVector() + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M(); + Double_t reco_mass_pimkm = (hp4v_sim->LorentzVector() + hm4v_sim->LorentzVector(K_MASS) + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M(); + h1_B_Mass_sim_unf->Fill(reco_mass_pipkp); + h1_B_Mass_sim_unf->Fill(reco_mass_pimkm); + if (B_BKGCAT == 30 && TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID) { TLorentzVector reconstructed_Kstar{}; @@ -202,7 +221,8 @@ int new_analysis_b02hphmmumu() } } - h1_B_Mass_sim->Fill(reconstructed_B_Mass); + B_Mass_sim_var = reconstructed_B_Mass; + tree_B_Mass_sim->Fill(); } } @@ -269,17 +289,22 @@ int new_analysis_b02hphmmumu() double mva_response = reader->EvaluateMVA("BDT"); h1_bdt_probs->Fill(mva_response); + h1_B_Mass_unf->Fill(reconstructed_B_Mass); + if (mva_response > mva_cut_value) { - if (TMath::Abs(reconstructed_Kstar.M() - KSTAR_MASS) < 50.) + h1_B_Mass_bdtf->Fill(reconstructed_B_Mass); + 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); + B_Mass_jpsi_var = reconstructed_B_Mass; + tree_B_Mass_jpsi->Fill(); } else if (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.) { - h1_B_Mass_psi2s->Fill(reconstructed_B_Mass); + B_Mass_psi2s_var = reconstructed_B_Mass; + tree_B_Mass_psi2s->Fill(); } } } @@ -294,14 +319,18 @@ int new_analysis_b02hphmmumu() 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_unf, analysis_name, 0.1); + DrawInDefaultCanvas(h1_B_Mass_sim_unf, analysis_name, 0.1); + DrawInDefaultCanvas(h1_B_Mass_bdtf, analysis_name, 0.1); + + DrawInDefaultCanvasStacked({h1_B_Mass_unf, h1_B_Mass_bdtf}, {kRed, kBlue}, {0, 3003}, analysis_name); - 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); - auto roofit_hist_psi2s_fitsum = CreateRooFitHistogramAndFitCB(h1_B_Mass_psi2s, true, true, roofit_hist_sim.shape_parameters); + auto roofit_hist_sim = CreateRooDataSetAndFitCB(tree_B_Mass_sim, B_Mass_sim_var_name, end_state_mass_literal, false, false, ShapeParamters{}); + auto roofit_hist_jpsi_fitsum = CreateRooDataSetAndFitCB(tree_B_Mass_jpsi, B_Mass_jpsi_var_name, end_state_mass_literal, true, true, roofit_hist_sim.shape_parameters); + auto roofit_hist_psi2s_fitsum = CreateRooDataSetAndFitCB(tree_B_Mass_psi2s, B_Mass_psi2s_var_name, end_state_mass_literal, true, true, roofit_hist_sim.shape_parameters); DrawInDefaultCanvas(roofit_hist_jpsi_fitsum, analysis_name); DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name); diff --git a/new_analysis_b02kppimmumu.cpp b/new_analysis_b02kppimmumu.cpp index f2b0c0c..b1d607f 100644 --- a/new_analysis_b02kppimmumu.cpp +++ b/new_analysis_b02kppimmumu.cpp @@ -78,15 +78,27 @@ int new_analysis_b02kppimmumu() sim_chain->SetBranchAddress("Pi_TRUEID", &Hm_TRUEID); sim_chain->SetBranchAddress("B0_BKGCAT", &B_BKGCAT); - 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); - TH1D *h1_B_Mass_sim = new TH1D("h1_B_Mass_sim", TString::Format("B Mass, Simualted (%s_noPID)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + Double_t B_Mass_jpsi_var, B_Mass_psi2s_var, B_Mass_sim_var; + TString B_Mass_jpsi_var_name = "B_Mass_jpsi_var"; + TString B_Mass_psi2s_var_name = "B_Mass_psi2s_var"; + TString B_Mass_sim_var_name = "B_Mass_sim_var"; + + TTree *tree_B_Mass_jpsi = new TTree("tree_B_Mass_jpsi", TString::Format("B Mass, J/#psi Mode (%s)", analysis_name)); + TTree *tree_B_Mass_psi2s = new TTree("tree_B_Mass_psi2s", TString::Format("B Mass, #psi(2S) Mode (%s)", analysis_name)); + TTree *tree_B_Mass_sim = new TTree("tree_B_Mass_sim", TString::Format("B Mass, Simualted (%s_noPID)", analysis_name)); + + tree_B_Mass_jpsi->Branch(B_Mass_jpsi_var_name, &B_Mass_jpsi_var); + tree_B_Mass_psi2s->Branch(B_Mass_psi2s_var_name, &B_Mass_psi2s_var); + tree_B_Mass_sim->Branch(B_Mass_sim_var_name, &B_Mass_sim_var); + + TH1D *h1_B_Mass_unf = new TH1D("h1_B_Mass_unf", TString::Format("B Mass (%s), Unfiltered", data_tree_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TH1D *h1_B_Mass_sim_unf = new TH1D("h1_B_Mass_sim_unf", TString::Format("B Mass, Simualted (%s), Unfiltered", sim_tree_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); - h1_B_Mass_sim->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_unf->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_sim_unf->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); @@ -101,9 +113,11 @@ int new_analysis_b02kppimmumu() { sim_chain->GetEntry(i); Double_t reconstructed_B_Mass = (hp4v_sim->LorentzVector() + hm4v_sim->LorentzVector() + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M(); + h1_B_Mass_sim_unf->Fill(reconstructed_B_Mass); if (B_BKGCAT == 0 && TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID && TMath::Abs(Hp_TRUEID) == PID_KAON && TMath::Abs(Hm_TRUEID) == PID_PION) { - h1_B_Mass_sim->Fill(reconstructed_B_Mass); + B_Mass_sim_var = reconstructed_B_Mass; + tree_B_Mass_sim->Fill(); } PrintProgress(TString::Format("%s SIG Collection", analysis_name), sim_entries, 10000, i); @@ -130,36 +144,36 @@ int new_analysis_b02kppimmumu() FillHlt1DecisionHistos(hlt1_decision_histos, reconstructed_B_Mass); } - if (B_MASS_VAR_MIN <= reconstructed_B_Mass && reconstructed_B_Mass <= B_MASS_VAR_MAX) + h1_B_Mass_unf->Fill(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.) { - if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) + if (Hp_Q == 1) { - 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); + kplus++; } - else if (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.) + else if (Hp_Q == -1) + { + kminus++; + } + + if (Hm_Q == 1) { - h1_B_Mass_psi2s->Fill(reconstructed_B_Mass); + piplus++; } + else if (Hm_Q == -1) + { + piminus++; + } + + B_Mass_jpsi_var = reconstructed_B_Mass; + tree_B_Mass_jpsi->Fill(); + } + else if (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.) + { + B_Mass_psi2s_var = reconstructed_B_Mass; + tree_B_Mass_psi2s->Fill(); } } @@ -172,14 +186,17 @@ 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(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); - auto roofit_hist_psi2s_fitsum = CreateRooFitHistogramAndFitCB(h1_B_Mass_psi2s, true, true, roofit_hist_sim.shape_parameters); + DrawInDefaultCanvas(h1_B_Mass_unf, analysis_name, 0.1); + DrawInDefaultCanvas(h1_B_Mass_sim_unf, analysis_name, 0.1); + + auto roofit_hist_sim = CreateRooDataSetAndFitCB(tree_B_Mass_sim, B_Mass_sim_var_name, end_state_mass_literal, false, false, ShapeParamters{}); + auto roofit_hist_jpsi_fitsum = CreateRooDataSetAndFitCB(tree_B_Mass_jpsi, B_Mass_jpsi_var_name, end_state_mass_literal, true, true, roofit_hist_sim.shape_parameters); + auto roofit_hist_psi2s_fitsum = CreateRooDataSetAndFitCB(tree_B_Mass_psi2s, B_Mass_psi2s_var_name, end_state_mass_literal, true, true, roofit_hist_sim.shape_parameters); DrawInDefaultCanvas(roofit_hist_jpsi_fitsum, analysis_name); DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name); @@ -187,7 +204,7 @@ int new_analysis_b02kppimmumu() // 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; + std::cout << "hist entries: " << tree_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 a4226e4..378d4af 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 = true; + 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"); @@ -66,22 +66,40 @@ int new_analysis_bu2hpmumu() FourVect *hp4v_sim = FourVect::Init(sim_chain, "Hp"); Int_t B_BKGCAT, L1_TRUEID, L2_TRUEID, Hp_TRUEID; + Double_t Hp_PID_K_sim; sim_chain->SetBranchAddress("L1_TRUEID", &L1_TRUEID); sim_chain->SetBranchAddress("L2_TRUEID", &L2_TRUEID); sim_chain->SetBranchAddress("Hp_TRUEID", &Hp_TRUEID); sim_chain->SetBranchAddress("B_BKGCAT", &B_BKGCAT); + sim_chain->SetBranchAddress("Hp_PID_K", &Hp_PID_K_sim); + + Double_t B_Mass_jpsi_var, B_Mass_psi2s_var, B_Mass_sim_var; + TString B_Mass_jpsi_var_name = "B_Mass_jpsi_var"; + TString B_Mass_psi2s_var_name = "B_Mass_psi2s_var"; + TString B_Mass_sim_var_name = "B_Mass_sim_var"; + + TTree *tree_B_Mass_jpsi = new TTree("tree_B_Mass_jpsi", TString::Format("B Mass, J/#psi Mode (%s)", analysis_name)); + TTree *tree_B_Mass_psi2s = new TTree("tree_B_Mass_psi2s", TString::Format("B Mass, #psi(2S) Mode (%s)", analysis_name)); + TTree *tree_B_Mass_sim = new TTree("tree_B_Mass_sim", TString::Format("B Mass, Simualted (%s_noPID)", analysis_name)); + + tree_B_Mass_jpsi->Branch(B_Mass_jpsi_var_name, &B_Mass_jpsi_var); + tree_B_Mass_psi2s->Branch(B_Mass_psi2s_var_name, &B_Mass_psi2s_var); + tree_B_Mass_sim->Branch(B_Mass_sim_var_name, &B_Mass_sim_var); + + TH1D *h1_B_Mass_unf = new TH1D("h1_B_Mass_unf", TString::Format("B Mass (%s), Unfiltered", data_tree_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TH1D *h1_B_Mass_bdtf = new TH1D("h1_B_Mass_bdtf", TString::Format("B Mass (%s), BDT Filter", data_tree_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TH1D *h1_B_Mass_sim_unf = new TH1D("h1_B_Mass_sim_unf", TString::Format("B Mass, Simualted (%s), Unfiltered", sim_tree_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); - 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); - TH1D *h1_B_Mass_sim = new TH1D("h1_B_Mass_sim", TString::Format("B Mass, Simualted (%s_noPID)", 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.); TH1D *h1_bdt_probs = new TH1D("h1_bdt_probs", "BDT Probabilities", 100, -2, 2); + TH1D *h1_Hp_PID_K_pref = new TH1D("h1_Hp_PID_K_pref", "H^{+} PID K, Before TrueID", 50, -10, 10); + TH1D *h1_Hp_PID_K_postf = new TH1D("h1_Hp_PID_K_postf", "H^{+} PID K, After TrueID", 50, -10, 10); - 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); + h1_B_Mass_unf->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_sim_unf->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_bdtf->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); @@ -149,8 +167,11 @@ int new_analysis_bu2hpmumu() { sim_chain->GetEntry(i); Double_t reconstructed_B_Mass = (hp4v_sim->LorentzVector(K_MASS) + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M(); + h1_Hp_PID_K_pref->Fill(Hp_PID_K_sim); + h1_B_Mass_sim_unf->Fill(reconstructed_B_Mass); if (B_BKGCAT == 30 && TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID && TMath::Abs(Hp_TRUEID) == PID_KAON) { + h1_Hp_PID_K_postf->Fill(Hp_PID_K_sim); if (sig_events < bkg_events) { if (retrain_bdt && std::all_of(vars.begin(), vars.end(), [](TV *v) @@ -161,7 +182,8 @@ int new_analysis_bu2hpmumu() } } - h1_B_Mass_sim->Fill(reconstructed_B_Mass); + B_Mass_sim_var = reconstructed_B_Mass; + tree_B_Mass_sim->Fill(); } PrintProgress(TString::Format("%s SIG Collection", analysis_name), sim_entries, 10000, i); @@ -178,7 +200,7 @@ int new_analysis_bu2hpmumu() Float_t *train_vars = new Float_t[vars.size()]; auto reader = SetupReader(vars, train_vars, "BuToHpMuMu"); - const double mva_cut_value = -0.03; + const double mva_cut_value = -0.02; for (unsigned int i = 0; i < data_entries; i++) { @@ -212,15 +234,20 @@ int new_analysis_bu2hpmumu() double mva_response = reader->EvaluateMVA("BDT"); h1_bdt_probs->Fill(mva_response); + h1_B_Mass_unf->Fill(reconstructed_B_Mass); + if (mva_response > mva_cut_value) { + h1_B_Mass_bdtf->Fill(reconstructed_B_Mass); if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) { - h1_B_Mass_jpsi->Fill(reconstructed_B_Mass); + B_Mass_jpsi_var = reconstructed_B_Mass; + tree_B_Mass_jpsi->Fill(); } else if (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.) { - h1_B_Mass_psi2s->Fill(reconstructed_B_Mass); + B_Mass_psi2s_var = reconstructed_B_Mass; + tree_B_Mass_psi2s->Fill(); } } @@ -233,14 +260,19 @@ 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(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_unf, analysis_name, 0.1); + DrawInDefaultCanvas(h1_B_Mass_sim_unf, analysis_name, 0.1); + DrawInDefaultCanvas(h1_B_Mass_bdtf, analysis_name, 0.1); + + DrawInDefaultCanvasStacked({h1_Hp_PID_K_pref, h1_Hp_PID_K_postf}, {kRed, kBlue}, {0, 3003},, analysis_name); + DrawInDefaultCanvasStacked({h1_B_Mass_unf, h1_B_Mass_bdtf}, {kRed, kBlue}, {0, 3003},, analysis_name); - 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); - auto roofit_hist_psi2s_fitsum = CreateRooFitHistogramAndFitCB(h1_B_Mass_psi2s, true, true, roofit_hist_sim.shape_parameters); + auto roofit_hist_sim = CreateRooDataSetAndFitCB(tree_B_Mass_sim, B_Mass_sim_var_name, end_state_mass_literal, false, false, ShapeParamters{}); + auto roofit_hist_jpsi_fitsum = CreateRooDataSetAndFitCB(tree_B_Mass_jpsi, B_Mass_jpsi_var_name, end_state_mass_literal, true, true, roofit_hist_sim.shape_parameters); + auto roofit_hist_psi2s_fitsum = CreateRooDataSetAndFitCB(tree_B_Mass_psi2s, B_Mass_psi2s_var_name, end_state_mass_literal, true, true, roofit_hist_sim.shape_parameters); DrawInDefaultCanvas(roofit_hist_jpsi_fitsum, analysis_name); DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name); @@ -252,92 +284,3 @@ int new_analysis_bu2hpmumu() return 0; } - -int analyze_bu2hpmumu_simulation() -{ - const char *analysis_name = "BuToHpMuMu_noPID"; - const char *end_state_mass_literal = "m(#pi^{+}_{(#rightarrow K^{+})}#mu^{+}#mu^{-})"; - - TChain *sim_chain = new TChain(TString::Format("%s/DecayTree", analysis_name)); - sim_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/rd_btoxll_simulation_fullstream_v0r0p6671378_BuToKpMuMu_12143001_magdown.root"); - - FourVect *l14v = FourVect::Init(sim_chain, "L1"); - FourVect *l24v = FourVect::Init(sim_chain, "L2"); - FourVect *hp4v = FourVect::Init(sim_chain, "Hp"); - - Int_t B_BKGCAT, L1_TRUEID, L2_TRUEID, Hp_TRUEID; - - sim_chain->SetBranchAddress("L1_TRUEID", &L1_TRUEID); - sim_chain->SetBranchAddress("L2_TRUEID", &L2_TRUEID); - sim_chain->SetBranchAddress("Hp_TRUEID", &Hp_TRUEID); - sim_chain->SetBranchAddress("B_BKGCAT", &B_BKGCAT); - - 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(sim_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 = sim_chain->GetEntries(); - for (unsigned int i = 0; i < entries; i++) - { - sim_chain->GetEntry(i); - - if (B_BKGCAT == 30 && TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID && TMath::Abs(Hp_TRUEID) == PID_KAON) - { - - TLorentzVector dimuon = l14v->LorentzVector() + l24v->LorentzVector(); - Double_t reconstructed_B_Mass = (hp4v->LorentzVector(K_MASS) + 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(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; -} diff --git a/new_analysis_bu2jpsikplus.cpp b/new_analysis_bu2jpsikplus.cpp new file mode 100644 index 0000000..18911a8 --- /dev/null +++ b/new_analysis_bu2jpsikplus.cpp @@ -0,0 +1,179 @@ +#include +#include +#include +#include +#include +#include + +#include "TH1D.h" +#include "TH2D.h" +#include "THStack.h" +#include "TGraph.h" +#include "TTree.h" +#include "TChain.h" +#include "TFile.h" +#include "TCanvas.h" +#include "TROOT.h" +#include "TStyle.h" +#include "TColor.h" +#include "TLorentzVector.h" +#include "TRandom3.h" +#include "TLegend.h" + +#include "RooDataHist.h" +#include "RooRealVar.h" +#include "RooPlot.h" +#include "RooGaussian.h" +#include "RooExponential.h" +#include "RooRealConstant.h" +#include "RooAddPdf.h" +#include "RooFitResult.h" +#include "RooProduct.h" +#include "RooCrystalBall.h" +#include "RooBreitWigner.h" +#include "RooArgSet.h" +#include "RooFFTConvPdf.h" +#include "RooNovosibirsk.h" + +#include "constants.h" +#include "basic_analysis.h" +#include "hlt1_decision_analysis.h" +#include "bdt_classification.h" + +// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Exclusive B2CC Bu To Kp Mu Mu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +int new_analysis_bu2jpsikplus() +{ + const char *analysis_name = "BuToJpsiKplus"; + const char *data_tree_name = "Hlt2B2CC_BuToJpsiKplus_JpsiToMuMu_Detached"; + const char *sim_tree_name = "BuToKpMuMu_noPID"; + const char *end_state_mass_literal = "m(K^{+}#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/b2cc_turbo_superliaison_collision23_butojpsikplus_turbopass_magdown_v0r0p6205916.root"); + + FourVect *l14v_data = FourVect::Init(data_chain, "muplus"); + FourVect *l24v_data = FourVect::Init(data_chain, "muminus"); + FourVect *hp4v_data = FourVect::Init(data_chain, "Kplus"); + Bool_t muplus_ISMUON, muminus_ISMUON, Hlt1TrackMVADecision, Hlt1TwoTrackMVADecision; + + data_chain->SetBranchAddress("muplus_ISMUON", &muplus_ISMUON); + data_chain->SetBranchAddress("muminus_ISMUON", &muminus_ISMUON); + data_chain->SetBranchAddress("Hlt1TrackMVADecision", &Hlt1TrackMVADecision); + data_chain->SetBranchAddress("Hlt1TwoTrackMVADecision", &Hlt1TwoTrackMVADecision); + + 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_BuToKpMuMu_12143001_magdown.root"); + + FourVect *l14v_sim = FourVect::Init(sim_chain, "L1"); + FourVect *l24v_sim = FourVect::Init(sim_chain, "L2"); + FourVect *hp4v_sim = FourVect::Init(sim_chain, "K"); + + Int_t B_BKGCAT, L1_TRUEID, L2_TRUEID, Hp_TRUEID; + + sim_chain->SetBranchAddress("L1_TRUEID", &L1_TRUEID); + sim_chain->SetBranchAddress("L2_TRUEID", &L2_TRUEID); + sim_chain->SetBranchAddress("K_TRUEID", &Hp_TRUEID); + sim_chain->SetBranchAddress("B_BKGCAT", &B_BKGCAT); + + Double_t B_Mass_jpsi_var, B_Mass_psi2s_var, B_Mass_sim_var; + TString B_Mass_jpsi_var_name = "B_Mass_jpsi_var"; + TString B_Mass_psi2s_var_name = "B_Mass_psi2s_var"; + TString B_Mass_sim_var_name = "B_Mass_sim_var"; + + TTree *tree_B_Mass_jpsi = new TTree("tree_B_Mass_jpsi", TString::Format("B Mass, J/#psi Mode (%s)", analysis_name)); + TTree *tree_B_Mass_psi2s = new TTree("tree_B_Mass_psi2s", TString::Format("B Mass, #psi(2S) Mode (%s)", analysis_name)); + TTree *tree_B_Mass_sim = new TTree("tree_B_Mass_sim", TString::Format("B Mass, Simualted (%s_noPID)", analysis_name)); + + tree_B_Mass_jpsi->Branch(B_Mass_jpsi_var_name, &B_Mass_jpsi_var); + tree_B_Mass_psi2s->Branch(B_Mass_psi2s_var_name, &B_Mass_psi2s_var); + tree_B_Mass_sim->Branch(B_Mass_sim_var_name, &B_Mass_sim_var); + + TH1D *h1_B_Mass_unf = new TH1D("h1_B_Mass_unf", TString::Format("B Mass (%s), Unfiltered", data_tree_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TH1D *h1_B_Mass_sim_unf = new TH1D("h1_B_Mass_sim_unf", TString::Format("B Mass, Simualted (%s), Unfiltered", sim_tree_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_unf->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_sim_unf->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 data_entries = data_chain->GetEntries(); + unsigned int sim_entries = sim_chain->GetEntries(); + + for (unsigned int i = 0; i < sim_entries; i++) + { + sim_chain->GetEntry(i); + Double_t reconstructed_B_Mass = (hp4v_sim->LorentzVector() + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M(); + h1_B_Mass_sim_unf->Fill(reconstructed_B_Mass); + if (B_BKGCAT == 0 && TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID && TMath::Abs(Hp_TRUEID) == PID_KAON) + { + B_Mass_sim_var = reconstructed_B_Mass; + tree_B_Mass_sim->Fill(); + } + + PrintProgress(TString::Format("%s SIG Collection", analysis_name), sim_entries, 10000, i); + } + + std::cout << "# Starting evaluation of data." << std::endl; + + for (unsigned int i = 0; i < data_entries; i++) + { + data_chain->GetEntry(i); + + TLorentzVector dimuon = l14v_data->LorentzVector() + l24v_data->LorentzVector(); + Double_t reconstructed_B_Mass = (hp4v_data->LorentzVector() + 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); + // } + + h1_B_Mass_unf->Fill(reconstructed_B_Mass); + + if (muplus_ISMUON && muminus_ISMUON && (Hlt1TrackMVADecision || Hlt1TwoTrackMVADecision)) + { + if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) + { + B_Mass_jpsi_var = reconstructed_B_Mass; + tree_B_Mass_jpsi->Fill(); + } + else if (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.) + { + B_Mass_psi2s_var = reconstructed_B_Mass; + tree_B_Mass_psi2s->Fill(); + } + } + + PrintProgress(TString::Format("%s BDT Evaluation", analysis_name), data_entries, 10000, i); + } + + // std::cout << "# Exclusive Hits" << std::endl; + // for (const auto &[line, hits] : exclusive_hits) + // { + // std::cout << line << ": " << hits << std::endl; + // } + + DrawInDefaultCanvas(h1_B_Mass_unf, analysis_name, 0.1); + DrawInDefaultCanvas(h1_B_Mass_sim_unf, analysis_name, 0.1); + + auto roofit_hist_sim = CreateRooDataSetAndFitCB(tree_B_Mass_sim, B_Mass_sim_var_name, end_state_mass_literal, false, false, ShapeParamters{}); + auto roofit_hist_jpsi_fitsum = CreateRooDataSetAndFitCB(tree_B_Mass_jpsi, B_Mass_jpsi_var_name, end_state_mass_literal, true, true, roofit_hist_sim.shape_parameters, 5140., 5460.); + auto roofit_hist_psi2s_fitsum = CreateRooDataSetAndFitCB(tree_B_Mass_psi2s, B_Mass_psi2s_var_name, end_state_mass_literal, true, true, roofit_hist_sim.shape_parameters, 5140., 5460.); + + DrawInDefaultCanvas(roofit_hist_jpsi_fitsum, analysis_name); + DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name); + DrawInDefaultCanvas(roofit_hist_sim, analysis_name); + + // DrawHlt1DecisionHistos(analysis_name, hlt1_decision_histos); + + return 0; +} \ No newline at end of file diff --git a/new_analysis_bu2kpmumu.cpp b/new_analysis_bu2kpmumu.cpp index e68443c..ed87c63 100644 --- a/new_analysis_bu2kpmumu.cpp +++ b/new_analysis_bu2kpmumu.cpp @@ -48,7 +48,7 @@ int new_analysis_bu2kpmumu() const char *data_tree_name = "Hlt2RD_BuToKpMuMu"; const char *sim_tree_name = "BuToKpMuMu_noPID"; const char *end_state_mass_literal = "m(K^{+}#mu^{+}#mu^{-})"; - const bool retrain_bdt = true; + 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/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_RD.root"); @@ -71,15 +71,27 @@ int new_analysis_bu2kpmumu() sim_chain->SetBranchAddress("K_TRUEID", &Hp_TRUEID); sim_chain->SetBranchAddress("B_BKGCAT", &B_BKGCAT); - 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); - TH1D *h1_B_Mass_sim = new TH1D("h1_B_Mass_sim", TString::Format("B Mass, Simualted (%s_noPID)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + Double_t B_Mass_jpsi_var, B_Mass_psi2s_var, B_Mass_sim_var; + TString B_Mass_jpsi_var_name = "B_Mass_jpsi_var"; + TString B_Mass_psi2s_var_name = "B_Mass_psi2s_var"; + TString B_Mass_sim_var_name = "B_Mass_sim_var"; + + TTree *tree_B_Mass_jpsi = new TTree("tree_B_Mass_jpsi", TString::Format("B Mass, J/#psi Mode (%s)", analysis_name)); + TTree *tree_B_Mass_psi2s = new TTree("tree_B_Mass_psi2s", TString::Format("B Mass, #psi(2S) Mode (%s)", analysis_name)); + TTree *tree_B_Mass_sim = new TTree("tree_B_Mass_sim", TString::Format("B Mass, Simualted (%s_noPID)", analysis_name)); + + tree_B_Mass_jpsi->Branch(B_Mass_jpsi_var_name, &B_Mass_jpsi_var); + tree_B_Mass_psi2s->Branch(B_Mass_psi2s_var_name, &B_Mass_psi2s_var); + tree_B_Mass_sim->Branch(B_Mass_sim_var_name, &B_Mass_sim_var); + + TH1D *h1_B_Mass_unf = new TH1D("h1_B_Mass_unf", TString::Format("B Mass (%s), Unfiltered", data_tree_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TH1D *h1_B_Mass_sim_unf = new TH1D("h1_B_Mass_sim_unf", TString::Format("B Mass, Simualted (%s), Unfiltered", sim_tree_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); - h1_B_Mass_sim->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_unf->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_sim_unf->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); @@ -94,9 +106,11 @@ int new_analysis_bu2kpmumu() { sim_chain->GetEntry(i); Double_t reconstructed_B_Mass = (hp4v_sim->LorentzVector() + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M(); + h1_B_Mass_sim_unf->Fill(reconstructed_B_Mass); if (B_BKGCAT == 0 && TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID && TMath::Abs(Hp_TRUEID) == PID_KAON) { - h1_B_Mass_sim->Fill(reconstructed_B_Mass); + B_Mass_sim_var = reconstructed_B_Mass; + tree_B_Mass_sim->Fill(); } PrintProgress(TString::Format("%s SIG Collection", analysis_name), sim_entries, 10000, i); @@ -117,13 +131,17 @@ int new_analysis_bu2kpmumu() FillHlt1DecisionHistos(hlt1_decision_histos, reconstructed_B_Mass); } + h1_B_Mass_unf->Fill(reconstructed_B_Mass); + if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) { - h1_B_Mass_jpsi->Fill(reconstructed_B_Mass); + B_Mass_jpsi_var = reconstructed_B_Mass; + tree_B_Mass_jpsi->Fill(); } else if (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.) { - h1_B_Mass_psi2s->Fill(reconstructed_B_Mass); + B_Mass_psi2s_var = reconstructed_B_Mass; + tree_B_Mass_psi2s->Fill(); } PrintProgress(TString::Format("%s BDT Evaluation", analysis_name), data_entries, 10000, i); @@ -135,14 +153,12 @@ 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(h1_B_Mass_unf, analysis_name, 0.1); + DrawInDefaultCanvas(h1_B_Mass_sim_unf, 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); - auto roofit_hist_psi2s_fitsum = CreateRooFitHistogramAndFitCB(h1_B_Mass_psi2s, true, true, roofit_hist_sim.shape_parameters); + auto roofit_hist_sim = CreateRooDataSetAndFitCB(tree_B_Mass_sim, B_Mass_sim_var_name, end_state_mass_literal, false, false, ShapeParamters{}); + auto roofit_hist_jpsi_fitsum = CreateRooDataSetAndFitCB(tree_B_Mass_jpsi, B_Mass_jpsi_var_name, end_state_mass_literal, true, true, roofit_hist_sim.shape_parameters); + auto roofit_hist_psi2s_fitsum = CreateRooDataSetAndFitCB(tree_B_Mass_psi2s, B_Mass_psi2s_var_name, end_state_mass_literal, true, true, roofit_hist_sim.shape_parameters); DrawInDefaultCanvas(roofit_hist_jpsi_fitsum, analysis_name); DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name);