subroutine guhadr
      implicit none
c
c*keep,gctrak.
*-- author :
      common/gctrak/vect(7),getot,gekin,vout(7),nmec,lmec(30),namec(30)
     + ,nstep ,maxnst,destep,destel,safety,sleng ,step  ,snext ,sfield
     + ,tofg  ,gekrat,upwght,ignext,inwvol,istop ,idecad,iekbin
     + , ilosl, imull,ingoto,nldown,nlevin,nlvsav,istory
      real*8 pout,px,py,pz
      real*8 xhatf,yhatf,zhatf,xhatr,yhatr,zhatr
      real*8 xhat1f,yhat1f,zhat1f,xhat1r,yhat1r,zhat1r
      real*8 psif,tthetaf,thetasc,phisc,psisc
      real*8 theta_sc,phi_sc,psi_sc,x_sc,y_sc,z_sc,r_sc
      real*8 phi_az,theta_az,pi,anorm
      real*8 a,b,c,d,ac,p,r,pp
      real*8 term1,term2,term3,term4
      real*8 prob,pmax,ratio,pmax1,pmax2
c
      real  vect,getot,gekin,vout,destep,destel,safety,sleng ,step
     + ,snext,sfield,tofg  ,gekrat,upwght,iii,ii
c end gctrak
      real rndm(3)
      integer nmec,lmec,namec,nstep ,maxnst,ignext,inwvol,istop
     + ,idecad,iekbin,ilosl, imull,ingoto,nldown,nlevin,nlvsav,istory

c      include 'geant_local.h'

      pi=3.14159265

c      write(*,*)'Calling hadronic interaction package ...'
c      pout=vect(7)*1000.
c      px=pout*vect(4)
c      py=pout*vect(5)
c      pz=pout*vect(6)
c      xhatf=px/pout
c      yhatf=py/pout
c      zhatf=pz/pout
c      anorm=sqrt(xhatf**2+yhatf**2+zhatf**2)
c      xhatf=xhatf/anorm
c      yhatf=yhatf/anorm
c      zhatf=zhatf/anorm
c     write(*,*)xhatf,yhatf,zhatf,sqrt(xhatf**2+yhatf**2+zhatf**2)

c      psif=dacos(sqrt(yhatf**2+zhatf**2))
c      if(xhatf.lt.0.0) psif=-1.0*psif
c      tthetaf=yhatf/zhatf

c      xhat1f=xhatf*dcos(psif)-
c     &     yhatf*dsin(datan(tthetaf))*dsin(psif)-
c     &     zhatf*dcos(datan(tthetaf))*dsin(psif)
c      yhat1f=
c     &     yhatf*dcos(datan(tthetaf))-
c     &     zhatf*dsin(datan(tthetaf))
c      zhat1f=xhatf*dsin(psif)+
c     &     yhatf*dsin(datan(tthetaf))*dcos(psif)+
c     &     zhatf*dcos(datan(tthetaf))*dcos(psif)
c      write(*,*)xhat1f,yhat1f,zhat1f,sqrt(xhat1f**2+yhat1f**2+zhat1f**2)

c      pp=pout/1000.0-1.4
c      r=pout/1000.0
      
 100  format(1x,3(f8.3,1x))
c
8876  call gheish
c
c      pout=vect(7)*1000.
c      px=pout*vect(4)
c      py=pout*vect(5)
c      pz=pout*vect(6)
c      xhatr=px/pout
c      yhatr=py/pout
c      zhatr=pz/pout
c      anorm=sqrt(xhatr**2+yhatr**2+zhatr**2)
c      xhatr=xhatr/anorm
c      yhatr=yhatr/anorm
c      zhatr=zhatr/anorm

c      xhat1r=xhatr*dcos(psif)-
c     &     yhatr*dsin(datan(tthetaf))*dsin(psif)-
c     &     zhatr*dcos(datan(tthetaf))*dsin(psif)
c      yhat1r=
c     &     yhatr*dcos(datan(tthetaf))-
c     &     zhatr*dsin(datan(tthetaf))
c      zhat1r=xhatr*dsin(psif)+
c     &     yhatr*dsin(datan(tthetaf))*dcos(psif)+
c     &     zhatr*dcos(datan(tthetaf))*dcos(psif)
c      write(*,*)xhat1r,yhat1r,zhat1r
c
 
