Belle II Software  release-05-01-25
B2JpsiKs_mu_qrBelleDataSplot.py
1 #!/usr/bin/env python
2 # -*- coding: utf-8 -*-
3 
4 # * ********************************** Flavor Tagging ************************************
5 # *
6 # * This validation script performs an Splot analysis of Belle Data using mbc as *
7 # * Splot variable in order to extract the qr distribution of signal B0s. *
8 # * The signal component of Mbc is first fitted on converted Belle MC using a *
9 # * 2 Gaussian function. The Mbc component is used together with an Argus function *
10 # * to fit the mbc distribution of converted Belle data. *
11 # * An Splot is performed and its output variable used to weigth the qr *
12 # * distribution. *
13 # * For more information see Sec. 4.8 in BELLE2-PTHESIS-2018-003 *
14 # * Usage: *
15 # * basf2 B2JpsiKs_mu_qrBelleDataSplot.py RootNtupleFilesMC RootNtupleFilesData treeName *
16 # *
17 # * Contributors: F. Abudinen (December 2018) *
18 # *
19 # ******************************************************************************************
20 
21 import ROOT
22 import numpy as np
23 import matplotlib as mpl
24 mpl.use('Agg')
25 mpl.rcParams['text.usetex'] = True
26 from defaultEvaluationParameters import categories, rbins, r_subsample, r_size
27 import matplotlib.pyplot as plt
28 from array import array
29 import pylab
30 import sys
31 import glob
32 import math
33 import random
34 
35 
36 PATH = "."
37 
38 workingFilesMC = str(sys.argv[1])
39 workingFilesData = str(sys.argv[2])
40 
41 workingFilesMC = glob.glob(str(workingFilesMC))
42 workingFilesData = glob.glob(str(workingFilesData))
43 
44 lim = 10
45 limZ = 0.03
46 limZsig = 0.05
47 
48 treenames = [str(sys.argv[3])]
49 
50 mc_total_notTagged = 0
51 
52 data_total_notTagged = 0
53 
54 data_total_Tagged = 0
55 
56 ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.WARNING)
57 
58 
59 def efficiencyCalculator(mc, data, method):
60 
61  # Calculate Efficiency for MC
62 
63  mc_dataB0 = mc[np.where(mc[:, 1] > 0)]
64  mc_dataB0bar = mc[np.where(mc[:, 1] < 0)]
65 
66  mc_totaldataB0 = np.absolute(mc_dataB0[:, 0])
67  mc_totaldataB0bar = np.absolute(mc_dataB0bar[:, 0])
68 
69  mc_rvalueB0 = (np.histogram(mc_totaldataB0, rbins, weights=mc_totaldataB0)[0] / np.histogram(mc_totaldataB0, rbins)[0])
70  mc_rvalueMSB0 = (
71  np.histogram(
72  mc_totaldataB0,
73  rbins,
74  weights=mc_totaldataB0 *
75  mc_totaldataB0)[0] /
76  np.histogram(
77  mc_totaldataB0,
78  rbins)[0])
79  mc_rvalueStdB0 = np.sqrt(abs(mc_rvalueMSB0 - mc_rvalueB0 * mc_rvalueB0))
80  mc_entriesB0 = np.histogram(mc_totaldataB0, rbins)[0]
81 
82  mc_rvalueB0bar = (
83  np.histogram(
84  mc_totaldataB0bar,
85  rbins,
86  weights=mc_totaldataB0bar)[0] /
87  np.histogram(
88  mc_totaldataB0bar,
89  rbins)[0])
90  mc_rvalueMSB0bar = (
91  np.histogram(
92  mc_totaldataB0bar,
93  rbins,
94  weights=mc_totaldataB0bar *
95  mc_totaldataB0bar)[0] /
96  np.histogram(
97  mc_totaldataB0bar,
98  rbins)[0])
99  mc_rvalueStdB0bar = np.sqrt(abs(mc_rvalueMSB0bar - mc_rvalueB0bar * mc_rvalueB0bar))
100  mc_entriesB0bar = np.histogram(mc_totaldataB0bar, rbins)[0]
101 
102  mc_total_entriesB0 = mc_totaldataB0.shape[0] + mc_total_notTagged / 2
103  mc_total_tagged_B0 = mc_totaldataB0.shape[0]
104  mc_event_fractionB0 = mc_entriesB0.astype(float) / mc_total_entriesB0
105 
106  mc_total_entriesB0bar = mc_totaldataB0bar.shape[0] + mc_total_notTagged / 2
107  mc_total_tagged_B0bar = mc_totaldataB0bar.shape[0]
108  mc_event_fractionB0bar = mc_entriesB0bar.astype(float) / mc_total_entriesB0bar
109 
110  mc_total_entries = mc_total_entriesB0 + mc_total_entriesB0bar
111  mc_total_tagged = mc_totaldataB0.shape[0] + mc_totaldataB0bar.shape[0]
112  mc_event_fractionTotal = (mc_entriesB0.astype(float) + mc_entriesB0bar.astype(float)) / mc_total_entries
113 
114  mc_tagging_eff = mc_total_tagged / (mc_total_tagged + mc_total_notTagged)
115  mc_DeltaTagging_eff = math.sqrt(mc_total_tagged * mc_total_notTagged**2 +
116  mc_total_notTagged * mc_total_tagged**2) / (mc_total_entries**2)
117  mc_tot_eff_effB0 = 0
118  mc_tot_eff_effB0bar = 0
119  mc_uncertainty_eff_effB0 = 0
120  mc_uncertainty_eff_effB0bar = 0
121  mc_uncertainty_eff_effAverage = 0
122  mc_diff_eff_Uncertainty = 0
123  mc_event_fractionDiff = array('f', [0] * r_size)
124  mc_rvalueB0Average = array('f', [0] * r_size)
125  mc_rvalueStdB0Average = array('f', [0] * r_size)
126  mc_wvalue = array('f', [0] * r_size)
127  mc_wvalueUncertainty = array('f', [0] * r_size)
128  mc_wvalueB0 = array('f', [0] * r_size)
129  mc_wvalueB0bar = array('f', [0] * r_size)
130  mc_wvalueDiff = array('f', [0] * r_size)
131  mc_wvalueDiffUncertainty = array('f', [0] * r_size)
132  mc_entries = array('f', [0] * r_size)
133  mc_iEffEfficiency = array('f', [0] * r_size)
134  mc_iEffEfficiencyUncertainty = array('f', [0] * r_size)
135  mc_iEffEfficiencyB0Uncertainty = array('f', [0] * r_size)
136  mc_iEffEfficiencyB0barUncertainty = array('f', [0] * r_size)
137  mc_iDeltaEffEfficiency = array('f', [0] * r_size)
138  mc_iDeltaEffEfficiencyUncertainty = array('f', [0] * r_size)
139 
140  print('****************** MEASURED EFFECTIVE EFFICIENCY FOR COMBINER USING ' + method + ' ***************************')
141  print('* *')
142  print('* --> DETERMINATION BASED ON MONTE CARLO *')
143  print('* *')
144  print('* ------------------------------------------------------------------------------------------------ *')
145  print('* r-interval <r> Efficiency Delta_Effcy w Delta_w *')
146  print('* ------------------------------------------------------------------------------------------------ *')
147 
148  for i in range(0, r_size - 1):
149  mc_rvalueB0Average[i] = (mc_rvalueB0[i] + mc_rvalueB0bar[i]) / 2
150  mc_rvalueStdB0Average[i] = math.sqrt(mc_rvalueStdB0[i]**2 / (mc_entriesB0[0] - 1) +
151  mc_rvalueStdB0bar[i]**2 / (mc_entriesB0bar[0] - 1)) / 2
152  mc_wvalueB0[i] = (1 - mc_rvalueB0[i]) / 2
153  mc_wvalueB0bar[i] = (1 - mc_rvalueB0bar[i]) / 2
154  mc_wvalueDiff[i] = mc_wvalueB0[i] - mc_wvalueB0bar[i]
155  mc_wvalueDiffUncertainty[i] = math.sqrt((mc_rvalueStdB0[i] / 2)**2 / (mc_entriesB0[0] - 1) +
156  (mc_rvalueStdB0bar[i] / 2)**2 / (mc_entriesB0bar[0] - 1))
157  mc_wvalue[i] = (mc_wvalueB0[i] + mc_wvalueB0bar[i]) / 2
158  mc_wvalueUncertainty[i] = mc_wvalueDiffUncertainty[i] / 2
159  # fraction of events/all events
160 
161  mc_event_fractionDiff[i] = (mc_entriesB0[i] - mc_entriesB0bar[i]) / mc_total_entries
162 
163  mc_iEffEfficiency[i] = mc_tagging_eff * (mc_event_fractionB0[i] * mc_rvalueB0[i] * mc_rvalueB0[i] +
164  mc_event_fractionB0bar[i] * mc_rvalueB0bar[i] * mc_rvalueB0bar[i]) / 2
165 
166  mc_iEffEfficiencyB0Uncertainty[i] = mc_rvalueB0[i] * \
167  math.sqrt((2 * mc_total_entriesB0 * mc_entriesB0[i] * mc_rvalueStdB0[i])**2 / (mc_entriesB0[0] - 1) +
168  mc_rvalueB0[i]**2 * mc_entriesB0[i] *
169  (mc_total_entriesB0 * (mc_total_entriesB0 - mc_entriesB0[i]) +
170  mc_entriesB0[i] * mc_total_notTagged)) / (mc_total_entriesB0**2)
171  mc_iEffEfficiencyB0barUncertainty[i] = mc_rvalueB0bar[i] * \
172  math.sqrt((2 * mc_total_entriesB0bar * mc_entriesB0bar[i] * mc_rvalueStdB0bar[i])**2 / (mc_entriesB0bar[0] - 1) +
173  mc_rvalueB0bar[i]**2 * mc_entriesB0bar[i] *
174  (mc_total_entriesB0bar * (mc_total_entriesB0bar - mc_entriesB0bar[i]) +
175  mc_entriesB0bar[i] * mc_total_notTagged)) / (mc_total_entriesB0bar**2)
176 
177  mc_iDeltaEffEfficiency[i] = mc_event_fractionB0[i] * mc_rvalueB0[i] * \
178  mc_rvalueB0[i] - mc_event_fractionB0bar[i] * mc_rvalueB0bar[i] * mc_rvalueB0bar[i]
179 
180  mc_iDeltaEffEfficiencyUncertainty[i] = math.sqrt(
181  mc_iEffEfficiencyB0Uncertainty[i]**2 +
182  mc_iEffEfficiencyB0barUncertainty[i]**2)
183 
184  mc_iEffEfficiencyUncertainty[i] = mc_iDeltaEffEfficiencyUncertainty[i] / 2
185 
186  mc_diff_eff_Uncertainty = mc_diff_eff_Uncertainty + mc_iDeltaEffEfficiencyUncertainty[i]**2
187 
188  # finally calculating the total effective efficiency
189  mc_tot_eff_effB0 = mc_tot_eff_effB0 + mc_event_fractionB0[i] * mc_rvalueB0[i] * mc_rvalueB0[i]
190  mc_tot_eff_effB0bar = mc_tot_eff_effB0bar + mc_event_fractionB0bar[i] * mc_rvalueB0bar[i] * mc_rvalueB0bar[i]
191  mc_uncertainty_eff_effAverage = mc_uncertainty_eff_effAverage + mc_iEffEfficiencyUncertainty[i]**2
192  mc_uncertainty_eff_effB0 = mc_uncertainty_eff_effB0 + mc_iEffEfficiencyB0Uncertainty[i]**2
193  mc_uncertainty_eff_effB0bar = mc_uncertainty_eff_effB0bar + mc_iEffEfficiencyB0barUncertainty[i]**2
194 
195  # intervallEff[i] = event_fractionTotal[i] * rvalueB0Average[i] * rvalueB0Average[i]
196  print('* ' + '{:.3f}'.format(r_subsample[i]) + ' - ' + '{:.3f}'.format(r_subsample[i + 1]) + ' ' +
197  '{:.3f}'.format(mc_rvalueB0Average[i]) + ' +- ' + '{:.4f}'.format(mc_rvalueStdB0Average[i]) + ' ' +
198  '{:.4f}'.format(mc_event_fractionTotal[i]) + ' ' +
199  '{: .4f}'.format(mc_event_fractionDiff[i]) + ' ' +
200  '{:.4f}'.format(mc_wvalue[i]) + ' +- ' + '{:.4f}'.format(mc_wvalueUncertainty[i]) + ' ' +
201  '{: .4f}'.format(mc_wvalueDiff[i]) + ' +- ' + '{:.4f}'.format(mc_wvalueDiffUncertainty[i]) + ' *')
202 
203  mc_average_eff_eff = (mc_tot_eff_effB0 + mc_tot_eff_effB0bar) / 2
204  mc_uncertainty_eff_effAverage = math.sqrt(mc_uncertainty_eff_effAverage)
205  mc_uncertainty_eff_effB0 = math.sqrt(mc_uncertainty_eff_effB0)
206  mc_uncertainty_eff_effB0bar = math.sqrt(mc_uncertainty_eff_effB0bar)
207  mc_diff_eff = mc_tot_eff_effB0 - mc_tot_eff_effB0bar
208  mc_diff_eff_Uncertainty = math.sqrt(mc_diff_eff_Uncertainty)
209 
210  print('* ------------------------------------------------------------------------------------------------ *')
211  print('* *')
212  print('* __________________________________________________________________________________________ *')
213  print('* | | *')
214  print('* | TOTAL NUMBER OF TAGGED EVENTS = ' +
215  '{:<24}'.format("%.0f" % mc_total_tagged) + '{:>36}'.format('| *'))
216  print('* | | *')
217  print(
218  '* | TOTAL AVERAGE EFFICIENCY (q=+-1)= ' +
219  '{:.2f}'.format(
220  mc_tagging_eff *
221  100) +
222  " +- " +
223  '{:.2f}'.format(
224  mc_DeltaTagging_eff *
225  100) +
226  ' % | *')
227  print('* | | *')
228  print(
229  '* | TOTAL AVERAGE EFFECTIVE EFFICIENCY (q=+-1)= ' +
230  '{:.6f}'.format(
231  mc_average_eff_eff *
232  100) +
233  " +- " +
234  '{:.6f}'.format(
235  mc_uncertainty_eff_effAverage *
236  100) +
237  ' % | *')
238  print('* | | *')
239  print(
240  '* | TOTAL AVERAGE EFFECTIVE EFFICIENCY ASYMMETRY (q=+-1)= ' +
241  '{:^9.6f}'.format(
242  mc_diff_eff *
243  100) +
244  " +- " +
245  '{:.6f}'.format(
246  mc_diff_eff_Uncertainty *
247  100) +
248  ' % | *')
249  print('* | | *')
250  print('* | B0-TAGGER TOTAL EFFECTIVE EFFICIENCIES: ' +
251  '{:.2f}'.format(mc_tot_eff_effB0 * 100) + " +-" + '{: 4.2f}'.format(mc_uncertainty_eff_effB0 * 100) +
252  ' % (q=+1) ' +
253  '{:.2f}'.format(mc_tot_eff_effB0bar * 100) + " +-" + '{: 4.2f}'.format(mc_uncertainty_eff_effB0bar * 100) +
254  ' % (q=-1) ' + ' | *')
255  print('* | | *')
256  print('* | FLAVOR PERCENTAGE (MC): ' +
257  '{:.2f}'.format(mc_total_tagged_B0 / mc_total_tagged * 100) + ' % (q=+1) ' +
258  '{:.2f}'.format(mc_total_tagged_B0bar / mc_total_tagged * 100) + ' % (q=-1) Diff=' +
259  '{:^5.2f}'.format((mc_total_tagged_B0 - mc_total_tagged_B0bar) / mc_total_tagged * 100) + ' % | *')
260  print('* |__________________________________________________________________________________________| *')
261  print('* *')
262  print('****************************************************************************************************')
263 
264  print(r'\begin{tabular}{ l r r r r r r r }')
265  print(r'\hline')
266  print(r'$r$- Interval & $\varepsilon_i\ $ & $w_i \pm \delta w_i\enskip\, $ ' +
267  r' & $\Delta w_i \pm \delta\Delta w_i $& $\varepsilon_{\text{eff}, i} \pm \delta\varepsilon_{\text{eff}, i}\enskip$ ' +
268  r' & & & $\Delta \varepsilon_{\text{eff}, i} \pm \delta\Delta \varepsilon_{\text{eff}, i} $\\ \hline\hline')
269  for i in range(0, r_size - 1):
270  print('$ ' + '{:.3f}'.format(r_subsample[i]) + ' - ' + '{:.3f}'.format(r_subsample[i + 1]) + '$ & $'
271  '{: 6.1f}'.format(mc_event_fractionTotal[i] * 100) + '$ & $' +
272  '{: 7.3f}'.format(mc_wvalue[i] * 100) + r" \pm " + '{:2.3f}'.format(mc_wvalueUncertainty[i] * 100) + r' $ & $' +
273  '{: 6.1f}'.format(mc_wvalueDiff[i] * 100) + r" \pm " +
274  '{:2.3f}'.format(mc_wvalueDiffUncertainty[i] * 100) + r'\, $ & $' +
275  '{: 8.4f}'.format(mc_iEffEfficiency[i] * 100) + # + '$ & $' +
276  r" \pm " + '{:2.4f}'.format(mc_iEffEfficiencyUncertainty[i] * 100) + r'\, $ & & & $' +
277  '{: 6.4f}'.format(mc_iDeltaEffEfficiency[i] * 100) + # +
278  r" \pm " + '{:2.4f}'.format(mc_iDeltaEffEfficiencyUncertainty[i] * 100) +
279  r'\enskip $ \\ ')
280  print(r'\hline\hline')
281  print(
282  r'\multicolumn{1}{r}{Total} & & \multicolumn{3}{r}{ $\varepsilon_\text{eff} = ' +
283  r'\sum_i \varepsilon_i \cdot \langle 1-2w_i\rangle^2 = ' +
284  '{: 6.2f}'.format(
285  mc_average_eff_eff *
286  100) +
287  r" \pm " +
288  '{: 6.2f}'.format(
289  mc_uncertainty_eff_effAverage *
290  100) +
291  r'\quad\,$ }')
292  print(r'& & \multicolumn{2}{r}{ $\Delta \varepsilon_\text{eff} = ' +
293  '{: 6.2f}'.format(mc_diff_eff * 100) + r" \pm " + '{: 6.2f}'.format(mc_diff_eff_Uncertainty * 100) + r'\ \quad\, $ }' +
294  r' \\')
295  print(r'\hline')
296  print(r'\end{tabular}')
297 
298  # Calculate Efficiency for Data
299 
300  data_totaldata = np.absolute(data[:, 0])
301  data_splotWeights = data[:, 1]
302 
303  data_entries = np.histogram(data_totaldata, rbins, weights=data_splotWeights)[0]
304  data_rvalueB0Average = (np.histogram(data_totaldata, rbins, weights=data_totaldata * data_splotWeights)[0] / data_entries)
305  data_rvalueMSB0Average = (
306  np.histogram(
307  data_totaldata,
308  rbins,
309  weights=(data_totaldata)**2 *
310  data_splotWeights)[0] /
311  data_entries)
312  data_rvalueStdB0Average = np.sqrt(abs(data_rvalueMSB0Average - data_rvalueB0Average * data_rvalueB0Average))
313 
314  data_total_tagged = 0
315  for iEntries in data_entries:
316  data_total_tagged += iEntries
317 
318  data_total_notTaggedWeighted = data_total_notTagged * data_total_tagged / data_totaldata.shape[0]
319 
320  data_total_entries = data_total_tagged + data_total_notTaggedWeighted
321  data_event_fractionTotal = data_entries.astype(float) / data_total_entries
322 
323  data_tagging_eff = data_total_tagged / data_total_entries
324  data_DeltaTagging_eff = math.sqrt(data_total_tagged * data_total_notTaggedWeighted**2 +
325  data_total_notTaggedWeighted * data_total_tagged**2) / (data_total_entries**2)
326  data_average_eff_eff = 0
327  data_uncertainty_eff_effAverage = 0
328  data_wvalue = array('f', [0] * r_size)
329  data_wvalueUncertainty = array('f', [0] * r_size)
330  data_iEffEfficiency = array('f', [0] * r_size)
331  data_iEffEfficiencyUncertainty = array('f', [0] * r_size)
332 
333  print('****************** MEASURED EFFECTIVE EFFICIENCY FOR COMBINER USING ' + method + ' ***************************')
334  print('* *')
335  print('* --> DETERMINATION BASED ON DATA *')
336  print('* *')
337  print('* ------------------------------------------------------------------------------------------------ *')
338  print('* r-interval <r> Efficiency w *')
339  print('* ------------------------------------------------------------------------------------------------ *')
340 
341  for i in range(0, r_size - 1):
342  data_wvalue[i] = (1 - data_rvalueB0Average[i]) / 2
343  data_wvalueUncertainty[i] = (data_rvalueStdB0Average[i] / math.sqrt(data_entries[i] - 1)) / 2
344  # fraction of events/all events
345 
346  data_iEffEfficiency[i] = data_tagging_eff * (data_event_fractionTotal[i] *
347  data_rvalueB0Average[i] * data_rvalueB0Average[i])
348 
349  data_iEffEfficiencyUncertainty[i] = data_rvalueB0Average[i] * \
350  math.sqrt((2 * data_total_entries * data_entries[i] *
351  (data_rvalueStdB0Average[i] / math.sqrt(data_entries[i] - 1)))**2 +
352  data_rvalueB0Average[i]**2 * data_entries[i] *
353  (data_total_entries * (data_total_entries - data_entries[i]) +
354  data_entries[i] * data_total_notTaggedWeighted)) / (data_total_entries**2)
355 
356  # finally calculating the total effective efficiency
357  data_average_eff_eff = data_average_eff_eff + \
358  data_event_fractionTotal[i] * data_rvalueB0Average[i] * data_rvalueB0Average[i]
359  data_uncertainty_eff_effAverage = data_uncertainty_eff_effAverage + data_iEffEfficiencyUncertainty[i]**2
360 
361  # intervallEff[i] = event_fractionTotal[i] * rvalueB0Average[i] * rvalueB0Average[i]
362  print('* ' + '{:.3f}'.format(r_subsample[i]) + ' - ' + '{:.3f}'.format(r_subsample[i + 1]) + ' ' +
363  '{:.3f}'.format(data_rvalueB0Average[i]) + ' +- ' + '{:.4f}'.format(data_rvalueStdB0Average[i]) + ' ' +
364  '{:.4f}'.format(data_event_fractionTotal[i]) + ' ' +
365  '{:.4f}'.format(data_wvalue[i]) + ' +- ' + '{:.4f}'.format(data_wvalueUncertainty[i]) + ' ' + ' *')
366 
367  data_uncertainty_eff_effAverage = math.sqrt(data_uncertainty_eff_effAverage)
368 
369  print('* ------------------------------------------------------------------------------------------------ *')
370  print('* *')
371  print('* __________________________________________________________________________________________ *')
372  print('* | | *')
373  print('* | TOTAL NUMBER OF TAGGED EVENTS = ' +
374  '{:<24}'.format("%.0f" % data_total_tagged) + '{:>36}'.format('| *'))
375  print('* | | *')
376  print(
377  '* | TOTAL AVERAGE EFFICIENCY (q=+-1)= ' +
378  '{:.2f}'.format(
379  data_tagging_eff *
380  100) +
381  " +- " +
382  '{:.2f}'.format(
383  data_DeltaTagging_eff *
384  100) +
385  ' % | *')
386  print('* | | *')
387  print(
388  '* | TOTAL AVERAGE EFFECTIVE EFFICIENCY (q=+-1)= ' +
389  '{:.6f}'.format(
390  data_average_eff_eff *
391  100) +
392  " +- " +
393  '{:.6f}'.format(
394  data_uncertainty_eff_effAverage *
395  100) +
396  ' % | *')
397  print('* |__________________________________________________________________________________________| *')
398  print('* *')
399  print('****************************************************************************************************')
400 
401  print(r'\begin{tabular}{ l r r r r}')
402  print(r'\hline')
403  print(r'$r$- Interval & $\varepsilon_i\ $ & $w_i \pm \delta w_i\ $ ' +
404  r' & $\varepsilon_{\text{eff}, i} \pm \delta\varepsilon_{\text{eff}, i}$ ' +
405  r'\\ \hline\hline')
406  for i in range(0, len(rbins) - 1):
407  print('$ ' + '{:.3f}'.format(r_subsample[i]) + ' - ' + '{:.3f}'.format(r_subsample[i + 1]) + '$ & $'
408  '{: 6.2f}'.format(data_event_fractionTotal[i] * 100) + '$ & $' +
409  '{: 6.2f}'.format(data_wvalue[i] * 100) + r" \pm " + '{:2.2f}'.format(data_wvalueUncertainty[i] * 100) + r'\, $ & $' +
410  '{: 7.3f}'.format(data_iEffEfficiency[i] * 100) + # + '$ & $' +
411  r" \pm " + '{:4.3f}'.format(data_iEffEfficiencyUncertainty[i] * 100) + r'\, $ \\ ')
412  print(r'\hline\hline')
413  print(
414  r'\multicolumn{1}{r}{Total} & \multicolumn{3}{r}{ $\varepsilon_\text{eff} = ' +
415  r'\sum_i \varepsilon_i \cdot \langle 1-2w_i\rangle^2 = ' +
416  '{: 6.1f}'.format(
417  data_average_eff_eff *
418  100) +
419  r" \pm " +
420  '{: 6.1f}'.format(
421  data_uncertainty_eff_effAverage *
422  100) +
423  r'\quad\,$ }' +
424  r' \\')
425  print(r'\hline')
426  print(r'\end{tabular}')
427 
428  # Plot of Signal for B2TIP
429 
430  # -----bins
431  bins = list(range(-25, 26, 1))
432  for i in range(0, len(bins)):
433  bins[i] = float(bins[i]) / 25
434  # ------
435 
436  fig1 = plt.figure(1, figsize=(11, 8))
437  ax1 = plt.axes([0.21, 0.15, 0.76, 0.8])
438 
439  qrDataPoints, qrDataBins = np.histogram(data[:, 0], bins=bins, weights=data_splotWeights)
440  qrDataPointsStd = np.sqrt(qrDataPoints)
441 
442  qrDataPointsNorm, qrDataBins = np.histogram(data[:, 0], bins=bins, weights=data_splotWeights, density=1)
443  qrDataPointsStdNorm = qrDataPointsStd * qrDataPointsNorm / qrDataPoints
444 
445  qrBincenters = 0.5 * (qrDataBins[1:] + qrDataBins[:-1])
446 
447  p3 = ax1.errorbar(qrBincenters, qrDataPointsNorm, xerr=0.02, yerr=qrDataPointsStdNorm, elinewidth=2, mew=2, ecolor='k',
448  fmt='o', mfc='k', markersize=6, label=r'${\rm Data}$')
449 
450  ax1.hist(mc[:, 0], bins=bins, density=1, histtype='step',
451  edgecolor='b', linewidth=4, linestyle='dashed', label=r'${\rm MC}$')
452 
453  p2, = ax1.plot([], label=r'${\rm MC}$', linewidth=4, linestyle='dashed', c='b')
454 
455  methodLabel = method
456  if method == "FANN":
457  methodLabel = "MLP"
458 
459  ax1.set_ylabel(r'${\rm Fraction\ \ of\ \ Events\ /\ (\, 0.04\, )}$', fontsize=35)
460  ax1.set_xlabel(r'$(q\cdot r)_{\rm ' + methodLabel + '}$', fontsize=35)
461  plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1],
462  [r'$-1$', r'', r'$-0.5$', r'', r'$0$', r'', r'$0.5$', r'', r'$1$'], rotation=0, size=40)
463  ax1.tick_params(axis='y', labelsize=35)
464  ax1.legend([p3, p2], [r'${\rm Data}$', r'${\rm MC}$'], prop={'size': 35}, loc=0, numpoints=1)
465  ax1.grid(True)
466 
467  if qrDataPointsNorm[np.argmax(qrDataPointsNorm)] < 1.4:
468  ax1.set_ylim(0, 1.4)
469 
470  ax1.set_xlim(-1.005, 1.005)
471  plt.savefig(PATH + '/QR_B0_B0bar' + methodLabel + '.pdf')
472  fig1.clear()
473 
474  # Plot of Continuum for B2TIP and Thesis
475 
476  fig2 = plt.figure(2, figsize=(11, 8))
477  ax1 = plt.axes([0.21, 0.15, 0.76, 0.8])
478 
479  qrDataPointsCont, qrDataBins = np.histogram(data[:, 0], bins=bins, weights=data[:, 2])
480  qrDataPointsContStd = np.sqrt(qrDataPointsCont)
481 
482  qrDataPointsContNorm, qrDataBins = np.histogram(data[:, 0], bins=bins, weights=data[:, 2], density=1)
483  qrDataPointsContStdNorm = qrDataPointsContStd * qrDataPointsContNorm / qrDataPointsCont
484 
485  qrBincenters = 0.5 * (qrDataBins[1:] + qrDataBins[:-1])
486 
487  p3 = ax1.errorbar(qrBincenters, qrDataPointsContNorm, xerr=0.02, yerr=qrDataPointsContStdNorm, elinewidth=2, mew=2, ecolor='k',
488  fmt='o', mfc='k', mec='#990000', markersize=6, label=r'${\rm Data}$')
489 
490  ax1.set_ylabel(r'${\rm Fraction\ \ of\ \ Events\ /\ (\, 0.04\, )}$', fontsize=35)
491  ax1.set_xlabel(r'$(q\cdot r)_{\rm ' + methodLabel + '}$', fontsize=35)
492  plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1],
493  [r'$-1$', r'', r'$-0.5$', r'', r'$0$', r'', r'$0.5$', r'', r'$1$'], rotation=0, size=40)
494  ax1.tick_params(axis='y', labelsize=35)
495  ax1.legend([p3], [r'${\rm Data}$'], prop={'size': 35}, loc=1, numpoints=1)
496  ax1.grid(True)
497 
498  if qrDataPointsContNorm[np.argmax(qrDataPointsContNorm)] < 1.4:
499  ax1.set_ylim(0, 1.4)
500 
501  ax1.set_xlim(-1.005, 1.005)
502  plt.savefig(PATH + '/QR_B0_B0bar' + methodLabel + 'Continuum.pdf')
503  fig2.clear()
504 
505  data_total_tagged_Continuum = 0
506  for iEntries in qrDataPointsCont:
507  data_total_tagged_Continuum += iEntries
508 
509  print('****************************************************************************************************')
510  print('* TOTAL NUMBER OF TAGGED CONTINUUM EVENTS = ' +
511  '{:<24}'.format("%.0f" % data_total_tagged_Continuum) + '{:>28}'.format(' *'))
512 
513  # Chi-squared between MC and Data
514 
515  mcHistW = ROOT.TH1F("mcHistW", "mcHistW", int(r_size - 2), r_subsample)
516  dataHistW = ROOT.TH1F("dataHistW", "dataHistW", int(r_size - 2), r_subsample)
517 
518  mcHistEff = ROOT.TH1F("mcHistEff", "mcHistEff", int(r_size - 2), r_subsample)
519  dataHistEff = ROOT.TH1F("dataHistEff", "dataHistEff", int(r_size - 2), r_subsample)
520 
521  print('****************************************************************************************************')
522  print('* *')
523  print('* --> TEST OF COMPATIBILITY BETWEEN DATA AND MC *')
524 
525  for i in range(0, len(rbins) - 1):
526 
527  mcHistW.SetBinContent(i + 1, mc_wvalue[i])
528  mcHistW.SetBinError(i + 1, mc_wvalueUncertainty[i])
529 
530  dataHistW.SetBinContent(i + 1, data_wvalue[i])
531  dataHistW.SetBinError(i + 1, data_wvalueUncertainty[i])
532 
533  mcHistEff.SetBinContent(i + 1, mc_iEffEfficiency[i])
534  mcHistEff.SetBinError(i + 1, mc_iEffEfficiencyUncertainty[i])
535 
536  dataHistEff.SetBinContent(i + 1, data_iEffEfficiency[i])
537  dataHistEff.SetBinError(i + 1, data_iEffEfficiencyUncertainty[i])
538 
539  residualsW = array('d', [0] * r_size)
540  residualsEff = array('d', [0] * r_size)
541  print('* *')
542  print('* TEST ON w *')
543  print('* *')
544 
545  chi2W = dataHistW.Chi2Test(mcHistW, "WW CHI2 P", residualsW)
546  print("1-CL for NDF = 7 is = ", ROOT.TMath.Prob(chi2W, 7))
547  print('* *')
548  print('* TEST ON EFFECTIVE EFFICIENCY *')
549  print('* *')
550  chi2Eff = dataHistEff.Chi2Test(mcHistEff, "WW P", residualsEff)
551  print("1-CL for NDF = 7 is = ", ROOT.TMath.Prob(chi2Eff, 7))
552 
553  print('****************************************************************************************************')
554 
555  # Plot of Signal with Residuals for Thesis
556 
557  qrMCPoints, qrDataBins = np.histogram(mc[:, 0], bins=bins)
558  qrMCPointsStd = np.sqrt(qrMCPoints)
559  qrMCPointsNorm, qrDataBins = np.histogram(mc[:, 0], bins=bins, density=1)
560  qrMCPointsStdNorm = qrMCPointsStd * qrMCPointsNorm / qrMCPoints
561 
562  qrDataMCResiduals = (qrDataPointsNorm - qrMCPointsNorm) \
563  / np.sqrt(qrDataPointsStdNorm * qrDataPointsStdNorm + qrMCPointsStdNorm * qrMCPointsStdNorm)
564 
565  fig3 = plt.figure(3, figsize=(11, 11))
566  ax1 = plt.axes([0.21, 0.37, 0.76, 0.60])
567  # print('Here Log Scale')
568  # ax1.set_yscale('log', nonposy='clip') # Only for Semileptonic
569  p3 = ax1.errorbar(qrBincenters, qrDataPointsNorm, xerr=0.02, yerr=qrDataPointsStdNorm,
570  elinewidth=2, mew=2, ecolor='k',
571  fmt='o', mfc='k', mec='#006600', markersize=6, label=r'${\rm Data}$')
572 
573  ax1.hist(mc[:, 0], bins=bins, density=1, histtype='step',
574  edgecolor='b', linewidth=4, linestyle='dashed', label=r'${\rm MC}$')
575 
576  p2, = ax1.plot([], label=r'${\rm MC}$', linewidth=4, linestyle='dashed', c='b')
577 
578  ax1.set_ylabel(r'${\rm Fraction\ \ of\ \ Events\ /\ (\, 0.04\, )}$', fontsize=35, labelpad=20)
579  # ax1.set_xlabel(r'', fontsize=35)
580  ax1.tick_params(axis='y', labelsize=35)
581  ax1.legend([p3, p2], [r'${\rm Data}$', r'${\rm MC}$'], prop={'size': 35}, loc=1, numpoints=1)
582  ax1.grid(True)
583 
584  if qrDataPointsNorm[np.argmax(qrDataPointsNorm)] < 1.4:
585  ax1.set_ylim(0, 1.4)
586  plt.yticks([0, 0.25, 0.5, 0.75, 1.0, 1.25],
587  [r'$0.00$', r'$0.25$', r'$0.5$', r'$0.75$', r'$1.00$', r'$1.25$'], rotation=0, size=35)
588  ax1.set_xlim(-1.005, 1.005)
589  plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1],
590  [r'', r'', r'', r'', r'', r'', r'', r'', r''], rotation=0, size=35)
591 
592  ax2 = plt.axes([0.21, 0.15, 0.76, 0.2])
593  ax2.errorbar(qrBincenters, qrDataMCResiduals, xerr=0.02, yerr=1, elinewidth=2, mew=2, ecolor='k',
594  fmt='o', mfc='k', mec='#006600', markersize=6, label=r'${\rm Data}$')
595  ax2.set_ylabel(r'${\rm Normalized}$' + '\n' + r'${\rm Residuals}$', fontsize=35, labelpad=20)
596  ax2.set_xlabel(r'$(q\cdot r)_{\rm ' + methodLabel + '}$', fontsize=35)
597  plt.xticks([-1, -0.75, -0.5, -0.25, 0, 0.25, 0.5, 0.75, 1],
598  [r'$-1$', r'', r'$-0.5$', r'', r'$0$', r'', r'$0.5$', r'', r'$1$'], rotation=0, size=35)
599  plt.yticks([-5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5],
600  [r'', r'$-4$', r'', r'$-2$', r'', r'$0$', r'', r'$2$', r'', r'$4$', r''], rotation=0, size=25)
601  ax2.set_ylim(-5, 5)
602  # ax2.tick_params(axis='y', labelsize=30)
603  ax2.set_xlim(-1.005, 1.005)
604  ax2.xaxis.grid(True) # linestyle='--'
605  # plt.axhline(y= 4, xmin=-1.005, xmax=1.005, linewidth=1, color = 'k', linestyle = '-')
606  plt.axhline(y=3, xmin=-1.005, xmax=1.005, linewidth=2, color='#006600', linestyle='-')
607  # plt.axhline(y= 2, xmin=-1.005, xmax=1.005, linewidth=1, color = 'k', linestyle = '-')
608  plt.axhline(y=1, xmin=-1.005, xmax=1.005, linewidth=2, color='b', linestyle='-')
609  plt.axhline(y=-1, xmin=-1.005, xmax=1.005, linewidth=2, color='b', linestyle='-')
610  # plt.axhline(y=-2, xmin=-1.005, xmax=1.005, linewidth=1, color = 'k', linestyle = '-')
611  plt.axhline(y=-3, xmin=-1.005, xmax=1.005, linewidth=2, color='#006600', linestyle='-')
612  # plt.axhline(y=-4, xmin=-1.005, xmax=1.005, linewidth=1, color = 'k', linestyle = '-')
613  plt.savefig(PATH + '/QR_B0_B0bar' + methodLabel + 'WithRes.pdf')
614  fig3.clear()
615 
616  return 1
617 
618 
619 def plotWithResiduals(rooFitFrame, rooRealVar, dots, modelCurve, units, nameOfPlot, legend, removeArtifacts):
620 
621  # plot for Thesis
622  rooFitFrame.Print()
623  rooFitFrame.GetXaxis().SetTitle("")
624  rooFitFrame.GetXaxis().SetLabelSize(0)
625 
626  rooFitFrame.GetYaxis().SetTitleSize(0.072)
627  rooFitFrame.GetYaxis().SetTitleOffset(0.98)
628  rooFitFrame.GetYaxis().SetLabelSize(0.055)
629 
630  yLabelBin = 0
631  yLabelBin = '{:0.0f}'.format(float(dots.GetErrorXlow(1) + dots.GetErrorXhigh(1)))
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  l = ROOT.TLegend(0.25, 0.43, 0.5, 0.89)
1036  l.AddEntry(dots, "Belle data")
1037  l.AddEntry(modelCurve, 'Fit result')
1038  l.AddEntry(signalCurve, "Signal")
1039  l.AddEntry(continuumCurve, 'Continuum')
1040 
1041  l.SetTextSize(0.055)
1042 
1043  plotWithResiduals(mbcFrame, B0_mbc, dots, modelCurve, "[GeV/#it{c}^{2}]", "Mbc", l, 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")