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'