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 :

      implicit none

      integer*4 mylast
      integer*4 nhu,nhv,nhx

      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
      logical idflag
      real*8 d1uetemp,d1xetemp,d1vetemp
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)
c      write(*,*)'In GUSTEP: ... inwvol = ',inwvol

      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 ( names(nlevel).eq."hch2" ) then
c	 write(*,*)'In hch2'
         if(inwvol.eq.1) then
c           write(6,*)'Coordinates at hch2'
c           write(6,*)'x=',vect(1),' y=',vect(2),' z=',vect(3)
           xahch2=vect(1)
           yahch2=vect(2)
           zahch2=vect(3)
         endif
         if(inwvol.eq.2) then
            xbhch2=vect(1)
            ybhch2=vect(2)
            zbhch2=vect(3)
         endif

         if ( istop.ne.0 ) then
            make_hist=0
         endif
      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,nhu,nhx,nhv,
     &           n1ua,n1xa,n1va,n1ub,n1xb,n1vb)

            nhu1=nhu
            nhx1=nhx
            nhv1=nhv

            if(nhu.eq.1.and.nhv.eq.1.and.nhx.eq.1) then
               idflag=.false.
               call get_drift_distance_ejb(x1a,y1a,z1a,x1b,y1b,z1b,
     &              n1ua,n1xa,n1va,d1ue,d1xe,d1ve,idflag)
               call get_drift_distance(x1a,y1a,z1a,x1b,y1b,z1b,
     &              n1ua,n1xa,n1va,d1u,d1x,d1v)
            else
               call get_drift_distance_ejb(x1a,y1a,z1a,x1b,y1b,z1b,
     &              n1ua,n1xa,n1va,d1ue,d1xe,d1ve,idflag)
               d1uetemp=d1ue
               d1xetemp=d1xe
               d1vetemp=d1ve
               call get_drift_distance_ejb(x1a,y1a,z1a,x1b,y1b,z1b,
     &              n1ub,n1xb,n1vb,d1ue,d1xe,d1ve,idflag)
               if(idflag)
     &             write(*,*)'Drift Distance 1a: ',
     &                 d1uetemp,d1xetemp,d1vetemp
               if(idflag)
     &             write(*,*)'Drift Distance 1b: ',d1ue,d1xe,d1ve
               call get_drift_distance(x1a,y1a,z1a,x1b,y1b,z1b,
     &              n1ua,n1xa,n1va,d1u,d1x,d1v)
               if(idflag)
     &             write(*,*)'Drift Distance 1c: ',d1u,d1x,d1v
               call get_drift_distance(x1a,y1a,z1a,x1b,y1b,z1b,
     &              n1ub,n1xb,n1vb,d1u,d1x,d1v)
               if(idflag)
     &             write(*,*)'Drift Distance 1d: ',d1u,d1x,d1v
            endif
          
c            write(*,*)'Hit in first chamber ...'
c            write(*,*)'Wire Numbers: ',n1ua,n1xa,n1va
               
         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,nhu,nhx,nhv,
     &           n2ua,n2xa,n2va,n2ub,n2xb,n2vb)

            nhu2=nhu
            nhx2=nhx
            nhv2=nhv

            if(nhu.eq.1.and.nhv.eq.1.and.nhx.eq.1) then
               call get_drift_distance_ejb(x2a,y2a,z2a,x2b,y2b,z2b,
     &              n2ua,n2xa,n2va,d2ue,d2xe,d2ve,idflag)
               call get_drift_distance(x2a,y2a,z2a,x2b,y2b,z2b,
     &              n2ua,n2xa,n2va,d2u,d2x,d2v)
            else
               call get_drift_distance_ejb(x2a,y2a,z2a,x2b,y2b,z2b,
     &              n2ua,n2xa,n2va,d2ue,d2xe,d2ve,idflag)
c               write(*,*)'Drift Distance 2a: ',d2ue,d2xe,d2ve
               call get_drift_distance_ejb(x2a,y2a,z2a,x2b,y2b,z2b,
     &              n2ub,n2xb,n2vb,d2ue,d2xe,d2ve,idflag)
