[NA] B2KKg decay

[mmikhasenko]

Reference

New LHCb paper on B2KKg is out

https://inspirehep.net/literature/2794646

Model content

Only resonances in KK,

Lineshapes

While the lineshape is tricky, it's kind of BW otimes Gauss`,

[mmikhasenko]

For validation, we'd need a few point (mKK, cos(theta)) for M for the original implementation.

mKK	cos(theta)	value
1.5	0.5	XXX
2.2	0.7	XXX
2.7	-0.1	XXX

[mmikhasenko]

A quick attempt to code fails to reproduce the paper plots

### A Pluto.jl notebook ###
# v0.19.42

using Markdown
using InteractiveUtils

# ╔═╡ 465dc2a2-2ce3-11ef-3449-59c4261d85ad
# ╠═╡ show_logs = false
begin
    using Pkg
    Pkg.activate()
    # 
    using DataFrames
    using ThreeBodyDecays
    using HadronicLineshapes
    using Plots
    using QuadGK
end

# ╔═╡ 757bf6cb-c729-4b70-91a6-2ba4212b4948
begin
    const mK = 0.458
    const mB = 5.4
end

# ╔═╡ b33bf712-90b3-4295-9a52-16a4ac68011a
begin
    ms = ThreeBodyMasses(mK,mK,0; m0=mB)
    two_js = ThreeBodySpins(0,0,2; two_h0=0)
    tbs = ThreeBodySystem(ms, two_js)
end; 

# ╔═╡ 7e257210-e8ec-47ae-8c20-3304c0398c7f
σs = x2σs([0.1, 0.2], ms; k=3)

# ╔═╡ e5b22100-c3e9-436e-b252-87ed18fb8f22
let two_j = 2
    dc = DecayChain(;
        k=3, two_j, tbs,
        Xlineshape=identity,
        Hij=NoRecoupling(0,0),
        HRk=ParityRecoupling(2,2,false))
    # 
    amplitude(dc, σs, (0,0,+2,0)) == amplitude(dc, σs, (0,0,-2,0))
end

# ╔═╡ 7c72fed9-f2b5-4de4-a5ff-2244c48687fb
df = [
    (j = 1, X = (m = 1.019, Γ = 4.249e-3), c=10*exp(0)), 
    (j = 2, X = (m = 1.5704, Γ = 86e-3), c=4.6*exp(0)),
    (j = 1, X = (m = 1.689, Γ = 211e-3), c=2.4*exp(137/180*π)),
    (j = 2, X = (m = 1.2755, Γ = 186.6e-3), c=1.07*exp(-55/180*π)),
    (j = 3, X = (m = 1.854, Γ = 78e-3), c=0.61*exp(-61/180*π)),
    (j = 2, X = (m = 2.011, Γ = 202e-3), c=0.74*exp(43/180*π)),
    (j = 1, X = σ->1.0, c=0.79)
] |> DataFrame

# ╔═╡ fa2964e8-61ff-42de-b264-10cd16efea32
df_l = transform(df, [:X, :j] => ByRow() do x, l
    !(x isa NamedTuple{(:m, :Γ)}) && return x
    scatt = BreitWigner(; x.m, x.Γ, ma=mK, mb=mK, l, d=3.0)
    ff_decay = BlattWeisskopf{l}(3.0)(breakup_ij(ms; k=3))
    ff_prod = BlattWeisskopf{l-1}(5.0)(breakup_Rk(ms; k=3))
    norm = ff_decay(x.m^2) * ff_prod(x.m^2) * sqrt(2l+1)
    scatt * ff_decay * ff_prod * (1/norm)
end => :X);

# ╔═╡ 9ddfaefb-1a39-4ba3-97dc-3f843e829e81
df_ch = transform(df_l, [:X, :j] => ByRow() do X, l
    DecayChain(;
        k=3, two_j=x2(l), tbs,
        Xlineshape=X,
        Hij=NoRecoupling(0,0),
        HRk=ParityRecoupling(2,2,false))
end => :chain);

# ╔═╡ 3b446004-609c-420c-b67e-31c2ff37242c
model = ThreeBodyDecay(
    "chain".*string.(1:size(df_ch,1)) .=> zip(df_ch.c, df_ch.chain));

# ╔═╡ f071e9ba-31ec-4927-8fbb-85b9c0bf4b52
unpolarized_intensity(model, σs)

# ╔═╡ 4c298f52-d0c4-4d25-be40-ef7082028dba
begin
    plot(1.05, 2.8) do mk
        I = Base.Fix1(unpolarized_intensity, model)
        integrand = projection_integrand(I, masses(model), mk^2; k = 3)
        mk * quadgk(integrand, 0, 1)[1]
    end
    for i in 1:length(model)
        plot!(1.05, 2.8) do mk
            I = Base.Fix1(unpolarized_intensity, model[i])
            integrand = projection_integrand(I, masses(model), mk^2; k = 3)
            mk * quadgk(integrand, 0, 1)[1]
        end
    end
    plot!(yscale=:log10)
end

# ╔═╡ Cell order:
# ╠═465dc2a2-2ce3-11ef-3449-59c4261d85ad
# ╠═757bf6cb-c729-4b70-91a6-2ba4212b4948
# ╠═b33bf712-90b3-4295-9a52-16a4ac68011a
# ╠═7e257210-e8ec-47ae-8c20-3304c0398c7f
# ╠═e5b22100-c3e9-436e-b252-87ed18fb8f22
# ╠═7c72fed9-f2b5-4de4-a5ff-2244c48687fb
# ╠═fa2964e8-61ff-42de-b264-10cd16efea32
# ╠═9ddfaefb-1a39-4ba3-97dc-3f843e829e81
# ╠═3b446004-609c-420c-b67e-31c2ff37242c
# ╠═f071e9ba-31ec-4927-8fbb-85b9c0bf4b52
# ╠═4c298f52-d0c4-4d25-be40-ef7082028dba