!----------------------------------------------------------------------
!
! Routines for ray tracing of photons inside the HERA-B RICH
!
! author: Marko Staric, IJS, Ljubljana
! version: 1.0, 16-dec-1999
! version: 1.1, 29-jan-2002, in (SINGLE_)SFM_HIT, (SINGLE_)PLM_HIT PATH not tested
!
!----------------------------------------------------------------------
!*** util
c      INTEGER FUNCTION KJE_HIT(Y)
c      SUBROUTINE NEW_POINT(R,DIR,PATH,RH)
!*** beam shrouds ***
c      LOGICAL FUNCTION IS_PBS_ONWAY(R,DIR,PATH)
c      LOGICAL FUNCTION IS_EBS_ONWAY(R,DIR,PATH)
!*** spherical mirror ***
c      SUBROUTINE SFM_HIT(R,DIR,K0,PATH,OK)
c      SUBROUTINE SINGLE_SFM_HIT(R,DIR,ID,PATH,OK)
c      SUBROUTINE WHICH_SFM(R,K0,ID)
c      LOGICAL FUNCTION INSIDE_SFM(R,ID)
c      SUBROUTINE SFM_REFLECT(R,DIR,K0,DIR_REFL,A)
c      SUBROUTINE SINGLE_SFM_REFLECT(R,DIR,ID,DIR_REFL,A)
!*** plane mirror ***
c      SUBROUTINE PLM_HIT(R,DIR,K0,PATH,OK)
c      SUBROUTINE SINGLE_PLM_HIT(R,DIR,ID,PATH,OK)
c      SUBROUTINE WHICH_PLM(R,K0,ID)
c      LOGICAL FUNCTION INSIDE_PLM(R,ID)
c      SUBROUTINE PLM_REFLECT(DIR,K0,DIR_REFL,A)
c      SUBROUTINE SINGLE_PLM_REFLECT(DIR,ID,DIR_REFL,A)
!*** detector entry window ***
c      SUBROUTINE DETWIN_HIT(R,DIR,K0,PATH,OK)
!*** detector ***
c      SUBROUTINE DETECTOR_HIT(R,DIR,K0,PATH,ALFA,BETA,ISM,OK)
c      SUBROUTINE SUPMOD_HIT(R,DIR,ISM,PATH,ALFA,BETA,OK)
c      SUBROUTINE DET_CHANNEL(ALFA,BETA,ISPM,ICH)
c      SUBROUTINE DET_RHIT(ALFA,BETA,ISM,RH)
!----------------------------------------------------------------------


      INTEGER FUNCTION KJE_HIT(Y)

      REAL Y                    ! [in]

      IF(Y.GT.0) THEN
         KJE_HIT=0
      ELSE
         KJE_HIT=1
      ENDIF
      RETURN
      END


      SUBROUTINE NEW_POINT(R,DIR,PATH,RH)

! moves the point R() along DIR() by PATH

      REAL R(3), DIR(3), PATH   ! [in]
      REAL RH(3)                ! [out]

      DO I=1,3
         RH(I)=R(I)+PATH*DIR(I)
      ENDDO

      RETURN
      END


!*** beam shrouds ***

      LOGICAL FUNCTION IS_PBS_ONWAY(R,DIR,PATH)

! returns .true., if ProtonBeamShroud is on the way (in forward dir)

#include "riter/RIC_VES.fi"
      REAL R(3), DIR(3)         ! photon [in]
      REAL PATH                 ! length from R() to hit, if true and >0 [out]

      X=R(1)-PBS_RC(1)
      Y=R(2)-PBS_RC(2)
      Z=R(3)-PBS_RC(3)
      C=X**2+Y**2-PBS_R**2
      IF(C.LT.0) THEN           ! photon inside the tube
         PATH=-1.
         IS_PBS_ONWAY=.TRUE.
      ELSE IF(C.EQ.0) THEN      ! photon on the tube wall
         PATH=0.
         B=DIR(1)*X+DIR(2)*Y
         IS_PBS_ONWAY=B.LT.0
      ELSE                      ! photon outside the tube
         IS_PBS_ONWAY=.FALSE.
         A=DIR(1)**2+DIR(2)**2
         IF(A.EQ.0) RETURN
         B=DIR(1)*X+DIR(2)*Y
         D=B**2-A*C
         IF(D.LT.0) RETURN
         D=SQRT(D)
         PATH=-(D+B)/A
         IF(PATH.LE.0) RETURN
         Z1=Z+PATH*DIR(3)
         IS_PBS_ONWAY=Z1.GT.0.AND.Z1.LT.PBS_LEN
         IF(.NOT.IS_PBS_ONWAY) THEN
            PATH=(D-B)/A
            Z1=Z+PATH*DIR(3)
            IS_PBS_ONWAY=Z1.GT.0.AND.Z1.LT.PBS_LEN
         ENDIF
      ENDIF
      
      RETURN
      END
   


      LOGICAL FUNCTION IS_EBS_ONWAY(R,DIR,PATH)

