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