C++ Interface to Tauola
tauola/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))
3246  & - hj(2)*aimag(hj(2)) - hj(1)*aimag(hj(1)) )
3247  RETURN
3248  END
3249  SUBROUTINE dampog(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3250 c ----------------------------------------------------------------------
3251 * calculates differential cross section and polarimeter vector
3252 * for tau decay into a1, a1 decays next into rho+pi and rho into pi+pi.
3253 * all spin effects in the full decay chain are taken into account.
3254 * calculations done in tau rest frame with z-axis along neutrino moment
3255 * the routine is writen for zero neutrino mass.
3256 c
3257 c called by : dphsaa
3258 c ----------------------------------------------------------------------
3259  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3260  * ,ampiz,ampi,amro,gamro,ama1,gama1
3261  * ,amk,amkz,amkst,gamkst
3262 c
3263  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3264  * ,ampiz,ampi,amro,gamro,ama1,gama1
3265  * ,amk,amkz,amkst,gamkst
3266  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3267  real*4 gfermi,gv,ga,ccabib,scabib,gamel
3268  COMMON /testa1/ keya1
3269  REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4)
3270  REAL PAA(4),VEC1(4),VEC2(4)
3271  REAL PIVEC(4),PIAKS(4),HVM(4)
3272  COMPLEX BWIGN,HADCUR(4),FNORM,FORMOM
3273  DATA icont /1/
3274 * this inline funct. calculates the scalar part of the propagator
3275 c ajwmod to satisfy compiler, comment out this unused function.
3276 c
3277 * four momentum of a1
3278  DO 10 i=1,4
3279  vec1(i)=0.0
3280  vec2(i)=0.0
3281  hv(i) =0.0
3282  10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3283  vec1(1)=1.0
3284 * masses of a1, and of two pi-pairs which may form rho
3285  xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3286  xmom =sqrt(abs( (pim2(4)+pipl(4))**2-(pim2(3)+pipl(3))**2
3287  $ -(pim2(2)+pipl(2))**2-(pim2(1)+pipl(1))**2 ))
3288  xmro2 =(pipl(1))**2 +(pipl(2))**2 +(pipl(3))**2
3289 * elements of hadron current
3290  prod1 =vec1(1)*pipl(1)
3291  prod2 =vec2(2)*pipl(2)
3292  p12 =pim1(4)*pim2(4)-pim1(1)*pim2(1)
3293  $ -pim1(2)*pim2(2)-pim1(3)*pim2(3)
3294  p1pl =pim1(4)*pipl(4)-pim1(1)*pipl(1)
3295  $ -pim1(2)*pipl(2)-pim1(3)*pipl(3)
3296  p2pl =pipl(4)*pim2(4)-pipl(1)*pim2(1)
3297  $ -pipl(2)*pim2(2)-pipl(3)*pim2(3)
3298  DO 40 i=1,3
3299  vec1(i)= (vec1(i)-prod1/xmro2*pipl(i))
3300  40 CONTINUE
3301  gnorm=sqrt(vec1(1)**2+vec1(2)**2+vec1(3)**2)
3302  DO 41 i=1,3
3303  vec1(i)= vec1(i)/gnorm
3304  41 CONTINUE
3305  vec2(1)=(vec1(2)*pipl(3)-vec1(3)*pipl(2))/sqrt(xmro2)
3306  vec2(2)=(vec1(3)*pipl(1)-vec1(1)*pipl(3))/sqrt(xmro2)
3307  vec2(3)=(vec1(1)*pipl(2)-vec1(2)*pipl(1))/sqrt(xmro2)
3308  p1vec1 =pim1(4)*vec1(4)-pim1(1)*vec1(1)
3309  $ -pim1(2)*vec1(2)-pim1(3)*vec1(3)
3310  p2vec1 =vec1(4)*pim2(4)-vec1(1)*pim2(1)
3311  $ -vec1(2)*pim2(2)-vec1(3)*pim2(3)
3312  p1vec2 =pim1(4)*vec2(4)-pim1(1)*vec2(1)
3313  $ -pim1(2)*vec2(2)-pim1(3)*vec2(3)
3314  p2vec2 =vec2(4)*pim2(4)-vec2(1)*pim2(1)
3315  $ -vec2(2)*pim2(2)-vec2(3)*pim2(3)
3316 * hadron current
3317  fnorm=formom(xmaa,xmom)
3318  brak=0.0
3319  DO 120 jj=1,2
3320  DO 45 i=1,4
3321  IF (jj.EQ.1) THEN
3322  hadcur(i) = fnorm *(
3323  $ vec1(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3324  $ -pim2(i)*(p2vec1*p1pl-p1vec1*p2pl)
3325  $ +pipl(i)*(p2vec1*p12 -p1vec1*(ampi**2+p2pl)) )
3326  ELSE
3327  hadcur(i) = fnorm *(
3328  $ vec2(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3329  $ -pim2(i)*(p2vec2*p1pl-p1vec2*p2pl)
3330  $ +pipl(i)*(p2vec2*p12 -p1vec2*(ampi**2+p2pl)) )
3331  ENDIF
3332  45 CONTINUE
3333 c
3334 * calculate pi-vectors: vector and axial
3335  CALL clvec(hadcur,pn,pivec)
3336  CALL claxi(hadcur,pn,piaks)
3337  CALL clnut(hadcur,brakm,hvm)
3338 * spin independent part of decay diff-cross-sect. in tau rest frame
3339  brak=brak+(gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3340  & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3341  DO 90 i=1,3
3342  hv(i)=hv(i)-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3343  & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3344  90 CONTINUE
3345 c hv is defined for tau- with gamma=b+hv*pol
3346  120 CONTINUE
3347  amplit=(gfermi*ccabib)**2*brak/2.
3348 c the statistical factor for identical pi-s was cancelled with
3349 c two, for two modes of a1 decay namelly pi+pi-pi- and pi-pi0pi0
3350 c polarimeter vector in tau rest frame
3351  DO 91 i=1,3
3352  hv(i)=-hv(i)/brak
3353  91 CONTINUE
3354 
3355  END
3356  SUBROUTINE damppk(MNUM,PT,PN,PIM1,PIM2,PIM3,AMPLIT,HV)
3357 c ----------------------------------------------------------------------
3358 * calculates differential cross section and polarimeter vector
3359 * for tau decay into k k pi, k pi pi.
3360 * all spin effects in the full decay chain are taken into account.
3361 * calculations done in tau rest frame with z-axis along neutrino moment
3362 c mnum decay mode identifier.
3363 c
3364 c called by : dphsaa
3365 c ----------------------------------------------------------------------
3366  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3367  * ,ampiz,ampi,amro,gamro,ama1,gama1
3368  * ,amk,amkz,amkst,gamkst
3369 c
3370  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3371  * ,ampiz,ampi,amro,gamro,ama1,gama1
3372  * ,amk,amkz,amkst,gamkst
3373  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3374  real*4 gfermi,gv,ga,ccabib,scabib,gamel
3375  REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4)
3376  REAL PAA(4),VEC1(4),VEC2(4),VEC3(4),VEC4(4),VEC5(4)
3377  REAL PIVEC(4),PIAKS(4),HVM(4)
3378  REAL FNORM(0:7),COEF(1:5,0:7)
3379  COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5,UROJ
3380  COMPLEX F1,F2,F3,F4,F5
3381  EXTERNAL form1,form2,form3,form4,form5
3382  DATA pi /3.141592653589793238462643/
3383  DATA icont /0/
3384 c
3385  DATA fpi /93.3e-3/
3386  IF (icont.EQ.0) THEN
3387  icont=1
3388  uroj=cmplx(0.0,1.0)
3389  dwapi0=sqrt(2.0)
3390  fnorm(0)=ccabib/fpi
3391  fnorm(1)=ccabib/fpi
3392  fnorm(2)=ccabib/fpi
3393  fnorm(3)=ccabib/fpi
3394  fnorm(4)=scabib/fpi/dwapi0
3395  fnorm(5)=scabib/fpi
3396  fnorm(6)=scabib/fpi
3397  fnorm(7)=ccabib/fpi
3398 c
3399  coef(1,0)= 2.0*sqrt(2.)/3.0
3400  coef(2,0)=-2.0*sqrt(2.)/3.0
3401 c ajw 2/98: add in the d-wave and i=0 3pi substructure:
3402  coef(3,0)= 2.0*sqrt(2.)/3.0
3403  coef(4,0)= fpi
3404  coef(5,0)= 0.0
3405 c
3406  coef(1,1)=-sqrt(2.)/3.0
3407  coef(2,1)= sqrt(2.)/3.0
3408  coef(3,1)= 0.0
3409  coef(4,1)= fpi
3410  coef(5,1)= sqrt(2.)
3411 c
3412  coef(1,2)=-sqrt(2.)/3.0
3413  coef(2,2)= sqrt(2.)/3.0
3414  coef(3,2)= 0.0
3415  coef(4,2)= 0.0
3416  coef(5,2)=-sqrt(2.)
3417 c
3418 c ajw 11/97: add in the k*-prim-s, ala finkemeier&mirkes
3419  coef(1,3)= 1./3.
3420  coef(2,3)=-2./3.
3421  coef(3,3)= 2./3.
3422  coef(4,3)= 0.0
3423  coef(5,3)= 0.0
3424 c
3425  coef(1,4)= 1.0/sqrt(2.)/3.0
3426  coef(2,4)=-1.0/sqrt(2.)/3.0
3427  coef(3,4)= 0.0
3428  coef(4,4)= 0.0
3429  coef(5,4)= 0.0
3430 c
3431  coef(1,5)=-sqrt(2.)/3.0
3432  coef(2,5)= sqrt(2.)/3.0
3433  coef(3,5)= 0.0
3434  coef(4,5)= 0.0
3435  coef(5,5)=-sqrt(2.)
3436 c
3437 c ajw 11/97: add in the k*-prim-s, ala finkemeier&mirkes
3438  coef(1,6)= 1./3.
3439  coef(2,6)=-2./3.
3440  coef(3,6)= 2./3.
3441  coef(4,6)= 0.0
3442  coef(5,6)=-2.0
3443 c
3444  coef(1,7)= 0.0
3445  coef(2,7)= 0.0
3446  coef(3,7)= 0.0
3447  coef(4,7)= 0.0
3448  coef(5,7)=-sqrt(2.0/3.0)
3449 c
3450  ENDIF
3451 c
3452  DO 10 i=1,4
3453  10 paa(i)=pim1(i)+pim2(i)+pim3(i)
3454  xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3455  xmro1 =sqrt(abs((pim3(4)+pim2(4))**2-(pim3(1)+pim2(1))**2
3456  $ -(pim3(2)+pim2(2))**2-(pim3(3)+pim2(3))**2))
3457  xmro2 =sqrt(abs((pim3(4)+pim1(4))**2-(pim3(1)+pim1(1))**2
3458  $ -(pim3(2)+pim1(2))**2-(pim3(3)+pim1(3))**2))
3459  xmro3 =sqrt(abs((pim1(4)+pim2(4))**2-(pim1(1)+pim2(1))**2
3460  $ -(pim1(2)+pim2(2))**2-(pim1(3)+pim2(3))**2))
3461 * elements of hadron current
3462  prod1 =paa(4)*(pim2(4)-pim3(4))-paa(1)*(pim2(1)-pim3(1))
3463  $ -paa(2)*(pim2(2)-pim3(2))-paa(3)*(pim2(3)-pim3(3))
3464  prod2 =paa(4)*(pim3(4)-pim1(4))-paa(1)*(pim3(1)-pim1(1))
3465  $ -paa(2)*(pim3(2)-pim1(2))-paa(3)*(pim3(3)-pim1(3))
3466  prod3 =paa(4)*(pim1(4)-pim2(4))-paa(1)*(pim1(1)-pim2(1))
3467  $ -paa(2)*(pim1(2)-pim2(2))-paa(3)*(pim1(3)-pim2(3))
3468  DO 40 i=1,4
3469  vec1(i)= pim2(i)-pim3(i) -paa(i)*prod1/xmaa**2
3470  vec2(i)= pim3(i)-pim1(i) -paa(i)*prod2/xmaa**2
3471  vec3(i)= pim1(i)-pim2(i) -paa(i)*prod3/xmaa**2
3472  40 vec4(i)= pim1(i)+pim2(i)+pim3(i)
3473  CALL prod5(pim1,pim2,pim3,vec5)
3474 * hadron current
3475 c be aware that sign of vec2 is opposite to sign of vec1 in a1 case
3476 c rationalize this code:
3477  f1 = cmplx(coef(1,mnum))*form1(mnum,xmaa**2,xmro1**2,xmro2**2)
3478  f2 = cmplx(coef(2,mnum))*form2(mnum,xmaa**2,xmro2**2,xmro1**2)
3479  f3 = cmplx(coef(3,mnum))*form3(mnum,xmaa**2,xmro3**2,xmro1**2)
3480  f4 = (-1.0*uroj)*
3481  $cmplx(coef(4,mnum))*form4(mnum,xmaa**2,xmro1**2,xmro2**2,xmro3**2)
3482  f5 = (-1.0)*uroj/4.0/pi**2/fpi**2*
3483  $ cmplx(coef(5,mnum))*form5(mnum,xmaa**2,xmro1**2,xmro2**2)
3484 
3485  DO 45 i=1,4
3486  hadcur(i)= cmplx(fnorm(mnum)) * (
3487  $ cmplx(vec1(i))*f1+cmplx(vec2(i))*f2+cmplx(vec3(i))*f3+
3488  $ cmplx(vec4(i))*f4+cmplx(vec5(i))*f5)
3489  45 CONTINUE
3490 c
3491 * calculate pi-vectors: vector and axial
3492  CALL clvec(hadcur,pn,pivec)
3493  CALL claxi(hadcur,pn,piaks)
3494  CALL clnut(hadcur,brakm,hvm)
3495 * spin independent part of decay diff-cross-sect. in tau rest frame
3496  brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3497  & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3498  amplit=(gfermi)**2*brak/2.
3499  IF (mnum.GE.9) THEN
3500  print *, 'MNUM=',mnum
3501  znak=-1.0
3502  xm1=0.0
3503  xm2=0.0
3504  xm3=0.0
3505  DO 77 k=1,4
3506  IF (k.EQ.4) znak=1.0
3507  xm1=znak*pim1(k)**2+xm1
3508  xm2=znak*pim2(k)**2+xm2
3509  xm3=znak*pim3(k)**2+xm3
3510  77 print *, 'PIM1=',pim1(k),'PIM2=',pim2(k),'PIM3=',pim3(k)
3511  print *, 'XM1=',sqrt(xm1),'XM2=',sqrt(xm2),'XM3=',sqrt(xm3)
3512  print *, '************************************************'
3513  ENDIF
3514 c polarimeter vector in tau rest frame
3515  DO 90 i=1,3
3516  hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3517  & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3518 c hv is defined for tau- with gamma=b+hv*pol
3519  hv(i)=-hv(i)/brak
3520  90 CONTINUE
3521  END
3522  SUBROUTINE prod5(P1,P2,P3,PIA)
3523 c ----------------------------------------------------------------------
3524 c external product of P1, P2, P3 4-momenta.
3525 c sign is chosen +/- for decay of tau +/- respectively
3526 c called by : dampaa, clnut
3527 c ----------------------------------------------------------------------
3528  COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3529  COMMON / idfc / idff
3530  REAL PIA(4),P1(4),P2(4),P3(4)
3531  det2(i,j)=p1(i)*p2(j)-p2(i)*p1(j)
3532 * -----------------------------------
3533  IF (ktom.EQ.1.OR.ktom.EQ.-1) THEN
3534  sign= idff/abs(idff)
3535  ELSEIF (ktom.EQ.2) THEN
3536  sign=-idff/abs(idff)
3537  ELSE
3538  print *, 'STOP IN PROD5: KTOM=',ktom
3539  stop
3540  ENDIF
3541 c
3542 c epsilon( p1(1), p2(2), p3(3), (4) ) = 1
3543 c
3544  pia(1)= -p3(3)*det2(2,4)+p3(4)*det2(2,3)+p3(2)*det2(3,4)
3545  pia(2)= -p3(4)*det2(1,3)+p3(3)*det2(1,4)-p3(1)*det2(3,4)
3546  pia(3)= p3(4)*det2(1,2)-p3(2)*det2(1,4)+p3(1)*det2(2,4)
3547  pia(4)= p3(3)*det2(1,2)-p3(2)*det2(1,3)+p3(1)*det2(2,3)
3548 c all four indices are up so pia(3) and pia(4) have same sign
3549  DO 20 i=1,4
3550  20 pia(i)=pia(i)*sign
3551  END
3552 
3553  SUBROUTINE dexnew(MODE,ISGN,POL,PNU,PAA,PNPI,JNPI)
3554 c ----------------------------------------------------------------------
3555 * this simulates tau decay in tau rest frame
3556 * into nu a1, next a1 decays into rho pi and finally rho into pi pi.
3557 * output four momenta: pnu tauneutrino,
3558 * paa a1
3559 * pim1 pion minus(or pi0) 1 (for tau minus)
3560 * pim2 pion minus(or pi0) 2
3561 * pipl pion plus(or pi-)
3562 * (pipl,pim1) form a rho
3563 c ----------------------------------------------------------------------
3564  COMMON / inout / inut,iout
3565  REAL POL(4),HV(4),PAA(4),PNU(4),PNPI(4,9),RN(1)
3566  DATA iwarm/0/
3567 c
3568  IF(mode.EQ.-1) THEN
3569 c ===================
3570  iwarm=1
3571  CALL dadnew( -1,isgn,hv,pnu,paa,pnpi,jdumm)
3572 cc CALL hbook1(816,'WEIGHT DISTRIBUTION DEXAA $',100,-2.,2.)
3573 c
3574  ELSEIF(mode.EQ. 0) THEN
3575 * =======================
3576  300 CONTINUE
3577  IF(iwarm.EQ.0) GOTO 902
3578  CALL dadnew( 0,isgn,hv,pnu,paa,pnpi,jnpi)
3579  wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
3580 cc CALL hfill(816,wt)
3581  CALL ranmar(rn,1)
3582  IF(rn(1).GT.wt) GOTO 300
3583 c
3584  ELSEIF(mode.EQ. 1) THEN
3585 * =======================
3586  CALL dadnew( 1,isgn,hv,pnu,paa,pnpi,jdumm)
3587 cc CALL hprint(816)
3588  ENDIF
3589 c =====
3590  RETURN
3591  902 WRITE(iout, 9020)
3592  9020 FORMAT(' ----- DEXNEW: LACK OF INITIALISATION')
3593  stop
3594  END
3595  SUBROUTINE dadnew(MODE,ISGN,HV,PNU,PWB,PNPI,JNPI)
3596 c ----------------------------------------------------------------------
3597  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3598  * ,ampiz,ampi,amro,gamro,ama1,gama1
3599  * ,amk,amkz,amkst,gamkst
3600 c
3601  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3602  * ,ampiz,ampi,amro,gamro,ama1,gama1
3603  * ,amk,amkz,amkst,gamkst
3604  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3605  real*4 gfermi,gv,ga,ccabib,scabib,gamel
3606  COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
3607  real*4 gampmc ,gamper
3608  COMMON / inout / inut,iout
3609  parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
3610  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3611  & ,names
3612  CHARACTER NAMES(NMODE)*31
3613 
3614  real*4 pnu(4),pwb(4),pnpi(4,9),hv(4),hhv(4)
3615  real*4 pdum1(4),pdum2(4),pdumi(4,9)
3616  real*4 rrr(3)
3617  real*4 wtmax(nmode)
3618  real*8 swt(nmode),sswt(nmode)
3619  dimension nevraw(nmode),nevovr(nmode),nevacc(nmode)
3620 c
3621  DATA pi /3.141592653589793238462643/
3622  DATA iwarm/0/
3623 c
3624  IF(mode.EQ.-1) THEN
3625 c ===================
3626 c -- at the moment only two decay modes of multipions have m. elem
3627  nmod=nmode
3628  iwarm=1
3629 c print 7003
3630  DO 1 jnpi=1,nmod
3631  nevraw(jnpi)=0
3632  nevacc(jnpi)=0
3633  nevovr(jnpi)=0
3634  swt(jnpi)=0
3635  sswt(jnpi)=0
3636  wtmax(jnpi)=-1.
3637 c for 4pi phase space, need lots more trials at initialization,
3638 c or use the wtmax determined with many trials for default model:
3639  ntrials = 500
3640  IF (jnpi.LE.nm4) THEN
3641  a = pkorb(3,37+jnpi)
3642  IF (a.LT.0.) THEN
3643  ntrials = 10000
3644  ELSE
3645  ntrials = 0
3646  wtmax(jnpi)=a
3647  END IF
3648  END IF
3649  DO i=1,ntrials
3650  IF (jnpi.LE.0) THEN
3651  GOTO 903
3652  ELSEIF(jnpi.LE.nm4) THEN
3653  CALL dph4pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3654  ELSEIF(jnpi.LE.nm4+nm5) THEN
3655  CALL dph5pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3656  ELSEIF(jnpi.LE.nm4+nm5+nm6) THEN
3657  CALL dphnpi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3658  ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3) THEN
3659  inum=jnpi-nm4-nm5-nm6
3660  CALL dphspk(wt,hv,pdum1,pdum2,pdumi,inum)
3661  ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2) THEN
3662  inum=jnpi-nm4-nm5-nm6-nm3
3663  CALL dphsrk(wt,hv,pdum1,pdum2,pdumi,inum)
3664  ELSE
3665  GOTO 903
3666  ENDIF
3667  IF(wt.GT.wtmax(jnpi)/1.2) wtmax(jnpi)=wt*1.2
3668  ENDDO
3669 c print *,' DADNEW JNPI,NTRIALS,WTMAX =',jnpi,ntrials,wtmax(jnpi)
3670 c CALL hbook1(801,'WEIGHT DISTRIBUTION DADNPI $',100,0.,2.,.0)
3671 c print 7004,wtmax(jnpi)
3672 1 CONTINUE
3673  WRITE(iout,7005)
3674 c
3675  ELSEIF(mode.EQ. 0) THEN
3676 c =======================
3677  IF(iwarm.EQ.0) GOTO 902
3678 c
3679 300 CONTINUE
3680  IF (jnpi.LE.0) THEN
3681  GOTO 903
3682  ELSEIF(jnpi.LE.nm4) THEN
3683  CALL dph4pi(wt,hhv,pnu,pwb,pnpi,jnpi)
3684  ELSEIF(jnpi.LE.nm4+nm5) THEN
3685  CALL dph5pi(wt,hhv,pnu,pwb,pnpi,jnpi)
3686  ELSEIF(jnpi.LE.nm4+nm5+nm6) THEN
3687  CALL dphnpi(wt,hhv,pnu,pwb,pnpi,jnpi)
3688  ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3) THEN
3689  inum=jnpi-nm4-nm5-nm6
3690  CALL dphspk(wt,hhv,pnu,pwb,pnpi,inum)
3691  ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2) THEN
3692  inum=jnpi-nm4-nm5-nm6-nm3
3693  CALL dphsrk(wt,hhv,pnu,pwb,pnpi,inum)
3694  ELSE
3695  GOTO 903
3696  ENDIF
3697  DO i=1,4
3698  hv(i)=-isgn*hhv(i)
3699  ENDDO
3700 c CALL hfill(801,wt/wtmax(jnpi))
3701  nevraw(jnpi)=nevraw(jnpi)+1
3702  swt(jnpi)=swt(jnpi)+wt
3703 cccm.s.>>>>>>
3704 cc sswt(jnpi)=sswt(jnpi)+wt**2
3705  sswt(jnpi)=sswt(jnpi)+dble(wt)**2
3706 cccm.s.<<<<<<
3707  CALL ranmar(rrr,3)
3708  rn=rrr(1)
3709  IF(wt.GT.wtmax(jnpi)) nevovr(jnpi)=nevovr(jnpi)+1
3710  IF(rn*wtmax(jnpi).GT.wt) GOTO 300
3711 c rotations to basic tau rest frame
3712  costhe=-1.+2.*rrr(2)
3713  thet=acos(costhe)
3714  phi =2*pi*rrr(3)
3715  CALL rotor2(thet,pnu,pnu)
3716  CALL rotor3( phi,pnu,pnu)
3717  CALL rotor2(thet,pwb,pwb)
3718  CALL rotor3( phi,pwb,pwb)
3719  CALL rotor2(thet,hv,hv)
3720  CALL rotor3( phi,hv,hv)
3721  nd=mulpik(jnpi)
3722  DO 301 i=1,nd
3723  CALL rotor2(thet,pnpi(1,i),pnpi(1,i))
3724  CALL rotor3( phi,pnpi(1,i),pnpi(1,i))
3725 301 CONTINUE
3726  nevacc(jnpi)=nevacc(jnpi)+1
3727 c
3728  ELSEIF(mode.EQ. 1) THEN
3729 c =======================
3730  DO 500 jnpi=1,nmod
3731  IF(nevraw(jnpi).EQ.0) GOTO 500
3732  pargam=swt(jnpi)/float(nevraw(jnpi)+1)
3733  error=0
3734  IF(nevraw(jnpi).NE.0)
3735  & error=sqrt(sswt(jnpi)/swt(jnpi)**2-1./float(nevraw(jnpi)))
3736  rat=pargam/gamel
3737  WRITE(iout, 7010) names(jnpi),
3738  & nevraw(jnpi),nevacc(jnpi),nevovr(jnpi),pargam,rat,error
3739 cc CALL hprint(801)
3740  gampmc(8+jnpi-1)=rat
3741  gamper(8+jnpi-1)=error
3742 cam nevdec(8+jnpi-1)=nevacc(jnpi)
3743  500 CONTINUE
3744  ENDIF
3745 c =====
3746  RETURN
3747  7003 FORMAT(///1x,15(5h*****)
3748  $ /,' *', 25x,'******** DADNEW INITIALISATION ********',9x,1h*
3749  $ )
3750  7004 FORMAT(' *',e20.5,5x,'WTMAX = MAXIMUM WEIGHT ',9x,1h*/)
3751  7005 FORMAT(
3752  $ /,1x,15(5h*****)/)
3753  7010 FORMAT(///1x,15(5h*****)
3754  $ /,' *', 25x,'******** DADNEW FINAL REPORT ******** ',9x,1h*
3755  $ /,' *', 25x,'CHANNEL:',a31 ,9x,1h*
3756  $ /,' *',i20 ,5x,'NEVRAW = NO. OF DECAYS TOTAL ',9x,1h*
3757  $ /,' *',i20 ,5x,'NEVACC = NO. OF DECAYS ACCEPTED ',9x,1h*
3758  $ /,' *',i20 ,5x,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
3759  $ /,' *',e20.5,5x,'PARTIAL WTDTH IN GEV UNITS ',9x,1h*
3760  $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
3761  $ /,' *',f20.8,5x,'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
3762  $ /,1x,15(5h*****)/)
3763  902 WRITE(iout, 9020)
3764  9020 FORMAT(' ----- DADNEW: LACK OF INITIALISATION')
3765  stop
3766  903 WRITE(iout, 9030) jnpi,mode
3767  9030 FORMAT(' ----- DADNEW: WRONG JNPI',2i5)
3768  stop
3769  END
3770 
3771 
3772  SUBROUTINE dph4pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
3773 c ----------------------------------------------------------------------
3774 * it simulates a1 decay in tau rest frame with
3775 * z-axis along a1 momentum
3776 c ----------------------------------------------------------------------
3777  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3778  * ,ampiz,ampi,amro,gamro,ama1,gama1
3779  * ,amk,amkz,amkst,gamkst
3780 c
3781  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3782  * ,ampiz,ampi,amro,gamro,ama1,gama1
3783  * ,amk,amkz,amkst,gamkst
3784  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3785  real*4 gfermi,gv,ga,ccabib,scabib,gamel
3786  REAL HV(4),PT(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4),PMULT(4,9)
3787  REAL PR(4),PIZ(4)
3788  real*4 rrr(9)
3789  real*8 uu,ff,ff1,ff2,ff3,ff4,gg1,gg2,gg3,gg4,rr
3790  DATA pi /3.141592653589793238462643/
3791  DATA icont /0/
3792  xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
3793 c amro, gamro is only a PARAMETER for geting hight efficiency
3794 c
3795 c three body phase space normalised as in bjorken-drell
3796 c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
3797  phspac=1./2**23/pi**11
3798  phsp=1./2**5/pi**2
3799  IF (jnpi.EQ.1) THEN
3800  prez=0.7
3801  amp1=ampi
3802  amp2=ampi
3803  amp3=ampi
3804  amp4=ampiz
3805  amrx=pkorb(1,14)
3806  gamrx=pkorb(2,14)
3807 c ajw: cant simply change amrop, etc, here!
3808 c choice is a by-hand tuning/optimization, no simple relationship
3809 c to actual resonance masses(accd to z.was).
3810 c what matters in the end is what you put in formf/curr .
3811  amrop =1.2
3812  gamrop=.46
3813  ELSE
3814  prez=0.0
3815  amp1=ampiz
3816  amp2=ampiz
3817  amp3=ampiz
3818  amp4=ampi
3819  amrx=1.4
3820  gamrx=.6
3821  amrop =amrx
3822  gamrop=gamrx
3823 
3824  ENDIF
3825  rrb=0.3
3826  CALL choice(100+jnpi,rrb,ichan,prob1,prob2,prob3,
3827  $ amrop,gamrop,amrx,gamrx,amrb,gamrb)
3828  prez=prob1+prob2
3829 c tau momentum
3830  pt(1)=0.
3831  pt(2)=0.
3832  pt(3)=0.
3833  pt(4)=amtau
3834 c
3835  CALL ranmar(rrr,9)
3836 c
3837 * masses of 4, 3 and 2 pi systems
3838 c 3 pi with sampling for resonance
3839 cam
3840  rr1=rrr(6)
3841  ams1=(amp1+amp2+amp3+amp4)**2
3842  ams2=(amtau-amnuta)**2
3843  alp1=atan((ams1-amrop**2)/amrop/gamrop)
3844  alp2=atan((ams2-amrop**2)/amrop/gamrop)
3845  alp=alp1+rr1*(alp2-alp1)
3846  am4sq =amrop**2+amrop*gamrop*tan(alp)
3847  am4 =sqrt(am4sq)
3848  phspac=phspac*
3849  $ ((am4sq-amrop**2)**2+(amrop*gamrop)**2)/(amrop*gamrop)
3850  phspac=phspac*(alp2-alp1)
3851 
3852 c
3853  rr1=rrr(1)
3854  ams1=(amp2+amp3+amp4)**2
3855  ams2=(am4-amp1)**2
3856  IF (rrr(9).GT.prez) THEN
3857  am3sq=ams1+ rr1*(ams2-ams1)
3858  am3 =sqrt(am3sq)
3859 c --- this part of jacobian will be recovered later
3860  ff1=ams2-ams1
3861  ELSE
3862 * phase space with sampling for omega resonance,
3863  alp1=atan((ams1-amrx**2)/amrx/gamrx)
3864  alp2=atan((ams2-amrx**2)/amrx/gamrx)
3865  alp=alp1+rr1*(alp2-alp1)
3866  am3sq =amrx**2+amrx*gamrx*tan(alp)
3867  am3 =sqrt(am3sq)
3868 c --- this part of the jacobian will be recovered later ---------------
3869  ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
3870  ff1=ff1*(alp2-alp1)
3871  ENDIF
3872 c mass of 2
3873  rr2=rrr(2)
3874  ams1=(amp3+amp4)**2
3875  ams2=(am3-amp2)**2
3876 * flat phase space;
3877  am2sq=ams1+ rr2*(ams2-ams1)
3878  am2 =sqrt(am2sq)
3879 c --- this part of jacobian will be recovered later
3880  ff2=(ams2-ams1)
3881 * 2 restframe, define piz and pipl
3882  enq1=(am2sq-amp3**2+amp4**2)/(2*am2)
3883  enq2=(am2sq+amp3**2-amp4**2)/(2*am2)
3884  ppi= enq1**2-amp4**2
3885  pppi=sqrt(abs(enq1**2-amp4**2))
3886  phspac=phspac*(4*pi)*(2*pppi/am2)
3887 * piz momentum in 2 rest frame
3888  CALL sphera(pppi,piz)
3889  piz(4)=enq1
3890 * pipl momentum in 2 rest frame
3891  DO 30 i=1,3
3892  30 pipl(i)=-piz(i)
3893  pipl(4)=enq2
3894 * 3 rest frame, define pim1
3895 * pr momentum
3896  pr(1)=0
3897  pr(2)=0
3898  pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
3899  pr(3)= sqrt(abs(pr(4)**2-am2**2))
3900  ppi = pr(4)**2-am2**2
3901 * pim1 momentum
3902  pim1(1)=0
3903  pim1(2)=0
3904  pim1(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
3905  pim1(3)=-pr(3)
3906 c --- this part of jacobian will be recovered later
3907  ff3=(4*pi)*(2*pr(3)/am3)
3908 * old pions boosted from 2 rest frame to 3 rest frame
3909  exe=(pr(4)+pr(3))/am2
3910  CALL bostr3(exe,piz,piz)
3911  CALL bostr3(exe,pipl,pipl)
3912  rr3=rrr(3)
3913  rr4=rrr(4)
3914  thet =acos(-1.+2*rr3)
3915  phi = 2*pi*rr4
3916  CALL rotpol(thet,phi,pipl)
3917  CALL rotpol(thet,phi,pim1)
3918  CALL rotpol(thet,phi,piz)
3919  CALL rotpol(thet,phi,pr)
3920 * 4 rest frame, define pim2
3921 * pr momentum
3922  pr(1)=0
3923  pr(2)=0
3924  pr(4)=1./(2*am4)*(am4**2+am3**2-amp1**2)
3925  pr(3)= sqrt(abs(pr(4)**2-am3**2))
3926  ppi = pr(4)**2-am3**2
3927 * pim2 momentum
3928  pim2(1)=0
3929  pim2(2)=0
3930  pim2(4)=1./(2*am4)*(am4**2-am3**2+amp1**2)
3931  pim2(3)=-pr(3)
3932 c --- this part of jacobian will be recovered later
3933  ff4=(4*pi)*(2*pr(3)/am4)
3934 * old pions boosted from 3 rest frame to 4 rest frame
3935  exe=(pr(4)+pr(3))/am3
3936  CALL bostr3(exe,piz,piz)
3937  CALL bostr3(exe,pipl,pipl)
3938  CALL bostr3(exe,pim1,pim1)
3939  rr3=rrr(7)
3940  rr4=rrr(8)
3941  thet =acos(-1.+2*rr3)
3942  phi = 2*pi*rr4
3943  CALL rotpol(thet,phi,pipl)
3944  CALL rotpol(thet,phi,pim1)
3945  CALL rotpol(thet,phi,pim2)
3946  CALL rotpol(thet,phi,piz)
3947  CALL rotpol(thet,phi,pr)
3948 c
3949 * now to the tau rest frame, define paa and neutrino momenta
3950 * paa momentum
3951  paa(1)=0
3952  paa(2)=0
3953  paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am4**2)
3954  paa(3)= sqrt(abs(paa(4)**2-am4**2))
3955  ppi = paa(4)**2-am4**2
3956  phspac=phspac*(4*pi)*(2*paa(3)/amtau)
3957  phsp=phsp*(4*pi)*(2*paa(3)/amtau)
3958 * tau-neutrino momentum
3959  pn(1)=0
3960  pn(2)=0
3961  pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am4**2)
3962  pn(3)=-paa(3)
3963 c zbw 20.12.2002 bug fix
3964  IF(rrr(9).LE.0.5*prez) THEN
3965  DO 72 i=1,4
3966  x=pim1(i)
3967  pim1(i)=pim2(i)
3968  72 pim2(i)=x
3969  ENDIF
3970 c end of bug fix
3971 c we include remaining part of the jacobian
3972 c --- flat channel
3973  am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
3974  $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
3975  ams2=(am4-amp2)**2
3976  ams1=(amp1+amp3+amp4)**2
3977  ff1=(ams2-ams1)
3978  ams1=(amp3+amp4)**2
3979  ams2=(sqrt(am3sq)-amp1)**2
3980  ff2=ams2-ams1
3981  ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
3982  ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
3983  uu=ff1*ff2*ff3*ff4
3984 c --- first channel
3985  am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
3986  $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
3987  ams2=(am4-amp2)**2
3988  ams1=(amp1+amp3+amp4)**2
3989  alp1=atan((ams1-amrx**2)/amrx/gamrx)
3990  alp2=atan((ams2-amrx**2)/amrx/gamrx)
3991  ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
3992  ff1=ff1*(alp2-alp1)
3993  ams1=(amp3+amp4)**2
3994  ams2=(sqrt(am3sq)-amp1)**2
3995  ff2=ams2-ams1
3996  ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
3997  ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
3998  ff=ff1*ff2*ff3*ff4
3999 c --- second channel
4000  am3sq=(pim2(4)+piz(4)+pipl(4))**2-(pim2(3)+piz(3)+pipl(3))**2
4001  $ -(pim2(2)+piz(2)+pipl(2))**2-(pim2(1)+piz(1)+pipl(1))**2
4002  ams2=(am4-amp1)**2
4003  ams1=(amp2+amp3+amp4)**2
4004  alp1=atan((ams1-amrx**2)/amrx/gamrx)
4005  alp2=atan((ams2-amrx**2)/amrx/gamrx)
4006  gg1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4007  gg1=gg1*(alp2-alp1)
4008  ams1=(amp3+amp4)**2
4009  ams2=(sqrt(am3sq)-amp2)**2
4010  gg2=ams2-ams1
4011  gg3=(4*pi)*(xlam(am2**2,amp2**2,am3sq)/am3sq)
4012  gg4=(4*pi)*(xlam(am3sq,amp1**2,am4**2)/am4**2)
4013  gg=gg1*gg2*gg3*gg4
4014 c --- jacobian averaged over the two
4015  IF ( ( (ff+gg)*uu+ff*gg ).GT.0.0d0) THEN
4016  rr=ff*gg*uu/(0.5*prez*(ff+gg)*uu+(1.0-prez)*ff*gg)
4017  phspac=phspac*rr
4018  ELSE
4019  phspac=0.0
4020  ENDIF
4021 * momenta of the two pi-minus are randomly symmetrised
4022  IF (jnpi.EQ.1) THEN
4023  rr5= rrr(5)
4024  IF(rr5.LE.0.5) THEN
4025  DO 70 i=1,4
4026  x=pim1(i)
4027  pim1(i)=pim2(i)
4028  70 pim2(i)=x
4029  ENDIF
4030  phspac=phspac/2.
4031  ELSE
4032 c momenta of pi0-s are generated uniformly only IF prez=0.0
4033  rr5= rrr(5)
4034  IF(rr5.LE.0.5) THEN
4035  DO 71 i=1,4
4036  x=pim1(i)
4037  pim1(i)=pim2(i)
4038  71 pim2(i)=x
4039  ENDIF
4040  phspac=phspac/6.
4041  ENDIF
4042 * all pions boosted from 4 rest frame to tau rest frame
4043 * z-axis antiparallel to neutrino momentum
4044  exe=(paa(4)+paa(3))/am4
4045  CALL bostr3(exe,piz,piz)
4046  CALL bostr3(exe,pipl,pipl)
4047  CALL bostr3(exe,pim1,pim1)
4048  CALL bostr3(exe,pim2,pim2)
4049  CALL bostr3(exe,pr,pr)
4050 c partial width consists of phase space and amplitude
4051 c check on consistency with dadnpi, THEN, code breakes uniform pion
4052 c distribution in hadronic system
4053 cam assume neutrino mass=0. and sum over final polarisation
4054 c amx2=am4**2
4055 c brak= 2*(amtau**2-amx2) * (amtau**2+2.*amx2)
4056 c amplit=ccabib**2*gfermi**2/2. * brak * amx2*sigee(amx2,1)
4057  IF (jnpi.EQ.1) THEN
4058  CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4059  ELSEIF (jnpi.EQ.2) THEN
4060  CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv)
4061  ENDIF
4062  dgamt=1/(2.*amtau)*amplit*phspac
4063 c phase space check
4064 c dgamt=phspac
4065  DO 77 k=1,4
4066  pmult(k,1)=pim1(k)
4067  pmult(k,2)=pim2(k)
4068  pmult(k,3)=pipl(k)
4069  pmult(k,4)=piz(k)
4070  77 CONTINUE
4071  END
4072  SUBROUTINE dam4pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,AMPLIT,HV)
4073 c ----------------------------------------------------------------------
4074 * calculates differential cross section and polarimeter vector
4075 * for tau decay into 4 pi modes
4076 * all spin effects in the full decay chain are taken into account.
4077 * calculations done in tau rest frame with z-axis along neutrino moment
4078 c mnum decay mode identifier.
4079 c
4080 c called by : dphsaa
4081 c ----------------------------------------------------------------------
4082  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4083  * ,ampiz,ampi,amro,gamro,ama1,gama1
4084  * ,amk,amkz,amkst,gamkst
4085 c
4086  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4087  * ,ampiz,ampi,amro,gamro,ama1,gama1
4088  * ,amk,amkz,amkst,gamkst
4089  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4090  real*4 gfermi,gv,ga,ccabib,scabib,gamel
4091  REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4),PIM4(4)
4092  REAL PIVEC(4),PIAKS(4),HVM(4)
4093  COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5
4094  EXTERNAL form1,form2,form3,form4,form5
4095  DATA pi /3.141592653589793238462643/
4096  DATA icont /0/
4097 c
4098  CALL curr_cleo(mnum,pim1,pim2,pim3,pim4,hadcur)
4099 c
4100 * calculate pi-vectors: vector and axial
4101  CALL clvec(hadcur,pn,pivec)
4102  CALL claxi(hadcur,pn,piaks)
4103  CALL clnut(hadcur,brakm,hvm)
4104 * spin independent part of decay diff-cross-sect. in tau rest frame
4105  brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4106  & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4107  amplit=(ccabib*gfermi)**2*brak/2.
4108 c polarimeter vector in tau rest frame
4109  DO 90 i=1,3
4110  hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4111  & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4112 c hv is defined for tau- with gamma=b+hv*pol
4113  IF (brak.NE.0.0)
4114  &hv(i)=-hv(i)/brak
4115  90 CONTINUE
4116  END
4117  SUBROUTINE dph5pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4118 c ----------------------------------------------------------------------
4119 * it simulates 5pi decay in tau rest frame with
4120 * z-axis along 5pi momentum
4121 c ----------------------------------------------------------------------
4122  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4123  * ,ampiz,ampi,amro,gamro,ama1,gama1
4124  * ,amk,amkz,amkst,gamkst
4125 c
4126  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4127  * ,ampiz,ampi,amro,gamro,ama1,gama1
4128 
4129 
4130  * ,amk,amkz,amkst,gamkst
4131  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4132  real*4 gfermi,gv,ga,ccabib,scabib,gamel
4133  parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4134  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4135  & ,names
4136  CHARACTER NAMES(NMODE)*31
4137  REAL HV(4),PT(4),PN(4),PAA(4),PMULT(4,9)
4138  real*4 pr(4),pi1(4),pi2(4),pi3(4),pi4(4),pi5(4)
4139  real*8 amp1,amp2,amp3,amp4,amp5,ams1,ams2,amom,gamom
4140  real*8 am5sq,am4sq,am3sq,am2sq,am5,am4,am3
4141  real*4 rrr(10)
4142  real*8 gg1,gg2,gg3,ff1,ff2,ff3,ff4,alp,alp1,alp2
4143  real*8 xm,am,gamma
4144 ccm.s.>>>>>>
4145  real*8 phspac
4146 ccm.s.<<<<<<
4147  DATA pi /3.141592653589793238462643/
4148  DATA icont /0/
4149  data fpi /93.3e-3/
4150 c
4151  COMPLEX BWIGN
4152 c
4153  bwign(xm,am,gamma)=xm**2/cmplx(xm**2-am**2,gamma*am)
4154 
4155 c
4156  amom=.782
4157  gamom=0.0085
4158 c
4159 c 6 body phase space normalised as in bjorken-drell
4160 c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
4161  phspac=1./2**29/pi**14
4162 c phspac=1./2**5/pi**2
4163 c init 5pi decay mode(jnpi)
4164  amp1=dcdmas(idffin(1,jnpi))
4165  amp2=dcdmas(idffin(2,jnpi))
4166  amp3=dcdmas(idffin(3,jnpi))
4167  amp4=dcdmas(idffin(4,jnpi))
4168  amp5=dcdmas(idffin(5,jnpi))
4169 c
4170 c tau momentum
4171  pt(1)=0.
4172  pt(2)=0.
4173  pt(3)=0.
4174  pt(4)=amtau
4175 c
4176  CALL ranmar(rrr,10)
4177 c
4178 c masses of 5, 4, 3 and 2 pi systems
4179 c 3 pi with sampling for omega resonance
4180 cam
4181 c mass of 5 (12345)
4182  rr1=rrr(10)
4183  ams1=(amp1+amp2+amp3+amp4+amp5)**2
4184  ams2=(amtau-amnuta)**2
4185  am5sq=ams1+ rr1*(ams2-ams1)
4186  am5 =sqrt(am5sq)
4187  phspac=phspac*(ams2-ams1)
4188 c
4189 c mass of 4 (2345)
4190 c flat phase space
4191  rr1=rrr(9)
4192  ams1=(amp2+amp3+amp4+amp5)**2
4193  ams2=(am5-amp1)**2
4194  am4sq=ams1+ rr1*(ams2-ams1)
4195  am4 =sqrt(am4sq)
4196  gg1=ams2-ams1
4197 c
4198 c mass of 3 (234)
4199 c phase space with sampling for omega resonance
4200  rr1=rrr(1)
4201  ams1=(amp2+amp3+amp4)**2
4202  ams2=(am4-amp5)**2
4203  alp1=atan((ams1-amom**2)/amom/gamom)
4204  alp2=atan((ams2-amom**2)/amom/gamom)
4205  alp=alp1+rr1*(alp2-alp1)
4206  am3sq =amom**2+amom*gamom*tan(alp)
4207  am3 =sqrt(am3sq)
4208 c --- this part of the jacobian will be recovered later ---------------
4209  gg2=((am3sq-amom**2)**2+(amom*gamom)**2)/(amom*gamom)
4210  gg2=gg2*(alp2-alp1)
4211 c flat phase space;
4212 c am3sq=ams1+ rr1*(ams2-ams1)
4213 c am3 =sqrt(am3sq)
4214 c --- this part of jacobian will be recovered later
4215 c gg2=ams2-ams1
4216 c
4217 c mass of 2 (34)
4218  rr2=rrr(2)
4219  ams1=(amp3+amp4)**2
4220  ams2=(am3-amp2)**2
4221 c flat phase space;
4222  am2sq=ams1+ rr2*(ams2-ams1)
4223  am2 =sqrt(am2sq)
4224 c --- this part of jacobian will be recovered later
4225  gg3=ams2-ams1
4226 c
4227 c(34) restframe, define pi3 and pi4
4228  enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4229  enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4230  ppi= enq1**2-amp3**2
4231  pppi=sqrt(abs(enq1**2-amp3**2))
4232  ff1=(4*pi)*(2*pppi/am2)
4233 c pi3 momentum in(34) rest frame
4234  call sphera(pppi,pi3)
4235  pi3(4)=enq1
4236 c pi4 momentum in(34) rest frame
4237  do 30 i=1,3
4238  30 pi4(i)=-pi3(i)
4239  pi4(4)=enq2
4240 c
4241 c(234) rest frame, define pi2
4242 c pr momentum
4243  pr(1)=0
4244  pr(2)=0
4245  pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4246  pr(3)= sqrt(abs(pr(4)**2-am2**2))
4247  ppi = pr(4)**2-am2**2
4248 c pi2 momentum
4249  pi2(1)=0
4250  pi2(2)=0
4251  pi2(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4252  pi2(3)=-pr(3)
4253 c --- this part of jacobian will be recovered later
4254  ff2=(4*pi)*(2*pr(3)/am3)
4255 c old pions boosted from 2 rest frame to 3 rest frame
4256  exe=(pr(4)+pr(3))/am2
4257  call bostr3(exe,pi3,pi3)
4258  call bostr3(exe,pi4,pi4)
4259  rr3=rrr(3)
4260  rr4=rrr(4)
4261  thet =acos(-1.+2*rr3)
4262  phi = 2*pi*rr4
4263  call rotpol(thet,phi,pi2)
4264  call rotpol(thet,phi,pi3)
4265  call rotpol(thet,phi,pi4)
4266 c
4267 c(2345) rest frame, define pi5
4268 c pr momentum
4269  pr(1)=0
4270  pr(2)=0
4271  pr(4)=1./(2*am4)*(am4**2+am3**2-amp5**2)
4272  pr(3)= sqrt(abs(pr(4)**2-am3**2))
4273  ppi = pr(4)**2-am3**2
4274 c pi5 momentum
4275  pi5(1)=0
4276  pi5(2)=0
4277  pi5(4)=1./(2*am4)*(am4**2-am3**2+amp5**2)
4278  pi5(3)=-pr(3)
4279 c --- this part of jacobian will be recovered later
4280  ff3=(4*pi)*(2*pr(3)/am4)
4281 c old pions boosted from 3 rest frame to 4 rest frame
4282  exe=(pr(4)+pr(3))/am3
4283  call bostr3(exe,pi2,pi2)
4284  call bostr3(exe,pi3,pi3)
4285  call bostr3(exe,pi4,pi4)
4286  rr3=rrr(5)
4287  rr4=rrr(6)
4288  thet =acos(-1.+2*rr3)
4289  phi = 2*pi*rr4
4290  call rotpol(thet,phi,pi2)
4291  call rotpol(thet,phi,pi3)
4292  call rotpol(thet,phi,pi4)
4293  call rotpol(thet,phi,pi5)
4294 c
4295 c(12345) rest frame, define pi1
4296 c pr momentum
4297  pr(1)=0
4298  pr(2)=0
4299  pr(4)=1./(2*am5)*(am5**2+am4**2-amp1**2)
4300  pr(3)= sqrt(abs(pr(4)**2-am4**2))
4301  ppi = pr(4)**2-am4**2
4302 c pi1 momentum
4303  pi1(1)=0
4304  pi1(2)=0
4305  pi1(4)=1./(2*am5)*(am5**2-am4**2+amp1**2)
4306  pi1(3)=-pr(3)
4307 c --- this part of jacobian will be recovered later
4308  ff4=(4*pi)*(2*pr(3)/am5)
4309 c old pions boosted from 4 rest frame to 5 rest frame
4310  exe=(pr(4)+pr(3))/am4
4311  call bostr3(exe,pi2,pi2)
4312  call bostr3(exe,pi3,pi3)
4313  call bostr3(exe,pi4,pi4)
4314  call bostr3(exe,pi5,pi5)
4315  rr3=rrr(7)
4316  rr4=rrr(8)
4317  thet =acos(-1.+2*rr3)
4318  phi = 2*pi*rr4
4319  call rotpol(thet,phi,pi1)
4320  call rotpol(thet,phi,pi2)
4321  call rotpol(thet,phi,pi3)
4322  call rotpol(thet,phi,pi4)
4323  call rotpol(thet,phi,pi5)
4324 c
4325 * now to the tau rest frame, define paa and neutrino momenta
4326 * paa momentum
4327  paa(1)=0
4328  paa(2)=0
4329 c paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5**2)
4330 c paa(3)= sqrt(abs(paa(4)**2-am5**2))
4331 c ppi = paa(4)**2-am5**2
4332  paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5sq)
4333  paa(3)= sqrt(abs(paa(4)**2-am5sq))
4334  ppi = paa(4)**2-am5sq
4335  phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4336 * tau-neutrino momentum
4337  pn(1)=0
4338  pn(2)=0
4339  pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am5**2)
4340  pn(3)=-paa(3)
4341 c
4342  phspac=phspac * gg1*gg2*gg3*ff1*ff2*ff3*ff4
4343 c
4344 c all pions boosted from 5 rest frame to tau rest frame
4345 c z-axis antiparallel to neutrino momentum
4346  exe=(paa(4)+paa(3))/am5
4347  call bostr3(exe,pi1,pi1)
4348  call bostr3(exe,pi2,pi2)
4349  call bostr3(exe,pi3,pi3)
4350  call bostr3(exe,pi4,pi4)
4351  call bostr3(exe,pi5,pi5)
4352 c
4353 c partial width consists of phase space and amplitude
4354 c amplitude(cf ys.tsai phys.rev.d4,2821(1971)
4355 c or f.gilman sh.rhie phys.rev.d31,1066(1985)
4356 c
4357  pxq=amtau*paa(4)
4358  pxn=amtau*pn(4)
4359  qxn=paa(4)*pn(4)-paa(1)*pn(1)-paa(2)*pn(2)-paa(3)*pn(3)
4360  brak=2*(gv**2+ga**2)*(2*pxq*qxn+am5sq*pxn)
4361  & -6*(gv**2-ga**2)*amtau*amnuta*am5sq
4362  fompp = cabs(bwign(am3,amom,gamom))**2
4363 c normalisation factor(to some numerical undimensioned factor;
4364 c cf r.fischer et al zphys c3, 313 (1980))
4365  fnorm = 1/fpi**6
4366 c amplit=ccabib**2*gfermi**2/2. * brak * am5sq*sigee(am5sq,jnpi)
4367  amplit=ccabib**2*gfermi**2/2. * brak
4368  amplit = amplit * fompp * fnorm
4369 c phase space test
4370 c amplit = amplit * fnorm
4371  dgamt=1/(2.*amtau)*amplit*phspac
4372 c ignore spin terms
4373  DO 40 i=1,3
4374  40 hv(i)=0.
4375 c
4376  do 77 k=1,4
4377  pmult(k,1)=pi1(k)
4378  pmult(k,2)=pi2(k)
4379  pmult(k,3)=pi3(k)
4380  pmult(k,4)=pi4(k)
4381  pmult(k,5)=pi5(k)
4382  77 continue
4383  return
4384 c missing: transposition of identical particles, startistical factors
4385 c for identical matrices, polarimetric vector. matrix element rather naive.
4386 c flat phase space in pion system + with breit wigner for omega
4387 c anyway it is better than nothing, and code is improvable.
4388  end
4389  SUBROUTINE dphsrk(DGAMT,HV,PN,PR,PMULT,INUM)
4390 c ----------------------------------------------------------------------
4391 c it simulates rho decay in tau rest frame with
4392 c z-axis along rho momentum
4393 c rho decays to k kbar
4394 c ----------------------------------------------------------------------
4395  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4396  * ,ampiz,ampi,amro,gamro,ama1,gama1
4397  * ,amk,amkz,amkst,gamkst
4398 c
4399  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4400  * ,ampiz,ampi,amro,gamro,ama1,gama1
4401  * ,amk,amkz,amkst,gamkst
4402  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4403  real*4 gfermi,gv,ga,ccabib,scabib,gamel
4404  REAL HV(4),PT(4),PN(4),PR(4),PKC(4),PKZ(4),QQ(4),PMULT(4,9)
4405  REAL RR1(1)
4406  DATA pi /3.141592653589793238462643/
4407  DATA icont /0/
4408 c
4409 c three body phase space normalised as in bjorken-drell
4410  phspac=1./2**11/pi**5
4411 c tau momentum
4412  pt(1)=0.
4413  pt(2)=0.
4414  pt(3)=0.
4415  pt(4)=amtau
4416 c mass of(real/virtual) rho
4417  ams1=(amk+amkz)**2
4418  ams2=(amtau-amnuta)**2
4419 c flat phase space
4420  CALL ranmar(rr1,1)
4421  amx2=ams1+ rr1(1)*(ams2-ams1)
4422  amx=sqrt(amx2)
4423  phspac=phspac*(ams2-ams1)
4424 c phase space with sampling for rho resonance
4425 c alp1=atan((ams1-amro**2)/amro/gamro)
4426 c alp2=atan((ams2-amro**2)/amro/gamro)
4427 cam
4428  100 CONTINUE
4429 c CALL ranmar(rr1,1)
4430 c alp=alp1+rr1(1)*(alp2-alp1)
4431 c amx2=amro**2+amro*gamro*tan(alp)
4432 c amx=sqrt(amx2)
4433 c IF(amx.LT.(amk+amkz)) GO TO 100
4434 cam
4435 c phspac=phspac*((amx2-amro**2)**2+(amro*gamro)**2)/(amro*gamro)
4436 c phspac=phspac*(alp2-alp1)
4437 c
4438 c tau-neutrino momentum
4439  pn(1)=0
4440  pn(2)=0
4441  pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
4442  pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
4443 c rho momentum
4444  pr(1)=0
4445  pr(2)=0
4446  pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
4447  pr(3)=-pn(3)
4448  phspac=phspac*(4*pi)*(2*pr(3)/amtau)
4449 c
4450 cam
4451  enq1=(amx2+amk**2-amkz**2)/(2.*amx)
4452  enq2=(amx2-amk**2+amkz**2)/(2.*amx)
4453  pppi=sqrt((enq1-amk)*(enq1+amk))
4454  phspac=phspac*(4*pi)*(2*pppi/amx)
4455 c charged pi momentum in rho rest frame
4456  CALL sphera(pppi,pkc)
4457  pkc(4)=enq1
4458 c neutral pi momentum in rho rest frame
4459  DO 20 i=1,3
4460 20 pkz(i)=-pkc(i)
4461  pkz(4)=enq2
4462  exe=(pr(4)+pr(3))/amx
4463 c pions boosted from rho rest frame to tau rest frame
4464  CALL bostr3(exe,pkc,pkc)
4465  CALL bostr3(exe,pkz,pkz)
4466  DO 30 i=1,4
4467  30 qq(i)=pkc(i)-pkz(i)
4468 c qq transverse to pr
4469  pksd =pr(4)*pr(4)-pr(3)*pr(3)-pr(2)*pr(2)-pr(1)*pr(1)
4470  qqpks=pr(4)* qq(4)-pr(3)* qq(3)-pr(2)* qq(2)-pr(1)* qq(1)
4471  DO 31 i=1,4
4472 31 qq(i)=qq(i)-pr(i)*qqpks/pksd
4473 c amplitude
4474  prodpq=pt(4)*qq(4)
4475  prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
4476  prodpn=pt(4)*pn(4)
4477  qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
4478  brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
4479  & +(gv**2-ga**2)*amtau*amnuta*qq2
4480  amplit=(gfermi*ccabib)**2*brak*2*fpirk(amx)
4481  dgamt=1/(2.*amtau)*amplit*phspac
4482  DO 40 i=1,3
4483  40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
4484  do 77 k=1,4
4485  pmult(k,1)=pkc(k)
4486  pmult(k,2)=pkz(k)
4487  77 continue
4488  RETURN
4489  END
4490  FUNCTION fpirk(W)
4491 c ----------------------------------------------------------
4492 c square of pion form factor
4493 c ----------------------------------------------------------
4494  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4495  * ,ampiz,ampi,amro,gamro,ama1,gama1
4496  * ,amk,amkz,amkst,gamkst
4497 c
4498  real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4499  * ,ampiz,ampi,amro,gamro,ama1,gama1
4500  * ,amk,amkz,amkst,gamkst
4501 c COMPLEX FPIKMK
4502  COMPLEX FPIKM
4503  fpirk=cabs(fpikm(w,amk,amkz))**2
4504 c fpirk=cabs(fpikmk(w,amk,amkz))**2
4505  END
4506  COMPLEX FUNCTION fpikmk(W,XM1,XM2)
4507 c **********************************************************
4508 c kaon form factor
4509 c **********************************************************
4510  COMPLEX BWIGM
4511  REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
4512  EXTERNAL bwig
4513  DATA init /0/
4514 c
4515 c ------------ parameters --------------------
4516  IF (init.EQ.0 ) THEN
4517  init=1
4518  pi=3.141592654
4519  pim=.140
4520  rom=0.773
4521  rog=0.145
4522  rom1=1.570
4523  rog1=0.510
4524 c beta1=-0.111
4525  beta1=-0.221
4526  ENDIF
4527 c -----------------------------------------------
4528  s=w**2
4529  fpikmk=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
4530  & /(1+beta1)
4531  RETURN
4532  END
4533  SUBROUTINE reslux
4534 c ****************
4535 c initialize lund COMMON
4536  nhep=0
4537  END
4538  SUBROUTINE dwrph(KTO,PHX)
4539 c
4540 c -------------------------
4541 c
4542  IMPLICIT REAL*8 (a-h,o-z)
4543  real*4 phx(4)
4544  real*4 qhot(4)
4545 c
4546  DO 9 k=1,4
4547  qhot(k) =0.0
4548  9 CONTINUE
4549 c CASE of tau radiative decays.
4550 c filling of the lund COMMON block.
4551  DO 1002 i=1,4
4552  1002 qhot(i)=phx(i)
4553  IF (qhot(4).GT.1.e-5) CALL dwluph(kto,qhot)
4554  RETURN
4555  END
4556  SUBROUTINE dwluph(KTO,PHOT)
4557 c---------------------------------------------------------------------
4558 c lorentz transformation to cmsystem and
4559 c updating of hepevt record
4560 c
4561 c called by : dexay1,(dekay1,dekay2)
4562 c
4563 c used when radiative corrections in decays are generated
4564 c---------------------------------------------------------------------
4565 c
4566  REAL PHOT(4)
4567  COMMON /taupos/ np1,np2
4568 c
4569 c check energy
4570  IF (phot(4).LE.0.0) RETURN
4571 c
4572 c position of decaying particle:
4573  IF((kto.EQ. 1).OR.(kto.EQ.11)) THEN
4574  nps=np1
4575  ELSE
4576  nps=np2
4577  ENDIF
4578 c
4579  ktos=kto
4580  IF(ktos.GT.10) ktos=ktos-10
4581 c boost and append photon(gamma is 22)
4582  CALL tralo4(ktos,phot,phot,am)
4583  CALL filhep(0,1,22,nps,nps,0,0,phot,0.0,.true.)
4584 c
4585  RETURN
4586  END
4587 
4588  SUBROUTINE dwluel(KTO,ISGN,PNU,PWB,PEL,PNE)
4589 c ----------------------------------------------------------------------
4590 c lorentz transformation to cmsystem and
4591 c updating of hepevt record
4592 c
4593 c isgn = 1/-1 for tau-/tau+
4594 c
4595 c called by : dexay,(dekay1,dekay2)
4596 c ----------------------------------------------------------------------
4597 c
4598  REAL PNU(4),PWB(4),PEL(4),PNE(4)
4599  COMMON /taupos/ np1,np2
4600 c
4601 c position of decaying particle:
4602  IF(kto.EQ. 1) THEN
4603  nps=np1
4604  ELSE
4605  nps=np2
4606  ENDIF
4607 c
4608 c tau neutrino(nu_tau is 16)
4609  CALL tralo4(kto,pnu,pnu,am)
4610  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4611 c
4612 c w boson(w+ is 24)
4613  CALL tralo4(kto,pwb,pwb,am)
4614 c CALL filhep(0,2,-24*isgn,nps,nps,0,0,pwb,am,.true.)
4615 c
4616 c electron(e- is 11)
4617  CALL tralo4(kto,pel,pel,am)
4618  CALL filhep(0,1,11*isgn,nps,nps,0,0,pel,am,.false.)
4619 c
4620 c anti electron neutrino(nu_e is 12)
4621  CALL tralo4(kto,pne,pne,am)
4622  CALL filhep(0,1,-12*isgn,nps,nps,0,0,pne,am,.true.)
4623 c
4624  RETURN
4625  END
4626  SUBROUTINE dwlumu(KTO,ISGN,PNU,PWB,PMU,PNM)
4627 c ----------------------------------------------------------------------
4628 c lorentz transformation to cmsystem and
4629 c updating of hepevt record
4630 c
4631 c isgn = 1/-1 for tau-/tau+
4632 c
4633 c called by : dexay,(dekay1,dekay2)
4634 c ----------------------------------------------------------------------
4635 c
4636  REAL PNU(4),PWB(4),PMU(4),PNM(4)
4637  COMMON /taupos/ np1,np2
4638 c
4639 c position of decaying particle:
4640  IF(kto.EQ. 1) THEN
4641  nps=np1
4642  ELSE
4643  nps=np2
4644  ENDIF
4645 c
4646 c tau neutrino(nu_tau is 16)
4647  CALL tralo4(kto,pnu,pnu,am)
4648  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4649 c
4650 c w boson(w+ is 24)
4651  CALL tralo4(kto,pwb,pwb,am)
4652 c CALL filhep(0,2,-24*isgn,nps,nps,0,0,pwb,am,.true.)
4653 c
4654 c muon(mu- is 13)
4655  CALL tralo4(kto,pmu,pmu,am)
4656  CALL filhep(0,1,13*isgn,nps,nps,0,0,pmu,am,.false.)
4657 c
4658 c anti muon neutrino(nu_mu is 14)
4659  CALL tralo4(kto,pnm,pnm,am)
4660  CALL filhep(0,1,-14*isgn,nps,nps,0,0,pnm,am,.true.)
4661 c
4662  RETURN
4663  END
4664  SUBROUTINE dwlupi(KTO,ISGN,PPI,PNU)
4665 c ----------------------------------------------------------------------
4666 c lorentz transformation to cmsystem and
4667 c updating of hepevt record
4668 c
4669 c isgn = 1/-1 for tau-/tau+
4670 c
4671 c called by : dexay,(dekay1,dekay2)
4672 c ----------------------------------------------------------------------
4673 c
4674  REAL PNU(4),PPI(4)
4675  COMMON /taupos/ np1,np2
4676 c
4677 c position of decaying particle:
4678  IF(kto.EQ. 1) THEN
4679  nps=np1
4680  ELSE
4681  nps=np2
4682  ENDIF
4683 c
4684 c tau neutrino(nu_tau is 16)
4685  CALL tralo4(kto,pnu,pnu,am)
4686  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4687 c
4688 c charged pi meson(pi+ is 211)
4689  CALL tralo4(kto,ppi,ppi,am)
4690  CALL filhep(0,1,-211*isgn,nps,nps,0,0,ppi,am,.true.)
4691 c
4692  RETURN
4693  END
4694  SUBROUTINE dwluro(KTO,ISGN,PNU,PRHO,PIC,PIZ)
4695 c ----------------------------------------------------------------------
4696 c lorentz transformation to cmsystem and
4697 c updating of hepevt record
4698 c
4699 c isgn = 1/-1 for tau-/tau+
4700 c
4701 c called by : dexay,(dekay1,dekay2)
4702 c ----------------------------------------------------------------------
4703 c
4704  REAL PNU(4),PRHO(4),PIC(4),PIZ(4)
4705  COMMON /taupos/ np1,np2
4706 c
4707 c position of decaying particle:
4708  IF(kto.EQ. 1) THEN
4709  nps=np1
4710  ELSE
4711  nps=np2
4712  ENDIF
4713 c
4714 c tau neutrino(nu_tau is 16)
4715  CALL tralo4(kto,pnu,pnu,am)
4716  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4717 c
4718 c charged rho meson(rho+ is 213)
4719  CALL tralo4(kto,prho,prho,am)
4720  CALL filhep(0,2,-213*isgn,nps,nps,0,0,prho,am,.true.)
4721 c
4722 c charged pi meson(pi+ is 211)
4723  CALL tralo4(kto,pic,pic,am)
4724  CALL filhep(0,1,-211*isgn,-1,-1,0,0,pic,am,.true.)
4725 c
4726 c pi0 meson(pi0 is 111)
4727  CALL tralo4(kto,piz,piz,am)
4728  CALL filhep(0,1,111,-2,-2,0,0,piz,am,.true.)
4729 c
4730  RETURN
4731  END
4732  SUBROUTINE dwluaa(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
4733 c ----------------------------------------------------------------------
4734 c lorentz transformation to cmsystem and
4735 c updating of hepevt record
4736 c
4737 c isgn = 1/-1 for tau-/tau+
4738 c jaa = 1 (2) for a_1- decay to pi+ 2pi- (pi- 2pi0)
4739 c
4740 c called by : dexay,(dekay1,dekay2)
4741 c ----------------------------------------------------------------------
4742 c
4743  REAL PNU(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
4744  COMMON /taupos/ np1,np2
4745 c
4746 c position of decaying particle:
4747  IF(kto.EQ. 1) THEN
4748  nps=np1
4749  ELSE
4750  nps=np2
4751  ENDIF
4752 c
4753 c tau neutrino(nu_tau is 16)
4754  CALL tralo4(kto,pnu,pnu,am)
4755  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4756 c
4757 c charged a_1 meson(a_1+ is 20213)
4758  CALL tralo4(kto,paa,paa,am)
4759  CALL filhep(0,1,-20213*isgn,nps,nps,0,0,paa,am,.true.)
4760 c
4761 c two possible decays of the charged a1 meson
4762  IF(jaa.EQ.1) THEN
4763 c
4764 c a1 --> pi+ pi- pi- (or charged conjugate)
4765 c
4766 c pi minus(or c.c.) (pi+ is 211)
4767  CALL tralo4(kto,pim2,pim2,am)
4768  CALL filhep(0,1,-211*isgn,-1,-1,0,0,pim2,am,.true.)
4769 c
4770 c pi minus(or c.c.) (pi+ is 211)
4771  CALL tralo4(kto,pim1,pim1,am)
4772  CALL filhep(0,1,-211*isgn,-2,-2,0,0,pim1,am,.true.)
4773 c
4774 c pi plus(or c.c.) (pi+ is 211)
4775  CALL tralo4(kto,pipl,pipl,am)
4776  CALL filhep(0,1, 211*isgn,-3,-3,0,0,pipl,am,.true.)
4777 c
4778  ELSE IF (jaa.EQ.2) THEN
4779 c
4780 c a1 --> pi- pi0 pi0(or charged conjugate)
4781 c
4782 c pi zero(pi0 is 111)
4783  CALL tralo4(kto,pim2,pim2,am)
4784  CALL filhep(0,1,111,-1,-1,0,0,pim2,am,.true.)
4785 c
4786 c pi zero(pi0 is 111)
4787  CALL tralo4(kto,pim1,pim1,am)
4788  CALL filhep(0,1,111,-2,-2,0,0,pim1,am,.true.)
4789 c
4790 c pi minus(or c.c.) (pi+ is 211)
4791  CALL tralo4(kto,pipl,pipl,am)
4792  CALL filhep(0,1,-211*isgn,-3,-3,0,0,pipl,am,.true.)
4793 c
4794  ENDIF
4795 c
4796  RETURN
4797  END
4798  SUBROUTINE dwlukk (KTO,ISGN,PKK,PNU)
4799 c ----------------------------------------------------------------------
4800 c lorentz transformation to cmsystem and
4801 c updating of hepevt record
4802 c
4803 c isgn = 1/-1 for tau-/tau+
4804 c
4805 c ----------------------------------------------------------------------
4806 c
4807  REAL PKK(4),PNU(4)
4808  COMMON /taupos/ np1,np2
4809 c
4810 c position of decaying particle
4811  IF (kto.EQ.1) THEN
4812  nps=np1
4813  ELSE
4814  nps=np2
4815  ENDIF
4816 c
4817 c tau neutrino(nu_tau is 16)
4818  CALL tralo4 (kto,pnu,pnu,am)
4819  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4820 c
4821 c k meson(k+ is 321)
4822  CALL tralo4 (kto,pkk,pkk,am)
4823  CALL filhep(0,1,-321*isgn,nps,nps,0,0,pkk,am,.true.)
4824 c
4825  RETURN
4826  END
4827  SUBROUTINE dwluks(KTO,ISGN,PNU,PKS,PKK,PPI,JKST)
4828  COMMON / taukle / bra1,brk0,brk0b,brks
4829  real*4 bra1,brk0,brk0b,brks
4830 c ----------------------------------------------------------------------
4831 c lorentz transformation to cmsystem and
4832 c updating of hepevt record
4833 c
4834 c isgn = 1/-1 for tau-/tau+
4835 c jkst=10 (20) corresponds to k0b pi- (k- pi0) decay
4836 c
4837 c ----------------------------------------------------------------------
4838 c
4839  REAL PNU(4),PKS(4),PKK(4),PPI(4),XIO(1)
4840  COMMON /taupos/ np1,np2
4841 c
4842 c position of decaying particle
4843  IF(kto.EQ. 1) THEN
4844  nps=np1
4845  ELSE
4846  nps=np2
4847  ENDIF
4848 c
4849 c tau neutrino(nu_tau is 16)
4850  CALL tralo4(kto,pnu,pnu,am)
4851  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4852 c
4853 c charged k* meson(k*+ is 323)
4854  CALL tralo4(kto,pks,pks,am)
4855  CALL filhep(0,1,-323*isgn,nps,nps,0,0,pks,am,.true.)
4856 c
4857 c two possible decay modes of charged k*
4858  IF(jkst.EQ.10) THEN
4859 c
4860 c k*- --> pi- k0b(or charged conjugate)
4861 c
4862 c charged pi meson(pi+ is 211)
4863  CALL tralo4(kto,ppi,ppi,am)
4864  CALL filhep(0,1,-211*isgn,-1,-1,0,0,ppi,am,.true.)
4865 c
4866  bran=brk0b
4867  IF (isgn.EQ.-1) bran=brk0
4868 c k0 --> k0_long(is 130) / k0_short(is 310) = 1/1
4869  CALL ranmar(xio,1)
4870  IF(xio(1).GT.bran) THEN
4871  k0type = 130
4872  ELSE
4873  k0type = 310
4874  ENDIF
4875 c
4876  CALL tralo4(kto,pkk,pkk,am)
4877  CALL filhep(0,1,k0type,-2,-2,0,0,pkk,am,.true.)
4878 c
4879  ELSE IF(jkst.EQ.20) THEN
4880 c
4881 c k*- --> pi0 k-
4882 c
4883 c pi zero(pi0 is 111)
4884  CALL tralo4(kto,ppi,ppi,am)
4885  CALL filhep(0,1,111,-1,-1,0,0,ppi,am,.true.)
4886 c
4887 c charged k meson(k+ is 321)
4888  CALL tralo4(kto,pkk,pkk,am)
4889  CALL filhep(0,1,-321*isgn,-2,-2,0,0,pkk,am,.true.)
4890 c
4891  ENDIF
4892 c
4893  RETURN
4894  END
4895  SUBROUTINE dwlnew(KTO,ISGN,PNU,PWB,PNPI,MODE)
4896 c ----------------------------------------------------------------------
4897 c lorentz transformation to cmsystem and
4898 c updating of hepevt record
4899 c
4900 c isgn = 1/-1 for tau-/tau+
4901 c
4902 c called by : dexay,(dekay1,dekay2)
4903 c ----------------------------------------------------------------------
4904 c
4905  parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4906  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4907  & ,names
4908  COMMON /taupos/ np1,np2
4909  CHARACTER NAMES(NMODE)*31
4910  REAL PNU(4),PWB(4),PNPI(4,9)
4911  REAL PPI(4)
4912 c
4913  jnpi=mode-7
4914 c position of decaying particle
4915  IF(kto.EQ. 1) THEN
4916  nps=np1
4917  ELSE
4918  nps=np2
4919  ENDIF
4920 c
4921 c tau neutrino(nu_tau is 16)
4922  CALL tralo4(kto,pnu,pnu,am)
4923  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4924 c
4925 c w boson(w+ is 24)
4926  CALL tralo4(kto,pwb,pwb,am)
4927  CALL filhep(0,1,-24*isgn,nps,nps,0,0,pwb,am,.true.)
4928 c
4929 c multi pi mode jnpi
4930 c
4931 c get multiplicity of mode jnpi
4932  nd=mulpik(jnpi)
4933  DO i=1,nd
4934  kfpi=lunpik(idffin(i,jnpi),-isgn)
4935 c for charged conjugate case, change charged pions only
4936 c IF(kfpi.NE.111)kfpi=kfpi*isgn
4937  DO j=1,4
4938  ppi(j)=pnpi(j,i)
4939  END DO
4940  CALL tralo4(kto,ppi,ppi,am)
4941  CALL filhep(0,1,kfpi,-i,-i,0,0,ppi,am,.true.)
4942  END DO
4943 c
4944  RETURN
4945  END
4946  FUNCTION amast(PP)
4947 c ----------------------------------------------------------------------
4948 c calculates mass of pp(double precision)
4949 c
4950 c used by : radkor
4951 c ----------------------------------------------------------------------
4952  IMPLICIT REAL*8 (a-h,o-z)
4953  real*8 pp(4)
4954  aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
4955 c
4956  IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
4957  amast=aaa
4958  RETURN
4959  END
4960  FUNCTION amas4(PP)
4961 c ******************
4962 c ----------------------------------------------------------------------
4963 c calculates mass of pp
4964 c
4965 c used by :
4966 c ----------------------------------------------------------------------
4967  REAL PP(4)
4968  aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
4969  IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
4970  amas4=aaa
4971  RETURN
4972  END
4973  FUNCTION angxy(X,Y)
4974 c ----------------------------------------------------------------------
4975 c
4976 c used by : koralz radkor
4977 c ----------------------------------------------------------------------
4978  IMPLICIT DOUBLE PRECISION (a-h,o-z)
4979  DATA pi /3.141592653589793238462643d0/
4980 c
4981  IF(abs(y).LT.abs(x)) THEN
4982  the=atan(abs(y/x))
4983  IF(x.LE.0d0) the=pi-the
4984  ELSE
4985  the=acos(x/sqrt(x**2+y**2))
4986  ENDIF
4987  angxy=the
4988  RETURN
4989  END
4990  FUNCTION angfi(X,Y)
4991 c ----------------------------------------------------------------------
4992 * calculates angle in(0,2*pi) range out of x-y
4993 c
4994 c used by : koralz radkor
4995 c ----------------------------------------------------------------------
4996  IMPLICIT DOUBLE PRECISION (a-h,o-z)
4997  DATA pi /3.141592653589793238462643d0/
4998 c
4999  IF(abs(y).LT.abs(x)) THEN
5000  the=atan(abs(y/x))
5001  IF(x.LE.0d0) the=pi-the
5002  ELSE
5003  the=acos(x/sqrt(x**2+y**2))
5004  ENDIF
5005  IF(y.LT.0d0) the=2d0*pi-the
5006  angfi=the
5007  END
5008  SUBROUTINE rotod1(PH1,PVEC,QVEC)
5009 c ----------------------------------------------------------------------
5010 c
5011 c used by : koralz
5012 c ----------------------------------------------------------------------
5013  IMPLICIT DOUBLE PRECISION (a-h,o-z)
5014  dimension pvec(4),qvec(4),rvec(4)
5015 c
5016  phi=ph1
5017  cs=cos(phi)
5018  sn=sin(phi)
5019  DO 10 i=1,4
5020  10 rvec(i)=pvec(i)
5021  qvec(1)=rvec(1)
5022  qvec(2)= cs*rvec(2)-sn*rvec(3)
5023  qvec(3)= sn*rvec(2)+cs*rvec(3)
5024  qvec(4)=rvec(4)
5025  RETURN
5026  END
5027  SUBROUTINE rotod2(PH1,PVEC,QVEC)
5028 c ----------------------------------------------------------------------
5029 c
5030 c used by : koralz radkor
5031 c ----------------------------------------------------------------------
5032  IMPLICIT DOUBLE PRECISION (a-h,o-z)
5033  dimension pvec(4),qvec(4),rvec(4)
5034 c
5035  phi=ph1
5036  cs=cos(phi)
5037  sn=sin(phi)
5038  DO 10 i=1,4
5039  10 rvec(i)=pvec(i)
5040  qvec(1)= cs*rvec(1)+sn*rvec(3)
5041  qvec(2)=rvec(2)
5042  qvec(3)=-sn*rvec(1)+cs*rvec(3)
5043  qvec(4)=rvec(4)
5044  RETURN
5045  END
5046  SUBROUTINE rotod3(PH1,PVEC,QVEC)
5047 c ----------------------------------------------------------------------
5048 c
5049 c used by : koralz radkor
5050 c ----------------------------------------------------------------------
5051  IMPLICIT DOUBLE PRECISION (a-h,o-z)
5052 c
5053  dimension pvec(4),qvec(4),rvec(4)
5054  phi=ph1
5055  cs=cos(phi)
5056  sn=sin(phi)
5057  DO 10 i=1,4
5058  10 rvec(i)=pvec(i)
5059  qvec(1)= cs*rvec(1)-sn*rvec(2)
5060  qvec(2)= sn*rvec(1)+cs*rvec(2)
5061  qvec(3)=rvec(3)
5062  qvec(4)=rvec(4)
5063  END
5064  SUBROUTINE bostr3(EXE,PVEC,QVEC)
5065 c ----------------------------------------------------------------------
5066 c boost along z axis, exe=exp(eta), eta= hiperbolic velocity.
5067 c
5068 c used by : tauola koralz(?)
5069 c ----------------------------------------------------------------------
5070  real*4 pvec(4),qvec(4),rvec(4)
5071 c
5072  DO 10 i=1,4
5073  10 rvec(i)=pvec(i)
5074  rpl=rvec(4)+rvec(3)
5075  rmi=rvec(4)-rvec(3)
5076  qpl=rpl*exe
5077  qmi=rmi/exe
5078  qvec(1)=rvec(1)
5079  qvec(2)=rvec(2)
5080  qvec(3)=(qpl-qmi)/2
5081  qvec(4)=(qpl+qmi)/2
5082  END
5083  SUBROUTINE bostd3(EXE,PVEC,QVEC)
5084 c ----------------------------------------------------------------------
5085 c boost along z axis, exe=exp(eta), eta= hiperbolic velocity.
5086 c
5087 c used by : koralz radkor
5088 c ----------------------------------------------------------------------
5089  IMPLICIT DOUBLE PRECISION (a-h,o-z)
5090  dimension pvec(4),qvec(4),rvec(4)
5091 c
5092  DO 10 i=1,4
5093  10 rvec(i)=pvec(i)
5094  rpl=rvec(4)+rvec(3)
5095  rmi=rvec(4)-rvec(3)
5096  qpl=rpl*exe
5097  qmi=rmi/exe
5098  qvec(1)=rvec(1)
5099  qvec(2)=rvec(2)
5100  qvec(3)=(qpl-qmi)/2
5101  qvec(4)=(qpl+qmi)/2
5102  RETURN
5103  END
5104  SUBROUTINE rotor1(PH1,PVEC,QVEC)
5105 c ----------------------------------------------------------------------
5106 c
5107 c called by :
5108 c ----------------------------------------------------------------------
5109  real*4 pvec(4),qvec(4),rvec(4)
5110 c
5111  phi=ph1
5112  cs=cos(phi)
5113  sn=sin(phi)
5114  DO 10 i=1,4
5115  10 rvec(i)=pvec(i)
5116  qvec(1)=rvec(1)
5117  qvec(2)= cs*rvec(2)-sn*rvec(3)
5118  qvec(3)= sn*rvec(2)+cs*rvec(3)
5119  qvec(4)=rvec(4)
5120  END
5121  SUBROUTINE rotor2(PH1,PVEC,QVEC)
5122 c ----------------------------------------------------------------------
5123 c
5124 c used by : tauola
5125 c ----------------------------------------------------------------------
5126  IMPLICIT REAL*4(a-h,o-z)
5127  real*4 pvec(4),qvec(4),rvec(4)
5128 c
5129  phi=ph1
5130  cs=cos(phi)
5131  sn=sin(phi)
5132  DO 10 i=1,4
5133  10 rvec(i)=pvec(i)
5134  qvec(1)= cs*rvec(1)+sn*rvec(3)
5135  qvec(2)=rvec(2)
5136  qvec(3)=-sn*rvec(1)+cs*rvec(3)
5137  qvec(4)=rvec(4)
5138  END
5139  SUBROUTINE rotor3(PHI,PVEC,QVEC)
5140 c ----------------------------------------------------------------------
5141 c
5142 c used by : tauola
5143 c ----------------------------------------------------------------------
5144  real*4 pvec(4),qvec(4),rvec(4)
5145 c
5146  cs=cos(phi)
5147  sn=sin(phi)
5148  DO 10 i=1,4
5149  10 rvec(i)=pvec(i)
5150  qvec(1)= cs*rvec(1)-sn*rvec(2)
5151  qvec(2)= sn*rvec(1)+cs*rvec(2)
5152  qvec(3)=rvec(3)
5153  qvec(4)=rvec(4)
5154  END
5155  SUBROUTINE spherd(R,X)
5156 c ----------------------------------------------------------------------
5157 c generates uniformly three-vector x on sphere of radius r
5158 c double precison version of sphera
5159 c ----------------------------------------------------------------------
5160  real*8 r,x(4),pi,costh,sinth
5161  real*4 rrr(2)
5162  DATA pi /3.141592653589793238462643d0/
5163 c
5164  CALL ranmar(rrr,2)
5165  costh=-1+2*rrr(1)
5166  sinth=sqrt(1 -costh**2)
5167  x(1)=r*sinth*cos(2*pi*rrr(2))
5168  x(2)=r*sinth*sin(2*pi*rrr(2))
5169  x(3)=r*costh
5170  RETURN
5171  END
5172  SUBROUTINE rotpox(THET,PHI,PP)
5173  IMPLICIT REAL*8 (a-h,o-z)
5174 c ----------------------------------------------------------------------
5175 c
5176 c ----------------------------------------------------------------------
5177  dimension pp(4)
5178 c
5179  CALL rotod2(thet,pp,pp)
5180  CALL rotod3( phi,pp,pp)
5181  RETURN
5182  END
5183  SUBROUTINE sphera(R,X)
5184 c ----------------------------------------------------------------------
5185 c generates uniformly three-vector x on sphere of radius r
5186 c
5187 c called by : dphsxx,dadmpi,dadmkk
5188 c ----------------------------------------------------------------------
5189  REAL X(4)
5190  real*4 rrr(2)
5191  DATA pi /3.141592653589793238462643/
5192 c
5193  CALL ranmar(rrr,2)
5194  costh=-1.+2.*rrr(1)
5195  sinth=sqrt(1.-costh**2)
5196  x(1)=r*sinth*cos(2*pi*rrr(2))
5197  x(2)=r*sinth*sin(2*pi*rrr(2))
5198  x(3)=r*costh
5199  RETURN
5200  END
5201  SUBROUTINE rotpol(THET,PHI,PP)
5202 c ----------------------------------------------------------------------
5203 c
5204 c called by : dadmaa,dphsaa
5205 c ----------------------------------------------------------------------
5206  REAL PP(4)
5207 c
5208  CALL rotor2(thet,pp,pp)
5209  CALL rotor3( phi,pp,pp)
5210  RETURN
5211  END
5212  SUBROUTINE ranmar(RVEC,LENV)
5213 c ----------------------------------------------------------------------
5214 c<<<<<FUNCTION ranmar(IDUMM)
5215 c cernlib v113, version with automatic default initialization
5216 c transformed to SUBROUTINE to be as in cernlib
5217 c am.lutz november 1988, feb. 1989
5218 c
5219 c!Universal random number generator proposed by Marsaglia and Zaman
5220 c in report fsu-scri-87-50
5221 c modified by f. james, 1988 and 1989, to generate a vector
5222 c of pseudorandom numbers rvec of length lenv, and to put in
5223 c the COMMON block everything needed to specify currrent state,
5224 c and to add input and output entry points rmarin, rmarut.
5225 c
5226 c unique random number used in the program
5227 c ----------------------------------------------------------------------
5228  COMMON / inout / inut,iout
5229  dimension rvec(*)
5230  common/raset1/u(97),c,i97,j97
5231  parameter(modcns=1000000000)
5232  DATA ntot,ntot2,ijkl/-1,0,0/
5233 c
5234  IF (ntot .GE. 0) GO TO 50
5235 c
5236 c default initialization. user has called ranmar without rmarin.
5237  ijkl = 54217137
5238  ntot = 0
5239  ntot2 = 0
5240  kalled = 0
5241  GO TO 1
5242 c
5243  entry rmarin(ijklin, ntotin,ntot2n)
5244 c initializing routine for ranmar, may be called before
5245 c generating pseudorandom numbers with ranmar. the input
5246 c values should be in the ranges: 0<=ijklin<=900 ooo ooo
5247 c 0<=ntotin<=999 999 999
5248 c 0<=ntot2n<<999 999 999!
5249 c to get the standard values in marsaglia-s paper, ijklin=54217137
5250 c ntotin,ntot2n=0
5251  ijkl = ijklin
5252  ntot = max(ntotin,0)
5253  ntot2= max(ntot2n,0)
5254  kalled = 1
5255 c always come here to initialize
5256  1 CONTINUE
5257  ij = ijkl/30082
5258  kl = ijkl - 30082*ij
5259  i = mod(ij/177, 177) + 2
5260  j = mod(ij, 177) + 2
5261  k = mod(kl/169, 178) + 1
5262  l = mod(kl, 169)
5263  WRITE(iout,201) ijkl,ntot,ntot2
5264  201 FORMAT(1x,' RANMAR INITIALIZED: ',i10,2x,2i10)
5265  DO 2 ii= 1, 97
5266  s = 0.
5267  t = .5
5268  DO 3 jj= 1, 24
5269  m = mod(mod(i*j,179)*k, 179)
5270  i = j
5271  j = k
5272  k = m
5273  l = mod(53*l+1, 169)
5274  IF (mod(l*m,64) .GE. 32) s = s+t
5275  3 t = 0.5*t
5276  2 u(ii) = s
5277  twom24 = 1.0
5278  DO 4 i24= 1, 24
5279  4 twom24 = 0.5*twom24
5280  c = 362436.*twom24
5281  cd = 7654321.*twom24
5282  cm = 16777213.*twom24
5283  i97 = 97
5284  j97 = 33
5285 c complete initialization by skipping
5286 c(ntot2*modcns + ntot) random numbers
5287  DO 45 loop2= 1, ntot2+1
5288  now = modcns
5289  IF (loop2 .EQ. ntot2+1) now=ntot
5290  IF (now .GT. 0) THEN
5291  WRITE (iout,'(A,I15)') ' RMARIN SKIPPING OVER ',now
5292  DO 40 idum = 1, ntot
5293  uni = u(i97)-u(j97)
5294  IF (uni .LT. 0.) uni=uni+1.
5295  u(i97) = uni
5296  i97 = i97-1
5297  IF (i97 .EQ. 0) i97=97
5298  j97 = j97-1
5299  IF (j97 .EQ. 0) j97=97
5300  c = c - cd
5301  IF (c .LT. 0.) c=c+cm
5302  40 CONTINUE
5303  ENDIF
5304  45 CONTINUE
5305  IF (kalled .EQ. 1) RETURN
5306 c
5307 c normal entry to generate lenv random numbers
5308  50 CONTINUE
5309  DO 100 ivec= 1, lenv
5310  uni = u(i97)-u(j97)
5311  IF (uni .LT. 0.) uni=uni+1.
5312  u(i97) = uni
5313  i97 = i97-1
5314  IF (i97 .EQ. 0) i97=97
5315  j97 = j97-1
5316  IF (j97 .EQ. 0) j97=97
5317  c = c - cd
5318  IF (c .LT. 0.) c=c+cm
5319  uni = uni-c
5320  IF (uni .LT. 0.) uni=uni+1.
5321 c replace exact zeroes by uniform distr. *2**-24
5322  IF (uni .EQ. 0.) THEN
5323  uni = twom24*u(2)
5324 c an exact zero here is very unlikely, but lets be safe.
5325  IF (uni .EQ. 0.) uni= twom24*twom24
5326  ENDIF
5327  rvec(ivec) = uni
5328  100 CONTINUE
5329  ntot = ntot + lenv
5330  IF (ntot .GE. modcns) THEN
5331  ntot2 = ntot2 + 1
5332  ntot = ntot - modcns
5333  ENDIF
5334  RETURN
5335 c entry to output current status
5336  entry rmarut(ijklut,ntotut,ntot2t)
5337  ijklut = ijkl
5338  ntotut = ntot
5339  ntot2t = ntot2
5340  RETURN
5341  END
5342  FUNCTION dilogt(X)
5343 c *****************
5344  IMPLICIT REAL*8(a-h,o-z)
5345 cern c304 version 29/07/71 dilog 59 c
5346  z=-1.64493406684822
5347  IF(x .LT.-1.0) GO TO 1
5348  IF(x .LE. 0.5) GO TO 2
5349  IF(x .EQ. 1.0) GO TO 3
5350  IF(x .LE. 2.0) GO TO 4
5351  z=3.2898681336964
5352  1 t=1.0/x
5353  s=-0.5
5354  z=z-0.5* log(abs(x))**2
5355  GO TO 5
5356  2 t=x
5357  s=0.5
5358  z=0.
5359  GO TO 5
5360  3 dilogt=1.64493406684822
5361  RETURN
5362  4 t=1.0-x
5363  s=-0.5
5364  z=1.64493406684822 - log(x)* log(abs(t))
5365  5 y=2.66666666666666 *t+0.66666666666666
5366  b= 0.00000 00000 00001
5367  a=y*b +0.00000 00000 00004
5368  b=y*a-b+0.00000 00000 00011
5369  a=y*b-a+0.00000 00000 00037
5370  b=y*a-b+0.00000 00000 00121
5371  a=y*b-a+0.00000 00000 00398
5372  b=y*a-b+0.00000 00000 01312
5373  a=y*b-a+0.00000 00000 04342
5374  b=y*a-b+0.00000 00000 14437
5375  a=y*b-a+0.00000 00000 48274
5376  b=y*a-b+0.00000 00001 62421
5377  a=y*b-a+0.00000 00005 50291
5378  b=y*a-b+0.00000 00018 79117
5379  a=y*b-a+0.00000 00064 74338
5380  b=y*a-b+0.00000 00225 36705
5381  a=y*b-a+0.00000 00793 87055
5382  b=y*a-b+0.00000 02835 75385
5383  a=y*b-a+0.00000 10299 04264
5384  b=y*a-b+0.00000 38163 29463
5385  a=y*b-a+0.00001 44963 00557
5386  b=y*a-b+0.00005 68178 22718
5387  a=y*b-a+0.00023 20021 96094
5388  b=y*a-b+0.00100 16274 96164
5389  a=y*b-a+0.00468 63619 59447
5390  b=y*a-b+0.02487 93229 24228
5391  a=y*b-a+0.16607 30329 27855
5392  a=y*a-b+1.93506 43008 6996
5393  dilogt=s*t*(a-b)+z
5394  RETURN
5395 c=======================================================================
5396 c===================END OF CPC PART ====================================
5397 c=======================================================================
5398  END
5399 
5400 c-----------------------------------------------------------------------
5401 c initialize rchl currents
5402 c dummy routine for compatibility with new updates of tauola
5403 c
5404 c added 25.jul.2012
5405 c-----------------------------------------------------------------------
5406  SUBROUTINE inirchl(FLAG)
5407  INTEGER FLAG
5408 
5409  IF(flag.NE.0) THEN
5410  WRITE(*,25) flag
5411  25 FORMAT(1x, "TAUOLA IniRChL: Fatal error, FLAG=",i2," but RChL currents missing")
5412  WRITE(*,*) " in loaded version of TAUOLA-FORTRAN library."
5413  stop
5414  ENDIF
5415 
5416  RETURN
5417  END
Generated by doxygen 1.9.1