1 /* copyright(c) 1991-2012 free software foundation, inc.
2 this file is part of the gnu c library.
4 the gnu c library is free software; you can redistribute it and/or
5 modify it under the terms of the gnu lesser general
Public
6 license as published by the free software foundation; either
7 version 2.1 of the license, or(at your option) any later version.
9 the gnu c library is distributed in the hope that it will be useful,
10 but without any warranty; without even the implied warranty of
11 merchantability or fitness for a particular purpose. see the gnu
12 lesser general
Public license for more details.
14 you should have received a copy of the gnu lesser general
Public
15 license along with the gnu c library;
if not, see
16 <http://www.gnu.org/licenses/>. */
19 /* this header is separate from features.h so that the compiler can
20 include it implicitly at the start of every compilation. it must
21 not itself include <features.h> or any other header that includes
22 <features.h> because the
implicit include comes before any feature
23 test macros that may be defined in a source file before it first
24 explicitly includes a system header. gcc knows the name of this
25 header in order to preinclude it. */
27 /* we
do support the iec 559 math functionality,
real and complex. */
29 /* wchar_t uses iso/iec 10646 (2nd ed., published 2011-03-15) /
32 /* we
do not support c11 <threads.h>. */
34 *ajw 1 version of curr from koralb.
35 SUBROUTINE curr_cleo(MNUM,PIM1,PIM2,PIM3,PIM4,HADCUR)
36 c ==================================================================
37 c ajw, 11/97 - based on original curr from tauola:
38 c hadronic current for 4 pi final state
39 c r. fisher, j. wess and f. wagner z. phys c3(1980) 313
40 c r. decker z. phys c36(1987) 487.
41 c m. gell-mann, d. sharp, w. wagner phys. rev. lett 8 (1962) 261.
42 c but, rewritten to be more general and less
"theoretical",
43 c using parameters tuned by vasia and dsc.
44 c ==================================================================
46 COMMON / parmas / amtau,amnuta,amel,amnue,ammu,amnumu
47 * ,ampiz,ampi,amro,gamro,ama1,gama1
48 * ,amk,amkz,amkst,gamkst
50 REAL*4 amtau,amnuta,amel,amnue,ammu,amnumu
51 * ,ampiz,ampi,amro,gamro,ama1,gama1
52 * ,amk,amkz,amkst,gamkst
54 REAL pim1(4),pim2(4),pim3(4),pim4(4)
57 INTEGER k,l,mnum,k1,k2,iro,i,j,kk
58 REAL pa(4),pb(4),paa(4)
60 REAL a,xm,xg,g1,g2,g,amro2,gamro2,amro3,gamro3,amom,gamom
61 REAL fro,coef1,fpi,coef2,qq,sk,denom,sig,qqa,ss23,ss24,ss34,qp1p2
62 REAL qp1p3,qp1p4,p1p2,p1p3,p1p4,sign
64 COMPLEX alf0,alf1,alf2,alf3
65 COMPLEX lam0,lam1,lam2,lam3
66 COMPLEX bet1,bet2,bet3
67 COMPLEX form1,form2,form3,form4,form2pi
68 COMPLEX bwigm,wigfor,fpikm,fpikmd
72 bwign(a,xm,xg)=1.0/cmplx(a-xm**2,xm*xg)
73 c*******************************************************************************
75 c --- masses and constants
76 IF (g1.NE.12.924)
THEN
82 coef1=2.0*sqrt(3.0)/fpi**2
86 c masses and widths for for rho-prim and rho-bis:
99 c amplitudes for(pi-pi-pi0pi+) -> ps, rho0, rho-, rho+, omega.
100 ampl(1) = cmplx(pkorb(3,31)*coef1,0.)
101 ampl(2) = cmplx(pkorb(3,32)*coef1,0.)*cexp(cmplx(0.,pkorb(3,42)))
102 ampl(3) = cmplx(pkorb(3,33)*coef1,0.)*cexp(cmplx(0.,pkorb(3,43)))
103 ampl(4) = cmplx(pkorb(3,34)*coef1,0.)*cexp(cmplx(0.,pkorb(3,44)))
104 ampl(5) = cmplx(pkorb(3,35)*coef2,0.)*cexp(cmplx(0.,pkorb(3,45)))
105 c amplitudes for(pi0pi0pi0pi-) -> ps, rho-.
106 ampl(6) = cmplx(pkorb(3,36)*coef1)
107 ampl(7) = cmplx(pkorb(3,37)*coef1)
109 c rho
' contributions to rho' -> pi-omega:
110 alf0 = cmplx(pkorb(3,51),0.0)
111 alf1 = cmplx(pkorb(3,52)*amro**2,0.0)
112 alf2 = cmplx(pkorb(3,53)*amro2**2,0.0)
113 alf3 = cmplx(pkorb(3,54)*amro3**2,0.0)
114 c rho
' contribtions to rho' -> rhopipi:
115 lam0 = cmplx(pkorb(3,55),0.0)
116 lam1 = cmplx(pkorb(3,56)*amro**2,0.0)
117 lam2 = cmplx(pkorb(3,57)*amro2**2,0.0)
118 lam3 = cmplx(pkorb(3,58)*amro3**2,0.0)
119 c rho contributions to rhopipi, rho -> 2pi:
120 bet1 = cmplx(pkorb(3,59)*amro**2,0.0)
121 bet2 = cmplx(pkorb(3,60)*amro2**2,0.0)
122 bet3 = cmplx(pkorb(3,61)*amro3**2,0.0)
125 c**************************************************
127 c --- initialization of four vectors
132 paa(k)=pim1(k)+pim2(k)+pim3(k)+pim4(k)
139 c ===================================================================
140 c pi- pi- p0 pi+
case ====
141 c ===================================================================
142 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
144 c add m(4pi)-dependence to rhopipi channels:
145 form4= lam0+lam1*bwign(qq,amro,gamro)
146 * +lam2*bwign(qq,amro2,gamro2)
147 * +lam3*bwign(qq,amro3,gamro3)
149 c --- loop over five contributions of the rho-pi-pi
155 ELSEIF (k2.EQ.3)
THEN
159 ELSEIF (k1.EQ.3)
THEN
169 sk=(pp(k1,4)+pp(k2,4))**2-(pp(k1,3)+pp(k2,3))**2
170 $ -(pp(k1,2)+pp(k2,2))**2-(pp(k1,1)+pp(k2,1))**2
172 c -- definition of aa matrix
178 c ... and the rest ...
180 IF (l.NE.k1.AND.l.NE.k2)
THEN
181 denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
182 $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
188 $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
193 c --- lets add something to hadcurr
194 c form1= fpikm(sqrt(sk),ampi,ampi) *fpikm(sqrt(qq),ampi,ampi)
195 c form1= ampl(1)+ampr*fpikm(sqrt(sk),ampi,ampi)
197 form2pi= bet1*bwigm(sk,amro,gamro,ampa,ampi)
198 1 +bet2*bwigm(sk,amro2,gamro2,ampa,ampi)
199 2 +bet3*bwigm(sk,amro3,gamro3,ampa,ampi)
200 form1= ampl(1)+ampr*form2pi
204 hadcur(i)=hadcur(i)+form1*form4*aa(i,j)*(pp(k1,j)-pp(k2,j))
206 c ---
end of the rho-pi-pi current(5 possibilities)
209 c ===================================================================
210 c now modify the coefficient for the omega-pi current: =
211 c ===================================================================
212 IF (ampl(5).EQ.cmplx(0.,0.)) goto 311
214 c overall rho+rhoprime for the 4pi system:
215 c form2=ampl(5)*(bwign(qq,amro,gamro)+elpha*bwign(qq,amrop,gamrop))
216 c modified m(4pi)-dependence:
217 form2=ampl(5)*(alf0+alf1*bwign(qq,amro,gamro)
218 * +alf2*bwign(qq,amro2,gamro2)
219 * +alf3*bwign(qq,amro3,gamro3))
221 c --- there are two possibilities for omega current
222 c --- pa pb are corresponding first and second pi-s
228 c --- lorentz invariants
241 IF (k.EQ.4) sign= 1.0
242 qqa=qqa+sign*(paa(k)-pa(k))**2
243 ss23=ss23+sign*(pb(k) +pim3(k))**2
244 ss24=ss24+sign*(pb(k) +pim4(k))**2
245 ss34=ss34+sign*(pim3(k)+pim4(k))**2
246 qp1p2=qp1p2+sign*(paa(k)-pa(k))*pb(k)
247 qp1p3=qp1p3+sign*(paa(k)-pa(k))*pim3(k)
248 qp1p4=qp1p4+sign*(paa(k)-pa(k))*pim4(k)
249 p1p2=p1p2+sign*pa(k)*pb(k)
250 p1p3=p1p3+sign*pa(k)*pim3(k)
251 p1p4=p1p4+sign*pa(k)*pim4(k)
254 c omega -> rho pi for the 3pi system:
255 c form3=bwign(qqa,amom,gamom)*(bwign(ss23,amro,gamro)+
256 c $ bwign(ss24,amro,gamro)+bwign(ss34,amro,gamro))
257 c no omega -> rho pi; just straight omega:
258 form3=bwign(qqa,amom,gamom)
261 hadcur(k)=hadcur(k)+form2*form3*(
262 $ pb(k)*(qp1p3*p1p4-qp1p4*p1p3)
263 $ +pim3(k)*(qp1p4*p1p2-qp1p2*p1p4)
264 $ +pim4(k)*(qp1p2*p1p3-qp1p3*p1p2) )
270 c ===================================================================
271 c pi0 pi0 p0 pi-
case ====
272 c ===================================================================
273 qq=paa(4)**2-paa(3)**2-paa(2)**2-paa(1)**2
275 c --- loop over three contribution of the non-omega current
277 sk=(pp(k,4)+pim4(4))**2-(pp(k,3)+pim4(3))**2
278 $ -(pp(k,2)+pim4(2))**2-(pp(k,1)+pim4(1))**2
280 c -- definition of aa matrix
287 c ... and the rest ...
290 denom=(paa(4)-pp(l,4))**2-(paa(3)-pp(l,3))**2
291 $ -(paa(2)-pp(l,2))**2-(paa(1)-pp(l,1))**2
297 $ -sig*(paa(i)-2.0*pp(l,i))*(paa(j)-pp(l,j))/denom
302 c --- lets add something to hadcurr
303 c form1= fpikm(sqrt(sk),ampi,ampi) *fpikmd(sqrt(qq),ampi,ampi)
304 ccccccccccccc form1=wigfor(sk,amro,gamro) (tests)
305 c form1= fpikm(sqrt(sk),ampi,ampi) *fpikm(sqrt(qq),ampi,ampi)
306 form1 = ampl(6)+ampl(7)*fpikm(sqrt(sk),ampi,ampi)
310 hadcur(i)=hadcur(i)+form1*aa(i,j)*(pp(k,j)-pp(4,j))
312 c ---
end of the non omega current(3 possibilities)