C++ Interface to Tauola
funct_rpt.f
1
2C **********************************************************
3C library of functions used in calculation of currents
4C references:
5C [1] arXiv:0911.4436 (hep-ph) D. Gomez Dumm et al. (tau -> 3pi nu)
6C [2] arXiv:0911.2640 (hep-ph) D. Gomezz Dumm et al. (tau -> KKpi nu)
7C [3] arXiv:0807.4883 (hep-ph) D. R. Boito et al. (tau -> Kpi nu)
8C [4] P. Roig, talk at (tau -> 2 pi nu)
9C [5] arXiv:0803.2039 (hep-ph) E. Arganda et al., Appendix B (tau -> KK nu)
10C **********************************************************
11 FUNCTION grho_rcht(XS,XMMM)
12 IMPLICIT NONE
13 REAL XS
14 DOUBLE PRECISION XMMM
15C **********************************************************
16C REAL FUNCTION Gamma Rho ; energy-dependent width of rho meson
17c in SU(2) limit mpi=mpi0, mk=mk0;
18C formula (14) of REF [1]
19C **********************************************************
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
48C **********************************************************
49C REAL FUNCTION Gamma Rho1 ; energy-dependent width of rho1 meson;
50c in SU(2) limit mpi=mpi0; only rho' -> 2pi loop;
51C formula (33) of REF [1]
52C **********************************************************
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
74C***********************************************************
75C DOUBLRE PRECISION FUNCTION
76c of two body phase space threshold: equal mass scalars
77C***********************************************************
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
99C***********************************************************
100C REAL FUNCTION LAMDA of three body phase space
101C***********************************************************
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
116C******************** FUNCTIONS FOR 2 SCALAR MODES ***********************
117
118
119
120C****************************************************************
121 FUNCTION r0scal_3pi(QX,SX)
122C****************************************************************
123C Complex Function: R_0 function (Pablo notes)
124C for the scalar contribution for three pion modes
125C Called by f3pi_rcht.f
126C **********************************************************
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'
132c$$$ a00_3piscal = 0.220
133c$$$ b00_3piscal = 0.268/mmpi_av**2
134c$$$ c00_3piscal = -0.0139/mmpi_av**4
135c$$$ d00_3piscal = -0.00139/mmpi_av**6
136c$$$ x00_3piscal = 36.77*mmpi_av**2
137c$$$c MMF0 = 0.98
138c$$$
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
169C****************************************************************
170 FUNCTION r2scal_3pi(QX,SX)
171C****************************************************************
172C Complex Function: R_2 function (Pablo notes)
173C for the scalar contribution for three pion modes
174C Called by f3pi_rcht.f
175C **********************************************************
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'
181c$$$ a02_3piscal = -0.0444
182c$$$ b02_3piscal = -0.0857/mmpi_av**2
183c$$$ c02_3piscal = -0.00221/mmpi_av**4
184c$$$ d02_3piscal = -0.000129/mmpi_av**6
185c$$$ x02_3piscal = -21.62*mmpi_av**2
186c$$$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
220C*************************************************************************
221C Functions for sigma contributions
222C*************************************************************************
223 FUNCTION ffsig(QX,XX)
224C **********************************************************
225C Complex Function:
226C Called by f3pi_rcht.f
227C **********************************************************
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
246C*************************************************************************
247 FUNCTION bwsig(XM,XG,XQ)
248C **********************************************************
249C Complex Function: S-wave Breit-Wigner
250C Called by f3pi_rcht.f
251C **********************************************************
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
267C*************************************************************************
268 FUNCTION decoul(mm1,mm3,ss2)
269C **********************************************************
270C Real Function: Coulomb interaction effects of two particles
271C with mm1 and mm2
272C Called by f3pi_rcht.f
273C **********************************************************
274 IMPLICIT NONE
275 REAL ss2
276 DOUBLE PRECISION mm1,mm3,betam1m3
277
278 include '../funct_declar.inc'
279 include '../parameter.inc'
280C*******************************************
281C COMMON block fixed in SUBROUTINE INIPHY
282c*******************************************
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
294C*************************************************************************
295 FUNCTION coul3part(mm1,mm2,mm3,ss1,ss2,ss3)
296C **********************************************************
297C Real Function: Coulomb interaction effects of two particles
298C with mm1 and mm2
299C Called by f3pi_rcht.f
300C **********************************************************
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
314C*************************************************************************
315 FUNCTION fattcoul(mm1,mm3,ss2)
316C **********************************************************
317C Real Function: Coulomb attraction of two particles
318C with mm1 and mm3
319C Called by f3pi_rcht.f
320C **********************************************************
321 IMPLICIT NONE
322 REAL ss2
323 DOUBLE PRECISION mm1,mm3,betam1m3
324
325 include '../funct_declar.inc'
326 include '../parameter.inc'
327C*******************************************
328C COMMON block fixed in SUBROUTINE INIPHY
329c*******************************************
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
346C*************************************************************************
347C*************************************************************************
348 FUNCTION frepcoul(mm1,mm3,ss2)
349C **********************************************************
350C Real Function: Coulomb repuslcion of two particles
351C with mm1 and mm3
352C Called by f3pi_rcht.f
353C **********************************************************
354 IMPLICIT NONE
355 REAL ss2
356 DOUBLE PRECISION mm1,mm3,betam1m3
357
358 include '../funct_declar.inc'
359 include '../parameter.inc'
360C*******************************************
361C COMMON block fixed in SUBROUTINE INIPHY
362c*******************************************
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
377C*************************************************************************