From 06a90384eb638e1c0589182dc30c97bd5d400fad Mon Sep 17 00:00:00 2001 From: Marius Pfeiffer Date: Wed, 3 Apr 2024 15:20:59 +0200 Subject: [PATCH] 20240403 --- B0ToHpHmMuMu_results.txt | 66 ++++++++++++-- B0ToKpPimMuMu_results.txt | 64 +++++++++++++ BuToHpMuMu_results.txt | 70 ++++++++++++-- BuToJpsiKplus_results.txt | 37 ++++++++ BuToKpMuMu_results.txt | 64 +++++++++++++ basic_analysis.h | 167 +++++++++++++++++++++++++--------- batch_run_all_anas.sh | 4 + bdt_classification.h | 4 +- constants.h | 2 + mapmc_b02hphmmumu.cpp | 85 +++++++++++++++-- mapmc_bu2hpmumu.cpp | 79 +++++++++++++--- new_analysis_b02hphmmumu.cpp | 93 +++++++++++++++---- new_analysis_b02kppimmumu.cpp | 77 +++++++++++++--- new_analysis_bu2hpmumu.cpp | 90 ++++++++++++++---- new_analysis_bu2jpsikplus.cpp | 82 +++++++++++++++-- new_analysis_bu2kpmumu.cpp | 81 +++++++++++++---- rounder.cpp | 78 ++++++++++++++++ tmva_plotting.cpp | 138 ++++++++++++++++++++++++++++ 18 files changed, 1123 insertions(+), 158 deletions(-) create mode 100644 B0ToKpPimMuMu_results.txt create mode 100644 BuToJpsiKplus_results.txt create mode 100644 BuToKpMuMu_results.txt create mode 100755 batch_run_all_anas.sh create mode 100644 rounder.cpp create mode 100644 tmva_plotting.cpp diff --git a/B0ToHpHmMuMu_results.txt b/B0ToHpHmMuMu_results.txt index 990e360..1dabbae 100644 --- a/B0ToHpHmMuMu_results.txt +++ b/B0ToHpHmMuMu_results.txt @@ -1,6 +1,60 @@ -#### Tue Mar 12 16:38:50 2024 #### -J/Psi Mode: 380.717 +- 23.9177 -Psi(2S) Mode: 14.7761 +- 5.0852 -Mode Yield Ratio: 0.0388112 +- 0.0135776 -Rel Br Frac MuMu: 7.7013 -Rel Br Frac: 0.298897 +#### B0ToHpHmMuMu @ Tue Apr 2 10:44:44 2024 #### +J/Psi Mode: 341 #pm 23 / 185 #pm 23 + Sig/Bkg: 1.85 #pm 0.26 +Psi(2S) Mode: 23 #pm 6 / 43 #pm 11 + Sig/Bkg: 0.54 #pm 0.19 +Mode Yield Ratio: 0.068 #pm 0.018 +Rel Br Frac MuMu: 7.701 #pm 0.804 +Rel Br Frac: 0.522 #pm 0.149 + + +Fitted Parameters: Simulation + +#mu = 5278.85 \pm 0.13 +#sigma_{LR} = 15.87 \pm 0.16 +#alpha_{L} = 1.72 \pm 0.06 +n_{L} = 3.34 \pm 0.30 +#alpha_{R} = 1.76 \pm 0.07 +n_{R} = 5.67 \pm 0.73 + +Fitted Parameters: J/PSI + +N_{Sig} = 341 \pm 23 +N_{Bkg} = 185 \pm 23 +N_{Bkg,2#sigma} = 97 \pm 6 +N_{Sig}/N_{Bkg} = 1.85 \pm 0.26 +N_{Sig}/N_{Bkg,2#sigma} = 3.53 \pm 0.32 +#lambda = -0.00217 \pm 0.00066 +#mu = 5270.03 \pm 1.92 +#sigma_{LR} = 26.26 \pm 1.98 +#alpha_{L} = 1.72 (c) +n_{L} = 3.34 (c) +#alpha_{R} = 1.76 (c) +n_{R} = 5.67 (c) + +Fitted Parameters: PSI(2S) + +N_{Sig} = 23 \pm 6 +N_{Bkg} = 43 \pm 10 +N_{Bkg,2#sigma} = 12 \pm 2 +N_{Sig}/N_{Bkg} = 0.54 \pm 0.19 +N_{Sig}/N_{Bkg,2#sigma} = 1.93 \pm 0.60 +#lambda = -0.00202 \pm 0.00120 +#mu = 5267.35 \pm 7.06 +#sigma_{LR} = 26.26 (c) +#alpha_{L} = 1.72 (c) +n_{L} = 3.34 (c) +#alpha_{R} = 1.76 (c) +n_{R} = 5.67 (c) + +# J/psi +\begin{tabular}{c|c} +$N_{Sig}$ & $N_{Bkg}$\\\hline +$341 \pm 23$ & $185 \pm 23$ +\end{tabular} + +# psi(2S) +\begin{tabular}{c|c} +$N_{Sig}$ & $N_{Bkg}$\\\hline +$23 \pm 6$ & $43 \pm 10$ +\end{tabular} diff --git a/B0ToKpPimMuMu_results.txt b/B0ToKpPimMuMu_results.txt new file mode 100644 index 0000000..7f24eca --- /dev/null +++ b/B0ToKpPimMuMu_results.txt @@ -0,0 +1,64 @@ +#### B0ToKpPimMuMu @ Tue Apr 2 10:46:02 2024 #### +J/Psi Mode: 305 #pm 21 / 175 #pm 21 + Sig/Bkg: 1.74 #pm 0.24 +Psi(2S) Mode: 26 #pm 7 / 68 #pm 19 + Sig/Bkg: 0.39 #pm 0.15 +Mode Yield Ratio: 0.086 #pm 0.022 +Rel Br Frac MuMu: 7.701 #pm 0.804 +Rel Br Frac: 0.662 #pm 0.179 + +Params from Sim: +aL: 1.74 +nL: 3.27 +aR: 1.82 +nR: 5.82 +S: 25.50 + +#mu = 5278.97 \pm 0.13 +#sigma_{LR} = 15.77 \pm 0.15 +#alpha_{L} = 1.74 \pm 0.06 +n_{L} = 3.27 \pm 0.30 +#alpha_{R} = 1.82 \pm 0.07 +n_{R} = 5.82 \pm 0.81 + +Fitted Parameters: J/PSI + +N_{Sig} = 305 \pm 21 +N_{Bkg} = 175 \pm 20 +N_{Bkg,2#sigma} = 98 \pm 6 +N_{Sig}/N_{Bkg} = 1.74 \pm 0.23 +N_{Sig}/N_{Bkg,2#sigma} = 3.12 \pm 0.28 +#lambda = -0.00332 \pm 0.00062 +#mu = 5269.84 \pm 1.85 +#sigma_{LR} = 25.50 \pm 1.77 +#alpha_{L} = 1.74 (c) +n_{L} = 3.27 (c) +#alpha_{R} = 1.82 (c) +n_{R} = 5.82 (c) + +Fitted Parameters: PSI(2S) + +N_{Sig} = 26 \pm 6 +N_{Bkg} = 68 \pm 19 +N_{Bkg,2#sigma} = 11 \pm 2 +N_{Sig}/N_{Bkg} = 0.39 \pm 0.14 +N_{Sig}/N_{Bkg,2#sigma} = 2.44 \pm 0.78 +#lambda = -0.00003 \pm 0.00398 +#mu = 5276.41 \pm 6.20 +#sigma_{LR} = 25.50 (c) +#alpha_{L} = 1.74 (c) +n_{L} = 3.27 (c) +#alpha_{R} = 1.82 (c) +n_{R} = 5.82 (c) + +# J/psi +\begin{tabular}{c|c} +$N_{Sig}$ & $N_{Bkg}$\\\hline +$305 \pm 21$ & $175 \pm 20$ +\end{tabular} + +# psi(2S) +\begin{tabular}{c|c} +$N_{Sig}$ & $N_{Bkg}$\\\hline +$26 \pm 6$ & $68 \pm 19$ +\end{tabular} diff --git a/BuToHpMuMu_results.txt b/BuToHpMuMu_results.txt index 8db82a9..368349e 100644 --- a/BuToHpMuMu_results.txt +++ b/BuToHpMuMu_results.txt @@ -1,6 +1,64 @@ -#### Tue Mar 12 16:35:42 2024 #### -J/Psi Mode: 950.052 +- 36.9322 -Psi(2S) Mode: 65.8457 +- 11.1206 -Mode Yield Ratio: 0.0693074 +- 0.0120113 -Rel Br Frac MuMu: 7.7013 -Rel Br Frac: 0.533757 +#### BuToHpMuMu @ Tue Apr 2 10:44:14 2024 #### +J/Psi Mode: 897 #pm 34 / 402 #pm 30 + Sig/Bkg: 2.23 #pm 0.19 +Psi(2S) Mode: 69 #pm 12 / 233 #pm 26 + Sig/Bkg: 0.30 #pm 0.06 +Mode Yield Ratio: 0.077 #pm 0.013 +Rel Br Frac MuMu: 7.701 #pm 0.804 +Rel Br Frac: 0.594 #pm 0.116 + +Params from Sim: +aL: 1.76 +nL: 4.02 +aR: 1.93 +nR: 4.33 +S: 23.48 + +#mu = 5278.16 \pm 0.05 +#sigma_{LR} = 16.43 \pm 0.05 +#alpha_{L} = 1.76 \pm 0.02 +n_{L} = 4.02 \pm 0.16 +#alpha_{R} = 1.93 \pm 0.03 +n_{R} = 4.33 \pm 0.19 + +Fitted Parameters: J/PSI + +N_{Sig} = 897 \pm 34 +N_{Bkg} = 402 \pm 30 +N_{Bkg,2#sigma} = 246 \pm 9 +N_{Sig}/N_{Bkg} = 2.23 \pm 0.19 +N_{Sig}/N_{Bkg,2#sigma} = 3.64 \pm 0.19 +#lambda = -0.00350 \pm 0.00041 +#mu = 5271.03 \pm 0.96 +#sigma_{LR} = 23.48 \pm 0.90 +#alpha_{L} = 1.76 (c) +n_{L} = 4.02 (c) +#alpha_{R} = 1.93 (c) +n_{R} = 4.33 (c) + +Fitted Parameters: PSI(2S) + +N_{Sig} = 69 \pm 11 +N_{Bkg} = 233 \pm 26 +N_{Bkg,2#sigma} = 44 \pm 4 +N_{Sig}/N_{Bkg} = 0.30 \pm 0.06 +N_{Sig}/N_{Bkg,2#sigma} = 1.56 \pm 0.29 +#lambda = -0.00142 \pm 0.00056 +#mu = 5263.32 \pm 4.14 +#sigma_{LR} = 23.48 (c) +#alpha_{L} = 1.76 (c) +n_{L} = 4.02 (c) +#alpha_{R} = 1.93 (c) +n_{R} = 4.33 (c) + +# J/psi +\begin{tabular}{c|c} +$N_{Sig}$ & $N_{Bkg}$\\\hline +$897 \pm 34$ & $402 \pm 30$ +\end{tabular} + +# psi(2S) +\begin{tabular}{c|c} +$N_{Sig}$ & $N_{Bkg}$\\\hline +$69 \pm 11$ & $233 \pm 26$ +\end{tabular} diff --git a/BuToJpsiKplus_results.txt b/BuToJpsiKplus_results.txt new file mode 100644 index 0000000..fad64e1 --- /dev/null +++ b/BuToJpsiKplus_results.txt @@ -0,0 +1,37 @@ +#### BuToJpsiKplus @ Wed Apr 3 13:28:02 2024 #### +J/Psi Mode: 892 #pm 61 / 1586 #pm 66 + Sig/Bkg: 0.56 #pm 0.05 +Params from Sim: +aL: 1.71 +nL: 4.54 +aR: 1.89 +nR: 5.30 +S: 26.83 + +#mu = 5278.10 \pm 0.05 +#sigma_{LR} = 16.27 \pm 0.06 +#alpha_{L} = 1.71 \pm 0.02 +n_{L} = 4.54 \pm 0.25 +#alpha_{R} = 1.89 \pm 0.03 +n_{R} = 5.30 \pm 0.33 + +Fitted Parameters: J/PSI + +N_{Sig} = 892 \pm 61 +N_{Bkg} = 1586 \pm 66 +N_{Bkg,2#sigma} = 916 \pm 33 +N_{Sig}/N_{Bkg} = 0.56 \pm 0.04 +N_{Sig}/N_{Bkg,2#sigma} = 0.97 \pm 0.08 +#lambda = -0.00140 \pm 0.00031 +#mu = 5271.73 \pm 1.57 +#sigma_{LR} = 26.83 \pm 2.04 +#alpha_{L} = 1.71 (c) +n_{L} = 4.54 (c) +#alpha_{R} = 1.89 (c) +n_{R} = 5.30 (c) + +# J/psi +\begin{tabular}{c|c} +$N_{Sig}$ & $N_{Bkg}$\\\hline +$892 \pm 61$ & $1586 \pm 66$ +\end{tabular} diff --git a/BuToKpMuMu_results.txt b/BuToKpMuMu_results.txt new file mode 100644 index 0000000..1d3f300 --- /dev/null +++ b/BuToKpMuMu_results.txt @@ -0,0 +1,64 @@ +#### BuToKpMuMu @ Tue Apr 2 10:45:35 2024 #### +J/Psi Mode: 815 #pm 35 / 743 #pm 45 + Sig/Bkg: 1.10 #pm 0.09 +Psi(2S) Mode: 69 #pm 12 / 290 #pm 30 + Sig/Bkg: 0.24 #pm 0.05 +Mode Yield Ratio: 0.085 #pm 0.015 +Rel Br Frac MuMu: 7.701 #pm 0.804 +Rel Br Frac: 0.654 #pm 0.132 + +Params from Sim: +aL: 1.77 +nL: 3.86 +aR: 1.96 +nR: 4.43 +S: 24.12 + +#mu = 5278.10 \pm 0.05 +#sigma_{LR} = 16.35 \pm 0.05 +#alpha_{L} = 1.77 \pm 0.02 +n_{L} = 3.86 \pm 0.16 +#alpha_{R} = 1.96 \pm 0.03 +n_{R} = 4.43 \pm 0.20 + +Fitted Parameters: J/PSI + +N_{Sig} = 815 \pm 34 +N_{Bkg} = 743 \pm 44 +N_{Bkg,2#sigma} = 259 \pm 9 +N_{Sig}/N_{Bkg} = 1.10 \pm 0.08 +N_{Sig}/N_{Bkg,2#sigma} = 3.15 \pm 0.18 +#lambda = -0.00207 \pm 0.00031 +#mu = 5271.47 \pm 1.08 +#sigma_{LR} = 24.12 \pm 1.02 +#alpha_{L} = 1.77 (c) +n_{L} = 3.86 (c) +#alpha_{R} = 1.96 (c) +n_{R} = 4.43 (c) + +Fitted Parameters: PSI(2S) + +N_{Sig} = 69 \pm 12 +N_{Bkg} = 290 \pm 30 +N_{Bkg,2#sigma} = 53 \pm 5 +N_{Sig}/N_{Bkg} = 0.24 \pm 0.05 +N_{Sig}/N_{Bkg,2#sigma} = 1.32 \pm 0.25 +#lambda = -0.00130 \pm 0.00051 +#mu = 5268.10 \pm 4.25 +#sigma_{LR} = 24.12 (c) +#alpha_{L} = 1.77 (c) +n_{L} = 3.86 (c) +#alpha_{R} = 1.96 (c) +n_{R} = 4.43 (c) + +# J/psi +\begin{tabular}{c|c} +$N_{Sig}$ & $N_{Bkg}$\\\hline +$815 \pm 34$ & $743 \pm 44$ +\end{tabular} + +# psi(2S) +\begin{tabular}{c|c} +$N_{Sig}$ & $N_{Bkg}$\\\hline +$69 \pm 12$ & $290 \pm 30$ +\end{tabular} diff --git a/basic_analysis.h b/basic_analysis.h index c75e8e0..b40c664 100644 --- a/basic_analysis.h +++ b/basic_analysis.h @@ -21,6 +21,50 @@ std::pair DivWithErr(double x, double dx, double y, double dy) { return std::make_pair(x / y, TMath::Sqrt((TMath::Sq(err_x) + TMath::Sq(err_y)))); } +std::pair MultWithErr(double x, double dx, double y, double dy) { + double err_x = y * dx; + double err_y = - x * dy; + return std::make_pair(x * y, TMath::Sqrt((TMath::Sq(err_x) + TMath::Sq(err_y)))); +} + +double RoundUp(double x, int n = 0) { + if(std::abs(x) > std::numeric_limits::max() / n) // v * p would overflow + throw std::overflow_error("rounding would overflow"); + const double multiplier = std::pow(10.0, n); + return std::ceil(x * multiplier) / multiplier; +} + +double RoundDown(double x, int n = 0) { + if(std::abs(x) > std::numeric_limits::max() / n) // v * p would overflow + throw std::overflow_error("rounding would overflow"); + const double multiplier = std::pow(10.0, n); + return std::floor(x * multiplier) / multiplier; +} + +double Round(double x, int n = 0) { + if(std::abs(x) > std::numeric_limits::max() / n) // v * p would overflow + throw std::overflow_error("rounding would overflow"); + const double multiplier = std::pow(10.0, n); + return std::round(x * multiplier) / multiplier; +} + +std::string ErrToStr(double val, double err, int n) { + return TString::Format("%.*f #pm %.*f", n, Round(val, n), n, RoundUp(err, n)).Data(); +} + +std::string ErrToStr(std::pair val_err, int n) { + return ErrToStr(val_err.first, val_err.second, n); +} + +std::pair CalcSignificance(double sig_val, double sig_err, double bkg_val, double bkg_err) { + double significance_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 significance_err = TMath::Sqrt(TMath::Sq(err_prop_sig * sig_err) + TMath::Sq(err_prop_bkg * bkg_err)); + + return std::make_pair(significance_val, significance_err); +} + class FourVect { private: @@ -84,14 +128,18 @@ struct FittedParam bool at_limit; bool constant; int decimals; + bool print_const; + bool ignore_print; - FittedParam(RooRealVar var, int decimals) + FittedParam(RooRealVar var, int decimals, bool print_const = false, bool ignore_print = false) { this->title = var.GetTitle(); this->name = var.GetName(); double val = var.getVal(); this->value = val; this->constant = var.isConstant(); + this->print_const = print_const; + this->ignore_print = ignore_print; if (!this->constant) { this->error = var.getError(); @@ -104,7 +152,7 @@ struct FittedParam this->decimals = decimals; } - FittedParam(std::string name, std::string title, double value, double err, int decimals) + FittedParam(std::string name, std::string title, double value, double err, int decimals, bool ignore_print = false) { this->title = title; this->name = name; @@ -113,17 +161,23 @@ struct FittedParam this->decimals = decimals; this->constant = false; this->at_limit = false; + this->ignore_print = ignore_print; } - std::string ToString() const + std::string ToString(bool latex = false) const { if (this->constant) { - return TString::Format("%s = %.*f", this->title.c_str(), this->decimals, this->value).Data(); + return TString::Format("%s = %.*f (c)", this->title.c_str(), this->decimals, this->value).Data(); } else { - return TString::Format("%s = %.*f #pm %.*f", this->title.c_str(), this->decimals, this->value, this->decimals, this->error).Data(); + if (this->error == 0) { + return TString::Format("%s = %.*f", this->title.c_str(), this->decimals, this->value).Data(); + } + else { + return TString::Format("%s = %.*f %s %.*f", this->title.c_str(), this->decimals, this->value, (latex ? "\\pm" : "#pm"), this->decimals, this->error).Data(); + } } } }; @@ -135,6 +189,10 @@ struct ShapeParamters double alpha_right; double n_right; double sigma_lr; + + std::string ToString() { + return TString::Format("aL: %.2f\nnL: %.2f\naR: %.2f\nnR: %.2f\nS: %.2f\n", this->alpha_left, this->n_left, this->alpha_right, this->n_right, this->sigma_lr).Data(); + } }; struct RooFitSummary @@ -160,21 +218,34 @@ 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 = "") +void DrawInDefaultCanvasStacked(std::vector histograms, std::vector colors, std::vector fill_style, const char *folder, const char* title, double margin_left = 0.12, 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(); + Double_t max = 0; + for (size_t i = 0; i < histograms.size(); i++) + { + if (histograms[i]->GetEntries() > max) { + max = histograms[i]->GetEntries(); + } + } + + TString stack_name = TString::Format("%s_stack", histograms[0]->GetName()); + TString stack_title = TString::Format("%s;%s;#Events normed to 1/%.0f", title, histograms[0]->GetXaxis()->GetTitle(), max); + + std::string drwopt_2 = std::string(option).empty() ? "HIST" : TString::Format("%s HIST", option).Data(); + auto stack = new THStack(stack_name, stack_title); 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->Scale(scaling_factor / max); hist_clone->SetLineColor(colors[i]); - hist_clone->SetMaximum(scaling_factor + (scaling_factor * 0.05)); + // hist_clone->SetMaximum(scaling_factor + (scaling_factor * 0.05)); + // hist_clone->SetMaximum(max + (max * 0.05)); hist_clone->SetMinimum(0.); hist_clone->SetStats(0); if (fill_style[i] != 0) @@ -185,15 +256,17 @@ void DrawInDefaultCanvasStacked(std::vector histograms, std::vector 0) { - hist_clone->Draw(drwopt_2.c_str()); + stack->Add(hist_clone, drwopt_2.c_str()); } else { - hist_clone->Draw(drwopt_2.c_str()); + stack->Add(hist_clone, drwopt_2.c_str()); } } - c->BuildLegend(0.55, 0.79, 0.87, 0.89); + stack->Draw("NOSTACK"); + + c->BuildLegend(0.55, 0.70, 0.87, 0.89); c->Draw(); c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data()); @@ -259,20 +332,8 @@ void DrawInDefaultCanvas(RooFitSummary fitSummary, const char *folder) Double_t fit_title_size = (pull_title_size * dwPad_area) / upPad_area; Double_t fit_label_size = (pull_label_size * 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->SetLabelOffset(0.02); @@ -298,7 +359,7 @@ void DrawInDefaultCanvas(RooFitSummary fitSummary, const char *folder) for (const auto &par : fitSummary.fitted_params) { - if (!par.constant) + if ((!par.ignore_print) && ((!par.constant) || par.print_const)) { leg1->AddEntry((TObject *)0, par.ToString().c_str(), ""); } @@ -373,7 +434,7 @@ void PrintProgress(const char *title, unsigned int total, unsigned int every, un RooPlot *CreateRooFitHistogram(TH1D *hist) { RooRealVar roo_var_mass("var_mass", "B Mass Variable", B_MASS_VAR_MIN, B_MASS_VAR_MAX); - roo_var_mass.setRange("fitting_range", B_MASS_VAR_MIN, B_MASS_VAR_MAX); + roo_var_mass.setRange("fitting_range", B_MASS_FIT_MIN, B_MASS_FIT_MAX); RooDataHist roohist_B_M("roohist_B_M", "B Mass Histogram", roo_var_mass, RooFit::Import(*hist)); @@ -391,13 +452,16 @@ RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString return TString::Format("%s_%s", text, name); }; - Double_t fitRangeUp = fit_low == 0 ? B_MASS_VAR_MIN : fit_low; - Double_t fitRangeLow = fit_up == 0 ? B_MASS_VAR_MAX : fit_up; + Double_t fitRangeUp = fit_low == 0 ? B_MASS_FIT_MIN : fit_low; + Double_t fitRangeLow = fit_up == 0 ? B_MASS_FIT_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, fitRangeUp, fitRangeLow); + const char *draw_range_name = "draw_range"; + roo_var_mass.setRange(draw_range_name, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + 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)); @@ -477,30 +541,45 @@ RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString 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(roodataset_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_pdf.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range(draw_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(draw_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(draw_range_name), RooFit::Name(signal_name)); - fitted_params.push_back(FittedParam(var_mass_nsig, 2)); - fitted_params.push_back(FittedParam(var_mass_nbkg, 2)); + fitted_params.push_back(FittedParam(var_mass_nsig, 0)); + fitted_params.push_back(FittedParam(var_mass_nbkg, 0)); + double mu_val = var_mass_x0.getVal(); + double sigma_val = var_mass_sigmaLR.getVal(); sig_val = var_mass_nsig.getVal(); sig_err = var_mass_nsig.getError(); bkg_val = var_mass_nbkg.getVal(); bkg_err = var_mass_nbkg.getError(); - double significance_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 significance_err = TMath::Sqrt(TMath::Sq(err_prop_sig * sig_err) + TMath::Sq(err_prop_bkg * bkg_err)); + double sig_bkg_sum_val = sig_val + bkg_val; + double sig_bkg_sum_err = TMath::Sqrt(TMath::Sq(sig_err) + TMath::Sq(bkg_err)); + + const int sigma_space = 2; + double int_lo = mu_val - sigma_space * sigma_val; + double int_up = mu_val + sigma_space * sigma_val; + roo_var_mass.setRange("peak_integral_range", int_lo, int_up); + std::cout << "#### Integral Over: " << int_lo << " -> " << int_up << std::endl; + auto bkg_int = bkg_exp.createIntegral(roo_var_mass, RooFit::NormSet(roo_var_mass), RooFit::Range("peak_integral_range")); + + auto bkg_yield_under_sig = std::make_pair(bkg_int->getVal() * sig_bkg_sum_val, bkg_int->getVal() * sig_bkg_sum_err); + + fitted_params.push_back(FittedParam("bkg_int_sig", TString::Format("N_{Bkg,%d#sigma}", sigma_space).Data(), bkg_yield_under_sig.first, bkg_yield_under_sig.second, 0)); + + auto significance = CalcSignificance(sig_val, sig_err, bkg_val, bkg_err); - auto sig_over_bkg = DivWithErr(sig_val, sig_err, bkg_val, bkg_err); + auto sig_over_bkg_3sig = DivWithErr(sig_val, sig_err, bkg_yield_under_sig.first, bkg_yield_under_sig.second); + auto sig_over_bkg_full = DivWithErr(sig_val, sig_err, bkg_val, bkg_err); - fitted_params.push_back(FittedParam("significance", "N_{Sig}/#sqrt{N_{Sig} + N_{Bkg}}", significance_val, significance_err, 2)); - fitted_params.push_back(FittedParam("sig_over_bkg", "N_{Sig}/N_{Bkg}", sig_over_bkg.first, sig_over_bkg.second, 2)); + // fitted_params.push_back(FittedParam("significance", "N_{Sig}/#sqrt{N_{Sig} + N_{Bkg}}", significance.first, significance.second, 2)); + fitted_params.push_back(FittedParam("sig_over_bkg", "N_{Sig}/N_{Bkg}", sig_over_bkg_full.first, sig_over_bkg_full.second, 2, true)); + fitted_params.push_back(FittedParam("sig_over_bkg", TString::Format("N_{Sig}/N_{Bkg,%d#sigma}", sigma_space).Data(), sig_over_bkg_3sig.first, sig_over_bkg_3sig.second, 2, true)); fitted_params.push_back(FittedParam(var_mass_bkg_c, 5)); @@ -508,14 +587,14 @@ RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString } else { - RooFitResult *fitres = sig_cb.fitTo(roodataset_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)); + sig_cb.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range(draw_range_name), RooFit::Name(signal_name)); pull_compare_name = signal_name; } fitted_params.push_back(FittedParam(var_mass_x0, 2)); - fitted_params.push_back(FittedParam(var_mass_sigmaLR, 2)); + fitted_params.push_back(FittedParam(var_mass_sigmaLR, 2, true)); fitted_params.push_back(FittedParam(var_mass_alphaL, 2)); fitted_params.push_back(FittedParam(var_mass_nL, 2)); fitted_params.push_back(FittedParam(var_mass_alphaR, 2)); diff --git a/batch_run_all_anas.sh b/batch_run_all_anas.sh new file mode 100755 index 0000000..7a4d5be --- /dev/null +++ b/batch_run_all_anas.sh @@ -0,0 +1,4 @@ +root -b -q ./new_analysis_bu2hpmumu.cpp +root -b -q ./new_analysis_b02hphmmumu.cpp +root -b -q ./new_analysis_bu2kpmumu.cpp +root -b -q ./new_analysis_b02kppimmumu.cpp \ No newline at end of file diff --git a/bdt_classification.h b/bdt_classification.h index 7e717b2..d5b9c3d 100644 --- a/bdt_classification.h +++ b/bdt_classification.h @@ -205,7 +205,7 @@ void TrainBDT(std::vector vars, const char* unique_id, TTree *sig_tree, TT TString outfile_name = TString::Format("%s_tmva_out.root", unique_id); TFile *output_file = TFile::Open(outfile_name, "RECREATE"); - TString factory_options("V:Silent:Color:DrawProgressBar:Transformations=I;D;P;G,D:AnalysisType=Auto"); + TString factory_options("V:!Silent:Color:DrawProgressBar:Transformations=I;D;P;G,D:AnalysisType=Auto"); TMVA::Factory *factory = new TMVA::Factory(TString::Format("%s_factory", unique_id), output_file, factory_options); TMVA::DataLoader *data_loader = new TMVA::DataLoader(TString::Format("%s_dataloader", unique_id)); @@ -227,7 +227,7 @@ void TrainBDT(std::vector vars, const char* unique_id, TTree *sig_tree, TT data_loader->AddSignalTree(sig_tree, signal_weight); data_loader->AddBackgroundTree(bkg_tree, background_weight); - data_loader->PrepareTrainingAndTestTree("", "", "nTrain_Signal=0:nTrain_Background=0:SplitMode=Random:NormMode=EqualNumEvents:!V"); + data_loader->PrepareTrainingAndTestTree("", "", "nTrain_Signal=0:nTrain_Background=0:SplitMode=Random:NormMode=NumEvents:V"); factory->BookMethod(data_loader, TMVA::Types::kBDT, "BDT", "!H:!V:NTrees=600:MinNodeSize=2.5%:CreateMVAPdfs:MaxDepth=3:BoostType=AdaBoost:AdaBoostBeta=0.5:UseBaggedBoost:BaggedSampleFraction=0.5:SeparationType=GiniIndex:nCuts=20"); diff --git a/constants.h b/constants.h index 0aa09ed..3daf1d8 100644 --- a/constants.h +++ b/constants.h @@ -5,6 +5,8 @@ const Double_t B_MASS_VAR_MIN = 5100.; const Double_t B_MASS_VAR_MAX = 6000.; +const Double_t B_MASS_FIT_MIN = 5100.; +const Double_t B_MASS_FIT_MAX = 5600.; const Int_t B_MASS_HIST_BINS = 70; const Double_t K_MASS = 493.677; diff --git a/mapmc_b02hphmmumu.cpp b/mapmc_b02hphmmumu.cpp index e281d61..d091db2 100644 --- a/mapmc_b02hphmmumu.cpp +++ b/mapmc_b02hphmmumu.cpp @@ -44,7 +44,7 @@ int mapmc_b02hphmmumu() const char *sim_tree_name = "B0ToHpHmMuMu_noPID"; 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_fullstream_v0r0p6671378_B0ToKpPimMuMu_11144002_magdown.root"); + sim_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/rd_btoxll_simulation_fullstream_v0r0p6919194_B0ToKpPimMuMu_11144002_magdown.root"); FourVect *l14v_sim = FourVect::Init(sim_chain, "L1"); FourVect *l24v_sim = FourVect::Init(sim_chain, "L2"); @@ -59,39 +59,77 @@ int mapmc_b02hphmmumu() sim_chain->SetBranchAddress("Hm_TRUEID", &Hm_TRUEID); sim_chain->SetBranchAddress("B0_BKGCAT", &B_BKGCAT); - Float_t B0_PT_in, B0_BPVFDCHI2_in, B0_BPVDIRA_in, Jpsi_BPVIPCHI2_in, Jpsi_PT_in, - Hp_BPVIPCHI2_in, Hp_PT_in, Hm_BPVIPCHI2_in, Hm_PT_in, L1_BPVIPCHI2_in, L2_BPVIPCHI2_in; - Double_t Hp_PROBNN_K_in, Hm_PROBNN_K_in; + Float_t B0_PT_in, B0_BPVFDCHI2_in, B0_BPVDIRA_in, + Jpsi_BPVIPCHI2_in, Jpsi_PT_in, Jpsi_BPVDIRA_in, Jpsi_BPVFDCHI2_in, + Res_BPVIPCHI2_in, Res_PT_in, Res_BPVFDCHI2_in, + Hp_BPVIPCHI2_in, Hp_PT_in, + Hm_BPVIPCHI2_in, Hm_PT_in, + L1_BPVIPCHI2_in, L1_PT_in, + L2_BPVIPCHI2_in, L2_PT_in; + Double_t B0_CHI2_in, B0_CHI2DOF_in, + Jpsi_CHI2_in, Jpsi_CHI2DOF_in, + Hp_PROBNN_K_in, Hp_PID_K_in, + Hm_PROBNN_K_in, Hm_PID_K_in; sim_chain->SetBranchAddress("B0_PT", &B0_PT_in); sim_chain->SetBranchAddress("B0_BPVFDCHI2", &B0_BPVFDCHI2_in); sim_chain->SetBranchAddress("B0_BPVDIRA", &B0_BPVDIRA_in); + sim_chain->SetBranchAddress("B0_CHI2", &B0_CHI2_in); + sim_chain->SetBranchAddress("B0_CHI2VXNDOF", &B0_CHI2DOF_in); + sim_chain->SetBranchAddress("Jpsi_BPVIPCHI2", &Jpsi_BPVIPCHI2_in); + sim_chain->SetBranchAddress("Jpsi_CHI2", &Jpsi_CHI2_in); + sim_chain->SetBranchAddress("Jpsi_CHI2DOF", &Jpsi_CHI2DOF_in); sim_chain->SetBranchAddress("Jpsi_PT", &Jpsi_PT_in); + sim_chain->SetBranchAddress("Jpsi_BPVFDCHI2", &Jpsi_BPVFDCHI2_in); + + sim_chain->SetBranchAddress("Res_BPVIPCHI2", &Res_BPVIPCHI2_in); + sim_chain->SetBranchAddress("Res_PT", &Res_PT_in); + sim_chain->SetBranchAddress("Res_BPVFDCHI2", &Res_BPVFDCHI2_in); + sim_chain->SetBranchAddress("Hp_BPVIPCHI2", &Hp_BPVIPCHI2_in); sim_chain->SetBranchAddress("Hp_PT", &Hp_PT_in); + sim_chain->SetBranchAddress("Hp_PROBNN_K", &Hp_PROBNN_K_in); + sim_chain->SetBranchAddress("Hp_PID_K", &Hp_PID_K_in); + sim_chain->SetBranchAddress("Hm_BPVIPCHI2", &Hm_BPVIPCHI2_in); sim_chain->SetBranchAddress("Hm_PT", &Hm_PT_in); + sim_chain->SetBranchAddress("Hm_PROBNN_K", &Hm_PROBNN_K_in); + sim_chain->SetBranchAddress("Hm_PID_K", &Hm_PID_K_in); + sim_chain->SetBranchAddress("L1_BPVIPCHI2", &L1_BPVIPCHI2_in); + sim_chain->SetBranchAddress("L1_PT", &L1_PT_in); + sim_chain->SetBranchAddress("L2_BPVIPCHI2", &L2_BPVIPCHI2_in); - sim_chain->SetBranchAddress("Hp_PROBNN_K", &Hp_PROBNN_K_in); - sim_chain->SetBranchAddress("Hm_PROBNN_K", &Hm_PROBNN_K_in); + sim_chain->SetBranchAddress("L2_PT", &L2_PT_in); TTree *output_tree = new TTree("DecayTree", "DecayTree"); - Float_t B0_PT_out, B0_BPVFDCHI2_out, B0_BPVDIRA_out, Jpsi_BPVIPCHI2_out, Jpsi_PT_out, - Kplus_BPVIPCHI2_out, Kplus_PT_out, piminus_BPVIPCHI2_out, piminus_PT_out, muplus_BPVIPCHI2_out, muminus_BPVIPCHI2_out, + Float_t B0_PT_out, B0_BPVFDCHI2_out, B0_BPVDIRA_out, + Jpsi_BPVIPCHI2_out, Jpsi_PT_out, Jpsi_BPVDIRA_out, Jpsi_BPVFDCHI2_out, + Kst0_BPVIPCHI2_out, Kst0_PT_out, Kst0_BPVFDCHI2_out, + Kplus_BPVIPCHI2_out, Kplus_PT_out, + piminus_BPVIPCHI2_out, piminus_PT_out, muplus_BPVIPCHI2_out, muplus_PT_out, + muminus_BPVIPCHI2_out, muminus_PT_out, B0_PX_out, B0_PY_out, B0_PZ_out, B0_ENERGY_out, muplus_PX_out, muplus_PY_out, muplus_PZ_out, muplus_ENERGY_out, muminus_PX_out, muminus_PY_out, muminus_PZ_out, muminus_ENERGY_out, Kplus_PX_out, Kplus_PY_out, Kplus_PZ_out, Kplus_ENERGY_out, piminus_PX_out, piminus_PY_out, piminus_PZ_out, piminus_ENERGY_out; - Double_t Kplus_PROBNN_K_out, piminus_PROBNN_K_out, B0_M_out, muplus_M_out, muminus_M_out, Kplus_M_out, piminus_M_out; + + Double_t B0_M_out, B0_CHI2_out, B0_CHI2DOF_out, + Jpsi_CHI2_out, Jpsi_CHI2DOF_out, + muplus_M_out, + muminus_M_out, + Kplus_M_out, Kplus_PROBNN_K_out, Kplus_PID_K_out, + piminus_M_out, piminus_PROBNN_K_out; Int_t muplus_ID_out, muminus_ID_out, Kplus_ID_out, piminus_ID_out; output_tree->Branch("B0_PT", &B0_PT_out, "B0_PT/F"); output_tree->Branch("B0_BPVFDCHI2", &B0_BPVFDCHI2_out, "B0_BPVFDCHI2/F"); output_tree->Branch("B0_BPVDIRA", &B0_BPVDIRA_out, "B0_BPVDIRA/F"); + output_tree->Branch("B0_CHI2", &B0_CHI2_out, "B0_CHI2/D"); + output_tree->Branch("B0_CHI2DOF", &B0_CHI2DOF_out, "B0_CHI2DOF/D"); output_tree->Branch("B0_PX", &B0_PX_out, "0_PX/F"); output_tree->Branch("B0_PY", &B0_PY_out, "0_PY/F"); output_tree->Branch("B0_PZ", &B0_PZ_out, "0_PZ/F"); @@ -100,9 +138,18 @@ int mapmc_b02hphmmumu() output_tree->Branch("Jpsi_BPVIPCHI2", &Jpsi_BPVIPCHI2_out, "Jpsi_BPVIPCHI2/F"); output_tree->Branch("Jpsi_PT", &Jpsi_PT_out, "Jpsi_PT/F"); + output_tree->Branch("Jpsi_BPVDIRA", &Jpsi_BPVDIRA_out, "Jpsi_BPVDIRA/F"); + output_tree->Branch("Jpsi_CHI2", &Jpsi_CHI2_out, "Jpsi_CHI2/D"); + output_tree->Branch("Jpsi_CHI2DOF", &Jpsi_CHI2DOF_out, "Jpsi_CHI2DOF/D"); + output_tree->Branch("Jpsi_BPVFDCHI2", &Jpsi_BPVFDCHI2_out, "Jpsi_BPVFDCHI2/F"); + + output_tree->Branch("Kst0_BPVIPCHI2", &Kst0_BPVIPCHI2_out, "Kst0_BPVIPCHI2/F"); + output_tree->Branch("Kst0_PT", &Kst0_PT_out, "Kst0_PT/F"); + output_tree->Branch("Kst0_BPVFDCHI2", &Kst0_BPVFDCHI2_out, "Kst0_BPVFDCHI2/F"); output_tree->Branch("Kplus_BPVIPCHI2", &Kplus_BPVIPCHI2_out, "Kplus_BPVIPCHI2/F"); output_tree->Branch("Kplus_PROBNN_K", &Kplus_PROBNN_K_out, "Kplus_PROBNN_K/D"); + output_tree->Branch("Kplus_PID_K", &Kplus_PID_K_out, "Kplus_PID_K/D"); output_tree->Branch("Kplus_PT", &Kplus_PT_out, "Kplus_PT/F"); output_tree->Branch("Kplus_PX", &Kplus_PX_out, "Kplus_PX/F"); output_tree->Branch("Kplus_PY", &Kplus_PY_out, "Kplus_PY/F"); @@ -122,6 +169,7 @@ int mapmc_b02hphmmumu() output_tree->Branch("piminus_ID", &piminus_ID_out, "piminus_ID/I"); output_tree->Branch("muminus_BPVIPCHI2", &muminus_BPVIPCHI2_out, "muminus_BPVIPCHI2/F"); + output_tree->Branch("muminus_PT", &muminus_PT_out, "muminus_PT/F"); output_tree->Branch("muminus_PX", &muminus_PX_out, "muminus_PX/F"); output_tree->Branch("muminus_PY", &muminus_PY_out, "muminus_PY/F"); output_tree->Branch("muminus_PZ", &muminus_PZ_out, "muminus_PZ/F"); @@ -130,6 +178,7 @@ int mapmc_b02hphmmumu() output_tree->Branch("muminus_ID", &muminus_ID_out, "muminus_ID/I"); output_tree->Branch("muplus_BPVIPCHI2", &muplus_BPVIPCHI2_out, "muplus_BPVIPCHI2/F"); + output_tree->Branch("muplus_PT", &muplus_PT_out, "muplus_PT/F"); output_tree->Branch("muplus_PX", &muplus_PX_out, "muplus_PX/F"); output_tree->Branch("muplus_PY", &muplus_PY_out, "muplus_PY/F"); output_tree->Branch("muplus_PZ", &muplus_PZ_out, "muplus_PZ/F"); @@ -138,6 +187,7 @@ int mapmc_b02hphmmumu() output_tree->Branch("muplus_ID", &muplus_ID_out, "muplus_ID/I"); unsigned int sim_entries = sim_chain->GetEntries(); + unsigned int selected_entries = 0; for (unsigned int i = 0; i < sim_entries; i++) { @@ -147,12 +197,23 @@ int mapmc_b02hphmmumu() B0_PT_out = B0_PT_in; B0_BPVFDCHI2_out = B0_BPVFDCHI2_in; B0_BPVDIRA_out = B0_BPVDIRA_in; + B0_CHI2_out = B0_CHI2_in; + B0_CHI2DOF_out = B0_CHI2DOF_in; Jpsi_BPVIPCHI2_out = Jpsi_BPVIPCHI2_in; Jpsi_PT_out = Jpsi_PT_in; + Jpsi_BPVDIRA_out = Jpsi_BPVDIRA_in; + Jpsi_CHI2_out = Jpsi_CHI2_in; + Jpsi_CHI2DOF_out = Jpsi_CHI2DOF_in; + Jpsi_BPVFDCHI2_out = Jpsi_BPVFDCHI2_in; + + Kst0_BPVIPCHI2_out = Res_BPVIPCHI2_in; + Kst0_BPVFDCHI2_out = Res_BPVFDCHI2_in; + Kst0_PT_out = Res_PT_in; muminus_BPVIPCHI2_out = L1_BPVIPCHI2_in; muminus_ID_out = L1_TRUEID; + muminus_PT_out = L1_PT_in; auto muminus_lv = l14v_sim->LorentzVector(); muminus_PX_out = muminus_lv.Px(); muminus_PY_out = muminus_lv.Py(); @@ -162,6 +223,7 @@ int mapmc_b02hphmmumu() muplus_BPVIPCHI2_out = L2_BPVIPCHI2_in; muplus_ID_out = L2_TRUEID; + muplus_PT_out = L2_PT_in; auto muplus_lv = l24v_sim->LorentzVector(); muplus_PX_out = muplus_lv.Px(); muplus_PY_out = muplus_lv.Py(); @@ -174,6 +236,7 @@ int mapmc_b02hphmmumu() Kplus_BPVIPCHI2_out = Hp_BPVIPCHI2_in; Kplus_PT_out = Hp_PT_in; Kplus_PROBNN_K_out = Hp_PROBNN_K_in; + Kplus_PID_K_out = Hp_PID_K_in; Kplus_ID_out = Hp_TRUEID; piminus_BPVIPCHI2_out = Hm_BPVIPCHI2_in; piminus_PT_out = Hm_PT_in; @@ -208,6 +271,7 @@ int mapmc_b02hphmmumu() Kplus_BPVIPCHI2_out = Hm_BPVIPCHI2_in; Kplus_PT_out = Hm_PROBNN_K_in; Kplus_PROBNN_K_out = Hm_PROBNN_K_in; + Kplus_PID_K_out = Hm_PID_K_in; Kplus_ID_out = Hm_TRUEID; piminus_BPVIPCHI2_out = Hp_BPVIPCHI2_in; piminus_PT_out = Hp_PT_in; @@ -239,6 +303,7 @@ int mapmc_b02hphmmumu() } output_tree->Fill(); + selected_entries++; } PrintProgress(TString::Format("%s Mapping", sim_tree_name), sim_entries, 10000, i); @@ -252,5 +317,7 @@ int mapmc_b02hphmmumu() output_tree->Write(); output_file->Close(); + std::cout << "Triggered Entries: " << sim_entries << " Selected Entries: " << selected_entries << std::endl; + return 0; } \ No newline at end of file diff --git a/mapmc_bu2hpmumu.cpp b/mapmc_bu2hpmumu.cpp index ed488d0..065e4c7 100644 --- a/mapmc_bu2hpmumu.cpp +++ b/mapmc_bu2hpmumu.cpp @@ -57,46 +57,80 @@ int mapmc_bu2hpmumu() sim_chain->SetBranchAddress("Hp_TRUEID", &Hp_TRUEID); sim_chain->SetBranchAddress("B_BKGCAT", &B_BKGCAT); - Float_t B_PT_in, B_BPVFDCHI2_in, B_BPVDIRA_in, Jpsi_BPVIPCHI2_in, Jpsi_PT_in, - Hp_BPVIPCHI2_in, Hp_PT_in, L1_BPVIPCHI2_in, L2_BPVIPCHI2_in; - Double_t Hp_PROBNN_K_in; + Float_t B_PT_in, B_BPVFDCHI2_in, B_BPVDIRA_in, + Jpsi_BPVIPCHI2_in, Jpsi_PT_in, Jpsi_BPVDIRA_in, + Hp_BPVIPCHI2_in, Hp_PT_in, + L1_BPVIPCHI2_in, L1_PT_in, + L2_BPVIPCHI2_in, L2_PT_in; + Double_t B_CHI2_in, B_CHI2DOF_in, + Jpsi_CHI2_in, Jpsi_CHI2DOF_in, + Hp_PROBNN_K_in, Hp_PID_K_in, + L1_PID_MU_in, + L2_PID_MU_in; sim_chain->SetBranchAddress("B_PT", &B_PT_in); sim_chain->SetBranchAddress("B_BPVFDCHI2", &B_BPVFDCHI2_in); sim_chain->SetBranchAddress("B_BPVDIRA", &B_BPVDIRA_in); + sim_chain->SetBranchAddress("B_CHI2", &B_CHI2_in); + sim_chain->SetBranchAddress("B_CHI2VXNDOF", &B_CHI2DOF_in); + sim_chain->SetBranchAddress("Jpsi_BPVIPCHI2", &Jpsi_BPVIPCHI2_in); + sim_chain->SetBranchAddress("Jpsi_BPVDIRA", &Jpsi_BPVDIRA_in); + sim_chain->SetBranchAddress("Jpsi_CHI2", &Jpsi_CHI2_in); + sim_chain->SetBranchAddress("Jpsi_CHI2DOF", &Jpsi_CHI2DOF_in); sim_chain->SetBranchAddress("Jpsi_PT", &Jpsi_PT_in); + sim_chain->SetBranchAddress("Hp_BPVIPCHI2", &Hp_BPVIPCHI2_in); sim_chain->SetBranchAddress("Hp_PT", &Hp_PT_in); + sim_chain->SetBranchAddress("Hp_PROBNN_K", &Hp_PROBNN_K_in); + sim_chain->SetBranchAddress("Hp_PID_K", &Hp_PID_K_in); + sim_chain->SetBranchAddress("L1_BPVIPCHI2", &L1_BPVIPCHI2_in); + sim_chain->SetBranchAddress("L1_PT", &L1_PT_in); + sim_chain->SetBranchAddress("L1_PID_MU", &L1_PID_MU_in); + sim_chain->SetBranchAddress("L2_BPVIPCHI2", &L2_BPVIPCHI2_in); - sim_chain->SetBranchAddress("Hp_PROBNN_K", &Hp_PROBNN_K_in); + sim_chain->SetBranchAddress("L2_PT", &L2_PT_in); + sim_chain->SetBranchAddress("L2_PID_MU", &L2_PID_MU_in); TTree *output_tree = new TTree("DecayTree", "DecayTree"); - Float_t B_PT_out, B_BPVFDCHI2_out, B_BPVDIRA_out, Jpsi_BPVIPCHI2_out, Jpsi_PT_out, - Kplus_BPVIPCHI2_out, Kplus_PT_out, muplus_BPVIPCHI2_out, muminus_BPVIPCHI2_out, - B_PX_out, B_PY_out, B_PZ_out, B_ENERGY_out, - muplus_PX_out, muplus_PY_out, muplus_PZ_out, muplus_ENERGY_out, - muminus_PX_out, muminus_PY_out, muminus_PZ_out, muminus_ENERGY_out, - Kplus_PX_out, Kplus_PY_out, Kplus_PZ_out, Kplus_ENERGY_out; - Double_t Kplus_PROBNN_K_out, B_M_out, muplus_M_out, muminus_M_out, Kplus_M_out; + Float_t B_PT_out, B_BPVFDCHI2_out, B_BPVDIRA_out, + Jpsi_BPVIPCHI2_out, Jpsi_PT_out, Jpsi_BPVDIRA_out, + Kplus_BPVIPCHI2_out, Kplus_PT_out, + muplus_BPVIPCHI2_out, muplus_PT_out, + muminus_BPVIPCHI2_out, muminus_PT_out, + B_PX_out, B_PY_out, B_PZ_out, B_ENERGY_out, + muplus_PX_out, muplus_PY_out, muplus_PZ_out, muplus_ENERGY_out, + muminus_PX_out, muminus_PY_out, muminus_PZ_out, muminus_ENERGY_out, + Kplus_PX_out, Kplus_PY_out, Kplus_PZ_out, Kplus_ENERGY_out; + Double_t B_M_out, B_CHI2_out, B_CHI2DOF_out, + Jpsi_CHI2_out, Jpsi_CHI2DOF_out, + Kplus_PROBNN_K_out, Kplus_M_out, Kplus_PID_K_out, + muplus_M_out, + muminus_M_out; Int_t muplus_ID_out, muminus_ID_out, Kplus_ID_out; - output_tree->Branch("B_PT", &B_PT_out, "B_PT/D"); + output_tree->Branch("B_PT", &B_PT_out, "B_PT/F"); output_tree->Branch("B_BPVFDCHI2", &B_BPVFDCHI2_out, "B_BPVFDCHI2/F"); output_tree->Branch("B_BPVDIRA", &B_BPVDIRA_out, "B_BPVDIRA/F"); - output_tree->Branch("B_PX", &B_PX_out, "0_PX/F"); - output_tree->Branch("B_PY", &B_PY_out, "0_PY/F"); - output_tree->Branch("B_PZ", &B_PZ_out, "0_PZ/F"); + output_tree->Branch("B_CHI2", &B_CHI2_out, "B_CHI2/D"); + output_tree->Branch("B_CHI2DOF", &B_CHI2DOF_out, "B_CHI2DOF/D"); + output_tree->Branch("B_PX", &B_PX_out, "B_PX/F"); + output_tree->Branch("B_PY", &B_PY_out, "B_PY/F"); + output_tree->Branch("B_PZ", &B_PZ_out, "B_PZ/F"); output_tree->Branch("B_ENERGY", &B_ENERGY_out, "B_ENERGY/F"); output_tree->Branch("B_M", &B_M_out, "B_M/D"); output_tree->Branch("Jpsi_BPVIPCHI2", &Jpsi_BPVIPCHI2_out, "Jpsi_BPVIPCHI2/F"); output_tree->Branch("Jpsi_PT", &Jpsi_PT_out, "Jpsi_PT/F"); + output_tree->Branch("Jpsi_BPVDIRA", &Jpsi_BPVDIRA_out, "Jpsi_BPVDIRA/F"); + output_tree->Branch("Jpsi_CHI2", &Jpsi_CHI2_out, "Jpsi_CHI2/D"); + output_tree->Branch("Jpsi_CHI2DOF", &Jpsi_CHI2DOF_out, "Jpsi_CHI2DOF/D"); output_tree->Branch("Kplus_BPVIPCHI2", &Kplus_BPVIPCHI2_out, "Kplus_BPVIPCHI2/F"); output_tree->Branch("Kplus_PROBNN_K", &Kplus_PROBNN_K_out, "Kplus_PROBNN_K/D"); + output_tree->Branch("Kplus_PID_K", &Kplus_PID_K_out, "Kplus_PID_K/D"); output_tree->Branch("Kplus_PT", &Kplus_PT_out, "Kplus_PT/F"); output_tree->Branch("Kplus_PX", &Kplus_PX_out, "Kplus_PX/F"); output_tree->Branch("Kplus_PY", &Kplus_PY_out, "Kplus_PY/F"); @@ -106,6 +140,7 @@ int mapmc_bu2hpmumu() output_tree->Branch("Kplus_ID", &Kplus_ID_out, "Kplus_ID/I"); output_tree->Branch("muminus_BPVIPCHI2", &muminus_BPVIPCHI2_out, "muminus_BPVIPCHI2/F"); + output_tree->Branch("muminus_PT", &muminus_PT_out, "muminus_PT/F"); output_tree->Branch("muminus_PX", &muminus_PX_out, "muminus_PX/F"); output_tree->Branch("muminus_PY", &muminus_PY_out, "muminus_PY/F"); output_tree->Branch("muminus_PZ", &muminus_PZ_out, "muminus_PZ/F"); @@ -114,6 +149,7 @@ int mapmc_bu2hpmumu() output_tree->Branch("muminus_ID", &muminus_ID_out, "muminus_ID/I"); output_tree->Branch("muplus_BPVIPCHI2", &muplus_BPVIPCHI2_out, "muplus_BPVIPCHI2/F"); + output_tree->Branch("muplus_PT", &muplus_PT_out, "muplus_PT/F"); output_tree->Branch("muplus_PX", &muplus_PX_out, "muplus_PX/F"); output_tree->Branch("muplus_PY", &muplus_PY_out, "muplus_PY/F"); output_tree->Branch("muplus_PZ", &muplus_PZ_out, "muplus_PZ/F"); @@ -122,6 +158,7 @@ int mapmc_bu2hpmumu() output_tree->Branch("muplus_ID", &muplus_ID_out, "muplus_ID/I"); unsigned int sim_entries = sim_chain->GetEntries(); + unsigned int selected_entries = 0; for (unsigned int i = 0; i < sim_entries; i++) { @@ -131,12 +168,18 @@ int mapmc_bu2hpmumu() B_PT_out = B_PT_in; B_BPVFDCHI2_out = B_BPVFDCHI2_in; B_BPVDIRA_out = B_BPVDIRA_in; + B_CHI2_out = B_CHI2_in; + B_CHI2DOF_out = B_CHI2DOF_in; Jpsi_BPVIPCHI2_out = Jpsi_BPVIPCHI2_in; Jpsi_PT_out = Jpsi_PT_in; + Jpsi_BPVDIRA_out = Jpsi_BPVDIRA_in; + Jpsi_CHI2_out = Jpsi_CHI2_in; + Jpsi_CHI2DOF_out = Jpsi_CHI2DOF_in; muminus_BPVIPCHI2_out = L1_BPVIPCHI2_in; muminus_ID_out = L1_TRUEID; + muminus_PT_out = L1_PT_in; auto muminus_lv = l14v_sim->LorentzVector(); muminus_PX_out = muminus_lv.Px(); muminus_PY_out = muminus_lv.Py(); @@ -146,6 +189,7 @@ int mapmc_bu2hpmumu() muplus_BPVIPCHI2_out = L2_BPVIPCHI2_in; muplus_ID_out = L2_TRUEID; + muplus_PT_out = L2_PT_in; auto muplus_lv = l24v_sim->LorentzVector(); muplus_PX_out = muplus_lv.Px(); muplus_PY_out = muplus_lv.Py(); @@ -156,6 +200,7 @@ int mapmc_bu2hpmumu() Kplus_BPVIPCHI2_out = Hp_BPVIPCHI2_in; Kplus_PT_out = Hp_PT_in; Kplus_PROBNN_K_out = Hp_PROBNN_K_in; + Kplus_PID_K_out = Hp_PID_K_in; Kplus_ID_out = Hp_TRUEID; auto Kplus_lv = hp4v_sim->LorentzVector(K_MASS); @@ -175,6 +220,8 @@ int mapmc_bu2hpmumu() B_M_out = B_lv.M(); output_tree->Fill(); + + selected_entries++; } PrintProgress(TString::Format("%s Mapping", sim_tree_name), sim_entries, 10000, i); @@ -188,5 +235,7 @@ int mapmc_bu2hpmumu() output_tree->Write(); output_file->Close(); + std::cout << "Triggered Entries: " << sim_entries << " Selected Entries: " << selected_entries << std::endl; + return 0; } \ No newline at end of file diff --git a/new_analysis_b02hphmmumu.cpp b/new_analysis_b02hphmmumu.cpp index 6bd3465..efd7798 100644 --- a/new_analysis_b02hphmmumu.cpp +++ b/new_analysis_b02hphmmumu.cpp @@ -54,12 +54,17 @@ int new_analysis_b02hphmmumu() 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_Sprucing23r1_90000000_RD.root"); + // data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Comm23_B6800GeV-VeCl-MgDo-Excl-UT_RealD_SprPass23r1_94000000_RD_45.root"); FourVect *l14v_data = FourVect::Init(data_chain, "muminus"); FourVect *l24v_data = FourVect::Init(data_chain, "muplus"); FourVect *hp4v_data = FourVect::Init(data_chain, "Kplus"); FourVect *hm4v_data = FourVect::Init(data_chain, "piminus"); + Double_t piminus_PID_K; + + data_chain->SetBranchAddress("piminus_PID_K", &piminus_PID_K); + TChain *sim_chain = new TChain(TString::Format("%s/DecayTree", sim_tree_name)); sim_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/B0ToHpHmMuMu_mapped_mc.root"); @@ -73,19 +78,19 @@ int new_analysis_b02hphmmumu() 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)); + TTree *tree_B_Mass_jpsi = new TTree("tree_B_Mass_jpsi", TString::Format("B^{0} Mass, J/#psi Mode (%s)", analysis_name)); + TTree *tree_B_Mass_psi2s = new TTree("tree_B_Mass_psi2s", TString::Format("B^{0} Mass, #psi(2S) Mode (%s)", analysis_name)); + TTree *tree_B_Mass_sim = new TTree("tree_B_Mass_sim", TString::Format("B^{0} Mass, Simualted (%s)", 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_unf = new TH1D("h1_B_Mass_unf", TString::Format("B^{0} Mass (%s), Unfiltered", analysis_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^{0} Mass (%s), BDT Filter", 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.); + TH2D *h2_Hlt1_flags_B_Mass = new TH2D("h2_Hlt1_flags_B_Mass", "Hlt1 Decision vs B^{0} 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^{0} Mass", 50, 5100, 5400, 13, 1., 14.); TH1D *h1_bdt_probs = new TH1D("h1_bdt_probs", "BDT Probabilities", 100, -1, 1); h1_B_Mass_unf->GetXaxis()->SetTitle(end_state_mass_literal); @@ -97,19 +102,25 @@ int new_analysis_b02hphmmumu() auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name); std::map exclusive_hits{}; + TV* kplus_pid_k_var = TV::Double("Kplus_PID_K", "Kplus_PID_K"); + std::vector vars{ TV::Float("B0_PT", "B0_PT"), TV::Float("B0_BPVFDCHI2", "B0_BPVFDCHI2"), TV::Float("B0_BPVDIRA", "B0_BPVDIRA"), TV::Float("Jpsi_BPVIPCHI2", "Jpsi_BPVIPCHI2"), TV::Float("Jpsi_PT", "Jpsi_PT"), + TV::Float("Kst0_BPVIPCHI2", "Kst0_BPVIPCHI2"), + TV::Float("Kst0_PT", "Kst0_PT"), TV::Float("Kplus_BPVIPCHI2", "Kplus_BPVIPCHI2"), TV::Float("Kplus_PT", "Kplus_PT"), TV::Float("piminus_BPVIPCHI2", "piminus_BPVIPCHI2"), TV::Float("piminus_PT", "piminus_PT"), - TV::Double("Kplus_PROBNN_K", "Kplus_PROBNN_K"), + kplus_pid_k_var, TV::Float("muminus_BPVIPCHI2", "muminus_BPVIPCHI2"), + // TV::Float("muminus_PT", "muminus_PT"), TV::Float("muplus_BPVIPCHI2", "muplus_BPVIPCHI2"), + // TV::Float("muplus_PT", "muplus_PT"), }; TTree *sig_tree = new TTree("TreeS", "tree containing signal data"); @@ -188,7 +199,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; + const double mva_cut_value = 0.0; for (unsigned int i = 0; i < data_entries; i++) { @@ -225,7 +236,7 @@ int new_analysis_b02hphmmumu() h1_B_Mass_unf->Fill(reconstructed_B_Mass); - if (mva_response > mva_cut_value) + if (mva_response > mva_cut_value && kplus_pid_k_var->GetDataDouble() > -3 && (kplus_pid_k_var->GetDataDouble() - piminus_PID_K) > 0) { h1_B_Mass_bdtf->Fill(reconstructed_B_Mass); if (TMath::Abs(reconstructed_Kstar.M() - KSTAR_MASS) < 100.) @@ -258,10 +269,13 @@ int new_analysis_b02hphmmumu() DrawInDefaultCanvas(h1_B_Mass_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); + DrawInDefaultCanvasStacked({h1_B_Mass_unf, h1_B_Mass_bdtf}, {kRed, kBlue}, {0, 3003}, analysis_name, TString::Format("B^{0} Mass (%s) Before/After BDT Selection", analysis_name)); 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); + + roofit_hist_sim.shape_parameters.sigma_lr = roofit_hist_jpsi_fitsum.shape_parameters.sigma_lr; + 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, true); DrawInDefaultCanvas(roofit_hist_jpsi_fitsum, analysis_name); @@ -274,18 +288,61 @@ int new_analysis_b02hphmmumu() auto signal_ratio = DivWithErr(roofit_hist_psi2s_fitsum.signal_yield.first, roofit_hist_psi2s_fitsum.signal_yield.second, roofit_hist_jpsi_fitsum.signal_yield.first, roofit_hist_jpsi_fitsum.signal_yield.second); - std::time_t t = std::time(nullptr); + std::time_t t = std::time(nullptr); std::tm tm = *std::localtime(&t); ofstream res_file; res_file.open(TString::Format("%s_results.txt", analysis_name).Data(), ios::out | ios::trunc); - res_file << "#### " << std::put_time(&tm, "%c") << " ####"<< std::endl; - res_file << "J/Psi Mode: " << roofit_hist_jpsi_fitsum.signal_yield.first << " +- " << roofit_hist_jpsi_fitsum.signal_yield.second << std::endl; - res_file << "Psi(2S) Mode: " << roofit_hist_psi2s_fitsum.signal_yield.first << " +- " << roofit_hist_psi2s_fitsum.signal_yield.second << std::endl; - res_file << "Mode Yield Ratio: " << signal_ratio.first << " +- " << signal_ratio.second << std::endl; - res_file << "Rel Br Frac MuMu: " << (BRF_JPSI_MUMU_VAL / BRF_PSI2S_MUMU_VAL) << std::endl; - res_file << "Rel Br Frac: " << signal_ratio.first * (BRF_JPSI_MUMU_VAL / BRF_PSI2S_MUMU_VAL) << std::endl; + auto jpsi_sigobkg = DivWithErr(roofit_hist_jpsi_fitsum.signal_yield.first, roofit_hist_jpsi_fitsum.signal_yield.second, roofit_hist_jpsi_fitsum.background_yield.first, roofit_hist_jpsi_fitsum.background_yield.second); + auto psi2s_sigobkg = DivWithErr(roofit_hist_psi2s_fitsum.signal_yield.first, roofit_hist_psi2s_fitsum.signal_yield.second, roofit_hist_psi2s_fitsum.background_yield.first, roofit_hist_psi2s_fitsum.background_yield.second); + + auto mumu_br_frac = DivWithErr(BRF_JPSI_MUMU_VAL, BRF_JPSI_MUMU_ERR, BRF_PSI2S_MUMU_VAL, BRF_PSI2S_MUMU_ERR); + + res_file << "#### " << analysis_name << " @ " << std::put_time(&tm, "%c") << " ####" << std::endl; + res_file << "J/Psi Mode: " << ErrToStr(roofit_hist_jpsi_fitsum.signal_yield, 0) << " / " << ErrToStr(roofit_hist_jpsi_fitsum.background_yield, 0) << std::endl; + res_file << " Sig/Bkg: " << ErrToStr(jpsi_sigobkg, 2) << std::endl; + res_file << "Psi(2S) Mode: " << ErrToStr(roofit_hist_psi2s_fitsum.signal_yield, 0) << " / " << ErrToStr(roofit_hist_psi2s_fitsum.background_yield, 0) << std::endl; + res_file << " Sig/Bkg: " << ErrToStr(psi2s_sigobkg, 2) << std::endl; + res_file << "Mode Yield Ratio: " << ErrToStr(signal_ratio, 3) << std::endl; + res_file << "Rel Br Frac MuMu: " << ErrToStr(mumu_br_frac, 3) << std::endl; + + auto br_frac = MultWithErr(signal_ratio.first, signal_ratio.second, mumu_br_frac.first, mumu_br_frac.second); + res_file << "Rel Br Frac: " << ErrToStr(br_frac, 3) << std::endl << std::endl; + + res_file << std::endl << "Fitted Parameters: Simulation" << std::endl << std::endl; + + for (const auto &par : roofit_hist_sim.fitted_params) + { + res_file << par.ToString(true).c_str() << std::endl; + } + + res_file << std::endl << "Fitted Parameters: J/PSI" << std::endl << std::endl; + + for (const auto &par : roofit_hist_jpsi_fitsum.fitted_params) + { + res_file << par.ToString(true).c_str() << std::endl; + } + + res_file << std::endl << "Fitted Parameters: PSI(2S)" << std::endl << std::endl; + + for (const auto &par : roofit_hist_psi2s_fitsum.fitted_params) + { + res_file << par.ToString(true).c_str() << std::endl; + } + + auto print_table = [&res_file](std::string name, std::pair sig, std::pair bkg) + { + res_file << std::endl; + res_file << "# " << name << std::endl; + res_file << "\\begin{tabular}{c|c}" << std::endl; + res_file << "$N_{Sig}$ & $N_{Bkg}$\\\\\\hline" << std::endl; + res_file << TString::Format("$%.0f \\pm %.0f$ & $%.0f \\pm %.0f$", sig.first, sig.second, bkg.first, bkg.second).Data() << std::endl; + res_file << "\\end{tabular}" << std::endl; + }; + + print_table("J/psi", roofit_hist_jpsi_fitsum.signal_yield, roofit_hist_jpsi_fitsum.background_yield); + print_table("psi(2S)", roofit_hist_psi2s_fitsum.signal_yield, roofit_hist_psi2s_fitsum.background_yield); res_file.close(); diff --git a/new_analysis_b02kppimmumu.cpp b/new_analysis_b02kppimmumu.cpp index 8033106..b95d68b 100644 --- a/new_analysis_b02kppimmumu.cpp +++ b/new_analysis_b02kppimmumu.cpp @@ -64,7 +64,7 @@ int new_analysis_b02kppimmumu() 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"); + sim_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/rd_btoxll_simulation_turbo_v0r0p6862057_B0ToKpPimMuMu_11144002_magdown.root"); FourVect *l14v_sim = FourVect::Init(sim_chain, "L1"); FourVect *l24v_sim = FourVect::Init(sim_chain, "L2"); @@ -84,19 +84,19 @@ int new_analysis_b02kppimmumu() 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)); + TTree *tree_B_Mass_jpsi = new TTree("tree_B_Mass_jpsi", TString::Format("B^{0} Mass, J/#psi Mode (%s)", analysis_name)); + TTree *tree_B_Mass_psi2s = new TTree("tree_B_Mass_psi2s", TString::Format("B^{0} Mass, #psi(2S) Mode (%s)", analysis_name)); + TTree *tree_B_Mass_sim = new TTree("tree_B_Mass_sim", TString::Format("B^{0} 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); + TH1D *h1_B_Mass_unf = new TH1D("h1_B_Mass_unf", TString::Format("B^{0} Mass (%s), Unfiltered", analysis_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^{0} 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.); + TH2D *h2_Hlt1_flags_B_Mass = new TH2D("h2_Hlt1_flags_B_Mass", "Hlt1 Decision vs B^{0} 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^{0} 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); @@ -174,7 +174,10 @@ int new_analysis_b02kppimmumu() 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); + + roofit_hist_sim.shape_parameters.sigma_lr = roofit_hist_jpsi_fitsum.shape_parameters.sigma_lr; + + 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, true); DrawInDefaultCanvas(roofit_hist_jpsi_fitsum, analysis_name); DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name); @@ -190,12 +193,56 @@ int new_analysis_b02kppimmumu() ofstream res_file; res_file.open(TString::Format("%s_results.txt", analysis_name).Data(), ios::out | ios::trunc); - res_file << "#### " << std::put_time(&tm, "%c") << " ####"<< std::endl; - res_file << "J/Psi Mode: " << roofit_hist_jpsi_fitsum.signal_yield.first << " +- " << roofit_hist_jpsi_fitsum.signal_yield.second << std::endl; - res_file << "Psi(2S) Mode: " << roofit_hist_psi2s_fitsum.signal_yield.first << " +- " << roofit_hist_psi2s_fitsum.signal_yield.second << std::endl; - res_file << "Mode Yield Ratio: " << signal_ratio.first << " +- " << signal_ratio.second << std::endl; - res_file << "Rel Br Frac MuMu: " << (BRF_JPSI_MUMU_VAL / BRF_PSI2S_MUMU_VAL) << std::endl; - res_file << "Rel Br Frac: " << signal_ratio.first * (BRF_JPSI_MUMU_VAL / BRF_PSI2S_MUMU_VAL) << std::endl; + auto jpsi_sigobkg = DivWithErr(roofit_hist_jpsi_fitsum.signal_yield.first, roofit_hist_jpsi_fitsum.signal_yield.second, roofit_hist_jpsi_fitsum.background_yield.first, roofit_hist_jpsi_fitsum.background_yield.second); + auto psi2s_sigobkg = DivWithErr(roofit_hist_psi2s_fitsum.signal_yield.first, roofit_hist_psi2s_fitsum.signal_yield.second, roofit_hist_psi2s_fitsum.background_yield.first, roofit_hist_psi2s_fitsum.background_yield.second); + + auto mumu_br_frac = DivWithErr(BRF_JPSI_MUMU_VAL, BRF_JPSI_MUMU_ERR, BRF_PSI2S_MUMU_VAL, BRF_PSI2S_MUMU_ERR); + + res_file << "#### " << analysis_name << " @ " << std::put_time(&tm, "%c") << " ####" << std::endl; + res_file << "J/Psi Mode: " << ErrToStr(roofit_hist_jpsi_fitsum.signal_yield, 0) << " / " << ErrToStr(roofit_hist_jpsi_fitsum.background_yield, 0) << std::endl; + res_file << " Sig/Bkg: " << ErrToStr(jpsi_sigobkg, 2) << std::endl; + res_file << "Psi(2S) Mode: " << ErrToStr(roofit_hist_psi2s_fitsum.signal_yield, 0) << " / " << ErrToStr(roofit_hist_psi2s_fitsum.background_yield, 0) << std::endl; + res_file << " Sig/Bkg: " << ErrToStr(psi2s_sigobkg, 2) << std::endl; + res_file << "Mode Yield Ratio: " << ErrToStr(signal_ratio, 3) << std::endl; + res_file << "Rel Br Frac MuMu: " << ErrToStr(mumu_br_frac, 3) << std::endl; + + auto br_frac = MultWithErr(signal_ratio.first, signal_ratio.second, mumu_br_frac.first, mumu_br_frac.second); + res_file << "Rel Br Frac: " << ErrToStr(br_frac, 3) << std::endl << std::endl; + + res_file << "Params from Sim:" << std::endl << roofit_hist_sim.shape_parameters.ToString() << std::endl; + + for (const auto &par : roofit_hist_sim.fitted_params) + { + res_file << par.ToString(true).c_str() << std::endl; + } + + res_file << std::endl << "Fitted Parameters: J/PSI" << std::endl << std::endl; + + for (const auto &par : roofit_hist_jpsi_fitsum.fitted_params) + { + res_file << par.ToString(true).c_str() << std::endl; + } + + res_file << std::endl << "Fitted Parameters: PSI(2S)" << std::endl << std::endl; + + for (const auto &par : roofit_hist_psi2s_fitsum.fitted_params) + { + res_file << par.ToString(true).c_str() << std::endl; + } + + + auto print_table = [&res_file](std::string name, std::pair sig, std::pair bkg) + { + res_file << std::endl; + res_file << "# " << name << std::endl; + res_file << "\\begin{tabular}{c|c}" << std::endl; + res_file << "$N_{Sig}$ & $N_{Bkg}$\\\\\\hline" << std::endl; + res_file << TString::Format("$%.0f \\pm %.0f$ & $%.0f \\pm %.0f$", sig.first, sig.second, bkg.first, bkg.second).Data() << std::endl; + res_file << "\\end{tabular}" << std::endl; + }; + + print_table("J/psi", roofit_hist_jpsi_fitsum.signal_yield, roofit_hist_jpsi_fitsum.background_yield); + print_table("psi(2S)", roofit_hist_psi2s_fitsum.signal_yield, roofit_hist_psi2s_fitsum.background_yield); res_file.close(); diff --git a/new_analysis_bu2hpmumu.cpp b/new_analysis_bu2hpmumu.cpp index b7cba27..48b29cd 100644 --- a/new_analysis_bu2hpmumu.cpp +++ b/new_analysis_bu2hpmumu.cpp @@ -54,11 +54,16 @@ int new_analysis_bu2hpmumu() 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_Sprucing23r1_90000000_RD.root"); + // data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Comm23_B6800GeV-VeCl-MgDo-Excl-UT_RealD_SprPass23r1_94000000_RD_45.root"); FourVect *l14v_data = FourVect::Init(data_chain, "muminus"); FourVect *l24v_data = FourVect::Init(data_chain, "muplus"); FourVect *hp4v_data = FourVect::Init(data_chain, "Kplus"); + // Double_t Kplus_PID_K; + + // data_chain->SetBranchAddress("Kplus_PID_K", &Kplus_PID_K); + TChain *sim_chain = new TChain(TString::Format("%s/DecayTree", sim_tree_name)); sim_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/BuToHpMuMu_mapped_mc.root"); @@ -71,19 +76,19 @@ int new_analysis_bu2hpmumu() 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)); + TTree *tree_B_Mass_jpsi = new TTree("tree_B_Mass_jpsi", TString::Format("B^{#pm} Mass, J/#psi Mode (%s)", analysis_name)); + TTree *tree_B_Mass_psi2s = new TTree("tree_B_Mass_psi2s", TString::Format("B^{#pm} Mass, #psi(2S) Mode (%s)", analysis_name)); + TTree *tree_B_Mass_sim = new TTree("tree_B_Mass_sim", TString::Format("B^{#pm} Mass, Simualted (%s)", 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_unf = new TH1D("h1_B_Mass_unf", TString::Format("B^{#pm} Mass (%s), Unfiltered", analysis_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^{#pm} Mass (%s), BDT Filter", 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.); + TH2D *h2_Hlt1_flags_B_Mass = new TH2D("h2_Hlt1_flags_B_Mass", "Hlt1 Decision vs B^{#pm} 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^{#pm} Mass", 50, 5100, 5400, 13, 1., 14.); TH1D *h1_bdt_probs = new TH1D("h1_bdt_probs", "BDT Probabilities", 100, -2, 2); h1_B_Mass_unf->GetXaxis()->SetTitle(end_state_mass_literal); @@ -95,6 +100,8 @@ int new_analysis_bu2hpmumu() auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name); std::map exclusive_hits{}; + TV* kplus_pid_k_var = TV::Double("Kplus_PID_K", "Kplus_PID_K"); + std::vector vars{ TV::Float("B_PT", "B_PT"), TV::Float("B_BPVFDCHI2", "B_BPVFDCHI2"), @@ -103,9 +110,11 @@ int new_analysis_bu2hpmumu() TV::Float("Jpsi_PT", "Jpsi_PT"), TV::Float("Kplus_BPVIPCHI2", "Kplus_BPVIPCHI2"), TV::Float("Kplus_PT", "Kplus_PT"), - TV::Double("Kplus_PROBNN_K", "Kplus_PROBNN_K"), + kplus_pid_k_var, TV::Float("muminus_BPVIPCHI2", "muminus_BPVIPCHI2"), + // TV::Float("muminus_PT", "muminus_PT"), TV::Float("muplus_BPVIPCHI2", "muplus_BPVIPCHI2"), + // TV::Float("muplus_PT", "muplus_PT"), }; TTree *sig_tree = new TTree("TreeS", "tree containing signal data"); @@ -182,7 +191,7 @@ int new_analysis_bu2hpmumu() Float_t *train_vars = new Float_t[vars.size()]; auto reader = SetupReader(vars, train_vars, analysis_name); - const double mva_cut_value = -0.05; + const double mva_cut_value = 0.0; for (unsigned int i = 0; i < data_entries; i++) { @@ -218,7 +227,7 @@ int new_analysis_bu2hpmumu() h1_B_Mass_unf->Fill(reconstructed_B_Mass); - if (mva_response > mva_cut_value) + if (mva_response > mva_cut_value && kplus_pid_k_var->GetDataDouble() > -3) { h1_B_Mass_bdtf->Fill(reconstructed_B_Mass); if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) @@ -248,10 +257,13 @@ int new_analysis_bu2hpmumu() DrawInDefaultCanvas(h1_B_Mass_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); + DrawInDefaultCanvasStacked({h1_B_Mass_unf, h1_B_Mass_bdtf}, {kRed, kBlue}, {0, 3003}, analysis_name, TString::Format("B^{#pm} Mass (%s) Before/After BDT Selection", analysis_name)); 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); + + roofit_hist_sim.shape_parameters.sigma_lr = roofit_hist_jpsi_fitsum.shape_parameters.sigma_lr; + 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, true); DrawInDefaultCanvas(roofit_hist_jpsi_fitsum, analysis_name); @@ -270,12 +282,56 @@ int new_analysis_bu2hpmumu() ofstream res_file; res_file.open(TString::Format("%s_results.txt", analysis_name).Data(), ios::out | ios::trunc); - res_file << "#### " << std::put_time(&tm, "%c") << " ####"<< std::endl; - res_file << "J/Psi Mode: " << roofit_hist_jpsi_fitsum.signal_yield.first << " +- " << roofit_hist_jpsi_fitsum.signal_yield.second << std::endl; - res_file << "Psi(2S) Mode: " << roofit_hist_psi2s_fitsum.signal_yield.first << " +- " << roofit_hist_psi2s_fitsum.signal_yield.second << std::endl; - res_file << "Mode Yield Ratio: " << signal_ratio.first << " +- " << signal_ratio.second << std::endl; - res_file << "Rel Br Frac MuMu: " << (BRF_JPSI_MUMU_VAL / BRF_PSI2S_MUMU_VAL) << std::endl; - res_file << "Rel Br Frac: " << signal_ratio.first * (BRF_JPSI_MUMU_VAL / BRF_PSI2S_MUMU_VAL) << std::endl; + auto jpsi_sigobkg = DivWithErr(roofit_hist_jpsi_fitsum.signal_yield.first, roofit_hist_jpsi_fitsum.signal_yield.second, roofit_hist_jpsi_fitsum.background_yield.first, roofit_hist_jpsi_fitsum.background_yield.second); + auto psi2s_sigobkg = DivWithErr(roofit_hist_psi2s_fitsum.signal_yield.first, roofit_hist_psi2s_fitsum.signal_yield.second, roofit_hist_psi2s_fitsum.background_yield.first, roofit_hist_psi2s_fitsum.background_yield.second); + + auto mumu_br_frac = DivWithErr(BRF_JPSI_MUMU_VAL, BRF_JPSI_MUMU_ERR, BRF_PSI2S_MUMU_VAL, BRF_PSI2S_MUMU_ERR); + + res_file << "#### " << analysis_name << " @ " << std::put_time(&tm, "%c") << " ####" << std::endl; + res_file << "J/Psi Mode: " << ErrToStr(roofit_hist_jpsi_fitsum.signal_yield, 0) << " / " << ErrToStr(roofit_hist_jpsi_fitsum.background_yield, 0) << std::endl; + res_file << " Sig/Bkg: " << ErrToStr(jpsi_sigobkg, 2) << std::endl; + res_file << "Psi(2S) Mode: " << ErrToStr(roofit_hist_psi2s_fitsum.signal_yield, 0) << " / " << ErrToStr(roofit_hist_psi2s_fitsum.background_yield, 0) << std::endl; + res_file << " Sig/Bkg: " << ErrToStr(psi2s_sigobkg, 2) << std::endl; + res_file << "Mode Yield Ratio: " << ErrToStr(signal_ratio, 3) << std::endl; + res_file << "Rel Br Frac MuMu: " << ErrToStr(mumu_br_frac, 3) << std::endl; + + auto br_frac = MultWithErr(signal_ratio.first, signal_ratio.second, mumu_br_frac.first, mumu_br_frac.second); + res_file << "Rel Br Frac: " << ErrToStr(br_frac, 3) << std::endl << std::endl; + + res_file << "Params from Sim:" << std::endl << roofit_hist_sim.shape_parameters.ToString() << std::endl; + + for (const auto &par : roofit_hist_sim.fitted_params) + { + res_file << par.ToString(true).c_str() << std::endl; + } + + res_file << std::endl << "Fitted Parameters: J/PSI" << std::endl << std::endl; + + for (const auto &par : roofit_hist_jpsi_fitsum.fitted_params) + { + res_file << par.ToString(true).c_str() << std::endl; + } + + res_file << std::endl << "Fitted Parameters: PSI(2S)" << std::endl << std::endl; + + for (const auto &par : roofit_hist_psi2s_fitsum.fitted_params) + { + res_file << par.ToString(true).c_str() << std::endl; + } + + + auto print_table = [&res_file](std::string name, std::pair sig, std::pair bkg) + { + res_file << std::endl; + res_file << "# " << name << std::endl; + res_file << "\\begin{tabular}{c|c}" << std::endl; + res_file << "$N_{Sig}$ & $N_{Bkg}$\\\\\\hline" << std::endl; + res_file << TString::Format("$%.0f \\pm %.0f$ & $%.0f \\pm %.0f$", sig.first, sig.second, bkg.first, bkg.second).Data() << std::endl; + res_file << "\\end{tabular}" << std::endl; + }; + + print_table("J/psi", roofit_hist_jpsi_fitsum.signal_yield, roofit_hist_jpsi_fitsum.background_yield); + print_table("psi(2S)", roofit_hist_psi2s_fitsum.signal_yield, roofit_hist_psi2s_fitsum.background_yield); res_file.close(); diff --git a/new_analysis_bu2jpsikplus.cpp b/new_analysis_bu2jpsikplus.cpp index 18911a8..83634fa 100644 --- a/new_analysis_bu2jpsikplus.cpp +++ b/new_analysis_bu2jpsikplus.cpp @@ -57,11 +57,13 @@ int new_analysis_bu2jpsikplus() FourVect *l24v_data = FourVect::Init(data_chain, "muminus"); FourVect *hp4v_data = FourVect::Init(data_chain, "Kplus"); Bool_t muplus_ISMUON, muminus_ISMUON, Hlt1TrackMVADecision, Hlt1TwoTrackMVADecision; + Double_t Kplus_PID_K; data_chain->SetBranchAddress("muplus_ISMUON", &muplus_ISMUON); data_chain->SetBranchAddress("muminus_ISMUON", &muminus_ISMUON); data_chain->SetBranchAddress("Hlt1TrackMVADecision", &Hlt1TrackMVADecision); data_chain->SetBranchAddress("Hlt1TwoTrackMVADecision", &Hlt1TwoTrackMVADecision); + data_chain->SetBranchAddress("Kplus_PID_K", &Kplus_PID_K); 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"); @@ -91,6 +93,7 @@ int new_analysis_bu2jpsikplus() 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_dimuon_mass = new TH1D("h1_dimuon_mass", TString::Format("#mu#mu Mass (%s)", data_tree_name), B_MASS_HIST_BINS, 2890., 3890.); 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.); @@ -139,8 +142,9 @@ int new_analysis_bu2jpsikplus() h1_B_Mass_unf->Fill(reconstructed_B_Mass); - if (muplus_ISMUON && muminus_ISMUON && (Hlt1TrackMVADecision || Hlt1TwoTrackMVADecision)) + if (muplus_ISMUON && muminus_ISMUON && Kplus_PID_K > -3 && (Hlt1TrackMVADecision || Hlt1TwoTrackMVADecision)) { + h1_dimuon_mass->Fill(dimuon.M()); if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) { B_Mass_jpsi_var = reconstructed_B_Mass; @@ -156,24 +160,82 @@ int new_analysis_bu2jpsikplus() 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); + DrawInDefaultCanvas(h1_dimuon_mass, 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.); + 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, false, 5150., 5460.); + + roofit_hist_sim.shape_parameters.sigma_lr = roofit_hist_jpsi_fitsum.shape_parameters.sigma_lr; + + //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, true, 5140., 5460.); DrawInDefaultCanvas(roofit_hist_jpsi_fitsum, analysis_name); - DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name); + //DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name); DrawInDefaultCanvas(roofit_hist_sim, analysis_name); // DrawHlt1DecisionHistos(analysis_name, hlt1_decision_histos); + //auto signal_ratio = DivWithErr(roofit_hist_psi2s_fitsum.signal_yield.first, roofit_hist_psi2s_fitsum.signal_yield.second, roofit_hist_jpsi_fitsum.signal_yield.first, roofit_hist_jpsi_fitsum.signal_yield.second); + + std::time_t t = std::time(nullptr); + std::tm tm = *std::localtime(&t); + + ofstream res_file; + res_file.open(TString::Format("%s_results.txt", analysis_name).Data(), ios::out | ios::trunc); + + auto jpsi_sigobkg = DivWithErr(roofit_hist_jpsi_fitsum.signal_yield.first, roofit_hist_jpsi_fitsum.signal_yield.second, roofit_hist_jpsi_fitsum.background_yield.first, roofit_hist_jpsi_fitsum.background_yield.second); + // auto psi2s_sigobkg = DivWithErr(roofit_hist_psi2s_fitsum.signal_yield.first, roofit_hist_psi2s_fitsum.signal_yield.second, roofit_hist_psi2s_fitsum.background_yield.first, roofit_hist_psi2s_fitsum.background_yield.second); + + auto mumu_br_frac = DivWithErr(BRF_JPSI_MUMU_VAL, BRF_JPSI_MUMU_ERR, BRF_PSI2S_MUMU_VAL, BRF_PSI2S_MUMU_ERR); + + res_file << "#### " << analysis_name << " @ " << std::put_time(&tm, "%c") << " ####" << std::endl; + res_file << "J/Psi Mode: " << ErrToStr(roofit_hist_jpsi_fitsum.signal_yield, 0) << " / " << ErrToStr(roofit_hist_jpsi_fitsum.background_yield, 0) << std::endl; + res_file << " Sig/Bkg: " << ErrToStr(jpsi_sigobkg, 2) << std::endl; + // res_file << "Psi(2S) Mode: " << ErrToStr(roofit_hist_psi2s_fitsum.signal_yield, 0) << " / " << ErrToStr(roofit_hist_psi2s_fitsum.background_yield, 0) << std::endl; + //res_file << " Sig/Bkg: " << ErrToStr(psi2s_sigobkg, 2) << std::endl; + // res_file << "Mode Yield Ratio: " << ErrToStr(signal_ratio, 3) << std::endl; + // res_file << "Rel Br Frac MuMu: " << ErrToStr(mumu_br_frac, 3) << std::endl; + + // auto br_frac = MultWithErr(signal_ratio.first, signal_ratio.second, mumu_br_frac.first, mumu_br_frac.second); + // res_file << "Rel Br Frac: " << ErrToStr(br_frac, 3) << std::endl << std::endl; + + res_file << "Params from Sim:" << std::endl << roofit_hist_sim.shape_parameters.ToString() << std::endl; + + for (const auto &par : roofit_hist_sim.fitted_params) + { + res_file << par.ToString(true).c_str() << std::endl; + } + + res_file << std::endl << "Fitted Parameters: J/PSI" << std::endl << std::endl; + + for (const auto &par : roofit_hist_jpsi_fitsum.fitted_params) + { + res_file << par.ToString(true).c_str() << std::endl; + } + + // res_file << std::endl << "Fitted Parameters: PSI(2S)" << std::endl << std::endl; + + // for (const auto &par : roofit_hist_psi2s_fitsum.fitted_params) + // { + // res_file << par.ToString(true).c_str() << std::endl; + // } + + + auto print_table = [&res_file](std::string name, std::pair sig, std::pair bkg) + { + res_file << std::endl; + res_file << "# " << name << std::endl; + res_file << "\\begin{tabular}{c|c}" << std::endl; + res_file << "$N_{Sig}$ & $N_{Bkg}$\\\\\\hline" << std::endl; + res_file << TString::Format("$%.0f \\pm %.0f$ & $%.0f \\pm %.0f$", sig.first, sig.second, bkg.first, bkg.second).Data() << std::endl; + res_file << "\\end{tabular}" << std::endl; + }; + + print_table("J/psi", roofit_hist_jpsi_fitsum.signal_yield, roofit_hist_jpsi_fitsum.background_yield); + // print_table("psi(2S)", roofit_hist_psi2s_fitsum.signal_yield, roofit_hist_psi2s_fitsum.background_yield); + + res_file.close(); return 0; } \ No newline at end of file diff --git a/new_analysis_bu2kpmumu.cpp b/new_analysis_bu2kpmumu.cpp index cd58590..b034d6d 100644 --- a/new_analysis_bu2kpmumu.cpp +++ b/new_analysis_bu2kpmumu.cpp @@ -3,7 +3,7 @@ #include #include #include -#include #include #include @@ -77,19 +77,19 @@ int new_analysis_bu2kpmumu() 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)); + TTree *tree_B_Mass_jpsi = new TTree("tree_B_Mass_jpsi", TString::Format("B^{#pm} Mass, J/#psi Mode (%s)", analysis_name)); + TTree *tree_B_Mass_psi2s = new TTree("tree_B_Mass_psi2s", TString::Format("B^{#pm} Mass, #psi(2S) Mode (%s)", analysis_name)); + TTree *tree_B_Mass_sim = new TTree("tree_B_Mass_sim", TString::Format("B^{#pm} 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); + TH1D *h1_B_Mass_unf = new TH1D("h1_B_Mass_unf", TString::Format("B^{#pm} Mass (%s), Unfiltered", analysis_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^{#pm} 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.); + TH2D *h2_Hlt1_flags_B_Mass = new TH2D("h2_Hlt1_flags_B_Mass", "Hlt1 Decision vs B^{#pm} 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^{#pm} 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); @@ -159,7 +159,10 @@ int new_analysis_bu2kpmumu() 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); + + roofit_hist_sim.shape_parameters.sigma_lr = roofit_hist_jpsi_fitsum.shape_parameters.sigma_lr; + + 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, true); DrawInDefaultCanvas(roofit_hist_jpsi_fitsum, analysis_name); DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name); @@ -168,20 +171,66 @@ int new_analysis_bu2kpmumu() // DrawHlt1DecisionHistos(analysis_name, hlt1_decision_histos); auto signal_ratio = DivWithErr(roofit_hist_psi2s_fitsum.signal_yield.first, roofit_hist_psi2s_fitsum.signal_yield.second, roofit_hist_jpsi_fitsum.signal_yield.first, roofit_hist_jpsi_fitsum.signal_yield.second); - std::time_t t = std::time(nullptr); + std::time_t t = std::time(nullptr); std::tm tm = *std::localtime(&t); ofstream res_file; res_file.open(TString::Format("%s_results.txt", analysis_name).Data(), ios::out | ios::trunc); - res_file << "#### " << std::put_time(&tm, "%c") << " ####"<< std::endl; - res_file << "J/Psi Mode: " << roofit_hist_jpsi_fitsum.signal_yield.first << " +- " << roofit_hist_jpsi_fitsum.signal_yield.second << std::endl; - res_file << "Psi(2S) Mode: " << roofit_hist_psi2s_fitsum.signal_yield.first << " +- " << roofit_hist_psi2s_fitsum.signal_yield.second << std::endl; - res_file << "Mode Yield Ratio: " << signal_ratio.first << " +- " << signal_ratio.second << std::endl; - res_file << "Rel Br Frac MuMu: " << (BRF_JPSI_MUMU_VAL / BRF_PSI2S_MUMU_VAL) << std::endl; - res_file << "Rel Br Frac: " << signal_ratio.first * (BRF_JPSI_MUMU_VAL / BRF_PSI2S_MUMU_VAL) << std::endl; + auto jpsi_sigobkg = DivWithErr(roofit_hist_jpsi_fitsum.signal_yield.first, roofit_hist_jpsi_fitsum.signal_yield.second, roofit_hist_jpsi_fitsum.background_yield.first, roofit_hist_jpsi_fitsum.background_yield.second); + auto psi2s_sigobkg = DivWithErr(roofit_hist_psi2s_fitsum.signal_yield.first, roofit_hist_psi2s_fitsum.signal_yield.second, roofit_hist_psi2s_fitsum.background_yield.first, roofit_hist_psi2s_fitsum.background_yield.second); + + auto mumu_br_frac = DivWithErr(BRF_JPSI_MUMU_VAL, BRF_JPSI_MUMU_ERR, BRF_PSI2S_MUMU_VAL, BRF_PSI2S_MUMU_ERR); + + res_file << "#### " << analysis_name << " @ " << std::put_time(&tm, "%c") << " ####" << std::endl; + res_file << "J/Psi Mode: " << ErrToStr(roofit_hist_jpsi_fitsum.signal_yield, 0) << " / " << ErrToStr(roofit_hist_jpsi_fitsum.background_yield, 0) << std::endl; + res_file << " Sig/Bkg: " << ErrToStr(jpsi_sigobkg, 2) << std::endl; + res_file << "Psi(2S) Mode: " << ErrToStr(roofit_hist_psi2s_fitsum.signal_yield, 0) << " / " << ErrToStr(roofit_hist_psi2s_fitsum.background_yield, 0) << std::endl; + res_file << " Sig/Bkg: " << ErrToStr(psi2s_sigobkg, 2) << std::endl; + res_file << "Mode Yield Ratio: " << ErrToStr(signal_ratio, 3) << std::endl; + res_file << "Rel Br Frac MuMu: " << ErrToStr(mumu_br_frac, 3) << std::endl; + + auto br_frac = MultWithErr(signal_ratio.first, signal_ratio.second, mumu_br_frac.first, mumu_br_frac.second); + res_file << "Rel Br Frac: " << ErrToStr(br_frac, 3) << std::endl << std::endl; + + res_file << "Params from Sim:" << std::endl << roofit_hist_sim.shape_parameters.ToString() << std::endl; + + for (const auto &par : roofit_hist_sim.fitted_params) + { + res_file << par.ToString(true).c_str() << std::endl; + } + + res_file << std::endl << "Fitted Parameters: J/PSI" << std::endl << std::endl; + + for (const auto &par : roofit_hist_jpsi_fitsum.fitted_params) + { + res_file << par.ToString(true).c_str() << std::endl; + } + + res_file << std::endl << "Fitted Parameters: PSI(2S)" << std::endl << std::endl; + + for (const auto &par : roofit_hist_psi2s_fitsum.fitted_params) + { + res_file << par.ToString(true).c_str() << std::endl; + } + + + auto print_table = [&res_file](std::string name, std::pair sig, std::pair bkg) + { + res_file << std::endl; + res_file << "# " << name << std::endl; + res_file << "\\begin{tabular}{c|c}" << std::endl; + res_file << "$N_{Sig}$ & $N_{Bkg}$\\\\\\hline" << std::endl; + res_file << TString::Format("$%.0f \\pm %.0f$ & $%.0f \\pm %.0f$", sig.first, sig.second, bkg.first, bkg.second).Data() << std::endl; + res_file << "\\end{tabular}" << std::endl; + }; + + print_table("J/psi", roofit_hist_jpsi_fitsum.signal_yield, roofit_hist_jpsi_fitsum.background_yield); + print_table("psi(2S)", roofit_hist_psi2s_fitsum.signal_yield, roofit_hist_psi2s_fitsum.background_yield); res_file.close(); + std::cout << "MC Truth Matched Sim Events:" << tree_B_Mass_sim->GetEntries() << std::endl; + return 0; } \ No newline at end of file diff --git a/rounder.cpp b/rounder.cpp new file mode 100644 index 0000000..f5717b3 --- /dev/null +++ b/rounder.cpp @@ -0,0 +1,78 @@ +#include + +#include "TMath.h" + +void println(std::string str) { + std::cout << str << std::endl; +} + +int get_ten_pow(double number) { + int sign = 0; + int pow = 0; + if (number >= 1) { + sign = 1; + while (number >= 1) { + // std::cout << number << " " << pow << std::endl; + number = number / 10; + pow++; + } + } else { + sign = -1; + while (number < 1) { + number = number * 10; + pow++; + } + } + + return sign * pow; +} + +int get_decimals(double number) { + int dec = 0; + while ((int) number != number) { + number = number * 10; + dec++; + } + + return dec; +} + +int rounder() { + // 2,543 ± 0,235 174,032 ± 82,543 + double val = 1.2456;// -> 4300 + double err = 0.2231; // -> 240 + + std::cout << get_ten_pow(231.4345) << " " << get_ten_pow(0.0045) << std::endl; + + int mv = 0; + if (err < 1) { + mv = 1; + } else { + mv = -1; + } + + int pot = 0; + while((mv == 1 && err < 3) || (mv == -1 && err > 30)) { + err = err * TMath::Power(10, mv); + pot++; + } + + int err_ten_pow = get_ten_pow(err); + int val_ten_pow = get_ten_pow(val); + + int diff = TMath::Abs(err_ten_pow - val_ten_pow); + + double input_err = err * TMath::Power(10, -mv * pot); + double round_err = TMath::Ceil(err) * TMath::Power(10, -mv * pot); + + int err_dec = get_decimals(round_err); + diff += err_dec; + + double mv_val = val * TMath::Power(10, mv * diff); + + double round_val = TMath::Ceil(mv_val) * TMath::Power(10, -mv * diff); + + std::cout << val << " " << mv_val << " " << err << " " << input_err<< " " << std::endl; + std::cout << val << "+-" << input_err << " -> " << round_val << "+-" << round_err << std::endl; + return 0; +} \ No newline at end of file diff --git a/tmva_plotting.cpp b/tmva_plotting.cpp new file mode 100644 index 0000000..ac369b4 --- /dev/null +++ b/tmva_plotting.cpp @@ -0,0 +1,138 @@ +#include + +#include "TMVA/efficiencies.h" + +#include "TH2F.h" +#include "TFile.h" +#include "TIterator.h" +#include "TKey.h" + +void set_tmva_style() +{ + TStyle *TMVAStyle = gROOT->GetStyle("TMVA"); + if (TMVAStyle != 0) + { + gROOT->SetStyle("TMVA"); + return; + } + + TMVAStyle = new TStyle(*gROOT->GetStyle("Plain")); // our style is based on Plain + TMVAStyle->SetName("TMVA"); + TMVAStyle->SetTitle("TMVA style based on \"Plain\" with modifications defined in tmvaglob.C"); + gROOT->GetListOfStyles()->Add(TMVAStyle); + gROOT->SetStyle("TMVA"); + + TMVAStyle->SetLineStyleString(5, "[52 12]"); + TMVAStyle->SetLineStyleString(6, "[22 12]"); + TMVAStyle->SetLineStyleString(7, "[22 10 7 10]"); + + // the pretty color palette of old + TMVAStyle->SetPalette((TMVA::gConfig().fVariablePlotting.fUsePaperStyle ? 18 : 1), 0); + + TMVAStyle->SetOptFile(0); + + // use plain black on white colors + TMVAStyle->SetFrameBorderMode(0); + TMVAStyle->SetCanvasBorderMode(0); + TMVAStyle->SetPadBorderMode(0); + TMVAStyle->SetPadColor(0); + TMVAStyle->SetFillStyle(0); + + TMVAStyle->SetLegendBorderSize(0); + + // title properties + // TMVAStyle->SetTitleW(.4); + // TMVAStyle->SetTitleH(.10); + // MVAStyle->SetTitleX(.5); + // TMVAStyle->SetTitleY(.9); + TMVAStyle->SetTitleFillColor(TColor::GetColor("#ffffff")); + TMVAStyle->SetTitleTextColor(TColor::GetColor("#000000")); + TMVAStyle->SetTitleBorderSize(0); + TMVAStyle->SetLineColor(TMVA::TMVAGlob::getTitleBorder()); + if (!TMVA::gConfig().fVariablePlotting.fUsePaperStyle) + { + TMVAStyle->SetFrameFillColor(TColor::GetColor("#ffffff")); + TMVAStyle->SetCanvasColor(TColor::GetColor("#ffffff")); + } + + // set the paper & margin sizes + TMVAStyle->SetPaperSize(20, 26); + TMVAStyle->SetPadTopMargin(0.08); + TMVAStyle->SetPadRightMargin(0.01); + TMVAStyle->SetPadBottomMargin(0.01); + TMVAStyle->SetPadLeftMargin(0.12); + + // use bold lines and markers + TMVAStyle->SetMarkerStyle(21); + TMVAStyle->SetMarkerSize(0.3); + TMVAStyle->SetHistLineWidth(2); + TMVAStyle->SetLineStyleString(2, "[12 12]"); // postscript dashes + + // do not display any of the standard histogram decorations + TMVAStyle->SetOptTitle(1); + TMVAStyle->SetTitleH(0.052); + + TMVAStyle->SetOptStat(0); + TMVAStyle->SetOptFit(0); + + // put tick marks on top and RHS of plots + TMVAStyle->SetPadTickX(1); + TMVAStyle->SetPadTickY(1); +} + +void set_title_vis(bool vis) +{ + TStyle *TMVAStyle = gROOT->GetStyle("TMVA"); + TMVAStyle->SetTitleTextColor(TColor::GetColor(vis ? "#000000" : "#ffffff")); +} + +void find_and_resize_canvas(const char *name, const char *save_folder_name, const char *save_file_name) +{ + auto canvs = (TCanvas *)gROOT->FindObject(name); + if (canvs != 0) + { + std::cout << "#### Resizing Canvas: " << canvs->GetTitle() << std::endl; + canvs->SetCanvasSize(canvs->GetWindowWidth() * 1.5, canvs->GetWindowHeight() * 1.5); + canvs->Update(); + canvs->SaveAs(TString::Format("output_files/tmva/%s/%s.pdf", save_folder_name, save_file_name)); + } +} + +void run(const char *trigger) { + std::filesystem::create_directory(TString::Format("output_files/tmva/%s", trigger).Data()); + + TString input_file = TString::Format("%s_tmva_out.root", trigger); + TString directory = TString::Format("%s_dataloader", trigger); + + set_title_vis(true); + TMVA::efficiencies(directory, input_file, 2, true); + find_and_resize_canvas("c", trigger, "efficiencies"); + + TMVA::mvas(directory, input_file, TMVA::kMVAType, true); + find_and_resize_canvas("canvas1", trigger, "class_output_testsamp"); + + TMVA::mvas(directory, input_file, TMVA::kCompareType, true); + find_and_resize_canvas("canvas1", trigger, "class_output_overtrain"); + + set_title_vis(false); + TMVA::variables(directory, input_file, "InputVariables_Id", "Input variables Id-transformed (training sample)", false, true); + for (int i = 0; i < 3; i++) + { + find_and_resize_canvas(TString::Format("canvas%d", i + 1), trigger, TString::Format("variables%d", i + 1)); + } +} + + +int tmva_plotting() +{ + set_tmva_style(); + run("BuToHpMuMu"); + run("B0ToHpHmMuMu"); + + return 0; +} + +/* +//title properties + +*/ \ No newline at end of file