GEANT321/gtrak/ggckov.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:41  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.23  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GGCKOV
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *   This routine is called for each tracking step of a charged   *
  17 C.    *   particle in a radiator. A Poisson-distributed number of      *
  18 C.    *   photons is generated according to the Cherenkov formula,     *
  19 C.    *   distributed evenly along the track segment and uniformly     *
  20 C.    *   azimuth w.r.t. the particle direction. The parameters are    *
  21 C.    *   then transformed into the Master Reference System, and they  *
  22 C.    *   are added to the particle stack.                             *
  23 C.    *                                                                *
  24 C.    *   ==>Called by : GTMUON, GTHADR, GTELEC                        *
  25 C.    *      Authors     R.Jones, F.Carminati ********                 *
  26 C.    *                                                                *
  27 C.    ******************************************************************
  28 C.
  29 #include "geant321/gcbank.inc"
  30 #include "geant321/gcjloc.inc"
  31 #include "geant321/gctmed.inc"
  32 #include "geant321/gcunit.inc"
  33 #include "geant321/gctrak.inc"
  34 #include "geant321/gckine.inc"
  35 #include "geant321/gcking.inc"
  36 #include "geant321/gconsp.inc"
  37 *
  38       REAL RPHOT(3)
  39       LOGICAL ROTATE
  40       PARAMETER (RFACT=369.81E9)
  41 *
  42 *   ------------------------------------------------------------------
  43 *
  44 * *** See whether we generate at least one photon
  45 *
  46 C     THRIND = GETOT/VECT(7)
  47 C     IF(Q(JINDEX+NPCKOV).LT.THRIND) THEN
  48 C        GO TO 999
  49 C     ELSEIF(Q(JINDEX+1).GE.THRIND) THEN
  50 C        PMIN  = Q(JTCKOV+2)
  51 C        DP    = Q(JTCKOV+NPCKOV+1)-PMIN
  52 C        GE    = Q(JCURIN+NPCKOV)
  53 C        JMIN  = 1
  54 C     ELSE
  55 C        JMIN = 1
  56 C        JMAX = NPCKOV
  57 C  10    JMED = (JMIN+JMAX)/2
  58 C        IF(Q(JINDEX+JMED).LT.THRIND) THEN
  59 C           JMIN = JMED
  60 C        ELSE
  61 C           JMAX = JMED
  62 C        ENDIF
  63 C        IF(JMAX-JMIN.GT.1) GO TO 10
  64 C        RATIO =
  65 C    +   (THRIND-Q(JINDEX+JMIN))/(Q(JINDEX+JMIN+1)-Q(JINDEX+JMIN))
  66 C        RATI1 = 1.-RATIO
  67 C        PMIN  = Q(JTCKOV+JMIN+1)*RATI1+Q(JTCKOV+JMIN+2)*RATIO
  68 C        DP    = Q(JTCKOV+NPCKOV+1)-PMIN
  69 C        GEMIN = Q(JCURIN+JMIN)*RATI1+Q(JCURIN+JMIN+1)*RATIO
  70 C        GE    = Q(JCURIN+NPCKOV)-GEMIN
  71 C     ENDIF
  72 C     DNDL = RFACT*(CHARGE**2)*(DP-GE*THRIND**2)
  73       IF(ITRTYP.NE.4.AND.ITRTYP.NE.8) CALL GNCKOV
  74       CALL GPOISS(DNDL*STEP,NGPHOT,1)
  75       IF(NGPHOT.EQ.0) THEN
  76          GO TO 999
  77       ELSEIF(NGPHOT.GT.MXPHOT) THEN
  78          WRITE(CHMAIL,10000) NGPHOT-MXPHOT
  79 10000   FORMAT(' **** GGCKOV Overflow in the photon stack, ',I10,
  80      +         ' photons are lost')
  81          CALL GMAIL(0,0)
  82          NGPHOT = MXPHOT
  83       ENDIF
  84 *
  85 * ***  Set up rotation to Particle frame
  86 *
  87       CALL GFANG(VECT(4),COSTH,SINTH,COSPH,SINPH,ROTATE)
  88 *
  89 * ***  Distribute the photons in origin, direction, momentum
  90       COSMX  = THRIND/Q(JINDEX+NPCKOV)
  91       SINMX2 = (1.-COSMX)*(1.+COSMX)
  92       DO 40 J=1,NGPHOT
  93          CALL GRNDM(RPHOT, 1)
  94          IF(IGNEXT.NE.0) THEN
  95             DS=(STEP-PREC)*RPHOT(1)+PREC
  96          ELSE
  97             DS = STEP*RPHOT(1)
  98          ENDIF
  99          XPHOT(1,J) = VECT(1)-VECT(4)*DS
 100          XPHOT(2,J) = VECT(2)-VECT(5)*DS
 101          XPHOT(3,J) = VECT(3)-VECT(6)*DS
 102          XPHOT(11,J)= TOFG+(STEP-DS)*GETOT/(VECT(7)*CLIGHT)
 103 *
 104 * *** Sample the momentum of the photon
 105    20    CALL GRNDM(RPHOT, 3)
 106          PPHOT=PMIN+RPHOT(1)*DP
 107 *
 108 * *** Find in which bin we are
 109          KMIN = JMIN
 110          KMAX = NPCKOV
 111    30    KMED = (KMIN+KMAX)/2
 112          IF(Q(JTCKOV+1+KMED).LT.PPHOT) THEN
 113             KMIN = KMED
 114          ELSE
 115             KMAX = KMED
 116          ENDIF
 117          IF(KMAX-KMIN.GT.1) GOTO 30
 118          RATIO = (PPHOT-Q(JTCKOV+1+KMIN))/
 119      +           (Q(JTCKOV+KMIN+2)-Q(JTCKOV+1+KMIN))
 120          RATI1  = (1.-RATIO)
 121 *
 122 * *** Find the density function value corresponding to the
 123 * *** momentum sampled
 124          RINDEX = Q(JINDEX+KMIN)*RATI1+Q(JINDEX+KMIN+1)*RATIO
 125          COST   = THRIND/RINDEX
 126          SINT2  = (1.-COST)*(1.+COST)
 127 *
 128 * *** Perform hit-and-miss
 129          IF(RPHOT(2)*SINMX2.GT.SINT2) GO TO 20
 130          SINT = SQRT(SINT2)
 131          PHI  = TWOPI*RPHOT(3)
 132          SINP = SIN(PHI)
 133          COSP = COS(PHI)
 134          XPHOT(4,J) = SINT*COSP
 135          XPHOT(5,J) = SINT*SINP
 136          XPHOT(6,J) = COST
 137          XPHOT(7,J) = PPHOT
 138          XPHOT(8,J) = COST*COSP
 139          XPHOT(9,J) = COST*SINP
 140          XPHOT(10,J) = -SINT
 141 *
 142          IF(ROTATE) THEN
 143             CALL GDROT(XPHOT(8,J),COSTH,SINTH,COSPH,SINPH)
 144             CALL GDROT(XPHOT(4,J),COSTH,SINTH,COSPH,SINPH)
 145          ENDIF
 146    40 CONTINUE
 147   999 END