2023-11-08-JuliaHEP-ThreeBodyDecay.jl

gistfile1.txt

### A Pluto.jl notebook ###
# v0.19.30

using Markdown
using InteractiveUtils

# ╔═╡ f4bbad10-7e14-11ee-0343-7da1e5ed0283
# ╠═╡ show_logs = false
begin
	cd(mktempdir())
	import Pkg
	Pkg.activate(".")
	Pkg.add([
		Pkg.PackageSpec(url="https://github.com/mmikhasenko/ThreeBodyDecay.jl"),
		Pkg.PackageSpec(url="https://github.com/mmikhasenko/SymbolicThreeBodyDecays.jl"),
		Pkg.PackageSpec("Plots")
	])
	# 
	using ThreeBodyDecay
	using SymbolicThreeBodyDecays
	using SymbolicThreeBodyDecays.SymPy
	# 
	using Plots
	import Plots.PlotMeasures: mm
end

# ╔═╡ a5d25226-f465-428f-931b-64825963c255
md"""
# Angular Analysis with Julia
"""

# ╔═╡ b4c00f7e-ce88-4516-a43a-d63e534bad6e
md"""
Hadronic decay modeling, Based on the **Dalitz Plot decomposition**, [PRD](https://inspirehep.net/literature/1758460)

In Julia:

- [ThreeBodyDecay.jl](https://github.com/mmikhasenko/ThreeBodyDecay.jl): core numerical implementation,
- [SymbolicThreeBodyDecay.jl](https://github.com/mmikhasenko/SymbolicThreeBodyDecays.jl): a nice small addition.

```
DPD: one that does spin right
```

Link to a [recent analysis](https://lc2pkpi-polarimetry.docs.cern.ch/index.html) [Python + Julia], worth checking out.

"""

# ╔═╡ 002b5c97-f66d-4ebd-b040-b3e341675575
theme(:wong2, frame=:box, grid=false, minorticks=true,
    guidefontvalign=:top, guidefonthalign=:right,
    xlim=(:auto, :auto), ylim=(:auto, :auto),
    lw=1.2, lab="", colorbar=false,
    bottom_margin=3mm, left_margin=3mm)

# ╔═╡ a2c951b4-7af4-48e2-ae41-70d9c3693d48
md"""
## Numerical
"""

# ╔═╡ 4c03daab-c976-4044-9023-d90e0dbb7a21
begin
	const mΛb = 5.61960
	const mπ = 0.14
	# intermediate resonances
	const mΣb = 5.81065;
	const ΓΣb = 0.005;
	const mΣb_x = 5.83032;
	const ΓΣb_x = 0.0094;
	
	const mΛb2S = 6.3; # just a peak of the plot
end

# ╔═╡ 2ab2ebe6-2177-4d73-af38-d27d79c30bf7
md"""
## Dynamics
"""

# ╔═╡ 6e1c0488-4cba-4e34-aeda-7d5c232243b4
full_a(dpp) = full_a(dpp.σs, dpp.two_λs)

# ╔═╡ d9c8a2e6-869e-4dbc-8c0d-8e88d69e1a4f
md"""
## Integration with Symbolic Computation
"""

# ╔═╡ cf4a76bd-a265-47be-a26c-71ebdc541838
begin
    @syms m1 m2 m3 m0
    @syms σ1 σ2 σ3
end;

# ╔═╡ dd425a12-1ffc-4eb3-9a05-64bdbfce4709
ms = ThreeBodyMasses(m1, m2, m3; m0)

# ╔═╡ d7b53398-05a8-4bb6-92d3-bfef3991ed2b
σs = Invariants(ms; σ1, σ2)

# ╔═╡ 2fe649e9-1f89-43ba-90f6-3243e15d85b0
function spinparity(p)
    pt = (p[2]..., p[1])
    jpv = str2jp.(pt)
    getproperty.(jpv, :j) .|> x2 |> ThreeBodyDecay.SpinTuple,
    getproperty.(jpv, :p) |> ThreeBodyDecay.ParityTuple
end

# ╔═╡ 74a39810-1a26-44f8-8525-a3c3836372ad
begin
	ifstate = "1/2+" => ("0-", "1/2+", "0-")
	two_js, Ps = ifstate |> spinparity
