mmikhasenko / ThreeBodyDecay.jl

Julia Implementation of the Dalitz plot decomposition
MIT License
9 stars 0 forks source link

Border fails when s1=s2=s3 is outside dalitz #22

Closed mmikhasenko closed 9 months ago

mmikhasenko commented 1 year ago

MWE

ms = ThreeBodyMasses(1,1,11; m0=14)
border(ms)

Fix

New border function

function fixedborder(ms::ThreeBodyDecay.MassTuple; Nx::Int=300)
    thresholds = ((ms.m2 + ms.m3), (ms.m3 + ms.m1), (ms.m1 + ms.m2))
    T0 = sum(abs2, ms) .- sum(abs2, thresholds)

    # 
    σs0 = alwaysin(ms)[end]
    Δσs0 = Tuple(σs0) .- thresholds .^ 2
    # Δσs = Tuple(σs) .- thresholds .^ 2
    # 
    Ps(θ) = 
        Polynomial([Δσs0[1], -cos(θ)]),
        Polynomial([Δσs0[2], cos(θ + π / 3)]),
        Polynomial([Δσs0[3], cos(θ - π / 3)])
    σs_P(θ) = thresholds .^ 2 .+ Ps(θ)
    # 
    ϕ = Base.Fix2(Kibble, ms^2)
    rborder(θ) = minimum(filter(x -> x > 0, roots(ϕ(σs_P(θ)))))
    σsborder(θ) = # evaluate the polynomials
        map(σs_P(θ)) do P
            P(rborder(θ))
        end
    θs = range(-π, π, length=Nx)
    return ThreeBodyDecay.MandestamTuple.(σsborder.(θs))
end

with the new center

function alwaysin(ms)
    @unpack m0, m1, m2, m3 = ms
    # 
    σ1 = ((m2+m3+m0-m1)/2)^2
    σ3 = σ3of1(0.1,σ1,ms^2)
    σs2 = Invariants(ms; σ1, σ3)
    # 
    σ2 = ((m3+m1+m0-m2)/2)^2
    σ1 = σ1of2(0,σ2,ms^2)
    σs3 = Invariants(ms; σ2,σ1)
    # 
    σ3 = ((m1+m2+m0-m3)/2)^2
    σ2 = σ2of3(0,σ3,ms^2)
    σs1 = Invariants(ms; σ3,σ2)

    σs = ThreeBodyDecay.MandestamTuple(sum((σs1, σs2, σs3) .|> collect) / 3)
    [σs1, σs2, σs3, σs]
end
mmikhasenko commented 9 months ago

Implemented in #32