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.7 ******',9x,1h*,
217 $ /,
' *', 25x,
'*DEC 1993; ALEPH fixes introd. dec 98 *',9x,1h*,
218 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
219 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
220 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
221 $ /,
' *', 25x,
'Physics initialization by ALEPH collab ',9x,1h*,
222 $ /,
' *', 25x,
'it is suggested to use this version ',9x,1h*,
223 $ /,
' *', 25x,
' with the help of the collab. advice ',9x,1h*,
224 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
225 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
227 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
228 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
229 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
230 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
231 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
232 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
233 $ /,
' *', 25x,
' Physics initialization by CLEO collab ',9x,1h*,
234 $ /,
' *', 25x,
' see Alain Weinstein www home page: ',9x,1h*,
235 $ /,
' *', 25x,
'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
236 $ /,
' *', 25x,
'/korb_doc.html#files ',9x,1h*,
237 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
238 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
240 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
241 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
242 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
243 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
244 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
245 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
246 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
247 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
248 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
249 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
251 $ /,
' *', 25x,
'****DEKAY ROUTINE: INITIALIZATION******',9x,1h*,
252 $ /,
' *',i20 ,5x,
'JAK1 = DECAY MODE TAU+ ',9x,1h*,
253 $ /,
' *',i20 ,5x,
'JAK2 = DECAY MODE TAU- ',9x,1h*,
255 7010
FORMAT(///1x,15(5h*****)
257 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
258 $ /,
' *', 25x,
'***********DECEMBER 1993***************',9x,1h*,
260 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
261 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
263 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
264 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
265 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
266 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
267 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
268 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
269 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
270 $ /,
' *', 25x,
'*****DEKAY ROUTINE: FINAL REPORT*******',9x,1h*,
271 $ /,
' *',i20 ,5x,
'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*,
272 $ /,
' *',i20 ,5x,
'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*,
273 $ /,
' *',i20 ,5x,
'NEVTOT = SUM ',9x,1h*,
275 $
' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
277 $ ,i10,2f12.7 ,
' DADMEL ELECTRON ',9x,1h*
278 $ /,
' *',i10,2f12.7 ,
' DADMMU MUON ',9x,1h*
279 $ /,
' *',i10,2f12.7 ,
' DADMPI PION ',9x,1h*
280 $ /,
' *',i10,2f12.7,
' DADMRO RHO (->2PI) ',9x,1h*
281 $ /,
' *',i10,2f12.7,
' DADMAA A1 (->3PI) ',9x,1h*
282 $ /,
' *',i10,2f12.7,
' DADMKK KAON ',9x,1h*
283 $ /,
' *',i10,2f12.7,
' DADMKS K* ',9x,1h*)
285 $ ,i10,2f12.7,a31 ,8x,1h*)
287 $ ,20x,
'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
288 $ /,
' *',20x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
291 9020
FORMAT(
' ----- DEKAY: LACK OF INITIALISATION')
294 9100
FORMAT(
' ----- DEKAY: WRONG VALUE OF KTO ')
297 SUBROUTINE dekay1(IMOD,HH,ISGN)
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.7 ******',9x,1h*,
611 $ /,
' *', 25x,
'*DEC 1993; ALEPH fixes introd. dec 98 *',9x,1h*,
612 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
613 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
614 $ /,
' *', 25x,
'Physics initialization by ALEPH collab ',9x,1h*,
615 $ /,
' *', 25x,
'it is suggested to use this version ',9x,1h*,
616 $ /,
' *', 25x,
' with the help of the collab. advice ',9x,1h*,
617 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
618 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
620 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
621 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
622 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
623 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
624 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
625 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
626 $ /,
' *', 25x,
' Physics initialization by CLEO collab ',9x,1h*,
627 $ /,
' *', 25x,
' see Alain Weinstein www home page: ',9x,1h*,
628 $ /,
' *', 25x,
'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
629 $ /,
' *', 25x,
'/korb_doc.html#files ',9x,1h*,
630 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
631 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
633 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
634 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
635 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
636 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
637 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
638 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
639 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
640 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
642 $ /,
' *', 25x,
'*******CERN-TH-6793 NOVEMBER 1992*****',9x,1h*,
643 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
644 $ /,
' *', 25x,
'******DEXAY ROUTINE: INITIALIZATION****',9x,1h*
645 $ /,
' *',i20 ,5x,
'JAK1 = DECAY MODE FERMION1 (TAU+) ',9x,1h*
646 $ /,
' *',i20 ,5x,
'JAK2 = DECAY MODE FERMION2 (TAU-) ',9x,1h*
650 7010
FORMAT(///1x,15(5h*****)
652 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
653 $ /,
' *', 25x,
'***********DECEMBER 1993***************',9x,1h*,
655 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.7 ******',9x,1h*,
656 $ /,
' *', 25x,
'***********August 1995***************',9x,1h*,
658 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
659 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
660 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
661 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
662 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
663 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
664 $ /,
' *', 25x,
'*******CERN-TH-6793 NOVEMBER 1992*****',9x,1h*,
665 $ /,
' *', 25x,
'******DEXAY ROUTINE: FINAL REPORT******',9x,1h*
666 $ /,
' *',i20 ,5x,
'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*
667 $ /,
' *',i20 ,5x,
'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*
668 $ /,
' *',i20 ,5x,
'NEVTOT = SUM ',9x,1h*
670 $
' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
672 $ ,i10,2f12.7 ,
' DADMEL ELECTRON ',9x,1h*
673 $ /,
' *',i10,2f12.7 ,
' DADMMU MUON ',9x,1h*
674 $ /,
' *',i10,2f12.7 ,
' DADMPI PION ',9x,1h*
675 $ /,
' *',i10,2f12.7,
' DADMRO RHO (->2PI) ',9x,1h*
676 $ /,
' *',i10,2f12.7,
' DADMAA A1 (->3PI) ',9x,1h*
677 $ /,
' *',i10,2f12.7,
' DADMKK KAON ',9x,1h*
678 $ /,
' *',i10,2f12.7,
' DADMKS K* ',9x,1h*)
680 $ ,i10,2f12.7,a31 ,8x,1h*)
682 $ ,20x,
'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
683 $ /,
' *',20x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
685 902
WRITE(iout, 9020)
686 9020
FORMAT(
' ----- DEXAY: LACK OF INITIALISATION')
688 910
WRITE(iout, 9100)
689 9100
FORMAT(
' ----- DEXAY: WRONG VALUE OF KTO ')
692 SUBROUTINE dexay1(KTO,JAKIN,JAK,POL,ISGN)
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)
4410 IF(rrr(9).LE.0.5*prez)
THEN
4419 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4420 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4422 ams1=(amp1+amp3+amp4)**2
4425 ams2=(sqrt(am3sq)-amp1)**2
4427 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4428 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4431 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4432 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4434 ams1=(amp1+amp3+amp4)**2
4435 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4436 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4437 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4440 ams2=(sqrt(am3sq)-amp1)**2
4442 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4443 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4446 am3sq=(pim2(4)+piz(4)+pipl(4))**2-(pim2(3)+piz(3)+pipl(3))**2
4447 $ -(pim2(2)+piz(2)+pipl(2))**2-(pim2(1)+piz(1)+pipl(1))**2
4449 ams1=(amp2+amp3+amp4)**2
4450 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4451 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4452 gg1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4455 ams2=(sqrt(am3sq)-amp2)**2
4457 gg3=(4*pi)*(xlam(am2**2,amp2**2,am3sq)/am3sq)
4458 gg4=(4*pi)*(xlam(am3sq,amp1**2,am4**2)/am4**2)
4461 IF ( ( (ff+gg)*uu+ff*gg ).GT.0.0d0)
THEN
4462 rr=ff*gg*uu/(0.5*prez*(ff+gg)*uu+(1.0-prez)*ff*gg)
4490 exe=(paa(4)+paa(3))/am4
4491 CALL bostr3(exe,piz,piz)
4492 CALL bostr3(exe,pipl,pipl)
4493 CALL bostr3(exe,pim1,pim1)
4494 CALL bostr3(exe,pim2,pim2)
4495 CALL bostr3(exe,pr,pr)
4500 CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4507 CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4508 ELSEIF (jnpi.EQ.2)
THEN
4509 CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv)
4512 dgamt=1/(2.*amtau)*amplit*phspac
4527 SUBROUTINE dam4pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,AMPLIT,HV)
4541 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4542 * ,ampiz,ampi,amro,gamro,ama1,gama1
4543 * ,amk,amkz,amkst,gamkst
4545 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4546 * ,ampiz,ampi,amro,gamro,ama1,gama1
4547 * ,amk,amkz,amkst,gamkst
4548 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4549 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4550 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4),PIM4(4)
4551 REAL PIVEC(4),PIAKS(4),HVM(4)
4552 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5
4553 EXTERNAL form1,form2,form3,form4,form5
4554 DATA pi /3.141592653589793238462643/
4558 CALL curr_cleo(mnum,pim1,pim2,pim3,pim4,hadcur)
4560 CALL curr(mnum,pim1,pim2,pim3,pim4,hadcur)
4564 CALL clvec(hadcur,pn,pivec)
4565 CALL claxi(hadcur,pn,piaks)
4566 CALL clnut(hadcur,brakm,hvm)
4568 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4569 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4570 amplit=(ccabib*gfermi)**2*brak/2.
4573 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4574 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4580 SUBROUTINE dph5pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4585 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4586 * ,ampiz,ampi,amro,gamro,ama1,gama1
4587 * ,amk,amkz,amkst,gamkst
4589 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4590 * ,ampiz,ampi,amro,gamro,ama1,gama1
4593 * ,amk,amkz,amkst,gamkst
4594 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4595 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4596 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4598 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4600 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4603 CHARACTER NAMES(NMODE)*31
4604 REAL HV(4),PT(4),PN(4),PAA(4),PMULT(4,9)
4605 real*4 pr(4),pi1(4),pi2(4),pi3(4),pi4(4),pi5(4)
4606 real*8 amp1,amp2,amp3,amp4,amp5,ams1,ams2,amom,gamom
4607 real*8 am5sq,am4sq,am3sq,am2sq,am5,am4,am3
4609 real*8 gg1,gg2,gg3,ff1,ff2,ff3,ff4,alp,alp1,alp2
4618 DATA pi /3.141592653589793238462643/
4625 bwign(xm,am,gammab)=xm**2/cmplx(xm**2-am**2,gammab*am)
4627 bwign(xm,am,gamma)=xm**2/cmplx(xm**2-am**2,gamma*am)
4636 phspac=1./2**29/pi**14
4639 amp1=dcdmas(idffin(1,jnpi))
4640 amp2=dcdmas(idffin(2,jnpi))
4641 amp3=dcdmas(idffin(3,jnpi))
4642 amp4=dcdmas(idffin(4,jnpi))
4643 amp5=dcdmas(idffin(5,jnpi))
4658 ams1=(amp1+amp2+amp3+amp4+amp5)**2
4659 ams2=(amtau-amnuta)**2
4660 am5sq=ams1+ rr1*(ams2-ams1)
4662 phspac=phspac*(ams2-ams1)
4667 ams1=(amp2+amp3+amp4+amp5)**2
4669 am4sq=ams1+ rr1*(ams2-ams1)
4676 ams1=(amp2+amp3+amp4)**2
4678 alp1=atan((ams1-amom**2)/amom/gamom)
4679 alp2=atan((ams2-amom**2)/amom/gamom)
4680 alp=alp1+rr1*(alp2-alp1)
4681 am3sq =amom**2+amom*gamom*tan(alp)
4684 gg2=((am3sq-amom**2)**2+(amom*gamom)**2)/(amom*gamom)
4697 am2sq=ams1+ rr2*(ams2-ams1)
4703 enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4704 enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4705 ppi= enq1**2-amp3**2
4706 pppi=sqrt(abs(enq1**2-amp3**2))
4707 ff1=(4*pi)*(2*pppi/am2)
4709 call sphera(pppi,pi3)
4720 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4721 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4722 ppi = pr(4)**2-am2**2
4726 pi2(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4729 ff2=(4*pi)*(2*pr(3)/am3)
4731 exe=(pr(4)+pr(3))/am2
4732 call bostr3(exe,pi3,pi3)
4733 call bostr3(exe,pi4,pi4)
4736 thet =acos(-1.+2*rr3)
4738 call rotpol(thet,phi,pi2)
4739 call rotpol(thet,phi,pi3)
4740 call rotpol(thet,phi,pi4)
4746 pr(4)=1./(2*am4)*(am4**2+am3**2-amp5**2)
4747 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4748 ppi = pr(4)**2-am3**2
4752 pi5(4)=1./(2*am4)*(am4**2-am3**2+amp5**2)
4755 ff3=(4*pi)*(2*pr(3)/am4)
4757 exe=(pr(4)+pr(3))/am3
4758 call bostr3(exe,pi2,pi2)
4759 call bostr3(exe,pi3,pi3)
4760 call bostr3(exe,pi4,pi4)
4763 thet =acos(-1.+2*rr3)
4765 call rotpol(thet,phi,pi2)
4766 call rotpol(thet,phi,pi3)
4767 call rotpol(thet,phi,pi4)
4768 call rotpol(thet,phi,pi5)
4774 pr(4)=1./(2*am5)*(am5**2+am4**2-amp1**2)
4775 pr(3)= sqrt(abs(pr(4)**2-am4**2))
4776 ppi = pr(4)**2-am4**2
4780 pi1(4)=1./(2*am5)*(am5**2-am4**2+amp1**2)
4783 ff4=(4*pi)*(2*pr(3)/am5)
4785 exe=(pr(4)+pr(3))/am4
4786 call bostr3(exe,pi2,pi2)
4787 call bostr3(exe,pi3,pi3)
4788 call bostr3(exe,pi4,pi4)
4789 call bostr3(exe,pi5,pi5)
4792 thet =acos(-1.+2*rr3)
4794 call rotpol(thet,phi,pi1)
4795 call rotpol(thet,phi,pi2)
4796 call rotpol(thet,phi,pi3)
4797 call rotpol(thet,phi,pi4)
4798 call rotpol(thet,phi,pi5)
4807 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5sq)
4808 paa(3)= sqrt(abs(paa(4)**2-am5sq))
4809 ppi = paa(4)**2-am5sq
4810 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4814 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am5**2)
4817 phspac=phspac * gg1*gg2*gg3*ff1*ff2*ff3*ff4
4821 exe=(paa(4)+paa(3))/am5
4822 call bostr3(exe,pi1,pi1)
4823 call bostr3(exe,pi2,pi2)
4824 call bostr3(exe,pi3,pi3)
4825 call bostr3(exe,pi4,pi4)
4826 call bostr3(exe,pi5,pi5)
4834 qxn=paa(4)*pn(4)-paa(1)*pn(1)-paa(2)*pn(2)-paa(3)*pn(3)
4835 brak=2*(gv**2+ga**2)*(2*pxq*qxn+am5sq*pxn)
4836 & -6*(gv**2-ga**2)*amtau*amnuta*am5sq
4837 fompp = cabs(bwign(am3,amom,gamom))**2
4842 amplit=ccabib**2*gfermi**2/2. * brak
4843 amplit = amplit * fompp * fnorm
4846 dgamt=1/(2.*amtau)*amplit*phspac
4869 SUBROUTINE dphsrk(DGAMT,HV,PN,PR,PMULT,INUM)
4875 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4876 * ,ampiz,ampi,amro,gamro,ama1,gama1
4877 * ,amk,amkz,amkst,gamkst
4879 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4880 * ,ampiz,ampi,amro,gamro,ama1,gama1
4881 * ,amk,amkz,amkst,gamkst
4882 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4883 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4884 REAL HV(4),PT(4),PN(4),PR(4),PKC(4),PKZ(4),QQ(4),PMULT(4,9)
4890 DATA pi /3.141592653589793238462643/
4894 phspac=1./2**11/pi**5
4902 ams2=(amtau-amnuta)**2
4905 amx2=ams1+ rr1(1)*(ams2-ams1)
4907 phspac=phspac*(ams2-ams1)
4925 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
4926 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
4930 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
4932 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
4935 enq1=(amx2+amk**2-amkz**2)/(2.*amx)
4936 enq2=(amx2-amk**2+amkz**2)/(2.*amx)
4937 pppi=sqrt((enq1-amk)*(enq1+amk))
4938 phspac=phspac*(4*pi)*(2*pppi/amx)
4940 CALL sphera(pppi,pkc)
4946 exe=(pr(4)+pr(3))/amx
4948 CALL bostr3(exe,pkc,pkc)
4949 CALL bostr3(exe,pkz,pkz)
4951 30 qq(i)=pkc(i)-pkz(i)
4953 pksd =pr(4)*pr(4)-pr(3)*pr(3)-pr(2)*pr(2)-pr(1)*pr(1)
4954 qqpks=pr(4)* qq(4)-pr(3)* qq(3)-pr(2)* qq(2)-pr(1)* qq(1)
4956 31 qq(i)=qq(i)-pr(i)*qqpks/pksd
4959 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
4961 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
4962 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
4963 & +(gv**2-ga**2)*amtau*amnuta*qq2
4964 amplit=(gfermi*ccabib)**2*brak*2*fpirk(amx)
4965 dgamt=1/(2.*amtau)*amplit*phspac
4967 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
4978 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4979 * ,ampiz,ampi,amro,gamro,ama1,gama1
4980 * ,amk,amkz,amkst,gamkst
4982 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4983 * ,ampiz,ampi,amro,gamro,ama1,gama1
4984 * ,amk,amkz,amkst,gamkst
4987 fpirk=cabs(fpikm(w,amk,amkz))**2
4990 COMPLEX FUNCTION fpikmk(W,XM1,XM2)
4995 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
5000 IF (init.EQ.0 )
THEN
5013 fpikmk=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
5022 parameter(nmxhep=2000)
5023 common/hepevtx/nevhep,nhep,isthep(nmxhep),idhep(nmxhep),
5024 &jmohep(2,nmxhep),jdahep(2,nmxhep),phep(5,nmxhep),vhep(4,nmxhep)
5029 SUBROUTINE dwrph(KTO,PHX)
5033 IMPLICIT REAL*8 (a-h,o-z)
5044 IF (qhot(4).GT.1.e-5)
CALL dwluph(kto,qhot)
5047 SUBROUTINE dwluph(KTO,PHOT)
5058 COMMON /taupos/ np1,np2
5064 COMMON /taupos/ np1,np2
5068 IF (phot(4).LE.0.0)
RETURN
5071 IF((kto.EQ. 1).OR.(kto.EQ.11))
THEN
5078 IF(ktos.GT.10) ktos=ktos-10
5080 CALL tralo4(ktos,phot,phot,am)
5081 CALL filhep(0,1,22,nps,nps,0,0,phot,0.0,.true.)
5086 SUBROUTINE dwluel(KTO,ISGN,PNU,PWB,PEL,PNE)
5097 COMMON /taupos/ np1,np2
5100 REAL PNU(4),PWB(4),PEL(4),PNE(4)
5103 COMMON /taupos/ np1,np2
5114 CALL tralo4(kto,pnu,pnu,am)
5115 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5118 CALL tralo4(kto,pwb,pwb,am)
5122 CALL tralo4(kto,pel,pel,am)
5123 CALL filhep(0,1,11*isgn,nps,nps,0,0,pel,am,.false.)
5126 CALL tralo4(kto,pne,pne,am)
5127 CALL filhep(0,1,-12*isgn,nps,nps,0,0,pne,am,.true.)
5131 SUBROUTINE dwlumu(KTO,ISGN,PNU,PWB,PMU,PNM)
5142 COMMON /taupos/ np1,np2
5145 REAL PNU(4),PWB(4),PMU(4),PNM(4)
5148 COMMON /taupos/ np1,np2
5159 CALL tralo4(kto,pnu,pnu,am)
5160 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5163 CALL tralo4(kto,pwb,pwb,am)
5167 CALL tralo4(kto,pmu,pmu,am)
5168 CALL filhep(0,1,13*isgn,nps,nps,0,0,pmu,am,.false.)
5171 CALL tralo4(kto,pnm,pnm,am)
5172 CALL filhep(0,1,-14*isgn,nps,nps,0,0,pnm,am,.true.)
5176 SUBROUTINE dwlupi(KTO,ISGN,PPI,PNU)
5187 COMMON /taupos/ np1,np2
5197 CALL tralo4(kto,pnu,pnu,am)
5198 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5201 CALL tralo4(kto,ppi,ppi,am)
5202 CALL filhep(0,1,-211*isgn,nps,nps,0,0,ppi,am,.true.)
5206 SUBROUTINE dwluro(KTO,ISGN,PNU,PRHO,PIC,PIZ)
5217 COMMON /taupos/ np1,np2
5220 REAL PNU(4),PRHO(4),PIC(4),PIZ(4)
5223 COMMON /taupos/ np1,np2
5234 CALL tralo4(kto,pnu,pnu,am)
5235 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5238 CALL tralo4(kto,prho,prho,am)
5239 CALL filhep(0,2,-213*isgn,nps,nps,0,0,prho,am,.true.)
5242 CALL tralo4(kto,pic,pic,am)
5243 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pic,am,.true.)
5246 CALL tralo4(kto,piz,piz,am)
5247 CALL filhep(0,1,111,-2,-2,0,0,piz,am,.true.)
5251 SUBROUTINE dwluaa(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
5263 COMMON /taupos/ np1,np2
5266 REAL PNU(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
5269 COMMON /taupos/ np1,np2
5280 CALL tralo4(kto,pnu,pnu,am)
5281 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5284 CALL tralo4(kto,paa,paa,am)
5285 CALL filhep(0,1,-20213*isgn,nps,nps,0,0,paa,am,.true.)
5293 CALL tralo4(kto,pim2,pim2,am)
5294 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pim2,am,.true.)
5297 CALL tralo4(kto,pim1,pim1,am)
5298 CALL filhep(0,1,-211*isgn,-2,-2,0,0,pim1,am,.true.)
5301 CALL tralo4(kto,pipl,pipl,am)
5302 CALL filhep(0,1, 211*isgn,-3,-3,0,0,pipl,am,.true.)
5304 ELSE IF (jaa.EQ.2)
THEN
5309 CALL tralo4(kto,pim2,pim2,am)
5310 CALL filhep(0,1,111,-1,-1,0,0,pim2,am,.true.)
5313 CALL tralo4(kto,pim1,pim1,am)
5314 CALL filhep(0,1,111,-2,-2,0,0,pim1,am,.true.)
5317 CALL tralo4(kto,pipl,pipl,am)
5318 CALL filhep(0,1,-211*isgn,-3,-3,0,0,pipl,am,.true.)
5324 SUBROUTINE dwlukk (KTO,ISGN,PKK,PNU)
5334 COMMON /taupos/ np1,np2
5348 CALL tralo4 (kto,pnu,pnu,am)
5349 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5352 CALL tralo4 (kto,pkk,pkk,am)
5353 CALL filhep(0,1,-321*isgn,nps,nps,0,0,pkk,am,.true.)
5357 SUBROUTINE dwluks(KTO,ISGN,PNU,PKS,PKK,PPI,JKST)
5358 COMMON / taukle / bra1,brk0,brk0b,brks
5359 real*4 bra1,brk0,brk0b,brks
5361 COMMON /taupos/ np1,np2
5374 REAL PNU(4),PKS(4),PKK(4),PPI(4)
5376 REAL PNU(4),PKS(4),PKK(4),PPI(4),XIO(1)
5377 COMMON /taupos/ np1,np2
5388 CALL tralo4(kto,pnu,pnu,am)
5389 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5392 CALL tralo4(kto,pks,pks,am)
5393 CALL filhep(0,1,-323*isgn,nps,nps,0,0,pks,am,.true.)
5401 CALL tralo4(kto,ppi,ppi,am)
5402 CALL filhep(0,1,-211*isgn,-1,-1,0,0,ppi,am,.true.)
5405 IF (isgn.EQ.-1) bran=brk0
5408 IF(xio(1).GT.bran)
THEN
5414 CALL tralo4(kto,pkk,pkk,am)
5415 CALL filhep(0,1,k0type,-2,-2,0,0,pkk,am,.true.)
5417 ELSE IF(jkst.EQ.20)
THEN
5422 CALL tralo4(kto,ppi,ppi,am)
5423 CALL filhep(0,1,111,-1,-1,0,0,ppi,am,.true.)
5426 CALL tralo4(kto,pkk,pkk,am)
5427 CALL filhep(0,1,-321*isgn,-2,-2,0,0,pkk,am,.true.)
5433 SUBROUTINE dwlnew(KTO,ISGN,PNU,PWB,PNPI,MODE)
5443 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
5445 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5447 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5450 COMMON /taupos/ np1,np2
5451 CHARACTER NAMES(NMODE)*31
5452 REAL PNU(4),PWB(4),PNPI(4,9)
5464 CALL tralo4(kto,pnu,pnu,am)
5465 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5468 CALL tralo4(kto,pwb,pwb,am)
5469 CALL filhep(0,1,-24*isgn,nps,nps,0,0,pwb,am,.true.)
5478 kfpi=lunpik(idffin(i,jnpi), isgn)
5480 kfpi=lunpik(idffin(i,jnpi),-isgn)
5487 CALL tralo4(kto,ppi,ppi,am)
5488 CALL filhep(0,1,kfpi,-i,-i,0,0,ppi,am,.true.)
5493 #if defined (CePeCe)
5502 IMPLICIT REAL*8 (a-h,o-z)
5504 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5506 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5518 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5519 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5528 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5529 DATA pi /3.141592653589793238462643d0/
5531 IF(abs(y).LT.abs(x))
THEN
5533 IF(x.LE.0d0) the=pi-the
5535 the=acos(x/sqrt(x**2+y**2))
5546 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5547 DATA pi /3.141592653589793238462643d0/
5549 IF(abs(y).LT.abs(x))
THEN
5551 IF(x.LE.0d0) the=pi-the
5553 the=acos(x/sqrt(x**2+y**2))
5555 IF(y.LT.0d0) the=2d0*pi-the
5558 SUBROUTINE rotod1(PH1,PVEC,QVEC)
5563 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5564 dimension pvec(4),qvec(4),rvec(4)
5572 qvec(2)= cs*rvec(2)-sn*rvec(3)
5573 qvec(3)= sn*rvec(2)+cs*rvec(3)
5577 SUBROUTINE rotod2(PH1,PVEC,QVEC)
5582 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5583 dimension pvec(4),qvec(4),rvec(4)
5590 qvec(1)= cs*rvec(1)+sn*rvec(3)
5592 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5596 SUBROUTINE rotod3(PH1,PVEC,QVEC)
5601 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5603 dimension pvec(4),qvec(4),rvec(4)
5609 qvec(1)= cs*rvec(1)-sn*rvec(2)
5610 qvec(2)= sn*rvec(1)+cs*rvec(2)
5614 SUBROUTINE bostr3(EXE,PVEC,QVEC)
5620 real*4 pvec(4),qvec(4),rvec(4)
5633 SUBROUTINE bostd3(EXE,PVEC,QVEC)
5639 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5640 dimension pvec(4),qvec(4),rvec(4)
5654 SUBROUTINE rotor1(PH1,PVEC,QVEC)
5659 real*4 pvec(4),qvec(4),rvec(4)
5667 qvec(2)= cs*rvec(2)-sn*rvec(3)
5668 qvec(3)= sn*rvec(2)+cs*rvec(3)
5671 SUBROUTINE rotor2(PH1,PVEC,QVEC)
5676 IMPLICIT REAL*4(a-h,o-z)
5677 real*4 pvec(4),qvec(4),rvec(4)
5684 qvec(1)= cs*rvec(1)+sn*rvec(3)
5686 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5689 SUBROUTINE rotor3(PHI,PVEC,QVEC)
5694 real*4 pvec(4),qvec(4),rvec(4)
5700 qvec(1)= cs*rvec(1)-sn*rvec(2)
5701 qvec(2)= sn*rvec(1)+cs*rvec(2)
5705 SUBROUTINE spherd(R,X)
5710 real*8 r,x(4),pi,costh,sinth
5712 DATA pi /3.141592653589793238462643d0/
5716 sinth=sqrt(1 -costh**2)
5717 x(1)=r*sinth*cos(2*pi*rrr(2))
5718 x(2)=r*sinth*sin(2*pi*rrr(2))
5722 SUBROUTINE rotpox(THET,PHI,PP)
5723 IMPLICIT REAL*8 (a-h,o-z)
5733 CALL rotod2(thet,pp,pp)
5734 CALL rotod3( phi,pp,pp)
5737 SUBROUTINE sphera(R,X)
5745 DATA pi /3.141592653589793238462643/
5749 sinth=sqrt(1.-costh**2)
5750 x(1)=r*sinth*cos(2*pi*rrr(2))
5751 x(2)=r*sinth*sin(2*pi*rrr(2))
5755 SUBROUTINE rotpol(THET,PHI,PP)
5762 CALL rotor2(thet,pp,pp)
5763 CALL rotor3( phi,pp,pp)
5766 #include "../randg/tauola-random.h"
5769 IMPLICIT REAL*8(a-h,o-z)
5772 IF(x .LT.-1.0)
GO TO 1
5773 IF(x .LE. 0.5)
GO TO 2
5774 IF(x .EQ. 1.0)
GO TO 3
5775 IF(x .LE. 2.0)
GO TO 4
5779 z=z-0.5* log(abs(x))**2
5785 3 dilogt=1.64493406684822
5789 z=1.64493406684822 - log(x)* log(abs(t))
5790 5 y=2.66666666666666 *t+0.66666666666666
5791 b= 0.00000 00000 00001
5792 a=y*b +0.00000 00000 00004
5793 b=y*a-b+0.00000 00000 00011
5794 a=y*b-a+0.00000 00000 00037
5795 b=y*a-b+0.00000 00000 00121
5796 a=y*b-a+0.00000 00000 00398
5797 b=y*a-b+0.00000 00000 01312
5798 a=y*b-a+0.00000 00000 04342
5799 b=y*a-b+0.00000 00000 14437
5800 a=y*b-a+0.00000 00000 48274
5801 b=y*a-b+0.00000 00001 62421
5802 a=y*b-a+0.00000 00005 50291
5803 b=y*a-b+0.00000 00018 79117
5804 a=y*b-a+0.00000 00064 74338
5805 b=y*a-b+0.00000 00225 36705
5806 a=y*b-a+0.00000 00793 87055
5807 b=y*a-b+0.00000 02835 75385
5808 a=y*b-a+0.00000 10299 04264
5809 b=y*a-b+0.00000 38163 29463
5810 a=y*b-a+0.00001 44963 00557
5811 b=y*a-b+0.00005 68178 22718
5812 a=y*b-a+0.00023 20021 96094
5813 b=y*a-b+0.00100 16274 96164
5814 a=y*b-a+0.00468 63619 59447
5815 b=y*a-b+0.02487 93229 24228
5816 a=y*b-a+0.16607 30329 27855
5817 a=y*a-b+1.93506 43008 6996