C++ Interface to Tauola
tauola-BBB/prod/f3pi.f
1*AJW 1 version of a1 form factor
2 COMPLEX FUNCTION f3pi(IFORM,QQ,SA,SB)
3C.......................................................................
4C.
5C. F3PI - 1 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
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*
61C. Hard-code the fit parameters:
62 IF (ifirst.EQ.0) THEN
63 ifirst = 1
64C 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
78C 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
93C Initialize to 0:
94 f3pi = cmplx(0.,0.)
95
96C. First determine whether we are doing pi-2pi0 or 3pi.
97C PKORB is set up to remember what flavor of 3pi it gave to KORALB,
98C since KORALB doesnt bother to remember!!
99 r = pkorb(4,11)
100 IF (r.EQ.0.) THEN
101C it is 2pi0pi-
102 idk = 1
103 m1 = mpiz
104 m2 = mpiz
105 m3 = mpic
106 ELSE
107C 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
117C. Then implement full form-factor from fit:
118C. [(rho-pi S-wave) + (rho-prim-pi S-wave) +
119C. (rho-pi D-wave) + (rho-prim-pi D-wave) +
120C. (f2 pi D-wave) + (sigmapi S-wave) + (f0pi S-wave)]
121C. based on fit to pi-2pi0 by M. Schmidler, CBX 97-64-Update (4/22/98)
122
123C Note that for FORM1, the arguments are S1, S2;
124C for FORM2, the arguments are S2, S1;
125C for FORM3, the arguments are S3, S1.
126C Here, we implement FORM1 and FORM2 at the same time,
127C 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
136C it is 2pi0pi-
137C 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
142C 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
156C F3PI = FPIKM(SQRT(S1),M2,M3)
157 ELSEIF (idk.EQ.2) THEN
158C it is 3pi
159C 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
165C 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
180C 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
190C it is 2pi0pi-
191C 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
196C 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
208C F3PI = CMPLX(0.,0.)
209 ELSEIF (idk.EQ.2) THEN
210C it is 3pi
211C 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
219C 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
238C F3PI = CMPLX(0.,0.)
239 END IF
240 END IF
241
242C Add overall a1/a1prime:
243 forma1 = fa1a1p(qq)
244 f3pi = f3pi*forma1
245
246 RETURN
247 END
248C **********************************************************
249 COMPLEX FUNCTION bwigml(S,M,G,M1,M2,L)
250C **********************************************************
251C L-WAVE BREIT-WIGNER
252C **********************************************************
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
273C=======================================================================
274 COMPLEX FUNCTION fa1a1p(XMSQ)
275C ==================================================================
276C complex form-factor for a1+a1prime. AJW 1/98
277C ==================================================================
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
288C 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.)
294C 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
311C=======================================================================
312 FUNCTION wga1(QQ)
313
314C mass-dependent M*Gamma of a1 through its decays to
315C. [(rho-pi S-wave) + (rho-pi D-wave) +
316C. (f2 pi D-wave) + (f0pi S-wave)]
317C. 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
323C-----------------------------------------------------------------------
324C
325 IF (ifirst.NE.987) THEN
326 ifirst = 987
327C
328C 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
333C coupling constants squared:
334 c3pi = 0.2384d0**2
335 ckst = 4.7621d0**2*c3pi
336 END IF
337
338C-----------------------------------------------------------------------
339C Parameterization of numerical integral of total width of a1 to 3pi.
340C From M. Schmidtler, CBX-97-64-Update.
341 s = dble(qq)
342 wg3pic = wga1c(s)
343 wg3pin = wga1n(s)
344
345C 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
353C=======================================================================
354 DOUBLE PRECISION FUNCTION wga1c(S)
355C
356C parameterization of m*Gamma(m^2) for pi-2pi0 system
357C
358 DOUBLE PRECISION S,STH,Q0,Q1,Q2,P0,P1,P2,P3,P4,G1_IM
359C
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)
363C
364 parameter(sth = 0.1753d0)
365C---------------------------------------------------------------------
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
378C=======================================================================
379 DOUBLE PRECISION FUNCTION wga1n(S)
380C
381C parameterization of m*Gamma(m^2) for pi-pi+pi- system
382C
383 DOUBLE PRECISION S,STH,Q0,Q1,Q2,P0,P1,P2,P3,P4,G1_IM
384C
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)
388C
389 parameter(sth = 0.1676d0)
390C---------------------------------------------------------------------
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