*------------------------------------------------------------------------ * * TRG_TRACK GEN Target Rracking routines * -=========- * * Raytracing of 3-d motion in polarized target field by solution * of differential equations of motion via 4th order Runge-Kutta. Field * orientation arbitrary. * * Note: - the HMS routines use a right handed coord. system with * x : pointing downwards * y : perpendicular to x,z, * pointing to the left (if seen in z-direction) * z : BEAM axis or HMS axis, pointing downstream or * from the target to the focal plane respectively * * - the B field map uses a cylindrical coordinate system * with z along the field axis and r perpendicular to it * * - all length (x,y,z,dl,l,...) are measured in [cm] * - all velocities are measured in [cm/ns] * - all angles are measured counter clock wise in [deg] * - time is measured in [ns] * - the B field is measured in [T] * * original devloped by ??? * widely modified by MM * - converted into subroutines * - rotation algorytm correced (at moment: phi==0 assumed) * - changed coordinate system * beam direction: z * horizontal plane: zy * out of plane: x (points downwards) * * Supplies: * trgInit (map,theta,phi) * load the target field map * trgTrack (u,E,dl,l) * tracks a single particle with given start parameter * trgXTrack (u,E,dl,l,id,Bdl) * tracks a single particle with given start parameter * storing the track in a PAW/hbook ntuple * trgTrackToPlane (u,E,dl,a,b,c,d) * track a single particle with given start parameters * and find the intersection of the particle track with * a given plane * * Note: - Before calling trgTrack,trgXTrack or trgTrackToPlane * the target field map has to be loaded by a call to * trgInit *------------------------------------------------------------------------ SUBROUTINE trgTrack (u,E,dl,l) IMPLICIT NONE REAL u(6),E,dl,l * -- track a single particle with given start parameters * * Parameter: * u IO : coordinate vector (initial/final) * u(1,2,3) : x, y, z [cm] * u(4,5,6) : dx/dt, dy/dt, dz/dt [cm/ns] * E I : particle energy [MeV] * sign of particle charge * (negative for electrons, positive for protons/deuterons) * dl I : step size [cm] * l I : tracking distance [cm] REAL factor COMMON /trgConversionFactor/factor REAL ts INTEGER i,n factor = 90./E ts = dl/30. n = ABS(l/dl) DO i=1,n !step thru time CALL trgRK4(u,u,ts) !solve diff. eqs. ENDDO RETURN END SUBROUTINE trgXTrack (u,E,dl,l,Bdl,Xfun,id) IMPLICIT NONE REAL u(6),E,dl,l,Bdl INTEGER id INTEGER Xfun EXTERNAL Xfun * -- track a single particle with given start parameters and * writes the track's coordinates to a hbook file * * Parameter: * u IO : coordinate vector (initial/final) * u(1,2,3) : x, y, z [cm] * u(4,5,6) : dx/dt, dy/dt, dz/dt [cm/ns] * E I : particle energy [MeV] * sign of particle charge * (negative for electrons, positive for protons/deuterons) * dl I : step size [cm] * l I : tracking distance [cm] * Bdl O : Integrated Bdl [Tcm] * Xfun I : external function called after every iteration * int Xfun (id,t,l,u) * int id as passed to trgXtrack * real t,l,u(9) time, pathlength and act. coordinates * id I : histogram id REAL factor COMMON /trgConversionFactor/factor REAL ts,uu(9) INTEGER i,n,x factor = 90./E ts = dl/30. n = l/dl ! book start location DO i=1,6 uu (i) = u(i) ENDDO DO i=7,9 uu (i) = 0. ENDDO x = Xfun(id,0.,0.,uu) ! track the particle and book location DO i=1,n CALL trgRK4Bdl(uu,uu,ts) ! solve diff. eqs. x = Xfun (id,i*ts,i*dl,uu) ENDDO DO i=1,6 u(i) = uu (i) ENDDO ! calculate Bdl ( B_x^2+B_y^2+B_z^2 ) Bdl = SQRT(uu(7)**2+uu(8)**2+uu(9)**2) RETURN END SUBROUTINE trgTrackToPlane (u0,E,dl,a,b,c,d) IMPLICIT NONE REAL u0(6),E,dl,a,b,c,d * -- track a single particle with given start parameters * and find the intersection of the particle track with a given plane * * Parameter: * u0 IO : coordinate vector (initial/final) * u0(1,2,3) : x, y, z [cm] * u0(4,5,6) : dx/dt, dy/dt, dz/dt [cm/ns] * E I : particle energy [MeV] * sign of particle charge * (negative for electrons, positive for protons/deuterons) * dl I : step size [cm] * a..d I : parameter of the intersection plane * 0 = a*x+b*y+c*z+d; REAL factor COMMON /trgConversionFactor/factor REAL ts,n,an,bn,cn,dn,dist0,dist1,u1(6) INTEGER i LOGICAL back n = 1/SQRT (a*a+b*b+c*c) an = a*n bn = b*n cn = c*n dn = d*n factor = 90./E ts = dl/30. dist0 = u0(1)*an + u0(2)*bn + u0(3)*cn + dn ! check for the tracking direction CALL trgRK4(u0,u1,ts) dist1 = u1(1)*an + u1(2)*bn + u1(3)*cn + dn IF ((SIGN(1.,dist0) .EQ. SIGN(1.,dist1)) .AND. > (ABS(dist0) .LT. ABS(dist1))) ts=-ts ! track through the intersection plane dist1 = dist0 DO WHILE (SIGN(1.,dist0) .EQ. SIGN(1.,dist1)) CALL trgRK4(u0,u1,ts) dist1 = u1(1)*an + u1(2)*bn + u1(3)*cn + dn IF (SIGN(1.,dist0) .EQ. SIGN(1.,dist1)) THEN CALL trgRK4(u1,u0,ts) dist0 = u0(1)*an + u0(2)*bn + u0(3)*cn + dn ENDIF ENDDO ! calculate the intersection point DO i=1,6 u0(i) = u0(i) + (u1(i)-u0(i)) * dist0/(dist0-dist1) ENDDO RETURN END *------------------------------------------------------------------------------ * load the field map and calculate the magnetic field strength * SUBROUTINE trgInit (map,theta,phi) IMPLICIT NONE CHARACTER map*(*) REAL theta,phi * -- read field map (for calculations in the LAB system) * * Parameter: * map I : filename of the fieldmap (=' ': uniform field test case) * theta,phi I : inplane (theta) and out of plane (phi) angle * * note: currently phi is always treated as 0 * INTEGER nz,nr PARAMETER (nz = 51) PARAMETER (nr = 51) REAL B_field_z(nz,nr),B_field_r(nz,nr),zz(nz),rr(nr) REAL B_theta,B_stheta,B_ctheta,B_phi,B_sphi,B_cphi COMMON /trgFieldStrength/ B_field_z,B_field_r,zz,rr COMMON /trgFieldAngles/ B_theta,B_stheta,B_ctheta, > B_phi, B_sphi, B_cphi REAL pi180 PARAMETER (pi180 = 3.141592653/180.) INTEGER ir,iz REAL xx B_theta = theta B_stheta = SIN(theta*pi180) B_ctheta = COS(theta*pi180) ! Note: for performance reasons B_phi is always treated 0 in trgField B_phi = phi B_sphi = SIN(phi*pi180) B_cphi = COS(phi*pi180) IF (map .NE. ' ') THEN !read in numerical field map OPEN (unit=1,file=map,status='old') DO ir=1,nr rr(ir) = 2.*float(ir-1) zz(ir) = 2.*float(ir-1) DO iz=1,nz READ (1,*)xx,xx,B_field_z(iz,ir),B_field_r(iz,ir),xx,xx,xx ENDDO ENDDO CLOSE (unit=1) ELSE DO ir=1,nr ! uniform 5T field over 26 cm in z rr(ir) = 2.*float(ir-1) ! and 16 cm in r zz(ir) = 2.*float(ir-1) DO iz=1,nz B_field_r(iz,ir) = 0. IF (rr(ir) .LE. 16. .and. zz(ir) .LE. 26.) THEN B_field_z(iz,ir) = 5.0 ELSE B_field_z(iz,ir) = 0.0 ENDIF ENDDO ENDDO ENDIF RETURN END SUBROUTINE trgField (x_,B_) IMPLICIT NONE REAL x_(3),B_(3) * -- calculate actual field * * Parameter: * x_ I : lab coordinates * B_ O : B field in lab coordinates * * Notes: * - 2-Dimensional Linear Interpolation: * Assumes uniform spacing of fieldmap in x,y * - for performance reasons B_phi is always treated 0 INTEGER nz,nr PARAMETER (nz = 51) PARAMETER (nr = 51) REAL B_field_z(nz,nr),B_field_r(nz,nr),zz(nz),rr(nr) REAL B_theta,B_stheta,B_ctheta,B_phi,B_sphi,B_cphi COMMON /trgFieldStrength/ B_field_z,B_field_r,zz,rr COMMON /trgFieldAngles/ B_theta,B_stheta,B_ctheta, > B_phi, B_sphi, B_cphi INTEGER i,j REAL x(3),B(3),z,r,az,ar,a0,a1 ! rotate to coordinates with z' along field direction x(1) = x_(1) x(2) = B_stheta*x_(3) + B_ctheta*x_(2) x(3) = B_ctheta*x_(3) - B_stheta*x_(2) ! compute zylinder coordinates z = ABS (x(3)) r = SQRT (x(1)**2 + x(2)**2) ! interpolate the field map i = INT((z-zz(1))/(zz(2)-zz(1))) + 1 j = INT((r-rr(1))/(rr(2)-rr(1))) + 1 IF ((i+1 .GT. nz) .OR. (i .LT. 1) .OR. > (j+1 .GT. nr) .OR. (j .LT. 1)) THEN B(1)=0. B(2)=0. B(3)=0. ELSE ! calculate the Bz component az = ((z-zz(i))/(zz(2)-zz(1))) ar = ((r-rr(j))/(rr(2)-rr(1))) a0=az*(B_field_z(i+1,j) -B_field_z(i,j)) +B_field_z(i,j) a1=az*(B_field_z(i+1,j+1)-B_field_z(i,j+1))+B_field_z(i,j+1) B(3) = (ar*(a1-a0)+a0) IF (r .gt. 0.) THEN ! calculate the Bx,By components a0=az*(B_field_r(i+1,j) -B_field_r(i,j)) +B_field_r(i,j) a1=az*(B_field_r(i+1,j+1)-B_field_r(i,j+1))+B_field_r(i,j+1) B(2) = (ar*(a1-a0)+a0)/r IF (x(3) .LT. 0.) B(2)= -B(2) B(1) = B(2)*x(1) B(2) = B(2)*x(2) ! transform B field to lab. system B_(1) = B(1) B_(2) = - B_stheta*B(3) + B_ctheta*B(2) B_(3) = B_ctheta*B(3) + B_stheta*B(2) ELSE B_(1) = 0. B_(2) = - B_stheta*B(3) B_(3) = B_ctheta*B(3) ENDIF ENDIF RETURN END *------------------------------------------------------------------------------ * solve the differential equation of the particle * SUBROUTINE trgDeriv(u,dudt) IMPLICIT NONE REAL u(9),dudt(9) * -- calculate the derivatives du(i)/dt for the runke kutta routine * * Parameter: * u I : actual coordinate vector * u(1,2,3) I : x, y, z * u(4,5,6) I : dx/dt, dy/dt, dz/dt * u(7,8,9) I : integral Bxdx, Bydy, Bzdz * dudt O : derivative du/dt * dudt(1,2,3) : dx/dt, dy/dt, dz/dt * dudt(4,5,6) : d^2xdt^2, d^2ydt^2, d^2zdt^2 * dudt(7,8,9) : B x v REAL factor COMMON /trgConversionFactor/factor REAL B(3) CALL trgField (u,B) ! These are just the velocities dudt(1) = u(4) dudt(2) = u(5) dudt(3) = u(6) ! This is just (v_vec X B_vec) dudt(7) = u(5)*B(3) - u(6)*B(2) dudt(8) = u(6)*B(1) - u(4)*B(3) dudt(9) = u(4)*B(2) - u(5)*B(1) ! This is just (v_vec X B_vec) * factor dudt(4) = dudt(7)*factor dudt(5) = dudt(8)*factor dudt(6) = dudt(9)*factor RETURN END SUBROUTINE trgRK4(u0,u1,h) IMPLICIT NONE REAL u0(6),u1(6),h * -- Fourth-order Runge-Kutta from Numerical Recipes book * for tracking through the target field * * Parameter: * u0 I : input coordinate vector * u1 O : output coordinate vector * u(1,2,3) : x, y, z * u(4,5,6) : dx/dt, dy/dt, dz/dt * h I : time step INTEGER i REAL ut(6),dudt(9),dut(9),dum(9),hh,h6 hh=h*0.5 h6=h/6. CALL trgDeriv(u0,dudt) DO i=1,6 ut(i) = u0(i) + hh*dudt(i) ENDDO CALL trgDeriv(ut,dut) DO i=1,6 ut(i) = u0(i) + hh*dut(i) ENDDO CALL trgDeriv(ut,dum) DO i=1,6 ut(i) = u0(i) +h*dum(i) dum(i)= dut(i) +dum(i) ENDDO CALL trgDeriv(ut,dut) DO i=1,6 u1(i)=u0(i)+h6*(dudt(i)+dut(i)+2.*dum(i)) ENDDO RETURN END SUBROUTINE trgRK4Bdl(u0,u1,h) IMPLICIT NONE REAL u0(9),u1(9),h * -- Fourth-order Runge-Kutta from Numerical Recipes book * for tracking through the target field (incl. B/dl calculation) * * Parameter: * u0 I : input coordinate vector * u1 O : output coordinate vector * u(1,2,3) : x, y, z * u(4,5,6) : dx/dt, dy/dt, dz/dt * u(7,8,9) : integral Bxdx, Bydy, Bzdz * h I : time step INTEGER i REAL ut(9),dudt(9),dut(9),dum(9),hh,h6 hh=h*0.5 h6=h/6. CALL trgDeriv(u0,dudt) DO i=1,9 ut(i) = u0(i) + hh*dudt(i) ENDDO CALL trgDeriv(ut,dut) DO i=1,9 ut(i) = u0(i) + hh*dut(i) ENDDO CALL trgDeriv(ut,dum) DO i=1,9 ut(i) = u0(i) +h*dum(i) dum(i)= dut(i) +dum(i) ENDDO CALL trgDeriv(ut,dut) DO i=1,9 u1(i)=u0(i)+h6*(dudt(i)+dut(i)+2.*dum(i)) ENDDO RETURN END