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

Merge pull request 'use_own_root_files' (#2) from use_own_root_files into master 

Reviewed-on: #2
master
Marius Pfeiffer 7 months ago
parent
84df7596bf bcb6249b55
commit
9a73561ed4
9 changed files with 264 additions and 251 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. 6
      B0ToHpHmMuMu_results.txt
  2. 6
      BuToHpMuMu_results.txt
  3. 85
      basic_analysis.h
  4. 6
      constants.h
  5. 2
      mapmc_b02hphmmumu.cpp
  6. 233
      new_analysis_b02hphmmumu.cpp
  7. 43
      new_analysis_b02kppimmumu.cpp
  8. 107
      new_analysis_bu2hpmumu.cpp
  9. 21
      new_analysis_bu2kpmumu.cpp

6
B0ToHpHmMuMu_results.txt
View File

@ -0,0 +1,6 @@
#### 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

6
BuToHpMuMu_results.txt
View File

@ -0,0 +1,6 @@
#### 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

85
basic_analysis.h
View File

@ -15,6 +15,12 @@
#include "constants.h" #include "constants.h"
std::pair<double, double> DivWithErr(double x, double dx, double y, double dy) {
double err_x = (1 / y) * dx;
double err_y = - (x / (y*y)) * dy;
return std::make_pair(x / y, TMath::Sqrt((TMath::Sq(err_x) + TMath::Sq(err_y))));
}
class FourVect class FourVect
{ {
private: private:
@ -98,7 +104,8 @@ 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)
{
this->title = title; this->title = title;
this->name = name; this->name = name;
this->value = value; this->value = value;
@ -127,6 +134,7 @@ struct ShapeParamters
double n_left; double n_left;
double alpha_right; double alpha_right;
double n_right; double n_right;
double sigma_lr;
}; };
struct RooFitSummary struct RooFitSummary
@ -134,6 +142,8 @@ struct RooFitSummary
RooPlot *fit_histogram; RooPlot *fit_histogram;
RooPlot *pull_histogram; RooPlot *pull_histogram;
std::map<std::string, std::string> pdf_names; std::map<std::string, std::string> pdf_names;
std::pair<double, double> signal_yield;
std::pair<double, double> background_yield;
std::vector<FittedParam> fitted_params; std::vector<FittedParam> fitted_params;
ShapeParamters shape_parameters; ShapeParamters shape_parameters;
}; };
@ -156,17 +166,19 @@ void DrawInDefaultCanvasStacked(std::vector<TH1D *> histograms, std::vector<Colo
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(); std::string drwopt_2 = std::string(option).empty() ? "SAME HIST" : TString::Format("%s SAME HIST", option).Data();
for (size_t i = 0; i < histograms.size(); i++) 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 / histograms[i]->GetMaximum());
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->SetMinimum(0.); hist_clone->SetMinimum(0.);
hist_clone->SetStats(0); hist_clone->SetStats(0);
if (fill_style[i] != 0) {
if (fill_style[i] != 0)
{
hist_clone->SetFillStyle(fill_style[i]); hist_clone->SetFillStyle(fill_style[i]);
hist_clone->SetFillColor(colors[i]); hist_clone->SetFillColor(colors[i]);
} }
@ -372,7 +384,7 @@ RooPlot *CreateRooFitHistogram(TH1D *hist)
return roo_frame_mass; return roo_frame_mass;
} }
RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString xAxis, bool hasExpBkg, bool useExtShape, ShapeParamters extShape, Double_t fit_low = 0, Double_t fit_up = 0)
RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString xAxis, bool hasExpBkg, bool useExtShape, ShapeParamters extShape, bool const_sigma = false, Double_t fit_low = 0, Double_t fit_up = 0)
{ {
auto suffix_name = [name = dataSet->GetName()](const char *text) auto suffix_name = [name = dataSet->GetName()](const char *text)
{ {
@ -387,7 +399,7 @@ RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString
roo_var_mass.setRange(fitting_range_name, fitRangeUp, fitRangeLow); roo_var_mass.setRange(fitting_range_name, fitRangeUp, fitRangeLow);
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));
RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title(dataSet->GetTitle()), RooFit::Name(TString::Format("%s_rplt", dataSet->GetName()))); RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title(dataSet->GetTitle()), RooFit::Name(TString::Format("%s_rplt", dataSet->GetName())));
@ -405,17 +417,35 @@ RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString
if (useExtShape) if (useExtShape)
{ {
var_mass_alphaL.setConstant(true);
var_mass_alphaL.setVal(extShape.alpha_left);
if (extShape.alpha_left != 0.)
{
var_mass_alphaL.setConstant(true);
var_mass_alphaL.setVal(extShape.alpha_left);
}
if (extShape.n_left != 0.)
{
var_mass_nL.setConstant(true);
var_mass_nL.setVal(extShape.n_left);
}
var_mass_nL.setConstant(true);
var_mass_nL.setVal(extShape.n_left);
if (extShape.alpha_right != 0.)
{
var_mass_alphaR.setConstant(true);
var_mass_alphaR.setVal(extShape.alpha_right);
}
var_mass_alphaR.setConstant(true);
var_mass_alphaR.setVal(extShape.alpha_right);
if (extShape.n_right != 0.)
{
var_mass_nR.setConstant(true);
var_mass_nR.setVal(extShape.n_right);
}
var_mass_nR.setConstant(true);
var_mass_nR.setVal(extShape.n_right);
if (extShape.sigma_lr != 0. && const_sigma)
{
var_mass_sigmaLR.setConstant(true);
var_mass_sigmaLR.setVal(extShape.sigma_lr);
}
} }
TString signal_name = suffix_name("sig_cb"); TString signal_name = suffix_name("sig_cb");
@ -427,6 +457,10 @@ RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString
pdf_names[signal_name.Data()] = sig_cb.getTitle().Data(); pdf_names[signal_name.Data()] = sig_cb.getTitle().Data();
TString pull_compare_name{}; TString pull_compare_name{};
double sig_val = 0.;
double sig_err = 0.;
double bkg_val = 0.;
double bkg_err = 0.;
if (hasExpBkg) if (hasExpBkg)
{ {
// Exponential for Background // Exponential for Background
@ -453,18 +487,20 @@ RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString
fitted_params.push_back(FittedParam(var_mass_nsig, 2)); 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_nbkg, 2));
sig_val = var_mass_nsig.getVal();
sig_err = var_mass_nsig.getError();
bkg_val = var_mass_nbkg.getVal();
bkg_err = var_mass_nbkg.getError();
double sig_val = var_mass_nsig.getVal();
double sig_err = var_mass_nsig.getError();
double bkg_val = var_mass_nbkg.getVal();
double bkg_err = var_mass_nbkg.getError();
double 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_over_bkg_val = sig_val / TMath::Sqrt(sig_val + bkg_val);
double err_prop_sig = (sig_val + 2 * bkg_val)/(2 * TMath::Power((sig_val + bkg_val), (3/2)));
double err_prop_bkg = -sig_val/(2 * TMath::Power((sig_val + bkg_val), (3/2)));
double sig_over_bkg_err = TMath::Sqrt(TMath::Sq(err_prop_sig * sig_err) + TMath::Sq(err_prop_bkg * bkg_err));
auto sig_over_bkg = DivWithErr(sig_val, sig_err, bkg_val, bkg_err);
fitted_params.push_back(FittedParam("sig_over_bkg", "N_{Sig}/#sqrt{N_{Sig} + N_{Bkg}}", sig_over_bkg_val, sig_over_bkg_err, 2));
fitted_params.push_back(FittedParam("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(var_mass_bkg_c, 5)); fitted_params.push_back(FittedParam(var_mass_bkg_c, 5));
@ -498,12 +534,15 @@ RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString
roo_frame_mass, roo_frame_mass,
roo_frame_pull, roo_frame_pull,
pdf_names, pdf_names,
std::make_pair(sig_val, sig_err),
std::make_pair(bkg_val, bkg_err),
fitted_params, fitted_params,
ShapeParamters{ ShapeParamters{
var_mass_alphaL.getVal(), var_mass_alphaL.getVal(),
var_mass_nL.getVal(), var_mass_nL.getVal(),
var_mass_alphaR.getVal(), var_mass_alphaR.getVal(),
var_mass_nR.getVal()}};
var_mass_nR.getVal(),
var_mass_sigmaLR.getVal()}};
} }
bool InRange(double value, double center, double low_intvl, double up_intvl) bool InRange(double value, double center, double low_intvl, double up_intvl)

