CDECK  ID>, GTHADR.
*CMZ :  3.21/02 03/07/94  17.57.24  by  S.Giani
*-- Author :
      SUBROUTINE GTHADR
C.
C.    ******************************************************************
C.    *                                                                *
C.    *   Charged hadron type track. Computes step size and propagates *
C.    *    particle through step.                                      *
C.    *                                                                *
C.    *   ==>Called by : GTRACK                                        *
C.    *       Authors    R.Brun, F.Bruyant, M.Maire ********           *
C.    *                                                                *
C.    ******************************************************************
C.
      PARAMETER (KWBANK=69000,KWWORK=5200)
      COMMON/GCBANK/NZEBRA,GVERSN,ZVERSN,IXSTOR,IXDIV,IXCONS,FENDQ(16)
     +             ,LMAIN,LR1,WS(KWBANK)
      DIMENSION IQ(2),Q(2),LQ(8000),IWS(2)
      EQUIVALENCE (Q(1),IQ(1),LQ(9)),(LQ(1),LMAIN),(IWS(1),WS(1))
      EQUIVALENCE (JCG,JGSTAT)
      COMMON/GCLINK/JDIGI ,JDRAW ,JHEAD ,JHITS ,JKINE ,JMATE ,JPART
     +      ,JROTM ,JRUNG ,JSET  ,JSTAK ,JGSTAT,JTMED ,JTRACK,JVERTX
     +      ,JVOLUM,JXYZ  ,JGPAR ,JGPAR2,JSKLT
C
      COMMON/GCCUTS/CUTGAM,CUTELE,CUTNEU,CUTHAD,CUTMUO,BCUTE,BCUTM
     +             ,DCUTE ,DCUTM ,PPCUTM,TOFMAX,GCUTS(5)
C
      COMMON/GCJLOC/NJLOC(2),JTM,JMA,JLOSS,JPROB,JMIXT,JPHOT,JANNI
     +                  ,JCOMP,JBREM,JPAIR,JDRAY,JPFIS,JMUNU,JRAYL
     +                 ,JMULOF,JCOEF,JRANG,JPHTA
C
      INTEGER       NJLOC   ,JTM,JMA,JLOSS,JPROB,JMIXT,JPHOT,JANNI
     +                  ,JCOMP,JBREM,JPAIR,JDRAY,JPFIS,JMUNU,JRAYL
     +                 ,JMULOF,JCOEF,JRANG,JPHTA
C
      COMMON/GCJLCK/NJLCK(2),JTCKOV,JABSCO,JEFFIC,JINDEX,JCURIN
     +                      ,JPOLAR,JTSTRA,JTSTCO,JTSTEN,JTASHO
C
      EQUIVALENCE (JLASTV,JTSTEN)
C
      INTEGER       NJLCK,JTCKOV,JABSCO,JEFFIC,JINDEX,JCURIN
     +                   ,JPOLAR,JLASTV,JTSTRA,JTSTCO,JTSTEN
     +                   ,JTASHO
C
      COMMON/GCKINE/IKINE,PKINE(10),ITRA,ISTAK,IVERT,IPART,ITRTYP
     +      ,NAPART(5),AMASS,CHARGE,TLIFE,VERT(3),PVERT(4),IPAOLD
C
      INTEGER MXGKIN
      PARAMETER (MXGKIN=100)
      COMMON/GCKING/KCASE,NGKINE,GKIN(5,MXGKIN),
     +                           TOFD(MXGKIN),IFLGK(MXGKIN)
      INTEGER       KCASE,NGKINE ,IFLGK,MXPHOT,NGPHOT
      REAL          GKIN,TOFD,XPHOT
C
      PARAMETER (MXPHOT=800)
      COMMON/GCKIN2/NGPHOT,XPHOT(11,MXPHOT)
C
      COMMON/GCKIN3/GPOS(3,MXGKIN)
      REAL          GPOS
C
      COMMON/GCMATE/NMAT,NAMATE(5),A,Z,DENS,RADL,ABSL
C
      INTEGER NMAT,NAMATE
      REAL A,Z,DENS,RADL,ABSL
