C++ Interface to Tauola
demo-factory/back/taumain.F
1  PROGRAM taudem
2 C **************
3 C NOTE THAT THE ROUTINES ARE NOT LIKE IN CPC DECK THIS IS HISTORICAL !!
4 C=======================================================================
5 C====================== DECTES : TEST OF TAU DECAY LIBRARY===========
6 C====================== KTORY = 1 : INTERFACE OF KORAL-Z TYPE ==========
7 C====================== KTORY = 2 : INTERFACE OF KORAL-B TYPE =========
8 C=======================================================================
9 C 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)
26 C ************************
27  REAL POL(4)
28  DOUBLE PRECISION HH(4)
29 C SWITCHES FOR TAUOLA;
30  COMMON / jaki / jak1,jak2,jakp,jakm,ktom
31  COMMON / idfc / idff
32 C I/O UNITS NUMBERS
33  COMMON / inout / inut,iout
34 C LUND TYPE IDENTIFIER FOR A1
35  COMMON / idpart / ia1
36 C /PTAU/ IS USED IN ROUTINE TRALO4
37  COMMON /ptau/ ptau
38  COMMON / taurad / xk0dec,itdkrc
39  real*8 xk0dec
40  COMMON /testa1/ keya1
41 C special switch for tests of dGamma/dQ**2 in a1 decay
42 C KEYA1=1 constant width of a1 and rho
43 C KEYA1=2 free choice of rho propagator (defined in function FPIK)
44 C and free choice of a1 mass and width. function g(Q**2)
45 C (see formula 3.48 in Comp. Phys. Comm. 64 (1991) 275)
46 C hard coded both in Monte Carlo and in testing distribution.
47 C KEYA1=3 function g(Q**2) hardcoded in the Monte Carlo
48 C (it is timy to calculate!), but appropriately adjusted in
49 C testing distribution.
50 C-----------------------------------------------------------------------
51 C INITIALIZATION
52 C-----------------------------------------------------------------------
53 C======================================
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
67 C======================================
68 C 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
78 C======================================
79  jak=0
80 C JAK1=5
81 C JAK2=5
82 C LUND IDENTIFIER (FOR TAU+) -15
83  IF (ktory.EQ.1) THEN
84  idff=-15
85  ELSE
86  idff= 15
87  ENDIF
88 C KTO=1 DENOTES TAU DEFINED BY IDFF (I.E. TAU+)
89 C 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
96 C TAU POLARIZATION IN ITS RESTFRAME;
97  pol(1)=0.
98  pol(2)=0.
99  pol(3)=.9
100 C TAU MOMENTUM IN GEV;
101 C PTAU=CMSENE/2.D0
102 C NUMBER OF EVENTS TO BE GENERATED;
103  nevtes=10
104  nevtes=nevt
105  print *, 'NEVTES= ',nevtes
106  WRITE(iout,7011) keya1
107 C
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
113 C INITIALISATION OF TAU DECAY PACKAGE TAUOLA
114 C ******************************************
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
125 C-----------------------------------------------------------------------
126 C GENERATION
127 C-----------------------------------------------------------------------
128  nev=0
129  DO 300 iev=1,nevtes
130  nev=nev+1
131 C RESLU INITIALISE THE LUND RECORD
132 #if defined (history)
133  CALL reslu
134 #else
135 #endif
136  CALL taufil
137 C 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
150 C 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
157 C-----------------------------------------------------------------------
158 C POSTGENERATION
159 C-----------------------------------------------------------------------
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
197 C
198  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
199  * ,ampiz,ampi,amro,gamro,ama1,gama1
200  * ,amk,amkz,amkst,gamkst
201 C
202  amrop=1.1
203  gamrop=0.36
204  amom=.782
205  gamom=0.0084
206 C 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
316 C
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
333 C
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
375 CAM normalised to e nu nutau channel
376 CAM enu munu pinu rhonu A1nu Knu K*nu pi
377 CAM DATA JLIST / 1, 2, 3, 4, 5, 6, 7,
378 #if defined (ALEPH)
379 CAM /0.1779,0.1731,0.1106,0.2530,0.1811,0.0072,0.0139
380 CAM DATA GAMPRT / 1.000,0.9732,0.6217,1.4221,1.0180,0.0405,0.0781
381 CAM DATA GAMPRT /1.000,0.9676,0.6154,1.3503,1.0225,0.0368,O.O758
382 CAM
383 C
384 C conventions of particles names
385 c
386 cam mode (JAK) 8 9
387 CAM channel pi- pi- pi0 pi+ 3pi0 pi-
388 cam particle code -1,-1, 2, 1, 0, 0, 2, 2, 2,-1, 0, 0,
389 CAM BR relative to electron .2414, .0601,
390 c
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
439 C
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)
468 C 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)
725 C THIS ROUTINE CAN BE CALLED BEFORE ANY TAU+ OR TAU- EVENT IS GENERATED
726 C IT CAN BE USED TO GENERATE TAU+ AND TAU- SAMPLES OF DIFFERENT
727 C 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
732 C ==================
733 C 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
739 C ====
740 C 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
746 C =====
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
864 C ----------------------------------------------------------------------
865 C INITIALISATION OF MASSES
866 C
867 C called by : KORALZ
868 C ----------------------------------------------------------------------
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
876 C
877 C IN-COMING / OUT-GOING FERMION MASSES
878  amtau = 1.7842
879 C --- 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
912 C
913 C
914 C IN-COMING / OUT-GOING FERMION MASSES
915 !! AMNUTA = PKORB(1,2)
916 !! AMNUE = PKORB(1,4)
917 !! AMNUMU = PKORB(1,6)
918 C
919 C 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)
930 C
931 #elif defined (ALEPH)
932  ampiz = 0.134964
933  ampi = 0.139568
934  amro = 0.7714
935  gamro = 0.153
936 cam AMRO = 0.773
937 cam 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
951 C *****************
952 C SUBSITUTE OF tau PRODUCTION GENERATOR
953 C
954  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
955  * ,ampiz,ampi,amro,gamro,ama1,gama1
956  * ,amk,amkz,amkst,gamkst
957 C
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
962 C positions of taus in the LUND common block
963 C it will be used by TAUOLA output routines.
964  COMMON /taupos / npa,npb
965  dimension xpb1(4),xpb2(4),aqf1(4),aqf2(4)
966 C
967 C --- 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
976 C --- TAU MOMENTA
977  CALL tralo4(1,aqf1,aqf1,am)
978  CALL tralo4(2,aqf2,aqf2,am)
979 C --- 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))
990 C --- Position of first and second tau in LUND common
991  npa=3
992  npb=4
993 C --- 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)
1000 C **************************
1001 C SUBSITUTE OF TRALO4
1002  REAL P(4),Q(4)
1003 C
1004  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1005  * ,ampiz,ampi,amro,gamro,ama1,gama1
1006  * ,amk,amkz,amkst,gamkst
1007 C
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)
1017 C ======================================================================
1018 C END OF THE TEST JOB
1019 C ======================================================================
1020  END
1021  SUBROUTINE filhep(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
1022 C ----------------------------------------------------------------------
1023 C this subroutine fills one entry into the HEPEVT common
1024 C and updates the information for affected mother entries
1025 C
1026 C written by Martin W. Gruenewald (91/01/28)
1027 C
1028 C called by : ZTOHEP,BTOHEP,DWLUxy
1029 C ----------------------------------------------------------------------
1030 C
1031 #include "../../include/HEPEVT.h"
1032 C PARAMETER (NMXHEP=2000)
1033 C COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
1034 C &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
1035 C SAVE /HEPEVT/
1036 C COMMON/PHOQED/QEDRAD(NMXHEP)
1037 C LOGICAL QEDRAD
1038 C SAVE /PHOQED/
1039  LOGICAL PHFLAG
1040 C
1041  real*4 p4(4)
1042 C
1043 C check address mode
1044  IF (n.EQ.0) THEN
1045 C
1046 C append mode
1047  ihep=nhep+1
1048  ELSE IF (n.GT.0) THEN
1049 C
1050 C absolute position
1051  ihep=n
1052  ELSE
1053 C
1054 C relative position
1055  ihep=nhep+n
1056  END IF
1057 C
1058 C check on IHEP
1059  IF ((ihep.LE.0).OR.(ihep.GT.nmxhep)) RETURN
1060 C
1061 C 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
1071 C
1072  DO i=1,4
1073  phep(i,ihep)=p4(i)
1074 C
1075 C 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
1079 C FLAG FOR PHOTOS...
1080  qedrad(ihep)=phflag
1081 C
1082 C update process:
1083  DO ip=jmohep(1,ihep),jmohep(2,ihep)
1084  IF(ip.GT.0)THEN
1085 C
1086 C if there is a daughter at IHEP, mother entry at IP has decayed
1087  IF(isthep(ip).EQ.1)isthep(ip)=2
1088 C
1089 C 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
1098 C
1099  RETURN
1100  END
1101 #if defined (ALEPH)
1102  FUNCTION dilogy(X)
1103 C *****************
1104  IMPLICIT REAL*8(a-h,o-z)
1105 CERN 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
1155 C=======================================================================
1156 C===================END OF CPC PART ====================================
1157 C=======================================================================
1158  END
1159 #endif