From 0990514af42bdfae0bdcaf00da4a9f46115a1149 Mon Sep 17 00:00:00 2001 From: Marius Pfeiffer Date: Wed, 8 Nov 2023 16:58:12 +0100 Subject: [PATCH] initial commit --- .gitignore | 1 + general_cuts.cpp | 170 +++++++++++++++++++++++++++++++ pre_selection_cuts.cpp | 225 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 396 insertions(+) create mode 100644 .gitignore create mode 100644 general_cuts.cpp create mode 100644 pre_selection_cuts.cpp diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..a608475 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +*.root \ No newline at end of file diff --git a/general_cuts.cpp b/general_cuts.cpp new file mode 100644 index 0000000..27bfaca --- /dev/null +++ b/general_cuts.cpp @@ -0,0 +1,170 @@ +#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 "TLorentzVector.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 +#include + +const int nBins = 200; + +template +class DrawVar +{ + private: + std::string _identifier; + TVar _min; + TVar _max; + TVar _value = TVar{}; + TH1 *_histogram = nullptr; + + public: + DrawVar(std::string identifier, TVar min, TVar max) + : _identifier { identifier }, _min { min }, _max { max } + { + } + + std::string getIdentifier() { + return _identifier; + } + + TVar getMin() { + return _min; + } + + TVar getMax() { + return _max; + } + + TVar* getValuePointer() { + return &_value; + } +}; + +int general_cuts() +{ + // just some plotting options + gROOT->SetStyle("Plain"); + TPad foo; + const int NRGBs = 5; + const int NCont = 250; + double stops[NRGBs] = {0.00, 0.34, 0.61, 0.84, 1.00}; + double red[NRGBs] = {0.00, 0.00, 0.87, 1.00, 0.51}; + double green[NRGBs] = {0.00, 0.81, 1.00, 0.20, 0.00}; + double blue[NRGBs] = {0.51, 1.00, 0.12, 0.00, 0.00}; + TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont); + gStyle->SetNumberContours(NCont); + gStyle->SetLabelFont(132, "xyz"); + gStyle->SetTitleFont(132, "xyz"); + gStyle->SetLegendFont(132); + gStyle->SetStatFont(132); + gStyle->SetEndErrorSize(10.0); + gStyle->SetOptStat(0); + gStyle->SetOptFit(0); + + std::vector> float_vars{ + DrawVar{"B0_PT", 0., 20000.} + }; + std::vector> double_vars{ + + }; + + // files to load + std::vector filenames = + { + "./B0ToKpPimMuMu_Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_Sprucing23r1_90000000_RD.root"}; + + TChain *chain = new TChain("SpruceRD_B0ToKpPimMuMu/DecayTree"); + for (unsigned int i = 0; i < filenames.size(); i++) + { + chain->Add(filenames.at(i).c_str()); + } + + Double_t B0_M; + chain->SetBranchAddress("B0_M", &B0_M); + TH1D *h1_B0_M = new TH1D("h1_B0_M", "B+ Mass", nBins, 4500, 7500); + + for (size_t i = 0; i < float_vars.size(); i++) + { + auto the_var = &float_vars[i]; + chain->SetBranchAddress(the_var->identifier.c_str(), &(the_var->value)); + the_var->histogram = new TH1D(TString::Format("h1_%s", the_var->identifier.c_str()).Data(), the_var->identifier.c_str(), nBins, the_var->min, the_var->max); + } + + for (size_t i = 0; i < double_vars.size(); i++) + { + auto the_var = &double_vars[i]; + chain->SetBranchAddress(the_var->identifier.c_str(), &(the_var->value)); + the_var->histogram = new TH1D(TString::Format("h1_%s", the_var->identifier.c_str()).Data(), the_var->identifier.c_str(), nBins, the_var->min, the_var->max); + } + + unsigned int entries = chain->GetEntries(); + // loop over all entries in the tree + for (unsigned int i = 0; i < entries; i++) + { + chain->GetEntry(i); + h1_B0_M->Fill(B0_M); + + for (size_t i = 0; i < float_vars.size(); i++) + { + float_vars[i].histogram->Fill(float_vars[i].value); + } + + for (size_t i = 0; i < double_vars.size(); i++) + { + double_vars[i].histogram->Fill(double_vars[i].value); + } + } + + TCanvas *c1 = new TCanvas("c1", "Canvas 1", 0, 0, 800, 600); + + h1_B0_M->Draw(); + + c1->Draw(); + + TCanvas *c2 = new TCanvas("c2", "Canvas 2", 0, 0, 1200, 600); + c2->Divide(4, 2); + + for (size_t i = 0; i < float_vars.size(); i++) + { + c2->cd(i+1); + float_vars[i].histogram->Draw(); + } + + c2->Draw(); + + TCanvas *c3 = new TCanvas("c3", "Canvas 3", 0, 0, 1200, 600); + c3->Divide(4, 2); + + for (size_t i = 0; i < double_vars.size(); i++) + { + c3->cd(i+1); + double_vars[i].histogram->Draw(); + } + + c3->Draw(); + + + return 0; +} \ No newline at end of file diff --git a/pre_selection_cuts.cpp b/pre_selection_cuts.cpp new file mode 100644 index 0000000..9e25d1f --- /dev/null +++ b/pre_selection_cuts.cpp @@ -0,0 +1,225 @@ +#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 "TLorentzVector.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 +#include + +const int nBins = 150; + +RooPlot* CreateRooFit(TH1D* hist); + +int pre_selection_cuts() +{ + // just some plotting options + gROOT->SetStyle("Plain"); + TPad foo; + const int NRGBs = 5; + const int NCont = 250; + double stops[NRGBs] = {0.00, 0.34, 0.61, 0.84, 1.00}; + double red[NRGBs] = {0.00, 0.00, 0.87, 1.00, 0.51}; + double green[NRGBs] = {0.00, 0.81, 1.00, 0.20, 0.00}; + double blue[NRGBs] = {0.51, 1.00, 0.12, 0.00, 0.00}; + TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont); + gStyle->SetNumberContours(NCont); + gStyle->SetLabelFont(132, "xyz"); + gStyle->SetTitleFont(132, "xyz"); + gStyle->SetLegendFont(132); + gStyle->SetStatFont(132); + gStyle->SetEndErrorSize(10.0); + gStyle->SetOptStat(0); + gStyle->SetOptFit(0); + + // files to load + std::vector filenames = + { + "./BuToKpMuMu_Collision23_Beam6800GeV-VeloClosed-MagDown-Excl-UT_RealData_Sprucing23r1_90000000_RD.root"}; + + TChain *chain = new TChain("SpruceRD_BuToKpMuMu/DecayTree"); + for (unsigned int i = 0; i < filenames.size(); i++) + { + chain->Add(filenames.at(i).c_str()); + } + + Double_t Bplus_M, // 4400 -> 8200 + JPsi_M, // 200 -> 6600 + Kplus_M; // 493.6769 -> 493.6779 + + Float_t Bplus_PT, // 0 -> 30 * 10^3 + Bplus_FDCHI2_OWNPV, // 0 -> 7000 + JPsi_IP_OWNPV, // 0 -> 3 + JPsi_PT, // 0 -> 26 * 10^3 + Kplus_IP_OWNPV, // 0 -> 7 + Kplus_PT; // 0 -> 11 * 10^3 + + chain->SetBranchAddress("Bplus_M", &Bplus_M); + chain->SetBranchAddress("Bplus_PT", &Bplus_PT); + chain->SetBranchAddress("Bplus_FDCHI2_OWNPV", &Bplus_FDCHI2_OWNPV); + chain->SetBranchAddress("JPsi_M", &JPsi_M); + chain->SetBranchAddress("JPsi_IP_OWNPV", &JPsi_IP_OWNPV); + chain->SetBranchAddress("JPsi_PT", &JPsi_PT); + chain->SetBranchAddress("Kplus_M", &Kplus_M); + chain->SetBranchAddress("Kplus_IP_OWNPV", &Kplus_IP_OWNPV); + chain->SetBranchAddress("Kplus_PT", &Kplus_PT); + + TH1D *h1_Bplus_M = new TH1D("h1_Bplus_M", "B+ Mass", nBins, 4500, 7500); + TH1D *h1_Bplus_PT = new TH1D("h1_Bplus_PT", "B+ P_{T}", nBins, 0, 15000); + TH1D *h1_Bplus_FDCHI2_OWNPV = new TH1D("h1_Bplus_FDCHI2_OWNPV", "B+ FD #chi^{2}", nBins, 0, 3000); + TH1D *h1_JPsi_M = new TH1D("h1_JPsi_M", "J/#psi Mass", nBins, 200, 6600); + TH1D *h1_JPsi_IP_OWNPV = new TH1D("h1_JPsi_IP_OWNPV", "J/#psi IP", nBins, 0, 1); + TH1D *h1_JPsi_PT = new TH1D("h1_JPsi_PT", "J/#psi P_{T}", nBins, 0, 15000); + TH1D *h1_Kplus_M = new TH1D("h1_Kplus_M", "K^{+} Mass", nBins, 493.676999999, 493.677000001); + TH1D *h1_Kplus_IP_OWNPV = new TH1D("h1_Kplus_IP_OWNPV", "K^{+} IP", nBins, 0, 3); + TH1D *h1_Kplus_PT = new TH1D("h1_Kplus_PT", "K^{+} P_{T}", nBins, 0, 6500); + + TH1D *h1_bkg_Bplus_M = new TH1D("h1_bkg_Bplus_M", "B+ Mass", nBins, 4500, 7500); + TH1D *h1_bkg_Bplus_PT = new TH1D("h1_bkg_Bplus_PT", "B+ P_{T}", nBins, 0, 15000); + TH1D *h1_bkg_Bplus_FDCHI2_OWNPV = new TH1D("h1_bkg_Bplus_FDCHI2_OWNPV", "B+ FD #chi^{2}", nBins, 0, 3000); + TH1D *h1_bkg_JPsi_M = new TH1D("h1_bkg_JPsi_M", "J/#psi Mass", nBins, 200, 6600); + TH1D *h1_bkg_JPsi_IP_OWNPV = new TH1D("h1_bkg_JPsi_IP_OWNPV", "J/#psi IP", nBins, 0, 1); + TH1D *h1_bkg_JPsi_PT = new TH1D("h1_bkg_JPsi_PT", "J/#psi P_{T}", nBins, 0, 15000); + TH1D *h1_bkg_Kplus_M = new TH1D("h1_bkg_Kplus_M", "K^{+} Mass", nBins, 493.676999999, 493.677000001); + TH1D *h1_bkg_Kplus_IP_OWNPV = new TH1D("h1_bkg_Kplus_IP_OWNPV", "K^{+} IP", nBins, 0, 3); + TH1D *h1_bkg_Kplus_PT = new TH1D("h1_bkg_Kplus_PT", "K^{+} P_{T}", nBins, 0, 6500); + + unsigned int entries = chain->GetEntries(); + // loop over all entries in the tree + for (unsigned int i = 0; i < entries; i++) + { + chain->GetEntry(i); + const double JSPI_MASS = 3096.; + + /* + B+ PT > 1500. + B+ FDCHI2 > 100. + J/Psi IP > 0.05 + J/Psi PT > 1000. + K+ IP > 0.2 + K+ PT > 650. + */ + if (Bplus_PT > 1200. && Bplus_FDCHI2_OWNPV > 90. && JPsi_IP_OWNPV > 0.05 && + JPsi_PT > 800. && Kplus_IP_OWNPV > 0.17 && Kplus_PT > 500. && (JSPI_MASS - 150. < JPsi_M && JPsi_M < JSPI_MASS + 150.)) { + h1_Bplus_M->Fill(Bplus_M); + } else { + h1_bkg_Bplus_M->Fill(Bplus_M); + } + + if (Bplus_M > 5700.) + { + h1_bkg_Bplus_PT->Fill(Bplus_PT); + h1_bkg_Bplus_FDCHI2_OWNPV->Fill(Bplus_FDCHI2_OWNPV); + h1_bkg_JPsi_M->Fill(JPsi_M); + h1_bkg_JPsi_IP_OWNPV->Fill(JPsi_IP_OWNPV); + h1_bkg_JPsi_PT->Fill(JPsi_PT); + h1_bkg_Kplus_M->Fill(Kplus_M); + h1_bkg_Kplus_IP_OWNPV->Fill(Kplus_IP_OWNPV); + h1_bkg_Kplus_PT->Fill(Kplus_PT); + } + else + { + h1_Bplus_PT->Fill(Bplus_PT); + h1_Bplus_FDCHI2_OWNPV->Fill(Bplus_FDCHI2_OWNPV); + h1_JPsi_M->Fill(JPsi_M); + h1_JPsi_IP_OWNPV->Fill(JPsi_IP_OWNPV); + h1_JPsi_PT->Fill(JPsi_PT); + h1_Kplus_M->Fill(Kplus_M); + h1_Kplus_IP_OWNPV->Fill(Kplus_IP_OWNPV); + h1_Kplus_PT->Fill(Kplus_PT); + } + } + + TCanvas *c1 = new TCanvas("c1", "Canvas 1", 0, 0, 800, 600); + + + //h1_bkg_Bplus_M->SetLineColor(kRed); + //h1_bkg_Bplus_M->Draw(); + h1_Bplus_M->Draw(); + + c1->Draw(); + + std::vector small_histos{h1_Bplus_PT, h1_Bplus_FDCHI2_OWNPV, h1_JPsi_M, h1_JPsi_IP_OWNPV, h1_JPsi_PT, h1_Kplus_M, h1_Kplus_IP_OWNPV, h1_Kplus_PT}; + std::vector small_bkg_histos{h1_bkg_Bplus_PT, h1_bkg_Bplus_FDCHI2_OWNPV, h1_bkg_JPsi_M, h1_bkg_JPsi_IP_OWNPV, h1_bkg_JPsi_PT, h1_bkg_Kplus_M, h1_bkg_Kplus_IP_OWNPV, h1_bkg_Kplus_PT}; + + TCanvas *c2 = new TCanvas("c2", "Canvas 2", 0, 0, 1200, 600); + c2->Divide(4, 2); + + for (size_t i = 0; i < small_histos.size(); i++) + { + c2->cd(i + 1); + + small_histos[i]->Draw(); + small_bkg_histos[i]->SetLineColor(kRed); + small_bkg_histos[i]->Draw("SAME"); + } + + c2->Draw(); + + TCanvas *c3 = new TCanvas("c3", "Canvas 3", 0, 0, 800, 600); + + auto fitFrame = CreateRooFit(h1_Bplus_M); + fitFrame->Draw(); + + c3->Draw(); + + return 0; +} + +RooPlot* CreateRooFit(TH1D* hist) { + RooRealVar roo_var_mass("roo_var_mass", "B+ Mass Variable", 4500., 7500.); + roo_var_mass.setRange("fitting_range", 4700., 6500.); + + TString hist_name = "roohist_bplus_M"; + RooDataHist roohist_bplus_M(hist_name, "B Plus Mass Histogram", roo_var_mass, RooFit::Import(*hist)); + + RooRealVar roo_sig_gauss_mean("roo_sig_gauss_mean", "Mass Gauss Mean", 5250., 5100., 5400.); + RooRealVar roo_sig_gauss_sigma("roo_sig_gauss_sigma", "Mass Gauss Sigma", 60., 0., 150.); + + RooGaussian roo_sig_gauss("roo_sig_gauss", "B+ Mass Signal Gaussian", roo_var_mass, roo_sig_gauss_mean, roo_sig_gauss_sigma); + + RooRealVar roo_bkg_exp_c("roo_bkg_exp_c", "Background C", -0.00117202, -0.003, -0.001); + RooExponential roo_bkg_exp("roo_bkg_exp", "B+ Mass Background Exp", roo_var_mass, roo_bkg_exp_c); + + RooRealVar roo_var_mass_nsig("roo_var_mass_nsig", "B+ Mass N Signal", 0., hist->GetEntries()); + RooRealVar roo_var_mass_nbkg("roo_var_mass_nbkg", "B+ Mass N Background", 0., hist->GetEntries()); + + TString pdf_name = "roo_pdf_sig_plus_bkg"; + RooAddPdf roo_pdf_sig_plus_bkg(pdf_name, "Sig + Bkg PDF", + RooArgList(roo_sig_gauss, roo_bkg_exp), + RooArgList(roo_var_mass_nsig, roo_var_mass_nbkg)); + + RooPlot *roo_frame_mass = roo_var_mass.frame(); + roohist_bplus_M.plotOn(roo_frame_mass, RooFit::Binning(nBins), RooFit::Name(hist_name)); + + RooFitResult *fitres = roo_pdf_sig_plus_bkg.fitTo(roohist_bplus_M, RooFit::Save(), RooFit::PrintLevel(1), RooFit::Range("fitting_range")); + + roo_pdf_sig_plus_bkg.plotOn(roo_frame_mass, RooFit::VisualizeError(*fitres, 1), RooFit::FillColor(kOrange + 1), RooFit::FillStyle(3144)); + roo_pdf_sig_plus_bkg.plotOn(roo_frame_mass, RooFit::LineColor(kRed), RooFit::LineStyle(kSolid), RooFit::Range("fitting_range"), RooFit::Name(pdf_name)); + roo_pdf_sig_plus_bkg.plotOn(roo_frame_mass, RooFit::Components(RooArgSet(roo_bkg_exp)), RooFit::LineColor(kBlue), RooFit::LineStyle(kDashed), RooFit::Range("fitting_range")); + + roo_sig_gauss.paramOn(roo_frame_mass, RooFit::Layout(0.60, 0.99, 0.90)); + roo_frame_mass->getAttText()->SetTextSize(0.027); + + return roo_frame_mass; +} \ No newline at end of file