      SUBROUTINE GUSTEP
C
C  $Id: gustep.F,v 1.21 2001/03/17 02:47:00 vavilov Exp $
C  $Author: vavilov $
C
C      Original Author:  Erik Ramberg,  April 1998
C
C  User routine called at the end of each tracking step
C
C  INWVOL is different from 0 when the track has reached
C       a volume boundary (=2 when exiting)
C  ISTOP is different from 0 if the track has stopped
C  NGKINE is different from 0 if secondaries generated
C

      Implicit none

#include "gtnum.inc"
#include "gcnum.inc"
#include "gttmed.inc"
#include "gctmed.inc"
#include "gtkine.inc"
#include "gckine.inc"
#include "gcking.inc"
#include "gtflag.inc"
#include "gcflag.inc"
#include "gttrak.inc"
#include "gctrak.inc"
#include "gtvolu.inc"
#include "gcvolu.inc"
#include "gtcuts.inc"
#include "gccuts.inc"
#include "gconsp.inc"
#include "gtscan.inc"
#include "gcscan.inc"
#include "gtlink.inc"
#include "gtbank.inc"
#include "gcbank.inc"

#include "ckm_para.inc"

C  Local Variables

      Logical LFirst,decay,first_kill,kill_track(10),kill_event
      Integer I, IT, IFL,j,mecdecay,norgpart,meccherenkov
      Integer nvert,nubuf,cherenkov_on_stack
#ifndef LINUX
      Integer ipart
#endif
      Real lvert(3),lpvert(4),ubuf(10)
      Integer lipart
      Data LFirst,First_Kill/.True.,.True./
      Save LFirst,mecdecay,norgpart,kill_track,first_kill,meccherenkov

      Character*4 NAMEC_TMP(MAXMEC)
      Equivalence (NAMEC_TMP,NAMEC)
      Character*4 NAMEC1_TMP(MAXME1)
      Equivalence (NAMEC1_TMP,NAMEC1)
      Character*4 NAMES_TMP(15)
      Equivalence (NAMES_TMP,NAMES)

      CHARACTER*4 CHCASE
      LOGICAL     LHADR, LDCAY, LCKOV
      INTEGER     IUCOMP, k
      REAL        MARK, EMAX, UPWGHT_SAVE
      PARAMETER  (MARK=99.9)

c      Write(*,*) 'In GUSTEP'

      If (LFirst) Then
        Do i = 1,MAXMEC ! find position of DCAY word in mechanism
          If (NAMEC_TMP(i).eq.'DCAY') 
     &         mecdecay = i
        EndDo
        Do i = 1,MAXME1 ! find position of DCAY word in mechanism
          If (NAMEC1_TMP(i).eq.'CKOV') 
     &         meccherenkov = 100+i
        EndDo
        CALL HBOOK1(1, 'N SEC', 50, -0.5, 49.5, 0.)
        CALL HBOOK1(2, 'N K+ ', 10, -0.5,  9.5, 0.)
        CALL HBOOK1(3, 'P KSEC', 100, 0., 25., 0.)
        LFirst = .False.
      EndIf

C
C         If SCAN is active...
C
      IF (SCANFL) THEN
          CALL GSCANU
          RETURN
      END IF

C
C  Put secondaries in stack 
C  Only save on gckine secondaries which are decay products of originals
C

      CALL UHTOC(KCASE,4,CHCASE,4)
      IF(ITRA.EQ.1.AND.ISTAK.EQ.0) UPWGHT=MARK

      IF(ITRTYP.NE.7 .AND. NGKINE.GT.0) THEN
         LCKOV = IUCOMP(105,LMEC,NMEC).NE.0
         LDCAY = IUCOMP(  5,LMEC,NMEC).NE.0
         LHADR = IUCOMP( 12,LMEC,NMEC).NE.0
         
c        IF(UPWGHT.EQ.MARK.AND.LDCAY.AND.NGKINE.NE.3) CALL DVDEBUG(3)

         IF(UPWGHT.EQ.MARK .AND. ISTOP.EQ.1 .AND. LDCAY) THEN
            DO i = 1, NGKINE
               IFLGK(i) = 1
               UPWGHT = 1.0
               CALL GSKING(i)

               NUBUF = 2
               UBUF(1) = ITRA
               UBUF(2) = tofg
               J=0
               CALL GSKINU(IFLGK(i),NUBUF,UBUF,J) ! add user info
            ENDDO
            UPWGHT = MARK
         ELSE

            IF(UPWGHT.EQ.MARK .AND. ISTOP.EQ.1) THEN
