gtgama.F

Go to the documentation of this file.

      SUBROUTINE gtgama
C.
C.    ******************************************************************
C.    *                                                                *
C.    *   Photon track. Computes step size and propagates particle     *
C.    *    through step.                                               *
C.    *                                                                *
C.    *   ==>Called by : GTRACK                                        *
C.    *      Authors    R.Brun, F.Bruyant L.Urban ********             *
C.    *                                                                *
C.    ******************************************************************
C.
#include "geant321/gcbank.inc"
#include "geant321/gccuts.inc"
#include "geant321/gcjloc.inc"
#include "geant321/gconsp.inc"
#include "geant321/gcphys.inc"
#include "geant321/gcstak.inc"
#include "geant321/gctmed.inc"
#include "geant321/gcmulo.inc"
#include "geant321/gctrak.inc"
#include "geant321/gcunit.inc"

      parameter(epsmac=1.e-6)
      DOUBLE PRECISION one,xcoef1,xcoef2,xcoef3,zero

      parameter(one=1,zero=0)
      parameter(epcut=1.022e-3)
C.
*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      iaban = nint(dphys1)
*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>      
C.    ------------------------------------------------------------------
*
* *** Particle below energy threshold ?  Short circuit
*
*
      IF (gekin.LE.cutgam) goto 998
*
* *** Update local pointers if medium has changed
*
      IF(iupd.EQ.0)THEN
         iupd  = 1
         jphot = lq(jma-6)
         jcomp = lq(jma-8)
         jpair = lq(jma-10)
         jpfis = lq(jma-12)
         jrayl = lq(jma-13)
      ENDIF
*
* *** Compute current step size
*
      iproc  = 103
      step   = stemax
      gekrt1 = 1 .-gekrat
*
*  **   Step limitation due to pair production ?
*
      IF (getot.GT.epcut) THEN
         IF (ipair.GT.0) THEN
            steppa = gekrt1*q(jpair+iekbin) +gekrat*q(jpair+iekbin+1)
            spair  = steppa*zintpa
            IF (spair.LT.step) THEN
               step  = spair
               iproc = 6
            ENDIF
         ENDIF
      ENDIF
*
*  **   Step limitation due to Compton scattering ?
*
      IF (icomp.GT.0) THEN
         stepco = gekrt1*q(jcomp+iekbin) +gekrat*q(jcomp+iekbin+1)
         scomp  = stepco*zintco
         IF (scomp.LT.step) THEN
            step  = scomp
            iproc = 7
         ENDIF
      ENDIF
*
*  **   Step limitation due to photo-electric effect ?
*
      IF (gekin.LT.0.4) THEN
         IF (iphot.GT.0) THEN
            stepph = gekrt1*q(jphot+iekbin) +gekrat*q(jphot+iekbin+1)
            sphot  = stepph*zintph
            IF (sphot.LT.step) THEN
               step  = sphot
               iproc = 8
            ENDIF
         ENDIF
      ENDIF
*
*  **   Step limitation due to photo-fission ?
*
      IF (jpfis.GT.0) THEN
         steppf = gekrt1*q(jpfis+iekbin) +gekrat*q(jpfis+iekbin+1)
         spfis  = steppf*zintpf
         IF (spfis.LT.step) THEN
            step  = spfis
            iproc = 23
         ENDIF
      ENDIF
*
*  **   Step limitation due to Rayleigh scattering ?
*
      IF (irayl.GT.0) THEN
         IF (gekin.LT.0.01) THEN
            stepra = gekrt1*q(jrayl+iekbin) +gekrat*q(jrayl+iekbin+1)
            srayl  = stepra*zintra
            IF (srayl.LT.step) THEN
               step  = srayl
               iproc = 25
            ENDIF
         ENDIF
      ENDIF
*
      IF (step.LT.0.) step = 0.
*
*  **   Step limitation due to geometry ?
*
      IF (step.GE.safety) THEN
         CALL gtnext
         IF (ignext.NE.0) THEN
            step   = snext + prec
            inwvol= 2
            iproc = 0
            nmec =  1
            lmec(1)=1
         ENDIF
