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))
37 COMPLEX FUNCTION fk1ab(XMSQ,INDX)
42 COMPLEX F1,F2,AMPA,AMPB
62 ampa = cmplx(pkorb(3,81),0.)
63 ampb = cmplx(pkorb(3,82),0.)
64 ELSE IF (indx.EQ.2)
THEN 65 ampa = cmplx(pkorb(3,83),0.)
66 ampb = cmplx(pkorb(3,84),0.)
67 ELSEIF (indx.EQ.3)
THEN 68 ampa = cmplx(pkorb(3,85),0.)
69 ampb = cmplx(pkorb(3,86),0.)
70 ELSEIF (indx.EQ.4)
THEN 71 ampa = cmplx(pkorb(3,87),0.)
72 ampb = cmplx(pkorb(3,88),0.)
78 f1 = cmplx(-xm1sq,0.0)/cmplx(xmsq-xm1sq,fg1)
79 f2 = cmplx(-xm2sq,0.0)/cmplx(xmsq-xm2sq,fg2)
80 fk1ab = ampa*f1+ampb*f2
85 FUNCTION form1(MNUM,QQ,S1,SDWA)
95 COMPLEX FORM1,WIGNER,WIGFOR,FPIKM,BWIGM
96 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
97 * ,ampiz,ampi,amro,gamro,ama1,gama1
98 * ,amk,amkz,amkst,gamkst
100 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
101 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
102 * ,AMK,AMKZ,AMKST,GAMKST
104 COMPLEX FORMA1,FORMK1,FORMRO,FORMKS
105 COMPLEX FA1A1P,FK1AB,F3PI
113 form1 = f3pi(1,qq,s1,sdwa)
116 ELSEIF (mnum.EQ.1)
THEN 118 formks = bwigm(s1,amkst,gamkst,ampi,amk)
120 form1 = forma1*formks
122 ELSEIF (mnum.EQ.2)
THEN 124 formks = bwigm(s1,amkst,gamkst,ampi,amk)
126 form1 = forma1*formks
128 ELSEIF (mnum.EQ.3)
THEN 130 formks = bwigm(s1,amkst,gamkst,ampi,amk)
132 form1 = forma1*formks
134 ELSEIF (mnum.EQ.4)
THEN 136 formks = bwigm(s1,amkst,gamkst,ampi,amk)
138 form1 = formk1*formks
140 ELSEIF (mnum.EQ.5)
THEN 143 formro = fpikm(sqrt(s1),ampi,ampi)
144 form1 = formk1*formro
146 ELSEIF (mnum.EQ.6)
THEN 149 formks = bwigm(s1,amkst,gamkst,amk,ampi)
150 form1 = formk1*formks
152 ELSEIF (mnum.EQ.7)
THEN 155 ELSEIF (mnum.EQ.9)
THEN 161 form1 = f3pi(1,qq,s1,sdwa)
163 ELSEIF (mnum.gt.9.and.mnum.lt.20 )
THEN 170 form1 = f3pi(1,qq,s1,sdwa)
175 FUNCTION form2(MNUM,QQ,S1,SDWA)
185 COMPLEX FORM2,WIGNER,WIGFOR,FPIKM,BWIGM
186 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
187 * ,ampiz,ampi,amro,gamro,ama1,gama1
188 * ,amk,amkz,amkst,gamkst
190 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
191 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
192 * ,AMK,AMKZ,AMKST,GAMKST
194 COMPLEX FORMA1,FORMK1,FORMRO,FORMKS
195 COMPLEX FA1A1P,FK1AB,F3PI
203 form2 = f3pi(2,qq,s1,sdwa)
204 ELSEIF (mnum.EQ.1)
THEN 206 formro = fpikm(sqrt(s1),amk,amk)
208 form2 = forma1*formro
210 ELSEIF (mnum.EQ.2)
THEN 212 formro = fpikm(sqrt(s1),amk,amk)
214 form2 = forma1*formro
216 ELSEIF (mnum.EQ.3)
THEN 218 formro = fpikm(sqrt(s1),amk,amk)
220 form2 = forma1*formro
222 ELSEIF (mnum.EQ.4)
THEN 224 formks = bwigm(s1,amkst,gamkst,ampi,amk)
226 form2 = formk1*formks
228 ELSEIF (mnum.EQ.5)
THEN 230 formks = bwigm(s1,amkst,gamkst,ampi,amk)
232 form2 = formk1*formks
234 ELSEIF (mnum.EQ.6)
THEN 236 formro = fpikm(sqrt(s1),ampi,ampi)
238 form2 = formk1*formro
240 ELSEIF (mnum.EQ.7)
THEN 243 ELSEIF (mnum.EQ.9)
THEN 249 form2 = f3pi(2,qq,s1,sdwa)
250 ELSEIF (mnum.gt.9.and.mnum.lt.20 )
THEN 256 form2 = f3pi(2,qq,s1,sdwa)
262 COMPLEX FUNCTION bwigm(S,M,G,XM1,XM2)
275 IF (s.GT.(xm1+xm2)**2)
THEN 276 qs=sqrt(abs((s -(xm1+xm2)**2)*(s -(xm1-xm2)**2)))/sqrt(s)
277 qm=sqrt(abs((m**2-(xm1+xm2)**2)*(m**2-(xm1-xm2)**2)))/m
279 gs=g*(m/w)**2*(qs/qm)**3
283 bwigm=m**2/cmplx(m**2-s,-sqrt(s)*gs)
286 COMPLEX FUNCTION fpikm(W,XM1,XM2)
291 REAL ROM,ROG,ROM1,ROG1,BETA1,PI,PIM,S,W
308 fpikm=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2))
312 COMPLEX FUNCTION fpikmd(W,XM1,XM2)
317 REAL ROM,ROG,ROM1,ROG1,PI,PIM,S,W
337 fpikmd=(delta*bwigm(s,rom,rog,xm1,xm2)
338 $ +beta*bwigm(s,rom1,rog1,xm1,xm2)
339 $ + bwigm(s,rom2,rog2,xm1,xm2))
344 FUNCTION form3(MNUM,QQ,S1,SDWA)
354 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
355 * ,ampiz,ampi,amro,gamro,ama1,gama1
356 * ,amk,amkz,amkst,gamkst
358 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
359 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
360 * ,AMK,AMKZ,AMKST,GAMKST
363 COMPLEX FORMA1,FORMK1,FORMRO,FORMKS
364 COMPLEX FA1A1P,FK1AB,F3PI
372 form3 = f3pi(3,qq,s1,sdwa)
373 ELSEIF (mnum.EQ.3)
THEN 375 formks = bwigm(s1,amkst,gamkst,ampiz,amk)
377 form3 = forma1*formks
379 ELSEIF (mnum.EQ.6)
THEN 381 formks = bwigm(s1,amkst,gamkst,amk,ampi)
383 form3 = formk1*formks
384 ELSEIF (mnum.EQ.9)
THEN 390 form3 = f3pi(3,qq,s1,sdwa)
391 ELSEIF (mnum.gt.9.and.mnum.lt.20 )
THEN 397 form3 = f3pi(3,qq,s1,sdwa)
405 FUNCTION form4(MNUM,QQ,S1,S2,S3)
413 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
414 * ,ampiz,ampi,amro,gamro,ama1,gama1
415 * ,amk,amkz,amkst,gamkst
417 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
418 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
419 * ,AMK,AMKZ,AMKST,GAMKST
420 COMPLEX FORM4,WIGNER,FPIKM
429 FUNCTION form5(MNUM,QQ,S1,S2)
438 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
439 * ,ampiz,ampi,amro,gamro,ama1,gama1
440 * ,amk,amkz,amkst,gamkst
442 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
443 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
444 * ,AMK,AMKZ,AMKST,GAMKST
445 COMPLEX FORM5,WIGNER,FPIKM,FPIKMD,BWIGM
450 ELSEIF (mnum.EQ.1)
THEN 453 form5=fpikmd(sqrt(qq),ampi,ampi)/(1+elpha)
454 $ *( fpikm(sqrt(s2),ampi,ampi)
455 $ +elpha*bwigm(s1,amkst,gamkst,ampi,amk))
456 ELSEIF (mnum.EQ.2)
THEN 459 form5=fpikmd(sqrt(qq),ampi,ampi)/(1+elpha)
460 $ *( fpikm(sqrt(s2),ampi,ampi)
461 $ +elpha*bwigm(s1,amkst,gamkst,ampi,amk))
462 ELSEIF (mnum.EQ.3)
THEN 465 ELSEIF (mnum.EQ.4)
THEN 468 ELSEIF (mnum.EQ.5)
THEN 471 form5=bwigm(qq,amkst,gamkst,ampi,amk)/(1+elpha)
472 $ *( fpikm(sqrt(s1),ampi,ampi)
473 $ +elpha*bwigm(s2,amkst,gamkst,ampi,amk))
474 ELSEIF (mnum.EQ.6)
THEN 477 form5=bwigm(qq,amkst,gamkst,ampi,amkz)/(1+elpha)
478 $ *( fpikm(sqrt(s2),ampi,ampi)
479 $ +elpha*bwigm(s1,amkst,gamkst,ampi,amk))
480 ELSEIF (mnum.EQ.7)
THEN 482 form5=fpikmd(sqrt(qq),ampi,ampi)*fpikm(sqrt(s1),ampi,ampi)
483 ELSEIF (mnum.EQ.9)
THEN 486 ELSEIF (mnum.gt.9.and.mnum.lt.20 )
THEN 493 SUBROUTINE currx(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
502 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
503 * ,ampiz,ampi,amro,gamro,ama1,gama1
504 * ,amk,amkz,amkst,gamkst
506 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
507 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
508 * ,AMK,AMKZ,AMKST,GAMKST
509 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
510 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
513 COMMON /arbit/ arflat,aromeg
515 REAL PIM1(4),PIM2(4),PIM3(4),PIM4(4),PAA(4)
517 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FPIKM
522 DATA pi /3.141592653589793238462643/
524 bwign(a,xm,xg)=1.0/cmplx(a-xm**2,xm*xg)
544 coef1=2.0*sqrt(3.0)/fpi**2*arflat
551 paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)
561 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
564 sk=(pp(k,4)+pim4(4))**2-(pp(k,3)+pim4(3))**2
565 $ -(pp(k,2)+pim4(2))**2-(pp(k,1)+pim4(1))**2
575 denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
576 $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
582 $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
588 form1= fpikm(sqrt(sk),ampi,ampi) *fpikm(sqrt(qq),ampi,ampi)
598 $ hadcur(i)+cmplx(fix*coef1)*form1*aa(i,j)*(pp(k,j)-pp(4,j))
624 IF (k.EQ.4) sign= 1.0
625 qqa=qqa+sign*(paa(k)-pa(k))**2
626 ss23=ss23+sign*(pb(k) +pim3(k))**2
627 ss24=ss24+sign*(pb(k) +pim4(k))**2
628 ss34=ss34+sign*(pim3(k)+pim4(k))**2
629 qp1p2=qp1p2+sign*(paa(k)-pa(k))*pb(k)
630 qp1p3=qp1p3+sign*(paa(k)-pa(k))*pim3(k)
631 qp1p4=qp1p4+sign*(paa(k)-pa(k))*pim4(k)
632 p1p2=p1p2+sign*pa(k)*pb(k)
633 p1p3=p1p3+sign*pa(k)*pim3(k)
634 p1p4=p1p4+sign*pa(k)*pim4(k)
637 form2=coef2*(bwign(qq,amro,gamro)+elpha*bwign(qq,amrop,gamrop))
640 form3=bwign(qqa,amom,gamom)
643 hadcur(k)=hadcur(k)+form2*form3*(
644 $ pb(k)*(qp1p3*p1p4-qp1p4*p1p3)
645 $ +pim3(k)*(qp1p4*p1p2-qp1p2*p1p4)
646 $ +pim4(k)*(qp1p2*p1p3-qp1p3*p1p2) )
654 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
657 sk=(pp(k,4)+pim4(4))**2-(pp(k,3)+pim4(3))**2
658 $ -(pp(k,2)+pim4(2))**2-(pp(k,1)+pim4(1))**2
669 denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
670 $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
676 $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
682 form1= fpikm(sqrt(sk),ampi,ampi) *fpikm(sqrt(qq),ampi,ampi)
689 $ hadcur(i)+cmplx(coef1)*form1*aa(i,j)*(pp(k,j)-pp(4,j))
695 FUNCTION wigfor(S,XM,XGAM)
696 COMPLEX WIGFOR,WIGNOR
697 wignor=cmplx(-xm**2,xm*xgam)
698 wigfor=wignor/cmplx(s-xm**2,xm*xgam)
704 COMMON /inout/ inut, iout
705 WRITE (unit = iout,fmt = 99)
706 WRITE (unit = iout,fmt = 98)
709 . /,
' *************************************************** ',
710 . /,
' YOU ARE USING THE 4 PION DECAY MODE FORM FACTORS ',
711 . /,
' WHICH HAVE BEEN DESCRIBED IN:',
712 . /,
' R. DECKER, M. FINKEMEIER, P. HEILIGER AND H.H. JONSSON',
713 ' "TAU DECAYS INTO FOUR PIONS" ',
714 . /,
' UNIVERSITAET KARLSRUHE PREPRINT TTP 94-13 (1994);',
715 . /,
' LNF-94/066(IR); HEP-PH/9410260 ',
717 . /,
' PLEASE NOTE THAT THIS ROUTINE IS USING PARAMETERS',
718 . /,
' RELATED TO THE 3 PION DECAY MODE (A1 MODE), SUCH AS',
719 . /,
' THE A1 MASS AND WIDTH (TAKEN FROM THE COMMON /PARMAS/)',
720 . /,
' AND THE 2 PION VECTOR RESONANCE FORM FACTOR (BY USING',
721 . /,
' THE ROUTINE FPIKM)' ,
722 . /,
' THUS IF YOU DECIDE TO CHANGE ANY OF THESE, YOU WILL' ,
723 . /,
' HAVE TO REFIT THE 4 PION PARAMETERS IN THE COMMON' )
725 .
' BLOCK /TAU4PI/, OR YOU MIGHT GET A BAD DISCRIPTION' ,
726 . /,
' OF TAU -> 4 PIONS' ,
727 . /,
' for these formfactors set in routine CHOICE for',
728 . /, .eq.
' mnum102 -- AMRX=1.42 and GAMRX=.21',
729 . /, .eq.
' mnum101 -- AMRX=1.3 and GAMRX=.46 PROB1,PROB2=0.2',
730 . /,
' to optimize phase space parametrization',
731 . /,
' *************************************************** ',
732 . /,
' coded by M. Finkemeier and P. Heiliger, 29. sept. 1994',
733 . /,
' incorporated to TAUOLA by Z. Was 17. jan. 1995',
736 . /,
' changed by: Z. Was on 17.01.95',
737 . /,
' changes by: M. Finkemeier on 30.01.95' )
741 COMMON /tau4pi/ gomega,gamma1,gamma2,rom1,rog1,beta1,
743 REAL*4 GOMEGA,GAMMA1,GAMMA2,ROM1,ROG1,BETA1,
755 COMPLEX FUNCTION bwiga1(QA)
760 COMMON / parmas/ amtau,amnuta,amel,amnue,ammu,amnumu,
761 % AMPIZ,ampi,amro,gamro,ama1,gama1,
762 % AMK,amkz,amkst,gamkst
765 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU,
766 % AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1,
767 % AMK,AMKZ,AMKST,GAMKST
769 wigner(a,b,c)=cmplx(1.0,0.0)/cmplx(a-b**2,b*c)
770 gamax=gama1*gfun(qa)/gfun(ama1**2)
771 bwiga1=-ama1**2*wigner(qa,ama1,gamax)
774 COMPLEX FUNCTION bwigeps(QEPS)
781 bwigeps=cmplx(meps**2,-meps*geps)/
782 % CMPLX(meps**2-qeps,-meps*geps)
785 COMPLEX FUNCTION frho4(W,XM1,XM2)
791 COMMON /tau4pi/ gomega,gamma1,gamma2,rom1,rog1,beta1,
793 REAL*4 GOMEGA,GAMMA1,GAMMA2,ROM1,ROG1,BETA1,
795 REAL ROM,ROG,PI,PIM,S,W
810 frho4=(bwigm(s,rom,rog,xm1,xm2)+beta1*bwigm(s,rom1,rog1,xm1,xm2)
811 & +beta2*bwigm(s,rom2,rog2,xm1,xm2))
815 SUBROUTINE curr(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
826 COMMON /tau4pi/ gomega,gamma1,gamma2,rom1,rog1,beta1,
828 REAL*4 GOMEGA,GAMMA1,GAMMA2,ROM1,ROG1,BETA1,
830 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
831 * ,ampiz,ampi,amro,gamro,ama1,gama1
832 * ,amk,amkz,amkst,gamkst
834 REAL*4 AMTAU,AMNUTA,AMEL,AMNUE,AMMU,AMNUMU
835 * ,AMPIZ,AMPI,AMRO,GAMRO,AMA1,GAMA1
836 * ,AMK,AMKZ,AMKST,GAMKST
837 COMMON / decpar / gfermi,gv,ga,ccabib,scabib,gamel
838 REAL*4 GFERMI,GV,GA,CCABIB,SCABIB,GAMEL
839 REAL PIM1(4),PIM2(4),PIM3(4),PIM4(4),PAA(4)
840 COMPLEX HADCUR(4),FORM1,FORM2,FORM3,FPIKM
842 COMPLEX BWIGEPS,BWIGA1
843 COMPLEX HCOMP1(4),HCOMP2(4),HCOMP3(4),HCOMP4(4)
845 COMPLEX T243,T213,T143,T123,T341,T342
846 COMPLEX T124,T134,T214,T234,T314,T324
847 COMPLEX S2413,S2314,S1423,S1324,S34
849 COMPLEX BRACK1,BRACK2,BRACK3,BRACK4A,BRACK4B,BRACK4
851 REAL QMP1,QMP2,QMP3,QMP4
852 REAL PS43,PS41,PS42,PS34,PS14,PS13,PS24,PS23
855 REAL PD243,PD241,PD213,PD143,PD142
856 REAL PD123,PD341,PD342,PD413,PD423
857 REAL PD124,PD134,PD214,PD234,PD314,PD324
862 DATA pi /3.141592653589793238462643/
865 bwign(a,xm,xg)=1.0/cmplx(a-xm**2,xm*xg)
886 coef1=2.0*sqrt(3.0)/fpi**2*arflat
893 paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)
903 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
909 qmp1=(pim2(4)+pim3(4)+pim4(4))**2-(pim2(3)+pim3(3)+pim4(3))**2
910 % -(pim2(2)+pim3(2)+pim4(2))**2-(pim2(1)+pim3(1)+pim4(1))**2
912 qmp2=(pim1(4)+pim3(4)+pim4(4))**2-(pim1(3)+pim3(3)+pim4(3))**2
913 % -(pim1(2)+pim3(2)+pim4(2))**2-(pim1(1)+pim3(1)+pim4(1))**2
915 qmp3=(pim1(4)+pim2(4)+pim4(4))**2-(pim1(3)+pim2(3)+pim4(3))**2
916 % -(pim1(2)+pim2(2)+pim4(2))**2-(pim1(1)+pim2(1)+pim4(1))**2
918 qmp4=(pim1(4)+pim2(4)+pim3(4))**2-(pim1(3)+pim2(3)+pim3(3))**2
919 % -(pim1(2)+pim2(2)+pim3(2))**2-(pim1(1)+pim2(1)+pim3(1))**2
924 ps43=(pim4(4)+pim3(4))**2-(pim4(3)+pim3(3))**2
925 % -(pim4(2)+pim3(2))**2-(pim4(1)+pim3(1))**2
927 ps41=(pim4(4)+pim1(4))**2-(pim4(3)+pim1(3))**2
928 % -(pim4(2)+pim1(2))**2-(pim4(1)+pim1(1))**2
930 ps42=(pim4(4)+pim2(4))**2-(pim4(3)+pim2(3))**2
931 % -(pim4(2)+pim2(2))**2-(pim4(1)+pim2(1))**2
937 ps13=(pim1(4)+pim3(4))**2-(pim1(3)+pim3(3))**2
938 % -(pim1(2)+pim3(2))**2-(pim1(1)+pim3(1))**2
942 ps23=(pim2(4)+pim3(4))**2-(pim2(3)+pim3(3))**2
943 % -(pim2(2)+pim3(2))**2-(pim2(1)+pim3(1))**2
945 pd243=pim2(4)*(pim4(4)-pim3(4))-pim2(3)*(pim4(3)-pim3(3))
946 % -pim2(2)*(pim4(2)-pim3(2))-pim2(1)*(pim4(1)-pim3(1))
948 pd241=pim2(4)*(pim4(4)-pim1(4))-pim2(3)*(pim4(3)-pim1(3))
949 % -pim2(2)*(pim4(2)-pim1(2))-pim2(1)*(pim4(1)-pim1(1))
951 pd213=pim2(4)*(pim1(4)-pim3(4))-pim2(3)*(pim1(3)-pim3(3))
952 % -pim2(2)*(pim1(2)-pim3(2))-pim2(1)*(pim1(1)-pim3(1))
954 pd143=pim1(4)*(pim4(4)-pim3(4))-pim1(3)*(pim4(3)-pim3(3))
955 % -pim1(2)*(pim4(2)-pim3(2))-pim1(1)*(pim4(1)-pim3(1))
957 pd142=pim1(4)*(pim4(4)-pim2(4))-pim1(3)*(pim4(3)-pim2(3))
958 % -pim1(2)*(pim4(2)-pim2(2))-pim1(1)*(pim4(1)-pim2(1))
960 pd123=pim1(4)*(pim2(4)-pim3(4))-pim1(3)*(pim2(3)-pim3(3))
961 % -pim1(2)*(pim2(2)-pim3(2))-pim1(1)*(pim2(1)-pim3(1))
963 pd341=pim3(4)*(pim4(4)-pim1(4))-pim3(3)*(pim4(3)-pim1(3))
964 % -pim3(2)*(pim4(2)-pim1(2))-pim3(1)*(pim4(1)-pim1(1))
966 pd342=pim3(4)*(pim4(4)-pim2(4))-pim3(3)*(pim4(3)-pim2(3))
967 % -pim3(2)*(pim4(2)-pim2(2))-pim3(1)*(pim4(1)-pim2(1))
969 pd413=pim4(4)*(pim1(4)-pim3(4))-pim4(3)*(pim1(3)-pim3(3))
970 % -pim4(2)*(pim1(2)-pim3(2))-pim4(1)*(pim1(1)-pim3(1))
972 pd423=pim4(4)*(pim2(4)-pim3(4))-pim4(3)*(pim2(3)-pim3(3))
973 % -pim4(2)*(pim2(2)-pim3(2))-pim4(1)*(pim2(1)-pim3(1))
977 qp1=pim1(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
978 % -pim1(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
979 % -pim1(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
980 % -pim1(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
982 qp2=pim2(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
983 % -pim2(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
984 % -pim2(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
985 % -pim2(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
987 qp3=pim3(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
988 % -pim3(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
989 % -pim3(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
990 % -pim3(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
992 qp4=pim4(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
993 % -pim4(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
994 % -pim4(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
995 % -pim4(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1001 t243=bwiga1(qmp2)*fpikm(sqrt(ps43),ampi,ampi)*gamma1
1002 t213=bwiga1(qmp2)*fpikm(sqrt(ps13),ampi,ampi)*gamma1
1003 t143=bwiga1(qmp1)*fpikm(sqrt(ps43),ampi,ampi)*gamma1
1004 t123=bwiga1(qmp1)*fpikm(sqrt(ps23),ampi,ampi)*gamma1
1005 t341=bwiga1(qmp3)*fpikm(sqrt(ps41),ampi,ampi)*gamma1
1006 t342=bwiga1(qmp3)*fpikm(sqrt(ps42),ampi,ampi)*gamma1
1010 s2413=frho4(sqrt(ps24),ampi,ampi)*gamma2
1011 s2314=frho4(sqrt(ps23),ampi,ampi)*bwigeps(ps14)*gamma2
1012 s1423=frho4(sqrt(ps14),ampi,ampi)*gamma2
1013 s1324=frho4(sqrt(ps13),ampi,ampi)*bwigeps(ps24)*gamma2
1014 s34=frho4(sqrt(ps34),ampi,ampi)*gamma2
1018 brack1=2.*t143+2.*t243+t123+t213
1019 % +t341*(pd241/qmp3-1.)+t342*(pd142/qmp3-1.)
1020 % +3./4.*(s1423+s2413-s2314-s1324)-3.*s34
1022 brack2=2.*t143*pd243/qmp1+3.*t213
1023 % +t123*(2.*pd423/qmp1+1.)+t341*(pd241/qmp3+3.)
1024 % +t342*(pd142/qmp3+1.)
1025 % -3./4.*(s2314+3.*s1324+3.*s1423+s2413)
1027 brack3=2.*t243*pd143/qmp2+3.*t123
1028 % +t213*(2.*pd413/qmp2+1.)+t341*(pd241/qmp3+1.)
1029 % +t342*(pd142/qmp3+3.)
1030 % -3./4.*(3.*s2314+s1324+s1423+3.*s2413)
1032 brack4a=2.*t143*(pd243/qq*(qp1/qmp1+1.)+pd143/qq)
1033 % +2.*t243*(pd143/qq*(qp2/qmp2+1.)+pd243/qq)
1035 % +2.*t123*(pd423/qq*(qp1/qmp1+1.)+pd123/qq)
1036 % +2.*t213*(pd413/qq*(qp2/qmp2+1.)+pd213/qq)
1037 % +t341*(pd241/qmp3+1.-2.*pd241/qq*(qp3/qmp3+1.)
1039 % +t342*(pd142/qmp3+1.-2.*pd142/qq*(qp3/qmp3+1.)
1042 brack4b=-3./4.*(s2314*(2.*(qp2-qp3)/qq+1.)
1043 % +s1324*(2.*(qp1-qp3)/qq+1.)
1044 % +s1423*(2.*(qp1-qp4)/qq+1.)
1045 % +s2413*(2.*(qp2-qp4)/qq+1.)
1046 % +4.*s34*(qp4-qp3)/qq)
1048 brack4=brack4a+brack4b
1052 hcomp1(k)=(pim3(k)-pim4(k))*brack1
1053 hcomp2(k)=pim1(k)*brack2
1054 hcomp3(k)=pim2(k)*brack3
1055 hcomp4(k)=(pim1(k)+pim2(k)+pim3(k)+pim4(k))*brack4
1061 hadcur(i)=hcomp1(i)-hcomp2(i)-hcomp3(i)+hcomp4(i)
1062 hadcur(i)=-coef1*frho4(sqrt(qq),ampi,ampi)*hadcur(i)
1091 IF (k.EQ.4) sign= 1.0
1092 qqa=qqa+sign*(paa(k)-pa(k))**2
1093 ss23=ss23+sign*(pb(k) +pim3(k))**2
1094 ss24=ss24+sign*(pb(k) +pim4(k))**2
1095 ss34=ss34+sign*(pim3(k)+pim4(k))**2
1096 qp1p2=qp1p2+sign*(paa(k)-pa(k))*pb(k)
1097 qp1p3=qp1p3+sign*(paa(k)-pa(k))*pim3(k)
1098 qp1p4=qp1p4+sign*(paa(k)-pa(k))*pim4(k)
1099 p1p2=p1p2+sign*pa(k)*pb(k)
1100 p1p3=p1p3+sign*pa(k)*pim3(k)
1101 p1p4=p1p4+sign*pa(k)*pim4(k)
1104 form2=coef2*(bwign(qq,amro,gamro)+elpha*bwign(qq,amrop,gamrop))
1107 form3=bwign(qqa,amom,gamom)
1110 hadcur(k)=hadcur(k)+form2*form3*(
1111 $ pb(k)*(qp1p3*p1p4-qp1p4*p1p3)
1112 $ +pim3(k)*(qp1p4*p1p2-qp1p2*p1p4)
1113 $ +pim4(k)*(qp1p2*p1p3-qp1p3*p1p2) )
1121 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
1128 qmp1=(pim2(4)+pim3(4)+pim4(4))**2-(pim2(3)+pim3(3)+pim4(3))**2
1129 % -(pim2(2)+pim3(2)+pim4(2))**2-(pim2(1)+pim3(1)+pim4(1))**2
1131 qmp2=(pim1(4)+pim3(4)+pim4(4))**2-(pim1(3)+pim3(3)+pim4(3))**2
1132 % -(pim1(2)+pim3(2)+pim4(2))**2-(pim1(1)+pim3(1)+pim4(1))**2
1134 qmp3=(pim1(4)+pim2(4)+pim4(4))**2-(pim1(3)+pim2(3)+pim4(3))**2
1135 % -(pim1(2)+pim2(2)+pim4(2))**2-(pim1(1)+pim2(1)+pim4(1))**2
1137 qmp4=(pim1(4)+pim2(4)+pim3(4))**2-(pim1(3)+pim2(3)+pim3(3))**2
1138 % -(pim1(2)+pim2(2)+pim3(2))**2-(pim1(1)+pim2(1)+pim3(1))**2
1143 ps14=(pim1(4)+pim4(4))**2-(pim1(3)+pim4(3))**2
1144 % -(pim1(2)+pim4(2))**2-(pim1(1)+pim4(1))**2
1146 ps21=(pim2(4)+pim1(4))**2-(pim2(3)+pim1(3))**2
1147 % -(pim2(2)+pim1(2))**2-(pim2(1)+pim1(1))**2
1149 ps23=(pim2(4)+pim3(4))**2-(pim2(3)+pim3(3))**2
1150 % -(pim2(2)+pim3(2))**2-(pim2(1)+pim3(1))**2
1152 ps24=(pim2(4)+pim4(4))**2-(pim2(3)+pim4(3))**2
1153 % -(pim2(2)+pim4(2))**2-(pim2(1)+pim4(1))**2
1155 ps31=(pim3(4)+pim1(4))**2-(pim3(3)+pim1(3))**2
1156 % -(pim3(2)+pim1(2))**2-(pim3(1)+pim1(1))**2
1158 ps34=(pim3(4)+pim4(4))**2-(pim3(3)+pim4(3))**2
1159 % -(pim3(2)+pim4(2))**2-(pim3(1)+pim4(1))**2
1163 pd324=pim3(4)*(pim2(4)-pim4(4))-pim3(3)*(pim2(3)-pim4(3))
1164 % -pim3(2)*(pim2(2)-pim4(2))-pim3(1)*(pim2(1)-pim4(1))
1166 pd314=pim3(4)*(pim1(4)-pim4(4))-pim3(3)*(pim1(3)-pim4(3))
1167 % -pim3(2)*(pim1(2)-pim4(2))-pim3(1)*(pim1(1)-pim4(1))
1169 pd234=pim2(4)*(pim3(4)-pim4(4))-pim2(3)*(pim3(3)-pim4(3))
1170 % -pim2(2)*(pim3(2)-pim4(2))-pim2(1)*(pim3(1)-pim4(1))
1172 pd214=pim2(4)*(pim1(4)-pim4(4))-pim2(3)*(pim1(3)-pim4(3))
1173 % -pim2(2)*(pim1(2)-pim4(2))-pim2(1)*(pim1(1)-pim4(1))
1175 pd134=pim1(4)*(pim3(4)-pim4(4))-pim1(3)*(pim3(3)-pim4(3))
1176 % -pim1(2)*(pim3(2)-pim4(2))-pim1(1)*(pim3(1)-pim4(1))
1178 pd124=pim1(4)*(pim2(4)-pim4(4))-pim1(3)*(pim2(3)-pim4(3))
1179 % -pim1(2)*(pim2(2)-pim4(2))-pim1(1)*(pim2(1)-pim4(1))
1183 qp1=pim1(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1184 % -pim1(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1185 % -pim1(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1186 % -pim1(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1188 qp2=pim2(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1189 % -pim2(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1190 % -pim2(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1191 % -pim2(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1193 qp3=pim3(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1194 % -pim3(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1195 % -pim3(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1196 % -pim3(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1198 qp4=pim4(4)*(pim1(4)+pim2(4)+pim3(4)+pim4(4))
1199 % -pim4(3)*(pim1(3)+pim2(3)+pim3(3)+pim4(3))
1200 % -pim4(2)*(pim1(2)+pim2(2)+pim3(2)+pim4(2))
1201 % -pim4(1)*(pim1(1)+pim2(1)+pim3(1)+pim4(1))
1206 t324=bwiga1(qmp3)*fpikm(sqrt(ps24),ampi,ampi)*gamma1
1207 t314=bwiga1(qmp3)*fpikm(sqrt(ps14),ampi,ampi)*gamma1
1208 t234=bwiga1(qmp2)*fpikm(sqrt(ps34),ampi,ampi)*gamma1
1209 t214=bwiga1(qmp2)*fpikm(sqrt(ps14),ampi,ampi)*gamma1
1210 t134=bwiga1(qmp1)*fpikm(sqrt(ps34),ampi,ampi)*gamma1
1211 t124=bwiga1(qmp1)*fpikm(sqrt(ps24),ampi,ampi)*gamma1
1215 s1423=frho4(sqrt(ps14),ampi,ampi)*bwigeps(ps23)*gamma2
1216 s2431=frho4(sqrt(ps24),ampi,ampi)*bwigeps(ps31)*gamma2
1217 s3421=frho4(sqrt(ps34),ampi,ampi)*bwigeps(ps21)*gamma2
1222 brack1=t234+t324+2.*t314+t134+2.*t214+t124
1223 % +t134*pd234/qmp1+t124*pd324/qmp1
1224 % -3./2.*(s3421+s2431+2.*s1423)
1227 brack2=t234*(1.+2.*pd134/qmp2)+3.*t324+3.*t124
1228 % +t134*(1.-pd234/qmp1)+2.*t214*pd314/qmp2
1230 % -3./2.*(s3421+3.*s2431)
1232 brack3=t324*(1.+2.*pd124/qmp3)+3.*t234+3.*t134
1233 % +t124*(1.-pd324/qmp1)+2.*t314*pd214/qmp3
1235 % -3./2.*(3.*s3421+s2431)
1237 brack4a=2.*t234*(1./2.+pd134/qq*(qp2/qmp2+1.)+pd234/qq)
1238 % +2.*t324*(1./2.+pd124/qq*(qp3/qmp3+1.)+pd324/qq)
1239 % +2.*t134*(1./2.+pd234/qq*(qp1/qmp1+1.)
1240 % -1./2.*pd234/qmp1+pd134/qq)
1241 % +2.*t124*(1./2.+pd324/qq*(qp1/qmp1+1.)
1242 % -1./2.*pd324/qmp1+pd124/qq)
1243 % +2.*t214*(pd314/qq*(qp2/qmp2+1.)+pd214/qq)
1244 % +2.*t314*(pd214/qq*(qp3/qmp3+1.)+pd314/qq)
1246 brack4b=-3./2.*(s3421*(2.*(qp3-qp4)/qq+1.)
1247 % +s2431*(2.*(qp2-qp4)/qq+1.)
1248 % +s1423*2.*(qp1-qp4)/qq)
1251 brack4=brack4a+brack4b
1255 hcomp1(k)=(pim1(k)-pim4(k))*brack1
1256 hcomp2(k)=pim2(k)*brack2
1257 hcomp3(k)=pim3(k)*brack3
1258 hcomp4(k)=(pim1(k)+pim2(k)+pim3(k)+pim4(k))*brack4
1264 hadcur(i)=hcomp1(i)+hcomp2(i)+hcomp3(i)-hcomp4(i)
1265 hadcur(i)=coef1*frho4(sqrt(qq),ampi,ampi)*hadcur(i)