How to symmetrize 3-body amplitude?

[JuanBSLeite]

Dear experts,

I'm trying to implement an example of the D+ -> pi+ pi+ pi- decay, but I'm struggling with the symmetrization of the amplitude (s12 -> s02) and with setting the coefficients with non-default values. I'm using the code below. As you can see, BW is set only on pair 12 for rho(770)0, pair 02 was ignored. Do you know what I'm missing?

import qrules
from IPython.display import Math
from sympy import latex

reaction = qrules.generate_transitions(
    initial_state=("D+",[0]),
    final_state=["pi+", "pi+", "pi-"],
    allowed_intermediate_particles=["rho(770)0"],
    allowed_interaction_types=["e", "w"],
    formalism="helicity",
    topology_building="isobar"
)

import ampform
from ampform.dynamics import PhaseSpaceFactor
from ampform.dynamics.builder import RelativisticBreitWignerBuilder, create_non_dynamic_with_ff,create_non_dynamic

bw_builder = RelativisticBreitWignerBuilder(
    energy_dependent_width=True,
    form_factor=True,
    phsp_factor=PhaseSpaceFactor,  # optional
)

model_builder = ampform.HelicityAmplitudeBuilder(reaction)
#model_builder.adapter.permutate_registered_topologies()
model_builder.dynamics.assign("D+", create_non_dynamic_with_ff)

for name in reaction.get_intermediate_particles().names:
    model_builder.dynamics.assign(name, bw_builder)

model = model_builder.formulate() # it builds the model

from ampform.io import aslatex
from IPython.display import Math

Math(aslatex(model.amplitudes))

Thanks,
Juan

[redeboer]

Hi Juan, thanks for reporting! This indeed looks like a bug 🐛 Your code can be run here:
https://compwa.github.io/ampform-475/lab?path=Dp-to-pip-pip-pim.ipynb
Will get back to this tomorrow — it would be good to have this channel in the unit tests.

[redeboer]

But I do realize now this may be a problem in amplitudes with indistinguishable fermions in the final state, due to Wigner rotations (for now better handled through ampform-dpd) 🤔 Hard to think of a three-body decay like that though.

[JuanBSLeite]

Maybe you could check how laura++ framework deals with it.

https://arxiv.org/abs/1711.09854

But I do realize now this may be a problem in amplitudes with indistinguishable fermions in the final state, due to Wigner rotations (for now better handled through ampform-dpd) 🤔 Hard to think of a three-body decay like that though.

[JuanBSLeite]

Got it.

Yes the naming is a bit unintuitive. QRules only generates the transition with ρ 0 → π + ( 1 ) π − ( 2 ) , because it only considers intrinsic, invariant particle properties, not kinematic properties (although it does generate spin projections and LS values, which are not invariant -- see ComPWA/qrules#219). The amplitude for ρ 0 → π + ( 0 ) π − ( 2 ) is generated by AmpForm when formulating the amplitude, and it is added as an additional term to A 12 . Also, there is indeed a sign flip in some cases, see #476 (comment), so you cannot simply take the sum of all amplitudes. Some contributions have a negative sign.

[redeboer]

Maybe you could check how laura++ framework deals with it.

Would be interesting to see, yes 👌

Fixed version available in AmpForm v0.15.9:

pip install ampform==0.15.9

Lite notebook with new version available here:
https://compwa.github.io/ampform-475/lab?path=Dp-to-pip-pip-pim-fixed.ipynb

[JuanBSLeite]

Now, I'm going to move to build a complete description, including detection efficiency and background. Our idea is to test it with one data sample published using Run 1 LHCb data.

[JuanBSLeite]

Hi @redeboer , I think this fix was not propagated to model.kinematic_variables

it is affecting this part of the example using Ds+->pi+pi+pi-:

from tensorwaves.data import (
    IntensityDistributionGenerator,
    TFWeightedPhaseSpaceGenerator,
)

weighted_phsp_generator = TFWeightedPhaseSpaceGenerator(
    initial_state_mass=reaction.initial_state[-1].mass,
    final_state_masses={i: p.mass for i, p in reaction.final_state.items()},
)
data_generator = IntensityDistributionGenerator(
    domain_generator=weighted_phsp_generator,
    function=intensity_func,
    domain_transformer=helicity_transformer,
)
data_momenta = data_generator.generate(1_000_000, rng)

output:

KeyError                                  Traceback (most recent call last)
[/tmp/ipython-input-159059717.py](https://localhost:8080/#) in <cell line: 0>()
     13     domain_transformer=helicity_transformer,
     14 )
---> 15 data_momenta = data_generator.generate(1_000_000, rng)
     16 # pd.DataFrame({
     17 #     (k, label): np.transpose(v)[i]

3 frames
[/usr/local/lib/python3.12/dist-packages/tensorwaves/function/__init__.py](https://localhost:8080/#) in __call__(self, data)
     88 
     89     def __call__(self, data: DataSample) -> np.ndarray:
---> 90         args = [data[var_name] for var_name in self.argument_order]
     91         return self.function(*args)
     92 

KeyError: 'm_02'

[redeboer]

Could be a bug then 😯

data_generator = IntensityDistributionGenerator(
    domain_generator=weighted_phsp_generator,
    function=intensity_func,
    domain_transformer=helicity_transformer,  # <-- ?
)

Which keys does the helicity_transformer.functions dictionary have?

[JuanBSLeite]

Here is:

{'m_0': PositionalArgumentFunction(function=<PjitFunction of <function _lambdifygenerated at 0x7ec2a4162e80>>, argument_order=('p0', 'p1', 'p2')),
 'm_1': PositionalArgumentFunction(function=<PjitFunction of <function _lambdifygenerated at 0x7ec2a41358a0>>, argument_order=('p0', 'p1', 'p2')),
 'm_2': PositionalArgumentFunction(function=<PjitFunction of <function _lambdifygenerated at 0x7ec2a41209a0>>, argument_order=('p0', 'p1', 'p2')),
 'm_012': PositionalArgumentFunction(function=<PjitFunction of <function _lambdifygenerated at 0x7ec2a4122ac0>>, argument_order=('p0', 'p1', 'p2')),
 'm_12': PositionalArgumentFunction(function=<PjitFunction of <function _lambdifygenerated at 0x7ec2a4135bc0>>, argument_order=('p0', 'p1', 'p2')),
 'phi_0': PositionalArgumentFunction(function=<PjitFunction of <function _lambdifygenerated at 0x7ec2a41342c0>>, argument_order=('p0', 'p1', 'p2')),
 'phi_1^12': PositionalArgumentFunction(function=<PjitFunction of <function _lambdifygenerated at 0x7ec2a4161940>>, argument_order=('p0', 'p1', 'p2')),
 'theta_0': PositionalArgumentFunction(function=<PjitFunction of <function _lambdifygenerated at 0x7ec2a4134c20>>, argument_order=('p0', 'p1', 'p2')),
 'theta_1^12': PositionalArgumentFunction(function=<PjitFunction of <function _lambdifygenerated at 0x7ec2a4160ae0>>, argument_order=('p0', 'p1', 'p2'))}