From 5153c4c4d68f08c88199e5cc1c1139fa35ac74c2 Mon Sep 17 00:00:00 2001 From: Marius Pfeiffer Date: Sun, 10 Mar 2024 20:24:41 +0100 Subject: [PATCH 1/2] alter ana macros to use pre-subbed root files, formatting, constant sigma --- basic_analysis.h | 58 +++++++--- mapmc_b02hphmmumu.cpp | 2 +- new_analysis_b02hphmmumu.cpp | 216 ++++++++++++----------------------- new_analysis_bu2hpmumu.cpp | 88 ++++++-------- 4 files changed, 147 insertions(+), 217 deletions(-) diff --git a/basic_analysis.h b/basic_analysis.h index 6aa0dae..d135dc0 100644 --- a/basic_analysis.h +++ b/basic_analysis.h @@ -98,7 +98,8 @@ struct FittedParam 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->name = name; this->value = value; @@ -127,6 +128,7 @@ struct ShapeParamters double n_left; double alpha_right; double n_right; + double sigma_lr; }; struct RooFitSummary @@ -156,17 +158,19 @@ void DrawInDefaultCanvasStacked(std::vector histograms, std::vectorGetName()); TCanvas *c = new TCanvas(name, histograms[0]->GetName(), 0, 0, 800, 600); c->SetLeftMargin(margin_left); + 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++) { 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->SetLineColor(colors[i]); hist_clone->SetMaximum(scaling_factor + (scaling_factor * 0.05)); hist_clone->SetMinimum(0.); hist_clone->SetStats(0); - if (fill_style[i] != 0) { + if (fill_style[i] != 0) + { hist_clone->SetFillStyle(fill_style[i]); hist_clone->SetFillColor(colors[i]); } @@ -372,7 +376,7 @@ RooPlot *CreateRooFitHistogram(TH1D *hist) 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) { @@ -381,13 +385,13 @@ RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString Double_t fitRangeUp = fit_low == 0 ? B_MASS_VAR_MIN : fit_low; Double_t fitRangeLow = fit_up == 0 ? B_MASS_VAR_MAX : fit_up; - + RooRealVar roo_var_mass(var_name, "B Mass Variable", B_MASS_VAR_MIN, B_MASS_VAR_MAX); const char *fitting_range_name = "fitting_range"; roo_var_mass.setRange(fitting_range_name, fitRangeUp, fitRangeLow); 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)); RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title(dataSet->GetTitle()), RooFit::Name(TString::Format("%s_rplt", dataSet->GetName()))); @@ -405,17 +409,35 @@ RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString 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); + } - var_mass_nL.setConstant(true); - var_mass_nL.setVal(extShape.n_left); + if (extShape.n_left != 0.) + { + var_mass_nL.setConstant(true); + var_mass_nL.setVal(extShape.n_left); + } - var_mass_alphaR.setConstant(true); - var_mass_alphaR.setVal(extShape.alpha_right); + if (extShape.alpha_right != 0.) + { + var_mass_alphaR.setConstant(true); + var_mass_alphaR.setVal(extShape.alpha_right); + } - var_mass_nR.setConstant(true); - var_mass_nR.setVal(extShape.n_right); + if (extShape.n_right != 0.) + { + 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"); @@ -452,7 +474,6 @@ 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(); @@ -460,8 +481,8 @@ RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString double bkg_err = var_mass_nbkg.getError(); 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 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("sig_over_bkg", "N_{Sig}/#sqrt{N_{Sig} + N_{Bkg}}", sig_over_bkg_val, sig_over_bkg_err, 2)); @@ -503,7 +524,8 @@ RooFitSummary CreateRooDataSetAndFitCB(TTree *dataSet, TString var_name, TString var_mass_alphaL.getVal(), var_mass_nL.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) diff --git a/mapmc_b02hphmmumu.cpp b/mapmc_b02hphmmumu.cpp index 24ba14d..e281d61 100644 --- a/mapmc_b02hphmmumu.cpp +++ b/mapmc_b02hphmmumu.cpp @@ -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; 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_BPVDIRA", &B0_BPVDIRA_out, "B0_BPVDIRA/F"); output_tree->Branch("B0_PX", &B0_PX_out, "0_PX/F"); diff --git a/new_analysis_b02hphmmumu.cpp b/new_analysis_b02hphmmumu.cpp index d6a0d95..edcb601 100644 --- a/new_analysis_b02hphmmumu.cpp +++ b/new_analysis_b02hphmmumu.cpp @@ -45,39 +45,26 @@ int new_analysis_b02hphmmumu() { 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 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/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)); - 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; TString B_Mass_jpsi_var_name = "B_Mass_jpsi_var"; @@ -94,14 +81,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_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_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); 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); h2_Hlt1_flags_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); h2_Hlt1_flags_excl_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); @@ -115,17 +100,14 @@ int new_analysis_b02hphmmumu() TV::Float("B0_BPVFDCHI2", "B0_BPVFDCHI2"), TV::Float("B0_BPVDIRA", "B0_BPVDIRA"), TV::Float("Jpsi_BPVIPCHI2", "Jpsi_BPVIPCHI2"), - // TV::Float("Jpsi_BPVDIRA", "Jpsi_BPVDIRA"), 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"); @@ -145,31 +127,17 @@ int new_analysis_b02hphmmumu() for (unsigned int i = 0; i < data_entries; 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 +155,28 @@ int new_analysis_b02hphmmumu() { 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) - { - 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) + if (retrain_bdt && std::all_of(vars.begin(), vars.end(), [](TV *v) + { return v->IsMCFinite(); })) { - 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); } if (retrain_bdt) { + std::cout << "# Added " << sig_events << " signal events." << std::endl; TrainBDT(vars, analysis_name, sig_tree, bkg_tree); std::cout << "# Finished BDT retrain." << std::endl; } @@ -240,72 +186,57 @@ int new_analysis_b02hphmmumu() Float_t *train_vars = new Float_t[vars.size()]; 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++) { data_chain->GetEntry(i); 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.) { - 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_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(); } } } @@ -321,16 +252,15 @@ int new_analysis_b02hphmmumu() DrawInDefaultCanvas(h2_Hlt1_flags_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_sim_unf, 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); 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_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_psi2s_fitsum, analysis_name); diff --git a/new_analysis_bu2hpmumu.cpp b/new_analysis_bu2hpmumu.cpp index d79a9b2..f098e03 100644 --- a/new_analysis_bu2hpmumu.cpp +++ b/new_analysis_bu2hpmumu.cpp @@ -45,34 +45,24 @@ int new_analysis_bu2hpmumu() { 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 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/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)); - 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; TString B_Mass_jpsi_var_name = "B_Mass_jpsi_var"; @@ -89,16 +79,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_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_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_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_sim_unf->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_excl_B_Mass->GetXaxis()->SetTitle(end_state_mass_literal); @@ -112,14 +98,12 @@ int new_analysis_bu2hpmumu() TV::Float("B_BPVFDCHI2", "B_BPVFDCHI2"), TV::Float("B_BPVDIRA", "B_BPVDIRA"), TV::Float("Jpsi_BPVIPCHI2", "Jpsi_BPVIPCHI2"), - // TV::Float("Jpsi_BPVDIRA", "Jpsi_BPVDIRA"), 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"); @@ -141,12 +125,13 @@ int new_analysis_bu2hpmumu() for (unsigned int i = 0; i < data_entries; 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) { 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_events++; @@ -166,48 +151,43 @@ int new_analysis_bu2hpmumu() for (unsigned int i = 0; i < sim_entries; 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); } 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 << "# Starting evaluation of data." << std::endl; 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++) { data_chain->GetEntry(i); 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) { return v->IsDataFinite(); })) @@ -264,15 +244,13 @@ int new_analysis_bu2hpmumu() 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_sim_unf, 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); 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_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_psi2s_fitsum, analysis_name); From bcb6249b551c08e651b7a16ffec478957a66904d Mon Sep 17 00:00:00 2001 From: Marius Pfeiffer Date: Tue, 12 Mar 2024 16:56:43 +0100 Subject: [PATCH 2/2] 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