c
#include "geant321/pilot.h"
c---------------------------------------------------------------------------
c
      subroutine gudcay
c
c     in case of muon started with partcode = 500, throw michel spectrum
c     and anisotropic angular distribution
c
c       T. Prokscha,    15-Sep-2000     PSI
c
c       Unix version
c
c       TP,     05-Apr-2001     PSI     Unix, NT, Linux
c
c---------------------------------------------------------------------------
c
        implicit none
c
#define CERNLIB_TYPE
#include  "geant321/gctrak.inc"
#include  "geant321/gckine.inc"
#include  "geant321/gcking.inc"
c
#include  "common.inc"
c
c     Michel spectrum
c
      real*4      pos_tot     ! total positron energy (including rest mass)
      real*4      p0          ! total positron momentum (MeV)
      real*4      px0, py0, pz0 ! components of p0
c
c     decay direction in muon rest system (Spin parallel to x-axis !!!)
c
      real*4      costheta, sintheta, cosphi, sinphi
c
c       4-momentum of muon particle (if decay in flight); if decay in flight
c       do a Lorentz-transformation to get momenta in real space
c
      real*4      decay_rest(4), beta(4), decay_lab(4), gamma, beta_tot
      real*4      decay_rot(3)  ! copy of momentum comp. of decay_lab(4) 
c
c----------------------------------------------------------------------------
c
c       get the Michel spectrum
c
      call michel(pos_tot, p0)
c
      call get_direction(costheta, sintheta, cosphi, sinphi, pos_tot)
c
c     get the single momentum components in the laboratory system
c     the muon spin is directed along the x-axis which points to the
c     LEFT detector.)
c
      px0 = p0 * costheta
      py0 = p0 * sintheta * cosphi
      pz0 = p0 * sintheta * sinphi
c
      decay_rest(1) = px0/1000.    ! GEANT uses GeV/c instead of MeV/c
      decay_rest(2) = py0/1000.
      decay_rest(3) = pz0/1000.
      decay_rest(4) = pos_tot/1000.
c
c     for the Lorentz transformation, maybe not important, because the
c     mu+ are slow
c
      gamma    = getot / MU_MASS
      beta_tot = vect(7) / getot    ! this is  p/E, vect(7) = p
      beta(1)  = -vect(4) * beta_tot ! vect(4) = px/p
      beta(2)  = -vect(5) * beta_tot ! vect(5) = py/p
      beta(3)  = -vect(6) * beta_tot ! vect(6) = pz/p
c
c     negative sign because we need the velocity of
c     the lab-system with respect to the rest-system
c
      beta(4)  = gamma
c
      call gloren(beta, decay_rest, decay_lab)
c        
c     if there is a B-field present we have to rotate the momentum vector
c     at the moment only Bz implemented; exit program if Bx or By are
c     given.
c       
c
      if ( bfield(1) .ne. 0. .or. bfield(2) .ne. 0. or. bfield(3) .ne. 
     1       0.) then
            decay_rot(1) = decay_lab(1)
            decay_rot(2) = decay_lab(2)
            decay_rot(3) = decay_lab(3)
c
            call rotate(decay_rot, tofg)
c
            decay_lab(1) = decay_rot(1)
            decay_lab(2) = decay_rot(2)
            decay_lab(3) = decay_rot(3)
      endif
c
      ngkine = 1
      tofd(1) = 0.
      gkin(1,1) = decay_lab(1)
      gkin(2,1) = decay_lab(2)
      gkin(3,1) = decay_lab(3)
      gkin(4,1) = decay_lab(4)
      gkin(5,1) = 2
      gpos(1,1) = vect(1)
      gpos(2,1) = vect(2)
      gpos(3,1) = vect(3)
c
      return
      end