P-vector fit comparison#
Studied decay#
FINAL_STATES: list[tuple[str, ...]] = [
["K0", "Sigma+", "p~"],
["eta", "p", "p~"],
]
REACTIONS: list[ReactionInfo] = [
qrules.generate_transitions(
initial_state="J/psi(1S)",
final_state=final_state,
allowed_intermediate_particles=["N"],
allowed_interaction_types=["strong"],
formalism="helicity",
particle_db=create_particle_database(),
)
for final_state in FINAL_STATES
]
Amplitude builder#
MODELS: list[HelicityModel] = []
for reaction in REACTIONS:
builder = ampform.get_builder(reaction)
builder.adapter.permutate_registered_topologies()
builder.config.scalar_initial_state_mass = True
builder.config.stable_final_state_ids = [0, 1, 2]
create_dynamics_symbol = DynamicsSymbolBuilder()
for resonance in reaction.get_intermediate_particles():
builder.set_dynamics(resonance.name, create_dynamics_symbol)
MODELS.append(builder.formulate())
del builder, reaction, resonance
\[\begin{split}\displaystyle \begin{array}{rcl}
A^{01}_{0, 0, - \frac{1}{2}, - \frac{1}{2}} &=& - C_{J/\psi(1S) \to {N_1(3/2^-)}_{+1/2} \overline{p}_{+1/2}; N_1(3/2^-) \to K^{0}_{0} \Sigma^{+}_{+1/2}} X_{J^P={\frac{3}{2}}^{+}, Q=+1} D^{1}_{0,0}\left(- \phi_{01},\theta_{01},0\right) D^{\frac{3}{2}}_{- \frac{1}{2},\frac{1}{2}}\left(- \phi^{01}_{0},\theta^{01}_{0},0\right) \\
&+& - C_{J/\psi(1S) \to {N_1(3/2^-)}_{+1/2} \overline{p}_{-1/2}; N_1(3/2^-) \to K^{0}_{0} \Sigma^{+}_{+1/2}} X_{J^P={\frac{3}{2}}^{+}, Q=+1} D^{1}_{0,1}\left(- \phi_{01},\theta_{01},0\right) D^{\frac{3}{2}}_{\frac{1}{2},\frac{1}{2}}\left(- \phi^{01}_{0},\theta^{01}_{0},0\right) \\
&+& - C_{J/\psi(1S) \to {N_1(3/2^-)}_{+3/2} \overline{p}_{+1/2}; N_1(3/2^-) \to K^{0}_{0} \Sigma^{+}_{+1/2}} X_{J^P={\frac{3}{2}}^{+}, Q=+1} D^{1}_{0,-1}\left(- \phi_{01},\theta_{01},0\right) D^{\frac{3}{2}}_{- \frac{3}{2},\frac{1}{2}}\left(- \phi^{01}_{0},\theta^{01}_{0},0\right) \\
\end{array}\end{split}\]
\[\begin{split}\displaystyle \begin{array}{cll}
X_{J^P={\frac{3}{2}}^{+}, Q=+1} \\
N_1(3/2^-) & m=1.82\text{ GeV} & \Gamma=0.6\text{ GeV} \\
\end{array}\end{split}\]
Dynamics parametrization#
Phasespace factor#
\[\begin{split}\displaystyle \begin{array}{rcl}
\rho^\mathrm{CM}_{m_{1},m_{2}}\left(s\right) &=& - 16 i \pi \Sigma\left(s\right) \\
\Sigma\left(s\right) &=& \frac{- \left(m_{1}^{2} - m_{2}^{2}\right) \left(- \frac{1}{\left(m_{1} + m_{2}\right)^{2}} + \frac{1}{s}\right) \log{\left(\frac{m_{1}}{m_{2}} \right)} + \frac{2 \log{\left(\frac{\left|{m_{1}^{2} + m_{2}^{2} + 2 \sqrt{s} q\left(s\right) - s}\right|}{2 m_{1} m_{2}} \right)} q\left(s\right)}{\sqrt{s}}}{16 \pi^{2}} \\
q\left(s\right) &=& \frac{\sqrt{\lambda\left(s, m_{1}^{2}, m_{2}^{2}\right)}}{2 \sqrt{s}} \\
\end{array}\end{split}\]
\(K\)-matrix formalism#
n_channels = len(REACTIONS)
I = sp.Identity(n_channels)
K = sp.MatrixSymbol("K", n_channels, n_channels)
P = sp.MatrixSymbol("P", n_channels, 1)
F = sp.MatrixSymbol("F", n_channels, 1)
rho = sp.MatrixSymbol("rho", n_channels, n_channels)
PARAMETERS_DEFAULTS = {}
for model in MODELS:
PARAMETERS_DEFAULTS.update(model.parameter_defaults)
del model
PARAMETERS_DEFAULTS = {
par: value
for par, value in PARAMETERS_DEFAULTS.items()
if not re.match(r"^m_\d+$", par.name)
}
\(K\)-matrix parametrization#
\[\begin{split}\displaystyle \left[\begin{matrix}\frac{\left(g_{N_1(3/2^-),0}\right)^{2}}{- m_{01}^{2} + \left(m_{N_1(3/2^-)}\right)^{2}} & \frac{g_{N_1(3/2^-),0} g_{N_1(3/2^-),1}}{- m_{01}^{2} + \left(m_{N_1(3/2^-)}\right)^{2}}\\\frac{g_{N_1(3/2^-),0} g_{N_1(3/2^-),1}}{- m_{01}^{2} + \left(m_{N_1(3/2^-)}\right)^{2}} & \frac{\left(g_{N_1(3/2^-),1}\right)^{2}}{- m_{01}^{2} + \left(m_{N_1(3/2^-)}\right)^{2}}\end{matrix}\right]\end{split}\]
\(P\)-vector parametrization#
\[\begin{split}\displaystyle \left[\begin{matrix}\frac{\beta_{N_1(3/2^-)} g_{N_1(3/2^-),0}}{- m_{01}^{2} + \left(m_{N_1(3/2^-)}\right)^{2}}\\\frac{\beta_{N_1(3/2^-)} g_{N_1(3/2^-),1}}{- m_{01}^{2} + \left(m_{N_1(3/2^-)}\right)^{2}}\end{matrix}\right]\end{split}\]
Phase space factor parametrization#
\[\begin{split}\displaystyle \left[\begin{matrix}\rho^\mathrm{CM}_{m_{0,0},m_{1,0}}\left(m_{01}^{2}\right) & 0\\0 & \rho^\mathrm{CM}_{m_{0,1},m_{1,1}}\left(m_{01}^{2}\right)\end{matrix}\right]\end{split}\]
\(F\)-vector construction#
Note
For some reason one has to leave out the multiplication of \(\rho\) by \(i\) within the calculation of the \(F\) vector
F = (I - sp.I * K * rho).inv() * P
F
\[\displaystyle \left(\mathbb{I} + - i K \rho\right)^{-1} P\]
F_vector = F.as_explicit()
parametrizations = {**K_expressions, **rho_expressions, **P_expressions}
F_exprs = F_vector.xreplace(parametrizations)
F_exprs[0].simplify(doit=False)
\[\displaystyle \frac{i \beta_{N_1(3/2^-)} g_{N_1(3/2^-),0}}{\left(g_{N_1(3/2^-),0}\right)^{2} \rho^\mathrm{CM}_{m_{0,0},m_{1,0}}\left(m_{01}^{2}\right) + \left(g_{N_1(3/2^-),1}\right)^{2} \rho^\mathrm{CM}_{m_{0,1},m_{1,1}}\left(m_{01}^{2}\right) - i m_{01}^{2} + i \left(m_{N_1(3/2^-)}\right)^{2}}\]
Create numerical functions#
FULL_EXPRESSIONS_FVECTOR = [
perform_cached_doit(MODELS_FVECTOR[i].expression).xreplace(
DYNAMICS_EXPRESSIONS_FVECTOR[i]
)
for i in range(n_channels)
]
INTENSITY_FUNCS_FVECTOR = [
create_parametrized_function(
expression=perform_cached_doit(FULL_EXPRESSIONS_FVECTOR[i]),
backend="jax",
parameters=MODELS_FVECTOR[i].parameter_defaults,
)
for i in range(n_channels)
]
Generate data#
Phase space sample#
HELICITY_TRANSFORMERS = [
SympyDataTransformer.from_sympy(model.kinematic_variables, backend="jax")
for model in MODELS_FVECTOR
]
PHSP = []
ε = 1e-8
for i in range(n_channels):
rng = TFUniformRealNumberGenerator(seed=0)
phsp_generator = TFPhaseSpaceGenerator(
initial_state_mass=REACTIONS[i].initial_state[-1].mass,
final_state_masses={it: p.mass for it, p in REACTIONS[i].final_state.items()},
)
phsp_momenta = phsp_generator.generate(100_000, rng)
phsp = HELICITY_TRANSFORMERS[i](phsp_momenta)
phsp = {k: v.real for k, v in phsp.items()}
phsp = {k: v + ε * 1j if re.match(r"^m_\d\d$", k) else v for k, v in phsp.items()}
PHSP.append(phsp)
Set parameters for toy model#
toy_parameters = {
R"\beta_{N_1(3/2^-)}": 1 + 0j,
R"m_{N_1(3/2^-)}": 1.71,
R"g_{N_1(3/2^-),0}": 3.2,
R"g_{N_1(3/2^-),1}": 1.5,
}
for func in INTENSITY_FUNCS_FVECTOR:
func.update_parameters(toy_parameters)
Toy data sample#
DATA = []
for i in range(n_channels):
weighted_phsp_generator = TFWeightedPhaseSpaceGenerator(
initial_state_mass=MODELS[i].reaction_info.initial_state[-1].mass,
final_state_masses={
i: p.mass for i, p in MODELS[i].reaction_info.final_state.items()
},
)
data_generator = IntensityDistributionGenerator(
domain_generator=weighted_phsp_generator,
function=INTENSITY_FUNCS_FVECTOR[i],
domain_transformer=HELICITY_TRANSFORMERS[i],
)
data_momenta = data_generator.generate(50_000, rng)
data = HELICITY_TRANSFORMERS[i](data_momenta)
DATA.append(data)
Perform fit#
Estimator definition#
class EstimatorSum(Estimator):
def __init__(self, estimators: Iterable[Estimator]) -> None:
self.__estimators = tuple(estimators)
def __call__(self, parameters: Mapping[str, ParameterValue]) -> float:
return sum(estimator(parameters) for estimator in self.__estimators)
def gradient(
self, parameters: Mapping[str, ParameterValue]
) -> dict[str, ParameterValue]:
raise NotImplementedError
combined_estimators = EstimatorSum(
UnbinnedNLL(
INTENSITY_FUNCS_FVECTOR[i],
data=DATA[i],
phsp=PHSP[i],
backend="jax",
)
for i in range(n_channels)
)
Initial parameters#
initial_parameters = {
R"m_{N_1(3/2^-)}": 1.9,
R"g_{N_1(3/2^-),0}": 2.8,
R"g_{N_1(3/2^-),1}": 1.6,
}
Optimize parameters#
minuit2 = Minuit2()
fit_result = minuit2.optimize(combined_estimators, initial_parameters)
assert fit_result.minimum_valid
fit_result
FitResult(
minimum_valid=True,
execution_time=2.1707212924957275,
function_calls=130,
estimator_value=-31725.142543423186,
parameter_values={
'm_{N_1(3/2^-)}': 1.7281071116110833,
'g_{N_1(3/2^-),0}': 3.3232991703411034,
'g_{N_1(3/2^-),1}': 1.4418247651741936,
},
parameter_errors={
'm_{N_1(3/2^-)}': 0.011169368106137329,
'g_{N_1(3/2^-),0}': 0.06518758826745548,
'g_{N_1(3/2^-),1}': 0.04383431405848117,
},
)
Fit quality check#
compute_aic_bic(fit_result)
(-63444.28508684637, -63417.82575199314)
initial |
fit result |
expected |
deviation |
|
---|---|---|---|---|
\(g_{N_1(3/2^-),0}\) |
2.8 |
3.32 |
3.2 |
3.9% |
\(g_{N_1(3/2^-),1}\) |
1.6 |
1.44 |
1.5 |
3.9% |
\(m_{N_1(3/2^-)}\) |
1.9 |
1.73 |
1.71 |
1.1% |