! returns .true., if ElectronBeamShroud is on the way

#include "riter/RIC_VES.fi"
      REAL R(3), DIR(3)         ! photon [in]
      REAL PATH                 ! length from R() to hit, if true and >0 [out]


      X=R(1)-EBS_RC(1)
      Y=R(2)-EBS_RC(2)
      Z=R(3)-EBS_RC(3)
      C=X**2+Y**2-EBS_R**2
      IF(C.LT.0) THEN           ! photon inside the tube
         PATH=-1.
         IS_EBS_ONWAY=.TRUE.
      ELSE IF(C.EQ.0) THEN      ! photon on the tube wall
         PATH=0.
         B=DIR(1)*X+DIR(2)*Y
         IS_EBS_ONWAY=B.LT.0
      ELSE                      ! photon outside the tube
         IS_EBS_ONWAY=.FALSE.
         A=DIR(1)**2+DIR(2)**2
         IF(A.EQ.0) RETURN
         B=DIR(1)*X+DIR(2)*Y
         D=B**2-A*C
         IF(D.LT.0) RETURN
         D=SQRT(D)
         PATH=-(D+B)/A
         IF(PATH.LE.0) RETURN
         Z1=Z+PATH*DIR(3)
         IS_EBS_ONWAY=Z1.GT.0.AND.Z1.LT.EBS_LEN
         IF(.NOT.IS_EBS_ONWAY) THEN
            PATH=(D-B)/A
            Z1=Z+PATH*DIR(3)
            IS_EBS_ONWAY=Z1.GT.0.AND.Z1.LT.EBS_LEN
         ENDIF
      ENDIF
      
      RETURN
      END


   
!*** spherical mirror ***

      SUBROUTINE SFM_HIT(R,DIR,K0,PATH,OK)

! returns length from R to photon hit on the sphere defined by nonsegmented 
! spherical mirror K0 (concave crossection)

#include "riter/RIC_SFM.fi"

      REAL R(3), DIR(3)         ! photon [in]
      INTEGER K0                ! 0 upper, 1 lower mirror [in]
      REAL PATH                 ! length from R to hit [out]
      LOGICAL OK                ! true, if solution exists  [out] 

      REAL R0(3)

      DO I=1,3
         R0(I)=R(I)-RCMI(I,K0)
      ENDDO
      RDIR=R0(1)*DIR(1)+R0(2)*DIR(2)+R0(3)*DIR(3)
      RR=R0(1)**2+R0(2)**2+R0(3)**2
      D=RDIR**2+R_MI(K0)**2-RR
      OK=D.GE.0
      IF(.NOT.OK) RETURN
      D=SQRT(D)
      PATH=D-RDIR

      RETURN
      END
  


      SUBROUTINE SINGLE_SFM_HIT(R,DIR,ID,PATH,OK)

! returns length from R to photon hit on the sphere defined by single
! spherical mirror ID (concave crossection)

#include "riter/RIC_SFM.fi"

      REAL R(3), DIR(3)         ! photon [in]
      INTEGER ID                ! mirror index [in]
      REAL PATH                 ! length from R to hit [out]
      LOGICAL OK                ! true, if solution exists [out] 

      REAL R0(3)

      OK=ID.GE.0.AND.ID.LT.80
      IF(.NOT.OK) RETURN
      DO I=1,3
         R0(I)=R(I)-SF_RC(I,ID)
      ENDDO
      RDIR=R0(1)*DIR(1)+R0(2)*DIR(2)+R0(3)*DIR(3)
      RR=R0(1)**2+R0(2)**2+R0(3)**2
      D=RDIR**2+SF_RMI(ID)**2-RR
      OK=D.GE.0
      IF(.NOT.OK) RETURN
      D=SQRT(D)
      PATH=D-RDIR

      RETURN
      END
  


      SUBROUTINE WHICH_SFM(R,K0,ID)

