C++InterfacetoTauola
modified/f3pi.f
1 
2 
3 *AJW 1 version of a1 form factor
4  COMPLEX FUNCTION f3pi(IFORM,QQ,SA,SB)
5 C.......................................................................
6 C.
7 C. F3PI - 1 version of a1 form factor, used in TAUOLA
8 C.
9 C. Inputs : None
10 C. :
11 C. Outputs : None
12 C.
13 C. COMMON : None
14 C.
15 C. Calls :
16 C. Called : by FORM1-FORM3 in $C_CVSSRC/korb/koralb/formf.F
17 C. Author : Alan Weinstein 2/98
18 C.
19 C. Detailed description
20 C. First determine whether we are doing pi-2pi0 or 3pi.
21 C. Then implement full form-factor from fit:
22 C. [(rho-pi S-wave) + (rho-prim-pi S-wave) +
23 C. (rho-pi D-wave) + (rho-prim-pi D-wave) +
24 C. (f2 pi D-wave) + (sigmapi S-wave) + (f0pi S-wave)]
25 C. based on fit to pi-2pi0 by M. Schmidler, CBX 97-64-Update (4/22/98)
26 C. All the parameters in this routine are hard-coded!!
27 C.
28 C.......................................................................
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 *
65 C. Hard-code the fit parameters:
66  IF (ifirst.EQ.0) THEN
67  ifirst = 1
68 C 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 
82 C 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 
97 C Initialize to 0:
98  f3pi = cmplx(0.,0.)
99 
100 C. First determine whether we are doing pi-2pi0 or 3pi.
101 C PKORB is set up to remember what flavor of 3pi it gave to KORALB,
102 C since KORALB doesnt bother to remember!!
103  r = pkorb(4,11)
104  IF (r.EQ.0.) THEN
105 C it is 2pi0pi-
106  idk = 1
107  m1 = mpiz
108  m2 = mpiz
109  m3 = mpic
110  ELSE
111 C 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 
121 C. Then implement full form-factor from fit:
122 C. [(rho-pi S-wave) + (rho-prim-pi S-wave) +
123 C. (rho-pi D-wave) + (rho-prim-pi D-wave) +
124 C. (f2 pi D-wave) + (sigmapi S-wave) + (f0pi S-wave)]
125 C. based on fit to pi-2pi0 by M. Schmidler, CBX 97-64-Update (4/22/98)
126 
127 C Note that for FORM1, the arguments are S1, S2;
128 C for FORM2, the arguments are S2, S1;
129 C for FORM3, the arguments are S3, S1.
130 C Here, we implement FORM1 and FORM2 at the same time,
131 C 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
140 C it is 2pi0pi-
141 C 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 
146 C 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 
160 C F3PI = FPIKM(SQRT(S1),M2,M3)
161  ELSEIF (idk.EQ.2) THEN
162 C it is 3pi
163 C 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 
169 C 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 
184 C 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
194 C it is 2pi0pi-
195 C 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 
200 C 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 
212 C F3PI = CMPLX(0.,0.)
213  ELSEIF (idk.EQ.2) THEN
214 C it is 3pi
215 C 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 
223 C 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 
242 C F3PI = CMPLX(0.,0.)
243  END IF
244  END IF
245 
246 C Add overall a1/a1prime:
247  forma1 = fa1a1p(qq)
248  f3pi = f3pi*forma1
249 
250  RETURN
251  END
252 C **********************************************************
253  COMPLEX FUNCTION bwigml(S,M,G,M1,M2,L)
254 C **********************************************************
255 C L-WAVE BREIT-WIGNER
256 C **********************************************************
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
277 C=======================================================================
278  COMPLEX FUNCTION fa1a1p(XMSQ)
279 C ==================================================================
280 C complex form-factor for a1+a1prime. AJW 1/98
281 C ==================================================================
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 
292 C 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.)
298 C 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
315 C=======================================================================
316  FUNCTION wga1(QQ)
317 
318 C mass-dependent M*Gamma of a1 through its decays to
319 C. [(rho-pi S-wave) + (rho-pi D-wave) +
320 C. (f2 pi D-wave) + (f0pi S-wave)]
321 C. 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
327 C-----------------------------------------------------------------------
328 C
329  IF (ifirst.NE.987) THEN
330  ifirst = 987
331 C
332 C 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
337 C coupling constants squared:
338  c3pi = 0.2384d0**2
339  ckst = 4.7621d0**2*c3pi
340  END IF
341 
342 C-----------------------------------------------------------------------
343 C Parameterization of numerical integral of total width of a1 to 3pi.
344 C From M. Schmidtler, CBX-97-64-Update.
345  s = dble(qq)
346  wg3pic = wga1c(s)
347  wg3pin = wga1n(s)
348 
349 C 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
357 C=======================================================================
358  DOUBLE PRECISION FUNCTION wga1c(S)
359 C
360 C parameterization of m*Gamma(m^2) for pi-2pi0 system
361 C
362  DOUBLE PRECISION s,sth,q0,q1,q2,p0,p1,p2,p3,p4,g1_im
363 C
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)
367 C
368  parameter(sth = 0.1753d0)
369 C---------------------------------------------------------------------
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
382 C=======================================================================
383  DOUBLE PRECISION FUNCTION wga1n(S)
384 C
385 C parameterization of m*Gamma(m^2) for pi-pi+pi- system
386 C
387  DOUBLE PRECISION s,sth,q0,q1,q2,p0,p1,p2,p3,p4,g1_im
388 C
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)
392 C
393  parameter(sth = 0.1676d0)
394 C---------------------------------------------------------------------
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
Generated on Tue Mar 12 2024 22:57:50 for C++InterfacetoTauola by   doxygen 1.8.7