diff --git a/BDT1_applyer.py b/BDT1_applyer.py new file mode 100644 index 0000000..ff34b12 --- /dev/null +++ b/BDT1_applyer.py @@ -0,0 +1,28 @@ +import os,sys +#from ROOT import TH1D,TH2D,TFile,TTree,TCanvas +import sys +import array + +# Import the other script + +# Script which creates the BDT_rootfile +import Branches_name +directory_in_str = './Cutted_files/' +#directory_in_str = './MC_Cutted_files/' +directory = os.fsencode(directory_in_str) + +path_to_BDT = '/tmva-branch-adder/src/' + + +#iterating over the files +for file in os.listdir(directory): + + filename = os.fsdecode(file) + +# if filename.endswith(".root"):# and filename.startswith("Bu2KpipiMM"): + if filename.endswith(".root") and filename.startswith("Bu2KpipiMM"): + print(filename) + + #Branches_name.main(filename,"mc_samples") + Branches_name.main(filename,"data_samples") + # break \ No newline at end of file diff --git a/BDT1_branches.py b/BDT1_branches.py new file mode 100644 index 0000000..a90ca72 --- /dev/null +++ b/BDT1_branches.py @@ -0,0 +1,33 @@ +from curses import echo +import os, sys + +def main(directory_name): + + path_to_directory = '/panfs/mleydier/Python_training/Project/' + directory = path_to_directory + directory_name # + '/' + + #iterating over the files + for file in os.listdir(directory): + filename = os.fsdecode(file) + + if filename.endswith(".root"):# and filename.startswith("Bu2KpipiMM"): +## if filename.endswith(".root") and filename.startswith("Cut_Bu2KpipiMM"): + print(filename) + #path_datafile = '/panfs/mleydier/Python_training/Project/Cutted_files/' + path_datafile ='/panfs/mleydier/Python_training/Project/MC_Cutted_files/' + #path_outputfile ='/panfs/mleydier/Python_training/Project/BDT_Cutted_files/' + path_outputfile ='/panfs/mleydier/Python_training/Project/BDT_MC_Cutted_files/' + + indatafile = path_datafile + filename + outdatafile = path_outputfile + "BDT1_" +filename + + # os.system allows to make the same commands in a script and in the shell + # spaces in the string are very important! + print("./main " + indatafile + " DecayTree factory_BDTG.weights.xml " + outdatafile + " BDT1") + print("******* Before os.system *******") + os.system("./main " + indatafile + " DecayTree factory_BDTG.weights.xml " + outdatafile + " BDT1") + print("******* After os.system *******") + +#end function +if(__name__=='__main__'): + sys.exit(main(sys.argv[1])) \ No newline at end of file diff --git a/Branches_name.py b/Branches_name.py new file mode 100644 index 0000000..5281ec9 --- /dev/null +++ b/Branches_name.py @@ -0,0 +1,108 @@ +import os,sys +from ROOT import TH1D,TH2D,TFile,TTree,TCanvas +import sys +import array +import numpy as np + +## MAIN FUNCTION ## +def main(FileName,datafile): + + inFileName = os.path.join(datafile,FileName) # Not relevant for the branches selection + inFile=TFile(inFileName) + intree = inFile.Get("DecayTree") + +# New cutted file +# cutted_sample_path = '/panfs/mleydier/Python_training/Project/MC_Cutted_files/' + cutted_sample_path = '/panfs/mleydier/Python_training/Project/Cutted_files/' + cutted_sample_name = 'Cut_' + FileName + + newfile = TFile(cutted_sample_path + cutted_sample_name,"RECREATE") + +#Branches to keep + + intree.SetBranchStatus("*",0) + branches_to_keep = ["K_PIDK","JPsi_MM","muplus_TRACK_GhostProb","muplus_ProbNNghost","muminus_ProbNNghost", + "muminus_TRACK_GhostProb","K_ProbNNghost","K_TRACK_GhostProb","B_L0MuonDecision_TOS","B_L0DiMuonDecision_TOS", + "B_DIRA_OWNPV","B_IPCHI2_OWNPV","B_DTF_PVandJpsi_CHI2NDOF","B_ENDVERTEX_CHI2","B_PT", "K1_1270_ENDVERTEX_CHI2", + "JPsi_FDCHI2_OWNPV","K_PT","muplus_PT","muminus_PT","JPsi_P","JPsi_PT","JPsi_CosTheta","B_MM", "piminus_PT", "piplus_PT"]#,"B_LOKI_DTF_CHI2NDOF"] + + for keepbranch in branches_to_keep: + intree.SetBranchStatus(keepbranch, 1) # activate branches of interest + +# Apply Sonia's cut : + LooseCut = "K_PIDK>-2 && muplus_TRACK_GhostProb<0.3 && muminus_TRACK_GhostProb<0.3 && K_TRACK_GhostProb<0.3 && (B_L0MuonDecision_TOS || B_L0DiMuonDe\ +cision_TOS) && B_DIRA_OWNPV > 0.9999 && B_DTF_PVandJpsi_CHI2NDOF<5" + + + newtree = intree.CopyTree(LooseCut) + + # pseudorapity eta = -ln() with JPsi_costheta + #Cos_Theta = intree.AsMatrix(["JPsi_CosTheta"]) + eta_arrqy = array.array('d',[0]) + phi_pt_qrray = array.array('d',[0]) + b_dtf_array = array.array('d',[0]) + kpipi_vtx_array = array.array('d',[0]) + jpsi_fd_array = array.array('d',[0]) + mu_plus_array = array.array('d',[0]) + mu_minus_array = array.array('d',[0]) + kaon_array = array.array('d',[0]) + + + #define all qrrqys for the brqnches to rename + #define all newtree.Branch for those branches + + #JPsi_eta = -1*np.log(np.tan(0.5*Cos_Theta)) + eta_branch = newtree.Branch("Jpsi_ETA", eta_arrqy , "Jpsi_ETA/D") + phi_pt_branch = newtree.Branch("phi_PT", phi_pt_qrray , "phi_PT/D") + b_dtf_array_branch = newtree.Branch("B_LOKI_DTF_CHI2NDOF", b_dtf_array, "B_LOKI_DTF_CHI2NDOF/D") + B_VxtChi2_branch = newtree.Branch("B_VtxChi2_234", kpipi_vtx_array, "B_VtxChi2_234/D") + jpsi_fd_a_branch = newtree.Branch("Jpsi_FDCHI2_OWNPV", jpsi_fd_array, "Jpsi_FDCHI2_OWNPV/D") + muplus_branch = newtree.Branch("mu1_PT", mu_plus_array, "mu1_PT/D") + mu_minus_branch = newtree.Branch("mu2_PT", mu_minus_array, "mu2_PT/D") + kaon_array_branch = newtree.Branch("Kaon_PT", kaon_array, "Kaon_PT/D") + + print(" Number of entries befor loop : ", newtree.GetEntries() ) + + for i in range(newtree.GetEntries()): + newtree.GetEntry(i) + eta_arrqy[0] = -1*np.log(np.tan(0.5*newtree.JPsi_CosTheta)) + phi_pt_qrray[0] = newtree.piplus_PT + newtree.piminus_PT + #print(eta_arrqy[0]) + #copy the value from old brqnch to array of new branch (element 0) + b_dtf_array[0] = newtree.B_DTF_PVandJpsi_CHI2NDOF + kpipi_vtx_array[0] = newtree.K1_1270_ENDVERTEX_CHI2 + jpsi_fd_array[0] = newtree.JPsi_FDCHI2_OWNPV + mu_plus_array[0] = newtree.muplus_PT + mu_minus_array[0] = newtree.muminus_PT + kaon_array[0] = newtree.K_PT + + + + eta_branch.Fill() + phi_pt_branch.Fill() + b_dtf_array_branch.Fill() + B_VxtChi2_branch.Fill() + jpsi_fd_a_branch.Fill() + muplus_branch.Fill() + mu_minus_branch.Fill() + kaon_array_branch.Fill() + + + print(" Number After loop : ", newtree.GetEntries() ) + + list_branches = [key.GetName() for key in intree.GetListOfBranches()] + for branch in list_branches: + if branch in branches_to_keep: + print(branch) + + +# Writing in the new tree + print("Tree written with ", newtree.GetEntries(), " entries copied") + print("File saved at location: ", cutted_sample_path+cutted_sample_name) + newtree.Write() + + inFile.Close() + print("File closed") +#end function +if(__name__=='__main__'): + sys.exit(main(sys.argv[1])) diff --git a/Final_plots.py b/Final_plots.py new file mode 100644 index 0000000..86cf666 --- /dev/null +++ b/Final_plots.py @@ -0,0 +1,249 @@ +from cmath import sqrt +import os,sys +import ROOT +from ROOT import TH1D,TH2D,TFile,TTree,TCanvas, TChain +import sys +import array +import numpy as np +import matplotlib.pyplot as plt +import scipy + +#ROOT.gROOT.SetBatch(1) + +## MAIN FUNCTION ## +def main(bins): + + #canvas1 = TCanvas("canvas1") + + bins = sys.argv[1] + #plotFileName2 = sys.argv[1] #"plotfit.pdf" + + BDT_files = "BDT_MC_Cutted_files" + Data_Cutted_files = "BDT_Cutted_files" # Muons cutted files and BDT added measured + Data_files = "data_samples" + + ########### Efficiency = measured / simulated ############## + + BDT_files_MC = TChain("DecayTree") + + BDT_files_MC.Add(os.path.join(BDT_files,"BDT1_Cut_Bu2KpipiJPsi-MM-MC-12145090-2018-MagDown-Sim09h-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) + BDT_files_MC.Add(os.path.join(BDT_files,"BDT1_Cut_Bu2KpipiJPsi-MM-MC-12145090-2018-MagUp-Sim09h-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) + + + Cutted_files_Data = TChain("DecayTree") + + Cutted_files_Data.Add(os.path.join(Data_Cutted_files,"BDT1_Cut_Bu2KpipiMM-Data-2018-MagDown-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) + Cutted_files_Data.Add(os.path.join(Data_Cutted_files,"BDT1_Cut_Bu2KpipiMM-Data-2018-MagUp-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) + + Raw_files_Data = TChain("DecayTree") + + Raw_files_Data.Add(os.path.join(Data_files,"Bu2KpipiMM-Data-2018-MagDown-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) + Raw_files_Data.Add(os.path.join(Data_files,"Bu2KpipiMM-Data-2018-MagUp-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) + + +################ ARRAY DECLARATION ################ + N_bins = int(bins) + BDT = np.empty(N_bins, dtype=object) + + +# *********** background *********** # + N_BG = np.empty(N_bins, dtype=object) + N_BG_BDT = np.empty(N_bins, dtype=object) + efficiency_B = np.empty(N_bins, dtype=object) + + significance = np.empty(N_bins, dtype=object) +# *********** signal *********** # + + efficiency_S_MC = np.empty(N_bins, dtype=object) + + N_S_BDT_MC = np.empty(N_bins, dtype=object) + N_S_MC = np.empty(N_bins, dtype=object) + Sqrt_B = np.empty(N_bins, dtype=object) +# *********** Errors *********** # + yerr = np.empty(N_bins, dtype=object) + N_BG_err = np.empty(N_bins, dtype=object) + N_S_err = np.empty(N_bins, dtype=object) + Efficiency_B_err = np.empty(N_bins, dtype=object) + Efficiency_S_err = np.empty(N_bins, dtype=object) + mass = np.empty(N_bins, dtype=object) + + + + +#=============== FOR LOOP ===============# + + for i in range(0,N_bins): + + BDT[i] = -1 + i*2/N_bins + + + #---------background cut without the peak---------# + + # BDT Cutted data | Background window + N_BG_BDT_i = Cutted_files_Data.Draw("B_MM>>histo_BG_BDT"," BDT1 > " + str(BDT[i]) + "& B_MM<5200 || B_MM>5350") + # Raw data | Background window + N_BG_i = Cutted_files_Data.Draw("B_MM>>histo_BG_BDT","B_MM<5200 || B_MM>5350") + + N_BG_BDT[i] = N_BG_BDT_i + N_BG[i] = N_BG_i + + efficiency_B[i] = N_BG_BDT[i] /N_BG[i] + + N_BG_err[i] = np.sqrt(N_BG_BDT_i) + Efficiency_B_err[i] = np.sqrt(N_BG_BDT_i)/N_BG_i +(N_BG_BDT_i*np.sqrt(N_BG_i))/(N_BG_i)**2 + + #------------------signal cut only the peak------------------# + + # BDT Cutted MC data | Signal window + N_S_BDT_MC_i = BDT_files_MC.Draw("B_MM>>histo_PEAK_BDT"," BDT1 > " + str(BDT[i]) + "& B_MM>5200 && B_MM<5350") + # BDT MC data | Signal window + N_S_MC_i = BDT_files_MC.Draw("B_MM>>histo_PEAK_BDT", "5200>histo_BG_BDT"," B_MM>5200 && B_MM<5350 ") + #S_new = N_BG_BDT_PEAK - B_new + # Significance ONLY OVER THE PEAK REGION + + S_BDT_PEAK = np.empty(N_bins, dtype=object) + B_BDT_PEAK = np.empty(N_bins, dtype=object) + + for i in range(0,N_bins): + BDT[i] = -1 + i*2/N_bins + #BDT_threshold=0.1 + Cutted_files_Data.Draw("B_MM>>histo_to_fit"," BDT1 > " + str(BDT[i]) + "& B_MM > 5000 && B_MM<5200 || B_MM>5350 && B_MM<6900") #B_MM > 3000 && B_MM<5200 || + + Histo_to_fit = ROOT.gDirectory.Get("histo_to_fit") + #Histo_to_fit.GetXaxis().SetTitle("Mass [MeV]") + #Histo_to_fit.GetYaxis().SetTitle("Number of events") + #Histo_to_fit.SetStats(0) + + myfitfunc = ROOT.TF1("matthieu", "[a]*exp([b]*x )", Histo_to_fit.GetXaxis().GetXmin(), Histo_to_fit.GetXaxis().GetXmax()) #[a]*exp([b]*x) [a]+[b]*x + myfitfunc.SetParameters(1,1) + Histo_to_fit.Fit("matthieu","", "", Histo_to_fit.GetXaxis().GetXmin(), Histo_to_fit.GetXaxis().GetXmax()) #,"B_MM<5200 || B_MM>5350" ) + A_i = myfitfunc.GetParameter(0) + B_i = myfitfunc.GetParameter(1) + + print("A_i = ", A_i ," and B_i= ", B_i) + plt.plot(Histo_to_fit) + #canvas1.Update() + #canvas1.SaveAs(plotFileName2) + + + # Histo_to_fit = ROOT.gDirectory.Get("histo_to_fit") + # Fit = ROOT.TF1("expfit","expo","10","-0.001") + # Histo_to_fit.Fit(Fit ,"E") + + # A_i = #(BDT) data extracted from the fit. + # B_i = #(BDT) + + B_BDT_PEAK[i] = (A_i/B_i)*( np.exp(B_i*M_max)- np.exp(B_i*M_min) ) + S_BDT_PEAK[i] = Cutted_files_Data.Draw("B_MM>>histo_PEAK"," BDT1 > " + str(BDT[i]) + "& B_MM<5200 && B_MM<5350")# - B_BDT_PEAK[i] + significance[i] = S_BDT_PEAK[i]/np.sqrt(B_BDT_PEAK[i]) + + # Plot Significance + + + #plt.errorbar(BDT, B_BDT_PEAK, xerr =None,yerr=None, fmt = 'o',color = 'navy', + # ecolor = 'cornflowerblue', elinewidth = 0, capsize=10,label='B') + #plt.errorbar(BDT, S_BDT_PEAK, xerr =None,yerr=None, fmt = 'o',color = 'red', + # ecolor = 'cornflowerblue', elinewidth = 0, capsize=10,label='S') + # Display graph + #plt.rc('ytick', labelsize=15) + #plt.rc('xtick', labelsize=15) + #plt.xlabel('BDT threshold') + #plt.ylabel('B and S') + #plt.legend(fontsize=15,loc="lower left") + #plt.show() + + + #plt.errorbar(BDT, significance, xerr =None,yerr=None, fmt = 'o',color = 'navy', + # ecolor = 'cornflowerblue', elinewidth = 0, capsize=10,label='Signification') + # Display graph + #plt.rc('ytick', labelsize=15) + #plt.rc('xtick', labelsize=15) + #plt.xlabel('BDT threshold') + #plt.ylabel('Significance') + #plt.legend(fontsize=15,loc="lower left") + #plt.show() + + + + + + + #================================= FITTING =================================# + + #plt.figure() + #plt.plot(BDT,efficiency_B) + #plt.plot(BDT,efficiency_S) + #plt.xlabel("BDT threshold ") + #plt.ylabel("Efficiency background") + #plt.show() + + + # plot of "significance" S/sqrt(B) + #plt.figure() + #plt.plot(BDT, N_cut_S/Sqrt_B) + #plt.xlabel("BDT threshold ") + #plt.ylabel("Significance ") + #plt.show() + +''' +#end function +if(__name__=='__main__'): + sys.exit(main(sys.argv[1])) diff --git a/Plot_B_mass.py b/Plot_B_mass.py new file mode 100644 index 0000000..ff1d841 --- /dev/null +++ b/Plot_B_mass.py @@ -0,0 +1,162 @@ +import ROOT +import os,sys +from ROOT import TH1D,TH2D,TFile,TTree,TCanvas,TChain +import sys +import array + +#histFileName = sys.argv [1] +plotFileName1 = sys.argv [1] +plotFileName2 = sys.argv [2] + +#histFile = TFile.Open(histFileName,"READ") +#dataHisto = histFile.Get("DecayTree") + +#if not dataHisto: +# print("Failed to get data histogram") +# sys.exit(1) + +canvas1 = TCanvas("canvas1") +#canvas1.Print(plotFileName1+"[") + +# Plotting # +# "B_MM" +data_files = "data_samples" +cutted_data_files = "Cutted_files" +# + +#Raw data +Run1 = TChain("DecayTree") + +Run1.Add(os.path.join(data_files,"Bu2KpipiMM-Data-2011-MagDown-StrippingBu2LLK_fiducialMM_preselPIDMM_selected.root")) +Run1.Add(os.path.join(data_files,"Bu2KpipiMM-Data-2011-MagUp-StrippingBu2LLK_fiducialMM_preselPIDMM_selected.root")) +Run1.Add(os.path.join(data_files,"Bu2KpipiMM-Data-2012-MagDown-StrippingBu2LLK_fiducialMM_preselPIDMM_selected.root")) +Run1.Add(os.path.join(data_files,"Bu2KpipiMM-Data-2012-MagUp-StrippingBu2LLK_fiducialMM_preselPIDMM_selected.root")) + +#Cutted data +Run1_Cut = TChain("DecayTree") + +Run1_Cut.Add(os.path.join(cutted_data_files,"Cut_Bu2KpipiMM-Data-2011-MagDown-StrippingBu2LLK_fiducialMM_preselPIDMM_selected.root")) +Run1_Cut.Add(os.path.join(cutted_data_files,"Cut_Bu2KpipiMM-Data-2011-MagUp-StrippingBu2LLK_fiducialMM_preselPIDMM_selected.root")) +Run1_Cut.Add(os.path.join(cutted_data_files,"Cut_Bu2KpipiMM-Data-2012-MagDown-StrippingBu2LLK_fiducialMM_preselPIDMM_selected.root")) +Run1_Cut.Add(os.path.join(cutted_data_files,"Cut_Bu2KpipiMM-Data-2012-MagUp-StrippingBu2LLK_fiducialMM_preselPIDMM_selected.root")) + + + +Run1.SetLineColor(ROOT.kBlack) +Run1_Cut.SetLineColor(ROOT.kRed) + +#Draw +#N_events = Run1.Draw("B_MM>>my_histogram1(50,3700,6900)","") + +Run1.Draw("B_MM>>my_histogram1(10,3700,6900)","") +Run1_Cut.Draw("B_MM>>my_histogram11(10,3700,6900)","B_MM > 5000 && B_MM<5200 || B_MM>5350","same") +#N_cut_events = Run1_Cut.Draw("B_MM>>my_histogram11(50,3700,6900)","B_MM<5000","same") +#print("N_events = ", N_events, "N_cut_events = ",N_cut_events) + +# Retrieve the histograms +my_histogram1 = ROOT.gDirectory.Get("my_histogram1") +my_histogram11 = ROOT.gDirectory.Get("my_histogram11") + +#Styling +my_histogram1.GetXaxis().SetTitle("Mass [MeV]") +my_histogram1.GetYaxis().SetTitle("Number of events") +my_histogram1.SetStats(0) +my_histogram11.SetStats(0) + +leg1 = ROOT.TLegend(.6,.6,.9,.7) +leg1.SetBorderSize(0) +leg1.SetFillColor(0) +leg1.SetFillStyle(0) +leg1.SetTextFont(42) +leg1.SetTextSize(0.035) +leg1.AddEntry(Run1,"Raw data","L") +leg1.AddEntry(Run1_Cut,"Cutted data","L") +leg1.Draw() + +latex = ROOT.TLatex() +latex.SetNDC() +latex.SetTextSize(0.045) +latex.DrawText(0.6,0.8,"LHCb Run 1") +latex.SetTextSize(0.02) +#latex.DrawText(0.7,0.77,"B \rightarrow Kpipi events") + +# Update the styling part +canvas1.Update() + + +canvas1.SaveAs(plotFileName1) + + + +canvas2 = TCanvas("canvas2") + +#Raw data +Run2 = TChain("DecayTree") + +#Run2.Add(os.path.join(data_files,"Bu2KpipiMM-Data-2015-MagDown-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) +#Run2.Add(os.path.join(data_files,"Bu2KpipiMM-Data-2015-MagUp-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) +#Run2.Add(os.path.join(data_files,"Bu2KpipiMM-Data-2016-MagDown-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) +#Run2.Add(os.path.join(data_files,"Bu2KpipiMM-Data-2016-MagUp-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) +#Run2.Add(os.path.join(data_files,"Bu2KpipiMM-Data-2017-MagDown-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) +#Run2.Add(os.path.join(data_files,"Bu2KpipiMM-Data-2017-MagUp-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) +Run2.Add(os.path.join(data_files,"Bu2KpipiMM-Data-2018-MagDown-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) +Run2.Add(os.path.join(data_files,"Bu2KpipiMM-Data-2018-MagUp-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) + +#Cutted data +Run2_Cut = TChain("DecayTree") + +#Run2_Cut.Add(os.path.join(cutted_data_files,"Cut_Bu2KpipiMM-Data-2015-MagDown-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) +#Run2_Cut.Add(os.path.join(cutted_data_files,"Cut_Bu2KpipiMM-Data-2015-MagUp-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) +#Run2_Cut.Add(os.path.join(cutted_data_files,"Cut_Bu2KpipiMM-Data-2016-MagDown-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) +#Run2_Cut.Add(os.path.join(cutted_data_files,"Cut_Bu2KpipiMM-Data-2016-MagUp-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) +#Run2_Cut.Add(os.path.join(cutted_data_files,"Cut_Bu2KpipiMM-Data-2017-MagDown-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) +#Run2_Cut.Add(os.path.join(cutted_data_files,"Cut_Bu2KpipiMM-Data-2017-MagUp-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) +Run2_Cut.Add(os.path.join(cutted_data_files,"Cut_Bu2KpipiMM-Data-2018-MagDown-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) +Run2_Cut.Add(os.path.join(cutted_data_files,"Cut_Bu2KpipiMM-Data-2018-MagUp-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) + +#Styling +Run2.SetLineColor(ROOT.kBlack) +Run2_Cut.SetLineColor(ROOT.kRed) + +# Draw +Run2.Draw("B_MM>>my_histogram2(100,4000,6900)","B_MM <5200 || B_MM>5350 ") +Run2_Cut.Draw("B_MM>>my_histogram22(100,4000,6900)","B_MM > 5000 && B_MM<5200 || B_MM>5350","same") +#Run2_Cut.fit() + + +# Retrieve the histograms +my_histogram2 = ROOT.gDirectory.Get("my_histogram2") +my_histogram22 = ROOT.gDirectory.Get("my_histogram22") + +#Styling +my_histogram2.GetXaxis().SetTitle("Mass [MeV]") +my_histogram2.GetYaxis().SetTitle("Number of events") +my_histogram2.SetStats(0) +my_histogram22.SetStats(0) + +leg2 = ROOT.TLegend(.6,.6,.9,.7) +leg2.SetBorderSize(0) +leg2.SetFillColor(0) +leg2.SetFillStyle(0) +leg2.SetTextFont(42) +leg2.SetTextSize(0.035) +leg2.AddEntry(Run2,"Raw data","L") +leg2.AddEntry(Run2_Cut,"Cutted data","L") +leg2.Draw() + +latex = ROOT.TLatex() +latex.SetNDC() +latex.SetTextSize(0.045) +latex.DrawText(0.6,0.8,"LHCb Run 2") +latex.SetTextSize(0.02) +#latex.DrawText(0.7,0.77,"B \rightarrow Kpipi events") + +# Update the styling part +canvas2.Update() + + +canvas2.SaveAs(plotFileName2) + + +# End of plotting + diff --git a/Plots_BDT.py b/Plots_BDT.py new file mode 100644 index 0000000..cf7790d --- /dev/null +++ b/Plots_BDT.py @@ -0,0 +1,108 @@ +import ROOT +import os,sys +from ROOT import TH1D,TH2D,TFile,TTree,TCanvas,TChain +import sys +import array + +#histFileName = sys.argv [1] +plotFileName1 = sys.argv [1] + + +canvas1 = TCanvas("canvas1") +#canvas1.Print(plotFileName1+"[") + +# Plotting # +# "B_MM" +data_files = "data_samples" +#compaired_data_files = "BDT_Cutted_files" +compaired_data_files = "Cutted_files" +muons_files = "Muons_samples" +# + +#Raw data +Run1 = TChain("DecayTree") + +Run1.Add(os.path.join(data_files,"Bu2KpipiMM-Data-2018-MagDown-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) +Run1.Add(os.path.join(data_files,"Bu2KpipiMM-Data-2018-MagUp-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) + + +#Cutted data +Run1_Cut = TChain("DecayTree") + +Run1_Cut.Add(os.path.join(compaired_data_files,"Cut_Bu2KpipiMM-Data-2018-MagDown-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) +Run1_Cut.Add(os.path.join(compaired_data_files,"Cut_Bu2KpipiMM-Data-2018-MagUp-StrippingBu2LLK_fiducialMM16_preselPIDMM16_selected.root")) + + +# Muons data +Run1_Muons = TChain("DecayTree") + +Run1_Muons.Add(os.path.join(muons_files,"BDT2mu4_added_Muon4Aliased_BDT2mu3_added_MuonAliased_BDT1added_preselectedDimuon_inclusiveJpsi_gangaAll_2018_Data_MD_211104.root")) +Run1_Muons.Add(os.path.join(muons_files,"BDT2mu4_added_Muon4Aliased_BDT2mu3_added_MuonAliased_BDT1added_preselectedDimuon_inclusiveJpsi_gangaAll_2018_Data_MU_211104.root")) + + + +#Styling +Run1.SetLineColor(ROOT.kBlue) +Run1_Cut.SetLineColor(ROOT.kRed) +Run1_Muons.SetLineColor(ROOT.kBlack) + +#Drawing + +#Run1.Draw("K_PT>>raw_histo(1000,0,7000)","") +#Run1_Cut.Draw("Kaon_PT>>cut_histo(1000,0,7000)","","same") +#Run1_Muons.Draw("Kaon_PT>>muon_histo(1000,0,7000)","","same") + +Run1.Draw("piplus_PT+piminus_PT>>raw_histo(1000,0,8000)","") +Run1_Cut.Draw("phi_PT>>cut_histo(1000,0,8000)","","same") +Run1_Muons.Draw("phi_PT>>muon_histo(1000,0,8000)","","same") + +#B_DTF_PVandJpsi_CHI2NDOF +#Run1.Draw("log10(abs(mu1_PT)>>raw_histo","") +#Run1_Cut.Draw("log10(abs(mu1_PT)>>cut_histo","","same") +#Run1_Muons.Draw("log10(abs(mu1_PT)>>muon_histo","","same") + +# Retrieve the histograms +raw_histo = ROOT.gDirectory.Get("raw_histo") +cut_histo = ROOT.gDirectory.Get("cut_histo") +muon_histo = ROOT.gDirectory.Get("muon_histo") + +#Normalization +factor = 1. +raw_histo.Scale(factor/raw_histo.Integral()) +cut_histo.Scale(factor/cut_histo.Integral()) +muon_histo.Scale(factor/muon_histo.Integral()) + +#Styling +raw_histo.GetXaxis().SetTitle("phi_PT") +raw_histo.GetYaxis().SetTitle("Events (Normalised) ") +raw_histo.SetStats(0) +#cut_histo.SetStats(0) + +leg1 = ROOT.TLegend(.6,.6,.9,.7) +leg1.SetBorderSize(0) +leg1.SetFillColor(0) +leg1.SetFillStyle(0) +leg1.SetTextFont(42) +leg1.SetTextSize(0.035) + +leg1.AddEntry(raw_histo,"Raw data","L") +leg1.AddEntry(cut_histo,"Cutted data","L") +leg1.AddEntry(muon_histo,"Muons data","L") + +leg1.Draw() + +latex = ROOT.TLatex() +latex.SetNDC() +latex.SetTextSize(0.045) +latex.DrawText(0.6,0.8,"LHCb 2018") +latex.SetTextSize(0.02) +#latex.DrawText(0.7,0.77,"B \rightarrow Kpipi events") + +# Update the styling part +canvas1.Update() + + +canvas1.SaveAs(plotFileName1) + + + diff --git a/applying_script.py b/applying_script.py new file mode 100644 index 0000000..3c59b7c --- /dev/null +++ b/applying_script.py @@ -0,0 +1,27 @@ +import os,sys +#from ROOT import TH1D,TH2D,TFile,TTree,TCanvas +import sys +import array +# Import the other script +import Branches_name +directory_in_str = './data_samples/' +#directory_in_str = './mc_samples/' +directory = os.fsencode(directory_in_str) + + +#iterating over the files +for file in os.listdir(directory): + + filename = os.fsdecode(file) + + # if filename.endswith(".root"):# and filename.startswith("Bu2KpipiMM"): + if filename.endswith(".root") and filename.startswith("Bu2KpipiMM"): + print(filename) + + #Branches_name.main(filename,"mc_samples") + Branches_name.main(filename,"data_samples") + #break + + + + diff --git a/efficiency.py b/efficiency.py new file mode 100644 index 0000000..0647d04 --- /dev/null +++ b/efficiency.py @@ -0,0 +1,81 @@ +from cmath import sqrt +import os,sys +import ROOT +from ROOT import TH1D,TH2D,TFile,TTree,TCanvas, TChain +import sys +import array +import numpy as np + +simulated_files = "mc_samples" +measured_files = "MC_Cutted_files" + +########### Efficiency = measured / simulated ############## + +#-----Run 1-----# + +measured_1 = TChain("DecayTree") +#add all MC cutted files for Run1 +measured_1.Add(os.path.join(measured_files,"Cut_Bu2KpipiJPsi-MM-MC-12145090-2011-MagDown-Sim09k-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +measured_1.Add(os.path.join(measured_files,"Cut_Bu2KpipiJPsi-MM-MC-12145090-2011-MagUp-Sim09k-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +measured_1.Add(os.path.join(measured_files,"Cut_Bu2KpipiJPsi-MM-MC-12145090-2012-MagDown-Sim09k-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +measured_1.Add(os.path.join(measured_files,"Cut_Bu2KpipiJPsi-MM-MC-12145090-2012-MagUp-Sim09k-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) + +simulated_1 = TChain("DecayTree") +#add all MC files for Run1 +simulated_1.Add(os.path.join(simulated_files,"Bu2KpipiJPsi-MM-MC-12145090-2011-MagDown-Sim09k-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +simulated_1.Add(os.path.join(simulated_files,"Bu2KpipiJPsi-MM-MC-12145090-2011-MagUp-Sim09k-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +simulated_1.Add(os.path.join(simulated_files,"Bu2KpipiJPsi-MM-MC-12145090-2012-MagDown-Sim09k-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +simulated_1.Add(os.path.join(simulated_files,"Bu2KpipiJPsi-MM-MC-12145090-2012-MagUp-Sim09k-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) + +N_cut_events_1 = measured_1.Draw("B_MM>>my_histogram1(50,3700,6900)","") +N_events_1 = simulated_1.Draw("B_MM>>my_histogram11(50,3700,6900)","","same") + +# efficiencies computation (estimator for probability) +efficiency_1 = N_cut_events_1/N_events_1 + +# Binomial uncertainty (with estimator P) : sigma = sqrt(P(1-P))/N +Proba_1 = efficiency_1 +uncertainty_1 = np.sqrt(N_events_1*Proba_1*(1-Proba_1))/N_events_1 + +print("---------LHC RUN 1---------") +print("N_events = ", N_events_1, "N_cut_events = ",N_cut_events_1) +print(" Efficiency for Run1 = ",efficiency_1, "±",uncertainty_1) + +#-----Run 2-----# + +measured_2 = TChain("DecayTree") +#add all MC cutted files for Run2 +measured_2.Add(os.path.join(measured_files,"Cut_Bu2KpipiJPsi-MM-MC-12145090-2015-MagDown-Sim09c-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +measured_2.Add(os.path.join(measured_files,"Cut_Bu2KpipiJPsi-MM-MC-12145090-2015-MagUp-Sim09c-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +measured_2.Add(os.path.join(measured_files,"Cut_Bu2KpipiJPsi-MM-MC-12145090-2016-MagDown-Sim09c-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +measured_2.Add(os.path.join(measured_files,"Cut_Bu2KpipiJPsi-MM-MC-12145090-2016-MagUp-Sim09c-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +measured_2.Add(os.path.join(measured_files,"Cut_Bu2KpipiJPsi-MM-MC-12145090-2017-MagDown-Sim09h-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +measured_2.Add(os.path.join(measured_files,"Cut_Bu2KpipiJPsi-MM-MC-12145090-2017-MagUp-Sim09h-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +measured_2.Add(os.path.join(measured_files,"Cut_Bu2KpipiJPsi-MM-MC-12145090-2018-MagDown-Sim09h-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +measured_2.Add(os.path.join(measured_files,"Cut_Bu2KpipiJPsi-MM-MC-12145090-2018-MagUp-Sim09h-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) + + +simulated_2 = TChain("DecayTree") +#add all MC files for Run2 +simulated_2.Add(os.path.join(simulated_files,"Bu2KpipiJPsi-MM-MC-12145090-2015-MagDown-Sim09c-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +simulated_2.Add(os.path.join(simulated_files,"Bu2KpipiJPsi-MM-MC-12145090-2015-MagUp-Sim09c-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +simulated_2.Add(os.path.join(simulated_files,"Bu2KpipiJPsi-MM-MC-12145090-2016-MagDown-Sim09c-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +simulated_2.Add(os.path.join(simulated_files,"Bu2KpipiJPsi-MM-MC-12145090-2016-MagUp-Sim09c-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +simulated_2.Add(os.path.join(simulated_files,"Bu2KpipiJPsi-MM-MC-12145090-2017-MagDown-Sim09h-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +simulated_2.Add(os.path.join(simulated_files,"Bu2KpipiJPsi-MM-MC-12145090-2017-MagUp-Sim09h-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +simulated_2.Add(os.path.join(simulated_files,"Bu2KpipiJPsi-MM-MC-12145090-2018-MagDown-Sim09h-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) +simulated_2.Add(os.path.join(simulated_files,"Bu2KpipiJPsi-MM-MC-12145090-2018-MagUp-Sim09h-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) + + +N_cut_events_2 = measured_2.Draw("B_MM>>my_histogram2(50,3700,6900)","") +N_events_2 = simulated_2.Draw("B_MM>>my_histogram22(50,3700,6900)","","same") + +efficiency_2 = N_cut_events_2/N_events_2 +Proba_2 = efficiency_2 +uncertainty_2 = np.sqrt(N_events_2*Proba_2*(1-Proba_2))/N_events_2 + +print("---------LHC RUN 2---------") +print("N_events = ", N_events_2, "N_cut_events = ",N_cut_events_2) +#print("Mean mu = ", Proba_2) +#print("sigma^2 = ",Proba_2*(1-Proba_2)) +print(" Efficiency for Run2 = ",efficiency_2, "±",uncertainty_2) \ No newline at end of file diff --git a/efficiency_BDT.py b/efficiency_BDT.py new file mode 100644 index 0000000..7e1ae5f --- /dev/null +++ b/efficiency_BDT.py @@ -0,0 +1,166 @@ +from cmath import sqrt +import os,sys +import ROOT +from ROOT import TH1D,TH2D,TFile,TTree,TCanvas, TChain +import sys +import array +import numpy as np +import matplotlib.pyplot as plt +import scipy + +ROOT.gROOT.SetBatch(1) + +## MAIN FUNCTION ## +def main(bins): #BDT_threshold + + #BDT_threshold = sys.argv [1] + bins = sys.argv[1] + + simulated_files = "mc_samples" + measured_files = "BDT_MC_Cutted_files" + + ########### Efficiency = measured / simulated ############## + + #-----Run 1-----# + + BDT_files = TChain("DecayTree") + + BDT_files.Add(os.path.join(measured_files,"BDT1_Cut_Bu2KpipiJPsi-MM-MC-12145090-2018-MagDown-Sim09h-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) + BDT_files.Add(os.path.join(measured_files,"BDT1_Cut_Bu2KpipiJPsi-MM-MC-12145090-2018-MagUp-Sim09h-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) + + + Raw_files = TChain("DecayTree") + + Raw_files.Add(os.path.join(simulated_files,"Bu2KpipiJPsi-MM-MC-12145090-2018-MagDown-Sim09h-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) + Raw_files.Add(os.path.join(simulated_files,"Bu2KpipiJPsi-MM-MC-12145090-2018-MagUp-Sim09h-StrippingBu2LLK_MCTruthMM_selected_resampled.root")) + + N_cut = BDT_files.Draw("B_MM>>histo_cut"," BDT1 > 0.1 & B_MM <5200") # + str(BDT_threshold) + N_events_1 = Raw_files.Draw("B_MM>>histo_raw"," B_MM <5200","same") + + # efficiencies computation (estimator for probability) + efficiency_1 = N_cut/N_events_1 + + # Binomial uncertainty (with estimator P) : sigma = sqrt(P(1-P))/N + Proba_1 = efficiency_1 + uncertainty_1 = np.sqrt(N_events_1*Proba_1*(1-Proba_1))/N_events_1 + + print("---------LHC 2018---------") + print("N_events = ", N_events_1, "N_cut_events = ",N_cut) + print(" Efficiency for 2018 = ",efficiency_1, "±",uncertainty_1) + + N_bins = int(bins) + BDT = np.empty(N_bins, dtype=object) + efficiency_B = np.empty(N_bins, dtype=object) + efficiency_S = np.empty(N_bins, dtype=object) + N_cut_B = np.empty(N_bins, dtype=object) + N_cut_S = np.empty(N_bins, dtype=object) + Sqrt_B = np.empty(N_bins, dtype=object) + + yerr = np.empty(N_bins, dtype=object) + mass = np.empty(N_bins, dtype=object) + + for i in range(0,N_bins): + + BDT[i] = -1 + i*2/N_bins + #---------background cut without the peak + + # BDT Cutted data | Background window + N_cut_BDT_i_1 = BDT_files.Draw("B_MM>>histo_cut"," BDT1 > " + str(BDT[i]) + "& B_MM<5200") # & B_MM>5400 + N_cut_BDT_i_2 = BDT_files.Draw("B_MM>>histo_cut"," BDT1 > " + str(BDT[i]) + "& B_MM>5400") # & B_MM>5400 + N_cut_BDT_i = N_cut_BDT_i_1 + N_cut_BDT_i_2 + + # Raw data | Background window + N_event_BDT_i_1 = Raw_files.Draw("B_MM>>histo_raw","B_MM<5200","same") + N_event_BDT_i_2 = Raw_files.Draw("B_MM>>histo_raw","B_MM>5400","same") + + N_event_BDT_i = N_event_BDT_i_1 + N_event_BDT_i_2 + + #---------signal cut only the peak + + # BDT Cutted data | Signal window + N_S_cut_i = BDT_files.Draw("B_MM>>histo_cut"," BDT1 > " + str(BDT[i]) + "& 5200>histo_cut", "5200