Resum f and g

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

@@ -3,6 +3,8 @@
 from pymule import *
 from pymule.plot import twopanel
 import arbuzov
+from scipy import optimize
+import scipy.stats
 
 gamma0 = Mmu**5/192/pi**3
 
@@ -232,5 +234,145 @@ ylim(-2e-5, 5e-5)
 xlabel("$E_e$")
 ylabel("$g_2$")
 ###########################################################}}}
+##########################################################################}}}
+## Beyond NNLO{{{
+### Prologue{{{
+# We now want to go beyond NNLO, especially towards the endpoint. For this we
+# have two options:
+#  * resum the endpoint using YFS and
+#  * include the knwon N$^3$LO-LL terms.
+# We start with the first point. Using Massimo's method, we expand $f_1/f_0$
+# around the endpoint. We do this in the notation of [1811.06461] with $x =
+# (p-q)^2/M^2$. We find
+# \begin{align}
+#  {\tt eplog} = \frac{f_1}{f_0} \sim \frac{g_1}{g_0} \sim
+#  \frac{\alpha \log(x)}{\pi}\Bigg(
+#    -2+\frac{1+z^2}{-1+z^2}\log\frac{m^2}{M^2}
+#  \Bigg) + \text{regular as $x\to0$}\,,\\
+# \end{align}
+
+z = Mel / Mmu
+eplogcoeff = alpha/pi*(
+    -2. + (1.+z**2)/(-1.+z**2) * np.log(z**2)
+)
+# The resummed endpoint is now just $e^{\tt eplog} = 1+{\tt eplog} + \frac{{\tt
+# eplog}^2}{2!}+\mathcal{O}((\alpha\log x)^3)$. Let us compare this at fixed
+# order.
+#################################################}}}
+### Compare endpoint at fixed order{{{
+# By re-expanding $e^{\tt eplog}$ we obtain only the $\log x$. To compare with
+# McMule we need to add higher-order terms. Hence, we fit the McMule end-point
+# as
+# \begin{align}
+#  f_1 &=            b\log x + c\ x^2+d\ x+e\,,\\
+#  f_2 &= a\log^2x + b\log x + c\ x^2+d\ x+e\,.
+# \end{align}
+#### Fitting{{{
+
+
+def fitfunc(p, Ee, o):
+    x = 1. - 2*Ee/Mmu + (Mel/Mmu)**2
+    ans = p[0] + p[1] * x + p[2] * x**2 + np.log(x) * p[3]
+    if o == 2:
+        ans += np.log(x)**2 * p[4]
+    return ans
 
+
+def fitep(data, o, lab):
+    xdata = data[2730:, 0]
+    ydata = data[2730:, 1]
+    edata = data[2730:, 2]
+
+    def errfunc(p, x, y, err):
+        return (y - fitfunc(p, x, o)) / err
+
+    pinit = [1.0] * (o+3)
+    out = optimize.leastsq(errfunc, pinit,
+                           args=(xdata, ydata, edata), full_output=1)
+
+    coeff = out[0]
+    covar = out[1]
+
+    def band(Ee):
+        cf = 0.95
+        x = 1. - 2*Ee/Mmu + (Mel/Mmu)**2
+        tval = scipy.stats.t.ppf((cf+1)/2, len(data) - len(coeff))
+        if o == 1:
+            xvec = np.array([1, x, x**2, np.log(x)])
+        elif o == 2:
+            xvec = np.array([1, x, x**2, np.log(x), np.log(x)**2])
+
+        return tval*np.sqrt(np.matmul(np.matmul(xvec, covar), xvec))
+
+    fitv = np.zeros(data[2730:].shape)
+    fitv[:, 0] = xdata
+    fitv[:, 1] = fitfunc(coeff, xdata, o)
+    fitv[:, 2] = [band(e) for e in fitv[:, 0]]
+    # One probably could calculate fitv[:,2] as well..
+
+    twopanel(
+        "$E_e$", labupleft=lab, labdownleft="res.",
+        upleft=[data[2735:], fitv],
+        downleft=[divideplots(data, fitv, offset=-1)]
+    )
+    return np.column_stack((coeff, np.sqrt(np.diag(covar))))
+
+
+# Let's perform the fit
+resf1 = fitep(divideplots(fnlo, flo), 1, "$f_1$")
+ylim(-3e-4, 5e-4)
+resg1 = fitep(divideplots(gnlo, glo), 1, "$g_1$")
+ylim(-2e-3, 2e-3)
+
+
+resf2 = fitep(divideplots(fnnlo, flo), 2, "$f_2$")
+ylim(-2e-3, 2e-3)
+resg2 = fitep(divideplots(gnnlo, glo), 2, "$g_2$")
+ylim(-2e-2, 2e-2)
+#################################################}}}
+#### Compare with resummation{{{
+print "Agreement for f1: ", printnumber(resf1[3]/eplogcoeff)
+print "Agreement for g1: ", printnumber(resg1[3]/eplogcoeff)
+print "Agreement for f2: ", printnumber(resf2[4]/(eplogcoeff**2/2))
+print "Agreement for g2: ", printnumber(resg2[4]/(eplogcoeff**2/2))
+#################################################}}}
+###########################################################}}}
+### Calculate YFS spectrum{{{
+# We can now calculate the spectrum, including soft logarithms by taking the
+# McMule data and adding the soft endpoint starting at $\alpha^3$.
+
+logx = np.log(1. - 2*Ee/Mmu + (Mel/Mmu)**2)
+eplog = logx * eplogcoeff
+epsub = np.exp(eplog) - (1. + eplog + eplog**2/2.)
+
+fLLYFS = flo.copy()
+gLLYFS = glo.copy()
+fLLYFS[:, 1] = fLLYFS[:, 1] * epsub
+fLLYFS[:, 2] = fLLYFS[:, 2] * epsub
+gLLYFS[:, 1] = gLLYFS[:, 1] * epsub
+gLLYFS[:, 2] = gLLYFS[:, 2] * epsub
+
+FYFS = addplots(F, fLLYFS)
+GYFS = addplots(G, gLLYFS)
+
+###########################################################}}}
+### Add Arbuzov's $(\alpha \log z)^3${{{
+# Arbuzov's $(\alpha\log z)^3$ is
+# \begin{align}
+#   f_3&= \Big(\frac{\alpha}{2\pi}L\Big)^3\frac1{6}\ f_3^{(0,\gamma)} + ...
+#   \,,\\
+#  f_3^{(0,\gamma)} &= {\tt f3LL} \sim 8\times\log^3x + ...\,.
+# \end{align}
+# This is the term we have to avoid to not double-count the soft-collinear
+# logarithm at $\alpha^3$.
+pref3loop = (alpha*L/(2*pi))**3 / 6.
+f3LLsub = pref3loop * (arbuzov.f3LL(xe) - 8 * logx**3)
+g3LLsub = pref3loop * (arbuzov.g3LL(xe) + 8 * logx**3)
+
+FFULL = FYFS.copy()
+GFULL = GYFS.copy()
+
+FFULL[:, 1] += f3LLsub
+GFULL[:, 1] += g3LLsub
+###########################################################}}}
 ##########################################################################}}}