Homec4science
Diffusion TP4b (master) Browse

R12329/efficiency_BDT.py
71d23c230be9master

Last Change
Raw File

efficiency_BDT.py
View Options

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<B_MM && B_MM<5400")
# BDT Raw data | Signal window
N_S_event_i = BDT_files.Draw("B_MM>>histo_cut", "5200<B_MM && B_MM<5400")
efficiency_B[i] = N_cut_BDT_i / N_event_BDT_i
efficiency_S[i] = N_S_cut_i / N_S_event_i
N_cut_B[i] = N_cut_BDT_i
N_cut_S[i] = N_S_cut_i
Sqrt_B[i] = np.sqrt(N_cut_BDT_i)
yerr[i] = 0.1
#plt.figure()
#plt.plot(BDT,N_cut_B,label=)
#plt.plot(BDT,N_cut_S)
#plt.xlabel("BDT threshold ")
#plt.ylabel("Nevents cutted background vs signal")
#plt.legend(fontsize=15, loc="lower left")
#plt.show()
# Plot error bar
plt.errorbar(BDT, N_cut_B, yerr=yerr, fmt = 'o',color = 'navy',
ecolor = 'cornflowerblue', elinewidth = 10, capsize=10,label='1$\sigma$')
plt.errorbar(BDT, N_cut_S, yerr, fmt = 'o',color = 'navy',
ecolor = 'crimson', elinewidth = 2, capsize=12,label='2$\sigma$')
# Display graph
plt.rc('ytick', labelsize=15)
plt.rc('xtick', labelsize=15)
plt.legend(fontsize=15,loc="lower left")
#plt.savefig('./plots/result_H0.png')
plt.show()
#================================= FITTING =================================#
#plt.figure()
#plt.plot(BDT,N_cut_B)
#plt.plot(BDT,N_cut_S)
#plt.xlabel("BDT threshold ")
#plt.ylabel("Nevents cutted background vs signal")
#plt.legend
#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()
#fig, ax = plt.subplots()
#ax.plot(BDT,N_cut_B, label='Background ')
#ax.plot(BDT,N_cut_S, label='Signal')
#ax.set_xlabel('BDT threshold') # Add an x-label to the axes.
#ax.set_ylabel('Nevents cutted background vs signal') # Add a y-label to the axes.
# ax.set_title("Simple Plot") # Add a title to the axes.
# ax.legend(); # Add a legend.
#plt.style.use('_mpl-gallery')
## plot:
#fig, ax = plt.subplots()
#ax.errorbar(BDT, N_cut_B, yerr, fmt='o', linewidth=2, capsize=6)
#ax.errorbar(BDT, N_cut_S, yerr, fmt='o', linewidth=2, capsize=6)
#ax.set(xlim=(0, 8), xticks=np.arange(1, 8),ylim=(0, 8), yticks=np.arange(1, 8))
#plt.show()
#end function
if(__name__=='__main__'):
sys.exit(main(sys.argv[1]))
c4science · Help