\(\Lambda_c^+ \to p K^- \pi^+\)

Model definition: lc2ppik-lhcb-2683025.json.

This page demonstrates deserialization and evaluation of an amplitude model for the decay \(\Lambda_c^+ \to p K^- \pi^+\). The amplitude analysis is performed based on roughly half a million of \(\Lambda_c^{\pm}\) decay candidates by the LHCb collaboration, INSPIRE-HEP 2683025. Details on the mapped of the amplitude model onto the standard helicity formalism can be found in appendix of INSPIRE-HEP 2623821.

import Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using ThreeBodyDecaysIO
using ThreeBodyDecaysIO.ThreeBodyDecays
using ThreeBodyDecaysIO.HadronicLineshapes
using ThreeBodyDecaysIO.Parameters
using ThreeBodyDecaysIO.DataFrames
using ThreeBodyDecaysIO.JSON
using Measurements
using Statistics
using QuadGK
using Plots

theme(:wong2, frame=:box, grid=false, minorticks=true,
    guidefontvalign=:top, guidefonthalign=:right,
    foreground_color_legend=nothing,
    xlim=(:auto, :auto), ylim=(:auto, :auto),
    lab="")

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.

input = open(joinpath(@__DIR__, "..", "..", "models", "lc2ppik-lhcb-2683025.json")) do io
    JSON.parse(io)
end

workspace = Dict{String,Any}()

@unpack functions = input
for fn in functions
    @unpack name, type = fn
    instance_type = eval(Symbol(type))
    workspace[name] = dict2instance(instance_type, fn)
end

@unpack distributions = input
for dist in distributions
    @unpack name, type = dist
    instance_type = eval(Symbol(type))
    workspace[name] = dict2instance(instance_type, distributions[1]; workspace)
end

┌ Warning: parameters are not used in HadronicUnpolarizedIntensity
└ @ ThreeBodyDecaysIO ~/.julia/packages/ThreeBodyDecaysIO/bNt1V/src/HadronicUnpolarizedIntensity.jl:18

Validation

The integrity of the model is checked by validating the value of distributions at a few phase space points. The table lists the validation checks and their status. The marks “🟢”, “🟡”, and “🔴” indicate an accuracy of \(<10^{-10}\), \(<10^{-2}\), or \(\ge10^{-2}\), respectively, for the difference between the reference and computed values.

let
    @unpack misc, parameter_points = input
    @unpack amplitude_model_checksums = misc

    map(amplitude_model_checksums) do check_point_info
        @unpack point, value, distribution = check_point_info
        #
        # pull distribution
        dist = workspace[distribution]

        # pull correct parameter point
        parameter_points_dict = array2dict(parameter_points; key="name")
        # find the point in the list of points
        parameter_point = parameter_points_dict[point]
        # compute, compare
        _parameters = array2dict(parameter_point["parameters"];
            key="name", apply=v -> v["value"])
        #
        computed_value = dist(_parameters)
        #
        tonumber(X::Number) = X
        tonumber(X::String) = string2complex(X)
        reference_value = tonumber(value)
        status = label_diff(reference_value - computed_value)
        #
        (; Distribution=distribution,
            Point=point, Status=status)
    end |> DataFrame
end

11×3 DataFrame

Row	Distribution	Point	Status
	String	String	Char
1	default_model	validation_point	🟡
2	L1405_Flatte	validation_point_m31sq	🟢
3	L1690_BW	validation_point_m31sq	🟢
4	D1232_BW	validation_point_m12sq	🟢
5	L1520_BW	validation_point_m31sq	🟢
6	L1600_BW	validation_point_m31sq	🟢
7	L2000_BW	validation_point_m31sq	🟢
8	D1600_BW	validation_point_m12sq	🟢
9	D1700_BW	validation_point_m12sq	🟢
10	K892_BW	validation_point_m23sq	🟢
11	L1670_BW	validation_point_m31sq	🟢

Visualization

The model describing the decay is fetched from the workspace

model_dist = [v for (k, v) in workspace if v isa HadronicUnpolarizedIntensity] |> first;

Dalitz plot

The Dalitz plot shows the probability distribution across two dimensional phase space of the decay.

let iσx = 2, iσy = 1
    xlab = ((i, j) = ij_from_k(iσx);
    "m²($i$j) [GeV²]")
    ylab = ((i, j) = ij_from_k(iσy);
    "m²($i$j) [GeV²]")
    model = model_dist.model

    plot(masses(model), Base.Fix1(unpolarized_intensity, model);
        iσx, iσy, xlab, ylab)
end

The projection of the model onto a mass variable is shown by black line. Contributions from individual resonances are shown by the colored lines.

let k = 2
    i, j = ij_from_k(k)
    xlab = "m($i$j) [GeV]"
    model = model_dist.model
    #
    mlims = sqrt.(lims(k, masses(model)))
    mv = range(mlims..., 150) |> shift_by_half
    #
    plot()
    plot!(mv, lab="Total") do m
        I = Base.Fix1(unpolarized_intensity, model)
        m * quadgk(projection_integrand(I, masses(model), m^2; k), 0, 1)[1]
    end
    chain_names = Set(model.names) |> collect |> sort
    for name in chain_names
        _model = model[model.names.==name]
        plot!(mv, lab=name) do m
            I = Base.Fix1(unpolarized_intensity, _model)
            m * quadgk(projection_integrand(I, masses(_model), m^2; k), 0, 1)[1]
        end
    end
    plot!(; xlab)
end

Fit Fractions

The contribution of different resonances to the overall decay process is quantitatively assessed using numerical evaluation on a finite-size sample. Statistical uncertainty in the reported values only reflects the precision of these calculations. For a understanding of model uncertainties, refer to the original publications.

let
    model = model_dist.model
    ms = masses(model)
    #
    x2 = rand(10000, 2)
    data = map(eachslice(x2; dims=1)) do (x, y)
        σ1 = lims1(ms)[1] + x * diff(lims1(ms) |> collect)[1]
        σ2 = lims2(ms)[1] + y * diff(lims2(ms) |> collect)[1]
        σs = Invariants(ms; σ1, σ2)
    end
    filter!(data) do σs
        Kibble(σs, ms^2) < 0
    end
    #
    chain_names = Set(model.names) |> collect |> sort
    _int_i = map(chain_names) do name
        _intensities = unpolarized_intensity.(model[model.names.==name] |> Ref, data)
        _value = mean(_intensities)
        _err = sqrt(cov(_intensities, _intensities) / length(data))
        _value ± _err
    end
    _int0 = sum(unpolarized_intensity.(model |> Ref, data)) / length(data)
    ff = round.(_int_i ./ _int0 .* 100; digits=2)
    DataFrame("Resonance" => chain_names, "Fit Fraction [%]" => ff)
end

12×2 DataFrame

Row	Resonance	Fit Fraction [%]
	String	Measurem…
1	D1232	28.99±0.6
2	D1600	4.5±0.06
3	D1700	3.89±0.06
4	K1430	14.7±0.22
5	K700	3.04±0.04
6	K892	22.6±0.56
7	L1405	7.75±0.22
8	L1520	1.9±0.09
9	L1600	5.13±0.06
10	L1670	1.16±0.04
11	L1690	1.2±0.03
12	L2000	9.5±0.12