! finds ID of the spherical mirror where the point R is
! R should be on the sphere defined by nonsegmented mirror K0

#include "riter/RIC_SFM.fi"

      REAL R(3)                 ! point on the spherical mirror [in]
      INTEGER K0                ! 0=upper, 1=lower mirror [in]
      INTEGER ID                ! ID of single mirror (0:79, -1 if out) [out]

! check, if point is inside nonsegmented mirror
      ID=-1
      IF(R(2).GT.0) THEN
         IF(K0.NE.0) RETURN
      ELSE
         IF(K0.NE.1) RETURN
      ENDIF
      IF(R(3).LE.0) RETURN
      TX=R(1)/R(3)
      IF(ABS(TX).GT.SF_TX.OR.ABS(R(2)).GT.SF_Y) RETURN

! find single mirror
      IF(R(1).GT.0) THEN
         K2=49-10*K0
         K1=K2-5
         KST=10*(1-2*K0)
         D_MIN=1.E30
         DO I=1,4
            D_MI=1.E30
            DO K=K1,K2
               D=(R(1)-SF_CEN(1,K))**2+(R(2)-SF_CEN(2,K))**2+
     &              (R(3)-SF_CEN(3,K))**2
               IF(D.GT.D_MI) GOTO 1
               D_MI=D
            ENDDO
 1          CONTINUE
            IF(D_MI.GT.D_MIN) GOTO 2
            D_MIN=D_MI
            ID=K-1
            K2=K2+KST
            K1=K2-5+MOD(I,2)
         ENDDO
 2       CONTINUE
      ELSE
         K2=40-10*K0
         K1=K2+4
         KST=10*(1-2*K0)
         D_MIN=1.E30
         DO I=1,4
            D_MI=1.E30
            DO K=K1,K2,-1
               D=(R(1)-SF_CEN(1,K))**2+(R(2)-SF_CEN(2,K))**2+
     &              (R(3)-SF_CEN(3,K))**2
               IF(D.GT.D_MI) GOTO 3
               D_MI=D
            ENDDO
 3          CONTINUE
            IF(D_MI.GT.D_MIN) GOTO 4
            D_MIN=D_MI
            ID=K+1
            K2=K2+KST
            K1=K2+4+MOD(I,2)
         ENDDO
 4       CONTINUE
      ENDIF

      RETURN
      END



      LOGICAL FUNCTION INSIDE_SFM(R,ID)

! returns true, if R is inside single spherical mirror ID
! R should be on the sphere defined by single mirror

#include "riter/RIC_SFM.fi"

      REAL R(3)                 ! point on the single mirror [in]
      INTEGER ID                ! ID of single mirror [in]

      INSIDE_SFM=.FALSE.
      IF(ID.LT.0.OR.ID.GT.79) RETURN
      DO I=1,SF_NVTC(ID)
         II=I-1
         IF(II.EQ.0) II=SF_NVTC(ID)
         X1=SF_VTCS(1,I,ID)-SF_VTCS(1,II,ID)
         Y1=SF_VTCS(2,I,ID)-SF_VTCS(2,II,ID)
         X2=R(1)-SF_VTCS(1,II,ID)
         Y2=R(2)-SF_VTCS(2,II,ID)
         Z=X1*Y2-Y1*X2
         INSIDE_SFM=Z.GT.0
         IF(.NOT.INSIDE_SFM) RETURN
      ENDDO
! hole for proton beam shroud
      IF(ID.EQ.34.OR.ID.EQ.44) THEN
         INSIDE_SFM=R(1).LT.0.OR.ABS(R(2)).GT.SF_PHOLE_R
         IF(.NOT.INSIDE_SFM) RETURN
         INSIDE_SFM=R(1)**2+R(2)**2.GT.SF_PHOLE_R**2
         IF(.NOT.INSIDE_SFM) RETURN
      ENDIF
! hole for electron beam shroud
      IF(R(2).LT.0) THEN
         RR=(R(1)-SF_EHOLE_X)**2+(R(2)-SF_EHOLE_Y)**2
         INSIDE_SFM=RR.GT.SF_EHOLE_R**2
      ENDIF

      RETURN
      END
 


      SUBROUTINE SFM_REFLECT(R,DIR,K0,DIR_REFL,A)

! reflection on nonsegmented spherical mirror K0
! R() must be photon hit on the mirror

