!  GeDecay3Body.f90
!
!  Takes care of the three body decay assuming isotropic phase space.
!  Parent particle and decay type are specified in old_tree.
!
!  Detlef Irmscher, Thomas Ullrich, Uni Heidelberg, September 1993
!  Last update: 17.3.95 tu
!
subroutine GeDecay3Body (oldindex)
    use f90_kind
    use genesis_m, self=>GeDecay3Body
    use mom_m
    implicit none
    integer, intent(in)         :: oldindex
    type(tree), pointer         :: old_tree, child1, child2, child3
    real(kind=double)           :: mp, m(3), m12, m13, e3, p3, e12
    real(kind=double)           :: m12lo, m12hi, t1
    real(kind=double)           :: costheta, sintheta, phi, cosphi, sinphi
    real(kind=double)           :: m13lo, m13hi, e1s, e3s, m13max, m13min
    integer                     :: i
    type(mom4)                  :: intmom(2)
    type(mom4)                  :: p4boost
    real                        :: rand(2)

    ! set up trees
    old_tree => treeList(oldindex)
    child1   => treeList(old_tree%child(1))
    child2   => treeList(old_tree%child(2))
    child3   => treeList(old_tree%child(3))

    ! check if really 3-body decay
    if (old_tree%decaytype%nBody /= 3) then
	write (errmess, *) 'Incorrect decay type:',         &
	old_tree%decaytype%nBody,                           &
	'-body decay ', old_tree%decaytype%name
        call GeErrorMessage ("GeDecay3Body", errmess, fatal = .true.)
    endif

    ! generate particle 3 in 123 cms:
    ! calculate mass limits
    mp    = old_tree%parttype%mass
    do i = 1, 3
        m(i)  = old_tree%decaytype%child(i)%mass
    enddo
    m12lo = (m(1) + m(2))**2
    m12hi = (mp   - m(3))**2
    m13lo = (m(1) + m(3))**2
    m13hi = (mp   - m(2))**2

    ! calculate particle 3
    ! get evenly distributed in phase space, then check kinematical limits
    do
        call random_number(rand)
        m12   = sqrt(m12lo + (m12hi - m12lo)*rand(1))
        m13   = sqrt(m13lo + (m13hi - m13lo)*rand(2))
        e1s = (m12**2 + m(1)**2 - m(2)**2)/(2.*m12)
        e3s = (mp**2  - m12**2  - m(3)**2)/(2.*m12)

        t1 = e1s**2 - m(1)**2         ! tmp to prevent sqrt domain errors
        if (t1 < 0) cycle

        m13max = sqrt((e1s + e3s)**2 - (sqrt(t1) -       &
                                        sqrt(e3s**2 - m(3)**2))**2)
        m13min = sqrt((e1s + e3s)**2 - (sqrt(t1) +       &
                                        sqrt(e3s**2 - m(3)**2))**2)
        if (m13 <= m13max .and. m13 >= m13min) exit
    enddo

!    m12   = m12lo + (m12hi - m12lo)*rndm()   wrong: not uniform phase space

    e3    = (mp**2 + m(3)**2 - m12**2)/(2.*mp)
    p3    = sqrt((e3 + m(3))*(e3 - m(3)))

    call random_number(rand)
    costheta = 2.*rand(1) - 1.
    sintheta = sqrt((1. + costheta)*(1. - costheta))
    phi      = TWOPI * rand(2)
    sinphi   = sin(phi)
    cosphi   = cos(phi)

    ! fill 4-vector of child 3
    child3%p4%E    = e3
    child3%p4%p%px = p3*sintheta*cosphi
    child3%p4%p%py = p3*sintheta*sinphi
    child3%p4%p%pz = p3*costheta

    ! generate 2 body decay in 12 rest frame
    call GeDeca2 (m12, m(1), m(2), intmom)

    ! boost back into 123 cms
    e12 = sqrt(p3**2 + m12**2)
    p4boost%E   = e12
    p4boost%p   = -child3%p4%p

    child1%p4 = GeBoost(intmom(1), p4boost)
    child2%p4 = GeBoost(intmom(2), p4boost)

    ! boost into lab system
    child1%p4  = GeBoost (child1%p4, old_tree%p4)
    child2%p4  = GeBoost (child2%p4, old_tree%p4)
    child3%p4  = GeBoost (child3%p4, old_tree%p4)

contains

   !
   !   2 body decay with isotropic angular distribution in cms.
   !
   subroutine GeDeca2 (mp, md1, md2, pout)
       implicit none
       real(kind=double), intent(in)   :: mp, md1, md2
       type(mom4), intent(out)         :: pout(2)
       real(kind=double)               :: e1, e2, p1
       real(kind=double)               :: costheta, sintheta, phi

       ! calculate CM energies and momenta of children
       e1  = (mp**2 + md1**2 - md2**2)/(2.*mp)
       e2  = mp - e1
       p1  = sqrt((e1 + md1)*(e1 - md1))

       ! orientation in decaying particle rest frame
       call random_number(rand)
       costheta = 2.*rand(1) - 1.
       sintheta = sqrt((1. + costheta)*(1. - costheta))
       phi      = TWOPI * rand(2)

       ! fill 4-vectors of daughters
       pout(1)%E    = e1
       pout(1)%p%px = p1*sintheta*cos(phi)
       pout(1)%p%py = p1*sintheta*sin(phi)
       pout(1)%p%pz = p1*costheta
       pout(2)%p    = -pout(1)%p
       pout(2)%E    = sqrt(pout(2)%p*pout(2)%p + md2**2)

   end subroutine GeDeca2

end subroutine GeDecay3Body