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)
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)
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,pdum,pdum1,pdum2)
141 CALL dadmro(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4)
142 CALL dadmaa(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5,jdum)
143 CALL dadmkk(-1,idum,pdum,pdum1,pdum2)
144 CALL dadmks(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,jdum)
145 CALL dadnew(-1,idum,hdum,pdum1,pdum2,pdumx,jdum)
151 ELSEIF(kto.EQ.1)
THEN
155 IF(iwarm.EQ.0)
GOTO 902
158 #elif defined (ALEPH)
163 CALL dekay1(0,h,isgn)
164 ELSEIF(kto.EQ.2)
THEN
168 IF(iwarm.EQ.0)
GOTO 902
171 #elif defined (ALEPH)
176 CALL dekay2(0,h,isgn)
177 ELSEIF(kto.EQ.11)
THEN
182 CALL dekay1(1,h,isgn)
183 ELSEIF(kto.EQ.12)
THEN
188 CALL dekay2(1,h,isgn)
189 ELSEIF(kto.EQ.100)
THEN
191 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
192 CALL dadmel( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
193 CALL dadmmu( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
194 CALL dadmpi( 1,idum,pdum,pdum1,pdum2)
195 CALL dadmro( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4)
196 CALL dadmaa( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5,jdum)
197 CALL dadmkk( 1,idum,pdum,pdum1,pdum2)
198 CALL dadmks( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,jdum)
199 CALL dadnew( 1,idum,hdum,pdum1,pdum2,pdumx,jdum)
200 WRITE(iout,7010) nev1,nev2,nevtot
201 WRITE(iout,7011) (nevdec(i),gampmc(i),gamper(i),i= 1,7)
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.5 ******',9x,1h*,
217 $ /,
' *', 25x,
'*DEC 1993; ALEPH fixes introd. dec 98 *',9x,1h*,
218 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
219 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
220 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
221 $ /,
' *', 25x,
'Physics initialization by ALEPH collab ',9x,1h*,
222 $ /,
' *', 25x,
'it is suggested to use this version ',9x,1h*,
223 $ /,
' *', 25x,
' with the help of the collab. advice ',9x,1h*,
224 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
225 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
227 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.6 ******',9x,1h*,
228 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
229 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
230 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
231 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
232 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
233 $ /,
' *', 25x,
' Physics initialization by CLEO collab ',9x,1h*,
234 $ /,
' *', 25x,
' see Alain Weinstein www home page: ',9x,1h*,
235 $ /,
' *', 25x,
'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
236 $ /,
' *', 25x,
'/korb_doc.html#files ',9x,1h*,
237 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
238 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
240 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.6 ******',9x,1h*,
241 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
242 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
243 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
244 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
245 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
246 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
247 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
248 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
249 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
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.5 ******',9x,1h*,
258 $ /,
' *', 25x,
'***********DECEMBER 1993***************',9x,1h*,
260 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.6 ******',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)
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
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)
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
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)
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)
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
574 IF(iwarm.EQ.0)
GOTO 902
577 CALL dexay1(kto,jak1,jakp,pol,isgn)
578 ELSEIF(kto.EQ.2)
THEN
583 IF(iwarm.EQ.0)
GOTO 902
584 isgn=-idff/iabs(idff)
586 CALL dexay1(kto,jak2,jakm,pol,isgn)
587 ELSEIF(kto.EQ.100)
THEN
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.5 ******',9x,1h*,
611 $ /,
' *', 25x,
'*DEC 1993; ALEPH fixes introd. dec 98 *',9x,1h*,
612 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
613 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
614 $ /,
' *', 25x,
'Physics initialization by ALEPH collab ',9x,1h*,
615 $ /,
' *', 25x,
'it is suggested to use this version ',9x,1h*,
616 $ /,
' *', 25x,
' with the help of the collab. advice ',9x,1h*,
617 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
618 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
620 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.6 ******',9x,1h*,
621 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
622 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
623 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
624 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
625 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
626 $ /,
' *', 25x,
' Physics initialization by CLEO collab ',9x,1h*,
627 $ /,
' *', 25x,
' see Alain Weinstein www home page: ',9x,1h*,
628 $ /,
' *', 25x,
'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
629 $ /,
' *', 25x,
'/korb_doc.html#files ',9x,1h*,
630 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
631 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
633 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.6 ******',9x,1h*,
634 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
635 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
636 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
637 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
638 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
639 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
640 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
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*
650 7010
FORMAT(///1x,15(5h*****)
652 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.5 ******',9x,1h*,
653 $ /,
' *', 25x,
'***********DECEMBER 1993***************',9x,1h*,
655 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.6 ******',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)
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)
757 REAL POL(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4),PH(4),RN(1)
763 CALL dadmel( -1,isgn,hv,pnu,pwb,q1,q2,ph)
766 ELSEIF(mode.EQ. 0)
THEN
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.
774 IF(rn(1).GT.wt)
GOTO 300
776 ELSEIF(mode.EQ. 1)
THEN
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)
796 COMMON / inout / inut,iout
797 REAL POL(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4),PH(4),RN(1)
803 CALL dadmmu( -1,isgn,hv,pnu,pwb,q1,q2,ph)
806 ELSEIF(mode.EQ. 0)
THEN
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.
814 IF(rn(1).GT.wt)
GOTO 300
816 ELSEIF(mode.EQ. 1)
THEN
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)
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/
869 CALL dphsel(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
870 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
874 ELSEIF(mode.EQ. 0)
THEN
877 IF(iwarm.EQ.0)
GOTO 902
879 CALL dphsel(wt,hv,pnu,pwb,q1,q2,phx)
885 IF(wt.GT.wtmax) nevovr=nevovr+1
886 IF(rn*wtmax.GT.wt)
GOTO 300
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
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)
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/
971 CALL dphsmu(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
972 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
976 ELSEIF(mode.EQ. 0)
THEN
979 IF(iwarm.EQ.0)
GOTO 902
981 CALL dphsmu(wt,hv,pnu,pwb,q1,q2,phx)
987 IF(wt.GT.wtmax) nevovr=nevovr+1
988 IF(rn*wtmax.GT.wt)
GOTO 300
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
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
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)
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)
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)
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/
1111 xlam(x,y,z)=sqrt((x-y-z)**2-4.0*y*z)
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)
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
1168 am2sq=ams1+ rr2*(ams2-ams1)
1170 phspac=phspac*(ams2-ams1)
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)
1178 CALL spherd(pppi,xn)
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)
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)
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)
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)
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)
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)
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
1329 rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(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)
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)
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)
1411 rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
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
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)
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 )
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
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)
1484 REAL POL(4),HV(4),PNU(4),PPI(4),RN(1)
1488 CALL dadmpi(-1,isgn,hv,ppi,pnu)
1491 ELSEIF(mode.EQ. 0)
THEN
1494 CALL dadmpi( 0,isgn,hv,ppi,pnu)
1495 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1498 IF(rn(1).GT.wt)
GOTO 300
1500 ELSEIF(mode.EQ. 1)
THEN
1502 CALL dadmpi( 1,isgn,hv,ppi,pnu)
1508 SUBROUTINE dadmpi(MODE,ISGN,HV,PPI,PNU)
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/
1528 ELSEIF(mode.EQ. 0)
THEN
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)
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
1552 IF(nevtot.EQ.0)
RETURN
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
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)
1588 COMMON / inout / inut,iout
1589 REAL POL(4),HV(4),PRO(4),PNU(4),PIC(4),PIZ(4),RN(1)
1595 CALL dadmro( -1,isgn,hv,pnu,pro,pic,piz)
1599 ELSEIF(mode.EQ. 0)
THEN
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.
1609 IF(rn(1).GT.wt)
GOTO 300
1611 ELSEIF(mode.EQ. 1)
THEN
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)
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/
1655 CALL dphsro(wt,hv,pdum1,pdum2,pdum3,pdum4)
1656 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1661 ELSEIF(mode.EQ. 0)
THEN
1664 IF(iwarm.EQ.0)
GOTO 902
1665 CALL dphsro(wt,hv,pnu,pro,pic,piz)
1672 IF(wt.GT.wtmax) nevovr=nevovr+1
1673 IF(rn*wtmax.GT.wt)
GOTO 300
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
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
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)
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/
1743 phspac=1./2**11/pi**5
1750 ams1=(ampi+ampiz)**2
1751 ams2=(amtau-amnuta)**2
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)
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)
1792 CALL sphera(pppi,pic)
1798 exe=(pr(4)+pr(3))/amx
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)
1828 COMMON / inout / inut,iout
1829 REAL POL(4),HV(4),PAA(4),PNU(4),PIM1(4),PIM2(4),PIPL(4),RN(1)
1835 CALL dadmaa( -1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1838 ELSEIF(mode.EQ. 0)
THEN
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.
1846 IF(rn(1).GT.wt)
GOTO 300
1848 ELSEIF(mode.EQ. 1)
THEN
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)
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/
1893 CALL dphsaa(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5,jaa)
1894 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1898 ELSEIF(mode.EQ. 0)
THEN
1901 IF(iwarm.EQ.0)
GOTO 902
1902 CALL dphsaa(wt,hv,pnu,paa,pim1,pim2,pipl,jaa)
1911 sswt=sswt+dble(wt)**2
1916 IF(wt.GT.wtmax) nevovr=nevovr+1
1917 IF(rn*wtmax.GT.wt)
GOTO 300
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
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
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)
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)
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)
2009 REAL POL(4),HV(4),PNU(4),PKK(4),RN(1)
2013 CALL dadmkk(-1,isgn,hv,pkk,pnu)
2016 ELSEIF(mode.EQ. 0)
THEN
2019 CALL dadmkk( 0,isgn,hv,pkk,pnu)
2020 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
2023 IF(rn(1).GT.wt)
GOTO 300
2025 ELSEIF(mode.EQ. 1)
THEN
2027 CALL dadmkk( 1,isgn,hv,pkk,pnu)
2033 SUBROUTINE dadmkk(MODE,ISGN,HV,PKK,PNU)
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/
2062 ELSEIF(mode.EQ. 0)
THEN
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)
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
2086 IF(nevtot.EQ.0)
RETURN
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
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)
2128 COMMON / inout / inut,iout
2129 REAL POL(4),HV(4),PKS(4),PNU(4),PKK(4),PPI(4),RN(1)
2136 CALL dadmks( -1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2140 ELSEIF(mode.EQ. 0)
THEN
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.
2150 IF(rn(1).GT.wt)
GOTO 300
2152 ELSEIF(mode.EQ. 1)
THEN
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)
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/
2200 CALL dphsks(wt,hv,pdum1,pdum2,pdum3,pdum4,jkst)
2201 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
2206 ELSEIF(mode.EQ. 0)
THEN
2208 IF(iwarm.EQ.0)
GOTO 902
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
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
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
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)
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/
2312 phspac=1./2**11/pi**5
2324 ams2=(amtau-amnuta)**2
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)
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)
2357 CALL sphera(pppi,ppi)
2362 pkk(4)=( amx**2+amkz**2-ampi**2 ) / ( 2*amx )
2363 exe=(pks(4)+pks(3))/amx
2365 CALL bostr3(exe,ppi,ppi)
2366 CALL bostr3(exe,pkk,pkk)
2368 30 qq(i)=ppi(i)-pkk(i)
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
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
2394 ELSEIF(jkst.EQ.20)
THEN
2398 ams2=(amtau-amnuta)**2
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)
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)
2434 CALL sphera(pppi,ppi)
2439 pkk(4)=( amx**2+amk**2-ampiz**2 ) / ( 2*amx )
2440 exe=(pks(4)+pks(3))/amx
2442 CALL bostr3(exe,ppi,ppi)
2443 CALL bostr3(exe,pkk,pkk)
2445 60 qq(i)=pkk(i)-ppi(i)
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
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)
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/
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/
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)
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)
2600 qxn=pr(4)*pn(4)-pr(1)*pn(1)-pr(2)*pn(2)-pr(3)*pn(3)
2606 brak=2*(gv**2+ga**2)*(2*pxq*qxn+amx2*pxn)
2607 & -6*(gv**2-ga**2)*amtau*amnuta*amx2
2610 amplit=ccabib**2*gfermi**2/2. * brak * amx2*sigee(amx2,jnpi)
2611 dgamt=1./(2.*amtau)*amplit*phspac
2615 phspac = 1./2.**(6*nd-7) /pi**(3*nd-4)
2616 phsmax = 1./dpar(nd-2)
2623 pv(5,nd)=ampik(nd,jnpi)
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))
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
2658 223 ams1=ams1+ampik(jl,jnpi)
2660 amx =(amx-ampik(il,jnpi))
2662 phsmax=phsmax * (ams2-ams1)
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)
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
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))
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)
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/
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)
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)
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
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)
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)
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)
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))
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
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)
3062 IF (ichan.LE.2)
THEN
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)
3083 am2sq=ams1+ rr2*(ams2-ams1)
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))
3097 phf1=(4*pi)*(2*pppi/am2)
3103 CALL sphera(pppi,pipl)
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)
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)
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)
3162 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am3**2)
3168 alp1=atan((ams1-amra**2)/amra/gamra)
3169 alp2=atan((ams2-amra**2)/amra/gamra)
3170 xpro = (pim1(3)+pipl(3))**2
3171 $ +(pim1(2)+pipl(2))**2+(pim1(1)+pipl(1))**2
3172 am2sq=-xpro+(pim1(4)+pipl(4))**2
3174 ff1 = ((am2sq-amra**2)**2+(amra*gamra)**2)/(amra*gamra)
3175 ff1 =ff1 *(alp2-alp1)
3177 gg1 = (4*pi)*(xlam(am2sq,amp2**2,amp3**2)/am2sq)
3179 gg1 =gg1 *(4*pi)*sqrt(4*xpro/am3sq)
3180 xjaje=gg1*(ams2-ams1)
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)
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)
3229 CALL dampog(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
3246 dgamt=1/(2.*amtau)*amplit*phspac
3248 SUBROUTINE dampaa(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
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
3278 bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
3282 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
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))
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
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))
3319 CALL clvec(hadcur,pn,pivec)
3320 CALL claxi(hadcur,pn,piaks)
3321 CALL clnut(hadcur,brakm,hvm)
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.
3330 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3331 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
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)
3361 REAL PI,PIM,QS,QM,W,GS,MK
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)
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)
3414 REAL PI,PIM,QS,QM,W,GS
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)
3439 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
3444 IF (init.EQ.0 )
THEN
3464 fpik= (bwig(s,rom,rog)+beta1*bwig(s,rom1,rog1))
3473 fpirho=cabs(fpik(w))**2
3475 SUBROUTINE clvec(HJ,PN,PIV)
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)
3499 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3500 COMMON / idfc / idff
3502 COMPLEX HJ(4),HJC(4)
3505 det2(i,j)=aimag(hjc(i)*hj(j)-hjc(j)*hj(i))
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)
3527 20 pia(i)=pia(i)*sign
3529 SUBROUTINE clnut(HJ,B,HV)
3541 b=real( hj(4)*aimag(hj(4)) - hj(3)*aimag(hj(3))
3542 & - hj(2)*aimag(hj(2)) - hj(1)*aimag(hj(1)) )
3545 SUBROUTINE dampog(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3559 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3560 * ,ampiz,ampi,amro,gamro,ama1,gama1
3561 * ,amk,amkz,amkst,gamkst
3563 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3564 * ,ampiz,ampi,amro,gamro,ama1,gama1
3565 * ,amk,amkz,amkst,gamkst
3566 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3567 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3568 COMMON /testa1/ keya1
3569 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4)
3570 REAL PAA(4),VEC1(4),VEC2(4)
3571 REAL PIVEC(4),PIAKS(4),HVM(4)
3572 COMPLEX BWIGN,HADCUR(4),FNORM,FORMOM
3578 bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
3586 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3589 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3590 xmom =sqrt(abs( (pim2(4)+pipl(4))**2-(pim2(3)+pipl(3))**2
3591 $ -(pim2(2)+pipl(2))**2-(pim2(1)+pipl(1))**2 ))
3592 xmro2 =(pipl(1))**2 +(pipl(2))**2 +(pipl(3))**2
3594 prod1 =vec1(1)*pipl(1)
3595 prod2 =vec2(2)*pipl(2)
3596 p12 =pim1(4)*pim2(4)-pim1(1)*pim2(1)
3597 $ -pim1(2)*pim2(2)-pim1(3)*pim2(3)
3598 p1pl =pim1(4)*pipl(4)-pim1(1)*pipl(1)
3599 $ -pim1(2)*pipl(2)-pim1(3)*pipl(3)
3600 p2pl =pipl(4)*pim2(4)-pipl(1)*pim2(1)
3601 $ -pipl(2)*pim2(2)-pipl(3)*pim2(3)
3603 vec1(i)= (vec1(i)-prod1/xmro2*pipl(i))
3605 gnorm=sqrt(vec1(1)**2+vec1(2)**2+vec1(3)**2)
3607 vec1(i)= vec1(i)/gnorm
3609 vec2(1)=(vec1(2)*pipl(3)-vec1(3)*pipl(2))/sqrt(xmro2)
3610 vec2(2)=(vec1(3)*pipl(1)-vec1(1)*pipl(3))/sqrt(xmro2)
3611 vec2(3)=(vec1(1)*pipl(2)-vec1(2)*pipl(1))/sqrt(xmro2)
3612 p1vec1 =pim1(4)*vec1(4)-pim1(1)*vec1(1)
3613 $ -pim1(2)*vec1(2)-pim1(3)*vec1(3)
3614 p2vec1 =vec1(4)*pim2(4)-vec1(1)*pim2(1)
3615 $ -vec1(2)*pim2(2)-vec1(3)*pim2(3)
3616 p1vec2 =pim1(4)*vec2(4)-pim1(1)*vec2(1)
3617 $ -pim1(2)*vec2(2)-pim1(3)*vec2(3)
3618 p2vec2 =vec2(4)*pim2(4)-vec2(1)*pim2(1)
3619 $ -vec2(2)*pim2(2)-vec2(3)*pim2(3)
3621 fnorm=formom(xmaa,xmom)
3626 hadcur(i) = fnorm *(
3627 $ vec1(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3628 $ -pim2(i)*(p2vec1*p1pl-p1vec1*p2pl)
3629 $ +pipl(i)*(p2vec1*p12 -p1vec1*(ampi**2+p2pl)) )
3631 hadcur(i) = fnorm *(
3632 $ vec2(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3633 $ -pim2(i)*(p2vec2*p1pl-p1vec2*p2pl)
3634 $ +pipl(i)*(p2vec2*p12 -p1vec2*(ampi**2+p2pl)) )
3639 CALL clvec(hadcur,pn,pivec)
3640 CALL claxi(hadcur,pn,piaks)
3641 CALL clnut(hadcur,brakm,hvm)
3643 brak=brak+(gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3644 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3646 hv(i)=hv(i)-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3647 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3651 amplit=(gfermi*ccabib)**2*brak/2.
3660 SUBROUTINE damppk(MNUM,PT,PN,PIM1,PIM2,PIM3,AMPLIT,HV)
3674 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3675 * ,ampiz,ampi,amro,gamro,ama1,gama1
3676 * ,amk,amkz,amkst,gamkst
3678 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3679 * ,ampiz,ampi,amro,gamro,ama1,gama1
3680 * ,amk,amkz,amkst,gamkst
3681 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3682 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3683 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4)
3684 REAL PAA(4),VEC1(4),VEC2(4),VEC3(4),VEC4(4),VEC5(4)
3685 REAL PIVEC(4),PIAKS(4),HVM(4)
3686 REAL FNORM(0:7),COEF(1:5,0:7)
3687 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5,UROJ
3689 COMPLEX F1,F2,F3,F4,F5
3691 EXTERNAL form1,form2,form3,form4,form5
3692 DATA pi /3.141592653589793238462643/
3696 IF (icont.EQ.0)
THEN
3704 fnorm(4)=scabib/fpi/dwapi0
3709 coef(1,0)= 2.0*sqrt(2.)/3.0
3710 coef(2,0)=-2.0*sqrt(2.)/3.0
3713 coef(3,0)= 2.0*sqrt(2.)/3.0
3720 coef(1,1)=-sqrt(2.)/3.0
3721 coef(2,1)= sqrt(2.)/3.0
3726 coef(1,2)=-sqrt(2.)/3.0
3727 coef(2,2)= sqrt(2.)/3.0
3745 coef(1,4)= 1.0/sqrt(2.)/3.0
3746 coef(2,4)=-1.0/sqrt(2.)/3.0
3751 coef(1,5)=-sqrt(2.)/3.0
3752 coef(2,5)= sqrt(2.)/3.0
3774 coef(5,7)=-sqrt(2.0/3.0)
3779 10 paa(i)=pim1(i)+pim2(i)+pim3(i)
3780 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3781 xmro1 =sqrt(abs((pim3(4)+pim2(4))**2-(pim3(1)+pim2(1))**2
3782 $ -(pim3(2)+pim2(2))**2-(pim3(3)+pim2(3))**2))
3783 xmro2 =sqrt(abs((pim3(4)+pim1(4))**2-(pim3(1)+pim1(1))**2
3784 $ -(pim3(2)+pim1(2))**2-(pim3(3)+pim1(3))**2))
3785 xmro3 =sqrt(abs((pim1(4)+pim2(4))**2-(pim1(1)+pim2(1))**2
3786 $ -(pim1(2)+pim2(2))**2-(pim1(3)+pim2(3))**2))
3788 prod1 =paa(4)*(pim2(4)-pim3(4))-paa(1)*(pim2(1)-pim3(1))
3789 $ -paa(2)*(pim2(2)-pim3(2))-paa(3)*(pim2(3)-pim3(3))
3790 prod2 =paa(4)*(pim3(4)-pim1(4))-paa(1)*(pim3(1)-pim1(1))
3791 $ -paa(2)*(pim3(2)-pim1(2))-paa(3)*(pim3(3)-pim1(3))
3792 prod3 =paa(4)*(pim1(4)-pim2(4))-paa(1)*(pim1(1)-pim2(1))
3793 $ -paa(2)*(pim1(2)-pim2(2))-paa(3)*(pim1(3)-pim2(3))
3795 vec1(i)= pim2(i)-pim3(i) -paa(i)*prod1/xmaa**2
3796 vec2(i)= pim3(i)-pim1(i) -paa(i)*prod2/xmaa**2
3797 vec3(i)= pim1(i)-pim2(i) -paa(i)*prod3/xmaa**2
3798 40 vec4(i)= pim1(i)+pim2(i)+pim3(i)
3799 CALL prod5(pim1,pim2,pim3,vec5)
3804 f1 = cmplx(coef(1,mnum))*form1(mnum,xmaa**2,xmro1**2,xmro2**2)
3805 f2 = cmplx(coef(2,mnum))*form2(mnum,xmaa**2,xmro2**2,xmro1**2)
3806 f3 = cmplx(coef(3,mnum))*form3(mnum,xmaa**2,xmro3**2,xmro1**2)
3808 $cmplx(coef(4,mnum))*form4(mnum,xmaa**2,xmro1**2,xmro2**2,xmro3**2)
3809 f5 = (-1.0)*uroj/4.0/pi**2/fpi**2*
3810 $ cmplx(coef(5,mnum))*form5(mnum,xmaa**2,xmro1**2,xmro2**2)
3813 hadcur(i)= cmplx(fnorm(mnum)) * (
3814 $ cmplx(vec1(i))*f1+cmplx(vec2(i))*f2+cmplx(vec3(i))*f3+
3815 $ cmplx(vec4(i))*f4+cmplx(vec5(i))*f5)
3819 hadcur(i)= cmplx(fnorm(mnum)) * (
3820 $cmplx(vec1(i)*coef(1,mnum))*form1(mnum,xmaa**2,xmro1**2,xmro2**2)+
3821 $cmplx(vec2(i)*coef(2,mnum))*form2(mnum,xmaa**2,xmro2**2,xmro1**2)+
3822 $cmplx(vec3(i)*coef(3,mnum))*form3(mnum,xmaa**2,xmro3**2,xmro1**2)+
3824 $cmplx(vec4(i)*coef(4,mnum))*form4(mnum,xmaa**2,xmro1**2,
3825 $ xmro2**2,xmro3**2) +
3826 $(-1.0)*uroj/4.0/pi**2/fpi**2*
3827 $cmplx(vec5(i)*coef(5,mnum))*form5(mnum,xmaa**2,xmro1**2,xmro2**2))
3832 CALL clvec(hadcur,pn,pivec)
3833 CALL claxi(hadcur,pn,piaks)
3834 CALL clnut(hadcur,brakm,hvm)
3836 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3837 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3838 amplit=(gfermi)**2*brak/2.
3840 print *,
'MNUM=',mnum
3846 IF (k.EQ.4) znak=1.0
3847 xm1=znak*pim1(k)**2+xm1
3848 xm2=znak*pim2(k)**2+xm2
3849 xm3=znak*pim3(k)**2+xm3
3850 77 print *,
'PIM1=',pim1(k),
'PIM2=',pim2(k),
'PIM3=',pim3(k)
3851 print *,
'XM1=',sqrt(xm1),
'XM2=',sqrt(xm2),
'XM3=',sqrt(xm3)
3852 print *,
'************************************************'
3856 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3857 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3862 SUBROUTINE prod5(P1,P2,P3,PIA)
3868 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3869 COMMON / idfc / idff
3870 REAL PIA(4),P1(4),P2(4),P3(4)
3871 det2(i,j)=p1(i)*p2(j)-p2(i)*p1(j)
3873 IF (ktom.EQ.1.OR.ktom.EQ.-1)
THEN
3874 sign= idff/abs(idff)
3875 ELSEIF (ktom.EQ.2)
THEN
3876 sign=-idff/abs(idff)
3878 print *,
'STOP IN PROD5: KTOM=',ktom
3884 pia(1)= -p3(3)*det2(2,4)+p3(4)*det2(2,3)+p3(2)*det2(3,4)
3885 pia(2)= -p3(4)*det2(1,3)+p3(3)*det2(1,4)-p3(1)*det2(3,4)
3886 pia(3)= p3(4)*det2(1,2)-p3(2)*det2(1,4)+p3(1)*det2(2,4)
3887 pia(4)= p3(3)*det2(1,2)-p3(2)*det2(1,3)+p3(1)*det2(2,3)
3890 20 pia(i)=pia(i)*sign
3893 SUBROUTINE dexnew(MODE,ISGN,POL,PNU,PAA,PNPI,JNPI)
3910 COMMON / inout / inut,iout
3911 REAL POL(4),HV(4),PAA(4),PNU(4),PNPI(4,9),RN(1)
3917 CALL dadnew( -1,isgn,hv,pnu,paa,pnpi,jdumm)
3924 ELSEIF(mode.EQ. 0)
THEN
3927 IF(iwarm.EQ.0)
GOTO 902
3928 CALL dadnew( 0,isgn,hv,pnu,paa,pnpi,jnpi)
3929 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
3932 IF(rn(1).GT.wt)
GOTO 300
3934 ELSEIF(mode.EQ. 1)
THEN
3936 CALL dadnew( 1,isgn,hv,pnu,paa,pnpi,jdumm)
3941 902
WRITE(iout, 9020)
3942 9020
FORMAT(
' ----- DEXNEW: LACK OF INITIALISATION')
3945 SUBROUTINE dadnew(MODE,ISGN,HV,PNU,PWB,PNPI,JNPI)
3947 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3948 * ,ampiz,ampi,amro,gamro,ama1,gama1
3949 * ,amk,amkz,amkst,gamkst
3951 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3952 * ,ampiz,ampi,amro,gamro,ama1,gama1
3953 * ,amk,amkz,amkst,gamkst
3954 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3955 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3956 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
3957 real*4 gampmc ,gamper
3960 COMMON / inout / inut,iout
3962 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
3964 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3966 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3969 CHARACTER NAMES(NMODE)*31
3971 COMMON / inout / inut,iout
3974 real*4 pnu(4),pwb(4),pnpi(4,9),hv(4),hhv(4)
3975 real*4 pdum1(4),pdum2(4),pdumi(4,9)
3978 real*8 swt(nmode),sswt(nmode)
3979 dimension nevraw(nmode),nevovr(nmode),nevacc(nmode)
3981 DATA pi /3.141592653589793238462643/
3997 #if defined (CePeCe)
3999 #elif defined (ALEPH)
4005 IF (jnpi.LE.nm4)
THEN
4006 a = pkorb(3,37+jnpi)
4018 ELSEIF(jnpi.LE.nm4)
THEN
4019 CALL dph4pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4020 ELSEIF(jnpi.LE.nm4+nm5)
THEN
4021 CALL dph5pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4022 ELSEIF(jnpi.LE.nm4+nm5+nm6)
THEN
4023 CALL dphnpi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4024 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3)
THEN
4025 inum=jnpi-nm4-nm5-nm6
4026 CALL dphspk(wt,hv,pdum1,pdum2,pdumi,inum)
4027 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2)
THEN
4028 inum=jnpi-nm4-nm5-nm6-nm3
4029 CALL dphsrk(wt,hv,pdum1,pdum2,pdumi,inum)
4033 IF(wt.GT.wtmax(jnpi)/1.2) wtmax(jnpi)=wt*1.2
4035 #if defined (CePeCe)
4036 #elif defined (ALEPH)
4045 ELSEIF(mode.EQ. 0)
THEN
4047 IF(iwarm.EQ.0)
GOTO 902
4052 ELSEIF(jnpi.LE.nm4)
THEN
4053 CALL dph4pi(wt,hhv,pnu,pwb,pnpi,jnpi)
4054 ELSEIF(jnpi.LE.nm4+nm5)
THEN
4055 CALL dph5pi(wt,hhv,pnu,pwb,pnpi,jnpi)
4056 ELSEIF(jnpi.LE.nm4+nm5+nm6)
THEN
4057 CALL dphnpi(wt,hhv,pnu,pwb,pnpi,jnpi)
4058 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3)
THEN
4059 inum=jnpi-nm4-nm5-nm6
4060 CALL dphspk(wt,hhv,pnu,pwb,pnpi,inum)
4061 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2)
THEN
4062 inum=jnpi-nm4-nm5-nm6-nm3
4063 CALL dphsrk(wt,hhv,pnu,pwb,pnpi,inum)
4071 nevraw(jnpi)=nevraw(jnpi)+1
4072 swt(jnpi)=swt(jnpi)+wt
4074 sswt(jnpi)=sswt(jnpi)+wt**2
4078 sswt(jnpi)=sswt(jnpi)+dble(wt)**2
4083 IF(wt.GT.wtmax(jnpi)) nevovr(jnpi)=nevovr(jnpi)+1
4084 IF(rn*wtmax(jnpi).GT.wt)
GOTO 300
4086 costhe=-1.+2.*rrr(2)
4089 CALL rotor2(thet,pnu,pnu)
4090 CALL rotor3( phi,pnu,pnu)
4091 CALL rotor2(thet,pwb,pwb)
4092 CALL rotor3( phi,pwb,pwb)
4093 CALL rotor2(thet,hv,hv)
4094 CALL rotor3( phi,hv,hv)
4097 CALL rotor2(thet,pnpi(1,i),pnpi(1,i))
4098 CALL rotor3( phi,pnpi(1,i),pnpi(1,i))
4100 nevacc(jnpi)=nevacc(jnpi)+1
4102 ELSEIF(mode.EQ. 1)
THEN
4105 IF(nevraw(jnpi).EQ.0)
GOTO 500
4106 pargam=swt(jnpi)/float(nevraw(jnpi)+1)
4108 IF(nevraw(jnpi).NE.0)
4109 & error=sqrt(sswt(jnpi)/swt(jnpi)**2-1./float(nevraw(jnpi)))
4111 WRITE(iout, 7010) names(jnpi),
4112 & nevraw(jnpi),nevacc(jnpi),nevovr(jnpi),pargam,rat,error
4114 gampmc(8+jnpi-1)=rat
4115 gamper(8+jnpi-1)=error
4121 7003
FORMAT(///1x,15(5h*****)
4122 $ /,
' *', 25x,
'******** DADNEW INITIALISATION ********',9x,1h*
4124 7004
FORMAT(
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*/)
4126 $ /,1x,15(5h*****)/)
4127 7010
FORMAT(///1x,15(5h*****)
4128 $ /,
' *', 25x,
'******** DADNEW FINAL REPORT ******** ',9x,1h*
4129 $ /,
' *', 25x,
'CHANNEL:',a31 ,9x,1h*
4130 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF DECAYS TOTAL ',9x,1h*
4131 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF DECAYS ACCEPTED ',9x,1h*
4132 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
4133 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH IN GEV UNITS ',9x,1h*
4134 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
4135 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
4136 $ /,1x,15(5h*****)/)
4137 902
WRITE(iout, 9020)
4138 9020
FORMAT(
' ----- DADNEW: LACK OF INITIALISATION')
4140 903
WRITE(iout, 9030) jnpi,mode
4141 9030
FORMAT(
' ----- DADNEW: WRONG JNPI',2i5)
4146 SUBROUTINE dph4pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4156 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4157 * ,ampiz,ampi,amro,gamro,ama1,gama1
4158 * ,amk,amkz,amkst,gamkst
4160 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4161 * ,ampiz,ampi,amro,gamro,ama1,gama1
4162 * ,amk,amkz,amkst,gamkst
4163 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4164 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4166 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4167 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4169 CHARACTER NAMES(NMODE)*31
4172 REAL HV(4),PT(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4),PMULT(4,9)
4175 real*8 uu,ff,ff1,ff2,ff3,ff4,gg1,gg2,gg3,gg4,rr
4176 DATA pi /3.141592653589793238462643/
4178 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
4183 phspac=1./2**23/pi**11
4187 amp1=dcdmas(idffin(1,jnpi))
4188 amp2=dcdmas(idffin(2,jnpi))
4189 amp3=dcdmas(idffin(3,jnpi))
4190 amp4=dcdmas(idffin(4,jnpi))
4194 #if defined (CePeCe)
4201 #elif defined (CLEO)
4236 CALL choice(100+jnpi,rrdum,ichan,prob1,prob2,prob3,
4239 CALL choice(100+jnpi,rrb,ichan,prob1,prob2,prob3,
4241 $ amrop,gamrop,amrx,gamrx,amrb,gamrb)
4255 ams1=(amp1+amp2+amp3+amp4)**2
4256 ams2=(amtau-amnuta)**2
4257 alp1=atan((ams1-amrop**2)/amrop/gamrop)
4258 alp2=atan((ams2-amrop**2)/amrop/gamrop)
4259 alp=alp1+rr1*(alp2-alp1)
4260 am4sq =amrop**2+amrop*gamrop*tan(alp)
4263 $ ((am4sq-amrop**2)**2+(amrop*gamrop)**2)/(amrop*gamrop)
4264 phspac=phspac*(alp2-alp1)
4268 ams1=(amp2+amp3+amp4)**2
4270 IF (rrr(9).GT.prez)
THEN
4271 am3sq=ams1+ rr1*(ams2-ams1)
4277 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4278 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4279 alp=alp1+rr1*(alp2-alp1)
4280 am3sq =amrx**2+amrx*gamrx*tan(alp)
4283 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4291 am2sq=ams1+ rr2*(ams2-ams1)
4297 enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4298 enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4299 ppi= enq1**2-amp3**2
4300 pppi=sqrt(abs(enq1**2-amp3**2))
4302 enq1=(am2sq-amp3**2+amp4**2)/(2*am2)
4303 enq2=(am2sq+amp3**2-amp4**2)/(2*am2)
4304 ppi= enq1**2-amp4**2
4305 pppi=sqrt(abs(enq1**2-amp4**2))
4307 phspac=phspac*(4*pi)*(2*pppi/am2)
4313 CALL sphera(pppi,piz)
4331 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4332 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4333 ppi = pr(4)**2-am2**2
4341 pim1(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4344 ff3=(4*pi)*(2*pr(3)/am3)
4346 exe=(pr(4)+pr(3))/am2
4347 CALL bostr3(exe,piz,piz)
4348 CALL bostr3(exe,pipl,pipl)
4351 thet =acos(-1.+2*rr3)
4353 CALL rotpol(thet,phi,pipl)
4354 CALL rotpol(thet,phi,pim1)
4355 CALL rotpol(thet,phi,piz)
4356 CALL rotpol(thet,phi,pr)
4366 pr(4)=1./(2*am4)*(am4**2+am3**2-amp1**2)
4367 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4368 ppi = pr(4)**2-am3**2
4376 pim2(4)=1./(2*am4)*(am4**2-am3**2+amp1**2)
4379 ff4=(4*pi)*(2*pr(3)/am4)
4381 exe=(pr(4)+pr(3))/am3
4382 CALL bostr3(exe,piz,piz)
4383 CALL bostr3(exe,pipl,pipl)
4384 CALL bostr3(exe,pim1,pim1)
4387 thet =acos(-1.+2*rr3)
4389 CALL rotpol(thet,phi,pipl)
4390 CALL rotpol(thet,phi,pim1)
4391 CALL rotpol(thet,phi,pim2)
4392 CALL rotpol(thet,phi,piz)
4393 CALL rotpol(thet,phi,pr)
4399 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am4**2)
4400 paa(3)= sqrt(abs(paa(4)**2-am4**2))
4401 ppi = paa(4)**2-am4**2
4402 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4403 phsp=phsp*(4*pi)*(2*paa(3)/amtau)
4407 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am4**2)
4411 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4412 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4414 ams1=(amp1+amp3+amp4)**2
4417 ams2=(sqrt(am3sq)-amp1)**2
4419 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4420 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4423 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4424 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4426 ams1=(amp1+amp3+amp4)**2
4427 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4428 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4429 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4432 ams2=(sqrt(am3sq)-amp1)**2
4434 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4435 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4438 am3sq=(pim2(4)+piz(4)+pipl(4))**2-(pim2(3)+piz(3)+pipl(3))**2
4439 $ -(pim2(2)+piz(2)+pipl(2))**2-(pim2(1)+piz(1)+pipl(1))**2
4441 ams1=(amp2+amp3+amp4)**2
4442 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4443 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4444 gg1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4447 ams2=(sqrt(am3sq)-amp2)**2
4449 gg3=(4*pi)*(xlam(am2**2,amp2**2,am3sq)/am3sq)
4450 gg4=(4*pi)*(xlam(am3sq,amp1**2,am4**2)/am4**2)
4453 IF ( ( (ff+gg)*uu+ff*gg ).GT.0.0d0)
THEN
4454 rr=ff*gg*uu/(0.5*prez*(ff+gg)*uu+(1.0-prez)*ff*gg)
4482 exe=(paa(4)+paa(3))/am4
4483 CALL bostr3(exe,piz,piz)
4484 CALL bostr3(exe,pipl,pipl)
4485 CALL bostr3(exe,pim1,pim1)
4486 CALL bostr3(exe,pim2,pim2)
4487 CALL bostr3(exe,pr,pr)
4492 CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4499 CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4500 ELSEIF (jnpi.EQ.2)
THEN
4501 CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv)
4504 dgamt=1/(2.*amtau)*amplit*phspac
4519 SUBROUTINE dam4pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,AMPLIT,HV)
4533 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4534 * ,ampiz,ampi,amro,gamro,ama1,gama1
4535 * ,amk,amkz,amkst,gamkst
4537 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4538 * ,ampiz,ampi,amro,gamro,ama1,gama1
4539 * ,amk,amkz,amkst,gamkst
4540 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4541 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4542 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4),PIM4(4)
4543 REAL PIVEC(4),PIAKS(4),HVM(4)
4544 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5
4545 EXTERNAL form1,form2,form3,form4,form5
4546 DATA pi /3.141592653589793238462643/
4550 CALL curr_cleo(mnum,pim1,pim2,pim3,pim4,hadcur)
4552 CALL curr(mnum,pim1,pim2,pim3,pim4,hadcur)
4556 CALL clvec(hadcur,pn,pivec)
4557 CALL claxi(hadcur,pn,piaks)
4558 CALL clnut(hadcur,brakm,hvm)
4560 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4561 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4562 amplit=(ccabib*gfermi)**2*brak/2.
4565 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4566 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4572 SUBROUTINE dph5pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4577 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4578 * ,ampiz,ampi,amro,gamro,ama1,gama1
4579 * ,amk,amkz,amkst,gamkst
4581 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4582 * ,ampiz,ampi,amro,gamro,ama1,gama1
4585 * ,amk,amkz,amkst,gamkst
4586 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4587 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4588 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4590 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4592 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4595 CHARACTER NAMES(NMODE)*31
4596 REAL HV(4),PT(4),PN(4),PAA(4),PMULT(4,9)
4597 real*4 pr(4),pi1(4),pi2(4),pi3(4),pi4(4),pi5(4)
4598 real*8 amp1,amp2,amp3,amp4,amp5,ams1,ams2,amom,gamom
4599 real*8 am5sq,am4sq,am3sq,am2sq,am5,am4,am3
4601 real*8 gg1,gg2,gg3,ff1,ff2,ff3,ff4,alp,alp1,alp2
4610 DATA pi /3.141592653589793238462643/
4617 bwign(xm,am,gammab)=xm**2/cmplx(xm**2-am**2,gammab*am)
4619 bwign(xm,am,gamma)=xm**2/cmplx(xm**2-am**2,gamma*am)
4628 phspac=1./2**29/pi**14
4631 amp1=dcdmas(idffin(1,jnpi))
4632 amp2=dcdmas(idffin(2,jnpi))
4633 amp3=dcdmas(idffin(3,jnpi))
4634 amp4=dcdmas(idffin(4,jnpi))
4635 amp5=dcdmas(idffin(5,jnpi))
4650 ams1=(amp1+amp2+amp3+amp4+amp5)**2
4651 ams2=(amtau-amnuta)**2
4652 am5sq=ams1+ rr1*(ams2-ams1)
4654 phspac=phspac*(ams2-ams1)
4659 ams1=(amp2+amp3+amp4+amp5)**2
4661 am4sq=ams1+ rr1*(ams2-ams1)
4668 ams1=(amp2+amp3+amp4)**2
4670 alp1=atan((ams1-amom**2)/amom/gamom)
4671 alp2=atan((ams2-amom**2)/amom/gamom)
4672 alp=alp1+rr1*(alp2-alp1)
4673 am3sq =amom**2+amom*gamom*tan(alp)
4676 gg2=((am3sq-amom**2)**2+(amom*gamom)**2)/(amom*gamom)
4689 am2sq=ams1+ rr2*(ams2-ams1)
4695 enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4696 enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4697 ppi= enq1**2-amp3**2
4698 pppi=sqrt(abs(enq1**2-amp3**2))
4699 ff1=(4*pi)*(2*pppi/am2)
4701 call sphera(pppi,pi3)
4712 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4713 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4714 ppi = pr(4)**2-am2**2
4718 pi2(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4721 ff2=(4*pi)*(2*pr(3)/am3)
4723 exe=(pr(4)+pr(3))/am2
4724 call bostr3(exe,pi3,pi3)
4725 call bostr3(exe,pi4,pi4)
4728 thet =acos(-1.+2*rr3)
4730 call rotpol(thet,phi,pi2)
4731 call rotpol(thet,phi,pi3)
4732 call rotpol(thet,phi,pi4)
4738 pr(4)=1./(2*am4)*(am4**2+am3**2-amp5**2)
4739 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4740 ppi = pr(4)**2-am3**2
4744 pi5(4)=1./(2*am4)*(am4**2-am3**2+amp5**2)
4747 ff3=(4*pi)*(2*pr(3)/am4)
4749 exe=(pr(4)+pr(3))/am3
4750 call bostr3(exe,pi2,pi2)
4751 call bostr3(exe,pi3,pi3)
4752 call bostr3(exe,pi4,pi4)
4755 thet =acos(-1.+2*rr3)
4757 call rotpol(thet,phi,pi2)
4758 call rotpol(thet,phi,pi3)
4759 call rotpol(thet,phi,pi4)
4760 call rotpol(thet,phi,pi5)
4766 pr(4)=1./(2*am5)*(am5**2+am4**2-amp1**2)
4767 pr(3)= sqrt(abs(pr(4)**2-am4**2))
4768 ppi = pr(4)**2-am4**2
4772 pi1(4)=1./(2*am5)*(am5**2-am4**2+amp1**2)
4775 ff4=(4*pi)*(2*pr(3)/am5)
4777 exe=(pr(4)+pr(3))/am4
4778 call bostr3(exe,pi2,pi2)
4779 call bostr3(exe,pi3,pi3)
4780 call bostr3(exe,pi4,pi4)
4781 call bostr3(exe,pi5,pi5)
4784 thet =acos(-1.+2*rr3)
4786 call rotpol(thet,phi,pi1)
4787 call rotpol(thet,phi,pi2)
4788 call rotpol(thet,phi,pi3)
4789 call rotpol(thet,phi,pi4)
4790 call rotpol(thet,phi,pi5)
4799 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5sq)
4800 paa(3)= sqrt(abs(paa(4)**2-am5sq))
4801 ppi = paa(4)**2-am5sq
4802 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4806 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am5**2)
4809 phspac=phspac * gg1*gg2*gg3*ff1*ff2*ff3*ff4
4813 exe=(paa(4)+paa(3))/am5
4814 call bostr3(exe,pi1,pi1)
4815 call bostr3(exe,pi2,pi2)
4816 call bostr3(exe,pi3,pi3)
4817 call bostr3(exe,pi4,pi4)
4818 call bostr3(exe,pi5,pi5)
4826 qxn=paa(4)*pn(4)-paa(1)*pn(1)-paa(2)*pn(2)-paa(3)*pn(3)
4827 brak=2*(gv**2+ga**2)*(2*pxq*qxn+am5sq*pxn)
4828 & -6*(gv**2-ga**2)*amtau*amnuta*am5sq
4829 fompp = cabs(bwign(am3,amom,gamom))**2
4834 amplit=ccabib**2*gfermi**2/2. * brak
4835 amplit = amplit * fompp * fnorm
4838 dgamt=1/(2.*amtau)*amplit*phspac
4861 SUBROUTINE dphsrk(DGAMT,HV,PN,PR,PMULT,INUM)
4867 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4868 * ,ampiz,ampi,amro,gamro,ama1,gama1
4869 * ,amk,amkz,amkst,gamkst
4871 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4872 * ,ampiz,ampi,amro,gamro,ama1,gama1
4873 * ,amk,amkz,amkst,gamkst
4874 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4875 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4876 REAL HV(4),PT(4),PN(4),PR(4),PKC(4),PKZ(4),QQ(4),PMULT(4,9)
4882 DATA pi /3.141592653589793238462643/
4886 phspac=1./2**11/pi**5
4894 ams2=(amtau-amnuta)**2
4897 amx2=ams1+ rr1(1)*(ams2-ams1)
4899 phspac=phspac*(ams2-ams1)
4917 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
4918 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
4922 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
4924 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
4927 enq1=(amx2+amk**2-amkz**2)/(2.*amx)
4928 enq2=(amx2-amk**2+amkz**2)/(2.*amx)
4929 pppi=sqrt((enq1-amk)*(enq1+amk))
4930 phspac=phspac*(4*pi)*(2*pppi/amx)
4932 CALL sphera(pppi,pkc)
4938 exe=(pr(4)+pr(3))/amx
4940 CALL bostr3(exe,pkc,pkc)
4941 CALL bostr3(exe,pkz,pkz)
4943 30 qq(i)=pkc(i)-pkz(i)
4945 pksd =pr(4)*pr(4)-pr(3)*pr(3)-pr(2)*pr(2)-pr(1)*pr(1)
4946 qqpks=pr(4)* qq(4)-pr(3)* qq(3)-pr(2)* qq(2)-pr(1)* qq(1)
4948 31 qq(i)=qq(i)-pr(i)*qqpks/pksd
4951 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
4953 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
4954 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
4955 & +(gv**2-ga**2)*amtau*amnuta*qq2
4956 amplit=(gfermi*ccabib)**2*brak*2*fpirk(amx)
4957 dgamt=1/(2.*amtau)*amplit*phspac
4959 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
4970 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4971 * ,ampiz,ampi,amro,gamro,ama1,gama1
4972 * ,amk,amkz,amkst,gamkst
4974 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4975 * ,ampiz,ampi,amro,gamro,ama1,gama1
4976 * ,amk,amkz,amkst,gamkst
4979 fpirk=cabs(fpikm(w,amk,amkz))**2
4982 COMPLEX FUNCTION fpikmk(W,XM1,XM2)
4987 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
4992 IF (init.EQ.0 )
THEN
5005 fpikmk=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
5014 parameter(nmxhep=2000)
5015 common/hepevtx/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
5016 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
5021 SUBROUTINE dwrph(KTO,PHX)
5025 IMPLICIT REAL*8 (a-h,o-z)
5036 IF (qhot(4).GT.1.e-5)
CALL dwluph(kto,qhot)
5039 SUBROUTINE dwluph(KTO,PHOT)
5050 COMMON /taupos/ np1,np2
5056 COMMON /taupos/ np1,np2
5060 IF (phot(4).LE.0.0)
RETURN
5063 IF((kto.EQ. 1).OR.(kto.EQ.11))
THEN
5070 IF(ktos.GT.10) ktos=ktos-10
5072 CALL tralo4(ktos,phot,phot,am)
5073 CALL filhep(0,1,22,nps,nps,0,0,phot,0.0,.true.)
5078 SUBROUTINE dwluel(KTO,ISGN,PNU,PWB,PEL,PNE)
5089 COMMON /taupos/ np1,np2
5092 REAL PNU(4),PWB(4),PEL(4),PNE(4)
5095 COMMON /taupos/ np1,np2
5106 CALL tralo4(kto,pnu,pnu,am)
5107 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5110 CALL tralo4(kto,pwb,pwb,am)
5114 CALL tralo4(kto,pel,pel,am)
5115 CALL filhep(0,1,11*isgn,nps,nps,0,0,pel,am,.false.)
5118 CALL tralo4(kto,pne,pne,am)
5119 CALL filhep(0,1,-12*isgn,nps,nps,0,0,pne,am,.true.)
5123 SUBROUTINE dwlumu(KTO,ISGN,PNU,PWB,PMU,PNM)
5134 COMMON /taupos/ np1,np2
5137 REAL PNU(4),PWB(4),PMU(4),PNM(4)
5140 COMMON /taupos/ np1,np2
5151 CALL tralo4(kto,pnu,pnu,am)
5152 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5155 CALL tralo4(kto,pwb,pwb,am)
5159 CALL tralo4(kto,pmu,pmu,am)
5160 CALL filhep(0,1,13*isgn,nps,nps,0,0,pmu,am,.false.)
5163 CALL tralo4(kto,pnm,pnm,am)
5164 CALL filhep(0,1,-14*isgn,nps,nps,0,0,pnm,am,.true.)
5168 SUBROUTINE dwlupi(KTO,ISGN,PPI,PNU)
5179 COMMON /taupos/ np1,np2
5189 CALL tralo4(kto,pnu,pnu,am)
5190 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5193 CALL tralo4(kto,ppi,ppi,am)
5194 CALL filhep(0,1,-211*isgn,nps,nps,0,0,ppi,am,.true.)
5198 SUBROUTINE dwluro(KTO,ISGN,PNU,PRHO,PIC,PIZ)
5209 COMMON /taupos/ np1,np2
5212 REAL PNU(4),PRHO(4),PIC(4),PIZ(4)
5215 COMMON /taupos/ np1,np2
5226 CALL tralo4(kto,pnu,pnu,am)
5227 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5230 CALL tralo4(kto,prho,prho,am)
5231 CALL filhep(0,2,-213*isgn,nps,nps,0,0,prho,am,.true.)
5234 CALL tralo4(kto,pic,pic,am)
5235 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pic,am,.true.)
5238 CALL tralo4(kto,piz,piz,am)
5239 CALL filhep(0,1,111,-2,-2,0,0,piz,am,.true.)
5243 SUBROUTINE dwluaa(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
5255 COMMON /taupos/ np1,np2
5258 REAL PNU(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
5261 COMMON /taupos/ np1,np2
5272 CALL tralo4(kto,pnu,pnu,am)
5273 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5276 CALL tralo4(kto,paa,paa,am)
5277 CALL filhep(0,1,-20213*isgn,nps,nps,0,0,paa,am,.true.)
5285 CALL tralo4(kto,pim2,pim2,am)
5286 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pim2,am,.true.)
5289 CALL tralo4(kto,pim1,pim1,am)
5290 CALL filhep(0,1,-211*isgn,-2,-2,0,0,pim1,am,.true.)
5293 CALL tralo4(kto,pipl,pipl,am)
5294 CALL filhep(0,1, 211*isgn,-3,-3,0,0,pipl,am,.true.)
5296 ELSE IF (jaa.EQ.2)
THEN
5301 CALL tralo4(kto,pim2,pim2,am)
5302 CALL filhep(0,1,111,-1,-1,0,0,pim2,am,.true.)
5305 CALL tralo4(kto,pim1,pim1,am)
5306 CALL filhep(0,1,111,-2,-2,0,0,pim1,am,.true.)
5309 CALL tralo4(kto,pipl,pipl,am)
5310 CALL filhep(0,1,-211*isgn,-3,-3,0,0,pipl,am,.true.)
5316 SUBROUTINE dwlukk (KTO,ISGN,PKK,PNU)
5326 COMMON /taupos/ np1,np2
5340 CALL tralo4 (kto,pnu,pnu,am)
5341 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5344 CALL tralo4 (kto,pkk,pkk,am)
5345 CALL filhep(0,1,-321*isgn,nps,nps,0,0,pkk,am,.true.)
5349 SUBROUTINE dwluks(KTO,ISGN,PNU,PKS,PKK,PPI,JKST)
5350 COMMON / taukle / bra1,brk0,brk0b,brks
5351 real*4 bra1,brk0,brk0b,brks
5353 COMMON /taupos/ np1,np2
5366 REAL PNU(4),PKS(4),PKK(4),PPI(4)
5368 REAL PNU(4),PKS(4),PKK(4),PPI(4),XIO(1)
5369 COMMON /taupos/ np1,np2
5380 CALL tralo4(kto,pnu,pnu,am)
5381 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5384 CALL tralo4(kto,pks,pks,am)
5385 CALL filhep(0,1,-323*isgn,nps,nps,0,0,pks,am,.true.)
5393 CALL tralo4(kto,ppi,ppi,am)
5394 CALL filhep(0,1,-211*isgn,-1,-1,0,0,ppi,am,.true.)
5397 IF (isgn.EQ.-1) bran=brk0
5400 IF(xio(1).GT.bran)
THEN
5406 CALL tralo4(kto,pkk,pkk,am)
5407 CALL filhep(0,1,k0type,-2,-2,0,0,pkk,am,.true.)
5409 ELSE IF(jkst.EQ.20)
THEN
5414 CALL tralo4(kto,ppi,ppi,am)
5415 CALL filhep(0,1,111,-1,-1,0,0,ppi,am,.true.)
5418 CALL tralo4(kto,pkk,pkk,am)
5419 CALL filhep(0,1,-321*isgn,-2,-2,0,0,pkk,am,.true.)
5425 SUBROUTINE dwlnew(KTO,ISGN,PNU,PWB,PNPI,MODE)
5435 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
5437 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5439 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5442 COMMON /taupos/ np1,np2
5443 CHARACTER NAMES(NMODE)*31
5444 REAL PNU(4),PWB(4),PNPI(4,9)
5456 CALL tralo4(kto,pnu,pnu,am)
5457 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5460 CALL tralo4(kto,pwb,pwb,am)
5461 CALL filhep(0,1,-24*isgn,nps,nps,0,0,pwb,am,.true.)
5470 kfpi=lunpik(idffin(i,jnpi), isgn)
5472 kfpi=lunpik(idffin(i,jnpi),-isgn)
5479 CALL tralo4(kto,ppi,ppi,am)
5480 CALL filhep(0,1,kfpi,-i,-i,0,0,ppi,am,.true.)
5485 #if defined (CePeCe)
5494 IMPLICIT REAL*8 (a-h,o-z)
5496 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5498 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5510 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5511 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5520 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5521 DATA pi /3.141592653589793238462643d0/
5523 IF(abs(y).LT.abs(x))
THEN
5525 IF(x.LE.0d0) the=pi-the
5527 the=acos(x/sqrt(x**2+y**2))
5538 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5539 DATA pi /3.141592653589793238462643d0/
5541 IF(abs(y).LT.abs(x))
THEN
5543 IF(x.LE.0d0) the=pi-the
5545 the=acos(x/sqrt(x**2+y**2))
5547 IF(y.LT.0d0) the=2d0*pi-the
5550 SUBROUTINE rotod1(PH1,PVEC,QVEC)
5555 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5556 dimension pvec(4),qvec(4),rvec(4)
5564 qvec(2)= cs*rvec(2)-sn*rvec(3)
5565 qvec(3)= sn*rvec(2)+cs*rvec(3)
5569 SUBROUTINE rotod2(PH1,PVEC,QVEC)
5574 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5575 dimension pvec(4),qvec(4),rvec(4)
5582 qvec(1)= cs*rvec(1)+sn*rvec(3)
5584 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5588 SUBROUTINE rotod3(PH1,PVEC,QVEC)
5593 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5595 dimension pvec(4),qvec(4),rvec(4)
5601 qvec(1)= cs*rvec(1)-sn*rvec(2)
5602 qvec(2)= sn*rvec(1)+cs*rvec(2)
5606 SUBROUTINE bostr3(EXE,PVEC,QVEC)
5612 real*4 pvec(4),qvec(4),rvec(4)
5625 SUBROUTINE bostd3(EXE,PVEC,QVEC)
5631 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5632 dimension pvec(4),qvec(4),rvec(4)
5646 SUBROUTINE rotor1(PH1,PVEC,QVEC)
5651 real*4 pvec(4),qvec(4),rvec(4)
5659 qvec(2)= cs*rvec(2)-sn*rvec(3)
5660 qvec(3)= sn*rvec(2)+cs*rvec(3)
5663 SUBROUTINE rotor2(PH1,PVEC,QVEC)
5668 IMPLICIT REAL*4(a-h,o-z)
5669 real*4 pvec(4),qvec(4),rvec(4)
5676 qvec(1)= cs*rvec(1)+sn*rvec(3)
5678 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5681 SUBROUTINE rotor3(PHI,PVEC,QVEC)
5686 real*4 pvec(4),qvec(4),rvec(4)
5692 qvec(1)= cs*rvec(1)-sn*rvec(2)
5693 qvec(2)= sn*rvec(1)+cs*rvec(2)
5697 SUBROUTINE spherd(R,X)
5702 real*8 r,x(4),pi,costh,sinth
5704 DATA pi /3.141592653589793238462643d0/
5708 sinth=sqrt(1 -costh**2)
5709 x(1)=r*sinth*cos(2*pi*rrr(2))
5710 x(2)=r*sinth*sin(2*pi*rrr(2))
5714 SUBROUTINE rotpox(THET,PHI,PP)
5715 IMPLICIT REAL*8 (a-h,o-z)
5725 CALL rotod2(thet,pp,pp)
5726 CALL rotod3( phi,pp,pp)
5729 SUBROUTINE sphera(R,X)
5737 DATA pi /3.141592653589793238462643/
5741 sinth=sqrt(1.-costh**2)
5742 x(1)=r*sinth*cos(2*pi*rrr(2))
5743 x(2)=r*sinth*sin(2*pi*rrr(2))
5747 SUBROUTINE rotpol(THET,PHI,PP)
5754 CALL rotor2(thet,pp,pp)
5755 CALL rotor3( phi,pp,pp)
5758 #include "../randg/tauola-random.h"
5761 IMPLICIT REAL*8(a-h,o-z)
5764 IF(x .LT.-1.0)
GO TO 1
5765 IF(x .LE. 0.5)
GO TO 2
5766 IF(x .EQ. 1.0)
GO TO 3
5767 IF(x .LE. 2.0)
GO TO 4
5771 z=z-0.5* log(abs(x))**2
5777 3 dilogt=1.64493406684822
5781 z=1.64493406684822 - log(x)* log(abs(t))
5782 5 y=2.66666666666666 *t+0.66666666666666
5783 b= 0.00000 00000 00001
5784 a=y*b +0.00000 00000 00004
5785 b=y*a-b+0.00000 00000 00011
5786 a=y*b-a+0.00000 00000 00037
5787 b=y*a-b+0.00000 00000 00121
5788 a=y*b-a+0.00000 00000 00398
5789 b=y*a-b+0.00000 00000 01312
5790 a=y*b-a+0.00000 00000 04342
5791 b=y*a-b+0.00000 00000 14437
5792 a=y*b-a+0.00000 00000 48274
5793 b=y*a-b+0.00000 00001 62421
5794 a=y*b-a+0.00000 00005 50291
5795 b=y*a-b+0.00000 00018 79117
5796 a=y*b-a+0.00000 00064 74338
5797 b=y*a-b+0.00000 00225 36705
5798 a=y*b-a+0.00000 00793 87055
5799 b=y*a-b+0.00000 02835 75385
5800 a=y*b-a+0.00000 10299 04264
5801 b=y*a-b+0.00000 38163 29463
5802 a=y*b-a+0.00001 44963 00557
5803 b=y*a-b+0.00005 68178 22718
5804 a=y*b-a+0.00023 20021 96094
5805 b=y*a-b+0.00100 16274 96164
5806 a=y*b-a+0.00468 63619 59447
5807 b=y*a-b+0.02487 93229 24228
5808 a=y*b-a+0.16607 30329 27855
5809 a=y*a-b+1.93506 43008 6996