6
constants.h
View File

@ -16,4 +16,10 @@ const int PID_KAON = 321;
const int PID_PION = 211; const int PID_PION = 211;
const int PID_MUON = 13; const int PID_MUON = 13;
const double BRF_JPSI_MUMU_VAL = 0.0593;
const double BRF_JPSI_MUMU_ERR = 0.0006;
const double BRF_PSI2S_MUMU_VAL = 0.0077;
const double BRF_PSI2S_MUMU_ERR = 0.0008;
#endif #endif

2
mapmc_b02hphmmumu.cpp
View File

@ -89,7 +89,7 @@ int mapmc_b02hphmmumu()
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 Kplus_PROBNN_K_out, piminus_PROBNN_K_out, B0_M_out, muplus_M_out, muminus_M_out, Kplus_M_out, piminus_M_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/D");
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_PX", &B0_PX_out, "0_PX/F"); output_tree->Branch("B0_PX", &B0_PX_out, "0_PX/F");

233
new_analysis_b02hphmmumu.cpp
View File

@ -4,6 +4,8 @@
#include <algorithm> #include <algorithm>
#include <filesystem> #include <filesystem>
#include <string_view> #include <string_view>
#include <ctime>
#include <fstream>
#include "TH1D.h" #include "TH1D.h"
#include "TH2D.h" #include "TH2D.h"
@ -45,39 +47,26 @@
int new_analysis_b02hphmmumu() int new_analysis_b02hphmmumu()
{ {
const char *analysis_name = "B0ToHpHmMuMu"; const char *analysis_name = "B0ToHpHmMuMu";
const char *data_tree_name = "B0ToHpHmMuMu";
const char *sim_tree_name = "B0ToHpHmMuMu_noPID";
const char *data_tree_name = "SpruceRD_B0ToHpHmMuMu";
const char *sim_tree_name = "B0ToHpHmMuMu_noPID_mapped";
const char *end_state_mass_literal = "m(#pi^{+}#pi^{-}_{(#rightarrow K^{-})}#mu^{+}#mu^{-} & #pi^{+}_{(#rightarrow K^{+})}#pi^{-}#mu^{+}#mu^{-})"; const char *end_state_mass_literal = "m(#pi^{+}#pi^{-}_{(#rightarrow K^{-})}#mu^{+}#mu^{-} & #pi^{+}_{(#rightarrow K^{+})}#pi^{-}#mu^{+}#mu^{-})";
const bool retrain_bdt = false; const bool retrain_bdt = false;
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/spruce_magdown_2023_v0_tuple_90000000_v0r0p6288631.root");
data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0r1_tuple_90000000_2023_v0r0p6288631.root");
data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_Sprucing23r1_90000000_RD.root");
Double_t Hp_PID_K, Hm_PID_K;
data_chain->SetBranchAddress("Hp_PID_K", &Hp_PID_K);
data_chain->SetBranchAddress("Hm_PID_K", &Hm_PID_K);
FourVect *l14v_data = FourVect::Init(data_chain, "L1");
FourVect *l24v_data = FourVect::Init(data_chain, "L2");
FourVect *hp4v_data = FourVect::Init(data_chain, "Hp");
FourVect *hm4v_data = FourVect::Init(data_chain, "Hm");
FourVect *l14v_data = FourVect::Init(data_chain, "muminus");
FourVect *l24v_data = FourVect::Init(data_chain, "muplus");
FourVect *hp4v_data = FourVect::Init(data_chain, "Kplus");
FourVect *hm4v_data = FourVect::Init(data_chain, "piminus");
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");
FourVect *l14v_sim = FourVect::Init(sim_chain, "L1");
FourVect *l24v_sim = FourVect::Init(sim_chain, "L2");
FourVect *hp4v_sim = FourVect::Init(sim_chain, "Hp");
FourVect *hm4v_sim = FourVect::Init(sim_chain, "Hm");
Int_t B_BKGCAT, L1_TRUEID, L2_TRUEID, Hp_TRUEID, Hm_TRUEID;
sim_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/B0ToHpHmMuMu_mapped_mc.root");
sim_chain->SetBranchAddress("L1_TRUEID", &L1_TRUEID);
sim_chain->SetBranchAddress("L2_TRUEID", &L2_TRUEID);
sim_chain->SetBranchAddress("Hp_TRUEID", &Hp_TRUEID);
sim_chain->SetBranchAddress("Hm_TRUEID", &Hm_TRUEID);
sim_chain->SetBranchAddress("B0_BKGCAT", &B_BKGCAT);
FourVect *l14v_sim = FourVect::Init(sim_chain, "muminus");
FourVect *l24v_sim = FourVect::Init(sim_chain, "muplus");
FourVect *hp4v_sim = FourVect::Init(sim_chain, "Kplus");
FourVect *hm4v_sim = FourVect::Init(sim_chain, "piminus");
Double_t B_Mass_jpsi_var, B_Mass_psi2s_var, B_Mass_sim_var; Double_t B_Mass_jpsi_var, B_Mass_psi2s_var, B_Mass_sim_var;
TString B_Mass_jpsi_var_name = "B_Mass_jpsi_var"; TString B_Mass_jpsi_var_name = "B_Mass_jpsi_var";
@ -94,14 +83,12 @@ int new_analysis_b02hphmmumu()
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_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_bdtf = new TH1D("h1_B_Mass_bdtf", TString::Format("B Mass (%s), BDT Filter", data_tree_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX);
TH1D *h1_B_Mass_sim_unf = new TH1D("h1_B_Mass_sim_unf", TString::Format("B Mass, Simualted (%s), Unfiltered", sim_tree_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX);
TH2D *h2_Hlt1_flags_B_Mass = new TH2D("h2_Hlt1_flags_B_Mass", "Hlt1 Decision vs B Mass", 50, 5100, 5400, 13, 1., 14.); TH2D *h2_Hlt1_flags_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_excl_B_Mass = new TH2D("h2_Hlt1_flags_excl_B_Mass", "Excl Hlt1 Decision vs B Mass", 50, 5100, 5400, 13, 1., 14.);
TH1D *h1_bdt_probs = new TH1D("h1_bdt_probs", "BDT Probabilities", 100, -1, 1); 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);
h1_B_Mass_sim_unf->GetXaxis()->SetTitle(end_state_mass_literal);
h1_B_Mass_bdtf->GetXaxis()->SetTitle(end_state_mass_literal); h1_B_Mass_bdtf->GetXaxis()->SetTitle(end_state_mass_literal);
h2_Hlt1_flags_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); h2_Hlt1_flags_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal);
h2_Hlt1_flags_excl_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); h2_Hlt1_flags_excl_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal);
@ -115,17 +102,14 @@ int new_analysis_b02hphmmumu()
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_BPVDIRA", "Jpsi_BPVDIRA"),
TV::Float("Jpsi_PT", "Jpsi_PT"), TV::Float("Jpsi_PT", "Jpsi_PT"),
TV::Float("Hp_BPVIPCHI2", "Hp_BPVIPCHI2"),
TV::Float("Hp_PT", "Hp_PT"),
TV::Float("Hm_BPVIPCHI2", "Hm_BPVIPCHI2"),
TV::Float("Hm_PT", "Hm_PT"),
// TV::Double("Kplus_PID_K", "K_PID_K"),
TV::Double("Hp_PROBNN_K", "Hp_PROBNN_K"),
TV::Double("Hm_PROBNN_K", "Hm_PROBNN_K"),
TV::Float("L1_BPVIPCHI2", "L1_BPVIPCHI2"),
TV::Float("L2_BPVIPCHI2", "L2_BPVIPCHI2"),
TV::Float("Kplus_BPVIPCHI2", "Kplus_BPVIPCHI2"),
TV::Float("Kplus_PT", "Kplus_PT"),
TV::Float("piminus_BPVIPCHI2", "piminus_BPVIPCHI2"),
TV::Float("piminus_PT", "piminus_PT"),
TV::Double("Kplus_PROBNN_K", "Kplus_PROBNN_K"),
TV::Float("muminus_BPVIPCHI2", "muminus_BPVIPCHI2"),
TV::Float("muplus_BPVIPCHI2", "muplus_BPVIPCHI2"),
}; };
TTree *sig_tree = new TTree("TreeS", "tree containing signal data"); TTree *sig_tree = new TTree("TreeS", "tree containing signal data");
@ -145,31 +129,17 @@ int new_analysis_b02hphmmumu()
for (unsigned int i = 0; i < data_entries; i++) for (unsigned int i = 0; i < data_entries; i++)
{ {
data_chain->GetEntry(i); data_chain->GetEntry(i);
TLorentzVector reconstructed_Kstar{};
bool found_k_star = false;
if (Hp_PID_K > 0 && Hm_PID_K < 0)
{
reconstructed_Kstar = hp4v_data->LorentzVector(K_MASS) + hm4v_data->LorentzVector();
found_k_star = true;
}
else if (Hm_PID_K > 0 && Hp_PID_K < 0)
{
reconstructed_Kstar = hp4v_data->LorentzVector() + hm4v_data->LorentzVector(K_MASS);
found_k_star = true;
}
TLorentzVector reconstructed_Kstar = hp4v_data->LorentzVector() + hm4v_data->LorentzVector();
TLorentzVector dimuon = l14v_data->LorentzVector() + l24v_data->LorentzVector();
Double_t reconstructed_B_Mass = (reconstructed_Kstar + dimuon).M();
if (found_k_star)
if (std::all_of(vars.begin(), vars.end(), [](TV *v)
{ return v->IsDataFinite(); }))
{ {
Double_t reconstructed_B_Mass = (reconstructed_Kstar + l14v_data->LorentzVector() + l24v_data->LorentzVector()).M();
if (std::all_of(vars.begin(), vars.end(), [](TV *v)
{ return v->IsDataFinite(); }))
if (reconstructed_B_Mass > 5500. && ((TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) || (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.)))
{ {
if (reconstructed_B_Mass > 5500.)
{
bkg_tree->Fill();
bkg_events++;
}
bkg_tree->Fill();
bkg_events++;
} }
} }
@ -187,50 +157,28 @@ int new_analysis_b02hphmmumu()
{ {
sim_chain->GetEntry(i); sim_chain->GetEntry(i);
Double_t reco_mass_pipkp = (hp4v_sim->LorentzVector(K_MASS) + hm4v_sim->LorentzVector() + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M();
Double_t reco_mass_pimkm = (hp4v_sim->LorentzVector() + hm4v_sim->LorentzVector(K_MASS) + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M();
h1_B_Mass_sim_unf->Fill(reco_mass_pipkp);
h1_B_Mass_sim_unf->Fill(reco_mass_pimkm);
TLorentzVector reconstructed_Kstar = hp4v_sim->LorentzVector() + hm4v_sim->LorentzVector();
Double_t reconstructed_B_Mass = (reconstructed_Kstar + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M();
if (B_BKGCAT == 30 && TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID)
if (sig_events < bkg_events)
{ {
TLorentzVector reconstructed_Kstar{};
bool found_k_star = false;
if (TMath::Abs(Hp_TRUEID) == PID_KAON && TMath::Abs(Hm_TRUEID) == PID_PION)
if (retrain_bdt && std::all_of(vars.begin(), vars.end(), [](TV *v)
{ return v->IsMCFinite(); }))
{ {
reconstructed_Kstar = hp4v_sim->LorentzVector(K_MASS) + hm4v_sim->LorentzVector();
found_k_star = true;
}
else if (TMath::Abs(Hp_TRUEID) == PID_PION && TMath::Abs(Hm_TRUEID) == PID_KAON)
{
reconstructed_Kstar = hp4v_sim->LorentzVector() + hm4v_sim->LorentzVector(K_MASS);
found_k_star = true;
}
if (found_k_star)
{
Double_t reconstructed_B_Mass = (reconstructed_Kstar + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M();
if (sig_events < bkg_events)
{
if (retrain_bdt && std::all_of(vars.begin(), vars.end(), [](TV *v)
{ return v->IsMCFinite(); }))
{
sig_tree->Fill();
sig_events++;
}
}
B_Mass_sim_var = reconstructed_B_Mass;
tree_B_Mass_sim->Fill();
sig_tree->Fill();
sig_events++;
} }
} }
B_Mass_sim_var = reconstructed_B_Mass;
tree_B_Mass_sim->Fill();
PrintProgress(TString::Format("%s SIG Collection", analysis_name), sim_entries, 10000, i); PrintProgress(TString::Format("%s SIG Collection", analysis_name), sim_entries, 10000, i);
} }
if (retrain_bdt) if (retrain_bdt)
{ {
std::cout << "# Added " << sig_events << " signal events." << std::endl;
TrainBDT(vars, analysis_name, sig_tree, bkg_tree); TrainBDT(vars, analysis_name, sig_tree, bkg_tree);
std::cout << "# Finished BDT retrain." << std::endl; std::cout << "# Finished BDT retrain." << std::endl;
} }
@ -240,72 +188,57 @@ 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.0508;
const double mva_cut_value = 0;
for (unsigned int i = 0; i < data_entries; i++) for (unsigned int i = 0; i < data_entries; i++)
{ {
data_chain->GetEntry(i); data_chain->GetEntry(i);
TLorentzVector dimuon = l14v_data->LorentzVector() + l24v_data->LorentzVector(); TLorentzVector dimuon = l14v_data->LorentzVector() + l24v_data->LorentzVector();
TLorentzVector reconstructed_Kstar{};
bool found_k_star = false;
if (Hp_PID_K > 0 && Hm_PID_K < 0)
{
reconstructed_Kstar = hp4v_data->LorentzVector(K_MASS) + hm4v_data->LorentzVector();
found_k_star = true;
}
else if (Hm_PID_K > 0 && Hp_PID_K < 0)
{
reconstructed_Kstar = hp4v_data->LorentzVector() + hm4v_data->LorentzVector(K_MASS);
found_k_star = true;
}
TLorentzVector reconstructed_Kstar = hp4v_data->LorentzVector() + hm4v_data->LorentzVector();
Double_t reconstructed_B_Mass = (reconstructed_Kstar + dimuon).M();
if (found_k_star)
if (std::all_of(vars.begin(), vars.end(), [](TV *v)
{ return v->IsDataFinite(); }))
{ {
Double_t reconstructed_B_Mass = (reconstructed_Kstar + dimuon).M();
if (std::all_of(vars.begin(), vars.end(), [](TV *v)
{ return v->IsDataFinite(); }))
for (size_t j = 0; j < vars.size(); j++)
{ {
for (size_t j = 0; j < vars.size(); j++)
if (vars[j]->IsDouble())
{
train_vars[j] = vars[j]->GetDataDouble();
}
else if (vars[j]->IsFloat())
{ {
if (vars[j]->IsDouble())
{
train_vars[j] = vars[j]->GetDataDouble();
}
else if (vars[j]->IsFloat())
{
train_vars[j] = vars[j]->GetDataFloat();
}
train_vars[j] = vars[j]->GetDataFloat();
} }
} }
}
if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.)
{
CheckHlt1Decisioins(h2_Hlt1_flags_B_Mass, h2_Hlt1_flags_excl_B_Mass, exclusive_hits, reconstructed_B_Mass);
FillHlt1DecisionHistos(hlt1_decision_histos, reconstructed_B_Mass);
}
if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.)
{
CheckHlt1Decisioins(h2_Hlt1_flags_B_Mass, h2_Hlt1_flags_excl_B_Mass, exclusive_hits, reconstructed_B_Mass);
FillHlt1DecisionHistos(hlt1_decision_histos, reconstructed_B_Mass);
}
double mva_response = reader->EvaluateMVA("BDT");
h1_bdt_probs->Fill(mva_response);
double mva_response = reader->EvaluateMVA("BDT");
h1_bdt_probs->Fill(mva_response);
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)
{
h1_B_Mass_bdtf->Fill(reconstructed_B_Mass);
if (TMath::Abs(reconstructed_Kstar.M() - KSTAR_MASS) < 100.)
{ {
h1_B_Mass_bdtf->Fill(reconstructed_B_Mass);
if (TMath::Abs(reconstructed_Kstar.M() - KSTAR_MASS) < 100.)
if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.)
{
B_Mass_jpsi_var = reconstructed_B_Mass;
tree_B_Mass_jpsi->Fill();
}
else if (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.)
{ {
if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.)
{
B_Mass_jpsi_var = reconstructed_B_Mass;
tree_B_Mass_jpsi->Fill();
}
else if (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.)
{
B_Mass_psi2s_var = reconstructed_B_Mass;
tree_B_Mass_psi2s->Fill();
}
B_Mass_psi2s_var = reconstructed_B_Mass;
tree_B_Mass_psi2s->Fill();
} }
} }
} }
@ -323,14 +256,13 @@ int new_analysis_b02hphmmumu()
DrawInDefaultCanvas(h2_Hlt1_flags_excl_B_Mass, analysis_name, 0.16, "COLZ"); DrawInDefaultCanvas(h2_Hlt1_flags_excl_B_Mass, analysis_name, 0.16, "COLZ");
DrawInDefaultCanvas(h1_B_Mass_unf, analysis_name, 0.1); DrawInDefaultCanvas(h1_B_Mass_unf, analysis_name, 0.1);
DrawInDefaultCanvas(h1_B_Mass_sim_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);
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);
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);
@ -340,5 +272,22 @@ int new_analysis_b02hphmmumu()
DrawBDTProbs(h1_bdt_probs, mva_cut_value, analysis_name); DrawBDTProbs(h1_bdt_probs, mva_cut_value, analysis_name);
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);
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;
res_file.close();
return 0; return 0;
} }

