GEANT321/erdecks/ertrch.F

   1 *CMZ :          24/11/95  16.28.16  by  S.Ravndal
   2 *CMZ :  3.21/02 29/03/94  15.41.49  by  S.Giani
   3 *-- Author :
   4       SUBROUTINE ERTRCH
   5 C.
   6 C.    ******************************************************************
   7 C.    *                                                                *
   8 C.    *    Average charged track is extrapolated by one step           *
   9 C.    *                                                                *
  10 C.    *    ==>Called by : ERTRGO                                       *
  11 C.    *       Original routine : GTHADR                                *
  12 C.    *       Authors   M.Maire, E.Nagy  *********                     *
  13 C.    *                                                                *
  14 C.    ******************************************************************
  15 C.
  16 #include "geant321/gcbank.inc"
  17 #include "geant321/gconst.inc"
  18 #include "geant321/gccuts.inc"
  19 #include "geant321/gcphys.inc"
  20 #include "geant321/gcjloc.inc"
  21 #include "geant321/gckine.inc"
  22 #include "geant321/gcmate.inc"
  23 #include "geant321/gcmulo.inc"
  24 #include "geant321/gctmed.inc"
  25 #include "geant321/gctrak.inc"
  26 #include "geant321/gcunit.inc"
  27 #include "geant321/ertrio.inc"
  28 #include "geant321/erwork.inc"
  29
  30 #if defined(CERNLIB_SINGLE)
  31       PARAMETER (EPSMAC=1.E-11)
  32 #endif
  33 #if !defined(CERNLIB_SINGLE)&&!defined(CERNLIB_IBM)
  34       PARAMETER (EPSMAC=5.E-6)
  35 #endif
  36 #if !defined(CERNLIB_SINGLE)&&defined(CERNLIB_IBM)
  37       PARAMETER (EPSMAC=5.E-5)
  38 #endif
  39 #if !defined(CERNLIB_SINGLE)
  40       DOUBLE PRECISION GKR,DEMEAN,STOPP1,STOPP2,STOPMX,STOPRG,STOPC
  41       DOUBLE PRECISION EKIPR
  42 #endif
  43       REAL VNEXT(6)
  44       SAVE CFLD,CHARG2,RMASS,CUTPRO,IKCUT,STOPC
  45 C.
  46 C.    ------------------------------------------------------------------
  47 *
  48 *
  49 * *** Update local pointers if medium has changed
  50 *
  51       IF (IUPD.EQ.0) THEN
  52          IUPD  = 1
  53          CHARG2 = CHARGE*CHARGE
  54          IF (IPART.LE.3) THEN
  55             CUTEK  = CUTELE
  56             RMASS  = 1.
  57             JRANG  = LQ(JMA-15)
  58          ELSE IF (IPART.LE.6) THEN
  59             CUTEK  = CUTMUO
  60             RMASS  = 1.
  61             JRANG  = LQ(JMA-16)
  62          ELSE
  63             CUTEK  = CUTHAD
  64             RMASS  = PMASS/AMASS
  65             JRANG  = LQ(JMA-16) + NEK1
  66          ENDIF
  67          CUTPRO = MAX(CUTEK*RMASS,ELOW(1))
  68          IKCUT = GEKA*LOG10(CUTPRO) + GEKB
  69          GKR   = (CUTPRO - ELOW(IKCUT))/(ELOW(IKCUT+1) - ELOW(IKCUT))
  70          STOPC = (1.-GKR)*Q(JRANG+IKCUT) + GKR*Q(JRANG+IKCUT+1)
  71          CFLD  = 0.
  72          IF (FIELDM.NE.0.) CFLD = 3333.*DEGRAD*TMAXFD/ABS(FIELDM*CHARGE)
  73       ENDIF
  74 *
  75 * *** Compute current step size
  76 *
  77       STEP   = BIG
  78       GEKRT1 = 1. - GEKRAT
  79 *
  80 * *** Step limitation due to energy loss (stopping range) ?
  81 *
  82       IF (ILOSS*DEEMAX.GT.0.) THEN
  83          STOPP1 = GEKRT1*Q(JRANG+IEKBIN) + GEKRAT*Q(JRANG+IEKBIN+1)
  84          STOPMX = (STOPP1 - STOPC)/(RMASS*CHARG2)
  85          EKF  = (1. - BACKTR*DEEMAX)*GEKIN*RMASS
  86          IF (EKF.LT.ELOW(1))    EKF = ELOW(1)
  87          IF (EKF.GE.ELOW(NEK1)) EKF = ELOW(NEK1)*0.99
  88          IKF=GEKA*LOG10(EKF)+GEKB
  89          GKR=(EKF-ELOW(IKF))/(ELOW(IKF+1)-ELOW(IKF))
  90          STOPP2 = (1.-GKR)*Q(JRANG+IKF) + GKR*Q(JRANG+IKF+1)
  91          SLOSP  = ABS (STOPP1 - STOPP2)
  92          STEP   = SLOSP/(RMASS*CHARG2)
  93       ENDIF
  94 *
  95 * *** Step limitation due to energy loss in magnetic field ?
  96 *
  97       IF (IFIELD*FIELDM.NE.0.) THEN
  98          SFIELD = CFLD*VECT(7)
  99          IF (SFIELD.LT.STEP) STEP = SFIELD
 100       ENDIF
 101 *
 102 * *** Compute point where to store error matrix
 103 *
 104       LERST  = 0
 105       STEPER = BIG
 106       ASCL1  = BIG
 107       DO 20 IPR = 1,NEPRED
 108          STEPE  = BIG
 109          IF (LELENG) STEPE = ERLENG(IPR) - SLENG
 110          IF (LEPLAN) THEN
 111             SCAL1 = 0.
 112             SCAL2 = 0.
 113             DO 18 I=1,3
 114                SCAL1 = SCAL1 + ERPLO(I,4,IPR)*(ERPLO(I,3,IPR)-VECT(I))
 115                SCAL2 = SCAL2 + ERPLO(I,4,IPR)*VECT(I+3)
 116    18       CONTINUE
 117             STEPE = SCAL1/SCAL2
 118          ENDIF
 119          IF (STEPE.LE.PREC) STEPE = BIG
 120          IF (STEPE.LT.STEPER) THEN
 121             STEPER = STEPE
 122             INLIST = IPR
 123             IF (LEPLAN) ASCL1  = ABS (SCAL1)
 124          ENDIF
 125    20 CONTINUE
 126       IF (STEPER.LE.STEP)  THEN
 127          STEP  = STEPER
 128          LERST = 1
 129       ENDIF
 130 *
 131 * *** Step limitation due to geometry ?
 132 *
 133       LNEXT  = 0
 134       IF (STEP.GE.0.95*SAFETY) THEN
 135          CALL GTNEXT
 136          IF (IGNEXT.NE.0) THEN
 137             STEP   = SNEXT + PREC
 138             LNEXT = 1
 139             IF ((STEPER-SNEXT).GT.(2*PREC)) LERST = 0
 140          ENDIF
 141       ENDIF
 142 *
 143 * *** Linear transport when no field or very short step
 144 *
 145       IF (IFIELD.EQ.0.OR.STEP.LE.2*PREC) THEN
 146         IF (IGNEXT.NE.0) THEN
 147           DO 25 I = 1,3
 148             VECTMP = VECT(I) +STEP*VECT(I+3)
 149             IF(VECTMP.EQ.VECT(I)) THEN
 150 *
 151 * *** Correct for machine precision
 152 *
 153                   IF(VECT(I+3).NE.0.) THEN
 154                      VECTMP =
 155      +               VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))*EPSMAC
 156 #if defined(DEBUG)
 157                      WRITE(CHMAIL, 10000)
 158                      CALL GMAIL(0,0)
 159                      WRITE(CHMAIL, 10100) GEKIN, NUMED, STEP, SNEXT
 160                      CALL GMAIL(0,0)
 161 10000 FORMAT(' Boundary correction in ERTRCH: ',
 162      +       '    GEKIN      NUMED       STEP      SNEXT')
 163 10100 FORMAT(31X,E10.3,1X,I10,1X,E10.3,1X,E10.3,1X)
 164 #endif
 165                   ENDIF
 166             ENDIF
 167             VOUT(I) = VECTMP
 168    25     CONTINUE
 169             INWVOL = 2
 170             NMEC = NMEC +1
 171             LMEC(NMEC) = 1
 172         ELSE
 173             DO 30 I = 1,3
 174                VOUT(I)  = VECT(I) +STEP*VECT(I+3)
 175    30       CONTINUE
 176         ENDIF
 177         DO 35 I = 4,6
 178           VOUT(I)  = VECT(I)
 179    35   CONTINUE
 180         GOTO 74
 181       END IF
 182 *
 183 * *** otherwise, swim particle in magnetic field
 184 *
 185       NMEC = NMEC +1
 186       LMEC(NMEC) = 4
 187 *
 188    50 LERST = 0
 189       LNEXT = 0
 190       CALL GUSWIM (CHTR , STEP, VECT, VOUT)
 191 *
 192 *     When near to boundary, take proper action (cut-step,crossing...)
 193       IF (STEP.GE.SAFETY) THEN
 194          INEAR = 0
 195          IF (IGNEXT.NE.0) THEN
 196            DO 51 I = 1,3
 197                VNEXT(I+3) = VECT(I+3)
 198                VNEXT(I) = VECT(I) +SNEXT*VECT(I+3)
 199    51      CONTINUE
 200            DO 52 I = 1,3
 201              IF (ABS(VOUT(I)-VNEXT(I)).GT.EPSIL) GO TO 55
 202    52      CONTINUE
 203            INEAR = 1
 204          ENDIF
 205 *
 206    55    CALL GINVOL (VOUT,ISAME)
 207          IF (ISAME.EQ.0) THEN
 208            IF ((INEAR.NE.0).OR.(STEP.LT.EPSIL)) THEN
 209              INWVOL = 2
 210              NMEC = NMEC +1
 211              LMEC(NMEC) = 1
 212              LNEXT = 1
 213            ELSE
 214 *              Cut step
 215              STEP = 0.5*STEP
 216              IF (LMEC(NMEC).NE.24) THEN
 217                NMEC = NMEC +1
 218                LMEC(NMEC) = 24
 219              ENDIF
 220              GOTO 50
 221            ENDIF
 222          ENDIF
 223       ENDIF
 224 *
 225 *
 226 *     preset plane reached  ?
 227    74 CONTINUE
 228       IF ((LEPLAN).AND.(STEP.GE.ASCL1)) THEN
 229          SCAL3 = 0.
 230          DO 28 I=1,3
 231             SCAL3=SCAL3+ERPLO(I,4,INLIST)*(ERPLO(I,3,INLIST)-VOUT(I))
 232    28    CONTINUE
 233          ASCL3 = ABS(SCAL3)
 234          SSCL1 = ASCL1/SCAL1
 235          IF (SCAL3*SSCL1.LT. -PREC) THEN
 236 *            Cut step
 237              STEP  = STEP*(ASCL1/(ASCL1+ASCL3))
 238              NMEC  = NMEC +1
 239              LMEC(NMEC) = 24
 240              GOTO 50
 241          ELSE
 242            IF(ASCL3.LE.PREC) LERST = 1
 243          ENDIF
 244       ENDIF
 245 *
 246       DO 75 I=1,6
 247            VECT(I) = VOUT(I)
 248    75 CONTINUE
 249 *
 250       IF (LELENG.AND.(STEP.GE.STEPER)) LERST = 1
 251 *
 252       SLENG = SLENG + STEP
 253 *
 254 * *** Now apply selected mechanisms
 255 *
 256       IF (LNEXT.EQ.1) THEN
 257           INWVOL = 2
 258           NMEC = NMEC + 1
 259           LMEC(NMEC) = 1
 260       ENDIF
 261 *
 262 * *** apply energy loss : find the kinetic energy corresponding
 263 *      to the new stopping range = stopmx -/+ step
 264 *      (take care of the back tracking !)
 265 *
 266       IF (ILOSS*DEEMAX.GT.0) THEN
 267          NMEC = NMEC +1
 268          LMEC(NMEC) = 3
 269          CALL ERLAND (STEP,Z,A,DENS,VECT(7),GETOT,AMASS,DEDX2)
 270          DEDX2  = DEDX2*CHARG2*CHARG2
 271          STOPRG = STOPP1 - BACKTR*STEP*RMASS*CHARG2
 272          IKF = IEKBIN
 273          IF (BACKTR.LE.0.) THEN
 274    95       IF (STOPRG.LT.Q(JRANG+IKF)) THEN
 275                IKF = IKF - 1
 276                IF (IKF.GT.1) GO TO 95
 277             ENDIF
 278          ELSE
 279    96       IF (STOPRG.GE.Q(JRANG+IKF+1)) THEN
 280                IKF = IKF + 1
 281                IF (IKF.LT.NEK1) GO TO 96
 282             ENDIF
 283          ENDIF
 284          GKR = (STOPRG - Q(JRANG+IKF)) / (Q(JRANG+IKF+1) - Q(JRANG+IKF))
 285          EKIPR = (1. -GKR)*ELOW(IKF) + GKR*ELOW(IKF+1)
 286          GEKINT =  EKIPR/RMASS
 287          IF (GEKINT.GT.CUTEK) THEN
 288             DESTEP = ABS (GEKIN - GEKINT)
 289             GEKIN  = GEKINT
 290             GETOT  = GEKIN + AMASS
 291             VECT(7)= SQRT((GETOT+AMASS)*GEKIN)
 292             CALL GEKBIN
 293          ELSE
 294             DESTEP = GEKINT
 295             GEKIN  = 0.
 296             GETOT  = AMASS
 297             VECT(7)= 0.
 298             INWVOL = 0
 299             ISTOP  = 2
 300             NMEC = NMEC + 1
 301             LMEC(NMEC) = 30
 302          ENDIF
 303       ENDIF
 304 *
 305 * *** Propagate error matrix
 306 *
 307       IF (.NOT. LEONLY) CALL ERPROP
 308 *
 309 * *** Store informations
 310 *
 311       IF(LERST.EQ.1) THEN
 312          NMEC = NMEC + 1
 313          LMEC(NMEC) = 27
 314          CALL ERSTOR
 315       ENDIF
 316 *
 317       END