Integrate over bins

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

@@ -4,6 +4,7 @@ from pymule import *
 from pymule.plot import twopanel
 import arbuzov
 from scipy import optimize
+import scipy.integrate as integrate
 import scipy.stats
 
 gamma0 = Mmu**5/192/pi**3
@@ -341,16 +342,43 @@ print "Agreement for g2: ", printnumber(resg2[4]/(eplogcoeff**2/2))
 # 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
+def dointegral(integrand):
+    diff = np.concatenate((
+        0.026*np.ones(981),
+        0.016*np.ones(1000),
+        0.008*np.ones(1000),
+        0.004*np.ones(706)
+    ))
+
+    return np.array([
+        np.insert(
+            np.array(integrate.quad(
+                integrand, i, j, epsrel=1e-5
+            )) / (j-i),
+            0,
+            (i + j)/2.
+        )
+        for i, j in np.column_stack((Ee-diff/2, Ee+diff/2))
+    ])
+
+
+def yfsintegrand(Ee, o):
+    xe = 2*Ee/Mmu
+
+    beta = sqrt(1-Mel**2/Ee**2)
+    fLO = 2/Mmu * xe**2*beta * (3-2*xe+xe/4*(3*xe-4)*(1-beta**2))
+    gLO = 2/Mmu * xe**2*beta * ((1-2*xe)*beta + 3*xe**2/4 * beta * (1-beta**2))
+
+    logx = np.log(1. - xe + (Mel/Mmu)**2)
+    eplog = logx * eplogcoeff
+    epsub = np.exp(eplog) - (1. + eplog + eplog**2/2.)
+    return [fLO * epsub, gLO * epsub][o]
+
+
+fLLYFS = dointegral(lambda e: yfsintegrand(e, 0))
+gLLYFS = dointegral(lambda e: yfsintegrand(e, 1))
+
 
 FYFS = addplots(F, fLLYFS)
 GYFS = addplots(G, gLLYFS)
@@ -365,14 +393,25 @@ GYFS = addplots(G, gLLYFS)
 # \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
+def collintegrand(Ee, o):
+    pref3loop = (alpha*L/(2*pi))**3 / 6.
+
+    xe = 2*Ee/Mmu
+    logx = np.log(1. - xe + (Mel/Mmu)**2)
+    if o == 0:
+        return pref3loop * (arbuzov.f3LL(xe) - 8 * logx**3)
+    else:
+        return pref3loop * (arbuzov.g3LL(xe) + 8 * logx**3)
+
+
+f3LLsub = dointegral(lambda e: collintegrand(e, 0))
+g3LLsub = dointegral(lambda e: collintegrand(e, 1))
+
+
+FFULL = addplots(FYFS, f3LLsub)
+GFULL = addplots(GYFS, g3LLsub)
+
 ###########################################################}}}
 ##########################################################################}}}