]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HERWIG/jimmy/cernmisc/jmdbska.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
In ESDtoAOD: updating instead of creating new AliAODHeader.
[u/mrichter/AliRoot.git] / HERWIG / jimmy / cernmisc / jmdbska.F
CommitLineData
ef94df36 1*
2* $Id$
3*
4* $Log$
a0a69b14 5* Revision 1.1 2006/10/02 10:45:20 morsch
6* JIMMY first commit.
7*
ef94df36 8* Revision 1.1 2005/07/29 16:15:46 jmb
9* Include the various CERNLIB functions jimmy needs
10*
11* Revision 1.1.1.1 1996/04/01 15:02:08 mclareni
12* Mathlib gen
13*
14*
15 SUBROUTINE JMDBSKA(X,IA,JA,NL,B)
16 IMPLICIT DOUBLE PRECISION (A-H,J,O-Z)
17 SAVE
18 LOGICAL LEX
a0a69b14 19 INTEGER JA,J
ef94df36 20
21 CHARACTER NAME*(*),ENAM*(*)
22 CHARACTER*80 ERRTXT
23 PARAMETER (NAME = 'BSKA/JMDBSKA', ENAM = 'EBSKA/DEBKA')
24
25 PARAMETER (Z1 = 1, Z2 = 2, Z3 = 3, Z4 = 4)
26 PARAMETER (Z12 = Z1/Z2, Z13 = Z1/Z3, Z14 = Z1/Z4, Z23 = Z2/Z3)
27 PARAMETER (Z34 = Z3/Z4)
28
29 DIMENSION B(0:*)
30
31 PARAMETER (PI = 3.14159 26535 89793D0, PIV = PI/4)
32
33 LEX=.FALSE.
34 GO TO 9
35
36 ENTRY JMDEBKA(X,IA,JA,NL,B)
37 LEX=.TRUE.
38
39 9 MODE=10*IA+JA
40 N=NL-1
41 U=2/X
42 IF(LEX) THEN
43 IF(X .LE. 0) THEN
44 N=0
45 WRITE(ERRTXT,101) X
46 CALL JMMLPT(ENAM,'C341.1',ERRTXT)
47 ELSEIF(NL .LT. 0 .OR. NL .GT. 100) THEN
48 N=0
49 WRITE(ERRTXT,103) NL
50 CALL JMMLPT(ENAM,'C341.3',ERRTXT)
51 ELSEIF(IA .EQ. 0) THEN
52 A=0
53 B(0)=JMDEBSK0(X)
54 B(1)=JMDEBSK1(X)
55 ELSEIF(MODE .EQ. 12) THEN
56 A=Z12
57 B(0)=SQRT(PIV*U)
58 B(1)=B(0)*(1+A*U)
59 ELSEIF(MODE .EQ. 13) THEN
60 A=Z13
61 B(0)=JMDEBKR3(X,1)
62 B(1)=JMDEBKR3(X,2)+A*U*B(0)
63 ELSEIF(MODE .EQ. 14) THEN
64 A=Z14
65 B(0)=JMDEBKR4(X,1)
66 B(1)=JMDEBKR4(X,3)+A*U*B(0)
67 ELSEIF(MODE .EQ. 23) THEN
68 A=Z23
69 B(0)=JMDEBKR3(X,2)
70 B(1)=JMDEBKR3(X,1)+A*U*B(0)
71 ELSEIF(MODE .EQ. 34) THEN
72 A=Z34
73 B(0)=JMDEBKR4(X,3)
74 B(1)=JMDEBKR4(X,1)+A*U*B(0)
75 ELSE
76 N=0
77 WRITE(ERRTXT,102) IA,JA
78 CALL JMMLPT(ENAM,'C341.2',ERRTXT)
79 ENDIF
80 ELSE
81 IF(X .LE. 0) THEN
82 N=0
83 WRITE(ERRTXT,101) X
84 CALL JMMLPT(NAME,'C341.1',ERRTXT)
85 ELSEIF(NL .LT. 0 .OR. NL .GT. 100) THEN
86 N=0
87 WRITE(ERRTXT,103) NL
88 CALL JMMLPT(NAME,'C341.3',ERRTXT)
89 ELSEIF(IA .EQ. 0) THEN
90 A=0
91 B(0)=JMDBESK0(X)
92 B(1)=JMDBESK1(X)
93 ELSEIF(MODE .EQ. 12) THEN
94 A=Z12
95 B(0)=EXP(-X)*SQRT(PIV*U)
96 B(1)=B(0)*(1+A*U)
97 ELSEIF(MODE .EQ. 13) THEN
98 A=Z13
99 B(0)=JMDBSKR3(X,1)
100 B(1)=JMDBSKR3(X,2)+A*U*B(0)
101 ELSEIF(MODE .EQ. 14) THEN
102 A=Z14
103 B(0)=JMDBSKR4(X,1)
104 B(1)=JMDBSKR4(X,3)+A*U*B(0)
105 ELSEIF(MODE .EQ. 23) THEN
106 A=Z23
107 B(0)=JMDBSKR3(X,2)
108 B(1)=JMDBSKR3(X,1)+A*U*B(0)
109 ELSEIF(MODE .EQ. 34) THEN
110 A=Z34
111 B(0)=JMDBSKR4(X,3)
112 B(1)=JMDBSKR4(X,1)+A*U*B(0)
113 ELSE
114 N=0
115 WRITE(ERRTXT,102) IA,JA
116 CALL JMMLPT(NAME,'C341.2',ERRTXT)
117 ENDIF
118 ENDIF
119 DO 1 J = 1,N
120 A=A+1
121 1 B(J+1)=B(J-1)+A*U*B(J)
122 RETURN
123 101 FORMAT('NON-POSITIVE ARGUMENT X = ',E15.6)
124 102 FORMAT('PAIR (IA,JA) = (',I5,I5,') ILLEGAL')
125 103 FORMAT('ILLEGAL NL =',I5)
126 END
127
128
129* $Id$
130*
131* $Log$
a0a69b14 132* Revision 1.1 2006/10/02 10:45:20 morsch
133* JIMMY first commit.
134*
ef94df36 135* Revision 1.1 2005/07/29 16:15:46 jmb
136* Include the various CERNLIB functions jimmy needs
137*
138* Revision 1.1.1.1 1996/04/01 15:02:00 mclareni
139* Mathlib gen
140*
141*
142 DOUBLE PRECISION FUNCTION JMDBESI0(X)
143 IMPLICIT DOUBLE PRECISION (A-H,J,O-Z)
144 SAVE
145
146 LOGICAL LEX
147 CHARACTER NAME0*(*),NAME1*(*),NAME0E*(*),NAME1E*(*)
148 CHARACTER*80 ERRTXT
149 DIMENSION CI(0:24,0:1),CK(0:16,0:1)
150
151 PARAMETER (NAME0 = 'JMBESK0/JMDBESK0', NAME0E =
152 $ 'JMEBESK0/JMDBESK0')
153 PARAMETER (NAME1 = 'JMBESK1/JMDBESK1', NAME1E =
154 $ 'JMEBESK1/JMDBESK1')
155
156 PARAMETER (EPS=1D-15)
157 PARAMETER (Z1 = 1, HF = Z1/2)
158 PARAMETER (PI = 3.14159 26535 89793 24D0)
159 PARAMETER (CE = 0.57721 56649 01532 86D0)
160 PARAMETER (PIH = PI/2, RPIH = 2/PI, RPI2 = 1/(2*PI))
161
162 DATA CI( 0,0) /+1.00827 92054 58740 032D0/
163 DATA CI( 1,0) /+0.00844 51226 24920 943D0/
164 DATA CI( 2,0) /+0.00017 27006 30777 567D0/
165 DATA CI( 3,0) /+0.00000 72475 91099 959D0/
166 DATA CI( 4,0) /+0.00000 05135 87726 878D0/
167 DATA CI( 5,0) /+0.00000 00568 16965 808D0/
168 DATA CI( 6,0) /+0.00000 00085 13091 223D0/
169 DATA CI( 7,0) /+0.00000 00012 38425 364D0/
170 DATA CI( 8,0) /+0.00000 00000 29801 672D0/
171 DATA CI( 9,0) /-0.00000 00000 78956 698D0/
172 DATA CI(10,0) /-0.00000 00000 33127 128D0/
173 DATA CI(11,0) /-0.00000 00000 04497 339D0/
174 DATA CI(12,0) /+0.00000 00000 01799 790D0/
175 DATA CI(13,0) /+0.00000 00000 00965 748D0/
176 DATA CI(14,0) /+0.00000 00000 00038 604D0/
177 DATA CI(15,0) /-0.00000 00000 00104 039D0/
178 DATA CI(16,0) /-0.00000 00000 00023 950D0/
179 DATA CI(17,0) /+0.00000 00000 00009 554D0/
180 DATA CI(18,0) /+0.00000 00000 00004 443D0/
181 DATA CI(19,0) /-0.00000 00000 00000 859D0/
182 DATA CI(20,0) /-0.00000 00000 00000 709D0/
183 DATA CI(21,0) /+0.00000 00000 00000 087D0/
184 DATA CI(22,0) /+0.00000 00000 00000 112D0/
185 DATA CI(23,0) /-0.00000 00000 00000 012D0/
186 DATA CI(24,0) /-0.00000 00000 00000 018D0/
187
188 DATA CI( 0,1) /+0.97580 06023 26285 926D0/
189 DATA CI( 1,1) /-0.02446 74429 63276 385D0/
190 DATA CI( 2,1) /-0.00027 72053 60763 829D0/
191 DATA CI( 3,1) /-0.00000 97321 46728 020D0/
192 DATA CI( 4,1) /-0.00000 06297 24238 640D0/
193 DATA CI( 5,1) /-0.00000 00659 61142 154D0/
194 DATA CI( 6,1) /-0.00000 00096 13872 919D0/
195 DATA CI( 7,1) /-0.00000 00014 01140 901D0/
196 DATA CI( 8,1) /-0.00000 00000 47563 167D0/
197 DATA CI( 9,1) /+0.00000 00000 81530 681D0/
198 DATA CI(10,1) /+0.00000 00000 35408 148D0/
199 DATA CI(11,1) /+0.00000 00000 05102 564D0/
200 DATA CI(12,1) /-0.00000 00000 01804 409D0/
201 DATA CI(13,1) /-0.00000 00000 01023 594D0/
202 DATA CI(14,1) /-0.00000 00000 00052 678D0/
203 DATA CI(15,1) /+0.00000 00000 00107 094D0/
204 DATA CI(16,1) /+0.00000 00000 00026 120D0/
205 DATA CI(17,1) /-0.00000 00000 00009 561D0/
206 DATA CI(18,1) /-0.00000 00000 00004 713D0/
207 DATA CI(19,1) /+0.00000 00000 00000 829D0/
208 DATA CI(20,1) /+0.00000 00000 00000 743D0/
209 DATA CI(21,1) /-0.00000 00000 00000 080D0/
210 DATA CI(22,1) /-0.00000 00000 00000 117D0/
211 DATA CI(23,1) /+0.00000 00000 00000 011D0/
212 DATA CI(24,1) /+0.00000 00000 00000 019D0/
213
214 DATA CK( 0,0) /+0.98840 81742 30825 800D0/
215 DATA CK( 1,0) /-0.01131 05046 46928 281D0/
216 DATA CK( 2,0) /+0.00026 95326 12762 724D0/
217 DATA CK( 3,0) /-0.00001 11066 85196 665D0/
218 DATA CK( 4,0) /+0.00000 06325 75108 500D0/
219 DATA CK( 5,0) /-0.00000 00450 47337 641D0/
220 DATA CK( 6,0) /+0.00000 00037 92996 456D0/
221 DATA CK( 7,0) /-0.00000 00003 64547 179D0/
222 DATA CK( 8,0) /+0.00000 00000 39043 756D0/
223 DATA CK( 9,0) /-0.00000 00000 04579 936D0/
224 DATA CK(10,0) /+0.00000 00000 00580 811D0/
225 DATA CK(11,0) /-0.00000 00000 00078 832D0/
226 DATA CK(12,0) /+0.00000 00000 00011 360D0/
227 DATA CK(13,0) /-0.00000 00000 00001 727D0/
228 DATA CK(14,0) /+0.00000 00000 00000 275D0/
229 DATA CK(15,0) /-0.00000 00000 00000 046D0/
230 DATA CK(16,0) /+0.00000 00000 00000 008D0/
231
232 DATA CK( 0,1) /+1.03595 08587 72358 331D0/
233 DATA CK( 1,1) /+0.03546 52912 43331 114D0/
234 DATA CK( 2,1) /-0.00046 84750 28166 889D0/
235 DATA CK( 3,1) /+0.00001 61850 63810 053D0/
236 DATA CK( 4,1) /-0.00000 08451 72048 124D0/
237 DATA CK( 5,1) /+0.00000 00571 32218 103D0/
238 DATA CK( 6,1) /-0.00000 00046 45554 607D0/
239 DATA CK( 7,1) /+0.00000 00004 35417 339D0/
240 DATA CK( 8,1) /-0.00000 00000 45757 297D0/
241 DATA CK( 9,1) /+0.00000 00000 05288 133D0/
242 DATA CK(10,1) /-0.00000 00000 00662 613D0/
243 DATA CK(11,1) /+0.00000 00000 00089 048D0/
244 DATA CK(12,1) /-0.00000 00000 00012 726D0/
245 DATA CK(13,1) /+0.00000 00000 00001 921D0/
246 DATA CK(14,1) /-0.00000 00000 00000 305D0/
247 DATA CK(15,1) /+0.00000 00000 00000 050D0/
248 DATA CK(16,1) /-0.00000 00000 00000 009D0/
249
250 NU=0
251 LEX=.FALSE.
252 GO TO 6
253
254 ENTRY JMDEBSI0(X)
255
256 NU=0
257 LEX=.TRUE.
258 GO TO 6
259
260 ENTRY JMDBESI1(X)
261
262 NU=1
263 LEX=.FALSE.
264 GO TO 6
265
266 ENTRY JMDEBSI1(X)
267
268 NU=1
269 LEX=.TRUE.
270
271 6 V=ABS(X)
272 IF(V .LT. 8) THEN
273 Y=(HF*V)**2
274 XL=NU+2
275 A0=1
276 A1=1+2*Y/((XL+1)*(XL-1))
277 A2=1+Y*(4+3*Y/((XL+2)*XL))/((XL+3)*(XL-1))
278 B0=1
279 B1=1-Y/(XL+1)
280 B2=1-Y*(1-Y/(2*(XL+2)))/(XL+3)
281 W1=3+XL
282 V1=3-XL
283 V3=XL-1
284 V2=V3+V3
285 C=0
286 DO 3 N = 3,30
287 C0=C
288 FN=N
289 W1=W1+2
290 W2=W1-1
291 W3=W2-1
292 W4=W3-1
293 W5=W4-1
294 W6=W5-1
295 V1=V1+1
296 V2=V2+1
297 V3=V3+1
298 U1=FN*W4
299 E=V3/(U1*W3)
300 U2=E*Y
301 F1=1+Y*V1/(U1*W1)
302 F2=(1+Y*V2/(V3*W2*W5))*U2
303 F3=-Y*Y*U2/(W4*W5*W5*W6)
304 A=F1*A2+F2*A1+F3*A0
305 B=F1*B2+F2*B1+F3*B0
306 C=A/B
307 IF(ABS(C0-C) .LT. EPS*ABS(C)) GO TO 4
308 A0=A1
309 A1=A2
310 A2=A
311 B0=B1
312 B1=B2
313 B2=B
314 3 CONTINUE
315 4 H=C
316 IF(NU .EQ. 1) H=HF*X*H
317 IF(LEX) H=EXP(-V)*H
318 ELSE
319 R=1/V
320 H=16*R-1
321 ALFA=H+H
322 B1=0
323 B2=0
324 DO 1 I = 24,0,-1
325 B0=CI(I,NU)+ALFA*B1-B2
326 B2=B1
327 1 B1=B0
328 H=SQRT(RPI2*R)*(B0-H*B2)
329 IF(NU*X .LT. 0) H=-H
330 IF(.NOT.LEX) H=EXP(V)*H
331 ENDIF
332 GO TO 9
333
334 ENTRY JMDBESK0(X)
335
336 NU=0
337 LEX=.FALSE.
338 GO TO 8
339
340 ENTRY JMDEBSK0(X)
341
342 NU=0
343 LEX=.TRUE.
344 GO TO 8
345
346 ENTRY JMDBESK1(X)
347
348 NU=1
349 LEX=.FALSE.
350 GO TO 8
351
352 ENTRY JMDEBSK1(X)
353
354 NU=1
355 LEX=.TRUE.
356
357 8 IF(X .LE. 0) THEN
358 H=0
359 WRITE(ERRTXT,101) X
360 IF(NU .EQ. 0 .AND. .NOT.LEX) CALL JMMLPT(NAME0 ,'C313.1',ERRTXT)
361 IF(NU .EQ. 0 .AND. LEX) CALL JMMLPT(NAME0E,'C313.1',ERRTXT)
362 IF(NU .EQ. 1 .AND. .NOT.LEX) CALL JMMLPT(NAME1 ,'C313.1',ERRTXT)
363 IF(NU .EQ. 1 .AND. LEX) CALL JMMLPT(NAME1E,'C313.1',ERRTXT)
364 ELSEIF(X .LT. 1) THEN
365 B=HF*X
366 BK=-(LOG(B)+CE)
367 F=BK
368 P=HF
369 Q=HF
370 C=1
371 D=B**2
372 BK1=P
373 DO 11 N = 1,15
374 FN=N
375 RFN=1/FN
376 P=P*RFN
377 Q=Q*RFN
378 F=(F+P+Q)*RFN
379 C=C*D*RFN
380 G=C*(P-FN*F)
381 H=C*F
382 BK=BK+H
383 BK1=BK1+G
384 IF(BK1*H+ABS(G)*BK .LE. EPS*BK*BK1) GO TO 12
385 11 CONTINUE
386 12 H=BK
387 IF(NU .EQ. 1) H=BK1/B
388 IF(LEX) H=EXP(X)*H
389 ELSEIF(X .LE. 5) THEN
390 XN=4*NU**2
391 A=9-XN
392 B=25-XN
393 C=768*X**2
394 C0=48*X
395 A0=1
396 A1=(16*X+7+XN)/A
397 A2=(C+C0*(XN+23)+XN*(XN+62)+129)/(A*B)
398 B0=1
399 B1=(16*X+9-XN)/A
400 B2=(C+C0*B)/(A*B)+1
401 C=0
402 DO 24 N = 3,30
403 C0=C
404 FN=N
405 FN2=FN+FN
406 FN1=FN2-1
407 FN3=FN1/(FN2-3)
408 FN4=12*FN**2-(1-XN)
409 FN5=16*FN1*X
410 RAN=1/((FN2+1)**2-XN)
411 F1=FN3*(FN4-20*FN)+FN5
412 F2=28*FN-FN4-8+FN5
413 F3=FN3*((FN2-5)**2-XN)
414 A=(F1*A2+F2*A1+F3*A0)*RAN
415 B=(F1*B2+F2*B1+F3*B0)*RAN
416 C=A/B
417 IF(ABS(C0-C) .LT. EPS*ABS(C)) GO TO 25
418 A0=A1
419 A1=A2
420 A2=A
421 B0=B1
422 B1=B2
423 B2=B
424 24 CONTINUE
425 25 H=C/SQRT(RPIH*X)
426 IF(.NOT.LEX) H=EXP(-X)*H
427 ELSE
428 R=1/X
429 H=10*R-1
430 ALFA=H+H
431 B1=0
432 B2=0
433 DO 23 I = 16,0,-1
434 B0=CK(I,NU)+ALFA*B1-B2
435 B2=B1
436 23 B1=B0
437 H=SQRT(PIH*R)*(B0-H*B2)
438 IF(.NOT.LEX) H=EXP(-X)*H
439 ENDIF
440 9 CONTINUE
441 JMDBESI0=H
442 RETURN
443 101 FORMAT(' NON-POSITIVE ARGUMENT X = ',1P,E15.6)
a0a69b14 444 END