GEANT321/gtrak/gthadr.F

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