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