43
new_analysis_b02kppimmumu.cpp
View File

@ -4,6 +4,8 @@
#include <algorithm> #include <algorithm>
#include <filesystem> #include <filesystem>
#include <string_view> #include <string_view>
#include <ctime>
#include <fstream>
#include "TH1D.h" #include "TH1D.h"
#include "TH2D.h" #include "TH2D.h"
@ -48,7 +50,6 @@ int new_analysis_b02kppimmumu()
const char *data_tree_name = "Hlt2RD_B0ToKpPimMuMu"; const char *data_tree_name = "Hlt2RD_B0ToKpPimMuMu";
const char *sim_tree_name = "B0ToKpPimMuMu_noPID"; const char *sim_tree_name = "B0ToKpPimMuMu_noPID";
const char *end_state_mass_literal = "m(K^{+}#pi^{-}#mu^{+}#mu^{-})"; const char *end_state_mass_literal = "m(K^{+}#pi^{-}#mu^{+}#mu^{-})";
const bool retrain_bdt = true;
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_SprucingPass23r1_94000000_RD.root"); data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_RD.root");
@ -125,11 +126,6 @@ int new_analysis_b02kppimmumu()
std::cout << "# Starting evaluation of data." << std::endl; std::cout << "# Starting evaluation of data." << std::endl;
int kplus = 0;
int kminus = 0;
int piplus = 0;
int piminus = 0;
for (unsigned int i = 0; i < data_entries; i++) for (unsigned int i = 0; i < data_entries; i++)
{ {
data_chain->GetEntry(i); data_chain->GetEntry(i);
@ -149,24 +145,6 @@ int new_analysis_b02kppimmumu()
{ {
if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.)
{ {
if (Hp_Q == 1)
{
kplus++;
}
else if (Hp_Q == -1)
{
kminus++;
}
if (Hm_Q == 1)
{
piplus++;
}
else if (Hm_Q == -1)
{
piminus++;
}
B_Mass_jpsi_var = reconstructed_B_Mass; B_Mass_jpsi_var = reconstructed_B_Mass;
tree_B_Mass_jpsi->Fill(); tree_B_Mass_jpsi->Fill();
} }
@ -204,7 +182,22 @@ int new_analysis_b02kppimmumu()
// DrawHlt1DecisionHistos(analysis_name, hlt1_decision_histos); // DrawHlt1DecisionHistos(analysis_name, hlt1_decision_histos);
std::cout << "hist entries: " << tree_B_Mass_jpsi->GetEntries() << ", kplus: " << kplus << ", kminus: " << kminus << ", piplus: " << piplus << ", piminus: " << piminus << std::endl;
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);
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;
res_file.close();
return 0; return 0;
} }

