* simulation of a1p experiment
*

      subroutine gen_evt

      implicit none

      include 'constants.inc'
      include 'beta_geom.inc'
      include 'sane_misc.inc'
      include 'sane.inc'

      real*8 u(6),u_vertex(6)

      integer*2 k
      integer*2 tgt
      logical ok
      real*8 E
      logical*1 check_accp
      real*8 tht_scat,phi_scat
      real*8 cross_section,calc_asym
      integer flg
      data flg /1/

C Now cycle through events
      if(flg.eq.1) then
      write(*,*) 'E_beam',E_beam,'ebeam2 = ',ebeam2
      E_beam =  ebeam2
      write(*,*) 'E_beam',E_beam,'ebeam2 = ',ebeam2
      flg = 0
      endif
      
 700  continue

      tgt = nint((num_tgts-1)*rand())+1
      zz = z(tgt)
      nn = n(tgt)

      call generate_event(u_vertex,E,part,tht_scat,phi_scat)

      xsn = cross_section(zz,nn,u_vertex,part)
      pp = sqrt(E**2-mass(part)**2)
      th = tht_scat
      ph = phi_scat
      EE = E

      if (part.eq.1.and.xsn.gt.0.d0) then
        asym = calc_asym(E_beam/1.d3,abs(E),tht_scat/d2r,theta_Bfield)
      else
        asym = 0.d0
      endif

      do k=1,6
        uu(k) = u_vertex(k)
      enddo

      return
      end

C------------------------------------------------------------------------------

      subroutine generate_event(u,E,part,tht,phi)

*       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 : points toward calorimeter from target
*
*             - 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]
* 
*       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 [GeV] * sign of particle charge
*                (negative for electrons, positive for protons/deuterons)
*
* Beam coordinate system
*    x points down
*    y points to beam left
*    z points down beam (downstream)
* this coordinate system is then simply related to the spectrometer 
* coordinate systems by a rotation about x, but it is different from
* the more typical where x points horizontal.

      implicit none

      include 'sane_accp.inc'
      include 'constants.inc'
      include 'sane.inc'

      real*8 c_lo,c_hi
      real*8 p,E,phi,c
      real*8 vel,u(6),v(6),tht,tht0
      real*8 cell_rad,cell_phi,cell_z
      integer*2 i
      integer*4 part


      tht0 = theta_0*d2r

      c_lo = cos(theta_max*d2r)
      c_hi = cos(theta_min*d2r)

C Throw momentum, phi and cos(theta)
          
      p = p_min + (p_max - p_min)*rand()
      phi = phi_min + (phi_max - phi_min)*rand()
      c = c_lo + (c_hi - c_lo)*rand()
      tht = acos(c)

      E = sqrt(p**2+mass(part)**2)
      vel = cc*p/E
      if (part_charge(part).ne.0) E = E*part_charge(part)

C     Generate vertex in target cell.  Assume distribution is a cylinder
C     of legnth cell_length and radius raster_radius.

      cell_z = cell_length*(rand()-0.5d0)
      cell_rad = raster_radius*sqrt(rand())  ! throw radius squared
      cell_phi = 360.d0*rand()*d2r
C      cell_z = 0.d0
C      cell_rad = 0.d0
C      cell_phi = 0.d0

      v(1) = cell_rad*cos(cell_phi)
      v(2) = cell_rad*sin(cell_phi)
      v(3) = cell_z
      
* generate angles in beam coordinate system

      v(4) = vel*sin(tht)*cos(phi*d2r)
      v(5) = vel*sin(tht)*sin(phi*d2r)
      v(6) = vel*cos(tht)

* Now, rotate angles so that z is pointing in direction of calorimeter

      u(1) = v(1)
      u(2) = cos(tht0)*v(2) - sin(tht0)*v(3)
      u(3) = sin(tht0)*v(2) + cos(tht0)*v(3)

      u(4) = v(4)
      u(5) = cos(tht0)*v(5) - sin(tht0)*v(6)
      u(6) = sin(tht0)*v(5) + cos(tht0)*v(6)

C      write(*,*) 'unrotated: ',v
C      write(*,*) '  rotated: ',u
      return
      end

C------------------------------------------------------------------------------

      logical*1 function check_accp(u)

      implicit none

      include 'beta_geom.inc'
      
      real*8 u(6)

      if ((abs(u(1)).le.cal_height/2.).and.
     1    (abs(u(2)).le.cal_width/2.)       ) then
        check_accp = .true.
      else
        check_accp = .false.
C        check_accp = .true.
      endif

      return
      end

C------------------------------------------------------------------------------

      subroutine write_setup_file(nevts)

      implicit none

      include 'sane.inc'
      include 'beta_geom.inc'
      include 'sane_accp.inc'

      integer*4 nevts

      open(29,file='a1p_model_setup.dat',status='unknown')

      write(29,'(f8.1,t30,"! ",a20)') E_beam,      "Beam Energy"
      write(29,'(f8.1,t30,"! ",a20)') theta_0,     "Theta_e"
      write(29,'(f8.1,t30,"! ",a20)') theta_Bfield,"Theta_B"
      write(29,'(f8.1,t30,"! ",a20)') theta_min,   "Phi Maximum"
      write(29,'(f8.1,t30,"! ",a20)') theta_max,   "Theta Minimum"
      write(29,'(f8.1,t30,"! ",a20)') phi_min,     "Phi Minimum"
      write(29,'(f8.1,t30,"! ",a20)') phi_max,     "Phi Maximum"
      write(29,'(f9.2,t30,"! ",a20)') p_min,       "Momentum Minimum"
      write(29,'(f9.2,t30,"! ",a20)') p_max,       "Momentum Maximum"
      write(29,'(f8.1,t30,"! ",a20)') cal_width,   "Cal Width"
      write(29,'(f8.1,t30,"! ",a20)') cal_height,  "Cal Height"
      write(29,'(f8.1,t30,"! ",a20)') cal_drift,   "Cal Drift"
      write(29,'(i8,t30,"! ",a20)')   nevts,       "Num of thrown events"

      close(29)

      return
      end
