Source code for gepard.wilson

"""Wilson coefficients and evolved Wilson coefficients."""

import numpy as np
from scipy.special import loggamma  # type: ignore

from . import adim, c1dvcs, c1dvmp, constants, evolution, qcd, theory, gegenbauer


[docs] def _fshu(j: np.ndarray) -> np.ndarray: """Shuvaev factor.""" # = 2^(J+1) Gamma(5/2+J) / Gamma(3/2) / Gamma(3+J) = # = 2^N Gamma(3/2+N) / Gamma(3/2) / Gamma(2+N) return (2**(j+1) * np.exp(loggamma(2.5 + j) - loggamma(3 + j) - loggamma(3/2)))
[docs] def calc_wc(m: theory.Theory, process_class: str): """Calculate Wilson coeffs. Args: m: instance of the Theory class process_class: 'DIS', 'DVCS' or 'DVMP' Returns: wc[s, k, p, a]: s in range(npws), k in range(npts), p in [LO, NLO], a in [Q, G, NSP, NSM] """ # Initialize array of correct shape with complex zeros nebas = len(m.evolution_basis) wc = np.zeros((m.npws, m.npts, m.p + 1, nebas), dtype=complex) # LO hard scattering coefs are equal to 1 for all flavors ... wc[:, :, 0, :] = np.ones_like(wc[:, :, 0, :]) # ... apart from gluons which are now set to zero: if 'G' in m.evolution_basis: wc[:, :, 0, m.evolution_basis.index('G')] = np.zeros_like(wc[:, :, 0, 0]) # Calculate NLO coefs if needed: if m.p == 1: # FIXME: C1 should calculate what's needed and we should not cut here wc[:, :, 1, :] = c1dvcs.C1(m, process_class)[..., :nebas] if process_class == 'DIS': prefac = np.ones_like(m.jpoints_pws) elif process_class == 'DVCS': prefac = _fshu(m.jpoints_pws) else: raise ValueError('process_class {} is not "DIS" or "DVCS"!'.format(process_class)) wc = np.einsum('sk,skpa->skpa', prefac, wc) return wc
[docs] def calc_wc_dvmp(m: theory.Theory, g: int): """Calculate Wilson coeffs for DVMP. Args: m: instance of the Theory class g: index of DA Gegenbauer moment (often k in published papers) Returns: wc[s, k, p, f]: s in range (npws), k in range(npts), p in [LO, NLO], f in [Q,G,NSP] """ nebas = len(m.evolution_basis) # GPD, not DA wc = np.zeros((m.npws, m.npts, m.p + 1, nebas), dtype=complex) j = m.jpoints_pws fshu = _fshu(j) # Normalizations. Factor 3 is from normalization of DA, so not NC # See p. 37, 39 of "Towards DVMP" paper. Eq. (3.62c) # FIXME: NSP normalizations for DVMP not checked yet quark_norm = 3 * fshu gluon_norm = 3 * fshu * 2 / constants.CF / (j + 3) # Normalizations are chosen so that LO is normalized to: wc[:, :, 0, 0] = quark_norm * np.ones_like(j) / m.nf # Q wc[:, :, 0, 1] = gluon_norm * np.ones_like(j) # G wc[:, :, 0, 2] = quark_norm * np.ones_like(j) # NSP ?! if m.p == 1: qp1, ps1, g1 = c1dvmp.c1dvmp(m, 1, j, g) q1 = qp1/m.nf + ps1 nsp1 = qp1 wc[:, :, 1, 0] = quark_norm * q1 # Q wc[:, :, 1, 1] = gluon_norm * g1 # Q wc[:, :, 1, 2] = quark_norm * nsp1 # NSP ?! return wc
[docs] def calc_pj(m: theory.Theory, x: float, eta: float): """Calculate gegenbauer p_j to transform GPD from j to x space. Args: m: instance of the Theory class x, eta: GPD kinematical variables Returns: p_j[s, k, a]: s in range(npws), k in range(npts), a in [Q, G, ...] i.e. m.evolution_basis Todo: Check conventions for various quark combinations. """ # FIXME: vectorize p_j more universally so this looping is not needed pj = [] for j in m.jpoints_pws: pjf = [] rest_ebas = 0 # index of next item in evolution basis pjQ = gegenbauer.p_j(j, x, eta, 3/2) if 'G' in m.evolution_basis: pjG = - gegenbauer.p_j(j-1, x, eta, 5/2) # FIXME: sign?! pjf.append(pjQ) pjf.append(pjG) rest_ebas = 2 # we processed first two items rest_of_basis = m.evolution_basis[rest_ebas:] # may be empty if only SI neded for k, item in enumerate(rest_of_basis): pjf.append(pjQ) pj.append(np.stack(pjf, axis=1)) return np.stack(pj, axis=0)
[docs] def calc_wce(m: theory.Theory, Q2: float, process_class: str): """Calculate evolved Wilson coeffs for given Q2, for all PWs. Args: m: instance of the Theory Q2: process scale (so renormalization scale is Q2/m.rr2 and factorization scale is Q2/m.rf2) process_class: 'DIS' or 'DVCS' Returns: wce[s,k,b]: s in range(npwmax), k in range(npts), b in [Q,G,NSP] """ # alpha_strong/2pi at renormalization, factorization and input scales: asmur2 = qcd.as2pf(m.p, m.nf, Q2/m.rr2, m.asp[m.p], m.r20) asmuf2 = qcd.as2pf(m.p, m.nf, Q2/m.rf2, m.asp[m.p], m.r20) asQ02 = qcd.as2pf(m.p, m.nf, m.Q02, m.asp[m.p], m.r20) R = asmuf2/asQ02 wc = calc_wc(m, process_class) evola = evolution.evolop(m, R, process_class) p_mat = np.array([[1,],]) # LO version of p_mat, see below if m.p == 1: # p_mat: matrix that combines (LO, NLO) evolution operator and Wilson coeffs # while canceling NNLO term NLO*NLO: p_mat = np.array([[1, asmuf2], [asmur2, 0]]) # 3. evolved Wilson coeff. wce = np.einsum('skpa,pq,skqab->skb', wc, p_mat, evola) return wce[..., :3]
[docs] def calc_wce_dvmp(m: theory.Theory, Q2: float): """Calculate evolved DVMP Wilson coeffs for given Q2, for all PWs. Args: m: instance of the Theory Q2: default renormalization and factorization scale (can be changed with m.rdvmpr2, m.rf2 and m.rdaf2) Returns: wce[s,k,g,j]: s in range(npwmax), k in range(npts), g in gegens, j in [Q,G,NSP] Todo: Code duplication, merge with calc_wce() """ process_class = 'DVMP' # alpha_strong/2pi at various scales: asmur2 = qcd.as2pf(m.p, m.nf, Q2/m.rdvmpr2, m.asp[m.p], m.r20) # GPD asmuf2 = qcd.as2pf(m.p, m.nf, Q2/m.rf2, m.asp[m.p], m.r20) asQ02 = qcd.as2pf(m.p, m.nf, m.Q02, m.asp[m.p], m.r20) R = asmuf2/asQ02 # DA asdamuf2 = qcd.as2pf(m.p, m.nf, Q2/m.rdaf2, m.asp[m.p], m.r20) asdaQ02 = qcd.as2pf(m.p, m.nf, m.daQ02, m.asp[m.p], m.r20) RDA = asdamuf2/asdaQ02 # unfortunately c1dvmp doesn't accept numpy array of Gegenbauers # so we have to build wc like this: wc = [] for g in m.gpoints: wc.append(calc_wc_dvmp(m, g)) wc = np.stack(wc, axis=0) # gegenbauer index is first # FIXME: next line is NS^+ to get agreement with tests but rho is C=-1 evola_da = evolution.evolop_da(m, RDA, prty=+1) # +1->-1 !! evola = evolution.evolop(m, R, process_class) # p_mat: matrix that combines (LO, NLO) evolution operator and Wilson coeffs # while canceling NNLO term NLO*NLO: p_mat = np.array([1]).reshape((1, 1, 1)) # LO if m.p == 1: p_mat = np.zeros((2, 2, 2)) p_mat[0, 0, 0] = 1. p_mat[1, 0, 0] = asmur2 p_mat[0, 1, 0] = asmuf2 p_mat[0, 0, 1] = asdamuf2 # 3. evolved Wilson coeff. (evola for GPD, evola_ns for DA) wce = np.einsum('gskpi,pqr,skqij,gr->skgj', wc, p_mat, evola, evola_da) return wce
[docs] def calc_j2x(m: theory.Theory, x: float, eta: float, Q2: float): """Calculate j2x coeffs, combined with evolution operator. Args: m: instance of the Theory x: long. momentum fraction argument of GPD eta: skewness Q2: final evolution scale Returns: wce[s,k,i,j]: s in range(npwmax); k in range(npts); i,j in [Q,G,NSP] Todo: Implement general (eta != x) j to x transform. """ asmuf2 = qcd.as2pf(m.p, m.nf, Q2/m.rf2, m.asp[m.p], m.r20) asQ02 = qcd.as2pf(m.p, m.nf, m.Q02, m.asp[m.p], m.r20) R = asmuf2/asQ02 nebas = len(m.evolution_basis) wc = np.zeros((m.npws, m.npts, m.p + 1, nebas), dtype=complex) if eta < 1e-8: # forward limit, PDF-like one = np.ones_like(wc[:, :, 0, 0]) rest_ebas = 0 if 'G' in m.evolution_basis: wc[:, :, 0, 0] = one # Q wc[:, :, 0, 1] = x * one # G rest_ebas = 2 rest_of_basis = m.evolution_basis[rest_ebas:] for k, item in enumerate(rest_of_basis): wc[:, :, 0, rest_ebas + k] = one # NS evola = evolution.evolop(m, R, 'DIS') elif abs(eta-x) < 1e-8: # cross-over, border eta=x limit j = m.jpoints_pws fshu = _fshu(j) wc[:, :, 0, 0] = fshu # Q wc[:, :, 0, 1] = fshu*2*x/(3+j) # G # NSP is left to be zero before normalization can be checked! evola = evolution.evolop(m, R, 'DVCS') else: # FIXME: a bit of a kludge follows pj = np.squeeze(calc_pj(m, x, eta)) if pj.ndim == 2: wc[:, :, 0, 0] = pj / np.sin(np.pi * (m.jpoints_pws + 1)) else: wc[:, :, 0, :] = np.einsum('sja,sj->sja', pj, 1 / np.sin(np.pi * (m.jpoints_pws + 1))) evola = evolution.evolop(m, R, 'DVCS') p_mat = np.array([[1,],]) # LO if m.p == 1: p_mat = np.array([[1, asmuf2], [0, 0]]) # Return evolved Wilson coeff. # Note the difference w.r.t. usual wce, where we have inner product of # wc and evola, while here we have element-wise product, so that we # are in gpd.py able to get separately quark and gluon GPDs in x-space wce = np.einsum('skpa,pq,skqab->skab', wc, p_mat, evola) return wce