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