GEANT321/gtrak/gtneut.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:44  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.24  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GTNEUT
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *   Neutral hadron type track. Computes step size and propagates *
  17 C.    *    particle through step.                                      *
  18 C.    *                                                                *
  19 C.    *   ==>Called by : GTRACK                                        *
  20 C.    *      Authors     R.Brun, F.Bruyant  *********                  *
  21 C.    *                                                                *
  22 C.    ******************************************************************
  23 C.
  24 #include "geant321/gcbank.inc"
  25 #include "geant321/gccuts.inc"
  26 #include "geant321/gckine.inc"
  27 #include "geant321/gcking.inc"
  28 #include "geant321/gconsp.inc"
  29 #include "geant321/gcphys.inc"
  30 #include "geant321/gcstak.inc"
  31 #include "geant321/gctmed.inc"
  32 #include "geant321/gctrak.inc"
  33 #if defined(CERNLIB_USRJMP)
  34 #include "geant321/gcjump.inc"
  35 #endif
  36 #if defined(CERNLIB_DEBUG)
  37 #include "geant321/gcunit.inc"
  38 #endif
  39 #if !defined(CERNLIB_SINGLE)
  40       PARAMETER (EPSMAC=1.E-6)
  41 #endif
  42 #if defined(CERNLIB_SINGLE)
  43       PARAMETER (EPSMAC=1.E-11)
  44 #endif
  45 C.
  46 C.    ------------------------------------------------------------------
  47 *
  48 * *** Particle below energy threshold ? Special treatment
  49 *
  50       IF (GEKIN.LE.CUTNEU) THEN
  51          GEKIN  = 0.
  52          GETOT  = AMASS
  53          VECT(7)= 0.
  54          ISTOP  = 2
  55          NMEC = NMEC + 1
  56          LMEC(NMEC) = 30
  57          IF (IHADR.EQ.0.AND.AMASS.GT.0.) THEN
  58             IF (IDCAY.NE.0) THEN
  59                GAMMA  = GETOT/AMASS
  60                TOFG = TOFG +GAMMA*SUMLIF/CLIGHT
  61                SUMLIF = 0.
  62                IF (TOFG.GE.TOFMAX) THEN
  63                   ISTOP = 4
  64                   NMEC = NMEC + 1
  65                   LMEC(NMEC) = 22
  66                   NGKINE = 0
  67                ELSE
  68                   NMEC = NMEC + 1
  69                   LMEC(NMEC) = 5
  70                   ISTOP =1
  71                   CALL GDECAY
  72                ENDIF
  73             ENDIF
  74             GO TO 999
  75          ENDIF
  76          IPROC = 12
  77          GO TO 40
  78       ENDIF
  79 *
  80 * *** Compute step size
  81 *
  82       IPROC  = 103
  83       STEP   = STEMAX
  84 *
  85 *  **   Step limitation due to hadron interaction ?
  86 *
  87       IF (IHADR.GT.0) THEN
  88 #if !defined(CERNLIB_USRJMP)
  89          CALL GUPHAD
  90 #endif
  91 #if defined(CERNLIB_USRJMP)
  92          CALL JUMPT0(JUPHAD)
  93 #endif
  94          IF (SHADR.LT.STEP) THEN
  95             IF (SHADR.LE.0.) SHADR = PREC
  96             STEP  = SHADR
  97             IPROC = 12
  98          ENDIF
  99       ENDIF
 100 *
 101 *  **   Step limitation due to decay ?
 102 *
 103       IF ((IDCAY.NE.0).AND.(AMASS.GT.0.)) THEN
 104          SDCAY = SUMLIF*VECT(7)/AMASS
 105          IF (SDCAY.LT.STEP) THEN
 106             STEP  = SDCAY
 107             IPROC = 5
 108          ENDIF
 109       ENDIF
 110 *
 111       IF (STEP.LT.0.) STEP = 0.
 112 *
 113 *  **   Step limitation due to geometry ?
 114 *
 115       IF (STEP.GE.SAFETY) THEN
 116          CALL GTNEXT
 117          IF (IGNEXT.NE.0) THEN
 118             STEP = SNEXT + PREC
 119             IPROC = 0
 120             INWVOL = 2
 121             NMEC = NMEC + 1
 122             LMEC(NMEC)  = 1
 123          ENDIF
 124 *
 125 *        Update SAFETY in stack companions, if any
 126          IF (IQ(JSTAK+3).NE.0) THEN
 127             DO 10 IST = IQ(JSTAK+3),IQ(JSTAK+1)
 128                JST = JSTAK +3 +(IST-1)*NWSTAK
 129                Q(JST+11) = SAFETY
 130    10       CONTINUE
 131             IQ(JSTAK+3) = 0
 132          ENDIF
 133 *
 134       ELSE
 135          IQ(JSTAK+3) = 0
 136       ENDIF
 137 *
 138 * *** Linear transport
 139 *
 140       IF (INWVOL.EQ.2) THEN
 141          DO 20 I = 1,3
 142             VECTMP  = VECT(I) +STEP*VECT(I+3)
 143             IF(VECTMP.EQ.VECT(I)) THEN
 144 *
 145 * *** Correct for machine precision
 146 *
 147                IF(VECT(I+3).NE.0.) THEN
 148                   VECTMP = VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))*
 149      +            EPSMAC
 150                   IF(NMEC.GT.0) THEN
 151                      IF(LMEC(NMEC).EQ.104) NMEC=NMEC-1
 152                   ENDIF
 153                   NMEC=NMEC+1
 154                   LMEC(NMEC)=104
 155 #if defined(CERNLIB_DEBUG)
 156                   WRITE(CHMAIL, 10000)
 157                   CALL GMAIL(0,0)
 158                   WRITE(CHMAIL, 10100) GEKIN, NUMED, STEP, SNEXT
 159                   CALL GMAIL(0,0)
 160 10000 FORMAT(' Boundary correction in GTNEUT: ',
 161      +       '    GEKIN      NUMED       STEP      SNEXT')
 162 10100 FORMAT(31X,E10.3,1X,I10,1X,E10.3,1X,E10.3,1X)
 163 #endif
 164                ENDIF
 165             ENDIF
 166             VECT(I) = VECTMP
 167    20    CONTINUE
 168       ELSE
 169          DO 30 I = 1,3
 170             VECT(I)  = VECT(I) +STEP*VECT(I+3)
 171    30    CONTINUE
 172       ENDIF
 173 *
 174       SLENG = SLENG +STEP
 175 *
 176 * *** Update time of flight
 177 *
 178       SUMLIF = SUMLIF -STEP*AMASS/VECT(7)
 179       TOFG = TOFG +STEP*GETOT/(VECT(7)*CLIGHT)
 180       IF (TOFG.GE.TOFMAX) THEN
 181          ISTOP = 4
 182          NMEC  = NMEC +1
 183          LMEC(NMEC) = 22
 184          GO TO 999
 185       ENDIF
 186 *
 187 * *** Update interaction probabilities
 188 *
 189       IF (IHADR.GT.0) ZINTHA = ZINTHA*(1.-STEP/SHADR)
 190 *
 191 * *** apply the selected process if any
 192 *
 193    40 IF (IPROC.EQ.0) GO TO 999
 194       NMEC = NMEC +1
 195       LMEC(NMEC) = IPROC
 196 *
 197 *  **   Hadron interaction ?
 198 *
 199       IF (IPROC.EQ.12) THEN
 200 #if !defined(CERNLIB_USRJMP)
 201          CALL GUHADR
 202 #endif
 203 #if defined(CERNLIB_USRJMP)
 204          CALL JUMPT0(JUHADR)
 205 #endif
 206 *          Check time cut-off for decays at rest
 207          IF (LMEC(NMEC).EQ.5) THEN
 208             TOFG   = TOFG +SUMLIF/CLIGHT
 209             SUMLIF = 0.
 210             IF (TOFG.GE.TOFMAX) THEN
 211                ISTOP = 4
 212                LMEC(NMEC) = 22
 213                NGKINE = 0
 214             ENDIF
 215          ENDIF
 216 *
 217 *  **   Decay ?
 218 *
 219       ELSE IF (IPROC.EQ.5) THEN
 220          ISTOP = 1
 221          CALL GDECAY
 222       ENDIF
 223 *                                                             END GTNEUT
 224   999 END