C++ Interface to Tauola
tauola-F/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.7 ******** *
11C **************DECEMBER 1993****************** *
12#else
13C *********TAUOLA LIBRARY: VERSION 2.7 ******** *
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,hdum,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,hdum,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,hdum,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,hdum,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.7 ******',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.7 ******',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.7 ******',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.7 ******',9x,1h*,
258 $ /,' *', 25x,'***********DECEMBER 1993***************',9x,1h*,
259#else
260 $ /,' *', 25x,'*****TAUOLA LIBRARY: VERSION 2.7 ******',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.7 ******',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.7 ******',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.7 ******',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.7 ******',9x,1h*,
653 $ /,' *', 25x,'***********DECEMBER 1993***************',9x,1h*,
654#else
655 $ /,' *', 25x,'*****TAUOLA LIBRARY: VERSION 2.7 ******',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 ZBW 20.12.2002 bug fix
4410 IF(rrr(9).LE.0.5*prez) THEN
4411 DO 72 i=1,4
4412 x=pim1(i)
4413 pim1(i)=pim2(i)
4414 72 pim2(i)=x
4415 ENDIF
4416C end of bug fix
4417C WE INCLUDE REMAINING PART OF THE JACOBIAN
4418C --- FLAT CHANNEL
4419 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4420 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4421 ams2=(am4-amp2)**2
4422 ams1=(amp1+amp3+amp4)**2
4423 ff1=(ams2-ams1)
4424 ams1=(amp3+amp4)**2
4425 ams2=(sqrt(am3sq)-amp1)**2
4426 ff2=ams2-ams1
4427 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4428 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4429 uu=ff1*ff2*ff3*ff4
4430C --- FIRST CHANNEL
4431 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4432 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4433 ams2=(am4-amp2)**2
4434 ams1=(amp1+amp3+amp4)**2
4435 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4436 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4437 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4438 ff1=ff1*(alp2-alp1)
4439 ams1=(amp3+amp4)**2
4440 ams2=(sqrt(am3sq)-amp1)**2
4441 ff2=ams2-ams1
4442 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4443 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4444 ff=ff1*ff2*ff3*ff4
4445C --- SECOND CHANNEL
4446 am3sq=(pim2(4)+piz(4)+pipl(4))**2-(pim2(3)+piz(3)+pipl(3))**2
4447 $ -(pim2(2)+piz(2)+pipl(2))**2-(pim2(1)+piz(1)+pipl(1))**2
4448 ams2=(am4-amp1)**2
4449 ams1=(amp2+amp3+amp4)**2
4450 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4451 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4452 gg1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4453 gg1=gg1*(alp2-alp1)
4454 ams1=(amp3+amp4)**2
4455 ams2=(sqrt(am3sq)-amp2)**2
4456 gg2=ams2-ams1
4457 gg3=(4*pi)*(xlam(am2**2,amp2**2,am3sq)/am3sq)
4458 gg4=(4*pi)*(xlam(am3sq,amp1**2,am4**2)/am4**2)
4459 gg=gg1*gg2*gg3*gg4
4460C --- JACOBIAN AVERAGED OVER THE TWO
4461 IF ( ( (ff+gg)*uu+ff*gg ).GT.0.0d0) THEN
4462 rr=ff*gg*uu/(0.5*prez*(ff+gg)*uu+(1.0-prez)*ff*gg)
4463 phspac=phspac*rr
4464 ELSE
4465 phspac=0.0
4466 ENDIF
4467* MOMENTA OF THE TWO PI-MINUS ARE RANDOMLY SYMMETRISED
4468 IF (jnpi.EQ.1) THEN
4469 rr5= rrr(5)
4470 IF(rr5.LE.0.5) THEN
4471 DO 70 i=1,4
4472 x=pim1(i)
4473 pim1(i)=pim2(i)
4474 70 pim2(i)=x
4475 ENDIF
4476 phspac=phspac/2.
4477 ELSE
4478C MOMENTA OF PI0-S ARE GENERATED UNIFORMLY ONLY IF PREZ=0.0
4479 rr5= rrr(5)
4480 IF(rr5.LE.0.5) THEN
4481 DO 71 i=1,4
4482 x=pim1(i)
4483 pim1(i)=pim2(i)
4484 71 pim2(i)=x
4485 ENDIF
4486 phspac=phspac/6.
4487 ENDIF
4488* ALL PIONS BOOSTED FROM 4 REST FRAME TO TAU REST FRAME
4489* Z-AXIS ANTIPARALLEL TO NEUTRINO MOMENTUM
4490 exe=(paa(4)+paa(3))/am4
4491 CALL bostr3(exe,piz,piz)
4492 CALL bostr3(exe,pipl,pipl)
4493 CALL bostr3(exe,pim1,pim1)
4494 CALL bostr3(exe,pim2,pim2)
4495 CALL bostr3(exe,pr,pr)
4496C PARTIAL WIDTH CONSISTS OF PHASE SPACE AND AMPLITUDE
4497C CHECK ON CONSISTENCY WITH DADNPI, THEN, CODE BREAKES UNIFORM PION
4498C DISTRIBUTION IN HADRONIC SYSTEM
4499#if defined (ALEPH)
4500 CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4501#else
4502CAM Assume neutrino mass=0. and sum over final polarisation
4503C AMX2=AM4**2
4504C BRAK= 2*(AMTAU**2-AMX2) * (AMTAU**2+2.*AMX2)
4505C AMPLIT=CCABIB**2*GFERMI**2/2. * BRAK * AMX2*SIGEE(AMX2,1)
4506 IF (jnpi.EQ.1) THEN
4507 CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4508 ELSEIF (jnpi.EQ.2) THEN
4509 CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv)
4510 ENDIF
4511#endif
4512 dgamt=1/(2.*amtau)*amplit*phspac
4513C PHASE SPACE CHECK
4514C DGAMT=PHSPAC
4515 DO 77 k=1,4
4516 pmult(k,1)=pim1(k)
4517 pmult(k,2)=pim2(k)
4518#if defined (ALEPH)
4519 pmult(k,3)=piz(k)
4520 pmult(k,4)=pipl(k)
4521#else
4522 pmult(k,3)=pipl(k)
4523 pmult(k,4)=piz(k)
4524#endif
4525 77 CONTINUE
4526 END
4527 SUBROUTINE dam4pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,AMPLIT,HV)
4528C ----------------------------------------------------------------------
4529* CALCULATES DIFFERENTIAL CROSS SECTION AND POLARIMETER VECTOR
4530* FOR TAU DECAY INTO 4 PI MODES
4531* ALL SPIN EFFECTS IN THE FULL DECAY CHAIN ARE TAKEN INTO ACCOUNT.
4532* CALCULATIONS DONE IN TAU REST FRAME WITH Z-AXIS ALONG NEUTRINO MOMENT
4533C MNUM DECAY MODE IDENTIFIER.
4534C
4535#if defined (ALEPH)
4536C called by : DPH4PI
4537#else
4538C called by : DPHSAA
4539#endif
4540C ----------------------------------------------------------------------
4541 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4542 * ,ampiz,ampi,amro,gamro,ama1,gama1
4543 * ,amk,amkz,amkst,gamkst
4544C
4545 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4546 * ,ampiz,ampi,amro,gamro,ama1,gama1
4547 * ,amk,amkz,amkst,gamkst
4548 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4549 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4550 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4),PIM4(4)
4551 REAL PIVEC(4),PIAKS(4),HVM(4)
4552 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5
4553 EXTERNAL form1,form2,form3,form4,form5
4554 DATA pi /3.141592653589793238462643/
4555 DATA icont /0/
4556C
4557#if defined (CLEO)
4558 CALL curr_cleo(mnum,pim1,pim2,pim3,pim4,hadcur)
4559#else
4560 CALL curr(mnum,pim1,pim2,pim3,pim4,hadcur)
4561#endif
4562C
4563* CALCULATE PI-VECTORS: VECTOR AND AXIAL
4564 CALL clvec(hadcur,pn,pivec)
4565 CALL claxi(hadcur,pn,piaks)
4566 CALL clnut(hadcur,brakm,hvm)
4567* SPIN INDEPENDENT PART OF DECAY DIFF-CROSS-SECT. IN TAU REST FRAME
4568 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4569 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4570 amplit=(ccabib*gfermi)**2*brak/2.
4571C POLARIMETER VECTOR IN TAU REST FRAME
4572 DO 90 i=1,3
4573 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4574 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4575C HV IS DEFINED FOR TAU- WITH GAMMA=B+HV*POL
4576 IF (brak.NE.0.0)
4577 &hv(i)=-hv(i)/brak
4578 90 CONTINUE
4579 END
4580 SUBROUTINE dph5pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4581C ----------------------------------------------------------------------
4582* IT SIMULATES 5pi DECAY IN TAU REST FRAME WITH
4583* Z-AXIS ALONG 5pi MOMENTUM
4584C ----------------------------------------------------------------------
4585 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4586 * ,ampiz,ampi,amro,gamro,ama1,gama1
4587 * ,amk,amkz,amkst,gamkst
4588C
4589 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4590 * ,ampiz,ampi,amro,gamro,ama1,gama1
4591
4592
4593 * ,amk,amkz,amkst,gamkst
4594 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4595 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4596 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4597#if defined (ALEPH)
4598 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4599#else
4600 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4601#endif
4602 & ,names
4603 CHARACTER NAMES(NMODE)*31
4604 REAL HV(4),PT(4),PN(4),PAA(4),PMULT(4,9)
4605 real*4 pr(4),pi1(4),pi2(4),pi3(4),pi4(4),pi5(4)
4606 real*8 amp1,amp2,amp3,amp4,amp5,ams1,ams2,amom,gamom
4607 real*8 am5sq,am4sq,am3sq,am2sq,am5,am4,am3
4608 real*4 rrr(10)
4609 real*8 gg1,gg2,gg3,ff1,ff2,ff3,ff4,alp,alp1,alp2
4610#if defined (ALEPH)
4611 real*8 xm,am,gammab
4612#else
4613 real*8 xm,am,gamma
4614ccM.S.>>>>>>
4615 real*8 phspac
4616ccM.S.<<<<<<
4617#endif
4618 DATA pi /3.141592653589793238462643/
4619 DATA icont /0/
4620 data fpi /93.3e-3/
4621c
4622 COMPLEX BWIGN
4623C
4624#if defined (ALEPH)
4625 bwign(xm,am,gammab)=xm**2/cmplx(xm**2-am**2,gammab*am)
4626#else
4627 bwign(xm,am,gamma)=xm**2/cmplx(xm**2-am**2,gamma*am)
4628#endif
4629
4630C
4631 amom=.782
4632 gamom=0.0085
4633c
4634C 6 BODY PHASE SPACE NORMALISED AS IN BJORKEN-DRELL
4635C D**3 P /2E/(2PI)**3 (2PI)**4 DELTA4(SUM P)
4636 phspac=1./2**29/pi**14
4637c PHSPAC=1./2**5/PI**2
4638C init 5pi decay mode (JNPI)
4639 amp1=dcdmas(idffin(1,jnpi))
4640 amp2=dcdmas(idffin(2,jnpi))
4641 amp3=dcdmas(idffin(3,jnpi))
4642 amp4=dcdmas(idffin(4,jnpi))
4643 amp5=dcdmas(idffin(5,jnpi))
4644c
4645C TAU MOMENTUM
4646 pt(1)=0.
4647 pt(2)=0.
4648 pt(3)=0.
4649 pt(4)=amtau
4650C
4651 CALL ranmar(rrr,10)
4652C
4653c masses of 5, 4, 3 and 2 pi systems
4654c 3 pi with sampling for omega resonance
4655cam
4656c mass of 5 (12345)
4657 rr1=rrr(10)
4658 ams1=(amp1+amp2+amp3+amp4+amp5)**2
4659 ams2=(amtau-amnuta)**2
4660 am5sq=ams1+ rr1*(ams2-ams1)
4661 am5 =sqrt(am5sq)
4662 phspac=phspac*(ams2-ams1)
4663c
4664c mass of 4 (2345)
4665c flat phase space
4666 rr1=rrr(9)
4667 ams1=(amp2+amp3+amp4+amp5)**2
4668 ams2=(am5-amp1)**2
4669 am4sq=ams1+ rr1*(ams2-ams1)
4670 am4 =sqrt(am4sq)
4671 gg1=ams2-ams1
4672c
4673c mass of 3 (234)
4674C phase space with sampling for omega resonance
4675 rr1=rrr(1)
4676 ams1=(amp2+amp3+amp4)**2
4677 ams2=(am4-amp5)**2
4678 alp1=atan((ams1-amom**2)/amom/gamom)
4679 alp2=atan((ams2-amom**2)/amom/gamom)
4680 alp=alp1+rr1*(alp2-alp1)
4681 am3sq =amom**2+amom*gamom*tan(alp)
4682 am3 =sqrt(am3sq)
4683c --- this part of the jacobian will be recovered later ---------------
4684 gg2=((am3sq-amom**2)**2+(amom*gamom)**2)/(amom*gamom)
4685 gg2=gg2*(alp2-alp1)
4686c flat phase space;
4687C am3sq=ams1+ rr1*(ams2-ams1)
4688C am3 =sqrt(am3sq)
4689c --- this part of jacobian will be recovered later
4690C gg2=ams2-ams1
4691c
4692C mass of 2 (34)
4693 rr2=rrr(2)
4694 ams1=(amp3+amp4)**2
4695 ams2=(am3-amp2)**2
4696c flat phase space;
4697 am2sq=ams1+ rr2*(ams2-ams1)
4698 am2 =sqrt(am2sq)
4699c --- this part of jacobian will be recovered later
4700 gg3=ams2-ams1
4701c
4702c (34) restframe, define pi3 and pi4
4703 enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4704 enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4705 ppi= enq1**2-amp3**2
4706 pppi=sqrt(abs(enq1**2-amp3**2))
4707 ff1=(4*pi)*(2*pppi/am2)
4708c pi3 momentum in (34) rest frame
4709 call sphera(pppi,pi3)
4710 pi3(4)=enq1
4711c pi4 momentum in (34) rest frame
4712 do 30 i=1,3
4713 30 pi4(i)=-pi3(i)
4714 pi4(4)=enq2
4715c
4716c (234) rest frame, define pi2
4717c pr momentum
4718 pr(1)=0
4719 pr(2)=0
4720 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4721 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4722 ppi = pr(4)**2-am2**2
4723c pi2 momentum
4724 pi2(1)=0
4725 pi2(2)=0
4726 pi2(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4727 pi2(3)=-pr(3)
4728c --- this part of jacobian will be recovered later
4729 ff2=(4*pi)*(2*pr(3)/am3)
4730c old pions boosted from 2 rest frame to 3 rest frame
4731 exe=(pr(4)+pr(3))/am2
4732 call bostr3(exe,pi3,pi3)
4733 call bostr3(exe,pi4,pi4)
4734 rr3=rrr(3)
4735 rr4=rrr(4)
4736 thet =acos(-1.+2*rr3)
4737 phi = 2*pi*rr4
4738 call rotpol(thet,phi,pi2)
4739 call rotpol(thet,phi,pi3)
4740 call rotpol(thet,phi,pi4)
4741C
4742C (2345) rest frame, define pi5
4743c pr momentum
4744 pr(1)=0
4745 pr(2)=0
4746 pr(4)=1./(2*am4)*(am4**2+am3**2-amp5**2)
4747 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4748 ppi = pr(4)**2-am3**2
4749c pi5 momentum
4750 pi5(1)=0
4751 pi5(2)=0
4752 pi5(4)=1./(2*am4)*(am4**2-am3**2+amp5**2)
4753 pi5(3)=-pr(3)
4754c --- this part of jacobian will be recovered later
4755 ff3=(4*pi)*(2*pr(3)/am4)
4756c old pions boosted from 3 rest frame to 4 rest frame
4757 exe=(pr(4)+pr(3))/am3
4758 call bostr3(exe,pi2,pi2)
4759 call bostr3(exe,pi3,pi3)
4760 call bostr3(exe,pi4,pi4)
4761 rr3=rrr(5)
4762 rr4=rrr(6)
4763 thet =acos(-1.+2*rr3)
4764 phi = 2*pi*rr4
4765 call rotpol(thet,phi,pi2)
4766 call rotpol(thet,phi,pi3)
4767 call rotpol(thet,phi,pi4)
4768 call rotpol(thet,phi,pi5)
4769C
4770C (12345) rest frame, define pi1
4771c pr momentum
4772 pr(1)=0
4773 pr(2)=0
4774 pr(4)=1./(2*am5)*(am5**2+am4**2-amp1**2)
4775 pr(3)= sqrt(abs(pr(4)**2-am4**2))
4776 ppi = pr(4)**2-am4**2
4777c pi1 momentum
4778 pi1(1)=0
4779 pi1(2)=0
4780 pi1(4)=1./(2*am5)*(am5**2-am4**2+amp1**2)
4781 pi1(3)=-pr(3)
4782c --- this part of jacobian will be recovered later
4783 ff4=(4*pi)*(2*pr(3)/am5)
4784c old pions boosted from 4 rest frame to 5 rest frame
4785 exe=(pr(4)+pr(3))/am4
4786 call bostr3(exe,pi2,pi2)
4787 call bostr3(exe,pi3,pi3)
4788 call bostr3(exe,pi4,pi4)
4789 call bostr3(exe,pi5,pi5)
4790 rr3=rrr(7)
4791 rr4=rrr(8)
4792 thet =acos(-1.+2*rr3)
4793 phi = 2*pi*rr4
4794 call rotpol(thet,phi,pi1)
4795 call rotpol(thet,phi,pi2)
4796 call rotpol(thet,phi,pi3)
4797 call rotpol(thet,phi,pi4)
4798 call rotpol(thet,phi,pi5)
4799c
4800* now to the tau rest frame, define paa and neutrino momenta
4801* paa momentum
4802 paa(1)=0
4803 paa(2)=0
4804c paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5**2)
4805c paa(3)= sqrt(abs(paa(4)**2-am5**2))
4806c ppi = paa(4)**2-am5**2
4807 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5sq)
4808 paa(3)= sqrt(abs(paa(4)**2-am5sq))
4809 ppi = paa(4)**2-am5sq
4810 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4811* tau-neutrino momentum
4812 pn(1)=0
4813 pn(2)=0
4814 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am5**2)
4815 pn(3)=-paa(3)
4816c
4817 phspac=phspac * gg1*gg2*gg3*ff1*ff2*ff3*ff4
4818c
4819C all pions boosted from 5 rest frame to tau rest frame
4820C z-axis antiparallel to neutrino momentum
4821 exe=(paa(4)+paa(3))/am5
4822 call bostr3(exe,pi1,pi1)
4823 call bostr3(exe,pi2,pi2)
4824 call bostr3(exe,pi3,pi3)
4825 call bostr3(exe,pi4,pi4)
4826 call bostr3(exe,pi5,pi5)
4827c
4828C partial width consists of phase space and amplitude
4829C AMPLITUDE (cf YS.Tsai Phys.Rev.D4,2821(1971)
4830C or F.Gilman SH.Rhie Phys.Rev.D31,1066(1985)
4831C
4832 pxq=amtau*paa(4)
4833 pxn=amtau*pn(4)
4834 qxn=paa(4)*pn(4)-paa(1)*pn(1)-paa(2)*pn(2)-paa(3)*pn(3)
4835 brak=2*(gv**2+ga**2)*(2*pxq*qxn+am5sq*pxn)
4836 & -6*(gv**2-ga**2)*amtau*amnuta*am5sq
4837 fompp = cabs(bwign(am3,amom,gamom))**2
4838c normalisation factor (to some numerical undimensioned factor;
4839c cf R.Fischer et al ZPhys C3, 313 (1980))
4840 fnorm = 1/fpi**6
4841c AMPLIT=CCABIB**2*GFERMI**2/2. * BRAK * AM5SQ*SIGEE(AM5SQ,JNPI)
4842 amplit=ccabib**2*gfermi**2/2. * brak
4843 amplit = amplit * fompp * fnorm
4844c phase space test
4845c amplit = amplit * fnorm
4846 dgamt=1/(2.*amtau)*amplit*phspac
4847c ignore spin terms
4848 DO 40 i=1,3
4849 40 hv(i)=0.
4850c
4851 do 77 k=1,4
4852 pmult(k,1)=pi1(k)
4853 pmult(k,2)=pi2(k)
4854 pmult(k,3)=pi3(k)
4855 pmult(k,4)=pi4(k)
4856 pmult(k,5)=pi5(k)
4857 77 continue
4858 return
4859#if defined (ALEPH)
4860C missing: transposition of identical particles, statistical factors
4861C for identical matrices, polarimetric vector. Matrix element rather nai
4862#else
4863C missing: transposition of identical particles, startistical factors
4864C for identical matrices, polarimetric vector. Matrix element rather naive.
4865#endif
4866C flat phase space in pion system + with breit wigner for omega
4867C anyway it is better than nothing, and code is improvable.
4868 end
4869 SUBROUTINE dphsrk(DGAMT,HV,PN,PR,PMULT,INUM)
4870C ----------------------------------------------------------------------
4871C IT SIMULATES RHO DECAY IN TAU REST FRAME WITH
4872C Z-AXIS ALONG RHO MOMENTUM
4873C Rho decays to K Kbar
4874C ----------------------------------------------------------------------
4875 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4876 * ,ampiz,ampi,amro,gamro,ama1,gama1
4877 * ,amk,amkz,amkst,gamkst
4878C
4879 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4880 * ,ampiz,ampi,amro,gamro,ama1,gama1
4881 * ,amk,amkz,amkst,gamkst
4882 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4883 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4884 REAL HV(4),PT(4),PN(4),PR(4),PKC(4),PKZ(4),QQ(4),PMULT(4,9)
4885#if defined (ALEPH)
4886 real*4 rr1(1)
4887#else
4888 REAL RR1(1)
4889#endif
4890 DATA pi /3.141592653589793238462643/
4891 DATA icont /0/
4892C
4893C THREE BODY PHASE SPACE NORMALISED AS IN BJORKEN-DRELL
4894 phspac=1./2**11/pi**5
4895C TAU MOMENTUM
4896 pt(1)=0.
4897 pt(2)=0.
4898 pt(3)=0.
4899 pt(4)=amtau
4900C MASS OF (REAL/VIRTUAL) RHO
4901 ams1=(amk+amkz)**2
4902 ams2=(amtau-amnuta)**2
4903C FLAT PHASE SPACE
4904 CALL ranmar(rr1,1)
4905 amx2=ams1+ rr1(1)*(ams2-ams1)
4906 amx=sqrt(amx2)
4907 phspac=phspac*(ams2-ams1)
4908C PHASE SPACE WITH SAMPLING FOR RHO RESONANCE
4909c ALP1=ATAN((AMS1-AMRO**2)/AMRO/GAMRO)
4910c ALP2=ATAN((AMS2-AMRO**2)/AMRO/GAMRO)
4911CAM
4912 100 CONTINUE
4913c CALL RANMAR(RR1,1)
4914c ALP=ALP1+RR1(1)*(ALP2-ALP1)
4915c AMX2=AMRO**2+AMRO*GAMRO*TAN(ALP)
4916c AMX=SQRT(AMX2)
4917c IF(AMX.LT.(AMK+AMKZ)) GO TO 100
4918CAM
4919c PHSPAC=PHSPAC*((AMX2-AMRO**2)**2+(AMRO*GAMRO)**2)/(AMRO*GAMRO)
4920c PHSPAC=PHSPAC*(ALP2-ALP1)
4921C
4922C TAU-NEUTRINO MOMENTUM
4923 pn(1)=0
4924 pn(2)=0
4925 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
4926 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
4927C RHO MOMENTUM
4928 pr(1)=0
4929 pr(2)=0
4930 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
4931 pr(3)=-pn(3)
4932 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
4933C
4934CAM
4935 enq1=(amx2+amk**2-amkz**2)/(2.*amx)
4936 enq2=(amx2-amk**2+amkz**2)/(2.*amx)
4937 pppi=sqrt((enq1-amk)*(enq1+amk))
4938 phspac=phspac*(4*pi)*(2*pppi/amx)
4939C CHARGED PI MOMENTUM IN RHO REST FRAME
4940 CALL sphera(pppi,pkc)
4941 pkc(4)=enq1
4942C NEUTRAL PI MOMENTUM IN RHO REST FRAME
4943 DO 20 i=1,3
494420 pkz(i)=-pkc(i)
4945 pkz(4)=enq2
4946 exe=(pr(4)+pr(3))/amx
4947C PIONS BOOSTED FROM RHO REST FRAME TO TAU REST FRAME
4948 CALL bostr3(exe,pkc,pkc)
4949 CALL bostr3(exe,pkz,pkz)
4950 DO 30 i=1,4
4951 30 qq(i)=pkc(i)-pkz(i)
4952C QQ transverse to PR
4953 pksd =pr(4)*pr(4)-pr(3)*pr(3)-pr(2)*pr(2)-pr(1)*pr(1)
4954 qqpks=pr(4)* qq(4)-pr(3)* qq(3)-pr(2)* qq(2)-pr(1)* qq(1)
4955 DO 31 i=1,4
495631 qq(i)=qq(i)-pr(i)*qqpks/pksd
4957C AMPLITUDE
4958 prodpq=pt(4)*qq(4)
4959 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
4960 prodpn=pt(4)*pn(4)
4961 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
4962 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
4963 & +(gv**2-ga**2)*amtau*amnuta*qq2
4964 amplit=(gfermi*ccabib)**2*brak*2*fpirk(amx)
4965 dgamt=1/(2.*amtau)*amplit*phspac
4966 DO 40 i=1,3
4967 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
4968 do 77 k=1,4
4969 pmult(k,1)=pkc(k)
4970 pmult(k,2)=pkz(k)
4971 77 continue
4972 RETURN
4973 END
4974 FUNCTION fpirk(W)
4975C ----------------------------------------------------------
4976c square of pion form factor
4977C ----------------------------------------------------------
4978 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4979 * ,ampiz,ampi,amro,gamro,ama1,gama1
4980 * ,amk,amkz,amkst,gamkst
4981C
4982 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4983 * ,ampiz,ampi,amro,gamro,ama1,gama1
4984 * ,amk,amkz,amkst,gamkst
4985c COMPLEX FPIKMK
4986 COMPLEX FPIKM
4987 fpirk=cabs(fpikm(w,amk,amkz))**2
4988c FPIRK=CABS(FPIKMK(W,AMK,AMKZ))**2
4989 END
4990 COMPLEX FUNCTION fpikmk(W,XM1,XM2)
4991C **********************************************************
4992C Kaon form factor
4993C **********************************************************
4994 COMPLEX BWIGM
4995 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
4996 EXTERNAL bwig
4997 DATA init /0/
4998C
4999C ------------ PARAMETERS --------------------
5000 IF (init.EQ.0 ) THEN
5001 init=1
5002 pi=3.141592654
5003 pim=.140
5004 rom=0.773
5005 rog=0.145
5006 rom1=1.570
5007 rog1=0.510
5008c BETA1=-0.111
5009 beta1=-0.221
5010 ENDIF
5011C -----------------------------------------------
5012 s=w**2
5013 fpikmk=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
5014 & /(1+beta1)
5015 RETURN
5016 END
5017 SUBROUTINE reslux
5018C ****************
5019C INITIALIZE LUND COMMON
5020#if defined (CLEO)
5021#else
5022 parameter(nmxhep=2000)
5023 common/hepevtx/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
5024 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
5025 SAVE /hepevtx/
5026#endif
5027 nhep=0
5028 END
5029 SUBROUTINE dwrph(KTO,PHX)
5030C
5031C -------------------------
5032C
5033 IMPLICIT REAL*8 (a-h,o-z)
5034 real*4 phx(4)
5035 real*4 qhot(4)
5036C
5037 DO 9 k=1,4
5038 qhot(k) =0.0
5039 9 CONTINUE
5040C CASE OF TAU RADIATIVE DECAYS.
5041C FILLING OF THE LUND COMMON BLOCK.
5042 DO 1002 i=1,4
5043 1002 qhot(i)=phx(i)
5044 IF (qhot(4).GT.1.e-5) CALL dwluph(kto,qhot)
5045 RETURN
5046 END
5047 SUBROUTINE dwluph(KTO,PHOT)
5048C---------------------------------------------------------------------
5049C Lorentz transformation to CMsystem and
5050C Updating of HEPEVT record
5051C
5052C called by : DEXAY1,(DEKAY1,DEKAY2)
5053C
5054C used when radiative corrections in decays are generated
5055C---------------------------------------------------------------------
5056C
5057#if defined (ALEPH)
5058 COMMON /taupos/ np1,np2
5059#else
5060#endif
5061 REAL PHOT(4)
5062#if defined (ALEPH)
5063#else
5064 COMMON /taupos/ np1,np2
5065#endif
5066C
5067C check energy
5068 IF (phot(4).LE.0.0) RETURN
5069C
5070C position of decaying particle:
5071 IF((kto.EQ. 1).OR.(kto.EQ.11)) THEN
5072 nps=np1
5073 ELSE
5074 nps=np2
5075 ENDIF
5076C
5077 ktos=kto
5078 IF(ktos.GT.10) ktos=ktos-10
5079C boost and append photon (gamma is 22)
5080 CALL tralo4(ktos,phot,phot,am)
5081 CALL filhep(0,1,22,nps,nps,0,0,phot,0.0,.true.)
5082C
5083 RETURN
5084 END
5085
5086 SUBROUTINE dwluel(KTO,ISGN,PNU,PWB,PEL,PNE)
5087C ----------------------------------------------------------------------
5088C Lorentz transformation to CMsystem and
5089C Updating of HEPEVT record
5090C
5091C ISGN = 1/-1 for tau-/tau+
5092C
5093C called by : DEXAY,(DEKAY1,DEKAY2)
5094C ----------------------------------------------------------------------
5095C
5096#if defined (ALEPH)
5097 COMMON /taupos/ np1,np2
5098#else
5099#endif
5100 REAL PNU(4),PWB(4),PEL(4),PNE(4)
5101#if defined (ALEPH)
5102#else
5103 COMMON /taupos/ np1,np2
5104#endif
5105C
5106C position of decaying particle:
5107 IF(kto.EQ. 1) THEN
5108 nps=np1
5109 ELSE
5110 nps=np2
5111 ENDIF
5112C
5113C tau neutrino (nu_tau is 16)
5114 CALL tralo4(kto,pnu,pnu,am)
5115 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5116C
5117C W boson (W+ is 24)
5118 CALL tralo4(kto,pwb,pwb,am)
5119C CALL FILHEP(0,2,-24*ISGN,NPS,NPS,0,0,PWB,AM,.TRUE.)
5120C
5121C electron (e- is 11)
5122 CALL tralo4(kto,pel,pel,am)
5123 CALL filhep(0,1,11*isgn,nps,nps,0,0,pel,am,.false.)
5124C
5125C anti electron neutrino (nu_e is 12)
5126 CALL tralo4(kto,pne,pne,am)
5127 CALL filhep(0,1,-12*isgn,nps,nps,0,0,pne,am,.true.)
5128C
5129 RETURN
5130 END
5131 SUBROUTINE dwlumu(KTO,ISGN,PNU,PWB,PMU,PNM)
5132C ----------------------------------------------------------------------
5133C Lorentz transformation to CMsystem and
5134C Updating of HEPEVT record
5135C
5136C ISGN = 1/-1 for tau-/tau+
5137C
5138C called by : DEXAY,(DEKAY1,DEKAY2)
5139C ----------------------------------------------------------------------
5140C
5141#if defined (ALEPH)
5142 COMMON /taupos/ np1,np2
5143#else
5144#endif
5145 REAL PNU(4),PWB(4),PMU(4),PNM(4)
5146#if defined (ALEPH)
5147#else
5148 COMMON /taupos/ np1,np2
5149#endif
5150C
5151C position of decaying particle:
5152 IF(kto.EQ. 1) THEN
5153 nps=np1
5154 ELSE
5155 nps=np2
5156 ENDIF
5157C
5158C tau neutrino (nu_tau is 16)
5159 CALL tralo4(kto,pnu,pnu,am)
5160 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5161C
5162C W boson (W+ is 24)
5163 CALL tralo4(kto,pwb,pwb,am)
5164C CALL FILHEP(0,2,-24*ISGN,NPS,NPS,0,0,PWB,AM,.TRUE.)
5165C
5166C muon (mu- is 13)
5167 CALL tralo4(kto,pmu,pmu,am)
5168 CALL filhep(0,1,13*isgn,nps,nps,0,0,pmu,am,.false.)
5169C
5170C anti muon neutrino (nu_mu is 14)
5171 CALL tralo4(kto,pnm,pnm,am)
5172 CALL filhep(0,1,-14*isgn,nps,nps,0,0,pnm,am,.true.)
5173C
5174 RETURN
5175 END
5176 SUBROUTINE dwlupi(KTO,ISGN,PPI,PNU)
5177C ----------------------------------------------------------------------
5178C Lorentz transformation to CMsystem and
5179C Updating of HEPEVT record
5180C
5181C ISGN = 1/-1 for tau-/tau+
5182C
5183C called by : DEXAY,(DEKAY1,DEKAY2)
5184C ----------------------------------------------------------------------
5185C
5186 REAL PNU(4),PPI(4)
5187 COMMON /taupos/ np1,np2
5188C
5189C position of decaying particle:
5190 IF(kto.EQ. 1) THEN
5191 nps=np1
5192 ELSE
5193 nps=np2
5194 ENDIF
5195C
5196C tau neutrino (nu_tau is 16)
5197 CALL tralo4(kto,pnu,pnu,am)
5198 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5199C
5200C charged pi meson (pi+ is 211)
5201 CALL tralo4(kto,ppi,ppi,am)
5202 CALL filhep(0,1,-211*isgn,nps,nps,0,0,ppi,am,.true.)
5203C
5204 RETURN
5205 END
5206 SUBROUTINE dwluro(KTO,ISGN,PNU,PRHO,PIC,PIZ)
5207C ----------------------------------------------------------------------
5208C Lorentz transformation to CMsystem and
5209C Updating of HEPEVT record
5210C
5211C ISGN = 1/-1 for tau-/tau+
5212C
5213C called by : DEXAY,(DEKAY1,DEKAY2)
5214C ----------------------------------------------------------------------
5215C
5216#if defined (ALEPH)
5217 COMMON /taupos/ np1,np2
5218#else
5219#endif
5220 REAL PNU(4),PRHO(4),PIC(4),PIZ(4)
5221#if defined (ALEPH)
5222#else
5223 COMMON /taupos/ np1,np2
5224#endif
5225C
5226C position of decaying particle:
5227 IF(kto.EQ. 1) THEN
5228 nps=np1
5229 ELSE
5230 nps=np2
5231 ENDIF
5232C
5233C tau neutrino (nu_tau is 16)
5234 CALL tralo4(kto,pnu,pnu,am)
5235 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5236C
5237C charged rho meson (rho+ is 213)
5238 CALL tralo4(kto,prho,prho,am)
5239 CALL filhep(0,2,-213*isgn,nps,nps,0,0,prho,am,.true.)
5240C
5241C charged pi meson (pi+ is 211)
5242 CALL tralo4(kto,pic,pic,am)
5243 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pic,am,.true.)
5244C
5245C pi0 meson (pi0 is 111)
5246 CALL tralo4(kto,piz,piz,am)
5247 CALL filhep(0,1,111,-2,-2,0,0,piz,am,.true.)
5248C
5249 RETURN
5250 END
5251 SUBROUTINE dwluaa(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
5252C ----------------------------------------------------------------------
5253C Lorentz transformation to CMsystem and
5254C Updating of HEPEVT record
5255C
5256C ISGN = 1/-1 for tau-/tau+
5257C JAA = 1 (2) FOR A_1- DECAY TO PI+ 2PI- (PI- 2PI0)
5258C
5259C called by : DEXAY,(DEKAY1,DEKAY2)
5260C ----------------------------------------------------------------------
5261C
5262#if defined (ALEPH)
5263 COMMON /taupos/ np1,np2
5264#else
5265#endif
5266 REAL PNU(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
5267#if defined (ALEPH)
5268#else
5269 COMMON /taupos/ np1,np2
5270#endif
5271C
5272C position of decaying particle:
5273 IF(kto.EQ. 1) THEN
5274 nps=np1
5275 ELSE
5276 nps=np2
5277 ENDIF
5278C
5279C tau neutrino (nu_tau is 16)
5280 CALL tralo4(kto,pnu,pnu,am)
5281 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5282C
5283C charged a_1 meson (a_1+ is 20213)
5284 CALL tralo4(kto,paa,paa,am)
5285 CALL filhep(0,1,-20213*isgn,nps,nps,0,0,paa,am,.true.)
5286C
5287C two possible decays of the charged a1 meson
5288 IF(jaa.EQ.1) THEN
5289C
5290C A1 --> PI+ PI- PI- (or charged conjugate)
5291C
5292C pi minus (or c.c.) (pi+ is 211)
5293 CALL tralo4(kto,pim2,pim2,am)
5294 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pim2,am,.true.)
5295C
5296C pi minus (or c.c.) (pi+ is 211)
5297 CALL tralo4(kto,pim1,pim1,am)
5298 CALL filhep(0,1,-211*isgn,-2,-2,0,0,pim1,am,.true.)
5299C
5300C pi plus (or c.c.) (pi+ is 211)
5301 CALL tralo4(kto,pipl,pipl,am)
5302 CALL filhep(0,1, 211*isgn,-3,-3,0,0,pipl,am,.true.)
5303C
5304 ELSE IF (jaa.EQ.2) THEN
5305C
5306C A1 --> PI- PI0 PI0 (or charged conjugate)
5307C
5308C pi zero (pi0 is 111)
5309 CALL tralo4(kto,pim2,pim2,am)
5310 CALL filhep(0,1,111,-1,-1,0,0,pim2,am,.true.)
5311C
5312C pi zero (pi0 is 111)
5313 CALL tralo4(kto,pim1,pim1,am)
5314 CALL filhep(0,1,111,-2,-2,0,0,pim1,am,.true.)
5315C
5316C pi minus (or c.c.) (pi+ is 211)
5317 CALL tralo4(kto,pipl,pipl,am)
5318 CALL filhep(0,1,-211*isgn,-3,-3,0,0,pipl,am,.true.)
5319C
5320 ENDIF
5321C
5322 RETURN
5323 END
5324 SUBROUTINE dwlukk (KTO,ISGN,PKK,PNU)
5325C ----------------------------------------------------------------------
5326C Lorentz transformation to CMsystem and
5327C Updating of HEPEVT record
5328C
5329C ISGN = 1/-1 for tau-/tau+
5330C
5331C ----------------------------------------------------------------------
5332C
5333 REAL PKK(4),PNU(4)
5334 COMMON /taupos/ np1,np2
5335C
5336C position of decaying particle
5337#if defined (ALEPH)
5338 IF(kto.EQ. 1) THEN
5339#else
5340 IF (kto.EQ.1) THEN
5341#endif
5342 nps=np1
5343 ELSE
5344 nps=np2
5345 ENDIF
5346C
5347C tau neutrino (nu_tau is 16)
5348 CALL tralo4 (kto,pnu,pnu,am)
5349 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5350C
5351C K meson (K+ is 321)
5352 CALL tralo4 (kto,pkk,pkk,am)
5353 CALL filhep(0,1,-321*isgn,nps,nps,0,0,pkk,am,.true.)
5354C
5355 RETURN
5356 END
5357 SUBROUTINE dwluks(KTO,ISGN,PNU,PKS,PKK,PPI,JKST)
5358 COMMON / taukle / bra1,brk0,brk0b,brks
5359 real*4 bra1,brk0,brk0b,brks
5360#if defined (ALEPH)
5361 COMMON /taupos/ np1,np2
5362 real*4 xio(1)
5363#endif
5364C ----------------------------------------------------------------------
5365C Lorentz transformation to CMsystem and
5366C Updating of HEPEVT record
5367C
5368C ISGN = 1/-1 for tau-/tau+
5369C JKST=10 (20) corresponds to K0B pi- (K- pi0) decay
5370C
5371C ----------------------------------------------------------------------
5372C
5373#if defined (ALEPH)
5374 REAL PNU(4),PKS(4),PKK(4),PPI(4)
5375#else
5376 REAL PNU(4),PKS(4),PKK(4),PPI(4),XIO(1)
5377 COMMON /taupos/ np1,np2
5378#endif
5379C
5380C position of decaying particle
5381 IF(kto.EQ. 1) THEN
5382 nps=np1
5383 ELSE
5384 nps=np2
5385 ENDIF
5386C
5387C tau neutrino (nu_tau is 16)
5388 CALL tralo4(kto,pnu,pnu,am)
5389 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5390C
5391C charged K* meson (K*+ is 323)
5392 CALL tralo4(kto,pks,pks,am)
5393 CALL filhep(0,1,-323*isgn,nps,nps,0,0,pks,am,.true.)
5394C
5395C two possible decay modes of charged K*
5396 IF(jkst.EQ.10) THEN
5397C
5398C K*- --> pi- K0B (or charged conjugate)
5399C
5400C charged pi meson (pi+ is 211)
5401 CALL tralo4(kto,ppi,ppi,am)
5402 CALL filhep(0,1,-211*isgn,-1,-1,0,0,ppi,am,.true.)
5403C
5404 bran=brk0b
5405 IF (isgn.EQ.-1) bran=brk0
5406C K0 --> K0_long (is 130) / K0_short (is 310) = 1/1
5407 CALL ranmar(xio,1)
5408 IF(xio(1).GT.bran) THEN
5409 k0type = 130
5410 ELSE
5411 k0type = 310
5412 ENDIF
5413C
5414 CALL tralo4(kto,pkk,pkk,am)
5415 CALL filhep(0,1,k0type,-2,-2,0,0,pkk,am,.true.)
5416C
5417 ELSE IF(jkst.EQ.20) THEN
5418C
5419C K*- --> pi0 K-
5420C
5421C pi zero (pi0 is 111)
5422 CALL tralo4(kto,ppi,ppi,am)
5423 CALL filhep(0,1,111,-1,-1,0,0,ppi,am,.true.)
5424C
5425C charged K meson (K+ is 321)
5426 CALL tralo4(kto,pkk,pkk,am)
5427 CALL filhep(0,1,-321*isgn,-2,-2,0,0,pkk,am,.true.)
5428C
5429 ENDIF
5430C
5431 RETURN
5432 END
5433 SUBROUTINE dwlnew(KTO,ISGN,PNU,PWB,PNPI,MODE)
5434C ----------------------------------------------------------------------
5435C Lorentz transformation to CMsystem and
5436C Updating of HEPEVT record
5437C
5438C ISGN = 1/-1 for tau-/tau+
5439C
5440C called by : DEXAY,(DEKAY1,DEKAY2)
5441C ----------------------------------------------------------------------
5442C
5443 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
5444#if defined (ALEPH)
5445 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5446#else
5447 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5448#endif
5449 & ,names
5450 COMMON /taupos/ np1,np2
5451 CHARACTER NAMES(NMODE)*31
5452 REAL PNU(4),PWB(4),PNPI(4,9)
5453 REAL PPI(4)
5454C
5455 jnpi=mode-7
5456C position of decaying particle
5457 IF(kto.EQ. 1) THEN
5458 nps=np1
5459 ELSE
5460 nps=np2
5461 ENDIF
5462C
5463C tau neutrino (nu_tau is 16)
5464 CALL tralo4(kto,pnu,pnu,am)
5465 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5466C
5467C W boson (W+ is 24)
5468 CALL tralo4(kto,pwb,pwb,am)
5469 CALL filhep(0,1,-24*isgn,nps,nps,0,0,pwb,am,.true.)
5470C
5471C multi pi mode JNPI
5472C
5473C get multiplicity of mode JNPI
5474 nd=mulpik(jnpi)
5475 DO i=1,nd
5476#if defined (ALEPH)
5477cam KFPI=LUNPIK(IDFFIN(I,JNPI),-ISGN)
5478 kfpi=lunpik(idffin(i,jnpi), isgn)
5479#else
5480 kfpi=lunpik(idffin(i,jnpi),-isgn)
5481#endif
5482C for charged conjugate case, change charged pions only
5483C IF(KFPI.NE.111)KFPI=KFPI*ISGN
5484 DO j=1,4
5485 ppi(j)=pnpi(j,i)
5486 END DO
5487 CALL tralo4(kto,ppi,ppi,am)
5488 CALL filhep(0,1,kfpi,-i,-i,0,0,ppi,am,.true.)
5489 END DO
5490C
5491 RETURN
5492 END
5493#if defined (CePeCe)
5494#else
5495#endif
5496 FUNCTION amast(PP)
5497C ----------------------------------------------------------------------
5498C CALCULATES MASS OF PP (DOUBLE PRECISION)
5499C
5500C USED BY : RADKOR
5501C ----------------------------------------------------------------------
5502 IMPLICIT REAL*8 (a-h,o-z)
5503 real*8 pp(4)
5504 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5505C
5506 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5507 amast=aaa
5508 RETURN
5509 END
5510 FUNCTION amas4(PP)
5511C ******************
5512C ----------------------------------------------------------------------
5513C CALCULATES MASS OF PP
5514C
5515C USED BY :
5516C ----------------------------------------------------------------------
5517 REAL PP(4)
5518 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5519 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5520 amas4=aaa
5521 RETURN
5522 END
5523 FUNCTION angxy(X,Y)
5524C ----------------------------------------------------------------------
5525C
5526C USED BY : KORALZ RADKOR
5527C ----------------------------------------------------------------------
5528 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5529 DATA pi /3.141592653589793238462643d0/
5530C
5531 IF(abs(y).LT.abs(x)) THEN
5532 the=atan(abs(y/x))
5533 IF(x.LE.0d0) the=pi-the
5534 ELSE
5535 the=acos(x/sqrt(x**2+y**2))
5536 ENDIF
5537 angxy=the
5538 RETURN
5539 END
5540 FUNCTION angfi(X,Y)
5541C ----------------------------------------------------------------------
5542* CALCULATES ANGLE IN (0,2*PI) RANGE OUT OF X-Y
5543C
5544C USED BY : KORALZ RADKOR
5545C ----------------------------------------------------------------------
5546 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5547 DATA pi /3.141592653589793238462643d0/
5548C
5549 IF(abs(y).LT.abs(x)) THEN
5550 the=atan(abs(y/x))
5551 IF(x.LE.0d0) the=pi-the
5552 ELSE
5553 the=acos(x/sqrt(x**2+y**2))
5554 ENDIF
5555 IF(y.LT.0d0) the=2d0*pi-the
5556 angfi=the
5557 END
5558 SUBROUTINE rotod1(PH1,PVEC,QVEC)
5559C ----------------------------------------------------------------------
5560C
5561C USED BY : KORALZ
5562C ----------------------------------------------------------------------
5563 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5564 dimension pvec(4),qvec(4),rvec(4)
5565C
5566 phi=ph1
5567 cs=cos(phi)
5568 sn=sin(phi)
5569 DO 10 i=1,4
5570 10 rvec(i)=pvec(i)
5571 qvec(1)=rvec(1)
5572 qvec(2)= cs*rvec(2)-sn*rvec(3)
5573 qvec(3)= sn*rvec(2)+cs*rvec(3)
5574 qvec(4)=rvec(4)
5575 RETURN
5576 END
5577 SUBROUTINE rotod2(PH1,PVEC,QVEC)
5578C ----------------------------------------------------------------------
5579C
5580C USED BY : KORALZ RADKOR
5581C ----------------------------------------------------------------------
5582 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5583 dimension pvec(4),qvec(4),rvec(4)
5584C
5585 phi=ph1
5586 cs=cos(phi)
5587 sn=sin(phi)
5588 DO 10 i=1,4
5589 10 rvec(i)=pvec(i)
5590 qvec(1)= cs*rvec(1)+sn*rvec(3)
5591 qvec(2)=rvec(2)
5592 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5593 qvec(4)=rvec(4)
5594 RETURN
5595 END
5596 SUBROUTINE rotod3(PH1,PVEC,QVEC)
5597C ----------------------------------------------------------------------
5598C
5599C USED BY : KORALZ RADKOR
5600C ----------------------------------------------------------------------
5601 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5602C
5603 dimension pvec(4),qvec(4),rvec(4)
5604 phi=ph1
5605 cs=cos(phi)
5606 sn=sin(phi)
5607 DO 10 i=1,4
5608 10 rvec(i)=pvec(i)
5609 qvec(1)= cs*rvec(1)-sn*rvec(2)
5610 qvec(2)= sn*rvec(1)+cs*rvec(2)
5611 qvec(3)=rvec(3)
5612 qvec(4)=rvec(4)
5613 END
5614 SUBROUTINE bostr3(EXE,PVEC,QVEC)
5615C ----------------------------------------------------------------------
5616C BOOST ALONG Z AXIS, EXE=EXP(ETA), ETA= HIPERBOLIC VELOCITY.
5617C
5618C USED BY : TAUOLA KORALZ (?)
5619C ----------------------------------------------------------------------
5620 real*4 pvec(4),qvec(4),rvec(4)
5621C
5622 DO 10 i=1,4
5623 10 rvec(i)=pvec(i)
5624 rpl=rvec(4)+rvec(3)
5625 rmi=rvec(4)-rvec(3)
5626 qpl=rpl*exe
5627 qmi=rmi/exe
5628 qvec(1)=rvec(1)
5629 qvec(2)=rvec(2)
5630 qvec(3)=(qpl-qmi)/2
5631 qvec(4)=(qpl+qmi)/2
5632 END
5633 SUBROUTINE bostd3(EXE,PVEC,QVEC)
5634C ----------------------------------------------------------------------
5635C BOOST ALONG Z AXIS, EXE=EXP(ETA), ETA= HIPERBOLIC VELOCITY.
5636C
5637C USED BY : KORALZ RADKOR
5638C ----------------------------------------------------------------------
5639 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5640 dimension pvec(4),qvec(4),rvec(4)
5641C
5642 DO 10 i=1,4
5643 10 rvec(i)=pvec(i)
5644 rpl=rvec(4)+rvec(3)
5645 rmi=rvec(4)-rvec(3)
5646 qpl=rpl*exe
5647 qmi=rmi/exe
5648 qvec(1)=rvec(1)
5649 qvec(2)=rvec(2)
5650 qvec(3)=(qpl-qmi)/2
5651 qvec(4)=(qpl+qmi)/2
5652 RETURN
5653 END
5654 SUBROUTINE rotor1(PH1,PVEC,QVEC)
5655C ----------------------------------------------------------------------
5656C
5657C called by :
5658C ----------------------------------------------------------------------
5659 real*4 pvec(4),qvec(4),rvec(4)
5660C
5661 phi=ph1
5662 cs=cos(phi)
5663 sn=sin(phi)
5664 DO 10 i=1,4
5665 10 rvec(i)=pvec(i)
5666 qvec(1)=rvec(1)
5667 qvec(2)= cs*rvec(2)-sn*rvec(3)
5668 qvec(3)= sn*rvec(2)+cs*rvec(3)
5669 qvec(4)=rvec(4)
5670 END
5671 SUBROUTINE rotor2(PH1,PVEC,QVEC)
5672C ----------------------------------------------------------------------
5673C
5674C USED BY : TAUOLA
5675C ----------------------------------------------------------------------
5676 IMPLICIT REAL*4(a-h,o-z)
5677 real*4 pvec(4),qvec(4),rvec(4)
5678C
5679 phi=ph1
5680 cs=cos(phi)
5681 sn=sin(phi)
5682 DO 10 i=1,4
5683 10 rvec(i)=pvec(i)
5684 qvec(1)= cs*rvec(1)+sn*rvec(3)
5685 qvec(2)=rvec(2)
5686 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5687 qvec(4)=rvec(4)
5688 END
5689 SUBROUTINE rotor3(PHI,PVEC,QVEC)
5690C ----------------------------------------------------------------------
5691C
5692C USED BY : TAUOLA
5693C ----------------------------------------------------------------------
5694 real*4 pvec(4),qvec(4),rvec(4)
5695C
5696 cs=cos(phi)
5697 sn=sin(phi)
5698 DO 10 i=1,4
5699 10 rvec(i)=pvec(i)
5700 qvec(1)= cs*rvec(1)-sn*rvec(2)
5701 qvec(2)= sn*rvec(1)+cs*rvec(2)
5702 qvec(3)=rvec(3)
5703 qvec(4)=rvec(4)
5704 END
5705 SUBROUTINE spherd(R,X)
5706C ----------------------------------------------------------------------
5707C GENERATES UNIFORMLY THREE-VECTOR X ON SPHERE OF RADIUS R
5708C DOUBLE PRECISON VERSION OF SPHERA
5709C ----------------------------------------------------------------------
5710 real*8 r,x(4),pi,costh,sinth
5711 real*4 rrr(2)
5712 DATA pi /3.141592653589793238462643d0/
5713C
5714 CALL ranmar(rrr,2)
5715 costh=-1+2*rrr(1)
5716 sinth=sqrt(1 -costh**2)
5717 x(1)=r*sinth*cos(2*pi*rrr(2))
5718 x(2)=r*sinth*sin(2*pi*rrr(2))
5719 x(3)=r*costh
5720 RETURN
5721 END
5722 SUBROUTINE rotpox(THET,PHI,PP)
5723 IMPLICIT REAL*8 (a-h,o-z)
5724C ----------------------------------------------------------------------
5725#if defined (ALEPH)
5726C double precison version of ROTPOL
5727#else
5728C
5729#endif
5730C ----------------------------------------------------------------------
5731 dimension pp(4)
5732C
5733 CALL rotod2(thet,pp,pp)
5734 CALL rotod3( phi,pp,pp)
5735 RETURN
5736 END
5737 SUBROUTINE sphera(R,X)
5738C ----------------------------------------------------------------------
5739C GENERATES UNIFORMLY THREE-VECTOR X ON SPHERE OF RADIUS R
5740C
5741C called by : DPHSxx,DADMPI,DADMKK
5742C ----------------------------------------------------------------------
5743 REAL X(4)
5744 real*4 rrr(2)
5745 DATA pi /3.141592653589793238462643/
5746C
5747 CALL ranmar(rrr,2)
5748 costh=-1.+2.*rrr(1)
5749 sinth=sqrt(1.-costh**2)
5750 x(1)=r*sinth*cos(2*pi*rrr(2))
5751 x(2)=r*sinth*sin(2*pi*rrr(2))
5752 x(3)=r*costh
5753 RETURN
5754 END
5755 SUBROUTINE rotpol(THET,PHI,PP)
5756C ----------------------------------------------------------------------
5757C
5758C called by : DADMAA,DPHSAA
5759C ----------------------------------------------------------------------
5760 REAL PP(4)
5761C
5762 CALL rotor2(thet,pp,pp)
5763 CALL rotor3( phi,pp,pp)
5764 RETURN
5765 END
5766#include "../randg/tauola-random.h"
5767 FUNCTION dilogt(X)
5768C *****************
5769 IMPLICIT REAL*8(a-h,o-z)
5770CERN C304 VERSION 29/07/71 DILOG 59 C
5771 z=-1.64493406684822
5772 IF(x .LT.-1.0) GO TO 1
5773 IF(x .LE. 0.5) GO TO 2
5774 IF(x .EQ. 1.0) GO TO 3
5775 IF(x .LE. 2.0) GO TO 4
5776 z=3.2898681336964
5777 1 t=1.0/x
5778 s=-0.5
5779 z=z-0.5* log(abs(x))**2
5780 GO TO 5
5781 2 t=x
5782 s=0.5
5783 z=0.
5784 GO TO 5
5785 3 dilogt=1.64493406684822
5786 RETURN
5787 4 t=1.0-x
5788 s=-0.5
5789 z=1.64493406684822 - log(x)* log(abs(t))
5790 5 y=2.66666666666666 *t+0.66666666666666
5791 b= 0.00000 00000 00001
5792 a=y*b +0.00000 00000 00004
5793 b=y*a-b+0.00000 00000 00011
5794 a=y*b-a+0.00000 00000 00037
5795 b=y*a-b+0.00000 00000 00121
5796 a=y*b-a+0.00000 00000 00398
5797 b=y*a-b+0.00000 00000 01312
5798 a=y*b-a+0.00000 00000 04342
5799 b=y*a-b+0.00000 00000 14437
5800 a=y*b-a+0.00000 00000 48274
5801 b=y*a-b+0.00000 00001 62421
5802 a=y*b-a+0.00000 00005 50291
5803 b=y*a-b+0.00000 00018 79117
5804 a=y*b-a+0.00000 00064 74338
5805 b=y*a-b+0.00000 00225 36705
5806 a=y*b-a+0.00000 00793 87055
5807 b=y*a-b+0.00000 02835 75385
5808 a=y*b-a+0.00000 10299 04264
5809 b=y*a-b+0.00000 38163 29463
5810 a=y*b-a+0.00001 44963 00557
5811 b=y*a-b+0.00005 68178 22718
5812 a=y*b-a+0.00023 20021 96094
5813 b=y*a-b+0.00100 16274 96164
5814 a=y*b-a+0.00468 63619 59447
5815 b=y*a-b+0.02487 93229 24228
5816 a=y*b-a+0.16607 30329 27855
5817 a=y*a-b+1.93506 43008 6996
5818 dilogt=s*t*(a-b)+z
5819 RETURN
5820C=======================================================================
5821C===================END OF CPC PART ====================================
5822C=======================================================================
5823 END
5824
5825C-----------------------------------------------------------------------
5826C Initialize RChL Currents
5827C Dummy routine for compatibility with new updates of TAUOLA
5828C
5829C Added 25.Jul.2012
5830C-----------------------------------------------------------------------
5831 SUBROUTINE inirchl(FLAG)
5832 INTEGER FLAG
5833
5834 IF(flag.NE.0) THEN
5835 WRITE(*,25) flag
5836 25 FORMAT(1x, "TAUOLA IniRChL: Fatal error, FLAG=",i2," but RChL currents missing")
5837 WRITE(*,*) " in loaded version of TAUOLA-FORTRAN library."
5838 stop
5839 ENDIF
5840
5841 RETURN
5842 END