GEANT321/gtrak/gtelec.F

   1 #include "geant321/pilot.h"
   2 *CMZ :  3.21/04 22/02/95  16.08.31  by  S.Giani
   3 *-- Author :
   4       SUBROUTINE GTELEC
   5 C.
   6 C.    ******************************************************************
   7 C.    *                                                                *
   8 C.    *   Electron type track. Computes step size and propagates       *
   9 C.    *    particle through step.                                      *
  10 C.    *                                                                *
  11 C.    *   ==> Called by : GTRACK                                       *
  12 C.    *       Authors    R.Brun, F.Bruyant, M.Maire L.Urban ********   *
  13 C.    *                                                                *
  14 C.    ******************************************************************
  15 C.
  16 #include "geant321/gcbank.inc"
  17 #include "geant321/gccuts.inc"
  18 #include "geant321/gcjloc.inc"
  19 #include "geant321/gckine.inc"
  20 #include "geant321/gcmate.inc"
  21 #include "geant321/gcmulo.inc"
  22 #include "geant321/gconsp.inc"
  23 #include "geant321/gcphys.inc"
  24 #include "geant321/gcstak.inc"
  25 #include "geant321/gctmed.inc"
  26 #include "geant321/gctrak.inc"
  27 #include "geant321/gcunit.inc"
  28 #include "geant321/gcking.inc"
  29 #if defined(CERNLIB_USRJMP)
  30 #include "geant321/gcjump.inc"
  31 #endif
  32
  33 #if !defined(CERNLIB_OLD)
  34 #include "geant321/gcvolu.inc"
  35 #include "geant321/gcvdma.inc"
  36 #endif
  37 #if !defined(CERNLIB_SINGLE)
  38       PARAMETER (EPSMAC=1.E-6)
  39       DOUBLE PRECISION DEMEAN,STOPRG,STOPMX,STOPC
  40       DOUBLE PRECISION ONE,XCOEF1,XCOEF2,XCOEF3,ZERO
  41 #endif
  42 #if defined(CERNLIB_SINGLE)
  43       PARAMETER (EPSMAC=1.E-11)
  44 #endif
  45       PARAMETER (ONE=1,ZERO=0.)
  46       REAL VNEXT(6)
  47       SAVE IKCUT,STOPC
  48 *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  49       DIMENSION RNDM(3)
  50       PARAMETER ( TLIM = 0.0002)
  51       IABAN = NINT(DPHYS1)
  52 *<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  53 C.
  54 C.    ------------------------------------------------------------------
  55 *
  56 * *** Particle below energy threshold ? short circuit
  57 *
  58       IF (GEKIN.LE.CUTELE) GO TO 100
  59 *
  60 * *** Update local pointers if medium or particle code has changed
  61       IF (IUPD.EQ.0) THEN
  62          IUPD  = 1
  63          JMULOF= LQ(JTM-1)
  64          IF (CHARGE.LT.0.) THEN
  65             JBREM = LQ(JMA-9)
  66             JLOSS = LQ(JMA-1)
  67             JDRAY = LQ(JMA-11)
  68             JRANG = LQ(JMA-15)
  69             JCOEF = LQ(JMA-17)
  70          ELSE
  71             JBREM = LQ(JMA-9)  +NEK1
  72             JLOSS = LQ(JMA-1)  +NEK1
  73             JDRAY = LQ(JMA-11) +NEK1
  74             JRANG = LQ(JMA-15) +NEK1
  75             JCOEF = LQ(JMA-17) +3*NEK1
  76             JANNI = LQ(JMA-7)
  77          ENDIF
  78          IF(IMCKOV.EQ.1) THEN
  79             JTCKOV = LQ(JTM-3)
  80             JABSCO = LQ(JTCKOV-1)
  81             JEFFIC = LQ(JTCKOV-2)
  82             JINDEX = LQ(JTCKOV-3)
  83             JCURIN = LQ(JTCKOV-4)
  84             NPCKOV = Q(JTCKOV+1)
  85          ENDIF
  86          OMCMOL = Q(JPROB+21)
  87          CHCMOL = Q(JPROB+25)
  88          IKCUT  = Q(JMULOF+NEK1+1)
  89          STOPC  = Q(JMULOF+NEK1+2)
  90          IF(ISTRA.GT.0) THEN
  91             JTSTRA = LQ(JMA-19)
  92             JTSTCO = LQ(JTSTRA-1)
  93             JTSTEN = LQ(JTSTRA-2)
  94 #if defined(CERNLIB_ASHO)
  95             IF(ISTRA.EQ.2) THEN
  96                JTASHO = LQ(JMA-20)
  97             ENDIF
  98 #endif
  99          ENDIF
 100       ENDIF
 101 *
 102 * *** Compute current step size
 103 *
 104       STEP   = STEMAX
 105       IPROC  = 103
 106       GEKRT1 = 1. -GEKRAT
 107 *
 108 *  **   Step limitation due to bremsstrahlung ?
 109 *
 110       IF (IBREM.GT.0) THEN
 111          STEPBR = GEKRT1*Q(JBREM+IEKBIN) +GEKRAT*Q(JBREM+IEKBIN+1)
 112          SBREM  = STEPBR*ZINTBR
 113          IF (SBREM.LT.STEP) THEN
 114             STEP  = SBREM
 115             IPROC = 9
 116          ENDIF
 117       ENDIF
 118 *
 119 *  **   Step limitation due to delta-ray production ?
 120 *
 121       IF (IDRAY.GT.0) THEN
 122          STEPDR = GEKRT1*Q(JDRAY+IEKBIN) +GEKRAT*Q(JDRAY+IEKBIN+1)
 123          SDRAY  = STEPDR*ZINTDR
 124          IF (SDRAY.LT.STEP) THEN
 125             STEP  = SDRAY
 126             IPROC = 10
 127          ENDIF
 128       ENDIF
 129 *
 130 *  **   Step limitation due to annihilation ?
 131 *
 132       IF (CHARGE.GT.0.) THEN
 133          IF (IANNI.GT.0) THEN
 134             STEPAN = GEKRT1*Q(JANNI+IEKBIN) +GEKRAT*Q(JANNI+IEKBIN+1)
 135             SANNI  = STEPAN*ZINTAN
 136             IF (SANNI.LT.STEP) THEN
 137                STEP  = SANNI
 138                IPROC = 11
 139             ENDIF
 140          ENDIF
 141       ENDIF
 142 *
 143       IF (STEP.LE.0.) THEN
 144          STEP = 0.
 145          GO TO 90
 146       ENDIF
 147 *
 148 *  **   Step limitation due to energy-loss,multiple scattering
 149 *             or magnetic field ?
 150 *
 151       IF (JMULOF.NE.0) THEN
 152          SMULOF  = GEKRT1*Q(JMULOF+IEKBIN) +GEKRAT*Q(JMULOF+IEKBIN+1)
 153          IF (SMULOF.LT.STEP) THEN
 154             STEP  = SMULOF
 155             IPROC = 0
 156          ENDIF
 157       ENDIF
 158 *
 159 *  **   Step limitation due to geometry ?
 160 *
 161       IF (STEP.GE.0.95*SAFETY) THEN
 162          CALL GTNEXT
 163          IF (IGNEXT.NE.0) THEN
 164             STEP   = SNEXT + PREC
 165             IPROC = 0
 166          ENDIF
 167 *
 168 *        Update SAFETY in stack companions, if any
 169          IF (IQ(JSTAK+3).NE.0) THEN
 170             DO 10 IST = IQ(JSTAK+3),IQ(JSTAK+1)
 171                JST    = JSTAK + 3 + (IST-1)*NWSTAK
 172                Q(JST+11) = SAFETY
 173    10       CONTINUE
 174             IQ(JSTAK+3) = 0
 175          ENDIF
 176       ELSE
 177          IQ(JSTAK+3) = 0
 178       ENDIF
 179 *
 180       IF (IFIELD.EQ.0.OR.STEP.LE.PREC) THEN
 181 *
 182 * *** Linear transport when no field or very short step
 183 *
 184          IF (IGNEXT.NE.0) THEN
 185 *
 186 * *** Particle is supposed to cross the boundary during step
 187 *
 188             DO 20 I = 1,3
 189                VECTMP  = VECT(I) +STEP*VECT(I+3)
 190                IF(VECTMP.EQ.VECT(I)) THEN
 191 *
 192 * *** Correct for machine precision
 193 *
 194                   IF(VECT(I+3).NE.0.) THEN
 195                      VECTMP =
 196      +               VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))*EPSMAC
 197                      IF(NMEC.GT.0) THEN
 198                         IF(LMEC(NMEC).EQ.104) NMEC=NMEC-1
 199                      ENDIF
 200                      NMEC=NMEC+1
 201                      LMEC(NMEC)=104
 202 #if defined(CERNLIB_DEBUG)
 203                      WRITE(CHMAIL, 10000)
 204                      CALL GMAIL(0,0)
 205                      WRITE(CHMAIL, 10100) GEKIN, NUMED, STEP, SNEXT
 206                      CALL GMAIL(0,0)
 207 10000 FORMAT(' Boundary correction in GTELEC: ',
 208      +       '    GEKIN      NUMED       STEP      SNEXT')
 209 10100 FORMAT(31X,E10.3,1X,I10,1X,E10.3,1X,E10.3,1X)
 210 #endif
 211                   ENDIF
 212                ENDIF
 213                VECT(I) = VECTMP
 214    20       CONTINUE
 215             INWVOL = 2
 216             NMEC = NMEC +1
 217             LMEC(NMEC) = 1
 218          ELSE
 219             DO 30 I = 1,3
 220                VECT(I)  = VECT(I) +STEP*VECT(I+3)
 221    30       CONTINUE
 222          ENDIF
 223       ELSE
 224 *
 225 * *** otherwise, swim particle in magnetic field
 226 *
 227 #if !defined(CERNLIB_OLD)
 228          if(mycoun.gt.1.and.nfmany.gt.0.and.step.ge.safety)then
 229            nlevel=manyle(nfmany)
 230            do 99 i=1,nlevel
 231              names(i)=manyna(nfmany,i)
 232              number(i)=manynu(nfmany,i)
 233  99        continue
 234            call glvolu(nlevel,names,number,ier)
 235            if(ier.ne.0)print *,'Fatal error in GLVOLU'
 236            ingoto=0
 237          endif
 238 #endif
 239          NMEC = NMEC +1
 240          LMEC(NMEC) = 4
 241 *
 242 #if !defined(CERNLIB_USRJMP)
 243    40    CALL GUSWIM (CHARGE, STEP, VECT, VOUT)
 244 #endif
 245 #if defined(CERNLIB_USRJMP)
 246    40    CALL JUMPT4(JUSWIM,CHARGE, STEP, VECT, VOUT)
 247 #endif
 248 *
 249 *  ** When near to boundary, take proper action (cut-step,crossing...)
 250 *
 251          IF(STEP.GE.SAFETY)THEN
 252             INEAR = 0
 253             IF (IGNEXT.NE.0) THEN
 254                DO 50 I = 1,3
 255                   VNEXT(I+3) = VECT(I+3)
 256                   VNEXT(I) = VECT(I) +SNEXT*VECT(I+3)
 257    50          CONTINUE
 258                DO 60 I = 1,3
 259                   IF (ABS(VOUT(I)-VNEXT(I)).GT.EPSIL) GO TO 70
 260    60          CONTINUE
 261                INEAR = 1
 262             ENDIF
 263 *
 264    70       CALL GINVOL (VOUT, ISAME)
 265             IF (ISAME.EQ.0) THEN
 266                IF ((INEAR.NE.0).OR.(STEP.LT.EPSIL)) THEN
 267                   INWVOL = 2
 268                   NMEC = NMEC +1
 269                   LMEC(NMEC) = 1
 270                ELSE
 271 *              Cut step
 272                   STEP = 0.5*STEP
 273                   IF (LMEC(NMEC).NE.24) THEN
 274                      NMEC = NMEC +1
 275                      LMEC(NMEC) = 24
 276                   ENDIF
 277                   GO TO 40
 278                ENDIF
 279             ENDIF
 280          ENDIF
 281 *
 282 *
 283          DO 80 I = 1,6
 284             VECT(I) = VOUT(I)
 285    80    CONTINUE
 286 *
 287       ENDIF
 288 *
 289 *
 290 * *** Correct the step due to multiple scattering
 291       IF (IMULL.NE.0) THEN
 292          STMULS = STEP
 293          CORR=0.0001*(STEP/RADL)*(GETOT/(VECT(7)*VECT(7)))**2
 294          IF (CORR.GT.0.25) CORR = 0.25
 295          STEP  = (1.+CORR)*STEP
 296       ENDIF
 297 *
 298       SLENG = SLENG + STEP
 299 *
 300 *  **  Apply synchrotron radiation if required
 301 *
 302       IF(ISYNC*IFIELD.NE.0) THEN
 303          CALL GSYNC
 304          NMEC = NMEC+1
 305          LMEC(NMEC) = 108
 306       ENDIF
 307 *
 308 * *** Generate Cherenkov photons if required
 309 *
 310       IF(IMCKOV.EQ.1) THEN
 311          CALL GGCKOV
 312          IF(NGPHOT.NE.0) THEN
 313             NMEC = NMEC+1
 314             LMEC(NMEC)=105
 315          ENDIF
 316       ENDIF
 317 *
 318 * *** Apply energy loss : find the kinetic energy corresponding
 319 *      to the new stopping range = stopmx - step
 320 *
 321       IF (ILOSL.NE.0) THEN
 322          NMEC = NMEC +1
 323          IF(GEKRAT.LT.0.7) THEN
 324             I1 = MAX(IEKBIN-1,1)
 325          ELSE
 326             I1 = MIN(IEKBIN,NEKBIN-1)
 327          ENDIF
 328          I1 = 3*(I1-1)+1
 329          XCOEF1 = Q(JCOEF+I1)
 330          XCOEF2 = Q(JCOEF+I1+1)
 331          XCOEF3 = Q(JCOEF+I1+2)
 332          IF(XCOEF1.NE.0.) THEN
 333             STOPMX = -XCOEF2+SIGN(ONE,XCOEF1)*SQRT(XCOEF2**2 - (XCOEF3-
 334      +      GEKIN/XCOEF1))
 335          ELSE
 336             STOPMX = - (XCOEF3-GEKIN)/XCOEF2
 337          ENDIF
 338 *
 339          STOPRG = STOPMX - STEP
 340          IF (STOPRG.LT.STOPC) THEN
 341             STEP = MAX(STOPMX - STOPC,ZERO)
 342             GO TO 100
 343          ENDIF
 344 *
 345          LMEC(NMEC) = 3
 346          IF(XCOEF1.NE.0.) THEN
 347             DEMEAN=GEKIN-XCOEF1*(XCOEF3+STOPRG*(2.*XCOEF2+STOPRG))
 348          ELSE
 349             DEMEAN=GEKIN-XCOEF2*STOPRG-XCOEF3
 350          ENDIF
 351          IF(DEMEAN.LE.5.*GEKIN*EPSMAC) THEN
 352             DEMEAN=(GEKRT1*Q(JLOSS+IEKBIN)+GEKRAT*Q(JLOSS+IEKBIN+1))
 353      +     *STEP
 354          ENDIF
 355          IF (ILOSS.EQ.4.OR.IEKBIN.LE.IKCUT+1) THEN
 356             DESTEP = DEMEAN
 357          ELSE
 358             DEMS = DEMEAN
 359             CALL GFLUCT (DEMS,DESTEP)
 360          ENDIF
 361          DESTEP=MAX(DESTEP,0.)
 362          GEKINT = GEKIN -DESTEP
 363          IF (GEKINT.LE.(1.01*CUTELE)) GO TO 100
 364          DESTEL = DESTEP
 365          GEKIN  = GEKINT
 366          GETOT  = GEKIN +AMASS
 367          VECT(7)= SQRT((GETOT+AMASS)*GEKIN)
 368          CALL GEKBIN
 369 *
 370 *        IABAN = 1 does not distinguish between sensitive and
 371 *        non-sensitive volumes, and can stop particles above the CUTE
 372 *
 373          IF(IABAN.EQ.1) THEN
 374 *
 375 *        STOP electron/positron if range < safety AND no brems
 376 *
 377            IF(STOPMX.LE.SAFETY) THEN
 378              IF(SBREM.GT.STOPMX) THEN
 379 *          + condition in case of Cherenkov generation:
 380 *          STOP if E/p >refractive index (i.e. e+/e- is below threshold)
 381                IF(IMCKOV.EQ.1) THEN
 382                  THRIND=GETOT/VECT(7)
 383 *****            IF(THRIND.GT.Q(JINDEX+1)) GOTO 100
 384                  IF(THRIND.GT.Q(JINDEX+1)) GOTO 98
 385                ELSE
 386                  GOTO 98
 387                ENDIF
 388              ENDIF
 389            ENDIF
 390            STOPRG = STOPMX - STEP
 391            IF (STOPRG.LT.STOPC) THEN
 392             STEP = MAX(STOPMX - STOPC,ZERO)
 393             GO TO 98
 394            ENDIF
 395          END IF
 396 *
 397 *        IABAN = 2 distinguishes between sensitive and non-sensitive
 398 *        volumes.
 399 *        In sensitive volumes additional tests are applied before
 400 *        the particle is stopped
 401 *
 402          IF(IABAN.EQ.2) THEN
 403            IF(ISVOL.LE.0) THEN
 404              IF(STOPMX.LE.SAFETY) THEN
 405 *
 406 *        test for brems and annihilation only
 407 *
 408                IF((IBREM.GT.0).AND.(SBREM.LE.STOPMX)) GOTO 97
 409                IF((CHARGE.GT.0.).AND.(IANNI.GT.0).and.
 410      +                                 (SANNI.LE.STOPMX)) GOTO 97
 411                GOTO 98
 412              END IF
 413            ELSE
 414 *
 415 *        sensitive volume ---> more tests !!
 416 *        is energy below TLIM (=200 keV ) ?
 417 *
 418              IF(GEKIN.LE.TLIM) THEN
 419 *
 420 *     range of the particle is the overestimated stopping range here
 421 *
 422                 TOSTOP=STOPMX-STEP*DESTEP/DEMEAN-STOPC
 423 *
 424 *        does the track remain in the actual volume?
 425 *
 426                 IF(TOSTOP.LE.SAFETY) THEN
 427 *
 428 *        is there no delta ray, brems, annihilation ?
 429 *
 430                   IF((IDRAY.GT.0).AND.(SDRAY.LE.TOSTOP)) GOTO 97
 431                   IF((IBREM.GT.0).AND.(SBREM.LE.TOSTOP)) GOTO 97
 432                   IF((CHARGE.GT.0.).AND.(IANNI.GT.0).and.
 433      +                                   (SANNI.LE.TOSTOP)) GOTO 97
 434 *
 435 *        extra condition in case of Cherenkov generation:
 436 *
 437                   IF(IMCKOV.EQ.1) THEN
 438                      THRIND=GETOT/VECT(7)
 439 *
 440 *        continue only if e+/e- below threshold
 441 *
 442                      IF(THRIND.LT.Q(JINDEX+1)) GOTO 97
 443                    ELSE
 444 *
 445 *        do not make transport if this estimated range negative...
 446 *
 447                      IF(TOSTOP.LE.0.) GOTO 98
 448 *
 449 *        estimate final position/direction of the particle
 450 *        from energy loss + multiple scattering
 451 *       ( multiple scattering with path length = range of the particle)
 452 *
 453                      ALFA=0.18*Z
 454                      ALFA1=ALFA+1.
 455                      TGTH2=ALFA*ALFA1/(1.2+ALFA)
 456                      S=TOSTOP*SQRT(1.+TGTH2)/ALFA1
 457                      TGTH=SQRT(TGTH2)
 458                      THET=ATAN(TGTH)
 459                      IF(THET.Lt.0.) THET=PI-THET
 460 *
 461 *        correct direction
 462 *
 463                      CT=COS(THET)
 464                      ST=SQRT(1.-CT*CT)
 465                      CALL GRNDM(RNDM,1)
 466                      PHI=TWOPI*RNDM(1)
 467 *
 468                      D1=ST*COS(PHI)
 469                      D2=ST*SIN(PHI)
 470                      D3=CT
 471                      VMM=SQRT(VECT(4)*VECT(4)+VECT(5)*VECT(5))
 472                      IF(VMM.NE.0.) THEN
 473                        PD1=VECT(4)/VMM
 474                        PD2=VECT(5)/VMM
 475                        V4=PD1*VECT(6)*D1-PD2*D2+VECT(4)*D3
 476                        V5=PD2*VECT(6)*D1+PD1*D2+VECT(5)*D3
 477                        V6=-VMM*D1+VECT(6)*D3
 478                      ELSE
 479                        V4=D1
 480                        V5=D2
 481                        V6=D3
 482                      ENDIF
 483                      VP=1./SQRT(V4*V4+V5*V5+V6*V6)
 484                      VECT(4)=V4*VP
 485                      VECT(5)=V5*VP
 486                      VECT(6)=V6*VP
 487 *
 488 *       transport particle ( assuming FIELD = 0. )
 489 *
 490                      DO 123 I=1,3
 491                      VECT(I)=VECT(I)+S*VECT(I+3)
 492 123                  CONTINUE
 493 *
 494 *       put back into GEKIN the original value in order to have
 495 *       a correct DESTEP
 496 *
 497                   GOTO 98
 498                 END IF
 499               END IF
 500             END IF
 501 *++++++++++++++++++++++++++++++++++++++++++++++++++++++
 502           END IF
 503 97        CONTINUE
 504         ENDIF
 505       ENDIF
 506 *<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 507 *
 508 * *** Apply multiple scattering.
 509 *
 510       IF (IMULL.NE.0) THEN
 511          NMEC = NMEC +1
 512          LMEC(NMEC) = 2
 513          CALL GMULTS
 514       ENDIF
 515 *
 516 * *** Update time of flight
 517 *
 518       TOFG = TOFG +STEP*GETOT/(VECT(7)*CLIGHT)
 519       IF (TOFG.GE.TOFMAX) THEN
 520          ISTOP = 4
 521          NMEC  = NMEC +1
 522          LMEC(NMEC) = 22
 523          GO TO 999
 524       ENDIF
 525 *
 526 * *** Update interaction probabilities
 527 *
 528       IF (IBREM.GT.0)    ZINTBR = ZINTBR -STEP/STEPBR
 529       IF (IDRAY.GT.0)    ZINTDR = ZINTDR -STEP/STEPDR
 530       IF (CHARGE.GT.0.) THEN
 531          IF (IANNI.GT.0) ZINTAN = ZINTAN -STEP/STEPAN
 532       ENDIF
 533 *
 534 * *** Apply the selected process if any
 535 *
 536    90 IF (IPROC.EQ.0) GO TO 999
 537       NMEC = NMEC +1
 538       LMEC(NMEC) = IPROC
 539 *
 540 *  **   Bremsstrahlung ?
 541 *
 542       IF (IPROC.EQ.9) THEN
 543          CALL GBREME
 544 *
 545 *  **   Delta ray ?
 546 *
 547       ELSE IF (IPROC.EQ.10) THEN
 548 *
 549        IF((IPART.EQ.2).OR.((IPART.EQ.3).AND.(GEKIN.GT.2.*DCUTE))) THEN
 550          CALL GDRAY
 551        ELSE
 552          GOTO 98
 553        ENDIF
 554 *
 555 *  **   Positron annihilation ?
 556 *
 557       ELSE IF (IPROC.EQ.11) THEN
 558          CALL GANNI
 559
 560       ENDIF
 561       GO TO 999
 562 *
 563 * *** Special treatment for overstopped tracks
 564 *
 565   98  GEKIN=GEKIN+DESTEP
 566   100 DESTEP = GEKIN
 567       DESTEL = DESTEP
 568       GEKIN  = 0.
 569       GETOT  = AMASS
 570       VECT(7)= 0.
 571       INWVOL = 0
 572       NMEC   = NMEC +1
 573       LMEC(NMEC) = 30
 574       IF ((CHARGE.LT.0.).OR.(IANNI.EQ.0)) THEN
 575          ISTOP = 2
 576       ELSE
 577          NMEC = NMEC +1
 578          LMEC(NMEC) = 11
 579          CALL GANNIR
 580       ENDIF
 581   999 IF(NGPHOT.GT.0) THEN
 582          IF(ITCKOV.EQ.2.AND.ISTOP.EQ.0) THEN
 583 *
 584 *  The electron has produced Cerenkov photons and it is still alive
 585 *  we put it in the stack and we let the photons to be tracked
 586             NGKINE = NGKINE+1
 587             GKIN(1,NGKINE) = VECT(4)*VECT(7)
 588             GKIN(2,NGKINE) = VECT(5)*VECT(7)
 589             GKIN(3,NGKINE) = VECT(6)*VECT(7)
 590             GKIN(4,NGKINE) = GETOT
 591             GKIN(5,NGKINE) = IPART
 592             TOFD(NGKINE) = 0.
 593             ISTOP = 1
 594 c----put position as well
 595             GPOS(1,NGKINE)=VECT(1)
 596             GPOS(2,NGKINE)=VECT(2)
 597             GPOS(3,NGKINE)=VECT(3)
 598          ENDIF
 599       ENDIF
 600       END
 601