pfeiffer
/
inclusive_detached_dilepton
1
0
Fork 0
Code Issues Pull Requests 1 Projects Releases Wiki Activity
Browse Source

20240403

master
Marius Pfeiffer 6 months ago
parent
9a73561ed4
commit
06a90384eb
18 changed files with 1123 additions and 158 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 66
      B0ToHpHmMuMu_results.txt
  2. 64
      B0ToKpPimMuMu_results.txt
  3. 70
      BuToHpMuMu_results.txt
  4. 37
      BuToJpsiKplus_results.txt
  5. 64
      BuToKpMuMu_results.txt
  6. 167
      basic_analysis.h
  7. 4
      batch_run_all_anas.sh
  8. 4
      bdt_classification.h
  9. 2
      constants.h
  10. 85
      mapmc_b02hphmmumu.cpp
  11. 71
      mapmc_bu2hpmumu.cpp
  12. 91
      new_analysis_b02hphmmumu.cpp
  13. 77
      new_analysis_b02kppimmumu.cpp
  14. 90
      new_analysis_bu2hpmumu.cpp
  15. 82
      new_analysis_bu2jpsikplus.cpp
  16. 79
      new_analysis_bu2kpmumu.cpp
  17. 78
      rounder.cpp
  18. 138
      tmva_plotting.cpp

66
B0ToHpHmMuMu_results.txt
View File

@ -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}

64
B0ToKpPimMuMu_results.txt
View File

@ -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}

70
BuToHpMuMu_results.txt
View File

@ -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}

37
BuToJpsiKplus_results.txt
View File

@ -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}

64
BuToKpMuMu_results.txt
View File

@ -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}

167
basic_analysis.h
View File

