GEANT321/gphys/gcomp.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:23  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.21  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GCOMP
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *  Simulates photon-electron COMPTON scattering.                 *
  17 C.    *                                                                *
  18 C.    *  The scattered photon energy is sampled using the quantum-     *
  19 C.    *  mechanical KLEIN-NISHINA formula. For this, the random       *
  20 C.    *  number techniques of BUTCHER and MESSEL(NUC. PHYS.20(1960),   *
  21 C.    *  15) are employed.                                             *
  22 C.    *  NOTE :                                                        *
  23 C.    *  (1) Effects due to binding of atomic electrons are            *
  24 C.    *      ignored(recoil electron energy assumed large compared     *
  25 C.    *      with binding energy).                                     *
  26 C.    *                                                                *
  27 C.    *    ==>Called by : GTGAMA                                       *
  28 C.    *       Authors    G.Patrick, L.Urban  *********                 *
  29 C.    *                                                                *
  30 C.    ******************************************************************
  31 C.
  32 #include "geant321/gcphys.inc"
  33 #include "geant321/gctrak.inc"
  34 #include "geant321/gcking.inc"
  35 #include "geant321/gconsp.inc"
  36 #include "geant321/gccuts.inc"
  37       DIMENSION PGAM(3)
  38       REAL RNDM(5)
  39       LOGICAL ROTATE
  40 C.
  41 C.    ------------------------------------------------------------------
  42 C.
  43       KCASE  = NAMEC(7)
  44       EGAM1  = VECT(7)
  45       EZERO  = EGAM1/EMASS
  46       EMINI  = 1.+2.*EZERO
  47       EMIN   = 1./EMINI
  48       DSIG1  = LOG(EMINI)
  49       DSIG2  = 0.5*(1.-EMIN*EMIN)
  50       DSIGT  = DSIG1+DSIG2
  51 C
  52 C             Decide which part of F(E)=1/E+E to sample from.
  53 C
  54    10 CALL GRNDM(RNDM,3)
  55       IF (DSIG1.LT.DSIGT*RNDM(1))THEN
  56 C
  57 C             Sample from F2(E) distribution.
  58 C
  59          BRD    = RNDM(2)
  60          CALL GRNDM(RNDM(4),1)
  61          IF (EZERO.GE.(EZERO+1.)*RNDM(4))THEN
  62             CALL GRNDM(RNDM(5),1)
  63             BRD    = MAX(BRD,RNDM(5))
  64          ENDIF
  65 C
  66          BR     = EMIN+(1.-EMIN)*BRD
  67       ELSE
  68          BR     = EMIN*EXP(DSIG1*RNDM(2))
  69       ENDIF
  70 C
  71 C             Scattered photon energy.
  72 C
  73       EGAM2  = BR*EGAM1
  74 C
  75 C             Calculate rejection function G(E).
  76 C
  77       T      = EMASS*(1.-BR)/EGAM2
  78       SINTH  = MAX(0.,T*(2.-T))
  79       REJ    = 1.0-(BR*SINTH)/(1.+BR*BR)
  80       IF (RNDM(3).GT.REJ)                        GO TO 10
  81 C
  82 C             Successful sampling of scattered photon.
  83 C
  84 C             CUTS ON ENERGY THRESHOLDS ?
  85 C
  86       TEL     = EGAM1-EGAM2
  87       IF((EGAM2.LE.CUTGAM).AND.(TEL.LE.CUTELE)) THEN
  88          ISTOP = 2
  89          VECT(7) = 0.
  90          GEKIN = 0.
  91          GETOT = 0.
  92          NGKINE = 0.
  93          DESTEP = DESTEP + EGAM1
  94          RETURN
  95       ENDIF
  96 C
  97 C             Generate photon angles with respect to a Z-axis
  98 C             defined along the parent photon.
  99 C             PHI is generated isotropically
 100 C
 101       SINTH  = SQRT(SINTH)
 102       COSTH  = 1.-T
 103       CALL GRNDM(RNDM,1)
 104       PHI    = TWOPI*RNDM(1)
 105       COSPHI = COS(PHI)
 106       SINPHI = SIN(PHI)
 107 C
 108 C             Polar co-ordinates to momentum components.
 109 C
 110       PGAM(1) = EGAM2*SINTH*COSPHI
 111       PGAM(2) = EGAM2*SINTH*SINPHI
 112       PGAM(3) = EGAM2*COSTH
 113 C
 114 C             Momentum vector of recoil electron.
 115 C
 116       NGKINE = 1
 117       EEL    = TEL + EMASS
 118       GKIN(1,1) = -PGAM(1)
 119       GKIN(2,1) = -PGAM(2)
 120       GKIN(3,1) = EGAM1-PGAM(3)
 121       GKIN(4,1)=EEL
 122       GKIN(5,1)=3
 123       TOFD(NGKINE)=0.
 124       GPOS(1,1) = VECT(1)
 125       GPOS(2,1) = VECT(2)
 126       GPOS(3,1) = VECT(3)
 127 C
 128 C             Rotate electron and scattered photon into GEANT system
 129 C
 130       CALL GFANG(VECT(4),COSTH,SINTH,COSPH,SINPH,ROTATE)
 131       IF(ROTATE) THEN
 132          CALL GDROT(PGAM(1),COSTH,SINTH,COSPH,SINPH)
 133          CALL GDROT(GKIN,COSTH,SINTH,COSPH,SINPH)
 134       ENDIF
 135 C
 136 C             Correct photon for energy lost and scattered angle
 137 C
 138       DO 60 I=1,3
 139    60 VECT(I+3) = PGAM(I)/EGAM2
 140       VECT(7) = EGAM2
 141       GETOT = EGAM2
 142       GEKIN = EGAM2
 143       CALL GEKBIN
 144 C
 145 C              Stop electron ?
 146 C
 147       IF((ICOMP.NE.1).OR.(TEL.LE.CUTELE)) THEN
 148          NGKINE = 0
 149          DESTEP = DESTEP + TEL
 150       ENDIF
 151 C
 152 C             Update probabilities
 153 C
 154       CALL GRNDM(RNDM,1)
 155       ZINTCO=-LOG(RNDM(1))
 156 C
 157       END