Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def helicity_amplitude(x, spin):
"""
Helicity amplitude for a resonance in scalar-scalar state
x : cos(helicity angle)
spin : spin of the resonance
"""
if spin == 0:
return tf.complex(tfext.constant(1.), tfext.constant(0.))
elif spin == 1:
return tf.complex(x, tfext.constant(0.))
elif spin == 2:
return tf.complex((3. * x ** 2 - 1.) / 2., tfext.constant(0.))
elif spin == 3:
return tf.complex((5. * x ** 3 - 3. * x) / 2., tfext.constant(0.))
elif spin == 4:
return tf.complex((35. * x ** 4 - 30. * x ** 2 + 3.) / 8., tfext.constant(0.))
else:
raise ValueError("Illegal spin number.")
def relativistic_breit_wigner(m2, mres, wres):
"""
Relativistic Breit-Wigner
"""
if wres.dtype is ctype:
return 1. / (tfext.to_complex(mres ** 2 - m2) - tf.complex(tfext.constant(0.), mres) * wres)
if wres.dtype is zfit.settings.fptype:
return 1. / tf.complex(mres ** 2 - m2, -mres * wres)
return None
def hankel1(x):
if l == 0:
return tfext.constant(1.)
if l == 1:
return 1 + x ** 2
if l == 2:
x2 = x ** 2
return 9 + x2 * (3. + x2)
if l == 3:
x2 = x ** 2
return 225 + x2 * (45 + x2 * (6 + x2))
if l == 4:
x2 = x ** 2
return 11025. + x2 * (1575. + x2 * (135. + x2 * (10. + x2)))
def ComplexTwoBodyMomentum(md, ma, mb):
"""
Momentum of two-body decay products D->AB in the D rest frame.
Output value is a complex number, analytic continuation for the
region below threshold.
"""
return tf.sqrt(tf.complex((md ** 2 - (ma + mb) ** 2) * (md ** 2 - (ma - mb) ** 2) /
(4 * md ** 2), tfext.constant(0.)))
def ZemachTensor(m2ab, m2ac, m2bc, m2d, m2a, m2b, m2c, spin, cache=False):
"""
Zemach tensor for 3-body D->ABC decay
"""
z = None
if spin == 0:
z = tf.complex(tfext.constant(1.), tfext.constant(0.))
if spin == 1:
z = tf.complex(m2ac - m2bc + (m2d - m2c) * (m2b - m2a) / m2ab, tfext.constant(0.))
if spin == 2:
z = tf.complex((m2bc - m2ac + (m2d - m2c) * (m2a - m2b) / m2ab) ** 2 - 1. / 3. * (
m2ab - 2. * (m2d + m2c) + (m2d - m2c) ** 2 / m2ab) * (
m2ab - 2. * (m2a + m2b) + (m2a - m2b) ** 2 / m2ab), tfext.constant(0.))
if cache:
optimization.cacheable_tensors += [z]
return z
def WignerD(phi, theta, psi, j2, m2_1, m2_2):
"""
Calculate Wigner capital-D function.
phi,
theta,
psi : Rotation angles
j : spin (in units of 1/2, e.g. 1 for spin=1/2)
m1 and m2 : spin projections (in units of 1/2, e.g. 1 for projection 1/2)
"""
i = tf.complex(tfext.constant(0), tfext.constant(1))
m1 = m2_1 / 2.
m2 = m2_2 / 2.
return tf.exp(-i * tfext.to_complex(m1 * phi)) * tfext.to_complex(
Wignerd(theta, j2, m2_1, m2_2)) * tf.exp(-i * tfext.to_complex(m2 * psi))
import tensorflow as tf
from zfit.core import tfext
from zfit.physics import constants as const
# from Sec.2.5 of A. Bharucha, D. Straub and R. Zwicky (BSZ2015)
t_plus = tf.square(const.MB + const.MKst)
t_minus = tf.square(const.MB - const.MKst)
t_zero = t_plus * (1. - tf.sqrt(1. - t_minus / t_plus))
# Mass resonances of Table 3 of BSZ2015
MR_A0 = tfext.constant(5.366)
MR_T1 = tfext.constant(5.415)
MR_V = MR_T1
MR_T2 = tfext.constant(5.829)
MR_T23 = MR_T2
MR_A1 = MR_T2
MR_A12 = MR_T2
# function for the transformation q2->z for the parametrization of FF and real parametrization of H
def z(t, t_plus, t_zero):
zeta = ((tf.sqrt(t_plus - t) - tf.sqrt(t_plus - t_zero)) /
(tf.sqrt(t_plus - t) + tf.sqrt(t_plus - t_zero)))
return zeta
def SpinRotationAngle(pa, pb, pc, bachelor=2):
"""
Calculate the angle between two spin-quantisation axes for the 3-body D->ABC decay
aligned along the particle B and particle A.
pa, pb, pc : 4-momenta of the final-state particles
bachelor : index of the "bachelor" particle (0=A, 1=B, or 2=C)
"""
if bachelor == 2:
return tfext.constant(0.)
pboost = LorentzVector(-SpatialComponents(pb) / Scalar(TimeComponent(pb)), TimeComponent(pb))
if bachelor == 0:
pa1 = SpatialComponents(LorentzBoost(pa, pboost))
pc1 = SpatialComponents(LorentzBoost(pc, pboost))
return tf.acos(ScalarProduct(pa1, pc1) / Norm(pa1) / Norm(pc1))
if bachelor == 1:
pac = pa + pc
pac1 = SpatialComponents(LorentzBoost(pac, pboost))
pa1 = SpatialComponents(LorentzBoost(pa, pboost))
return tf.acos(ScalarProduct(pac1, pa1) / Norm(pac1) / Norm(pa1))
return None