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 a1 form factor
35 COMPLEX FUNCTION f3pi(IFORM,QQ,SA,SB)
36 c.......................................................................
38 c. f3pi - 1 version of a1 form factor, used in tauola
47 c. called : by form1-form3 in $c_cvssrc/korb/koralb/formf.f
48 c. author : alan weinstein 2/98
50 c. detailed description
51 c. first determine whether we are doing pi-2pi0 or 3pi.
52 c.
Then implement full form-factor from fit:
53 c. [(rho-pi s-wave) + (rho-prim-pi s-wave) +
54 c. (rho-pi d-wave) + (rho-prim-pi d-wave) +
55 c. (f2 pi d-wave) + (sigmapi s-wave) + (f0pi s-wave)]
56 c. based on fit to pi-2pi0 by m. schmidler, cbx 97-64-update(4/22/98)
57 c. all the parameters in this routine are hard-coded
59 c.......................................................................
60 * -------------------- argument declarations ---------------
64 * --------------------
EXTERNAL declarations ---------------
68 * -------------------- sequence declarations ---------------
70 * -------------------- local declarations ---------------
73 parameter( crname =
'F3PI' )
76 REAL mro,gro,mrp,grp,mf2,gf2,mf0,gf0,msg,gsg
77 REAL m1,m2,m3,m1sq,m2sq,m3sq,mpiz,mpic
79 REAL f134,f150,f15a,f15b,f167
80 REAL f34a,f34b,f35,f35a,f35b,f36a,f36b
81 COMPLEX bt1,bt2,bt3,bt4,bt5,bt6,bt7
82 COMPLEX fro1,fro2,frp1,frp2
83 COMPLEX ff21,ff22,ff23,fsg1,fsg2,fsg3,ff01,ff02,ff03
86 * --------------------
SAVE declarations ---------------
88 * --------------------
DATA initializations ---------------
91 * ----------------- executable code starts here ------------
93 c. hard-code the fit parameters:
96 c rho, rhoprime, f2(1275), f0(1186), sigma(made up
110 c fit coefficients for each of the contributions:
113 bt2 = cmplx(0.12,0.)*cexp(cmplx(0., 0.99*pi))
114 bt3 = cmplx(0.37,0.)*cexp(cmplx(0.,-0.15*pi))
115 bt4 = cmplx(0.87,0.)*cexp(cmplx(0., 0.53*pi))
116 bt5 = cmplx(0.71,0.)*cexp(cmplx(0., 0.56*pi))
117 bt6 = cmplx(2.10,0.)*cexp(cmplx(0., 0.23*pi))
118 bt7 = cmplx(0.77,0.)*cexp(cmplx(0.,-0.54*pi))
120 print *,
' In F3pi: add (rho-pi S-wave) + (rhop-pi S-wave) +'
121 print *,
' (rho-pi D-wave) + (rhop-pi D-wave) +'
122 print *,
' (f2 pi D-wave) + (sigmapi S-wave) + (f0pi S-wave)'
128 c. first determine whether we are doing pi-2pi0 or 3pi.
129 c pkorb is set up to remember what flavor of 3pi it gave to koralb,
130 c since koralb doesnt bother to remember
149 c.
Then implement full form-factor from fit:
150 c. [(rho-pi s-wave) + (rho-prim-pi s-wave) +
151 c. (rho-pi d-wave) + (rho-prim-pi d-wave) +
152 c. (f2 pi d-wave) + (sigmapi s-wave) + (f0pi s-wave)]
153 c. based on fit to pi-2pi0 by m. schmidler, cbx 97-64-update(4/22/98)
155 c note that for form1, the arguments are s1, s2;
156 c for form2, the arguments are s2, s1;
157 c for form3, the arguments are s3, s1.
158 c here, we implement form1 and form2 at the same time,
159 c so the above switch is just what we need
161 IF (iform.EQ.1.OR.iform.EQ.2)
THEN
164 s3 = qq-sa-sb+m1sq+m2sq+m3sq
165 IF (s3.LE.0..OR.s2.LE.0.)
RETURN
169 c lorentz invariants for all the contributions:
170 f134 = -(1./3.)*((s3-m3sq)-(s1-m1sq))
171 f150 = (1./18.)*(qq-m3sq+s3)*(2.*m1sq+2.*m2sq-s3)/s3
174 c breit wigners for all the contributions:
175 fro1 = bwigml(s1,mro,gro,m2,m3,1)
176 frp1 = bwigml(s1,mrp,grp,m2,m3,1)
177 fro2 = bwigml(s2,mro,gro,m3,m1,1)
178 frp2 = bwigml(s2,mrp,grp,m3,m1,1)
179 ff23 = bwigml(s3,mf2,gf2,m1,m2,2)
180 fsg3 = bwigml(s3,msg,gsg,m1,m2,0)
181 ff03 = bwigml(s3,mf0,gf0,m1,m2,0)
183 f3pi = bt1*fro1+bt2*frp1+
184 1 bt3*cmplx(f134,0.)*fro2+bt4*cmplx(f134,0.)*frp2+
185 1 bt5*cmplx(f150,0.)*ff23+
186 1 bt6*cmplx(f167,0.)*fsg3+bt7*cmplx(f167,0.)*ff03
188 c f3pi = fpikm(sqrt(s1),m2,m3)
189 ELSEIF (idk.EQ.2)
THEN
191 c lorentz invariants for all the contributions:
192 f134 = -(1./3.)*((s3-m3sq)-(s1-m1sq))
193 f15a = -(1./2.)*((s2-m2sq)-(s3-m3sq))
194 f15b = -(1./18.)*(qq-m2sq+s2)*(2.*m1sq+2.*m3sq-s2)/s2
197 c breit wigners for all the contributions:
198 fro1 = bwigml(s1,mro,gro,m2,m3,1)
199 frp1 = bwigml(s1,mrp,grp,m2,m3,1)
200 fro2 = bwigml(s2,mro,gro,m3,m1,1)
201 frp2 = bwigml(s2,mrp,grp,m3,m1,1)
202 ff21 = bwigml(s1,mf2,gf2,m2,m3,2)
203 ff22 = bwigml(s2,mf2,gf2,m3,m1,2)
204 fsg2 = bwigml(s2,msg,gsg,m3,m1,0)
205 ff02 = bwigml(s2,mf0,gf0,m3,m1,0)
207 f3pi = bt1*fro1+bt2*frp1+
208 1 bt3*cmplx(f134,0.)*fro2+bt4*cmplx(f134,0.)*frp2
209 1 -bt5*cmplx(f15a,0.)*ff21-bt5*cmplx(f15b,0.)*ff22
210 1 -bt6*cmplx(f167,0.)*fsg2-bt7*cmplx(f167,0.)*ff02
212 c f3pi = fpikm(sqrt(s1),m2,m3)
215 ELSE IF (iform.EQ.3)
THEN
218 s2 = qq-sa-sb+m1sq+m2sq+m3sq
219 IF (s1.LE.0..OR.s2.LE.0.)
RETURN
223 c lorentz invariants for all the contributions:
224 f34a = (1./3.)*((s2-m2sq)-(s3-m3sq))
225 f34b = (1./3.)*((s3-m3sq)-(s1-m1sq))
226 f35 =-(1./2.)*((s1-m1sq)-(s2-m2sq))
228 c breit wigners for all the contributions:
229 fro1 = bwigml(s1,mro,gro,m2,m3,1)
230 frp1 = bwigml(s1,mrp,grp,m2,m3,1)
231 fro2 = bwigml(s2,mro,gro,m3,m1,1)
232 frp2 = bwigml(s2,mrp,grp,m3,m1,1)
233 ff23 = bwigml(s3,mf2,gf2,m1,m2,2)
236 1 bt3*(cmplx(f34a,0.)*fro1+cmplx(f34b,0.)*fro2)+
237 1 bt4*(cmplx(f34a,0.)*frp1+cmplx(f34b,0.)*frp2)+
238 1 bt5*cmplx(f35,0.)*ff23
240 c f3pi = cmplx(0.,0.)
241 ELSEIF (idk.EQ.2)
THEN
243 c lorentz invariants for all the contributions:
244 f34a = (1./3.)*((s2-m2sq)-(s3-m3sq))
245 f34b = (1./3.)*((s3-m3sq)-(s1-m1sq))
246 f35a = -(1./18.)*(qq-m1sq+s1)*(2.*m2sq+2.*m3sq-s1)/s1
247 f35b = (1./18.)*(qq-m2sq+s2)*(2.*m3sq+2.*m1sq-s2)/s2
251 c breit wigners for all the contributions:
252 fro1 = bwigml(s1,mro,gro,m2,m3,1)
253 frp1 = bwigml(s1,mrp,grp,m2,m3,1)
254 fro2 = bwigml(s2,mro,gro,m3,m1,1)
255 frp2 = bwigml(s2,mrp,grp,m3,m1,1)
256 ff21 = bwigml(s1,mf2,gf2,m2,m3,2)
257 ff22 = bwigml(s2,mf2,gf2,m3,m1,2)
258 fsg1 = bwigml(s1,msg,gsg,m2,m3,0)
259 fsg2 = bwigml(s2,msg,gsg,m3,m1,0)
260 ff01 = bwigml(s1,mf0,gf0,m2,m3,0)
261 ff02 = bwigml(s2,mf0,gf0,m3,m1,0)
264 1 bt3*(cmplx(f34a,0.)*fro1+cmplx(f34b,0.)*fro2)+
265 1 bt4*(cmplx(f34a,0.)*frp1+cmplx(f34b,0.)*frp2)
266 1 -bt5*(cmplx(f35a,0.)*ff21+cmplx(f35b,0.)*ff22)
267 1 -bt6*(cmplx(f36a,0.)*fsg1+cmplx(f36b,0.)*fsg2)
268 1 -bt7*(cmplx(f36a,0.)*ff01+cmplx(f36b,0.)*ff02)
270 c f3pi = cmplx(0.,0.)
274 c add overall a1/a1prime:
280 c **********************************************************
281 COMPLEX FUNCTION bwigml(S,M,G,M1,M2,L)
282 c **********************************************************
283 c l-wave breit-wigner
284 c **********************************************************
287 REAL msq,w,wgs,mp,mm,qs,qm
294 IF (w.GT.(m1+m2))
THEN
295 qs=sqrt(abs((s -mp)*(s -mm)))/w
296 qm=sqrt(abs((msq -mp)*(msq -mm)))/m
298 wgs=g*(msq/w)*(qs/qm)**ipow
301 bwigml=cmplx(msq,0.)/cmplx(msq-s,-wgs)
305 c=======================================================================
306 COMPLEX FUNCTION fa1a1p(XMSQ)
307 c ==================================================================
308 c
complex form-factor for a1+a1prime. ajw 1/98
309 c ==================================================================
313 REAL xm1,xg1,xm2,xg2,xm1sq,xm2sq,gg1,gg2,gf,fg1,fg2
317 IF (ifirst.EQ.0)
THEN
320 c the user may choose masses and widths that differ from nominal:
325 bet = cmplx(pkorb(3,17),0.)
326 c scale factors relative to nominal:
327 gg1 = xm1*xg1/(1.3281*0.806)
328 gg2 = xm2*xg2/(1.3281*0.806)
337 f1 = cmplx(-xm1sq,0.0)/cmplx(xmsq-xm1sq,fg1)
338 f2 = cmplx(-xm2sq,0.0)/cmplx(xmsq-xm2sq,fg2)
343 c=======================================================================
346 c mass-dependent m*gamma of a1 through its decays to
347 c. [(rho-pi s-wave) + (rho-pi d-wave) +
348 c. (f2 pi d-wave) + (f0pi s-wave)]
349 c. and simple k*k s-wave
352 DOUBLE PRECISION mkst,mk,mk1sq,mk2sq,c3pi,ckst
353 DOUBLE PRECISION s,wga1c,wga1n,wg3pic,wg3pin,gkst
355 c-----------------------------------------------------------------------
357 IF (ifirst.NE.987)
THEN
360 c contribution to m*gamma(m(3pi)^2) from s-wave k*k:
365 c coupling constants squared:
367 ckst = 4.7621d0**2*c3pi
370 c-----------------------------------------------------------------------
371 c parameterization of numerical integral of total width of a1 to 3pi.
372 c from m. schmidtler, cbx-97-64-update.
377 c contribution to m*gamma(m(3pi)^2) from s-wave k*k,
if above threshold
379 IF (s.GT.mk1sq) gkst = sqrt((s-mk1sq)*(s-mk2sq))/(2.*s)
381 wga1 = sngl(c3pi*(wg3pic+wg3pin)+ckst*gkst)
385 c=======================================================================
386 DOUBLE PRECISION FUNCTION wga1c(S)
388 c parameterization of m*gamma(m^2) for pi-2pi0 system
390 DOUBLE PRECISION s,sth,q0,q1,q2,p0,p1,p2,p3,p4,g1_im
392 parameter(q0 = 5.80900d0,q1 = -3.00980d0,q2 = 4.57920d0,
393 1 p0 = -13.91400d0,p1 = 27.67900d0,p2 = -13.39300d0,
394 2 p3 = 3.19240d0,p4 = -0.10487d0)
396 parameter(sth = 0.1753d0)
397 c---------------------------------------------------------------------
401 ELSEIF((s.GT.sth).AND.(s.LT.0.823d0))
THEN
402 g1_im = q0*(s-sth)**3*(1. + q1*(s-sth) + q2*(s-sth)**2)
404 g1_im = p0 + p1*s + p2*s**2+ p3*s**3 + p4*s**4
410 c=======================================================================
411 DOUBLE PRECISION FUNCTION wga1n(S)
413 c parameterization of m*gamma(m^2) for pi-pi+pi- system
415 DOUBLE PRECISION s,sth,q0,q1,q2,p0,p1,p2,p3,p4,g1_im
417 parameter(q0 = 6.28450d0,q1 = -2.95950d0,q2 = 4.33550d0,
418 1 p0 = -15.41100d0,p1 = 32.08800d0,p2 = -17.66600d0,
419 2 p3 = 4.93550d0,p4 = -0.37498d0)
421 parameter(sth = 0.1676d0)
422 c---------------------------------------------------------------------
426 ELSEIF((s.GT.sth).AND.(s.LT.0.823d0))
THEN
427 g1_im = q0*(s-sth)**3*(1. + q1*(s-sth) + q2*(s-sth)**2)
429 g1_im = p0 + p1*s + p2*s**2+ p3*s**3 + p4*s**4