#include "riter/RIC_SFM.fi"

      REAL R(3), DIR(3)         ! photon (point on the mirror) [in]
      INTEGER K0                ! 0 upper, 1 lower mirror [in]
      REAL DIR_REFL(3)          ! reflected ray [out]
      REAL A                    ! if A<0, hit into surface with refl. coating  

      REAL NORM(3)
      
      DO I=1,3
         NORM(I)=(RCMI(I,K0)-R(I))/R_MI(K0)
      ENDDO
      A=DIR(1)*NORM(1)+DIR(2)*NORM(2)+DIR(3)*NORM(3)
      DO I=1,3
         DIR_REFL(I)=DIR(I)-2.*A*NORM(I)
      ENDDO

      RETURN
      END


      SUBROUTINE SINGLE_SFM_REFLECT(R,DIR,ID,DIR_REFL,A)

! reflection on single spherical mirror ID
! R() must be photon hit on the mirror

#include "riter/RIC_SFM.fi"

      REAL R(3), DIR(3)         ! photon (point on the mirror) [in]
      INTEGER ID                ! valid mirror index (0:79) [in]
      REAL DIR_REFL(3)          ! reflected ray [out]
      REAL A                    ! A<0, hit into surface with refl.coating [out]

      REAL NORM(3)
      
      DO I=1,3
         NORM(I)=(SF_RC(I,ID)-R(I))/SF_RMI(ID)
      ENDDO
      A=DIR(1)*NORM(1)+DIR(2)*NORM(2)+DIR(3)*NORM(3)
      DO I=1,3
         DIR_REFL(I)=DIR(I)-2.*A*NORM(I)
      ENDDO

      RETURN
      END


!*** plane mirror ***

      SUBROUTINE PLM_HIT(R,DIR,K0,PATH,OK)

! returns the length from R() to photon hit on the plane defined by 
! nonsegmented plane mirror K0 (crossection in forward direction)

#include "riter/RIC_PLM.fi"

      REAL R(3), DIR(3)         ! photon [in]
      INTEGER K0                ! 0 upper, 1 lower mirror [in]
      REAL PATH                 ! length from R to hit [out]
      LOGICAL OK                ! true, if solution exists [out]

      REAL R0(3)

      DO I=1,3
         R0(I)=R(I)-PL_R0(I,K0)
      ENDDO
      CALL SKAL_PROD(R0,PL_VN(1,K0),S1)
      CALL SKAL_PROD(DIR,PL_VN(1,K0),S2)
      OK=S2.NE.0
      IF(.NOT.OK) RETURN
      PATH=-S1/S2

      RETURN
      END



      SUBROUTINE SINGLE_PLM_HIT(R,DIR,ID,PATH,OK)

! returns the length from R() to photon hit on the plane defined by 
! single plane mirror K0 (crossection in forward direction)

#include "riter/RIC_PLM.fi"

      REAL R(3), DIR(3)         ! photon [in]
      INTEGER ID                ! mirror index [in]
      REAL PATH                 ! length from R to hit [out]
      LOGICAL OK                ! true, if solution exists [out]

      REAL R0(3)

      OK=ID.GE.0.AND.ID.LT.36
      IF(.NOT.OK) RETURN
      DO I=1,3
         R0(I)=R(I)-PL_ORIG(I,ID)
      ENDDO
      CALL SKAL_PROD(R0,PL_NOR(1,ID),S1)
      CALL SKAL_PROD(DIR,PL_NOR(1,ID),S2)
      OK=S2.NE.0
      IF(.NOT.OK) RETURN
      PATH=-S1/S2

      RETURN
      END



      SUBROUTINE WHICH_PLM(R,K0,ID)

! finds ID of plane mirror where the point R is
! R should be on the plane defined by nonsegmented mirror K0

#include "riter/RIC_PLM.fi"

      REAL R(3)                 ! point on the plane mirror [in]
      INTEGER K0                ! 0=upper, 1=lower mirror [in]
      INTEGER ID                ! ID of single mirror (0:35, -1 if out) [out]

      REAL R0(3)

      DO I=1,3
         R0(I)=R(I)-PL_R0(I,K0)
      ENDDO
      CALL SKAL_PROD(R0,PL_VX(1,K0),X)
      X=X/PL_DX(K0)
      CALL SKAL_PROD(R0,PL_VY(1,K0),Y)
      Y=Y/PL_DY(K0)
      IF(X.LT.0.OR.X.GE.1.OR.Y.LT.0.OR.Y.GE.1) THEN
         ID=-1
         RETURN
      ENDIF
      X=X*6.
      IX=IFIX(X)
      Y=Y*3.
      IY=IFIX(Y)
      ID=IX+6*IY+(1-K0)*18

      RETURN
      END
 


      LOGICAL FUNCTION INSIDE_PLM(R,ID)