c
c     Initial particle stopped. Why ?
c     Reasonably there can be two possibilities 

               IF    (LCKOV.AND.(  ! 1) Artificial stop by CKOV process
     +                CHCASE.EQ.'NULL' .AND. NGKINE.EQ.1 .OR.
     +                CHCASE.EQ.'DRAY' .AND. NGKINE.EQ.2 .OR.
     +                CHCASE.EQ.'HADR' .AND. NGKINE.EQ.1 .AND.
     +                                       IPART.EQ.GKIN(5,1))) THEN

                  DO i = 1, NGKINE
                     IFLGK(i) = 0
                     UPWGHT = 1.0
                     IF(GKIN(5,i).EQ.IPART) UPWGHT = MARK
                     CALL GSKING(i)
                  ENDDO
                  UPWGHT = MARK
               ELSEIF(LHADR) THEN ! 2) Hadronic interatcion

                  IF(LCKOV.AND.NGKINE.EQ.1) CALL DVDEBUG(2)

                  k = 0
                  EMAX = -1.0 ! Find the fastest K
                  DO i = 1, NGKINE
                     IF(GKIN(5,i).EQ.IPART.AND.GKIN(4,i).GT.EMAX)THEN
                        EMAX = GKIN(4,i)
                        k = i
                     ENDIF
                  ENDDO

                  j = 0
                  CALL HF1(1, REAL(NGKINE), 1.)
                  DO i = 1, NGKINE
                     IFLGK(i) = 0
                     UPWGHT = 1.0
                     IF(GKIN(5,i).EQ.IPART) THEN
                        j = j + 1
                        CALL HF1(3,SQRT(GKIN(1,i)**2+GKIN(2,i)**2+
     +                                             GKIN(3,i)**2),1.)
                        IF(i.EQ.k.AND.GKIN(4,i).GT.ISWIT(7)) UPWGHT=MARK
                     ENDIF
                     CALL GSKING(i)
                  ENDDO
                  CALL HF1(2, REAL(j), 1.)
                  UPWGHT = MARK
               ELSE ! 3) Check if we cover all possibilities
                  CALL DVDEBUG(3)
               ENDIF

            ELSE ! Not an initial particle or not stop
               UPWGHT_SAVE = UPWGHT
               DO i = 1, NGKINE
                  IFLGK(i) = 0
                  UPWGHT =  1.0
                  CALL GSKING(i)
               ENDDO
               UPWGHT = UPWGHT_SAVE
            ENDIF ! Not an initial particle

         ENDIF ! Decay of an initial particle

      ENDIF ! There are secondaries

      cherenkov_on_stack = -1
      IF(ITRTYP.NE.7) THEN
        If (NGKINE.eq.1 .and. 
     &   (NAMES_TMP(NLEVEL).EQ.'KRGA' .or. NAMES_TMP(NLEVEL).EQ.'PRGA')
     &       ) Then
          If (ISTOP.eq.1) Then                    ! particle disappeared
            Do j = 1,nmec
              If (lmec(j).eq.meccherenkov) 
     &             cherenkov_on_stack = IQ(JSTAK+1)
            EndDo
          EndIf
        EndIf
      ENDIF
            

C     RICH
C  Put Cherenkov photons on stack
c
      IF (NAMES_TMP(NLEVEL).EQ.'PRGA' .or. NAMES_TMP(NLEVEL).eq.'KRGA') 
     &     Call GSKPHO(0) 

C call subroutine gustep for sub packages

      IF(NAMES_TMP(NLEVEL).EQ.'PRGA') 
     &     CALL PRCH_GUSTEP(cherenkov_on_stack)

      IF(NAMES_TMP(NLEVEL).EQ.'KRGA') 
     &     CALL KRCH_GUSTEP(cherenkov_on_stack)

