ATLAS Offline Software
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
5 def 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
65 def 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
126 def 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
146 def 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 
166 def 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
189 def 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
234 def 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 
336 def globalStrawLayer(wheel,strawlayer):
337  if wheel < 6:
338  return wheel*16 + strawlayer
339 
340  return 6*16 + (wheel - 6)*8 + strawlayer
341 
342 def 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
comparisonUtils.DrawTRTEndcapALines
def DrawTRTEndcapALines(yRange)
Definition: comparisonUtils.py:304
comparisonUtils.fillBarrelHists
def fillBarrelHists(x1, y1, z1, x2, y2, z2, theHists)
Definition: comparisonUtils.py:5
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
TH2F
Definition: rootspy.cxx:420
comparisonUtils.makeSubSysGeometryHist
def makeSubSysGeometryHist(name, title, color)
Definition: comparisonUtils.py:126
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
comparisonUtils.DrawTRTEndcapCLines
def DrawTRTEndcapCLines(yRange)
Definition: comparisonUtils.py:270
comparisonUtils.makeSubSysDiffHist
def makeSubSysDiffHist(name, title, xName, xMins, xMaxes, yName, yRange, color, nBins=120)
Definition: comparisonUtils.py:146
comparisonUtils.make1D_SubSysDiffHist
def make1D_SubSysDiffHist(name, title, xName, xMins, xMaxes, yName, color)
Definition: comparisonUtils.py:166
comparisonUtils.readInData
def readInData(inputfile)
Definition: comparisonUtils.py:234
comparisonUtils.fillEndcapHists
def fillEndcapHists(x1, y1, z1, x2, y2, z2, theHists, side)
Definition: comparisonUtils.py:65
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:195
comparisonUtils.WriteHist
def WriteHist(title, hists, drawEndcapLines=False, drawTRTFirst=False)
Definition: comparisonUtils.py:189
TH1F
Definition: rootspy.cxx:320
comparisonUtils.globalStrawLayer
def globalStrawLayer(wheel, strawlayer)
Definition: comparisonUtils.py:336
calibdata.cd
cd
Definition: calibdata.py:51
readCCLHist.float
float
Definition: readCCLHist.py:83
comparisonUtils.drawStrawPlaneTransCan
def drawStrawPlaneTransCan(bec, wheel, strawlayer, x1, y1, x2, y2)
Definition: comparisonUtils.py:342