"""Deeply Virtual Compton Scattering (DVCS) observables."""
# Unlike DVMP and DIS, DVCS theory is separated into three files
# dvcs.py - cross-section and derived observables
# bmk.py - BMK formulas for interference, DVCS-squared and BH-squared terms
# cff.py - models for Compton Form Factors
from numpy import cos, pi, sin, sqrt
from . import cff, data, quadrature, theory
from .constants import GeV2nb, Mp2, alpha
from .kinematics import weight_BH
[docs]
class DVCS(theory.Theory):
"""DVCS observables base class.
Implements cross-section for electroproduction of polarized real photon,
on polarized nucleon target, as well as derived cross-sections and asymmetries.
Subclasses should, as their methods, implement expressions for squared
BH, interference and squared DVCS amplitudes:
- TBH2unp, TINTunp, TDVCS2unp (for unpolarized target)
- TBH2LP, TINTLP, TDVCS2LP (for longitudinally polarized target)
- TBH2TP, TINTTP, TDVCS2TP (for transversally polarized target)
Implemented subclasses are:
- BMK - From Belitsky, Mueller and Kirchner arXiv:hep-ph/0112108
- hotfixedBMK - simple valence xB region improvement by Dieter
- BM10ex
- BM10 - Belitsky, Mueller and Ji
- BM10tw2 - BM10 with higher twists set to zero
"""
[docs]
def PreFacSigma(self, pt: data.DataPoint) -> float:
"""Overall prefactor of 4-fold XS.
Args:
pt: instance of DataPoint
Returns:
Overal prefactor from Eq. (22) of BMK, times 2pi because of
proton Phi integration and times y/Q2 because of dy -> dQ2.
Converted to nanobarns.
"""
return alpha**3 * pt.xB * pt.y**2 / (8. * pi * pt.Q2**2 *
sqrt(1.+pt.eps2)) * GeV2nb
[docs]
def XS(self, pt: data.DataPoint, **kwargs) -> float:
"""Differential 4-fold e p --> e p gamma cross section.
Args:
pt: DataPoint instance
**kwargs: keyword arguments
Keyword args:
vars (dict): dict which can override pt kinematics, e.g., vars={'phi': 0}
flip (string, list): flip the sign of pt attribute,
e.g., flip='in1polarization'
Returns:
Cross section.
This is for calculation of fully beam and target polarized cross-section.
Most of the usually measured cross-sections and asymmetries are derived
from this one.
Todo:
* transversal target BMK/Trento switching
"""
# We create new instance of data.Datapoint, which can be manipulated with
# without messing up the original pt.
# Overriding pt kinematics with those from 'vars' keyword attribute
if 'vars' in kwargs:
kin = data.DataPoint(kindict=kwargs['vars'])
data._fill_kinematics(kin, old=pt)
else:
kin = pt.copy()
# Copy non-kinematical info
for atr in ['in1charge', 'in1polarizationvector', 'in1polarization',
'in2polarizationvector', 'in2polarization', 'in2particle',
'process', 'exptype', 'in1energy', 'in2energy', 'observable']:
if atr in pt:
setattr(kin, atr, getattr(pt, atr))
# Pre-calculate some kinematical functions
kin.prepare()
# Flipping spins and/or charges for asymmetries, if asked for
if 'flip' in kwargs and kwargs['flip']:
if isinstance(kwargs['flip'], list):
for item in kwargs['flip']:
setattr(kin, item, - getattr(pt, item))
else:
setattr(kin, kwargs['flip'], - getattr(pt, kwargs['flip']))
# Weighting the integrand by BH propagators, if asked for
if 'weighted' in kwargs and kwargs['weighted']:
wgh = weight_BH(kin)
else:
wgh = 1
# Finally, we build up the cross-section
# 1. -- unpolarized target part --
aux = self.TBH2unp(kin) + self.TINTunp(kin) + self.TDVCS2unp(kin)
if hasattr(pt, 'in2polarizationvector'):
# 2. -- longitudinally polarized target part --
if pt.in2polarizationvector == 'L':
aux += kin.in2polarization*(
self.TBH2LP(kin) + self.TINTLP(kin) + self.TDVCS2LP(kin))
elif pt.in2polarizationvector == 'T':
# 3. -- transversally polarized target part --
# We directly take cos(varphi) or sin(varphi) terms depending if
# varFTn is specified
if hasattr(pt, 'varFTn'):
# FIXME: should deal properly with Trento/BKM differences
# Also, this gives result in Trento convention now,
# and is completely O.K. for Trento asymmetries
# kin.varphi = (3*pt.varFTn+1)*pi/4. # pi for cos, -pi/2 for sin
# And this is for BMK
kin.varphi = (1-pt.varFTn)*pi/4. # 0 for cos, pi/2 for sin
aux += kin.in2polarization*(
self.TBH2TP(kin) + self.TINTTP(kin) + self.TDVCS2TP(kin))
elif pt.in2polarizationvector == 'U':
pass
else:
raise ValueError('in2polarizationvector must be L, T, or U!')
return wgh * self.PreFacSigma(kin) * aux
[docs]
def _XGAMMA_DVCS_t_Approx(self, pt):
"""Partial DVCS (gamma* p -> gamma p) cross section differential in t.
Args:
pt: instance of DataPoint
Returns:
Cross section. Approx. formula used in NPB10 paper.
"""
W2 = pt.W * pt.W
# Simplified formula used also in Fortran gepard code
ReH, ImH, ReE, ImE, ReHt, ImHt, ReEt, ImEt = self.m.cff(pt)
# 260.5... = 4 pi alpha^2 GeV^2 nb
res = 260.5633976788416 * W2 * (
(ImH**2 + ReH**2)
- pt.t/(4.*Mp2)*(ReE**2 + ImE**2)) / (
(W2 + pt.Q2) * (2.0 * W2 + pt.Q2)**2)
return res
[docs]
def _XGAMMA_DVCS_t_Ex(self, pt):
"""Partial DVCS (gamma* p -> gamma p) cross section differential in t.
Args:
pt: instance of DataPoint
**kwargs: keyword arguments
Returns:
Cross section
Todo:
Htilde is missing here. Small at low-x, but bigger than some terms that are retained.
"""
eps2 = 4. * pt.xB**2 * Mp2 / pt.Q2
if cff.HybridCFF in self.__class__.mro():
# For hybrid models we cannot ask for cff() since self gets misinterpreted
ReH = self.m.ReH(pt)
ImH = self.m.ImH(pt)
ReE = self.m.ReE(pt)
ImE = self.m.ImE(pt)
else:
ReH, ImH, ReE, ImE, ReHt, ImHt, ReEt, ImEt = self.m.cff(pt)
# 65.1... = pi alpha^2 GeV^2 nb
res = 65.14079453579676 * (pt.xB**2 / pt.Q2**2 / (1-pt.xB) / (2-pt.xB)**2 /
sqrt(1 + eps2) * (
4 * (1 - pt.xB) * (ImH**2 + ReH**2)
- pt.xB**2 * (ReE**2+ImE**2 + 2*ReE*ReH + 2*ImE*ImH)
- (2-pt.xB)**2 * pt.t/4/Mp2*(ImE**2+ReE**2)))
return res
# Choice of what is actually used:
_XGAMMA_DVCS_t = _XGAMMA_DVCS_t_Ex
[docs]
def XSintphi(self, pt: data.DataPoint, **kwargs) -> float:
"""XS integrated over azimuthal angle.
Args:
pt: instance of DataPoint
**kwargs: keyword arguments
Returns:
Cross section (unpolarized)
Todo:
This is unnecessary function. Should be changed to XUU.
"""
mem = pt.__dict__.pop('FTn', None)
pt.FTn = 0
res = self.XUU(pt, **kwargs)
# restore old value if needed
if mem:
pt.in2polarization = mem
return 2*pi*res
[docs]
def _phiharmonic(self, fun, pt, **kwargs):
"""Return fun evaluated for phi=pt.phi, or proper harmonic."""
if 'phi' in pt or ('vars' in kwargs
and 'phi' in kwargs['vars']):
return fun(pt, **kwargs)
elif 'FTn' in pt:
if pt.FTn < 0:
res = quadrature.Hquadrature(
lambda phi: fun(pt, vars={'phi': phi}, **kwargs)
* sin(-pt.FTn*phi), 0, 2*pi)
elif pt.FTn > 0:
res = quadrature.Hquadrature(
lambda phi: fun(pt, vars={'phi': phi}, **kwargs)
* cos(pt.FTn*phi), 0, 2*pi)
elif pt.FTn == 0:
res = quadrature.Hquadrature(
lambda phi: fun(pt, vars={'phi': phi}, **kwargs), 0, 2*pi)/2
else:
raise ValueError('FTn = {} looks wrong!'.format(str(pt.FTn)))
return res / pi
else:
raise ValueError(
'{} has neither azimuthal angle phi nor harmonic FTn defined!'.format(pt))
# Observables: cross-sections
[docs]
def _XUU(self, pt, **kwargs):
"""4-fold beam helicity-independent cross section."""
R = kwargs.copy()
R.update({'flip': 'in1polarization'})
return (self.XS(pt, **kwargs) + self.XS(pt, **R)) / 2
[docs]
def XUU(self, pt, **kwargs):
"""4-fold beam helicity-independent cross section."""
return self._phiharmonic(self._XUU, pt, **kwargs)
[docs]
def XUUw(self, pt, **kwargs):
"""Weighted 4-fold beam helicity-independent cross section."""
kwargs['weighted'] = True
return self._phiharmonic(self.XUU, pt, **kwargs)
[docs]
def _XLU(self, pt, **kwargs):
"""4-fold beam helicity-dependent cross section."""
R = kwargs.copy()
R.update({'flip': 'in1polarization'})
return (self.XS(pt, **kwargs) - self.XS(pt, **R)) / 2
[docs]
def XLU(self, pt, **kwargs):
"""4-fold beam helicity-dependent cross section."""
return self._phiharmonic(self._XLU, pt, **kwargs)
[docs]
def XLUw(self, pt, **kwargs):
"""Weighted 4-fold beam helicity-dependent cross section."""
kwargs['weighted'] = True
return self._phiharmonic(self.XLU, pt, **kwargs)
[docs]
def _XUD(self, pt, **kwargs):
"""Cross-section difference on (long. or trans.) polarized target."""
pol = kwargs.copy()
pol.update({'flip': 'in2polarization'})
o = self.XS(pt, **kwargs)
f = self.XS(pt, **pol)
return (o-f)/2
[docs]
def XUL(self, pt, **kwargs):
"""Cross-section difference on longitudinally polarized target."""
assert pt.in2polarizationvector == 'L'
return self._XUD(pt, **kwargs)
[docs]
def XUT(self, pt, **kwargs):
"""Cross-section difference on transversally polarized target."""
assert pt.in2polarizationvector == 'T'
return self._XUD(pt, **kwargs)
[docs]
def _XCLU(self, pt, **kwargs):
"""4-fold beam charge-spin cross section difference.
Args:
pt: instance of DataPoint
**kwargs: keyword arguments
Returns:
D_{CS,U}/2, as measured by COMPASS.
"""
R = kwargs.copy()
R.update({'flip': ['in1polarization', 'in1charge']})
return (self.XS(pt, **kwargs) - self.XS(pt, **R)) / 2
[docs]
def XCLU(self, pt, **kwargs):
"""4-fold beam charge-spin cross section difference.
Args:
pt: instance of DataPoint
**kwargs: keyword arguments
Returns:
D_{CS,U}/2, as measured by COMPASS. (NOT TESTED.)
"""
return self._phiharmonic(self._XCLU, pt, **kwargs)
[docs]
def _XCUU(self, pt, **kwargs):
"""4-fold beam charge-spin cross section sum.
Args:
pt: instance of DataPoint
**kwargs: keyword arguments
Returns:
S_{CS,U}/2, as measured by COMPASS.
"""
R = kwargs.copy()
R.update({'flip': ['in1polarization', 'in1charge']})
return (self.XS(pt, **kwargs) + self.XS(pt, **R)) / 2
[docs]
def XCUU(self, pt, **kwargs):
"""4-fold beam charge-spin cross section sum.
Args:
pt: instance of DataPoint
**kwargs: keyword arguments
Returns:
S_{CS,U}/2, as measured by COMPASS. (NOT TESTED.)
"""
return self._phiharmonic(self._XCUU, pt, **kwargs)
# Observables: asymmetries
[docs]
def _AC(self, pt, **kwargs):
"""Calculate beam charge asymmtery."""
chg = kwargs.copy()
chg.update({'flip': 'in1charge'})
o = self.XS(pt, **kwargs)
c = self.XS(pt, **chg)
return (o - c) / (o + c)
# optimized formula (don't calculate parts which cancel anyway)
# return self.TINTunp(pt, phi, 0, 1) / (
# self.TBH2unp(pt, phi) + self.TDVCS2unp(pt, phi) )
[docs]
def AC(self, pt, **kwargs):
"""Calculate beam charge asymmetry."""
res = self._phiharmonic(self._AC, pt, **kwargs)
return res
[docs]
def _ALU(self, pt, **kwargs):
"""Calculate beam spin asymmetry."""
return self._XLU(pt, **kwargs) / self._XUU(pt, **kwargs)
[docs]
def _ALUapprox(self, pt, **kwargs):
"""Calculate beam spin asymmetry.
Args:
pt: instance of DataPoint
**kwargs: keyword arguments
Returns:
Beam spin asymmetry A_LU.
For harmonics, fast approximate formula is used.
"""
if 'phi' in pt:
return self._ALU(pt, **kwargs)
elif 'FTn' in pt and pt.FTn == -1:
if 'vars' in kwargs:
kwargs['vars'].update({'phi': pi/2.})
else:
kwargs['vars'] = {'phi': pi/2.}
return self._ALU(pt, **kwargs)
else:
raise ValueError('[%s] has neither azimuthal angle phi\
nor harmonic FTn = -1 defined!' % pt)
[docs]
def _ALUexact(self, pt, **kwargs):
"""Calculate beam spin asymmetry."""
res = self._phiharmonic(self._ALU, pt, **kwargs)
return res
# Make a choice
ALU = _ALUexact
[docs]
def _TSA(self, pt, **kwargs):
"""Target spin asymmetry (transversal or longitudinal)."""
pol = kwargs.copy()
pol.update({'flip': 'in2polarization'})
o = self.XS(pt, **kwargs)
p = self.XS(pt, **pol)
return (o-p)/(o+p)
[docs]
def TSA(self, pt, **kwargs):
"""Calculate target spin asymmetry (AUL or AUT).
Args:
pt: instance of DataPoint
**kwargs: keyword arguments
Returns:
Target spin asymmetry.
Whether AUL or AUT is calculated is determined by the
pt.in2polarizationvector being 'L' or 'T'.
"""
return self._phiharmonic(self._TSA, pt, **kwargs)
[docs]
def AUL(self, pt, **kwargs):
"""Calculate longitudinal target spin asymmetry."""
assert pt.in2polarizationvector == 'L'
return self.TSA(pt, **kwargs)
[docs]
def AUT(self, pt, **kwargs):
"""Calculate transversal target spin asymmetry."""
assert pt.in2polarizationvector == 'T'
return self.TSA(pt, **kwargs)
[docs]
def _BTSA(self, pt, **kwargs):
"""Calculate beam-target spin asymmetry (BTSA)."""
# According to 1004.0177 Eq. (1.8)
bpol = kwargs.copy()
bpol.update({'flip': 'in1polarization'})
tpol = kwargs.copy()
tpol.update({'flip': 'in2polarization'})
both = kwargs.copy()
both.update({'flip': ['in1polarization', 'in2polarization']})
o = self.XS(pt, **kwargs)
p = self.XS(pt, **bpol)
t = self.XS(pt, **tpol)
b = self.XS(pt, **both)
return ((o+b) - (p+t)) / ((o+b) + (p+t))
[docs]
def BTSA(self, pt, **kwargs):
"""Calculate beam-target spin asymmetry (ALL or ALT).
Args:
pt: instance of DataPoint
**kwargs: keyword arguments
Returns:
Beam-target spin asymmetry.
Whether ALL or ALT is calculated is determined by the
pt.in2polarizationvector being 'L' or 'T'.
"""
return self._phiharmonic(self._BTSA, pt, **kwargs)
[docs]
def ALL(self, pt, **kwargs):
"""Calculate beam and longitudinal target spin asymmetry."""
assert pt.in2polarizationvector == 'L'
return self.BTSA(pt, **kwargs)
[docs]
def ALT(self, pt, **kwargs):
"""Calculate beam and transversal target spin asymmetry."""
assert pt.in2polarizationvector == 'T'
return self.BTSA(pt, **kwargs)
[docs]
def _CBTSA(self, pt, chargepar=-1, **kwargs):
"""Calculate double charge and spin asymmetries."""
# According to 1106.2990 Eq. (18)(chargepar=1) and (19)(chargepar=-1)
T_args = kwargs.copy()
T_args.update({'flip': 'in2polarization'})
B_args = kwargs.copy()
B_args.update({'flip': 'in1polarization'})
BT_args = kwargs.copy()
BT_args.update({'flip': ['in1polarization', 'in2polarization']})
C_args = kwargs.copy()
C_args.update({'flip': 'in1charge'})
CT_args = kwargs.copy()
CT_args.update({'flip': ['in1charge', 'in2polarization']})
CB_args = kwargs.copy()
CB_args.update({'flip': ['in1charge', 'in1polarization']})
CBT_args = kwargs.copy()
CBT_args.update({'flip': ['in1charge', 'in1polarization', 'in2polarization']})
o = self.XS(pt, **kwargs)
t = self.XS(pt, **T_args)
b = self.XS(pt, **B_args)
bt = self.XS(pt, **BT_args)
c = self.XS(pt, **C_args)
ct = self.XS(pt, **CT_args)
cb = self.XS(pt, **CB_args)
cbt = self.XS(pt, **CBT_args)
return ((o-t-b+bt) + chargepar*(c-ct-cb+cbt)) / (o+t+b+bt+c+ct+cb+cbt)
[docs]
def ALTI(self, pt, **kwargs):
"""Calculate {A_LT,I} as defined by HERMES 1106.2990 Eq. (19)."""
return self._phiharmonic(self._CBTSA, pt, **kwargs)
[docs]
def ALTBHDVCS(self, pt, **kwargs):
"""Calculate {A_LT,BHDVCS} as defined by HERMES 1106.2990 Eq. (18)."""
return self._phiharmonic(self._CBTSA, pt, chargepar=+1, **kwargs)
[docs]
def _ALUI(self, pt, **kwargs):
"""Beam spin asymmetry, interference part."""
pol = kwargs.copy()
pol.update({'flip': 'in1polarization'})
chg = kwargs.copy()
chg.update({'flip': 'in1charge'})
both = kwargs.copy()
both.update({'flip': ['in1polarization', 'in1charge']})
o = self.XS(pt, **kwargs)
p = self.XS(pt, **pol)
c = self.XS(pt, **chg)
b = self.XS(pt, **both)
return ((o-p) - (c-b)) / ((o+p) + (c+b))
[docs]
def ALUI(self, pt, **kwargs):
"""Calculate beam spin asymmetry, interference part."""
# As defined by HERMES 0909.3587 Eq. (2.2).
return self._phiharmonic(self._ALUI, pt, **kwargs)
[docs]
def _ALUDVCS(self, pt, **kwargs):
"""Calculate ALU as defined by HERMES 0909.3587 Eq. (2.3)."""
pol = kwargs.copy()
pol.update({'flip': 'in1polarization'})
chg = kwargs.copy()
chg.update({'flip': 'in1charge'})
both = kwargs.copy()
both.update({'flip': ['in1polarization', 'in1charge']})
o = self.XS(pt, **kwargs)
p = self.XS(pt, **pol)
c = self.XS(pt, **chg)
b = self.XS(pt, **both)
return ((o-p) + (c-b)) / ((o+p) + (c+b))
[docs]
def ALUDVCS(self, pt, **kwargs):
"""Calculate beam spin asymmetry, DVCS part."""
# As defined by HERMES 0909.3587 Eq. (2.3).
return self._phiharmonic(self._ALUDVCS, pt, **kwargs)
[docs]
def _AUTI(self, pt, **kwargs):
"""Calculate trans. target asymmetry, defined by HERMES 0802.2499 Eq. (15)."""
pol = kwargs.copy()
pol.update({'flip': 'in2polarization'})
chg = kwargs.copy()
chg.update({'flip': 'in1charge'})
both = kwargs.copy()
both.update({'flip': ['in2polarization', 'in1charge']})
o = self.XS(pt, **kwargs)
p = self.XS(pt, **pol)
c = self.XS(pt, **chg)
b = self.XS(pt, **both)
return ((o-p) - (c-b)) / ((o+p) + (c+b))
[docs]
def AUTI(self, pt, **kwargs):
"""Calculate trans. target asymmetry, defined by HERMES 0802.2499 Eq. (15)."""
return self._phiharmonic(self._AUTI, pt, **kwargs)
[docs]
def _AUTDVCS(self, pt, **kwargs):
"""Calculate trans. target asymmetry, defined by HERMES 0802.2499 Eq. (14)."""
pol = kwargs.copy()
pol.update({'flip': 'in2polarization'})
chg = kwargs.copy()
chg.update({'flip': 'in1charge'})
both = kwargs.copy()
both.update({'flip': ['in2polarization', 'in1charge']})
o = self.XS(pt, **kwargs)
p = self.XS(pt, **pol)
c = self.XS(pt, **chg)
b = self.XS(pt, **both)
return ((o-p) + (c-b)) / ((o+p) + (c+b))
[docs]
def AUTDVCS(self, pt, **kwargs):
"""Calculate trans. target asymmetry, defined by HERMES 0802.2499 Eq. (14)."""
return self._phiharmonic(self._AUTDVCS, pt, **kwargs)
# Observables: ad-hoc, one-off stuff
[docs]
def AC0minusr1(self, pt):
"""Return some combination liked by Dieter."""
return self.ACcos0(pt) - pt.r * self.ACcos1(pt)
[docs]
def XLUw2C(self, pt):
"""Im(C^I) as defined by HALL A."""
return self.ImCCALINTunp(pt)
[docs]
def XUUw2C(self, pt):
"""Re(C^I) or Re(C^I + Del C^I) as defined by HALL A.
Args:
pt: instance of DataPoint
Returns:
Harmonic coefficients.
Todo:
Although it is attributed to FTn=0, Re(C^I + Del C^I)
is only a part of zeroth harmonic.
"""
if pt.FTn == 0:
return self.ReCCALINTunp(pt) + self.ReDELCCALINTunp(pt)
elif pt.FTn == 1:
return self.ReCCALINTunp(pt)
[docs]
def XwA(self, pt):
"""Ratio of first two cos harmonics of w-weighted cross section."""
# In BMK, not Trento?
b0 = quadrature.Hquadrature(
lambda phi: self.XUU(pt, vars={'phi': phi}, weighted=True),
0, 2.0*pi) / (2.0*pi)
b1 = quadrature.Hquadrature(
lambda phi: self.XUU(pt, vars={'phi': phi}, weighted=True) * cos(phi),
0, 2.0*pi) / pi
return b1/b0
[docs]
def BCSA(self, pt, **kwargs):
"""Beam charge-spin asymmetry as measured by COMPASS."""
return self.BCSD(pt, **kwargs) / self.BCSS(pt, **kwargs)
# This PEP8 violating end-of-file import serves just to bring BMK into dvcs namespace
from .bmk import BM10, BMK, BM10ex, BM10tw2, hotfixedBMK # noqa: F401, E402