2 c=============================================================
5 c *********************
7 c **********************************************************************
10 c *********tauola library: version 2.7 ******** *
11 c **************december 1993****************** *
13 c *********tauola library: version 2.7 ******** *
14 c **************august 1995****************** *
16 c ** authors: s.jadach, z.was ***** *
17 c ** r. decker, m. jezabek, j.h.kuehn, ***** *
18 c ********available from: wasm at cernvm ****** *
19 c *******published in comp. phys. comm.******** *
20 c *** preprint cern-th-5856 september 1990 **** *
21 c *** preprint cern-th-6195 october 1991 **** *
22 c *** preprint cern-th-6793 november 1992 **** *
23 c **********************************************************************
25 c ----------------------------------------------------------------------
27 c chooses decay mode according to list of branching ratios
37 c jak=14-19 kkpi & kpipi modes
38 c jak=20-21 eta pi pi; gamma pi pi modes
44 c ----------------------------------------------------------------------
45 COMMON / taubra / gamprt(30),jlist(30),nchan
52 IF(nchan.LE.0.OR.nchan.GT.30) goto 902
59 IF(rrr(1).LT.cumul(i)/cumul(nchan)) ji=i
64 9020
FORMAT(
' ----- JAKER: WRONG NCHAN')
67 SUBROUTINE dekay(KTO,HX)
68 c ***********************
69 c this dekay is in spirit of the
'DECAY' which
70 c was included in koral-b program, comp. phys. commun.
71 c vol. 36 (1985) 191, see comments on general philosophy there.
72 c kto=0 initialisation(obligatory)
73 c kto=1,11 denotes tau+ and kto=2,12 tau-
74 c dekay(1,h) and dekay(2,h) is called internally by mc generator.
75 c h denotes the polarimetric vector, used by the host
PROGRAM for
76 c calculation of the spin weight.
77 c user may optionally CALL dekay(11,h) dekay(12,h) in order
78 c to transform decay products to cms and
WRITE lund record in /lujets/.
79 c kto=100, print final report(
OPTIONAL).
81 c jak=1 electron decay
90 c jak=14-19 kkpi & kpipi modes
91 c jak=20-21 eta pi pi; gamma pi pi modes
92 c jak=0 inclusive: jak=1-21
95 c jak=0 inclusive: jak=1,2,3,4,5,6,7,8
99 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
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)
110 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
112 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
115 CHARACTER names(nmode)*31
116 COMMON / inout / inut,iout
117 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4),hdum(4)
126 c initialisation or reinitialisation
127 c first or second tau positions in hepevt as in koralb/z
131 IF (iwarm.EQ.1) x=5/(iwarm-1)
133 WRITE(iout,7001) jak1,jak2
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)
151 ELSEIF(kto.EQ.1)
THEN
152 c =====================
153 c decay of tau+ in the tau rest frame
155 IF(iwarm.EQ.0) goto 902
158 #elif defined (ALEPH)
160 c ajwmod to change brs depending on sign:
163 CALL dekay1(0,h,isgn)
164 ELSEIF(kto.EQ.2)
THEN
165 c =================================
166 c decay of tau- in the tau rest frame
168 IF(iwarm.EQ.0) goto 902
171 #elif defined (ALEPH)
173 c ajwmod to change brs depending on sign:
176 CALL dekay2(0,h,isgn)
177 ELSEIF(kto.EQ.11)
THEN
178 c ======================
179 c rest of decay procedure for accepted tau+ decay
182 CALL dekay1(1,h,isgn)
183 ELSEIF(kto.EQ.12)
THEN
184 c ======================
185 c rest of decay procedure for accepted tau- decay
188 CALL dekay2(1,h,isgn)
189 ELSEIF(kto.EQ.100)
THEN
190 c =======================
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)
203 $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
214 7001
FORMAT(///1x,15(5h*****)
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*,
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*,
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*,
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*,
255 7010
FORMAT(///1x,15(5h*****)
257 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
258 $ /,
' *', 25x,
'***********DECEMBER 1993***************',9x,1h*,
260 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
261 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
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*,
275 $
' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
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*)
285 $ ,i10,2f12.7,a31 ,8x,1h*)
287 $ ,20x,
'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
288 $ /,
' *',20x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
291 9020
FORMAT(
' ----- DEKAY: LACK OF INITIALISATION')
294 9100
FORMAT(
' ----- DEKAY: WRONG VALUE OF KTO ')
297 SUBROUTINE dekay1(IMOD,HH,ISGN)
298 c *******************************
299 c this routine simulates tau+ decay
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
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
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)
324 IF(jak1.EQ.-1)
RETURN
329 IF(jak1.EQ.0) CALL jaker(jak)
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)
345 CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
351 ELSEIF(imd.EQ.1)
THEN
352 c =====================
355 nevdec(jak)=nevdec(jak)+1
360 CALL dwluel(1,isgn,pnu,pwb,pmu,pnm)
361 CALL dwrph(ktom,phot)
365 ELSEIF(jak.EQ.2)
THEN
366 CALL dwlumu(1,isgn,pnu,pwb,pmu,pnm)
367 CALL dwrph(ktom,phot)
371 ELSEIF(jak.EQ.3)
THEN
372 CALL dwlupi(1,isgn,ppi,pnu)
376 ELSEIF(jak.EQ.4)
THEN
377 CALL dwluro(1,isgn,pnu,prho,pic,piz)
381 ELSEIF(jak.EQ.5)
THEN
382 CALL dwluaa(1,isgn,pnu,paa,pim1,pim2,pipl,jaa)
385 ELSEIF(jak.EQ.6)
THEN
386 CALL dwlukk(1,isgn,pkk,pnu)
389 ELSEIF(jak.EQ.7)
THEN
390 CALL dwluks(1,isgn,pnu,pks,pkk,ppi,jkst)
395 CALL dwlnew(1,isgn,pnu,pwb,pnpi,jak)
403 SUBROUTINE dekay2(IMOD,HH,ISGN)
404 c *******************************
405 c this routine simulates tau- decay
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
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
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)
430 IF(jak2.EQ.-1)
RETURN
435 IF(jak2.EQ.0) CALL jaker(jak)
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)
451 CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
456 ELSEIF(imd.EQ.1)
THEN
457 c =====================
460 nevdec(jak)=nevdec(jak)+1
465 CALL dwluel(2,isgn,pnu,pwb,pmu,pnm)
466 CALL dwrph(ktom,phot)
470 ELSEIF(jak.EQ.2)
THEN
471 CALL dwlumu(2,isgn,pnu,pwb,pmu,pnm)
472 CALL dwrph(ktom,phot)
476 ELSEIF(jak.EQ.3)
THEN
477 CALL dwlupi(2,isgn,ppi,pnu)
481 ELSEIF(jak.EQ.4)
THEN
482 CALL dwluro(2,isgn,pnu,prho,pic,piz)
486 ELSEIF(jak.EQ.5)
THEN
487 CALL dwluaa(2,isgn,pnu,paa,pim1,pim2,pipl,jaa)
490 ELSEIF(jak.EQ.6)
THEN
491 CALL dwlukk(2,isgn,pkk,pnu)
494 ELSEIF(jak.EQ.7)
THEN
495 CALL dwluks(2,isgn,pnu,pks,pkk,ppi,jkst)
500 CALL dwlnew(2,isgn,pnu,pwb,pnpi,jak)
508 SUBROUTINE dexay(KTO,POL)
509 c ----------------------------------------------------------------------
510 c this
'DEXAY' is a routine which generates decay of the single
511 c polarized tau, pol is a polarization vector(not a polarimeter
512 c vector as in dekay) of the tau and it is an input
PARAMETER.
513 c kto=0 initialisation(obligatory)
514 c kto=1 denotes tau+ and kto=2 tau-
515 c dexay(1,pol) and dexay(2,pol) are called internally by mc generator.
516 c decay products are transformed readily
517 c to cms and writen in the lund record in /lujets/
518 c kto=100, print final report(
OPTIONAL).
521 c ----------------------------------------------------------------------
522 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
523 REAL*4 gampmc ,gamper
524 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
526 COMMON /taupos/ np1,np2
527 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
529 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
531 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
534 CHARACTER names(nmode)*31
535 COMMON / inout / inut,iout
537 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4)
546 c initialisation or reinitialisation
547 c first or second tau positions in hepevt as in koralb/z
551 WRITE(iout, 7001) jak1,jak2
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)
569 ELSEIF(kto.EQ.1)
THEN
570 c =====================
571 c decay of tau+ in the tau rest frame
574 IF(iwarm.EQ.0) goto 902
576 cam CALL dexay1(pol,isgn)
577 CALL dexay1(kto,jak1,jakp,pol,isgn)
578 ELSEIF(kto.EQ.2)
THEN
579 c =================================
580 c decay of tau- in the tau rest frame
583 IF(iwarm.EQ.0) goto 902
584 isgn=-idff/iabs(idff)
585 cam CALL dexay2(pol,isgn)
586 CALL dexay1(kto,jak2,jakm,pol,isgn)
587 ELSEIF(kto.EQ.100)
THEN
588 c =======================
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)
601 $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
608 7001
FORMAT(///1x,15(5h*****)
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*,
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*,
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*,
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*
648 chbu
format 7010 had more than 19 continuation lines
650 7010
FORMAT(///1x,15(5h*****)
652 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
653 $ /,
' *', 25x,
'***********DECEMBER 1993***************',9x,1h*,
655 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
656 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
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*
670 $
' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
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*)
680 $ ,i10,2f12.7,a31 ,8x,1h*)
682 $ ,20x,
'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
683 $ /,
' *',20x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
685 902
WRITE(iout, 9020)
686 9020
FORMAT(
' ----- DEXAY: LACK OF INITIALISATION')
688 910
WRITE(iout, 9100)
689 9100
FORMAT(
' ----- DEXAY: WRONG VALUE OF KTO ')
692 SUBROUTINE dexay1(KTO,JAKIN,JAK,POL,ISGN)
693 c ---------------------------------------------------------------------
694 c this routine simulates tau+- decay
697 c ---------------------------------------------------------------------
698 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
699 REAL*4 gampmc ,gamper
700 COMMON / inout / inut,iout
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)
711 IF(jakin.EQ.-1)
RETURN
718 IF(jak.EQ.0) CALL jaker(jak)
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)
745 CALL dexnew(0, isgn,polar,pnu,pwb,pnpi,jnpi)
746 CALL dwlnew(kto,isgn,pnu,pwb,pnpi,jak)
748 nevdec(jak)=nevdec(jak)+1
750 SUBROUTINE dexel(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
751 c ----------------------------------------------------------------------
752 c this simulates tau decay in tau rest frame
753 c into electron and two neutrinos
755 c called by : dexay,dexay1
756 c ----------------------------------------------------------------------
757 REAL pol(4),hv(4),pwb(4),pnu(4),q1(4),q2(4),ph(4),rn(1)
761 c ===================
763 CALL dadmel( -1,isgn,hv,pnu,pwb,q1,q2,ph)
764 cc CALL hbook1(813,
'WEIGHT DISTRIBUTION DEXEL $',100,0,2)
766 ELSEIF(mode.EQ. 0)
THEN
767 c =======================
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.
772 cc CALL hfill(813,wt)
774 IF(rn(1).GT.wt) goto 300
776 ELSEIF(mode.EQ. 1)
THEN
777 c =======================
778 CALL dadmel( 1,isgn,hv,pnu,pwb,q1,q2,ph)
784 9020
FORMAT(
' ----- DEXEL: LACK OF INITIALISATION')
787 SUBROUTINE dexmu(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
788 c ----------------------------------------------------------------------
789 c this simulates tau decay in its rest frame
790 c into muon and two neutrinos
791 c output four momenta: pnu tauneutrino,
795 c ----------------------------------------------------------------------
796 COMMON / inout / inut,iout
797 REAL pol(4),hv(4),pwb(4),pnu(4),q1(4),q2(4),ph(4),rn(1)
801 c ===================
803 CALL dadmmu( -1,isgn,hv,pnu,pwb,q1,q2,ph)
804 cc CALL hbook1(814,
'WEIGHT DISTRIBUTION DEXMU $',100,0,2)
806 ELSEIF(mode.EQ. 0)
THEN
807 c =======================
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.
812 cc CALL hfill(814,wt)
814 IF(rn(1).GT.wt) goto 300
816 ELSEIF(mode.EQ. 1)
THEN
817 c =======================
818 CALL dadmmu( 1,isgn,hv,pnu,pwb,q1,q2,ph)
823 902
WRITE(iout, 9020)
824 9020
FORMAT(
' ----- DEXMU: LACK OF INITIALISATION')
827 SUBROUTINE dadmel(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
828 c ----------------------------------------------------------------------
830 c called by : dexel,(dekay,dekay1)
831 c ----------------------------------------------------------------------
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
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
847 COMMON / inout / inut,iout
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)
856 DATA pi /3.141592653589793238462643/
860 c ===================
869 CALL dphsel(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
870 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
872 cc CALL hbook1(803,
'WEIGHT DISTRIBUTION DADMEL $',100,0,2)
874 ELSEIF(mode.EQ. 0)
THEN
875 c =======================
877 IF(iwarm.EQ.0) goto 902
879 CALL dphsel(wt,hv,pnu,pwb,q1,q2,phx)
880 cc CALL hfill(803,wt/wtmax)
885 IF(wt.GT.wtmax) nevovr=nevovr+1
886 IF(rn*wtmax.GT.wt) goto 300
887 c rotations to basic tau rest frame
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)
906 44 hhv(i)=-isgn*hv(i)
909 ELSEIF(mode.EQ. 1)
THEN
910 c =======================
911 IF(nevraw.EQ.0)
RETURN
912 pargam=swt/float(nevraw+1)
914 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
916 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
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*
935 902
WRITE(iout, 9020)
936 9020
FORMAT(
' ----- DADMEL: LACK OF INITIALISATION')
939 SUBROUTINE dadmmu(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
940 c ----------------------------------------------------------------------
941 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
942 * ,ampiz,ampi,amro,gamro,ama1,gama1
943 * ,amk,amkz,amkst,gamkst
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
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)
958 DATA pi /3.141592653589793238462643/
962 c ===================
971 CALL dphsmu(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
972 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
974 cc CALL hbook1(802,
'WEIGHT DISTRIBUTION DADMMU $',100,0,2)
976 ELSEIF(mode.EQ. 0)
THEN
977 c =======================
979 IF(iwarm.EQ.0) goto 902
981 CALL dphsmu(wt,hv,pnu,pwb,q1,q2,phx)
982 cc CALL hfill(802,wt/wtmax)
987 IF(wt.GT.wtmax) nevovr=nevovr+1
988 IF(rn*wtmax.GT.wt) goto 300
989 c rotations to basic tau rest frame
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)
1006 44 hhv(i)=-isgn*hv(i)
1009 ELSEIF(mode.EQ. 1)
THEN
1010 c =======================
1011 IF(nevraw.EQ.0)
RETURN
1012 pargam=swt/float(nevraw+1)
1014 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1016 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1020 cam nevdec(2)=nevacc
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')
1039 SUBROUTINE dphsel(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
1040 c xnx,xna was flipped in parameters of dphsel and dphsmu
1041 c *********************************************************************
1042 c * electron decay mode *
1043 c *********************************************************************
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)
1049 CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
1060 SUBROUTINE dphsmu(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
1061 c xnx,xna was flipped in parameters of dphsel and dphsmu
1062 c *********************************************************************
1063 c * muon decay mode *
1064 c *********************************************************************
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)
1070 CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
1081 SUBROUTINE drcmu(DGAMT,HV,PH,PAA,XA,QP,XN,IELMU)
1082 IMPLICIT REAL*8 (a-h,o-z)
1083 c ----------------------------------------------------------------------
1084 * it simulates e,mu channels of tau decay in its rest frame with
1085 * qed order alpha corrections
1086 c ----------------------------------------------------------------------
1087 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1088 * ,ampiz,ampi,amro,gamro,ama1,gama1
1089 * ,amk,amkz,amkst,gamkst
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
1097 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1098 REAL*4 gampmc ,gamper
1100 COMMON / inout / inut,iout
1101 COMMON / taurad / xk0dec,itdkrc
1103 REAL*8 hv(4),pt(4),ph(4),paa(4),xa(4),qp(4),xn(4)
1107 DATA pi /3.141592653589793238462643d0/
1109 c ajwmod to satisfy compiler, comment out this unused function.
1111 xlam(x,y,z)=sqrt((x-y-z)**2-4.0*y*z)
1113 c amro, gamro is only a
PARAMETER for geting hight efficiency
1115 c three body phase space normalised as in bjorken-drell
1116 c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
1117 phspac=1./2**17/pi**8
1127 IF (ielmu.EQ.1)
THEN
1134 IF ( itdkrc.EQ.0) prhard=0d0
1136 IF(prsoft.LT.0.1)
THEN
1137 print *,
'ERROR IN DRCMU; PRSOFT=',prsoft
1142 ihard=(rr5.GT.prsoft)
1144 c tau decay to
'TAU+photon'
1146 ams1=(amu+amnuta)**2
1149 xl1=log(xk1/2/xk0dec)
1154 phspac=phspac*ams2*xl1*xk
1155 phspac=phspac/prhard
1158 phspac=phspac*2**6*pi**3
1159 phspac=phspac/prsoft
1161 c mass of neutrina system
1168 am2sq=ams1+ rr2*(ams2-ams1)
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)
1180 * nu light in nunu rest frame
1184 * tau-prim rest frame, define qp(muon
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
1194 qp(4)=1.d0/(2*am3)*(am3**2-am2**2+amu**2)
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)
1204 eps=4*(amu/amtax)**2
1205 xl1=log((2+eps)/eps)
1207 eta =exp(xl1*rr3+xl0)
1210 phspac=phspac*xl1/2*eta
1212 CALL rotpox(thet,phi,xn)
1213 CALL rotpox(thet,phi,xa)
1214 CALL rotpox(thet,phi,qp)
1215 CALL rotpox(thet,phi,pr)
1217 * now to the tau rest frame, define tau-prim and gamma momenta
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)
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)
1238 thet =acos(-1.+2*rr3)
1240 CALL rotpox(thet,phi,xn)
1241 CALL rotpox(thet,phi,xa)
1242 CALL rotpox(thet,phi,qp)
1243 CALL rotpox(thet,phi,pr)
1245 * now to the tau rest frame, define tau-prim and gamma momenta
1257 c 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
1261 SUBROUTINE dampry(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1262 IMPLICIT REAL*8 (a-h,o-z)
1263 c ----------------------------------------------------------------------
1264 c it calculates matrix element for the
1265 c tau --> mu(e) nu nubar decay mode
1266 c including complete order alpha qed corrections.
1267 c ----------------------------------------------------------------------
1268 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1269 * ,ampiz,ampi,amro,gamro,ama1,gama1
1270 * ,amk,amkz,amkst,gamkst
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)
1279 IF(xk(4).LT.0.1d0*ak0)
THEN
1280 amplit=thb(itdkrc,qp,xn,xa,ak0,hv)
1282 amplit=sqm2(itdkrc,qp,xn,xa,xk,ak0,hv)
1286 FUNCTION sqm2(ITDKRC,QP,XN,XA,XK,AK0,HV)
1288 c **********************************************************************
1289 c
REAL photon matrix element squared *
1291 c hv- polarimetric four-vector of tau *
1292 c qp,xn,xa,xk - 4-momenta of electron(muon), nu, nubar and photon *
1293 c all four-vectors in tau rest frame(in gev) *
1294 c ak0 - infrared cutoff, minimal energy of hard photons(gev) *
1295 c sqm2 - value for s=0 *
1296 c see eqs. (2.9)-(2.10) from cjk( nucl.phys.b(1991) ) *
1297 c **********************************************************************
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
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)
1314 REAL*8 s0(3),rxa(3),rxk(3),rqp(3)
1315 DATA pi /3.141592653589793238462643d0/
1321 emass2=qp(4)**2-qp(1)**2-qp(2)**2-qp(3)**2
1323 c scalar products of four-momenta
1329 rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1330 c 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)
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)
1337 c 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)
1349 s1= qpxn*txa*( -emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
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
1357 s0(i) = qpxn*rxa(i)*(-emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
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
1364 FUNCTION thb(ITDKRC,QP,XN,XA,AK0,HV)
1366 c **********************************************************************
1367 c born +virtual+soft photon matrix element**2 o(alpha) *
1369 c hv- polarimetric four-vector of tau *
1370 c qp,xn,xa - four-momenta of electron(muon), nu and nubar in gev *
1371 c all four-vectors in tau rest frame *
1372 c ak0 - infrared cutoff, minimal energy of hard photons *
1373 c thb - value for s=0 *
1374 c see eqs. (2.2),(2.4)-(2.5) from cjk(nucl.phys.b351(1991)70 *
1375 c and(c.2) from jk(nucl.phys.b320(1991)20 ) *
1376 c **********************************************************************
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
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)
1393 REAL*8 rxa(3),rxn(3),rqp(3)
1394 REAL*8 bornpl(3),am3pol(3),xm3pol(3)
1395 DATA pi /3.141592653589793238462643d0/
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)
1410 c 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)
1413 c quasi two-body variables
1415 u3=sqrt(qp(1)**2+qp(2)**2+qp(3)**2)/tmass
1417 w0=(xn(4)+xa(4))/tmass
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
1436 c 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)
1443 c 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 )
1449 c 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
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
1458 c v-a and v+a couplings, but in the born part only
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) )*
1465 5 hv(i)=(bornpl(i)+am3pol(i))/(born+am3)-1.d0
1467 IF (thb/born.LT.0.1d0)
THEN
1468 print *,
'ERROR IN THB, THB/BORN=',thb/born
1477 SUBROUTINE dexpi(MODE,ISGN,POL,PPI,PNU)
1478 c ----------------------------------------------------------------------
1479 c tau decay into pion and tau-neutrino
1481 c output four momenta: pnu tauneutrino,
1483 c ----------------------------------------------------------------------
1484 REAL pol(4),hv(4),pnu(4),ppi(4),rn(1)
1487 c ===================
1488 CALL dadmpi(-1,isgn,hv,ppi,pnu)
1489 cc CALL hbook1(815,
'WEIGHT DISTRIBUTION DEXPI $',100,0,2)
1491 ELSEIF(mode.EQ. 0)
THEN
1492 c =======================
1494 CALL dadmpi( 0,isgn,hv,ppi,pnu)
1495 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1496 cc CALL hfill(815,wt)
1498 IF(rn(1).GT.wt) goto 300
1500 ELSEIF(mode.EQ. 1)
THEN
1501 c =======================
1502 CALL dadmpi( 1,isgn,hv,ppi,pnu)
1508 SUBROUTINE dadmpi(MODE,ISGN,HV,PPI,PNU)
1509 c ----------------------------------------------------------------------
1510 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1511 * ,ampiz,ampi,amro,gamro,ama1,gama1
1512 * ,amk,amkz,amkst,gamkst
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/
1526 c ===================
1528 ELSEIF(mode.EQ. 0)
THEN
1529 c =======================
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)
1535 CALL sphera(xpi,ppi)
1537 c tau-neutrino momentum
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
1547 40 hv(i)=-isgn*2*ga*gv*amtau*(2*ppi(i)*qxn-pnu(i)*ampi**2)/brak
1550 ELSEIF(mode.EQ. 1)
THEN
1551 c =======================
1552 IF(nevtot.EQ.0)
RETURN
1554 c gamm=(gfermi*fpi)**2/(16.*pi)*amtau**3*
1555 c * (brak/amtau**4)**2
1556 czw 7.02.93 here was an error affecting non standard model
1557 c configurations only
1558 gamm=(gfermi*fpi)**2/(16.*pi)*amtau**3*
1560 $ sqrt((amtau**2-ampi**2-amnuta**2)**2
1561 $ -4*ampi**2*amnuta**2 )/amtau**2
1564 WRITE(iout, 7010) nevtot,gamm,rat,error
1567 cam nevdec(3)=nevtot
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*****)/)
1579 SUBROUTINE dexro(MODE,ISGN,POL,PNU,PRO,PIC,PIZ)
1580 c ----------------------------------------------------------------------
1581 c this simulates tau decay in tau rest frame
1582 c into nu rho, next rho decays into pion pair.
1583 c output four momenta: pnu tauneutrino,
1587 c ----------------------------------------------------------------------
1588 COMMON / inout / inut,iout
1589 REAL pol(4),hv(4),pro(4),pnu(4),pic(4),piz(4),rn(1)
1593 c ===================
1595 CALL dadmro( -1,isgn,hv,pnu,pro,pic,piz)
1596 cc CALL hbook1(816,
'WEIGHT DISTRIBUTION DEXRO $',100,0,2)
1597 cc CALL hbook1(916,
'ABS2 OF HV IN ROUTINE DEXRO $',100,0,2)
1599 ELSEIF(mode.EQ. 0)
THEN
1600 c =======================
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.
1605 cc CALL hfill(816,wt)
1606 cc xhelp=hv(1)**2+hv(2)**2+hv(3)**2
1607 cc CALL hfill(916,xhelp)
1609 IF(rn(1).GT.wt) goto 300
1611 ELSEIF(mode.EQ. 1)
THEN
1612 c =======================
1613 CALL dadmro( 1,isgn,hv,pnu,pro,pic,piz)
1619 902
WRITE(iout, 9020)
1620 9020
FORMAT(
' ----- DEXRO: LACK OF INITIALISATION')
1623 SUBROUTINE dadmro(MODE,ISGN,HHV,PNU,PRO,PIC,PIZ)
1624 c ----------------------------------------------------------------------
1625 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1626 * ,ampiz,ampi,amro,gamro,ama1,gama1
1627 * ,amk,amkz,amkst,gamkst
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
1638 REAL hv(4),pro(4),pnu(4),pic(4),piz(4)
1639 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4)
1642 DATA pi /3.141592653589793238462643/
1646 c ===================
1655 CALL dphsro(wt,hv,pdum1,pdum2,pdum3,pdum4)
1656 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1658 cc CALL hbook1(801,
'WEIGHT DISTRIBUTION DADMRO $',100,0,2)
1661 ELSEIF(mode.EQ. 0)
THEN
1662 c =======================
1664 IF(iwarm.EQ.0) goto 902
1665 CALL dphsro(wt,hv,pnu,pro,pic,piz)
1666 cc CALL hfill(801,wt/wtmax)
1672 IF(wt.GT.wtmax) nevovr=nevovr+1
1673 IF(rn*wtmax.GT.wt) goto 300
1674 c rotations to basic tau rest frame
1675 costhe=-1.+2.*rrr(2)
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)
1689 44 hhv(i)=-isgn*hv(i)
1692 ELSEIF(mode.EQ. 1)
THEN
1693 c =======================
1694 IF(nevraw.EQ.0)
RETURN
1695 pargam=swt/float(nevraw+1)
1697 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1699 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1703 cam nevdec(4)=nevacc
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')
1724 SUBROUTINE dphsro(DGAMT,HV,PN,PR,PIC,PIZ)
1725 c ----------------------------------------------------------------------
1726 c it simulates rho decay in tau rest frame with
1727 c z-axis along rho momentum
1728 c ----------------------------------------------------------------------
1729 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1730 * ,ampiz,ampi,amro,gamro,ama1,gama1
1731 * ,amk,amkz,amkst,gamkst
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/
1742 c three body phase space normalised as in bjorken-drell
1743 phspac=1./2**11/pi**5
1749 c mass of(real/virtual) rho
1750 ams1=(ampi+ampiz)**2
1751 ams2=(amtau-amnuta)**2
1754 c amx2=ams1+ rr1(1)*(ams2-ams1)
1756 c amx2=ams1+ rr1*(ams2-ams1)
1759 c phspac=phspac*(ams2-ams1)
1760 c phase space with sampling for rho resonance
1761 alp1=atan((ams1-amro**2)/amro/gamro)
1762 alp2=atan((ams2-amro**2)/amro/gamro)
1766 alp=alp1+rr1(1)*(alp2-alp1)
1767 amx2=amro**2+amro*gamro*tan(alp)
1769 IF(amx.LT.2.*ampi) go to 100
1771 phspac=phspac*((amx2-amro**2)**2+(amro*gamro)**2)/(amro*gamro)
1772 phspac=phspac*(alp2-alp1)
1774 c tau-neutrino momentum
1777 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
1778 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
1782 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
1784 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
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)
1791 c charged pi momentum in rho rest frame
1792 CALL sphera(pppi,pic)
1794 c neutral pi momentum in rho rest frame
1798 exe=(pr(4)+pr(3))/amx
1799 c pions boosted from rho rest frame to tau rest frame
1800 CALL bostr3(exe,pic,pic)
1801 CALL bostr3(exe,piz,piz)
1803 30 qq(i)=pic(i)-piz(i)
1806 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
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
1814 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
1817 SUBROUTINE dexaa(MODE,ISGN,POL,PNU,PAA,PIM1,PIM2,PIPL,JAA)
1818 c ----------------------------------------------------------------------
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,
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
1827 c ----------------------------------------------------------------------
1828 COMMON / inout / inut,iout
1829 REAL pol(4),hv(4),paa(4),pnu(4),pim1(4),pim2(4),pipl(4),rn(1)
1833 c ===================
1835 CALL dadmaa( -1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1836 cc CALL hbook1(816,
'WEIGHT DISTRIBUTION DEXAA $',100,-2.,2.)
1838 ELSEIF(mode.EQ. 0)
THEN
1839 * =======================
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.
1844 cc CALL hfill(816,wt)
1846 IF(rn(1).GT.wt) goto 300
1848 ELSEIF(mode.EQ. 1)
THEN
1849 * =======================
1850 CALL dadmaa( 1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1855 902
WRITE(iout, 9020)
1856 9020
FORMAT(
' ----- DEXAA: LACK OF INITIALISATION')
1859 SUBROUTINE dadmaa(MODE,ISGN,HHV,PNU,PAA,PIM1,PIM2,PIPL,JAA)
1860 c ----------------------------------------------------------------------
1861 * a1 decay unweighted events
1862 c ----------------------------------------------------------------------
1863 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1864 * ,ampiz,ampi,amro,gamro,ama1,gama1
1865 * ,amk,amkz,amkst,gamkst
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
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)
1880 DATA pi /3.141592653589793238462643/
1884 c ===================
1893 CALL dphsaa(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5,jaa)
1894 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1896 cc CALL hbook1(801,
'WEIGHT DISTRIBUTION DADMAA $',100,0,2)
1898 ELSEIF(mode.EQ. 0)
THEN
1899 c =======================
1901 IF(iwarm.EQ.0) goto 902
1902 CALL dphsaa(wt,hv,pnu,paa,pim1,pim2,pipl,jaa)
1903 cc CALL hfill(801,wt/wtmax)
1911 sswt=sswt+dble(wt)**2
1916 IF(wt.GT.wtmax) nevovr=nevovr+1
1917 IF(rn*wtmax.GT.wt) goto 300
1918 c rotations to basic tau rest frame
1919 costhe=-1.+2.*rrr(2)
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)
1929 44 hhv(i)=-isgn*hv(i)
1932 ELSEIF(mode.EQ. 1)
THEN
1933 c =======================
1934 IF(nevraw.EQ.0)
RETURN
1935 pargam=swt/float(nevraw+1)
1937 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1939 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1943 cam nevdec(5)=nevacc
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')
1964 SUBROUTINE dphsaa(DGAMT,HV,PN,PAA,PIM1,PIM2,PIPL,JAA)
1965 c ----------------------------------------------------------------------
1966 * it simulates a1 decay in tau rest frame with
1967 * z-axis along a1 momentum
1968 c ----------------------------------------------------------------------
1969 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1970 * ,ampiz,ampi,amro,gamro,ama1,gama1
1971 * ,amk,amkz,amkst,gamkst
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)
1982 c matrix element number:
1984 c
TYPE of the generation:
1988 IF (rmod.LT.bra1)
THEN
2000 $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
2002 SUBROUTINE dexkk(MODE,ISGN,POL,PKK,PNU)
2003 c ----------------------------------------------------------------------
2004 c tau decay into kaon and tau-neutrino
2006 c output four momenta: pnu tauneutrino,
2008 c ----------------------------------------------------------------------
2009 REAL pol(4),hv(4),pnu(4),pkk(4),rn(1)
2012 c ===================
2013 CALL dadmkk(-1,isgn,hv,pkk,pnu)
2014 cc CALL hbook1(815,
'WEIGHT DISTRIBUTION DEXPI $',100,0,2)
2016 ELSEIF(mode.EQ. 0)
THEN
2017 c =======================
2019 CALL dadmkk( 0,isgn,hv,pkk,pnu)
2020 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
2021 cc CALL hfill(815,wt)
2023 IF(rn(1).GT.wt) goto 300
2025 ELSEIF(mode.EQ. 1)
THEN
2026 c =======================
2027 CALL dadmkk( 1,isgn,hv,pkk,pnu)
2033 SUBROUTINE dadmkk(MODE,ISGN,HV,PKK,PNU)
2034 c ----------------------------------------------------------------------
2038 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2039 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2041 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2042 * ,ampiz,ampi,amro,gamro,ama1,gama1
2043 * ,amk,amkz,amkst,gamkst
2045 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2046 * ,ampiz,ampi,amro,gamro,ama1,gama1
2047 * ,amk,amkz,amkst,gamkst
2049 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2050 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
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/
2060 c ===================
2062 ELSEIF(mode.EQ. 0)
THEN
2063 c =======================
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)
2069 CALL sphera(xkk,pkk)
2071 c tau-neutrino momentum
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
2081 40 hv(i)=-isgn*2*ga*gv*amtau*(2*pkk(i)*qxn-pnu(i)*amk**2)/brak
2084 ELSEIF(mode.EQ. 1)
THEN
2085 c =======================
2086 IF(nevtot.EQ.0)
RETURN
2088 cfz there was brak/amtau**4 before
2089 c gamm=(gfermi*fkk)**2/(16.*pi)*amtau**3*
2090 c * (brak/amtau**4)**2
2091 czw 7.02.93 here was an error affecting non standard model
2092 c configurations only
2093 gamm=(gfermi*fkk)**2/(16.*pi)*amtau**3*
2095 $ sqrt((amtau**2-amk**2-amnuta**2)**2
2096 $ -4*amk**2*amnuta**2 )/amtau**2
2101 WRITE(iout, 7010) nevtot,gamm,rat,error
2104 cam nevdec(6)=nevtot
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*****)/)
2116 SUBROUTINE dexks(MODE,ISGN,POL,PNU,PKS,PKK,PPI,JKST)
2117 c ----------------------------------------------------------------------
2118 c this simulates tau decay in tau rest frame
2119 c into nu k*,
THEN k* decays into pi0,k+-(jkst=20)
2120 c or pi+-,k0(jkst=10).
2121 c output four momenta: pnu tauneutrino,
2127 c ----------------------------------------------------------------------
2128 COMMON / inout / inut,iout
2129 REAL pol(4),hv(4),pks(4),pnu(4),pkk(4),ppi(4),rn(1)
2133 c ===================
2135 cfz initialisation done with the gharged pion neutral kaon mode(jkst=10
2136 CALL dadmks( -1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2137 cc CALL hbook1(816,
'WEIGHT DISTRIBUTION DEXKS $',100,0,2)
2138 cc CALL hbook1(916,
'ABS2 OF HV IN ROUTINE DEXKS $',100,0,2)
2140 ELSEIF(mode.EQ. 0)
THEN
2141 c =======================
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.
2146 cc CALL hfill(816,wt)
2147 cc xhelp=hv(1)**2+hv(2)**2+hv(3)**2
2148 cc CALL hfill(916,xhelp)
2150 IF(rn(1).GT.wt) goto 300
2152 ELSEIF(mode.EQ. 1)
THEN
2153 c ======================================
2154 CALL dadmks( 1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2160 902
WRITE(iout, 9020)
2161 9020
FORMAT(
' ----- DEXKS: LACK OF INITIALISATION')
2164 SUBROUTINE dadmks(MODE,ISGN,HHV,PNU,PKS,PKK,PPI,JKST)
2165 c ----------------------------------------------------------------------
2166 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2167 * ,ampiz,ampi,amro,gamro,ama1,gama1
2168 * ,amk,amkz,amkst,gamkst
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
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)
2185 DATA pi /3.141592653589793238462643/
2189 c ===================
2198 c the initialisation is done with the 66.7% MODE
2200 CALL dphsks(wt,hv,pdum1,pdum2,pdum3,pdum4,jkst)
2201 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
2203 cc CALL hbook1(801,
'WEIGHT DISTRIBUTION DADMKS $',100,0,2)
2205 cc CALL hbook1(112,
'-------- K* MASS -------- $',100,0.,2.)
2206 ELSEIF(mode.EQ. 0)
THEN
2207 c =====================================
2208 IF(iwarm.EQ.0) goto 902
2209 c here we choose randomly between k0 pi+_(66.7%)
2210 c and k+_ pi0(33.3%)
2214 IF(rmod(1).LT.dec1)
THEN
2219 CALL dphsks(wt,hv,pnu,pks,pkk,ppi,jkst)
2222 IF(wt.GT.wtmax) nevovr=nevovr+1
2226 IF(rn*wtmax.GT.wt) goto 400
2227 c rotations to basic tau rest frame
2228 costhe=-1.+2.*rrr(2)
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)
2242 44 hhv(i)=-isgn*hv(i)
2245 ELSEIF(mode.EQ. 1)
THEN
2246 c =======================
2247 IF(nevraw.EQ.0)
RETURN
2248 pargam=swt/float(nevraw+1)
2250 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
2252 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
2256 cam nevdec(7)=nevacc
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')
2277 SUBROUTINE dphsks(DGAMT,HV,PN,PKS,PKK,PPI,JKST)
2278 c ----------------------------------------------------------------------
2279 c it simulates kaon* decay in tau rest frame with
2280 c z-axis along kaon* momentum
2281 c jkst=10 for k* --->k0 + pi+-
2282 c jkst=20 for k* --->k+- + pi0
2283 c ----------------------------------------------------------------------
2286 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2287 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2289 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2290 * ,ampiz,ampi,amro,gamro,ama1,gama1
2291 * ,amk,amkz,amkst,gamkst
2293 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2294 * ,ampiz,ampi,amro,gamro,ama1,gama1
2295 * ,amk,amkz,amkst,gamkst
2297 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2298 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2301 REAL hv(4),pt(4),pn(4),pks(4),pkk(4),ppi(4),qq(4),rr1(1)
2308 DATA pi /3.141592653589793238462643/
2311 c three body phase space normalised as in bjorken-drell
2312 phspac=1./2**11/pi**5
2319 c here begin the k0,pi+_ decay
2321 c ==================
2322 c mass of(real/virtual) k*
2324 ams2=(amtau-amnuta)**2
2326 c amx2=ams1+ rr1(1)*(ams2-ams1)
2328 c phspac=phspac*(ams2-ams1)
2329 c 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)
2335 phspac=phspac*((amx2-amkst**2)**2+(amkst*gamkst)**2)
2337 phspac=phspac*(alp2-alp1)
2339 c tau-neutrino momentum
2342 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2343 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2348 pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2350 phspac=phspac*(4*pi)*(2*pks(3)/amtau)
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)
2356 c charged pi momentum in kaon* rest frame
2357 CALL sphera(pppi,ppi)
2359 c neutral kaon momentum in k* rest frame
2362 pkk(4)=( amx**2+amkz**2-ampi**2 ) / ( 2*amx )
2363 exe=(pks(4)+pks(3))/amx
2364 c pion and k boosted from k* rest frame to tau rest frame
2365 CALL bostr3(exe,ppi,ppi)
2366 CALL bostr3(exe,pkk,pkk)
2368 30 qq(i)=ppi(i)-pkk(i)
2369 c 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)
2373 31 qq(i)=qq(i)-pks(i)*qqpks/pksd
2376 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
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
2381 c a simple breit-wigner is chosen for k* resonance
2383 cam fks=cabs(bwigs(amx2,amkst,gamkst))**2
2384 fks=cabs(bwigm(amx2,amkst,gamkst,ampi,amkz))**2
2386 fks=cabs(bwigs(amx2,amkst,gamkst))**2
2388 amplit=(gfermi*scabib)**2*brak*2*fks
2389 dgamt=1/(2.*amtau)*amplit*phspac
2391 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
2393 c here begin the k+-,pi0 decay
2394 ELSEIF(jkst.EQ.20)
THEN
2395 c ======================
2396 c mass of(real/virtual) k*
2398 ams2=(amtau-amnuta)**2
2401 c amx2=ams1+ rr1(1)*(ams2-ams1)
2403 c amx2=ams1+ rr1*(ams2-ams1)
2406 c phspac=phspac*(ams2-ams1)
2407 c 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)
2413 phspac=phspac*((amx2-amkst**2)**2+(amkst*gamkst)**2)
2415 phspac=phspac*(alp2-alp1)
2417 c tau-neutrino momentum
2420 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2421 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2425 pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2427 phspac=phspac*(4*pi)*(2*pks(3)/amtau)
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)
2433 c neutral pi momentum in k* rest frame
2434 CALL sphera(pppi,ppi)
2436 c charged kaon momentum in k* rest frame
2439 pkk(4)=( amx**2+amk**2-ampiz**2 ) / ( 2*amx )
2440 exe=(pks(4)+pks(3))/amx
2441 c pion and k boosted from k* rest frame to tau rest frame
2442 CALL bostr3(exe,ppi,ppi)
2443 CALL bostr3(exe,pkk,pkk)
2445 60 qq(i)=pkk(i)-ppi(i)
2446 c 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)
2450 61 qq(i)=qq(i)-pks(i)*qqpks/pksd
2453 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
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
2458 c a simple breit-wigner is chosen for the k* resonance
2460 cam fks=cabs(bwigs(amx2,amkst,gamkst))**2
2461 fks=cabs(bwigm(amx2,amkst,gamkst,amk,ampiz))**2
2463 fks=cabs(bwigs(amx2,amkst,gamkst))**2
2465 amplit=(gfermi*scabib)**2*brak*2*fks
2466 dgamt=1/(2.*amtau)*amplit*phspac
2468 70 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
2476 SUBROUTINE dphnpi(DGAMT,HV,PN,PR,PPI,JNPI)
2478 SUBROUTINE dphnpi(DGAMT,HVX,PNX,PRX,PPIX,JNPI)
2480 c ----------------------------------------------------------------------
2481 c it simulates multipi decay in tau rest frame with
2482 c z-axis opposite to neutrino momentum
2483 c ----------------------------------------------------------------------
2484 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2485 * ,ampiz,ampi,amro,gamro,ama1,gama1
2486 * ,amk,amkz,amkst,gamkst
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)
2495 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2497 CHARACTER names(nmode)*31
2500 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2502 CHARACTER names(nmode)*31
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)
2512 DATA pi /3.141592653589793238462643/
2513 DATA dpar/2.,5.,15.,60.,250.,1500.,1.2e4,1.2e5/
2515 c 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)
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
2525 REAL*8 gam,bep,phi,a,b,c
2527 REAL*4 rrr(9),rrx(2),rn(1),rr2(1)
2529 DATA pi /3.141592653589793238462643/
2530 DATA wetmax /20*1d-15/
2532 cc-- pawt(a,b,c)=sqrt((a**2-(b+c)**2)*(a**2-(b-c)**2))/(2.*a)
2535 $ sqrt(max(0.d0,(a**2-(b+c)**2)*(a**2-(b-c)**2)))/(2.d0*a)
2538 ampik(i,j)=dcdmas(idffin(i,j))
2543 IF ((jnpi.LE.0).OR.jnpi.GT.20)
THEN
2544 WRITE(6,*)
'JNPI OUTSIDE RANGE DEFINED BY WETMAX; JNPI=',jnpi
2562 phspac = 1./2.**5 /pi**2
2564 4 ps =ps+ampik(i,jnpi)
2571 ams2=(amtau-amnuta)**2
2575 amx2=ams1+ rr1(1)*(ams2-ams1)
2577 amx2=ams1+ rr2(1)*(ams2-ams1)
2581 phspac=phspac * (ams2-ams1)
2583 c tau-neutrino momentum
2586 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx2)
2587 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2591 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx2)
2593 phspac=phspac * (4.*pi) * (2.*pr(3)/amtau)
2595 c amplitude(cf ys.tsai phys.rev.d4,2821(1971)
2596 c or f.gilman sh.rhie phys.rev.d31,1066(1985)
2600 qxn=pr(4)*pn(4)-pr(1)*pn(1)-pr(2)*pn(2)-pr(3)*pn(3)
2603 c here was an error. 20.10.91 (zw)
2604 c brak=2*(gv**2+ga**2)*(2*pxq*pxn+amx2*qxn)
2606 brak=2*(gv**2+ga**2)*(2*pxq*qxn+amx2*pxn)
2607 & -6*(gv**2-ga**2)*amtau*amnuta*amx2
2608 cam assume neutrino mass=0. and sum over final polarisation
2609 c 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
2613 c isotropic w decay in w rest frame
2615 phspac = 1./2.**(6*nd-7) /pi**(3*nd-4)
2616 phsmax = 1./dpar(nd-2)
2623 pv(5,nd)=ampik(nd,jnpi)
2624 c compute max. phase space factor
2625 pmax=amw-ps+ampik(nd,jnpi)
2628 pmax=pmax+ampik(il,jnpi)
2629 pmin=pmin+ampik(il+1,jnpi)
2631 220 phsmax=phsmax*pawt(pmax,pmin,ampik(il,jnpi))
2632 cam generate nd-2 effective masses(cf ludecy)
2633 phspac = 1./2.**(6*nd-7) /pi**(3*nd-4)
2635 CALL ranmar(rrr,nd-1)
2639 IF(rsav.LE.rord(jl)) goto 260
2640 250 rord(jl+1)=rord(jl)
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))
2649 IF(phs.LT.rn*phsmax) goto 240
2651 220 phsmax=phsmax*pawt(pmax,pmin,ampik(il,jnpi))/pmax
2653 c --- 2.02.94 zw 9 lines
2658 223 ams1=ams1+ampik(jl,jnpi)
2660 amx =(amx-ampik(il,jnpi))
2662 phsmax=phsmax * (ams2-ams1)
2667 cam generate nd-2 effective masses
2669 phspac = 1./2.**(6*nd-7) /pi**(3*nd-4)
2671 CALL ranmar(rrr,nd-2)
2675 231 ams1=ams1+ampik(jl,jnpi)
2677 ams2=(amx-ampik(il,jnpi))**2
2679 amx2=ams1+ rr1*(ams2-ams1)
2682 phspac=phspac * (ams2-ams1)
2683 c --- 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)
2688 pa=pawt(pv(5,nd-1),ampik(nd-1,jnpi),ampik(nd,jnpi))
2689 phs =phs *pa/pv(5,nd-1)
2691 wetmax(jnpi)=1.2d0*max(wetmax(jnpi)/1.2d0,phs/phsmax)
2692 IF (ncont.EQ.500 000)
THEN
2695 xnpi=xnpi+ampik(kk,jnpi)
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
2703 IF(rn(1)*phsmax*wetmax(jnpi).GT.phs) go to 100
2705 c...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))
2712 ue(1)=sqrt(1.-ue(3)**2)*cos(phi)
2713 ue(2)=sqrt(1.-ue(3)**2)*sin(phi)
2718 ue(1)=sqrt(1.d0-ue(3)**2)*cos(phi)
2719 ue(2)=sqrt(1.d0-ue(3)**2)*sin(phi)
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))
2728 c...lorentz transform decay products to tau frame
2730 310 ppi(j,nd)=pv(j,nd)
2733 320 be(j)=pv(j,il)/pv(4,il)
2734 gam=pv(4,il)/pv(5,il)
2736 bep=be(1)*ppi(1,i)+be(2)*ppi(2,i)+be(3)*ppi(3,i)
2739 330 ppi(j,i)=ppi(j,i)+gam*(gam*bep/(1.+gam)+ppi(4,i))*be(j)
2741 330 ppi(j,i)=ppi(j,i)+gam*(gam*bep/(1.d0+gam)+ppi(4,i))*be(j)
2743 ppi(4,i)=gam*(ppi(4,i)+bep)
2763 FUNCTION sigee(Q2,JNP)
2764 c ----------------------------------------------------------------------
2765 c e+e- cross section in the(1.gev2,amtau**2) region
2766 c normalised to sig0 = 4/3 pi alfa2
2767 c used in matrix element for multipion tau decays
2768 c cf ys.tsai phys.rev d4 ,2821(1971)
2769 c f.gilman et al phys.rev d17,1846(1978)
2770 c c.kiesling, to be pub. in high energy e+e- physics(1988)
2771 c datsig(*,1) = e+e- -> pi+pi-2pi0
2772 c datsig(*,2) = e+e- -> 2pi+2pi-
2773 c datsig(*,3) = 5-pion contribution(a la tn.pham et al)
2774 c(phys lett 78b,623(1978)
2775 c datsig(*,5) = e+e- -> 6pi
2777 c 4- and 6-pion cross sections from
data
2778 c 5-pion contribution related to 4-pion cross section
2781 c ----------------------------------------------------------------------
2782 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2783 * ,ampiz,ampi,amro,gamro,ama1,gama1
2784 * ,amk,amkz,amkst,gamkst
2786 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2787 * ,ampiz,ampi,amro,gamro,ama1,gama1
2788 * ,amk,amkz,amkst,gamkst
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,
2798 7 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25,
2801 DATA pi /3.141592653589793238462643/
2810 c ajwmod: initialize
if called from outside qq:
2811 IF (ampi.LT.0.139) ampi = 0.1395675
2817 datsig(i,2) = datsig(i,2)/2.
2818 datsig(i,1) = datsig(i,1) + datsig(i,2)
2824 IF(t . gt. s-ampi ) go to 201
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)
2833 c
WRITE(6,1000) datsig
2834 1000
FORMAT(///1x,
' EE SIGMA USED IN MULTIPI DECAYS'/
2840 sigee=datsig(1,jnpi)+
2841 & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
2842 ELSEIF(q.LT.1.8)
THEN
2845 IF(q.LT.qmax) go to 2
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
2854 IF(sigee.LT..0) sigee=0.
2856 sigee = sigee/(6.*pi**2*sig0)
2861 FUNCTION sigold(Q2,JNPI)
2862 c ----------------------------------------------------------------------
2863 c e+e- cross section in the(1.gev2,amtau**2) region
2864 c normalised to sig0 = 4/3 pi alfa2
2865 c used in matrix element for multipion tau decays
2866 c cf ys.tsai phys.rev d4 ,2821(1971)
2867 c f.gilman et al phys.rev d17,1846(1978)
2868 c c.kiesling, to be pub. in high energy e+e- physics(1988)
2869 c datsig(*,1) = e+e- -> pi+pi-2pi0
2870 c datsig(*,2) = e+e- -> 2pi+2pi-
2871 c datsig(*,3) = 5-pion contribution(a la tn.pham et al)
2872 c(phys lett 78b,623(1978)
2873 c datsig(*,4) = e+e- -> 6pi
2875 c 4- and 6-pion cross sections from
data
2876 c 5-pion contribution related to 4-pion cross section
2879 c ----------------------------------------------------------------------
2880 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2881 * ,ampiz,ampi,amro,gamro,ama1,gama1
2882 * ,amk,amkz,amkst,gamkst
2884 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2885 * ,ampiz,ampi,amro,gamro,ama1,gama1
2886 * ,amk,amkz,amkst,gamkst
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,
2895 6 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25/
2897 DATA pi /3.141592653589793238462643/
2905 datsig(i,2) = datsig(i,2)/2.
2906 datsig(i,1) = datsig(i,1) + datsig(i,2)
2912 IF(t . gt. s-ampi ) go to 201
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
2919 c
WRITE(6,1000) datsig
2920 1000
FORMAT(///1x,
' EE SIGMA USED IN MULTIPI DECAYS'/
2926 sigee=datsig(1,jnpi)+
2927 & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
2928 ELSEIF(q.LT.1.8)
THEN
2931 IF(q.LT.qmax) go to 2
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
2940 IF(sigee.LT..0) sigee=0.
2942 sigee = sigee/(6.*pi**2*sig0)
2947 SUBROUTINE dphspk(DGAMT,HV,PN,PAA,PNPI,JAA)
2948 c ----------------------------------------------------------------------
2949 * it simulates three pi(k) decay in the tau rest frame
2950 * z-axis along hadronic system
2951 c ----------------------------------------------------------------------
2952 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
2954 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2956 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2959 CHARACTER names(nmode)*31
2961 REAL hv(4),pn(4),paa(4),pim1(4),pim2(4),pipl(4),pnpi(4,9)
2962 c matrix element number:
2964 c
TYPE of the generation:
2967 c --- 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))
2972 $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
2984 $ dphtre(dgamt,hv,pn,paa,pim1,ampa,pim2,ampb,pipl,amp3,keyt,mnum)
2985 c ----------------------------------------------------------------------
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.
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
3000 c ----------------------------------------------------------------------
3001 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3002 * ,ampiz,ampi,amro,gamro,ama1,gama1
3003 * ,amk,amkz,amkst,gamkst
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)
3013 DATA pi /3.141592653589793238462643/
3015 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
3016 c amro, gamro is only a
PARAMETER for geting hight efficiency
3018 c three body phase space normalised as in bjorken-drell
3019 c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
3020 phspac=1./2**17/pi**8
3030 CALL choice(mnum,rr,ichan,prob1,prob2,prob3,
3031 $ amrx,gamrx,amra,gamra,amrb,gamrb)
3032 IF (ichan.EQ.1)
THEN
3035 ELSEIF (ichan.EQ.2)
THEN
3044 ams1=(amp1+amp2+amp3)**2
3045 ams2=(amtau-amnuta)**2
3047 c phase space with sampling for a1 resonance
3049 * phase space with sampling for a1 resonance
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)
3056 phspac=phspac*((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
3057 phspac=phspac*(alp2-alp1)
3058 c mass of(real/virtual) rho -
3062 IF (ichan.LE.2)
THEN
3064 c phase space with sampling for rho resonance,
3066 * phase space with sampling for rho resonance,
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)
3073 c --- this part of the jacobian will be recovered later ---------------
3074 c phspac=phspac*(alp2-alp1)
3075 c phspac=phspac*((am2sq-amra**2)**2+(amra*gamra)**2)/(amra*gamra)
3076 c----------------------------------------------------------------------
3083 am2sq=ams1+ rr2*(ams2-ams1)
3088 c rho restframe, define pipl and pim1
3090 * rho restframe, define pipl and pim1
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))
3096 c --- this part of jacobian will be recovered later
3097 phf1=(4*pi)*(2*pppi/am2)
3099 c pi minus momentum in rho rest frame
3101 * pi minus momentum in rho rest frame
3103 CALL sphera(pppi,pipl)
3106 c pi0 1 momentum in rho rest frame
3108 * pi0 1 momentum in rho rest frame
3114 c a1 rest frame, define pim2
3116 * a1 rest frame, define pim2
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
3127 pim2(4)=1./(2*am3)*(am3**2-am2**2+amp1**2)
3129 phf2=(4*pi)*(2*pr(3)/am3)
3131 c old pions boosted from rho rest frame to a1 rest frame
3133 * old pions boosted from rho rest frame to a1 rest frame
3135 exe=(pr(4)+pr(3))/am2
3136 CALL bostr3(exe,pipl,pipl)
3137 CALL bostr3(exe,pim1,pim1)
3144 thet =acos(-1.+2*rr3)
3146 CALL rotpol(thet,phi,pipl)
3147 CALL rotpol(thet,phi,pim1)
3148 CALL rotpol(thet,phi,pim2)
3149 CALL rotpol(thet,phi,pr)
3151 * now to the tau rest frame, define a1 and neutrino momenta
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
3162 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am3**2)
3164 c here we correct for the jacobians of the two chains
3165 c ---first channel ------- pim1+pipl
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
3173 c jacobian of speeding
3174 ff1 = ((am2sq-amra**2)**2+(amra*gamra)**2)/(amra*gamra)
3175 ff1 =ff1 *(alp2-alp1)
3176 c lambda of rho decay
3177 gg1 = (4*pi)*(xlam(am2sq,amp2**2,amp3**2)/am2sq)
3178 c lambda of a1 decay
3179 gg1 =gg1 *(4*pi)*sqrt(4*xpro/am3sq)
3180 xjaje=gg1*(ams2-ams1)
3181 c ---second channel ------ pim2+pipl
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)
3200 IF (ichan.EQ.2)
THEN
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
3209 IF (a1+a2+a3.NE.0.0)
THEN
3210 phspac=phspac/(a1+a2+a3)
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)
3227 c partial width consists of phase space and amplitude
3229 CALL dampog(pt,pn,pim1,pim2,pipl,amplit,hv)
3230 c
ELSEIF (mnum.EQ.0)
THEN
3231 c CALL dampaa(pt,pn,pim1,pim2,pipl,amplit,hv)
3233 CALL damppk(mnum,pt,pn,pim1,pim2,pipl,amplit,hv)
3235 IF (keyt.EQ.1.OR.keyt.EQ.2)
THEN
3236 c the statistical factor for identical pi-s is cancelled with
3237 c two, for two modes of a1 decay namelly pi+pi-pi- and pi-pi0pi0
3239 cam phspac=phspac*2.0
3240 cam phspac=phspac/2.
3246 dgamt=1/(2.*amtau)*amplit*phspac
3248 SUBROUTINE dampaa(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3249 c ----------------------------------------------------------------------
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.
3256 c called by : dphsaa
3257 c ----------------------------------------------------------------------
3258 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3259 * ,ampiz,ampi,amro,gamro,ama1,gama1
3260 * ,amk,amkz,amkst,gamkst
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
3274 * f constants for a1, a1-rho-pi, and rho-pi-pi
3277 * this inline funct. calculates the scalar part of the propagator
3278 bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
3280 * four momentum of a1
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))
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
3302 fnorm=fa1/sqrt(2.)*faropi*fro2pi
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))
3309 fnorm=2.0*sqrt(2.)/3.0/fpi
3310 gamax=gama1*gfun(xmaa**2)/gfun(ama1**2)
3312 hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gamax)
3313 $ *(cmplx(vec1(i))*fpik(xmro1)
3314 $ +cmplx(vec2(i))*fpik(xmro2))
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.
3326 c the statistical factor for identical pi-s was cancelled with
3327 c two, for two modes of a1 decay namelly pi+pi-pi- and pi-pi0pi0
3328 c polarimeter vector in tau rest frame
3330 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3331 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3332 c hv is defined for tau- with gamma=b+hv*pol
3338 c ****************************************************************
3339 c g-
FUNCTION used to inroduce energy dependence in a1 width
3340 c ****************************************************************
3341 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3342 * ,ampiz,ampi,amro,gamro,ama1,gama1
3343 * ,amk,amkz,amkst,gamkst
3345 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3346 * ,ampiz,ampi,amro,gamro,ama1,gama1
3347 * ,amk,amkz,amkst,gamkst
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)
3353 gfun=qkwa*(1.623+10.38/qkwa-9.32/qkwa**2+0.65/qkwa**3)
3356 COMPLEX FUNCTION bwigs(S,M,G)
3357 c **********************************************************
3358 c p-wave breit-wigner for k*
3359 c **********************************************************
3361 REAL pi,pim,qs,qm,w,gs,mk
3363 c ajw: add k*-prim possibility:
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)
3371 c ------------ parameters --------------------
3378 c ajw: add k*-prim possibility:
3384 c ------- breit-wigner -----------------------
3387 IF (s.GT.(pim+mk)**2)
THEN
3389 qs=p(s,pim**2,mk**2)
3390 qm=p(m**2,pim**2,mk**2)
3392 gs=g*(m/w)*(qs/qm)**3
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)
3403 bwigs=m**2/cmplx(m**2-s,-m*gs)
3405 bwigs=m**2/cmplx(m**2-s,-m*gs)
3409 COMPLEX FUNCTION bwig(S,M,G)
3410 c **********************************************************
3411 c p-wave breit-wigner for rho
3412 c **********************************************************
3414 REAL pi,pim,qs,qm,w,gs
3416 c ------------ parameters --------------------
3421 c ------- breit-wigner -----------------------
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)
3427 gs=g*(m/w)*(qs/qm)**3
3431 bwig=m**2/cmplx(m**2-s,-m*gs)
3434 COMPLEX FUNCTION fpik(W)
3435 c **********************************************************
3437 c **********************************************************
3439 REAL rom,rog,rom1,rog1,beta1,pi,pim,s,w
3443 c ------------ parameters --------------------
3444 IF (init.EQ.0 )
THEN
3462 c -----------------------------------------------
3464 fpik= (bwig(s,rom,rog)+beta1*bwig(s,rom1,rog1))
3469 c **********************************************************
3470 c square of pion form factor
3471 c **********************************************************
3473 fpirho=cabs(fpik(w))**2
3475 SUBROUTINE clvec(HJ,PN,PIV)
3476 c ----------------------------------------------------------------------
3477 * calculates the
"VECTOR TYPE" pi-vector piv
3478 * note that the neutrino mom. pn is assumed to be along z-axis
3480 c called by : dampaa
3481 c ----------------------------------------------------------------------
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)))
3489 10 piv(i)=4.*
REAL(hn*conjg(hj(i)))-2.*hh*pn(i)
3492 SUBROUTINE claxi(HJ,PN,PIA)
3493 c ----------------------------------------------------------------------
3494 * calculates the
"AXIAL TYPE" pi-vector pia
3495 * note that the neutrino mom. pn is assumed to be along z-axis
3496 c sign is chosen +/- for decay of tau +/- respectively
3497 c called by : dampaa, clnut
3498 c ----------------------------------------------------------------------
3499 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3500 COMMON / idfc / idff
3502 COMPLEX hj(4),hjc(4)
3503 c det2(i,j)=aimag(hj(i)*hjc(j)-hj(j)*hjc(i))
3504 c -- here was an error(zw, 21.11.1991)
3505 det2(i,j)=aimag(hjc(i)*hj(j)-hjc(j)*hj(i))
3506 c -- it was affecting sign of a_lr asymmetry in a1 decay.
3507 c -- note also collision of notation of gamma_va as defined in
3508 c -- 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)
3515 print *,
'STOP IN CLAXI: KTOM=',ktom
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)
3525 c all four indices are up so pia(3) and pia(4) have same sign
3527 20 pia(i)=pia(i)*sign
3529 SUBROUTINE clnut(HJ,B,HV)
3530 c ----------------------------------------------------------------------
3531 * calculates the contribution by neutrino mass
3532 * note the tau is assumed to be at rest
3534 c called by : dampaa
3535 c ----------------------------------------------------------------------
3541 b=
REAL( HJ(4)*AIMAG(HJ(4)) - HJ(3)*AIMAG(HJ(3))
& - HJ(2)*AIMAG(HJ(2)) - HJ(1)*AIMAG(HJ(1)) )
3544 SUBROUTINE dampog(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3545 c ----------------------------------------------------------------------
3546 * calculates differential cross section and polarimeter vector
3547 * for tau decay into a1, a1 decays next into rho+pi and rho into pi+pi.
3548 * all spin effects in the full decay chain are taken into account.
3549 * calculations done in tau rest frame with z-axis along neutrino moment
3550 * the routine is writen for zero neutrino mass.
3553 c called by : dphtre
3555 c called by : dphsaa
3557 c ----------------------------------------------------------------------
3558 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3559 * ,ampiz,ampi,amro,gamro,ama1,gama1
3560 * ,amk,amkz,amkst,gamkst
3562 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3563 * ,ampiz,ampi,amro,gamro,ama1,gama1
3564 * ,amk,amkz,amkst,gamkst
3565 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3566 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3567 COMMON /testa1/ keya1
3568 REAL hv(4),pt(4),pn(4),pim1(4),pim2(4),pipl(4)
3569 REAL paa(4),vec1(4),vec2(4)
3570 REAL pivec(4),piaks(4),hvm(4)
3571 COMPLEX bwign,hadcur(4),fnorm,formom
3573 * this inline funct. calculates the scalar part of the propagator
3575 c ajwmod to satisfy compiler, comment out this unused function.
3577 bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
3580 * four momentum of a1
3585 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3587 * masses of a1, and of two pi-pairs which may form rho
3588 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3589 xmom =sqrt(abs( (pim2(4)+pipl(4))**2-(pim2(3)+pipl(3))**2
3590 $ -(pim2(2)+pipl(2))**2-(pim2(1)+pipl(1))**2 ))
3591 xmro2 =(pipl(1))**2 +(pipl(2))**2 +(pipl(3))**2
3592 * elements of hadron current
3593 prod1 =vec1(1)*pipl(1)
3594 prod2 =vec2(2)*pipl(2)
3595 p12 =pim1(4)*pim2(4)-pim1(1)*pim2(1)
3596 $ -pim1(2)*pim2(2)-pim1(3)*pim2(3)
3597 p1pl =pim1(4)*pipl(4)-pim1(1)*pipl(1)
3598 $ -pim1(2)*pipl(2)-pim1(3)*pipl(3)
3599 p2pl =pipl(4)*pim2(4)-pipl(1)*pim2(1)
3600 $ -pipl(2)*pim2(2)-pipl(3)*pim2(3)
3602 vec1(i)= (vec1(i)-prod1/xmro2*pipl(i))
3604 gnorm=sqrt(vec1(1)**2+vec1(2)**2+vec1(3)**2)
3606 vec1(i)= vec1(i)/gnorm
3608 vec2(1)=(vec1(2)*pipl(3)-vec1(3)*pipl(2))/sqrt(xmro2)
3609 vec2(2)=(vec1(3)*pipl(1)-vec1(1)*pipl(3))/sqrt(xmro2)
3610 vec2(3)=(vec1(1)*pipl(2)-vec1(2)*pipl(1))/sqrt(xmro2)
3611 p1vec1 =pim1(4)*vec1(4)-pim1(1)*vec1(1)
3612 $ -pim1(2)*vec1(2)-pim1(3)*vec1(3)
3613 p2vec1 =vec1(4)*pim2(4)-vec1(1)*pim2(1)
3614 $ -vec1(2)*pim2(2)-vec1(3)*pim2(3)
3615 p1vec2 =pim1(4)*vec2(4)-pim1(1)*vec2(1)
3616 $ -pim1(2)*vec2(2)-pim1(3)*vec2(3)
3617 p2vec2 =vec2(4)*pim2(4)-vec2(1)*pim2(1)
3618 $ -vec2(2)*pim2(2)-vec2(3)*pim2(3)
3620 fnorm=formom(xmaa,xmom)
3625 hadcur(i) = fnorm *(
3626 $ vec1(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3627 $ -pim2(i)*(p2vec1*p1pl-p1vec1*p2pl)
3628 $ +pipl(i)*(p2vec1*p12 -p1vec1*(ampi**2+p2pl)) )
3630 hadcur(i) = fnorm *(
3631 $ vec2(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3632 $ -pim2(i)*(p2vec2*p1pl-p1vec2*p2pl)
3633 $ +pipl(i)*(p2vec2*p12 -p1vec2*(ampi**2+p2pl)) )
3637 * calculate pi-vectors: vector and axial
3638 CALL clvec(hadcur,pn,pivec)
3639 CALL claxi(hadcur,pn,piaks)
3640 CALL clnut(hadcur,brakm,hvm)
3641 * spin independent part of decay diff-cross-sect. in tau rest frame
3642 brak=brak+(gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3643 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3645 hv(i)=hv(i)-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3646 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3648 c hv is defined for tau- with gamma=b+hv*pol
3650 amplit=(gfermi*ccabib)**2*brak/2.
3651 c the statistical factor for identical pi-s was cancelled with
3652 c two, for two modes of a1 decay namelly pi+pi-pi- and pi-pi0pi0
3653 c polarimeter vector in tau rest frame
3659 SUBROUTINE damppk(MNUM,PT,PN,PIM1,PIM2,PIM3,AMPLIT,HV)
3660 c ----------------------------------------------------------------------
3661 * calculates differential cross section and polarimeter vector
3662 * for tau decay into k k pi, k pi pi.
3663 * all spin effects in the full decay chain are taken into account.
3664 * calculations done in tau rest frame with z-axis along neutrino moment
3665 c mnum decay mode identifier.
3668 c called by : dphtre
3670 c called by : dphsaa
3672 c ----------------------------------------------------------------------
3673 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3674 * ,ampiz,ampi,amro,gamro,ama1,gama1
3675 * ,amk,amkz,amkst,gamkst
3677 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3678 * ,ampiz,ampi,amro,gamro,ama1,gama1
3679 * ,amk,amkz,amkst,gamkst
3680 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3681 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3682 REAL hv(4),pt(4),pn(4),pim1(4),pim2(4),pim3(4)
3683 REAL paa(4),vec1(4),vec2(4),vec3(4),vec4(4),vec5(4)
3684 REAL pivec(4),piaks(4),hvm(4)
3685 REAL fnorm(0:7),coef(1:5,0:7)
3686 COMPLEX hadcur(4),form1,form2,form3,form4,form5,uroj
3688 COMPLEX f1,f2,f3,f4,f5
3690 EXTERNAL form1,form2,form3,form4,form5
3691 DATA pi /3.141592653589793238462643/
3695 IF (icont.EQ.0)
THEN
3703 fnorm(4)=scabib/fpi/dwapi0
3708 coef(1,0)= 2.0*sqrt(2.)/3.0
3709 coef(2,0)=-2.0*sqrt(2.)/3.0
3711 c ajw 2/98: add in the d-wave and i=0 3pi substructure:
3712 coef(3,0)= 2.0*sqrt(2.)/3.0
3719 coef(1,1)=-sqrt(2.)/3.0
3720 coef(2,1)= sqrt(2.)/3.0
3725 coef(1,2)=-sqrt(2.)/3.0
3726 coef(2,2)= sqrt(2.)/3.0
3732 c ajw 11/97: add in the k*-prim-s, ala finkemeier&mirkes
3744 coef(1,4)= 1.0/sqrt(2.)/3.0
3745 coef(2,4)=-1.0/sqrt(2.)/3.0
3750 coef(1,5)=-sqrt(2.)/3.0
3751 coef(2,5)= sqrt(2.)/3.0
3757 c ajw 11/97: add in the k*-prim-s, ala finkemeier&mirkes
3773 coef(5,7)=-sqrt(2.0/3.0)
3778 10 paa(i)=pim1(i)+pim2(i)+pim3(i)
3779 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3780 xmro1 =sqrt(abs((pim3(4)+pim2(4))**2-(pim3(1)+pim2(1))**2
3781 $ -(pim3(2)+pim2(2))**2-(pim3(3)+pim2(3))**2))
3782 xmro2 =sqrt(abs((pim3(4)+pim1(4))**2-(pim3(1)+pim1(1))**2
3783 $ -(pim3(2)+pim1(2))**2-(pim3(3)+pim1(3))**2))
3784 xmro3 =sqrt(abs((pim1(4)+pim2(4))**2-(pim1(1)+pim2(1))**2
3785 $ -(pim1(2)+pim2(2))**2-(pim1(3)+pim2(3))**2))
3786 * elements of hadron current
3787 prod1 =paa(4)*(pim2(4)-pim3(4))-paa(1)*(pim2(1)-pim3(1))
3788 $ -paa(2)*(pim2(2)-pim3(2))-paa(3)*(pim2(3)-pim3(3))
3789 prod2 =paa(4)*(pim3(4)-pim1(4))-paa(1)*(pim3(1)-pim1(1))
3790 $ -paa(2)*(pim3(2)-pim1(2))-paa(3)*(pim3(3)-pim1(3))
3791 prod3 =paa(4)*(pim1(4)-pim2(4))-paa(1)*(pim1(1)-pim2(1))
3792 $ -paa(2)*(pim1(2)-pim2(2))-paa(3)*(pim1(3)-pim2(3))
3794 vec1(i)= pim2(i)-pim3(i) -paa(i)*prod1/xmaa**2
3795 vec2(i)= pim3(i)-pim1(i) -paa(i)*prod2/xmaa**2
3796 vec3(i)= pim1(i)-pim2(i) -paa(i)*prod3/xmaa**2
3797 40 vec4(i)= pim1(i)+pim2(i)+pim3(i)
3798 CALL prod5(pim1,pim2,pim3,vec5)
3800 c be aware that sign of vec2 is opposite to sign of vec1 in a1
case
3802 c rationalize this code:
3803 f1 = cmplx(coef(1,mnum))*form1(mnum,xmaa**2,xmro1**2,xmro2**2)
3804 f2 = cmplx(coef(2,mnum))*form2(mnum,xmaa**2,xmro2**2,xmro1**2)
3805 f3 = cmplx(coef(3,mnum))*form3(mnum,xmaa**2,xmro3**2,xmro1**2)
3807 $cmplx(coef(4,mnum))*form4(mnum,xmaa**2,xmro1**2,xmro2**2,xmro3**2)
3808 f5 = (-1.0)*uroj/4.0/pi**2/fpi**2*
3809 $ cmplx(coef(5,mnum))*form5(mnum,xmaa**2,xmro1**2,xmro2**2)
3812 hadcur(i)= cmplx(fnorm(mnum)) * (
3813 $ cmplx(vec1(i))*f1+cmplx(vec2(i))*f2+cmplx(vec3(i))*f3+
3814 $ cmplx(vec4(i))*f4+cmplx(vec5(i))*f5)
3818 hadcur(i)= cmplx(fnorm(mnum)) * (
3819 $cmplx(vec1(i)*coef(1,mnum))*form1(mnum,xmaa**2,xmro1**2,xmro2**2)+
3820 $cmplx(vec2(i)*coef(2,mnum))*form2(mnum,xmaa**2,xmro2**2,xmro1**2)+
3821 $cmplx(vec3(i)*coef(3,mnum))*form3(mnum,xmaa**2,xmro3**2,xmro1**2)+
3823 $cmplx(vec4(i)*coef(4,mnum))*form4(mnum,xmaa**2,xmro1**2,
3824 $ xmro2**2,xmro3**2) +
3825 $(-1.0)*uroj/4.0/pi**2/fpi**2*
3826 $cmplx(vec5(i)*coef(5,mnum))*form5(mnum,xmaa**2,xmro1**2,xmro2**2))
3830 * calculate pi-vectors: vector and axial
3831 CALL clvec(hadcur,pn,pivec)
3832 CALL claxi(hadcur,pn,piaks)
3833 CALL clnut(hadcur,brakm,hvm)
3834 * spin independent part of decay diff-cross-sect. in tau rest frame
3835 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3836 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3837 amplit=(gfermi)**2*brak/2.
3839 print *,
'MNUM=',mnum
3845 IF (k.EQ.4) znak=1.0
3846 xm1=znak*pim1(k)**2+xm1
3847 xm2=znak*pim2(k)**2+xm2
3848 xm3=znak*pim3(k)**2+xm3
3849 77 print *,
'PIM1=',pim1(k),
'PIM2=',pim2(k),
'PIM3=',pim3(k)
3850 print *,
'XM1=',sqrt(xm1),
'XM2=',sqrt(xm2),
'XM3=',sqrt(xm3)
3851 print *,
'************************************************'
3853 c polarimeter vector in tau rest frame
3855 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3856 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3857 c hv is defined for tau- with gamma=b+hv*pol
3861 SUBROUTINE prod5(P1,P2,P3,PIA)
3862 c ----------------------------------------------------------------------
3863 c
external product of p1, p2, p3 4-momenta.
3864 c sign is chosen +/- for decay of tau +/- respectively
3865 c called by : dampaa, clnut
3866 c ----------------------------------------------------------------------
3867 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3868 COMMON / idfc / idff
3869 REAL pia(4),p1(4),p2(4),p3(4)
3870 det2(i,j)=p1(i)*p2(j)-p2(i)*p1(j)
3871 * -----------------------------------
3872 IF (ktom.EQ.1.OR.ktom.EQ.-1)
THEN
3873 sign= idff/abs(idff)
3874 ELSEIF (ktom.EQ.2)
THEN
3875 sign=-idff/abs(idff)
3877 print *,
'STOP IN PROD5: KTOM=',ktom
3881 c epsilon( p1(1), p2(2), p3(3), (4) ) = 1
3883 pia(1)= -p3(3)*det2(2,4)+p3(4)*det2(2,3)+p3(2)*det2(3,4)
3884 pia(2)= -p3(4)*det2(1,3)+p3(3)*det2(1,4)-p3(1)*det2(3,4)
3885 pia(3)= p3(4)*det2(1,2)-p3(2)*det2(1,4)+p3(1)*det2(2,4)
3886 pia(4)= p3(3)*det2(1,2)-p3(2)*det2(1,3)+p3(1)*det2(2,3)
3887 c all four indices are up so pia(3) and pia(4) have same sign
3889 20 pia(i)=pia(i)*sign
3892 SUBROUTINE dexnew(MODE,ISGN,POL,PNU,PAA,PNPI,JNPI)
3893 c ----------------------------------------------------------------------
3894 * this simulates tau decay in tau rest frame
3895 * into nu a1, next a1 decays into rho pi and finally rho into pi pi.
3896 * output four momenta: pnu tauneutrino,
3898 * paa hadron 4-vector
3899 * pnpi final state particles
3903 * pim1 pion minus(or pi0) 1 (for tau minus)
3904 * pim2 pion minus(or pi0) 2
3905 * pipl pion plus(or pi-)
3906 * (pipl,pim1) form a rho
3908 c ----------------------------------------------------------------------
3909 COMMON / inout / inut,iout
3910 REAL pol(4),hv(4),paa(4),pnu(4),pnpi(4,9),rn(1)
3914 c ===================
3916 CALL dadnew( -1,isgn,hv,pnu,paa,pnpi,jdumm)
3918 cc CALL hbook1(816,
'WEIGHT DISTRIBUTION DEXNEW $',100,-2.,2.)
3920 cc CALL hbook1(816,
'WEIGHT DISTRIBUTION DEXAA $',100,-2.,2.)
3923 ELSEIF(mode.EQ. 0)
THEN
3924 * =======================
3926 IF(iwarm.EQ.0) goto 902
3927 CALL dadnew( 0,isgn,hv,pnu,paa,pnpi,jnpi)
3928 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
3929 cc CALL hfill(816,wt)
3931 IF(rn(1).GT.wt) goto 300
3933 ELSEIF(mode.EQ. 1)
THEN
3934 * =======================
3935 CALL dadnew( 1,isgn,hv,pnu,paa,pnpi,jdumm)
3940 902
WRITE(iout, 9020)
3941 9020
FORMAT(
' ----- DEXNEW: LACK OF INITIALISATION')
3944 SUBROUTINE dadnew(MODE,ISGN,HV,PNU,PWB,PNPI,JNPI)
3945 c ----------------------------------------------------------------------
3946 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3947 * ,ampiz,ampi,amro,gamro,ama1,gama1
3948 * ,amk,amkz,amkst,gamkst
3950 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3951 * ,ampiz,ampi,amro,gamro,ama1,gama1
3952 * ,amk,amkz,amkst,gamkst
3953 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3954 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3955 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
3956 REAL*4 gampmc ,gamper
3959 COMMON / inout / inut,iout
3961 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
3963 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3965 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3968 CHARACTER names(nmode)*31
3970 COMMON / inout / inut,iout
3973 REAL*4 pnu(4),pwb(4),pnpi(4,9),hv(4),hhv(4)
3974 REAL*4 pdum1(4),pdum2(4),pdumi(4,9)
3977 REAL*8 swt(nmode),sswt(nmode)
3978 dimension nevraw(nmode),nevovr(nmode),nevacc(nmode)
3980 DATA pi /3.141592653589793238462643/
3984 c ===================
3985 c -- at the moment only two decay modes of multipions have m. elem
3996 #if defined (CePeCe)
3998 #elif defined (ALEPH)
4001 c for 4pi phase space, need lots more trials at initialization,
4002 c or
use the wtmax determined with many trials for default model:
4004 IF (jnpi.LE.nm4)
THEN
4005 a = pkorb(3,37+jnpi)
4017 ELSEIF(jnpi.LE.nm4)
THEN
4018 CALL dph4pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4019 ELSEIF(jnpi.LE.nm4+nm5)
THEN
4020 CALL dph5pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4021 ELSEIF(jnpi.LE.nm4+nm5+nm6)
THEN
4022 CALL dphnpi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4023 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3)
THEN
4024 inum=jnpi-nm4-nm5-nm6
4025 CALL dphspk(wt,hv,pdum1,pdum2,pdumi,inum)
4026 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2)
THEN
4027 inum=jnpi-nm4-nm5-nm6-nm3
4028 CALL dphsrk(wt,hv,pdum1,pdum2,pdumi,inum)
4032 IF(wt.GT.wtmax(jnpi)/1.2) wtmax(jnpi)=wt*1.2
4034 #if defined (CePeCe)
4035 #elif defined (ALEPH)
4037 c print *,
' DADNEW JNPI,NTRIALS,WTMAX =',jnpi,ntrials,wtmax(jnpi)
4039 c CALL hbook1(801,
'WEIGHT DISTRIBUTION DADNPI $',100,0.,2.,.0)
4040 c print 7004,wtmax(jnpi)
4044 ELSEIF(mode.EQ. 0)
THEN
4045 c =======================
4046 IF(iwarm.EQ.0) goto 902
4051 ELSEIF(jnpi.LE.nm4)
THEN
4052 CALL dph4pi(wt,hhv,pnu,pwb,pnpi,jnpi)
4053 ELSEIF(jnpi.LE.nm4+nm5)
THEN
4054 CALL dph5pi(wt,hhv,pnu,pwb,pnpi,jnpi)
4055 ELSEIF(jnpi.LE.nm4+nm5+nm6)
THEN
4056 CALL dphnpi(wt,hhv,pnu,pwb,pnpi,jnpi)
4057 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3)
THEN
4058 inum=jnpi-nm4-nm5-nm6
4059 CALL dphspk(wt,hhv,pnu,pwb,pnpi,inum)
4060 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2)
THEN
4061 inum=jnpi-nm4-nm5-nm6-nm3
4062 CALL dphsrk(wt,hhv,pnu,pwb,pnpi,inum)
4069 c CALL hfill(801,wt/wtmax(jnpi))
4070 nevraw(jnpi)=nevraw(jnpi)+1
4071 swt(jnpi)=swt(jnpi)+wt
4073 sswt(jnpi)=sswt(jnpi)+wt**2
4076 cc sswt(jnpi)=sswt(jnpi)+wt**2
4077 sswt(jnpi)=sswt(jnpi)+dble(wt)**2
4082 IF(wt.GT.wtmax(jnpi)) nevovr(jnpi)=nevovr(jnpi)+1
4083 IF(rn*wtmax(jnpi).GT.wt) goto 300
4084 c rotations to basic tau rest frame
4085 costhe=-1.+2.*rrr(2)
4088 CALL rotor2(thet,pnu,pnu)
4089 CALL rotor3( phi,pnu,pnu)
4090 CALL rotor2(thet,pwb,pwb)
4091 CALL rotor3( phi,pwb,pwb)
4092 CALL rotor2(thet,hv,hv)
4093 CALL rotor3( phi,hv,hv)
4096 CALL rotor2(thet,pnpi(1,i),pnpi(1,i))
4097 CALL rotor3( phi,pnpi(1,i),pnpi(1,i))
4099 nevacc(jnpi)=nevacc(jnpi)+1
4101 ELSEIF(mode.EQ. 1)
THEN
4102 c =======================
4104 IF(nevraw(jnpi).EQ.0) goto 500
4105 pargam=swt(jnpi)/float(nevraw(jnpi)+1)
4107 IF(nevraw(jnpi).NE.0)
4108 & error=sqrt(sswt(jnpi)/swt(jnpi)**2-1./float(nevraw(jnpi)))
4110 WRITE(iout, 7010) names(jnpi),
4111 & nevraw(jnpi),nevacc(jnpi),nevovr(jnpi),pargam,rat,error
4113 gampmc(8+jnpi-1)=rat
4114 gamper(8+jnpi-1)=error
4115 cam nevdec(8+jnpi-1)=nevacc(jnpi)
4120 7003
FORMAT(///1x,15(5h*****)
4121 $ /,
' *', 25x,
'******** DADNEW INITIALISATION ********',9x,1h*
4123 7004
FORMAT(
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*/)
4125 $ /,1x,15(5h*****)/)
4126 7010
FORMAT(///1x,15(5h*****)
4127 $ /,
' *', 25x,
'******** DADNEW FINAL REPORT ******** ',9x,1h*
4128 $ /,
' *', 25x,
'CHANNEL:',a31 ,9x,1h*
4129 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF DECAYS TOTAL ',9x,1h*
4130 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF DECAYS ACCEPTED ',9x,1h*
4131 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
4132 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH IN GEV UNITS ',9x,1h*
4133 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
4134 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
4135 $ /,1x,15(5h*****)/)
4136 902
WRITE(iout, 9020)
4137 9020
FORMAT(
' ----- DADNEW: LACK OF INITIALISATION')
4139 903
WRITE(iout, 9030) jnpi,mode
4140 9030
FORMAT(
' ----- DADNEW: WRONG JNPI',2i5)
4145 SUBROUTINE dph4pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4146 c ----------------------------------------------------------------------
4148 * it simulates 4pi decay in tau rest frame with
4149 * z-axis along 4pi momentum
4151 * it simulates a1 decay in tau rest frame with
4152 * z-axis along a1 momentum
4154 c ----------------------------------------------------------------------
4155 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4156 * ,ampiz,ampi,amro,gamro,ama1,gama1
4157 * ,amk,amkz,amkst,gamkst
4159 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4160 * ,ampiz,ampi,amro,gamro,ama1,gama1
4161 * ,amk,amkz,amkst,gamkst
4162 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4163 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
4165 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4166 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4168 CHARACTER names(nmode)*31
4171 REAL hv(4),pt(4),pn(4),paa(4),pim1(4),pim2(4),pipl(4),pmult(4,9)
4174 REAL*8 uu,ff,ff1,ff2,ff3,ff4,gg1,gg2,gg3,gg4,rr
4175 DATA pi /3.141592653589793238462643/
4177 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
4178 c amro, gamro is only a
PARAMETER for geting hight efficiency
4180 c three body phase space normalised as in bjorken-drell
4181 c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
4182 phspac=1./2**23/pi**11
4185 c init decay mode jnpi
4186 amp1=dcdmas(idffin(1,jnpi))
4187 amp2=dcdmas(idffin(2,jnpi))
4188 amp3=dcdmas(idffin(3,jnpi))
4189 amp4=dcdmas(idffin(4,jnpi))
4193 #if defined (CePeCe)
4200 #elif defined (CLEO)
4207 c ajw: cant simply change amrop, etc, here
4208 c choice is a by-hand tuning/optimization, no simple relationship
4209 c to actual resonance masses(accd to z.was).
4210 c what matters in the
end is what you put in formf/curr .
4235 CALL choice(100+jnpi,rrdum,ichan,prob1,prob2,prob3,
4238 CALL choice(100+jnpi,rrb,ichan,prob1,prob2,prob3,
4240 $ amrop,gamrop,amrx,gamrx,amrb,gamrb)
4250 * masses of 4, 3 and 2 pi systems
4251 c 3 pi with sampling for resonance
4254 ams1=(amp1+amp2+amp3+amp4)**2
4255 ams2=(amtau-amnuta)**2
4256 alp1=atan((ams1-amrop**2)/amrop/gamrop)
4257 alp2=atan((ams2-amrop**2)/amrop/gamrop)
4258 alp=alp1+rr1*(alp2-alp1)
4259 am4sq =amrop**2+amrop*gamrop*tan(alp)
4262 $ ((am4sq-amrop**2)**2+(amrop*gamrop)**2)/(amrop*gamrop)
4263 phspac=phspac*(alp2-alp1)
4267 ams1=(amp2+amp3+amp4)**2
4269 IF (rrr(9).GT.prez)
THEN
4270 am3sq=ams1+ rr1*(ams2-ams1)
4272 c --- this part of jacobian will be recovered later
4275 * phase space with sampling for omega resonance,
4276 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4277 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4278 alp=alp1+rr1*(alp2-alp1)
4279 am3sq =amrx**2+amrx*gamrx*tan(alp)
4281 c --- this part of the jacobian will be recovered later ---------------
4282 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4290 am2sq=ams1+ rr2*(ams2-ams1)
4292 c --- this part of jacobian will be recovered later
4294 * 2 restframe, define piz and pipl
4296 enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4297 enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4298 ppi= enq1**2-amp3**2
4299 pppi=sqrt(abs(enq1**2-amp3**2))
4301 enq1=(am2sq-amp3**2+amp4**2)/(2*am2)
4302 enq2=(am2sq+amp3**2-amp4**2)/(2*am2)
4303 ppi= enq1**2-amp4**2
4304 pppi=sqrt(abs(enq1**2-amp4**2))
4306 phspac=phspac*(4*pi)*(2*pppi/am2)
4308 * piz momentum in 2 rest frame(piz is 3rd pi)
4310 * piz momentum in 2 rest frame
4312 CALL sphera(pppi,piz)
4315 c pipl momentum in 2 rest frame(pipl is 4th pi)
4317 * pipl momentum in 2 rest frame
4322 * 3 rest frame, define pim1
4330 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4331 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4332 ppi = pr(4)**2-am2**2
4340 pim1(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4342 c --- this part of jacobian will be recovered later
4343 ff3=(4*pi)*(2*pr(3)/am3)
4344 * old pions boosted from 2 rest frame to 3 rest frame
4345 exe=(pr(4)+pr(3))/am2
4346 CALL bostr3(exe,piz,piz)
4347 CALL bostr3(exe,pipl,pipl)
4350 thet =acos(-1.+2*rr3)
4352 CALL rotpol(thet,phi,pipl)
4353 CALL rotpol(thet,phi,pim1)
4354 CALL rotpol(thet,phi,piz)
4355 CALL rotpol(thet,phi,pr)
4357 c 4 rest frame, define pim2
4360 * 4 rest frame, define pim2
4365 pr(4)=1./(2*am4)*(am4**2+am3**2-amp1**2)
4366 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4367 ppi = pr(4)**2-am3**2
4375 pim2(4)=1./(2*am4)*(am4**2-am3**2+amp1**2)
4377 c --- this part of jacobian will be recovered later
4378 ff4=(4*pi)*(2*pr(3)/am4)
4379 * old pions boosted from 3 rest frame to 4 rest frame
4380 exe=(pr(4)+pr(3))/am3
4381 CALL bostr3(exe,piz,piz)
4382 CALL bostr3(exe,pipl,pipl)
4383 CALL bostr3(exe,pim1,pim1)
4386 thet =acos(-1.+2*rr3)
4388 CALL rotpol(thet,phi,pipl)
4389 CALL rotpol(thet,phi,pim1)
4390 CALL rotpol(thet,phi,pim2)
4391 CALL rotpol(thet,phi,piz)
4392 CALL rotpol(thet,phi,pr)
4394 * now to the tau rest frame, define paa and neutrino momenta
4398 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am4**2)
4399 paa(3)= sqrt(abs(paa(4)**2-am4**2))
4400 ppi = paa(4)**2-am4**2
4401 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4402 phsp=phsp*(4*pi)*(2*paa(3)/amtau)
4403 * tau-neutrino momentum
4406 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am4**2)
4408 c zbw 20.12.2002 bug fix
4409 IF(rrr(9).LE.0.5*prez)
THEN
4416 c we include remaining part of the jacobian
4418 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4419 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4421 ams1=(amp1+amp3+amp4)**2
4424 ams2=(sqrt(am3sq)-amp1)**2
4426 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4427 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4430 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4431 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4433 ams1=(amp1+amp3+amp4)**2
4434 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4435 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4436 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4439 ams2=(sqrt(am3sq)-amp1)**2
4441 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4442 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4444 c --- second channel
4445 am3sq=(pim2(4)+piz(4)+pipl(4))**2-(pim2(3)+piz(3)+pipl(3))**2
4446 $ -(pim2(2)+piz(2)+pipl(2))**2-(pim2(1)+piz(1)+pipl(1))**2
4448 ams1=(amp2+amp3+amp4)**2
4449 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4450 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4451 gg1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4454 ams2=(sqrt(am3sq)-amp2)**2
4456 gg3=(4*pi)*(xlam(am2**2,amp2**2,am3sq)/am3sq)
4457 gg4=(4*pi)*(xlam(am3sq,amp1**2,am4**2)/am4**2)
4459 c --- jacobian averaged over the two
4460 IF ( ( (ff+gg)*uu+ff*gg ).GT.0.0d0)
THEN
4461 rr=ff*gg*uu/(0.5*prez*(ff+gg)*uu+(1.0-prez)*ff*gg)
4466 * momenta of the two pi-minus are randomly symmetrised
4477 c momenta of pi0-s are generated uniformly only
IF prez=0.0
4487 * all pions boosted from 4 rest frame to tau rest frame
4488 * z-axis antiparallel to neutrino momentum
4489 exe=(paa(4)+paa(3))/am4
4490 CALL bostr3(exe,piz,piz)
4491 CALL bostr3(exe,pipl,pipl)
4492 CALL bostr3(exe,pim1,pim1)
4493 CALL bostr3(exe,pim2,pim2)
4494 CALL bostr3(exe,pr,pr)
4495 c partial width consists of phase space and amplitude
4496 c check on consistency with dadnpi,
THEN, code breakes uniform pion
4497 c distribution in hadronic system
4499 CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4501 cam assume neutrino mass=0. and sum over final polarisation
4503 c brak= 2*(amtau**2-amx2) * (amtau**2+2.*amx2)
4504 c amplit=ccabib**2*gfermi**2/2. * brak * amx2*sigee(amx2,1)
4506 CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4507 ELSEIF (jnpi.EQ.2)
THEN
4508 CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv)
4511 dgamt=1/(2.*amtau)*amplit*phspac
4526 SUBROUTINE dam4pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,AMPLIT,HV)
4527 c ----------------------------------------------------------------------
4528 * calculates differential cross section and polarimeter vector
4529 * for tau decay into 4 pi modes
4530 * all spin effects in the full decay chain are taken into account.
4531 * calculations done in tau rest frame with z-axis along neutrino moment
4532 c mnum decay mode identifier.
4535 c called by : dph4pi
4537 c called by : dphsaa
4539 c ----------------------------------------------------------------------
4540 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4541 * ,ampiz,ampi,amro,gamro,ama1,gama1
4542 * ,amk,amkz,amkst,gamkst
4544 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4545 * ,ampiz,ampi,amro,gamro,ama1,gama1
4546 * ,amk,amkz,amkst,gamkst
4547 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4548 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
4549 REAL hv(4),pt(4),pn(4),pim1(4),pim2(4),pim3(4),pim4(4)
4550 REAL pivec(4),piaks(4),hvm(4)
4551 COMPLEX hadcur(4),form1,form2,form3,form4,form5
4552 EXTERNAL form1,form2,form3,form4,form5
4553 DATA pi /3.141592653589793238462643/
4557 CALL curr_cleo(mnum,pim1,pim2,pim3,pim4,hadcur)
4559 CALL curr(mnum,pim1,pim2,pim3,pim4,hadcur)
4562 * calculate pi-vectors: vector and axial
4563 CALL clvec(hadcur,pn,pivec)
4564 CALL claxi(hadcur,pn,piaks)
4565 CALL clnut(hadcur,brakm,hvm)
4566 * spin independent part of decay diff-cross-sect. in tau rest frame
4567 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4568 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4569 amplit=(ccabib*gfermi)**2*brak/2.
4570 c polarimeter vector in tau rest frame
4572 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4573 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4574 c hv is defined for tau- with gamma=b+hv*pol
4579 SUBROUTINE dph5pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4580 c ----------------------------------------------------------------------
4581 * it simulates 5pi decay in tau rest frame with
4582 * z-axis along 5pi momentum
4583 c ----------------------------------------------------------------------
4584 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4585 * ,ampiz,ampi,amro,gamro,ama1,gama1
4586 * ,amk,amkz,amkst,gamkst
4588 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4589 * ,ampiz,ampi,amro,gamro,ama1,gama1
4592 * ,amk,amkz,amkst,gamkst
4593 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4594 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
4595 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4597 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4599 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4602 CHARACTER names(nmode)*31
4603 REAL hv(4),pt(4),pn(4),paa(4),pmult(4,9)
4604 REAL*4 pr(4),pi1(4),pi2(4),pi3(4),pi4(4),pi5(4)
4605 REAL*8 amp1,amp2,amp3,amp4,amp5,ams1,ams2,amom,gamom
4606 REAL*8 am5sq,am4sq,am3sq,am2sq,am5,am4,am3
4608 REAL*8 gg1,gg2,gg3,ff1,ff2,ff3,ff4,alp,alp1,alp2
4617 DATA pi /3.141592653589793238462643/
4624 bwign(xm,am,gammab)=xm**2/cmplx(xm**2-am**2,gammab*am)
4626 bwign(xm,am,gamma)=xm**2/cmplx(xm**2-am**2,gamma*am)
4633 c 6 body phase space normalised as in bjorken-drell
4634 c d**3 p /2e/(2pi)**3 (2pi)**4 delta4(sum p)
4635 phspac=1./2**29/pi**14
4636 c phspac=1./2**5/pi**2
4637 c init 5pi decay mode(jnpi)
4638 amp1=dcdmas(idffin(1,jnpi))
4639 amp2=dcdmas(idffin(2,jnpi))
4640 amp3=dcdmas(idffin(3,jnpi))
4641 amp4=dcdmas(idffin(4,jnpi))
4642 amp5=dcdmas(idffin(5,jnpi))
4652 c masses of 5, 4, 3 and 2 pi systems
4653 c 3 pi with sampling for omega resonance
4657 ams1=(amp1+amp2+amp3+amp4+amp5)**2
4658 ams2=(amtau-amnuta)**2
4659 am5sq=ams1+ rr1*(ams2-ams1)
4661 phspac=phspac*(ams2-ams1)
4666 ams1=(amp2+amp3+amp4+amp5)**2
4668 am4sq=ams1+ rr1*(ams2-ams1)
4673 c phase space with sampling for omega resonance
4675 ams1=(amp2+amp3+amp4)**2
4677 alp1=atan((ams1-amom**2)/amom/gamom)
4678 alp2=atan((ams2-amom**2)/amom/gamom)
4679 alp=alp1+rr1*(alp2-alp1)
4680 am3sq =amom**2+amom*gamom*tan(alp)
4682 c --- this part of the jacobian will be recovered later ---------------
4683 gg2=((am3sq-amom**2)**2+(amom*gamom)**2)/(amom*gamom)
4686 c am3sq=ams1+ rr1*(ams2-ams1)
4688 c --- this part of jacobian will be recovered later
4696 am2sq=ams1+ rr2*(ams2-ams1)
4698 c --- this part of jacobian will be recovered later
4701 c(34) restframe, define pi3 and pi4
4702 enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4703 enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4704 ppi= enq1**2-amp3**2
4705 pppi=sqrt(abs(enq1**2-amp3**2))
4706 ff1=(4*pi)*(2*pppi/am2)
4707 c pi3 momentum in(34) rest frame
4708 call sphera(pppi,pi3)
4710 c pi4 momentum in(34) rest frame
4715 c(234) rest frame, define pi2
4719 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4720 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4721 ppi = pr(4)**2-am2**2
4725 pi2(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4727 c --- this part of jacobian will be recovered later
4728 ff2=(4*pi)*(2*pr(3)/am3)
4729 c old pions boosted from 2 rest frame to 3 rest frame
4730 exe=(pr(4)+pr(3))/am2
4731 call bostr3(exe,pi3,pi3)
4732 call bostr3(exe,pi4,pi4)
4735 thet =acos(-1.+2*rr3)
4737 call rotpol(thet,phi,pi2)
4738 call rotpol(thet,phi,pi3)
4739 call rotpol(thet,phi,pi4)
4741 c(2345) rest frame, define pi5
4745 pr(4)=1./(2*am4)*(am4**2+am3**2-amp5**2)
4746 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4747 ppi = pr(4)**2-am3**2
4751 pi5(4)=1./(2*am4)*(am4**2-am3**2+amp5**2)
4753 c --- this part of jacobian will be recovered later
4754 ff3=(4*pi)*(2*pr(3)/am4)
4755 c old pions boosted from 3 rest frame to 4 rest frame
4756 exe=(pr(4)+pr(3))/am3
4757 call bostr3(exe,pi2,pi2)
4758 call bostr3(exe,pi3,pi3)
4759 call bostr3(exe,pi4,pi4)
4762 thet =acos(-1.+2*rr3)
4764 call rotpol(thet,phi,pi2)
4765 call rotpol(thet,phi,pi3)
4766 call rotpol(thet,phi,pi4)
4767 call rotpol(thet,phi,pi5)
4769 c(12345) rest frame, define pi1
4773 pr(4)=1./(2*am5)*(am5**2+am4**2-amp1**2)
4774 pr(3)= sqrt(abs(pr(4)**2-am4**2))
4775 ppi = pr(4)**2-am4**2
4779 pi1(4)=1./(2*am5)*(am5**2-am4**2+amp1**2)
4781 c --- this part of jacobian will be recovered later
4782 ff4=(4*pi)*(2*pr(3)/am5)
4783 c old pions boosted from 4 rest frame to 5 rest frame
4784 exe=(pr(4)+pr(3))/am4
4785 call bostr3(exe,pi2,pi2)
4786 call bostr3(exe,pi3,pi3)
4787 call bostr3(exe,pi4,pi4)
4788 call bostr3(exe,pi5,pi5)
4791 thet =acos(-1.+2*rr3)
4793 call rotpol(thet,phi,pi1)
4794 call rotpol(thet,phi,pi2)
4795 call rotpol(thet,phi,pi3)
4796 call rotpol(thet,phi,pi4)
4797 call rotpol(thet,phi,pi5)
4799 * now to the tau rest frame, define paa and neutrino momenta
4803 c paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5**2)
4804 c paa(3)= sqrt(abs(paa(4)**2-am5**2))
4805 c ppi = paa(4)**2-am5**2
4806 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5sq)
4807 paa(3)= sqrt(abs(paa(4)**2-am5sq))
4808 ppi = paa(4)**2-am5sq
4809 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4810 * tau-neutrino momentum
4813 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am5**2)
4816 phspac=phspac * gg1*gg2*gg3*ff1*ff2*ff3*ff4
4818 c all pions boosted from 5 rest frame to tau rest frame
4819 c z-axis antiparallel to neutrino momentum
4820 exe=(paa(4)+paa(3))/am5
4821 call bostr3(exe,pi1,pi1)
4822 call bostr3(exe,pi2,pi2)
4823 call bostr3(exe,pi3,pi3)
4824 call bostr3(exe,pi4,pi4)
4825 call bostr3(exe,pi5,pi5)
4827 c partial width consists of phase space and amplitude
4828 c amplitude(cf ys.tsai phys.rev.d4,2821(1971)
4829 c or f.gilman sh.rhie phys.rev.d31,1066(1985)
4833 qxn=paa(4)*pn(4)-paa(1)*pn(1)-paa(2)*pn(2)-paa(3)*pn(3)
4834 brak=2*(gv**2+ga**2)*(2*pxq*qxn+am5sq*pxn)
4835 & -6*(gv**2-ga**2)*amtau*amnuta*am5sq
4836 fompp = cabs(bwign(am3,amom,gamom))**2
4837 c normalisation factor(to some numerical undimensioned factor;
4838 c cf r.fischer et al zphys c3, 313 (1980))
4840 c amplit=ccabib**2*gfermi**2/2. * brak * am5sq*sigee(am5sq,jnpi)
4841 amplit=ccabib**2*gfermi**2/2. * brak
4842 amplit = amplit * fompp * fnorm
4844 c amplit = amplit * fnorm
4845 dgamt=1/(2.*amtau)*amplit*phspac
4859 c missing: transposition of identical particles, statistical factors
4860 c for identical matrices, polarimetric vector. matrix element rather nai
4862 c missing: transposition of identical particles, startistical factors
4863 c for identical matrices, polarimetric vector. matrix element rather naive.
4865 c flat phase space in pion system + with breit wigner for omega
4866 c anyway it is better than nothing, and code is improvable.
4868 SUBROUTINE dphsrk(DGAMT,HV,PN,PR,PMULT,INUM)
4869 c ----------------------------------------------------------------------
4870 c it simulates rho decay in tau rest frame with
4871 c z-axis along rho momentum
4872 c rho decays to k kbar
4873 c ----------------------------------------------------------------------
4874 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4875 * ,ampiz,ampi,amro,gamro,ama1,gama1
4876 * ,amk,amkz,amkst,gamkst
4878 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4879 * ,ampiz,ampi,amro,gamro,ama1,gama1
4880 * ,amk,amkz,amkst,gamkst
4881 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4882 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
4883 REAL hv(4),pt(4),pn(4),pr(4),pkc(4),pkz(4),qq(4),pmult(4,9)
4889 DATA pi /3.141592653589793238462643/
4892 c three body phase space normalised as in bjorken-drell
4893 phspac=1./2**11/pi**5
4899 c mass of(real/virtual) rho
4901 ams2=(amtau-amnuta)**2
4904 amx2=ams1+ rr1(1)*(ams2-ams1)
4906 phspac=phspac*(ams2-ams1)
4907 c phase space with sampling for rho resonance
4908 c alp1=atan((ams1-amro**2)/amro/gamro)
4909 c alp2=atan((ams2-amro**2)/amro/gamro)
4912 c CALL ranmar(rr1,1)
4913 c alp=alp1+rr1(1)*(alp2-alp1)
4914 c amx2=amro**2+amro*gamro*tan(alp)
4916 c
IF(amx.LT.(amk+amkz)) go to 100
4918 c phspac=phspac*((amx2-amro**2)**2+(amro*gamro)**2)/(amro*gamro)
4919 c phspac=phspac*(alp2-alp1)
4921 c tau-neutrino momentum
4924 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
4925 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
4929 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
4931 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
4934 enq1=(amx2+amk**2-amkz**2)/(2.*amx)
4935 enq2=(amx2-amk**2+amkz**2)/(2.*amx)
4936 pppi=sqrt((enq1-amk)*(enq1+amk))
4937 phspac=phspac*(4*pi)*(2*pppi/amx)
4938 c charged pi momentum in rho rest frame
4939 CALL sphera(pppi,pkc)
4941 c neutral pi momentum in rho rest frame
4945 exe=(pr(4)+pr(3))/amx
4946 c pions boosted from rho rest frame to tau rest frame
4947 CALL bostr3(exe,pkc,pkc)
4948 CALL bostr3(exe,pkz,pkz)
4950 30 qq(i)=pkc(i)-pkz(i)
4951 c qq transverse to pr
4952 pksd =pr(4)*pr(4)-pr(3)*pr(3)-pr(2)*pr(2)-pr(1)*pr(1)
4953 qqpks=pr(4)* qq(4)-pr(3)* qq(3)-pr(2)* qq(2)-pr(1)* qq(1)
4955 31 qq(i)=qq(i)-pr(i)*qqpks/pksd
4958 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
4960 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
4961 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
4962 & +(gv**2-ga**2)*amtau*amnuta*qq2
4963 amplit=(gfermi*ccabib)**2*brak*2*fpirk(amx)
4964 dgamt=1/(2.*amtau)*amplit*phspac
4966 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
4974 c ----------------------------------------------------------
4975 c square of pion form factor
4976 c ----------------------------------------------------------
4977 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4978 * ,ampiz,ampi,amro,gamro,ama1,gama1
4979 * ,amk,amkz,amkst,gamkst
4981 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4982 * ,ampiz,ampi,amro,gamro,ama1,gama1
4983 * ,amk,amkz,amkst,gamkst
4986 fpirk=cabs(fpikm(w,amk,amkz))**2
4987 c fpirk=cabs(fpikmk(w,amk,amkz))**2
4989 COMPLEX FUNCTION fpikmk(W,XM1,XM2)
4990 c **********************************************************
4992 c **********************************************************
4994 REAL rom,rog,rom1,rog1,beta1,pi,pim,s,w
4998 c ------------ parameters --------------------
4999 IF (init.EQ.0 )
THEN
5010 c -----------------------------------------------
5012 fpikmk=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
5018 c initialize lund
COMMON
5021 parameter(nmxhep=2000)
5022 common/hepevtx/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
5023 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
5028 SUBROUTINE dwrph(KTO,PHX)
5030 c -------------------------
5032 IMPLICIT REAL*8 (a-h,o-z)
5039 c
CASE of tau radiative decays.
5040 c filling of the lund
COMMON block.
5043 IF (qhot(4).GT.1.e-5) CALL dwluph(kto,qhot)
5046 SUBROUTINE dwluph(KTO,PHOT)
5047 c---------------------------------------------------------------------
5048 c lorentz transformation to cmsystem and
5049 c updating of hepevt record
5051 c called by : dexay1,(dekay1,dekay2)
5053 c used when radiative corrections in decays are generated
5054 c---------------------------------------------------------------------
5057 COMMON /taupos/ np1,np2
5063 COMMON /taupos/ np1,np2
5067 IF (phot(4).LE.0.0)
RETURN
5069 c position of decaying particle:
5070 IF((kto.EQ. 1).OR.(kto.EQ.11))
THEN
5077 IF(ktos.GT.10) ktos=ktos-10
5078 c boost and append photon(gamma is 22)
5079 CALL tralo4(ktos,phot,phot,am)
5080 CALL filhep(0,1,22,nps,nps,0,0,phot,0.0,.true.)
5085 SUBROUTINE dwluel(KTO,ISGN,PNU,PWB,PEL,PNE)
5086 c ----------------------------------------------------------------------
5087 c lorentz transformation to cmsystem and
5088 c updating of hepevt record
5090 c isgn = 1/-1 for tau-/tau+
5092 c called by : dexay,(dekay1,dekay2)
5093 c ----------------------------------------------------------------------
5096 COMMON /taupos/ np1,np2
5099 REAL pnu(4),pwb(4),pel(4),pne(4)
5102 COMMON /taupos/ np1,np2
5105 c position of decaying particle:
5112 c tau neutrino(nu_tau is 16)
5113 CALL tralo4(kto,pnu,pnu,am)
5114 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5117 CALL tralo4(kto,pwb,pwb,am)
5118 c CALL filhep(0,2,-24*isgn,nps,nps,0,0,pwb,am,.true.)
5120 c electron(e- is 11)
5121 CALL tralo4(kto,pel,pel,am)
5122 CALL filhep(0,1,11*isgn,nps,nps,0,0,pel,am,.false.)
5124 c anti electron neutrino(nu_e is 12)
5125 CALL tralo4(kto,pne,pne,am)
5126 CALL filhep(0,1,-12*isgn,nps,nps,0,0,pne,am,.true.)
5130 SUBROUTINE dwlumu(KTO,ISGN,PNU,PWB,PMU,PNM)
5131 c ----------------------------------------------------------------------
5132 c lorentz transformation to cmsystem and
5133 c updating of hepevt record
5135 c isgn = 1/-1 for tau-/tau+
5137 c called by : dexay,(dekay1,dekay2)
5138 c ----------------------------------------------------------------------
5141 COMMON /taupos/ np1,np2
5144 REAL pnu(4),pwb(4),pmu(4),pnm(4)
5147 COMMON /taupos/ np1,np2
5150 c position of decaying particle:
5157 c tau neutrino(nu_tau is 16)
5158 CALL tralo4(kto,pnu,pnu,am)
5159 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5162 CALL tralo4(kto,pwb,pwb,am)
5163 c CALL filhep(0,2,-24*isgn,nps,nps,0,0,pwb,am,.true.)
5166 CALL tralo4(kto,pmu,pmu,am)
5167 CALL filhep(0,1,13*isgn,nps,nps,0,0,pmu,am,.false.)
5169 c anti muon neutrino(nu_mu is 14)
5170 CALL tralo4(kto,pnm,pnm,am)
5171 CALL filhep(0,1,-14*isgn,nps,nps,0,0,pnm,am,.true.)
5175 SUBROUTINE dwlupi(KTO,ISGN,PPI,PNU)
5176 c ----------------------------------------------------------------------
5177 c lorentz transformation to cmsystem and
5178 c updating of hepevt record
5180 c isgn = 1/-1 for tau-/tau+
5182 c called by : dexay,(dekay1,dekay2)
5183 c ----------------------------------------------------------------------
5186 COMMON /taupos/ np1,np2
5188 c position of decaying particle:
5195 c tau neutrino(nu_tau is 16)
5196 CALL tralo4(kto,pnu,pnu,am)
5197 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5199 c charged pi meson(pi+ is 211)
5200 CALL tralo4(kto,ppi,ppi,am)
5201 CALL filhep(0,1,-211*isgn,nps,nps,0,0,ppi,am,.true.)
5205 SUBROUTINE dwluro(KTO,ISGN,PNU,PRHO,PIC,PIZ)
5206 c ----------------------------------------------------------------------
5207 c lorentz transformation to cmsystem and
5208 c updating of hepevt record
5210 c isgn = 1/-1 for tau-/tau+
5212 c called by : dexay,(dekay1,dekay2)
5213 c ----------------------------------------------------------------------
5216 COMMON /taupos/ np1,np2
5219 REAL pnu(4),prho(4),pic(4),piz(4)
5222 COMMON /taupos/ np1,np2
5225 c position of decaying particle:
5232 c tau neutrino(nu_tau is 16)
5233 CALL tralo4(kto,pnu,pnu,am)
5234 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5236 c charged rho meson(rho+ is 213)
5237 CALL tralo4(kto,prho,prho,am)
5238 CALL filhep(0,2,-213*isgn,nps,nps,0,0,prho,am,.true.)
5240 c charged pi meson(pi+ is 211)
5241 CALL tralo4(kto,pic,pic,am)
5242 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pic,am,.true.)
5244 c pi0 meson(pi0 is 111)
5245 CALL tralo4(kto,piz,piz,am)
5246 CALL filhep(0,1,111,-2,-2,0,0,piz,am,.true.)
5250 SUBROUTINE dwluaa(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
5251 c ----------------------------------------------------------------------
5252 c lorentz transformation to cmsystem and
5253 c updating of hepevt record
5255 c isgn = 1/-1 for tau-/tau+
5256 c jaa = 1 (2) for a_1- decay to pi+ 2pi- (pi- 2pi0)
5258 c called by : dexay,(dekay1,dekay2)
5259 c ----------------------------------------------------------------------
5262 COMMON /taupos/ np1,np2
5265 REAL pnu(4),paa(4),pim1(4),pim2(4),pipl(4)
5268 COMMON /taupos/ np1,np2
5271 c position of decaying particle:
5278 c tau neutrino(nu_tau is 16)
5279 CALL tralo4(kto,pnu,pnu,am)
5280 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5282 c charged a_1 meson(a_1+ is 20213)
5283 CALL tralo4(kto,paa,paa,am)
5284 CALL filhep(0,1,-20213*isgn,nps,nps,0,0,paa,am,.true.)
5286 c two possible decays of the charged a1 meson
5289 c a1 --> pi+ pi- pi- (or charged conjugate)
5291 c pi minus(or c.c.) (pi+ is 211)
5292 CALL tralo4(kto,pim2,pim2,am)
5293 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pim2,am,.true.)
5295 c pi minus(or c.c.) (pi+ is 211)
5296 CALL tralo4(kto,pim1,pim1,am)
5297 CALL filhep(0,1,-211*isgn,-2,-2,0,0,pim1,am,.true.)
5299 c pi plus(or c.c.) (pi+ is 211)
5300 CALL tralo4(kto,pipl,pipl,am)
5301 CALL filhep(0,1, 211*isgn,-3,-3,0,0,pipl,am,.true.)
5303 ELSE IF (jaa.EQ.2)
THEN
5305 c a1 --> pi- pi0 pi0(or charged conjugate)
5307 c pi zero(pi0 is 111)
5308 CALL tralo4(kto,pim2,pim2,am)
5309 CALL filhep(0,1,111,-1,-1,0,0,pim2,am,.true.)
5311 c pi zero(pi0 is 111)
5312 CALL tralo4(kto,pim1,pim1,am)
5313 CALL filhep(0,1,111,-2,-2,0,0,pim1,am,.true.)
5315 c pi minus(or c.c.) (pi+ is 211)
5316 CALL tralo4(kto,pipl,pipl,am)
5317 CALL filhep(0,1,-211*isgn,-3,-3,0,0,pipl,am,.true.)
5323 SUBROUTINE dwlukk (KTO,ISGN,PKK,PNU)
5324 c ----------------------------------------------------------------------
5325 c lorentz transformation to cmsystem and
5326 c updating of hepevt record
5328 c isgn = 1/-1 for tau-/tau+
5330 c ----------------------------------------------------------------------
5333 COMMON /taupos/ np1,np2
5335 c position of decaying particle
5346 c tau neutrino(nu_tau is 16)
5347 CALL tralo4(kto,pnu,pnu,am)
5348 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5350 c k meson(k+ is 321)
5351 CALL tralo4(kto,pkk,pkk,am)
5352 CALL filhep(0,1,-321*isgn,nps,nps,0,0,pkk,am,.true.)
5356 SUBROUTINE dwluks(KTO,ISGN,PNU,PKS,PKK,PPI,JKST)
5357 COMMON / taukle / bra1,brk0,brk0b,brks
5358 REAL*4 bra1,brk0,brk0b,brks
5360 COMMON /taupos/ np1,np2
5363 c ----------------------------------------------------------------------
5364 c lorentz transformation to cmsystem and
5365 c updating of hepevt record
5367 c isgn = 1/-1 for tau-/tau+
5368 c jkst=10 (20) corresponds to k0b pi- (k- pi0) decay
5370 c ----------------------------------------------------------------------
5373 REAL pnu(4),pks(4),pkk(4),ppi(4)
5375 REAL pnu(4),pks(4),pkk(4),ppi(4),xio(1)
5376 COMMON /taupos/ np1,np2
5379 c position of decaying particle
5386 c tau neutrino(nu_tau is 16)
5387 CALL tralo4(kto,pnu,pnu,am)
5388 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5390 c charged k* meson(k*+ is 323)
5391 CALL tralo4(kto,pks,pks,am)
5392 CALL filhep(0,1,-323*isgn,nps,nps,0,0,pks,am,.true.)
5394 c two possible decay modes of charged k*
5397 c k*- --> pi- k0b(or charged conjugate)
5399 c charged pi meson(pi+ is 211)
5400 CALL tralo4(kto,ppi,ppi,am)
5401 CALL filhep(0,1,-211*isgn,-1,-1,0,0,ppi,am,.true.)
5404 IF (isgn.EQ.-1) bran=brk0
5405 c k0 --> k0_long(is 130) / k0_short(is 310) = 1/1
5407 IF(xio(1).GT.bran)
THEN
5413 CALL tralo4(kto,pkk,pkk,am)
5414 CALL filhep(0,1,k0type,-2,-2,0,0,pkk,am,.true.)
5416 ELSE IF(jkst.EQ.20)
THEN
5420 c pi zero(pi0 is 111)
5421 CALL tralo4(kto,ppi,ppi,am)
5422 CALL filhep(0,1,111,-1,-1,0,0,ppi,am,.true.)
5424 c charged k meson(k+ is 321)
5425 CALL tralo4(kto,pkk,pkk,am)
5426 CALL filhep(0,1,-321*isgn,-2,-2,0,0,pkk,am,.true.)
5432 SUBROUTINE dwlnew(KTO,ISGN,PNU,PWB,PNPI,MODE)
5433 c ----------------------------------------------------------------------
5434 c lorentz transformation to cmsystem and
5435 c updating of hepevt record
5437 c isgn = 1/-1 for tau-/tau+
5439 c called by : dexay,(dekay1,dekay2)
5440 c ----------------------------------------------------------------------
5442 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
5444 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5446 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5449 COMMON /taupos/ np1,np2
5450 CHARACTER names(nmode)*31
5451 REAL pnu(4),pwb(4),pnpi(4,9)
5455 c position of decaying particle
5462 c tau neutrino(nu_tau is 16)
5463 CALL tralo4(kto,pnu,pnu,am)
5464 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5467 CALL tralo4(kto,pwb,pwb,am)
5468 CALL filhep(0,1,-24*isgn,nps,nps,0,0,pwb,am,.true.)
5470 c multi pi mode jnpi
5472 c get multiplicity of mode jnpi
5476 cam kfpi=lunpik(idffin(i,jnpi),-isgn)
5477 kfpi=lunpik(idffin(i,jnpi), isgn)
5479 kfpi=lunpik(idffin(i,jnpi),-isgn)
5481 c for charged conjugate
case, change charged pions only
5482 c
IF(kfpi.NE.111)kfpi=kfpi*isgn
5486 CALL tralo4(kto,ppi,ppi,am)
5487 CALL filhep(0,1,kfpi,-i,-i,0,0,ppi,am,.true.)
5492 #if defined (CePeCe)
5496 c ----------------------------------------------------------------------
5497 c calculates mass of pp(double precision)
5500 c ----------------------------------------------------------------------
5501 IMPLICIT REAL*8 (a-h,o-z)
5503 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5505 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5510 c ******************
5511 c ----------------------------------------------------------------------
5512 c calculates mass of pp
5515 c ----------------------------------------------------------------------
5517 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5518 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5523 c ----------------------------------------------------------------------
5525 c used by : koralz radkor
5526 c ----------------------------------------------------------------------
5527 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5528 DATA pi /3.141592653589793238462643d0/
5530 IF(abs(y).LT.abs(x))
THEN
5532 IF(x.LE.0d0) the=pi-the
5534 the=acos(x/sqrt(x**2+y**2))
5540 c ----------------------------------------------------------------------
5541 * calculates angle in(0,2*pi) range out of x-y
5543 c used by : koralz radkor
5544 c ----------------------------------------------------------------------
5545 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5546 DATA pi /3.141592653589793238462643d0/
5548 IF(abs(y).LT.abs(x))
THEN
5550 IF(x.LE.0d0) the=pi-the
5552 the=acos(x/sqrt(x**2+y**2))
5554 IF(y.LT.0d0) the=2d0*pi-the
5557 SUBROUTINE rotod1(PH1,PVEC,QVEC)
5558 c ----------------------------------------------------------------------
5561 c ----------------------------------------------------------------------
5562 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5563 dimension pvec(4),qvec(4),rvec(4)
5571 qvec(2)= cs*rvec(2)-sn*rvec(3)
5572 qvec(3)= sn*rvec(2)+cs*rvec(3)
5576 SUBROUTINE rotod2(PH1,PVEC,QVEC)
5577 c ----------------------------------------------------------------------
5579 c used by : koralz radkor
5580 c ----------------------------------------------------------------------
5581 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5582 dimension pvec(4),qvec(4),rvec(4)
5589 qvec(1)= cs*rvec(1)+sn*rvec(3)
5591 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5595 SUBROUTINE rotod3(PH1,PVEC,QVEC)
5596 c ----------------------------------------------------------------------
5598 c used by : koralz radkor
5599 c ----------------------------------------------------------------------
5600 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5602 dimension pvec(4),qvec(4),rvec(4)
5608 qvec(1)= cs*rvec(1)-sn*rvec(2)
5609 qvec(2)= sn*rvec(1)+cs*rvec(2)
5613 SUBROUTINE bostr3(EXE,PVEC,QVEC)
5614 c ----------------------------------------------------------------------
5615 c boost along z axis, exe=exp(eta), eta= hiperbolic velocity.
5617 c used by : tauola koralz(?)
5618 c ----------------------------------------------------------------------
5619 REAL*4 pvec(4),qvec(4),rvec(4)
5632 SUBROUTINE bostd3(EXE,PVEC,QVEC)
5633 c ----------------------------------------------------------------------
5634 c boost along z axis, exe=exp(eta), eta= hiperbolic velocity.
5636 c used by : koralz radkor
5637 c ----------------------------------------------------------------------
5638 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5639 dimension pvec(4),qvec(4),rvec(4)
5653 SUBROUTINE rotor1(PH1,PVEC,QVEC)
5654 c ----------------------------------------------------------------------
5657 c ----------------------------------------------------------------------
5658 REAL*4 pvec(4),qvec(4),rvec(4)
5666 qvec(2)= cs*rvec(2)-sn*rvec(3)
5667 qvec(3)= sn*rvec(2)+cs*rvec(3)
5670 SUBROUTINE rotor2(PH1,PVEC,QVEC)
5671 c ----------------------------------------------------------------------
5674 c ----------------------------------------------------------------------
5675 IMPLICIT REAL*4(a-h,o-z)
5676 REAL*4 pvec(4),qvec(4),rvec(4)
5683 qvec(1)= cs*rvec(1)+sn*rvec(3)
5685 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5688 SUBROUTINE rotor3(PHI,PVEC,QVEC)
5689 c ----------------------------------------------------------------------
5692 c ----------------------------------------------------------------------
5693 REAL*4 pvec(4),qvec(4),rvec(4)
5699 qvec(1)= cs*rvec(1)-sn*rvec(2)
5700 qvec(2)= sn*rvec(1)+cs*rvec(2)
5704 SUBROUTINE spherd(R,X)
5705 c ----------------------------------------------------------------------
5706 c generates uniformly three-vector x on sphere of radius r
5707 c double precison version of sphera
5708 c ----------------------------------------------------------------------
5709 REAL*8 r,x(4),pi,costh,sinth
5711 DATA pi /3.141592653589793238462643d0/
5715 sinth=sqrt(1 -costh**2)
5716 x(1)=r*sinth*cos(2*pi*rrr(2))
5717 x(2)=r*sinth*sin(2*pi*rrr(2))
5721 SUBROUTINE rotpox(THET,PHI,PP)
5722 IMPLICIT REAL*8 (a-h,o-z)
5723 c ----------------------------------------------------------------------
5725 c double precison version of rotpol
5729 c ----------------------------------------------------------------------
5732 CALL rotod2(thet,pp,pp)
5733 CALL rotod3( phi,pp,pp)
5736 SUBROUTINE sphera(R,X)
5737 c ----------------------------------------------------------------------
5738 c generates uniformly three-vector x on sphere of radius r
5740 c called by : dphsxx,dadmpi,dadmkk
5741 c ----------------------------------------------------------------------
5744 DATA pi /3.141592653589793238462643/
5748 sinth=sqrt(1.-costh**2)
5749 x(1)=r*sinth*cos(2*pi*rrr(2))
5750 x(2)=r*sinth*sin(2*pi*rrr(2))
5754 SUBROUTINE rotpol(THET,PHI,PP)
5755 c ----------------------------------------------------------------------
5757 c called by : dadmaa,dphsaa
5758 c ----------------------------------------------------------------------
5761 CALL rotor2(thet,pp,pp)
5762 CALL rotor3( phi,pp,pp)
5765 #include "../randg/tauola-random.h"
5768 IMPLICIT REAL*8(a-h,o-z)
5769 cern c304 version 29/07/71 dilog 59 c
5771 IF(x .LT.-1.0) go to 1
5772 IF(x .LE. 0.5) go to 2
5773 IF(x .EQ. 1.0) go to 3
5774 IF(x .LE. 2.0) go to 4
5778 z=z-0.5* log(abs(x))**2
5784 3 dilogt=1.64493406684822
5788 z=1.64493406684822 - log(x)* log(abs(t))
5789 5 y=2.66666666666666 *t+0.66666666666666
5790 b= 0.00000 00000 00001
5791 a=y*b +0.00000 00000 00004
5792 b=y*a-b+0.00000 00000 00011
5793 a=y*b-a+0.00000 00000 00037
5794 b=y*a-b+0.00000 00000 00121
5795 a=y*b-a+0.00000 00000 00398
5796 b=y*a-b+0.00000 00000 01312
5797 a=y*b-a+0.00000 00000 04342
5798 b=y*a-b+0.00000 00000 14437
5799 a=y*b-a+0.00000 00000 48274
5800 b=y*a-b+0.00000 00001 62421
5801 a=y*b-a+0.00000 00005 50291
5802 b=y*a-b+0.00000 00018 79117
5803 a=y*b-a+0.00000 00064 74338
5804 b=y*a-b+0.00000 00225 36705
5805 a=y*b-a+0.00000 00793 87055
5806 b=y*a-b+0.00000 02835 75385
5807 a=y*b-a+0.00000 10299 04264
5808 b=y*a-b+0.00000 38163 29463
5809 a=y*b-a+0.00001 44963 00557
5810 b=y*a-b+0.00005 68178 22718
5811 a=y*b-a+0.00023 20021 96094
5812 b=y*a-b+0.00100 16274 96164
5813 a=y*b-a+0.00468 63619 59447
5814 b=y*a-b+0.02487 93229 24228
5815 a=y*b-a+0.16607 30329 27855
5816 a=y*a-b+1.93506 43008 6996
5819 c=======================================================================
5820 c===================
END of cpc part ====================================
5821 c=======================================================================
5824 c-----------------------------------------------------------------------
5825 c initialize rchl currents
5826 c dummy routine for compatibility with new updates of tauola
5829 c-----------------------------------------------------------------------
5830 SUBROUTINE inirchl(FLAG)
5835 25
FORMAT(1x,
"TAUOLA IniRChL: Fatal error, FLAG=",i2,
" but RChL currents missing")
5836 WRITE(*,*)
" in loaded version of TAUOLA-FORTRAN library."