ATLAS Offline Software
h6prod_getxy.py
Go to the documentation of this file.
1 # Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
2 
3 
4 from math import sqrt
5 import ROOT
6 from array import array
7 
8 
9 def getXYInPolygon(s,nev0=1500):
10  xset,yset,nset = getXYSetInPolygon(nev0)
11  ss = s%(len(xset))
12  return xset[ss], yset[ss], nset[ss]
13 
14 
15 
16 def drawXYInPolygon(nev0=1500):
17  nevtot = 0
18  nevinwin = 0
19  xset,yset,nset = getXYSetInPolygon(nev0)
20  print ("Number of points ", len(xset), len(yset))
21  #histXY = ROOT.TH2F( 'clusE', 'clusE', 500, -500.0, 500.0, 500, -500., 500. )
22  histXY = ROOT.TH2F( 'clusE', 'clusE', 2000, -800.0, 800.0, 2000, -800., 800. )
23  histXY.SetStats(0)
24  histXY.SetTitle("")
25  for i in range(0,len(xset)):
26  #print ("nrun:",i+1," xCryo:",xset[i], "yTable:",yset[i])
27  print ('nrun %3d xCryo %5d yTable %5d nev %5d' % (i+1,xset[i], yset[i], nset[i]))
28  nevtot = nevtot + nset[i]
29  dist = sqrt(xset[i]*xset[i] + (70-yset[i])*(70-yset[i]))
30  if dist < 300:
31  nevinwin = nevinwin + nset[i]
32 
33  histXY.Fill(xset[i], yset[i])
34  c0 = ROOT.TCanvas( 'c0', 'c0', 800, 800 )
35  c0.SetFillColor(10)
36  c0.SetFillStyle(0)
37  c0.cd()
38  ROOT.gPad.SetGrid()
39  ROOT.gPad.SetLeftMargin(0.12)
40  ROOT.gPad.SetRightMargin(0.12)
41  ROOT.gPad.SetTopMargin(0.12)
42  ROOT.gPad.SetBottomMargin(0.12)
43  tc = ROOT.TEllipse(0,70,300)
44  tc.SetFillStyle(0)
45  tc.SetLineWidth(2)
46  histXY.SetMarkerStyle(20)
47  histXY.SetLineColor(2)
48  histXY.SetMarkerColor(2)
49  histXY.Draw()
50  tc.Draw()
51  #poly_np = 9
52  #poly_x = array('d',[0, 300, 150, 250, 0, -250, -150, -300, 0])
53  #poly_y = array('d',[320, 220, -80, -130, -280, -130, -80, 220, 320])
54  poly_np = 8
55  poly_x = array('d',[0, 500, 700, 60, -60, -700, -500, 0])
56  poly_y = array('d',[700, 600, 400, -400, -400, 400, 600, 700])
57  myPoly = ROOT.TPolyLine(poly_np, poly_x, poly_y)
58  myPoly.SetLineColor(4)
59  myPoly.SetLineWidth(2)
60  myPoly.Draw()
61 
62  #for i in range(0,500):
63  #x,y = getXYInPolygon(i)
64  #print ("nrun:",i+1," xCryo:",x, "yTable:",y)
65 
66  print ("nevtot:",nevtot, "in win:",nevinwin)
67  return c0, histXY, tc, myPoly
68 
69 
70 #def getXYInPolygon(s):
71 
72 
73 
77 def getXYSetInPolygon(nev0=1500):
78  # this is our polygon drawn on xCryo, yTable plane in which we are going to place our cryostat and table
79  #poly_np = 9
80  #poly_x = array('i',[0, 300, 150, 250, 0, -250, -150, -300, 0])
81  #poly_y = array('i',[320, 220, -80, -130, -280, -130, -80, 220, 320])
82 
83  poly_np = 8
84  poly_x = array('i',[0, 500, 700, 60, -60, -700, -500, 0])
85  poly_y = array('i',[700, 600, 400, -400, -400, 400, 600, 700])
86 
87  # this is starting point of grid
88  x0 = 0
89  y0 = 0
90  # this is grid step
91  dx = 25
92  dy = 25
93  # this is number of events to simulate in points Y=-400 and Y=+400
94  #n_Ytop = 2000
95  #n_Ybot = 4000
96  n_Ytop = nev0
97  y_Ytop = 700
98  n_Ybot = nev0*2
99  y_Ybot = -400
100  print ("Initial number of events:", nev0, "y_Ytop:",y_Ytop, "n_Ytop:",n_Ytop, "y_Ybot:",y_Ybot, "n_Ybot:",n_Ybot)
101  # - - - - - -
102  # - 2 3 4 - -
103  # - 1 + 5 - -
104  # - 0 7 6 - -
105  # - - - - - -
106  #
107  xset = []
108  yset = []
109  nset = []
110  xset.append(x0)
111  yset.append(y0)
112  #nset.append((n_Ytop+(y0-400)*(n_Ytop-n_Ybot)/800)/100*100)
113  nset.append( (n_Ybot+(y0-y_Ybot)*(n_Ytop-n_Ybot)/(y_Ytop-y_Ybot))/10*10 )
114  for iRound in range(1,1000):
115  boxSide = 1+iRound*2
116  nPoints = (boxSide-1)*4
117  #print ("----------------------------------------")
118  #print ("iRound:", iRound, " boxSide ", boxSide, "nPoints ", nPoints)
119  nGoodPoints = 0
120  for iPoint in range(0, nPoints):
121  boxSideNum = iPoint/(boxSide-1)
122  boxSideStep = iPoint%(boxSide-1)
123  if boxSideNum == 0:
124  iX = x0 - iRound*dx
125  iY = y0 - iRound*dy + dy*boxSideStep
126  if boxSideNum == 1:
127  iX = x0 - iRound*dx + dx*boxSideStep
128  iY = y0 + iRound*dy
129  if boxSideNum == 2:
130  iX = x0 + iRound*dx
131  iY = y0 + iRound*dy - dy*boxSideStep
132  if boxSideNum == 3:
133  iX = x0 + iRound*dx - dx*boxSideStep
134  iY = y0 - iRound*dy
135  mt = ROOT.TMath
136  result = mt.IsInside(iX, iY, poly_np, poly_x, poly_y)
137  #print ("iPoint ", iPoint," from nPoints ",nPoints, " boxSideNum:",boxSideNum, " mod:", boxSideStep, " iX:", iX, " iY:", iY, result)
138  if result == 1:
139  nGoodPoints += 1
140  xset.append(iX)
141  yset.append(iY)
142  nset.append( (n_Ybot+(iY-y_Ybot)*(n_Ytop-n_Ybot)/(y_Ytop-y_Ybot))/10*10 )
143  #print ("nGoodPoints:", nGoodPoints, " nPoints:", nPoints)
144  if nGoodPoints == 0:
145  break
146 
147  return xset, yset, nset
148 
149 
150 #def getXYInRectangle(s):
151  #xmin=-240
152  #xmax=240
153  #ymin=-170
154  #ymax=181
155  #dx=100
156  #dy=100
157  #nx=(xmax-xmin)/dx + 1
158  #ny=(ymax-ymin)/dy + 1
159  #ss = s%(nx*ny)
160  #iy=ss/nx
161  #ix=ss%nx
162 
167 
168 
169 #def getXYInCircle(s):
170 
182 
183  #nx=(xmax-xmin)/dx + 1
184  #ny=(ymax-ymin)/dy + 1
185  #ntot = nx*ny
186  #xset = []
187  #yset = []
188  #for ix in range(0, int(nx)):
189  #for iy in range(0,int(ny)):
190  #x=xmin+ix*dx
191  #y=ymin+iy*dy
192  #dist = sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0))
193 
201 
python.h6prod_getxy.getXYInPolygon
def getXYInPolygon(s, nev0=1500)
Definition: h6prod_getxy.py:9
python.h6prod_getxy.drawXYInPolygon
def drawXYInPolygon(nev0=1500)
Definition: h6prod_getxy.py:16
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:194
python.h6prod_getxy.getXYSetInPolygon
def getXYSetInPolygon(nev0=1500)
Function returns two arrays with xCryo, yTable positions inside user defined shape (polygon) Points a...
Definition: h6prod_getxy.py:77
array