end

# ╔═╡ 85ed6757-c4f2-4cb7-8f73-a97aaa262bc7
tbs = ThreeBodySystem(mπ, mΛb, mπ; m0=mΛb2S, two_js)

# ╔═╡ 6232f908-2621-4907-ac76-4e7df350af0d
dpp = randomPoint(tbs)

# ╔═╡ c8b1f1a2-1590-45ed-a70a-3600b8dd9e2e
data = flatDalitzPlotSample(tbs.ms)

# ╔═╡ 1fbf018b-8005-49b7-a51b-87d19ef9fd14
histogram2d(getindex.(data, :σ3), getindex.(data, :σ2), bins=50)

# ╔═╡ 04cc49ff-ee7c-40c2-aa04-0f26a80e96b4
begin
	# lineshape
	dc_Σb1 = DecayChainLS(1, σ -> BW(σ, mΣb, ΓΣb); tbs, Ps, two_s=1, parity='+')
	dc_Σb3 = DecayChainLS(3, σ -> BW(σ, mΣb, ΓΣb); tbs, Ps, two_s=1, parity='+')
	dc_Σb_x1 = DecayChainLS(1, σ -> BW(σ, mΣb_x, ΓΣb_x); tbs, Ps, two_s=3, parity='+')
	dc_Σb_x3 = DecayChainLS(3, σ -> BW(σ, mΣb_x, ΓΣb_x); tbs, Ps, two_s=3, parity='+')
	dc_f0 = DecayChainLS(2, σ -> 1.0; tbs, Ps, two_s=0, parity='+')
end;

# ╔═╡ 0b2ae68c-abd7-4dc9-81d5-ccab8ca3338b
full_a(σs, two_λs) = 
	sum(amplitude(ch, σs, two_λs)
		for ch in [dc_Σb1, dc_Σb_x1, dc_Σb3, dc_Σb_x3])

# ╔═╡ df9db828-b0e2-401f-bdcd-a81720dd3ec1
total_I(σs) = sum(abs2(full_a(σs, two_λs)) for two_λs in itr(two_js))

# ╔═╡ cce113c3-71b4-4b10-a329-ae9c1d73971d
let
    σ3v = range(lims3(tbs.ms)..., length=151)
    σ2v = range(lims2(tbs.ms)..., length=150)
    cali = [
        let σs = Invariants(tbs.ms; σ2, σ3)
            Kibble(σs, tbs.ms^2) < 0 ? total_I(σs) : NaN
		end for (σ2, σ3) in Iterators.product(σ2v, σ3v)
    ]
    heatmap(σ3v, σ2v, cali)
end

# ╔═╡ 1e450d31-9393-4655-8909-96732007e963
cali = total_I.(data)

# ╔═╡ 3bd38333-9efe-4a9d-8b5f-61482e1172d6
histogram2d(getindex.(data, :σ3), getindex.(data, :σ2), weights=cali, bins=50)

# ╔═╡ 4caea52a-2ae8-46a7-acd5-7a15c629ae8f
# Ωb (1/2+) -> Ωc(JP) [-> Lc(1/2+) K(0+)] π(0+) 
# 
function decay_via_xic2645(jpR)
    # 
    ifstate = "1/2+" => ("1/2+", "0-", "0-")
    js, Ps = ifstate |> spinparity
    tbs = ThreeBodySystem(ms, js)
    # 
    # resonance
    Rjp = str2jp(jpR)
    R(σ) = Sym("R")
    # 
    # decay chain
    dc_all = DecayChainsLS(3, R; two_s=Rjp.j |> x2, parity=Rjp.p, Ps, tbs)
	size(dc_all) != (1,1) && @info "More couplings"
    return dc_all |> first
end

# ╔═╡ a3024546-3529-41ce-81da-b43bceeb6bd3
dc_test = let
    dc = decay_via_xic2645("1/2+")
    @show dc.HRk.two_ls
    @show dc.Hij.two_ls
    dc
end;

# ╔═╡ 45e58167-a78e-4852-8af9-c62349583723
function I(jp0)
    dc = decay_via_xic2645(jp0)
    full_amplitude = sum(itr(dc.tbs.two_js)) do two_λs
        A = amplitude(dc, σs |> StickySymTuple, two_λs .|> Sym; refζs=(3, 1, 1, 3)).doit()
		abs2(A)
    end
    full_amplitude |> simplify |> simplify