C
      COMMON/GCMULO/SINMUL(101),COSMUL(101),SQRMUL(101),OMCMOL,CHCMOL
     +  ,EKMIN,EKMAX,NEKBIN,NEK1,EKINV,GEKA,GEKB,EKBIN(200),ELOW(200)
C
      REAL SINMUL,COSMUL,SQRMUL,OMCMOL,CHCMOL,EKMIN,EKMAX,ELOW,EKINV
      REAL GEKA,GEKB,EKBIN
      INTEGER NEKBIN,NEK1
C
      DOUBLE PRECISION PI,TWOPI,PIBY2,DEGRAD,RADDEG,CLIGHT,BIG,EMASS
      DOUBLE PRECISION EMMU,PMASS,AVO
*
      PARAMETER (PI=3.14159265358979324D0)
      PARAMETER (TWOPI=6.28318530717958648D0)
      PARAMETER (PIBY2=1.57079632679489662D0)
      PARAMETER (DEGRAD=0.0174532925199432958D0)
      PARAMETER (RADDEG=57.2957795130823209D0)
      PARAMETER (CLIGHT=29979245800.D0)
      PARAMETER (BIG=10000000000.D0)
      PARAMETER (EMASS=0.0005109990615D0)
      PARAMETER (EMMU=0.105658387D0)
      PARAMETER (PMASS=0.9382723128D0)
      PARAMETER (AVO=0.60221367D0)
*
      COMMON/GCPHYS/IPAIR,SPAIR,SLPAIR,ZINTPA,STEPPA
     +             ,ICOMP,SCOMP,SLCOMP,ZINTCO,STEPCO
     +             ,IPHOT,SPHOT,SLPHOT,ZINTPH,STEPPH
     +             ,IPFIS,SPFIS,SLPFIS,ZINTPF,STEPPF
     +             ,IDRAY,SDRAY,SLDRAY,ZINTDR,STEPDR
     +             ,IANNI,SANNI,SLANNI,ZINTAN,STEPAN
     +             ,IBREM,SBREM,SLBREM,ZINTBR,STEPBR
     +             ,IHADR,SHADR,SLHADR,ZINTHA,STEPHA
     +             ,IMUNU,SMUNU,SLMUNU,ZINTMU,STEPMU
     +             ,IDCAY,SDCAY,SLIFE ,SUMLIF,DPHYS1
     +             ,ILOSS,SLOSS,SOLOSS,STLOSS,DPHYS2
     +             ,IMULS,SMULS,SOMULS,STMULS,DPHYS3
     +             ,IRAYL,SRAYL,SLRAYL,ZINTRA,STEPRA
     +             ,IPHTA,SPHTA,SLPHTA,ZINTAP,STEPAP
      COMMON/GCPHLT/ILABS,SLABS,SLLABS,ZINTLA,STEPLA
     +             ,ISYNC
     +             ,ISTRA
*
      PARAMETER (NWSTAK=12,NWINT=11,NWREAL=12,NWTRAC=NWINT+NWREAL+6)
      COMMON /GCSTAK/ NJTMAX, NJTMIN, NTSTKP, NTSTKS, NDBOOK, NDPUSH,
     +                NJFREE, NJGARB, NJINVO, LINSAV(15), LMXSAV(15)
      EQUIVALENCE (ISTORD,NJTMIN)
C
      COMMON/GCTMED/NUMED,NATMED(5),ISVOL,IFIELD,FIELDM,TMAXFD,STEMAX
     +      ,DEEMAX,EPSIL,STMIN,CFIELD,PREC,IUPD,ISTPAR,NUMOLD
      COMMON/GCTLIT/THRIND,PMIN,DP,DNDL,JMIN,ITCKOV,IMCKOV,NPCKOV
C
      PARAMETER (MAXMEC=30)
      COMMON/GCTRAK/VECT(7),GETOT,GEKIN,VOUT(7),NMEC,LMEC(MAXMEC)
     + ,NAMEC(MAXMEC),NSTEP ,MAXNST,DESTEP,DESTEL,SAFETY,SLENG
     + ,STEP  ,SNEXT ,SFIELD,TOFG  ,GEKRAT,UPWGHT,IGNEXT,INWVOL
     + ,ISTOP ,IGAUTO,IEKBIN, ILOSL, IMULL,INGOTO,NLDOWN,NLEVIN
     + ,NLVSAV,ISTORY
      PARAMETER (MAXME1=30)
      COMMON/GCTPOL/POLAR(3), NAMEC1(MAXME1)
