C++InterfacetoTauola
F/prod/f3pi.f
1 /* copyright(c) 1991-2012 free software foundation, inc.
2  this file is part of the gnu c library.
3 
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.
8 
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.
13 
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/>. */
17 
18 
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. */
26 
27 /* we do support the iec 559 math functionality, real and complex. */
28 
29 /* wchar_t uses iso/iec 10646 (2nd ed., published 2011-03-15) /
30  unicode 6.0. */
31 
32 /* we do not support c11 <threads.h>. */
33 
34 *ajw 1 version of a1 form factor
35  COMPLEX FUNCTION f3pi(IFORM,QQ,SA,SB)
36 c.......................................................................
37 c.
38 c. f3pi - 1 version of a1 form factor, used in tauola
39 c.
40 c. inputs : none
41 c. :
42 c. outputs : none
43 c.
44 c. COMMON : none
45 c.
46 c. calls :
47 c. called : by form1-form3 in $c_cvssrc/korb/koralb/formf.f
48 c. author : alan weinstein 2/98
49 c.
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!!
58 c.
59 c.......................................................................
60 * -------------------- argument declarations ---------------
61 
62  INTEGER iform
63  REAL qq,sa,sb
64 * -------------------- EXTERNAL declarations ---------------
65 *
66  REAL pkorb
67  COMPLEX bwigml
68 * -------------------- sequence declarations ---------------
69 *
70 * -------------------- local declarations ---------------
71 *
72  CHARACTER*(*) crname
73  parameter( crname = 'F3PI' )
74 *
75  INTEGER ifirst,idk
76  REAL mro,gro,mrp,grp,mf2,gf2,mf0,gf0,msg,gsg
77  REAL m1,m2,m3,m1sq,m2sq,m3sq,mpiz,mpic
78  REAL s1,s2,s3,r,pi
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
84  COMPLEX fa1a1p,forma1
85 
86 * -------------------- SAVE declarations ---------------
87 *
88 * -------------------- DATA initializations ---------------
89 *
90  DATA ifirst/0/
91 * ----------------- executable code starts here ------------
92 *
93 c. hard-code the fit parameters:
94  IF (ifirst.EQ.0) THEN
95  ifirst = 1
96 c rho, rhoprime, f2(1275), f0(1186), sigma(made up!)
97  mro = 0.7743
98  gro = 0.1491
99  mrp = 1.370
100  grp = 0.386
101  mf2 = 1.275
102  gf2 = 0.185
103  mf0 = 1.186
104  gf0 = 0.350
105  msg = 0.860
106  gsg = 0.880
107  mpiz = pkorb(1,7)
108  mpic = pkorb(1,8)
109 
110 c fit coefficients for each of the contributions:
111  pi = 3.14159
112  bt1 = cmplx(1.,0.)
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))
119 
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)'
123  END IF
124 
125 c initialize to 0:
126  f3pi = cmplx(0.,0.)
127 
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!!
131  r = pkorb(4,11)
132  IF (r.EQ.0.) THEN
133 c it is 2pi0pi-
134  idk = 1
135  m1 = mpiz
136  m2 = mpiz
137  m3 = mpic
138  ELSE
139 c it is 3pi
140  idk = 2
141  m1 = mpic
142  m2 = mpic
143  m3 = mpic
144  END IF
145  m1sq = m1*m1
146  m2sq = m2*m2
147  m3sq = m3*m3
148 
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)
154 
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!
160 
161  IF (iform.EQ.1.OR.iform.EQ.2) THEN
162  s1 = sa
163  s2 = sb
164  s3 = qq-sa-sb+m1sq+m2sq+m3sq
165  IF (s3.LE.0..OR.s2.LE.0.) RETURN
166 
167  IF (idk.EQ.1) THEN
168 c it is 2pi0pi-
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
172  f167 = (2./3.)
173 
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)
182 
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
187 
188 c f3pi = fpikm(sqrt(s1),m2,m3)
189  ELSEIF (idk.EQ.2) THEN
190 c it is 3pi
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
195  f167 = -(2./3.)
196 
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)
206 
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
211 
212 c f3pi = fpikm(sqrt(s1),m2,m3)
213  END IF
214 
215  ELSE IF (iform.EQ.3) THEN
216  s3 = sa
217  s1 = sb
218  s2 = qq-sa-sb+m1sq+m2sq+m3sq
219  IF (s1.LE.0..OR.s2.LE.0.) RETURN
220 
221  IF (idk.EQ.1) THEN
222 c it is 2pi0pi-
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))
227 
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)
234 
235  f3pi =
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
239 
240 c f3pi = cmplx(0.,0.)
241  ELSEIF (idk.EQ.2) THEN
242 c it is 3pi
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
248  f36a = -(2./3.)
249  f36b = (2./3.)
250 
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)
262 
263  f3pi =
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)
269 
270 c f3pi = cmplx(0.,0.)
271  END IF
272  END IF
273 
274 c add overall a1/a1prime:
275  forma1 = fa1a1p(qq)
276  f3pi = f3pi*forma1
277 
278  RETURN
279  END
280 c **********************************************************
281  COMPLEX FUNCTION bwigml(S,M,G,M1,M2,L)
282 c **********************************************************
283 c l-wave breit-wigner
284 c **********************************************************
285  REAL s,m,g,m1,m2
286  INTEGER l,ipow
287  REAL msq,w,wgs,mp,mm,qs,qm
288 
289  mp = (m1+m2)**2
290  mm = (m1-m2)**2
291  msq = m*m
292  w = sqrt(s)
293  wgs = 0.0
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
297  ipow = 2*l+1
298  wgs=g*(msq/w)*(qs/qm)**ipow
299  ENDIF
300 
301  bwigml=cmplx(msq,0.)/cmplx(msq-s,-wgs)
302 
303  RETURN
304  END
305 c=======================================================================
306  COMPLEX FUNCTION fa1a1p(XMSQ)
307 c ==================================================================
308 c complex form-factor for a1+a1prime. ajw 1/98
309 c ==================================================================
310 
311  REAL xmsq
312  REAL pkorb,wga1
313  REAL xm1,xg1,xm2,xg2,xm1sq,xm2sq,gg1,gg2,gf,fg1,fg2
314  COMPLEX bet,f1,f2
315  INTEGER ifirst/0/
316 
317  IF (ifirst.EQ.0) THEN
318  ifirst = 1
319 
320 c the user may choose masses and widths that differ from nominal:
321  xm1 = pkorb(1,10)
322  xg1 = pkorb(2,10)
323  xm2 = pkorb(1,17)
324  xg2 = pkorb(2,17)
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)
329 
330  xm1sq = xm1*xm1
331  xm2sq = xm2*xm2
332  END IF
333 
334  gf = wga1(xmsq)
335  fg1 = gg1*gf
336  fg2 = gg2*gf
337  f1 = cmplx(-xm1sq,0.0)/cmplx(xmsq-xm1sq,fg1)
338  f2 = cmplx(-xm2sq,0.0)/cmplx(xmsq-xm2sq,fg2)
339  fa1a1p = f1+bet*f2
340 
341  RETURN
342  END
343 c=======================================================================
344  FUNCTION wga1(QQ)
345 
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
350 
351  REAL qq,wga1
352  DOUBLE PRECISION mkst,mk,mk1sq,mk2sq,c3pi,ckst
353  DOUBLE PRECISION s,wga1c,wga1n,wg3pic,wg3pin,gkst
354  INTEGER ifirst
355 c-----------------------------------------------------------------------
356 c
357  IF (ifirst.NE.987) THEN
358  ifirst = 987
359 c
360 c contribution to m*gamma(m(3pi)^2) from s-wave k*k:
361  mkst = 0.894d0
362  mk = 0.496d0
363  mk1sq = (mkst+mk)**2
364  mk2sq = (mkst-mk)**2
365 c coupling constants squared:
366  c3pi = 0.2384d0**2
367  ckst = 4.7621d0**2*c3pi
368  END IF
369 
370 c-----------------------------------------------------------------------
371 c parameterization of numerical integral of total width of a1 to 3pi.
372 c from m. schmidtler, cbx-97-64-update.
373  s = dble(qq)
374  wg3pic = wga1c(s)
375  wg3pin = wga1n(s)
376 
377 c contribution to m*gamma(m(3pi)^2) from s-wave k*k, if above threshold
378  gkst = 0.d0
379  IF (s.GT.mk1sq) gkst = sqrt((s-mk1sq)*(s-mk2sq))/(2.*s)
380 
381  wga1 = sngl(c3pi*(wg3pic+wg3pin)+ckst*gkst)
382 
383  RETURN
384  END
385 c=======================================================================
386  DOUBLE PRECISION FUNCTION wga1c(S)
387 c
388 c parameterization of m*gamma(m^2) for pi-2pi0 system
389 c
390  DOUBLE PRECISION s,sth,q0,q1,q2,p0,p1,p2,p3,p4,g1_im
391 c
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)
395 c
396  parameter(sth = 0.1753d0)
397 c---------------------------------------------------------------------
398 
399  IF(s.LT.sth) THEN
400  g1_im = 0.d0
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)
403  ELSE
404  g1_im = p0 + p1*s + p2*s**2+ p3*s**3 + p4*s**4
405  ENDIF
406 
407  wga1c = g1_im
408  RETURN
409  END
410 c=======================================================================
411  DOUBLE PRECISION FUNCTION wga1n(S)
412 c
413 c parameterization of m*gamma(m^2) for pi-pi+pi- system
414 c
415  DOUBLE PRECISION s,sth,q0,q1,q2,p0,p1,p2,p3,p4,g1_im
416 c
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)
420 c
421  parameter(sth = 0.1676d0)
422 c---------------------------------------------------------------------
423 
424  IF(s.LT.sth) THEN
425  g1_im = 0.d0
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)
428  ELSE
429  g1_im = p0 + p1*s + p2*s**2+ p3*s**3 + p4*s**4
430  ENDIF
431 
432  wga1n = g1_im
433  RETURN
434  END