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