Belle II Software light-2406-ragdoll
B2JpsiKs_mu_qrBelleDataSplot.py
1#!/usr/bin/env python
2
3
10
11# ************************************ Flavor Tagging ***********************************
12# * *
13# * This validation script performs an Splot analysis of Belle Data using mbc as *
14# * Splot variable in order to extract the qr distribution of signal B0s. *
15# * The signal component of Mbc is first fitted on converted Belle MC using a *
16# * 2 Gaussian function. The Mbc component is used together with an Argus function *
17# * to fit the mbc distribution of converted Belle data. *
18# * An Splot is performed and its output variable used to weight the qr *
19# * distribution. *
20# * For more information see Sec. 4.8 in BELLE2-PTHESIS-2018-003 *
21# * Usage: *
22# * basf2 B2JpsiKs_mu_qrBelleDataSplot.py RootNtupleFilesMC RootNtupleFilesData treeName *
23# * *
24# ******************************************************************************************
25
26import math
27import glob
28import sys
29from array import array
30import matplotlib.pyplot as plt
31from defaultEvaluationParameters import r_size, r_subsample, rbins
32import ROOT
33import numpy as np
34import matplotlib as mpl
35mpl.use('Agg')
36mpl.rcParams['text.usetex'] = True
37
38
39PATH = "."
40
41workingFilesMC = str(sys.argv[1])
42workingFilesData = str(sys.argv[2])
43
44workingFilesMC = glob.glob(str(workingFilesMC))
45workingFilesData = glob.glob(str(workingFilesData))
46
47lim = 10
48limZ = 0.03
49limZsig = 0.05
50
51treenames = [str(sys.argv[3])]
52
53mc_total_notTagged = 0
54
55data_total_notTagged = 0
56
57data_total_Tagged = 0
58
59ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.WARNING)
60
61
62def efficiencyCalculator(mc, data, method):
63
64 # Calculate Efficiency for MC
65
66 mc_dataB0 = mc[np.where(mc[:, 1] > 0)]
67 mc_dataB0bar = mc[np.where(mc[:, 1] < 0)]
68
69 mc_totaldataB0 = np.absolute(mc_dataB0[:, 0])
70 mc_totaldataB0bar = np.absolute(mc_dataB0bar[:, 0])
71
72 mc_rvalueB0 = (np.histogram(mc_totaldataB0, rbins, weights=mc_totaldataB0)[0] / np.histogram(mc_totaldataB0, rbins)[0])
73 mc_rvalueMSB0 = (
74 np.histogram(
75 mc_totaldataB0,
76 rbins,
77 weights=mc_totaldataB0 *
78 mc_totaldataB0)[0] /
79 np.histogram(
80 mc_totaldataB0,
81 rbins)[0])
82 mc_rvalueStdB0 = np.sqrt(abs(mc_rvalueMSB0 - mc_rvalueB0 * mc_rvalueB0))
83 mc_entriesB0 = np.histogram(mc_totaldataB0, rbins)[0]
84
85 mc_rvalueB0bar = (
86 np.histogram(
87 mc_totaldataB0bar,
88 rbins,
89 weights=mc_totaldataB0bar)[0] /
90 np.histogram(
91 mc_totaldataB0bar,
92 rbins)[0])
93 mc_rvalueMSB0bar = (
94 np.histogram(
95 mc_totaldataB0bar,
96 rbins,
97 weights=mc_totaldataB0bar *
98 mc_totaldataB0bar)[0] /
99 np.histogram(
100 mc_totaldataB0bar,
101 rbins)[0])
102 mc_rvalueStdB0bar = np.sqrt(abs(mc_rvalueMSB0bar - mc_rvalueB0bar * mc_rvalueB0bar))
103 mc_entriesB0bar = np.histogram(mc_totaldataB0bar, rbins)[0]
104
105 mc_total_entriesB0 = mc_totaldataB0.shape[0] + mc_total_notTagged / 2
106 mc_total_tagged_B0 = mc_totaldataB0.shape[0]
107 mc_event_fractionB0 = mc_entriesB0.astype(float) / mc_total_entriesB0
108
109 mc_total_entriesB0bar = mc_totaldataB0bar.shape[0] + mc_total_notTagged / 2
110 mc_total_tagged_B0bar = mc_totaldataB0bar.shape[0]
111 mc_event_fractionB0bar = mc_entriesB0bar.astype(float) / mc_total_entriesB0bar
112
113 mc_total_entries = mc_total_entriesB0 + mc_total_entriesB0bar
114 mc_total_tagged = mc_totaldataB0.shape[0] + mc_totaldataB0bar.shape[0]
115 mc_event_fractionTotal = (mc_entriesB0.astype(float) + mc_entriesB0bar.astype(float)) / mc_total_entries
116
117 mc_tagging_eff = mc_total_tagged / (mc_total_tagged + mc_total_notTagged)
118 mc_DeltaTagging_eff = math.sqrt(mc_total_tagged * mc_total_notTagged**2 +
119 mc_total_notTagged * mc_total_tagged**2) / (mc_total_entries**2)
120 mc_tot_eff_effB0 = 0
121 mc_tot_eff_effB0bar = 0
122 mc_uncertainty_eff_effB0 = 0
123 mc_uncertainty_eff_effB0bar = 0
124 mc_uncertainty_eff_effAverage = 0
125 mc_diff_eff_Uncertainty = 0
126 mc_event_fractionDiff = array('f', [0] * r_size)
127 mc_rvalueB0Average = array('f', [0] * r_size)
128 mc_rvalueStdB0Average = array('f', [0] * r_size)
129 mc_wvalue = array('f', [0] * r_size)
130 mc_wvalueUncertainty = array('f', [0] * r_size)
131 mc_wvalueB0 = array('f', [0] * r_size)
132 mc_wvalueB0bar = array('f', [0] * r_size)
133 mc_wvalueDiff = array('f', [0] * r_size)
134 mc_wvalueDiffUncertainty = array('f', [0] * r_size)
135 mc_iEffEfficiency = array('f', [0] * r_size)
136 mc_iEffEfficiencyUncertainty = array('f', [0] * r_size)
137 mc_iEffEfficiencyB0Uncertainty = array('f', [0] * r_size)
138 mc_iEffEfficiencyB0barUncertainty = array('f', [0] * r_size)
139 mc_iDeltaEffEfficiency = array('f', [0] * r_size)
140 mc_iDeltaEffEfficiencyUncertainty = array('f', [0] * r_size)
141
142 print('****************** MEASURED EFFECTIVE EFFICIENCY FOR COMBINER USING ' + method + ' ***************************')
143 print('* *')
144 print('* --> DETERMINATION BASED ON MONTE CARLO *')
145 print('* *')
146 print('* ------------------------------------------------------------------------------------------------ *')
147 print('* r-interval <r> Efficiency Delta_Effcy w Delta_w *')
148 print('* ------------------------------------------------------------------------------------------------ *')
149
150 for i in range(0, r_size - 1):
151 mc_rvalueB0Average[i] = (mc_rvalueB0[i] + mc_rvalueB0bar[i]) / 2
152 mc_rvalueStdB0Average[i] = math.sqrt(mc_rvalueStdB0[i]**2 / (mc_entriesB0[0] - 1) +
153 mc_rvalueStdB0bar[i]**2 / (mc_entriesB0bar[0] - 1)) / 2
154 mc_wvalueB0[i] = (1 - mc_rvalueB0[i]) / 2
155 mc_wvalueB0bar[i] = (1 - mc_rvalueB0bar[i]) / 2
156 mc_wvalueDiff[i] = mc_wvalueB0[i] - mc_wvalueB0bar[i]
157 mc_wvalueDiffUncertainty[i] = math.sqrt((mc_rvalueStdB0[i] / 2)**2 / (mc_entriesB0[0] - 1) +
158 (mc_rvalueStdB0bar[i] / 2)**2 / (mc_entriesB0bar[0] - 1))
159 mc_wvalue[i] = (mc_wvalueB0[i] + mc_wvalueB0bar[i]) / 2
160 mc_wvalueUncertainty[i] = mc_wvalueDiffUncertainty[i] / 2
161 # fraction of events/all events
162
163 mc_event_fractionDiff[i] = (mc_entriesB0[i] - mc_entriesB0bar[i]) / mc_total_entries
164
165 mc_iEffEfficiency[i] = mc_tagging_eff * (mc_event_fractionB0[i] * mc_rvalueB0[i] * mc_rvalueB0[i] +
166 mc_event_fractionB0bar[i] * mc_rvalueB0bar[i] * mc_rvalueB0bar[i]) / 2
167
168 mc_iEffEfficiencyB0Uncertainty[i] = mc_rvalueB0[i] * \
169 math.sqrt((2 * mc_total_entriesB0 * mc_entriesB0[i] * mc_rvalueStdB0[i])**2 / (mc_entriesB0[0] - 1) +
170 mc_rvalueB0[i]**2 * mc_entriesB0[i] *
171 (mc_total_entriesB0 * (mc_total_entriesB0 - mc_entriesB0[i]) +
172 mc_entriesB0[i] * mc_total_notTagged)) / (mc_total_entriesB0**2)
173 mc_iEffEfficiencyB0barUncertainty[i] = mc_rvalueB0bar[i] * \
174 math.sqrt((2 * mc_total_entriesB0bar * mc_entriesB0bar[i] * mc_rvalueStdB0bar[i])**2 / (mc_entriesB0bar[0] - 1) +
175 mc_rvalueB0bar[i]**2 * mc_entriesB0bar[i] *
176 (mc_total_entriesB0bar * (mc_total_entriesB0bar - mc_entriesB0bar[i]) +
177 mc_entriesB0bar[i] * mc_total_notTagged)) / (mc_total_entriesB0bar**2)
178
179 mc_iDeltaEffEfficiency[i] = mc_event_fractionB0[i] * mc_rvalueB0[i] * \
180 mc_rvalueB0[i] - mc_event_fractionB0bar[i] * mc_rvalueB0bar[i] * mc_rvalueB0bar[i]
181
182 mc_iDeltaEffEfficiencyUncertainty[i] = math.sqrt(
183 mc_iEffEfficiencyB0Uncertainty[i]**2 +
184 mc_iEffEfficiencyB0barUncertainty[i]**2)
185
186 mc_iEffEfficiencyUncertainty[i] = mc_iDeltaEffEfficiencyUncertainty[i] / 2
187
188 mc_diff_eff_Uncertainty = mc_diff_eff_Uncertainty + mc_iDeltaEffEfficiencyUncertainty[i]**2
189
190 # finally calculating the total effective efficiency
191 mc_tot_eff_effB0 = mc_tot_eff_effB0 + mc_event_fractionB0[i] * mc_rvalueB0[i] * mc_rvalueB0[i]
192 mc_tot_eff_effB0bar = mc_tot_eff_effB0bar + mc_event_fractionB0bar[i] * mc_rvalueB0bar[i] * mc_rvalueB0bar[i]
193 mc_uncertainty_eff_effAverage = mc_uncertainty_eff_effAverage + mc_iEffEfficiencyUncertainty[i]**2
194 mc_uncertainty_eff_effB0 = mc_uncertainty_eff_effB0 + mc_iEffEfficiencyB0Uncertainty[i]**2
195 mc_uncertainty_eff_effB0bar = mc_uncertainty_eff_effB0bar + mc_iEffEfficiencyB0barUncertainty[i]**2
196
197 # intervallEff[i] = event_fractionTotal[i] * rvalueB0Average[i] * rvalueB0Average[i]
198 print('* ' + f'{r_subsample[i]:.3f}' + ' - ' + f'{r_subsample[i + 1]:.3f}' + ' ' +
199 f'{mc_rvalueB0Average[i]:.3f}' + ' +- ' + f'{mc_rvalueStdB0Average[i]:.4f}' + ' ' +
200 f'{mc_event_fractionTotal[i]:.4f}' + ' ' +
201 f'{mc_event_fractionDiff[i]: .4f}' + ' ' +
202 f'{mc_wvalue[i]:.4f}' + ' +- ' + f'{mc_wvalueUncertainty[i]:.4f}' + ' ' +
203 f'{mc_wvalueDiff[i]: .4f}' + ' +- ' + f'{mc_wvalueDiffUncertainty[i]:.4f}' + ' *')
204
205 mc_average_eff_eff = (mc_tot_eff_effB0 + mc_tot_eff_effB0bar) / 2
206 mc_uncertainty_eff_effAverage = math.sqrt(mc_uncertainty_eff_effAverage)
207 mc_uncertainty_eff_effB0 = math.sqrt(mc_uncertainty_eff_effB0)
208 mc_uncertainty_eff_effB0bar = math.sqrt(mc_uncertainty_eff_effB0bar)
209 mc_diff_eff = mc_tot_eff_effB0 - mc_tot_eff_effB0bar
210 mc_diff_eff_Uncertainty = math.sqrt(mc_diff_eff_Uncertainty)
211
212 print('* ------------------------------------------------------------------------------------------------ *')
213 print('* *')
214 print('* __________________________________________________________________________________________ *')
215 print('* | | *')
216 print('* | TOTAL NUMBER OF TAGGED EVENTS = ' +
217 f"{f'{mc_total_tagged:.0f}':<24}" + f"{'| *':>36}")
218 print('* | | *')
219 print(
220 '* | TOTAL AVERAGE EFFICIENCY (q=+-1)= ' +
221 f'{mc_tagging_eff * 100:.2f}' +
222 " +- " +
223 f'{mc_DeltaTagging_eff * 100:.2f}' +
224 ' % | *')
225 print('* | | *')
226 print(
227 '* | TOTAL AVERAGE EFFECTIVE EFFICIENCY (q=+-1)= ' +
228 f'{mc_average_eff_eff * 100:.6f}' +
229 " +- " +
230 f'{mc_uncertainty_eff_effAverage * 100:.6f}' +
231 ' % | *')
232 print('* | | *')
233 print(
234 '* | TOTAL AVERAGE EFFECTIVE EFFICIENCY ASYMMETRY (q=+-1)= ' +
235 f'{mc_diff_eff * 100:^9.6f}' +
236 " +- " +
237 f'{mc_diff_eff_Uncertainty * 100:.6f}' +
238 ' % | *')
239 print('* | | *')
240 print('* | B0-TAGGER TOTAL EFFECTIVE EFFICIENCIES: ' +
241 f'{mc_tot_eff_effB0 * 100:.2f}' + " +-" + f'{mc_uncertainty_eff_effB0 * 100: 4.2f}' +
242 ' % (q=+1) ' +
243 f'{mc_tot_eff_effB0bar * 100:.2f}' + " +-" + f'{mc_uncertainty_eff_effB0bar * 100: 4.2f}' +
244 ' % (q=-1) ' + ' | *')
245 print('* | | *')
246 print('* | FLAVOR PERCENTAGE (MC): ' +
247 f'{mc_total_tagged_B0 / mc_total_tagged * 100:.2f}' + ' % (q=+1) ' +
248 f'{mc_total_tagged_B0bar / mc_total_tagged * 100:.2f}' + ' % (q=-1) Diff=' +
249 f'{(mc_total_tagged_B0 - mc_total_tagged_B0bar) / mc_total_tagged * 100:^5.2f}' + ' % | *')
250 print('* |__________________________________________________________________________________________| *')
251 print('* *')
252 print('****************************************************************************************************')
253
254 print(r'\begin{tabular}{ l r r r r r r r }')
255 print(r'\hline')
256 print(r'$r$- Interval & $\varepsilon_i\ $ & $w_i \pm \delta w_i\enskip\, $ ' +
257 r' & $\Delta w_i \pm \delta\Delta w_i $& $\varepsilon_{\text{eff}, i} \pm \delta\varepsilon_{\text{eff}, i}\enskip$ ' +
258 r' & & & $\Delta \varepsilon_{\text{eff}, i} \pm \delta\Delta \varepsilon_{\text{eff}, i} $\\ \hline\hline')
259 for i in range(0, r_size - 1):
260 print('$ ' + f'{r_subsample[i]:.3f}' + ' - ' + f'{r_subsample[i + 1]:.3f}' + '$ & $' +
261 f'{mc_event_fractionTotal[i] * 100: 6.1f}' + '$ & $' +
262 f'{mc_wvalue[i] * 100: 7.3f}' + r" \pm " + f'{mc_wvalueUncertainty[i] * 100:2.3f}' + r' $ & $' +
263 f'{mc_wvalueDiff[i] * 100: 6.1f}' + r" \pm " +
264 f'{mc_wvalueDiffUncertainty[i] * 100:2.3f}' + r'\, $ & $' +
265 f'{mc_iEffEfficiency[i] * 100: 8.4f}' + # + '$ & $' +
266 r" \pm " + f'{mc_iEffEfficiencyUncertainty[i] * 100:2.4f}' + r'\, $ & & & $' +
267 f'{mc_iDeltaEffEfficiency[i] * 100: 6.4f}' + # +
268 r" \pm " + f'{mc_iDeltaEffEfficiencyUncertainty[i] * 100:2.4f}' +
269 r'\enskip $ \\ ')
270 print(r'\hline\hline')
271 print(
272 r'\multicolumn{1}{r}{Total} & & \multicolumn{3}{r}{ $\varepsilon_\text{eff} = ' +
273 r'\sum_i \varepsilon_i \cdot \langle 1-2w_i\rangle^2 = ' +
274 f'{mc_average_eff_eff * 100: 6.2f}' +
275 r" \pm " +
276 f'{mc_uncertainty_eff_effAverage * 100: 6.2f}' +
277 r'\quad\,$ }')
278 print(r'& & \multicolumn{2}{r}{ $\Delta \varepsilon_\text{eff} = ' +
279 f'{mc_diff_eff * 100: 6.2f}' + r" \pm " + f'{mc_diff_eff_Uncertainty * 100: 6.2f}' + r'\ \quad\, $ }' +
280 r' \\')
281 print(r'\hline')
282 print(r'\end{tabular}')
283
284 # Calculate Efficiency for Data
285
286 data_totaldata = np.absolute(data[:, 0])
287 data_splotWeights = data[:, 1]
288
289 data_entries = np.histogram(data_totaldata, rbins, weights=data_splotWeights)[0]
290 data_rvalueB0Average = (np.histogram(data_totaldata, rbins, weights=data_totaldata * data_splotWeights)[0] / data_entries)
291 data_rvalueMSB0Average = (
292 np.histogram(
293 data_totaldata,
294 rbins,
295 weights=(data_totaldata)**2 *
296 data_splotWeights)[0] /
297 data_entries)
298 data_rvalueStdB0Average = np.sqrt(abs(data_rvalueMSB0Average - data_rvalueB0Average * data_rvalueB0Average))
299
300 data_total_tagged = 0
301 for iEntries in data_entries:
302 data_total_tagged += iEntries
303
304 data_total_notTaggedWeighted = data_total_notTagged * data_total_tagged / data_totaldata.shape[0]
305
306 data_total_entries = data_total_tagged + data_total_notTaggedWeighted
307 data_event_fractionTotal = data_entries.astype(float) / data_total_entries
308
309 data_tagging_eff = data_total_tagged / data_total_entries
310 data_DeltaTagging_eff = math.sqrt(data_total_tagged * data_total_notTaggedWeighted**2 +
311 data_total_notTaggedWeighted * data_total_tagged**2) / (data_total_entries**2)
312 data_average_eff_eff = 0
313 data_uncertainty_eff_effAverage = 0
314 data_wvalue = array('f', [0] * r_size)
315 data_wvalueUncertainty = array('f', [0] * r_size)
316 data_iEffEfficiency = array('f', [0] * r_size)
317 data_iEffEfficiencyUncertainty = array('f', [0] * r_size)
318
319 print('****************** MEASURED EFFECTIVE EFFICIENCY FOR COMBINER USING ' + method + ' ***************************')
320 print('* *')
321 print('* --> DETERMINATION BASED ON DATA *')
322 print('* *')
323 print('* ------------------------------------------------------------------------------------------------ *')
324 print('* r-interval <r> Efficiency w *')
325 print('* ------------------------------------------------------------------------------------------------ *')
326
327 for i in range(0, r_size - 1):
328 data_wvalue[i] = (1 - data_rvalueB0Average[i]) / 2
329 data_wvalueUncertainty[i] = (data_rvalueStdB0Average[i] / math.sqrt(data_entries[i] - 1)) / 2
330 # fraction of events/all events
331
332 data_iEffEfficiency[i] = data_tagging_eff * (data_event_fractionTotal[i] *
333 data_rvalueB0Average[i] * data_rvalueB0Average[i])
334
335 data_iEffEfficiencyUncertainty[i] = data_rvalueB0Average[i] * \
336 math.sqrt((2 * data_total_entries * data_entries[i] *
337 (data_rvalueStdB0Average[i] / math.sqrt(data_entries[i] - 1)))**2 +
338 data_rvalueB0Average[i]**2 * data_entries[i] *
339 (data_total_entries * (data_total_entries - data_entries[i]) +
340 data_entries[i] * data_total_notTaggedWeighted)) / (data_total_entries**2)
341
342 # finally calculating the total effective efficiency
343 data_average_eff_eff = data_average_eff_eff + \
344 data_event_fractionTotal[i] * data_rvalueB0Average[i] * data_rvalueB0Average[i]
345 data_uncertainty_eff_effAverage = data_uncertainty_eff_effAverage + data_iEffEfficiencyUncertainty[i]**2
346
347 # intervallEff[i] = event_fractionTotal[i] * rvalueB0Average[i] * rvalueB0Average[i]
348 print('* ' + f'{r_subsample[i]:.3f}' + ' - ' + f'{r_subsample[i + 1]:.3f}' + ' ' +
349 f'{data_rvalueB0Average[i]:.3f}' + ' +- ' + f'{data_rvalueStdB0Average[i]:.4f}' + ' ' +
350 f'{data_event_fractionTotal[i]:.4f}' + ' ' +
351 f'{data_wvalue[i]:.4f}' + ' +- ' + f'{data_wvalueUncertainty[i]:.4f}' + ' ' + ' *')
352
353 data_uncertainty_eff_effAverage = math.sqrt(data_uncertainty_eff_effAverage)
354
355 print('* ------------------------------------------------------------------------------------------------ *')
356 print('* *')
357 print('* __________________________________________________________________________________________ *')
358 print('* | | *')
359 print('* | TOTAL NUMBER OF TAGGED EVENTS = ' +
360 f"{f'{data_total_tagged:.0f}':<24}" + f"{'| *':>36}")
361 print('* | | *')
362 print(
363 '* | TOTAL AVERAGE EFFICIENCY (q=+-1)= ' +
364 f'{data_tagging_eff * 100:.2f}' +
365 " +- " +
366 f'{data_DeltaTagging_eff * 100:.2f}' +
367 ' % | *')
368 print('* | | *')
369 print(
370 '* | TOTAL AVERAGE EFFECTIVE EFFICIENCY (q=+-1)= ' +
371 f'{data_average_eff_eff * 100:.6f}' +
372 " +- " +
373 f'{data_uncertainty_eff_effAverage * 100:.6f}' +
374 ' % | *')
375 print('* |__________________________________________________________________________________________| *')
376 print('* *')
377 print('****************************************************************************************************')
378
379 print(r'\begin{tabular}{ l r r r r}')
380 print(r'\hline')
381 print(r'$r$- Interval & $\varepsilon_i\ $ & $w_i \pm \delta w_i\ $ ' +
382 r' & $\varepsilon_{\text{eff}, i} \pm \delta\varepsilon_{\text{eff}, i}$ ' +
383 r'\\ \hline\hline')
384 for i in range(0, len(rbins) - 1):
385 print('$ ' + f'{r_subsample[i]:.3f}' + ' - ' + f'{r_subsample[i + 1]:.3f}' + '$ & $' +
386 f'{data_event_fractionTotal[i] * 100: 6.2f}' + '$ & $' +
387 f'{data_wvalue[i] * 100: 6.2f}' + r" \pm " + f'{data_wvalueUncertainty[i] * 100:2.2f}' + r'\, $ & $' +
388 f'{data_iEffEfficiency[i] * 100: 7.3f}' + # + '$ & $' +
389 r" \pm " + f'{data_iEffEfficiencyUncertainty[i] * 100:4.3f}' + r'\, $ \\ ')
390 print(r'\hline\hline')
391 print(
392 r'\multicolumn{1}{r}{Total} & \multicolumn{3}{r}{ $\varepsilon_\text{eff} = ' +
393 r'\sum_i \varepsilon_i \cdot \langle 1-2w_i\rangle^2 = ' +
394 f'{data_average_eff_eff * 100: 6.1f}' +
395 r" \pm " +
396 f'{data_uncertainty_eff_effAverage * 100: 6.1f}' +
397 r'\quad\,$ }' +
398 r' \\')
399 print(r'\hline')
400 print(r'\end{tabular}')
401
402 # Plot of Signal for B2TIP
403
404 # -----bins
405 bins = list(range(-25, 26, 1))
406 for i in range(0, len(bins)):
407 bins[i] = float(bins[i]) / 25
408 # ------
409
410 fig1 = plt.figure(1, figsize=(11, 8))
411 ax1 = plt.axes([0.21, 0.15, 0.76, 0.8])
412
413 qrDataPoints, qrDataBins = np.histogram(data[:, 0], bins=bins, weights=data_splotWeights)
414 qrDataPointsStd = np.sqrt(qrDataPoints)
415
416 qrDataPointsNorm, qrDataBins = np.histogram(data[:, 0], bins=bins, weights=data_splotWeights, density=1)
417 qrDataPointsStdNorm = qrDataPointsStd * qrDataPointsNorm / qrDataPoints
418
419 qrBincenters = 0.5 * (qrDataBins[1:] + qrDataBins[:-1])
420
421 p3 = ax1.errorbar(qrBincenters, qrDataPointsNorm, xerr=0.02, yerr=qrDataPointsStdNorm, elinewidth=2, mew=2, ecolor='k',
422 fmt='o', mfc='k', markersize=6, label=r'${\rm Data}$')
423
424 ax1.hist(mc[:, 0], bins=bins, density=1, histtype='step',
425 edgecolor='b', linewidth=4, linestyle='dashed', label=r'${\rm MC}$')
426
427 p2, = ax1.plot([], label=r'${\rm MC}$', linewidth=4, linestyle='dashed', c='b')
428
429 methodLabel = method
430 if method == "FANN":
431 methodLabel = "MLP"
432
433 ax1.set_ylabel(r'${\rm Fraction\ \ of\ \ Events\ /\ (\, 0.04\, )}$', fontsize=35)
434 ax1.set_xlabel(r'$(q\cdot r)_{\rm ' + methodLabel + '}$', fontsize=35)
435 plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1],
436 [r'$-1$', r'', r'$-0.5$', r'', r'$0$', r'', r'$0.5$', r'', r'$1$'], rotation=0, size=40)
437 ax1.tick_params(axis='y', labelsize=35)
438 ax1.legend([p3, p2], [r'${\rm Data}$', r'${\rm MC}$'], prop={'size': 35}, loc=0, numpoints=1)
439 ax1.grid(True)
440
441 if qrDataPointsNorm[np.argmax(qrDataPointsNorm)] < 1.4:
442 ax1.set_ylim(0, 1.4)
443
444 ax1.set_xlim(-1.005, 1.005)
445 plt.savefig(PATH + '/QR_B0_B0bar' + methodLabel + '.pdf')
446 fig1.clear()
447
448 # Plot of Continuum for B2TIP and Thesis
449
450 fig2 = plt.figure(2, figsize=(11, 8))
451 ax1 = plt.axes([0.21, 0.15, 0.76, 0.8])
452
453 qrDataPointsCont, qrDataBins = np.histogram(data[:, 0], bins=bins, weights=data[:, 2])
454 qrDataPointsContStd = np.sqrt(qrDataPointsCont)
455
456 qrDataPointsContNorm, qrDataBins = np.histogram(data[:, 0], bins=bins, weights=data[:, 2], density=1)
457 qrDataPointsContStdNorm = qrDataPointsContStd * qrDataPointsContNorm / qrDataPointsCont
458
459 qrBincenters = 0.5 * (qrDataBins[1:] + qrDataBins[:-1])
460
461 p3 = ax1.errorbar(qrBincenters, qrDataPointsContNorm, xerr=0.02, yerr=qrDataPointsContStdNorm, elinewidth=2, mew=2, ecolor='k',
462 fmt='o', mfc='k', mec='#990000', markersize=6, label=r'${\rm Data}$')
463
464 ax1.set_ylabel(r'${\rm Fraction\ \ of\ \ Events\ /\ (\, 0.04\, )}$', fontsize=35)
465 ax1.set_xlabel(r'$(q\cdot r)_{\rm ' + methodLabel + '}$', fontsize=35)
466 plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1],
467 [r'$-1$', r'', r'$-0.5$', r'', r'$0$', r'', r'$0.5$', r'', r'$1$'], rotation=0, size=40)
468 ax1.tick_params(axis='y', labelsize=35)
469 ax1.legend([p3], [r'${\rm Data}$'], prop={'size': 35}, loc=1, numpoints=1)
470 ax1.grid(True)
471
472 if qrDataPointsContNorm[np.argmax(qrDataPointsContNorm)] < 1.4:
473 ax1.set_ylim(0, 1.4)
474
475 ax1.set_xlim(-1.005, 1.005)
476 plt.savefig(PATH + '/QR_B0_B0bar' + methodLabel + 'Continuum.pdf')
477 fig2.clear()
478
479 data_total_tagged_Continuum = 0
480 for iEntries in qrDataPointsCont:
481 data_total_tagged_Continuum += iEntries
482
483 print('****************************************************************************************************')
484 print('* TOTAL NUMBER OF TAGGED CONTINUUM EVENTS = ' +
485 f"{f'{data_total_tagged_Continuum:.0f}':<24}" + f"{' *':>28}")
486
487 # Chi-squared between MC and Data
488
489 mcHistW = ROOT.TH1F("mcHistW", "mcHistW", int(r_size - 2), r_subsample)
490 dataHistW = ROOT.TH1F("dataHistW", "dataHistW", int(r_size - 2), r_subsample)
491
492 mcHistEff = ROOT.TH1F("mcHistEff", "mcHistEff", int(r_size - 2), r_subsample)
493 dataHistEff = ROOT.TH1F("dataHistEff", "dataHistEff", int(r_size - 2), r_subsample)
494
495 print('****************************************************************************************************')
496 print('* *')
497 print('* --> TEST OF COMPATIBILITY BETWEEN DATA AND MC *')
498
499 for i in range(0, len(rbins) - 1):
500
501 mcHistW.SetBinContent(i + 1, mc_wvalue[i])
502 mcHistW.SetBinError(i + 1, mc_wvalueUncertainty[i])
503
504 dataHistW.SetBinContent(i + 1, data_wvalue[i])
505 dataHistW.SetBinError(i + 1, data_wvalueUncertainty[i])
506
507 mcHistEff.SetBinContent(i + 1, mc_iEffEfficiency[i])
508 mcHistEff.SetBinError(i + 1, mc_iEffEfficiencyUncertainty[i])
509
510 dataHistEff.SetBinContent(i + 1, data_iEffEfficiency[i])
511 dataHistEff.SetBinError(i + 1, data_iEffEfficiencyUncertainty[i])
512
513 residualsW = array('d', [0] * r_size)
514 residualsEff = array('d', [0] * r_size)
515 print('* *')
516 print('* TEST ON w *')
517 print('* *')
518
519 chi2W = dataHistW.Chi2Test(mcHistW, "WW CHI2 P", residualsW)
520 print("1-CL for NDF = 7 is = ", ROOT.TMath.Prob(chi2W, 7))
521 print('* *')
522 print('* TEST ON EFFECTIVE EFFICIENCY *')
523 print('* *')
524 chi2Eff = dataHistEff.Chi2Test(mcHistEff, "WW P", residualsEff)
525 print("1-CL for NDF = 7 is = ", ROOT.TMath.Prob(chi2Eff, 7))
526
527 print('****************************************************************************************************')
528
529 # Plot of Signal with Residuals for Thesis
530
531 qrMCPoints, qrDataBins = np.histogram(mc[:, 0], bins=bins)
532 qrMCPointsStd = np.sqrt(qrMCPoints)
533 qrMCPointsNorm, qrDataBins = np.histogram(mc[:, 0], bins=bins, density=1)
534 qrMCPointsStdNorm = qrMCPointsStd * qrMCPointsNorm / qrMCPoints
535
536 qrDataMCResiduals = (qrDataPointsNorm - qrMCPointsNorm) \
537 / np.sqrt(qrDataPointsStdNorm * qrDataPointsStdNorm + qrMCPointsStdNorm * qrMCPointsStdNorm)
538
539 fig3 = plt.figure(3, figsize=(11, 11))
540 ax1 = plt.axes([0.21, 0.37, 0.76, 0.60])
541 # print('Here Log Scale')
542 # ax1.set_yscale('log', nonposy='clip') # Only for Semileptonic
543 p3 = ax1.errorbar(qrBincenters, qrDataPointsNorm, xerr=0.02, yerr=qrDataPointsStdNorm,
544 elinewidth=2, mew=2, ecolor='k',
545 fmt='o', mfc='k', mec='#006600', markersize=6, label=r'${\rm Data}$')
546
547 ax1.hist(mc[:, 0], bins=bins, density=1, histtype='step',
548 edgecolor='b', linewidth=4, linestyle='dashed', label=r'${\rm MC}$')
549
550 p2, = ax1.plot([], label=r'${\rm MC}$', linewidth=4, linestyle='dashed', c='b')
551
552 ax1.set_ylabel(r'${\rm Fraction\ \ of\ \ Events\ /\ (\, 0.04\, )}$', fontsize=35, labelpad=20)
553 # ax1.set_xlabel(r'', fontsize=35)
554 ax1.tick_params(axis='y', labelsize=35)
555 ax1.legend([p3, p2], [r'${\rm Data}$', r'${\rm MC}$'], prop={'size': 35}, loc=1, numpoints=1)
556 ax1.grid(True)
557
558 if qrDataPointsNorm[np.argmax(qrDataPointsNorm)] < 1.4:
559 ax1.set_ylim(0, 1.4)
560 plt.yticks([0, 0.25, 0.5, 0.75, 1.0, 1.25],
561 [r'$0.00$', r'$0.25$', r'$0.5$', r'$0.75$', r'$1.00$', r'$1.25$'], rotation=0, size=35)
562 ax1.set_xlim(-1.005, 1.005)
563 plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1],
564 [r'', r'', r'', r'', r'', r'', r'', r'', r''], rotation=0, size=35)
565
566 ax2 = plt.axes([0.21, 0.15, 0.76, 0.2])
567 ax2.errorbar(qrBincenters, qrDataMCResiduals, xerr=0.02, yerr=1, elinewidth=2, mew=2, ecolor='k',
568 fmt='o', mfc='k', mec='#006600', markersize=6, label=r'${\rm Data}$')
569 ax2.set_ylabel(r'${\rm Normalized}$' + '\n' + r'${\rm Residuals}$', fontsize=35, labelpad=20)
570 ax2.set_xlabel(r'$(q\cdot r)_{\rm ' + methodLabel + '}$', fontsize=35)
571 plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1],
572 [r'$-1$', r'', r'$-0.5$', r'', r'$0$', r'', r'$0.5$', r'', r'$1$'], rotation=0, size=35)
573 plt.yticks([-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5],
574 [r'', r'$-4$', r'', r'$-2$', r'', r'$0$', r'', r'$2$', r'', r'$4$', r''], rotation=0, size=25)
575 ax2.set_ylim(-5, 5)
576 # ax2.tick_params(axis='y', labelsize=30)
577 ax2.set_xlim(-1.005, 1.005)
578 ax2.xaxis.grid(True) # linestyle='--'
579 # plt.axhline(y= 4, xmin=-1.005, xmax=1.005, linewidth=1, color = 'k', linestyle = '-')
580 plt.axhline(y=3, xmin=-1.005, xmax=1.005, linewidth=2, color='#006600', linestyle='-')
581 # plt.axhline(y= 2, xmin=-1.005, xmax=1.005, linewidth=1, color = 'k', linestyle = '-')
582 plt.axhline(y=1, xmin=-1.005, xmax=1.005, linewidth=2, color='b', linestyle='-')
583 plt.axhline(y=-1, xmin=-1.005, xmax=1.005, linewidth=2, color='b', linestyle='-')
584 # plt.axhline(y=-2, xmin=-1.005, xmax=1.005, linewidth=1, color = 'k', linestyle = '-')
585 plt.axhline(y=-3, xmin=-1.005, xmax=1.005, linewidth=2, color='#006600', linestyle='-')
586 # plt.axhline(y=-4, xmin=-1.005, xmax=1.005, linewidth=1, color = 'k', linestyle = '-')
587 plt.savefig(PATH + '/QR_B0_B0bar' + methodLabel + 'WithRes.pdf')
588 fig3.clear()
589
590 return 1
591
592
593def plotWithResiduals(rooFitFrame, rooRealVar, dots, modelCurve, units, nameOfPlot, legend, removeArtifacts):
594
595 # plot for Thesis
596 rooFitFrame.Print()
597 rooFitFrame.GetXaxis().SetTitle("")
598 rooFitFrame.GetXaxis().SetLabelSize(0)
599
600 rooFitFrame.GetYaxis().SetTitleSize(0.072)
601 rooFitFrame.GetYaxis().SetTitleOffset(0.98)
602 rooFitFrame.GetYaxis().SetLabelSize(0.055)
603
604 pointsHist = ROOT.RooHist()
605
606 iBin = 0
607
608 # xValModel = ROOT.Double(-1.E30)
609 # yValModel = ROOT.Double(-1.E30)
610 xValDot = ROOT.Double(-1.E30)
611 yValDot = ROOT.Double(-1.E30)
612
613 iDotPoint = ROOT.RooRealVar("iDotPoint", "", 0.)
614 iModelPoint = ROOT.RooRealVar("iModelPoint", "", 0.)
615 iDotError = ROOT.RooRealVar("iDotError", "", 0.)
616 # yComb = ROOT.RooFormulaVar("yComb", "", "-1*@0", ROOT.RooArgList(y1RV))
617 iResidual = ROOT.RooFormulaVar("yComb", "", "(@0 - @1)/@2", ROOT.RooArgList(iDotPoint, iModelPoint, iDotError))
618
619 while iBin < dots.GetN() - 1:
620 dots.GetPoint(iBin, xValDot, yValDot)
621
622 iDotPoint.setVal(yValDot)
623 iModelPoint.setVal(modelCurve.interpolate(xValDot))
624 iDotError.setVal(float(dots.GetErrorYlow(iBin) + dots.GetErrorYhigh(iBin)))
625
626 if np.isnan(iResidual.getVal()) is not True:
627
628 residualValue = iResidual.getVal()
629
630 if removeArtifacts:
631 if abs(residualValue) > 3:
632 residualValue = 0
633
634 # print("xval = ", xValDot)
635 # print(iBin, " = ", iResidual.getVal(), ", @0 = ", yValDot, ", @1 = ",
636 # modelCurve.interpolate(xValDot), ", @2 = ",
637 # float(dots.GetErrorYlow(iBin) + dots.GetErrorYhigh(iBin)))
638 pointsHist.addBinWithXYError(xValDot, residualValue,
639 dots.GetErrorXlow(iBin),
640 dots.GetErrorXhigh(iBin),
641 1, 1,)
642
643 iBin += 1
644
645 pointsHist.SetMarkerStyle(dots.GetMarkerStyle())
646 pointsHist.SetMarkerSize(dots.GetMarkerSize())
647
648 rooFitFrameRes = rooRealVar.frame()
649 # rooFitFrameRes.addObject(mbcFrame.pullHist())
650 # rooFitFrameRes.addPlotable(pointsHist, "P")
651
652 rooFitFrameRes.SetTitle("")
653 rooFitFrameRes.GetXaxis().SetTitle(rooRealVar.GetTitle() + " " + units)
654 rooFitFrameRes.GetXaxis().SetTitleSize(0.2)
655 rooFitFrameRes.GetXaxis().SetTitleOffset(0.9)
656
657 rooFitFrameRes.GetXaxis().SetTickSize(0.07)
658 rooFitFrameRes.GetYaxis().SetTickSize(0.024)
659
660 rooFitFrameRes.GetYaxis().SetTitle("#splitline{Normalized}{ Residuals}")
661 rooFitFrameRes.GetYaxis().SetTitleSize(0.16)
662 rooFitFrameRes.GetYaxis().SetTitleOffset(0.4)
663
664 rooFitFrameRes.GetXaxis().SetLabelSize(0.125)
665 rooFitFrameRes.GetYaxis().SetLabelSize(0.120)
666
667 rooFitFrameRes.SetAxisRange(-5, 5, "Y")
668 rooFitFrameRes.GetYaxis().SetNdivisions(10)
669 rooFitFrameRes.GetYaxis().ChangeLabel(1, -1, 0.)
670 rooFitFrameRes.GetYaxis().ChangeLabel(3, -1, 0.)
671 rooFitFrameRes.GetYaxis().ChangeLabel(5, -1, 0.)
672 rooFitFrameRes.GetYaxis().ChangeLabel(7, -1, 0.)
673 rooFitFrameRes.GetYaxis().ChangeLabel(9, -1, 0.)
674 rooFitFrameRes.GetYaxis().ChangeLabel(11, -1, 0.)
675
676 gLine1 = ROOT.TLine(5.2, 3, 5.295, 3)
677 gLine2 = ROOT.TLine(5.2, 1, 5.295, 1)
678 gLine3 = ROOT.TLine(5.2, -1, 5.295, -1)
679 gLine4 = ROOT.TLine(5.2, -3, 5.295, -3)
680 gLine1.SetLineColor(ROOT.kGreen + 3)
681 gLine2.SetLineColor(ROOT.kBlue)
682 gLine3.SetLineColor(ROOT.kBlue)
683 gLine4.SetLineColor(ROOT.kGreen + 3)
684
685 gLine1.SetLineWidth(2)
686 gLine2.SetLineWidth(2)
687 gLine3.SetLineWidth(2)
688 gLine4.SetLineWidth(2)
689
690 # print("Number of Y axis bins = ", rooFitFrameRes.GetYaxis().GetNbins() )
691
692 c1 = ROOT.TCanvas("c1", "c1", 700, 640)
693 c1.SetBottomMargin(0)
694 c1.cd()
695
696 Pad1 = ROOT.TPad("p1", "p1", 0, 0.277, 1, 1, 0)
697 Pad2 = ROOT.TPad("p2", "p2", 0, 0, 1, 0.276, 0)
698 Pad1.Draw()
699 Pad2.Draw()
700
701 Pad1.SetLeftMargin(0.15)
702 Pad1.SetBottomMargin(0.02)
703 Pad1.SetTopMargin(0.06)
704
705 Pad2.SetLeftMargin(0.15)
706 Pad2.SetBottomMargin(0.4)
707 Pad2.SetTopMargin(0.01)
708
709 Pad2.cd()
710 rooFitFrameRes.Draw()
711 gLine1.Draw("SAME")
712 gLine2.Draw("SAME")
713 gLine3.Draw("SAME")
714 gLine4.Draw("SAME")
715 pointsHist.Draw("P SAME")
716 # Pad2.Update()
717
718 Pad1.cd()
719
720 rooFitFrame.Draw()
721
722 if legend != "":
723 legend.Draw()
724 # Pad1.Update()
725 nameOfPlot = nameOfPlot + "WithResiduals.pdf"
726 c1.SaveAs(nameOfPlot)
727 c1.Destructor()
728
729
730for treename in treenames:
731
732 print(treename)
733
734 tdat = ROOT.TChain(treename)
735 tMC = ROOT.TChain(treename)
736
737 for iFile in workingFilesData:
738 tdat.AddFile(iFile)
739
740 for iFile in workingFilesMC:
741 tMC.AddFile(iFile)
742
743 histo_notTaggedEventsMC = ROOT.TH1F('notTaggedEventsMC',
744 'Histogram for not tagged events',
745 1, -3.0, -1.0)
746
747 histo_notTaggedEventsData = ROOT.TH1F('notTaggedEventsData',
748 'Histogram for not tagged events',
749 1, -3.0, -1.0)
750
751 tMC.Draw(
752 'FBDT_qrCombined>>notTaggedEventsMC',
753 'abs(qrMC) == 0 && isSignal == 1 && FBDT_qrCombined < -1')
754
755 tdat.Draw(
756 'FBDT_qrCombined>>notTaggedEventsData',
757 'FBDT_qrCombined < -1')
758
759 mc_total_notTagged = histo_notTaggedEventsMC.GetEntries()
760 data_total_notTagged = histo_notTaggedEventsData.GetEntries()
761
762 B0_mbc = ROOT.RooRealVar("Mbc", "#it{m}_{bc}", 0, 5.2, 5.295, "GeV/#it{c}^{2}")
763 B0_deltaE = ROOT.RooRealVar("deltaE", "#Delta#it{E}", 0, -0.15, 0.15, "GeV/#it{c}^{2}")
764
765 B0_FBDT_qrCombined = ROOT.RooRealVar("FBDT_qrCombined", "FBDT_qrCombined", 0, -1.01, 1.01)
766 B0_FANN_qrCombined = ROOT.RooRealVar("FANN_qrCombined", "FANN_qrCombined", 0, -1.01, 1.01)
767
768 B0_qrMC = ROOT.RooRealVar("qrMC", "qrMC", 0., -100, 100)
769
770 argSet = ROOT.RooArgSet(B0_mbc, B0_deltaE)
771 argSet.add(B0_FBDT_qrCombined)
772 argSet.add(B0_FANN_qrCombined)
773 argSet.add(B0_qrMC)
774
775 cutData = "abs(FBDT_qrCombined) < 1.01"
776 cutMC = "abs(qrMC) ==1 && abs(FBDT_qrCombined) < 1.01"
777
778 data = ROOT.RooDataSet(
779 "data",
780 "data",
781 tdat,
782 argSet,
783 cutData)
784
785 mc = ROOT.RooDataSet(
786 "mc",
787 "mc",
788 tMC,
789 argSet,
790 cutMC)
791
792 fitMC = ROOT.RooDataSet("fitMC", "fitMC", ROOT.RooArgSet(B0_mbc, B0_FBDT_qrCombined, B0_FANN_qrCombined, B0_qrMC), "")
793 fitData = ROOT.RooDataSet("fitData", "fitData", ROOT.RooArgSet(B0_mbc, B0_FBDT_qrCombined, B0_FANN_qrCombined), "")
794
795 fitMC.append(mc)
796 fitData.append(data)
797
798 data_total_Tagged = fitData.numEntries()
799 data_total_notTagged = histo_notTaggedEventsData.GetEntries()
800
801 maxBo_mbc = tdat.GetMaximum("Mbc") # maximum B0_mbc = 5.289213180541992
802 print("maximum B0_mbc = ", maxBo_mbc)
803
804 B0_mbc.setRange("fitRange", 5.2, maxBo_mbc)
805
806# Definition of Signal Components and Fit on MC
807
808 mbcMuGauss1 = ROOT.RooRealVar("mbcMuGauss1", "mbcMuGauss1", 5.27943e+00, 5.27, 5.2893, "GeV")
809 mbcSigmaGauss1 = ROOT.RooRealVar("mbcSigmaGauss1", "mbcSigmaGauss1", 2.62331e-03, 0, 0.1, "GeV")
810 mbcGauss1 = ROOT.RooGaussModel("mbcGauss1", "mbcGauss1", B0_mbc, mbcMuGauss1, mbcSigmaGauss1)
811
812 mbcMu2Shift = ROOT.RooRealVar("mbcMu2Shift", "mbcMu2Shift", 7.62270e-03, -0.01, 0.01)
813 mbcSig2Factor = ROOT.RooRealVar("mbcSig2Factor", "mbcSig2Factor", 2.83547e-01, 0, 1)
814 mbcMuGauss2 = ROOT.RooFormulaVar("mbcMuGauss2", "mbcMuGauss2", "@0-@1", ROOT.RooArgList(mbcMuGauss1, mbcMu2Shift))
815 mbcSigmaGauss2 = ROOT.RooFormulaVar("mbcSigmaGauss2", "mbcSigmaGauss2", "@0/@1", ROOT.RooArgList(mbcSigmaGauss1, mbcSig2Factor))
816 mbcGauss2 = ROOT.RooGaussModel("mbcGauss2", "mbcGauss2", B0_mbc, mbcMuGauss2, mbcSigmaGauss2)
817
818 mbcFracGauss2 = ROOT.RooRealVar("mbcFracGauss2", "mbcFracGauss2", 9.93351e-01, 0.0, 1.)
819
820 mbcSignalModel = ROOT.RooAddModel(
821 "mbcSignalModel", "mbcSignalModel", ROOT.RooArgList(
822 mbcGauss1, mbcGauss2), ROOT.RooArgList(mbcFracGauss2))
823
824 mcFit = mbcSignalModel.fitTo(
825 fitMC,
826 ROOT.RooFit.Minos(ROOT.kFALSE), ROOT.RooFit.Extended(ROOT.kFALSE),
827 ROOT.RooFit.NumCPU(8), ROOT.RooFit.Save())
828
829 mbcMu2Shift.setConstant()
830 mbcSig2Factor.setConstant()
831 mbcFracGauss2.setConstant()
832
833 mbcSignalArgset1 = ROOT.RooArgSet(mbcGauss1)
834 mbcSignalArgset2 = ROOT.RooArgSet(mbcGauss2)
835
836 mbcLogFrameMC = B0_mbc.frame()
837 fitMC.plotOn(mbcLogFrameMC, ROOT.RooFit.LineWidth(4))
838 mbcSignalModel.plotOn(mbcLogFrameMC)
839 mbcSignalModel.plotOn(
840 mbcLogFrameMC,
841 ROOT.RooFit.Components(mbcSignalArgset1),
842 ROOT.RooFit.LineColor(ROOT.kRed + 2), ROOT.RooFit.LineWidth(4))
843 mbcSignalModel.plotOn(
844 mbcLogFrameMC,
845 ROOT.RooFit.Components(mbcSignalArgset2),
846 ROOT.RooFit.LineColor(ROOT.kGreen + 3), ROOT.RooFit.LineWidth(4))
847
848 mbcLogFrameMC.SetTitle("")
849 mbcXtitle = B0_mbc.GetTitle() + " [GeV/#it{c}^{2}]"
850 mbcLogFrameMC.GetXaxis().SetTitle(mbcXtitle)
851 mbcLogFrameMC.GetXaxis().SetTitleSize(0.05)
852 mbcLogFrameMC.GetXaxis().SetLabelSize(0.045)
853 mbcLogFrameMC.GetYaxis().SetTitleSize(0.05)
854 mbcLogFrameMC.GetYaxis().SetTitleOffset(1.5)
855 mbcLogFrameMC.GetYaxis().SetLabelSize(0.045)
856
857 c1 = ROOT.TCanvas("c1", "c1", 1400, 1100)
858 c1.cd()
859 Pad = ROOT.TPad("p1", "p1", 0, 0, 1, 1, 0, 0, 0)
860 Pad.SetLeftMargin(0.15)
861 Pad.SetBottomMargin(0.15)
862 Pad.SetLogy()
863 Pad.Draw()
864 Pad.cd()
865 mbcLogFrameMC.Draw()
866 Pad.Update()
867 nPlot = PATH + "/B0_mbcMC.pdf"
868 c1.SaveAs(nPlot)
869 c1.Destructor()
870
871# Definition of Continuum component and fit of Signal Model and Continuum component on Data
872
873 mbcSignYield = ROOT.RooRealVar("mbcSignYield", "mbcSignYield", 1.53208e+03, 0.0, 20000)
874 mbcContYield = ROOT.RooRealVar("mbcContYield", "mbcContYield", 1.13456e+03, 0.0, 20000)
875
876 mbcArgusMax = ROOT.RooRealVar("mbcArgusMax", "mbcArgusMax", maxBo_mbc, "GeV")
877 mbcArgusC = ROOT.RooRealVar("mbcArgusC", "mbcArgusC", -4.65996e+01, -100, 0)
878 mbcArgusPow = ROOT.RooRealVar("mbcArgusPow", "mbcArgusPow", 5.13442e-01, -2, 2, "GeV")
879 mbcArgus = ROOT.RooArgusBG("mbcArgus", "mbcArgus", B0_mbc, mbcArgusMax, mbcArgusC, mbcArgusPow)
880
881 mbcModel = ROOT.RooAddPdf(
882 "mbcModel", "mbcModel", ROOT.RooArgList(
883 mbcSignalModel, mbcArgus), ROOT.RooArgList(
884 mbcSignYield, mbcContYield))
885
886 dataFit = mbcModel.fitTo(
887 fitData,
888 ROOT.RooFit.Minos(ROOT.kFALSE), ROOT.RooFit.Extended(ROOT.kFALSE),
889 ROOT.RooFit.NumCPU(8), ROOT.RooFit.Save())
890
891 mbcMuGauss1.setConstant()
892 mbcSigmaGauss1.setConstant()
893 mbcArgusC.setConstant()
894 mbcArgusPow.setConstant()
895
896# Now we use the SPlot class to add SWeights to our data set
897
898 sData = ROOT.RooStats.SPlot("sData", "SPlot using B0_mbc", fitData, mbcModel, ROOT.RooArgList(mbcSignYield, mbcContYield))
899
900
901# Check that our weights have the desired properties
902
903 print("Check SWeights:")
904
905 print("Yield of mbcSignYield is ", mbcSignYield.getVal(), ". From sWeights it is ",
906 sData.GetYieldFromSWeight("mbcSignYield"))
907
908 print("Yield of mbcContYield is ", mbcContYield.getVal(), ". From sWeights it is ",
909 sData.GetYieldFromSWeight("mbcContYield"))
910
911 # for i in range(0, 10):
912 # print("mbcSignYield Weight ", sData.GetSWeight(i,"mbcSignYield"),
913 # " mbcContYield Weight ", sData.GetSWeight(i,"mbcContYield"),
914 # " Total Weight ", sData.GetSumOfEventSWeight(i))
915
916 dataFitExtended = mbcModel.fitTo(
917 fitData, ROOT.RooFit.Extended(),
918 ROOT.RooFit.Minos(ROOT.kFALSE), ROOT.RooFit.Extended(ROOT.kFALSE),
919 ROOT.RooFit.NumCPU(8), ROOT.RooFit.Save())
920
921 mbcModel.Print()
922
923 mbcArgset1 = ROOT.RooArgSet(mbcArgus)
924 mbcArgset2 = ROOT.RooArgSet(mbcSignalModel)
925
926 B0_mbc.setBins(38)
927 mbcFrame = B0_mbc.frame()
928 fitData.plotOn(mbcFrame)
929 mbcModel.plotOn(mbcFrame, ROOT.RooFit.LineWidth(4))
930 mbcModel.plotOn(
931 mbcFrame,
932 ROOT.RooFit.Components(mbcArgset1),
933 ROOT.RooFit.LineColor(
934 ROOT.kRed + 2),
935 ROOT.RooFit.LineStyle(5),
936 ROOT.RooFit.LineWidth(4))
937 mbcModel.plotOn(
938 mbcFrame,
939 ROOT.RooFit.Components(mbcArgset2),
940 ROOT.RooFit.LineColor(
941 ROOT.kGreen + 3),
942 ROOT.RooFit.LineStyle(7),
943 ROOT.RooFit.LineWidth(4))
944
945 mbcFrame.SetTitle("")
946 mbcFrame.GetXaxis().SetTitle(mbcXtitle)
947 mbcFrame.SetTitleOffset(0.9, "X")
948 mbcFrame.GetXaxis().SetTitleSize(0.075)
949 mbcFrame.GetXaxis().SetLabelSize(0.045)
950 mbcFrame.GetYaxis().SetTitleSize(0.06)
951 mbcFrame.GetYaxis().SetTitleOffset(0.95)
952 mbcFrame.GetYaxis().SetLabelSize(0.045)
953
954 c1 = ROOT.TCanvas("c1", "c1", 10, 15, 700, 500)
955 c1.SetLeftMargin(0.15)
956 c1.SetRightMargin(0.04)
957 c1.SetTopMargin(0.03)
958 c1.SetBottomMargin(0.15)
959 c1.cd()
960 # Pad = ROOT.TPad("p1", "p1", 0, 0, 1, 1, 0, 0, 0)
961 # Pad.SetLeftMargin(0.15)
962 # Pad.SetBottomMargin(0.15)
963 # Pad.Draw()
964 # Pad.cd()
965 mbcFrame.Draw()
966 # Pad.Update()
967 nPlot = PATH + "/B0_mbc.pdf"
968 c1.SaveAs(nPlot)
969 c1.Destructor()
970
971 # Plot Mbc MC with residuals
972 mbcFrameMC = B0_mbc.frame()
973 mbcFrameMC.SetTitle("")
974
975 fitMC.plotOn(mbcFrameMC)
976 mbcSignalModel.plotOn(mbcFrameMC, ROOT.RooFit.LineColor(ROOT.kGreen + 3), ROOT.RooFit.LineStyle(7), ROOT.RooFit.LineWidth(4))
977
978 dotsMC = mbcFrameMC.getHist("h_fitMC")
979 modelCurveMC = mbcFrameMC.getCurve("mbcSignalModel_Norm[Mbc]")
980
981 dotsMC.SetFillColor(ROOT.kWhite)
982 modelCurveMC.SetFillColor(ROOT.kWhite)
983
984 mbcFrameMC.Print()
985
986 lMC = ROOT.TLegend(0.25, 0.63, 0.5, 0.89)
987 lMC.AddEntry(dotsMC, "Belle MC")
988 lMC.AddEntry(modelCurveMC, 'Signal')
989
990 lMC.SetTextSize(0.055)
991
992 plotWithResiduals(mbcFrameMC, B0_mbc, dotsMC, modelCurveMC, "[GeV/#it{c}^{2}]", "mbcMC", lMC, True)
993
994 # Plot of Model with Mbc and residuals
995
996 dots = mbcFrame.getHist("h_fitData")
997 modelCurve = mbcFrame.getCurve("mbcModel_Norm[Mbc]")
998 signalCurve = mbcFrame.getCurve("mbcModel_Norm[Mbc]_Comp[mbcSignalModel]")
999 continuumCurve = mbcFrame.getCurve("mbcModel_Norm[Mbc]_Comp[mbcArgus]")
1000
1001 dots.SetFillColor(ROOT.kWhite)
1002 modelCurve.SetFillColor(ROOT.kWhite)
1003 signalCurve.SetFillColor(ROOT.kWhite)
1004 continuumCurve.SetFillColor(ROOT.kWhite)
1005
1006 legend = ROOT.TLegend(0.25, 0.43, 0.5, 0.89)
1007 legend.AddEntry(dots, "Belle data")
1008 legend.AddEntry(modelCurve, 'Fit result')
1009 legend.AddEntry(signalCurve, "Signal")
1010 legend.AddEntry(continuumCurve, 'Continuum')
1011
1012 legend.SetTextSize(0.055)
1013
1014 plotWithResiduals(mbcFrame, B0_mbc, dots, modelCurve, "[GeV/#it{c}^{2}]", "Mbc", legend, False)
1015
1016 signalData = ROOT.RooDataSet("signalData", "signalData", fitData, fitData.get(), "", "mbcSignYield_sw")
1017
1018 print("Variable Names")
1019
1020 qrSigFrame = B0_FBDT_qrCombined.frame()
1021 signalData.plotOn(qrSigFrame, ROOT.RooFit.DataError(ROOT.RooAbsData.SumW2))
1022
1023 qrSigFrame.SetTitle("")
1024 mbcXtitle = "qr"
1025 qrSigFrame.GetXaxis().SetTitle(mbcXtitle)
1026 qrSigFrame.GetXaxis().SetTitleSize(0.05)
1027 qrSigFrame.GetXaxis().SetLabelSize(0.045)
1028 qrSigFrame.GetYaxis().SetTitleSize(0.05)
1029 qrSigFrame.GetYaxis().SetTitleOffset(1.5)
1030 qrSigFrame.GetYaxis().SetLabelSize(0.045)
1031
1032 c1 = ROOT.TCanvas("c1", "c1", 1400, 1100)
1033 c1.cd()
1034 Pad = ROOT.TPad("p1", "p1", 0, 0, 1, 1, 0, 0, 0)
1035 Pad.SetLeftMargin(0.15)
1036 Pad.SetBottomMargin(0.15)
1037 Pad.Draw()
1038 Pad.cd()
1039 qrSigFrame.Draw()
1040
1041 Pad.Update()
1042 nPlot = PATH + "/B0_qrSignal.pdf"
1043 c1.SaveAs(nPlot)
1044 c1.Destructor()
1045
1046 continuumData = ROOT.RooDataSet("continuumData", "continuumData", fitData, fitData.get(), "", "mbcContYield_sw")
1047
1048 qrContFrame = B0_FBDT_qrCombined.frame()
1049 continuumData.plotOn(qrContFrame, ROOT.RooFit.DataError(ROOT.RooAbsData.SumW2))
1050
1051 qrContFrame.SetTitle("")
1052 mbcXtitle = "qr"
1053 qrContFrame.GetXaxis().SetTitle(mbcXtitle)
1054 qrContFrame.GetXaxis().SetTitleSize(0.05)
1055 qrContFrame.GetXaxis().SetLabelSize(0.045)
1056 qrContFrame.GetYaxis().SetTitleSize(0.05)
1057 qrContFrame.GetYaxis().SetTitleOffset(1.5)
1058 qrContFrame.GetYaxis().SetLabelSize(0.045)
1059
1060 c1 = ROOT.TCanvas("c1", "c1", 1400, 1100)
1061 c1.cd()
1062 Pad = ROOT.TPad("p1", "p1", 0, 0, 1, 1, 0, 0, 0)
1063 Pad.SetLeftMargin(0.15)
1064 Pad.SetBottomMargin(0.15)
1065 Pad.Draw()
1066 Pad.cd()
1067 qrContFrame.Draw()
1068 Pad.Update()
1069 nPlot = PATH + "/B0_qrContinuum.pdf"
1070 c1.SaveAs(nPlot)
1071 c1.Destructor()
1072
1073 fitData.Print()
1074
1075 fitMC.Print()
1076
1077 sPlotResultsFBDT = np.zeros((fitData.numEntries(), 3))
1078 sPlotResultsFANN = np.zeros((fitData.numEntries(), 3))
1079 for i in range(fitData.numEntries()):
1080 row = fitData.get(i)
1081 sPlotResultsFBDT[i][0] = row.getRealValue('FBDT_qrCombined', 0, ROOT.kTRUE)
1082 sPlotResultsFBDT[i][1] = row.getRealValue('mbcSignYield_sw', 0, ROOT.kTRUE)
1083 sPlotResultsFBDT[i][2] = row.getRealValue('mbcContYield_sw', 0, ROOT.kTRUE)
1084 sPlotResultsFANN[i][0] = row.getRealValue('FANN_qrCombined', 0, ROOT.kTRUE)
1085 sPlotResultsFANN[i][1] = row.getRealValue('mbcSignYield_sw', 0, ROOT.kTRUE)
1086 sPlotResultsFANN[i][2] = row.getRealValue('mbcContYield_sw', 0, ROOT.kTRUE)
1087
1088 mcSPlotFBDT = np.zeros((fitMC.numEntries(), 2))
1089 mcSPlotFANN = np.zeros((fitMC.numEntries(), 2))
1090 for i in range(fitMC.numEntries()):
1091 row = fitMC.get(i)
1092 mcSPlotFBDT[i][0] = row.getRealValue('FBDT_qrCombined', 0, ROOT.kTRUE)
1093 mcSPlotFBDT[i][1] = row.getRealValue('qrMC', 0, ROOT.kTRUE)
1094 mcSPlotFANN[i][0] = row.getRealValue('FANN_qrCombined', 0, ROOT.kTRUE)
1095 mcSPlotFANN[i][1] = row.getRealValue('qrMC', 0, ROOT.kTRUE)
1096
1097 efficiencyCalculator(mcSPlotFBDT, sPlotResultsFBDT, "FBDT")
1098 efficiencyCalculator(mcSPlotFANN, sPlotResultsFANN, "FANN")