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