feat: Add cosθ projection integrand for three-body decay analysis (#71)

* feat: Add cosθ projection integrand for three-body decay analysis

- Implement `project_cosθij_intergand` function to project onto cos(θij)
- Add new function to ThreeBodyDecays module exports
- Include test case to verify cosθ projection integral matches phase space integral

* formatting and cspell

Author: mmikhasenko

.cspell/project.txt

@@ -18,6 +18,7 @@ guidefontvalign
 helicities
 helicity
 imag
+intergand
 issequential
 jldoctest
 JPAC

src/ThreeBodyDecays.jl

@@ -74,7 +74,8 @@ include("decay_model.jl")
 export change_basis_3from1, change_basis_1from2, change_basis_2from3
 include("cross_channel.jl")
 
-export phase_space_integrand, projection_integrand
+export phase_space_integrand
+export projection_integrand, project_cosθij_intergand
 include("integrand.jl")
 
 end  # module ThreeBodyDecays

src/integrand.jl

@@ -1,5 +1,5 @@
 """
-	phase_space_integrand(function_σs, ms; k::Int)
+    phase_space_integrand(function_σs, ms; k::Int)
 
 Calculate the phase space integrand for a given function `function_σs`,
 and mass tuple `ms`. The key argument `k` specifies the mapping: `σk->[0,1]`, zk->[0,1].
@@ -38,7 +38,7 @@ function phase_space_integrand(function_σs, ms; k::Int)
 end
 
 """
-	projection_integrand(function_σs, ms, σk; k)
+projection_integrand(function_σs, ms, σk; k)
 
 Calculate the projection integrand for a given function `function_σs`,
 mass tuple `ms`, and Mandelstam variable `σk`, with `k` specified by a keyword argument.
@@ -53,9 +53,9 @@ It returns an integrand function of x, x ∈ [0,1] to pass to a numerical integr
 # Usage
 ```julia
 plot(4.2, 4.6) do e1
-	I = Base.Fix1(unpolarized_intensity, model)
-	integrand = projection_integrand(I, masses(model), e1^2; k = 3)
-	e1 * quadgk(integrand, 0, 1)[1]
+I = Base.Fix1(unpolarized_intensity, model)
+integrand = projection_integrand(I, masses(model), e1^2; k = 3)
+e1 * quadgk(integrand, 0, 1)[1]
 end
 ```
 """
@@ -72,3 +72,42 @@ function projection_integrand(function_σs, ms, σk; k)
     end
     return integrand
 end
+
+"""
+    project_cosθij_intergand(fs, ms, z; k)
+
+Calculate the integrand for projecting onto cos(θij) for a given function `fs`,
+mass tuple `ms`, and cosine value `z`, with `k` specifying the spectator.
+Returns an integrand function of x ∈ [0,1] to pass to a numerical integrator.
+
+# Arguments
+- `fs`: A function that takes a MandelstamTuple and returns a scalar.
+- `ms`: A MassTuple containing the masses.
+- `z`: The cosine value to project onto.
+- `k`: The spectator index.
+
+# Returns
+- An integrand function for numerical integration.
+
+# Example
+```julia
+let Nb = 100
+    [sum(range(-1, 1, Nb)) do z
+        integrand = project_cosθij_intergand(I, ms, z; k)
+        quadgk(integrand, 0, 1)[1] * 2/Nb
+    end for k in 1:3]
+end
+```
+"""
+function project_cosθij_intergand(fs, ms, z; k)
+    i, j = ij_from_k(k)
+    mi, mj, mk, m0 = ms[i], ms[j], ms[k], ms[4]
+
+    function integrand(xσk)
+        σs = x2σs([xσk, (z + 1) / 2], ms; k)
+        σk = σs[k]
+        jac = sqrt(Kallen(σk, mk^2, m0^2) * Kallen(σk, mi^2, mj^2)) / σk
+        return fs(σs) / 2 * jac * ((m0 - mk)^2 - (mi + mj)^2)
+    end
+    return integrand
+end

test/test_integrals.jl

@@ -1,7 +1,7 @@
 using ThreeBodyDecays
+using QuadGK
 using Test
 using Cuba
-using QuadGK
 
 @testset "phase space mapping" begin
     ms = ThreeBodyMasses(0.141, 0.142, 0.143; m0 = 3.09)
@@ -36,3 +36,31 @@ end
     @test quadgk(projection_integrand(σs -> 1, ms, 0.7^2; k = 3), 0, 1)[1] ≈
           8.26409477583501
 end
+
+@testset "cosθ projection integral" begin
+    ms = ThreeBodyMasses(0.141, 0.142, 0.143; m0 = 3.09)
+    I = σs -> σs[1] * σs[2] - σs[3]  # Unit amplitude for testing
+
+    # Test integration over cosθ matches phase space integral
+    Nb = 500
+    # Integration over cosθ
+    cosθ_integral = [
+        sum(range(-1, 1, Nb)) do z
+            integrand = project_cosθij_intergand(I, ms, z; k)
+            quadgk(integrand, 0, 1)[1] * 2 / Nb
+        end for k = 1:3
+    ]
+
+    # Integration over σk
+    σk_integral = [
+        sum(range(lims(ms; k)..., Nb)) do σk
+            integrand = projection_integrand(I, ms, σk; k)
+            quadgk(integrand, 0, 1)[1] * diff(collect(lims(ms; k)))[1] / Nb
+        end for k = 1:3
+    ]
+
+    # Both methods should give approximately the same result
+    @test isapprox(cosθ_integral[1], σk_integral[1], rtol = 1e-2)
+    @test isapprox(cosθ_integral[2], σk_integral[2], rtol = 1e-2)
+    @test isapprox(cosθ_integral[3], σk_integral[3], rtol = 1e-2)
+end