!-----------------------------------------------------------------------
! Likelihood calculation
! M. Staric, jan-2000, changed: nov-2000
!
! SUBROUTINE LIKELIHOOD - this does all
! SUBROUTINE LOG_LIKE(X,SIG,WT,N,X0,STFS,STFB,A,PLKH,SFOT)
! REAL FUNCTION POIS_LOG(MEAN,N)
! SUBROUTINE NORMALIZE_LKH(PLKH,NHYP,NORM,IDP)                                      
!
!-----------------------------------------------------------------------


      SUBROUTINE LIKELIHOOD(NITER)

      INTEGER NITER

#include "riter/RIC_TRA.fi"

! returns normalised log Likelihoods in /TRACKS/
! input: expected values from /TRACKS/, 
!        cerenkov angles, sigmas and weights from ring_iter
! photons with negative or zero sigma are ignored

! reset
      DO ITRA=1,NTRA
         DO K=0,NHYP-1
            PLKH(K,ITRA)=0.
         ENDDO
         SFOT(ITRA)=0.
      ENDDO
! fill
      CALL GET_RING_TAB_ROWS(NT)
      IF(NITER.GT.0) THEN
         DO K=1,NT
            CALL GET_RING_TAB(K,ITRA,IFOT,SIG,TH_CER,FI_CER,WT) 
            IF(SIG.GT.0) THEN
               CALL LOG_LIKE(TH_CER,SIG,WT,
     &              NHYP,THC(0,ITRA),STFS(0,ITRA),STFB(ITRA),AX3(ITRA),
     &              PLKH(0,ITRA),SFOT(ITRA))
            ENDIF
         ENDDO
      ELSE
         DO K=1,NT
            CALL GET_RING_TAB(K,ITRA,IFOT,SIG,TH_CER,FI_CER,WT) 
            IF(SIG.GT.0) THEN
               CALL LOG_LIKE_0(TH_CER,SIG,WT,
     &              NHYP,THC(0,ITRA),STFS(0,ITRA),STFB(ITRA),AX1(ITRA),
     &              PLKH(0,ITRA),SFOT(ITRA))
            ENDIF
         ENDDO
      ENDIF
! add poissonian part, normalize it
      DO ITRA=1,NTRA
         IF(RECO_FLAG(ITRA).NE.0) THEN
            DO K=0,NHYP-1
               STFE=STFS(K,ITRA)+STFB(ITRA)
               PLKH(K,ITRA)=PLKH(K,ITRA)+POIS_LOG(STFE,SFOT(ITRA))
            ENDDO
            CALL NORMALIZE_LKH(PLKH(0,ITRA),NHYP,NHYP,PID(ITRA))
         ENDIF
      ENDDO
      RETURN
      END



      SUBROUTINE LOG_LIKE(X,SIG,WT,N,X0,STFS,STFB,A,PLKH,SFOT)

      REAL X, SIG, WT           ! cerenkov angle, sigma, weight [in]
      INTEGER N                 ! number of hypotheses [in]
      REAL X0(0:N-1)            ! cerenkov angles for other,e,mu,pi,K,p [in]
      REAL STFS(0:N-1)          ! expected number of cerenkov photons [in]
      REAL STFB                 ! expected number of background photons [in]
      REAL A                    ! background normalisation constant [in]
      REAL PLKH(0:N-1)          ! logLikelihoods (spatial part) [in,out]
      REAL SFOT                 ! number of photons (weighted sum) [in,out]

      PARAMETER (SQPI=2.50662827463)

      IF(X.LE.0) RETURN

      IF(STFB.EQ.0) THEN
         CONST=LOG(SQPI*SIG)
         DO K=0,N-1
            IF(STFS(K).GT.0) THEN
               Z=(X-X0(K))/SIG
               F=-(Z**2/2+CONST)
               PLKH(K)=PLKH(K)+WT*F
            ENDIF
         ENDDO
      ELSE
         B=A*X**3
         FB=LOG(B)
         DO K=0,N-1
            IF(STFS(K).GT.0) THEN
               Z=(X-X0(K))/SIG
               IF(ABS(Z).LT.3) THEN
                  G=EXP(-Z**2/2)/SQPI/SIG
                  F=LOG((STFS(K)*G+STFB*B)/(STFS(K)+STFB))
               ELSE
                  F=LOG(STFB*B/(STFS(K)+STFB))
               ENDIF
            ELSE
               F=FB
            ENDIF
            PLKH(K)=PLKH(K)+WT*F
         ENDDO
      ENDIF
      SFOT=SFOT+WT

      RETURN
      END
      
      
      SUBROUTINE LOG_LIKE_0(X,SIG,WT,N,X0,STFS,STFB,A,PLKH,SFOT)

