"""Mellin-Barnes integrations.
Todo:
Various MB routines should be integrated into one.
"""
from math import log, pi
import numpy as np
[docs]
class MellinBarnes(object):
"""Base class for models built by Mellin-Barnes integration."""
[docs]
def __init__(self, **kwargs) -> None:
self.tgj = np.tan(pi*self.jpoints/2)
# print('MellinBarnes init done.')
[docs]
def _mellin_barnes_integral(self, xi, wce, gpd):
"""Return convolution of evolved Wilson coefs and GPDs."""
eph = np.exp(self.phi*1j)
cfacj = eph * xi**(-self.jpoints-1)
cch = np.einsum('j,sa,sja,ja->j', cfacj,
self.pw_strengths(), wce, gpd)
imh = np.dot(self.wg, cch.imag)
np.multiply(cch, self.tgj, out=cch)
reh = np.dot(self.wg, cch.imag)
return reh, imh
[docs]
def _mellin_barnes_integral_HE(self, xi, wce, H, E):
"""Return convolution of evolved Wilson coefs and GPDs H and E."""
eph = np.exp(self.phi*1j)
cfacj = eph * xi**(-self.jpoints-1)
cch_H = np.einsum('j,sa,sja,ja->j', cfacj,
self.pw_strengths(), wce, H)
cch_E = np.einsum('j,sa,sja,ja->j', cfacj,
self.pw_strengths_E(), wce, E)
imh = np.dot(self.wg, cch_H.imag)
ime = np.dot(self.wg, cch_E.imag)
np.multiply(cch_H, self.tgj, out=cch_H)
np.multiply(cch_E, self.tgj, out=cch_E)
reh = np.dot(self.wg, cch_H.imag)
ree = np.dot(self.wg, cch_E.imag)
return reh, imh, ree, ime
[docs]
def _dis_mellin_barnes_integral(self, xi, wce, pdf):
"""Return convolution of evolved Wilson coefs and PDFs."""
eph = np.exp(self.phi*1j)
cfacj = eph * np.exp((self.jpoints) * log(1/xi)) # eph/xi**j
cch = np.einsum('j,ja,ja->j', cfacj, wce, pdf)
mb_int = np.dot(self.wg, cch.imag)
return mb_int
[docs]
def _j2x_mellin_barnes_integral(self, x, eta, wce, gpd):
"""Return convolution of j->x coef, evolution operator and GPD."""
# difference wrt above integrations is that here we do NOT sum over flavors
eph = np.exp(self.phi*1j)
cfacj = eph * x**(-self.jpoints-1) / np.pi
if eta < 1e-8 and x>0:
# forward limit, PDF-like, so only zero-th PW is taken
cch = np.einsum('k,ia,kab,kb->ki', cfacj, self.antifrot_pdf, wce[0, :, :, :], gpd)
# in evol basis:
# cch = np.einsum('j,jab,jb->ja', cfacj, wce[0, :, :, :], gpd)
elif abs(eta-x) < 1e-8 and x>0:
# cross-over, border eta=x limit
cch = np.einsum('j,sa,sjab,jb->ja', cfacj,
self.pw_strengths(), wce, gpd)
else:
cch = eph * np.einsum('si,ia,skab,kb->ki', self.pw_strengths(),
self.antifrot_pdf, wce, gpd)
# If you want result in evolution basis:
# cch = eph * np.einsum('sa,skab,kb->ka', self.pw_strengths(),
# wce, gpd)
#raise Exception('eta has to be either 0 or equal to x')
mb_int_flav = np.dot(self.wg, cch.imag)
return mb_int_flav
[docs]
def _j2x_mellin_barnes_integral_E(self, x, eta, wce, gpd):
"""Return convolution of j->x coef, evolution operator and GPD."""
# difference wrt above integrations is that here we do NOT sum over flavors
eph = np.exp(self.phi*1j)
cfacj = eph * x**(-self.jpoints-1) / np.pi
if eta < 1e-8 and x>0:
# forward limit, PDF-like, so only zero-th PW is taken
cch = np.einsum('j,jab,jb->ja', cfacj, wce[0, :, :, :], gpd)
elif abs(eta-x) < 1e-8 and x>0:
# cross-over, border eta=x limit
cch = np.einsum('j,sa,sjab,jb->ja', cfacj,
self.pw_strengths_E(), wce, gpd)
else:
cch = eph * np.einsum('sa,sjab,jb->ja',
self.pw_strengths_E(), wce, gpd)
#raise Exception('eta has to be either 0 or equal to x')
mb_int_flav = np.dot(self.wg, cch.imag)
return mb_int_flav