### A Pluto.jl notebook ###
# v0.20.21

using Markdown
using InteractiveUtils

# ╔═╡ 430d4cf8-59dc-11f1-069f-f3a576bc61c8
# ╠═╡ show_logs = false
begin
	using Pkg
	Pkg.activate(mktempdir())
	Pkg.add([
		Pkg.PackageSpec(url="https://github.com/JuliaHEP/FourVectors.jl.git"),
		Pkg.PackageSpec(url="https://github.com/mmikhasenko/RamboOnDiet.jl.git"),
		Pkg.PackageSpec("HadronicLineshapes"),
		Pkg.PackageSpec("DataFrames"),
		Pkg.PackageSpec("Plots"),
	])
	# 
	using FourVectors
	using RemboOnDiet
	using HadronicLineshapes
	using Plots
	using DataFrames
end

# ╔═╡ 565b2bf7-48e4-4783-8376-100e4166361d
md"""
# Example-03: phase space
"""

# ╔═╡ 6d1685cb-b1b3-4921-a69b-043fac7de7f7
theme(:boxed)

# ╔═╡ 92300ccc-0b51-4f32-8323-db3fccc6f54e
default_pars = (
	m1 = 0.8, Γ1 = 0.15, c1 = 2.00*cis(0.3),
	m2 = 1.2, Γ2 = 0.05, c2 = 1.0,
	m3 = 1.6, Γ3 = 0.05, c3 = 3.0)

# ╔═╡ 5a850db8-eceb-4a72-99c1-90b5604b5bd9
amplitude_model(σ1::Real, σ2::Real, σ3::Real; pars) = let
	(; m1, Γ1, c1) = pars
	(; m2, Γ2, c2) = pars
	(; m3, Γ3, c3) = pars
	# 
	c1 * BreitWigner(m1,Γ1)(σ1) +
	c2 * BreitWigner(m2,Γ2)(σ2) + 
	c3 * BreitWigner(m3,Γ3)(σ3)
end

# ╔═╡ 7f26d8fe-2534-4d1a-a6ab-467ff260ca0a
function amplitude(x::Array{<:FourVector}; pars)
	p1, p2, p3 = x
	σ1 = mass2(p2+p3)
	σ2 = mass2(p3+p1)
	σ3 = mass2(p1+p2)
	# 
	amplitude_model(σ1,σ2,σ3; pars)
end

# ╔═╡ 28b0605a-68c6-41a5-a6bc-19ec060b8826
# x ∈ R¹², p ∈ R⁶×C³
unnpormalized_density(x; pars) = abs2(amplitude(x; pars))

# ╔═╡ 375ed5f1-5862-4c4f-8840-cf6dfa316df5
md"""
## What is $\vec x$
`unnpormalized_density` is defined on a `x` which is 3x4 dimentions
```julia
 p1 = [p1x,p1y,p1z,E1],
 p2 = [p2x,p2y,p2z,E2],
 p3 = [p3x,p3y,p3z,E3]
```
"""

# ╔═╡ b978c053-b61c-4b1b-bfdf-ecd566e55b89
md"""
## Evaluating the density on MC sample
"""

# ╔═╡ fa25432c-5dc6-4df2-8606-fddcef1cf9e0
data = let N = 100_000
	masses = [0.93827208816, 0.493677, 0.13957039]
	generator = PhaseSpaceGenerator(masses, 2.28646)
	df = [rand(generator) for _ in 1:N] |> DataFrame
	# 
	select(df, :momenta => ByRow() do (p1, p2, p3)
		σ1, σ2 = mass2(p2+p3), mass2(p3+p1)
		(; p1, p2, p3, σ1, σ2)
		end => AsTable, :weight)
end;

# ╔═╡ 74cedaae-1959-461b-a0e9-e75442837fc3
density = select(data, [:p1,:p2,:p3] => ByRow() do p1, p2, p3
	unnpormalized_density([p1,p2,p3]; pars=default_pars)
end => :weight);

# ╔═╡ c4c65450-a11f-4611-8de9-2910400b009d
histogram2d(data.σ1, data.σ2, bins=50; weights=data.weight)

# ╔═╡ f2f0bc51-30ea-441e-be2b-7bd0e116600d
histogram2d(data.σ1, data.σ2, bins=50;
			weights=data.weight .* density.weight)

# ╔═╡ Cell order:
# ╟─565b2bf7-48e4-4783-8376-100e4166361d
# ╠═430d4cf8-59dc-11f1-069f-f3a576bc61c8
# ╠═6d1685cb-b1b3-4921-a69b-043fac7de7f7
# ╠═7f26d8fe-2534-4d1a-a6ab-467ff260ca0a
# ╠═92300ccc-0b51-4f32-8323-db3fccc6f54e
# ╠═5a850db8-eceb-4a72-99c1-90b5604b5bd9
# ╠═28b0605a-68c6-41a5-a6bc-19ec060b8826
# ╟─375ed5f1-5862-4c4f-8840-cf6dfa316df5
# ╟─b978c053-b61c-4b1b-bfdf-ecd566e55b89
# ╠═fa25432c-5dc6-4df2-8606-fddcef1cf9e0
# ╠═74cedaae-1959-461b-a0e9-e75442837fc3
# ╠═c4c65450-a11f-4611-8de9-2910400b009d
# ╠═f2f0bc51-30ea-441e-be2b-7bd0e116600d