@ -21,6 +21,50 @@ std::pair<double, double> 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)))); return std::make_pair(x / y, TMath::Sqrt((TMath::Sq(err_x) + TMath::Sq(err_y))));
} }
std::pair<double, double> 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<double>::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<double>::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<double>::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<double, double> val_err, int n) {
return ErrToStr(val_err.first, val_err.second, n);
}
std::pair<double, double> 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 class FourVect
{ {
private: private:
@ -84,14 +128,18 @@ struct FittedParam
bool at_limit; bool at_limit;
bool constant; bool constant;
int decimals; 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->title = var.GetTitle();
this->name = var.GetName(); this->name = var.GetName();
double val = var.getVal(); double val = var.getVal();
this->value = val; this->value = val;
this->constant = var.isConstant(); this->constant = var.isConstant();
this->print_const = print_const;
this->ignore_print = ignore_print;
if (!this->constant) if (!this->constant)
{ {
this->error = var.getError(); this->error = var.getError();
@ -104,7 +152,7 @@ struct FittedParam
this->decimals = decimals; 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->title = title;
this->name = name; this->name = name;
@ -113,17 +161,23 @@ struct FittedParam
this->decimals = decimals; this->decimals = decimals;
this->constant = false; this->constant = false;
this->at_limit = false; this->at_limit = false;
this->ignore_print = ignore_print;
} }
std::string ToString() const
std::string ToString(bool latex = false) const
{ {
if (this->constant) 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 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 alpha_right;
double n_right; double n_right;
double sigma_lr; 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 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()); c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data());
} }
void DrawInDefaultCanvasStacked(std::vector<TH1D *> histograms, std::vector<Color_t> colors, std::vector<Style_t> fill_style, const char *folder, double margin_left = 0.1, Option_t *option = "")
void DrawInDefaultCanvasStacked(std::vector<TH1D *> histograms, std::vector<Color_t> colors, std::vector<Style_t> 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()); std::filesystem::create_directory(TString::Format("output_files/analysis/%s", folder).Data());
TString name = TString::Format("%s_stack_canvas", histograms[0]->GetName()); TString name = TString::Format("%s_stack_canvas", histograms[0]->GetName());
TCanvas *c = new TCanvas(name, histograms[0]->GetName(), 0, 0, 800, 600); TCanvas *c = new TCanvas(name, histograms[0]->GetName(), 0, 0, 800, 600);
c->SetLeftMargin(margin_left); 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++) for (size_t i = 0; i < histograms.size(); i++)
{ {
const Double_t scaling_factor = 1.; const Double_t scaling_factor = 1.;
auto hist_clone = (TH1 *)histograms[i]->Clone(TString::Format("%s_clone", histograms[i]->GetName())); 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->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->SetMinimum(0.);
hist_clone->SetStats(0); hist_clone->SetStats(0);
if (fill_style[i] != 0) if (fill_style[i] != 0)
@ -185,15 +256,17 @@ void DrawInDefaultCanvasStacked(std::vector<TH1D *> histograms, std::vector<Colo
if (i > 0) if (i > 0)
{ {
hist_clone->Draw(drwopt_2.c_str());
stack->Add(hist_clone, drwopt_2.c_str());
} }
else 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->Draw();
c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data()); 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_title_size = (pull_title_size * dwPad_area) / upPad_area;
Double_t fit_label_size = (pull_label_size * dwPad_area) / upPad_area; Double_t fit_label_size = (pull_label_size * dwPad_area) / upPad_area;
// canvas->Update();
c->cd(1); 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_Xaxis = fitSummary.fit_histogram->GetXaxis();
auto fit_Yaxis = fitSummary.fit_histogram->GetYaxis(); auto fit_Yaxis = fitSummary.fit_histogram->GetYaxis();
fit_Xaxis->SetLabelOffset(0.02); fit_Xaxis->SetLabelOffset(0.02);
@ -298,7 +359,7 @@ void DrawInDefaultCanvas(RooFitSummary fitSummary, const char *folder)
for (const auto &par : fitSummary.fitted_params) 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(), ""); 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) RooPlot *CreateRooFitHistogram(TH1D *hist)
{ {
RooRealVar roo_var_mass("var_mass", "B Mass Variable", B_MASS_VAR_MIN, B_MASS_VAR_MAX); 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)); 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); 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); RooRealVar roo_var_mass(var_name, "B Mass Variable", B_MASS_VAR_MIN, B_MASS_VAR_MAX);
const char *fitting_range_name = "fitting_range"; const char *fitting_range_name = "fitting_range";
roo_var_mass.setRange(fitting_range_name, fitRangeUp, fitRangeLow); 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"); TString dataset_name = suffix_name("roodataset_B_M");
// RooDataHist roodataset_B_M(hist_name, "B Mass Histogram", roo_var_mass, RooFit::Import(*hist)); // 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)); 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)); 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(); 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)); // 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_val = var_mass_nsig.getVal();
sig_err = var_mass_nsig.getError(); sig_err = var_mass_nsig.getError();
bkg_val = var_mass_nbkg.getVal(); bkg_val = var_mass_nbkg.getVal();
bkg_err = var_mass_nbkg.getError(); 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)); fitted_params.push_back(FittedParam(var_mass_bkg_c, 5));
@ -508,14 +587,14 @@ RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString
} }
else 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)); // 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; pull_compare_name = signal_name;
} }
fitted_params.push_back(FittedParam(var_mass_x0, 2)); 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_alphaL, 2));
fitted_params.push_back(FittedParam(var_mass_nL, 2)); fitted_params.push_back(FittedParam(var_mass_nL, 2));
fitted_params.push_back(FittedParam(var_mass_alphaR, 2)); fitted_params.push_back(FittedParam(var_mass_alphaR, 2));

4
batch_run_all_anas.sh
View File

@ -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

4
bdt_classification.h
View File

@ -205,7 +205,7 @@ void TrainBDT(std::vector<TV *> vars, const char* unique_id, TTree *sig_tree, TT
TString outfile_name = TString::Format("%s_tmva_out.root", unique_id); TString outfile_name = TString::Format("%s_tmva_out.root", unique_id);
TFile *output_file = TFile::Open(outfile_name, "RECREATE"); 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::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)); TMVA::DataLoader *data_loader = new TMVA::DataLoader(TString::Format("%s_dataloader", unique_id));
@ -227,7 +227,7 @@ void TrainBDT(std::vector<TV *> vars, const char* unique_id, TTree *sig_tree, TT
data_loader->AddSignalTree(sig_tree, signal_weight); data_loader->AddSignalTree(sig_tree, signal_weight);
data_loader->AddBackgroundTree(bkg_tree, background_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"); 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");

2
constants.h
View File

@ -5,6 +5,8 @@
const Double_t B_MASS_VAR_MIN = 5100.; const Double_t B_MASS_VAR_MIN = 5100.;
const Double_t B_MASS_VAR_MAX = 6000.; 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 Int_t B_MASS_HIST_BINS = 70;
const Double_t K_MASS = 493.677; const Double_t K_MASS = 493.677;

85
mapmc_b02hphmmumu.cpp
View File

@ -44,7 +44,7 @@ int mapmc_b02hphmmumu()
const char *sim_tree_name = "B0ToHpHmMuMu_noPID"; const char *sim_tree_name = "B0ToHpHmMuMu_noPID";
TChain *sim_chain = new TChain(TString::Format("%s/DecayTree", sim_tree_name)); 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 *l14v_sim = FourVect::Init(sim_chain, "L1");
FourVect *l24v_sim = FourVect::Init(sim_chain, "L2"); 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("Hm_TRUEID", &Hm_TRUEID);
sim_chain->SetBranchAddress("B0_BKGCAT", &B_BKGCAT); 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_PT", &B0_PT_in);
sim_chain->SetBranchAddress("B0_BPVFDCHI2", &B0_BPVFDCHI2_in); sim_chain->SetBranchAddress("B0_BPVFDCHI2", &B0_BPVFDCHI2_in);
sim_chain->SetBranchAddress("B0_BPVDIRA", &B0_BPVDIRA_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_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_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_BPVIPCHI2", &Hp_BPVIPCHI2_in);
sim_chain->SetBranchAddress("Hp_PT", &Hp_PT_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_BPVIPCHI2", &Hm_BPVIPCHI2_in);
sim_chain->SetBranchAddress("Hm_PT", &Hm_PT_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_BPVIPCHI2", &L1_BPVIPCHI2_in);
sim_chain->SetBranchAddress("L1_PT", &L1_PT_in);
sim_chain->SetBranchAddress("L2_BPVIPCHI2", &L2_BPVIPCHI2_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"); 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, B0_PX_out, B0_PY_out, B0_PZ_out, B0_ENERGY_out,
muplus_PX_out, muplus_PY_out, muplus_PZ_out, muplus_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, muminus_PX_out, muminus_PY_out, muminus_PZ_out, muminus_ENERGY_out,
Kplus_PX_out, Kplus_PY_out, Kplus_PZ_out, Kplus_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; 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; 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_PT", &B0_PT_out, "B0_PT/F");
output_tree->Branch("B0_BPVFDCHI2", &B0_BPVFDCHI2_out, "B0_BPVFDCHI2/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_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_PX", &B0_PX_out, "0_PX/F");
output_tree->Branch("B0_PY", &B0_PY_out, "0_PY/F"); output_tree->Branch("B0_PY", &B0_PY_out, "0_PY/F");
output_tree->Branch("B0_PZ", &B0_PZ_out, "0_PZ/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_BPVIPCHI2", &Jpsi_BPVIPCHI2_out, "Jpsi_BPVIPCHI2/F");
output_tree->Branch("Jpsi_PT", &Jpsi_PT_out, "Jpsi_PT/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_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_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_PT", &Kplus_PT_out, "Kplus_PT/F");
output_tree->Branch("Kplus_PX", &Kplus_PX_out, "Kplus_PX/F"); output_tree->Branch("Kplus_PX", &Kplus_PX_out, "Kplus_PX/F");
output_tree->Branch("Kplus_PY", &Kplus_PY_out, "Kplus_PY/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("piminus_ID", &piminus_ID_out, "piminus_ID/I");
output_tree->Branch("muminus_BPVIPCHI2", &muminus_BPVIPCHI2_out, "muminus_BPVIPCHI2/F"); 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_PX", &muminus_PX_out, "muminus_PX/F");
output_tree->Branch("muminus_PY", &muminus_PY_out, "muminus_PY/F"); output_tree->Branch("muminus_PY", &muminus_PY_out, "muminus_PY/F");
output_tree->Branch("muminus_PZ", &muminus_PZ_out, "muminus_PZ/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("muminus_ID", &muminus_ID_out, "muminus_ID/I");
output_tree->Branch("muplus_BPVIPCHI2", &muplus_BPVIPCHI2_out, "muplus_BPVIPCHI2/F"); 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_PX", &muplus_PX_out, "muplus_PX/F");
output_tree->Branch("muplus_PY", &muplus_PY_out, "muplus_PY/F"); output_tree->Branch("muplus_PY", &muplus_PY_out, "muplus_PY/F");
output_tree->Branch("muplus_PZ", &muplus_PZ_out, "muplus_PZ/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"); output_tree->Branch("muplus_ID", &muplus_ID_out, "muplus_ID/I");
unsigned int sim_entries = sim_chain->GetEntries(); unsigned int sim_entries = sim_chain->GetEntries();
unsigned int selected_entries = 0;
for (unsigned int i = 0; i < sim_entries; i++) for (unsigned int i = 0; i < sim_entries; i++)
{ {
@ -147,12 +197,23 @@ int mapmc_b02hphmmumu()
B0_PT_out = B0_PT_in; B0_PT_out = B0_PT_in;
B0_BPVFDCHI2_out = B0_BPVFDCHI2_in; B0_BPVFDCHI2_out = B0_BPVFDCHI2_in;
B0_BPVDIRA_out = B0_BPVDIRA_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_BPVIPCHI2_out = Jpsi_BPVIPCHI2_in;
Jpsi_PT_out = Jpsi_PT_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_BPVIPCHI2_out = L1_BPVIPCHI2_in;
muminus_ID_out = L1_TRUEID; muminus_ID_out = L1_TRUEID;
muminus_PT_out = L1_PT_in;
auto muminus_lv = l14v_sim->LorentzVector(); auto muminus_lv = l14v_sim->LorentzVector();
muminus_PX_out = muminus_lv.Px(); muminus_PX_out = muminus_lv.Px();
muminus_PY_out = muminus_lv.Py(); muminus_PY_out = muminus_lv.Py();
@ -162,6 +223,7 @@ int mapmc_b02hphmmumu()
muplus_BPVIPCHI2_out = L2_BPVIPCHI2_in; muplus_BPVIPCHI2_out = L2_BPVIPCHI2_in;
muplus_ID_out = L2_TRUEID; muplus_ID_out = L2_TRUEID;
muplus_PT_out = L2_PT_in;
auto muplus_lv = l24v_sim->LorentzVector(); auto muplus_lv = l24v_sim->LorentzVector();
muplus_PX_out = muplus_lv.Px(); muplus_PX_out = muplus_lv.Px();
muplus_PY_out = muplus_lv.Py(); muplus_PY_out = muplus_lv.Py();
@ -174,6 +236,7 @@ int mapmc_b02hphmmumu()
Kplus_BPVIPCHI2_out = Hp_BPVIPCHI2_in; Kplus_BPVIPCHI2_out = Hp_BPVIPCHI2_in;
Kplus_PT_out = Hp_PT_in; Kplus_PT_out = Hp_PT_in;
Kplus_PROBNN_K_out = Hp_PROBNN_K_in; Kplus_PROBNN_K_out = Hp_PROBNN_K_in;
Kplus_PID_K_out = Hp_PID_K_in;
Kplus_ID_out = Hp_TRUEID; Kplus_ID_out = Hp_TRUEID;
piminus_BPVIPCHI2_out = Hm_BPVIPCHI2_in; piminus_BPVIPCHI2_out = Hm_BPVIPCHI2_in;
piminus_PT_out = Hm_PT_in; piminus_PT_out = Hm_PT_in;
@ -208,6 +271,7 @@ int mapmc_b02hphmmumu()
Kplus_BPVIPCHI2_out = Hm_BPVIPCHI2_in; Kplus_BPVIPCHI2_out = Hm_BPVIPCHI2_in;
Kplus_PT_out = Hm_PROBNN_K_in; Kplus_PT_out = Hm_PROBNN_K_in;
Kplus_PROBNN_K_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; Kplus_ID_out = Hm_TRUEID;
piminus_BPVIPCHI2_out = Hp_BPVIPCHI2_in; piminus_BPVIPCHI2_out = Hp_BPVIPCHI2_in;
piminus_PT_out = Hp_PT_in; piminus_PT_out = Hp_PT_in;
@ -239,6 +303,7 @@ int mapmc_b02hphmmumu()
} }
output_tree->Fill(); output_tree->Fill();
selected_entries++;
} }
PrintProgress(TString::Format("%s Mapping", sim_tree_name), sim_entries, 10000, i); PrintProgress(TString::Format("%s Mapping", sim_tree_name), sim_entries, 10000, i);
@ -252,5 +317,7 @@ int mapmc_b02hphmmumu()
output_tree->Write(); output_tree->Write();
output_file->Close(); output_file->Close();
std::cout << "Triggered Entries: " << sim_entries << " Selected Entries: " << selected_entries << std::endl;
return 0; return 0;
} }

71
mapmc_bu2hpmumu.cpp
View File

@ -57,46 +57,80 @@ int mapmc_bu2hpmumu()
sim_chain->SetBranchAddress("Hp_TRUEID", &Hp_TRUEID); sim_chain->SetBranchAddress("Hp_TRUEID", &Hp_TRUEID);
sim_chain->SetBranchAddress("B_BKGCAT", &B_BKGCAT); 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_PT", &B_PT_in);
sim_chain->SetBranchAddress("B_BPVFDCHI2", &B_BPVFDCHI2_in); sim_chain->SetBranchAddress("B_BPVFDCHI2", &B_BPVFDCHI2_in);
sim_chain->SetBranchAddress("B_BPVDIRA", &B_BPVDIRA_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_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("Jpsi_PT", &Jpsi_PT_in);
sim_chain->SetBranchAddress("Hp_BPVIPCHI2", &Hp_BPVIPCHI2_in); sim_chain->SetBranchAddress("Hp_BPVIPCHI2", &Hp_BPVIPCHI2_in);
sim_chain->SetBranchAddress("Hp_PT", &Hp_PT_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_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("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"); 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,
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, B_PX_out, B_PY_out, B_PZ_out, B_ENERGY_out,
muplus_PX_out, muplus_PY_out, muplus_PZ_out, muplus_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, muminus_PX_out, muminus_PY_out, muminus_PZ_out, muminus_ENERGY_out,
Kplus_PX_out, Kplus_PY_out, Kplus_PZ_out, Kplus_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;
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; 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_BPVFDCHI2", &B_BPVFDCHI2_out, "B_BPVFDCHI2/F");
output_tree->Branch("B_BPVDIRA", &B_BPVDIRA_out, "B_BPVDIRA/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_ENERGY", &B_ENERGY_out, "B_ENERGY/F");
output_tree->Branch("B_M", &B_M_out, "B_M/D"); 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_BPVIPCHI2", &Jpsi_BPVIPCHI2_out, "Jpsi_BPVIPCHI2/F");
output_tree->Branch("Jpsi_PT", &Jpsi_PT_out, "Jpsi_PT/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_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_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_PT", &Kplus_PT_out, "Kplus_PT/F");
output_tree->Branch("Kplus_PX", &Kplus_PX_out, "Kplus_PX/F"); output_tree->Branch("Kplus_PX", &Kplus_PX_out, "Kplus_PX/F");
output_tree->Branch("Kplus_PY", &Kplus_PY_out, "Kplus_PY/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("Kplus_ID", &Kplus_ID_out, "Kplus_ID/I");
output_tree->Branch("muminus_BPVIPCHI2", &muminus_BPVIPCHI2_out, "muminus_BPVIPCHI2/F"); 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_PX", &muminus_PX_out, "muminus_PX/F");
output_tree->Branch("muminus_PY", &muminus_PY_out, "muminus_PY/F"); output_tree->Branch("muminus_PY", &muminus_PY_out, "muminus_PY/F");
output_tree->Branch("muminus_PZ", &muminus_PZ_out, "muminus_PZ/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("muminus_ID", &muminus_ID_out, "muminus_ID/I");
output_tree->Branch("muplus_BPVIPCHI2", &muplus_BPVIPCHI2_out, "muplus_BPVIPCHI2/F"); 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_PX", &muplus_PX_out, "muplus_PX/F");
output_tree->Branch("muplus_PY", &muplus_PY_out, "muplus_PY/F"); output_tree->Branch("muplus_PY", &muplus_PY_out, "muplus_PY/F");
output_tree->Branch("muplus_PZ", &muplus_PZ_out, "muplus_PZ/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"); output_tree->Branch("muplus_ID", &muplus_ID_out, "muplus_ID/I");
unsigned int sim_entries = sim_chain->GetEntries(); unsigned int sim_entries = sim_chain->GetEntries();
unsigned int selected_entries = 0;
for (unsigned int i = 0; i < sim_entries; i++) for (unsigned int i = 0; i < sim_entries; i++)
{ {
@ -131,12 +168,18 @@ int mapmc_bu2hpmumu()
B_PT_out = B_PT_in; B_PT_out = B_PT_in;
B_BPVFDCHI2_out = B_BPVFDCHI2_in; B_BPVFDCHI2_out = B_BPVFDCHI2_in;
B_BPVDIRA_out = B_BPVDIRA_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_BPVIPCHI2_out = Jpsi_BPVIPCHI2_in;
Jpsi_PT_out = Jpsi_PT_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_BPVIPCHI2_out = L1_BPVIPCHI2_in;
muminus_ID_out = L1_TRUEID; muminus_ID_out = L1_TRUEID;
muminus_PT_out = L1_PT_in;
auto muminus_lv = l14v_sim->LorentzVector(); auto muminus_lv = l14v_sim->LorentzVector();
muminus_PX_out = muminus_lv.Px(); muminus_PX_out = muminus_lv.Px();
muminus_PY_out = muminus_lv.Py(); muminus_PY_out = muminus_lv.Py();
@ -146,6 +189,7 @@ int mapmc_bu2hpmumu()
muplus_BPVIPCHI2_out = L2_BPVIPCHI2_in; muplus_BPVIPCHI2_out = L2_BPVIPCHI2_in;
muplus_ID_out = L2_TRUEID; muplus_ID_out = L2_TRUEID;
muplus_PT_out = L2_PT_in;
auto muplus_lv = l24v_sim->LorentzVector(); auto muplus_lv = l24v_sim->LorentzVector();
muplus_PX_out = muplus_lv.Px(); muplus_PX_out = muplus_lv.Px();
muplus_PY_out = muplus_lv.Py(); muplus_PY_out = muplus_lv.Py();
@ -156,6 +200,7 @@ int mapmc_bu2hpmumu()
Kplus_BPVIPCHI2_out = Hp_BPVIPCHI2_in; Kplus_BPVIPCHI2_out = Hp_BPVIPCHI2_in;
Kplus_PT_out = Hp_PT_in; Kplus_PT_out = Hp_PT_in;
Kplus_PROBNN_K_out = Hp_PROBNN_K_in; Kplus_PROBNN_K_out = Hp_PROBNN_K_in;
Kplus_PID_K_out = Hp_PID_K_in;
Kplus_ID_out = Hp_TRUEID; Kplus_ID_out = Hp_TRUEID;
auto Kplus_lv = hp4v_sim->LorentzVector(K_MASS); auto Kplus_lv = hp4v_sim->LorentzVector(K_MASS);
@ -175,6 +220,8 @@ int mapmc_bu2hpmumu()
B_M_out = B_lv.M(); B_M_out = B_lv.M();
output_tree->Fill(); output_tree->Fill();
selected_entries++;
} }
PrintProgress(TString::Format("%s Mapping", sim_tree_name), sim_entries, 10000, i); PrintProgress(TString::Format("%s Mapping", sim_tree_name), sim_entries, 10000, i);
@ -188,5 +235,7 @@ int mapmc_bu2hpmumu()
output_tree->Write(); output_tree->Write();
output_file->Close(); output_file->Close();
std::cout << "Triggered Entries: " << sim_entries << " Selected Entries: " << selected_entries << std::endl;
return 0; return 0;
} }

91
new_analysis_b02hphmmumu.cpp
View File

@ -54,12 +54,17 @@ int new_analysis_b02hphmmumu()
TChain *data_chain = new TChain(TString::Format("%s/DecayTree", data_tree_name)); 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/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 *l14v_data = FourVect::Init(data_chain, "muminus");
FourVect *l24v_data = FourVect::Init(data_chain, "muplus"); FourVect *l24v_data = FourVect::Init(data_chain, "muplus");
FourVect *hp4v_data = FourVect::Init(data_chain, "Kplus"); FourVect *hp4v_data = FourVect::Init(data_chain, "Kplus");
FourVect *hm4v_data = FourVect::Init(data_chain, "piminus"); 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)); 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"); 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_psi2s_var_name = "B_Mass_psi2s_var";
TString B_Mass_sim_var_name = "B_Mass_sim_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_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_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); 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); TH1D *h1_bdt_probs = new TH1D("h1_bdt_probs", "BDT Probabilities", 100, -1, 1);
h1_B_Mass_unf->GetXaxis()->SetTitle(end_state_mass_literal); 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); auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name);
std::map<std::string, int> exclusive_hits{}; std::map<std::string, int> exclusive_hits{};
TV* kplus_pid_k_var = TV::Double("Kplus_PID_K", "Kplus_PID_K");
std::vector<TV *> vars{ std::vector<TV *> vars{
TV::Float("B0_PT", "B0_PT"), TV::Float("B0_PT", "B0_PT"),
TV::Float("B0_BPVFDCHI2", "B0_BPVFDCHI2"), TV::Float("B0_BPVFDCHI2", "B0_BPVFDCHI2"),
TV::Float("B0_BPVDIRA", "B0_BPVDIRA"), TV::Float("B0_BPVDIRA", "B0_BPVDIRA"),
TV::Float("Jpsi_BPVIPCHI2", "Jpsi_BPVIPCHI2"), TV::Float("Jpsi_BPVIPCHI2", "Jpsi_BPVIPCHI2"),
TV::Float("Jpsi_PT", "Jpsi_PT"), 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_BPVIPCHI2", "Kplus_BPVIPCHI2"),
TV::Float("Kplus_PT", "Kplus_PT"), TV::Float("Kplus_PT", "Kplus_PT"),
TV::Float("piminus_BPVIPCHI2", "piminus_BPVIPCHI2"), TV::Float("piminus_BPVIPCHI2", "piminus_BPVIPCHI2"),
TV::Float("piminus_PT", "piminus_PT"), 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_BPVIPCHI2", "muminus_BPVIPCHI2"),
// TV::Float("muminus_PT", "muminus_PT"),
TV::Float("muplus_BPVIPCHI2", "muplus_BPVIPCHI2"), TV::Float("muplus_BPVIPCHI2", "muplus_BPVIPCHI2"),
// TV::Float("muplus_PT", "muplus_PT"),
}; };
TTree *sig_tree = new TTree("TreeS", "tree containing signal data"); 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()]; Float_t *train_vars = new Float_t[vars.size()];
auto reader = SetupReader(vars, train_vars, analysis_name); 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++) 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); 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); h1_B_Mass_bdtf->Fill(reconstructed_B_Mass);
if (TMath::Abs(reconstructed_Kstar.M() - KSTAR_MASS) < 100.) 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_unf, analysis_name, 0.1);
DrawInDefaultCanvas(h1_B_Mass_bdtf, 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_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_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); 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_jpsi_fitsum, analysis_name);
@ -280,12 +294,55 @@ int new_analysis_b02hphmmumu()
ofstream res_file; ofstream res_file;
res_file.open(TString::Format("%s_results.txt", analysis_name).Data(), ios::out | ios::trunc); 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<double, double> sig, std::pair<double, double> 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(); res_file.close();

77
new_analysis_b02kppimmumu.cpp
View File

@ -64,7 +64,7 @@ int new_analysis_b02kppimmumu()
data_chain->SetBranchAddress("piminus_Q", &Hm_Q); data_chain->SetBranchAddress("piminus_Q", &Hm_Q);
TChain *sim_chain = new TChain(TString::Format("%s/DecayTree", sim_tree_name)); 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 *l14v_sim = FourVect::Init(sim_chain, "L1");
FourVect *l24v_sim = FourVect::Init(sim_chain, "L2"); 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_psi2s_var_name = "B_Mass_psi2s_var";
TString B_Mass_sim_var_name = "B_Mass_sim_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_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_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); 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_unf->GetXaxis()->SetTitle(end_state_mass_literal);
h1_B_Mass_sim_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_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_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_jpsi_fitsum, analysis_name);
DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name); DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name);
@ -190,12 +193,56 @@ int new_analysis_b02kppimmumu()
ofstream res_file; ofstream res_file;
res_file.open(TString::Format("%s_results.txt", analysis_name).Data(), ios::out | ios::trunc); 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<double, double> sig, std::pair<double, double> 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(); res_file.close();

90
new_analysis_bu2hpmumu.cpp
View File

@ -54,11 +54,16 @@ int new_analysis_bu2hpmumu()
TChain *data_chain = new TChain(TString::Format("%s/DecayTree", data_tree_name)); 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/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 *l14v_data = FourVect::Init(data_chain, "muminus");
FourVect *l24v_data = FourVect::Init(data_chain, "muplus"); FourVect *l24v_data = FourVect::Init(data_chain, "muplus");
FourVect *hp4v_data = FourVect::Init(data_chain, "Kplus"); 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)); 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"); 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_psi2s_var_name = "B_Mass_psi2s_var";
TString B_Mass_sim_var_name = "B_Mass_sim_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_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_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); 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); TH1D *h1_bdt_probs = new TH1D("h1_bdt_probs", "BDT Probabilities", 100, -2, 2);
h1_B_Mass_unf->GetXaxis()->SetTitle(end_state_mass_literal); 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); auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name);
std::map<std::string, int> exclusive_hits{}; std::map<std::string, int> exclusive_hits{};
TV* kplus_pid_k_var = TV::Double("Kplus_PID_K", "Kplus_PID_K");
std::vector<TV *> vars{ std::vector<TV *> vars{
TV::Float("B_PT", "B_PT"), TV::Float("B_PT", "B_PT"),
TV::Float("B_BPVFDCHI2", "B_BPVFDCHI2"), TV::Float("B_BPVFDCHI2", "B_BPVFDCHI2"),
@ -103,9 +110,11 @@ int new_analysis_bu2hpmumu()
TV::Float("Jpsi_PT", "Jpsi_PT"), TV::Float("Jpsi_PT", "Jpsi_PT"),
TV::Float("Kplus_BPVIPCHI2", "Kplus_BPVIPCHI2"), TV::Float("Kplus_BPVIPCHI2", "Kplus_BPVIPCHI2"),
TV::Float("Kplus_PT", "Kplus_PT"), 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_BPVIPCHI2", "muminus_BPVIPCHI2"),
// TV::Float("muminus_PT", "muminus_PT"),
TV::Float("muplus_BPVIPCHI2", "muplus_BPVIPCHI2"), TV::Float("muplus_BPVIPCHI2", "muplus_BPVIPCHI2"),
// TV::Float("muplus_PT", "muplus_PT"),
}; };
TTree *sig_tree = new TTree("TreeS", "tree containing signal data"); 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()]; Float_t *train_vars = new Float_t[vars.size()];
auto reader = SetupReader(vars, train_vars, analysis_name); 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++) 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); 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); h1_B_Mass_bdtf->Fill(reconstructed_B_Mass);
if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) 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_unf, analysis_name, 0.1);
DrawInDefaultCanvas(h1_B_Mass_bdtf, 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_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_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); 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_jpsi_fitsum, analysis_name);
@ -270,12 +282,56 @@ int new_analysis_bu2hpmumu()
ofstream res_file; ofstream res_file;
res_file.open(TString::Format("%s_results.txt", analysis_name).Data(), ios::out | ios::trunc); 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<double, double> sig, std::pair<double, double> 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(); res_file.close();

82
new_analysis_bu2jpsikplus.cpp
View File

@ -57,11 +57,13 @@ int new_analysis_bu2jpsikplus()
FourVect *l24v_data = FourVect::Init(data_chain, "muminus"); FourVect *l24v_data = FourVect::Init(data_chain, "muminus");
FourVect *hp4v_data = FourVect::Init(data_chain, "Kplus"); FourVect *hp4v_data = FourVect::Init(data_chain, "Kplus");
Bool_t muplus_ISMUON, muminus_ISMUON, Hlt1TrackMVADecision, Hlt1TwoTrackMVADecision; Bool_t muplus_ISMUON, muminus_ISMUON, Hlt1TrackMVADecision, Hlt1TwoTrackMVADecision;
Double_t Kplus_PID_K;
data_chain->SetBranchAddress("muplus_ISMUON", &muplus_ISMUON); data_chain->SetBranchAddress("muplus_ISMUON", &muplus_ISMUON);
data_chain->SetBranchAddress("muminus_ISMUON", &muminus_ISMUON); data_chain->SetBranchAddress("muminus_ISMUON", &muminus_ISMUON);
data_chain->SetBranchAddress("Hlt1TrackMVADecision", &Hlt1TrackMVADecision); data_chain->SetBranchAddress("Hlt1TrackMVADecision", &Hlt1TrackMVADecision);
data_chain->SetBranchAddress("Hlt1TwoTrackMVADecision", &Hlt1TwoTrackMVADecision); 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)); 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"); 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); 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_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); 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_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); 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.) if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.)
{ {
B_Mass_jpsi_var = reconstructed_B_Mass; 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); 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_unf, analysis_name, 0.1);
DrawInDefaultCanvas(h1_B_Mass_sim_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_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_jpsi_fitsum, analysis_name);
DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name);
//DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name);
DrawInDefaultCanvas(roofit_hist_sim, analysis_name); DrawInDefaultCanvas(roofit_hist_sim, analysis_name);
// DrawHlt1DecisionHistos(analysis_name, hlt1_decision_histos); // 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<double, double> sig, std::pair<double, double> 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; return 0;
} }

