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