Homec4science
Diffusion TP4b (master) Browse

R12329/Final_plots.py
71d23c230be9master

Last Change
Raw File

Final_plots.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):
#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<B_MM && B_MM<5350")
N_S_BDT_MC[i] = N_S_BDT_MC_i
N_S_MC[i] = N_S_MC_i
efficiency_S_MC[i] = N_S_BDT_MC[i] / N_S_MC[i]
N_S_err[i] = np.sqrt(N_S_BDT_MC_i)
Efficiency_S_err[i] = (np.sqrt(N_S_BDT_MC_i))/N_S_MC_i +(N_S_BDT_MC_i*np.sqrt(N_S_MC_i))/((N_S_MC_i)**2)
#print("Efficiency_S_err[i] = ", Efficiency_S_err[i])
# Plot Background VS Signal
plt.errorbar(BDT, N_BG_BDT,xerr=None, yerr=N_BG_err, fmt = '.',color = 'navy',
ecolor = 'cornflowerblue', elinewidth = 1, capsize=1,label='BDT Cutted Muons')
plt.errorbar(BDT, N_S_BDT_MC,xerr=None,yerr=N_S_err, fmt = '.',color = 'red',
ecolor = 'crimson', elinewidth = 1, capsize=1,label='BDT cutted MC Signal')
# Display graph
plt.rc('ytick', labelsize=15)
plt.rc('xtick', labelsize=15)
plt.xlabel('BDT threshold')
plt.ylabel('Events')
plt.legend(fontsize=15,loc="lower left")
plt.show()
# Plot Efficiency Signal / Background
plt.errorbar(BDT, efficiency_B, xerr =None,yerr=Efficiency_B_err, fmt = '.',color = 'navy',
ecolor = 'cornflowerblue', elinewidth =1, capsize=1,label='Efficiency Background')
plt.errorbar(BDT, efficiency_S_MC, xerr=None,yerr=Efficiency_S_err, fmt = '.',color = 'red',
ecolor = 'crimson', elinewidth =1, capsize=1,label='Efficiency MC signal')
# Display graph
plt.rc('ytick', labelsize=15)
plt.rc('xtick', labelsize=15)
plt.xlabel('BDT threshold')
plt.ylabel('Efficiency')
plt.legend(fontsize=15,loc="lower left")
plt.show()
#================================= FITTING =================================#
'''
# FIT(mass) = A•exp(B•mass)
#A= 1.43828*10**1 # ± 1.08816e-2
#B= -1.12837*10**-3 # ± 1.91257e-6
M_min = 5200 #MeV
M_max = 5350
#B_new = (A/B)*( np.exp(B*M_max)- np.exp(B*M_min) )
#N _BG_BDT_PEAK = Cutted_files_Data.Draw("B_MM>>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]))
c4science · Help