      SUBROUTINE GUKINE
C
C Generates Kinematics for primary track
C
      implicit none

C Include files

       include 'gcbank.inc'
       include 'gcflag.inc'
       include 'gconsp.inc'
       include 'mcflags.inc'
       include 'vols.inc'
       include 'part.inc'
       include 'ntuple.inc'
       include 'spec.inc'

C Local variables

      integer retry, iprim, nvert, nktrk, nktrk2, I
      real ran, angle, displacement, x, y, z
      real momentum_tot(5), pmom(3)

******************************************************************************

C Reset the event data

100   call resetevt
C
C Here put code to generate input particles and
C attach them to the GEANT structures

C Put vertex in middle of Cs source

      VERTEX(1)=0.
      VERTEX(2)=0.
      VERTEX(3)=-(2.*LIDPAR(3) + SSPAR(3))


C Random radius and angle produced within Cs source

      CALL RANLUX (ran, 1)
      displacement = ((CSPAR(2)**2.)*ran)**0.5
      

      CALL RANLUX (ran, 1)
      angle = ran*2.*Pi
      IF (angle.lt.Pi/2.) THEN
         VERTEX(1) = VERTEX(1) + displacement*cos(angle)
         VERTEX(2) = VERTEX(2) + displacement*sin(angle)
      ELSE IF (angle.lt.Pi) THEN
         angle = Pi - angle
         VERTEX(1) = VERTEX(1) - displacement*cos(angle)
         VERTEX(2) = VERTEX(2) + displacement*sin(angle)
      ELSE IF (angle.lt.3.*Pi/4.) THEN
         angle = angle - Pi
         VERTEX(1) = VERTEX(1) - displacement*cos(angle)
         VERTEX(2) = VERTEX(2) - displacement*sin(angle)
      ELSE 
         angle = 2.*Pi - angle
         VERTEX(1) = VERTEX(1) + displacement*cos(angle)
         VERTEX(2) = VERTEX(2) - displacement*sin(angle)
      END IF

      CALL RANLUX (ran,1)
      VERTEX(3) = VERTEX(3)+((2.*ran)-1.)*(CSPAR(3))

C Generate a random direction for the initial particles and assign
C momentum accordingly
      IF (PSWITCH .EQ. 1) THEN
         
         DO I=1,5
         CALL RANLUX (ran,1)
         z = 2.*ran-1.
         CALL RANLUX (ran,1)
         angle = ran * 2. * Pi
         IF (angle.lt.Pi/2.) THEN
            x = ((1-z**2.)**0.5)*cos(angle)
            y = ((1-z**2.)**0.5)*sin(angle)
         ELSE IF (angle.lt.Pi) THEN
            angle = Pi - angle
            x = -1.*((1-z**2.)**0.5)*cos(angle)
            y = ((1-z**2.)**0.5)*sin(angle)
         ELSE IF (angle.lt.3.*Pi/4.) THEN
            angle = angle - Pi
            x = -1.*((1-z**2.)**0.5)*cos(angle)
            y = -1.*((1-z**2.)**0.5)*sin(angle)
         ELSE 
            angle = 2.*Pi - angle
            x = ((1-z**2.)**0.5)*cos(angle)
            y = -1.*((1-z**2.)**0.5)*sin(angle)
         END IF

          momentum_tot(I) = (PX(I)**2.+PY(I)**2.+PZ(I)**2.)**0.5
          PX(I) = momentum_tot(I) * x
          PY(I) = momentum_tot(I) * y
          PZ(I) = momentum_tot(I) * z
         END DO
      ELSE IF (PSWITCH.EQ.2) THEN
C
C Generate a random direction
C
       CALL RANLUX (ran,1)
       z = 2.*ran-1.
       CALL RANLUX (ran,1)
       angle = ran * 2. * Pi
       IF (angle.lt.Pi/2.) THEN
          x = ((1-z**2.)**0.5)*cos(angle)
          y = ((1-z**2.)**0.5)*sin(angle)
       ELSE IF (angle.lt.Pi) THEN
          angle = Pi - angle
          x = -1.*((1-z**2.)**0.5)*cos(angle)
          y = ((1-z**2.)**0.5)*sin(angle)
       ELSE IF (angle.lt.3.*Pi/4.) THEN
          angle = angle - Pi
          x = -1.*((1-z**2.)**0.5)*cos(angle)
          y = -1.*((1-z**2.)**0.5)*sin(angle)
       ELSE
          angle = 2.*Pi - angle
          x = ((1-z**2.)**0.5)*cos(angle)
          y = -1.*((1-z**2.)**0.5)*sin(angle)
       END IF

       VERTEX(3) = -(2.*LIDPAR(3) + SSPAR(3)) + SSPAR(3) + 0.01
       call rnhran(beta_spectrum,550,0.,1.,momentum_tot(1))
       momentum_tot(1)=momentum_tot(1)/1E6
       px(1) = x*momentum_tot(1)
       py(1) = y*momentum_tot(1)
       pz(1) = z*momentum_tot(1)
       ptype(1) = 3 ! electron
      END IF

C This call 'UCTOH' takes a string 'PRIM' and converts it to an
C integer (IPRIM) which GEANT uses to keep track of what vertex it was
C There is an additional call 'UHTOC' which goes the other way around..

      CALL UCTOH('PRIM',IPRIM,4,4)

C Attach the VERTEX 'PRIM' to the event structure
C This returns the vertex id : NVERT 

      CALL GSVERT(VERTEX,0,0,IPRIM,1,NVERT)

C Now attach a photon and electron - with GEANT particle IDs of , 
C and momentum vector PMOMPH and PMOMEL to vertex NVERT
C This returns the value of NKTRK which is the number of tracks stored
C on this vertex
      NPART = 0
      DO I = 1,5
       IF (PTYPE(I) .EQ. 0) CYCLE
       NPART = NPART + 1
       pmom(1) = PX(I)
       pmom(2) = PY(I)
       pmom(3) = PZ(I)
       CALL GSKINE(pmom,PTYPE(I),NVERT,0,0,NKTRK)
      END DO

      CALL GPVERT(0)
      CALL GPKINE(0)

      WRITE(*,10) IRUN,IEVENT

10    FORMAT(6X,'--------- Run : ',I8,' Event : ',I8,'---------') 
      RETURN
      END
