This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / gphys / gvaviv.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:21:35  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       FUNCTION GVAVIV(RKAPPA,BETA2,RAN)
13 C.
14 C.    ******************************************************************
15 C.    *                                                                *
16 C.    *                                                                *
17 C.    *  GVAVIV                                                        *
18 C.    *  ------                                                        *
19 C.    *                                                                *
20 C.    *  Initializes the parameters for a Vavilov distribution         *
21 C.    *                                                                *
22 C.    *  BETA2 is the particle velocity squared                        *
23 C.    *  RKAPPA the usual straggling parameter                         *
24 C.    *                                                                *
25 C.    *  This routine has been extracted from a proposed CERNLIB       *
26 C.    *  set of routines (VAVCOE,VAVFSM).                              *
27 C.    *  The authors of these routines have submitted an article       *
28 C.    *  to NIM:                                                       *
29 C.    *                                                                *
30 C.    *    Alberto Rotondi, Paolo Montagna                             *
31 C.    *    Fast calculation of Vavilov distribution NIMB 1990          *
32 C.    *                                                                *
33 C.    *       Author    K.S.Koelbig     ********                       *
34 C.    *                                                                *
35 C.    *   ==> Called by : GLANDO                                       *
36 C.    *                                                                *
37 C.    ******************************************************************
38 C.
39       DIMENSION AC(0:12),HC(0:8),H(9)
40  
41       PARAMETER
42      1(BKMNX1 = 0.02, BKMNY1 = 0.05, BKMNX2 = 0.12, BKMNY2 = 0.05,
43      2 BKMNX3 = 0.22, BKMNY3 = 0.05, BKMXX1 = 0.1 , BKMXY1 = 1   ,
44      3 BKMXX2 = 0.2 , BKMXY2 = 1   , BKMXX3 = 0.3 , BKMXY3 = 1   )
45       PARAMETER
46      1(FBKX1 = 2/(BKMXX1-BKMNX1), FBKX2 = 2/(BKMXX2-BKMNX2),
47      2 FBKX3 = 2/(BKMXX3-BKMNX3), FBKY1 = 2/(BKMXY1-BKMNY1),
48      3 FBKY2 = 2/(BKMXY2-BKMNY2), FBKY3 = 2/(BKMXY3-BKMNY3))
49  
50       DIMENSION EDGEC(2:7),FNINV(5),DRK(5),DSIGM(5),ALFA(2:5)
51       DIMENSION U1(13),U2(13),U3(13),U4(12),U5(13),U6(13),U7( 8),U8(13)
52       DIMENSION V1(12),V2(12),V3(13),V4(12),V5(13),V6(13),V7(11),V8(11)
53       DIMENSION W1(13),W2(11),W3(13),W4(13),W5(13),W6(13),       W8( 8)
54  
55       DATA FNINV /1, 0.5, 0.33333333, 0.25, 0.2/
56  
57       DATA (EDGEC(J),J=2,7)
58      1/ 0.16666667E+0, 0.41666667E-1, 0.83333333E-2,
59      2  0.13888889E-1, 0.69444444E-2, 0.77160493E-3/
60  
61       DATA (U1(K),K=1,13)
62      1/ 0.25850868E+0,  0.32477982E-1, -0.59020496E-2,
63      2  0.            , 0.24880692E-1,  0.47404356E-2,
64      3 -0.74445130E-3,  0.73225731E-2,  0.           ,
65      4  0.11668284E-2,  0.           , -0.15727318E-2,-0.11210142E-2/
66  
67       DATA (U2(K),K=1,13)
68      1/ 0.43142611E+0,  0.40797543E-1, -0.91490215E-2,
69      2  0.           ,  0.42127077E-1,  0.73167928E-2,
70      3 -0.14026047E-2,  0.16195241E-1,  0.24714789E-2,
71      4  0.20751278E-2,  0.           , -0.25141668E-2,-0.14064022E-2/
72  
73       DATA (U3(K),K=1,13)
74      1/ 0.25225955E+0,  0.64820468E-1, -0.23615759E-1,
75      2  0.           ,  0.23834176E-1,  0.21624675E-2,
76      3 -0.26865597E-2, -0.54891384E-2,  0.39800522E-2,
77      4  0.48447456E-2, -0.89439554E-2, -0.62756944E-2,-0.24655436E-2/
78  
79       DATA (U4(K),K=1,12)
80      1/ 0.12593231E+1, -0.20374501E+0,  0.95055662E-1,
81      2 -0.20771531E-1, -0.46865180E-1, -0.77222986E-2,
82      3  0.32241039E-2,  0.89882920E-2, -0.67167236E-2,
83      4 -0.13049241E-1,  0.18786468E-1,  0.14484097E-1/
84  
85       DATA (U5(K),K=1,13)
86      1/-0.24864376E-1, -0.10368495E-2,  0.14330117E-2,
87      2  0.20052730E-3,  0.18751903E-2,  0.12668869E-2,
88      3  0.48736023E-3,  0.34850854E-2,  0.           ,
89      4 -0.36597173E-3,  0.19372124E-2,  0.70761825E-3, 0.46898375E-3/
90  
91       DATA (U6(K),K=1,13)
92      1/ 0.35855696E-1, -0.27542114E-1,  0.12631023E-1,
93      2 -0.30188807E-2, -0.84479939E-3,  0.           ,
94      3  0.45675843E-3, -0.69836141E-2,  0.39876546E-2,
95      4 -0.36055679E-2,  0.           ,  0.15298434E-2, 0.19247256E-2/
96  
97       DATA (U7(K),K=1,8)
98      1/ 0.10234691E+2, -0.35619655E+1,  0.69387764E+0,
99      2 -0.14047599E+0, -0.19952390E+1, -0.45679694E+0,
100      3  0.           ,  0.50505298E+0/
101  
102       DATA (U8(K),K=1,13)
103      1/ 0.21487518E+2, -0.11825253E+2,  0.43133087E+1,
104      2 -0.14500543E+1, -0.34343169E+1, -0.11063164E+1,
105      3 -0.21000819E+0,  0.17891643E+1, -0.89601916E+0,
106      4  0.39120793E+0,  0.73410606E+0,  0.           ,-0.32454506E+0/
107  
108       DATA (V1(K),K=1,12)
109      1/ 0.27827257E+0, -0.14227603E-2,  0.24848327E-2,
110      2  0.           ,  0.45091424E-1,  0.80559636E-2,
111      3 -0.38974523E-2,  0.           , -0.30634124E-2,
112      4  0.75633702E-3,  0.54730726E-2,  0.19792507E-2/
113  
114       DATA (V2(K),K=1,12)
115      1/ 0.41421789E+0, -0.30061649E-1,  0.52249697E-2,
116      2  0.           ,  0.12693873E+0,  0.22999801E-1,
117      3 -0.86792801E-2,  0.31875584E-1, -0.61757928E-2,
118      4  0.           ,  0.19716857E-1,  0.32596742E-2/
119  
120       DATA (V3(K),K=1,13)
121      1/ 0.20191056E+0, -0.46831422E-1,  0.96777473E-2,
122      2 -0.17995317E-2,  0.53921588E-1,  0.35068740E-2,
123      3 -0.12621494E-1, -0.54996531E-2, -0.90029985E-2,
124      4  0.34958743E-2,  0.18513506E-1,  0.68332334E-2,-0.12940502E-2/
125  
126       DATA (V4(K),K=1,12)
127      1/ 0.13206081E+1,  0.10036618E+0, -0.22015201E-1,
128      2  0.61667091E-2, -0.14986093E+0, -0.12720568E-1,
129      3  0.24972042E-1, -0.97751962E-2,  0.26087455E-1,
130      4 -0.11399062E-1, -0.48282515E-1, -0.98552378E-2/
131  
132       DATA (V5(K),K=1,13)
133      1/ 0.16435243E-1,  0.36051400E-1,  0.23036520E-2,
134      2 -0.61666343E-3, -0.10775802E-1,  0.51476061E-2,
135      3  0.56856517E-2, -0.13438433E-1,  0.           ,
136      4  0.           , -0.25421507E-2,  0.20169108E-2,-0.15144931E-2/
137  
138       DATA (V6(K),K=1,13)
139      1/ 0.33432405E-1,  0.60583916E-2, -0.23381379E-2,
140      2  0.83846081E-3, -0.13346861E-1, -0.17402116E-2,
141      3  0.21052496E-2,  0.15528195E-2,  0.21900670E-2,
142      4 -0.13202847E-2, -0.45124157E-2, -0.15629454E-2, 0.22499176E-3/
143  
144       DATA (V7(K),K=1,11)
145      1/ 0.54529572E+1, -0.90906096E+0,  0.86122438E-1,
146      2  0.           , -0.12218009E+1, -0.32324120E+0,
147      3 -0.27373591E-1,  0.12173464E+0,  0.           ,
148      4  0.           ,  0.40917471E-1/
149  
150       DATA (V8(K),K=1,11)
151      1/ 0.93841352E+1, -0.16276904E+1,  0.16571423E+0,
152      2  0.           , -0.18160479E+1, -0.50919193E+0,
153      3 -0.51384654E-1,  0.21413992E+0,  0.           ,
154      4  0.           ,  0.66596366E-1/
155  
156       DATA (W1(K),K=1,13)
157      1/ 0.29712951E+0,  0.97572934E-2,  0.           ,
158      2 -0.15291686E-2,  0.35707399E-1,  0.96221631E-2,
159      3 -0.18402821E-2, -0.49821585E-2,  0.18831112E-2,
160      4  0.43541673E-2,  0.20301312E-2, -0.18723311E-2,-0.73403108E-3/
161  
162       DATA (W2(K),K=1,11)
163      1/ 0.40882635E+0,  0.14474912E-1,  0.25023704E-2,
164      2 -0.37707379E-2,  0.18719727E+0,  0.56954987E-1,
165      3  0.           ,  0.23020158E-1,  0.50574313E-2,
166      4  0.94550140E-2,  0.19300232E-1/
167  
168       DATA (W3(K),K=1,13)
169      1/ 0.16861629E+0,  0.           ,  0.36317285E-2,
170      2 -0.43657818E-2,  0.30144338E-1,  0.13891826E-1,
171      3 -0.58030495E-2, -0.38717547E-2,  0.85359607E-2,
172      4  0.14507659E-1,  0.82387775E-2, -0.10116105E-1,-0.55135670E-2/
173  
174       DATA (W4(K),K=1,13)
175      1/ 0.13493891E+1, -0.26863185E-2, -0.35216040E-2,
176      2  0.24434909E-1, -0.83447911E-1, -0.48061360E-1,
177      3  0.76473951E-2,  0.24494430E-1, -0.16209200E-1,
178      4 -0.37768479E-1, -0.47890063E-1,  0.17778596E-1, 0.13179324E-1/
179  
180       DATA (W5(K),K=1,13)
181      1/ 0.10264945E+0,  0.32738857E-1,  0.           ,
182      2  0.43608779E-2, -0.43097757E-1, -0.22647176E-2,
183      3  0.94531290E-2, -0.12442571E-1, -0.32283517E-2,
184      4 -0.75640352E-2, -0.88293329E-2,  0.52537299E-2, 0.13340546E-2/
185  
186       DATA (W6(K),K=1,13)
187      1/ 0.29568177E-1, -0.16300060E-2, -0.21119745E-3,
188      2  0.23599053E-2, -0.48515387E-2, -0.40797531E-2,
189      3  0.40403265E-3,  0.18200105E-2, -0.14346306E-2,
190      4 -0.39165276E-2, -0.37432073E-2,  0.19950380E-2, 0.12222675E-2/
191  
192       DATA (W8(K),K=1,8)
193      1/ 0.66184645E+1, -0.73866379E+0,  0.44693973E-1,
194      2  0.           , -0.14540925E+1, -0.39529833E+0,
195      3 -0.44293243E-1,  0.88741049E-1/
196  
197       GVAVIV=0
198       IF(RKAPPA .LT. 0.01 .OR. RKAPPA .GT. 12) RETURN
199       IF(RKAPPA .GE. 0.29) THEN
200        ITYPE=1
201        NPT=100
202        WK=1/SQRT(RKAPPA)
203        AC(0)=(-0.032227*BETA2-0.074275)*RKAPPA+
204      1   (0.24533*BETA2+0.070152)*WK+(-0.55610*BETA2-3.1579)
205        AC(8)=(-0.013483*BETA2-0.048801)*RKAPPA+
206      1   (-1.6921*BETA2+8.3656)*WK+(-0.73275*BETA2-3.5226)
207        DRK(1)=WK**2
208        DSIGM(1)=SQRT(RKAPPA/(1-0.5*BETA2))
209        DO 1 J = 1,4
210        DRK(J+1)=DRK(1)*DRK(J)
211        DSIGM(J+1)=DSIGM(1)*DSIGM(J)
212     1  ALFA(J+1)=(FNINV(J)-BETA2*FNINV(J+1))*DRK(J)
213  
214        HC(0)=LOG(RKAPPA)+BETA2+0.42278434
215        HC(1)=DSIGM(1)
216        HC(2)=ALFA(3)*DSIGM(3)
217        HC(3)=(3*ALFA(2)**2+ALFA(4))*DSIGM(4)-3
218        HC(4)=(10*ALFA(2)*ALFA(3)+ALFA(5))*DSIGM(5)-10*HC(2)
219        HC(5)=HC(2)**2
220        HC(6)=HC(2)*HC(3)
221        HC(7)=HC(2)*HC(5)
222        DO 2 J = 2,7
223     2  HC(J)=EDGEC(J)*HC(J)
224        HC(8)=0.39894228*HC(1)
225       ELSEIF(RKAPPA .GE. 0.22) THEN
226        ITYPE=2
227        NPT=150
228        X=1+(RKAPPA-BKMXX3)*FBKX3
229        Y=1+(SQRT(BETA2)-BKMXY3)*FBKY3
230        XX=2*X
231        YY=2*Y
232        X2=XX*X-1
233        X3=XX*X2-X
234        Y2=YY*Y-1
235        Y3=YY*Y2-Y
236        XY=X*Y
237        P2=X2*Y
238        P3=X3*Y
239        Q2=Y2*X
240        Q3=Y3*X
241        PQ=X2*Y2
242        AC(1)=W1(1)+W1(2)*X+W1(4)*X3+W1(5)*Y+W1(6)*Y2+W1(7)*Y3+
243      1  W1(8)*XY+W1(9)*P2+W1(10)*P3+W1(11)*Q2+W1(12)*Q3+W1(13)*PQ
244        AC(2)=W2(1)+W2(2)*X+W2(3)*X2+W2(4)*X3+W2(5)*Y+W2(6)*Y2+
245      1  W2(8)*XY+W2(9)*P2+W2(10)*P3+W2(11)*Q2
246        AC(3)=W3(1)+W3(3)*X2+W3(4)*X3+W3(5)*Y+W3(6)*Y2+W3(7)*Y3+
247      1  W3(8)*XY+W3(9)*P2+W3(10)*P3+W3(11)*Q2+W3(12)*Q3+W3(13)*PQ
248        AC(4)=W4(1)+W4(2)*X+W4(3)*X2+W4(4)*X3+W4(5)*Y+W4(6)*Y2+W4(7)*Y3+
249      1  W4(8)*XY+W4(9)*P2+W4(10)*P3+W4(11)*Q2+W4(12)*Q3+W4(13)*PQ
250        AC(5)=W5(1)+W5(2)*X+W5(4)*X3+W5(5)*Y+W5(6)*Y2+W5(7)*Y3+
251      1  W5(8)*XY+W5(9)*P2+W5(10)*P3+W5(11)*Q2+W5(12)*Q3+W5(13)*PQ
252        AC(6)=W6(1)+W6(2)*X+W6(3)*X2+W6(4)*X3+W6(5)*Y+W6(6)*Y2+W6(7)*Y3+
253      1  W6(8)*XY+W6(9)*P2+W6(10)*P3+W6(11)*Q2+W6(12)*Q3+W6(13)*PQ
254        AC(8)=W8(1)+W8(2)*X+W8(3)*X2+W8(5)*Y+W8(6)*Y2+W8(7)*Y3+W8(8)*XY
255        AC(0)=-3.05
256       ELSEIF(RKAPPA .GE. 0.12) THEN
257        ITYPE=3
258        NPT=200
259        X=1+(RKAPPA-BKMXX2)*FBKX2
260        Y=1+(SQRT(BETA2)-BKMXY2)*FBKY2
261        XX=2*X
262        YY=2*Y
263        X2=XX*X-1
264        X3=XX*X2-X
265        Y2=YY*Y-1
266        Y3=YY*Y2-Y
267        XY=X*Y
268        P2=X2*Y
269        P3=X3*Y
270        Q2=Y2*X
271        Q3=Y3*X
272        PQ=X2*Y2
273        AC(1)=V1(1)+V1(2)*X+V1(3)*X2+V1(5)*Y+V1(6)*Y2+V1(7)*Y3+
274      1  V1(9)*P2+V1(10)*P3+V1(11)*Q2+V1(12)*Q3
275        AC(2)=V2(1)+V2(2)*X+V2(3)*X2+V2(5)*Y+V2(6)*Y2+V2(7)*Y3+
276      1  V2(8)*XY+V2(9)*P2+V2(11)*Q2+V2(12)*Q3
277        AC(3)=V3(1)+V3(2)*X+V3(3)*X2+V3(4)*X3+V3(5)*Y+V3(6)*Y2+V3(7)*Y3+
278      1  V3(8)*XY+V3(9)*P2+V3(10)*P3+V3(11)*Q2+V3(12)*Q3+V3(13)*PQ
279        AC(4)=V4(1)+V4(2)*X+V4(3)*X2+V4(4)*X3+V4(5)*Y+V4(6)*Y2+V4(7)*Y3+
280      1  V4(8)*XY+V4(9)*P2+V4(10)*P3+V4(11)*Q2+V4(12)*Q3
281        AC(5)=V5(1)+V5(2)*X+V5(3)*X2+V5(4)*X3+V5(5)*Y+V5(6)*Y2+V5(7)*Y3+
282      1  V5(8)*XY+V5(11)*Q2+V5(12)*Q3+V5(13)*PQ
283        AC(6)=V6(1)+V6(2)*X+V6(3)*X2+V6(4)*X3+V6(5)*Y+V6(6)*Y2+V6(7)*Y3+
284      1  V6(8)*XY+V6(9)*P2+V6(10)*P3+V6(11)*Q2+V6(12)*Q3+V6(13)*PQ
285        AC(7)=V7(1)+V7(2)*X+V7(3)*X2+V7(5)*Y+V7(6)*Y2+V7(7)*Y3+
286      1  V7(8)*XY+V7(11)*Q2
287        AC(8)=V8(1)+V8(2)*X+V8(3)*X2+V8(5)*Y+V8(6)*Y2+V8(7)*Y3+
288      1  V8(8)*XY+V8(11)*Q2
289        AC(0)=-3.04
290       ELSE
291        ITYPE=4
292        IF(RKAPPA .GE. 0.02) ITYPE=3
293        NPT=200
294        X=1+(RKAPPA-BKMXX1)*FBKX1
295        Y=1+(SQRT(BETA2)-BKMXY1)*FBKY1
296        XX=2*X
297        YY=2*Y
298        X2=XX*X-1
299        X3=XX*X2-X
300        Y2=YY*Y-1
301        Y3=YY*Y2-Y
302        XY=X*Y
303        P2=X2*Y
304        P3=X3*Y
305        Q2=Y2*X
306        Q3=Y3*X
307        PQ=X2*Y2
308        IF(ITYPE .EQ. 4) GO TO 4
309        AC(1)=U1(1)+U1(2)*X+U1(3)*X2+U1(5)*Y+U1(6)*Y2+U1(7)*Y3+
310      1  U1(8)*XY+U1(10)*P3+U1(12)*Q3+U1(13)*PQ
311        AC(2)=U2(1)+U2(2)*X+U2(3)*X2+U2(5)*Y+U2(6)*Y2+U2(7)*Y3+
312      1  U2(8)*XY+U2(9)*P2+U2(10)*P3+U2(12)*Q3+U2(13)*PQ
313        AC(3)=U3(1)+U3(2)*X+U3(3)*X2+U3(5)*Y+U3(6)*Y2+U3(7)*Y3+
314      1  U3(8)*XY+U3(9)*P2+U3(10)*P3+U3(11)*Q2+U3(12)*Q3+U3(13)*PQ
315        AC(4)=U4(1)+U4(2)*X+U4(3)*X2+U4(4)*X3+U4(5)*Y+U4(6)*Y2+U4(7)*Y3+
316      1  U4(8)*XY+U4(9)*P2+U4(10)*P3+U4(11)*Q2+U4(12)*Q3
317        AC(5)=U5(1)+U5(2)*X+U5(3)*X2+U5(4)*X3+U5(5)*Y+U5(6)*Y2+U5(7)*Y3+
318      1  U5(8)*XY+U5(10)*P3+U5(11)*Q2+U5(12)*Q3+U5(13)*PQ
319        AC(6)=U6(1)+U6(2)*X+U6(3)*X2+U6(4)*X3+U6(5)*Y+U6(7)*Y3+
320      1  U6(8)*XY+U6(9)*P2+U6(10)*P3+U6(12)*Q3+U6(13)*PQ
321     4  AC(7)=U7(1)+U7(2)*X+U7(3)*X2+U7(4)*X3+U7(5)*Y+U7(6)*Y2+U7(8)*XY
322        AC(8)=U8(1)+U8(2)*X+U8(3)*X2+U8(4)*X3+U8(5)*Y+U8(6)*Y2+U8(7)*Y3+
323      1  U8(8)*XY+U8(9)*P2+U8(10)*P3+U8(11)*Q2+U8(13)*PQ
324        AC(0)=-3.03
325       ENDIF
326       AC(9)=(AC(8)-AC(0))/NPT
327       IF(ITYPE .EQ. 3) THEN
328        X=(AC(7)-AC(8))/(AC(7)*AC(8))
329        Y=1/LOG(AC(8)/AC(7))
330        P2=AC(7)**2
331        AC(11)=P2*(AC(1)*EXP(-AC(2)*(AC(7)+AC(5)*P2)-
332      1            AC(3)*EXP(-AC(4)*(AC(7)+AC(6)*P2)))-0.045*Y/AC(7))/
333      2            (1+X*Y*AC(7))
334        AC(12)=(0.045+X*AC(11))*Y
335       ENDIF
336       IF(ITYPE .EQ. 4) AC(10)=0.995/GLANDS(AC(8))
337  
338       T=2*RAN/AC(9)
339       RLAM=AC(0)
340       FL=0
341       S=0
342       DO 21 N = 1,NPT
343       RLAM=RLAM+AC(9)
344       IF(ITYPE .EQ. 1) THEN
345        FN=1
346        X=(RLAM+HC(0))*HC(1)
347        H(1)=X
348        H(2)=X**2-1
349        DO 31 K = 2,8
350        FN=FN+1
351    31  H(K+1)=X*H(K)-FN*H(K-1)
352        Y=1+HC(7)*H(9)
353        DO 32 K = 2,6
354    32  Y=Y+HC(K)*H(K+1)
355        FU=HC(8)*EXP(-0.5*X**2)*MAX(Y,0.)
356       ELSEIF(ITYPE .EQ. 2) THEN
357        X=RLAM**2
358        FU=AC(1)*EXP(-AC(2)*(RLAM+AC(5)*X)-
359      1    AC(3)*EXP(-AC(4)*(RLAM+AC(6)*X)))
360       ELSEIF(ITYPE .EQ. 3) THEN
361        IF(RLAM .LT. AC(7)) THEN
362         X=RLAM**2
363         FU=AC(1)*EXP(-AC(2)*(RLAM+AC(5)*X)-
364      1     AC(3)*EXP(-AC(4)*(RLAM+AC(6)*X)))
365        ELSE
366         X=1/RLAM
367         FU=(AC(11)*X+AC(12))*X
368        ENDIF
369       ELSE
370        FU=AC(10)*GLANDE(RLAM)
371       ENDIF
372       S=S+FL+FU
373       IF(S .GT. T) GO TO 22
374    21 FL=FU
375  
376    22 S0=S-FL-FU
377       GVAVIV=RLAM-AC(9)
378       IF(S .GT. S0) GVAVIV=GVAVIV+AC(9)*(T-S0)/(S-S0)
379       RETURN
380       END