Add dalitzprojection recipe for 1D Dalitz plot projections (#83)

* Initial plan
* Add dalitzprojection recipe and visualization tutorial
Co-authored-by: mmikhasenko <[email protected]>
* update visualization notebook
* referenced visualization
* update quad to the last argument
* test integrals
* added ;

---------

Co-authored-by: copilot-swe-agent[bot] <[email protected]>
Co-authored-by: mmikhasenko <[email protected]>
Co-authored-by: Mikhail Mikhasenko <[email protected]>

Author: Copilot

.cspell/project.txt

@@ -10,8 +10,10 @@ cosζj
 cosζk
 Dalitz
 dalitzplot
+dalitzprojection
 doublearg
 eachslice
+fillalpha
 findmin
 getindex
 getproperty
@@ -59,6 +61,7 @@ vcat
 Weisskopf
 wignerd
 xbins
+xguide
 xlab
 xlim
 Xlineshape

docs/src/10-energy-dependence.md

@@ -17,7 +17,7 @@ When building the model, these dependencies are specified as follows:
 ```julia
 tbs = ThreeBodySystem(
     ThreeBodyMasses(1.0, 2.0, 3.0; m0=16.0),
-    ThreeBodySpins(0, 0, 0; two_h0=0))
+    ThreeBodySpins(0, 0, 0; two_h0=0));
 
 L = 1  # orbital angular momentum for production
 l = 1  # orbital angular momentum for decay
@@ -32,7 +32,7 @@ chain1 = DecayChain(;
     Hij = VertexFunction(  # decay form-factor
         RecouplingLS((2l, 0)),  # (two_l, two_s)
         BlattWeisskopf{l}(1.5)),  # {l}
-    tbs)
+    tbs);
 ```
 
 Many common models for the propagators and form factors are available in `HadronicLineshapes.jl`.
@@ -43,17 +43,17 @@ Alternatively, you can use lambda functions for simple cases, define any custom
 struct MyCoupling <: Recoupling
     # custom fields
 end
-amplitude(cs::MyCoupling, (two_λa, two_λb), (two_j, two_ja, two_jb))  # to defined
+amplitude(cs::MyCoupling, (two_λa, two_λb), (two_j, two_ja, two_jb));  # to defined
 
 # Custom form factor
 struct MyFormFactor
     # custom fields
 end
-(ff::MyFormFactor)(m0², m1², m2²)  # to defined
+(ff::MyFormFactor)(m0², m1², m2²);  # to defined
 
 # Custom propagator
 struct MyPropagator
     # custom fields
 end
-(prop::MyPropagator)(σk)  # to defined
+(prop::MyPropagator)(σk);  # to defined
 ```

docs/src/30-tutorial.jl

@@ -26,7 +26,7 @@ theme(
 )
 
 # decay Λb ⟶ Jψ p K
-constants = Dict("mJψ" => 3.09, "mp" => 0.938, "mK" => 0.49367, "mLb" => 5.62) # masses of the particles
+constants = Dict("mJψ" => 3.09, "mp" => 0.938, "mK" => 0.49367, "mLb" => 5.62);  # masses of the particles
 
 # `ThreeBodySystem` creates an immutable structure that describes the setup.
 # Two work with particles with non-integer spin, the doubled quantum numbers are stored.
@@ -36,13 +36,13 @@ ms = ThreeBodyMasses(   # masses m1,m2,m3,m0
     constants["mp"],
     constants["mK"];
     m0 = constants["mLb"],
-)
+);
 
 # create two-body system
-tbs = ThreeBodySystem(; ms, two_js = ThreeBodySpins(2, 1, 0; two_h0 = 1)) # twice spin
+tbs = ThreeBodySystem(; ms, two_js = ThreeBodySpins(2, 1, 0; two_h0 = 1));  # twice spin
 
-Conserving = ThreeBodyParities('-', '+', '-'; P0 = '+')
-Violating = ThreeBodyParities('-', '+', '-'; P0 = '-')
+Conserving = ThreeBodyParities('-', '+', '-'; P0 = '+');
+Violating = ThreeBodyParities('-', '+', '-'; P0 = '-');
 
 # - the invariant variables, `σs = [σ₁,σ₂,σ₃]`,
 # - helicities `two_λs = [λ₁,λ₂,λ₃,λ₀]`
@@ -60,7 +60,7 @@ struct BW
     m::Float64
     Γ::Float64
 end
-(bw::BW)(σ::Number) = 1 / (bw.m^2 - σ - 1im * bw.m * bw.Γ)
+(bw::BW)(σ::Number) = 1 / (bw.m^2 - σ - 1im * bw.m * bw.Γ);
 
 
 # chains-1, i.e. (2+3): Λs with the lowest ls, LS
@@ -70,21 +70,21 @@ end
     jp = jp"3/2+",
     Ps = Conserving,
     tbs = tbs,
-)
+);
 Λ1690 = DecayChainLS(
     k = 1,
     Xlineshape = BW(1.685, 0.050),
     jp = jp"1/2+",
     Ps = Conserving,
     tbs = tbs,
-)
+);
 Λ1810 = DecayChainLS(
     k = 1,
     Xlineshape = BW(1.80, 0.090),
     jp = jp"5/2+",
     Ps = Conserving,
     tbs = tbs,
-)
+);
 Λs = (Λ1520, Λ1690, Λ1810);
 
 # chains-3, i.e. (1+2): Pentaquarks with the lowest ls, LS
@@ -94,21 +94,21 @@ Pc4312 = DecayChainLS(;
     jp = jp"1/2+",
     Ps = Conserving,
     tbs = tbs,
-)
+);
 Pc4440 = DecayChainLS(;
     k = 3,
     Xlineshape = BW(4.440, 0.010),
     jp = jp"1/2+",
     Ps = Conserving,
     tbs = tbs,
-)
+);
 Pc4457 = DecayChainLS(;
     k = 3,
     Xlineshape = BW(4.457, 0.020),
     jp = jp"3/2+",
     Ps = Conserving,
     tbs = tbs,
-)
+);
 Pcs = (Pc4312, Pc4440, Pc4457);
 
 # ## Unpolarized intensity
@@ -121,17 +121,17 @@ const model = ThreeBodyDecay(
         [2, 2.1, 1.4im, 0.4, 0.3im, -0.8im],
         [Λ1520, Λ1690, Λ1810, Pc4312, Pc4440, Pc4457],
     ),
-)
+);
 
 # just a random point of the Dalitz Plot
-σs = randomPoint(tbs.ms)
-two_λs = randomPoint(tbs.two_js)
+σs = randomPoint(tbs.ms);
+two_λs = randomPoint(tbs.two_js);
 
 # Model is build, one can compute unpolarized intensity with it
-@show amplitude(model, σs, two_λs)
+@show amplitude(model, σs, two_λs);
 
 # gives a real number - probability
-@show unpolarized_intensity(model, σs)
+@show unpolarized_intensity(model, σs);
 
 # ## Plotting API
 #
@@ -151,20 +151,12 @@ plot(
 # Visualize the matrix element as a function of kinematic variables:
 
 plot(masses(model), σs -> abs2(amplitude(Pc4312, σs, (2, -1, 0, 1))))
-dalitzplot(
-    masses(model),
-    Base.Fix1(unpolarized_intensity, model);
-    iσx = 1,
-    iσy = 3,
-    xlim = (2.0, 4.0),
-    ylim = (20.0, 27.0),
-    xbins = 100,
-    ybins = 120,
-)
 
 # Compute Dalitz plot projections numerically:
 
 plot(4.2, 4.6) do e1
     I = Base.Fix1(unpolarized_intensity, model)
     e1 * quadgk(projection_integrand(I, masses(model), e1^2; k = 3), 0, 1)[1]
 end
+
+# See more in the [Visualization Tutorial](32-visualization.md).

docs/src/32-visualization.jl

@@ -0,0 +1,180 @@
+#md # # Visualization Tutorial
+#nb # # ThreeBodyDecays.jl Visualization Tutorial
+#jl # # ThreeBodyDecays.jl Visualization Tutorial
+
+# This tutorial demonstrates how to visualize three-body decay models using
+# Dalitz plots and Dalitz plot projections. We'll use the same Λb ⟶ J/ψ p K
+# decay example from the main tutorial.
+
+using ThreeBodyDecays
+using Plots
+using QuadGK
+#
+theme(
+    :wong,
+    frame = :box,
+    lab = "",
+    minorticks = true,
+    guidefontvalign = :top,
+    guidefonthalign = :right,
+    xlim = (:auto, :auto),
+    ylim = (:auto, :auto),
+    grid = false,
+);
+
+# ## Setting up Kinematics
+
+# First, let's define the masses and spin of particles:
+
+ms = ThreeBodyMasses(3.09, 0.938, 0.49367; m0 = 5.62);
+tbs = ThreeBodySystem(; ms, two_js = ThreeBodySpins(2, 1, 0; two_h0 = 1));
+
+# ## Phase Space Boundaries
+
+# You can also visualize the phase space boundaries:
+
+plot(
+    border31(ms),
+    xlab = "σ₃ [GeV²]",
+    ylab = "σ₁ [GeV²]",
+    title = "Phase Space Boundary",
+    aspect_ratio = 1,
+);
+
+# The `dalitzplot` and `dalitzprojection` recipes provide a flexible way to
+# visualize three-body decay models. The key parameters are:
+#
+# **For `dalitzplot`:**
+# - `iσx`, `iσy`: Choose which invariants to plot (1, 2, or 3)
+# - `xbins`, `ybins`: Resolution of the plot
+# - `xlims`, `ylims`: Custom axis limits (use `:auto` for automatic)
+#
+# **For `dalitzprojection`:**
+# - First arguments: Mass tuple and intensity function
+# - Last argument: Integration function (e.g., `quadgk` from QuadGK.jl)
+# - `k`: Which invariant to project onto (1, 2, or 3)
+# - `bins`: Number of bins in the projection
+# - `xlims`: Custom axis limits (use `:auto` for automatic)
+
+
+
+
+# ## Setting up the model
+
+# For a simple model, we will use a Breit-Wigner lineshape:
+struct BW
+    m::Float64
+    Γ::Float64
+end
+(bw::BW)(σ::Number) = 1 / (bw.m^2 - σ - 1im * bw.m * bw.Γ);
+
+# Create decay chains for Lambda resonances:
+
+# parities, needed only to identify available couplings
+Ps = ThreeBodyParities('-', '+', '-'; P0 = '+');
+
+const model = ThreeBodyDecay(
+    ["Λ1520", "Λ1690"] .=> zip(
+        [1.0, 1.5],
+        [
+            DecayChainLS(; k = 1, Xlineshape = BW(1.5195, 0.0156), jp = jp"3/2+", Ps, tbs),
+            DecayChainLS(; k = 1, Xlineshape = BW(1.685, 0.050), jp = jp"1/2+", Ps, tbs),
+        ],
+    ),
+);
+
+# ## Dalitz Plot Visualization
+
+# The Dalitz plot shows the intensity distribution across the phase space.
+# We can create a Dalitz plot using the `dalitzplot` recipe:
+
+dalitzplot(
+    masses(model),
+    Base.Fix1(unpolarized_intensity, model);
+    iσx = 1,
+    iσy = 3,
+    xbins = 100,
+    ybins = 100,
+    xlab = "σ₁ = m²(J/ψ p) [GeV²]",
+    ylab = "σ₃ = m²(p K) [GeV²]",
+    title = "Dalitz Plot: Λb → J/ψ p K",
+)
+
+# You can also specify custom limits:
+
+dalitzplot(
+    masses(model),
+    Base.Fix1(unpolarized_intensity, model);
+    iσx = 2,
+    iσy = 1,
+    ylims = (2.0, 3.5),
+    xbins = 120,
+    ybins = 120,
+    xlab = "σ₂ = m²(J/ψ K) [GeV²]",
+    ylab = "σ₁ = m²(J/ψ p) [GeV²]",
+)
+
+# ## Dalitz Plot Projections
+
+# Projections are 1D plots obtained by integrating the Dalitz plot intensity
+# over one invariant mass coordinate. This is useful for visualizing features
+# that may be clearer in 1D.
+
+# Project onto σ₁ (m²(J/ψ p)):
+
+p1 = dalitzprojection(
+    masses(model),
+    Base.Fix1(unpolarized_intensity, model),
+    quadgk;
+    k = 1,
+    bins = 300,
+    xlims = (2.0, 3.5),
+    fill = true,
+    fillalpha = 0.5,
+    xlab = "σ₁ = m²(J/ψ p) [GeV²]",
+    ylab = "Integrated Intensity",
+    title = "Projection onto σ₁ (zoomed)",
+)
+
+# Project onto σ₂ (m²(J/ψ K)):
+
+p2 = dalitzprojection(
+    masses(model),
+    Base.Fix1(unpolarized_intensity, model),
+    quadgk;
+    k = 2,
+    bins = 100,
+    fill = true,
+    fillalpha = 0.5,
+    xlab = "σ₂ = m²(J/ψ K) [GeV²]",
+    ylab = "Integrated Intensity",
+    title = "Projection onto σ₂",
+)
+
+# Project onto σ₃ (m²(p K)) with custom limits:
+
+p3 = dalitzprojection(
+    masses(model),
+    Base.Fix1(unpolarized_intensity, model),
+    quadgk;
+    k = 3,
+    fill = true,
+    fillalpha = 0.5,
+    bins = 150,
+    xlab = "σ₃ = m²(p K) [GeV²]",
+    ylab = "Integrated Intensity",
+    title = "Projection onto σ₃",
+)
+
+# ## Comparing Multiple Projections
+
+# You can compare projections from different models or resonances:
+
+plot(
+    p1,
+    p2,
+    p3,
+    layout = (1, 3),
+    size = (1200, 400),
+    bottom_margin = 6Plots.PlotMeasures.mm,
+)

src/ThreeBodyDecays.jl

@@ -42,6 +42,8 @@ include("phase_space.jl")
 export border31, border12, border23, border13, border21, border32, border
 export DalitzPlot # just a docstring
 export dalitzplot # user plot function
+export DalitzProjection # just a docstring
+export dalitzprojection # user plot function
 include("dalitz.jl")
 
 # Wigner Rotations - zeta angles