end

# ╔═╡ ff32e136-eeda-4de5-bb93-7ff000dd760e
md"""
## Test different $J^P$ of the resonance
"""

# ╔═╡ b194ffb0-805d-48b1-99a8-285dabfc698a
I("1/2+") == I("1/2-") && I("1/2+")

# ╔═╡ dcf6a109-971e-4532-b48f-d548188c8dd4
I("3/2+") == I("3/2-") && I("3/2+")

# ╔═╡ 1b3dce79-2526-4c16-90b0-62a326396c46
I("5/2+") == I("5/2-") && I("5/2+")

# ╔═╡ 2942c655-c55e-4ef0-a1a1-c584d724a849
begin
    plot(title="\$\\varOmega_b^-(1/2^\\pm) \\to \\varOmega_c^{**}(J^P) \\pi^-\$")
    for jp in vec(string.(1:2:5) .* "/2+")
		I_jp = I(jp)
		@syms θ12::real=>"θ_12" cosθ
        expr = I_jp.subs(Dict(
			dc_test.Xlineshape(0) => Sym(1),
			θ12 => acos(cosθ)
		), simultaneous=true) + 1e-7 * cosθ
        plot!(expr, -1, 1, lab=jp, lw=1.5)
    end
    plot!(ylim=(0, :auto), legend_title="JP",
		xlab="cosine of \$\\varOmega_c^{**0}\$ helicity angle")
end

# ╔═╡ Cell order:
# ╟─a5d25226-f465-428f-931b-64825963c255
# ╟─b4c00f7e-ce88-4516-a43a-d63e534bad6e
# ╠═f4bbad10-7e14-11ee-0343-7da1e5ed0283
# ╠═002b5c97-f66d-4ebd-b040-b3e341675575
# ╟─a2c951b4-7af4-48e2-ae41-70d9c3693d48
# ╠═4c03daab-c976-4044-9023-d90e0dbb7a21
# ╠═85ed6757-c4f2-4cb7-8f73-a97aaa262bc7
# ╠═6232f908-2621-4907-ac76-4e7df350af0d
# ╠═c8b1f1a2-1590-45ed-a70a-3600b8dd9e2e
# ╠═1fbf018b-8005-49b7-a51b-87d19ef9fd14
# ╟─2ab2ebe6-2177-4d73-af38-d27d79c30bf7
# ╠═74a39810-1a26-44f8-8525-a3c3836372ad
# ╠═04cc49ff-ee7c-40c2-aa04-0f26a80e96b4
# ╠═0b2ae68c-abd7-4dc9-81d5-ccab8ca3338b
# ╠═6e1c0488-4cba-4e34-aeda-7d5c232243b4
# ╠═df9db828-b0e2-401f-bdcd-a81720dd3ec1
# ╠═cce113c3-71b4-4b10-a329-ae9c1d73971d
# ╠═1e450d31-9393-4655-8909-96732007e963
# ╠═3bd38333-9efe-4a9d-8b5f-61482e1172d6
# ╟─d9c8a2e6-869e-4dbc-8c0d-8e88d69e1a4f
# ╠═cf4a76bd-a265-47be-a26c-71ebdc541838
# ╠═dd425a12-1ffc-4eb3-9a05-64bdbfce4709
# ╠═d7b53398-05a8-4bb6-92d3-bfef3991ed2b
# ╠═2fe649e9-1f89-43ba-90f6-3243e15d85b0
# ╠═4caea52a-2ae8-46a7-acd5-7a15c629ae8f
# ╠═a3024546-3529-41ce-81da-b43bceeb6bd3
# ╠═45e58167-a78e-4852-8af9-c62349583723
# ╟─ff32e136-eeda-4de5-bb93-7ff000dd760e
# ╠═b194ffb0-805d-48b1-99a8-285dabfc698a
# ╠═dcf6a109-971e-4532-b48f-d548188c8dd4
# ╠═1b3dce79-2526-4c16-90b0-62a326396c46
# ╠═2942c655-c55e-4ef0-a1a1-c584d724a849