*------------------------------------------------------------------------ * * HMS_TRACK HMS Tracking Routines * -=========- * * Forward and Backward Tracking of electrons in the Jlab HMS hall * C spectrometer * * Note: - the HMS routines use a lab (HMS) coord. system * and the corresponding COSY coord. system, both * right handed with * x : pointing downwards * y : perpendicular to x,z, * pointing to the left (if seen in z-direction) * z : HMS axis, pointing from the target to the focal plane * * - all lengths (x,y,z,l,...) are measured in [m] * - all angles are measured as dx/dz,dy/dz (lab coords.) * or as A,B (COSY coords.) * - the momentum is measured in delta (relative momentum * deviation = 1-p0/pHMS) * * PART 1: Forward trackinh using COSY transport matrices * * * PART 2: Reconstruction (backward tracking) using reconstruction * and COSY transport matrices (including the * effects of a vertical beam offset (out-of plane)) * * * written by Markus Muehlbauer for the GEN Experiment *------------------------------------------------------------------------ *------------------------------------------------------------------------ * * PART 1: HMS Forward Tracking (Target to Focal Plane) * -=======- * * Forward tracking in the Jlab HMS hall C spectrometer * using COSY transport matrices * * developed by Cris Cothran * modified by Markus Muehlbauer * - CCs orignal program converted into subroutines * - mad additions for pure tracking, without checking the acceptance * - and changed the innermost loops applying the matrix * (which speeds up the whole thing by a factor of about 30) * * Supplies: * hmsInitForward (map) * load the forward transport maps * hmsForward (uT,zT,u,z) * make a single step transport calculation * (without treating the acceptance) * hmsAccept (uT,zT,u,z) * make a multi step transport calculation * (also treating the acceptance) * * Note: - Before calling hmsForward or hmsAccept the forward * transport maps have to be loaded by a call to hmsInitForward *------------------------------------------------------------------------ SUBROUTINE hmsInitForward (filename,maxorder,path,p0) IMPLICIT NONE CHARACTER filename*(*) INTEGER maxorder LOGICAL path REAL p0 * -- load the HMS forward (cosy) map * * Parameter: * filename I : name of the forward map * maxorder I : maximal order to take into account * path I : calculate the path length/TOF variation * p0 I : momentum the spectrometer is set to * * reconInitFormat uses the same map-file as the simulation * routines and digs out the necessary matrix elements to * perform the focal plane corrections * * Map Line Ax Ax' Ay Ay' Al Adelta XxYyld * * Map : Map number 0-20 * 20: full HMS (target to focal plane) * Line : Line number * Ax : coeff's applied to x * Ax' : coeff's applied to x' * Ay : coeff's applied to y * Ay' : coeff's applied to y' * Al : coeff's applied to path length deviation * X : exponent of the actual x * x : exponent of the actual x' * Y : exponent of the actual y * y : exponent of the actual y' * l : exponent of the actual path length * d : exponent of the delta * REAL fac,me PARAMETER (fac = 0.9904599) PARAMETER (me = 0.00051099906) ! transport map common block INTEGER NMAPS,NLINES PARAMETER (NMAPS = 20) PARAMETER (NLINES = 10000) INTEGER first(NMAPS),last(NMAPS),order INTEGER e1(NLINES),e2(NLINES),e3(NLINES) INTEGER e4(NLINES),e5(NLINES),e6(NLINES) REAL c1(NLINES),c2(NLINES),c3(NLINES) REAL c4(NLINES),c5(NLINES) COMMON /hmsTrack/first,last,order, > e1,e2,e3,e4,e5,e6,c1,c2,c3,c4,c5 REAL delta,mep0 COMMON /hmsTrackLabCOSY/delta,mep0 ! other variables INTEGER i,m,n,num,eof ! read the necessary coeffs from the map data file 10 FORMAT (X,I2,X,I5,X,5(E14.7E2,X),6I1) mep0 = me/(p0*fac) DO i=1,NMAPS first(i) = 0 last (i) = 0 ENDDO OPEN (UNIT=97,FILE=filename,STATUS='OLD') num = 1 order = 0 READ (97,10,IOSTAT=eof) m,n, > c1(num), c2(num), c3(num), c4(num), c5(num), > e1(num), e2(num), e3(num), e4(num), e5(num), e6(num) IF (.not. path) c5(num) = 0. DO WHILE (eof .GE. 0) i = e1(num)+e2(num)+e3(num)+e4(num)+e5(num)+e6(num) IF ((i .LE. maxorder) .AND. ((c1(num) .NE. 0.0) .OR. > (c2(num) .NE. 0.0) .OR. (c3(num) .NE. 0.0) .OR. > (c4(num) .NE. 0.0) .OR. (c5(num) .NE. 0.0))) THEN IF (i .GT. order) order = i IF (0 .EQ. first(m)) first(m) = num last(m) = num num = num+1 ENDIF READ (97,10,IOSTAT=eof) m,n, > c1(num), c2(num), c3(num), c4(num), c5(num), > e1(num), e2(num), e3(num), e4(num), e5(num), e6(num) IF (.not. path) c5(num) = 0. ENDDO CLOSE (97) RETURN END SUBROUTINE hmsApplyCOSY (u,num) IMPLICIT NONE REAL u(6) INTEGER num * -- apply a COSY matrix on the COSY vector u * (used in the forward tracking only) * * Parameter * u IO : coordinate vector (COSY) * u(1,2) : x [m], A = out of plane coordinates (downwards) * u(3,4) : y [m], B = inplane coordinates (perp. on x,z) * u(5) : l [m] = [?] = path length deviation (?) * u(6) : delta = relative deviation of the particle * momentum from p0 * num I : cosy matrix to apply * ! transport map common block INTEGER NMAPS,NLINES PARAMETER (NMAPS = 20) PARAMETER (NLINES = 10000) INTEGER first(NMAPS),last(NMAPS),order INTEGER e1(NLINES),e2(NLINES),e3(NLINES) INTEGER e4(NLINES),e5(NLINES),e6(NLINES) REAL c1(NLINES),c2(NLINES),c3(NLINES) REAL c4(NLINES),c5(NLINES) COMMON /hmsTrack/first,last,order, > e1,e2,e3,e4,e5,e6,c1,c2,c3,c4,c5 ! other variables REAL a,uu1(0:10),uu2(0:10),uu3(0:10) REAL uu4(0:10),uu5(0:10),uu6(0:10) INTEGER i ! calculate the powers of the focal plane coordinates uu1(0) = 1. uu2(0) = 1. uu3(0) = 1. uu4(0) = 1. uu5(0) = 1. uu6(0) = 1. uu1(1) = u(1) uu2(1) = u(2) uu3(1) = u(3) uu4(1) = u(4) uu5(1) = u(5) uu6(1) = u(6) DO i=2,order uu1(i)=uu1(i-1)*uu1(1) uu2(i)=uu2(i-1)*uu2(1) uu3(i)=uu3(i-1)*uu3(1) uu4(i)=uu4(i-1)*uu4(1) uu5(i)=uu5(i-1)*uu5(1) uu6(i)=uu6(i-1)*uu6(1) ENDDO DO i=1,5 u(i) = 0. ENDDO ! apply the cosy matrix DO i = first(num),last(num) a = uu1(e1(i)) * uu2(e2(i)) * uu3(e3(i)) * > uu4(e4(i)) * uu5(e5(i)) * uu6(e6(i)) u(1) = u(1) + c1(i) * a u(2) = u(2) + c2(i) * a u(3) = u(3) + c3(i) * a u(4) = u(4) + c4(i) * a u(5) = u(5) + c5(i) * a ENDDO RETURN END SUBROUTINE hmsLab2COSY (u) IMPLICIT NONE REAL u(6) * -- convert from LAB coordinates to COSY coordinates * * Parameter * u IO : coordinate vector (lab -> COSY) * u(1,2) : x [m], dx/dy -> A = out of plane coord. (downwards) * u(3,4) : y [m], dy/dz -> B = inplane coord. (perp. on x,z) * u(5) : l [m] = [?] = path length deviation (?) * u(6) : delta = relative deviation of the particle * momentum from p0 REAL delta,mep0 COMMON /hmsTrackLabCOSY/delta,mep0 REAL a delta = u(6) a = (1.+delta)/SQRT(1.+u(2)**2+u(4)**2) u(2) = a*u(2) u(4) = a*u(4) u(6) = (SQRT((1.+delta)**2+mep0**2)-SQRT(1.+mep0**2))/ > (SQRT(1.+mep0**2)-mep0) RETURN END SUBROUTINE hmsCOSY2Lab (u) IMPLICIT NONE REAL u(6) * -- convert from COSY coordinates to LAB coordinates * * Parameter * u IO : coordinates vector (COSY -> lab) * u(1,2) : x [m], A -> dx/dy = out of plane coord. (downwards) * u(3,4) : y [m], B -> dy/dz = inplane coord. (perp. on x,z) * u(5) : l [m] = [?] = path length deviation (?) * u(6) : delta = relative deviation of the particle * momentum from p0 REAL delta,mep0 COMMON /hmsTrackLabCOSY/delta,mep0 REAL a a = 1./SQRT((1.+delta)**2-u(2)**2-u(4)**2) u(2) = a*u(2) u(4) = a*u(4) u(6) = delta RETURN END SUBROUTINE hmsForward (uT,zT,u,zz,lost) IMPLICIT none REAL uT(6),zT,u(6),zz LOGICAL lost * -- make a single step transport calculation (without treating the acceptance) * * Parameter: * uT,zT I : target coordinates * uT(1,2) : x [m], dx/dz = out of plane coord. (downwards) * uT(3,4) : y [m], dy/dz = inplane coord. (perp. on x,z) * uT(5) : l [?] = path length deviation (?) * uT(6) : delta; relative deviation of the particle * momentum from p0 * zT : z [m]; in axis coordinate (towards HMS) * u,zz O : focal plane coordinates * u(1,2) : x [m], dx/dz; out of plane coord. (downwards) * u(3,4) : y [m], dy/dz; inplane coord. (perp. on x,z) * u(5) : l [?]; path length deviation (?) * u(6) : delta; relative deviation of the particle * momentum from p0 * zz : z [m]; in axis coordinate (towards HMS) * position where the particle stops inside HMS * (here always z-calorimeter) * lost O : set to .true. if the particle is lost in the HMS, * otherwise .false. (here always false) ! other variables INTEGER i lost = .FALSE. zz = 26.44743 ! copy the target coordinates to the u vector DO i=1,6 u(i)=uT(i) ENDDO ! drift backwards to z=0 u(1)=u(1)-u(2)*zT u(3)=u(3)-u(4)*zT !go up to the hut CALL hmsLab2COSY (u) CALL hmsApplyCOSY (u,20) CALL hmsCOSY2Lab (u) RETURN END LOGICAL FUNCTION hmsCheckDipole (u,zz,dz) IMPLICIT NONE REAL u(6),zz,dz * -- check for the dipole aperture * * Parameter * u I : target coordinates (COSY or LAB) * u(1,2) : x [m], A or dx/dz = out of plane coord. (downwards) * u(3,4) : y [m], A or dy/dz = inplane coord. (perp. on x,z) * u(5) : l [?] = path length deviation (?) * u(6) : delta; relative deviation of the particle * momentum from p0 * zz IO : z [m]; in axis coordinate (towards HMS) * dz I : distance from last check point [m] zz = zz+dz hmsCheckDipole = .TRUE. IF (ABS(u(1)).GT.0.27940.OR.ABS(u(3)).GT.0.18415) THEN IF (ABS(u(1)).GT.0.34290.OR.ABS(u(3)).GT.0.20320) RETURN IF (ABS(u(1)).GT.0.27940.AND.ABS(u(3)).GT.0.12065) THEN IF (((ABS(u(1))-0.27940)**2 + > (ABS(u(3))-0.12065)**2).GT.0.06350**2) RETURN ENDIF IF (ABS(u(1)).GT.0.13970 .OR. > (ABS(u(1))-10.1852*ABS(u(3))).GT.2.069633) RETURN ENDIF hmsCheckDipole = .FALSE. RETURN END LOGICAL FUNCTION hmsCheckQuad (u,zz,dz,r) IMPLICIT NONE REAL u(6),zz,dz,r * -- check for the dipoleaperture * * Parameter * u I : target coordinates (COSY or LAB) * u(1,2) : x [m], A or dx/dz = out of plane coord. (downwards) * u(3,4) : y [m], A or dy/dz = inplane coord. (perp. on x,z) * u(5) : l [?] = path length deviation (?) * u(6) : delta; relative deviation of the particle * momentum from p0 * zz IO : z [m]; in axis coordinate (towards HMS) * dz I : distance from last check point [m] * r I : aperture radius [m] zz = zz+dz hmsCheckQuad = ((u(1)**2+u(3)**2) .GT. r**2) RETURN END LOGICAL FUNCTION hmsDriftOcta (u,zz,dz,x,y,m,b) IMPLICIT NONE REAL u(6),zz,dz,x,y,m,b * -- drift electron and check for the octagon * * Parameter * u I : target coordinates (COSY or LAB) * u(1,2) : x [m], A or dx/dz = out of plane coord. (downwards) * u(3,4) : y [m], A or dy/dz = inplane coord. (perp. on x,z) * u(5) : l [?] = path length deviation (?) * u(6) : delta; relative deviation of the particle * momentum from p0 * zz IO : z [m]; in axis coordinate (towards HMS) * dz I : distance from last check point [m] * x,y,m,b I : octagon coordinates (x,y,b [m], m [1]) zz=zz+dz u(1)=u(1)+u(2)*dz u(3)=u(3)+u(4)*dz hmsDriftOcta = (ABS(u(1)) .GT. x) .OR. (ABS(u(3)) .GT. y) .OR. > ((ABS(u(1))+m*ABS(u(3))) .GT. b) RETURN END LOGICAL FUNCTION hmsDriftTPlate (u,zz,dz,r,y) IMPLICIT NONE REAL u(6),zz,dz,r,y * -- drift electron and check for the dipole transition plate * * Parameter * u IO : coordinate vector * zz IO : z position [m] * dz I : drift distance * r,y I : aperture zz=zz+dz u(1)=u(1)+u(2)*dz u(3)=u(3)+u(4)*dz hmsDriftTPlate = ((u(1)**2+u(3)**2) .GT. r**2) .OR. > (ABS(u(3)) .GT. y) RETURN END LOGICAL FUNCTION hmsDriftCirc (u,zz,dz,r) IMPLICIT NONE REAL u(6),zz,dz,r * -- drift electron and check for circular aperture * * Parameter * u I : target coordinates (COSY or LAB) * u(1,2) : x [m], A or dx/dz = out of plane coord. (downwards) * u(3,4) : y [m], A or dy/dz = inplane coord. (perp. on x,z) * u(5) : l [?] = path length deviation (?) * u(6) : delta; relative deviation of the particle * momentum from p0 * zz IO : z [m]; in axis coordinate (towards HMS) * dz I : distance from last check point [m] * r,y I : transition plate aperture radius [m] zz=zz+dz u(1)=u(1)+u(2)*dz u(3)=u(3)+u(4)*dz hmsDriftCirc = ((u(1)**2+u(3)**2) .GT. r**2) RETURN END LOGICAL FUNCTION hmsDriftRect (u,zz,dz,x0,x,y0,y) IMPLICIT NONE REAL u(6),zz,dz,x0,x,y0,y * -- drift electron and check for rectangular aperture * * Parameter * u I : target coordinates (COSY or LAB) * u(1,2) : x [m], A or dx/dz = out of plane coord. (downwards) * u(3,4) : y [m], A or dy/dz = inplane coord. (perp. on x,z) * u(5) : l [?] = path length deviation (?) * u(6) : delta; relative deviation of the particle * momentum from p0 * zz IO : z [m]; in axis coordinate (towards HMS) * dz I : distance from last check point [m] * x,x0, * y,y0 I : rectangular aperture [m] (half size and offset) zz=zz+dz u(1)=u(1)+u(2)*dz u(3)=u(3)+u(4)*dz hmsDriftRect = (ABS(u(1)-x0) .GT. x) .OR. > (ABS(u(3)-y0) .GT. y) RETURN END SUBROUTINE hmsAccept (uT,zT,u,zz,lost) IMPLICIT none REAL uT(6),zT,u(6),zz LOGICAL lost * -- make a transport calculation to find the acceptance * * Parameter: * uT,zT I : target coordinates * uT(1,2) : x [m], dx/dz; out of plane coord. (downwards) * uT(3,4) : y [m], dy/dz; inplane coord. (perp. on x,z) * uT(5) : l [?]; path length deviation (?) * uT(6) : delta; relative deviation of the particle * momentum from p0 * zT : z [m]; in axis coordinate (towards HMS) * u,zz O : focal plane coordinates * u(1,2) : x [m], dx/dz; out of plane coord. (downwards) * u(3,4) : y [m], dy/dz; inplane coord. (perp. on x,z) * u(5) : l [?]; path length deviation (?) * u(6) : delta; relative deviation of the particle * momentum from p0 * zz : z [m]; in axis coordinate (towards HMS) * position where the particle stops inside HMS * lost O : set to .true. if the particle is lost in the HMS, * otherwise .false. LOGICAL hmsCheckQuad LOGICAL hmsCheckDipole LOGICAL hmsDriftOcta LOGICAL hmsDriftTPlate LOGICAL hmsDriftCirc LOGICAL hmsDriftRect REAL uS(6) INTEGER i lost = .TRUE. zz = zT ! copy the target coordinates to the u vector DO i=1,6 u(i)=uT(i) ENDDO ! ----------------------------------------------------------- sive slit ! drift to sieve slit: z=1.27636m, ! octagon edge m=2.545977, b=0.13502640m IF (hmsDriftOcta(u,zz,1.27636-zT,0.0900176,0.0353568, > 2.5459770,0.1350264)) RETURN ! drift to back of sieve slit: z=1.33986m, dz=0.0635m ! octagon edge m=2.546569, b=0.141655189 m IF (hmsDriftOcta(u,zz,0.0635,0.0944368,0.0370839, > 2.546569,0.141655189)) RETURN ! ------------------------------------------------------------------ Q1 ! drift to mechanical entrance of Q1: z=1.4960m, dz=0.15614m IF (hmsDriftCirc(u,zz,0.15614,0.202575)) RETURN ! drift to Q1 entrance EFB: z=1.775805635m, dz=0.279805635m u(1)=u(1)+u(2)*0.279805635 u(3)=u(3)+u(4)*0.279805635 ! convert to COSY coordinates: CALL hmsLab2COSY (u) ! and save values DO i=1,6 uS(i)=u(i) ENDDO ! transport through Q1 fringe fields: CALL hmsApplyCOSY (u,1) IF (hmsCheckQuad(u,zz,0.279805635,0.202575)) RETURN ! transport through Q1, 1/5 at a time: DO i=1,4 CALL hmsApplyCOSY (u,2) IF (hmsCheckQuad(u,zz,1.87838873/5.,0.202575)) RETURN ENDDO ! restore values: DO i=1,6 u(i)=uS(i) ENDDO ! transport to Q1 exit EFB: CALL hmsApplyCOSY (u,3) IF (hmsCheckQuad(u,zz,1.87838873/5.,0.202575)) RETURN ! convert back LAB coordinates: CALL hmsCOSY2Lab (u) ! drift to Q1 mechanical exit: z=3.9340m, dz=0.279805635m IF (hmsDriftCirc(u,zz,0.279805635,0.202575)) RETURN ! ------------------------------------------------------------------ Q2 ! drift to mechanical entrance of Q2: z=4.5610m, dz=0.6270m IF (hmsDriftCirc(u,zz,0.6270,0.29840)) RETURN ! drift to Q2 entrance EFB: z=4.887021890m, dz=0.326021890m u(1)=u(1)+u(2)*0.326021890 u(3)=u(3)+u(4)*0.326021890 ! convert to COSY coordinates: CALL hmsLab2COSY (u) ! and save values DO i=1,6 uS(i)=u(i) ENDDO ! transport through Q2 fringe fields: CALL hmsApplyCOSY (u,4) IF (hmsCheckQuad(u,zz,0.326021890,0.29840)) RETURN ! transport through Q2, 1/5 at a time: DO i=1,4 CALL hmsApplyCOSY (u,5) IF (hmsCheckQuad(u,zz,2.15595622/5.,0.29840)) RETURN ENDDO ! restore values: DO i=1,6 u(i)=uS(i) ENDDO ! transport to Q2 exit EFB: CALL hmsApplyCOSY (u,6) IF (hmsCheckQuad(u,zz,2.15595622/5.,0.29840)) RETURN ! convert back LAB coordinates: CALL hmsCOSY2LAB (u) ! drift to Q2 mechanical exit: z=7.3690m, dz=0.326021890m IF (hmsDriftCirc(u,zz,0.326021890,0.29840)) RETURN ! ------------------------------------------------------------------ Q3 ! drift to mechanical entrance of Q3: z=7.6610m, dz=0.2920m IF (hmsDriftCirc(u,zz,0.2920,0.29840)) RETURN ! drift to Q3 entrance EFB: z=7.990200290m, dz=0.329200290m u(1)=u(1)+u(2)*0.329200290 u(3)=u(3)+u(4)*0.329200290 ! convert to COSY coordinates: CALL hmsLab2COSY (u) ! save values DO i=1,6 uS(i)=u(i) ENDDO ! transport through Q3 fringe fields: CALL hmsApplyCOSY (u,7) IF (hmsCheckQuad(u,zz,0.329200290,0.29840)) RETURN DO i=1,4 CALL hmsApplyCOSY (u,8) IF (hmsCheckQuad(u,zz,2.14959942/5.,0.29840)) RETURN ENDDO ! and restore values: DO i=1,6 u(i)=uS(i) ENDDO ! transport to Q3 exit EFB: CALL hmsApplyCOSY (u,9) IF (hmsCheckQuad(u,zz,2.14959942/5.,0.29840)) RETURN ! convert back LAB coordinates: CALL hmsCOSY2Lab (u) ! drift to Q3 mechanical exit: z=10.4690m, dz=0.329200290m IF (hmsDriftCirc(u,zz,0.329200290,0.29840)) RETURN ! -------------------------------------------------------------- Dipole ! drift to transition plate: z=11.058002m, dz=0.589002m IF (hmsDriftTPlate(u,zz,0.589002,0.30480,0.205232)) RETURN ! drift to opposite side of transition plate: z=11.092800m, dz=0.034798m IF (hmsDriftTPlate(u,zz,0.034798,0.30480,0.205232)) RETURN IF (hmsCheckDipole(u,zz,0.)) RETURN ! drift to D magnetic entrance: z=11.55m, dz=0.4572m u(1)=u(1)+u(2)*0.4572 u(3)=u(3)+u(4)*0.4572 IF (hmsCheckDipole(u,zz,0.4572)) RETURN ! convert to COSY coordinates: CALL hmsLab2COSY (u) ! save values: DO i=1,6 uS(i)=u(i) ENDDO ! transport through 1/5 D with rotated entrance face: CALL hmsApplyCOSY (u,10) IF (hmsCheckDipole(u,zz,5.26053145/5.)) RETURN ! transport through 3/5 D with sector segments: DO i=1,3 CALL hmsApplyCOSY (u,11) IF (hmsCheckDipole(u,zz,5.26053145/5.)) RETURN ENDDO ! restore values: DO i=1,6 u(i)=uS(i) ENDDO ! transport through D (entrance to exit, fringe fields included): CALL hmsApplyCOSY (u,13) IF (hmsCheckDipole(u,zz,5.26053145/5.)) RETURN ! convert back LAB coordinates: CALL hmsCOSY2Lab (u) ! drift to transition plate: z=0.4572m, dz=0.4572m IF (hmsDriftTPlate(u,zz,0.457200,0.34290,0.205232)) RETURN IF (hmsCheckDipole(u,zz,0.)) RETURN ! drift to opposite side of transition plate: z=0.491998m,dz=0.034798m IF (hmsDriftTPlate(u,zz,0.034798,0.34290,0.205232)) RETURN ! drift to end of first piece of telescope: z=1.119378m, dz=0.62738m IF (hmsDriftCirc(u,zz,0.62738,0.338450)) RETURN ! drift to end of second piece of telescope: z=4.086098m, dz=2.96672m IF (hmsDriftCirc(u,zz,2.96672,0.384175)) RETURN ! drift to end of third piece of telescope: z=5.578398m, dz=1.4923m IF (hmsDriftCirc(u,zz,1.49230,0.460375)) RETURN ! ----------------------------------------------------------------- hut ! drift to focal plane. This is the reference point for detector positions. u(1)=u(1)+u(2)*0.671602 u(3)=u(3)+u(4)*0.671602 zz = zz + 0.671602 DO i=1,6 uS(i)=u(i) ENDDO ! drift to DC1 entrance: z=-0.51923-0.036=-0.55523m, dz=-0.55523m IF(hmsDriftRect(u,zz,-0.55523,-0.01670,0.565,-0.00343,0.26))RETURN ! drift to DC1 exit: z=-0.51923+0.054=-0.46523m, dz=0.090m IF(hmsDriftRect(u,zz, 0.09000,-0.01670,0.565,-0.00343,0.26))RETURN ! drift to DC2 entrance: z=0.29299-0.036=0.25699m, dz=0.72222m IF(hmsDriftRect(u,zz, 0.72222,-0.02758,0.565,-0.01653,0.26))RETURN ! drift to DC2 exit: z=0.29299+0.054=0.34699m, dz=0.090m IF(hmsDriftRect(u,zz, 0.09000,-0.02758,0.565,-0.01653,0.26))RETURN ! drift to S1X: z=0.7783m, dz=0.43131m IF (hmsDriftRect(u,zz,0.43131,0.015,0.6025,0.000,0.3775)) RETURN ! drift to S1Y: z=0.9752m, dz=0.1969m IF (hmsDriftRect(u,zz,0.19690,0.000,0.6025,0.001,0.3775)) RETURN ! skip CK - no survey information ! drift to S2X: z=2.9882m, dz=2.013m IF (hmsDriftRect(u,zz,2.01300,0.004,0.6025,0.000,0.3775)) RETURN ! drift to S2Y: z=3.1851m, dz=0.1969m IF (hmsDriftRect(u,zz,0.19690,0.000,0.6025,0.013,0.3775)) RETURN ! drift to CAL: z=3.3869m, dz=0.2018m IF (hmsDriftRect(u,zz,0.20180,-0.134,0.6000,0.000,0.3000)) RETURN ! --------------------------------------------------------------- done lost = .FALSE. DO i=1,6 u(i)=uS(i) ENDDO RETURN END *------------------------------------------------------------------------ * * PART 2: HMS Reconstruction (Backward Tracking; Focal Plane to Target) * -=======- * * Reconstruction (backward tracking) in the Jlab HMS hall C * spectrometer using reconstruction and forward COSY matrices * (including the effects of beam offsets (out-of plane)) * * Both the normal in-plane scattering and the more special * out-of-plane scattering are handeled. The later makes use * of the forward COSY matrices. The algorithm was tested for * beam offsets in the range of cm (up or below the * nominal scattering plane) * * Supplies: * hmsInitRecon (map) * load the reconstruction maps * hmsInPlane (u,uT) * reconstruction of the target coordinates * (delta, dx/dz, y, dy/dz) at z=0 * hmsOutOfPlane (u,x,uT) * reconstruction of the target coordinates * (delta, dx/dz, y, dy/dz) at z=0 including the * vertical beam offset * * Note: - Before calling hmsReconInPlane or hmsReconOutOfPlaneAccept * the reconstruction map has to be loaded by a call to * hmsInitRecon * - Before calling hmsReconOutOfPlane the forward transport * maps have to be loaded by a call to hmsInitForward *------------------------------------------------------------------------ SUBROUTINE hmsInitRecon (filename) IMPLICIT NONE CHARACTER filename*(*) * -- load the HMS reconstruction map * * Parameter: * filename I : name of the reconstruction map * * File format: * Ax' Ay Ay' Ad XxYy * * Ax' : coeff's giving target x' * Ay : coeff's giving target y * Ay' : coeff's giving target y' * Ad : coeff's giving the delta * X : exponent of the focal plane x * x : exponent of the focal plane x' * Y : exponent of the focal plane y * y : exponent of the focal plane y' * ! matrix elements for the reconstruction INTEGER NTERMS PARAMETER (NTERMS=1000) INTEGER e1(NTERMS),e2(NTERMS),e3(NTERMS),e4(NTERMS),num,order REAL c1(NTERMS),c2(NTERMS),c3(NTERMS),c4(NTERMS) COMMON /hmsRecon/num,order,e1,e2,e3,e4,c1,c2,c3,c4 ! other variables INTEGER i,eof ! read the necessary coeffs from the map data file 10 FORMAT (4(E14.7E2,X),4I1) OPEN (UNIT=97,FILE=filename,STATUS='OLD') num = 1 order = 0 READ (97,10,IOSTAT=eof) > c2(num), c3(num), c4(num), c1(num), > e1(num), e2(num), e3(num), e4(num) DO WHILE (eof .GE. 0) IF ((c1(num) .NE. 0.0) .OR. (c2(num) .NE. 0.0) .OR. > (c3(num) .NE. 0.0) .OR. (c4(num) .NE. 0.0)) THEN i = e1(num) + e2(num) + e3(num) + e4(num) IF (i .GT. order) order = i num = num+1 ENDIF READ (97,10,IOSTAT=eof) > c2(num), c3(num), c4(num), c1(num), > e1(num), e2(num), e3(num), e4(num) ENDDO num = num-1 CLOSE (97) RETURN END SUBROUTINE hmsReconOffset (uT,du) IMPLICIT NONE REAL uT(6), du(4) * -- calculates the focal plane correction for a given * set of target coordinates * * Parameter * uT I : target coordinates (lab) * uT(1,2) : x [m], dx/dz = out of plane coord. (downwards) * uT(3,4) : y [m], dy/dz = inplane coord. (perp. on x,z) * uT(5) : z [m] = in axis coordinate (towards HMS) * uT(6) : delta = relative deviation of the particle * momentum from p0 * du O : correction values for the focal plane quantities * du(1,2) : correction in x [m], dx/dz * du(3,4) : correction in y [m], dy/dz ! transport map common block INTEGER NMAPS,NLINES PARAMETER (NMAPS = 20) PARAMETER (NLINES = 10000) INTEGER first(NMAPS),last(NMAPS),order INTEGER e1(NLINES),e2(NLINES),e3(NLINES) INTEGER e4(NLINES),e5(NLINES),e6(NLINES) REAL c1(NLINES),c2(NLINES),c3(NLINES) REAL c4(NLINES),c5(NLINES) COMMON /hmsTrack/first,last,order, > e1,e2,e3,e4,e5,e6,c1,c2,c3,c4,c5 REAL delta,mep0 COMMON /hmsTrackLabCOSY/delta,mep0 ! other variables REAL a,b,uu1(0:10),uu2(0:10),uu3(0:10),uu4(0:10),uu6(0:10),u(4) INTEGER i uu1(0) = 1 uu2(0) = 1 uu3(0) = 1 uu4(0) = 1 uu6(0) = 1 ! drift backwards to z=0 uu1(1)=uT(1)-uT(2)*uT(5) uu3(3)=uT(3)-uT(4)*uT(5) ! convert to COSY coordinates a = (1.+uT(6))/SQRT(1.+uT(2)**2+uT(4)**2) uu2(1) = a*uT(2) uu4(1) = a*uT(4) uu6(1) = (SQRT((1.+uT(6))**2+mep0**2)-SQRT(1.+mep0**2))/ > (SQRT(1.+mep0**2)-mep0) ! calculate the powers of the COSY coordinates DO i=1,4 du(i) = 0 u (i) = 0 ENDDO DO i=2,order uu1(i)=uu1(i)*uu1(1) uu2(i)=uu2(i)*uu2(1) uu3(i)=uu3(i)*uu3(1) uu4(i)=uu4(i)*uu4(1) uu6(i)=uu6(i)*uu6(1) ENDDO ! calculate the focal plane offsets DO i = first(20), last(20) a = uu2(e2(i)) * uu3(e3(i)) * uu4(e4(i)) * uu6(e6(i)) IF (e1(i) .EQ. 0) THEN u(2) = u (2)+c2(i) * a u(4) = u (4)+c4(i) * a ELSE a = uu1(e1(i)) * a du(1) = du(1)+c1(i) * a du(2) = du(2)+c2(i) * a du(3) = du(3)+c3(i) * a du(4) = du(4)+c4(i) * a ENDIF ENDDO ! convert COSY coordinates to differences in lab coordinates du(2) = du(2)+u(2) du(4) = du(4)+u(4) a = 1/SQRT((1.+uT(6))**2-du(2)**2-du(4)**2) b = 1/SQRT((1.+uT(6))**2-u (2)**2-u (4)**2) du(2) = du(2)*a-u(2)*b du(4) = du(4)*a-u(4)*b RETURN END SUBROUTINE hmsReconInPlane (u,uT) IMPLICIT NONE REAL u(4),uT(6) * -- performs the reconstruction of the target coordinates * (delta, dx/dz, y, dy/dz) at z=0 * * Parameter: * u I : focal plane coordinates (lab) * u(1,2) : x [m], dx/dz = out of plane coord. (downwards) * u(3,4) : y [m], dy/dz = inplane coord. (perp. on x,z) * uT O : target coordinates (lab) * uT(1,2) : x [m], dx/dz = out of plane coord. (downwards) * uT(3,4) : y [m], dy/dz = inplane coord. (perp. on x,z) * uT(5) : z [m] = in axis coordinate (towards HMS) * uT(6) : delta (relative deviation of the particle * momentum from p0) ! matrix elemnts needed for calculating the focal plane offset INTEGER NTERMS PARAMETER (NTERMS=1000) INTEGER e1(NTERMS),e2(NTERMS),e3(NTERMS),e4(NTERMS),num,order REAL c1(NTERMS),c2(NTERMS),c3(NTERMS),c4(NTERMS) COMMON /hmsRecon/num,order,e1,e2,e3,e4,c1,c2,c3,c4 ! other variables REAL a,uu1(0:10),uu2(0:10),uu3(0:10),uu4(0:10) INTEGER i ! calculate the powers of the focal plane coordinates DO i=1,6 uT(i) = 0 ENDDO uu1(0) = 1. uu2(0) = 1. uu3(0) = 1. uu4(0) = 1. uu1(1) = u(1) uu2(1) = u(2) uu3(1) = u(3) uu4(1) = u(4) DO i=2,order uu1(i)=uu1(i-1)*uu1(1) uu2(i)=uu2(i-1)*uu2(1) uu3(i)=uu3(i-1)*uu3(1) uu4(i)=uu4(i-1)*uu4(1) ENDDO ! calculate the target coordinates DO i = 1, num a = uu1(e1(i)) * uu2(e2(i)) * uu3(e3(i)) * uu4(e4(i)) uT(6) = uT(6)+c1(i) * a uT(2) = uT(2)+c2(i) * a uT(3) = uT(3)+c3(i) * a uT(4) = uT(4)+c4(i) * a ENDDO RETURN END SUBROUTINE hmsReconOutOfPlane (u,x,uT) IMPLICIT NONE REAL u(4),x,uT(6) * -- performs the reconstruction of the target coordinates * (delta, dx/dz, y, dy/dz) for the out of plane case * * Parameter: * u I : focal plane coordinates (lab) * u(1,2) : x [m], dx/dz = out of plane coord. (downwards) * u(3,4) : y [m], dy/dz = inplane coord. (perp. on x,z) * x I : x-offset at the target [m] (z=0) (out of plane; downwards) * uT O : target coordinates (lab) * uT(1,2) : x [m], dx/dz = out of plane coord. (downwards) * uT(3,4) : y [m], dy/dz = inplane coord. (perp. on x,z) * uT(5) : z [m] = in axis coordinate (towards HMS) * uT(6) : delta (relative deviation of the particle * momentum from p0) REAL eps PARAMETER (eps = 0.0005) REAL dd,du(4),u0(4) INTEGER n CALL hmsReconInPlane (u,uT) dd = 1 n = 0 DO WHILE ((abs(dd) .GT. eps) .AND. (n .LT. 10)) uT(1) = x CALL hmsReconOffset (uT,du) u0(1) = u(1)-du(1) u0(2) = u(2)-du(2) u0(3) = u(3)-du(3) u0(4) = u(4)-du(4) dd = uT(6) CALL hmsReconInPlane (u0,uT) dd = dd-uT(6) n = n+1 ENDDO uT(1) = x RETURN END