C++ Interface to Tauola
demo-factory/prod/taumain.F
1 PROGRAM taudem
2C **************
3C NOTE THAT THE ROUTINES ARE NOT LIKE IN CPC DECK THIS IS HISTORICAL !!
4C=======================================================================
5C====================== DECTES : TEST OF TAU DECAY LIBRARY===========
6C====================== KTORY = 1 : INTERFACE OF KORAL-Z TYPE ==========
7C====================== KTORY = 2 : INTERFACE OF KORAL-B TYPE =========
8C=======================================================================
9C COMMON /PAWC/ BLAN(10000)
10 COMMON / / blan(10000)
11 CHARACTER*7 DNAME
12 COMMON / inout / inut,iout
13 dname='KKPI'
14! CALL GLIMIT(20000)
15! CALL GOUTPU(16)
16 inut=5
17 iout=6
18 OPEN(iout,file="./tauola.output")
19 OPEN(inut,file="./dane.dat")
20 ktory=1
21 CALL dectes(ktory)
22 ktory=2
23 CALL dectes(ktory)
24 END
25 SUBROUTINE dectes(KTORY)
26C ************************
27 REAL POL(4)
28 DOUBLE PRECISION HH(4)
29C SWITCHES FOR TAUOLA;
30 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
31 COMMON / idfc / idff
32C I/O UNITS NUMBERS
33 COMMON / inout / inut,iout
34C LUND TYPE IDENTIFIER FOR A1
35 COMMON / idpart / ia1
36C /PTAU/ IS USED IN ROUTINE TRALO4
37 COMMON /ptau/ ptau
38 COMMON / taurad / xk0dec,itdkrc
39 real*8 xk0dec
40 COMMON /testa1/ keya1
41C special switch for tests of dGamma/dQ**2 in a1 decay
42C KEYA1=1 constant width of a1 and rho
43C KEYA1=2 free choice of rho propagator (defined in function FPIK)
44C and free choice of a1 mass and width. function g(Q**2)
45C (see formula 3.48 in Comp. Phys. Comm. 64 (1991) 275)
46C hard coded both in Monte Carlo and in testing distribution.
47C KEYA1=3 function g(Q**2) hardcoded in the Monte Carlo
48C (it is timy to calculate!), but appropriately adjusted in
49C testing distribution.
50C-----------------------------------------------------------------------
51C INITIALIZATION
52C-----------------------------------------------------------------------
53C======================================
54 ninp=inut
55 nout=iout
56 3000 FORMAT(a80)
57 3001 FORMAT(8i2)
58 3002 FORMAT(i10)
59 3003 FORMAT(f10.0)
60 IF (ktory.EQ.1) THEN
61 READ( ninp,3000) testit
62 WRITE(nout,3000) testit
63 READ( ninp,3001) kat1,kat2,kat3,kat4,kat5,kat6
64 READ( ninp,3002) nevt,jak1,jak2,itdkrc
65 READ( ninp,3003) ptau,xk0dec
66 ENDIF
67C======================================
68C control output
69 WRITE(nout,'(6A6/6I6)')
70 $ 'KAT1','KAT2','KAT3','KAT4','KAT5','KAT6',
71 $ kat1 , kat2 , kat3 , kat4 , kat5 , kat6
72 WRITE(nout,'(4A12/4I12)')
73 $ 'NEVT','JAK1','JAK2','ITDKRC',
74 $ nevt, jak1 , jak2 , itdkrc
75 WRITE(nout,'(2A12/2F12.6)')
76 $ 'PTAU','XK0DEC',
77 $ ptau , xk0dec
78C======================================
79 jak=0
80C JAK1=5
81C JAK2=5
82C LUND IDENTIFIER (FOR TAU+) -15
83 IF (ktory.EQ.1) THEN
84 idff=-15
85 ELSE
86 idff= 15
87 ENDIF
88C KTO=1 DENOTES TAU DEFINED BY IDFF (I.E. TAU+)
89C KTO=2 DENOTES THE OPPOSITE (I.E. TAU-)
90 kto=2
91 IF (kto.NE.2) THEN
92 print *, 'for the sake of these tests KTO has to be 2'
93 print *, 'to change tau- to tau+ change IDFF from -15 to 15'
94 stop
95 ENDIF
96C TAU POLARIZATION IN ITS RESTFRAME;
97 pol(1)=0.
98 pol(2)=0.
99 pol(3)=.9
100C TAU MOMENTUM IN GEV;
101C PTAU=CMSENE/2.D0
102C NUMBER OF EVENTS TO BE GENERATED;
103 nevtes=10
104 nevtes=nevt
105 print *, 'NEVTES= ',nevtes
106 WRITE(iout,7011) keya1
107C
108 IF (ktory.EQ.1) THEN
109 WRITE(iout,7001) jak,idff,pol(3),ptau
110 ELSE
111 WRITE(iout,7004) jak,idff,pol(3),ptau
112 ENDIF
113C INITIALISATION OF TAU DECAY PACKAGE TAUOLA
114C ******************************************
115 CALL inimas
116 CALL initdk
117
118
119 CALL iniphy(0.1d0)
120 IF (ktory.EQ.1) THEN
121 CALL dexay(-1,pol)
122 ELSE
123 CALL dekay(-1,hh)
124 ENDIF
125C-----------------------------------------------------------------------
126C GENERATION
127C-----------------------------------------------------------------------
128 nev=0
129 DO 300 iev=1,nevtes
130 nev=nev+1
131C RESLU INITIALISE THE LUND RECORD
132#if defined (history)
133 CALL reslu
134#else
135#endif
136 CALL taufil
137C DECAY....
138 IF (ktory.EQ.1) THEN
139 CALL dexay(kto,pol)
140 ELSE
141 CALL dekay(kto,hh)
142 CALL dekay(kto+10,hh)
143 ENDIF
144 CALL luhepc(2)
145 IF(iev.LE.44) THEN
146 WRITE(iout,7002) iev
147 IF (ktory.NE.1) THEN
148 WRITE(iout,7003) hh
149 ENDIF
150C CALL LULIST(11)
151 CALL lulist(2)
152 ENDIF
153 ipri=mod(nev,1000)
154 IF(ipri.EQ.1) print *, ' event no: ',nev,' NEVTES: ',nevtes
155 300 CONTINUE
156 301 CONTINUE
157C-----------------------------------------------------------------------
158C POSTGENERATION
159C-----------------------------------------------------------------------
160 IF (ktory.EQ.1) THEN
161 CALL dexay(100,pol)
162 ELSE
163 CALL dekay(100,hh)
164 ENDIF
165 RETURN
166 7001 FORMAT(//4(/1x,15(5h=====))
167 $ /,' ', 19x,' TEST OF RAD. CORR IN ELECTRON DECAY ',9x,1h ,
168 $ /,' ', 19x,' TESTS OF TAU DECAY ROUTINES ',9x,1h ,
169 $ /,' ', 19x,' INTERFACE OF THE KORAL-Z TYPE ',9x,1h ,
170 $ 2(/,1x,15(5h=====)),
171 $ /,5x ,'JAK =',i7 ,' KEY DEFINING DECAY TYPE ',9x,1h ,
172 $ /,5x ,'IDFF =',i7 ,' LUND IDENTIFIER FOR FIRST TAU ',9x,1h ,
173 $ /,5x ,'POL(3)=',f7.2,' THIRD COMPONENT OF TAU POLARIZ. ',9x,1h ,
174 $ /,5x ,'PTAU =',f7.2,' THIRD COMPONENT OF TAU MOM. GEV ',9x,1h ,
175 $ 2(/,1x,15(5h=====))/)
176 7002 FORMAT(///1x, '===== EVENT NO.',i4,1x,5h=====)
177 7003 FORMAT(5x,'POLARIMETRIC VECTOR: ',
178 $ 7x,'HH(1)',7x,'HH(2)',7x,'HH(3)',7x,'HH(4)',
179 $ /, 5x,' ', 4(1x,f11.8) )
180 7004 FORMAT(//4(/1x,15(5h=====))
181 $ /,' ', 19x,' TEST OF RAD. CORR IN ELECTRON DECAY ',9x,1h ,
182 $ /,' ', 19x,' TESTS OF TAU DECAY ROUTINES ',9x,1h ,
183 $ /,' ', 19x,' INTERFACE OF THE KORAL-B TYPE ',9x,1h ,
184 $ 2(/,1x,15(5h=====)),
185 $ /,5x ,'JAK =',i7 ,' KEY DEFINING DECAY TYPE ',9x,1h ,
186 $ /,5x ,'IDFF =',i7 ,' LUND IDENTIFIER FOR FIRST TAU ',9x,1h ,
187 $ /,5x ,'POL(3)=',f7.2,' THIRD COMPONENT OF TAU POLARIZ. ',9x,1h ,
188 $ /,5x ,'PTAU =',f7.2,' THIRD COMPONENT OF TAU MOM. GEV ',9x,1h ,
189 $ 2(/,1x,15(5h=====))/)
190 7011 FORMAT(///1x, '===== TYPE OF CURRENT',i4,1x,5h=====)
191 END
192 SUBROUTINE choice(MNUM,RR,ICHAN,PROB1,PROB2,PROB3,
193 $ AMRX,GAMRX,AMRA,GAMRA,AMRB,GAMRB)
194 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
195 * ,ampiz,ampi,amro,gamro,ama1,gama1
196 * ,amk,amkz,amkst,gamkst
197C
198 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
199 * ,ampiz,ampi,amro,gamro,ama1,gama1
200 * ,amk,amkz,amkst,gamkst
201C
202 amrop=1.1
203 gamrop=0.36
204 amom=.782
205 gamom=0.0084
206C XXXXA CORRESPOND TO S2 CHANNEL !
207 IF(mnum.EQ.0) THEN
208 prob1=0.5
209 prob2=0.5
210 amrx =ama1
211 gamrx=gama1
212 amra =amro
213 gamra=gamro
214 amrb =amro
215 gamrb=gamro
216 ELSEIF(mnum.EQ.1) THEN
217 prob1=0.5
218 prob2=0.5
219 amrx =1.57
220 gamrx=0.9
221 amrb =amkst
222 gamrb=gamkst
223 amra =amro
224 gamra=gamro
225 ELSEIF(mnum.EQ.2) THEN
226 prob1=0.5
227 prob2=0.5
228 amrx =1.57
229 gamrx=0.9
230 amrb =amkst
231 gamrb=gamkst
232 amra =amro
233 gamra=gamro
234 ELSEIF(mnum.EQ.3) THEN
235 prob1=0.5
236 prob2=0.5
237 amrx =1.27
238 gamrx=0.3
239 amra =amkst
240 gamra=gamkst
241 amrb =amkst
242 gamrb=gamkst
243 ELSEIF(mnum.EQ.4) THEN
244 prob1=0.5
245 prob2=0.5
246 amrx =1.27
247 gamrx=0.3
248 amra =amkst
249 gamra=gamkst
250 amrb =amkst
251 gamrb=gamkst
252 ELSEIF(mnum.EQ.5) THEN
253 prob1=0.5
254 prob2=0.5
255 amrx =1.27
256 gamrx=0.3
257 amra =amkst
258 gamra=gamkst
259 amrb =amro
260 gamrb=gamro
261 ELSEIF(mnum.EQ.6) THEN
262 prob1=0.4
263 prob2=0.4
264 amrx =1.27
265 gamrx=0.3
266 amra =amro
267 gamra=gamro
268 amrb =amkst
269 gamrb=gamkst
270 ELSEIF(mnum.EQ.7) THEN
271 prob1=0.0
272 prob2=1.0
273 amrx =1.27
274 gamrx=0.9
275 amra =amro
276 gamra=gamro
277 amrb =amro
278 gamrb=gamro
279 ELSEIF(mnum.EQ.8) THEN
280 prob1=0.0
281 prob2=1.0
282 amrx =amrop
283 gamrx=gamrop
284 amrb =amom
285 gamrb=gamom
286 amra =amro
287 gamra=gamro
288 ELSEIF(mnum.EQ.101) THEN
289 prob1=.35
290 prob2=.35
291 amrx =1.2
292 gamrx=.46
293 amrb =amom
294 gamrb=gamom
295 amra =amom
296 gamra=gamom
297 ELSEIF(mnum.EQ.102) THEN
298 prob1=0.0
299 prob2=0.0
300 amrx =1.4
301 gamrx=.6
302 amrb =amom
303 gamrb=gamom
304 amra =amom
305 gamra=gamom
306 ELSE
307 prob1=0.0
308 prob2=0.0
309 amrx =ama1
310 gamrx=gama1
311 amra =amro
312 gamra=gamro
313 amrb =amro
314 gamrb=gamro
315 ENDIF
316C
317 IF (rr.LE.prob1) THEN
318 ichan=1
319 ELSEIF(rr.LE.(prob1+prob2)) THEN
320 ichan=2
321 ax =amra
322 gx =gamra
323 amra =amrb
324 gamra=gamrb
325 amrb =ax
326 gamrb=gx
327 px =prob1
328 prob1=prob2
329 prob2=px
330 ELSE
331 ichan=3
332 ENDIF
333C
334 prob3=1.0-prob1-prob2
335 END
336 SUBROUTINE initdk
337* ----------------------------------------------------------------------
338* INITIALISATION OF TAU DECAY PARAMETERS and routines
339*
340* called by : KORALZ
341* ----------------------------------------------------------------------
342
343 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
344 real*4 gfermi,gv,ga,ccabib,scabib,gamel
345 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
346 * ,ampiz,ampi,amro,gamro,ama1,gama1
347 * ,amk,amkz,amkst,gamkst
348*
349 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
350 * ,ampiz,ampi,amro,gamro,ama1,gama1
351 * ,amk,amkz,amkst,gamkst
352 COMMON / taubra / gamprt(30),jlist(30),nchan
353 COMMON / taukle / bra1,brk0,brk0b,brks
354 real*4 bra1,brk0,brk0b,brks
355#if defined (ALEPH)
356 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
357 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
358 & ,names
359 CHARACTER NAMES(NMODE)*31
360#else
361 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
362 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
363 & ,names
364 CHARACTER NAMES(NMODE)*31
365#endif
366 CHARACTER OLDNAMES(7)*31
367 CHARACTER*80 bxINIT
368 parameter(
369 $ bxinit ='(1x,1h*,g17.8, 16x, a31,a4,a4, 1x,1h*)'
370 $ )
371 real*4 pi,pol1(4)
372*
373*
374* LIST OF BRANCHING RATIOS
375CAM normalised to e nu nutau channel
376CAM enu munu pinu rhonu A1nu Knu K*nu pi
377CAM DATA JLIST / 1, 2, 3, 4, 5, 6, 7,
378#if defined (ALEPH)
379CAM /0.1779,0.1731,0.1106,0.2530,0.1811,0.0072,0.0139
380CAM DATA GAMPRT / 1.000,0.9732,0.6217,1.4221,1.0180,0.0405,0.0781
381CAM DATA GAMPRT /1.000,0.9676,0.6154,1.3503,1.0225,0.0368,O.O758
382CAM
383C
384C conventions of particles names
385c
386cam mode (JAK) 8 9
387CAM channel pi- pi- pi0 pi+ 3pi0 pi-
388cam particle code -1,-1, 2, 1, 0, 0, 2, 2, 2,-1, 0, 0,
389CAM BR relative to electron .2414, .0601,
390c
391* 10 11
392* 1 3pi+- 2pi0 5pi+-
393* 1 -1,-1, 1, 2, 2, 0, -1,-1,-1, 1, 1, 0,
394* 1 .0281, .0045,
395
396* 12 13
397* 2 5pi+- pi0 3pi+- 3pi0
398* 2 -1,-1,-1, 1, 1, 2, -1,-1, 1, 2, 2, 2,
399* 2 .0010, .0062,
400
401* 14 15
402* 3 K- pi- K+ K0 pi- KB
403* 3 -3,-1, 3, 0, 0, 0, 4,-1,-4, 0, 0, 0,
404* 3 .0096, .0169,
405
406* 16 17
407* 4 K- pi0 K0 2pi0 K-
408* 4 -3, 2, 4, 0, 0, 0, 2, 2,-3, 0, 0, 0,
409* 4 .0056, .0045,
410
411* 18 19
412* 5 K- pi- pi+ pi- KB pi0
413* 5 -3,-1, 1, 0, 0, 0, -1,-4, 2, 0, 0, 0,
414* 5 .0219, .0180,
415
416* 20 21
417* 6 eta pi- pi0 pi- pi0 gamma
418* 6 9,-1, 2, 0, 0, 0, -1, 2, 8, 0, 0, 0,
419* 6 .0096, .0088,
420
421* 22 /
422* 7 K- K0 /
423* 7 -3, 4 /
424* 7 .0146 /
425#else
426*AM DATA GAMPRT /1.000,0.9730,0.6054,1.2432,0.8432,0.0432,O.O811,0.616
427*AM
428*AM multipion decays
429*
430* conventions of particles names
431* K-,P-,K+, K0,P-,KB, K-,P0,K0
432* 3, 1,-3 , 4, 1,-4 , 3, 2, 4 ,
433* P0,P0,K-, K-,P-,P+, P-,KB,P0
434* 2, 2, 3 , 3, 1,-1 , 1,-4, 2 ,
435* ET,P-,P0 P-,P0,GM
436* 9, 1, 2 , 1, 2, 8
437*
438#endif
439C
440 dimension nopik(6,nmode),npik(nmode)
441*AM outgoing multiplicity and flavors of multi-pion /multi-K modes
442 DATA npik / 4, 4,
443 1 5, 5,
444 2 6, 6,
445 3 3, 3,
446 4 3, 3,
447 5 3, 3,
448 6 3, 3,
449 7 2 /
450#if defined (ALEPH)
451 DATA nopik / -1,-1, 2, 1, 0, 0, 2, 2, 2,-1, 0, 0,
452 1 -1,-1, 1, 2, 2, 0, -1,-1,-1, 1, 1, 0,
453 2 -1,-1,-1, 1, 1, 2, -1,-1, 1, 2, 2, 2,
454 3 -3,-1, 3, 0, 0, 0, 4,-1,-4, 0, 0, 0,
455 4 -3, 2, 4, 0, 0, 0, 2, 2,-3, 0, 0, 0,
456 5 -3,-1, 1, 0, 0, 0, -1,-4, 2, 0, 0, 0,
457 6 9,-1, 2, 0, 0, 0, -1, 2, 8, 0, 0, 0,
458#else
459 DATA nopik / -1,-1, 1, 2, 0, 0, 2, 2, 2,-1, 0, 0,
460 1 -1,-1, 1, 2, 2, 0, -1,-1,-1, 1, 1, 0,
461 2 -1,-1,-1, 1, 1, 2, -1,-1, 1, 2, 2, 2,
462 3 -3,-1, 3, 0, 0, 0, -4,-1, 4, 0, 0, 0,
463 4 -3, 2,-4, 0, 0, 0, 2, 2,-3, 0, 0, 0,
464 5 -3,-1, 1, 0, 0, 0, -1, 4, 2, 0, 0, 0,
465 6 9,-1, 2, 0, 0, 0, -1, 2, 8, 0, 0, 0,
466#endif
467#if defined (CLEO)
468C AJWMOD fix sign bug, 2/22/99
469 7 -3,-4, 0, 0, 0, 0 /
470#else
471 7 -3, 4, 0, 0, 0, 0 /
472#endif
473* LIST OF BRANCHING RATIOS
474 nchan = nmode + 7
475 DO 1 i = 1,30
476 IF (i.LE.nchan) THEN
477 jlist(i) = i
478#if defined (CePeCe)
479 IF(i.EQ. 1) gamprt(i) = 1.0000
480 IF(i.EQ. 2) gamprt(i) = 1.0000
481 IF(i.EQ. 3) gamprt(i) = 1.0000
482 IF(i.EQ. 4) gamprt(i) = 1.0000
483 IF(i.EQ. 5) gamprt(i) = 1.0000
484 IF(i.EQ. 6) gamprt(i) = 1.0000
485 IF(i.EQ. 7) gamprt(i) = 1.0000
486 IF(i.EQ. 8) gamprt(i) = 1.0000
487 IF(i.EQ. 9) gamprt(i) = 1.0000
488 IF(i.EQ.10) gamprt(i) = 1.0000
489 IF(i.EQ.11) gamprt(i) = 1.0000
490 IF(i.EQ.12) gamprt(i) = 1.0000
491 IF(i.EQ.13) gamprt(i) = 1.0000
492 IF(i.EQ.14) gamprt(i) = 1.0000
493 IF(i.EQ.15) gamprt(i) = 1.0000
494 IF(i.EQ.16) gamprt(i) = 1.0000
495 IF(i.EQ.17) gamprt(i) = 1.0000
496 IF(i.EQ.18) gamprt(i) = 1.0000
497 IF(i.EQ.19) gamprt(i) = 1.0000
498 IF(i.EQ.20) gamprt(i) = 1.0000
499 IF(i.EQ.21) gamprt(i) = 1.0000
500 IF(i.EQ.22) gamprt(i) = 1.0000
501#elif defined (CLEO)
502 IF(i.EQ. 1) gamprt(i) =0.1800
503 IF(i.EQ. 2) gamprt(i) =0.1751
504 IF(i.EQ. 3) gamprt(i) =0.1110
505 IF(i.EQ. 4) gamprt(i) =0.2515
506 IF(i.EQ. 5) gamprt(i) =0.1790
507 IF(i.EQ. 6) gamprt(i) =0.0071
508 IF(i.EQ. 7) gamprt(i) =0.0134
509 IF(i.EQ. 8) gamprt(i) =0.0450
510 IF(i.EQ. 9) gamprt(i) =0.0100
511 IF(i.EQ.10) gamprt(i) =0.0009
512 IF(i.EQ.11) gamprt(i) =0.0004
513 IF(i.EQ.12) gamprt(i) =0.0003
514 IF(i.EQ.13) gamprt(i) =0.0005
515 IF(i.EQ.14) gamprt(i) =0.0015
516 IF(i.EQ.15) gamprt(i) =0.0015
517 IF(i.EQ.16) gamprt(i) =0.0015
518 IF(i.EQ.17) gamprt(i) =0.0005
519 IF(i.EQ.18) gamprt(i) =0.0050
520 IF(i.EQ.19) gamprt(i) =0.0055
521 IF(i.EQ.20) gamprt(i) =0.0017
522 IF(i.EQ.21) gamprt(i) =0.0013
523 IF(i.EQ.22) gamprt(i) =0.0010
524#elif defined (ALEPH)
525 IF(i.EQ. 1) gamprt(i) = 1.0000
526 IF(i.EQ. 2) gamprt(i) = .9732
527 IF(i.EQ. 3) gamprt(i) = .6217
528 IF(i.EQ. 4) gamprt(i) = 1.4221
529 IF(i.EQ. 5) gamprt(i) = 1.0180
530 IF(i.EQ. 6) gamprt(i) = .0405
531 IF(i.EQ. 7) gamprt(i) = .0781
532 IF(i.EQ. 8) gamprt(i) = .2414
533 IF(i.EQ. 9) gamprt(i) = .0601
534 IF(i.EQ.10) gamprt(i) = .0281
535 IF(i.EQ.11) gamprt(i) = .0045
536 IF(i.EQ.12) gamprt(i) = .0010
537 IF(i.EQ.13) gamprt(i) = .0062
538 IF(i.EQ.14) gamprt(i) = .0096
539 IF(i.EQ.15) gamprt(i) = .0169
540 IF(i.EQ.16) gamprt(i) = .0056
541 IF(i.EQ.17) gamprt(i) = .0045
542 IF(i.EQ.18) gamprt(i) = .0219
543 IF(i.EQ.19) gamprt(i) = .0180
544 IF(i.EQ.20) gamprt(i) = .0096
545 IF(i.EQ.21) gamprt(i) = .0088
546 IF(i.EQ.22) gamprt(i) = .0146
547#else
548#endif
549 IF(i.EQ. 1) oldnames(i)=' TAU- --> E- '
550 IF(i.EQ. 2) oldnames(i)=' TAU- --> MU- '
551 IF(i.EQ. 3) oldnames(i)=' TAU- --> PI- '
552 IF(i.EQ. 4) oldnames(i)=' TAU- --> PI-, PI0 '
553 IF(i.EQ. 5) oldnames(i)=' TAU- --> A1- (two subch) '
554 IF(i.EQ. 6) oldnames(i)=' TAU- --> K- '
555 IF(i.EQ. 7) oldnames(i)=' TAU- --> K*- (two subch) '
556 IF(i.EQ. 8) names(i-7)=' TAU- --> 2PI-, PI0, PI+ '
557 IF(i.EQ. 9) names(i-7)=' TAU- --> 3PI0, PI- '
558 IF(i.EQ.10) names(i-7)=' TAU- --> 2PI-, PI+, 2PI0 '
559 IF(i.EQ.11) names(i-7)=' TAU- --> 3PI-, 2PI+, '
560 IF(i.EQ.12) names(i-7)=' TAU- --> 3PI-, 2PI+, PI0 '
561 IF(i.EQ.13) names(i-7)=' TAU- --> 2PI-, PI+, 3PI0 '
562 IF(i.EQ.14) names(i-7)=' TAU- --> K-, PI-, K+ '
563 IF(i.EQ.15) names(i-7)=' TAU- --> K0, PI-, K0B '
564#if defined (ALEPH)
565 IF(i.EQ.16) names(i-7)=' TAU- --> K- PI0 K0 '
566#else
567 IF(i.EQ.16) names(i-7)=' TAU- --> K-, K0, PI0 '
568#endif
569 IF(i.EQ.17) names(i-7)=' TAU- --> PI0 PI0 K- '
570 IF(i.EQ.18) names(i-7)=' TAU- --> K- PI- PI+ '
571 IF(i.EQ.19) names(i-7)=' TAU- --> PI- K0B PI0 '
572 IF(i.EQ.20) names(i-7)=' TAU- --> ETA PI- PI0 '
573 IF(i.EQ.21) names(i-7)=' TAU- --> PI- PI0 GAM '
574 IF(i.EQ.22) names(i-7)=' TAU- --> K- K0 '
575 ELSE
576 jlist(i) = 0
577 gamprt(i) = 0.
578 ENDIF
579 1 CONTINUE
580 DO i=1,nmode
581 mulpik(i)=npik(i)
582 DO j=1,mulpik(i)
583 idffin(j,i)=nopik(j,i)
584 ENDDO
585 ENDDO
586*
587*
588* --- COEFFICIENTS TO FIX RATIO OF:
589* --- A1 3CHARGED/ A1 1CHARGED 2 NEUTRALS MATRIX ELEMENTS (MASLESS LIM.)
590* --- PROBABILITY OF K0 TO BE KS
591* --- PROBABILITY OF K0B TO BE KS
592* --- RATIO OF COEFFICIENTS FOR K*--> K0 PI-
593* --- ALL COEFFICENTS SHOULD BE IN THE RANGE (0.0,1.0)
594* --- THEY MEANING IS PROBABILITY OF THE FIRST CHOICE ONLY IF ONE
595* --- NEGLECTS MASS-PHASE SPACE EFFECTS
596 bra1=0.5
597 brk0=0.5
598 brk0b=0.5
599 brks=0.6667
600*
601
602 gfermi = 1.16637e-5
603 ccabib = 0.975
604 gv = 1.0
605 ga =-1.0
606
607
608
609* ZW 13.04.89 HERE WAS AN ERROR
610 scabib = sqrt(1.-ccabib**2)
611 pi =4.*atan(1.)
612 gamel = gfermi**2*amtau**5/(192*pi**3)
613*
614* CALL DEXAY(-1,pol1)
615*
616 RETURN
617 END
618 FUNCTION dcdmas(IDENT)
619 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
620 * ,ampiz,ampi,amro,gamro,ama1,gama1
621 * ,amk,amkz,amkst,gamkst
622*
623 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
624 * ,ampiz,ampi,amro,gamro,ama1,gama1
625 * ,amk,amkz,amkst,gamkst
626 IF (ident.EQ. 1) THEN
627 apkmas=ampi
628 ELSEIF (ident.EQ.-1) THEN
629 apkmas=ampi
630 ELSEIF (ident.EQ. 2) THEN
631 apkmas=ampiz
632 ELSEIF (ident.EQ.-2) THEN
633 apkmas=ampiz
634 ELSEIF (ident.EQ. 3) THEN
635 apkmas=amk
636 ELSEIF (ident.EQ.-3) THEN
637 apkmas=amk
638 ELSEIF (ident.EQ. 4) THEN
639 apkmas=amkz
640 ELSEIF (ident.EQ.-4) THEN
641 apkmas=amkz
642 ELSEIF (ident.EQ. 8) THEN
643 apkmas=0.0001
644 ELSEIF (ident.EQ.-8) THEN
645 apkmas=0.0001
646 ELSEIF (ident.EQ. 9) THEN
647 apkmas=0.5488
648 ELSEIF (ident.EQ.-9) THEN
649 apkmas=0.5488
650 ELSE
651 print *, 'STOP IN APKMAS, WRONG IDENT=',ident
652 stop
653 ENDIF
654 dcdmas=apkmas
655 END
656 FUNCTION lunpik(ID,ISGN)
657 COMMON / taukle / bra1,brk0,brk0b,brks
658 real*4 bra1,brk0,brk0b,brks
659 real*4 xio(1)
660 ident=id*isgn
661#if defined (ALEPH)
662 IF (ident.EQ. 1) THEN
663 ipkdef= 211
664 ELSEIF (ident.EQ.-1) THEN
665 ipkdef=-211
666 ELSEIF (ident.EQ. 2) THEN
667 ipkdef= 111
668 ELSEIF (ident.EQ.-2) THEN
669 ipkdef= 111
670 ELSEIF (ident.EQ. 3) THEN
671 ipkdef= 321
672 ELSEIF (ident.EQ.-3) THEN
673 ipkdef=-321
674#else
675 IF (ident.EQ. 1) THEN
676 ipkdef=-211
677 ELSEIF (ident.EQ.-1) THEN
678 ipkdef= 211
679 ELSEIF (ident.EQ. 2) THEN
680 ipkdef=111
681 ELSEIF (ident.EQ.-2) THEN
682 ipkdef=111
683 ELSEIF (ident.EQ. 3) THEN
684 ipkdef=-321
685 ELSEIF (ident.EQ.-3) THEN
686 ipkdef= 321
687#endif
688 ELSEIF (ident.EQ. 4) THEN
689*
690* K0 --> K0_LONG (IS 130) / K0_SHORT (IS 310) = 1/1
691 CALL ranmar(xio,1)
692 IF (xio(1).GT.brk0) THEN
693 ipkdef= 130
694 ELSE
695 ipkdef= 310
696 ENDIF
697 ELSEIF (ident.EQ.-4) THEN
698*
699* K0B--> K0_LONG (IS 130) / K0_SHORT (IS 310) = 1/1
700 CALL ranmar(xio,1)
701 IF (xio(1).GT.brk0b) THEN
702 ipkdef= 130
703 ELSE
704 ipkdef= 310
705 ENDIF
706 ELSEIF (ident.EQ. 8) THEN
707 ipkdef= 22
708 ELSEIF (ident.EQ.-8) THEN
709 ipkdef= 22
710 ELSEIF (ident.EQ. 9) THEN
711 ipkdef= 221
712 ELSEIF (ident.EQ.-9) THEN
713 ipkdef= 221
714 ELSE
715 print *, 'STOP IN IPKDEF, WRONG IDENT=',ident
716 stop
717 ENDIF
718 lunpik=ipkdef
719 END
720
721
722#if defined (CLEO)
723
724 SUBROUTINE taurdf(KTO)
725C THIS ROUTINE CAN BE CALLED BEFORE ANY TAU+ OR TAU- EVENT IS GENERATED
726C IT CAN BE USED TO GENERATE TAU+ AND TAU- SAMPLES OF DIFFERENT
727C CONTENTS
728 COMMON / taukle / bra1,brk0,brk0b,brks
729 real*4 bra1,brk0,brk0b,brks
730 COMMON / taubra / gamprt(30),jlist(30),nchan
731 IF (kto.EQ.1) THEN
732C ==================
733C AJWMOD: Set the BRs for (A1+ -> rho+ pi0) and (K*+ -> K0 pi+)
734 bra1 = pkorb(4,1)
735 brks = pkorb(4,3)
736 brk0 = pkorb(4,5)
737 brk0b = pkorb(4,6)
738 ELSE
739C ====
740C AJWMOD: Set the BRs for (A1+ -> rho+ pi0) and (K*+ -> K0 pi+)
741 bra1 = pkorb(4,2)
742 brks = pkorb(4,4)
743 brk0 = pkorb(4,5)
744 brk0b = pkorb(4,6)
745 ENDIF
746C =====
747 END
748#else
749
750 SUBROUTINE taurdf(KTO)
751* THIS ROUTINE CAN BE CALLED BEFORE ANY TAU+ OR TAU- EVENT IS GENERATED
752* IT CAN BE USED TO GENERATE TAU+ AND TAU- SAMPLES OF DIFFERENT
753* CONTENTS
754 COMMON / taukle / bra1,brk0,brk0b,brks
755 real*4 bra1,brk0,brk0b,brks
756 COMMON / taubra / gamprt(30),jlist(30),nchan
757 IF (kto.EQ.1) THEN
758* ==================
759* LIST OF BRANCHING RATIOS
760 nchan = 19
761 DO 1 i = 1,30
762 IF (i.LE.nchan) THEN
763 jlist(i) = i
764 IF(i.EQ. 1) gamprt(i) = .0000
765 IF(i.EQ. 2) gamprt(i) = .0000
766 IF(i.EQ. 3) gamprt(i) = .0000
767 IF(i.EQ. 4) gamprt(i) = .0000
768 IF(i.EQ. 5) gamprt(i) = .0000
769 IF(i.EQ. 6) gamprt(i) = .0000
770 IF(i.EQ. 7) gamprt(i) = .0000
771 IF(i.EQ. 8) gamprt(i) = 1.0000
772 IF(i.EQ. 9) gamprt(i) = 1.0000
773 IF(i.EQ.10) gamprt(i) = 1.0000
774 IF(i.EQ.11) gamprt(i) = 1.0000
775 IF(i.EQ.12) gamprt(i) = 1.0000
776 IF(i.EQ.13) gamprt(i) = 1.0000
777 IF(i.EQ.14) gamprt(i) = 1.0000
778 IF(i.EQ.15) gamprt(i) = 1.0000
779 IF(i.EQ.16) gamprt(i) = 1.0000
780 IF(i.EQ.17) gamprt(i) = 1.0000
781 IF(i.EQ.18) gamprt(i) = 1.0000
782 IF(i.EQ.19) gamprt(i) = 1.0000
783 ELSE
784 jlist(i) = 0
785 gamprt(i) = 0.
786 ENDIF
787 1 CONTINUE
788* --- COEFFICIENTS TO FIX RATIO OF:
789* --- A1 3CHARGED/ A1 1CHARGED 2 NEUTRALS MATRIX ELEMENTS (MASLESS LIM.)
790* --- PROBABILITY OF K0 TO BE KS
791* --- PROBABILITY OF K0B TO BE KS
792* --- RATIO OF COEFFICIENTS FOR K*--> K0 PI-
793* --- ALL COEFFICENTS SHOULD BE IN THE RANGE (0.0,1.0)
794* --- THEY MEANING IS PROBABILITY OF THE FIRST CHOICE ONLY IF ONE
795* --- NEGLECTS MASS-PHASE SPACE EFFECTS
796 bra1=0.5
797 brk0=0.5
798 brk0b=0.5
799 brks=0.6667
800 ELSE
801* ====
802* LIST OF BRANCHING RATIOS
803 nchan = 19
804 DO 2 i = 1,30
805 IF (i.LE.nchan) THEN
806 jlist(i) = i
807 IF(i.EQ. 1) gamprt(i) = .0000
808 IF(i.EQ. 2) gamprt(i) = .0000
809 IF(i.EQ. 3) gamprt(i) = .0000
810 IF(i.EQ. 4) gamprt(i) = .0000
811 IF(i.EQ. 5) gamprt(i) = .0000
812 IF(i.EQ. 6) gamprt(i) = .0000
813 IF(i.EQ. 7) gamprt(i) = .0000
814 IF(i.EQ. 8) gamprt(i) = 1.0000
815 IF(i.EQ. 9) gamprt(i) = 1.0000
816 IF(i.EQ.10) gamprt(i) = 1.0000
817 IF(i.EQ.11) gamprt(i) = 1.0000
818 IF(i.EQ.12) gamprt(i) = 1.0000
819 IF(i.EQ.13) gamprt(i) = 1.0000
820 IF(i.EQ.14) gamprt(i) = 1.0000
821 IF(i.EQ.15) gamprt(i) = 1.0000
822 IF(i.EQ.16) gamprt(i) = 1.0000
823 IF(i.EQ.17) gamprt(i) = 1.0000
824 IF(i.EQ.18) gamprt(i) = 1.0000
825 IF(i.EQ.19) gamprt(i) = 1.0000
826 ELSE
827 jlist(i) = 0
828 gamprt(i) = 0.
829 ENDIF
830 2 CONTINUE
831* --- COEFFICIENTS TO FIX RATIO OF:
832* --- A1 3CHARGED/ A1 1CHARGED 2 NEUTRALS MATRIX ELEMENTS (MASLESS LIM.)
833* --- PROBABILITY OF K0 TO BE KS
834* --- PROBABILITY OF K0B TO BE KS
835* --- RATIO OF COEFFICIENTS FOR K*--> K0 PI-
836* --- ALL COEFFICENTS SHOULD BE IN THE RANGE (0.0,1.0)
837* --- THEY MEANING IS PROBABILITY OF THE FIRST CHOICE ONLY IF ONE
838* --- NEGLECTS MASS-PHASE SPACE EFFECTS
839 bra1=0.5
840 brk0=0.5
841 brk0b=0.5
842 brks=0.6667
843 ENDIF
844* =====
845 END
846#endif
847
848 SUBROUTINE iniphy(XK00)
849* ----------------------------------------------------------------------
850* INITIALISATION OF PARAMETERS
851* USED IN QED and/or GSW ROUTINES
852* ----------------------------------------------------------------------
853 COMMON / qedprm /alfinv,alfpi,xk0
854 real*8 alfinv,alfpi,xk0
855 real*8 pi8,xk00
856*
857 pi8 = 4.d0*datan(1.d0)
858 alfinv = 137.03604d0
859 alfpi = 1d0/(alfinv*pi8)
860 xk0=xk00
861 END
862
863 SUBROUTINE inimas
864C ----------------------------------------------------------------------
865C INITIALISATION OF MASSES
866C
867C called by : KORALZ
868C ----------------------------------------------------------------------
869 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
870 * ,ampiz,ampi,amro,gamro,ama1,gama1
871 * ,amk,amkz,amkst,gamkst
872*
873 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
874 * ,ampiz,ampi,amro,gamro,ama1,gama1
875 * ,amk,amkz,amkst,gamkst
876C
877C IN-COMING / OUT-GOING FERMION MASSES
878 amtau = 1.7842
879C --- let us update tau mass ...
880 amtau = 1.777
881 amnuta = 0.010
882 amel = 0.0005111
883 amnue = 0.0
884 ammu = 0.105659
885 amnumu = 0.0
886*
887* MASSES USED IN TAU DECAYS
888#if defined (CePeCe)
889 ampiz = 0.134964
890 ampi = 0.139568
891 amro = 0.773
892 gamro = 0.145
893*C GAMRO = 0.666
894 ama1 = 1.251
895 gama1 = 0.599
896 amk = 0.493667
897 amkz = 0.49772
898 amkst = 0.8921
899 gamkst = 0.0513
900#elif defined (CLEO)
901 ampiz = 0.134964
902 ampi = 0.139568
903 amro = 0.773
904 gamro = 0.145
905*C GAMRO = 0.666
906 ama1 = 1.251
907 gama1 = 0.599
908 amk = 0.493667
909 amkz = 0.49772
910 amkst = 0.8921
911 gamkst = 0.0513
912C
913C
914C IN-COMING / OUT-GOING FERMION MASSES
915!! AMNUTA = PKORB(1,2)
916!! AMNUE = PKORB(1,4)
917!! AMNUMU = PKORB(1,6)
918C
919C MASSES USED IN TAU DECAYS Cleo settings
920!! AMPIZ = PKORB(1,7)
921!! AMPI = PKORB(1,8)
922!! AMRO = PKORB(1,9)
923!! GAMRO = PKORB(2,9)
924 ama1 = 1.275 !! PKORB(1,10)
925 gama1 = 0.615 !! PKORB(2,10)
926!! AMK = PKORB(1,11)
927!! AMKZ = PKORB(1,12)
928!! AMKST = PKORB(1,13)
929!! GAMKST = PKORB(2,13)
930C
931#elif defined (ALEPH)
932 ampiz = 0.134964
933 ampi = 0.139568
934 amro = 0.7714
935 gamro = 0.153
936cam AMRO = 0.773
937cam GAMRO = 0.145
938 ama1 = 1.251! PMAS(LUCOMP(ia1),1) ! AMA1 = 1.251
939 gama1 = 0.599! PMAS(LUCOMP(ia1),2) ! GAMA1 = 0.599
940 print *,'INIMAS a1 mass= ',ama1,gama1
941 amk = 0.493667
942 amkz = 0.49772
943 amkst = 0.8921
944 gamkst = 0.0513
945#else
946#endif
947
948 RETURN
949 END
950 SUBROUTINE taufil
951C *****************
952C SUBSITUTE OF tau PRODUCTION GENERATOR
953C
954 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
955 * ,ampiz,ampi,amro,gamro,ama1,gama1
956 * ,amk,amkz,amkst,gamkst
957C
958 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
959 * ,ampiz,ampi,amro,gamro,ama1,gama1
960 * ,amk,amkz,amkst,gamkst
961 COMMON / idfc / idff
962C positions of taus in the LUND common block
963C it will be used by TAUOLA output routines.
964 COMMON /taupos / npa,npb
965 dimension xpb1(4),xpb2(4),aqf1(4),aqf2(4)
966C
967C --- DEFINING DUMMY EVENTS MOMENTA
968 DO 4 k=1,3
969 xpb1(k)=0.0
970 xpb2(k)=0.0
971 aqf1(k)=0.0
972 aqf2(k)=0.0
973 4 CONTINUE
974 aqf1(4)=amtau
975 aqf2(4)=amtau
976C --- TAU MOMENTA
977 CALL tralo4(1,aqf1,aqf1,am)
978 CALL tralo4(2,aqf2,aqf2,am)
979C --- BEAMS MOMENTA AND IDENTIFIERS
980 kfb1= 11*idff/iabs(idff)
981 kfb2=-11*idff/iabs(idff)
982 xpb1(4)= aqf1(4)
983 xpb1(3)= aqf1(4)
984 IF(aqf1(3).NE.0.0)
985 $ xpb1(3)= aqf1(4)*aqf1(3)/abs(aqf1(3))
986 xpb2(4)= aqf2(4)
987 xpb2(3)=-aqf2(4)
988 IF(aqf2(3).NE.0.0)
989 $ xpb2(3)= aqf2(4)*aqf2(3)/abs(aqf2(3))
990C --- Position of first and second tau in LUND common
991 npa=3
992 npb=4
993C --- FILL TO LUND COMMON
994 CALL filhep( 1,3, kfb1,0,0,0,0,xpb1, amel,.true.)
995 CALL filhep( 2,3, kfb2,0,0,0,0,xpb2, amel,.true.)
996 CALL filhep(npa,1, idff,1,2,0,0,aqf1,amtau,.true.)
997 CALL filhep(npb,1,-idff,1,2,0,0,aqf2,amtau,.true.)
998 END
999 SUBROUTINE tralo4(KTO,P,Q,AM)
1000C **************************
1001C SUBSITUTE OF TRALO4
1002 REAL P(4),Q(4)
1003C
1004 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1005 * ,ampiz,ampi,amro,gamro,ama1,gama1
1006 * ,amk,amkz,amkst,gamkst
1007C
1008 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1009 * ,ampiz,ampi,amro,gamro,ama1,gama1
1010 * ,amk,amkz,amkst,gamkst
1011 COMMON /ptau/ ptau
1012 am=amas4(p)
1013 etau=sqrt(ptau**2+amtau**2)
1014 exe=(etau+ptau)/amtau
1015 IF(kto.EQ.2) exe=(etau-ptau)/amtau
1016 CALL bostr3(exe,p,q)
1017C ======================================================================
1018C END OF THE TEST JOB
1019C ======================================================================
1020 END
1021 SUBROUTINE filhep(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
1022C ----------------------------------------------------------------------
1023C this subroutine fills one entry into the HEPEVT common
1024C and updates the information for affected mother entries
1025C
1026C written by Martin W. Gruenewald (91/01/28)
1027C
1028C called by : ZTOHEP,BTOHEP,DWLUxy
1029C ----------------------------------------------------------------------
1030C
1031#include "../../include/HEPEVT.h"
1032C PARAMETER (NMXHEP=2000)
1033C COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
1034C &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
1035C SAVE /HEPEVT/
1036C COMMON/PHOQED/QEDRAD(NMXHEP)
1037C LOGICAL QEDRAD
1038C SAVE /PHOQED/
1039 LOGICAL PHFLAG
1040C
1041 real*4 p4(4)
1042C
1043C check address mode
1044 IF (n.EQ.0) THEN
1045C
1046C append mode
1047 ihep=nhep+1
1048 ELSE IF (n.GT.0) THEN
1049C
1050C absolute position
1051 ihep=n
1052 ELSE
1053C
1054C relative position
1055 ihep=nhep+n
1056 END IF
1057C
1058C check on IHEP
1059 IF ((ihep.LE.0).OR.(ihep.GT.nmxhep)) RETURN
1060C
1061C add entry
1062 nhep=ihep
1063 isthep(ihep)=ist
1064 idhep(ihep)=id
1065 jmohep(1,ihep)=jmo1
1066 IF(jmo1.LT.0)jmohep(1,ihep)=jmohep(1,ihep)+ihep
1067 jmohep(2,ihep)=jmo2
1068 IF(jmo2.LT.0)jmohep(2,ihep)=jmohep(2,ihep)+ihep
1069 jdahep(1,ihep)=jda1
1070 jdahep(2,ihep)=jda2
1071C
1072 DO i=1,4
1073 phep(i,ihep)=p4(i)
1074C
1075C KORAL-B and KORAL-Z do not provide vertex and/or lifetime informations
1076 vhep(i,ihep)=0.0
1077 END DO
1078 phep(5,ihep)=pinv
1079C FLAG FOR PHOTOS...
1080 qedrad(ihep)=phflag
1081C
1082C update process:
1083 DO ip=jmohep(1,ihep),jmohep(2,ihep)
1084 IF(ip.GT.0)THEN
1085C
1086C if there is a daughter at IHEP, mother entry at IP has decayed
1087 IF(isthep(ip).EQ.1)isthep(ip)=2
1088C
1089C and daughter pointers of mother entry must be updated
1090 IF(jdahep(1,ip).EQ.0)THEN
1091 jdahep(1,ip)=ihep
1092 jdahep(2,ip)=ihep
1093 ELSE
1094 jdahep(2,ip)=max(ihep,jdahep(2,ip))
1095 END IF
1096 END IF
1097 END DO
1098C
1099 RETURN
1100 END
1101#if defined (ALEPH)
1102 FUNCTION dilogy(X)
1103C *****************
1104 IMPLICIT REAL*8(a-h,o-z)
1105CERN C304 VERSION 29/07/71 DILOG 59 C
1106 z=-1.64493406684822
1107 IF(x .LT.-1.0) GO TO 1
1108 IF(x .LE. 0.5) GO TO 2
1109 IF(x .EQ. 1.0) GO TO 3
1110 IF(x .LE. 2.0) GO TO 4
1111 z=3.2898681336964
1112 1 t=1.0/x
1113 s=-0.5
1114 z=z-0.5* log(abs(x))**2
1115 GO TO 5
1116 2 t=x
1117 s=0.5
1118 z=0.
1119 GO TO 5
1120 3 dilogy=1.64493406684822
1121 RETURN
1122 4 t=1.0-x
1123 s=-0.5
1124 z=1.64493406684822 - log(x)* log(abs(t))
1125 5 y=2.66666666666666 *t+0.66666666666666
1126 b= 0.00000 00000 00001
1127 a=y*b +0.00000 00000 00004
1128 b=y*a-b+0.00000 00000 00011
1129 a=y*b-a+0.00000 00000 00037
1130 b=y*a-b+0.00000 00000 00121
1131 a=y*b-a+0.00000 00000 00398
1132 b=y*a-b+0.00000 00000 01312
1133 a=y*b-a+0.00000 00000 04342
1134 b=y*a-b+0.00000 00000 14437
1135 a=y*b-a+0.00000 00000 48274
1136 b=y*a-b+0.00000 00001 62421
1137 a=y*b-a+0.00000 00005 50291
1138 b=y*a-b+0.00000 00018 79117
1139 a=y*b-a+0.00000 00064 74338
1140 b=y*a-b+0.00000 00225 36705
1141 a=y*b-a+0.00000 00793 87055
1142 b=y*a-b+0.00000 02835 75385
1143 a=y*b-a+0.00000 10299 04264
1144 b=y*a-b+0.00000 38163 29463
1145 a=y*b-a+0.00001 44963 00557
1146 b=y*a-b+0.00005 68178 22718
1147 a=y*b-a+0.00023 20021 96094
1148 b=y*a-b+0.00100 16274 96164
1149 a=y*b-a+0.00468 63619 59447
1150 b=y*a-b+0.02487 93229 24228
1151 a=y*b-a+0.16607 30329 27855
1152 a=y*a-b+1.93506 43008 6996
1153 dilogy=s*t*(a-b)+z
1154 RETURN
1155C=======================================================================
1156C===================END OF CPC PART ====================================
1157C=======================================================================
1158 END
1159#endif