C
      COMMON/GCUNIT/LIN,LOUT,NUNITS,LUNITS(5)
      INTEGER LIN,LOUT,NUNITS,LUNITS
      COMMON/GCMAIL/CHMAIL
      CHARACTER*132 CHMAIL
C
      PARAMETER    (MAXJMP=30)
      COMMON/GCJUMP/JUDCAY, JUDIGI, JUDTIM, JUFLD , JUHADR, JUIGET,
     +              JUINME, JUINTI, JUKINE, JUNEAR, JUOUT , JUPHAD,
     +              JUSKIP, JUSTEP, JUSWIM, JUTRAK, JUTREV, JUVIEW,
     +              JUPARA
      DIMENSION     JMPADR(MAXJMP)
      EQUIVALENCE  (JMPADR(1), JUDCAY)
*
 
      COMMON/GCVOLU/NLEVEL,NAMES(15),NUMBER(15),
     +LVOLUM(15),LINDEX(15),INFROM,NLEVMX,NLDEV(15),LINMX(15),
     +GTRAN(3,15),GRMAT(10,15),GONLY(15),GLX(3)
C
      COMMON/GCVDMA/NVMANY,MANYLE(20),MANYNA(20,15),
     +MANYNU(20,15),NFMANY,MYCOUN,IMYSE,RAYTRA,VECCOS(3)
C
C
      PARAMETER (EPSMAC=1.E-6)
      DOUBLE PRECISION GKR,DEMEAN,STOPP,STOPMX,STOPRG,STOPC,EKIPR
      DOUBLE PRECISION ONE,XCOEF1,XCOEF2,XCOEF3,YCOEF1,YCOEF2,YCOEF3
      PARAMETER (THRESH=0.7,ONE=1)
      REAL VNEXT(6)
      SAVE CFLD,CHARG2,RMASS,CUTPRO,IKCUT,STOPC
C.
C.    ------------------------------------------------------------------
*
* *** Particle below energy threshold ? short circuit
*
      IF (GEKIN.LE.CUTHAD) GO TO 100
*
* *** Update local pointers if medium has changed
*
      IF (IUPD.EQ.0) THEN
         IUPD  = 1
         JLOSS = LQ(JMA-3)
         JRANG = LQ(JMA-16) + NEK1
         JCOEF = LQ(JMA-18) + 3*NEK1
         CHARG2 = CHARGE*CHARGE
         RMASS  = PMASS/AMASS
         OMCMOL = Q(JPROB+21)*CHARG2
         CHCMOL = Q(JPROB+25)*ABS(CHARGE)
         CUTPRO = MAX(CUTHAD*RMASS,ELOW(1))
         IKCUT = GEKA*LOG10(CUTPRO) + GEKB
         GKR   = (CUTPRO - ELOW(IKCUT))/(ELOW(IKCUT+1) - ELOW(IKCUT))
         STOPC = (1.-GKR)*Q(JRANG+IKCUT) + GKR*Q(JRANG+IKCUT+1)
         IF (FIELDM.NE.0.) THEN
            CFLD = 3333.*DEGRAD*TMAXFD/ABS(FIELDM*CHARGE)
         ELSE
            CFLD = BIG
         ENDIF
         IF(IMCKOV.EQ.1) THEN
            JTCKOV = LQ(JTM-3)
            JABSCO = LQ(JTCKOV-1)
            JEFFIC = LQ(JTCKOV-2)
            JINDEX = LQ(JTCKOV-3)
            JCURIN = LQ(JTCKOV-4)
            NPCKOV = Q(JTCKOV+1)
         ENDIF
         IF(ISTRA.GT.0) THEN
            JTSTRA = LQ(JMA-19)
            JTSTCO = LQ(JTSTRA-1)
            JTSTEN = LQ(JTSTRA-2)
         ENDIF
      ENDIF
