Added comparison with [hep-ph/0205172] and [hep-ph/0505069]

michel-decay/validation/validation.py

+# vim: foldmethod=marker
+from pymule import *
+from pymule.maths import Li2, Li3, zeta2, zeta3
+import pymule.compress
+
+gamma0 = Mmu**5/192/pi**3
+
+## Validation{{{
+## Load data{{{
+setup(folder='out.tar.bz2')
+
+openLeptons = mergefks(sigma('m2ennee0', merge={'xh':2,'xs':2}))
+openLeptons = addplots(openLeptons['xh'], openLeptons['xs'], sa=0.5,sb=0.5)
+
+nnlo = mergefks(
+    sigma('m2ennRF', merge={'xe': 5}),
+    sigma('m2ennRR', merge={'xe': 2}),
+    sigma('m2ennFF', merge={'xe': 2}, folder='ff.tar.bz2', flavour='mu-e')
+)
+xev = np.array([1,1/gamma0,1/gamma0])*addplots(nnlo['xe'], openLeptons)
+##########################################}}}
+## Load logarithms from [hep-ph/0205172] (Arbuzov et al.){{{
+# see also meg2d/arbuzov.py
+
+def fLL(xe):
+    gam = (
+        -32*xe**3/9 
+        + 4*xe**2*(3 - 2*xe)*(
+            -Li2(1 - xe) - zeta2 + np.log(xe)**2/2 
+            - 2*np.log(xe)*np.log(1 - xe) + np.log(1 - xe)**2
+        )
+        + 8*xe**2/3 + 17*xe/6 
+        + (-32*xe**3/3 + 8*xe**2 - 2*xe - 5./6)*np.log(xe) 
+        + (32*xe**3/3 - 16*xe**2 + 8*xe + 10./3)*np.log(1 - xe) + 11./36
+    )
+    NS = (
+        8*xe**3/3 + 2*xe**2*(3 - 2*xe)*np.log((1 - xe)/xe) 
+        - 4*xe**2 + 2*xe + 5./6
+    )
+    S = (
+        -8*xe**3/9 - 14*xe**2/3 + 3*xe 
+        + (4*xe**2 + 4*xe + 5./3)*np.log(xe) + 17./9 + 2/(3*xe)
+    )
+    return gam/2 + NS/3 + S/2
+
+
+def fNLL(xe):
+    def Sot(z):
+        return 0.5*np.log(1-z)**2*np.log(z) + np.log(1-z)*Li2(1-z) - Li3(1-z)+zeta3
+    
+    gam = (
+        -559*xe**3/54 
+        + 2*xe**2*(3 - 2*xe)*(
+            -2*Sot(xe) + 2*Li2(xe)*np.log(xe) + 2*Li2(xe)*np.log(1 - xe)
+            - 2*Li3(xe) - 2*zeta2*np.log(xe) - 2*zeta2*np.log(1 - xe)
+            + 7*zeta3 + 2*np.log(xe)**3 - 5*np.log(xe)**2*np.log(1 - xe) 
+            + 5*np.log(xe)*np.log(1 - xe)**2
+        )
+        + 83*xe**2/3 - 287*xe/12 + (8*xe**3 - 4*xe**2 - 12*xe)*np.log(1 - xe)**2 
+        + (-98*xe**3/3 + 35*xe**2 - 2*xe - 10./3)*zeta2 
+        + (-70*xe**3/3 + 22*xe**2 - 5*xe - 25./12)*np.log(xe)**2 
+        + (-12*xe**3 + 64*xe**2/3 - 53*xe/3 - 17./3)*np.log(1 - xe) 
+        + (44*xe**3/9 + 4*xe**2/3 + 37*xe/6 - 3./4)*np.log(xe) 
+        + (92*xe**3/3 - 54*xe**2 + 32*xe + 25./3)*np.log(xe)*np.log(1 - xe) 
+        + (92*xe**3/3 - 40*xe**2 + 14*xe + 10./3)*Li2(xe) + 211./216
+    )
+    NS = (
+        -64*xe**3/9
+        + 2*xe**2*(3 - 2*xe)*(
+            -2*Li2(1 - xe) - 2*zeta2/3 - np.log(xe)**2 
+            - 2*np.log(xe)*np.log(1 - xe)/3 + 2*np.log(1 - xe)**2/3
+        ) 
+        + 100*xe**2/9 - 19*xe/3 
+        + (-76*xe**3/9 + 8*xe**2 + 4*xe/3 + 5./9)*np.log(xe) 
+        + (12*xe**3 - 46*xe**2/3 - 4*xe/3 + 10./9)*np.log(1 - xe) 
+        - 11./6
+    )
+    S = (
+        10*xe**3/9 + 77*xe**2/18 + 43*xe/18 +
+        (Li2(1 - xe) + np.log(xe)*np.log(1 - xe))*(4*xe**2 + 4*xe + 5./3)
+        + (4*xe**2 + 6*xe + 5./2)*np.log(xe)**2 + (-19*xe**2/3 - 5*xe/6 + 8./9 + 4/(3*xe))*np.log(xe)
+        + (-8*xe**3/9 - 14*xe**2/3 + 3*xe + 17./9 + 2/(3*xe))*np.log(1 - xe)
+        - 67./9 - 1/(3*xe)
+    )
+    int = (
+        104*xe**3/9 
+        + 2*xe**2*(3 - 2*xe)*(-4*Sot(1 - xe) - 2*Li2(1 - xe)*np.log(xe) + 2*Li3(1 - xe))
+        - 55*xe**2/3 + 41*xe/3 + (26*xe**3/3 - 9*xe**2)*np.log(xe)**2 
+        + (-28*xe**2/3 - 5*xe/3 - 5./3)*np.log(xe) 
+        + (52*xe**3/3 - 26*xe**2 + 4*xe + 5./3)*Li2(1 - xe) - 62./9
+    )
+    return gam+NS+S+int
+
+
+def f2loop(xe):
+    z = Mel / Mmu
+    f0 = -np.sqrt(xe**2 - 4*z**2)*(2*xe**2 + 4*z**2 - 3*xe*(1 + z**2))
+    f00 = -xe**2*(-3 + 2*xe)
+    L = np.log(Mmu**2/Mel**2)
+    return f0 / f00 * (L**2 * fLL(xe) + L * fNLL(xe) ) / 2 / pi / pi
+
+###########################################################}}}
+## Load data from [hep-ph/0505069] (Anastasiou et al.){{{
+
+const = np.array(pymule.compress.uncompress(
+    "1:eJx11WtQU0cUAOAQ8AFqLRGFgAqCxAfYagBFCnNT5SEWTGtIHB21QuFCGST0BsSqgHU"
+    "oBWsrAuWRUARtaUu0gkAtFguCRYQEgwgBEiCIQsBoRBRx1LQn/ursTs6Pndn5ZnfP2d3Z"
+    "XREh3ElyaTSayAKakM9FiST9/z3KoCrSmbO7CKrGpWwm22qUQLwjz578+zz4BKtiuV4wg"
+    "rrTa3mc0vYuQWWwrkZrL2hQpwmZDvNzwX1D8sKbZ6lRP7I94lumfTdBzfdTpDn39qD+jp"
+    "gT5VADbueyYg7vZifqiqvM5o6oHoJa/XjN0SXyW6i7p7Hv/zQIPndVv1ow1Ix6Tmmlo/d"
+    "eJUFxqtZbSizrUc+Qt9DHNOCVL+4nhJdVYvXNTvquKq6XoBo/qpmqrcpHPVOSyd78xujP"
+    "VNvN+bkNiJ8rMkYfQZXyAm3cbX5H/V27FP24Qz/4e7p5pw1/oJ6afyDApwTcybWzKJHZg"
+    "Pr0F7/VBrmqYP+WNnzTKGpGPSb0qSv/gtFrvty4x/Mm6rGeXX6pPmqC+pjHOnttuA31Q4"
+    "eZCZFycM5CiT5oSwfq4dvavAWhAwQVR1fMyU5XoN7ecvH78tvgd463Pzsw2In61nPThmr"
+    "BIEHFzK6TvFrZhbqvO/XaUQOuyD/Td2PHXdT9FvFsguOHCGrj5qCYjG3dqN/4laxTGMBL"
+    "H/UFcRg9qC9vUtrot2gIil4TfPipFPPe91/LkjLAS76SRsSzlaj3PKtrd2sHv7Nz9lruJ"
+    "czv1CW3DS8cJihvN570pV0v6h2pXrcKwsBlVjmjB4WYy7ZO3QwtAhdPNqmklzFvnVvVYj"
+    "4ITlPrLu7QYZ5+qDv15Ip7BEW2LDk+6dqH+vEPclM4JLiskrMz71PMj9E+OTZdAa4o+cy"
+    "s+DTm3vFJnimPjX4qo/J8A+aeHsYYgfxOSJhJasw3zDzZMJ4MvjG58kTgS8x1UYsnhH+B"
+    "k3H/PLRd0o/6+NourSvtPkEVh/eHjbMxH9Nnj6n8wWP3aMWOgZiX7as9G5EJLhPor13lY"
+    "l7iLCqxk4OLuW/MJyMwLx5b/2OH9QOTLuAbA5zNm9bwd2MeZr8obJ74gcn1F2Rump8zAG"
+    "6xfzg53Rdzq9DnVtucR036HEa1pYEcNVl/U5ryenrFqMn8G/1/aPTRgyv8GQeb0jCvt+Q"
+    "1PPEYI6hCP6fqmRzMk98GuILttbTtEuaJvl6H1tWDG1Zt/1oiwzzBbCpxhKaF82WS1nOn"
+    "MHcT2q6L9QeXvXP0SN9CFeqrPbvdnLLA6RbZIxXumI+UvRwpl4EXPpcKbXmYa8gr9/Zaj"
+    "0P+E429WiHmA25JwwwBOKlWR/LzMS94+/6DF3dOylnVmOft31XIHQSPbPt30b4WzHNcFh"
+    "fMcpkw6SFnNnNPkeAG+cyT6/WYBwtehG6RTphcP9ChNmRGb3p+s6w+izSPh7B/9Q6lt3S"
+    "YG0ILzL0Og8eUB1BxLDXqrxh8+sN68MKslJMBIsxfKBlmZWY6uH9RzdqwLswvB2yqWxUA"
+    "XvihZVCkxwDql6yeXxnIMo5fxv1ZlIt57Nv/C5w2LWEtf4N5tJ8y3p7xCOrv0P7SSg6i7"
+    "ix9xSrjg7eWe60RKTB3TKhz3SUG90w7UbHMZwj1pV7JKxcMgYt33XZrPY95zzmJOtf5MU"
+    "F55FtrM+ZpUO+K3q0KjgaPpP+5ei9H0/Affiqv0w=="
+))
+
+const[:,1] = (
+    const[:,1] * (
+        -np.sqrt(const[:,0]**2 - 4*(Mel/Mmu)**2)*(
+            2*const[:,0]**2 + 4*(Mel/Mmu)**2 - 3*const[:,0]*(1 + (Mel/Mmu)**2)
+        )
+    ) / alpha**2 *2e-4
+)
+###########################################################}}}
+## Make plot{{{
+xe = np.linspace(0,1,200)[2:-1]
+
+ax = gca()
+errorband(xev, ax=ax, col='C3')
+ax.plot(xe,f2loop(xe), 'C2', linestyle='dashed')
+ax.plot(const[:,0],const[:,1] + f2loop(const[:,0]), 'C1')
+ax.set_ylim(-15,25)
+ax.set_xlabel('$x_e$')
+ax.set_ylabel("$\\D\\sigma^{(2)} / \\D x_e / \\sigma_0$")
+mulify(ax.figure)
+
+ax.legend([
+    '$\\textsc{McMule}$',
+    '\\cite{Arbuzov:2002pp,Arbuzov:2002cn}',
+    '\\cite{Anastasiou2005The-electron}'
+], loc='upper center')
+
+##########################################}}}
+#####################################################################}}}