36 COMMON / taubra / gamprt(500),jlist(500),nchan
40 REAL cumul(500),rrr(1)
42 IF(nchan.LE.0.OR.nchan.GT.500) goto 902
49 IF(rrr(1).LT.cumul(i)/cumul(nchan)) ji=i
54 9020
FORMAT(
' ----- JAKER: WRONG NCHAN')
57 SUBROUTINE dekay(KTO,HX)
84 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
88 COMMON /taupos/ np1,np2
89 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
91 parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
93 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
96 CHARACTER names(nmode)*31
97 COMMON / inout / inut,iout
98 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4),hdum(4)
110 IF (iwarm.EQ.1) x=5/(iwarm-1)
112 WRITE(iout,7001) jak1,jak2
116 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
117 CALL dadmel(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
118 CALL dadmmu(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
119 CALL dadmpi(-1,idum,pdum,pdum1,pdum2)
120 CALL dadmro(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4)
121 CALL dadmaa(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5,jdum)
122 CALL dadmkk(-1,idum,pdum,pdum1,pdum2)
123 CALL dadmks(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,jdum)
124 CALL dadnew(-1,idum,hdum,pdum1,pdum2,pdumx,jdum)
130 ELSEIF(kto.EQ.1)
THEN
134 IF(iwarm.EQ.0) goto 902
140 CALL dekay1(0,h,isgn)
141 ELSEIF(kto.EQ.2)
THEN
145 IF(iwarm.EQ.0) goto 902
151 CALL dekay2(0,h,isgn)
152 ELSEIF(kto.EQ.11)
THEN
157 CALL dekay1(1,h,isgn)
158 ELSEIF(kto.EQ.12)
THEN
163 CALL dekay2(1,h,isgn)
164 ELSEIF(kto.EQ.100)
THEN
166 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
167 CALL dadmel( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
168 CALL dadmmu( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
169 CALL dadmpi( 1,idum,pdum,pdum1,pdum2)
170 CALL dadmro( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4)
171 CALL dadmaa( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5,jdum)
172 CALL dadmkk( 1,idum,pdum,pdum1,pdum2)
173 CALL dadmks( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,jdum)
174 CALL dadnew( 1,idum,hdum,pdum1,pdum2,pdumx,jdum)
175 WRITE(iout,7010) nev1,nev2,nevtot
176 WRITE(iout,7011) (nevdec(i),gampmc(i),gamper(i),i= 1,7)
178 $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
189 7001
FORMAT(///1x,15(5h*****)
191 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.9 ******',9x,1h*,
192 $ /,
' *', 25x,
'***********Sep 2005***************',9x,1h*,
193 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
194 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
195 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
196 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
197 $ /,
' *', 25x,
' Physics initialization by CLEO collab ',9x,1h*,
198 $ /,
' *', 25x,
' see Alain Weinstein www home page: ',9x,1h*,
199 $ /,
' *', 25x,
'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
200 $ /,
' *', 25x,
'/korb_doc.html#files ',9x,1h*,
201 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
202 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
204 $ /,
' *', 25x,
'****DEKAY ROUTINE: INITIALIZATION******',9x,1h*,
205 $ /,
' *',i20 ,5x,
'JAK1 = DECAY MODE TAU+ ',9x,1h*,
206 $ /,
' *',i20 ,5x,
'JAK2 = DECAY MODE TAU- ',9x,1h*,
208 7010
FORMAT(///1x,15(5h*****)
210 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.9 ******',9x,1h*,
211 $ /,
' *', 25x,
'***********Sep 2005***************',9x,1h*,
213 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
214 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
215 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
216 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
217 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
218 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
219 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
220 $ /,
' *', 25x,
'*****DEKAY ROUTINE: FINAL REPORT*******',9x,1h*,
221 $ /,
' *',i20 ,5x,
'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*,
222 $ /,
' *',i20 ,5x,
'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*,
223 $ /,
' *',i20 ,5x,
'NEVTOT = SUM ',9x,1h*,
225 $
' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
227 $ ,i10,2f12.7 ,
' DADMEL ELECTRON ',9x,1h*
228 $ /,
' *',i10,2f12.7 ,
' DADMMU MUON ',9x,1h*
229 $ /,
' *',i10,2f12.7 ,
' DADMPI PION ',9x,1h*
230 $ /,
' *',i10,2f12.7,
' DADMRO RHO (->2PI) ',9x,1h*
231 $ /,
' *',i10,2f12.7,
' DADMAA A1 (->3PI) ',9x,1h*
232 $ /,
' *',i10,2f12.7,
' DADMKK KAON ',9x,1h*
233 $ /,
' *',i10,2f12.7,
' DADMKS K* ',9x,1h*)
235 $ ,i10,2f12.7,a31 ,8x,1h*)
237 $ ,20x,
'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
238 $ /,
' *',20x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
241 9020
FORMAT(
' ----- DEKAY: LACK OF INITIALISATION')
244 9100
FORMAT(
' ----- DEKAY: WRONG VALUE OF KTO ')
247 SUBROUTINE dekay1(IMOD,HH,ISGN)
251 COMMON / decp4 / pp1(4),pp2(4),kf1,kf2
252 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
253 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
254 REAL*4 gampmc ,gamper
257 REAL hv(4),pnu(4),ppi(4)
258 REAL pwb(4),pmu(4),pnm(4)
259 REAL prho(4),pic(4),piz(4)
260 REAL paa(4),pim1(4),pim2(4),pipl(4)
267 IF(jak1.EQ.-1)
RETURN
272 IF(jak1.EQ.0) CALL jaker(jak)
274 CALL dadmel(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
275 ELSEIF(jak.EQ.2)
THEN
276 CALL dadmmu(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
277 ELSEIF(jak.EQ.3)
THEN
278 CALL dadmpi(0, isgn,hv,ppi,pnu)
279 ELSEIF(jak.EQ.4)
THEN
280 CALL dadmro(0, isgn,hv,pnu,prho,pic,piz)
281 ELSEIF(jak.EQ.5)
THEN
282 CALL dadmaa(0, isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
283 ELSEIF(jak.EQ.6)
THEN
284 CALL dadmkk(0, isgn,hv,pkk,pnu)
285 ELSEIF(jak.EQ.7)
THEN
286 CALL dadmks(0, isgn,hv,pnu,pks ,pkk,ppi,jkst)
288 CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
294 ELSEIF(imd.EQ.1)
THEN
298 nevdec(jak)=nevdec(jak)+1
303 CALL dwluel(1,isgn,pnu,pwb,pmu,pnm)
304 CALL dwrph(ktom,phot)
308 ELSEIF(jak.EQ.2)
THEN
309 CALL dwlumu(1,isgn,pnu,pwb,pmu,pnm)
310 CALL dwrph(ktom,phot)
314 ELSEIF(jak.EQ.3)
THEN
315 CALL dwlupi(1,isgn,ppi,pnu)
319 ELSEIF(jak.EQ.4)
THEN
320 CALL dwluro(1,isgn,pnu,prho,pic,piz)
324 ELSEIF(jak.EQ.5)
THEN
325 CALL dwluaa(1,isgn,pnu,paa,pim1,pim2,pipl,jaa)
328 ELSEIF(jak.EQ.6)
THEN
329 CALL dwlukk(1,isgn,pkk,pnu)
332 ELSEIF(jak.EQ.7)
THEN
333 CALL dwluks(1,isgn,pnu,pks,pkk,ppi,jkst)
338 CALL dwlnew(1,isgn,pnu,pwb,pnpi,jak)
346 SUBROUTINE dekay2(IMOD,HH,ISGN)
350 COMMON / decp4 / pp1(4),pp2(4),kf1,kf2
351 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
352 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
353 REAL*4 gampmc ,gamper
356 REAL hv(4),pnu(4),ppi(4)
357 REAL pwb(4),pmu(4),pnm(4)
358 REAL prho(4),pic(4),piz(4)
359 REAL paa(4),pim1(4),pim2(4),pipl(4)
366 IF(jak2.EQ.-1)
RETURN
371 IF(jak2.EQ.0) CALL jaker(jak)
373 CALL dadmel(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
374 ELSEIF(jak.EQ.2)
THEN
375 CALL dadmmu(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
376 ELSEIF(jak.EQ.3)
THEN
377 CALL dadmpi(0, isgn,hv,ppi,pnu)
378 ELSEIF(jak.EQ.4)
THEN
379 CALL dadmro(0, isgn,hv,pnu,prho,pic,piz)
380 ELSEIF(jak.EQ.5)
THEN
381 CALL dadmaa(0, isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
382 ELSEIF(jak.EQ.6)
THEN
383 CALL dadmkk(0, isgn,hv,pkk,pnu)
384 ELSEIF(jak.EQ.7)
THEN
385 CALL dadmks(0, isgn,hv,pnu,pks ,pkk,ppi,jkst)
387 CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
392 ELSEIF(imd.EQ.1)
THEN
396 nevdec(jak)=nevdec(jak)+1
401 CALL dwluel(2,isgn,pnu,pwb,pmu,pnm)
402 CALL dwrph(ktom,phot)
406 ELSEIF(jak.EQ.2)
THEN
407 CALL dwlumu(2,isgn,pnu,pwb,pmu,pnm)
408 CALL dwrph(ktom,phot)
412 ELSEIF(jak.EQ.3)
THEN
413 CALL dwlupi(2,isgn,ppi,pnu)
417 ELSEIF(jak.EQ.4)
THEN
418 CALL dwluro(2,isgn,pnu,prho,pic,piz)
422 ELSEIF(jak.EQ.5)
THEN
423 CALL dwluaa(2,isgn,pnu,paa,pim1,pim2,pipl,jaa)
426 ELSEIF(jak.EQ.6)
THEN
427 CALL dwlukk(2,isgn,pkk,pnu)
430 ELSEIF(jak.EQ.7)
THEN
431 CALL dwluks(2,isgn,pnu,pks,pkk,ppi,jkst)
436 CALL dwlnew(2,isgn,pnu,pwb,pnpi,jak)
444 SUBROUTINE dexay(KTO,POL)
458 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
459 REAL*4 gampmc ,gamper
460 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
462 COMMON /taupos/ np1,np2
463 parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
465 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
468 CHARACTER names(nmode)*31
469 COMMON / inout / inut,iout
471 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4)
485 WRITE(iout, 7001) jak1,jak2
489 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
490 CALL dexel(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
491 CALL dexmu(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
492 CALL dexpi(-1,idum,pdum,pdum1,pdum2)
493 CALL dexro(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4)
494 CALL dexaa(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5,idum)
495 CALL dexkk(-1,idum,pdum,pdum1,pdum2)
496 CALL dexks(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,idum)
497 CALL dexnew(-1,idum,pdum,pdum1,pdum2,pdumi,idum)
503 ELSEIF(kto.EQ.1)
THEN
508 IF(iwarm.EQ.0) goto 902
511 CALL dexay1(kto,jak1,jakp,pol,isgn)
512 ELSEIF(kto.EQ.2)
THEN
517 IF(iwarm.EQ.0) goto 902
518 isgn=-idff/iabs(idff)
520 CALL dexay1(kto,jak2,jakm,pol,isgn)
521 ELSEIF(kto.EQ.100)
THEN
523 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
524 CALL dexel( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
525 CALL dexmu( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
526 CALL dexpi( 1,idum,pdum,pdum1,pdum2)
527 CALL dexro( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4)
528 CALL dexaa( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5,idum)
529 CALL dexkk( 1,idum,pdum,pdum1,pdum2)
530 CALL dexks( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,idum)
531 CALL dexnew( 1,idum,pdum,pdum1,pdum2,pdumi,idum)
532 WRITE(iout,7010) nev1,nev2,nevtot
533 WRITE(iout,7011) (nevdec(i),gampmc(i),gamper(i),i= 1,7)
535 $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
542 7001
FORMAT(///1x,15(5h*****)
544 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.9 ******',9x,1h*,
545 $ /,
' *', 25x,
'*********** Sep 2005***************',9x,1h*,
546 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
547 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
548 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
549 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
550 $ /,
' *', 25x,
' Physics initialization by CLEO collab ',9x,1h*,
551 $ /,
' *', 25x,
' see Alain Weinstein www home page: ',9x,1h*,
552 $ /,
' *', 25x,
'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
553 $ /,
' *', 25x,
'/korb_doc.html#files ',9x,1h*,
554 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
555 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
557 $ /,
' *', 25x,
'*******CERN-TH-6793 NOVEMBER 1992*****',9x,1h*,
558 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
559 $ /,
' *', 25x,
'******DEXAY ROUTINE: INITIALIZATION****',9x,1h*
560 $ /,
' *',i20 ,5x,
'JAK1 = DECAY MODE FERMION1 (TAU+) ',9x,1h*
561 $ /,
' *',i20 ,5x,
'JAK2 = DECAY MODE FERMION2 (TAU-) ',9x,1h*
565 7010
FORMAT(///1x,15(5h*****)
567 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.9 ******',9x,1h*,
568 $ /,
' *', 25x,
'***********Sep 2005***************',9x,1h*,
569 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
570 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
571 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
572 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
573 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
574 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
575 $ /,
' *', 25x,
'*******CERN-TH-6793 NOVEMBER 1992*****',9x,1h*,
576 $ /,
' *', 25x,
'******DEXAY ROUTINE: FINAL REPORT******',9x,1h*
577 $ /,
' *',i20 ,5x,
'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*
578 $ /,
' *',i20 ,5x,
'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*
579 $ /,
' *',i20 ,5x,
'NEVTOT = SUM ',9x,1h*
581 $
' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
583 $ ,i10,2f12.7 ,
' DADMEL ELECTRON ',9x,1h*
584 $ /,
' *',i10,2f12.7 ,
' DADMMU MUON ',9x,1h*
585 $ /,
' *',i10,2f12.7 ,
' DADMPI PION ',9x,1h*
586 $ /,
' *',i10,2f12.7,
' DADMRO RHO (->2PI) ',9x,1h*
587 $ /,
' *',i10,2f12.7,
' DADMAA A1 (->3PI) ',9x,1h*
588 $ /,
' *',i10,2f12.7,
' DADMKK KAON ',9x,1h*
589 $ /,
' *',i10,2f12.7,
' DADMKS K* ',9x,1h*)
591 $ ,i10,2f12.7,a31 ,8x,1h*)
593 $ ,20x,
'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
594 $ /,
' *',20x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
596 902
WRITE(iout, 9020)
597 9020
FORMAT(
' ----- DEXAY: LACK OF INITIALISATION')
599 910
WRITE(iout, 9100)
600 9100
FORMAT(
' ----- DEXAY: WRONG VALUE OF KTO ')
603 SUBROUTINE dexay1(KTO,JAKIN,JAK,POL,ISGN)
609 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
610 REAL*4 gampmc ,gamper
611 COMMON / inout / inut,iout
614 REAL prho(4),pic(4),piz(4)
615 REAL pwb(4),pmu(4),pnm(4)
616 REAL paa(4),pim1(4),pim2(4),pipl(4)
622 IF(jakin.EQ.-1)
RETURN
629 IF(jak.EQ.0) CALL jaker(jak)
632 CALL dexel(0, isgn,polar,pnu,pwb,pmu,pnm,phot)
633 CALL dwluel(kto,isgn,pnu,pwb,pmu,pnm)
634 CALL dwrph(kto,phot )
635 ELSEIF(jak.EQ.2)
THEN
636 CALL dexmu(0, isgn,polar,pnu,pwb,pmu,pnm,phot)
637 CALL dwlumu(kto,isgn,pnu,pwb,pmu,pnm)
638 CALL dwrph(kto,phot )
639 ELSEIF(jak.EQ.3)
THEN
640 CALL dexpi(0, isgn,polar,ppi,pnu)
641 CALL dwlupi(kto,isgn,ppi,pnu)
642 ELSEIF(jak.EQ.4)
THEN
643 CALL dexro(0, isgn,polar,pnu,prho,pic,piz)
644 CALL dwluro(kto,isgn,pnu,prho,pic,piz)
645 ELSEIF(jak.EQ.5)
THEN
646 CALL dexaa(0, isgn,polar,pnu,paa,pim1,pim2,pipl,jaa)
647 CALL dwluaa(kto,isgn,pnu,paa,pim1,pim2,pipl,jaa)
648 ELSEIF(jak.EQ.6)
THEN
649 CALL dexkk(0, isgn,polar,pkk,pnu)
650 CALL dwlukk(kto,isgn,pkk,pnu)
651 ELSEIF(jak.EQ.7)
THEN
652 CALL dexks(0, isgn,polar,pnu,pks,pkk,ppi,jkst)
653 CALL dwluks(kto,isgn,pnu,pks,pkk,ppi,jkst)
656 CALL dexnew(0, isgn,polar,pnu,pwb,pnpi,jnpi)
657 CALL dwlnew(kto,isgn,pnu,pwb,pnpi,jak)
659 nevdec(jak)=nevdec(jak)+1
661 SUBROUTINE dexel(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
668 REAL pol(4),hv(4),pwb(4),pnu(4),q1(4),q2(4),ph(4),rn(1)
674 CALL dadmel( -1,isgn,hv,pnu,pwb,q1,q2,ph)
677 ELSEIF(mode.EQ. 0)
THEN
680 IF(iwarm.EQ.0) goto 902
681 CALL dadmel( 0,isgn,hv,pnu,pwb,q1,q2,ph)
682 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
685 IF(rn(1).GT.wt) goto 300
687 ELSEIF(mode.EQ. 1)
THEN
689 CALL dadmel( 1,isgn,hv,pnu,pwb,q1,q2,ph)
695 9020
FORMAT(
' ----- DEXEL: LACK OF INITIALISATION')
698 SUBROUTINE dexmu(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
707 COMMON / inout / inut,iout
708 REAL pol(4),hv(4),pwb(4),pnu(4),q1(4),q2(4),ph(4),rn(1)
714 CALL dadmmu( -1,isgn,hv,pnu,pwb,q1,q2,ph)
717 ELSEIF(mode.EQ. 0)
THEN
720 IF(iwarm.EQ.0) goto 902
721 CALL dadmmu( 0,isgn,hv,pnu,pwb,q1,q2,ph)
722 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
725 IF(rn(1).GT.wt) goto 300
727 ELSEIF(mode.EQ. 1)
THEN
729 CALL dadmmu( 1,isgn,hv,pnu,pwb,q1,q2,ph)
734 902
WRITE(iout, 9020)
735 9020
FORMAT(
' ----- DEXMU: LACK OF INITIALISATION')
738 SUBROUTINE dadmel(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
743 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
744 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
745 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
746 * ,ampiz,ampi,amro,gamro,ama1,gama1
747 * ,amk,amkz,amkst,gamkst
749 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
750 * ,ampiz,ampi,amro,gamro,ama1,gama1
751 * ,amk,amkz,amkst,gamkst
752 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
753 REAL*4 gampmc ,gamper
757 COMMON / inout / inut,iout
760 REAL hhv(4),hv(4),pwb(4),pnu(4),q1(4),q2(4)
761 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4)
764 DATA pi /3.141592653589793238462643/
777 CALL dphsel(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
778 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
782 ELSEIF(mode.EQ. 0)
THEN
785 IF(iwarm.EQ.0) goto 902
787 CALL dphsel(wt,hv,pnu,pwb,q1,q2,phx)
793 IF(wt.GT.wtmax) nevovr=nevovr+1
794 IF(rn*wtmax.GT.wt) goto 300
801 CALL rotor2(thet,pnu,pnu)
802 CALL rotor3( phi,pnu,pnu)
803 CALL rotor2(thet,pwb,pwb)
804 CALL rotor3( phi,pwb,pwb)
805 CALL rotor2(thet,q1,q1)
806 CALL rotor3( phi,q1,q1)
807 CALL rotor2(thet,q2,q2)
808 CALL rotor3( phi,q2,q2)
809 CALL rotor2(thet,hv,hv)
810 CALL rotor3( phi,hv,hv)
811 CALL rotor2(thet,phx,phx)
812 CALL rotor3( phi,phx,phx)
814 44 hhv(i)=-isgn*hv(i)
817 ELSEIF(mode.EQ. 1)
THEN
819 IF(nevraw.EQ.0)
RETURN
820 pargam=swt/float(nevraw+1)
822 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
824 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
832 7010
FORMAT(///1x,15(5h*****)
833 $ /,
' *', 25x,
'******** DADMEL FINAL REPORT ******** ',9x,1h*
834 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF EL DECAYS TOTAL ',9x,1h*
835 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF EL DECS. ACCEPTED ',9x,1h*
836 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
837 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( ELECTRON) IN GEV UNITS ',9x,1h*
838 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
839 $ /,
' *',f20.9,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
840 $ /,
' *',25x,
'COMPLETE QED CORRECTIONS INCLUDED ',9x,1h*
841 $ /,
' *',25x,
'BUT ONLY V-A CUPLINGS ',9x,1h*
843 902
WRITE(iout, 9020)
844 9020
FORMAT(
' ----- DADMEL: LACK OF INITIALISATION')
847 SUBROUTINE dadmmu(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
849 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
850 * ,ampiz,ampi,amro,gamro,ama1,gama1
851 * ,amk,amkz,amkst,gamkst
853 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
854 * ,ampiz,ampi,amro,gamro,ama1,gama1
855 * ,amk,amkz,amkst,gamkst
856 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
857 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
858 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
859 REAL*4 gampmc ,gamper
860 COMMON / inout / inut,iout
862 REAL hhv(4),hv(4),pnu(4),pwb(4),q1(4),q2(4)
863 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4)
866 DATA pi /3.141592653589793238462643/
879 CALL dphsmu(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
880 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
884 ELSEIF(mode.EQ. 0)
THEN
887 IF(iwarm.EQ.0) goto 902
889 CALL dphsmu(wt,hv,pnu,pwb,q1,q2,phx)
895 IF(wt.GT.wtmax) nevovr=nevovr+1
896 IF(rn*wtmax.GT.wt) goto 300
901 CALL rotor2(thet,pnu,pnu)
902 CALL rotor3( phi,pnu,pnu)
903 CALL rotor2(thet,pwb,pwb)
904 CALL rotor3( phi,pwb,pwb)
905 CALL rotor2(thet,q1,q1)
906 CALL rotor3( phi,q1,q1)
907 CALL rotor2(thet,q2,q2)
908 CALL rotor3( phi,q2,q2)
909 CALL rotor2(thet,hv,hv)
910 CALL rotor3( phi,hv,hv)
911 CALL rotor2(thet,phx,phx)
912 CALL rotor3( phi,phx,phx)
914 44 hhv(i)=-isgn*hv(i)
917 ELSEIF(mode.EQ. 1)
THEN
919 IF(nevraw.EQ.0)
RETURN
920 pargam=swt/float(nevraw+1)
922 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
924 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
932 7010
FORMAT(///1x,15(5h*****)
933 $ /,
' *', 25x,
'******** DADMMU FINAL REPORT ******** ',9x,1h*
934 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF MU DECAYS TOTAL ',9x,1h*
935 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF MU DECS. ACCEPTED ',9x,1h*
936 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
937 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (MU DECAY) IN GEV UNITS ',9x,1h*
938 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
939 $ /,
' *',f20.9,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
940 $ /,
' *',25x,
'COMPLETE QED CORRECTIONS INCLUDED ',9x,1h*
941 $ /,
' *',25x,
'BUT ONLY V-A CUPLINGS ',9x,1h*
943 902
WRITE(iout, 9020)
944 9020
FORMAT(
' ----- DADMMU: LACK OF INITIALISATION')
947 SUBROUTINE dphsel(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
953 REAL*4 hvx(4),paax(4),xax(4),qpx(4),xnx(4)
954 REAL*8 hv(4),ph(4),paa(4),xa(4),qp(4),xn(4)
957 CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
968 SUBROUTINE dphsmu(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
974 REAL*4 hvx(4),paax(4),xax(4),qpx(4),xnx(4)
975 REAL*8 hv(4),ph(4),paa(4),xa(4),qp(4),xn(4)
978 CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
989 SUBROUTINE drcmu(DGAMT,HV,PH,PAA,XA,QP,XN,IELMU)
990 IMPLICIT REAL*8 (a-h,o-z)
995 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
996 * ,ampiz,ampi,amro,gamro,ama1,gama1
997 * ,amk,amkz,amkst,gamkst
999 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1000 * ,ampiz,ampi,amro,gamro,ama1,gama1
1001 * ,amk,amkz,amkst,gamkst
1002 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1003 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1005 COMMON / inout / inut,iout
1006 COMMON / taurad / xk0dec,itdkrc
1008 REAL*8 hv(4),pt(4),ph(4),paa(4),xa(4),qp(4),xn(4)
1012 DATA pi /3.141592653589793238462643d0/
1020 phspac=1./2**17/pi**8
1030 IF (ielmu.EQ.1)
THEN
1037 IF ( itdkrc.EQ.0) prhard=0d0
1039 IF(prsoft.LT.0.1)
THEN
1040 print *,
'ERROR IN DRCMU; PRSOFT=',prsoft
1045 ihard=(rr5.GT.prsoft)
1049 ams1=(amu+amnuta)**2
1052 xl1=log(xk1/2/xk0dec)
1057 phspac=phspac*ams2*xl1*xk
1058 phspac=phspac/prhard
1061 phspac=phspac*2**6*pi**3
1062 phspac=phspac/prsoft
1071 am2sq=ams1+ rr2*(ams2-ams1)
1073 phspac=phspac*(ams2-ams1)
1075 enq1=(am2sq+amnuta**2)/(2*am2)
1076 enq2=(am2sq-amnuta**2)/(2*am2)
1077 ppi= enq1**2-amnuta**2
1078 pppi=sqrt(abs(enq1**2-amnuta**2))
1079 phspac=phspac*(4*pi)*(2*pppi/am2)
1081 CALL spherd(pppi,xn)
1091 pr(4)=1.d0/(2*am3)*(am3**2+am2**2-amu**2)
1092 pr(3)= sqrt(abs(pr(4)**2-am2**2))
1093 ppi = pr(4)**2-am2**2
1097 qp(4)=1.d0/(2*am3)*(am3**2-am2**2+amu**2)
1099 phspac=phspac*(4*pi)*(2*pr(3)/am3)
1101 exe=(pr(4)+pr(3))/am2
1102 CALL bostd3(exe,xn,xn)
1103 CALL bostd3(exe,xa,xa)
1107 eps=4*(amu/amtax)**2
1108 xl1=log((2+eps)/eps)
1110 eta =exp(xl1*rr3+xl0)
1113 phspac=phspac*xl1/2*eta
1115 CALL rotpox(thet,phi,xn)
1116 CALL rotpox(thet,phi,xa)
1117 CALL rotpox(thet,phi,qp)
1118 CALL rotpox(thet,phi,pr)
1124 paa(4)=1/(2*amtax)*(amtax**2+am3**2)
1125 paa(3)= sqrt(abs(paa(4)**2-am3**2))
1126 ppi = paa(4)**2-am3**2
1127 phspac=phspac*(4*pi)*(2*paa(3)/amtax)
1135 exe=(paa(4)+paa(3))/am3
1136 CALL bostd3(exe,xn,xn)
1137 CALL bostd3(exe,xa,xa)
1138 CALL bostd3(exe,qp,qp)
1139 CALL bostd3(exe,pr,pr)
1141 thet =acos(-1.+2*rr3)
1143 CALL rotpox(thet,phi,xn)
1144 CALL rotpox(thet,phi,xa)
1145 CALL rotpox(thet,phi,qp)
1146 CALL rotpox(thet,phi,pr)
1161 CALL dampry(itdkrc,xk0dec,ph,xa,qp,xn,amplit,hv)
1162 dgamt=1/(2.*amtax)*amplit*phspac
1164 SUBROUTINE dampry(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1165 IMPLICIT REAL*8 (a-h,o-z)
1171 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1172 * ,ampiz,ampi,amro,gamro,ama1,gama1
1173 * ,amk,amkz,amkst,gamkst
1175 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1176 * ,ampiz,ampi,amro,gamro,ama1,gama1
1177 * ,amk,amkz,amkst,gamkst
1178 REAL*8 hv(4),qp(4),xn(4),xa(4),xk(4)
1182 IF(xk(4).LT.0.1d0*ak0)
THEN
1183 amplit=thb(itdkrc,qp,xn,xa,ak0,hv)
1185 amplit=sqm2(itdkrc,qp,xn,xa,xk,ak0,hv)
1189 FUNCTION sqm2(ITDKRC,QP,XN,XA,XK,AK0,HV)
1202 IMPLICIT REAL*8(a-h,o-z)
1203 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1204 * ,ampiz,ampi,amro,gamro,ama1,gama1
1205 * ,amk,amkz,amkst,gamkst
1207 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1208 * ,ampiz,ampi,amro,gamro,ama1,gama1
1209 * ,amk,amkz,amkst,gamkst
1210 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1211 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1212 COMMON / qedprm /alfinv,alfpi,xk0
1213 REAL*8 alfinv,alfpi,xk0
1214 REAL*8 qp(4),xn(4),xa(4),xk(4)
1217 REAL*8 s0(3),rxa(3),rxk(3),rqp(3)
1218 DATA pi /3.141592653589793238462643d0/
1224 emass2=qp(4)**2-qp(1)**2-qp(2)**2-qp(3)**2
1232 rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1234 rxk(i)=r(4)*xk(4)-r(1)*xk(1)-r(2)*xk(2)-r(3)*xk(3)
1235 rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
1237 qpxn=qp(4)*xn(4)-qp(1)*xn(1)-qp(2)*xn(2)-qp(3)*xn(3)
1238 qpxa=qp(4)*xa(4)-qp(1)*xa(1)-qp(2)*xa(2)-qp(3)*xa(3)
1239 qpxk=qp(4)*xk(4)-qp(1)*xk(1)-qp(2)*xk(2)-qp(3)*xk(3)
1241 xnxk=xn(4)*xk(4)-xn(1)*xk(1)-xn(2)*xk(2)-xn(3)*xk(3)
1242 xaxk=xa(4)*xk(4)-xa(1)*xk(1)-xa(2)*xk(2)-xa(3)*xk(3)
1252 s1= qpxn*txa*( -emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
1254 $qpxn/txk**2* ( tmass2*xaxk - txa*txk+ xaxk*txk) -
1255 $txa*txn/txk - qpxn/(qpxk*txk)* (tqp*xaxk-txk*qpxa)
1256 const4=256*pi/alphai*gf**2
1257 IF (itdkrc.EQ.0) const4=0d0
1260 s0(i) = qpxn*rxa(i)*(-emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
1262 $ qpxn/txk**2* (tmass2*xaxk - txa*rxk(i)+ xaxk*rxk(i))-
1263 $ rxa(i)*txn/txk - qpxn/(qpxk*txk)*(rqp(i)*xaxk- rxk(i)*qpxa)
1264 5 hv(i)=s0(i)/s1-1.d0
1267 FUNCTION thb(ITDKRC,QP,XN,XA,AK0,HV)
1281 IMPLICIT REAL*8(a-h,o-z)
1282 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1283 * ,ampiz,ampi,amro,gamro,ama1,gama1
1284 * ,amk,amkz,amkst,gamkst
1286 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1287 * ,ampiz,ampi,amro,gamro,ama1,gama1
1288 * ,amk,amkz,amkst,gamkst
1289 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1290 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1291 COMMON / qedprm /alfinv,alfpi,xk0
1292 REAL*8 alfinv,alfpi,xk0
1293 dimension qp(4),xn(4),xa(4)
1296 REAL*8 rxa(3),rxn(3),rqp(3)
1297 REAL*8 bornpl(3),am3pol(3),xm3pol(3)
1298 DATA pi /3.141592653589793238462643d0/
1311 rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1312 rxn(i)=r(4)*xn(4)-r(1)*xn(1)-r(2)*xn(2)-r(3)*xn(3)
1314 rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
1318 u3=sqrt(qp(1)**2+qp(2)**2+qp(3)**2)/tmass
1320 w0=(xn(4)+xa(4))/tmass
1332 f0=2*u0/u3*( dilogt(1-(um*wm/(up*wp)))- dilogt(1-wm/wp) +
1333 $dilogt(1-um/up) -2*yu+ 2*log(up)*(yw+yu) ) +
1334 $1/y* ( 2*u3*yu + (1-eps2- 2*y)*log(eps) ) +
1335 $ 2 - 4*(u0/u3*yu -1)* log(2*al)
1336 fp= yu/(2*u3)*(1 + (1-eps2)/y ) + log(eps)/y
1337 fm= yu/(2*u3)*(1 - (1-eps2)/y ) - log(eps)/y
1340 qpxn=qp(4)*xn(4)-qp(1)*xn(1)-qp(2)*xn(2)-qp(3)*xn(3)
1341 qpxa=qp(4)*xa(4)-qp(1)*xa(1)-qp(2)*xa(2)-qp(3)*xa(3)
1342 xnxa=xn(4)*xa(4)-xn(1)*xa(1)-xn(2)*xa(2)-xn(3)*xa(3)
1347 const3=1/(2*alphai*pi)*64*gf**2
1348 IF (itdkrc.EQ.0) const3=0d0
1349 xm3= -( f0* qpxn*txa + fp*eps2* txn*txa +
1350 $fm* qpxn*qpxa + f3* tmass2*xnxa )
1353 brak= (gv+ga)**2*tqp*xnxa+(gv-ga)**2*txa*qpxn
1354 & -(gv**2-ga**2)*tmass*amnuta*qpxa
1355 born= 32*(gfermi**2/2.)*brak
1357 xm3pol(i)= -( f0* qpxn*rxa(i) + fp*eps2* txn*rxa(i) +
1358 $ fm* qpxn* (qpxa + (rxa(i)*tqp-txa*rqp(i))/tmass2 ) +
1359 $ f3* (tmass2*xnxa +txn*rxa(i) -rxn(i)*txa) )
1360 am3pol(i)=xm3pol(i)*const3
1363 & (gv+ga)**2*tmass*xnxa*qp(i)
1364 & -(gv-ga)**2*tmass*qpxn*xa(i)
1365 & +(gv**2-ga**2)*amnuta*txa*qp(i)
1366 & -(gv**2-ga**2)*amnuta*tqp*xa(i) )*
1368 5 hv(i)=(bornpl(i)+am3pol(i))/(born+am3)-1.d0
1370 IF (thb/born.LT.0.1d0)
THEN
1371 print *,
'ERROR IN THB, THB/BORN=',thb/born
1378 SUBROUTINE dexpi(MODE,ISGN,POL,PPI,PNU)
1385 REAL pol(4),hv(4),pnu(4),ppi(4),rn(1)
1389 CALL dadmpi(-1,isgn,hv,ppi,pnu)
1392 ELSEIF(mode.EQ. 0)
THEN
1395 CALL dadmpi( 0,isgn,hv,ppi,pnu)
1396 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1399 IF(rn(1).GT.wt) goto 300
1401 ELSEIF(mode.EQ. 1)
THEN
1403 CALL dadmpi( 1,isgn,hv,ppi,pnu)
1409 SUBROUTINE dadmpi(MODE,ISGN,HV,PPI,PNU)
1411 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1412 * ,ampiz,ampi,amro,gamro,ama1,gama1
1413 * ,amk,amkz,amkst,gamkst
1415 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1416 * ,ampiz,ampi,amro,gamro,ama1,gama1
1417 * ,amk,amkz,amkst,gamkst
1418 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1419 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1420 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
1421 REAL*4 gampmc ,gamper
1422 COMMON / inout / inut,iout
1423 REAL ppi(4),pnu(4),hv(4)
1424 DATA pi /3.141592653589793238462643/
1429 ELSEIF(mode.EQ. 0)
THEN
1432 epi= (amtau**2+ampi**2-amnuta**2)/(2*amtau)
1433 enu= (amtau**2-ampi**2+amnuta**2)/(2*amtau)
1434 xpi= sqrt(epi**2-ampi**2)
1436 CALL sphera(xpi,ppi)
1444 qxn=ppi(4)*pnu(4)-ppi(1)*pnu(1)-ppi(2)*pnu(2)-ppi(3)*pnu(3)
1445 brak=(gv**2+ga**2)*(2*pxq*qxn-ampi**2*pxn)
1446 & +(gv**2-ga**2)*amtau*amnuta*ampi**2
1448 40 hv(i)=-isgn*2*ga*gv*amtau*(2*ppi(i)*qxn-pnu(i)*ampi**2)/brak
1451 ELSEIF(mode.EQ. 1)
THEN
1453 IF(nevtot.EQ.0)
RETURN
1459 gamm=(gfermi*fpi)**2/(16.*pi)*amtau**3*
1461 $ sqrt((amtau**2-ampi**2-amnuta**2)**2
1462 $ -4*ampi**2*amnuta**2 )/amtau**2
1465 WRITE(iout, 7010) nevtot,gamm,rat,error
1472 7010
FORMAT(///1x,15(5h*****)
1473 $ /,
' *', 25x,
'******** DADMPI FINAL REPORT ******** ',9x,1h*
1474 $ /,
' *',i20 ,5x,
'NEVTOT = NO. OF PI DECAYS TOTAL ',9x,1h*
1475 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( PI DECAY) IN GEV UNITS ',9x,1h*
1476 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1477 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9x,1h*
1478 $ /,1x,15(5h*****)/)
1480 SUBROUTINE dexro(MODE,ISGN,POL,PNU,PRO,PIC,PIZ)
1489 COMMON / inout / inut,iout
1490 REAL pol(4),hv(4),pro(4),pnu(4),pic(4),piz(4),rn(1)
1496 CALL dadmro( -1,isgn,hv,pnu,pro,pic,piz)
1500 ELSEIF(mode.EQ. 0)
THEN
1503 IF(iwarm.EQ.0) goto 902
1504 CALL dadmro( 0,isgn,hv,pnu,pro,pic,piz)
1505 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1510 IF(rn(1).GT.wt) goto 300
1512 ELSEIF(mode.EQ. 1)
THEN
1514 CALL dadmro( 1,isgn,hv,pnu,pro,pic,piz)
1520 902
WRITE(iout, 9020)
1521 9020
FORMAT(
' ----- DEXRO: LACK OF INITIALISATION')
1524 SUBROUTINE dadmro(MODE,ISGN,HHV,PNU,PRO,PIC,PIZ)
1526 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1527 * ,ampiz,ampi,amro,gamro,ama1,gama1
1528 * ,amk,amkz,amkst,gamkst
1530 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1531 * ,ampiz,ampi,amro,gamro,ama1,gama1
1532 * ,amk,amkz,amkst,gamkst
1533 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1534 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1535 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
1536 REAL*4 gampmc ,gamper
1537 COMMON / inout / inut,iout
1539 REAL hv(4),pro(4),pnu(4),pic(4),piz(4)
1540 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4)
1543 DATA pi /3.141592653589793238462643/
1556 CALL dphsro(wt,hv,pdum1,pdum2,pdum3,pdum4)
1557 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1562 ELSEIF(mode.EQ. 0)
THEN
1565 IF(iwarm.EQ.0) goto 902
1566 CALL dphsro(wt,hv,pnu,pro,pic,piz)
1573 IF(wt.GT.wtmax) nevovr=nevovr+1
1574 IF(rn*wtmax.GT.wt) goto 300
1576 costhe=-1.+2.*rrr(2)
1579 CALL rotor2(thet,pnu,pnu)
1580 CALL rotor3( phi,pnu,pnu)
1581 CALL rotor2(thet,pro,pro)
1582 CALL rotor3( phi,pro,pro)
1583 CALL rotor2(thet,pic,pic)
1584 CALL rotor3( phi,pic,pic)
1585 CALL rotor2(thet,piz,piz)
1586 CALL rotor3( phi,piz,piz)
1587 CALL rotor2(thet,hv,hv)
1588 CALL rotor3( phi,hv,hv)
1590 44 hhv(i)=-isgn*hv(i)
1593 ELSEIF(mode.EQ. 1)
THEN
1595 IF(nevraw.EQ.0)
RETURN
1596 pargam=swt/float(nevraw+1)
1598 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1600 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1608 7003
FORMAT(///1x,15(5h*****)
1609 $ /,
' *', 25x,
'******** DADMRO INITIALISATION ********',9x,1h*
1610 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
1611 $ /,1x,15(5h*****)/)
1612 7010
FORMAT(///1x,15(5h*****)
1613 $ /,
' *', 25x,
'******** DADMRO FINAL REPORT ******** ',9x,1h*
1614 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF RHO DECAYS TOTAL ',9x,1h*
1615 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF RHO DECS. ACCEPTED ',9x,1h*
1616 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
1617 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (RHO DECAY) IN GEV UNITS ',9x,1h*
1618 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1619 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
1620 $ /,1x,15(5h*****)/)
1621 902
WRITE(iout, 9020)
1622 9020
FORMAT(
' ----- DADMRO: LACK OF INITIALISATION')
1625 SUBROUTINE dphsro(DGAMT,HV,PN,PR,PIC,PIZ)
1630 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1631 * ,ampiz,ampi,amro,gamro,ama1,gama1
1632 * ,amk,amkz,amkst,gamkst
1634 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1635 * ,ampiz,ampi,amro,gamro,ama1,gama1
1636 * ,amk,amkz,amkst,gamkst
1637 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1638 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1639 REAL hv(4),pt(4),pn(4),pr(4),pic(4),piz(4),qq(4),rr1(1)
1640 DATA pi /3.141592653589793238462643/
1644 phspac=1./2**11/pi**5
1651 ams1=(ampi+ampiz)**2
1652 ams2=(amtau-amnuta)**2
1660 alp1=atan((ams1-amro**2)/amro/gamro)
1661 alp2=atan((ams2-amro**2)/amro/gamro)
1665 alp=alp1+rr1(1)*(alp2-alp1)
1666 amx2=amro**2+amro*gamro*tan(alp)
1668 IF(amx.LT.2.*ampi) go to 100
1670 phspac=phspac*((amx2-amro**2)**2+(amro*gamro)**2)/(amro*gamro)
1671 phspac=phspac*(alp2-alp1)
1676 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
1677 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
1681 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
1683 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
1686 enq1=(amx2+ampi**2-ampiz**2)/(2.*amx)
1687 enq2=(amx2-ampi**2+ampiz**2)/(2.*amx)
1688 pppi=sqrt((enq1-ampi)*(enq1+ampi))
1689 phspac=phspac*(4*pi)*(2*pppi/amx)
1691 CALL sphera(pppi,pic)
1697 exe=(pr(4)+pr(3))/amx
1699 CALL bostr3(exe,pic,pic)
1700 CALL bostr3(exe,piz,piz)
1702 30 qq(i)=pic(i)-piz(i)
1705 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
1707 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
1708 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
1709 & +(gv**2-ga**2)*amtau*amnuta*qq2
1710 amplit=(gfermi*ccabib)**2*brak*2*fpirho(amx)
1711 dgamt=1/(2.*amtau)*amplit*phspac
1713 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
1716 SUBROUTINE dexaa(MODE,ISGN,POL,PNU,PAA,PIM1,PIM2,PIPL,JAA)
1727 COMMON / inout / inut,iout
1728 REAL pol(4),hv(4),paa(4),pnu(4),pim1(4),pim2(4),pipl(4),rn(1)
1734 CALL dadmaa( -1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1737 ELSEIF(mode.EQ. 0)
THEN
1740 IF(iwarm.EQ.0) goto 902
1741 CALL dadmaa( 0,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1742 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1745 IF(rn(1).GT.wt) goto 300
1747 ELSEIF(mode.EQ. 1)
THEN
1749 CALL dadmaa( 1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1754 902
WRITE(iout, 9020)
1755 9020
FORMAT(
' ----- DEXAA: LACK OF INITIALISATION')
1758 SUBROUTINE dadmaa(MODE,ISGN,HHV,PNU,PAA,PIM1,PIM2,PIPL,JAA)
1762 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1763 * ,ampiz,ampi,amro,gamro,ama1,gama1
1764 * ,amk,amkz,amkst,gamkst
1766 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1767 * ,ampiz,ampi,amro,gamro,ama1,gama1
1768 * ,amk,amkz,amkst,gamkst
1769 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1770 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1771 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
1772 REAL*4 gampmc ,gamper
1773 COMMON / inout / inut,iout
1775 REAL hv(4),paa(4),pnu(4),pim1(4),pim2(4),pipl(4)
1776 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4),pdum5(4)
1779 DATA pi /3.141592653589793238462643/
1792 CALL dphsaa(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5,jaa)
1793 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1797 ELSEIF(mode.EQ. 0)
THEN
1800 IF(iwarm.EQ.0) goto 902
1801 CALL dphsaa(wt,hv,pnu,paa,pim1,pim2,pipl,jaa)
1808 sswt=sswt+dble(wt)**2
1813 IF(wt.GT.wtmax) nevovr=nevovr+1
1814 IF(rn*wtmax.GT.wt) goto 300
1816 costhe=-1.+2.*rrr(2)
1819 CALL rotpol(thet,phi,pnu)
1820 CALL rotpol(thet,phi,paa)
1821 CALL rotpol(thet,phi,pim1)
1822 CALL rotpol(thet,phi,pim2)
1823 CALL rotpol(thet,phi,pipl)
1824 CALL rotpol(thet,phi,hv)
1826 44 hhv(i)=-isgn*hv(i)
1829 ELSEIF(mode.EQ. 1)
THEN
1831 IF(nevraw.EQ.0)
RETURN
1832 pargam=swt/float(nevraw+1)
1834 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1836 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1844 7003
FORMAT(///1x,15(5h*****)
1845 $ /,
' *', 25x,
'******** DADMAA INITIALISATION ********',9x,1h*
1846 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
1847 $ /,1x,15(5h*****)/)
1848 7010
FORMAT(///1x,15(5h*****)
1849 $ /,
' *', 25x,
'******** DADMAA FINAL REPORT ******** ',9x,1h*
1850 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF A1 DECAYS TOTAL ',9x,1h*
1851 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF A1 DECS. ACCEPTED ',9x,1h*
1852 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
1853 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (A1 DECAY) IN GEV UNITS ',9x,1h*
1854 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1855 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
1856 $ /,1x,15(5h*****)/)
1857 902
WRITE(iout, 9020)
1858 9020
FORMAT(
' ----- DADMAA: LACK OF INITIALISATION')
1861 SUBROUTINE dphsaa(DGAMT,HV,PN,PAA,PIM1,PIM2,PIPL,JAA)
1866 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1867 * ,ampiz,ampi,amro,gamro,ama1,gama1
1868 * ,amk,amkz,amkst,gamkst
1870 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1871 * ,ampiz,ampi,amro,gamro,ama1,gama1
1872 * ,amk,amkz,amkst,gamkst
1873 COMMON / taukle / bra1,brk0,brk0b,brks
1874 REAL*4 bra1,brk0,brk0b,brks
1875 REAL hv(4),pn(4),paa(4),pim1(4),pim2(4),pipl(4)
1885 IF (rmod.LT.bra1)
THEN
1897 $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
1899 SUBROUTINE dexkk(MODE,ISGN,POL,PKK,PNU)
1906 REAL pol(4),hv(4),pnu(4),pkk(4),rn(1)
1910 CALL dadmkk(-1,isgn,hv,pkk,pnu)
1913 ELSEIF(mode.EQ. 0)
THEN
1916 CALL dadmkk( 0,isgn,hv,pkk,pnu)
1917 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1920 IF(rn(1).GT.wt) goto 300
1922 ELSEIF(mode.EQ. 1)
THEN
1924 CALL dadmkk( 1,isgn,hv,pkk,pnu)
1930 SUBROUTINE dadmkk(MODE,ISGN,HV,PKK,PNU)
1934 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1935 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
1937 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1938 * ,ampiz,ampi,amro,gamro,ama1,gama1
1939 * ,amk,amkz,amkst,gamkst
1941 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1942 * ,ampiz,ampi,amro,gamro,ama1,gama1
1943 * ,amk,amkz,amkst,gamkst
1946 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
1947 REAL*4 gampmc ,gamper
1948 COMMON / inout / inut,iout
1949 REAL pkk(4),pnu(4),hv(4)
1950 DATA pi /3.141592653589793238462643/
1955 ELSEIF(mode.EQ. 0)
THEN
1958 ekk= (amtau**2+amk**2-amnuta**2)/(2*amtau)
1959 enu= (amtau**2-amk**2+amnuta**2)/(2*amtau)
1960 xkk= sqrt(ekk**2-amk**2)
1962 CALL sphera(xkk,pkk)
1970 qxn=pkk(4)*pnu(4)-pkk(1)*pnu(1)-pkk(2)*pnu(2)-pkk(3)*pnu(3)
1971 brak=(gv**2+ga**2)*(2*pxq*qxn-amk**2*pxn)
1972 & +(gv**2-ga**2)*amtau*amnuta*amk**2
1974 40 hv(i)=-isgn*2*ga*gv*amtau*(2*pkk(i)*qxn-pnu(i)*amk**2)/brak
1977 ELSEIF(mode.EQ. 1)
THEN
1979 IF(nevtot.EQ.0)
RETURN
1986 gamm=(gfermi*fkk)**2/(16.*pi)*amtau**3*
1988 $ sqrt((amtau**2-amk**2-amnuta**2)**2
1989 $ -4*amk**2*amnuta**2 )/amtau**2
1994 WRITE(iout, 7010) nevtot,gamm,rat,error
2001 7010
FORMAT(///1x,15(5h*****)
2002 $ /,
' *', 25x,
'******** DADMKK FINAL REPORT ********',9x,1h*
2003 $ /,
' *',i20 ,5x,
'NEVTOT = NO. OF K DECAYS TOTAL ',9x,1h*,
2004 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( K DECAY) IN GEV UNITS ',9x,1h*,
2005 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2006 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9x,1h*
2007 $ /,1x,15(5h*****)/)
2009 SUBROUTINE dexks(MODE,ISGN,POL,PNU,PKS,PKK,PPI,JKST)
2021 COMMON / inout / inut,iout
2022 REAL pol(4),hv(4),pks(4),pnu(4),pkk(4),ppi(4),rn(1)
2029 CALL dadmks( -1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2033 ELSEIF(mode.EQ. 0)
THEN
2036 IF(iwarm.EQ.0) goto 902
2037 CALL dadmks( 0,isgn,hv,pnu,pks,pkk,ppi,jkst)
2038 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
2043 IF(rn(1).GT.wt) goto 300
2045 ELSEIF(mode.EQ. 1)
THEN
2047 CALL dadmks( 1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2053 902
WRITE(iout, 9020)
2054 9020
FORMAT(
' ----- DEXKS: LACK OF INITIALISATION')
2057 SUBROUTINE dadmks(MODE,ISGN,HHV,PNU,PKS,PKK,PPI,JKST)
2059 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2060 * ,ampiz,ampi,amro,gamro,ama1,gama1
2061 * ,amk,amkz,amkst,gamkst
2063 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2064 * ,ampiz,ampi,amro,gamro,ama1,gama1
2065 * ,amk,amkz,amkst,gamkst
2066 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2067 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2068 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
2069 REAL*4 gampmc ,gamper
2070 COMMON / taukle / bra1,brk0,brk0b,brks
2071 REAL*4 bra1,brk0,brk0b,brks
2072 COMMON / inout / inut,iout
2074 REAL hv(4),pks(4),pnu(4),pkk(4),ppi(4)
2075 REAL pdum1(4),pdum2(4),pdum3(4),pdum4(4)
2076 REAL*4 rrr(3),rmod(1)
2078 DATA pi /3.141592653589793238462643/
2093 CALL dphsks(wt,hv,pdum1,pdum2,pdum3,pdum4,jkst)
2094 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
2099 ELSEIF(mode.EQ. 0)
THEN
2101 IF(iwarm.EQ.0) goto 902
2107 IF(rmod(1).LT.dec1)
THEN
2112 CALL dphsks(wt,hv,pnu,pks,pkk,ppi,jkst)
2115 IF(wt.GT.wtmax) nevovr=nevovr+1
2119 IF(rn*wtmax.GT.wt) goto 400
2121 costhe=-1.+2.*rrr(2)
2124 CALL rotor2(thet,pnu,pnu)
2125 CALL rotor3( phi,pnu,pnu)
2126 CALL rotor2(thet,pks,pks)
2127 CALL rotor3( phi,pks,pks)
2128 CALL rotor2(thet,pkk,pkk)
2129 CALL rotor3(phi,pkk,pkk)
2130 CALL rotor2(thet,ppi,ppi)
2131 CALL rotor3( phi,ppi,ppi)
2132 CALL rotor2(thet,hv,hv)
2133 CALL rotor3( phi,hv,hv)
2135 44 hhv(i)=-isgn*hv(i)
2138 ELSEIF(mode.EQ. 1)
THEN
2140 IF(nevraw.EQ.0)
RETURN
2141 pargam=swt/float(nevraw+1)
2143 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
2145 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
2153 7003
FORMAT(///1x,15(5h*****)
2154 $ /,
' *', 25x,
'******** DADMKS INITIALISATION ********',9x,1h*
2155 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
2156 $ /,1x,15(5h*****)/)
2157 7010
FORMAT(///1x,15(5h*****)
2158 $ /,
' *', 25x,
'******** DADMKS FINAL REPORT ********',9x,1h*
2159 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF K* DECAYS TOTAL ',9x,1h*,
2160 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF K* DECS. ACCEPTED ',9x,1h*,
2161 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
2162 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (K* DECAY) IN GEV UNITS ',9x,1h*,
2163 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2164 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
2165 $ /,1x,15(5h*****)/)
2166 902
WRITE(iout, 9020)
2167 9020
FORMAT(
' ----- DADMKS: LACK OF INITIALISATION')
2170 SUBROUTINE dphsks(DGAMT,HV,PN,PKS,PKK,PPI,JKST)
2178 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2179 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2181 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2182 * ,ampiz,ampi,amro,gamro,ama1,gama1
2183 * ,amk,amkz,amkst,gamkst
2185 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2186 * ,ampiz,ampi,amro,gamro,ama1,gama1
2187 * ,amk,amkz,amkst,gamkst
2190 REAL hv(4),pt(4),pn(4),pks(4),pkk(4),ppi(4),qq(4),rr1(1)
2194 DATA pi /3.141592653589793238462643/
2198 phspac=1./2**11/pi**5
2210 ams2=(amtau-amnuta)**2
2216 alp1=atan((ams1-amkst**2)/amkst/gamkst)
2217 alp2=atan((ams2-amkst**2)/amkst/gamkst)
2218 alp=alp1+rr1(1)*(alp2-alp1)
2219 amx2=amkst**2+amkst*gamkst*tan(alp)
2221 phspac=phspac*((amx2-amkst**2)**2+(amkst*gamkst)**2)
2223 phspac=phspac*(alp2-alp1)
2228 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2229 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2234 pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2236 phspac=phspac*(4*pi)*(2*pks(3)/amtau)
2239 enpi=( amx**2+ampi**2-amkz**2 ) / ( 2*amx )
2240 pppi=sqrt((enpi-ampi)*(enpi+ampi))
2241 phspac=phspac*(4*pi)*(2*pppi/amx)
2243 CALL sphera(pppi,ppi)
2248 pkk(4)=( amx**2+amkz**2-ampi**2 ) / ( 2*amx )
2249 exe=(pks(4)+pks(3))/amx
2251 CALL bostr3(exe,ppi,ppi)
2252 CALL bostr3(exe,pkk,pkk)
2254 30 qq(i)=ppi(i)-pkk(i)
2256 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
2257 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
2259 31 qq(i)=qq(i)-pks(i)*qqpks/pksd
2262 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
2264 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
2265 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
2266 & +(gv**2-ga**2)*amtau*amnuta*qq2
2269 fks=cabs(bwigs(amx2,amkst,gamkst))**2
2271 amplit=(gfermi*scabib)**2*brak*2*fks
2272 dgamt=1/(2.*amtau)*amplit*phspac
2274 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
2277 ELSEIF(jkst.EQ.20)
THEN
2281 ams2=(amtau-amnuta)**2
2289 alp1=atan((ams1-amkst**2)/amkst/gamkst)
2290 alp2=atan((ams2-amkst**2)/amkst/gamkst)
2291 alp=alp1+rr1(1)*(alp2-alp1)
2292 amx2=amkst**2+amkst*gamkst*tan(alp)
2294 phspac=phspac*((amx2-amkst**2)**2+(amkst*gamkst)**2)
2296 phspac=phspac*(alp2-alp1)
2301 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2302 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2306 pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2308 phspac=phspac*(4*pi)*(2*pks(3)/amtau)
2311 enpi=( amx**2+ampiz**2-amk**2 ) / ( 2*amx )
2312 pppi=sqrt((enpi-ampiz)*(enpi+ampiz))
2313 phspac=phspac*(4*pi)*(2*pppi/amx)
2315 CALL sphera(pppi,ppi)
2320 pkk(4)=( amx**2+amk**2-ampiz**2 ) / ( 2*amx )
2321 exe=(pks(4)+pks(3))/amx
2323 CALL bostr3(exe,ppi,ppi)
2324 CALL bostr3(exe,pkk,pkk)
2326 60 qq(i)=pkk(i)-ppi(i)
2328 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
2329 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
2331 61 qq(i)=qq(i)-pks(i)*qqpks/pksd
2334 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
2336 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
2337 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
2338 & +(gv**2-ga**2)*amtau*amnuta*qq2
2341 fks=cabs(bwigs(amx2,amkst,gamkst))**2
2343 amplit=(gfermi*scabib)**2*brak*2*fks
2344 dgamt=1/(2.*amtau)*amplit*phspac
2346 70 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
2354 SUBROUTINE dphnpi(DGAMT,HVX,PNX,PRX,PPIX,JNPI)
2360 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2361 * ,ampiz,ampi,amro,gamro,ama1,gama1
2362 * ,amk,amkz,amkst,gamkst
2364 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2365 * ,ampiz,ampi,amro,gamro,ama1,gama1
2366 * ,amk,amkz,amkst,gamkst
2367 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2368 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2369 parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
2371 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2373 CHARACTER names(nmode)*31
2378 REAL*8 pn(4),pr(4),ppi(4,9),hv(4)
2379 REAL*4 pnx(4),prx(4),ppix(4,9),hvx(4)
2380 REAL*8 pv(5,9),pt(4),ue(3),be(3)
2381 REAL*8 pawt,amx,ams1,ams2,pa,phs,phsmax,pmin,pmax
2385 REAL*8 gam,bep,phi,a,b,c
2387 REAL*4 rrr(9),rrx(2),rn(1),rr2(1)
2389 DATA pi /3.141592653589793238462643/
2390 DATA wetmax /500*1d-15/
2395 $ sqrt(max(0.d0,(a**2-(b+c)**2)*(a**2-(b-c)**2)))/(2.d0*a)
2398 ampik(i,j)=dcdmas(idffin(i,j))
2402 IF ((jnpi.LE.0).OR.jnpi.GT.100)
THEN
2403 WRITE(6,*)
'JNPI OUTSIDE RANGE DEFINED BY WETMAX; JNPI=',jnpi
2420 phspac = 1./2.**5 /pi**2
2422 4 ps =ps+ampik(i,jnpi)
2427 ams2=(amtau-amnuta)**2
2431 amx2=ams1+ rr2(1)*(ams2-ams1)
2435 phspac=phspac * (ams2-ams1)
2440 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx2)
2441 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2445 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx2)
2447 phspac=phspac * (4.*pi) * (2.*pr(3)/amtau)
2454 qxn=pr(4)*pn(4)-pr(1)*pn(1)-pr(2)*pn(2)-pr(3)*pn(3)
2459 brak=2*(gv**2+ga**2)*(2*pxq*qxn+amx2*pxn)
2460 & -6*(gv**2-ga**2)*amtau*amnuta*amx2
2466 amplit=ccabib**2*gfermi**2/2.* brak*amx2*sigee(amx2,jn)
2467 dgamt=1./(2.*amtau)*amplit*phspac
2476 pv(5,nd)=ampik(nd,jnpi)
2478 pmax=amw-ps+ampik(nd,jnpi)
2481 pmax=pmax+ampik(il,jnpi)
2482 pmin=pmin+ampik(il+1,jnpi)
2484 220 phsmax=phsmax*pawt(pmax,pmin,ampik(il,jnpi))/pmax
2491 223 ams1=ams1+ampik(jl,jnpi)
2493 amx =(amx-ampik(il,jnpi))
2495 phsmax=phsmax * (ams2-ams1)
2502 phspac = 1./2.**(6*nd-7) /pi**(3*nd-4)
2504 CALL ranmar(rrr,nd-2)
2508 231 ams1=ams1+ampik(jl,jnpi)
2510 ams2=(amx-ampik(il,jnpi))**2
2512 amx2=ams1+ rr1*(ams2-ams1)
2515 phspac=phspac * (ams2-ams1)
2517 phs=phs* (ams2-ams1)
2518 pa=pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2519 phs =phs *pa/pv(5,il)
2521 pa=pawt(pv(5,nd-1),ampik(nd-1,jnpi),ampik(nd,jnpi))
2522 phs =phs *pa/pv(5,nd-1)
2524 wetmax(jnpi)=1.2d0*max(wetmax(jnpi)/1.2d0,phs/phsmax)
2525 IF (ncont.EQ.500 000)
THEN
2528 xnpi=xnpi+ampik(kk,jnpi)
2530 WRITE(6,*)
'ROUNDING INSTABILITY IN DPHNPI ?'
2531 WRITE(6,*)
'AMW=',amw,
'XNPI=',xnpi
2532 WRITE(6,*)
'IF =AMW= IS NEARLY EQUAL =XNPI= THAT IS IT'
2533 WRITE(6,*)
'PHS=',phs,
'PHSMAX=',phsmax
2536 IF(rn(1)*phsmax*wetmax(jnpi).GT.phs) go to 100
2539 280
DO 300 il=1,nd-1
2540 pa=pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2545 ue(1)=sqrt(1.d0-ue(3)**2)*cos(phi)
2546 ue(2)=sqrt(1.d0-ue(3)**2)*sin(phi)
2550 290 pv(j,il+1)=-pa*ue(j)
2551 ppi(4,il)=sqrt(pa**2+ampik(il,jnpi)**2)
2552 pv(4,il+1)=sqrt(pa**2+pv(5,il+1)**2)
2553 phspac=phspac *(4.*pi)*(2.*pa/pv(5,il))
2557 310 ppi(j,nd)=pv(j,nd)
2560 320 be(j)=pv(j,il)/pv(4,il)
2561 gam=pv(4,il)/pv(5,il)
2563 bep=be(1)*ppi(1,i)+be(2)*ppi(2,i)+be(3)*ppi(3,i)
2566 330 ppi(j,i)=ppi(j,i)+gam*(gam*bep/(1.d0+gam)+ppi(4,i))*be(j)
2568 ppi(4,i)=gam*(ppi(4,i)+bep)
2587 FUNCTION sigee(Q2,JNP)
2606 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2607 * ,ampiz,ampi,amro,gamro,ama1,gama1
2608 * ,amk,amkz,amkst,gamkst
2610 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2611 * ,ampiz,ampi,amro,gamro,ama1,gama1
2612 * ,amk,amkz,amkst,gamkst
2616 1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
2617 2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
2618 3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
2619 4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
2622 7 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25,
2625 DATA pi /3.141592653589793238462643/
2630 IF(jnpi.GT.6) jnpi=6
2643 datsig(i,2) = datsig(i,2)/2.
2644 datsig(i,1) = datsig(i,1) + datsig(i,2)
2650 IF(t . gt. s-ampi ) go to 201
2652 fact=(t2/s2)**2*sqrt((s2-t2-ampi2)**2-4.*t2*ampi2)/s2 *2.*t*.05
2653 fact = fact * (datsig(j,1)+datsig(j+1,1))
2654 200 datsig(i,3) = datsig(i,3) + fact
2655 201 datsig(i,3) = datsig(i,3) /(2*pi*fpi)**2
2656 datsig(i,4) = datsig(i,3)
2657 datsig(i,6) = datsig(i,5)
2660 1000
FORMAT(///1x,
' EE SIGMA USED IN MULTIPI DECAYS'/
2666 sigee=datsig(1,jnpi)+
2667 & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
2668 ELSEIF(q.LT.1.8)
THEN
2671 IF(q.LT.qmax) go to 2
2674 2 sigee=datsig(i,jnpi)+
2675 & (datsig(i+1,jnpi)-datsig(i,jnpi)) * (q-qmin)/.05
2676 ELSEIF(q.GT.1.8)
THEN
2677 sigee=datsig(17,jnpi)+
2678 & (datsig(17,jnpi)-datsig(16,jnpi)) * (q-1.8)/.05
2680 IF(sigee.LT..0) sigee=0.
2682 sigee = sigee/(6.*pi**2*sig0)
2687 FUNCTION sigold(Q2,JNPI)
2706 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2707 * ,ampiz,ampi,amro,gamro,ama1,gama1
2708 * ,amk,amkz,amkst,gamkst
2710 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2711 * ,ampiz,ampi,amro,gamro,ama1,gama1
2712 * ,amk,amkz,amkst,gamkst
2716 1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
2717 2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
2718 3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
2719 4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
2721 6 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25/
2723 DATA pi /3.141592653589793238462643/
2731 datsig(i,2) = datsig(i,2)/2.
2732 datsig(i,1) = datsig(i,1) + datsig(i,2)
2738 IF(t . gt. s-ampi ) go to 201
2740 fact=(t2/s2)**2*sqrt((s2-t2-ampi2)**2-4.*t2*ampi2)/s2 *2.*t*.05
2741 fact = fact * (datsig(j,1)+datsig(j+1,1))
2742 200 datsig(i,3) = datsig(i,3) + fact
2743 201 datsig(i,3) = datsig(i,3) /(2*pi*fpi)**2
2746 1000
FORMAT(///1x,
' EE SIGMA USED IN MULTIPI DECAYS'/
2752 sigee=datsig(1,jnpi)+
2753 & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
2754 ELSEIF(q.LT.1.8)
THEN
2757 IF(q.LT.qmax) go to 2
2760 2 sigee=datsig(i,jnpi)+
2761 & (datsig(i+1,jnpi)-datsig(i,jnpi)) * (q-qmin)/.05
2762 ELSEIF(q.GT.1.8)
THEN
2763 sigee=datsig(17,jnpi)+
2764 & (datsig(17,jnpi)-datsig(16,jnpi)) * (q-1.8)/.05
2766 IF(sigee.LT..0) sigee=0.
2768 sigee = sigee/(6.*pi**2*sig0)
2773 SUBROUTINE dphspk(DGAMT,HV,PN,PAA,PNPI,JAA)
2778 parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
2780 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2783 CHARACTER names(nmode)*31
2785 REAL hv(4),pn(4),paa(4),pim1(4),pim2(4),pipl(4),pnpi(4,9)
2793 amp1=dcdmas(idffin(1,jaa+nm4+nm5+nm6))
2794 amp2=dcdmas(idffin(2,jaa+nm4+nm5+nm6))
2795 amp3=dcdmas(idffin(3,jaa+nm4+nm5+nm6))
2797 $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
2809 $ dphtre(dgamt,hv,pn,paa,pim1,ampa,pim2,ampb,pipl,amp3,keyt,mnum)
2826 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2827 * ,ampiz,ampi,amro,gamro,ama1,gama1
2828 * ,amk,amkz,amkst,gamkst
2830 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
2831 * ,ampiz,ampi,amro,gamro,ama1,gama1
2832 * ,amk,amkz,amkst,gamkst
2833 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2834 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
2835 REAL hv(4),pt(4),pn(4),paa(4),pim1(4),pim2(4),pipl(4)
2838 DATA pi /3.141592653589793238462643/
2840 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
2845 phspac=1./2**17/pi**8
2855 CALL choice(mnum,rr,ichan,prob1,prob2,prob3,
2856 $ amrx,gamrx,amra,gamra,amrb,gamrb)
2857 IF (ichan.EQ.1)
THEN
2860 ELSEIF (ichan.EQ.2)
THEN
2869 ams1=(amp1+amp2+amp3)**2
2870 ams2=(amtau-amnuta)**2
2874 alp1=atan((ams1-amrx**2)/amrx/gamrx)
2875 alp2=atan((ams2-amrx**2)/amrx/gamrx)
2876 alp=alp1+rr1*(alp2-alp1)
2877 am3sq =amrx**2+amrx*gamrx*tan(alp)
2879 phspac=phspac*((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
2880 phspac=phspac*(alp2-alp1)
2885 IF (ichan.LE.2)
THEN
2889 alp1=atan((ams1-amra**2)/amra/gamra)
2890 alp2=atan((ams2-amra**2)/amra/gamra)
2891 alp=alp1+rr2*(alp2-alp1)
2892 am2sq =amra**2+amra*gamra*tan(alp)
2902 am2sq=ams1+ rr2*(ams2-ams1)
2909 enq1=(am2sq-amp2**2+amp3**2)/(2*am2)
2910 enq2=(am2sq+amp2**2-amp3**2)/(2*am2)
2911 ppi= enq1**2-amp3**2
2912 pppi=sqrt(abs(enq1**2-amp3**2))
2914 phf1=(4*pi)*(2*pppi/am2)
2918 CALL sphera(pppi,pipl)
2932 pr(4)=1./(2*am3)*(am3**2+am2**2-amp1**2)
2933 pr(3)= sqrt(abs(pr(4)**2-am2**2))
2934 ppi = pr(4)**2-am2**2
2938 pim2(4)=1./(2*am3)*(am3**2-am2**2+amp1**2)
2940 phf2=(4*pi)*(2*pr(3)/am3)
2944 exe=(pr(4)+pr(3))/am2
2945 CALL bostr3(exe,pipl,pipl)
2946 CALL bostr3(exe,pim1,pim1)
2952 thet =acos(-1.+2*rr3)
2954 CALL rotpol(thet,phi,pipl)
2955 CALL rotpol(thet,phi,pim1)
2956 CALL rotpol(thet,phi,pim2)
2957 CALL rotpol(thet,phi,pr)
2963 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am3**2)
2964 paa(3)= sqrt(abs(paa(4)**2-am3**2))
2965 ppi = paa(4)**2-am3**2
2966 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
2970 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am3**2)
2976 alp1=atan((ams1-amra**2)/amra/gamra)
2977 alp2=atan((ams2-amra**2)/amra/gamra)
2978 xpro = (pim1(3)+pipl(3))**2
2979 $ +(pim1(2)+pipl(2))**2+(pim1(1)+pipl(1))**2
2980 am2sq=-xpro+(pim1(4)+pipl(4))**2
2982 ff1 = ((am2sq-amra**2)**2+(amra*gamra)**2)/(amra*gamra)
2983 ff1 =ff1 *(alp2-alp1)
2985 gg1 = (4*pi)*(xlam(am2sq,amp2**2,amp3**2)/am2sq)
2987 gg1 =gg1 *(4*pi)*sqrt(4*xpro/am3sq)
2988 xjaje=gg1*(ams2-ams1)
2992 alp1=atan((ams1-amrb**2)/amrb/gamrb)
2993 alp2=atan((ams2-amrb**2)/amrb/gamrb)
2994 xpro = (pim2(3)+pipl(3))**2
2995 $ +(pim2(2)+pipl(2))**2+(pim2(1)+pipl(1))**2
2996 am2sq=-xpro+(pim2(4)+pipl(4))**2
2997 ff2 = ((am2sq-amrb**2)**2+(amrb*gamrb)**2)/(amrb*gamrb)
2998 ff2 =ff2 *(alp2-alp1)
2999 gg2 = (4*pi)*(xlam(am2sq,amp1**2,amp3**2)/am2sq)
3000 gg2 =gg2 *(4*pi)*sqrt(4*xpro/am3sq)
3001 xjadw=gg2*(ams2-ams1)
3008 IF (ichan.EQ.2)
THEN
3013 IF (xjac1.NE.0.0) a1=prob1/xjac1
3014 IF (xjac2.NE.0.0) a2=prob2/xjac2
3015 IF (xjac3.NE.0.0) a3=prob3/xjac3
3017 IF (a1+a2+a3.NE.0.0)
THEN
3018 phspac=phspac/(a1+a2+a3)
3030 exe=(paa(4)+paa(3))/am3
3031 CALL bostr3(exe,pipl,pipl)
3032 CALL bostr3(exe,pim1,pim1)
3033 CALL bostr3(exe,pim2,pim2)
3034 CALL bostr3(exe,pr,pr)
3037 CALL dampog(pt,pn,pim1,pim2,pipl,amplit,hv)
3042 CALL damppk(mnum,pt,pn,pim1,pim2,pipl,amplit,hv)
3047 IF (keyt.EQ.1.OR.keyt.EQ.2)
THEN
3054 dgamt=1/(2.*amtau)*amplit*phspac
3056 SUBROUTINE dampaa(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3066 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3067 * ,ampiz,ampi,amro,gamro,ama1,gama1
3068 * ,amk,amkz,amkst,gamkst
3070 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3071 * ,ampiz,ampi,amro,gamro,ama1,gama1
3072 * ,amk,amkz,amkst,gamkst
3073 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3074 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3075 COMMON /testa1/ keya1
3076 REAL hv(4),pt(4),pn(4),pim1(4),pim2(4),pipl(4)
3077 REAL paa(4),vec1(4),vec2(4)
3078 REAL pivec(4),piaks(4),hvm(4)
3079 COMPLEX bwign,hadcur(4),fpik
3086 bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
3090 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3092 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3093 xmro1 =sqrt(abs((pipl(4)+pim1(4))**2-(pipl(1)+pim1(1))**2
3094 $ -(pipl(2)+pim1(2))**2-(pipl(3)+pim1(3))**2))
3095 xmro2 =sqrt(abs((pipl(4)+pim2(4))**2-(pipl(1)+pim2(1))**2
3096 $ -(pipl(2)+pim2(2))**2-(pipl(3)+pim2(3))**2))
3098 prod1 =paa(4)*(pim1(4)-pipl(4))-paa(1)*(pim1(1)-pipl(1))
3099 $ -paa(2)*(pim1(2)-pipl(2))-paa(3)*(pim1(3)-pipl(3))
3100 prod2 =paa(4)*(pim2(4)-pipl(4))-paa(1)*(pim2(1)-pipl(1))
3101 $ -paa(2)*(pim2(2)-pipl(2))-paa(3)*(pim2(3)-pipl(3))
3103 vec1(i)= pim1(i)-pipl(i) -paa(i)*prod1/xmaa**2
3104 40 vec2(i)= pim2(i)-pipl(i) -paa(i)*prod2/xmaa**2
3106 IF (keya1.EQ.1)
THEN
3110 fnorm=fa1/sqrt(2.)*faropi*fro2pi
3112 hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gama1)
3113 $ *(cmplx(vec1(i))*amro**2*bwign(xmro1,amro,gamro)
3114 $ +cmplx(vec2(i))*amro**2*bwign(xmro2,amro,gamro))
3117 fnorm=2.0*sqrt(2.)/3.0/fpi
3118 gamax=gama1*gfun(xmaa**2)/gfun(ama1**2)
3120 hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gamax)
3121 $ *(cmplx(vec1(i))*fpik(xmro1)
3122 $ +cmplx(vec2(i))*fpik(xmro2))
3127 CALL clvec(hadcur,pn,pivec)
3128 CALL claxi(hadcur,pn,piaks)
3129 CALL clnut(hadcur,brakm,hvm)
3131 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3132 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3133 amplit=(gfermi*ccabib)**2*brak/2.
3138 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3139 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3149 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3150 * ,ampiz,ampi,amro,gamro,ama1,gama1
3151 * ,amk,amkz,amkst,gamkst
3153 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3154 * ,ampiz,ampi,amro,gamro,ama1,gama1
3155 * ,amk,amkz,amkst,gamkst
3157 IF (qkwa.LT.(amro+ampi)**2)
THEN
3158 gfun=4.1*(qkwa-9*ampiz**2)**3
3159 $ *(1.-3.3*(qkwa-9*ampiz**2)+5.8*(qkwa-9*ampiz**2)**2)
3161 gfun=qkwa*(1.623+10.38/qkwa-9.32/qkwa**2+0.65/qkwa**3)
3164 COMPLEX FUNCTION bwigs(S,M,G)
3169 REAL pi,pim,qs,qm,w,gs,mk
3176 p(a,b,c)=sqrt(abs(abs(((a+b-c)**2-4.*a*b)/4./a)
3177 $ +(((a+b-c)**2-4.*a*b)/4./a))/2.0)
3193 qs=p(s,pim**2,mk**2)
3194 qm=p(m**2,pim**2,mk**2)
3196 gs=g*(m/w)*(qs/qm)**3
3198 bw=m**2/cmplx(m**2-s,-m*gs)
3199 qpm=p(pm**2,pim**2,mk**2)
3200 g1=pg*(pm/w)*(qs/qpm)**3
3201 bwp=pm**2/cmplx(pm**2-s,-pm*g1)
3202 bwigs= (bw+pbeta*bwp)/(1+pbeta)
3206 COMPLEX FUNCTION bwig(S,M,G)
3211 REAL pi,pim,qs,qm,w,gs
3220 IF (s.GT.4.*pim**2)
THEN
3221 qs=sqrt(abs(abs(s/4.-pim**2)+(s/4.-pim**2))/2.0)
3222 qm=sqrt(m**2/4.-pim**2)
3224 gs=g*(m/w)*(qs/qm)**3
3228 bwig=m**2/cmplx(m**2-s,-m*gs)
3231 COMPLEX FUNCTION fpik(W)
3236 REAL rom,rog,rom1,rog1,beta1,pi,pim,s,w
3241 IF (init.EQ.0 )
THEN
3255 fpik= (bwig(s,rom,rog)+beta1*bwig(s,rom1,rog1))
3264 fpirho=cabs(fpik(w))**2
3266 SUBROUTINE clvec(HJ,PN,PIV)
3276 hn= hj(4)*cmplx(pn(4))-hj(3)*cmplx(pn(3))
3277 hh=
REAL(hj(4)*conjg(hj(4))-hj(3)*conjg(hj(3))
3278 $ -hj(2)*conjg(hj(2))-hj(1)*conjg(hj(1)))
3280 10 piv(i)=4.*
REAL(hn*conjg(hj(i)))-2.*hh*pn(i)
3283 SUBROUTINE claxi(HJ,PN,PIA)
3290 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3291 COMMON / idfc / idff
3293 COMPLEX hj(4),hjc(4)
3296 det2(i,j)=aimag(hjc(i)*hj(j)-hjc(j)*hj(i))
3301 IF (ktom.EQ.1.OR.ktom.EQ.-1)
THEN
3302 sign= idff/abs(idff)
3303 ELSEIF (ktom.EQ.2)
THEN
3304 sign=-idff/abs(idff)
3306 print *,
'STOP IN CLAXI: KTOM=',ktom
3311 10 hjc(i)=conjg(hj(i))
3312 pia(1)= -2.*pn(3)*det2(2,4)+2.*pn(4)*det2(2,3)
3313 pia(2)= -2.*pn(4)*det2(1,3)+2.*pn(3)*det2(1,4)
3314 pia(3)= 2.*pn(4)*det2(1,2)
3315 pia(4)= 2.*pn(3)*det2(1,2)
3318 20 pia(i)=pia(i)*sign
3320 SUBROUTINE clnut(HJ,B,HV)
3332 b=
REAL( HJ(4)*AIMAG(HJ(4)) - HJ(3)*AIMAG(HJ(3))
& - HJ(2)*AIMAG(HJ(2)) - HJ(1)*AIMAG(HJ(1)) )
3335 SUBROUTINE dampog(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3347 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3348 * ,ampiz,ampi,amro,gamro,ama1,gama1
3349 * ,amk,amkz,amkst,gamkst
3351 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3352 * ,ampiz,ampi,amro,gamro,ama1,gama1
3353 * ,amk,amkz,amkst,gamkst
3354 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3355 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3356 COMMON /testa1/ keya1
3357 REAL hv(4),pt(4),pn(4),pim1(4),pim2(4),pipl(4)
3358 REAL paa(4),vec1(4),vec2(4)
3359 REAL pivec(4),piaks(4),hvm(4)
3360 COMPLEX bwign,hadcur(4),fnorm,formom
3372 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3375 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3376 xmom =sqrt(abs( (pim2(4)+pipl(4))**2-(pim2(3)+pipl(3))**2
3377 $ -(pim2(2)+pipl(2))**2-(pim2(1)+pipl(1))**2 ))
3378 xmro2 =(pipl(1))**2 +(pipl(2))**2 +(pipl(3))**2
3380 prod1 =vec1(1)*pipl(1)
3381 prod2 =vec2(2)*pipl(2)
3382 p12 =pim1(4)*pim2(4)-pim1(1)*pim2(1)
3383 $ -pim1(2)*pim2(2)-pim1(3)*pim2(3)
3384 p1pl =pim1(4)*pipl(4)-pim1(1)*pipl(1)
3385 $ -pim1(2)*pipl(2)-pim1(3)*pipl(3)
3386 p2pl =pipl(4)*pim2(4)-pipl(1)*pim2(1)
3387 $ -pipl(2)*pim2(2)-pipl(3)*pim2(3)
3389 vec1(i)= (vec1(i)-prod1/xmro2*pipl(i))
3391 gnorm=sqrt(vec1(1)**2+vec1(2)**2+vec1(3)**2)
3393 vec1(i)= vec1(i)/gnorm
3395 vec2(1)=(vec1(2)*pipl(3)-vec1(3)*pipl(2))/sqrt(xmro2)
3396 vec2(2)=(vec1(3)*pipl(1)-vec1(1)*pipl(3))/sqrt(xmro2)
3397 vec2(3)=(vec1(1)*pipl(2)-vec1(2)*pipl(1))/sqrt(xmro2)
3398 p1vec1 =pim1(4)*vec1(4)-pim1(1)*vec1(1)
3399 $ -pim1(2)*vec1(2)-pim1(3)*vec1(3)
3400 p2vec1 =vec1(4)*pim2(4)-vec1(1)*pim2(1)
3401 $ -vec1(2)*pim2(2)-vec1(3)*pim2(3)
3402 p1vec2 =pim1(4)*vec2(4)-pim1(1)*vec2(1)
3403 $ -pim1(2)*vec2(2)-pim1(3)*vec2(3)
3404 p2vec2 =vec2(4)*pim2(4)-vec2(1)*pim2(1)
3405 $ -vec2(2)*pim2(2)-vec2(3)*pim2(3)
3407 fnorm=formom(xmaa,xmom)
3412 hadcur(i) = fnorm *(
3413 $ vec1(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3414 $ -pim2(i)*(p2vec1*p1pl-p1vec1*p2pl)
3415 $ +pipl(i)*(p2vec1*p12 -p1vec1*(ampi**2+p2pl)) )
3417 hadcur(i) = fnorm *(
3418 $ vec2(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3419 $ -pim2(i)*(p2vec2*p1pl-p1vec2*p2pl)
3420 $ +pipl(i)*(p2vec2*p12 -p1vec2*(ampi**2+p2pl)) )
3425 CALL clvec(hadcur,pn,pivec)
3426 CALL claxi(hadcur,pn,piaks)
3427 CALL clnut(hadcur,brakm,hvm)
3429 brak=brak+(gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3430 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3432 hv(i)=hv(i)-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3433 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3437 amplit=(gfermi*ccabib)**2*brak/2.
3446 SUBROUTINE damppk(MNUM,PT,PN,PIM1,PIM2,PIM3,AMPLIT,HV)
3458 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3459 * ,ampiz,ampi,amro,gamro,ama1,gama1
3460 * ,amk,amkz,amkst,gamkst
3462 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3463 * ,ampiz,ampi,amro,gamro,ama1,gama1
3464 * ,amk,amkz,amkst,gamkst
3465 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3466 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3467 REAL hv(4),pt(4),pn(4),pim1(4),pim2(4),pim3(4)
3468 REAL paa(4),vec1(4),vec2(4),vec3(4),vec4(4),vec5(4)
3469 REAL pivec(4),piaks(4),hvm(4)
3470 REAL fnorm(0:19),coef(1:5,0:19)
3471 COMPLEX hadcur(4),form1,form2,form3,form4,form5,uroj
3473 COMPLEX f1,f2,f3,f4,f5
3475 EXTERNAL form1,form2,form3,form4,form5
3476 DATA pi /3.141592653589793238462643/
3480 IF (icont.EQ.0)
THEN
3488 fnorm(4)=scabib/fpi/dwapi0
3498 coef(1,0)= 2.0*sqrt(2.)/3.0
3499 coef(2,0)=-2.0*sqrt(2.)/3.0
3502 coef(3,0)= 2.0*sqrt(2.)/3.0
3507 coef(1,1)=-sqrt(2.)/3.0
3508 coef(2,1)= sqrt(2.)/3.0
3513 coef(1,2)=-sqrt(2.)/3.0
3514 coef(2,2)= sqrt(2.)/3.0
3528 coef(1,4)= 1.0/sqrt(2.)/3.0
3529 coef(2,4)=-1.0/sqrt(2.)/3.0
3534 coef(1,5)=-sqrt(2.)/3.0
3535 coef(2,5)= sqrt(2.)/3.0
3553 coef(5,7)=-sqrt(2.0/3.0)
3562 coef(1,9)= 2.0*sqrt(2.)/3.0
3563 coef(2,9)=-2.0*sqrt(2.)/3.0
3564 coef(3,9)= 2.0*sqrt(2.)/3.0
3570 coef(1,k)= 2.0*sqrt(2.)/3.0
3571 coef(2,k)=-2.0*sqrt(2.)/3.0
3572 coef(3,k)= 2.0*sqrt(2.)/3.0
3583 10 paa(i)=pim1(i)+pim2(i)+pim3(i)
3584 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3585 xmro1 =sqrt(abs((pim3(4)+pim2(4))**2-(pim3(1)+pim2(1))**2
3586 $ -(pim3(2)+pim2(2))**2-(pim3(3)+pim2(3))**2))
3587 xmro2 =sqrt(abs((pim3(4)+pim1(4))**2-(pim3(1)+pim1(1))**2
3588 $ -(pim3(2)+pim1(2))**2-(pim3(3)+pim1(3))**2))
3589 xmro3 =sqrt(abs((pim1(4)+pim2(4))**2-(pim1(1)+pim2(1))**2
3590 $ -(pim1(2)+pim2(2))**2-(pim1(3)+pim2(3))**2))
3592 prod1 =paa(4)*(pim2(4)-pim3(4))-paa(1)*(pim2(1)-pim3(1))
3593 $ -paa(2)*(pim2(2)-pim3(2))-paa(3)*(pim2(3)-pim3(3))
3594 prod2 =paa(4)*(pim3(4)-pim1(4))-paa(1)*(pim3(1)-pim1(1))
3595 $ -paa(2)*(pim3(2)-pim1(2))-paa(3)*(pim3(3)-pim1(3))
3596 prod3 =paa(4)*(pim1(4)-pim2(4))-paa(1)*(pim1(1)-pim2(1))
3597 $ -paa(2)*(pim1(2)-pim2(2))-paa(3)*(pim1(3)-pim2(3))
3599 vec1(i)= pim2(i)-pim3(i) -paa(i)*prod1/xmaa**2
3600 vec2(i)= pim3(i)-pim1(i) -paa(i)*prod2/xmaa**2
3601 vec3(i)= pim1(i)-pim2(i) -paa(i)*prod3/xmaa**2
3602 40 vec4(i)= pim1(i)+pim2(i)+pim3(i)
3603 CALL prod5(pim1,pim2,pim3,vec5)
3608 f1 = cmplx(coef(1,mnum))*form1(mnum,xmaa**2,xmro1**2,xmro2**2)
3609 f2 = cmplx(coef(2,mnum))*form2(mnum,xmaa**2,xmro2**2,xmro1**2)
3610 f3 = cmplx(coef(3,mnum))*form3(mnum,xmaa**2,xmro3**2,xmro1**2)
3612 $cmplx(coef(4,mnum))*form4(mnum,xmaa**2,xmro1**2,xmro2**2,xmro3**2)
3613 f5 = (-1.0)*uroj/4.0/pi**2/fpi**2*
3614 $ cmplx(coef(5,mnum))*form5(mnum,xmaa**2,xmro1**2,xmro2**2)
3617 hadcur(i)= cmplx(fnorm(mnum)) * (
3618 $ cmplx(vec1(i))*f1+cmplx(vec2(i))*f2+cmplx(vec3(i))*f3+
3619 $ cmplx(vec4(i))*f4+cmplx(vec5(i))*f5)
3624 CALL clvec(hadcur,pn,pivec)
3625 CALL claxi(hadcur,pn,piaks)
3626 CALL clnut(hadcur,brakm,hvm)
3628 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3629 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3630 amplit=(gfermi)**2*brak/2.
3631 IF (mnum.GE.20)
THEN
3632 print *,
'MNUM=',mnum
3638 IF (k.EQ.4) znak=1.0
3639 xm1=znak*pim1(k)**2+xm1
3640 xm2=znak*pim2(k)**2+xm2
3641 xm3=znak*pim3(k)**2+xm3
3642 77 print *,
'PIM1=',pim1(k),
'PIM2=',pim2(k),
'PIM3=',pim3(k)
3643 print *,
'XM1=',sqrt(xm1),
'XM2=',sqrt(xm2),
'XM3=',sqrt(xm3)
3644 print *,
'************************************************'
3648 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3649 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3654 SUBROUTINE prod5(P1,P2,P3,PIA)
3660 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3661 COMMON / idfc / idff
3662 REAL pia(4),p1(4),p2(4),p3(4)
3663 det2(i,j)=p1(i)*p2(j)-p2(i)*p1(j)
3665 IF (ktom.EQ.1.OR.ktom.EQ.-1)
THEN
3666 sign= idff/abs(idff)
3667 ELSEIF (ktom.EQ.2)
THEN
3668 sign=-idff/abs(idff)
3670 print *,
'STOP IN PROD5: KTOM=',ktom
3676 pia(1)= -p3(3)*det2(2,4)+p3(4)*det2(2,3)+p3(2)*det2(3,4)
3677 pia(2)= -p3(4)*det2(1,3)+p3(3)*det2(1,4)-p3(1)*det2(3,4)
3678 pia(3)= p3(4)*det2(1,2)-p3(2)*det2(1,4)+p3(1)*det2(2,4)
3679 pia(4)= p3(3)*det2(1,2)-p3(2)*det2(1,3)+p3(1)*det2(2,3)
3682 20 pia(i)=pia(i)*sign
3685 SUBROUTINE dexnew(MODE,ISGN,POL,PNU,PAA,PNPI,JNPI)
3698 COMMON / inout / inut,iout
3699 REAL pol(4),hv(4),paa(4),pnu(4),pnpi(4,9),rn(1)
3705 CALL dadnew( -1,isgn,hv,pnu,paa,pnpi,jdumm)
3710 ELSEIF(mode.EQ. 0)
THEN
3713 IF(iwarm.EQ.0) goto 902
3714 CALL dadnew( 0,isgn,hv,pnu,paa,pnpi,jnpi)
3715 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
3718 IF(rn(1).GT.wt) goto 300
3720 ELSEIF(mode.EQ. 1)
THEN
3722 CALL dadnew( 1,isgn,hv,pnu,paa,pnpi,jdumm)
3727 902
WRITE(iout, 9020)
3728 9020
FORMAT(
' ----- DEXNEW: LACK OF INITIALISATION')
3731 SUBROUTINE dadnew(MODE,ISGN,HV,PNU,PWB,PNPI,JNPI)
3733 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3734 * ,ampiz,ampi,amro,gamro,ama1,gama1
3735 * ,amk,amkz,amkst,gamkst
3737 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3738 * ,ampiz,ampi,amro,gamro,ama1,gama1
3739 * ,amk,amkz,amkst,gamkst
3740 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3741 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3742 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
3743 REAL*4 gampmc ,gamper
3745 COMMON / inout / inut,iout
3747 parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
3749 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3752 CHARACTER names(nmode)*31
3755 REAL*4 pnu(4),pwb(4),pnpi(4,9),hv(4),hhv(4)
3756 REAL*4 pdum1(4),pdum2(4),pdumi(4,9)
3759 REAL*8 swt(nmode),sswt(nmode)
3760 dimension nevraw(nmode),nevovr(nmode),nevacc(nmode)
3762 DATA pi /3.141592653589793238462643/
3783 a = pkorb(3,37+jnpi)
3796 ELSEIF(jnpi.LE.nm4)
THEN
3797 CALL dph4pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3799 ELSEIF(jnpi.LE.nm4+nm5)
THEN
3800 CALL dph5pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3801 ELSEIF(jnpi.LE.nm4+nm5+nm6)
THEN
3802 CALL dphnpi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3803 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3)
THEN
3804 inum=jnpi-nm4-nm5-nm6
3805 CALL dphspk(wt,hv,pdum1,pdum2,pdumi,inum)
3806 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2)
THEN
3807 inum=jnpi-nm4-nm5-nm6-nm3
3808 CALL dphsrk(wt,hv,pdum1,pdum2,pdumi,inum)
3812 IF(wt.GT.wtmax(jnpi)/1.2) wtmax(jnpi)=wt*1.2
3822 ELSEIF(mode.EQ. 0)
THEN
3824 IF(iwarm.EQ.0) goto 902
3829 ELSEIF(jnpi.LE.nm4)
THEN
3830 CALL dph4pi(wt,hhv,pnu,pwb,pnpi,jnpi)
3831 ELSEIF(jnpi.LE.nm4+nm5)
THEN
3832 CALL dph5pi(wt,hhv,pnu,pwb,pnpi,jnpi)
3833 ELSEIF(jnpi.LE.nm4+nm5+nm6)
THEN
3834 CALL dphnpi(wt,hhv,pnu,pwb,pnpi,jnpi)
3835 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3)
THEN
3836 inum=jnpi-nm4-nm5-nm6
3837 CALL dphspk(wt,hhv,pnu,pwb,pnpi,inum)
3838 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2)
THEN
3839 inum=jnpi-nm4-nm5-nm6-nm3
3840 CALL dphsrk(wt,hhv,pnu,pwb,pnpi,inum)
3848 nevraw(jnpi)=nevraw(jnpi)+1
3849 swt(jnpi)=swt(jnpi)+wt
3853 sswt(jnpi)=sswt(jnpi)+dble(wt)**2
3858 IF(wt.GT.wtmax(jnpi)) nevovr(jnpi)=nevovr(jnpi)+1
3859 IF(rn*wtmax(jnpi).GT.wt) goto 300
3861 costhe=-1.+2.*rrr(2)
3864 CALL rotor2(thet,pnu,pnu)
3865 CALL rotor3( phi,pnu,pnu)
3866 CALL rotor2(thet,pwb,pwb)
3867 CALL rotor3( phi,pwb,pwb)
3868 CALL rotor2(thet,hv,hv)
3869 CALL rotor3( phi,hv,hv)
3872 CALL rotor2(thet,pnpi(1,i),pnpi(1,i))
3873 CALL rotor3( phi,pnpi(1,i),pnpi(1,i))
3875 nevacc(jnpi)=nevacc(jnpi)+1
3877 ELSEIF(mode.EQ. 1)
THEN
3880 IF(nevraw(jnpi).EQ.0) goto 500
3881 pargam=swt(jnpi)/float(nevraw(jnpi)+1)
3883 IF(nevraw(jnpi).NE.0)
3884 & error=sqrt(sswt(jnpi)/swt(jnpi)**2-1./float(nevraw(jnpi)))
3886 WRITE(iout, 7010) names(jnpi),
3887 & nevraw(jnpi),nevacc(jnpi),nevovr(jnpi),pargam,rat,error
3889 gampmc(8+jnpi-1)=rat
3890 gamper(8+jnpi-1)=error
3896 7003
FORMAT(///1x,15(5h*****)
3897 $ /,
' *', 25x,
'******** DADNEW INITIALISATION ********',9x,1h*
3899 7004
FORMAT(
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*/)
3901 $ /,1x,15(5h*****)/)
3902 7010
FORMAT(///1x,15(5h*****)
3903 $ /,
' *', 25x,
'******** DADNEW FINAL REPORT ******** ',9x,1h*
3904 $ /,
' *', 25x,
'CHANNEL:',a31 ,9x,1h*
3905 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF DECAYS TOTAL ',9x,1h*
3906 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF DECAYS ACCEPTED ',9x,1h*
3907 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
3908 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH IN GEV UNITS ',9x,1h*
3909 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
3910 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
3911 $ /,1x,15(5h*****)/)
3912 902
WRITE(iout, 9020)
3913 9020
FORMAT(
' ----- DADNEW: LACK OF INITIALISATION')
3915 903
WRITE(iout, 9030) jnpi,mode
3916 9030
FORMAT(
' ----- DADNEW: WRONG JNPI',2i5)
3921 SUBROUTINE dph4pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
3928 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3929 * ,ampiz,ampi,amro,gamro,ama1,gama1
3930 * ,amk,amkz,amkst,gamkst
3932 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
3933 * ,ampiz,ampi,amro,gamro,ama1,gama1
3934 * ,amk,amkz,amkst,gamkst
3935 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3936 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
3939 REAL hv(4),pt(4),pn(4),paa(4),pim1(4),pim2(4),pipl(4),pmult(4,9)
3942 REAL*8 uu,ff,ff1,ff2,ff3,ff4,gg1,gg2,gg3,gg4,rr
3943 DATA pi /3.141592653589793238462643/
3945 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
3950 phspac=1./2**23/pi**11
3985 CALL choice(100+jnpi,rrb,ichan,prob1,prob2,prob3,
3987 $ amrop,gamrop,amrx,gamrx,amrb,gamrb)
4001 ams1=(amp1+amp2+amp3+amp4)**2
4002 ams2=(amtau-amnuta)**2
4003 alp1=atan((ams1-amrop**2)/amrop/gamrop)
4004 alp2=atan((ams2-amrop**2)/amrop/gamrop)
4005 alp=alp1+rr1*(alp2-alp1)
4006 am4sq =amrop**2+amrop*gamrop*tan(alp)
4009 $ ((am4sq-amrop**2)**2+(amrop*gamrop)**2)/(amrop*gamrop)
4010 phspac=phspac*(alp2-alp1)
4014 ams1=(amp2+amp3+amp4)**2
4016 IF (rrr(9).GT.prez)
THEN
4017 am3sq=ams1+ rr1*(ams2-ams1)
4023 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4024 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4025 alp=alp1+rr1*(alp2-alp1)
4026 am3sq =amrx**2+amrx*gamrx*tan(alp)
4029 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4037 am2sq=ams1+ rr2*(ams2-ams1)
4043 enq1=(am2sq-amp3**2+amp4**2)/(2*am2)
4044 enq2=(am2sq+amp3**2-amp4**2)/(2*am2)
4045 ppi= enq1**2-amp4**2
4046 pppi=sqrt(abs(enq1**2-amp4**2))
4048 phspac=phspac*(4*pi)*(2*pppi/am2)
4052 CALL sphera(pppi,piz)
4066 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4067 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4068 ppi = pr(4)**2-am2**2
4074 pim1(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4077 ff3=(4*pi)*(2*pr(3)/am3)
4079 exe=(pr(4)+pr(3))/am2
4080 CALL bostr3(exe,piz,piz)
4081 CALL bostr3(exe,pipl,pipl)
4084 thet =acos(-1.+2*rr3)
4086 CALL rotpol(thet,phi,pipl)
4087 CALL rotpol(thet,phi,pim1)
4088 CALL rotpol(thet,phi,piz)
4089 CALL rotpol(thet,phi,pr)
4096 pr(4)=1./(2*am4)*(am4**2+am3**2-amp1**2)
4097 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4098 ppi = pr(4)**2-am3**2
4104 pim2(4)=1./(2*am4)*(am4**2-am3**2+amp1**2)
4107 ff4=(4*pi)*(2*pr(3)/am4)
4109 exe=(pr(4)+pr(3))/am3
4110 CALL bostr3(exe,piz,piz)
4111 CALL bostr3(exe,pipl,pipl)
4112 CALL bostr3(exe,pim1,pim1)
4115 thet =acos(-1.+2*rr3)
4117 CALL rotpol(thet,phi,pipl)
4118 CALL rotpol(thet,phi,pim1)
4119 CALL rotpol(thet,phi,pim2)
4120 CALL rotpol(thet,phi,piz)
4121 CALL rotpol(thet,phi,pr)
4127 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am4**2)
4128 paa(3)= sqrt(abs(paa(4)**2-am4**2))
4129 ppi = paa(4)**2-am4**2
4130 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4131 phsp=phsp*(4*pi)*(2*paa(3)/amtau)
4135 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am4**2)
4138 IF(rrr(9).LE.0.5*prez)
THEN
4147 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4148 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4150 ams1=(amp1+amp3+amp4)**2
4153 ams2=(sqrt(am3sq)-amp1)**2
4155 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4156 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4159 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4160 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4162 ams1=(amp1+amp3+amp4)**2
4163 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4164 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4165 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4168 ams2=(sqrt(am3sq)-amp1)**2
4170 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4171 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4174 am3sq=(pim2(4)+piz(4)+pipl(4))**2-(pim2(3)+piz(3)+pipl(3))**2
4175 $ -(pim2(2)+piz(2)+pipl(2))**2-(pim2(1)+piz(1)+pipl(1))**2
4177 ams1=(amp2+amp3+amp4)**2
4178 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4179 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4180 gg1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4183 ams2=(sqrt(am3sq)-amp2)**2
4185 gg3=(4*pi)*(xlam(am2**2,amp2**2,am3sq)/am3sq)
4186 gg4=(4*pi)*(xlam(am3sq,amp1**2,am4**2)/am4**2)
4189 IF ( ( (ff+gg)*uu+ff*gg ).GT.0.0d0)
THEN
4190 rr=ff*gg*uu/(0.5*prez*(ff+gg)*uu+(1.0-prez)*ff*gg)
4218 exe=(paa(4)+paa(3))/am4
4219 CALL bostr3(exe,piz,piz)
4220 CALL bostr3(exe,pipl,pipl)
4221 CALL bostr3(exe,pim1,pim1)
4222 CALL bostr3(exe,pim2,pim2)
4223 CALL bostr3(exe,pr,pr)
4233 CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4234 ELSEIF (jnpi.EQ.2)
THEN
4235 CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv)
4237 CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv)
4240 dgamt=1/(2.*amtau)*amplit*phspac
4252 SUBROUTINE dam4pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,AMPLIT,HV)
4264 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4265 * ,ampiz,ampi,amro,gamro,ama1,gama1
4266 * ,amk,amkz,amkst,gamkst
4268 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4269 * ,ampiz,ampi,amro,gamro,ama1,gama1
4270 * ,amk,amkz,amkst,gamkst
4271 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4272 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
4273 REAL hv(4),pt(4),pn(4),pim1(4),pim2(4),pim3(4),pim4(4)
4274 REAL pivec(4),piaks(4),hvm(4)
4275 COMPLEX hadcur(4),form1,form2,form3,form4,form5
4276 EXTERNAL form1,form2,form3,form4,form5
4277 DATA pi /3.141592653589793238462643/
4281 CALL curr_cleo(mnum,pim1,pim2,pim3,pim4,hadcur)
4285 CALL clvec(hadcur,pn,pivec)
4286 CALL claxi(hadcur,pn,piaks)
4287 CALL clnut(hadcur,brakm,hvm)
4289 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4290 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4291 amplit=(ccabib*gfermi)**2*brak/2.
4294 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4295 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4301 SUBROUTINE dam5pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,PIM5,AMPLIT,HV)
4313 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4314 * ,ampiz,ampi,amro,gamro,ama1,gama1
4315 * ,amk,amkz,amkst,gamkst
4317 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4318 * ,ampiz,ampi,amro,gamro,ama1,gama1
4319 * ,amk,amkz,amkst,gamkst
4320 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4321 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
4322 REAL hv(4),pt(4),pn(4),pim1(4),pim2(4),pim3(4),pim4(4),pim5(4)
4323 REAL pivec(4),piaks(4),hvm(4)
4324 COMPLEX hadcur(4),form1,form2,form3,form4,form5
4325 EXTERNAL form1,form2,form3,form4,form5
4326 DATA pi /3.141592653589793238462643/
4330 CALL curr5(mnum,pim1,pim2,pim3,pim4,pim5,hadcur)
4334 CALL clvec(hadcur,pn,pivec)
4335 CALL claxi(hadcur,pn,piaks)
4336 CALL clnut(hadcur,brakm,hvm)
4338 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4339 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4340 amplit=(ccabib*gfermi)**2*brak/2.
4343 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4344 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4350 SUBROUTINE dph5pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4355 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4356 * ,ampiz,ampi,amro,gamro,ama1,gama1
4357 * ,amk,amkz,amkst,gamkst
4359 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4360 * ,ampiz,ampi,amro,gamro,ama1,gama1
4363 * ,amk,amkz,amkst,gamkst
4364 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4365 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
4366 parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
4368 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4371 CHARACTER names(nmode)*31
4372 REAL hv(4),pt(4),pn(4),paa(4),pmult(4,9)
4373 REAL*4 pr(4),pi1(4),pi2(4),pi3(4),pi4(4),pi5(4)
4374 REAL*8 amp1,amp2,amp3,amp4,amp5,ams1,ams2,amom,gamom
4375 REAL*8 am5sq,am4sq,am3sq,am2sq,am5,am4,am3
4377 REAL*8 gg1,gg2,gg3,ff1,ff2,ff3,ff4,alp,alp1,alp2
4384 DATA pi /3.141592653589793238462643/
4391 bwign(xm,am,gamma)=xm**2/cmplx(xm**2-am**2,gamma*am)
4402 IF (jnpi.EQ.23.OR.jnpi.EQ.24) gamom= 0.0085
4407 phspac=1./2**29/pi**14
4410 amp1=dcdmas(idffin(1,jnpi))
4411 amp2=dcdmas(idffin(2,jnpi))
4412 amp3=dcdmas(idffin(3,jnpi))
4413 amp4=dcdmas(idffin(4,jnpi))
4414 amp5=dcdmas(idffin(5,jnpi))
4428 IF (rrr(11).GT.proba2)
THEN
4431 ams1=(amp1+amp2+amp3+amp4+amp5)**2
4432 ams2=(amtau-amnuta)**2
4433 alp1=atan((ams1-ama2**2)/ama2/gama2)
4434 alp2=atan((ams2-ama2**2)/ama2/gama2)
4435 am5sq=ams1+ rr1*(ams2-ams1)
4441 ams1=(amp1+amp2+amp3+amp4+amp5)**2
4442 ams2=(amtau-amnuta)**2
4443 alp1=atan((ams1-ama2**2)/ama2/gama2)
4444 alp2=atan((ams2-ama2**2)/ama2/gama2)
4445 alp=alp1+rr1*(alp2-alp1)
4446 am5sq =ama2**2+ama2*gama2*tan(alp)
4450 gg5=((am5sq-ama2**2)**2+(ama2*gama2)**2)/(ama2*gama2)
4452 phspac=phspac/(proba2/gg5+(1d0-proba2)/(ams2-ams1) )
4457 ams1=(amp2+amp3+amp4+amp5)**2
4459 am4sq=ams1+ rr1*(ams2-ams1)
4465 IF (rrr(12).LT.probom)
THEN
4468 ams1=(amp2+amp3+amp4)**2
4470 alp1=atan((ams1-amom**2)/amom/gamom)
4471 alp2=atan((ams2-amom**2)/amom/gamom)
4472 alp=alp1+rr1*(alp2-alp1)
4473 am3sq =amom**2+amom*gamom*tan(alp)
4478 ams1=(amp2+amp3+amp4)**2
4480 alp1=atan((ams1-amom**2)/amom/gamom)
4481 alp2=atan((ams2-amom**2)/amom/gamom)
4483 am3sq=ams1+ rr1*(ams2-ams1)
4488 gg2=((am3sq-amom**2)**2+(amom*gamom)**2)/(amom*gamom)
4490 gg2=1d0/(probom/gg2+(1d0-probom)/(ams2-ams1))
4497 am2sq=ams1+ rr2*(ams2-ams1)
4503 enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4504 enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4505 ppi= enq1**2-amp3**2
4506 pppi=sqrt(abs(enq1**2-amp3**2))
4507 ff1=(4*pi)*(2*pppi/am2)
4509 call sphera(pppi,pi3)
4520 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4521 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4522 ppi = pr(4)**2-am2**2
4526 pi2(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4529 ff2=(4*pi)*(2*pr(3)/am3)
4531 exe=(pr(4)+pr(3))/am2
4532 call bostr3(exe,pi3,pi3)
4533 call bostr3(exe,pi4,pi4)
4536 thet =acos(-1.+2*rr3)
4538 call rotpol(thet,phi,pi2)
4539 call rotpol(thet,phi,pi3)
4540 call rotpol(thet,phi,pi4)
4546 pr(4)=1./(2*am4)*(am4**2+am3**2-amp5**2)
4547 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4548 ppi = pr(4)**2-am3**2
4552 pi5(4)=1./(2*am4)*(am4**2-am3**2+amp5**2)
4555 ff3=(4*pi)*(2*pr(3)/am4)
4557 exe=(pr(4)+pr(3))/am3
4558 call bostr3(exe,pi2,pi2)
4559 call bostr3(exe,pi3,pi3)
4560 call bostr3(exe,pi4,pi4)
4563 thet =acos(-1.+2*rr3)
4565 call rotpol(thet,phi,pi2)
4566 call rotpol(thet,phi,pi3)
4567 call rotpol(thet,phi,pi4)
4568 call rotpol(thet,phi,pi5)
4574 pr(4)=1./(2*am5)*(am5**2+am4**2-amp1**2)
4575 pr(3)= sqrt(abs(pr(4)**2-am4**2))
4576 ppi = pr(4)**2-am4**2
4580 pi1(4)=1./(2*am5)*(am5**2-am4**2+amp1**2)
4583 ff4=(4*pi)*(2*pr(3)/am5)
4585 exe=(pr(4)+pr(3))/am4
4586 call bostr3(exe,pi2,pi2)
4587 call bostr3(exe,pi3,pi3)
4588 call bostr3(exe,pi4,pi4)
4589 call bostr3(exe,pi5,pi5)
4592 thet =acos(-1.+2*rr3)
4594 call rotpol(thet,phi,pi1)
4595 call rotpol(thet,phi,pi2)
4596 call rotpol(thet,phi,pi3)
4597 call rotpol(thet,phi,pi4)
4598 call rotpol(thet,phi,pi5)
4607 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5sq)
4608 paa(3)= sqrt(abs(paa(4)**2-am5sq))
4609 ppi = paa(4)**2-am5sq
4610 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4614 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am5**2)
4617 phspac=phspac * gg1*gg2*gg3*ff1*ff2*ff3*ff4
4621 exe=(paa(4)+paa(3))/am5
4622 call bostr3(exe,pi1,pi1)
4623 call bostr3(exe,pi2,pi2)
4624 call bostr3(exe,pi3,pi3)
4625 call bostr3(exe,pi4,pi4)
4626 call bostr3(exe,pi5,pi5)
4634 qxn=paa(4)*pn(4)-paa(1)*pn(1)-paa(2)*pn(2)-paa(3)*pn(3)
4635 brak=2*(gv**2+ga**2)*(2*pxq*qxn+am5sq*pxn)
4636 & -6*(gv**2-ga**2)*amtau*amnuta*am5sq
4637 fompp = cabs(bwign(am3,amom,gamom))**2
4642 amplit=ccabib**2*gfermi**2/2. * brak
4643 amplit = amplit * fompp * fnorm
4656 $ CALL dam5pi(jnpi,pt,pn,pi1,pi2,pi3,pi4,pi5,amplit,hv)
4657 dgamt=1/(2.*amtau)*amplit*phspac
4674 SUBROUTINE dphsrk(DGAMT,HV,PN,PR,PMULT,INUM)
4680 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4681 * ,ampiz,ampi,amro,gamro,ama1,gama1
4682 * ,amk,amkz,amkst,gamkst
4684 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4685 * ,ampiz,ampi,amro,gamro,ama1,gama1
4686 * ,amk,amkz,amkst,gamkst
4687 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4688 REAL*4 gfermi,gv,ga,ccabib,scabib,gamel
4689 REAL hv(4),pt(4),pn(4),pr(4),pkc(4),pkz(4),qq(4),pmult(4,9)
4693 DATA pi /3.141592653589793238462643/
4697 phspac=1./2**11/pi**5
4705 ams2=(amtau-amnuta)**2
4708 amx2=ams1+ rr1(1)*(ams2-ams1)
4710 phspac=phspac*(ams2-ams1)
4728 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
4729 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
4733 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
4735 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
4738 enq1=(amx2+amk**2-amkz**2)/(2.*amx)
4739 enq2=(amx2-amk**2+amkz**2)/(2.*amx)
4740 pppi=sqrt((enq1-amk)*(enq1+amk))
4741 phspac=phspac*(4*pi)*(2*pppi/amx)
4743 CALL sphera(pppi,pkc)
4749 exe=(pr(4)+pr(3))/amx
4751 CALL bostr3(exe,pkc,pkc)
4752 CALL bostr3(exe,pkz,pkz)
4754 30 qq(i)=pkc(i)-pkz(i)
4756 pksd =pr(4)*pr(4)-pr(3)*pr(3)-pr(2)*pr(2)-pr(1)*pr(1)
4757 qqpks=pr(4)* qq(4)-pr(3)* qq(3)-pr(2)* qq(2)-pr(1)* qq(1)
4759 31 qq(i)=qq(i)-pr(i)*qqpks/pksd
4762 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
4764 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
4765 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
4766 & +(gv**2-ga**2)*amtau*amnuta*qq2
4767 amplit=(gfermi*ccabib)**2*brak*2*fpirk(amx)
4768 dgamt=1/(2.*amtau)*amplit*phspac
4770 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
4781 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4782 * ,ampiz,ampi,amro,gamro,ama1,gama1
4783 * ,amk,amkz,amkst,gamkst
4785 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
4786 * ,ampiz,ampi,amro,gamro,ama1,gama1
4787 * ,amk,amkz,amkst,gamkst
4790 fpirk=cabs(fpikm(w,amk,amkz))**2
4793 COMPLEX FUNCTION fpikmk(W,XM1,XM2)
4798 REAL rom,rog,rom1,rog1,beta1,pi,pim,s,w
4803 IF (init.EQ.0 )
THEN
4816 fpikmk=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
4827 SUBROUTINE dwrph(KTO,PHX)
4831 IMPLICIT REAL*8 (a-h,o-z)
4842 IF (qhot(4).GT.1.e-5) CALL dwluph(kto,qhot)
4845 SUBROUTINE dwluph(KTO,PHOT)
4859 COMMON /taupos/ np1,np2
4863 IF (phot(4).LE.0.0)
RETURN
4866 IF((kto.EQ. 1).OR.(kto.EQ.11))
THEN
4873 IF(ktos.GT.10) ktos=ktos-10
4875 CALL tralo4(ktos,phot,phot,am)
4876 CALL filhep(0,1,22,nps,nps,0,0,phot,0.0,.true.)
4881 SUBROUTINE dwluel(KTO,ISGN,PNU,PWB,PEL,PNE)
4893 REAL pnu(4),pwb(4),pel(4),pne(4)
4895 COMMON /taupos/ np1,np2
4906 CALL tralo4(kto,pnu,pnu,am)
4907 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4910 CALL tralo4(kto,pwb,pwb,am)
4914 CALL tralo4(kto,pel,pel,am)
4915 CALL filhep(0,1,11*isgn,nps,nps,0,0,pel,am,.false.)
4918 CALL tralo4(kto,pne,pne,am)
4919 CALL filhep(0,1,-12*isgn,nps,nps,0,0,pne,am,.true.)
4923 SUBROUTINE dwlumu(KTO,ISGN,PNU,PWB,PMU,PNM)
4935 REAL pnu(4),pwb(4),pmu(4),pnm(4)
4937 COMMON /taupos/ np1,np2
4948 CALL tralo4(kto,pnu,pnu,am)
4949 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4952 CALL tralo4(kto,pwb,pwb,am)
4956 CALL tralo4(kto,pmu,pmu,am)
4957 CALL filhep(0,1,13*isgn,nps,nps,0,0,pmu,am,.false.)
4960 CALL tralo4(kto,pnm,pnm,am)
4961 CALL filhep(0,1,-14*isgn,nps,nps,0,0,pnm,am,.true.)
4965 SUBROUTINE dwlupi(KTO,ISGN,PPI,PNU)
4976 COMMON /taupos/ np1,np2
4986 CALL tralo4(kto,pnu,pnu,am)
4987 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4990 CALL tralo4(kto,ppi,ppi,am)
4991 CALL filhep(0,1,-211*isgn,nps,nps,0,0,ppi,am,.true.)
4995 SUBROUTINE dwluro(KTO,ISGN,PNU,PRHO,PIC,PIZ)
5007 REAL pnu(4),prho(4),pic(4),piz(4)
5009 COMMON /taupos/ np1,np2
5020 CALL tralo4(kto,pnu,pnu,am)
5021 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5024 CALL tralo4(kto,prho,prho,am)
5025 CALL filhep(0,2,-213*isgn,nps,nps,0,0,prho,am,.true.)
5028 CALL tralo4(kto,pic,pic,am)
5029 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pic,am,.true.)
5032 CALL tralo4(kto,piz,piz,am)
5033 CALL filhep(0,1,111,-2,-2,0,0,piz,am,.true.)
5037 SUBROUTINE dwluaa(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
5050 REAL pnu(4),paa(4),pim1(4),pim2(4),pipl(4)
5052 COMMON /taupos/ np1,np2
5063 CALL tralo4(kto,pnu,pnu,am)
5064 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5067 CALL tralo4(kto,paa,paa,am)
5068 CALL filhep(0,1,-20213*isgn,nps,nps,0,0,paa,am,.true.)
5076 CALL tralo4(kto,pim2,pim2,am)
5077 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pim2,am,.true.)
5080 CALL tralo4(kto,pim1,pim1,am)
5081 CALL filhep(0,1,-211*isgn,-2,-2,0,0,pim1,am,.true.)
5084 CALL tralo4(kto,pipl,pipl,am)
5085 CALL filhep(0,1, 211*isgn,-3,-3,0,0,pipl,am,.true.)
5087 ELSE IF (jaa.EQ.2)
THEN
5092 CALL tralo4(kto,pim2,pim2,am)
5093 CALL filhep(0,1,111,-1,-1,0,0,pim2,am,.true.)
5096 CALL tralo4(kto,pim1,pim1,am)
5097 CALL filhep(0,1,111,-2,-2,0,0,pim1,am,.true.)
5100 CALL tralo4(kto,pipl,pipl,am)
5101 CALL filhep(0,1,-211*isgn,-3,-3,0,0,pipl,am,.true.)
5107 SUBROUTINE dwlukk (KTO,ISGN,PKK,PNU)
5117 COMMON /taupos/ np1,np2
5129 CALL tralo4(kto,pnu,pnu,am)
5130 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5133 CALL tralo4(kto,pkk,pkk,am)
5134 CALL filhep(0,1,-321*isgn,nps,nps,0,0,pkk,am,.true.)
5138 SUBROUTINE dwluks(KTO,ISGN,PNU,PKS,PKK,PPI,JKST)
5139 COMMON / taukle / bra1,brk0,brk0b,brks
5140 REAL*4 bra1,brk0,brk0b,brks
5152 REAL pnu(4),pks(4),pkk(4),ppi(4),xio(1)
5153 COMMON /taupos/ np1,np2
5164 CALL tralo4(kto,pnu,pnu,am)
5165 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5168 CALL tralo4(kto,pks,pks,am)
5169 CALL filhep(0,1,-323*isgn,nps,nps,0,0,pks,am,.true.)
5177 CALL tralo4(kto,ppi,ppi,am)
5178 CALL filhep(0,1,-211*isgn,-1,-1,0,0,ppi,am,.true.)
5181 IF (isgn.EQ.-1) bran=brk0
5184 IF(xio(1).GT.bran)
THEN
5190 CALL tralo4(kto,pkk,pkk,am)
5191 CALL filhep(0,1,k0type,-2,-2,0,0,pkk,am,.true.)
5193 ELSE IF(jkst.EQ.20)
THEN
5198 CALL tralo4(kto,ppi,ppi,am)
5199 CALL filhep(0,1,111,-1,-1,0,0,ppi,am,.true.)
5202 CALL tralo4(kto,pkk,pkk,am)
5203 CALL filhep(0,1,-321*isgn,-2,-2,0,0,pkk,am,.true.)
5209 SUBROUTINE dwlnew(KTO,ISGN,PNU,PWB,PNPI,MODE)
5219 parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
5221 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5224 COMMON /taupos/ np1,np2
5225 CHARACTER names(nmode)*31
5226 REAL pnu(4),pwb(4),pnpi(4,9)
5238 CALL tralo4(kto,pnu,pnu,am)
5239 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5242 CALL tralo4(kto,pwb,pwb,am)
5243 CALL filhep(0,1,-24*isgn,nps,nps,0,0,pwb,am,.true.)
5251 kfpi=lunpik(idffin(i,jnpi),-isgn)
5258 CALL tralo4(kto,ppi,ppi,am)
5259 CALL filhep(0,1,kfpi,-i,-i,0,0,ppi,am,.true.)
5272 IMPLICIT REAL*8 (a-h,o-z)
5274 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5276 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5288 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5289 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5298 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5299 DATA pi /3.141592653589793238462643d0/
5301 IF(abs(y).LT.abs(x))
THEN
5303 IF(x.LE.0d0) the=pi-the
5305 the=acos(x/sqrt(x**2+y**2))
5316 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5317 DATA pi /3.141592653589793238462643d0/
5319 IF(abs(y).LT.abs(x))
THEN
5321 IF(x.LE.0d0) the=pi-the
5323 the=acos(x/sqrt(x**2+y**2))
5325 IF(y.LT.0d0) the=2d0*pi-the
5328 SUBROUTINE rotod1(PH1,PVEC,QVEC)
5333 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5334 dimension pvec(4),qvec(4),rvec(4)
5342 qvec(2)= cs*rvec(2)-sn*rvec(3)
5343 qvec(3)= sn*rvec(2)+cs*rvec(3)
5347 SUBROUTINE rotod2(PH1,PVEC,QVEC)
5352 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5353 dimension pvec(4),qvec(4),rvec(4)
5360 qvec(1)= cs*rvec(1)+sn*rvec(3)
5362 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5366 SUBROUTINE rotod3(PH1,PVEC,QVEC)
5371 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5373 dimension pvec(4),qvec(4),rvec(4)
5379 qvec(1)= cs*rvec(1)-sn*rvec(2)
5380 qvec(2)= sn*rvec(1)+cs*rvec(2)
5384 SUBROUTINE bostr3(EXE,PVEC,QVEC)
5390 REAL*4 pvec(4),qvec(4),rvec(4)
5403 SUBROUTINE bostd3(EXE,PVEC,QVEC)
5409 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5410 dimension pvec(4),qvec(4),rvec(4)
5424 SUBROUTINE rotor1(PH1,PVEC,QVEC)
5429 REAL*4 pvec(4),qvec(4),rvec(4)
5437 qvec(2)= cs*rvec(2)-sn*rvec(3)
5438 qvec(3)= sn*rvec(2)+cs*rvec(3)
5441 SUBROUTINE rotor2(PH1,PVEC,QVEC)
5446 IMPLICIT REAL*4(a-h,o-z)
5447 REAL*4 pvec(4),qvec(4),rvec(4)
5454 qvec(1)= cs*rvec(1)+sn*rvec(3)
5456 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5459 SUBROUTINE rotor3(PHI,PVEC,QVEC)
5464 REAL*4 pvec(4),qvec(4),rvec(4)
5470 qvec(1)= cs*rvec(1)-sn*rvec(2)
5471 qvec(2)= sn*rvec(1)+cs*rvec(2)
5475 SUBROUTINE spherd(R,X)
5480 REAL*8 r,x(4),pi,costh,sinth
5482 DATA pi /3.141592653589793238462643d0/
5486 sinth=sqrt(1 -costh**2)
5487 x(1)=r*sinth*cos(2*pi*rrr(2))
5488 x(2)=r*sinth*sin(2*pi*rrr(2))
5492 SUBROUTINE rotpox(THET,PHI,PP)
5493 IMPLICIT REAL*8 (a-h,o-z)
5501 CALL rotod2(thet,pp,pp)
5502 CALL rotod3( phi,pp,pp)
5505 SUBROUTINE sphera(R,X)
5513 DATA pi /3.141592653589793238462643/
5517 sinth=sqrt(1.-costh**2)
5518 x(1)=r*sinth*cos(2*pi*rrr(2))
5519 x(2)=r*sinth*sin(2*pi*rrr(2))
5523 SUBROUTINE rotpol(THET,PHI,PP)
5530 CALL rotor2(thet,pp,pp)
5531 CALL rotor3( phi,pp,pp)
5534 #include "../randg/tauola-random.h"
5537 IMPLICIT REAL*8(a-h,o-z)
5540 IF(x .LT.-1.0) go to 1
5541 IF(x .LE. 0.5) go to 2
5542 IF(x .EQ. 1.0) go to 3
5543 IF(x .LE. 2.0) go to 4
5547 z=z-0.5* log(abs(x))**2
5553 3 dilogt=1.64493406684822
5557 z=1.64493406684822 - log(x)* log(abs(t))
5558 5 y=2.66666666666666 *t+0.66666666666666
5559 b= 0.00000 00000 00001
5560 a=y*b +0.00000 00000 00004
5561 b=y*a-b+0.00000 00000 00011
5562 a=y*b-a+0.00000 00000 00037
5563 b=y*a-b+0.00000 00000 00121
5564 a=y*b-a+0.00000 00000 00398
5565 b=y*a-b+0.00000 00000 01312
5566 a=y*b-a+0.00000 00000 04342
5567 b=y*a-b+0.00000 00000 14437
5568 a=y*b-a+0.00000 00000 48274
5569 b=y*a-b+0.00000 00001 62421
5570 a=y*b-a+0.00000 00005 50291
5571 b=y*a-b+0.00000 00018 79117
5572 a=y*b-a+0.00000 00064 74338
5573 b=y*a-b+0.00000 00225 36705
5574 a=y*b-a+0.00000 00793 87055
5575 b=y*a-b+0.00000 02835 75385
5576 a=y*b-a+0.00000 10299 04264
5577 b=y*a-b+0.00000 38163 29463
5578 a=y*b-a+0.00001 44963 00557
5579 b=y*a-b+0.00005 68178 22718
5580 a=y*b-a+0.00023 20021 96094
5581 b=y*a-b+0.00100 16274 96164
5582 a=y*b-a+0.00468 63619 59447
5583 b=y*a-b+0.02487 93229 24228
5584 a=y*b-a+0.16607 30329 27855
5585 a=y*a-b+1.93506 43008 6996