GEANT321/gphys/gmults.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:28  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.22  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GMULTS
  13 C
  14 C     ******************************************************************
  15 C     *                                                                *
  16 C     *       Steering routine for the multiple scattering.            *
  17 C     *       select Moliere theory , Gaussian approximation           *
  18 C     *       or single Coulomb scattering depending of their range    *
  19 C     *       of validity.                                             *
  20 C     *                                                                *
  21 C     *    ==>Called by : GTELEC , GTHADR , GTMUON                     *
  22 C.    *       Author     M.Maire  *********                            *
  23 C     *                                                                *
  24 C     ******************************************************************
  25 *
  26       PARAMETER (THRMOL=50.,THRGAU=0.01)
  27 #include "geant321/gcbank.inc"
  28 #include "geant321/gcjloc.inc"
  29 #include "geant321/gctrak.inc"
  30 #include "geant321/gckine.inc"
  31 #include "geant321/gcphys.inc"
  32 #include "geant321/gcmate.inc"
  33 #include "geant321/gcmulo.inc"
  34       DIMENSION DIN(3)
  35 *
  36       BETA2 = (VECT(7)/GETOT)**2
  37       IF(BETA2.LE.0.) RETURN
  38       IF(IMULS.EQ.3) THEN
  39          CALL GMGAUS(BETA2,DIN)
  40 *
  41       ELSE
  42          CHARG2=CHARGE**2
  43 *
  44 *     Here we decide whether we have to recalculate OMCMOL or not.
  45 *     OMCMOL has been calculated in GMULOF according to a formula
  46 *     which contains the term (1 + 3.34*(Alpha*z*Z/Beta)**2) where
  47 *     z is the incident charge (CHARGE), setting Beta=1 and z=1.
  48 *     We do not recalculate OMCMOL if:
  49 *
  50 *              3.34*(Z*Alpha)**2*((z/Beta)**2-1) << 1
  51 *                   Z**2*(z**2-Beta**2)/Beta**2  << 3.34/Alpha**2 = 5500
  52 *                   Z**2*(z**2-Beta**2)  <= 50.*Beta**2
  53 *
  54 *     We further multiply the first term for the number of elements
  55 *     in the mixture because the approssimation is particularly bad
  56 *     for mixtures so we make the condition more restrictive.
  57 *     All this thanks to Gerry Lynch.
  58 *
  59          IF(Z**2*(CHARG2-BETA2)*Q(JMA+11).GT.THRMOL*BETA2) THEN
  60             NLM=Q(JMA+11)
  61 *
  62             IF(NLM.EQ.1)THEN
  63                CALL GMOLIO(A,Z,1.,1,DENS,BETA2,CHARG2,OMCMOL)
  64 *
  65             ELSE
  66                CALL GMOLIO(Q(JMIXT+1),Q(JMIXT+NLM+1),Q(JMIXT+2*NLM+1),
  67      +         NLM,DENS,BETA2,CHARG2,OMCMOL)
  68 *
  69             ENDIF
  70          ELSE
  71             JPROB = LQ(JMA-4)
  72             OMCMOL = Q(JPROB+21)*CHARG2
  73 *
  74          ENDIF
  75 *
  76          OMEGA = OMCMOL*STMULS/BETA2
  77 *
  78          IF (OMEGA.LE.20.) THEN
  79             CALL GMCOUL(OMEGA,DIN)
  80 *
  81          ELSE
  82             CALL GMOLIE(OMEGA,BETA2,DIN)
  83 *
  84          ENDIF
  85       ENDIF
  86 *
  87 * *** Computes rotation matrix around particle direction
  88 * *** Compute new direction cosines
  89 *
  90       VMM = SQRT(VECT(4)*VECT(4)+VECT(5)*VECT(5))
  91       IF (VMM.NE.0.) THEN
  92          PD1=VECT(4)/VMM
  93          PD2=VECT(5)/VMM
  94          V4= PD1*VECT(6)*DIN(1) -PD2*DIN(2) +VECT(4)*DIN(3)
  95          V5= PD2*VECT(6)*DIN(1) +PD1*DIN(2) +VECT(5)*DIN(3)
  96          V6= -VMM*DIN(1) +VECT(6)*DIN(3)
  97       ELSE
  98          V4= DIN(1)
  99          V5= DIN(2)
 100          V6= DIN(3)*SIGN(1.,VECT(6))
 101       ENDIF
 102 *
 103 * *** Renormalize direction cosines
 104 *
 105       VP = 1./SQRT(V4*V4+V5*V5+V6*V6)
 106       VECT(4) = V4*VP
 107       VECT(5) = V5*VP
 108       VECT(6) = V6*VP
 109 *
 110       END