TEvtGen/Tauola/tauola-fortran/new-currents/RChL-currents/rcht_common/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*************************************************************************