54 COMMON / taubra / gamprt(30),jlist(30),nchan
61 IF(nchan.LE.0.OR.nchan.GT.30)
GOTO 902
68 IF(rrr(1).LT.cumul(i)/cumul(nchan)) ji=i
73 9020
FORMAT(
' ----- JAKER: WRONG NCHAN')
76 SUBROUTINE dekay(KTO,HX)
108 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
114 COMMON /taupos/ np1,np2
115 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
116 real*4 gampmc ,gamper
117 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
121 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
124 CHARACTER NAMES(NMODE)*31
125 COMMON / inout / inut,iout
129 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4),HDUM(4)
143 IF (iwarm.EQ.1) x=5/(iwarm-1)
145 WRITE(iout,7001) jak1,jak
146 WRITE(iout,7002) iver
150 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
151 CALL dadmel(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
152 CALL dadmmu(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
153 CALL dadmpi(-1,idum,hdum,pdum1,pdum2)
154 CALL dadmro(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4)
155 CALL dadmaa(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5,jdum)
156 CALL dadmkk(-1,idum,hdum,pdum1,pdum2)
157 CALL dadmks(-1,idum,hdum,pdum1,pdum2,pdum3,pdum4,jdum)
158 CALL dadnew(-1,idum,hdum,pdum1,pdum2,pdumx,jdum)
164 ELSEIF(kto.EQ.1)
THEN
168 IF(iwarm.EQ.0)
GOTO 902
176 CALL dekay1(0,h,isgn)
177 ELSEIF(kto.EQ.2)
THEN
181 IF(iwarm.EQ.0)
GOTO 902
189 CALL dekay2(0,h,isgn)
190 ELSEIF(kto.EQ.11)
THEN
195 CALL dekay1(1,h,isgn)
196 ELSEIF(kto.EQ.12)
THEN
201 CALL dekay2(1,h,isgn)
202 ELSEIF(kto.EQ.100)
THEN
204 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
205 CALL dadmel( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
206 CALL dadmmu( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5)
207 CALL dadmpi( 1,idum,hdum,pdum1,pdum2)
208 CALL dadmro( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4)
209 CALL dadmaa( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,pdum5,jdum)
210 CALL dadmkk( 1,idum,hdum,pdum1,pdum2)
211 CALL dadmks( 1,idum,hdum,pdum1,pdum2,pdum3,pdum4,jdum)
212 CALL dadnew( 1,idum,hdum,pdum1,pdum2,pdumx,jdum)
213 WRITE(iout,7010) nev1,nev2,nevtot
214 WRITE(iout,7011) (nevdec(i),gampmc(i),gamper(i),i= 1,7)
216 $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
227 7001
FORMAT(///1x,15(5h*****)
228 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.9 ******',9x,1h*,
229 $ /,
' *', 25x,
'***********October 2011 ***************',9x,1h*,
230 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
231 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
232 $ /,
' *', 25x,
'**AVAILABLE FROM: www.cern.ch/wasm**** ',9x,1h*,
233 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
234 $ /,
' *', 25x,
'0: Physics initialization CLEO collab ',9x,1h*,
235 $ /,
' *', 25x,
' see Alain Weinstein www home page: ',9x,1h*,
236 $ /,
' *', 25x,
'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
237 $ /,
' *', 25x,
'/korb_doc.html#files ',9x,1h*,
238 $ /,
' *', 25x,
'1: Physics initialization RChL of: ',9x,1h*,
239 $ /,
' *', 25x,
' O. Shekhovtsova, T. Przedzinski, ',9x,1h*,
240 $ /,
' *', 25x,
' P. Roig and Z. Was ',9x,1h*,
241 $ /,
' *', 25x,
' IFJPAN-2013-5, UAB-FT-731 ',9x,1h*,
242 $ /,
' *', 25x,
'*******CPC 76 (1993) 361 *****',9x,1h*,
243 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
244 $ /,
' *', 25x,
'****DEKAY ROUTINE: INITIALIZATION******',9x,1h*,
245 $ /,
' *',i20 ,5x,
'JAK1 = DECAY MODE TAU+ ',9x,1h*,
246 $ /,
' *',i20 ,5x,
'JAK2 = DECAY MODE TAU- ',9x,1h*,
248 7002
FORMAT(
' *',i20 ,5x,
'IVER = hadronic current version ',9x,1h*)
249 7010
FORMAT(///1x,15(5h*****)
250 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.9 ******',9x,1h*,
251 $ /,
' *', 25x,
'***********October 2011 ***************',9x,1h*,
252 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
253 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
254 $ /,
' *', 25x,
'**AVAILABLE FROM: www.cern.ch/wasm ****',9x,1h*,
255 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
256 $ /,
' *', 25x,
'******* 64 (1990) 275 *****',9x,1h*,
257 $ /,
' *', 25x,
'******* 70 (1992) 69 *****',9x,1h*,
258 $ /,
' *', 25x,
'******* 76 (1993) 361 *****',9x,1h*,
259 $ /,
' *', 25x,
'******* IFJPAN-2013-5, UAB-FT-731 **',9x,1h*,
260 $ /,
' *', 25x,
'*****DEKAY ROUTINE: FINAL REPORT*******',9x,1h*,
261 $ /,
' *',i20 ,5x,
'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*,
262 $ /,
' *',i20 ,5x,
'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*,
263 $ /,
' *',i20 ,5x,
'NEVTOT = SUM ',9x,1h*,
265 $
' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
267 $ ,i10,2f12.7 ,
' DADMEL ELECTRON ',9x,1h*
268 $ /,
' *',i10,2f12.7 ,
' DADMMU MUON ',9x,1h*
269 $ /,
' *',i10,2f12.7 ,
' DADMPI PION ',9x,1h*
270 $ /,
' *',i10,2f12.7,
' DADMRO RHO (->2PI) ',9x,1h*
271 $ /,
' *',i10,2f12.7,
' DADMAA A1 (->3PI) ',9x,1h*
272 $ /,
' *',i10,2f12.7,
' DADMKK KAON ',9x,1h*
273 $ /,
' *',i10,2f12.7,
' DADMKS K* ',9x,1h*)
275 $ ,i10,2f12.7,a31 ,8x,1h*)
277 $ ,20x,
'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
278 $ /,
' *',20x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
281 9020
FORMAT(
' ----- DEKAY: LACK OF INITIALISATION')
284 9100
FORMAT(
' ----- DEKAY: WRONG VALUE OF KTO ')
287 SUBROUTINE dekay1(IMOD,HH,ISGN)
290 COMMON / decp4 / pp1(4),pp2(4),kf1,kf2
291 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
292 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
293 real*4 gampmc ,gamper
295 REAL HV(4),PNU(4),PPI(4)
296 REAL PWB(4),PMU(4),PNM(4)
297 REAL PRHO(4),PIC(4),PIZ(4)
298 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
305 IF(jak1.EQ.-1)
RETURN
310 IF(jak1.EQ.0)
CALL jaker(jak)
312 CALL dadmel(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
313 ELSEIF(jak.EQ.2)
THEN
314 CALL dadmmu(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
315 ELSEIF(jak.EQ.3)
THEN
316 CALL dadmpi(0, isgn,hv,ppi,pnu)
317 ELSEIF(jak.EQ.4)
THEN
318 CALL dadmro(0, isgn,hv,pnu,prho,pic,piz)
319 ELSEIF(jak.EQ.5)
THEN
320 CALL dadmaa(0, isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
321 ELSEIF(jak.EQ.6)
THEN
322 CALL dadmkk(0, isgn,hv,pkk,pnu)
323 ELSEIF(jak.EQ.7)
THEN
324 CALL dadmks(0, isgn,hv,pnu,pks ,pkk,ppi,jkst)
326 CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
332 ELSEIF(imd.EQ.1)
THEN
336 nevdec(jak)=nevdec(jak)+1
341 CALL dwluel(1,isgn,pnu,pwb,pmu,pnm)
342 CALL dwrph(ktom,phot)
346 ELSEIF(jak.EQ.2)
THEN
347 CALL dwlumu(1,isgn,pnu,pwb,pmu,pnm)
348 CALL dwrph(ktom,phot)
352 ELSEIF(jak.EQ.3)
THEN
353 CALL dwlupi(1,isgn,ppi,pnu)
357 ELSEIF(jak.EQ.4)
THEN
358 CALL dwluro(1,isgn,pnu,prho,pic,piz)
362 ELSEIF(jak.EQ.5)
THEN
363 CALL dwluaa(1,isgn,pnu,paa,pim1,pim2,pipl,jaa)
366 ELSEIF(jak.EQ.6)
THEN
367 CALL dwlukk(1,isgn,pkk,pnu)
370 ELSEIF(jak.EQ.7)
THEN
371 CALL dwluks(1,isgn,pnu,pks,pkk,ppi,jkst)
376 CALL dwlnew(1,isgn,pnu,pwb,pnpi,jak)
384 SUBROUTINE dekay2(IMOD,HH,ISGN)
387 COMMON / decp4 / pp1(4),pp2(4),kf1,kf2
388 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
389 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
390 real*4 gampmc ,gamper
392 REAL HV(4),PNU(4),PPI(4)
393 REAL PWB(4),PMU(4),PNM(4)
394 REAL PRHO(4),PIC(4),PIZ(4)
395 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
402 IF(jak2.EQ.-1)
RETURN
407 IF(jak2.EQ.0)
CALL jaker(jak)
409 CALL dadmel(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
410 ELSEIF(jak.EQ.2)
THEN
411 CALL dadmmu(0, isgn,hv,pnu,pwb,pmu,pnm,phot)
412 ELSEIF(jak.EQ.3)
THEN
413 CALL dadmpi(0, isgn,hv,ppi,pnu)
414 ELSEIF(jak.EQ.4)
THEN
415 CALL dadmro(0, isgn,hv,pnu,prho,pic,piz)
416 ELSEIF(jak.EQ.5)
THEN
417 CALL dadmaa(0, isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
418 ELSEIF(jak.EQ.6)
THEN
419 CALL dadmkk(0, isgn,hv,pkk,pnu)
420 ELSEIF(jak.EQ.7)
THEN
421 CALL dadmks(0, isgn,hv,pnu,pks ,pkk,ppi,jkst)
423 CALL dadnew(0, isgn,hv,pnu,pwb,pnpi,jak-7)
428 ELSEIF(imd.EQ.1)
THEN
432 nevdec(jak)=nevdec(jak)+1
437 CALL dwluel(2,isgn,pnu,pwb,pmu,pnm)
438 CALL dwrph(ktom,phot)
442 ELSEIF(jak.EQ.2)
THEN
443 CALL dwlumu(2,isgn,pnu,pwb,pmu,pnm)
444 CALL dwrph(ktom,phot)
448 ELSEIF(jak.EQ.3)
THEN
449 CALL dwlupi(2,isgn,ppi,pnu)
453 ELSEIF(jak.EQ.4)
THEN
454 CALL dwluro(2,isgn,pnu,prho,pic,piz)
458 ELSEIF(jak.EQ.5)
THEN
459 CALL dwluaa(2,isgn,pnu,paa,pim1,pim2,pipl,jaa)
462 ELSEIF(jak.EQ.6)
THEN
463 CALL dwlukk(2,isgn,pkk,pnu)
466 ELSEIF(jak.EQ.7)
THEN
467 CALL dwluks(2,isgn,pnu,pks,pkk,ppi,jkst)
472 CALL dwlnew(2,isgn,pnu,pwb,pnpi,jak)
480 SUBROUTINE dexay(KTO,POL)
494 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
495 real*4 gampmc ,gamper
496 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
498 COMMON /taupos/ np1,np2
499 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
500 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
502 CHARACTER NAMES(NMODE)*31
503 COMMON / inout / inut,iout
508 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
522 WRITE(iout, 7001) jak1,jak2
523 WRITE(iout,7002) iver
527 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
528 CALL dexel(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
529 CALL dexmu(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
530 CALL dexpi(-1,idum,pdum,pdum1,pdum2)
531 CALL dexro(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4)
532 CALL dexaa(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5,idum)
533 CALL dexkk(-1,idum,pdum,pdum1,pdum2)
534 CALL dexks(-1,idum,pdum,pdum1,pdum2,pdum3,pdum4,idum)
535 CALL dexnew(-1,idum,pdum,pdum1,pdum2,pdumi,idum)
541 ELSEIF(kto.EQ.1)
THEN
546 IF(iwarm.EQ.0)
GOTO 902
549 CALL dexay1(kto,jak1,jakp,pol,isgn)
550 ELSEIF(kto.EQ.2)
THEN
555 IF(iwarm.EQ.0)
GOTO 902
556 isgn=-idff/iabs(idff)
558 CALL dexay1(kto,jak2,jakm,pol,isgn)
559 ELSEIF(kto.EQ.100)
THEN
561 IF(jak1.NE.-1.OR.jak2.NE.-1)
THEN
562 CALL dexel( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
563 CALL dexmu( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5)
564 CALL dexpi( 1,idum,pdum,pdum1,pdum2)
565 CALL dexro( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4)
566 CALL dexaa( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,pdum5,idum)
567 CALL dexkk( 1,idum,pdum,pdum1,pdum2)
568 CALL dexks( 1,idum,pdum,pdum1,pdum2,pdum3,pdum4,idum)
569 CALL dexnew( 1,idum,pdum,pdum1,pdum2,pdumi,idum)
570 WRITE(iout,7010) nev1,nev2,nevtot
571 WRITE(iout,7011) (nevdec(i),gampmc(i),gamper(i),i= 1,7)
573 $ (nevdec(i),gampmc(i),gamper(i),names(i-7),i=8,7+nmode)
580 7001
FORMAT(///1x,15(5h*****)
581 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.9 ******',9x,1h*,
582 $ /,
' *', 25x,
'***********October 2011***************',9x,1h*,
583 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
584 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
585 $ /,
' *', 25x,
'**AVAILABLE FROM: www.cern.ch/wasm****',9x,1h*,
586 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
587 $ /,
' *', 25x,
'0: Physics initialization CLEO collab ',9x,1h*,
588 $ /,
' *', 25x,
' see Alain Weinstein www home page: ',9x,1h*,
589 $ /,
' *', 25x,
'http://www.cithep.caltech.edu/~ajw/ ',9x,1h*,
590 $ /,
' *', 25x,
'/korb_doc.html#files ',9x,1h*,
591 $ /,
' *', 25x,
'1: Physics initialization RChL of: ',9x,1h*,
592 $ /,
' *', 25x,
' O. Shekhovtsova, T. Przedzinski, ',9x,1h*,
593 $ /,
' *', 25x,
' P. Roig and Z. Was ',9x,1h*,
594 $ /,
' *', 25x,
' IFJPAN-2013-5, UAB-FT-731 ',9x,1h*,
595 $ /,
' *', 25x,
'****** CPC 76 (1993) 361 ******',9x,1h*,
596 $ /,
' *', 25x,
'**5 or more pi dec.: precision limited ',9x,1h*,
597 $ /,
' *', 25x,
'******DEXAY ROUTINE: INITIALIZATION****',9x,1h*
598 $ /,
' *',i20 ,5x,
'JAK1 = DECAY MODE FERMION1 (TAU+) ',9x,1h*
599 $ /,
' *',i20 ,5x,
'JAK2 = DECAY MODE FERMION2 (TAU-) ',9x,1h*
601 7002
FORMAT(
' *',i20 ,5x,
'IVER = hadronic current version ',9x,1h*)
604 7010
FORMAT(///1x,15(5h*****)
605 $ /,
' *', 25x,
'*****TAUOLA LIBRARY: VERSION 2.9 ******',9x,1h*,
606 $ /,
' *', 25x,
'***********October 2011***************',9x,1h*
607 $ /,
' *', 25x,
'**AUTHORS: S.JADACH, Z.WAS*************',9x,1h*,
608 $ /,
' *', 25x,
'**R. DECKER, M. JEZABEK, J.H.KUEHN*****',9x,1h*,
609 $ /,
' *', 25x,
'**AVAILABLE FROM: www.cern.ch/wasm ****',9x,1h*,
610 $ /,
' *', 25x,
'***** PUBLISHED IN COMP. PHYS. COMM.***',9x,1h*,
611 $ /,
' *', 25x,
'******* 64 (1990) 275 *****',9x,1h*,
612 $ /,
' *', 25x,
'******* 70 (1992) 69 *****',9x,1h*,
613 $ /,
' *', 25x,
'******* 76 (1993) 361 *****',9x,1h*,
614 $ /,
' *', 25x,
'******* IFJPAN-2013-5, UAB-FT-731 **',9x,1h*,
615 $ /,
' *', 25x,
'******DEXAY ROUTINE: FINAL REPORT******',9x,1h*
616 $ /,
' *',i20 ,5x,
'NEV1 = NO. OF TAU+ DECS. ACCEPTED ',9x,1h*
617 $ /,
' *',i20 ,5x,
'NEV2 = NO. OF TAU- DECS. ACCEPTED ',9x,1h*
618 $ /,
' *',i20 ,5x,
'NEVTOT = SUM ',9x,1h*
620 $
' PART.WIDTH ERROR ROUTINE DECAY MODE ',9x,1h*)
622 $ ,i10,2f12.7 ,
' DADMEL ELECTRON ',9x,1h*
623 $ /,
' *',i10,2f12.7 ,
' DADMMU MUON ',9x,1h*
624 $ /,
' *',i10,2f12.7 ,
' DADMPI PION ',9x,1h*
625 $ /,
' *',i10,2f12.7,
' DADMRO RHO (->2PI) ',9x,1h*
626 $ /,
' *',i10,2f12.7,
' DADMAA A1 (->3PI) ',9x,1h*
627 $ /,
' *',i10,2f12.7,
' DADMKK KAON ',9x,1h*
628 $ /,
' *',i10,2f12.7,
' DADMKS K* ',9x,1h*)
630 $ ,i10,2f12.7,a31 ,8x,1h*)
632 $ ,20x,
'THE ERROR IS RELATIVE AND PART.WIDTH ',10x,1h*
633 $ /,
' *',20x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',10x,1h*
635 902
WRITE(iout, 9020)
636 9020
FORMAT(
' ----- DEXAY: LACK OF INITIALISATION')
638 910
WRITE(iout, 9100)
639 9100
FORMAT(
' ----- DEXAY: WRONG VALUE OF KTO ')
642 SUBROUTINE dexay1(KTO,JAKIN,JAK,POL,ISGN)
648 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
649 real*4 gampmc ,gamper
650 COMMON / inout / inut,iout
653 REAL PRHO(4),PIC(4),PIZ(4)
654 REAL PWB(4),PMU(4),PNM(4)
655 REAL PAA(4),PIM1(4),PIM2(4),PIPL(4)
661 IF(jakin.EQ.-1)
RETURN
668 IF(jak.EQ.0)
CALL jaker(jak)
671 CALL dexel(0, isgn,polar,pnu,pwb,pmu,pnm,phot)
672 CALL dwluel(kto,isgn,pnu,pwb,pmu,pnm)
673 CALL dwrph(kto,phot )
674 ELSEIF(jak.EQ.2)
THEN
675 CALL dexmu(0, isgn,polar,pnu,pwb,pmu,pnm,phot)
676 CALL dwlumu(kto,isgn,pnu,pwb,pmu,pnm)
677 CALL dwrph(kto,phot )
678 ELSEIF(jak.EQ.3)
THEN
679 CALL dexpi(0, isgn,polar,ppi,pnu)
680 CALL dwlupi(kto,isgn,ppi,pnu)
681 ELSEIF(jak.EQ.4)
THEN
682 CALL dexro(0, isgn,polar,pnu,prho,pic,piz)
683 CALL dwluro(kto,isgn,pnu,prho,pic,piz)
684 ELSEIF(jak.EQ.5)
THEN
685 CALL dexaa(0, isgn,polar,pnu,paa,pim1,pim2,pipl,jaa)
686 CALL dwluaa(kto,isgn,pnu,paa,pim1,pim2,pipl,jaa)
687 ELSEIF(jak.EQ.6)
THEN
688 CALL dexkk(0, isgn,polar,pkk,pnu)
689 CALL dwlukk(kto,isgn,pkk,pnu)
690 ELSEIF(jak.EQ.7)
THEN
691 CALL dexks(0, isgn,polar,pnu,pks,pkk,ppi,jkst)
692 CALL dwluks(kto,isgn,pnu,pks,pkk,ppi,jkst)
695 CALL dexnew(0, isgn,polar,pnu,pwb,pnpi,jnpi)
696 CALL dwlnew(kto,isgn,pnu,pwb,pnpi,jak)
698 nevdec(jak)=nevdec(jak)+1
700 SUBROUTINE dexel(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
707 REAL POL(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4),PH(4),RN(1)
713 CALL dadmel( -1,isgn,hv,pnu,pwb,q1,q2,ph)
716 ELSEIF(mode.EQ. 0)
THEN
719 IF(iwarm.EQ.0)
GOTO 902
720 CALL dadmel( 0,isgn,hv,pnu,pwb,q1,q2,ph)
721 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
724 IF(rn(1).GT.wt)
GOTO 300
726 ELSEIF(mode.EQ. 1)
THEN
728 CALL dadmel( 1,isgn,hv,pnu,pwb,q1,q2,ph)
734 9020
FORMAT(
' ----- DEXEL: LACK OF INITIALISATION')
737 SUBROUTINE dexmu(MODE,ISGN,POL,PNU,PWB,Q1,Q2,PH)
746 COMMON / inout / inut,iout
747 REAL POL(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4),PH(4),RN(1)
753 CALL dadmmu( -1,isgn,hv,pnu,pwb,q1,q2,ph)
756 ELSEIF(mode.EQ. 0)
THEN
759 IF(iwarm.EQ.0)
GOTO 902
760 CALL dadmmu( 0,isgn,hv,pnu,pwb,q1,q2,ph)
761 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
764 IF(rn(1).GT.wt)
GOTO 300
766 ELSEIF(mode.EQ. 1)
THEN
768 CALL dadmmu( 1,isgn,hv,pnu,pwb,q1,q2,ph)
773 902
WRITE(iout, 9020)
774 9020
FORMAT(
' ----- DEXMU: LACK OF INITIALISATION')
777 SUBROUTINE dadmel(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
782 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
783 real*4 gfermi,gv,ga,ccabib,scabib,gamel
784 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
785 * ,ampiz,ampi,amro,gamro,ama1,gama1
786 * ,amk,amkz,amkst,gamkst
788 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
789 * ,ampiz,ampi,amro,gamro,ama1,gama1
790 * ,amk,amkz,amkst,gamkst
791 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
792 real*4 gampmc ,gamper
794 COMMON / inout / inut,iout
795 REAL HHV(4),HV(4),PWB(4),PNU(4),Q1(4),Q2(4)
796 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
799 DATA pi /3.141592653589793238462643/
812 CALL dphsel(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
813 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
817 ELSEIF(mode.EQ. 0)
THEN
820 IF(iwarm.EQ.0)
GOTO 902
822 CALL dphsel(wt,hv,pnu,pwb,q1,q2,phx)
828 IF(wt.GT.wtmax) nevovr=nevovr+1
829 IF(rn*wtmax.GT.wt)
GOTO 300
836 CALL rotor2(thet,pnu,pnu)
837 CALL rotor3( phi,pnu,pnu)
838 CALL rotor2(thet,pwb,pwb)
839 CALL rotor3( phi,pwb,pwb)
840 CALL rotor2(thet,q1,q1)
841 CALL rotor3( phi,q1,q1)
842 CALL rotor2(thet,q2,q2)
843 CALL rotor3( phi,q2,q2)
844 CALL rotor2(thet,hv,hv)
845 CALL rotor3( phi,hv,hv)
846 CALL rotor2(thet,phx,phx)
847 CALL rotor3( phi,phx,phx)
849 44 hhv(i)=-isgn*hv(i)
852 ELSEIF(mode.EQ. 1)
THEN
854 IF(nevraw.EQ.0)
RETURN
855 pargam=swt/float(nevraw+1)
857 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
859 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
867 7010
FORMAT(///1x,15(5h*****)
868 $ /,
' *', 25x,
'******** DADMEL FINAL REPORT ******** ',9x,1h*
869 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF EL DECAYS TOTAL ',9x,1h*
870 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF EL DECS. ACCEPTED ',9x,1h*
871 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
872 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( ELECTRON) IN GEV UNITS ',9x,1h*
873 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
874 $ /,
' *',f20.9,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
875 $ /,
' *',25x,
'COMPLETE QED CORRECTIONS INCLUDED ',9x,1h*
876 $ /,
' *',25x,
'BUT ONLY V-A CUPLINGS ',9x,1h*
878 902
WRITE(iout, 9020)
879 9020
FORMAT(
' ----- DADMEL: LACK OF INITIALISATION')
882 SUBROUTINE dadmmu(MODE,ISGN,HHV,PNU,PWB,Q1,Q2,PHX)
884 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
885 * ,ampiz,ampi,amro,gamro,ama1,gama1
886 * ,amk,amkz,amkst,gamkst
888 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
889 * ,ampiz,ampi,amro,gamro,ama1,gama1
890 * ,amk,amkz,amkst,gamkst
891 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
892 real*4 gfermi,gv,ga,ccabib,scabib,gamel
893 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
894 real*4 gampmc ,gamper
895 COMMON / inout / inut,iout
897 REAL HHV(4),HV(4),PNU(4),PWB(4),Q1(4),Q2(4)
898 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
901 DATA pi /3.141592653589793238462643/
914 CALL dphsmu(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5)
915 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
919 ELSEIF(mode.EQ. 0)
THEN
922 IF(iwarm.EQ.0)
GOTO 902
924 CALL dphsmu(wt,hv,pnu,pwb,q1,q2,phx)
930 IF(wt.GT.wtmax) nevovr=nevovr+1
931 IF(rn*wtmax.GT.wt)
GOTO 300
936 CALL rotor2(thet,pnu,pnu)
937 CALL rotor3( phi,pnu,pnu)
938 CALL rotor2(thet,pwb,pwb)
939 CALL rotor3( phi,pwb,pwb)
940 CALL rotor2(thet,q1,q1)
941 CALL rotor3( phi,q1,q1)
942 CALL rotor2(thet,q2,q2)
943 CALL rotor3( phi,q2,q2)
944 CALL rotor2(thet,hv,hv)
945 CALL rotor3( phi,hv,hv)
946 CALL rotor2(thet,phx,phx)
947 CALL rotor3( phi,phx,phx)
949 44 hhv(i)=-isgn*hv(i)
952 ELSEIF(mode.EQ. 1)
THEN
954 IF(nevraw.EQ.0)
RETURN
955 pargam=swt/float(nevraw+1)
957 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
959 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
967 7010
FORMAT(///1x,15(5h*****)
968 $ /,
' *', 25x,
'******** DADMMU FINAL REPORT ******** ',9x,1h*
969 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF MU DECAYS TOTAL ',9x,1h*
970 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF MU DECS. ACCEPTED ',9x,1h*
971 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
972 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (MU DECAY) IN GEV UNITS ',9x,1h*
973 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
974 $ /,
' *',f20.9,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
975 $ /,
' *',25x,
'COMPLETE QED CORRECTIONS INCLUDED ',9x,1h*
976 $ /,
' *',25x,
'BUT ONLY V-A CUPLINGS ',9x,1h*
978 902
WRITE(iout, 9020)
979 9020
FORMAT(
' ----- DADMMU: LACK OF INITIALISATION')
982 SUBROUTINE dphsel(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
988 real*4 hvx(4),paax(4),xax(4),qpx(4),xnx(4)
989 real*8 hv(4),ph(4),paa(4),xa(4),qp(4),xn(4)
992 CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
1003 SUBROUTINE dphsmu(DGAMX,HVX,XNX,PAAX,QPX,XAX,PHX)
1009 real*4 hvx(4),paax(4),xax(4),qpx(4),xnx(4)
1010 real*8 hv(4),ph(4),paa(4),xa(4),qp(4),xn(4)
1013 CALL drcmu(dgamt,hv,ph,paa,xa,qp,xn,ielmu)
1024 SUBROUTINE drcmu(DGAMT,HV,PH,PAA,XA,QP,XN,IELMU)
1025 IMPLICIT REAL*8 (a-h,o-z)
1030 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1031 * ,ampiz,ampi,amro,gamro,ama1,gama1
1032 * ,amk,amkz,amkst,gamkst
1034 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1035 * ,ampiz,ampi,amro,gamro,ama1,gama1
1036 * ,amk,amkz,amkst,gamkst
1037 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1038 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1039 COMMON / inout / inut,iout
1040 COMMON / taurad / xk0dec,itdkrc
1042 real*8 hv(4),pt(4),ph(4),paa(4),xa(4),qp(4),xn(4)
1046 DATA pi /3.141592653589793238462643d0/
1052 phspac=1./2**17/pi**8
1062 IF (ielmu.EQ.1)
THEN
1069 IF ( itdkrc.EQ.0) prhard=0d0
1071 IF(prsoft.LT.0.1)
THEN
1072 print *,
'ERROR IN DRCMU; PRSOFT=',prsoft
1077 ihard=(rr5.GT.prsoft)
1081 ams1=(amu+amnuta)**2
1084 xl1=log(xk1/2/xk0dec)
1089 phspac=phspac*ams2*xl1*xk
1090 phspac=phspac/prhard
1093 phspac=phspac*2**6*pi**3
1094 phspac=phspac/prsoft
1103 am2sq=ams1+ rr2*(ams2-ams1)
1105 phspac=phspac*(ams2-ams1)
1107 enq1=(am2sq+amnuta**2)/(2*am2)
1108 enq2=(am2sq-amnuta**2)/(2*am2)
1109 ppi= enq1**2-amnuta**2
1110 pppi=sqrt(abs(enq1**2-amnuta**2))
1111 phspac=phspac*(4*pi)*(2*pppi/am2)
1113 CALL spherd(pppi,xn)
1123 pr(4)=1.d0/(2*am3)*(am3**2+am2**2-amu**2)
1124 pr(3)= sqrt(abs(pr(4)**2-am2**2))
1125 ppi = pr(4)**2-am2**2
1129 qp(4)=1.d0/(2*am3)*(am3**2-am2**2+amu**2)
1131 phspac=phspac*(4*pi)*(2*pr(3)/am3)
1133 exe=(pr(4)+pr(3))/am2
1134 CALL bostd3(exe,xn,xn)
1135 CALL bostd3(exe,xa,xa)
1139 eps=4*(amu/amtax)**2
1140 xl1=log((2+eps)/eps)
1142 eta =exp(xl1*rr3+xl0)
1145 phspac=phspac*xl1/2*eta
1147 CALL rotpox(thet,phi,xn)
1148 CALL rotpox(thet,phi,xa)
1149 CALL rotpox(thet,phi,qp)
1150 CALL rotpox(thet,phi,pr)
1156 paa(4)=1/(2*amtax)*(amtax**2+am3**2)
1157 paa(3)= sqrt(abs(paa(4)**2-am3**2))
1158 ppi = paa(4)**2-am3**2
1159 phspac=phspac*(4*pi)*(2*paa(3)/amtax)
1167 exe=(paa(4)+paa(3))/am3
1168 CALL bostd3(exe,xn,xn)
1169 CALL bostd3(exe,xa,xa)
1170 CALL bostd3(exe,qp,qp)
1171 CALL bostd3(exe,pr,pr)
1173 thet =acos(-1.+2*rr3)
1175 CALL rotpox(thet,phi,xn)
1176 CALL rotpox(thet,phi,xa)
1177 CALL rotpox(thet,phi,qp)
1178 CALL rotpox(thet,phi,pr)
1193 CALL dampry(itdkrc,xk0dec,ph,xa,qp,xn,amplit,hv)
1194 dgamt=1/(2.*amtax)*amplit*phspac
1196 SUBROUTINE dampry(ITDKRC,XK0DEC,XK,XA,QP,XN,AMPLIT,HV)
1197 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 real*8 hv(4),qp(4),xn(4),xa(4),xk(4)
1214 IF(xk(4).LT.0.1d0*ak0)
THEN
1215 amplit=thb(itdkrc,qp,xn,xa,ak0,hv)
1217 amplit=sqm2(itdkrc,qp,xn,xa,xk,ak0,hv)
1221 FUNCTION sqm2(ITDKRC,QP,XN,XA,XK,AK0,HV)
1234 IMPLICIT REAL*8(a-h,o-z)
1235 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1236 * ,ampiz,ampi,amro,gamro,ama1,gama1
1237 * ,amk,amkz,amkst,gamkst
1239 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1240 * ,ampiz,ampi,amro,gamro,ama1,gama1
1241 * ,amk,amkz,amkst,gamkst
1242 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1243 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1244 COMMON / qedprm /alfinv,alfpi,xk0
1245 real*8 alfinv,alfpi,xk0
1246 real*8 qp(4),xn(4),xa(4),xk(4)
1249 real*8 s0(3),rxa(3),rxk(3),rqp(3)
1250 DATA pi /3.141592653589793238462643d0/
1256 emass2=qp(4)**2-qp(1)**2-qp(2)**2-qp(3)**2
1264 rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1266 rxk(i)=r(4)*xk(4)-r(1)*xk(1)-r(2)*xk(2)-r(3)*xk(3)
1267 rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
1269 qpxn=qp(4)*xn(4)-qp(1)*xn(1)-qp(2)*xn(2)-qp(3)*xn(3)
1270 qpxa=qp(4)*xa(4)-qp(1)*xa(1)-qp(2)*xa(2)-qp(3)*xa(3)
1271 qpxk=qp(4)*xk(4)-qp(1)*xk(1)-qp(2)*xk(2)-qp(3)*xk(3)
1273 xnxk=xn(4)*xk(4)-xn(1)*xk(1)-xn(2)*xk(2)-xn(3)*xk(3)
1274 xaxk=xa(4)*xk(4)-xa(1)*xk(1)-xa(2)*xk(2)-xa(3)*xk(3)
1284 s1= qpxn*txa*( -emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
1286 $qpxn/txk**2* ( tmass2*xaxk - txa*txk+ xaxk*txk) -
1287 $txa*txn/txk - qpxn/(qpxk*txk)* (tqp*xaxk-txk*qpxa)
1288 const4=256*pi/alphai*gf**2
1289 IF (itdkrc.EQ.0) const4=0d0
1292 s0(i) = qpxn*rxa(i)*(-emass2/qpxk**2*a + 2*tqp/(qpxk*txk)*b-
1294 $ qpxn/txk**2* (tmass2*xaxk - txa*rxk(i)+ xaxk*rxk(i))-
1295 $ rxa(i)*txn/txk - qpxn/(qpxk*txk)*(rqp(i)*xaxk- rxk(i)*qpxa)
1296 5 hv(i)=s0(i)/s1-1.d0
1299 FUNCTION thb(ITDKRC,QP,XN,XA,AK0,HV)
1313 IMPLICIT REAL*8(a-h,o-z)
1314 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1315 * ,ampiz,ampi,amro,gamro,ama1,gama1
1316 * ,amk,amkz,amkst,gamkst
1318 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1319 * ,ampiz,ampi,amro,gamro,ama1,gama1
1320 * ,amk,amkz,amkst,gamkst
1321 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1322 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1323 COMMON / qedprm /alfinv,alfpi,xk0
1324 real*8 alfinv,alfpi,xk0
1325 dimension qp(4),xn(4),xa(4)
1328 real*8 rxa(3),rxn(3),rqp(3)
1329 real*8 bornpl(3),am3pol(3),xm3pol(3)
1330 DATA pi /3.141592653589793238462643d0/
1343 rxa(i)=r(4)*xa(4)-r(1)*xa(1)-r(2)*xa(2)-r(3)*xa(3)
1344 rxn(i)=r(4)*xn(4)-r(1)*xn(1)-r(2)*xn(2)-r(3)*xn(3)
1346 rqp(i)=r(4)*qp(4)-r(1)*qp(1)-r(2)*qp(2)-r(3)*qp(3)
1350 u3=sqrt(qp(1)**2+qp(2)**2+qp(3)**2)/tmass
1352 w0=(xn(4)+xa(4))/tmass
1364 f0=2*u0/u3*( dilogt(1-(um*wm/(up*wp)))- dilogt(1-wm/wp) +
1365 $dilogt(1-um/up) -2*yu+ 2*log(up)*(yw+yu) ) +
1366 $1/y* ( 2*u3*yu + (1-eps2- 2*y)*log(eps) ) +
1367 $ 2 - 4*(u0/u3*yu -1)* log(2*al)
1368 fp= yu/(2*u3)*(1 + (1-eps2)/y ) + log(eps)/y
1369 fm= yu/(2*u3)*(1 - (1-eps2)/y ) - log(eps)/y
1372 qpxn=qp(4)*xn(4)-qp(1)*xn(1)-qp(2)*xn(2)-qp(3)*xn(3)
1373 qpxa=qp(4)*xa(4)-qp(1)*xa(1)-qp(2)*xa(2)-qp(3)*xa(3)
1374 xnxa=xn(4)*xa(4)-xn(1)*xa(1)-xn(2)*xa(2)-xn(3)*xa(3)
1379 const3=1/(2*alphai*pi)*64*gf**2
1380 IF (itdkrc.EQ.0) const3=0d0
1381 xm3= -( f0* qpxn*txa + fp*eps2* txn*txa +
1382 $fm* qpxn*qpxa + f3* tmass2*xnxa )
1385 brak= (gv+ga)**2*tqp*xnxa+(gv-ga)**2*txa*qpxn
1386 & -(gv**2-ga**2)*tmass*amnuta*qpxa
1387 born= 32*(gfermi**2/2.)*brak
1389 xm3pol(i)= -( f0* qpxn*rxa(i) + fp*eps2* txn*rxa(i) +
1390 $ fm* qpxn* (qpxa + (rxa(i)*tqp-txa*rqp(i))/tmass2 ) +
1391 $ f3* (tmass2*xnxa +txn*rxa(i) -rxn(i)*txa) )
1392 am3pol(i)=xm3pol(i)*const3
1395 & (gv+ga)**2*tmass*xnxa*qp(i)
1396 & -(gv-ga)**2*tmass*qpxn*xa(i)
1397 & +(gv**2-ga**2)*amnuta*txa*qp(i)
1398 & -(gv**2-ga**2)*amnuta*tqp*xa(i) )*
1400 5 hv(i)=(bornpl(i)+am3pol(i))/(born+am3)-1.d0
1402 IF (thb/born.LT.0.1d0)
THEN
1403 print *,
'ERROR IN THB, THB/BORN=',thb/born
1408 SUBROUTINE dexpi(MODE,ISGN,POL,PPI,PNU)
1415 REAL POL(4),HV(4),PNU(4),PPI(4),RN(1)
1419 CALL dadmpi(-1,isgn,hv,ppi,pnu)
1422 ELSEIF(mode.EQ. 0)
THEN
1425 CALL dadmpi( 0,isgn,hv,ppi,pnu)
1426 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1429 IF(rn(1).GT.wt)
GOTO 300
1431 ELSEIF(mode.EQ. 1)
THEN
1433 CALL dadmpi( 1,isgn,hv,ppi,pnu)
1439 SUBROUTINE dadmpi(MODE,ISGN,HV,PPI,PNU)
1441 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1442 * ,ampiz,ampi,amro,gamro,ama1,gama1
1443 * ,amk,amkz,amkst,gamkst
1445 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1446 * ,ampiz,ampi,amro,gamro,ama1,gama1
1447 * ,amk,amkz,amkst,gamkst
1448 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1449 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1450 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1451 real*4 gampmc ,gamper
1452 COMMON / inout / inut,iout
1453 REAL PPI(4),PNU(4),HV(4)
1454 DATA pi /3.141592653589793238462643/
1459 ELSEIF(mode.EQ. 0)
THEN
1462 epi= (amtau**2+ampi**2-amnuta**2)/(2*amtau)
1463 enu= (amtau**2-ampi**2+amnuta**2)/(2*amtau)
1464 xpi= sqrt(epi**2-ampi**2)
1466 CALL sphera(xpi,ppi)
1474 qxn=ppi(4)*pnu(4)-ppi(1)*pnu(1)-ppi(2)*pnu(2)-ppi(3)*pnu(3)
1475 brak=(gv**2+ga**2)*(2*pxq*qxn-ampi**2*pxn)
1476 & +(gv**2-ga**2)*amtau*amnuta*ampi**2
1478 40 hv(i)=-isgn*2*ga*gv*amtau*(2*ppi(i)*qxn-pnu(i)*ampi**2)/brak
1481 ELSEIF(mode.EQ. 1)
THEN
1483 IF(nevtot.EQ.0)
RETURN
1489 gamm=(gfermi*fpi)**2/(16.*pi)*amtau**3*
1491 $ sqrt((amtau**2-ampi**2-amnuta**2)**2
1492 $ -4*ampi**2*amnuta**2 )/amtau**2
1495 WRITE(iout, 7010) nevtot,gamm,rat,error
1502 7010
FORMAT(///1x,15(5h*****)
1503 $ /,
' *', 25x,
'******** DADMPI FINAL REPORT ******** ',9x,1h*
1504 $ /,
' *',i20 ,5x,
'NEVTOT = NO. OF PI DECAYS TOTAL ',9x,1h*
1505 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( PI DECAY) IN GEV UNITS ',9x,1h*
1506 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1507 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9x,1h*
1508 $ /,1x,15(5h*****)/)
1510 SUBROUTINE dexro(MODE,ISGN,POL,PNU,PRO,PIC,PIZ)
1519 COMMON / inout / inut,iout
1520 REAL POL(4),HV(4),PRO(4),PNU(4),PIC(4),PIZ(4),RN(1)
1526 CALL dadmro( -1,isgn,hv,pnu,pro,pic,piz)
1530 ELSEIF(mode.EQ. 0)
THEN
1533 IF(iwarm.EQ.0)
GOTO 902
1534 CALL dadmro( 0,isgn,hv,pnu,pro,pic,piz)
1535 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
1540 IF(rn(1).GT.wt)
GOTO 300
1542 ELSEIF(mode.EQ. 1)
THEN
1544 CALL dadmro( 1,isgn,hv,pnu,pro,pic,piz)
1550 902
WRITE(iout, 9020)
1551 9020
FORMAT(
' ----- DEXRO: LACK OF INITIALISATION')
1554 SUBROUTINE dadmro(MODE,ISGN,HHV,PNU,PRO,PIC,PIZ)
1556 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1557 * ,ampiz,ampi,amro,gamro,ama1,gama1
1558 * ,amk,amkz,amkst,gamkst
1560 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1561 * ,ampiz,ampi,amro,gamro,ama1,gama1
1562 * ,amk,amkz,amkst,gamkst
1563 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1564 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1565 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
1566 real*4 gampmc ,gamper
1567 COMMON / inout / inut,iout
1569 REAL HV(4),PRO(4),PNU(4),PIC(4),PIZ(4)
1570 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4)
1573 DATA pi /3.141592653589793238462643/
1586 CALL dphsro(wt,hv,pdum1,pdum2,pdum3,pdum4)
1587 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
1592 ELSEIF(mode.EQ. 0)
THEN
1595 IF(iwarm.EQ.0)
GOTO 902
1596 CALL dphsro(wt,hv,pnu,pro,pic,piz)
1603 IF(wt.GT.wtmax) nevovr=nevovr+1
1604 IF(rn*wtmax.GT.wt)
GOTO 300
1606 costhe=-1.+2.*rrr(2)
1609 CALL rotor2(thet,pnu,pnu)
1610 CALL rotor3( phi,pnu,pnu)
1611 CALL rotor2(thet,pro,pro)
1612 CALL rotor3( phi,pro,pro)
1613 CALL rotor2(thet,pic,pic)
1614 CALL rotor3( phi,pic,pic)
1615 CALL rotor2(thet,piz,piz)
1616 CALL rotor3( phi,piz,piz)
1617 CALL rotor2(thet,hv,hv)
1618 CALL rotor3( phi,hv,hv)
1620 44 hhv(i)=-isgn*hv(i)
1623 ELSEIF(mode.EQ. 1)
THEN
1625 IF(nevraw.EQ.0)
RETURN
1626 pargam=swt/float(nevraw+1)
1628 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
1630 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
1638 7003
FORMAT(///1x,15(5h*****)
1639 $ /,
' *', 25x,
'******** DADMRO INITIALISATION ********',9x,1h*
1640 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
1641 $ /,1x,15(5h*****)/)
1642 7010
FORMAT(///1x,15(5h*****)
1643 $ /,
' *', 25x,
'******** DADMRO FINAL REPORT ******** ',9x,1h*
1644 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF RHO DECAYS TOTAL ',9x,1h*
1645 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF RHO DECS. ACCEPTED ',9x,1h*
1646 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
1647 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (RHO DECAY) IN GEV UNITS ',9x,1h*
1648 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
1649 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
1650 $ /,1x,15(5h*****)/)
1651 902
WRITE(iout, 9020)
1652 9020
FORMAT(
' ----- DADMRO: LACK OF INITIALISATION')
1655 SUBROUTINE dphsro(DGAMT,HV,PN,PR,PIC,PIZ)
1660 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1661 * ,ampiz,ampi,amro,gamro,ama1,gama1
1662 * ,amk,amkz,amkst,gamkst
1664 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1665 * ,ampiz,ampi,amro,gamro,ama1,gama1
1666 * ,amk,amkz,amkst,gamkst
1667 REAL HV(4),PT(4),PN(4),PR(4),PIC(4),PIZ(4),QQ(4),RR1(1)
1668 DATA pi /3.141592653589793238462643/
1672 phspac=1./2**11/pi**5
1679 ams1=(ampi+ampiz)**2
1680 ams2=(amtau-amnuta)**2
1686 alp1=atan((ams1-amro**2)/amro/gamro)
1687 alp2=atan((ams2-amro**2)/amro/gamro)
1691 alp=alp1+rr1(1)*(alp2-alp1)
1692 amx2=amro**2+amro*gamro*tan(alp)
1694 IF(amx.LT.2.*ampi)
GO TO 100
1696 phspac=phspac*((amx2-amro**2)**2+(amro*gamro)**2)/(amro*gamro)
1697 phspac=phspac*(alp2-alp1)
1702 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
1703 pn(3)=-sqrt(abs((pn(4)-amnuta)*(pn(4)+amnuta)))
1707 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
1709 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
1712 enq1=(amx2+ampi**2-ampiz**2)/(2.*amx)
1713 enq2=(amx2-ampi**2+ampiz**2)/(2.*amx)
1714 pppi=sqrt((enq1-ampi)*(enq1+ampi))
1715 phspac=phspac*(4*pi)*(2*pppi/amx)
1717 CALL sphera(pppi,pic)
1723 exe=(pr(4)+pr(3))/amx
1725 CALL bostr3(exe,pic,pic)
1726 CALL bostr3(exe,piz,piz)
1728 CALL dam2pi(0,pt,pn,pic,piz,amplit,hv)
1729 dgamt=1/(2.*amtau)*amplit*phspac
1744 SUBROUTINE dam2pi(MNUM,PT,PN,PIM1,PIM2,AMPLIT,HV)
1754 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1755 * ,ampiz,ampi,amro,gamro,ama1,gama1
1756 * ,amk,amkz,amkst,gamkst
1758 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1759 * ,ampiz,ampi,amro,gamro,ama1,gama1
1760 * ,amk,amkz,amkst,gamkst
1761 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
1762 real*4 gfermi,gv,ga,ccabib,scabib,gamel
1763 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4)
1764 REAL PIVEC(4),PIAKS(4),HVM(4)
1766 DATA pi /3.141592653589793238462643/
1770 CALL curr_pipi0(pim1,pim2,hadcur)
1771 ELSEIF (mnum.EQ.1)
THEN
1772 CALL curr_pik0(pim1,pim2,hadcur)
1773 ELSEIF (mnum.EQ.2)
THEN
1774 CALL curr_kpi0(pim1,pim2,hadcur)
1775 ELSEIF (mnum.EQ.3)
THEN
1776 CALL curr_kk0(pim1,pim2,hadcur)
1778 write(*,*)
'DAM2PI: wrong MNUM= ',mnum
1784 CALL clvec(hadcur,pn,pivec)
1785 CALL claxi(hadcur,pn,piaks)
1786 CALL clnut(hadcur,brakm,hvm)
1788 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
1789 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
1790 IF (mnum.EQ.0.OR.mnum.EQ.3)
THEN
1791 amplit=(ccabib*gfermi)**2*brak
1793 amplit=(scabib*gfermi)**2*brak
1797 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
1798 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
1805 SUBROUTINE curr_pipi0(PC,PN,HADCUR)
1814 COMPLEX BWIGS,HADCUR(4),FKPIPL,FRHO_PI
1816 REAL PC(4),PN(4),QQ(4),PKS(4),FPIRHO
1817 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1818 * ,ampiz,ampi,amro,gamro,ama1,gama1
1819 * ,amk,amkz,amkst,gamkst
1821 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1822 * ,ampiz,ampi,amro,gamro,ama1,gama1
1823 * ,amk,amkz,amkst,gamkst
1832 pks(ik)=pc(ik)+ pn(ik)
1833 qq(ik)=pc(ik)- pn(ik)
1836 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
1837 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
1839 31 qq(ik)=qq(ik)-pks(ik)*qqpks/pksd
1842 CALL getff2pirho(ff2pirho)
1843 IF (ff2pirho.EQ.2)
THEN
1846 hadcur(k)=qq(k)* fpibel(sqrt(pksd),0)
1848 ELSEIF (ff2pirho.EQ.3)
THEN
1853 hadcur(k)=qq(k)* fpibel(sqrt(pksd),1)
1856 write(*,*)
'problem in 2-scalars current FF2PIRHO=',ff2pirho
1859 ELSEIF (iver.EQ.0)
THEN
1861 hadcur(k)=qq(k)* sqrt(fpirho(sqrt(pksd)))
1865 write(*,*)
'problem in 2-scalars current IVER=',iver
1872 SUBROUTINE curr_pik0(PC,PN,HADCUR)
1881 COMPLEX BWIGS,HADCUR(4),FKPIPL
1882 REAL PC(4),PN(4),QQ(4),PKS(4),FKPISC,PKSD,QQPKS
1883 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1884 & ,ampiz,ampi,amro,gamro,ama1,gama1
1885 & ,amk,amkz,amkst,gamkst
1887 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1888 & ,ampiz,ampi,amro,gamro,ama1,gama1
1889 & ,amk,amkz,amkst,gamkst,fact_k0pi
1890 COMMON / taukle / bra1,brk0,brk0b,brks
1891 real*4 bra1,brk0,brk0b,brks
1900 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
1901 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
1903 31 qq(i)=qq(i)-pks(i)*qqpks/pksd
1907 hadcur(k)=qq(k)*bwigs(pksd,amkst,gamkst)
1913 SUBROUTINE curr_kpi0(PC,PN,HADCUR)
1922 COMPLEX BWIGS,HADCUR(4),FKPIPL
1923 REAL PC(4),PN(4),QQ(4),PKS(4),FKPISC,PKSD,QQPKS
1924 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1925 * ,ampiz,ampi,amro,gamro,ama1,gama1
1926 * ,amk,amkz,amkst,gamkst
1928 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1929 * ,ampiz,ampi,amro,gamro,ama1,gama1
1930 * ,amk,amkz,amkst,gamkst,fact_kpi0
1931 COMMON / taukle / bra1,brk0,brk0b,brks
1932 real*4 bra1,brk0,brk0b,brks
1937 30 qq(i)=pc(i)- pn(i)
1939 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
1940 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
1942 31 qq(i)=qq(i)-pks(i)*qqpks/pksd
1944 hadcur(k)=qq(k)*bwigs(pksd,amkst,gamkst)
1950 SUBROUTINE curr_kk0(PC,PN,HADCUR)
1959 COMPLEX BWIGS,HADCUR(4),FKK0_RCHT
1960 REAL PC(4),PN(4),QQ(4),PKS(4),PKSD,QQPKS,FPIRK
1962 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1963 * ,ampiz,ampi,amro,gamro,ama1,gama1
1964 * ,amk,amkz,amkst,gamkst
1966 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
1967 * ,ampiz,ampi,amro,gamro,ama1,gama1
1968 * ,amk,amkz,amkst,gamkst
1974 pksd =pks(4)*pks(4)-pks(3)*pks(3)-pks(2)*pks(2)-pks(1)*pks(1)
1975 qqpks=pks(4)* qq(4)-pks(3)* qq(3)-pks(2)* qq(2)-pks(1)* qq(1)
1977 31 qq(i)=qq(i)-pks(i)*qqpks/pksd
1980 hadcur(k)=qq(k)*sqrt(fpirk(sqrt(pksd)))
1996 DATA pi /3.141592653589793238462643/
2014 IF (icont.EQ.0)
THEN
2019 coefc(1,0)= 2.0*sqrt(2.)/3.0
2020 coefc(2,0)=-2.0*sqrt(2.)/3.0
2022 coefc(3,0)= 2.0*sqrt(2.)/3.0
2026 coefc(1,1)=-sqrt(2.)/3.0
2027 coefc(2,1)= sqrt(2.)/3.0
2030 coefc(5,1)= sqrt(2.)
2033 coefc(1,2)=-sqrt(2.)/3.0
2034 coefc(2,2)= sqrt(2.)/3.0
2038 coefc(5,2)=-sqrt(2.)
2048 coefc(1,4)= 1.0/sqrt(2.)/3.0
2049 coefc(2,4)=-1.0/sqrt(2.)/3.0
2054 coefc(1,5)=-sqrt(2.)/3.0
2055 coefc(2,5)= sqrt(2.)/3.0
2058 coefc(5,5)=-sqrt(2.)
2071 coefc(5,7)=-sqrt(2.0/3.0)
2074 IF (iver.EQ.0.OR.j.NE.0)
THEN
2076 ELSEIF (iver.EQ.1)
THEN
2079 write(*,*)
'wrong IVER=',iver
2085 SUBROUTINE inirchl(IVERI)
2098 CALL rchl_parameters(1)
2103 SUBROUTINE inirchlget(I)
2122 SUBROUTINE dexaa(MODE,ISGN,POL,PNU,PAA,PIM1,PIM2,PIPL,JAA)
2133 COMMON / inout / inut,iout
2134 REAL POL(4),HV(4),PAA(4),PNU(4),PIM1(4),PIM2(4),PIPL(4),RN(1)
2140 CALL dadmaa( -1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
2143 ELSEIF(mode.EQ. 0)
THEN
2146 IF(iwarm.EQ.0)
GOTO 902
2147 CALL dadmaa( 0,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
2148 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
2151 IF(rn(1).GT.wt)
GOTO 300
2153 ELSEIF(mode.EQ. 1)
THEN
2155 CALL dadmaa( 1,isgn,hv,pnu,paa,pim1,pim2,pipl,jaa)
2160 902
WRITE(iout, 9020)
2161 9020
FORMAT(
' ----- DEXAA: LACK OF INITIALISATION')
2164 SUBROUTINE dadmaa(MODE,ISGN,HHV,PNU,PAA,PIM1,PIM2,PIPL,JAA)
2168 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2169 * ,ampiz,ampi,amro,gamro,ama1,gama1
2170 * ,amk,amkz,amkst,gamkst
2172 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2173 * ,ampiz,ampi,amro,gamro,ama1,gama1
2174 * ,amk,amkz,amkst,gamkst
2175 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2176 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2177 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
2178 real*4 gampmc ,gamper
2179 COMMON / inout / inut,iout
2181 REAL HV(4),PAA(4),PNU(4),PIM1(4),PIM2(4),PIPL(4)
2182 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4),PDUM5(4)
2185 DATA pi /3.141592653589793238462643/
2198 CALL dphsaa(wt,hv,pdum1,pdum2,pdum3,pdum4,pdum5,jaa)
2199 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
2203 ELSEIF(mode.EQ. 0)
THEN
2206 IF(iwarm.EQ.0)
GOTO 902
2207 CALL dphsaa(wt,hv,pnu,paa,pim1,pim2,pipl,jaa)
2213 sswt=sswt+dble(wt)**2
2217 IF(wt.GT.wtmax) nevovr=nevovr+1
2218 IF(rn*wtmax.GT.wt)
GOTO 300
2220 costhe=-1.+2.*rrr(2)
2223 CALL rotpol(thet,phi,pnu)
2224 CALL rotpol(thet,phi,paa)
2225 CALL rotpol(thet,phi,pim1)
2226 CALL rotpol(thet,phi,pim2)
2227 CALL rotpol(thet,phi,pipl)
2228 CALL rotpol(thet,phi,hv)
2230 44 hhv(i)=-isgn*hv(i)
2233 ELSEIF(mode.EQ. 1)
THEN
2235 IF(nevraw.EQ.0)
RETURN
2236 pargam=swt/float(nevraw+1)
2238 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
2240 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
2248 7003
FORMAT(///1x,15(5h*****)
2249 $ /,
' *', 25x,
'******** DADMAA INITIALISATION ********',9x,1h*
2250 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
2251 $ /,1x,15(5h*****)/)
2252 7010
FORMAT(///1x,15(5h*****)
2253 $ /,
' *', 25x,
'******** DADMAA FINAL REPORT ******** ',9x,1h*
2254 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF A1 DECAYS TOTAL ',9x,1h*
2255 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF A1 DECS. ACCEPTED ',9x,1h*
2256 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
2257 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (A1 DECAY) IN GEV UNITS ',9x,1h*
2258 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2259 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
2260 $ /,1x,15(5h*****)/)
2261 902
WRITE(iout, 9020)
2262 9020
FORMAT(
' ----- DADMAA: LACK OF INITIALISATION')
2265 SUBROUTINE dphsaa(DGAMT,HV,PN,PAA,PIM1,PIM2,PIPL,JAA)
2270 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2271 * ,ampiz,ampi,amro,gamro,ama1,gama1
2272 * ,amk,amkz,amkst,gamkst
2274 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2275 * ,ampiz,ampi,amro,gamro,ama1,gama1
2276 * ,amk,amkz,amkst,gamkst
2277 COMMON / taukle / bra1,brk0,brk0b,brks
2278 real*4 bra1,brk0,brk0b,brks
2279 REAL HV(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
2289 IF (rmod.LT.bra1)
THEN
2302 $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
2304 SUBROUTINE dexkk(MODE,ISGN,POL,PKK,PNU)
2311 REAL POL(4),HV(4),PNU(4),PKK(4),RN(1)
2315 CALL dadmkk(-1,isgn,hv,pkk,pnu)
2318 ELSEIF(mode.EQ. 0)
THEN
2321 CALL dadmkk( 0,isgn,hv,pkk,pnu)
2322 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
2325 IF(rn(1).GT.wt)
GOTO 300
2327 ELSEIF(mode.EQ. 1)
THEN
2329 CALL dadmkk( 1,isgn,hv,pkk,pnu)
2335 SUBROUTINE dadmkk(MODE,ISGN,HV,PKK,PNU)
2338 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2339 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2340 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2341 * ,ampiz,ampi,amro,gamro,ama1,gama1
2342 * ,amk,amkz,amkst,gamkst
2344 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2345 * ,ampiz,ampi,amro,gamro,ama1,gama1
2346 * ,amk,amkz,amkst,gamkst
2347 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
2348 real*4 gampmc ,gamper
2349 COMMON / inout / inut,iout
2350 REAL PKK(4),PNU(4),HV(4)
2351 DATA pi /3.141592653589793238462643/
2356 ELSEIF(mode.EQ. 0)
THEN
2359 ekk= (amtau**2+amk**2-amnuta**2)/(2*amtau)
2360 enu= (amtau**2-amk**2+amnuta**2)/(2*amtau)
2361 xkk= sqrt(ekk**2-amk**2)
2363 CALL sphera(xkk,pkk)
2371 qxn=pkk(4)*pnu(4)-pkk(1)*pnu(1)-pkk(2)*pnu(2)-pkk(3)*pnu(3)
2372 brak=(gv**2+ga**2)*(2*pxq*qxn-amk**2*pxn)
2373 & +(gv**2-ga**2)*amtau*amnuta*amk**2
2375 40 hv(i)=-isgn*2*ga*gv*amtau*(2*pkk(i)*qxn-pnu(i)*amk**2)/brak
2378 ELSEIF(mode.EQ. 1)
THEN
2380 IF(nevtot.EQ.0)
RETURN
2387 gamm=(gfermi*fkk)**2/(16.*pi)*amtau**3*
2389 $ sqrt((amtau**2-amk**2-amnuta**2)**2
2390 $ -4*amk**2*amnuta**2 )/amtau**2
2395 WRITE(iout, 7010) nevtot,gamm,rat,error
2402 7010
FORMAT(///1x,15(5h*****)
2403 $ /,
' *', 25x,
'******** DADMKK FINAL REPORT ********',9x,1h*
2404 $ /,
' *',i20 ,5x,
'NEVTOT = NO. OF K DECAYS TOTAL ',9x,1h*,
2405 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH ( K DECAY) IN GEV UNITS ',9x,1h*,
2406 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2407 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH (STAT.)',9x,1h*
2408 $ /,1x,15(5h*****)/)
2410 SUBROUTINE dexks(MODE,ISGN,POL,PNU,PKS,PKK,PPI,JKST)
2422 COMMON / inout / inut,iout
2423 REAL POL(4),HV(4),PKS(4),PNU(4),PKK(4),PPI(4),RN(1)
2430 CALL dadmks( -1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2434 ELSEIF(mode.EQ. 0)
THEN
2437 IF(iwarm.EQ.0)
GOTO 902
2438 CALL dadmks( 0,isgn,hv,pnu,pks,pkk,ppi,jkst)
2439 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
2444 IF(rn(1).GT.wt)
GOTO 300
2446 ELSEIF(mode.EQ. 1)
THEN
2448 CALL dadmks( 1,isgn,hv,pnu,pks,pkk,ppi,jkst)
2454 902
WRITE(iout, 9020)
2455 9020
FORMAT(
' ----- DEXKS: LACK OF INITIALISATION')
2458 SUBROUTINE dadmks(MODE,ISGN,HHV,PNU,PKS,PKK,PPI,JKST)
2460 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2461 * ,ampiz,ampi,amro,gamro,ama1,gama1
2462 * ,amk,amkz,amkst,gamkst
2464 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2465 * ,ampiz,ampi,amro,gamro,ama1,gama1
2466 * ,amk,amkz,amkst,gamkst
2467 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2468 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2469 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
2470 real*4 gampmc ,gamper
2471 COMMON / taukle / bra1,brk0,brk0b,brks
2472 real*4 bra1,brk0,brk0b,brks
2473 COMMON / inout / inut,iout
2475 REAL HV(4),PKS(4),PNU(4),PKK(4),PPI(4)
2476 REAL PDUM1(4),PDUM2(4),PDUM3(4),PDUM4(4)
2477 real*4 rrr(3),rmod(1)
2479 DATA pi /3.141592653589793238462643/
2494 CALL dphsks(wt,hv,pdum1,pdum2,pdum3,pdum4,jkst)
2495 IF(wt.GT.wtmax/1.2) wtmax=wt*1.2
2500 ELSEIF(mode.EQ. 0)
THEN
2502 IF(iwarm.EQ.0)
GOTO 902
2508 IF(rmod(1).LT.dec1)
THEN
2513 CALL dphsks(wt,hv,pnu,pks,pkk,ppi,jkst)
2516 IF(wt.GT.wtmax) nevovr=nevovr+1
2520 IF(rn*wtmax.GT.wt)
GOTO 400
2522 costhe=-1.+2.*rrr(2)
2525 CALL rotor2(thet,pnu,pnu)
2526 CALL rotor3( phi,pnu,pnu)
2527 CALL rotor2(thet,pks,pks)
2528 CALL rotor3( phi,pks,pks)
2529 CALL rotor2(thet,pkk,pkk)
2530 CALL rotor3(phi,pkk,pkk)
2531 CALL rotor2(thet,ppi,ppi)
2532 CALL rotor3( phi,ppi,ppi)
2533 CALL rotor2(thet,hv,hv)
2534 CALL rotor3( phi,hv,hv)
2536 44 hhv(i)=-isgn*hv(i)
2539 ELSEIF(mode.EQ. 1)
THEN
2541 IF(nevraw.EQ.0)
RETURN
2542 pargam=swt/float(nevraw+1)
2544 IF(nevraw.NE.0) error=sqrt(sswt/swt**2-1./float(nevraw))
2546 WRITE(iout, 7010) nevraw,nevacc,nevovr,pargam,rat,error
2554 7003
FORMAT(///1x,15(5h*****)
2555 $ /,
' *', 25x,
'******** DADMKS INITIALISATION ********',9x,1h*
2556 $ /,
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*
2557 $ /,1x,15(5h*****)/)
2558 7010
FORMAT(///1x,15(5h*****)
2559 $ /,
' *', 25x,
'******** DADMKS FINAL REPORT ********',9x,1h*
2560 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF K* DECAYS TOTAL ',9x,1h*,
2561 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF K* DECS. ACCEPTED ',9x,1h*,
2562 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
2563 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH (K* DECAY) IN GEV UNITS ',9x,1h*,
2564 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
2565 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
2566 $ /,1x,15(5h*****)/)
2567 902
WRITE(iout, 9020)
2568 9020
FORMAT(
' ----- DADMKS: LACK OF INITIALISATION')
2571 SUBROUTINE dphsks(DGAMT,HV,PN,PKS,PKK,PPI,JKST)
2578 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2579 * ,ampiz,ampi,amro,gamro,ama1,gama1
2580 * ,amk,amkz,amkst,gamkst
2582 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2583 * ,ampiz,ampi,amro,gamro,ama1,gama1
2584 * ,amk,amkz,amkst,gamkst
2585 REAL HV(4),PT(4),PN(4),PKS(4),PKK(4),PPI(4),QQ(4),RR1(1),RR2(1)
2586 DATA pi /3.141592653589793238462643/
2590 phspac=1./2**11/pi**5
2602 ams2=(amtau-amnuta)**2
2603 alp1=atan((ams1-amkst**2)/amkst/gamkst)
2604 alp2=atan((ams2-amkst**2)/amkst/gamkst)
2607 IF (rr2(1).LT.prob1)
THEN
2609 amx2=ams1+ rr1(1)*(ams2-ams1)
2613 alp=alp1+rr1(1)*(alp2-alp1)
2614 amx2=amkst**2+amkst*gamkst*tan(alp)
2619 phspac2=((amx2-amkst**2)**2+(amkst*gamkst)**2)
2621 phspac2=phspac2*(alp2-alp1)
2624 IF (phspac1.NE.0.0) a1=prob1 /phspac1
2625 IF (phspac2.NE.0.0) a2=(1-prob1)/phspac2
2629 IF (a1+a2.NE.0.0)
THEN
2630 phspac=phspac/(a1+a2)
2639 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2640 pn(3)=-sqrt(abs((pn(4)-amnuta)*(pn(4)+amnuta)))
2645 pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2647 phspac=phspac*(4*pi)*(2*pks(3)/amtau)
2650 enpi=( amx**2+ampi**2-amkz**2 ) / ( 2*amx )
2651 pppi=sqrt(abs(enpi-ampi)*(enpi+ampi))
2652 phspac=phspac*(4*pi)*(2*pppi/amx)
2654 CALL sphera(pppi,ppi)
2659 pkk(4)=( amx**2+amkz**2-ampi**2 ) / ( 2*amx )
2660 exe=(pks(4)+pks(3))/amx
2662 CALL bostr3(exe,ppi,ppi)
2663 CALL bostr3(exe,pkk,pkk)
2665 CALL dam2pi(1,pt,pn,ppi,pkk,amplit,hv)
2666 dgamt=1/(2.*amtau)*amplit*phspac
2669 ELSEIF(jkst.EQ.20)
THEN
2673 ams2=(amtau-amnuta)**2
2674 alp1=atan((ams1-amkst**2)/amkst/gamkst)
2675 alp2=atan((ams2-amkst**2)/amkst/gamkst)
2679 IF (rr2(1).LT.prob1)
THEN
2681 amx2=ams1+ rr1(1)*(ams2-ams1)
2685 alp=alp1+rr1(1)*(alp2-alp1)
2686 amx2=amkst**2+amkst*gamkst*tan(alp)
2691 phspac2=((amx2-amkst**2)**2+(amkst*gamkst)**2)
2693 phspac2=phspac2*(alp2-alp1)
2696 IF (phspac1.NE.0.0) a1=prob1 /phspac1
2697 IF (phspac2.NE.0.0) a2=(1-prob1)/phspac2
2699 IF (a1+a2.NE.0.0)
THEN
2700 phspac=phspac/(a1+a2)
2709 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
2710 pn(3)=-sqrt(abs((pn(4)-amnuta)*(pn(4)+amnuta)))
2714 pks(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
2716 phspac=phspac*(4*pi)*(2*pks(3)/amtau)
2719 enpi=( amx**2+ampiz**2-amk**2 ) / ( 2*amx )
2720 pppi=sqrt(abs(enpi-ampiz)*(enpi+ampiz))
2721 phspac=phspac*(4*pi)*(2*pppi/amx)
2723 CALL sphera(pppi,ppi)
2728 pkk(4)=( amx**2+amk**2-ampiz**2 ) / ( 2*amx )
2729 exe=(pks(4)+pks(3))/amx
2731 CALL bostr3(exe,ppi,ppi)
2732 CALL bostr3(exe,pkk,pkk)
2734 CALL dam2pi(2,pt,pn,pkk,ppi,amplit,hv)
2735 dgamt=1/(2.*amtau)*amplit*phspac
2744 SUBROUTINE dphnpi(DGAMT,HVX,PNX,PRX,PPIX,JNPI)
2749 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2750 * ,ampiz,ampi,amro,gamro,ama1,gama1
2751 * ,amk,amkz,amkst,gamkst
2753 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2754 * ,ampiz,ampi,amro,gamro,ama1,gama1
2755 * ,amk,amkz,amkst,gamkst
2756 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
2757 real*4 gfermi,gv,ga,ccabib,scabib,gamel
2758 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
2759 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
2761 CHARACTER NAMES(NMODE)*31
2764 real*8 pn(4),pr(4),ppi(4,9),hv(4)
2765 real*4 pnx(4),prx(4),ppix(4,9),hvx(4)
2766 real*8 pv(5,9),pt(4),ue(3),be(3)
2767 real*8 pawt,amx,ams1,ams2,pa,phs,phsmax,pmin,pmax
2771 real*8 gam,bep,phi,a,b,c
2773 real*4 rrr(9),rrx(2),rn(1),rr2(1)
2775 DATA pi /3.141592653589793238462643/
2776 DATA wetmax /20*1d-15/
2781 $ sqrt(max(0.d0,(a**2-(b+c)**2)*(a**2-(b-c)**2)))/(2.d0*a)
2783 ampik(i,j)=dcdmas(idffin(i,j))
2786 IF ((jnpi.LE.0).OR.jnpi.GT.20)
THEN
2787 WRITE(6,*)
'JNPI OUTSIDE RANGE DEFINED BY WETMAX; JNPI=',jnpi
2801 phspac = 1./2.**5 /pi**2
2803 4 ps =ps+ampik(i,jnpi)
2806 ams2=(amtau-amnuta)**2
2809 amx2=ams1+ rr2(1)*(ams2-ams1)
2812 phspac=phspac * (ams2-ams1)
2817 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx2)
2818 pn(3)=-sqrt(abs((pn(4)-amnuta)*(pn(4)+amnuta)))
2822 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx2)
2824 phspac=phspac * (4.*pi) * (2.*pr(3)/amtau)
2831 qxn=pr(4)*pn(4)-pr(1)*pn(1)-pr(2)*pn(2)-pr(3)*pn(3)
2834 brak=2*(gv**2+ga**2)*(2*pxq*qxn+amx2*pxn)
2835 & -6*(gv**2-ga**2)*amtau*amnuta*amx2
2838 amplit=ccabib**2*gfermi**2/2. * brak * amx2*sigee(amx2,jnpi)
2839 dgamt=1./(2.*amtau)*amplit*phspac
2846 pv(5,nd)=ampik(nd,jnpi)
2848 pmax=amw-ps+ampik(nd,jnpi)
2851 pmax=pmax+ampik(il,jnpi)
2852 pmin=pmin+ampik(il+1,jnpi)
2853 220 phsmax=phsmax*pawt(pmax,pmin,ampik(il,jnpi))/pmax
2860 223 ams1=ams1+ampik(jl,jnpi)
2862 amx =(amx-ampik(il,jnpi))
2864 phsmax=phsmax * (ams2-ams1)
2871 phspac = 1./2.**(6*nd-7) /pi**(3*nd-4)
2873 CALL ranmar(rrr,nd-2)
2877 231 ams1=ams1+ampik(jl,jnpi)
2879 ams2=(amx-ampik(il,jnpi))**2
2881 amx2=ams1+ rr1*(ams2-ams1)
2884 phspac=phspac * (ams2-ams1)
2886 phs=phs* (ams2-ams1)
2887 pa=pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2888 phs =phs *pa/pv(5,il)
2890 pa=pawt(pv(5,nd-1),ampik(nd-1,jnpi),ampik(nd,jnpi))
2891 phs =phs *pa/pv(5,nd-1)
2894 IF(phsmax.NE.0.0)
THEN
2896 wetmax(jnpi)=1.2d0*max(wetmax(jnpi)/1.2d0,phs/phsmax)
2898 wetmax(jnpi)=1.2d0*wetmax(jnpi)/1.2d0
2901 IF (ncont.EQ.500 000)
THEN
2904 xnpi=xnpi+ampik(kk,jnpi)
2906 WRITE(6,*)
'ROUNDING INSTABILITY IN DPHNPI ?'
2907 WRITE(6,*)
'AMW=',amw,
'XNPI=',xnpi
2908 WRITE(6,*)
'IF =AMW= IS NEARLY EQUAL =XNPI= THAT IS IT'
2909 WRITE(6,*)
'PHS=',phs,
'PHSMAX=',phsmax
2912 IF(rn(1)*phsmax*wetmax(jnpi).GT.phs)
GO TO 100
2914 280
DO 300 il=1,nd-1
2915 pa=pawt(pv(5,il),pv(5,il+1),ampik(il,jnpi))
2919 ue(1)=sqrt(1.d0-ue(3)**2)*cos(phi)
2920 ue(2)=sqrt(1.d0-ue(3)**2)*sin(phi)
2923 290 pv(j,il+1)=-pa*ue(j)
2924 ppi(4,il)=sqrt(pa**2+ampik(il,jnpi)**2)
2925 pv(4,il+1)=sqrt(pa**2+pv(5,il+1)**2)
2926 phspac=phspac *(4.*pi)*(2.*pa/pv(5,il))
2930 310 ppi(j,nd)=pv(j,nd)
2933 320 be(j)=pv(j,il)/pv(4,il)
2934 gam=pv(4,il)/pv(5,il)
2936 bep=be(1)*ppi(1,i)+be(2)*ppi(2,i)+be(3)*ppi(3,i)
2938 330 ppi(j,i)=ppi(j,i)+gam*(gam*bep/(1.d0+gam)+ppi(4,i))*be(j)
2939 ppi(4,i)=gam*(ppi(4,i)+bep)
2956 FUNCTION sigee(Q2,JNP)
2975 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
2976 * ,ampiz,ampi,amro,gamro,ama1,gama1
2977 * ,amk,amkz,amkst,gamkst
2979 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
2980 * ,ampiz,ampi,amro,gamro,ama1,gama1
2981 * ,amk,amkz,amkst,gamkst
2985 1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
2986 2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
2987 3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
2988 4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
2991 7 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25,
2994 DATA pi /3.141592653589793238462643/
3007 datsig(i,2) = datsig(i,2)/2.
3008 datsig(i,1) = datsig(i,1) + datsig(i,2)
3014 IF(t . gt. s-ampi )
GO TO 201
3016 fact=(t2/s2)**2*sqrt((s2-t2-ampi2)**2-4.*t2*ampi2)/s2 *2.*t*.05
3017 fact = fact * (datsig(j,1)+datsig(j+1,1))
3018 200 datsig(i,3) = datsig(i,3) + fact
3019 201 datsig(i,3) = datsig(i,3) /(2*pi*fpi)**2
3020 datsig(i,4) = datsig(i,3)
3021 datsig(i,6) = datsig(i,5)
3024 1000
FORMAT(///1x,
' EE SIGMA USED IN MULTIPI DECAYS'/
3030 sigee=datsig(1,jnpi)+
3031 & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
3032 ELSEIF(q.LT.1.8)
THEN
3035 IF(q.LT.qmax)
GO TO 2
3038 2 sigee=datsig(i,jnpi)+
3039 & (datsig(i+1,jnpi)-datsig(i,jnpi)) * (q-qmin)/.05
3040 ELSEIF(q.GT.1.8)
THEN
3041 sigee=datsig(17,jnpi)+
3042 & (datsig(17,jnpi)-datsig(16,jnpi)) * (q-1.8)/.05
3044 IF(sigee.LT..0) sigee=0.
3046 sigee = sigee/(6.*pi**2*sig0)
3051 FUNCTION sigold(Q2,JNPI)
3070 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3071 * ,ampiz,ampi,amro,gamro,ama1,gama1
3072 * ,amk,amkz,amkst,gamkst
3074 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3075 * ,ampiz,ampi,amro,gamro,ama1,gama1
3076 * ,amk,amkz,amkst,gamkst
3080 1 7.40,12.00,16.15,21.25,24.90,29.55,34.15,37.40,37.85,37.40,
3081 2 36.00,33.25,30.50,27.70,24.50,21.25,18.90,
3082 3 1.24, 2.50, 3.70, 5.40, 7.45,10.75,14.50,18.20,22.30,28.90,
3083 4 29.35,25.60,22.30,18.60,14.05,11.60, 9.10,
3085 6 9*.0,.65,1.25,2.20,3.15,5.00,5.75,7.80,8.25/
3087 DATA pi /3.141592653589793238462643/
3095 datsig(i,2) = datsig(i,2)/2.
3096 datsig(i,1) = datsig(i,1) + datsig(i,2)
3102 IF(t . gt. s-ampi )
GO TO 201
3104 fact=(t2/s2)**2*sqrt((s2-t2-ampi2)**2-4.*t2*ampi2)/s2 *2.*t*.05
3105 fact = fact * (datsig(j,1)+datsig(j+1,1))
3106 200 datsig(i,3) = datsig(i,3) + fact
3107 201 datsig(i,3) = datsig(i,3) /(2*pi*fpi)**2
3110 1000
FORMAT(///1x,
' EE SIGMA USED IN MULTIPI DECAYS'/
3116 sigee=datsig(1,jnpi)+
3117 & (datsig(2,jnpi)-datsig(1,jnpi))*(q-1.)/.05
3118 ELSEIF(q.LT.1.8)
THEN
3121 IF(q.LT.qmax)
GO TO 2
3124 2 sigee=datsig(i,jnpi)+
3125 & (datsig(i+1,jnpi)-datsig(i,jnpi)) * (q-qmin)/.05
3126 ELSEIF(q.GT.1.8)
THEN
3127 sigee=datsig(17,jnpi)+
3128 & (datsig(17,jnpi)-datsig(16,jnpi)) * (q-1.8)/.05
3130 IF(sigee.LT..0) sigee=0.
3132 sigee = sigee/(6.*pi**2*sig0)
3137 SUBROUTINE dphspk(DGAMT,HV,PN,PAA,PNPI,JAA)
3142 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
3143 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3145 CHARACTER NAMES(NMODE)*31
3147 REAL HV(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4),PNPI(4,9)
3154 amp1=dcdmas(idffin(1,jaa+nm4+nm5+nm6))
3155 amp2=dcdmas(idffin(2,jaa+nm4+nm5+nm6))
3156 amp3=dcdmas(idffin(3,jaa+nm4+nm5+nm6))
3158 $ dphtre(dgamt,hv,pn,paa,pim1,amp1,pim2,amp2,pipl,amp3,keyt,mnum)
3170 $ dphtre(dgamt,hv,pn,paa,pim1,ampa,pim2,ampb,pipl,amp3,keyt,mnum)
3187 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3188 * ,ampiz,ampi,amro,gamro,ama1,gama1
3189 * ,amk,amkz,amkst,gamkst
3191 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3192 * ,ampiz,ampi,amro,gamro,ama1,gama1
3193 * ,amk,amkz,amkst,gamkst
3194 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3195 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3196 REAL HV(4),PT(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
3199 DATA pi /3.141592653589793238462643/
3201 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
3206 phspac=1./2**17/pi**8
3216 CALL choice(mnum,rr,ichan,prob1,prob2,prob3,
3217 $ amrx,gamrx,amra,gamra,amrb,gamrb)
3218 IF (ichan.EQ.1)
THEN
3221 ELSEIF (ichan.EQ.2)
THEN
3230 ams1=(amp1+amp2+amp3)**2
3231 ams2=(amtau-amnuta)**2
3233 alp1=atan((ams1-amrx**2)/amrx/gamrx)
3234 alp2=atan((ams2-amrx**2)/amrx/gamrx)
3235 alp=alp1+rr1*(alp2-alp1)
3236 am3sq =amrx**2+amrx*gamrx*tan(alp)
3238 phspac=phspac*((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
3239 phspac=phspac*(alp2-alp1)
3244 IF (ichan.LE.2)
THEN
3246 alp1=atan((ams1-amra**2)/amra/gamra)
3247 alp2=atan((ams2-amra**2)/amra/gamra)
3248 alp=alp1+rr2*(alp2-alp1)
3249 am2sq =amra**2+amra*gamra*tan(alp)
3257 am2sq=ams1+ rr2*(ams2-ams1)
3262 enq1=(am2sq-amp2**2+amp3**2)/(2*am2)
3263 enq2=(am2sq+amp2**2-amp3**2)/(2*am2)
3264 ppi= enq1**2-amp3**2
3265 pppi=sqrt(abs(enq1**2-amp3**2))
3267 phf1=(4*pi)*(2*pppi/am2)
3269 CALL sphera(pppi,pipl)
3279 pr(4)=1./(2*am3)*(am3**2+am2**2-amp1**2)
3280 pr(3)= sqrt(abs(pr(4)**2-am2**2))
3281 ppi = pr(4)**2-am2**2
3285 pim2(4)=1./(2*am3)*(am3**2-am2**2+amp1**2)
3287 phf2=(4*pi)*(2*pr(3)/am3)
3289 exe=(pr(4)+pr(3))/am2
3290 CALL bostr3(exe,pipl,pipl)
3291 CALL bostr3(exe,pim1,pim1)
3295 thet =acos(-1.+2*rr3)
3297 CALL rotpol(thet,phi,pipl)
3298 CALL rotpol(thet,phi,pim1)
3299 CALL rotpol(thet,phi,pim2)
3300 CALL rotpol(thet,phi,pr)
3306 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am3**2)
3307 paa(3)= sqrt(abs(paa(4)**2-am3**2))
3308 ppi = paa(4)**2-am3**2
3309 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
3313 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am3**2)
3319 alp1=atan((ams1-amra**2)/amra/gamra)
3320 alp2=atan((ams2-amra**2)/amra/gamra)
3321 xpro = (pim1(3)+pipl(3))**2
3322 $ +(pim1(2)+pipl(2))**2+(pim1(1)+pipl(1))**2
3323 am2sq=-xpro+(pim1(4)+pipl(4))**2
3325 ff1 = ((am2sq-amra**2)**2+(amra*gamra)**2)/(amra*gamra)
3326 ff1 =ff1 *(alp2-alp1)
3328 gg1 = (4*pi)*(xlam(am2sq,amp2**2,amp3**2)/am2sq)
3330 gg1 =gg1 *(4*pi)*sqrt(4*xpro/am3sq)
3331 xjaje=gg1*(ams2-ams1)
3335 alp1=atan((ams1-amrb**2)/amrb/gamrb)
3336 alp2=atan((ams2-amrb**2)/amrb/gamrb)
3337 xpro = (pim2(3)+pipl(3))**2
3338 $ +(pim2(2)+pipl(2))**2+(pim2(1)+pipl(1))**2
3339 am2sq=-xpro+(pim2(4)+pipl(4))**2
3340 ff2 = ((am2sq-amrb**2)**2+(amrb*gamrb)**2)/(amrb*gamrb)
3341 ff2 =ff2 *(alp2-alp1)
3342 gg2 = (4*pi)*(xlam(am2sq,amp1**2,amp3**2)/am2sq)
3343 gg2 =gg2 *(4*pi)*sqrt(4*xpro/am3sq)
3344 xjadw=gg2*(ams2-ams1)
3351 IF (ichan.EQ.2)
THEN
3356 IF (xjac1.NE.0.0) a1=prob1/xjac1
3357 IF (xjac2.NE.0.0) a2=prob2/xjac2
3358 IF (xjac3.NE.0.0) a3=prob3/xjac3
3360 IF (a1+a2+a3.NE.0.0)
THEN
3361 phspac=phspac/(a1+a2+a3)
3373 exe=(paa(4)+paa(3))/am3
3374 CALL bostr3(exe,pipl,pipl)
3375 CALL bostr3(exe,pim1,pim1)
3376 CALL bostr3(exe,pim2,pim2)
3377 CALL bostr3(exe,pr,pr)
3380 CALL dampog(pt,pn,pim1,pim2,pipl,amplit,hv)
3384 CALL damppk(mnum,pt,pn,pim1,pim2,pipl,amplit,hv)
3386 IF (keyt.EQ.1.OR.keyt.EQ.2)
THEN
3392 dgamt=1/(2.*amtau)*amplit*phspac
3394 SUBROUTINE dampaa(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3404 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3405 * ,ampiz,ampi,amro,gamro,ama1,gama1
3406 * ,amk,amkz,amkst,gamkst
3408 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3409 * ,ampiz,ampi,amro,gamro,ama1,gama1
3410 * ,amk,amkz,amkst,gamkst
3411 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3412 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3413 COMMON /testa1/ keya1
3414 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4)
3415 REAL PAA(4),VEC1(4),VEC2(4)
3416 REAL PIVEC(4),PIAKS(4),HVM(4)
3417 COMPLEX BWIGN,HADCUR(4),FPIK
3424 bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
3428 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3430 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3431 xmro1 =sqrt(abs((pipl(4)+pim1(4))**2-(pipl(1)+pim1(1))**2
3432 $ -(pipl(2)+pim1(2))**2-(pipl(3)+pim1(3))**2))
3433 xmro2 =sqrt(abs((pipl(4)+pim2(4))**2-(pipl(1)+pim2(1))**2
3434 $ -(pipl(2)+pim2(2))**2-(pipl(3)+pim2(3))**2))
3436 prod1 =paa(4)*(pim1(4)-pipl(4))-paa(1)*(pim1(1)-pipl(1))
3437 $ -paa(2)*(pim1(2)-pipl(2))-paa(3)*(pim1(3)-pipl(3))
3438 prod2 =paa(4)*(pim2(4)-pipl(4))-paa(1)*(pim2(1)-pipl(1))
3439 $ -paa(2)*(pim2(2)-pipl(2))-paa(3)*(pim2(3)-pipl(3))
3441 vec1(i)= pim1(i)-pipl(i) -paa(i)*prod1/xmaa**2
3442 40 vec2(i)= pim2(i)-pipl(i) -paa(i)*prod2/xmaa**2
3444 IF (keya1.EQ.1)
THEN
3448 fnorm=fa1/sqrt(2.)*faropi*fro2pi
3450 hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gama1)
3451 $ *(cmplx(vec1(i))*amro**2*bwign(xmro1,amro,gamro)
3452 $ +cmplx(vec2(i))*amro**2*bwign(xmro2,amro,gamro))
3455 fnorm=2.0*sqrt(2.)/3.0/fpi
3456 gamax=gama1*gfun(xmaa**2)/gfun(ama1**2)
3458 hadcur(i)= cmplx(fnorm) *ama1**2*bwign(xmaa,ama1,gamax)
3459 $ *(cmplx(vec1(i))*fpik(xmro1)
3460 $ +cmplx(vec2(i))*fpik(xmro2))
3465 CALL clvec(hadcur,pn,pivec)
3466 CALL claxi(hadcur,pn,piaks)
3467 CALL clnut(hadcur,brakm,hvm)
3469 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3470 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3471 amplit=(gfermi*ccabib)**2*brak/2.
3476 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3477 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3487 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3488 * ,ampiz,ampi,amro,gamro,ama1,gama1
3489 * ,amk,amkz,amkst,gamkst
3491 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3492 * ,ampiz,ampi,amro,gamro,ama1,gama1
3493 * ,amk,amkz,amkst,gamkst
3495 IF (qkwa.LT.(amro+ampi)**2)
THEN
3496 gfun=4.1*(qkwa-9*ampiz**2)**3
3497 $ *(1.-3.3*(qkwa-9*ampiz**2)+5.8*(qkwa-9*ampiz**2)**2)
3499 gfun=qkwa*(1.623+10.38/qkwa-9.32/qkwa**2+0.65/qkwa**3)
3502 COMPLEX FUNCTION bwigs(S,M,G)
3507 REAL PI,PIM,QS,QM,W,GS,MK
3512 p(a,b,c)=sqrt(abs(abs(((a+b-c)**2-4.*a*b)/4./a)
3513 $ +(((a+b-c)**2-4.*a*b)/4./a))/2.0)
3526 qs=p(s,pim**2,mk**2)
3527 qm=p(m**2,pim**2,mk**2)
3529 gs=g*(m/w)*(qs/qm)**3
3530 bw=m**2/cmplx(m**2-s,-m*gs)
3531 qpm=p(pm**2,pim**2,mk**2)
3532 g1=pg*(pm/w)*(qs/qpm)**3
3533 bwp=pm**2/cmplx(pm**2-s,-pm*g1)
3534 bwigs= (bw+pbeta*bwp)/(1+pbeta)
3537 COMPLEX FUNCTION bwig(S,M,G)
3542 REAL PI,PIM,QS,QM,W,GS
3551 IF (s.GT.4.*pim**2)
THEN
3552 qs=sqrt(abs(abs(s/4.-pim**2)+(s/4.-pim**2))/2.0)
3553 qm=sqrt(m**2/4.-pim**2)
3555 gs=g*(m/w)*(qs/qm)**3
3559 bwig=m**2/cmplx(m**2-s,-m*gs)
3562 COMPLEX FUNCTION fpik(W)
3567 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
3572 IF (init.EQ.0 )
THEN
3584 fpik= (bwig(s,rom,rog)+beta1*bwig(s,rom1,rog1))
3593 fpirho=cabs(fpik(w))**2
3595 SUBROUTINE clvec(HJ,PN,PIV)
3605 hn= hj(4)*cmplx(pn(4))-hj(3)*cmplx(pn(3))
3606 hh= real(hj(4)*conjg(hj(4))-hj(3)*conjg(hj(3))
3607 $ -hj(2)*conjg(hj(2))-hj(1)*conjg(hj(1)))
3609 10 piv(i)=4.*real(hn*conjg(hj(i)))-2.*hh*pn(i)
3612 SUBROUTINE claxi(HJ,PN,PIA)
3619 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3620 COMMON / idfc / idff
3622 COMPLEX HJ(4),HJC(4)
3625 det2(i,j)=aimag(hjc(i)*hj(j)-hjc(j)*hj(i))
3630 IF (ktom.EQ.1.OR.ktom.EQ.-1)
THEN
3631 sign= idff/abs(idff)
3632 ELSEIF (ktom.EQ.2)
THEN
3633 sign=-idff/abs(idff)
3635 print *,
'STOP IN CLAXI: KTOM=',ktom
3640 10 hjc(i)=conjg(hj(i))
3641 pia(1)= -2.*pn(3)*det2(2,4)+2.*pn(4)*det2(2,3)
3642 pia(2)= -2.*pn(4)*det2(1,3)+2.*pn(3)*det2(1,4)
3643 pia(3)= 2.*pn(4)*det2(1,2)
3644 pia(4)= 2.*pn(3)*det2(1,2)
3647 20 pia(i)=pia(i)*sign
3649 SUBROUTINE clnut(HJ,B,HV)
3661 b=real( hj(4)*aimag(hj(4)) - hj(3)*aimag(hj(3))
3662 & - hj(2)*aimag(hj(2)) - hj(1)*aimag(hj(1)) )
3665 SUBROUTINE dampog(PT,PN,PIM1,PIM2,PIPL,AMPLIT,HV)
3675 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3676 * ,ampiz,ampi,amro,gamro,ama1,gama1
3677 * ,amk,amkz,amkst,gamkst
3679 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3680 * ,ampiz,ampi,amro,gamro,ama1,gama1
3681 * ,amk,amkz,amkst,gamkst
3682 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3683 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3684 COMMON /testa1/ keya1
3685 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIPL(4)
3686 REAL PAA(4),VEC1(4),VEC2(4)
3687 REAL PIVEC(4),PIAKS(4),HVM(4)
3688 COMPLEX BWIGN,HADCUR(4),FNORM,FORMOM
3698 10 paa(i)=pim1(i)+pim2(i)+pipl(i)
3701 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3702 xmom =sqrt(abs( (pim2(4)+pipl(4))**2-(pim2(3)+pipl(3))**2
3703 $ -(pim2(2)+pipl(2))**2-(pim2(1)+pipl(1))**2 ))
3704 xmro2 =(pipl(1))**2 +(pipl(2))**2 +(pipl(3))**2
3706 prod1 =vec1(1)*pipl(1)
3707 prod2 =vec2(2)*pipl(2)
3708 p12 =pim1(4)*pim2(4)-pim1(1)*pim2(1)
3709 $ -pim1(2)*pim2(2)-pim1(3)*pim2(3)
3710 p1pl =pim1(4)*pipl(4)-pim1(1)*pipl(1)
3711 $ -pim1(2)*pipl(2)-pim1(3)*pipl(3)
3712 p2pl =pipl(4)*pim2(4)-pipl(1)*pim2(1)
3713 $ -pipl(2)*pim2(2)-pipl(3)*pim2(3)
3715 vec1(i)= (vec1(i)-prod1/xmro2*pipl(i))
3717 gnorm=sqrt(vec1(1)**2+vec1(2)**2+vec1(3)**2)
3719 vec1(i)= vec1(i)/gnorm
3721 vec2(1)=(vec1(2)*pipl(3)-vec1(3)*pipl(2))/sqrt(xmro2)
3722 vec2(2)=(vec1(3)*pipl(1)-vec1(1)*pipl(3))/sqrt(xmro2)
3723 vec2(3)=(vec1(1)*pipl(2)-vec1(2)*pipl(1))/sqrt(xmro2)
3724 p1vec1 =pim1(4)*vec1(4)-pim1(1)*vec1(1)
3725 $ -pim1(2)*vec1(2)-pim1(3)*vec1(3)
3726 p2vec1 =vec1(4)*pim2(4)-vec1(1)*pim2(1)
3727 $ -vec1(2)*pim2(2)-vec1(3)*pim2(3)
3728 p1vec2 =pim1(4)*vec2(4)-pim1(1)*vec2(1)
3729 $ -pim1(2)*vec2(2)-pim1(3)*vec2(3)
3730 p2vec2 =vec2(4)*pim2(4)-vec2(1)*pim2(1)
3731 $ -vec2(2)*pim2(2)-vec2(3)*pim2(3)
3733 fnorm=formom(xmaa,xmom)
3738 hadcur(i) = fnorm *(
3739 $ vec1(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3740 $ -pim2(i)*(p2vec1*p1pl-p1vec1*p2pl)
3741 $ +pipl(i)*(p2vec1*p12 -p1vec1*(ampi**2+p2pl)) )
3743 hadcur(i) = fnorm *(
3744 $ vec2(i)*(ampi**2*p1pl-p2pl*(p12-p1pl))
3745 $ -pim2(i)*(p2vec2*p1pl-p1vec2*p2pl)
3746 $ +pipl(i)*(p2vec2*p12 -p1vec2*(ampi**2+p2pl)) )
3751 CALL clvec(hadcur,pn,pivec)
3752 CALL claxi(hadcur,pn,piaks)
3753 CALL clnut(hadcur,brakm,hvm)
3755 brak=brak+(gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3756 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3758 hv(i)=hv(i)-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3759 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3763 amplit=(gfermi*ccabib)**2*brak/2.
3772 SUBROUTINE damppk(MNUM,PT,PN,PIM1,PIM2,PIM3,AMPLIT,HV)
3782 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3783 * ,ampiz,ampi,amro,gamro,ama1,gama1
3784 * ,amk,amkz,amkst,gamkst
3786 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3787 * ,ampiz,ampi,amro,gamro,ama1,gama1
3788 * ,amk,amkz,amkst,gamkst
3789 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3790 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3793 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4)
3794 REAL PAA(4),VEC1(4),VEC2(4),VEC3(4),VEC4(4),VEC5(4)
3795 REAL PIVEC(4),PIAKS(4),HVM(4)
3796 REAL FNORM(0:7),COEF
3797 DOUBLE PRECISION GETFPIRPT
3798 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5,UROJ
3799 COMPLEX F1,F2,F3,F4,F5
3800 EXTERNAL form1,form2,form3,form4,form5
3801 DATA pi /3.141592653589793238462643/
3805 IF (icont.EQ.0)
THEN
3811 IF (iver.EQ.0.OR.mnum.NE.0)
THEN
3813 ELSEIF (iver.EQ.1)
THEN
3816 write(*,*)
'wrong IVER=',iver
3823 fnorm(4)=scabib/fpi/dwapi0
3831 10 paa(i)=pim1(i)+pim2(i)+pim3(i)
3832 xmaa =sqrt(abs(paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2))
3833 xmro1 =sqrt(abs((pim3(4)+pim2(4))**2-(pim3(1)+pim2(1))**2
3834 $ -(pim3(2)+pim2(2))**2-(pim3(3)+pim2(3))**2))
3835 xmro2 =sqrt(abs((pim3(4)+pim1(4))**2-(pim3(1)+pim1(1))**2
3836 $ -(pim3(2)+pim1(2))**2-(pim3(3)+pim1(3))**2))
3837 xmro3 =sqrt(abs((pim1(4)+pim2(4))**2-(pim1(1)+pim2(1))**2
3838 $ -(pim1(2)+pim2(2))**2-(pim1(3)+pim2(3))**2))
3840 prod1 =paa(4)*(pim2(4)-pim3(4))-paa(1)*(pim2(1)-pim3(1))
3841 $ -paa(2)*(pim2(2)-pim3(2))-paa(3)*(pim2(3)-pim3(3))
3842 prod2 =paa(4)*(pim3(4)-pim1(4))-paa(1)*(pim3(1)-pim1(1))
3843 $ -paa(2)*(pim3(2)-pim1(2))-paa(3)*(pim3(3)-pim1(3))
3844 prod3 =paa(4)*(pim1(4)-pim2(4))-paa(1)*(pim1(1)-pim2(1))
3845 $ -paa(2)*(pim1(2)-pim2(2))-paa(3)*(pim1(3)-pim2(3))
3847 vec1(i)= pim2(i)-pim3(i) -paa(i)*prod1/xmaa**2
3848 vec2(i)= pim3(i)-pim1(i) -paa(i)*prod2/xmaa**2
3849 vec3(i)= pim1(i)-pim2(i) -paa(i)*prod3/xmaa**2
3850 40 vec4(i)= pim1(i)+pim2(i)+pim3(i)
3851 CALL prod5(pim1,pim2,pim3,vec5)
3855 f1 = cmplx(coef(1,mnum))*form1(mnum,xmaa**2,xmro1**2,xmro2**2)
3856 f2 = cmplx(coef(2,mnum))*form2(mnum,xmaa**2,xmro2**2,xmro1**2)
3857 f3 = cmplx(coef(3,mnum))*form3(mnum,xmaa**2,xmro3**2,xmro1**2)
3859 $cmplx(coef(4,mnum))*form4(mnum,xmaa**2,xmro1**2,xmro2**2,xmro3**2)
3860 f5 = (-1.0)*uroj/4.0/pi**2/fpi**2*
3861 $ cmplx(coef(5,mnum))*form5(mnum,xmaa**2,xmro1**2,xmro2**2)
3864 hadcur(i)= cmplx(fnorm(mnum)) * (
3865 $ cmplx(vec1(i))*f1+cmplx(vec2(i))*f2+cmplx(vec3(i))*f3+
3866 $ cmplx(vec4(i))*f4+cmplx(vec5(i))*f5)
3870 CALL clvec(hadcur,pn,pivec)
3871 CALL claxi(hadcur,pn,piaks)
3872 CALL clnut(hadcur,brakm,hvm)
3874 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
3875 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
3876 amplit=(gfermi)**2*brak/2.
3878 print *,
'MNUM=',mnum
3884 IF (k.EQ.4) znak=1.0
3885 xm1=znak*pim1(k)**2+xm1
3886 xm2=znak*pim2(k)**2+xm2
3887 xm3=znak*pim3(k)**2+xm3
3888 77 print *,
'PIM1=',pim1(k),
'PIM2=',pim2(k),
'PIM3=',pim3(k)
3889 print *,
'XM1=',sqrt(xm1),
'XM2=',sqrt(xm2),
'XM3=',sqrt(xm3)
3890 print *,
'************************************************'
3894 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
3895 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
3900 SUBROUTINE prod5(P1,P2,P3,PIA)
3906 COMMON / jaki / jak1,jak2,jakp,jakm,ktom
3907 COMMON / idfc / idff
3908 REAL PIA(4),P1(4),P2(4),P3(4)
3909 det2(i,j)=p1(i)*p2(j)-p2(i)*p1(j)
3911 IF (ktom.EQ.1.OR.ktom.EQ.-1)
THEN
3912 sign= idff/abs(idff)
3913 ELSEIF (ktom.EQ.2)
THEN
3914 sign=-idff/abs(idff)
3916 print *,
'STOP IN PROD5: KTOM=',ktom
3922 pia(1)= -p3(3)*det2(2,4)+p3(4)*det2(2,3)+p3(2)*det2(3,4)
3923 pia(2)= -p3(4)*det2(1,3)+p3(3)*det2(1,4)-p3(1)*det2(3,4)
3924 pia(3)= p3(4)*det2(1,2)-p3(2)*det2(1,4)+p3(1)*det2(2,4)
3925 pia(4)= p3(3)*det2(1,2)-p3(2)*det2(1,3)+p3(1)*det2(2,3)
3928 20 pia(i)=pia(i)*sign
3931 SUBROUTINE dexnew(MODE,ISGN,POL,PNU,PAA,PNPI,JNPI)
3942 COMMON / inout / inut,iout
3943 REAL POL(4),HV(4),PAA(4),PNU(4),PNPI(4,9),RN(1)
3949 CALL dadnew( -1,isgn,hv,pnu,paa,pnpi,jdumm)
3952 ELSEIF(mode.EQ. 0)
THEN
3955 IF(iwarm.EQ.0)
GOTO 902
3956 CALL dadnew( 0,isgn,hv,pnu,paa,pnpi,jnpi)
3957 wt=(1+pol(1)*hv(1)+pol(2)*hv(2)+pol(3)*hv(3))/2.
3960 IF(rn(1).GT.wt)
GOTO 300
3962 ELSEIF(mode.EQ. 1)
THEN
3964 CALL dadnew( 1,isgn,hv,pnu,paa,pnpi,jdumm)
3969 902
WRITE(iout, 9020)
3970 9020
FORMAT(
' ----- DEXNEW: LACK OF INITIALISATION')
3973 SUBROUTINE dadnew(MODE,ISGN,HV,PNU,PWB,PNPI,JNPI)
3975 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
3976 * ,ampiz,ampi,amro,gamro,ama1,gama1
3977 * ,amk,amkz,amkst,gamkst
3979 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
3980 * ,ampiz,ampi,amro,gamro,ama1,gama1
3981 * ,amk,amkz,amkst,gamkst
3982 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
3983 real*4 gfermi,gv,ga,ccabib,scabib,gamel
3984 COMMON / taubmc / gampmc(30),gamper(30),nevdec(30)
3985 real*4 gampmc ,gamper
3986 COMMON / inout / inut,iout
3987 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
3988 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
3990 CHARACTER NAMES(NMODE)*31
3992 real*4 pnu(4),pwb(4),pnpi(4,9),hv(4),hhv(4)
3993 real*4 pdum1(4),pdum2(4),pdumi(4,9)
3996 real*8 swt(nmode),sswt(nmode)
3997 dimension nevraw(nmode),nevovr(nmode),nevacc(nmode)
3999 DATA pi /3.141592653589793238462643/
4018 IF (jnpi.LE.nm4)
THEN
4020 wtmax(jnpi) = pkorb(3,37+jnpi)
4026 ELSEIF(jnpi.LE.nm4)
THEN
4027 CALL dph4pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4028 ELSEIF(jnpi.LE.nm4+nm5)
THEN
4029 CALL dph5pi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4030 ELSEIF(jnpi.LE.nm4+nm5+nm6)
THEN
4031 CALL dphnpi(wt,hv,pdum1,pdum2,pdumi,jnpi)
4032 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3)
THEN
4033 inum=jnpi-nm4-nm5-nm6
4034 CALL dphspk(wt,hv,pdum1,pdum2,pdumi,inum)
4035 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2)
THEN
4036 inum=jnpi-nm4-nm5-nm6-nm3
4037 CALL dphsrk(wt,hv,pdum1,pdum2,pdumi,inum)
4041 IF(wt.GT.wtmax(jnpi)/1.2) wtmax(jnpi)=wt*1.2
4049 ELSEIF(mode.EQ. 0)
THEN
4051 IF(iwarm.EQ.0)
GOTO 902
4056 ELSEIF(jnpi.LE.nm4)
THEN
4057 CALL dph4pi(wt,hhv,pnu,pwb,pnpi,jnpi)
4058 ELSEIF(jnpi.LE.nm4+nm5)
THEN
4059 CALL dph5pi(wt,hhv,pnu,pwb,pnpi,jnpi)
4060 ELSEIF(jnpi.LE.nm4+nm5+nm6)
THEN
4061 CALL dphnpi(wt,hhv,pnu,pwb,pnpi,jnpi)
4062 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3)
THEN
4063 inum=jnpi-nm4-nm5-nm6
4064 CALL dphspk(wt,hhv,pnu,pwb,pnpi,inum)
4065 ELSEIF(jnpi.LE.nm4+nm5+nm6+nm3+nm2)
THEN
4066 inum=jnpi-nm4-nm5-nm6-nm3
4067 CALL dphsrk(wt,hhv,pnu,pwb,pnpi,inum)
4075 nevraw(jnpi)=nevraw(jnpi)+1
4076 swt(jnpi)=swt(jnpi)+wt
4079 sswt(jnpi)=sswt(jnpi)+dble(wt)**2
4083 IF(wt.GT.wtmax(jnpi)) nevovr(jnpi)=nevovr(jnpi)+1
4084 IF(rn*wtmax(jnpi).GT.wt)
GOTO 300
4086 costhe=-1.+2.*rrr(2)
4089 CALL rotor2(thet,pnu,pnu)
4090 CALL rotor3( phi,pnu,pnu)
4091 CALL rotor2(thet,pwb,pwb)
4092 CALL rotor3( phi,pwb,pwb)
4093 CALL rotor2(thet,hv,hv)
4094 CALL rotor3( phi,hv,hv)
4097 CALL rotor2(thet,pnpi(1,i),pnpi(1,i))
4098 CALL rotor3( phi,pnpi(1,i),pnpi(1,i))
4100 nevacc(jnpi)=nevacc(jnpi)+1
4102 ELSEIF(mode.EQ. 1)
THEN
4105 IF(nevraw(jnpi).EQ.0)
GOTO 500
4106 pargam=swt(jnpi)/float(nevraw(jnpi)+1)
4108 IF(nevraw(jnpi).NE.0)
4109 & error=sqrt(sswt(jnpi)/swt(jnpi)**2-1./float(nevraw(jnpi)))
4111 WRITE(iout, 7010) names(jnpi),
4112 & nevraw(jnpi),nevacc(jnpi),nevovr(jnpi),pargam,rat,error
4114 gampmc(8+jnpi-1)=rat
4115 gamper(8+jnpi-1)=error
4121 7003
FORMAT(///1x,15(5h*****)
4122 $ /,
' *', 25x,
'******** DADNEW INITIALISATION ********',9x,1h*
4124 7004
FORMAT(
' *',e20.5,5x,
'WTMAX = MAXIMUM WEIGHT ',9x,1h*/)
4126 $ /,1x,15(5h*****)/)
4127 7010
FORMAT(///1x,15(5h*****)
4128 $ /,
' *', 25x,
'******** DADNEW FINAL REPORT ******** ',9x,1h*
4129 $ /,
' *', 25x,
'CHANNEL:',a31 ,9x,1h*
4130 $ /,
' *',i20 ,5x,
'NEVRAW = NO. OF DECAYS TOTAL ',9x,1h*
4131 $ /,
' *',i20 ,5x,
'NEVACC = NO. OF DECAYS ACCEPTED ',9x,1h*
4132 $ /,
' *',i20 ,5x,
'NEVOVR = NO. OF OVERWEIGHTED EVENTS ',9x,1h*
4133 $ /,
' *',e20.5,5x,
'PARTIAL WTDTH IN GEV UNITS ',9x,1h*
4134 $ /,
' *',f20.9,5x,
'IN UNITS GFERMI**2*MASS**5/192/PI**3 ',9x,1h*
4135 $ /,
' *',f20.8,5x,
'RELATIVE ERROR OF PARTIAL WIDTH ',9x,1h*
4136 $ /,1x,15(5h*****)/)
4137 902
WRITE(iout, 9020)
4138 9020
FORMAT(
' ----- DADNEW: LACK OF INITIALISATION')
4140 903
WRITE(iout, 9030) jnpi,mode
4141 9030
FORMAT(
' ----- DADNEW: WRONG JNPI',2i5)
4146 SUBROUTINE dph4pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4151 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4152 * ,ampiz,ampi,amro,gamro,ama1,gama1
4153 * ,amk,amkz,amkst,gamkst
4155 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4156 * ,ampiz,ampi,amro,gamro,ama1,gama1
4157 * ,amk,amkz,amkst,gamkst
4158 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4159 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4160 REAL HV(4),PT(4),PN(4),PAA(4),PIM1(4),PIM2(4),PIPL(4),PMULT(4,9)
4163 real*8 uu,ff,ff1,ff2,ff3,ff4,gg1,gg2,gg3,gg4,rr
4164 DATA pi /3.141592653589793238462643/
4166 xlam(x,y,z)=sqrt(abs((x-y-z)**2-4.0*y*z))
4171 phspac=1./2**23/pi**11
4200 CALL choice(100+jnpi,rrb,ichan,prob1,prob2,prob3,
4201 $ amrop,gamrop,amrx,gamrx,amrb,gamrb)
4215 ams1=(amp1+amp2+amp3+amp4)**2
4216 ams2=(amtau-amnuta)**2
4217 alp1=atan((ams1-amrop**2)/amrop/gamrop)
4218 alp2=atan((ams2-amrop**2)/amrop/gamrop)
4219 alp=alp1+rr1*(alp2-alp1)
4220 am4sq =amrop**2+amrop*gamrop*tan(alp)
4223 $ ((am4sq-amrop**2)**2+(amrop*gamrop)**2)/(amrop*gamrop)
4224 phspac=phspac*(alp2-alp1)
4228 ams1=(amp2+amp3+amp4)**2
4230 IF (rrr(9).GT.prez)
THEN
4231 am3sq=ams1+ rr1*(ams2-ams1)
4237 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4238 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4239 alp=alp1+rr1*(alp2-alp1)
4240 am3sq =amrx**2+amrx*gamrx*tan(alp)
4243 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4251 am2sq=ams1+ rr2*(ams2-ams1)
4256 enq1=(am2sq-amp3**2+amp4**2)/(2*am2)
4257 enq2=(am2sq+amp3**2-amp4**2)/(2*am2)
4258 ppi= enq1**2-amp4**2
4259 pppi=sqrt(abs(enq1**2-amp4**2))
4260 phspac=phspac*(4*pi)*(2*pppi/am2)
4262 CALL sphera(pppi,piz)
4272 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4273 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4274 ppi = pr(4)**2-am2**2
4278 pim1(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4281 ff3=(4*pi)*(2*pr(3)/am3)
4283 exe=(pr(4)+pr(3))/am2
4284 CALL bostr3(exe,piz,piz)
4285 CALL bostr3(exe,pipl,pipl)
4288 thet =acos(-1.+2*rr3)
4290 CALL rotpol(thet,phi,pipl)
4291 CALL rotpol(thet,phi,pim1)
4292 CALL rotpol(thet,phi,piz)
4293 CALL rotpol(thet,phi,pr)
4298 pr(4)=1./(2*am4)*(am4**2+am3**2-amp1**2)
4299 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4300 ppi = pr(4)**2-am3**2
4304 pim2(4)=1./(2*am4)*(am4**2-am3**2+amp1**2)
4307 ff4=(4*pi)*(2*pr(3)/am4)
4309 exe=(pr(4)+pr(3))/am3
4310 CALL bostr3(exe,piz,piz)
4311 CALL bostr3(exe,pipl,pipl)
4312 CALL bostr3(exe,pim1,pim1)
4315 thet =acos(-1.+2*rr3)
4317 CALL rotpol(thet,phi,pipl)
4318 CALL rotpol(thet,phi,pim1)
4319 CALL rotpol(thet,phi,pim2)
4320 CALL rotpol(thet,phi,piz)
4321 CALL rotpol(thet,phi,pr)
4327 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am4**2)
4328 paa(3)= sqrt(abs(paa(4)**2-am4**2))
4329 ppi = paa(4)**2-am4**2
4330 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4331 phsp=phsp*(4*pi)*(2*paa(3)/amtau)
4335 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am4**2)
4338 IF(rrr(9).LE.0.5*prez)
THEN
4347 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4348 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4350 ams1=(amp1+amp3+amp4)**2
4353 ams2=(sqrt(am3sq)-amp1)**2
4355 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4356 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4359 am3sq=(pim1(4)+piz(4)+pipl(4))**2-(pim1(3)+piz(3)+pipl(3))**2
4360 $ -(pim1(2)+piz(2)+pipl(2))**2-(pim1(1)+piz(1)+pipl(1))**2
4362 ams1=(amp1+amp3+amp4)**2
4363 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4364 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4365 ff1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4368 ams2=(sqrt(am3sq)-amp1)**2
4370 ff3=(4*pi)*(xlam(am2**2,amp1**2,am3sq)/am3sq)
4371 ff4=(4*pi)*(xlam(am3sq,amp2**2,am4**2)/am4**2)
4374 am3sq=(pim2(4)+piz(4)+pipl(4))**2-(pim2(3)+piz(3)+pipl(3))**2
4375 $ -(pim2(2)+piz(2)+pipl(2))**2-(pim2(1)+piz(1)+pipl(1))**2
4377 ams1=(amp2+amp3+amp4)**2
4378 alp1=atan((ams1-amrx**2)/amrx/gamrx)
4379 alp2=atan((ams2-amrx**2)/amrx/gamrx)
4380 gg1=((am3sq-amrx**2)**2+(amrx*gamrx)**2)/(amrx*gamrx)
4383 ams2=(sqrt(am3sq)-amp2)**2
4385 gg3=(4*pi)*(xlam(am2**2,amp2**2,am3sq)/am3sq)
4386 gg4=(4*pi)*(xlam(am3sq,amp1**2,am4**2)/am4**2)
4390 IF ( (0.5*prez*(ff+gg)*uu+(1.0-prez)*ff*gg).GT.0.0d0)
THEN
4391 rr=ff*gg*uu/(0.5*prez*(ff+gg)*uu+(1.0-prez)*ff*gg)
4419 exe=(paa(4)+paa(3))/am4
4420 CALL bostr3(exe,piz,piz)
4421 CALL bostr3(exe,pipl,pipl)
4422 CALL bostr3(exe,pim1,pim1)
4423 CALL bostr3(exe,pim2,pim2)
4424 CALL bostr3(exe,pr,pr)
4433 CALL dam4pi(jnpi,pt,pn,pim1,pim2,piz,pipl,amplit,hv)
4434 ELSEIF (jnpi.EQ.2)
THEN
4435 CALL dam4pi(jnpi,pt,pn,pim1,pim2,pipl,piz,amplit,hv)
4437 dgamt=1/(2.*amtau)*amplit*phspac
4447 SUBROUTINE dam4pi(MNUM,PT,PN,PIM1,PIM2,PIM3,PIM4,AMPLIT,HV)
4457 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4458 * ,ampiz,ampi,amro,gamro,ama1,gama1
4459 * ,amk,amkz,amkst,gamkst
4461 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4462 * ,ampiz,ampi,amro,gamro,ama1,gama1
4463 * ,amk,amkz,amkst,gamkst
4464 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4465 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4466 REAL HV(4),PT(4),PN(4),PIM1(4),PIM2(4),PIM3(4),PIM4(4)
4467 REAL PIVEC(4),PIAKS(4),HVM(4)
4468 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FORM4,FORM5
4469 EXTERNAL form1,form2,form3,form4,form5
4470 DATA pi /3.141592653589793238462643/
4473 CALL curr_cleo(mnum,pim1,pim2,pim3,pim4,hadcur)
4476 CALL clvec(hadcur,pn,pivec)
4477 CALL claxi(hadcur,pn,piaks)
4478 CALL clnut(hadcur,brakm,hvm)
4480 brak= (gv**2+ga**2)*pt(4)*pivec(4) +2.*gv*ga*pt(4)*piaks(4)
4481 & +2.*(gv**2-ga**2)*amnuta*amtau*brakm
4482 amplit=(ccabib*gfermi)**2*brak/2.
4485 hv(i)=-(amtau*((gv**2+ga**2)*piaks(i)+2.*gv*ga*pivec(i)))
4486 & +(gv**2-ga**2)*amnuta*amtau*hvm(i)
4491 SUBROUTINE dph5pi(DGAMT,HV,PN,PAA,PMULT,JNPI)
4496 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4497 * ,ampiz,ampi,amro,gamro,ama1,gama1
4498 * ,amk,amkz,amkst,gamkst
4500 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4501 * ,ampiz,ampi,amro,gamro,ama1,gama1
4504 * ,amk,amkz,amkst,gamkst
4505 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
4506 real*4 gfermi,gv,ga,ccabib,scabib,gamel
4507 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
4508 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
4510 CHARACTER NAMES(NMODE)*31
4511 REAL HV(4),PT(4),PN(4),PAA(4),PMULT(4,9)
4512 real*4 pr(4),pi1(4),pi2(4),pi3(4),pi4(4),pi5(4)
4513 real*8 amp1,amp2,amp3,amp4,amp5,ams1,ams2,amom,gamom
4514 real*8 am5sq,am4sq,am3sq,am2sq,am5,am4,am3
4516 real*8 gg1,gg2,gg3,ff1,ff2,ff3,ff4,alp,alp1,alp2
4521 DATA pi /3.141592653589793238462643/
4527 bwign(xm,am,gamma)=xm**2/cmplx(xm**2-am**2,gamma*am)
4535 phspac=1./2**29/pi**14
4538 amp1=dcdmas(idffin(1,jnpi))
4539 amp2=dcdmas(idffin(2,jnpi))
4540 amp3=dcdmas(idffin(3,jnpi))
4541 amp4=dcdmas(idffin(4,jnpi))
4542 amp5=dcdmas(idffin(5,jnpi))
4557 ams1=(amp1+amp2+amp3+amp4+amp5)**2
4558 ams2=(amtau-amnuta)**2
4559 am5sq=ams1+ rr1*(ams2-ams1)
4561 phspac=phspac*(ams2-ams1)
4566 ams1=(amp2+amp3+amp4+amp5)**2
4568 am4sq=ams1+ rr1*(ams2-ams1)
4575 ams1=(amp2+amp3+amp4)**2
4577 alp1=atan((ams1-amom**2)/amom/gamom)
4578 alp2=atan((ams2-amom**2)/amom/gamom)
4579 alp=alp1+rr1*(alp2-alp1)
4580 am3sq =amom**2+amom*gamom*tan(alp)
4583 gg2=((am3sq-amom**2)**2+(amom*gamom)**2)/(amom*gamom)
4596 am2sq=ams1+ rr2*(ams2-ams1)
4602 enq1=(am2sq+amp3**2-amp4**2)/(2*am2)
4603 enq2=(am2sq-amp3**2+amp4**2)/(2*am2)
4604 ppi= enq1**2-amp3**2
4605 pppi=sqrt(abs(enq1**2-amp3**2))
4606 ff1=(4*pi)*(2*pppi/am2)
4608 call sphera(pppi,pi3)
4619 pr(4)=1./(2*am3)*(am3**2+am2**2-amp2**2)
4620 pr(3)= sqrt(abs(pr(4)**2-am2**2))
4621 ppi = pr(4)**2-am2**2
4625 pi2(4)=1./(2*am3)*(am3**2-am2**2+amp2**2)
4628 ff2=(4*pi)*(2*pr(3)/am3)
4630 exe=(pr(4)+pr(3))/am2
4631 call bostr3(exe,pi3,pi3)
4632 call bostr3(exe,pi4,pi4)
4635 thet =acos(-1.+2*rr3)
4637 call rotpol(thet,phi,pi2)
4638 call rotpol(thet,phi,pi3)
4639 call rotpol(thet,phi,pi4)
4645 pr(4)=1./(2*am4)*(am4**2+am3**2-amp5**2)
4646 pr(3)= sqrt(abs(pr(4)**2-am3**2))
4647 ppi = pr(4)**2-am3**2
4651 pi5(4)=1./(2*am4)*(am4**2-am3**2+amp5**2)
4654 ff3=(4*pi)*(2*pr(3)/am4)
4656 exe=(pr(4)+pr(3))/am3
4657 call bostr3(exe,pi2,pi2)
4658 call bostr3(exe,pi3,pi3)
4659 call bostr3(exe,pi4,pi4)
4662 thet =acos(-1.+2*rr3)
4664 call rotpol(thet,phi,pi2)
4665 call rotpol(thet,phi,pi3)
4666 call rotpol(thet,phi,pi4)
4667 call rotpol(thet,phi,pi5)
4673 pr(4)=1./(2*am5)*(am5**2+am4**2-amp1**2)
4674 pr(3)= sqrt(abs(pr(4)**2-am4**2))
4675 ppi = pr(4)**2-am4**2
4679 pi1(4)=1./(2*am5)*(am5**2-am4**2+amp1**2)
4682 ff4=(4*pi)*(2*pr(3)/am5)
4684 exe=(pr(4)+pr(3))/am4
4685 call bostr3(exe,pi2,pi2)
4686 call bostr3(exe,pi3,pi3)
4687 call bostr3(exe,pi4,pi4)
4688 call bostr3(exe,pi5,pi5)
4691 thet =acos(-1.+2*rr3)
4693 call rotpol(thet,phi,pi1)
4694 call rotpol(thet,phi,pi2)
4695 call rotpol(thet,phi,pi3)
4696 call rotpol(thet,phi,pi4)
4697 call rotpol(thet,phi,pi5)
4706 paa(4)=1./(2*amtau)*(amtau**2-amnuta**2+am5sq)
4707 paa(3)= sqrt(abs(paa(4)**2-am5sq))
4708 ppi = paa(4)**2-am5sq
4709 phspac=phspac*(4*pi)*(2*paa(3)/amtau)
4713 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-am5**2)
4716 phspac=phspac * gg1*gg2*gg3*ff1*ff2*ff3*ff4
4720 exe=(paa(4)+paa(3))/am5
4721 call bostr3(exe,pi1,pi1)
4722 call bostr3(exe,pi2,pi2)
4723 call bostr3(exe,pi3,pi3)
4724 call bostr3(exe,pi4,pi4)
4725 call bostr3(exe,pi5,pi5)
4733 qxn=paa(4)*pn(4)-paa(1)*pn(1)-paa(2)*pn(2)-paa(3)*pn(3)
4734 brak=2*(gv**2+ga**2)*(2*pxq*qxn+am5sq*pxn)
4735 & -6*(gv**2-ga**2)*amtau*amnuta*am5sq
4736 fompp = cabs(bwign(am3,amom,gamom))**2
4741 amplit=ccabib**2*gfermi**2/2. * brak
4742 amplit = amplit * fompp * fnorm
4745 dgamt=1/(2.*amtau)*amplit*phspac
4763 SUBROUTINE dphsrk(DGAMT,HV,PN,PR,PMULT,INUM)
4769 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4770 * ,ampiz,ampi,amro,gamro,ama1,gama1
4771 * ,amk,amkz,amkst,gamkst
4773 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4774 * ,ampiz,ampi,amro,gamro,ama1,gama1
4775 * ,amk,amkz,amkst,gamkst
4777 REAL HV(4),PT(4),PN(4),PR(4),PKC(4),PKZ(4),QQ(4),PMULT(4,9)
4779 DATA pi /3.141592653589793238462643/
4783 phspac=1./2**11/pi**5
4791 ams2=(amtau-amnuta)**2
4794 amx2=ams1+ rr1(1)*(ams2-ams1)
4796 phspac=phspac*(ams2-ams1)
4814 pn(4)=1./(2*amtau)*(amtau**2+amnuta**2-amx**2)
4815 pn(3)=-sqrt(abs((pn(4)-amnuta)*(pn(4)+amnuta)))
4819 pr(4)=1./(2*amtau)*(amtau**2-amnuta**2+amx**2)
4821 phspac=phspac*(4*pi)*(2*pr(3)/amtau)
4824 enq1=(amx2+amk**2-amkz**2)/(2.*amx)
4825 enq2=(amx2-amk**2+amkz**2)/(2.*amx)
4826 pppi=sqrt(abs(enq1-amk)*(enq1+amk))
4827 phspac=phspac*(4*pi)*(2*pppi/amx)
4829 CALL sphera(pppi,pkc)
4835 exe=(pr(4)+pr(3))/amx
4837 CALL bostr3(exe,pkc,pkc)
4838 CALL bostr3(exe,pkz,pkz)
4841 CALL dam2pi(3,pt,pn,pkc,pkz,amplit,hv)
4842 dgamt=1/(2.*amtau)*amplit*phspac
4854 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
4855 * ,ampiz,ampi,amro,gamro,ama1,gama1
4856 * ,amk,amkz,amkst,gamkst
4858 real*4 amtau,amnuta,amel,amnue,ammu,amnumu
4859 * ,ampiz,ampi,amro,gamro,ama1,gama1
4860 * ,amk,amkz,amkst,gamkst
4863 fpirk=cabs(fpikm(w,amk,amkz))**2
4866 COMPLEX FUNCTION fpikmk(W,XM1,XM2)
4871 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
4876 IF (init.EQ.0 )
THEN
4889 fpikmk=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
4898 SUBROUTINE dwrph(KTO,PHX)
4902 IMPLICIT REAL*8 (a-h,o-z)
4913 IF (qhot(4).GT.1.e-5)
CALL dwluph(kto,qhot)
4916 SUBROUTINE dwluph(KTO,PHOT)
4927 COMMON /taupos/ np1,np2
4930 IF (phot(4).LE.0.0)
RETURN
4933 IF((kto.EQ. 1).OR.(kto.EQ.11))
THEN
4940 IF(ktos.GT.10) ktos=ktos-10
4942 CALL tralo4(ktos,phot,phot,am)
4943 CALL filhep(0,1,22,nps,nps,0,0,phot,0.0,.true.)
4948 SUBROUTINE dwluel(KTO,ISGN,PNU,PWB,PEL,PNE)
4958 REAL PNU(4),PWB(4),PEL(4),PNE(4)
4959 COMMON /taupos/ np1,np2
4969 CALL tralo4(kto,pnu,pnu,am)
4970 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
4973 CALL tralo4(kto,pwb,pwb,am)
4977 CALL tralo4(kto,pel,pel,am)
4978 CALL filhep(0,1,11*isgn,nps,nps,0,0,pel,am,.false.)
4981 CALL tralo4(kto,pne,pne,am)
4982 CALL filhep(0,1,-12*isgn,nps,nps,0,0,pne,am,.true.)
4986 SUBROUTINE dwlumu(KTO,ISGN,PNU,PWB,PMU,PNM)
4996 REAL PNU(4),PWB(4),PMU(4),PNM(4)
4997 COMMON /taupos/ np1,np2
5007 CALL tralo4(kto,pnu,pnu,am)
5008 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5011 CALL tralo4(kto,pwb,pwb,am)
5015 CALL tralo4(kto,pmu,pmu,am)
5016 CALL filhep(0,1,13*isgn,nps,nps,0,0,pmu,am,.false.)
5019 CALL tralo4(kto,pnm,pnm,am)
5020 CALL filhep(0,1,-14*isgn,nps,nps,0,0,pnm,am,.true.)
5024 SUBROUTINE dwlupi(KTO,ISGN,PPI,PNU)
5035 COMMON /taupos/ np1,np2
5045 CALL tralo4(kto,pnu,pnu,am)
5046 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5049 CALL tralo4(kto,ppi,ppi,am)
5050 CALL filhep(0,1,-211*isgn,nps,nps,0,0,ppi,am,.true.)
5054 SUBROUTINE dwluro(KTO,ISGN,PNU,PRHO,PIC,PIZ)
5064 REAL PNU(4),PRHO(4),PIC(4),PIZ(4)
5065 COMMON /taupos/ np1,np2
5075 CALL tralo4(kto,pnu,pnu,am)
5076 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5079 CALL tralo4(kto,prho,prho,am)
5080 CALL filhep(0,2,-213*isgn,nps,nps,0,0,prho,am,.true.)
5083 CALL tralo4(kto,pic,pic,am)
5084 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pic,am,.true.)
5087 CALL tralo4(kto,piz,piz,am)
5088 CALL filhep(0,1,111,-2,-2,0,0,piz,am,.true.)
5092 SUBROUTINE dwluaa(KTO,ISGN,PNU,PAA,PIM1,PIM2,PIPL,JAA)
5103 REAL PNU(4),PAA(4),PIM1(4),PIM2(4),PIPL(4)
5104 COMMON /taupos/ np1,np2
5114 CALL tralo4(kto,pnu,pnu,am)
5115 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5118 CALL tralo4(kto,paa,paa,am)
5119 CALL filhep(0,1,-20213*isgn,nps,nps,0,0,paa,am,.true.)
5127 CALL tralo4(kto,pim2,pim2,am)
5128 CALL filhep(0,1,-211*isgn,-1,-1,0,0,pim2,am,.true.)
5131 CALL tralo4(kto,pim1,pim1,am)
5132 CALL filhep(0,1,-211*isgn,-2,-2,0,0,pim1,am,.true.)
5135 CALL tralo4(kto,pipl,pipl,am)
5136 CALL filhep(0,1, 211*isgn,-3,-3,0,0,pipl,am,.true.)
5138 ELSE IF (jaa.EQ.2)
THEN
5143 CALL tralo4(kto,pim2,pim2,am)
5144 CALL filhep(0,1,111,-1,-1,0,0,pim2,am,.true.)
5147 CALL tralo4(kto,pim1,pim1,am)
5148 CALL filhep(0,1,111,-2,-2,0,0,pim1,am,.true.)
5151 CALL tralo4(kto,pipl,pipl,am)
5152 CALL filhep(0,1,-211*isgn,-3,-3,0,0,pipl,am,.true.)
5158 SUBROUTINE dwlukk (KTO,ISGN,PKK,PNU)
5168 COMMON /taupos/ np1,np2
5178 CALL tralo4 (kto,pnu,pnu,am)
5179 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5182 CALL tralo4 (kto,pkk,pkk,am)
5183 CALL filhep(0,1,-321*isgn,nps,nps,0,0,pkk,am,.true.)
5187 SUBROUTINE dwluks(KTO,ISGN,PNU,PKS,PKK,PPI,JKST)
5188 COMMON / taukle / bra1,brk0,brk0b,brks
5189 real*4 bra1,brk0,brk0b,brks
5199 REAL PNU(4),PKS(4),PKK(4),PPI(4),XIO(1)
5200 COMMON /taupos/ np1,np2
5210 CALL tralo4(kto,pnu,pnu,am)
5211 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5214 CALL tralo4(kto,pks,pks,am)
5215 CALL filhep(0,1,-323*isgn,nps,nps,0,0,pks,am,.true.)
5223 CALL tralo4(kto,ppi,ppi,am)
5224 CALL filhep(0,1,-211*isgn,-1,-1,0,0,ppi,am,.true.)
5227 IF (isgn.EQ.-1) bran=brk0
5230 IF(xio(1).GT.bran)
THEN
5236 CALL tralo4(kto,pkk,pkk,am)
5237 CALL filhep(0,1,k0type,-2,-2,0,0,pkk,am,.true.)
5239 ELSE IF(jkst.EQ.20)
THEN
5244 CALL tralo4(kto,ppi,ppi,am)
5245 CALL filhep(0,1,111,-1,-1,0,0,ppi,am,.true.)
5248 CALL tralo4(kto,pkk,pkk,am)
5249 CALL filhep(0,1,-321*isgn,-2,-2,0,0,pkk,am,.true.)
5255 SUBROUTINE dwlnew(KTO,ISGN,PNU,PWB,PNPI,MODE)
5265 parameter(nmode=15,nm1=0,nm2=1,nm3=8,nm4=2,nm5=1,nm6=3)
5266 COMMON / taudcd /idffin(9,nmode),mulpik(nmode)
5268 COMMON /taupos/ np1,np2
5269 CHARACTER NAMES(NMODE)*31
5270 REAL PNU(4),PWB(4),PNPI(4,9)
5282 CALL tralo4(kto,pnu,pnu,am)
5283 CALL filhep(0,1,16*isgn,nps,nps,0,0,pnu,am,.true.)
5286 CALL tralo4(kto,pwb,pwb,am)
5287 CALL filhep(0,1,-24*isgn,nps,nps,0,0,pwb,am,.true.)
5294 kfpi=lunpik(idffin(i,jnpi),-isgn)
5300 CALL tralo4(kto,ppi,ppi,am)
5301 CALL filhep(0,1,kfpi,-i,-i,0,0,ppi,am,.true.)
5312 IMPLICIT REAL*8 (a-h,o-z)
5314 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5316 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5328 aaa=pp(4)**2-pp(3)**2-pp(2)**2-pp(1)**2
5329 IF(aaa.NE.0.0) aaa=aaa/sqrt(abs(aaa))
5338 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5339 DATA pi /3.141592653589793238462643d0/
5341 IF(abs(y).LT.abs(x))
THEN
5343 IF(x.LE.0d0) the=pi-the
5345 the=acos(x/sqrt(x**2+y**2))
5356 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5357 DATA pi /3.141592653589793238462643d0/
5359 IF(abs(y).LT.abs(x))
THEN
5361 IF(x.LE.0d0) the=pi-the
5363 the=acos(x/sqrt(x**2+y**2))
5365 IF(y.LT.0d0) the=2d0*pi-the
5368 SUBROUTINE rotod1(PH1,PVEC,QVEC)
5373 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5374 dimension pvec(4),qvec(4),rvec(4)
5382 qvec(2)= cs*rvec(2)-sn*rvec(3)
5383 qvec(3)= sn*rvec(2)+cs*rvec(3)
5387 SUBROUTINE rotod2(PH1,PVEC,QVEC)
5392 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5393 dimension pvec(4),qvec(4),rvec(4)
5400 qvec(1)= cs*rvec(1)+sn*rvec(3)
5402 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5406 SUBROUTINE rotod3(PH1,PVEC,QVEC)
5411 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5413 dimension pvec(4),qvec(4),rvec(4)
5419 qvec(1)= cs*rvec(1)-sn*rvec(2)
5420 qvec(2)= sn*rvec(1)+cs*rvec(2)
5424 SUBROUTINE bostr3(EXE,PVEC,QVEC)
5430 real*4 pvec(4),qvec(4),rvec(4)
5443 SUBROUTINE bostd3(EXE,PVEC,QVEC)
5449 IMPLICIT DOUBLE PRECISION (a-h,o-z)
5450 dimension pvec(4),qvec(4),rvec(4)
5464 SUBROUTINE rotor1(PH1,PVEC,QVEC)
5469 real*4 pvec(4),qvec(4),rvec(4)
5477 qvec(2)= cs*rvec(2)-sn*rvec(3)
5478 qvec(3)= sn*rvec(2)+cs*rvec(3)
5481 SUBROUTINE rotor2(PH1,PVEC,QVEC)
5486 IMPLICIT REAL*4(a-h,o-z)
5487 real*4 pvec(4),qvec(4),rvec(4)
5494 qvec(1)= cs*rvec(1)+sn*rvec(3)
5496 qvec(3)=-sn*rvec(1)+cs*rvec(3)
5499 SUBROUTINE rotor3(PHI,PVEC,QVEC)
5504 real*4 pvec(4),qvec(4),rvec(4)
5510 qvec(1)= cs*rvec(1)-sn*rvec(2)
5511 qvec(2)= sn*rvec(1)+cs*rvec(2)
5515 SUBROUTINE spherd(R,X)
5520 real*8 r,x(4),pi,costh,sinth
5522 DATA pi /3.141592653589793238462643d0/
5526 sinth=sqrt(1 -costh**2)
5527 x(1)=r*sinth*cos(2*pi*rrr(2))
5528 x(2)=r*sinth*sin(2*pi*rrr(2))
5532 SUBROUTINE rotpox(THET,PHI,PP)
5533 IMPLICIT REAL*8 (a-h,o-z)
5539 CALL rotod2(thet,pp,pp)
5540 CALL rotod3( phi,pp,pp)
5543 SUBROUTINE sphera(R,X)
5551 DATA pi /3.141592653589793238462643/
5555 sinth=sqrt(1.-costh**2)
5556 x(1)=r*sinth*cos(2*pi*rrr(2))
5557 x(2)=r*sinth*sin(2*pi*rrr(2))
5561 SUBROUTINE rotpol(THET,PHI,PP)
5568 CALL rotor2(thet,pp,pp)
5569 CALL rotor3( phi,pp,pp)
5572 SUBROUTINE ranmar(RVEC,LENV)
5588 COMMON / inout / inut,iout
5590 common/raset1/u(97),c,i97,j97
5591 parameter(modcns=1000000000)
5592 DATA ntot,ntot2,ijkl/-1,0,0/
5594 IF (ntot .GE. 0)
GO TO 50
5603 entry rmarin(ijklin, ntotin,ntot2n)
5612 ntot = max(ntotin,0)
5613 ntot2= max(ntot2n,0)
5618 kl = ijkl - 30082*ij
5619 i = mod(ij/177, 177) + 2
5620 j = mod(ij, 177) + 2
5621 k = mod(kl/169, 178) + 1
5623 WRITE(iout,201) ijkl,ntot,ntot2
5624 201
FORMAT(1x,
' RANMAR INITIALIZED: ',i10,2x,2i10)
5629 m = mod(mod(i*j,179)*k, 179)
5633 l = mod(53*l+1, 169)
5634 IF (mod(l*m,64) .GE. 32) s = s+t
5639 4 twom24 = 0.5*twom24
5641 cd = 7654321.*twom24
5642 cm = 16777213.*twom24
5647 DO 45 loop2= 1, ntot2+1
5649 IF (loop2 .EQ. ntot2+1) now=ntot
5650 IF (now .GT. 0)
THEN
5651 WRITE (iout,
'(A,I15)')
' RMARIN SKIPPING OVER ',now
5652 DO 40 idum = 1, ntot
5654 IF (uni .LT. 0.) uni=uni+1.
5657 IF (i97 .EQ. 0) i97=97
5659 IF (j97 .EQ. 0) j97=97
5661 IF (c .LT. 0.) c=c+cm
5665 IF (kalled .EQ. 1)
RETURN
5669 DO 100 ivec= 1, lenv
5671 IF (uni .LT. 0.) uni=uni+1.
5674 IF (i97 .EQ. 0) i97=97
5676 IF (j97 .EQ. 0) j97=97
5678 IF (c .LT. 0.) c=c+cm
5680 IF (uni .LT. 0.) uni=uni+1.
5682 IF (uni .EQ. 0.)
THEN
5685 IF (uni .EQ. 0.) uni= twom24*twom24
5690 IF (ntot .GE. modcns)
THEN
5692 ntot = ntot - modcns
5696 entry rmarut(ijklut,ntotut,ntot2t)
5704 IMPLICIT REAL*8(a-h,o-z)
5707 IF(x .LT.-1.0)
GO TO 1
5708 IF(x .LE. 0.5)
GO TO 2
5709 IF(x .EQ. 1.0)
GO TO 3
5710 IF(x .LE. 2.0)
GO TO 4
5714 z=z-0.5* log(abs(x))**2
5720 3 dilogt=1.64493406684822
5724 z=1.64493406684822 - log(x)* log(abs(t))
5725 5 y=2.66666666666666 *t+0.66666666666666
5726 b= 0.00000 00000 00001
5727 a=y*b +0.00000 00000 00004
5728 b=y*a-b+0.00000 00000 00011
5729 a=y*b-a+0.00000 00000 00037
5730 b=y*a-b+0.00000 00000 00121
5731 a=y*b-a+0.00000 00000 00398
5732 b=y*a-b+0.00000 00000 01312
5733 a=y*b-a+0.00000 00000 04342
5734 b=y*a-b+0.00000 00000 14437
5735 a=y*b-a+0.00000 00000 48274
5736 b=y*a-b+0.00000 00001 62421
5737 a=y*b-a+0.00000 00005 50291
5738 b=y*a-b+0.00000 00018 79117
5739 a=y*b-a+0.00000 00064 74338
5740 b=y*a-b+0.00000 00225 36705
5741 a=y*b-a+0.00000 00793 87055
5742 b=y*a-b+0.00000 02835 75385
5743 a=y*b-a+0.00000 10299 04264
5744 b=y*a-b+0.00000 38163 29463
5745 a=y*b-a+0.00001 44963 00557
5746 b=y*a-b+0.00005 68178 22718
5747 a=y*b-a+0.00023 20021 96094
5748 b=y*a-b+0.00100 16274 96164
5749 a=y*b-a+0.00468 63619 59447
5750 b=y*a-b+0.02487 93229 24228
5751 a=y*b-a+0.16607 30329 27855
5752 a=y*a-b+1.93506 43008 6996