/
angles.py
170 lines (159 loc) · 7.32 KB
/
angles.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
"""Formulate expressions for scattering and alignment angles."""
from __future__ import annotations
import sympy as sp
from ampform.kinematics.phasespace import Kallen
def formulate_scattering_angle(
state_id: int, sibling_id: int
) -> tuple[sp.Symbol, sp.acos]:
r"""Formulate the scattering angle in the rest frame of the resonance.
Compute the :math:`\theta_{ij}` scattering angle as formulated in `Eq (A1) in the
DPD paper <https://journals.aps.org/prd/pdf/10.1103/PhysRevD.101.034033#page=9>`_.
The angle is that between particle :math:`i` and spectator particle :math:`k` in the
rest frame of the isobar resonance :math:`(ij)`.
"""
if not {state_id, sibling_id} <= {1, 2, 3}:
msg = "Child IDs need to be one of 1, 2, 3"
raise ValueError(msg)
if {state_id, sibling_id} in { # pyright: ignore[reportUnnecessaryContains]
(2, 1),
(3, 2),
(1, 3),
}:
msg = f"Cannot compute scattering angle θ{state_id}{sibling_id}"
raise NotImplementedError(msg)
if state_id == sibling_id:
msg = f"IDs of the decay products cannot be equal: {state_id}"
raise ValueError(msg)
symbol = sp.Symbol(Rf"theta_{state_id}{sibling_id}", real=True)
spectator_id = next(iter({1, 2, 3} - {state_id, sibling_id}))
m0 = sp.Symbol("m0", nonnegative=True)
mi = sp.Symbol(f"m{state_id}", nonnegative=True)
mj = sp.Symbol(f"m{sibling_id}", nonnegative=True)
mk = sp.Symbol(f"m{spectator_id}", nonnegative=True)
σj = sp.Symbol(f"sigma{sibling_id}", nonnegative=True)
σk = sp.Symbol(f"sigma{spectator_id}", nonnegative=True)
theta = sp.acos(
(2 * σk * (σj - mk**2 - mi**2) - (σk + mi**2 - mj**2) * (m0**2 - σk - mk**2))
/ (sp.sqrt(Kallen(m0**2, mk**2, σk)) * sp.sqrt(Kallen(σk, mi**2, mj**2)))
)
return symbol, theta
def formulate_theta_hat_angle(
isobar_id: int, aligned_subsystem: int
) -> tuple[sp.Symbol, sp.acos]:
r"""Formulate an expression for :math:`\hat\theta_{i(j)}`."""
allowed_ids = {1, 2, 3}
if not {isobar_id, aligned_subsystem} <= allowed_ids:
msg = f"Child IDs need to be one of {', '.join(map(str, allowed_ids))}"
raise ValueError(msg)
symbol = sp.Symbol(Rf"\hat\theta_{isobar_id}({aligned_subsystem})", real=True)
if isobar_id == aligned_subsystem:
return symbol, sp.S.Zero
if (isobar_id, aligned_subsystem) in {(3, 1), (1, 2), (2, 3)}:
remaining_id = next(iter(allowed_ids - {isobar_id, aligned_subsystem}))
m0 = sp.Symbol("m0", nonnegative=True)
mi = sp.Symbol(f"m{isobar_id}", nonnegative=True)
mj = sp.Symbol(f"m{aligned_subsystem}", nonnegative=True)
σi = sp.Symbol(f"sigma{isobar_id}", nonnegative=True)
σj = sp.Symbol(f"sigma{aligned_subsystem}", nonnegative=True)
σk = sp.Symbol(f"sigma{remaining_id}", nonnegative=True)
theta = sp.acos(
(
(m0**2 + mi**2 - σi) * (m0**2 + mj**2 - σj)
- 2 * m0**2 * (σk - mi**2 - mj**2)
)
/ (sp.sqrt(Kallen(m0**2, mj**2, σj)) * sp.sqrt(Kallen(m0**2, σi, mi**2)))
)
return symbol, theta
_, theta = formulate_theta_hat_angle(aligned_subsystem, isobar_id)
return symbol, -theta
def formulate_zeta_angle( # noqa: C901, PLR0911, PLR0914
rotated_state: int,
aligned_subsystem: int,
reference_subsystem: int,
) -> tuple[sp.Symbol, sp.acos]:
r"""Formulate an expression for the alignment angle :math:`\zeta^i_{j(k)}`."""
zeta_symbol = sp.Symbol(
Rf"\zeta^{rotated_state}_{{{aligned_subsystem}({reference_subsystem})}}",
real=True,
)
if rotated_state == 0:
_, theta = formulate_theta_hat_angle(aligned_subsystem, reference_subsystem)
return zeta_symbol, theta
if reference_subsystem == 0:
_, zeta = formulate_zeta_angle(rotated_state, aligned_subsystem, rotated_state)
return zeta_symbol, zeta
if aligned_subsystem == reference_subsystem:
return zeta_symbol, sp.S.Zero
m0, m1, m2, m3 = sp.symbols("m:4", nonnegative=True)
σ1, σ2, σ3 = sp.symbols("sigma1:4", nonnegative=True)
if (rotated_state, aligned_subsystem, reference_subsystem) == (1, 1, 3):
cos_zeta_expr = (
2 * m1**2 * (σ2 - m0**2 - m2**2)
+ (m0**2 + m1**2 - σ1) * (σ3 - m1**2 - m2**2)
) / (sp.sqrt(Kallen(m0**2, m1**2, σ1)) * sp.sqrt(Kallen(σ3, m1**2, m2**2)))
return zeta_symbol, sp.acos(cos_zeta_expr)
if (rotated_state, aligned_subsystem, reference_subsystem) == (1, 2, 1):
cos_zeta_expr = (
2 * m1**2 * (σ3 - m0**2 - m3**2)
+ (m0**2 + m1**2 - σ1) * (σ2 - m1**2 - m3**2)
) / (sp.sqrt(Kallen(m0**2, m1**2, σ1)) * sp.sqrt(Kallen(σ2, m1**2, m3**2)))
return zeta_symbol, sp.acos(cos_zeta_expr)
if (rotated_state, aligned_subsystem, reference_subsystem) == (2, 2, 1):
cos_zeta_expr = (
2 * m2**2 * (σ3 - m0**2 - m3**2)
+ (m0**2 + m2**2 - σ2) * (σ1 - m2**2 - m3**2)
) / (sp.sqrt(Kallen(m0**2, m2**2, σ2)) * sp.sqrt(Kallen(σ1, m2**2, m3**2)))
return zeta_symbol, sp.acos(cos_zeta_expr)
if (rotated_state, aligned_subsystem, reference_subsystem) == (2, 3, 2):
cos_zeta_expr = (
2 * m2**2 * (σ1 - m0**2 - m1**2)
+ (m0**2 + m2**2 - σ2) * (σ3 - m2**2 - m1**2)
) / (sp.sqrt(Kallen(m0**2, m2**2, σ2)) * sp.sqrt(Kallen(σ3, m2**2, m1**2)))
return zeta_symbol, sp.acos(cos_zeta_expr)
if (rotated_state, aligned_subsystem, reference_subsystem) == (3, 3, 2):
cos_zeta_expr = (
2 * m3**2 * (σ1 - m0**2 - m1**2)
+ (m0**2 + m3**2 - σ3) * (σ2 - m3**2 - m1**2)
) / (sp.sqrt(Kallen(m0**2, m3**2, σ3)) * sp.sqrt(Kallen(σ2, m3**2, m1**2)))
return zeta_symbol, sp.acos(cos_zeta_expr)
if (rotated_state, aligned_subsystem, reference_subsystem) == (3, 1, 3):
cos_zeta_expr = (
2 * m3**2 * (σ2 - m0**2 - m2**2)
+ (m0**2 + m3**2 - σ3) * (σ1 - m3**2 - m2**2)
) / (sp.sqrt(Kallen(m0**2, m3**2, σ3)) * sp.sqrt(Kallen(σ1, m3**2, m2**2)))
return zeta_symbol, sp.acos(cos_zeta_expr)
if (rotated_state, aligned_subsystem, reference_subsystem) in { # Eq (A10)
(1, 2, 3),
(2, 3, 1),
(3, 1, 2),
}:
def create_symbols(i):
return sp.symbols(f"m{i} sigma{i}", nonnegative=True)
mi, σi = create_symbols(rotated_state)
mj, σj = create_symbols(aligned_subsystem)
mk, σk = create_symbols(reference_subsystem)
cos_zeta_expr = (
2 * mi**2 * (mj**2 + mk**2 - σi)
+ (σj - mi**2 - mk**2) * (σk - mi**2 - mj**2)
) / (sp.sqrt(Kallen(σj, mk**2, mi**2)) * sp.sqrt(Kallen(σk, mi**2, mj**2)))
return zeta_symbol, sp.acos(cos_zeta_expr)
if (rotated_state, aligned_subsystem, reference_subsystem) in {
(1, 3, 1),
(2, 1, 2),
(3, 2, 3),
(1, 1, 2),
(2, 2, 3),
(3, 3, 1),
(1, 3, 2),
(2, 1, 3),
(3, 2, 1),
}:
_, zeta = formulate_zeta_angle(
rotated_state, reference_subsystem, aligned_subsystem
)
return zeta_symbol, -zeta
msg = (
"No expression for"
f" ζ^{rotated_state}_{aligned_subsystem}({reference_subsystem})"
)
raise NotImplementedError(msg)