C     IF ( NAMES_TMP(NLEVEL).EQ.'VMSC' .or.
C    &     NAMES_TMP(NLEVEL).EQ.'VMPB' .or.
C    &     NAMES_TMP(NLEVEL).EQ.'VMFE' .or.
C    &     NAMES_TMP(NLEVEL).EQ.'VMV1' .or.
C    &     NAMES_TMP(NLEVEL).EQ.'VMV2' .or.
C    &     NAMES_TMP(NLEVEL).EQ.'VMV3' .or.
C    &     NAMES_TMP(NLEVEL).EQ.'VVIO' .or.
C    &     NAMES_TMP(NLEVEL).EQ.'VVWN' .or.
C    &     NAMES_TMP(NLEVEL).EQ.'VDV1' .or.
C    &     NAMES_TMP(NLEVEL).EQ.'VDV2' )       
C    &     CALL VVS_GUSTEP
C     IF ( NAMES_TMP(NLEVEL).EQ.'BISC' .or.
C    &     NAMES_TMP(NLEVEL).EQ.'BIPB' .or.
C    &     NAMES_TMP(NLEVEL).EQ.'BIDV' )
C    &     CALL BIVS_GUSTEP

      IF ( NAMES_TMP(NLEVEL).EQ.'PLXX' .or. 
     &     NAMES_TMP(NLEVEL).EQ.'PLX1' .or.
     &     NAMES_TMP(NLEVEL).EQ.'PLUU' .or. 
     &     NAMES_TMP(NLEVEL).EQ.'PLU1' .or.
     &     NAMES_TMP(NLEVEL).EQ.'PLVV' .or. 
     &     NAMES_TMP(NLEVEL).EQ.'PLV1'      ) 
     &  call ums_gustep

      IF ( NAMES_TMP(NLEVEL).EQ.'KTXX' .or. 
     &     NAMES_TMP(NLEVEL).EQ.'KTX1' .or.
     &     NAMES_TMP(NLEVEL).EQ.'KTUU' .or. 
     &     NAMES_TMP(NLEVEL).EQ.'KTU1' .or.
     &     NAMES_TMP(NLEVEL).EQ.'KTVV' .or. 
     &     NAMES_TMP(NLEVEL).EQ.'KTV1'      ) 
     &  call keat_gustep

      IF(NAMES_TMP(NLEVEL).EQ.'STGA') 
     &   call straw_gustep

      IF(NAMES_TMP(NLEVEL).EQ.'FVSC') 
     &   call fvs_gustep

      IF(NAMES_TMP(NLEVEL).EQ.'MVSC') 
     &   call mvs_gustep
c
c     Stop a track if enter any VETO if ISWIT(4) = 1 (DV)
c
      IF(ISWIT(4).EQ.1 .AND. INWVOL.EQ.1 .AND. (
     &   NAMES_TMP(NLEVEL).EQ.'BIPB' .OR.
     &   NAMES_TMP(NLEVEL).EQ.'BISC' .OR.
     &   NAMES_TMP(NLEVEL).EQ.'VMPB' .OR.
     &   NAMES_TMP(NLEVEL).EQ.'VMSC' .OR.
     &   NAMES_TMP(NLEVEL).EQ.'FVSC' .OR.
     &   NAMES_TMP(NLEVEL).EQ.'FVPB' .OR.
     &   NAMES_TMP(NLEVEL).EQ.'MVSC' .OR.
     &   NAMES_TMP(NLEVEL).EQ.'MVSC' )) THEN
         ISTOP = 2
      ENDIF
C
C     Call user step routine
C
      CALL UUSTEP

      IF (ISWIT(1).GE.1) Then
        If (ISWIT(2)/100.ge.1 .or. ITRTYP.ne.7) Then !avoid plotting Ch-Photons
          CALL GSXYZ
          IF (ISWIT(1).GE.2) CALL GDCXYZ
          if (ISWIT(3).eq.1) then
            CALL GPCXYZ                           ! C.M.-5Apr-2000
            if(inwvol.eq.1 .or. inwvol.eq.2) print *,'inwvol=',inwvol
          endif
        EndIf
      EndIf

      CALL UCTOH('NULL',KCASE,4,4) ! Don't write after this line. DV.

      RETURN
      
      END
      
      SUBROUTINE DVDEBUG(n)
      IMPLICIT NONE
#include "gtnum.inc"
#include "gcnum.inc"
#include "gttmed.inc"
#include "gctmed.inc"
#include "gtkine.inc"
#include "gckine.inc"
#include "gcking.inc"
#include "gtflag.inc"
#include "gcflag.inc"
#include "gttrak.inc"
#include "gctrak.inc"
#include "gtvolu.inc"
#include "gcvolu.inc"
#include "gtcuts.inc"
#include "gccuts.inc"
#include "gconsp.inc"
#include "gtscan.inc"
#include "gcscan.inc"
#include "gtlink.inc"
#include "gtbank.inc"
#include "gcbank.inc"
#include "ckm_para.inc"

      INTEGER n, i
      CHARACTER*4 CHCASE
      REAL        ub(10)
      INTEGER     nb, lpart
      REAL        lvert(3), lpvert(4)

      CALL UHTOC(KCASE,4,CHCASE,4)
      
      PRINT *,'CASE ',N,' ISTOP,NGKINE ',ISTOP,NGKINE
      PRINT *,'IEVENT', ievent
      CALL GFKINE(1, lvert, lpvert, lpart, i, ub, nb)
      PRINT 100, lpart, lvert, (lpvert(i),i=1,3)
 100  FORMAT(' KINE -1 ',I2,2(X,F7.3),X,F9.3,3(X,F7.3))
      PRINT *,'RNDM ', NRNDM(1), NRNDM(2)

      PRINT *,'PARTICLES:', ITRA, IPART
      DO i = 1, NGKINE
         print *,GKIN(5,i),GKIN(4,i)
      ENDDO

      PRINT *,'STACK: ',IQ(JSTAK+1)
      DO i = 1, IQ(JSTAK+1)
         PRINT *,IQ(JSTAK+(i-1)*12+5),Q(JSTAK+(i-1)*12+15)
      ENDDO

      PRINT *,'MECHANISMS: ',CHCASE
      DO i = 1, NMEC
         IF(LMEC(i).LT.100) THEN
            WRITE(*,'(A4)') NAMEC(LMEC(i))
         ELSE
            WRITE(*,'(A4)') NAMEC1(LMEC(i)-100)
         ENDIF
      ENDDO

