GEANT321/gtrak/gthion.F

   1 #include "geant321/pilot.h"
   2 *CMZ :  3.21/02 03/07/94  17.58.49  by  S.Giani
   3 *-- Author :
   4       SUBROUTINE GTHION
   5 C.
   6 C.    ******************************************************************
   7 C.    *                                                                *
   8 C.    *   Heavy  ion  type track. Computes step size and propagates    *
   9 C.    *    particle through step.                                      *
  10 C.    *                                                                *
  11 C.    *   The ionisation energy loss is calculated here (mean +        *
  12 C.    *       fluctuations)                                            *
  13 C.    *   The fluctuations are the same for ILOSS=1,2,3 and            *
  14 C.    *       there is no fluctuation for ILOSS=4.                     *
  15 C.    *                                                                *
  16 C.    *   ==>Called by : GTRACK                                        *
  17 C.    *       Authors    R.Brun, F.Bruyant, M.Maire, L.Urban ***       *
  18 C.    *                                                                *
  19 C.    ******************************************************************
  20 C.
  21 #include "geant321/gcbank.inc"
  22 #include "geant321/gccuts.inc"
  23 #include "geant321/gcjloc.inc"
  24 #include "geant321/gckine.inc"
  25 #include "geant321/gcking.inc"
  26 #include "geant321/gcmate.inc"
  27 #include "geant321/gcmulo.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 #include "geant321/gcunit.inc"
  34 #if defined(CERNLIB_USRJMP)
  35 #include "geant321/gcjump.inc"
  36 #endif
  37
  38 #if !defined(CERNLIB_OLD)
  39 #include "geant321/gcvolu.inc"
  40 #include "geant321/gcvdma.inc"
  41 #endif
  42 #if !defined(CERNLIB_SINGLE)
  43       PARAMETER (EPSMAC=1.E-6)
  44       DOUBLE PRECISION GKR,DEMEAN,STOPP,STOPMX,STOPRG,STOPC,EKIPR
  45       DOUBLE PRECISION ONE,XCOEF1,XCOEF2,XCOEF3,YCOEF1,YCOEF2,YCOEF3
  46 #endif
  47 #if defined(CERNLIB_SINGLE)
  48       PARAMETER (EPSMAC=1.E-11)
  49 #endif
  50       PARAMETER (THRESH=0.7,ONE=1)
  51       PARAMETER (TWOTHR=2*ONE/3,AMU=0.9314943)
  52       PARAMETER (DME=7.84572E-8,CNORM=2.5)
  53       REAL VNEXT(6)
  54       DIMENSION RNDM(2)
  55       SAVE RMASS,CUTPRO,IKCUT,STOPC,FACFLU,CHAR23
  56 C.
  57 C.    ------------------------------------------------------------------
  58 *
  59 * *** Particle below energy threshold ? short circuit
  60 *
  61       IF (GEKIN.LE.CUTHAD) GO TO 100
  62 *
  63 * *** Update local pointers if medium has changed
  64 *
  65       IF (IUPD.EQ.0) THEN
  66          IUPD  = 1
  67          JLOSS = LQ(JMA-3)
  68          JRANG = LQ(JMA-16) + NEK1
  69          JCOEF = LQ(JMA-18) + 3*NEK1
  70          RMASS  = PMASS/AMASS
  71          CUTPRO = MAX(CUTHAD*RMASS,ELOW(1))
  72          IKCUT = GEKA*LOG10(CUTPRO) + GEKB
  73          GKR   = (CUTPRO - ELOW(IKCUT))/(ELOW(IKCUT+1) - ELOW(IKCUT))
  74          STOPC = (1.-GKR)*Q(JRANG+IKCUT) + GKR*Q(JRANG+IKCUT+1)
  75          FACFLU = DME*(Z*DENS/A)
  76          CHAR23 = ONE/CHARGE**TWOTHR
  77          IF(IMCKOV.EQ.1) THEN
  78             JTCKOV = LQ(JTM-3)
  79             JABSCO = LQ(JTCKOV-1)
  80             JEFFIC = LQ(JTCKOV-2)
  81             JINDEX = LQ(JTCKOV-3)
  82             JCURIN = LQ(JTCKOV-4)
  83             NPCKOV = Q(JTCKOV+1)
  84          ENDIF
  85       ENDIF
  86 *
  87 * *** Compute energy dependent parameters
  88 *
  89       GAMASS=GETOT+AMASS
  90       BET2=GEKIN*GAMASS/(GETOT*GETOT)
  91       BET=SQRT(BET2)
  92       W1=1.034-0.1777*EXP(-0.08114*CHARGE)
  93       W2=BET*CHAR23
  94       W3=121.4139*W2+0.0378*SIN(190.7165*W2)
  95       CHARG1=CHARGE*(1.-W1*EXP(-W3))
  96 *
  97 *              the effective charge  CHARG1
  98 *            can be negative only for very low energy and
  99 *     for CHARGE > 20 ( very low energy : T/A < 20 keV/nucleon)
 100 *              in this case short circuit
 101 *
 102       IF(CHARG1.LT.0.) GOTO 100
 103       CHARG2=CHARG1**2
 104 *
 105       OMCMOL=Q(JPROB+21)*CHARG2
 106       CHCMOL=Q(JPROB+25)*ABS(CHARG1)
 107       IF(FIELDM.NE.0.) THEN
 108          CFLD=3333.*DEGRAD*TMAXFD/ABS(FIELDM*CHARG1)
 109       ELSE
 110          CFLD=BIG
 111       ENDIF
 112 *
 113 * *** Compute current step size
 114 *
 115       STEP   = STEMAX
 116       IPROC  = 103
 117       GEKRT1 = 1. -GEKRAT
 118 *
 119 *  **   Step limitation due to hadron interaction ?
 120 *
 121       IF (IHADR.GT.0) THEN
 122 #if !defined(CERNLIB_USRJMP)
 123          CALL GUPHAD
 124 #endif
 125 #if defined(CERNLIB_USRJMP)
 126          CALL JUMPT0(JUPHAD)
 127 #endif
 128          IF (SHADR.LT.STEP) THEN
 129             IF (SHADR.LE.0.) SHADR = PREC
 130             STEP  = SHADR
 131             IPROC = 12
 132          ENDIF
 133       ENDIF
 134 *
 135 *  ** Step limitation due to delta-ray production ?
 136 *       (Cannot be tabulated easily because dependent on AMASS)
 137 *
 138       IF (IDRAY.GT.0) THEN
 139          STEPDR = BIG
 140          IF (GEKIN.GT.DCUTM) THEN
 141             TMAX   = EMASS*GEKIN*GAMASS/(0.5*AMASS*AMASS+EMASS*GETOT)
 142             IF (TMAX.GT.DCUTM) THEN
 143                Y    = DCUTM/TMAX
 144                SIG  = (1.-Y+BET2*Y*LOG(Y))/DCUTM
 145 *              extra term for spin 1/2
 146                IF (AMASS.GT.0.9) SIG=SIG+0.5*(TMAX-DCUTM)/(GETOT*GETOT)
 147                SIG = SIG*Q(JPROB+17)*CHARG2*EMASS/BET2
 148 *
 149                IF (SIG.GT.0.) THEN
 150                   STEPDR = 1./SIG
 151                   SDRAY  = STEPDR*ZINTDR
 152                   IF (SDRAY.LE.STEP) THEN
 153                      STEP  = SDRAY
 154                      IPROC = 10
 155                   ENDIF
 156                ENDIF
 157             ENDIF
 158          ENDIF
 159       ENDIF
 160 *
 161       IF (STEP.LE.0.) THEN
 162          STEP  = 0.
 163          GO TO 110
 164       ENDIF
 165 *
 166 *  **   Step limitation due to energy loss (stopping range) ?
 167 *
 168       IF (ILOSL.GT.0) THEN
 169          IF(GEKRAT.LT.THRESH) THEN
 170             I1 = MAX(IEKBIN-1,1)
 171          ELSE
 172             I1 = MIN(IEKBIN,NEKBIN-1)
 173          ENDIF
 174          I1 = 3*(I1-1)+1
 175          XCOEF1 = Q(JCOEF+I1)
 176          XCOEF2 = Q(JCOEF+I1+1)
 177          XCOEF3 = Q(JCOEF+I1+2)
 178          IF(XCOEF1.NE.0) THEN
 179             STOPP = -XCOEF2+SIGN(ONE,XCOEF1)* SQRT(XCOEF2
 180      +      **2 -(XCOEF3-GEKIN*RMASS/XCOEF1))
 181          ELSE
 182             STOPP = - (XCOEF3-GEKIN*RMASS)/XCOEF2
 183          ENDIF
 184          STOPMX = (STOPP - STOPC)/(RMASS*CHARG2)
 185          IF (STOPMX.LT.MIN(STEP,STMIN)) THEN
 186             STEP = STOPMX
 187             IPROC = 0
 188             IF(STEP.LE.0.)THEN
 189                GO TO 100
 190             ENDIF
 191             GO TO 10
 192          ENDIF
 193          EKF = (1. - DEEMAX)*GEKIN*RMASS
 194          IF (EKF.LT.ELOW(1)) THEN
 195             EKF = ELOW(1)
 196          ELSEIF (EKF.GE.ELOW(NEK1)) THEN
 197             EKF = ELOW(NEK1)*0.99
 198          ENDIF
 199          IKF=GEKA*LOG10(EKF)+GEKB
 200          GKR=(EKF-ELOW(IKF))/(ELOW(IKF+1)-ELOW(IKF))
 201          IF(GKR.LT.THRESH) THEN
 202             IK1 = MAX(IKF-1,1)
 203          ELSE
 204             IK1 = MIN(IKF,NEKBIN-1)
 205          ENDIF
 206          IK1 = 3*(IK1-1)+1
 207          YCOEF1=Q(JCOEF+IK1)
 208          YCOEF2=Q(JCOEF+IK1+1)
 209          YCOEF3=Q(JCOEF+IK1+2)
 210          IF(YCOEF1.NE.0.) THEN
 211             SLOSP = -YCOEF2+SIGN(ONE,YCOEF1)*SQRT(YCOEF2**2- (YCOEF3-
 212      +      EKF/YCOEF1))
 213          ELSE
 214             SLOSP = - (YCOEF3-EKF)/YCOEF2
 215          ENDIF
 216          SLOSP = STOPP - SLOSP
 217          SLOSS = MAX(STMIN, SLOSP/(RMASS*CHARG2) )
 218          IF (SLOSS.LT.STEP) THEN
 219             STEP = SLOSS
 220             IPROC = 0
 221          ENDIF
 222       ENDIF
 223 *
 224 *  **   Step limitation due to energy loss in magnetic field ?
 225 *
 226       IF (IFIELD.NE.0) THEN
 227          SFIELD = CFLD*VECT(7)
 228          SFIELD=MAX(SFIELD, STMIN)
 229          IF (SFIELD.LT.STEP) THEN
 230             STEP  = SFIELD
 231             IPROC = 0
 232          ENDIF
 233       ENDIF
 234 *
 235 *  **   Step limitation due to multiple scattering ?
 236 *
 237       IF (IMULL.GT.0) THEN
 238          SMULS=MIN(2232.*RADL*((VECT(7)**2)/(GETOT*CHARG1))**2,10.*RADL)
 239          SMULS  = MAX(STMIN, SMULS )
 240          IF (SMULS.LT.STEP) THEN
 241             STEP  = SMULS
 242             IPROC = 0
 243          ENDIF
 244       ENDIF
 245 *
 246    10 CONTINUE
 247 *
 248 *  **   Step limitation due to Cerenkov production ?
 249 *
 250       IF (IMCKOV.GT.0) THEN
 251          CALL GNCKOV
 252          STCKOV = MXPHOT/MAX(3.*DNDL,1E-10)
 253          SMULS  = MAX(STMIN, STCKOV)
 254          IF (SMULS.LT.STEP) THEN
 255             STEP  = STCKOV
 256             IPROC = 0
 257          ENDIF
 258       ENDIF
 259 *
 260 *  **   Step limitation due to geometry ?
 261 *
 262       IF (STEP.GE.0.95*SAFETY) THEN
 263          CALL GTNEXT
 264          IF (IGNEXT.NE.0) THEN
 265             STEP   = SNEXT + PREC
 266             IPROC = 0
 267          ENDIF
 268 *
 269 *        Update SAFETY in stack companions, if any
 270          IF (IQ(JSTAK+3).NE.0) THEN
 271             DO 20 IST = IQ(JSTAK+3),IQ(JSTAK+1)
 272                JST    = JSTAK + 3 + (IST-1)*NWSTAK
 273                Q(JST+11) = SAFETY
 274    20       CONTINUE
 275             IQ(JSTAK+3) = 0
 276          ENDIF
 277       ELSE
 278          IQ(JSTAK+3) = 0
 279       ENDIF
 280 *
 281 * *** Linear transport when no field or very short step
 282 *
 283       IF (IFIELD.EQ.0.OR.STEP.LE.PREC) THEN
 284 *
 285          IF (IGNEXT.NE.0) THEN
 286             DO 30 I = 1,3
 287                VECTMP  = VECT(I) +STEP*VECT(I+3)
 288                IF(VECTMP.EQ.VECT(I)) THEN
 289 *
 290 * *** Correct for machine precision
 291 *
 292                   IF(VECT(I+3).NE.0.) THEN
 293                      VECTMP =
 294      +               VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))*EPSMAC
 295                      IF(NMEC.GT.0) THEN
 296                         IF(LMEC(NMEC).EQ.104) NMEC=NMEC-1
 297                      ENDIF
 298                      NMEC=NMEC+1
 299                      LMEC(NMEC)=104
 300 #if defined(CERNLIB_DEBUG)
 301                      WRITE(CHMAIL, 10000)
 302                      CALL GMAIL(0,0)
 303                      WRITE(CHMAIL, 10100) GEKIN, NUMED, STEP, SNEXT
 304                      CALL GMAIL(0,0)
 305 10000 FORMAT(' Boundary correction in GTHION: ',
 306      +       '    GEKIN      NUMED       STEP      SNEXT')
 307 10100 FORMAT(31X,E10.3,1X,I10,1X,E10.3,1X,E10.3,1X)
 308 #endif
 309                   ENDIF
 310                ENDIF
 311                VECT(I) = VECTMP
 312    30       CONTINUE
 313             INWVOL = 2
 314             NMEC = NMEC +1
 315             LMEC(NMEC) = 1
 316          ELSE
 317             DO 40 I = 1,3
 318                VECT(I)  = VECT(I) +STEP*VECT(I+3)
 319    40       CONTINUE
 320          ENDIF
 321       ELSE
 322 *
 323 * ***   otherwise, swim particle in magnetic field
 324 *
 325 #if !defined(CERNLIB_OLD)
 326          if(mycoun.gt.1.and.nfmany.gt.0.and.step.ge.safety)then
 327            nlevel=manyle(nfmany)
 328            do 99 i=1,nlevel
 329              names(i)=manyna(nfmany,i)
 330              number(i)=manynu(nfmany,i)
 331  99        continue
 332            call glvolu(nlevel,names,number,ier)
 333            if(ier.ne.0)print *,'Fatal error in GLVOLU'
 334            ingoto=0
 335          endif
 336 #endif
 337          NMEC = NMEC +1
 338          LMEC(NMEC) = 4
 339 *
 340 #if !defined(CERNLIB_USRJMP)
 341    50    CALL GUSWIM (CHARG1, STEP, VECT, VOUT)
 342 #endif
 343 #if defined(CERNLIB_USRJMP)
 344    50    CALL JUMPT4(JUSWIM, CHARG1, STEP, VECT, VOUT)
 345 #endif
 346 *
 347 *  ** When near to boundary, take proper action (cut-step,crossing...)
 348 *
 349          IF(STEP.GE.SAFETY)THEN
 350             INEAR = 0
 351             IF (IGNEXT.NE.0) THEN
 352                DO 60 I = 1,3
 353                   VNEXT(I+3) = VECT(I+3)
 354                   VNEXT(I) = VECT(I) +SNEXT*VECT(I+3)
 355    60          CONTINUE
 356                DO 70 I = 1,3
 357                   IF (ABS(VOUT(I)-VNEXT(I)).GT.EPSIL) GO TO 80
 358    70          CONTINUE
 359                INEAR = 1
 360             ENDIF
 361 *
 362    80       CALL GINVOL (VOUT, ISAME)
 363             IF (ISAME.EQ.0)THEN
 364                IF ((INEAR.NE.0).OR.(STEP.LT.EPSIL)) THEN
 365                   INWVOL = 2
 366                   NMEC = NMEC +1
 367                   LMEC(NMEC) = 1
 368                ELSE
 369 *              Cut step
 370                   STEP = 0.5*STEP
 371                   IF (LMEC(NMEC).NE.24) THEN
 372                      NMEC = NMEC +1
 373                      LMEC(NMEC) = 24
 374                   ENDIF
 375                   GO TO 50
 376                ENDIF
 377             ENDIF
 378          ENDIF
 379 *
 380          DO 90 I = 1,6
 381             VECT(I) = VOUT(I)
 382    90    CONTINUE
 383 *
 384       ENDIF
 385 *
 386 * *** Correct the step due to multiple scattering
 387       IF (IMULL.NE.0) THEN
 388          STMULS = STEP
 389          CORR=0.0001*CHARG2*(STEP/RADL)*(GETOT/(VECT(7)*VECT(7)))**2
 390          IF (CORR.GT.0.25) CORR = 0.25
 391          STEP  = (1.+CORR)*STEP
 392       ENDIF
 393 *
 394       SLENG = SLENG + STEP
 395 *
 396 * *** Generate Cherenkov photons if required
 397 *
 398       IF(IMCKOV.EQ.1) THEN
 399          CALL GGCKOV
 400          NMEC=NMEC+1
 401          LMEC(NMEC)=105
 402       ENDIF
 403 *
 404 * *** apply energy loss : find the kinetic energy corresponding
 405 *      to the new stopping range = stopmx - step
 406 *
 407       IF (ILOSL.NE.0) THEN
 408          NMEC = NMEC +1
 409          LMEC(NMEC) = 3
 410          STOPRG = STOPP - STEP*RMASS*CHARG2
 411          IF (STOPRG.LE.STOPC) THEN
 412             STEP = STOPMX
 413             GO TO 100
 414          ENDIF
 415          IF(XCOEF1.NE.0.) THEN
 416             EKIPR = XCOEF1*(XCOEF3+STOPRG*(2.*XCOEF2+STOPRG))
 417          ELSE
 418             EKIPR = XCOEF2*STOPRG+XCOEF3
 419          ENDIF
 420          DEMEAN=GEKIN - EKIPR/RMASS
 421          IF(DEMEAN.LE.5.*GEKIN*EPSMAC) THEN
 422             DEMEAN=(GEKRT1*Q(JLOSS+IEKBIN)+GEKRAT*Q(JLOSS+IEKBIN+1))
 423      +             *STEP*CHARG2
 424          ENDIF
 425 *
 426 *        fluctuations : differ from that of 'ordinary' hadrons
 427 *
 428          IF (ILOSS.EQ.4.OR.IEKBIN.LE.IKCUT+1) THEN
 429             DESTEP = DEMEAN
 430          ELSE
 431 *
 432 *     Charge exchange fluctuations + Gaussian 'Landau' fluctuations
 433 *           (it is the same for ILOSS=1,2,3 !)
 434 *
 435             SIGMA2=CNORM*CHARG1*(1.-CHARG1/CHARGE)
 436             SIGMA2=MAX(SIGMA2,0.)
 437             TA = RMASS*GEKIN
 438             TAM=TA/AMU
 439             SIGMA2=SIGMA2+2.+TAM*(2.+TAM)
 440 *
 441             SIGMA2=FACFLU*CHARG2*STEP*SIGMA2
 442             IF(SIGMA2.GT.0.0) THEN
 443                 SIGMA=SQRT(SIGMA2)
 444             ELSE
 445                 SIGMA= 0.0
 446             END IF
 447 *
 448 *     Check if we are in 'Gaussian' regime ...
 449 *
 450             CAPPA=(1.+TAM)/(TAM*(2.+TAM)*EMASS)
 451             CAPPA=0.5*CAPPA**2*FACFLU*CHARG2*STEP
 452 *
 453 *     ... if not , correct SIGMA !
 454
 455             IF( (CAPPA.LT.10.) .AND. (CAPPA.GT.0.0) ) THEN
 456                SIGMA=SIGMA/(0.97+0.03*SQRT(10./CAPPA))
 457             ENDIF
 458 *
 459             CALL GRNDM(RNDM,2)
 460             DEFLUC=SIGMA*SIN(TWOPI*RNDM(1))*SQRT(-2.*LOG(RNDM(2)))
 461             DESTEP=DEMEAN+DEFLUC
 462          ENDIF
 463 *
 464 *     protection against negative destep
 465 *
 466          IF(DESTEP.LT.0.) DESTEP=DEMEAN
 467 *                          IF (DESTEP.LT.0.) DESTEP = 0.
 468          GEKINT = GEKIN -DESTEP
 469          IF (GEKINT.LE.(1.01*CUTHAD)) GO TO 100
 470          DESTEL = DESTEP
 471          GEKIN  = GEKINT
 472          GETOT  = GEKIN +AMASS
 473          VECT(7)= SQRT((GETOT+AMASS)*GEKIN)
 474          CALL GEKBIN
 475       ENDIF
 476 *
 477 * *** Apply multiple scattering.
 478 *
 479       IF (IMULL.NE.0) THEN
 480          NMEC = NMEC +1
 481          LMEC(NMEC) = 2
 482 *      check charge dependence ...........!!!!!!!  (later..)
 483          CALL GMULTS
 484       ENDIF
 485 *
 486 * *** Update time of flight
 487 *
 488       SUMLIF = SUMLIF -STEP*AMASS/VECT(7)
 489       TOFG   = TOFG +STEP*GETOT/(VECT(7)*CLIGHT)
 490       IF (TOFG.GE.TOFMAX) THEN
 491          ISTOP = 4
 492          NMEC  = NMEC +1
 493          LMEC(NMEC) = 22
 494          GO TO 999
 495       ENDIF
 496 *
 497 * *** Update interaction probabilities
 498 *
 499       IF (IHADR.GT.0) ZINTHA = ZINTHA*(1.-STEP/SHADR)
 500       IF (IDRAY.GT.0) ZINTDR = ZINTDR -STEP/STEPDR
 501 *
 502       GO TO 110
 503 *
 504 *  **   Special treatment for overstopped tracks
 505 *
 506   100 DESTEP = GEKIN
 507       DESTEL = DESTEP
 508       GEKIN  = 0.
 509       GETOT  = AMASS
 510       VECT(7)= 0.
 511       INWVOL = 0
 512       ISTOP  = 2
 513       NMEC = NMEC + 1
 514       LMEC(NMEC) = 30
 515       IF (IHADR.EQ.0) GO TO 999
 516       IPROC = 12
 517 *
 518 * *** apply slected process if any
 519 *
 520   110 IF (IPROC.EQ.0) GO TO 999
 521       NMEC = NMEC +1
 522       LMEC(NMEC) = IPROC
 523 *
 524 *  **   Hadron interaction ?
 525 *
 526       IF (IPROC.EQ.12) THEN
 527 #if !defined(CERNLIB_USRJMP)
 528          CALL GUHADR
 529 #endif
 530 #if defined(CERNLIB_USRJMP)
 531          CALL JUMPT0(JUHADR)
 532 #endif
 533 *   *   Check time cut-off for decays at rest
 534          IF (LMEC(NMEC).EQ.5) THEN
 535             TOFG   = TOFG +SUMLIF/CLIGHT
 536             SUMLIF = 0.
 537             IF (TOFG.GE.TOFMAX) THEN
 538                NGKINE = 0
 539                ISTOP  = 4
 540                LMEC(NMEC) = 22
 541             ENDIF
 542          ENDIF
 543 *
 544 *  **   Delta-ray ?
 545 *
 546       ELSE IF (IPROC.EQ.10) THEN
 547          CALL GDRAY
 548       ENDIF
 549   999 END