c**********************************************************************
c
      subroutine richstep
c
c   design by W. Przygoda (TUM)
c   last modified on: 01/09/2001  by L. Fabbietti
c**********************************************************************

      implicit none

#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"
#include "richtups.inc"
#include "geamedia.inc"

      integer i, ihit, numhit, k
      parameter (numhit=11)
      real hitrich(numhit)
      real xm(3), xd(3), pm(3), pd(3)
      real effi, fate, rndm, divdif
      logical ishit


c --------------- tracking of photons ---------------------------------------
     
      if (ngphot.gt.0 .and.
     +    (numed.eq.med_c4f10 .or.
     +     numed.eq.med_caf2  .or.
     +     numed.eq.med_ch4   .or.
     +     numed.eq.med_det   .or.
     +     numed.eq.med_sio2a .or.
     +     numed.eq.med_sio2b .or.
     +     numed.eq.med_mgf2a .or.
     +     numed.eq.med_mgf2b .or.
     +     numed.eq.med_n2    .or. 
     +     numed.eq.med_mirror)) then
         call gskpho(0)   ! store all Cerenkov photons
c         do i=1,ngphot   ! loop over generated Cerenkov photons and store only
c                         ! those that are detected
c            fate = rndm(i)
c -----  linear interpolation of detection efficiency in table
c -----  (PHOTEFF,PHOTLEN) using CERNlib function divdif (E105)
c            effi = divdif(PHOTEFF,PHOTLEN,OPTBIN,xphot(7,i),1)
c            if(fate.le.effi) call gskpho(i) ! store this Cerenkov photon with
c                                            ! track nb = parent track nb
c         enddo

      endif
c --------------- store hits ------------------------------------------------

      if (numed.eq.med_det) then    ! hit on the pad plane volume

         if(ipart.ne.50) then
            if (inwvol.eq.1) then   ! pad detector volume entered, reset enlos
               enlos = 0.0
            else
               enlos = enlos + destep
            endif
         endif

         if ((ipart.eq.50 .and. inwvol.eq.2) .or.  ! Cerenkov reached pad plane
     +       (ipart.ne.50 .and. 
     +       (inwvol.eq.1 .or. inwvol.eq.2 .or. istop.ge.1)) ) then
                   ! particle entered, left or stopped in detector

            do i=1,3                 ! transform to pad-plane reference frame
               xm(i) = vect(i)
               pm(i) = vect(i+3)
            enddo
            call gmtod(xm,xd,1)
            call gmtod(pm,pd,2)
      
            hitrich(1) = xd(1)          ! x, y z in pad-plane frame
            hitrich(2) = xd(2)
            hitrich(3) = xd(3)
            hitrich(4) = pd(1)          ! momentum components in detector frame
            hitrich(5) = pd(2)
            hitrich(6) = pd(3)
            hitrich(7) = vect(7)        ! total momentum
            hitrich(8) = float(ipart)   ! particle id of hit
            hitrich(9) = 1.e3*enlos     ! energy loss in MeV
            hitrich(10) = float(itra)   ! track number (= parent for Cerenkovs)
            hitrich(11) = float(inwvol) ! volume flag
c           hitrich(12) = 0             ! number of Photons reflected by mirror

            if(ipart.eq.50) then        ! for Cerenkov photons go to eV
               hitrich(7) = 1.e9*vect(7)
               hitrich(9) = 0.0
            end if

            call gsahit(iset,idet,itra,numbv,hitrich,ihit)
            if(ihit.eq.0) write(6,*) 'RICHSTEP: hit not stored'
         endif
      endif
  
c     call gprint('HITS',0)


c---------------counting photons reflected by mirror---------------
        
         if (evtnow.lt.idevt) then 
           do i=1,nmircount
             numphotall(i) = 0
             leptrackall(i) = 0
           enddo
           nmircount = 0
           evtnow = idevt
         endif
         
         k = 1 
         if ((numed .eq. med_c4f10) .or. (numed .eq. med_N2)) then
           do i=1,nmec  
             if (lmec(i) .eq. 106) then
               if(vect(3).gt.0) then

                 ishit = .false.
                 do k=1, nmircount
                   if (leptrackall(k) .eq. itra) then      
                     numphotall(k) = numphotall(k) + 1 ! variable from
                                                       ! common block richtups 
c                    write(6,*) 'number of photons and track ',
c     +              numphotall(k),leptrackall(k)
                     ishit = .true. 
                   endif
                 enddo               
                 if (.not. ishit) then
                   if (nmircount.lt.maxmirr) then 
                     nmircount = nmircount + 1
                     leptrackall(nmircount) = itra
                     numphotall(nmircount) = 1
                   else
                     write(6,*) 'MAXMIRR exceeded in richstep: ',maxmirr
                   endif
                 endif
               endif 
             endif   

           enddo !loop on mechanism active in the present  step      
         endif
       
c ------------------ store hit on mirror ------------------

      if (numed.eq.med_mirror .or. numed.eq.med_glassmirror) then

        if (((ipart.eq.2 .and. inwvol.eq.2) .or.
     +            (ipart.eq.3 .and. inwvol.eq.2)) .and.
     +            (vect(6).gt.0)) then

      ! lepton crosses the mirror
      
            hitrich(1) = vect(1)        ! x, y z in mirror frame
            hitrich(2) = vect(2)
            hitrich(3) = vect(3)
            hitrich(4) = vect(4)        ! momentum components in detector frame
            hitrich(5) = vect(5)
            hitrich(6) = vect(6)
            hitrich(7) = vect(7)        ! total momentum
            hitrich(8) = float(ipart)   ! particle id of hit
c           write(6,*)'Track Number for leptons',itra
c           write(6,*)'Id of lepton',ipart  
            hitrich(9) = 0              ! energy loss in MeV
            hitrich(10) = float(itra)   ! track number (= parent for Cerenkovs)
            hitrich(11) = float(inwvol) ! volume flag
            call gsahit(iset,idet,itra,numbv,hitrich,ihit)
            if(ihit.eq.0) write(6,*) 'RICHSTEP: mirror hit not stored'
        endif
      endif    ! end condition on mirror hit


      return
      end