diff --git a/simc/event.f b/simc/event.f index f425a33900e0d7645c374b10f01d665ec9cd8952..40f6c4c901102ad4e962a5444423cdf487740062 100644 --- a/simc/event.f +++ b/simc/event.f @@ -299,7 +299,7 @@ C modified 5/15/06 for poinct endif !not (doing_hyd_elast) if ( doing_2pi) then ok_2pi = .true. - call get_file_event(spec%e%theta,spec%p%theta, + call get_file_event(electron_arm,spec%e%theta,spec%p%theta, > vertex%e%xptar,vertex%e%yptar,vertex%e%p,vertex%e%E, > vertex%p%xptar,vertex%p%yptar,vertex%p%p,vertex%p%E, > main%gen_weight) diff --git a/simc/get_file_event.f b/simc/get_file_event.f index 765da16aa48aacbe93266bb38660464c2ce276c6..1eceead3c33146ff5ae84c2552fe057b2176705e 100644 --- a/simc/get_file_event.f +++ b/simc/get_file_event.f @@ -1,7 +1,9 @@ - subroutine get_file_event(th_spec_e,th_spec_p, + subroutine get_file_event(e_arm,th_spec_e,th_spec_p, > dxdz,dydz,e_mom,e_E,dxdzp,dydzp,p_mom,p_E,weight) c c input variables: +c electron_arm : 1 means HMS is electron arm +c 5 means SHMS is the electron arm c th_spec_e : central spec angle for electron (rad) c th_spec_p : central spec angle for proton (rad) c output variables: @@ -20,12 +22,14 @@ c implicit none include 'simulate.inc' c - real*8 e_4v(4), e_4vr(4) - real*8 p_4v(4), p_4vr(4) + integer*4 e_arm,ii + real*8 SHMS_4v(4), e_4vr(4), e_4v(4) + real*8 HMS_4v(4), p_4vr(4), p_4v(4) real*8 cth_hms, sth_hms, cth_shms, sth_shms real*8 th_spec_e real*8 th_spec_p - real*8 e_mom,e_E,dxdz,dydz,e_vz + real*8 sth_elec,cth_elec,sth_prot,cth_prot + real*8 e_mom,e_E,dxdz,dydz,e_vz,SHMS_vz,HMS_vz real*8 p_mom,p_E,dxdzp,dydzp,p_vz real*8 w, weight character*80 multpifile @@ -47,40 +51,69 @@ c c c end_of_2pi_file = .false. - read(51,*,end=999,err=999) e_4v,e_vz,p_4v,p_vz,w + 888 continue + read(51,*,end=999,err=999) HMS_4v,HMS_vz,SHMS_4v,SHMS_vz,w count = count + 1 +c if (count .gt. 100) goto 999 + +c Rotate 4-vectors to their respective frames +c This is a rotation around the x-axis, which points downwards +c --> the HMS is at negative angles, and the SHMS at positive angles + if (e_arm .eq. 1) then + cth_elec = cos(th_spec_e) + sth_elec = sin(-th_spec_e) + cth_prot= cos(th_spec_p) + sth_prot = sin(th_spec_p) + e_vz = HMS_vz * 1. + p_vz = SHMS_vz * 1. + e_4v(1) = HMS_4v(1) + p_4v(1) = SHMS_4v(1) + e_4v(2) = HMS_4v(2) + p_4v(2) = SHMS_4v(2) + e_4v(3) = HMS_4v(3) + p_4v(3) = SHMS_4v(3) + e_4v(4) = HMS_4v(4) + p_4v(4) = SHMS_4v(4) + else + cth_prot = cos(th_spec_p) + sth_prot = sin(-th_spec_p) + cth_elec = cos(th_spec_e) + sth_elec = sin(th_spec_e) + e_vz = SHMS_vz * 1. + p_vz = HMS_vz * 1. + e_4v(1) = SHMS_4v(1) + p_4v(1) = HMS_4v(1) + e_4v(2) = SHMS_4v(2) + p_4v(2) = HMS_4v(2) + e_4v(3) = SHMS_4v(3) + p_4v(3) = HMS_4v(3) + e_4v(4) = SHMS_4v(4) + p_4v(4) = HMS_4v(4) + endif if(debug(5)) then write(*,*) ' ' write(*,*) ' NEW EVENT: ',count, weight - write(*,*) ' HMS: ',e_4v(1),e_4v(2),e_4v(3),e_4v(4) - write(*,*) ' SHMS: ',p_4v(1),p_4v(2),p_4v(3),p_4v(4) + write(*,*) ' e: ',e_4v(1),e_4v(2),e_4v(3),e_4v(4) + write(*,*) ' p: ',p_4v(1),p_4v(2),p_4v(3),p_4v(4) write(*,*) ' vertex: ',e_vz,p_vz endif !debug -c -c Rotate 4-vectors to their respective frames -c This is a rotation around the x-axis, which points downwards -c --> the HMS is at negative angles, and the SHMS at positive angles - cth_hms = cos(th_spec_e) - sth_hms = sin(-th_spec_e) - cth_shms = cos(th_spec_p) - sth_shms = sin(th_spec_p) -c Rotatation about the x-axis --> only y, and z change +cc Rotatation about the x-axis --> only y, and z change e_4vr(1) = e_4v(1) - e_4vr(2) = cth_hms * e_4v(2) - sth_hms * e_4v(3) - e_4vr(3) = sth_hms * e_4v(2) + cth_hms * e_4v(3) + e_4vr(2) = cth_elec * e_4v(2) - sth_elec * e_4v(3) + e_4vr(3) = sth_elec * e_4v(2) + cth_elec * e_4v(3) e_4vr(4) = e_4v(4) p_4vr(1) = p_4v(1) - p_4vr(2) = cth_shms * p_4v(2) - sth_shms * p_4v(3) - p_4vr(3) = sth_shms * p_4v(2) + cth_shms * p_4v(3) + p_4vr(2) = cth_prot * p_4v(2) - sth_prot * p_4v(3) + p_4vr(3) = sth_prot * p_4v(2) + cth_prot * p_4v(3) p_4vr(4) = p_4v(4) if(debug(5)) then write(*,*) ' ' write(*,*) ' SPEC:' - write(*,*) ' HMS: ',-th_spec_e/pi*180 - write(*,*) ' SHMS: ',th_spec_p/pi*180 + write(*,*) ' e: ',th_spec_e/pi*180 + write(*,*) ' p: ',th_spec_p/pi*180 write(*,*) ' ROTATED: ' - write(*,*) ' HMS: ',e_4vr(1),e_4vr(2),e_4vr(3),e_4vr(4) - write(*,*) ' SHMS: ',p_4vr(1),p_4vr(2),p_4vr(3),p_4vr(4) + write(*,*) ' e: ',e_4vr(1),e_4vr(2),e_4vr(3),e_4vr(4) + write(*,*) ' p: ',p_4vr(1),p_4vr(2),p_4vr(3),p_4vr(4) endif !debug c Calculate dxdz and dydz e_mom = sqrt(e_4vr(1)**2 + e_4vr(2)**2 + e_4vr(3)**2) @@ -96,16 +129,18 @@ c Convert momenta and energies into MeV p_E = p_4vr(4) * 1000. c Weight and vertex positions are read in correctly. c Vertez z is read in cm - e_vz = e_vz * 1. - p_vz = p_vz * 1. weight = w if(debug(5)) then write(*,*) ' ' write(*,*) 'SIMC INPUT: ' - write(*,*) ' HMS: ',dxdz,dydz,e_mom,e_E,e_vz - write(*,*) ' SHMS: ',dxdzp,dydzp,p_mom,p_E,p_vz + write(*,*) ' e: ',dxdz,dydz,e_mom,e_E,e_vz + write(*,*) ' p: ',dxdzp,dydzp,p_mom,p_E,p_vz endif !debug c + if ( abs(atan(e_4v(1)/e_4v(3))*57.3-th_spec_e)<3. + c .and. abs(atan(p_4v(1)/p_4v(3))*57.3-th_spec_p)<3. ) then + count_miss = count_miss+1 + endif return c 999 write(*,*) ' reached end of file'