|
@ -50,7 +50,8 @@ int new_analysis_bu2hpmumu() |
|
|
const char *data_tree_name = "SpruceRD_BuToHpMuMu"; |
|
|
const char *data_tree_name = "SpruceRD_BuToHpMuMu"; |
|
|
const char *sim_tree_name = "BuToHpMuMu_noPID_mapped"; |
|
|
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 = true; |
|
|
|
|
|
const bool skip_fit = 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/Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_Sprucing23r1_90000000_RD.root"); |
|
@ -60,9 +61,9 @@ int new_analysis_bu2hpmumu() |
|
|
FourVect *l24v_data = FourVect::Init(data_chain, "muplus"); |
|
|
FourVect *l24v_data = FourVect::Init(data_chain, "muplus"); |
|
|
FourVect *hp4v_data = FourVect::Init(data_chain, "Kplus"); |
|
|
FourVect *hp4v_data = FourVect::Init(data_chain, "Kplus"); |
|
|
|
|
|
|
|
|
// Double_t Kplus_PID_K;
|
|
|
|
|
|
|
|
|
Double_t Kplus_PID_K; |
|
|
|
|
|
|
|
|
// data_chain->SetBranchAddress("Kplus_PID_K", &Kplus_PID_K);
|
|
|
|
|
|
|
|
|
data_chain->SetBranchAddress("Kplus_PID_K", &Kplus_PID_K); |
|
|
|
|
|
|
|
|
TChain *sim_chain = new TChain(TString::Format("%s/DecayTree", sim_tree_name)); |
|
|
TChain *sim_chain = new TChain(TString::Format("%s/DecayTree", sim_tree_name)); |
|
|
sim_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/BuToHpMuMu_mapped_mc.root"); |
|
|
sim_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/MC/BuToHpMuMu_mapped_mc.root"); |
|
@ -100,17 +101,19 @@ int new_analysis_bu2hpmumu() |
|
|
auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name); |
|
|
auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name); |
|
|
std::map<std::string, int> exclusive_hits{}; |
|
|
std::map<std::string, int> exclusive_hits{}; |
|
|
|
|
|
|
|
|
TV* kplus_pid_k_var = TV::Double("Kplus_PID_K", "Kplus_PID_K"); |
|
|
|
|
|
|
|
|
// TV* kplus_pid_k_var = TV::Double("Kplus_PID_K", "Kplus_PID_K");
|
|
|
|
|
|
|
|
|
std::vector<TV *> vars{ |
|
|
std::vector<TV *> vars{ |
|
|
TV::Float("B_PT", "B_PT"), |
|
|
TV::Float("B_PT", "B_PT"), |
|
|
TV::Float("B_BPVFDCHI2", "B_BPVFDCHI2"), |
|
|
TV::Float("B_BPVFDCHI2", "B_BPVFDCHI2"), |
|
|
TV::Float("B_BPVDIRA", "B_BPVDIRA"), |
|
|
TV::Float("B_BPVDIRA", "B_BPVDIRA"), |
|
|
|
|
|
TV::Double("B_CHI2", "B_CHI2"), |
|
|
TV::Float("Jpsi_BPVIPCHI2", "Jpsi_BPVIPCHI2"), |
|
|
TV::Float("Jpsi_BPVIPCHI2", "Jpsi_BPVIPCHI2"), |
|
|
TV::Float("Jpsi_PT", "Jpsi_PT"), |
|
|
TV::Float("Jpsi_PT", "Jpsi_PT"), |
|
|
|
|
|
TV::Double("Jpsi_CHI2", "Jpsi_CHI2"), |
|
|
TV::Float("Kplus_BPVIPCHI2", "Kplus_BPVIPCHI2"), |
|
|
TV::Float("Kplus_BPVIPCHI2", "Kplus_BPVIPCHI2"), |
|
|
TV::Float("Kplus_PT", "Kplus_PT"), |
|
|
TV::Float("Kplus_PT", "Kplus_PT"), |
|
|
kplus_pid_k_var, |
|
|
|
|
|
|
|
|
// kplus_pid_k_var,
|
|
|
TV::Float("muminus_BPVIPCHI2", "muminus_BPVIPCHI2"), |
|
|
TV::Float("muminus_BPVIPCHI2", "muminus_BPVIPCHI2"), |
|
|
// TV::Float("muminus_PT", "muminus_PT"),
|
|
|
// TV::Float("muminus_PT", "muminus_PT"),
|
|
|
TV::Float("muplus_BPVIPCHI2", "muplus_BPVIPCHI2"), |
|
|
TV::Float("muplus_BPVIPCHI2", "muplus_BPVIPCHI2"), |
|
@ -142,7 +145,7 @@ int new_analysis_bu2hpmumu() |
|
|
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. && ((TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) || (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.)) && Kplus_PID_K > -3) |
|
|
{ |
|
|
{ |
|
|
bkg_tree->Fill(); |
|
|
bkg_tree->Fill(); |
|
|
bkg_events++; |
|
|
bkg_events++; |
|
@ -186,6 +189,12 @@ int new_analysis_bu2hpmumu() |
|
|
std::cout << "# Finished BDT retrain." << std::endl; |
|
|
std::cout << "# Finished BDT retrain." << std::endl; |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
if (skip_fit) |
|
|
|
|
|
{ |
|
|
|
|
|
std::cout << "# Skipping evaluation of data." << std::endl; |
|
|
|
|
|
return 0; |
|
|
|
|
|
} |
|
|
|
|
|
|
|
|
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()]; |
|
@ -222,23 +231,26 @@ int new_analysis_bu2hpmumu() |
|
|
FillHlt1DecisionHistos(hlt1_decision_histos, reconstructed_B_Mass); |
|
|
FillHlt1DecisionHistos(hlt1_decision_histos, reconstructed_B_Mass); |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
double mva_response = reader->EvaluateMVA("BDT"); |
|
|
|
|
|
h1_bdt_probs->Fill(mva_response); |
|
|
|
|
|
|
|
|
if (Kplus_PID_K > -3 && ((TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) || (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.))) |
|
|
|
|
|
{ |
|
|
|
|
|
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 && kplus_pid_k_var->GetDataDouble() > -3) |
|
|
|
|
|
{ |
|
|
|
|
|
h1_B_Mass_bdtf->Fill(reconstructed_B_Mass); |
|
|
|
|
|
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 (mva_response > mva_cut_value) |
|
|
{ |
|
|
{ |
|
|
B_Mass_psi2s_var = reconstructed_B_Mass; |
|
|
|
|
|
tree_B_Mass_psi2s->Fill(); |
|
|
|
|
|
|
|
|
h1_B_Mass_bdtf->Fill(reconstructed_B_Mass); |
|
|
|
|
|
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(); |
|
|
|
|
|
} |
|
|
} |
|
|
} |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
@ -290,36 +302,41 @@ int new_analysis_bu2hpmumu() |
|
|
res_file << "#### " << analysis_name << " @ " << std::put_time(&tm, "%c") << " ####" << std::endl; |
|
|
res_file << "#### " << analysis_name << " @ " << std::put_time(&tm, "%c") << " ####" << std::endl; |
|
|
res_file << "J/Psi Mode: " << ErrToStr(roofit_hist_jpsi_fitsum.signal_yield, 0) << " / " << ErrToStr(roofit_hist_jpsi_fitsum.background_yield, 0) << std::endl; |
|
|
res_file << "J/Psi Mode: " << ErrToStr(roofit_hist_jpsi_fitsum.signal_yield, 0) << " / " << ErrToStr(roofit_hist_jpsi_fitsum.background_yield, 0) << std::endl; |
|
|
res_file << " Sig/Bkg: " << ErrToStr(jpsi_sigobkg, 2) << std::endl; |
|
|
res_file << " Sig/Bkg: " << ErrToStr(jpsi_sigobkg, 2) << std::endl; |
|
|
res_file << "Psi(2S) Mode: " << ErrToStr(roofit_hist_psi2s_fitsum.signal_yield, 0) << " / " << ErrToStr(roofit_hist_psi2s_fitsum.background_yield, 0) << std::endl; |
|
|
|
|
|
|
|
|
res_file << "Psi(2S) Mode: " << ErrToStr(roofit_hist_psi2s_fitsum.signal_yield, 0) << " / " << ErrToStr(roofit_hist_psi2s_fitsum.background_yield, 0) << std::endl; |
|
|
res_file << " Sig/Bkg: " << ErrToStr(psi2s_sigobkg, 2) << std::endl; |
|
|
res_file << " Sig/Bkg: " << ErrToStr(psi2s_sigobkg, 2) << std::endl; |
|
|
res_file << "Mode Yield Ratio: " << ErrToStr(signal_ratio, 3) << std::endl; |
|
|
res_file << "Mode Yield Ratio: " << ErrToStr(signal_ratio, 3) << std::endl; |
|
|
res_file << "Rel Br Frac MuMu: " << ErrToStr(mumu_br_frac, 3) << std::endl; |
|
|
res_file << "Rel Br Frac MuMu: " << ErrToStr(mumu_br_frac, 3) << std::endl; |
|
|
|
|
|
|
|
|
auto br_frac = MultWithErr(signal_ratio.first, signal_ratio.second, mumu_br_frac.first, mumu_br_frac.second); |
|
|
auto br_frac = MultWithErr(signal_ratio.first, signal_ratio.second, mumu_br_frac.first, mumu_br_frac.second); |
|
|
res_file << "Rel Br Frac: " << ErrToStr(br_frac, 3) << std::endl << std::endl; |
|
|
|
|
|
|
|
|
res_file << "Rel Br Frac: " << ErrToStr(br_frac, 3) << std::endl |
|
|
|
|
|
<< std::endl; |
|
|
|
|
|
|
|
|
res_file << "Params from Sim:" << std::endl << roofit_hist_sim.shape_parameters.ToString() << std::endl; |
|
|
|
|
|
|
|
|
res_file << "Params from Sim:" << std::endl |
|
|
|
|
|
<< roofit_hist_sim.shape_parameters.ToString() << std::endl; |
|
|
|
|
|
|
|
|
for (const auto &par : roofit_hist_sim.fitted_params) |
|
|
for (const auto &par : roofit_hist_sim.fitted_params) |
|
|
{ |
|
|
{ |
|
|
res_file << par.ToString(true).c_str() << std::endl; |
|
|
res_file << par.ToString(true).c_str() << std::endl; |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
res_file << std::endl << "Fitted Parameters: J/PSI" << std::endl << std::endl; |
|
|
|
|
|
|
|
|
res_file << std::endl |
|
|
|
|
|
<< "Fitted Parameters: J/PSI" << std::endl |
|
|
|
|
|
<< std::endl; |
|
|
|
|
|
|
|
|
for (const auto &par : roofit_hist_jpsi_fitsum.fitted_params) |
|
|
for (const auto &par : roofit_hist_jpsi_fitsum.fitted_params) |
|
|
{ |
|
|
{ |
|
|
res_file << par.ToString(true).c_str() << std::endl; |
|
|
res_file << par.ToString(true).c_str() << std::endl; |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
res_file << std::endl << "Fitted Parameters: PSI(2S)" << std::endl << std::endl; |
|
|
|
|
|
|
|
|
res_file << std::endl |
|
|
|
|
|
<< "Fitted Parameters: PSI(2S)" << std::endl |
|
|
|
|
|
<< std::endl; |
|
|
|
|
|
|
|
|
for (const auto &par : roofit_hist_psi2s_fitsum.fitted_params) |
|
|
for (const auto &par : roofit_hist_psi2s_fitsum.fitted_params) |
|
|
{ |
|
|
{ |
|
|
res_file << par.ToString(true).c_str() << std::endl; |
|
|
res_file << par.ToString(true).c_str() << std::endl; |
|
|
} |
|
|
} |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
auto print_table = [&res_file](std::string name, std::pair<double, double> sig, std::pair<double, double> bkg) |
|
|
auto print_table = [&res_file](std::string name, std::pair<double, double> sig, std::pair<double, double> bkg) |
|
|
{ |
|
|
{ |
|
|
res_file << std::endl; |
|
|
res_file << std::endl; |
|
|