subroutine randdir(cosx,dir,thmin,thmax)
      real*4 cosx(3),pi,dir(3),thmin,thmax
      real*4 locphi,locxaxis(3),loctheta,loccosx(3)
      COMMON/GCONST/PI,TWOPI,PIBY2,DEGRAD,RADDEG,CLIGHT,BIG,EMASS

      loctheta = acos(cos(thmin)-(cos(thmin)-cos(thmax))*ran(dum))
      locphi = 2.*pi*ran(dum)
      call dezero(dir(3),1.e-10)
      call vecfill(locxaxis,-dir(3)/sqrt(dir(3)**2+dir(1)**2),0.,
     .                  dir(1)/sqrt(dir(3)**2+dir(1)**2) )
      call vecfill(loccosx,sin(loctheta)*cos(locphi),
     .                     sin(loctheta)*sin(locphi),cos(loctheta))
      call trans(loccosx,cosx,dir,locxaxis)
      end
                    
 
      subroutine vecfill(vec,val1,val2,val3)
      real*4 vec(3),val1,val2,val3
      vec(1) = val1
      vec(2) = val2
      vec(3) = val3
      end

      SUBROUTINE TRANS(in,out,zaxis,xaxis)
C       SUBROUTINE TO PERFORM A ROTATION OF THREE MUTUALLY                    
C       ORTHOGONAL AXES.  In is the vector in frame A, out is
C       the vector in frame B (desired result) and zaxis and
C       xaxis are the coordinates of those A axes in the B frame.
C       T is a rotation matrix.
      real*4 in(3),out(3),zaxis(3),xaxis(3)
      real*4 t(3,3)
      integer i,j,k                 
      do 900 i=1,3
         out(i) = 0.
         t(1,i) = xaxis(i)
         t(3,i) = zaxis(i)        
         j = mod(i,3)+1
         k = mod(j,3)+1
         t(2,i) = xaxis(k)*zaxis(j) - xaxis(j)*zaxis(k)
900   continue
      do 950 i=1,3
         do 930 j=1,3
            out(i) = out(i) + t(j,i)*in(j)
930      continue
950   continue
      return
      end

      subroutine dezero(x,y)
      real*4 x,y
      if (x .lt. y .and. x .ge. 0.) x = y
      if (x .gt. -y .and. x .lt. 0.) x = -y
      end