\(\Lambda_b^0 \to p K^- \gamma\)

Model definition: lb2pkg-lhcb-2765817.json.

This page demonstrates deserialization and evaluation of an amplitude model for the decay \(\Lambda_b^0 \to p K^- \gamma\). The resonant structure is studied using proton-proton collision data recorded at centre-of-mass energies of \(7\), \(8\), and \(13\) TeV collected with the LHCb detector, INSPIRE-HEP 2765817.

Activate environment 
using Pkg: 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
using LaTeXStrings

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

Function definitions

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.

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", "lb2pkg-lhcb-2765817.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, dist; workspace)
end

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^{-12}\), \(<10^{-2}\), or \(\ge10^{-2}\), respectively, for the difference between the reference and computed values.

A loop over validation points 
df = 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, Reference=value, Computed_value=computed_value,
            Point=point, Status=status)
    end |> DataFrame
end
17Γ—5 DataFrame
Row Distribution Reference Computed_value Point Status
String Any Complex… String Char
1 default_model 645.688 645.688+0.0im validation_point1 🟑
2 default_model 33.3131 33.3131+0.0im validation_point2 🟒
3 default_model 38.2867 38.2867+0.0im validation_point3 🟒
4 L1405_Flatte -0.740197550450631 + 0.23636822774611727i -0.740198+0.236368im validation_point_m12sq 🟒
5 L1890_BW 2.198693289424333 + 1.0366413279995204i 2.19869+1.03664im validation_point_m12sq 🟒
6 L1800_BW 0.3117344649326182 + 2.7741995686483376i 0.311734+2.7742im validation_point_m12sq 🟒
7 L1520_BW -0.6431712763866598 + 0.5539481728042333i -0.643171+0.553948im validation_point_m12sq 🟒
8 L2110_BW 0.8540683809737195 + 0.0415635681035755i 0.854068+0.0415636im validation_point_m12sq 🟒
9 L1680_BW -2.4613302416589518 + 0.36109808787606623i -2.46133+0.361098im validation_point_m12sq 🟒
10 L2350_BW 0.43057041569422266 + 0.00019710073244917788i 0.43057+0.000197101im validation_point_m12sq 🟒
11 L1600_BW -0.6293313194817324 + 0.7663369213404696i -0.629331+0.766337im validation_point_m12sq 🟒
12 L1810_BW 0.1064360891979254 + 5.093985975099388i 0.106436+5.09399im validation_point_m12sq 🟒
13 L1820_BW 4.647408096395898 + 4.443944768616481i 4.64741+4.44394im validation_point_m12sq 🟒
14 L1830_BW 3.824528421578999 + 3.827749839464702i 3.82453+3.82775im validation_point_m12sq 🟒
15 L2100_BW 0.8262246387007841 + 0.013532417439695257i 0.826225+0.0135324im validation_point_m12sq 🟒
16 L1690_BW -1.8597288296334098 + 1.396126484908448i -1.85973+1.39613im validation_point_m12sq 🟒
17 LNR30_NR 1.0 + 0.0i 1.0+0.0im validation_point_m12sq 🟒

Visualization

The model describing the decay is fetched from the workspace

model_dist = [v for (k, v) in workspace if v isa HadronicUnpolarizedIntensity] |> first;
model = model_dist.model
ThreeBodyDecay{74, DecayChain{X, ThreeBodySystem{@NamedTuple{m1::Float64, m2::Float64, m3::Float64, m0::Float64}, @NamedTuple{two_h1::Int64, two_h2::Int64, two_h3::Int64, two_h0::Int64}}} where X, ComplexF64, String}
  chains: StaticArraysCore.SVector{74, DecayChain{X, ThreeBodySystem{@NamedTuple{m1::Float64, m2::Float64, m3::Float64, m0::Float64}, @NamedTuple{two_h1::Int64, two_h2::Int64, two_h3::Int64, two_h0::Int64}}} where X}
  couplings: StaticArraysCore.SVector{74, ComplexF64}
  names: StaticArraysCore.SVector{74, String}

Dalitz plot

The Dalitz plot shows the probability distribution across two dimensional phase space of the decay in the range of mass \(m_{pK^{-}}\) from \(1.5\) to \(2.5 GeV/c^2\).

Dalitz plot plotting 
let
    ms = masses(model)
    #
    Οƒ3v = range(1.9, 2.5^2, 100)
    Οƒ2v = range(5, 27, 80)
    #
    _model = model
    f(Οƒs) = Kibble(Οƒs, ms^2) > 0 ? NaN : unpolarized_intensity(model, Οƒs)

    values = [f(Invariants(ms; Οƒ3, Οƒ2)) for Οƒ2 in Οƒ2v, Οƒ3 in Οƒ3v]
    heatmap(Οƒ3v, Οƒ2v, values, colorbar=false, c=:speed)
    plot!(border32(masses(model)),
        l=(3, :black),
        xlim=(Οƒ3v[1], Οƒ3v[end]), ylim=(Οƒ2v[1], Οƒ2v[end]),
        xlab=L"m_{pK^-}^2\,\, [\textrm{GeV}^2]",
        ylab=L"m_{p\gamma}^2\,\, [\textrm{GeV}^2]")
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.

Computation of projections 
let k = 3
    i, j = ij_from_k(k)
    model = model_dist.model
    #
    mlims = (sqrt(lims(3, masses(model))[1]), 2.5)
    mv = range(mlims..., 200) |> shift_by_half
    #
    plot()
    plot!(mv, lab=L"\textrm{Total}") do m
        I = Base.Fix1(unpolarized_intensity, model)
        m * quadgk(projection_integrand(I, masses(model), m^2; k), 0, 1)[1] / 1e3
    end
    chain_names = Set(model.names) |> collect |> sort
    for name in chain_names
        _model = model[model.names.==name]
        lab = replace(name, "L" => "\\varLambda(\\textrm{") * "})"
        plot!(mv, lab=latexstring(lab)) do m
            I = Base.Fix1(unpolarized_intensity, _model)
            m * quadgk(projection_integrand(I, masses(_model), m^2; k), 0, 1)[1] / 1e3
        end
    end
    plot!(;
        xlab=L"m_{pK^-}\,\, [\textrm{GeV}]",
        ylab=L"\textrm{d} N /\textrm{d} m_{pK^-}\,\, [\textrm{GeV}^{-1}]")
end