*********************************************************
*
      real function zz4l()
*
******************************************************************
*     Sample program of Higgs -> 4 lepton analysis               *
*     (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    *
*                                                                *
*  TWO PARTS OF THIS PROGRAM MUST BE ADDED BY YOU                *
*  2002.8.20 modified to exclude Z->tautau cases                 *
*                           (pointed by Kawachi (Kyoto Ed U)     *
******************************************************************
*
      implicit none      
*
      include 'rantpl3.inc'

*--   parameters
      real ptmucut, ptelcut
      parameter (ptmucut = 20., ptelcut = 20.)
      real etamucut, etaelcut
      parameter (etamucut = 2.5, etaelcut = 2.5)
      real dmzcut
      parameter (dmzcut = 10.)
      
*---  local variables ----------------
      integer   kmuon

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

*--
      real xmll

      integer ilep1,ilep2
      integer ilep3,ilep4
      real dmmin

*--   4-vectors
      real pz1(4),pz2(4)
      real xmz1, xmz2

*--   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

*      itagmu(i) : tag number of selected i-th muon
*      itag2(1,j), itag2(2,j) : tag number of selected j-th muon pair

*-------------- 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 pttop(1,plep(1,nlep), 
     &                 pmu(1,i),pmu(2,i),pmu(3,i))
               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
         

PLEASE ADD YOUR PROGRAM FOR THIS PART FOR ELECTRONS.
YOU CAN COPY MOST OF THE PART FOR MUONS JUST ABOVE HERE.

         enddo
         
      endif

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

*------------check if there are 4 or more muons 
      if (nlep.lt.4) 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(41,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 .or. ilep2.eq.0) return
      call psum(plep(1,ilep1),plep(1,ilep2), pz1)
      xmz1 = amass(pz1)
      call hf1(42, xmz1, 1.)
      
*---  find the 2nd pair that is closest to mz mass------------
      ilep3 = 0
      ilep4 = 0
      dmmin = 0.
      do i=1,nlep-1



PLEASE ADD YOUR PROGRAM FOR THIS PART FOR 
FINDING THE 2ND PAIR THAT IS CLOSEST TO MZ MASS.
TAKE A LOOK OF xmass.f AS A GOOD GUIDLINE FOR THIS PART.



      enddo
      if(ilep3.eq.0 .or. ilep4.eq.0) return
      call psum(plep(1,ilep3),plep(1,ilep4), pz2)
      xmz2 = amass(pz2)
      call hf1(44, xmz2, 1.)
      call hf2(45, xmz1, xmz2, 1.)

*---  cut onf delta_mz
      if (abs(xmz1-mz).gt.dmzcut) return
      if (abs(xmz2-mz).gt.dmzcut) return
      
*--- select 1st and 2nd pairs for higgs candidates---------------
      call psum(pz1,pz2, ph)
      xmh = amass(ph)
      call pttop(2,ph, etah,phih,pth)

      call hf1(51, xmh, 1.)
      call hf1(52, pth, 1.)
      call hf1(53, etah, 1.)
      
      return
      end

      include 'library/subroutines.inc'
