C++InterfacetoTauola
tauola-factory/Standart_Tauola/tauola.F
1 #if defined (ALEPH)
2 c=============================================================
3 #endif
4  SUBROUTINE jaker(JAK)
5 c *********************
6 c
7 c **********************************************************************
8 c *
9 #if defined (ALEPH)
10 c *********tauola library: version 2.5 ******** *
11 c **************december 1993****************** *
12 #else
13 c *********tauola library: version 2.6 ******** *
14 c **************august 1995****************** *
15 #endif
16 c ** authors: s.jadach, z.was ***** *
17 c ** r. decker, m. jezabek, j.h.kuehn, ***** *
18 c ********available from: wasm at cernvm ****** *
19 c *******published in comp. phys. comm.******** *
20 c *** preprint cern-th-5856 september 1990 **** *
21 c *** preprint cern-th-6195 october 1991 **** *
22 c *** preprint cern-th-6793 november 1992 **** *
23 c **********************************************************************
24 c
25 c ----------------------------------------------------------------------
26 c SUBROUTINE jaker,
27 c chooses decay mode according to list of branching ratios
28 c jak=1 electron mode
29 c jak=2 muon mode
30 c jak=3 pion mode
31 c jak=4 rho mode
32 c jak=5 a1 mode
33 c jak=6 k mode
34 c jak=7 k* mode
35 #if defined (ALEPH)
36 c jak=8-13 npi modes
37 c jak=14-19 kkpi & kpipi modes
38 c jak=20-21 eta pi pi; gamma pi pi modes
39 #else
40 c jak=8 npi mode
41 #endif
42 c
43 c called by : dexay
44 c ----------------------------------------------------------------------
45  COMMON / taubra / gamprt(30),jlist(30),nchan
46 #if defined (ALEPH)
47 #else
48 c REAL cumul(20)
49 #endif
50  REAL cumul(30),rrr(1)
51 c
52  IF(nchan.LE.0.OR.nchan.GT.30) goto 902
53  CALL ranmar(rrr,1)
54  sum=0
55  DO 20 i=1,nchan
56  sum=sum+gamprt(i)
57  20 cumul(i)=sum
58  DO 25 i=nchan,1,-1
59  IF(rrr(1).LT.cumul(i)/cumul(nchan)) ji=i
60  25 CONTINUE
61  jak=jlist(ji)
62  RETURN
63  902 print 9020
64  9020 FORMAT(' ----- JAKER: WRONG NCHAN')
65  stop
66  END
67  SUBROUTINE dekay(KTO,HX)
68 c ***********************
69 c this dekay is in spirit of the 'DECAY' which
70 c was included in koral-b program, comp. phys. commun.
71 c vol. 36 (1985) 191, see comments on general philosophy there.
72 c kto=0 initialisation(obligatory)
73 c kto=1,11 denotes tau+ and kto=2,12 tau-
74 c dekay(1,h) and dekay(2,h) is called internally by mc generator.
75 c h denotes the polarimetric vector, used by the host PROGRAM for
76 c calculation of the spin weight.
77 c user may optionally CALL dekay(11,h) dekay(12,h) in order
78 c to transform decay products to cms and WRITE lund record in /lujets/.
79 c kto=100, print final report(OPTIONAL).
80 c decay modes:
81 c jak=1 electron decay
82 c jak=2 mu decay
83 c jak=3 pi decay
84 c jak=4 rho decay
85 c jak=5 a1 decay
86 c jak=6 k decay
87 c jak=7 k* decay
88 #if defined (ALEPH)
89 c jak= 8-13 npi modes
90 c jak=14-19 kkpi & kpipi modes
91 c jak=20-21 eta pi pi; gamma pi pi modes
92 c jak=0 inclusive: jak=1-21
93 #else
94 c jak=8 npi decay
95 c jak=0 inclusive: jak=1,2,3,4,5,6,7,8
96 #endif
97  REAL h(4)
98  REAL*8 hx(4)
99  COMMON / jaki / jak1,jak2,jakp,jakm,ktom
100 #if defined (ALEPH)
101  COMMON / idfc / idff
102 #else
103  COMMON / idfc / idf
104 #endif
105  COMMON /taupos/ np1,np2
106  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
107  REAL*4 gampmc ,gamper
108  parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
109 #if defined (ALEPH)
110  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
111 #else
112  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
113 #endif
114  & ,names
115  CHARACTER names(nmode)*31
116  COMMON / inout / inut,iout
117  REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4),hdum(4)
118  REAL pdumx(4,9)
119  DATA iwarm/0/
120  ktom=kto
121 #if defined (ALEPH)
122  idf =idff
123 #endif
124  IF(kto.EQ.-1) THEN
125 c ==================
126 c initialisation or reinitialisation
127 c first or second tau positions in hepevt as in koralb/z
128  np1=3
129  np2=4
130  ktom=1
131  IF (iwarm.EQ.1) x=5/(iwarm-1)
132  iwarm=1
133  WRITE(iout,7001) jak1,jak2
134  nevtot=0
135  nev1=0
136  nev2=0
137  IF(jak1.NE.-1.OR.jak2.NE.-1) THEN
138  CALL dadmel(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
139  CALL dadmmu(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
140  CALL dadmpi(-1,idum,pdum,pdum1,pdum2)
141  CALL dadmro(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4)
142  CALL dadmaa(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5,jdum)
143  CALL dadmkk(-1,idum,pdum,pdum1,pdum2)
144  CALL dadmks(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,jdum)
145  CALL dadnew(-1,idum,hdum,pdum1,pdum2,pdumx,jdum)
146  ENDIF
147  DO 21 i=1,30
148  nevdec(i)=0
149  gampmc(i)=0
150  21 gamper(i)=0
151  ELSEIF(kto.EQ.1) THEN
152 c =====================
153 c decay of tau+ in the tau rest frame
154  nevtot=nevtot+1
155  IF(iwarm.EQ.0) goto 902
156  isgn= idf/iabs(idf)
157 #if defined (CePeCe)
158 #elif defined (ALEPH)
159 #else
160 c ajwmod to change brs depending on sign:
161  CALL taurdf(kto)
162 #endif
163  CALL dekay1(0,h,isgn)
164  ELSEIF(kto.EQ.2) THEN
165 c =================================
166 c decay of tau- in the tau rest frame
167  nevtot=nevtot+1
168  IF(iwarm.EQ.0) goto 902
169  isgn=-idf/iabs(idf)
170 #if defined (CePeCe)
171 #elif defined (ALEPH)
172 #else
173 c ajwmod to change brs depending on sign:
174  CALL taurdf(kto)
175 #endif
176  CALL dekay2(0,h,isgn)
177  ELSEIF(kto.EQ.11) THEN
178 c ======================
179 c rest of decay procedure for accepted tau+ decay
180  nev1=nev1+1
181  isgn= idf/iabs(idf)
182  CALL dekay1(1,h,isgn)
183  ELSEIF(kto.EQ.12) THEN
184 c ======================
185 c rest of decay procedure for accepted tau- decay
186  nev2=nev2+1
187  isgn=-idf/iabs(idf)
188  CALL dekay2(1,h,isgn)
189  ELSEIF(kto.EQ.100) THEN
190 c =======================
191  IF(jak1.NE.-1.OR.jak2.NE.-1) THEN
192  CALL dadmel( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
193  CALL dadmmu( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
194  CALL dadmpi( 1,idum,pdum,pdum1,pdum2)
195  CALL dadmro( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4)
196  CALL dadmaa( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5,jdum)
197  CALL dadmkk( 1,idum,pdum,pdum1,pdum2)
198  CALL dadmks( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,jdum)
199  CALL dadnew( 1,idum,hdum,pdum1,pdum2,pdumx,jdum)
200  WRITE(iout,7010) nev1,nev2,nevtot
201  WRITE(iout,7011) (nevdec(i),gampmc(i),gamper(i),i= 1,7)
202  WRITE(iout,7012)
203  $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
204  WRITE(iout,7013)
205  ENDIF
206  ELSE
207 c ====
208  goto 910
209  ENDIF
210 c =====
211  DO 78 k=1,4
212  78 hx(k)=h(k)
213  RETURN
214  7001 FORMAT(///1x,15(5h*****)
215 #if defined (ALEPH)
216  $ /,' *', 25x,'*****TAUOLA LIBRARY: VERSION 2.5 ******',9x,1h*,
217  $ /,' *', 25x,'*DEC 1993; ALEPH fixes introd. dec 98 *',9x,1h*,
218  $ /,' *', 25x,'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
219  $ /,' *', 25x,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
220  $ /,' *', 25x,'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
221  $ /,' *', 25x,'Physics initialization by ALEPH collab ',9x,1h*,
222  $ /,' *', 25x,'it is suggested to use this version ',9x,1h*,
223  $ /,' *', 25x,' with the help of the collab. advice ',9x,1h*,
224  $ /,' *', 25x,'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
225  $ /,' *', 25x,'**5 or more pi dec.: precision limited ',9x,1h*,
226 #elif defined (CLEO)
227  $ /,' *', 25x,'*****TAUOLA LIBRARY: VERSION 2.6 ******',9x,1h*,
228  $ /,' *', 25x,'***********August 1995***************',9x,1h*,
229  $ /,' *', 25x,'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
230  $ /,' *', 25x,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
231  $ /,' *', 25x,'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
232  $ /,' *', 25x,'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
233  $ /,' *', 25x,' Physics initialization by CLEO collab ',9x,1h*,
234  $ /,' *', 25x,' see Alain Weinstein www home page: ',9x,1h*,
235  $ /,' *', 25x,'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
236  $ /,' *', 25x,'/korb_doc.html#files ',9x,1h*,
237  $ /,' *', 25x,'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
238  $ /,' *', 25x,'**5 or more pi dec.: precision limited ',9x,1h*,
239 #else
240  $ /,' *', 25x,'*****TAUOLA LIBRARY: VERSION 2.6 ******',9x,1h*,
241  $ /,' *', 25x,'***********August 1995***************',9x,1h*,
242  $ /,' *', 25x,'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
243  $ /,' *', 25x,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
244  $ /,' *', 25x,'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
245  $ /,' *', 25x,'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
246  $ /,' *', 25x,'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
247  $ /,' *', 25x,'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
248  $ /,' *', 25x,'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
249  $ /,' *', 25x,'**5 or more pi dec.: precision limited ',9x,1h*,
250 #endif
251  $ /,' *', 25x,'****DEKAY ROUTINE: INITIALIZATION******',9x,1h*,
252  $ /,' *',i20 ,5x,'JAK1 = DECAY MODE TAU+ ',9x,1h*,
253  $ /,' *',i20 ,5x,'JAK2 = DECAY MODE TAU- ',9x,1h*,
254  $ /,1x,15(5h*****)/)
255  7010 FORMAT(///1x,15(5h*****)
256 #if defined (ALEPH)
257  $ /,' *', 25x,'*****TAUOLA LIBRARY: VERSION 2.5 ******',9x,1h*,
258  $ /,' *', 25x,'***********DECEMBER 1993***************',9x,1h*,
259 #else
260  $ /,' *', 25x,'*****TAUOLA LIBRARY: VERSION 2.6 ******',9x,1h*,
261  $ /,' *', 25x,'***********August 1995***************',9x,1h*,
262 #endif
263  $ /,' *', 25x,'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
264  $ /,' *', 25x,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
265  $ /,' *', 25x,'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
266  $ /,' *', 25x,'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
267  $ /,' *', 25x,'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
268  $ /,' *', 25x,'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
269  $ /,' *', 25x,'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
270  $ /,' *', 25x,'*****DEKAY ROUTINE: FINAL REPORT*******',9x,1h*,
271  $ /,' *',i20 ,5x,'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*,
272  $ /,' *',i20 ,5x,'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*,
273  $ /,' *',i20 ,5x,'NEVTOT = SUM ',9x,1h*,
274  $ /,' *',' NOEVTS ',
275  $ ' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
276  7011 FORMAT(1x,'*'
277  $ ,i10,2f12.7 ,' DADMEL ELECTRON ',9x,1h*
278  $ /,' *',i10,2f12.7 ,' DADMMU MUON ',9x,1h*
279  $ /,' *',i10,2f12.7 ,' DADMPI PION ',9x,1h*
280  $ /,' *',i10,2f12.7, ' DADMRO RHO (->2PI) ',9x,1h*
281  $ /,' *',i10,2f12.7, ' DADMAA A1 (->3PI) ',9x,1h*
282  $ /,' *',i10,2f12.7, ' DADMKK KAON ',9x,1h*
283  $ /,' *',i10,2f12.7, ' DADMKS K* ',9x,1h*)
284  7012 FORMAT(1x,'*'
285  $ ,i10,2f12.7,a31 ,8x,1h*)
286  7013 FORMAT(1x,'*'
287  $ ,20x,'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
288  $ /,' *',20x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
289  $ /,1x,15(5h*****)/)
290  902 print 9020
291  9020 FORMAT(' ----- DEKAY: LACK OF INITIALISATION')
292  stop
293  910 print 9100
294  9100 FORMAT(' ----- DEKAY: WRONG VALUE OF KTO ')
295  stop
296  END
297  SUBROUTINE dekay1(IMOD,HH,ISGN)
298 c *******************************
299 c this routine simulates tau+ decay
300 #if defined (ALEPH)
301  COMMON / jaki / jak1,jak2,jakp,jakm,ktom
302  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
303  REAL*4 gampmc ,gamper
304  COMMON / decp4 / pp1(4),pp2(4),kff1,kff2
305  REAL*4 pp1 ,pp2
306  INTEGER kff1,kff2
307 #else
308  COMMON / decp4 / pp1(4),pp2(4),kf1,kf2
309  COMMON / jaki / jak1,jak2,jakp,jakm,ktom
310  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
311  REAL*4 gampmc ,gamper
312 #endif
313  REAL hh(4)
314  REAL hv(4),pnu(4),ppi(4)
315  REAL pwb(4),pmu(4),pnm(4)
316  REAL prho(4),pic(4),piz(4)
317  REAL paa(4),pim1(4),pim2(4),pipl(4)
318  REAL pkk(4),pks(4)
319  REAL pnpi(4,9)
320  REAL phot(4)
321  REAL pdum(4)
322  DATA nev,nprin/0,10/
323  kto=1
324  IF(jak1.EQ.-1) RETURN
325  imd=imod
326  IF(imd.EQ.0) THEN
327 c =================
328  jak=jak1
329  IF(jak1.EQ.0) CALL jaker(jak)
330  IF(jak.EQ.1) THEN
331  CALL dadmel(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
332  ELSEIF(jak.EQ.2) THEN
333  CALL dadmmu(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
334  ELSEIF(jak.EQ.3) THEN
335  CALL dadmpi(0, isgn,hv,ppi,pnu)
336  ELSEIF(jak.EQ.4) THEN
337  CALL dadmro(0, isgn,hv,pnu,prho,pic,piz)
338  ELSEIF(jak.EQ.5) THEN
339  CALL dadmaa(0, isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
340  ELSEIF(jak.EQ.6) THEN
341  CALL dadmkk(0, isgn,hv,pkk,pnu)
342  ELSEIF(jak.EQ.7) THEN
343  CALL dadmks(0, isgn,hv,pnu,pks ,pkk,ppi,jkst)
344  ELSE
345  CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
346  ENDIF
347  DO 33 i=1,3
348  33 hh(i)=hv(i)
349  hh(4)=1.0
350 
351  ELSEIF(imd.EQ.1) THEN
352 c =====================
353  nev=nev+1
354  IF (jak.LT.31) THEN
355  nevdec(jak)=nevdec(jak)+1
356  ENDIF
357  DO 34 i=1,4
358  34 pdum(i)=.0
359  IF(jak.EQ.1) THEN
360  CALL dwluel(1,isgn,pnu,pwb,pmu,pnm)
361  CALL dwrph(ktom,phot)
362  DO 10 i=1,4
363  10 pp1(i)=pmu(i)
364 
365  ELSEIF(jak.EQ.2) THEN
366  CALL dwlumu(1,isgn,pnu,pwb,pmu,pnm)
367  CALL dwrph(ktom,phot)
368  DO 20 i=1,4
369  20 pp1(i)=pmu(i)
370 
371  ELSEIF(jak.EQ.3) THEN
372  CALL dwlupi(1,isgn,ppi,pnu)
373  DO 30 i=1,4
374  30 pp1(i)=ppi(i)
375 
376  ELSEIF(jak.EQ.4) THEN
377  CALL dwluro(1,isgn,pnu,prho,pic,piz)
378  DO 40 i=1,4
379  40 pp1(i)=prho(i)
380 
381  ELSEIF(jak.EQ.5) THEN
382  CALL dwluaa(1,isgn,pnu,paa,pim1,pim2,pipl,jaa)
383  DO 50 i=1,4
384  50 pp1(i)=paa(i)
385  ELSEIF(jak.EQ.6) THEN
386  CALL dwlukk(1,isgn,pkk,pnu)
387  DO 60 i=1,4
388  60 pp1(i)=pkk(i)
389  ELSEIF(jak.EQ.7) THEN
390  CALL dwluks(1,isgn,pnu,pks,pkk,ppi,jkst)
391  DO 70 i=1,4
392  70 pp1(i)=pks(i)
393  ELSE
394 cam multipion decay
395  CALL dwlnew(1,isgn,pnu,pwb,pnpi,jak)
396  DO 80 i=1,4
397  80 pp1(i)=pwb(i)
398  ENDIF
399 
400  ENDIF
401 c =====
402  END
403  SUBROUTINE dekay2(IMOD,HH,ISGN)
404 c *******************************
405 c this routine simulates tau- decay
406 #if defined (ALEPH)
407  COMMON / jaki / jak1,jak2,jakp,jakm,ktom
408  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
409  REAL*4 gampmc ,gamper
410  COMMON / decp4 / pp1(4),pp2(4),kff1,kff2
411  REAL*4 pp1 ,pp2
412  INTEGER kff1,kff2
413 #else
414  COMMON / decp4 / pp1(4),pp2(4),kf1,kf2
415  COMMON / jaki / jak1,jak2,jakp,jakm,ktom
416  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
417  REAL*4 gampmc ,gamper
418 #endif
419  REAL hh(4)
420  REAL hv(4),pnu(4),ppi(4)
421  REAL pwb(4),pmu(4),pnm(4)
422  REAL prho(4),pic(4),piz(4)
423  REAL paa(4),pim1(4),pim2(4),pipl(4)
424  REAL pkk(4),pks(4)
425  REAL pnpi(4,9)
426  REAL phot(4)
427  REAL pdum(4)
428  DATA nev,nprin/0,10/
429  kto=2
430  IF(jak2.EQ.-1) RETURN
431  imd=imod
432  IF(imd.EQ.0) THEN
433 c =================
434  jak=jak2
435  IF(jak2.EQ.0) CALL jaker(jak)
436  IF(jak.EQ.1) THEN
437  CALL dadmel(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
438  ELSEIF(jak.EQ.2) THEN
439  CALL dadmmu(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
440  ELSEIF(jak.EQ.3) THEN
441  CALL dadmpi(0, isgn,hv,ppi,pnu)
442  ELSEIF(jak.EQ.4) THEN
443  CALL dadmro(0, isgn,hv,pnu,prho,pic,piz)
444  ELSEIF(jak.EQ.5) THEN
445  CALL dadmaa(0, isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
446  ELSEIF(jak.EQ.6) THEN
447  CALL dadmkk(0, isgn,hv,pkk,pnu)
448  ELSEIF(jak.EQ.7) THEN
449  CALL dadmks(0, isgn,hv,pnu,pks ,pkk,ppi,jkst)
450  ELSE
451  CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
452  ENDIF
453  DO 33 i=1,3
454  33 hh(i)=hv(i)
455  hh(4)=1.0
456  ELSEIF(imd.EQ.1) THEN
457 c =====================
458  nev=nev+1
459  IF (jak.LT.31) THEN
460  nevdec(jak)=nevdec(jak)+1
461  ENDIF
462  DO 34 i=1,4
463  34 pdum(i)=.0
464  IF(jak.EQ.1) THEN
465  CALL dwluel(2,isgn,pnu,pwb,pmu,pnm)
466  CALL dwrph(ktom,phot)
467  DO 10 i=1,4
468  10 pp2(i)=pmu(i)
469 
470  ELSEIF(jak.EQ.2) THEN
471  CALL dwlumu(2,isgn,pnu,pwb,pmu,pnm)
472  CALL dwrph(ktom,phot)
473  DO 20 i=1,4
474  20 pp2(i)=pmu(i)
475 
476  ELSEIF(jak.EQ.3) THEN
477  CALL dwlupi(2,isgn,ppi,pnu)
478  DO 30 i=1,4
479  30 pp2(i)=ppi(i)
480 
481  ELSEIF(jak.EQ.4) THEN
482  CALL dwluro(2,isgn,pnu,prho,pic,piz)
483  DO 40 i=1,4
484  40 pp2(i)=prho(i)
485 
486  ELSEIF(jak.EQ.5) THEN
487  CALL dwluaa(2,isgn,pnu,paa,pim1,pim2,pipl,jaa)
488  DO 50 i=1,4
489  50 pp2(i)=paa(i)
490  ELSEIF(jak.EQ.6) THEN
491  CALL dwlukk(2,isgn,pkk,pnu)
492  DO 60 i=1,4
493  60 pp1(i)=pkk(i)
494  ELSEIF(jak.EQ.7) THEN
495  CALL dwluks(2,isgn,pnu,pks,pkk,ppi,jkst)
496  DO 70 i=1,4
497  70 pp1(i)=pks(i)
498  ELSE
499 cam multipion decay
500  CALL dwlnew(2,isgn,pnu,pwb,pnpi,jak)
501  DO 80 i=1,4
502  80 pp1(i)=pwb(i)
503  ENDIF
504 c
505  ENDIF
506 c =====
507  END
508  SUBROUTINE dexay(KTO,POL)
509 c ----------------------------------------------------------------------
510 c this 'DEXAY' is a routine which generates decay of the single
511 c polarized tau, pol is a polarization vector(not a polarimeter
512 c vector as in dekay) of the tau and it is an input PARAMETER.
513 c kto=0 initialisation(obligatory)
514 c kto=1 denotes tau+ and kto=2 tau-
515 c dexay(1,pol) and dexay(2,pol) are called internally by mc generator.
516 c decay products are transformed readily
517 c to cms and writen in the lund record in /lujets/
518 c kto=100, print final report(OPTIONAL).
519 c
520 c called by : koralz
521 c ----------------------------------------------------------------------
522  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
523  REAL*4 gampmc ,gamper
524  COMMON / jaki / jak1,jak2,jakp,jakm,ktom
525  COMMON / idfc / idff
526  COMMON /taupos/ np1,np2
527  parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
528 #if defined (ALEPH)
529  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
530 #else
531  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
532 #endif
533  & ,names
534  CHARACTER names(nmode)*31
535  COMMON / inout / inut,iout
536  REAL pol(4)
537  REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4)
538  REAL pdum(4)
539  REAL pdumi(4,9)
540  DATA iwarm/0/
541  ktom=kto
542 c
543  IF(kto.EQ.-1) THEN
544 c ==================
545 
546 c initialisation or reinitialisation
547 c first or second tau positions in hepevt as in koralb/z
548  np1=3
549  np2=4
550  iwarm=1
551  WRITE(iout, 7001) jak1,jak2
552  nevtot=0
553  nev1=0
554  nev2=0
555  IF(jak1.NE.-1.OR.jak2.NE.-1) THEN
556  CALL dexel(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
557  CALL dexmu(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
558  CALL dexpi(-1,idum,pdum,pdum1,pdum2)
559  CALL dexro(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4)
560  CALL dexaa(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5,idum)
561  CALL dexkk(-1,idum,pdum,pdum1,pdum2)
562  CALL dexks(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,idum)
563  CALL dexnew(-1,idum,pdum,pdum1,pdum2,pdumi,idum)
564  ENDIF
565  DO 21 i=1,30
566  nevdec(i)=0
567  gampmc(i)=0
568  21 gamper(i)=0
569  ELSEIF(kto.EQ.1) THEN
570 c =====================
571 c decay of tau+ in the tau rest frame
572  nevtot=nevtot+1
573  nev1=nev1+1
574  IF(iwarm.EQ.0) goto 902
575  isgn=idff/iabs(idff)
576 cam CALL dexay1(pol,isgn)
577  CALL dexay1(kto,jak1,jakp,pol,isgn)
578  ELSEIF(kto.EQ.2) THEN
579 c =================================
580 c decay of tau- in the tau rest frame
581  nevtot=nevtot+1
582  nev2=nev2+1
583  IF(iwarm.EQ.0) goto 902
584  isgn=-idff/iabs(idff)
585 cam CALL dexay2(pol,isgn)
586  CALL dexay1(kto,jak2,jakm,pol,isgn)
587  ELSEIF(kto.EQ.100) THEN
588 c =======================
589  IF(jak1.NE.-1.OR.jak2.NE.-1) THEN
590  CALL dexel( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
591  CALL dexmu( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
592  CALL dexpi( 1,idum,pdum,pdum1,pdum2)
593  CALL dexro( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4)
594  CALL dexaa( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5,idum)
595  CALL dexkk( 1,idum,pdum,pdum1,pdum2)
596  CALL dexks( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,idum)
597  CALL dexnew( 1,idum,pdum,pdum1,pdum2,pdumi,idum)
598  WRITE(iout,7010) nev1,nev2,nevtot
599  WRITE(iout,7011) (nevdec(i),gampmc(i),gamper(i),i= 1,7)
600  WRITE(iout,7012)
601  $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
602  WRITE(iout,7013)
603  ENDIF
604  ELSE
605  goto 910
606  ENDIF
607  RETURN
608  7001 FORMAT(///1x,15(5h*****)
609 #if defined (ALEPH)
610  $ /,' *', 25x,'*****TAUOLA LIBRARY: VERSION 2.5 ******',9x,1h*,
611  $ /,' *', 25x,'*DEC 1993; ALEPH fixes introd. dec 98 *',9x,1h*,
612  $ /,' *', 25x,'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
613  $ /,' *', 25x,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
614  $ /,' *', 25x,'Physics initialization by ALEPH collab ',9x,1h*,
615  $ /,' *', 25x,'it is suggested to use this version ',9x,1h*,
616  $ /,' *', 25x,' with the help of the collab. advice ',9x,1h*,
617  $ /,' *', 25x,'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
618  $ /,' *', 25x,'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
619 #elif defined (CLEO)
620  $ /,' *', 25x,'*****TAUOLA LIBRARY: VERSION 2.6 ******',9x,1h*,
621  $ /,' *', 25x,'***********August 1995***************',9x,1h*,
622  $ /,' *', 25x,'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
623  $ /,' *', 25x,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
624  $ /,' *', 25x,'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
625  $ /,' *', 25x,'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
626  $ /,' *', 25x,' Physics initialization by CLEO collab ',9x,1h*,
627  $ /,' *', 25x,' see Alain Weinstein www home page: ',9x,1h*,
628  $ /,' *', 25x,'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
629  $ /,' *', 25x,'/korb_doc.html#files ',9x,1h*,
630  $ /,' *', 25x,'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
631  $ /,' *', 25x,'**5 or more pi dec.: precision limited ',9x,1h*,
632 #else
633  $ /,' *', 25x,'*****TAUOLA LIBRARY: VERSION 2.6 ******',9x,1h*,
634  $ /,' *', 25x,'***********August 1995***************',9x,1h*,
635  $ /,' *', 25x,'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
636  $ /,' *', 25x,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
637  $ /,' *', 25x,'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
638  $ /,' *', 25x,'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
639  $ /,' *', 25x,'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
640  $ /,' *', 25x,'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
641 #endif
642  $ /,' *', 25x,'*******CERN-TH-6793 NOVEMBER 1992*****',9x,1h*,
643  $ /,' *', 25x,'**5 or more pi dec.: precision limited ',9x,1h*,
644  $ /,' *', 25x,'******DEXAY ROUTINE: INITIALIZATION****',9x,1h*
645  $ /,' *',i20 ,5x,'JAK1 = DECAY MODE FERMION1 (TAU+) ',9x,1h*
646  $ /,' *',i20 ,5x,'JAK2 = DECAY MODE FERMION2 (TAU-) ',9x,1h*
647  $ /,1x,15(5h*****)/)
648 chbu format 7010 had more than 19 continuation lines
649 chbu split into two
650  7010 FORMAT(///1x,15(5h*****)
651 #if defined (ALEPH)
652  $ /,' *', 25x,'*****TAUOLA LIBRARY: VERSION 2.5 ******',9x,1h*,
653  $ /,' *', 25x,'***********DECEMBER 1993***************',9x,1h*,
654 #else
655  $ /,' *', 25x,'*****TAUOLA LIBRARY: VERSION 2.6 ******',9x,1h*,
656  $ /,' *', 25x,'***********August 1995***************',9x,1h*,
657 #endif
658  $ /,' *', 25x,'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
659  $ /,' *', 25x,'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
660  $ /,' *', 25x,'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
661  $ /,' *', 25x,'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
662  $ /,' *', 25x,'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
663  $ /,' *', 25x,'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
664  $ /,' *', 25x,'*******CERN-TH-6793 NOVEMBER 1992*****',9x,1h*,
665  $ /,' *', 25x,'******DEXAY ROUTINE: FINAL REPORT******',9x,1h*
666  $ /,' *',i20 ,5x,'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*
667  $ /,' *',i20 ,5x,'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*
668  $ /,' *',i20 ,5x,'NEVTOT = SUM ',9x,1h*
669  $ /,' *',' NOEVTS ',
670  $ ' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
671  7011 FORMAT(1x,'*'
672  $ ,i10,2f12.7 ,' DADMEL ELECTRON ',9x,1h*
673  $ /,' *',i10,2f12.7 ,' DADMMU MUON ',9x,1h*
674  $ /,' *',i10,2f12.7 ,' DADMPI PION ',9x,1h*
675  $ /,' *',i10,2f12.7, ' DADMRO RHO (->2PI) ',9x,1h*
676  $ /,' *',i10,2f12.7, ' DADMAA A1 (->3PI) ',9x,1h*
677  $ /,' *',i10,2f12.7, ' DADMKK KAON ',9x,1h*
678  $ /,' *',i10,2f12.7, ' DADMKS K* ',9x,1h*)
679  7012 FORMAT(1x,'*'
680  $ ,i10,2f12.7,a31 ,8x,1h*)
681  7013 FORMAT(1x,'*'
682  $ ,20x,'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
683  $ /,' *',20x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
684  $ /,1x,15(5h*****)/)
685  902 WRITE(iout, 9020)
686  9020 FORMAT(' ----- DEXAY: LACK OF INITIALISATION')
687  stop
688  910 WRITE(iout, 9100)
689  9100 FORMAT(' ----- DEXAY: WRONG VALUE OF KTO ')
690  stop
691  END
692  SUBROUTINE dexay1(KTO,JAKIN,JAK,POL,ISGN)
693 c ---------------------------------------------------------------------
694 c this routine simulates tau+- decay
695 c
696 c called by : dexay
697 c ---------------------------------------------------------------------
698  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
699  REAL*4 gampmc ,gamper
700  COMMON / inout / inut,iout
701  REAL pol(4),polar(4)
702  REAL pnu(4),ppi(4)
703  REAL prho(4),pic(4),piz(4)
704  REAL pwb(4),pmu(4),pnm(4)
705  REAL paa(4),pim1(4),pim2(4),pipl(4)
706  REAL pkk(4),pks(4)
707  REAL pnpi(4,9)
708  REAL phot(4)
709  REAL pdum(4)
710 c
711  IF(jakin.EQ.-1) RETURN
712  DO 33 i=1,3
713  33 polar(i)=pol(i)
714  polar(4)=0.
715  DO 34 i=1,4
716  34 pdum(i)=.0
717  jak=jakin
718  IF(jak.EQ.0) CALL jaker(jak)
719 cam
720  IF(jak.EQ.1) THEN
721  CALL dexel(0, isgn,polar,pnu,pwb,pmu,pnm,phot)
722  CALL dwluel(kto,isgn,pnu,pwb,pmu,pnm)
723  CALL dwrph(kto,phot )
724  ELSEIF(jak.EQ.2) THEN
725  CALL dexmu(0, isgn,polar,pnu,pwb,pmu,pnm,phot)
726  CALL dwlumu(kto,isgn,pnu,pwb,pmu,pnm)
727  CALL dwrph(kto,phot )
728  ELSEIF(jak.EQ.3) THEN
729  CALL dexpi(0, isgn,polar,ppi,pnu)
730  CALL dwlupi(kto,isgn,ppi,pnu)
731  ELSEIF(jak.EQ.4) THEN
732  CALL dexro(0, isgn,polar,pnu,prho,pic,piz)
733  CALL dwluro(kto,isgn,pnu,prho,pic,piz)
734  ELSEIF(jak.EQ.5) THEN
735  CALL dexaa(0, isgn,polar,pnu,paa,pim1,pim2,pipl,jaa)
736  CALL dwluaa(kto,isgn,pnu,paa,pim1,pim2,pipl,jaa)
737  ELSEIF(jak.EQ.6) THEN
738  CALL dexkk(0, isgn,polar,pkk,pnu)
739  CALL dwlukk(kto,isgn,pkk,pnu)
740  ELSEIF(jak.EQ.7) THEN
741  CALL dexks(0, isgn,polar,pnu,pks,pkk,ppi,jkst)
742  CALL dwluks(kto,isgn,pnu,pks,pkk,ppi,jkst)
743  ELSE
744  jnpi=jak-7
745  CALL dexnew(0, isgn,polar,pnu,pwb,pnpi,jnpi)
746  CALL dwlnew(kto,isgn,pnu,pwb,pnpi,jak)
747  ENDIF
748  nevdec(jak)=nevdec(jak)+1
749  END
750  SUBROUTINE dexel(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
751 c ----------------------------------------------------------------------
752 c this simulates tau decay in tau rest frame
753 c into electron and two neutrinos
754 c
755 c called by : dexay,dexay1
756 c ----------------------------------------------------------------------
757  REAL pol(4),hv(4),pwb(4),pnu(4),q1(4),q2(4),ph(4),rn(1)
758  DATA iwarm/0/
759 c
760  IF(mode.EQ.-1) THEN
761 c ===================
762  iwarm=1
763  CALL dadmel( -1,isgn,hv,pnu,pwb,q1,q2,ph)
764 cc CALL hbook1(813,'WEIGHT DISTRIBUTION DEXEL $',100,0,2)
765 c
766  ELSEIF(mode.EQ. 0) THEN
767 c =======================
768 300 CONTINUE
769  IF(iwarm.EQ.0) goto 902
770  CALL dadmel( 0,isgn,hv,pnu,pwb,q1,q2,ph)
771  wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
772 cc CALL hfill(813,wt)
773  CALL ranmar(rn,1)
774  IF(rn(1).GT.wt) goto 300
775 c
776  ELSEIF(mode.EQ. 1) THEN
777 c =======================
778  CALL dadmel( 1,isgn,hv,pnu,pwb,q1,q2,ph)
779 cc CALL hprint(813)
780  ENDIF
781 c =====
782  RETURN
783  902 print 9020
784  9020 FORMAT(' ----- DEXEL: LACK OF INITIALISATION')
785  stop
786  END
787  SUBROUTINE dexmu(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
788 c ----------------------------------------------------------------------
789 c this simulates tau decay in its rest frame
790 c into muon and two neutrinos
791 c output four momenta: pnu tauneutrino,
792 c pwb w-boson
793 c q1 muon
794 c q2 muon-neutrino
795 c ----------------------------------------------------------------------
796  COMMON / inout / inut,iout
797  REAL pol(4),hv(4),pwb(4),pnu(4),q1(4),q2(4),ph(4),rn(1)
798  DATA iwarm/0/
799 c
800  IF(mode.EQ.-1) THEN
801 c ===================
802  iwarm=1
803  CALL dadmmu( -1,isgn,hv,pnu,pwb,q1,q2,ph)
804 cc CALL hbook1(814,'WEIGHT DISTRIBUTION DEXMU $',100,0,2)
805 c
806  ELSEIF(mode.EQ. 0) THEN
807 c =======================
808 300 CONTINUE
809  IF(iwarm.EQ.0) goto 902
810  CALL dadmmu( 0,isgn,hv,pnu,pwb,q1,q2,ph)
811  wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
812 cc CALL hfill(814,wt)
813  CALL ranmar(rn,1)
814  IF(rn(1).GT.wt) goto 300
815 c
816  ELSEIF(mode.EQ. 1) THEN
817 c =======================
818  CALL dadmmu( 1,isgn,hv,pnu,pwb,q1,q2,ph)
819 cc CALL hprint(814)
820  ENDIF
821 c =====
822  RETURN
823  902 WRITE(iout, 9020)
824  9020 FORMAT(' ----- DEXMU: LACK OF INITIALISATION')
825  stop
826  END
827  SUBROUTINE dadmel(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
828 c ----------------------------------------------------------------------
829 c
830 c called by : dexel,(dekay,dekay1)
831 c ----------------------------------------------------------------------
832  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
833  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
834  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
835  * ,ampiz,ampi,amro,gamro,ama1,gama1
836  * ,amk,amkz,amkst,gamkst
837 c
838  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
839  * ,ampiz,ampi,amro,gamro,ama1,gama1
840  * ,amk,amkz,amkst,gamkst
841  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
842  REAL*4 gampmc ,gamper
843 #if defined (ALEPH)
844 #else
845  REAL*4 phx(4)
846 #endif
847  COMMON / inout / inut,iout
848 #if defined (ALEPH)
849  REAL*4 phx(4)
850 #else
851 #endif
852  REAL hhv(4),hv(4),pwb(4),pnu(4),q1(4),q2(4)
853  REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4)
854  REAL*4 rrr(3)
855  REAL*8 swt, sswt
856  DATA pi /3.141592653589793238462643/
857  DATA iwarm/0/
858 c
859  IF(mode.EQ.-1) THEN
860 c ===================
861  iwarm=1
862  nevraw=0
863  nevacc=0
864  nevovr=0
865  swt=0
866  sswt=0
867  wtmax=1e-20
868  DO 15 i=1,500
869  CALL dphsel(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
870  IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
871 15 CONTINUE
872 cc CALL hbook1(803,'WEIGHT DISTRIBUTION DADMEL $',100,0,2)
873 c
874  ELSEIF(mode.EQ. 0) THEN
875 c =======================
876 300 CONTINUE
877  IF(iwarm.EQ.0) goto 902
878  nevraw=nevraw+1
879  CALL dphsel(wt,hv,pnu,pwb,q1,q2,phx)
880 cc CALL hfill(803,wt/wtmax)
881  swt=swt+wt
882  sswt=sswt+wt**2
883  CALL ranmar(rrr,3)
884  rn=rrr(1)
885  IF(wt.GT.wtmax) nevovr=nevovr+1
886  IF(rn*wtmax.GT.wt) goto 300
887 c rotations to basic tau rest frame
888  rr2=rrr(2)
889  costhe=-1.+2.*rr2
890  thet=acos(costhe)
891  rr3=rrr(3)
892  phi =2*pi*rr3
893  CALL rotor2(thet,pnu,pnu)
894  CALL rotor3( phi,pnu,pnu)
895  CALL rotor2(thet,pwb,pwb)
896  CALL rotor3( phi,pwb,pwb)
897  CALL rotor2(thet,q1,q1)
898  CALL rotor3( phi,q1,q1)
899  CALL rotor2(thet,q2,q2)
900  CALL rotor3( phi,q2,q2)
901  CALL rotor2(thet,hv,hv)
902  CALL rotor3( phi,hv,hv)
903  CALL rotor2(thet,phx,phx)
904  CALL rotor3( phi,phx,phx)
905  DO 44,i=1,3
906  44 hhv(i)=-isgn*hv(i)
907  nevacc=nevacc+1
908 c
909  ELSEIF(mode.EQ. 1) THEN
910 c =======================
911  IF(nevraw.EQ.0) RETURN
912  pargam=swt/float(nevraw+1)
913  error=0
914  IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
915  rat=pargam/gamel
916  WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
917 cc CALL hprint(803)
918  gampmc(1)=rat
919  gamper(1)=error
920 cam nevdec(1)=nevacc
921  ENDIF
922 c =====
923  RETURN
924  7010 FORMAT(///1x,15(5h*****)
925  $ /,' *', 25x,'******** DADMEL FINAL REPORT ******** ',9x,1h*
926  $ /,' *',i20 ,5x,'NEVRAW = NO. OF EL DECAYS TOTAL ',9x,1h*
927  $ /,' *',i20 ,5x,'NEVACC = NO. OF EL DECS. ACCEPTED ',9x,1h*
928  $ /,' *',i20 ,5x,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
929  $ /,' *',e20.5,5x,'PARTIAL WTDTH ( ELECTRON) IN GEV UNITS ',9x,1h*
930  $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
931  $ /,' *',f20.9,5x,'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
932  $ /,' *',25x, 'COMPLETE QED CORRECTIONS INCLUDED ',9x,1h*
933  $ /,' *',25x, 'BUT ONLY V-A CUPLINGS ',9x,1h*
934  $ /,1x,15(5h*****)/)
935  902 WRITE(iout, 9020)
936  9020 FORMAT(' ----- DADMEL: LACK OF INITIALISATION')
937  stop
938  END
939  SUBROUTINE dadmmu(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
940 c ----------------------------------------------------------------------
941  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
942  * ,ampiz,ampi,amro,gamro,ama1,gama1
943  * ,amk,amkz,amkst,gamkst
944 c
945  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
946  * ,ampiz,ampi,amro,gamro,ama1,gama1
947  * ,amk,amkz,amkst,gamkst
948  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
949  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
950  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
951  REAL*4 gampmc ,gamper
952  COMMON / inout / inut,iout
953  REAL*4 phx(4)
954  REAL hhv(4),hv(4),pnu(4),pwb(4),q1(4),q2(4)
955  REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4)
956  REAL*4 rrr(3)
957  REAL*8 swt, sswt
958  DATA pi /3.141592653589793238462643/
959  DATA iwarm /0/
960 c
961  IF(mode.EQ.-1) THEN
962 c ===================
963  iwarm=1
964  nevraw=0
965  nevacc=0
966  nevovr=0
967  swt=0
968  sswt=0
969  wtmax=1e-20
970  DO 15 i=1,500
971  CALL dphsmu(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
972  IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
973 15 CONTINUE
974 cc CALL hbook1(802,'WEIGHT DISTRIBUTION DADMMU $',100,0,2)
975 c
976  ELSEIF(mode.EQ. 0) THEN
977 c =======================
978 300 CONTINUE
979  IF(iwarm.EQ.0) goto 902
980  nevraw=nevraw+1
981  CALL dphsmu(wt,hv,pnu,pwb,q1,q2,phx)
982 cc CALL hfill(802,wt/wtmax)
983  swt=swt+wt
984  sswt=sswt+wt**2
985  CALL ranmar(rrr,3)
986  rn=rrr(1)
987  IF(wt.GT.wtmax) nevovr=nevovr+1
988  IF(rn*wtmax.GT.wt) goto 300
989 c rotations to basic tau rest frame
990  costhe=-1.+2.*rrr(2)
991  thet=acos(costhe)
992  phi =2*pi*rrr(3)
993  CALL rotor2(thet,pnu,pnu)
994  CALL rotor3( phi,pnu,pnu)
995  CALL rotor2(thet,pwb,pwb)
996  CALL rotor3( phi,pwb,pwb)
997  CALL rotor2(thet,q1,q1)
998  CALL rotor3( phi,q1,q1)
999  CALL rotor2(thet,q2,q2)
1000  CALL rotor3( phi,q2,q2)
1001  CALL rotor2(thet,hv,hv)
1002  CALL rotor3( phi,hv,hv)
1003  CALL rotor2(thet,phx,phx)
1004  CALL rotor3( phi,phx,phx)
1005  DO 44,i=1,3
1006  44 hhv(i)=-isgn*hv(i)
1007  nevacc=nevacc+1
1008 c
1009  ELSEIF(mode.EQ. 1) THEN
1010 c =======================
1011  IF(nevraw.EQ.0) RETURN
1012  pargam=swt/float(nevraw+1)
1013  error=0
1014  IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1015  rat=pargam/gamel
1016  WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1017 cc CALL hprint(802)
1018  gampmc(2)=rat
1019  gamper(2)=error
1020 cam nevdec(2)=nevacc
1021  ENDIF
1022 c =====
1023  RETURN
1024  7010 FORMAT(///1x,15(5h*****)
1025  $ /,' *', 25x,'******** DADMMU FINAL REPORT ******** ',9x,1h*
1026  $ /,' *',i20 ,5x,'NEVRAW = NO. OF MU DECAYS TOTAL ',9x,1h*
1027  $ /,' *',i20 ,5x,'NEVACC = NO. OF MU DECS. ACCEPTED ',9x,1h*
1028  $ /,' *',i20 ,5x,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
1029  $ /,' *',e20.5,5x,'PARTIAL WTDTH (MU DECAY) IN GEV UNITS ',9x,1h*
1030  $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1031  $ /,' *',f20.9,5x,'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
1032  $ /,' *',25x, 'COMPLETE QED CORRECTIONS INCLUDED ',9x,1h*
1033  $ /,' *',25x, 'BUT ONLY V-A CUPLINGS ',9x,1h*
1034  $ /,1x,15(5h*****)/)
1035  902 WRITE(iout, 9020)
1036  9020 FORMAT(' ----- DADMMU: LACK OF INITIALISATION')
1037  stop
1038  END
1039  SUBROUTINE dphsel(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
1040 c xnx,xna was flipped in parameters of dphsel and dphsmu
1041 c *********************************************************************
1042 c * electron decay mode *
1043 c *********************************************************************
1044  REAL*4 phx(4)
1045  REAL*4 hvx(4),paax(4),xax(4),qpx(4),xnx(4)
1046  REAL*8 hv(4),ph(4),paa(4),xa(4),qp(4),xn(4)
1047  REAL*8 dgamt
1048  ielmu=1
1049  CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
1050  DO 7 k=1,4
1051  hvx(k)=hv(k)
1052  phx(k)=ph(k)
1053  paax(k)=paa(k)
1054  xax(k)=xa(k)
1055  qpx(k)=qp(k)
1056  xnx(k)=xn(k)
1057  7 CONTINUE
1058  dgamx=dgamt
1059  END
1060  SUBROUTINE dphsmu(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
1061 c xnx,xna was flipped in parameters of dphsel and dphsmu
1062 c *********************************************************************
1063 c * muon decay mode *
1064 c *********************************************************************
1065  REAL*4 phx(4)
1066  REAL*4 hvx(4),paax(4),xax(4),qpx(4),xnx(4)
1067  REAL*8 hv(4),ph(4),paa(4),xa(4),qp(4),xn(4)
1068  REAL*8 dgamt
1069  ielmu=2
1070  CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
1071  DO 7 k=1,4
1072  hvx(k)=hv(k)
1073  phx(k)=ph(k)
1074  paax(k)=paa(k)
1075  xax(k)=xa(k)
1076  qpx(k)=qp(k)
1077  xnx(k)=xn(k)
1078  7 CONTINUE
1079  dgamx=dgamt
1080  END
1081  SUBROUTINE drcmu(DGAMT,HV,PH,PAA,XA,QP,XN,IELMU)
1082  IMPLICIT REAL*8 (a-h,o-z)
1083 c ----------------------------------------------------------------------
1084 * it simulates e,mu channels of tau decay in its rest frame with
1085 * qed order alpha corrections
1086 c ----------------------------------------------------------------------
1087  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1088  * ,ampiz,ampi,amro,gamro,ama1,gama1
1089  * ,amk,amkz,amkst,gamkst
1090 c
1091  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1092  * ,ampiz,ampi,amro,gamro,ama1,gama1
1093  * ,amk,amkz,amkst,gamkst
1094  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1095  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1096 #if defined (ALEPH)
1097  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1098  REAL*4 gampmc ,gamper
1099 #endif
1100  COMMON / inout / inut,iout
1101  COMMON / taurad / xk0dec,itdkrc
1102  REAL*8 xk0dec
1103  REAL*8 hv(4),pt(4),ph(4),paa(4),xa(4),qp(4),xn(4)
1104  REAL*8 pr(4)
1105  REAL*4 rrr(6)
1106  LOGICAL ihard
1107  DATA pi /3.141592653589793238462643d0/
1108 #if defined (CLEO)
1109 c ajwmod to satisfy compiler, comment out this unused function.
1110 #else
1111  xlam(x,y,z)=sqrt((x-y-z)**2-4.0*y*z)
1112 #endif
1113 c amro, gamro is only a PARAMETER for geting hight efficiency
1114 c
1115 c three body phase space normalised as in bjorken-drell
1116 c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
1117  phspac=1./2**17/pi**8
1118  amtax=amtau
1119 c tau momentum
1120  pt(1)=0.d0
1121  pt(2)=0.d0
1122  pt(3)=0.d0
1123  pt(4)=amtax
1124 c
1125  CALL ranmar(rrr,6)
1126 c
1127  IF (ielmu.EQ.1) THEN
1128  amu=amel
1129  ELSE
1130  amu=ammu
1131  ENDIF
1132 c
1133  prhard=0.30d0
1134  IF ( itdkrc.EQ.0) prhard=0d0
1135  prsoft=1.-prhard
1136  IF(prsoft.LT.0.1) THEN
1137  print *, 'ERROR IN DRCMU; PRSOFT=',prsoft
1138  stop
1139  ENDIF
1140 c
1141  rr5=rrr(5)
1142  ihard=(rr5.GT.prsoft)
1143  IF (ihard) THEN
1144 c tau decay to 'TAU+photon'
1145  rr1=rrr(1)
1146  ams1=(amu+amnuta)**2
1147  ams2=(amtax)**2
1148  xk1=1-ams1/ams2
1149  xl1=log(xk1/2/xk0dec)
1150  xl0=log(2*xk0dec)
1151  xk=exp(xl1*rr1+xl0)
1152  am3sq=(1-xk)*ams2
1153  am3 =sqrt(am3sq)
1154  phspac=phspac*ams2*xl1*xk
1155  phspac=phspac/prhard
1156  ELSE
1157  am3=amtax
1158  phspac=phspac*2**6*pi**3
1159  phspac=phspac/prsoft
1160  ENDIF
1161 c mass of neutrina system
1162  rr2=rrr(2)
1163  ams1=(amnuta)**2
1164  ams2=(am3-amu)**2
1165 cam
1166 cam
1167 * flat phase space;
1168  am2sq=ams1+ rr2*(ams2-ams1)
1169  am2 =sqrt(am2sq)
1170  phspac=phspac*(ams2-ams1)
1171 * neutrina rest frame, define xn and xa
1172  enq1=(am2sq+amnuta**2)/(2*am2)
1173  enq2=(am2sq-amnuta**2)/(2*am2)
1174  ppi= enq1**2-amnuta**2
1175  pppi=sqrt(abs(enq1**2-amnuta**2))
1176  phspac=phspac*(4*pi)*(2*pppi/am2)
1177 * nu tau in nunu rest frame
1178  CALL spherd(pppi,xn)
1179  xn(4)=enq1
1180 * nu light in nunu rest frame
1181  DO 30 i=1,3
1182  30 xa(i)=-xn(i)
1183  xa(4)=enq2
1184 * tau-prim rest frame, define qp(muon
1185 * nunu momentum
1186  pr(1)=0
1187  pr(2)=0
1188  pr(4)=1.d0/(2*am3)*(am3**2+am2**2-amu**2)
1189  pr(3)= sqrt(abs(pr(4)**2-am2**2))
1190  ppi = pr(4)**2-am2**2
1191 * muon momentum
1192  qp(1)=0
1193  qp(2)=0
1194  qp(4)=1.d0/(2*am3)*(am3**2-am2**2+amu**2)
1195  qp(3)=-pr(3)
1196  phspac=phspac*(4*pi)*(2*pr(3)/am3)
1197 * neutrina boosted from their frame to tau-prim rest frame
1198  exe=(pr(4)+pr(3))/am2
1199  CALL bostd3(exe,xn,xn)
1200  CALL bostd3(exe,xa,xa)
1201  rr3=rrr(3)
1202  rr4=rrr(4)
1203  IF (ihard) THEN
1204  eps=4*(amu/amtax)**2
1205  xl1=log((2+eps)/eps)
1206  xl0=log(eps)
1207  eta =exp(xl1*rr3+xl0)
1208  cthet=1+eps-eta
1209  thet =acos(cthet)
1210  phspac=phspac*xl1/2*eta
1211  phi = 2*pi*rr4
1212  CALL rotpox(thet,phi,xn)
1213  CALL rotpox(thet,phi,xa)
1214  CALL rotpox(thet,phi,qp)
1215  CALL rotpox(thet,phi,pr)
1216 c
1217 * now to the tau rest frame, define tau-prim and gamma momenta
1218 * tau-prim momentum
1219  paa(1)=0
1220  paa(2)=0
1221  paa(4)=1/(2*amtax)*(amtax**2+am3**2)
1222  paa(3)= sqrt(abs(paa(4)**2-am3**2))
1223  ppi = paa(4)**2-am3**2
1224  phspac=phspac*(4*pi)*(2*paa(3)/amtax)
1225 * gamma momentum
1226  ph(1)=0
1227  ph(2)=0
1228  ph(4)=paa(3)
1229  ph(3)=-paa(3)
1230 * all momenta boosted from tau-prim rest frame to tau rest frame
1231 * z-axis antiparallel to photon momentum
1232  exe=(paa(4)+paa(3))/am3
1233  CALL bostd3(exe,xn,xn)
1234  CALL bostd3(exe,xa,xa)
1235  CALL bostd3(exe,qp,qp)
1236  CALL bostd3(exe,pr,pr)
1237  ELSE
1238  thet =acos(-1.+2*rr3)
1239  phi = 2*pi*rr4
1240  CALL rotpox(thet,phi,xn)
1241  CALL rotpox(thet,phi,xa)
1242  CALL rotpox(thet,phi,qp)
1243  CALL rotpox(thet,phi,pr)
1244 c
1245 * now to the tau rest frame, define tau-prim and gamma momenta
1246 * tau-prim momentum
1247  paa(1)=0
1248  paa(2)=0
1249  paa(4)=amtax
1250  paa(3)=0
1251 * gamma momentum
1252  ph(1)=0
1253  ph(2)=0
1254  ph(4)=0
1255  ph(3)=0
1256  ENDIF
1257 c partial width consists of phase space and amplitude
1258  CALL dampry(itdkrc,xk0dec,ph,xa,qp,xn,amplit,hv)
1259  dgamt=1/(2.*amtax)*amplit*phspac
1260  END
1261  SUBROUTINE dampry(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1262  IMPLICIT REAL*8 (a-h,o-z)
1263 c ----------------------------------------------------------------------
1264 c it calculates matrix element for the
1265 c tau --> mu(e) nu nubar decay mode
1266 c including complete order alpha qed corrections.
1267 c ----------------------------------------------------------------------
1268  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1269  * ,ampiz,ampi,amro,gamro,ama1,gama1
1270  * ,amk,amkz,amkst,gamkst
1271 c
1272  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1273  * ,ampiz,ampi,amro,gamro,ama1,gama1
1274  * ,amk,amkz,amkst,gamkst
1275  REAL*8 hv(4),qp(4),xn(4),xa(4),xk(4)
1276 c
1277  hv(4)=1.d0
1278  ak0=xk0dec*amtau
1279  IF(xk(4).LT.0.1d0*ak0) THEN
1280  amplit=thb(itdkrc,qp,xn,xa,ak0,hv)
1281  ELSE
1282  amplit=sqm2(itdkrc,qp,xn,xa,xk,ak0,hv)
1283  ENDIF
1284  RETURN
1285  END
1286  FUNCTION sqm2(ITDKRC,QP,XN,XA,XK,AK0,HV)
1287 c
1288 c **********************************************************************
1289 c REAL photon matrix element squared *
1290 c parameters: *
1291 c hv- polarimetric four-vector of tau *
1292 c qp,xn,xa,xk - 4-momenta of electron(muon), nu, nubar and photon *
1293 c all four-vectors in tau rest frame(in gev) *
1294 c ak0 - infrared cutoff, minimal energy of hard photons(gev) *
1295 c sqm2 - value for s=0 *
1296 c see eqs. (2.9)-(2.10) from cjk( nucl.phys.b(1991) ) *
1297 c **********************************************************************
1298 c
1299  IMPLICIT REAL*8(a-h,o-z)
1300  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1301  * ,ampiz,ampi,amro,gamro,ama1,gama1
1302  * ,amk,amkz,amkst,gamkst
1303 c
1304  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1305  * ,ampiz,ampi,amro,gamro,ama1,gama1
1306  * ,amk,amkz,amkst,gamkst
1307  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1308  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1309  COMMON / qedprm /alfinv,alfpi,xk0
1310  REAL*8 alfinv,alfpi,xk0
1311  REAL*8 qp(4),xn(4),xa(4),xk(4)
1312  REAL*8 r(4)
1313  REAL*8 hv(4)
1314  REAL*8 s0(3),rxa(3),rxk(3),rqp(3)
1315  DATA pi /3.141592653589793238462643d0/
1316 c
1317  tmass=amtau
1318  gf=gfermi
1319  alphai=alfinv
1320  tmass2=tmass**2
1321  emass2=qp(4)**2-qp(1)**2-qp(2)**2-qp(3)**2
1322  r(4)=tmass
1323 c scalar products of four-momenta
1324  DO 7 i=1,3
1325  r(1)=0.d0
1326  r(2)=0.d0
1327  r(3)=0.d0
1328  r(i)=tmass
1329  rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1330 c rxn(i)=r(4)*xn(4)-r(1)*xn(1)-r(2)*xn(2)-r(3)*xn(3)
1331  rxk(i)=r(4)*xk(4)-r(1)*xk(1)-r(2)*xk(2)-r(3)*xk(3)
1332  rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
1333  7 CONTINUE
1334  qpxn=qp(4)*xn(4)-qp(1)*xn(1)-qp(2)*xn(2)-qp(3)*xn(3)
1335  qpxa=qp(4)*xa(4)-qp(1)*xa(1)-qp(2)*xa(2)-qp(3)*xa(3)
1336  qpxk=qp(4)*xk(4)-qp(1)*xk(1)-qp(2)*xk(2)-qp(3)*xk(3)
1337 c xnxa=xn(4)*xa(4)-xn(1)*xa(1)-xn(2)*xa(2)-xn(3)*xa(3)
1338  xnxk=xn(4)*xk(4)-xn(1)*xk(1)-xn(2)*xk(2)-xn(3)*xk(3)
1339  xaxk=xa(4)*xk(4)-xa(1)*xk(1)-xa(2)*xk(2)-xa(3)*xk(3)
1340  txn=tmass*xn(4)
1341  txa=tmass*xa(4)
1342  tqp=tmass*qp(4)
1343  txk=tmass*xk(4)
1344 c
1345  x= xnxk/qpxn
1346  z= txk/tqp
1347  a= 1+x
1348  b= 1+ x*(1+z)/2+z/2
1349  s1= qpxn*txa*( -emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
1350  $tmass2/txk**2) +
1351  $qpxn/txk**2* ( tmass2*xaxk - txa*txk+ xaxk*txk) -
1352  $txa*txn/txk - qpxn/(qpxk*txk)* (tqp*xaxk-txk*qpxa)
1353  const4=256*pi/alphai*gf**2
1354  IF (itdkrc.EQ.0) const4=0d0
1355  sqm2=s1*const4
1356  DO 5 i=1,3
1357  s0(i) = qpxn*rxa(i)*(-emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
1358  $ tmass2/txk**2) +
1359  $ qpxn/txk**2* (tmass2*xaxk - txa*rxk(i)+ xaxk*rxk(i))-
1360  $ rxa(i)*txn/txk - qpxn/(qpxk*txk)*(rqp(i)*xaxk- rxk(i)*qpxa)
1361  5 hv(i)=s0(i)/s1-1.d0
1362  RETURN
1363  END
1364  FUNCTION thb(ITDKRC,QP,XN,XA,AK0,HV)
1365 c
1366 c **********************************************************************
1367 c born +virtual+soft photon matrix element**2 o(alpha) *
1368 c parameters: *
1369 c hv- polarimetric four-vector of tau *
1370 c qp,xn,xa - four-momenta of electron(muon), nu and nubar in gev *
1371 c all four-vectors in tau rest frame *
1372 c ak0 - infrared cutoff, minimal energy of hard photons *
1373 c thb - value for s=0 *
1374 c see eqs. (2.2),(2.4)-(2.5) from cjk(nucl.phys.b351(1991)70 *
1375 c and(c.2) from jk(nucl.phys.b320(1991)20 ) *
1376 c **********************************************************************
1377 c
1378  IMPLICIT REAL*8(a-h,o-z)
1379  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1380  * ,ampiz,ampi,amro,gamro,ama1,gama1
1381  * ,amk,amkz,amkst,gamkst
1382 c
1383  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1384  * ,ampiz,ampi,amro,gamro,ama1,gama1
1385  * ,amk,amkz,amkst,gamkst
1386  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1387  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1388  COMMON / qedprm /alfinv,alfpi,xk0
1389  REAL*8 alfinv,alfpi,xk0
1390  dimension qp(4),xn(4),xa(4)
1391  REAL*8 hv(4)
1392  dimension r(4)
1393  REAL*8 rxa(3),rxn(3),rqp(3)
1394  REAL*8 bornpl(3),am3pol(3),xm3pol(3)
1395  DATA pi /3.141592653589793238462643d0/
1396 c
1397  tmass=amtau
1398  gf=gfermi
1399  alphai=alfinv
1400 c
1401  tmass2=tmass**2
1402  r(4)=tmass
1403  DO 7 i=1,3
1404  r(1)=0.d0
1405  r(2)=0.d0
1406  r(3)=0.d0
1407  r(i)=tmass
1408  rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1409  rxn(i)=r(4)*xn(4)-r(1)*xn(1)-r(2)*xn(2)-r(3)*xn(3)
1410 c rxk(i)=r(4)*xk(4)-r(1)*xk(1)-r(2)*xk(2)-r(3)*xk(3)
1411  rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
1412  7 CONTINUE
1413 c quasi two-body variables
1414  u0=qp(4)/tmass
1415  u3=sqrt(qp(1)**2+qp(2)**2+qp(3)**2)/tmass
1416  w3=u3
1417  w0=(xn(4)+xa(4))/tmass
1418  up=u0+u3
1419  um=u0-u3
1420  wp=w0+w3
1421  wm=w0-w3
1422  yu=log(up/um)/2
1423  yw=log(wp/wm)/2
1424  eps2=u0**2-u3**2
1425  eps=sqrt(eps2)
1426  y=w0**2-w3**2
1427  al=ak0/tmass
1428 c formfactors
1429  f0=2*u0/u3*( dilogt(1-(um*wm/(up*wp)))- dilogt(1-wm/wp) +
1430  $dilogt(1-um/up) -2*yu+ 2*log(up)*(yw+yu) ) +
1431  $1/y* ( 2*u3*yu + (1-eps2- 2*y)*log(eps) ) +
1432  $ 2 - 4*(u0/u3*yu -1)* log(2*al)
1433  fp= yu/(2*u3)*(1 + (1-eps2)/y ) + log(eps)/y
1434  fm= yu/(2*u3)*(1 - (1-eps2)/y ) - log(eps)/y
1435  f3= eps2*(fp+fm)/2
1436 c scalar products of four-momenta
1437  qpxn=qp(4)*xn(4)-qp(1)*xn(1)-qp(2)*xn(2)-qp(3)*xn(3)
1438  qpxa=qp(4)*xa(4)-qp(1)*xa(1)-qp(2)*xa(2)-qp(3)*xa(3)
1439  xnxa=xn(4)*xa(4)-xn(1)*xa(1)-xn(2)*xa(2)-xn(3)*xa(3)
1440  txn=tmass*xn(4)
1441  txa=tmass*xa(4)
1442  tqp=tmass*qp(4)
1443 c decay differential width without and with polarization
1444  const3=1/(2*alphai*pi)*64*gf**2
1445  IF (itdkrc.EQ.0) const3=0d0
1446  xm3= -( f0* qpxn*txa + fp*eps2* txn*txa +
1447  $fm* qpxn*qpxa + f3* tmass2*xnxa )
1448  am3=xm3*const3
1449 c v-a and v+a couplings, but in the born part only
1450  brak= (gv+ga)**2*tqp*xnxa+(gv-ga)**2*txa*qpxn
1451  & -(gv**2-ga**2)*tmass*amnuta*qpxa
1452  born= 32*(gfermi**2/2.)*brak
1453  DO 5 i=1,3
1454  xm3pol(i)= -( f0* qpxn*rxa(i) + fp*eps2* txn*rxa(i) +
1455  $ fm* qpxn* (qpxa + (rxa(i)*tqp-txa*rqp(i))/tmass2 ) +
1456  $ f3* (tmass2*xnxa +txn*rxa(i) -rxn(i)*txa) )
1457  am3pol(i)=xm3pol(i)*const3
1458 c v-a and v+a couplings, but in the born part only
1459  bornpl(i)=born+(
1460  & (gv+ga)**2*tmass*xnxa*qp(i)
1461  & -(gv-ga)**2*tmass*qpxn*xa(i)
1462  & +(gv**2-ga**2)*amnuta*txa*qp(i)
1463  & -(gv**2-ga**2)*amnuta*tqp*xa(i) )*
1464  & 32*(gfermi**2/2.)
1465  5 hv(i)=(bornpl(i)+am3pol(i))/(born+am3)-1.d0
1466  thb=born+am3
1467  IF (thb/born.LT.0.1d0) THEN
1468  print *, 'ERROR IN THB, THB/BORN=',thb/born
1469 #if defined (CLEO)
1470  thb=0.d0
1471 #else
1472  stop
1473 #endif
1474  ENDIF
1475  RETURN
1476  END
1477  SUBROUTINE dexpi(MODE,ISGN,POL,PPI,PNU)
1478 c ----------------------------------------------------------------------
1479 c tau decay into pion and tau-neutrino
1480 c in tau rest frame
1481 c output four momenta: pnu tauneutrino,
1482 c ppi pion charged
1483 c ----------------------------------------------------------------------
1484  REAL pol(4),hv(4),pnu(4),ppi(4),rn(1)
1485 cc
1486  IF(mode.EQ.-1) THEN
1487 c ===================
1488  CALL dadmpi(-1,isgn,hv,ppi,pnu)
1489 cc CALL hbook1(815,'WEIGHT DISTRIBUTION DEXPI $',100,0,2)
1490 
1491  ELSEIF(mode.EQ. 0) THEN
1492 c =======================
1493 300 CONTINUE
1494  CALL dadmpi( 0,isgn,hv,ppi,pnu)
1495  wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1496 cc CALL hfill(815,wt)
1497  CALL ranmar(rn,1)
1498  IF(rn(1).GT.wt) goto 300
1499 c
1500  ELSEIF(mode.EQ. 1) THEN
1501 c =======================
1502  CALL dadmpi( 1,isgn,hv,ppi,pnu)
1503 cc CALL hprint(815)
1504  ENDIF
1505 c =====
1506  RETURN
1507  END
1508  SUBROUTINE dadmpi(MODE,ISGN,HV,PPI,PNU)
1509 c ----------------------------------------------------------------------
1510  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1511  * ,ampiz,ampi,amro,gamro,ama1,gama1
1512  * ,amk,amkz,amkst,gamkst
1513 c
1514  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1515  * ,ampiz,ampi,amro,gamro,ama1,gama1
1516  * ,amk,amkz,amkst,gamkst
1517  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1518  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1519  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1520  REAL*4 gampmc ,gamper
1521  COMMON / inout / inut,iout
1522  REAL ppi(4),pnu(4),hv(4)
1523  DATA pi /3.141592653589793238462643/
1524 c
1525  IF(mode.EQ.-1) THEN
1526 c ===================
1527  nevtot=0
1528  ELSEIF(mode.EQ. 0) THEN
1529 c =======================
1530  nevtot=nevtot+1
1531  epi= (amtau**2+ampi**2-amnuta**2)/(2*amtau)
1532  enu= (amtau**2-ampi**2+amnuta**2)/(2*amtau)
1533  xpi= sqrt(epi**2-ampi**2)
1534 c pi momentum
1535  CALL sphera(xpi,ppi)
1536  ppi(4)=epi
1537 c tau-neutrino momentum
1538  DO 30 i=1,3
1539 30 pnu(i)=-ppi(i)
1540  pnu(4)=enu
1541  pxq=amtau*epi
1542  pxn=amtau*enu
1543  qxn=ppi(4)*pnu(4)-ppi(1)*pnu(1)-ppi(2)*pnu(2)-ppi(3)*pnu(3)
1544  brak=(gv**2+ga**2)*(2*pxq*qxn-ampi**2*pxn)
1545  & +(gv**2-ga**2)*amtau*amnuta*ampi**2
1546  DO 40 i=1,3
1547 40 hv(i)=-isgn*2*ga*gv*amtau*(2*ppi(i)*qxn-pnu(i)*ampi**2)/brak
1548  hv(4)=1
1549 c
1550  ELSEIF(mode.EQ. 1) THEN
1551 c =======================
1552  IF(nevtot.EQ.0) RETURN
1553  fpi=0.1284
1554 c gamm=(gfermi*fpi)**2/(16.*pi)*amtau**3*
1555 c * (brak/amtau**4)**2
1556 czw 7.02.93 here was an error affecting non standard model
1557 c configurations only
1558  gamm=(gfermi*fpi)**2/(16.*pi)*amtau**3*
1559  $ (brak/amtau**4)*
1560  $ sqrt((amtau**2-ampi**2-amnuta**2)**2
1561  $ -4*ampi**2*amnuta**2 )/amtau**2
1562  error=0
1563  rat=gamm/gamel
1564  WRITE(iout, 7010) nevtot,gamm,rat,error
1565  gampmc(3)=rat
1566  gamper(3)=error
1567 cam nevdec(3)=nevtot
1568  ENDIF
1569 c =====
1570  RETURN
1571  7010 FORMAT(///1x,15(5h*****)
1572  $ /,' *', 25x,'******** DADMPI FINAL REPORT ******** ',9x,1h*
1573  $ /,' *',i20 ,5x,'NEVTOT = NO. OF PI DECAYS TOTAL ',9x,1h*
1574  $ /,' *',e20.5,5x,'PARTIAL WTDTH ( PI DECAY) IN GEV UNITS ',9x,1h*
1575  $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1576  $ /,' *',f20.8,5x,'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9x,1h*
1577  $ /,1x,15(5h*****)/)
1578  END
1579  SUBROUTINE dexro(MODE,ISGN,POL,PNU,PRO,PIC,PIZ)
1580 c ----------------------------------------------------------------------
1581 c this simulates tau decay in tau rest frame
1582 c into nu rho, next rho decays into pion pair.
1583 c output four momenta: pnu tauneutrino,
1584 c pro rho
1585 c pic pion charged
1586 c piz pion zero
1587 c ----------------------------------------------------------------------
1588  COMMON / inout / inut,iout
1589  REAL pol(4),hv(4),pro(4),pnu(4),pic(4),piz(4),rn(1)
1590  DATA iwarm/0/
1591 c
1592  IF(mode.EQ.-1) THEN
1593 c ===================
1594  iwarm=1
1595  CALL dadmro( -1,isgn,hv,pnu,pro,pic,piz)
1596 cc CALL hbook1(816,'WEIGHT DISTRIBUTION DEXRO $',100,0,2)
1597 cc CALL hbook1(916,'ABS2 OF HV IN ROUTINE DEXRO $',100,0,2)
1598 c
1599  ELSEIF(mode.EQ. 0) THEN
1600 c =======================
1601 300 CONTINUE
1602  IF(iwarm.EQ.0) goto 902
1603  CALL dadmro( 0,isgn,hv,pnu,pro,pic,piz)
1604  wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1605 cc CALL hfill(816,wt)
1606 cc xhelp=hv(1)**2+hv(2)**2+hv(3)**2
1607 cc CALL hfill(916,xhelp)
1608  CALL ranmar(rn,1)
1609  IF(rn(1).GT.wt) goto 300
1610 c
1611  ELSEIF(mode.EQ. 1) THEN
1612 c =======================
1613  CALL dadmro( 1,isgn,hv,pnu,pro,pic,piz)
1614 cc CALL hprint(816)
1615 cc CALL hprint(916)
1616  ENDIF
1617 c =====
1618  RETURN
1619  902 WRITE(iout, 9020)
1620  9020 FORMAT(' ----- DEXRO: LACK OF INITIALISATION')
1621  stop
1622  END
1623  SUBROUTINE dadmro(MODE,ISGN,HHV,PNU,PRO,PIC,PIZ)
1624 c ----------------------------------------------------------------------
1625  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1626  * ,ampiz,ampi,amro,gamro,ama1,gama1
1627  * ,amk,amkz,amkst,gamkst
1628 c
1629  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1630  * ,ampiz,ampi,amro,gamro,ama1,gama1
1631  * ,amk,amkz,amkst,gamkst
1632  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1633  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1634  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1635  REAL*4 gampmc ,gamper
1636  COMMON / inout / inut,iout
1637  REAL hhv(4)
1638  REAL hv(4),pro(4),pnu(4),pic(4),piz(4)
1639  REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4)
1640  REAL*4 rrr(3)
1641  REAL*8 swt, sswt
1642  DATA pi /3.141592653589793238462643/
1643  DATA iwarm/0/
1644 c
1645  IF(mode.EQ.-1) THEN
1646 c ===================
1647  iwarm=1
1648  nevraw=0
1649  nevacc=0
1650  nevovr=0
1651  swt=0
1652  sswt=0
1653  wtmax=1e-20
1654  DO 15 i=1,500
1655  CALL dphsro(wt,hv,pdum1,pdum2,pdum3,pdum4)
1656  IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1657 15 CONTINUE
1658 cc CALL hbook1(801,'WEIGHT DISTRIBUTION DADMRO $',100,0,2)
1659 cc print 7003,wtmax
1660 c
1661  ELSEIF(mode.EQ. 0) THEN
1662 c =======================
1663 300 CONTINUE
1664  IF(iwarm.EQ.0) goto 902
1665  CALL dphsro(wt,hv,pnu,pro,pic,piz)
1666 cc CALL hfill(801,wt/wtmax)
1667  nevraw=nevraw+1
1668  swt=swt+wt
1669  sswt=sswt+wt**2
1670  CALL ranmar(rrr,3)
1671  rn=rrr(1)
1672  IF(wt.GT.wtmax) nevovr=nevovr+1
1673  IF(rn*wtmax.GT.wt) goto 300
1674 c rotations to basic tau rest frame
1675  costhe=-1.+2.*rrr(2)
1676  thet=acos(costhe)
1677  phi =2*pi*rrr(3)
1678  CALL rotor2(thet,pnu,pnu)
1679  CALL rotor3( phi,pnu,pnu)
1680  CALL rotor2(thet,pro,pro)
1681  CALL rotor3( phi,pro,pro)
1682  CALL rotor2(thet,pic,pic)
1683  CALL rotor3( phi,pic,pic)
1684  CALL rotor2(thet,piz,piz)
1685  CALL rotor3( phi,piz,piz)
1686  CALL rotor2(thet,hv,hv)
1687  CALL rotor3( phi,hv,hv)
1688  DO 44 i=1,3
1689  44 hhv(i)=-isgn*hv(i)
1690  nevacc=nevacc+1
1691 c
1692  ELSEIF(mode.EQ. 1) THEN
1693 c =======================
1694  IF(nevraw.EQ.0) RETURN
1695  pargam=swt/float(nevraw+1)
1696  error=0
1697  IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1698  rat=pargam/gamel
1699  WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1700 cc CALL hprint(801)
1701  gampmc(4)=rat
1702  gamper(4)=error
1703 cam nevdec(4)=nevacc
1704  ENDIF
1705 c =====
1706  RETURN
1707  7003 FORMAT(///1x,15(5h*****)
1708  $ /,' *', 25x,'******** DADMRO INITIALISATION ********',9x,1h*
1709  $ /,' *',e20.5,5x,'WTMAX = MAXIMUM WEIGHT ',9x,1h*
1710  $ /,1x,15(5h*****)/)
1711  7010 FORMAT(///1x,15(5h*****)
1712  $ /,' *', 25x,'******** DADMRO FINAL REPORT ******** ',9x,1h*
1713  $ /,' *',i20 ,5x,'NEVRAW = NO. OF RHO DECAYS TOTAL ',9x,1h*
1714  $ /,' *',i20 ,5x,'NEVACC = NO. OF RHO DECS. ACCEPTED ',9x,1h*
1715  $ /,' *',i20 ,5x,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
1716  $ /,' *',e20.5,5x,'PARTIAL WTDTH (RHO DECAY) IN GEV UNITS ',9x,1h*
1717  $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1718  $ /,' *',f20.8,5x,'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
1719  $ /,1x,15(5h*****)/)
1720  902 WRITE(iout, 9020)
1721  9020 FORMAT(' ----- DADMRO: LACK OF INITIALISATION')
1722  stop
1723  END
1724  SUBROUTINE dphsro(DGAMT,HV,PN,PR,PIC,PIZ)
1725 c ----------------------------------------------------------------------
1726 c it simulates rho decay in tau rest frame with
1727 c z-axis along rho momentum
1728 c ----------------------------------------------------------------------
1729  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1730  * ,ampiz,ampi,amro,gamro,ama1,gama1
1731  * ,amk,amkz,amkst,gamkst
1732 c
1733  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1734  * ,ampiz,ampi,amro,gamro,ama1,gama1
1735  * ,amk,amkz,amkst,gamkst
1736  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1737  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1738  REAL hv(4),pt(4),pn(4),pr(4),pic(4),piz(4),qq(4),rr1(1)
1739  DATA pi /3.141592653589793238462643/
1740  DATA icont /0/
1741 c
1742 c three body phase space normalised as in bjorken-drell
1743  phspac=1./2**11/pi**5
1744 c tau momentum
1745  pt(1)=0.
1746  pt(2)=0.
1747  pt(3)=0.
1748  pt(4)=amtau
1749 c mass of(real/virtual) rho
1750  ams1=(ampi+ampiz)**2
1751  ams2=(amtau-amnuta)**2
1752 c flat phase space
1753 #if defined (ALEPH)
1754 c amx2=ams1+ rr1(1)*(ams2-ams1)
1755 #else
1756 c amx2=ams1+ rr1*(ams2-ams1)
1757 #endif
1758 c amx=sqrt(amx2)
1759 c phspac=phspac*(ams2-ams1)
1760 c phase space with sampling for rho resonance
1761  alp1=atan((ams1-amro**2)/amro/gamro)
1762  alp2=atan((ams2-amro**2)/amro/gamro)
1763 cam
1764  100 CONTINUE
1765  CALL ranmar(rr1,1)
1766  alp=alp1+rr1(1)*(alp2-alp1)
1767  amx2=amro**2+amro*gamro*tan(alp)
1768  amx=sqrt(amx2)
1769  IF(amx.LT.2.*ampi) go to 100
1770 cam
1771  phspac=phspac*((amx2-amro**2)**2+(amro*gamro)**2)/(amro*gamro)
1772  phspac=phspac*(alp2-alp1)
1773 c
1774 c tau-neutrino momentum
1775  pn(1)=0
1776  pn(2)=0
1777  pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
1778  pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
1779 c rho momentum
1780  pr(1)=0
1781  pr(2)=0
1782  pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
1783  pr(3)=-pn(3)
1784  phspac=phspac*(4*pi)*(2*pr(3)/amtau)
1785 c
1786 cam
1787  enq1=(amx2+ampi**2-ampiz**2)/(2.*amx)
1788  enq2=(amx2-ampi**2+ampiz**2)/(2.*amx)
1789  pppi=sqrt((enq1-ampi)*(enq1+ampi))
1790  phspac=phspac*(4*pi)*(2*pppi/amx)
1791 c charged pi momentum in rho rest frame
1792  CALL sphera(pppi,pic)
1793  pic(4)=enq1
1794 c neutral pi momentum in rho rest frame
1795  DO 20 i=1,3
1796 20 piz(i)=-pic(i)
1797  piz(4)=enq2
1798  exe=(pr(4)+pr(3))/amx
1799 c pions boosted from rho rest frame to tau rest frame
1800  CALL bostr3(exe,pic,pic)
1801  CALL bostr3(exe,piz,piz)
1802  DO 30 i=1,4
1803 30 qq(i)=pic(i)-piz(i)
1804 c amplitude
1805  prodpq=pt(4)*qq(4)
1806  prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
1807  prodpn=pt(4)*pn(4)
1808  qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
1809  brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
1810  & +(gv**2-ga**2)*amtau*amnuta*qq2
1811  amplit=(gfermi*ccabib)**2*brak*2*fpirho(amx)
1812  dgamt=1/(2.*amtau)*amplit*phspac
1813  DO 40 i=1,3
1814  40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
1815  RETURN
1816  END
1817  SUBROUTINE dexaa(MODE,ISGN,POL,PNU,PAA,PIM1,PIM2,PIPL,JAA)
1818 c ----------------------------------------------------------------------
1819 * this simulates tau decay in tau rest frame
1820 * into nu a1, next a1 decays into rho pi and finally rho into pi pi.
1821 * output four momenta: pnu tauneutrino,
1822 * paa a1
1823 * pim1 pion minus(or pi0) 1 (for tau minus)
1824 * pim2 pion minus(or pi0) 2
1825 * pipl pion plus(or pi-)
1826 * (pipl,pim1) form a rho
1827 c ----------------------------------------------------------------------
1828  COMMON / inout / inut,iout
1829  REAL pol(4),hv(4),paa(4),pnu(4),pim1(4),pim2(4),pipl(4),rn(1)
1830  DATA iwarm/0/
1831 c
1832  IF(mode.EQ.-1) THEN
1833 c ===================
1834  iwarm=1
1835  CALL dadmaa( -1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1836 cc CALL hbook1(816,'WEIGHT DISTRIBUTION DEXAA $',100,-2.,2.)
1837 c
1838  ELSEIF(mode.EQ. 0) THEN
1839 * =======================
1840  300 CONTINUE
1841  IF(iwarm.EQ.0) goto 902
1842  CALL dadmaa( 0,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1843  wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1844 cc CALL hfill(816,wt)
1845  CALL ranmar(rn,1)
1846  IF(rn(1).GT.wt) goto 300
1847 c
1848  ELSEIF(mode.EQ. 1) THEN
1849 * =======================
1850  CALL dadmaa( 1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1851 cc CALL hprint(816)
1852  ENDIF
1853 c =====
1854  RETURN
1855  902 WRITE(iout, 9020)
1856  9020 FORMAT(' ----- DEXAA: LACK OF INITIALISATION')
1857  stop
1858  END
1859  SUBROUTINE dadmaa(MODE,ISGN,HHV,PNU,PAA,PIM1,PIM2,PIPL,JAA)
1860 c ----------------------------------------------------------------------
1861 * a1 decay unweighted events
1862 c ----------------------------------------------------------------------
1863  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1864  * ,ampiz,ampi,amro,gamro,ama1,gama1
1865  * ,amk,amkz,amkst,gamkst
1866 c
1867  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1868  * ,ampiz,ampi,amro,gamro,ama1,gama1
1869  * ,amk,amkz,amkst,gamkst
1870  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1871  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1872  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1873  REAL*4 gampmc ,gamper
1874  COMMON / inout / inut,iout
1875  REAL hhv(4)
1876  REAL hv(4),paa(4),pnu(4),pim1(4),pim2(4),pipl(4)
1877  REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4)
1878  REAL*4 rrr(3)
1879  REAL*8 swt, sswt
1880  DATA pi /3.141592653589793238462643/
1881  DATA iwarm/0/
1882 c
1883  IF(mode.EQ.-1) THEN
1884 c ===================
1885  iwarm=1
1886  nevraw=0
1887  nevacc=0
1888  nevovr=0
1889  swt=0
1890  sswt=0
1891  wtmax=1e-20
1892  DO 15 i=1,500
1893  CALL dphsaa(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5,jaa)
1894  IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1895 15 CONTINUE
1896 cc CALL hbook1(801,'WEIGHT DISTRIBUTION DADMAA $',100,0,2)
1897 c
1898  ELSEIF(mode.EQ. 0) THEN
1899 c =======================
1900 300 CONTINUE
1901  IF(iwarm.EQ.0) goto 902
1902  CALL dphsaa(wt,hv,pnu,paa,pim1,pim2,pipl,jaa)
1903 cc CALL hfill(801,wt/wtmax)
1904  nevraw=nevraw+1
1905  swt=swt+wt
1906 #if defined (ALEPH)
1907  sswt=sswt+wt**2
1908 #else
1909 ccm.s.>>>>>>
1910 cc sswt=sswt+wt**2
1911  sswt=sswt+dble(wt)**2
1912 ccm.s.<<<<<<
1913 #endif
1914  CALL ranmar(rrr,3)
1915  rn=rrr(1)
1916  IF(wt.GT.wtmax) nevovr=nevovr+1
1917  IF(rn*wtmax.GT.wt) goto 300
1918 c rotations to basic tau rest frame
1919  costhe=-1.+2.*rrr(2)
1920  thet=acos(costhe)
1921  phi =2*pi*rrr(3)
1922  CALL rotpol(thet,phi,pnu)
1923  CALL rotpol(thet,phi,paa)
1924  CALL rotpol(thet,phi,pim1)
1925  CALL rotpol(thet,phi,pim2)
1926  CALL rotpol(thet,phi,pipl)
1927  CALL rotpol(thet,phi,hv)
1928  DO 44 i=1,3
1929  44 hhv(i)=-isgn*hv(i)
1930  nevacc=nevacc+1
1931 c
1932  ELSEIF(mode.EQ. 1) THEN
1933 c =======================
1934  IF(nevraw.EQ.0) RETURN
1935  pargam=swt/float(nevraw+1)
1936  error=0
1937  IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1938  rat=pargam/gamel
1939  WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1940 cc CALL hprint(801)
1941  gampmc(5)=rat
1942  gamper(5)=error
1943 cam nevdec(5)=nevacc
1944  ENDIF
1945 c =====
1946  RETURN
1947  7003 FORMAT(///1x,15(5h*****)
1948  $ /,' *', 25x,'******** DADMAA INITIALISATION ********',9x,1h*
1949  $ /,' *',e20.5,5x,'WTMAX = MAXIMUM WEIGHT ',9x,1h*
1950  $ /,1x,15(5h*****)/)
1951  7010 FORMAT(///1x,15(5h*****)
1952  $ /,' *', 25x,'******** DADMAA FINAL REPORT ******** ',9x,1h*
1953  $ /,' *',i20 ,5x,'NEVRAW = NO. OF A1 DECAYS TOTAL ',9x,1h*
1954  $ /,' *',i20 ,5x,'NEVACC = NO. OF A1 DECS. ACCEPTED ',9x,1h*
1955  $ /,' *',i20 ,5x,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
1956  $ /,' *',e20.5,5x,'PARTIAL WTDTH (A1 DECAY) IN GEV UNITS ',9x,1h*
1957  $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1958  $ /,' *',f20.8,5x,'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
1959  $ /,1x,15(5h*****)/)
1960  902 WRITE(iout, 9020)
1961  9020 FORMAT(' ----- DADMAA: LACK OF INITIALISATION')
1962  stop
1963  END
1964  SUBROUTINE dphsaa(DGAMT,HV,PN,PAA,PIM1,PIM2,PIPL,JAA)
1965 c ----------------------------------------------------------------------
1966 * it simulates a1 decay in tau rest frame with
1967 * z-axis along a1 momentum
1968 c ----------------------------------------------------------------------
1969  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1970  * ,ampiz,ampi,amro,gamro,ama1,gama1
1971  * ,amk,amkz,amkst,gamkst
1972 c
1973  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1974  * ,ampiz,ampi,amro,gamro,ama1,gama1
1975  * ,amk,amkz,amkst,gamkst
1976  COMMON / taukle / bra1,brk0,brk0b,brks
1977  REAL*4 bra1,brk0,brk0b,brks
1978  REAL hv(4),pn(4),paa(4),pim1(4),pim2(4),pipl(4)
1979 
1980 
1981  REAL*4 rrr(1)
1982 c matrix element number:
1983  mnum=0
1984 c TYPE of the generation:
1985  keyt=1
1986  CALL ranmar(rrr,1)
1987  rmod=rrr(1)
1988  IF (rmod.LT.bra1) THEN
1989  jaa=1
1990  amp1=ampi
1991  amp2=ampi
1992  amp3=ampi
1993  ELSE
1994  jaa=2
1995  amp1=ampiz
1996  amp2=ampiz
1997  amp3=ampi
1998  ENDIF
1999  CALL
2000  $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
2001  END
2002  SUBROUTINE dexkk(MODE,ISGN,POL,PKK,PNU)
2003 c ----------------------------------------------------------------------
2004 c tau decay into kaon and tau-neutrino
2005 c in tau rest frame
2006 c output four momenta: pnu tauneutrino,
2007 c pkk kaon charged
2008 c ----------------------------------------------------------------------
2009  REAL pol(4),hv(4),pnu(4),pkk(4),rn(1)
2010 c
2011  IF(mode.EQ.-1) THEN
2012 c ===================
2013  CALL dadmkk(-1,isgn,hv,pkk,pnu)
2014 cc CALL hbook1(815,'WEIGHT DISTRIBUTION DEXPI $',100,0,2)
2015 c
2016  ELSEIF(mode.EQ. 0) THEN
2017 c =======================
2018 300 CONTINUE
2019  CALL dadmkk( 0,isgn,hv,pkk,pnu)
2020  wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
2021 cc CALL hfill(815,wt)
2022  CALL ranmar(rn,1)
2023  IF(rn(1).GT.wt) goto 300
2024 c
2025  ELSEIF(mode.EQ. 1) THEN
2026 c =======================
2027  CALL dadmkk( 1,isgn,hv,pkk,pnu)
2028 cc CALL hprint(815)
2029  ENDIF
2030 c =====
2031  RETURN
2032  END
2033  SUBROUTINE dadmkk(MODE,ISGN,HV,PKK,PNU)
2034 c ----------------------------------------------------------------------
2035 c fz
2036 #if defined (ALEPH)
2037 #else
2038  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2039  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2040 #endif
2041  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2042  * ,ampiz,ampi,amro,gamro,ama1,gama1
2043  * ,amk,amkz,amkst,gamkst
2044 c
2045  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2046  * ,ampiz,ampi,amro,gamro,ama1,gama1
2047  * ,amk,amkz,amkst,gamkst
2048 #if defined (ALEPH)
2049  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2050  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2051 #else
2052 #endif
2053  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
2054  REAL*4 gampmc ,gamper
2055  COMMON / inout / inut,iout
2056  REAL pkk(4),pnu(4),hv(4)
2057  DATA pi /3.141592653589793238462643/
2058 c
2059  IF(mode.EQ.-1) THEN
2060 c ===================
2061  nevtot=0
2062  ELSEIF(mode.EQ. 0) THEN
2063 c =======================
2064  nevtot=nevtot+1
2065  ekk= (amtau**2+amk**2-amnuta**2)/(2*amtau)
2066  enu= (amtau**2-amk**2+amnuta**2)/(2*amtau)
2067  xkk= sqrt(ekk**2-amk**2)
2068 c k momentum
2069  CALL sphera(xkk,pkk)
2070  pkk(4)=ekk
2071 c tau-neutrino momentum
2072  DO 30 i=1,3
2073 30 pnu(i)=-pkk(i)
2074  pnu(4)=enu
2075  pxq=amtau*ekk
2076  pxn=amtau*enu
2077  qxn=pkk(4)*pnu(4)-pkk(1)*pnu(1)-pkk(2)*pnu(2)-pkk(3)*pnu(3)
2078  brak=(gv**2+ga**2)*(2*pxq*qxn-amk**2*pxn)
2079  & +(gv**2-ga**2)*amtau*amnuta*amk**2
2080  DO 40 i=1,3
2081 40 hv(i)=-isgn*2*ga*gv*amtau*(2*pkk(i)*qxn-pnu(i)*amk**2)/brak
2082  hv(4)=1
2083 c
2084  ELSEIF(mode.EQ. 1) THEN
2085 c =======================
2086  IF(nevtot.EQ.0) RETURN
2087  fkk=0.0354
2088 cfz there was brak/amtau**4 before
2089 c gamm=(gfermi*fkk)**2/(16.*pi)*amtau**3*
2090 c * (brak/amtau**4)**2
2091 czw 7.02.93 here was an error affecting non standard model
2092 c configurations only
2093  gamm=(gfermi*fkk)**2/(16.*pi)*amtau**3*
2094  $ (brak/amtau**4)*
2095  $ sqrt((amtau**2-amk**2-amnuta**2)**2
2096  $ -4*amk**2*amnuta**2 )/amtau**2
2097  error=0
2098 
2099  error=0
2100  rat=gamm/gamel
2101  WRITE(iout, 7010) nevtot,gamm,rat,error
2102  gampmc(6)=rat
2103  gamper(6)=error
2104 cam nevdec(6)=nevtot
2105  ENDIF
2106 c =====
2107  RETURN
2108  7010 FORMAT(///1x,15(5h*****)
2109  $ /,' *', 25x,'******** DADMKK FINAL REPORT ********',9x,1h*
2110  $ /,' *',i20 ,5x,'NEVTOT = NO. OF K DECAYS TOTAL ',9x,1h*,
2111  $ /,' *',e20.5,5x,'PARTIAL WTDTH ( K DECAY) IN GEV UNITS ',9x,1h*,
2112  $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2113  $ /,' *',f20.8,5x,'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9x,1h*
2114  $ /,1x,15(5h*****)/)
2115  END
2116  SUBROUTINE dexks(MODE,ISGN,POL,PNU,PKS,PKK,PPI,JKST)
2117 c ----------------------------------------------------------------------
2118 c this simulates tau decay in tau rest frame
2119 c into nu k*, THEN k* decays into pi0,k+-(jkst=20)
2120 c or pi+-,k0(jkst=10).
2121 c output four momenta: pnu tauneutrino,
2122 c pks k* charged
2123 c pk0 k zero
2124 c pkc k charged
2125 c pic pion charged
2126 c piz pion zero
2127 c ----------------------------------------------------------------------
2128  COMMON / inout / inut,iout
2129  REAL pol(4),hv(4),pks(4),pnu(4),pkk(4),ppi(4),rn(1)
2130  DATA iwarm/0/
2131 c
2132  IF(mode.EQ.-1) THEN
2133 c ===================
2134  iwarm=1
2135 cfz initialisation done with the gharged pion neutral kaon mode(jkst=10
2136  CALL dadmks( -1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2137 cc CALL hbook1(816,'WEIGHT DISTRIBUTION DEXKS $',100,0,2)
2138 cc CALL hbook1(916,'ABS2 OF HV IN ROUTINE DEXKS $',100,0,2)
2139 c
2140  ELSEIF(mode.EQ. 0) THEN
2141 c =======================
2142 300 CONTINUE
2143  IF(iwarm.EQ.0) goto 902
2144  CALL dadmks( 0,isgn,hv,pnu,pks,pkk,ppi,jkst)
2145  wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
2146 cc CALL hfill(816,wt)
2147 cc xhelp=hv(1)**2+hv(2)**2+hv(3)**2
2148 cc CALL hfill(916,xhelp)
2149  CALL ranmar(rn,1)
2150  IF(rn(1).GT.wt) goto 300
2151 c
2152  ELSEIF(mode.EQ. 1) THEN
2153 c ======================================
2154  CALL dadmks( 1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2155 cc CALL hprint(816)
2156 cc CALL hprint(916)
2157  ENDIF
2158 c =====
2159  RETURN
2160  902 WRITE(iout, 9020)
2161  9020 FORMAT(' ----- DEXKS: LACK OF INITIALISATION')
2162  stop
2163  END
2164  SUBROUTINE dadmks(MODE,ISGN,HHV,PNU,PKS,PKK,PPI,JKST)
2165 c ----------------------------------------------------------------------
2166  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2167  * ,ampiz,ampi,amro,gamro,ama1,gama1
2168  * ,amk,amkz,amkst,gamkst
2169 c
2170  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2171  * ,ampiz,ampi,amro,gamro,ama1,gama1
2172  * ,amk,amkz,amkst,gamkst
2173  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2174  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2175  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
2176  REAL*4 gampmc ,gamper
2177  COMMON / taukle / bra1,brk0,brk0b,brks
2178  REAL*4 bra1,brk0,brk0b,brks
2179  COMMON / inout / inut,iout
2180  REAL hhv(4)
2181  REAL hv(4),pks(4),pnu(4),pkk(4),ppi(4)
2182  REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4)
2183  REAL*4 rrr(3),rmod(1)
2184  REAL*8 swt, sswt
2185  DATA pi /3.141592653589793238462643/
2186  DATA iwarm/0/
2187 c
2188  IF(mode.EQ.-1) THEN
2189 c ===================
2190  iwarm=1
2191  nevraw=0
2192  nevacc=0
2193  nevovr=0
2194  swt=0
2195  sswt=0
2196  wtmax=1e-20
2197  DO 15 i=1,500
2198 c the initialisation is done with the 66.7% MODE
2199  jkst=10
2200  CALL dphsks(wt,hv,pdum1,pdum2,pdum3,pdum4,jkst)
2201  IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
2202 15 CONTINUE
2203 cc CALL hbook1(801,'WEIGHT DISTRIBUTION DADMKS $',100,0,2)
2204 cc print 7003,wtmax
2205 cc CALL hbook1(112,'-------- K* MASS -------- $',100,0.,2.)
2206  ELSEIF(mode.EQ. 0) THEN
2207 c =====================================
2208  IF(iwarm.EQ.0) goto 902
2209 c here we choose randomly between k0 pi+_(66.7%)
2210 c and k+_ pi0(33.3%)
2211  dec1=brks
2212 400 CONTINUE
2213  CALL ranmar(rmod,1)
2214  IF(rmod(1).LT.dec1) THEN
2215  jkst=10
2216  ELSE
2217  jkst=20
2218  ENDIF
2219  CALL dphsks(wt,hv,pnu,pks,pkk,ppi,jkst)
2220  CALL ranmar(rrr,3)
2221  rn=rrr(1)
2222  IF(wt.GT.wtmax) nevovr=nevovr+1
2223  nevraw=nevraw+1
2224  swt=swt+wt
2225  sswt=sswt+wt**2
2226  IF(rn*wtmax.GT.wt) goto 400
2227 c rotations to basic tau rest frame
2228  costhe=-1.+2.*rrr(2)
2229  thet=acos(costhe)
2230  phi =2*pi*rrr(3)
2231  CALL rotor2(thet,pnu,pnu)
2232  CALL rotor3( phi,pnu,pnu)
2233  CALL rotor2(thet,pks,pks)
2234  CALL rotor3( phi,pks,pks)
2235  CALL rotor2(thet,pkk,pkk)
2236  CALL rotor3(phi,pkk,pkk)
2237  CALL rotor2(thet,ppi,ppi)
2238  CALL rotor3( phi,ppi,ppi)
2239  CALL rotor2(thet,hv,hv)
2240  CALL rotor3( phi,hv,hv)
2241  DO 44 i=1,3
2242  44 hhv(i)=-isgn*hv(i)
2243  nevacc=nevacc+1
2244 c
2245  ELSEIF(mode.EQ. 1) THEN
2246 c =======================
2247  IF(nevraw.EQ.0) RETURN
2248  pargam=swt/float(nevraw+1)
2249  error=0
2250  IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
2251  rat=pargam/gamel
2252  WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
2253 cc CALL hprint(801)
2254  gampmc(7)=rat
2255  gamper(7)=error
2256 cam nevdec(7)=nevacc
2257  ENDIF
2258 c =====
2259  RETURN
2260  7003 FORMAT(///1x,15(5h*****)
2261  $ /,' *', 25x,'******** DADMKS INITIALISATION ********',9x,1h*
2262  $ /,' *',e20.5,5x,'WTMAX = MAXIMUM WEIGHT ',9x,1h*
2263  $ /,1x,15(5h*****)/)
2264  7010 FORMAT(///1x,15(5h*****)
2265  $ /,' *', 25x,'******** DADMKS FINAL REPORT ********',9x,1h*
2266  $ /,' *',i20 ,5x,'NEVRAW = NO. OF K* DECAYS TOTAL ',9x,1h*,
2267  $ /,' *',i20 ,5x,'NEVACC = NO. OF K* DECS. ACCEPTED ',9x,1h*,
2268  $ /,' *',i20 ,5x,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
2269  $ /,' *',e20.5,5x,'PARTIAL WTDTH (K* DECAY) IN GEV UNITS ',9x,1h*,
2270  $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2271  $ /,' *',f20.8,5x,'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
2272  $ /,1x,15(5h*****)/)
2273  902 WRITE(iout, 9020)
2274  9020 FORMAT(' ----- DADMKS: LACK OF INITIALISATION')
2275  stop
2276  END
2277  SUBROUTINE dphsks(DGAMT,HV,PN,PKS,PKK,PPI,JKST)
2278 c ----------------------------------------------------------------------
2279 c it simulates kaon* decay in tau rest frame with
2280 c z-axis along kaon* momentum
2281 c jkst=10 for k* --->k0 + pi+-
2282 c jkst=20 for k* --->k+- + pi0
2283 c ----------------------------------------------------------------------
2284 #if defined (ALEPH)
2285 #else
2286  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2287  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2288 #endif
2289  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2290  * ,ampiz,ampi,amro,gamro,ama1,gama1
2291  * ,amk,amkz,amkst,gamkst
2292 c
2293  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2294  * ,ampiz,ampi,amro,gamro,ama1,gama1
2295  * ,amk,amkz,amkst,gamkst
2296 #if defined (ALEPH)
2297  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2298  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2299 #else
2300 #endif
2301  REAL hv(4),pt(4),pn(4),pks(4),pkk(4),ppi(4),qq(4),rr1(1)
2302 #if defined (ALEPH)
2303 cam COMPLEX bwigs
2304  COMPLEX bwigm
2305 #else
2306  COMPLEX bwigs
2307 #endif
2308  DATA pi /3.141592653589793238462643/
2309 c
2310  DATA icont /0/
2311 c three body phase space normalised as in bjorken-drell
2312  phspac=1./2**11/pi**5
2313 c tau momentum
2314  pt(1)=0.
2315  pt(2)=0.
2316  pt(3)=0.
2317  pt(4)=amtau
2318  CALL ranmar(rr1,1)
2319 c here begin the k0,pi+_ decay
2320  IF(jkst.EQ.10)THEN
2321 c ==================
2322 c mass of(real/virtual) k*
2323  ams1=(ampi+amkz)**2
2324  ams2=(amtau-amnuta)**2
2325 c flat phase space
2326 c amx2=ams1+ rr1(1)*(ams2-ams1)
2327 c amx=sqrt(amx2)
2328 c phspac=phspac*(ams2-ams1)
2329 c phase space with sampling for k* resonance
2330  alp1=atan((ams1-amkst**2)/amkst/gamkst)
2331  alp2=atan((ams2-amkst**2)/amkst/gamkst)
2332  alp=alp1+rr1(1)*(alp2-alp1)
2333  amx2=amkst**2+amkst*gamkst*tan(alp)
2334  amx=sqrt(amx2)
2335  phspac=phspac*((amx2-amkst**2)**2+(amkst*gamkst)**2)
2336  & /(amkst*gamkst)
2337  phspac=phspac*(alp2-alp1)
2338 c
2339 c tau-neutrino momentum
2340  pn(1)=0
2341  pn(2)=0
2342  pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2343  pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2344 c
2345 c k* momentum
2346  pks(1)=0
2347  pks(2)=0
2348  pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2349  pks(3)=-pn(3)
2350  phspac=phspac*(4*pi)*(2*pks(3)/amtau)
2351 c
2352 cam
2353  enpi=( amx**2+ampi**2-amkz**2 ) / ( 2*amx )
2354  pppi=sqrt((enpi-ampi)*(enpi+ampi))
2355  phspac=phspac*(4*pi)*(2*pppi/amx)
2356 c charged pi momentum in kaon* rest frame
2357  CALL sphera(pppi,ppi)
2358  ppi(4)=enpi
2359 c neutral kaon momentum in k* rest frame
2360  DO 20 i=1,3
2361 20 pkk(i)=-ppi(i)
2362  pkk(4)=( amx**2+amkz**2-ampi**2 ) / ( 2*amx )
2363  exe=(pks(4)+pks(3))/amx
2364 c pion and k boosted from k* rest frame to tau rest frame
2365  CALL bostr3(exe,ppi,ppi)
2366  CALL bostr3(exe,pkk,pkk)
2367  DO 30 i=1,4
2368 30 qq(i)=ppi(i)-pkk(i)
2369 c qq transverse to pks
2370  pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
2371  qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
2372  DO 31 i=1,4
2373 31 qq(i)=qq(i)-pks(i)*qqpks/pksd
2374 c amplitude
2375  prodpq=pt(4)*qq(4)
2376  prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
2377  prodpn=pt(4)*pn(4)
2378  qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
2379  brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
2380  & +(gv**2-ga**2)*amtau*amnuta*qq2
2381 c a simple breit-wigner is chosen for k* resonance
2382 #if defined (ALEPH)
2383 cam fks=cabs(bwigs(amx2,amkst,gamkst))**2
2384  fks=cabs(bwigm(amx2,amkst,gamkst,ampi,amkz))**2
2385 #else
2386  fks=cabs(bwigs(amx2,amkst,gamkst))**2
2387 #endif
2388  amplit=(gfermi*scabib)**2*brak*2*fks
2389  dgamt=1/(2.*amtau)*amplit*phspac
2390  DO 40 i=1,3
2391  40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
2392 c
2393 c here begin the k+-,pi0 decay
2394  ELSEIF(jkst.EQ.20)THEN
2395 c ======================
2396 c mass of(real/virtual) k*
2397  ams1=(ampiz+amk)**2
2398  ams2=(amtau-amnuta)**2
2399 c flat phase space
2400 #if defined (ALEPH)
2401 c amx2=ams1+ rr1(1)*(ams2-ams1)
2402 #else
2403 c amx2=ams1+ rr1*(ams2-ams1)
2404 #endif
2405 c amx=sqrt(amx2)
2406 c phspac=phspac*(ams2-ams1)
2407 c phase space with sampling for k* resonance
2408  alp1=atan((ams1-amkst**2)/amkst/gamkst)
2409  alp2=atan((ams2-amkst**2)/amkst/gamkst)
2410  alp=alp1+rr1(1)*(alp2-alp1)
2411  amx2=amkst**2+amkst*gamkst*tan(alp)
2412  amx=sqrt(amx2)
2413  phspac=phspac*((amx2-amkst**2)**2+(amkst*gamkst)**2)
2414  & /(amkst*gamkst)
2415  phspac=phspac*(alp2-alp1)
2416 c
2417 c tau-neutrino momentum
2418  pn(1)=0
2419  pn(2)=0
2420  pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2421  pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2422 c kaon* momentum
2423  pks(1)=0
2424  pks(2)=0
2425  pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2426  pks(3)=-pn(3)
2427  phspac=phspac*(4*pi)*(2*pks(3)/amtau)
2428 c
2429 cam
2430  enpi=( amx**2+ampiz**2-amk**2 ) / ( 2*amx )
2431  pppi=sqrt((enpi-ampiz)*(enpi+ampiz))
2432  phspac=phspac*(4*pi)*(2*pppi/amx)
2433 c neutral pi momentum in k* rest frame
2434  CALL sphera(pppi,ppi)
2435  ppi(4)=enpi
2436 c charged kaon momentum in k* rest frame
2437  DO 50 i=1,3
2438 50 pkk(i)=-ppi(i)
2439  pkk(4)=( amx**2+amk**2-ampiz**2 ) / ( 2*amx )
2440  exe=(pks(4)+pks(3))/amx
2441 c pion and k boosted from k* rest frame to tau rest frame
2442  CALL bostr3(exe,ppi,ppi)
2443  CALL bostr3(exe,pkk,pkk)
2444  DO 60 i=1,4
2445 60 qq(i)=pkk(i)-ppi(i)
2446 c qq transverse to pks
2447  pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
2448  qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
2449  DO 61 i=1,4
2450 61 qq(i)=qq(i)-pks(i)*qqpks/pksd
2451 c amplitude
2452  prodpq=pt(4)*qq(4)
2453  prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
2454  prodpn=pt(4)*pn(4)
2455  qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
2456  brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
2457  & +(gv**2-ga**2)*amtau*amnuta*qq2
2458 c a simple breit-wigner is chosen for the k* resonance
2459 #if defined (ALEPH)
2460 cam fks=cabs(bwigs(amx2,amkst,gamkst))**2
2461  fks=cabs(bwigm(amx2,amkst,gamkst,amk,ampiz))**2
2462 #else
2463  fks=cabs(bwigs(amx2,amkst,gamkst))**2
2464 #endif
2465  amplit=(gfermi*scabib)**2*brak*2*fks
2466  dgamt=1/(2.*amtau)*amplit*phspac
2467  DO 70 i=1,3
2468  70 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
2469  ENDIF
2470  RETURN
2471  END
2472 
2473 
2474 
2475 #if defined (ALEPH)
2476  SUBROUTINE dphnpi(DGAMT,HV,PN,PR,PPI,JNPI)
2477 #else
2478  SUBROUTINE dphnpi(DGAMT,HVX,PNX,PRX,PPIX,JNPI)
2479 #endif
2480 c ----------------------------------------------------------------------
2481 c it simulates multipi decay in tau rest frame with
2482 c z-axis opposite to neutrino momentum
2483 c ----------------------------------------------------------------------
2484  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2485  * ,ampiz,ampi,amro,gamro,ama1,gama1
2486  * ,amk,amkz,amkst,gamkst
2487 c
2488  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2489  * ,ampiz,ampi,amro,gamro,ama1,gama1
2490  * ,amk,amkz,amkst,gamkst
2491  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2492  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2493  parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
2494 #if defined (ALEPH)
2495  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2496  & ,names
2497  CHARACTER names(nmode)*31
2498 c
2499 #else
2500  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2501  & ,names
2502  CHARACTER names(nmode)*31
2503  REAL*8 wetmax(20)
2504 c
2505 #endif
2506 #if defined (ALEPH)
2507  REAL pn(4),pr(4),ppi(4,9),hv(4)
2508  REAL pv(5,9),pt(4),ue(3),be(3)
2509  REAL*4 rrr(9),rord(9),rr1(1)
2510  real dpar(8)
2511 c
2512  DATA pi /3.141592653589793238462643/
2513  DATA dpar/2.,5.,15.,60.,250.,1500.,1.2e4,1.2e5/
2514 c
2515 c pawt(a,b,c)=sqrt((a**2-(b+c)**2)*(a**2-(b-c)**2))/(2.*a)
2516  pawt(a,b,c)=sqrt(max(0.,(a**2-(b+c)**2)*(a**2-(b-c)**2)))/(2.*a)
2517 #else
2518  REAL*8 pn(4),pr(4),ppi(4,9),hv(4)
2519  REAL*4 pnx(4),prx(4),ppix(4,9),hvx(4)
2520  REAL*8 pv(5,9),pt(4),ue(3),be(3)
2521  REAL*8 pawt,amx,ams1,ams2,pa,phs,phsmax,pmin,pmax
2522 !!! M.S. to fix underflow >>>
2523  REAL*8 phspac
2524 !!! M.S. to fix underflow <<<
2525  REAL*8 gam,bep,phi,a,b,c
2526  REAL*8 ampik
2527  REAL*4 rrr(9),rrx(2),rn(1),rr2(1)
2528 c
2529  DATA pi /3.141592653589793238462643/
2530  DATA wetmax /20*1d-15/
2531 c
2532 cc-- pawt(a,b,c)=sqrt((a**2-(b+c)**2)*(a**2-(b-c)**2))/(2.*a)
2533 c
2534  pawt(a,b,c)=
2535  $ sqrt(max(0.d0,(a**2-(b+c)**2)*(a**2-(b-c)**2)))/(2.d0*a)
2536 #endif
2537 c
2538  ampik(i,j)=dcdmas(idffin(i,j))
2539 c
2540 c
2541 #if defined (ALEPH)
2542 #else
2543  IF ((jnpi.LE.0).OR.jnpi.GT.20) THEN
2544  WRITE(6,*) 'JNPI OUTSIDE RANGE DEFINED BY WETMAX; JNPI=',jnpi
2545  stop
2546  ENDIF
2547 
2548 #endif
2549 c tau momentum
2550  pt(1)=0.
2551  pt(2)=0.
2552  pt(3)=0.
2553  pt(4)=amtau
2554 c
2555 #if defined (ALEPH)
2556 #else
2557  500 CONTINUE
2558 #endif
2559 c mass of virtual w
2560  nd=mulpik(jnpi)
2561  ps=0.
2562  phspac = 1./2.**5 /pi**2
2563  DO 4 i=1,nd
2564 4 ps =ps+ampik(i,jnpi)
2565 #if defined (ALEPH)
2566  CALL ranmar(rr1,1)
2567 #else
2568  CALL ranmar(rr2,1)
2569 #endif
2570  ams1=ps**2
2571  ams2=(amtau-amnuta)**2
2572 c
2573 c
2574 #if defined (ALEPH)
2575  amx2=ams1+ rr1(1)*(ams2-ams1)
2576 #else
2577  amx2=ams1+ rr2(1)*(ams2-ams1)
2578 #endif
2579  amx =sqrt(amx2)
2580  amw =amx
2581  phspac=phspac * (ams2-ams1)
2582 c
2583 c tau-neutrino momentum
2584  pn(1)=0
2585  pn(2)=0
2586  pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx2)
2587  pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2588 c w momentum
2589  pr(1)=0
2590  pr(2)=0
2591  pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx2)
2592  pr(3)=-pn(3)
2593  phspac=phspac * (4.*pi) * (2.*pr(3)/amtau)
2594 c
2595 c amplitude(cf ys.tsai phys.rev.d4,2821(1971)
2596 c or f.gilman sh.rhie phys.rev.d31,1066(1985)
2597 c
2598  pxq=amtau*pr(4)
2599  pxn=amtau*pn(4)
2600  qxn=pr(4)*pn(4)-pr(1)*pn(1)-pr(2)*pn(2)-pr(3)*pn(3)
2601 #if defined (ALEPH)
2602 #else
2603 c here was an error. 20.10.91 (zw)
2604 c brak=2*(gv**2+ga**2)*(2*pxq*pxn+amx2*qxn)
2605 #endif
2606  brak=2*(gv**2+ga**2)*(2*pxq*qxn+amx2*pxn)
2607  & -6*(gv**2-ga**2)*amtau*amnuta*amx2
2608 cam assume neutrino mass=0. and sum over final polarisation
2609 c brak= 2*(amtau**2-amx2) * (amtau**2+2.*amx2)
2610  amplit=ccabib**2*gfermi**2/2. * brak * amx2*sigee(amx2,jnpi)
2611  dgamt=1./(2.*amtau)*amplit*phspac
2612 c
2613 c isotropic w decay in w rest frame
2614 #if defined (ALEPH)
2615  phspac = 1./2.**(6*nd-7) /pi**(3*nd-4)
2616  phsmax = 1./dpar(nd-2)
2617 #else
2618  phsmax = 1.
2619 #endif
2620  DO 200 i=1,4
2621  200 pv(i,1)=pr(i)
2622  pv(5,1)=amw
2623  pv(5,nd)=ampik(nd,jnpi)
2624 c compute max. phase space factor
2625  pmax=amw-ps+ampik(nd,jnpi)
2626  pmin=.0
2627  DO 220 il=nd-1,1,-1
2628  pmax=pmax+ampik(il,jnpi)
2629  pmin=pmin+ampik(il+1,jnpi)
2630 #if defined (ALEPH)
2631  220 phsmax=phsmax*pawt(pmax,pmin,ampik(il,jnpi))
2632 cam generate nd-2 effective masses(cf ludecy)
2633  phspac = 1./2.**(6*nd-7) /pi**(3*nd-4)
2634  240 rord(1)=1.
2635  CALL ranmar(rrr,nd-1)
2636  DO 260 il=2,nd-1
2637  rsav=rrr(il)
2638  DO 250 jl=il-1,1,-1
2639  IF(rsav.LE.rord(jl)) goto 260
2640  250 rord(jl+1)=rord(jl)
2641  260 rord(jl+1)=rsav
2642  rord(nd)=0.
2643  phs=1.
2644  DO 270 il=nd-1,1,-1
2645  pv(5,il)=pv(5,il+1)+ampik(il,jnpi)
2646  & +(rord(il)-rord(il+1))*(pv(5,1)-ps)
2647  270 phs=phs*pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2648  rn = rrr(1)
2649  IF(phs.LT.rn*phsmax) goto 240
2650 #else
2651  220 phsmax=phsmax*pawt(pmax,pmin,ampik(il,jnpi))/pmax
2652 
2653 c --- 2.02.94 zw 9 lines
2654  amx=amw
2655  DO 222 il=1,nd-2
2656  ams1=.0
2657  DO 223 jl=il+1,nd
2658  223 ams1=ams1+ampik(jl,jnpi)
2659  ams1=ams1**2
2660  amx =(amx-ampik(il,jnpi))
2661  ams2=(amx)**2
2662  phsmax=phsmax * (ams2-ams1)
2663  222 CONTINUE
2664  ncont=0
2665  100 CONTINUE
2666  ncont=ncont+1
2667 cam generate nd-2 effective masses
2668  phs=1.d0
2669  phspac = 1./2.**(6*nd-7) /pi**(3*nd-4)
2670  amx=amw
2671  CALL ranmar(rrr,nd-2)
2672  DO 230 il=1,nd-2
2673  ams1=.0d0
2674  DO 231 jl=il+1,nd
2675  231 ams1=ams1+ampik(jl,jnpi)
2676  ams1=ams1**2
2677  ams2=(amx-ampik(il,jnpi))**2
2678  rr1=rrr(il)
2679  amx2=ams1+ rr1*(ams2-ams1)
2680  amx=sqrt(amx2)
2681  pv(5,il+1)=amx
2682  phspac=phspac * (ams2-ams1)
2683 c --- 2.02.94 zw 1 line
2684  phs=phs* (ams2-ams1)
2685  pa=pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2686  phs =phs *pa/pv(5,il)
2687  230 CONTINUE
2688  pa=pawt(pv(5,nd-1),ampik(nd-1,jnpi),ampik(nd,jnpi))
2689  phs =phs *pa/pv(5,nd-1)
2690  CALL ranmar(rn,1)
2691  wetmax(jnpi)=1.2d0*max(wetmax(jnpi)/1.2d0,phs/phsmax)
2692  IF (ncont.EQ.500 000) THEN
2693  xnpi=0.0
2694  DO kk=1,nd
2695  xnpi=xnpi+ampik(kk,jnpi)
2696  ENDDO
2697  WRITE(6,*) 'ROUNDING INSTABILITY IN DPHNPI ?'
2698  WRITE(6,*) 'AMW=',amw,'XNPI=',xnpi
2699  WRITE(6,*) 'IF =AMW= IS NEARLY EQUAL =XNPI= THAT IS IT'
2700  WRITE(6,*) 'PHS=',phs,'PHSMAX=',phsmax
2701  goto 500
2702  ENDIF
2703  IF(rn(1)*phsmax*wetmax(jnpi).GT.phs) go to 100
2704 #endif
2705 c...perform successive two-particle decays in respective cm frame
2706  280 DO 300 il=1,nd-1
2707  pa=pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2708 #if defined (ALEPH)
2709  CALL ranmar(rrr,2)
2710  ue(3)=2.*rrr(1)-1.
2711  phi=2.*pi*rrr(2)
2712  ue(1)=sqrt(1.-ue(3)**2)*cos(phi)
2713  ue(2)=sqrt(1.-ue(3)**2)*sin(phi)
2714 #else
2715  CALL ranmar(rrx,2)
2716  ue(3)=2.*rrx(1)-1.
2717  phi=2.*pi*rrx(2)
2718  ue(1)=sqrt(1.d0-ue(3)**2)*cos(phi)
2719  ue(2)=sqrt(1.d0-ue(3)**2)*sin(phi)
2720 #endif
2721  DO 290 j=1,3
2722  ppi(j,il)=pa*ue(j)
2723  290 pv(j,il+1)=-pa*ue(j)
2724  ppi(4,il)=sqrt(pa**2+ampik(il,jnpi)**2)
2725  pv(4,il+1)=sqrt(pa**2+pv(5,il+1)**2)
2726  phspac=phspac *(4.*pi)*(2.*pa/pv(5,il))
2727  300 CONTINUE
2728 c...lorentz transform decay products to tau frame
2729  DO 310 j=1,4
2730  310 ppi(j,nd)=pv(j,nd)
2731  DO 340 il=nd-1,1,-1
2732  DO 320 j=1,3
2733  320 be(j)=pv(j,il)/pv(4,il)
2734  gam=pv(4,il)/pv(5,il)
2735  DO 340 i=il,nd
2736  bep=be(1)*ppi(1,i)+be(2)*ppi(2,i)+be(3)*ppi(3,i)
2737  DO 330 j=1,3
2738 #if defined (ALEPH)
2739  330 ppi(j,i)=ppi(j,i)+gam*(gam*bep/(1.+gam)+ppi(4,i))*be(j)
2740 #else
2741  330 ppi(j,i)=ppi(j,i)+gam*(gam*bep/(1.d0+gam)+ppi(4,i))*be(j)
2742 #endif
2743  ppi(4,i)=gam*(ppi(4,i)+bep)
2744  340 CONTINUE
2745 c
2746  hv(4)=1.
2747  hv(3)=0.
2748  hv(2)=0.
2749  hv(1)=0.
2750 #if defined (ALEPH)
2751 #else
2752  DO k=1,4
2753  pnx(k)=pn(k)
2754  prx(k)=pr(k)
2755  hvx(k)=hv(k)
2756  DO l=1,nd
2757  ppix(k,l)=ppi(k,l)
2758  ENDDO
2759  ENDDO
2760 #endif
2761  RETURN
2762  END
2763  FUNCTION sigee(Q2,JNP)
2764 c ----------------------------------------------------------------------
2765 c e+e- cross section in the(1.gev2,amtau**2) region
2766 c normalised to sig0 = 4/3 pi alfa2
2767 c used in matrix element for multipion tau decays
2768 c cf ys.tsai phys.rev d4 ,2821(1971)
2769 c f.gilman et al phys.rev d17,1846(1978)
2770 c c.kiesling, to be pub. in high energy e+e- physics(1988)
2771 c datsig(*,1) = e+e- -> pi+pi-2pi0
2772 c datsig(*,2) = e+e- -> 2pi+2pi-
2773 c datsig(*,3) = 5-pion contribution(a la tn.pham et al)
2774 c(phys lett 78b,623(1978)
2775 c datsig(*,5) = e+e- -> 6pi
2776 c
2777 c 4- and 6-pion cross sections from data
2778 c 5-pion contribution related to 4-pion cross section
2779 c
2780 c called by dphnpi
2781 c ----------------------------------------------------------------------
2782  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2783  * ,ampiz,ampi,amro,gamro,ama1,gama1
2784  * ,amk,amkz,amkst,gamkst
2785 c
2786  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2787  * ,ampiz,ampi,amro,gamro,ama1,gama1
2788  * ,amk,amkz,amkst,gamkst
2789  REAL*4 datsig(17,6)
2790 c
2791  DATA datsig/
2792  1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
2793  2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
2794  3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
2795  4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
2796  5 17*.0,
2797  6 17*.0,
2798  7 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25,
2799  8 17*.0/
2800  DATA sig0 / 86.8 /
2801  DATA pi /3.141592653589793238462643/
2802  DATA init / 0 /
2803 c
2804  jnpi=jnp
2805  IF(jnp.EQ.4) jnpi=3
2806  IF(jnp.EQ.3) jnpi=4
2807  IF(init.EQ.0) THEN
2808  init=1
2809 #if defined (CLEO)
2810 c ajwmod: initialize if called from outside qq:
2811  IF (ampi.LT.0.139) ampi = 0.1395675
2812 #else
2813 #endif
2814  ampi2=ampi**2
2815  fpi = .943*ampi
2816  DO 100 i=1,17
2817  datsig(i,2) = datsig(i,2)/2.
2818  datsig(i,1) = datsig(i,1) + datsig(i,2)
2819  s = 1.025+(i-1)*.05
2820  fact=0.
2821  s2=s**2
2822  DO 200 j=1,17
2823  t= 1.025+(j-1)*.05
2824  IF(t . gt. s-ampi ) go to 201
2825  t2=t**2
2826  fact=(t2/s2)**2*sqrt((s2-t2-ampi2)**2-4.*t2*ampi2)/s2 *2.*t*.05
2827  fact = fact * (datsig(j,1)+datsig(j+1,1))
2828  200 datsig(i,3) = datsig(i,3) + fact
2829  201 datsig(i,3) = datsig(i,3) /(2*pi*fpi)**2
2830  datsig(i,4) = datsig(i,3)
2831  datsig(i,6) = datsig(i,5)
2832  100 CONTINUE
2833 c WRITE(6,1000) datsig
2834  1000 FORMAT(///1x,' EE SIGMA USED IN MULTIPI DECAYS'/
2835  % (17f7.2/))
2836  ENDIF
2837  q=sqrt(q2)
2838  qmin=1.
2839  IF(q.LT.qmin) THEN
2840  sigee=datsig(1,jnpi)+
2841  & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
2842  ELSEIF(q.LT.1.8) THEN
2843  DO 1 i=1,16
2844  qmax = qmin + .05
2845  IF(q.LT.qmax) go to 2
2846  qmin = qmin + .05
2847  1 CONTINUE
2848  2 sigee=datsig(i,jnpi)+
2849  & (datsig(i+1,jnpi)-datsig(i,jnpi)) * (q-qmin)/.05
2850  ELSEIF(q.GT.1.8) THEN
2851  sigee=datsig(17,jnpi)+
2852  & (datsig(17,jnpi)-datsig(16,jnpi)) * (q-1.8)/.05
2853  ENDIF
2854  IF(sigee.LT..0) sigee=0.
2855 c
2856  sigee = sigee/(6.*pi**2*sig0)
2857 c
2858  RETURN
2859  END
2860 
2861  FUNCTION sigold(Q2,JNPI)
2862 c ----------------------------------------------------------------------
2863 c e+e- cross section in the(1.gev2,amtau**2) region
2864 c normalised to sig0 = 4/3 pi alfa2
2865 c used in matrix element for multipion tau decays
2866 c cf ys.tsai phys.rev d4 ,2821(1971)
2867 c f.gilman et al phys.rev d17,1846(1978)
2868 c c.kiesling, to be pub. in high energy e+e- physics(1988)
2869 c datsig(*,1) = e+e- -> pi+pi-2pi0
2870 c datsig(*,2) = e+e- -> 2pi+2pi-
2871 c datsig(*,3) = 5-pion contribution(a la tn.pham et al)
2872 c(phys lett 78b,623(1978)
2873 c datsig(*,4) = e+e- -> 6pi
2874 c
2875 c 4- and 6-pion cross sections from data
2876 c 5-pion contribution related to 4-pion cross section
2877 c
2878 c called by dphnpi
2879 c ----------------------------------------------------------------------
2880  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2881  * ,ampiz,ampi,amro,gamro,ama1,gama1
2882  * ,amk,amkz,amkst,gamkst
2883 c
2884  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2885  * ,ampiz,ampi,amro,gamro,ama1,gama1
2886  * ,amk,amkz,amkst,gamkst
2887  REAL*4 datsig(17,4)
2888 c
2889  DATA datsig/
2890  1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
2891  2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
2892  3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
2893  4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
2894  5 17*.0,
2895  6 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25/
2896  DATA sig0 / 86.8 /
2897  DATA pi /3.141592653589793238462643/
2898  DATA init / 0 /
2899 c
2900  IF(init.EQ.0) THEN
2901  init=1
2902  ampi2=ampi**2
2903  fpi = .943*ampi
2904  DO 100 i=1,17
2905  datsig(i,2) = datsig(i,2)/2.
2906  datsig(i,1) = datsig(i,1) + datsig(i,2)
2907  s = 1.025+(i-1)*.05
2908  fact=0.
2909  s2=s**2
2910  DO 200 j=1,17
2911  t= 1.025+(j-1)*.05
2912  IF(t . gt. s-ampi ) go to 201
2913  t2=t**2
2914  fact=(t2/s2)**2*sqrt((s2-t2-ampi2)**2-4.*t2*ampi2)/s2 *2.*t*.05
2915  fact = fact * (datsig(j,1)+datsig(j+1,1))
2916  200 datsig(i,3) = datsig(i,3) + fact
2917  201 datsig(i,3) = datsig(i,3) /(2*pi*fpi)**2
2918  100 CONTINUE
2919 c WRITE(6,1000) datsig
2920  1000 FORMAT(///1x,' EE SIGMA USED IN MULTIPI DECAYS'/
2921  % (17f7.2/))
2922  ENDIF
2923  q=sqrt(q2)
2924  qmin=1.
2925  IF(q.LT.qmin) THEN
2926  sigee=datsig(1,jnpi)+
2927  & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
2928  ELSEIF(q.LT.1.8) THEN
2929  DO 1 i=1,16
2930  qmax = qmin + .05
2931  IF(q.LT.qmax) go to 2
2932  qmin = qmin + .05
2933  1 CONTINUE
2934  2 sigee=datsig(i,jnpi)+
2935  & (datsig(i+1,jnpi)-datsig(i,jnpi)) * (q-qmin)/.05
2936  ELSEIF(q.GT.1.8) THEN
2937  sigee=datsig(17,jnpi)+
2938  & (datsig(17,jnpi)-datsig(16,jnpi)) * (q-1.8)/.05
2939  ENDIF
2940  IF(sigee.LT..0) sigee=0.
2941 c
2942  sigee = sigee/(6.*pi**2*sig0)
2943  sigold=sigee
2944 c
2945  RETURN
2946  END
2947  SUBROUTINE dphspk(DGAMT,HV,PN,PAA,PNPI,JAA)
2948 c ----------------------------------------------------------------------
2949 * it simulates three pi(k) decay in the tau rest frame
2950 * z-axis along hadronic system
2951 c ----------------------------------------------------------------------
2952  parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
2953 #if defined (ALEPH)
2954  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2955 #else
2956  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2957 #endif
2958  & ,names
2959  CHARACTER names(nmode)*31
2960 
2961  REAL hv(4),pn(4),paa(4),pim1(4),pim2(4),pipl(4),pnpi(4,9)
2962 c matrix element number:
2963  mnum=jaa
2964 c TYPE of the generation:
2965  keyt=4
2966  IF(jaa.EQ.7) keyt=3
2967 c --- masses of the decay products
2968  amp1=dcdmas(idffin(1,jaa+nm4+nm5+nm6))
2969  amp2=dcdmas(idffin(2,jaa+nm4+nm5+nm6))
2970  amp3=dcdmas(idffin(3,jaa+nm4+nm5+nm6))
2971  CALL
2972  $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
2973  DO i=1,4
2974  pnpi(i,1)=pim1(i)
2975  pnpi(i,2)=pim2(i)
2976  pnpi(i,3)=pipl(i)
2977  ENDDO
2978  END
2979 
2980 
2981 
2982 
2983  subroutine
2984  $ dphtre(dgamt,hv,pn,paa,pim1,ampa,pim2,ampb,pipl,amp3,keyt,mnum)
2985 c ----------------------------------------------------------------------
2986 * it simulates a1 decay in tau rest frame with
2987 * z-axis along a1 momentum
2988 * it can be also used to generate k k pi and k pi pi tau decays.
2989 * input parameters
2990 * keyt - algorithm controlling switch
2991 * 2 - flat phase space pim1 pim2 symmetrized statistical factor 1/2
2992 * 1 - like 1 but peaked around a1 and rho(two channels) masses.
2993 * 3 - peaked around omega, all particles different
2994 * other- flat phase space, all particles different
2995 * amp1 - mass of first pi, etc. (1-3)
2996 * mnum - matrix element type
2997 * 0 - a1 matrix element
2998 * 1-6 - matrix element for k pi pi, k k pi decay modes
2999 * 7 - pi- pi0 gamma matrix element
3000 c ----------------------------------------------------------------------
3001  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3002  * ,ampiz,ampi,amro,gamro,ama1,gama1
3003  * ,amk,amkz,amkst,gamkst
3004 c
3005  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3006  * ,ampiz,ampi,amro,gamro,ama1,gama1
3007  * ,amk,amkz,amkst,gamkst
3008  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3009  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3010  REAL hv(4),pt(4),pn(4),paa(4),pim1(4),pim2(4),pipl(4)
3011  REAL pr(4)
3012  REAL*4 rrr(5)
3013  DATA pi /3.141592653589793238462643/
3014  DATA icont /0/
3015  xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
3016 c amro, gamro is only a PARAMETER for geting hight efficiency
3017 c
3018 c three body phase space normalised as in bjorken-drell
3019 c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
3020  phspac=1./2**17/pi**8
3021 c tau momentum
3022  pt(1)=0.
3023  pt(2)=0.
3024  pt(3)=0.
3025  pt(4)=amtau
3026 c
3027  CALL ranmar(rrr,5)
3028  rr=rrr(5)
3029 c
3030  CALL choice(mnum,rr,ichan,prob1,prob2,prob3,
3031  $ amrx,gamrx,amra,gamra,amrb,gamrb)
3032  IF (ichan.EQ.1) THEN
3033  amp1=ampb
3034  amp2=ampa
3035  ELSEIF (ichan.EQ.2) THEN
3036  amp1=ampa
3037  amp2=ampb
3038  ELSE
3039  amp1=ampb
3040  amp2=ampa
3041  ENDIF
3042 cam
3043  rr1=rrr(1)
3044  ams1=(amp1+amp2+amp3)**2
3045  ams2=(amtau-amnuta)**2
3046 #if defined (ALEPH)
3047 c phase space with sampling for a1 resonance
3048 #else
3049 * phase space with sampling for a1 resonance
3050 #endif
3051  alp1=atan((ams1-amrx**2)/amrx/gamrx)
3052  alp2=atan((ams2-amrx**2)/amrx/gamrx)
3053  alp=alp1+rr1*(alp2-alp1)
3054  am3sq =amrx**2+amrx*gamrx*tan(alp)
3055  am3 =sqrt(am3sq)
3056  phspac=phspac*((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
3057  phspac=phspac*(alp2-alp1)
3058 c mass of(real/virtual) rho -
3059  rr2=rrr(2)
3060  ams1=(amp2+amp3)**2
3061  ams2=(am3-amp1)**2
3062  IF (ichan.LE.2) THEN
3063 #if defined (ALEPH)
3064 c phase space with sampling for rho resonance,
3065 #else
3066 * phase space with sampling for rho resonance,
3067 #endif
3068  alp1=atan((ams1-amra**2)/amra/gamra)
3069  alp2=atan((ams2-amra**2)/amra/gamra)
3070  alp=alp1+rr2*(alp2-alp1)
3071  am2sq =amra**2+amra*gamra*tan(alp)
3072  am2 =sqrt(am2sq)
3073 c --- this part of the jacobian will be recovered later ---------------
3074 c phspac=phspac*(alp2-alp1)
3075 c phspac=phspac*((am2sq-amra**2)**2+(amra*gamra)**2)/(amra*gamra)
3076 c----------------------------------------------------------------------
3077  ELSE
3078 #if defined (ALEPH)
3079 c flat phase space;
3080 #else
3081 * flat phase space;
3082 #endif
3083  am2sq=ams1+ rr2*(ams2-ams1)
3084  am2 =sqrt(am2sq)
3085  phf0=(ams2-ams1)
3086  ENDIF
3087 #if defined (ALEPH)
3088 c rho restframe, define pipl and pim1
3089 #else
3090 * rho restframe, define pipl and pim1
3091 #endif
3092  enq1=(am2sq-amp2**2+amp3**2)/(2*am2)
3093  enq2=(am2sq+amp2**2-amp3**2)/(2*am2)
3094  ppi= enq1**2-amp3**2
3095  pppi=sqrt(abs(enq1**2-amp3**2))
3096 c --- this part of jacobian will be recovered later
3097  phf1=(4*pi)*(2*pppi/am2)
3098 #if defined (ALEPH)
3099 c pi minus momentum in rho rest frame
3100 #else
3101 * pi minus momentum in rho rest frame
3102 #endif
3103  CALL sphera(pppi,pipl)
3104  pipl(4)=enq1
3105 #if defined (ALEPH)
3106 c pi0 1 momentum in rho rest frame
3107 #else
3108 * pi0 1 momentum in rho rest frame
3109 #endif
3110  DO 30 i=1,3
3111  30 pim1(i)=-pipl(i)
3112  pim1(4)=enq2
3113 #if defined (ALEPH)
3114 c a1 rest frame, define pim2
3115 #else
3116 * a1 rest frame, define pim2
3117 #endif
3118 * rho momentum
3119  pr(1)=0
3120  pr(2)=0
3121  pr(4)=1./(2*am3)*(am3**2+am2**2-amp1**2)
3122  pr(3)= sqrt(abs(pr(4)**2-am2**2))
3123  ppi = pr(4)**2-am2**2
3124 * pi0 2 momentum
3125  pim2(1)=0
3126  pim2(2)=0
3127  pim2(4)=1./(2*am3)*(am3**2-am2**2+amp1**2)
3128  pim2(3)=-pr(3)
3129  phf2=(4*pi)*(2*pr(3)/am3)
3130 #if defined (ALEPH)
3131 c old pions boosted from rho rest frame to a1 rest frame
3132 #else
3133 * old pions boosted from rho rest frame to a1 rest frame
3134 #endif
3135  exe=(pr(4)+pr(3))/am2
3136  CALL bostr3(exe,pipl,pipl)
3137  CALL bostr3(exe,pim1,pim1)
3138  rr3=rrr(3)
3139  rr4=rrr(4)
3140 #if defined (ALEPH)
3141 #else
3142 cam thet =pi*rr3
3143 #endif
3144  thet =acos(-1.+2*rr3)
3145  phi = 2*pi*rr4
3146  CALL rotpol(thet,phi,pipl)
3147  CALL rotpol(thet,phi,pim1)
3148  CALL rotpol(thet,phi,pim2)
3149  CALL rotpol(thet,phi,pr)
3150 c
3151 * now to the tau rest frame, define a1 and neutrino momenta
3152 * a1 momentum
3153  paa(1)=0
3154  paa(2)=0
3155  paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am3**2)
3156  paa(3)= sqrt(abs(paa(4)**2-am3**2))
3157  ppi = paa(4)**2-am3**2
3158  phspac=phspac*(4*pi)*(2*paa(3)/amtau)
3159 * tau-neutrino momentum
3160  pn(1)=0
3161  pn(2)=0
3162  pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am3**2)
3163  pn(3)=-paa(3)
3164 c here we correct for the jacobians of the two chains
3165 c ---first channel ------- pim1+pipl
3166  ams1=(amp2+amp3)**2
3167  ams2=(am3-amp1)**2
3168  alp1=atan((ams1-amra**2)/amra/gamra)
3169  alp2=atan((ams2-amra**2)/amra/gamra)
3170  xpro = (pim1(3)+pipl(3))**2
3171  $ +(pim1(2)+pipl(2))**2+(pim1(1)+pipl(1))**2
3172  am2sq=-xpro+(pim1(4)+pipl(4))**2
3173 c jacobian of speeding
3174  ff1 = ((am2sq-amra**2)**2+(amra*gamra)**2)/(amra*gamra)
3175  ff1 =ff1 *(alp2-alp1)
3176 c lambda of rho decay
3177  gg1 = (4*pi)*(xlam(am2sq,amp2**2,amp3**2)/am2sq)
3178 c lambda of a1 decay
3179  gg1 =gg1 *(4*pi)*sqrt(4*xpro/am3sq)
3180  xjaje=gg1*(ams2-ams1)
3181 c ---second channel ------ pim2+pipl
3182  ams1=(amp1+amp3)**2
3183  ams2=(am3-amp2)**2
3184  alp1=atan((ams1-amrb**2)/amrb/gamrb)
3185  alp2=atan((ams2-amrb**2)/amrb/gamrb)
3186  xpro = (pim2(3)+pipl(3))**2
3187  $ +(pim2(2)+pipl(2))**2+(pim2(1)+pipl(1))**2
3188  am2sq=-xpro+(pim2(4)+pipl(4))**2
3189  ff2 = ((am2sq-amrb**2)**2+(amrb*gamrb)**2)/(amrb*gamrb)
3190  ff2 =ff2 *(alp2-alp1)
3191  gg2 = (4*pi)*(xlam(am2sq,amp1**2,amp3**2)/am2sq)
3192  gg2 =gg2 *(4*pi)*sqrt(4*xpro/am3sq)
3193  xjadw=gg2*(ams2-ams1)
3194 c
3195  a1=0.0
3196  a2=0.0
3197  a3=0.0
3198  xjac1=ff1*gg1
3199  xjac2=ff2*gg2
3200  IF (ichan.EQ.2) THEN
3201  xjac3=xjadw
3202  ELSE
3203  xjac3=xjaje
3204  ENDIF
3205  IF (xjac1.NE.0.0) a1=prob1/xjac1
3206  IF (xjac2.NE.0.0) a2=prob2/xjac2
3207  IF (xjac3.NE.0.0) a3=prob3/xjac3
3208 c
3209  IF (a1+a2+a3.NE.0.0) THEN
3210  phspac=phspac/(a1+a2+a3)
3211  ELSE
3212  phspac=0.0
3213  ENDIF
3214  IF(ichan.EQ.2) THEN
3215  DO 70 i=1,4
3216  x=pim1(i)
3217  pim1(i)=pim2(i)
3218  70 pim2(i)=x
3219  ENDIF
3220 * all pions boosted from a1 rest frame to tau rest frame
3221 * z-axis antiparallel to neutrino momentum
3222  exe=(paa(4)+paa(3))/am3
3223  CALL bostr3(exe,pipl,pipl)
3224  CALL bostr3(exe,pim1,pim1)
3225  CALL bostr3(exe,pim2,pim2)
3226  CALL bostr3(exe,pr,pr)
3227 c partial width consists of phase space and amplitude
3228  IF (mnum.EQ.8) THEN
3229  CALL dampog(pt,pn,pim1,pim2,pipl,amplit,hv)
3230 c ELSEIF (mnum.EQ.0) THEN
3231 c CALL dampaa(pt,pn,pim1,pim2,pipl,amplit,hv)
3232  ELSE
3233  CALL damppk(mnum,pt,pn,pim1,pim2,pipl,amplit,hv)
3234  ENDIF
3235  IF (keyt.EQ.1.OR.keyt.EQ.2) THEN
3236 c the statistical factor for identical pi-s is cancelled with
3237 c two, for two modes of a1 decay namelly pi+pi-pi- and pi-pi0pi0
3238 #if defined (ALEPH)
3239 cam phspac=phspac*2.0
3240 cam phspac=phspac/2.
3241 #else
3242  phspac=phspac*2.0
3243  phspac=phspac/2.
3244 #endif
3245  ENDIF
3246  dgamt=1/(2.*amtau)*amplit*phspac
3247  END
3248  SUBROUTINE dampaa(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3249 c ----------------------------------------------------------------------
3250 * calculates differential cross section and polarimeter vector
3251 * for tau decay into a1, a1 decays next into rho+pi and rho into pi+pi.
3252 * all spin effects in the full decay chain are taken into account.
3253 * calculations done in tau rest frame with z-axis along neutrino moment
3254 * the routine is writen for zero neutrino mass.
3255 c
3256 c called by : dphsaa
3257 c ----------------------------------------------------------------------
3258  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3259  * ,ampiz,ampi,amro,gamro,ama1,gama1
3260  * ,amk,amkz,amkst,gamkst
3261 c
3262  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3263  * ,ampiz,ampi,amro,gamro,ama1,gama1
3264  * ,amk,amkz,amkst,gamkst
3265  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3266  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3267  COMMON /testa1/ keya1
3268  REAL hv(4),pt(4),pn(4),pim1(4),pim2(4),pipl(4)
3269  REAL paa(4),vec1(4),vec2(4)
3270  REAL pivec(4),piaks(4),hvm(4)
3271  COMPLEX bwign,hadcur(4),fpik
3272  DATA icont /1/
3273 c
3274 * f constants for a1, a1-rho-pi, and rho-pi-pi
3275 *
3276  DATA fpi /93.3e-3/
3277 * this inline funct. calculates the scalar part of the propagator
3278  bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
3279 c
3280 * four momentum of a1
3281  DO 10 i=1,4
3282  10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3283 * masses of a1, and of two pi-pairs which may form rho
3284  xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3285  xmro1 =sqrt(abs((pipl(4)+pim1(4))**2-(pipl(1)+pim1(1))**2
3286  $ -(pipl(2)+pim1(2))**2-(pipl(3)+pim1(3))**2))
3287  xmro2 =sqrt(abs((pipl(4)+pim2(4))**2-(pipl(1)+pim2(1))**2
3288  $ -(pipl(2)+pim2(2))**2-(pipl(3)+pim2(3))**2))
3289 * elements of hadron current
3290  prod1 =paa(4)*(pim1(4)-pipl(4))-paa(1)*(pim1(1)-pipl(1))
3291  $ -paa(2)*(pim1(2)-pipl(2))-paa(3)*(pim1(3)-pipl(3))
3292  prod2 =paa(4)*(pim2(4)-pipl(4))-paa(1)*(pim2(1)-pipl(1))
3293  $ -paa(2)*(pim2(2)-pipl(2))-paa(3)*(pim2(3)-pipl(3))
3294  DO 40 i=1,4
3295  vec1(i)= pim1(i)-pipl(i) -paa(i)*prod1/xmaa**2
3296  40 vec2(i)= pim2(i)-pipl(i) -paa(i)*prod2/xmaa**2
3297 * hadron current saturated with a1 and rho resonances
3298  IF (keya1.EQ.1) THEN
3299  fa1=9.87
3300  faropi=1.0
3301  fro2pi=1.0
3302  fnorm=fa1/sqrt(2.)*faropi*fro2pi
3303  DO 45 i=1,4
3304  hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gama1)
3305  $ *(cmplx(vec1(i))*amro**2*bwign(xmro1,amro,gamro)
3306  $ +cmplx(vec2(i))*amro**2*bwign(xmro2,amro,gamro))
3307  45 CONTINUE
3308  ELSE
3309  fnorm=2.0*sqrt(2.)/3.0/fpi
3310  gamax=gama1*gfun(xmaa**2)/gfun(ama1**2)
3311  DO 46 i=1,4
3312  hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gamax)
3313  $ *(cmplx(vec1(i))*fpik(xmro1)
3314  $ +cmplx(vec2(i))*fpik(xmro2))
3315  46 CONTINUE
3316  ENDIF
3317 c
3318 * calculate pi-vectors: vector and axial
3319  CALL clvec(hadcur,pn,pivec)
3320  CALL claxi(hadcur,pn,piaks)
3321  CALL clnut(hadcur,brakm,hvm)
3322 * spin independent part of decay diff-cross-sect. in tau rest frame
3323  brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3324  & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3325  amplit=(gfermi*ccabib)**2*brak/2.
3326 c the statistical factor for identical pi-s was cancelled with
3327 c two, for two modes of a1 decay namelly pi+pi-pi- and pi-pi0pi0
3328 c polarimeter vector in tau rest frame
3329  DO 90 i=1,3
3330  hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3331  & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3332 c hv is defined for tau- with gamma=b+hv*pol
3333  hv(i)=-hv(i)/brak
3334  90 CONTINUE
3335  END
3336 
3337  FUNCTION gfun(QKWA)
3338 c ****************************************************************
3339 c g-FUNCTION used to inroduce energy dependence in a1 width
3340 c ****************************************************************
3341  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3342  * ,ampiz,ampi,amro,gamro,ama1,gama1
3343  * ,amk,amkz,amkst,gamkst
3344 c
3345  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3346  * ,ampiz,ampi,amro,gamro,ama1,gama1
3347  * ,amk,amkz,amkst,gamkst
3348 c
3349  IF (qkwa.LT.(amro+ampi)**2) THEN
3350  gfun=4.1*(qkwa-9*ampiz**2)**3
3351  $ *(1.-3.3*(qkwa-9*ampiz**2)+5.8*(qkwa-9*ampiz**2)**2)
3352  ELSE
3353  gfun=qkwa*(1.623+10.38/qkwa-9.32/qkwa**2+0.65/qkwa**3)
3354  ENDIF
3355  END
3356  COMPLEX FUNCTION bwigs(S,M,G)
3357 c **********************************************************
3358 c p-wave breit-wigner for k*
3359 c **********************************************************
3360  REAL s,m,g
3361  REAL pi,pim,qs,qm,w,gs,mk
3362 #if defined (CLEO)
3363 c ajw: add k*-prim possibility:
3364  REAL pm, pg, pbeta
3365  COMPLEX bw,bwp
3366 #else
3367 #endif
3368  DATA init /0/
3369  p(a,b,c)=sqrt(abs(abs(((a+b-c)**2-4.*a*b)/4./a)
3370  $ +(((a+b-c)**2-4.*a*b)/4./a))/2.0)
3371 c ------------ parameters --------------------
3372  IF (init.EQ.0) THEN
3373  init=1
3374  pi=3.141592654
3375  pim=.139
3376  mk=.493667
3377 #if defined (CLEO)
3378 c ajw: add k*-prim possibility:
3379  pm = pkorb(1,16)
3380  pg = pkorb(2,16)
3381  pbeta = pkorb(3,16)
3382 #else
3383 #endif
3384 c ------- breit-wigner -----------------------
3385  ENDIF
3386 #if defined (ALEPH)
3387  IF (s.GT.(pim+mk)**2) THEN
3388 #endif
3389  qs=p(s,pim**2,mk**2)
3390  qm=p(m**2,pim**2,mk**2)
3391  w=sqrt(s)
3392  gs=g*(m/w)*(qs/qm)**3
3393 #if defined (CLEO)
3394  bw=m**2/cmplx(m**2-s,-m*gs)
3395  qpm=p(pm**2,pim**2,mk**2)
3396  g1=pg*(pm/w)*(qs/qpm)**3
3397  bwp=pm**2/cmplx(pm**2-s,-pm*g1)
3398  bwigs= (bw+pbeta*bwp)/(1+pbeta)
3399 #elif defined (ALEPH)
3400  ELSE
3401  gs=0.0
3402  ENDIF
3403  bwigs=m**2/cmplx(m**2-s,-m*gs)
3404 #else
3405  bwigs=m**2/cmplx(m**2-s,-m*gs)
3406 #endif
3407  RETURN
3408  END
3409  COMPLEX FUNCTION bwig(S,M,G)
3410 c **********************************************************
3411 c p-wave breit-wigner for rho
3412 c **********************************************************
3413  REAL s,m,g
3414  REAL pi,pim,qs,qm,w,gs
3415  DATA init /0/
3416 c ------------ parameters --------------------
3417  IF (init.EQ.0) THEN
3418  init=1
3419  pi=3.141592654
3420  pim=.139
3421 c ------- breit-wigner -----------------------
3422  ENDIF
3423  IF (s.GT.4.*pim**2) THEN
3424  qs=sqrt(abs(abs(s/4.-pim**2)+(s/4.-pim**2))/2.0)
3425  qm=sqrt(m**2/4.-pim**2)
3426  w=sqrt(s)
3427  gs=g*(m/w)*(qs/qm)**3
3428  ELSE
3429  gs=0.0
3430  ENDIF
3431  bwig=m**2/cmplx(m**2-s,-m*gs)
3432  RETURN
3433  END
3434  COMPLEX FUNCTION fpik(W)
3435 c **********************************************************
3436 c pion form factor
3437 c **********************************************************
3438  COMPLEX bwig
3439  REAL rom,rog,rom1,rog1,beta1,pi,pim,s,w
3440  EXTERNAL bwig
3441  DATA init /0/
3442 c
3443 c ------------ parameters --------------------
3444  IF (init.EQ.0 ) THEN
3445  init=1
3446  pi=3.141592654
3447  pim=.140
3448 #if defined (CLEO)
3449  rom=pkorb(1,9)
3450  rog=pkorb(2,9)
3451  rom1=pkorb(1,15)
3452  rog1=pkorb(2,15)
3453  beta1=pkorb(3,15)
3454 #else
3455  rom=0.773
3456  rog=0.145
3457  rom1=1.370
3458  rog1=0.510
3459  beta1=-0.145
3460 #endif
3461  ENDIF
3462 c -----------------------------------------------
3463  s=w**2
3464  fpik= (bwig(s,rom,rog)+beta1*bwig(s,rom1,rog1))
3465  & /(1+beta1)
3466  RETURN
3467  END
3468  FUNCTION fpirho(W)
3469 c **********************************************************
3470 c square of pion form factor
3471 c **********************************************************
3472  COMPLEX fpik
3473  fpirho=cabs(fpik(w))**2
3474  END
3475  SUBROUTINE clvec(HJ,PN,PIV)
3476 c ----------------------------------------------------------------------
3477 * calculates the "VECTOR TYPE" pi-vector piv
3478 * note that the neutrino mom. pn is assumed to be along z-axis
3479 c
3480 c called by : dampaa
3481 c ----------------------------------------------------------------------
3482  REAL piv(4),pn(4)
3483  COMPLEX hj(4),hn
3484 c
3485  hn= hj(4)*cmplx(pn(4))-hj(3)*cmplx(pn(3))
3486  hh= REAL(hj(4)*conjg(hj(4))-hj(3)*conjg(hj(3))
3487  $ -hj(2)*conjg(hj(2))-hj(1)*conjg(hj(1)))
3488  DO 10 i=1,4
3489  10 piv(i)=4.*REAL(hn*conjg(hj(i)))-2.*hh*pn(i)
3490  RETURN
3491  END
3492  SUBROUTINE claxi(HJ,PN,PIA)
3493 c ----------------------------------------------------------------------
3494 * calculates the "AXIAL TYPE" pi-vector pia
3495 * note that the neutrino mom. pn is assumed to be along z-axis
3496 c sign is chosen +/- for decay of tau +/- respectively
3497 c called by : dampaa, clnut
3498 c ----------------------------------------------------------------------
3499  COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3500  COMMON / idfc / idff
3501  REAL pia(4),pn(4)
3502  COMPLEX hj(4),hjc(4)
3503 c det2(i,j)=aimag(hj(i)*hjc(j)-hj(j)*hjc(i))
3504 c -- here was an error(zw, 21.11.1991)
3505  det2(i,j)=aimag(hjc(i)*hj(j)-hjc(j)*hj(i))
3506 c -- it was affecting sign of a_lr asymmetry in a1 decay.
3507 c -- note also collision of notation of gamma_va as defined in
3508 c -- tauola paper and j.h. kuhn and santamaria z. phys c 48 (1990) 445
3509 * -----------------------------------
3510  IF (ktom.EQ.1.OR.ktom.EQ.-1) THEN
3511  sign= idff/abs(idff)
3512  ELSEIF (ktom.EQ.2) THEN
3513  sign=-idff/abs(idff)
3514  ELSE
3515  print *, 'STOP IN CLAXI: KTOM=',ktom
3516  stop
3517  ENDIF
3518 c
3519  DO 10 i=1,4
3520  10 hjc(i)=conjg(hj(i))
3521  pia(1)= -2.*pn(3)*det2(2,4)+2.*pn(4)*det2(2,3)
3522  pia(2)= -2.*pn(4)*det2(1,3)+2.*pn(3)*det2(1,4)
3523  pia(3)= 2.*pn(4)*det2(1,2)
3524  pia(4)= 2.*pn(3)*det2(1,2)
3525 c all four indices are up so pia(3) and pia(4) have same sign
3526  DO 20 i=1,4
3527  20 pia(i)=pia(i)*sign
3528  END
3529  SUBROUTINE clnut(HJ,B,HV)
3530 c ----------------------------------------------------------------------
3531 * calculates the contribution by neutrino mass
3532 * note the tau is assumed to be at rest
3533 c
3534 c called by : dampaa
3535 c ----------------------------------------------------------------------
3536  COMPLEX hj(4)
3537  REAL hv(4),p(4)
3538  DATA p /3*0.,1.0/
3539 c
3540  CALL claxi(hj,p,hv)
3541  b=REAL( HJ(4)*AIMAG(HJ(4)) - HJ(3)*AIMAG(HJ(3)) & - HJ(2)*AIMAG(HJ(2)) - HJ(1)*AIMAG(HJ(1)) )
3542  RETURN
3543  END
3544  SUBROUTINE dampog(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3545 c ----------------------------------------------------------------------
3546 * calculates differential cross section and polarimeter vector
3547 * for tau decay into a1, a1 decays next into rho+pi and rho into pi+pi.
3548 * all spin effects in the full decay chain are taken into account.
3549 * calculations done in tau rest frame with z-axis along neutrino moment
3550 * the routine is writen for zero neutrino mass.
3551 c
3552 #if defined (ALEPH)
3553 c called by : dphtre
3554 #else
3555 c called by : dphsaa
3556 #endif
3557 c ----------------------------------------------------------------------
3558  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3559  * ,ampiz,ampi,amro,gamro,ama1,gama1
3560  * ,amk,amkz,amkst,gamkst
3561 c
3562  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3563  * ,ampiz,ampi,amro,gamro,ama1,gama1
3564  * ,amk,amkz,amkst,gamkst
3565  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3566  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3567  COMMON /testa1/ keya1
3568  REAL hv(4),pt(4),pn(4),pim1(4),pim2(4),pipl(4)
3569  REAL paa(4),vec1(4),vec2(4)
3570  REAL pivec(4),piaks(4),hvm(4)
3571  COMPLEX bwign,hadcur(4),fnorm,formom
3572  DATA icont /1/
3573 * this inline funct. calculates the scalar part of the propagator
3574 #if defined (CLEO)
3575 c ajwmod to satisfy compiler, comment out this unused function.
3576 #else
3577  bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
3578 #endif
3579 c
3580 * four momentum of a1
3581  DO 10 i=1,4
3582  vec1(i)=0.0
3583  vec2(i)=0.0
3584  hv(i) =0.0
3585  10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3586  vec1(1)=1.0
3587 * masses of a1, and of two pi-pairs which may form rho
3588  xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3589  xmom =sqrt(abs( (pim2(4)+pipl(4))**2-(pim2(3)+pipl(3))**2
3590  $ -(pim2(2)+pipl(2))**2-(pim2(1)+pipl(1))**2 ))
3591  xmro2 =(pipl(1))**2 +(pipl(2))**2 +(pipl(3))**2
3592 * elements of hadron current
3593  prod1 =vec1(1)*pipl(1)
3594  prod2 =vec2(2)*pipl(2)
3595  p12 =pim1(4)*pim2(4)-pim1(1)*pim2(1)
3596  $ -pim1(2)*pim2(2)-pim1(3)*pim2(3)
3597  p1pl =pim1(4)*pipl(4)-pim1(1)*pipl(1)
3598  $ -pim1(2)*pipl(2)-pim1(3)*pipl(3)
3599  p2pl =pipl(4)*pim2(4)-pipl(1)*pim2(1)
3600  $ -pipl(2)*pim2(2)-pipl(3)*pim2(3)
3601  DO 40 i=1,3
3602  vec1(i)= (vec1(i)-prod1/xmro2*pipl(i))
3603  40 CONTINUE
3604  gnorm=sqrt(vec1(1)**2+vec1(2)**2+vec1(3)**2)
3605  DO 41 i=1,3
3606  vec1(i)= vec1(i)/gnorm
3607  41 CONTINUE
3608  vec2(1)=(vec1(2)*pipl(3)-vec1(3)*pipl(2))/sqrt(xmro2)
3609  vec2(2)=(vec1(3)*pipl(1)-vec1(1)*pipl(3))/sqrt(xmro2)
3610  vec2(3)=(vec1(1)*pipl(2)-vec1(2)*pipl(1))/sqrt(xmro2)
3611  p1vec1 =pim1(4)*vec1(4)-pim1(1)*vec1(1)
3612  $ -pim1(2)*vec1(2)-pim1(3)*vec1(3)
3613  p2vec1 =vec1(4)*pim2(4)-vec1(1)*pim2(1)
3614  $ -vec1(2)*pim2(2)-vec1(3)*pim2(3)
3615  p1vec2 =pim1(4)*vec2(4)-pim1(1)*vec2(1)
3616  $ -pim1(2)*vec2(2)-pim1(3)*vec2(3)
3617  p2vec2 =vec2(4)*pim2(4)-vec2(1)*pim2(1)
3618  $ -vec2(2)*pim2(2)-vec2(3)*pim2(3)
3619 * hadron current
3620  fnorm=formom(xmaa,xmom)
3621  brak=0.0
3622  DO 120 jj=1,2
3623  DO 45 i=1,4
3624  IF (jj.EQ.1) THEN
3625  hadcur(i) = fnorm *(
3626  $ vec1(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3627  $ -pim2(i)*(p2vec1*p1pl-p1vec1*p2pl)
3628  $ +pipl(i)*(p2vec1*p12 -p1vec1*(ampi**2+p2pl)) )
3629  ELSE
3630  hadcur(i) = fnorm *(
3631  $ vec2(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3632  $ -pim2(i)*(p2vec2*p1pl-p1vec2*p2pl)
3633  $ +pipl(i)*(p2vec2*p12 -p1vec2*(ampi**2+p2pl)) )
3634  ENDIF
3635  45 CONTINUE
3636 c
3637 * calculate pi-vectors: vector and axial
3638  CALL clvec(hadcur,pn,pivec)
3639  CALL claxi(hadcur,pn,piaks)
3640  CALL clnut(hadcur,brakm,hvm)
3641 * spin independent part of decay diff-cross-sect. in tau rest frame
3642  brak=brak+(gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3643  & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3644  DO 90 i=1,3
3645  hv(i)=hv(i)-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3646  & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3647  90 CONTINUE
3648 c hv is defined for tau- with gamma=b+hv*pol
3649  120 CONTINUE
3650  amplit=(gfermi*ccabib)**2*brak/2.
3651 c the statistical factor for identical pi-s was cancelled with
3652 c two, for two modes of a1 decay namelly pi+pi-pi- and pi-pi0pi0
3653 c polarimeter vector in tau rest frame
3654  DO 91 i=1,3
3655  hv(i)=-hv(i)/brak
3656  91 CONTINUE
3657 
3658  END
3659  SUBROUTINE damppk(MNUM,PT,PN,PIM1,PIM2,PIM3,AMPLIT,HV)
3660 c ----------------------------------------------------------------------
3661 * calculates differential cross section and polarimeter vector
3662 * for tau decay into k k pi, k pi pi.
3663 * all spin effects in the full decay chain are taken into account.
3664 * calculations done in tau rest frame with z-axis along neutrino moment
3665 c mnum decay mode identifier.
3666 c
3667 #if defined (ALEPH)
3668 c called by : dphtre
3669 #else
3670 c called by : dphsaa
3671 #endif
3672 c ----------------------------------------------------------------------
3673  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3674  * ,ampiz,ampi,amro,gamro,ama1,gama1
3675  * ,amk,amkz,amkst,gamkst
3676 c
3677  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3678  * ,ampiz,ampi,amro,gamro,ama1,gama1
3679  * ,amk,amkz,amkst,gamkst
3680  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3681  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3682  REAL hv(4),pt(4),pn(4),pim1(4),pim2(4),pim3(4)
3683  REAL paa(4),vec1(4),vec2(4),vec3(4),vec4(4),vec5(4)
3684  REAL pivec(4),piaks(4),hvm(4)
3685  REAL fnorm(0:7),coef(1:5,0:7)
3686  COMPLEX hadcur(4),form1,form2,form3,form4,form5,uroj
3687 #if defined (CLEO)
3688  COMPLEX f1,f2,f3,f4,f5
3689 #endif
3690  EXTERNAL form1,form2,form3,form4,form5
3691  DATA pi /3.141592653589793238462643/
3692  DATA icont /0/
3693 c
3694  DATA fpi /93.3e-3/
3695  IF (icont.EQ.0) THEN
3696  icont=1
3697  uroj=cmplx(0.0,1.0)
3698  dwapi0=sqrt(2.0)
3699  fnorm(0)=ccabib/fpi
3700  fnorm(1)=ccabib/fpi
3701  fnorm(2)=ccabib/fpi
3702  fnorm(3)=ccabib/fpi
3703  fnorm(4)=scabib/fpi/dwapi0
3704  fnorm(5)=scabib/fpi
3705  fnorm(6)=scabib/fpi
3706  fnorm(7)=ccabib/fpi
3707 c
3708  coef(1,0)= 2.0*sqrt(2.)/3.0
3709  coef(2,0)=-2.0*sqrt(2.)/3.0
3710 #if defined (CLEO)
3711 c ajw 2/98: add in the d-wave and i=0 3pi substructure:
3712  coef(3,0)= 2.0*sqrt(2.)/3.0
3713 #else
3714  coef(3,0)= 0.0
3715 #endif
3716  coef(4,0)= fpi
3717  coef(5,0)= 0.0
3718 c
3719  coef(1,1)=-sqrt(2.)/3.0
3720  coef(2,1)= sqrt(2.)/3.0
3721  coef(3,1)= 0.0
3722  coef(4,1)= fpi
3723  coef(5,1)= sqrt(2.)
3724 c
3725  coef(1,2)=-sqrt(2.)/3.0
3726  coef(2,2)= sqrt(2.)/3.0
3727  coef(3,2)= 0.0
3728  coef(4,2)= 0.0
3729  coef(5,2)=-sqrt(2.)
3730 c
3731 #if defined (CLEO)
3732 c ajw 11/97: add in the k*-prim-s, ala finkemeier&mirkes
3733  coef(1,3)= 1./3.
3734  coef(2,3)=-2./3.
3735  coef(3,3)= 2./3.
3736 #else
3737  coef(1,3)= 0.0
3738  coef(2,3)=-1.0
3739  coef(3,3)= 0.0
3740 #endif
3741  coef(4,3)= 0.0
3742  coef(5,3)= 0.0
3743 c
3744  coef(1,4)= 1.0/sqrt(2.)/3.0
3745  coef(2,4)=-1.0/sqrt(2.)/3.0
3746  coef(3,4)= 0.0
3747  coef(4,4)= 0.0
3748  coef(5,4)= 0.0
3749 c
3750  coef(1,5)=-sqrt(2.)/3.0
3751  coef(2,5)= sqrt(2.)/3.0
3752  coef(3,5)= 0.0
3753  coef(4,5)= 0.0
3754  coef(5,5)=-sqrt(2.)
3755 c
3756 #if defined (CLEO)
3757 c ajw 11/97: add in the k*-prim-s, ala finkemeier&mirkes
3758  coef(1,6)= 1./3.
3759  coef(2,6)=-2./3.
3760  coef(3,6)= 2./3.
3761 #else
3762  coef(1,6)= 0.0
3763  coef(2,6)=-1.0
3764  coef(3,6)= 0.0
3765 #endif
3766  coef(4,6)= 0.0
3767  coef(5,6)=-2.0
3768 c
3769  coef(1,7)= 0.0
3770  coef(2,7)= 0.0
3771  coef(3,7)= 0.0
3772  coef(4,7)= 0.0
3773  coef(5,7)=-sqrt(2.0/3.0)
3774 c
3775  ENDIF
3776 c
3777  DO 10 i=1,4
3778  10 paa(i)=pim1(i)+pim2(i)+pim3(i)
3779  xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3780  xmro1 =sqrt(abs((pim3(4)+pim2(4))**2-(pim3(1)+pim2(1))**2
3781  $ -(pim3(2)+pim2(2))**2-(pim3(3)+pim2(3))**2))
3782  xmro2 =sqrt(abs((pim3(4)+pim1(4))**2-(pim3(1)+pim1(1))**2
3783  $ -(pim3(2)+pim1(2))**2-(pim3(3)+pim1(3))**2))
3784  xmro3 =sqrt(abs((pim1(4)+pim2(4))**2-(pim1(1)+pim2(1))**2
3785  $ -(pim1(2)+pim2(2))**2-(pim1(3)+pim2(3))**2))
3786 * elements of hadron current
3787  prod1 =paa(4)*(pim2(4)-pim3(4))-paa(1)*(pim2(1)-pim3(1))
3788  $ -paa(2)*(pim2(2)-pim3(2))-paa(3)*(pim2(3)-pim3(3))
3789  prod2 =paa(4)*(pim3(4)-pim1(4))-paa(1)*(pim3(1)-pim1(1))
3790  $ -paa(2)*(pim3(2)-pim1(2))-paa(3)*(pim3(3)-pim1(3))
3791  prod3 =paa(4)*(pim1(4)-pim2(4))-paa(1)*(pim1(1)-pim2(1))
3792  $ -paa(2)*(pim1(2)-pim2(2))-paa(3)*(pim1(3)-pim2(3))
3793  DO 40 i=1,4
3794  vec1(i)= pim2(i)-pim3(i) -paa(i)*prod1/xmaa**2
3795  vec2(i)= pim3(i)-pim1(i) -paa(i)*prod2/xmaa**2
3796  vec3(i)= pim1(i)-pim2(i) -paa(i)*prod3/xmaa**2
3797  40 vec4(i)= pim1(i)+pim2(i)+pim3(i)
3798  CALL prod5(pim1,pim2,pim3,vec5)
3799 * hadron current
3800 c be aware that sign of vec2 is opposite to sign of vec1 in a1 case
3801 #if defined (CLEO)
3802 c rationalize this code:
3803  f1 = cmplx(coef(1,mnum))*form1(mnum,xmaa**2,xmro1**2,xmro2**2)
3804  f2 = cmplx(coef(2,mnum))*form2(mnum,xmaa**2,xmro2**2,xmro1**2)
3805  f3 = cmplx(coef(3,mnum))*form3(mnum,xmaa**2,xmro3**2,xmro1**2)
3806  f4 = (-1.0*uroj)*
3807  $cmplx(coef(4,mnum))*form4(mnum,xmaa**2,xmro1**2,xmro2**2,xmro3**2)
3808  f5 = (-1.0)*uroj/4.0/pi**2/fpi**2*
3809  $ cmplx(coef(5,mnum))*form5(mnum,xmaa**2,xmro1**2,xmro2**2)
3810 
3811  DO 45 i=1,4
3812  hadcur(i)= cmplx(fnorm(mnum)) * (
3813  $ cmplx(vec1(i))*f1+cmplx(vec2(i))*f2+cmplx(vec3(i))*f3+
3814  $ cmplx(vec4(i))*f4+cmplx(vec5(i))*f5)
3815  45 CONTINUE
3816 #else
3817  DO 45 i=1,4
3818  hadcur(i)= cmplx(fnorm(mnum)) * (
3819  $cmplx(vec1(i)*coef(1,mnum))*form1(mnum,xmaa**2,xmro1**2,xmro2**2)+
3820  $cmplx(vec2(i)*coef(2,mnum))*form2(mnum,xmaa**2,xmro2**2,xmro1**2)+
3821  $cmplx(vec3(i)*coef(3,mnum))*form3(mnum,xmaa**2,xmro3**2,xmro1**2)+
3822  *(-1.0*uroj)*
3823  $cmplx(vec4(i)*coef(4,mnum))*form4(mnum,xmaa**2,xmro1**2,
3824  $ xmro2**2,xmro3**2) +
3825  $(-1.0)*uroj/4.0/pi**2/fpi**2*
3826  $cmplx(vec5(i)*coef(5,mnum))*form5(mnum,xmaa**2,xmro1**2,xmro2**2))
3827  45 CONTINUE
3828 #endif
3829 c
3830 * calculate pi-vectors: vector and axial
3831  CALL clvec(hadcur,pn,pivec)
3832  CALL claxi(hadcur,pn,piaks)
3833  CALL clnut(hadcur,brakm,hvm)
3834 * spin independent part of decay diff-cross-sect. in tau rest frame
3835  brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3836  & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3837  amplit=(gfermi)**2*brak/2.
3838  IF (mnum.GE.9) THEN
3839  print *, 'MNUM=',mnum
3840  znak=-1.0
3841  xm1=0.0
3842  xm2=0.0
3843  xm3=0.0
3844  DO 77 k=1,4
3845  IF (k.EQ.4) znak=1.0
3846  xm1=znak*pim1(k)**2+xm1
3847  xm2=znak*pim2(k)**2+xm2
3848  xm3=znak*pim3(k)**2+xm3
3849  77 print *, 'PIM1=',pim1(k),'PIM2=',pim2(k),'PIM3=',pim3(k)
3850  print *, 'XM1=',sqrt(xm1),'XM2=',sqrt(xm2),'XM3=',sqrt(xm3)
3851  print *, '************************************************'
3852  ENDIF
3853 c polarimeter vector in tau rest frame
3854  DO 90 i=1,3
3855  hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3856  & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3857 c hv is defined for tau- with gamma=b+hv*pol
3858  hv(i)=-hv(i)/brak
3859  90 CONTINUE
3860  END
3861  SUBROUTINE prod5(P1,P2,P3,PIA)
3862 c ----------------------------------------------------------------------
3863 c external product of p1, p2, p3 4-momenta.
3864 c sign is chosen +/- for decay of tau +/- respectively
3865 c called by : dampaa, clnut
3866 c ----------------------------------------------------------------------
3867  COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3868  COMMON / idfc / idff
3869  REAL pia(4),p1(4),p2(4),p3(4)
3870  det2(i,j)=p1(i)*p2(j)-p2(i)*p1(j)
3871 * -----------------------------------
3872  IF (ktom.EQ.1.OR.ktom.EQ.-1) THEN
3873  sign= idff/abs(idff)
3874  ELSEIF (ktom.EQ.2) THEN
3875  sign=-idff/abs(idff)
3876  ELSE
3877  print *, 'STOP IN PROD5: KTOM=',ktom
3878  stop
3879  ENDIF
3880 c
3881 c epsilon( p1(1), p2(2), p3(3), (4) ) = 1
3882 c
3883  pia(1)= -p3(3)*det2(2,4)+p3(4)*det2(2,3)+p3(2)*det2(3,4)
3884  pia(2)= -p3(4)*det2(1,3)+p3(3)*det2(1,4)-p3(1)*det2(3,4)
3885  pia(3)= p3(4)*det2(1,2)-p3(2)*det2(1,4)+p3(1)*det2(2,4)
3886  pia(4)= p3(3)*det2(1,2)-p3(2)*det2(1,3)+p3(1)*det2(2,3)
3887 c all four indices are up so pia(3) and pia(4) have same sign
3888  DO 20 i=1,4
3889  20 pia(i)=pia(i)*sign
3890  END
3891 
3892  SUBROUTINE dexnew(MODE,ISGN,POL,PNU,PAA,PNPI,JNPI)
3893 c ----------------------------------------------------------------------
3894 * this simulates tau decay in tau rest frame
3895 * into nu a1, next a1 decays into rho pi and finally rho into pi pi.
3896 * output four momenta: pnu tauneutrino,
3897 #if defined (ALEPH)
3898 * paa hadron 4-vector
3899 * pnpi final state particles
3900 * jnpi decay type
3901 #else
3902 * paa a1
3903 * pim1 pion minus(or pi0) 1 (for tau minus)
3904 * pim2 pion minus(or pi0) 2
3905 * pipl pion plus(or pi-)
3906 * (pipl,pim1) form a rho
3907 #endif
3908 c ----------------------------------------------------------------------
3909  COMMON / inout / inut,iout
3910  REAL pol(4),hv(4),paa(4),pnu(4),pnpi(4,9),rn(1)
3911  DATA iwarm/0/
3912 c
3913  IF(mode.EQ.-1) THEN
3914 c ===================
3915  iwarm=1
3916  CALL dadnew( -1,isgn,hv,pnu,paa,pnpi,jdumm)
3917 #if defined (ALEPH)
3918 cc CALL hbook1(816,'WEIGHT DISTRIBUTION DEXNEW $',100,-2.,2.)
3919 #else
3920 cc CALL hbook1(816,'WEIGHT DISTRIBUTION DEXAA $',100,-2.,2.)
3921 #endif
3922 c
3923  ELSEIF(mode.EQ. 0) THEN
3924 * =======================
3925  300 CONTINUE
3926  IF(iwarm.EQ.0) goto 902
3927  CALL dadnew( 0,isgn,hv,pnu,paa,pnpi,jnpi)
3928  wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
3929 cc CALL hfill(816,wt)
3930  CALL ranmar(rn,1)
3931  IF(rn(1).GT.wt) goto 300
3932 c
3933  ELSEIF(mode.EQ. 1) THEN
3934 * =======================
3935  CALL dadnew( 1,isgn,hv,pnu,paa,pnpi,jdumm)
3936 cc CALL hprint(816)
3937  ENDIF
3938 c =====
3939  RETURN
3940  902 WRITE(iout, 9020)
3941  9020 FORMAT(' ----- DEXNEW: LACK OF INITIALISATION')
3942  stop
3943  END
3944  SUBROUTINE dadnew(MODE,ISGN,HV,PNU,PWB,PNPI,JNPI)
3945 c ----------------------------------------------------------------------
3946  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3947  * ,ampiz,ampi,amro,gamro,ama1,gama1
3948  * ,amk,amkz,amkst,gamkst
3949 c
3950  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3951  * ,ampiz,ampi,amro,gamro,ama1,gama1
3952  * ,amk,amkz,amkst,gamkst
3953  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3954  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3955  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
3956  REAL*4 gampmc ,gamper
3957 #if defined (ALEPH)
3958 #else
3959  COMMON / inout / inut,iout
3960 #endif
3961  parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
3962 #if defined (ALEPH)
3963  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3964 #else
3965  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3966 #endif
3967  & ,names
3968  CHARACTER names(nmode)*31
3969 #if defined (ALEPH)
3970  COMMON / inout / inut,iout
3971 #endif
3972 
3973  REAL*4 pnu(4),pwb(4),pnpi(4,9),hv(4),hhv(4)
3974  REAL*4 pdum1(4),pdum2(4),pdumi(4,9)
3975  REAL*4 rrr(3)
3976  REAL*4 wtmax(nmode)
3977  REAL*8 swt(nmode),sswt(nmode)
3978  dimension nevraw(nmode),nevovr(nmode),nevacc(nmode)
3979 c
3980  DATA pi /3.141592653589793238462643/
3981  DATA iwarm/0/
3982 c
3983  IF(mode.EQ.-1) THEN
3984 c ===================
3985 c -- at the moment only two decay modes of multipions have m. elem
3986  nmod=nmode
3987  iwarm=1
3988 c print 7003
3989  DO 1 jnpi=1,nmod
3990  nevraw(jnpi)=0
3991  nevacc(jnpi)=0
3992  nevovr(jnpi)=0
3993  swt(jnpi)=0
3994  sswt(jnpi)=0
3995  wtmax(jnpi)=-1.
3996 #if defined (CePeCe)
3997  DO i=1,500
3998 #elif defined (ALEPH)
3999  DO i=1,500
4000 #else
4001 c for 4pi phase space, need lots more trials at initialization,
4002 c or use the wtmax determined with many trials for default model:
4003  ntrials = 500
4004  IF (jnpi.LE.nm4) THEN
4005  a = pkorb(3,37+jnpi)
4006  IF (a.LT.0.) THEN
4007  ntrials = 10000
4008  ELSE
4009  ntrials = 0
4010  wtmax(jnpi)=a
4011  END IF
4012  END IF
4013  DO i=1,ntrials
4014 #endif
4015  IF (jnpi.LE.0) THEN
4016  goto 903
4017  ELSEIF(jnpi.LE.nm4) THEN
4018  CALL dph4pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4019  ELSEIF(jnpi.LE.nm4+nm5) THEN
4020  CALL dph5pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4021  ELSEIF(jnpi.LE.nm4+nm5+nm6) THEN
4022  CALL dphnpi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4023  ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3) THEN
4024  inum=jnpi-nm4-nm5-nm6
4025  CALL dphspk(wt,hv,pdum1,pdum2,pdumi,inum)
4026  ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2) THEN
4027  inum=jnpi-nm4-nm5-nm6-nm3
4028  CALL dphsrk(wt,hv,pdum1,pdum2,pdumi,inum)
4029  ELSE
4030  goto 903
4031  ENDIF
4032  IF(wt.GT.wtmax(jnpi)/1.2) wtmax(jnpi)=wt*1.2
4033  ENDDO
4034 #if defined (CePeCe)
4035 #elif defined (ALEPH)
4036 #else
4037 c print *,' DADNEW JNPI,NTRIALS,WTMAX =',jnpi,ntrials,wtmax(jnpi)
4038 #endif
4039 c CALL hbook1(801,'WEIGHT DISTRIBUTION DADNPI $',100,0.,2.,.0)
4040 c print 7004,wtmax(jnpi)
4041 1 CONTINUE
4042  WRITE(iout,7005)
4043 c
4044  ELSEIF(mode.EQ. 0) THEN
4045 c =======================
4046  IF(iwarm.EQ.0) goto 902
4047 c
4048 300 CONTINUE
4049  IF (jnpi.LE.0) THEN
4050  goto 903
4051  ELSEIF(jnpi.LE.nm4) THEN
4052  CALL dph4pi(wt,hhv,pnu,pwb,pnpi,jnpi)
4053  ELSEIF(jnpi.LE.nm4+nm5) THEN
4054  CALL dph5pi(wt,hhv,pnu,pwb,pnpi,jnpi)
4055  ELSEIF(jnpi.LE.nm4+nm5+nm6) THEN
4056  CALL dphnpi(wt,hhv,pnu,pwb,pnpi,jnpi)
4057  ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3) THEN
4058  inum=jnpi-nm4-nm5-nm6
4059  CALL dphspk(wt,hhv,pnu,pwb,pnpi,inum)
4060  ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2) THEN
4061  inum=jnpi-nm4-nm5-nm6-nm3
4062  CALL dphsrk(wt,hhv,pnu,pwb,pnpi,inum)
4063  ELSE
4064  goto 903
4065  ENDIF
4066  DO i=1,4
4067  hv(i)=-isgn*hhv(i)
4068  ENDDO
4069 c CALL hfill(801,wt/wtmax(jnpi))
4070  nevraw(jnpi)=nevraw(jnpi)+1
4071  swt(jnpi)=swt(jnpi)+wt
4072 #if defined (ALEPH)
4073  sswt(jnpi)=sswt(jnpi)+wt**2
4074 #else
4075 cccm.s.>>>>>>
4076 cc sswt(jnpi)=sswt(jnpi)+wt**2
4077  sswt(jnpi)=sswt(jnpi)+dble(wt)**2
4078 cccm.s.<<<<<<
4079 #endif
4080  CALL ranmar(rrr,3)
4081  rn=rrr(1)
4082  IF(wt.GT.wtmax(jnpi)) nevovr(jnpi)=nevovr(jnpi)+1
4083  IF(rn*wtmax(jnpi).GT.wt) goto 300
4084 c rotations to basic tau rest frame
4085  costhe=-1.+2.*rrr(2)
4086  thet=acos(costhe)
4087  phi =2*pi*rrr(3)
4088  CALL rotor2(thet,pnu,pnu)
4089  CALL rotor3( phi,pnu,pnu)
4090  CALL rotor2(thet,pwb,pwb)
4091  CALL rotor3( phi,pwb,pwb)
4092  CALL rotor2(thet,hv,hv)
4093  CALL rotor3( phi,hv,hv)
4094  nd=mulpik(jnpi)
4095  DO 301 i=1,nd
4096  CALL rotor2(thet,pnpi(1,i),pnpi(1,i))
4097  CALL rotor3( phi,pnpi(1,i),pnpi(1,i))
4098 301 CONTINUE
4099  nevacc(jnpi)=nevacc(jnpi)+1
4100 c
4101  ELSEIF(mode.EQ. 1) THEN
4102 c =======================
4103  DO 500 jnpi=1,nmod
4104  IF(nevraw(jnpi).EQ.0) goto 500
4105  pargam=swt(jnpi)/float(nevraw(jnpi)+1)
4106  error=0
4107  IF(nevraw(jnpi).NE.0)
4108  & error=sqrt(sswt(jnpi)/swt(jnpi)**2-1./float(nevraw(jnpi)))
4109  rat=pargam/gamel
4110  WRITE(iout, 7010) names(jnpi),
4111  & nevraw(jnpi),nevacc(jnpi),nevovr(jnpi),pargam,rat,error
4112 cc CALL hprint(801)
4113  gampmc(8+jnpi-1)=rat
4114  gamper(8+jnpi-1)=error
4115 cam nevdec(8+jnpi-1)=nevacc(jnpi)
4116  500 CONTINUE
4117  ENDIF
4118 c =====
4119  RETURN
4120  7003 FORMAT(///1x,15(5h*****)
4121  $ /,' *', 25x,'******** DADNEW INITIALISATION ********',9x,1h*
4122  $ )
4123  7004 FORMAT(' *',e20.5,5x,'WTMAX = MAXIMUM WEIGHT ',9x,1h*/)
4124  7005 FORMAT(
4125  $ /,1x,15(5h*****)/)
4126  7010 FORMAT(///1x,15(5h*****)
4127  $ /,' *', 25x,'******** DADNEW FINAL REPORT ******** ',9x,1h*
4128  $ /,' *', 25x,'CHANNEL:',a31 ,9x,1h*
4129  $ /,' *',i20 ,5x,'NEVRAW = NO. OF DECAYS TOTAL ',9x,1h*
4130  $ /,' *',i20 ,5x,'NEVACC = NO. OF DECAYS ACCEPTED ',9x,1h*
4131  $ /,' *',i20 ,5x,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
4132  $ /,' *',e20.5,5x,'PARTIAL WTDTH IN GEV UNITS ',9x,1h*
4133  $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
4134  $ /,' *',f20.8,5x,'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
4135  $ /,1x,15(5h*****)/)
4136  902 WRITE(iout, 9020)
4137  9020 FORMAT(' ----- DADNEW: LACK OF INITIALISATION')
4138  stop
4139  903 WRITE(iout, 9030) jnpi,mode
4140  9030 FORMAT(' ----- DADNEW: WRONG JNPI',2i5)
4141  stop
4142  END
4143 
4144 
4145  SUBROUTINE dph4pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4146 c ----------------------------------------------------------------------
4147 #if defined (ALEPH)
4148 * it simulates 4pi decay in tau rest frame with
4149 * z-axis along 4pi momentum
4150 #else
4151 * it simulates a1 decay in tau rest frame with
4152 * z-axis along a1 momentum
4153 #endif
4154 c ----------------------------------------------------------------------
4155  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4156  * ,ampiz,ampi,amro,gamro,ama1,gama1
4157  * ,amk,amkz,amkst,gamkst
4158 c
4159  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4160  * ,ampiz,ampi,amro,gamro,ama1,gama1
4161  * ,amk,amkz,amkst,gamkst
4162  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4163  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
4164 #if defined (ALEPH)
4165  parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4166  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4167  & ,names
4168  CHARACTER names(nmode)*31
4169 #else
4170 #endif
4171  REAL hv(4),pt(4),pn(4),paa(4),pim1(4),pim2(4),pipl(4),pmult(4,9)
4172  REAL pr(4),piz(4)
4173  REAL*4 rrr(9)
4174  REAL*8 uu,ff,ff1,ff2,ff3,ff4,gg1,gg2,gg3,gg4,rr
4175  DATA pi /3.141592653589793238462643/
4176  DATA icont /0/
4177  xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
4178 c amro, gamro is only a PARAMETER for geting hight efficiency
4179 c
4180 c three body phase space normalised as in bjorken-drell
4181 c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
4182  phspac=1./2**23/pi**11
4183  phsp=1./2**5/pi**2
4184 #if defined (ALEPH)
4185 c init decay mode jnpi
4186  amp1=dcdmas(idffin(1,jnpi))
4187  amp2=dcdmas(idffin(2,jnpi))
4188  amp3=dcdmas(idffin(3,jnpi))
4189  amp4=dcdmas(idffin(4,jnpi))
4190 #endif
4191  IF (jnpi.EQ.1) THEN
4192  prez=0.7
4193 #if defined (CePeCe)
4194  amp1=ampi
4195  amp2=ampi
4196  amp3=ampi
4197  amp4=ampiz
4198  amrx=0.782
4199  gamrx=0.0084
4200 #elif defined (CLEO)
4201  amp1=ampi
4202  amp2=ampi
4203  amp3=ampi
4204  amp4=ampiz
4205  amrx=pkorb(1,14)
4206  gamrx=pkorb(2,14)
4207 c ajw: cant simply change amrop, etc, here!
4208 c choice is a by-hand tuning/optimization, no simple relationship
4209 c to actual resonance masses(accd to z.was).
4210 c what matters in the end is what you put in formf/curr .
4211 #else
4212  amrx=0.782
4213  gamrx=0.0084
4214 #endif
4215  amrop =1.2
4216  gamrop=.46
4217  ELSE
4218  prez=0.0
4219 #if defined (ALEPH)
4220 #else
4221  amp1=ampiz
4222  amp2=ampiz
4223  amp3=ampiz
4224  amp4=ampi
4225 #endif
4226  amrx=1.4
4227  gamrx=.6
4228  amrop =amrx
4229  gamrop=gamrx
4230 
4231  ENDIF
4232 #if defined (ALEPH)
4233 ! 07.06.96 here was an error in the type of variable.
4234  rrdum=0.3
4235  CALL choice(100+jnpi,rrdum,ichan,prob1,prob2,prob3,
4236 #else
4237  rrb=0.3
4238  CALL choice(100+jnpi,rrb,ichan,prob1,prob2,prob3,
4239 #endif
4240  $ amrop,gamrop,amrx,gamrx,amrb,gamrb)
4241  prez=prob1+prob2
4242 c tau momentum
4243  pt(1)=0.
4244  pt(2)=0.
4245  pt(3)=0.
4246  pt(4)=amtau
4247 c
4248  CALL ranmar(rrr,9)
4249 c
4250 * masses of 4, 3 and 2 pi systems
4251 c 3 pi with sampling for resonance
4252 cam
4253  rr1=rrr(6)
4254  ams1=(amp1+amp2+amp3+amp4)**2
4255  ams2=(amtau-amnuta)**2
4256  alp1=atan((ams1-amrop**2)/amrop/gamrop)
4257  alp2=atan((ams2-amrop**2)/amrop/gamrop)
4258  alp=alp1+rr1*(alp2-alp1)
4259  am4sq =amrop**2+amrop*gamrop*tan(alp)
4260  am4 =sqrt(am4sq)
4261  phspac=phspac*
4262  $ ((am4sq-amrop**2)**2+(amrop*gamrop)**2)/(amrop*gamrop)
4263  phspac=phspac*(alp2-alp1)
4264 
4265 c
4266  rr1=rrr(1)
4267  ams1=(amp2+amp3+amp4)**2
4268  ams2=(am4-amp1)**2
4269  IF (rrr(9).GT.prez) THEN
4270  am3sq=ams1+ rr1*(ams2-ams1)
4271  am3 =sqrt(am3sq)
4272 c --- this part of jacobian will be recovered later
4273  ff1=ams2-ams1
4274  ELSE
4275 * phase space with sampling for omega resonance,
4276  alp1=atan((ams1-amrx**2)/amrx/gamrx)
4277  alp2=atan((ams2-amrx**2)/amrx/gamrx)
4278  alp=alp1+rr1*(alp2-alp1)
4279  am3sq =amrx**2+amrx*gamrx*tan(alp)
4280  am3 =sqrt(am3sq)
4281 c --- this part of the jacobian will be recovered later ---------------
4282  ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4283  ff1=ff1*(alp2-alp1)
4284  ENDIF
4285 c mass of 2
4286  rr2=rrr(2)
4287  ams1=(amp3+amp4)**2
4288  ams2=(am3-amp2)**2
4289 * flat phase space;
4290  am2sq=ams1+ rr2*(ams2-ams1)
4291  am2 =sqrt(am2sq)
4292 c --- this part of jacobian will be recovered later
4293  ff2=(ams2-ams1)
4294 * 2 restframe, define piz and pipl
4295 #if defined (ALEPH)
4296  enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4297  enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4298  ppi= enq1**2-amp3**2
4299  pppi=sqrt(abs(enq1**2-amp3**2))
4300 #else
4301  enq1=(am2sq-amp3**2+amp4**2)/(2*am2)
4302  enq2=(am2sq+amp3**2-amp4**2)/(2*am2)
4303  ppi= enq1**2-amp4**2
4304  pppi=sqrt(abs(enq1**2-amp4**2))
4305 #endif
4306  phspac=phspac*(4*pi)*(2*pppi/am2)
4307 #if defined (ALEPH)
4308 * piz momentum in 2 rest frame(piz is 3rd pi)
4309 #else
4310 * piz momentum in 2 rest frame
4311 #endif
4312  CALL sphera(pppi,piz)
4313  piz(4)=enq1
4314 #if defined (ALEPH)
4315 c pipl momentum in 2 rest frame(pipl is 4th pi)
4316 #else
4317 * pipl momentum in 2 rest frame
4318 #endif
4319  DO 30 i=1,3
4320  30 pipl(i)=-piz(i)
4321  pipl(4)=enq2
4322 * 3 rest frame, define pim1
4323 #if defined (ALEPH)
4324 c pr momentum
4325 #else
4326 * pr momentum
4327 #endif
4328  pr(1)=0
4329  pr(2)=0
4330  pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4331  pr(3)= sqrt(abs(pr(4)**2-am2**2))
4332  ppi = pr(4)**2-am2**2
4333 #if defined (ALEPH)
4334 c pim1 momentum
4335 #else
4336 * pim1 momentum
4337 #endif
4338  pim1(1)=0
4339  pim1(2)=0
4340  pim1(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4341  pim1(3)=-pr(3)
4342 c --- this part of jacobian will be recovered later
4343  ff3=(4*pi)*(2*pr(3)/am3)
4344 * old pions boosted from 2 rest frame to 3 rest frame
4345  exe=(pr(4)+pr(3))/am2
4346  CALL bostr3(exe,piz,piz)
4347  CALL bostr3(exe,pipl,pipl)
4348  rr3=rrr(3)
4349  rr4=rrr(4)
4350  thet =acos(-1.+2*rr3)
4351  phi = 2*pi*rr4
4352  CALL rotpol(thet,phi,pipl)
4353  CALL rotpol(thet,phi,pim1)
4354  CALL rotpol(thet,phi,piz)
4355  CALL rotpol(thet,phi,pr)
4356 #if defined (ALEPH)
4357 c 4 rest frame, define pim2
4358 c pr momentum
4359 #else
4360 * 4 rest frame, define pim2
4361 * pr momentum
4362 #endif
4363  pr(1)=0
4364  pr(2)=0
4365  pr(4)=1./(2*am4)*(am4**2+am3**2-amp1**2)
4366  pr(3)= sqrt(abs(pr(4)**2-am3**2))
4367  ppi = pr(4)**2-am3**2
4368 #if defined (ALEPH)
4369 c pim2 momentum
4370 #else
4371 * pim2 momentum
4372 #endif
4373  pim2(1)=0
4374  pim2(2)=0
4375  pim2(4)=1./(2*am4)*(am4**2-am3**2+amp1**2)
4376  pim2(3)=-pr(3)
4377 c --- this part of jacobian will be recovered later
4378  ff4=(4*pi)*(2*pr(3)/am4)
4379 * old pions boosted from 3 rest frame to 4 rest frame
4380  exe=(pr(4)+pr(3))/am3
4381  CALL bostr3(exe,piz,piz)
4382  CALL bostr3(exe,pipl,pipl)
4383  CALL bostr3(exe,pim1,pim1)
4384  rr3=rrr(7)
4385  rr4=rrr(8)
4386  thet =acos(-1.+2*rr3)
4387  phi = 2*pi*rr4
4388  CALL rotpol(thet,phi,pipl)
4389  CALL rotpol(thet,phi,pim1)
4390  CALL rotpol(thet,phi,pim2)
4391  CALL rotpol(thet,phi,piz)
4392  CALL rotpol(thet,phi,pr)
4393 c
4394 * now to the tau rest frame, define paa and neutrino momenta
4395 * paa momentum
4396  paa(1)=0
4397  paa(2)=0
4398  paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am4**2)
4399  paa(3)= sqrt(abs(paa(4)**2-am4**2))
4400  ppi = paa(4)**2-am4**2
4401  phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4402  phsp=phsp*(4*pi)*(2*paa(3)/amtau)
4403 * tau-neutrino momentum
4404  pn(1)=0
4405  pn(2)=0
4406  pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am4**2)
4407  pn(3)=-paa(3)
4408 c we include remaining part of the jacobian
4409 c --- flat channel
4410  am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4411  $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4412  ams2=(am4-amp2)**2
4413  ams1=(amp1+amp3+amp4)**2
4414  ff1=(ams2-ams1)
4415  ams1=(amp3+amp4)**2
4416  ams2=(sqrt(am3sq)-amp1)**2
4417  ff2=ams2-ams1
4418  ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4419  ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4420  uu=ff1*ff2*ff3*ff4
4421 c --- first channel
4422  am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4423  $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4424  ams2=(am4-amp2)**2
4425  ams1=(amp1+amp3+amp4)**2
4426  alp1=atan((ams1-amrx**2)/amrx/gamrx)
4427  alp2=atan((ams2-amrx**2)/amrx/gamrx)
4428  ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4429  ff1=ff1*(alp2-alp1)
4430  ams1=(amp3+amp4)**2
4431  ams2=(sqrt(am3sq)-amp1)**2
4432  ff2=ams2-ams1
4433  ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4434  ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4435  ff=ff1*ff2*ff3*ff4
4436 c --- second channel
4437  am3sq=(pim2(4)+piz(4)+pipl(4))**2-(pim2(3)+piz(3)+pipl(3))**2
4438  $ -(pim2(2)+piz(2)+pipl(2))**2-(pim2(1)+piz(1)+pipl(1))**2
4439  ams2=(am4-amp1)**2
4440  ams1=(amp2+amp3+amp4)**2
4441  alp1=atan((ams1-amrx**2)/amrx/gamrx)
4442  alp2=atan((ams2-amrx**2)/amrx/gamrx)
4443  gg1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4444  gg1=gg1*(alp2-alp1)
4445  ams1=(amp3+amp4)**2
4446  ams2=(sqrt(am3sq)-amp2)**2
4447  gg2=ams2-ams1
4448  gg3=(4*pi)*(xlam(am2**2,amp2**2,am3sq)/am3sq)
4449  gg4=(4*pi)*(xlam(am3sq,amp1**2,am4**2)/am4**2)
4450  gg=gg1*gg2*gg3*gg4
4451 c --- jacobian averaged over the two
4452  IF ( ( (ff+gg)*uu+ff*gg ).GT.0.0d0) THEN
4453  rr=ff*gg*uu/(0.5*prez*(ff+gg)*uu+(1.0-prez)*ff*gg)
4454  phspac=phspac*rr
4455  ELSE
4456  phspac=0.0
4457  ENDIF
4458 * momenta of the two pi-minus are randomly symmetrised
4459  IF (jnpi.EQ.1) THEN
4460  rr5= rrr(5)
4461  IF(rr5.LE.0.5) THEN
4462  DO 70 i=1,4
4463  x=pim1(i)
4464  pim1(i)=pim2(i)
4465  70 pim2(i)=x
4466  ENDIF
4467  phspac=phspac/2.
4468  ELSE
4469 c momenta of pi0-s are generated uniformly only IF prez=0.0
4470  rr5= rrr(5)
4471  IF(rr5.LE.0.5) THEN
4472  DO 71 i=1,4
4473  x=pim1(i)
4474  pim1(i)=pim2(i)
4475  71 pim2(i)=x
4476  ENDIF
4477  phspac=phspac/6.
4478  ENDIF
4479 * all pions boosted from 4 rest frame to tau rest frame
4480 * z-axis antiparallel to neutrino momentum
4481  exe=(paa(4)+paa(3))/am4
4482  CALL bostr3(exe,piz,piz)
4483  CALL bostr3(exe,pipl,pipl)
4484  CALL bostr3(exe,pim1,pim1)
4485  CALL bostr3(exe,pim2,pim2)
4486  CALL bostr3(exe,pr,pr)
4487 c partial width consists of phase space and amplitude
4488 c check on consistency with dadnpi, THEN, code breakes uniform pion
4489 c distribution in hadronic system
4490 #if defined (ALEPH)
4491  CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4492 #else
4493 cam assume neutrino mass=0. and sum over final polarisation
4494 c amx2=am4**2
4495 c brak= 2*(amtau**2-amx2) * (amtau**2+2.*amx2)
4496 c amplit=ccabib**2*gfermi**2/2. * brak * amx2*sigee(amx2,1)
4497  IF (jnpi.EQ.1) THEN
4498  CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4499  ELSEIF (jnpi.EQ.2) THEN
4500  CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv)
4501  ENDIF
4502 #endif
4503  dgamt=1/(2.*amtau)*amplit*phspac
4504 c phase space check
4505 c dgamt=phspac
4506  DO 77 k=1,4
4507  pmult(k,1)=pim1(k)
4508  pmult(k,2)=pim2(k)
4509 #if defined (ALEPH)
4510  pmult(k,3)=piz(k)
4511  pmult(k,4)=pipl(k)
4512 #else
4513  pmult(k,3)=pipl(k)
4514  pmult(k,4)=piz(k)
4515 #endif
4516  77 CONTINUE
4517  END
4518  SUBROUTINE dam4pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,AMPLIT,HV)
4519 c ----------------------------------------------------------------------
4520 * calculates differential cross section and polarimeter vector
4521 * for tau decay into 4 pi modes
4522 * all spin effects in the full decay chain are taken into account.
4523 * calculations done in tau rest frame with z-axis along neutrino moment
4524 c mnum decay mode identifier.
4525 c
4526 #if defined (ALEPH)
4527 c called by : dph4pi
4528 #else
4529 c called by : dphsaa
4530 #endif
4531 c ----------------------------------------------------------------------
4532  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4533  * ,ampiz,ampi,amro,gamro,ama1,gama1
4534  * ,amk,amkz,amkst,gamkst
4535 c
4536  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4537  * ,ampiz,ampi,amro,gamro,ama1,gama1
4538  * ,amk,amkz,amkst,gamkst
4539  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4540  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
4541  REAL hv(4),pt(4),pn(4),pim1(4),pim2(4),pim3(4),pim4(4)
4542  REAL pivec(4),piaks(4),hvm(4)
4543  COMPLEX hadcur(4),form1,form2,form3,form4,form5
4544  EXTERNAL form1,form2,form3,form4,form5
4545  DATA pi /3.141592653589793238462643/
4546  DATA icont /0/
4547 c
4548 #if defined (CLEO)
4549  CALL curr_cleo(mnum,pim1,pim2,pim3,pim4,hadcur)
4550 #else
4551  CALL curr(mnum,pim1,pim2,pim3,pim4,hadcur)
4552 #endif
4553 c
4554 * calculate pi-vectors: vector and axial
4555  CALL clvec(hadcur,pn,pivec)
4556  CALL claxi(hadcur,pn,piaks)
4557  CALL clnut(hadcur,brakm,hvm)
4558 * spin independent part of decay diff-cross-sect. in tau rest frame
4559  brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4560  & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4561  amplit=(ccabib*gfermi)**2*brak/2.
4562 c polarimeter vector in tau rest frame
4563  DO 90 i=1,3
4564  hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4565  & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4566 c hv is defined for tau- with gamma=b+hv*pol
4567  IF (brak.NE.0.0)
4568  &hv(i)=-hv(i)/brak
4569  90 CONTINUE
4570  END
4571  SUBROUTINE dph5pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4572 c ----------------------------------------------------------------------
4573 * it simulates 5pi decay in tau rest frame with
4574 * z-axis along 5pi momentum
4575 c ----------------------------------------------------------------------
4576  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4577  * ,ampiz,ampi,amro,gamro,ama1,gama1
4578  * ,amk,amkz,amkst,gamkst
4579 c
4580  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4581  * ,ampiz,ampi,amro,gamro,ama1,gama1
4582 
4583 
4584  * ,amk,amkz,amkst,gamkst
4585  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4586  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
4587  parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4588 #if defined (ALEPH)
4589  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4590 #else
4591  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4592 #endif
4593  & ,names
4594  CHARACTER names(nmode)*31
4595  REAL hv(4),pt(4),pn(4),paa(4),pmult(4,9)
4596  REAL*4 pr(4),pi1(4),pi2(4),pi3(4),pi4(4),pi5(4)
4597  REAL*8 amp1,amp2,amp3,amp4,amp5,ams1,ams2,amom,gamom
4598  REAL*8 am5sq,am4sq,am3sq,am2sq,am5,am4,am3
4599  REAL*4 rrr(10)
4600  REAL*8 gg1,gg2,gg3,ff1,ff2,ff3,ff4,alp,alp1,alp2
4601 #if defined (ALEPH)
4602  REAL*8 xm,am,gammab
4603 #else
4604  REAL*8 xm,am,gamma
4605 ccm.s.>>>>>>
4606  real*8 phspac
4607 ccm.s.<<<<<<
4608 #endif
4609  DATA pi /3.141592653589793238462643/
4610  DATA icont /0/
4611  data fpi /93.3e-3/
4612 c
4613  COMPLEX bwign
4614 c
4615 #if defined (ALEPH)
4616  bwign(xm,am,gammab)=xm**2/cmplx(xm**2-am**2,gammab*am)
4617 #else
4618  bwign(xm,am,gamma)=xm**2/cmplx(xm**2-am**2,gamma*am)
4619 #endif
4620 
4621 c
4622  amom=.782
4623  gamom=0.0085
4624 c
4625 c 6 body phase space normalised as in bjorken-drell
4626 c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
4627  phspac=1./2**29/pi**14
4628 c phspac=1./2**5/pi**2
4629 c init 5pi decay mode(jnpi)
4630  amp1=dcdmas(idffin(1,jnpi))
4631  amp2=dcdmas(idffin(2,jnpi))
4632  amp3=dcdmas(idffin(3,jnpi))
4633  amp4=dcdmas(idffin(4,jnpi))
4634  amp5=dcdmas(idffin(5,jnpi))
4635 c
4636 c tau momentum
4637  pt(1)=0.
4638  pt(2)=0.
4639  pt(3)=0.
4640  pt(4)=amtau
4641 c
4642  CALL ranmar(rrr,10)
4643 c
4644 c masses of 5, 4, 3 and 2 pi systems
4645 c 3 pi with sampling for omega resonance
4646 cam
4647 c mass of 5 (12345)
4648  rr1=rrr(10)
4649  ams1=(amp1+amp2+amp3+amp4+amp5)**2
4650  ams2=(amtau-amnuta)**2
4651  am5sq=ams1+ rr1*(ams2-ams1)
4652  am5 =sqrt(am5sq)
4653  phspac=phspac*(ams2-ams1)
4654 c
4655 c mass of 4 (2345)
4656 c flat phase space
4657  rr1=rrr(9)
4658  ams1=(amp2+amp3+amp4+amp5)**2
4659  ams2=(am5-amp1)**2
4660  am4sq=ams1+ rr1*(ams2-ams1)
4661  am4 =sqrt(am4sq)
4662  gg1=ams2-ams1
4663 c
4664 c mass of 3 (234)
4665 c phase space with sampling for omega resonance
4666  rr1=rrr(1)
4667  ams1=(amp2+amp3+amp4)**2
4668  ams2=(am4-amp5)**2
4669  alp1=atan((ams1-amom**2)/amom/gamom)
4670  alp2=atan((ams2-amom**2)/amom/gamom)
4671  alp=alp1+rr1*(alp2-alp1)
4672  am3sq =amom**2+amom*gamom*tan(alp)
4673  am3 =sqrt(am3sq)
4674 c --- this part of the jacobian will be recovered later ---------------
4675  gg2=((am3sq-amom**2)**2+(amom*gamom)**2)/(amom*gamom)
4676  gg2=gg2*(alp2-alp1)
4677 c flat phase space;
4678 c am3sq=ams1+ rr1*(ams2-ams1)
4679 c am3 =sqrt(am3sq)
4680 c --- this part of jacobian will be recovered later
4681 c gg2=ams2-ams1
4682 c
4683 c mass of 2 (34)
4684  rr2=rrr(2)
4685  ams1=(amp3+amp4)**2
4686  ams2=(am3-amp2)**2
4687 c flat phase space;
4688  am2sq=ams1+ rr2*(ams2-ams1)
4689  am2 =sqrt(am2sq)
4690 c --- this part of jacobian will be recovered later
4691  gg3=ams2-ams1
4692 c
4693 c(34) restframe, define pi3 and pi4
4694  enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4695  enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4696  ppi= enq1**2-amp3**2
4697  pppi=sqrt(abs(enq1**2-amp3**2))
4698  ff1=(4*pi)*(2*pppi/am2)
4699 c pi3 momentum in(34) rest frame
4700  call sphera(pppi,pi3)
4701  pi3(4)=enq1
4702 c pi4 momentum in(34) rest frame
4703  do 30 i=1,3
4704  30 pi4(i)=-pi3(i)
4705  pi4(4)=enq2
4706 c
4707 c(234) rest frame, define pi2
4708 c pr momentum
4709  pr(1)=0
4710  pr(2)=0
4711  pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4712  pr(3)= sqrt(abs(pr(4)**2-am2**2))
4713  ppi = pr(4)**2-am2**2
4714 c pi2 momentum
4715  pi2(1)=0
4716  pi2(2)=0
4717  pi2(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4718  pi2(3)=-pr(3)
4719 c --- this part of jacobian will be recovered later
4720  ff2=(4*pi)*(2*pr(3)/am3)
4721 c old pions boosted from 2 rest frame to 3 rest frame
4722  exe=(pr(4)+pr(3))/am2
4723  call bostr3(exe,pi3,pi3)
4724  call bostr3(exe,pi4,pi4)
4725  rr3=rrr(3)
4726  rr4=rrr(4)
4727  thet =acos(-1.+2*rr3)
4728  phi = 2*pi*rr4
4729  call rotpol(thet,phi,pi2)
4730  call rotpol(thet,phi,pi3)
4731  call rotpol(thet,phi,pi4)
4732 c
4733 c(2345) rest frame, define pi5
4734 c pr momentum
4735  pr(1)=0
4736  pr(2)=0
4737  pr(4)=1./(2*am4)*(am4**2+am3**2-amp5**2)
4738  pr(3)= sqrt(abs(pr(4)**2-am3**2))
4739  ppi = pr(4)**2-am3**2
4740 c pi5 momentum
4741  pi5(1)=0
4742  pi5(2)=0
4743  pi5(4)=1./(2*am4)*(am4**2-am3**2+amp5**2)
4744  pi5(3)=-pr(3)
4745 c --- this part of jacobian will be recovered later
4746  ff3=(4*pi)*(2*pr(3)/am4)
4747 c old pions boosted from 3 rest frame to 4 rest frame
4748  exe=(pr(4)+pr(3))/am3
4749  call bostr3(exe,pi2,pi2)
4750  call bostr3(exe,pi3,pi3)
4751  call bostr3(exe,pi4,pi4)
4752  rr3=rrr(5)
4753  rr4=rrr(6)
4754  thet =acos(-1.+2*rr3)
4755  phi = 2*pi*rr4
4756  call rotpol(thet,phi,pi2)
4757  call rotpol(thet,phi,pi3)
4758  call rotpol(thet,phi,pi4)
4759  call rotpol(thet,phi,pi5)
4760 c
4761 c(12345) rest frame, define pi1
4762 c pr momentum
4763  pr(1)=0
4764  pr(2)=0
4765  pr(4)=1./(2*am5)*(am5**2+am4**2-amp1**2)
4766  pr(3)= sqrt(abs(pr(4)**2-am4**2))
4767  ppi = pr(4)**2-am4**2
4768 c pi1 momentum
4769  pi1(1)=0
4770  pi1(2)=0
4771  pi1(4)=1./(2*am5)*(am5**2-am4**2+amp1**2)
4772  pi1(3)=-pr(3)
4773 c --- this part of jacobian will be recovered later
4774  ff4=(4*pi)*(2*pr(3)/am5)
4775 c old pions boosted from 4 rest frame to 5 rest frame
4776  exe=(pr(4)+pr(3))/am4
4777  call bostr3(exe,pi2,pi2)
4778  call bostr3(exe,pi3,pi3)
4779  call bostr3(exe,pi4,pi4)
4780  call bostr3(exe,pi5,pi5)
4781  rr3=rrr(7)
4782  rr4=rrr(8)
4783  thet =acos(-1.+2*rr3)
4784  phi = 2*pi*rr4
4785  call rotpol(thet,phi,pi1)
4786  call rotpol(thet,phi,pi2)
4787  call rotpol(thet,phi,pi3)
4788  call rotpol(thet,phi,pi4)
4789  call rotpol(thet,phi,pi5)
4790 c
4791 * now to the tau rest frame, define paa and neutrino momenta
4792 * paa momentum
4793  paa(1)=0
4794  paa(2)=0
4795 c paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5**2)
4796 c paa(3)= sqrt(abs(paa(4)**2-am5**2))
4797 c ppi = paa(4)**2-am5**2
4798  paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5sq)
4799  paa(3)= sqrt(abs(paa(4)**2-am5sq))
4800  ppi = paa(4)**2-am5sq
4801  phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4802 * tau-neutrino momentum
4803  pn(1)=0
4804  pn(2)=0
4805  pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am5**2)
4806  pn(3)=-paa(3)
4807 c
4808  phspac=phspac * gg1*gg2*gg3*ff1*ff2*ff3*ff4
4809 c
4810 c all pions boosted from 5 rest frame to tau rest frame
4811 c z-axis antiparallel to neutrino momentum
4812  exe=(paa(4)+paa(3))/am5
4813  call bostr3(exe,pi1,pi1)
4814  call bostr3(exe,pi2,pi2)
4815  call bostr3(exe,pi3,pi3)
4816  call bostr3(exe,pi4,pi4)
4817  call bostr3(exe,pi5,pi5)
4818 c
4819 c partial width consists of phase space and amplitude
4820 c amplitude(cf ys.tsai phys.rev.d4,2821(1971)
4821 c or f.gilman sh.rhie phys.rev.d31,1066(1985)
4822 c
4823  pxq=amtau*paa(4)
4824  pxn=amtau*pn(4)
4825  qxn=paa(4)*pn(4)-paa(1)*pn(1)-paa(2)*pn(2)-paa(3)*pn(3)
4826  brak=2*(gv**2+ga**2)*(2*pxq*qxn+am5sq*pxn)
4827  & -6*(gv**2-ga**2)*amtau*amnuta*am5sq
4828  fompp = cabs(bwign(am3,amom,gamom))**2
4829 c normalisation factor(to some numerical undimensioned factor;
4830 c cf r.fischer et al zphys c3, 313 (1980))
4831  fnorm = 1/fpi**6
4832 c amplit=ccabib**2*gfermi**2/2. * brak * am5sq*sigee(am5sq,jnpi)
4833  amplit=ccabib**2*gfermi**2/2. * brak
4834  amplit = amplit * fompp * fnorm
4835 c phase space test
4836 c amplit = amplit * fnorm
4837  dgamt=1/(2.*amtau)*amplit*phspac
4838 c ignore spin terms
4839  DO 40 i=1,3
4840  40 hv(i)=0.
4841 c
4842  do 77 k=1,4
4843  pmult(k,1)=pi1(k)
4844  pmult(k,2)=pi2(k)
4845  pmult(k,3)=pi3(k)
4846  pmult(k,4)=pi4(k)
4847  pmult(k,5)=pi5(k)
4848  77 continue
4849  return
4850 #if defined (ALEPH)
4851 c missing: transposition of identical particles, statistical factors
4852 c for identical matrices, polarimetric vector. matrix element rather nai
4853 #else
4854 c missing: transposition of identical particles, startistical factors
4855 c for identical matrices, polarimetric vector. matrix element rather naive.
4856 #endif
4857 c flat phase space in pion system + with breit wigner for omega
4858 c anyway it is better than nothing, and code is improvable.
4859  end
4860  SUBROUTINE dphsrk(DGAMT,HV,PN,PR,PMULT,INUM)
4861 c ----------------------------------------------------------------------
4862 c it simulates rho decay in tau rest frame with
4863 c z-axis along rho momentum
4864 c rho decays to k kbar
4865 c ----------------------------------------------------------------------
4866  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4867  * ,ampiz,ampi,amro,gamro,ama1,gama1
4868  * ,amk,amkz,amkst,gamkst
4869 c
4870  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4871  * ,ampiz,ampi,amro,gamro,ama1,gama1
4872  * ,amk,amkz,amkst,gamkst
4873  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4874  REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
4875  REAL hv(4),pt(4),pn(4),pr(4),pkc(4),pkz(4),qq(4),pmult(4,9)
4876 #if defined (ALEPH)
4877  REAL*4 rr1(1)
4878 #else
4879  REAL rr1(1)
4880 #endif
4881  DATA pi /3.141592653589793238462643/
4882  DATA icont /0/
4883 c
4884 c three body phase space normalised as in bjorken-drell
4885  phspac=1./2**11/pi**5
4886 c tau momentum
4887  pt(1)=0.
4888  pt(2)=0.
4889  pt(3)=0.
4890  pt(4)=amtau
4891 c mass of(real/virtual) rho
4892  ams1=(amk+amkz)**2
4893  ams2=(amtau-amnuta)**2
4894 c flat phase space
4895  CALL ranmar(rr1,1)
4896  amx2=ams1+ rr1(1)*(ams2-ams1)
4897  amx=sqrt(amx2)
4898  phspac=phspac*(ams2-ams1)
4899 c phase space with sampling for rho resonance
4900 c alp1=atan((ams1-amro**2)/amro/gamro)
4901 c alp2=atan((ams2-amro**2)/amro/gamro)
4902 cam
4903  100 CONTINUE
4904 c CALL ranmar(rr1,1)
4905 c alp=alp1+rr1(1)*(alp2-alp1)
4906 c amx2=amro**2+amro*gamro*tan(alp)
4907 c amx=sqrt(amx2)
4908 c IF(amx.LT.(amk+amkz)) go to 100
4909 cam
4910 c phspac=phspac*((amx2-amro**2)**2+(amro*gamro)**2)/(amro*gamro)
4911 c phspac=phspac*(alp2-alp1)
4912 c
4913 c tau-neutrino momentum
4914  pn(1)=0
4915  pn(2)=0
4916  pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
4917  pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
4918 c rho momentum
4919  pr(1)=0
4920  pr(2)=0
4921  pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
4922  pr(3)=-pn(3)
4923  phspac=phspac*(4*pi)*(2*pr(3)/amtau)
4924 c
4925 cam
4926  enq1=(amx2+amk**2-amkz**2)/(2.*amx)
4927  enq2=(amx2-amk**2+amkz**2)/(2.*amx)
4928  pppi=sqrt((enq1-amk)*(enq1+amk))
4929  phspac=phspac*(4*pi)*(2*pppi/amx)
4930 c charged pi momentum in rho rest frame
4931  CALL sphera(pppi,pkc)
4932  pkc(4)=enq1
4933 c neutral pi momentum in rho rest frame
4934  DO 20 i=1,3
4935 20 pkz(i)=-pkc(i)
4936  pkz(4)=enq2
4937  exe=(pr(4)+pr(3))/amx
4938 c pions boosted from rho rest frame to tau rest frame
4939  CALL bostr3(exe,pkc,pkc)
4940  CALL bostr3(exe,pkz,pkz)
4941  DO 30 i=1,4
4942  30 qq(i)=pkc(i)-pkz(i)
4943 c qq transverse to pr
4944  pksd =pr(4)*pr(4)-pr(3)*pr(3)-pr(2)*pr(2)-pr(1)*pr(1)
4945  qqpks=pr(4)* qq(4)-pr(3)* qq(3)-pr(2)* qq(2)-pr(1)* qq(1)
4946  DO 31 i=1,4
4947 31 qq(i)=qq(i)-pr(i)*qqpks/pksd
4948 c amplitude
4949  prodpq=pt(4)*qq(4)
4950  prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
4951  prodpn=pt(4)*pn(4)
4952  qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
4953  brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
4954  & +(gv**2-ga**2)*amtau*amnuta*qq2
4955  amplit=(gfermi*ccabib)**2*brak*2*fpirk(amx)
4956  dgamt=1/(2.*amtau)*amplit*phspac
4957  DO 40 i=1,3
4958  40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
4959  do 77 k=1,4
4960  pmult(k,1)=pkc(k)
4961  pmult(k,2)=pkz(k)
4962  77 continue
4963  RETURN
4964  END
4965  FUNCTION fpirk(W)
4966 c ----------------------------------------------------------
4967 c square of pion form factor
4968 c ----------------------------------------------------------
4969  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4970  * ,ampiz,ampi,amro,gamro,ama1,gama1
4971  * ,amk,amkz,amkst,gamkst
4972 c
4973  REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4974  * ,ampiz,ampi,amro,gamro,ama1,gama1
4975  * ,amk,amkz,amkst,gamkst
4976 c COMPLEX fpikmk
4977  COMPLEX fpikm
4978  fpirk=cabs(fpikm(w,amk,amkz))**2
4979 c fpirk=cabs(fpikmk(w,amk,amkz))**2
4980  END
4981  COMPLEX FUNCTION fpikmk(W,XM1,XM2)
4982 c **********************************************************
4983 c kaon form factor
4984 c **********************************************************
4985  COMPLEX bwigm
4986  REAL rom,rog,rom1,rog1,beta1,pi,pim,s,w
4987  EXTERNAL bwig
4988  DATA init /0/
4989 c
4990 c ------------ parameters --------------------
4991  IF (init.EQ.0 ) THEN
4992  init=1
4993  pi=3.141592654
4994  pim=.140
4995  rom=0.773
4996  rog=0.145
4997  rom1=1.570
4998  rog1=0.510
4999 c beta1=-0.111
5000  beta1=-0.221
5001  ENDIF
5002 c -----------------------------------------------
5003  s=w**2
5004  fpikmk=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
5005  & /(1+beta1)
5006  RETURN
5007  END
5008  SUBROUTINE reslux
5009 c ****************
5010 c initialize lund COMMON
5011 #if defined (CLEO)
5012 #else
5013  parameter(nmxhep=2000)
5014  common/hepevtx/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
5015  &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
5016  SAVE /hepevtx/
5017 #endif
5018  nhep=0
5019  END
5020  SUBROUTINE dwrph(KTO,PHX)
5021 c
5022 c -------------------------
5023 c
5024  IMPLICIT REAL*8 (a-h,o-z)
5025  REAL*4 phx(4)
5026  REAL*4 qhot(4)
5027 c
5028  DO 9 k=1,4
5029  qhot(k) =0.0
5030  9 CONTINUE
5031 c CASE of tau radiative decays.
5032 c filling of the lund COMMON block.
5033  DO 1002 i=1,4
5034  1002 qhot(i)=phx(i)
5035  IF (qhot(4).GT.1.e-5) CALL dwluph(kto,qhot)
5036  RETURN
5037  END
5038  SUBROUTINE dwluph(KTO,PHOT)
5039 c---------------------------------------------------------------------
5040 c lorentz transformation to cmsystem and
5041 c updating of hepevt record
5042 c
5043 c called by : dexay1,(dekay1,dekay2)
5044 c
5045 c used when radiative corrections in decays are generated
5046 c---------------------------------------------------------------------
5047 c
5048 #if defined (ALEPH)
5049  COMMON /taupos/ np1,np2
5050 #else
5051 #endif
5052  REAL phot(4)
5053 #if defined (ALEPH)
5054 #else
5055  COMMON /taupos/ np1,np2
5056 #endif
5057 c
5058 c check energy
5059  IF (phot(4).LE.0.0) RETURN
5060 c
5061 c position of decaying particle:
5062  IF((kto.EQ. 1).OR.(kto.EQ.11)) THEN
5063  nps=np1
5064  ELSE
5065  nps=np2
5066  ENDIF
5067 c
5068  ktos=kto
5069  IF(ktos.GT.10) ktos=ktos-10
5070 c boost and append photon(gamma is 22)
5071  CALL tralo4(ktos,phot,phot,am)
5072  CALL filhep(0,1,22,nps,nps,0,0,phot,0.0,.true.)
5073 c
5074  RETURN
5075  END
5076 
5077  SUBROUTINE dwluel(KTO,ISGN,PNU,PWB,PEL,PNE)
5078 c ----------------------------------------------------------------------
5079 c lorentz transformation to cmsystem and
5080 c updating of hepevt record
5081 c
5082 c isgn = 1/-1 for tau-/tau+
5083 c
5084 c called by : dexay,(dekay1,dekay2)
5085 c ----------------------------------------------------------------------
5086 c
5087 #if defined (ALEPH)
5088  COMMON /taupos/ np1,np2
5089 #else
5090 #endif
5091  REAL pnu(4),pwb(4),pel(4),pne(4)
5092 #if defined (ALEPH)
5093 #else
5094  COMMON /taupos/ np1,np2
5095 #endif
5096 c
5097 c position of decaying particle:
5098  IF(kto.EQ. 1) THEN
5099  nps=np1
5100  ELSE
5101  nps=np2
5102  ENDIF
5103 c
5104 c tau neutrino(nu_tau is 16)
5105  CALL tralo4(kto,pnu,pnu,am)
5106  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5107 c
5108 c w boson(w+ is 24)
5109  CALL tralo4(kto,pwb,pwb,am)
5110 c CALL filhep(0,2,-24*isgn,nps,nps,0,0,pwb,am,.true.)
5111 c
5112 c electron(e- is 11)
5113  CALL tralo4(kto,pel,pel,am)
5114  CALL filhep(0,1,11*isgn,nps,nps,0,0,pel,am,.false.)
5115 c
5116 c anti electron neutrino(nu_e is 12)
5117  CALL tralo4(kto,pne,pne,am)
5118  CALL filhep(0,1,-12*isgn,nps,nps,0,0,pne,am,.true.)
5119 c
5120  RETURN
5121  END
5122  SUBROUTINE dwlumu(KTO,ISGN,PNU,PWB,PMU,PNM)
5123 c ----------------------------------------------------------------------
5124 c lorentz transformation to cmsystem and
5125 c updating of hepevt record
5126 c
5127 c isgn = 1/-1 for tau-/tau+
5128 c
5129 c called by : dexay,(dekay1,dekay2)
5130 c ----------------------------------------------------------------------
5131 c
5132 #if defined (ALEPH)
5133  COMMON /taupos/ np1,np2
5134 #else
5135 #endif
5136  REAL pnu(4),pwb(4),pmu(4),pnm(4)
5137 #if defined (ALEPH)
5138 #else
5139  COMMON /taupos/ np1,np2
5140 #endif
5141 c
5142 c position of decaying particle:
5143  IF(kto.EQ. 1) THEN
5144  nps=np1
5145  ELSE
5146  nps=np2
5147  ENDIF
5148 c
5149 c tau neutrino(nu_tau is 16)
5150  CALL tralo4(kto,pnu,pnu,am)
5151  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5152 c
5153 c w boson(w+ is 24)
5154  CALL tralo4(kto,pwb,pwb,am)
5155 c CALL filhep(0,2,-24*isgn,nps,nps,0,0,pwb,am,.true.)
5156 c
5157 c muon(mu- is 13)
5158  CALL tralo4(kto,pmu,pmu,am)
5159  CALL filhep(0,1,13*isgn,nps,nps,0,0,pmu,am,.false.)
5160 c
5161 c anti muon neutrino(nu_mu is 14)
5162  CALL tralo4(kto,pnm,pnm,am)
5163  CALL filhep(0,1,-14*isgn,nps,nps,0,0,pnm,am,.true.)
5164 c
5165  RETURN
5166  END
5167  SUBROUTINE dwlupi(KTO,ISGN,PPI,PNU)
5168 c ----------------------------------------------------------------------
5169 c lorentz transformation to cmsystem and
5170 c updating of hepevt record
5171 c
5172 c isgn = 1/-1 for tau-/tau+
5173 c
5174 c called by : dexay,(dekay1,dekay2)
5175 c ----------------------------------------------------------------------
5176 c
5177  REAL pnu(4),ppi(4)
5178  COMMON /taupos/ np1,np2
5179 c
5180 c position of decaying particle:
5181  IF(kto.EQ. 1) THEN
5182  nps=np1
5183  ELSE
5184  nps=np2
5185  ENDIF
5186 c
5187 c tau neutrino(nu_tau is 16)
5188  CALL tralo4(kto,pnu,pnu,am)
5189  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5190 c
5191 c charged pi meson(pi+ is 211)
5192  CALL tralo4(kto,ppi,ppi,am)
5193  CALL filhep(0,1,-211*isgn,nps,nps,0,0,ppi,am,.true.)
5194 c
5195  RETURN
5196  END
5197  SUBROUTINE dwluro(KTO,ISGN,PNU,PRHO,PIC,PIZ)
5198 c ----------------------------------------------------------------------
5199 c lorentz transformation to cmsystem and
5200 c updating of hepevt record
5201 c
5202 c isgn = 1/-1 for tau-/tau+
5203 c
5204 c called by : dexay,(dekay1,dekay2)
5205 c ----------------------------------------------------------------------
5206 c
5207 #if defined (ALEPH)
5208  COMMON /taupos/ np1,np2
5209 #else
5210 #endif
5211  REAL pnu(4),prho(4),pic(4),piz(4)
5212 #if defined (ALEPH)
5213 #else
5214  COMMON /taupos/ np1,np2
5215 #endif
5216 c
5217 c position of decaying particle:
5218  IF(kto.EQ. 1) THEN
5219  nps=np1
5220  ELSE
5221  nps=np2
5222  ENDIF
5223 c
5224 c tau neutrino(nu_tau is 16)
5225  CALL tralo4(kto,pnu,pnu,am)
5226  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5227 c
5228 c charged rho meson(rho+ is 213)
5229  CALL tralo4(kto,prho,prho,am)
5230  CALL filhep(0,2,-213*isgn,nps,nps,0,0,prho,am,.true.)
5231 c
5232 c charged pi meson(pi+ is 211)
5233  CALL tralo4(kto,pic,pic,am)
5234  CALL filhep(0,1,-211*isgn,-1,-1,0,0,pic,am,.true.)
5235 c
5236 c pi0 meson(pi0 is 111)
5237  CALL tralo4(kto,piz,piz,am)
5238  CALL filhep(0,1,111,-2,-2,0,0,piz,am,.true.)
5239 c
5240  RETURN
5241  END
5242  SUBROUTINE dwluaa(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
5243 c ----------------------------------------------------------------------
5244 c lorentz transformation to cmsystem and
5245 c updating of hepevt record
5246 c
5247 c isgn = 1/-1 for tau-/tau+
5248 c jaa = 1 (2) for a_1- decay to pi+ 2pi- (pi- 2pi0)
5249 c
5250 c called by : dexay,(dekay1,dekay2)
5251 c ----------------------------------------------------------------------
5252 c
5253 #if defined (ALEPH)
5254  COMMON /taupos/ np1,np2
5255 #else
5256 #endif
5257  REAL pnu(4),paa(4),pim1(4),pim2(4),pipl(4)
5258 #if defined (ALEPH)
5259 #else
5260  COMMON /taupos/ np1,np2
5261 #endif
5262 c
5263 c position of decaying particle:
5264  IF(kto.EQ. 1) THEN
5265  nps=np1
5266  ELSE
5267  nps=np2
5268  ENDIF
5269 c
5270 c tau neutrino(nu_tau is 16)
5271  CALL tralo4(kto,pnu,pnu,am)
5272  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5273 c
5274 c charged a_1 meson(a_1+ is 20213)
5275  CALL tralo4(kto,paa,paa,am)
5276  CALL filhep(0,1,-20213*isgn,nps,nps,0,0,paa,am,.true.)
5277 c
5278 c two possible decays of the charged a1 meson
5279  IF(jaa.EQ.1) THEN
5280 c
5281 c a1 --> pi+ pi- pi- (or charged conjugate)
5282 c
5283 c pi minus(or c.c.) (pi+ is 211)
5284  CALL tralo4(kto,pim2,pim2,am)
5285  CALL filhep(0,1,-211*isgn,-1,-1,0,0,pim2,am,.true.)
5286 c
5287 c pi minus(or c.c.) (pi+ is 211)
5288  CALL tralo4(kto,pim1,pim1,am)
5289  CALL filhep(0,1,-211*isgn,-2,-2,0,0,pim1,am,.true.)
5290 c
5291 c pi plus(or c.c.) (pi+ is 211)
5292  CALL tralo4(kto,pipl,pipl,am)
5293  CALL filhep(0,1, 211*isgn,-3,-3,0,0,pipl,am,.true.)
5294 c
5295  ELSE IF (jaa.EQ.2) THEN
5296 c
5297 c a1 --> pi- pi0 pi0(or charged conjugate)
5298 c
5299 c pi zero(pi0 is 111)
5300  CALL tralo4(kto,pim2,pim2,am)
5301  CALL filhep(0,1,111,-1,-1,0,0,pim2,am,.true.)
5302 c
5303 c pi zero(pi0 is 111)
5304  CALL tralo4(kto,pim1,pim1,am)
5305  CALL filhep(0,1,111,-2,-2,0,0,pim1,am,.true.)
5306 c
5307 c pi minus(or c.c.) (pi+ is 211)
5308  CALL tralo4(kto,pipl,pipl,am)
5309  CALL filhep(0,1,-211*isgn,-3,-3,0,0,pipl,am,.true.)
5310 c
5311  ENDIF
5312 c
5313  RETURN
5314  END
5315  SUBROUTINE dwlukk (KTO,ISGN,PKK,PNU)
5316 c ----------------------------------------------------------------------
5317 c lorentz transformation to cmsystem and
5318 c updating of hepevt record
5319 c
5320 c isgn = 1/-1 for tau-/tau+
5321 c
5322 c ----------------------------------------------------------------------
5323 c
5324  REAL pkk(4),pnu(4)
5325  COMMON /taupos/ np1,np2
5326 c
5327 c position of decaying particle
5328 #if defined (ALEPH)
5329  IF(kto.EQ. 1) THEN
5330 #else
5331  IF (kto.EQ.1) THEN
5332 #endif
5333  nps=np1
5334  ELSE
5335  nps=np2
5336  ENDIF
5337 c
5338 c tau neutrino(nu_tau is 16)
5339  CALL tralo4(kto,pnu,pnu,am)
5340  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5341 c
5342 c k meson(k+ is 321)
5343  CALL tralo4(kto,pkk,pkk,am)
5344  CALL filhep(0,1,-321*isgn,nps,nps,0,0,pkk,am,.true.)
5345 c
5346  RETURN
5347  END
5348  SUBROUTINE dwluks(KTO,ISGN,PNU,PKS,PKK,PPI,JKST)
5349  COMMON / taukle / bra1,brk0,brk0b,brks
5350  REAL*4 bra1,brk0,brk0b,brks
5351 #if defined (ALEPH)
5352  COMMON /taupos/ np1,np2
5353  REAL*4 xio(1)
5354 #endif
5355 c ----------------------------------------------------------------------
5356 c lorentz transformation to cmsystem and
5357 c updating of hepevt record
5358 c
5359 c isgn = 1/-1 for tau-/tau+
5360 c jkst=10 (20) corresponds to k0b pi- (k- pi0) decay
5361 c
5362 c ----------------------------------------------------------------------
5363 c
5364 #if defined (ALEPH)
5365  REAL pnu(4),pks(4),pkk(4),ppi(4)
5366 #else
5367  REAL pnu(4),pks(4),pkk(4),ppi(4),xio(1)
5368  COMMON /taupos/ np1,np2
5369 #endif
5370 c
5371 c position of decaying particle
5372  IF(kto.EQ. 1) THEN
5373  nps=np1
5374  ELSE
5375  nps=np2
5376  ENDIF
5377 c
5378 c tau neutrino(nu_tau is 16)
5379  CALL tralo4(kto,pnu,pnu,am)
5380  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5381 c
5382 c charged k* meson(k*+ is 323)
5383  CALL tralo4(kto,pks,pks,am)
5384  CALL filhep(0,1,-323*isgn,nps,nps,0,0,pks,am,.true.)
5385 c
5386 c two possible decay modes of charged k*
5387  IF(jkst.EQ.10) THEN
5388 c
5389 c k*- --> pi- k0b(or charged conjugate)
5390 c
5391 c charged pi meson(pi+ is 211)
5392  CALL tralo4(kto,ppi,ppi,am)
5393  CALL filhep(0,1,-211*isgn,-1,-1,0,0,ppi,am,.true.)
5394 c
5395  bran=brk0b
5396  IF (isgn.EQ.-1) bran=brk0
5397 c k0 --> k0_long(is 130) / k0_short(is 310) = 1/1
5398  CALL ranmar(xio,1)
5399  IF(xio(1).GT.bran) THEN
5400  k0type = 130
5401  ELSE
5402  k0type = 310
5403  ENDIF
5404 c
5405  CALL tralo4(kto,pkk,pkk,am)
5406  CALL filhep(0,1,k0type,-2,-2,0,0,pkk,am,.true.)
5407 c
5408  ELSE IF(jkst.EQ.20) THEN
5409 c
5410 c k*- --> pi0 k-
5411 c
5412 c pi zero(pi0 is 111)
5413  CALL tralo4(kto,ppi,ppi,am)
5414  CALL filhep(0,1,111,-1,-1,0,0,ppi,am,.true.)
5415 c
5416 c charged k meson(k+ is 321)
5417  CALL tralo4(kto,pkk,pkk,am)
5418  CALL filhep(0,1,-321*isgn,-2,-2,0,0,pkk,am,.true.)
5419 c
5420  ENDIF
5421 c
5422  RETURN
5423  END
5424  SUBROUTINE dwlnew(KTO,ISGN,PNU,PWB,PNPI,MODE)
5425 c ----------------------------------------------------------------------
5426 c lorentz transformation to cmsystem and
5427 c updating of hepevt record
5428 c
5429 c isgn = 1/-1 for tau-/tau+
5430 c
5431 c called by : dexay,(dekay1,dekay2)
5432 c ----------------------------------------------------------------------
5433 c
5434  parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
5435 #if defined (ALEPH)
5436  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5437 #else
5438  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5439 #endif
5440  & ,names
5441  COMMON /taupos/ np1,np2
5442  CHARACTER names(nmode)*31
5443  REAL pnu(4),pwb(4),pnpi(4,9)
5444  REAL ppi(4)
5445 c
5446  jnpi=mode-7
5447 c position of decaying particle
5448  IF(kto.EQ. 1) THEN
5449  nps=np1
5450  ELSE
5451  nps=np2
5452  ENDIF
5453 c
5454 c tau neutrino(nu_tau is 16)
5455  CALL tralo4(kto,pnu,pnu,am)
5456  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5457 c
5458 c w boson(w+ is 24)
5459  CALL tralo4(kto,pwb,pwb,am)
5460  CALL filhep(0,1,-24*isgn,nps,nps,0,0,pwb,am,.true.)
5461 c
5462 c multi pi mode jnpi
5463 c
5464 c get multiplicity of mode jnpi
5465  nd=mulpik(jnpi)
5466  DO i=1,nd
5467 #if defined (ALEPH)
5468 cam kfpi=lunpik(idffin(i,jnpi),-isgn)
5469  kfpi=lunpik(idffin(i,jnpi), isgn)
5470 #else
5471  kfpi=lunpik(idffin(i,jnpi),-isgn)
5472 #endif
5473 c for charged conjugate case, change charged pions only
5474 c IF(kfpi.NE.111)kfpi=kfpi*isgn
5475  DO j=1,4
5476  ppi(j)=pnpi(j,i)
5477  END DO
5478  CALL tralo4(kto,ppi,ppi,am)
5479  CALL filhep(0,1,kfpi,-i,-i,0,0,ppi,am,.true.)
5480  END DO
5481 c
5482  RETURN
5483  END
5484 #if defined (CePeCe)
5485 #else
5486 #endif
5487  FUNCTION amast(PP)
5488 c ----------------------------------------------------------------------
5489 c calculates mass of pp(double precision)
5490 c
5491 c used by : radkor
5492 c ----------------------------------------------------------------------
5493  IMPLICIT REAL*8 (a-h,o-z)
5494  REAL*8 pp(4)
5495  aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5496 c
5497  IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5498  amast=aaa
5499  RETURN
5500  END
5501  FUNCTION amas4(PP)
5502 c ******************
5503 c ----------------------------------------------------------------------
5504 c calculates mass of pp
5505 c
5506 c used by :
5507 c ----------------------------------------------------------------------
5508  REAL pp(4)
5509  aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5510  IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5511  amas4=aaa
5512  RETURN
5513  END
5514  FUNCTION angxy(X,Y)
5515 c ----------------------------------------------------------------------
5516 c
5517 c used by : koralz radkor
5518 c ----------------------------------------------------------------------
5519  IMPLICIT DOUBLE PRECISION (a-h,o-z)
5520  DATA pi /3.141592653589793238462643d0/
5521 c
5522  IF(abs(y).LT.abs(x)) THEN
5523  the=atan(abs(y/x))
5524  IF(x.LE.0d0) the=pi-the
5525  ELSE
5526  the=acos(x/sqrt(x**2+y**2))
5527  ENDIF
5528  angxy=the
5529  RETURN
5530  END
5531  FUNCTION angfi(X,Y)
5532 c ----------------------------------------------------------------------
5533 * calculates angle in(0,2*pi) range out of x-y
5534 c
5535 c used by : koralz radkor
5536 c ----------------------------------------------------------------------
5537  IMPLICIT DOUBLE PRECISION (a-h,o-z)
5538  DATA pi /3.141592653589793238462643d0/
5539 c
5540  IF(abs(y).LT.abs(x)) THEN
5541  the=atan(abs(y/x))
5542  IF(x.LE.0d0) the=pi-the
5543  ELSE
5544  the=acos(x/sqrt(x**2+y**2))
5545  ENDIF
5546  IF(y.LT.0d0) the=2d0*pi-the
5547  angfi=the
5548  END
5549  SUBROUTINE rotod1(PH1,PVEC,QVEC)
5550 c ----------------------------------------------------------------------
5551 c
5552 c used by : koralz
5553 c ----------------------------------------------------------------------
5554  IMPLICIT DOUBLE PRECISION (a-h,o-z)
5555  dimension pvec(4),qvec(4),rvec(4)
5556 c
5557  phi=ph1
5558  cs=cos(phi)
5559  sn=sin(phi)
5560  DO 10 i=1,4
5561  10 rvec(i)=pvec(i)
5562  qvec(1)=rvec(1)
5563  qvec(2)= cs*rvec(2)-sn*rvec(3)
5564  qvec(3)= sn*rvec(2)+cs*rvec(3)
5565  qvec(4)=rvec(4)
5566  RETURN
5567  END
5568  SUBROUTINE rotod2(PH1,PVEC,QVEC)
5569 c ----------------------------------------------------------------------
5570 c
5571 c used by : koralz radkor
5572 c ----------------------------------------------------------------------
5573  IMPLICIT DOUBLE PRECISION (a-h,o-z)
5574  dimension pvec(4),qvec(4),rvec(4)
5575 c
5576  phi=ph1
5577  cs=cos(phi)
5578  sn=sin(phi)
5579  DO 10 i=1,4
5580  10 rvec(i)=pvec(i)
5581  qvec(1)= cs*rvec(1)+sn*rvec(3)
5582  qvec(2)=rvec(2)
5583  qvec(3)=-sn*rvec(1)+cs*rvec(3)
5584  qvec(4)=rvec(4)
5585  RETURN
5586  END
5587  SUBROUTINE rotod3(PH1,PVEC,QVEC)
5588 c ----------------------------------------------------------------------
5589 c
5590 c used by : koralz radkor
5591 c ----------------------------------------------------------------------
5592  IMPLICIT DOUBLE PRECISION (a-h,o-z)
5593 c
5594  dimension pvec(4),qvec(4),rvec(4)
5595  phi=ph1
5596  cs=cos(phi)
5597  sn=sin(phi)
5598  DO 10 i=1,4
5599  10 rvec(i)=pvec(i)
5600  qvec(1)= cs*rvec(1)-sn*rvec(2)
5601  qvec(2)= sn*rvec(1)+cs*rvec(2)
5602  qvec(3)=rvec(3)
5603  qvec(4)=rvec(4)
5604  END
5605  SUBROUTINE bostr3(EXE,PVEC,QVEC)
5606 c ----------------------------------------------------------------------
5607 c boost along z axis, exe=exp(eta), eta= hiperbolic velocity.
5608 c
5609 c used by : tauola koralz(?)
5610 c ----------------------------------------------------------------------
5611  REAL*4 pvec(4),qvec(4),rvec(4)
5612 c
5613  DO 10 i=1,4
5614  10 rvec(i)=pvec(i)
5615  rpl=rvec(4)+rvec(3)
5616  rmi=rvec(4)-rvec(3)
5617  qpl=rpl*exe
5618  qmi=rmi/exe
5619  qvec(1)=rvec(1)
5620  qvec(2)=rvec(2)
5621  qvec(3)=(qpl-qmi)/2
5622  qvec(4)=(qpl+qmi)/2
5623  END
5624  SUBROUTINE bostd3(EXE,PVEC,QVEC)
5625 c ----------------------------------------------------------------------
5626 c boost along z axis, exe=exp(eta), eta= hiperbolic velocity.
5627 c
5628 c used by : koralz radkor
5629 c ----------------------------------------------------------------------
5630  IMPLICIT DOUBLE PRECISION (a-h,o-z)
5631  dimension pvec(4),qvec(4),rvec(4)
5632 c
5633  DO 10 i=1,4
5634  10 rvec(i)=pvec(i)
5635  rpl=rvec(4)+rvec(3)
5636  rmi=rvec(4)-rvec(3)
5637  qpl=rpl*exe
5638  qmi=rmi/exe
5639  qvec(1)=rvec(1)
5640  qvec(2)=rvec(2)
5641  qvec(3)=(qpl-qmi)/2
5642  qvec(4)=(qpl+qmi)/2
5643  RETURN
5644  END
5645  SUBROUTINE rotor1(PH1,PVEC,QVEC)
5646 c ----------------------------------------------------------------------
5647 c
5648 c called by :
5649 c ----------------------------------------------------------------------
5650  REAL*4 pvec(4),qvec(4),rvec(4)
5651 c
5652  phi=ph1
5653  cs=cos(phi)
5654  sn=sin(phi)
5655  DO 10 i=1,4
5656  10 rvec(i)=pvec(i)
5657  qvec(1)=rvec(1)
5658  qvec(2)= cs*rvec(2)-sn*rvec(3)
5659  qvec(3)= sn*rvec(2)+cs*rvec(3)
5660  qvec(4)=rvec(4)
5661  END
5662  SUBROUTINE rotor2(PH1,PVEC,QVEC)
5663 c ----------------------------------------------------------------------
5664 c
5665 c used by : tauola
5666 c ----------------------------------------------------------------------
5667  IMPLICIT REAL*4(a-h,o-z)
5668  REAL*4 pvec(4),qvec(4),rvec(4)
5669 c
5670  phi=ph1
5671  cs=cos(phi)
5672  sn=sin(phi)
5673  DO 10 i=1,4
5674  10 rvec(i)=pvec(i)
5675  qvec(1)= cs*rvec(1)+sn*rvec(3)
5676  qvec(2)=rvec(2)
5677  qvec(3)=-sn*rvec(1)+cs*rvec(3)
5678  qvec(4)=rvec(4)
5679  END
5680  SUBROUTINE rotor3(PHI,PVEC,QVEC)
5681 c ----------------------------------------------------------------------
5682 c
5683 c used by : tauola
5684 c ----------------------------------------------------------------------
5685  REAL*4 pvec(4),qvec(4),rvec(4)
5686 c
5687  cs=cos(phi)
5688  sn=sin(phi)
5689  DO 10 i=1,4
5690  10 rvec(i)=pvec(i)
5691  qvec(1)= cs*rvec(1)-sn*rvec(2)
5692  qvec(2)= sn*rvec(1)+cs*rvec(2)
5693  qvec(3)=rvec(3)
5694  qvec(4)=rvec(4)
5695  END
5696  SUBROUTINE spherd(R,X)
5697 c ----------------------------------------------------------------------
5698 c generates uniformly three-vector x on sphere of radius r
5699 c double precison version of sphera
5700 c ----------------------------------------------------------------------
5701  REAL*8 r,x(4),pi,costh,sinth
5702  REAL*4 rrr(2)
5703  DATA pi /3.141592653589793238462643d0/
5704 c
5705  CALL ranmar(rrr,2)
5706  costh=-1+2*rrr(1)
5707  sinth=sqrt(1 -costh**2)
5708  x(1)=r*sinth*cos(2*pi*rrr(2))
5709  x(2)=r*sinth*sin(2*pi*rrr(2))
5710  x(3)=r*costh
5711  RETURN
5712  END
5713  SUBROUTINE rotpox(THET,PHI,PP)
5714  IMPLICIT REAL*8 (a-h,o-z)
5715 c ----------------------------------------------------------------------
5716 #if defined (ALEPH)
5717 c double precison version of rotpol
5718 #else
5719 c
5720 #endif
5721 c ----------------------------------------------------------------------
5722  dimension pp(4)
5723 c
5724  CALL rotod2(thet,pp,pp)
5725  CALL rotod3( phi,pp,pp)
5726  RETURN
5727  END
5728  SUBROUTINE sphera(R,X)
5729 c ----------------------------------------------------------------------
5730 c generates uniformly three-vector x on sphere of radius r
5731 c
5732 c called by : dphsxx,dadmpi,dadmkk
5733 c ----------------------------------------------------------------------
5734  REAL x(4)
5735  REAL*4 rrr(2)
5736  DATA pi /3.141592653589793238462643/
5737 c
5738  CALL ranmar(rrr,2)
5739  costh=-1.+2.*rrr(1)
5740  sinth=sqrt(1.-costh**2)
5741  x(1)=r*sinth*cos(2*pi*rrr(2))
5742  x(2)=r*sinth*sin(2*pi*rrr(2))
5743  x(3)=r*costh
5744  RETURN
5745  END
5746  SUBROUTINE rotpol(THET,PHI,PP)
5747 c ----------------------------------------------------------------------
5748 c
5749 c called by : dadmaa,dphsaa
5750 c ----------------------------------------------------------------------
5751  REAL pp(4)
5752 c
5753  CALL rotor2(thet,pp,pp)
5754  CALL rotor3( phi,pp,pp)
5755  RETURN
5756  END
5757 #include "../randg/tauola-random.h"
5758  FUNCTION dilogt(X)
5759 c *****************
5760  IMPLICIT REAL*8(a-h,o-z)
5761 cern c304 version 29/07/71 dilog 59 c
5762  z=-1.64493406684822
5763  IF(x .LT.-1.0) go to 1
5764  IF(x .LE. 0.5) go to 2
5765  IF(x .EQ. 1.0) go to 3
5766  IF(x .LE. 2.0) go to 4
5767  z=3.2898681336964
5768  1 t=1.0/x
5769  s=-0.5
5770  z=z-0.5* log(abs(x))**2
5771  go to 5
5772  2 t=x
5773  s=0.5
5774  z=0.
5775  go to 5
5776  3 dilogt=1.64493406684822
5777  RETURN
5778  4 t=1.0-x
5779  s=-0.5
5780  z=1.64493406684822 - log(x)* log(abs(t))
5781  5 y=2.66666666666666 *t+0.66666666666666
5782  b= 0.00000 00000 00001
5783  a=y*b +0.00000 00000 00004
5784  b=y*a-b+0.00000 00000 00011
5785  a=y*b-a+0.00000 00000 00037
5786  b=y*a-b+0.00000 00000 00121
5787  a=y*b-a+0.00000 00000 00398
5788  b=y*a-b+0.00000 00000 01312
5789  a=y*b-a+0.00000 00000 04342
5790  b=y*a-b+0.00000 00000 14437
5791  a=y*b-a+0.00000 00000 48274
5792  b=y*a-b+0.00000 00001 62421
5793  a=y*b-a+0.00000 00005 50291
5794  b=y*a-b+0.00000 00018 79117
5795  a=y*b-a+0.00000 00064 74338
5796  b=y*a-b+0.00000 00225 36705
5797  a=y*b-a+0.00000 00793 87055
5798  b=y*a-b+0.00000 02835 75385
5799  a=y*b-a+0.00000 10299 04264
5800  b=y*a-b+0.00000 38163 29463
5801  a=y*b-a+0.00001 44963 00557
5802  b=y*a-b+0.00005 68178 22718
5803  a=y*b-a+0.00023 20021 96094
5804  b=y*a-b+0.00100 16274 96164
5805  a=y*b-a+0.00468 63619 59447
5806  b=y*a-b+0.02487 93229 24228
5807  a=y*b-a+0.16607 30329 27855
5808  a=y*a-b+1.93506 43008 6996
5809  dilogt=s*t*(a-b)+z
5810  RETURN
5811 c=======================================================================
5812 c===================END of cpc part ====================================
5813 c=======================================================================
5814  END
5815