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