!-----------------------------------------------------------------------
! track points inside the RICH needed for cerenkov angle calculation
! M. Staric, jan-2000
!
!--> treba se izpiliti!!!!! 
!
! SUBROUTINE TRACK_IN_RICH(ITRA)
! SUBROUTINE TO_PBEAM_SHROUD(R,DIR)
!
!-----------------------------------------------------------------------


      SUBROUTINE TRACK_IN_RICH(ITRA)

! fills /RIC_TRIC/ and sets /RIC_TRA/RECO_FLAG

      INTEGER ITRA              ! valid track index

#include "riter/RIC_TRA.fi"
#include "riter/RIC_TRIC.fi"
#include "riter/RIC_PAR.fi"
#include "riter/RIC_VES.fi"
#include "riter/RIC_SUPM.fi"
#include "riter/RIC_SFM.fi"

      REAL PATH(0:1)
      LOGICAL OK
      REAL LAMBDA
      REAL PHI_LIM /250.E-3/, LAM_LIM /180.E-3/ ! to nekam v commom!!!

! set reco flag to false
      ITRACK=0
      RECO_FLAG(ITRA)=0
! momentum cut (1 GeV/c)
      IF(T_MOM(ITRA).LT.1) RETURN
! lambda - phi acceptance cut
      DO I=1,3
         TDIR(I)=T_DIR(I,ITRA)
      ENDDO
      IF(TDIR(3).LT.(.5)) RETURN
      LAMBDA=ASIN(TDIR(2))
      IF(ABS(LAMBDA).GT.LAM_LIM) RETURN
      PHI=ATAN(TDIR(1)/TDIR(3))
      IF(ABS(PHI).GT.PHI_LIM) RETURN
! entrance point (entrance window or plane mirror or proton shroud)
      TR(3)=Z_IN
      CALL PREMICA(T_R(1,ITRA),TDIR,TR,PATHA,3)
      CALL TO_PBEAM_SHROUD(TR,TDIR)
      CALL TO_PLANE_MIRROR(TR,TDIR)
      CALL NEW_POINT(TR,TDIR,1.,TR) ! move point for 1 mm (numeric errors!)
! exit point (hits on nonsegmented spherical mirrors), track length
      DO K=0,1
         CALL SFM_HIT(TR,TDIR,K,PATH(K),OK)
         IF(.NOT.OK) RETURN
         CALL NEW_POINT(TR,TDIR,PATH(K),TR1(1,K))
      ENDDO      
      KJE_T=KJE_HIT(TR1(2,0))
      TLEN=PATH(KJE_T)
      IF(TLEN.LE.0) RETURN
      DO K=0,1
         CALL WHICH_SFM(TR1(1,K),K,ID_SPH(K))
      ENDDO
! mean emmission point
      DO I=1,3
         TR0(I)=(TR(I)+TR1(I,KJE_T))/2.
      ENDDO
! reflected track on spherical mirror
      DO K=0,1
         CALL SFM_REFLECT(TR1(1,K),TDIR,K,TDIR1(1,K),AA)
      ENDDO
! plane mirror hit and reflected track
      DO K=0,1
         CALL PLM_HIT(TR1(1,K),TDIR1(1,K),K,PATHA,OK)
         IF(.NOT.OK) RETURN
         CALL NEW_POINT(TR1(1,K),TDIR1(1,K),PATHA,TR2(1,K))
         CALL WHICH_PLM(TR2(1,K),K,ID_PLM(K))
         CALL PLM_REFLECT(TDIR1(1,K),K,TDIR2(1,K),AA)
      ENDDO
! ring center
      IF(.NOT.SM_VIRTUAL) CALL TO_VIRT_POS
      DO K=0,1
         CALL DETECTOR_HIT(TR1(1,K),TDIR1(1,K),K,PATHA,A,B,ISM,OK)
         IF(.NOT.OK) RETURN
         CALL DET_RHIT(A,B,ISM,TRC(1,K))
      ENDDO
! misc
      SIG_T2(1)=T_SIG(1,ITRA)**2
      SIG_T2(2)=T_SIG(2,ITRA)**2
      RMAX2=(THC_LIM*(R_MI(0)+R_MI(1))/4.)**2
! set reco flag to true
      ITRACK=ITRA
      RECO_FLAG(ITRA)=1

      RETURN
      END


      SUBROUTINE TO_PBEAM_SHROUD(R,DIR)

! move point to proton beam shroud, if point is inside shroud

      REAL R(3), DIR(3)
#include "riter/RIC_VES.fi"

      X=R(1)-PBS_RC(1)
      Y=R(2)-PBS_RC(2)
      C=X**2+Y**2-PBS_R**2
      IF(C.GE.0) RETURN
      B=DIR(1)**2+DIR(2)**2
      IF(B.EQ.0) THEN
         R(3)=Z_OUT
         RETURN 
      ENDIF
      A=DIR(1)*X+DIR(2)*Y
      D=A**2-B*C
      IF(D.LT.0) RETURN
      D=SQRT(D)
      PATH=(D-A)/B
      IF(PATH.LE.0) RETURN
      DO I=1,3
         R(I)=R(I)+PATH*DIR(I)
      ENDDO

      RETURN
      END



      SUBROUTINE TO_PLANE_MIRROR(R,DIR)

! move point to plane mirror, if track crosses it

      REAL R(3), DIR(3)
#include "riter/RIC_PLM.fi"

      LOGICAL OK

      K0=KJE_HIT(R(2))
      CALL PLM_HIT(R,DIR,K0,PATH,OK)
      IF(.NOT.OK.OR.PATH.LT.0) RETURN
      Y=R(2)+PATH*DIR(2)
      IF(ABS(Y).GT.PL_YMIN(K0)) THEN
         CALL NEW_POINT(R,DIR,PATH,R)
      ENDIF
      RETURN
      END

