Belle II Software  release-08-01-10
VXDOverlaps.py
1 
8 
9 import os
10 import ROOT
11 import numpy
12 from array import array
13 from ROOT import TCanvas, TH1F, TLine
14 from hist_utils import hist2array
15 
16 
17 '''Takes the output of combined OverlapResiduals and HistoManager modules,
18  as input, providess hit-maps for overlapping VXD hits in Layer:Sensor
19  plots and computes statistics for 2D monitoring plots.'''
20 
21 
22 def Median_plots_phi(filename, lyr_num, phi_bins, phi_inf, phi_sup):
23  """
24  Function to compute the median of the projected DeltaResU(V)
25  distributions for each azimuthal overlap in DeltaResU(V) vs phi plots
26  """
27  # Accesses the TDirectory containing VXD histograms of residuals differences for overlapping hits
28  f = ROOT.TFile.Open(filename, 'read')
29  mn = f.Get('Monitoring_VXDOverlaps')
30  # Gets 2D DeltaRes_u vs phi histograms stored in histofile.root
31  h_PhiU = mn.Get('h_DeltaResUPhi_Lyr' + str(lyr_num))
32  # Gets 2D DeltaRes_v vs phi histograms stored in histofile.root
33  h_PhiV = mn.Get('h_DeltaResVPhi_Lyr' + str(lyr_num))
34  # Defines histograms of projected DeltaRes_u distributions
35  h_UMedians = TH1F(
36  'h_UMedians_Lyr' +
37  str(lyr_num),
38  'Layer' +
39  str(lyr_num) +
40  ': medians of #Deltares_{u} for each overlap',
41  phi_bins,
42  phi_inf,
43  phi_sup)
44  # Defines histograms of projected DeltaRes_u distributions
45  h_VMedians = TH1F(
46  'h_VMedians_Lyr' +
47  str(lyr_num),
48  'Layer' +
49  str(lyr_num) +
50  ': medians of #Deltares_{v} for each overlap',
51  phi_bins,
52  phi_inf,
53  phi_sup)
54  # Sets the number of resampled samples for bootstrap
55  Nrs = 100
56  # Initializes parameters to compute the median of ROOT TH1
57  q_U = array('d', [0])
58  p_U = array('d', [0.5])
59  q_V = array('d', [0])
60  p_V = array('d', [0.5])
61  # Lists storing medians of toy MC distributions for statistical bootstrap
62  l_U_median = []
63  l_V_median = []
64  # Lists storing estimated medians of projected DeltaRes distribution for every azimuthal overlap
65  l_U_median_pos = []
66  l_V_median_pos = []
67  # Computes medians and related uncertainties (statistical bootstrapping) and produces plots
68  c_PhiU = TCanvas('c_PhiU_' + str(lyr_num), 'DeltaResUPhi_' + str(lyr_num), 700, 500)
69  c_PhiV = TCanvas('c_PhiV_' + str(lyr_num), 'DeltaResVPhi_' + str(lyr_num), 700, 500)
70  if(lyr_num == 1 or lyr_num == 3):
71  c_PhiU.Divide(4, 2)
72  c_PhiV.Divide(4, 2)
73  if(lyr_num == 4):
74  c_PhiU.Divide(5, 2)
75  c_PhiV.Divide(5, 2)
76  if(lyr_num == 5):
77  c_PhiU.Divide(5, 3)
78  c_PhiV.Divide(5, 3)
79  if(lyr_num == 6):
80  c_PhiU.Divide(5, 4)
81  c_PhiV.Divide(5, 4)
82  c_PhiUMedians = TCanvas('c_PhiUMedians_' + str(lyr_num), 'UMedians_' + str(lyr_num), 700, 500)
83  c_PhiVMedians = TCanvas('c_PhiVMedians_' + str(lyr_num), 'VMedians_' + str(lyr_num), 700, 500)
84  for i in range(0, phi_bins):
85  xinf = phi_inf + i * (phi_sup - phi_inf) / phi_bins
86  xsup = phi_inf + (i + 1) * (phi_sup - phi_inf) / phi_bins
87  h_PhiU.GetXaxis().SetRangeUser(xinf, xsup)
88  h_PhiV.GetXaxis().SetRangeUser(xinf, xsup)
89  h_PhiU.ProjectionY().GetQuantiles(1, q_U, p_U)
90  h_PhiV.ProjectionY().GetQuantiles(1, q_V, p_V)
91  h_U = h_PhiU.ProjectionY()
92  h_V = h_PhiV.ProjectionY()
93  h_U.SetTitle(str(round(xinf, 3)) + ' < #phi < ' + str(round(xsup, 3)))
94  h_V.SetTitle(str(round(xinf, 3)) + ' < #phi < ' + str(round(xsup, 3)))
95  h_U.GetXaxis().SetRangeUser(-200, 200)
96  h_V.GetXaxis().SetRangeUser(-200, 200)
97  if(lyr_num == 1):
98  h_U.GetYaxis().SetRangeUser(0, 300)
99  h_V.GetYaxis().SetRangeUser(0, 300)
100  median_pos_U = TLine(q_U[0], 0, q_U[0], 300)
101  median_pos_V = TLine(q_V[0], 0, q_V[0], 300)
102  if(lyr_num == 3):
103  h_U.GetYaxis().SetRangeUser(0, 500)
104  h_V.GetYaxis().SetRangeUser(0, 500)
105  median_pos_U = TLine(q_U[0], 0, q_U[0], 500)
106  median_pos_V = TLine(q_V[0], 0, q_V[0], 500)
107  if(lyr_num == 4):
108  h_U.GetYaxis().SetRangeUser(0, 2000)
109  h_V.GetYaxis().SetRangeUser(0, 2000)
110  median_pos_U = TLine(q_U[0], 0, q_U[0], 2000)
111  median_pos_V = TLine(q_V[0], 0, q_V[0], 2000)
112  if(lyr_num == 5):
113  h_U.GetYaxis().SetRangeUser(0, 1000)
114  h_V.GetYaxis().SetRangeUser(0, 1000)
115  median_pos_U = TLine(q_U[0], 0, q_U[0], 1000)
116  median_pos_V = TLine(q_V[0], 0, q_V[0], 1000)
117  if(lyr_num == 6):
118  h_U.GetYaxis().SetRangeUser(0, 1000)
119  h_V.GetYaxis().SetRangeUser(0, 1000)
120  median_pos_U = TLine(q_U[0], 0, q_U[0], 1000)
121  median_pos_V = TLine(q_V[0], 0, q_V[0], 1000)
122  median_pos_U.SetLineWidth(2)
123  median_pos_V.SetLineWidth(2)
124  median_pos_U.SetLineColor(2)
125  median_pos_V.SetLineColor(2)
126  l_U_median_pos.append(median_pos_U)
127  l_V_median_pos.append(median_pos_V)
128  h_U.GetYaxis().SetTitle('counts')
129  h_V.GetYaxis().SetTitle('counts')
130  meas_U = hist2array(h_U)
131  meas_V = hist2array(h_V)
132  bs_U = numpy.random.poisson(1., (len(meas_U), Nrs))
133  bs_V = numpy.random.poisson(1., (len(meas_V), Nrs))
134  for j in range(Nrs):
135  toy_U = numpy.repeat(meas_U, bs_U[:, j])
136  toy_V = numpy.repeat(meas_V, bs_V[:, j])
137  median_U_toy = numpy.median(toy_U)
138  median_V_toy = numpy.median(toy_V)
139  l_U_median.append(median_U_toy)
140  l_V_median.append(median_V_toy)
141  median_U_rs = numpy.array(l_U_median)
142  median_V_rs = numpy.array(l_V_median)
143  median_U_dev = numpy.std(median_U_rs)
144  median_V_dev = numpy.std(median_V_rs)
145  h_UMedians.SetBinContent(i + 1, q_U[0])
146  h_UMedians.SetBinError(i + 1, median_U_dev)
147  h_VMedians.SetBinContent(i + 1, q_V[0])
148  h_VMedians.SetBinError(i + 1, median_V_dev)
149  c_PhiU.cd(i + 1)
150  h_U.DrawCopy()
151  l_U_median_pos[i].Draw("SAME")
152  c_PhiV.cd(i + 1)
153  h_V.DrawCopy()
154  l_V_median_pos[i].Draw("SAME")
155  c_PhiUMedians.cd()
156  h_UMedians.GetXaxis().SetTitle('#phi (rad)')
157  h_UMedians.GetYaxis().SetTitle('Median of #Deltares_{u} (#mum)')
158  h_UMedians.Draw()
159  c_PhiVMedians.cd()
160  h_VMedians.GetXaxis().SetTitle('#phi (rad)')
161  h_VMedians.GetYaxis().SetTitle('Median of #Deltares_{V} (#mum)')
162  h_VMedians.Draw()
163  # If not existing, creates a dedicated folder
164  if not os.path.exists('Median_plots_OverlapsPhi'):
165  os.mkdir('Median_plots_OverlapsPhi')
166  c_PhiU.SaveAs('Median_plots_OverlapsPhi/Median_and_DeltaResUPhi_Lyr' + str(lyr_num) + '.root')
167  c_PhiU.SaveAs('Median_plots_OverlapsPhi/Median_and_DeltaResUPhi_Lyr' + str(lyr_num) + '.pdf')
168  c_PhiV.SaveAs('Median_plots_OverlapsPhi/Median_and_DeltaResVPhi_Lyr' + str(lyr_num) + '.root')
169  c_PhiV.SaveAs('Median_plots_OverlapsPhi/Median_and_DeltaResVPhi_Lyr' + str(lyr_num) + '.pdf')
170  c_PhiUMedians.SaveAs('Median_plots_OverlapsPhi/Lyr' + str(lyr_num) + '_DeltaResUMedians_vs_phi.root')
171  c_PhiUMedians.SaveAs('Median_plots_OverlapsPhi/Lyr' + str(lyr_num) + '_DeltaResUMedians_vs_phi.pdf')
172  c_PhiVMedians.SaveAs('Median_plots_OverlapsPhi/Lyr' + str(lyr_num) + '_DeltaResVMedians_vs_phi.root')
173  c_PhiVMedians.SaveAs('Median_plots_OverlapsPhi/Lyr' + str(lyr_num) + '_DeltaResVMedians_vs_phi.pdf')
174  return
175 
176 
177 def Median_plots_z(filename, lyr_num, z_bins, z_inf, z_sup):
178  """
179  Function to compute the median of the projected DeltaResU(V) distributions
180  for each sensor in DeltaResU(V) vs z plots
181  """
182  # Accesses the TDirectory containing VXD histograms of residuals differences for overlapping hits
183  f = ROOT.TFile.Open(filename, 'read')
184  mn = f.Get('Monitoring_VXDOverlaps')
185  # Gets 2D DeltaRes_u vs z histograms stored in histofile.root
186  h_ZU = mn.Get('h_DeltaResUz_Lyr' + str(lyr_num))
187  # Gets 2D DeltaRes_v vs z histograms stored in histofile.root
188  h_ZV = mn.Get('h_DeltaResVz_Lyr' + str(lyr_num))
189  # Defines histograms of projected DeltaRes_u distributions
190  h_UMedians = TH1F('h_UMedians_Lyr' + str(lyr_num), 'Layer' + str(lyr_num) +
191  ': medians of #Deltares_{u} for each sensor', z_bins, z_inf, z_sup)
192  # Defines histograms of projected DeltaRes_u distributions
193  h_VMedians = TH1F('h_VMedians_Lyr' + str(lyr_num), 'Layer' + str(lyr_num) +
194  ': medians of #Deltares_{v} for each sensor', z_bins, z_inf, z_sup)
195  # Sets the number of resampled samples for bootstrap
196  Nrs = 100
197  # Initializes parameters to compute the median of ROOT TH1
198  q_U = array('d', [0])
199  p_U = array('d', [0.5])
200  q_V = array('d', [0])
201  p_V = array('d', [0.5])
202  # Lists storing medians of toy MC distributions for statistical bootstrap
203  l_U_median = []
204  l_V_median = []
205  # Lists storing estimated medians of projected DeltaRes distribution for every sensor along z
206  l_U_median_pos = []
207  l_V_median_pos = []
208  # Computes medians and related uncertainties (statistical bootstrapping) and produces plots
209  c_ZU = TCanvas('c_ZU_' + str(lyr_num), 'DeltaResUZ_' + str(lyr_num), 700, 500)
210  c_ZV = TCanvas('c_ZV_' + str(lyr_num), 'DeltaResVZ_' + str(lyr_num), 700, 500)
211  if(lyr_num == 1 or lyr_num == 3):
212  c_ZU.Divide(2, 1)
213  c_ZV.Divide(2, 1)
214  if(lyr_num == 4):
215  c_ZU.Divide(3, 1)
216  c_ZV.Divide(3, 1)
217  if(lyr_num == 5 or lyr_num == 6):
218  c_ZU.Divide(3, 2)
219  c_ZV.Divide(3, 2)
220  c_ZUMedians = TCanvas('c_ZUMedians_' + str(lyr_num), 'UMedians_' + str(lyr_num), 700, 500)
221  c_ZVMedians = TCanvas('c_ZVMedians_' + str(lyr_num), 'VMedians_' + str(lyr_num), 700, 500)
222  for i in range(0, z_bins):
223  xinf = z_inf + i * (z_sup - z_inf) / z_bins
224  xsup = z_inf + (i + 1) * (z_sup - z_inf) / z_bins
225  h_ZU.GetXaxis().SetRangeUser(xinf, xsup)
226  h_ZV.GetXaxis().SetRangeUser(xinf, xsup)
227  h_ZU.ProjectionY().GetQuantiles(1, q_U, p_U)
228  h_ZV.ProjectionY().GetQuantiles(1, q_V, p_V)
229  h_U = h_ZU.ProjectionY()
230  h_V = h_ZV.ProjectionY()
231  h_U.SetTitle(str(round(xinf, 3)) + ' (cm) < z < ' + str(round(xsup, 3)) + ' (cm)')
232  h_V.SetTitle(str(round(xinf, 3)) + ' (cm) < z < ' + str(round(xsup, 3)) + ' (cm)')
233  h_U.GetXaxis().SetRangeUser(-200, 200)
234  h_V.GetXaxis().SetRangeUser(-200, 200)
235  if(lyr_num == 1):
236  h_U.GetYaxis().SetRangeUser(0, 500)
237  h_V.GetYaxis().SetRangeUser(0, 500)
238  median_pos_U = TLine(q_U[0], 0, q_U[0], 500)
239  median_pos_V = TLine(q_V[0], 0, q_V[0], 500)
240  if(lyr_num == 3):
241  h_U.GetYaxis().SetRangeUser(0, 2000)
242  h_V.GetYaxis().SetRangeUser(0, 2000)
243  median_pos_U = TLine(q_U[0], 0, q_U[0], 2000)
244  median_pos_V = TLine(q_V[0], 0, q_V[0], 2000)
245  if(lyr_num == 4):
246  h_U.GetYaxis().SetRangeUser(0, 7000)
247  h_V.GetYaxis().SetRangeUser(0, 7000)
248  median_pos_U = TLine(q_U[0], 0, q_U[0], 7000)
249  median_pos_V = TLine(q_V[0], 0, q_V[0], 7000)
250  if(lyr_num == 5):
251  h_U.GetYaxis().SetRangeUser(0, 3000)
252  h_V.GetYaxis().SetRangeUser(0, 3000)
253  median_pos_U = TLine(q_U[0], 0, q_U[0], 3000)
254  median_pos_V = TLine(q_V[0], 0, q_V[0], 3000)
255  if(lyr_num == 6):
256  h_U.GetYaxis().SetRangeUser(0, 3000)
257  h_V.GetYaxis().SetRangeUser(0, 3000)
258  median_pos_U = TLine(q_U[0], 0, q_U[0], 3000)
259  median_pos_V = TLine(q_V[0], 0, q_V[0], 3000)
260  median_pos_U.SetLineWidth(2)
261  median_pos_V.SetLineWidth(2)
262  median_pos_U.SetLineColor(2)
263  median_pos_V.SetLineColor(2)
264  l_U_median_pos.append(median_pos_U)
265  l_V_median_pos.append(median_pos_V)
266  h_U.GetYaxis().SetTitle('counts')
267  h_V.GetYaxis().SetTitle('counts')
268  meas_U = hist2array(h_U)
269  meas_V = hist2array(h_V)
270  bs_U = numpy.random.poisson(1., (len(meas_U), Nrs))
271  bs_V = numpy.random.poisson(1., (len(meas_V), Nrs))
272  for j in range(Nrs):
273  toy_U = numpy.repeat(meas_U, bs_U[:, j])
274  toy_V = numpy.repeat(meas_V, bs_V[:, j])
275  median_U_toy = numpy.median(toy_U)
276  median_V_toy = numpy.median(toy_V)
277  l_U_median.append(median_U_toy)
278  l_V_median.append(median_V_toy)
279  median_U_rs = numpy.array(l_U_median)
280  median_V_rs = numpy.array(l_V_median)
281  median_U_dev = numpy.std(median_U_rs)
282  median_V_dev = numpy.std(median_V_rs)
283  h_UMedians.SetBinContent(i + 1, q_U[0])
284  h_UMedians.SetBinError(i + 1, median_U_dev)
285  h_VMedians.SetBinContent(i + 1, q_V[0])
286  h_VMedians.SetBinError(i + 1, median_V_dev)
287  c_ZU.cd(i + 1)
288  h_U.DrawCopy()
289  l_U_median_pos[i].Draw("SAME")
290  c_ZV.cd(i + 1)
291  h_V.DrawCopy()
292  l_V_median_pos[i].Draw("SAME")
293  c_ZUMedians.cd()
294  h_UMedians.GetXaxis().SetTitle('z (cm)')
295  h_UMedians.GetYaxis().SetTitle('Median of #Deltares_{u} (#mum)')
296  h_UMedians.Draw()
297  c_ZVMedians.cd()
298  h_VMedians.GetXaxis().SetTitle('z (cm)')
299  h_VMedians.GetYaxis().SetTitle('Median of #Deltares_{V} (#mum)')
300  h_VMedians.Draw()
301  # If not existing, creates a dedicated folder
302  if not os.path.exists('Median_plots_OverlapsZ'):
303  os.mkdir('Median_plots_OverlapsZ')
304  c_ZU.SaveAs('Median_plots_OverlapsZ/Median_and_DeltaResUZ_Lyr' + str(lyr_num) + '.root')
305  c_ZU.SaveAs('Median_plots_OverlapsZ/Median_and_DeltaResUZ_Lyr' + str(lyr_num) + '.pdf')
306  c_ZV.SaveAs('Median_plots_OverlapsZ/Median_and_DeltaResVZ_Lyr' + str(lyr_num) + '.root')
307  c_ZV.SaveAs('Median_plots_OverlapsZ/Median_and_DeltaResVZ_Lyr' + str(lyr_num) + '.pdf')
308  c_ZUMedians.SaveAs('Median_plots_OverlapsZ/Lyr' + str(lyr_num) + '_DeltaResUMedians_vs_z.root')
309  c_ZUMedians.SaveAs('Median_plots_OverlapsZ/Lyr' + str(lyr_num) + '_DeltaResUMedians_vs_z.pdf')
310  c_ZVMedians.SaveAs('Median_plots_OverlapsZ/Lyr' + str(lyr_num) + '_DeltaResVMedians_vs_z.root')
311  c_ZVMedians.SaveAs('Median_plots_OverlapsZ/Lyr' + str(lyr_num) + '_DeltaResVMedians_vs_z.pdf')
312  return
313 
314 
315 def LayerSensorPlots(filename, lyr_num, lddr_num, snsr_num):
316  """
317  Creates and saves Layer.Sensor plots for overlapping hits hitmaps
318  """
319  # Accesses the TDirectory containing VXD hit-maps for overlapping hits
320  f = ROOT.TFile.Open(filename, 'read')
321  hm = f.Get('HitMaps_VXDOverlaps')
322  # If not existing, creates a dedicated directory
323  if not os.path.exists('HitMaps_plots_Overlaps'):
324  os.mkdir('HitMaps_plots_Overlaps')
325  # Produces Layer.Sensor plots containing hit-maps for overlapping hits
326  # in all the ladders of a specific layer
327  for i in range(1, snsr_num + 1):
328  c_Meas = TCanvas('c_Meas_' + str(lyr_num) + ':' + str(i), 'Layer:Sensor = ' + str(lyr_num) + ':' + str(i), 500, 700)
329  if(lyr_num == 6):
330  c_Meas.Divide(4, 4)
331  elif(lyr_num == 5):
332  c_Meas.Divide(4, 3)
333  elif(lyr_num == 4):
334  c_Meas.Divide(5, 2)
335  elif(lyr_num == 3):
336  c_Meas.Divide(7, 1)
337  elif(lyr_num == 2):
338  c_Meas.Divide(6, 2)
339  elif(lyr_num == 1):
340  c_Meas.Divide(4, 2)
341  for k in range(1, lddr_num + 1):
342  histo = hm.Get('h_' + str(lyr_num) + str(k) + str(i))
343  c_Meas.cd(k)
344  histo.Draw('COLZ')
345  c_Meas.SaveAs('HitMaps_plots_Overlaps/c_Layer:Sensor_' + str(lyr_num) + str(i) + '.root')
346  c_Meas.SaveAs('HitMaps_plots_Overlaps/c_Layer:Sensor_' + str(lyr_num) + str(i) + '.pdf')
347  return
348 
349 
350 # Root output of module OverlapResiduals
351 filename = 'histofile.root'
352 
353 if __name__ == "__main__":
354  # Dictionary for VXD layers with overlaps
355  VXDLayers = {1: {'Layer': 1, 'Ladders': 8, 'Sensors': 2, 'Phi_bins': 8,
356  'Phi_inf': -3.2, 'Phi_sup': 3.2, 'Z_bins': 2, 'Z_inf': -3.2, 'Z_sup': 5.9},
357  3: {'Layer': 3, 'Ladders': 7, 'Sensors': 2, 'Phi_bins': 7,
358  'Phi_inf': -3.0, 'Phi_sup': 3.0, 'Z_bins': 2, 'Z_inf': -9.5, 'Z_sup': 15.5},
359  4: {'Layer': 4, 'Ladders': 10, 'Sensors': 3, 'Phi_bins': 10,
360  'Phi_inf': -3.0, 'Phi_sup': 3.0, 'Z_bins': 3, 'Z_inf': -16.5, 'Z_sup': 21.5},
361  5: {'Layer': 5, 'Ladders': 12, 'Sensors': 4, 'Phi_bins': 13,
362  'Phi_inf': -3.2, 'Phi_sup': 3.2, 'Z_bins': 4, 'Z_inf': -20.5, 'Z_sup': 29.5},
363  6: {'Layer': 6, 'Ladders': 16, 'Sensors': 5, 'Phi_bins': 17,
364  'Phi_inf': -3.3, 'Phi_sup': 3.3, 'Z_bins': 5, 'Z_inf': -25.5, 'Z_sup': 36.5}}
365  # Calls the defined functions
366  for i in range(1, 7):
367  if(i == 2):
368  continue # No overlaps for layer 2 in Phase3
369  else:
370  LayerSensorPlots(
371  filename,
372  lyr_num=VXDLayers[i]['Layer'],
373  lddr_num=VXDLayers[i]['Ladders'],
374  snsr_num=VXDLayers[i]['Sensors'])
375  Median_plots_phi(
376  filename,
377  lyr_num=VXDLayers[i]['Layer'],
378  phi_bins=VXDLayers[i]['Phi_bins'],
379  phi_inf=VXDLayers[i]['Phi_inf'],
380  phi_sup=VXDLayers[i]['Phi_sup'])
381  Median_plots_z(
382  filename,
383  lyr_num=VXDLayers[i]['Layer'],
384  z_bins=VXDLayers[i]['Z_bins'],
385  z_inf=VXDLayers[i]['Z_inf'],
386  z_sup=VXDLayers[i]['Z_sup'])