GEANT321/gtrak/gtmuon.F

   1 #include "geant321/pilot.h"
   2 *CMZ :  3.21/02 29/03/94  15.41.24  by  S.Giani
   3 *-- Author :
   4       SUBROUTINE GTMUON
   5 C.
   6 C.    ******************************************************************
   7 C.    *                                                                *
   8 C.    *   Muon track. Computes step size and propagates particle       *
   9 C.    *    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/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
  41 #endif
  42 #if defined(CERNLIB_SINGLE)
  43       PARAMETER (EPSMAC=1.E-11)
  44 #endif
  45       PARAMETER (ONE=1)
  46       REAL VNEXT(6)
  47       SAVE IKCUT,STOPC
  48 C.
  49 C.    ------------------------------------------------------------------
  50 *
  51 * *** Particle below energy threshold ? short circuit
  52 *
  53       IF (GEKIN.LE.CUTMUO) GO TO 100
  54 *
  55 * *** Update local pointers if medium has changed
  56       IF (IUPD.EQ.0) THEN
  57          IUPD  = 1
  58          JLOSS = LQ(JMA-2)
  59          JBREM = LQ(JMA-9)
  60          JPAIR = LQ(JMA-10)
  61          JDRAY = LQ(JMA-11)
  62          JMUNU = LQ(JMA-14)
  63          JRANG = LQ(JMA-16)
  64          JCOEF = LQ(JMA-18)
  65          JMULOF= LQ(JTM-2)
  66          IF(IMCKOV.EQ.1) THEN
  67             JTCKOV = LQ(JTM-3)
  68             JABSCO = LQ(JTCKOV-1)
  69             JEFFIC = LQ(JTCKOV-2)
  70             JINDEX = LQ(JTCKOV-3)
  71             JCURIN = LQ(JTCKOV-4)
  72             NPCKOV = Q(JTCKOV+1)
  73          ENDIF
  74          OMCMOL= Q(JPROB+21)
  75          CHCMOL= Q(JPROB+25)
  76          IKCUT = Q(JMULOF+NEK1+1)
  77          STOPC = Q(JMULOF+NEK1+2)
  78          IF(ISTRA.GT.0) THEN
  79             JTSTRA = LQ(JMA-19)
  80             JTSTCO = LQ(JTSTRA-1)
  81             JTSTEN = LQ(JTSTRA-2)
  82 #if defined(CERNLIB_ASHO)
  83             IF(ISTRA.EQ.2) THEN
  84                JTASHO = LQ(JMA-20)
  85             ENDIF
  86 #endif
  87          ENDIF
  88       ENDIF
  89 *
  90 * *** Compute current step size
  91 *
  92       STEP   = STEMAX
  93       IPROC  = 103
  94       GEKRT1 = 1. -GEKRAT
  95       IEK1   = IEKBIN+NEK1
  96       IEK2   = IEKBIN+2*NEK1
  97 *
  98 *  **   Step limitation due to bremsstrahlung ?
  99 *
 100       IF (IBREM.GT.0) THEN
 101          STEPBR = GEKRT1*Q(JBREM+IEK2) +GEKRAT*Q(JBREM+IEK2+1)
 102          SBREM  = STEPBR*ZINTBR
 103          IF (SBREM.LT.STEP) THEN
 104             STEP  = SBREM
 105             IPROC = 9
 106          ENDIF
 107       ENDIF
 108 *
 109 *  **   Step limitation due to pair production ?
 110 *
 111       IF (IPAIR.GT.0) THEN
 112          STEPPA = GEKRT1*Q(JPAIR+IEK1) +GEKRAT*Q(JPAIR+IEK1+1)
 113          SPAIR  = STEPPA*ZINTPA
 114          IF (SPAIR.LT.STEP) THEN
 115             STEP  = SPAIR
 116             IPROC = 6
 117          ENDIF
 118       ENDIF
 119 *
 120 *  **   Step limitation due to decay ?
 121 *
 122       IF (IDCAY.NE.0) THEN
 123          SDCAY = SUMLIF*VECT(7)/AMASS
 124          IF (SDCAY.LT.STEP) THEN
 125             STEP  = SDCAY
 126             IPROC = 5
 127          ENDIF
 128       ENDIF
 129 *
 130 *  **   Step limitation due to delta-ray ?
 131 *
 132       IF (IDRAY.GT.0) THEN
 133          STEPDR = GEKRT1*Q(JDRAY+IEK2) +GEKRAT*Q(JDRAY+IEK2+1)
 134          SDRAY  = STEPDR*ZINTDR
 135          IF (SDRAY.LT.STEP) THEN
 136             STEP  = SDRAY
 137             IPROC = 10
 138          ENDIF
 139       ENDIF
 140 *
 141 *  **   Step limitation due to nuclear interaction ?
 142 *
 143       IF (IMUNU.GT.0) THEN
 144          IF(GEKIN.GE.5.)THEN
 145             STEPMU = GEKRT1*Q(JMUNU+IEKBIN) +GEKRAT*Q(JMUNU+IEKBIN+1)
 146             SMUNU = STEPMU*ZINTMU
 147             IF (SMUNU.LT.STEP) THEN
 148                STEP  = SMUNU
 149                IPROC = 21
 150             ENDIF
 151          ELSE
 152             STEPMU = BIG
 153          ENDIF
 154       ENDIF
 155 *
 156       IF (STEP.LE.0.) THEN
 157          STEP = 0.
 158          GO TO 90
 159       ENDIF
 160 *
 161 *  **   Step limitation due to energy-loss,multiple scattering
 162 *             or magnetic field ?
 163 *
 164       IF (JMULOF.NE.0) THEN
 165          SMULOF  = GEKRT1*Q(JMULOF+IEKBIN) +GEKRAT*Q(JMULOF+IEKBIN+1)
 166          IF (SMULOF.LT.STEP) THEN
 167             STEP  = SMULOF
 168             IPROC = 0
 169          ENDIF
 170       ENDIF
 171 *
 172 *  **   Step limitation due to geometry ?
 173 *
 174       IF (STEP.GE.0.95*SAFETY) THEN
 175          CALL GTNEXT
 176          IF (IGNEXT.NE.0) THEN
 177             STEP  = SNEXT + PREC
 178             IPROC = 0
 179          ENDIF
 180 *
 181 *        Update SAFETY in stack companions, if any
 182          IF (IQ(JSTAK+3).NE.0) THEN
 183             DO 10 IST = IQ(JSTAK+3),IQ(JSTAK+1)
 184                JST    = JSTAK + 3 + (IST-1)*NWSTAK
 185                Q(JST+11) = SAFETY
 186    10       CONTINUE
 187             IQ(JSTAK+3) = 0
 188          ENDIF
 189       ELSE
 190          IQ(JSTAK+3) = 0
 191       ENDIF
 192 *
 193 * *** Linear transport when no field or very short step
 194 *
 195       IF (IFIELD.EQ.0.OR.STEP.LE.PREC) THEN
 196 *
 197          IF (IGNEXT.NE.0) THEN
 198             DO 20 I = 1,3
 199                VECTMP  = VECT(I) +STEP*VECT(I+3)
 200                IF(VECTMP.EQ.VECT(I)) THEN
 201 *
 202 * *** Correct for machine precision
 203 *
 204                   IF(VECT(I+3).NE.0.) THEN
 205                      VECTMP =
 206      +               VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))*EPSMAC
 207                      IF(NMEC.GT.0) THEN
 208                         IF(LMEC(NMEC).EQ.104) NMEC=NMEC-1
 209                      ENDIF
 210                      NMEC=NMEC+1
 211                      LMEC(NMEC)=104
 212 #if defined(CERNLIB_DEBUG)
 213                      WRITE(CHMAIL, 10000)
 214                      CALL GMAIL(0,0)
 215                      WRITE(CHMAIL, 10100) GEKIN, NUMED, STEP, SNEXT
 216                      CALL GMAIL(0,0)
 217 10000 FORMAT(' Boundary correction in GTMUON: ',
 218      +       '    GEKIN      NUMED       STEP      SNEXT')
 219 10100 FORMAT(31X,E10.3,1X,I10,1X,E10.3,1X,E10.3,1X)
 220 #endif
 221                   ENDIF
 222                ENDIF
 223                VECT(I) = VECTMP
 224    20       CONTINUE
 225             INWVOL = 2
 226             NMEC = NMEC +1
 227             LMEC(NMEC) = 1
 228          ELSE
 229             DO 30 I = 1,3
 230                VECT(I)  = VECT(I) +STEP*VECT(I+3)
 231    30       CONTINUE
 232          ENDIF
 233       ELSE
 234 *
 235 * ***   otherwise, swim particle in magnetic field
 236 *
 237 #if !defined(CERNLIB_OLD)
 238          if(mycoun.gt.1.and.nfmany.gt.0.and.step.ge.safety)then
 239            nlevel=manyle(nfmany)
 240            do 99 i=1,nlevel
 241              names(i)=manyna(nfmany,i)
 242              number(i)=manynu(nfmany,i)
 243  99        continue
 244            call glvolu(nlevel,names,number,ier)
 245            if(ier.ne.0)print *,'Fatal error in GLVOLU'
 246            ingoto=0
 247          endif
 248 #endif
 249          NMEC = NMEC +1
 250          LMEC(NMEC) = 4
 251 *
 252 #if !defined(CERNLIB_USRJMP)
 253    40    CALL GUSWIM (CHARGE, STEP, VECT, VOUT)
 254 #endif
 255 #if defined(CERNLIB_USRJMP)
 256    40    CALL JUMPT4(JUSWIM, CHARGE, STEP, VECT, VOUT)
 257 #endif
 258 *
 259 *  ** When near to boundary, take proper action (cut-step,crossing...)
 260 *
 261          IF(STEP.GE.SAFETY)THEN
 262             INEAR = 0
 263             IF (IGNEXT.NE.0) THEN
 264                DO 50 I = 1,3
 265                   VNEXT(I+3) = VECT(I+3)
 266                   VNEXT(I) = VECT(I) +SNEXT*VECT(I+3)
 267    50          CONTINUE
 268                DO 60 I = 1,3
 269                   IF (ABS(VOUT(I)-VNEXT(I)).GT.EPSIL) GO TO 70
 270    60          CONTINUE
 271                INEAR = 1
 272             ENDIF
 273 *
 274    70       CALL GINVOL (VOUT, ISAME)
 275             IF (ISAME.EQ.0)THEN
 276                IF ((INEAR.NE.0).OR.(STEP.LT.EPSIL)) THEN
 277                   INWVOL = 2
 278                   NMEC = NMEC +1
 279                   LMEC(NMEC) = 1
 280                ELSE
 281 *              Cut step
 282                   STEP = 0.5*STEP
 283                   IF (LMEC(NMEC).NE.24) THEN
 284                      NMEC = NMEC +1
 285                      LMEC(NMEC) = 24
 286                   ENDIF
 287                   GO TO 40
 288                ENDIF
 289             ENDIF
 290          ENDIF
 291 *
 292          DO 80 I = 1,6
 293             VECT(I) = VOUT(I)
 294    80    CONTINUE
 295 *
 296       ENDIF
 297 *
 298 * *** Correct the step due to multiple scattering
 299       IF (IMULL.NE.0) THEN
 300          STMULS = STEP
 301          CORR=0.0001*(STEP/RADL)*(GETOT/(VECT(7)*VECT(7)))**2
 302          IF (CORR.GT.0.25) CORR = 0.25
 303          STEP  = (1.+CORR)*STEP
 304       ENDIF
 305 *
 306       SLENG = SLENG + STEP
 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          LMEC(NMEC) = 3
 324          IF(GEKRAT.LT.0.7) THEN
 325             I1 = MAX(IEKBIN-1,1)
 326          ELSE
 327             I1 = MIN(IEKBIN,NEKBIN-1)
 328          ENDIF
 329          I1 = 3*(I1-1)+1
 330          XCOEF1 = Q(JCOEF+I1)
 331          XCOEF2 = Q(JCOEF+I1+1)
 332          XCOEF3 = Q(JCOEF+I1+2)
 333          IF(XCOEF1.NE.0.) THEN
 334             STOPMX = -XCOEF2+SIGN(ONE,XCOEF1)*SQRT(XCOEF2**2 - (XCOEF3-
 335      +      GEKIN/XCOEF1))
 336          ELSE
 337             STOPMX = - (XCOEF3-GEKIN)/XCOEF2
 338          ENDIF
 339          STOPRG = STOPMX - STEP
 340          IF (STOPRG.LT.STOPC) THEN
 341             STEP = STOPMX - STOPC
 342             GO TO 100
 343          ENDIF
 344 *
 345          IF(XCOEF1.NE.0.) THEN
 346             DEMEAN=GEKIN-XCOEF1*(XCOEF3+STOPRG*(2.*XCOEF2+STOPRG))
 347          ELSE
 348             DEMEAN=GEKIN-XCOEF2*STOPRG-XCOEF3
 349          ENDIF
 350          IF(DEMEAN.LE.5.*GEKIN*EPSMAC) THEN
 351             DEMEAN=(GEKRT1*Q(JLOSS+IEKBIN)+GEKRAT*Q(JLOSS+IEKBIN+1))
 352      +      *STEP
 353          ENDIF
 354          IF (ILOSS.EQ.4.OR.IEKBIN.LE.IKCUT+1) THEN
 355             DESTEP = DEMEAN
 356          ELSE
 357             DEMS = DEMEAN
 358             CALL GFLUCT(DEMS,DESTEP)
 359          ENDIF
 360          IF (DESTEP.LT.0.) DESTEP = 0.
 361          GEKINT = GEKIN -DESTEP
 362          IF (GEKINT.LE.(1.01*CUTMUO)) GO TO 100
 363          DESTEL = DESTEP
 364          GEKIN  = GEKINT
 365          GETOT  = GEKIN +AMASS
 366          VECT(7)= SQRT((GETOT+AMASS)*GEKIN)
 367          CALL GEKBIN
 368       ENDIF
 369 *
 370 * *** Apply multiple scattering.
 371 *
 372       IF (IMULL.NE.0) THEN
 373          NMEC   = NMEC +1
 374          LMEC(NMEC) = 2
 375          CALL GMULTS
 376       ENDIF
 377 *
 378 * *** Update time of flight
 379 *
 380       SUMLIF = SUMLIF -STEP*AMASS/VECT(7)
 381       TOFG = TOFG +STEP*GETOT/(VECT(7)*CLIGHT)
 382       IF (TOFG.GE.TOFMAX) THEN
 383          ISTOP = 4
 384          NMEC  = NMEC +1
 385          LMEC(NMEC) = 22
 386          GO TO 999
 387       ENDIF
 388 *
 389 * *** Update interaction probabilities
 390 *
 391       IF (IBREM.GT.0) ZINTBR = ZINTBR -STEP/STEPBR
 392       IF (IPAIR.GT.0) ZINTPA = ZINTPA -STEP/STEPPA
 393       IF (IDRAY.GT.0) ZINTDR = ZINTDR -STEP/STEPDR
 394       IF (IMUNU.GT.0) ZINTMU = ZINTMU -STEP/STEPMU
 395 *
 396 * ***   otherwise, apply the selected process if any
 397 *
 398    90 IF (IPROC.EQ.0) GO TO 999
 399       NMEC = NMEC +1
 400       LMEC(NMEC) = IPROC
 401 *
 402 *  **   Bremsstrahlung ?
 403 *
 404       IF (IPROC.EQ.9) THEN
 405          CALL GBREMM
 406 *
 407 *  **   Pair production ?
 408 *
 409       ELSE IF (IPROC.EQ.6) THEN
 410          CALL GPAIRM
 411 *
 412 *  **   Decay ?
 413 *
 414       ELSE IF (IPROC.EQ.5) THEN
 415          ISTOP = 1
 416          CALL GDECAY
 417 *
 418 *  **   Delta-ray ?
 419 *
 420       ELSE IF (IPROC.EQ.10) THEN
 421          CALL GDRAY
 422 *
 423 *  **   Nuclear interaction ?
 424 *
 425       ELSE IF (IPROC.EQ.21) THEN
 426          CALL GMUNU
 427       ENDIF
 428       GO TO 999
 429 *
 430 * *** Special treatment for overstopped tracks
 431 *
 432   100 DESTEP = GEKIN
 433       DESTEL = DESTEP
 434       GEKIN  = 0.
 435       GETOT  = AMASS
 436       VECT(7)= 0.
 437       INWVOL = 0
 438       NMEC   = NMEC +1
 439       LMEC(NMEC) = 30
 440       IF (IDCAY.EQ.0) THEN
 441          ISTOP = 2
 442       ELSE
 443          NMEC   = NMEC +1
 444          TOFG   = TOFG +SUMLIF/CLIGHT
 445          SUMLIF = 0.
 446          IF (TOFG.GE.TOFMAX) THEN
 447             ISTOP = 4
 448             LMEC(NMEC) = 22
 449             GO TO 999
 450          ENDIF
 451          LMEC(NMEC) = 5
 452          ISTOP = 1
 453          CALL GDECAY
 454       ENDIF
 455   999 IF(NGPHOT.GT.0) THEN
 456          IF(ITCKOV.EQ.2.AND.ISTOP.EQ.0) THEN
 457 *
 458 *  The muon has produced Cerenkov photons and it is still alive
 459 *  we put it in the stack and we let the photons to be tracked
 460             NGKINE = NGKINE+1
 461             GKIN(1,NGKINE) = VECT(4)*VECT(7)
 462             GKIN(2,NGKINE) = VECT(5)*VECT(7)
 463             GKIN(3,NGKINE) = VECT(6)*VECT(7)
 464             GKIN(4,NGKINE) = GETOT
 465             GKIN(5,NGKINE) = IPART
 466             TOFD(NGKINE) = 0.
 467             ISTOP = 1
 468 c----put position as well
 469             GPOS(1,NGKINE)=VECT(1)
 470             GPOS(2,NGKINE)=VECT(2)
 471             GPOS(3,NGKINE)=VECT(3)
 472          ENDIF
 473       ENDIF
 474       END