From bcb6249b551c08e651b7a16ffec478957a66904d Mon Sep 17 00:00:00 2001 From: Marius Pfeiffer Date: Tue, 12 Mar 2024 16:56:43 +0100 Subject: [PATCH] sig over bkg, br frac calc --- B0ToHpHmMuMu_results.txt | 6 +++++ BuToHpMuMu_results.txt | 6 +++++ basic_analysis.h | 31 +++++++++++++++++++------ constants.h | 6 +++++ new_analysis_b02hphmmumu.cpp | 19 ++++++++++++++++ new_analysis_b02kppimmumu.cpp | 43 +++++++++++++++-------------------- new_analysis_bu2hpmumu.cpp | 19 ++++++++++++++++ new_analysis_bu2kpmumu.cpp | 21 +++++++++++++++-- 8 files changed, 117 insertions(+), 34 deletions(-) create mode 100644 B0ToHpHmMuMu_results.txt create mode 100644 BuToHpMuMu_results.txt diff --git a/B0ToHpHmMuMu_results.txt b/B0ToHpHmMuMu_results.txt new file mode 100644 index 0000000..990e360 --- /dev/null +++ b/B0ToHpHmMuMu_results.txt @@ -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 diff --git a/BuToHpMuMu_results.txt b/BuToHpMuMu_results.txt new file mode 100644 index 0000000..8db82a9 --- /dev/null +++ b/BuToHpMuMu_results.txt @@ -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 diff --git a/basic_analysis.h b/basic_analysis.h index d135dc0..c75e8e0 100644 --- a/basic_analysis.h +++ b/basic_analysis.h @@ -15,6 +15,12 @@ #include "constants.h" +std::pair 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 { private: @@ -136,6 +142,8 @@ struct RooFitSummary RooPlot *fit_histogram; RooPlot *pull_histogram; std::map pdf_names; + std::pair signal_yield; + std::pair background_yield; std::vector fitted_params; ShapeParamters shape_parameters; }; @@ -449,6 +457,10 @@ RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString pdf_names[signal_name.Data()] = sig_cb.getTitle().Data(); TString pull_compare_name{}; + double sig_val = 0.; + double sig_err = 0.; + double bkg_val = 0.; + double bkg_err = 0.; if (hasExpBkg) { // Exponential for Background @@ -475,17 +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_nbkg, 2)); - 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(); + 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_over_bkg_val = sig_val / TMath::Sqrt(sig_val + bkg_val); + 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 sig_over_bkg_err = TMath::Sqrt(TMath::Sq(err_prop_sig * sig_err) + TMath::Sq(err_prop_bkg * bkg_err)); + double significance_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)); @@ -519,6 +534,8 @@ RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString roo_frame_mass, roo_frame_pull, pdf_names, + std::make_pair(sig_val, sig_err), + std::make_pair(bkg_val, bkg_err), fitted_params, ShapeParamters{ var_mass_alphaL.getVal(), diff --git a/constants.h b/constants.h index 9e481e8..0aa09ed 100644 --- a/constants.h +++ b/constants.h @@ -16,4 +16,10 @@ const int PID_KAON = 321; const int PID_PION = 211; 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 \ No newline at end of file diff --git a/new_analysis_b02hphmmumu.cpp b/new_analysis_b02hphmmumu.cpp index edcb601..6bd3465 100644 --- a/new_analysis_b02hphmmumu.cpp +++ b/new_analysis_b02hphmmumu.cpp @@ -4,6 +4,8 @@ #include #include #include +#include +#include #include "TH1D.h" #include "TH2D.h" @@ -270,5 +272,22 @@ int new_analysis_b02hphmmumu() 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; } \ No newline at end of file diff --git a/new_analysis_b02kppimmumu.cpp b/new_analysis_b02kppimmumu.cpp index b1d607f..8033106 100644 --- a/new_analysis_b02kppimmumu.cpp +++ b/new_analysis_b02kppimmumu.cpp @@ -4,6 +4,8 @@ #include #include #include +#include +#include #include "TH1D.h" #include "TH2D.h" @@ -48,7 +50,6 @@ int new_analysis_b02kppimmumu() const char *data_tree_name = "Hlt2RD_B0ToKpPimMuMu"; const char *sim_tree_name = "B0ToKpPimMuMu_noPID"; 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)); 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; - int kplus = 0; - int kminus = 0; - int piplus = 0; - int piminus = 0; - for (unsigned int i = 0; i < data_entries; i++) { data_chain->GetEntry(i); @@ -149,24 +145,6 @@ int new_analysis_b02kppimmumu() { 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; tree_B_Mass_jpsi->Fill(); } @@ -204,7 +182,22 @@ int new_analysis_b02kppimmumu() // 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; } \ No newline at end of file diff --git a/new_analysis_bu2hpmumu.cpp b/new_analysis_bu2hpmumu.cpp index f098e03..b7cba27 100644 --- a/new_analysis_bu2hpmumu.cpp +++ b/new_analysis_bu2hpmumu.cpp @@ -4,6 +4,8 @@ #include #include #include +#include +#include #include "TH1D.h" #include "TH2D.h" @@ -260,5 +262,22 @@ int new_analysis_bu2hpmumu() 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; } diff --git a/new_analysis_bu2kpmumu.cpp b/new_analysis_bu2kpmumu.cpp index ed87c63..cd58590 100644 --- a/new_analysis_bu2kpmumu.cpp +++ b/new_analysis_bu2kpmumu.cpp @@ -3,7 +3,9 @@ #include #include #include -#include +#include +#include #include "TH1D.h" #include "TH2D.h" @@ -48,7 +50,6 @@ int new_analysis_bu2kpmumu() const char *data_tree_name = "Hlt2RD_BuToKpMuMu"; const char *sim_tree_name = "BuToKpMuMu_noPID"; 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)); 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); // 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; } \ No newline at end of file