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