Belle II Software development
VXDOverlaps.py
1
8
9import os
10import ROOT
11import numpy
12from array import array
13from ROOT import TCanvas, TH1F, TLine
14from 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
22def 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
177def 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
315def 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
351filename = 'histofile.root'
352
353if __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'])