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