C----------------------------------------------------------------------
C       requested default input cards:
C
C KTJET_etamax  3.2
C KTJET_ptmin   2.5
C KTJET_IMODE : 3111
C KTJET_IRECO :    1
C KTJET_RCUT1 : 0.1
C KTJET_RCUT2 : 0.5
C KTJET_RMAC  : 1.0
C KTJET_SCALE : -3.
C KTJET_SCVAR : 0.8
C
C----------------------------------------------------------------------
C----------------------------------------------------------------------
      subroutine ktjzes(vtx,sira_flag,sira_prob,ebeam,pbeam)
*
*     'KT algorithm' type jetfinder to run in ZES for DIS jets
*                    interface to KTCLUS package by M.SEYMOUR
*
*     input: tables:
*            SIDAT95 : SINISTRA electron
*            CALTRU  : calorimeter
*            PARTAP
*
*            via parameter:
*            vtx(3)    : vertex 3-vector
*            sira_flag : has Sinistra run allready? (1/0) (yes/no)
*            sira_prob : probability, if SIRA still has to be run
*                        (recommended is: 0.9)
*            ebeam     : momentum of incoming lepton in -z direction
*            pbeam     : momentum of incoming hadron in +z direction
*                        (this depends on the year...)
*                        (note: only |ebeam| and |pbeam| are needed!)
*
*            via input cards:       (recommended:)
*            CUTS ON JETS:
*            eta_cut : KTJET_etamax :  3.2
*            pt_cut  : KTJET_ptmin  :  3.0
*            STEERING OF KTCLUS:
*            Imode  : KTJET_IMODE : 3111
*            Ireco  : KTJET_IRECO :    1
*            rcut   : KTJET_RCUT1 :    0.1      resolution parameter
*                     KTJET_RCUT2 :    0.5      for jetsearch
*            rmac   : KTJET_RMAC  :    1.0      resolution parameter
*                                               for macro jet separation
*            scale  : KTJET_SCALE :             denominator of kt measure
*                                   >  0.  -> use value of "scale" FIXED
*                                   =  0.  -> use KTCLUS default (E_tot)
*                                   = -1.  -> use (E_tot-E_el)*KTJET_SCVAR
*                                   = -2.  -> use E_t_el*KTJET_SCVAR
*                                   = -3.  -> use Q_el*KTJET_SCVAR
*                                                 from SINISTRA electron
*            scvar  : KTJET_SCVAR :    0.8      scale correspondence to Q_da
*
*     output: common /ktjjets/
*
*     called routines:  sira95 (optional)
*                       cccapo
*                       cccxyz
*                       vfill
*                       KTBREI )
*                       KTCLUS ) from KTCLUS package
*                       KTRECO )
*
*     Author: St.Eisenhardt
*     Date  : 06/11/97
*
C----------------------------------------------------------------------

      Implicit NONE

#include "partap.inc"
#include "sidat95.inc"
#include "caltru.inc"
#include "ktjjets.inc"
#include "ktcut.inc"
#include "zescommon.inc"
#include "zestmpcommon.inc"

C-------------------------------------------------------------------
C-------  variable & parameter definition
C-------------------------------------------------------------------
C-------  GLOBAL STEERING FLAGS/PARMETERS
C-------  ===============================
C------- max_njin is the maximum number of 4-vectors KTCLUS can handle

      Integer           nj,j,i,tt
      Integer           vdim
      Integer           Ierr

      Integer           max_njin
      Parameter         (max_njin=512)

      Integer           sira_flag
      Real              sira_prob
      Real              ktpi,vtx(3)
      Parameter         (ktpi=3.141592654)
      Double Precision  dktpi


C-------  VARIABLES for electron handling
C-------  ===============================

      Real              eesira,e_r,e_theta,e_phi
      Real              siraHitPosVTXcor(3)
      Integer           FindisOpt
      Integer           neposer
      Integer           ecand,eposer(100)


