subroutine gustep
c
c     this subroutine was written by jpsullivan april 21-22, 1993
c     the tracking realated part is relatively simple -- if the
c     particle leave the volume called 'targ', throw it away.
c     it also makes a bunch of histograms
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
c     
      integer nmec,lmec,namec,nstep ,maxnst,ignext,inwvol,istop
     +     ,idecad,iekbin,ilosl, imull,ingoto,nldown,nlevin
     +     ,nlvsav,istory
      real  vect,getot,gekin,vout,destep,destel,safety,sleng ,step
     +     ,snext,sfield,tofg  ,gekrat,upwght
c     end gctrak
*     keep,gcvolu.
*--   author :
      common/gcvolu/nlevel,names(15),number(15),
     +     lvolum(15),lindex(15),infrom,nlevmx,nldev(15),linmx(15),
     +     gtran(3,15),grmat(10,15),gonly(15),glx(3)
c     
      integer nlevel,number,lvolum,lindex,infrom,nlevmx,
     +     nldev,linmx
      character*4 names
      real gtran,grmat,gonly,glx
c     end gcvolu
c     
*     keep,gcbank.
*--   author :
      integer iq,lq,nzebra,ixstor,ixdiv,ixcons,lmain,lr1
      integer kwbank,kwwork,iws
      real gversn,zversn,fendq,ws,q
c     
      parameter (kwbank=69000,kwwork=5200)
      common/gcbank/nzebra,gversn,zversn,ixstor,ixdiv,ixcons,fendq(16)
     +     ,lmain,lr1,ws(kwbank)
      dimension iq(2),q(2),lq(8000),iws(2)
      equivalence (q(1),iq(1),lq(9)),(lq(1),lmain),(iws(1),ws(1))
      common/gclink/jdigi ,jdraw ,jhead ,jhits ,jkine ,jmate ,jpart
     +     ,jrotm ,jrung ,jset  ,jstak ,jgstat,jtmed ,jtrack,jvertx
     +     ,jvolum,jxyz  ,jgpar ,jgpar2,jsklt
c     
      integer       jdigi ,jdraw ,jhead ,jhits ,jkine ,jmate ,jpart
     +     ,jrotm ,jrung ,jset  ,jstak ,jgstat,jtmed ,jtrack,jvertx
     +     ,jvolum,jxyz  ,jgpar,jgpar2 ,jsklt
c     
*     keep,gcking.
*--   author :
      common/gcking/kcase,ngkine,gkin(5,100),tofd(100),iflgk(100)
      integer       kcase,ngkine ,iflgk
      real          gkin,tofd
c     end gcking
c     
*     keep,gckine.
*--   author :
*--   author :
      integer       ikine,itra,istak,ivert,ipart,itrtyp,napart,ipaold
      real          pkine,amass,charge,tlife,vert,pvert
      common/gckine/ikine,pkine(10),itra,istak,ivert,ipart,itrtyp
     +     ,napart(5),amass,charge,tlife,vert(3),pvert(4),ipaold
c     end gckine
c     
      integer ihset,ihdet,iset,idet,idtype,nvname,numbv
      common/gcsets/ihset,ihdet,iset,idet,idtype,nvname,numbv(20)
      real x1,y1,z1,lpar,v1,v2,v3,newdist,x1new,y1new,z1new
      real xstr,ystr,zstr,xstrnew,ystrnew,zstrnew
      real rotmat2,rotmat3,rotmat4,rotmat1
c
      common/geomstep/rotmat1(3,3),rotmat2(3,3),rotmat3(3,3),
     &     rotmat4(3,3)
c     
      include 'fpp_local.h'
      include 'geant_local.h'
c      include 'parameter.h'
c      include 'espace_type.h'
c      include 'detector.h'
c      include 'transport.h'
c      include 'option.h'
c     
c     
      integer i,make_hist,ioff,ihit,ieffcheck
      real pt_pi,ppar_pi,arg1,arg2,rapid_pi,pchmb,hits(6)
      real a,b,c,beta,z,y,straw,ypath,rndm(3),ycompare
c      write(6,*)'entering gustep'
c      write(6,*)'inwvol =',inwvol
c      write(6,*)'position =',vect(1),vect(2),vect(3)
c      write(6,*)'names =',names(nlevel)
c     
c
c
      if ( ngkine.gt.0. ) then
         do i=1,nmec
            if ( lmec(i).eq.12 ) then
c               write ( 6,* ) ' gustep: hadronic interaction'
c               write ( 6,* ) ' nevent=',nevent
            end if
         end do
         mylast = min(100,ngkine)
         do i=1,mylast
            iflgk(i) = 1
            if ( gkin(5,i).eq.4 ) iflgk(i) = 0
            if ( gkin(4,i).gt.0.001 ) iflgk(i)=0
         end do
c         n_2nd = n_2nd + ngkine
      endif
c     
c  the following call makes sure all of the secondary particles
c  get tracked too (provided the flag iflgk(i) for that particle
c  was set in the loop above -- this point is not correctly or
c  clearly documented in the version of the geant manual i have).
c
c      if(sectrack) then
c         call gsking ( 0 )
c      endif
c     
c      write(*,*)'In GUSTEP: ... names(nlevel) = ',names(nlevel)

      make_hist = 0
      if ( inwvol.eq.1 .and. names(nlevel).eq."hall" ) then
         make_hist=1
c     
c     if we get here, the tracking is done for this track
c     (istop=1) but make a bunch of histograms before exitting
c     note that ipart=8 means a pi+ and 9 is a pi-
c     
c     istop = 1
      else if ( istop.ne.0.and. names(nlevel).eq."aira" ) then
         make_hist=0
      else if ( istop.ne.0.and. names(nlevel).eq."airb" ) then
         make_hist=0
      else if ( istop.ne.0.and. names(nlevel).eq."airc" ) then
         make_hist=0
      else if ( istop.ne.0.and. names(nlevel).eq."aird" ) then
         make_hist=0
      else if ( istop.ne.0.and. names(nlevel).eq."aire" ) then
         make_hist=0
      else if ( istop.ne.0.and. names(nlevel).eq."airf" ) then
         make_hist=0
      else if ( names(nlevel).eq."airg" ) then
c         write(*,*)'In airg ... inwvol = ',inwvol
c         write(*,*)'Z-value = ',vect(3)
         if (istop.ne.0) then
            make_hist=0
         endif
      else if ( names(nlevel).eq."HALL" ) then
c         write(*,*)'In HALL ... inwvol = ',inwvol
c         write(*,*)'Z-value = ',vect(3)
         if (istop.ne.0) then
            make_hist=0
         endif
      else if ( istop.ne.0.and. names(nlevel).eq."airh" ) then
         make_hist=0
      else if ( istop.ne.0.and. names(nlevel).eq."hch1" ) then
         make_hist=0
      else if ( istop.ne.0.and. names(nlevel).eq."hch2" ) then
         make_hist=0
      else if ( names(nlevel).eq."fch1" ) then
c         write(*,*)'In fch1 ... inwvol = ',inwvol
         if(inwvol.eq.1) then
c           write(6,*)'Coordinates at fch1'
c           write(6,*)'x=',vect(1),' y=',vect(2),' z=',vect(3)
           x1a=vect(1)
           y1a=vect(2)
           z1a=vect(3)
         endif
         if(inwvol.eq.2) then
            x1b=vect(1)
            y1b=vect(2)
            z1b=vect(3)

            call get_wire_numbers(x1a,y1a,z1a,x1b,y1b,z1b,n1u,n1x,n1v)

            call get_drift_distance(x1a,y1a,z1a,x1b,y1b,z1b,
     &           n1u,n1x,n1v,d1u,d1x,d1v)

c            write(*,*)'Hit in first chamber ...'
c            write(*,*)'Wire Numbers: ',n1u,n1x,n1v

         endif

         if ( istop.ne.0 ) then
            make_hist=0
         endif
      else if ( names(nlevel).eq."fch2" ) then
c         write(*,*)'In fch2 ... inwvol = ',inwvol
c         write(*,*)'Z-value = ',vect(3)
         if(inwvol.eq.1) then
           x2a=vect(1)
           y2a=vect(2)
           z2a=vect(3)
         endif
         if(inwvol.eq.2) then
            x2b=vect(1)
            y2b=vect(2)
            z2b=vect(3)

            call get_wire_numbers(x2a,y2a,z2a,x2b,y2b,z2b,n2u,n2x,n2v)

            call get_drift_distance(x2a,y2a,z2a,x2b,y2b,z2b,
     &           n2u,n2x,n2v,d2u,d2x,d2v)

