!-----------------------------------------------------------------------
! cerenkov angle calculation for track-photon pair
! M. Staric, jan-2000, jan-2002
!
!-----------------------------------------------------------------------

      SUBROUTINE CALC_THCER(RF,SIG_F,TH_CER,FI_CER,SIG,IMIR)

      IMPLICIT NONE
      REAL RF(3), SIG_F
      REAL TH_CER, FI_CER
      REAL SIG
      INTEGER IMIR

#include "riter/RIC_PAR.fi"
#include "riter/RIC_TRIC.fi"
#include "riter/RIC_SFM.fi"

      LOGICAL IN_NARROW_WIN, IN_WIDE_WIN, IN_NARROW_WIN2, IN_WIDE_WIN2
      LOGICAL IS_PBS_ONWAY, IS_EBS_ONWAY
      INTEGER KJE_HIT
      LOGICAL OK
      REAL DIRF(3), DIRCH(3), R1(3), DIRF1(3)
      INTEGER K0,ID,KK,ISFM,IPLM
      REAL X,Y,Z,RADIJ2,AA,PATH,PATH1,CFI2

      SIG=0.
! which detector
      K0=KJE_HIT(RF(2))
! coarse limitation
      X=RF(1)-TRC(1,K0)
      Y=RF(2)-TRC(2,K0)
      Z=RF(3)-TRC(3,K0)
      RADIJ2=X**2+Y**2+Z**2
      IF(RADIJ2.GT.RMAX2) RETURN
! approximate cerenkov angle**2
      RADIJ2=RADIJ2/(R_MI(K0)/2.)**2
      IF(.NOT.IN_WIDE_WIN2(RADIJ2)) RETURN
! accurate cerenkov angle
      ID=ID_SPH(K0)
      IF(ID.GE.0) THEN
         CALL CERENKOV_KOT(TR0,TDIR,RF,SF_RC(1,ID),SF_RMI(ID),
     &        DIRF,DIRCH,OK)         
      ELSE
         CALL CERENKOV_KOT(TR0,TDIR,RF,RCMI(1,K0),R_MI(K0),
     &        DIRF,DIRCH,OK)         
      ENDIF
      IF(.NOT.OK) RETURN
! was photon emitted by this track (emit it from entrance to RICH)
      CALL SFM_HIT(TR,DIRF,K0,PATH,OK)
      IF(.NOT.OK.OR.PATH.LT.0) RETURN
      CALL NEW_POINT(TR,DIRF,PATH,R1)
      IF(K0.EQ.0) THEN
         OK=KJE_T.EQ.0.OR.R1(2).GE.0
      ELSE
         OK=KJE_T.EQ.1.OR.R1(2).LE.0
      ENDIF
      IF(.NOT.OK) RETURN
! screening of beam shroud (too restrictive: some photons are lost!!)
! before reflection
      IF(IS_PBS_ONWAY(TR,DIRF,PATH1)) THEN
         IF(PATH1.LT.PATH) RETURN
      ENDIF
! after reflection
      CALL SFM_REFLECT(R1,DIRF,K0,DIRF1,AA)
      IF(IS_PBS_ONWAY(R1,DIRF1,PATH)) RETURN
! cerenkov theta, phi
      TH_CER=SQRT(DIRCH(1)**2+DIRCH(2)**2)
      IF(TH_CER.GT.0) THEN
         FI_CER=ATAN2(DIRCH(2),DIRCH(1))
      ELSE
         FI_CER=0.
      ENDIF
! corrections for misaligned mirrors, if enabled
      IF(LOC_AL) THEN
         CALL WHICH_MIRRORS(TR0,DIRF,K0,ISFM,IPLM,OK)
         IF(OK) CALL CORRECT_THC(ISFM,IPLM,TH_CER,FI_CER)
      ENDIF
! set angular sigma and sign flag
      CFI2=COS(FI_CER)**2
      SIG=SQRT(SIG_F**2+SIG_T2(1)*CFI2+SIG_T2(2)*(1.-CFI2))
      IF(.NOT.IN_NARROW_WIN(TH_CER)) SIG=-SIG
! mirror pairs - for alignment
      IMIR=0
      CALL WHICH_MIRRORS(TR,DIRF,K0,ISFM,IPLM,OK)
      IF(OK) THEN
         IMIR=IOR(ISHFT(IMIR,8),ISFM+1)
         IMIR=IOR(ISHFT(IMIR,8),IPLM+1)
      ENDIF
      CALL WHICH_MIRRORS(TR1(1,K0),DIRF,K0,ISFM,IPLM,OK)
      IF(OK) THEN
         IMIR=IOR(ISHFT(IMIR,8),ISFM+1)
         IMIR=IOR(ISHFT(IMIR,8),IPLM+1)
      ENDIF
  
      RETURN
      END



      SUBROUTINE WHICH_MIRRORS(R,DIR,K0,ISFM,IPLM,OK)

      IMPLICIT NONE
      REAL R(3), DIR(3)         ! reconstructed photon [in]
      INTEGER K0                ! detected on: 0 upper, 1 lower detector [in]
      INTEGER ISFM, IPLM        ! spherical, planar mirror index [out]
      LOGICAL OK                ! status