C-------  VARIABLES for CAL EXTRACTION
C-------  ============================

      Real              cell_x,cell_y,cell_z,cell_e
      Real              cell_r,cell_theta,cell_phi
      Real              cell_E_x,cell_E_y,cell_E_z

C-------  JET VARIABLES & jet algorithm variables
C-------  =======================================

      Integer           njin,njout,nsubjet,njet
      Integer           Iboost,Imode,Ireco
      Real              scale,scvar,kt_etot,redge(max_njin)
      Real              rmac
      Integer           subjet1(max_njin)
      Double Precision  djin_l(4,max_njin)
      Double Precision  djin_b(4,max_njin)
      Double Precision  djout_l(4,max_njin)
      Double Precision  djout_b(4,max_njin)
      Real              ebeam,pbeam
      Double Precision  debeam,dpbeam,descat(4)
      real              rcut(2)

C-------  INTERMEDIATE KINEMATIC VARIABLES
C-------  ================================

*      Double Precision  dene(max_njin)
*      Double Precision  det(max_njin)
      Double Precision  dpt(max_njin)
      Double Precision  dpz(max_njin)
*      Double Precision  dmj(max_njin)
*      Double Precision  depzj(max_njin)
      Double Precision  deta(max_njin)
      Double Precision  dth(max_njin)
      Double Precision  dphi(max_njin)

C-------------------------------------------------------------------
C-------  initialization
C-------------------------------------------------------------------

C-------  COMMON block variables
C-------  ======================
C------- ktj_maxjets is the maximum number of jets allowed
C------- to be stored in the common block KTJJETS

      call vfill(njkt,ktj_maxcuts,0)
      vdim = ktj_maxcuts * ktj_maxjets * 4
      call vfill(pjkt,vdim,-999.99)


C-------  VARIABLES for electron handling
C-------  ===============================

      eesira = 0.0
      e_r = 0.0
      e_theta = 0.0
      e_phi = 0.0
      call vfill(siraHitPosVTXcor,3,0.0)

C-------  JET VARIABLES & jet algorithm variables
C-------  =======================================

      njin = 0
      njout = 0
      njet = 0
      nsubjet = 0
      dktpi = dble(ktpi)

      scale = 0.
      scvar = 0.
      kt_etot = 0.

      vdim = 4 * max_njin * 2           ! * 2 for double variables
      call vfill(djin_l,vdim,0.0)
      call vfill(djin_b,vdim,0.0)
      call vfill(djout_l,vdim,0.0)
      call vfill(djout_b,vdim,0.0)
      vdim = 4 * 2                      ! * 2 for double variables
      call vfill(descat,vdim,0.0)

      call vfill(redge,max_njin,0.0)
      call vfill(subjet1,max_njin,0.0)

      vdim = max_njin * 2               ! * 2 for double variables
*      call vfill(dene,vdim,0.0)
*      call vfill(det,vdim,0.0)
      call vfill(dpt,vdim,0.0)
*      call vfill(dmj,vdim,0.0)
*      call vfill(depzj,vdim,0.0)
      call vfill(deta,vdim,0.0)
      call vfill(dth,vdim,0.0)
      call vfill(dphi,vdim,0.0)

C-------------------------------------------------------------------
C-------------------------------------------------------------------
C-------  start
C-------------------------------------------------------------------
C-------------------------------------------------------------------

C-------  check for events with # of cell > 512

      if ( coutab(CALTRU).gt.512 ) then
         write(*,*) 'KTJZES: too many cells in CALTRU =',coutab(CALTRU)
         RETURN
      endif

C-------------------------------------------------------------------
C-------------------------------------------------------------------
C-------  get 4-vector of scattered electron
C-------------------------------------------------------------------
C-------------------------------------------------------------------

C-------  check if SINISTRA95 has allready been run

