      SUBROUTINE GUSTEP
      IMPLICIT NONE

       include 'gcsets.inc'
       include 'gcking.inc'
       include 'gckine.inc'
       include 'gcflag.inc'
       include 'gctrak.inc'
       include 'gctmed.inc'
       include 'mcflags.inc'
       include 'sets.inc'
       include 'gcvolu.inc'
       include 'ntuple.inc'

      INTEGER neutrino_code, trTypeGeantino, i
      PARAMETER (neutrino_code = 4)
      PARAMETER (trTypeGeantino = 6)
      CHARACTER*4 vname
      LOGICAL inHOLD, inSCIN


      INTEGER historyHOLD, historyWRLD, historyNEW
      INTEGER historySCIN
      PARAMETER (historyNEW = 0)
      PARAMETER (historyHOLD = 1)
      PARAMETER (historySCIN = 2)
      PARAMETER (historyWRLD = 3)
      INTEGER sameVolume, leavingVolume, enteringVolume
      PARAMETER (sameVolume = 0)
      PARAMETER (enteringVolume = 1)
      PARAMETER (leavingVolume = 2)

C Subroutine to handle how we deal with secondaries in this 
C tracking step for this particular particle

      CALL doSecondaries

C If the particle is a neutrino or is unspecified, then just ignore
C IPart if contained in the common block gckine.inc

      IF (iPart .EQ. neutrino_code .OR. 
     +    iTrTyp.EQ.trTypeGeantino) RETURN

C If we are in the scintillator volume (Set=1) then make a hit
C IsVol is contained in the common block gctmed.inc - it is 1 if we are in a sensitive volume
C and 0 if not. We only record hits for sensitive volumes
C
C ISet is contained in the common block gcsets.inc - identifies what set the current sensitive
C detector is in

      IF (isVol.eq.1 .and. iSet.EQ.iSetScint) THEN
       CALL doScintHit
       N_S_HITS = N_S_HITS + 1
      ENDIF

      IF (isVol.eq.1 .and. iSet.EQ.iSetBody) THEN
       CALL doBodyHit
       N_B_HITS = N_B_HITS + 1
      ENDIF

C
C Routine for identifying if particle is leaving holder or entering
C scintillator
C
      inHOLD = .false.
      inSCIN = .false.
      DO i = 1,nlevel
         CALL UHTOC(NAMES(i),4,vname,4)
         IF (vname.eq.'HOLD') inHOLD = .true.
         IF (vname.eq.'SCIN') inSCIN = .true.
      END DO
      IF (istory .NE. historyNEW) THEN
         IF ((inwvol.EQ.enteringVolume) .AND. 
     +  .NOT. inHOLD .AND.
     + (istory.EQ.historyHOLD)) CALL doExitHold
         IF ((inwvol.EQ.enteringVolume) .AND. inSCIN .AND.
     + (istory.NE.historySCIN)) CALL doEnterScin
      END IF


C Updates particle history for next run
      IF (istory .EQ. historyNEW .OR. 
     +    inwvol .NE. sameVolume) THEN
        IF (inHOLD) THEN
          istory = historyHOLD
        ELSE IF (inSCIN) THEN
          istory = historySCIN
        ELSE
          istory = historyWRLD
        END IF
      END IF

C Store space points for debugging information

      CALL GSXYZ
      CALL GDEBUG
      RETURN
      END
C-------------------------------------------------------------------------------
      SUBROUTINE doSecondaries
      IMPLICIT NONE
       include 'gcking.inc'
       include 'gckine.inc'
       include 'gcbank.inc'
      LOGICAL generatingVertex
      INTEGER jk, iOrigVtx, i, j, iAdrDum
      CHARACTER*4 mechanism

C Ngkine is > 0 if secondary particle have been generated,
C store them in the particle stack
C Secondary neutrinos are not stored

      if (ngkine .GT. 0) then

          generatingVertex = .FALSE.
          call uhtoc(kcase,4,mechanism,4)
 
C Loop over all secondary particles from this interaction
 
          do i = 1,ngkine
              j = gkin(5,i)                !Particle type
              if (j.eq.4) then
C  Ignore secondary neutrinos
                  iflgk(i) = -1
              else if (mechanism .EQ. 'DCAY') then
C  Store all other decay products
                  iflgk(i) = 1
              else if((mechanism .EQ. 'HADR' .OR.
     &                 mechanism .EQ. 'MUNU' .OR. 
     +                 mechanism.EQ.'EINC' .OR.
     &                 mechanism .EQ. 'ECOH' .OR. 
     +                 mechanism.EQ.'ANNH' .OR.
     &                 mechanism .EQ. 'INHE') .AND. (j.NE.1)) then
C  Load secondary hadrons (but not photons) into the main particle (KINE) list
                  iflgk(i) = 1
              else
C  Load shower products onto the temporary stack; ITRA remains the KINE particle number
                  iflgk(i) = 0
              endif
              if (iflgk(i).EQ.1) generatingVertex = .TRUE.
          enddo

C Load everything onto the stack

          call gsking(0)

C If we are adding a new vertex to JKINE, then add the mechanism that created
C this vertex (kcase) to the vertex datastructure

          if (generatingVertex .and. iflgk(ngkine+1).gt.0) then
              iadrDum = 0
              call gsveru(iflgk(ngkine+1),1,kcase,iadrDum)
          endif
      ENDIF


      RETURN
      END
