GEANT321/gphys/ganni.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.2  1996/02/27 10:08:11  ravndal
   6 * Precision problem in cos(theta) solved
   7 *
   8 * Revision 1.1.1.1  1995/10/24 10:21:21  cernlib
   9 * Geant
  10 *
  11 *
  12 #include "geant321/pilot.h"
  13 *CMZ :  3.21/02 29/03/94  15.41.21  by  S.Giani
  14 *-- Author :
  15       SUBROUTINE GANNI
  16 C.
  17 C.    ******************************************************************
  18 C.    *                                                                *
  19 C.    *  Generates positron annihilation                               *
  20 C.    *                                                                *
  21 C.    *    ==>Called by : GTELEC                                       *
  22 C.    *       Author    L.Urban *********                              *
  23 C.    *       10/06/93: modified by Georges Azuelos (Vancouver)        *
  24 C     *                 to include 1-quantum annihilation              *
  25 C.    *                                                                *
  26 C.    ******************************************************************
  27 C.
  28 #include "geant321/gcphys.inc"
  29 #include "geant321/gctrak.inc"
  30 #include "geant321/gccuts.inc"
  31 #include "geant321/gcking.inc"
  32 #include "geant321/gconsp.inc"
  33 #include "geant321/gcbank.inc"
  34 #include "geant321/gcmulo.inc"
  35 #include "geant321/gcjloc.inc"
  36 #include "geant321/gcmate.inc"
  37       DIMENSION PGAM(3),RNDM(2)
  38       LOGICAL ROTATE
  39       PARAMETER (ALFA=7.29735E-3)
  40 C.
  41 C.    ------------------------------------------------------------------
  42 C.
  43       KCASE = NAMEC(11)
  44       IF((IANNI.NE.1).OR.((GETOT+EMASS).LE.CUTGAM)) THEN
  45          ISTOP = 2
  46          DESTEP = DESTEP + GETOT + EMASS
  47          GEKIN = 0.
  48          GETOT = EMASS
  49          VECT(7)= 0.
  50          GO TO 999
  51       ENDIF
  52 C
  53       XE=GETOT
  54       GAM=XE/EMASS
  55       GAM2=GAM**2
  56       GAM1=MAX(GAM2-1.,0.)
  57       GAMP1=GAM+1.
  58       C=SQRT(GAM1)
  59 *
  60       SIG=(GAM2+4.*GAM+1.)*LOG(GAM+C)/GAM1-(GAM+3.)/C
  61       SIG=0.5*Q(JPROB+17)*SIG/GAMP1
  62 *
  63       BIND=0.5*(Z*ALFA)**2*EMASS
  64       E1Q=XE+EMASS-BIND
  65       IF(E1Q .GT. 0.)THEN
  66          GVE=VECT(7)/EMASS
  67          SIG1=GAM2+2.*(GAM+2.)/3.-(GAM+2.)/GVE*LOG(GAM+GVE)
  68          SIG1=2.*Q(JPROB+18)*SIG1/(GVE*GAMP1**2)
  69       ELSE
  70          SIG1=0.
  71       END IF
  72 *
  73       SIG=SIG+SIG1
  74       CALL GRNDM(RNDM,1)
  75 C
  76       IF(RNDM(1).GE.SIG1/SIG)THEN
  77          GAMP12=GAMP1**2
  78          P=VECT(7)
  79          E0=1./(GAMP1+C)
  80 C
  81    10    CALL GRNDM(RNDM,2)
  82          E=E0*((1.-E0)/E0)**RNDM(1)
  83 C
  84          SCREJ=(GAMP12+2.*GAMP1-2.-GAMP12*E-1./E)/(GAMP12-2.)
  85          IF(RNDM(2).GT.SCREJ) GOTO 10
  86 C
  87          EPHOT1=(XE+EMASS)*E
  88 C
  89          COSTH=(GEKIN+EMASS*(2.*E-1.)/E)/P
  90 C
  91 C restrict COSTH to [-1.,+1.]
  92 C
  93          COSTH = MIN( 1. , MAX( -1. , COSTH ) )
  94          SINTH=SQRT((1.-COSTH)*(1.+COSTH))
  95          CALL GRNDM(RNDM,1)
  96          PHI = TWOPI * RNDM(1)
  97          COSPHI = COS(PHI)
  98          SINPHI = SIN(PHI)
  99 C
 100          PGAM(1) = EPHOT1* SINTH * COSPHI
 101          PGAM(2) = EPHOT1* SINTH * SINPHI
 102          PGAM(3) = EPHOT1* COSTH
 103 C
 104 C             Rotate tracks into GEANT system and store.
 105 C
 106          CALL GFANG(VECT(4),COSAL,SINAL,COSBT,SINBT,ROTATE)
 107 C
 108 C            Polar co-ordinates to momentum components.
 109 C
 110          NGGAMM = 0
 111          IF(EPHOT1.GT.CUTGAM) THEN
 112             NGGAMM = NGGAMM + 1
 113             NGKINE = NGKINE + 1
 114             GKIN(1,NGKINE) = PGAM(1)
 115             GKIN(2,NGKINE) = PGAM(2)
 116             GKIN(3,NGKINE) = PGAM(3)
 117             GKIN(4,NGKINE) = EPHOT1
 118             GKIN(5,NGKINE) = 1
 119             TOFD(NGKINE)=0.
 120             GPOS(1,NGKINE) = VECT(1)
 121             GPOS(2,NGKINE) = VECT(2)
 122             GPOS(3,NGKINE) = VECT(3)
 123             IF(ROTATE)
 124      +      CALL GDROT(GKIN(1,NGKINE),COSAL,SINAL,COSBT,SINBT)
 125          ELSE
 126             DESTEP = DESTEP + EPHOT1
 127          ENDIF
 128 C
 129 C             Momentum vector of second photon.
 130 C
 131          EPHOT2 = GETOT + EMASS - EPHOT1
 132          IF(EPHOT2.GT.CUTGAM) THEN
 133             NGGAMM = NGGAMM + 1
 134             NGKINE = NGKINE + 1
 135             GKIN(1,NGKINE) = - PGAM(1)
 136             GKIN(2,NGKINE) = - PGAM(2)
 137             GKIN(3,NGKINE) = P - PGAM(3)
 138             GKIN(4,NGKINE) = EPHOT2
 139             GKIN(5,NGKINE) = 1
 140             TOFD(NGKINE)=0.
 141             GPOS(1,NGKINE) = VECT(1)
 142             GPOS(2,NGKINE) = VECT(2)
 143             GPOS(3,NGKINE) = VECT(3)
 144             IF(ROTATE)
 145      +      CALL GDROT(GKIN(1,NGKINE),COSAL,SINAL,COSBT,SINBT)
 146          ELSE
 147             DESTEP = DESTEP + EPHOT2
 148          ENDIF
 149       ELSE
 150 C 1-quantum annihilation
 151          P=VECT(7)
 152          EPHOT=E1Q
 153 C Assume photon collinear with positron
 154          PGAM(1) = 0.
 155          PGAM(2) = 0.
 156          PGAM(3) = EPHOT
 157 C
 158 C             Rotate tracks into GEANT system and store.
 159 C
 160          CALL GFANG(VECT(4),COSAL,SINAL,COSBT,SINBT,ROTATE)
 161 C
 162 C            Polar co-ordinates to momentum components.
 163 C
 164          NGGAMM = 0
 165          IF(EPHOT.GT.CUTGAM) THEN
 166             NGGAMM = NGGAMM + 1
 167             NGKINE = NGKINE + 1
 168             GKIN(1,NGKINE) = PGAM(1)
 169             GKIN(2,NGKINE) = PGAM(2)
 170             GKIN(3,NGKINE) = PGAM(3)
 171             GKIN(4,NGKINE) = EPHOT
 172             GKIN(5,NGKINE) = 1
 173             TOFD(NGKINE)=0.
 174             GPOS(1,NGKINE) = VECT(1)
 175             GPOS(2,NGKINE) = VECT(2)
 176             GPOS(3,NGKINE) = VECT(3)
 177             IF(ROTATE)
 178      +      CALL GDROT(GKIN(1,NGKINE),COSAL,SINAL,COSBT,SINBT)
 179          ELSE
 180             DESTEP = DESTEP + EPHOT
 181          ENDIF
 182       END IF
 183 C
 184       IF(NGGAMM.GT.0) THEN
 185          ISTOP = 1
 186       ELSE
 187          ISTOP = 2
 188       ENDIF
 189 C
 190  999  END