C      if ( sira_flag.ne.1 ) then
C         call sira95(vtx,sira_prob,Ierr)
C      endif

C-------  check if at least one electron has been found .or. SIRA error

C      if ( Ncand.lt.1 .or. Ierr.ne.0 ) then
C         RETURN
C      endif

C-------  recommended electron finding for DIS:
C-------            highest prop. and E -> ( FindisOpt = 1 )

C      FindisOpt = 1

C-------  find most propable electron candidate ( ecand )
C-------  and corresponding number of CAL cells ( neposer )

C      call findis95(FindisOpt,sira_prob,ecand,Ierr)

C-------  there seems to be a bug in FINDIS95 ( ecand=0 gives no error... )

C      if ( Ierr.ne.0 .or. ecand.le.0 ) then
C         RETURN
C      endif

C-------  electron energy stays uncorrected in ZES as for TLT
C-------  for the hit position of the electron on the CAL surface
C-------  one has to correct for the vertex
      eesira = ZES_Eeu_si
      siraHitPosVTXcor(1) = ZES_Xeu_si-vtx(1)
      siraHitPosVTXcor(2) = ZES_Yeu_si-vtx(2)
      siraHitPosVTXcor(3) = ZES_Zeu_si-vtx(3)
      neposer = ZESTMP_Ncell_si

      call cccapo(siraHitPosVTXcor(1),
     &            siraHitPosVTXcor(2),
     &            siraHitPosVTXcor(3),
     &            e_r,e_theta,e_phi)

      descat(1) = dble( eesira * sin(e_theta) * cos(e_phi) )
      descat(2) = dble( eesira * sin(e_theta) * sin(e_phi) )
      descat(3) = dble( eesira * cos(e_theta) )
      descat(4) = dble( eesira )

C-------  make list of electron cells to be excluded from CALTRU
C-------  get poser numbers of electron related CAL cells

	 do j = 1,ZESTMP_Ncell_si
            eposer(j) = ZESTMP_CellList_si(j)
	 enddo

C-------------------------------------------------------------------
C-------------------------------------------------------------------
C-------  get 4-vectors from CAL
C-------------------------------------------------------------------
C-------------------------------------------------------------------

C-------  extract 4-vectors from CALTRU

      do 100 j = 1,coutab(CalTru)

         call fettab(CalTru,id,j)

C-------  skip cells beloning to "found electron cells"

	 do i = 1,neposer
	    if ( CalTru_CellNr.eq.eposer(i) ) goto 100
	 enddo

C-------  get cell geometry: x, y, z, r, theta, phi and correct for vertex

         call cccxyz(CalTru_CellNr,cell_x,cell_y,cell_z,Ierr)

         if (Ierr.ne.0) then
	    write(*,*) 'KTJZES: ERROR in CCCXYZ : Ierr = ',Ierr
	    RETURN
         endif

         cell_x = cell_x - vtx(1)
         cell_y = cell_y - vtx(2)
         cell_z = cell_z - vtx(3)
         call cccapo(cell_x,cell_y,cell_z,
     &        cell_r,cell_theta,cell_phi)

C-------  get projections of cell energy

         cell_E_x = 0.0
         cell_E_y = 0.0
         cell_E_z = 0.0
         cell_E = 0.0

         cell_E_x = CalTru_E * sin(cell_theta) * cos(cell_phi)
         cell_E_y = CalTru_E * sin(cell_theta) * sin(cell_phi)
         cell_E_z = CalTru_E * cos(cell_theta)
         cell_E = CalTru_E

C-------  fill 4-vectors for jet finding

         njin = njin + 1

         djin_l(1,njin) = dble(cell_E_x)
         djin_l(2,njin) = dble(cell_E_y)
         djin_l(3,njin) = dble(cell_E_z)
         djin_l(4,njin) = dble(cell_E)

         kt_etot = kt_etot + djin_l(4,njin)

 100  continue

