\(X\to \pi^-\pi^+\pi^-\)

Model definition: x2pipipi-compass-1391643.json.

This page demonstrates deserialization and evaluation the a partial wave analysis of diffractively produced \(3\pi\) system. The analysis is performed independently in bins of \(3\pi\) mass, and transferred momentum. The total number of bins is 100. For every bin, the model includes about 170 decay chains (about 88 waves, symmetrized for two \(\pi^+\pi^-\) pairs). The main paper describing details of the analysis is INSPIRE-HEP 1391643.

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", "x2pipipi-compass-1391643.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^{-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, computed_value, value, Status=status)
    end |> DataFrame
end

21×5 DataFrame

Row	Distribution	Point	computed_value	value	Status
	String	String	Complex…	Any	Char
1	compass_3pi_JP=1+_M=0_1540_1560	validation_point1	0.557753+0.0im	0.557753	🟢
2	compass_3pi_JP=1+_M=0_1540_1560	validation_point2	1.43625+0.0im	1.43625	🟢
3	compass_3pi_JP=1+_M=0_1540_1560	validation_point3	0.108395+0.0im	0.108395	🟢
4	compass_3pi_JP=1+_M=0_1540_1560	validation_point4	0.791847+0.0im	0.791847	🟢
5	compass_3pi_JP=1-_M=1_1540_1560	validation_point1	0.0452047+0.0im	0.0452047	🟢
6	compass_3pi_JP=1-_M=1_1540_1560	validation_point2	0.010651+0.0im	0.010651	🟢
7	compass_3pi_JP=1-_M=1_1540_1560	validation_point3	0.01857+0.0im	0.01857	🟢
8	compass_3pi_JP=1-_M=1_1540_1560	validation_point4	0.0900839+0.0im	0.0900839	🟢
9	compass_3pi_JP=2+_M=1_1540_1560	validation_point1	0.0170822+0.0im	0.0170822	🟢
10	compass_3pi_JP=2+_M=1_1540_1560	validation_point2	0.0579535+0.0im	0.0579535	🟢
11	compass_3pi_JP=2+_M=1_1540_1560	validation_point3	0.035384+0.0im	0.035384	🟢
12	compass_3pi_JP=2+_M=1_1540_1560	validation_point4	0.112218+0.0im	0.112218	🟢
13	compass_3pi_JP=4+_M=1_1540_1560	validation_point1	0.0205558+0.0im	0.0205558	🟢
14	compass_3pi_JP=4+_M=1_1540_1560	validation_point2	0.0118184+0.0im	0.0118184	🟢
15	compass_3pi_JP=4+_M=1_1540_1560	validation_point3	0.0175504+0.0im	0.0175504	🟢
16	compass_3pi_JP=4+_M=1_1540_1560	validation_point4	0.105833+0.0im	0.105833	🟢
17	R(1274)	validation_point_m12sq	1.68268+0.620869im	1.682682255119889 + 0.62086927105196i	🟢
18	R(1690)	validation_point_m12sq	0.564308+0.103183im	0.5643082293246228 + 0.10318283266108788i	🟢
19	R(600)	validation_point_m12sq	0.101724+1.02734im	0.10172436035046228 + 1.0273440332286132i	🟢
20	R(768)	validation_point_m12sq	-1.73218+0.632399im	-1.7321826562432676 + 0.6323990884902008i	🟢
21	R(965)	validation_point_m12sq	-0.921258+2.14704im	-0.9212576583634395 + 2.1470398931994485i	🟢

Visualization

The model describing the decay is fetched from the workspace

model_names = [k for (k, v) in workspace if v isa HadronicUnpolarizedIntensity]

4-element Vector{String}:
 "compass_3pi_JP=1+_M=0_1540_1560"
 "compass_3pi_JP=4+_M=1_1540_1560"
 "compass_3pi_JP=2+_M=1_1540_1560"
 "compass_3pi_JP=1-_M=1_1540_1560"

Dalitz plot

The Dalitz plot shows the probability distribution across two dimensional phase space of the decay. Below the distribution is shown for all models in the file

let iσx = 3, iσy = 1
    n_models = length(model_names)
    Nx, Ny = 2, div(n_models + 1, 2)
    plot(
        map(model_names) do model_name
            model = workspace[model_name].model
            #
            xlab = ((i, j) = ij_from_k(iσx);
            "m²($i$j) [GeV²]")
            ylab = ((i, j) = ij_from_k(iσy);
            "m²($i$j) [GeV²]")

            plot(masses(model), Base.Fix1(unpolarized_intensity, model);
                iσx, iσy, xlab, ylab, title=model_name[13:end], aspect_ratio=1)
        end...,
        layout=grid(Ny, Nx),
        size=(350 * Nx, 350 * Ny))
end