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