pfeiffer/inclusive_detached_dilepton
1
0
Fork 0
Code Issues Pull Requests 1 repo.project_board Releases Wiki Activity

Compare commits

base: pfeiffer:9a73561ed4f0cecdf5c86a27f9f6f2ef1d24489b
Branches Tags
pfeiffer:master
pfeiffer:training_tryouts
..
compare: pfeiffer:84df7596bf1312ecca658b3b9a96c35db964f888
Branches Tags
pfeiffer:training_tryouts
pfeiffer:master

No commits in common. "9a73561ed4f0cecdf5c86a27f9f6f2ef1d24489b" and "84df7596bf1312ecca658b3b9a96c35db964f888" have entirely different histories.
9a73561ed4 ... 84df7596bf

9 changed files with 256 additions and 269 deletions
Show Stats Expand all files Collapse all files

6
B0ToHpHmMuMu_results.txt
View File

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

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

69
basic_analysis.h
View File

@ -15,12 +15,6 @@
#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:
@ -104,8 +98,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) {
{
this->title = title; this->title = title;
this->name = name; this->name = name;
this->value = value; this->value = value;
@ -134,7 +127,6 @@ 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
@ -142,8 +134,6 @@ 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;
}; };
@ -166,19 +156,17 @@ 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]);
} }
@ -384,7 +372,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, bool const_sigma = false, 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, 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)
{ {
@ -399,7 +387,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())));
@ -416,38 +404,20 @@ RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString
RooRealVar var_mass_nR(suffix_name("var_mass_nR"), "n_{R}", 5., 0., 15.); RooRealVar var_mass_nR(suffix_name("var_mass_nR"), "n_{R}", 5., 0., 15.);
if (useExtShape) if (useExtShape)
{
if (extShape.alpha_left != 0.)
{ {
var_mass_alphaL.setConstant(true); var_mass_alphaL.setConstant(true);
var_mass_alphaL.setVal(extShape.alpha_left); var_mass_alphaL.setVal(extShape.alpha_left);
}
if (extShape.n_left != 0.)
{
var_mass_nL.setConstant(true); var_mass_nL.setConstant(true);
var_mass_nL.setVal(extShape.n_left); var_mass_nL.setVal(extShape.n_left);
}
if (extShape.alpha_right != 0.)
{
var_mass_alphaR.setConstant(true); var_mass_alphaR.setConstant(true);
var_mass_alphaR.setVal(extShape.alpha_right); var_mass_alphaR.setVal(extShape.alpha_right);
}
if (extShape.n_right != 0.)
{
var_mass_nR.setConstant(true); var_mass_nR.setConstant(true);
var_mass_nR.setVal(extShape.n_right); 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");
RooCrystalBall sig_cb(signal_name, "CB Signal", roo_var_mass, var_mass_x0, var_mass_sigmaLR, var_mass_alphaL, var_mass_nL, var_mass_alphaR, var_mass_nR); RooCrystalBall sig_cb(signal_name, "CB Signal", roo_var_mass, var_mass_x0, var_mass_sigmaLR, var_mass_alphaL, var_mass_nL, var_mass_alphaR, var_mass_nR);
@ -457,10 +427,6 @@ 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
@ -487,20 +453,18 @@ 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 significance_val = sig_val / TMath::Sqrt(sig_val + bkg_val); double sig_val = var_mass_nsig.getVal();
double err_prop_sig = (sig_val + 2 * bkg_val) / (2 * TMath::Power((sig_val + bkg_val), (3 / 2))); double sig_err = var_mass_nsig.getError();
double err_prop_bkg = -sig_val / (2 * TMath::Power((sig_val + bkg_val), (3 / 2))); double bkg_val = var_mass_nbkg.getVal();
double significance_err = TMath::Sqrt(TMath::Sq(err_prop_sig * sig_err) + TMath::Sq(err_prop_bkg * bkg_err)); double bkg_err = var_mass_nbkg.getError();
auto sig_over_bkg = DivWithErr(sig_val, sig_err, bkg_val, 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));
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}/#sqrt{N_{Sig} + N_{Bkg}}", sig_over_bkg_val, sig_over_bkg_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));
@ -534,15 +498,12 @@ 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,10 +16,4 @@ 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/F"); output_tree->Branch("B0_PT", &B0_PT_out, "B0_PT/D");
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");

145
new_analysis_b02hphmmumu.cpp
View File

@ -4,8 +4,6 @@
#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"
@ -47,26 +45,39 @@
int new_analysis_b02hphmmumu() int new_analysis_b02hphmmumu()
{ {
const char *analysis_name = "B0ToHpHmMuMu"; const char *analysis_name = "B0ToHpHmMuMu";
const char *data_tree_name = "SpruceRD_B0ToHpHmMuMu"; const char *data_tree_name = "B0ToHpHmMuMu";
const char *sim_tree_name = "B0ToHpHmMuMu_noPID_mapped"; const char *sim_tree_name = "B0ToHpHmMuMu_noPID";
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/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_Sprucing23r1_90000000_RD.root"); // 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");
FourVect *l14v_data = FourVect::Init(data_chain, "muminus"); Double_t Hp_PID_K, Hm_PID_K;
FourVect *l24v_data = FourVect::Init(data_chain, "muplus"); data_chain->SetBranchAddress("Hp_PID_K", &Hp_PID_K);
FourVect *hp4v_data = FourVect::Init(data_chain, "Kplus"); data_chain->SetBranchAddress("Hm_PID_K", &Hm_PID_K);
FourVect *hm4v_data = FourVect::Init(data_chain, "piminus");
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");
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/rd_btoxll_simulation_fullstream_v0r0p6671378_B0ToKpPimMuMu_11144002_magdown.root");
FourVect *l14v_sim = FourVect::Init(sim_chain, "muminus"); FourVect *l14v_sim = FourVect::Init(sim_chain, "L1");
FourVect *l24v_sim = FourVect::Init(sim_chain, "muplus"); FourVect *l24v_sim = FourVect::Init(sim_chain, "L2");
FourVect *hp4v_sim = FourVect::Init(sim_chain, "Kplus"); FourVect *hp4v_sim = FourVect::Init(sim_chain, "Hp");
FourVect *hm4v_sim = FourVect::Init(sim_chain, "piminus"); FourVect *hm4v_sim = FourVect::Init(sim_chain, "Hm");
Int_t B_BKGCAT, L1_TRUEID, L2_TRUEID, Hp_TRUEID, Hm_TRUEID;
sim_chain->SetBranchAddress("L1_TRUEID", &L1_TRUEID);
sim_chain->SetBranchAddress("L2_TRUEID", &L2_TRUEID);
sim_chain->SetBranchAddress("Hp_TRUEID", &Hp_TRUEID);
sim_chain->SetBranchAddress("Hm_TRUEID", &Hm_TRUEID);
sim_chain->SetBranchAddress("B0_BKGCAT", &B_BKGCAT);
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";
@ -83,12 +94,14 @@ 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);
@ -102,14 +115,17 @@ 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("Kplus_BPVIPCHI2", "Kplus_BPVIPCHI2"), TV::Float("Hp_BPVIPCHI2", "Hp_BPVIPCHI2"),
TV::Float("Kplus_PT", "Kplus_PT"), TV::Float("Hp_PT", "Hp_PT"),
TV::Float("piminus_BPVIPCHI2", "piminus_BPVIPCHI2"), TV::Float("Hm_BPVIPCHI2", "Hm_BPVIPCHI2"),
TV::Float("piminus_PT", "piminus_PT"), TV::Float("Hm_PT", "Hm_PT"),
TV::Double("Kplus_PROBNN_K", "Kplus_PROBNN_K"), // TV::Double("Kplus_PID_K", "K_PID_K"),
TV::Float("muminus_BPVIPCHI2", "muminus_BPVIPCHI2"), TV::Double("Hp_PROBNN_K", "Hp_PROBNN_K"),
TV::Float("muplus_BPVIPCHI2", "muplus_BPVIPCHI2"), TV::Double("Hm_PROBNN_K", "Hm_PROBNN_K"),
TV::Float("L1_BPVIPCHI2", "L1_BPVIPCHI2"),
TV::Float("L2_BPVIPCHI2", "L2_BPVIPCHI2"),
}; };
TTree *sig_tree = new TTree("TreeS", "tree containing signal data"); TTree *sig_tree = new TTree("TreeS", "tree containing signal data");
@ -129,19 +145,33 @@ 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 = hp4v_data->LorentzVector() + hm4v_data->LorentzVector(); TLorentzVector reconstructed_Kstar{};
TLorentzVector dimuon = l14v_data->LorentzVector() + l24v_data->LorentzVector(); bool found_k_star = false;
Double_t reconstructed_B_Mass = (reconstructed_Kstar + dimuon).M(); 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;
}
if (found_k_star)
{
Double_t reconstructed_B_Mass = (reconstructed_Kstar + l14v_data->LorentzVector() + l24v_data->LorentzVector()).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. && ((TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) || (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.))) if (reconstructed_B_Mass > 5500.)
{ {
bkg_tree->Fill(); bkg_tree->Fill();
bkg_events++; bkg_events++;
} }
} }
}
PrintProgress(TString::Format("%s BKG Collection", analysis_name), data_entries, 10000, i); PrintProgress(TString::Format("%s BKG Collection", analysis_name), data_entries, 10000, i);
} }
@ -157,7 +187,28 @@ int new_analysis_b02hphmmumu()
{ {
sim_chain->GetEntry(i); sim_chain->GetEntry(i);
TLorentzVector reconstructed_Kstar = hp4v_sim->LorentzVector() + hm4v_sim->LorentzVector(); Double_t reco_mass_pipkp = (hp4v_sim->LorentzVector(K_MASS) + hm4v_sim->LorentzVector() + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M();
Double_t reco_mass_pimkm = (hp4v_sim->LorentzVector() + hm4v_sim->LorentzVector(K_MASS) + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M();
h1_B_Mass_sim_unf->Fill(reco_mass_pipkp);
h1_B_Mass_sim_unf->Fill(reco_mass_pimkm);
if (B_BKGCAT == 30 && TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID)
{
TLorentzVector reconstructed_Kstar{};
bool found_k_star = false;
if (TMath::Abs(Hp_TRUEID) == PID_KAON && TMath::Abs(Hm_TRUEID) == PID_PION)
{
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(); Double_t reconstructed_B_Mass = (reconstructed_Kstar + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M();
if (sig_events < bkg_events) if (sig_events < bkg_events)
@ -172,13 +223,14 @@ int new_analysis_b02hphmmumu()
B_Mass_sim_var = reconstructed_B_Mass; B_Mass_sim_var = reconstructed_B_Mass;
tree_B_Mass_sim->Fill(); 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;
} }
@ -188,14 +240,28 @@ 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.0508;
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 = hp4v_data->LorentzVector() + hm4v_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;
}
if (found_k_star)
{
Double_t reconstructed_B_Mass = (reconstructed_Kstar + dimuon).M(); Double_t reconstructed_B_Mass = (reconstructed_Kstar + dimuon).M();
if (std::all_of(vars.begin(), vars.end(), [](TV *v) if (std::all_of(vars.begin(), vars.end(), [](TV *v)
@ -242,6 +308,7 @@ int new_analysis_b02hphmmumu()
} }
} }
} }
}
PrintProgress(TString::Format("%s BDT Evaluation", analysis_name), data_entries, 10000, i); PrintProgress(TString::Format("%s BDT Evaluation", analysis_name), data_entries, 10000, i);
} }
@ -256,13 +323,14 @@ 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, 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);
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);
@ -272,22 +340,5 @@ 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,8 +4,6 @@
#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"
@ -50,6 +48,7 @@ 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");
@ -126,6 +125,11 @@ 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);
@ -145,6 +149,24 @@ 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();
} }
@ -182,22 +204,7 @@ int new_analysis_b02kppimmumu()
// 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::cout << "hist entries: " << tree_B_Mass_jpsi->GetEntries() << ", kplus: " << kplus << ", kminus: " << kminus << ", piplus: " << piplus << ", piminus: " << piminus << std::endl;
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;
} }

