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