36 COMMON / taubra / gamprt(500),jlist(500),nchan
40 REAL CUMUL(500),RRR(1)
42 IF(nchan.LE.0.OR.nchan.GT.500)
GOTO 902
49 IF(rrr(1).LT.cumul(i)/cumul(nchan)) ji=i
54 9020
FORMAT(
' ----- JAKER: WRONG NCHAN')
57 SUBROUTINE dekay(KTO,HX)
84 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
88 COMMON /taupos/ np1,np2
89 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
91 parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
93 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
96 CHARACTER NAMES(NMODE)*31
97 COMMON / inout / inut,iout
98 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4),HDUM(4)
110 IF (iwarm.EQ.1) x=5/(iwarm-1)
112 WRITE(iout,7001) jak1,jak2
116 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
117 CALL dadmel(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
118 CALL dadmmu(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
119 CALL dadmpi(-1,idum,pdum,pdum1,pdum2)
120 CALL dadmro(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4)
121 CALL dadmaa(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5,jdum)
122 CALL dadmkk(-1,idum,pdum,pdum1,pdum2)
123 CALL dadmks(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,jdum)
124 CALL dadnew(-1,idum,hdum,pdum1,pdum2,pdumx,jdum)
130 ELSEIF(kto.EQ.1)
THEN
134 IF(iwarm.EQ.0)
GOTO 902
140 CALL dekay1(0,h,isgn)
141 ELSEIF(kto.EQ.2)
THEN
145 IF(iwarm.EQ.0)
GOTO 902
151 CALL dekay2(0,h,isgn)
152 ELSEIF(kto.EQ.11)
THEN
157 CALL dekay1(1,h,isgn)
158 ELSEIF(kto.EQ.12)
THEN
163 CALL dekay2(1,h,isgn)
164 ELSEIF(kto.EQ.100)
THEN
166 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
167 CALL dadmel( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
168 CALL dadmmu( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
169 CALL dadmpi( 1,idum,pdum,pdum1,pdum2)
170 CALL dadmro( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4)
171 CALL dadmaa( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5,jdum)
172 CALL dadmkk( 1,idum,pdum,pdum1,pdum2)
173 CALL dadmks( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,jdum)
174 CALL dadnew( 1,idum,hdum,pdum1,pdum2,pdumx,jdum)
175 WRITE(iout,7010) nev1,nev2,nevtot
176 WRITE(iout,7011) (nevdec(i),gampmc(i),gamper(i),i= 1,7)
178 $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
189 7001
FORMAT(///1x,15(5h*****)
191 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.9 ******',9x,1h*,
192 $ /,
' *', 25x,
'***********Sep 2005***************',9x,1h*,
193 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
194 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
195 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
196 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
197 $ /,
' *', 25x,
' Physics initialization by CLEO collab ',9x,1h*,
198 $ /,
' *', 25x,
' see Alain Weinstein www home page: ',9x,1h*,
199 $ /,
' *', 25x,
'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
200 $ /,
' *', 25x,
'/korb_doc.html#files ',9x,1h*,
201 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
202 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
204 $ /,
' *', 25x,
'****DEKAY ROUTINE: INITIALIZATION******',9x,1h*,
205 $ /,
' *',i20 ,5x,
'JAK1 = DECAY MODE TAU+ ',9x,1h*,
206 $ /,
' *',i20 ,5x,
'JAK2 = DECAY MODE TAU- ',9x,1h*,
208 7010
FORMAT(///1x,15(5h*****)
210 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.9 ******',9x,1h*,
211 $ /,
' *', 25x,
'***********Sep 2005***************',9x,1h*,
213 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
214 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
215 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
216 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
217 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
218 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
219 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
220 $ /,
' *', 25x,
'*****DEKAY ROUTINE: FINAL REPORT*******',9x,1h*,
221 $ /,
' *',i20 ,5x,
'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*,
222 $ /,
' *',i20 ,5x,
'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*,
223 $ /,
' *',i20 ,5x,
'NEVTOT = SUM ',9x,1h*,
225 $
' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
227 $ ,i10,2f12.7 ,
' DADMEL ELECTRON ',9x,1h*
228 $ /,
' *',i10,2f12.7 ,
' DADMMU MUON ',9x,1h*
229 $ /,
' *',i10,2f12.7 ,
' DADMPI PION ',9x,1h*
230 $ /,
' *',i10,2f12.7,
' DADMRO RHO (->2PI) ',9x,1h*
231 $ /,
' *',i10,2f12.7,
' DADMAA A1 (->3PI) ',9x,1h*
232 $ /,
' *',i10,2f12.7,
' DADMKK KAON ',9x,1h*
233 $ /,
' *',i10,2f12.7,
' DADMKS K* ',9x,1h*)
235 $ ,i10,2f12.7,a31 ,8x,1h*)
237 $ ,20x,
'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
238 $ /,
' *',20x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
241 9020
FORMAT(
' ----- DEKAY: LACK OF INITIALISATION')
244 9100
FORMAT(
' ----- DEKAY: WRONG VALUE OF KTO ')
247 SUBROUTINE dekay1(IMOD,HH,ISGN)
251 COMMON / decp4 / pp1(4),pp2(4),kf1,kf2
252 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
253 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
254 real*4 gampmc ,gamper
257 REAL HV(4),PNU(4),PPI(4)
258 REAL PWB(4),PMU(4),PNM(4)
259 REAL PRHO(4),PIC(4),PIZ(4)
260 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
267 IF(jak1.EQ.-1)
RETURN
272 IF(jak1.EQ.0)
CALL jaker(jak)
274 CALL dadmel(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
275 ELSEIF(jak.EQ.2)
THEN
276 CALL dadmmu(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
277 ELSEIF(jak.EQ.3)
THEN
278 CALL dadmpi(0, isgn,hv,ppi,pnu)
279 ELSEIF(jak.EQ.4)
THEN
280 CALL dadmro(0, isgn,hv,pnu,prho,pic,piz)
281 ELSEIF(jak.EQ.5)
THEN
282 CALL dadmaa(0, isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
283 ELSEIF(jak.EQ.6)
THEN
284 CALL dadmkk(0, isgn,hv,pkk,pnu)
285 ELSEIF(jak.EQ.7)
THEN
286 CALL dadmks(0, isgn,hv,pnu,pks ,pkk,ppi,jkst)
288 CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
294 ELSEIF(imd.EQ.1)
THEN
298 nevdec(jak)=nevdec(jak)+1
303 CALL dwluel(1,isgn,pnu,pwb,pmu,pnm)
304 CALL dwrph(ktom,phot)
308 ELSEIF(jak.EQ.2)
THEN
309 CALL dwlumu(1,isgn,pnu,pwb,pmu,pnm)
310 CALL dwrph(ktom,phot)
314 ELSEIF(jak.EQ.3)
THEN
315 CALL dwlupi(1,isgn,ppi,pnu)
319 ELSEIF(jak.EQ.4)
THEN
320 CALL dwluro(1,isgn,pnu,prho,pic,piz)
324 ELSEIF(jak.EQ.5)
THEN
325 CALL dwluaa(1,isgn,pnu,paa,pim1,pim2,pipl,jaa)
328 ELSEIF(jak.EQ.6)
THEN
329 CALL dwlukk(1,isgn,pkk,pnu)
332 ELSEIF(jak.EQ.7)
THEN
333 CALL dwluks(1,isgn,pnu,pks,pkk,ppi,jkst)
338 CALL dwlnew(1,isgn,pnu,pwb,pnpi,jak)
346 SUBROUTINE dekay2(IMOD,HH,ISGN)
350 COMMON / decp4 / pp1(4),pp2(4),kf1,kf2
351 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
352 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
353 real*4 gampmc ,gamper
356 REAL HV(4),PNU(4),PPI(4)
357 REAL PWB(4),PMU(4),PNM(4)
358 REAL PRHO(4),PIC(4),PIZ(4)
359 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
366 IF(jak2.EQ.-1)
RETURN
371 IF(jak2.EQ.0)
CALL jaker(jak)
373 CALL dadmel(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
374 ELSEIF(jak.EQ.2)
THEN
375 CALL dadmmu(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
376 ELSEIF(jak.EQ.3)
THEN
377 CALL dadmpi(0, isgn,hv,ppi,pnu)
378 ELSEIF(jak.EQ.4)
THEN
379 CALL dadmro(0, isgn,hv,pnu,prho,pic,piz)
380 ELSEIF(jak.EQ.5)
THEN
381 CALL dadmaa(0, isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
382 ELSEIF(jak.EQ.6)
THEN
383 CALL dadmkk(0, isgn,hv,pkk,pnu)
384 ELSEIF(jak.EQ.7)
THEN
385 CALL dadmks(0, isgn,hv,pnu,pks ,pkk,ppi,jkst)
387 CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
392 ELSEIF(imd.EQ.1)
THEN
396 nevdec(jak)=nevdec(jak)+1
401 CALL dwluel(2,isgn,pnu,pwb,pmu,pnm)
402 CALL dwrph(ktom,phot)
406 ELSEIF(jak.EQ.2)
THEN
407 CALL dwlumu(2,isgn,pnu,pwb,pmu,pnm)
408 CALL dwrph(ktom,phot)
412 ELSEIF(jak.EQ.3)
THEN
413 CALL dwlupi(2,isgn,ppi,pnu)
417 ELSEIF(jak.EQ.4)
THEN
418 CALL dwluro(2,isgn,pnu,prho,pic,piz)
422 ELSEIF(jak.EQ.5)
THEN
423 CALL dwluaa(2,isgn,pnu,paa,pim1,pim2,pipl,jaa)
426 ELSEIF(jak.EQ.6)
THEN
427 CALL dwlukk(2,isgn,pkk,pnu)
430 ELSEIF(jak.EQ.7)
THEN
431 CALL dwluks(2,isgn,pnu,pks,pkk,ppi,jkst)
436 CALL dwlnew(2,isgn,pnu,pwb,pnpi,jak)
444 SUBROUTINE dexay(KTO,POL)
458 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
459 real*4 gampmc ,gamper
460 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
462 COMMON /taupos/ np1,np2
463 parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
465 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
468 CHARACTER NAMES(NMODE)*31
469 COMMON / inout / inut,iout
471 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
485 WRITE(iout, 7001) jak1,jak2
489 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
490 CALL dexel(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
491 CALL dexmu(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
492 CALL dexpi(-1,idum,pdum,pdum1,pdum2)
493 CALL dexro(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4)
494 CALL dexaa(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5,idum)
495 CALL dexkk(-1,idum,pdum,pdum1,pdum2)
496 CALL dexks(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,idum)
497 CALL dexnew(-1,idum,pdum,pdum1,pdum2,pdumi,idum)
503 ELSEIF(kto.EQ.1)
THEN
508 IF(iwarm.EQ.0)
GOTO 902
511 CALL dexay1(kto,jak1,jakp,pol,isgn)
512 ELSEIF(kto.EQ.2)
THEN
517 IF(iwarm.EQ.0)
GOTO 902
518 isgn=-idff/iabs(idff)
520 CALL dexay1(kto,jak2,jakm,pol,isgn)
521 ELSEIF(kto.EQ.100)
THEN
523 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
524 CALL dexel( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
525 CALL dexmu( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
526 CALL dexpi( 1,idum,pdum,pdum1,pdum2)
527 CALL dexro( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4)
528 CALL dexaa( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5,idum)
529 CALL dexkk( 1,idum,pdum,pdum1,pdum2)
530 CALL dexks( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,idum)
531 CALL dexnew( 1,idum,pdum,pdum1,pdum2,pdumi,idum)
532 WRITE(iout,7010) nev1,nev2,nevtot
533 WRITE(iout,7011) (nevdec(i),gampmc(i),gamper(i),i= 1,7)
535 $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
542 7001
FORMAT(///1x,15(5h*****)
544 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.9 ******',9x,1h*,
545 $ /,
' *', 25x,
'*********** Sep 2005***************',9x,1h*,
546 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
547 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
548 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
549 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
550 $ /,
' *', 25x,
' Physics initialization by CLEO collab ',9x,1h*,
551 $ /,
' *', 25x,
' see Alain Weinstein www home page: ',9x,1h*,
552 $ /,
' *', 25x,
'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
553 $ /,
' *', 25x,
'/korb_doc.html#files ',9x,1h*,
554 $ /,
' *', 25x,
'*******CERN TH-6793 NOVEMBER 1992*****',9x,1h*,
555 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
557 $ /,
' *', 25x,
'*******CERN-TH-6793 NOVEMBER 1992*****',9x,1h*,
558 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
559 $ /,
' *', 25x,
'******DEXAY ROUTINE: INITIALIZATION****',9x,1h*
560 $ /,
' *',i20 ,5x,
'JAK1 = DECAY MODE FERMION1 (TAU+) ',9x,1h*
561 $ /,
' *',i20 ,5x,
'JAK2 = DECAY MODE FERMION2 (TAU-) ',9x,1h*
565 7010
FORMAT(///1x,15(5h*****)
567 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.9 ******',9x,1h*,
568 $ /,
' *', 25x,
'***********Sep 2005***************',9x,1h*,
569 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
570 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
571 $ /,
' *', 25x,
'**AVAILABLE FROM: WASM AT CERNVM ******',9x,1h*,
572 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
573 $ /,
' *', 25x,
'*******CERN-TH-5856 SEPTEMBER 1990*****',9x,1h*,
574 $ /,
' *', 25x,
'*******CERN-TH-6195 SEPTEMBER 1991*****',9x,1h*,
575 $ /,
' *', 25x,
'*******CERN-TH-6793 NOVEMBER 1992*****',9x,1h*,
576 $ /,
' *', 25x,
'******DEXAY ROUTINE: FINAL REPORT******',9x,1h*
577 $ /,
' *',i20 ,5x,
'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*
578 $ /,
' *',i20 ,5x,
'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*
579 $ /,
' *',i20 ,5x,
'NEVTOT = SUM ',9x,1h*
581 $
' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
583 $ ,i10,2f12.7 ,
' DADMEL ELECTRON ',9x,1h*
584 $ /,
' *',i10,2f12.7 ,
' DADMMU MUON ',9x,1h*
585 $ /,
' *',i10,2f12.7 ,
' DADMPI PION ',9x,1h*
586 $ /,
' *',i10,2f12.7,
' DADMRO RHO (->2PI) ',9x,1h*
587 $ /,
' *',i10,2f12.7,
' DADMAA A1 (->3PI) ',9x,1h*
588 $ /,
' *',i10,2f12.7,
' DADMKK KAON ',9x,1h*
589 $ /,
' *',i10,2f12.7,
' DADMKS K* ',9x,1h*)
591 $ ,i10,2f12.7,a31 ,8x,1h*)
593 $ ,20x,
'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
594 $ /,
' *',20x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
596 902
WRITE(iout, 9020)
597 9020
FORMAT(
' ----- DEXAY: LACK OF INITIALISATION')
599 910
WRITE(iout, 9100)
600 9100
FORMAT(
' ----- DEXAY: WRONG VALUE OF KTO ')
603 SUBROUTINE dexay1(KTO,JAKIN,JAK,POL,ISGN)
609 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
610 real*4 gampmc ,gamper
611 COMMON / inout / inut,iout
614 REAL PRHO(4),PIC(4),PIZ(4)
615 REAL PWB(4),PMU(4),PNM(4)
616 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
622 IF(jakin.EQ.-1)
RETURN
629 IF(jak.EQ.0)
CALL jaker(jak)
632 CALL dexel(0, isgn,polar,pnu,pwb,pmu,pnm,phot)
633 CALL dwluel(kto,isgn,pnu,pwb,pmu,pnm)
634 CALL dwrph(kto,phot )
635 ELSEIF(jak.EQ.2)
THEN
636 CALL dexmu(0, isgn,polar,pnu,pwb,pmu,pnm,phot)
637 CALL dwlumu(kto,isgn,pnu,pwb,pmu,pnm)
638 CALL dwrph(kto,phot )
639 ELSEIF(jak.EQ.3)
THEN
640 CALL dexpi(0, isgn,polar,ppi,pnu)
641 CALL dwlupi(kto,isgn,ppi,pnu)
642 ELSEIF(jak.EQ.4)
THEN
643 CALL dexro(0, isgn,polar,pnu,prho,pic,piz)
644 CALL dwluro(kto,isgn,pnu,prho,pic,piz)
645 ELSEIF(jak.EQ.5)
THEN
646 CALL dexaa(0, isgn,polar,pnu,paa,pim1,pim2,pipl,jaa)
647 CALL dwluaa(kto,isgn,pnu,paa,pim1,pim2,pipl,jaa)
648 ELSEIF(jak.EQ.6)
THEN
649 CALL dexkk(0, isgn,polar,pkk,pnu)
650 CALL dwlukk(kto,isgn,pkk,pnu)
651 ELSEIF(jak.EQ.7)
THEN
652 CALL dexks(0, isgn,polar,pnu,pks,pkk,ppi,jkst)
653 CALL dwluks(kto,isgn,pnu,pks,pkk,ppi,jkst)
656 CALL dexnew(0, isgn,polar,pnu,pwb,pnpi,jnpi)
657 CALL dwlnew(kto,isgn,pnu,pwb,pnpi,jak)
659 nevdec(jak)=nevdec(jak)+1
661 SUBROUTINE dexel(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
668 REAL POL(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4),PH(4),RN(1)
674 CALL dadmel( -1,isgn,hv,pnu,pwb,q1,q2,ph)
677 ELSEIF(mode.EQ. 0)
THEN
680 IF(iwarm.EQ.0)
GOTO 902
681 CALL dadmel( 0,isgn,hv,pnu,pwb,q1,q2,ph)
682 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
685 IF(rn(1).GT.wt)
GOTO 300
687 ELSEIF(mode.EQ. 1)
THEN
689 CALL dadmel( 1,isgn,hv,pnu,pwb,q1,q2,ph)
695 9020
FORMAT(
' ----- DEXEL: LACK OF INITIALISATION')
698 SUBROUTINE dexmu(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
707 COMMON / inout / inut,iout
708 REAL POL(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4),PH(4),RN(1)
714 CALL dadmmu( -1,isgn,hv,pnu,pwb,q1,q2,ph)
717 ELSEIF(mode.EQ. 0)
THEN
720 IF(iwarm.EQ.0)
GOTO 902
721 CALL dadmmu( 0,isgn,hv,pnu,pwb,q1,q2,ph)
722 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
725 IF(rn(1).GT.wt)
GOTO 300
727 ELSEIF(mode.EQ. 1)
THEN
729 CALL dadmmu( 1,isgn,hv,pnu,pwb,q1,q2,ph)
734 902
WRITE(iout, 9020)
735 9020
FORMAT(
' ----- DEXMU: LACK OF INITIALISATION')
738 SUBROUTINE dadmel(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
743 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
744 real*4 gfermi,gv,ga,ccabib,scabib,gamel
745 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
746 * ,ampiz,ampi,amro,gamro,ama1,gama1
747 * ,amk,amkz,amkst,gamkst
749 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
750 * ,ampiz,ampi,amro,gamro,ama1,gama1
751 * ,amk,amkz,amkst,gamkst
752 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
753 real*4 gampmc ,gamper
757 COMMON / inout / inut,iout
760 REAL HHV(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4)
761 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
764 DATA pi /3.141592653589793238462643/
777 CALL dphsel(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
778 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
782 ELSEIF(mode.EQ. 0)
THEN
785 IF(iwarm.EQ.0)
GOTO 902
787 CALL dphsel(wt,hv,pnu,pwb,q1,q2,phx)
793 IF(wt.GT.wtmax) nevovr=nevovr+1
794 IF(rn*wtmax.GT.wt)
GOTO 300
801 CALL rotor2(thet,pnu,pnu)
802 CALL rotor3( phi,pnu,pnu)
803 CALL rotor2(thet,pwb,pwb)
804 CALL rotor3( phi,pwb,pwb)
805 CALL rotor2(thet,q1,q1)
806 CALL rotor3( phi,q1,q1)
807 CALL rotor2(thet,q2,q2)
808 CALL rotor3( phi,q2,q2)
809 CALL rotor2(thet,hv,hv)
810 CALL rotor3( phi,hv,hv)
811 CALL rotor2(thet,phx,phx)
812 CALL rotor3( phi,phx,phx)
814 44 hhv(i)=-isgn*hv(i)
817 ELSEIF(mode.EQ. 1)
THEN
819 IF(nevraw.EQ.0)
RETURN
820 pargam=swt/float(nevraw+1)
822 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
824 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
832 7010
FORMAT(///1x,15(5h*****)
833 $ /,
' *', 25x,
'******** DADMEL FINAL REPORT ******** ',9x,1h*
834 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF EL DECAYS TOTAL ',9x,1h*
835 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF EL DECS. ACCEPTED ',9x,1h*
836 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
837 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( ELECTRON) IN GEV UNITS ',9x,1h*
838 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
839 $ /,
' *',f20.9,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
840 $ /,
' *',25x,
'COMPLETE QED CORRECTIONS INCLUDED ',9x,1h*
841 $ /,
' *',25x,
'BUT ONLY V-A CUPLINGS ',9x,1h*
843 902
WRITE(iout, 9020)
844 9020
FORMAT(
' ----- DADMEL: LACK OF INITIALISATION')
847 SUBROUTINE dadmmu(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
849 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
850 * ,ampiz,ampi,amro,gamro,ama1,gama1
851 * ,amk,amkz,amkst,gamkst
853 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
854 * ,ampiz,ampi,amro,gamro,ama1,gama1
855 * ,amk,amkz,amkst,gamkst
856 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
857 real*4 gfermi,gv,ga,ccabib,scabib,gamel
858 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
859 real*4 gampmc ,gamper
860 COMMON / inout / inut,iout
862 REAL HHV(4),HV(4),PNU(4),PWB(4),Q1(4),Q2(4)
863 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
866 DATA pi /3.141592653589793238462643/
879 CALL dphsmu(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
880 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
884 ELSEIF(mode.EQ. 0)
THEN
887 IF(iwarm.EQ.0)
GOTO 902
889 CALL dphsmu(wt,hv,pnu,pwb,q1,q2,phx)
895 IF(wt.GT.wtmax) nevovr=nevovr+1
896 IF(rn*wtmax.GT.wt)
GOTO 300
901 CALL rotor2(thet,pnu,pnu)
902 CALL rotor3( phi,pnu,pnu)
903 CALL rotor2(thet,pwb,pwb)
904 CALL rotor3( phi,pwb,pwb)
905 CALL rotor2(thet,q1,q1)
906 CALL rotor3( phi,q1,q1)
907 CALL rotor2(thet,q2,q2)
908 CALL rotor3( phi,q2,q2)
909 CALL rotor2(thet,hv,hv)
910 CALL rotor3( phi,hv,hv)
911 CALL rotor2(thet,phx,phx)
912 CALL rotor3( phi,phx,phx)
914 44 hhv(i)=-isgn*hv(i)
917 ELSEIF(mode.EQ. 1)
THEN
919 IF(nevraw.EQ.0)
RETURN
920 pargam=swt/float(nevraw+1)
922 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
924 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
932 7010
FORMAT(///1x,15(5h*****)
933 $ /,
' *', 25x,
'******** DADMMU FINAL REPORT ******** ',9x,1h*
934 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF MU DECAYS TOTAL ',9x,1h*
935 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF MU DECS. ACCEPTED ',9x,1h*
936 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
937 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (MU DECAY) IN GEV UNITS ',9x,1h*
938 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
939 $ /,
' *',f20.9,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
940 $ /,
' *',25x,
'COMPLETE QED CORRECTIONS INCLUDED ',9x,1h*
941 $ /,
' *',25x,
'BUT ONLY V-A CUPLINGS ',9x,1h*
943 902
WRITE(iout, 9020)
944 9020
FORMAT(
' ----- DADMMU: LACK OF INITIALISATION')
947 SUBROUTINE dphsel(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
953 real*4 hvx(4),paax(4),xax(4),qpx(4),xnx(4)
954 real*8 hv(4),ph(4),paa(4),xa(4),qp(4),xn(4)
957 CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
968 SUBROUTINE dphsmu(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
974 real*4 hvx(4),paax(4),xax(4),qpx(4),xnx(4)
975 real*8 hv(4),ph(4),paa(4),xa(4),qp(4),xn(4)
978 CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
989 SUBROUTINE drcmu(DGAMT,HV,PH,PAA,XA,QP,XN,IELMU)
990 IMPLICIT REAL*8 (a-h,o-z)
995 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
996 * ,ampiz,ampi,amro,gamro,ama1,gama1
997 * ,amk,amkz,amkst,gamkst
999 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1000 * ,ampiz,ampi,amro,gamro,ama1,gama1
1001 * ,amk,amkz,amkst,gamkst
1002 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1003 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1005 COMMON / inout / inut,iout
1006 COMMON / taurad / xk0dec,itdkrc
1008 real*8 hv(4),pt(4),ph(4),paa(4),xa(4),qp(4),xn(4)
1012 DATA pi /3.141592653589793238462643d0/
1020 phspac=1./2**17/pi**8
1030 IF (ielmu.EQ.1)
THEN
1037 IF ( itdkrc.EQ.0) prhard=0d0
1039 IF(prsoft.LT.0.1)
THEN
1040 print *,
'ERROR IN DRCMU; PRSOFT=',prsoft
1045 ihard=(rr5.GT.prsoft)
1049 ams1=(amu+amnuta)**2
1052 xl1=log(xk1/2/xk0dec)
1057 phspac=phspac*ams2*xl1*xk
1058 phspac=phspac/prhard
1061 phspac=phspac*2**6*pi**3
1062 phspac=phspac/prsoft
1071 am2sq=ams1+ rr2*(ams2-ams1)
1073 phspac=phspac*(ams2-ams1)
1075 enq1=(am2sq+amnuta**2)/(2*am2)
1076 enq2=(am2sq-amnuta**2)/(2*am2)
1077 ppi= enq1**2-amnuta**2
1078 pppi=sqrt(abs(enq1**2-amnuta**2))
1079 phspac=phspac*(4*pi)*(2*pppi/am2)
1081 CALL spherd(pppi,xn)
1091 pr(4)=1.d0/(2*am3)*(am3**2+am2**2-amu**2)
1092 pr(3)= sqrt(abs(pr(4)**2-am2**2))
1093 ppi = pr(4)**2-am2**2
1097 qp(4)=1.d0/(2*am3)*(am3**2-am2**2+amu**2)
1099 phspac=phspac*(4*pi)*(2*pr(3)/am3)
1101 exe=(pr(4)+pr(3))/am2
1102 CALL bostd3(exe,xn,xn)
1103 CALL bostd3(exe,xa,xa)
1107 eps=4*(amu/amtax)**2
1108 xl1=log((2+eps)/eps)
1110 eta =exp(xl1*rr3+xl0)
1113 phspac=phspac*xl1/2*eta
1115 CALL rotpox(thet,phi,xn)
1116 CALL rotpox(thet,phi,xa)
1117 CALL rotpox(thet,phi,qp)
1118 CALL rotpox(thet,phi,pr)
1124 paa(4)=1/(2*amtax)*(amtax**2+am3**2)
1125 paa(3)= sqrt(abs(paa(4)**2-am3**2))
1126 ppi = paa(4)**2-am3**2
1127 phspac=phspac*(4*pi)*(2*paa(3)/amtax)
1135 exe=(paa(4)+paa(3))/am3
1136 CALL bostd3(exe,xn,xn)
1137 CALL bostd3(exe,xa,xa)
1138 CALL bostd3(exe,qp,qp)
1139 CALL bostd3(exe,pr,pr)
1141 thet =acos(-1.+2*rr3)
1143 CALL rotpox(thet,phi,xn)
1144 CALL rotpox(thet,phi,xa)
1145 CALL rotpox(thet,phi,qp)
1146 CALL rotpox(thet,phi,pr)
1161 CALL dampry(itdkrc,xk0dec,ph,xa,qp,xn,amplit,hv)
1162 dgamt=1/(2.*amtax)*amplit*phspac
1164 SUBROUTINE dampry(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1165 IMPLICIT REAL*8 (a-h,o-z)
1171 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1172 * ,ampiz,ampi,amro,gamro,ama1,gama1
1173 * ,amk,amkz,amkst,gamkst
1175 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1176 * ,ampiz,ampi,amro,gamro,ama1,gama1
1177 * ,amk,amkz,amkst,gamkst
1178 real*8 hv(4),qp(4),xn(4),xa(4),xk(4)
1182 IF(xk(4).LT.0.1d0*ak0)
THEN
1183 amplit=thb(itdkrc,qp,xn,xa,ak0,hv)
1185 amplit=sqm2(itdkrc,qp,xn,xa,xk,ak0,hv)
1189 FUNCTION sqm2(ITDKRC,QP,XN,XA,XK,AK0,HV)
1202 IMPLICIT REAL*8(a-h,o-z)
1203 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1204 * ,ampiz,ampi,amro,gamro,ama1,gama1
1205 * ,amk,amkz,amkst,gamkst
1207 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1208 * ,ampiz,ampi,amro,gamro,ama1,gama1
1209 * ,amk,amkz,amkst,gamkst
1210 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1211 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1212 COMMON / qedprm /alfinv,alfpi,xk0
1213 real*8 alfinv,alfpi,xk0
1214 real*8 qp(4),xn(4),xa(4),xk(4)
1217 real*8 s0(3),rxa(3),rxk(3),rqp(3)
1218 DATA pi /3.141592653589793238462643d0/
1224 emass2=qp(4)**2-qp(1)**2-qp(2)**2-qp(3)**2
1232 rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1234 rxk(i)=r(4)*xk(4)-r(1)*xk(1)-r(2)*xk(2)-r(3)*xk(3)
1235 rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
1237 qpxn=qp(4)*xn(4)-qp(1)*xn(1)-qp(2)*xn(2)-qp(3)*xn(3)
1238 qpxa=qp(4)*xa(4)-qp(1)*xa(1)-qp(2)*xa(2)-qp(3)*xa(3)
1239 qpxk=qp(4)*xk(4)-qp(1)*xk(1)-qp(2)*xk(2)-qp(3)*xk(3)
1241 xnxk=xn(4)*xk(4)-xn(1)*xk(1)-xn(2)*xk(2)-xn(3)*xk(3)
1242 xaxk=xa(4)*xk(4)-xa(1)*xk(1)-xa(2)*xk(2)-xa(3)*xk(3)
1252 s1= qpxn*txa*( -emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
1254 $qpxn/txk**2* ( tmass2*xaxk - txa*txk+ xaxk*txk) -
1255 $txa*txn/txk - qpxn/(qpxk*txk)* (tqp*xaxk-txk*qpxa)
1256 const4=256*pi/alphai*gf**2
1257 IF (itdkrc.EQ.0) const4=0d0
1260 s0(i) = qpxn*rxa(i)*(-emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
1262 $ qpxn/txk**2* (tmass2*xaxk - txa*rxk(i)+ xaxk*rxk(i))-
1263 $ rxa(i)*txn/txk - qpxn/(qpxk*txk)*(rqp(i)*xaxk- rxk(i)*qpxa)
1264 5 hv(i)=s0(i)/s1-1.d0
1267 FUNCTION thb(ITDKRC,QP,XN,XA,AK0,HV)
1281 IMPLICIT REAL*8(a-h,o-z)
1282 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1283 * ,ampiz,ampi,amro,gamro,ama1,gama1
1284 * ,amk,amkz,amkst,gamkst
1286 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1287 * ,ampiz,ampi,amro,gamro,ama1,gama1
1288 * ,amk,amkz,amkst,gamkst
1289 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1290 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1291 COMMON / qedprm /alfinv,alfpi,xk0
1292 real*8 alfinv,alfpi,xk0
1293 dimension qp(4),xn(4),xa(4)
1296 real*8 rxa(3),rxn(3),rqp(3)
1297 real*8 bornpl(3),am3pol(3),xm3pol(3)
1298 DATA pi /3.141592653589793238462643d0/
1311 rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1312 rxn(i)=r(4)*xn(4)-r(1)*xn(1)-r(2)*xn(2)-r(3)*xn(3)
1314 rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
1318 u3=sqrt(qp(1)**2+qp(2)**2+qp(3)**2)/tmass
1320 w0=(xn(4)+xa(4))/tmass
1332 f0=2*u0/u3*( dilogt(1-(um*wm/(up*wp)))- dilogt(1-wm/wp) +
1333 $dilogt(1-um/up) -2*yu+ 2*log(up)*(yw+yu) ) +
1334 $1/y* ( 2*u3*yu + (1-eps2- 2*y)*log(eps) ) +
1335 $ 2 - 4*(u0/u3*yu -1)* log(2*al)
1336 fp= yu/(2*u3)*(1 + (1-eps2)/y ) + log(eps)/y
1337 fm= yu/(2*u3)*(1 - (1-eps2)/y ) - log(eps)/y
1340 qpxn=qp(4)*xn(4)-qp(1)*xn(1)-qp(2)*xn(2)-qp(3)*xn(3)
1341 qpxa=qp(4)*xa(4)-qp(1)*xa(1)-qp(2)*xa(2)-qp(3)*xa(3)
1342 xnxa=xn(4)*xa(4)-xn(1)*xa(1)-xn(2)*xa(2)-xn(3)*xa(3)
1347 const3=1/(2*alphai*pi)*64*gf**2
1348 IF (itdkrc.EQ.0) const3=0d0
1349 xm3= -( f0* qpxn*txa + fp*eps2* txn*txa +
1350 $fm* qpxn*qpxa + f3* tmass2*xnxa )
1353 brak= (gv+ga)**2*tqp*xnxa+(gv-ga)**2*txa*qpxn
1354 & -(gv**2-ga**2)*tmass*amnuta*qpxa
1355 born= 32*(gfermi**2/2.)*brak
1357 xm3pol(i)= -( f0* qpxn*rxa(i) + fp*eps2* txn*rxa(i) +
1358 $ fm* qpxn* (qpxa + (rxa(i)*tqp-txa*rqp(i))/tmass2 ) +
1359 $ f3* (tmass2*xnxa +txn*rxa(i) -rxn(i)*txa) )
1360 am3pol(i)=xm3pol(i)*const3
1363 & (gv+ga)**2*tmass*xnxa*qp(i)
1364 & -(gv-ga)**2*tmass*qpxn*xa(i)
1365 & +(gv**2-ga**2)*amnuta*txa*qp(i)
1366 & -(gv**2-ga**2)*amnuta*tqp*xa(i) )*
1368 5 hv(i)=(bornpl(i)+am3pol(i))/(born+am3)-1.d0
1370 IF (thb/born.LT.0.1d0)
THEN
1371 print *,
'ERROR IN THB, THB/BORN=',thb/born
1378 SUBROUTINE dexpi(MODE,ISGN,POL,PPI,PNU)
1385 REAL POL(4),HV(4),PNU(4),PPI(4),RN(1)
1389 CALL dadmpi(-1,isgn,hv,ppi,pnu)
1392 ELSEIF(mode.EQ. 0)
THEN
1395 CALL dadmpi( 0,isgn,hv,ppi,pnu)
1396 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1399 IF(rn(1).GT.wt)
GOTO 300
1401 ELSEIF(mode.EQ. 1)
THEN
1403 CALL dadmpi( 1,isgn,hv,ppi,pnu)
1409 SUBROUTINE dadmpi(MODE,ISGN,HV,PPI,PNU)
1411 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1412 * ,ampiz,ampi,amro,gamro,ama1,gama1
1413 * ,amk,amkz,amkst,gamkst
1415 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1416 * ,ampiz,ampi,amro,gamro,ama1,gama1
1417 * ,amk,amkz,amkst,gamkst
1418 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1419 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1420 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
1421 real*4 gampmc ,gamper
1422 COMMON / inout / inut,iout
1423 REAL PPI(4),PNU(4),HV(4)
1424 DATA pi /3.141592653589793238462643/
1429 ELSEIF(mode.EQ. 0)
THEN
1432 epi= (amtau**2+ampi**2-amnuta**2)/(2*amtau)
1433 enu= (amtau**2-ampi**2+amnuta**2)/(2*amtau)
1434 xpi= sqrt(epi**2-ampi**2)
1436 CALL sphera(xpi,ppi)
1444 qxn=ppi(4)*pnu(4)-ppi(1)*pnu(1)-ppi(2)*pnu(2)-ppi(3)*pnu(3)
1445 brak=(gv**2+ga**2)*(2*pxq*qxn-ampi**2*pxn)
1446 & +(gv**2-ga**2)*amtau*amnuta*ampi**2
1448 40 hv(i)=-isgn*2*ga*gv*amtau*(2*ppi(i)*qxn-pnu(i)*ampi**2)/brak
1451 ELSEIF(mode.EQ. 1)
THEN
1453 IF(nevtot.EQ.0)
RETURN
1459 gamm=(gfermi*fpi)**2/(16.*pi)*amtau**3*
1461 $ sqrt((amtau**2-ampi**2-amnuta**2)**2
1462 $ -4*ampi**2*amnuta**2 )/amtau**2
1465 WRITE(iout, 7010) nevtot,gamm,rat,error
1472 7010
FORMAT(///1x,15(5h*****)
1473 $ /,
' *', 25x,
'******** DADMPI FINAL REPORT ******** ',9x,1h*
1474 $ /,
' *',i20 ,5x,
'NEVTOT = NO. OF PI DECAYS TOTAL ',9x,1h*
1475 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( PI DECAY) IN GEV UNITS ',9x,1h*
1476 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1477 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9x,1h*
1478 $ /,1x,15(5h*****)/)
1480 SUBROUTINE dexro(MODE,ISGN,POL,PNU,PRO,PIC,PIZ)
1489 COMMON / inout / inut,iout
1490 REAL POL(4),HV(4),PRO(4),PNU(4),PIC(4),PIZ(4),RN(1)
1496 CALL dadmro( -1,isgn,hv,pnu,pro,pic,piz)
1500 ELSEIF(mode.EQ. 0)
THEN
1503 IF(iwarm.EQ.0)
GOTO 902
1504 CALL dadmro( 0,isgn,hv,pnu,pro,pic,piz)
1505 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1510 IF(rn(1).GT.wt)
GOTO 300
1512 ELSEIF(mode.EQ. 1)
THEN
1514 CALL dadmro( 1,isgn,hv,pnu,pro,pic,piz)
1520 902
WRITE(iout, 9020)
1521 9020
FORMAT(
' ----- DEXRO: LACK OF INITIALISATION')
1524 SUBROUTINE dadmro(MODE,ISGN,HHV,PNU,PRO,PIC,PIZ)
1526 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1527 * ,ampiz,ampi,amro,gamro,ama1,gama1
1528 * ,amk,amkz,amkst,gamkst
1530 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1531 * ,ampiz,ampi,amro,gamro,ama1,gama1
1532 * ,amk,amkz,amkst,gamkst
1533 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1534 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1535 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
1536 real*4 gampmc ,gamper
1537 COMMON / inout / inut,iout
1539 REAL HV(4),PRO(4),PNU(4),PIC(4),PIZ(4)
1540 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4)
1543 DATA pi /3.141592653589793238462643/
1556 CALL dphsro(wt,hv,pdum1,pdum2,pdum3,pdum4)
1557 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1562 ELSEIF(mode.EQ. 0)
THEN
1565 IF(iwarm.EQ.0)
GOTO 902
1566 CALL dphsro(wt,hv,pnu,pro,pic,piz)
1573 IF(wt.GT.wtmax) nevovr=nevovr+1
1574 IF(rn*wtmax.GT.wt)
GOTO 300
1576 costhe=-1.+2.*rrr(2)
1579 CALL rotor2(thet,pnu,pnu)
1580 CALL rotor3( phi,pnu,pnu)
1581 CALL rotor2(thet,pro,pro)
1582 CALL rotor3( phi,pro,pro)
1583 CALL rotor2(thet,pic,pic)
1584 CALL rotor3( phi,pic,pic)
1585 CALL rotor2(thet,piz,piz)
1586 CALL rotor3( phi,piz,piz)
1587 CALL rotor2(thet,hv,hv)
1588 CALL rotor3( phi,hv,hv)
1590 44 hhv(i)=-isgn*hv(i)
1593 ELSEIF(mode.EQ. 1)
THEN
1595 IF(nevraw.EQ.0)
RETURN
1596 pargam=swt/float(nevraw+1)
1598 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1600 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1608 7003
FORMAT(///1x,15(5h*****)
1609 $ /,
' *', 25x,
'******** DADMRO INITIALISATION ********',9x,1h*
1610 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
1611 $ /,1x,15(5h*****)/)
1612 7010
FORMAT(///1x,15(5h*****)
1613 $ /,
' *', 25x,
'******** DADMRO FINAL REPORT ******** ',9x,1h*
1614 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF RHO DECAYS TOTAL ',9x,1h*
1615 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF RHO DECS. ACCEPTED ',9x,1h*
1616 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
1617 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (RHO DECAY) IN GEV UNITS ',9x,1h*
1618 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1619 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
1620 $ /,1x,15(5h*****)/)
1621 902
WRITE(iout, 9020)
1622 9020
FORMAT(
' ----- DADMRO: LACK OF INITIALISATION')
1625 SUBROUTINE dphsro(DGAMT,HV,PN,PR,PIC,PIZ)
1630 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1631 * ,ampiz,ampi,amro,gamro,ama1,gama1
1632 * ,amk,amkz,amkst,gamkst
1634 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1635 * ,ampiz,ampi,amro,gamro,ama1,gama1
1636 * ,amk,amkz,amkst,gamkst
1637 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1638 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1639 REAL HV(4),PT(4),PN(4),PR(4),PIC(4),PIZ(4),QQ(4),RR1(1)
1640 DATA pi /3.141592653589793238462643/
1644 phspac=1./2**11/pi**5
1651 ams1=(ampi+ampiz)**2
1652 ams2=(amtau-amnuta)**2
1660 alp1=atan((ams1-amro**2)/amro/gamro)
1661 alp2=atan((ams2-amro**2)/amro/gamro)
1665 alp=alp1+rr1(1)*(alp2-alp1)
1666 amx2=amro**2+amro*gamro*tan(alp)
1668 IF(amx.LT.2.*ampi)
GO TO 100
1670 phspac=phspac*((amx2-amro**2)**2+(amro*gamro)**2)/(amro*gamro)
1671 phspac=phspac*(alp2-alp1)
1676 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
1677 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
1681 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
1683 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
1686 enq1=(amx2+ampi**2-ampiz**2)/(2.*amx)
1687 enq2=(amx2-ampi**2+ampiz**2)/(2.*amx)
1688 pppi=sqrt((enq1-ampi)*(enq1+ampi))
1689 phspac=phspac*(4*pi)*(2*pppi/amx)
1691 CALL sphera(pppi,pic)
1697 exe=(pr(4)+pr(3))/amx
1699 CALL bostr3(exe,pic,pic)
1700 CALL bostr3(exe,piz,piz)
1702 30 qq(i)=pic(i)-piz(i)
1705 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
1707 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
1708 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
1709 & +(gv**2-ga**2)*amtau*amnuta*qq2
1710 amplit=(gfermi*ccabib)**2*brak*2*fpirho(amx)
1711 dgamt=1/(2.*amtau)*amplit*phspac
1713 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
1716 SUBROUTINE dexaa(MODE,ISGN,POL,PNU,PAA,PIM1,PIM2,PIPL,JAA)
1727 COMMON / inout / inut,iout
1728 REAL POL(4),HV(4),PAA(4),PNU(4),PIM1(4),PIM2(4),PIPL(4),RN(1)
1734 CALL dadmaa( -1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1737 ELSEIF(mode.EQ. 0)
THEN
1740 IF(iwarm.EQ.0)
GOTO 902
1741 CALL dadmaa( 0,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1742 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1745 IF(rn(1).GT.wt)
GOTO 300
1747 ELSEIF(mode.EQ. 1)
THEN
1749 CALL dadmaa( 1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
1754 902
WRITE(iout, 9020)
1755 9020
FORMAT(
' ----- DEXAA: LACK OF INITIALISATION')
1758 SUBROUTINE dadmaa(MODE,ISGN,HHV,PNU,PAA,PIM1,PIM2,PIPL,JAA)
1762 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1763 * ,ampiz,ampi,amro,gamro,ama1,gama1
1764 * ,amk,amkz,amkst,gamkst
1766 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1767 * ,ampiz,ampi,amro,gamro,ama1,gama1
1768 * ,amk,amkz,amkst,gamkst
1769 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1770 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1771 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
1772 real*4 gampmc ,gamper
1773 COMMON / inout / inut,iout
1775 REAL HV(4),PAA(4),PNU(4),PIM1(4),PIM2(4),PIPL(4)
1776 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
1779 DATA pi /3.141592653589793238462643/
1792 CALL dphsaa(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5,jaa)
1793 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1797 ELSEIF(mode.EQ. 0)
THEN
1800 IF(iwarm.EQ.0)
GOTO 902
1801 CALL dphsaa(wt,hv,pnu,paa,pim1,pim2,pipl,jaa)
1808 sswt=sswt+dble(wt)**2
1813 IF(wt.GT.wtmax) nevovr=nevovr+1
1814 IF(rn*wtmax.GT.wt)
GOTO 300
1816 costhe=-1.+2.*rrr(2)
1819 CALL rotpol(thet,phi,pnu)
1820 CALL rotpol(thet,phi,paa)
1821 CALL rotpol(thet,phi,pim1)
1822 CALL rotpol(thet,phi,pim2)
1823 CALL rotpol(thet,phi,pipl)
1824 CALL rotpol(thet,phi,hv)
1826 44 hhv(i)=-isgn*hv(i)
1829 ELSEIF(mode.EQ. 1)
THEN
1831 IF(nevraw.EQ.0)
RETURN
1832 pargam=swt/float(nevraw+1)
1834 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1836 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1844 7003
FORMAT(///1x,15(5h*****)
1845 $ /,
' *', 25x,
'******** DADMAA INITIALISATION ********',9x,1h*
1846 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
1847 $ /,1x,15(5h*****)/)
1848 7010
FORMAT(///1x,15(5h*****)
1849 $ /,
' *', 25x,
'******** DADMAA FINAL REPORT ******** ',9x,1h*
1850 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF A1 DECAYS TOTAL ',9x,1h*
1851 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF A1 DECS. ACCEPTED ',9x,1h*
1852 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
1853 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (A1 DECAY) IN GEV UNITS ',9x,1h*
1854 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1855 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
1856 $ /,1x,15(5h*****)/)
1857 902
WRITE(iout, 9020)
1858 9020
FORMAT(
' ----- DADMAA: LACK OF INITIALISATION')
1861 SUBROUTINE dphsaa(DGAMT,HV,PN,PAA,PIM1,PIM2,PIPL,JAA)
1866 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1867 * ,ampiz,ampi,amro,gamro,ama1,gama1
1868 * ,amk,amkz,amkst,gamkst
1870 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1871 * ,ampiz,ampi,amro,gamro,ama1,gama1
1872 * ,amk,amkz,amkst,gamkst
1873 COMMON / taukle / bra1,brk0,brk0b,brks
1874 real*4 bra1,brk0,brk0b,brks
1875 REAL HV(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
1885 IF (rmod.LT.bra1)
THEN
1897 $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
1899 SUBROUTINE dexkk(MODE,ISGN,POL,PKK,PNU)
1906 REAL POL(4),HV(4),PNU(4),PKK(4),RN(1)
1910 CALL dadmkk(-1,isgn,hv,pkk,pnu)
1913 ELSEIF(mode.EQ. 0)
THEN
1916 CALL dadmkk( 0,isgn,hv,pkk,pnu)
1917 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1920 IF(rn(1).GT.wt)
GOTO 300
1922 ELSEIF(mode.EQ. 1)
THEN
1924 CALL dadmkk( 1,isgn,hv,pkk,pnu)
1930 SUBROUTINE dadmkk(MODE,ISGN,HV,PKK,PNU)
1934 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1935 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1937 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1938 * ,ampiz,ampi,amro,gamro,ama1,gama1
1939 * ,amk,amkz,amkst,gamkst
1941 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1942 * ,ampiz,ampi,amro,gamro,ama1,gama1
1943 * ,amk,amkz,amkst,gamkst
1946 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
1947 real*4 gampmc ,gamper
1948 COMMON / inout / inut,iout
1949 REAL PKK(4),PNU(4),HV(4)
1950 DATA pi /3.141592653589793238462643/
1955 ELSEIF(mode.EQ. 0)
THEN
1958 ekk= (amtau**2+amk**2-amnuta**2)/(2*amtau)
1959 enu= (amtau**2-amk**2+amnuta**2)/(2*amtau)
1960 xkk= sqrt(ekk**2-amk**2)
1962 CALL sphera(xkk,pkk)
1970 qxn=pkk(4)*pnu(4)-pkk(1)*pnu(1)-pkk(2)*pnu(2)-pkk(3)*pnu(3)
1971 brak=(gv**2+ga**2)*(2*pxq*qxn-amk**2*pxn)
1972 & +(gv**2-ga**2)*amtau*amnuta*amk**2
1974 40 hv(i)=-isgn*2*ga*gv*amtau*(2*pkk(i)*qxn-pnu(i)*amk**2)/brak
1977 ELSEIF(mode.EQ. 1)
THEN
1979 IF(nevtot.EQ.0)
RETURN
1986 gamm=(gfermi*fkk)**2/(16.*pi)*amtau**3*
1988 $ sqrt((amtau**2-amk**2-amnuta**2)**2
1989 $ -4*amk**2*amnuta**2 )/amtau**2
1994 WRITE(iout, 7010) nevtot,gamm,rat,error
2001 7010
FORMAT(///1x,15(5h*****)
2002 $ /,
' *', 25x,
'******** DADMKK FINAL REPORT ********',9x,1h*
2003 $ /,
' *',i20 ,5x,
'NEVTOT = NO. OF K DECAYS TOTAL ',9x,1h*,
2004 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( K DECAY) IN GEV UNITS ',9x,1h*,
2005 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2006 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9x,1h*
2007 $ /,1x,15(5h*****)/)
2009 SUBROUTINE dexks(MODE,ISGN,POL,PNU,PKS,PKK,PPI,JKST)
2021 COMMON / inout / inut,iout
2022 REAL POL(4),HV(4),PKS(4),PNU(4),PKK(4),PPI(4),RN(1)
2029 CALL dadmks( -1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2033 ELSEIF(mode.EQ. 0)
THEN
2036 IF(iwarm.EQ.0)
GOTO 902
2037 CALL dadmks( 0,isgn,hv,pnu,pks,pkk,ppi,jkst)
2038 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
2043 IF(rn(1).GT.wt)
GOTO 300
2045 ELSEIF(mode.EQ. 1)
THEN
2047 CALL dadmks( 1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2053 902
WRITE(iout, 9020)
2054 9020
FORMAT(
' ----- DEXKS: LACK OF INITIALISATION')
2057 SUBROUTINE dadmks(MODE,ISGN,HHV,PNU,PKS,PKK,PPI,JKST)
2059 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2060 * ,ampiz,ampi,amro,gamro,ama1,gama1
2061 * ,amk,amkz,amkst,gamkst
2063 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2064 * ,ampiz,ampi,amro,gamro,ama1,gama1
2065 * ,amk,amkz,amkst,gamkst
2066 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2067 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2068 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
2069 real*4 gampmc ,gamper
2070 COMMON / taukle / bra1,brk0,brk0b,brks
2071 real*4 bra1,brk0,brk0b,brks
2072 COMMON / inout / inut,iout
2074 REAL HV(4),PKS(4),PNU(4),PKK(4),PPI(4)
2075 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4)
2076 real*4 rrr(3),rmod(1)
2078 DATA pi /3.141592653589793238462643/
2093 CALL dphsks(wt,hv,pdum1,pdum2,pdum3,pdum4,jkst)
2094 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
2099 ELSEIF(mode.EQ. 0)
THEN
2101 IF(iwarm.EQ.0)
GOTO 902
2107 IF(rmod(1).LT.dec1)
THEN
2112 CALL dphsks(wt,hv,pnu,pks,pkk,ppi,jkst)
2115 IF(wt.GT.wtmax) nevovr=nevovr+1
2119 IF(rn*wtmax.GT.wt)
GOTO 400
2121 costhe=-1.+2.*rrr(2)
2124 CALL rotor2(thet,pnu,pnu)
2125 CALL rotor3( phi,pnu,pnu)
2126 CALL rotor2(thet,pks,pks)
2127 CALL rotor3( phi,pks,pks)
2128 CALL rotor2(thet,pkk,pkk)
2129 CALL rotor3(phi,pkk,pkk)
2130 CALL rotor2(thet,ppi,ppi)
2131 CALL rotor3( phi,ppi,ppi)
2132 CALL rotor2(thet,hv,hv)
2133 CALL rotor3( phi,hv,hv)
2135 44 hhv(i)=-isgn*hv(i)
2138 ELSEIF(mode.EQ. 1)
THEN
2140 IF(nevraw.EQ.0)
RETURN
2141 pargam=swt/float(nevraw+1)
2143 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
2145 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
2153 7003
FORMAT(///1x,15(5h*****)
2154 $ /,
' *', 25x,
'******** DADMKS INITIALISATION ********',9x,1h*
2155 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
2156 $ /,1x,15(5h*****)/)
2157 7010
FORMAT(///1x,15(5h*****)
2158 $ /,
' *', 25x,
'******** DADMKS FINAL REPORT ********',9x,1h*
2159 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF K* DECAYS TOTAL ',9x,1h*,
2160 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF K* DECS. ACCEPTED ',9x,1h*,
2161 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
2162 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (K* DECAY) IN GEV UNITS ',9x,1h*,
2163 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2164 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
2165 $ /,1x,15(5h*****)/)
2166 902
WRITE(iout, 9020)
2167 9020
FORMAT(
' ----- DADMKS: LACK OF INITIALISATION')
2170 SUBROUTINE dphsks(DGAMT,HV,PN,PKS,PKK,PPI,JKST)
2178 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2179 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2181 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2182 * ,ampiz,ampi,amro,gamro,ama1,gama1
2183 * ,amk,amkz,amkst,gamkst
2185 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2186 * ,ampiz,ampi,amro,gamro,ama1,gama1
2187 * ,amk,amkz,amkst,gamkst
2190 REAL HV(4),PT(4),PN(4),PKS(4),PKK(4),PPI(4),QQ(4),RR1(1)
2194 DATA pi /3.141592653589793238462643/
2198 phspac=1./2**11/pi**5
2210 ams2=(amtau-amnuta)**2
2216 alp1=atan((ams1-amkst**2)/amkst/gamkst)
2217 alp2=atan((ams2-amkst**2)/amkst/gamkst)
2218 alp=alp1+rr1(1)*(alp2-alp1)
2219 amx2=amkst**2+amkst*gamkst*tan(alp)
2221 phspac=phspac*((amx2-amkst**2)**2+(amkst*gamkst)**2)
2223 phspac=phspac*(alp2-alp1)
2228 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2229 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2234 pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2236 phspac=phspac*(4*pi)*(2*pks(3)/amtau)
2239 enpi=( amx**2+ampi**2-amkz**2 ) / ( 2*amx )
2240 pppi=sqrt((enpi-ampi)*(enpi+ampi))
2241 phspac=phspac*(4*pi)*(2*pppi/amx)
2243 CALL sphera(pppi,ppi)
2248 pkk(4)=( amx**2+amkz**2-ampi**2 ) / ( 2*amx )
2249 exe=(pks(4)+pks(3))/amx
2251 CALL bostr3(exe,ppi,ppi)
2252 CALL bostr3(exe,pkk,pkk)
2254 30 qq(i)=ppi(i)-pkk(i)
2256 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
2257 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
2259 31 qq(i)=qq(i)-pks(i)*qqpks/pksd
2262 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
2264 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
2265 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
2266 & +(gv**2-ga**2)*amtau*amnuta*qq2
2269 fks=cabs(bwigs(amx2,amkst,gamkst))**2
2271 amplit=(gfermi*scabib)**2*brak*2*fks
2272 dgamt=1/(2.*amtau)*amplit*phspac
2274 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
2277 ELSEIF(jkst.EQ.20)
THEN
2281 ams2=(amtau-amnuta)**2
2289 alp1=atan((ams1-amkst**2)/amkst/gamkst)
2290 alp2=atan((ams2-amkst**2)/amkst/gamkst)
2291 alp=alp1+rr1(1)*(alp2-alp1)
2292 amx2=amkst**2+amkst*gamkst*tan(alp)
2294 phspac=phspac*((amx2-amkst**2)**2+(amkst*gamkst)**2)
2296 phspac=phspac*(alp2-alp1)
2301 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2302 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2306 pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2308 phspac=phspac*(4*pi)*(2*pks(3)/amtau)
2311 enpi=( amx**2+ampiz**2-amk**2 ) / ( 2*amx )
2312 pppi=sqrt((enpi-ampiz)*(enpi+ampiz))
2313 phspac=phspac*(4*pi)*(2*pppi/amx)
2315 CALL sphera(pppi,ppi)
2320 pkk(4)=( amx**2+amk**2-ampiz**2 ) / ( 2*amx )
2321 exe=(pks(4)+pks(3))/amx
2323 CALL bostr3(exe,ppi,ppi)
2324 CALL bostr3(exe,pkk,pkk)
2326 60 qq(i)=pkk(i)-ppi(i)
2328 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
2329 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
2331 61 qq(i)=qq(i)-pks(i)*qqpks/pksd
2334 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
2336 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
2337 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
2338 & +(gv**2-ga**2)*amtau*amnuta*qq2
2341 fks=cabs(bwigs(amx2,amkst,gamkst))**2
2343 amplit=(gfermi*scabib)**2*brak*2*fks
2344 dgamt=1/(2.*amtau)*amplit*phspac
2346 70 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
2354 SUBROUTINE dphnpi(DGAMT,HVX,PNX,PRX,PPIX,JNPI)
2360 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2361 * ,ampiz,ampi,amro,gamro,ama1,gama1
2362 * ,amk,amkz,amkst,gamkst
2364 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2365 * ,ampiz,ampi,amro,gamro,ama1,gama1
2366 * ,amk,amkz,amkst,gamkst
2367 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2368 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2369 parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
2371 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2373 CHARACTER NAMES(NMODE)*31
2378 real*8 pn(4),pr(4),ppi(4,9),hv(4)
2379 real*4 pnx(4),prx(4),ppix(4,9),hvx(4)
2380 real*8 pv(5,9),pt(4),ue(3),be(3)
2381 real*8 pawt,amx,ams1,ams2,pa,phs,phsmax,pmin,pmax
2385 real*8 gam,bep,phi,a,b,c
2387 real*4 rrr(9),rrx(2),rn(1),rr2(1)
2389 DATA pi /3.141592653589793238462643/
2390 DATA wetmax /500*1d-15/
2395 $ sqrt(max(0.d0,(a**2-(b+c)**2)*(a**2-(b-c)**2)))/(2.d0*a)
2398 ampik(i,j)=dcdmas(idffin(i,j))
2402 IF ((jnpi.LE.0).OR.jnpi.GT.100)
THEN
2403 WRITE(6,*)
'JNPI OUTSIDE RANGE DEFINED BY WETMAX; JNPI=',jnpi
2420 phspac = 1./2.**5 /pi**2
2422 4 ps =ps+ampik(i,jnpi)
2427 ams2=(amtau-amnuta)**2
2431 amx2=ams1+ rr2(1)*(ams2-ams1)
2435 phspac=phspac * (ams2-ams1)
2440 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx2)
2441 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
2445 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx2)
2447 phspac=phspac * (4.*pi) * (2.*pr(3)/amtau)
2454 qxn=pr(4)*pn(4)-pr(1)*pn(1)-pr(2)*pn(2)-pr(3)*pn(3)
2459 brak=2*(gv**2+ga**2)*(2*pxq*qxn+amx2*pxn)
2460 & -6*(gv**2-ga**2)*amtau*amnuta*amx2
2466 amplit=ccabib**2*gfermi**2/2.* brak*amx2*sigee(amx2,jn)
2467 dgamt=1./(2.*amtau)*amplit*phspac
2476 pv(5,nd)=ampik(nd,jnpi)
2478 pmax=amw-ps+ampik(nd,jnpi)
2481 pmax=pmax+ampik(il,jnpi)
2482 pmin=pmin+ampik(il+1,jnpi)
2484 220 phsmax=phsmax*pawt(pmax,pmin,ampik(il,jnpi))/pmax
2491 223 ams1=ams1+ampik(jl,jnpi)
2493 amx =(amx-ampik(il,jnpi))
2495 phsmax=phsmax * (ams2-ams1)
2502 phspac = 1./2.**(6*nd-7) /pi**(3*nd-4)
2504 CALL ranmar(rrr,nd-2)
2508 231 ams1=ams1+ampik(jl,jnpi)
2510 ams2=(amx-ampik(il,jnpi))**2
2512 amx2=ams1+ rr1*(ams2-ams1)
2515 phspac=phspac * (ams2-ams1)
2517 phs=phs* (ams2-ams1)
2518 pa=pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2519 phs =phs *pa/pv(5,il)
2521 pa=pawt(pv(5,nd-1),ampik(nd-1,jnpi),ampik(nd,jnpi))
2522 phs =phs *pa/pv(5,nd-1)
2524 wetmax(jnpi)=1.2d0*max(wetmax(jnpi)/1.2d0,phs/phsmax)
2525 IF (ncont.EQ.500 000)
THEN
2528 xnpi=xnpi+ampik(kk,jnpi)
2530 WRITE(6,*)
'ROUNDING INSTABILITY IN DPHNPI ?'
2531 WRITE(6,*)
'AMW=',amw,
'XNPI=',xnpi
2532 WRITE(6,*)
'IF =AMW= IS NEARLY EQUAL =XNPI= THAT IS IT'
2533 WRITE(6,*)
'PHS=',phs,
'PHSMAX=',phsmax
2536 IF(rn(1)*phsmax*wetmax(jnpi).GT.phs)
GO TO 100
2539 280
DO 300 il=1,nd-1
2540 pa=pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2545 ue(1)=sqrt(1.d0-ue(3)**2)*cos(phi)
2546 ue(2)=sqrt(1.d0-ue(3)**2)*sin(phi)
2550 290 pv(j,il+1)=-pa*ue(j)
2551 ppi(4,il)=sqrt(pa**2+ampik(il,jnpi)**2)
2552 pv(4,il+1)=sqrt(pa**2+pv(5,il+1)**2)
2553 phspac=phspac *(4.*pi)*(2.*pa/pv(5,il))
2557 310 ppi(j,nd)=pv(j,nd)
2560 320 be(j)=pv(j,il)/pv(4,il)
2561 gam=pv(4,il)/pv(5,il)
2563 bep=be(1)*ppi(1,i)+be(2)*ppi(2,i)+be(3)*ppi(3,i)
2566 330 ppi(j,i)=ppi(j,i)+gam*(gam*bep/(1.d0+gam)+ppi(4,i))*be(j)
2568 ppi(4,i)=gam*(ppi(4,i)+bep)
2587 FUNCTION sigee(Q2,JNP)
2606 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2607 * ,ampiz,ampi,amro,gamro,ama1,gama1
2608 * ,amk,amkz,amkst,gamkst
2610 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2611 * ,ampiz,ampi,amro,gamro,ama1,gama1
2612 * ,amk,amkz,amkst,gamkst
2616 1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
2617 2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
2618 3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
2619 4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
2622 7 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25,
2625 DATA pi /3.141592653589793238462643/
2630 IF(jnpi.GT.6) jnpi=6
2643 datsig(i,2) = datsig(i,2)/2.
2644 datsig(i,1) = datsig(i,1) + datsig(i,2)
2650 IF(t . gt. s-ampi )
GO TO 201
2652 fact=(t2/s2)**2*sqrt((s2-t2-ampi2)**2-4.*t2*ampi2)/s2 *2.*t*.05
2653 fact = fact * (datsig(j,1)+datsig(j+1,1))
2654 200 datsig(i,3) = datsig(i,3) + fact
2655 201 datsig(i,3) = datsig(i,3) /(2*pi*fpi)**2
2656 datsig(i,4) = datsig(i,3)
2657 datsig(i,6) = datsig(i,5)
2660 1000
FORMAT(///1x,
' EE SIGMA USED IN MULTIPI DECAYS'/
2666 sigee=datsig(1,jnpi)+
2667 & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
2668 ELSEIF(q.LT.1.8)
THEN
2671 IF(q.LT.qmax)
GO TO 2
2674 2 sigee=datsig(i,jnpi)+
2675 & (datsig(i+1,jnpi)-datsig(i,jnpi)) * (q-qmin)/.05
2676 ELSEIF(q.GT.1.8)
THEN
2677 sigee=datsig(17,jnpi)+
2678 & (datsig(17,jnpi)-datsig(16,jnpi)) * (q-1.8)/.05
2680 IF(sigee.LT..0) sigee=0.
2682 sigee = sigee/(6.*pi**2*sig0)
2687 FUNCTION sigold(Q2,JNPI)
2706 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2707 * ,ampiz,ampi,amro,gamro,ama1,gama1
2708 * ,amk,amkz,amkst,gamkst
2710 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2711 * ,ampiz,ampi,amro,gamro,ama1,gama1
2712 * ,amk,amkz,amkst,gamkst
2716 1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
2717 2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
2718 3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
2719 4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
2721 6 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25/
2723 DATA pi /3.141592653589793238462643/
2731 datsig(i,2) = datsig(i,2)/2.
2732 datsig(i,1) = datsig(i,1) + datsig(i,2)
2738 IF(t . gt. s-ampi )
GO TO 201
2740 fact=(t2/s2)**2*sqrt((s2-t2-ampi2)**2-4.*t2*ampi2)/s2 *2.*t*.05
2741 fact = fact * (datsig(j,1)+datsig(j+1,1))
2742 200 datsig(i,3) = datsig(i,3) + fact
2743 201 datsig(i,3) = datsig(i,3) /(2*pi*fpi)**2
2746 1000
FORMAT(///1x,
' EE SIGMA USED IN MULTIPI DECAYS'/
2752 sigee=datsig(1,jnpi)+
2753 & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
2754 ELSEIF(q.LT.1.8)
THEN
2757 IF(q.LT.qmax)
GO TO 2
2760 2 sigee=datsig(i,jnpi)+
2761 & (datsig(i+1,jnpi)-datsig(i,jnpi)) * (q-qmin)/.05
2762 ELSEIF(q.GT.1.8)
THEN
2763 sigee=datsig(17,jnpi)+
2764 & (datsig(17,jnpi)-datsig(16,jnpi)) * (q-1.8)/.05
2766 IF(sigee.LT..0) sigee=0.
2768 sigee = sigee/(6.*pi**2*sig0)
2773 SUBROUTINE dphspk(DGAMT,HV,PN,PAA,PNPI,JAA)
2778 parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
2780 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2783 CHARACTER NAMES(NMODE)*31
2785 REAL HV(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4),PNPI(4,9)
2793 amp1=dcdmas(idffin(1,jaa+nm4+nm5+nm6))
2794 amp2=dcdmas(idffin(2,jaa+nm4+nm5+nm6))
2795 amp3=dcdmas(idffin(3,jaa+nm4+nm5+nm6))
2797 $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
2809 $ dphtre(dgamt,hv,pn,paa,pim1,ampa,pim2,ampb,pipl,amp3,keyt,mnum)
2826 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2827 * ,ampiz,ampi,amro,gamro,ama1,gama1
2828 * ,amk,amkz,amkst,gamkst
2830 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2831 * ,ampiz,ampi,amro,gamro,ama1,gama1
2832 * ,amk,amkz,amkst,gamkst
2833 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2834 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2835 REAL HV(4),PT(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
2838 DATA pi /3.141592653589793238462643/
2840 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
2845 phspac=1./2**17/pi**8
2855 CALL choice(mnum,rr,ichan,prob1,prob2,prob3,
2856 $ amrx,gamrx,amra,gamra,amrb,gamrb)
2857 IF (ichan.EQ.1)
THEN
2860 ELSEIF (ichan.EQ.2)
THEN
2869 ams1=(amp1+amp2+amp3)**2
2870 ams2=(amtau-amnuta)**2
2874 alp1=atan((ams1-amrx**2)/amrx/gamrx)
2875 alp2=atan((ams2-amrx**2)/amrx/gamrx)
2876 alp=alp1+rr1*(alp2-alp1)
2877 am3sq =amrx**2+amrx*gamrx*tan(alp)
2879 phspac=phspac*((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
2880 phspac=phspac*(alp2-alp1)
2885 IF (ichan.LE.2)
THEN
2889 alp1=atan((ams1-amra**2)/amra/gamra)
2890 alp2=atan((ams2-amra**2)/amra/gamra)
2891 alp=alp1+rr2*(alp2-alp1)
2892 am2sq =amra**2+amra*gamra*tan(alp)
2902 am2sq=ams1+ rr2*(ams2-ams1)
2909 enq1=(am2sq-amp2**2+amp3**2)/(2*am2)
2910 enq2=(am2sq+amp2**2-amp3**2)/(2*am2)
2911 ppi= enq1**2-amp3**2
2912 pppi=sqrt(abs(enq1**2-amp3**2))
2914 phf1=(4*pi)*(2*pppi/am2)
2918 CALL sphera(pppi,pipl)
2932 pr(4)=1./(2*am3)*(am3**2+am2**2-amp1**2)
2933 pr(3)= sqrt(abs(pr(4)**2-am2**2))
2934 ppi = pr(4)**2-am2**2
2938 pim2(4)=1./(2*am3)*(am3**2-am2**2+amp1**2)
2940 phf2=(4*pi)*(2*pr(3)/am3)
2944 exe=(pr(4)+pr(3))/am2
2945 CALL bostr3(exe,pipl,pipl)
2946 CALL bostr3(exe,pim1,pim1)
2952 thet =acos(-1.+2*rr3)
2954 CALL rotpol(thet,phi,pipl)
2955 CALL rotpol(thet,phi,pim1)
2956 CALL rotpol(thet,phi,pim2)
2957 CALL rotpol(thet,phi,pr)
2963 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am3**2)
2964 paa(3)= sqrt(abs(paa(4)**2-am3**2))
2965 ppi = paa(4)**2-am3**2
2966 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
2970 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am3**2)
2976 alp1=atan((ams1-amra**2)/amra/gamra)
2977 alp2=atan((ams2-amra**2)/amra/gamra)
2978 xpro = (pim1(3)+pipl(3))**2
2979 $ +(pim1(2)+pipl(2))**2+(pim1(1)+pipl(1))**2
2980 am2sq=-xpro+(pim1(4)+pipl(4))**2
2982 ff1 = ((am2sq-amra**2)**2+(amra*gamra)**2)/(amra*gamra)
2983 ff1 =ff1 *(alp2-alp1)
2985 gg1 = (4*pi)*(xlam(am2sq,amp2**2,amp3**2)/am2sq)
2987 gg1 =gg1 *(4*pi)*sqrt(4*xpro/am3sq)
2988 xjaje=gg1*(ams2-ams1)
2992 alp1=atan((ams1-amrb**2)/amrb/gamrb)
2993 alp2=atan((ams2-amrb**2)/amrb/gamrb)
2994 xpro = (pim2(3)+pipl(3))**2
2995 $ +(pim2(2)+pipl(2))**2+(pim2(1)+pipl(1))**2
2996 am2sq=-xpro+(pim2(4)+pipl(4))**2
2997 ff2 = ((am2sq-amrb**2)**2+(amrb*gamrb)**2)/(amrb*gamrb)
2998 ff2 =ff2 *(alp2-alp1)
2999 gg2 = (4*pi)*(xlam(am2sq,amp1**2,amp3**2)/am2sq)
3000 gg2 =gg2 *(4*pi)*sqrt(4*xpro/am3sq)
3001 xjadw=gg2*(ams2-ams1)
3008 IF (ichan.EQ.2)
THEN
3013 IF (xjac1.NE.0.0) a1=prob1/xjac1
3014 IF (xjac2.NE.0.0) a2=prob2/xjac2
3015 IF (xjac3.NE.0.0) a3=prob3/xjac3
3017 IF (a1+a2+a3.NE.0.0)
THEN
3018 phspac=phspac/(a1+a2+a3)
3030 exe=(paa(4)+paa(3))/am3
3031 CALL bostr3(exe,pipl,pipl)
3032 CALL bostr3(exe,pim1,pim1)
3033 CALL bostr3(exe,pim2,pim2)
3034 CALL bostr3(exe,pr,pr)
3037 CALL dampog(pt,pn,pim1,pim2,pipl,amplit,hv)
3042 CALL damppk(mnum,pt,pn,pim1,pim2,pipl,amplit,hv)
3047 IF (keyt.EQ.1.OR.keyt.EQ.2)
THEN
3054 dgamt=1/(2.*amtau)*amplit*phspac
3056 SUBROUTINE dampaa(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3066 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3067 * ,ampiz,ampi,amro,gamro,ama1,gama1
3068 * ,amk,amkz,amkst,gamkst
3070 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3071 * ,ampiz,ampi,amro,gamro,ama1,gama1
3072 * ,amk,amkz,amkst,gamkst
3073 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3074 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3075 COMMON /testa1/ keya1
3076 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4)
3077 REAL PAA(4),VEC1(4),VEC2(4)
3078 REAL PIVEC(4),PIAKS(4),HVM(4)
3079 COMPLEX BWIGN,HADCUR(4),FPIK
3086 bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
3090 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3092 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3093 xmro1 =sqrt(abs((pipl(4)+pim1(4))**2-(pipl(1)+pim1(1))**2
3094 $ -(pipl(2)+pim1(2))**2-(pipl(3)+pim1(3))**2))
3095 xmro2 =sqrt(abs((pipl(4)+pim2(4))**2-(pipl(1)+pim2(1))**2
3096 $ -(pipl(2)+pim2(2))**2-(pipl(3)+pim2(3))**2))
3098 prod1 =paa(4)*(pim1(4)-pipl(4))-paa(1)*(pim1(1)-pipl(1))
3099 $ -paa(2)*(pim1(2)-pipl(2))-paa(3)*(pim1(3)-pipl(3))
3100 prod2 =paa(4)*(pim2(4)-pipl(4))-paa(1)*(pim2(1)-pipl(1))
3101 $ -paa(2)*(pim2(2)-pipl(2))-paa(3)*(pim2(3)-pipl(3))
3103 vec1(i)= pim1(i)-pipl(i) -paa(i)*prod1/xmaa**2
3104 40 vec2(i)= pim2(i)-pipl(i) -paa(i)*prod2/xmaa**2
3106 IF (keya1.EQ.1)
THEN
3110 fnorm=fa1/sqrt(2.)*faropi*fro2pi
3112 hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gama1)
3113 $ *(cmplx(vec1(i))*amro**2*bwign(xmro1,amro,gamro)
3114 $ +cmplx(vec2(i))*amro**2*bwign(xmro2,amro,gamro))
3117 fnorm=2.0*sqrt(2.)/3.0/fpi
3118 gamax=gama1*gfun(xmaa**2)/gfun(ama1**2)
3120 hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gamax)
3121 $ *(cmplx(vec1(i))*fpik(xmro1)
3122 $ +cmplx(vec2(i))*fpik(xmro2))
3127 CALL clvec(hadcur,pn,pivec)
3128 CALL claxi(hadcur,pn,piaks)
3129 CALL clnut(hadcur,brakm,hvm)
3131 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3132 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3133 amplit=(gfermi*ccabib)**2*brak/2.
3138 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3139 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3149 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3150 * ,ampiz,ampi,amro,gamro,ama1,gama1
3151 * ,amk,amkz,amkst,gamkst
3153 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3154 * ,ampiz,ampi,amro,gamro,ama1,gama1
3155 * ,amk,amkz,amkst,gamkst
3157 IF (qkwa.LT.(amro+ampi)**2)
THEN
3158 gfun=4.1*(qkwa-9*ampiz**2)**3
3159 $ *(1.-3.3*(qkwa-9*ampiz**2)+5.8*(qkwa-9*ampiz**2)**2)
3161 gfun=qkwa*(1.623+10.38/qkwa-9.32/qkwa**2+0.65/qkwa**3)
3164 COMPLEX FUNCTION bwigs(S,M,G)
3169 REAL PI,PIM,QS,QM,W,GS,MK
3176 p(a,b,c)=sqrt(abs(abs(((a+b-c)**2-4.*a*b)/4./a)
3177 $ +(((a+b-c)**2-4.*a*b)/4./a))/2.0)
3193 qs=p(s,pim**2,mk**2)
3194 qm=p(m**2,pim**2,mk**2)
3196 gs=g*(m/w)*(qs/qm)**3
3198 bw=m**2/cmplx(m**2-s,-m*gs)
3199 qpm=p(pm**2,pim**2,mk**2)
3200 g1=pg*(pm/w)*(qs/qpm)**3
3201 bwp=pm**2/cmplx(pm**2-s,-pm*g1)
3202 bwigs= (bw+pbeta*bwp)/(1+pbeta)
3206 COMPLEX FUNCTION bwig(S,M,G)
3211 REAL PI,PIM,QS,QM,W,GS
3220 IF (s.GT.4.*pim**2)
THEN
3221 qs=sqrt(abs(abs(s/4.-pim**2)+(s/4.-pim**2))/2.0)
3222 qm=sqrt(m**2/4.-pim**2)
3224 gs=g*(m/w)*(qs/qm)**3
3228 bwig=m**2/cmplx(m**2-s,-m*gs)
3231 COMPLEX FUNCTION fpik(W)
3236 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
3241 IF (init.EQ.0 )
THEN
3255 fpik= (bwig(s,rom,rog)+beta1*bwig(s,rom1,rog1))
3264 fpirho=cabs(fpik(w))**2
3266 SUBROUTINE clvec(HJ,PN,PIV)
3276 hn= hj(4)*cmplx(pn(4))-hj(3)*cmplx(pn(3))
3277 hh= real(hj(4)*conjg(hj(4))-hj(3)*conjg(hj(3))
3278 $ -hj(2)*conjg(hj(2))-hj(1)*conjg(hj(1)))
3280 10 piv(i)=4.*real(hn*conjg(hj(i)))-2.*hh*pn(i)
3283 SUBROUTINE claxi(HJ,PN,PIA)
3290 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3291 COMMON / idfc / idff
3293 COMPLEX HJ(4),HJC(4)
3296 det2(i,j)=aimag(hjc(i)*hj(j)-hjc(j)*hj(i))
3301 IF (ktom.EQ.1.OR.ktom.EQ.-1)
THEN
3302 sign= idff/abs(idff)
3303 ELSEIF (ktom.EQ.2)
THEN
3304 sign=-idff/abs(idff)
3306 print *,
'STOP IN CLAXI: KTOM=',ktom
3311 10 hjc(i)=conjg(hj(i))
3312 pia(1)= -2.*pn(3)*det2(2,4)+2.*pn(4)*det2(2,3)
3313 pia(2)= -2.*pn(4)*det2(1,3)+2.*pn(3)*det2(1,4)
3314 pia(3)= 2.*pn(4)*det2(1,2)
3315 pia(4)= 2.*pn(3)*det2(1,2)
3318 20 pia(i)=pia(i)*sign
3320 SUBROUTINE clnut(HJ,B,HV)
3332 b=real( hj(4)*aimag(hj(4)) - hj(3)*aimag(hj(3))
3333 & - hj(2)*aimag(hj(2)) - hj(1)*aimag(hj(1)) )
3336 SUBROUTINE dampog(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3348 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3349 * ,ampiz,ampi,amro,gamro,ama1,gama1
3350 * ,amk,amkz,amkst,gamkst
3352 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3353 * ,ampiz,ampi,amro,gamro,ama1,gama1
3354 * ,amk,amkz,amkst,gamkst
3355 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3356 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3357 COMMON /testa1/ keya1
3358 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4)
3359 REAL PAA(4),VEC1(4),VEC2(4)
3360 REAL PIVEC(4),PIAKS(4),HVM(4)
3361 COMPLEX BWIGN,HADCUR(4),FNORM,FORMOM
3373 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3376 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3377 xmom =sqrt(abs( (pim2(4)+pipl(4))**2-(pim2(3)+pipl(3))**2
3378 $ -(pim2(2)+pipl(2))**2-(pim2(1)+pipl(1))**2 ))
3379 xmro2 =(pipl(1))**2 +(pipl(2))**2 +(pipl(3))**2
3381 prod1 =vec1(1)*pipl(1)
3382 prod2 =vec2(2)*pipl(2)
3383 p12 =pim1(4)*pim2(4)-pim1(1)*pim2(1)
3384 $ -pim1(2)*pim2(2)-pim1(3)*pim2(3)
3385 p1pl =pim1(4)*pipl(4)-pim1(1)*pipl(1)
3386 $ -pim1(2)*pipl(2)-pim1(3)*pipl(3)
3387 p2pl =pipl(4)*pim2(4)-pipl(1)*pim2(1)
3388 $ -pipl(2)*pim2(2)-pipl(3)*pim2(3)
3390 vec1(i)= (vec1(i)-prod1/xmro2*pipl(i))
3392 gnorm=sqrt(vec1(1)**2+vec1(2)**2+vec1(3)**2)
3394 vec1(i)= vec1(i)/gnorm
3396 vec2(1)=(vec1(2)*pipl(3)-vec1(3)*pipl(2))/sqrt(xmro2)
3397 vec2(2)=(vec1(3)*pipl(1)-vec1(1)*pipl(3))/sqrt(xmro2)
3398 vec2(3)=(vec1(1)*pipl(2)-vec1(2)*pipl(1))/sqrt(xmro2)
3399 p1vec1 =pim1(4)*vec1(4)-pim1(1)*vec1(1)
3400 $ -pim1(2)*vec1(2)-pim1(3)*vec1(3)
3401 p2vec1 =vec1(4)*pim2(4)-vec1(1)*pim2(1)
3402 $ -vec1(2)*pim2(2)-vec1(3)*pim2(3)
3403 p1vec2 =pim1(4)*vec2(4)-pim1(1)*vec2(1)
3404 $ -pim1(2)*vec2(2)-pim1(3)*vec2(3)
3405 p2vec2 =vec2(4)*pim2(4)-vec2(1)*pim2(1)
3406 $ -vec2(2)*pim2(2)-vec2(3)*pim2(3)
3408 fnorm=formom(xmaa,xmom)
3413 hadcur(i) = fnorm *(
3414 $ vec1(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3415 $ -pim2(i)*(p2vec1*p1pl-p1vec1*p2pl)
3416 $ +pipl(i)*(p2vec1*p12 -p1vec1*(ampi**2+p2pl)) )
3418 hadcur(i) = fnorm *(
3419 $ vec2(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3420 $ -pim2(i)*(p2vec2*p1pl-p1vec2*p2pl)
3421 $ +pipl(i)*(p2vec2*p12 -p1vec2*(ampi**2+p2pl)) )
3426 CALL clvec(hadcur,pn,pivec)
3427 CALL claxi(hadcur,pn,piaks)
3428 CALL clnut(hadcur,brakm,hvm)
3430 brak=brak+(gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3431 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3433 hv(i)=hv(i)-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3434 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3438 amplit=(gfermi*ccabib)**2*brak/2.
3447 SUBROUTINE damppk(MNUM,PT,PN,PIM1,PIM2,PIM3,AMPLIT,HV)
3459 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3460 * ,ampiz,ampi,amro,gamro,ama1,gama1
3461 * ,amk,amkz,amkst,gamkst
3463 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3464 * ,ampiz,ampi,amro,gamro,ama1,gama1
3465 * ,amk,amkz,amkst,gamkst
3466 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3467 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3468 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4)
3469 REAL PAA(4),VEC1(4),VEC2(4),VEC3(4),VEC4(4),VEC5(4)
3470 REAL PIVEC(4),PIAKS(4),HVM(4)
3471 REAL FNORM(0:19),COEF(1:5,0:19)
3472 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5,UROJ
3474 COMPLEX F1,F2,F3,F4,F5
3476 EXTERNAL form1,form2,form3,form4,form5
3477 DATA pi /3.141592653589793238462643/
3481 IF (icont.EQ.0)
THEN
3489 fnorm(4)=scabib/fpi/dwapi0
3499 coef(1,0)= 2.0*sqrt(2.)/3.0
3500 coef(2,0)=-2.0*sqrt(2.)/3.0
3503 coef(3,0)= 2.0*sqrt(2.)/3.0
3508 coef(1,1)=-sqrt(2.)/3.0
3509 coef(2,1)= sqrt(2.)/3.0
3514 coef(1,2)=-sqrt(2.)/3.0
3515 coef(2,2)= sqrt(2.)/3.0
3529 coef(1,4)= 1.0/sqrt(2.)/3.0
3530 coef(2,4)=-1.0/sqrt(2.)/3.0
3535 coef(1,5)=-sqrt(2.)/3.0
3536 coef(2,5)= sqrt(2.)/3.0
3554 coef(5,7)=-sqrt(2.0/3.0)
3563 coef(1,9)= 2.0*sqrt(2.)/3.0
3564 coef(2,9)=-2.0*sqrt(2.)/3.0
3565 coef(3,9)= 2.0*sqrt(2.)/3.0
3571 coef(1,k)= 2.0*sqrt(2.)/3.0
3572 coef(2,k)=-2.0*sqrt(2.)/3.0
3573 coef(3,k)= 2.0*sqrt(2.)/3.0
3584 10 paa(i)=pim1(i)+pim2(i)+pim3(i)
3585 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3586 xmro1 =sqrt(abs((pim3(4)+pim2(4))**2-(pim3(1)+pim2(1))**2
3587 $ -(pim3(2)+pim2(2))**2-(pim3(3)+pim2(3))**2))
3588 xmro2 =sqrt(abs((pim3(4)+pim1(4))**2-(pim3(1)+pim1(1))**2
3589 $ -(pim3(2)+pim1(2))**2-(pim3(3)+pim1(3))**2))
3590 xmro3 =sqrt(abs((pim1(4)+pim2(4))**2-(pim1(1)+pim2(1))**2
3591 $ -(pim1(2)+pim2(2))**2-(pim1(3)+pim2(3))**2))
3593 prod1 =paa(4)*(pim2(4)-pim3(4))-paa(1)*(pim2(1)-pim3(1))
3594 $ -paa(2)*(pim2(2)-pim3(2))-paa(3)*(pim2(3)-pim3(3))
3595 prod2 =paa(4)*(pim3(4)-pim1(4))-paa(1)*(pim3(1)-pim1(1))
3596 $ -paa(2)*(pim3(2)-pim1(2))-paa(3)*(pim3(3)-pim1(3))
3597 prod3 =paa(4)*(pim1(4)-pim2(4))-paa(1)*(pim1(1)-pim2(1))
3598 $ -paa(2)*(pim1(2)-pim2(2))-paa(3)*(pim1(3)-pim2(3))
3600 vec1(i)= pim2(i)-pim3(i) -paa(i)*prod1/xmaa**2
3601 vec2(i)= pim3(i)-pim1(i) -paa(i)*prod2/xmaa**2
3602 vec3(i)= pim1(i)-pim2(i) -paa(i)*prod3/xmaa**2
3603 40 vec4(i)= pim1(i)+pim2(i)+pim3(i)
3604 CALL prod5(pim1,pim2,pim3,vec5)
3609 f1 = cmplx(coef(1,mnum))*form1(mnum,xmaa**2,xmro1**2,xmro2**2)
3610 f2 = cmplx(coef(2,mnum))*form2(mnum,xmaa**2,xmro2**2,xmro1**2)
3611 f3 = cmplx(coef(3,mnum))*form3(mnum,xmaa**2,xmro3**2,xmro1**2)
3613 $cmplx(coef(4,mnum))*form4(mnum,xmaa**2,xmro1**2,xmro2**2,xmro3**2)
3614 f5 = (-1.0)*uroj/4.0/pi**2/fpi**2*
3615 $ cmplx(coef(5,mnum))*form5(mnum,xmaa**2,xmro1**2,xmro2**2)
3618 hadcur(i)= cmplx(fnorm(mnum)) * (
3619 $ cmplx(vec1(i))*f1+cmplx(vec2(i))*f2+cmplx(vec3(i))*f3+
3620 $ cmplx(vec4(i))*f4+cmplx(vec5(i))*f5)
3625 CALL clvec(hadcur,pn,pivec)
3626 CALL claxi(hadcur,pn,piaks)
3627 CALL clnut(hadcur,brakm,hvm)
3629 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3630 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3631 amplit=(gfermi)**2*brak/2.
3632 IF (mnum.GE.20)
THEN
3633 print *,
'MNUM=',mnum
3639 IF (k.EQ.4) znak=1.0
3640 xm1=znak*pim1(k)**2+xm1
3641 xm2=znak*pim2(k)**2+xm2
3642 xm3=znak*pim3(k)**2+xm3
3643 77 print *,
'PIM1=',pim1(k),
'PIM2=',pim2(k),
'PIM3=',pim3(k)
3644 print *,
'XM1=',sqrt(xm1),
'XM2=',sqrt(xm2),
'XM3=',sqrt(xm3)
3645 print *,
'************************************************'
3649 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3650 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3655 SUBROUTINE prod5(P1,P2,P3,PIA)
3661 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3662 COMMON / idfc / idff
3663 REAL PIA(4),P1(4),P2(4),P3(4)
3664 det2(i,j)=p1(i)*p2(j)-p2(i)*p1(j)
3666 IF (ktom.EQ.1.OR.ktom.EQ.-1)
THEN
3667 sign= idff/abs(idff)
3668 ELSEIF (ktom.EQ.2)
THEN
3669 sign=-idff/abs(idff)
3671 print *,
'STOP IN PROD5: KTOM=',ktom
3677 pia(1)= -p3(3)*det2(2,4)+p3(4)*det2(2,3)+p3(2)*det2(3,4)
3678 pia(2)= -p3(4)*det2(1,3)+p3(3)*det2(1,4)-p3(1)*det2(3,4)
3679 pia(3)= p3(4)*det2(1,2)-p3(2)*det2(1,4)+p3(1)*det2(2,4)
3680 pia(4)= p3(3)*det2(1,2)-p3(2)*det2(1,3)+p3(1)*det2(2,3)
3683 20 pia(i)=pia(i)*sign
3686 SUBROUTINE dexnew(MODE,ISGN,POL,PNU,PAA,PNPI,JNPI)
3699 COMMON / inout / inut,iout
3700 REAL POL(4),HV(4),PAA(4),PNU(4),PNPI(4,9),RN(1)
3706 CALL dadnew( -1,isgn,hv,pnu,paa,pnpi,jdumm)
3711 ELSEIF(mode.EQ. 0)
THEN
3714 IF(iwarm.EQ.0)
GOTO 902
3715 CALL dadnew( 0,isgn,hv,pnu,paa,pnpi,jnpi)
3716 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
3719 IF(rn(1).GT.wt)
GOTO 300
3721 ELSEIF(mode.EQ. 1)
THEN
3723 CALL dadnew( 1,isgn,hv,pnu,paa,pnpi,jdumm)
3728 902
WRITE(iout, 9020)
3729 9020
FORMAT(
' ----- DEXNEW: LACK OF INITIALISATION')
3732 SUBROUTINE dadnew(MODE,ISGN,HV,PNU,PWB,PNPI,JNPI)
3734 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3735 * ,ampiz,ampi,amro,gamro,ama1,gama1
3736 * ,amk,amkz,amkst,gamkst
3738 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3739 * ,ampiz,ampi,amro,gamro,ama1,gama1
3740 * ,amk,amkz,amkst,gamkst
3741 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3742 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3743 COMMON / taubmc / gampmc(500),gamper(500),nevdec(500)
3744 real*4 gampmc ,gamper
3746 COMMON / inout / inut,iout
3748 parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
3750 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3753 CHARACTER NAMES(NMODE)*31
3756 real*4 pnu(4),pwb(4),pnpi(4,9),hv(4),hhv(4)
3757 real*4 pdum1(4),pdum2(4),pdumi(4,9)
3760 real*8 swt(nmode),sswt(nmode)
3761 dimension nevraw(nmode),nevovr(nmode),nevacc(nmode)
3763 DATA pi /3.141592653589793238462643/
3784 a = pkorb(3,37+jnpi)
3797 ELSEIF(jnpi.LE.nm4)
THEN
3798 CALL dph4pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3800 ELSEIF(jnpi.LE.nm4+nm5)
THEN
3801 CALL dph5pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3802 ELSEIF(jnpi.LE.nm4+nm5+nm6)
THEN
3803 CALL dphnpi(wt,hv,pdum1,pdum2,pdumi,jnpi)
3804 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3)
THEN
3805 inum=jnpi-nm4-nm5-nm6
3806 CALL dphspk(wt,hv,pdum1,pdum2,pdumi,inum)
3807 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2)
THEN
3808 inum=jnpi-nm4-nm5-nm6-nm3
3809 CALL dphsrk(wt,hv,pdum1,pdum2,pdumi,inum)
3813 IF(wt.GT.wtmax(jnpi)/1.2) wtmax(jnpi)=wt*1.2
3823 ELSEIF(mode.EQ. 0)
THEN
3825 IF(iwarm.EQ.0)
GOTO 902
3830 ELSEIF(jnpi.LE.nm4)
THEN
3831 CALL dph4pi(wt,hhv,pnu,pwb,pnpi,jnpi)
3832 ELSEIF(jnpi.LE.nm4+nm5)
THEN
3833 CALL dph5pi(wt,hhv,pnu,pwb,pnpi,jnpi)
3834 ELSEIF(jnpi.LE.nm4+nm5+nm6)
THEN
3835 CALL dphnpi(wt,hhv,pnu,pwb,pnpi,jnpi)
3836 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3)
THEN
3837 inum=jnpi-nm4-nm5-nm6
3838 CALL dphspk(wt,hhv,pnu,pwb,pnpi,inum)
3839 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2)
THEN
3840 inum=jnpi-nm4-nm5-nm6-nm3
3841 CALL dphsrk(wt,hhv,pnu,pwb,pnpi,inum)
3849 nevraw(jnpi)=nevraw(jnpi)+1
3850 swt(jnpi)=swt(jnpi)+wt
3854 sswt(jnpi)=sswt(jnpi)+dble(wt)**2
3859 IF(wt.GT.wtmax(jnpi)) nevovr(jnpi)=nevovr(jnpi)+1
3860 IF(rn*wtmax(jnpi).GT.wt)
GOTO 300
3862 costhe=-1.+2.*rrr(2)
3865 CALL rotor2(thet,pnu,pnu)
3866 CALL rotor3( phi,pnu,pnu)
3867 CALL rotor2(thet,pwb,pwb)
3868 CALL rotor3( phi,pwb,pwb)
3869 CALL rotor2(thet,hv,hv)
3870 CALL rotor3( phi,hv,hv)
3873 CALL rotor2(thet,pnpi(1,i),pnpi(1,i))
3874 CALL rotor3( phi,pnpi(1,i),pnpi(1,i))
3876 nevacc(jnpi)=nevacc(jnpi)+1
3878 ELSEIF(mode.EQ. 1)
THEN
3881 IF(nevraw(jnpi).EQ.0)
GOTO 500
3882 pargam=swt(jnpi)/float(nevraw(jnpi)+1)
3884 IF(nevraw(jnpi).NE.0)
3885 & error=sqrt(sswt(jnpi)/swt(jnpi)**2-1./float(nevraw(jnpi)))
3887 WRITE(iout, 7010) names(jnpi),
3888 & nevraw(jnpi),nevacc(jnpi),nevovr(jnpi),pargam,rat,error
3890 gampmc(8+jnpi-1)=rat
3891 gamper(8+jnpi-1)=error
3897 7003
FORMAT(///1x,15(5h*****)
3898 $ /,
' *', 25x,
'******** DADNEW INITIALISATION ********',9x,1h*
3900 7004
FORMAT(
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*/)
3902 $ /,1x,15(5h*****)/)
3903 7010
FORMAT(///1x,15(5h*****)
3904 $ /,
' *', 25x,
'******** DADNEW FINAL REPORT ******** ',9x,1h*
3905 $ /,
' *', 25x,
'CHANNEL:',a31 ,9x,1h*
3906 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF DECAYS TOTAL ',9x,1h*
3907 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF DECAYS ACCEPTED ',9x,1h*
3908 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
3909 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH IN GEV UNITS ',9x,1h*
3910 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
3911 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
3912 $ /,1x,15(5h*****)/)
3913 902
WRITE(iout, 9020)
3914 9020
FORMAT(
' ----- DADNEW: LACK OF INITIALISATION')
3916 903
WRITE(iout, 9030) jnpi,mode
3917 9030
FORMAT(
' ----- DADNEW: WRONG JNPI',2i5)
3922 SUBROUTINE dph4pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
3929 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3930 * ,ampiz,ampi,amro,gamro,ama1,gama1
3931 * ,amk,amkz,amkst,gamkst
3933 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3934 * ,ampiz,ampi,amro,gamro,ama1,gama1
3935 * ,amk,amkz,amkst,gamkst
3936 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3937 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3940 REAL HV(4),PT(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4),PMULT(4,9)
3943 real*8 uu,ff,ff1,ff2,ff3,ff4,gg1,gg2,gg3,gg4,rr
3944 DATA pi /3.141592653589793238462643/
3946 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
3951 phspac=1./2**23/pi**11
3986 CALL choice(100+jnpi,rrb,ichan,prob1,prob2,prob3,
3988 $ amrop,gamrop,amrx,gamrx,amrb,gamrb)
4002 ams1=(amp1+amp2+amp3+amp4)**2
4003 ams2=(amtau-amnuta)**2
4004 alp1=atan((ams1-amrop**2)/amrop/gamrop)
4005 alp2=atan((ams2-amrop**2)/amrop/gamrop)
4006 alp=alp1+rr1*(alp2-alp1)
4007 am4sq =amrop**2+amrop*gamrop*tan(alp)
4010 $ ((am4sq-amrop**2)**2+(amrop*gamrop)**2)/(amrop*gamrop)
4011 phspac=phspac*(alp2-alp1)
4015 ams1=(amp2+amp3+amp4)**2
4017 IF (rrr(9).GT.prez)
THEN
4018 am3sq=ams1+ rr1*(ams2-ams1)
4024 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4025 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4026 alp=alp1+rr1*(alp2-alp1)
4027 am3sq =amrx**2+amrx*gamrx*tan(alp)
4030 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4038 am2sq=ams1+ rr2*(ams2-ams1)
4044 enq1=(am2sq-amp3**2+amp4**2)/(2*am2)
4045 enq2=(am2sq+amp3**2-amp4**2)/(2*am2)
4046 ppi= enq1**2-amp4**2
4047 pppi=sqrt(abs(enq1**2-amp4**2))
4049 phspac=phspac*(4*pi)*(2*pppi/am2)
4053 CALL sphera(pppi,piz)
4067 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4068 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4069 ppi = pr(4)**2-am2**2
4075 pim1(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4078 ff3=(4*pi)*(2*pr(3)/am3)
4080 exe=(pr(4)+pr(3))/am2
4081 CALL bostr3(exe,piz,piz)
4082 CALL bostr3(exe,pipl,pipl)
4085 thet =acos(-1.+2*rr3)
4087 CALL rotpol(thet,phi,pipl)
4088 CALL rotpol(thet,phi,pim1)
4089 CALL rotpol(thet,phi,piz)
4090 CALL rotpol(thet,phi,pr)
4097 pr(4)=1./(2*am4)*(am4**2+am3**2-amp1**2)
4098 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4099 ppi = pr(4)**2-am3**2
4105 pim2(4)=1./(2*am4)*(am4**2-am3**2+amp1**2)
4108 ff4=(4*pi)*(2*pr(3)/am4)
4110 exe=(pr(4)+pr(3))/am3
4111 CALL bostr3(exe,piz,piz)
4112 CALL bostr3(exe,pipl,pipl)
4113 CALL bostr3(exe,pim1,pim1)
4116 thet =acos(-1.+2*rr3)
4118 CALL rotpol(thet,phi,pipl)
4119 CALL rotpol(thet,phi,pim1)
4120 CALL rotpol(thet,phi,pim2)
4121 CALL rotpol(thet,phi,piz)
4122 CALL rotpol(thet,phi,pr)
4128 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am4**2)
4129 paa(3)= sqrt(abs(paa(4)**2-am4**2))
4130 ppi = paa(4)**2-am4**2
4131 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4132 phsp=phsp*(4*pi)*(2*paa(3)/amtau)
4136 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am4**2)
4139 IF(rrr(9).LE.0.5*prez)
THEN
4148 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4149 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4151 ams1=(amp1+amp3+amp4)**2
4154 ams2=(sqrt(am3sq)-amp1)**2
4156 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4157 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4160 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4161 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4163 ams1=(amp1+amp3+amp4)**2
4164 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4165 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4166 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4169 ams2=(sqrt(am3sq)-amp1)**2
4171 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4172 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4175 am3sq=(pim2(4)+piz(4)+pipl(4))**2-(pim2(3)+piz(3)+pipl(3))**2
4176 $ -(pim2(2)+piz(2)+pipl(2))**2-(pim2(1)+piz(1)+pipl(1))**2
4178 ams1=(amp2+amp3+amp4)**2
4179 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4180 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4181 gg1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4184 ams2=(sqrt(am3sq)-amp2)**2
4186 gg3=(4*pi)*(xlam(am2**2,amp2**2,am3sq)/am3sq)
4187 gg4=(4*pi)*(xlam(am3sq,amp1**2,am4**2)/am4**2)
4190 IF ( ( (ff+gg)*uu+ff*gg ).GT.0.0d0)
THEN
4191 rr=ff*gg*uu/(0.5*prez*(ff+gg)*uu+(1.0-prez)*ff*gg)
4219 exe=(paa(4)+paa(3))/am4
4220 CALL bostr3(exe,piz,piz)
4221 CALL bostr3(exe,pipl,pipl)
4222 CALL bostr3(exe,pim1,pim1)
4223 CALL bostr3(exe,pim2,pim2)
4224 CALL bostr3(exe,pr,pr)
4234 CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4235 ELSEIF (jnpi.EQ.2)
THEN
4236 CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv)
4238 CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv)
4241 dgamt=1/(2.*amtau)*amplit*phspac
4253 SUBROUTINE dam4pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,AMPLIT,HV)
4265 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4266 * ,ampiz,ampi,amro,gamro,ama1,gama1
4267 * ,amk,amkz,amkst,gamkst
4269 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4270 * ,ampiz,ampi,amro,gamro,ama1,gama1
4271 * ,amk,amkz,amkst,gamkst
4272 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4273 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4274 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4),PIM4(4)
4275 REAL PIVEC(4),PIAKS(4),HVM(4)
4276 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5
4277 EXTERNAL form1,form2,form3,form4,form5
4278 DATA pi /3.141592653589793238462643/
4282 CALL curr_cleo(mnum,pim1,pim2,pim3,pim4,hadcur)
4286 CALL clvec(hadcur,pn,pivec)
4287 CALL claxi(hadcur,pn,piaks)
4288 CALL clnut(hadcur,brakm,hvm)
4290 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4291 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4292 amplit=(ccabib*gfermi)**2*brak/2.
4295 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4296 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4302 SUBROUTINE dam5pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,PIM5,AMPLIT,HV)
4314 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4315 * ,ampiz,ampi,amro,gamro,ama1,gama1
4316 * ,amk,amkz,amkst,gamkst
4318 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4319 * ,ampiz,ampi,amro,gamro,ama1,gama1
4320 * ,amk,amkz,amkst,gamkst
4321 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4322 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4323 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4),PIM4(4),PIM5(4)
4324 REAL PIVEC(4),PIAKS(4),HVM(4)
4325 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5
4326 EXTERNAL form1,form2,form3,form4,form5
4327 DATA pi /3.141592653589793238462643/
4331 CALL curr5(mnum,pim1,pim2,pim3,pim4,pim5,hadcur)
4335 CALL clvec(hadcur,pn,pivec)
4336 CALL claxi(hadcur,pn,piaks)
4337 CALL clnut(hadcur,brakm,hvm)
4339 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4340 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4341 amplit=(ccabib*gfermi)**2*brak/2.
4344 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4345 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4351 SUBROUTINE dph5pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4356 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4357 * ,ampiz,ampi,amro,gamro,ama1,gama1
4358 * ,amk,amkz,amkst,gamkst
4360 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4361 * ,ampiz,ampi,amro,gamro,ama1,gama1
4364 * ,amk,amkz,amkst,gamkst
4365 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4366 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4367 parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
4369 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4372 CHARACTER NAMES(NMODE)*31
4373 REAL HV(4),PT(4),PN(4),PAA(4),PMULT(4,9)
4374 real*4 pr(4),pi1(4),pi2(4),pi3(4),pi4(4),pi5(4)
4375 real*8 amp1,amp2,amp3,amp4,amp5,ams1,ams2,amom,gamom
4376 real*8 am5sq,am4sq,am3sq,am2sq,am5,am4,am3
4378 real*8 gg1,gg2,gg3,ff1,ff2,ff3,ff4,alp,alp1,alp2
4385 DATA pi /3.141592653589793238462643/
4392 bwign(xm,am,gamma)=xm**2/cmplx(xm**2-am**2,gamma*am)
4403 IF (jnpi.EQ.23.OR.jnpi.EQ.24) gamom= 0.0085
4408 phspac=1./2**29/pi**14
4411 amp1=dcdmas(idffin(1,jnpi))
4412 amp2=dcdmas(idffin(2,jnpi))
4413 amp3=dcdmas(idffin(3,jnpi))
4414 amp4=dcdmas(idffin(4,jnpi))
4415 amp5=dcdmas(idffin(5,jnpi))
4429 IF (rrr(11).GT.proba2)
THEN
4432 ams1=(amp1+amp2+amp3+amp4+amp5)**2
4433 ams2=(amtau-amnuta)**2
4434 alp1=atan((ams1-ama2**2)/ama2/gama2)
4435 alp2=atan((ams2-ama2**2)/ama2/gama2)
4436 am5sq=ams1+ rr1*(ams2-ams1)
4442 ams1=(amp1+amp2+amp3+amp4+amp5)**2
4443 ams2=(amtau-amnuta)**2
4444 alp1=atan((ams1-ama2**2)/ama2/gama2)
4445 alp2=atan((ams2-ama2**2)/ama2/gama2)
4446 alp=alp1+rr1*(alp2-alp1)
4447 am5sq =ama2**2+ama2*gama2*tan(alp)
4451 gg5=((am5sq-ama2**2)**2+(ama2*gama2)**2)/(ama2*gama2)
4453 phspac=phspac/(proba2/gg5+(1d0-proba2)/(ams2-ams1) )
4458 ams1=(amp2+amp3+amp4+amp5)**2
4460 am4sq=ams1+ rr1*(ams2-ams1)
4466 IF (rrr(12).LT.probom)
THEN
4469 ams1=(amp2+amp3+amp4)**2
4471 alp1=atan((ams1-amom**2)/amom/gamom)
4472 alp2=atan((ams2-amom**2)/amom/gamom)
4473 alp=alp1+rr1*(alp2-alp1)
4474 am3sq =amom**2+amom*gamom*tan(alp)
4479 ams1=(amp2+amp3+amp4)**2
4481 alp1=atan((ams1-amom**2)/amom/gamom)
4482 alp2=atan((ams2-amom**2)/amom/gamom)
4484 am3sq=ams1+ rr1*(ams2-ams1)
4489 gg2=((am3sq-amom**2)**2+(amom*gamom)**2)/(amom*gamom)
4491 gg2=1d0/(probom/gg2+(1d0-probom)/(ams2-ams1))
4498 am2sq=ams1+ rr2*(ams2-ams1)
4504 enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4505 enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4506 ppi= enq1**2-amp3**2
4507 pppi=sqrt(abs(enq1**2-amp3**2))
4508 ff1=(4*pi)*(2*pppi/am2)
4510 call sphera(pppi,pi3)
4521 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4522 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4523 ppi = pr(4)**2-am2**2
4527 pi2(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4530 ff2=(4*pi)*(2*pr(3)/am3)
4532 exe=(pr(4)+pr(3))/am2
4533 call bostr3(exe,pi3,pi3)
4534 call bostr3(exe,pi4,pi4)
4537 thet =acos(-1.+2*rr3)
4539 call rotpol(thet,phi,pi2)
4540 call rotpol(thet,phi,pi3)
4541 call rotpol(thet,phi,pi4)
4547 pr(4)=1./(2*am4)*(am4**2+am3**2-amp5**2)
4548 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4549 ppi = pr(4)**2-am3**2
4553 pi5(4)=1./(2*am4)*(am4**2-am3**2+amp5**2)
4556 ff3=(4*pi)*(2*pr(3)/am4)
4558 exe=(pr(4)+pr(3))/am3
4559 call bostr3(exe,pi2,pi2)
4560 call bostr3(exe,pi3,pi3)
4561 call bostr3(exe,pi4,pi4)
4564 thet =acos(-1.+2*rr3)
4566 call rotpol(thet,phi,pi2)
4567 call rotpol(thet,phi,pi3)
4568 call rotpol(thet,phi,pi4)
4569 call rotpol(thet,phi,pi5)
4575 pr(4)=1./(2*am5)*(am5**2+am4**2-amp1**2)
4576 pr(3)= sqrt(abs(pr(4)**2-am4**2))
4577 ppi = pr(4)**2-am4**2
4581 pi1(4)=1./(2*am5)*(am5**2-am4**2+amp1**2)
4584 ff4=(4*pi)*(2*pr(3)/am5)
4586 exe=(pr(4)+pr(3))/am4
4587 call bostr3(exe,pi2,pi2)
4588 call bostr3(exe,pi3,pi3)
4589 call bostr3(exe,pi4,pi4)
4590 call bostr3(exe,pi5,pi5)
4593 thet =acos(-1.+2*rr3)
4595 call rotpol(thet,phi,pi1)
4596 call rotpol(thet,phi,pi2)
4597 call rotpol(thet,phi,pi3)
4598 call rotpol(thet,phi,pi4)
4599 call rotpol(thet,phi,pi5)
4608 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5sq)
4609 paa(3)= sqrt(abs(paa(4)**2-am5sq))
4610 ppi = paa(4)**2-am5sq
4611 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4615 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am5**2)
4618 phspac=phspac * gg1*gg2*gg3*ff1*ff2*ff3*ff4
4622 exe=(paa(4)+paa(3))/am5
4623 call bostr3(exe,pi1,pi1)
4624 call bostr3(exe,pi2,pi2)
4625 call bostr3(exe,pi3,pi3)
4626 call bostr3(exe,pi4,pi4)
4627 call bostr3(exe,pi5,pi5)
4635 qxn=paa(4)*pn(4)-paa(1)*pn(1)-paa(2)*pn(2)-paa(3)*pn(3)
4636 brak=2*(gv**2+ga**2)*(2*pxq*qxn+am5sq*pxn)
4637 & -6*(gv**2-ga**2)*amtau*amnuta*am5sq
4638 fompp = cabs(bwign(am3,amom,gamom))**2
4643 amplit=ccabib**2*gfermi**2/2. * brak
4644 amplit = amplit * fompp * fnorm
4657 $
CALL dam5pi(jnpi,pt,pn,pi1,pi2,pi3,pi4,pi5,amplit,hv)
4658 dgamt=1/(2.*amtau)*amplit*phspac
4675 SUBROUTINE dphsrk(DGAMT,HV,PN,PR,PMULT,INUM)
4681 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4682 * ,ampiz,ampi,amro,gamro,ama1,gama1
4683 * ,amk,amkz,amkst,gamkst
4685 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4686 * ,ampiz,ampi,amro,gamro,ama1,gama1
4687 * ,amk,amkz,amkst,gamkst
4688 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4689 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4690 REAL HV(4),PT(4),PN(4),PR(4),PKC(4),PKZ(4),QQ(4),PMULT(4,9)
4694 DATA pi /3.141592653589793238462643/
4698 phspac=1./2**11/pi**5
4706 ams2=(amtau-amnuta)**2
4709 amx2=ams1+ rr1(1)*(ams2-ams1)
4711 phspac=phspac*(ams2-ams1)
4729 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
4730 pn(3)=-sqrt((pn(4)-amnuta)*(pn(4)+amnuta))
4734 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
4736 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
4739 enq1=(amx2+amk**2-amkz**2)/(2.*amx)
4740 enq2=(amx2-amk**2+amkz**2)/(2.*amx)
4741 pppi=sqrt((enq1-amk)*(enq1+amk))
4742 phspac=phspac*(4*pi)*(2*pppi/amx)
4744 CALL sphera(pppi,pkc)
4750 exe=(pr(4)+pr(3))/amx
4752 CALL bostr3(exe,pkc,pkc)
4753 CALL bostr3(exe,pkz,pkz)
4755 30 qq(i)=pkc(i)-pkz(i)
4757 pksd =pr(4)*pr(4)-pr(3)*pr(3)-pr(2)*pr(2)-pr(1)*pr(1)
4758 qqpks=pr(4)* qq(4)-pr(3)* qq(3)-pr(2)* qq(2)-pr(1)* qq(1)
4760 31 qq(i)=qq(i)-pr(i)*qqpks/pksd
4763 prodnq=pn(4)*qq(4)-pn(1)*qq(1)-pn(2)*qq(2)-pn(3)*qq(3)
4765 qq2= qq(4)**2-qq(1)**2-qq(2)**2-qq(3)**2
4766 brak=(gv**2+ga**2)*(2*prodpq*prodnq-prodpn*qq2)
4767 & +(gv**2-ga**2)*amtau*amnuta*qq2
4768 amplit=(gfermi*ccabib)**2*brak*2*fpirk(amx)
4769 dgamt=1/(2.*amtau)*amplit*phspac
4771 40 hv(i)=2*gv*ga*amtau*(2*prodnq*qq(i)-qq2*pn(i))/brak
4782 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4783 * ,ampiz,ampi,amro,gamro,ama1,gama1
4784 * ,amk,amkz,amkst,gamkst
4786 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4787 * ,ampiz,ampi,amro,gamro,ama1,gama1
4788 * ,amk,amkz,amkst,gamkst
4791 fpirk=cabs(fpikm(w,amk,amkz))**2
4794 COMPLEX FUNCTION fpikmk(W,XM1,XM2)
4799 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
4804 IF (init.EQ.0 )
THEN
4817 fpikmk=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
4828 SUBROUTINE dwrph(KTO,PHX)
4832 IMPLICIT REAL*8 (a-h,o-z)
4843 IF (qhot(4).GT.1.e-5)
CALL dwluph(kto,qhot)
4846 SUBROUTINE dwluph(KTO,PHOT)
4860 COMMON /taupos/ np1,np2
4864 IF (phot(4).LE.0.0)
RETURN
4867 IF((kto.EQ. 1).OR.(kto.EQ.11))
THEN
4874 IF(ktos.GT.10) ktos=ktos-10
4876 CALL tralo4(ktos,phot,phot,am)
4877 CALL filhep(0,1,22,nps,nps,0,0,phot,0.0,.true.)
4882 SUBROUTINE dwluel(KTO,ISGN,PNU,PWB,PEL,PNE)
4894 REAL PNU(4),PWB(4),PEL(4),PNE(4)
4896 COMMON /taupos/ np1,np2
4907 CALL tralo4(kto,pnu,pnu,am)
4908 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4911 CALL tralo4(kto,pwb,pwb,am)
4915 CALL tralo4(kto,pel,pel,am)
4916 CALL filhep(0,1,11*isgn,nps,nps,0,0,pel,am,.false.)
4919 CALL tralo4(kto,pne,pne,am)
4920 CALL filhep(0,1,-12*isgn,nps,nps,0,0,pne,am,.true.)
4924 SUBROUTINE dwlumu(KTO,ISGN,PNU,PWB,PMU,PNM)
4936 REAL PNU(4),PWB(4),PMU(4),PNM(4)
4938 COMMON /taupos/ np1,np2
4949 CALL tralo4(kto,pnu,pnu,am)
4950 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4953 CALL tralo4(kto,pwb,pwb,am)
4957 CALL tralo4(kto,pmu,pmu,am)
4958 CALL filhep(0,1,13*isgn,nps,nps,0,0,pmu,am,.false.)
4961 CALL tralo4(kto,pnm,pnm,am)
4962 CALL filhep(0,1,-14*isgn,nps,nps,0,0,pnm,am,.true.)
4966 SUBROUTINE dwlupi(KTO,ISGN,PPI,PNU)
4977 COMMON /taupos/ np1,np2
4987 CALL tralo4(kto,pnu,pnu,am)
4988 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4991 CALL tralo4(kto,ppi,ppi,am)
4992 CALL filhep(0,1,-211*isgn,nps,nps,0,0,ppi,am,.true.)
4996 SUBROUTINE dwluro(KTO,ISGN,PNU,PRHO,PIC,PIZ)
5008 REAL PNU(4),PRHO(4),PIC(4),PIZ(4)
5010 COMMON /taupos/ np1,np2
5021 CALL tralo4(kto,pnu,pnu,am)
5022 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5025 CALL tralo4(kto,prho,prho,am)
5026 CALL filhep(0,2,-213*isgn,nps,nps,0,0,prho,am,.true.)
5029 CALL tralo4(kto,pic,pic,am)
5030 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pic,am,.true.)
5033 CALL tralo4(kto,piz,piz,am)
5034 CALL filhep(0,1,111,-2,-2,0,0,piz,am,.true.)
5038 SUBROUTINE dwluaa(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
5051 REAL PNU(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
5053 COMMON /taupos/ np1,np2
5064 CALL tralo4(kto,pnu,pnu,am)
5065 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5068 CALL tralo4(kto,paa,paa,am)
5069 CALL filhep(0,1,-20213*isgn,nps,nps,0,0,paa,am,.true.)
5077 CALL tralo4(kto,pim2,pim2,am)
5078 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pim2,am,.true.)
5081 CALL tralo4(kto,pim1,pim1,am)
5082 CALL filhep(0,1,-211*isgn,-2,-2,0,0,pim1,am,.true.)
5085 CALL tralo4(kto,pipl,pipl,am)
5086 CALL filhep(0,1, 211*isgn,-3,-3,0,0,pipl,am,.true.)
5088 ELSE IF (jaa.EQ.2)
THEN
5093 CALL tralo4(kto,pim2,pim2,am)
5094 CALL filhep(0,1,111,-1,-1,0,0,pim2,am,.true.)
5097 CALL tralo4(kto,pim1,pim1,am)
5098 CALL filhep(0,1,111,-2,-2,0,0,pim1,am,.true.)
5101 CALL tralo4(kto,pipl,pipl,am)
5102 CALL filhep(0,1,-211*isgn,-3,-3,0,0,pipl,am,.true.)
5108 SUBROUTINE dwlukk (KTO,ISGN,PKK,PNU)
5118 COMMON /taupos/ np1,np2
5130 CALL tralo4 (kto,pnu,pnu,am)
5131 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5134 CALL tralo4 (kto,pkk,pkk,am)
5135 CALL filhep(0,1,-321*isgn,nps,nps,0,0,pkk,am,.true.)
5139 SUBROUTINE dwluks(KTO,ISGN,PNU,PKS,PKK,PPI,JKST)
5140 COMMON / taukle / bra1,brk0,brk0b,brks
5141 real*4 bra1,brk0,brk0b,brks
5153 REAL PNU(4),PKS(4),PKK(4),PPI(4),XIO(1)
5154 COMMON /taupos/ np1,np2
5165 CALL tralo4(kto,pnu,pnu,am)
5166 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5169 CALL tralo4(kto,pks,pks,am)
5170 CALL filhep(0,1,-323*isgn,nps,nps,0,0,pks,am,.true.)
5178 CALL tralo4(kto,ppi,ppi,am)
5179 CALL filhep(0,1,-211*isgn,-1,-1,0,0,ppi,am,.true.)
5182 IF (isgn.EQ.-1) bran=brk0
5185 IF(xio(1).GT.bran)
THEN
5191 CALL tralo4(kto,pkk,pkk,am)
5192 CALL filhep(0,1,k0type,-2,-2,0,0,pkk,am,.true.)
5194 ELSE IF(jkst.EQ.20)
THEN
5199 CALL tralo4(kto,ppi,ppi,am)
5200 CALL filhep(0,1,111,-1,-1,0,0,ppi,am,.true.)
5203 CALL tralo4(kto,pkk,pkk,am)
5204 CALL filhep(0,1,-321*isgn,-2,-2,0,0,pkk,am,.true.)
5210 SUBROUTINE dwlnew(KTO,ISGN,PNU,PWB,PNPI,MODE)
5220 parameter(nmode=86,nm1=0,nm2=11,nm3=19,nm4=22,nm5=21,nm6=13)
5222 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5225 COMMON /taupos/ np1,np2
5226 CHARACTER NAMES(NMODE)*31
5227 REAL PNU(4),PWB(4),PNPI(4,9)
5239 CALL tralo4(kto,pnu,pnu,am)
5240 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5243 CALL tralo4(kto,pwb,pwb,am)
5244 CALL filhep(0,1,-24*isgn,nps,nps,0,0,pwb,am,.true.)
5252 kfpi=lunpik(idffin(i,jnpi),-isgn)
5259 CALL tralo4(kto,ppi,ppi,am)
5260 CALL filhep(0,1,kfpi,-i,-i,0,0,ppi,am,.true.)
5273 IMPLICIT REAL*8 (a-h,o-z)
5275 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5277 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5289 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5290 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5299 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5300 DATA pi /3.141592653589793238462643d0/
5302 IF(abs(y).LT.abs(x))
THEN
5304 IF(x.LE.0d0) the=pi-the
5306 the=acos(x/sqrt(x**2+y**2))
5317 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5318 DATA pi /3.141592653589793238462643d0/
5320 IF(abs(y).LT.abs(x))
THEN
5322 IF(x.LE.0d0) the=pi-the
5324 the=acos(x/sqrt(x**2+y**2))
5326 IF(y.LT.0d0) the=2d0*pi-the
5329 SUBROUTINE rotod1(PH1,PVEC,QVEC)
5334 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5335 dimension pvec(4),qvec(4),rvec(4)
5343 qvec(2)= cs*rvec(2)-sn*rvec(3)
5344 qvec(3)= sn*rvec(2)+cs*rvec(3)
5348 SUBROUTINE rotod2(PH1,PVEC,QVEC)
5353 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5354 dimension pvec(4),qvec(4),rvec(4)
5361 qvec(1)= cs*rvec(1)+sn*rvec(3)
5363 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5367 SUBROUTINE rotod3(PH1,PVEC,QVEC)
5372 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5374 dimension pvec(4),qvec(4),rvec(4)
5380 qvec(1)= cs*rvec(1)-sn*rvec(2)
5381 qvec(2)= sn*rvec(1)+cs*rvec(2)
5385 SUBROUTINE bostr3(EXE,PVEC,QVEC)
5391 real*4 pvec(4),qvec(4),rvec(4)
5404 SUBROUTINE bostd3(EXE,PVEC,QVEC)
5410 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5411 dimension pvec(4),qvec(4),rvec(4)
5425 SUBROUTINE rotor1(PH1,PVEC,QVEC)
5430 real*4 pvec(4),qvec(4),rvec(4)
5438 qvec(2)= cs*rvec(2)-sn*rvec(3)
5439 qvec(3)= sn*rvec(2)+cs*rvec(3)
5442 SUBROUTINE rotor2(PH1,PVEC,QVEC)
5447 IMPLICIT REAL*4(a-h,o-z)
5448 real*4 pvec(4),qvec(4),rvec(4)
5455 qvec(1)= cs*rvec(1)+sn*rvec(3)
5457 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5460 SUBROUTINE rotor3(PHI,PVEC,QVEC)
5465 real*4 pvec(4),qvec(4),rvec(4)
5471 qvec(1)= cs*rvec(1)-sn*rvec(2)
5472 qvec(2)= sn*rvec(1)+cs*rvec(2)
5476 SUBROUTINE spherd(R,X)
5481 real*8 r,x(4),pi,costh,sinth
5483 DATA pi /3.141592653589793238462643d0/
5487 sinth=sqrt(1 -costh**2)
5488 x(1)=r*sinth*cos(2*pi*rrr(2))
5489 x(2)=r*sinth*sin(2*pi*rrr(2))
5493 SUBROUTINE rotpox(THET,PHI,PP)
5494 IMPLICIT REAL*8 (a-h,o-z)
5502 CALL rotod2(thet,pp,pp)
5503 CALL rotod3( phi,pp,pp)
5506 SUBROUTINE sphera(R,X)
5514 DATA pi /3.141592653589793238462643/
5518 sinth=sqrt(1.-costh**2)
5519 x(1)=r*sinth*cos(2*pi*rrr(2))
5520 x(2)=r*sinth*sin(2*pi*rrr(2))
5524 SUBROUTINE rotpol(THET,PHI,PP)
5531 CALL rotor2(thet,pp,pp)
5532 CALL rotor3( phi,pp,pp)
5535 SUBROUTINE ranmar(RVEC,LENV)
5551 COMMON / inout / inut,iout
5553 common/raset1/u(97),c,i97,j97
5554 parameter(modcns=1000000000)
5555 DATA ntot,ntot2,ijkl/-1,0,0/
5557 IF (ntot .GE. 0)
GO TO 50
5566 entry rmarin(ijklin, ntotin,ntot2n)
5575 ntot = max(ntotin,0)
5576 ntot2= max(ntot2n,0)
5581 kl = ijkl - 30082*ij
5582 i = mod(ij/177, 177) + 2
5583 j = mod(ij, 177) + 2
5584 k = mod(kl/169, 178) + 1
5586 WRITE(iout,201) ijkl,ntot,ntot2
5587 201
FORMAT(1x,
' RANMAR INITIALIZED: ',i10,2x,2i10)
5592 m = mod(mod(i*j,179)*k, 179)
5596 l = mod(53*l+1, 169)
5597 IF (mod(l*m,64) .GE. 32) s = s+t
5602 4 twom24 = 0.5*twom24
5604 cd = 7654321.*twom24
5605 cm = 16777213.*twom24
5610 DO 45 loop2= 1, ntot2+1
5612 IF (loop2 .EQ. ntot2+1) now=ntot
5613 IF (now .GT. 0)
THEN
5614 WRITE (iout,
'(A,I15)')
' RMARIN SKIPPING OVER ',now
5615 DO 40 idum = 1, ntot
5617 IF (uni .LT. 0.) uni=uni+1.
5620 IF (i97 .EQ. 0) i97=97
5622 IF (j97 .EQ. 0) j97=97
5624 IF (c .LT. 0.) c=c+cm
5628 IF (kalled .EQ. 1)
RETURN
5632 DO 100 ivec= 1, lenv
5634 IF (uni .LT. 0.) uni=uni+1.
5637 IF (i97 .EQ. 0) i97=97
5639 IF (j97 .EQ. 0) j97=97
5641 IF (c .LT. 0.) c=c+cm
5643 IF (uni .LT. 0.) uni=uni+1.
5645 IF (uni .EQ. 0.)
THEN
5648 IF (uni .EQ. 0.) uni= twom24*twom24
5653 IF (ntot .GE. modcns)
THEN
5655 ntot = ntot - modcns
5659 entry rmarut(ijklut,ntotut,ntot2t)
5667 IMPLICIT REAL*8(a-h,o-z)
5670 IF(x .LT.-1.0)
GO TO 1
5671 IF(x .LE. 0.5)
GO TO 2
5672 IF(x .EQ. 1.0)
GO TO 3
5673 IF(x .LE. 2.0)
GO TO 4
5677 z=z-0.5* log(abs(x))**2
5683 3 dilogt=1.64493406684822
5687 z=1.64493406684822 - log(x)* log(abs(t))
5688 5 y=2.66666666666666 *t+0.66666666666666
5689 b= 0.00000 00000 00001
5690 a=y*b +0.00000 00000 00004
5691 b=y*a-b+0.00000 00000 00011
5692 a=y*b-a+0.00000 00000 00037
5693 b=y*a-b+0.00000 00000 00121
5694 a=y*b-a+0.00000 00000 00398
5695 b=y*a-b+0.00000 00000 01312
5696 a=y*b-a+0.00000 00000 04342
5697 b=y*a-b+0.00000 00000 14437
5698 a=y*b-a+0.00000 00000 48274
5699 b=y*a-b+0.00000 00001 62421
5700 a=y*b-a+0.00000 00005 50291
5701 b=y*a-b+0.00000 00018 79117
5702 a=y*b-a+0.00000 00064 74338
5703 b=y*a-b+0.00000 00225 36705
5704 a=y*b-a+0.00000 00793 87055
5705 b=y*a-b+0.00000 02835 75385
5706 a=y*b-a+0.00000 10299 04264
5707 b=y*a-b+0.00000 38163 29463
5708 a=y*b-a+0.00001 44963 00557
5709 b=y*a-b+0.00005 68178 22718
5710 a=y*b-a+0.00023 20021 96094
5711 b=y*a-b+0.00100 16274 96164
5712 a=y*b-a+0.00468 63619 59447
5713 b=y*a-b+0.02487 93229 24228
5714 a=y*b-a+0.16607 30329 27855
5715 a=y*a-b+1.93506 43008 6996