Skip to content

Commit

Permalink
Lineshape cleanup (#69)
Browse files Browse the repository at this point in the history 
* replaced Bugg in json, cleaned qmd (#68)
* got rid of constant (#67)
  • Loading branch information
@mmikhasenko
mmikhasenko authored Jun 5, 2024
1 parent 50fabc0 commit d93a407
Show file tree
Hide file tree
Showing 6 changed files with 18 additions and 312 deletions.
10 changes: 5 additions & 5 deletions docs/julia/Manifest.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -276,11 +276,11 @@ version = "1.10.8"

[[deps.HadronicLineshapes]]
deps = ["Parameters", "QuadGK", "StaticArrays"]
git-tree-sha1 = "bad0c37beb3c481eaafbffc3df166852a9eb462f"
git-tree-sha1 = "f06145c91f2316bb9cc69dd4afa8ad51b3c2bb34"
repo-rev = "main"
repo-url = "https://github.com/mmikhasenko/HadronicLineshapes.jl"
uuid = "49c9d978-1f9d-4e96-a984-0a9783c0b9bf"
version = "0.3.3"
version = "0.3.4"

[[deps.HarfBuzz_jll]]
deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"]
Expand Down Expand Up @@ -898,12 +898,12 @@ uuid = "e6563dab-9ca1-5843-bde3-2ccf38d63843"
version = "0.11.3"

[[deps.ThreeBodyDecaysIO]]
deps = ["DataFrames", "HadronicLineshapes", "JSON", "MacroTools", "OrderedCollections", "Parameters", "ThreeBodyDecays"]
git-tree-sha1 = "84fbac13959f426213ae1bc5f49c46b17fab7e30"
deps = ["DataFrames", "HadronicLineshapes", "JSON", "MacroTools", "OrderedCollections", "Parameters", "Polynomials", "ThreeBodyDecays"]
git-tree-sha1 = "7a8326221b833ad250c8e50d497b5de190510f3d"
repo-rev = "main"
repo-url = "https://github.com/mmikhasenko/ThreeBodyDecaysIO.jl"
uuid = "418e7ecf-680e-4cb5-ad61-5e2f006aefac"
version = "0.4.1"
version = "0.4.2"

[[deps.TranscodingStreams]]
git-tree-sha1 = "a947ea21087caba0a798c5e494d0bb78e3a1a3a0"
Expand Down
11 changes: 0 additions & 11 deletions docs/julia/lb2pkg-lhcb-2765817.qmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,17 +42,6 @@ theme(:wong2, frame=:box, grid=false, minorticks=true,
Non-standard lineshapes are used to model resonances that do not conform to a simple `BreitWigner` distributions, or a `MultichannelBreitWigner` has to be defined explicitly.
The code below defines a new `NonResonant` lineshape, and its deserialization method. In this case this is just a constant.

```{julia}
struct Constant <: HadronicLineshapes.AbstractFlexFunc
end
(X::Constant)(σ::Number) = 1
import ThreeBodyDecaysIO: dict2instance
function dict2instance(::Type{Constant}, dict)
return NamedArgFunc(Constant(), [""])
end
```

## Deserialization of Objects to a Workspace

Model components are deserialized from a JSON file into computational objects within a workspace for further manipulation. First, functions representing lineshapes and form-factors are built. Following this, distributions are processed and added to the workspace.
Expand Down
30 changes: 0 additions & 30 deletions docs/julia/lc2ppik-lhcb-2683025.qmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,30 +36,6 @@ theme(:wong2, frame=:box, grid=false, minorticks=true,
lab="")
```

## Function definitions

The model contains non-standard lineshapes that do not conform to a simple `BreitWigner` or `MultichannelBreitWigner` distribution, so they have to be defined explicitly.
The code below defines a new lineshape and implements its deserialization method (specializing `dict2instance()` using [multiple dispatch](https://docs.julialang.org/en/v1/manual/methods/#man-method-specializations)).

```{julia}
struct BreitWignerWidthExpLikeBugg <: HadronicLineshapes.AbstractFlexFunc
m::Float64
Γ::Float64
γ::Float64
end
function (BW::BreitWignerWidthExpLikeBugg)(σ::Number)
mK = 0.493677
mπ = 0.13957018
σA = mK^2 - mπ^2 / 2
@unpack m, Γ, γ = BW
Γt = (σ - σA) / (m^2 - σA) * Γ * exp(-γ * σ)
1 / (m^2 - σ - 1im * m * Γt)
end
function ThreeBodyDecaysIO.dict2instance(::Type{BreitWignerWidthExpLikeBugg}, dict)
@unpack mass, width, slope, x = dict
return NamedArgFunc(BreitWignerWidthExpLikeBugg(mass, width, slope), [x])
end
```

## Deserialization of Objects to a Workspace

Expand Down Expand Up @@ -114,12 +90,6 @@ let
_parameters = array2dict(parameter_point["parameters"];
key="name", apply=v -> v["value"])
#
function label_diff(diff; levels=[1e-2, 1e-10])
_diff = abs(diff)
_diff < levels[2] && return '🟢'
_diff < levels[1] && return '🟡'
return '🔴'
end
computed_value = dist(_parameters)
#
tonumber(X::Number) = X
Expand Down
255 changes: 1 addition & 254 deletions docs/python/lc2ppik-lhcb-2683025.qmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
jupyter: python3
---

# $\Lambda_c^+ \to p K^- \pi^+$c
# $\Lambda_c^+ \to p K^- \pi^+$

::: {.callout-note appearance="simple"}
Model definition: [`lc2ppik-lhcb-2683025.json`](https://github.com/RUB-EP1/amplitude-serialization/blob/main/models/lc2ppik-lhcb-2683025.json).
Expand Down Expand Up @@ -38,7 +38,6 @@ from ampform_dpd import DefinedExpression
from ampform_dpd.decay import FinalStateID, State, ThreeBodyDecay
from ampform_dpd.dynamics import (
BreitWigner,
BuggBreitWigner,
ChannelArguments,
EnergyDependentWidth,
MultichannelBreitWigner,
Expand Down Expand Up @@ -211,79 +210,6 @@ L1405_Flatte = formulate_multichannel_breit_wigner(
Math(aslatex(L1405_Flatte))
```

#### Breit-Wigner with exponential

The model contains one lineshape function that is not standard, so we have to implement a custom propagator dynamics builder for this.

```{python}
#| code-fold: true
s, m0, Γ0, m1, m2, γ = sp.symbols("s m0 Gamma0 m1 m2 gamma", nonnegative=True)
expr = BuggBreitWigner(s, m0, Γ0, m1, m2, γ)
Math(aslatex({expr: expr.doit(deep=False)}))
```

```{python}
CHAIN_DEFS[18]
```

```{python}
get_function_definition("K700_BuggBW", MODEL_DEFINITION)
```

```{python}
def formulate_bugg_breit_wigner(
propagator: Propagator, resonance: str, model: ModelDefinition
) -> DefinedExpression:
function_definition = get_function_definition(propagator["parametrization"], model)
node = propagator["node"]
i, j = node
s = to_mandelstam_symbol(node)
mass = sp.Symbol(f"m_{{{resonance}}}", nonnegative=True)
width = sp.Symbol(Rf"\Gamma_{{{resonance}}}", nonnegative=True)
γ = sp.Symbol(Rf"\gamma_{{{resonance}}}", nonnegative=True)
m1 = to_mass_symbol(i)
m2 = to_mass_symbol(j)
final_state = get_final_state(model)
return DefinedExpression(
expression=BuggBreitWigner(s, mass, width, m1, m2, γ),
definitions={
mass: function_definition["mass"],
width: function_definition["width"],
m1: final_state[i].mass,
m2: final_state[j].mass,
γ: function_definition["slope"],
},
)
```

```{python}
CHAIN_18 = CHAIN_DEFS[18]
K700_BuggBW = formulate_bugg_breit_wigner(
propagator=CHAIN_18["propagators"][0],
resonance=to_latex(CHAIN_18["name"]),
model=MODEL_DEFINITION,
)
Math(aslatex(K700_BuggBW))
```

#### General propagator dynamics builder

```{python}
DYNAMICS_BUILDERS = {
"BreitWignerWidthExpLikeBugg": formulate_bugg_breit_wigner,
}
```

```{python}
#| code-fold: true
exprs = [
formulate_dynamics(CHAIN_DEFS[0], MODEL_DEFINITION, to_latex, DYNAMICS_BUILDERS),
formulate_dynamics(CHAIN_DEFS[18], MODEL_DEFINITION, to_latex, DYNAMICS_BUILDERS),
formulate_dynamics(CHAIN_DEFS[20], MODEL_DEFINITION, to_latex, DYNAMICS_BUILDERS),
]
Math(aslatex(exprs))
```

## Construct `AmplitudeModel`

### Unpolarized intensity
Expand Down Expand Up @@ -335,66 +261,6 @@ definitions = formulate_chain_amplitude(λ0, λ1, λ2, λ3, MODEL_DEFINITION, ch
Math(aslatex(definitions))
```

### Full amplitude model

```{python}
MODEL = formulate(
MODEL_DEFINITION,
additional_builders=DYNAMICS_BUILDERS,
cleanup_summations=True,
to_latex=to_latex,
)
MODEL.intensity
```

```{python}
#| code-fold: true
if "EXECUTE_NB" in os.environ:
selected_amplitudes = MODEL.amplitudes
else:
selected_amplitudes = {
k: v for i, (k, v) in enumerate(MODEL.amplitudes.items()) if i < 2
}
Math(aslatex(selected_amplitudes, terms_per_line=1))
```

```{python}
#| code-fold: true
Math(aslatex(MODEL.variables))
```

```{python}
#| code-fold: true
Math(aslatex({**MODEL.invariants, **MODEL.masses}))
```

## Numeric results

```{python}
intensity_expr = MODEL.full_expression.xreplace(MODEL.variables)
intensity_expr = intensity_expr.xreplace(MODEL.parameter_defaults)
```

```{python}
#| echo: false
free_symbols = intensity_expr.free_symbols
assert len(free_symbols) == 3
assert str(sorted(free_symbols, key=str)) == "[sigma1, sigma2, sigma3]"
```

```{python}
#| code-summary: Lambdify to numeric function
#| code-fold: true
intensity_funcs = {}
for s, s_expr in tqdm(MODEL.invariants.items()):
k = int(str(s)[-1])
s_expr = s_expr.xreplace(MODEL.masses).doit()
expr = perform_cached_doit(intensity_expr.xreplace({s: s_expr}))
func = perform_cached_lambdify(expr, backend="jax")
assert len(func.argument_order) == 2, func.argument_order
intensity_funcs[k] = func
```

### Validation

::: {.callout-warning}
Expand All @@ -417,122 +283,3 @@ checksum_points = {
}
checksum_points
```
```{python}
def sigma1_of_sigma2_cos_theta(sigma2, cos_theta, masses):
# Define the masses
m0, m1, m2, m3 = masses
def Kallen(x, y, z):
return (x**2 + y**2 + z**2 - 2*(x*y + y*z + z*x))
# Calculations
EE4sigma = (sigma2 + m1**2 - m3**2) * (sigma2 + m0**2 - m2**2)
p2q24sigma = Kallen(sigma2, m3**2, m1**2) * Kallen(m0**2, sigma2, m2**2)
# Calculate σi
sigma1_expr = m0**2 + m1**2 - (EE4sigma - sp.sqrt(p2q24sigma) * cos_theta) / (2 * sigma2)
return sigma1_expr
```
```{python}
#| code-fold: true
array = []
for distribution_name, checksum in checksums.items():
for point_name, expected in checksum.items():
parameters = checksum_points[point_name]
if len(parameters) < 6:
continue
s2 = parameters["m_31"] ** 2
cos_theta_31 = parameters["cos_theta_31"]
s1 = float(sigma1_of_sigma2_cos_theta(
s2,
cos_theta_31,
list(MODEL.masses.values())
))
print({"point_name": point_name, "s1": s1, "s2": s2})
# checked correct here, julia gives 0.7980703453578917
computed = intensity_funcs[3]({"sigma1": s1, "sigma2": s2})
status = "🟢" if computed == expected else "🔴"
array.append((distribution_name, point_name, computed, expected, status))
pd.DataFrame(array, columns=["Distribution", "Point", "Computed", "Expected", "Status"])
```

::: {.callout-warning}
See [ComPWA/ampform-dpd#133](https://github.com/ComPWA/ampform-dpd/issues/133).
:::

### Dalitz plot

```{python}
#| code-fold: true
i, j = (2, 1)
k, *_ = {1, 2, 3} - {i, j}
σk, σk_expr = list(MODEL.invariants.items())[k - 1]
Math(aslatex({σk: σk_expr}))
```

```{python}
#| code-fold: true
#| code-summary: Define meshgrid for Dalitz plot
resolution = 1_000
m = sorted(MODEL.masses, key=str)
x_min = float(((m[j] + m[k]) ** 2).xreplace(MODEL.masses))
x_max = float(((m[0] - m[i]) ** 2).xreplace(MODEL.masses))
y_min = float(((m[i] + m[k]) ** 2).xreplace(MODEL.masses))
y_max = float(((m[0] - m[j]) ** 2).xreplace(MODEL.masses))
x_diff = x_max - x_min
y_diff = y_max - y_min
x_min -= 0.05 * x_diff
x_max += 0.05 * x_diff
y_min -= 0.05 * y_diff
y_max += 0.05 * y_diff
X, Y = jnp.meshgrid(
jnp.linspace(x_min, x_max, num=resolution),
jnp.linspace(y_min, y_max, num=resolution),
)
dalitz_data = {
f"sigma{i}": X,
f"sigma{j}": Y,
}
```

```{python}
#| code-summary: Prepare parametrized numerical function
intensities = intensity_funcs[k](dalitz_data)
```

```{python}
#| echo: false
assert not jnp.all(jnp.isnan(intensities)), "All intensities are NaN"
```

```{python}
#| code-fold: true
#| code-summary: Dalitz plot is not yet correct
#| output: false
def get_decay_products(
decay: ThreeBodyDecay, subsystem_id: FinalStateID
) -> tuple[State, State]:
if subsystem_id not in decay.final_state:
msg = f"Subsystem ID {subsystem_id} is not a valid final state ID"
raise ValueError(msg)
return tuple(s for s in decay.final_state.values() if s.index != subsystem_id)
plt.rc("font", size=18)
I_tot = jnp.nansum(intensities)
normalized_intensities = intensities / I_tot
fig, ax = plt.subplots(figsize=(14, 10))
mesh = ax.pcolormesh(X, Y, normalized_intensities)
ax.set_aspect("equal")
c_bar = plt.colorbar(mesh, ax=ax, pad=0.01)
c_bar.ax.set_ylabel("Normalized intensity (a.u.)")
sigma_labels = {
i: Rf"$\sigma_{i} = M^2\left({' '.join(p.latex for p in get_decay_products(DECAY, i))}\right)$"
for i in (1, 2, 3)
}
ax.set_xlabel(sigma_labels[i])
ax.set_ylabel(sigma_labels[j])
plt.show()
```
Loading

0 comments on commit d93a407

Please sign in to comment.