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 brash 1.1.2.9 SUBROUTINE h_fpp_FP2DC(iSet,iCham,iLay,Slope,FPcoords,DCcoords)
|
70 frw 1.1.2.1 *--------------------------------------------------------
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 * 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 brash 1.1.2.9 integer*4 iCham
90 integer*4 iLay
91 integer*4 iLocalPlane
|
92 frw 1.1.2.1 logical*4 Slope
93 real*4 FPcoords(3), DCcoords(3)
94
95 integer*4 i,j
96 real*4 MYcoords(3)
97
|
98 brash 1.1.2.9 iLocalPlane=(iSet-1)*6+(iCham-1)*3+iLay
|
99 frw 1.1.2.1
100 if (Slope) then
101 * * for slopes, we can ignore any position offset
102 MYcoords(1) = FPcoords(1)
103 MYcoords(2) = FPcoords(2)
104 MYcoords(3) = FPcoords(3)
105 else
106 * * for coordinates, we need to subtract the offset
|
107 brash 1.1.2.9 MYcoords(1) = FPcoords(1) - HFPP_Xoff(iSet) -
108 & HFPP_Xoff_fine(iLocalPlane)
109 MYcoords(2) = FPcoords(2) - HFPP_Yoff(iSet) -
110 & HFPP_Yoff_fine(iLocalPlane)
|
111 frw 1.1.2.1 MYcoords(3) = FPcoords(3) - HFPP_Zoff(iSet)
112 endif
113
114 * * use rotation matrix to rotate from focal plane coords to DCset coords
115 do i=1,3 !x,y,z for DC
116 DCcoords(i) = 0.0
117 do j=1,3 !x,y,z for FP
118 DCcoords(i) = DCcoords(i) + HFPP_Mrotation(iSet,i,j) * MYcoords(j)
119 enddo
120 enddo
121
122 if (slope) then
123 * * for slopes, we need to renormalize to dz=1
124 if (DCcoords(3).eq.0.0) then
125 DCcoords(1) = H_FPP_BAD_COORD
126 DCcoords(2) = H_FPP_BAD_COORD
127 else
128 DCcoords(1) = DCcoords(1) / DCcoords(3)
129 DCcoords(2) = DCcoords(2) / DCcoords(3)
130 DCcoords(3) = 1.0
131 endif
132 frw 1.1.2.1 endif
133
134
135 RETURN
136 END
137
138
139 c==============================================================================
140 c==============================================================================
141 c==============================================================================
142 c==============================================================================
143
144
145 SUBROUTINE h_fpp_DC2FP(iSet,Slope,DCcoords,FPcoords)
146 *--------------------------------------------------------
147 * Hall C HMS Focal Plane Polarimeter Code
148 *
149 * Purpose: transforms coordinates from the coord system
150 * of the specified set of FPP drift chambers to
151 * the the HMS focal plane system
152 * alternatively transforms SLOPES
153 frw 1.1.2.1 *
154 * Created by Frank R. Wesselmann, February 2004
155 *
156 *--------------------------------------------------------
157
158 IMPLICIT NONE
159
160 INCLUDE 'hms_data_structures.cmn'
161 INCLUDE 'hms_geometry.cmn'
162 c INCLUDE 'hms_fpp_event.cmn'
163
164 integer*4 iSet
165 logical*4 Slope
166 real*4 DCcoords(3), FPcoords(3)
167
168 integer*4 i,j
169
170
171 * * use INVERSE rotation matrix to rotate from DCset coords to focal plane coords
172 do i=1,3 !x,y,z for FP
173 FPcoords(i) = 0.0
174 frw 1.1.2.1 do j=1,3 !x,y,z for DC
175 FPcoords(i) = FPcoords(i) + HFPP_Irotation(iSet,i,j) * DCcoords(j)
176 enddo
177 enddo
178
179 if (Slope) then
180 * * for slopes, we need to renormalize to dz=1 if possible
181 if (FPcoords(3).eq.0.0) then
182 FPcoords(1) = H_FPP_BAD_COORD
183 FPcoords(2) = H_FPP_BAD_COORD
184 else
185 FPcoords(1) = FPcoords(1) / FPcoords(3)
186 FPcoords(2) = FPcoords(2) / FPcoords(3)
187 FPcoords(3) = 1.0
188 endif
189
190 else
191 * * for coordinates, we need to add the offset
192 FPcoords(1) = FPcoords(1) + HFPP_Xoff(iSet)
193 FPcoords(2) = FPcoords(2) + HFPP_Yoff(iSet)
194 FPcoords(3) = FPcoords(3) + HFPP_Zoff(iSet)
195 frw 1.1.2.1 endif
196
197 RETURN
198 END
199
200
201 c==============================================================================
202 c==============================================================================
203 c==============================================================================
204 c==============================================================================
205
206
|
207 cdaq 1.1.2.8 SUBROUTINE h_fpp_closest(Track1,Track2,sclose,zclose)
|
208 frw 1.1.2.1 *--------------------------------------------------------
209 * Hall C HMS Focal Plane Polarimeter Code
210 *
211 * Purpose: given two lines (tracks) in space, determine
212 * the distance of closest approach and the
213 * average z-coordinate of the closest approach
214 *
215 * Created by Frank R. Wesselmann, February 2004
216 *
217 *--------------------------------------------------------
218
219 IMPLICIT NONE
220
221 INCLUDE 'hms_data_structures.cmn'
222
223 real*4 Track1(4), Track2(4) ! IN mx,bx,my,by of two tracks
224 real*4 sclose ! OUT distance at closest approach
225 real*4 zclose ! OUT average z-coordinate at c.a.
226
|
227 brash 1.1.2.4 real*8 mx1,my1,bx1,by1
228 real*8 mx2,my2,bx2,by2
229 real*8 a1,a2,b,c1,c2
230 real*8 x1,x2,y1,y2,z1,z2
231 real*8 denom
|
232 frw 1.1.2.1
|
233 cdaq 1.1.2.8 mx1 = dble(Track1(1))
234 bx1 = dble(Track1(2))
235 my1 = dble(Track1(3))
236 by1 = dble(Track1(4))
|
237 frw 1.1.2.1
|
238 cdaq 1.1.2.8 mx2 = dble(Track2(1))
239 bx2 = dble(Track2(2))
240 my2 = dble(Track2(3))
241 by2 = dble(Track2(4))
|
242 frw 1.1.2.1
243 a1 = mx1*mx1 + my1*my1 + 1
244 a2 = mx2*mx2 + my2*my2 + 1
245 b = mx1*mx2 + my1*my2 + 1
|
246 cdaq 1.1.2.8 c1 = (bx1 - bx2)*mx1 + (by1 - by2)*my1
247 c2 = (bx1 - bx2)*mx2 + (by1 - by2)*my2
|
248 frw 1.1.2.1
249 denom = b*b - a1*a2
|
250 brash 1.1.2.4 if (denom.eq.0.0d0) then
|
251 frw 1.1.2.1 zclose = H_FPP_BAD_COORD
252 sclose = H_FPP_BAD_COORD
253 else
254 z1 = (a2*c1 - b*c2) / denom
255 z2 = ( b*c1 - a1*c2) / denom
256
257 x1 = z1 * mx1 + bx1
258 y1 = z1 * my1 + by1
259
260 x2 = z2 * mx2 + bx2
261 y2 = z2 * my2 + by2
262
|
263 cdaq 1.1.2.8 zclose = sngl(0.5d0*(z1 + z2))
264 sclose = sngl(sqrt( (x1-x2)**2 + (y1-y2)**2 + (z1-z2)**2 ))
|
265 frw 1.1.2.1 endif
266
|
267 brash 1.1.2.4 c write(*,*)'Zclose calculation 1: ',mx1,mx2,my1,my2
268 c write(*,*)'Zclose calculation 2: ',bx1,bx2,by1,by2
269 c write(*,*)'Zclose calculation 2a:',b,denom
270 c write(*,*)'Zclose calculation 2b:',a1,a2,c1,c2
271 c write(*,*)'Zclose calculation 3: ',ztrack2,z1,z2
272 c write(*,*)'Zclose calculation 4: ',x1,x2,y1,y2
273 c write(*,*)'Zclose calculation 5: ',zclose,sclose
274
|
275 frw 1.1.2.1 RETURN
276 END
277
|
278 brash 1.1.2.6 c==============================================================================
279 c==============================================================================
280 c==============================================================================
281 c==============================================================================
282
283
284 SUBROUTINE h_fpp_conetest(Track1,DCset,zclose,theta,icone)
285 *--------------------------------------------------------
286 * Hall C HMS Focal Plane Polarimeter Code
287 *
288 * Purpose: Calculate FPP conetest variable - assumes elliptical projection
289 * onto the last layer of the 2nd chamber in a set.
290 *
291 * Created by Edward J. Brash, September 2007
292 *
293 *--------------------------------------------------------
294
295 IMPLICIT NONE
296
297 INCLUDE 'hms_data_structures.cmn'
298 INCLUDE 'hms_geometry.cmn'
299 brash 1.1.2.6
300 real*4 Track1(4) ! IN mx,bx,my,by of front track
301 integer*4 DCset ! which chamber set we are in
302 real*4 zclose ! previously calculated z of closest approach
303 real*4 theta ! polar scattering angle
304 integer*4 icone ! Cone-test variable (1=pass, 0=fail)
305
306 real*8 mx1,my1,bx1,by1
307 real*4 ztrack2 ! central z-position of FPP set
308 real*4 zback_off ! offset from cntrl. pos. of last layer
309 real*8 zback,xfront,yfront,ttheta
310 real*8 r1x,r1y,r2x,r2y,xmin,xmax,ymin,ymax
|
311 cdaq 1.1.2.7 real*8 xpt1,xpt2,xpt3,xpt4,ypt1,ypt2,ypt3,ypt4
312 integer*4 iSet,iChamber,iLayer,iPlane
|
313 brash 1.1.2.6
314 mx1 = Track1(1)*1.0d0
315 bx1 = Track1(2)*1.0d0
316 my1 = Track1(3)*1.0d0
317 by1 = Track1(4)*1.0d0
318
319 xmin=HFPP_Xoff(DCset)-HFPP_Xsize(DCset)/2.0
320 xmax=HFPP_Xoff(DCset)+HFPP_Xsize(DCset)/2.0
321 ymin=HFPP_Yoff(DCset)-HFPP_Ysize(DCset)/2.0
322 ymax=HFPP_Yoff(DCset)+HFPP_Ysize(DCset)/2.0
323
324 icone=1
325
326 iSet=DCset
327 do iChamber=1, H_FPP_N_DCINSET
328 do iLayer=1, H_FPP_N_DCLAYERS
329
330 iPlane = H_FPP_N_DCLAYERS * H_FPP_N_DCINSET * (iSet-1)
331 > + H_FPP_N_DCLAYERS * (iChamber-1)
332 > + iLayer
333
334 brash 1.1.2.6 zback_off = HFPP_layerZ(iSet,iChamber,iLayer)
335 ztrack2 = HFPP_Zoff(iSet)
336 zback = ztrack2 + zback_off
337
338 xfront = bx1 + mx1*zback
339 yfront = by1 + my1*zback
340
341 ttheta=tan(theta)
342
343 r1x = (zback-zclose)*(mx1 + (ttheta-mx1)/(1.0+ttheta*mx1))
344 r2x = (zback-zclose)*((ttheta+mx1)/(1.0-ttheta*mx1) - mx1)
345 r1y = (zback-zclose)*(my1 +
346 & (ttheta-my1)/(1.0+ttheta*my1))
347 r2y = (zback-zclose)*((ttheta+my1)/(1.0-ttheta*my1)
348 & - my1)
349
350 xpt1=xfront-abs(r1x)
351 ypt1=yfront
352 xpt2=xfront+abs(r2x)
353 ypt2=yfront
354 xpt3=xfront
355 brash 1.1.2.6 ypt3=yfront-abs(r1y)
356 xpt4=xfront
357 ypt4=yfront+abs(r2y)
358
359 if(xpt1.lt.xmin)icone=0
360 if(xpt2.gt.xmax)icone=0
361 if(ypt3.lt.ymin)icone=0
362 if(ypt4.gt.ymax)icone=0
363
364 c if(icone.eq.0) then
365 c write(*,*)'chamber limits ',xmin,xmax,ymin,ymax
366 c write(*,*)'(',xpt1,',',ypt1,')'
367 c write(*,*)'(',xpt2,',',ypt2,')'
368 c write(*,*)'(',xpt3,',',ypt3,')'
369 c write(*,*)'(',xpt4,',',ypt4,')'
370 c endif
371
372 enddo ! iLayer
373 enddo ! iChamber
374
375 RETURN
376 brash 1.1.2.6 END
|
377 frw 1.1.2.1
378 c==============================================================================
379 c==============================================================================
380 c==============================================================================
381 c==============================================================================
382
383
384 SUBROUTINE h_fpp_relative_angles(mx_ref,my_ref,mx_new,my_new,theta,phi)
385 *--------------------------------------------------------
386 * Hall C HMS Focal Plane Polarimeter Code
387 *
388 * Purpose: find the POLAR angles between two tracks
389 * reference track is likely the incident HMS track, and the
390 * new track is probably the FPP track
391 * tracks are only identified by horizontal and vertical slopes
392 * we rotate both tracks s.th. the ref track is the new z axis
393 * then the simple polar angles of the rotated "new" track are
394 * the ones we want
395 *
396 * Created by Frank R. Wesselmann, February 2004
397 *
398 frw 1.1.2.1 *--------------------------------------------------------
399
400 IMPLICIT NONE
401
402 INCLUDE 'hms_data_structures.cmn'
403
404 real*4 mx_ref,my_ref ! IN slopes of reference track (incident?)
405 real*4 mx_new,my_new ! IN slopes of interesting track (analyzer?)
406 real*4 theta,phi ! OUT _polar_ angles of new track relative to ref
407
408 real*8 alpha, beta ! horizontal and vertical angle of ref track
409 real*8 M(3,3) ! rotation matrix: (row,column)
410 real*8 r_i(3) ! unit vector along "new" track before rotation
411 real*8 r_f(3) ! unit vector along "new" track after rotation
|
412 brash 1.1.2.5 real*8 r_in(3) ! unit vector along "in" track before rotation
413 real*8 r_fin(3) ! unit vector along "in" track after rotation
|
414 frw 1.1.2.1 real*8 magnitude
415 real*8 dtheta,dphi ! for convenience, double precision versions of OUT
|
416 brash 1.1.2.5 real*8 x,y,z,xin,yin,zin
|
417 frw 1.1.2.1
418 integer i,j
419
420
421 * * figure out rotation matrix
422
|
423 brash 1.1.2.5 c write(*,*)'Theta calculation 1: ',mx_ref,mx_new,my_ref,my_new
|
424 brash 1.1.2.4
|
425 frw 1.1.2.1 beta = datan(dble(my_ref))
|
426 brash 1.1.2.5 alpha = datan(dble(mx_ref)*dcos(beta)) ! this ought to be safe as the negative angle works
|
427 frw 1.1.2.1
428 M(1,1) = dcos(alpha)
429 M(1,2) = -1.d0*dsin(alpha)*dsin(beta)
430 M(1,3) = -1.d0*dsin(alpha)*dcos(beta)
431
432 M(2,1) = 0.d0
433 M(2,2) = dcos(beta)
434 M(2,3) = -1.d0* dsin(beta)
435
436 M(3,1) = dsin(alpha)
437 M(3,2) = dcos(alpha)*dsin(beta)
438 M(3,3) = dcos(alpha)*dcos(beta)
439
|
440 brash 1.1.2.5 * * normalize incoming vector
441
442 xin = dble(mx_ref)
443 yin = dble(my_ref)
444 zin = 1.d0
445 magnitude = dsqrt(xin*xin+yin*yin+zin*zin)
446 r_in(1)=xin/magnitude
447 r_in(2)=yin/magnitude
448 r_in(3)=zin/magnitude
449
450 do i=1,3
451 r_fin(i) = 0.d0
452 do j=1,3
453 r_fin(i) = r_fin(i) + M(i,j)*r_in(j)
454 enddo !j
455 enddo !i
456
457 c write(*,*)r_in(1),r_in(2),r_in(3)
458 c write(*,*)r_fin(1),r_fin(2),r_fin(3)
459
|
460 frw 1.1.2.1 * * normalize direction vector
461
462 x = dble(mx_new)
|
463 brash 1.1.2.4 y = dble(my_new)
|
464 frw 1.1.2.1 z = 1.d0
465 magnitude = dsqrt(x*x + y*y + z*z)
466 r_i(1) = x / magnitude
467 r_i(2) = y / magnitude
468 r_i(3) = z / magnitude
469
|
470 brash 1.1.2.4
|
471 frw 1.1.2.1 * * rotate vector of interest
472
473 do i=1,3
474 r_f(i) = 0.d0
475 do j=1,3
476 r_f(i) = r_f(i) + M(i,j)*r_i(j)
477 enddo !j
478 enddo !i
479
|
480 brash 1.1.2.5 c write(*,*)r_i(1),r_i(2),r_i(3)
|
481 brash 1.1.2.4 c write(*,*)r_f(1),r_f(2),r_f(3)
482
|
483 frw 1.1.2.1 * * find polar angles
484
485 dtheta = dacos(r_f(3)) ! z = cos(theta)
486
487 if (r_f(1).ne.0.d0) then
|
488 brash 1.1.2.3 if (r_f(1).gt.0.d0) then
489 if (r_f(2).gt.0.d0) then
490 dphi = datan( r_f(2)/r_f(1) ) ! y/x = tan(phi)
491 else
492 dphi = datan( r_f(2)/r_f(1) ) ! y/x = tan(phi)
493 dphi = dphi + 6.28318d0
494 endif
495 else
496 dphi = datan( r_f(2)/r_f(1) ) ! y/x = tan(phi)
497 dphi = dphi + 3.14159d0
498 endif
|
499 frw 1.1.2.1 elseif (r_f(2).gt.0.d0) then
|
500 brash 1.1.2.3 dphi = 1.57080d0 ! phi = +90
|
501 frw 1.1.2.1 elseif (r_f(2).lt.0.d0) then
|
502 brash 1.1.2.3 phi = 4.71239d0 ! phi = +270
|
503 frw 1.1.2.1 else
504 dphi = 0.d0 ! phi undefined if theta=0 or r=0
505 endif
506
507 theta = sngl(dtheta)
508 phi = sngl(dphi)
|
509 brash 1.1.2.4
510 c write(*,*)'Theta, phi = ',theta*180.0/3.14159265,phi*180.0/3.14159265
|
511 frw 1.1.2.1
512
513 RETURN
514 END
|
Return to
h_fpp_geometry.f
CVS log
Up to [HallC] / Analyzer / HTRACKING