Added f and g analysis

michel-decay/f-and-g/arbuzov.py

+import requests
+import zlib
+import sympy as sp
+import sympy.parsing.latex
+import scipy.special
+import numpy as np
+import StringIO
+import tarfile
+import pymule.maths
+
+
+def dilog(z):
+    return pymule.maths.Li2(z)
+
+
+def trilog(z):
+    return pymule.maths.Li3(z)
+
+
+def Sot(z):
+    return (
+        0.5*np.log(1-z)**2*np.log(z) + np.log(1-z)*dilog(1-z)
+        - trilog(1-z) + pymule.maths.zeta3
+    )
+
+
+def Phi(z):
+    return dilog(-(1-z)/z) + np.log((1-z)/z)**2-np.pi**2/6.
+
+
+def phi(z):
+    return dilog(-(1-z)/z) + sp.log((1-z)/z) - pi**2/6
+
+
+def parse(lines, s):
+    return sympy.parsing.latex.parse_latex(
+        ''.join(lines)
+        .replace("f^{\\mathrm{Born}}(x)", "f*")
+        .replace("g^{\\mathrm{Born}}(x)", "g*")
+        .replace("\\frac{q^2}{m_\\mu^2}", "\\qq")
+        .replace("\\nonumber", "")
+        .replace("\\biggr", "")
+        .replace("\\biggl", "")
+        .replace("\\left", "")
+        .replace("\\right", "")
+        .replace("\\\\", "")
+        .replace("&", "")
+        .replace("\\{", "(")
+        .replace('\\}', ")")
+        .replace('\\Sot{x}', "\\Sot(x)")
+        .replace('\\Sot{1-x}', "\\Sot(1-x)")
+        .replace('\\Li_2', "\\dilog")
+        .replace('\\Li_3', "\\trilog")
+        .replace('\\Li_3', "\\trilog")
+        .replace('\\zeta_2', "\\zeta(2)")
+        .replace('\\zeta_3', "\\zeta(3)")
+        .replace(
+            r'\Li{2}{\frac{(1-\beta)(2-x(1+\beta))}{(1+\beta)(2-x(1-\beta))}}',
+            r'\dilog(\frac{(1-\beta)(2-x(1+\beta))}{(1+\beta)(2-x(1-\beta))})'
+        ).replace(
+            r'\Li{2}{\frac{2-x(1+\beta)}{2-x(1-\beta)}}',
+            r'\dilog(\frac{2-x(1+\beta)}{2-x(1-\beta)})'
+        ).replace('x(', 'x*(')
+        .replace('L (', 'L*(')
+        .replace('beta(', 'beta*(')[s:-1]).subs('x', 'xe')
+
+
+x, xe, beta, z = sp.var('x xe beta z')
+
+r = requests.get('https://arxiv.org/e-print/hep-ph/0110047')
+tex = zlib.decompress(r.content, 16+zlib.MAX_WBITS)
+
+flo = 3-2*xe+xe/4*(3*xe-4)*(1-beta**2)
+glo = (1-2*xe)*beta+3*xe**2/4*(1-beta**2)*beta
+
+fnlo = parse(tex.splitlines()[231:262], 9)
+Adef = parse(tex.splitlines()[265:276], 3)
+gnlo = parse(tex.splitlines()[279:326], 9)
+
+subs = {
+    'f': flo, 'g': glo, 'A': Adef,
+    'beta': sp.sqrt(1 - 4*z**2/xe**2),
+    'L': -2*sp.log(z), 'qq': (1-xe)+z**2
+}
+
+fnlo = fnlo.subs(subs).subs(subs)
+gnlo = gnlo.subs(subs).subs(subs)
+
+f1 = sp.lambdify((xe, z), fnlo, modules=['numpy', {'dilog': dilog}])
+g1 = sp.lambdify((xe, z), gnlo, modules=['numpy', {'dilog': dilog}])
+
+
+r = requests.get('https://arxiv.org/e-print/hep-ph/0202102')
+fp = StringIO.StringIO(r.content)
+tar = tarfile.open(fileobj=fp)
+tex = tar.extractfile('mu_lla_2.tex').read()
+
+fLL = parse(tex.splitlines()[322:327], 1)\
+      .subs('z', 'xe/(1+m**2/M**2)').subs('m', 'z*M')
+gLL = parse(tex.splitlines()[330:335], 1)\
+      .subs('z', 'xe/(1+m**2/M**2)').subs('m', 'z*M')
+
+f2LL = sp.lambdify((xe, z), fLL, modules=['numpy', {'Phi': Phi}])
+g2LL = sp.lambdify((xe, z), gLL, modules=['numpy', {'Phi': Phi}])
+
+
+r = requests.get('https://arxiv.org/e-print/hep-ph/0205172')
+fp = StringIO.StringIO(r.content)
+tar = tarfile.open(fileobj=fp)
+tex = tar.extractfile('mu_ff_k4.tex').read()
+
+fLL2 = parse(tex.splitlines()[542:550], 21)
+fNLL = parse(tex.replace("\\ln^2x \\biggl(", "\\ln(x)^2 * (")
+                .replace('\\ln x \n\\biggl(', '\\ln(x)*\n(')
+                .splitlines()[557:577], 0)
+
+fLLns = parse(tex.splitlines()[598:600], 20)
+fNLLns = parse(tex.replace('\\ln x \\biggl(', '\\ln (x)* \\biggl(')
+                  .splitlines()[602:610], 0)
+fLLs = parse(tex.splitlines()[625:628], 0)
+fNLLs = parse(tex.replace("\\ln x\\biggl(", "\\ln(x)\\biggl(")
+                 .replace("\\ln^2x \\biggr(", "\\ln(x)^2* \\biggr(")
+                 .splitlines()[630:641], 0)
+fNLLint = parse(tex.replace("\\ln^2x \\biggl(", "\\ln(x)^2 * (")
+                   .replace('\\ln x \\biggl(', '\\ln(x)*(')
+                   .splitlines()[655:662], 0)
+
+# print "{" + ",\n".join(
+#     str(i) for i in [fLL2, fNLL, fLLns, fNLLns, fLLs, fNLLs, fNLLint]
+# ) + "}"

