1 frw 1.1.2.1 *--------------------------------------------------------
2 *--------------------------------------------------------
3 *--------------------------------------------------------
4 *
5 * Hall C HMS Focal Plane Polarimeter Code
6 *
7 * Created by Frank R. Wesselmann, February 2004
8 *
9 *--------------------------------------------------------
10 *
11 * this file contains several small geometry related routines
12 *
13 *--------------------------------------------------------
14
15
16
17 SUBROUTINE h_fpp_uTrack(iSet,iCham,iLay,iTrack,uCoord)
18 *--------------------------------------------------------
19 * Hall C HMS Focal Plane Polarimeter Code
20 *
21 * Purpose: determine in-layer coordinate of intersection
22 frw 1.1.2.1 * of given track and given drift chamber layer
23 *
24 * Created by Frank R. Wesselmann, February 2004
25 *
26 *--------------------------------------------------------
27
28 IMPLICIT NONE
29
30 INCLUDE 'hms_data_structures.cmn'
31 INCLUDE 'hms_geometry.cmn'
|
32 frw 1.1.2.2 include 'gen_detectorids.par'
|
33 frw 1.1.2.1 include 'gen_decode_common.cmn'
34 INCLUDE 'hms_fpp_event.cmn'
35
36 integer*4 iSet, iCham, iLay, iTrack
37 real*4 uCoord
38
39 real*4 x,y,z
40
41 uCoord = H_FPP_BAD_COORD
42
43 if (HFPP_N_tracks(iSet).le.0) RETURN
44 if (iTrack.le.0.or.iTrack.gt.HFPP_N_tracks(iSet)) RETURN
45
46 if (iSet.le.0.or.iSet.gt.H_FPP_N_DCSETS) RETURN
47 if (iCham.le.0.or.iCham.gt.H_FPP_N_DCINSET) RETURN
48 if (iLay.le.0.or.iLay.gt.H_FPP_N_DCLAYERS) RETURN
49
50 z = HFPP_layerZ(iSet,iCham,iLay)
51
52 x = HFPP_track_fine(iSet,iTrack,2) + HFPP_track_fine(iSet,iTrack,1)*z
53 y = HFPP_track_fine(iSet,iTrack,4) + HFPP_track_fine(iSet,iTrack,3)*z
54 frw 1.1.2.1
55 * * determine generalized in-layer coordinate
56 uCoord = HFPP_direction(iSet,iCham,iLay,1) * x
57 > + HFPP_direction(iSet,iCham,iLay,2) * y
58
59 RETURN
60 END
61
62
63 c==============================================================================
64 c==============================================================================
65 c==============================================================================
66 c==============================================================================
67
68
69 SUBROUTINE h_fpp_FP2DC(iSet,Slope,FPcoords,DCcoords)
70 *--------------------------------------------------------
71 * Hall C HMS Focal Plane Polarimeter Code
72 *
73 * Purpose: transforms coordinates from HMS focal plane
74 * system to the coord system of the specified
75 frw 1.1.2.1 * set of FPP drift chambers
76 * alternatively transforms SLOPES
77 *
78 * Created by Frank R. Wesselmann, February 2004
79 *
80 *--------------------------------------------------------
81
82 IMPLICIT NONE
83
84 INCLUDE 'hms_data_structures.cmn'
85 INCLUDE 'hms_geometry.cmn'
86 c INCLUDE 'hms_fpp_event.cmn'
87
88 integer*4 iSet
89 logical*4 Slope
90 real*4 FPcoords(3), DCcoords(3)
91
92 integer*4 i,j
93 real*4 MYcoords(3)
94
95
96 frw 1.1.2.1 if (Slope) then
97 * * for slopes, we can ignore any position offset
98 MYcoords(1) = FPcoords(1)
99 MYcoords(2) = FPcoords(2)
100 MYcoords(3) = FPcoords(3)
101 else
102 * * for coordinates, we need to subtract the offset
103 MYcoords(1) = FPcoords(1) - HFPP_Xoff(iSet)
104 MYcoords(2) = FPcoords(2) - HFPP_Yoff(iSet)
105 MYcoords(3) = FPcoords(3) - HFPP_Zoff(iSet)
106 endif
107
108 * * use rotation matrix to rotate from focal plane coords to DCset coords
109 do i=1,3 !x,y,z for DC
110 DCcoords(i) = 0.0
111 do j=1,3 !x,y,z for FP
112 DCcoords(i) = DCcoords(i) + HFPP_Mrotation(iSet,i,j) * MYcoords(j)
113 enddo
114 enddo
115
116 if (slope) then
117 frw 1.1.2.1 * * for slopes, we need to renormalize to dz=1
118 if (DCcoords(3).eq.0.0) then
119 DCcoords(1) = H_FPP_BAD_COORD
120 DCcoords(2) = H_FPP_BAD_COORD
121 else
122 DCcoords(1) = DCcoords(1) / DCcoords(3)
123 DCcoords(2) = DCcoords(2) / DCcoords(3)
124 DCcoords(3) = 1.0
125 endif
126 endif
127
128
129 RETURN
130 END
131
132
133 c==============================================================================
134 c==============================================================================
135 c==============================================================================
136 c==============================================================================
137
138 frw 1.1.2.1
139 SUBROUTINE h_fpp_DC2FP(iSet,Slope,DCcoords,FPcoords)
140 *--------------------------------------------------------
141 * Hall C HMS Focal Plane Polarimeter Code
142 *
143 * Purpose: transforms coordinates from the coord system
144 * of the specified set of FPP drift chambers to
145 * the the HMS focal plane system
146 * alternatively transforms SLOPES
147 *
148 * Created by Frank R. Wesselmann, February 2004
149 *
150 *--------------------------------------------------------
151
152 IMPLICIT NONE
153
154 INCLUDE 'hms_data_structures.cmn'
155 INCLUDE 'hms_geometry.cmn'
156 c INCLUDE 'hms_fpp_event.cmn'
157
158 integer*4 iSet
159 frw 1.1.2.1 logical*4 Slope
160 real*4 DCcoords(3), FPcoords(3)
161
162 integer*4 i,j
163
164
165 * * use INVERSE rotation matrix to rotate from DCset coords to focal plane coords
166 do i=1,3 !x,y,z for FP
167 FPcoords(i) = 0.0
168 do j=1,3 !x,y,z for DC
169 FPcoords(i) = FPcoords(i) + HFPP_Irotation(iSet,i,j) * DCcoords(j)
170 enddo
171 enddo
172
173 if (Slope) then
174 * * for slopes, we need to renormalize to dz=1 if possible
175 if (FPcoords(3).eq.0.0) then
176 FPcoords(1) = H_FPP_BAD_COORD
177 FPcoords(2) = H_FPP_BAD_COORD
178 else
179 FPcoords(1) = FPcoords(1) / FPcoords(3)
180 frw 1.1.2.1 FPcoords(2) = FPcoords(2) / FPcoords(3)
181 FPcoords(3) = 1.0
182 endif
183
184 else
185 * * for coordinates, we need to add the offset
186 FPcoords(1) = FPcoords(1) + HFPP_Xoff(iSet)
187 FPcoords(2) = FPcoords(2) + HFPP_Yoff(iSet)
188 FPcoords(3) = FPcoords(3) + HFPP_Zoff(iSet)
189 endif
190
191 RETURN
192 END
193
194
195 c==============================================================================
196 c==============================================================================
197 c==============================================================================
198 c==============================================================================
199
200
|
201 brash 1.1.2.4 SUBROUTINE h_fpp_closest(Track1,Track2,ztrack2,sclose,zclose)
|
202 frw 1.1.2.1 *--------------------------------------------------------
203 * Hall C HMS Focal Plane Polarimeter Code
204 *
205 * Purpose: given two lines (tracks) in space, determine
206 * the distance of closest approach and the
207 * average z-coordinate of the closest approach
208 *
209 * Created by Frank R. Wesselmann, February 2004
210 *
211 *--------------------------------------------------------
212
213 IMPLICIT NONE
214
215 INCLUDE 'hms_data_structures.cmn'
216
217 real*4 Track1(4), Track2(4) ! IN mx,bx,my,by of two tracks
|
218 brash 1.1.2.4 real*4 ztrack2 ! central z-position of FPP set
|
219 frw 1.1.2.1 real*4 sclose ! OUT distance at closest approach
220 real*4 zclose ! OUT average z-coordinate at c.a.
221
|
222 brash 1.1.2.4 real*8 mx1,my1,bx1,by1
223 real*8 mx2,my2,bx2,by2
224 real*8 a1,a2,b,c1,c2
225 real*8 x1,x2,y1,y2,z1,z2
226 real*8 denom
|
227 frw 1.1.2.1
228
|
229 brash 1.1.2.4 mx1 = Track1(1)*1.0d0
230 bx1 = Track1(2)*1.0d0
231 my1 = Track1(3)*1.0d0
232 by1 = Track1(4)*1.0d0
|
233 frw 1.1.2.1
|
234 brash 1.1.2.4 mx2 = Track2(1)*1.0d0
235 bx2 = Track2(2)*1.0d0
236 my2 = Track2(3)*1.0d0
237 by2 = Track2(4)*1.0d0
|
238 frw 1.1.2.1
239 a1 = mx1*mx1 + my1*my1 + 1
240 a2 = mx2*mx2 + my2*my2 + 1
241 b = mx1*mx2 + my1*my2 + 1
|
242 brash 1.1.2.4 c1 = (bx1 - bx2)*mx1 + (by1 - by2)*my1 - ztrack2
243 c2 = (bx1 - bx2)*mx2 + (by1 - by2)*my2 - ztrack2
|
244 frw 1.1.2.1
245 denom = b*b - a1*a2
|
246 brash 1.1.2.4 if (denom.eq.0.0d0) then
|
247 frw 1.1.2.1 zclose = H_FPP_BAD_COORD
248 sclose = H_FPP_BAD_COORD
249 else
250 z1 = (a2*c1 - b*c2) / denom
251 z2 = ( b*c1 - a1*c2) / denom
252
253 x1 = z1 * mx1 + bx1
254 y1 = z1 * my1 + by1
255
256 x2 = z2 * mx2 + bx2
257 y2 = z2 * my2 + by2
258
|
259 brash 1.1.2.4 zclose = 0.5d0*(z1 + z2 + ztrack2)
260 sclose = sqrt( (x1-x2)**2 + (y1-y2)**2 + (z1-z2-ztrack2)**2 )
|
261 frw 1.1.2.1 endif
262
|
263 brash 1.1.2.4 c write(*,*)'Zclose calculation 1: ',mx1,mx2,my1,my2
264 c write(*,*)'Zclose calculation 2: ',bx1,bx2,by1,by2
265 c write(*,*)'Zclose calculation 2a:',b,denom
266 c write(*,*)'Zclose calculation 2b:',a1,a2,c1,c2
267 c write(*,*)'Zclose calculation 3: ',ztrack2,z1,z2
268 c write(*,*)'Zclose calculation 4: ',x1,x2,y1,y2
269 c write(*,*)'Zclose calculation 5: ',zclose,sclose
270
|
271 frw 1.1.2.1 RETURN
272 END
273
274
275 c==============================================================================
276 c==============================================================================
277 c==============================================================================
278 c==============================================================================
279
280
281 SUBROUTINE h_fpp_relative_angles(mx_ref,my_ref,mx_new,my_new,theta,phi)
282 *--------------------------------------------------------
283 * Hall C HMS Focal Plane Polarimeter Code
284 *
285 * Purpose: find the POLAR angles between two tracks
286 * reference track is likely the incident HMS track, and the
287 * new track is probably the FPP track
288 * tracks are only identified by horizontal and vertical slopes
289 * we rotate both tracks s.th. the ref track is the new z axis
290 * then the simple polar angles of the rotated "new" track are
291 * the ones we want
292 frw 1.1.2.1 *
293 * Created by Frank R. Wesselmann, February 2004
294 *
295 *--------------------------------------------------------
296
297 IMPLICIT NONE
298
299 INCLUDE 'hms_data_structures.cmn'
300
301 real*4 mx_ref,my_ref ! IN slopes of reference track (incident?)
302 real*4 mx_new,my_new ! IN slopes of interesting track (analyzer?)
303 real*4 theta,phi ! OUT _polar_ angles of new track relative to ref
304
305 real*8 alpha, beta ! horizontal and vertical angle of ref track
306 real*8 M(3,3) ! rotation matrix: (row,column)
307 real*8 r_i(3) ! unit vector along "new" track before rotation
308 real*8 r_f(3) ! unit vector along "new" track after rotation
309 real*8 magnitude
310 real*8 dtheta,dphi ! for convenience, double precision versions of OUT
311 real*8 x,y,z
312
313 frw 1.1.2.1 integer i,j
314
315
316 * * figure out rotation matrix
317
|
318 brash 1.1.2.4 c write(*,*)mx_ref,my_ref,mx_new,my_new
319
|
320 frw 1.1.2.1 alpha = datan(dble(mx_ref)) ! this ought to be safe as the negative angle works
321 beta = datan(dble(my_ref))
322
323 M(1,1) = dcos(alpha)
324 M(1,2) = -1.d0*dsin(alpha)*dsin(beta)
325 M(1,3) = -1.d0*dsin(alpha)*dcos(beta)
326
327 M(2,1) = 0.d0
328 M(2,2) = dcos(beta)
329 M(2,3) = -1.d0* dsin(beta)
330
331 M(3,1) = dsin(alpha)
332 M(3,2) = dcos(alpha)*dsin(beta)
333 M(3,3) = dcos(alpha)*dcos(beta)
334
335 * * normalize direction vector
336
337 x = dble(mx_new)
|
338 brash 1.1.2.4 y = dble(my_new)
|
339 frw 1.1.2.1 z = 1.d0
340 magnitude = dsqrt(x*x + y*y + z*z)
341 r_i(1) = x / magnitude
342 r_i(2) = y / magnitude
343 r_i(3) = z / magnitude
344
|
345 brash 1.1.2.4 c write(*,*)r_i(1),r_i(2),r_i(3)
346
|
347 frw 1.1.2.1 * * rotate vector of interest
348
349 do i=1,3
350 r_f(i) = 0.d0
351 do j=1,3
352 r_f(i) = r_f(i) + M(i,j)*r_i(j)
353 enddo !j
354 enddo !i
355
|
356 brash 1.1.2.4 c write(*,*)r_f(1),r_f(2),r_f(3)
357
|
358 frw 1.1.2.1 * * find polar angles
359
360 dtheta = dacos(r_f(3)) ! z = cos(theta)
361
362 if (r_f(1).ne.0.d0) then
|
363 brash 1.1.2.3 if (r_f(1).gt.0.d0) then
364 if (r_f(2).gt.0.d0) then
365 dphi = datan( r_f(2)/r_f(1) ) ! y/x = tan(phi)
366 else
367 dphi = datan( r_f(2)/r_f(1) ) ! y/x = tan(phi)
368 dphi = dphi + 6.28318d0
369 endif
370 else
371 dphi = datan( r_f(2)/r_f(1) ) ! y/x = tan(phi)
372 dphi = dphi + 3.14159d0
373 endif
|
374 frw 1.1.2.1 elseif (r_f(2).gt.0.d0) then
|
375 brash 1.1.2.3 dphi = 1.57080d0 ! phi = +90
|
376 frw 1.1.2.1 elseif (r_f(2).lt.0.d0) then
|
377 brash 1.1.2.3 phi = 4.71239d0 ! phi = +270
|
378 frw 1.1.2.1 else
379 dphi = 0.d0 ! phi undefined if theta=0 or r=0
380 endif
381
382 theta = sngl(dtheta)
383 phi = sngl(dphi)
|
384 brash 1.1.2.4
385 c write(*,*)'Theta, phi = ',theta*180.0/3.14159265,phi*180.0/3.14159265
|
386 frw 1.1.2.1
387
388 RETURN
389 END
|
Return to
h_fpp_geometry.f
CVS log
Up to [HallC] / Analyzer / HTRACKING