#include "riter/RIC_TRIC.fi"

      INTEGER KJE_HIT
      REAL PATH, R1(3), DIR1(3), R2(3), AA
      INTEGER KK

      CALL SFM_HIT(R,DIR,K0,PATH,OK)
      IF(.NOT.OK) RETURN
      CALL NEW_POINT(R,DIR,PATH,R1)
      KK=KJE_HIT(R1(2))
      IF(KK.EQ.K0) THEN
         CALL WHICH_SFM(R1,K0,ISFM)
      ELSE
         ISFM=ID_SPH(K0)
      ENDIF
      CALL SFM_REFLECT(R1,DIR,K0,DIR1,AA)
      CALL PLM_HIT(R1,DIR1,K0,PATH,OK)
      IF(.NOT.OK) RETURN
      CALL NEW_POINT(R1,DIR1,PATH,R2)
      CALL WHICH_PLM(R2,K0,IPLM)

      RETURN
      END



      SUBROUTINE CERENKOV_KOT(R0,DIR,RH,RC,R_MI,DIRF,DIRCH,OK)

! calculation of cerenkov angle from photon hit on detector and part. track
! resi_sin: linear approximation

      REAL R0(3)                ! point of photon emission (mid point on track)
      REAL DIR(3)               ! particle direction (unit vector)
      REAL RH(3)                ! photon hit on detector
      REAL RC(3), R_MI          ! center of curvature of mirror

      REAL DIRF(3)              ! result (photon direction vector in HERA-B) if
      REAL DIRCH(3)             ! result (photon direction vector in DIR) if
      LOGICAL OK                ! OK = .true.

      REAL R1(3), R2(3), RR(3), R1N, R2N, RRN
      REAL A, B, ALFA, TH, FI
      REAL C1, C2, S_ALFA

      DO I=1,3
         R1(I)=R0(I)-RC(I)
         R2(I)=RH(I)-RC(I)
      ENDDO
      R1N=SQRT(R1(1)**2+R1(2)**2+R1(3)**2)	
      R2N=SQRT(R2(1)**2+R2(2)**2+R2(3)**2)	
      A=R1N/R_MI
      B=R2N/R1N
      ALFA=(R1(1)*R2(1)+R1(2)*R2(2)+R1(3)*R2(3))/R1N/R2N
      ALFA=MAX(MIN(ALFA,1.),-1.)
      ALFA=ACOS(ALFA)
      CALL RESI_SIN(A,B,ALFA,TH,FI,OK)
      IF(.NOT.OK) RETURN
      S_ALFA=SIN(ALFA)
      C1=R_MI/R1N*SIN(ALFA-FI)/S_ALFA
      C2=R_MI/R2N*SIN(FI)/S_ALFA
      DO I=1,3
         RR(I)=(C1-1.)*R1(I)+C2*R2(I)
      ENDDO
      RRN=SQRT(RR(1)**2+RR(2)**2+RR(3)**2)	
      DO I=1,3
         DIRF(I)=RR(I)/RRN
      ENDDO
      CALL TRANSFORM_FROM_FIXED(DIR,DIRF,DIRCH)
      
      RETURN
      END


      SUBROUTINE RESI_SIN(A,B,C,X,Y,OK)
      
! find solution of equations:	sin(x) = a*sin(x+y)
!				sin(x+y)=b*sin(c+x-y)
! approximative solution (second order)

      REAL A,B,C                ! coeficients
      REAL X,Y                  ! result
      LOGICAL OK                ! if all OK
      
      ABB=1.-2*A*B+B
      IF(ABB.EQ.0) THEN
         OK=.FALSE.
         RETURN
      ENDIF

      X=A*B*C/ABB
      Y=(1.-A)*B*C/ABB

      SX=SIN(X)
      CX=COS(X)
      SXY=SIN(X+Y)
      CXY=COS(X+Y)
      SCXY=SIN(C+X-Y)
      CCXY=COS(C+X-Y)
      
      A11=CX-A*CXY
      A12=-A*CXY
      B1=A*SXY-SX
      A21=CXY-B*CCXY
      A22=CXY+B*CCXY
      B2=B*SCXY-SXY
         
      D=A11*A22-A12*A21
      IF(D.EQ.0) THEN
         OK=.FALSE.
         RETURN
      ENDIF
      
      X=X+(B1*A22-B2*A12)/D
      Y=Y+(B2*A11-B1*A21)/D

      OK=.TRUE.
      RETURN
      END
	  	