C-------------------------------------------------------------------
C-------------------------------------------------------------------
C-------  jet finding
C-------------------------------------------------------------------
C-------------------------------------------------------------------

C-------  run KT algorithm
C-------  ================

C-------  scale for jet finding:
C-------        scale >  0.  -> use "scale" FIXED
C-------              =  0.  -> use KTCLUS default (E_tot)
C-------              = -1.  -> use (E_tot-E_el)*KTJET_SCVAR
C-------              = -2.  -> use E_t_el*KTJET_SCVAR
C-------              = -3.  -> use Q_el*KTJET_SCVAR
C-------                            from SINISTRA electron

      scale = KTJET_SCALE
      scvar = KTJET_SCVAR

      if ( scale.lt.0. ) then

         if ( scale.eq.-1. ) then
            scale = ( kt_etot-sngl(descat(4)) )*scvar

         elseif ( scale.eq.-2. ) then
            scale = sngl(dsqrt( descat(1)**2 + descat(2)**2 )) * scvar

         elseif ( scale.eq.-3. ) then
            scale = ( 2.*abs(ebeam)*sngl(descat(4)) *
     &                (1.+cos(e_theta)) ) *                ! Q_el
     &              scvar
         else
            scale = 0.
            write(*,*)
     &      'KTJZES: Warning: wrong scale -> KTCLUS default choosen!'
         endif

      endif

C-------  z components of beam vectors

      debeam = dble( ebeam )
      dpbeam = dble( pbeam )

C-------  boost input vectors to BREIT frame : djin_l(k,j) -> djin_b(k,j)
C-------    input : Iboost      = 0/2=to BREIT frame, 1/3=from BREIT frame
C-------                          0/1=no azimuthal rotation afterwards
C-------                          2/3=lepton plane rotated into the x-z plane
C-------            debeam      = momentum of incoming lepton in +z direction
C-------            dpbeam      = momentum of incoming hadron in +z direction
C-------            descat(k)   = 4-momentum of outgoing lepton
C-------            njin        = number of particles in djin(k,j)
C-------            djin_l(k,j) = 4-momentum of jth particle
C-------    output: djin_b(k,j) = 4-momentum of jth particle

      Iboost = 0                        ! LAB -> BREIT

      call KTBREI(Iboost,debeam,dpbeam,descat,njin,
     &            djin_l,djin_b,*501)

C-------  do clustering :
C-------    input : Imode    = 3111 
C-------                     = z-> : p-><-e , angular kt def., derive angles
C-------                       relative to proton dir., E recombination scheme
C-------            djin_b(k,j) = 4-momentum of jth particle
C-------            njin     = number of particles
C-------            scale    = denominator of kt measure. etot if ecut=0
C-------    output: redge(j) = rcut edge for j -> j-1 jets

      Imode = KTJET_IMODE
      call KTCLUS(Imode,djin_b,njin,scale,redge,*502)

C-------  reconstruct jet momenta at rcut
C-------    input : Ireco    = recombination scheme (may differ to KTCLUS)
C-------            djin_b(k,j) = 4-momentum of jth particle
C-------            njin     = number of particles
C-------            scale    = denominator of kt measure. etot if ecut=0
C-------            rcut     = rcut value where to reconstruct jet momenta
C-------            rmac     = rcut value to determine which jets are sub-jets
C-------    output: djout_b(k,nj) = 4-momentum of njth jet at scale rcut
C-------            subjet1(j) = macro-jet containing jth jet, 0 if macro-jet
C-------            njet     = number of jets
C-------            nsubjet  = number of sub-jets

C------- At this point we can reconstruct jets at more than one value
C------- of rcut. Currently the only input variable that changes is rcut.
C------- For ZES presently two different rcuts are used.

      Ireco = KTJET_IRECO
      rmac  = KTJET_RMAC
      rcut(1) = KTJET_RCUT1
      rcut(2) = KTJET_RCUT2
     
      do i = 1,2

