C++ Interface to Tauola
tauola-F/f3pi.F
1*AJW CLEO version of a1 form factor
2 COMPLEX FUNCTION f3pi(IFORM,QQ,SA,SB)
3C.......................................................................
4C.
5C. F3PI - CLEO version of a1 form factor, used in TAUOLA
6C.
7C. Inputs : None
8C. :
9C. Outputs : None
10C.
11C. COMMON : None
12C.
13C. Calls :
14C. Called : by FORM1-FORM3 in $C_CVSSRC/korb/koralb/formf.F
15C. Author : Alan Weinstein 2/98
16C.
17C. Detailed description
18C. First determine whether we are doing pi-2pi0 or 3pi.
19C. Then implement full form-factor from fit:
20C. [(rho-pi S-wave) + (rho-prim-pi S-wave) +
21C. (rho-pi D-wave) + (rho-prim-pi D-wave) +
22C. (f2 pi D-wave) + (sigmapi S-wave) + (f0pi S-wave)]
23C. based on fit to pi-2pi0 by M. Schmidler, CBX 97-64-Update (4/22/98)
24C. All the parameters in this routine are hard-coded!!
25C.
26C.......................................................................
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*
63C. Hard-code the fit parameters:
64 IF (ifirst.EQ.0) THEN
65 ifirst = 1
66C 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
80C 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
95C Initialize to 0:
96 f3pi = cmplx(0.,0.)
97
98C. First determine whether we are doing pi-2pi0 or 3pi.
99C PKORB is set up to remember what flavor of 3pi it gave to KORALB,
100C since KORALB doesnt bother to remember!!
101 r = pkorb(4,11)
102 IF (r.EQ.0.) THEN
103C it is 2pi0pi-
104 idk = 1
105 m1 = mpiz
106 m2 = mpiz
107 m3 = mpic
108 ELSE
109C 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
119C. Then implement full form-factor from fit:
120C. [(rho-pi S-wave) + (rho-prim-pi S-wave) +
121C. (rho-pi D-wave) + (rho-prim-pi D-wave) +
122C. (f2 pi D-wave) + (sigmapi S-wave) + (f0pi S-wave)]
123C. based on fit to pi-2pi0 by M. Schmidler, CBX 97-64-Update (4/22/98)
124
125C Note that for FORM1, the arguments are S1, S2;
126C for FORM2, the arguments are S2, S1;
127C for FORM3, the arguments are S3, S1.
128C Here, we implement FORM1 and FORM2 at the same time,
129C 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
138C it is 2pi0pi-
139C 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
144C 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
158C F3PI = FPIKM(SQRT(S1),M2,M3)
159 ELSEIF (idk.EQ.2) THEN
160C it is 3pi
161C 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
167C 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
182C 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
192C it is 2pi0pi-
193C 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
198C 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
210C F3PI = CMPLX(0.,0.)
211 ELSEIF (idk.EQ.2) THEN
212C it is 3pi
213C 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
221C 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
240C F3PI = CMPLX(0.,0.)
241 END IF
242 END IF
243
244C Add overall a1/a1prime:
245 forma1 = fa1a1p(qq)
246 f3pi = f3pi*forma1
247
248 RETURN
249 END
250C **********************************************************
251 COMPLEX FUNCTION bwigml(S,M,G,M1,M2,L)
252C **********************************************************
253C L-WAVE BREIT-WIGNER
254C **********************************************************
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
275C=======================================================================
276 COMPLEX FUNCTION fa1a1p(XMSQ)
277C ==================================================================
278C complex form-factor for a1+a1prime. AJW 1/98
279C ==================================================================
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
290C 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.)
296C 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
313C=======================================================================
314 FUNCTION wga1(QQ)
315
316C mass-dependent M*Gamma of a1 through its decays to
317C. [(rho-pi S-wave) + (rho-pi D-wave) +
318C. (f2 pi D-wave) + (f0pi S-wave)]
319C. 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
325C-----------------------------------------------------------------------
326C
327 IF (ifirst.NE.987) THEN
328 ifirst = 987
329C
330C 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
335C coupling constants squared:
336 c3pi = 0.2384d0**2
337 ckst = 4.7621d0**2*c3pi
338 END IF
339
340C-----------------------------------------------------------------------
341C Parameterization of numerical integral of total width of a1 to 3pi.
342C From M. Schmidtler, CBX-97-64-Update.
343 s = dble(qq)
344 wg3pic = wga1c(s)
345 wg3pin = wga1n(s)
346
347C 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
355C=======================================================================
356 DOUBLE PRECISION FUNCTION wga1c(S)
357C
358C parameterization of m*Gamma(m^2) for pi-2pi0 system
359C
360 DOUBLE PRECISION S,STH,Q0,Q1,Q2,P0,P1,P2,P3,P4,G1_IM
361C
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)
365C
366 parameter(sth = 0.1753d0)
367C---------------------------------------------------------------------
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
380C=======================================================================
381 DOUBLE PRECISION FUNCTION wga1n(S)
382C
383C parameterization of m*Gamma(m^2) for pi-pi+pi- system
384C
385 DOUBLE PRECISION S,STH,Q0,Q1,Q2,P0,P1,P2,P3,P4,G1_IM
386C
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)
390C
391 parameter(sth = 0.1676d0)
392C---------------------------------------------------------------------
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