This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / gphys / glando.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:21:25  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/02 29/03/94  15.41.21  by  S.Giani
11 *-- Author :
12       SUBROUTINE GLANDO(IMODE,STEP,Z,A,RHO,P,E,XMASS,DE,IFLAG)
13 C.
14 C.    ******************************************************************
15 C.    *                                                                *
16 C.    *                                                                *
17 C.    *  GLANDO                                                        *
18 C.    *  ------                                                        *
19 C.    *                                                                *
20 C.    *  Energy straggling using Gaussian, Landau & Vavilov theories.  *
21 C.    *                                                                *
22 C.    *  Input                                                         *
23 C.    *  -----                                                         *
24 C.    *  IMODE  =  3   Landau(RANLAN)  sampling                        *
25 C.    *         =  4   Landau(GENLAN)  sampling                        *
26 C.    *         =  5   Vavilov(DINVAV) sampling                        *
27 C.    *         =  6   Gaussian sampling                               *
28 C.    *         =  2   Automatic selection of relevant distribution    *
29 C.    *                                                                *
30 C.    *  STEP   =  current step-length (cm)                            *
31 C.    *                                                                *
32 C.    *  Output                                                        *
33 C.    *  ------                                                        *
34 C.    *  DE     =  DE/DX - <DE/DX>     (GeV)                           *
35 C.    *                                                                *
36 C.    *  IFLAG  =  3   Landau(RANLAN) sampling used                    *
37 C.    *         =  4   Landau(GENLAN) sampling used                    *
38 C.    *         =  5   Vavilov(GVAVIV) sampling used                   *
39 C.    *         =  6   Gaussian sampling used                          *
40 C.    *                                                                *
41 C.    *  Warning                                                       *
42 C.    *  -------                                                       *
43 C.    *  Only Landau sampling should be used since this has been well  *
44 C.    *  tested whereas both Vavilov and Gaussian sampling are being   *
45 C.    *  developed.                                                    *
46 C.    *                                                                *
47 C.    *  Author      : G.N. Patrick                                    *
48 C.    *  Date        : 03.05.1985                                      *
49 C.    *  Last update : 09.09.1985                                      *
50 C.    *                                                                *
51 C.    ******************************************************************
52 C.
53       PARAMETER (EULER=0.577215)
54       PARAMETER (P1=.60715,P2=.67794,P3=.52382E-1,P4=.94753,
55      +           P5=.74442,P6=1.1934)
56 #include "geant321/gconsp.inc"
57       DIMENSION RNDM(2)
58       FLAND(X) = P1+P6*X+(P2+P3*X)*EXP(P4*X+P5)
59 C.
60 C.    ------------------------------------------------------------------
61 C.
62       IF(STEP.LT.1.E-7)THEN
63          DE=0.
64          IFLAG=0
65          RETURN
66       ENDIF
67 C
68 C     Calculate xi factor (KeV).
69 C
70       BETA   = P/E
71       GAMMA  = E/XMASS
72       XI     = (153.5*Z*STEP*RHO)/(A*BETA*BETA)
73 C
74 C     Maximum energy transfer to atomic electron (KeV).
75 C
76       ETA    = BETA*GAMMA
77       ETASQ  = ETA*ETA
78       RATIO  = EMASS/XMASS
79       F1     = 2.*EMASS*ETASQ
80       F2     = 1.+2.*RATIO*GAMMA+RATIO*RATIO
81       EMAX   = F1*1.E+6/F2
82 C
83 C     Calculate Kappa significance ratio.
84 C
85       CAPPA  = XI/EMAX
86 C
87 C     Choose correct function if IMODE set to automatic selection.
88 C
89       IMODEI = IMODE
90       IF (IMODEI.EQ.2)    THEN
91          IF (CAPPA.LT.0.01)                              IMODEI = 4
92 *****    IF (CAPPA.GE.0.01.AND.CAPPA.LE.10.)             IMODEI = 3
93          IF (CAPPA.GE.0.01.AND.CAPPA.LE.10.)             IMODEI = 5
94          IF (CAPPA.GT.10.)                               IMODEI = 6
95       ENDIF
96 C
97 C     +---------------------------------------------------------------+
98 C     I Sample lambda variable from Kolbig/Schorr Landau distribution I
99 C     +---------------------------------------------------------------+
100 C
101    10 CONTINUE
102       IF (IMODEI.EQ.3) THEN
103          XMEAN = -BETA**2-LOG(XI/EMAX)+EULER-1.
104          XLAMX = FLAND(XMEAN)
105    20    CALL GRNDM(RNDM,1)
106          IF( RNDM(1) .GT. 0.980 ) GO TO 20
107          XLAMB = GLANDR(RNDM(1))
108          IF(XLAMB.GT.XLAMX) GO TO 20
109          IFLAG  = 3
110 C
111 C            +----------------------------------------------+
112 C            I  Sample lambda variable from James & Hancock I
113 C            I  Landau distribution                         I
114 C            +----------------------------------------------+
115 C
116  
117       ELSEIF (IMODEI.EQ.4) THEN
118          XMEAN = -BETA**2-LOG(XI/EMAX)+EULER-1.
119          XLAMX = FLAND(XMEAN)
120   25     CALL GLANDG(XLAMB)
121          IF(XLAMB.GT.XLAMX) GO TO 25
122          IFLAG = 4
123 C
124 C            +---------------------------------------------------------+
125 C            I Sample lambda variable (Landau not Vavilov) from        I
126 C            I Rotondi&Montagna&Kolbig Vavilov distribution            I
127 C            +---------------------------------------------------------+
128 C
129       ELSEIF (IMODEI.EQ.5) THEN
130 C
131 C            Keep within computation bounds.
132 C
133          RKA = CAPPA
134          BE2 = BETA*BETA
135          CALL GRNDM(RNDM,1)
136          XLAMB = GVAVIV(RKA,BE2,RNDM(1))
137          IFLAG = 5
138 C
139 C     +-----------------------------------+
140 C     I Sample from Gaussian distribution I
141 C     +-----------------------------------+
142 C
143       ELSEIF (IMODEI.EQ.6) THEN
144          SIGMA  = XI*EMAX*(1.-(BETA*BETA/2.))
145          SIGMA  = SQRT(SIGMA)
146    30    CALL GRNDM(RNDM,2)
147          IF(RNDM(1).LE.0.)GO TO 30
148          F1     = -2.*LOG(RNDM(1))
149          DEKEV  = SIGMA*SQRT(F1)*COS(2.*PI*RNDM(2))
150          IFLAG  = 6
151          GOTO 40
152       ENDIF
153 C
154 C     Calculate DE/DX - <DE/DX>
155 C
156       DEKEV  = XI*(XLAMB+BETA*BETA+LOG(XI/EMAX)-EULER+1.)
157 C
158    40 DE=DEKEV*1.E-6
159       END