diff --git a/.gitignore b/.gitignore index 38ba3e8..d04eb15 100644 --- a/.gitignore +++ b/.gitignore @@ -4,4 +4,5 @@ status_report/** branchnames.txt images/** output_files/** -BuToHpMuMu_dataloader/** \ No newline at end of file +BuToHpMuMu_dataloader/** +B0ToHpHmMuMu_dataloader/** \ No newline at end of file diff --git a/basic_analysis.h b/basic_analysis.h index 5edba4b..5aacbbd 100644 --- a/basic_analysis.h +++ b/basic_analysis.h @@ -65,6 +65,65 @@ public: } }; +struct FittedParam +{ + std::string title; + std::string name; + double value; + double error; + bool at_limit; + bool constant; + int decimals; + + FittedParam(RooRealVar var, int decimals) + { + this->title = var.GetTitle(); + this->name = var.GetName(); + double val = var.getVal(); + this->value = val; + this->constant = var.isConstant(); + if (!this->constant) + { + this->error = var.getError(); + this->at_limit = ((val == var.getMin()) || (val == var.getMax())); + } + else + { + this->error = 0.; + } + this->decimals = decimals; + } + + std::string ToString() const + { + if (this->constant) + { + return TString::Format("%s = %.*f", this->title.c_str(), this->decimals, this->value).Data(); + } + else + { + return TString::Format("%s = %.*f #pm %.*f", this->title.c_str(), this->decimals, this->value, this->decimals, this->error).Data(); + } + } +}; + +struct ShapeParamters +{ + double alpha_left; + double n_left; + double alpha_right; + double n_right; +}; + +struct RooFitSummary +{ + RooPlot *fit_histogram; + RooPlot *pull_histogram; + std::map pdf_names; + std::vector fitted_params; + ShapeParamters shape_parameters; +}; + void DrawInDefaultCanvas(TH1 *histogram, const char *folder, double margin_left = 0, Option_t *option = "") { std::filesystem::create_directory(TString::Format("output_files/analysis/%s", folder).Data()); @@ -88,6 +147,52 @@ void DrawInDefaultCanvas(RooPlot *rooPlot, const char *folder, double margin_lef c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data()); } +void DrawInDefaultCanvas(RooFitSummary fitSummary, const char *folder) +{ + std::filesystem::create_directory(TString::Format("output_files/analysis/%s", folder).Data()); + TString name = TString::Format("%s_canvas", fitSummary.fit_histogram->GetName()); + TCanvas *c = new TCanvas(name, fitSummary.fit_histogram->GetTitle(), 0, 0, 800, 600); + + TPad *pull_pad = new TPad(TString::Format("%s_canv_pull", fitSummary.fit_histogram->GetName()), "Pull Pad", 0., 0., 1., 0.3); + pull_pad->Draw(); + pull_pad->SetTopMargin(0.001); + pull_pad->SetBottomMargin(0.3); + pull_pad->SetGrid(); + + TPad *fit_pad = new TPad(TString::Format("%s_canv_fit", fitSummary.fit_histogram->GetName()), "Fit Pad", 0., 0.3, 1., 1.); + fit_pad->Draw(); + fit_pad->SetBottomMargin(0.001); + fit_pad->cd(); + + fitSummary.fit_histogram->Draw(); + + TLegend *leg1 = new TLegend(0.60, 0.55, 0.87, 0.89); + leg1->SetFillColor(kWhite); + leg1->SetLineColor(kBlack); + + for (const auto &[name, title] : fitSummary.pdf_names) + { + leg1->AddEntry(fitSummary.fit_histogram->findObject(name.c_str()), title.c_str(), "LP"); + } + + for (const auto &par : fitSummary.fitted_params) + { + if (!par.constant) + { + leg1->AddEntry((TObject *)0, par.ToString().c_str(), ""); + } + } + + leg1->Draw(); + + pull_pad->cd(); + + fitSummary.pull_histogram->Draw(); + + c->Draw(); + c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data()); +} + void PrintProgress(const char *title, unsigned int total, unsigned int every, unsigned int current) { if ((current + 1) % every == 0 || current + 1 == total) @@ -112,21 +217,6 @@ RooPlot *CreateRooFitHistogram(TH1D *hist) return roo_frame_mass; } -struct ShapeParamters -{ - double alpha_left; - double n_left; - double alpha_right; - double n_right; -}; - -struct RooFitSummary -{ - RooPlot *fit_histogram; - RooPlot *pull_histogram; - ShapeParamters shape_parameters; -}; - RooFitSummary CreateRooFitHistogramAndFitCB(TH1D *hist, bool hasExpBkg, bool useExtShape, ShapeParamters extShape) { auto suffix_name = [name = hist->GetName()](const char *text) @@ -154,7 +244,8 @@ RooFitSummary CreateRooFitHistogramAndFitCB(TH1D *hist, bool hasExpBkg, bool use RooRealVar var_mass_alphaR(suffix_name("var_mass_alphaR"), "#alpha_{R}", 2., 0., 4.); RooRealVar var_mass_nR(suffix_name("var_mass_nR"), "n_{R}", 5., 0., 15.); - if (useExtShape) { + if (useExtShape) + { var_mass_alphaL.setConstant(true); var_mass_alphaL.setVal(extShape.alpha_left); @@ -165,11 +256,16 @@ RooFitSummary CreateRooFitHistogramAndFitCB(TH1D *hist, bool hasExpBkg, bool use var_mass_alphaR.setVal(extShape.alpha_right); var_mass_nR.setConstant(true); - var_mass_nR.setVal(extShape.n_left); + var_mass_nR.setVal(extShape.n_right); } TString signal_name = suffix_name("sig_cb"); - RooCrystalBall sig_cb(signal_name, "Signal Crystal Ball", 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); + + std::vector fitted_params{}; + std::map pdf_names{}; + + pdf_names[signal_name.Data()] = sig_cb.getTitle().Data(); TString pull_compare_name{}; if (hasExpBkg) @@ -179,12 +275,14 @@ RooFitSummary CreateRooFitHistogramAndFitCB(TH1D *hist, bool hasExpBkg, bool use TString background_name = suffix_name("bgk_exp"); RooExponential bkg_exp(background_name, "Exp Background", roo_var_mass, var_mass_bkg_c); + pdf_names[background_name.Data()] = bkg_exp.getTitle().Data(); RooRealVar var_mass_nsig(suffix_name("nsig"), "N_{Sig}", 0., hist->GetEntries()); RooRealVar var_mass_nbkg(suffix_name("nbkg"), "N_{Bkg}", 0., hist->GetEntries()); TString fitted_pdf_name = suffix_name("sigplusbkg"); - RooAddPdf fitted_pdf(fitted_pdf_name, "Sig and Bkg PDF", RooArgList(sig_cb, bkg_exp), RooArgList(var_mass_nsig, var_mass_nbkg)); + RooAddPdf fitted_pdf(fitted_pdf_name, "Sig + Bkg", RooArgList(sig_cb, bkg_exp), RooArgList(var_mass_nsig, var_mass_nbkg)); + pdf_names[fitted_pdf_name.Data()] = fitted_pdf.getTitle().Data(); RooFitResult *fitres = fitted_pdf.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range(fitting_range_name)); @@ -192,26 +290,36 @@ RooFitSummary CreateRooFitHistogramAndFitCB(TH1D *hist, bool hasExpBkg, bool use fitted_pdf.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range(fitting_range_name), RooFit::Name(fitted_pdf_name)); fitted_pdf.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(bkg_exp)), RooFit::LineColor(kBlue - 7), RooFit::LineStyle(kDashed), RooFit::Range(fitting_range_name), RooFit::Name(background_name)); fitted_pdf.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(sig_cb)), RooFit::LineColor(kRed - 7), RooFit::LineStyle(kDashed), RooFit::Range(fitting_range_name), RooFit::Name(signal_name)); - fitted_pdf.paramOn(roo_frame_mass); + + fitted_params.push_back(FittedParam(var_mass_bkg_c, 5)); + fitted_params.push_back(FittedParam(var_mass_nsig, 2)); + fitted_params.push_back(FittedParam(var_mass_nbkg, 2)); pull_compare_name = fitted_pdf_name; } else { RooFitResult *fitres = sig_cb.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range(fitting_range_name)); - // sigplusbkg.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144)); sig_cb.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range(fitting_range_name), RooFit::Name(signal_name)); pull_compare_name = signal_name; - sig_cb.paramOn(roo_frame_mass); } + fitted_params.push_back(FittedParam(var_mass_x0, 2)); + fitted_params.push_back(FittedParam(var_mass_sigmaLR, 2)); + fitted_params.push_back(FittedParam(var_mass_alphaL, 2)); + fitted_params.push_back(FittedParam(var_mass_nL, 2)); + fitted_params.push_back(FittedParam(var_mass_alphaR, 2)); + fitted_params.push_back(FittedParam(var_mass_nR, 2)); + RooPlot *roo_frame_pull = roo_var_mass.frame(RooFit::Title("Pull Distribution"), RooFit::Name(TString::Format("%s_pull_rplt", hist->GetName()))); roo_frame_pull->addPlotable(roo_frame_mass->pullHist(hist_name, pull_compare_name), "P"); return RooFitSummary{ roo_frame_mass, roo_frame_pull, + pdf_names, + fitted_params, ShapeParamters{ var_mass_alphaL.getVal(), var_mass_nL.getVal(), diff --git a/bdt_classification.h b/bdt_classification.h index b83fad0..7e717b2 100644 --- a/bdt_classification.h +++ b/bdt_classification.h @@ -202,7 +202,7 @@ void ConnectVarsToData(std::vector vars, TChain *data_chain, TChain *mc_ch void TrainBDT(std::vector vars, const char* unique_id, TTree *sig_tree, TTree *bkg_tree) { - TString outfile_name = TString::Format("%s_out.root", unique_id); + TString outfile_name = TString::Format("%s_tmva_out.root", unique_id); TFile *output_file = TFile::Open(outfile_name, "RECREATE"); TString factory_options("V:Silent:Color:DrawProgressBar:Transformations=I;D;P;G,D:AnalysisType=Auto"); diff --git a/new_analysis_b02hphmmumu.cpp b/new_analysis_b02hphmmumu.cpp new file mode 100644 index 0000000..f681d80 --- /dev/null +++ b/new_analysis_b02hphmmumu.cpp @@ -0,0 +1,309 @@ +#include +#include +#include +#include +#include +#include + +#include "TH1D.h" +#include "TH2D.h" +#include "THStack.h" +#include "TGraph.h" +#include "TTree.h" +#include "TChain.h" +#include "TFile.h" +#include "TCanvas.h" +#include "TROOT.h" +#include "TStyle.h" +#include "TColor.h" +#include "TLorentzVector.h" +#include "TRandom3.h" +#include "TLegend.h" + +#include "RooDataHist.h" +#include "RooRealVar.h" +#include "RooPlot.h" +#include "RooGaussian.h" +#include "RooExponential.h" +#include "RooRealConstant.h" +#include "RooAddPdf.h" +#include "RooFitResult.h" +#include "RooProduct.h" +#include "RooCrystalBall.h" +#include "RooBreitWigner.h" +#include "RooArgSet.h" +#include "RooFFTConvPdf.h" +#include "RooNovosibirsk.h" + +#include "constants.h" +#include "basic_analysis.h" +#include "hlt1_decision_analysis.h" +#include "bdt_classification.h" + +// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Inclusive :: B0 To Hp Hm Mu Mu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +int new_analysis_b02hphmmumu() +{ + const char *analysis_name = "B0ToHpHmMuMu"; + const char *data_tree_name = "B0ToHpHmMuMu"; + const char *sim_tree_name = "B0ToHpHmMuMu_noPID"; + const char *end_state_mass_literal = "m(#pi^{+}_{(#rightarrow K^{+})}#pi^{-}_{(#rightarrow K^{-})}#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/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"); + + 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"); + + 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->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); + + TH1D *h1_B_Mass_jpsi = new TH1D("h1_B_Mass_jpsi", TString::Format("B Mass, J/#psi Mode (%s)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TH1D *h1_B_Mass_psi2s = new TH1D("h1_B_Mass_psi2s", TString::Format("B Mass, #psi(2S) Mode (%s)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TH1D *h1_B_Mass_sim = new TH1D("h1_B_Mass_sim", TString::Format("B Mass, Simualted (%s_noPID)", analysis_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_jpsi->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_psi2s->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); + + ConnectHlt1Decisions(data_chain, h2_Hlt1_flags_B_Mass, h2_Hlt1_flags_excl_B_Mass); + auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name); + std::map exclusive_hits{}; + + std::vector vars{ + TV::Float("B0_PT", "B0_PT"), + 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"), + }; + + TTree *sig_tree = new TTree("TreeS", "tree containing signal data"); + TTree *bkg_tree = new TTree("TreeB", "tree containing background data"); + + ConnectVarsToData(vars, data_chain, sim_chain, sig_tree, bkg_tree); + + unsigned int data_entries = data_chain->GetEntries(); + unsigned int sim_entries = 100000; // sim_chain->GetEntries(); + + unsigned int bkg_events = 0; + + if (retrain_bdt) + { + std::cout << "# Starting with BDT retrain." << std::endl; + 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; + } + + 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) + { return v->IsDataFinite(); })) + { + if (reconstructed_B_Mass > 5500.) + { + bkg_tree->Fill(); + bkg_events++; + } + } + } + + PrintProgress(TString::Format("%s BKG Collection", analysis_name), data_entries, 10000, i); + } + } + else + { + std::cout << "# Starting without BDT retrain." << std::endl; + } + + for (unsigned int i = 0; i < sim_entries; i++) + { + sim_chain->GetEntry(i); + + 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(); + + if (std::all_of(vars.begin(), vars.end(), [](TV *v) + { return v->IsMCFinite(); })) + { + if (retrain_bdt) + { + sig_tree->Fill(); + } + h1_B_Mass_sim->Fill(reconstructed_B_Mass); + } + } + } + + PrintProgress(TString::Format("%s SIG Collection", analysis_name), sim_entries, 10000, i); + } + + if (retrain_bdt) + { + 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, analysis_name); + + const double mva_cut_value = -0.12; // 0.09; // -0.02; + + 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; + } + + if (found_k_star) + { + 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++) + { + if (vars[j]->IsDouble()) + { + train_vars[j] = vars[j]->GetDataDouble(); + } + else if (vars[j]->IsFloat()) + { + 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); + } + + double mva_response = reader->EvaluateMVA("BDT"); + h1_bdt_probs->Fill(mva_response); + + if (mva_response > mva_cut_value) + { + if (TMath::Abs(reconstructed_Kstar.M() - KSTAR_MASS) < 100.) + { + if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) + { + h1_B_Mass_jpsi->Fill(reconstructed_B_Mass); + } + else if (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.) + { + h1_B_Mass_psi2s->Fill(reconstructed_B_Mass); + } + } + } + } + + PrintProgress(TString::Format("%s BDT Evaluation", analysis_name), data_entries, 10000, i); + } + + std::cout << "# Exclusive Hits" << std::endl; + for (const auto &[line, hits] : exclusive_hits) + { + std::cout << line << ": " << hits << std::endl; + } + + // 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_jpsi, analysis_name, 0.1); + DrawInDefaultCanvas(h1_B_Mass_psi2s, analysis_name, 0.1); + + auto roofit_hist_sim = CreateRooFitHistogramAndFitCB(h1_B_Mass_sim, false, false, ShapeParamters{}); + auto roofit_hist_jpsi_fitsum = CreateRooFitHistogramAndFitCB(h1_B_Mass_jpsi, true, true, roofit_hist_sim.shape_parameters); + auto roofit_hist_psi2s_fitsum = CreateRooFitHistogramAndFitCB(h1_B_Mass_psi2s, true, true, roofit_hist_sim.shape_parameters); + + DrawInDefaultCanvas(roofit_hist_jpsi_fitsum, analysis_name); + DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name); + DrawInDefaultCanvas(roofit_hist_sim, analysis_name); + + // DrawHlt1DecisionHistos(analysis_name, hlt1_decision_histos); + + DrawBDTProbs(h1_bdt_probs, mva_cut_value, analysis_name); + + return 0; +} \ No newline at end of file diff --git a/new_analysis_b02kppimmumu.cpp b/new_analysis_b02kppimmumu.cpp new file mode 100644 index 0000000..3921230 --- /dev/null +++ b/new_analysis_b02kppimmumu.cpp @@ -0,0 +1,160 @@ +#include +#include +#include +#include +#include +#include + +#include "TH1D.h" +#include "TH2D.h" +#include "THStack.h" +#include "TGraph.h" +#include "TTree.h" +#include "TChain.h" +#include "TFile.h" +#include "TCanvas.h" +#include "TROOT.h" +#include "TStyle.h" +#include "TColor.h" +#include "TLorentzVector.h" +#include "TRandom3.h" +#include "TLegend.h" + +#include "RooDataHist.h" +#include "RooRealVar.h" +#include "RooPlot.h" +#include "RooGaussian.h" +#include "RooExponential.h" +#include "RooRealConstant.h" +#include "RooAddPdf.h" +#include "RooFitResult.h" +#include "RooProduct.h" +#include "RooCrystalBall.h" +#include "RooBreitWigner.h" +#include "RooArgSet.h" +#include "RooFFTConvPdf.h" +#include "RooNovosibirsk.h" + +#include "constants.h" +#include "basic_analysis.h" +#include "hlt1_decision_analysis.h" +#include "bdt_classification.h" + +// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Exclusive B0 To Kp Pim Mu Mu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +int new_analysis_b02kppimmumu() +{ + const char *analysis_name = "B0ToKpPimMuMu"; + 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"); + + FourVect *l14v_data = FourVect::Init(data_chain, "muplus"); + FourVect *l24v_data = FourVect::Init(data_chain, "muminus"); + 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_turbo_v0r0p6657752_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, "K"); + FourVect *hm4v_sim = FourVect::Init(sim_chain, "Pi"); + + 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("K_TRUEID", &Hp_TRUEID); + sim_chain->SetBranchAddress("Pi_TRUEID", &Hm_TRUEID); + sim_chain->SetBranchAddress("B0_BKGCAT", &B_BKGCAT); + + TH1D *h1_B_Mass_jpsi = new TH1D("h1_B_Mass_jpsi", TString::Format("B Mass, J/#psi Mode (%s)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TH1D *h1_B_Mass_psi2s = new TH1D("h1_B_Mass_psi2s", TString::Format("B Mass, #psi(2S) Mode (%s)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TH1D *h1_B_Mass_sim = new TH1D("h1_B_Mass_sim", TString::Format("B Mass, Simualted (%s_noPID)", analysis_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.); + + h1_B_Mass_jpsi->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_psi2s->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); + + ConnectHlt1Decisions(data_chain, h2_Hlt1_flags_B_Mass, h2_Hlt1_flags_excl_B_Mass); + auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name); + std::map exclusive_hits{}; + + unsigned int data_entries = data_chain->GetEntries(); + unsigned int sim_entries = sim_chain->GetEntries(); + + for (unsigned int i = 0; i < sim_entries; i++) + { + sim_chain->GetEntry(i); + Double_t reconstructed_B_Mass = (hp4v_sim->LorentzVector() + hm4v_sim->LorentzVector() + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M(); + if (B_BKGCAT == 0 && TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID && TMath::Abs(Hp_TRUEID) == PID_KAON && TMath::Abs(Hm_TRUEID) == PID_PION) + { + h1_B_Mass_sim->Fill(reconstructed_B_Mass); + } + + PrintProgress(TString::Format("%s SIG Collection", analysis_name), sim_entries, 10000, i); + } + + std::cout << "# Starting evaluation of data." << std::endl; + + for (unsigned int i = 0; i < data_entries; i++) + { + data_chain->GetEntry(i); + + TLorentzVector dimuon = l14v_data->LorentzVector() + l24v_data->LorentzVector(); + TLorentzVector reconstructed_Kstar = hp4v_data->LorentzVector() + hm4v_data->LorentzVector(); + Double_t reconstructed_B_Mass = (reconstructed_Kstar + dimuon).M(); + + 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(reconstructed_Kstar.M() - KSTAR_MASS) < 100.) + { + if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) + { + h1_B_Mass_jpsi->Fill(reconstructed_B_Mass); + } + else if (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.) + { + h1_B_Mass_psi2s->Fill(reconstructed_B_Mass); + } + } + + PrintProgress(TString::Format("%s Data Evaluation", analysis_name), data_entries, 10000, i); + } + + std::cout << "# Exclusive Hits" << std::endl; + for (const auto &[line, hits] : exclusive_hits) + { + std::cout << line << ": " << hits << std::endl; + } + + 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_jpsi, analysis_name, 0.1); + DrawInDefaultCanvas(h1_B_Mass_psi2s, analysis_name, 0.1); + + auto roofit_hist_sim = CreateRooFitHistogramAndFitCB(h1_B_Mass_sim, false, false, ShapeParamters{}); + auto roofit_hist_jpsi_fitsum = CreateRooFitHistogramAndFitCB(h1_B_Mass_jpsi, true, true, roofit_hist_sim.shape_parameters); + auto roofit_hist_psi2s_fitsum = CreateRooFitHistogramAndFitCB(h1_B_Mass_psi2s, true, true, roofit_hist_sim.shape_parameters); + + DrawInDefaultCanvas(roofit_hist_jpsi_fitsum, analysis_name); + DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name); + DrawInDefaultCanvas(roofit_hist_sim, analysis_name); + + DrawHlt1DecisionHistos(analysis_name, hlt1_decision_histos); + + return 0; +} \ No newline at end of file diff --git a/new_analysis.cpp b/new_analysis_bu2hpmumu.cpp similarity index 96% rename from new_analysis.cpp rename to new_analysis_bu2hpmumu.cpp index d2ba7a1..ae260fd 100644 --- a/new_analysis.cpp +++ b/new_analysis_bu2hpmumu.cpp @@ -42,21 +42,23 @@ // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Bu To Hp Mu Mu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -int analyze_bu2hpmumu_data() +int new_analysis_bu2hpmumu() { const char *analysis_name = "BuToHpMuMu"; + const char *data_tree_name = "BuToHpMuMu"; + const char *sim_tree_name = "BuToHpMuMu_noPID"; 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", analysis_name)); - data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0_tuple_90000000_v0r0p6288631.root"); + 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"); FourVect *l14v_data = FourVect::Init(data_chain, "L1"); FourVect *l24v_data = FourVect::Init(data_chain, "L2"); FourVect *hp4v_data = FourVect::Init(data_chain, "Hp"); - TChain *sim_chain = new TChain(TString::Format("%s_noPID/DecayTree", analysis_name)); + 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"); FourVect *l14v_sim = FourVect::Init(sim_chain, "L1"); @@ -107,7 +109,7 @@ int analyze_bu2hpmumu_data() ConnectVarsToData(vars, data_chain, sim_chain, sig_tree, bkg_tree); unsigned int data_entries = data_chain->GetEntries(); - unsigned int sim_entries = sim_chain->GetEntries(); + unsigned int sim_entries = data_entries;// sim_chain->GetEntries(); unsigned int bkg_events = 0; @@ -167,7 +169,7 @@ int analyze_bu2hpmumu_data() Float_t *train_vars = new Float_t[vars.size()]; auto reader = SetupReader(vars, train_vars, "BuToHpMuMu"); - const double mva_cut_value = 0; // 0.09; // -0.02; + const double mva_cut_value = -0.03; // 0.09; // -0.02; for (unsigned int i = 0; i < data_entries; i++) { @@ -230,20 +232,15 @@ int analyze_bu2hpmumu_data() auto roofit_hist_sim = CreateRooFitHistogramAndFitCB(h1_B_Mass_sim, false, false, ShapeParamters{}); auto roofit_hist_jpsi_fitsum = CreateRooFitHistogramAndFitCB(h1_B_Mass_jpsi, true, true, roofit_hist_sim.shape_parameters); auto roofit_hist_psi2s_fitsum = CreateRooFitHistogramAndFitCB(h1_B_Mass_psi2s, true, true, roofit_hist_sim.shape_parameters); - - DrawInDefaultCanvas(roofit_hist_jpsi_fitsum.fit_histogram, analysis_name, 0.1); - DrawInDefaultCanvas(roofit_hist_psi2s_fitsum.fit_histogram, analysis_name, 0.1); - DrawInDefaultCanvas(roofit_hist_sim.fit_histogram, analysis_name, 0.1); + + DrawInDefaultCanvas(roofit_hist_jpsi_fitsum, analysis_name); + DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name); + DrawInDefaultCanvas(roofit_hist_sim, analysis_name); DrawHlt1DecisionHistos(analysis_name, hlt1_decision_histos); DrawBDTProbs(h1_bdt_probs, mva_cut_value, analysis_name); - // auto test_bkg_tree_file = TFile::Open("test_bkg_tree.root", "RECREATE"); - // bkg_tree->Write(); - // sig_tree->Write(); - // test_bkg_tree_file->Close(); - return 0; } @@ -621,7 +618,7 @@ int analyze_Hlt2RD_B0ToKpPimMuMu() return 0; } -int new_analysis() -{ - return analyze_bu2hpmumu_data(); -} \ No newline at end of file +// int new_analysis() +// { +// return analyze_bu2hpmumu_data(); +// } \ No newline at end of file diff --git a/new_analysis_bu2kpmumu.cpp b/new_analysis_bu2kpmumu.cpp new file mode 100644 index 0000000..0adf66a --- /dev/null +++ b/new_analysis_bu2kpmumu.cpp @@ -0,0 +1,153 @@ +#include +#include +#include +#include +#include +#include + +#include "TH1D.h" +#include "TH2D.h" +#include "THStack.h" +#include "TGraph.h" +#include "TTree.h" +#include "TChain.h" +#include "TFile.h" +#include "TCanvas.h" +#include "TROOT.h" +#include "TStyle.h" +#include "TColor.h" +#include "TLorentzVector.h" +#include "TRandom3.h" +#include "TLegend.h" + +#include "RooDataHist.h" +#include "RooRealVar.h" +#include "RooPlot.h" +#include "RooGaussian.h" +#include "RooExponential.h" +#include "RooRealConstant.h" +#include "RooAddPdf.h" +#include "RooFitResult.h" +#include "RooProduct.h" +#include "RooCrystalBall.h" +#include "RooBreitWigner.h" +#include "RooArgSet.h" +#include "RooFFTConvPdf.h" +#include "RooNovosibirsk.h" + +#include "constants.h" +#include "basic_analysis.h" +#include "hlt1_decision_analysis.h" +#include "bdt_classification.h" + +// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Exclusive Bu To Kp Mu Mu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +int new_analysis_bu2kpmumu() +{ + const char *analysis_name = "BuToKpMuMu"; + 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 = 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"); + + FourVect *l14v_data = FourVect::Init(data_chain, "muplus"); + FourVect *l24v_data = FourVect::Init(data_chain, "muminus"); + 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_turbo_v0r0p6657752_BuToKpMuMu_12143001_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, "K"); + + Int_t B_BKGCAT, L1_TRUEID, L2_TRUEID, Hp_TRUEID; + + sim_chain->SetBranchAddress("L1_TRUEID", &L1_TRUEID); + sim_chain->SetBranchAddress("L2_TRUEID", &L2_TRUEID); + sim_chain->SetBranchAddress("K_TRUEID", &Hp_TRUEID); + sim_chain->SetBranchAddress("B_BKGCAT", &B_BKGCAT); + + TH1D *h1_B_Mass_jpsi = new TH1D("h1_B_Mass_jpsi", TString::Format("B Mass, J/#psi Mode (%s)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TH1D *h1_B_Mass_psi2s = new TH1D("h1_B_Mass_psi2s", TString::Format("B Mass, #psi(2S) Mode (%s)", analysis_name), B_MASS_HIST_BINS, B_MASS_VAR_MIN, B_MASS_VAR_MAX); + TH1D *h1_B_Mass_sim = new TH1D("h1_B_Mass_sim", TString::Format("B Mass, Simualted (%s_noPID)", analysis_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.); + + h1_B_Mass_jpsi->GetXaxis()->SetTitle(end_state_mass_literal); + h1_B_Mass_psi2s->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); + + ConnectHlt1Decisions(data_chain, h2_Hlt1_flags_B_Mass, h2_Hlt1_flags_excl_B_Mass); + auto hlt1_decision_histos = CreateHlt1DecisionHistos(analysis_name); + std::map exclusive_hits{}; + + unsigned int data_entries = data_chain->GetEntries(); + unsigned int sim_entries = sim_chain->GetEntries(); + + for (unsigned int i = 0; i < sim_entries; i++) + { + sim_chain->GetEntry(i); + Double_t reconstructed_B_Mass = (hp4v_sim->LorentzVector() + l14v_sim->LorentzVector() + l24v_sim->LorentzVector()).M(); + if (B_BKGCAT == 0 && TMath::Abs(L1_TRUEID) == PID_MUON && L2_TRUEID == -L1_TRUEID && TMath::Abs(Hp_TRUEID) == PID_KAON) + { + h1_B_Mass_sim->Fill(reconstructed_B_Mass); + } + + PrintProgress(TString::Format("%s SIG Collection", analysis_name), sim_entries, 10000, i); + } + + std::cout << "# Starting evaluation of data." << std::endl; + + 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() + dimuon).M(); + + 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.) + { + h1_B_Mass_jpsi->Fill(reconstructed_B_Mass); + } + else if (TMath::Abs(dimuon.M() - PSI2S_MASS) < 100.) + { + h1_B_Mass_psi2s->Fill(reconstructed_B_Mass); + } + + PrintProgress(TString::Format("%s BDT Evaluation", analysis_name), data_entries, 10000, i); + } + + std::cout << "# Exclusive Hits" << std::endl; + for (const auto &[line, hits] : exclusive_hits) + { + std::cout << line << ": " << hits << std::endl; + } + + 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_jpsi, analysis_name, 0.1); + DrawInDefaultCanvas(h1_B_Mass_psi2s, analysis_name, 0.1); + + auto roofit_hist_sim = CreateRooFitHistogramAndFitCB(h1_B_Mass_sim, false, false, ShapeParamters{}); + auto roofit_hist_jpsi_fitsum = CreateRooFitHistogramAndFitCB(h1_B_Mass_jpsi, true, true, roofit_hist_sim.shape_parameters); + auto roofit_hist_psi2s_fitsum = CreateRooFitHistogramAndFitCB(h1_B_Mass_psi2s, true, true, roofit_hist_sim.shape_parameters); + + DrawInDefaultCanvas(roofit_hist_jpsi_fitsum, analysis_name); + DrawInDefaultCanvas(roofit_hist_psi2s_fitsum, analysis_name); + DrawInDefaultCanvas(roofit_hist_sim, analysis_name); + + DrawHlt1DecisionHistos(analysis_name, hlt1_decision_histos); + + return 0; +} \ No newline at end of file