*
* *** Compute current step size
*
      STEP   = STEMAX
      IPROC  = 103
      GEKRT1 = 1. -GEKRAT
*
*  **   Step limitation due to hadron interaction ?
*
      IF (IHADR.GT.0) THEN
         CALL JUMPT0(JUPHAD)
         IF (SHADR.LT.STEP) THEN
            IF (SHADR.LE.0.) SHADR = PREC
            STEP  = SHADR
            IPROC = 12
         ENDIF
      ENDIF
*
*  **   Step limitation due to decay ?
*
      IF (IDCAY.GT.0) THEN
         SDCAY = SUMLIF*VECT(7)/AMASS
         IF (SDCAY.LT.STEP) THEN
            STEP  = SDCAY
            IPROC = 5
         ENDIF
      ENDIF
*
*  ** Step limitation due to delta-ray production ?
*       (Cannot be tabulated easily because dependent on AMASS)
*
      IF (IDRAY.GT.0) THEN
         STEPDR = BIG
         IF (GEKIN.GT.DCUTM) THEN
            GAMASS = GETOT +AMASS
            TMAX   = EMASS*GEKIN*GAMASS/(0.5*AMASS*AMASS+EMASS*GETOT)
            IF (TMAX.GT.DCUTM) THEN
               BET2 = GEKIN*GAMASS/(GETOT*GETOT)
               Y    = DCUTM/TMAX
               SIG  = (1.-Y+BET2*Y*LOG(Y))/DCUTM
*              extra term for spin 1/2
               IF (AMASS.GT.0.9) SIG=SIG+0.5*(TMAX-DCUTM)/(GETOT*GETOT)
               SIG = SIG*Q(JPROB+17)*CHARG2*EMASS/BET2
*
               IF (SIG.GT.0.) THEN
                  STEPDR = 1./SIG
                  SDRAY  = STEPDR*ZINTDR
                  IF (SDRAY.LE.STEP) THEN
                     STEP  = SDRAY
                     IPROC = 10
                  ENDIF
               ENDIF
            ENDIF
         ENDIF
      ENDIF
*
      IF (STEP.LE.0.) THEN
         STEP  = 0.
         GO TO 110
      ENDIF
*
*  **   Step limitation due to energy loss (stopping range) ?
*
      IF (ILOSL.GT.0) THEN
         IF(GEKRAT.LT.THRESH) THEN
            I1 = MAX(IEKBIN-1,1)
         ELSE
            I1 = MIN(IEKBIN,NEKBIN-1)
         ENDIF
         I1 = 3*(I1-1)+1
         XCOEF1 = Q(JCOEF+I1)
         XCOEF2 = Q(JCOEF+I1+1)
         XCOEF3 = Q(JCOEF+I1+2)
         IF(XCOEF1.NE.0) THEN
            STOPP = -XCOEF2+SIGN(ONE,XCOEF1)* SQRT(XCOEF2
     +      **2 -(XCOEF3-GEKIN*RMASS/XCOEF1))
         ELSE
            STOPP = - (XCOEF3-GEKIN*RMASS)/XCOEF2
         ENDIF
         STOPMX = (STOPP - STOPC)/(RMASS*CHARG2)
         IF (STOPMX.LT.MIN(STEP,STMIN)) THEN
            STEP = STOPMX
            IPROC = 0
            IF(STEP.LE.0.)THEN
               GO TO 100
            ENDIF
            GO TO 10
         ENDIF
         EKF = (1. - DEEMAX)*GEKIN*RMASS
         IF (EKF.LT.ELOW(1)) THEN
            EKF = ELOW(1)
         ELSEIF (EKF.GE.ELOW(NEK1)) THEN
            EKF = ELOW(NEK1)*0.99
         ENDIF
         IKF=GEKA*LOG10(EKF)+GEKB
         GKR=(EKF-ELOW(IKF))/(ELOW(IKF+1)-ELOW(IKF))
         IF(GKR.LT.THRESH) THEN
            IK1 = MAX(IKF-1,1)
         ELSE
            IK1 = MIN(IKF,NEKBIN-1)
         ENDIF
         IK1 = 3*(IK1-1)+1
         YCOEF1=Q(JCOEF+IK1)
         YCOEF2=Q(JCOEF+IK1+1)
         YCOEF3=Q(JCOEF+IK1+2)
         IF(YCOEF1.NE.0.) THEN
            SLOSP = -YCOEF2+SIGN(ONE,YCOEF1)*SQRT(YCOEF2**2- (YCOEF3-
     +      EKF/YCOEF1))
         ELSE
            SLOSP = - (YCOEF3-EKF)/YCOEF2
         ENDIF
         SLOSP = STOPP - SLOSP
         SLOSS = MAX(STMIN, SLOSP/(RMASS*CHARG2) )
         IF (SLOSS.LT.STEP) THEN
            STEP = SLOSS
            IPROC = 0
         ENDIF
      ENDIF
