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