"""Definitions of Compton Form Factors."""
from __future__ import annotations
from math import log, pi
from typing import Dict, Union
import numpy as np
from . import data, gpd, model, quadrature, wilson # noqa: F401
[docs]
def _flatten(T):
"""Flatten the tuple."""
if not isinstance(T, tuple):
return (T,)
elif len(T) == 0:
return ()
else:
return _flatten(T[0]) + _flatten(T[1:])
[docs]
class CFF(model.ParameterModel):
"""Base class for CFFs.
CFFs are set to be zero here. Actual models are built by subclassing this.
"""
allCFFs = ['ImH', 'ReH', 'ImE', 'ReE', 'ImHt', 'ReHt', 'ImEt', 'ReEt']
allCFFsb = ['ImH', 'ReH', 'ImE', 'ReE', 'ImHt', 'ReHt', 'ImEt', 'ReEb']
allCFFeffs = ['ImHeff', 'ReHeff', 'ImEeff', 'ReEeff',
'ImHteff', 'ReHteff', 'ImEteff', 'ReEteff']
# allGPDs = []
[docs]
def __init__(self, **kwargs):
self.nf = kwargs.setdefault('nf', 4)
self.rr2 = kwargs.setdefault('rr2', 1) # ratio Q2/(DVCS renorm. scale squared)
# squared DVCS charge factors
if self.nf == 3:
qs = 2/9
qns = 1/9
else: # nf = 4
qs = 5/18
qns = 1/6
self.dvcs_charges = (qs, qs, qns)
super().__init__(**kwargs)
[docs]
def print_CFFs(self, pt: data.DataPoint, format: str = None):
"""Print values of CFFs at given kinematic point.
Args:
pt: DataPoint instance
format: 'mma' to get output pastable into Mathematica
"""
vals = [getattr(self, cff)(pt) for cff in self.allCFFs]
if format == 'mma':
s = "{" + 8*"%s -> %f, "
s = s[:-2] + "}"
else:
s = 8*"%4s = %5.2f\n"
print(s % _flatten(tuple(zip(self.allCFFs, vals))))
# Initial definition of all CFFs. All just return zero.
for name in allCFFs:
exec('def %s(self, pt): return 0.' % name)
for name in allCFFeffs:
exec('def %s(self, pt): return 0.' % name)
[docs]
def ReEb(self, pt: data.DataPoint):
"""Return Re(E-bar).
Args:
pt: DataPoint instance
Returns:
xi*ReEt
"""
return (pt.xi)*self.ReEt(pt)
[docs]
def cff(self, pt: data.DataPoint) -> np.ndarray:
"""Return array(ReH, ImH, ReE, ...) for kinematic point.
Args:
pt: DataPoint instance
Returns:
Eight CFFs (ReH, ImH, ReE, ..., ImEt).
This is generic top-class method that only gathers results
of dedicated ReH(), ImH(), etc. functions.
The idea is that this method can be implemented more efficiently
in subclases where components of CFFs can be calculated at the
same time or in parallel.
"""
return np.array([self.ReH(pt), self.ImH(pt), self.ReE(pt), self.ImE(pt),
self.ReHt(pt), self.ImHt(pt), self.ReEt(pt), self.ImEt(pt)])
[docs]
def is_within_model_kinematics(self, pt: data.DataPoint):
"""Check if kinematics is sensible.
Args:
pt: DataPoint instance
Returns:
True if pt kinematics is inside the model validity range.
"""
return ((1.5 <= pt.Q2 <= 5.) and
(pt.tm < min(1., pt.Q2/4)) and
(1e-3 < pt.xB < 0.5))
[docs]
class MellinBarnesCFF(CFF):
"""Class of models built by Mellin-Barnes integration."""
[docs]
def __init__(self, **kwargs) -> None:
self.tgj = np.tan(pi*self.jpoints/2.)
# wce[Q2] = wce[spw, j, a] - Wilson coeffs evolved; local to model instance
self.wce: Dict[float, np.ndarray] = {} # DVCS Wilson coefs.
super().__init__(**kwargs)
[docs]
def cff(self, pt: data.DataPoint) -> np.ndarray:
"""Return array(ReH, ImH, ReE, ...) for kinematic point."""
try:
wce_ar = self.wce[pt.Q2]
except KeyError:
# calculate it
wce_ar = wilson.calc_wce(self, pt.Q2, 'DVCS')
# memorize it for future
self.wce[pt.Q2] = wce_ar
# Evaluations depending on model parameters:
h_prerot = self.H(pt.xi, pt.t)
h = np.einsum('f,fa,ja->jf', self.dvcs_charges, self.frot, h_prerot)
e_prerot = self.E(pt.xi, pt.t)
e = np.einsum('f,fa,ja->jf', self.dvcs_charges, self.frot, e_prerot)
reh, imh, ree, ime = self._mellin_barnes_integral_HE(pt.xi, wce_ar, h, e)
return np.array([reh, imh, ree, ime, 0, 0, 0, 0])
[docs]
def ReH(self, pt: data.DataPoint) -> float:
"""Return Re(CFF H) for kinematic point."""
return self.cff(pt)[0]
[docs]
def ImH(self, pt: data.DataPoint) -> float:
"""Return Im(CFF H) for kinematic point."""
return self.cff(pt)[1]
[docs]
def ReE(self, pt: data.DataPoint) -> float:
"""Return Re(CFF E) for kinematic point."""
return self.cff(pt)[2]
[docs]
def ImE(self, pt: data.DataPoint) -> float:
"""Return Im(CFF E) for kinematic point."""
return self.cff(pt)[3]
[docs]
class DispersionCFF(CFF):
"""Use dispersion relations for real parts of CFFs.
methods: ReH, ReE, ReHt, ReEt, subtraction
Subclass should implement ansaetze for ImH, ImE, ImHt, ImEt
and subtraction. This class implements just dispersion integrals.
"""
[docs]
def __init__(self, **kwargs) -> None:
super().__init__(**kwargs)
[docs]
def dispargV(self, x, fun, pt):
"""Integrand of the dispersion integral (vector case).
Args:
x: integration variable, longitudinal momentum fraction
fun: function providing imaginary part of CFF
pt: DataPoint instance
Returns:
Integrand of dispersion integral.
With variable change x->x^(1/(1-ga))=u in
order to tame the singularity at x=0.
"""
ga = 0.9 # Nice value obtained by experimentation in Mathematica
u = x**(1./(1.-ga))
if hasattr(fun, '__self__'):
res = u**ga * (fun(pt, u) - fun(pt))
else: # unbound method
res = u**ga * (fun(self, pt, u) - fun(self, pt))
return (2.*u) / (pt.xi**2 - u**2) * res / (1.-ga)
[docs]
def dispargA(self, x, fun, pt):
"""Integrand of the dispersion integral (axial-vector case).
Args:
x: integration variable, longitudinal momentum fraction
fun: function providing imaginary part of CFF
pt: DataPoint instance
Returns:
Integrand of dispersion integral.
With variable change x->x^(1/(1-ga))=u in
order to tame the singularity at x=0.
"""
ga = 0.9 # Value same as for V-case (TODO: is this the best choice?)
u = x**(1./(1.-ga))
if hasattr(fun, '__self__'):
res = u**ga * (fun(pt, u) - fun(pt))
else: # unbound method
res = u**ga * (fun(self, pt, u) - fun(self, pt))
return (2. * pt.xi) / (pt.xi**2 - u**2) * res / (1.-ga)
[docs]
def subtraction(self, pt):
"""Subtraction constant."""
return 0 # default
[docs]
def _ReV(self, pt, imfun, subsign):
"""Real part of vector CFFs H and E.
Args:
pt: DataPoint instance
imfun: function providing imaginary part of CFF
subsign: sign of subtraction constant (-1 for H, +1 for E)
Returns:
Dispersion integral over imfun minus/plus subtraction constant.
"""
res = quadrature.PVquadrature(self.dispargV, 0, 1, (imfun, pt))
if hasattr(imfun, '__self__'):
pvpi = (res + log(pt.xi**2 / (1.-pt.xi**2)) * imfun(pt)) / pi
else:
pvpi = (res + log(pt.xi**2 / (1.-pt.xi**2)) * imfun(self, pt)) / pi
# P.V./pi - subtraction constant C/(1-t/MC^2)^2
return pvpi + subsign * self.subtraction(pt)
[docs]
def _ReA(self, pt, imfun):
"""Real part of axial vector CFFs Ht and Et.
Args:
pt: DataPoint instance
imfun: function providing imaginary part of CFF
Returns:
Dispersion integral over imfun.
"""
res = quadrature.PVquadrature(self.dispargA, 0, 1, (imfun, pt))
if hasattr(imfun, '__self__'):
pvpi = (res + log((1.+pt.xi)/(1.-pt.xi)) * imfun(pt))/pi
else: # unbound method
pvpi = (res + log((1.+pt.xi)/(1.-pt.xi)) * imfun(self, pt))/pi
return pvpi # this is P.V./pi
[docs]
def ReH(self, pt, imfun=None):
"""Real part of CFF H."""
if imfun is None:
imfun = self.ImH
return self._ReV(pt, imfun, -1)
[docs]
def ReE(self, pt, imfun=None):
"""Real part of CFF E."""
if imfun is None:
imfun = self.ImE
return self._ReV(pt, imfun, +1)
[docs]
def ReHt(self, pt, imfun=None):
"""Real part of CFF Ht."""
if imfun is None:
imfun = self.ImHt
return self._ReA(pt, imfun)
[docs]
def ReEt(self, pt, imfun=None):
"""Real part of CFF Et."""
if imfun is None:
imfun = self.ImEt
return self._ReA(pt, imfun)
[docs]
class PionPole(object):
"""Various options for pion-pole contribution."""
[docs]
def DMfixpole(self, pt):
"""Return fixed pion-pole as used by Dieter."""
pole = (2.2390424 * (1. - (1.7*(0.0196 - pt.t))/(1.
- pt.t/2.)**2))/((0.0196 - pt.t)*pt.xi)
if 'in2particle' in pt and pt.in2particle == 'n':
return -pole # neutron
else:
return pole # proton
[docs]
def DMfreepole(self, pt):
"""Return free pion-pole as proposed by Dieter."""
pole = (self.parameters['rpi'] * 2.16444 / (0.0196 - pt.t) / (1.
- pt.t/self.parameters['mpi2'])**2 / pt.xi)
if 'in2particle' in pt and pt.in2particle == 'n':
return -pole # neutron
else:
return pole # proton
[docs]
class DispersionFixedPoleCFF(DispersionCFF, PionPole):
"""Model for CFFs as in arXiv:0904.0458."""
[docs]
def __init__(self, **kwargs):
self.nf = kwargs.setdefault('nf', 4)
# initial values of parameters and limits on their values
self.add_parameters({'Nsea': 1.5, 'alS': 1.13, 'alpS': 0.15,
'mS2': 0.499849, 'rS': 1.0, 'bS': 2.0,
'Nv': 1.35, 'alv': 0.43, 'alpv': 0.85,
'mv2': 1.0, 'rv': 0.5, 'bv': 2.2,
'C': 7.0, 'mC2': 1.69,
'tNv': 0.0, 'tal': 0.43, 'talp': 0.85,
'tmv2': 7.29, 'trv': 6.0, 'tbv': 3.0})
self.add_parameters_limits({'bS': (0.4, 5.0), 'mv2': (0.16, 2.25),
'rv': (0., 8.), 'bv': (0.4, 5.),
'C': (-10., 10.), 'mC2': (0.16, 4.),
'tmv2': (0.16, 4.), 'trv': (0., 8.),
'tbv': (0.4, 5.)})
super().__init__(**kwargs)
[docs]
def subtraction(self, pt):
"""Dispersion relations subtraction constant."""
return self.parameters['C']/(1.-pt.t/self.parameters['mC2'])**2
[docs]
def ImH(self, pt: data.DataPoint, xi: Union[float, np.ndarray] = 0):
"""Imaginary part of CFF H."""
p = self.parameters # just a shortcut
# FIXME: The following solution is not elegant
if isinstance(xi, np.ndarray):
# function was called with third argument that is xi nd array
x = xi
elif xi != 0:
# function was called with third argument that is xi number
x = xi
else:
# xi should be taken from pt object
x = pt.xi
t = pt.t
twox = 2.*x / (1.+x)
onex = (1.-x) / (1.+x)
if 'in2particle' in pt and pt.in2particle == 'n':
chgfac = (1.*4./9. + 2./9.) # neutron
else:
chgfac = (2.*4./9. + 1./9.) # proton
val = (chgfac * p['Nv'] * p['rv'] * twox**(-p['alv']-p['alpv']*t) *
onex**p['bv'] / (1. - onex*t/(p['mv2'])))
sea = ((2./9.) * p['Nsea'] * p['rS'] * twox**(-p['alS']-p['alpS']*t) *
onex**p['bS'] / (1. - onex*t/(p['mS2']))**2)
return pi * (val + sea) / (1.+x)
[docs]
def ImHt(self, pt: data.DataPoint, xi: Union[float, np.ndarray] = 0):
"""Imaginary part of CFF Ht i.e. tilde{H}."""
p = self.parameters # just a shortcut
# FIXME: The following solution is not elegant
if isinstance(xi, np.ndarray):
# function was called with third argument that is xi nd array
x = xi
elif xi != 0:
# function was called with third argument that is xi number
x = xi
else:
# xi should be taken from pt object
x = pt.xi
t = pt.t
twox = 2.*x / (1.+x)
onex = (1.-x) / (1.+x)
try:
regge = (-p['tal']-p['talp']*t)
except KeyError:
# Old models take Regge trajectory params from H:
regge = (-p['alv']-p['alpv']*t)
if 'in2particle' in pt and pt.in2particle == 'n':
chgfac = (1.*4./9. + 2./9.) # neutron
else:
chgfac = (2.*4./9. + 1./9.) # proton
val = (chgfac * p['tNv'] * p['trv']
* twox**regge * onex**p['tbv'] / (1. - onex*t/(p['tmv2'])))
return pi * val / (1.+x)
[docs]
def ImE(self, pt: data.DataPoint, xi: Union[float, np.ndarray] = 0):
"""Imaginary part of CFF E."""
# Just changing function signature w.r.t. CFF
# to make it compatible for dispersion integral
return 0
[docs]
def ReEt(self, pt: data.DataPoint):
"""Instead of disp. rel. use pole formula."""
return self.DMfixpole(pt)
[docs]
class DispersionFreePoleCFF(DispersionFixedPoleCFF):
"""Model for CFFs as in arXiv:0904.0458. + free pion pole."""
[docs]
def __init__(self, **kwargs):
# Adding two extra parameters:
self.add_parameters({'rpi': 1.0, 'mpi2': 1.0})
self.add_parameters_limits({
'rpi': (-8, 8.),
'mpi2': (0.16, 16.)})
super().__init__(**kwargs)
[docs]
def ReEt(self, pt):
"""Instead of disp. rel. use pole formula."""
return self.DMfreepole(pt)
[docs]
class HybridCFF(MellinBarnesCFF, DispersionFixedPoleCFF, CFF):
"""Combines MB model for small xB and DR model for valence xB.
Todo:
Check general consistency for xi != pt.xi for imag CFFs.
"""
[docs]
def is_within_model_kinematics(self, pt):
"""Test kinematics of datapoint."""
# relaxing xBmin and removing Q2max
return ((1.5 <= pt.Q2) and
(pt.tm < min(1., pt.Q2/4)) and
(1e-5 < pt.xB < 0.65))
# FIXME: this below looks inconsistent generally for xi != pt.xi !!
[docs]
def ImH(self, pt, xi=0):
"""Imaginary part of CFF H."""
return MellinBarnesCFF.ImH(self, pt) + DispersionFixedPoleCFF.ImH(self, pt, xi)
[docs]
def ReH(self, pt):
"""Real part of CFF H."""
return MellinBarnesCFF.ReH(self, pt) + DispersionFixedPoleCFF.ReH(
self, pt, imfun=DispersionFixedPoleCFF.ImH)
[docs]
def ImE(self, pt, xi=0):
"""Imaginary part of CFF E."""
return MellinBarnesCFF.ImE(self, pt) + DispersionFixedPoleCFF.ImE(self, pt, xi)
[docs]
def ReE(self, pt):
"""Real part of CFF E."""
return MellinBarnesCFF.ReE(self, pt) + DispersionFixedPoleCFF.ReE(
self, pt, imfun=DispersionFixedPoleCFF.ImE)
# Tildes are not yet provided by MellinBarnesCFF model
[docs]
def ImHt(self, pt, xi=0):
"""Imaginary part of CFF Ht."""
return DispersionFixedPoleCFF.ImHt(self, pt, xi)
[docs]
def ReHt(self, pt):
"""Real part of CFF Ht."""
return DispersionFixedPoleCFF.ReHt(self, pt, imfun=DispersionFixedPoleCFF.ImHt)
[docs]
def ImEt(self, pt):
"""Imaginary part of CFF Et."""
return DispersionFixedPoleCFF.ImEt(self, pt)
[docs]
def ReEt(self, pt):
"""Real part of CFF Et. Provided by subclasses."""
return 0
[docs]
class HybridFixedPoleCFF(HybridCFF):
"""HybridCFF model with fixed pion pole model for ReEt."""
[docs]
def ReEt(self, pt):
"""Instead of disp. rel. use fixed pion pole formula."""
return self.DMfixpole(pt)
[docs]
class HybridFreePoleCFF(HybridCFF):
"""HybridCFF model with free pion pole model for ReEt."""
[docs]
def ReEt(self, pt):
"""Instead of disp. rel. use free pion pole formula."""
return self.DMfreepole(pt)
# This PEP8 violating end-of-file import serves just to bring GK code into cff namespace
from .gk import GoloskokovKrollCFF # noqa: F401, E402