1 jones 1.1 *------------------------------------------------------------------------
2 *
3 * HMS_TRACK HMS Tracking Routines
4 * -=========-
5 *
6 * Forward and Backward Tracking of electrons in the Jlab HMS hall
7 * C spectrometer
8 *
9 * Note: - the HMS routines use a lab (HMS) coord. system
10 * and the corresponding COSY coord. system, both
11 * right handed with
12 * x : pointing downwards
13 * y : perpendicular to x,z,
14 * pointing to the left (if seen in z-direction)
15 * z : HMS axis, pointing from the target to the focal plane
16 *
17 * - all lengths (x,y,z,l,...) are measured in [m]
18 * - all angles are measured as dx/dz,dy/dz (lab coords.)
19 * or as A,B (COSY coords.)
20 * - the momentum is measured in delta (relative momentum
21 * deviation = 1-p0/pHMS)
22 jones 1.1 *
23 * PART 1: Forward trackinh using COSY transport matrices
24 *
25 *
26 * PART 2: Reconstruction (backward tracking) using reconstruction
27 * and COSY transport matrices (including the
28 * effects of a vertical beam offset (out-of plane))
29 *
30 *
31 * written by Markus Muehlbauer for the GEN Experiment
32 *
33 * frw 9/2000
34 * changes made in the course of the migration to g77 compiler:
35 * - fixed some typos (old ones, too)
36 * - all COSY conversion code was already commented out (why?
37 * by who?) so I removed it and made the code more readable
38 * - various variables were not initialized prior to reading
39 * from file. This may or may not be an issue, but fixed it
40 * anyway
41 *
42 *------------------------------------------------------------------------
43 jones 1.1
44 *------------------------------------------------------------------------
45 *
46 * PART 1: HMS Forward Tracking (Target to Focal Plane)
47 * -=======-
48 *
49 * Forward tracking in the Jlab HMS hall C spectrometer
50 * using COSY transport matrices
51 *
52 * developed by Cris Cothran
53 * modified by Markus Muehlbauer
54 * - CCs orignal program converted into subroutines
55 * - mad additions for pure tracking, without checking the acceptance
56 * - and changed the innermost loops applying the matrix
57 * (which speeds up the whole thing by a factor of about 30)
58 *
59 * Supplies:
60 * hmsInitForward (map)
61 * load the forward transport maps
62 * hmsForward (uT,zT,u,z)
63 * make a single step transport calculation
64 jones 1.1 * (without treating the acceptance)
65 * hmsAccept (uT,zT,u,z)
66 * make a multi step transport calculation
67 * (also treating the acceptance)
68 *
69 * Note: - Before calling hmsForward or hmsAccept the forward
70 * transport maps have to be loaded by a call to hmsInitForward
71 *------------------------------------------------------------------------
72
73 SUBROUTINE hmsInitForward (filename,maxorder,path,p0)
74 IMPLICIT NONE
75 CHARACTER filename*(*)
76 INTEGER maxorder
77 LOGICAL path
78 REAL p0
79
80 * -- load the HMS forward (cosy) map
81 *
82 * Parameter:
83 * filename I : name of the forward map
84 * maxorder I : maximal order to take into account
85 jones 1.1 * path I : calculate the path length/TOF variation
86 * p0 I : momentum the spectrometer is set to
87 *
88 * reconInitFormat uses the same map-file as the simulation
89 * routines and digs out the necessary matrix elements to
90 * perform the focal plane corrections
91 *
92 * Map Line Ax Ax' Ay Ay' Al Adelta XxYyld
93 *
94 * Map : Map number 0-20
95 * 20: full HMS (target to focal plane)
96 * Line : Line number
97 * Ax : coeff's applied to x
98 * Ax' : coeff's applied to x'
99 * Ay : coeff's applied to y
100 * Ay' : coeff's applied to y'
101 * Al : coeff's applied to path length deviation
102 * X : exponent of the actual x
103 * x : exponent of the actual x'
104 * Y : exponent of the actual y
105 * y : exponent of the actual y'
106 jones 1.1 * l : exponent of the actual path length
107 * d : exponent of the delta
108
109 REAL fac,me
110 PARAMETER (fac = 0.9904599)
111 PARAMETER (me = 0.00051099906)
112
113 ! transport map common block
114 INTEGER NMAPS,NLINES
115 PARAMETER (NMAPS = 20)
116 PARAMETER (NLINES = 10000)
117
118 INTEGER first(NMAPS),last(NMAPS),order
119 INTEGER e1(NLINES),e2(NLINES),e3(NLINES)
120 INTEGER e4(NLINES),e5(NLINES),e6(NLINES)
121 REAL c1(NLINES),c2(NLINES),c3(NLINES)
122 REAL c4(NLINES),c5(NLINES)
123
124 COMMON /hmsTrack/first,last,order,
125 > e1,e2,e3,e4,e5,e6,c1,c2,c3,c4,c5
126
127 jones 1.1 REAL delta,mep0
128 COMMON /hmsTrackLabCOSY/delta,mep0
129
130 ! other variables
131 INTEGER i,m,n,num,eof
132
133 ! read the necessary coeffs from the map data file
134 10 FORMAT (1X,I2,1X,I5,1X,5(E14.7E2,1X),6I1)
135
136 mep0 = me/(p0*fac)
137
138 DO i=1,NMAPS
139 first(i) = 0
140 last (i) = 0
141 ENDDO
142
143 DO i=1,NLINES
144 c1(i) = 0.
145 c2(i) = 0.
146 c3(i) = 0.
147 c4(i) = 0.
148 jones 1.1 c5(i) = 0.
149 e1(i) = 0
150 e2(i) = 0
151 e3(i) = 0
152 e4(i) = 0
153 e5(i) = 0
154 e6(i) = 0
155 ENDDO
156
157 num = 1
158 order = 0
159
160 OPEN (UNIT=97,FILE=filename,STATUS='OLD')
161
162 READ (97,10,IOSTAT=eof) m,n,
163 > c1(num), c2(num), c3(num), c4(num), c5(num),
164 > e1(num), e2(num), e3(num), e4(num), e5(num), e6(num)
165 IF (.not. path) c5(num) = 0.
166
167 DO WHILE (eof .GE. 0)
168
169 jones 1.1 i = e1(num)+e2(num)+e3(num)+e4(num)+e5(num)+e6(num)
170 IF ((i .LE. maxorder) .AND. ((c1(num) .NE. 0.0) .OR.
171 > (c2(num) .NE. 0.0) .OR. (c3(num) .NE. 0.0) .OR.
172 > (c4(num) .NE. 0.0) .OR. (c5(num) .NE. 0.0))) THEN
173 IF (i .GT. order) order = i
174 IF (0 .EQ. first(m)) first(m) = num
175 last(m) = num
176 num = num+1
177 ENDIF
178 READ (97,10,IOSTAT=eof) m,n,
179 > c1(num), c2(num), c3(num), c4(num), c5(num),
180 > e1(num), e2(num), e3(num), e4(num), e5(num), e6(num)
181 IF (.not. path) c5(num) = 0.
182
183 ENDDO
184
185 CLOSE (97)
186
187 RETURN
188 END
189
190 jones 1.1
191
192 *------------------------------------------------------------------------
193
194 SUBROUTINE hmsApplyCOSY (u,num)
195 IMPLICIT NONE
196 REAL u(6)
197 INTEGER num
198
199 * -- apply a COSY matrix on the COSY vector u
200 * (used in the forward tracking only)
201 *
202 * Parameter
203 * u IO : coordinate vector (COSY)
204 * u(1,2) : x [m], A = out of plane coordinates (downwards)
205 * u(3,4) : y [m], B = inplane coordinates (perp. on x,z)
206 * u(5) : l [m] = [?] = path length deviation (?)
207 * u(6) : delta = relative deviation of the particle
208 * momentum from p0
209 * num I : cosy matrix to apply
210
211 jones 1.1 ! transport map common block
212 INTEGER NMAPS,NLINES
213 PARAMETER (NMAPS = 20)
214 PARAMETER (NLINES = 10000)
215
216 INTEGER first(NMAPS),last(NMAPS),order
217 INTEGER e1(NLINES),e2(NLINES),e3(NLINES)
218 INTEGER e4(NLINES),e5(NLINES),e6(NLINES)
219 REAL c1(NLINES),c2(NLINES),c3(NLINES)
220 REAL c4(NLINES),c5(NLINES)
221
222 COMMON /hmsTrack/first,last,order,
223 > e1,e2,e3,e4,e5,e6,c1,c2,c3,c4,c5
224
225 ! other variables
226 REAL a,uu1(0:10),uu2(0:10),uu3(0:10)
227 REAL uu4(0:10),uu5(0:10),uu6(0:10)
228 INTEGER i
229
230 ! calculate the powers of the focal plane coordinates
231 uu1(0) = 1.
232 jones 1.1 uu2(0) = 1.
233 uu3(0) = 1.
234 uu4(0) = 1.
235 uu5(0) = 1.
236 uu6(0) = 1.
237 uu1(1) = u(1)
238 uu2(1) = u(2)
239 uu3(1) = u(3)
240 uu4(1) = u(4)
241 uu5(1) = u(5)
242 uu6(1) = u(6)
243 DO i=2,order
244 uu1(i)=uu1(i-1)*uu1(1)
245 uu2(i)=uu2(i-1)*uu2(1)
246 uu3(i)=uu3(i-1)*uu3(1)
247 uu4(i)=uu4(i-1)*uu4(1)
248 uu5(i)=uu5(i-1)*uu5(1)
249 uu6(i)=uu6(i-1)*uu6(1)
250 ENDDO
251
252 DO i=1,5
253 jones 1.1 u(i) = 0.
254 ENDDO
255
256 ! apply the cosy matrix
257 DO i = first(num),last(num)
258 a = uu1(e1(i)) * uu2(e2(i)) * uu3(e3(i)) *
259 > uu4(e4(i)) * uu5(e5(i)) * uu6(e6(i))
260 u(1) = u(1) + c1(i) * a
261 u(2) = u(2) + c2(i) * a
262 u(3) = u(3) + c3(i) * a
263 u(4) = u(4) + c4(i) * a
264 u(5) = u(5) + c5(i) * a
265 ENDDO
266
267 RETURN
268 END
269
270 *------------------------------------------------------------------------
271
272 SUBROUTINE hmsForward (uT,zT,u,zz,lost)
273 IMPLICIT none
274 jones 1.1 REAL uT(6),zT,u(6),zz
275 LOGICAL lost
276
277 * -- make a single step transport calculation (without treating the acceptance)
278 *
279 * Parameter:
280 * uT,zT I : target coordinates
281 * uT(1,2) : x [m], dx/dz = out of plane coord. (downwards)
282 * uT(3,4) : y [m], dy/dz = inplane coord. (perp. on x,z)
283 * uT(5) : l [?] = path length deviation (?)
284 * uT(6) : delta; relative deviation of the particle
285 * momentum from p0
286 * zT : z [m]; in axis coordinate (towards HMS)
287 * u,zz O : focal plane coordinates
288 * u(1,2) : x [m], dx/dz; out of plane coord. (downwards)
289 * u(3,4) : y [m], dy/dz; inplane coord. (perp. on x,z)
290 * u(5) : l [?]; path length deviation (?)
291 * u(6) : delta; relative deviation of the particle
292 * momentum from p0
293 * zz : z [m]; in axis coordinate (towards HMS)
294 * position where the particle stops inside HMS
295 jones 1.1 * (here always z-calorimeter)
296 * lost O : set to .true. if the particle is lost in the HMS,
297 * otherwise .false. (here always false)
298
299 ! other variables
300 INTEGER i
301
302 lost = .FALSE.
303 zz = 26.44743
304
305 ! copy the target coordinates to the u vector
306 DO i=1,6
307 u(i)=uT(i)
308 ENDDO
309
310 ! drift backwards to z=0
311 u(1)=u(1)-u(2)*zT
312 u(3)=u(3)-u(4)*zT
313
314 !go up to the hut
315 CALL hmsApplyCOSY (u,20)
316 jones 1.1
317 RETURN
318 END
319
320 *------------------------------------------------------------------------
321
322 LOGICAL FUNCTION hmsCheckDipole (u,zz,dz)
323 IMPLICIT NONE
324 REAL u(6),zz,dz
325
326 * -- check for the dipole aperture
327 *
328 * Parameter
329 * u I : target coordinates (COSY or LAB)
330 * u(1,2) : x [m], A or dx/dz = out of plane coord. (downwards)
331 * u(3,4) : y [m], A or dy/dz = inplane coord. (perp. on x,z)
332 * u(5) : l [?] = path length deviation (?)
333 * u(6) : delta; relative deviation of the particle
334 * momentum from p0
335 * zz IO : z [m]; in axis coordinate (towards HMS)
336 * dz I : distance from last check point [m]
337 jones 1.1
338 zz = zz+dz
339
340 hmsCheckDipole = .TRUE.
341
342 IF (ABS(u(1)).GT.0.27940.OR.ABS(u(3)).GT.0.18415) THEN
343 IF (ABS(u(1)).GT.0.34290.OR.ABS(u(3)).GT.0.20320) RETURN
344 IF (ABS(u(1)).GT.0.27940.AND.ABS(u(3)).GT.0.12065) THEN
345 IF (((ABS(u(1))-0.27940)**2 +
346 > (ABS(u(3))-0.12065)**2).GT.0.06350**2) RETURN
347 ENDIF
348 IF (ABS(u(1)).GT.0.13970 .OR.
349 > (ABS(u(1))-10.1852*ABS(u(3))).GT.2.069633) RETURN
350 ENDIF
351
352 hmsCheckDipole = .FALSE.
353
354 RETURN
355 END
356
357 *------------------------------------------------------------------------
358 jones 1.1
359 LOGICAL FUNCTION hmsCheckQuad (u,zz,dz,r)
360 IMPLICIT NONE
361 REAL u(6),zz,dz,r
362
363 * -- check for the quadrupole aperture
364 *
365 * Parameter
366 * u I : target coordinates (COSY or LAB)
367 * u(1,2) : x [m], A or dx/dz = out of plane coord. (downwards)
368 * u(3,4) : y [m], A or dy/dz = inplane coord. (perp. on x,z)
369 * u(5) : l [?] = path length deviation (?)
370 * u(6) : delta; relative deviation of the particle
371 * momentum from p0
372 * zz IO : z [m]; in axis coordinate (towards HMS)
373 * dz I : distance from last check point [m]
374 * r I : aperture radius [m]
375
376 zz = zz+dz
377 hmsCheckQuad = ((u(1)**2+u(3)**2) .GT. r**2)
378
379 jones 1.1 RETURN
380 END
381
382 *------------------------------------------------------------------------
383
384 LOGICAL FUNCTION hmsDriftOcta (u,zz,dz,x,y,m,b)
385 IMPLICIT NONE
386 REAL u(6),zz,dz,x,y,m,b
387
388 * -- drift electron and check for the octagon
389 *
390 * Parameter
391 * u I : target coordinates (COSY or LAB)
392 * u(1,2) : x [m], A or dx/dz = out of plane coord. (downwards)
393 * u(3,4) : y [m], A or dy/dz = inplane coord. (perp. on x,z)
394 * u(5) : l [?] = path length deviation (?)
395 * u(6) : delta; relative deviation of the particle
396 * momentum from p0
397 * zz IO : z [m]; in axis coordinate (towards HMS)
398 * dz I : distance from last check point [m]
399 * x,y,m,b I : octagon coordinates (x,y,b [m], m [1])
400 jones 1.1
401 zz=zz+dz
402 u(1)=u(1)+u(2)*dz
403 u(3)=u(3)+u(4)*dz
404 hmsDriftOcta = (ABS(u(1)) .GT. x) .OR. (ABS(u(3)) .GT. y) .OR.
405 > ((ABS(u(1))+m*ABS(u(3))) .GT. b)
406
407 RETURN
408 END
409
410 *------------------------------------------------------------------------
411
412 LOGICAL FUNCTION hmsDriftTPlate (u,zz,dz,r,y)
413 IMPLICIT NONE
414 REAL u(6),zz,dz,r,y
415
416 * -- drift electron and check for the dipole transition plate
417 *
418 * Parameter
419 * u IO : coordinate vector
420 * zz IO : z position [m]
421 jones 1.1 * dz I : drift distance
422 * r,y I : aperture
423
424 zz=zz+dz
425 u(1)=u(1)+u(2)*dz
426 u(3)=u(3)+u(4)*dz
427 hmsDriftTPlate = ((u(1)**2+u(3)**2) .GT. r**2) .OR.
428 > (ABS(u(3)) .GT. y)
429
430 RETURN
431 END
432
433 *------------------------------------------------------------------------
434
435 LOGICAL FUNCTION hmsDriftCirc (u,zz,dz,r)
436 IMPLICIT NONE
437 REAL u(6),zz,dz,r
438
439 * -- drift electron and check for circular aperture
440 *
441 * Parameter
442 jones 1.1 * u I : target coordinates (COSY or LAB)
443 * u(1,2) : x [m], A or dx/dz = out of plane coord. (downwards)
444 * u(3,4) : y [m], A or dy/dz = inplane coord. (perp. on x,z)
445 * u(5) : l [?] = path length deviation (?)
446 * u(6) : delta; relative deviation of the particle
447 * momentum from p0
448 * zz IO : z [m]; in axis coordinate (towards HMS)
449 * dz I : distance from last check point [m]
450 * r,y I : transition plate aperture radius [m]
451
452 zz=zz+dz
453 u(1)=u(1)+u(2)*dz
454 u(3)=u(3)+u(4)*dz
455 hmsDriftCirc = ((u(1)**2+u(3)**2) .GT. r**2)
456
457 RETURN
458 END
459
460 *------------------------------------------------------------------------
461
462 LOGICAL FUNCTION hmsDriftRect (u,zz,dz,x0,x,y0,y)
463 jones 1.1 IMPLICIT NONE
464 REAL u(6),zz,dz,x0,x,y0,y
465
466 * -- drift electron and check for rectangular aperture
467 *
468 * Parameter
469 * u I : target coordinates (COSY or LAB)
470 * u(1,2) : x [m], A or dx/dz = out of plane coord. (downwards)
471 * u(3,4) : y [m], A or dy/dz = inplane coord. (perp. on x,z)
472 * u(5) : l [?] = path length deviation (?)
473 * u(6) : delta; relative deviation of the particle
474 * momentum from p0
475 * zz IO : z [m]; in axis coordinate (towards HMS)
476 * dz I : distance from last check point [m]
477 * x,x0,
478 * y,y0 I : rectangular aperture [m] (half size and offset)
479
480 zz=zz+dz
481 u(1)=u(1)+u(2)*dz
482 u(3)=u(3)+u(4)*dz
483 hmsDriftRect = (ABS(u(1)-x0) .GT. x) .OR.
484 jones 1.1 > (ABS(u(3)-y0) .GT. y)
485 RETURN
486 END
487
488 *------------------------------------------------------------------------
489
490 SUBROUTINE hmsAccept (uT,zT,u,zz,lost)
491 IMPLICIT none
492 REAL uT(6),zT,u(6),zz
493 LOGICAL lost
494
495 * -- make a transport calculation to find the acceptance
496 *
497 * Parameter:
498 * uT,zT I : target coordinates
499 * uT(1,2) : x [m], dx/dz; out of plane coord. (downwards)
500 * uT(3,4) : y [m], dy/dz; inplane coord. (perp. on x,z)
501 * uT(5) : l [?]; path length deviation (?)
502 * uT(6) : delta; relative deviation of the particle
503 * momentum from p0
504 * zT : z [m]; in axis coordinate (towards HMS)
505 jones 1.1 * u,zz O : focal plane coordinates
506 * u(1,2) : x [m], dx/dz; out of plane coord. (downwards)
507 * u(3,4) : y [m], dy/dz; inplane coord. (perp. on x,z)
508 * u(5) : l [?]; path length deviation (?)
509 * u(6) : delta; relative deviation of the particle
510 * momentum from p0
511 * zz : z [m]; in axis coordinate (towards HMS)
512 * position where the particle stops inside HMS
513 * lost O : set to .true. if the particle is lost in the HMS,
514 * otherwise .false.
515
516 LOGICAL hmsCheckQuad
517 LOGICAL hmsCheckDipole
518 LOGICAL hmsDriftOcta
519 LOGICAL hmsDriftTPlate
520 LOGICAL hmsDriftCirc
521 LOGICAL hmsDriftRect
522
523 REAL uS(6)
524 INTEGER i
525
526 jones 1.1 lost = .TRUE.
527 zz = zT
528
529 ! copy the target coordinates to the u vector
530 DO i=1,6
531 u(i)=uT(i)
532 ENDDO
533
534 ! ----------------------------------------------------------- sive slit
535 ! drift to sieve slit: z=1.27636m,
536 ! octagon edge m=2.545977, b=0.13502640m
537 IF (hmsDriftOcta(u,zz,1.27636-zT,0.0900176,0.0353568,
538 > 2.5459770,0.1350264)) RETURN
539
540 ! drift to back of sieve slit: z=1.33986m, dz=0.0635m
541 ! octagon edge m=2.546569, b=0.141655189 m
542 IF (hmsDriftOcta(u,zz,0.0635,0.0944368,0.0370839,
543 > 2.546569,0.141655189)) RETURN
544
545 ! ------------------------------------------------------------------ Q1
546 ! drift to mechanical entrance of Q1: z=1.4960m, dz=0.15614m
547 jones 1.1 IF (hmsDriftCirc(u,zz,0.15614,0.202575)) RETURN
548
549 ! drift to Q1 entrance EFB: z=1.775805635m, dz=0.279805635m
550 u(1)=u(1)+u(2)*0.279805635
551 u(3)=u(3)+u(4)*0.279805635
552
553 ! and save values
554 DO i=1,6
555 uS(i)=u(i)
556 ENDDO
557
558 ! transport through Q1 fringe fields:
559 CALL hmsApplyCOSY (u,1)
560 IF (hmsCheckQuad(u,zz,0.279805635,0.202575)) RETURN
561
562 ! transport through Q1, 1/5 at a time:
563 DO i=1,4
564 CALL hmsApplyCOSY (u,2)
565 IF (hmsCheckQuad(u,zz,1.87838873/5.,0.202575)) RETURN
566 ENDDO
567
568 jones 1.1 ! restore values:
569 DO i=1,6
570 u(i)=uS(i)
571 ENDDO
572
573 ! transport to Q1 exit EFB:
574 CALL hmsApplyCOSY (u,3)
575 IF (hmsCheckQuad(u,zz,1.87838873/5.,0.202575)) RETURN
576
577 ! drift to Q1 mechanical exit: z=3.9340m, dz=0.279805635m
578 IF (hmsDriftCirc(u,zz,0.279805635,0.202575)) RETURN
579
580 ! ------------------------------------------------------------------ Q2
581 ! drift to mechanical entrance of Q2: z=4.5610m, dz=0.6270m
582 IF (hmsDriftCirc(u,zz,0.6270,0.29840)) RETURN
583
584 ! drift to Q2 entrance EFB: z=4.887021890m, dz=0.326021890m
585 u(1)=u(1)+u(2)*0.326021890
586 u(3)=u(3)+u(4)*0.326021890
587
588 ! and save values
589 jones 1.1 DO i=1,6
590 uS(i)=u(i)
591 ENDDO
592
593 ! transport through Q2 fringe fields:
594 CALL hmsApplyCOSY (u,4)
595 IF (hmsCheckQuad(u,zz,0.326021890,0.29840)) RETURN
596
597 ! transport through Q2, 1/5 at a time:
598 DO i=1,4
599 CALL hmsApplyCOSY (u,5)
600 IF (hmsCheckQuad(u,zz,2.15595622/5.,0.29840)) RETURN
601 ENDDO
602
603 ! restore values:
604 DO i=1,6
605 u(i)=uS(i)
606 ENDDO
607
608 ! transport to Q2 exit EFB:
609 CALL hmsApplyCOSY (u,6)
610 jones 1.1 IF (hmsCheckQuad(u,zz,2.15595622/5.,0.29840)) RETURN
611
612 ! drift to Q2 mechanical exit: z=7.3690m, dz=0.326021890m
613 IF (hmsDriftCirc(u,zz,0.326021890,0.29840)) RETURN
614
615 ! ------------------------------------------------------------------ Q3
616 ! drift to mechanical entrance of Q3: z=7.6610m, dz=0.2920m
617 IF (hmsDriftCirc(u,zz,0.2920,0.29840)) RETURN
618
619 ! drift to Q3 entrance EFB: z=7.990200290m, dz=0.329200290m
620 u(1)=u(1)+u(2)*0.329200290
621 u(3)=u(3)+u(4)*0.329200290
622
623 ! save values
624 DO i=1,6
625 uS(i)=u(i)
626 ENDDO
627
628 ! transport through Q3 fringe fields:
629 CALL hmsApplyCOSY (u,7)
630 IF (hmsCheckQuad(u,zz,0.329200290,0.29840)) RETURN
631 jones 1.1
632 DO i=1,4
633 CALL hmsApplyCOSY (u,8)
634 IF (hmsCheckQuad(u,zz,2.14959942/5.,0.29840)) RETURN
635 ENDDO
636
637 ! and restore values:
638 DO i=1,6
639 u(i)=uS(i)
640 ENDDO
641
642 ! transport to Q3 exit EFB:
643 CALL hmsApplyCOSY (u,9)
644 IF (hmsCheckQuad(u,zz,2.14959942/5.,0.29840)) RETURN
645
646 ! drift to Q3 mechanical exit: z=10.4690m, dz=0.329200290m
647 IF (hmsDriftCirc(u,zz,0.329200290,0.29840)) RETURN
648
649 ! -------------------------------------------------------------- Dipole
650 ! drift to transition plate: z=11.058002m, dz=0.589002m
651 IF (hmsDriftTPlate(u,zz,0.589002,0.30480,0.205232)) RETURN
652 jones 1.1
653 ! drift to opposite side of transition plate: z=11.092800m, dz=0.034798m
654 IF (hmsDriftTPlate(u,zz,0.034798,0.30480,0.205232)) RETURN
655 IF (hmsCheckDipole(u,zz,0.)) RETURN
656
657 ! drift to D magnetic entrance: z=11.55m, dz=0.4572m
658 u(1)=u(1)+u(2)*0.4572
659 u(3)=u(3)+u(4)*0.4572
660 IF (hmsCheckDipole(u,zz,0.4572)) RETURN
661
662 ! save values:
663 DO i=1,6
664 uS(i)=u(i)
665 ENDDO
666
667 ! transport through 1/5 D with rotated entrance face:
668 CALL hmsApplyCOSY (u,10)
669 IF (hmsCheckDipole(u,zz,5.26053145/5.)) RETURN
670
671 ! transport through 3/5 D with sector segments:
672 DO i=1,3
673 jones 1.1 CALL hmsApplyCOSY (u,11)
674 IF (hmsCheckDipole(u,zz,5.26053145/5.)) RETURN
675 ENDDO
676
677 ! restore values:
678 DO i=1,6
679 u(i)=uS(i)
680 ENDDO
681
682 ! transport through D (entrance to exit, fringe fields included):
683 CALL hmsApplyCOSY (u,13)
684 IF (hmsCheckDipole(u,zz,5.26053145/5.)) RETURN
685
686 ! drift to transition plate: z=0.4572m, dz=0.4572m
687 IF (hmsDriftTPlate(u,zz,0.457200,0.34290,0.205232)) RETURN
688 IF (hmsCheckDipole(u,zz,0.)) RETURN
689
690 ! drift to opposite side of transition plate: z=0.491998m,dz=0.034798m
691 IF (hmsDriftTPlate(u,zz,0.034798,0.34290,0.205232)) RETURN
692
693 ! drift to end of first piece of telescope: z=1.119378m, dz=0.62738m
694 jones 1.1 IF (hmsDriftCirc(u,zz,0.62738,0.338450)) RETURN
695
696 ! drift to end of second piece of telescope: z=4.086098m, dz=2.96672m
697 IF (hmsDriftCirc(u,zz,2.96672,0.384175)) RETURN
698
699 ! drift to end of third piece of telescope: z=5.578398m, dz=1.4923m
700 IF (hmsDriftCirc(u,zz,1.49230,0.460375)) RETURN
701
702 ! ----------------------------------------------------------------- hut
703 ! drift to focal plane. This is the reference point for detector positions.
704 u(1)=u(1)+u(2)*0.671602
705 u(3)=u(3)+u(4)*0.671602
706 zz = zz + 0.671602
707 DO i=1,6
708 uS(i)=u(i)
709 ENDDO
710
711 ! drift to DC1 entrance: z=-0.51923-0.036=-0.55523m, dz=-0.55523m
712 IF(hmsDriftRect(u,zz,-0.55523,-0.01670,0.565,-0.00343,0.26))RETURN
713 ! drift to DC1 exit: z=-0.51923+0.054=-0.46523m, dz=0.090m
714 IF(hmsDriftRect(u,zz, 0.09000,-0.01670,0.565,-0.00343,0.26))RETURN
715 jones 1.1 ! drift to DC2 entrance: z=0.29299-0.036=0.25699m, dz=0.72222m
716 IF(hmsDriftRect(u,zz, 0.72222,-0.02758,0.565,-0.01653,0.26))RETURN
717 ! drift to DC2 exit: z=0.29299+0.054=0.34699m, dz=0.090m
718 IF(hmsDriftRect(u,zz, 0.09000,-0.02758,0.565,-0.01653,0.26))RETURN
719
720 ! drift to S1X: z=0.7783m, dz=0.43131m
721 IF (hmsDriftRect(u,zz,0.43131,0.015,0.6025,0.000,0.3775)) RETURN
722 ! drift to S1Y: z=0.9752m, dz=0.1969m
723 IF (hmsDriftRect(u,zz,0.19690,0.000,0.6025,0.001,0.3775)) RETURN
724 ! skip CK - no survey information
725 ! drift to S2X: z=2.9882m, dz=2.013m
726 IF (hmsDriftRect(u,zz,2.01300,0.004,0.6025,0.000,0.3775)) RETURN
727 ! drift to S2Y: z=3.1851m, dz=0.1969m
728 IF (hmsDriftRect(u,zz,0.19690,0.000,0.6025,0.013,0.3775)) RETURN
729
730 ! drift to CAL: z=3.3869m, dz=0.2018m
731 IF (hmsDriftRect(u,zz,0.20180,-0.134,0.6000,0.000,0.3000)) RETURN
732
733 ! --------------------------------------------------------------- done
734 lost = .FALSE.
735 DO i=1,6
736 jones 1.1 u(i)=uS(i)
737 ENDDO
738
739 RETURN
740 END
741
742 *------------------------------------------------------------------------
743 *------------------------------------------------------------------------
744 *------------------------------------------------------------------------
745 *
746 * PART 2: HMS Reconstruction (Backward Tracking; Focal Plane to Target)
747 * -=======-
748 *
749 * Reconstruction (backward tracking) in the Jlab HMS hall C
750 * spectrometer using reconstruction and forward COSY matrices
751 * (including the effects of beam offsets (out-of plane))
752 *
753 * Both the normal in-plane scattering and the more special
754 * out-of-plane scattering are handeled. The later makes use
755 * of the forward COSY matrices. The algorithm was tested for
756 * beam offsets in the range of cm (up or below the
757 jones 1.1 * nominal scattering plane)
758 *
759 * Supplies:
760 * hmsInitRecon (map,p0)
761 * load the reconstruction maps
762 * hmsInPlane (u,uT,ok)
763 * reconstruction of the target coordinates
764 * (delta, dx/dz, y, dy/dz) at z=0
765 * hmsOutOfPlane (u,x,uT,ok)
766 * reconstruction of the target coordinates
767 * (delta, dx/dz, y, dy/dz) at z=0 including the
768 * vertical beam offset
769 *
770 * Note: - Before calling hmsReconInPlane or hmsReconOutOfPlaneAccept
771 * the reconstruction map has to be loaded by a call to
772 * hmsInitRecon
773 * - Before calling hmsReconOutOfPlane the forward transport
774 * maps have to be loaded by a call to hmsInitForward
775 *------------------------------------------------------------------------
776
777 SUBROUTINE hmsInitRecon (filename)
778 jones 1.1 IMPLICIT NONE
779 CHARACTER filename*(*)
780
781 * -- load the HMS reconstruction map
782 *
783 * Parameter:
784 * filename I : name of the reconstruction map
785 * p0 I : momentum the spectrometer is set to
786 *
787 * File format:
788 * Ax' Ay Ay' Ad XxYy
789 *
790 * Ax' : coeff's giving target x'
791 * Ay : coeff's giving target y
792 * Ay' : coeff's giving target y'
793 * Ad : coeff's giving the delta
794 * X : exponent of the focal plane x
795 * x : exponent of the focal plane x'
796 * Y : exponent of the focal plane y
797 * y : exponent of the focal plane y'
798
799 jones 1.1 ! matrix elements for the reconstruction
800 INTEGER NTERMS
801 PARAMETER (NTERMS=1000)
802
803 INTEGER e1(NTERMS),e2(NTERMS),e3(NTERMS),e4(NTERMS),num,order
804 REAL c1(NTERMS),c2(NTERMS),c3(NTERMS),c4(NTERMS)
805 COMMON /hmsRecon/num,order,e1,e2,e3,e4,c1,c2,c3,c4
806
807 ! other variables
808 CHARACTER line*256
809 INTEGER i,eof,e5
810
811 ! read the necessary coeffs from the map data file
812 10 FORMAT (A)
813 20 FORMAT (1X,4G16.9,1X,5I1)
814
815 OPEN (UNIT=97,FILE=filename,STATUS='OLD')
816
817 num = 1
818 order = 0
819
820 jones 1.1 READ (97,10,IOSTAT=eof) line
821
822 DO WHILE (eof .EQ. 0)
823 IF ((line(1:1) .NE. '!') .AND.
824 > (line(1:4) .NE. ' ---') .AND.
825 > (line(1:2) .NE. 'h_')) THEN
826 ! read the reconstruction coefficents
827 READ (line,20,IOSTAT=eof)
828 > c2(num), c3(num), c4(num), c1(num),
829 > e1(num), e2(num), e3(num), e4(num), e5
830 IF (((c1(num) .NE. 0.0) .OR. (c2(num) .NE. 0.0) .OR.
831 > (c3(num) .NE. 0.0) .OR. (c4(num) .NE. 0.0)) .AND.
832 > (e5 .EQ. 0) .AND. (eof .EQ. 0)) THEN
833 i = e1(num) + e2(num) + e3(num) + e4(num)
834 IF (i .GT. order) order = i
835 num = num+1
836 ENDIF
837 ENDIF
838
839 IF (eof .EQ. 0) READ (97,10,IOSTAT=eof) line
840 ENDDO
841 jones 1.1 num = num-1
842
843 CLOSE (97)
844
845 RETURN
846 END
847
848 *------------------------------------------------------------------------
849
850 SUBROUTINE hmsReconOffset (uT,du)
851
852 IMPLICIT NONE
853
854 REAL uT(6), du(4)
855
856 * -- calculates the focal plane correction for a given
857 * set of target coordinates
858 *
859 * Parameter
860 * uT I : target coordinates (lab)
861 * uT(1,2) : x [m], dx/dz = out of plane coord. (downwards)
862 jones 1.1 * uT(3,4) : y [m], dy/dz = inplane coord. (perp. on x,z)
863 * uT(5) : z [m] = in axis coordinate (towards HMS)
864 * uT(6) : delta = relative deviation of the particle
865 * momentum from p0
866 * du O : correction values for the focal plane quantities
867 * du(1,2) : correction in x [m], dx/dz
868 * du(3,4) : correction in y [m], dy/dz
869 *
870 * frw 9/2000
871 * adjusted and repaired code in course of migration to g77
872 *
873 * some notes on how this routine works, since it does not obviously
874 * match the description given on MM's web page:
875 *
876 * Based on the web page documentation, this subroutine is supposed to transform
877 * the current guess at the reconstructed target coordinates back to the focal
878 * plane twice, once with x_target = beam value and once for x_target=0 and
879 * take the difference between the results as the correction du
880 * This code takes advantage of numerous cancellation and thereby reduces the
881 * calculation effort:
882 * the sum below (DO i = first(20), last(20) ...) would be evaluated for both
883 jones 1.1 * values, uT(1)=x and uT(1)=0. The difference would be du.
884 * if e1(i)=0 then the result is independent of the value of uT(1) and since
885 * that's the only difference between the 2 evaluations, these terms cancel
886 * in the difference giving du, thus they are omitted.
887 * similarly, if e1(i)<>0, the evaluation of the terms in the sum for uT(1)=0
888 * would result in all terms being 0, thus subtracting them from the sum for
889 * uT(1)=x changes nothing. Therefore, the difference between the transform
890 * for uT(1)=x and uT(1)=0 is simply equal to the transform for uT(1)=x if the
891 * elements for which e1(i)=0 are skipped.
892
893 ! transport map common block
894 INTEGER NMAPS,NLINES
895 PARAMETER (NMAPS = 20)
896 PARAMETER (NLINES = 10000)
897
898 INTEGER first(NMAPS),last(NMAPS),order
899 INTEGER e1(NLINES),e2(NLINES),e3(NLINES)
900 INTEGER e4(NLINES),e5(NLINES),e6(NLINES)
901 REAL c1(NLINES),c2(NLINES),c3(NLINES)
902 REAL c4(NLINES),c5(NLINES)
903
904 jones 1.1 COMMON /hmsTrack/first,last,order,
905 > e1,e2,e3,e4,e5,e6,c1,c2,c3,c4,c5
906
907 REAL delta,mep0
908 COMMON /hmsTrackLabCOSY/delta,mep0
909
910 ! other variables
911 REAL a,uu1(0:10),uu2(0:10),uu3(0:10),uu4(0:10),uu6(0:10)
912 INTEGER i
913
914 uu1(0) = 1.
915 uu2(0) = 1.
916 uu3(0) = 1.
917 uu4(0) = 1.
918 uu6(0) = 1.
919
920 ! drift backwards to z=0
921 uu1(1)=uT(1)-uT(2)*uT(5)
922 uu3(1)=uT(3)-uT(4)*uT(5)
923
924 uu2(1) = uT(2)
925 jones 1.1 uu4(1) = uT(4)
926 uu6(1) = uT(6)
927
928
929 ! calculate the powers of the COSY coordinates
930 DO i=1,4
931 du(i) = 0.
932 ENDDO
933
934 DO i=2,order
935 uu1(i)=uu1(i-1)*uu1(1)
936 uu2(i)=uu2(i-1)*uu2(1)
937 uu3(i)=uu3(i-1)*uu3(1)
938 uu4(i)=uu4(i-1)*uu4(1)
939 uu6(i)=uu6(i-1)*uu6(1)
940 ENDDO
941
942 ! calculate the focal plane offsets
943
944 DO i = first(20), last(20)
945
946 jones 1.1 IF (e1(i) .NE. 0) THEN
947
948 a = uu1(e1(i)) * uu2(e2(i)) * uu3(e3(i))
949 > * uu4(e4(i)) * uu6(e6(i))
950
951 du(1) = du(1) + c1(i) * a
952 du(2) = du(2) + c2(i) * a
953 du(3) = du(3) + c3(i) * a
954 du(4) = du(4) + c4(i) * a
955
956 ENDIF
957
958 ENDDO
959
960 RETURN
961 END
962
963 *------------------------------------------------------------------------
964
965 SUBROUTINE hmsReconInPlane (u,uT,ok)
966 IMPLICIT NONE
967 jones 1.1 REAL u(4),uT(6)
968 LOGICAL ok
969
970 * -- performs the reconstruction of the target coordinates
971 * (delta, dx/dz, y, dy/dz) at z=0
972 *
973 * Parameter:
974 * u I : focal plane coordinates (lab)
975 * u(1,2) : x [m], dx/dz = out of plane coord. (downwards)
976 * u(3,4) : y [m], dy/dz = inplane coord. (perp. on x,z)
977 * uT O : target coordinates (lab)
978 * uT(1,2) : x [m], dx/dz = out of plane coord. (downwards)
979 * uT(3,4) : y [m], dy/dz = inplane coord. (perp. on x,z)
980 * uT(5) : z [m] = in axis coordinate (towards HMS)
981 * uT(6) : delta (relative deviation of the particle
982 * momentum from p0)
983 * ok IO : status variable
984 * - if false no action is taken
985 * - set to false when no reconstruction is found
986
987 ! matrix elemnts needed for calculating the focal plane offset
988 jones 1.1 INTEGER NTERMS
989 PARAMETER (NTERMS=1000)
990
991 INTEGER e1(NTERMS),e2(NTERMS),e3(NTERMS),e4(NTERMS),num,order
992 REAL c1(NTERMS),c2(NTERMS),c3(NTERMS),c4(NTERMS)
993 COMMON /hmsRecon/num,order,e1,e2,e3,e4,c1,c2,c3,c4
994
995 ! other variables
996 REAL a,uu1(0:10),uu2(0:10),uu3(0:10),uu4(0:10)
997 INTEGER i
998
999
1000 DO i=1,6
1001 uT(i) = 0.
1002 ENDDO
1003
1004 ! calculate the powers of the focal plane coordinates
1005
1006 uu1(0) = 1.
1007 uu2(0) = 1.
1008 uu3(0) = 1.
1009 jones 1.1 uu4(0) = 1.
1010
1011 uu1(1) = u(1)
1012 uu2(1) = u(2)
1013 uu3(1) = u(3)
1014 uu4(1) = u(4)
1015
1016 DO i=2,order
1017 uu1(i)=uu1(i-1)*uu1(1)
1018 uu2(i)=uu2(i-1)*uu2(1)
1019 uu3(i)=uu3(i-1)*uu3(1)
1020 uu4(i)=uu4(i-1)*uu4(1)
1021 ENDDO
1022
1023 ! calculate the target coordinates
1024 DO i = 1, num
1025 a = uu1(e1(i)) * uu2(e2(i)) * uu3(e3(i)) * uu4(e4(i))
1026
1027 uT(6) = uT(6) + c1(i) * a
1028 uT(2) = uT(2) + c2(i) * a
1029 uT(3) = uT(3) + c3(i) * a
1030 jones 1.1 uT(4) = uT(4) + c4(i) * a
1031 ENDDO
1032
1033 ok = ((ABS(uT(2)) .LT. 1.) .AND. (ABS(uT(3)) .LT. 1.) .AND.
1034 > (ABS(uT(4)) .LT. 1.) .AND. (ABS(uT(6)) .LT. 1.))
1035
1036 RETURN
1037 END
1038
1039 *------------------------------------------------------------------------
1040
1041 SUBROUTINE hmsReconOutOfPlane (u,x,uT,ok)
1042 IMPLICIT NONE
1043 REAL u(4),x,uT(6)
1044 LOGICAL ok
1045
1046 * -- performs the reconstruction of the target coordinates
1047 * (delta, dx/dz, y, dy/dz) for the out of plane case
1048 *
1049 * Parameter:
1050 * u I : focal plane coordinates (lab)
1051 jones 1.1 * u(1,2) : x [m], dx/dz = out of plane coord. (downwards)
1052 * u(3,4) : y [m], dy/dz = inplane coord. (perp. on x,z)
1053 * x I : x-offset at the target [m] (z=0) (out of plane; downwards)
1054 * uT O : target coordinates (lab)
1055 * uT(1,2) : x [m], dx/dz = out of plane coord. (downwards)
1056 * uT(3,4) : y [m], dy/dz = inplane coord. (perp. on x,z)
1057 * uT(5) : z [m] = in axis coordinate (towards HMS)
1058 * uT(6) : delta (relative deviation of the particle
1059 * ok IO : status variable
1060 * - if false no action is taken
1061 * - set to false when no reconstruction is found
1062 * momentum from p0)
1063
1064 REAL eps
1065 PARAMETER (eps = 0.0005) ! accuracy in delta
1066
1067 REAL dd,du(4),u0(4)
1068 INTEGER n
1069
1070
1071
1072 jones 1.1 CALL hmsReconInPlane (u,uT,ok) ! first guess
1073 uT(1) = x
1074
1075 du(1) = 0.
1076 du(2) = 0.
1077 du(3) = 0.
1078 du(4) = 0.
1079
1080 dd = 1.
1081 n = 0
1082
1083 DO WHILE ((ABS(dd) .GT. eps) .AND. (n .LT. 10) .AND. ok)
1084
1085 CALL hmsReconOffset (uT,du)
1086
1087 u0(1) = u(1)-du(1)
1088 u0(2) = u(2)-du(2)
1089 u0(3) = u(3)-du(3)
1090 u0(4) = u(4)-du(4)
1091
1092 dd = uT(6)
1093 jones 1.1
1094 CALL hmsReconInPlane (u0,uT,ok)
1095 uT(1) = x
1096
1097 dd = dd-uT(6)
1098 n = n+1
1099
1100 ENDDO
1101
1102 IF (ABS(dd) .GT. eps) ok = .FALSE. ! not converged
1103
1104 RETURN
1105 END
1106
|
Return to
hms_track.f
CVS log
Up to [HallC] / pol_hms_single