c
c----------------------------------------------------------------------------
c
      subroutine get_direction(costheta, sintheta, cosphi, sinphi, 
     1                         pos_tot)
c
c
c       T. Prokscha,    15-Sep-2000     PSI
c
c       Unix Version
c
c------------------------------------------------------------------------------
c
c       generates the initial direction of the decay-e+
c
      implicit none
c
#define CERNLIB_TYPE
#include "geant321/gclist.inc"
#include "geant321/gctrak.inc"
#include        "common.inc"
#include        "cwn.inc"
c
      real*4 pos_tot                  ! total e+ energy including mass
      real*4 p                        ! polarization == bfield(4)
      real*4 costheta,phi             ! decay direction in rotated system
      real*4 sintheta
      real*4 cosphi, sinphi
      real*4 eps                      ! pos_tot/MAX_TOTAL_ENERGY
      real*4 Dw                       ! asymmetry factor
      real*4 random(3)
      real*4 x,y,z,r
      real*4 t_in_sample,tof,beta,gamma_sqr ! to add relaxation in sample
      real*4 mass, mass_sqr, pz_sqr, z0
      real*4 arg
c
c-----------------------------------------------------------------------------
c   
        if ( lsets(2) .eq. 0 ) then ! positron created in MCP2 or sample
          t_in_sample = tofg*1.e9   ! time in nsec
          arg = beam_parameter(4)*t_in_sample
          if (arg > 10. ) then
           p = 0.
          else
           p = bfield(5) / 100. * dexp(-dble(arg)) ! expo relaxation
          endif
        else
          x = vect(1)
          y = vect(2)
          z = vect(3)
          r = sqrt(x**2 + y**2)
          if (1.4 .le. z .and. z .le. 1.41 .and. r .le. 3.5) then
            mass     = MU_MASS*1000. ! muon mass in MeV/c
            mass_sqr = mass*mass
            pz_sqr   = nt_pz0*nt_pz0
            gamma_sqr = (mass_sqr + pz_sqr)/mass_sqr  
            beta = sqrt(gamma_sqr - 1.)
            if ( lsets(2) .lt. 0 ) z0 = -lsets(2)
            tof         = z0/beta/C_LIGHT
            t_in_sample = (tofg - tof)*1e9 ! time in nsec 
            p = bfield(5) / 100. * exp(-beam_parameter(4)*t_in_sample) ! expo relaxation
          else
            p = bfield(4) / 100.
          endif
        endif
c        print*,' lsets(2)               = ',lsets(2)
c        print*,' nt_pz0                 = ',nt_pz0
c        print*,' beta                   = ',beta
c        print*,' tofg                   = ',tofg
c        print*,' tof                    = ',tof
c        print*,' t_in_sample            = ',t_in_sample
c        print*,' Polarization           = ',p
c
#if defined (OS_UNIX)
        call ranlux(random, 3)	! random generator from Mathlib
#else
        random(1) = ran(ix1)
        random(2) = ran(ix1)
        random(3) = ran(ix1)
#endif
        if (p.eq.0.0) then
c
          costheta = 1. - 2.*random(1)
c 
        else
c
c           see TP Logbook 7, p. 167 for getting the equations
c
          eps = pos_tot / MAX_TOTAL_ENERGY
          Dw = (2. * eps - 1.) / (3. - 2. * eps) * p

          costheta = 1/Dw*(-1. + 
     1               sqrt( 1. - 2*Dw*( 2.*random(2) - 1. ) + Dw**2))
        endif
c
        sintheta = sqrt(1. - costheta**2)
c
        phi=random(3) * TWO_PI
c
        cosphi = cos(phi)
        sinphi = sin(phi)
c
c       now we have the decay direction in the rest system of the muon
c       return.
c
        return
        end