ATLAS Offline Software
Loading...
Searching...
No Matches
comparisonUtils.py
Go to the documentation of this file.
1# Functions needed by CompareGeometries
2# Author John Alison <johnda@hep.upenn.edu>
3
4# Function to fill the Barrel Histograms
5def fillBarrelHists(x1,y1,z1
6 ,x2,y2,z2
7 ,theHists
8 ):
9
10 #wkd2
11 #print x1,y1,z1,x2,y2,z2
12
13 # Get the nicknames for the hists
14 h_geo1 = theHists.geo1[1]
15 h_geo2 = theHists.geo2[1]
16 h_XvsR = theHists.DxVsR[1]
17 h_XvsZ = theHists.DxVsZ[1]
18 h_YvsR = theHists.DyVsR[1]
19 h_YvsZ = theHists.DyVsZ[1]
20 h_PhivsZ = theHists.DphiVsZ[1]
21 h_XYvsR = theHists.DxyVsR[1]
22 h_PhivsR = theHists.DphiVsR[1]
23 h_RvsR = theHists.DrVsR[1]
24 h_ZvsR = theHists.DzVsR[1]
25 h_XYvsZ = theHists.DxyVsZ[1]
26 h_RvsPhi = theHists.DrVsPhi[1]
27 h_X = theHists.Dx[1]
28 h_Y = theHists.Dy[1]
29 h_Z = theHists.Dz[1]
30
31 h_geo1.Fill(x1,y1)
32 h_geo2.Fill(x2,y2)
33
34 h_XvsR.Fill(sqrt(pow(x2,2) + pow(y2,2)),x1 - x2)
35 h_XvsZ.Fill(z2,x1 - x2)
36
37 h_YvsR.Fill(sqrt(pow(x2,2) + pow(y2,2)),y1 - y2)
38 h_YvsZ.Fill(z2,y1 - y2)
39
40 h_PhivsZ.Fill(z2,atan2(y1,x1)-atan2(y2,x2))
41
42 h_XYvsR.Fill(sqrt(pow(x2,2) + pow(y2,2)),sqrt(pow(x1 - x2,2)+pow(y1 - y2 ,2)))
43
44 #wkd2 - div by zero check
45 if not sqrt(x1*x1+y1*y1) == 0:
46 h_PhivsR.Fill(sqrt(pow(x2,2) + pow(y2,2)),1.0/sqrt(x1*x1+y1*y1)*((x1 - x2)*y1 - (y1 - y2)* x1))
47
48 h_RvsR.Fill(sqrt(pow(x2,2) + pow(y2,2)), (1.0/sqrt(x1*x1+y1*y1)*((x1 - x2)*x1 + (y1 - y2)* y1)))
49
50 h_X.Fill(x1 - x2)
51 h_Y.Fill(y1 - y2)
52 h_Z.Fill(z1 - z2)
53
54 if z1 > 0:
55 h_ZvsR.Fill(sqrt(pow(x2,2) + pow(y2,2)),sqrt(pow(z1 - z2,2)))
56
57 h_XYvsZ.Fill(z2,sqrt(pow(x1 - x2,2)+pow(y1 - y2 ,2)))
58
59 if x2:
60 h_RvsPhi.Fill(atan2(y2,x2),sqrt(pow(x1,2)+pow(y1,2)) -sqrt(pow(x2,2)+pow(y2,2)))
61
62 return
63
64# Function to fill the Endcap Histograms
65def fillEndcapHists(x1,y1,z1
66 ,x2,y2,z2
67 ,theHists
68 ,side
69 ):
70
71#---wkd2
72# print x1, y1, z1, x2, y2, z2
73#wkd2---
74
75 h_geo1 = theHists.geo1[side]
76 h_geo2 = theHists.geo2[side]
77 h_XvsR = theHists.DxVsR[side]
78 h_XvsZ = theHists.DxVsZ[side]
79 h_YvsR = theHists.DyVsR[side]
80 h_YvsZ = theHists.DyVsZ[side]
81 h_PhivsZ = theHists.DphiVsZ[side]
82 h_XYvsR = theHists.DxyVsR[side]
83 h_PhivsR = theHists.DphiVsR[side]
84 h_RvsR = theHists.DrVsR[side]
85 h_ZvsR = theHists.DzVsR[side]
86 h_XYvsZ = theHists.DxyVsZ[side]
87 h_RvsPhi = theHists.DrVsPhi[side]
88 h_X = theHists.Dx[side]
89 h_Y = theHists.Dy[side]
90 h_Z = theHists.Dz[side]
91
92 h_geo1.Fill(z1,atan2(y1,x1))
93 h_geo2.Fill(z2,atan2(y2,x2))
94
95 h_XvsR.Fill(sqrt(pow(x2,2) + pow(y2,2)),x1 - x2)
96 h_XvsZ.Fill(z2,x1 - x2)
97
98 h_YvsR.Fill(sqrt(pow(x2,2) + pow(y2,2)),y1 - y2)
99 h_YvsZ.Fill(z2,y1 - y2)
100
101 h_PhivsZ.Fill(z2,atan2(y1,x1)-atan2(y2,x2))
102
103 h_XYvsR.Fill(sqrt(pow(x2,2) + pow(y2,2)),sqrt(pow(x1 - x2,2)+pow(y1 - y2 ,2)))
104
105 #wkd2 - div by 0 check
106 if not sqrt(x1*x1+y1*y1) == 0:
107 h_PhivsR.Fill(sqrt(pow(x2,2) + pow(y2,2)),1.0/sqrt(x1*x1+y1*y1)*((x1 - x2)*y1 - (y1 - y2)* x1))
108
109 h_RvsR.Fill(sqrt(pow(x2,2) + pow(y2,2)), (1.0/sqrt(x1*x1+y1*y1)*((x1 - x2)*x1 + (y1 - y2)* y1)))
110
111 h_ZvsR.Fill(sqrt(pow(x2,2) + pow(y2,2)),sqrt(pow(z1 - z2,2)))
112
113 h_XYvsZ.Fill(z2,sqrt(pow(x1 - x2,2)+pow(y1 - y2 ,2)))
114
115
116 h_X.Fill(x1 - x2)
117 h_Y.Fill(y1 - y2)
118 h_Z.Fill(z1 - z2)
119
120 if x2:
121 h_RvsPhi.Fill(atan2(y2,x2),sqrt(pow(x1,2)+pow(y1,2)) -sqrt(pow(x2,2)+pow(y2,2)))
122
123 return
124
125# Function to make the geometry histograms
126def makeSubSysGeometryHist(name,title,color):
127
128 h_Ecc = TH2F(name+"_Ecc",title+" (Endcap C)",2400,-2900,-500,2400,-3.14,3.14)
129 h_Ecc.SetMarkerColor(color+1)
130 h_Ecc.GetXaxis().SetTitle("z [mm]")
131 h_Ecc.GetYaxis().SetTitle("phi")
132
133 h_b = TH2F(name,title,2400,-1200,1200,2400,-1200,1200)
134 h_b.SetMarkerColor(color+1)
135 h_b.GetXaxis().SetTitle("x [mm]")
136 h_b.GetYaxis().SetTitle("y [mm]")
137
138 h_Eca = TH2F(name+"_Eca",title+" (Endcap A)",2400,400,2900,2400,-3.14,3.14)
139 h_Eca.SetMarkerColor(color+1)
140 h_Eca.GetXaxis().SetTitle("z [mm]")
141 h_Eca.GetYaxis().SetTitle("phi")
142
143 return h_Ecc,h_b,h_Eca,
144
145# Make the difference hists for a particular sub system
146def makeSubSysDiffHist(name,title,xName,xMins,xMaxes,yName,yRange,color,nBins = 120):
147
148
149 h_Ecc = TH2F(name+"_Ecc",title+" (Endcap C)", nBins,xMins[0],xMaxes[0],nBins,-1*yRange,yRange)
150 h_Ecc.SetMarkerColor(color+1)
151 h_Ecc.GetXaxis().SetTitle(xName)
152 h_Ecc.GetYaxis().SetTitle(yName)
153
154 h_b = TH2F(name,title,nBins,xMins[1],xMaxes[1],nBins,-1*yRange,yRange)
155 h_b.SetMarkerColor(color+1)
156 h_b.GetXaxis().SetTitle(xName)
157 h_b.GetYaxis().SetTitle(yName)
158
159 h_Eca = TH2F(name+"_Eca",title+" (Endcap A)", nBins,xMins[2],xMaxes[2],nBins,-1*yRange,yRange)
160 h_Eca.SetMarkerColor(color+1)
161 h_Eca.GetXaxis().SetTitle(xName)
162 h_Eca.GetYaxis().SetTitle(yName)
163
164 return h_Ecc,h_b,h_Eca
165
166def make1D_SubSysDiffHist(name,title,xName,xMins,xMaxes,yName,color):
167 nBins = 120
168 h_Ecc = TH1F(name+"_Ecc",title+" (Endcap C)", nBins,xMins[0],xMaxes[0])
169 h_Ecc.SetMarkerColor(color+1)
170 h_Ecc.SetLineColor(color+1)
171 h_Ecc.GetXaxis().SetTitle(xName)
172 h_Ecc.GetYaxis().SetTitle(yName)
173
174 h_b = TH1F(name,title,nBins,xMins[1],xMaxes[1])
175 h_b.SetMarkerColor(color+1)
176 h_b.SetLineColor(color+1)
177 h_b.GetXaxis().SetTitle(xName)
178 h_b.GetYaxis().SetTitle(yName)
179
180 h_Eca = TH1F(name+"_Eca",title+" (Endcap A)", nBins,xMins[2],xMaxes[2])
181 h_Eca.SetMarkerColor(color+1)
182 h_Eca.SetLineColor(color+1)
183 h_Eca.GetXaxis().SetTitle(xName)
184 h_Eca.GetYaxis().SetTitle(yName)
185
186 return h_Ecc,h_b,h_Eca
187
188# Function to write out all the histograms
189def WriteHist(title
190 ,hists
191 ,drawEndcapLines=False
192 ,drawTRTFirst=False):
193
194 can = TCanvas(title,title, 700, 500)
195 can.Divide(3)
196
197 can.cd(1)
198 for i in range(len(hists)):
199 if i:
200 hists[i][0].Draw("same")
201 else:
202 hists[i][0].Draw("")
203
204
205 if drawEndcapLines:
206 lineLength = 0.0008
207 endcapLinesC_ForGeo = DrawTRTEndcapCLines(lineLength)
208
209
210 can.cd(2)
211 for i in range(len(hists)):
212 if i:
213 hists[i][1].Draw("same")
214 else:
215 hists[i][1].Draw("")
216
217 can.cd(3)
218 for i in range(len(hists)):
219 if i:
220 hists[i][2].Draw("same")
221 else:
222 hists[i][2].Draw("")
223
224
225 if drawEndcapLines:
226 lineLength = 0.0008
227 endcapLinesA_ForGeo = DrawTRTEndcapALines(lineLength)
228
229 can.Write()
230
231 return
232
233# Function to read in the data from the input text files
234def readInData(inputfile):
235 m_iblElements = [] # IBL separate because it has more eta modules than the PIX.
236 m_pixelElements = []
237 m_sctElements = []
238 m_trtElements = []
239 for line in inputfile:
240 words = line.split()
241
242 m_thisPosition = []
243 for i in range(len(words)-1):
244 m_thisPosition.append(float(words[i+1]))
245
246 # Pixel
247 if words[0] == "1":
248 if words[1] == "4" or words[1] == "-4":
249 # DBM, we're gonna lump it in with the IBL
250 m_iblElements.append(m_thisPosition)
251
252 elif words[1] == "0" and words[2] == "0": # IBL
253 m_iblElements.append(m_thisPosition)
254
255 else:
256 m_pixelElements.append(m_thisPosition)
257
258 # SCT
259 if words[0] == "2":
260 m_sctElements.append(m_thisPosition)
261
262 # TRT
263 if words[0] == "3":
264 m_trtElements.append(m_thisPosition)
265
266 return m_iblElements,m_pixelElements,m_sctElements,m_trtElements
267
268
269# Functions to draw lines on the TRT endcap plots
271
272 #B Wheels
273 bLines =[]
274 for i in range(9):
275 bWheels = -2710 + 122*i
276 bLines.append(TLine(bWheels,-1*yRange,bWheels,yRange))
277 bLines[i].SetLineStyle(2)
278 bLines[i].Draw()
279
280 seperator = TLine(-1725,-1*yRange,-1725,yRange)
281 seperator.SetLineColor(kRed)
282 seperator.Draw()
283
284 #A Wheels
285 aLines = []
286 do16plane = False
287 if do16plane:
288 for i in range(7):
289 # For the 16-plane wheel
290 aWheels = -1705 + 142 *i;
291 aLines.append(TLine(aWheels,-1*yRange,aWheels,yRange))
292 aLines[i].SetLineStyle(2)
293 aLines[i].Draw()
294 else:
295 for i in range(13):
296 # For the 8-plane wheel
297 aWheels = -1705 + 71 *i;
298 aLines.append(TLine(aWheels,-1*yRange,aWheels,yRange))
299 aLines[i].SetLineStyle(2)
300 aLines[i].Draw()
301
302 return bLines, seperator, aLines
303
305 #B Wheels
306 bLines =[]
307 for i in range(9):
308 bWheels = 2710 - 122*i
309 bLines.append(TLine(bWheels,-1*yRange,bWheels,yRange))
310 bLines[i].SetLineStyle(2)
311 bLines[i].Draw()
312
313 seperator = TLine(1725,-1*yRange,1725,yRange)
314 seperator.SetLineColor(kRed)
315 seperator.Draw()
316
317 #A Wheels
318 aLines = []
319 do16plane = False
320 if do16plane:
321 for i in range(7):
322 aWheels = 1705 - 142 *i;
323 aLines.append(TLine(aWheels,-1*yRange,aWheels,yRange))
324 aLines[i].SetLineStyle(2)
325 aLines[i].Draw()
326 else:
327 for i in range(13):
328 aWheels = 1705 - 71 *i;
329 aLines.append(TLine(aWheels,-1*yRange,aWheels,yRange))
330 aLines[i].SetLineStyle(2)
331 aLines[i].Draw()
332 return bLines, seperator, aLines
333
334
335
336def globalStrawLayer(wheel,strawlayer):
337 if wheel < 6:
338 return wheel*16 + strawlayer
339
340 return 6*16 + (wheel - 6)*8 + strawlayer
341
342def drawStrawPlaneTransCan(bec,wheel,strawlayer,x1,y1,x2,y2):
343
344 # 0 - 160
345 strawPlane = int(globalStrawLayer(wheel,strawlayer))
346
347 # Fill TRT Trans can
348 if bec == 2:
349 if strawPlane >= len(trtEndcapAStrawPlanes):
350 return
351
352 trtEndcapAStrawPlanes[strawPlane].cd()
353
354 line.SetLineColor(kGreen+1)
355 line.DrawArrow(x1,y1,x1+TRASL_FACTOR*(x1-x2),y1+TRASL_FACTOR*(y1-y2),0.01,"")
356
357 if bec == -2:
358 if strawPlane >= len(trtEndcapCStrawPlanes):
359 return
360
361 trtEndcapCStrawPlanes[strawPlane].cd()
362 line.SetLineColor(kGreen+1)
363 line.DrawArrow(x1,y1,x1+TRASL_FACTOR*(x1-x2),y1+TRASL_FACTOR*(y1-y2),0.01,"")
364 return
constexpr int pow(int base, int exp) noexcept
makeSubSysGeometryHist(name, title, color)
DrawTRTEndcapALines(yRange)
readInData(inputfile)
DrawTRTEndcapCLines(yRange)
globalStrawLayer(wheel, strawlayer)
makeSubSysDiffHist(name, title, xName, xMins, xMaxes, yName, yRange, color, nBins=120)
fillBarrelHists(x1, y1, z1, x2, y2, z2, theHists)
WriteHist(title, hists, drawEndcapLines=False, drawTRTFirst=False)
drawStrawPlaneTransCan(bec, wheel, strawlayer, x1, y1, x2, y2)
fillEndcapHists(x1, y1, z1, x2, y2, z2, theHists, side)
make1D_SubSysDiffHist(name, title, xName, xMins, xMaxes, yName, color)