***********************************************
*                                             *
        SUBROUTINE MDCDIGI
*                                             *
*                                             *
*  last updated on: 28/11/2008 by I. Koenig   *
***********************************************

      IMPLICIT NONE
#include "drtups.inc"
#include "geant321/gcflag.inc"
#include "user.inc"
#include "kinetups.inc"

      INTEGER*4 NVDIM, NHDIM, NHMAX, NSEC, NMOD, NLAY, NTR, NH, NTL
      PARAMETER (NVDIM=4, NHDIM=16, NHMAX=600, NSEC=6, NMOD=4, NLAY=6+1)
      CHARACTER*4 IUDET
      INTEGER*4 ITRACK(NHMAX), NUMBV(NVDIM,NHMAX), NHITS
      REAL*4 HITS(NHDIM,NHMAX)
      INTEGER*4 NUMVS(NVDIM)
      INTEGER*4 JSEC, JMOD, JLAY, JHIT, JTRK, I, J, K, L
      INTEGER*4 TRK(MAXHIT)
      INTEGER*4 TRKHIT(MAXHIT)

*       extract drift-chamber hits from 'JHITS'
      DO 20 JSEC=1,NSEC   ! loop over sectors  (=6)

         NUMVS(1) = JSEC
         NUMVS(2) = 0
         NUMVS(3) = 0
         NUMVS(4) = 0    ! ==> select all hits on sector JSEC in gfhits(...)

         DO 25 JMOD=1,NMOD  ! loop over drift chambers  (=4)
 
            NTR = 0
            NH = 0

            DO 30 JLAY=1,NLAY   ! loop over wire planes (6x sense + 1 cathode)

               NHITS = 0
               NTL = 0          ! number of hits in a single layer
               IUDET = 'D'//CHAR(JMOD+48)//'S'//CHAR(JLAY+48)    ! ==> select det
               IF(JLAY.EQ.NLAY) IUDET = 'D'//CHAR(JMOD+48)//'C4' ! 4th cathode
               CALL GFHITS('DRIF',IUDET,NVDIM,NHDIM,NHMAX,0,NUMVS,
     &               ITRACK,NUMBV,HITS,NHITS)                    ! retrieve hits from Geant store
         
               IF (NHITS.GT.NHMAX) THEN  ! overflow
                 WRITE(6,*) 'ERROR: overflow in mdcdigi.F (nHits>',
     &                      NHMAX,') in MDC layer',JSEC,JMOD,JLAY
                 RETURN
               ENDIF

               DO 40 JHIT=1,NHITS   ! loop over hits in Geant store

                  JTRK = ITRACK(JHIT)
                  IF(JTRK.LE.0) THEN
                     WRITE(6,*) 'ERROR: ITRA=0 in mdcdigi.F'    ! this should never happen...
                     RETURN
                  ENDIF

                  DO K=1,NTR   ! loop over tracks seen so far (in all layers)
                    IF(TRK(K).EQ.JTRK) THEN
                      TRKHIT(K) = TRKHIT(K)+1   ! increment track hit counter
                      IF(TRKHIT(K).GT.MAXMUL) THEN
                        IF(TRKHIT(K).EQ.MAXMUL+1) THEN  ! print message once only
                           WRITE(6,*) 'Track with too many hits (>',
     &                                MAXMUL, ') in MDC',JSEC,JMOD
                        ENDIF
                        GOTO 40
                      ENDIF
                      GOTO 401
                    ENDIF
                  ENDDO
                  IF(NTR.LT.MAXHIT) THEN
                    NTR = NTR + 1
                    TRK(NTR) = JTRK  ! new track, add to list
                    TRKHIT(NTR) = 1
                  ENDIF

 401              NH = NH + 1

                  NTL = NTL + 1
                  IF(NTL.GT.MAXTRKMDC) THEN  ! cannot handle this many...
                     WRITE(6,*) 'Too many Tracks in MDC layer',
     &                          JSEC,JMOD,JLAY
                     RETURN
                  ENDIF

c                 Here come the goodies that go into the tree
c
                  MDCTRAK(NH) = JTRK              ! track number
                  LAYER(NH) = JLAY                ! wire plane number
                  MDCX(NH) = LUNIT*HITS(7,JHIT)   ! y distance (in mm)
                  MDCY(NH) = LUNIT*HITS(8,JHIT)   ! x distance (in mm)
                  MDCTHETA(NH) = HITS(15,JHIT)    ! angle of incidence
                  MDCPHI(NH) = HITS(16,JHIT)      ! azimuthal angle
                  MDCTOF(NH) = HITS(11,JHIT)*1.e9 ! tof (in ns)
                  MDCMOM(NH) = 1.e3*SQRT(HITS(4,JHIT)**2+
     &                           HITS(5,JHIT)**2+HITS(6,JHIT)**2)
                                                  ! momentum (in MeV/c)

c      if (eventid.eq.154) write(6,4111) MDCTRAK(NH),JSEC,JMOD,JLAY,
c     &   MDCTOF(NH), MDCMOM(NH)
c 4111 format('Trk:',i2,' (s,m,l)=',3i2,'  tof= ',f10.6,' p= ',f10.6)
c---------------------------------------------------------------------------
c                  write(6,*) NH,MDCTRAK(NH),JSEC,JMOD,LAYER(NH),NHITS

 40            CONTINUE                  ! end loop over hits

 30         CONTINUE   ! end loop over wire planes

            NHIT = NH

            IF(IMDC.GE.1) THEN
               CALL FILLMDC(JSEC,JMOD)
            ENDIF

 25      CONTINUE               ! end loop over drift chambers

 20   CONTINUE   ! end loop over sectors

      RETURN
      END