C++ Interface to Tauola
demo-factory/back/attic/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( 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)
27 C ************************
28  REAL POL(4)
29  DOUBLE PRECISION HH(4)
30 C SWITCHES FOR TAUOLA;
31  COMMON / jaki / jak1,jak2,jakp,jakm,ktom
32  COMMON / idfc / idff
33 C I/O UNITS NUMBERS
34  COMMON / inout / inut,iout
35 C LUND TYPE IDENTIFIER FOR A1
36  COMMON / idpart / ia1
37 C /PTAU/ IS USED IN ROUTINE TRALO4
38  COMMON /ptau/ ptau
39  COMMON / taurad / xk0dec,itdkrc
40  real*8 xk0dec
41  COMMON /testa1/ keya1
42 C special switch for tests of dGamma/dQ**2 in a1 decay
43 C KEYA1=1 constant width of a1 and rho
44 C KEYA1=2 free choice of rho propagator (defined in function FPIK)
45 C and free choice of a1 mass and width. function g(Q**2)
46 C (see formula 3.48 in Comp. Phys. Comm. 64 (1991) 275)
47 C hard coded both in Monte Carlo and in testing distribution.
48 C KEYA1=3 function g(Q**2) hardcoded in the Monte Carlo
49 C (it is timy to calculate!), but appropriately adjusted in
50 C testing distribution.
51 C-----------------------------------------------------------------------
52 C INITIALIZATION
53 C-----------------------------------------------------------------------
54 C======================================
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
68 C======================================
69 C 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
79 C======================================
80  jak=0
81 C JAK1=5
82 C JAK2=5
83 C LUND IDENTIFIER (FOR TAU+) -15
84  IF (ktory.EQ.1) THEN
85  idff=-15
86  ELSE
87  idff= 15
88  ENDIF
89 C KTO=1 DENOTES TAU DEFINED BY IDFF (I.E. TAU+)
90 C 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
97 C TAU POLARIZATION IN ITS RESTFRAME;
98  pol(1)=0.
99  pol(2)=0.
100  pol(3)=.9
101 C TAU MOMENTUM IN GEV;
102 C PTAU=CMSENE/2.D0
103 C NUMBER OF EVENTS TO BE GENERATED;
104  nevtes=10
105  nevtes=nevt
106  print *, 'NEVTES= ',nevtes
107  WRITE(iout,7011) keya1
108 C
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
114 C INITIALISATION OF TAU DECAY PACKAGE TAUOLA
115 C ******************************************
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
126 C-----------------------------------------------------------------------
127 C GENERATION
128 C-----------------------------------------------------------------------
129  nev=0
130  DO 300 iev=1,nevtes
131  nev=nev+1
132 C RESLU INITIALISE THE LUND RECORD
133 #if defined (history)
134  CALL reslu
135 #else
136 #endif
137  CALL taufil
138 C 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
151 C 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
158 C-----------------------------------------------------------------------
159 C POSTGENERATION
160 C-----------------------------------------------------------------------
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
198 C
199  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
200  * ,ampiz,ampi,amro,gamro,ama1,gama1
201  * ,amk,amkz,amkst,gamkst
202 C
203  amrop=1.1
204  gamrop=0.36
205  amom=.782
206  gamom=0.0084
207 C 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
317 C
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
334 C
335  prob3=1.0-prob1-prob2
336  END
337  SUBROUTINE initdk
338 C ----------------------------------------------------------------------
339 C INITIALISATION OF TAU DECAY PARAMETERS and routines
340 C
341 C called by : KORALZ
342 C ----------------------------------------------------------------------
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 C
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
367 C
368 C LIST OF BRANCHING RATIOS
369 CAM normalised to e nu nutau channel
370 CAM enu munu pinu rhonu A1nu Knu K*nu pi
371 CAM DATA JLIST / 1, 2, 3, 4, 5, 6, 7,
372 #if defined (ALEPH)
373 CAM /0.1779,0.1731,0.1106,0.2530,0.1811,0.0072,0.0139
374 CAM DATA GAMPRT / 1.000,0.9732,0.6217,1.4221,1.0180,0.0405,0.0781
375 CAM DATA GAMPRT /1.000,0.9676,0.6154,1.3503,1.0225,0.0368,O.O758
376 CAM
377 C
378 C conventions of particles names
379 c
380 cam mode (JAK) 8 9
381 CAM channel pi- pi- pi0 pi+ 3pi0 pi-
382 cam particle code -1,-1, 2, 1, 0, 0, 2, 2, 2,-1, 0, 0,
383 CAM BR relative to electron .2414, .0601,
384 c
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
433 C
434  dimension nopik(6,nmode),npik(nmode)
435 CAM 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)
462 C 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
467 C 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
595 C
596 C
597 C --- COEFFICIENTS TO FIX RATIO OF:
598 C --- A1 3CHARGED/ A1 1CHARGED 2 NEUTRALS MATRIX ELEMENTS (MASLESS LIM.)
599 C --- PROBABILITY OF K0 TO BE KS
600 C --- PROBABILITY OF K0B TO BE KS
601 C --- RATIO OF COEFFICIENTS FOR K*--> K0 PI-
602 C --- ALL COEFFICENTS SHOULD BE IN THE RANGE (0.0,1.0)
603 C --- THEY MEANING IS PROBABILITY OF THE FIRST CHOICE ONLY IF ONE
604 C --- NEGLECTS MASS-PHASE SPACE EFFECTS
605  bra1=0.5
606  brk0=0.5
607  brk0b=0.5
608  brks=0.6667
609 C
610 C --- remaining constants
611  pi =4.*atan(1.)
612  gfermi = 1.16637e-5
613  ccabib = 0.975
614  gv = 1.0
615  ga =-1.0
616 C ZW 13.04.89 HERE WAS AN ERROR
617  scabib = sqrt(1.-ccabib**2)
618  gamel = gfermi**2*amtau**5/(192*pi**3)
619 C
620 C CALL DEXAY(-1)
621 C
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
628 C
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
696 C
697 C 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
705 C
706 C 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)
730 C THIS ROUTINE CAN BE CALLED BEFORE ANY TAU+ OR TAU- EVENT IS GENERATED
731 C IT CAN BE USED TO GENERATE TAU+ AND TAU- SAMPLES OF DIFFERENT
732 C 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
737 C ==================
738 C 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
744 C ====
745 C 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
751 C =====
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)
853 C ----------------------------------------------------------------------
854 C INITIALISATION OF PARAMETERS
855 C USED IN QED and/or GSW ROUTINES
856 C ----------------------------------------------------------------------
857  COMMON / qedprm /alfinv,alfpi,xk0
858  real*8 alfinv,alfpi,xk0
859  real*8 pi8,xk00
860 C
861  pi8 = 4.d0*datan(1.d0)
862  alfinv = 137.03604d0
863  alfpi = 1d0/(alfinv*pi8)
864  xk0=xk00
865  END
866  SUBROUTINE inimas
867 C ----------------------------------------------------------------------
868 C INITIALISATION OF MASSES
869 C
870 C called by : KORALZ
871 C ----------------------------------------------------------------------
872  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
873  * ,ampiz,ampi,amro,gamro,ama1,gama1
874  * ,amk,amkz,amkst,gamkst
875 C
876  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
877  * ,ampiz,ampi,amro,gamro,ama1,gama1
878  * ,amk,amkz,amkst,gamkst
879 C
880 C 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
887 C
888 C MASSES USED IN TAU DECAYS
889  ampiz = 0.134964
890  ampi = 0.139568
891  amro = 0.773
892  gamro = 0.145
893 CC 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 C
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
925 C
926 C
927 C IN-COMING / OUT-GOING FERMION MASSES
928 !! AMNUTA = PKORB(1,2)
929 !! AMNUE = PKORB(1,4)
930 !! AMNUMU = PKORB(1,6)
931 C
932 C 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)
943 C
944 #elif defined (ALEPH)
945  ampiz = 0.134964
946  ampi = 0.139568
947  amro = 0.7714
948  gamro = 0.153
949 cam AMRO = 0.773
950 cam 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
964 C *****************
965 C SUBSITUTE OF tau PRODUCTION GENERATOR
966 C
967  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
968  * ,ampiz,ampi,amro,gamro,ama1,gama1
969  * ,amk,amkz,amkst,gamkst
970 C
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
975 C positions of taus in the LUND common block
976 C it will be used by TAUOLA output routines.
977  COMMON /taupos / npa,npb
978  dimension xpb1(4),xpb2(4),aqf1(4),aqf2(4)
979 C
980 C --- 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
989 C --- TAU MOMENTA
990  CALL tralo4(1,aqf1,aqf1,am)
991  CALL tralo4(2,aqf2,aqf2,am)
992 C --- 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))
1003 C --- Position of first and second tau in LUND common
1004  npa=3
1005  npb=4
1006 C --- 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)
1013 C **************************
1014 C SUBSITUTE OF TRALO4
1015  REAL P(4),Q(4)
1016 C
1017  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1018  * ,ampiz,ampi,amro,gamro,ama1,gama1
1019  * ,amk,amkz,amkst,gamkst
1020 C
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)
1030 C ======================================================================
1031 C END OF THE TEST JOB
1032 C ======================================================================
1033  END
1034  SUBROUTINE filhep(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
1035 C ----------------------------------------------------------------------
1036 C this subroutine fills one entry into the HEPEVT common
1037 C and updates the information for affected mother entries
1038 C
1039 C written by Martin W. Gruenewald (91/01/28)
1040 C
1041 C called by : ZTOHEP,BTOHEP,DWLUxy
1042 C ----------------------------------------------------------------------
1043 C
1044 #include "../../include/HEPEVT.h"
1045 C PARAMETER (NMXHEP=2000)
1046 C COMMON/HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
1047 C &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
1048 C SAVE /HEPEVT/
1049 C COMMON/PHOQED/QEDRAD(NMXHEP)
1050 C LOGICAL QEDRAD
1051 C SAVE /PHOQED/
1052  LOGICAL PHFLAG
1053 C
1054  real*4 p4(4)
1055 C
1056 C check address mode
1057  IF (n.EQ.0) THEN
1058 C
1059 C append mode
1060  ihep=nhep+1
1061  ELSE IF (n.GT.0) THEN
1062 C
1063 C absolute position
1064  ihep=n
1065  ELSE
1066 C
1067 C relative position
1068  ihep=nhep+n
1069  END IF
1070 C
1071 C check on IHEP
1072  IF ((ihep.LE.0).OR.(ihep.GT.nmxhep)) RETURN
1073 C
1074 C 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
1084 C
1085  DO i=1,4
1086  phep(i,ihep)=p4(i)
1087 C
1088 C 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
1092 C FLAG FOR PHOTOS...
1093  qedrad(ihep)=phflag
1094 C
1095 C update process:
1096  DO ip=jmohep(1,ihep),jmohep(2,ihep)
1097  IF(ip.GT.0)THEN
1098 C
1099 C if there is a daughter at IHEP, mother entry at IP has decayed
1100  IF(isthep(ip).EQ.1)isthep(ip)=2
1101 C
1102 C 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
1111 C
1112  RETURN
1113  END
1114 #if defined (ALEPH)
1115  FUNCTION dilogy(X)
1116 C *****************
1117  IMPLICIT REAL*8(a-h,o-z)
1118 CERN 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
1168 C=======================================================================
1169 C===================END OF CPC PART ====================================
1170 C=======================================================================
1171  END
1172 #endif