1 FUNCTION formom(XMAA,XMOM)
6 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
7 * ,ampiz,ampi,amro,gamro,ama1,gama1
8 * ,amk,amkz,amkst,gamkst
10 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
11 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
12 * ,AMK,AMKZ,AMKST,GAMKST
13 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
14 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
19 bwign(xm,am,gamma)=1./cmplx(xm**2-am**2,gamma*am)
30 fqed =sqrt(4.0*3.1415926535/137.03604)
31 formom=fqed*fro**2/sqrt(2.0)*gcoup**2*bwign(xmom,amom,gamom)
32 $ *(bwign(xmaa,amro,gamro)+elpha*bwign(xmaa,amrop,gamrop))
33 $ *(bwign( 0.0,amro,gamro)+elpha*bwign( 0.0,amrop,gamrop))
39 COMPLEX FUNCTION fk1ab(XMSQ,INDX)
44 COMPLEX F1,F2,AMPA,AMPB
64 ampa = cmplx(pkorb(3,81),0.)
65 ampb = cmplx(pkorb(3,82),0.)
66 ELSE IF (indx.EQ.2)
THEN 67 ampa = cmplx(pkorb(3,83),0.)
68 ampb = cmplx(pkorb(3,84),0.)
69 ELSEIF (indx.EQ.3)
THEN 70 ampa = cmplx(pkorb(3,85),0.)
71 ampb = cmplx(pkorb(3,86),0.)
72 ELSEIF (indx.EQ.4)
THEN 73 ampa = cmplx(pkorb(3,87),0.)
74 ampb = cmplx(pkorb(3,88),0.)
80 f1 = cmplx(-xm1sq,0.0)/cmplx(xmsq-xm1sq,fg1)
81 f2 = cmplx(-xm2sq,0.0)/cmplx(xmsq-xm2sq,fg2)
82 fk1ab = ampa*f1+ampb*f2
87 FUNCTION form1(MNUM,QQ,S1,SDWA)
97 COMPLEX FORM1,WIGNER,WIGFOR,FPIKM,BWIGM
98 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
99 * ,ampiz,ampi,amro,gamro,ama1,gama1
100 * ,amk,amkz,amkst,gamkst
102 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
103 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
104 * ,AMK,AMKZ,AMKST,GAMKST
106 COMPLEX FORMA1,FORMK1,FORMRO,FORMKS
107 COMPLEX FA1A1P,FK1AB,F3PI
115 form1 = f3pi(1,qq,s1,sdwa)
117 ELSEIF (mnum.EQ.1)
THEN 119 formks = bwigm(s1,amkst,gamkst,ampi,amk)
121 form1 = forma1*formks
123 ELSEIF (mnum.EQ.2)
THEN 125 formks = bwigm(s1,amkst,gamkst,ampi,amk)
127 form1 = forma1*formks
129 ELSEIF (mnum.EQ.3)
THEN 131 formks = bwigm(s1,amkst,gamkst,ampi,amk)
133 form1 = forma1*formks
135 ELSEIF (mnum.EQ.4)
THEN 137 formks = bwigm(s1,amkst,gamkst,ampi,amk)
139 form1 = formk1*formks
141 ELSEIF (mnum.EQ.5)
THEN 144 formro = fpikm(sqrt(s1),ampi,ampi)
145 form1 = formk1*formro
147 ELSEIF (mnum.EQ.6)
THEN 150 formks = bwigm(s1,amkst,gamkst,amk,ampi)
151 form1 = formk1*formks
153 ELSEIF (mnum.EQ.7)
THEN 158 wigner(a,b,c)= cmplx(1.0,0.0)/cmplx(a-b**2,b*c)
161 gamax=gama1*gfun(qq)/gfun(ama1**2)
162 form1=ama1**2*wigner(qq,ama1,gamax)*fpikm(sqrt(s1),ampi,ampi)
163 ELSEIF (mnum.EQ.1)
THEN 166 form1=bwigm(s1,amkst,gamkst,ampi,amkz)
168 form1=bwigm(s1,amkst,gamkst,ampi,amk)
170 gamax=gama1*gfun(qq)/gfun(ama1**2)
171 form1=ama1**2*wigner(qq,ama1,gamax)*form1
172 ELSEIF (mnum.EQ.2)
THEN 175 form1=bwigm(s1,amkst,gamkst,ampi,amkz)
177 form1=bwigm(s1,amkst,gamkst,ampi,amk)
179 gamax=gama1*gfun(qq)/gfun(ama1**2)
180 form1=ama1**2*wigner(qq,ama1,gamax)*form1
181 ELSEIF (mnum.EQ.3)
THEN 184 gamax=gama1*gfun(qq)/gfun(ama1**2)
185 form1=ama1**2*wigner(qq,ama1,gamax)*form1
186 ELSEIF (mnum.EQ.4)
THEN 191 form1=bwigm(s1,amkst,gamkst,amk,ampiz)
193 form1=bwigm(s1,amkst,gamkst,amk,ampi)
195 form1=wigfor(qq,xm2,gam2)*form1
196 ELSEIF (mnum.EQ.5)
THEN 200 form1=wigfor(qq,xm2,gam2)*fpikm(sqrt(s1),ampi,ampi)
201 ELSEIF (mnum.EQ.6)
THEN 203 ELSEIF (mnum.EQ.7)
THEN 209 FUNCTION form2(MNUM,QQ,S1,SDWA)
219 COMPLEX FORM2,WIGNER,WIGFOR,FPIKM,BWIGM
220 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
221 * ,ampiz,ampi,amro,gamro,ama1,gama1
222 * ,amk,amkz,amkst,gamkst
224 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
225 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
226 * ,AMK,AMKZ,AMKST,GAMKST
228 COMPLEX FORMA1,FORMK1,FORMRO,FORMKS
229 COMPLEX FA1A1P,FK1AB,F3PI
237 form2 = f3pi(2,qq,s1,sdwa)
239 ELSEIF (mnum.EQ.1)
THEN 241 formro = fpikm(sqrt(s1),amk,amk)
243 form2 = forma1*formro
245 ELSEIF (mnum.EQ.2)
THEN 247 formro = fpikm(sqrt(s1),amk,amk)
249 form2 = forma1*formro
251 ELSEIF (mnum.EQ.3)
THEN 253 formro = fpikm(sqrt(s1),amk,amk)
255 form2 = forma1*formro
257 ELSEIF (mnum.EQ.4)
THEN 259 formks = bwigm(s1,amkst,gamkst,ampi,amk)
261 form2 = formk1*formks
263 ELSEIF (mnum.EQ.5)
THEN 265 formks = bwigm(s1,amkst,gamkst,ampi,amk)
267 form2 = formk1*formks
269 ELSEIF (mnum.EQ.6)
THEN 271 formro = fpikm(sqrt(s1),ampi,ampi)
273 form2 = formk1*formro
275 ELSEIF (mnum.EQ.7)
THEN 281 wigner(a,b,c)= cmplx(1.0,0.0)/cmplx(a-b**2,b*c)
284 gamax=gama1*gfun(qq)/gfun(ama1**2)
285 form2=ama1**2*wigner(qq,ama1,gamax)*fpikm(sqrt(s1),ampi,ampi)
286 ELSEIF (mnum.EQ.1)
THEN 288 gamax=gama1*gfun(qq)/gfun(ama1**2)
289 form2=ama1**2*wigner(qq,ama1,gamax)*fpikm(sqrt(s1),ampi,ampi)
290 ELSEIF (mnum.EQ.2)
THEN 292 gamax=gama1*gfun(qq)/gfun(ama1**2)
293 form2=ama1**2*wigner(qq,ama1,gamax)*fpikm(sqrt(s1),ampi,ampi)
294 ELSEIF (mnum.EQ.3)
THEN 296 gamax=gama1*gfun(qq)/gfun(ama1**2)
297 form2=ama1**2*wigner(qq,ama1,gamax)*fpikm(sqrt(s1),ampi,ampi)
298 ELSEIF (mnum.EQ.4)
THEN 303 form2=bwigm(s1,amkst,gamkst,amk,ampiz)
305 form2=bwigm(s1,amkst,gamkst,amk,ampi)
307 form2=wigfor(qq,xm2,gam2)*form2
308 ELSEIF (mnum.EQ.5)
THEN 313 form2=bwigm(s1,amkst,gamkst,amk,ampiz)
315 form2=bwigm(s1,amkst,gamkst,amk,ampi)
317 form2=wigfor(qq,xm2,gam2)*form2
319 ELSEIF (mnum.EQ.6)
THEN 322 form2=wigfor(qq,xm2,gam2)*fpikm(sqrt(s1),ampi,ampi)
324 ELSEIF (mnum.EQ.7)
THEN 331 COMPLEX FUNCTION bwigm(S,M,G,XM1,XM2)
344 IF (s.GT.(xm1+xm2)**2)
THEN 345 qs=sqrt(abs((s -(xm1+xm2)**2)*(s -(xm1-xm2)**2)))/sqrt(s)
346 qm=sqrt(abs((m**2-(xm1+xm2)**2)*(m**2-(xm1-xm2)**2)))/m
348 gs=g*(m/w)**2*(qs/qm)**3
352 bwigm=m**2/cmplx(m**2-s,-sqrt(s)*gs)
355 COMPLEX FUNCTION fpikm(W,XM1,XM2)
360 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
377 fpikm=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
381 COMPLEX FUNCTION fpikmd(W,XM1,XM2)
386 REAL ROM,ROG,ROM1,ROG1,PI,PIM,S,W
406 fpikmd=(delta*bwigm(s,rom,rog,xm1,xm2)
407 $ +beta*bwigm(s,rom1,rog1,xm1,xm2)
408 $ + bwigm(s,rom2,rog2,xm1,xm2))
413 FUNCTION form3(MNUM,QQ,S1,SDWA)
423 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
424 * ,ampiz,ampi,amro,gamro,ama1,gama1
425 * ,amk,amkz,amkst,gamkst
427 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
428 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
429 * ,AMK,AMKZ,AMKST,GAMKST
432 COMPLEX FORMA1,FORMK1,FORMRO,FORMKS
433 COMPLEX FA1A1P,FK1AB,F3PI
441 form3 = f3pi(3,qq,s1,sdwa)
443 ELSEIF (mnum.EQ.3)
THEN 445 formks = bwigm(s1,amkst,gamkst,ampiz,amk)
447 form3 = forma1*formks
449 ELSEIF (mnum.EQ.6)
THEN 451 formks = bwigm(s1,amkst,gamkst,amk,ampi)
453 form3 = formk1*formks
469 FUNCTION form4(MNUM,QQ,S1,S2,S3)
477 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
478 * ,ampiz,ampi,amro,gamro,ama1,gama1
479 * ,amk,amkz,amkst,gamkst
481 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
482 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
483 * ,AMK,AMKZ,AMKST,GAMKST
484 COMPLEX FORM4,WIGNER,FPIKM
488 wigner(a,b,c)=cmplx(1.0,0.0) /cmplx(a-b**2,b*c)
501 IF (s.GT.(xm1+xm2)**2)
THEN 502 qs=sqrt(abs((s -(xm1+xm2)**2)*(s -(xm1-xm2)**2)))/sqrt(s)
503 qm=sqrt(abs((m**2-(xm1+xm2)**2)*(m**2-(xm1-xm2)**2)))/m
505 gs=g*(m/w)**2*(qs/qm)**5
510 form4=g1*g2*fpip/amro**4/ampip**2
511 $ *ampip**2*wigner(qq,ampip,gamx)
512 $ *( s1*(s2-s3)*fpikm(sqrt(s1),ampiz,ampiz)
513 $ +s2*(s1-s3)*fpikm(sqrt(s2),ampiz,ampiz) )
514 ELSEIF (mnum.EQ.1)
THEN 526 IF (s.GT.(xm1+xm2)**2)
THEN 527 qs=sqrt(abs((s -(xm1+xm2)**2)*(s -(xm1-xm2)**2)))/sqrt(s)
528 qm=sqrt(abs((m**2-(xm1+xm2)**2)*(m**2-(xm1-xm2)**2)))/m
530 gs=g*(m/w)**2*(qs/qm)**5
535 form4=g1*g2*fpip/amro**4/ampip**2
536 $ *ampip**2*wigner(qq,ampip,gamx)
537 $ *( s1*(s2-s3)*fpikm(sqrt(s1),ampiz,ampiz)
538 $ +s2*(s1-s3)*fpikm(sqrt(s2),ampiz,ampiz) )
551 FUNCTION form5(MNUM,QQ,S1,S2)
560 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
561 * ,ampiz,ampi,amro,gamro,ama1,gama1
562 * ,amk,amkz,amkst,gamkst
564 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
565 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
566 * ,AMK,AMKZ,AMKST,GAMKST
567 COMPLEX FORM5,WIGNER,FPIKM,FPIKMD,BWIGM
570 wigner(a,b,c)=cmplx(1.0,0.0)/cmplx(a-b**2,b*c)
571 #elif defined (ALEPH) 572 wigner(a,b,c)=cmplx(1.0,0.0)/cmplx(a-b**2,b*c)
577 ELSEIF (mnum.EQ.1)
THEN 580 form5=fpikmd(sqrt(qq),ampi,ampi)/(1+elpha)
581 $ *( fpikm(sqrt(s2),ampi,ampi)
582 $ +elpha*bwigm(s1,amkst,gamkst,ampi,amk))
583 ELSEIF (mnum.EQ.2)
THEN 586 form5=fpikmd(sqrt(qq),ampi,ampi)/(1+elpha)
587 $ *( fpikm(sqrt(s2),ampi,ampi)
588 $ +elpha*bwigm(s1,amkst,gamkst,ampi,amk))
589 ELSEIF (mnum.EQ.3)
THEN 592 ELSEIF (mnum.EQ.4)
THEN 595 ELSEIF (mnum.EQ.5)
THEN 598 form5=bwigm(qq,amkst,gamkst,ampi,amk)/(1+elpha)
599 $ *( fpikm(sqrt(s1),ampi,ampi)
600 $ +elpha*bwigm(s2,amkst,gamkst,ampi,amk))
601 ELSEIF (mnum.EQ.6)
THEN 604 form5=bwigm(qq,amkst,gamkst,ampi,amkz)/(1+elpha)
605 $ *( fpikm(sqrt(s2),ampi,ampi)
606 $ +elpha*bwigm(s1,amkst,gamkst,ampi,amk))
607 ELSEIF (mnum.EQ.7)
THEN 609 form5=fpikmd(sqrt(qq),ampi,ampi)*fpikm(sqrt(s1),ampi,ampi)
614 SUBROUTINE curr(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
616 SUBROUTINE currx(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
625 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
626 * ,ampiz,ampi,amro,gamro,ama1,gama1
627 * ,amk,amkz,amkst,gamkst
629 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
630 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
631 * ,AMK,AMKZ,AMKST,GAMKST
632 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
633 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
637 COMMON /arbit/ arflat,aromeg
639 REAL PIM1(4),PIM2(4),PIM3(4),PIM4(4),PAA(4)
642 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,WIGFOR
644 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FPIKM
649 DATA pi /3.141592653589793238462643/
651 bwign(a,xm,xg)=1.0/cmplx(a-xm**2,xm*xg)
686 coef1=2.0*sqrt(3.0)/fpi**2*arflat
693 paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)
703 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
706 sk=(pp(k,4)+pim4(4))**2-(pp(k,3)+pim4(3))**2
707 $ -(pp(k,2)+pim4(2))**2-(pp(k,1)+pim4(1))**2
717 denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
718 $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
724 $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
732 form1=wigfor(sk,amro,gamro)
734 form1= fpikm(sqrt(sk),ampi,ampi) *fpikm(sqrt(qq),ampi,ampi)
744 $ hadcur(i)+cmplx(fix*coef1)*form1*aa(i,j)*(pp(k,j)-pp(4,j))
770 IF (k.EQ.4) sign= 1.0
771 qqa=qqa+sign*(paa(k)-pa(k))**2
772 ss23=ss23+sign*(pb(k) +pim3(k))**2
773 ss24=ss24+sign*(pb(k) +pim4(k))**2
774 ss34=ss34+sign*(pim3(k)+pim4(k))**2
775 qp1p2=qp1p2+sign*(paa(k)-pa(k))*pb(k)
776 qp1p3=qp1p3+sign*(paa(k)-pa(k))*pim3(k)
777 qp1p4=qp1p4+sign*(paa(k)-pa(k))*pim4(k)
778 p1p2=p1p2+sign*pa(k)*pb(k)
779 p1p3=p1p3+sign*pa(k)*pim3(k)
780 p1p4=p1p4+sign*pa(k)*pim4(k)
783 form2=coef2*(bwign(qq,amro,gamro)+elpha*bwign(qq,amrop,gamrop))
786 form3=bwign(qqa,amom,gamom)
789 hadcur(k)=hadcur(k)+form2*form3*(
790 $ pb(k)*(qp1p3*p1p4-qp1p4*p1p3)
791 $ +pim3(k)*(qp1p4*p1p2-qp1p2*p1p4)
792 $ +pim4(k)*(qp1p2*p1p3-qp1p3*p1p2) )
800 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
803 sk=(pp(k,4)+pim4(4))**2-(pp(k,3)+pim4(3))**2
804 $ -(pp(k,2)+pim4(2))**2-(pp(k,1)+pim4(1))**2
815 denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
816 $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
822 $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
830 form1=wigfor(sk,amro,gamro)
832 form1= fpikm(sqrt(sk),ampi,ampi) *fpikm(sqrt(qq),ampi,ampi)
839 $ hadcur(i)+cmplx(coef1)*form1*aa(i,j)*(pp(k,j)-pp(4,j))
845 FUNCTION wigfor(S,XM,XGAM)
846 COMPLEX WIGFOR,WIGNOR
847 wignor=cmplx(-xm**2,xm*xgam)
848 wigfor=wignor/cmplx(s-xm**2,xm*xgam)
855 COMMON /inout/ inut, iout
856 WRITE (unit = iout,fmt = 99)
857 WRITE (unit = iout,fmt = 98)
860 . /,
' *************************************************** ',
861 . /,
' YOU ARE USING THE 4 PION DECAY MODE FORM FACTORS ',
862 . /,
' WHICH HAVE BEEN DESCRIBED IN:',
863 . /,
' R. DECKER, M. FINKEMEIER, P. HEILIGER AND H.H. JONSSON',
864 . /,
' "TAU DECAYS INTO FOUR PIONS" ',
865 . /,
' UNIVERSITAET KARLSRUHE PREPRINT TTP 94-13 (1994);',
866 . /,
' LNF-94/066(IR); HEP-PH/9410260 ',
868 . /,
' PLEASE NOTE THAT THIS ROUTINE IS USING PARAMETERS',
869 . /,
' RELATED TO THE 3 PION DECAY MODE (A1 MODE), SUCH AS',
870 . /,
' THE A1 MASS AND WIDTH (TAKEN FROM THE COMMON /PARMAS/)',
871 . /,
' AND THE 2 PION VECTOR RESONANCE FORM FACTOR (BY USING',
872 . /,
' THE ROUTINE FPIKM)' ,
873 . /,
' THUS IF YOU DECIDE TO CHANGE ANY OF THESE, YOU WILL' ,
874 . /,
' HAVE TO REFIT THE 4 PION PARAMETERS IN THE COMMON' )
876 .
' BLOCK /TAU4PI/, OR YOU MIGHT GET A BAD DISCRIPTION' ,
877 . /,
' OF TAU -> 4 PIONS' ,
878 . /,
' for these formfactors set in routine CHOICE for',
879 . /, .eq.
' mnum102 -- AMRX=1.42 and GAMRX=.21',
880 . /, .eq.
' mnum101 -- AMRX=1.3 and GAMRX=.46 PROB1,PROB2=0.2',
881 . /,
' to optimize phase space parametrization',
882 . /,
' *************************************************** ',
883 . /,
' coded by M. Finkemeier and P. Heiliger, 29. sept. 1994',
884 . /,
' incorporated to TAUOLA by Z. Was 17. jan. 1995',
887 . /,
' changed by: Z. Was on 17.01.95',
888 . /,
' changes by: M. Finkemeier on 30.01.95' )
892 COMMON /tau4pi/ gomega,gamma1,gamma2,rom1,rog1,beta1,
894 REAL*4 GOMEGA,GAMMA1,GAMMA2,ROM1,ROG1,BETA1,
906 COMPLEX FUNCTION bwiga1(QA)
911 COMMON / parmas/ amtau,amnuta,amel,amnue,ammu,amnumu,
912 % AMPIZ,ampi,amro,gamro,ama1,gama1,
913 % AMK,amkz,amkst,gamkst
916 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU,
917 % AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1,
918 % AMK,AMKZ,AMKST,GAMKST
920 wigner(a,b,c)=cmplx(1.0,0.0)/cmplx(a-b**2,b*c)
921 gamax=gama1*gfun(qa)/gfun(ama1**2)
922 bwiga1=-ama1**2*wigner(qa,ama1,gamax)
925 COMPLEX FUNCTION bwigeps(QEPS)
932 bwigeps=cmplx(meps**2,-meps*geps)/
933 % CMPLX(meps**2-qeps,-meps*geps)
936 COMPLEX FUNCTION frho4(W,XM1,XM2)
942 COMMON /tau4pi/ gomega,gamma1,gamma2,rom1,rog1,beta1,
944 REAL*4 GOMEGA,GAMMA1,GAMMA2,ROM1,ROG1,BETA1,
946 REAL ROM,ROG,PI,PIM,S,W
961 frho4=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2)
962 & +beta2*bwigm(s,rom2,rog2,xm1,xm2))
966 SUBROUTINE curr(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
977 COMMON /tau4pi/ gomega,gamma1,gamma2,rom1,rog1,beta1,
979 REAL*4 GOMEGA,GAMMA1,GAMMA2,ROM1,ROG1,BETA1,
981 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
982 * ,ampiz,ampi,amro,gamro,ama1,gama1
983 * ,amk,amkz,amkst,gamkst
985 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
986 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
987 * ,AMK,AMKZ,AMKST,GAMKST
988 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
989 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
990 REAL PIM1(4),PIM2(4),PIM3(4),PIM4(4),PAA(4)
991 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FPIKM
993 COMPLEX BWIGEPS,BWIGA1
994 COMPLEX HCOMP1(4),HCOMP2(4),HCOMP3(4),HCOMP4(4)
996 COMPLEX T243,T213,T143,T123,T341,T342
997 COMPLEX T124,T134,T214,T234,T314,T324
998 COMPLEX S2413,S2314,S1423,S1324,S34
1000 COMPLEX BRACK1,BRACK2,BRACK3,BRACK4A,BRACK4B,BRACK4
1002 REAL QMP1,QMP2,QMP3,QMP4
1003 REAL PS43,PS41,PS42,PS34,PS14,PS13,PS24,PS23
1006 REAL PD243,PD241,PD213,PD143,PD142
1007 REAL PD123,PD341,PD342,PD413,PD423
1008 REAL PD124,PD134,PD214,PD234,PD314,PD324
1009 REAL QP1,QP2,QP3,QP4
1012 REAL AA(4,4),PP(4,4)
1013 DATA pi /3.141592653589793238462643/
1016 bwign(a,xm,xg)=1.0/cmplx(a-xm**2,xm*xg)
1037 coef1=2.0*sqrt(3.0)/fpi**2*arflat
1043 hadcur(k)=cmplx(0.0)
1044 paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)
1054 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
1060 qmp1=(pim2(4)+pim3(4)+pim4(4))**2-(pim2(3)+pim3(3)+pim4(3))**2
1061 % -(pim2(2)+pim3(2)+pim4(2))**2-(pim2(1)+pim3(1)+pim4(1))**2
1063 qmp2=(pim1(4)+pim3(4)+pim4(4))**2-(pim1(3)+pim3(3)+pim4(3))**2
1064 % -(pim1(2)+pim3(2)+pim4(2))**2-(pim1(1)+pim3(1)+pim4(1))**2
1066 qmp3=(pim1(4)+pim2(4)+pim4(4))**2-(pim1(3)+pim2(3)+pim4(3))**2
1067 % -(pim1(2)+pim2(2)+pim4(2))**2-(pim1(1)+pim2(1)+pim4(1))**2
1069 qmp4=(pim1(4)+pim2(4)+pim3(4))**2-(pim1(3)+pim2(3)+pim3(3))**2
1070 % -(pim1(2)+pim2(2)+pim3(2))**2-(pim1(1)+pim2(1)+pim3(1))**2
1075 ps43=(pim4(4)+pim3(4))**2-(pim4(3)+pim3(3))**2
1076 % -(pim4(2)+pim3(2))**2-(pim4(1)+pim3(1))**2
1078 ps41=(pim4(4)+pim1(4))**2-(pim4(3)+pim1(3))**2
1079 % -(pim4(2)+pim1(2))**2-(pim4(1)+pim1(1))**2
1081 ps42=(pim4(4)+pim2(4))**2-(pim4(3)+pim2(3))**2
1082 % -(pim4(2)+pim2(2))**2-(pim4(1)+pim2(1))**2
1088 ps13=(pim1(4)+pim3(4))**2-(pim1(3)+pim3(3))**2
1089 % -(pim1(2)+pim3(2))**2-(pim1(1)+pim3(1))**2
1093 ps23=(pim2(4)+pim3(4))**2-(pim2(3)+pim3(3))**2
1094 % -(pim2(2)+pim3(2))**2-(pim2(1)+pim3(1))**2
1096 pd243=pim2(4)*(pim4(4)-pim3(4))-pim2(3)*(pim4(3)-pim3(3))
1097 % -pim2(2)*(pim4(2)-pim3(2))-pim2(1)*(pim4(1)-pim3(1))
1099 pd241=pim2(4)*(pim4(4)-pim1(4))-pim2(3)*(pim4(3)-pim1(3))
1100 % -pim2(2)*(pim4(2)-pim1(2))-pim2(1)*(pim4(1)-pim1(1))
1102 pd213=pim2(4)*(pim1(4)-pim3(4))-pim2(3)*(pim1(3)-pim3(3))
1103 % -pim2(2)*(pim1(2)-pim3(2))-pim2(1)*(pim1(1)-pim3(1))
1105 pd143=pim1(4)*(pim4(4)-pim3(4))-pim1(3)*(pim4(3)-pim3(3))
1106 % -pim1(2)*(pim4(2)-pim3(2))-pim1(1)*(pim4(1)-pim3(1))
1108 pd142=pim1(4)*(pim4(4)-pim2(4))-pim1(3)*(pim4(3)-pim2(3))
1109 % -pim1(2)*(pim4(2)-pim2(2))-pim1(1)*(pim4(1)-pim2(1))
1111 pd123=pim1(4)*(pim2(4)-pim3(4))-pim1(3)*(pim2(3)-pim3(3))
1112 % -pim1(2)*(pim2(2)-pim3(2))-pim1(1)*(pim2(1)-pim3(1))
1114 pd341=pim3(4)*(pim4(4)-pim1(4))-pim3(3)*(pim4(3)-pim1(3))
1115 % -pim3(2)*(pim4(2)-pim1(2))-pim3(1)*(pim4(1)-pim1(1))
1117 pd342=pim3(4)*(pim4(4)-pim2(4))-pim3(3)*(pim4(3)-pim2(3))
1118 % -pim3(2)*(pim4(2)-pim2(2))-pim3(1)*(pim4(1)-pim2(1))
1120 pd413=pim4(4)*(pim1(4)-pim3(4))-pim4(3)*(pim1(3)-pim3(3))
1121 % -pim4(2)*(pim1(2)-pim3(2))-pim4(1)*(pim1(1)-pim3(1))
1123 pd423=pim4(4)*(pim2(4)-pim3(4))-pim4(3)*(pim2(3)-pim3(3))
1124 % -pim4(2)*(pim2(2)-pim3(2))-pim4(1)*(pim2(1)-pim3(1))
1128 qp1=pim1(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1129 % -pim1(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1130 % -pim1(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1131 % -pim1(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1133 qp2=pim2(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1134 % -pim2(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1135 % -pim2(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1136 % -pim2(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1138 qp3=pim3(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1139 % -pim3(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1140 % -pim3(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1141 % -pim3(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1143 qp4=pim4(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1144 % -pim4(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1145 % -pim4(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1146 % -pim4(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1152 t243=bwiga1(qmp2)*fpikm(sqrt(ps43),ampi,ampi)*gamma1
1153 t213=bwiga1(qmp2)*fpikm(sqrt(ps13),ampi,ampi)*gamma1
1154 t143=bwiga1(qmp1)*fpikm(sqrt(ps43),ampi,ampi)*gamma1
1155 t123=bwiga1(qmp1)*fpikm(sqrt(ps23),ampi,ampi)*gamma1
1156 t341=bwiga1(qmp3)*fpikm(sqrt(ps41),ampi,ampi)*gamma1
1157 t342=bwiga1(qmp3)*fpikm(sqrt(ps42),ampi,ampi)*gamma1
1161 s2413=frho4(sqrt(ps24),ampi,ampi)*gamma2
1162 s2314=frho4(sqrt(ps23),ampi,ampi)*bwigeps(ps14)*gamma2
1163 s1423=frho4(sqrt(ps14),ampi,ampi)*gamma2
1164 s1324=frho4(sqrt(ps13),ampi,ampi)*bwigeps(ps24)*gamma2
1165 s34=frho4(sqrt(ps34),ampi,ampi)*gamma2
1169 brack1=2.*t143+2.*t243+t123+t213
1170 % +t341*(pd241/qmp3-1.)+t342*(pd142/qmp3-1.)
1171 % +3./4.*(s1423+s2413-s2314-s1324)-3.*s34
1173 brack2=2.*t143*pd243/qmp1+3.*t213
1174 % +t123*(2.*pd423/qmp1+1.)+t341*(pd241/qmp3+3.)
1175 % +t342*(pd142/qmp3+1.)
1176 % -3./4.*(s2314+3.*s1324+3.*s1423+s2413)
1178 brack3=2.*t243*pd143/qmp2+3.*t123
1179 % +t213*(2.*pd413/qmp2+1.)+t341*(pd241/qmp3+1.)
1180 % +t342*(pd142/qmp3+3.)
1181 % -3./4.*(3.*s2314+s1324+s1423+3.*s2413)
1183 brack4a=2.*t143*(pd243/qq*(qp1/qmp1+1.)+pd143/qq)
1184 % +2.*t243*(pd143/qq*(qp2/qmp2+1.)+pd243/qq)
1186 % +2.*t123*(pd423/qq*(qp1/qmp1+1.)+pd123/qq)
1187 % +2.*t213*(pd413/qq*(qp2/qmp2+1.)+pd213/qq)
1188 % +t341*(pd241/qmp3+1.-2.*pd241/qq*(qp3/qmp3+1.)
1190 % +t342*(pd142/qmp3+1.-2.*pd142/qq*(qp3/qmp3+1.)
1193 brack4b=-3./4.*(s2314*(2.*(qp2-qp3)/qq+1.)
1194 % +s1324*(2.*(qp1-qp3)/qq+1.)
1195 % +s1423*(2.*(qp1-qp4)/qq+1.)
1196 % +s2413*(2.*(qp2-qp4)/qq+1.)
1197 % +4.*s34*(qp4-qp3)/qq)
1199 brack4=brack4a+brack4b
1203 hcomp1(k)=(pim3(k)-pim4(k))*brack1
1204 hcomp2(k)=pim1(k)*brack2
1205 hcomp3(k)=pim2(k)*brack3
1206 hcomp4(k)=(pim1(k)+pim2(k)+pim3(k)+pim4(k))*brack4
1212 hadcur(i)=hcomp1(i)-hcomp2(i)-hcomp3(i)+hcomp4(i)
1213 hadcur(i)=-coef1*frho4(sqrt(qq),ampi,ampi)*hadcur(i)
1242 IF (k.EQ.4) sign= 1.0
1243 qqa=qqa+sign*(paa(k)-pa(k))**2
1244 ss23=ss23+sign*(pb(k) +pim3(k))**2
1245 ss24=ss24+sign*(pb(k) +pim4(k))**2
1246 ss34=ss34+sign*(pim3(k)+pim4(k))**2
1247 qp1p2=qp1p2+sign*(paa(k)-pa(k))*pb(k)
1248 qp1p3=qp1p3+sign*(paa(k)-pa(k))*pim3(k)
1249 qp1p4=qp1p4+sign*(paa(k)-pa(k))*pim4(k)
1250 p1p2=p1p2+sign*pa(k)*pb(k)
1251 p1p3=p1p3+sign*pa(k)*pim3(k)
1252 p1p4=p1p4+sign*pa(k)*pim4(k)
1255 form2=coef2*(bwign(qq,amro,gamro)+elpha*bwign(qq,amrop,gamrop))
1258 form3=bwign(qqa,amom,gamom)
1261 hadcur(k)=hadcur(k)+form2*form3*(
1262 $ pb(k)*(qp1p3*p1p4-qp1p4*p1p3)
1263 $ +pim3(k)*(qp1p4*p1p2-qp1p2*p1p4)
1264 $ +pim4(k)*(qp1p2*p1p3-qp1p3*p1p2) )
1272 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
1279 qmp1=(pim2(4)+pim3(4)+pim4(4))**2-(pim2(3)+pim3(3)+pim4(3))**2
1280 % -(pim2(2)+pim3(2)+pim4(2))**2-(pim2(1)+pim3(1)+pim4(1))**2
1282 qmp2=(pim1(4)+pim3(4)+pim4(4))**2-(pim1(3)+pim3(3)+pim4(3))**2
1283 % -(pim1(2)+pim3(2)+pim4(2))**2-(pim1(1)+pim3(1)+pim4(1))**2
1285 qmp3=(pim1(4)+pim2(4)+pim4(4))**2-(pim1(3)+pim2(3)+pim4(3))**2
1286 % -(pim1(2)+pim2(2)+pim4(2))**2-(pim1(1)+pim2(1)+pim4(1))**2
1288 qmp4=(pim1(4)+pim2(4)+pim3(4))**2-(pim1(3)+pim2(3)+pim3(3))**2
1289 % -(pim1(2)+pim2(2)+pim3(2))**2-(pim1(1)+pim2(1)+pim3(1))**2
1294 ps14=(pim1(4)+pim4(4))**2-(pim1(3)+pim4(3))**2
1295 % -(pim1(2)+pim4(2))**2-(pim1(1)+pim4(1))**2
1297 ps21=(pim2(4)+pim1(4))**2-(pim2(3)+pim1(3))**2
1298 % -(pim2(2)+pim1(2))**2-(pim2(1)+pim1(1))**2
1300 ps23=(pim2(4)+pim3(4))**2-(pim2(3)+pim3(3))**2
1301 % -(pim2(2)+pim3(2))**2-(pim2(1)+pim3(1))**2
1303 ps24=(pim2(4)+pim4(4))**2-(pim2(3)+pim4(3))**2
1304 % -(pim2(2)+pim4(2))**2-(pim2(1)+pim4(1))**2
1306 ps31=(pim3(4)+pim1(4))**2-(pim3(3)+pim1(3))**2
1307 % -(pim3(2)+pim1(2))**2-(pim3(1)+pim1(1))**2
1309 ps34=(pim3(4)+pim4(4))**2-(pim3(3)+pim4(3))**2
1310 % -(pim3(2)+pim4(2))**2-(pim3(1)+pim4(1))**2
1314 pd324=pim3(4)*(pim2(4)-pim4(4))-pim3(3)*(pim2(3)-pim4(3))
1315 % -pim3(2)*(pim2(2)-pim4(2))-pim3(1)*(pim2(1)-pim4(1))
1317 pd314=pim3(4)*(pim1(4)-pim4(4))-pim3(3)*(pim1(3)-pim4(3))
1318 % -pim3(2)*(pim1(2)-pim4(2))-pim3(1)*(pim1(1)-pim4(1))
1320 pd234=pim2(4)*(pim3(4)-pim4(4))-pim2(3)*(pim3(3)-pim4(3))
1321 % -pim2(2)*(pim3(2)-pim4(2))-pim2(1)*(pim3(1)-pim4(1))
1323 pd214=pim2(4)*(pim1(4)-pim4(4))-pim2(3)*(pim1(3)-pim4(3))
1324 % -pim2(2)*(pim1(2)-pim4(2))-pim2(1)*(pim1(1)-pim4(1))
1326 pd134=pim1(4)*(pim3(4)-pim4(4))-pim1(3)*(pim3(3)-pim4(3))
1327 % -pim1(2)*(pim3(2)-pim4(2))-pim1(1)*(pim3(1)-pim4(1))
1329 pd124=pim1(4)*(pim2(4)-pim4(4))-pim1(3)*(pim2(3)-pim4(3))
1330 % -pim1(2)*(pim2(2)-pim4(2))-pim1(1)*(pim2(1)-pim4(1))
1334 qp1=pim1(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1335 % -pim1(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1336 % -pim1(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1337 % -pim1(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1339 qp2=pim2(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1340 % -pim2(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1341 % -pim2(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1342 % -pim2(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1344 qp3=pim3(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1345 % -pim3(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1346 % -pim3(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1347 % -pim3(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1349 qp4=pim4(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1350 % -pim4(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1351 % -pim4(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1352 % -pim4(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1357 t324=bwiga1(qmp3)*fpikm(sqrt(ps24),ampi,ampi)*gamma1
1358 t314=bwiga1(qmp3)*fpikm(sqrt(ps14),ampi,ampi)*gamma1
1359 t234=bwiga1(qmp2)*fpikm(sqrt(ps34),ampi,ampi)*gamma1
1360 t214=bwiga1(qmp2)*fpikm(sqrt(ps14),ampi,ampi)*gamma1
1361 t134=bwiga1(qmp1)*fpikm(sqrt(ps34),ampi,ampi)*gamma1
1362 t124=bwiga1(qmp1)*fpikm(sqrt(ps24),ampi,ampi)*gamma1
1366 s1423=frho4(sqrt(ps14),ampi,ampi)*bwigeps(ps23)*gamma2
1367 s2431=frho4(sqrt(ps24),ampi,ampi)*bwigeps(ps31)*gamma2
1368 s3421=frho4(sqrt(ps34),ampi,ampi)*bwigeps(ps21)*gamma2
1373 brack1=t234+t324+2.*t314+t134+2.*t214+t124
1374 % +t134*pd234/qmp1+t124*pd324/qmp1
1375 % -3./2.*(s3421+s2431+2.*s1423)
1378 brack2=t234*(1.+2.*pd134/qmp2)+3.*t324+3.*t124
1379 % +t134*(1.-pd234/qmp1)+2.*t214*pd314/qmp2
1381 % -3./2.*(s3421+3.*s2431)
1383 brack3=t324*(1.+2.*pd124/qmp3)+3.*t234+3.*t134
1384 % +t124*(1.-pd324/qmp1)+2.*t314*pd214/qmp3
1386 % -3./2.*(3.*s3421+s2431)
1388 brack4a=2.*t234*(1./2.+pd134/qq*(qp2/qmp2+1.)+pd234/qq)
1389 % +2.*t324*(1./2.+pd124/qq*(qp3/qmp3+1.)+pd324/qq)
1390 % +2.*t134*(1./2.+pd234/qq*(qp1/qmp1+1.)
1391 % -1./2.*pd234/qmp1+pd134/qq)
1392 % +2.*t124*(1./2.+pd324/qq*(qp1/qmp1+1.)
1393 % -1./2.*pd324/qmp1+pd124/qq)
1394 % +2.*t214*(pd314/qq*(qp2/qmp2+1.)+pd214/qq)
1395 % +2.*t314*(pd214/qq*(qp3/qmp3+1.)+pd314/qq)
1397 brack4b=-3./2.*(s3421*(2.*(qp3-qp4)/qq+1.)
1398 % +s2431*(2.*(qp2-qp4)/qq+1.)
1399 % +s1423*2.*(qp1-qp4)/qq)
1402 brack4=brack4a+brack4b
1406 hcomp1(k)=(pim1(k)-pim4(k))*brack1
1407 hcomp2(k)=pim2(k)*brack2
1408 hcomp3(k)=pim3(k)*brack3
1409 hcomp4(k)=(pim1(k)+pim2(k)+pim3(k)+pim4(k))*brack4
1415 hadcur(i)=hcomp1(i)+hcomp2(i)+hcomp3(i)-hcomp4(i)
1416 hadcur(i)=coef1*frho4(sqrt(qq),ampi,ampi)*hadcur(i)