*
*  **   Step limitation due to bending in magnetic field ?
*
      IF (IFIELD.NE.0) THEN
         SFIELD = CFLD*VECT(7)
         SFIELD=MAX(SFIELD, STMIN)
         IF (SFIELD.LT.STEP) THEN
            STEP  = SFIELD
            IPROC = 0
         ENDIF
      ENDIF
*
*  **   Step limitation due to multiple scattering ?
*
      IF (IMULL.GT.0) THEN
         SMULS=MIN(2232.*RADL*((VECT(7)**2)/(GETOT*CHARGE))**2,10.*RADL)
         SMULS  = MAX(STMIN, SMULS )
         IF (SMULS.LT.STEP) THEN
            STEP  = SMULS
            IPROC = 0
         ENDIF
      ENDIF
*
   10 CONTINUE
*
*  **   Step limitation due to Cerenkov production ?
*
      IF (IMCKOV.GT.0) THEN
         CALL GNCKOV
         STCKOV = MXPHOT/MAX(3.*DNDL,1E-10)
         SMULS  = MAX(STMIN, STCKOV)
         IF (SMULS.LT.STEP) THEN
            STEP  = STCKOV
            IPROC = 0
         ENDIF
      ENDIF
*
*  **   Step limitation due to geometry ?
*
      IF (STEP.GE.0.95*SAFETY) THEN
         CALL GTNEXT
         IF (IGNEXT.NE.0) THEN
            STEP   = SNEXT + PREC
            IPROC = 0
         ENDIF
*
*        Update SAFETY in stack companions, if any
         IF (IQ(JSTAK+3).NE.0) THEN
            DO 20 IST = IQ(JSTAK+3),IQ(JSTAK+1)
               JST    = JSTAK + 3 + (IST-1)*NWSTAK
               Q(JST+11) = SAFETY
   20       CONTINUE
            IQ(JSTAK+3) = 0
         ENDIF
      ELSE
         IQ(JSTAK+3) = 0
      ENDIF
*
* *** Linear transport when no field or very short step
*
      IF (IFIELD.EQ.0.OR.STEP.LE.PREC) THEN
*
         IF (IGNEXT.NE.0) THEN
            DO 30 I = 1,3
               VECTMP  = VECT(I) +STEP*VECT(I+3)
               IF(VECTMP.EQ.VECT(I)) THEN
*
* *** Correct for machine precision
*
                  IF(VECT(I+3).NE.0.) THEN
                     VECTMP =
     +               VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))*EPSMAC
                     IF(NMEC.GT.0) THEN
                        IF(LMEC(NMEC).EQ.104) NMEC=NMEC-1
                     ENDIF
                     NMEC=NMEC+1
                     LMEC(NMEC)=104
                  ENDIF
               ENDIF
               VECT(I) = VECTMP
   30       CONTINUE
            INWVOL = 2
            NMEC = NMEC +1
            LMEC(NMEC) = 1
         ELSE
            DO 40 I = 1,3
               VECT(I)  = VECT(I) +STEP*VECT(I+3)
   40       CONTINUE
         ENDIF
      ELSE
*
* ***   otherwise, swim particle in magnetic field
*
         if(mycoun.gt.1.and.nfmany.gt.0.and.step.ge.safety)then
           nlevel=manyle(nfmany)
           do 99 i=1,nlevel
             names(i)=manyna(nfmany,i)
             number(i)=manynu(nfmany,i)
 99        continue
           call glvolu(nlevel,names,number,ier)
           if(ier.ne.0)print *,'Fatal error in GLVOLU'
           ingoto=0
         endif
         NMEC = NMEC +1
         LMEC(NMEC) = 4
