*******************************************************
      SUBROUTINE ZUFODIF(VPX,VPY,VPZ,ETAF,ETAB,
     >                     ETAMAX,HPX1,HPY1,HPZ1,HE1)

*******************************************************
*
* INPUT: VERTEX X,Y,Z: VPX,VPY,VPZ
*
* OUTPUT:
*
*       ZES name:      Local name:
*      ------------    ----------
*      ZES_Fgap=       ETAF
*      ZES_Bgap=       ETAB
*      ZES_Etamax_zu = ETAMAX
*      ZES_V_H_px_zu = HPX1  	
*      ZES_V_H_py_zu = HPY1  	
*      ZES_V_H_pz_zu = HPZ1  
*      ZES_V_H_E_zu  = HE1  	
*
***********************************
        IMPLICIT NONE

#include "zisles.inc"
#include "zdrecgb.inc"


      REAL VP(3),ETA(5000),VPX,VPY,VPZ
      REAL ZTHE,GAP,DIF2,PIPPO
      INTEGER I,NETA,J,IND2,IERR
      REAL ETAF,ETAB,ETAMAX
      REAL HE1,HPX1,HPY1,HPZ1



*-->  Zero Hadronic Info

      HE1  = 0.
      HPX1 = 0.
      HPY1 = 0.
      HPZ1 = 0.
      ETAF = -1000
      ETAB = -1000
      ETAMAX = -1000
      NETA=0.      
      CALL VZERO(ETA,5000)

      VP(1)=VPX
      VP(2)=VPY
      VP(3)=VPZ

* I assume ZUFO has already been called, so I can just do: 

      CALL Z_RECGB(1,1,VP,IERR)


      IF (IERR.LT.0) THEN
       PRINT *,'*-2-> Z_RECGB(EVENT REJECTED) ERROR : ',IERR
       RETURN
      ENDIF

* Loop over ZUFOS and  calculate the rapidity and x,y,z for each ZUFO

      DO I=1,NZUFOS
       ZTHE = ATAN2(SQRT(ZUFO(1,I)**2+ZUFO(2,I)**2),ZUFO(3,I))
       ETA(I)=-LOG(TAN(ZTHE/2))
       HPX1 = HPX1 + ZUFO(1,I)
       HPY1 = HPY1 + ZUFO(2,I)
       HPZ1 = HPZ1 + ZUFO(3,I)
       HE1  = HE1  + ZUFO(4,I)
      ENDDO

C Here I order the eta array in order 
C i=1 (max) --> i - nzufo (min) pseudorap.

      NETA=NZUFOS

      Do  I = 1,NETA
         DO  J=I+1, NETA
            IF(ETA(J).GT.ETA(I)) THEN
               PIPPO=ETA(I)
               ETA(I)=ETA(J)
               ETA(J)=PIPPO
            ENDIF
         ENDDO
      ENDDO
c
c Compute the forward edge of the gap and put it in IND2
c

      GAP = ETA(1)-ETA(2)
      IND2 = 1
      DO I=2,NETA-1
         j=i+1
         DIF2=ETA(I)-ETA(j)
         IF(DIF2.GT.GAP) THEN
            GAP=DIF2
            IND2=I
         ENDIF
      ENDDO

* Fill the GAP block

      ETAF = ETA(IND2)
      IF (NETA.GT.1) ETAB = ETA(IND2+1)
      ETAMAX = ETA(1)

      return
      end


**********************************************************************