*********************************************************
*
      real function xmass(ptcut)
*
******************************************************************
*     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    *
******************************************************************
*
      implicit none      
*
      include 'rantpl3.inc'

*--   arguments
      real ptcut
*
*---  local variables ----------------
      integer   kmuon

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

*--
      real xm

      integer imin1,jmin1
      real dmmin1
      integer imin2,jmin2
      real dmmin2

*--   4-vectors
      integer k1,k2,k3,k4
      real plep1(4),plep2(4),plep3(4),plep4(4)
      real pz1(4),pz2(4)

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

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

*--   functions
      real pmass2,amass2,amass

*--   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 ------------------------
*
      kmuon = 0
      
      call hf1(11,real(nisomu), 1.)
      if (nisomu.le.0) return

*--   list of isolated muons with 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(pmu(3,i).ge.ptcut) then
            kmuon = kmuon + 1
            if (kmuon.le.nlepmax) then
               iplep(kmuon) = i
               idlep(kmuon) = idmu
               chglep(kmuon) = chgmu(i)
               do k=1,3
                  plep(k,kmuon) = pmu(k,i)
               enddo
            endif
         endif
      enddo
      if (kmuon.lt.2) return

*--   invariant mass plot of muon pairs
      do i=1,kmuon-1
         ii = iplep(i)
         do j=i+1, kmuon
            jj = iplep(j)
            xm = pmass2(pmu(1,ii),pmu(1,jj))
*--         for all muon pairs
            call hf1(21, xm, 1.)
*--         for oppositely charged pairs
            if(chglep(i).ne.chglep(j)) then  
               call hf1(22, xm, 1.)
            endif
         enddo
      enddo
      
*--- end of muon-pair reconstruction -----------------
      
*------------check if there are 4 or more muons 
      if (kmuon.lt.4) return
      
*---  find a pair with invariant mass closest to mz mass
      imin1 = 0
      jmin1 = 0
      do i=1,kmuon-1
         ii = iplep(i)
         do j=i+1, kmuon
            jj = iplep(j)
            xm = pmass2(pmu(1,ii),pmu(1,jj))
            call hf1(23,xm,1.)
            if (idlep(i).eq.idlep(j) .and.
     &           chglep(i).ne.chglep(j)) then
               call hf1(24,xm,1.)
               if (imin1.eq.0 .or. 
     &              abs(xm-mz).lt.dmmin1) then
                  imin1 = i
                  jmin1 = j
                  dmmin1 = abs(xm-mz)
               endif
            endif
         enddo
      enddo
      k1 = iplep(imin1)
      call pttop(1,plep1, pmu(1,k1),pmu(2,k1),pmu(3,k1))
      k2 = iplep(jmin1)
      call pttop(1,plep2, pmu(1,k2),pmu(2,k2),pmu(3,k2))
      call psum(plep1,plep2, pz1)
      call hf1(25, amass2(plep1,plep2),1.)
      call hf1(26, amass2(plep1,plep2),1.)
      
*---  find the 2nd pair that is closest to mz mass------------
      imin2 = 0
      jmin2 = 0
      dmmin2 = 0.
      do i=1,kmuon-1
         if (i.ne.imin1 .and. i.ne.jmin1) then
            ii = iplep(i)
            do j=i+1, kmuon
               if (j.ne.imin1 .and. j.ne.jmin1) then
                  if (idlep(i).eq.idlep(j) .and.
     &                 chglep(i).ne.chglep(j)) then
                     jj = iplep(j)
                     xm = pmass2(pmu(1,ii),pmu(1,jj))
                     if (imin2.eq.0 .or. 
     &                    abs(xm-mz).lt.dmmin2) then
                        imin2 = i
                        jmin2 = j
                        dmmin2 = abs(xm-mz)
                     endif
                  endif
               endif
            enddo
         endif
      enddo
      k3 = iplep(imin2)
      call pttop(1,plep3, pmu(1,k3),pmu(2,k3),pmu(3,k3))
      k4 = iplep(jmin2)
      call pttop(1,plep4, pmu(1,k4),pmu(2,k4),pmu(3,k4))
      call psum(plep3,plep4, pz2)
      call hf1(27, amass2(plep3,plep4),1.)
      call hf1(28, amass2(plep3,plep4),1.)
      call hf2(29, amass2(plep1,plep2),amass2(plep3,plep4),1.)
      
*--- select 1st and 2nd pairs for higgs candidates---------------
      call psum(pz1,pz2, ph)
      xhmass = amass(ph)
      call pttop(2,ph, etah,phih,pth)
c      pth = sqrt(ph(1)**2 + ph(2)**2)
c      thetah = atan2(pth, ph(3))
c      etah = -alog(tan(thetah*0.5))
      call hf1(31, xhmass, 1.)
      call hf1(32, pth, 1.)
      call hf1(33, etah, 1.)
      
      return
      end

      include 'library/amass.f'
      include 'library/amass2.f'
      include 'library/pmass2.f'
      include 'library/psum.f'
      include 'library/pttop.f'



