!  GeLetDecay.f90
!
!  Takes care of the decay of a particle according to the decays in
!  decayList. Produced particles are added to the end of treeList.
!
!  Detlef Irmscher, Thomas Ullrich, Uni Heidelberg, September 1993
!  Last update: 17.3.95 tu
!
subroutine GeLetDecay(currentIndex)
    use genesis_m, self=>GeLetDecay
    use f90_kind
    use FormDecay_m
    implicit none
    integer, intent(in)  :: currentIndex
    type(tree), pointer  :: currentTree
    real(kind=double)    :: branchsum, decayrndm
    integer              :: idec, ichild, use_decay, which_decay(30)
    real                 :: rand
    integer              :: temp_index

    ! set up tree
    if (currentIndex < 1 .or. currentIndex > lasttree) then
        write (errmess, *) 'funny index: ', currentIndex
        call GeErrorMessage ("GeLetDecay", errmess)
        return
    endif
    currentTree => treeList(currentIndex)

    ! find out which decay to use
    if (is_mass_weight == 0) then               ! use branching ratios
        ! random number determining the decay
        call random_number(rand)
        decayrndm = rand
        branchsum = 0_double

        ! loop over decays
        use_decay = 0
        do idec = 1, size(decayList)
            ! check if decay for this particle
            if (decayList(idec)%parent /= currentTree%parttype)  cycle

            ! sum branching ratio
            branchsum = branchsum + decayList(idec)%branchingRatio

            ! use this decay?
            if (decayrndm > branchsum .and. is_mass_weight .eq. 0) cycle

            ! decay enabled? if not, exit loop
            if (.not.decayList(idec)%enabled) exit

            use_decay = idec
            exit
        enddo
    else                                          ! each decay equally
        ! get number of enabled decays for this particle
        use_decay = 0
        do idec = 1, size(decayList)
            ! check if decay for this particle
            if (decayList(idec)%parent /= currentTree%parttype)  cycle

            ! decay enabled?
            if (decayList(idec)%enabled) then
                use_decay = use_decay + 1
                which_decay(use_decay) = idec
            endif
        enddo

        ! get number of decay to use
        if (use_decay /= 0) then
            decayrndm = use_decay                 ! to ensure equal weighting
            call random_number(rand)
            idec = int(rand*use_decay) + 1
            use_decay = which_decay(idec)
        endif
    endif

    ! now decay, if decaymode identified
    if (use_decay .ne. 0) then
        idec = use_decay

        ! mark decay
        currentTree%decaytype => decayList(idec)
        currentTree%decayed   = .true.

        ! set up pointer to children
        do ichild = 1, decayList(idec)%nBody
            temp_index =                                  &
               GeAllocTree (currentIndex, decayList(idec)%child(ichild)%name)
            treeList(currentIndex)%child(ichild) = temp_index
        enddo

        ! reset current tree in case of realloc by allocTree
        currentTree => treeList(currentIndex)

        ! do the actual decay
        if (decayList(idec)%nBody.eq.2) then
            call GeDecay2Body (currentIndex)
        else if (decayList(idec)%nBody.eq.3) then
            if (.not.GeIsDalitzDecay(decayList(idec))) then
                call GeDecay3Body (currentIndex)
            else
                call GeDalitzDecay (currentIndex)
            endif
        else
            write (errmess, *) "decay ", decayList(idec)%name,              &
                               "has funny nBody: ", decayList(idec)%nBody
            call GeErrorMessage ("GeLetDecay", errmess)
        endif

        ! set children to stable, if indicated in decay
        if (.not.decayList(idec)%decayChildren) then
            do ichild = 1, decayList(idec)%nBody
                treeList(currentTree%child(ichild))%isStable = .true.
            enddo
        endif

        ! fill in mass weight
        if (is_mass_weight .eq. 1) then
            do ichild = 1, decayList(idec)%nBody
                treeList(currentTree%child(ichild))%mass_weight =           &
                    currentTree%mass_weight*decayList(idec)%branchingRatio  &
                    *decayrndm
            enddo
        endif

        ! fill in pt weight
        if (is_pt_weight .eq. 1) then
            do ichild = 1, decayList(idec)%nBody
                treeList(currentTree%child(ichild))%pt_weight =           &
                         currentTree%pt_weight
            enddo
        endif

    endif

end subroutine GeLetDecay