91
new_analysis_bu2hpmumu.cpp
View File

@ -4,8 +4,6 @@
#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"
@ -47,24 +45,34 @@
int new_analysis_bu2hpmumu() int new_analysis_bu2hpmumu()
{ {
const char *analysis_name = "BuToHpMuMu"; const char *analysis_name = "BuToHpMuMu";
const char *data_tree_name = "SpruceRD_BuToHpMuMu"; const char *data_tree_name = "BuToHpMuMu";
const char *sim_tree_name = "BuToHpMuMu_noPID_mapped"; const char *sim_tree_name = "BuToHpMuMu_noPID";
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/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_Sprucing23r1_90000000_RD.root"); // 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");
FourVect *l14v_data = FourVect::Init(data_chain, "muminus"); FourVect *l14v_data = FourVect::Init(data_chain, "L1");
FourVect *l24v_data = FourVect::Init(data_chain, "muplus"); FourVect *l24v_data = FourVect::Init(data_chain, "L2");
FourVect *hp4v_data = FourVect::Init(data_chain, "Kplus"); FourVect *hp4v_data = FourVect::Init(data_chain, "Hp");
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/rd_btoxll_simulation_fullstream_v0r0p6671378_BuToKpMuMu_12143001_magdown.root");
FourVect *l14v_sim = FourVect::Init(sim_chain, "muminus"); FourVect *l14v_sim = FourVect::Init(sim_chain, "L1");
FourVect *l24v_sim = FourVect::Init(sim_chain, "muplus"); FourVect *l24v_sim = FourVect::Init(sim_chain, "L2");
FourVect *hp4v_sim = FourVect::Init(sim_chain, "Kplus"); FourVect *hp4v_sim = FourVect::Init(sim_chain, "Hp");
Int_t B_BKGCAT, L1_TRUEID, L2_TRUEID, Hp_TRUEID;
Double_t Hp_PID_K_sim;
sim_chain->SetBranchAddress("L1_TRUEID", &L1_TRUEID);
sim_chain->SetBranchAddress("L2_TRUEID", &L2_TRUEID);
sim_chain->SetBranchAddress("Hp_TRUEID", &Hp_TRUEID);
sim_chain->SetBranchAddress("B_BKGCAT", &B_BKGCAT);
sim_chain->SetBranchAddress("Hp_PID_K", &Hp_PID_K_sim);
Double_t B_Mass_jpsi_var, B_Mass_psi2s_var, B_Mass_sim_var; 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";
@ -81,12 +89,16 @@ 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);
@ -100,12 +112,14 @@ 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("Kplus_BPVIPCHI2", "Kplus_BPVIPCHI2"), TV::Float("Hp_BPVIPCHI2", "Hp_BPVIPCHI2"),
TV::Float("Kplus_PT", "Kplus_PT"), TV::Float("Hp_PT", "Hp_PT"),
TV::Double("Kplus_PROBNN_K", "Kplus_PROBNN_K"), // TV::Double("Kplus_PID_K", "K_PID_K"),
TV::Float("muminus_BPVIPCHI2", "muminus_BPVIPCHI2"), TV::Double("Hp_PROBNN_K", "Hp_PROBNN_K"),
TV::Float("muplus_BPVIPCHI2", "muplus_BPVIPCHI2"), TV::Float("L1_BPVIPCHI2", "L1_BPVIPCHI2"),
TV::Float("L2_BPVIPCHI2", "L2_BPVIPCHI2"),
}; };
TTree *sig_tree = new TTree("TreeS", "tree containing signal data"); TTree *sig_tree = new TTree("TreeS", "tree containing signal data");
@ -127,13 +141,12 @@ 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);
TLorentzVector dimuon = l14v_data->LorentzVector() + l24v_data->LorentzVector(); Double_t reconstructed_B_Mass = (hp4v_data->LorentzVector(K_MASS) + l14v_data->LorentzVector() + l24v_data->LorentzVector()).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(); }))
{ {
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_tree->Fill();
bkg_events++; bkg_events++;
@ -153,7 +166,12 @@ 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() + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M(); Double_t reconstructed_B_Mass = (hp4v_sim->LorentzVector(K_MASS) + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M();
h1_Hp_PID_K_pref->Fill(Hp_PID_K_sim);
h1_B_Mass_sim_unf->Fill(reconstructed_B_Mass);
if (B_BKGCAT == 30 && TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID && TMath::Abs(Hp_TRUEID) == PID_KAON)
{
h1_Hp_PID_K_postf->Fill(Hp_PID_K_sim);
if (sig_events < bkg_events) if (sig_events < bkg_events)
{ {
if (retrain_bdt && std::all_of(vars.begin(), vars.end(), [](TV *v) if (retrain_bdt && std::all_of(vars.begin(), vars.end(), [](TV *v)
@ -166,30 +184,30 @@ int new_analysis_bu2hpmumu()
B_Mass_sim_var = reconstructed_B_Mass; B_Mass_sim_var = reconstructed_B_Mass;
tree_B_Mass_sim->Fill(); 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, "BuToHpMuMu", 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;
} }
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, analysis_name); auto reader = SetupReader(vars, train_vars, "BuToHpMuMu");
const double mva_cut_value = -0.05; const double mva_cut_value = -0.02;
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() + dimuon).M(); Double_t reconstructed_B_Mass = (hp4v_data->LorentzVector(K_MASS) + 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(); }))
@ -246,13 +264,15 @@ 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, 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);
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);
@ -262,22 +282,5 @@ 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,9 +3,7 @@
#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"
@ -50,6 +48,7 @@ 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");
@ -166,22 +165,6 @@ 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;
} }