Belle II Software development
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
18import basf2
19import ROOT
20from ROOT import Belle2
21from ROOT import TNamed
22from ROOT.Math import XYZVector
23import math
24import numpy
25
26
27class 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.evtgen = evtgen
35
36 self.check_eklm = check_eklm
37
38 self.output_file = ROOT.TFile(output_file, 'recreate')
39 contact = 'Leo Piilonen (piilonen@vt.edu)'
40
41 self.hist_nkl = ROOT.TH1F('k0l_number',
42 'Number of KLM clusters per 1 MC particle',
43 5, -0.5, 4.5)
44 self.hist_nkl.SetXTitle('KLM clusters')
45 self.hist_nkl.SetYTitle('Events')
46 functions = self.hist_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_xres = ROOT.TH1F('k0l_xres',
54 'KLM K0L decay vertex X resolution',
55 100, -50, 50)
56 self.hist_xres.SetXTitle('cm')
57 self.hist_xres.SetYTitle('Events')
58 functions = self.hist_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_yres = ROOT.TH1F('k0l_yres',
65 'KLM K0L decay vertex Y resolution',
66 100, -50, 50)
67 self.hist_yres.SetXTitle('cm')
68 self.hist_yres.SetYTitle('Events')
69 functions = self.hist_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_zres = ROOT.TH1F('k0l_zres',
76 'KLM K0L decay vertex Z resolution',
77 100, -50, 50)
78 self.hist_zres.SetXTitle('cm')
79 self.hist_zres.SetYTitle('Events')
80 functions = self.hist_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_tres = ROOT.TH1F('k0l_tres',
87 'KLM K0L decay time resolution',
88 100, -15., 15.)
89 self.hist_tres.SetXTitle('ns')
90 self.hist_tres.SetYTitle('Events')
91 functions = self.hist_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_pres = ROOT.TH1F('k0l_pres',
98 'KLM K0L momentum resolution',
99 100, -3., 3.)
100 self.hist_pres.SetXTitle('GeV')
101 self.hist_pres.SetYTitle('Events')
102 functions = self.hist_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_ptres = ROOT.TH1F('k0l_ptres',
109 'KLM K0L momentum theta resolution',
110 100, -0.2, 0.2)
111 self.hist_ptres.SetXTitle('rad')
112 self.hist_ptres.SetYTitle('Events')
113 functions = self.hist_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_ppres = ROOT.TH1F('k0l_ppres',
120 'KLM K0L momentum phi resolution',
121 100, -0.2, 0.2)
122 self.hist_ppres.SetXTitle('rad')
123 self.hist_ppres.SetYTitle('Events')
124 functions = self.hist_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_covmat = ROOT.TH1F('k0l_covmat',
131 'KLM K0L coordinates covariance matrix',
132 6, 0, 1)
133 self.hist_covmat.GetXaxis().SetBinLabel(1, 'xx')
134 self.hist_covmat.GetXaxis().SetBinLabel(2, 'xy')
135 self.hist_covmat.GetXaxis().SetBinLabel(3, 'xz')
136 self.hist_covmat.GetXaxis().SetBinLabel(4, 'yy')
137 self.hist_covmat.GetXaxis().SetBinLabel(5, 'yz')
138 self.hist_covmat.GetXaxis().SetBinLabel(6, 'zz')
139 self.hist_covmat.SetYTitle('Covariance, cm^{2}')
140 functions = self.hist_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_corrmat = ROOT.TH1F('k0l_corrmat',
147 'KLM K0L correlation matrix',
148 10, 0, 1)
149 self.hist_corrmat.GetXaxis().SetBinLabel(1, 'xx')
150 self.hist_corrmat.GetXaxis().SetBinLabel(2, 'xy')
151 self.hist_corrmat.GetXaxis().SetBinLabel(3, 'xz')
152 self.hist_corrmat.GetXaxis().SetBinLabel(4, 'xp')
153 self.hist_corrmat.GetXaxis().SetBinLabel(5, 'yy')
154 self.hist_corrmat.GetXaxis().SetBinLabel(6, 'yz')
155 self.hist_corrmat.GetXaxis().SetBinLabel(7, 'yp')
156 self.hist_corrmat.GetXaxis().SetBinLabel(8, 'zz')
157 self.hist_corrmat.GetXaxis().SetBinLabel(9, 'zp')
158 self.hist_corrmat.GetXaxis().SetBinLabel(10, 'pp')
159 self.hist_corrmat.SetYTitle('Correlation coefficient')
160 functions = self.hist_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_av = XYZVector(0, 0, 0)
167
168 self.vertex = []
169
170 self.momentum_av = 0
171
172 self.momentum = []
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.evtgen):
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_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_nkl.Fill(len(klm_clusters))
206 for klm_cluster in klm_clusters:
207 vertex_k = XYZVector(klm_cluster.getClusterPosition()) - vertex
208 self.vertex.append(vertex_k)
209 self.vertex_k_av = self.vertex_k_av + vertex_k
210 self.momentum.append(klm_cluster.getMomentumMag())
211 self.momentum_av = self.momentum_av + \
212 klm_cluster.getMomentumMag()
213 time_k = klm_cluster.getTime()
214 four_momentum_k = klm_cluster.getMomentum()
215 self.hist_xres.Fill(vertex_k.X())
216 self.hist_yres.Fill(vertex_k.Y())
217 self.hist_zres.Fill(vertex_k.Z())
218 self.hist_tres.Fill(time_k - time)
219 self.hist_pres.Fill(four_momentum_k.P() - momentum.R())
220 self.hist_ptres.Fill(four_momentum_k.Theta() - momentum.Theta())
221 self.hist_ppres.Fill(four_momentum_k.Phi() - momentum.Phi())
222
223 def terminate(self):
224 """ Termination function. """
225 self.vertex_k_av = self.vertex_k_av * (1.0 / len(self.vertex))
226 self.momentum_av = self.momentum_av * (1.0 / len(self.vertex))
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.vertex)):
233 cov_mat[0][0] = cov_mat[0][0] + \
234 (self.vertex[i].X() - self.vertex_k_av.X()) * \
235 (self.vertex[i].X() - self.vertex_k_av.X())
236 cov_mat[0][1] = cov_mat[0][1] + \
237 (self.vertex[i].X() - self.vertex_k_av.X()) * \
238 (self.vertex[i].Y() - self.vertex_k_av.Y())
239 cov_mat[0][2] = cov_mat[0][2] + \
240 (self.vertex[i].X() - self.vertex_k_av.X()) * \
241 (self.vertex[i].Z() - self.vertex_k_av.Z())
242 cov_mat[0][3] = cov_mat[0][3] + \
243 (self.vertex[i].X() - self.vertex_k_av.X()) * \
244 (self.momentum[i] - self.momentum_av)
245 cov_mat[1][1] = cov_mat[1][1] + \
246 (self.vertex[i].Y() - self.vertex_k_av.Y()) * \
247 (self.vertex[i].Y() - self.vertex_k_av.Y())
248 cov_mat[1][2] = cov_mat[1][2] + \
249 (self.vertex[i].Y() - self.vertex_k_av.Y()) * \
250 (self.vertex[i].Z() - self.vertex_k_av.Z())
251 cov_mat[1][3] = cov_mat[1][3] + \
252 (self.vertex[i].Y() - self.vertex_k_av.Y()) * \
253 (self.momentum[i] - self.momentum_av)
254 cov_mat[2][2] = cov_mat[2][2] + \
255 (self.vertex[i].Z() - self.vertex_k_av.Z()) * \
256 (self.vertex[i].Z() - self.vertex_k_av.Z())
257 cov_mat[2][3] = cov_mat[2][3] + \
258 (self.vertex[i].Z() - self.vertex_k_av.Z()) * \
259 (self.momentum[i] - self.momentum_av)
260 cov_mat[3][3] = cov_mat[3][3] + \
261 (self.momentum[i] - self.momentum_av) * \
262 (self.momentum[i] - self.momentum_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.vertex) - 1)
266 for i in range(len(self.vertex)):
267 cov_mat_err[0][0] = cov_mat_err[0][0] + \
268 pow((self.vertex[i].X() - self.vertex_k_av.X()) *
269 (self.vertex[i].X() - self.vertex_k_av.X()) - cov_mat[0][0], 2)
270 cov_mat_err[0][1] = cov_mat_err[0][1] + \
271 pow((self.vertex[i].X() - self.vertex_k_av.X()) *
272 (self.vertex[i].Y() - self.vertex_k_av.Y()) - cov_mat[0][1], 2)
273 cov_mat_err[0][2] = cov_mat_err[0][2] + \
274 pow((self.vertex[i].X() - self.vertex_k_av.X()) *
275 (self.vertex[i].Z() - self.vertex_k_av.Z()) - cov_mat[0][2], 2)
276 cov_mat_err[0][3] = cov_mat_err[0][3] + \
277 pow((self.vertex[i].X() - self.vertex_k_av.X()) *
278 (self.momentum[i] - self.momentum_av) - cov_mat[0][3], 2)
279 cov_mat_err[1][1] = cov_mat_err[1][1] + \
280 pow((self.vertex[i].Y() - self.vertex_k_av.Y()) *
281 (self.vertex[i].Y() - self.vertex_k_av.Y()) - cov_mat[1][1], 2)
282 cov_mat_err[1][2] = cov_mat_err[1][2] + \
283 pow((self.vertex[i].Y() - self.vertex_k_av.Y()) *
284 (self.vertex[i].Z() - self.vertex_k_av.Z()) - cov_mat[1][2], 2)
285 cov_mat_err[1][3] = cov_mat_err[1][3] + \
286 pow((self.vertex[i].Y() - self.vertex_k_av.Y()) *
287 (self.momentum[i] - self.momentum_av) - cov_mat[1][3], 2)
288 cov_mat_err[2][2] = cov_mat_err[2][2] + \
289 pow((self.vertex[i].Z() - self.vertex_k_av.Z()) *
290 (self.vertex[i].Z() - self.vertex_k_av.Z()) - cov_mat[2][2], 2)
291 cov_mat_err[2][3] = cov_mat_err[2][3] + \
292 pow((self.vertex[i].Z() - self.vertex_k_av.Z()) *
293 (self.momentum[i] - self.momentum_av) - cov_mat[2][3], 2)
294 cov_mat_err[3][3] = cov_mat_err[3][3] + \
295 pow((self.momentum[i] - self.momentum_av) *
296 (self.momentum[i] - self.momentum_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.vertex) - 2) * len(self.vertex)))
302 self.hist_covmat.SetBinContent(1, cov_mat[0][0])
303 self.hist_covmat.SetBinContent(2, cov_mat[0][1])
304 self.hist_covmat.SetBinContent(3, cov_mat[0][2])
305 self.hist_covmat.SetBinContent(4, cov_mat[1][1])
306 self.hist_covmat.SetBinContent(5, cov_mat[1][2])
307 self.hist_covmat.SetBinContent(6, cov_mat[2][2])
308 self.hist_covmat.SetBinError(1, cov_mat_err[0][0])
309 self.hist_covmat.SetBinError(2, cov_mat_err[0][1])
310 self.hist_covmat.SetBinError(3, cov_mat_err[0][2])
311 self.hist_covmat.SetBinError(4, cov_mat_err[1][1])
312 self.hist_covmat.SetBinError(5, cov_mat_err[1][2])
313 self.hist_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_corrmat.SetBinContent(1, corr_mat[0][0])
322 self.hist_corrmat.SetBinContent(2, corr_mat[0][1])
323 self.hist_corrmat.SetBinContent(3, corr_mat[0][2])
324 self.hist_corrmat.SetBinContent(4, corr_mat[0][3])
325 self.hist_corrmat.SetBinContent(5, corr_mat[1][1])
326 self.hist_corrmat.SetBinContent(6, corr_mat[1][2])
327 self.hist_corrmat.SetBinContent(7, corr_mat[1][3])
328 self.hist_corrmat.SetBinContent(8, corr_mat[2][2])
329 self.hist_corrmat.SetBinContent(9, corr_mat[2][3])
330 self.hist_corrmat.SetBinContent(10, corr_mat[3][3])
331 self.hist_corrmat.SetBinError(1, corr_mat_err[0][0])
332 self.hist_corrmat.SetBinError(2, corr_mat_err[0][1])
333 self.hist_corrmat.SetBinError(3, corr_mat_err[0][2])
334 self.hist_corrmat.SetBinError(4, corr_mat_err[0][3])
335 self.hist_corrmat.SetBinError(5, corr_mat_err[1][1])
336 self.hist_corrmat.SetBinError(6, corr_mat_err[1][2])
337 self.hist_corrmat.SetBinError(7, corr_mat_err[1][3])
338 self.hist_corrmat.SetBinError(8, corr_mat_err[2][2])
339 self.hist_corrmat.SetBinError(9, corr_mat_err[2][3])
340 self.hist_corrmat.SetBinError(10, corr_mat_err[3][3])
341 self.output_file.cd()
342 self.hist_nkl.Write()
343 self.hist_xres.Write()
344 self.hist_yres.Write()
345 self.hist_zres.Write()
346 self.hist_tres.Write()
347 self.hist_pres.Write()
348 self.hist_ptres.Write()
349 self.hist_ppres.Write()
350 self.hist_covmat.Write()
351 self.hist_corrmat.Write()
352 self.output_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