Belle II Software  release-08-01-10
KLMK0LPlotModule.py
1 #!/usr/bin/env python3
2 
3 
10 
11 """
12 <header>
13  <noexecute>Definition of plotting class</noexecute>
14  <description>Creation of KLM K0L validation plots.</description>
15 </header>
16 """
17 
18 import basf2
19 import ROOT
20 from ROOT import Belle2
21 from ROOT import TNamed
22 from ROOT.Math import XYZVector
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().__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 = 'Leo Piilonen (piilonen@vt.edu)'
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, -15., 15.)
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 = XYZVector(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 = XYZVector(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  four_momentum_k = klm_cluster.getMomentum()
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(four_momentum_k.P() - momentum.R())
220  self.hist_ptreshist_ptres.Fill(four_momentum_k.Theta() - momentum.Theta())
221  self.hist_ppreshist_ppres.Fill(four_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:72
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