c            write(*,*)'Hit in second chamber ...'
c            write(*,*)'Wire Numbers: ',n2u,n2x,n2v

         endif
c
         if ( istop.ne.0 ) then
            make_hist=0
         endif
      else if ( names(nlevel).eq."fch3" ) then
         if(inwvol.eq.1) then
           x3a=vect(1)
           y3a=vect(2)
           z3a=vect(3)
         endif
         if(inwvol.eq.2) then
            x3b=vect(1)
            y3b=vect(2)
            z3b=vect(3)

            call get_wire_numbers(x3a,y3a,z3a,x3b,y3b,z3b,n3u,n3x,n3v)

            call get_drift_distance(x3a,y3a,z3a,x3b,y3b,z3b,
     &           n3u,n3x,n3v,d3u,d3x,d3v)

c            write(*,*)'Hit in third chamber ...'
c            write(*,*)'Wire Numbers: ',n3u,n3x,n3v

         endif
c
         if ( istop.ne.0 ) then
            make_hist=0
         endif     
      else if ( names(nlevel).eq."fch4" ) then
         if(inwvol.eq.1) then
           x4a=vect(1)
           y4a=vect(2)
           z4a=vect(3)
         endif
         if(inwvol.eq.2) then
            x4b=vect(1)
            y4b=vect(2)
            z4b=vect(3)

            call get_wire_numbers(x4a,y4a,z4a,x4b,y4b,z4b,n4u,n4x,n4v)

            call get_drift_distance(x4a,y4a,z4a,x4b,y4b,z4b,
     &           n4u,n4x,n4v,d4u,d4x,d4v)

c            write(*,*)'Hit in fourth chamber ...'
c            write(*,*)'Wire Numbers: ',n4u,n4x,n4v

         endif
c
         if ( istop.ne.0 ) then
            make_hist=0
         endif
      else if ( istop.ne.0.and. names(nlevel).eq."sci1" ) then
         make_hist=0
      else if ( names(nlevel).eq."anl1" ) then
         if(inwvol.eq.1) then
c           write(6,*)'We have a hit in the fist analyzer at'
c           write(6,*)'x=',vect(1),' y=',vect(2),' z=',vect(3)
           xdet=vect(1)
           ydet=vect(2)
           zdet=vect(3)
         endif
         if ( istop.ne.0 ) then
            make_hist=0
         endif
      end if
c     
c     
c     store current track parameters (including position ) in jxyz structure.
c     
      call gsxyz
c     
c     moved histograming stuff to gulast
c     
c     write(6,*)'done in gustep'
 9999 return
      end

      subroutine get_wire_numbers(xa,ya,za,xb,yb,zb,nu,nx,nv)

      implicit none

      real*8 xa,ya,za,xb,yb,zb
      integer*4 nu,nx,nv

      real*8 zc,zu,zx,zv,zt,xp,yp,xu,xx,xv,yu,yx,yv,uw,xw,vw
      real*8 anu,anx,anv

c We have the (x,y,z) coordinates of the entrance (a) and exit (b) points
c of the track.  We can use this information to calculate the wire numbers that
c were hit in each plane.
c
      zc=(zb-za)/2.0+za
      zu=zc-1.60
      zx=zc
      zv=zc+1.60
      zt=(zb-za)
      xp=(xb-xa)/zt
      yp=(yb-ya)/zt
      xu=xa+xp*(zu-za)
      yu=ya+yp*(zu-za)
      xx=xa+xp*(zx-za)
      yx=ya+yp*(zx-za)
      xv=xa+xp*(zv-za)
      yv=ya+yp*(zv-za)
c
      uw=(xu+yu)/sqrt(2.0)
      xw=xx
      vw=(-xv+yv)/sqrt(2.0)
c
c      write(*,*)'********************'
c      write(*,*)'A: ',xa,ya,za
c      write(*,*)'B: ',xb,yb,zb
c      write(*,*)'U: ',xu,yu,zu
c      write(*,*)'X: ',xx,yx,zx
c      write(*,*)'V: ',xv,yv,zv
c      write(*,*)'W: ',uw,xw,vw
c      write(*,*)'********************'
c     
      anu=(-uw-3.592+104.0)/2.0
      anx=(-xw-5.080+84.0)/2.0
      anv=(vw-3.592+104.0)/2.0
