***********************************************************************
c     
      SUBROUTINE RPCDIGI
c     
c     last modified on: 05/12/2008  by I. Koenig
c     
***********************************************************************
      IMPLICIT NONE
#include "geant321/gcflag.inc"
#include "rpctuple.inc"
#include "user.inc"

      INTEGER NHDIM, NHMAX, NVDIM 
      PARAMETER (NHDIM=14,NHMAX=MAXTRKRPC,NVDIM=6)
      INTEGER ITRACK(NHMAX), NUMVS(NVDIM), NUMBV(NVDIM,NHMAX)
      DATA NUMVS /0, 0, 0, 0, 0, 0/ !gas box, col, vicol, tube, inside, gap
      INTEGER NHITS, COL
      REAL HITS(NHDIM,NHMAX)
      CHARACTER CCOL*1, CCELL*1, GGAP*1, NAME*2, IUDET*4
      INTEGER ISECT, I, J, L, N_T, CUR_TR
      REAL THETA, PHI, PXY, RAD2DEG
      
      RAD2DEG = 57.29577

c     loop over all the geant hits
     
      DO ISECT=1, NSECT         

         NTRK = 0                ! NTRK starts from zero for each sector
         
         NUMVS(1) = ISECT        ! 1st index for 'copy volumes'

         DO COL=1, NCOL

            WRITE(CCOL,'(I1)') COL

            DO I=1, NCELLS      

               WRITE(CCELL,'(I1)') I
               IF (I.GT.9) WRITE(CCELL,'(A1)') I+55 !ASCII CODE for A,B,C...

               DO J=1, NGAPS
                  
                  NUMVS(6) = J   ! 2nd index for 'copy volumes'
                  WRITE(NAME,'(A1, A1)') CCOL, CCELL 
                  IUDET = 'EG'//NAME                   ! Name of the active volume

                  CALL GFHITS('RPCG',IUDET,NVDIM,NHDIM,NHMAX,0,NUMVS,
     &                 ITRACK,NUMBV,HITS,NHITS)

                  IF (NHITS.GT.NHMAX) THEN  ! overflow
                    WRITE(6,*) 'ERROR: overflow in rpcdigi.F (nHits>',
     &                         NHMAX,') in RPC gap',ISECT,COL,I,J
                    RETURN
                  ENDIF

c                 loop over all the hits per gap
                  DO L=1,NHITS 
                     IF (NTRK.LT.MAXTRKRPC) THEN

c                       calculate the local angles

                        THETA = RAD2DEG*acos(HITS(13,L))

                        PXY   = sqrt(HITS(11,L)**2 + HITS(12,L)**2)                    

                        IF (PXY.NE.0) THEN
                           IF (HITS(12,L).GT.0.0) THEN
                              PHI = RAD2DEG*acos(HITS(11,L)/PXY) 
                           ELSE 
                              PHI = 360 - RAD2DEG*acos(HITS(11,L)/PXY) 
                           ENDIF
                        ELSE
                           PHI = 0.0
                        ENDIF
                     
c                       fill hits in the common block 
                  
                        NTRK           =   NTRK+1            ! track index [1:MAXTRKRPC]
                        RPCTRK(NTRK)   =   HITS(5,L)         ! geant track number
                        RPCDET(NTRK)   =   4*HITS(7,L) + J-1 ! detector number (begins at zero)
                        RPCE(NTRK)     =   HITS(6,L)*1.0e3   ! energy loss (in MeV)
                        RPCX(NTRK)     =   LUNIT*HITS(1,L)   ! x in mm
                        RPCY(NTRK)     =   LUNIT*HITS(2,L)   ! y in mm
                        RPCZ(NTRK)     =   LUNIT*HITS(3,L)   ! z in mm
                        RPCTHETA(NTRK) =   THETA             ! polar angle in deg
                        RPCPHI(NTRK)   =   PHI               ! azimuthal angle in deg
                        RPCTOF(NTRK)   =   HITS(4,L)*1.0e9   ! tof in ns
                        RPCMOM(NTRK)   =   HITS(8,L)*1.0e3   ! momentum in MeV/c
                        RPCLEN(NTRK)   =   LUNIT*HITS(9,L)   ! total track length in mm
                        RPCLOCLEN(NTRK)=   LUNIT*HITS(10,L)  ! local track length in mm
                     ELSE
                        WRITE(6,*) 'ERROR: too many hits (>',
     &                             MAXTRKRPC,') in RPC sector',ISECT
                        RETURN
                     ENDIF

                  ENDDO         ! end of hit loop

               ENDDO            ! end of gap loop

            ENDDO               ! end of cell loop

         ENDDO                  ! end of column loop


c     Now read the hits of the virtual plane (detector number == -1)

         IUDET    = 'EBOX'      ! Name of the active volume
         NUMVS(6) = 0           ! Set the gap copy number to zero

         
         CALL GFHITS('RPCG',IUDET,NVDIM,NHDIM,NHMAX,0,NUMVS,
     &        ITRACK,NUMBV,HITS,NHITS)
         
         IF (NHITS.GT.NHMAX) THEN  ! overflow
            WRITE(6,*) 'ERROR: overflow in rpcdigi.F (nHits>',
     &                 NHMAX,') in RPC virtual plane',ISECT,COL,I,J
            RETURN
         ENDIF

c     loop over all the hits per gap
         DO L=1,NHITS 
            
            IF (NTRK.LT.MAXTRKRPC) THEN
c     calculate the local angles

               THETA = RAD2DEG*acos(HITS(13,L))

               PXY   = sqrt(HITS(11,L)**2 + HITS(12,L)**2)                    
            
               IF (PXY.NE.0) THEN
                  IF (HITS(12,L).GT.0.0) THEN
                     PHI = RAD2DEG*acos(HITS(11,L)/PXY) 
                  ELSE 
                     PHI = 360 - RAD2DEG*acos(HITS(11,L)/PXY) 
                  ENDIF
               ELSE
                  PHI = 0.0
               ENDIF
            
c     fill hits in the common block 
            
               NTRK           =   NTRK+1           ! track index [1:MAXTRKRPC]
               RPCTRK(NTRK)   =   HITS(5,L)        ! geant track number
               RPCDET(NTRK)   =   -ISECT           ! detector number for EBOX
               RPCE(NTRK)     =   HITS(6,L)*1.0e3  ! energy loss (in MeV)
               RPCX(NTRK)     =   LUNIT*HITS(1,L)  ! x in mm
               RPCY(NTRK)     =   LUNIT*HITS(2,L)  ! y in mm
               RPCZ(NTRK)     =   LUNIT*HITS(3,L)  ! z in mm
               RPCTHETA(NTRK) =   THETA            ! polar angle in deg
               RPCPHI(NTRK)   =   PHI              ! azimuthal angle in deg
               RPCTOF(NTRK)   =   HITS(4,L)*1.0e9  ! tof in ns
               RPCMOM(NTRK)   =   HITS(8,L)*1.0e3  ! momentum in MeV/c
               RPCLEN(NTRK)   =   LUNIT*HITS(9,L)  ! total track length in mm
               RPCLOCLEN(NTRK)=   LUNIT*HITS(10,L) ! local track length in mm
            ELSE
               WRITE(6,*) 'ERROR: too many hits (>',
     &                    MAXTRKRPC,') in RPC virtual plane',ISECT
               RETURN
            ENDIF

         ENDDO

         CALL FILLRPC(ISECT)    ! fill hgeantrpc

      ENDDO                     ! end of sector loop


      RETURN
      END