pfeiffer
/
inclusive_detached_dilepton
1
0
Fork 0
Code Issues Pull Requests 1 Projects Releases Wiki Activity
Browse Source

move old macros, add new one with hlt1 check and psi2s

pull/1/head
Marius Pfeiffer 8 months ago
parent
b1709c1a44
commit
8fd431b64d
9 changed files with 1047 additions and 174 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 1
      .gitignore
  2. 732
      new_analysis.cpp
  3. 115
      old_analysis_macros/analysis_fullstream_b02hphmmumu.cpp
  4. 216
      old_analysis_macros/analysis_fullstream_bu2hpmumu.cpp
  5. 108
      old_analysis_macros/analysis_turbo.cpp
  6. 41
      old_analysis_macros/do_all.cpp
  7. 0
      old_analysis_macros/simulation.cpp
  8. 0
      old_analysis_macros/simulation_fullstream_B02KpPimMuMu.cpp
  9. 0
      old_analysis_macros/simulation_fullstream_Bu2KpMuMu.cpp

1
.gitignore
View File

@ -3,3 +3,4 @@ dataloader/**
status_report/** status_report/**
branchnames.txt branchnames.txt
images/** images/**
output_files/**

732
new_analysis.cpp
View File

@ -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 <string>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <filesystem>
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<Hlt1Decision> 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<std::string> &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<std::string> &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<std::string, int> 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<std::string, int> 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<std::string, int> 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<std::string, int> 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();
}

115
analysis_fullstream_b02hphmmumu.cpp → old_analysis_macros/analysis_fullstream_b02hphmmumu.cpp
View File

@ -41,7 +41,7 @@ const Double_t MASS_HIST_FIT_MIN = 5100.;
const Double_t MASS_HIST_FIT_MAX = 5700.; const Double_t MASS_HIST_FIT_MAX = 5700.;
const char *TITLE = "SpruceRD_B0ToHpHmMuMu (#pi^{-} #rightarrow K^{-})"; 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 char *MASS_LITERAL = "m(#pi^{+}#pi^{-}_{(#rightarrow K^{-})}#mu^{+}#mu^{-})";
const Double_t K_MASS = 493.677; 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 JPSI_MASS = 3096.9;
const Double_t KSTAR_MASS = 891.76; const Double_t KSTAR_MASS = 891.76;
struct EventOccurance {
int occurances;
double initialCalcP;
double initialFileP;
};
void CreateRooFitAndDraw(TH1D *hist); void CreateRooFitAndDraw(TH1D *hist);
int analysis_fullstream_b02hphmmumu() int analysis_fullstream_b02hphmmumu()
@ -72,7 +78,8 @@ int analysis_fullstream_b02hphmmumu()
Hm_P, Hm_P,
Res_PT, Res_PT,
Res_BPVIPCHI2, Res_BPVIPCHI2,
Res_P;
Res_P,
B0_P;
Double_t L1_PID_MU, Double_t L1_PID_MU,
L2_PID_MU, L2_PID_MU,
@ -96,7 +103,14 @@ int analysis_fullstream_b02hphmmumu()
Hlt1TrackMVADecision, Hlt1TrackMVADecision,
Hlt1TwoTrackMVADecision; 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_M", &B0_M);
data_chain->SetBranchAddress("B0_P", &B0_P);
data_chain->SetBranchAddress("B0_BPVFDCHI2", &B0_BPVFDCHI2); 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_PZ", &Hm_PZ);
data_chain->SetBranchAddress("Hm_ENERGY", &Hm_ENERGY); 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_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_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.); TH1D *h1_Hp_PID_K = new TH1D("h1_Hp_PID_K", "H^{+} PID K", nBins, -10., 10.);
std::map<std::pair<unsigned int, unsigned long>, EventOccurance> run_event_num_dict;
unsigned int entries = data_chain->GetEntries(); unsigned int entries = data_chain->GetEntries();
std::cout << "#### Total Entries From Files: " << entries << std::endl;
unsigned int fitting_entries = 0; unsigned int fitting_entries = 0;
for (unsigned int i = 0; i < entries; i++) 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_Km_Pip_Mass = (Km_4v + Pip_4v).M();
Double_t reconstructed_B0_Km_Pip_Mass = (Km_4v + Pip_4v + l1_4v + l2_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 (( if ((
(B0_BPVFDCHI2 > 64) & (B0_CHI2VXNDOF < 9) & (B0_BPVIPCHI2 < 25) & (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) & (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) & (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)))) (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)
{
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)
{ {
h1_B0_M->Fill(reconstructed_B0_Km_Pip_Mass);
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) if ((i + 1) % 10000 == 0 || i + 1 == entries)
{ {
std::cout << "[" 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->Draw();
c1->SaveAs(TString::Format("images/root_hist_%s_bmass.pdf", FILE_NAME).Data()); c1->SaveAs(TString::Format("images/root_hist_%s_bmass.pdf", FILE_NAME).Data());
@ -255,7 +303,7 @@ int analysis_fullstream_b02hphmmumu()
c3->BuildLegend(); c3->BuildLegend();
c3->Draw(); c3->Draw();
CreateRooFitAndDraw(h1_B0_M);
CreateRooFitAndDraw(h1_B0_M_jpsi);
return 0; 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")); 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); TCanvas *c = new TCanvas("roofit_c", "roofit_c", 0, 0, 800, 600);
roo_frame_mass->Draw(); 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->SetFillColor(kWhite);
leg1->SetLineColor(kBlack); 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(); leg1->Draw();

216
analysis_fullstream_bu2hpmumu.cpp → old_analysis_macros/analysis_fullstream_bu2hpmumu.cpp
View File

@ -46,8 +46,8 @@ const Double_t KSTAR_MASS = 891.760;
void CreateRooFitAndDraw(TH1D *hist, const char *name); 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^{-})"; const char *MASS_LITERAL = "m(#pi^{+}_{(#rightarrow K^{+})}#mu^{+}#mu^{-})";
struct Hlt1Decision struct Hlt1Decision
@ -113,8 +113,8 @@ int CombineHlt1DecisionFlags()
int analysis_fullstream_bu2hpmumu() int analysis_fullstream_bu2hpmumu()
{ {
TChain *data_chain = new TChain("BuToHpMuMu/DecayTree"); 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_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, Float_t B_BPVFDCHI2,
B_BPVIPCHI2, B_BPVIPCHI2,
@ -145,6 +145,12 @@ int analysis_fullstream_bu2hpmumu()
Hlt1TrackMVADecision, Hlt1TrackMVADecision,
Hlt1TwoTrackMVADecision; 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_M", &B_M);
data_chain->SetBranchAddress("B_BPVFDCHI2", &B_BPVFDCHI2); 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_4BodyDecision", &Hlt2_InclDetDiMuon_4BodyDecision);
data_chain->SetBranchAddress("Hlt2_InclDetDiMuon_3BodyDecision", &Hlt2_InclDetDiMuon_3BodyDecision); data_chain->SetBranchAddress("Hlt2_InclDetDiMuon_3BodyDecision", &Hlt2_InclDetDiMuon_3BodyDecision);
data_chain->SetBranchAddress("Hlt2_InclDetDiMuonDecision", &Hlt2_InclDetDiMuonDecision); 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; 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); 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_PZ", &Hp_PZ);
data_chain->SetBranchAddress("Hp_ENERGY", &Hp_ENERGY); 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_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); 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<std::pair<unsigned int, unsigned long>, int> run_event_num_dict;
unsigned int entries = data_chain->GetEntries(); unsigned int entries = data_chain->GetEntries();
unsigned int fitting_entries = 0;
unsigned int selected_entries = 0;
for (unsigned int i = 0; i < entries; i++) for (unsigned int i = 0; i < entries; i++)
{ {
data_chain->GetEntry(i); data_chain->GetEntry(i);
@ -218,43 +229,65 @@ int analysis_fullstream_bu2hpmumu()
TLorentzVector l2_4v(L2_PX, L2_PY, L2_PZ, L2_ENERGY); TLorentzVector l2_4v(L2_PX, L2_PY, L2_PZ, L2_ENERGY);
Double_t reconstructed_B_Mass = (K_4v + l1_4v + l2_4v).M(); 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);
}
}
}
// 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);
// }
// }
// for (const auto &var : Hlt1Decisions)
// {
// if (var.value)
// {
// h2_Hlt1_flags_excl_B_Mass->Fill(reconstructed_B_Mass, var.name.c_str(), 1);
// break;
// }
// }
// }
// std::cout << "HLT1 Decision Flags: " << combinedFlags << std::endl; // std::cout << "HLT1 Decision Flags: " << combinedFlags << std::endl;
if (trackMVADec)
// 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)))
{ {
if (TMath::Abs(Jpsi_M - JPSI_MASS) < 100)
// 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)
{ {
h1_B_M_jpsi->Fill(reconstructed_B_Mass);
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 if (TMath::Abs(Jpsi_M - PSI2S_MASS) < 100)
else
{ {
h1_B_M_psi2s->Fill(reconstructed_B_Mass);
}
run_event_unique_num = 1;
h1_B_M_jpsi->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)))))
// {
// if ((TMath::Abs(Jpsi_M - 3096.9) < 100) & ((Hlt1TrackMVADecision == 1) | (Hlt1TwoTrackMVADecision == 1)))
// {
// h1_B_M->Fill(reconstructed_B_Mass);
// if (MASS_HIST_FIT_MIN <= reconstructed_B_Mass && reconstructed_B_Mass <= MASS_HIST_FIT_MAX)
// {
// fitting_entries++;
// }
// }
}
// } // }
if ((i + 1) % 10000 == 0 || i + 1 == entries) 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(); h1_B_M_jpsi->Draw();
c1->cd(2);
h1_B_M_psi2s->Draw();
c1->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); 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->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"); // 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; return 0;
} }
@ -291,11 +341,9 @@ void CreateRooFitAndDraw(TH1D *hist, const char *name)
{ {
auto suffix = [name](const char *text) 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); 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); 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"); 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)); 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::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); 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(); 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->SetFillColor(kWhite);
leg1->SetLineColor(kBlack); 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(); leg1->Draw();
p2->cd();
// p2->cd();
roo_frame_pull->Draw();
// roo_frame_pull->Draw();
c->Draw(); c->Draw();
c->SaveAs(TString::Format("images/data/roofit_hist_%s_bmass.pdf", suffix(FILE_NAME).Data()).Data()); c->SaveAs(TString::Format("images/data/roofit_hist_%s_bmass.pdf", suffix(FILE_NAME).Data()).Data());

108
analysis_turbo.cpp → old_analysis_macros/analysis_turbo.cpp
View File

@ -37,12 +37,15 @@ const int nBins = 70;
const Double_t MASS_HIST_MIN = 5150.; const Double_t MASS_HIST_MIN = 5150.;
const Double_t MASS_HIST_MAX = 6000.; 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 *B0_DECAY = "Hlt2RD_B0ToKpPimMuMu";
const char *Bu_DECAY = "Hlt2RD_BuToKpMuMu"; 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); void CreateRooFitAndDraw(TH1D *hist, int fitting_entries);
@ -64,6 +67,12 @@ int analysis_turbo()
Bool_t Hlt1TrackMVADecision, Hlt1TwoTrackMVADecision; Bool_t Hlt1TrackMVADecision, Hlt1TwoTrackMVADecision;
UInt_t runNumber;
ULong64_t eventNumber;
data_chain->SetBranchAddress("RUNNUMBER", &runNumber);
data_chain->SetBranchAddress("EVENTNUMBER", &eventNumber);
if (TITLE == B0_DECAY) if (TITLE == B0_DECAY)
{ {
data_chain->SetBranchAddress("B0_M", &B0_M); data_chain->SetBranchAddress("B0_M", &B0_M);
@ -78,34 +87,77 @@ int analysis_turbo()
data_chain->SetBranchAddress("Hlt1TrackMVADecision", &Hlt1TrackMVADecision); data_chain->SetBranchAddress("Hlt1TrackMVADecision", &Hlt1TrackMVADecision);
data_chain->SetBranchAddress("Hlt1TwoTrackMVADecision", &Hlt1TwoTrackMVADecision); 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<std::pair<unsigned int, unsigned long>, int> run_event_num_dict;
int fitting_entries = 0; int fitting_entries = 0;
unsigned int entries = data_chain->GetEntries(); unsigned int entries = data_chain->GetEntries();
unsigned int selected_entries = 0;
for (unsigned int i = 0; i < entries; i++) for (unsigned int i = 0; i < entries; i++)
{ {
data_chain->GetEntry(i); data_chain->GetEntry(i);
if (TITLE == B0_DECAY) 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)))
{
// 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->Fill(B0_M);
if (MASS_HIST_MIN <= B0_M && B0_M <= MASS_HIST_MAX)
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 ((B_M > 4500.) && (B_M < 6000.) && ((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 = " << B_M << " MeV"
<< std::endl;
continue;
} }
else else
{ {
if ((TMath::Abs(Jpsi_M - 3096.9) < 100.) && (B_M > 4500.) && (B_M < 6000.) && ((Hlt1TrackMVADecision) || (Hlt1TwoTrackMVADecision)))
run_event_unique_num = 1;
}
if (TMath::Abs(Jpsi_M - JPSI_MASS) < 100)
{ {
h1_B_M->Fill(B_M);
if (MASS_HIST_MIN <= B_M && B_M <= MASS_HIST_MAX)
h1_B_M_jpsi->Fill(B_M);
}
else if (TMath::Abs(Jpsi_M - PSI2S_MASS) < 100)
{ {
fitting_entries++;
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->Draw();
c1->SaveAs(TString::Format("images/root_hist_%s_bmass.pdf", TITLE).Data()); 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; return 0;
} }
@ -174,7 +233,7 @@ void CreateRooFitAndDraw(TH1D *hist, int fitting_entries)
// TString pdf_name = TString::Format("%s_sigplusbkg", var_id.c_str()); // 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)); 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::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")); // 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("sig_cb"), "Signal", "LP");
// leg1->AddEntry(roo_frame_mass->findObject("bkg_exp"), "Background", "LP"); // leg1->AddEntry(roo_frame_mass->findObject("bkg_exp"), "Background", "LP");
// leg1->AddEntry((TObject *)0, "", ""); // 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(); leg1->Draw();

41
do_all.cpp → old_analysis_macros/do_all.cpp
View File

@ -275,9 +275,13 @@ public:
bools[name] = Variable<Bool_t>(name, cuts); bools[name] = Variable<Bool_t>(name, cuts);
} }
void AddTriggerCut(const std::vector<std::string> &triggers)
void AddTriggerCut(const std::vector<std::string> &triggers_list)
{ {
trigger_cuts.push_back(triggers);
for (const auto &trig : triggers_list)
{
bools[trig.c_str()] = Variable<Bool_t>(trig.c_str(), {});
}
trigger_cuts.push_back(triggers_list);
} }
void AddFourMomentumFor(const char *name) void AddFourMomentumFor(const char *name)
@ -418,7 +422,7 @@ public:
} }
}; };
int do_all()
void BuToHpMuMu()
{ {
TChain *data_chain = new TChain("BuToHpMuMu/DecayTree"); 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"); data_chain->Add("/auto/data/pfeiffer/inclusive_detached_dilepton/data_samples/spruce_magdown_2023_v0_tuple_90000000_v0r0p6288631.root");
VariableCollection col{}; VariableCollection col{};
col.AddFloat("B_BPVFDCHI2", {new GT<float>(36.)});
col.AddDouble("B_CHI2VXNDOF", {new LT(16.)});
col.AddFloat("B_BPVIPCHI2", {new LT<float>(25.)});
col.AddDouble("Jpsi_MAXDOCACHI2", {new LT(36.)});
col.AddFloat("L1_BPVIPCHI2", {new GT<float>(9.)});
col.AddFloat("L2_BPVIPCHI2", {new GT<float>(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<float>(350.)});
col.AddFloat("L2_PT", {new GT<float>(350.)});
col.AddDouble("Jpsi_CHI2DOF", {new LT(9.)});
col.AddFloat("Jpsi_BPVFDCHI2", {new GT<float>(16.)});
col.AddDouble("Jpsi_M", {new LT(5500.), new BT(3096.9 - 100., 3096.9 + 100.)});
col.AddFloat("Hp_PT", {new GT<float>(400.)});
col.AddFloat("Hp_BPVIPCHI2", {new GT<float>(6.)});
col.AddFloat("Hp_P", {new GT<float>(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("L1");
col.AddFourMomentumFor("L2"); col.AddFourMomentumFor("L2");
col.AddFourMomentumFor("Hp"); col.AddFourMomentumFor("Hp");
col.AddTriggerCut({"Hlt2_InclDetDiMuonDecision", "Hlt2_InclDetDiMuon_3BodyDecision", "Hlt2_InclDetDiMuon_4BodyDecision"}); col.AddTriggerCut({"Hlt2_InclDetDiMuonDecision", "Hlt2_InclDetDiMuon_3BodyDecision", "Hlt2_InclDetDiMuon_4BodyDecision"});
col.AddTriggerCut({"Hlt1TrackMVADecision", "Hlt1TwoTrackMVADecision"}); col.AddTriggerCut({"Hlt1TrackMVADecision", "Hlt1TwoTrackMVADecision"});
col.AddTriggerCut({"Hlt2RD_BuToKpMuMuDecision"});
col.Connect(data_chain); col.Connect(data_chain);
@ -494,3 +472,8 @@ int do_all()
return 0; return 0;
} }
int do_all() {
BuToHpMuMu();
return 0;
}

0
simulation.cpp → old_analysis_macros/simulation.cpp
View File

0
simulation_fullstream_B02KpPimMuMu.cpp → old_analysis_macros/simulation_fullstream_B02KpPimMuMu.cpp
View File

0
simulation_fullstream_Bu2KpMuMu.cpp → old_analysis_macros/simulation_fullstream_Bu2KpMuMu.cpp
View File

Write Preview
Loading…
Cancel
Save