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