This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / 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