! no iterations, linear background assumed

      REAL X, SIG, WT           ! cerenkov angle, sigma, weight [in]
      INTEGER N                 ! number of hypotheses [in]
      REAL X0(0:N-1)            ! cerenkov angles for other,e,mu,pi,K,p [in]
      REAL STFS(0:N-1)          ! expected number of cerenkov photons [in]
      REAL STFB                 ! expected number of background photons [in]
      REAL A                    ! background normalisation constant [in]
      REAL PLKH(0:N-1)          ! logLikelihoods (spatial part) [in,out]
      REAL SFOT                 ! number of photons (weighted sum) [in,out]

      PARAMETER (SQPI=2.50662827463)

      IF(X.LE.0) RETURN

      IF(STFB.EQ.0) THEN
         CONST=LOG(SQPI*SIG)
         DO K=0,N-1
            IF(STFS(K).GT.0) THEN
               Z=(X-X0(K))/SIG
               F=-(Z**2/2+CONST)
               PLKH(K)=PLKH(K)+WT*F
            ENDIF
         ENDDO
      ELSE
         B=A*X
         FB=LOG(B)
         DO K=0,N-1
            IF(STFS(K).GT.0) THEN
               Z=(X-X0(K))/SIG
               IF(ABS(Z).LT.3) THEN
                  G=EXP(-Z**2/2)/SQPI/SIG
                  F=LOG((STFS(K)*G+STFB*B)/(STFS(K)+STFB))
               ELSE
                  F=LOG(STFB*B/(STFS(K)+STFB))
               ENDIF
            ELSE
               F=FB
            ENDIF
            PLKH(K)=PLKH(K)+WT*F
         ENDDO
      ENDIF
      SFOT=SFOT+WT

      RETURN
      END
      
      
      REAL FUNCTION POIS_LOG(MEAN,N)
 
! log(Poisson(mean,n))
 
      REAL MEAN
      REAL N
      
      IF(MEAN.EQ.0) THEN
         IF(N.EQ.0) THEN
            POIS_LOG=0.
         ELSE
            POIS_LOG=-1.E10
         ENDIF
      ELSE
         POIS_LOG=N*LOG(MEAN)-MEAN-ALOGAM(N+1.)
      ENDIF
      RETURN
      END



      SUBROUTINE NORMALIZE_LKH(PLKH,NHYP,NORM,IDP)
 
      REAL PLKH(0:NHYP-1)       ! log Likelihoods (on input and output) [in,out]
      INTEGER NHYP              ! number of all hypotheses [in]
      INTEGER NORM              ! number of hypotheses to normalise them to sum=1 [in]
      INTEGER IDP               ! most probable particle ID [out]
 
      N=MIN(NORM,NHYP)
      IDP=0
      IF(N.LE.0) RETURN
      PMA=PLKH(IDP)
      DO I=0,N-1
         IF(PLKH(I).GT.PMA) THEN
            PMA=PLKH(I)
            IDP=I
         ENDIF
      ENDDO
      DO I=0,NHYP-1
         PLKH(I)=PLKH(I)-PMA
      ENDDO
      C=0.
      DO I=0,N-1
         C=C+EXP(PLKH(I))
      ENDDO
      C=LOG(C)
      DO I=0,NHYP-1
         PLKH(I)=PLKH(I)-C      ! EXP(PLKH()) is probability
      ENDDO
      RETURN
      END
 
 
   
