1 SUBROUTINE tauola(MODE,KEYPOL)
9 #include "../../include/HEPEVT.h"
10 COMMON /taupos/ np1, np2
11 REAL*4 phoi(4),phof(4)
12 double precision q1(4),q2(4),p1(4),p2(4),p3(4),p4(4)
13 COMMON / momdec / q1,q2,p1,p2,p3,p4
15 COMMON /libra/ jak1,jak2,itdkrc,ifphot,ifhadm,ifhadp
18 common /pseudocoup/ csc,ssc
21 COMMON / inout / inut,iout
24 dimension pol1(4), pol2(4)
25 double precision pol1x(4), pol2x(4)
27 DATA pol1 /0.0,0.0,0.0,0.0/
28 DATA pol2 /0.0,0.0,0.0,0.0/
29 DATA pi /3.141592653589793238462643d0/
37 COMMON /isons_tau/ison(2)
63 betah=sqrt(1d0-4*xmtau**2/xmh**2)
67 IF (ifphot.EQ.1) CALL phoini
68 CALL inietc(jak1,jak2,itdkrc,ifphot)
78 WRITE(iout,7001) pol,psi,ion(1),ion(2),ion(3)
79 ELSEIF(mode.EQ.0)
THEN
84 CALL phyfix(nstop,nstart)
97 IF(abs(idhep(i)).EQ.kftau.AND.isthep(i).EQ.1.AND.
98 $ (isthep(i).GE.125.OR.isthep(i).LT.120))
THEN
100 DO WHILE (abs(idhep(imoth)).EQ.kftau)
101 imoth=jmohep(1,imoth)
103 IF (isthep(imoth).EQ.3.OR.
104 $ (isthep(imoth).GE.120.AND.isthep(imoth).LE.125))
THEN
106 IF (idhep(j).EQ.idhep(imoth).AND.
107 $ jmohep(1,j).EQ.imoth.AND.
108 $ isthep(j).EQ.2)
THEN
118 IF(jmoth.EQ.imother(ii)) goto 9999
138 IF (idhep(jmohep(1,im0)).EQ.idhep(im0)) im0=jmohep(1,im0)
141 IF (isthep(i).EQ.3.OR.
142 $ (isthep(i).GE.120.AND.isthep(i).LE.125))
THEN
146 DO WHILE (idhep(imoth).EQ.idhep(i).OR.abs(idhep(imoth)).EQ.kftau)
147 imoth=jmohep(1,imoth)
149 IF ((imoth.EQ.im0.OR.imoth.EQ.im).AND.isel.EQ.-1)
THEN
152 ELSEIF ((imoth.EQ.im0.OR.imoth.EQ.im).AND.isel.EQ.0)
THEN
154 ELSEIF ((imoth.NE.im0.AND.imoth.NE.im).AND.isel.EQ.0)
THEN
179 IF(idhep(i).EQ.-kftau.OR.idhep(i).EQ.-kfnue) ncount=ncount+1
180 IF(idhep(i).EQ. kftau.OR.idhep(i).EQ. kfnue) ncount=ncount+1
187 DO i=max(np1+1,ison(1)),ison(2)
189 IF(idhep(i).EQ.-kftau.OR.idhep(i).EQ.-kfnue) np1=i
192 DO i=max(np2+1,ison(1)),ison(2)
194 IF(idhep(i).EQ. kftau.OR.idhep(i).EQ. kfnue) np2=i
209 IF(idhep(im).EQ.kfhiggs(1).OR.idhep(im).EQ.kfhiggs(2).OR.
210 $ idhep(im).EQ.kfhiggs(3))
THEN
211 IF(rrr(1).LT.0.5)
THEN
218 ELSEIF((idhep(im).EQ.kfz0).OR.(idhep(im).EQ.kfgam))
THEN
226 polz0=plzapx(.true.,im,np1,np2)
227 IF(rrr(1).LT.polz0)
THEN
234 ELSEIF(idhep(np1).EQ.-idhep(np2))
THEN
235 polz0=plzapx(.true.,im,np1,np2)
236 IF(rrr(1).LT.polz0)
THEN
243 ELSEIF(abs(idhep(im)).EQ.kfhigch)
THEN
253 IF(idhep(im).EQ.kfhiggs(1).OR.idhep(im).EQ.kfhiggs(2).OR.
254 $ idhep(im).EQ.kfhiggs(3))
THEN
255 IF(idhep(np1).EQ.-kftau .AND.
256 $ (jdahep(1,np1).LE.np1.OR.jdahep(1,np1).GT.nhep) .AND.
257 $ idhep(np2).EQ. kftau .AND.
258 $ (jdahep(1,np2).LE.np2.OR.jdahep(1,np2).GT.nhep)
260 IF (idhep(im).EQ.kfhiggs(1))
THEN
262 ELSEIF (idhep(im).EQ.kfhiggs(2))
THEN
264 ELSEIF (idhep(im).EQ.kfhiggs(3))
THEN
267 WRITE(*,*)
'Warning from TAUOLA:'
268 WRITE(*,*)
'I stop this run, wrong IDHEP(IM)=',idhep(im)
271 CALL spinhiggs(im,np1,np2,ifpseudo,pol1,pol2)
272 IF (ifphot.EQ.1) CALL photos(im)
278 IF(idhep(np1).EQ.-kftau.AND.
279 $ (jdahep(1,np1).LE.np1.OR.jdahep(1,np1).GT.nhep))
THEN
282 IF (ifphot.EQ.1) CALL photos(np1)
286 IF(idhep(np2).EQ. kftau.AND.
287 $ (jdahep(1,np2).LE.np2.OR.jdahep(1,np2).GT.nhep))
THEN
290 IF (ifphot.EQ.1) CALL photos(np2)
295 IF (ncount.GT.0) goto 666
298 ELSEIF(mode.EQ.1)
THEN
302 CALL dekay(100,pol1x)
306 7001
FORMAT(///1x,15(5h*****)
307 $ /,
' *', 25x,
'*****TAUOLA UNIVERSAL INTERFACE: ******',9x,1h*,
308 $ /,
' *', 25x,
'*****VERSION 1.22, April 2009 (gfort)**',9x,1h*,
309 $ /,
' *', 25x,
'**AUTHORS: P. Golonka, B. Kersevan, ***',9x,1h*,
310 $ /,
' *', 25x,
'**T. Pierzchala, E. Richter-Was, ******',9x,1h*,
311 $ /,
' *', 25x,
'****** Z. Was, M. Worek ***************',9x,1h*,
312 $ /,
' *', 25x,
'**USEFUL DISCUSSIONS, IN PARTICULAR ***',9x,1h*,
313 $ /,
' *', 25x,
'*WITH C. Biscarat and S. Slabospitzky**',9x,1h*,
314 $ /,
' *', 25x,
'****** are warmly acknowledged ********',9x,1h*,
315 $ /,
' *', 25x,
' ',9x,1h*,
316 $ /,
' *', 25x,
'********** INITIALIZATION ************',9x,1h*,
317 $ /,
' *',f20.5,5x,
'tau polarization switch must be 1 or 0 ',9x,1h*,
318 $ /,
' *',f20.5,5x,
'Higs scalar/pseudo mix CERN-TH/2003-166',9x,1h*,
319 $ /,
' *',i10, 15x,
'PI0 decay switch must be 1 or 0 ',9x,1h*,
320 $ /,
' *',i10, 15x,
'ETA decay switch must be 1 or 0 ',9x,1h*,
321 $ /,
' *',i10, 15x,
'K0S decay switch must be 1 or 0 ',9x,1h*,
324 7002
FORMAT(///1x,15(5h*****)
325 $ /,
' *', 25x,
'*****TAUOLA UNIVERSAL INTERFACE: ******',9x,1h*,
326 $ /,
' *', 25x,
'*****VERSION 1.22, April 2009 (gfort)**',9x,1h*,
327 $ /,
' *', 25x,
'**AUTHORS: P. Golonka, B. Kersevan, ***',9x,1h*,
328 $ /,
' *', 25x,
'**T. Pierzchala, E. Richter-Was, ******',9x,1h*,
329 $ /,
' *', 25x,
'****** Z. Was, M. Worek ***************',9x,1h*,
330 $ /,
' *', 25x,
'**USEFUL DISCUSSIONS, IN PARTICULAR ***',9x,1h*,
331 $ /,
' *', 25x,
'*WITH C. Biscarat and S. Slabospitzky**',9x,1h*,
332 $ /,
' *', 25x,
'****** are warmly acknowledged ********',9x,1h*,
333 $ /,
' *', 25x,
'****** END OF MODULE OPERATION ********',9x,1h*,
338 SUBROUTINE spinhiggs(IM,NP1,NP2,IFPSEUDO,Pol1,Pol2)
340 REAL*8 hh1,hh2,wthiggs
341 dimension pol1(4), pol2(4),hh1(4),hh2(4), rrr(1)
349 wt=wthiggs(ifpseudo,hh1,hh2)
350 IF (rrr(1).GT.wt) goto 10
356 FUNCTION wthiggs(IFPSEUDO,HH1,HH2)
358 common /pseudocoup/ csc,ssc
361 REAL*8 hh1(4),hh2(4),r(4,4),wthiggs
375 r(1,1)=(csc**2-ssc**2)/(csc**2+ssc**2)
376 r(2,2)=(csc**2-ssc**2)/(csc**2+ssc**2)
377 r(1,2)=2*csc*ssc/(csc**2+ssc**2)
378 r(2,1)=-2*csc*ssc/(csc**2+ssc**2)
388 wthiggs=wthiggs+r(k,l)*hh1(k)*hh2(l)
394 FUNCTION plzapx(HOPEin,IM0,NP1,NP2)
399 REAL*8 plzap0,svar,costhe,sini,sfin,zprop2,
400 $ p1(4),p2(4),q1(4),q2(4),qq(4),ph(4),pd1(4),pd2(4),
401 $ pq1(4),pq2(4),pb(4),pa(4)
405 #include "../../include/HEPEVT.h"
409 COMMON /isons_tau/ison(2)
420 IF (idhep(im0).EQ.idhep(im00)) im=jmohep(1,im0)
429 IF (isthep(im00).EQ.120)
THEN
437 IF (imo1.EQ.0.AND.imo2.EQ.0)
THEN
440 IF (idhep(imo1).EQ.idhep(np1)) imo1=jmohep(1,imo1)
444 IF (idhep(imo2).EQ.idhep(np2)) imo2=jmohep(1,imo2)
449 ELSEIF (idhep(im).NE.22.AND.idhep(im).NE.23)
THEN
452 IF (idhep(imo1).EQ.idhep(np1)) imo1=jmohep(1,imo1)
456 IF (idhep(imo2).EQ.idhep(np2)) imo2=jmohep(1,imo2)
463 IF (imo1.EQ.0) hope=.false.
464 IF (imo2.EQ.0) hope=.false.
465 IF (imo1.EQ.imo2) hope=.false.
475 IF (im.EQ.jmohep(1,im0).AND.
476 $ (idhep(im).EQ.22.OR.idhep(im).EQ.23))
THEN
483 CALL bostdq( 1,pa, q1, q1)
484 CALL bostdq( 1,pa, q2, q2)
485 CALL bostdq(-1,pb, q1, q1)
486 CALL bostdq(-1,pb, q2, q2)
492 IF (hope) p1(i)=phep(i,imo1)
493 IF (hope) p2(i)=phep(i,imo2)
501 IF (hope) idfp1=idhep(imo1)
502 IF (hope) idfp2=idhep(imo2)
504 svar=qq(4)**2-qq(3)**2-qq(2)**2-qq(1)**2
524 IF (iffull.EQ.1) inbr=np1
525 IF (idhep(jmohep(1,inbr)).EQ.idhep(inbr)) inbr=jmohep(1,inbr)
527 IF(nx1.EQ.0.OR.nx2.EQ.0)
THEN
531 IF(jmohep(1,inbr-k).EQ.jmohep(1,inbr))
THEN
540 IF(jmohep(1,k).EQ.jmohep(1,inbr))
THEN
550 IF (abs(idfp1).GE.20.AND.abs(idfp2).GE.20) hope=.false.
552 IF (abs(idfp1).LE.20.AND.abs(idfp2).LE.20.AND.idfp1+idfp2.NE.0)
562 IF (abs(idfp1).LT.20) ide= idfp1
563 IF (abs(idfp2).LT.20) ide=-idfp2
583 iflav=min(abs(idfp1),abs(idfp2))
593 IF (abs(idfp1).GE.20)
THEN
596 IF (abs(idp).EQ.iflav)
THEN
604 IF (abs(idfp2).GE.20)
THEN
607 IF (abs(idp).EQ.iflav)
THEN
617 IF (abs(idfp1).GE.20)
THEN
621 xm1=abs((pd1(4)+ph(4))**2-(pd1(3)+ph(3))**2
622 $ -(pd1(2)+ph(2))**2-(pd1(1)+ph(1))**2)
623 xm2=abs((pd2(4)+ph(4))**2-(pd2(3)+ph(3))**2
624 $ -(pd2(2)+ph(2))**2-(pd2(1)+ph(1))**2)
639 IF (abs(idfp2).GE.20)
THEN
643 xm1=abs((pd1(4)+ph(4))**2-(pd1(3)+ph(3))**2
644 $ -(pd1(2)+ph(2))**2-(pd1(1)+ph(1))**2)
645 xm2=abs((pd2(4)+ph(4))**2-(pd2(3)+ph(3))**2
646 $ -(pd2(2)+ph(2))**2-(pd2(1)+ph(1))**2)
663 IF (idhep(jmohep(1,np1)).EQ.idhep(np1)) nph1=jmohep(1,np1)
664 IF (idhep(jmohep(1,np2)).EQ.idhep(np2)) nph2=jmohep(1,np2)
668 IF (abs(idhep(k)).NE.iflav.AND.k.NE.im.AND.
670 $ k.NE.nph1.AND.k.NE.nph2)
THEN
672 IF(idhep(k).EQ.22.AND.iffull.EQ.1)
THEN
676 xm1=abs((pd1(4)-ph(4))**2-(pd1(3)-ph(3))**2
677 $ -(pd1(2)-ph(2))**2-(pd1(1)-ph(1))**2)
678 xm2=abs((pd2(4)-ph(4))**2-(pd2(3)-ph(3))**2
679 $ -(pd2(2)-ph(2))**2-(pd2(1)-ph(1))**2)
680 xm3=abs((pq1(4)+ph(4))**2-(pq1(3)+ph(3))**2
681 $ -(pq1(2)+ph(2))**2-(pq1(1)+ph(1))**2)
682 xm4=abs((pq2(4)+ph(4))**2-(pq2(3)+ph(3))**2
683 $ -(pq2(2)+ph(2))**2-(pq2(1)+ph(1))**2)
686 sini=abs((pd1(4)+pd2(4)-ph(4))**2-(pd1(3)+pd2(3)-ph(3))**2
687 $ -(pd1(2)+pd2(2)-ph(2))**2-(pd1(1)+pd2(1)-ph(1))**2)
688 sfin=abs((pd1(4)+pd2(4) )**2-(pd1(3)+pd2(3) )**2
689 $ -(pd1(2)+pd2(2) )**2-(pd1(1)+pd2(1) )**2)
699 xm=min(xm1,xm2,xm3,xm4)
704 ELSEIF (xm2.EQ.xm)
THEN
708 ELSEIF (xm3.EQ.xm)
THEN
721 xm1=abs((pd1(4)-ph(4))**2-(pd1(3)-ph(3))**2
722 $ -(pd1(2)-ph(2))**2-(pd1(1)-ph(1))**2)
723 xm2=abs((pd2(4)-ph(4))**2-(pd2(3)-ph(3))**2
724 $ -(pd2(2)-ph(2))**2-(pd2(1)-ph(1))**2)
754 IF (abs(idhep(im0)).EQ.22.OR.abs(idhep(im0)).EQ.23)
THEN
756 IF(abs(idhep(k)).EQ.22)
THEN
763 xm3=abs((pq1(4)+ph(4))**2-(pq1(3)+ph(3))**2
764 $ -(pq1(2)+ph(2))**2-(pq1(1)+ph(1))**2)
765 xm4=abs((pq2(4)+ph(4))**2-(pq2(3)+ph(3))**2
766 $ -(pq2(2)+ph(2))**2-(pq2(1)+ph(1))**2)
789 CALL angulu(pd1,pd2,q1,q2,costhe)
791 plzapx=plzap0(ide,idfq1,svar,costhe)
794 SUBROUTINE angulu(PD1,PD2,Q1,Q2,COSTHE)
795 REAL*8 pd1(4),pd2(4),q1(4),q2(4),costhe,p(4),qq(4),qt(4)
801 xm1=abs(pd1(4)**2-pd1(3)**2-pd1(2)**2-pd1(1)**2)
802 xm2=abs(pd2(4)**2-pd2(3)**2-pd2(2)**2-pd2(1)**2)
820 xmqq=sqrt(qq(4)**2-qq(3)**2-qq(2)**2-qq(1)**2)
822 qtxqq=qt(4)*qq(4)-qt(3)*qq(3)-qt(2)*qq(2)-qt(1)*qq(1)
824 qt(k)=qt(k)-qq(k)*qtxqq/xmqq**2
827 pxqq=p(4)*qq(4)-p(3)*qq(3)-p(2)*qq(2)-p(1)*qq(1)
829 p(k)=p(k)-qq(k)*pxqq/xmqq**2
832 pxp =sqrt(p(1)**2+p(2)**2+p(3)**2-p(4)**2)
833 qtxqt=sqrt(qt(3)**2+qt(2)**2+qt(1)**2-qt(4)**2)
834 pxqt =p(3)*qt(3)+p(2)*qt(2)+p(1)*qt(1)-p(4)*qt(4)
835 costhe=pxqt/pxp/qtxqt
839 FUNCTION plzap0(IDE,IDF,SVAR,COSTH0)
842 REAL*8 plzap0,svar,costhe,costh0,t_born
849 CALL initwk(ide,idf,svar)
851 CALL initwk(-ide,-idf,svar)
853 plzap0=t_born(0,svar,costhe,1d0,1d0)
854 $ /(t_born(0,svar,costhe,1d0,1d0)+t_born(0,svar,costhe,-1d0,-1d0))
858 FUNCTION t_born(MODE,SVAR,COSTHE,TA,TB)
869 IMPLICIT REAL*8(a-h,o-z)
870 COMMON / t_beampm / ene ,amin,amfin,ide,idf
871 REAL*8 ene ,amin,amfin
872 COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
873 & ,xupgi ,xupzi ,xupgf ,xupzf
874 & ,ndiag0,ndiaga,keya,keyz
875 & ,itce,jtce,itcf,jtcf,kolor
876 REAL*8 ss,poln,t3e,qe,t3f,qf
877 & ,xupgi(2),xupzi(2),xupgf(2),xupzf(2)
880 COMMON / t_gswprm /swsq,amw,amz,amh,amtop,gammz
881 REAL*8 swsq,amw,amz,amh,amtop,gammz
887 COMPLEX*16 aborn(2,2),aphot(2,2),azett(2,2)
888 COMPLEX*16 xupzfp(2),xupzip(2)
889 COMPLEX*16 abornm(2,2),aphotm(2,2),azettm(2,2)
890 COMPLEX*16 propa,propz
892 COMPLEX*16 xupf,xupi,xff(4),xfem,xfota,xrho,xke,xkf,xkef
893 COMPLEX*16 xthing,xve,xvf,xvef
894 DATA xi/(0.d0,1.d0)/,xr/(1.d0,0.d0)/
897 DATA svar0,cost0 /-5.d0,-6.d0/
898 DATA pi /3.141592653589793238462643d0/
899 DATA seps1,seps2 /0d0,0d0/
902 IF ( mode.NE.mode0.OR.svar.NE.svar0.OR.costhe.NE.cost0
903 $ .OR.ide0.NE.ide)
THEN
911 sinthe=sqrt(1.d0-costhe**2)
912 beta=sqrt(max(0d0,1d0-4d0*amfin**2/svar))
914 xupzfp(1)=0.5d0*(xupzf(1)+xupzf(2))+0.5*beta*(xupzf(1)-xupzf(2))
915 xupzfp(2)=0.5d0*(xupzf(1)+xupzf(2))-0.5*beta*(xupzf(1)-xupzf(2))
916 xupzip(1)=0.5d0*(xupzi(1)+xupzi(2))+0.5*(xupzi(1)-xupzi(2))
917 xupzip(2)=0.5d0*(xupzi(1)+xupzi(2))-0.5*(xupzi(1)-xupzi(2))
919 xupf =0.5d0*(xupzf(1)+xupzf(2))
920 xupi =0.5d0*(xupzi(1)+xupzi(2))
924 propz =1d0/dcmplx(svar-amz**2,svar/amz*gammz)
925 IF (keygsw.EQ.0) propz=0.d0
928 regula= (3-2*i)*(3-2*j) + costhe
929 regulm=-(3-2*i)*(3-2*j) * sinthe *2.d0*amfin/sqrt(svar)
930 aphot(i,j)=propa*(xupgi(i)*xupgf(j)*regula)
931 azett(i,j)=propz*(xupzip(i)*xupzfp(j)+xthing)*regula
932 aborn(i,j)=aphot(i,j)+azett(i,j)
933 aphotm(i,j)=propa*dcmplx(0d0,1d0)*xupgi(i)*xupgf(j)*regulm
934 azettm(i,j)=propz*dcmplx(0d0,1d0)*(xupzip(i)*xupf+xthing)*regulm
935 abornm(i,j)=aphotm(i,j)+azettm(i,j)
950 factor=kolor*(1d0+helic*polar1)*(1d0-helic*polar2)/4d0
951 factom=factor*(1+helit*ta)*(1-helit*tb)
952 factor=factor*(1+helit*ta)*(1+helit*tb)
954 born=born+cdabs(aborn(i,j))**2*factor
957 born=born+cdabs(abornm(i,j))**2*factom
963 IF(funt.LT.0.d0) funt=born
966 IF (svar.GT.4d0*amfin**2)
THEN
968 thresh=sqrt(1-4d0*amfin**2/svar)
969 t_born= funt*svar**2*thresh
979 SUBROUTINE initwk(IDEX,IDFX,SVAR)
981 IMPLICIT REAL*8 (a-h,o-z)
982 COMMON / t_beampm / ene ,amin,amfin,ide,idf
983 REAL*8 ene ,amin,amfin
984 COMMON / t_gauspm /ss,poln,t3e,qe,t3f,qf
985 & ,xupgi ,xupzi ,xupgf ,xupzf
986 & ,ndiag0,ndiaga,keya,keyz
987 & ,itce,jtce,itcf,jtcf,kolor
988 REAL*8 ss,poln,t3e,qe,t3f,qf
989 & ,xupgi(2),xupzi(2),xupgf(2),xupzf(2)
990 COMMON / t_gswprm /swsq,amw,amz,amh,amtop,gammz
991 REAL*8 swsq,amw,amz,amh,amtop,gammz
1003 IF (idfx.EQ. 15)
then
1006 ELSEIF (idfx.EQ.-15)
then
1010 WRITE(*,*)
'INITWK: WRONG IDFX'
1014 IF (idex.EQ. 11)
then
1017 ELSEIF (idex.EQ.-11)
then
1020 ELSEIF (idex.EQ. 13)
then
1023 ELSEIF (idex.EQ.-13)
then
1026 ELSEIF (idex.EQ. 1)
then
1029 ELSEIF (idex.EQ.- 1)
then
1032 ELSEIF (idex.EQ. 2)
then
1035 ELSEIF (idex.EQ.- 2)
then
1038 ELSEIF (idex.EQ. 3)
then
1041 ELSEIF (idex.EQ.- 3)
then
1044 ELSEIF (idex.EQ. 4)
then
1047 ELSEIF (idex.EQ.- 4)
then
1050 ELSEIF (idex.EQ. 5)
then
1053 ELSEIF (idex.EQ.- 5)
then
1056 ELSEIF (idex.EQ. 12)
then
1059 ELSEIF (idex.EQ.- 12)
then
1062 ELSEIF (idex.EQ. 14)
then
1065 ELSEIF (idex.EQ.- 14)
then
1068 ELSEIF (idex.EQ. 16)
then
1071 ELSEIF (idex.EQ.- 16)
then
1076 WRITE(*,*)
'INITWK: WRONG IDEX'
1090 CALL t_givizo( ide, 1,aizor,qe,kdumm)
1091 CALL t_givizo( ide,-1,aizol,qe,kdumm)
1095 xupzi(1)=(aizor-qe*swsq)/sqrt(swsq*(1-swsq))
1096 xupzi(2)=(aizol-qe*swsq)/sqrt(swsq*(1-swsq))
1097 CALL t_givizo( idf, 1,aizor,qf,kolor)
1098 CALL t_givizo( idf,-1,aizol,qf,kolor)
1102 xupzf(1)=(aizor-qf*swsq)/sqrt(swsq*(1-swsq))
1103 xupzf(2)=(aizol-qf*swsq)/sqrt(swsq*(1-swsq))
1114 SUBROUTINE t_givizo(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR)
1126 IMPLICIT REAL*8(a-h,o-z)
1128 IF(idferm.EQ.0.OR.iabs(idferm).GT.4) goto 901
1129 IF(iabs(ihelic).NE.1) goto 901
1131 idtype =iabs(idferm)
1133 lepqua=int(idtype*0.4999999d0)
1134 iupdow=idtype-2*lepqua-1
1135 charge =(-iupdow+2d0/3d0*lepqua)*ic
1136 sizo3 =0.25d0*(ic-ih)*(1-2*iupdow)
1141 901 print *,
' STOP IN GIVIZO: WRONG PARAMS.'
1144 #include "../../include/phyfix.h"
1145 SUBROUTINE filhep(N,IST,ID,JMO1,JMO2,JDA1,JDA2,P4,PINV,PHFLAG)
1156 #include "../../include/HEPEVT.h"
1166 ELSE IF (n.GT.0)
THEN
1177 IF ((ihep.LE.0).OR.(ihep.GT.nmxhep))
RETURN
1184 IF(jmo1.LT.0)jmohep(1,ihep)=jmohep(1,ihep)+ihep
1186 IF(jmo2.LT.0)jmohep(2,ihep)=jmohep(2,ihep)+ihep
1201 DO ip=jmohep(1,ihep),jmohep(2,ihep)
1205 IF(isthep(ip).EQ.1)isthep(ip)=2
1208 IF(jdahep(1,ip).EQ.0)
THEN
1212 jdahep(2,ip)=max(ihep,jdahep(2,ip))
1221 FUNCTION ihepdim(DUM)
1223 #include "../../include/HEPEVT.h"
1227 IMPLICIT REAL*8(a-h,o-z)
1228 COMPLEX*16 cprz0,cprz0m
1231 cprz0=dcmplx((s-amz**2),s/amz*gammz)
1233 zprop2=(abs(cprz0m))**2
1236 SUBROUTINE taupi0(MODE,JAK,ION)
1248 #include "../../include/HEPEVT.h"
1251 COMMON /taupos/ np1,np2
1253 REAL phot1(4),phot2(4)
1254 REAL*8 r,x(4),y(4),pi0(4)
1255 INTEGER jezeli(3),ion(3)
1258 IF (mode.EQ.-1)
THEN
1264 IF (jezeli(1).EQ.0)
RETURN
1265 IF (jezeli(2).EQ.1) CALL taueta(jak)
1266 IF (jezeli(3).EQ.1) CALL tauk0s(jak)
1268 IF((kto.EQ. 1).OR.(kto.EQ.11))
THEN
1274 DO k=jdahep(1,nps),nhepm
1275 IF (idhep(k).EQ.111.AND.jdahep(1,k).LE.k)
THEN
1280 r=sqrt(pi0(4)**2-pi0(3)**2-pi0(2)**2-pi0(1)**2)/2d0
1289 CALL bostdq(-1,pi0,x,x)
1290 CALL bostdq(-1,pi0,y,y)
1296 CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1297 CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1302 SUBROUTINE taueta(JAK)
1308 #include "../../include/HEPEVT.h"
1309 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1310 * ,ampiz,ampi,amro,gamro,ama1,gama1
1311 * ,amk,amkz,amkst,gamkst
1313 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1314 * ,ampiz,ampi,amro,gamro,ama1,gama1
1315 * ,amk,amkz,amkst,gamkst
1318 COMMON /taupos/ np1,np2
1320 REAL rrr(1),brsum(3), rr(2)
1321 REAL phot1(4),phot2(4),phot3(4)
1322 REAL*8 x(4), y(4), z(4)
1324 REAL*8 r,ru,peta(4),xm1,xm2,xm3,xm,xlam
1326 xlam(a,b,c)=sqrt(abs((a-b-c)**2-4.0*b*c))
1328 IF((kto.EQ. 1).OR.(kto.EQ.11))
THEN
1334 DO k=jdahep(1,nps),nhepm
1335 IF (idhep(k).EQ.221.AND.jdahep(1,k).LE.k)
THEN
1341 brsum(2)=brsum(1)+0.319
1342 brsum(3)=brsum(2)+0.237
1345 IF (rrr(1).LT.brsum(1))
THEN
1347 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1356 CALL bostdq(-1,peta,x,x)
1357 CALL bostdq(-1,peta,y,y)
1363 CALL filhep(0,1,22,k,k,0,0,phot1,0.0,.true.)
1364 CALL filhep(0,1,22,k,k,0,0,phot2,0.0,.true.)
1366 IF(rrr(1).LT.brsum(2))
THEN
1373 ELSEIF(rrr(1).LT.brsum(3))
THEN
1390 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)
1393 am2=sqrt(amin**2+rr(1)*(amax**2-amin**2))
1395 wt=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)
1396 & *xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)
1398 IF (rr(2).GT.wt) goto 7
1400 ru=xlam(1d0*am2**2,1d0*xm1**2,1d0*xm2**2)/am2/2
1404 x(4)=sqrt(ru**2+xm1**2)
1405 y(4)=sqrt(ru**2+xm2**2)
1411 ru=xlam(1d0*r**2,1d0*am2**2,1d0*xm3**2)/r/2
1413 z(4)=sqrt(ru**2+am2**2)
1415 CALL bostdq(-1,z,x,x)
1416 CALL bostdq(-1,z,y,y)
1421 z(4)=sqrt(ru**2+xm3**2)
1423 CALL bostdq(-1,peta,x,x)
1424 CALL bostdq(-1,peta,y,y)
1425 CALL bostdq(-1,peta,z,z)
1435 CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1436 CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)
1437 CALL filhep(0,1,id3,k,k,0,0,phot3,ym3,.true.)
1444 SUBROUTINE tauk0s(JAK)
1450 #include "../../include/HEPEVT.h"
1452 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
1453 * ,ampiz,ampi,amro,gamro,ama1,gama1
1454 * ,amk,amkz,amkst,gamkst
1456 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
1457 * ,ampiz,ampi,amro,gamro,ama1,gama1
1458 * ,amk,amkz,amkst,gamkst
1461 COMMON /taupos/ np1,np2
1463 REAL rrr(1),brsum(3), rr(2)
1464 REAL phot1(4),phot2(4),phot3(4)
1465 REAL*8 x(4), y(4), z(4)
1467 REAL*8 r,ru,peta(4),xm1,xm2,xm3,xm,xlam
1469 xlam(a,b,c)=sqrt(abs((a-b-c)**2-4.0*b*c))
1471 IF((kto.EQ. 1).OR.(kto.EQ.11))
THEN
1477 DO k=jdahep(1,nps),nhepm
1478 IF (idhep(k).EQ.310.AND.jdahep(1,k).LE.k)
THEN
1487 brsum(3)=brsum(2)+0.237
1490 IF(rrr(1).LT.brsum(1))
THEN
1495 ELSEIF(rrr(1).LT.brsum(2))
THEN
1508 r=sqrt(peta(4)**2-peta(3)**2-peta(2)**2-peta(1)**2)/2d0
1510 r=sqrt(abs(r**2-xm1**2))
1519 CALL bostdq(-1,peta,x,x)
1520 CALL bostdq(-1,peta,y,y)
1529 CALL filhep(0,1,id1,k,k,0,0,phot1,ym1,.true.)
1530 CALL filhep(0,1,id2,k,k,0,0,phot2,ym2,.true.)