sig over bkg, br frac calc

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

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

basic_analysis.h

@@ -15,6 +15,12 @@
 
 #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
 {
 private:
@@ -136,6 +142,8 @@ struct RooFitSummary
   RooPlot *fit_histogram;
   RooPlot *pull_histogram;
   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;
   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();
+    sig_val = var_mass_nsig.getVal();
-    double sig_err = var_mass_nsig.getError();
+    sig_err = var_mass_nsig.getError();
-    double bkg_val = var_mass_nbkg.getVal();
+    bkg_val = var_mass_nbkg.getVal();
-    double bkg_err = var_mass_nbkg.getError();
+    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));
 
-    fitted_params.push_back(FittedParam("sig_over_bkg", "N_{Sig}/#sqrt{N_{Sig} + N_{Bkg}}", sig_over_bkg_val, sig_over_bkg_err, 2));
+    auto sig_over_bkg = DivWithErr(sig_val, sig_err, bkg_val, bkg_err);
+
+    fitted_params.push_back(FittedParam("significance", "N_{Sig}/#sqrt{N_{Sig} + N_{Bkg}}", significance_val, significance_err, 2));
+    fitted_params.push_back(FittedParam("sig_over_bkg", "N_{Sig}/N_{Bkg}", sig_over_bkg.first, sig_over_bkg.second, 2));
 
     fitted_params.push_back(FittedParam(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(),

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

new_analysis_b02hphmmumu.cpp

@@ -4,6 +4,8 @@
 #include <algorithm>
 #include <filesystem>
 #include <string_view>
+#include <ctime>
+#include <fstream>
 
 #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;
 }

new_analysis_b02kppimmumu.cpp

@@ -4,6 +4,8 @@
 #include <algorithm>
 #include <filesystem>
 #include <string_view>
+#include <ctime>
+#include <fstream>
 
 #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;
 }

new_analysis_bu2hpmumu.cpp

@@ -4,6 +4,8 @@
 #include <algorithm>
 #include <filesystem>
 #include <string_view>
+#include <ctime>
+#include <fstream>
 
 #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;
 }

new_analysis_bu2kpmumu.cpp

@@ -3,7 +3,9 @@
 #include <cmath>
 #include <algorithm>
 #include <filesystem>
-#include <string_view>
+#include <string_view
+#include <ctime>
+#include <fstream>
 
 #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;
 }