#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; }