*
*        Update SAFETY in stack companions, if any
         IF (iq(jstak+3).NE.0) THEN
            DO 10 ist = iq(jstak+3),iq(jstak+1)
               jst = jstak +3 +(ist-1)*nwstak
               q(jst+11) = safety
   10       CONTINUE
            iq(jstak+3) = 0
         ENDIF
*
      ELSE
         iq(jstak+3) = 0
      ENDIF
*
* *** Linear transport
*
      IF (inwvol.EQ.2) THEN
         DO 20 i = 1,3
            vectmp  = vect(i) +step*vect(i+3)
            IF(vectmp.EQ.vect(i)) THEN
*
* *** Correct for machine precision
*
               IF(vect(i+3).NE.0.) THEN
                  vectmp = vect(i)+abs(vect(i))*sign(1.,vect(i+3))*
     +            epsmac
                  IF(nmec.GT.0) THEN
                     IF(lmec(nmec).EQ.104) nmec=nmec-1
                  ENDIF
                  nmec=nmec+1
                  lmec(nmec)=104
#ifdef G3DEBUG
                  WRITE(chmail, 10000)
                  CALL gmail(0,0)
                  WRITE(chmail, 10100) gekin, numed, step, snext
                  CALL gmail(0,0)
10000 FORMAT(' Boundary correction in GTGAMA: ',
     +       '    GEKIN      NUMED       STEP      SNEXT')
10100 FORMAT(31x,e10.3,1x,i10,1x,e10.3,1x,e10.3,1x)
#endif
               ENDIF
            ENDIF
            vect(i) = vectmp
   20    CONTINUE
      ELSE
         DO 30 i = 1,3
            vect(i)  = vect(i) +step*vect(i+3)
   30    CONTINUE
      ENDIF
*
      sleng = sleng +step
*
* *** Update time of flight
*
      tofg = tofg +step/clight
*
* *** Update interaction probabilities
*
      IF (getot.GT.epcut) THEN
         IF (ipair.GT.0) zintpa = zintpa -step/steppa
      ENDIF
      IF (icomp.GT.0)    zintco = zintco -step/stepco
      IF (gekin.LT.0.4) THEN
         IF (iphot.GT.0) zintph = zintph -step/stepph
      ENDIF
      IF (jpfis.GT.0)    zintpf = zintpf -step/steppf
      IF (irayl.GT.0) THEN
         IF (gekin.LT.0.01) zintra = zintra -step/stepra
      ENDIF
*
      IF (iproc.EQ.0) go to 999
      nmec = 1
      lmec(1) = iproc
*
*  ** Pair production ?
*
      IF (iproc.EQ.6) THEN
         CALL gpairg
*
*  ** Compton scattering ?
*
      ELSE IF (iproc.EQ.7) THEN
         CALL gcomp
*
*  ** Photo-electric effect ?
*
      ELSE IF (iproc.EQ.8) THEN
*
        IF ((iaban.NE.0).AND.(gekin.LE.0.001))  THEN
*           Calculate range of the photoelectron ( with kin. energy Ephot)        
            jcoef = lq(jma-17)
         IF(gekrat.LT.0.7) THEN
            i1 = max(iekbin-1,1)
         ELSE
            i1 = min(iekbin,nekbin-1)
         ENDIF
         i1 = 3*(i1-1)+1
         xcoef1 = q(jcoef+i1)
         xcoef2 = q(jcoef+i1+1)
         xcoef3 = q(jcoef+i1+2)
         IF(xcoef1.NE.0.) THEN
            stopmx = -xcoef2+sign(one,xcoef1)*sqrt(xcoef2**2 - (xcoef3-
     +      gekin/xcoef1))
         ELSE
            stopmx = - (xcoef3-gekin)/xcoef2
         ENDIF
*
*        DO NOT call GPHOT if this (overestimated) range is smaller
*                than SAFETY
*
         IF (stopmx.LE.safety) goto 998
        ENDIF
 
         CALL gphot
*
*  ** Rayleigh effect ?
*
      ELSE IF (iproc.EQ.25) THEN
         CALL grayl
*
*  ** Photo-fission ?
*
      ELSE IF (iproc.EQ.23) THEN
         CALL gpfis
*
      ENDIF
*
         goto 999
998      destep = gekin
         gekin  = 0.
         getot  = 0.
         vect(7)= 0.
         istop  = 2
         nmec   = 1
         lmec(1)= 30
  999 END