79
new_analysis_bu2kpmumu.cpp
View File

@ -3,7 +3,7 @@
#include <cmath> #include <cmath>
#include <algorithm> #include <algorithm>
#include <filesystem> #include <filesystem>
#include <string_view
#include <string_view>
#include <ctime> #include <ctime>
#include <fstream> #include <fstream>
@ -77,19 +77,19 @@ int new_analysis_bu2kpmumu()
TString B_Mass_psi2s_var_name = "B_Mass_psi2s_var"; TString B_Mass_psi2s_var_name = "B_Mass_psi2s_var";
TString B_Mass_sim_var_name = "B_Mass_sim_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_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_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); 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_unf->GetXaxis()->SetTitle(end_state_mass_literal);
h1_B_Mass_sim_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_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_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_jpsi_fitsum, analysis_name);
DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name); DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name);
@ -174,14 +177,60 @@ int new_analysis_bu2kpmumu()
ofstream res_file; ofstream res_file;
res_file.open(TString::Format("%s_results.txt", analysis_name).Data(), ios::out | ios::trunc); 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<double, double> sig, std::pair<double, double> 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(); res_file.close();
std::cout << "MC Truth Matched Sim Events:" << tree_B_Mass_sim->GetEntries() << std::endl;
return 0; return 0;
} }

78
rounder.cpp
View File

@ -0,0 +1,78 @@
#include <iostream>
#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;
}

138
tmva_plotting.cpp
View File

@ -0,0 +1,138 @@
#include <filesystem>
#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
*/
Write Preview
Loading…
Cancel
Save