*------------------------------------------------------------------------ * * 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,ok) * 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 * * 99/11/24 GAW: in trgInit - fixed bug in uniform field * - added zero field option * in trgTrackToPlane - fixed bug in which routine took * an extra step if crossed the targeted plane * in the first step * 99/12/06 GAW: numerous modifications to work with 2 spectrometers *------------------------------------------------------------------------ SUBROUTINE trgTrack (u,E,dl,l,spect) IMPLICIT NONE REAL u(6),E,dl,l INTEGER spect * -- 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/sqrt(u(4)**2+u(5)**2+u(6)**2) n = ABS(l/dl) DO i=1,n !step thru time CALL trgRK4(u,u,ts,spect) !solve diff. eqs. ENDDO RETURN END SUBROUTINE trgXTrack (u,E,dl,l,Bdl,Xfun,id,spect) IMPLICIT NONE REAL u(6),E,dl,l,Bdl INTEGER id,spect 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/sqrt(u(4)**2+u(5)**2+u(6)**2) 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,spect) ! 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 (u,E,dl,a,b,c,d,ok,spect) IMPLICIT NONE REAL u(6),E,dl,a,b,c,d INTEGER spect LOGICAL ok * -- track a single particle with given start parameters * and find the intersection of the particle track with a given plane * * Parameter: * u 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; * ok IO : status variable * - if false no action is taken * - set to false when no intersection point is found * REAL factor COMMON /trgConversionFactor/factor REAL ts,n,an,bn,cn,dn,maxdist,dist0,dist1,u0(6),u1(6) INTEGER i,steps,max_steps IF (.NOT. OK) RETURN 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/sqrt(u(4)**2+u(5)**2+u(6)**2) dist0 = u(1)*an + u(2)*bn + u(3)*cn + dn maxdist = max(ABS(dist0)*4.,1.0) ! check for the tracking direction CALL trgRK4(u,u1,ts,spect) 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 ! GAW 99/11/22 ! Previously, if dist1 and dist0 had different signs, it move the track one ! extra step so that interpolation in end is wrong. The added if prevents that. steps = 0 max_steps = int(max(dist0,10.*dl)/dl)*10 IF (SIGN(1.,dist0).eq.SIGN(1.,dist1)) THEN dist1 = dist0 DO WHILE ((SIGN(1.,dist0) .EQ. SIGN(1.,dist1)) .AND. ok) CALL trgRK4(u1,u0,ts,spect) dist0 = u0(1)*an + u0(2)*bn + u0(3)*cn + dn IF (SIGN(1.,dist0) .EQ. SIGN(1.,dist1)) THEN CALL trgRK4(u0,u1,ts,spect) dist1 = u1(1)*an + u1(2)*bn + u1(3)*cn + dn ENDIF ok = (ABS(dist1) .LT. maxdist).and.steps.lt.max_steps C write(*,*) dist0,dist1 steps = steps+1 ENDDO ELSE DO i=1,6 u0(i) = u(i) ENDDO ENDIF IF (ok) THEN ! calculate the intersection point DO i=1,6 u(i) = u0(i) + (u1(i)-u0(i)) * dist0/(dist0-dist1) ENDDO ENDIF RETURN END *------------------------------------------------------------------------------ * load the field map and calculate the magnetic field strength * SUBROUTINE trgInit (map,theta_e,phi_e,theta_p,phi_p) IMPLICIT NONE CHARACTER map*(*) REAL theta_e,phi_e,theta_p,phi_p * -- read field map (for calculations in the LAB system) * * Parameter: * map I : filename of the fieldmap (=' ': uniform field test case) * theta_e,phi_e I : inplane(theta) & out of plane(phi) angle for e spect * theta_p,phi_p I : inplane(theta) & out of plane(phi) angle for p spect * * note: currently phi is always treated as 0 * * GAW 99/12/06 modified to work with 2 spectrometers 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_e,B_stheta_e,B_ctheta_e,B_phi_e,B_sphi_e,B_cphi_e REAL B_theta_p,B_stheta_p,B_ctheta_p,B_phi_p,B_sphi_p,B_cphi_p COMMON /trgFieldStrength/ B_field_z,B_field_r,zz,rr COMMON /trgFieldAngles_e/ B_theta_e,B_stheta_e,B_ctheta_e, > B_phi_e, B_sphi_e, B_cphi_e COMMON /trgFieldAngles_p/ B_theta_p,B_stheta_p,B_ctheta_p, > B_phi_p, B_sphi_p, B_cphi_p REAL pi180 PARAMETER (pi180 = 3.141592653/180.) INTEGER ir,iz REAL xx B_theta_e = theta_e B_stheta_e = SIN(theta_e*pi180) B_ctheta_e = COS(theta_e*pi180) B_theta_p = theta_p B_stheta_p = SIN(theta_p*pi180) B_ctheta_p = COS(theta_p*pi180) ! Note: for performance reasons B_phi is always treated 0 in trgField B_phi_e = phi_e B_sphi_e = SIN(phi_e*pi180) B_cphi_e = COS(phi_e*pi180) B_phi_p = phi_p B_sphi_p = SIN(phi_p*pi180) B_cphi_p = COS(phi_p*pi180) CGAW write(*,*) 'trginit',theta_e,theta_p ! GAW 99/11/22: Add zero field option IF (map.EQ.'0') THEN DO ir=1,nr rr(ir) = 2.*float(ir-1) zz(ir) = 2.*float(ir-1) DO iz=1,nz B_field_r(iz,ir) = 0. B_field_z(iz,ir) = 0. ENDDO ENDDO ! GAW 99/11/22: End ELSEIF (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 ! GAW 99/11/19: Must initialize rr and zz before going through loop since ! do tests on them. 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) ENDDO DO ir=1,nr ! uniform 5T field over 26 cm in z DO iz=1,nz B_field_r(iz,ir) = 0. IF (rr(ir) .LE. 16. .and. zz(iz) .LE. 26.) THEN CGAW B_field_z(iz,ir) = 0.0 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_,spect) IMPLICIT NONE REAL x_(3),B_(3) * -- calculate actual field * * Parameter: * x_ I : lab coordinates * B_ O : B field in lab coordinates * spect I: id for spectrometer (-1=e-, +1=p) * * Notes: * - 2-Dimensional Linear Interpolation: * Assumes uniform spacing of fieldmap in x,y * - for performance reasons B_phi is always treated 0 * * GAW 99/12/06 modified to work with 2 spectrometers 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_e,B_stheta_e,B_ctheta_e,B_phi_e,B_sphi_e,B_cphi_e REAL B_theta_p,B_stheta_p,B_ctheta_p,B_phi_p,B_sphi_p,B_cphi_p COMMON /trgFieldStrength/ B_field_z,B_field_r,zz,rr COMMON /trgFieldAngles_e/ B_theta_e,B_stheta_e,B_ctheta_e, > B_phi_e, B_sphi_e, B_cphi_e COMMON /trgFieldAngles_p/ B_theta_p,B_stheta_p,B_ctheta_p, > B_phi_p, B_sphi_p, B_cphi_p REAL B_tht, B_stht, B_ctht, B_ph, B_sph, B_cph INTEGER i,j,spect REAL x(3),B(3),z,r,az,ar,a0,a1 ! set up angles for electron or proton spectrometer if (spect.eq.-1) then B_tht = B_theta_e B_stht = B_stheta_e B_ctht = B_ctheta_e B_ph = B_phi_e B_sph = B_sphi_e B_cph = B_cphi_e else if (spect.eq.1) then B_tht = B_theta_p B_stht = B_stheta_p B_ctht = B_ctheta_p B_ph = B_phi_p B_sph = B_sphi_p B_cph = B_cphi_p endif ! rotate to coordinates with z' along field direction x(1) = x_(1) x(2) = B_stht*x_(3) + B_ctht*x_(2) x(3) = B_ctht*x_(3) - B_stht*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_stht*B(3) + B_ctht*B(2) B_(3) = B_ctht*B(3) + B_stht*B(2) ELSE B_(1) = 0. B_(2) = - B_stht*B(3) B_(3) = B_ctht*B(3) ENDIF ENDIF RETURN END *------------------------------------------------------------------------------ * solve the differential equation of the particle * SUBROUTINE trgDeriv(u,dudt,spect) 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 * spect I : -1 for e spectrometer, +1 for p spectrometer REAL factor COMMON /trgConversionFactor/factor INTEGER spect REAL B(3) CALL trgField (u,B,spect) ! 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,spect) 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 * spect I: -1 for e spectrometer, +1 for p spectrometer INTEGER i,spect REAL ut(6),dudt(9),dut(9),dum(9),hh,h6 hh=h*0.5 h6=h/6. CALL trgDeriv(u0,dudt,spect) DO i=1,6 ut(i) = u0(i) + hh*dudt(i) ENDDO CALL trgDeriv(ut,dut,spect) DO i=1,6 ut(i) = u0(i) + hh*dut(i) ENDDO CALL trgDeriv(ut,dum,spect) DO i=1,6 ut(i) = u0(i) +h*dum(i) dum(i)= dut(i) +dum(i) ENDDO CALL trgDeriv(ut,dut,spect) DO i=1,6 u1(i)=u0(i)+h6*(dudt(i)+dut(i)+2.*dum(i)) ENDDO RETURN END SUBROUTINE trgRK4Bdl(u0,u1,h,spect) 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 * spect I: -1 for e spectrometer, +1 for p spectrometer INTEGER i,spect REAL ut(9),dudt(9),dut(9),dum(9),hh,h6 hh=h*0.5 h6=h/6. CALL trgDeriv(u0,dudt,spect) DO i=1,9 ut(i) = u0(i) + hh*dudt(i) ENDDO CALL trgDeriv(ut,dut,spect) DO i=1,9 ut(i) = u0(i) + hh*dut(i) ENDDO CALL trgDeriv(ut,dum,spect) DO i=1,9 ut(i) = u0(i) +h*dum(i) dum(i)= dut(i) +dum(i) ENDDO CALL trgDeriv(ut,dut,spect) DO i=1,9 u1(i)=u0(i)+h6*(dudt(i)+dut(i)+2.*dum(i)) ENDDO RETURN END C------------------------------------------------------------------------------ ! ! Tracking of charged particles through the target field. ! ! Author: Glen Warren, December 1999 ! ! Much of the code is taken from Markus Muehlbauer's work on the Hall C replay ! engine. ! C------------------------------------------------------------------------------ subroutine track_from_tgt(x,y,z,dx,dy,mom,mass,spect,ok) C Given vertex coordinates, momentum and mass, tracks particle through a C field to field-free region 100 cm from target. It is assumed that the C code that calls this routine will reconstruct the track to z=0. C C x,y,z,dx,dy coordinates follow COSY a la Hall C: C z = into spectrometer C x = down (direction of increasing momentum) C y = z cross x. C dx = dx/dz C dy = dy/dz C C coordinates for tracking are: C vT(1,2,3) are the position in X,Y,Z [cm] C vT(4,5,6) are the velocity in the X,Y,Z direction [cm/ns]. implicit none real*8 x,y,z,dx,dy ! in: initial coords. out: coords of image track real*8 mom ! momentum (MeV). (mom<0 for e-, mom>0 for p,d) real*8 mass ! mass of particle (MeV) integer spect logical ok real*8 cc ! speed of light in cm/ns parameter (cc = 29.9792458) real*8 vel ! velocity of particle [cm/ns] real*8 eng ! energy of particle real vT(6) C write(*,*) 'from target',spect C write(*,*) mom,mass call print_coord2('init track, beam: ',x,y,z,dx,dy) vel = abs(mom)/sqrt(mom**2+mass**2)*cc eng = sign(1.,mom)*sqrt(mom**2+mass**2) vT(1) = x vT(2) = y vT(3) = z vT(6) = vel/sqrt(1+dx**2+dy**2) vT(4) = dx*vT(6) vT(5) = dy*vT(6) ! for debugging, run track first to z=0. ok = .true. call trgTrackToPlane(vT,eng,1.,0.,0.,1.,0.,ok,spect) C write(*,*) ok,spect call print_coord3('init track, z=0:',vt) ! track through magnetic field to z=100 cm plane ok = .true. call trgTrackToPlane(vT,eng,1.,0.,0.,1.,-100.,ok,spect) C write(*,*) ok,spect call print_coord3('init track, z=100:',vt) ! translate back into SIMC variables x = vT(1) y = vT(2) z = vT(3) dx = vT(4)/vT(6) dy = vT(5)/vt(6) return end C------------------------------------------------------------------------------ C Tracks proton through away from target through a target field without C considering additional optics-- appropriate for E93-026 - GAW C C This routine is a mostly a copy of what is done for the full Monte C Carlo of the proton arm. It would be aesthetically pleasing to rewrite C that part of simc.f so that the repeated code is not repeated. But, C that is aesthetics. subroutine transport_proton(main,P_arm_z) implicit none include 'simulate.inc' real*8 x_P_arm,y_P_arm,z_P_arm,dx_P_arm,dy_P_arm,delta_P_arm real*8 P_arm_z logical ok_P_arm record /event_main/ main delta_P_arm = main.SP.p.delta x_P_arm = -(main.target.y+spec.p.offset.y) if(hms_e_flag)then y_P_arm = -(main.target.z+spec.p.offset.z)*spec.p.sin_th+ > (main.target.x+spec.p.offset.x)*spec.p.cos_th z_P_arm = (main.target.x+spec.p.offset.x)*spec.p.sin_th+ > (main.target.z+spec.p.offset.z)*spec.p.cos_th else y_P_arm = (main.target.z+spec.p.offset.z)*spec.p.sin_th+ > (main.target.x+spec.p.offset.x)*spec.p.cos_th z_P_arm = -(main.target.x+spec.p.offset.x)*spec.p.sin_th+ > (main.target.z+spec.p.offset.z)*spec.p.cos_th endif dx_P_arm = main.SP.p.xptar - spec.p.offset.xptar dy_P_arm = main.SP.p.yptar - spec.p.offset.xptar ! project vertex coordinates to z=0 to be used to compare with reconstructed ! target positions P_arm_z = y_P_arm - z_P_arm*dy_P_arm call track_from_tgt(x_P_arm,y_P_arm,z_P_arm,dx_P_arm,dy_P_arm, > spec.p.P*(1+delta_P_arm/100.),Mp,+1,ok_P_arm) if (.not.ok_P_arm) dx_P_arm = -1000. ! ........ drift this to z=0 x_P_arm = x_P_arm - z_P_arm*dx_P_arm y_P_arm = y_P_arm - z_P_arm*dy_P_arm z_P_arm = 0.0 main.RECON.p.delta = delta_P_arm main.RECON.p.yptar = dy_P_arm main.RECON.p.xptar = dx_P_arm main.RECON.p.z = y_P_arm ! construct vectors at first plane of bars for the focal plane coordinates main.FP.p.x = x_P_arm + drift_to_ndet*dx_P_arm main.FP.p.y = y_P_arm + drift_to_ndet*dy_P_arm main.FP.p.dx = dx_P_arm main.FP.p.dy = dy_P_arm return end C------------------------------------------------------------------------------ subroutine track_to_tgt(delta,y,dx,dy,frx,fry,mom,mass,ctheta, > stheta,spect,ok) implicit none include 'simulate.inc' real*8 delta,y,dx,dy ! in: first guess reconstructed coords. ! out: final reconstructed coords. real*8 frx ! raster horizontal position (points left) real*8 fry ! raster vertical position (points up) real*8 mom ! momentum (MeV). (mom<0 for e-, mom>0 for p,d) real*8 mass ! mass of particle (MeV) real ctheta,stheta ! cosine and sine of central spectrometer angle integer spect logical ok real vT(6),vTx(6) real*4 xx,delx real*8 xxd integer*2 i,n real*8 vel,cc,eng,mom_0 parameter (cc=29.9792458) ! do reconstruction considering target field. Taken from gen_track.f from ! Markus Muehlbauer ! C write(*,*) 'to target',spect call print_coord1('after first mc_hms_recon',y,delta,dx,dy,fry,0.) ! ! copy vertical offset into another variable xx = -fry ! use first call to mc_hms_recon as a first guess. Next trace back 100 cm ! to enter field free region and calculate vector for field tracking program. vel = abs(mom)/sqrt(mom**2+mass**2)*cc eng = sign(1.,mom)*sqrt(mom**2+mass**2) mom_0 = mom/(1.d0+delta/100.d0) vT(1) = -fry + 100.*dx vT(2) = y + 100.*dy vT(3) = 100. vT(6) = vel/SQRT(1+dy**2+dx**2) vT(4) = dx*vT(6) vT(5) = dy*vT(6) call print_coord3('first track at z=100',vT) ! and track into the magnetic field to the beam plane (perp. to y) CALL trgTrackToPlane (vT,eng,1.,0.,-ctheta,stheta,frx,ok,spect) call print_coord3( 'first track on beam',vT) n = 0 delx = 1. DO WHILE ((delx .GT. .0001) .AND. (n .LT. 10) .AND. ok) delx = abs(-fry-vT(1)) if (debug(6)) write(*,*) '----------------' if (debug(6)) write(*,*) 'do while: ',n,delx ! track to the z=0 plane to find a correction for the x-offset vTx(1) = -fry DO i=2,6 vTx(i) = vT(i) ENDDO call print_coord3( 'vT, beam',vT) call print_coord3( 'vTx, beam',vTx) CALL trgTrackToPlane (vT, eng,1.,0.,0.,1.,0.,ok,spect) CALL trgTrackToPlane (vTx,eng,1.,0.,0.,1.,0.,ok,spect) call print_coord3( 'vT, z=0',vT) call print_coord3( 'vTx, z=0',vTx) xx = xx+min(1.,max(-1.,(vTx(1)-vT(1)))) xxd = xx ! now find a better approximation call print_coord1('before mc_hms_recon',y,delta,dx,dy,-xxd,0.) call mc_hms_recon(delta,dy,dx,y,xxd) mom = mom_0*(1.d0+delta/100.d0) vel = abs(mom)/sqrt(mom**2+mass**2)*cc eng = sign(1.,mom)*sqrt(mom**2+mass**2) call print_coord1('after mc_hms_recon',y,delta,dx,dy,-xxd,0.) ! drift to a field free region and calculate the velocities vT(1) = xx + 100.*dx vT(2) = y + 100.*dy vT(3) = 100. vT(6) = vel/SQRT(1+dy**2+dx**2) vT(4) = dx*vT(6) vT(5) = dy*vT(6) ! and track into the magnetic field to the beam plane (perp. to y) call print_coord3( 'before last track',vT) CALL trgTrackToPlane (vT,eng,1.,0.,-ctheta,stheta,frx,ok,spect) call print_coord3( 'after last track:',vT) n = n+1 ENDDO IF (delx .GT. .2) ok = .FALSE. ! calculate the result in HMS coordinates dy = vT(5)/vT(6) dx = vT(4)/vT(6) y = vT(2) return end ! routines to print coordinates given different input - GAW subroutine print_coord1(txt,y,delta,dxdz,dydz,x,z) implicit none include 'simulate.inc' character*(*) txt real*8 y,delta,dxdz,dydz,x,z real*8 vT(6) vt(1) = -x vt(2) = y vt(3) = z vt(6) = 30/sqrt(1.+dxdz**2+dydz**2) vt(4) = dxdz vt(5) = dydz if (debug(6)) write(*,100) txt,vt(1),vt(2),vt(3),vt(4), > vt(5),sqrt(vt(4)**2+vt(5)**2+vt(6)**2) 100 format(a20,6f10.4) return end subroutine print_coord2(txt,x,y,z,dxdz,dydz) implicit none include 'simulate.inc' character*(*) txt real*8 y,dxdz,dydz,x,z real*8 vT(6) vt(1) = x vt(2) = y vt(3) = z vt(6) = 30/sqrt(1.+dxdz**2+dydz**2) vt(4) = dxdz vt(5) = dydz if (debug(6)) write(*,100) txt,vt(1),vt(2),vt(3),vt(4), > vt(5),sqrt(vt(4)**2+vt(5)**2+vt(6)**2) 100 format(a20,6f10.4) return end subroutine print_coord3(txt,vt) implicit none include 'simulate.inc' character*(*) txt real*8 vT(6) if (debug(6)) write(*,100) txt,vt(1),vt(2),vt(3),vt(4)/vt(6), > vt(5)/vt(6),sqrt(vt(4)**2+vt(5)**2+vt(6)**2) 100 format(a20,6f10.4) return end