C ----------------------------------------------------------------------------------------
      SUBROUTINE doScintHit
      IMPLICIT NONE
       include 'gcsets.inc'
       include 'gckine.inc'
       include 'gctrak.inc'
       include 'mcflags.inc'
       include 'gctmed.inc'
       include 'gcvolu.inc'

C Make a scintillator hit - basically this means that I record some
C information about the hit in an array and store it in the HIT data 
C structure.

      REAL HIT(13)
      INTEGER IHIT

C If no energy has been lost then there wasn't a hit...
 
      IF (DESTEP.LE.0) RETURN
C
C Define the hit : Usually this will just be the position (X,Y,Z), the
C track identifier (ITRA), the time of the hit (TOFG), the momentum
C and energy of the particle after the hit (PX,PY,PZ,E), and the amount
C of energy deposited in the step, DESTEP
C
C All the variables (ISET,IDET,NUMBV,IHIT) are already defined in
C the common block gcsets
C
      HIT(1) = VECT(1)   ! X
      HIT(2) = VECT(2)   ! Y
      HIT(3) = VECT(3)   ! Z
      HIT(4) = ITRA      ! Track id
      HIT(5) = IPART     ! Time
      HIT(6) = VECT(4)   ! P_X
      HIT(7) = VECT(5)   ! P_Y
      HIT(8) = VECT(6)   ! P_Z
      HIT(9) = VECT(7)   ! E
      HIT(10) = DESTEP   ! Delta(E)
      CALL GSAHIT(ISET,IDET,ITRA,NUMBV,HIT,IHIT)


      RETURN
      END

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

      SUBROUTINE doBodyHit
      IMPLICIT NONE
       include 'gcsets.inc'
       include 'gckine.inc'
       include 'gctrak.inc'
       include 'mcflags.inc'
       include 'gctmed.inc'
       include 'gcvolu.inc'

C Make a body hit - basically this means that I record some
C information about the hit in an array and store it in the HIT data 
C structure.

      REAL HIT(13)
      INTEGER IHIT

C If no energy has been lost then there wasn't a hit...
 
      IF (DESTEP.LE.0) RETURN
C
C Define the hit : Usually this will just be the position (X,Y,Z), the
C track identifier (ITRA), the time of the hit (TOFG), the momentum
C and energy of the particle after the hit (PX,PY,PZ,E), and the amount
C of energy deposited in the step, DESTEP
C
C All the variables (ISET,IDET,NUMBV,IHIT) are already defined in
C the common block gcsets
C
      HIT(1) = VECT(1)   ! X
      HIT(2) = VECT(2)   ! Y
      HIT(3) = VECT(3)   ! Z
      HIT(4) = ITRA      ! Track id
      HIT(5) = IPART     ! Time
      HIT(6) = VECT(4)   ! P_X
      HIT(7) = VECT(5)   ! P_Y
      HIT(8) = VECT(6)   ! P_Z
      HIT(9) = VECT(7)   ! E
      HIT(10) = DESTEP   ! Delta(E)
      CALL GSAHIT(ISET,IDET,ITRA,NUMBV,HIT,IHIT)

      RETURN
      END

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

      SUBROUTINE doExitHold
      IMPLICIT NONE

       include 'gcsets.inc'
       include 'gckine.inc'
       include 'gctrak.inc'
       include 'mcflags.inc'
       include 'gctmed.inc'
       include 'gcvolu.inc'
       include 'ntuple.inc'

C Gets data in the event of particles leaving the holder

      INTEGER geant2pdg
      EXTERNAL geant2pdg

      NESCAPE=NESCAPE+1
      IF (NESCAPE.LT.100)THEN
C         NTESCAPE(NESCAPE) = ITRA
         PTYP_ESC(NESCAPE) = geant2PDG(IPART)
         X_ESC(NESCAPE) = vect(1)
         Y_ESC(NESCAPE) = vect(2)
         Z_ESC(NESCAPE) = vect(3)
C         NPEX(NESCAPE) = vect(7)*vect(4)
C         NPEY(NESCAPE) = vect(7)*vect(5)
C         NPEZ(NESCAPE) = vect(7)*vect(6)
         PE_ESC(NESCAPE) = getot
      END IF

      RETURN
      END

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

      SUBROUTINE doEnterScin
      IMPLICIT NONE
       include 'gcsets.inc'
       include 'gckine.inc'
       include 'gctrak.inc'
       include 'mcflags.inc'
       include 'gctmed.inc'
       include 'gcvolu.inc'
       include 'ntuple.inc'

C Gets data in the event of particles entering the scintillator

      INTEGER geant2pdg
      EXTERNAL geant2pdg

      N_ENT_SCINT=N_ENT_SCINT+1
      IF (N_ENT_SCINT.LT.100)THEN
C         NTESCAPE(NESCAPE) = ITRA
         PTYP_ENT_SCINT(N_ENT_SCINT) = geant2PDG(IPART)
         X_ENT_SCINT(N_ENT_SCINT) = vect(1)
         Y_ENT_SCINT(N_ENT_SCINT) = vect(2)
         Z_ENT_SCINT(N_ENT_SCINT) = vect(3)
C         NPEX(NESCAPE) = vect(7)*vect(4)
C         NPEY(NESCAPE) = vect(7)*vect(5)
C         NPEZ(NESCAPE) = vect(7)*vect(6)
         PE_ENT_SCINT(N_ENT_SCINT) = getot
      END IF

      RETURN
      END
