*------------------------------------------------------------------------ * TARGET tracking an electron through the Gen target field * -======- * * Raytracing of 3-d motion in polarized target field by solution * of differential equations of motion via 4th order Runge-Kutta. Field * orientation arbitrary. * - reads the start parameters form STDIN and * - writes the results to an PAW/hbook ntuple file * * written by Markus Muehlbauer *------------------------------------------------------------------------ IMPLICIT NONE INTEGER nwpawc,icycle,istat PARAMETER (nwpawc = 1000000) REAL paw(nwpawc) COMMON /pawc/paw CHARACTER*8 tags(9) DATA tags/'N','t','l','x','y','z','Vx','Vy','Vz'/ REAL u(6),vi(3),v0(3) LOGICAL uniform,recoils,backwards CHARACTER ans*1 REAL B_theta,B_phi REAL e0_e,ef_e,e_tot,beta,qvec,t_max,tmass REAL tht_e,tht_q,phi,phi_max,phi_step REAL Bdl,theta_bdl,theta_if,arg1,arg2,arg3,ratio INTEGER n_events,n_steps,i,j REAL c,pmass,dmass,pi180 PARAMETER (c = 30.0) ! speed of light (cm/nanosec) PARAMETER (pmass = 938.28) ! proton mass in MeV PARAMETER (dmass = 1875.6) ! deuteron mass in MeV PARAMETER (pi180 = 3.141592653/180.) INTEGER bookTrack EXTERNAL bookTrack uniform = .false. recoils = .true. backwards = .false. ! Ask questions WRITE (6,'(a)')'Gen 3D target raytraying:' WRITE (6,'(a)')'=========================' WRITE (6,'(a,$)')'Enter E_0 [GeV], theta_e [deg]> ' READ (5,*)e0_e,tht_e WRITE (6,'(a,$)')'Enter phi_max [deg], '// > 'and # of steps in phi_max> ' READ (5,*)phi_max,n_events WRITE (6,'(a,$)')'Enter t_max [ns] and # of time steps> ' READ (5,*)t_max,n_steps WRITE (6,'(a,$)')'Enter spherical angles (theta,phi) of B [deg]> ' READ (5,*)B_theta,B_phi WRITE (6,'(a,$)')'Uniform field test? (y/n)> ' READ (5,'(a)')ans IF (ans .eq. 'y' .OR. ans .EQ. 'Y') uniform = .true. WRITE (6,'(a,$)')'Track electrons or recoil particles? (e/r)> ' READ (5,'(a)')ans IF(ans .eq. 'e' .OR. ans .EQ. 'E') recoils = .false. WRITE (6,'(a,$)')'Is target proton or deuteron? (p/d) >' READ (5,'(a,$)')ans IF(ans .EQ. 'p' .OR. ans .EQ. 'P') THEN tmass = pmass ELSE tmass = dmass ENDIF IF (.not. recoils) THEN WRITE (6,'(a,$)')'Track electrons backwards? (y/n) >' READ (5,'(a)')ans IF (ans .EQ. 'y' .OR. ans .EQ. 'Y') backwards = .true. ENDIF WRITE(6,*) ! Use Oxford field map or do uniform field test case IF (uniform) THEN CALL trgInit (' ',B_theta,B_phi) ELSE CALL trgInit ('field_map.dat',B_theta,B_phi) ENDIF ! calculate kinematic quantities e0_e = e0_e*1000. !to MeV ef_e = e0_e/(1.+2.*e0_e*SIN(0.5*tht_e*pi180)**2/tmass) qvec = SQRT((e0_e-ef_e)*2.*tmass + (e0_e-ef_e)**2) tht_q = -ASIN(ef_e*SIN(tht_e*pi180)/qvec)/pi180 phi_step = 2.*phi_max/n_events ! set up initial "velocity" vector IF (recoils) THEN ! tracking recoils (protons or deuterons) e_tot = sqrt(qvec**2 + tmass**2) beta = qvec/e_tot vi(1) = beta*c*SIN(tht_q*pi180) vi(2) = beta*c*SIN(tht_q*pi180) vi(3) = beta*c*COS(tht_q*pi180) ELSEIF (backwards) THEN ! tracking incident beam beta = 1.0 e_tot = -e0_e vi(1) = 0. vi(2) = 0. vi(3) = -c ELSE ! tracking scattered e- beta = 1.0 e_tot = -ef_e vi(1) = c*SIN(tht_e*pi180) vi(2) = c*SIN(tht_e*pi180) vi(3) = c*COS(tht_e*pi180) ENDIF OPEN (unit=10,file='bdl.out',status='unknown') CALL hlimit(nwpawc) CALL hropen(20,'Tracks','tracks.hbook','NX',4096,ISTAT) CALL hbookn(1,'Particle Tracks',9,'//Tracks',4096,tags) DO i=1,n_events !azimuthal angle loop phi = -phi_max+phi_step*i ! set up initial conditions, including actual velocity vector u(1) = 0. ! x u(2) = 0. ! z u(3) = 0. ! y u(4) = vi(1)*SIN(phi*pi180) ! dx/dt u(5) = vi(2)*COS(phi*pi180) ! dy/dt u(6) = vi(3) ! dz/dt DO j=1,3 ! save initial velocity v0(j) = u(3+j) ENDDO ! track particle CALL trgXTrack (u,e_tot,t_max*30./n_steps,t_max*30.,Bdl, > bookTrack,i) ! calculate total deflection angle based on integral B-dl theta_Bdl = abs(171.887*Bdl/e_tot/beta) ! also calculate it from v_i-dot-v_f arg1 = v0(1)*u(4) + v0(2)*u(5) + v0(3)*u(6) arg2 = v0(1)**2 + v0(2)**2 + v0(3)**2 arg3 = u(4)**2 + u(5)**2 + u(6)**2 theta_if = ACOS(arg1/SQRT(arg2*arg3))/pi180 ratio = 0 IF (theta_if .NE. 0.) ratio = theta_bdl/theta_if WRITE (10,1000) i,Bdl,theta_if,theta_Bdl,ratio ENDDO !end of 'phi' loop CALL hrout(0,ICYCLE,' ') CALL hrend('Tracks') 1000 FORMAT(/,t10,'Track #: ',i3,t40,'Integral B_dl = ',f8.3, > ' T-cm',/,t10,'angle(v_i,v_f) = ',f9.3, > t40,'angle_Bdl = ',f9.3,' degrees',/, > t10,'Bdl/v_if = ',f9.4) STOP END INTEGER FUNCTION bookTrack (id,t,l,u) IMPLICIT NONE INTEGER id REAL t,l,u(9) * -- books the particle location in the n-tuple file REAL event(9) INTEGER i event (1) = id event (2) = t event (3) = l DO i=1,6 event (i+3) = u(i) ENDDO CALL hfn(1,event) bookTrack = id RETURN END