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