c
c      write(*,*)'Wires: ',anu,anx,anv
c      write(*,*)'********************'
      nu=anu
      nv=anv
      nx=anx
      if((anu-nu).ge.0.500)nu=nu+1
      if((anx-nx).ge.0.500)nx=nx+1
      if((anv-nv).ge.0.500)nv=nv+1
c
      return
      end

      subroutine get_drift_distance(xa,ya,za,xb,yb,zb,
     &           nu,nx,nv,du,dx,dv)

      implicit none

      real*8 xa,ya,za,xb,yb,zb,du,dx,dv
      integer*4 nu,nx,nv
      real*8 tphi,ttheta
      real*8 uw,xw,vw
      real*8 c11,c12,c21,c22,d1,d2,a,zt,xt,yt
      real*8 xu,yu,zu
      real*8 xv,yv,zv
      real*8 xx,yx,zx

      tphi=(xb-xa)/(zb-za)
      ttheta=(yb-ya)/(zb-za)

      uw=-2.0*nu-3.592+104.0
      xw=-2.0*nx-5.080+84.0
      vw=2.0*nv+3.592-104.0
      zu=(zb+za)/2.0-1.60
      zv=(zb+za)/2.0+1.60
      zx=(zb+za)/2.0

c      write(*,*)'W: ',uw,xw,vw
c      write(*,*)'********************'

      c11=tphi-ttheta
      c12=sqrt(2.0)
      d1=-xa+ya
      c21=tphi*tphi+ttheta*ttheta+1.0
      c22=-ttheta/sqrt(2.0)+tphi/sqrt(2.0)
      d2=uw*(ttheta+tphi)/sqrt(2.0)-xa*tphi-ya*ttheta+zu-za

      a=(d1*c21-d2*c11)/(c21*c12-c22*c11)
      zt=(d1-c12*a)/c11
c      write(*,*)'U Plane zt,a: ',zt,a

      xt=xa+zt*tphi
      yt=ya+zt*ttheta
      xu=(uw-a)/sqrt(2.0)
      yu=(uw+a)/sqrt(2.0)

      zt=zt+za

      du=sqrt((xt-xu)**2+(yt-yu)**2+(zt-zu)**2)      
       
      c11=tphi+ttheta
      c12=-sqrt(2.0)
      d1=-xa-ya
      c21=tphi*tphi+ttheta*ttheta+1.0
      c22=-ttheta/sqrt(2.0)-tphi/sqrt(2.0)
      d2=vw*(ttheta+tphi)/sqrt(2.0)-xa*tphi-ya*ttheta+zv-za

      a=(d1*c21-d2*c11)/(c21*c12-c22*c11)
      zt=(d1-c12*a)/c11
c      write(*,*)'V Plane zt,a: ',zt,a

      xt=xa+zt*tphi
      yt=ya+zt*ttheta
      xv=-(vw-a)/sqrt(2.0)
      yv=(vw+a)/sqrt(2.0)

      zt=zt+za
      
      dv=sqrt((xt-xv)**2+(yt-yv)**2+(zt-zv)**2)
     
      c11=tphi
      c12=-1.0
      d1=-ya
      c21=tphi*tphi+ttheta*ttheta+1.0
      c22=-ttheta
      d2=xw*tphi-xa*tphi-ya*ttheta+zx-za

      a=(d1*c21-d2*c11)/(c21*c12-c22*c11)
      zt=(d1-c12*a)/c11
c      write(*,*)'X Plane zt,a: ',zt,a

      xt=xa+zt*tphi
      yt=ya+zt*ttheta
      xx=xw
      yx=a

      zt=zt+za
      
      dx=sqrt((xt-xx)**2+(yt-yx)**2+(zt-zx)**2)

      if((du.gt.1.00)) then
         write(*,*)'Calculating Drift Distance ...'
         write(*,*)'A: ',xa,ya,za
         write(*,*)'B: ',xb,yb,zb
         write(*,*)'U: ',xu,yu,zu
         write(*,*)'V: ',xv,yv,zv
         write(*,*)'X: ',xx,yx,zx
         write(*,*)'U Drift distance = ',du
         write(*,*)'V Drift distance = ',dv
         write(*,*)'X Drift distance = ',dx
      endif

      return
      end