! returns true, if R is inside single plane mirror ID
! R should be on the plane defined by single mirror

#include "riter/RIC_PLM.fi"

      REAL R(3)                 ! point on the single mirror [in]
      INTEGER ID                ! ID of single mirror [in]
 
      REAL R0(3)

      INSIDE_PLM=.FALSE.
      IF(ID.LT.0.OR.ID.GT.35) RETURN
      DO I=1,3
         R0(I)=R(I)-PL_ORIG(I,ID)
      ENDDO
      CALL SKAL_PROD(R0,PL_X_VEC(1,ID),X)
      X=X/PL_DXVEC(ID)
      CALL SKAL_PROD(R0,PL_Y_VEC(1,ID),Y)
      Y=Y/PL_DYVEC(ID)
      INSIDE_PLM=X.GT.0.AND.X.LT.1.AND.Y.GT.0.AND.Y.LT.1
! hole for electron beam shroud
      IF(R(2).LT.0) THEN
         X=R(1)-PL_HOLE_X
         Y=R(2)-PL_HOLE_Y
         INSIDE_PLM=INSIDE_PLM.AND.X**2+Y**2.GT.PL_HOLE_R**2
      ENDIF

      RETURN
      END
 


      SUBROUTINE PLM_REFLECT(DIR,K0,DIR_REFL,A)

! reflection on nonsegmented plane mirror K0

#include "riter/RIC_PLM.fi"

      REAL DIR(3)               ! photon direction [in]
      INTEGER K0                ! 0 upper, 1 lower mirror [in]
      REAL DIR_REFL(3)          ! reflected ray [out]
      REAL A                    ! A<0, hit into surface with refl.coating [out]

      A=DIR(1)*PL_VN(1,K0)+DIR(2)*PL_VN(2,K0)+DIR(3)*PL_VN(3,K0)
      DO I=1,3
         DIR_REFL(I)=DIR(I)-2.*A*PL_VN(I,K0)
      ENDDO

      RETURN
      END



      SUBROUTINE SINGLE_PLM_REFLECT(DIR,ID,DIR_REFL,A)

! reflection on single plane mirror ID

#include "riter/RIC_PLM.fi"

      REAL DIR(3)               ! photon direction [in]
      INTEGER ID                ! valid mirror index (0:35) [in]
      REAL DIR_REFL(3)          ! reflected ray [out]
      REAL A                    ! A<0, hit into surface with refl.coating [out]

      A=DIR(1)*PL_NOR(1,ID)+DIR(2)*PL_NOR(2,ID)+DIR(3)*PL_NOR(3,ID)
      DO I=1,3
         DIR_REFL(I)=DIR(I)-2.*A*PL_NOR(I,ID)
      ENDDO

      RETURN
      END


!*** detector entry window ***

      SUBROUTINE DETWIN_HIT(R,DIR,K0,PATH,OK)

! returns length from R to photon hit on the detector entry window K0
! (crossection inside the window in forward direction)

#include "riter/RIC_DETW.fi"

      REAL R(3), DIR(3)         ! photon [in]
      INTEGER K0                ! 0=upper, 1=lower [in]
      REAL PATH                 ! length from R to hit, if OK=.true. [out]
      LOGICAL OK                ! if true, PATH is returned [out]

      REAL  P(3)

      SN=S_PROD(DIR,DW_VN(1,K0))
      OK=SN.NE.0
      IF(.NOT.OK) RETURN
      DO I=1,3
         P(I)=R(I)-DW_R0(I,K0)
      ENDDO
      S=S_PROD(P,DW_VN(1,K0))
      PATH=-S/SN
      OK=PATH.GT.0
      IF(.NOT.OK) RETURN
      DO I=1,3
         P(I)=P(I)+PATH*DIR(I)
      ENDDO
      A=S_PROD(P,DW_VX(1,K0))/DW_DX(K0)
      OK=A.GT.0.AND.A.LT.1
      IF(.NOT.OK) RETURN
      A=S_PROD(P,DW_VY(1,K0))/DW_DY(K0)
      OK=A.GT.0.AND.A.LT.1

      RETURN
      END