*
   50    CALL JUMPT4(JUSWIM, CHARGE, STEP, VECT, VOUT)
c
c       Update the total energy, if the muon has been accelrated
C       Stolen from P. Lebrun's DPgeant  23 Sep 98
c
        if( abs(vout(7) - vect(7)) .gt. 1.e-6) then
           getot = sqrt(vout(7)**2 + amass**2)
           gekin = getot - amass
           call gekbin
        end if
*
*  ** When near to boundary, take proper action (cut-step,crossing...)
*
         IF(STEP.GE.SAFETY)THEN
            INEAR = 0
            IF (IGNEXT.NE.0) THEN
               DO 60 I = 1,3
                  VNEXT(I+3) = VECT(I+3)
                  VNEXT(I) = VECT(I) +SNEXT*VECT(I+3)
   60          CONTINUE
               DO 70 I = 1,3
                  IF (ABS(VOUT(I)-VNEXT(I)).GT.EPSIL) GO TO 80
   70          CONTINUE
               INEAR = 1
            ENDIF
*
   80       CALL GINVOL (VOUT, ISAME)
            IF (ISAME.EQ.0)THEN
               IF ((INEAR.NE.0).OR.(STEP.LT.EPSIL)) THEN
                  INWVOL = 2
                  NMEC = NMEC +1
                  LMEC(NMEC) = 1
               ELSE
*              Cut step
                  STEP = 0.5*STEP
                  IF (LMEC(NMEC).NE.24) THEN
                     NMEC = NMEC +1
                     LMEC(NMEC) = 24
                  ENDIF
                  GO TO 50
               ENDIF
            ENDIF
         ENDIF
*
         DO 90 I = 1,6
            VECT(I) = VOUT(I)
   90    CONTINUE
*
      ENDIF
*
* *** Correct the step due to multiple scattering
      IF (IMULL.NE.0) THEN
         STMULS = STEP
         CORR=0.0001*CHARG2*(STEP/RADL)*(GETOT/(VECT(7)*VECT(7)))**2
         IF (CORR.GT.0.25) CORR = 0.25
         STEP  = (1.+CORR)*STEP
      ENDIF
*
      SLENG = SLENG + STEP
*
* *** Generate Cherenkov photons if required
*
      IF(IMCKOV.EQ.1) THEN
         CALL GGCKOV
         IF(NGPHOT.NE.0) THEN
            NMEC=NMEC+1
            LMEC(NMEC)=105
         ENDIF
      ENDIF
*
* *** apply energy loss : find the kinetic energy corresponding
*      to the new stopping range = stopmx - step
*
      IF (ILOSL.NE.0) THEN
         NMEC = NMEC +1
         LMEC(NMEC) = 3
         STOPRG = STOPP - STEP*RMASS*CHARG2
         IF (STOPRG.LE.STOPC) THEN
            STEP = STOPMX
            GO TO 100
         ENDIF
         IF(XCOEF1.NE.0.) THEN
            EKIPR = XCOEF1*(XCOEF3+STOPRG*(2.*XCOEF2+STOPRG))
         ELSE
            EKIPR = XCOEF2*STOPRG+XCOEF3
         ENDIF
         DEMEAN=GEKIN - EKIPR/RMASS
         IF(DEMEAN.LE.5.*GEKIN*EPSMAC) THEN
            DEMEAN=(GEKRT1*Q(JLOSS+IEKBIN)+GEKRAT*Q(JLOSS+IEKBIN+1))
     +      *STEP*CHARG2
         ENDIF
         IF (ILOSS.EQ.4.OR.IEKBIN.LE.IKCUT+1) THEN
            DESTEP = DEMEAN
         ELSE
            DEMS   = DEMEAN/CHARG2
            CALL GFLUCT(DEMS,DESTEP)
            DESTEP = DESTEP*CHARG2
         ENDIF
         DESTEP=MAX(DESTEP,0.)
         GEKINT = GEKIN -DESTEP
         IF (GEKINT.LE.(1.01*CUTHAD)) GO TO 100
         DESTEL = DESTEP
         GEKIN  = GEKINT
         GETOT  = GEKIN +AMASS
         VECT(7)= SQRT((GETOT+AMASS)*GEKIN)
         CALL GEKBIN
      ENDIF
