C++InterfacetoTauola
BBB/prod/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 ** 1: 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)) & - HJ(2)*AIMAG(HJ(2)) - HJ(1)*AIMAG(HJ(1)) )
3333  RETURN
3334  END
3335  SUBROUTINE dampog(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3336 C ----------------------------------------------------------------------
3337 * CALCULATES DIFFERENTIAL CROSS SECTION AND POLARIMETER VECTOR
3338 * FOR TAU DECAY INTO A1, A1 DECAYS NEXT INTO RHO+PI AND RHO INTO PI+PI.
3339 * ALL SPIN EFFECTS IN THE FULL DECAY CHAIN ARE TAKEN INTO ACCOUNT.
3340 * CALCULATIONS DONE IN TAU REST FRAME WITH Z-AXIS ALONG NEUTRINO MOMENT
3341 * THE ROUTINE IS WRITEN FOR ZERO NEUTRINO MASS.
3342 C
3343 
3344 C called by : DPHSAA
3345 
3346 C ----------------------------------------------------------------------
3347  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3348  * ,ampiz,ampi,amro,gamro,ama1,gama1
3349  * ,amk,amkz,amkst,gamkst
3350 C
3351  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
3352  * ,ampiz,ampi,amro,gamro,ama1,gama1
3353  * ,amk,amkz,amkst,gamkst
3354  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3355  REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
3356  COMMON /testa1/ keya1
3357  REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4)
3358  REAL PAA(4),VEC1(4),VEC2(4)
3359  REAL PIVEC(4),PIAKS(4),HVM(4)
3360  COMPLEX BWIGN,HADCUR(4),FNORM,FORMOM
3361  DATA icont /1/
3362 * THIS INLINE FUNCT. CALCULATES THE SCALAR PART OF THE PROPAGATOR
3363 
3364 C AJWMOD to satisfy compiler, comment out this unused function.
3365 
3366 C
3367 * FOUR MOMENTUM OF A1
3368  DO 10 i=1,4
3369  vec1(i)=0.0
3370  vec2(i)=0.0
3371  hv(i) =0.0
3372  10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3373  vec1(1)=1.0
3374 * MASSES OF A1, AND OF TWO PI-PAIRS WHICH MAY FORM RHO
3375  xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3376  xmom =sqrt(abs( (pim2(4)+pipl(4))**2-(pim2(3)+pipl(3))**2
3377  $ -(pim2(2)+pipl(2))**2-(pim2(1)+pipl(1))**2 ))
3378  xmro2 =(pipl(1))**2 +(pipl(2))**2 +(pipl(3))**2
3379 * ELEMENTS OF HADRON CURRENT
3380  prod1 =vec1(1)*pipl(1)
3381  prod2 =vec2(2)*pipl(2)
3382  p12 =pim1(4)*pim2(4)-pim1(1)*pim2(1)
3383  $ -pim1(2)*pim2(2)-pim1(3)*pim2(3)
3384  p1pl =pim1(4)*pipl(4)-pim1(1)*pipl(1)
3385  $ -pim1(2)*pipl(2)-pim1(3)*pipl(3)
3386  p2pl =pipl(4)*pim2(4)-pipl(1)*pim2(1)
3387  $ -pipl(2)*pim2(2)-pipl(3)*pim2(3)
3388  DO 40 i=1,3
3389  vec1(i)= (vec1(i)-prod1/xmro2*pipl(i))
3390  40 CONTINUE
3391  gnorm=sqrt(vec1(1)**2+vec1(2)**2+vec1(3)**2)
3392  DO 41 i=1,3
3393  vec1(i)= vec1(i)/gnorm
3394  41 CONTINUE
3395  vec2(1)=(vec1(2)*pipl(3)-vec1(3)*pipl(2))/sqrt(xmro2)
3396  vec2(2)=(vec1(3)*pipl(1)-vec1(1)*pipl(3))/sqrt(xmro2)
3397  vec2(3)=(vec1(1)*pipl(2)-vec1(2)*pipl(1))/sqrt(xmro2)
3398  p1vec1 =pim1(4)*vec1(4)-pim1(1)*vec1(1)
3399  $ -pim1(2)*vec1(2)-pim1(3)*vec1(3)
3400  p2vec1 =vec1(4)*pim2(4)-vec1(1)*pim2(1)
3401  $ -vec1(2)*pim2(2)-vec1(3)*pim2(3)
3402  p1vec2 =pim1(4)*vec2(4)-pim1(1)*vec2(1)
3403  $ -pim1(2)*vec2(2)-pim1(3)*vec2(3)
3404  p2vec2 =vec2(4)*pim2(4)-vec2(1)*pim2(1)
3405  $ -vec2(2)*pim2(2)-vec2(3)*pim2(3)
3406 * HADRON CURRENT
3407  fnorm=formom(xmaa,xmom)
3408  brak=0.0
3409  DO 120 jj=1,2
3410  DO 45 i=1,4
3411  IF (jj.EQ.1) THEN
3412  hadcur(i) = fnorm *(
3413  $ vec1(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3414  $ -pim2(i)*(p2vec1*p1pl-p1vec1*p2pl)
3415  $ +pipl(i)*(p2vec1*p12 -p1vec1*(ampi**2+p2pl)) )
3416  ELSE
3417  hadcur(i) = fnorm *(
3418  $ vec2(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3419  $ -pim2(i)*(p2vec2*p1pl-p1vec2*p2pl)
3420  $ +pipl(i)*(p2vec2*p12 -p1vec2*(ampi**2+p2pl)) )
3421  ENDIF
3422  45 CONTINUE
3423 C
3424 * CALCULATE PI-VECTORS: VECTOR AND AXIAL
3425  CALL clvec(hadcur,pn,pivec)
3426  CALL claxi(hadcur,pn,piaks)
3427  CALL clnut(hadcur,brakm,hvm)
3428 * SPIN INDEPENDENT PART OF DECAY DIFF-CROSS-SECT. IN TAU REST FRAME
3429  brak=brak+(gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3430  & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3431  DO 90 i=1,3
3432  hv(i)=hv(i)-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3433  & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3434  90 CONTINUE
3435 C HV IS DEFINED FOR TAU- WITH GAMMA=B+HV*POL
3436  120 CONTINUE
3437  amplit=(gfermi*ccabib)**2*brak/2.
3438 C THE STATISTICAL FACTOR FOR IDENTICAL PI-S WAS CANCELLED WITH
3439 C TWO, FOR TWO MODES OF A1 DECAY NAMELLY PI+PI-PI- AND PI-PI0PI0
3440 C POLARIMETER VECTOR IN TAU REST FRAME
3441  DO 91 i=1,3
3442  hv(i)=-hv(i)/brak
3443  91 CONTINUE
3444 
3445  END
3446  SUBROUTINE damppk(MNUM,PT,PN,PIM1,PIM2,PIM3,AMPLIT,HV)
3447 C ----------------------------------------------------------------------
3448 * CALCULATES DIFFERENTIAL CROSS SECTION AND POLARIMETER VECTOR
3449 * FOR TAU DECAY INTO K K pi, K pi pi.
3450 * ALL SPIN EFFECTS IN THE FULL DECAY CHAIN ARE TAKEN INTO ACCOUNT.
3451 * CALCULATIONS DONE IN TAU REST FRAME WITH Z-AXIS ALONG NEUTRINO MOMENT
3452 C MNUM DECAY MODE IDENTIFIER.
3453 C
3454 
3455 C called by : DPHSAA
3456 
3457 C ----------------------------------------------------------------------
3458  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3459  * ,ampiz,ampi,amro,gamro,ama1,gama1
3460  * ,amk,amkz,amkst,gamkst
3461 C
3462  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
3463  * ,ampiz,ampi,amro,gamro,ama1,gama1
3464  * ,amk,amkz,amkst,gamkst
3465  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3466  REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
3467  REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4)
3468  REAL PAA(4),VEC1(4),VEC2(4),VEC3(4),VEC4(4),VEC5(4)
3469  REAL PIVEC(4),PIAKS(4),HVM(4)
3470  REAL FNORM(0:19),COEF(1:5,0:19)
3471  COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5,UROJ
3472 
3473  COMPLEX F1,F2,F3,F4,F5
3474 
3475  EXTERNAL form1,form2,form3,form4,form5
3476  DATA pi /3.141592653589793238462643/
3477  DATA icont /0/
3478 C
3479  DATA fpi /93.3e-3/
3480  IF (icont.EQ.0) THEN
3481  icont=1
3482  uroj=cmplx(0.0,1.0)
3483  dwapi0=sqrt(2.0)
3484  fnorm(0)=ccabib/fpi
3485  fnorm(1)=ccabib/fpi
3486  fnorm(2)=ccabib/fpi
3487  fnorm(3)=ccabib/fpi
3488  fnorm(4)=scabib/fpi/dwapi0
3489  fnorm(5)=scabib/fpi
3490  fnorm(6)=scabib/fpi
3491  fnorm(7)=ccabib/fpi
3492  fnorm(8)=0.0 ! this chanel is dead
3493  fnorm(9)=ccabib/fpi
3494  DO k=10,19
3495  fnorm(k)=fnorm(9) ! these chanells are not initialized
3496  ENDDO
3497 C
3498  coef(1,0)= 2.0*sqrt(2.)/3.0
3499  coef(2,0)=-2.0*sqrt(2.)/3.0
3500 
3501 C AJW 2/98: Add in the D-wave and I=0 3pi substructure:
3502  coef(3,0)= 2.0*sqrt(2.)/3.0
3503 
3504  coef(4,0)= fpi
3505  coef(5,0)= 0.0
3506 C
3507  coef(1,1)=-sqrt(2.)/3.0
3508  coef(2,1)= sqrt(2.)/3.0
3509  coef(3,1)= 0.0
3510  coef(4,1)= fpi
3511  coef(5,1)= sqrt(2.)
3512 C
3513  coef(1,2)=-sqrt(2.)/3.0
3514  coef(2,2)= sqrt(2.)/3.0
3515  coef(3,2)= 0.0
3516  coef(4,2)= 0.0
3517  coef(5,2)=-sqrt(2.)
3518 C
3519 
3520 C AJW 11/97: Add in the K*-prim-s, ala Finkemeier&Mirkes
3521  coef(1,3)= 1./3.
3522  coef(2,3)=-2./3.
3523  coef(3,3)= 2./3.
3524 
3525  coef(4,3)= 0.0
3526  coef(5,3)= 0.0
3527 C
3528  coef(1,4)= 1.0/sqrt(2.)/3.0
3529  coef(2,4)=-1.0/sqrt(2.)/3.0
3530  coef(3,4)= 0.0
3531  coef(4,4)= 0.0
3532  coef(5,4)= 0.0
3533 C
3534  coef(1,5)=-sqrt(2.)/3.0
3535  coef(2,5)= sqrt(2.)/3.0
3536  coef(3,5)= 0.0
3537  coef(4,5)= 0.0
3538  coef(5,5)=-sqrt(2.)
3539 C
3540 
3541 C AJW 11/97: Add in the K*-prim-s, ala Finkemeier&Mirkes
3542  coef(1,6)= 1./3.
3543  coef(2,6)=-2./3.
3544  coef(3,6)= 2./3.
3545 
3546  coef(4,6)= 0.0
3547  coef(5,6)=-2.0
3548 C
3549  coef(1,7)= 0.0
3550  coef(2,7)= 0.0
3551  coef(3,7)= 0.0
3552  coef(4,7)= 0.0
3553  coef(5,7)=-sqrt(2.0/3.0)
3554 C
3555  coef(1,8)= 0.0
3556  coef(2,8)= 0.0
3557  coef(3,8)= 0.0
3558  coef(4,8)= 0.0
3559  coef(5,8)= 0.0
3560 C
3561 C
3562  coef(1,9)= 2.0*sqrt(2.)/3.0
3563  coef(2,9)=-2.0*sqrt(2.)/3.0
3564  coef(3,9)= 2.0*sqrt(2.)/3.0
3565 
3566  coef(4,9)= fpi
3567  coef(5,9)= 0.0
3568 
3569  DO k=10,19 ! these chanells are not initialized
3570  coef(1,k)= 2.0*sqrt(2.)/3.0
3571  coef(2,k)=-2.0*sqrt(2.)/3.0
3572  coef(3,k)= 2.0*sqrt(2.)/3.0
3573 
3574  coef(4,k)= fpi
3575  coef(5,k)= 0.0
3576 
3577  ENDDO
3578 
3579 
3580  ENDIF
3581 C
3582  DO 10 i=1,4
3583  10 paa(i)=pim1(i)+pim2(i)+pim3(i)
3584  xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3585  xmro1 =sqrt(abs((pim3(4)+pim2(4))**2-(pim3(1)+pim2(1))**2
3586  $ -(pim3(2)+pim2(2))**2-(pim3(3)+pim2(3))**2))
3587  xmro2 =sqrt(abs((pim3(4)+pim1(4))**2-(pim3(1)+pim1(1))**2
3588  $ -(pim3(2)+pim1(2))**2-(pim3(3)+pim1(3))**2))
3589  xmro3 =sqrt(abs((pim1(4)+pim2(4))**2-(pim1(1)+pim2(1))**2
3590  $ -(pim1(2)+pim2(2))**2-(pim1(3)+pim2(3))**2))
3591 * ELEMENTS OF HADRON CURRENT
3592  prod1 =paa(4)*(pim2(4)-pim3(4))-paa(1)*(pim2(1)-pim3(1))
3593  $ -paa(2)*(pim2(2)-pim3(2))-paa(3)*(pim2(3)-pim3(3))
3594  prod2 =paa(4)*(pim3(4)-pim1(4))-paa(1)*(pim3(1)-pim1(1))
3595  $ -paa(2)*(pim3(2)-pim1(2))-paa(3)*(pim3(3)-pim1(3))
3596  prod3 =paa(4)*(pim1(4)-pim2(4))-paa(1)*(pim1(1)-pim2(1))
3597  $ -paa(2)*(pim1(2)-pim2(2))-paa(3)*(pim1(3)-pim2(3))
3598  DO 40 i=1,4
3599  vec1(i)= pim2(i)-pim3(i) -paa(i)*prod1/xmaa**2
3600  vec2(i)= pim3(i)-pim1(i) -paa(i)*prod2/xmaa**2
3601  vec3(i)= pim1(i)-pim2(i) -paa(i)*prod3/xmaa**2
3602  40 vec4(i)= pim1(i)+pim2(i)+pim3(i)
3603  CALL prod5(pim1,pim2,pim3,vec5)
3604 * HADRON CURRENT
3605 C be aware that sign of vec2 is opposite to sign of vec1 in a1 case
3606 
3607 C Rationalize this code:
3608  f1 = cmplx(coef(1,mnum))*form1(mnum,xmaa**2,xmro1**2,xmro2**2)
3609  f2 = cmplx(coef(2,mnum))*form2(mnum,xmaa**2,xmro2**2,xmro1**2)
3610  f3 = cmplx(coef(3,mnum))*form3(mnum,xmaa**2,xmro3**2,xmro1**2)
3611  f4 = (-1.0*uroj)*
3612  $cmplx(coef(4,mnum))*form4(mnum,xmaa**2,xmro1**2,xmro2**2,xmro3**2)
3613  f5 = (-1.0)*uroj/4.0/pi**2/fpi**2*
3614  $ cmplx(coef(5,mnum))*form5(mnum,xmaa**2,xmro1**2,xmro2**2)
3615 ! if (mnum.eq.9) write(*,*) 'effy=', mnum,'>>',f1,f2,f3,f4,f5
3616  DO 45 i=1,4
3617  hadcur(i)= cmplx(fnorm(mnum)) * (
3618  $ cmplx(vec1(i))*f1+cmplx(vec2(i))*f2+cmplx(vec3(i))*f3+
3619  $ cmplx(vec4(i))*f4+cmplx(vec5(i))*f5)
3620  45 CONTINUE
3621 
3622 C
3623 * CALCULATE PI-VECTORS: VECTOR AND AXIAL
3624  CALL clvec(hadcur,pn,pivec)
3625  CALL claxi(hadcur,pn,piaks)
3626  CALL clnut(hadcur,brakm,hvm)
3627 * SPIN INDEPENDENT PART OF DECAY DIFF-CROSS-SECT. IN TAU REST FRAME
3628  brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3629  & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3630  amplit=(gfermi)**2*brak/2.
3631  IF (mnum.GE.20) THEN
3632  print *, 'MNUM=',mnum
3633  znak=-1.0
3634  xm1=0.0
3635  xm2=0.0
3636  xm3=0.0
3637  DO 77 k=1,4
3638  IF (k.EQ.4) znak=1.0
3639  xm1=znak*pim1(k)**2+xm1
3640  xm2=znak*pim2(k)**2+xm2
3641  xm3=znak*pim3(k)**2+xm3
3642  77 print *, 'PIM1=',pim1(k),'PIM2=',pim2(k),'PIM3=',pim3(k)
3643  print *, 'XM1=',sqrt(xm1),'XM2=',sqrt(xm2),'XM3=',sqrt(xm3)
3644  print *, '************************************************'
3645  ENDIF
3646 C POLARIMETER VECTOR IN TAU REST FRAME
3647  DO 90 i=1,3
3648  hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3649  & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3650 C HV IS DEFINED FOR TAU- WITH GAMMA=B+HV*POL
3651  hv(i)=-hv(i)/brak
3652  90 CONTINUE
3653  END
3654  SUBROUTINE prod5(P1,P2,P3,PIA)
3655 C ----------------------------------------------------------------------
3656 C external product of P1, P2, P3 4-momenta.
3657 C SIGN is chosen +/- for decay of TAU +/- respectively
3658 C called by : DAMPAA, CLNUT
3659 C ----------------------------------------------------------------------
3660  COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3661  COMMON / idfc / idff
3662  REAL PIA(4),P1(4),P2(4),P3(4)
3663  det2(i,j)=p1(i)*p2(j)-p2(i)*p1(j)
3664 * -----------------------------------
3665  IF (ktom.EQ.1.OR.ktom.EQ.-1) THEN
3666  sign= idff/abs(idff)
3667  ELSEIF (ktom.EQ.2) THEN
3668  sign=-idff/abs(idff)
3669  ELSE
3670  print *, 'STOP IN PROD5: KTOM=',ktom
3671  stop
3672  ENDIF
3673 C
3674 C EPSILON( p1(1), p2(2), p3(3), (4) ) = 1
3675 C
3676  pia(1)= -p3(3)*det2(2,4)+p3(4)*det2(2,3)+p3(2)*det2(3,4)
3677  pia(2)= -p3(4)*det2(1,3)+p3(3)*det2(1,4)-p3(1)*det2(3,4)
3678  pia(3)= p3(4)*det2(1,2)-p3(2)*det2(1,4)+p3(1)*det2(2,4)
3679  pia(4)= p3(3)*det2(1,2)-p3(2)*det2(1,3)+p3(1)*det2(2,3)
3680 C ALL FOUR INDICES ARE UP SO PIA(3) AND PIA(4) HAVE SAME SIGN
3681  DO 20 i=1,4
3682  20 pia(i)=pia(i)*sign
3683  END
3684 
3685  SUBROUTINE dexnew(MODE,ISGN,POL,PNU,PAA,PNPI,JNPI)
3686 C ----------------------------------------------------------------------
3687 * THIS SIMULATES TAU DECAY IN TAU REST FRAME
3688 * INTO NU A1, NEXT A1 DECAYS INTO RHO PI AND FINALLY RHO INTO PI PI.
3689 * OUTPUT FOUR MOMENTA: PNU TAUNEUTRINO,
3690 
3691 * PAA A1
3692 * PIM1 PION MINUS (OR PI0) 1 (FOR TAU MINUS)
3693 * PIM2 PION MINUS (OR PI0) 2
3694 * PIPL PION PLUS (OR PI-)
3695 * (PIPL,PIM1) FORM A RHO
3696 
3697 C ----------------------------------------------------------------------
3698  COMMON / inout / inut,iout
3699  REAL POL(4),HV(4),PAA(4),PNU(4),PNPI(4,9),RN(1)
3700  DATA iwarm/0/
3701 C
3702  IF(mode.EQ.-1) THEN
3703 C ===================
3704  iwarm=1
3705  CALL dadnew( -1,isgn,hv,pnu,paa,pnpi,jdumm)
3706 
3707 CC CALL HBOOK1(816,'WEIGHT DISTRIBUTION DEXAA $',100,-2.,2.)
3708 
3709 C
3710  ELSEIF(mode.EQ. 0) THEN
3711 * =======================
3712  300 CONTINUE
3713  IF(iwarm.EQ.0) GOTO 902
3714  CALL dadnew( 0,isgn,hv,pnu,paa,pnpi,jnpi)
3715  wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
3716 CC CALL HFILL(816,WT)
3717  CALL ranmar(rn,1)
3718  IF(rn(1).GT.wt) GOTO 300
3719 C
3720  ELSEIF(mode.EQ. 1) THEN
3721 * =======================
3722  CALL dadnew( 1,isgn,hv,pnu,paa,pnpi,jdumm)
3723 CC CALL HPRINT(816)
3724  ENDIF
3725 C =====
3726  RETURN
3727  902 WRITE(iout, 9020)
3728  9020 FORMAT(' ----- DEXNEW: LACK OF INITIALISATION')
3729  stop
3730  END
3731  SUBROUTINE dadnew(MODE,ISGN,HV,PNU,PWB,PNPI,JNPI)
3732 C ----------------------------------------------------------------------
3733  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3734  * ,ampiz,ampi,amro,gamro,ama1,gama1
3735  * ,amk,amkz,amkst,gamkst
3736 C
3737  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
3738  * ,ampiz,ampi,amro,gamro,ama1,gama1
3739  * ,amk,amkz,amkst,gamkst
3740  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3741  REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
3742  COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
3743  REAL*4 GAMPMC ,GAMPER
3744 
3745  COMMON / inout / inut,iout
3746 
3747  parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
3748 
3749  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3750 
3751  & ,names
3752  CHARACTER NAMES(NMODE)*31
3753 
3754 
3755  REAL*4 PNU(4),PWB(4),PNPI(4,9),HV(4),HHV(4)
3756  REAL*4 PDUM1(4),PDUM2(4),PDUMI(4,9)
3757  REAL*4 RRR(3)
3758  REAL*4 WTMAX(NMODE)
3759  REAL*8 SWT(NMODE),SSWT(NMODE)
3760  dimension nevraw(nmode),nevovr(nmode),nevacc(nmode)
3761 C
3762  DATA pi /3.141592653589793238462643/
3763  DATA iwarm/0/
3764 C
3765  IF(mode.EQ.-1) THEN
3766 C ===================
3767 C -- AT THE MOMENT ONLY TWO DECAY MODES OF MULTIPIONS HAVE M. ELEM
3768  nmod=nmode
3769  iwarm=1
3770 C PRINT 7003
3771  DO 1 jnpi=1,nmod
3772  nevraw(jnpi)=0
3773  nevacc(jnpi)=0
3774  nevovr(jnpi)=0
3775  swt(jnpi)=0
3776  sswt(jnpi)=0
3777  wtmax(jnpi)=-1.
3778 
3779 C for 4pi phase space, need lots more trials at initialization,
3780 C or use the WTMAX determined with many trials for default model:
3781  ntrials = 45000
3782  IF (jnpi.LE.2) THEN
3783  a = pkorb(3,37+jnpi)
3784  IF (a.LT.0.) THEN
3785  ntrials = 30000
3786  ELSE
3787  ntrials = 0
3788  wtmax(jnpi)=a
3789  END IF
3790  END IF
3791  ! write(*,*) 'jnpi=',jnpi, ntrials, a
3792  DO i=1,ntrials
3793 
3794  IF (jnpi.LE.0) THEN
3795  GOTO 903
3796  ELSEIF(jnpi.LE.nm4) THEN
3797  CALL dph4pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3798 ! IF (I.eq.1) write(*,*) '4 pi jnpi=',jnpi
3799  ELSEIF(jnpi.LE.nm4+nm5) THEN
3800  CALL dph5pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3801  ELSEIF(jnpi.LE.nm4+nm5+nm6) THEN
3802  CALL dphnpi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3803  ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3) THEN
3804  inum=jnpi-nm4-nm5-nm6
3805  CALL dphspk(wt,hv,pdum1,pdum2,pdumi,inum)
3806  ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2) THEN
3807  inum=jnpi-nm4-nm5-nm6-nm3
3808  CALL dphsrk(wt,hv,pdum1,pdum2,pdumi,inum)
3809  ELSE
3810  GOTO 903
3811  ENDIF
3812  IF(wt.GT.wtmax(jnpi)/1.2) wtmax(jnpi)=wt*1.2
3813  ENDDO
3814 
3815 C PRINT *,' DADNEW JNPI,NTRIALS,WTMAX =',JNPI,NTRIALS,WTMAX(JNPI)
3816 
3817 C CALL HBOOK1(801,'WEIGHT DISTRIBUTION DADNPI $',100,0.,2.,.0)
3818 C PRINT 7004,WTMAX(JNPI)
3819 1 CONTINUE
3820  WRITE(iout,7005)
3821 C
3822  ELSEIF(mode.EQ. 0) THEN
3823 C =======================
3824  IF(iwarm.EQ.0) GOTO 902
3825 C
3826 300 CONTINUE
3827  IF (jnpi.LE.0) THEN
3828  GOTO 903
3829  ELSEIF(jnpi.LE.nm4) THEN
3830  CALL dph4pi(wt,hhv,pnu,pwb,pnpi,jnpi)
3831  ELSEIF(jnpi.LE.nm4+nm5) THEN
3832  CALL dph5pi(wt,hhv,pnu,pwb,pnpi,jnpi)
3833  ELSEIF(jnpi.LE.nm4+nm5+nm6) THEN
3834  CALL dphnpi(wt,hhv,pnu,pwb,pnpi,jnpi)
3835  ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3) THEN
3836  inum=jnpi-nm4-nm5-nm6
3837  CALL dphspk(wt,hhv,pnu,pwb,pnpi,inum)
3838  ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2) THEN
3839  inum=jnpi-nm4-nm5-nm6-nm3
3840  CALL dphsrk(wt,hhv,pnu,pwb,pnpi,inum)
3841  ELSE
3842  GOTO 903
3843  ENDIF
3844  DO i=1,4
3845  hv(i)=-isgn*hhv(i)
3846  ENDDO
3847 C CALL HFILL(801,WT/WTMAX(JNPI))
3848  nevraw(jnpi)=nevraw(jnpi)+1
3849  swt(jnpi)=swt(jnpi)+wt
3850 
3851 cccM.S.>>>>>>
3852 cc SSWT(JNPI)=SSWT(JNPI)+WT**2
3853  sswt(jnpi)=sswt(jnpi)+dble(wt)**2
3854 cccM.S.<<<<<<
3855 
3856  CALL ranmar(rrr,3)
3857  rn=rrr(1)
3858  IF(wt.GT.wtmax(jnpi)) nevovr(jnpi)=nevovr(jnpi)+1
3859  IF(rn*wtmax(jnpi).GT.wt) GOTO 300
3860 C ROTATIONS TO BASIC TAU REST FRAME
3861  costhe=-1.+2.*rrr(2)
3862  thet=acos(costhe)
3863  phi =2*pi*rrr(3)
3864  CALL rotor2(thet,pnu,pnu)
3865  CALL rotor3( phi,pnu,pnu)
3866  CALL rotor2(thet,pwb,pwb)
3867  CALL rotor3( phi,pwb,pwb)
3868  CALL rotor2(thet,hv,hv)
3869  CALL rotor3( phi,hv,hv)
3870  nd=mulpik(jnpi)
3871  DO 301 i=1,nd
3872  CALL rotor2(thet,pnpi(1,i),pnpi(1,i))
3873  CALL rotor3( phi,pnpi(1,i),pnpi(1,i))
3874 301 CONTINUE
3875  nevacc(jnpi)=nevacc(jnpi)+1
3876 C
3877  ELSEIF(mode.EQ. 1) THEN
3878 C =======================
3879  DO 500 jnpi=1,nmod
3880  IF(nevraw(jnpi).EQ.0) GOTO 500
3881  pargam=swt(jnpi)/float(nevraw(jnpi)+1)
3882  error=0
3883  IF(nevraw(jnpi).NE.0)
3884  & error=sqrt(sswt(jnpi)/swt(jnpi)**2-1./float(nevraw(jnpi)))
3885  rat=pargam/gamel
3886  WRITE(iout, 7010) names(jnpi),
3887  & nevraw(jnpi),nevacc(jnpi),nevovr(jnpi),pargam,rat,error
3888 CC CALL HPRINT(801)
3889  gampmc(8+jnpi-1)=rat
3890  gamper(8+jnpi-1)=error
3891 CAM NEVDEC(8+JNPI-1)=NEVACC(JNPI)
3892  500 CONTINUE
3893  ENDIF
3894 C =====
3895  RETURN
3896  7003 FORMAT(///1x,15(5h*****)
3897  $ /,' *', 25x,'******** DADNEW INITIALISATION ********',9x,1h*
3898  $ )
3899  7004 FORMAT(' *',e20.5,5x,'WTMAX = MAXIMUM WEIGHT ',9x,1h*/)
3900  7005 FORMAT(
3901  $ /,1x,15(5h*****)/)
3902  7010 FORMAT(///1x,15(5h*****)
3903  $ /,' *', 25x,'******** DADNEW FINAL REPORT ******** ',9x,1h*
3904  $ /,' *', 25x,'CHANNEL:',a31 ,9x,1h*
3905  $ /,' *',i20 ,5x,'NEVRAW = NO. OF DECAYS TOTAL ',9x,1h*
3906  $ /,' *',i20 ,5x,'NEVACC = NO. OF DECAYS ACCEPTED ',9x,1h*
3907  $ /,' *',i20 ,5x,'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
3908  $ /,' *',e20.5,5x,'PARTIAL WTDTH IN GEV UNITS ',9x,1h*
3909  $ /,' *',f20.9,5x,'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
3910  $ /,' *',f20.8,5x,'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
3911  $ /,1x,15(5h*****)/)
3912  902 WRITE(iout, 9020)
3913  9020 FORMAT(' ----- DADNEW: LACK OF INITIALISATION')
3914  stop
3915  903 WRITE(iout, 9030) jnpi,mode
3916  9030 FORMAT(' ----- DADNEW: WRONG JNPI',2i5)
3917  stop
3918  END
3919 
3920 
3921  SUBROUTINE dph4pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
3922 C ----------------------------------------------------------------------
3923 
3924 * IT SIMULATES A1 DECAY IN TAU REST FRAME WITH
3925 * Z-AXIS ALONG A1 MOMENTUM
3926 
3927 C ----------------------------------------------------------------------
3928  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3929  * ,ampiz,ampi,amro,gamro,ama1,gama1
3930  * ,amk,amkz,amkst,gamkst
3931 C
3932  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
3933  * ,ampiz,ampi,amro,gamro,ama1,gama1
3934  * ,amk,amkz,amkst,gamkst
3935  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3936  REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
3937 
3938 
3939  REAL HV(4),PT(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4),PMULT(4,9)
3940  REAL PR(4),PIZ(4)
3941  REAL*4 RRR(9)
3942  REAL*8 UU,FF,FF1,FF2,FF3,FF4,GG1,GG2,GG3,GG4,RR
3943  DATA pi /3.141592653589793238462643/
3944  DATA icont /0/
3945  xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
3946 C AMRO, GAMRO IS ONLY A PARAMETER FOR GETING HIGHT EFFICIENCY
3947 C
3948 C THREE BODY PHASE SPACE NORMALISED AS IN BJORKEN-DRELL
3949 C D**3 P /2E/(2PI)**3 (2PI)**4 DELTA4(SUM P)
3950  phspac=1./2**23/pi**11
3951  phsp=1./2**5/pi**2
3952 
3953  IF (jnpi.EQ.1) THEN
3954  prez=0.7
3955 
3956  amp1=ampi
3957  amp2=ampi
3958  amp3=ampi
3959  amp4=ampiz
3960  amrx=pkorb(1,14)
3961  gamrx=pkorb(2,14)
3962 C AJW: cant simply change AMROP, etc, here!
3963 C CHOICE is a by-hand tuning/optimization, no simple relationship
3964 C to actual resonance masses (accd to Z.Was).
3965 C What matters in the end is what you put in formf/curr .
3966 
3967  amrop =1.2
3968  gamrop=.46
3969  ELSE
3970  prez=0.0
3971 
3972  amp1=ampiz
3973  amp2=ampiz
3974  amp3=ampiz
3975  amp4=ampi
3976 
3977  amrx=1.4
3978  gamrx=.6
3979  amrop =amrx
3980  gamrop=gamrx
3981 
3982  ENDIF
3983 
3984  rrb=0.3
3985  CALL choice(100+jnpi,rrb,ichan,prob1,prob2,prob3,
3986 
3987  $ amrop,gamrop,amrx,gamrx,amrb,gamrb)
3988  prez=prob1+prob2
3989 C TAU MOMENTUM
3990  pt(1)=0.
3991  pt(2)=0.
3992  pt(3)=0.
3993  pt(4)=amtau
3994 C
3995  CALL ranmar(rrr,9)
3996 C
3997 * MASSES OF 4, 3 AND 2 PI SYSTEMS
3998 C 3 PI WITH SAMPLING FOR RESONANCE
3999 CAM
4000  rr1=rrr(6)
4001  ams1=(amp1+amp2+amp3+amp4)**2
4002  ams2=(amtau-amnuta)**2
4003  alp1=atan((ams1-amrop**2)/amrop/gamrop)
4004  alp2=atan((ams2-amrop**2)/amrop/gamrop)
4005  alp=alp1+rr1*(alp2-alp1)
4006  am4sq =amrop**2+amrop*gamrop*tan(alp)
4007  am4 =sqrt(am4sq)
4008  phspac=phspac*
4009  $ ((am4sq-amrop**2)**2+(amrop*gamrop)**2)/(amrop*gamrop)
4010  phspac=phspac*(alp2-alp1)
4011 
4012 C
4013  rr1=rrr(1)
4014  ams1=(amp2+amp3+amp4)**2
4015  ams2=(am4-amp1)**2
4016  IF (rrr(9).GT.prez) THEN
4017  am3sq=ams1+ rr1*(ams2-ams1)
4018  am3 =sqrt(am3sq)
4019 C --- this part of jacobian will be recovered later
4020  ff1=ams2-ams1
4021  ELSE
4022 * PHASE SPACE WITH SAMPLING FOR OMEGA RESONANCE,
4023  alp1=atan((ams1-amrx**2)/amrx/gamrx)
4024  alp2=atan((ams2-amrx**2)/amrx/gamrx)
4025  alp=alp1+rr1*(alp2-alp1)
4026  am3sq =amrx**2+amrx*gamrx*tan(alp)
4027  am3 =sqrt(am3sq)
4028 C --- THIS PART OF THE JACOBIAN WILL BE RECOVERED LATER ---------------
4029  ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4030  ff1=ff1*(alp2-alp1)
4031  ENDIF
4032 C MASS OF 2
4033  rr2=rrr(2)
4034  ams1=(amp3+amp4)**2
4035  ams2=(am3-amp2)**2
4036 * FLAT PHASE SPACE;
4037  am2sq=ams1+ rr2*(ams2-ams1)
4038  am2 =sqrt(am2sq)
4039 C --- this part of jacobian will be recovered later
4040  ff2=(ams2-ams1)
4041 * 2 RESTFRAME, DEFINE PIZ AND PIPL
4042 
4043  enq1=(am2sq-amp3**2+amp4**2)/(2*am2)
4044  enq2=(am2sq+amp3**2-amp4**2)/(2*am2)
4045  ppi= enq1**2-amp4**2
4046  pppi=sqrt(abs(enq1**2-amp4**2))
4047 
4048  phspac=phspac*(4*pi)*(2*pppi/am2)
4049 
4050 * PIZ MOMENTUM IN 2 REST FRAME
4051 
4052  CALL sphera(pppi,piz)
4053  piz(4)=enq1
4054 
4055 * PIPL MOMENTUM IN 2 REST FRAME
4056 
4057  DO 30 i=1,3
4058  30 pipl(i)=-piz(i)
4059  pipl(4)=enq2
4060 * 3 REST FRAME, DEFINE PIM1
4061 
4062 * PR MOMENTUM
4063 
4064  pr(1)=0
4065  pr(2)=0
4066  pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4067  pr(3)= sqrt(abs(pr(4)**2-am2**2))
4068  ppi = pr(4)**2-am2**2
4069 
4070 * PIM1 MOMENTUM
4071 
4072  pim1(1)=0
4073  pim1(2)=0
4074  pim1(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4075  pim1(3)=-pr(3)
4076 C --- this part of jacobian will be recovered later
4077  ff3=(4*pi)*(2*pr(3)/am3)
4078 * OLD PIONS BOOSTED FROM 2 REST FRAME TO 3 REST FRAME
4079  exe=(pr(4)+pr(3))/am2
4080  CALL bostr3(exe,piz,piz)
4081  CALL bostr3(exe,pipl,pipl)
4082  rr3=rrr(3)
4083  rr4=rrr(4)
4084  thet =acos(-1.+2*rr3)
4085  phi = 2*pi*rr4
4086  CALL rotpol(thet,phi,pipl)
4087  CALL rotpol(thet,phi,pim1)
4088  CALL rotpol(thet,phi,piz)
4089  CALL rotpol(thet,phi,pr)
4090 
4091 * 4 REST FRAME, DEFINE PIM2
4092 * PR MOMENTUM
4093 
4094  pr(1)=0
4095  pr(2)=0
4096  pr(4)=1./(2*am4)*(am4**2+am3**2-amp1**2)
4097  pr(3)= sqrt(abs(pr(4)**2-am3**2))
4098  ppi = pr(4)**2-am3**2
4099 
4100 * PIM2 MOMENTUM
4101 
4102  pim2(1)=0
4103  pim2(2)=0
4104  pim2(4)=1./(2*am4)*(am4**2-am3**2+amp1**2)
4105  pim2(3)=-pr(3)
4106 C --- this part of jacobian will be recovered later
4107  ff4=(4*pi)*(2*pr(3)/am4)
4108 * OLD PIONS BOOSTED FROM 3 REST FRAME TO 4 REST FRAME
4109  exe=(pr(4)+pr(3))/am3
4110  CALL bostr3(exe,piz,piz)
4111  CALL bostr3(exe,pipl,pipl)
4112  CALL bostr3(exe,pim1,pim1)
4113  rr3=rrr(7)
4114  rr4=rrr(8)
4115  thet =acos(-1.+2*rr3)
4116  phi = 2*pi*rr4
4117  CALL rotpol(thet,phi,pipl)
4118  CALL rotpol(thet,phi,pim1)
4119  CALL rotpol(thet,phi,pim2)
4120  CALL rotpol(thet,phi,piz)
4121  CALL rotpol(thet,phi,pr)
4122 C
4123 * NOW TO THE TAU REST FRAME, DEFINE PAA AND NEUTRINO MOMENTA
4124 * PAA MOMENTUM
4125  paa(1)=0
4126  paa(2)=0
4127  paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am4**2)
4128  paa(3)= sqrt(abs(paa(4)**2-am4**2))
4129  ppi = paa(4)**2-am4**2
4130  phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4131  phsp=phsp*(4*pi)*(2*paa(3)/amtau)
4132 * TAU-NEUTRINO MOMENTUM
4133  pn(1)=0
4134  pn(2)=0
4135  pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am4**2)
4136  pn(3)=-paa(3)
4137 C ZBW 20.12.2002 bug fix
4138  IF(rrr(9).LE.0.5*prez) THEN
4139  DO 72 i=1,4
4140  x=pim1(i)
4141  pim1(i)=pim2(i)
4142  72 pim2(i)=x
4143  ENDIF
4144 C end of bug fix
4145 C WE INCLUDE REMAINING PART OF THE JACOBIAN
4146 C --- FLAT CHANNEL
4147  am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4148  $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4149  ams2=(am4-amp2)**2
4150  ams1=(amp1+amp3+amp4)**2
4151  ff1=(ams2-ams1)
4152  ams1=(amp3+amp4)**2
4153  ams2=(sqrt(am3sq)-amp1)**2
4154  ff2=ams2-ams1
4155  ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4156  ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4157  uu=ff1*ff2*ff3*ff4
4158 C --- FIRST CHANNEL
4159  am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4160  $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4161  ams2=(am4-amp2)**2
4162  ams1=(amp1+amp3+amp4)**2
4163  alp1=atan((ams1-amrx**2)/amrx/gamrx)
4164  alp2=atan((ams2-amrx**2)/amrx/gamrx)
4165  ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4166  ff1=ff1*(alp2-alp1)
4167  ams1=(amp3+amp4)**2
4168  ams2=(sqrt(am3sq)-amp1)**2
4169  ff2=ams2-ams1
4170  ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4171  ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4172  ff=ff1*ff2*ff3*ff4
4173 C --- SECOND CHANNEL
4174  am3sq=(pim2(4)+piz(4)+pipl(4))**2-(pim2(3)+piz(3)+pipl(3))**2
4175  $ -(pim2(2)+piz(2)+pipl(2))**2-(pim2(1)+piz(1)+pipl(1))**2
4176  ams2=(am4-amp1)**2
4177  ams1=(amp2+amp3+amp4)**2
4178  alp1=atan((ams1-amrx**2)/amrx/gamrx)
4179  alp2=atan((ams2-amrx**2)/amrx/gamrx)
4180  gg1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4181  gg1=gg1*(alp2-alp1)
4182  ams1=(amp3+amp4)**2
4183  ams2=(sqrt(am3sq)-amp2)**2
4184  gg2=ams2-ams1
4185  gg3=(4*pi)*(xlam(am2**2,amp2**2,am3sq)/am3sq)
4186  gg4=(4*pi)*(xlam(am3sq,amp1**2,am4**2)/am4**2)
4187  gg=gg1*gg2*gg3*gg4
4188 C --- JACOBIAN AVERAGED OVER THE TWO
4189  IF ( ( (ff+gg)*uu+ff*gg ).GT.0.0d0) THEN
4190  rr=ff*gg*uu/(0.5*prez*(ff+gg)*uu+(1.0-prez)*ff*gg)
4191  phspac=phspac*rr
4192  ELSE
4193  phspac=0.0
4194  ENDIF
4195 * MOMENTA OF THE TWO PI-MINUS ARE RANDOMLY SYMMETRISED
4196  IF (jnpi.EQ.1) THEN
4197  rr5= rrr(5)
4198  IF(rr5.LE.0.5) THEN
4199  DO 70 i=1,4
4200  x=pim1(i)
4201  pim1(i)=pim2(i)
4202  70 pim2(i)=x
4203  ENDIF
4204  phspac=phspac/2.
4205  ELSE
4206 C MOMENTA OF PI0-S ARE GENERATED UNIFORMLY ONLY IF PREZ=0.0
4207  rr5= rrr(5)
4208  IF(rr5.LE.0.5) THEN
4209  DO 71 i=1,4
4210  x=pim1(i)
4211  pim1(i)=pim2(i)
4212  71 pim2(i)=x
4213  ENDIF
4214  phspac=phspac/6.
4215  ENDIF
4216 * ALL PIONS BOOSTED FROM 4 REST FRAME TO TAU REST FRAME
4217 * Z-AXIS ANTIPARALLEL TO NEUTRINO MOMENTUM
4218  exe=(paa(4)+paa(3))/am4
4219  CALL bostr3(exe,piz,piz)
4220  CALL bostr3(exe,pipl,pipl)
4221  CALL bostr3(exe,pim1,pim1)
4222  CALL bostr3(exe,pim2,pim2)
4223  CALL bostr3(exe,pr,pr)
4224 C PARTIAL WIDTH CONSISTS OF PHASE SPACE AND AMPLITUDE
4225 C CHECK ON CONSISTENCY WITH DADNPI, THEN, CODE BREAKES UNIFORM PION
4226 C DISTRIBUTION IN HADRONIC SYSTEM
4227 
4228 CAM Assume neutrino mass=0. and sum over final polarisation
4229 C AMX2=AM4**2
4230 C BRAK= 2*(AMTAU**2-AMX2) * (AMTAU**2+2.*AMX2)
4231 C AMPLIT=CCABIB**2*GFERMI**2/2. * BRAK * AMX2*SIGEE(AMX2,1)
4232  IF (jnpi.EQ.1) THEN
4233  CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4234  ELSEIF (jnpi.EQ.2) THEN
4235  CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv)
4236  ELSE
4237  CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv) ! temporarily
4238  ENDIF
4239 
4240  dgamt=1/(2.*amtau)*amplit*phspac
4241 C PHASE SPACE CHECK
4242 C DGAMT=PHSPAC
4243  DO 77 k=1,4
4244  pmult(k,1)=pim1(k)
4245  pmult(k,2)=pim2(k)
4246 
4247  pmult(k,3)=pipl(k)
4248  pmult(k,4)=piz(k)
4249 
4250  77 CONTINUE
4251  END
4252  SUBROUTINE dam4pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,AMPLIT,HV)
4253 C ----------------------------------------------------------------------
4254 * CALCULATES DIFFERENTIAL CROSS SECTION AND POLARIMETER VECTOR
4255 * FOR TAU DECAY INTO 4 PI MODES
4256 * ALL SPIN EFFECTS IN THE FULL DECAY CHAIN ARE TAKEN INTO ACCOUNT.
4257 * CALCULATIONS DONE IN TAU REST FRAME WITH Z-AXIS ALONG NEUTRINO MOMENT
4258 C MNUM DECAY MODE IDENTIFIER.
4259 C
4260 
4261 C called by : DPHSAA
4262 
4263 C ----------------------------------------------------------------------
4264  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4265  * ,ampiz,ampi,amro,gamro,ama1,gama1
4266  * ,amk,amkz,amkst,gamkst
4267 C
4268  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
4269  * ,ampiz,ampi,amro,gamro,ama1,gama1
4270  * ,amk,amkz,amkst,gamkst
4271  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4272  REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
4273  REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4),PIM4(4)
4274  REAL PIVEC(4),PIAKS(4),HVM(4)
4275  COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5
4276  EXTERNAL form1,form2,form3,form4,form5
4277  DATA pi /3.141592653589793238462643/
4278  DATA icont /0/
4279 C
4280 ! write(*,*) 'falanti',mnum
4281  CALL curr_cleo(mnum,pim1,pim2,pim3,pim4,hadcur)
4282 
4283 C
4284 * CALCULATE PI-VECTORS: VECTOR AND AXIAL
4285  CALL clvec(hadcur,pn,pivec)
4286  CALL claxi(hadcur,pn,piaks)
4287  CALL clnut(hadcur,brakm,hvm)
4288 * SPIN INDEPENDENT PART OF DECAY DIFF-CROSS-SECT. IN TAU REST FRAME
4289  brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4290  & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4291  amplit=(ccabib*gfermi)**2*brak/2.
4292 C POLARIMETER VECTOR IN TAU REST FRAME
4293  DO 90 i=1,3
4294  hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4295  & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4296 C HV IS DEFINED FOR TAU- WITH GAMMA=B+HV*POL
4297  IF (brak.NE.0.0)
4298  &hv(i)=-hv(i)/brak
4299  90 CONTINUE
4300  END
4301  SUBROUTINE dam5pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,PIM5,AMPLIT,HV)
4302 C ----------------------------------------------------------------------
4303 * CALCULATES DIFFERENTIAL CROSS SECTION AND POLARIMETER VECTOR
4304 * FOR TAU DECAY INTO 4 PI MODES
4305 * ALL SPIN EFFECTS IN THE FULL DECAY CHAIN ARE TAKEN INTO ACCOUNT.
4306 * CALCULATIONS DONE IN TAU REST FRAME WITH Z-AXIS ALONG NEUTRINO MOMENT
4307 C MNUM DECAY MODE IDENTIFIER.
4308 C
4309 
4310 C called by : DPHSAA
4311 
4312 C ----------------------------------------------------------------------
4313  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4314  * ,ampiz,ampi,amro,gamro,ama1,gama1
4315  * ,amk,amkz,amkst,gamkst
4316 C
4317  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
4318  * ,ampiz,ampi,amro,gamro,ama1,gama1
4319  * ,amk,amkz,amkst,gamkst
4320  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4321  REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
4322  REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4),PIM4(4),PIM5(4)
4323  REAL PIVEC(4),PIAKS(4),HVM(4)
4324  COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5
4325  EXTERNAL form1,form2,form3,form4,form5
4326  DATA pi /3.141592653589793238462643/
4327  DATA icont /0/
4328 C
4329 ! write(*,*) 'falanti',mnum
4330  CALL curr5(mnum,pim1,pim2,pim3,pim4,pim5,hadcur)
4331 
4332 C
4333 * CALCULATE PI-VECTORS: VECTOR AND AXIAL
4334  CALL clvec(hadcur,pn,pivec)
4335  CALL claxi(hadcur,pn,piaks)
4336  CALL clnut(hadcur,brakm,hvm)
4337 * SPIN INDEPENDENT PART OF DECAY DIFF-CROSS-SECT. IN TAU REST FRAME
4338  brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4339  & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4340  amplit=(ccabib*gfermi)**2*brak/2.
4341 C POLARIMETER VECTOR IN TAU REST FRAME
4342  DO 90 i=1,3
4343  hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4344  & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4345 C HV IS DEFINED FOR TAU- WITH GAMMA=B+HV*POL
4346  IF (brak.NE.0.0)
4347  &hv(i)=-hv(i)/brak
4348  90 CONTINUE
4349  END
4350  SUBROUTINE dph5pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4351 C ----------------------------------------------------------------------
4352 * IT SIMULATES 5pi DECAY IN TAU REST FRAME WITH
4353 * Z-AXIS ALONG 5pi MOMENTUM
4354 C ----------------------------------------------------------------------
4355  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4356  * ,ampiz,ampi,amro,gamro,ama1,gama1
4357  * ,amk,amkz,amkst,gamkst
4358 C
4359  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
4360  * ,ampiz,ampi,amro,gamro,ama1,gama1
4361 
4362 
4363  * ,amk,amkz,amkst,gamkst
4364  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4365  REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
4366  parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
4367 
4368  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4369 
4370  & ,names
4371  CHARACTER NAMES(NMODE)*31
4372  REAL HV(4),PT(4),PN(4),PAA(4),PMULT(4,9)
4373  REAL*4 PR(4),PI1(4),PI2(4),PI3(4),PI4(4),PI5(4)
4374  REAL*8 AMP1,AMP2,AMP3,AMP4,AMP5,ams1,ams2,amom,gamom
4375  REAL*8 AM5SQ,AM4SQ,AM3SQ,AM2SQ,AM5,AM4,AM3
4376  REAL*4 RRR(12)
4377  REAL*8 gg1,gg2,gg3,ff1,ff2,ff3,ff4,alp,alp1,alp2
4378 
4379  REAL*8 XM,AM,GAMMA
4380 ccM.S.>>>>>>
4381  real*8 phspac
4382 ccM.S.<<<<<<
4383 
4384  DATA pi /3.141592653589793238462643/
4385  DATA icont /0/
4386  data fpi /93.3e-3/
4387 c
4388  COMPLEX BWIGN
4389 C
4390 
4391  bwign(xm,am,gamma)=xm**2/cmplx(xm**2-am**2,gamma*am)
4392 
4393 
4394  proba2=0.7
4395  probom=0.7
4396  ama2=1.260
4397  gama2=0.400
4398 
4399  amom=.782
4400  gamom=.7
4401 
4402  IF (jnpi.EQ.23.OR.jnpi.EQ.24) gamom= 0.0085 ! it is too early for
4403  ! centralized control of presampling
4404 c
4405 C 6 BODY PHASE SPACE NORMALISED AS IN BJORKEN-DRELL
4406 C D**3 P /2E/(2PI)**3 (2PI)**4 DELTA4(SUM P)
4407  phspac=1./2**29/pi**14
4408 c PHSPAC=1./2**5/PI**2
4409 C init 5pi decay mode (JNPI)
4410  amp1=dcdmas(idffin(1,jnpi))
4411  amp2=dcdmas(idffin(2,jnpi))
4412  amp3=dcdmas(idffin(3,jnpi))
4413  amp4=dcdmas(idffin(4,jnpi))
4414  amp5=dcdmas(idffin(5,jnpi))
4415 c
4416 C TAU MOMENTUM
4417  pt(1)=0.
4418  pt(2)=0.
4419  pt(3)=0.
4420  pt(4)=amtau
4421 C
4422  CALL ranmar(rrr,12)
4423 C
4424 c masses of 5, 4, 3 and 2 pi systems
4425 c 3 pi with sampling for omega resonance
4426 cam
4427 c mass of 5 (12345)
4428  IF (rrr(11).GT.proba2) THEN
4429 c flat phase space:
4430  rr1=rrr(10)
4431  ams1=(amp1+amp2+amp3+amp4+amp5)**2
4432  ams2=(amtau-amnuta)**2
4433  alp1=atan((ams1-ama2**2)/ama2/gama2)
4434  alp2=atan((ams2-ama2**2)/ama2/gama2)
4435  am5sq=ams1+ rr1*(ams2-ams1)
4436  am5 =sqrt(am5sq)
4437 C phspac=phspac*(ams2-ams1)
4438 c or peaked phase space for a1(?) resonance:
4439  ELSE
4440  rr1=rrr(10)
4441  ams1=(amp1+amp2+amp3+amp4+amp5)**2
4442  ams2=(amtau-amnuta)**2
4443  alp1=atan((ams1-ama2**2)/ama2/gama2)
4444  alp2=atan((ams2-ama2**2)/ama2/gama2)
4445  alp=alp1+rr1*(alp2-alp1)
4446  am5sq =ama2**2+ama2*gama2*tan(alp)
4447  am5 =sqrt(am5sq)
4448  ENDIF
4449 c --- these are two parts of jacobian, plugged here ---------------
4450  gg5=((am5sq-ama2**2)**2+(ama2*gama2)**2)/(ama2*gama2)
4451  gg5=gg5*(alp2-alp1)
4452  phspac=phspac/(proba2/gg5+(1d0-proba2)/(ams2-ams1) )
4453 c
4454 c mass of 4 (2345)
4455 c flat phase space
4456  rr1=rrr(9)
4457  ams1=(amp2+amp3+amp4+amp5)**2
4458  ams2=(am5-amp1)**2
4459  am4sq=ams1+ rr1*(ams2-ams1)
4460  am4 =sqrt(am4sq)
4461  gg1=ams2-ams1
4462 c
4463 c mass of 3 (234)
4464 
4465  IF (rrr(12).LT.probom) THEN
4466 C phase space with sampling for omega resonance
4467  rr1=rrr(1)
4468  ams1=(amp2+amp3+amp4)**2
4469  ams2=(am4-amp5)**2
4470  alp1=atan((ams1-amom**2)/amom/gamom)
4471  alp2=atan((ams2-amom**2)/amom/gamom)
4472  alp=alp1+rr1*(alp2-alp1)
4473  am3sq =amom**2+amom*gamom*tan(alp)
4474  am3 =sqrt(am3sq)
4475  ELSE
4476 c flat phase space;
4477  rr1=rrr(1)
4478  ams1=(amp2+amp3+amp4)**2
4479  ams2=(am4-amp5)**2
4480  alp1=atan((ams1-amom**2)/amom/gamom)
4481  alp2=atan((ams2-amom**2)/amom/gamom)
4482 
4483  am3sq=ams1+ rr1*(ams2-ams1)
4484  am3 =sqrt(am3sq)
4485 c --- this part of jacobian will be recovered later
4486  ENDIF
4487 c --- this part of the jacobian will be recovered later ---------------
4488  gg2=((am3sq-amom**2)**2+(amom*gamom)**2)/(amom*gamom)
4489  gg2=gg2*(alp2-alp1)
4490  gg2=1d0/(probom/gg2+(1d0-probom)/(ams2-ams1))
4491 c
4492 C mass of 2 (34)
4493  rr2=rrr(2)
4494  ams1=(amp3+amp4)**2
4495  ams2=(am3-amp2)**2
4496 c flat phase space;
4497  am2sq=ams1+ rr2*(ams2-ams1)
4498  am2 =sqrt(am2sq)
4499 c --- this part of jacobian will be recovered later
4500  gg3=ams2-ams1
4501 c
4502 c (34) restframe, define pi3 and pi4
4503  enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4504  enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4505  ppi= enq1**2-amp3**2
4506  pppi=sqrt(abs(enq1**2-amp3**2))
4507  ff1=(4*pi)*(2*pppi/am2)
4508 c pi3 momentum in (34) rest frame
4509  call sphera(pppi,pi3)
4510  pi3(4)=enq1
4511 c pi4 momentum in (34) rest frame
4512  do 30 i=1,3
4513  30 pi4(i)=-pi3(i)
4514  pi4(4)=enq2
4515 c
4516 c (234) rest frame, define pi2
4517 c pr momentum
4518  pr(1)=0
4519  pr(2)=0
4520  pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4521  pr(3)= sqrt(abs(pr(4)**2-am2**2))
4522  ppi = pr(4)**2-am2**2
4523 c pi2 momentum
4524  pi2(1)=0
4525  pi2(2)=0
4526  pi2(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4527  pi2(3)=-pr(3)
4528 c --- this part of jacobian will be recovered later
4529  ff2=(4*pi)*(2*pr(3)/am3)
4530 c old pions boosted from 2 rest frame to 3 rest frame
4531  exe=(pr(4)+pr(3))/am2
4532  call bostr3(exe,pi3,pi3)
4533  call bostr3(exe,pi4,pi4)
4534  rr3=rrr(3)
4535  rr4=rrr(4)
4536  thet =acos(-1.+2*rr3)
4537  phi = 2*pi*rr4
4538  call rotpol(thet,phi,pi2)
4539  call rotpol(thet,phi,pi3)
4540  call rotpol(thet,phi,pi4)
4541 C
4542 C (2345) rest frame, define pi5
4543 c pr momentum
4544  pr(1)=0
4545  pr(2)=0
4546  pr(4)=1./(2*am4)*(am4**2+am3**2-amp5**2)
4547  pr(3)= sqrt(abs(pr(4)**2-am3**2))
4548  ppi = pr(4)**2-am3**2
4549 c pi5 momentum
4550  pi5(1)=0
4551  pi5(2)=0
4552  pi5(4)=1./(2*am4)*(am4**2-am3**2+amp5**2)
4553  pi5(3)=-pr(3)
4554 c --- this part of jacobian will be recovered later
4555  ff3=(4*pi)*(2*pr(3)/am4)
4556 c old pions boosted from 3 rest frame to 4 rest frame
4557  exe=(pr(4)+pr(3))/am3
4558  call bostr3(exe,pi2,pi2)
4559  call bostr3(exe,pi3,pi3)
4560  call bostr3(exe,pi4,pi4)
4561  rr3=rrr(5)
4562  rr4=rrr(6)
4563  thet =acos(-1.+2*rr3)
4564  phi = 2*pi*rr4
4565  call rotpol(thet,phi,pi2)
4566  call rotpol(thet,phi,pi3)
4567  call rotpol(thet,phi,pi4)
4568  call rotpol(thet,phi,pi5)
4569 C
4570 C (12345) rest frame, define pi1
4571 c pr momentum
4572  pr(1)=0
4573  pr(2)=0
4574  pr(4)=1./(2*am5)*(am5**2+am4**2-amp1**2)
4575  pr(3)= sqrt(abs(pr(4)**2-am4**2))
4576  ppi = pr(4)**2-am4**2
4577 c pi1 momentum
4578  pi1(1)=0
4579  pi1(2)=0
4580  pi1(4)=1./(2*am5)*(am5**2-am4**2+amp1**2)
4581  pi1(3)=-pr(3)
4582 c --- this part of jacobian will be recovered later
4583  ff4=(4*pi)*(2*pr(3)/am5)
4584 c old pions boosted from 4 rest frame to 5 rest frame
4585  exe=(pr(4)+pr(3))/am4
4586  call bostr3(exe,pi2,pi2)
4587  call bostr3(exe,pi3,pi3)
4588  call bostr3(exe,pi4,pi4)
4589  call bostr3(exe,pi5,pi5)
4590  rr3=rrr(7)
4591  rr4=rrr(8)
4592  thet =acos(-1.+2*rr3)
4593  phi = 2*pi*rr4
4594  call rotpol(thet,phi,pi1)
4595  call rotpol(thet,phi,pi2)
4596  call rotpol(thet,phi,pi3)
4597  call rotpol(thet,phi,pi4)
4598  call rotpol(thet,phi,pi5)
4599 c
4600 * now to the tau rest frame, define paa and neutrino momenta
4601 * paa momentum
4602  paa(1)=0
4603  paa(2)=0
4604 c paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5**2)
4605 c paa(3)= sqrt(abs(paa(4)**2-am5**2))
4606 c ppi = paa(4)**2-am5**2
4607  paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5sq)
4608  paa(3)= sqrt(abs(paa(4)**2-am5sq))
4609  ppi = paa(4)**2-am5sq
4610  phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4611 * tau-neutrino momentum
4612  pn(1)=0
4613  pn(2)=0
4614  pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am5**2)
4615  pn(3)=-paa(3)
4616 c
4617  phspac=phspac * gg1*gg2*gg3*ff1*ff2*ff3*ff4
4618 c
4619 C all pions boosted from 5 rest frame to tau rest frame
4620 C z-axis antiparallel to neutrino momentum
4621  exe=(paa(4)+paa(3))/am5
4622  call bostr3(exe,pi1,pi1)
4623  call bostr3(exe,pi2,pi2)
4624  call bostr3(exe,pi3,pi3)
4625  call bostr3(exe,pi4,pi4)
4626  call bostr3(exe,pi5,pi5)
4627 c
4628 C partial width consists of phase space and amplitude
4629 C AMPLITUDE (cf YS.Tsai Phys.Rev.D4,2821(1971)
4630 C or F.Gilman SH.Rhie Phys.Rev.D31,1066(1985)
4631 C
4632  pxq=amtau*paa(4)
4633  pxn=amtau*pn(4)
4634  qxn=paa(4)*pn(4)-paa(1)*pn(1)-paa(2)*pn(2)-paa(3)*pn(3)
4635  brak=2*(gv**2+ga**2)*(2*pxq*qxn+am5sq*pxn)
4636  & -6*(gv**2-ga**2)*amtau*amnuta*am5sq
4637  fompp = cabs(bwign(am3,amom,gamom))**2
4638 c normalisation factor (to some numerical undimensioned factor;
4639 c cf R.Fischer et al ZPhys C3, 313 (1980))
4640  fnorm = 1/fpi**6
4641 c AMPLIT=CCABIB**2*GFERMI**2/2. * BRAK * AM5SQ*SIGEE(AM5SQ,JNPI)
4642  amplit=ccabib**2*gfermi**2/2. * brak !* (1D0*(jnpi-12))
4643  amplit = amplit * fompp * fnorm
4644 c phase space test
4645 c amplit = amplit * fnorm
4646 
4647 ! write(*,*) '5pi jnpi=',jnpi
4648 c ignore spin terms
4649  DO 40 i=1,3
4650  40 hv(i)=0.
4651 
4652 ! write(*,*) jnpi
4653 ! stop
4654 
4655  if (jnpi.gt.23) ! for the time being we want to keep old wrong m.e.
4656  $ CALL dam5pi(jnpi,pt,pn,pi1,pi2,pi3,pi4,pi5,amplit,hv)
4657  dgamt=1/(2.*amtau)*amplit*phspac
4658 c
4659  do 77 k=1,4
4660  pmult(k,1)=pi1(k)
4661  pmult(k,2)=pi2(k)
4662  pmult(k,3)=pi3(k)
4663  pmult(k,4)=pi4(k)
4664  pmult(k,5)=pi5(k)
4665  77 continue
4666  return
4667 
4668 C missing: transposition of identical particles, startistical factors
4669 C for identical matrices, polarimetric vector. Matrix element rather naive.
4670 
4671 C flat phase space in pion system + with breit wigner for omega
4672 C anyway it is better than nothing, and code is improvable.
4673  end
4674  SUBROUTINE dphsrk(DGAMT,HV,PN,PR,PMULT,INUM)
4675 C ----------------------------------------------------------------------
4676 C IT SIMULATES RHO DECAY IN TAU REST FRAME WITH
4677 C Z-AXIS ALONG RHO MOMENTUM
4678 C Rho decays to K Kbar
4679 C ----------------------------------------------------------------------
4680  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4681  * ,ampiz,ampi,amro,gamro,ama1,gama1
4682  * ,amk,amkz,amkst,gamkst
4683 C
4684  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
4685  * ,ampiz,ampi,amro,gamro,ama1,gama1
4686  * ,amk,amkz,amkst,gamkst
4687  COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4688  REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
4689  REAL HV(4),PT(4),PN(4),PR(4),PKC(4),PKZ(4),QQ(4),PMULT(4,9)
4690 
4691  REAL RR1(1)
4692 
4693  DATA pi /3.141592653589793238462643/
4694  DATA icont /0/
4695 C
4696 C THREE BODY PHASE SPACE NORMALISED AS IN BJORKEN-DRELL
4697  phspac=1./2**11/pi**5
4698 C TAU MOMENTUM
4699  pt(1)=0.
4700  pt(2)=0.
4701  pt(3)=0.
4702  pt(4)=amtau
4703 C MASS OF (REAL/VIRTUAL) RHO
4704  ams1=(amk+amkz)**2
4705  ams2=(amtau-amnuta)**2
4706 C FLAT PHASE SPACE
4707  CALL ranmar(rr1,1)
4708  amx2=ams1+ rr1(1)*(ams2-ams1)
4709  amx=sqrt(amx2)
4710  phspac=phspac*(ams2-ams1)
4711 C PHASE SPACE WITH SAMPLING FOR RHO RESONANCE
4712 c ALP1=ATAN((AMS1-AMRO**2)/AMRO/GAMRO)
4713 c ALP2=ATAN((AMS2-AMRO**2)/AMRO/GAMRO)
4714 CAM
4715  100 CONTINUE
4716 c CALL RANMAR(RR1,1)
4717 c ALP=ALP1+RR1(1)*(ALP2-ALP1)
4718 c AMX2=AMRO**2+AMRO*GAMRO*TAN(ALP)
4719 c AMX=SQRT(AMX2)
4720 c IF(AMX.LT.(AMK+AMKZ)) GO TO 100
4721 CAM
4722 c PHSPAC=PHSPAC*((AMX2-AMRO**2)**2+(AMRO*GAMRO)**2)/(AMRO*GAMRO)
4723 c PHSPAC=PHSPAC*(ALP2-ALP1)
4724 C
4725 C TAU-NEUTRINO MOMENTUM
4726  pn(1)=0
4727  pn(2)=0
4728  pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
4729  pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
4730 C RHO MOMENTUM
4731  pr(1)=0
4732  pr(2)=0
4733  pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
4734  pr(3)=-pn(3)
4735  phspac=phspac*(4*pi)*(2*pr(3)/amtau)
4736 C
4737 CAM
4738  enq1=(amx2+amk**2-amkz**2)/(2.*amx)
4739  enq2=(amx2-amk**2+amkz**2)/(2.*amx)
4740  pppi=sqrt((enq1-amk)*(enq1+amk))
4741  phspac=phspac*(4*pi)*(2*pppi/amx)
4742 C CHARGED PI MOMENTUM IN RHO REST FRAME
4743  CALL sphera(pppi,pkc)
4744  pkc(4)=enq1
4745 C NEUTRAL PI MOMENTUM IN RHO REST FRAME
4746  DO 20 i=1,3
4747 20 pkz(i)=-pkc(i)
4748  pkz(4)=enq2
4749  exe=(pr(4)+pr(3))/amx
4750 C PIONS BOOSTED FROM RHO REST FRAME TO TAU REST FRAME
4751  CALL bostr3(exe,pkc,pkc)
4752  CALL bostr3(exe,pkz,pkz)
4753  DO 30 i=1,4
4754  30 qq(i)=pkc(i)-pkz(i)
4755 C QQ transverse to PR
4756  pksd =pr(4)*pr(4)-pr(3)*pr(3)-pr(2)*pr(2)-pr(1)*pr(1)
4757  qqpks=pr(4)* qq(4)-pr(3)* qq(3)-pr(2)* qq(2)-pr(1)* qq(1)
4758  DO 31 i=1,4
4759 31 qq(i)=qq(i)-pr(i)*qqpks/pksd
4760 C AMPLITUDE
4761  prodpq=pt(4)*qq(4)
4762  prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
4763  prodpn=pt(4)*pn(4)
4764  qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
4765  brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
4766  & +(gv**2-ga**2)*amtau*amnuta*qq2
4767  amplit=(gfermi*ccabib)**2*brak*2*fpirk(amx)
4768  dgamt=1/(2.*amtau)*amplit*phspac
4769  DO 40 i=1,3
4770  40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
4771  do 77 k=1,4
4772  pmult(k,1)=pkc(k)
4773  pmult(k,2)=pkz(k)
4774  77 continue
4775  RETURN
4776  END
4777  FUNCTION fpirk(W)
4778 C ----------------------------------------------------------
4779 c square of pion form factor
4780 C ----------------------------------------------------------
4781  COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4782  * ,ampiz,ampi,amro,gamro,ama1,gama1
4783  * ,amk,amkz,amkst,gamkst
4784 C
4785  REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
4786  * ,ampiz,ampi,amro,gamro,ama1,gama1
4787  * ,amk,amkz,amkst,gamkst
4788 c COMPLEX FPIKMK
4789  COMPLEX FPIKM
4790  fpirk=cabs(fpikm(w,amk,amkz))**2
4791 c FPIRK=CABS(FPIKMK(W,AMK,AMKZ))**2
4792  END
4793  COMPLEX FUNCTION fpikmk(W,XM1,XM2)
4794 C **********************************************************
4795 C Kaon form factor
4796 C **********************************************************
4797  COMPLEX BWIGM
4798  REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
4799  EXTERNAL bwig
4800  DATA init /0/
4801 C
4802 C ------------ PARAMETERS --------------------
4803  IF (init.EQ.0 ) THEN
4804  init=1
4805  pi=3.141592654
4806  pim=.140
4807  rom=0.773
4808  rog=0.145
4809  rom1=1.570
4810  rog1=0.510
4811 c BETA1=-0.111
4812  beta1=-0.221
4813  ENDIF
4814 C -----------------------------------------------
4815  s=w**2
4816  fpikmk=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
4817  & /(1+beta1)
4818  RETURN
4819  END
4820  SUBROUTINE reslux
4821 C ****************
4822 C INITIALIZE LUND COMMON
4823 
4824 
4825  nhep=0
4826  END
4827  SUBROUTINE dwrph(KTO,PHX)
4828 C
4829 C -------------------------
4830 C
4831  IMPLICIT REAL*8 (a-h,o-z)
4832  REAL*4 PHX(4)
4833  REAL*4 QHOT(4)
4834 C
4835  DO 9 k=1,4
4836  qhot(k) =0.0
4837  9 CONTINUE
4838 C CASE OF TAU RADIATIVE DECAYS.
4839 C FILLING OF THE LUND COMMON BLOCK.
4840  DO 1002 i=1,4
4841  1002 qhot(i)=phx(i)
4842  IF (qhot(4).GT.1.e-5) CALL dwluph(kto,qhot)
4843  RETURN
4844  END
4845  SUBROUTINE dwluph(KTO,PHOT)
4846 C---------------------------------------------------------------------
4847 C Lorentz transformation to CMsystem and
4848 C Updating of HEPEVT record
4849 C
4850 C called by : DEXAY1,(DEKAY1,DEKAY2)
4851 C
4852 C used when radiative corrections in decays are generated
4853 C---------------------------------------------------------------------
4854 C
4855 
4856 
4857  REAL PHOT(4)
4858 
4859  COMMON /taupos/ np1,np2
4860 
4861 C
4862 C check energy
4863  IF (phot(4).LE.0.0) RETURN
4864 C
4865 C position of decaying particle:
4866  IF((kto.EQ. 1).OR.(kto.EQ.11)) THEN
4867  nps=np1
4868  ELSE
4869  nps=np2
4870  ENDIF
4871 C
4872  ktos=kto
4873  IF(ktos.GT.10) ktos=ktos-10
4874 C boost and append photon (gamma is 22)
4875  CALL tralo4(ktos,phot,phot,am)
4876  CALL filhep(0,1,22,nps,nps,0,0,phot,0.0,.true.)
4877 C
4878  RETURN
4879  END
4880 
4881  SUBROUTINE dwluel(KTO,ISGN,PNU,PWB,PEL,PNE)
4882 C ----------------------------------------------------------------------
4883 C Lorentz transformation to CMsystem and
4884 C Updating of HEPEVT record
4885 C
4886 C ISGN = 1/-1 for tau-/tau+
4887 C
4888 C called by : DEXAY,(DEKAY1,DEKAY2)
4889 C ----------------------------------------------------------------------
4890 C
4891 
4892 
4893  REAL PNU(4),PWB(4),PEL(4),PNE(4)
4894 
4895  COMMON /taupos/ np1,np2
4896 
4897 C
4898 C position of decaying particle:
4899  IF(kto.EQ. 1) THEN
4900  nps=np1
4901  ELSE
4902  nps=np2
4903  ENDIF
4904 C
4905 C tau neutrino (nu_tau is 16)
4906  CALL tralo4(kto,pnu,pnu,am)
4907  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4908 C
4909 C W boson (W+ is 24)
4910  CALL tralo4(kto,pwb,pwb,am)
4911 C CALL FILHEP(0,2,-24*ISGN,NPS,NPS,0,0,PWB,AM,.TRUE.)
4912 C
4913 C electron (e- is 11)
4914  CALL tralo4(kto,pel,pel,am)
4915  CALL filhep(0,1,11*isgn,nps,nps,0,0,pel,am,.false.)
4916 C
4917 C anti electron neutrino (nu_e is 12)
4918  CALL tralo4(kto,pne,pne,am)
4919  CALL filhep(0,1,-12*isgn,nps,nps,0,0,pne,am,.true.)
4920 C
4921  RETURN
4922  END
4923  SUBROUTINE dwlumu(KTO,ISGN,PNU,PWB,PMU,PNM)
4924 C ----------------------------------------------------------------------
4925 C Lorentz transformation to CMsystem and
4926 C Updating of HEPEVT record
4927 C
4928 C ISGN = 1/-1 for tau-/tau+
4929 C
4930 C called by : DEXAY,(DEKAY1,DEKAY2)
4931 C ----------------------------------------------------------------------
4932 C
4933 
4934 
4935  REAL PNU(4),PWB(4),PMU(4),PNM(4)
4936 
4937  COMMON /taupos/ np1,np2
4938 
4939 C
4940 C position of decaying particle:
4941  IF(kto.EQ. 1) THEN
4942  nps=np1
4943  ELSE
4944  nps=np2
4945  ENDIF
4946 C
4947 C tau neutrino (nu_tau is 16)
4948  CALL tralo4(kto,pnu,pnu,am)
4949  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4950 C
4951 C W boson (W+ is 24)
4952  CALL tralo4(kto,pwb,pwb,am)
4953 C CALL FILHEP(0,2,-24*ISGN,NPS,NPS,0,0,PWB,AM,.TRUE.)
4954 C
4955 C muon (mu- is 13)
4956  CALL tralo4(kto,pmu,pmu,am)
4957  CALL filhep(0,1,13*isgn,nps,nps,0,0,pmu,am,.false.)
4958 C
4959 C anti muon neutrino (nu_mu is 14)
4960  CALL tralo4(kto,pnm,pnm,am)
4961  CALL filhep(0,1,-14*isgn,nps,nps,0,0,pnm,am,.true.)
4962 C
4963  RETURN
4964  END
4965  SUBROUTINE dwlupi(KTO,ISGN,PPI,PNU)
4966 C ----------------------------------------------------------------------
4967 C Lorentz transformation to CMsystem and
4968 C Updating of HEPEVT record
4969 C
4970 C ISGN = 1/-1 for tau-/tau+
4971 C
4972 C called by : DEXAY,(DEKAY1,DEKAY2)
4973 C ----------------------------------------------------------------------
4974 C
4975  REAL PNU(4),PPI(4)
4976  COMMON /taupos/ np1,np2
4977 C
4978 C position of decaying particle:
4979  IF(kto.EQ. 1) THEN
4980  nps=np1
4981  ELSE
4982  nps=np2
4983  ENDIF
4984 C
4985 C tau neutrino (nu_tau is 16)
4986  CALL tralo4(kto,pnu,pnu,am)
4987  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4988 C
4989 C charged pi meson (pi+ is 211)
4990  CALL tralo4(kto,ppi,ppi,am)
4991  CALL filhep(0,1,-211*isgn,nps,nps,0,0,ppi,am,.true.)
4992 C
4993  RETURN
4994  END
4995  SUBROUTINE dwluro(KTO,ISGN,PNU,PRHO,PIC,PIZ)
4996 C ----------------------------------------------------------------------
4997 C Lorentz transformation to CMsystem and
4998 C Updating of HEPEVT record
4999 C
5000 C ISGN = 1/-1 for tau-/tau+
5001 C
5002 C called by : DEXAY,(DEKAY1,DEKAY2)
5003 C ----------------------------------------------------------------------
5004 C
5005 
5006 
5007  REAL PNU(4),PRHO(4),PIC(4),PIZ(4)
5008 
5009  COMMON /taupos/ np1,np2
5010 
5011 C
5012 C position of decaying particle:
5013  IF(kto.EQ. 1) THEN
5014  nps=np1
5015  ELSE
5016  nps=np2
5017  ENDIF
5018 C
5019 C tau neutrino (nu_tau is 16)
5020  CALL tralo4(kto,pnu,pnu,am)
5021  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5022 C
5023 C charged rho meson (rho+ is 213)
5024  CALL tralo4(kto,prho,prho,am)
5025  CALL filhep(0,2,-213*isgn,nps,nps,0,0,prho,am,.true.)
5026 C
5027 C charged pi meson (pi+ is 211)
5028  CALL tralo4(kto,pic,pic,am)
5029  CALL filhep(0,1,-211*isgn,-1,-1,0,0,pic,am,.true.)
5030 C
5031 C pi0 meson (pi0 is 111)
5032  CALL tralo4(kto,piz,piz,am)
5033  CALL filhep(0,1,111,-2,-2,0,0,piz,am,.true.)
5034 C
5035  RETURN
5036  END
5037  SUBROUTINE dwluaa(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
5038 C ----------------------------------------------------------------------
5039 C Lorentz transformation to CMsystem and
5040 C Updating of HEPEVT record
5041 C
5042 C ISGN = 1/-1 for tau-/tau+
5043 C JAA = 1 (2) FOR A_1- DECAY TO PI+ 2PI- (PI- 2PI0)
5044 C
5045 C called by : DEXAY,(DEKAY1,DEKAY2)
5046 C ----------------------------------------------------------------------
5047 C
5048 
5049 
5050  REAL PNU(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
5051 
5052  COMMON /taupos/ np1,np2
5053 
5054 C
5055 C position of decaying particle:
5056  IF(kto.EQ. 1) THEN
5057  nps=np1
5058  ELSE
5059  nps=np2
5060  ENDIF
5061 C
5062 C tau neutrino (nu_tau is 16)
5063  CALL tralo4(kto,pnu,pnu,am)
5064  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5065 C
5066 C charged a_1 meson (a_1+ is 20213)
5067  CALL tralo4(kto,paa,paa,am)
5068  CALL filhep(0,1,-20213*isgn,nps,nps,0,0,paa,am,.true.)
5069 C
5070 C two possible decays of the charged a1 meson
5071  IF(jaa.EQ.1) THEN
5072 C
5073 C A1 --> PI+ PI- PI- (or charged conjugate)
5074 C
5075 C pi minus (or c.c.) (pi+ is 211)
5076  CALL tralo4(kto,pim2,pim2,am)
5077  CALL filhep(0,1,-211*isgn,-1,-1,0,0,pim2,am,.true.)
5078 C
5079 C pi minus (or c.c.) (pi+ is 211)
5080  CALL tralo4(kto,pim1,pim1,am)
5081  CALL filhep(0,1,-211*isgn,-2,-2,0,0,pim1,am,.true.)
5082 C
5083 C pi plus (or c.c.) (pi+ is 211)
5084  CALL tralo4(kto,pipl,pipl,am)
5085  CALL filhep(0,1, 211*isgn,-3,-3,0,0,pipl,am,.true.)
5086 C
5087  ELSE IF (jaa.EQ.2) THEN
5088 C
5089 C A1 --> PI- PI0 PI0 (or charged conjugate)
5090 C
5091 C pi zero (pi0 is 111)
5092  CALL tralo4(kto,pim2,pim2,am)
5093  CALL filhep(0,1,111,-1,-1,0,0,pim2,am,.true.)
5094 C
5095 C pi zero (pi0 is 111)
5096  CALL tralo4(kto,pim1,pim1,am)
5097  CALL filhep(0,1,111,-2,-2,0,0,pim1,am,.true.)
5098 C
5099 C pi minus (or c.c.) (pi+ is 211)
5100  CALL tralo4(kto,pipl,pipl,am)
5101  CALL filhep(0,1,-211*isgn,-3,-3,0,0,pipl,am,.true.)
5102 C
5103  ENDIF
5104 C
5105  RETURN
5106  END
5107  SUBROUTINE dwlukk (KTO,ISGN,PKK,PNU)
5108 C ----------------------------------------------------------------------
5109 C Lorentz transformation to CMsystem and
5110 C Updating of HEPEVT record
5111 C
5112 C ISGN = 1/-1 for tau-/tau+
5113 C
5114 C ----------------------------------------------------------------------
5115 C
5116  REAL PKK(4),PNU(4)
5117  COMMON /taupos/ np1,np2
5118 C
5119 C position of decaying particle
5120 
5121  IF (kto.EQ.1) THEN
5122 
5123  nps=np1
5124  ELSE
5125  nps=np2
5126  ENDIF
5127 C
5128 C tau neutrino (nu_tau is 16)
5129  CALL tralo4 (kto,pnu,pnu,am)
5130  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5131 C
5132 C K meson (K+ is 321)
5133  CALL tralo4 (kto,pkk,pkk,am)
5134  CALL filhep(0,1,-321*isgn,nps,nps,0,0,pkk,am,.true.)
5135 C
5136  RETURN
5137  END
5138  SUBROUTINE dwluks(KTO,ISGN,PNU,PKS,PKK,PPI,JKST)
5139  COMMON / taukle / bra1,brk0,brk0b,brks
5140  REAL*4 BRA1,BRK0,BRK0B,BRKS
5141 
5142 C ----------------------------------------------------------------------
5143 C Lorentz transformation to CMsystem and
5144 C Updating of HEPEVT record
5145 C
5146 C ISGN = 1/-1 for tau-/tau+
5147 C JKST=10 (20) corresponds to K0B pi- (K- pi0) decay
5148 C
5149 C ----------------------------------------------------------------------
5150 C
5151 
5152  REAL PNU(4),PKS(4),PKK(4),PPI(4),XIO(1)
5153  COMMON /taupos/ np1,np2
5154 
5155 C
5156 C position of decaying particle
5157  IF(kto.EQ. 1) THEN
5158  nps=np1
5159  ELSE
5160  nps=np2
5161  ENDIF
5162 C
5163 C tau neutrino (nu_tau is 16)
5164  CALL tralo4(kto,pnu,pnu,am)
5165  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5166 C
5167 C charged K* meson (K*+ is 323)
5168  CALL tralo4(kto,pks,pks,am)
5169  CALL filhep(0,1,-323*isgn,nps,nps,0,0,pks,am,.true.)
5170 C
5171 C two possible decay modes of charged K*
5172  IF(jkst.EQ.10) THEN
5173 C
5174 C K*- --> pi- K0B (or charged conjugate)
5175 C
5176 C charged pi meson (pi+ is 211)
5177  CALL tralo4(kto,ppi,ppi,am)
5178  CALL filhep(0,1,-211*isgn,-1,-1,0,0,ppi,am,.true.)
5179 C
5180  bran=brk0b
5181  IF (isgn.EQ.-1) bran=brk0
5182 C K0 --> K0_long (is 130) / K0_short (is 310) = 1/1
5183  CALL ranmar(xio,1)
5184  IF(xio(1).GT.bran) THEN
5185  k0type = 130
5186  ELSE
5187  k0type = 310
5188  ENDIF
5189 C
5190  CALL tralo4(kto,pkk,pkk,am)
5191  CALL filhep(0,1,k0type,-2,-2,0,0,pkk,am,.true.)
5192 C
5193  ELSE IF(jkst.EQ.20) THEN
5194 C
5195 C K*- --> pi0 K-
5196 C
5197 C pi zero (pi0 is 111)
5198  CALL tralo4(kto,ppi,ppi,am)
5199  CALL filhep(0,1,111,-1,-1,0,0,ppi,am,.true.)
5200 C
5201 C charged K meson (K+ is 321)
5202  CALL tralo4(kto,pkk,pkk,am)
5203  CALL filhep(0,1,-321*isgn,-2,-2,0,0,pkk,am,.true.)
5204 C
5205  ENDIF
5206 C
5207  RETURN
5208  END
5209  SUBROUTINE dwlnew(KTO,ISGN,PNU,PWB,PNPI,MODE)
5210 C ----------------------------------------------------------------------
5211 C Lorentz transformation to CMsystem and
5212 C Updating of HEPEVT record
5213 C
5214 C ISGN = 1/-1 for tau-/tau+
5215 C
5216 C called by : DEXAY,(DEKAY1,DEKAY2)
5217 C ----------------------------------------------------------------------
5218 C
5219  parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
5220 
5221  COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5222 
5223  & ,names
5224  COMMON /taupos/ np1,np2
5225  CHARACTER NAMES(NMODE)*31
5226  REAL PNU(4),PWB(4),PNPI(4,9)
5227  REAL PPI(4)
5228 C
5229  jnpi=mode-7
5230 C position of decaying particle
5231  IF(kto.EQ. 1) THEN
5232  nps=np1
5233  ELSE
5234  nps=np2
5235  ENDIF
5236 C
5237 C tau neutrino (nu_tau is 16)
5238  CALL tralo4(kto,pnu,pnu,am)
5239  CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5240 C
5241 C W boson (W+ is 24)
5242  CALL tralo4(kto,pwb,pwb,am)
5243  CALL filhep(0,1,-24*isgn,nps,nps,0,0,pwb,am,.true.)
5244 C
5245 C multi pi mode JNPI
5246 C
5247 C get multiplicity of mode JNPI
5248  nd=mulpik(jnpi)
5249  DO i=1,nd
5250 
5251  kfpi=lunpik(idffin(i,jnpi),-isgn)
5252 
5253 C for charged conjugate case, change charged pions only
5254 C IF(KFPI.NE.111)KFPI=KFPI*ISGN
5255  DO j=1,4
5256  ppi(j)=pnpi(j,i)
5257  END DO
5258  CALL tralo4(kto,ppi,ppi,am)
5259  CALL filhep(0,1,kfpi,-i,-i,0,0,ppi,am,.true.)
5260  END DO
5261 C
5262  RETURN
5263  END
5264 
5265 
5266  FUNCTION amast(PP)
5267 C ----------------------------------------------------------------------
5268 C CALCULATES MASS OF PP (DOUBLE PRECISION)
5269 C
5270 C USED BY : RADKOR
5271 C ----------------------------------------------------------------------
5272  IMPLICIT REAL*8 (a-h,o-z)
5273  REAL*8 PP(4)
5274  aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5275 C
5276  IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5277  amast=aaa
5278  RETURN
5279  END
5280  FUNCTION amas4(PP)
5281 C ******************
5282 C ----------------------------------------------------------------------
5283 C CALCULATES MASS OF PP
5284 C
5285 C USED BY :
5286 C ----------------------------------------------------------------------
5287  REAL PP(4)
5288  aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5289  IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5290  amas4=aaa
5291  RETURN
5292  END
5293  FUNCTION angxy(X,Y)
5294 C ----------------------------------------------------------------------
5295 C
5296 C USED BY : KORALZ RADKOR
5297 C ----------------------------------------------------------------------
5298  IMPLICIT DOUBLE PRECISION (a-h,o-z)
5299  DATA pi /3.141592653589793238462643d0/
5300 C
5301  IF(abs(y).LT.abs(x)) THEN
5302  the=atan(abs(y/x))
5303  IF(x.LE.0d0) the=pi-the
5304  ELSE
5305  the=acos(x/sqrt(x**2+y**2))
5306  ENDIF
5307  angxy=the
5308  RETURN
5309  END
5310  FUNCTION angfi(X,Y)
5311 C ----------------------------------------------------------------------
5312 * CALCULATES ANGLE IN (0,2*PI) RANGE OUT OF X-Y
5313 C
5314 C USED BY : KORALZ RADKOR
5315 C ----------------------------------------------------------------------
5316  IMPLICIT DOUBLE PRECISION (a-h,o-z)
5317  DATA pi /3.141592653589793238462643d0/
5318 C
5319  IF(abs(y).LT.abs(x)) THEN
5320  the=atan(abs(y/x))
5321  IF(x.LE.0d0) the=pi-the
5322  ELSE
5323  the=acos(x/sqrt(x**2+y**2))
5324  ENDIF
5325  IF(y.LT.0d0) the=2d0*pi-the
5326  angfi=the
5327  END
5328  SUBROUTINE rotod1(PH1,PVEC,QVEC)
5329 C ----------------------------------------------------------------------
5330 C
5331 C USED BY : KORALZ
5332 C ----------------------------------------------------------------------
5333  IMPLICIT DOUBLE PRECISION (a-h,o-z)
5334  dimension pvec(4),qvec(4),rvec(4)
5335 C
5336  phi=ph1
5337  cs=cos(phi)
5338  sn=sin(phi)
5339  DO 10 i=1,4
5340  10 rvec(i)=pvec(i)
5341  qvec(1)=rvec(1)
5342  qvec(2)= cs*rvec(2)-sn*rvec(3)
5343  qvec(3)= sn*rvec(2)+cs*rvec(3)
5344  qvec(4)=rvec(4)
5345  RETURN
5346  END
5347  SUBROUTINE rotod2(PH1,PVEC,QVEC)
5348 C ----------------------------------------------------------------------
5349 C
5350 C USED BY : KORALZ RADKOR
5351 C ----------------------------------------------------------------------
5352  IMPLICIT DOUBLE PRECISION (a-h,o-z)
5353  dimension pvec(4),qvec(4),rvec(4)
5354 C
5355  phi=ph1
5356  cs=cos(phi)
5357  sn=sin(phi)
5358  DO 10 i=1,4
5359  10 rvec(i)=pvec(i)
5360  qvec(1)= cs*rvec(1)+sn*rvec(3)
5361  qvec(2)=rvec(2)
5362  qvec(3)=-sn*rvec(1)+cs*rvec(3)
5363  qvec(4)=rvec(4)
5364  RETURN
5365  END
5366  SUBROUTINE rotod3(PH1,PVEC,QVEC)
5367 C ----------------------------------------------------------------------
5368 C
5369 C USED BY : KORALZ RADKOR
5370 C ----------------------------------------------------------------------
5371  IMPLICIT DOUBLE PRECISION (a-h,o-z)
5372 C
5373  dimension pvec(4),qvec(4),rvec(4)
5374  phi=ph1
5375  cs=cos(phi)
5376  sn=sin(phi)
5377  DO 10 i=1,4
5378  10 rvec(i)=pvec(i)
5379  qvec(1)= cs*rvec(1)-sn*rvec(2)
5380  qvec(2)= sn*rvec(1)+cs*rvec(2)
5381  qvec(3)=rvec(3)
5382  qvec(4)=rvec(4)
5383  END
5384  SUBROUTINE bostr3(EXE,PVEC,QVEC)
5385 C ----------------------------------------------------------------------
5386 C BOOST ALONG Z AXIS, EXE=EXP(ETA), ETA= HIPERBOLIC VELOCITY.
5387 C
5388 C USED BY : TAUOLA KORALZ (?)
5389 C ----------------------------------------------------------------------
5390  REAL*4 PVEC(4),QVEC(4),RVEC(4)
5391 C
5392  DO 10 i=1,4
5393  10 rvec(i)=pvec(i)
5394  rpl=rvec(4)+rvec(3)
5395  rmi=rvec(4)-rvec(3)
5396  qpl=rpl*exe
5397  qmi=rmi/exe
5398  qvec(1)=rvec(1)
5399  qvec(2)=rvec(2)
5400  qvec(3)=(qpl-qmi)/2
5401  qvec(4)=(qpl+qmi)/2
5402  END
5403  SUBROUTINE bostd3(EXE,PVEC,QVEC)
5404 C ----------------------------------------------------------------------
5405 C BOOST ALONG Z AXIS, EXE=EXP(ETA), ETA= HIPERBOLIC VELOCITY.
5406 C
5407 C USED BY : KORALZ RADKOR
5408 C ----------------------------------------------------------------------
5409  IMPLICIT DOUBLE PRECISION (a-h,o-z)
5410  dimension pvec(4),qvec(4),rvec(4)
5411 C
5412  DO 10 i=1,4
5413  10 rvec(i)=pvec(i)
5414  rpl=rvec(4)+rvec(3)
5415  rmi=rvec(4)-rvec(3)
5416  qpl=rpl*exe
5417  qmi=rmi/exe
5418  qvec(1)=rvec(1)
5419  qvec(2)=rvec(2)
5420  qvec(3)=(qpl-qmi)/2
5421  qvec(4)=(qpl+qmi)/2
5422  RETURN
5423  END
5424  SUBROUTINE rotor1(PH1,PVEC,QVEC)
5425 C ----------------------------------------------------------------------
5426 C
5427 C called by :
5428 C ----------------------------------------------------------------------
5429  REAL*4 PVEC(4),QVEC(4),RVEC(4)
5430 C
5431  phi=ph1
5432  cs=cos(phi)
5433  sn=sin(phi)
5434  DO 10 i=1,4
5435  10 rvec(i)=pvec(i)
5436  qvec(1)=rvec(1)
5437  qvec(2)= cs*rvec(2)-sn*rvec(3)
5438  qvec(3)= sn*rvec(2)+cs*rvec(3)
5439  qvec(4)=rvec(4)
5440  END
5441  SUBROUTINE rotor2(PH1,PVEC,QVEC)
5442 C ----------------------------------------------------------------------
5443 C
5444 C USED BY : TAUOLA
5445 C ----------------------------------------------------------------------
5446  IMPLICIT REAL*4(a-h,o-z)
5447  REAL*4 PVEC(4),QVEC(4),RVEC(4)
5448 C
5449  phi=ph1
5450  cs=cos(phi)
5451  sn=sin(phi)
5452  DO 10 i=1,4
5453  10 rvec(i)=pvec(i)
5454  qvec(1)= cs*rvec(1)+sn*rvec(3)
5455  qvec(2)=rvec(2)
5456  qvec(3)=-sn*rvec(1)+cs*rvec(3)
5457  qvec(4)=rvec(4)
5458  END
5459  SUBROUTINE rotor3(PHI,PVEC,QVEC)
5460 C ----------------------------------------------------------------------
5461 C
5462 C USED BY : TAUOLA
5463 C ----------------------------------------------------------------------
5464  REAL*4 PVEC(4),QVEC(4),RVEC(4)
5465 C
5466  cs=cos(phi)
5467  sn=sin(phi)
5468  DO 10 i=1,4
5469  10 rvec(i)=pvec(i)
5470  qvec(1)= cs*rvec(1)-sn*rvec(2)
5471  qvec(2)= sn*rvec(1)+cs*rvec(2)
5472  qvec(3)=rvec(3)
5473  qvec(4)=rvec(4)
5474  END
5475  SUBROUTINE spherd(R,X)
5476 C ----------------------------------------------------------------------
5477 C GENERATES UNIFORMLY THREE-VECTOR X ON SPHERE OF RADIUS R
5478 C DOUBLE PRECISON VERSION OF SPHERA
5479 C ----------------------------------------------------------------------
5480  REAL*8 R,X(4),PI,COSTH,SINTH
5481  REAL*4 RRR(2)
5482  DATA pi /3.141592653589793238462643d0/
5483 C
5484  CALL ranmar(rrr,2)
5485  costh=-1+2*rrr(1)
5486  sinth=sqrt(1 -costh**2)
5487  x(1)=r*sinth*cos(2*pi*rrr(2))
5488  x(2)=r*sinth*sin(2*pi*rrr(2))
5489  x(3)=r*costh
5490  RETURN
5491  END
5492  SUBROUTINE rotpox(THET,PHI,PP)
5493  IMPLICIT REAL*8 (a-h,o-z)
5494 C ----------------------------------------------------------------------
5495 
5496 C
5497 
5498 C ----------------------------------------------------------------------
5499  dimension pp(4)
5500 C
5501  CALL rotod2(thet,pp,pp)
5502  CALL rotod3( phi,pp,pp)
5503  RETURN
5504  END
5505  SUBROUTINE sphera(R,X)
5506 C ----------------------------------------------------------------------
5507 C GENERATES UNIFORMLY THREE-VECTOR X ON SPHERE OF RADIUS R
5508 C
5509 C called by : DPHSxx,DADMPI,DADMKK
5510 C ----------------------------------------------------------------------
5511  REAL X(4)
5512  REAL*4 RRR(2)
5513  DATA pi /3.141592653589793238462643/
5514 C
5515  CALL ranmar(rrr,2)
5516  costh=-1.+2.*rrr(1)
5517  sinth=sqrt(1.-costh**2)
5518  x(1)=r*sinth*cos(2*pi*rrr(2))
5519  x(2)=r*sinth*sin(2*pi*rrr(2))
5520  x(3)=r*costh
5521  RETURN
5522  END
5523  SUBROUTINE rotpol(THET,PHI,PP)
5524 C ----------------------------------------------------------------------
5525 C
5526 C called by : DADMAA,DPHSAA
5527 C ----------------------------------------------------------------------
5528  REAL PP(4)
5529 C
5530  CALL rotor2(thet,pp,pp)
5531  CALL rotor3( phi,pp,pp)
5532  RETURN
5533  END
5534  SUBROUTINE ranmar(RVEC,LENV)
5535 C ----------------------------------------------------------------------
5536 C<<<<<FUNCTION RANMAR(IDUMM)
5537 C CERNLIB V113, VERSION WITH AUTOMATIC DEFAULT INITIALIZATION
5538 C Transformed to SUBROUTINE to be as in CERNLIB
5539 C AM.Lutz November 1988, Feb. 1989
5540 C
5541 C!Universal random number generator proposed by Marsaglia and Zaman
5542 C in report FSU-SCRI-87-50
5543 C modified by F. James, 1988 and 1989, to generate a vector
5544 C of pseudorandom numbers RVEC of length LENV, and to put in
5545 C the COMMON block everything needed to specify currrent state,
5546 C and to add input and output entry points RMARIN, RMARUT.
5547 C
5548 C Unique random number used in the program
5549 C ----------------------------------------------------------------------
5550  COMMON / inout / inut,iout
5551  dimension rvec(*)
5552  common/raset1/u(97),c,i97,j97
5553  parameter(modcns=1000000000)
5554  DATA ntot,ntot2,ijkl/-1,0,0/
5555 C
5556  IF (ntot .GE. 0) GO TO 50
5557 C
5558 C Default initialization. User has called RANMAR without RMARIN.
5559  ijkl = 54217137
5560  ntot = 0
5561  ntot2 = 0
5562  kalled = 0
5563  GO TO 1
5564 C
5565  entry rmarin(ijklin, ntotin,ntot2n)
5566 C Initializing routine for RANMAR, may be called before
5567 C generating pseudorandom numbers with RANMAR. The input
5568 C values should be in the ranges: 0<=IJKLIN<=900 OOO OOO
5569 C 0<=NTOTIN<=999 999 999
5570 C 0<=NTOT2N<<999 999 999!
5571 C To get the standard values in Marsaglia-s paper, IJKLIN=54217137
5572 C NTOTIN,NTOT2N=0
5573  ijkl = ijklin
5574  ntot = max(ntotin,0)
5575  ntot2= max(ntot2n,0)
5576  kalled = 1
5577 C always come here to initialize
5578  1 CONTINUE
5579  ij = ijkl/30082
5580  kl = ijkl - 30082*ij
5581  i = mod(ij/177, 177) + 2
5582  j = mod(ij, 177) + 2
5583  k = mod(kl/169, 178) + 1
5584  l = mod(kl, 169)
5585  WRITE(iout,201) ijkl,ntot,ntot2
5586  201 FORMAT(1x,' RANMAR INITIALIZED: ',i10,2x,2i10)
5587  DO 2 ii= 1, 97
5588  s = 0.
5589  t = .5
5590  DO 3 jj= 1, 24
5591  m = mod(mod(i*j,179)*k, 179)
5592  i = j
5593  j = k
5594  k = m
5595  l = mod(53*l+1, 169)
5596  IF (mod(l*m,64) .GE. 32) s = s+t
5597  3 t = 0.5*t
5598  2 u(ii) = s
5599  twom24 = 1.0
5600  DO 4 i24= 1, 24
5601  4 twom24 = 0.5*twom24
5602  c = 362436.*twom24
5603  cd = 7654321.*twom24
5604  cm = 16777213.*twom24
5605  i97 = 97
5606  j97 = 33
5607 C Complete initialization by skipping
5608 C (NTOT2*MODCNS + NTOT) random numbers
5609  DO 45 loop2= 1, ntot2+1
5610  now = modcns
5611  IF (loop2 .EQ. ntot2+1) now=ntot
5612  IF (now .GT. 0) THEN
5613  WRITE (iout,'(A,I15)') ' RMARIN SKIPPING OVER ',now
5614  DO 40 idum = 1, ntot
5615  uni = u(i97)-u(j97)
5616  IF (uni .LT. 0.) uni=uni+1.
5617  u(i97) = uni
5618  i97 = i97-1
5619  IF (i97 .EQ. 0) i97=97
5620  j97 = j97-1
5621  IF (j97 .EQ. 0) j97=97
5622  c = c - cd
5623  IF (c .LT. 0.) c=c+cm
5624  40 CONTINUE
5625  ENDIF
5626  45 CONTINUE
5627  IF (kalled .EQ. 1) RETURN
5628 C
5629 C Normal entry to generate LENV random numbers
5630  50 CONTINUE
5631  DO 100 ivec= 1, lenv
5632  uni = u(i97)-u(j97)
5633  IF (uni .LT. 0.) uni=uni+1.
5634  u(i97) = uni
5635  i97 = i97-1
5636  IF (i97 .EQ. 0) i97=97
5637  j97 = j97-1
5638  IF (j97 .EQ. 0) j97=97
5639  c = c - cd
5640  IF (c .LT. 0.) c=c+cm
5641  uni = uni-c
5642  IF (uni .LT. 0.) uni=uni+1.
5643 C Replace exact zeroes by uniform distr. *2**-24
5644  IF (uni .EQ. 0.) THEN
5645  uni = twom24*u(2)
5646 C An exact zero here is very unlikely, but lets be safe.
5647  IF (uni .EQ. 0.) uni= twom24*twom24
5648  ENDIF
5649  rvec(ivec) = uni
5650  100 CONTINUE
5651  ntot = ntot + lenv
5652  IF (ntot .GE. modcns) THEN
5653  ntot2 = ntot2 + 1
5654  ntot = ntot - modcns
5655  ENDIF
5656  RETURN
5657 C Entry to output current status
5658  entry rmarut(ijklut,ntotut,ntot2t)
5659  ijklut = ijkl
5660  ntotut = ntot
5661  ntot2t = ntot2
5662  RETURN
5663  END
5664  FUNCTION dilogt(X)
5665 C *****************
5666  IMPLICIT REAL*8(a-h,o-z)
5667 CERN C304 VERSION 29/07/71 DILOG 59 C
5668  z=-1.64493406684822
5669  IF(x .LT.-1.0) GO TO 1
5670  IF(x .LE. 0.5) GO TO 2
5671  IF(x .EQ. 1.0) GO TO 3
5672  IF(x .LE. 2.0) GO TO 4
5673  z=3.2898681336964
5674  1 t=1.0/x
5675  s=-0.5
5676  z=z-0.5* log(abs(x))**2
5677  GO TO 5
5678  2 t=x
5679  s=0.5
5680  z=0.
5681  GO TO 5
5682  3 dilogt=1.64493406684822
5683  RETURN
5684  4 t=1.0-x
5685  s=-0.5
5686  z=1.64493406684822 - log(x)* log(abs(t))
5687  5 y=2.66666666666666 *t+0.66666666666666
5688  b= 0.00000 00000 00001
5689  a=y*b +0.00000 00000 00004
5690  b=y*a-b+0.00000 00000 00011
5691  a=y*b-a+0.00000 00000 00037
5692  b=y*a-b+0.00000 00000 00121
5693  a=y*b-a+0.00000 00000 00398
5694  b=y*a-b+0.00000 00000 01312
5695  a=y*b-a+0.00000 00000 04342
5696  b=y*a-b+0.00000 00000 14437
5697  a=y*b-a+0.00000 00000 48274
5698  b=y*a-b+0.00000 00001 62421
5699  a=y*b-a+0.00000 00005 50291
5700  b=y*a-b+0.00000 00018 79117
5701  a=y*b-a+0.00000 00064 74338
5702  b=y*a-b+0.00000 00225 36705
5703  a=y*b-a+0.00000 00793 87055
5704  b=y*a-b+0.00000 02835 75385
5705  a=y*b-a+0.00000 10299 04264
5706  b=y*a-b+0.00000 38163 29463
5707  a=y*b-a+0.00001 44963 00557
5708  b=y*a-b+0.00005 68178 22718
5709  a=y*b-a+0.00023 20021 96094
5710  b=y*a-b+0.00100 16274 96164
5711  a=y*b-a+0.00468 63619 59447
5712  b=y*a-b+0.02487 93229 24228
5713  a=y*b-a+0.16607 30329 27855
5714  a=y*a-b+1.93506 43008 6996
5715  dilogt=s*t*(a-b)+z
5716  RETURN
5717 C=======================================================================
5718 C===================END OF CPC PART ====================================
5719 C=======================================================================
5720  END
5721