!*** detector ***

      SUBROUTINE DETECTOR_HIT(R,DIR,K0,PATH,ALFA,BETA,ISM,OK)

! returns supermodule ID and normalised projections of a hit into 
! supermodule x and y axis. For hit inside the active area ALFA and BETA
! have values between 0 and 1. Values are valid, if OK=.true.
 
      REAL R(3), DIR(3)         ! photon [in]
      INTEGER K0                ! 0=upper, 1=lower [in]
      REAL PATH                 ! length from R to hit, if OK=.true. [out]
      REAL ALFA, BETA           ! normalised projections (x,y) [out]
      INTEGER ISM               ! supermodule ID [out]
      LOGICAL OK                ! true, if hit in forward direction [out]

      ISM=3+K0*7
      CALL SUPMOD_HIT(R,DIR,ISM,PATH,ALFA,BETA,OK)
      IF(.NOT.OK) RETURN
      IF(ALFA.LT.0) THEN
         DO WHILE(ALFA.LT.0.AND.ISM.GT.K0*7)
            ISM=ISM-1
            CALL SUPMOD_HIT(R,DIR,ISM,PATH,ALFA,BETA,OK)
            IF(.NOT.OK) RETURN
         ENDDO
      ELSE IF(ALFA.GT.1) THEN
         DO WHILE(ALFA.GT.1.AND.ISM.LT.6+K0*7)
            ISM=ISM+1
            CALL SUPMOD_HIT(R,DIR,ISM,PATH,ALFA,BETA,OK)
            IF(.NOT.OK) RETURN
         ENDDO
      ENDIF

      RETURN
      END



      SUBROUTINE SUPMOD_HIT(R,DIR,ISM,PATH,ALFA,BETA,OK)

#include "riter/RIC_SUPM.fi"

      REAL R(3), DIR(3)         ! photon [in]
      INTEGER ISM               ! valid supermodule ID (0:13) [in]
      REAL PATH                 ! length from R to hit, if OK=.true. [out]
      REAL ALFA, BETA           ! normalised projections (x,y) [out]
      LOGICAL OK                ! if true, ALFA and BETA returned [out]

      REAL  P(3)

      SN=S_PROD(DIR,SM_VN(1,ISM))
      OK=SN.NE.0
      IF(.NOT.OK) RETURN
      DO I=1,3
         P(I)=R(I)-SM_R0(I,ISM)
      ENDDO
      S=S_PROD(P,SM_VN(1,ISM))
      PATH=-S/SN
      OK=PATH.GT.0
      IF(.NOT.OK) RETURN
      DO I=1,3
         P(I)=P(I)+PATH*DIR(I)
      ENDDO
      ALFA=S_PROD(P,SM_VX(1,ISM))/SM_DX(ISM)
      BETA=S_PROD(P,SM_VY(1,ISM))/SM_DY(ISM)

      RETURN
      END



      SUBROUTINE DET_CHANNEL(ALFA,BETA,ISPM,ICH)

! returns detector channel number or -1 (-2), if it doesn't exist

#include "riter/RIC_CHA.fi"

      REAL ALFA, BETA           ! normalised projections into supmod axes [in]
      INTEGER ISPM              ! valid supmod ID (0:13) [in]
      INTEGER ICH               ! channel number, or -1 (-2) [out]

      ICH=-2
      IF(BETA.LT.0.OR.BETA.GE.1) RETURN
      IF(ALFA.LT.0.OR.ALFA.GE.1) RETURN
      IROW=IFIX(BETA*128)
      ICOL=IFIX(ALFA*48)
      ICH=CH_DET(IROW,ICOL,ISPM)
      RETURN
      END


  
      SUBROUTINE DET_RHIT(ALFA,BETA,ISM,RH)

! returns coordinates of a hit

#include "riter/RIC_SUPM.fi"

      REAL ALFA, BETA           ! normalised projections into supmod axes [in]
      INTEGER ISM               ! valid supmod ID (0:13) [in]
      REAL RH(3)                ! hit in HERA-B coordinates [out]

      DO I=1,3
         RH(I)=ALFA*SM_VX(I,ISM)+BETA*SM_VY(I,ISM)+SM_R0(I,ISM)
      ENDDO
      RETURN
      END





