1 SUBROUTINE tauola(MODE,KEYPOL)
11 parameter(nmxhep=4000)
13 INTEGER nevhep,nhep,isthep,idhep,jmohep,
32 COMMON /taupos/ np1, np2
33 REAL*4 PHOI(4),PHOF(4)
34 double precision Q1(4),Q2(4),P1(4),P2(4),P3(4),P4(4)
35 COMMON / momdec / q1,q2,p1,p2,p3,p4
37 COMMON /libra/ jak1,jak2,itdkrc,ifphot,ifhadm,ifhadp
40 common /pseudocoup/ csc,ssc
43 COMMON / inout / inut,iout
46 dimension pol1(4), pol2(4)
47 double precision POL1x(4), POL2x(4)
49 DATA pol1 /0.0,0.0,0.0,0.0/
50 DATA pol2 /0.0,0.0,0.0,0.0/
51 DATA pi /3.141592653589793238462643d0/
59 COMMON /isons_tau/ison(2)
85 betah=sqrt(1d0-4*xmtau**2/xmh**2)
89 IF (ifphot.EQ.1)
CALL phoini
90 CALL inietc(jak1,jak2,itdkrc,ifphot)
100 WRITE(iout,7001) pol,psi,ion(1),ion(2),ion(3)
101 ELSEIF(mode.EQ.0)
THEN 106 CALL phyfix(nstop,nstart)
119 IF(abs(idhep(i)).EQ.kftau.AND.isthep(i).EQ.1.AND.
120 $ (isthep(i).GE.125.OR.isthep(i).LT.120))
THEN 122 DO WHILE (abs(idhep(imoth)).EQ.kftau)
123 imoth=jmohep(1,imoth)
125 IF (isthep(imoth).EQ.3.OR.
126 $ (isthep(imoth).GE.120.AND.isthep(imoth).LE.125))
THEN 128 IF (idhep(j).EQ.idhep(imoth).AND.
129 $ jmohep(1,j).EQ.imoth.AND.
130 $ isthep(j).EQ.2)
THEN 140 IF(jmoth.EQ.imother(ii))
GOTO 9999
160 IF (idhep(jmohep(1,im0)).EQ.idhep(im0)) im0=jmohep(1,im0)
163 IF (isthep(i).EQ.3.OR.
164 $ (isthep(i).GE.120.AND.isthep(i).LE.125))
THEN 168 DO WHILE (idhep(imoth).EQ.idhep(i).OR.abs(idhep(imoth)).EQ.kftau
169 imoth=jmohep(1,imoth)
171 IF ((imoth.EQ.im0.OR.imoth.EQ.im).AND.isel.EQ.-1)
THEN 174 ELSEIF ((imoth.EQ.im0.OR.imoth.EQ.im).AND.isel.EQ.0)
THEN 176 ELSEIF ((imoth.NE.im0.AND.imoth.NE.im).AND.isel.EQ.0)
THEN 201 IF(idhep(i).EQ.-kftau.OR.idhep(i).EQ.-kfnue) ncount=ncount+1
202 IF(idhep(i).EQ. kftau.OR.idhep(i).EQ. kfnue) ncount=ncount+1
209 DO i=max(np1+1,ison(1)),ison(2)
211 IF(idhep(i).EQ.-kftau.OR.idhep(i).EQ.-kfnue) np1=i
214 DO i=max(np2+1,ison(1)),ison(2)
216 IF(idhep(i).EQ. kftau.OR.idhep(i).EQ. kfnue) np2=i
231 IF(idhep(im).EQ.kfhiggs(1).OR.idhep(im).EQ.kfhiggs(2).OR.
232 $ idhep(im).EQ.kfhiggs(3))
THEN 233 IF(rrr(1).LT.0.5)
THEN 240 ELSEIF((idhep(im).EQ.kfz0).OR.(idhep(im).EQ.kfgam))
THEN 248 polz0=plzapx(.true.,im,np1,np2)
249 IF(rrr(1).LT.polz0)
THEN 256 ELSEIF(idhep(np1).EQ.-idhep(np2))
THEN 257 polz0=plzapx(.true.,im,np1,np2)
258 IF(rrr(1).LT.polz0)
THEN 265 ELSEIF(abs(idhep(im)).EQ.kfhigch)
THEN 275 IF(idhep(im).EQ.kfhiggs(1).OR.idhep(im).EQ.kfhiggs(2).OR.
276 $ idhep(im).EQ.kfhiggs(3))
THEN 277 IF(idhep(np1).EQ.-kftau .AND.
278 $ (jdahep(1,np1).LE.np1.OR.jdahep(1,np1).GT.nhep) .AND.
279 $ idhep(np2).EQ. kftau .AND.
280 $ (jdahep(1,np2).LE.np2.OR.jdahep(1,np2).GT.nhep)
282 IF (idhep(im).EQ.kfhiggs(1))
THEN 284 ELSEIF (idhep(im).EQ.kfhiggs(2))
THEN 286 ELSEIF (idhep(im).EQ.kfhiggs(3))
THEN 289 WRITE(*,*)
'Warning from TAUOLA:' 290 WRITE(*,*)
'I stop this run, wrong IDHEP(IM)=',idhep(im
293 CALL spinhiggs(im,np1,np2,ifpseudo,pol1,pol2)
294 IF (ifphot.EQ.1)
CALL photos(im)
300 IF(idhep(np1).EQ.-kftau.AND.
301 $ (jdahep(1,np1).LE.np1.OR.jdahep(1,np1).GT.nhep))
THEN 304 IF (ifphot.EQ.1)
CALL photos(np1)
308 IF(idhep(np2).EQ. kftau.AND.
309 $ (jdahep(1,np2).LE.np2.OR.jdahep(1,np2).GT.nhep))
THEN 312 IF (ifphot.EQ.1)
CALL photos(np2)
317 IF (ncount.GT.0)
GOTO 666
320 ELSEIF(mode.EQ.1)
THEN 324 CALL dekay(100,pol1x)
328 7001
FORMAT(///1x,15(5h*****)
329 $ /,
' *', 25x,
'*****TAUOLA UNIVERSAL INTERFACE: ******',9x,1h*,
330 $ /,
' *', 25x,
'*****VERSION 1.21, September 2005******',9x,1h*,
331 $ /,
' *', 25x,
'**AUTHORS: P. Golonka, B. Kersevan, ***',9x,1h*,
332 $ /,
' *', 25x,
'**T. Pierzchala, E. Richter-Was, ******',9x,1h*,
333 $ /,
' *', 25x,
'****** Z. Was, M. Worek ***************',9x,1h*,
334 $ /,
' *', 25x,
'**USEFUL DISCUSSIONS, IN PARTICULAR ***',9x,1h*,
335 $ /,
' *', 25x,
'*WITH C. Biscarat and S. Slabospitzky**',9x,1h*,
336 $ /,
' *', 25x,
'****** are warmly acknowledged ********',9x,1h*,
337 $ /,
' *', 25x,
' ',9x,1h*,
338 $ /,
' *', 25x,
'********** INITIALIZATION ************',9x,1h*,
339 $ /,
' *',f20.5,5x,
'tau polarization switch must be 1 or 0 ',9x,1h*,
340 $ /,
' *',f20.5,5x,
'Higs scalar/pseudo mix CERN-TH/2003-166',9x,1h*,
341 $ /,
' *',i10, 15x,
'PI0 decay switch must be 1 or 0 ',9x,1h*,
342 $ /,
' *',i10, 15x,
'ETA decay switch must be 1 or 0 ',9x,1h*,
343 $ /,
' *',i10, 15x,
'K0S decay switch must be 1 or 0 ',9x,1h*,
346 7002
FORMAT(///1x,15(5h*****)
347 $ /,
' *', 25x,
'*****TAUOLA UNIVERSAL INTERFACE: ******',9x,1h*,
348 $ /,
' *', 25x,
'*****VERSION 1.21, September2005 ******',9x,1h*,
349 $ /,
' *', 25x,
'**AUTHORS: P. Golonka, B. Kersevan, ***',9x,1h*,
350 $ /,
' *', 25x,
'**T. Pierzchala, E. Richter-Was, ******',9x,1h*,
351 $ /,
' *', 25x,
'****** Z. Was, M. Worek ***************',9x,1h*,
352 $ /,
' *', 25x,
'**USEFUL DISCUSSIONS, IN PARTICULAR ***',9x,1h*,
353 $ /,
' *', 25x,
'*WITH C. Biscarat and S. Slabospitzky**',9x,1h*,
354 $ /,
' *', 25x,
'****** are warmly acknowledged ********',9x,1h*,
355 $ /,
' *', 25x,
'****** END OF MODULE OPERATION ********',9x,1h*,
360 SUBROUTINE spinhiggs(IM,NP1,NP2,IFPSEUDO,Pol1,Pol2)
362 REAL*8 HH1,HH2,wthiggs
363 dimension pol1(4), pol2(4),hh1(4),hh2(4), rrr(1)
371 wt=wthiggs(ifpseudo,hh1,hh2)
372 IF (rrr(1).GT.wt)
GOTO 10
378 FUNCTION wthiggs(IFPSEUDO,HH1,HH2)
380 common /pseudocoup/ csc,ssc
383 REAL*8 HH1(4),HH2(4),R(4,4),wthiggs
397 r(1,1)=(csc**2-ssc**2)/(csc**2+ssc**2)
398 r(2,2)=(csc**2-ssc**2)/(csc**2+ssc**2)
399 r(1,2)=2*csc*ssc/(csc**2+ssc**2)
400 r(2,1)=-2*csc*ssc/(csc**2+ssc**2)
410 wthiggs=wthiggs+r(k,l)*hh1(k)*hh2(l)
416 FUNCTION plzapx(HOPEin,IM0,NP1,NP2)
421 REAL*8 PLZAPX,PLZAP0,SVAR,COSTHE,sini,sfin,ZPROP2,
422 $ P1(4),P2(4),Q1(4),Q2(4),QQ(4),PH(4),PD1(4),PD2(4),
423 $ PQ1(4),PQ2(4),PB(4),PA(4)
428 parameter(nmxhep=4000)
430 INTEGER nevhep,nhep,isthep,idhep,jmohep,
452 COMMON /isons_tau/ison(2)
463 IF (idhep(im0).EQ.idhep(im00)) im=jmohep(1,im0)
472 IF (isthep(im00).EQ.120)
THEN 480 IF (imo1.EQ.0.AND.imo2.EQ.0)
THEN 483 IF (idhep(imo1).EQ.idhep(np1)) imo1=jmohep(1,imo1)
487 IF (idhep(imo2).EQ.idhep(np2)) imo2=jmohep(1,imo2)
492 ELSEIF (idhep(im).NE.22.AND.idhep(im).NE.23)
THEN 495 IF (idhep(imo1).EQ.idhep(np1)) imo1=jmohep(1,imo1)
499 IF (idhep(imo2).EQ.idhep(np2)) imo2=jmohep(1,imo2)
506 IF (imo1.EQ.0) hope=.false.
507 IF (imo2.EQ.0) hope=.false.
508 IF (imo1.EQ.imo2) hope=.false.
518 IF (im.EQ.jmohep(1,im0).AND.
519 $ (idhep(im).EQ.22.OR.idhep(im).EQ.23))
THEN 526 CALL bostdq( 1,pa, q1, q1)
527 CALL bostdq( 1,pa, q2, q2)
528 CALL bostdq(-1,pb, q1, q1)
529 CALL bostdq(-1,pb, q2, q2)
535 IF (hope) p1(i)=phep(i,imo1)
536 IF (hope) p2(i)=phep(i,imo2)
544 IF (hope) idfp1=idhep(imo1)
545 IF (hope) idfp2=idhep(imo2)
547 svar=qq(4)**2-qq(3)**2-qq(2)**2-qq(1)**2
567 IF (iffull.EQ.1) inbr=np1
568 IF (idhep(jmohep(1,inbr)).EQ.idhep(inbr)) inbr=jmohep(1,inbr)
570 IF(nx1.EQ.0.OR.nx2.EQ.0)
THEN 574 IF(jmohep(1,inbr-k).EQ.jmohep(1,inbr))
THEN 583 IF(jmohep(1,k).EQ.jmohep(1,inbr))
THEN 593 IF (abs(idfp1).GE.20.AND.abs(idfp2).GE.20) hope=.false.
595 IF (abs(idfp1).LE.20.AND.abs(idfp2).LE.20.AND.idfp1+idfp2.NE.0)
605 IF (abs(idfp1).LT.20) ide= idfp1
606 IF (abs(idfp2).LT.20) ide=-idfp2
626 iflav=min(abs(idfp1),abs(idfp2))
636 IF (abs(idfp1).GE.20)
THEN 639 IF (abs(idp).EQ.iflav)
THEN 647 IF (abs(idfp2).GE.20)
THEN 650 IF (abs(idp).EQ.iflav)
THEN 660 IF (abs(idfp1).GE.20)
THEN 664 xm1=abs((pd1(4)+ph(4))**2-(pd1(3)+ph(3))**2
665 $ -(pd1(2)+ph(2))**2-(pd1(1)+ph(1))**2)
666 xm2=abs((pd2(4)+ph(4))**2-(pd2(3)+ph(3))**2
667 $ -(pd2(2)+ph(2))**2-(pd2(1)+ph(1))**2)
682 IF (abs(idfp2).GE.20)
THEN 686 xm1=abs((pd1(4)+ph(4))**2-(pd1(3)+ph(3))**2
687 $ -(pd1(2)+ph(2))**2-(pd1(1)+ph(1))**2)
688 xm2=abs((pd2(4)+ph(4))**2-(pd2(3)+ph(3))**2
689 $ -(pd2(2)+ph(2))**2-(pd2(1)+ph(1))**2)
706 IF (idhep(jmohep(1,np1)).EQ.idhep(np1)) nph1=jmohep(1,np1)
707 IF (idhep(jmohep(1,np2)).EQ.idhep(np2)) nph2=jmohep(1,np2)
711 IF (abs(idhep(k)).NE.iflav.AND.k.NE.im.AND.
713 $ k.NE.nph1.AND.k.NE.nph2)
THEN 715 IF(idhep(k).EQ.22.AND.iffull.EQ.1)
THEN 719 xm1=abs((pd1(4)-ph(4))**2-(pd1(3)-ph(3))**2
720 $ -(pd1(2)-ph(2))**2-(pd1(1)-ph(1))**2)
721 xm2=abs((pd2(4)-ph(4))**2-(pd2(3)-ph(3))**2
722 $ -(pd2(2)-ph(2))**2-(pd2(1)-ph(1))**2)
723 xm3=abs((pq1(4)+ph(4))**2-(pq1(3)+ph(3))**2
724 $ -(pq1(2)+ph(2))**2-(pq1(1)+ph(1))**2)
725 xm4=abs((pq2(4)+ph(4))**2-(pq2(3)+ph(3))**2
726 $ -(pq2(2)+ph(2))**2-(pq2(1)+ph(1))**2)
729 sini=abs((pd1(4)+pd2(4)-ph(4))**2-(pd1(3)+pd2(3)-ph(3))**2
730 $ -(pd1(2)+pd2(2)-ph(2))**2-(pd1(1)+pd2(1)-ph(1))**2)
731 sfin=abs((pd1(4)+pd2(4) )**2-(pd1(3)+pd2(3) )**2
732 $ -(pd1(2)+pd2(2) )**2-(pd1(1)+pd2(1) )**2)
742 xm=min(xm1,xm2,xm3,xm4)
747 ELSEIF (xm2.EQ.xm)
THEN 751 ELSEIF (xm3.EQ.xm)
THEN 764 xm1=abs((pd1(4)-ph(4))**2-(pd1(3)-ph(3))**2
765 $ -(pd1(2)-ph(2))**2-(pd1(1)-ph(1))**2)
766 xm2=abs((pd2(4)-ph(4))**2-(pd2(3)-ph(3))**2
767 $ -(pd2(2)-ph(2))**2-(pd2(1)-ph(1))**2)
797 IF (abs(idhep(im0)).EQ.22.OR.abs(idhep(im0)).EQ.23)
THEN 799 IF(abs(idhep(k)).EQ.22)
THEN 806 xm3=abs((pq1(4)+ph(4))**2-(pq1(3)+ph(3))**2
807 $ -(pq1(2)+ph(2))**2-(pq1(1)+ph(1))**2)
808 xm4=abs((pq2(4)+ph(4))**2-(pq2(3)+ph(3))**2
809 $ -(pq2(2)+ph(2))**2-(pq2(1)+ph(1))**2)
832 CALL angulu(pd1,pd2,q1,q2,costhe)
834 plzapx=plzap0(ide,idfq1,svar,costhe)
837 SUBROUTINE angulu(PD1,PD2,Q1,Q2,COSTHE)
838 REAL*8 PD1(4),PD2(4),Q1(4),Q2(4),COSTHE,P(4),QQ(4),QT(4)
844 xm1=abs(pd1(4)**2-pd1(3)**2-pd1(2)**2-pd1(1)**2)
845 xm2=abs(pd2(4)**2-pd2(3)**2-pd2(2)**2-pd2(1)**2)
863 xmqq=sqrt(qq(4)**2-qq(3)**2-qq(2)**2-qq(1)**2)
865 qtxqq=qt(4)*qq(4)-qt(3)*qq(3)-qt(2)*qq(2)-qt(1)*qq(1)
867 qt(k)=qt(k)-qq(k)*qtxqq/xmqq**2
870 pxqq=p(4)*qq(4)-p(3)*qq(3)-p(2)*qq(2)-p(1)*qq(1)
872 p(k)=p(k)-qq(k)*pxqq/xmqq**2
875 pxp =sqrt(p(1)**2+p(2)**2+p(3)**2-p(4)**2)
876 qtxqt=sqrt(qt(3)**2+qt(2)**2+qt(1)**2-qt(4)**2)
877 pxqt =p(3)*qt(3)+p(2)*qt(2)+p(1)*qt(1)-p(4)*qt(4)
878 costhe=pxqt/pxp/qtxqt
882 FUNCTION plzap0(IDE,IDF,SVAR,COSTH0)
885 REAL*8 PLZAP0,SVAR,COSTHE,COSTH0
892 CALL initwk(ide,idf,svar)
894 CALL initwk(-ide,-idf,svar)
896 plzap0=t_born(0,svar,costhe,1d0,1d0)
897 $ /(t_born(0,svar,costhe,1d0,1d0)+t_born(0,svar,costhe,-1d0,-1d0))
901 FUNCTION t_born(MODE,SVAR,COSTHE,TA,TB)
912 IMPLICIT REAL*8(a-h,o-z)
913 COMMON / t_beampm / ene ,amin,amfin,ide,idf
914 REAL*8 ENE ,AMIN,AMFIN
915 COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
916 & ,xupgi ,xupzi ,xupgf ,xupzf
917 & ,ndiag0,ndiaga,keya,keyz
918 & ,itce,jtce,itcf,jtcf,kolor
919 REAL*8 SS,POLN,T3E,QE,T3F,QF
920 & ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2)
923 COMMON / t_gswprm /swsq,amw,amz,amh,amtop,gammz
924 REAL*8 SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
930 COMPLEX*16 ABORN(2,2),APHOT(2,2),AZETT(2,2)
931 COMPLEX*16 XUPZFP(2),XUPZIP(2)
932 COMPLEX*16 ABORNM(2,2),APHOTM(2,2),AZETTM(2,2)
933 COMPLEX*16 PROPA,PROPZ
935 COMPLEX*16 XUPF,XUPI,XFF(4),XFEM,XFOTA,XRHO,XKE,XKF,XKEF
936 COMPLEX*16 XTHING,XVE,XVF,XVEF
937 DATA xi/(0.d0,1.d0)/,xr/(1.d0,0.d0)/
940 DATA svar0,cost0 /-5.d0,-6.d0/
941 DATA pi /3.141592653589793238462643d0/
942 DATA seps1,seps2 /0d0,0d0/
945 IF ( mode.NE.mode0.OR.svar.NE.svar0.OR.costhe.NE.cost0
946 $ .OR.ide0.NE.ide)
THEN 954 sinthe=sqrt(1.d0-costhe**2)
955 beta=sqrt(max(0d0,1d0-4d0*amfin**2/svar))
957 xupzfp(1)=0.5d0*(xupzf(1)+xupzf(2))+0.5*beta*(xupzf(1)-xupzf(2))
958 xupzfp(2)=0.5d0*(xupzf(1)+xupzf(2))-0.5*beta*(xupzf(1)-xupzf(2))
959 xupzip(1)=0.5d0*(xupzi(1)+xupzi(2))+0.5*(xupzi(1)-xupzi(2))
960 xupzip(2)=0.5d0*(xupzi(1)+xupzi(2))-0.5*(xupzi(1)-xupzi(2))
962 xupf =0.5d0*(xupzf(1)+xupzf(2))
963 xupi =0.5d0*(xupzi(1)+xupzi(2))
967 propz =1d0/dcmplx(svar-amz**2,svar/amz*gammz)
968 IF (keygsw.EQ.0) propz=0.d0
971 regula= (3-2*i)*(3-2*j) + costhe
972 regulm=-(3-2*i)*(3-2*j) * sinthe *2.d0*amfin/sqrt(svar)
973 aphot(i,j)=propa*(xupgi(i)*xupgf(j)*regula)
974 azett(i,j)=propz*(xupzip(i)*xupzfp(j)+xthing)*regula
975 aborn(i,j)=aphot(i,j)+azett(i,j)
976 aphotm(i,j)=propa*dcmplx(0d0,1d0)*xupgi(i)*xupgf(j)*regulm
977 azettm(i,j)=propz*dcmplx(0d0,1d0)*(xupzip(i)*xupf+xthing)*regulm
978 abornm(i,j)=aphotm(i,j)+azettm(i,j)
993 factor=kolor*(1d0+helic*polar1)*(1d0-helic*polar2)/4d0
994 factom=factor*(1+helit*ta)*(1-helit*tb)
995 factor=factor*(1+helit*ta)*(1+helit*tb)
997 born=born+cdabs(aborn(i,j))**2*factor
1000 born=born+cdabs(abornm(i,j))**2*factom
1006 IF(funt.LT.0.d0) funt=born
1009 IF (svar.GT.4d0*amfin**2)
THEN 1011 thresh=sqrt(1-4d0*amfin**2/svar)
1012 t_born= funt*svar**2*thresh
1022 SUBROUTINE initwk(IDEX,IDFX,SVAR)
1024 IMPLICIT REAL*8 (a-h,o-z)
1025 COMMON / t_beampm / ene ,amin,amfin,ide,idf
1026 REAL*8 ENE ,AMIN,AMFIN
1027 COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
1028 & ,xupgi ,xupzi ,xupgf ,xupzf
1029 & ,ndiag0,ndiaga,keya,keyz
1030 & ,itce,jtce,itcf,jtcf,kolor
1031 REAL*8 SS,POLN,T3E,QE,T3F,QF
1032 & ,XUPGI(2),XUPZI(2),XUPGF(2),XUPZF(2)
1033 COMMON / t_gswprm /swsq,amw,amz,amh,amtop,gammz
1034 REAL*8 SWSQ,AMW,AMZ,AMH,AMTOP,GAMMZ
1046 IF (idfx.EQ. 15)
then 1049 ELSEIF (idfx.EQ.-15)
then 1053 WRITE(*,*)
'INITWK: WRONG IDFX' 1057 IF (idex.EQ. 11)
then 1060 ELSEIF (idex.EQ.-11)
then 1063 ELSEIF (idex.EQ. 13)
then 1066 ELSEIF (idex.EQ.-13)
then 1069 ELSEIF (idex.EQ. 1)
then 1072 ELSEIF (idex.EQ.- 1)
then 1075 ELSEIF (idex.EQ. 2)
then 1078 ELSEIF (idex.EQ.- 2)
then 1081 ELSEIF (idex.EQ. 3)
then 1084 ELSEIF (idex.EQ.- 3)
then 1087 ELSEIF (idex.EQ. 4)
then 1090 ELSEIF (idex.EQ.- 4)
then 1093 ELSEIF (idex.EQ. 5)
then 1096 ELSEIF (idex.EQ.- 5)
then 1099 ELSEIF (idex.EQ. 12)
then 1102 ELSEIF (idex.EQ.- 12)
then 1105 ELSEIF (idex.EQ. 14)
then 1108 ELSEIF (idex.EQ.- 14)
then 1111 ELSEIF (idex.EQ. 16)
then 1114 ELSEIF (idex.EQ.- 16)
then 1119 WRITE(*,*)
'INITWK: WRONG IDEX' 1133 CALL t_givizo( ide, 1,aizor,qe,kdumm)
1134 CALL t_givizo( ide,-1,aizol,qe,kdumm)
1138 xupzi(1)=(aizor-qe*swsq)/sqrt(swsq*(1-swsq))
1139 xupzi(2)=(aizol-qe*swsq)/sqrt(swsq*(1-swsq))
1140 CALL t_givizo( idf, 1,aizor,qf,kolor)
1141 CALL t_givizo( idf,-1,aizol,qf,kolor)
1145 xupzf(1)=(aizor-qf*swsq)/sqrt(swsq*(1-swsq))
1146 xupzf(2)=(aizol-qf*swsq)/sqrt(swsq*(1-swsq))
1157 SUBROUTINE t_givizo(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR)
1169 IMPLICIT REAL*8(a-h,o-z)
1171 IF(idferm.EQ.0.OR.iabs(idferm).GT.4)
GOTO 901
1172 IF(iabs(ihelic).NE.1)
GOTO 901
1174 idtype =iabs(idferm)
1176 lepqua=int(idtype*0.4999999d0)
1177 iupdow=idtype-2*lepqua-1
1178 charge =(-iupdow+2d0/3d0*lepqua)*ic
1179 sizo3 =0.25d0*(ic-ih)*(1-2*iupdow)
1184 901 print *,
' STOP IN GIVIZO: WRONG PARAMS.' 1187 SUBROUTINE phyfix(NSTOP,NSTART)
1188 common/lujets/n,k(4000,5),p(4000,5),v(4000,5)
1194 IF(k(i,1).NE.21)
THEN 1202 SUBROUTINE filhep(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
1215 parameter(nmxhep=4000)
1217 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1245 ELSE IF (n.GT.0)
THEN 1256 IF ((ihep.LE.0).OR.(ihep.GT.nmxhep))
RETURN 1263 IF(jmo1.LT.0)jmohep(1,ihep)=jmohep(1,ihep)+ihep
1265 IF(jmo2.LT.0)jmohep(2,ihep)=jmohep(2,ihep)+ihep
1280 DO ip=jmohep(1,ihep),jmohep(2,ihep)
1284 IF(isthep(ip).EQ.1)isthep(ip)=2
1287 IF(jdahep(1,ip).EQ.0)
THEN 1291 jdahep(2,ip)=max(ihep,jdahep(2,ip))
1300 FUNCTION ihepdim(DUM)
1304 parameter(nmxhep=4000)
1306 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1328 IMPLICIT REAL*8(a-h,o-z)
1329 COMPLEX*16 CPRZ0,CPRZ0M
1332 cprz0=dcmplx((s-amz**2),s/amz*gammz)
1334 zprop2=(abs(cprz0m))**2
1337 SUBROUTINE taupi0(MODE,JAK,ION)
1351 parameter(nmxhep=4000)
1353 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1374 COMMON /taupos/ np1,np2
1376 REAL PHOT1(4),PHOT2(4)
1377 REAL*8 R,X(4),Y(4),PI0(4)
1378 INTEGER JEZELI(3),ION(3)
1381 IF (mode.EQ.-1)
THEN 1387 IF (jezeli(1).EQ.0)
RETURN 1388 IF (jezeli(2).EQ.1)
CALL taueta(jak)
1389 IF (jezeli(3).EQ.1)
CALL tauk0s(jak)
1391 IF((kto.EQ. 1).OR.(kto.EQ.11))
THEN 1397 DO k=jdahep(1,nps),nhepm
1398 IF (idhep(k).EQ.111.AND.jdahep(1,k).LE.k)
THEN 1403 r=sqrt(pi0(4)**2-pi0(3)**2-pi0(2)**2-pi0(1)**2)/2d0
1412 CALL bostdq(-1,pi0,x,x)
1413 CALL bostdq(-1,pi0,y,y)
1419 CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1420 CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1425 SUBROUTINE taueta(JAK)
1433 parameter(nmxhep=4000)
1435 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1454 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1455 * ,ampiz,ampi,amro,gamro,ama1,gama1
1456 * ,amk,amkz,amkst,gamkst
1458 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1459 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1460 * ,AMK,AMKZ,AMKST,GAMKST
1463 COMMON /taupos/ np1,np2
1465 REAL RRR(1),BRSUM(3), RR(2)
1466 REAL PHOT1(4),PHOT2(4),PHOT3(4)
1467 REAL*8 X(4), Y(4), Z(4)
1469 REAL*8 R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM
1471 xlam(a,b,c)=sqrt(abs((a-b-c)**2-4.0*b*c))
1473 IF((kto.EQ. 1).OR.(kto.EQ.11))
THEN 1479 DO k=jdahep(1,nps),nhepm
1480 IF (idhep(k).EQ.221.AND.jdahep(1,k).LE.k)
THEN 1486 brsum(2)=brsum(1)+0.319
1487 brsum(3)=brsum(2)+0.237
1490 IF (rrr(1).LT.brsum(1))
THEN 1492 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1501 CALL bostdq(-1,peta,x,x)
1502 CALL bostdq(-1,peta,y,y)
1508 CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1509 CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1511 IF(rrr(1).LT.brsum(2))
THEN 1518 ELSEIF(rrr(1).LT.brsum(3))
THEN 1535 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)
1538 am2=sqrt(amin**2+rr(1)*(amax**2-amin**2))
1540 wt=xlam(r**2,am2**2,xm3**2)*xlam(am2**2,xm1**2,xm2**2)
1542 IF (rr(2).GT.wt)
GOTO 7
1544 ru=xlam(am2**2,xm1**2,xm2**2)/am2/2
1548 x(4)=sqrt(ru**2+xm1**2)
1549 y(4)=sqrt(ru**2+xm2**2)
1555 ru=xlam(r**2,am2**2,xm3**2)/r/2
1557 z(4)=sqrt(ru**2+am2**2)
1559 CALL bostdq(-1,z,x,x)
1560 CALL bostdq(-1,z,y,y)
1565 z(4)=sqrt(ru**2+xm3**2)
1567 CALL bostdq(-1,peta,x,x)
1568 CALL bostdq(-1,peta,y,y)
1569 CALL bostdq(-1,peta,z,z)
1579 CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1580 CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1581 CALL filhep(0,1,id3,k,k,0,0,phot3,ym3,.true.)
1588 SUBROUTINE tauk0s(JAK)
1596 parameter(nmxhep=4000)
1598 INTEGER nevhep,nhep,isthep,idhep,jmohep,
1618 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1619 * ,ampiz,ampi,amro,gamro,ama1,gama1
1620 * ,amk,amkz,amkst,gamkst
1622 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
1623 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
1624 * ,AMK,AMKZ,AMKST,GAMKST
1627 COMMON /taupos/ np1,np2
1629 REAL RRR(1),BRSUM(3), RR(2)
1630 REAL PHOT1(4),PHOT2(4),PHOT3(4)
1631 REAL*8 X(4), Y(4), Z(4)
1633 REAL*8 R,RU,PETA(4),XM1,XM2,XM3,XM,XLAM
1635 xlam(a,b,c)=sqrt(abs((a-b-c)**2-4.0*b*c))
1637 IF((kto.EQ. 1).OR.(kto.EQ.11))
THEN 1643 DO k=jdahep(1,nps),nhepm
1644 IF (idhep(k).EQ.310.AND.jdahep(1,k).LE.k)
THEN 1653 brsum(3)=brsum(2)+0.237
1656 IF(rrr(1).LT.brsum(1))
THEN 1661 ELSEIF(rrr(1).LT.brsum(2))
THEN 1674 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1676 r=sqrt(abs(r**2-xm1**2))
1685 CALL bostdq(-1,peta,x,x)
1686 CALL bostdq(-1,peta,y,y)
1695 CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1696 CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)