C++InterfacetoTauola
funct_rpt.f
1 
2 C **********************************************************
3 C library of functions used in calculation of currents
4 C references:
5 C [1] arXiv:0911.4436 (hep-ph) D. Gomez Dumm et al. (tau -> 3pi nu)
6 C [2] arXiv:0911.2640 (hep-ph) D. Gomezz Dumm et al. (tau -> KKpi nu)
7 C [3] arXiv:0807.4883 (hep-ph) D. R. Boito et al. (tau -> Kpi nu)
8 C [4] P. Roig, talk at (tau -> 2 pi nu)
9 C [5] arXiv:0803.2039 (hep-ph) E. Arganda et al., Appendix B (tau -> KK nu)
10 C **********************************************************
11  FUNCTION grho_rcht(XS,XMMM)
12  IMPLICIT NONE
13  REAL xs
14  DOUBLE PRECISION xmmm
15 C **********************************************************
16 C REAL FUNCTION Gamma Rho ; energy-dependent width of rho meson
17 c in SU(2) limit mpi=mpi0, mk=mk0;
18 C formula (14) of REF [1]
19 C **********************************************************
20  include '../parameter.inc'
21  include '../funct_declar.inc'
22 
23  REAL mmpi_av2,mmk_2
24 
25  mmpi_av2 = mmpi_av**2
26  mmk_2 = mmk**2
27 
28  IF(xs.GE.(4.*mmk_2)) THEN
29  grho_rcht=xmmm*xs*((1.-4.*mmpi_av2/xs)**1.5
30  $ +0.5*(1.-4.*mmk_2/xs)**1.5)
31  $ /(96.*pi*fpi_rpt**2)
32  ELSE IF((xs.GE.(4.*mmpi_av2)).AND.(xs.LE.(4.*mmk_2))) THEN
33  grho_rcht=xmmm*xs*(1.-4.*mmpi_av2/xs)**1.5
34  $ /(96.*pi*fpi_rpt**2)
35  ELSE
36  grho_rcht = 0.
37  ENDIF
38 
39  RETURN
40 
41  END
42 
43 
44  FUNCTION grho1_rcht(XS,XMMM,XGGG)
45  IMPLICIT NONE
46  REAL xs
47  DOUBLE PRECISION xmmm,xggg
48 C **********************************************************
49 C REAL FUNCTION Gamma Rho1 ; energy-dependent width of rho1 meson;
50 c in SU(2) limit mpi=mpi0; only rho' -> 2pi loop;
51 C formula (33) of REF [1]
52 C **********************************************************
53  include '../parameter.inc'
54  include '../funct_declar.inc'
55 
56  REAL mmpi_av2
57  mmpi_av2 = mmpi_av**2
58 
59  IF (xs.GE.(4.*mmpi_av2)) THEN
60  grho1_rcht=xggg*sqrt(xmmm**2/xs)*
61  & ((xs-4.*mmpi_av2)/(xmmm**2-4.*mmpi_av2))**1.5
62  ELSE
63  grho1_rcht = 0.
64  ENDIF
65 
66  RETURN
67  END
68 
69 
70 
71  FUNCTION sigp(SS)
72  IMPLICIT NONE
73  DOUBLE PRECISION ss
74 C***********************************************************
75 C DOUBLRE PRECISION FUNCTION
76 c of two body phase space threshold: equal mass scalars
77 C***********************************************************
78  REAL tt
79  include '../parameter.inc'
80  include '../funct_declar.inc'
81 
82  tt = 1. - 4.*mmpi_av**2/ss
83 
84  IF (tt.GE.0) THEN
85  sigp = sqrt(tt)
86  ELSE
87  sigp = 0.
88  ENDIF
89 
90  RETURN
91  END
92 
93 
94 
95 
96  FUNCTION lamb_rcht(X1,X2,X3)
97  IMPLICIT NONE
98  REAL x1,x2,x3,arg_rcht
99 C***********************************************************
100 C REAL FUNCTION LAMDA of three body phase space
101 C***********************************************************
102  include '../funct_declar.inc'
103 
104  arg_rcht = (x1-x2-x3)**2 - 4.*x3*x2
105  IF(arg_rcht.GE.0.) THEN
106  lamb_rcht = arg_rcht
107  ELSE
108  lamb_rcht = 0.
109  ENDIF
110 
111  RETURN
112  END
113 
114 
115 
116 C******************** FUNCTIONS FOR 2 SCALAR MODES ***********************
117 
118 
119 
120 C****************************************************************
121  FUNCTION r0scal_3pi(QX,SX)
122 C****************************************************************
123 C Complex Function: R_0 function (Pablo notes)
124 C for the scalar contribution for three pion modes
125 C Called by f3pi_rcht.f
126 C **********************************************************
127  IMPLICIT NONE
128  REAL qx,sx,delta0_3piscal,xsx,xst
129  DOUBLE PRECISION xm1,xm2,dsx
130  include '../funct_declar.inc'
131  include '../parameter.inc'
132 c$$$ a00_3piscal = 0.220
133 c$$$ b00_3piscal = 0.268/mmpi_av**2
134 c$$$ c00_3piscal = -0.0139/mmpi_av**4
135 c$$$ d00_3piscal = -0.00139/mmpi_av**6
136 c$$$ x00_3piscal = 36.77*mmpi_av**2
137 c$$$c MMF0 = 0.98
138 c$$$
139  dsx = sx
140  xsx = sx/4.*sigp(dsx)**2
141  xst = sqrt(sx)
142 
143  if(sx.le.0.7) then
144  delta0_3piscal = sigp(dsx)*(a00_3piscal + b00_3piscal*xsx +
145  & c00_3piscal*xsx**2 + d00_3piscal*xsx**3)
146 
147  delta0_3piscal = delta0_3piscal*
148  & (4*mmpi_av**2 - x00_3piscal)/(sx-x00_3piscal)
149  else if (xst.le.1.21) then
150  delta0_3piscal = -10572.0+50658.0*xst-87903.0*xst**2+66886.0*xst**3
151  & -18699.0*xst**4
152  delta0_3piscal = delta0_3piscal*pi/180.0
153  else
154  delta0_3piscal = 255.0*pi/180.0
155  endif
156 
157 
158  delta0_3piscal = atan(delta0_3piscal)
159 
160  r0scal_3pi = alpha0_3pi/qx +
161  & alpha1_3pi/qx**2*(sx - mmf0**2)
162 
163  r0scal_3pi = r0scal_3pi*(cos(delta0_3piscal) +
164  & i*sin(delta0_3piscal))
165 
166  RETURN
167  END
168 
169 C****************************************************************
170  FUNCTION r2scal_3pi(QX,SX)
171 C****************************************************************
172 C Complex Function: R_2 function (Pablo notes)
173 C for the scalar contribution for three pion modes
174 C Called by f3pi_rcht.f
175 C **********************************************************
176  IMPLICIT NONE
177  REAL qx,sx,delta2_3piscal,xsx,xst
178  DOUBLE PRECISION xm1,xm2,dsx
179  include '../funct_declar.inc'
180  include '../parameter.inc'
181 c$$$ a02_3piscal = -0.0444
182 c$$$ b02_3piscal = -0.0857/mmpi_av**2
183 c$$$ c02_3piscal = -0.00221/mmpi_av**4
184 c$$$ d02_3piscal = -0.000129/mmpi_av**6
185 c$$$ x02_3piscal = -21.62*mmpi_av**2
186 c$$$c MMF0 = 0.98
187 
188  dsx = sx
189  xsx = sx/4.*sigp(dsx)**2
190  xst = sqrt(sx)
191 
192  if(sx.le.0.7) then
193  delta2_3piscal = sigp(dsx)*(a02_3piscal + b02_3piscal*xsx +
194  & c02_3piscal*xsx**2 + d02_3piscal*xsx**3)
195 
196  delta2_3piscal = delta2_3piscal*
197  & (4*mmpi_av**2 - x02_3piscal)/(sx-x02_3piscal)
198  else if(xst.le.1.21) then
199  delta2_3piscal = 282.9-1314.9*xst+2153.4*xst**2-1574.5d0*xst**3+
200  & 428.06d0*xst**4
201  delta2_3piscal = delta2_3piscal*pi/180.0
202  else
203  delta2_3piscal = -27.0*pi/180.0
204  endif
205 
206  delta2_3piscal = atan(delta2_3piscal)
207 
208  r2scal_3pi = gamma0_3pi/qx +
209  & gamma1_3pi/qx**2*(sx - mmf0**2)
210 
211  r2scal_3pi = r2scal_3pi*(cos(delta2_3piscal) +
212  & i*sin(delta2_3piscal))
213 
214  RETURN
215  END
216 
217 
218 
219 
220 C*************************************************************************
221 C Functions for sigma contributions
222 C*************************************************************************
223  FUNCTION ffsig(QX,XX)
224 C **********************************************************
225 C Complex Function:
226 C Called by f3pi_rcht.f
227 C **********************************************************
228  IMPLICIT NONE
229  REAL qx,xx,mm2
230  DOUBLE PRECISION xm1,xm2,xphi
231  include '../funct_declar.inc'
232  include '../parameter.inc'
233 
234  mm2 = mmpi_av**2
235 
236  xphi = - rsigma**2* lamb_rcht(qx,xx,mm2)/(8.*qx)
237 
238  ffsig = dexp(xphi)
239 
240 
241  RETURN
242  END
243 
244 
245 
246 C*************************************************************************
247  FUNCTION bwsig(XM,XG,XQ)
248 C **********************************************************
249 C Complex Function: S-wave Breit-Wigner
250 C Called by f3pi_rcht.f
251 C **********************************************************
252  IMPLICIT NONE
253  REAL xq
254  DOUBLE PRECISION xm,xg,xm2,xxq,gamma
255  include '../funct_declar.inc'
256  include '../parameter.inc'
257 
258  xxq = xq
259  xm2 = xm**2
260  gamma = xg*sigp(xxq)/sigp(xm2)
261  bwsig = xm*xm/cmplx(xm*xm-xq, -xm*gamma)
262 
263 
264  RETURN
265  END
266 
267 C*************************************************************************
268  FUNCTION decoul(mm1,mm3,ss2)
269 C **********************************************************
270 C Real Function: Coulomb interaction effects of two particles
271 C with mm1 and mm2
272 C Called by f3pi_rcht.f
273 C **********************************************************
274  IMPLICIT NONE
275  REAL ss2
276  DOUBLE PRECISION mm1,mm3,betam1m3
277 
278  include '../funct_declar.inc'
279  include '../parameter.inc'
280 C*******************************************
281 C COMMON block fixed in SUBROUTINE INIPHY
282 c*******************************************
283  COMMON / qedprm /alfinv,alfpi,xk0
284  REAL*8 alfinv,alfpi,xk0
285  betam1m3 = 1.d0 - (mm1 +mm3)**2/ss2
286  betam1m3 = dsqrt(betam1m3)
287 
288  if(ss2.gt.(mm1+mm3)**2)
289  & decoul = pi/2.d0/betam1m3/alfinv
290 
291  RETURN
292  END
293 
294 C*************************************************************************
295  FUNCTION coul3part(mm1,mm2,mm3,ss1,ss2,ss3)
296 C **********************************************************
297 C Real Function: Coulomb interaction effects of two particles
298 C with mm1 and mm2
299 C Called by f3pi_rcht.f
300 C **********************************************************
301  IMPLICIT NONE
302  REAL ss3,ss2,ss1
303  DOUBLE PRECISION mm1,mm2,mm3
304  include '../funct_declar.inc'
305  include '../parameter.inc'
306 
307  coul3part = decoul(mm1,mm3,ss2) + decoul(mm2,mm3,ss1)
308  & - decoul(mm1,mm2,ss3)
309 
310  coul3part = exp(coul3part)
311 
312  RETURN
313  END
314 C*************************************************************************
315  FUNCTION fattcoul(mm1,mm3,ss2)
316 C **********************************************************
317 C Real Function: Coulomb attraction of two particles
318 C with mm1 and mm3
319 C Called by f3pi_rcht.f
320 C **********************************************************
321  IMPLICIT NONE
322  REAL ss2
323  DOUBLE PRECISION mm1,mm3,betam1m3
324 
325  include '../funct_declar.inc'
326  include '../parameter.inc'
327 C*******************************************
328 C COMMON block fixed in SUBROUTINE INIPHY
329 c*******************************************
330  COMMON / qedprm /alfinv,alfpi,xk0
331  REAL*8 alfinv,alfpi,xk0
332  if(ss2.gt.(mm1+mm3)**2) then
333  betam1m3 = 2.*dsqrt(1.d0 - (mm1 +mm3)**2/ss2)
334  & /(1.+ (1.d0 - (mm1 +mm3)**2/ss2))
335 
336 
337  fattcoul = 2.*pi/betam1m3/alfinv
338  & /(1.-exp(-2.*pi/betam1m3/alfinv))
339  else
340  fattcoul = 1
341  endif
342 
343  RETURN
344  END
345 
346 C*************************************************************************
347 C*************************************************************************
348  FUNCTION frepcoul(mm1,mm3,ss2)
349 C **********************************************************
350 C Real Function: Coulomb repuslcion of two particles
351 C with mm1 and mm3
352 C Called by f3pi_rcht.f
353 C **********************************************************
354  IMPLICIT NONE
355  REAL ss2
356  DOUBLE PRECISION mm1,mm3,betam1m3
357 
358  include '../funct_declar.inc'
359  include '../parameter.inc'
360 C*******************************************
361 C COMMON block fixed in SUBROUTINE INIPHY
362 c*******************************************
363  COMMON / qedprm /alfinv,alfpi,xk0
364  REAL*8 alfinv,alfpi,xk0
365  if(ss2.gt.(mm1+mm3)**2) then
366  betam1m3 = 2.*dsqrt(1.d0 - (mm1 +mm3)**2/ss2)
367  & /(1.+ (1.d0 - (mm1 +mm3)**2/ss2))
368 
369  frepcoul = 2.*pi/betam1m3/alfinv
370  & /(-1.+exp(2.*pi/betam1m3/alfinv))
371  else
372  frepcoul = 1
373  endif
374  RETURN
375  END
376 
377 C*************************************************************************