C++InterfacetoTauola
BBB/prod/f3pi.f
1 *AJW 1 version of a1 form factor
2  COMPLEX FUNCTION f3pi(IFORM,QQ,SA,SB)
3 C.......................................................................
4 C.
5 C. F3PI - 1 version of a1 form factor, used in TAUOLA
6 C.
7 C. Inputs : None
8 C. :
9 C. Outputs : None
10 C.
11 C. COMMON : None
12 C.
13 C. Calls :
14 C. Called : by FORM1-FORM3 in $C_CVSSRC/korb/koralb/formf.F
15 C. Author : Alan Weinstein 2/98
16 C.
17 C. Detailed description
18 C. First determine whether we are doing pi-2pi0 or 3pi.
19 C. Then implement full form-factor from fit:
20 C. [(rho-pi S-wave) + (rho-prim-pi S-wave) +
21 C. (rho-pi D-wave) + (rho-prim-pi D-wave) +
22 C. (f2 pi D-wave) + (sigmapi S-wave) + (f0pi S-wave)]
23 C. based on fit to pi-2pi0 by M. Schmidler, CBX 97-64-Update (4/22/98)
24 C. All the parameters in this routine are hard-coded!!
25 C.
26 C.......................................................................
27 
28 * -------------------- Argument declarations ---------------
29 
30  INTEGER iform
31  REAL qq,sa,sb
32 * -------------------- EXTERNAL declarations ---------------
33 *
34  REAL pkorb
35  COMPLEX bwigml
36 * -------------------- SEQUENCE declarations ---------------
37 *
38 * -------------------- Local declarations ---------------
39 *
40  CHARACTER*(*) crname
41  parameter( crname = 'F3PI' )
42 *
43  INTEGER ifirst,idk
44  REAL mro,gro,mrp,grp,mf2,gf2,mf0,gf0,msg,gsg
45  REAL m1,m2,m3,m1sq,m2sq,m3sq,mpiz,mpic
46  REAL s1,s2,s3,r,pi
47  REAL f134,f150,f15a,f15b,f167
48  REAL f34a,f34b,f35,f35a,f35b,f36a,f36b
49  COMPLEX bt1,bt2,bt3,bt4,bt5,bt6,bt7
50  COMPLEX fro1,fro2,frp1,frp2
51  COMPLEX ff21,ff22,ff23,fsg1,fsg2,fsg3,ff01,ff02,ff03
52  COMPLEX fa1a1p,forma1
53 
54 * -------------------- SAVE declarations ---------------
55 *
56 * -------------------- DATA initializations ---------------
57 *
58  DATA ifirst/0/
59 * ----------------- Executable code starts here ------------
60 *
61 C. Hard-code the fit parameters:
62  IF (ifirst.EQ.0) THEN
63  ifirst = 1
64 C rho, rhoprime, f2(1275), f0(1186), sigma(made up!)
65  mro = 0.7743
66  gro = 0.1491
67  mrp = 1.370
68  grp = 0.386
69  mf2 = 1.275
70  gf2 = 0.185
71  mf0 = 1.186
72  gf0 = 0.350
73  msg = 0.860
74  gsg = 0.880
75  mpiz = pkorb(1,7)
76  mpic = pkorb(1,8)
77 
78 C Fit coefficients for each of the contributions:
79  pi = 3.14159
80  bt1 = cmplx(1.,0.)
81  bt2 = cmplx(0.12,0.)*cexp(cmplx(0., 0.99*pi))
82  bt3 = cmplx(0.37,0.)*cexp(cmplx(0.,-0.15*pi))
83  bt4 = cmplx(0.87,0.)*cexp(cmplx(0., 0.53*pi))
84  bt5 = cmplx(0.71,0.)*cexp(cmplx(0., 0.56*pi))
85  bt6 = cmplx(2.10,0.)*cexp(cmplx(0., 0.23*pi))
86  bt7 = cmplx(0.77,0.)*cexp(cmplx(0.,-0.54*pi))
87 
88  print *,' In F3pi: add (rho-pi S-wave) + (rhop-pi S-wave) +'
89  print *,' (rho-pi D-wave) + (rhop-pi D-wave) +'
90  print *,' (f2 pi D-wave) + (sigmapi S-wave) + (f0pi S-wave)'
91  END IF
92 
93 C Initialize to 0:
94  f3pi = cmplx(0.,0.)
95 
96 C. First determine whether we are doing pi-2pi0 or 3pi.
97 C PKORB is set up to remember what flavor of 3pi it gave to KORALB,
98 C since KORALB doesnt bother to remember!!
99  r = pkorb(4,11)
100  IF (r.EQ.0.) THEN
101 C it is 2pi0pi-
102  idk = 1
103  m1 = mpiz
104  m2 = mpiz
105  m3 = mpic
106  ELSE
107 C it is 3pi
108  idk = 2
109  m1 = mpic
110  m2 = mpic
111  m3 = mpic
112  END IF
113  m1sq = m1*m1
114  m2sq = m2*m2
115  m3sq = m3*m3
116 
117 C. Then implement full form-factor from fit:
118 C. [(rho-pi S-wave) + (rho-prim-pi S-wave) +
119 C. (rho-pi D-wave) + (rho-prim-pi D-wave) +
120 C. (f2 pi D-wave) + (sigmapi S-wave) + (f0pi S-wave)]
121 C. based on fit to pi-2pi0 by M. Schmidler, CBX 97-64-Update (4/22/98)
122 
123 C Note that for FORM1, the arguments are S1, S2;
124 C for FORM2, the arguments are S2, S1;
125 C for FORM3, the arguments are S3, S1.
126 C Here, we implement FORM1 and FORM2 at the same time,
127 C so the above switch is just what we need!
128 
129  IF (iform.EQ.1.OR.iform.EQ.2) THEN
130  s1 = sa
131  s2 = sb
132  s3 = qq-sa-sb+m1sq+m2sq+m3sq
133  IF (s3.LE.0..OR.s2.LE.0.) RETURN
134 
135  IF (idk.EQ.1) THEN
136 C it is 2pi0pi-
137 C Lorentz invariants for all the contributions:
138  f134 = -(1./3.)*((s3-m3sq)-(s1-m1sq))
139  f150 = (1./18.)*(qq-m3sq+s3)*(2.*m1sq+2.*m2sq-s3)/s3
140  f167 = (2./3.)
141 
142 C Breit Wigners for all the contributions:
143  fro1 = bwigml(s1,mro,gro,m2,m3,1)
144  frp1 = bwigml(s1,mrp,grp,m2,m3,1)
145  fro2 = bwigml(s2,mro,gro,m3,m1,1)
146  frp2 = bwigml(s2,mrp,grp,m3,m1,1)
147  ff23 = bwigml(s3,mf2,gf2,m1,m2,2)
148  fsg3 = bwigml(s3,msg,gsg,m1,m2,0)
149  ff03 = bwigml(s3,mf0,gf0,m1,m2,0)
150 
151  f3pi = bt1*fro1+bt2*frp1+
152  1 bt3*cmplx(f134,0.)*fro2+bt4*cmplx(f134,0.)*frp2+
153  1 bt5*cmplx(f150,0.)*ff23+
154  1 bt6*cmplx(f167,0.)*fsg3+bt7*cmplx(f167,0.)*ff03
155 
156 C F3PI = FPIKM(SQRT(S1),M2,M3)
157  ELSEIF (idk.EQ.2) THEN
158 C it is 3pi
159 C Lorentz invariants for all the contributions:
160  f134 = -(1./3.)*((s3-m3sq)-(s1-m1sq))
161  f15a = -(1./2.)*((s2-m2sq)-(s3-m3sq))
162  f15b = -(1./18.)*(qq-m2sq+s2)*(2.*m1sq+2.*m3sq-s2)/s2
163  f167 = -(2./3.)
164 
165 C Breit Wigners for all the contributions:
166  fro1 = bwigml(s1,mro,gro,m2,m3,1)
167  frp1 = bwigml(s1,mrp,grp,m2,m3,1)
168  fro2 = bwigml(s2,mro,gro,m3,m1,1)
169  frp2 = bwigml(s2,mrp,grp,m3,m1,1)
170  ff21 = bwigml(s1,mf2,gf2,m2,m3,2)
171  ff22 = bwigml(s2,mf2,gf2,m3,m1,2)
172  fsg2 = bwigml(s2,msg,gsg,m3,m1,0)
173  ff02 = bwigml(s2,mf0,gf0,m3,m1,0)
174 
175  f3pi = bt1*fro1+bt2*frp1+
176  1 bt3*cmplx(f134,0.)*fro2+bt4*cmplx(f134,0.)*frp2
177  1 -bt5*cmplx(f15a,0.)*ff21-bt5*cmplx(f15b,0.)*ff22
178  1 -bt6*cmplx(f167,0.)*fsg2-bt7*cmplx(f167,0.)*ff02
179 
180 C F3PI = FPIKM(SQRT(S1),M2,M3)
181  END IF
182 
183  ELSE IF (iform.EQ.3) THEN
184  s3 = sa
185  s1 = sb
186  s2 = qq-sa-sb+m1sq+m2sq+m3sq
187  IF (s1.LE.0..OR.s2.LE.0.) RETURN
188 
189  IF (idk.EQ.1) THEN
190 C it is 2pi0pi-
191 C Lorentz invariants for all the contributions:
192  f34a = (1./3.)*((s2-m2sq)-(s3-m3sq))
193  f34b = (1./3.)*((s3-m3sq)-(s1-m1sq))
194  f35 =-(1./2.)*((s1-m1sq)-(s2-m2sq))
195 
196 C Breit Wigners for all the contributions:
197  fro1 = bwigml(s1,mro,gro,m2,m3,1)
198  frp1 = bwigml(s1,mrp,grp,m2,m3,1)
199  fro2 = bwigml(s2,mro,gro,m3,m1,1)
200  frp2 = bwigml(s2,mrp,grp,m3,m1,1)
201  ff23 = bwigml(s3,mf2,gf2,m1,m2,2)
202 
203  f3pi =
204  1 bt3*(cmplx(f34a,0.)*fro1+cmplx(f34b,0.)*fro2)+
205  1 bt4*(cmplx(f34a,0.)*frp1+cmplx(f34b,0.)*frp2)+
206  1 bt5*cmplx(f35,0.)*ff23
207 
208 C F3PI = CMPLX(0.,0.)
209  ELSEIF (idk.EQ.2) THEN
210 C it is 3pi
211 C Lorentz invariants for all the contributions:
212  f34a = (1./3.)*((s2-m2sq)-(s3-m3sq))
213  f34b = (1./3.)*((s3-m3sq)-(s1-m1sq))
214  f35a = -(1./18.)*(qq-m1sq+s1)*(2.*m2sq+2.*m3sq-s1)/s1
215  f35b = (1./18.)*(qq-m2sq+s2)*(2.*m3sq+2.*m1sq-s2)/s2
216  f36a = -(2./3.)
217  f36b = (2./3.)
218 
219 C Breit Wigners for all the contributions:
220  fro1 = bwigml(s1,mro,gro,m2,m3,1)
221  frp1 = bwigml(s1,mrp,grp,m2,m3,1)
222  fro2 = bwigml(s2,mro,gro,m3,m1,1)
223  frp2 = bwigml(s2,mrp,grp,m3,m1,1)
224  ff21 = bwigml(s1,mf2,gf2,m2,m3,2)
225  ff22 = bwigml(s2,mf2,gf2,m3,m1,2)
226  fsg1 = bwigml(s1,msg,gsg,m2,m3,0)
227  fsg2 = bwigml(s2,msg,gsg,m3,m1,0)
228  ff01 = bwigml(s1,mf0,gf0,m2,m3,0)
229  ff02 = bwigml(s2,mf0,gf0,m3,m1,0)
230 
231  f3pi =
232  1 bt3*(cmplx(f34a,0.)*fro1+cmplx(f34b,0.)*fro2)+
233  1 bt4*(cmplx(f34a,0.)*frp1+cmplx(f34b,0.)*frp2)
234  1 -bt5*(cmplx(f35a,0.)*ff21+cmplx(f35b,0.)*ff22)
235  1 -bt6*(cmplx(f36a,0.)*fsg1+cmplx(f36b,0.)*fsg2)
236  1 -bt7*(cmplx(f36a,0.)*ff01+cmplx(f36b,0.)*ff02)
237 
238 C F3PI = CMPLX(0.,0.)
239  END IF
240  END IF
241 
242 C Add overall a1/a1prime:
243  forma1 = fa1a1p(qq)
244  f3pi = f3pi*forma1
245 
246  RETURN
247  END
248 C **********************************************************
249  COMPLEX FUNCTION bwigml(S,M,G,M1,M2,L)
250 C **********************************************************
251 C L-WAVE BREIT-WIGNER
252 C **********************************************************
253  REAL s,m,g,m1,m2
254  INTEGER l,ipow
255  REAL msq,w,wgs,mp,mm,qs,qm
256 
257  mp = (m1+m2)**2
258  mm = (m1-m2)**2
259  msq = m*m
260  w = sqrt(s)
261  wgs = 0.0
262  IF (w.GT.(m1+m2)) THEN
263  qs=sqrt(abs((s -mp)*(s -mm)))/w
264  qm=sqrt(abs((msq -mp)*(msq -mm)))/m
265  ipow = 2*l+1
266  wgs=g*(msq/w)*(qs/qm)**ipow
267  ENDIF
268 
269  bwigml=cmplx(msq,0.)/cmplx(msq-s,-wgs)
270 
271  RETURN
272  END
273 C=======================================================================
274  COMPLEX FUNCTION fa1a1p(XMSQ)
275 C ==================================================================
276 C complex form-factor for a1+a1prime. AJW 1/98
277 C ==================================================================
278 
279  REAL xmsq
280  REAL pkorb,wga1
281  REAL xm1,xg1,xm2,xg2,xm1sq,xm2sq,gg1,gg2,gf,fg1,fg2
282  COMPLEX bet,f1,f2
283  INTEGER ifirst/0/
284 
285  IF (ifirst.EQ.0) THEN
286  ifirst = 1
287 
288 C The user may choose masses and widths that differ from nominal:
289  xm1 = pkorb(1,10)
290  xg1 = pkorb(2,10)
291  xm2 = pkorb(1,17)
292  xg2 = pkorb(2,17)
293  bet = cmplx(pkorb(3,17),0.)
294 C scale factors relative to nominal:
295  gg1 = xm1*xg1/(1.3281*0.806)
296  gg2 = xm2*xg2/(1.3281*0.806)
297 
298  xm1sq = xm1*xm1
299  xm2sq = xm2*xm2
300  END IF
301 
302  gf = wga1(xmsq)
303  fg1 = gg1*gf
304  fg2 = gg2*gf
305  f1 = cmplx(-xm1sq,0.0)/cmplx(xmsq-xm1sq,fg1)
306  f2 = cmplx(-xm2sq,0.0)/cmplx(xmsq-xm2sq,fg2)
307  fa1a1p = f1+bet*f2
308 
309  RETURN
310  END
311 C=======================================================================
312  FUNCTION wga1(QQ)
313 
314 C mass-dependent M*Gamma of a1 through its decays to
315 C. [(rho-pi S-wave) + (rho-pi D-wave) +
316 C. (f2 pi D-wave) + (f0pi S-wave)]
317 C. AND simple K*K S-wave
318 
319  REAL qq,wga1
320  DOUBLE PRECISION mkst,mk,mk1sq,mk2sq,c3pi,ckst
321  DOUBLE PRECISION s,wga1c,wga1n,wg3pic,wg3pin,gkst
322  INTEGER ifirst
323 C-----------------------------------------------------------------------
324 C
325  IF (ifirst.NE.987) THEN
326  ifirst = 987
327 C
328 C Contribution to M*Gamma(m(3pi)^2) from S-wave K*K:
329  mkst = 0.894d0
330  mk = 0.496d0
331  mk1sq = (mkst+mk)**2
332  mk2sq = (mkst-mk)**2
333 C coupling constants squared:
334  c3pi = 0.2384d0**2
335  ckst = 4.7621d0**2*c3pi
336  END IF
337 
338 C-----------------------------------------------------------------------
339 C Parameterization of numerical integral of total width of a1 to 3pi.
340 C From M. Schmidtler, CBX-97-64-Update.
341  s = dble(qq)
342  wg3pic = wga1c(s)
343  wg3pin = wga1n(s)
344 
345 C Contribution to M*Gamma(m(3pi)^2) from S-wave K*K, if above threshold
346  gkst = 0.d0
347  IF (s.GT.mk1sq) gkst = sqrt((s-mk1sq)*(s-mk2sq))/(2.*s)
348 
349  wga1 = sngl(c3pi*(wg3pic+wg3pin)+ckst*gkst)
350 
351  RETURN
352  END
353 C=======================================================================
354  DOUBLE PRECISION FUNCTION wga1c(S)
355 C
356 C parameterization of m*Gamma(m^2) for pi-2pi0 system
357 C
358  DOUBLE PRECISION s,sth,q0,q1,q2,p0,p1,p2,p3,p4,g1_im
359 C
360  parameter(q0 = 5.80900d0,q1 = -3.00980d0,q2 = 4.57920d0,
361  1 p0 = -13.91400d0,p1 = 27.67900d0,p2 = -13.39300d0,
362  2 p3 = 3.19240d0,p4 = -0.10487d0)
363 C
364  parameter(sth = 0.1753d0)
365 C---------------------------------------------------------------------
366 
367  IF(s.LT.sth) THEN
368  g1_im = 0.d0
369  ELSEIF((s.GT.sth).AND.(s.LT.0.823d0)) THEN
370  g1_im = q0*(s-sth)**3*(1. + q1*(s-sth) + q2*(s-sth)**2)
371  ELSE
372  g1_im = p0 + p1*s + p2*s**2+ p3*s**3 + p4*s**4
373  ENDIF
374 
375  wga1c = g1_im
376  RETURN
377  END
378 C=======================================================================
379  DOUBLE PRECISION FUNCTION wga1n(S)
380 C
381 C parameterization of m*Gamma(m^2) for pi-pi+pi- system
382 C
383  DOUBLE PRECISION s,sth,q0,q1,q2,p0,p1,p2,p3,p4,g1_im
384 C
385  parameter(q0 = 6.28450d0,q1 = -2.95950d0,q2 = 4.33550d0,
386  1 p0 = -15.41100d0,p1 = 32.08800d0,p2 = -17.66600d0,
387  2 p3 = 4.93550d0,p4 = -0.37498d0)
388 C
389  parameter(sth = 0.1676d0)
390 C---------------------------------------------------------------------
391 
392  IF(s.LT.sth) THEN
393  g1_im = 0.d0
394  ELSEIF((s.GT.sth).AND.(s.LT.0.823d0)) THEN
395  g1_im = q0*(s-sth)**3*(1. + q1*(s-sth) + q2*(s-sth)**2)
396  ELSE
397  g1_im = p0 + p1*s + p2*s**2+ p3*s**3 + p4*s**4
398  ENDIF
399 
400  wga1n = g1_im
401  RETURN
402  END