c               write(*,*)'Drift Distance 2b: ',d2ue,d2xe,d2ve
               call get_drift_distance(x2a,y2a,z2a,x2b,y2b,z2b,
     &              n2ua,n2xa,n2va,d2u,d2x,d2v)
c               write(*,*)'Drift Distance 2c: ',d2u,d2x,d2v
               call get_drift_distance(x2a,y2a,z2a,x2b,y2b,z2b,
     &              n2ub,n2xb,n2vb,d2u,d2x,d2v)
c               write(*,*)'Drift Distance 2c: ',d2u,d2x,d2v
            endif

           call calc_theta_phi(xahch2,yahch2,zahch2,
     &            xbhch2,ybhch2,zbhch2,
     &           x1a,y1a,z1a,x2b,y2b,z2b,
     &            theta_front,phi_front)

c            write(*,*)'Hit in first chamber ...'
c            write(*,*)'Wire Numbers: ',n2ua,n2xa,n2va

         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,nhu,nhx,nhv,
     &           n3ua,n3xa,n3va,n3ub,n3xb,n3vb)

            nhu3=nhu
            nhx3=nhx
            nhv3=nhv

            if(nhu.eq.1.and.nhv.eq.1.and.nhx.eq.1) then
               call get_drift_distance_ejb(x3a,y3a,z3a,x3b,y3b,z3b,
     &              n3ua,n3xa,n3va,d3ue,d3xe,d3ve,idflag)
               call get_drift_distance(x3a,y3a,z3a,x3b,y3b,z3b,
     &              n3ua,n3xa,n3va,d3u,d3x,d3v)
            else
               call get_drift_distance_ejb(x3a,y3a,z3a,x3b,y3b,z3b,
     &              n3ua,n3xa,n3va,d3ue,d3xe,d3ve,idflag)
c               write(*,*)'Drift Distance 3a: ',d3ue,d3xe,d3ve
               call get_drift_distance_ejb(x3a,y3a,z3a,x3b,y3b,z3b,
     &              n3ub,n3xb,n3vb,d3ue,d3xe,d3ve,idflag)
c               write(*,*)'Drift Distance 3b: ',d3ue,d3xe,d3ve
               call get_drift_distance(x3a,y3a,z3a,x3b,y3b,z3b,
     &              n3ua,n3xa,n3va,d3u,d3x,d3v)
c               write(*,*)'Drift Distance 3c: ',d3u,d3x,d3v
            endif
          

c            write(*,*)'Hit in first chamber ...'
c            write(*,*)'Wire Numbers: ',n3ua,n3xa,n3va

         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,nhu,nhx,nhv,
     &           n4ua,n4xa,n4va,n4ub,n4xb,n4vb)

            nhu4=nhu
            nhx4=nhx
            nhv4=nhv

            if(nhu.eq.1.and.nhv.eq.1.and.nhx.eq.1) then
               call get_drift_distance_ejb(x4a,y4a,z4a,x4b,y4b,z4b,
     &              n4ua,n4xa,n4va,d4ue,d4xe,d4ve,idflag)
               call get_drift_distance(x4a,y4a,z4a,x4b,y4b,z4b,
     &              n4ua,n4xa,n4va,d4u,d4x,d4v)
            else
               call get_drift_distance_ejb(x4a,y4a,z4a,x4b,y4b,z4b,
     &              n4ua,n4xa,n4va,d4ue,d4xe,d4ve,idflag)
c               write(*,*)'Drift Distance 4a: ',d4ue,d4xe,d4ve
               call get_drift_distance_ejb(x4a,y4a,z4a,x4b,y4b,z4b,
     &              n4ub,n4xb,n4vb,d4ue,d4xe,d4ve,idflag)
c               write(*,*)'Drift Distance 4b: ',d4ue,d4xe,d4ve
               call get_drift_distance(x4a,y4a,z4a,x4b,y4b,z4b,
     &              n4ua,n4xa,n4va,d4u,d4x,d4v)
c               write(*,*)'Drift Distance 4c: ',d4u,d4x,d4v
            endif
          

c            write(*,*)'Hit in first chamber ...'
c            write(*,*)'Wire Numbers: ',n3ua,n3xa,n3va

         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,nhu,nhx,nhv,
     %     nu1,nx1,nv1,nu2,nx2,nv2)

      implicit none

      real*8 xa,ya,za,xb,yb,zb
      integer*4 nu1,nx1,nv1
      integer*4 nu2,nx2,nv2
      integer*4 nhu,nhx,nhv

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

      nu2=0
      nx2=0
      nv2=0
c
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
c
c Project to the FRONT of the "cell" associated with each plane.
c
      xu=xa+xp*(zu-za-0.8)
      yu=ya+yp*(zu-za-0.8)
      xx=xa+xp*(zx-za-0.8)
      yx=ya+yp*(zx-za-0.8)
      xv=xa+xp*(zv-za-0.8)
      yv=ya+yp*(zv-za-0.8)      
c
c      xu=xa
c      yu=ya
c      xx=xa
c      yx=ya
c      xv=xa
c      yv=ya
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(*,*)'********************'
      nu1=anu
      nv1=anv
      nx1=anx
      if((anu-nu1).ge.0.500)nu1=nu1+1
      if((anx-nx1).ge.0.500)nx1=nx1+1
      if((anv-nv1).ge.0.500)nv1=nv1+1
c      write(*,*)'using front of chamber: ',nu1,nx1,nv1
c
c Now project to the BACK of the "cell" associated with each plane.
c
      xu=xa+xp*(zu-za+0.8)
      yu=ya+yp*(zu-za+0.8)
      xx=xa+xp*(zx-za+0.8)
      yx=ya+yp*(zx-za+0.8)
      xv=xa+xp*(zv-za+0.8)
      yv=ya+yp*(zv-za+0.8)      
c
c      xu=xb
c      yu=yb
c      xx=xb
c      yx=yb
c      xv=xb
c      yv=yb
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(*,*)'********************'
      nu2=anu
      nv2=anv
      nx2=anx
      if((anu-nu2).ge.0.500)nu2=nu2+1
      if((anx-nx2).ge.0.500)nx2=nx2+1
      if((anv-nv2).ge.0.500)nv2=nv2+1
      nhu=1
      nhx=1
      nhv=1
      if(nu1.ne.nu2)nhu=2
      if(nx1.ne.nx2)nhx=2
      if(nv1.ne.nv2)nhv=2

      return
      end

c      subroutine get_drift_distance_ejb(xa,ya,za,xb,yb,zb,
c     &           nu,nx,nv,du,dx,dv)
c
c      implicit none
c
c      real*8 xa,ya,za,xb,yb,zb,du,dx,dv
c      integer*4 nu,nx,nv
c      real*8 tphi,ttheta
c      real*8 uw,xw,vw
c      real*8 c11,c12,c21,c22,d1,d2,a,zt,xt,yt
c      real*8 xu,yu,zu
c      real*8 xv,yv,zv
c      real*8 xx,yx,zx
c
c      tphi=(xb-xa)/(zb-za)
c      ttheta=(yb-ya)/(zb-za)
c
c      uw=-2.0*nu-3.592+104.0
c      xw=-2.0*nx-5.080+84.0
c      vw=2.0*nv+3.592-104.0
c      zu=(zb+za)/2.0-1.60
c      zv=(zb+za)/2.0+1.60
c      zx=(zb+za)/2.0
c
cc      write(*,*)'W: ',uw,xw,vw
cc      write(*,*)'********************'
c
c      c11=tphi-ttheta
c      c12=sqrt(2.0)
c      d1=-xa+ya
c      c21=tphi*tphi+ttheta*ttheta+1.0
c      c22=-ttheta/sqrt(2.0)+tphi/sqrt(2.0)
c      d2=uw*(ttheta+tphi)/sqrt(2.0)-xa*tphi-ya*ttheta+zu-za
c
c      a=(d1*c21-d2*c11)/(c21*c12-c22*c11)
c      zt=(d1-c12*a)/c11
cc      write(*,*)'U Plane zt,a: ',zt,a
c
c      xt=xa+zt*tphi
c      yt=ya+zt*ttheta
c      xu=(uw-a)/sqrt(2.0)
c      yu=(uw+a)/sqrt(2.0)
c
c      zt=zt+za
c
c      du=sqrt((xt-xu)**2+(yt-yu)**2+(zt-zu)**2)      
c       
c      c11=tphi+ttheta
c      c12=-sqrt(2.0)
c      d1=-xa-ya
c      c21=tphi*tphi+ttheta*ttheta+1.0
c      c22=-ttheta/sqrt(2.0)-tphi/sqrt(2.0)
c      d2=vw*(ttheta+tphi)/sqrt(2.0)-xa*tphi-ya*ttheta+zv-za
c
c      a=(d1*c21-d2*c11)/(c21*c12-c22*c11)
c      zt=(d1-c12*a)/c11
cc      write(*,*)'V Plane zt,a: ',zt,a
c
c      xt=xa+zt*tphi
c      yt=ya+zt*ttheta
c      xv=-(vw-a)/sqrt(2.0)
c      yv=(vw+a)/sqrt(2.0)
c
c      zt=zt+za
c      
c      dv=sqrt((xt-xv)**2+(yt-yv)**2+(zt-zv)**2)
c     
c      c11=tphi
c      c12=-1.0
c      d1=-ya
c      c21=tphi*tphi+ttheta*ttheta+1.0
c      c22=-ttheta
c      d2=xw*tphi-xa*tphi-ya*ttheta+zx-za
c
c      a=(d1*c21-d2*c11)/(c21*c12-c22*c11)
c      zt=(d1-c12*a)/c11
cc      write(*,*)'X Plane zt,a: ',zt,a
c
c      xt=xa+zt*tphi
c      yt=ya+zt*ttheta
c      xx=xw
c      yx=a
c
c      zt=zt+za
c      
c      dx=sqrt((xt-xx)**2+(yt-yx)**2+(zt-zx)**2)
c
cc      if((du.gt.1.00)) then
c         write(*,*)'Calculating Drift Distance ...'
c         write(*,*)'A: ',xa,ya,za
c         write(*,*)'B: ',xb,yb,zb
cc         write(*,*)'U: ',xu,yu,zu
cc         write(*,*)'V: ',xv,yv,zv
cc         write(*,*)'X: ',xx,yx,zx
c         write(*,*)'U Drift distance = ',du
c         write(*,*)'V Drift distance = ',dv
c         write(*,*)'X Drift distance = ',dx
cc      endif
c
c      return
c      end


      subroutine get_drift_distance(xa,ya,za,xb,yb,zb,
     &           nu,nx,nv,distu,distx,distv)
c
c     This subroutine uses the entrance (a) and exit (b)
c     points of a chamber to define the line of the track.
c     It uses the wire number and wire direction to define
c     the wire line.  It then uses these to build parallel
c     planes by computing a normal vector to both of the lines.
c     Finally, it calculates the distance between these two
c     planes which is the distance of shortest approach.
c      
      implicit none
c
      real*8 xa,ya,za,xb,yb,zb
      real*8 distu,distx,distv,xp,yp,zc
      integer*4 nu,nx,nv
      real*8 uw,xw,vw
      real*8 zt,xt,yt
      real*8 xu,yu,zu
      real*8 xv,yv,zv
      real*8 xx,yx,zx
c
c     the direction of each of the lines
c     vect1=track   vectu=u wire
      real*8 vect1(1:3), vectu(1:3)
      real*8 vectx(1:3), vectv(1:3)
c   
c     the normal vector to both lines
      real*8 normu(1:3)
      real*8 normx(1:3)
      real*8 normv(1:3)
c
c     the coefficients of the plane
      real*8 pvectu(1:4)
      real*8 pvectx(1:4)
      real*8 pvectv(1:4)
c
      vect1(1)=xb-xa
      vect1(2)=yb-ya
      vect1(3)=zb-za
c
      zu=(zb+za)/2.0-1.60
      zv=(zb+za)/2.0+1.60
      zx=(zb+za)/2.0
c
c     use line number to calculate distance relative to
c     wire plane, then convert to x and y
c      write(*,*)'nu nx nv',nu,nx,nv
      uw=-2.0*nu-3.592+104.0
      xw=-2.0*nx-5.080+84.0
      vw=2.0*nv+3.592-104.0
      xu=uw/sqrt(2.0)
      yu=uw/sqrt(2.0)      
      xx=xw
      yx=0
      xv=-vw/sqrt(2.0)
      yv=vw/sqrt(2.0)
c      write(*,*)'uw xu yu zu',uw,xu,yu,zu
c      write(*,*)'xw xx yx zx',xw,xx,yx,zx
c      write(*,*)'vw xv yv zv',vw,xv,yv,zv
c
c     define direction vector for wires, will be the same
c     for each wire in a given plane, and is known
c     for each plane
      vectu(1)=1.0/sqrt(2.0)
      vectu(2)=1.0/sqrt(2.0)
      vectu(3)=0.0
      vectx(1)=0.0
      vectx(2)=1.0
      vectx(3)=0.0
      vectv(1)=-1.0/sqrt(2.0)
      vectv(2)=1.0/sqrt(2.0)
      vectv(3)=0.0
c
c      write(*,*)'distance calculations .....'
c      write(*,*)xa,ya,za
c      write(*,*)vect1(1),vect1(2),vect1(3)
c      write(*,*)xu,yu,zu
c      write(*,*)vectu(1),vectu(2),vectu(3)
c      write(*,*)xx,yx,zx
c      write(*,*)vectx(1),vectx(2),vectx(3)
c      write(*,*)xv,yv,zv
c      write(*,*)vectv(1),vectv(2),vectv(3)
c      write(*,*)'distance calculations .....'
c
c     cross product
      normu(1)=vect1(2)*vectu(3)-vect1(3)*vectu(2)
      normu(2)=vect1(1)*vectu(3)-vect1(3)*vectu(1)
      normu(3)=vect1(1)*vectu(2)-vect1(2)*vectu(1)
      normx(1)=vect1(2)*vectx(3)-vect1(3)*vectx(2)
      normx(2)=vect1(1)*vectx(3)-vect1(3)*vectx(1)
      normx(3)=vect1(1)*vectx(2)-vect1(2)*vectx(1)
      normv(1)=vect1(2)*vectv(3)-vect1(3)*vectv(2)
      normv(2)=vect1(1)*vectv(3)-vect1(3)*vectv(1)
      normv(3)=vect1(1)*vectv(2)-vect1(2)*vectv(1)
c      
      pvectu(1)=normu(1)
      pvectu(2)=normu(2)
      pvectu(3)=normu(3)
      pvectu(4)=normu(1)*(-xu)+normu(2)*(-yu)+normu(3)*(-zu)
      pvectx(1)=normx(1)
      pvectx(2)=normx(2)
      pvectx(3)=normx(3)
      pvectx(4)=normx(1)*(-xx)+normx(2)*(-yx)+normx(3)*(-zx)
      pvectv(1)=normv(1)
      pvectv(2)=normv(2)
      pvectv(3)=normv(3)
      pvectv(4)=normv(1)*(-xv)+normv(2)*(-yv)+normv(3)*(-zv)
c
c     distance formula
      distu=(pvectu(1)*xa+pvectu(2)*ya+pvectu(3)*za+pvectu(4))
     &     /sqrt(normu(1)**2+normu(2)**2+normu(3)**2)
      distx=(pvectx(1)*xa+pvectx(2)*ya+pvectx(3)*za+pvectx(4))
     &     /sqrt(normx(1)**2+normx(2)**2+normx(3)**2)
      distv=(pvectv(1)*xa+pvectv(2)*ya+pvectv(3)*za+pvectv(4))
     &     /sqrt(normv(1)**2+normv(2)**2+normv(3)**2)
c      write(*,*)'Drift distance: ',distu, distx, distv
c     
c        write(*,*)'Brads routine ....' 
c        write(*,*)xa,ya,za
c        write(*,*)xb,yb,zb
c        write(*,*)nu,nx,nv
c   	write(*,*)uw,xw,vw
c        write(*,*)'Drift distance: ',distu, distx, distv
      
      return
      end

      subroutine get_drift_distance_ejb(xa,ya,za,xb,yb,zb,
     &           nu,nx,nv,distu,distx,distv,idflag)
c
c Author: Ed Brash - December 15th, 2005
c Yet another attempt at a full drift distance calculation
c
      implicit none
c
      real*8 xa,ya,za,xb,yb,zb
      real*8 distu,distx,distv,xp,yp,zc
      integer*4 nu,nx,nv
      real*8 uw,xw,vw
      real*8 zt,xt,yt
      real*8 xu,yu,zu
      real*8 xv,yv,zv
      real*8 xx,yx,zx
      logical idflag
c
c     the direction of each of the lines
c     vect1=track   vectu=u wire
      real*8 vect1(1:3), vectu(1:3)
      real*8 vectx(1:3), vectv(1:3)

c     the difference vector between the defining points
      real*8 du(1:3),dx(1:3),dv(1:3)
c   
c     the normal vector to both lines
      real*8 normu(1:3)
      real*8 normx(1:3)
      real*8 normv(1:3)
      real*8 normumag,normxmag,normvmag
c
c     the coefficients of the distance vector
      real*8 dvectu(1:4)
      real*8 dvectx(1:4)
      real*8 dvectv(1:4)
c
      idflag=.false.
      vect1(1)=xb-xa
      vect1(2)=yb-ya
      vect1(3)=zb-za
c
      zu=(zb+za)/2.0-1.60
      zv=(zb+za)/2.0+1.60
      zx=(zb+za)/2.0
c
c     use line number to calculate distance relative to
c     wire plane, then convert to x and y
c      write(*,*)'nu nx nv',nu,nx,nv
      uw=-2.0*nu-3.592+104.0
      xw=-2.0*nx-5.080+84.0
      vw=2.0*nv+3.592-104.0
      xu=uw/sqrt(2.0)
      yu=uw/sqrt(2.0)      
      xx=xw
      yx=0
      xv=-vw/sqrt(2.0)
      yv=vw/sqrt(2.0)
c      write(*,*)'uw xu yu zu',uw,xu,yu,zu
c      write(*,*)'xw xx yx zx',xw,xx,yx,zx
c      write(*,*)'vw xv yv zv',vw,xv,yv,zv
c
c     define direction vector for wires, will be the same
c     for each wire in a given plane, and is known
c     for each plane
      vectu(1)=1.0/sqrt(2.0)
      vectu(2)=-1.0/sqrt(2.0)
      vectu(3)=0.0
      vectx(1)=0.0
      vectx(2)=1.0
      vectx(3)=0.0
      vectv(1)=1.0/sqrt(2.0)
      vectv(2)=1.0/sqrt(2.0)
      vectv(3)=0.0
c
c      write(*,*)'distance calculations .....'
c      write(*,*)xa,ya,za
c      write(*,*)vect1(1),vect1(2),vect1(3)
c      write(*,*)xu,yu,zu
c      write(*,*)vectu(1),vectu(2),vectu(3)
c      write(*,*)xx,yx,zx
c      write(*,*)vectx(1),vectx(2),vectx(3)
c      write(*,*)xv,yv,zv
c      write(*,*)vectv(1),vectv(2),vectv(3)
c      write(*,*)'distance calculations .....'
c
c     cross product
      normu(1)=vect1(2)*vectu(3)-vect1(3)*vectu(2)
      normu(2)=vect1(3)*vectu(1)-vect1(1)*vectu(3)
      normu(3)=vect1(1)*vectu(2)-vect1(2)*vectu(1)
      normx(1)=vect1(2)*vectx(3)-vect1(3)*vectx(2)
      normx(2)=vect1(3)*vectx(1)-vect1(1)*vectx(3)
      normx(3)=vect1(1)*vectx(2)-vect1(2)*vectx(1)
      normv(1)=vect1(2)*vectv(3)-vect1(3)*vectv(2)
      normv(2)=vect1(3)*vectv(1)-vect1(1)*vectv(3)
      normv(3)=vect1(1)*vectv(2)-vect1(2)*vectv(1)
c      write(*,*)normu(1),normu(2),normu(3)
c
      normumag=sqrt(normu(1)**2+normu(2)**2+normu(3)**2)
      normxmag=sqrt(normx(1)**2+normx(2)**2+normx(3)**2)
      normvmag=sqrt(normv(1)**2+normv(2)**2+normv(3)**2)
      normu(1)=normu(1)/normumag
      normu(2)=normu(2)/normumag
      normu(3)=normu(3)/normumag
      normx(1)=normx(1)/normxmag
      normx(2)=normx(2)/normxmag
      normx(3)=normx(3)/normxmag
      normv(1)=normv(1)/normvmag
      normv(2)=normv(2)/normvmag
      normv(3)=normv(3)/normvmag
c      write(*,*)normumag
c
      du(1)=xa-xu
      du(2)=ya-yu
      du(3)=za-zu      
      dx(1)=xa-xx
      dx(2)=ya-yx
      dx(3)=za-zx      
      dv(1)=xa-xv
      dv(2)=ya-yv
      dv(3)=za-zv
c
c
c     distance formula
      distu=du(1)*normu(1)+du(2)*normu(2)+du(3)*normu(3)     
      distx=dx(1)*normx(1)+dx(2)*normx(2)+dx(3)*normx(3)     
      distv=dv(1)*normv(1)+dv(2)*normv(2)+dv(3)*normv(3)     
c
      if(distu.gt.1.0.or.distx.gt.1.0.or.distv.gt.1.0)
     &            idflag=.true. 
      if(distu.gt.1.28.or.distx.gt.1.28.or.distv.gt.1.28)then
        write(*,*)'Problem Child !!!'
        write(*,*)'distance calculations .....'
        write(*,*)xa,ya,za
c        write(*,*)vect1(1),vect1(2),vect1(3)
        write(*,*)xu,yu,zu
c        write(*,*)vectu(1),vectu(2),vectu(3)
        write(*,*)xx,yx,zx
c        write(*,*)vectx(1),vectx(2),vectx(3)
        write(*,*)xv,yv,zv
c        write(*,*)vectv(1),vectv(2),vectv(3)
        write(*,*)'normalization factors'
        write(*,*)normu(1),normu(2),normu(3)
        write(*,*)normumag
        write(*,*)normx(1),normx(2),normx(3)
        write(*,*)normxmag
        write(*,*)normv(1),normv(2),normv(3)
        write(*,*)normvmag
        write(*,*)'Drift distance: ',distu, distx, distv
      endif
c      
      return
      end


      subroutine calc_theta_phi(xin1,yin1,zin1,xin2,yin2,zin2,
     &     xsc1,ysc1,zsc1,xsc2,ysc2,zsc2,theta,phi)
c
      implicit none
      include 'fpp_local.h'
      include 'geant_local.h'
c
      real*8 xin1,yin1,zin1,xin2,yin2,zin2
      real*8 xsc1,ysc1,zsc1,xsc2,ysc2,zsc2,theta,phi
      real*8 ftheta, fphi, fpsi
      real*8 lin,lout,theta_ejb,phi_ejb
c
      real invect(1:3)
      real scvect(1:3)
      real scvect2(1:3)
      real in(1:3)
      real out(1:3)
      real scat(1:3)
c
      invect(1)=xin2-xin1
      invect(2)=yin2-yin1
      invect(3)=zin2-zin1
      scvect(1)=xsc2-xsc1
      scvect(2)=ysc2-ysc1
      scvect(3)=zsc2-zsc1
c      write(*,*)'INCOMING: ',invect(1),invect(2),invect(3)
c      write(*,*)'SCATTERED: ',scvect(1),scvect(2),scvect(3)
c
c EJB calculation of theta and phi
c
      in(1)=invect(1)/invect(3)
      in(2)=invect(2)/invect(3)
      in(3)=invect(3)/invect(3)
      out(1)=scvect(1)/scvect(3)
      out(2)=scvect(2)/scvect(3)
      out(3)=scvect(3)/scvect(3)
      lin=sqrt(in(1)**2+in(2)**2+in(3)**2)
      lout=sqrt(out(1)**2+out(2)**2+out(3)**2)
      scat(1)=out(1)-in(1)
      scat(2)=out(2)-in(2)
      scat(3)=out(3)
      x_ejb=scat(1)
      y_ejb=scat(2)
      z_ejb=scat(3)
      if(scat(1).ge.0.0.and.scat(2).gt.0.0)then
         phi_ejb=atan(scat(1)/scat(2))*57.2957795
      else if(scat(1).ge.0.0.and.scat(2).lt.0.0)then
         phi_ejb=atan(scat(1)/scat(2))*57.2957795+180.00
      else if(scat(1).le.0.0.and.scat(2).lt.0.0)then
         phi_ejb=atan(scat(1)/scat(2))*57.2957795+180.00
      else if(scat(1).le.0.0.and.scat(2).gt.0.0)then
         phi_ejb=atan(scat(1)/scat(2))*57.2957795+360.00
      endif
c
      theta_ejb=acos((in(1)*out(1)+in(2)*out(2)+in(3)*out(3))/
     &     (lin*lout))*57.2957795
c      write(*,*)'EJB Incoming Vector = ',in(1),in(2),in(3)
c      write(*,*)'EJB Outgoing Vector = ',out(1),out(2),out(3)
c      write(*,*)'EJB Scattered Vector = ',scat(1),scat(2),scat(3)
c      write(*,*)'EJB Thetas = ',theta_ejb,phi_ejb
c
c end EJB calculation
c      


      ftheta=acos(invect(3)/sqrt(invect(1)**2+invect(3)**2))
      fphi=acos(invect(3)/sqrt(invect(2)**2+invect(3)**2))
      fpsi=acos(sqrt(invect(2)**2+invect(3)**2)/sqrt(invect(1)**2
     &     +invect(2)**2+invect(3)**2))
c
c      write(*,*)'ftheta, fphi, fpsi',
c     &    ftheta*57.296,fphi*57.296,fpsi*57.296
      scvect2(1)=scvect(1)*cos(fpsi)-sin(fpsi)*(scvect(2)*sin(fphi)
     &     +scvect(3)*cos(fphi))
      scvect2(2)=scvect(2)*cos(fphi)-scvect(3)*sin(fphi)
      scvect2(3)=scvect(1)*sin(fpsi)+cos(fpsi)*(scvect(2)*sin(fphi)
     &     +scvect(3)*cos(fphi))
c
c      write(*,*)'SCATTERED 2: ',scvect2(1),scvect2(2),scvect2(3)
      theta=atan(sqrt(scvect2(1)**2+scvect2(2)**2)/scvect2(3))*57.2957795
      phi=atan(scvect2(1)/scvect2(2))*57.2957795
      if (scvect2(1).lt.0.0.and.scvect2(2).gt.0.0) 
     &  phi=phi+360.00
      if (scvect2(1).lt.0.0.and.scvect2(2).lt.0.0) 
     &  phi=phi+180.00
      if (scvect2(1).gt.0.0.and.scvect2(2).lt.0.0) 
     &  phi=phi+180.00

c      write(*,*)'Theta,phi =',theta,phi
      theta=theta_ejb
      phi=phi_ejb
c
      return
      end