c      psisc=dasin(xhat1r)
c      if(dcos(psisc).ne.0.0D0) then 
c         thetasc=dasin(yhat1r/dcos(psisc))
c      else
cc         write(*,*)'Zero problem with psisc'
cc         theta_az=100.0
cc         phi_az=400.0
c         goto 3344
c      endif
c      if(dcos(thetasc).ne.0.0D0) then 
c         phisc=datan(dtan(psisc)/dcos(thetasc))
c      else
cc         write(*,*)'Zero problem with thetasc'
cc         theta_az=100.0
cc         phi_az=400.0
cc         goto 3344
c      endif

c      theta_sc=thetasc
c      phi_sc=phisc

cc
cc Another method of changing to spherical coordinates
cc
c      psi_sc=datan(dtan(phi_sc)*dcos(theta_sc))
c      x_sc=dsin(psi_sc)
c      y_sc=dcos(psi_sc)*dsin(theta_sc)
c      z_sc=dcos(psi_sc)*dcos(theta_sc)
c      r_sc=sqrt(x_sc**2+y_sc**2)
c      if(z_sc.ne.0.0.and.y_sc.ne.0.0) then
c         theta_az=datan(r_sc/z_sc)
c         phi_az=datan(x_sc/y_sc)
c      else
cc         write(*,*)'Zero problem with zsc or ysc'
cc         theta_az=100.0
cc         phi_az=400.0
c         goto 3344         
c      endif
c      theta_az=theta_az*180.0/pi
c      phi_az=phi_az*180.0/pi
c      if(y_sc.lt.0.0) phi_az=phi_az+180.0
c      if(y_sc.ge.0.0.and.x_sc.lt.0.0) phi_az=phi_az+360.0
c      geant_theta=theta_az
c      geant_phi=phi_az
c
c      a=acoeff(1,1)+acoeff(1,2)*pp+acoeff(1,3)*pp**2+
c     >   acoeff(1,4)*pp**3+acoeff(1,5)*pp**4+acoeff(1,6)*pp**5
c      b=acoeff(2,1)+acoeff(2,2)*pp+acoeff(2,3)*pp**2+
c     >   acoeff(2,4)*pp**3+acoeff(2,5)*pp**4+acoeff(2,6)*pp**5
c      c=acoeff(3,1)+acoeff(3,2)*pp+acoeff(3,3)*pp**2+
c     >   acoeff(3,4)*pp**3+acoeff(3,5)*pp**4+acoeff(3,6)*pp**5
c      d=acoeff(4,1)+acoeff(4,2)*pp+acoeff(4,3)*pp**2+
c     >   acoeff(4,4)*pp**3+acoeff(4,5)*pp**4+acoeff(4,6)*pp**5

c      p=r
c      r=r*dsin(theta_az/180.0*pi)
c      term1=a*r
c      term2=b*r**2
c      term3=c*r**4
c      term4=d*p*dsin(5.0*theta_az/180.0*pi)
c      ac=term1/(1+term2+term3)+term4
c
cc Now we have the values of phi, ac(theta) and pn,pt.  Now calculate
cc the relative probability.
cc
c      pmax1=1+pn*ac*dcos(datan(+1.0*pt/pn))+
c     >     pt*ac*dsin(datan(+1.0*pt/pn))
c      pmax2=1+pn*ac*dcos(datan(+1.0*pt/pn)+pi)+
c     >     pt*ac*dsin(datan(+1.0*pt/pn)+pi)
c      pmax=max(pmax1,pmax2)
c      prob=1+pn*ac*dcos(phi_az*pi/180.0)+
c     >	         pt*ac*dsin(phi_az*pi/180.0)
c      ratio=prob/pmax
cc      write(*,*)pn/0.384/(-0.87524),pt/(0.384)
c      call grndm(rndm,3)
cc      write(*,*)rndm(1)
c      if(rndm(1).gt.ratio) then
cc	  write(*,*)'Regeneration ...'
c          goto 8876
c      endif

c 3344 continue
cc      write(*,*)'Returning from guhadr ...' 
      return
      end