*********************************************************
*
      real function zzllbb()
*
******************************************************************
*     Sample program of Higgs -> l l b b                         *
*     (1) Form invariant mass of muon pairs                      *
*     (2) Select opposite charged pairs and store                *
*     (3) Select 2 pairs closest to mass(Z)                      *
*     (4) Form the invariant masss of selected two muon pairs    *
******************************************************************
*
      implicit none      
*
      include 'rantpl3.inc'

*--   parameters

*--   for leptons
      real ptmucut, etamucut
      parameter (ptmucut = 20., etamucut = 2.5)
      real ptelcut, etaelcut
      parameter (ptelcut = 20., etaelcut = 2.5)
      
      integer nlepcut, nbjetcut
      parameter (nlepcut = 2, nbjetcut = 2)

*--   for ets
      real ptbjcut, etabjcut
      parameter (ptbjcut = 20., etabjcut = 2.5)
      real ptqjcut, etaqjcut
      parameter (ptqjcut = 20., etaqjcut = 5.0)
            
      real dmzllcut, dmzbbcut
      parameter (dmzllcut = 10., dmzbbcut = 25.)

*--   list of leptons      
      integer nlep,nlepmax
      integer nmuo,nele
      parameter (nlepmax = 40)
      integer iplep(nlepmax),idlep(nlepmax),chglep(nlepmax)
      real plep(7,nlepmax)
      integer idmu,idel
      parameter (idmu=1,idel=2)

*--   list of jets
      integer njetmax
      parameter (njetmax = 100)
      integer nqjet,ipqjet(njetmax)
      real pqjet(7,njetmax)
      integer nbjet,ipbjet(njetmax)
      real pbjet(7,njetmax)
      integer idbj,idtauj,idqj
      parameter (idbj=5,idtauj=15,idqj=98)

*--
      real xmll, xmbb

      integer ilep1,ilep2
      integer ibjet1,ibjet2
      real dmmin

*--   4-vectors
      real pzll(4),pzbb(4)
      real xmzll, xmzbb

*--   for higgs
      real ph(4)
      real xmh, pth, thetah, etah, phih

*--   functions
      real amass2,amass

*--   constants
      real mz
      parameter (mz = 91.187)

*--   temporary variables
      integer i,j,k
      integer ii,jj

      integer nsetbje
      parameter (nsetbje = 5)
      integer iftau,ifveto
      parameter (iftau=0, ifveto=0)
      real epsitau
      parameter (epsitau = 0.5)
      integer indveto
      parameter (indveto = 2)


      call atlfastb(nsetbje,iftau,epsitau,ifveto,indveto)
      call jetcorr

*-------------- initialization ------------------------
*
      nlep = 0
      nmuo = 0
      nele = 0
      
      call hf1(11,real(nisomu), 1.)
      
      if (nisomu.gt.0) then

*--   list of isolated muons with |eta|<etacut and pt>ptcut
         do i=1,nisomu
            call hf1(12, pmu(3,i), 1.)
            call hf1(13, pmu(1,i), 1.) 
            call hf1(14, pmu(2,i), 1.)
            if(abs(pmu(1,i)).lt.etamucut .and.
     &           pmu(3,i).gt.ptmucut) then
               nlep = nlep + 1
               nmuo = nmuo + 1
               if (nlep.le.nlepmax) then
                  iplep(nlep) = i
                  idlep(nlep) = idmu
                  chglep(nlep) = chgmu(i)
                  call pcopy37(pmu(1,i),plep(1,nlep))
               endif
            endif
         enddo
         
      endif
         
      call hf1(21,real(nel), 1.)
      
      if (nel.gt.0) then
      
*--   list of electrons with |eta|<etacut and pt>ptcut
         do i=1,nel
            call hf1(22, pel(3,i), 1.)
            call hf1(23, pel(1,i), 1.) 
            call hf1(24, pel(2,i), 1.)
            if(abs(pel(1,i)).lt.etaelcut .and.
     &           pel(3,i).gt.ptelcut) then
               nlep = nlep + 1
               nele = nele + 1
               if (nlep.le.nlepmax) then
                  iplep(nlep) = i
                  idlep(nlep) = idel
                  chglep(nlep) = chgel(i)
                  call pcopy37(pel(1,i),plep(1,nlep))
               endif
            endif
         enddo
         
      endif

      call hf1(31,real(nlep),1.)
      call hf1(32,real(nmuo),1.)
      call hf1(33,real(nele),1.)