*
* *** Apply multiple scattering.
*
      IF (IMULL.NE.0) THEN
         NMEC = NMEC +1
         LMEC(NMEC) = 2
         CALL GMULTS
      ENDIF
*
* *** Update time of flight
*
      SUMLIF = SUMLIF -STEP*AMASS/VECT(7)
      TOFG   = TOFG +STEP*GETOT/(VECT(7)*CLIGHT)
      IF (TOFG.GE.TOFMAX) THEN
         ISTOP = 4
         NMEC  = NMEC +1
         LMEC(NMEC) = 22
         GO TO 999
      ENDIF
*
* *** Update interaction probabilities
*
      IF (IHADR.GT.0) ZINTHA = ZINTHA*(1.-STEP/SHADR)
      IF (IDRAY.GT.0) ZINTDR = ZINTDR -STEP/STEPDR
*
      GO TO 110
*
*  **   Special treatment for overstopped tracks
*
  100 DESTEP = GEKIN
      DESTEL = DESTEP
      GEKIN  = 0.
      GAMMA  = GETOT/AMASS
      GETOT  = AMASS
      VECT(7)= 0.
      INWVOL = 0
      ISTOP  = 2
      NMEC = NMEC + 1
      LMEC(NMEC) = 30
      IF (IHADR.EQ.0) THEN
         IF (IDCAY.NE.0) THEN
            TOFG = TOFG +0.5*(1+GAMMA)*SUMLIF/CLIGHT
            SUMLIF = 0.
            IF (TOFG.GE.TOFMAX) THEN
               ISTOP = 4
               NMEC = NMEC + 1
               LMEC(NMEC) = 22
               NGKINE = 0
            ELSE
               NMEC = NMEC + 1
               LMEC(NMEC) = 5
               ISTOP =1
               CALL GDECAY
            ENDIF
         ENDIF
         GO TO 999
      ENDIF
      IPROC = 12
*
* *** apply slected process if any
*
  110 IF (IPROC.EQ.0) GO TO 999
      NMEC = NMEC +1
      LMEC(NMEC) = IPROC
*
*  **   Hadron interaction ?
*
      IF (IPROC.EQ.12) THEN
         CALL JUMPT0(JUHADR)
*   *   Check time cut-off for decays at rest
         IF (LMEC(NMEC).EQ.5) THEN
            TOFG   = TOFG +SUMLIF/CLIGHT
            SUMLIF = 0.
            IF (TOFG.GE.TOFMAX) THEN
               NGKINE = 0
               ISTOP  = 4
               LMEC(NMEC) = 22
            ENDIF
         ENDIF
*
*  **   Decay ?
*
      ELSE IF (IPROC.EQ.5) THEN
         ISTOP = 1
         CALL GDECAY
*
*  **   Delta-ray ?
*
      ELSE IF (IPROC.EQ.10) THEN
         CALL GDRAY
      ENDIF
  999 IF(NGPHOT.GT.0) THEN
         IF(ITCKOV.EQ.2.AND.ISTOP.EQ.0) THEN
*
*  The hadron has produced Cerenkov photons and it is still alive
*  we put it in the stack and we let the photons to be tracked
            NGKINE = NGKINE+1
            GKIN(1,NGKINE) = VECT(4)*VECT(7)
            GKIN(2,NGKINE) = VECT(5)*VECT(7)
            GKIN(3,NGKINE) = VECT(6)*VECT(7)
            GKIN(4,NGKINE) = GETOT
            GKIN(5,NGKINE) = IPART
            TOFD(NGKINE) = 0.
            ISTOP = 1
c----put position as well
            GPOS(1,NGKINE)=VECT(1)
            GPOS(2,NGKINE)=VECT(2)
            GPOS(3,NGKINE)=VECT(3)
         ENDIF
      ENDIF
      END