107
new_analysis_bu2hpmumu.cpp
View File

@ -4,6 +4,8 @@
#include <algorithm> #include <algorithm>
#include <filesystem> #include <filesystem>
#include <string_view> #include <string_view>
#include <ctime>
#include <fstream>
#include "TH1D.h" #include "TH1D.h"
#include "TH2D.h" #include "TH2D.h"
@ -45,34 +47,24 @@
int new_analysis_bu2hpmumu() int new_analysis_bu2hpmumu()
{ {
const char *analysis_name = "BuToHpMuMu"; const char *analysis_name = "BuToHpMuMu";
const char *data_tree_name = "BuToHpMuMu";
const char *sim_tree_name = "BuToHpMuMu_noPID";
const char *data_tree_name = "SpruceRD_BuToHpMuMu";
const char *sim_tree_name = "BuToHpMuMu_noPID_mapped";
const char *end_state_mass_literal = "m(#pi^{+}_{(#rightarrow K^{+})}#mu^{+}#mu^{-})"; const char *end_state_mass_literal = "m(#pi^{+}_{(#rightarrow K^{+})}#mu^{+}#mu^{-})";
const bool retrain_bdt = false; const bool retrain_bdt = false;
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/spruce_magdown_2023_v0_tuple_90000000_v0r0p6288631.root");
data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0r1_tuple_90000000_2023_v0r0p6288631.root");
data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_Sprucing23r1_90000000_RD.root");
FourVect *l14v_data = FourVect::Init(data_chain, "L1");
FourVect *l24v_data = FourVect::Init(data_chain, "L2");
FourVect *hp4v_data = FourVect::Init(data_chain, "Hp");
FourVect *l14v_data = FourVect::Init(data_chain, "muminus");
FourVect *l24v_data = FourVect::Init(data_chain, "muplus");
FourVect *hp4v_data = FourVect::Init(data_chain, "Kplus");
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_BuToKpMuMu_12143001_magdown.root");
sim_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/BuToHpMuMu_mapped_mc.root");
FourVect *l14v_sim = FourVect::Init(sim_chain, "L1");
FourVect *l24v_sim = FourVect::Init(sim_chain, "L2");
FourVect *hp4v_sim = FourVect::Init(sim_chain, "Hp");
Int_t B_BKGCAT, L1_TRUEID, L2_TRUEID, Hp_TRUEID;
Double_t Hp_PID_K_sim;
sim_chain->SetBranchAddress("L1_TRUEID", &L1_TRUEID);
sim_chain->SetBranchAddress("L2_TRUEID", &L2_TRUEID);
sim_chain->SetBranchAddress("Hp_TRUEID", &Hp_TRUEID);
sim_chain->SetBranchAddress("B_BKGCAT", &B_BKGCAT);
sim_chain->SetBranchAddress("Hp_PID_K", &Hp_PID_K_sim);
FourVect *l14v_sim = FourVect::Init(sim_chain, "muminus");
FourVect *l24v_sim = FourVect::Init(sim_chain, "muplus");
FourVect *hp4v_sim = FourVect::Init(sim_chain, "Kplus");
Double_t B_Mass_jpsi_var, B_Mass_psi2s_var, B_Mass_sim_var; Double_t B_Mass_jpsi_var, B_Mass_psi2s_var, B_Mass_sim_var;
TString B_Mass_jpsi_var_name = "B_Mass_jpsi_var"; TString B_Mass_jpsi_var_name = "B_Mass_jpsi_var";
@ -89,16 +81,12 @@ int new_analysis_bu2hpmumu()
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_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_bdtf = new TH1D("h1_B_Mass_bdtf", TString::Format("B Mass (%s), BDT Filter", data_tree_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX);
TH1D *h1_B_Mass_sim_unf = new TH1D("h1_B_Mass_sim_unf", TString::Format("B Mass, Simualted (%s), Unfiltered", sim_tree_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX);
TH2D *h2_Hlt1_flags_B_Mass = new TH2D("h2_Hlt1_flags_B_Mass", "Hlt1 Decision vs B Mass", 50, 5100, 5400, 13, 1., 14.); TH2D *h2_Hlt1_flags_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_excl_B_Mass = new TH2D("h2_Hlt1_flags_excl_B_Mass", "Excl Hlt1 Decision vs B Mass", 50, 5100, 5400, 13, 1., 14.);
TH1D *h1_bdt_probs = new TH1D("h1_bdt_probs", "BDT Probabilities", 100, -2, 2); TH1D *h1_bdt_probs = new TH1D("h1_bdt_probs", "BDT Probabilities", 100, -2, 2);
TH1D *h1_Hp_PID_K_pref = new TH1D("h1_Hp_PID_K_pref", "H^{+} PID K, Before TrueID", 50, -10, 10);
TH1D *h1_Hp_PID_K_postf = new TH1D("h1_Hp_PID_K_postf", "H^{+} PID K, After TrueID", 50, -10, 10);
h1_B_Mass_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_bdtf->GetXaxis()->SetTitle(end_state_mass_literal); h1_B_Mass_bdtf->GetXaxis()->SetTitle(end_state_mass_literal);
h2_Hlt1_flags_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); h2_Hlt1_flags_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal);
h2_Hlt1_flags_excl_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); h2_Hlt1_flags_excl_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal);
@ -112,14 +100,12 @@ int new_analysis_bu2hpmumu()
TV::Float("B_BPVFDCHI2", "B_BPVFDCHI2"), TV::Float("B_BPVFDCHI2", "B_BPVFDCHI2"),
TV::Float("B_BPVDIRA", "B_BPVDIRA"), TV::Float("B_BPVDIRA", "B_BPVDIRA"),
TV::Float("Jpsi_BPVIPCHI2", "Jpsi_BPVIPCHI2"), TV::Float("Jpsi_BPVIPCHI2", "Jpsi_BPVIPCHI2"),
// TV::Float("Jpsi_BPVDIRA", "Jpsi_BPVDIRA"),
TV::Float("Jpsi_PT", "Jpsi_PT"), TV::Float("Jpsi_PT", "Jpsi_PT"),
TV::Float("Hp_BPVIPCHI2", "Hp_BPVIPCHI2"),
TV::Float("Hp_PT", "Hp_PT"),
// TV::Double("Kplus_PID_K", "K_PID_K"),
TV::Double("Hp_PROBNN_K", "Hp_PROBNN_K"),
TV::Float("L1_BPVIPCHI2", "L1_BPVIPCHI2"),
TV::Float("L2_BPVIPCHI2", "L2_BPVIPCHI2"),
TV::Float("Kplus_BPVIPCHI2", "Kplus_BPVIPCHI2"),
TV::Float("Kplus_PT", "Kplus_PT"),
TV::Double("Kplus_PROBNN_K", "Kplus_PROBNN_K"),
TV::Float("muminus_BPVIPCHI2", "muminus_BPVIPCHI2"),
TV::Float("muplus_BPVIPCHI2", "muplus_BPVIPCHI2"),
}; };
TTree *sig_tree = new TTree("TreeS", "tree containing signal data"); TTree *sig_tree = new TTree("TreeS", "tree containing signal data");
@ -141,12 +127,13 @@ int new_analysis_bu2hpmumu()
for (unsigned int i = 0; i < data_entries; i++) for (unsigned int i = 0; i < data_entries; i++)
{ {
data_chain->GetEntry(i); data_chain->GetEntry(i);
Double_t reconstructed_B_Mass = (hp4v_data->LorentzVector(K_MASS) + l14v_data->LorentzVector() + l24v_data->LorentzVector()).M();
TLorentzVector dimuon = l14v_data->LorentzVector() + l24v_data->LorentzVector();
Double_t reconstructed_B_Mass = (hp4v_data->LorentzVector() + dimuon).M();
if (std::all_of(vars.begin(), vars.end(), [](TV *v) if (std::all_of(vars.begin(), vars.end(), [](TV *v)
{ return v->IsDataFinite(); })) { return v->IsDataFinite(); }))
{ {
if (reconstructed_B_Mass > 5500.)
if (reconstructed_B_Mass > 5500. && ((TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) || (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.)))
{ {
bkg_tree->Fill(); bkg_tree->Fill();
bkg_events++; bkg_events++;
@ -166,48 +153,43 @@ int new_analysis_bu2hpmumu()
for (unsigned int i = 0; i < sim_entries; i++) for (unsigned int i = 0; i < sim_entries; i++)
{ {
sim_chain->GetEntry(i); sim_chain->GetEntry(i);
Double_t reconstructed_B_Mass = (hp4v_sim->LorentzVector(K_MASS) + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M();
h1_Hp_PID_K_pref->Fill(Hp_PID_K_sim);
h1_B_Mass_sim_unf->Fill(reconstructed_B_Mass);
if (B_BKGCAT == 30 && TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID && TMath::Abs(Hp_TRUEID) == PID_KAON)
Double_t reconstructed_B_Mass = (hp4v_sim->LorentzVector() + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M();
if (sig_events < bkg_events)
{ {
h1_Hp_PID_K_postf->Fill(Hp_PID_K_sim);
if (sig_events < bkg_events)
if (retrain_bdt && std::all_of(vars.begin(), vars.end(), [](TV *v)
{ return v->IsMCFinite(); }))
{ {
if (retrain_bdt && std::all_of(vars.begin(), vars.end(), [](TV *v)
{ return v->IsMCFinite(); }))
{
sig_tree->Fill();
sig_events++;
}
sig_tree->Fill();
sig_events++;
} }
B_Mass_sim_var = reconstructed_B_Mass;
tree_B_Mass_sim->Fill();
} }
B_Mass_sim_var = reconstructed_B_Mass;
tree_B_Mass_sim->Fill();
PrintProgress(TString::Format("%s SIG Collection", analysis_name), sim_entries, 10000, i); PrintProgress(TString::Format("%s SIG Collection", analysis_name), sim_entries, 10000, i);
} }
if (retrain_bdt) if (retrain_bdt)
{ {
TrainBDT(vars, "BuToHpMuMu", sig_tree, bkg_tree);
std::cout << "# Added " << sig_events << " signal events." << std::endl;
TrainBDT(vars, analysis_name, sig_tree, bkg_tree);
std::cout << "# Finished BDT retrain." << std::endl; std::cout << "# Finished BDT retrain." << std::endl;
} }
std::cout << "# Starting evaluation of data." << std::endl; std::cout << "# Starting evaluation of data." << std::endl;
Float_t *train_vars = new Float_t[vars.size()]; Float_t *train_vars = new Float_t[vars.size()];
auto reader = SetupReader(vars, train_vars, "BuToHpMuMu");
auto reader = SetupReader(vars, train_vars, analysis_name);
const double mva_cut_value = -0.02;
const double mva_cut_value = -0.05;
for (unsigned int i = 0; i < data_entries; i++) for (unsigned int i = 0; i < data_entries; i++)
{ {
data_chain->GetEntry(i); data_chain->GetEntry(i);
TLorentzVector dimuon = l14v_data->LorentzVector() + l24v_data->LorentzVector(); TLorentzVector dimuon = l14v_data->LorentzVector() + l24v_data->LorentzVector();
Double_t reconstructed_B_Mass = (hp4v_data->LorentzVector(K_MASS) + dimuon).M();
Double_t reconstructed_B_Mass = (hp4v_data->LorentzVector() + dimuon).M();
if (std::all_of(vars.begin(), vars.end(), [](TV *v) if (std::all_of(vars.begin(), vars.end(), [](TV *v)
{ return v->IsDataFinite(); })) { return v->IsDataFinite(); }))
@ -264,15 +246,13 @@ int new_analysis_bu2hpmumu()
DrawInDefaultCanvas(h2_Hlt1_flags_excl_B_Mass, analysis_name, 0.16, "COLZ"); DrawInDefaultCanvas(h2_Hlt1_flags_excl_B_Mass, analysis_name, 0.16, "COLZ");
DrawInDefaultCanvas(h1_B_Mass_unf, analysis_name, 0.1); DrawInDefaultCanvas(h1_B_Mass_unf, analysis_name, 0.1);
DrawInDefaultCanvas(h1_B_Mass_sim_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_Hp_PID_K_pref, h1_Hp_PID_K_postf}, {kRed, kBlue}, {0, 3003}, analysis_name);
DrawInDefaultCanvasStacked({h1_B_Mass_unf, h1_B_Mass_bdtf}, {kRed, kBlue}, {0, 3003}, analysis_name); DrawInDefaultCanvasStacked({h1_B_Mass_unf, h1_B_Mass_bdtf}, {kRed, kBlue}, {0, 3003}, 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);
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);
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);
@ -282,5 +262,22 @@ int new_analysis_bu2hpmumu()
DrawBDTProbs(h1_bdt_probs, mva_cut_value, analysis_name); DrawBDTProbs(h1_bdt_probs, mva_cut_value, analysis_name);
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);
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;
res_file.close();
return 0; return 0;
} }

21
new_analysis_bu2kpmumu.cpp
View File

@ -3,7 +3,9 @@
#include <cmath> #include <cmath>
#include <algorithm> #include <algorithm>
#include <filesystem> #include <filesystem>
#include <string_view>
#include <string_view
#include <ctime>
#include <fstream>
#include "TH1D.h" #include "TH1D.h"
#include "TH2D.h" #include "TH2D.h"
@ -48,7 +50,6 @@ int new_analysis_bu2kpmumu()
const char *data_tree_name = "Hlt2RD_BuToKpMuMu"; const char *data_tree_name = "Hlt2RD_BuToKpMuMu";
const char *sim_tree_name = "BuToKpMuMu_noPID"; const char *sim_tree_name = "BuToKpMuMu_noPID";
const char *end_state_mass_literal = "m(K^{+}#mu^{+}#mu^{-})"; const char *end_state_mass_literal = "m(K^{+}#mu^{+}#mu^{-})";
const bool retrain_bdt = false;
TChain *data_chain = new TChain(TString::Format("%s/DecayTree", data_tree_name)); TChain *data_chain = new TChain(TString::Format("%s/DecayTree", data_tree_name));
data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_RD.root"); data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_SprucingPass23r1_94000000_RD.root");
@ -165,6 +166,22 @@ int new_analysis_bu2kpmumu()
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);
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;
res_file.close();
return 0; return 0;
} }
Write Preview
Loading…
Cancel
Save