GEANT321/gphys/gdray.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 02/07/94  18.28.12  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GDRAY
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *  Generates Delta rays                                          *
  17 C.    *                                                                *
  18 C.    *    ==>Called by : GTELEC,GTHADR,GTMUON                         *
  19 C.    *       Authors    D.Ward, L.Urban  ********                     *
  20 C.    *                                                                *
  21 C.    ******************************************************************
  22 C.
  23 #include "geant321/gcphys.inc"
  24 #include "geant321/gctrak.inc"
  25 #include "geant321/gconsp.inc"
  26 #include "geant321/gckine.inc"
  27 #include "geant321/gcking.inc"
  28 #include "geant321/gccuts.inc"
  29       DIMENSION PELS(3)
  30       DIMENSION RNDM(2)
  31       LOGICAL ROTATE
  32 C.
  33 C.    ------------------------------------------------------------------
  34 C.
  35       P=VECT(7)
  36       XE=GETOT
  37       TE=GEKIN
  38       GAM=XE/AMASS
  39       GAM2=GAM*GAM
  40       T=GAM-1.
  41       X=DCUTE/(T*EMASS)
  42 C
  43       KCASE = NAMEC(10)
  44       IF(IPART.EQ.3)   THEN
  45 C
  46 C======>       Moller scattering
  47 C
  48         IF(X.GE.0.5) GO TO 90
  49         CC=1.-2.*X
  50 C
  51   10    CALL GRNDM(RNDM,2)
  52         E=X/(1.-CC*RNDM(1))
  53 C
  54         B1=4./(9.*GAM2-10.*GAM+5.)
  55         B2=T*T*B1
  56         B3=(2.*GAM2+2.*GAM-1.)*B1
  57         E1=1.-E
  58 C
  59         SCREJ=B2*E*E-B3*E/E1+B1*GAM2/(E1*E1)
  60 C
  61         IF(RNDM(2).GT.SCREJ) GOTO 10
  62 C
  63       ELSEIF(IPART.EQ.2)THEN
  64 C
  65 C======>       Bhabha scattering
  66 C
  67         IF(X.GE.1.) GO TO 90
  68         X1=1.-X
  69   20    CALL GRNDM(RNDM,2)
  70         E=X/(1.-X1*RNDM(1))
  71 C
  72         Y=1./(GAM+1.)
  73         Y2=Y*Y
  74         C=1.-2.*Y
  75         B1=2.-Y2
  76         B2=C*(3.+Y2)
  77         C2=C*C
  78         B4=C2*C
  79         B3=C2+B4
  80         B0=GAM2/(GAM2-1.)
  81 C
  82         SCREJ=(((B4*E-B3)*E+B2)*E-B1)*E+B0
  83         SCREJ=SCREJ/((((B4*X-B3)*X+B2)*X-B1)*X+B0)
  84         IF(RNDM(2).GT.SCREJ) GOTO 20
  85 C
  86       ELSE
  87 C
  88 C======>     Heavy particle.
  89 C
  90         TMAX=2.*EMASS*(GAM2-1.)/
  91      +  (1.+2.*GAM*EMASS/AMASS+(EMASS/AMASS)**2)
  92         IF(TMAX.LE.DCUTM)  GOTO 90
  93   40    CALL GRNDM(RNDM,2)
  94         E=1./DCUTM+RNDM(1)*(1./TMAX-1./DCUTM)
  95         E=1./E
  96         BET2=1.-1./GAM2
  97         SCREJ=1.-BET2*(E/TMAX)
  98 C ---         extra term for spin 1/2 parent.
  99         IF(AMASS.GT.0.9 .OR. AMASS.LT.0.12)
 100      +  SCREJ=SCREJ+0.5*(E/GETOT)**2
 101         IF(RNDM(2).GT.SCREJ) GO TO 40
 102         E=E/(T*EMASS)
 103 C
 104       ENDIF
 105 C
 106       EEL=(T*E+1.)*EMASS
 107       TEL=EEL-EMASS
 108       PEL=SQRT(ABS((EEL+EMASS)*TEL))
 109       COSTH=(XE*EEL+EMASS*(TEL-XE))/(P*PEL)
 110       IF(COSTH.GE.1.) THEN
 111          COSTH=1.
 112          SINTH=0.
 113       ELSEIF(COSTH.LE.-1.) THEN
 114          COSTH=-1.
 115          SINTH=0.
 116       ELSE
 117          SINTH=SQRT((1.+COSTH)*(1.-COSTH))
 118       ENDIF
 119       CALL GRNDM(RNDM,1)
 120       PHI    = TWOPI*RNDM(1)
 121       COSPHI = COS(PHI)
 122       SINPHI = SIN(PHI)
 123 C
 124 C             Polar co-ordinates to momentum components.
 125 C
 126       NGKINE = 1
 127       GKIN(1,1)=PEL*SINTH*COSPHI
 128       GKIN(2,1)=PEL*SINTH*SINPHI
 129       GKIN(3,1)=PEL*COSTH
 130       GKIN(4,1)=EEL
 131       GKIN(5,1)=3
 132       TOFD(NGKINE)=0.
 133       GPOS(1,NGKINE) = VECT(1)
 134       GPOS(2,NGKINE) = VECT(2)
 135       GPOS(3,NGKINE) = VECT(3)
 136 C
 137       PELS(1)=-GKIN(1,1)
 138       PELS(2)=-GKIN(2,1)
 139       PELS(3)=P-GKIN(3,1)
 140 C
 141       CALL GFANG(VECT(4),COSTH,SINTH,COSPH,SINPH,ROTATE)
 142       IF(ROTATE) THEN
 143          CALL GDROT(PELS(1),COSTH,SINTH,COSPH,SINPH)
 144          CALL GDROT(GKIN,COSTH,SINTH,COSPH,SINPH)
 145       ENDIF
 146 C
 147 C             Correct track vector for lost energy and scattered angles
 148 C
 149       TELS=TE-TEL
 150       EELS=TELS+AMASS
 151       PEELS=SQRT(ABS((EELS+AMASS)*TELS))
 152       IF(PEELS.GT.0.)THEN
 153          DO 55 I=1,3
 154             VECT(I+3) = PELS(I)/PEELS
 155   55     CONTINUE
 156       ENDIF
 157       VECT(7) = PEELS
 158       GEKIN=TELS
 159       GETOT=EELS
 160       CALL GEKBIN
 161       IF((IDRAY.NE.1).OR.(TEL.LE.CUTELE)) THEN
 162           NGKINE = 0
 163           DESTEP = DESTEP + TEL
 164       ENDIF
 165 C
 166 C             Update probabilities
 167 C
 168   90  CALL GRNDM(RNDM,1)
 169       ZINTDR=-LOG(RNDM(1))
 170       SLDRAY=SLENG
 171       STEPDR=BIG
 172 C
 173       END