changed input file format to use simple 4-vectors

README.md

@@ -86,17 +86,17 @@ simc reads in your events from a plain text file. The format for this file is
 one line per event where the following HMS and SHMS variables are given.
 Note that at this point, we don't pass the vertex information yet **TODO**.
 ```bash
-hms_xptar hms_yptar hms_particle_mom shms_xptar shms_yptar shms_particle_mom weight
+px_HMS py_HMS pz_HMS E_HMS vz_HMS px_SHMS py_SHMS pz_SHMS E_SHMS vz_SHMS weight
 ```
 Definitions:
 ```
-hms_* :  HMS particle
-shms_*:  SHMS particle
-*_xptar: dx/dz relative to the central ray
-*_yptar: dy/dz relative to the central ray
-*_particle_momentum: Particle momentum in MeV
+*_HMS :  HMS particle
+*_SHMS:  SHMS particle
+px,py,pz: paricle 3-momentum components in GeV
+E: particle Energy in GeV
 weight: event weight
 ```
+Make sure you use the right coordinate system!
 The coordinate system is at the target center with z pointing along the central
 spectrometer angle, +x pointing vertical down and +y pointing to the left
 (i.e. smaller HMS angles and larger SHMS angles).

simc/event.f

@@ -300,14 +300,12 @@ C modified 5/15/06 for poinct
 	if ( doing_2pi) then
 	   ok_2pi = .true.
 	   call get_file_event(spec%e%theta,spec%p%theta,
-     >       vertex%e%xptar,vertex%e%yptar,vertex%e%p,
-     >       vertex%p%xptar,vertex%p%yptar,vertex%p%p,
+     >       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)
 	     if ( .not. ok_2pi) goto 100
 	     vertex%e%delta = 100.*(vertex%e%P-spec%e%P)/spec%e%P
 	     vertex%p%delta = 100.*(vertex%p%P-spec%p%P)/spec%p%P
-	     vertex%e%E = sqrt(vertex%e%P*vertex%e%P + 0.511*0.511)
-	     vertex%p%E = sqrt(vertex%p%P*vertex%p%P + 0.511*0.511)
        if(debug(5)) then
           write(*,*) ' E and Delta: '
           write(*,*) '     HMS: ',vertex%e%E,vertex%e%delta

simc/get_file_event.f

       subroutine get_file_event(th_spec_e,th_spec_p,
-     >    dxdz,dydz,e_mom,dxdzp,dydzp,p_mom,weight)
+     >    dxdz,dydz,e_mom,e_E,dxdzp,dydzp,p_mom,p_E,weight)
 c
 c  input variables:
 c        th_spec_e : central spec angle for electron (rad)
@@ -8,19 +8,25 @@ c  output variables:
 c        dxdz : xptar for electron
 c        dydz : yptar for electron
 c        e_mom : electron momentum ( MeV)
+c        e_E : electron Energy ( MeV)
 c        dxdzp : xptar for proton
 c        dydzp : yptar for proton
 c        p_mom : proton momentum ( MeV)
+c        p_E : proton Energy ( MeV)
+c        weight: event weight
+c
+c        NOTE: vertex also read in but currently ignored
 c
       implicit none
 	include 'simulate.inc'
 c
-         real*8 e_mom
-         real*8 p_mom
+         real*8 e_4v(4), e_4vr(4)
+         real*8 p_4v(4), p_4vr(4)
+         real*8 cth_hms, sth_hms, cth_shms, sth_shms
          real*8 th_spec_e
          real*8 th_spec_p
-         real*8 dxdz,dydz
-         real*8 dxdzp,dydzp
+         real*8 e_mom,e_E,dxdz,dydz,e_vz
+         real*8 p_mom,p_E,dxdzp,dydzp,p_vz
          real*8 w, weight
          character*80 multpifile
          integer count,count_miss
@@ -41,23 +47,63 @@ c
 c
 c
          end_of_2pi_file = .false. 
-         read(51,*,end=999,err=999) dxdz,dydz,e_mom,dxdzp,dydzp,p_mom,w
+         read(51,*,end=999,err=999) e_4v,e_vz,p_4v,p_vz,w
          count = count + 1
+	       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(*,*) '    vertex: ',e_vz,p_vz
+	       endif !debug
 c
-         p_mom = p_mom*1000.
+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
+         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(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(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(*,*) ' 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)
+	       endif !debug
+c Calculate dxdz and dydz
+         e_mom = sqrt(e_4vr(1)**2 + e_4vr(2)**2 + e_4vr(3)**2)
+         p_mom = sqrt(p_4vr(1)**2 + p_4vr(2)**2 + p_4vr(3)**2)
+         dxdz = e_4vr(1)/e_mom
+         dydz = e_4vr(2)/e_mom
+         dxdzp = p_4vr(1)/p_mom
+         dydzp = p_4vr(2)/p_mom
+c Convert momenta and energies into MeV
          e_mom = e_mom * 1000.
-         dxdz = dxdz/1.
-         dydz = dydz/1.
-
-         dxdzp = dxdzp/1.
-         dydzp = dydzp/1.
+         p_mom = p_mom * 1000.
+         e_E = e_4vr(4) * 1000.
+         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(*,*) ' NEW EVENT: ',count
-             write(*,*) '       HMS: ',dxdz,dydz,e_mom
-             write(*,*) '      SHMS: ',dxdzp,dydzp,p_mom
+             write(*,*) 'SIMC INPUT: '
+             write(*,*) '       HMS: ',dxdz,dydz,e_mom,e_E,e_vz
+             write(*,*) '      SHMS: ',dxdzp,dydzp,p_mom,p_E,p_vz
 	       endif !debug
 c        
          return