C------- Add following three lines for intialization though it may make no
C------- difference. 

         vdim = 4 * max_njin * 2        ! * 2 for double variables
         call vfill(djout_l,vdim,0.0)
         call vfill(djout_b,vdim,0.0)

         call KTRECO(Ireco,djin_b,njin,scale,rcut(i),rmac,
     &               djout_b,subjet1,njet,nsubjet,*503)


C-------  boost output vectors to LAB frame : djout_b(k,j) -> djout_l(k,j)

         Iboost = 1             ! BREIT -> LAB

         call KTBREI(Iboost,debeam,dpbeam,descat,njet,
     &               djout_b,djout_l,*501)

C-------  calculate jet properties for LAB frame

         njout = 0

         do nj = 1,njet

C-------  energy
*         dene(nj) = djout_l(4,nj)
C-------  pz
          dpz(nj) = djout_l(3,nj)


C-------  pt
            dpt(nj) = dsqrt ( djout_l(1,nj)**2 +
     &                        djout_l(2,nj)**2 )

C-------  mass
*         dmj(nj) = dsqrt ( djout_l(4,nj)**2 -
*     &                     ( djout_l(1,nj)**2 +
*     &                       djout_l(2,nj)**2 +
*     &                       djout_l(3,nj)**2 ) )

C-------  theta
            if ( dpt(nj).eq.0.d0 .and. 
     &           djout_l(3,nj).eq.0.d0 ) then
               dth(nj) = dktpi / 2.d0
            else
               dth(nj) = datan2( dpt(nj),djout_l(3,nj) )
               if (dth(nj).lt.0.d0 ) then
                  dth(nj) = dth(nj) + dktpi
               endif
            endif

C-------  eta = pseudo rapidity
            if ( dth(nj).eq.0.d0 ) then
               deta(nj) = 999.
            else
               deta(nj) = -1.d0 * dlog(dtan( dth(nj)/2.d0 ))
            endif

C-------  Et
*            det(nj) = djout_l(4,nj) * dsin(dth(nj))

C-------  E - pz
*           depzj(nj) = djout_l(4,nj) * ( 1.d0 - dcos(dth(nj)) )

C-------  phi
            if ( djout_l(1,nj).eq.0.d0 .and. 
     &           djout_l(2,nj).eq.0.d0 ) then
               dphi(nj) = 0.d0
            else
               dphi(nj) = datan2( djout_l(2,nj),djout_l(1,nj) )
               if (dphi(nj).lt.0.d0 ) then
                  dphi(nj) = dphi(nj) + 2.d0 * dktpi
               endif
            endif

C-------  fill common block variables

C-------  check cuts to accept the valid jets

            if ( sngl(deta(nj)).le.KTJET_etamax .and.
     &           sngl(dpt(nj)).ge.KTJET_ptmin ) then
               
               njout = njout + 1

C-------  check for events with njout > ktj_maxjets
C-------        -> throw away left over jets

               if ( njout.le.ktj_maxjets ) then
                  njkt(i) = njout                  
                  pjkt(i,njout,1) = sngl( dpt(nj) )
                  pjkt(i,njout,2) = sngl( deta(nj) )
                  pjkt(i,njout,3) = sngl( dphi(nj) )
                  pjkt(i,njout,4) = sngl( dpz(nj) )
               endif

            endif

*         write(*,*)
*     &   'KTJZES:',nj,njout,njkt

         enddo

      enddo

      RETURN

 501  write(25,'(a26,i4)')
     &     ' ERROR in KTBREI : njin = ',njin
      RETURN

 502  write(25,'(a26,2i4)')
     &     ' ERROR in KTCLUS : njin = ',njin
      RETURN

 503  write(25,'(a26,2i4)')
     &     ' ERROR in KTRECO : njin = ',njin
      RETURN

 510  continue
      RETURN

      END
