***********************************************************************

      SUBROUTINE MDCSTEP
c
c  store hit info for MDC set DRIF as:
c
c  HITS(1)  = x of current particle in cave system
c  HITS(2)  = y               "
c  HITS(3)  = z               "
c  HITS(4)  = px of current particle in cave system
c  HITS(5)  = py               "
c  HITS(6)  = pz               "
c  HITS(7)  = x in MDC system  (TRD1 shape used for drift volumes, ie x changes
c  HITS(8)  = y in MDC system    along z  --> y and z are permuted) 
c  HITS(9)  = kinetic energy of current particle (in GeV)
c  HITS(10) = current particle ID
c  HITS(11) = TOF of particle (in sec)
c  HITS(12) = distance along wire (= xprime)
c  HITS(13) = perpendicular distance to central wire in plane (= 0-wire)
c  HITS(14) = volume copy number (of level=3)
c  HITS(15) = angle of incidence on plane (0-180 deg)
c  HITS(16) = azimuthal angle (0-360 deg)
c             ->needed to compute the shortest distance to nearest sense wire
c
c  HITS(12),HITS(13), HITS(16) are not used 
c
c  last modified on: 29/11/06  by R.Holzmann
***********************************************************************
      IMPLICIT NONE
#include "drtups.inc"
#include "mdcdconst.inc"
#include "geant321/gcflag.inc"
#include "geant321/gckine.inc"
#include "geant321/gcking.inc"
#include "geant321/gclist.inc"
#include "geant321/gcsets.inc"
#include "geant321/gctmed.inc"
#include "geant321/gctrak.inc"
#include "geant321/gcvolu.inc"
#ifdef CALC_REF_MATRIX
#include "guref.inc"
#endif
#include "kinetups.inc"
      INTEGER NH
      PARAMETER (NH=16)
      REAL HITS(NH)

      CHARACTER*4 CHVOL(10), CHSET
      EQUIVALENCE (IHSET,CHSET)
      EQUIVALENCE (NAMES,CHVOL)

      REAL XM(3), XD(3), PM(3), PD(3), PDZ
      REAL DEG, ANG, FPZ, PXY
      INTEGER I, J, IHIT, ID1, ID2, IP3
      LOGICAL LVOL

      INTEGER ISTEP
      SAVE ISTEP

      IF((INWVOL.EQ.1 .OR. INWVOL.EQ.2)
     +          .AND. CHSET.EQ.'DRIF' .AND. CHARGE.NE.0.) THEN
                      ! charged track has entered or left MDC drift volume

         IF(VECT(7)<0.0002) RETURN  ! momentum is below 0.2 MeV, abort

         ID1 = IDTYPE/10           ! identify wire plane via IDTYPE... (gcsets)
         ID2 = IDTYPE - 10*ID1
         ANG = PHI(ID1,ID2)        ! ...and select appropriate wire angle

         IF(MDCSHAPE(1:4).NE.'TRD1')
     +           WRITE(6,*) ' MDCSTEP: Wrong shape: ', mdcshape

         DO I = 1,3
            HITS(I) = VECT(I)      ! position of hit and particle momentum
            HITS(I+3) = VECT(I+3)*VECT(7)
            XM(I) = VECT(I)
            PM(I) = VECT(I+3)
         ENDDO

         CALL GMTOD(XM,XD,1)       ! transform coordinates into MDC system
         CALL GMTOD(PM,PD,2)       ! transform momenta into MDC system

         PDZ =  PD(2)        ! pz is along GEANT y axis of drift volumes
         FPZ = PDZ/ABS(PDZ)  ! particle goes forwards or backwards ?

         IF((INWVOL.EQ.1 .AND. FPZ.LT.0) .OR.      ! entered backwards
     +      (INWVOL.EQ.2 .AND. FPZ.GT.0)) RETURN   ! or left forwards

         IF (INCWIRES.EQ.0) THEN        ! geometry without wires
            IF (NSTEP-ISTEP.EQ.1) THEN  ! next hit follows too close to previous
               ISTEP = NSTEP
               RETURN
            ENDIF
         ELSE
            IF (NSTEP-ISTEP.EQ.1) THEN  ! no wire hit or too close to previous
               ISTEP = NSTEP
               RETURN
            ELSE
               IF ((FPZ.GT.0 .AND. INFROM.GT.0) .OR.
     +             ( FPZ.LT.0 .AND. INGOTO.GT.0)) THEN
                  RETURN
               ENDIF
            ENDIF
         ENDIF
 
         ISTEP = NSTEP

         HITS(7) = -XD(1)+MDCDX(ID2,ID1)
         HITS(8) = XD(3)+MDCDY(ID2,ID1)
         HITS(9) = GEKIN
         HITS(10) = FLOAT(IPART)
         HITS(11) = TOFG
         HITS(14) = FLOAT(NUMBER(3))
         HITS(15) = 57.29577*ACOS(PD(2))   ! angle of incidence  (with y <-> z)
         PXY = SQRT(PD(1)**2+PD(3)**2)     ! with y <-> z
         IF(PXY.NE.0.0) THEN
            HITS(16) = 57.29577*ACOS(-PD(1)/PXY)   ! x axis flipped
         ELSE
            HITS(16) = 0.0
         ENDIF
         IF(PD(3).LT.0.0) HITS(16) = 360.0 - HITS(16)

         CALL GSAHIT(ISET,IDET,ITRA,NUMBV,HITS,IHIT)
      ENDIF

      RETURN
      END