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.)