Belle II Software  release-06-02-00
KLMK0LPlotModule.py
1 #!/usr/bin/env python3
2 # -*- coding: utf-8 -*-
3 
4 
11 
12 """
13 <header>
14  <contact>Kirill Chilikin (chilikin@lebedev.ru)</contact>
15  <description>Creation of KLM K0L validation plots.</description>
16 </header>
17 """
18 
19 import basf2
20 import ROOT
21 from ROOT import Belle2
22 from ROOT import TNamed
23 import math
24 import numpy
25 
26 
27 class KLMK0LPlotModule(basf2.Module):
28  """ Class for creation of KLM K0L validation plots. """
29 
30  def __init__(self, output_file, evtgen, check_eklm):
31  """Initialization."""
32  super(KLMK0LPlotModule, self).__init__()
33 
34  self.evtgenevtgen = evtgen
35 
36  self.check_eklmcheck_eklm = check_eklm
37 
38  self.output_fileoutput_file = ROOT.TFile(output_file, 'recreate')
39  contact = 'Kirill Chilikin (chilikin@lebedev.ru)'
40 
41  self.hist_nklhist_nkl = ROOT.TH1F('k0l_number',
42  'Number of KLM clusters per 1 MC particle',
43  5, -0.5, 4.5)
44  self.hist_nklhist_nkl.SetXTitle('KLM clusters')
45  self.hist_nklhist_nkl.SetYTitle('Events')
46  functions = self.hist_nklhist_nkl.GetListOfFunctions()
47  functions.Add(TNamed('Description', 'Number of KLM clusters per 1 MC particle'))
48  functions.Add(TNamed('Check', 'No efficiency decrease or multiple candidates \
49  increase'))
50  functions.Add(TNamed('Contact', contact))
51  functions.Add(TNamed('MetaOptions', 'shifter'))
52 
53  self.hist_xreshist_xres = ROOT.TH1F('k0l_xres',
54  'KLM K0L decay vertex X resolution',
55  100, -50, 50)
56  self.hist_xreshist_xres.SetXTitle('cm')
57  self.hist_xreshist_xres.SetYTitle('Events')
58  functions = self.hist_xreshist_xres.GetListOfFunctions()
59  functions.Add(TNamed('Description', 'X resolution'))
60  functions.Add(TNamed('Check', 'No bias, resolution ~ 16 cm.'))
61  functions.Add(TNamed('Contact', contact))
62  functions.Add(TNamed('MetaOptions', 'shifter'))
63 
64  self.hist_yreshist_yres = ROOT.TH1F('k0l_yres',
65  'KLM K0L decay vertex Y resolution',
66  100, -50, 50)
67  self.hist_yreshist_yres.SetXTitle('cm')
68  self.hist_yreshist_yres.SetYTitle('Events')
69  functions = self.hist_yreshist_yres.GetListOfFunctions()
70  functions.Add(TNamed('Description', 'Y resolution'))
71  functions.Add(TNamed('Check', 'No bias, resolution ~ 16 cm.'))
72  functions.Add(TNamed('Contact', contact))
73  functions.Add(TNamed('MetaOptions', 'shifter'))
74 
75  self.hist_zreshist_zres = ROOT.TH1F('k0l_zres',
76  'KLM K0L decay vertex Z resolution',
77  100, -50, 50)
78  self.hist_zreshist_zres.SetXTitle('cm')
79  self.hist_zreshist_zres.SetYTitle('Events')
80  functions = self.hist_zreshist_zres.GetListOfFunctions()
81  functions.Add(TNamed('Description', 'Z resolution'))
82  functions.Add(TNamed('Check', 'No bias, resolution ~ 16 cm.'))
83  functions.Add(TNamed('Contact', contact))
84  functions.Add(TNamed('MetaOptions', 'shifter'))
85 
86  self.hist_treshist_tres = ROOT.TH1F('k0l_tres',
87  'KLM K0L decay time resolution',
88  100, -20., 10.)
89  self.hist_treshist_tres.SetXTitle('ns')
90  self.hist_treshist_tres.SetYTitle('Events')
91  functions = self.hist_treshist_tres.GetListOfFunctions()
92  functions.Add(TNamed('Description', 'Time resolution'))
93  functions.Add(TNamed('Check', 'No bias.'))
94  functions.Add(TNamed('Contact', contact))
95  functions.Add(TNamed('MetaOptions', 'shifter'))
96 
97  self.hist_preshist_pres = ROOT.TH1F('k0l_pres',
98  'KLM K0L momentum resolution',
99  100, -3., 3.)
100  self.hist_preshist_pres.SetXTitle('GeV')
101  self.hist_preshist_pres.SetYTitle('Events')
102  functions = self.hist_preshist_pres.GetListOfFunctions()
103  functions.Add(TNamed('Description', 'Momentum resolution'))
104  functions.Add(TNamed('Check', 'No bias.'))
105  functions.Add(TNamed('Contact', contact))
106  functions.Add(TNamed('MetaOptions', 'shifter'))
107 
108  self.hist_ptreshist_ptres = ROOT.TH1F('k0l_ptres',
109  'KLM K0L momentum theta resolution',
110  100, -0.2, 0.2)
111  self.hist_ptreshist_ptres.SetXTitle('rad')
112  self.hist_ptreshist_ptres.SetYTitle('Events')
113  functions = self.hist_ptreshist_ptres.GetListOfFunctions()
114  functions.Add(TNamed('Description', 'Momentum theta resolution'))
115  functions.Add(TNamed('Check', 'No bias, resolution ~ 0.06'))
116  functions.Add(TNamed('Contact', contact))
117  functions.Add(TNamed('MetaOptions', 'shifter'))
118 
119  self.hist_ppreshist_ppres = ROOT.TH1F('k0l_ppres',
120  'KLM K0L momentum phi resolution',
121  100, -0.2, 0.2)
122  self.hist_ppreshist_ppres.SetXTitle('rad')
123  self.hist_ppreshist_ppres.SetYTitle('Events')
124  functions = self.hist_ppreshist_ppres.GetListOfFunctions()
125  functions.Add(TNamed('Description', 'Momentum phi resolution'))
126  functions.Add(TNamed('Check', 'No bias, resolution ~ 0.07'))
127  functions.Add(TNamed('Contact', contact))
128  functions.Add(TNamed('MetaOptions', 'shifter'))
129 
130  self.hist_covmathist_covmat = ROOT.TH1F('k0l_covmat',
131  'KLM K0L coordinates covariance matrix',
132  6, 0, 1)
133  self.hist_covmathist_covmat.GetXaxis().SetBinLabel(1, 'xx')
134  self.hist_covmathist_covmat.GetXaxis().SetBinLabel(2, 'xy')
135  self.hist_covmathist_covmat.GetXaxis().SetBinLabel(3, 'xz')
136  self.hist_covmathist_covmat.GetXaxis().SetBinLabel(4, 'yy')
137  self.hist_covmathist_covmat.GetXaxis().SetBinLabel(5, 'yz')
138  self.hist_covmathist_covmat.GetXaxis().SetBinLabel(6, 'zz')
139  self.hist_covmathist_covmat.SetYTitle('Covariance, cm^{2}')
140  functions = self.hist_covmathist_covmat.GetListOfFunctions()
141  functions.Add(TNamed('Description', 'Momentum phi resolution'))
142  functions.Add(TNamed('Check', 'No large off-diagonal elements.'))
143  functions.Add(TNamed('Contact', contact))
144  functions.Add(TNamed('MetaOptions', 'shifter'))
145 
146  self.hist_corrmathist_corrmat = ROOT.TH1F('k0l_corrmat',
147  'KLM K0L correlation matrix',
148  10, 0, 1)
149  self.hist_corrmathist_corrmat.GetXaxis().SetBinLabel(1, 'xx')
150  self.hist_corrmathist_corrmat.GetXaxis().SetBinLabel(2, 'xy')
151  self.hist_corrmathist_corrmat.GetXaxis().SetBinLabel(3, 'xz')
152  self.hist_corrmathist_corrmat.GetXaxis().SetBinLabel(4, 'xp')
153  self.hist_corrmathist_corrmat.GetXaxis().SetBinLabel(5, 'yy')
154  self.hist_corrmathist_corrmat.GetXaxis().SetBinLabel(6, 'yz')
155  self.hist_corrmathist_corrmat.GetXaxis().SetBinLabel(7, 'yp')
156  self.hist_corrmathist_corrmat.GetXaxis().SetBinLabel(8, 'zz')
157  self.hist_corrmathist_corrmat.GetXaxis().SetBinLabel(9, 'zp')
158  self.hist_corrmathist_corrmat.GetXaxis().SetBinLabel(10, 'pp')
159  self.hist_corrmathist_corrmat.SetYTitle('Correlation coefficient')
160  functions = self.hist_corrmathist_corrmat.GetListOfFunctions()
161  functions.Add(TNamed('Description', 'Momentum phi resolution'))
162  functions.Add(TNamed('Check', 'No large off-diagonal elements.'))
163  functions.Add(TNamed('Contact', contact))
164  functions.Add(TNamed('MetaOptions', 'shifter'))
165 
166  self.vertex_k_avvertex_k_av = ROOT.TVector3(0, 0, 0)
167 
168  self.vertexvertex = []
169 
170  self.momentum_avmomentum_av = 0
171 
172  self.momentummomentum = []
173 
174  def event(self):
175  """ Event function. """
176  mc_particles = Belle2.PyStoreArray('MCParticles')
177  for mc_particle in mc_particles:
178  # Select K_L0.
179  if (mc_particle.getPDG() != 130):
180  continue
181  # Select primary K_L0.
182  if (self.evtgenevtgen):
183  b_pdg = abs(mc_particle.getMother().getPDG())
184  if (not(b_pdg == 511)):
185  continue
186  else:
187  if (mc_particle.getProductionTime() > 0):
188  continue
189  vertex = mc_particle.getDecayVertex()
190  time = mc_particle.getDecayTime()
191  momentum = mc_particle.getMomentum()
192  if (self.check_eklmcheck_eklm):
193  x = vertex.x()
194  y = vertex.y()
195  z = vertex.z()
196  r = math.sqrt(x * x + y * y)
197  if (r < 132.5 or r > 329.0):
198  continue
199  if (abs(x) < 8.2 or abs(y) < 8.2):
200  continue
201  if (not((z > -315.1 and z < -183.0) or
202  (z > 277.0 and z < 409.1))):
203  continue
204  klm_clusters = mc_particle.getRelationsFrom('KLMClusters')
205  self.hist_nklhist_nkl.Fill(len(klm_clusters))
206  for klm_cluster in klm_clusters:
207  vertex_k = klm_cluster.getClusterPosition() - vertex
208  self.vertexvertex.append(vertex_k)
209  self.vertex_k_avvertex_k_av = self.vertex_k_avvertex_k_av + vertex_k
210  self.momentummomentum.append(klm_cluster.getMomentumMag())
211  self.momentum_avmomentum_av = self.momentum_avmomentum_av + \
212  klm_cluster.getMomentumMag()
213  time_k = klm_cluster.getTime()
214  momentum_k = klm_cluster.getMomentum().Vect()
215  self.hist_xreshist_xres.Fill(vertex_k.x())
216  self.hist_yreshist_yres.Fill(vertex_k.y())
217  self.hist_zreshist_zres.Fill(vertex_k.z())
218  self.hist_treshist_tres.Fill(time_k - time)
219  self.hist_preshist_pres.Fill(momentum_k.Mag() - momentum.Mag())
220  self.hist_ptreshist_ptres.Fill(momentum_k.Theta() - momentum.Theta())
221  self.hist_ppreshist_ppres.Fill(momentum_k.Phi() - momentum.Phi())
222 
223  def terminate(self):
224  """ Termination function. """
225  self.vertex_k_avvertex_k_av = self.vertex_k_avvertex_k_av * (1.0 / len(self.vertexvertex))
226  self.momentum_avmomentum_av = self.momentum_avmomentum_av * (1.0 / len(self.vertexvertex))
227  # x, y, z, e
228  cov_mat = numpy.zeros((4, 4))
229  cov_mat_err = numpy.zeros((4, 4))
230  corr_mat = numpy.zeros((4, 4))
231  corr_mat_err = numpy.zeros((4, 4))
232  for i in range(len(self.vertexvertex)):
233  cov_mat[0][0] = cov_mat[0][0] + \
234  (self.vertexvertex[i].x() - self.vertex_k_avvertex_k_av.x()) * \
235  (self.vertexvertex[i].x() - self.vertex_k_avvertex_k_av.x())
236  cov_mat[0][1] = cov_mat[0][1] + \
237  (self.vertexvertex[i].x() - self.vertex_k_avvertex_k_av.x()) * \
238  (self.vertexvertex[i].y() - self.vertex_k_avvertex_k_av.y())
239  cov_mat[0][2] = cov_mat[0][2] + \
240  (self.vertexvertex[i].x() - self.vertex_k_avvertex_k_av.x()) * \
241  (self.vertexvertex[i].z() - self.vertex_k_avvertex_k_av.z())
242  cov_mat[0][3] = cov_mat[0][3] + \
243  (self.vertexvertex[i].x() - self.vertex_k_avvertex_k_av.x()) * \
244  (self.momentummomentum[i] - self.momentum_avmomentum_av)
245  cov_mat[1][1] = cov_mat[1][1] + \
246  (self.vertexvertex[i].y() - self.vertex_k_avvertex_k_av.y()) * \
247  (self.vertexvertex[i].y() - self.vertex_k_avvertex_k_av.y())
248  cov_mat[1][2] = cov_mat[1][2] + \
249  (self.vertexvertex[i].y() - self.vertex_k_avvertex_k_av.y()) * \
250  (self.vertexvertex[i].z() - self.vertex_k_avvertex_k_av.z())
251  cov_mat[1][3] = cov_mat[1][3] + \
252  (self.vertexvertex[i].y() - self.vertex_k_avvertex_k_av.y()) * \
253  (self.momentummomentum[i] - self.momentum_avmomentum_av)
254  cov_mat[2][2] = cov_mat[2][2] + \
255  (self.vertexvertex[i].z() - self.vertex_k_avvertex_k_av.z()) * \
256  (self.vertexvertex[i].z() - self.vertex_k_avvertex_k_av.z())
257  cov_mat[2][3] = cov_mat[2][3] + \
258  (self.vertexvertex[i].z() - self.vertex_k_avvertex_k_av.z()) * \
259  (self.momentummomentum[i] - self.momentum_avmomentum_av)
260  cov_mat[3][3] = cov_mat[3][3] + \
261  (self.momentummomentum[i] - self.momentum_avmomentum_av) * \
262  (self.momentummomentum[i] - self.momentum_avmomentum_av)
263  for i in range(0, 4):
264  for j in range(i, 4):
265  cov_mat[i][j] = cov_mat[i][j] / (len(self.vertexvertex) - 1)
266  for i in range(len(self.vertexvertex)):
267  cov_mat_err[0][0] = cov_mat_err[0][0] + \
268  pow((self.vertexvertex[i].x() - self.vertex_k_avvertex_k_av.x()) *
269  (self.vertexvertex[i].x() - self.vertex_k_avvertex_k_av.x()) - cov_mat[0][0], 2)
270  cov_mat_err[0][1] = cov_mat_err[0][1] + \
271  pow((self.vertexvertex[i].x() - self.vertex_k_avvertex_k_av.x()) *
272  (self.vertexvertex[i].y() - self.vertex_k_avvertex_k_av.y()) - cov_mat[0][1], 2)
273  cov_mat_err[0][2] = cov_mat_err[0][2] + \
274  pow((self.vertexvertex[i].x() - self.vertex_k_avvertex_k_av.x()) *
275  (self.vertexvertex[i].z() - self.vertex_k_avvertex_k_av.z()) - cov_mat[0][2], 2)
276  cov_mat_err[0][3] = cov_mat_err[0][3] + \
277  pow((self.vertexvertex[i].x() - self.vertex_k_avvertex_k_av.x()) *
278  (self.momentummomentum[i] - self.momentum_avmomentum_av) - cov_mat[0][3], 2)
279  cov_mat_err[1][1] = cov_mat_err[1][1] + \
280  pow((self.vertexvertex[i].y() - self.vertex_k_avvertex_k_av.y()) *
281  (self.vertexvertex[i].y() - self.vertex_k_avvertex_k_av.y()) - cov_mat[1][1], 2)
282  cov_mat_err[1][2] = cov_mat_err[1][2] + \
283  pow((self.vertexvertex[i].y() - self.vertex_k_avvertex_k_av.y()) *
284  (self.vertexvertex[i].z() - self.vertex_k_avvertex_k_av.z()) - cov_mat[1][2], 2)
285  cov_mat_err[1][3] = cov_mat_err[1][3] + \
286  pow((self.vertexvertex[i].y() - self.vertex_k_avvertex_k_av.y()) *
287  (self.momentummomentum[i] - self.momentum_avmomentum_av) - cov_mat[1][3], 2)
288  cov_mat_err[2][2] = cov_mat_err[2][2] + \
289  pow((self.vertexvertex[i].z() - self.vertex_k_avvertex_k_av.z()) *
290  (self.vertexvertex[i].z() - self.vertex_k_avvertex_k_av.z()) - cov_mat[2][2], 2)
291  cov_mat_err[2][3] = cov_mat_err[2][3] + \
292  pow((self.vertexvertex[i].z() - self.vertex_k_avvertex_k_av.z()) *
293  (self.momentummomentum[i] - self.momentum_avmomentum_av) - cov_mat[2][3], 2)
294  cov_mat_err[3][3] = cov_mat_err[3][3] + \
295  pow((self.momentummomentum[i] - self.momentum_avmomentum_av) *
296  (self.momentummomentum[i] - self.momentum_avmomentum_av) - cov_mat[3][3], 2)
297  for i in range(0, 4):
298  for j in range(i, 4):
299  cov_mat_err[i][j] = \
300  math.sqrt(cov_mat_err[i][j] /
301  ((len(self.vertexvertex) - 2) * len(self.vertexvertex)))
302  self.hist_covmathist_covmat.SetBinContent(1, cov_mat[0][0])
303  self.hist_covmathist_covmat.SetBinContent(2, cov_mat[0][1])
304  self.hist_covmathist_covmat.SetBinContent(3, cov_mat[0][2])
305  self.hist_covmathist_covmat.SetBinContent(4, cov_mat[1][1])
306  self.hist_covmathist_covmat.SetBinContent(5, cov_mat[1][2])
307  self.hist_covmathist_covmat.SetBinContent(6, cov_mat[2][2])
308  self.hist_covmathist_covmat.SetBinError(1, cov_mat_err[0][0])
309  self.hist_covmathist_covmat.SetBinError(2, cov_mat_err[0][1])
310  self.hist_covmathist_covmat.SetBinError(3, cov_mat_err[0][2])
311  self.hist_covmathist_covmat.SetBinError(4, cov_mat_err[1][1])
312  self.hist_covmathist_covmat.SetBinError(5, cov_mat_err[1][2])
313  self.hist_covmathist_covmat.SetBinError(6, cov_mat_err[2][2])
314  for i in range(0, 4):
315  for j in range(i, 4):
316  corr_mat[i][j] = cov_mat[i][j] / \
317  math.sqrt(cov_mat[i][i]) / math.sqrt(cov_mat[j][j])
318  # Normalization error is not taken into account.
319  corr_mat_err[i][j] = cov_mat_err[i][j] / cov_mat[i][j] * \
320  corr_mat[i][j]
321  self.hist_corrmathist_corrmat.SetBinContent(1, corr_mat[0][0])
322  self.hist_corrmathist_corrmat.SetBinContent(2, corr_mat[0][1])
323  self.hist_corrmathist_corrmat.SetBinContent(3, corr_mat[0][2])
324  self.hist_corrmathist_corrmat.SetBinContent(4, corr_mat[0][3])
325  self.hist_corrmathist_corrmat.SetBinContent(5, corr_mat[1][1])
326  self.hist_corrmathist_corrmat.SetBinContent(6, corr_mat[1][2])
327  self.hist_corrmathist_corrmat.SetBinContent(7, corr_mat[1][3])
328  self.hist_corrmathist_corrmat.SetBinContent(8, corr_mat[2][2])
329  self.hist_corrmathist_corrmat.SetBinContent(9, corr_mat[2][3])
330  self.hist_corrmathist_corrmat.SetBinContent(10, corr_mat[3][3])
331  self.hist_corrmathist_corrmat.SetBinError(1, corr_mat_err[0][0])
332  self.hist_corrmathist_corrmat.SetBinError(2, corr_mat_err[0][1])
333  self.hist_corrmathist_corrmat.SetBinError(3, corr_mat_err[0][2])
334  self.hist_corrmathist_corrmat.SetBinError(4, corr_mat_err[0][3])
335  self.hist_corrmathist_corrmat.SetBinError(5, corr_mat_err[1][1])
336  self.hist_corrmathist_corrmat.SetBinError(6, corr_mat_err[1][2])
337  self.hist_corrmathist_corrmat.SetBinError(7, corr_mat_err[1][3])
338  self.hist_corrmathist_corrmat.SetBinError(8, corr_mat_err[2][2])
339  self.hist_corrmathist_corrmat.SetBinError(9, corr_mat_err[2][3])
340  self.hist_corrmathist_corrmat.SetBinError(10, corr_mat_err[3][3])
341  self.output_fileoutput_file.cd()
342  self.hist_nklhist_nkl.Write()
343  self.hist_xreshist_xres.Write()
344  self.hist_yreshist_yres.Write()
345  self.hist_zreshist_zres.Write()
346  self.hist_treshist_tres.Write()
347  self.hist_preshist_pres.Write()
348  self.hist_ptreshist_ptres.Write()
349  self.hist_ppreshist_ppres.Write()
350  self.hist_covmathist_covmat.Write()
351  self.hist_corrmathist_corrmat.Write()
352  self.output_fileoutput_file.Close()
Belle2::PyStoreArray
a (simplified) python wrapper for StoreArray.
Definition: PyStoreArray.h:56
KLMK0LPlotModule.KLMK0LPlotModule.hist_corrmat
hist_corrmat
Correlation matrix histogram.
Definition: KLMK0LPlotModule.py:146
KLMK0LPlotModule.KLMK0LPlotModule.hist_tres
hist_tres
Time resolution histogram.
Definition: KLMK0LPlotModule.py:86
KLMK0LPlotModule.KLMK0LPlotModule.hist_nkl
hist_nkl
Number of K0L histogram.
Definition: KLMK0LPlotModule.py:41
KLMK0LPlotModule.KLMK0LPlotModule.check_eklm
check_eklm
Whether to check if cluster is in EKLM.
Definition: KLMK0LPlotModule.py:36
KLMK0LPlotModule.KLMK0LPlotModule.evtgen
evtgen
True for evtgen events, false for particle gun.
Definition: KLMK0LPlotModule.py:34
KLMK0LPlotModule.KLMK0LPlotModule.hist_pres
hist_pres
Momentum resolution histogram.
Definition: KLMK0LPlotModule.py:97
KLMK0LPlotModule.KLMK0LPlotModule.hist_xres
hist_xres
X resolution histogram.
Definition: KLMK0LPlotModule.py:53
KLMK0LPlotModule.KLMK0LPlotModule.hist_zres
hist_zres
Z resolution histogram.
Definition: KLMK0LPlotModule.py:75
KLMK0LPlotModule.KLMK0LPlotModule.hist_yres
hist_yres
Y resolution histogram.
Definition: KLMK0LPlotModule.py:64
KLMK0LPlotModule.KLMK0LPlotModule.__init__
def __init__(self, output_file, evtgen, check_eklm)
Definition: KLMK0LPlotModule.py:30
KLMK0LPlotModule.KLMK0LPlotModule.hist_covmat
hist_covmat
Covariance matrix histogram.
Definition: KLMK0LPlotModule.py:130
KLMK0LPlotModule.KLMK0LPlotModule.hist_ptres
hist_ptres
Momentum theta resolution histogram.
Definition: KLMK0LPlotModule.py:108
KLMK0LPlotModule.KLMK0LPlotModule.hist_ppres
hist_ppres
Momentum phi resolution histogram.
Definition: KLMK0LPlotModule.py:119