Complex phsp factor produces nonsensical intensity

[sebastianJaeger]

If the complex phasespace factor is used to model a resonance below the kinematic threshold, the resulting intensity has an extremely high and sharp peak at a (seemingly) random position within the kinematic borders. A minimal working example to generate a model with this problem is:

import qrules
from ampform.dynamics import PhaseSpaceFactorComplex
from ampform.dynamics.builder import RelativisticBreitWignerBuilder
from ampform import get_builder

reaction = qrules.generate_transitions(
    initial_state=("J/psi(1S)", [-1, +1]),
    final_state=["p", "p~", "eta"],
    allowed_intermediate_particles=[
        "N(1440)",                            
    ],
    allowed_interaction_types=["strong", "EM"],
    formalism="canonical-helicity",
    mass_conservation_factor=1,
)

model_builder = get_builder(reaction)

complex_bw_builder = RelativisticBreitWignerBuilder(form_factor=True, phsp_factor=PhaseSpaceFactorComplex)
for name in reaction.get_intermediate_particles().names:
    model_builder.set_dynamics(name, complex_bw_builder)
        
model = model_builder.formulate()

Attached is a simple plot of the distribution of events generated with this model.
Complex_PHSP_Factor_example.pdf

[redeboer]

Hi @sebastianJaeger, thanks for your report. I adapted your code snippet a bit so that it has exactly one resonance and extended it so that it generates momenta. (Note: it's not necessary to generate a hit-and-miss data sample to view the intensity distribution.)

Extended code snippet

Running on cfca442 and ComPWA/tensorwaves@e422bce:

import ampform
import matplotlib.pyplot as plt
import numpy as np
import qrules
from ampform.dynamics import PhaseSpaceFactorComplex
from ampform.dynamics.builder import RelativisticBreitWignerBuilder
from matplotlib import cm
from tensorwaves.data import (
    SympyDataTransformer,
    TFPhaseSpaceGenerator,
    TFUniformRealNumberGenerator,
)
from tensorwaves.function.sympy import create_parametrized_function

# Generate model
reaction = qrules.generate_transitions(
    initial_state=("J/psi(1S)", [-1, +1]),
    final_state=["p", "p~", "eta"],
    allowed_intermediate_particles=[
        "N(1440)+",
    ],
    allowed_interaction_types=["strong", "EM"],
    formalism="canonical-helicity",
    mass_conservation_factor=1,
)

model_builder = ampform.get_builder(reaction)
breit_wigner_builder = RelativisticBreitWignerBuilder(
    form_factor=True, phsp_factor=PhaseSpaceFactorComplex
)
for name in reaction.get_intermediate_particles().names:
    model_builder.set_dynamics(name, breit_wigner_builder)
model = model_builder.formulate()

# Generate phase space sample
rng = TFUniformRealNumberGenerator(seed=0)
phsp_generator = TFPhaseSpaceGenerator(
    initial_state_mass=reaction.initial_state[-1].mass,
    final_state_masses={
        i: p.mass for i, p in reaction.final_state.items()
    },
)
phsp_momenta = phsp_generator.generate(1_000_000, rng)
helicity_transformer = SympyDataTransformer.from_sympy(
    model.kinematic_variables, backend="jax"
)
phsp = helicity_transformer(phsp_momenta)

# Visualize intensity distribution
unfolded_expression = model.expression.doit()
intensity_func = create_parametrized_function(
    expression=unfolded_expression,
    parameters=model.parameter_defaults,
    backend="jax",
)

fig, ax = plt.subplots(figsize=(9, 4))
intensities = intensity_func(phsp)
ax.hist(
    np.array(phsp["m_02"].real),
    weights=np.array(intensities),
    bins=200,
    alpha=0.5,
    density=True,
)
ax.set_yscale("log")
ax.set_xlabel(R"$m_{p\eta}$ [GeV]")

resonances = sorted(
    reaction.get_intermediate_particles(),
    key=lambda p: p.mass,
)
evenly_spaced_interval = np.linspace(0, 1, len(resonances))
colors = [cm.rainbow(x) for x in evenly_spaced_interval]
for p, color in zip(resonances, colors):
    ax.axvline(x=p.mass, linestyle="dotted", label=p.name, color=color)
ax.legend()
plt.show()

Initially, I thought that something was wrong with the new Breit-Wigner builder (#206), but if the same script is run on a resonance above threshold (where PhaseSpaceFactorComplex is equivalent to the default PhaseSpaceFactor), the distribution look smooth.

Snippet to generate resonance above threshold

reaction = qrules.generate_transitions(
    initial_state=("J/psi(1S)", [-1, +1]),
    final_state=["gamma", "pi0", "pi0"],
    allowed_intermediate_particles=["f(0)(980)"],
    allowed_interaction_types=["strong", "EM"],
)

My guess is that this phase space factor simply does not make sense for decay products with unequal masses (here: N(1440)⁺ → proton eta), see #151. PhaseSpaceFactorComplex only makes sure that numbers are produced for the negative values under its square root, but the physics probably does not make sense. An idea may be to to investigate the model interactively as shown in this tutorial.

[redeboer]

Another observation: the peak seems to originate from the form factor. Same code as in #235 (comment) with RelativisticBreitWignerBuilder(form_factor=False, ...) results in:

Versions: cfca442 and ComPWA/tensorwaves@e422bce

%config InlineBackend.figure_formats = ['svg']
import ampform
import matplotlib.pyplot as plt
import numpy as np
import qrules
from ampform.dynamics import PhaseSpaceFactorComplex
from ampform.dynamics.builder import RelativisticBreitWignerBuilder
from matplotlib import cm
from tensorwaves.data import (
    SympyDataTransformer,
    TFPhaseSpaceGenerator,
    TFUniformRealNumberGenerator,
)
from tensorwaves.function.sympy import create_parametrized_function

# Generate model
reaction = qrules.generate_transitions(
    initial_state=("J/psi(1S)", [-1, +1]),
    final_state=["p", "p~", "eta"],
    allowed_intermediate_particles=["N(1440)+"],
    allowed_interaction_types=["strong", "EM"],
    mass_conservation_factor=1,
)

model_builder = ampform.get_builder(reaction)
breit_wigner_builder = RelativisticBreitWignerBuilder(
    form_factor=False,
    phsp_factor=PhaseSpaceFactorComplex,
)
for name in reaction.get_intermediate_particles().names:
    model_builder.set_dynamics(name, breit_wigner_builder)
model = model_builder.formulate()

# Generate phase space sample
rng = TFUniformRealNumberGenerator(seed=0)
phsp_generator = TFPhaseSpaceGenerator(
    initial_state_mass=reaction.initial_state[-1].mass,
    final_state_masses={i: p.mass for i, p in reaction.final_state.items()},
)
phsp_momenta = phsp_generator.generate(1_000_000, rng)
helicity_transformer = SympyDataTransformer.from_sympy(
    model.kinematic_variables, backend="jax"
)
phsp = helicity_transformer(phsp_momenta)

# Visualize intensity distribution
unfolded_expression = model.expression.doit()
intensity_func = create_parametrized_function(
    expression=unfolded_expression,
    parameters=model.parameter_defaults,
    backend="jax",
)

fig, ax = plt.subplots(figsize=(9, 4))
intensities = intensity_func(phsp)
ax.hist(
    np.array(phsp["m_02"].real),
    weights=np.array(intensities),
    bins=200,
    alpha=0.5,
    density=True,
)
ax.set_yscale("log")
ax.set_xlabel(R"$m_{p\eta}$ [GeV]")

resonances = sorted(
    reaction.get_intermediate_particles(),
    key=lambda p: p.mass,
)
evenly_spaced_interval = np.linspace(0, 1, len(resonances))
colors = [cm.rainbow(x) for x in evenly_spaced_interval]
for p, color in zip(resonances, colors):
    ax.axvline(x=p.mass, linestyle="dotted", label=p.name, color=color)
ax.legend()
plt.show()

This seems to be a bug, though, because setting form_factor=False also switches off the energy dependent width (where the phase space factor is of relevance), letting it fall back to a standard relativistic_breit_wigner().

model.expression.args[0].args[0].args[0].args[0].doit()

[redeboer]

Another observation: the peak seems to originate from the form factor.

Okay, I take that back.
#236 makes it possible to independently add a form factor and/or set an energy-dependent width. Playing around a bit with that and also making the final state masses equal (proton gets eta mass) gives me the impression that

1. the peak comes from the energy dependent width.
2. the peak becomes more extreme the closer the sub-threshold resonance is to threshold.
3. with PhaseSpaceFactorComplex, it doesn't matter whether the masses of the decay product are equal or not.

Source code

Note: this code could be simplified with scalar masses, which turns the final state masses into parameters.

pip install \
  git+https://github.com/ComPWA/ampform@3ed3ed5 \
  "tensorwaves[jax,pwa] @ git+https://github.com/ComPWA/tensorwaves@bea8cd9"

# %config InlineBackend.figure_formats = ['svg']
import ampform
import attrs
import matplotlib.pyplot as plt
import numpy as np
import qrules
from ampform.dynamics import PhaseSpaceFactorComplex
from ampform.dynamics.builder import RelativisticBreitWignerBuilder
from matplotlib import cm

from tensorwaves.data import (
    SympyDataTransformer,
    TFPhaseSpaceGenerator,
    TFUniformRealNumberGenerator,
)
from tensorwaves.function.sympy import create_parametrized_function

# Modify particle masses
PDG = qrules.load_pdg()
eta_mass = PDG["eta"].mass
proton_mass = PDG["eta"].mass
resonance_mass = 0.97 * (proton_mass + eta_mass)

particle_set = set(PDG)

particle = PDG["p"]
modified_particle = attrs.evolve(particle, mass=proton_mass)
particle_set.remove(particle)
particle_set.add(modified_particle)

particle = PDG["N(1440)+"]
modified_particle = attrs.evolve(particle, mass=resonance_mass)
particle_set.remove(particle)
particle_set.add(modified_particle)

PDG = qrules.ParticleCollection(particle_set)

# Create amplitude model
reaction = qrules.generate_transitions(
    initial_state=("J/psi(1S)", [+1]),
    final_state=["p", "p~", "eta"],
    allowed_intermediate_particles=["N(1440)+"],
    allowed_interaction_types=["strong", "EM"],
    particle_db=PDG,
    mass_conservation_factor=1,
    # max_angular_momentum=1,  # filter out L=2
)
assert len(reaction.get_intermediate_particles()) == 1

model_builder = ampform.get_builder(reaction)
breit_wigner_builder = RelativisticBreitWignerBuilder(
    phsp_factor=PhaseSpaceFactorComplex
)
for name in reaction.get_intermediate_particles().names:
    model_builder.set_dynamics(name, breit_wigner_builder)

breit_wigner_builder.energy_dependent_width = False
breit_wigner_builder.form_factor = False
model_standard_bw = model_builder.formulate()

breit_wigner_builder.energy_dependent_width = False
breit_wigner_builder.form_factor = True
model_standard_bw_with_ff = model_builder.formulate()

breit_wigner_builder.energy_dependent_width = True
breit_wigner_builder.form_factor = False
model_energy_dependent = model_builder.formulate()

breit_wigner_builder.energy_dependent_width = True
breit_wigner_builder.form_factor = True
model_energy_dependent_with_ff = model_builder.formulate()

# Generate phase space sample
rng = TFUniformRealNumberGenerator(seed=0)
phsp_generator = TFPhaseSpaceGenerator(
    initial_state_mass=reaction.initial_state[-1].mass,
    final_state_masses={i: p.mass for i, p in reaction.final_state.items()},
)
phsp_momenta = phsp_generator.generate(1_000_000, rng)
helicity_transformer = SympyDataTransformer.from_sympy(
    model_standard_bw.kinematic_variables, backend="jax"
)
phsp = helicity_transformer(phsp_momenta)
phsp = {k: v.real for k, v in phsp.items()}  # important for sqrt of Blatt-Weisskopf!

# Generate intensity distributions
func_standard_bw = create_parametrized_function(
    expression=model_standard_bw.expression.doit(),
    parameters=model_standard_bw.parameter_defaults,
    backend="jax",
)
func_standard_bw_with_ff = create_parametrized_function(
    expression=model_standard_bw_with_ff.expression.doit(),
    parameters=model_standard_bw_with_ff.parameter_defaults,
    backend="jax",
)
func_energy_dependent = create_parametrized_function(
    expression=model_energy_dependent.expression.doit(),
    parameters=model_energy_dependent.parameter_defaults,
    backend="jax",
)
func_energy_dependent_with_ff = create_parametrized_function(
    expression=model_energy_dependent_with_ff.expression.doit(),
    parameters=model_energy_dependent_with_ff.parameter_defaults,
    backend="jax",
)

# Plot 'em!
fig, ax = plt.subplots(figsize=(9, 4))

hist_kwargs = dict(
    bins=200,
    histtype="step",
    density=True,
)
ax.hist(
    np.array(phsp["m_02"].real),
    weights=np.array(func_standard_bw(phsp)),
    label="Standard Breit-Wigner",
    **hist_kwargs,
)
ax.hist(
    np.array(phsp["m_02"].real),
    weights=np.array(func_standard_bw_with_ff(phsp)),
    label="Standard Breit-Wigner with form factor",
    linestyle="dotted",
    **hist_kwargs,
)
ax.hist(
    np.array(phsp["m_02"].real),
    weights=np.array(func_energy_dependent(phsp)),
    label="Energy-dependent Breit-Wigner",
    **hist_kwargs,
)
ax.hist(
    np.array(phsp["m_02"].real),
    weights=np.array(func_energy_dependent_with_ff(phsp)),
    label="Energy-dependent Breit-Wigner with form factor",
    linestyle="dotted",
    **hist_kwargs,
)

ax.set_yscale("log")
ax.set_xlabel(R"$m_{p\eta}$ [GeV]")

resonances = sorted(
    reaction.get_intermediate_particles(),
    key=lambda p: p.mass,
)
evenly_spaced_interval = np.linspace(0, 1, len(resonances))
colors = [cm.rainbow(x) for x in evenly_spaced_interval]
for p, color in zip(resonances, colors):
    ax.axvline(x=p.mass, linestyle="dotted", label=p.name, color=color)
ax.legend()
plt.show()

N(1440)⁺ mass at 97% of proton+eta mass

N(1440)⁺ mass at 99% of proton+eta mass

N(1440)⁺ mass at 95% of 2x eta mass (proton mass set to eta mass)

[redeboer]

@sebastianJaeger, it may be worth looking into the new Resonances section of the PDG. Existing implementation is based on PDG 2020. For the energy dependent width, compare in particular:

• PDG2020, §Resonances, Eq. (49.21)
• PDG2021, §Resonances, Eq. (50.25)

See also existing implementation/formula for EnergyDependentWidth.

[redeboer]

This problem is probably fixed through #265: the narrow peak is probably an artifact of the phase space factor behaviour below m₁-m₂. Figure below shows the phase space factors that were implemented as of posting this bug report. Chew-Mandelstam (lower right) displays the proper behaviour below m₁-m₂ and was implemented in #265.

[redeboer]

PhaseSpaceFactorSWave #235 (comment) should display the correct behaviour. There is still a narrow peak, but this is expected behaviour (after some discussion with @mmikhasenko).

Sub-threshold Breit-Wigner comparison form this page, without phase space.

Note that, technically speaking, one should compute the phase space factor from the dispersion integral (TR-003) in the case of higher angular momenta (here: L=1), but using this phase space factor for S-waves at least ensures better behaviour of the phase space factor below threshold.

---

Modified snippet of #235 (comment)

Source code

pip install ampform==0.14.0 tensorwaves[jax,pwa]==0.4.4

# %config InlineBackend.figure_formats = ['svg']
import ampform
import attrs
import matplotlib.pyplot as plt
import numpy as np
import qrules
from ampform.dynamics import PhaseSpaceFactorSWave
from ampform.dynamics.builder import RelativisticBreitWignerBuilder
from matplotlib import cm

from tensorwaves.data import (
    SympyDataTransformer,
    TFPhaseSpaceGenerator,
    TFUniformRealNumberGenerator,
)
from tensorwaves.function.sympy import create_parametrized_function

# Modify particle masses
PDG = qrules.load_pdg()
eta_mass = PDG["eta"].mass
proton_mass = PDG["eta"].mass
resonance_mass = 0.97 * (proton_mass + eta_mass)

particle_set = set(PDG)

particle = PDG["p"]
modified_particle = attrs.evolve(particle, mass=proton_mass)
particle_set.remove(particle)
particle_set.add(modified_particle)

particle = PDG["N(1440)+"]
modified_particle = attrs.evolve(particle, mass=resonance_mass)
particle_set.remove(particle)
particle_set.add(modified_particle)

PDG = qrules.ParticleCollection(particle_set)

# Create amplitude model
reaction = qrules.generate_transitions(
    initial_state=("J/psi(1S)", [+1]),
    final_state=["p", "p~", "eta"],
    allowed_intermediate_particles=["N(1440)+"],
    allowed_interaction_types=["strong", "EM"],
    particle_db=PDG,
    mass_conservation_factor=1,
)
assert len(reaction.get_intermediate_particles()) == 1

model_builder = ampform.get_builder(reaction)
model_builder.scalar_initial_state_mass = True
model_builder.stable_final_state_ids = [0, 1, 2]
breit_wigner_builder = RelativisticBreitWignerBuilder(
    phsp_factor=PhaseSpaceFactorSWave   # <-- fixes the problem
)
for name in reaction.get_intermediate_particles().names:
    model_builder.set_dynamics(name, breit_wigner_builder)

breit_wigner_builder.energy_dependent_width = False
breit_wigner_builder.form_factor = False
model_standard_bw = model_builder.formulate()

breit_wigner_builder.energy_dependent_width = False
breit_wigner_builder.form_factor = True
model_standard_bw_with_ff = model_builder.formulate()

breit_wigner_builder.energy_dependent_width = True
breit_wigner_builder.form_factor = False
model_energy_dependent = model_builder.formulate()

breit_wigner_builder.energy_dependent_width = True
breit_wigner_builder.form_factor = True
model_energy_dependent_with_ff = model_builder.formulate()

# Generate phase space sample
rng = TFUniformRealNumberGenerator(seed=0)
phsp_generator = TFPhaseSpaceGenerator(
    initial_state_mass=reaction.initial_state[-1].mass,
    final_state_masses={i: p.mass for i, p in reaction.final_state.items()},
)
phsp_momenta = phsp_generator.generate(1_000_000, rng)
helicity_transformer = SympyDataTransformer.from_sympy(
    model_standard_bw.kinematic_variables, backend="jax"
)
phsp = helicity_transformer(phsp_momenta)
phsp = {k: v.real for k, v in phsp.items()}  # important for sqrt of Blatt-Weisskopf!

# Generate intensity distributions
func_standard_bw = create_parametrized_function(
    expression=model_standard_bw.expression.doit(),
    parameters=model_standard_bw.parameter_defaults,
    backend="jax",
)
func_standard_bw_with_ff = create_parametrized_function(
    expression=model_standard_bw_with_ff.expression.doit(),
    parameters=model_standard_bw_with_ff.parameter_defaults,
    backend="jax",
)
func_energy_dependent = create_parametrized_function(
    expression=model_energy_dependent.expression.doit(),
    parameters=model_energy_dependent.parameter_defaults,
    backend="jax",
)
func_energy_dependent_with_ff = create_parametrized_function(
    expression=model_energy_dependent_with_ff.expression.doit(),
    parameters=model_energy_dependent_with_ff.parameter_defaults,
    backend="jax",
)

# Plot 'em!
fig, ax = plt.subplots(figsize=(9, 4))

hist_kwargs = dict(
    bins=200,
    histtype="step",
    density=True,
)
ax.hist(
    np.array(phsp["m_02"].real),
    weights=np.array(func_standard_bw(phsp)),
    label="Standard Breit-Wigner",
    **hist_kwargs,
)
ax.hist(
    np.array(phsp["m_02"].real),
    weights=np.array(func_standard_bw_with_ff(phsp)),
    label="Standard Breit-Wigner with form factor",
    linestyle="dotted",
    **hist_kwargs,
)
ax.hist(
    np.array(phsp["m_02"].real),
    weights=np.array(func_energy_dependent(phsp)),
    label="Energy-dependent Breit-Wigner",
    **hist_kwargs,
)
ax.hist(
    np.array(phsp["m_02"].real),
    weights=np.array(func_energy_dependent_with_ff(phsp)),
    label="Energy-dependent Breit-Wigner with form factor",
    linestyle="dotted",
    **hist_kwargs,
)

ax.set_yscale("log")
ax.set_xlabel(R"$m_{p\eta}$ [GeV]")

resonances = sorted(
    reaction.get_intermediate_particles(),
    key=lambda p: p.mass,
)
evenly_spaced_interval = np.linspace(0, 1, len(resonances))
colors = [cm.rainbow(x) for x in evenly_spaced_interval]
for p, color in zip(resonances, colors):
    ax.axvline(x=p.mass, linestyle="dotted", label=p.name, color=color)
ax.legend()
plt.show()

N(1440)⁺ mass set to 97% of p+η mass

N(1440)⁺ mass at 95% of 2x eta mass (proton mass set to eta mass)

[redeboer]

• @sebastianJaeger First attempt: fall back to normal Breit-Wigner with create_relativistic_breit_wigner()
• @sebastianJaeger Second attempt: PhaseSpaceFactorSWave with RelativisticBreitWignerBuilder. Example in code snippet above Complex phsp factor produces nonsensical intensity #235 (comment)
• @redeboer Test dispersion integral for higher angular momenta in above distribution

[mmikhasenko]

Let me add to the discussion.

• Once you want to describe the state with the mass below the threshold you have to be careful. The standard single channel Breit-Wigner won't be correct approach.
• The most common case (including N1440) is when you describe a state predominately decaying to the channel A with the lower threshold (N->ppi), while you want to describe the channel B with the higher threshold (N->p eta). In that case, you should
  • Call the standard Breit-Wigner with the masses of p and pi.
  • Plot against m(ppi) make sure it looks right, the mass and width
  • multiple this |BW|^2 by the phase space of the channel B.

[redeboer]

Thanks for your feedback!

Once you want to describe the state with the mass below the threshold you have to be careful. The standard single channel Breit-Wigner won't be correct approach.

The analysis will for now remain an efficiency study, so it was decided that a normal BW suffices for now just to describe the distribution (without the resonance positions necessarily making sense).

The most common case (including N1440) is when you describe a state predominately decaying to the channel A with the lower threshold (N->ppi), while you want to describe the channel B with the higher threshold (N->p eta).

Yes, you mentioned. From an internal note:
N(1440) predominantly decays to pπ, not pη. This has effect on the energy-dependent width (see e.g. #235), through PDG2021 Equation (50.24-25): the coupling g for pπ dominates in this sum over the g for pη.
The way I understand it, is that you can use a pion mass instead of an eta mass in Eq. (50.28), because the coupling for pη in Eq. (50.25) is much smaller than that for pπ. This removes the need for a phase space factor with analytic continuation.

Your recommendations can be addressed once the analysis can move to an actual PWA stage.

[mmikhasenko]

Just one more practical thing:

• once you add p eta channel with whatever fraction, the first thing to look at is m(ppi) again and check the mass and the width, renormalize them. It is needed because the higher threshold will impact the peak position, as you well see.