diff --git a/.gitignore b/.gitignore index d7f4219..df8967a 100644 --- a/.gitignore +++ b/.gitignore @@ -2,4 +2,5 @@ dataloader/** status_report/** branchnames.txt -images/** \ No newline at end of file +images/** +output_files/** \ No newline at end of file diff --git a/new_analysis.cpp b/new_analysis.cpp new file mode 100644 index 0000000..7126da6 --- /dev/null +++ b/new_analysis.cpp @@ -0,0 +1,732 @@ +#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 +#include +#include +#include +#include + +const Double_t B_MASS_VAR_MIN = 5100.; +const Double_t B_MASS_VAR_MAX = 6000.; +const Int_t B_MASS_HIST_BINS = 70; + +const Double_t K_MASS = 493.677; +const Double_t PSI2S_MASS = 3686.093; +const Double_t JPSI_MASS = 3096.900; +const Double_t KSTAR_MASS = 891.760; + +class FourVect +{ +private: + Float_t px, py, pz, energy; + + Float_t *PtrPX() + { + return &px; + } + + Float_t *PtrPY() + { + return &py; + } + + Float_t *PtrPZ() + { + return &pz; + } + + Float_t *PtrEnergy() + { + return &energy; + } + +public: + static FourVect *Init(TTree *tree, const char *particle) + { + FourVect *v = new FourVect{}; + tree->SetBranchAddress(TString::Format("%s_PX", particle).Data(), v->PtrPX()); + tree->SetBranchAddress(TString::Format("%s_PY", particle).Data(), v->PtrPY()); + tree->SetBranchAddress(TString::Format("%s_PZ", particle).Data(), v->PtrPZ()); + tree->SetBranchAddress(TString::Format("%s_ENERGY", particle).Data(), v->PtrEnergy()); + return v; + } + + TLorentzVector LorentzVector() + { + return TLorentzVector(px, py, pz, energy); + } + + TLorentzVector LorentzVector(double sub_mass) + { + TVector3 momentum(px, py, pz); + float energy = TMath::Sqrt(TMath::Sq(sub_mass) + momentum.Mag2()); + return TLorentzVector(momentum, energy); + } + + std::string ToString() + { + return TString::Format("(%f, %f, %f, %f)", px, py, pz, energy).Data(); + } +}; + +struct Hlt1Decision +{ + std::string name; + int index; + Bool_t value; + + std::string GetName() const + { + return TString::Format("Hlt1%sDecision", name.c_str()).Data(); + } + + Bool_t *GetValuePointer() + { + return &value; + } +}; + +std::vector Hlt1Decisions{ + Hlt1Decision{"DiMuonHighMass", 1}, + Hlt1Decision{"DiMuonLowMass", 2}, + Hlt1Decision{"DiMuonNoIP", 3}, + Hlt1Decision{"DiMuonSoft", 4}, + Hlt1Decision{"DisplacedLeptons", 5}, + Hlt1Decision{"LowPtDiMuon", 6}, + Hlt1Decision{"LowPtMuon", 7}, + Hlt1Decision{"OneMuonTrackLine", 8}, + Hlt1Decision{"SingleHighEt", 9}, + Hlt1Decision{"SingleHighPtMuon", 10}, + Hlt1Decision{"TrackMVA", 11}, + Hlt1Decision{"TrackMuonMVA", 12}, + Hlt1Decision{"TwoTrackMVA", 13}, +}; + +bool CutHlt1DecisionsAnd(const std::set &decision_list) +{ + bool okay = true; + for (const auto &var : Hlt1Decisions) + { + if (decision_list.find(var.name) != decision_list.end()) + { + okay = okay && var.value; + } + } + return okay; +} + +bool CutHlt1DecisionsOr(const std::set &decision_list) +{ + bool okay = false; + for (const auto &var : Hlt1Decisions) + { + if (decision_list.find(var.name) != decision_list.end()) + { + okay = okay || var.value; + } + } + return okay; +} + +void ConnectHlt1Decisions(TChain *chain) +{ + for (auto &var : Hlt1Decisions) + { + if (chain->FindBranch(var.GetName().c_str())) + chain->SetBranchAddress(var.GetName().c_str(), var.GetValuePointer()); + } +} + +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()); + TString name = TString::Format("%s_canvas", histogram->GetName()); + TCanvas *c = new TCanvas(name, histogram->GetName(), 0, 0, 800, 600); + c->SetLeftMargin(margin_left); + histogram->SetStats(0); + histogram->Draw(option); + c->Draw(); + c->SaveAs(TString::Format("output_files/analysis/%s/%s.pdf", folder, name.Data()).Data()); +} + +void DrawInDefaultCanvas(RooPlot *rooPlot, const char *folder, double margin_left = 0, Option_t *option = "") +{ + std::filesystem::create_directory(TString::Format("output_files/analysis/%s", folder).Data()); + TString name = TString::Format("%s_canvas", rooPlot->GetName()); + TCanvas *c = new TCanvas(name, rooPlot->GetName(), 0, 0, 800, 600); + c->SetLeftMargin(margin_left); + rooPlot->Draw(option); + 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) + { + std::cout << "[" + << title + << "] Processed event: " << current + 1 << " (" << TString::Format("%.2f", ((double)(current + 1) / (double)total) * 100.) << "%)" << std::endl; + } +} + +RooPlot* CreateRooFitHistogram(TH1D* hist) { + RooRealVar roo_var_mass("var_mass", "B Mass Variable", B_MASS_VAR_MIN, B_MASS_VAR_MAX); + roo_var_mass.setRange("fitting_range", B_MASS_VAR_MIN, B_MASS_VAR_MAX); + + RooDataHist roohist_B_M("roohist_B_M", "B Mass Histogram", roo_var_mass, RooFit::Import(*hist)); + + RooPlot *roo_frame_mass = roo_var_mass.frame(RooFit::Title(hist->GetTitle())); + roohist_B_M.plotOn(roo_frame_mass, RooFit::Binning(B_MASS_HIST_BINS), RooFit::Name("B Mass Distribution")); + roo_frame_mass->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle()); + + return roo_frame_mass; +} + +// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Bu To Hp Mu Mu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +int analyze_bu2hpmumu() +{ + const char *analysis_name = "BuToHpMuMu"; + const char *end_state_mass_literal = "m(#pi^{+}_{(#rightarrow K^{+})}#mu^{+}#mu^{-})"; + + 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"); + data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0r1_tuple_90000000_2023_v0r0p6288631.root"); + + FourVect *l14v = FourVect::Init(data_chain, "L1"); + FourVect *l24v = FourVect::Init(data_chain, "L2"); + FourVect *hp4v = FourVect::Init(data_chain, "Hp"); + + 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); + 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); + + for (const auto &var : Hlt1Decisions) + { + if (data_chain->FindBranch(var.GetName().c_str())) + { + h2_Hlt1_flags_B_Mass->Fill(0., var.name.c_str(), 0.); + h2_Hlt1_flags_excl_B_Mass->Fill(0., var.name.c_str(), 0.); + } + } + + std::map exclusive_hits{}; + + unsigned int entries = data_chain->GetEntries(); + for (unsigned int i = 0; i < entries; i++) + { + data_chain->GetEntry(i); + + TLorentzVector dimuon = l14v->LorentzVector() + l24v->LorentzVector(); + Double_t reconstructed_B_Mass = (hp4v->LorentzVector(K_MASS) + dimuon).M(); + + if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) + { + for (const auto &var : Hlt1Decisions) + { + if (var.value) + { + h2_Hlt1_flags_B_Mass->Fill(reconstructed_B_Mass, var.name.c_str(), 1); + } + } + + bool exclusive = true; + std::string line{}; + for (const auto &var : Hlt1Decisions) + { + if (var.value) + { + if (!line.empty()) + { + exclusive = false; + } + + line = var.name; + } + } + + if (!line.empty() && exclusive) + { + int &hits = exclusive_hits[line]; + if (hits) + { + hits += 1; + } + else + { + hits = 1; + } + h2_Hlt1_flags_excl_B_Mass->Fill(reconstructed_B_Mass, line.c_str(), 1); + } + } + + if (CutHlt1DecisionsOr({"TwoTrackMuonMVA", "TrackMuonMVA"})) + { + 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(analysis_name, 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_jpsi = CreateRooFitHistogram(h1_B_Mass_jpsi); + auto roofit_hist_psi2s = CreateRooFitHistogram(h1_B_Mass_psi2s); + DrawInDefaultCanvas(roofit_hist_jpsi, analysis_name, 0.1); + DrawInDefaultCanvas(roofit_hist_psi2s, analysis_name, 0.1); + + return 0; +} + +// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% B0 To Hp Hm Mu Mu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +int analyze_b02hphmmumu() +{ + const char *analysis_name = "B0ToHpHmMuMu"; + const char *end_state_mass_literal = "m(#pi^{+}_{(#rightarrow K^{+})} #pi^{-} #mu^{+}#mu^{-})"; + + 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"); + 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, Hp_PID_PI, Hm_PID_PI; + + data_chain->SetBranchAddress("Hp_PID_K", &Hp_PID_K); + data_chain->SetBranchAddress("Hm_PID_K", &Hm_PID_K); + data_chain->SetBranchAddress("Hp_PID_PI", &Hp_PID_PI); + data_chain->SetBranchAddress("Hm_PID_PI", &Hm_PID_PI); + + FourVect *l14v = FourVect::Init(data_chain, "L1"); + FourVect *l24v = FourVect::Init(data_chain, "L2"); + FourVect *hp4v = FourVect::Init(data_chain, "Hp"); + FourVect *hm4v = FourVect::Init(data_chain, "Hm"); + + 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); + 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); + + for (const auto &var : Hlt1Decisions) + { + h2_Hlt1_flags_B_Mass->Fill(0., var.name.c_str(), 0.); + h2_Hlt1_flags_excl_B_Mass->Fill(0., var.name.c_str(), 0.); + } + + std::map exclusive_hits{}; + + unsigned int entries = data_chain->GetEntries(); + for (unsigned int i = 0; i < entries; i++) + { + data_chain->GetEntry(i); + + TLorentzVector dimuon = l14v->LorentzVector() + l24v->LorentzVector(); + + TLorentzVector lv_hp{}; + TLorentzVector lv_hm{}; + + if (Hp_PID_K > 0 && Hm_PID_K < 0) + { + lv_hp = hp4v->LorentzVector(K_MASS); + lv_hm = hm4v->LorentzVector(); + } + else if (Hp_PID_K < 0 && Hm_PID_K > 0) + { + continue; + lv_hp = hp4v->LorentzVector(); + lv_hm = hm4v->LorentzVector(K_MASS); + } + else + { + continue; + } + + Double_t reconstructed_Kstar_Mass = (lv_hp + lv_hm).M(); + Double_t reconstructed_B_Mass = (lv_hp + lv_hm + dimuon).M(); + + if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) + { + for (const auto &var : Hlt1Decisions) + { + if (var.value) + { + h2_Hlt1_flags_B_Mass->Fill(reconstructed_B_Mass, var.name.c_str(), 1); + } + } + + bool exclusive = true; + std::string line{}; + for (const auto &var : Hlt1Decisions) + { + if (var.value) + { + if (!line.empty()) + { + exclusive = false; + } + + line = var.name; + } + } + + if (!line.empty() && exclusive) + { + int &hits = exclusive_hits[line]; + if (hits) + { + hits += 1; + } + else + { + hits = 1; + } + h2_Hlt1_flags_excl_B_Mass->Fill(reconstructed_B_Mass, line.c_str(), 1); + } + } + + if (CutHlt1DecisionsOr({"TwoTrackMuonMVA", "TrackMuonMVA"})) + { + if (TMath::Abs(reconstructed_Kstar_Mass - KSTAR_MASS) < 50.) + { + 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(analysis_name, 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_jpsi = CreateRooFitHistogram(h1_B_Mass_jpsi); + auto roofit_hist_psi2s = CreateRooFitHistogram(h1_B_Mass_psi2s); + DrawInDefaultCanvas(roofit_hist_jpsi, analysis_name, 0.1); + DrawInDefaultCanvas(roofit_hist_psi2s, analysis_name, 0.1); + + return 0; +} + +// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Hlt2RD_BuToKpMuMu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +int analyze_Hlt2RD_BuToKpMuMu() +{ + const char *analysis_name = "Hlt2RD_BuToKpMuMu"; + const char *end_state_mass_literal = "m(K^{+} #mu^{+}#mu^{-})"; + + TChain *data_chain = new TChain(TString::Format("%s/DecayTree", analysis_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 = FourVect::Init(data_chain, "muplus"); + FourVect *l24v = FourVect::Init(data_chain, "muminus"); + FourVect *hp4v = FourVect::Init(data_chain, "Kplus"); + + 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); + 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); + + for (const auto &var : Hlt1Decisions) + { + h2_Hlt1_flags_B_Mass->Fill(0., var.name.c_str(), 0.); + h2_Hlt1_flags_excl_B_Mass->Fill(0., var.name.c_str(), 0.); + } + + std::map exclusive_hits{}; + + unsigned int entries = data_chain->GetEntries(); + for (unsigned int i = 0; i < entries; i++) + { + data_chain->GetEntry(i); + + TLorentzVector dimuon = l14v->LorentzVector() + l24v->LorentzVector(); + + TLorentzVector lv_hp = hp4v->LorentzVector(); + + Double_t reconstructed_B_Mass = (lv_hp + dimuon).M(); + + if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) + { + for (const auto &var : Hlt1Decisions) + { + if (var.value) + { + h2_Hlt1_flags_B_Mass->Fill(reconstructed_B_Mass, var.name.c_str(), 1); + } + } + + bool exclusive = true; + std::string line{}; + for (const auto &var : Hlt1Decisions) + { + if (var.value) + { + if (!line.empty()) + { + exclusive = false; + } + + line = var.name; + } + } + + if (!line.empty() && exclusive) + { + int &hits = exclusive_hits[line]; + if (hits) + { + hits += 1; + } + else + { + hits = 1; + } + h2_Hlt1_flags_excl_B_Mass->Fill(reconstructed_B_Mass, line.c_str(), 1); + } + } + + if (CutHlt1DecisionsOr({"TwoTrackMuonMVA", "TrackMuonMVA"})) + { + 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("Hlt2RD_BuToKpMuMu", 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_jpsi = CreateRooFitHistogram(h1_B_Mass_jpsi); + auto roofit_hist_psi2s = CreateRooFitHistogram(h1_B_Mass_psi2s); + DrawInDefaultCanvas(roofit_hist_jpsi, analysis_name, 0.1); + DrawInDefaultCanvas(roofit_hist_psi2s, analysis_name, 0.1); + + return 0; +} + +// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Hlt2RD_B0ToKpPimMuMu %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +int analyze_Hlt2RD_B0ToKpPimMuMu() +{ + const char *analysis_name = "Hlt2RD_B0ToKpPimMuMu"; + const char *end_state_mass_literal = "m(K^{+} #pi^{-} #mu^{+}#mu^{-})"; + + TChain *data_chain = new TChain(TString::Format("%s/DecayTree", analysis_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 = FourVect::Init(data_chain, "muplus"); + FourVect *l24v = FourVect::Init(data_chain, "muminus"); + FourVect *hp4v = FourVect::Init(data_chain, "Kplus"); + FourVect *hm4v = FourVect::Init(data_chain, "piminus"); + + 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); + 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); + + for (const auto &var : Hlt1Decisions) + { + h2_Hlt1_flags_B_Mass->Fill(0., var.name.c_str(), 0.); + h2_Hlt1_flags_excl_B_Mass->Fill(0., var.name.c_str(), 0.); + } + + std::map exclusive_hits{}; + + unsigned int entries = data_chain->GetEntries(); + for (unsigned int i = 0; i < entries; i++) + { + data_chain->GetEntry(i); + + TLorentzVector dimuon = l14v->LorentzVector() + l24v->LorentzVector(); + + TLorentzVector lv_hp = hp4v->LorentzVector(); + TLorentzVector lv_hm = hm4v->LorentzVector(); + + Double_t reconstructed_Kstar_Mass = (lv_hp + lv_hm).M(); + Double_t reconstructed_B_Mass = (lv_hp + lv_hm + dimuon).M(); + + if (TMath::Abs(dimuon.M() - JPSI_MASS) < 100.) + { + for (const auto &var : Hlt1Decisions) + { + if (var.value) + { + h2_Hlt1_flags_B_Mass->Fill(reconstructed_B_Mass, var.name.c_str(), 1); + } + } + + bool exclusive = true; + std::string line{}; + for (const auto &var : Hlt1Decisions) + { + if (var.value) + { + if (!line.empty()) + { + exclusive = false; + } + + line = var.name; + } + } + + if (!line.empty() && exclusive) + { + int &hits = exclusive_hits[line]; + if (hits) + { + hits += 1; + } + else + { + hits = 1; + } + h2_Hlt1_flags_excl_B_Mass->Fill(reconstructed_B_Mass, line.c_str(), 1); + } + } + + if (CutHlt1DecisionsOr({"TwoTrackMuonMVA", "TrackMuonMVA"})) + { + if (TMath::Abs(reconstructed_Kstar_Mass - KSTAR_MASS) < 50.) + { + 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(analysis_name, 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_jpsi = CreateRooFitHistogram(h1_B_Mass_jpsi); + auto roofit_hist_psi2s = CreateRooFitHistogram(h1_B_Mass_psi2s); + DrawInDefaultCanvas(roofit_hist_jpsi, analysis_name, 0.1); + DrawInDefaultCanvas(roofit_hist_psi2s, analysis_name, 0.1); + + return 0; +} + +int new_analysis() +{ + return analyze_Hlt2RD_B0ToKpPimMuMu(); +} \ No newline at end of file diff --git a/analysis_fullstream_b02hphmmumu.cpp b/old_analysis_macros/analysis_fullstream_b02hphmmumu.cpp similarity index 68% rename from analysis_fullstream_b02hphmmumu.cpp rename to old_analysis_macros/analysis_fullstream_b02hphmmumu.cpp index b51e6d5..979c5cd 100644 --- a/analysis_fullstream_b02hphmmumu.cpp +++ b/old_analysis_macros/analysis_fullstream_b02hphmmumu.cpp @@ -41,7 +41,7 @@ const Double_t MASS_HIST_FIT_MIN = 5100.; const Double_t MASS_HIST_FIT_MAX = 5700.; const char *TITLE = "SpruceRD_B0ToHpHmMuMu (#pi^{-} #rightarrow K^{-})"; -const char *FILE_NAME = "SpruceRD_B0ToHpHmMuMu_Pip2Kp"; +const char *FILE_NAME = "SpruceRD_B0ToHpHmMuMu_Pim2Km"; const char *MASS_LITERAL = "m(#pi^{+}#pi^{-}_{(#rightarrow K^{-})}#mu^{+}#mu^{-})"; const Double_t K_MASS = 493.677; @@ -49,6 +49,12 @@ const Double_t PSI2S_MASS = 3686.093; const Double_t JPSI_MASS = 3096.9; const Double_t KSTAR_MASS = 891.76; +struct EventOccurance { + int occurances; + double initialCalcP; + double initialFileP; +}; + void CreateRooFitAndDraw(TH1D *hist); int analysis_fullstream_b02hphmmumu() @@ -56,7 +62,7 @@ int analysis_fullstream_b02hphmmumu() TChain *data_chain = new TChain("B0ToHpHmMuMu/DecayTree"); 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/spruce_magdown_2023_v0_tuple_90000000_v0r0p6288631.root"); - + Float_t B0_BPVFDCHI2, B0_BPVIPCHI2, L1_BPVIPCHI2, @@ -72,7 +78,8 @@ int analysis_fullstream_b02hphmmumu() Hm_P, Res_PT, Res_BPVIPCHI2, - Res_P; + Res_P, + B0_P; Double_t L1_PID_MU, L2_PID_MU, @@ -96,7 +103,14 @@ int analysis_fullstream_b02hphmmumu() Hlt1TrackMVADecision, Hlt1TwoTrackMVADecision; + UInt_t runNumber; + ULong64_t eventNumber; + + data_chain->SetBranchAddress("RUNNUMBER", &runNumber); + data_chain->SetBranchAddress("EVENTNUMBER", &eventNumber); + data_chain->SetBranchAddress("B0_M", &B0_M); + data_chain->SetBranchAddress("B0_P", &B0_P); data_chain->SetBranchAddress("B0_BPVFDCHI2", &B0_BPVFDCHI2); @@ -164,12 +178,17 @@ int analysis_fullstream_b02hphmmumu() data_chain->SetBranchAddress("Hm_PZ", &Hm_PZ); data_chain->SetBranchAddress("Hm_ENERGY", &Hm_ENERGY); - TH1D *h1_B0_M = new TH1D("h1_B0_M", TITLE, nBins, MASS_HIST_MIN, MASS_HIST_MAX); + TH1D *h1_B0_M_jpsi = new TH1D("h1_B0_M_jpsi", TString::Format("%s (J/#psi)", TITLE), nBins, MASS_HIST_MIN, MASS_HIST_MAX); + TH1D *h1_B0_M_jpsi_dupl = new TH1D("h1_B0_M_jpsi_dupl", TString::Format("%s (J/#psi)", TITLE), nBins, MASS_HIST_MIN, MASS_HIST_MAX); + TH1D *h1_B0_M_psi2s = new TH1D("h1_B0_M_psi2s", TString::Format("%s (#psi(2S))", TITLE), nBins, MASS_HIST_MIN, MASS_HIST_MAX); TH1D *h1_Res_M = new TH1D("h1_Res_M", "K^{*0} Mass", nBins, 0., 3500.); TH1D *h1_Hm_PID_K = new TH1D("h1_Hm_PID_K", "H^{-} PID K", nBins, -10., 10.); TH1D *h1_Hp_PID_K = new TH1D("h1_Hp_PID_K", "H^{+} PID K", nBins, -10., 10.); + std::map, EventOccurance> run_event_num_dict; + unsigned int entries = data_chain->GetEntries(); + std::cout << "#### Total Entries From Files: " << entries << std::endl; unsigned int fitting_entries = 0; for (unsigned int i = 0; i < entries; i++) { @@ -194,11 +213,6 @@ int analysis_fullstream_b02hphmmumu() Double_t reconstructed_Km_Pip_Mass = (Km_4v + Pip_4v).M(); Double_t reconstructed_B0_Km_Pip_Mass = (Km_4v + Pip_4v + l1_4v + l2_4v).M(); - h1_Res_M->Fill(reconstructed_Kp_Pim_Mass); - - h1_Hm_PID_K->Fill(Hm_PID_K); - h1_Hp_PID_K->Fill(Hp_PID_K); - if (( (B0_BPVFDCHI2 > 64) & (B0_CHI2VXNDOF < 9) & (B0_BPVIPCHI2 < 25) & // @@ -208,21 +222,43 @@ int analysis_fullstream_b02hphmmumu() // (Jpsi_CHI2DOF < 9) & (Jpsi_BPVFDCHI2 > 16) & (Jpsi_M < 5500) & // - (Res_PT > 400) & (Res_MAXDOCACHI2 < 36) & (reconstructed_Kp_Pim_Mass < 2600) & (reconstructed_Km_Pip_Mass < 2600) & (Res_CHI2DOF < 25) & (Res_BPVIPCHI2 > 4) & + (Res_PT > 400) & (Res_MAXDOCACHI2 < 36) & (reconstructed_Km_Pip_Mass < 2600) & (Res_CHI2DOF < 25) & (Res_BPVIPCHI2 > 4) & // (Hm_BPVIPCHI2 > 6) & (Hm_PT > 250) & (Hm_P > 2000) & (Hm_PID_K > 0) & // - (Hp_BPVIPCHI2 > 6) & (Hp_PT > 250) & (Hp_P > 2000) & (Hp_PID_K < 0) - ) & + (Hp_BPVIPCHI2 > 6) & (Hp_PT > 250) & (Hp_P > 2000) & (Hp_PID_K < 0)) & // (Hlt2RD_B0ToKpPimMuMuDecision) & (((Hlt2_InclDetDiMuon_4BodyDecision) | (Hlt2_InclDetDiMuon_3BodyDecision) | (Hlt2_InclDetDiMuonDecision)))) { - if ((TMath::Abs(Jpsi_M - JPSI_MASS) < 100) & ((Hlt1TrackMVADecision) | (Hlt1TwoTrackMVADecision))) + if (TMath::Abs(Jpsi_M - 3096.9) < 100 && ((Hlt1TrackMVADecision) | (Hlt1TwoTrackMVADecision)) + && MASS_HIST_FIT_MIN <= reconstructed_B0_Km_Pip_Mass && reconstructed_B0_Km_Pip_Mass <= MASS_HIST_FIT_MAX) { - h1_B0_M->Fill(reconstructed_B0_Km_Pip_Mass); + auto calcP = (Km_4v.Vect() + Pip_4v.Vect() + l1_4v.Vect() + l2_4v.Vect()).Mag(); + // check for unique run/event number + auto run_event_pair = std::make_pair((unsigned int)runNumber, (unsigned long)eventNumber); + EventOccurance &run_event_unique_num = run_event_num_dict[run_event_pair]; + if (run_event_unique_num.occurances > 0) + { + std::cout << "R: " << runNumber << " / E: " << eventNumber + << ", P Calc = " << calcP << ", P File = " << B0_P + << " Init: P Calc = " << run_event_unique_num.initialCalcP << " P File = " << run_event_unique_num.initialFileP + << " M = " << reconstructed_B0_Km_Pip_Mass + << std::endl; + + run_event_num_dict[run_event_pair] = EventOccurance{run_event_unique_num.occurances + 1, run_event_unique_num.initialCalcP, run_event_unique_num.initialFileP}; + h1_B0_M_jpsi_dupl->Fill(reconstructed_B0_Km_Pip_Mass); + } + else + { + run_event_num_dict[run_event_pair] = EventOccurance{1, calcP, B0_P}; + h1_B0_M_jpsi->Fill(reconstructed_B0_Km_Pip_Mass); + } } } + h1_Hm_PID_K->Fill(Hm_PID_K); + h1_Hp_PID_K->Fill(Hp_PID_K); + if ((i + 1) % 10000 == 0 || i + 1 == entries) { std::cout << "[" @@ -231,10 +267,22 @@ int analysis_fullstream_b02hphmmumu() } } - h1_B0_M->GetXaxis()->SetTitle(MASS_LITERAL); + std::cout << "#### Run/Event Num Dict Size: " << run_event_num_dict.size() << std::endl; + + h1_B0_M_jpsi->GetXaxis()->SetTitle(MASS_LITERAL); + h1_B0_M_psi2s->GetXaxis()->SetTitle(MASS_LITERAL); + h1_B0_M_jpsi_dupl->GetXaxis()->SetTitle(MASS_LITERAL); + h1_B0_M_jpsi_dupl->SetLineColor(kRed); + + - TCanvas *c1 = new TCanvas("c1", "c1", 0, 0, 800, 600); - h1_B0_M->Draw(); + TCanvas *c1 = new TCanvas("c1", "c1", 0, 0, 1600, 600); + c1->Divide(2, 1); + c1->cd(1); + h1_B0_M_jpsi->Draw(); + h1_B0_M_jpsi_dupl->Draw("SAME"); + c1->cd(2); + h1_B0_M_psi2s->Draw(); c1->Draw(); c1->SaveAs(TString::Format("images/root_hist_%s_bmass.pdf", FILE_NAME).Data()); @@ -255,7 +303,7 @@ int analysis_fullstream_b02hphmmumu() c3->BuildLegend(); c3->Draw(); - CreateRooFitAndDraw(h1_B0_M); + CreateRooFitAndDraw(h1_B0_M_jpsi); return 0; } @@ -301,30 +349,31 @@ void CreateRooFitAndDraw(TH1D *hist) RooFitResult *fitres = sigplusbkg.fitTo(roohist_B0_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range("fitting_range")); - sigplusbkg.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144)); - sigplusbkg.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range("fitting_range"), RooFit::Name("sigplusbkg")); - sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(bkg_exp)), RooFit::LineColor(kBlue - 7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("bkg_exp")); - sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(sig_cb)), RooFit::LineColor(kRed - 7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("sig_cb")); + // sigplusbkg.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144)); + // sigplusbkg.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range("fitting_range"), RooFit::Name("sigplusbkg")); + // sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(bkg_exp)), RooFit::LineColor(kBlue - 7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("bkg_exp")); + // sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(sig_cb)), RooFit::LineColor(kRed - 7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("sig_cb")); TCanvas *c = new TCanvas("roofit_c", "roofit_c", 0, 0, 800, 600); roo_frame_mass->Draw(); - TLegend *leg1 = new TLegend(0.50, 0.58, 0.87, 0.89); + // TLegend *leg1 = new TLegend(0.50, 0.58, 0.87, 0.89); + TLegend *leg1 = new TLegend(0.50, 0.80, 0.87, 0.89); leg1->SetFillColor(kWhite); leg1->SetLineColor(kBlack); - leg1->AddEntry(roo_frame_mass->findObject("sigplusbkg"), "Signal + Background", "LP"); - leg1->AddEntry(roo_frame_mass->findObject("sig_cb"), "Signal", "LP"); - leg1->AddEntry(roo_frame_mass->findObject("bkg_exp"), "Background", "LP"); - leg1->AddEntry((TObject *)0, "", ""); - leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_x0.getTitle().Data(), var_mass_x0.getVal(), var_mass_x0.getError()).Data(), ""); - leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_sigmaLR.getTitle().Data(), var_mass_sigmaLR.getVal(), var_mass_sigmaLR.getError()).Data(), ""); - leg1->AddEntry((TObject *)0, TString::Format("%s = %.6f #pm %.6f", var_mass_bkg_c.getTitle().Data(), var_mass_bkg_c.getVal(), var_mass_bkg_c.getError()).Data(), ""); - leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "S", var_mass_nsig.getVal(), var_mass_nsig.getError()).Data(), ""); - leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "B", var_mass_nbkg.getVal(), var_mass_nbkg.getError()).Data(), ""); - - leg1->AddEntry((TObject *)0, TString::Format("Entries: %.1f", var_mass_nsig.getVal() + var_mass_nbkg.getVal()), ""); + // leg1->AddEntry(roo_frame_mass->findObject("sigplusbkg"), "Signal + Background", "LP"); + // leg1->AddEntry(roo_frame_mass->findObject("sig_cb"), "Signal", "LP"); + // leg1->AddEntry(roo_frame_mass->findObject("bkg_exp"), "Background", "LP"); + // leg1->AddEntry((TObject *)0, "", ""); + // leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_x0.getTitle().Data(), var_mass_x0.getVal(), var_mass_x0.getError()).Data(), ""); + // leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_sigmaLR.getTitle().Data(), var_mass_sigmaLR.getVal(), var_mass_sigmaLR.getError()).Data(), ""); + // leg1->AddEntry((TObject *)0, TString::Format("%s = %.6f #pm %.6f", var_mass_bkg_c.getTitle().Data(), var_mass_bkg_c.getVal(), var_mass_bkg_c.getError()).Data(), ""); + // leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "S", var_mass_nsig.getVal(), var_mass_nsig.getError()).Data(), ""); + // leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "B", var_mass_nbkg.getVal(), var_mass_nbkg.getError()).Data(), ""); + + leg1->AddEntry((TObject *)0, TString::Format("Entries: %.0f", roohist_B0_M.sumEntries()), ""); leg1->Draw(); diff --git a/analysis_fullstream_bu2hpmumu.cpp b/old_analysis_macros/analysis_fullstream_bu2hpmumu.cpp similarity index 58% rename from analysis_fullstream_bu2hpmumu.cpp rename to old_analysis_macros/analysis_fullstream_bu2hpmumu.cpp index 748004d..0d6f729 100644 --- a/analysis_fullstream_bu2hpmumu.cpp +++ b/old_analysis_macros/analysis_fullstream_bu2hpmumu.cpp @@ -46,8 +46,8 @@ const Double_t KSTAR_MASS = 891.760; void CreateRooFitAndDraw(TH1D *hist, const char *name); -const char *TITLE = "SpruceRD_BuToHpMuMu (#pi^{+} #rightarrow K^{+}, w/o excl Hlt2 decision cut)"; -const char *FILE_NAME = "SpruceRD_BuToHpMuMu_Pip2Kp_wohlt2cut"; +const char *TITLE = "SpruceRD_BuToHpMuMu (#pi^{+} #rightarrow K^{+}, Small Cut)"; +const char *FILE_NAME = "SpruceRD_BuToHpMuMu_Pip2Kp_smallcut"; const char *MASS_LITERAL = "m(#pi^{+}_{(#rightarrow K^{+})}#mu^{+}#mu^{-})"; struct Hlt1Decision @@ -113,8 +113,8 @@ int CombineHlt1DecisionFlags() int analysis_fullstream_bu2hpmumu() { TChain *data_chain = new TChain("BuToHpMuMu/DecayTree"); - 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/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"); Float_t B_BPVFDCHI2, B_BPVIPCHI2, @@ -145,6 +145,12 @@ int analysis_fullstream_bu2hpmumu() Hlt1TrackMVADecision, Hlt1TwoTrackMVADecision; + UInt_t runNumber; + ULong64_t eventNumber; + + data_chain->SetBranchAddress("RUNNUMBER", &runNumber); + data_chain->SetBranchAddress("EVENTNUMBER", &eventNumber); + data_chain->SetBranchAddress("B_M", &B_M); data_chain->SetBranchAddress("B_BPVFDCHI2", &B_BPVFDCHI2); @@ -176,8 +182,8 @@ int analysis_fullstream_bu2hpmumu() data_chain->SetBranchAddress("Hlt2_InclDetDiMuon_4BodyDecision", &Hlt2_InclDetDiMuon_4BodyDecision); data_chain->SetBranchAddress("Hlt2_InclDetDiMuon_3BodyDecision", &Hlt2_InclDetDiMuon_3BodyDecision); data_chain->SetBranchAddress("Hlt2_InclDetDiMuonDecision", &Hlt2_InclDetDiMuonDecision); - // data_chain->SetBranchAddress("Hlt1TrackMVADecision", &Hlt1TrackMVADecision); - // data_chain->SetBranchAddress("Hlt1TwoTrackMVADecision", &Hlt1TwoTrackMVADecision); + data_chain->SetBranchAddress("Hlt1TrackMVADecision", &Hlt1TrackMVADecision); + data_chain->SetBranchAddress("Hlt1TwoTrackMVADecision", &Hlt1TwoTrackMVADecision); Float_t L1_PX, L1_PY, L1_PZ, L1_ENERGY, L2_PX, L2_PY, L2_PZ, L2_ENERGY, Hp_PX, Hp_PY, Hp_PZ, Hp_ENERGY; data_chain->SetBranchAddress("L1_PX", &L1_PX); @@ -193,20 +199,25 @@ int analysis_fullstream_bu2hpmumu() data_chain->SetBranchAddress("Hp_PZ", &Hp_PZ); data_chain->SetBranchAddress("Hp_ENERGY", &Hp_ENERGY); - ConnectHlt1Decisions(data_chain); + // ConnectHlt1Decisions(data_chain); TH1D *h1_B_M_jpsi = new TH1D("h1_B_M_jpsi", "B Mass, J/#psi Mode", nBins, MASS_HIST_MIN, MASS_HIST_MAX); TH1D *h1_B_M_psi2s = new TH1D("h1_B_M_psi2s", "B Mass, #psi(2S) Mode", nBins, MASS_HIST_MIN, MASS_HIST_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.); - h2_Hlt1_flags_B_Mass->SetCanExtend(TH1::kYaxis); + TH1D *h1_JPsi_M = new TH1D("h1_JPsi_M", "J/#psi", nBins, JPSI_MASS - 100., PSI2S_MASS + 100.); + // 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.); + // h2_Hlt1_flags_B_Mass->SetCanExtend(TH1::kYaxis); - for (const auto &var : Hlt1Decisions) - { - h2_Hlt1_flags_B_Mass->Fill(0., var.name.c_str(), 0.); - } + // for (const auto &var : Hlt1Decisions) + // { + // h2_Hlt1_flags_B_Mass->Fill(0., var.name.c_str(), 0.); + // h2_Hlt1_flags_excl_B_Mass->Fill(0., var.name.c_str(), 0.); + // } + + std::map, int> run_event_num_dict; unsigned int entries = data_chain->GetEntries(); - unsigned int fitting_entries = 0; + unsigned int selected_entries = 0; for (unsigned int i = 0; i < entries; i++) { data_chain->GetEntry(i); @@ -218,44 +229,66 @@ int analysis_fullstream_bu2hpmumu() TLorentzVector l2_4v(L2_PX, L2_PY, L2_PZ, L2_ENERGY); Double_t reconstructed_B_Mass = (K_4v + l1_4v + l2_4v).M(); - int combinedFlags = CombineHlt1DecisionFlags(); - bool trackMVADec = (((combinedFlags & FLAG_TrackMVA) == FLAG_TrackMVA) | ((combinedFlags & FLAG_TwoTrackMVA) == FLAG_TwoTrackMVA)); - if (TMath::Abs(Jpsi_M - JPSI_MASS) < 100 && trackMVADec) - { - for (const auto &var : Hlt1Decisions) - { - if (var.value) - { - h2_Hlt1_flags_B_Mass->Fill(reconstructed_B_Mass, var.name.c_str(), 1); - } - } - } - // std::cout << "HLT1 Decision Flags: " << combinedFlags << std::endl; - - if (trackMVADec) - { - if (TMath::Abs(Jpsi_M - JPSI_MASS) < 100) - { - h1_B_M_jpsi->Fill(reconstructed_B_Mass); - } - else if (TMath::Abs(Jpsi_M - PSI2S_MASS) < 100) - { - h1_B_M_psi2s->Fill(reconstructed_B_Mass); - } - } - - // if ((((B_BPVFDCHI2 > 36) & (B_CHI2VXNDOF < 16) & (B_BPVIPCHI2 < 25) & (Jpsi_MAXDOCACHI2 < 36) & (L1_BPVIPCHI2 > 9) & (L2_BPVIPCHI2 > 9) & (L1_PID_MU > -3) & (L2_PID_MU > -3) & (L1_ISMUON == 1) & (L2_ISMUON == 1) & (L1_PT > 350) & (L2_PT > 350) & (Jpsi_CHI2DOF < 9) & (Jpsi_BPVFDCHI2 > 16) & (Jpsi_M < 5500) & (Hp_PT > 400) & (Hp_BPVIPCHI2 > 6) & (Hp_PT > 400) & (Hp_BPVIPCHI2 > 6) & (Hp_P > 2000) & (Hp_PID_K > -4)) /*& (Hlt2RD_BuToKpMuMuDecision == 1)*/ & (((Hlt2_InclDetDiMuon_4BodyDecision == 1) | (Hlt2_InclDetDiMuon_3BodyDecision == 1) | (Hlt2_InclDetDiMuonDecision == 1))))) + // int combinedFlags = CombineHlt1DecisionFlags(); + // bool trackMVADec = (((combinedFlags & FLAG_TrackMVA) == FLAG_TrackMVA) | ((combinedFlags & FLAG_TwoTrackMVA) == FLAG_TwoTrackMVA)); + // if (TMath::Abs(Jpsi_M - JPSI_MASS) < 100 && trackMVADec) // { - // if ((TMath::Abs(Jpsi_M - 3096.9) < 100) & ((Hlt1TrackMVADecision == 1) | (Hlt1TwoTrackMVADecision == 1))) + // for (const auto &var : Hlt1Decisions) // { - // h1_B_M->Fill(reconstructed_B_Mass); + // if (var.value) + // { + // h2_Hlt1_flags_B_Mass->Fill(reconstructed_B_Mass, var.name.c_str(), 1); + // } + // } - // if (MASS_HIST_FIT_MIN <= reconstructed_B_Mass && reconstructed_B_Mass <= MASS_HIST_FIT_MAX) + // for (const auto &var : Hlt1Decisions) + // { + // if (var.value) // { - // fitting_entries++; + // h2_Hlt1_flags_excl_B_Mass->Fill(reconstructed_B_Mass, var.name.c_str(), 1); + // break; // } // } // } + // std::cout << "HLT1 Decision Flags: " << combinedFlags << std::endl; + + // if (trackMVADec) + // { + // h1_JPsi_M->Fill(Jpsi_M); + // if (TMath::Abs(Jpsi_M - JPSI_MASS) < 100) + // { + // h1_B_M_jpsi->Fill(reconstructed_B_Mass); + // } + // else if (TMath::Abs(Jpsi_M - PSI2S_MASS) < 100) + // { + // h1_B_M_psi2s->Fill(reconstructed_B_Mass); + // } + // } + + // if ( + // (((B_BPVFDCHI2 > 36) & (B_CHI2VXNDOF < 16) & (B_BPVIPCHI2 < 25) & (Jpsi_MAXDOCACHI2 < 36) & (L1_BPVIPCHI2 > 9) & (L2_BPVIPCHI2 > 9) & (L1_PID_MU > -3) & (L2_PID_MU > -3) & (L1_ISMUON) & (L2_ISMUON) & (L1_PT > 350) & (L2_PT > 350) & (Jpsi_CHI2DOF < 9) & (Jpsi_BPVFDCHI2 > 16) & (Jpsi_M < 5500) & (Hp_PT > 400) & (Hp_BPVIPCHI2 > 6) & (Hp_PT > 400) & (Hp_BPVIPCHI2 > 6) & (Hp_P > 2000) & (Hp_PID_K > -4)) & (Hlt2RD_BuToKpMuMuDecision) & (((Hlt2_InclDetDiMuon_4BodyDecision) | (Hlt2_InclDetDiMuon_3BodyDecision) | (Hlt2_InclDetDiMuonDecision))))) + // { + if ((TMath::Abs(Jpsi_M - 3096.9) < 100) & ((Hlt1TrackMVADecision) | (Hlt1TwoTrackMVADecision))) + { + // check for unique run/event number + selected_entries++; + auto run_event_pair = std::make_pair((unsigned int)runNumber, (unsigned long)eventNumber); + int &run_event_unique_num = run_event_num_dict[run_event_pair]; + if (run_event_unique_num) + { + run_event_unique_num = run_event_unique_num + 1; + std::cout << "R: " << runNumber << " / E: " << eventNumber << " / # " << run_event_unique_num + << ", M = " << reconstructed_B_Mass << " MeV" + << std::endl; + } + else + { + run_event_unique_num = 1; + h1_B_M_jpsi->Fill(reconstructed_B_Mass); + } + + } + // } if ((i + 1) % 10000 == 0 || i + 1 == entries) { @@ -265,24 +298,41 @@ int analysis_fullstream_bu2hpmumu() } } - TCanvas *c1 = new TCanvas("c1", "c1", 0, 0, 800, 600); + std::cout << "#### " << run_event_num_dict.size() << " out of " << selected_entries << " unique (" + << TString::Format("%.2f", ((double)(run_event_num_dict.size()) / (double)selected_entries) * 100.) << "%)." << std::endl; + + TCanvas *c1 = new TCanvas("c1", "c1", 0, 0, 1400, 600); + c1->Divide(2, 1); + c1->cd(1); h1_B_M_jpsi->Draw(); + c1->cd(2); + h1_B_M_psi2s->Draw(); c1->Draw(); - c1->SaveAs(TString::Format("images/root_hist_%s_%s_bmass.pdf", FILE_NAME, "jpsi").Data()); + c1->SaveAs(TString::Format("images/root_hist_%s_%s_bmass.pdf", FILE_NAME, "jpsi_psi2s").Data()); TCanvas *c2 = new TCanvas("c2", "c2", 0, 0, 800, 600); - h1_B_M_psi2s->Draw(); + h1_JPsi_M->SetStats(0); + h1_JPsi_M->Draw(); c2->Draw(); - c2->SaveAs(TString::Format("images/root_hist_%s_%s_bmass.pdf", FILE_NAME, "psi2s").Data()); - - TCanvas *c3 = new TCanvas("c3", "c3", 0, 0, 800, 600); - c3->SetLeftMargin(0.20); - h2_Hlt1_flags_B_Mass->SetStats(0); - h2_Hlt1_flags_B_Mass->Draw("COLZ"); - c3->Draw(); - // CreateRooFitAndDraw(h1_B_M_jpsi, "jpsi"); + // TCanvas *c3 = new TCanvas("c3", "c3", 0, 0, 1400, 600); + // c3->Divide(2, 1, 0.04); + // c3->cd(1); + // //c3->SetLeftMargin(0.40); + // h2_Hlt1_flags_B_Mass->SetStats(0); + // h2_Hlt1_flags_B_Mass->Draw("COLZ"); + // c3->cd(2); + // //c3->SetLeftMargin(0.40); + // h2_Hlt1_flags_excl_B_Mass->SetStats(0); + // h2_Hlt1_flags_excl_B_Mass->Draw("COLZ"); + // c3->Draw(); + // c3->SaveAs(TString::Format("images/root_hist_%s_hlt1dec.pdf", FILE_NAME).Data()); + + CreateRooFitAndDraw(h1_B_M_jpsi, "jpsi"); // CreateRooFitAndDraw(h1_B_M_psi2s, "psi2s"); + // zweiter Plot, welche exklusiv von *einer* HLT1 Line getriggert werden + // DecProd Cut in Decfile von MonteCarlo Production + // --> Lumi * (B Prod Crosssection für B-Prouktion) * (Branching Frcation B->J/Psi K) * (Brranch Fraction J/Psi->mumu) * Effizienz return 0; } @@ -291,11 +341,9 @@ void CreateRooFitAndDraw(TH1D *hist, const char *name) { auto suffix = [name](const char *text) { - return TString::Format("%s_%s", text, name); + return TString::Format("%s%s", text, ""); }; - std::cout << " ###### " << suffix("AABB") << std::endl; - RooRealVar roo_var_mass(suffix("var_mass"), "B Mass Variable", MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX); roo_var_mass.setRange("fitting_range", MASS_HIST_FIT_MIN, MASS_HIST_FIT_MAX); @@ -334,52 +382,52 @@ void CreateRooFitAndDraw(TH1D *hist, const char *name) TString pdf_name = suffix("sigplusbkg"); RooAddPdf sigplusbkg(pdf_name, "Sig and Bkg PDF", RooArgList(sig_cb, bkg_exp), RooArgList(var_mass_nsig, var_mass_nbkg)); - RooFitResult *fitres = sigplusbkg.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range("fitting_range")); + // RooFitResult *fitres = sigplusbkg.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range("fitting_range")); // sigplusbkg.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144)); - sigplusbkg.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range("fitting_range"), RooFit::Name(pdf_name)); - sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(bkg_exp)), RooFit::LineColor(kBlue - 7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("bkg_exp")); - sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(sig_cb)), RooFit::LineColor(kRed - 7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("sig_cb")); + // sigplusbkg.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range("fitting_range"), RooFit::Name(pdf_name)); + // sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(bkg_exp)), RooFit::LineColor(kBlue - 7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("bkg_exp")); + // sigplusbkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(sig_cb)), RooFit::LineColor(kRed - 7), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range"), RooFit::Name("sig_cb")); - RooPlot *roo_frame_pull = roo_var_mass.frame(RooFit::Title("Pull Distribution")); - roo_frame_pull->addPlotable(roo_frame_mass->pullHist(hist_name, pdf_name), "P"); + // RooPlot *roo_frame_pull = roo_var_mass.frame(RooFit::Title("Pull Distribution")); + // roo_frame_pull->addPlotable(roo_frame_mass->pullHist(hist_name, pdf_name), "P"); TCanvas *c = new TCanvas(suffix("roofit_c"), "roofit_c", 0, 0, 800, 600); - auto *p2 = new TPad(suffix("p2"), "Pull", 0., 0., 1., 0.3); - p2->Draw(); - p2->SetTopMargin(0.001); - p2->SetBottomMargin(0.3); - p2->SetGrid(); + // auto *p2 = new TPad(suffix("p2"), "Pull", 0., 0., 1., 0.3); + // p2->Draw(); + // p2->SetTopMargin(0.001); + // p2->SetBottomMargin(0.3); + // p2->SetGrid(); - auto *p1 = new TPad(suffix("p1"), "Fit", 0., 0.32, 1., 1.); - p1->Draw(); - p1->SetBottomMargin(0.001); - p1->cd(); + // auto *p1 = new TPad(suffix("p1"), "Fit", 0., 0.32, 1., 1.); + // p1->Draw(); + // p1->SetBottomMargin(0.001); + // p1->cd(); roo_frame_mass->Draw(); - TLegend *leg1 = new TLegend(0.50, 0.50, 0.87, 0.89); + TLegend *leg1 = new TLegend(0.50, 0.80, 0.87, 0.89); leg1->SetFillColor(kWhite); leg1->SetLineColor(kBlack); - leg1->AddEntry(roo_frame_mass->findObject(pdf_name), "Signal + Background", "LP"); - leg1->AddEntry(roo_frame_mass->findObject("sig_cb"), "Signal", "LP"); - leg1->AddEntry(roo_frame_mass->findObject("bkg_exp"), "Background", "LP"); - leg1->AddEntry((TObject *)0, "", ""); - leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_x0.getTitle().Data(), var_mass_x0.getVal(), var_mass_x0.getError()).Data(), ""); - leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_sigmaLR.getTitle().Data(), var_mass_sigmaLR.getVal(), var_mass_sigmaLR.getError()).Data(), ""); - leg1->AddEntry((TObject *)0, TString::Format("%s = %.6f #pm %.6f", var_mass_bkg_c.getTitle().Data(), var_mass_bkg_c.getVal(), var_mass_bkg_c.getError()).Data(), ""); - leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "S", var_mass_nsig.getVal(), var_mass_nsig.getError()).Data(), ""); - leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "B", var_mass_nbkg.getVal(), var_mass_nbkg.getError()).Data(), ""); + // leg1->AddEntry(roo_frame_mass->findObject(pdf_name), "Signal + Background", "LP"); + // leg1->AddEntry(roo_frame_mass->findObject("sig_cb"), "Signal", "LP"); + // leg1->AddEntry(roo_frame_mass->findObject("bkg_exp"), "Background", "LP"); + // leg1->AddEntry((TObject *)0, "", ""); + // leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_x0.getTitle().Data(), var_mass_x0.getVal(), var_mass_x0.getError()).Data(), ""); + // leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_sigmaLR.getTitle().Data(), var_mass_sigmaLR.getVal(), var_mass_sigmaLR.getError()).Data(), ""); + // leg1->AddEntry((TObject *)0, TString::Format("%s = %.6f #pm %.6f", var_mass_bkg_c.getTitle().Data(), var_mass_bkg_c.getVal(), var_mass_bkg_c.getError()).Data(), ""); + // leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "S", var_mass_nsig.getVal(), var_mass_nsig.getError()).Data(), ""); + // leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "B", var_mass_nbkg.getVal(), var_mass_nbkg.getError()).Data(), ""); - leg1->AddEntry((TObject *)0, TString::Format("Entries: %.1f", var_mass_nsig.getVal() + var_mass_nbkg.getVal()), ""); + leg1->AddEntry((TObject *)0, TString::Format("Entries: %.0f", roohist_B_M.sumEntries()), ""); leg1->Draw(); - p2->cd(); + // p2->cd(); - roo_frame_pull->Draw(); + // roo_frame_pull->Draw(); c->Draw(); c->SaveAs(TString::Format("images/data/roofit_hist_%s_bmass.pdf", suffix(FILE_NAME).Data()).Data()); diff --git a/analysis_turbo.cpp b/old_analysis_macros/analysis_turbo.cpp similarity index 57% rename from analysis_turbo.cpp rename to old_analysis_macros/analysis_turbo.cpp index 0afbf56..d8b6151 100644 --- a/analysis_turbo.cpp +++ b/old_analysis_macros/analysis_turbo.cpp @@ -37,12 +37,15 @@ const int nBins = 70; const Double_t MASS_HIST_MIN = 5150.; const Double_t MASS_HIST_MAX = 6000.; +const Double_t PSI2S_MASS = 3686.093; +const Double_t JPSI_MASS = 3096.9; + const char *B0_DECAY = "Hlt2RD_B0ToKpPimMuMu"; const char *Bu_DECAY = "Hlt2RD_BuToKpMuMu"; -const char *TITLE = B0_DECAY; -const char *HIST_TITLE = "Hlt2RD_BuToKpMuMu"; -const char *X_AXIS = "m(K^{+}#mu^{+}#mu^{-})"; +const char *TITLE = Bu_DECAY; +const char *HIST_TITLE = Bu_DECAY; // "Hlt2RD_BuToKpMuMu"; +const char *X_AXIS = "m(K^{+}#mu^{+}#mu^{-})"; // "m(K^{+}#mu^{+}#mu^{-})"; void CreateRooFitAndDraw(TH1D *hist, int fitting_entries); @@ -64,6 +67,12 @@ int analysis_turbo() Bool_t Hlt1TrackMVADecision, Hlt1TwoTrackMVADecision; + UInt_t runNumber; + ULong64_t eventNumber; + + data_chain->SetBranchAddress("RUNNUMBER", &runNumber); + data_chain->SetBranchAddress("EVENTNUMBER", &eventNumber); + if (TITLE == B0_DECAY) { data_chain->SetBranchAddress("B0_M", &B0_M); @@ -78,34 +87,77 @@ int analysis_turbo() data_chain->SetBranchAddress("Hlt1TrackMVADecision", &Hlt1TrackMVADecision); data_chain->SetBranchAddress("Hlt1TwoTrackMVADecision", &Hlt1TwoTrackMVADecision); - TH1D *h1_B_M = new TH1D("h1_B_M", HIST_TITLE, nBins, MASS_HIST_MIN, MASS_HIST_MAX); - h1_B_M->GetXaxis()->SetTitle(X_AXIS); + TH1D *h1_B_M_jpsi = new TH1D("h1_B_M_jpsi", TString::Format("%s (J/#psi)", HIST_TITLE), nBins, MASS_HIST_MIN, MASS_HIST_MAX); + TH1D *h1_B_M_psi2s = new TH1D("h1_B_M_psi2s", TString::Format("%s (#psi(2S))", HIST_TITLE), nBins, MASS_HIST_MIN, MASS_HIST_MAX); + h1_B_M_jpsi->GetXaxis()->SetTitle(X_AXIS); + h1_B_M_psi2s->GetXaxis()->SetTitle(X_AXIS); + + std::map, int> run_event_num_dict; int fitting_entries = 0; unsigned int entries = data_chain->GetEntries(); + unsigned int selected_entries = 0; for (unsigned int i = 0; i < entries; i++) { data_chain->GetEntry(i); if (TITLE == B0_DECAY) { - if ((TMath::Abs(Jpsi_M - 3096.9) < 100) && (B0_M > 4500) && (B0_M < 6000) && (TMath::Abs(Kst_M - 895.55) < 50) && ((Hlt1TrackMVADecision) || (Hlt1TwoTrackMVADecision))) + if ((B0_M > 4500) && (B0_M < 6000) && (TMath::Abs(Kst_M - 895.55) < 50) && ((Hlt1TrackMVADecision) || (Hlt1TwoTrackMVADecision))) { - h1_B_M->Fill(B0_M); - if (MASS_HIST_MIN <= B0_M && B0_M <= MASS_HIST_MAX) + // check for unique run/event number + selected_entries++; + auto run_event_pair = std::make_pair((unsigned int)runNumber, (unsigned long)eventNumber); + int &run_event_unique_num = run_event_num_dict[run_event_pair]; + if (run_event_unique_num) + { + run_event_unique_num = run_event_unique_num + 1; + std::cout << "R: " << runNumber << " / E: " << eventNumber << " / # " << run_event_unique_num + << ", M = " << B0_M << " MeV" + << std::endl; + continue; + } + else + { + run_event_unique_num = 1; + } + if (TMath::Abs(Jpsi_M - JPSI_MASS) < 100) + { + h1_B_M_jpsi->Fill(B0_M); + } + else if (TMath::Abs(Jpsi_M - PSI2S_MASS) < 100) { - fitting_entries++; + h1_B_M_psi2s->Fill(B0_M); } } } else { - if ((TMath::Abs(Jpsi_M - 3096.9) < 100.) && (B_M > 4500.) && (B_M < 6000.) && ((Hlt1TrackMVADecision) || (Hlt1TwoTrackMVADecision))) + if ((B_M > 4500.) && (B_M < 6000.) && ((Hlt1TrackMVADecision) || (Hlt1TwoTrackMVADecision))) { - h1_B_M->Fill(B_M); - if (MASS_HIST_MIN <= B_M && B_M <= MASS_HIST_MAX) + // check for unique run/event number + selected_entries++; + auto run_event_pair = std::make_pair((unsigned int)runNumber, (unsigned long)eventNumber); + int &run_event_unique_num = run_event_num_dict[run_event_pair]; + if (run_event_unique_num) + { + run_event_unique_num = run_event_unique_num + 1; + std::cout << "R: " << runNumber << " / E: " << eventNumber << " / # " << run_event_unique_num + << ", M = " << B_M << " MeV" + << std::endl; + continue; + } + else + { + run_event_unique_num = 1; + } + if (TMath::Abs(Jpsi_M - JPSI_MASS) < 100) { - fitting_entries++; + h1_B_M_jpsi->Fill(B_M); + } + else if (TMath::Abs(Jpsi_M - PSI2S_MASS) < 100) + { + h1_B_M_psi2s->Fill(B_M); } } } @@ -118,12 +170,19 @@ int analysis_turbo() } } - TCanvas *c1 = new TCanvas("c1", "c1", 0, 0, 800, 600); - h1_B_M->Draw(); + std::cout << "#### " << run_event_num_dict.size() << " out of " << selected_entries << " unique (" + << TString::Format("%.2f", ((double)(run_event_num_dict.size()) / (double)selected_entries) * 100.) << "%)." << std::endl; + + TCanvas *c1 = new TCanvas("c1", "c1", 0, 0, 1400, 600); + c1->Divide(2, 1); + c1->cd(1); + h1_B_M_jpsi->Draw(); + c1->cd(2); + h1_B_M_psi2s->Draw(); c1->Draw(); c1->SaveAs(TString::Format("images/root_hist_%s_bmass.pdf", TITLE).Data()); - CreateRooFitAndDraw(h1_B_M, fitting_entries); + CreateRooFitAndDraw(h1_B_M_jpsi, fitting_entries); return 0; } @@ -174,7 +233,7 @@ void CreateRooFitAndDraw(TH1D *hist, int fitting_entries) // TString pdf_name = TString::Format("%s_sigplusbkg", var_id.c_str()); RooAddPdf sigplusbkg("sigplusbkg", "Sig and Bkg PDF", RooArgList(sig_cb, bkg_exp), RooArgList(var_mass_nsig, var_mass_nbkg)); - RooFitResult *fitres = sigplusbkg.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range("fitting_range")); + // RooFitResult *fitres = sigplusbkg.fitTo(roohist_B_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range("fitting_range")); // sigplusbkg.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kYellow -7), RooFit::FillStyle(3144)); // sigplusbkg.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range("fitting_range"), RooFit::Name("sigplusbkg")); @@ -193,13 +252,14 @@ void CreateRooFitAndDraw(TH1D *hist, int fitting_entries) // leg1->AddEntry(roo_frame_mass->findObject("sig_cb"), "Signal", "LP"); // leg1->AddEntry(roo_frame_mass->findObject("bkg_exp"), "Background", "LP"); // leg1->AddEntry((TObject *)0, "", ""); - leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_x0.getTitle().Data(), var_mass_x0.getVal(), var_mass_x0.getError()).Data(), ""); - leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_sigmaLR.getTitle().Data(), var_mass_sigmaLR.getVal(), var_mass_sigmaLR.getError()).Data(), ""); - leg1->AddEntry((TObject *)0, TString::Format("%s = %.6f #pm %.6f", var_mass_bkg_c.getTitle().Data(), var_mass_bkg_c.getVal(), var_mass_bkg_c.getError()).Data(), ""); - leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "S", var_mass_nsig.getVal(), var_mass_nsig.getError()).Data(), ""); - leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "B", var_mass_nbkg.getVal(), var_mass_nbkg.getError()).Data(), ""); - - leg1->AddEntry((TObject *)0, TString::Format("Entries: %d", fitting_entries), ""); + // leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_x0.getTitle().Data(), var_mass_x0.getVal(), var_mass_x0.getError()).Data(), ""); + // leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", var_mass_sigmaLR.getTitle().Data(), var_mass_sigmaLR.getVal(), var_mass_sigmaLR.getError()).Data(), ""); + // leg1->AddEntry((TObject *)0, TString::Format("%s = %.6f #pm %.6f", var_mass_bkg_c.getTitle().Data(), var_mass_bkg_c.getVal(), var_mass_bkg_c.getError()).Data(), ""); + // leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "S", var_mass_nsig.getVal(), var_mass_nsig.getError()).Data(), ""); + // leg1->AddEntry((TObject *)0, TString::Format("%s = %.3f #pm %.3f", "B", var_mass_nbkg.getVal(), var_mass_nbkg.getError()).Data(), ""); + + // leg1->AddEntry((TObject *)0, TString::Format("Entries: %.2f", var_mass_nsig.getVal() + var_mass_nbkg.getVal()), ""); + leg1->AddEntry((TObject *)0, TString::Format("Entries: %.0f", roohist_B_M.sumEntries()), ""); leg1->Draw(); diff --git a/do_all.cpp b/old_analysis_macros/do_all.cpp similarity index 87% rename from do_all.cpp rename to old_analysis_macros/do_all.cpp index 339e5f4..237fc03 100644 --- a/do_all.cpp +++ b/old_analysis_macros/do_all.cpp @@ -275,9 +275,13 @@ public: bools[name] = Variable(name, cuts); } - void AddTriggerCut(const std::vector &triggers) + void AddTriggerCut(const std::vector &triggers_list) { - trigger_cuts.push_back(triggers); + for (const auto &trig : triggers_list) + { + bools[trig.c_str()] = Variable(trig.c_str(), {}); + } + trigger_cuts.push_back(triggers_list); } void AddFourMomentumFor(const char *name) @@ -418,7 +422,7 @@ public: } }; -int do_all() +void BuToHpMuMu() { TChain *data_chain = new TChain("BuToHpMuMu/DecayTree"); @@ -426,38 +430,12 @@ int do_all() data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0_tuple_90000000_v0r0p6288631.root"); VariableCollection col{}; - col.AddFloat("B_BPVFDCHI2", {new GT(36.)}); - col.AddDouble("B_CHI2VXNDOF", {new LT(16.)}); - col.AddFloat("B_BPVIPCHI2", {new LT(25.)}); - col.AddDouble("Jpsi_MAXDOCACHI2", {new LT(36.)}); - col.AddFloat("L1_BPVIPCHI2", {new GT(9.)}); - col.AddFloat("L2_BPVIPCHI2", {new GT(9.)}); - col.AddDouble("L1_PID_MU", {new GT(-3.)}); - col.AddDouble("L2_PID_MU", {new GT(-3.)}); - col.AddBool("L1_ISMUON", {new EQ(true)}); - col.AddBool("L2_ISMUON", {new EQ(true)}); - col.AddFloat("L1_PT", {new GT(350.)}); - col.AddFloat("L2_PT", {new GT(350.)}); - col.AddDouble("Jpsi_CHI2DOF", {new LT(9.)}); - col.AddFloat("Jpsi_BPVFDCHI2", {new GT(16.)}); - col.AddDouble("Jpsi_M", {new LT(5500.), new BT(3096.9 - 100., 3096.9 + 100.)}); - col.AddFloat("Hp_PT", {new GT(400.)}); - col.AddFloat("Hp_BPVIPCHI2", {new GT(6.)}); - col.AddFloat("Hp_P", {new GT(2000.)}); - col.AddDouble("Hp_PID_K", {new GT(-4.)}); - col.AddBool("Hlt2_InclDetDiMuon_4BodyDecision", {}); - col.AddBool("Hlt2_InclDetDiMuon_3BodyDecision", {}); - col.AddBool("Hlt2_InclDetDiMuonDecision", {}); - col.AddBool("Hlt2RD_BuToKpMuMuDecision", {}); - col.AddBool("Hlt1TrackMVADecision", {}); - col.AddBool("Hlt1TwoTrackMVADecision", {}); col.AddFourMomentumFor("L1"); col.AddFourMomentumFor("L2"); col.AddFourMomentumFor("Hp"); col.AddTriggerCut({"Hlt2_InclDetDiMuonDecision", "Hlt2_InclDetDiMuon_3BodyDecision", "Hlt2_InclDetDiMuon_4BodyDecision"}); col.AddTriggerCut({"Hlt1TrackMVADecision", "Hlt1TwoTrackMVADecision"}); - col.AddTriggerCut({"Hlt2RD_BuToKpMuMuDecision"}); col.Connect(data_chain); @@ -492,5 +470,10 @@ int do_all() h1_B_M->Draw(); c1->Draw(); + return 0; +} + +int do_all() { + BuToHpMuMu(); return 0; } \ No newline at end of file diff --git a/simulation.cpp b/old_analysis_macros/simulation.cpp similarity index 100% rename from simulation.cpp rename to old_analysis_macros/simulation.cpp diff --git a/simulation_fullstream_B02KpPimMuMu.cpp b/old_analysis_macros/simulation_fullstream_B02KpPimMuMu.cpp similarity index 100% rename from simulation_fullstream_B02KpPimMuMu.cpp rename to old_analysis_macros/simulation_fullstream_B02KpPimMuMu.cpp diff --git a/simulation_fullstream_Bu2KpMuMu.cpp b/old_analysis_macros/simulation_fullstream_Bu2KpMuMu.cpp similarity index 100% rename from simulation_fullstream_Bu2KpMuMu.cpp rename to old_analysis_macros/simulation_fullstream_Bu2KpMuMu.cpp