-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathFunctions.py
More file actions
190 lines (154 loc) · 8.53 KB
/
Functions.py
File metadata and controls
190 lines (154 loc) · 8.53 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# Known result Predicted results Want to plot? For Multianalysis for pdf save
def predictionY_interval(binsTab, y_train, y_test, y_pred_onTrain, y_pred_onTest, plotResult, computeInterval, pdfName='yPred'): #Function form to use in Multianalysis
interval=yTrainSHist=yTrainBHist=yTestSHist=yTestBHist=0 # For the return if not computed
yTestS=y_pred_onTest[y_test==1]
yTestB=y_pred_onTest[y_test==0]
if plotResult:
yTrainS=y_pred_onTrain[y_train==1] # Take values at index depending of BDT y output
yTrainB=y_pred_onTrain[y_train==0]
fig, ax = plt.subplots(figsize=(8, 4), dpi=300) # Want to stack histograms on top of same axes
ax.hist(yTrainS, bins=binsTab, color='r',histtype='bar',density=True, alpha=0.3) # Adding transparency for histogram columns
ax.hist(yTrainB, bins=binsTab, color='b', histtype='bar',density=True, alpha=0.3) # collect the results in varibles for later KS test
yTrainSHist = ax.hist(yTrainS, bins=binsTab, color='r', histtype='step',density=True, label='Signal (train)') # Histogram values
yTrainBHist = ax.hist(yTrainB, bins=binsTab, color='b', histtype='step',density=True, label='Background (train)')
yTestSHist = plt.hist(yTestS, bins=binsTab, density=True,alpha = 0.0) # Invisible to collect histogram values
yTestBHist = plt.hist(yTestB, bins=binsTab, density=True,alpha = 0.0)
bin_centers = 0.5*(binsTab[1:] + binsTab[:-1])
ax.scatter(bin_centers, yTestSHist[0], marker='o', c='r', s=20, alpha=1,label='Signal (test)') # Display as dot in histogram
ax.scatter(bin_centers, yTestBHist[0], marker='o', c='b', s=20, alpha=1,label='Background (test)') # Must take first array for data, second one is bins values
plt.xlabel('BDT signal response')
plt.ylabel('Normalized number of events')
plt.legend()
plt.savefig(f'plots/{pdfName}.pdf')
plt.close()
if computeInterval:
numberBins=len(binsTab)
errorMax=2*(1/(numberBins-1))
binsTab=np.linspace(0, 1, num=numberBins)
yTestSHist=np.histogram(yTestS,bins=binsTab,density=False) #Density=True to normalise, not necessary here we just want the max
yTestBHist=np.histogram(yTestB,bins=binsTab,density=False)
maxS=np.argmax(yTestSHist[0]) #first tab is the histogram values, second is the bins values on axis
maxB=np.argmax(yTestBHist[0])
interval=abs(binsTab[maxS]-binsTab[maxB])
#print(f'Interval between Signal and backgroung max in histogram: {interval} pm {errorMax}')
return interval,yTrainSHist,yTrainBHist,yTestSHist,yTestBHist
def bgFromExp(blindMin,blindMax,massWidth,paramExpA,paramExpB):
return (1/massWidth)*( (paramExpA/paramExpB)*(np.exp(-paramExpB*blindMin)-np.exp(-paramExpB*blindMax)) )
def monoExp(x, a, b):
return (a * np.exp(-b * x) )
def fitBackground(XselectedBg,printPlotFit):
# C,D,E Normalisation factors for easier fit:
C=100
E=5500
bgMin=5500
bgMax=6500
steps=21
massTab=np.linspace(bgMin,bgMax,steps) #X axis
binCenters=[(0.5*(massTab[i]+massTab[i+1])-E)/C for i in range (len(massTab)-1)] # X axis transformaion for easier fit
massB0=XselectedBg['B_s0_DTF_M']
massHist,binEdges=np.histogram(massB0,bins=massTab,density=False) # Histogram computation
D=massHist[0]+1 # +1 to avoid dividing by 0
massHist=massHist/D # Histogram normilised to be ~1 at x=0
params, paramCov = curve_fit(monoExp, binCenters, massHist,bounds=([0.0,0.0], [2.0, 2.0])) # Do the easier exponential fit
a, b=params
A=D*float(a)*np.exp(float(b)*E/C)
B=float(b)/C
dA=np.sqrt(paramCov[0][0])*D*np.exp(B*E/C)
dB=np.sqrt(paramCov[1][1])/C
#print(f'a: {a}, b: {b}, c:{c}')
print(f'paramCov: {paramCov}')
#dA, dB, dC=paramCov # Estimated covarience on parameters
if printPlotFit:
fig, ax = plt.subplots(figsize=(8, 4), dpi=300)
ax.plot(binCenters,massHist,'k.',label=f'Mass Histogram')
ax.plot(binCenters,monoExp(np.array(binCenters),*params),'b--', # * in call to expands the tuple into separate elements
label=f'fit: N={"{0:.4f}".format(a)}*exp(-{"{0:.4f}".format(b)}*M)')
plt.xlabel('B_0 mass')
plt.ylabel('Number of events')
plt.legend()
plt.show()
plt.close()
blindingMin=5100
blindingMax=5500
nBgOutBlinding=len(massB0)
massWidth=(bgMax-bgMin)/(steps-1)
nBgInBlinding=bgFromExp(blindMin=blindingMin,blindMax=blindingMax,massWidth=massWidth,paramExpA=A,paramExpB=B)#+ c*(blindingMax-blindingMin) )
nBgcheck1=bgFromExp(blindMin=blindingMin,blindMax=blindingMax,massWidth=massWidth,paramExpA=A+dA,paramExpB=B+dB)
nBgcheck2=bgFromExp(blindMin=blindingMin,blindMax=blindingMax,massWidth=massWidth,paramExpA=A-dA,paramExpB=B-dB)
upperNBg=bgFromExp(blindMin=blindingMin,blindMax=blindingMax,massWidth=massWidth,paramExpA=A+dA,paramExpB=B-dB)
lowerNBg=bgFromExp(blindMin=blindingMin,blindMax=blindingMax,massWidth=massWidth,paramExpA=A-dA,paramExpB=B+dB)
#print(f'nBgInBlinding:{nBgInBlinding}')
#print(f'upperNBg:{upperNBg}')
#print(f'lowerNBg:{lowerNBg}')
#print(f'nBgcheck1:{nBgcheck1}')
#print(f'nBgcheck2:{nBgcheck2}')
if( (upperNBg<nBgcheck1)|(upperNBg<nBgcheck2)|(lowerNBg>nBgcheck1)|(lowerNBg>nBgcheck2) ):
print("ERROR in uncertainty computation, min/max are not extremums.")
#print(f'Rapport {nBgInBlinding/nBgOutBlinding} between Inside {nBgInBlinding} and outside bliding: {nBgOutBlinding}')
if nBgOutBlinding < 10:
nBgInBlinding=100000
print('Too few event to fit background (<10)')
errNBg=np.maximum(abs(upperNBg-nBgInBlinding),abs(nBgInBlinding-lowerNBg))
#print(f'Background: {nBgInBlinding} pm {errNBg}')
return nBgInBlinding,errNBg
# Try without scaling function to reduce uncertainty
def fitBackgroundTry(XselectedBg,printPlotFit):
# C,D,E Normalisation factors for easier fit:
C=100
E=5500
bgMin=5500
bgMax=6500
steps=21
massTab=np.linspace(bgMin,bgMax,steps) #X axis
binCenters=[(0.5*(massTab[i]+massTab[i+1])-E)/C for i in range (len(massTab)-1)] # X axis transformaion for easier fit
massB0=XselectedBg['B_s0_DTF_M']
massHist,binEdges=np.histogram(massB0,bins=massTab,density=False) # Histogram computation
D=massHist[0]+1 # +1 to avoid dividing by 0
massHist=massHist/D # Histogram normilised to be ~1 at x=0
params, paramCov = curve_fit(monoExp, binCenters, massHist,bounds=([0.0,0.0], [2.0, 2.0])) # Do the easier exponential fit
a, b=params
A=D*float(a)#*np.exp(float(b)*E/C)
B=float(b)#/C
dA=np.sqrt(paramCov[0][0])#*D*np.exp(B*E/C)
dB=np.sqrt(paramCov[1][1])#/C
#print(f'a: {a}, b: {b}, c:{c}')
#print(f'paramCov: {paramCov}')
#dA, dB, dC=paramCov # Estimated covarience on parameters
if printPlotFit:
fig, ax = plt.subplots(figsize=(8, 4), dpi=300)
ax.plot(binCenters,massHist,'k.',label=f'Mass Histogram')
ax.plot(binCenters,monoExp(np.array(binCenters),*params),'b--', # * in call to expands the tuple into separate elements
label=f'fit: N={"{0:.4f}".format(a)}*exp(-{"{0:.4f}".format(b)}*M)')
plt.xlabel('B_0 mass')
plt.ylabel('Number of events')
plt.legend()
plt.show()
plt.close()
blindingMin=(5100-E)/C
blindingMax=(5500-E)/C
nBgOutBlinding=len(massB0)
massWidth=(bgMax-bgMin)/(steps-1)/C
nBgInBlinding=bgFromExp(blindMin=blindingMin,blindMax=blindingMax,massWidth=massWidth,paramExpA=A,paramExpB=B)#+ c*(blindingMax-blindingMin) )
nBgExtremum1=bgFromExp(blindMin=blindingMin,blindMax=blindingMax,massWidth=massWidth,paramExpA=A+dA,paramExpB=B+dB)
nBgExtremum2=bgFromExp(blindMin=blindingMin,blindMax=blindingMax,massWidth=massWidth,paramExpA=A-dA,paramExpB=B-dB)
nBgExtremum3=bgFromExp(blindMin=blindingMin,blindMax=blindingMax,massWidth=massWidth,paramExpA=A+dA,paramExpB=B-dB)
nBgExtremum4=bgFromExp(blindMin=blindingMin,blindMax=blindingMax,massWidth=massWidth,paramExpA=A-dA,paramExpB=B+dB)
print(f'nBgInBlinding:{nBgInBlinding}')
print(f'nBgExtremum1:{nBgExtremum1}')
print(f'nBgExtremum2:{nBgExtremum2}')
print(f'nBgExtremum3:{nBgExtremum3}')
print(f'nBgExtremum4:{nBgExtremum4}')
#if( (upperNBg<nBgcheck1)|(upperNBg<nBgcheck2)|(lowerNBg>nBgcheck1)|(lowerNBg>nBgcheck2) ):
# print("ERROR in uncertainty computation, min/max are not extremums.")
#print(f'Rapport {nBgInBlinding/nBgOutBlinding} between Inside {nBgInBlinding} and outside bliding: {nBgOutBlinding}')
if nBgOutBlinding < 10:
nBgInBlinding=100000
print('Too few event to fit background (<10)')
upperNBg=np.amax([nBgExtremum1,nBgExtremum2,nBgExtremum3,nBgExtremum4])
lowerNBg=np.amin([nBgExtremum1,nBgExtremum2,nBgExtremum3,nBgExtremum4])
errNBg=np.maximum(abs(upperNBg-nBgInBlinding),abs(nBgInBlinding-lowerNBg))
print(f'Background: {nBgInBlinding} pm {errNBg}')
return nBgInBlinding,errNBg