michel-decay/f-and-g/f-and-g.py

+# vim: foldmethod=marker
+## Init{{{
+from pymule import *
+from pymule.plot import twopanel
+from arbuzov import f1, g1, f2LL, g2LL
+
+gamma0 = Mmu**5/192/pi**3
+
+# We *define* $f$ and $g$ through
+# \begin{align}
+# \frac{\mathrm{d}^2\Gamma}{\mathrm{d}E_e\,\mathrm{d}(\cos\theta)}
+#   = \Gamma_0\Big( f(E_e) + P\,\cos\theta\ g(E_e)\Big)\,.
+# \end{align}
+# At LO these are for $m=0$
+# \begin{align}
+#  f^{(0)}(E_e) &= (\tfrac{2E_e}{M})^2\big(3-2(\tfrac{2E_e}{M})\big)\,,\\
+#  g^{(0)}(E_e) &= (\tfrac{2E_e}{M})^2\big(1-2(\tfrac{2E_e}{M})\big)\,.
+# \end{align}
+# We now calculate $f$ and $g$ by calculating with $P=0.85$
+# \begin{align}
+# \frac{\mathrm{d}\Gamma_-}{\mathrm{d} E_e}
+#        =\int_{-1}^0\mathrm{d}(\cos\theta)\ \frac{\mathrm{d}^2\Gamma}
+#                                    {\mathrm{d} E_e\,\mathrm{d}(\cos\theta)}
+# \qquad\text{and}\qquad
+# \frac{\mathrm{d}\Gamma_+}{\mathrm{d} E_e}
+#        =\int_{ 0}^1\mathrm{d}(\cos\theta)\ \frac{\mathrm{d}^2\Gamma}
+#                                    {\mathrm{d} E_e\,\mathrm{d}(\cos\theta)}
+# \,.
+# \end{align}
+# Allowing us to obtain $f$ and $g$.
+#
+#
+# To get a good enough sampling, we split the numerical $E_e$ integration in
+# four integrals
+# \begin{align}
+#  \int_m^{E_\text{max}}\mathrm{d}E_e =
+#  \int_{ 0}^{26}\mathrm{d}E_e
+# +\int_{26}^{42}\mathrm{d}E_e
+# +\int_{42}^{50}\mathrm{d}E_e
+# +\int_{50}^{54}\mathrm{d}E_e
+# \end{align}
+# This is implemented in `user.f95` as follows
+# ```fortran
+#   ! 1st digit is cos sign (1 = -, 2 = +)
+#   ! 2nd digit is the energy range
+#
+#   whatplot = mod(p_encrypted,10)
+#
+#   if(p_encrypted/10 == 1) then
+#     whatplot = - whatplot
+#   endif
+#  ...
+#   if(abs(whatplot) == 1) then
+#      if(q2(4) > 26.) set_zero = 0
+#    elseif(abs(whatplot) == 2) then
+#      if(q2(4) < 26. .or. q2(4) > 42.) set_zero = 0
+#    elseif(abs(whatplot) == 3) then
+#      if(q2(4) < 42. .or. q2(4) > 50.) set_zero = 0
+#    elseif(abs(whatplot) == 4) then
+#      if(q2(4) < 50.) set_zero = 0
+#    endif
+#    if (whatplot > 0) then
+#      ! We are in the cos > 0 region
+#      if (cos_th(q2,ez) < 0) set_zero = 0
+#    else
+#      ! We are in the cos < 0 region
+#      if (cos_th(q2,ez) > 0) set_zero = 0
+#    endif
+# ```
+##########################################################################}}}
+## Load data{{{
+setup(folder='meg-2d/out.tar.bz2')
+
+
+### Helper functions{{{
+# These functions load the full dataset for a given sign
+def lfunclo(s):
+    return [
+        scaleset(mergefks(
+            sigma('m2enn0', obs=s+obs)
+        ), alpha**0/gamma0)
+
+        for obs in ['1', '2', '3', '4']
+    ]
+
+
+def lfuncnlo(s):
+    return [
+        scaleset(mergefks(
+            sigma('m2ennR', obs=s+obs),
+            sigma('m2ennF', obs=s+obs)
+        ), alpha**1/gamma0)
+
+        for obs in ['1', '2', '3', '4']
+    ]
+
+
+def lfuncnnlo(s):
+    return [
+        scaleset(mergefks(
+            sigma(
+                'm2ennRF', obs=s+obs,
+                sanitycheck=lambda v:abs(v['value'])[0] < 1e10
+            ),
+            sigma('m2ennRR', obs=s+obs),
+            sigma('m2ennFF', obs=s+obs)
+        ), alpha**2/gamma0)
+
+        for obs in ['1', '2', '3', '4']
+    ]
+
+
+# This functions merges the different runs for one sign into one large run that
+# has full coverage of the spectrum
+def mergesign(lfunc, s, name):
+    sets = lfunc(s)
+
+    dic = {
+        'value': plusnumbers([i['value'] for i in sets]),
+        'chi2a': sum(
+            i['chi2a'][0][0] if type(i['chi2a'][0]) == list else i['chi2a'][0]
+            for i in sets
+        )/len(sets),
+        'time': sum(i['time'] for i in sets)
+    }
+    for i in range(4):
+        dic['Ee%d' % (i+1)] = sets[i]['Ee%d' % (i+1)][1:-1]
+
+    return dic
+
+
+###########################################################}}}
+### Assemble $f$ and $g${{{
+# We can now build $f$ and $g$ as
+# \begin{align}
+# f &= \frac{\Gamma_- + \Gamma_+}{2}\,,\\
+# g &= \frac{-\Gamma_- + \Gamma_+}{0.85}
+# \end{align}
+def getfandg(lfunc, name):
+    r1 = mergesign(lfunc, '1', name)
+    r2 = mergesign(lfunc, '2', name)
+
+    espec1 = np.concatenate((r1['Ee1'], r1['Ee2'], r1['Ee3'], r1['Ee4']))
+    espec2 = np.concatenate((r2['Ee1'], r2['Ee2'], r2['Ee3'], r2['Ee4']))
+
+    return (
+        addplots(espec1, espec2, sa=0.5, sb=0.5),
+        addplots(espec1, espec2, sa=-1/0.85, sb=1/0.85)
+    )
+
+
+flo, glo = getfandg(lfunclo, 'lo')
+fnlo, gnlo = getfandg(lfuncnlo, 'nlo')
+fnnlo, gnnlo = getfandg(lfuncnnlo, 'nnlo')
+
+
+F = addplots(flo, addplots(fnlo, fnnlo))
+G = addplots(glo, addplots(gnlo, gnnlo))
+###########################################################}}}
+##########################################################################}}}
+## Compare at LO and NLO{{{
+# Here we compare the results with Arbuzov [??]
+Ee = flo[:, 0]
+xe = 2*Ee/Mmu
+beta = sqrt(1-Mel**2/Ee**2)
+
+### Compare at LO{{{
+fLOref = 2/Mmu * xe**2*beta * (3-2*xe+xe/4*(3*xe-4)*(1-beta**2))
+gLOref = 2/Mmu * xe**2*beta * ((1-2*xe)*beta + 3*xe**2/4 * beta * (1-beta**2))
+fLOref[Ee < Mel] = 0.
+gLOref[Ee < Mel] = 0.
+fLOref[Ee > Mmu/2 + Mel**2/Mmu**2] = 0.
+gLOref[Ee > Mmu/2 + Mel**2/Mmu**2] = 0.
+
+
+twopanel(
+    "$E_e$", labupleft="$f_0$",
+    upleft=[flo, np.column_stack((Ee, fLOref, np.zeros((4000))))],
+    downleft=[np.column_stack((Ee, 1-flo[:, 1]/fLOref, flo[:, 2]/fLOref))]
+)
+ylim(-1e-3, 1e-3)
+
+
+twopanel(
+    "$E_e$", labupleft="$g_0$",
+    upleft=[glo, np.column_stack((Ee, gLOref, np.zeros((4000))))],
+    downleft=[np.column_stack((Ee, 1-glo[:, 1]/gLOref, glo[:, 2]/gLOref))]
+)
+ylim(-5e-2, 5e-2)
+###########################################################}}}
+### Compare at NLO{{{
+fNLOref = alpha/(2*pi) * 2/Mmu * xe**2*beta * f1(xe, Mel/Mmu)
+gNLOref = alpha/(2*pi) * 2/Mmu * xe**2*beta * g1(xe, Mel/Mmu)
+fNLOref[Ee < Mel] = 0.
+gNLOref[Ee < Mel] = 0.
+fNLOref[Ee > Mmu/2 + Mel**2/Mmu**2] = 0.
+gNLOref[Ee > Mmu/2 + Mel**2/Mmu**2] = 0.
+
+twopanel(
+    "$E_e$", labupleft="$f_1$",
+    upleft=[fnlo, np.column_stack((Ee, fNLOref, np.zeros((4000))))],
+    downleft=[np.column_stack((Ee, 1-fnlo[:, 1]/fNLOref, fnlo[:, 2]/fNLOref))]
+)
+ylim(-2e-3, 2e-3)
+twopanel(
+    "$E_e$", labupleft="$g_1$",
+    upleft=[gnlo, np.column_stack((Ee, gNLOref, np.zeros((4000))))],
+    downleft=[np.column_stack((Ee, 1-gnlo[:, 1]/gNLOref, gnlo[:, 2]/gNLOref))]
+)
+ylim(-2e-3, 2e-3)
+###########################################################}}}
+### Compare at NNLO{{{
+fNNLOref = (alpha/(2*pi))**2 * 2/Mmu * (np.log(Mmu**2/Mel**2)-1)**2/2 * f2LL(xe, Mel/Mmu)  # nopep8
+gNNLOref = (alpha/(2*pi))**2 * 2/Mmu * (np.log(Mmu**2/Mel**2)-1)**2/2 * g2LL(xe, Mel/Mmu)  # nopep8
+fNNLOref[Ee < Mel] = 0.
+gNNLOref[Ee < Mel] = 0.
+fNNLOref[Ee > Mmu/2 + Mel**2/Mmu**2] = 0.
+gNNLOref[Ee > Mmu/2 + Mel**2/Mmu**2] = 0.
+###########################################################}}}
+
+##########################################################################}}}