*------------check if there are 2 or more leptons
      if (nlep.lt.nlepcut) return

      nqjet = 0
      nbjet = 0

      call hf1(41,real(njt),1.)

      do i=1,njt
         if (idjt(i).eq.idqj) then
            call hf1(42, pjt(3,i), 1.)
            call hf1(43, pjt(1,i), 1.) 
            call hf1(44, pjt(2,i), 1.)
            if(abs(pjt(1,i)).lt.etaqjcut .and.
     &           pjt(3,i).gt.ptqjcut) then
               nqjet = nqjet + 1
               if (nqjet.le.njetmax) then
                  ipqjet(nqjet) = i
                  call pcopy37(pjt(1,i),pqjet(1,nqjet))
               endif
            endif
         elseif (idjt(i).eq.idbj) then
            call hf1(45, pjt(3,i), 1.)
            call hf1(46, pjt(1,i), 1.) 
            call hf1(47, pjt(2,i), 1.)
            if(abs(pjt(1,i)).lt.etabjcut .and.
     &           pjt(3,i).gt.ptbjcut) then
               nbjet = nbjet + 1
               if (nbjet.le.njetmax) then
                  ipbjet(nbjet) = i
                  call pcopy37(pjt(1,i),pbjet(1,nbjet))
               endif
            endif
         endif
      enddo

      call hf1(51,real(nqjet),1.)
      call hf1(52,real(nbjet),1.)

      if (nbjet.lt.nbjetcut) return

*---  find a pair with invariant mass closest to mz mass
      ilep1 = 0
      ilep2 = 0
      dmmin = 0.
      do i=1,nlep-1
         do j=i+1, nlep
            xmll = amass2(plep(1,i),plep(1,j))
            if (idlep(i).eq.idlep(j) .and.
     &           chglep(i).ne.chglep(j)) then
               call hf1(61,xmll,1.)
               if (ilep1.eq.0 .or. 
     &              abs(xmll-mz).lt.dmmin) then
                  ilep1 = i
                  ilep2 = j
                  dmmin = abs(xmll-mz)
               endif
            endif
         enddo
      enddo

      if (ilep1.eq.0) return

      call psum(plep(1,ilep1),plep(1,ilep2), pzll)
      xmzll = amass(pzll)
      call hf1(62, xmzll, 1.)

*---  cut onf delta_mz
      if (abs(xmzll-mz).gt.dmzllcut) return
      
*---  find the 2nd pair that is closest to mz mass------------
      ibjet1 = 0
      ibjet2 = 0
      dmmin = 0.
      do i=1,nbjet-1
         do j=i+1, nbjet
            xmbb = amass2(pbjet(1,i),pbjet(1,j))
            call hf1(63,xmbb,1.)
            if (ibjet1.eq.0 .or. 
     &           abs(xmbb-mz).lt.dmmin) then
               ibjet1 = i
               ibjet2 = j
               dmmin = abs(xmbb-mz)
            endif
         enddo
      enddo

      if (ibjet1.eq.0) return
      
      call psum(pbjet(1,ibjet1),pbjet(1,ibjet2), pzbb)
      xmzbb = amass(pzbb)
      call hf1(64, xmzbb, 1.)
      call hf2(65, xmzll, xmzbb, 1.)

      if (abs(xmzbb-mz).gt.dmzbbcut) return
      
*--- select 1st and 2nd pairs for higgs candidates---------------
      call psum(pzll,pzbb, ph)
      xmh = amass(ph)
      call pttop(2,ph, etah,phih,pth)

      call hf1(71, xmh, 1.)
      call hf1(72, pth, 1.)
      call hf1(73, etah, 1.)
      
      return
      end
      include 'library/subroutines.inc'