c     PAUSE

      RETURN
      END
C
C  $Log: gustep.F,v $
C  Revision 1.21  2001/03/17 02:47:00  vavilov
C  Bug fixed which prevented an initial particle to be forced to decay.
C
C  Revision 1.20  2001/03/16 02:56:24  vavilov
C  Substantial changes aimed to introduce a new way of keeping track on an
C  initial particle. Initail particle first marked by putting the variable
C  UPWGHT (/GCTRAK/) to some value (MARK) diffrent from default, then keeping this
C  value with the particle through its stops (it turned out that GEANT can stop
C  particle not only for the reason of decay or hadronic interaction but also
C  when it emmits cherenkov photons).
C  Now particle is forced to decay and its decay procuducts are saved to JKINE
C  structure only if it is the initial particle (UPWGHT=MARK)!
C
C  Revision 1.19  2001/03/02 18:37:57  jurgen
C  make names linux compatible
C
C  Revision 1.18  2001/02/28 22:51:01  vavilov
C  Add code that stops any particle enetering any of 4 VETOs: BIVS, VVS, FVS, MVS
C  if ISWIT(4) = 1 (this increases a simulation speed).
C
C  Revision 1.17  2001/02/12 18:23:07  jurgen
C  Fix overwriting of ipart in gfkine calls
C
C  Revision 1.16  2001/02/09 21:32:33  godang
C  fixed KEAT layer from KTYY and KTY1 to KTVV and KTV1
C
C  Revision 1.15  2001/01/29 21:30:21  caroline
C        call bivs_gustep the stepping routines in BIVS
C
C  Revision 1.14  2001/01/17 00:34:43  caroline
C       Updatedd for VVS Endcaps and its windows
C
C  Revision 1.13  2001/01/11 15:41:50  ramberg
C  Fixed UMS layers from PLYY and PLY1 to PLVV and PLV1
C
C  Revision 1.12  2001/01/04 02:23:29  jurgen
C  avoid plotting of cherenkov photons (only for debugging)
C
C  Revision 1.11  2000/12/27 18:46:19  jurgen
C  changed vvs names (Caroline)
C
C  Revision 1.10  2000/12/19 18:00:08  ramberg
C  Added call to KEAT gustep routine
C
C  Revision 1.9  2000/12/04 06:44:49  ramberg
C  - Took out hooks for Linux system.  Now all systems compare names to
C    the equivalenced variables
C  - Added call to ums_gustep
C
C  Revision 1.8  2000/10/06 03:45:32  ramberg
C
C  Substantial revision to the geometry, with new straw geometry from Hogan
C  and new FVS and MVS geometry from Erik.  All three have hit output.
C
C  Revision 1.7  2000/09/12 17:57:49  jurgen
C  Add changes by Caroline concerning vvs configs.
C
C  Revision 1.6  2000/08/17 01:00:03  ckm
C  add linux compatibility.
C  Fix decay point cut.
C
C  Revision 1.5  2000/08/02 17:40:13  ckm
C  put decays products of original particles also into gckine
C  put cherenkov photons on stack
C  add handling of decay volume
C  call package gustep routines
C
C  Revision 1.4  1999/10/28 00:15:53  jurgen
C  Substantial changes by Dima Vavilov
C
c Revision 1.3  1998/05/20  19:14:39  jurgen
c Declare local variable
c
c Revision 1.2  1998/05/20  19:02:12  jurgen
c also include gt...
c
c Revision 1.1  1998/05/13  22:56:43  jurgen
c First time to cvs
c
C
