Skip to content

Projection recipe #42

@mmikhasenko

Description

@mmikhasenko

It would be nice to introduce the projection recipe.

  1. The recipe should fill natular and have the similar kw as dalitzplot recipe (bins for the number of bins, xlims to the ranges, k to control to which exes it's projected.)

One problem is that the function needs to know how to integrate, while the package does not depend on any integrator
The integrator method should also be provided, (quadgk would be the standard thing to pass).

  1. Here is a prototype to
Prototype on how to use recipe
# ╔═╡ d9037098-a76c-11f0-30fc-531a7805bbd4
begin
	using Plots
	using RecipesBase
end

# ╔═╡ 710536c7-d95c-47db-8226-aae5230606f2
using QuadGK

# ╔═╡ cd0f5776-2901-4704-b0f3-10256a1850d7
theme(:boxed)

# ╔═╡ e3d92a28-7cea-4fba-9a6f-06f93aa4393a
@userplot DalitzPlot

# ╔═╡ 6c17d220-307e-425d-8287-e4840fdda881
@recipe function f(d::DalitzPlot, args...)
	(d.args...,)
end

# ╔═╡ 3b05df27-ad78-4dbb-8b07-450ec2ffda5c
dalitzplot(s0, r->1/(1+r^2); xlim=(0,4), ylim=(0,4), aspect_ratio=1)

# ╔═╡ cd97427d-f0bc-40aa-ac45-4408454a31ee
struct Settings
	r::Float64
	xy::Tuple{Float64,Float64}
end

# ╔═╡ ac092645-7af1-47cd-947a-3d7e7139abff
@recipe function f(s::Settings, func, args...)
	x0,y0 = s.xy
	r0 = s.r
	xv = x0 .+ range(-r0, r0, 100)
	yv = y0 .+ range(-r0, r0, 100)
	Z = map(Iterators.product(xv,yv)) do (x,y)
		r′ = sqrt((x-x0)^2+(y-y0)^2)
		r′ < r0 ? func(r′) : NaN
	end
	seriestype --> :heatmap
	(xv, yv, Z)
end

# ╔═╡ 172049ae-dc71-48fd-9561-f2dea46c7bfa
s0 = Settings(1.1, (2.0, 1.3))

# ╔═╡ da17f7c5-5eac-4eab-b4b8-8441477c80e1
plot(s0, r->1/(1+r^2); xlim=(0,4), ylim=(0,4), aspect_ratio=1)

# ╔═╡ 7e5677de-bf86-4d9c-99e6-f67c0f850f1f
plot(r->r, s0; xlim=(0,4), ylim=(0,4), aspect_ratio=1)

# ╔═╡ b62319c4-0615-4eda-87da-526b8ac2dc00
begin
	struct ProjectionIntegrand{F,X}
		x::X
		s::Settings
		func::F
		k::Int
	end
	ProjectionIntegrand(x, s::Settings, func; k) = ProjectionIntegrand(x, s::Settings, func, k)
	function (pri::ProjectionIntegrand)(var)
		Δx = pri.x-pri.s.xy[1]
		abs(Δx) >= pri.s.r && return 0.0
		ylim1, ylim2 = sqrt(pri.s.r^2-Δx^2) .* [-1, 1]
		Δy = ylim1 + var * (ylim2-ylim1)
		r = sqrt(Δx^2+Δy^2)
		pri.func(r)*(ylim2-ylim1)
	end
end

# ╔═╡ 67ff49c5-fc3e-4c84-b435-793d5cdefed3
begin
	dalits_projection(quadgk, s0, r->r^2;
					  fill=0, fillalpha = 0.3, xlim=(-1.0, 5.0), k=1)
	# plot!(xlim=(-1,3))
end

# ╔═╡ fc9a9276-ee59-4522-92cb-0e2686c14192
@userplot Dalits_Projection

# ╔═╡ b7bdf152-f8aa-4da4-bd47-7e58d227253f
@recipe function f(dp::Dalits_Projection, args...; k::Int, xlims=(:auto, :auto), xpoints=100)
	tool = dp.args[1]
	s = dp.args[2]
	#
	x1 = xlims[1] == :auto ? s.xy[1]-s.r : xlims[1]
	x2 = xlims[2] == :auto ? s.xy[1]+s.r : xlims[2]
	xv = range(x1, x2, xpoints)
	# 
	plotting(x) = tool(ProjectionIntegrand(x, dp.args[2:end]...; k), 0, 1)[1]
	(plotting, xv)
end
  1. Here is an old prototype that handles k, and calls projection_integrand correctly.
begin
	@shorthands dalitz_plot_projection
	# 
	@recipe function f(::Type{Val{:dalitz_plot_projection}}, x, y, z)
		d = plotattributes
		# 
		function_on_dalitz = get(d, :function_on_dalitz, nothing)
		if function_on_dalitz == nothing
			error("what do you want me to plot? Provide it with `f` argument.")
		end
		ms = get(d, :masses, nothing)
		if ms == nothing
			error("Provide masses with the `masses` argument.")
		end
		axis = get(d, :axis_k, nothing)
		if axis == nothing
			error("Provide axis (k) with the `axis` argument.")
		end
		
		function projection_function(σk)
			integrand = projection_integrand(function_on_dalitz, ms, σk; k=axis)
			quadgk(integrand, 0, 1)[1]
		end
		# 
	    seriestype := :path
		x := y
		y := projection_function.(y)
		# 
		()
	end
end

Metadata

Metadata

Labels

No labels
No labels

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions