!-----------------------------------------------------------------------
! estimation of signal
! M. Staric, jan-2000, rev. jul-2000, nov-2000, jan-2002
!
! SUBROUTINE SIGNAL_ESTIMATE(ITRA) - this does all
! SUBROUTINE RICH_ACC(THC,ACC)
! SUBROUTINE RAY_TRACE_0(R,DIR,K0,Y,ICH)
! SUBROUTINE RAY_TRACE_1(R,DIR,K0,ICH)
! REAL FUNCTION CHER_ANGLE(P,MASS,REF_IND)                                  
!
!-----------------------------------------------------------------------

      SUBROUTINE SIGNAL_ESTIMATE(ITRA)

! must be called after TRACK_IN_RICH(ITRA)
! fills /RIC_TRA/THC(), /RIC_TRA/STFS()
! resets RECO_FLAG, if there is no expected photons for any mass hyp.,
! or sets it to 100*acc(electron)

      INTEGER ITRA              ! valid track index
 
#include "riter/RIC_PAR.fi"
#include "riter/RIC_TRA.fi"
#include "riter/RIC_TRIC.fi"
      
      integer NMASS
      PARAMETER (NMASS=5)
      REAL P_MASS(NMASS)
      DATA P_MASS /.511E-3,105.66E-3,139.57E-3,493.68E-3,938.27E-3/
      SAVE P_MASS
      DATA DTHC /2.E-3/
      SAVE DTHC

      IF(ITRA.LT.1) RETURN
      IF(ITRA.GT.NTRA) RETURN
! cerenkov angle of hypetheses
      THC(0,ITRA)=0.
      DO I=1,MIN(NHYP-1,NMASS)
         THC(I,ITRA)=CHER_ANGLE(T_MOM(ITRA),P_MASS(I),REF_IND)
      ENDDO
! geometrical acceptance and number of photons
      S=0.
      STFS(0,ITRA)=0.
      ACC=0.
      DO I=1,MIN(NHYP-1,NMASS)
         IF(ABS(THC(I,ITRA)-THC(I-1,ITRA)).GT.DTHC) THEN
            CALL RICH_ACC(THC(I,ITRA),ACC)
         ENDIF
         STFS(I,ITRA)=N_ZERO*TLEN*ETA*ACC*THC(I,ITRA)**2
         S=MAX(S,STFS(I,ITRA))
         IF(I.EQ.1) ACCE=ACC
      ENDDO
      IF(S.EQ.0) THEN
         RECO_FLAG(ITRA)=0
      ELSE
         RECO_FLAG(ITRA)=NINT(100.*ACCE)
      ENDIF
      DO I=NMASS+1,NHYP-1
         THC(I,ITRA)=THC(1,ITRA)
         STFS(I,ITRA)=STFS(1,ITRA)*(I-NMASS+1)
      ENDDO

      RETURN
      END


      SUBROUTINE RICH_ACC(THC,ACC)

! calculates geometrical acceptance for given track using ray_tracing routines

      REAL THC                  ! cerenkov angle [in]
      REAL ACC                  ! geometrical acceptance [out]

#include "riter/RIC_TRIC.fi"
#include "riter/RIC_SUPM.fi"

      PARAMETER (NRAY=50)
      REAL DIRF(3), DIRCH(3)
      REAL DFI, CFI, SFI
      LOGICAL START /.TRUE./
      SAVE PI, DFI, CFI, SFI, START

      IF(START) THEN
         START=.FALSE.
         PI=4.*ATAN(1.)
         DFI=2.*PI/NRAY
         CFI=COS(DFI)
         SFI=SIN(DFI)
      ENDIF
      ACC=0.
      IF(THC.LE.0) RETURN
      DIRCH(1)=SIN(THC)
      DIRCH(2)=0.
      DIRCH(3)=COS(THC)
      IF(.NOT.SM_VIRTUAL) CALL TO_VIRT_POS

      K0=KJE_T
      DO IRAY=1,NRAY
         XX=DIRCH(1)
         YY=DIRCH(2)
         DIRCH(1)=XX*CFI-YY*SFI
         DIRCH(2)=XX*SFI+YY*CFI
         CALL ROTATE_TO_FIXED(TDIR,DIRCH,DIRF)
         CALL RAY_TRACE_0(TR,DIRF,K0,Y,ICH)
         IF(ICH.GE.-2) THEN
            ZN=0.
            IF(K0.EQ.KJE_T) THEN
               IF(ICH.GE.0) ZN=1.
            ELSE
               DD=Y-TR1(2,KJE_T)
               IF(DD.NE.0) THEN
                  PP=ABS(Y/DD)
               ELSE
                  PP=0.
               ENDIF
               IF(ICH.GE.0) ZN=PP
               CALL RAY_TRACE_1(TR1(1,KJE_T),DIRF,KJE_T,ICH)
               IF(ICH.GE.0) ZN=ZN+1.-PP
            ENDIF
            ACC=ACC+ZN
         ENDIF
      ENDDO
      ACC=ACC/NRAY

      RETURN
      END



      SUBROUTINE RAY_TRACE_0(R,DIR,K0,Y,ICH)

! from inside the vessel, toward spherical mirror -->
! reflection or screening --> detector hit (virtual) or screening -->
! channel number 

      REAL R(3), DIR(3)
      INTEGER K0
      REAL Y
      INTEGER ICH
#include "riter/RIC_SUPM.fi"

      LOGICAL OK
      REAL R1(3), DIR1(3), R2(3), DIR2(3)
      LOGICAL IS_PBS_ONWAY, IS_EBS_ONWAY

      ICH=-3
      CALL SFM_HIT(R,DIR,K0,PATH,OK)
      OK=OK.AND.PATH.GT.0
      K=K0
      IF(OK) THEN
         CALL NEW_POINT(R,DIR,PATH,R1)
         K=KJE_HIT(R1(2))
      ENDIF
      IF(.NOT.OK.OR.K.NE.K0) THEN
         K0=MOD(K0+1,2)
         CALL SFM_HIT(R,DIR,K0,PATH,OK)
         OK=OK.AND.PATH.GT.0
         IF(.NOT.OK) RETURN
         CALL NEW_POINT(R,DIR,PATH,R1)
      ENDIF
      Y=R1(2)
      IF(IS_PBS_ONWAY(R,DIR,PATH1)) THEN
         IF(PATH1.LT.PATH) RETURN
      ENDIF
      CALL SFM_REFLECT(R1,DIR,K0,DIR1,AA)
      IF(IS_PBS_ONWAY(R1,DIR1,PATH)) RETURN
      IF(R1(1).GT.0.AND.R1(2).LT.0) THEN
         IF(IS_EBS_ONWAY(R1,DIR1,PATH)) RETURN
      ENDIF
      CALL DETECTOR_HIT(R1,DIR1,K0,PATH,ALFA,BETA,ISM,OK)
      IF(.NOT.OK) RETURN
      CALL DET_CHANNEL(ALFA,BETA,ISM,ICH)
! lens acceptance
      IF(ICH.LT.0) RETURN
      CTH=ABS(S_PROD(DIR1,SM_VN(1,ISM)))
      IF(CTH.LT.(.9877)) ICH=-1

      RETURN
      END


      SUBROUTINE RAY_TRACE_1(R,DIR,K0,ICH)

! point on spherical mirror, toward mirror -->
! reflection --> detector hit --> channel number 
! (no screening)

      REAL R(3), DIR(3)
      INTEGER K0
      INTEGER ICH
#include "riter/RIC_SUPM.fi"

      LOGICAL OK
      REAL DIR1(3)
      LOGICAL IS_PBS_ONWAY, IS_EBS_ONWAY

      ICH=-3
      CALL SFM_REFLECT(R,DIR,K0,DIR1,AA)
      CALL DETECTOR_HIT(R,DIR1,K0,PATH,ALFA,BETA,ISM,OK)
      IF(.NOT.OK) RETURN
      CALL DET_CHANNEL(ALFA,BETA,ISM,ICH)
! lens acceptance
      IF(ICH.LT.0) RETURN
      CTH=ABS(S_PROD(DIR1,SM_VN(1,ISM)))
      IF(CTH.LT.(.9877)) ICH=-1

      RETURN
      END



      REAL FUNCTION CHER_ANGLE(P,MASS,REF_IND)
 
! calculates Cerenkov angle of a particle of mass MASS and momentum P
 
      REAL P, MASS, REF_IND
 
      IF(MASS.LE.0.OR.P.LT.0) THEN
         CHER_ANGLE=0.
         RETURN
      ENDIF
      A=P/MASS
      B=A/SQRT(1.+A**2)*REF_IND
      IF(B.GT.1) THEN
         CHER_ANGLE=ACOS(1./B)
      ELSE
         CHER_ANGLE=0.
      ENDIF
      RETURN
      END
                                                                                
