Belle II Software development
plotVXDAlignmentPayload.py
1#!/usr/bin/env python3
2
3
10
11import numpy as np
12import pandas as pd
13import math
14from pylab import savefig, subplot
15import matplotlib.pyplot as plt
16from matplotlib import ticker
17import matplotlib as mpl
18import sys
19
20
21import ROOT
22from ROOT import Belle2
23if len(sys.argv) < 2:
24 sys.exit("No input .root file specified!")
25
26inputroot = sys.argv[1]
27file = ROOT.TFile(inputroot, "OPEN")
28vxd = file.Get("VXDAlignment")
29
30fileName = inputroot + '.txt'
31text_file = open(fileName, "w")
32text_file.write("layer ladder sensor param value\n")
33for entry in vxd.getMap():
34 element_parameter = entry.first
35 value = entry.second
36 element = element_parameter.first
37 param = element_parameter.second
38 vxdid = Belle2.VxdID(element)
39 layer = vxdid.getLayerNumber()
40 ladder = vxdid.getLadderNumber()
41 sensor = vxdid.getSensorNumber()
42 if sensor != 0:
43 text_file.write("{} {} {} {} {}\n".format(layer, ladder, sensor, param, value))
44text_file.close()
45
46# set to True, if automatic scaling is to be used, otherwise the values in the unit_scale dict will be used.
47autoscale = True
48
49# datafile structure
50components = {'value'}
51comp_map = {'value': 0}
52parameters = {'du', 'dv', 'dw', 'alpha', 'beta', 'gamma', 'P_20', 'P_11', 'P_02', 'P_30', 'P_21', 'P_12', 'P_03'}
53param_map = {1: 'du', 2: 'dv', 3: 'dw', 4: 'alpha', 5: 'beta', 6: 'gamma',
54 31: 'P_20', 32: 'P_11', 33: 'P_02',
55 41: 'P_30', 42: 'P_21', 43: 'P_12', 44: 'P_03'}
56param_nums = {u: v for (v, u) in param_map.items()}
57layers = {1, 2, 3, 4, 5, 6}
58
59n_data = 10 # number of independent samples
60
61# OK, we could make a dataframe for this, but then - why?
62layer_phi0 = {
63 1: 0.0,
64 2: 0.0,
65 3: 33.624 * np.pi / 180,
66 4: 8 * np.pi / 180,
67 5: -8 * np.pi / 180,
68 6: -4 * np.pi / 180,
69}
70layer_nladders = {
71 1: 8,
72 2: 12,
73 3: 7,
74 4: 10,
75 5: 12,
76 6: 16,
77}
78layer_nsensors = {
79 1: 2,
80 2: 2,
81 3: 2,
82 4: 3,
83 5: 4,
84 6: 5,
85}
86
87# If you want to use preset scales, enter them here and set autoscale to False.
88unit_scale = {'value': {'shift': (-10.0,
89 10.,
90 plt.cm.bwr,
91 ' [um]'),
92 'rotation': (-1.0,
93 1.0,
94 plt.cm.bwr,
95 ' [mrad]'),
96 'surface': (-10.0,
97 10.,
98 plt.cm.bwr,
99 ' [um]')}}
100
101# Read data into a pandas DataFrame
102print('reading from file ' + fileName)
103
104dataframe = pd.read_table(
105 fileName,
106 sep=' ',
107 skipinitialspace=True,
108 skiprows=1,
109 names=[
110 'layer',
111 'ladder',
112 'sensor',
113 'parameter',
114 'value'])
115
116# Scale things properly
117scale = np.where((dataframe.parameter < 4) | (dataframe.parameter > 6), 1.0e4, 1.0e3)
118for colname in ['value']:
119 dataframe[colname] *= scale
120
121# If bias is average of 10 samples, scale t_mille by sqrt(10)
122# dataframe['t_mille'] *= np.sqrt(10)
123
124# Calculate and print scales.
125# Maxima can be calculated for series of columns
126data_columns = ['value']
127min_shift = dataframe[data_columns][dataframe.parameter < 4].min().min()
128max_shift = dataframe[data_columns][dataframe.parameter < 4].max().max()
129min_shift = min(min_shift, -max_shift)
130max_shift = max(-min_shift, max_shift)
131min_rot = dataframe[data_columns][(3 < dataframe.parameter) & (dataframe.parameter < 7)].min().min()
132max_rot = dataframe[data_columns][(3 < dataframe.parameter) & (dataframe.parameter < 7)].max().max()
133min_rot = min(min_rot, -max_rot)
134max_rot = max(-min_rot, max_rot)
135min_sur2 = dataframe[data_columns][(30 < dataframe.parameter) & (dataframe.parameter < 34)].min().min()
136max_sur2 = dataframe[data_columns][(30 < dataframe.parameter) & (dataframe.parameter < 34)].max().max()
137min_sur2 = min(min_sur2, -max_sur2)
138max_sur2 = max(-min_sur2, max_sur2)
139min_sur3 = dataframe[data_columns][(40 < dataframe.parameter) & (dataframe.parameter < 45)].min().min()
140max_sur3 = dataframe[data_columns][(40 < dataframe.parameter) & (dataframe.parameter < 45)].max().max()
141min_sur3 = min(min_sur3, -max_sur3)
142max_sur3 = max(-min_sur3, max_sur3)
143
144if math.isnan(min_rot):
145 min_rot = 0.
146if math.isnan(max_rot):
147 max_rot = 0.
148if math.isnan(min_shift):
149 min_shift = 0.
150if math.isnan(max_shift):
151 max_shift = 0.
152if math.isnan(min_sur2):
153 min_sur2 = 0.
154if math.isnan(max_sur2):
155 max_sur2 = 0.
156if math.isnan(min_sur3):
157 min_sur3 = 0.
158if math.isnan(max_sur3):
159 max_sur3 = 0.
160
161print('Symmetrized scales for value: ')
162print('Shifts: ' + str(min_shift) + ' to ' + str(max_shift))
163print('Rotations: ' + str(min_rot) + ' to ' + str(max_rot))
164print('Surface^2 ' + str(min_sur2) + ' to ' + str(max_sur2))
165print('Surface^3 ' + str(min_sur3) + ' to ' + str(max_sur3))
166
167if autoscale:
168 unit_scale['value']['shift'] = (min_shift, max_shift, plt.cm.bwr, ' [um]')
169 unit_scale['value']['rotation'] = (min_rot, max_rot, plt.cm.bwr, ' [mrad]')
170 unit_scale['value']['surface2'] = (min_sur2, max_sur2, plt.cm.bwr, ' [um]')
171 unit_scale['value']['surface3'] = (min_sur3, max_sur3, plt.cm.bwr, ' [um]')
172
173fig = plt.figure(figsize=(10, 6.5))
174for value in components:
175 for shiftrot, i_params in zip(['shift', 'rotation'], [[1, 2, 3], [4, 5, 6]]):
176 cmin, cmax, cmap, unit = unit_scale[value][shiftrot]
177 for i_param in i_params:
178 parameter = param_map[i_param]
179 # Make a subframe
180 val_param = dataframe[['layer', 'ladder', 'sensor', value]][dataframe['parameter'] == i_param]
181 val_param[value] = val_param[value]
182 # Here comes the real thing
183 ax = subplot(2, 3, i_param, projection='polar')
184 layer_gap = 24
185 bottom = 24 # starts at this, increases by sensor height in layers and layer_gap in between.
186 for layer in layers:
187 ladder_width = 2 * np.pi / layer_nladders[layer]
188 ladder_angles = np.linspace(layer_phi0[layer], layer_phi0[layer] + 2 * np.pi, layer_nladders[layer], endpoint=False)
189
190 for sensor in range(1, 1 + layer_nsensors[layer]):
191 height = (layer > 3 and sensor == 1) * 10.0 + 28.0 / (1.0 + 0.3 * layer + 0.5 * sensor)
192 content = height * np.ones(layer_nladders[layer])
193 bars = ax.bar(
194 ladder_angles,
195 content,
196 width=ladder_width,
197 bottom=bottom,
198 label=str(sensor + 1),
199 align='center',
200 # edgecolor='k',
201 )
202 bottom += height
203 for (ladder, bar) in zip(range(1, 1 + layer_nladders[layer]), bars):
204 cond = (val_param.layer == layer) & (val_param.ladder == ladder) & (val_param.sensor == sensor)
205 try:
206 c = float(val_param[value][cond])
207 except BaseException:
208 c = 0.
209 if (cmax - cmin) > 0.:
210 color_value = (c - cmin) / (cmax - cmin)
211 else:
212 color_value = 0.
213 bar.set_facecolor(cmap(color_value))
214 bar.set_linewidth(0.2)
215
216 label_r = bottom + 8
217 for i in range(layer_nladders[layer]):
218 ax.text(
219 ladder_angles[i],
220 label_r,
221 str(i + 1),
222 ha='center',
223 va='center',
224 fontsize=6,
225 # color='k',
226 )
227 bottom += layer_gap * (1 + 0.5 * (layer == 2))
228
229 ax.grid(False)
230 ax.axis('off')
231
232 plt.title(parameter)
233
234 # if we are done with shifts/rotations, collect some data
235 # for the colorbar
236 if shiftrot == 'shift':
237 par_shifts = (cmap, cmin, cmax, unit)
238 elif shiftrot == 'rotation':
239 par_rots = (cmap, cmin, cmax, unit)
240
241 # Plot colorbars after rings are in place
242 axb1 = fig.add_axes([0.01, 0.55, 0.02, 0.33])
243 cmap, cmin, cmax, unit = par_shifts
244 norm1 = mpl.colors.Normalize(vmin=cmin, vmax=cmax)
245 cb1 = mpl.colorbar.ColorbarBase(axb1, cmap=cmap, norm=norm1, orientation='vertical')
246 cb1.set_label('Shifts ' + unit)
247 cb1.locator = ticker.MaxNLocator(nbins=11)
248 cb1.update_ticks()
249 # rotations
250 axb2 = fig.add_axes([0.01, 0.11, 0.02, 0.33])
251 cmap, cmin, cmax, unit = par_rots
252 norm2 = mpl.colors.Normalize(vmin=cmin, vmax=cmax)
253 cb2 = mpl.colorbar.ColorbarBase(axb2, cmap=cmap, norm=norm2, orientation='vertical')
254 cb2.set_label('Rotations ' + unit)
255 cb2.locator = ticker.MaxNLocator(nbins=11)
256 cb2.update_ticks()
257
258 fig.suptitle(value, fontsize=20, x=0.01, y=0.95, horizontalalignment='left')
259 figname = inputroot + '.rigid.png'
260 savefig(figname)
261 print('Saved figure ' + figname)
262
263fig = plt.figure(figsize=(14, 6.5))
264for value in components:
265 for shiftrot, i_params in zip(['surface2', 'surface3'], [[31, 32, 33], [41, 42, 43, 44]]):
266 cmin, cmax, cmap, unit = unit_scale[value][shiftrot]
267 for i_param in i_params:
268 parameter = param_map[i_param]
269 # Make a subframe
270 val_param = dataframe[['layer', 'ladder', 'sensor', value]][dataframe['parameter'] == i_param]
271 val_param[value] = val_param[value]
272 # Here comes the real thing
273 if shiftrot == 'surface2':
274 ax = subplot(2, 4, i_param-30, projection='polar')
275 if shiftrot == 'surface3':
276 ax = subplot(2, 4, i_param-36, projection='polar')
277 layer_gap = 24
278 bottom = 24 # starts at this, increases by sensor height in layers and layer_gap in between.
279 for layer in layers:
280 ladder_width = 2 * np.pi / layer_nladders[layer]
281 ladder_angles = np.linspace(layer_phi0[layer], layer_phi0[layer] + 2 * np.pi, layer_nladders[layer], endpoint=False)
282
283 for sensor in range(1, 1 + layer_nsensors[layer]):
284 height = (layer > 3 and sensor == 1) * 10.0 + 28.0 / (1.0 + 0.3 * layer + 0.5 * sensor)
285 content = height * np.ones(layer_nladders[layer])
286 bars = ax.bar(
287 ladder_angles,
288 content,
289 width=ladder_width,
290 bottom=bottom,
291 label=str(sensor + 1),
292 align='center',
293 # edgecolor='k',
294 )
295 bottom += height
296 for (ladder, bar) in zip(range(1, 1 + layer_nladders[layer]), bars):
297 cond = (val_param.layer == layer) & (val_param.ladder == ladder) & (val_param.sensor == sensor)
298 try:
299 c = float(val_param[value][cond])
300 except BaseException:
301 c = 0.
302 if (cmax - cmin) > 0.:
303 color_value = (c - cmin) / (cmax - cmin)
304 else:
305 color_value = 0.
306 bar.set_facecolor(cmap(color_value))
307 bar.set_linewidth(0.2)
308
309 label_r = bottom + 8
310 for i in range(layer_nladders[layer]):
311 ax.text(
312 ladder_angles[i],
313 label_r,
314 str(i + 1),
315 ha='center',
316 va='center',
317 fontsize=6,
318 color='k',
319 )
320 bottom += layer_gap * (1 + 0.5 * (layer == 2))
321
322 ax.grid(False)
323 ax.axis('off')
324
325 plt.title(parameter)
326
327 # if we are done with shifts/rotations, collect some data
328 # for the colorbar
329 if shiftrot == 'surface2':
330 par_sur2 = (cmap, cmin, cmax, unit)
331 elif shiftrot == 'surface3':
332 par_sur3 = (cmap, cmin, cmax, unit)
333
334 # Plot colorbars after rings are in place
335 axb1 = fig.add_axes([0.01, 0.55, 0.02, 0.33])
336 cmap, cmin, cmax, unit = par_sur2
337 norm1 = mpl.colors.Normalize(vmin=cmin, vmax=cmax)
338 cb1 = mpl.colorbar.ColorbarBase(axb1, cmap=cmap, norm=norm1, orientation='vertical')
339 cb1.set_label('Quadratic parameters ' + unit)
340 cb1.locator = ticker.MaxNLocator(nbins=11)
341 cb1.update_ticks()
342 # rotations
343 axb2 = fig.add_axes([0.01, 0.11, 0.02, 0.33])
344 cmap, cmin, cmax, unit = par_sur3
345 norm2 = mpl.colors.Normalize(vmin=cmin, vmax=cmax)
346 cb2 = mpl.colorbar.ColorbarBase(axb2, cmap=cmap, norm=norm2, orientation='vertical')
347 cb2.set_label('Cubic parameters ' + unit)
348 cb2.locator = ticker.MaxNLocator(nbins=11)
349 cb2.update_ticks()
350
351 fig.suptitle(value, fontsize=20, x=0.01, y=0.95, horizontalalignment='left')
352
353 figname = inputroot + '.surface.png'
354 savefig(figname)
355 print('Saved figure ' + figname)
Belle2::VxdID
Class to uniquely identify a any structure of the PXD and SVD.
Definition: VxdID.h:33