Subroutine MFH123(iopt,nn)
C    ===========================
C
      Implicit None
C
#include "partap.inc"
#include "zdskey.inc"
#include "fmckin.inc"
#include "fwmip.inc"
#include "fwtrk.inc"
#include "fwhmu.inc"
#include "mfsigm.inc"
#include "mfcuts.inc"
#include "mfmtch.inc"
#include "mfflg.inc"
c
      integer iucomp

      logical ok
c
      integer i1,i2,i3,j,nn,nloop(3),iopt,imask(4),i
      real x1,y1,z1,x2,y2,z2,x3,y3,z3,x12,y12,z12,x13,y13,z13,x23,y23,
     .    z23,cx12,cy12,cz12,cx13,cy13,cz13,cx23,cy23,cz23,r12,r13,r23,
     .    cst,snt,csp,snp,aaa,s2x,s2y,sig2xctd,sig2yctd,r3sq,pmu2,cms,
     .    cms12,cms23,cms13,swt,swp,chip,chit,phav,thav,theta(4),phi(4),
     .    err2p(4),err2t(4),th,ph,chisq,resp
      integer offset,kk,ii(4)
      real prob
      real probtp,probt,probp
c
      real pi,twopi,pi4
      save pi,twopi,pi4
c
      logical first/.true./
      save    first
c
      real delta_cut /0.25/
      save delta_cut
      integer ithz
      real sigt(4,4,3),sigp(4,4,3),scale
      save sigt,sigp,scale
      data scale /2./
      data sigt/
c iopt=1
     & 0.020,0.020,0.020,0.020,
     & 0.020,0.020,0.020,0.020,
     & 0.110,0.065,0.040,0.030,
     & 0.060,0.030,0.015,0.010,
c iopt=2
     & 0.015,0.015,0.015,0.015,
     & 0.020,0.020,0.020,0.020,
     & 4*0,4*0,
c iopt=3
     & 4*0,
     & 0.050,0.030,0.020,0.020,
     & 4*0,
     & 0.085,0.040,0.020,0.010/
      data sigp/
c iopt=1
     & 0.060,0.060,0.060,0.060,
     & 0.025,0.025,0.025,0.025,
     & 0.150,0.110,0.080,0.050,
     & 0.095,0.050,0.030,0.030,
c iopt=2
     & 0.050,0.050,0.050,0.050,
     & 0.050,0.050,0.050,0.050,
     & 4*0,4*0,
c iopt=3
     & 4*0,
     & 0.070,0.040,0.025,0.025,
     & 4*0,
     & 0.070,0.040,0.025,0.025/
c
c
      if(first) then
       first=.false.
       pi=acos(-1.)
       twopi=2.*pi
       pi4=pi/4.
      endif
c
      call vzero(theta,4)
      call vzero(phi,4)

      nn=0
c
      if(Mdebug.gt.0) then
        write(6,*) ' NUMBER OF FMU ELEMENTS ',nhmf
        write(6,*) ' NUMBER OF MIP ELEMENTS ',nmip(1)
        write(6,*) ' NUMBER OF CTD ELEMENTS ',nctd(1)
      endif
c
      if(iopt.eq.1) then
        nloop(1)=nhmf
        nloop(2)=nmip(1)
        nloop(3)=nctd(1)
        imask(1)=1
        imask(2)=1
        imask(3)=1
        imask(4)=1
        kk=4
        ii(1)=1
        ii(2)=2
        ii(3)=3
        ii(4)=4
      elseif(iopt.eq.2) then
        nloop(1)=nhmf
        nloop(2)=nmip(1)
        nloop(3)=1
        imask(1)=1
        imask(2)=1
        imask(3)=0
        imask(4)=0
        kk=2
        ii(1)=1
        ii(2)=2
      elseif(iopt.eq.3) then
        nloop(1)=nhmf
        nloop(2)=1
        nloop(3)=nctd(1)
        imask(1)=0
        imask(2)=1
        imask(3)=0
        imask(4)=1
        kk=2
        ii(1)=2
        ii(2)=4
      endif
c
c---> loop over mu candidates
      do i1=1,nloop(1)
         x1=xhmf(1,i1)
         y1=xhmf(2,i1)
         z1=xhmf(3,i1)
         if(iopt.eq.2) then
          theta(1)=xhmf(4,i1)
          phi(1)=xhmf(5,i1)
          err2t(1)=ehmf(3,i1)
          err2p(1)=ehmf(4,i1)
         endif
         do i2=1,nloop(2)
c------> loop over mip candidates
           if(nmip(1).gt.0.and.iopt.ne.3.and.imipused(i2,1).eq.0) then
ccc           if(nmip(1).gt.0.and.iopt.ne.3) then
            x2=xmip(i2,1)
            y2=ymip(i2,1)
            z2=zmip(i2,1)
            x12=x1-x2
            y12=y1-y2
            z12=z1-z2
            Call mfCaPo(x12,y12,z12,r12,th,ph)
            if(iopt.ne.2) then
             cx12=x12/r12
             cy12=y12/r12
             cz12=z12/r12
             theta(1)=th
             phi(1)=ph
            else
             theta(2)=thmip(i2,1)
             phi(2)=phmip(i2,1)
            endif
           endif
           do i3=1,nloop(3)
c----------> loop over ctd candidates
             if(nctd(1).gt.0.and.iopt.ne.2.and.ictdused(i3,1).eq.0)then
ccc             if(nctd(1).gt.0.and.iopt.ne.2) then
               x3=xctd(i3,1)
               y3=yctd(i3,1)
               z3=zctd(i3,1)
               theta(4)=thctd(i3,1)
               phi(4)=psctd(i3,1)
               if(imask(2).gt.0) then
                x13=x1-x3
                y13=y1-y3
                z13=z1-z3
                Call mfCaPo(x13,y13,z13,r13,th,ph)
                theta(2)=th
                phi(2)=ph
                if(phi(2).lt.0) phi(2)=phi(2)+twopi
               endif
               if(imask(3).gt.0) then
                x23=x2-x3
                y23=y2-y3
                z23=z2-z3
                Call mfCaPo(x23,y23,z23,r23,th,ph)
                cx23=x23/r23
                cy23=y23/r23
                cz23=z23/r23
                theta(3)=th
                phi(3)=ph
               endif
             endif
c
             if(iopt.eq.2) then
                ithz=1
             else
                if(ndfctd(i3,1).lt.10) then
                   ithz=1
                elseif(ndfctd(i3,1).lt.20) then
                   ithz=2
                elseif(ndfctd(i3,1).lt.30) then
                   ithz=3
                else
                   ithz=4
                endif
             endif
             do j=1,kk
                err2t(ii(j))=sigt(ithz,ii(j),iopt)**2
                err2p(ii(j))=sigp(ithz,ii(j),iopt)**2
                if(iopt.eq.1) then
                   if(ii(j).eq.3.or.ii(j).eq.4) then
                      err2t(ii(j))=err2t(ii(j))*scale
                      err2p(ii(j))=err2p(ii(j))*scale
                   endif
                elseif(iopt.eq.3) then
                   err2t(ii(j))=err2t(ii(j))*scale
                   err2p(ii(j))=err2p(ii(j))*scale
                endif
             enddo
c
c----> average values and chisquare
              thav=0.
              phav=0.
              swt=0.
              swp=0.
              chip=0.
              chit=0.
              offset=0
              do j=2,kk
                 resp=abs(phi(ii(j))-phi(ii(1)))
                 if(resp.gt.pi) then
                  offset=1
                  resp=twopi-resp
                 endif
                 if(resp.gt.pi4) goto 100
              enddo
              if(offset.eq.1) then
                 resp=pi-phi(ii(1))
                 phi(ii(1))=pi
                 do j=2,kk
                   phi(ii(j))=phi(ii(j))+resp
                   if(phi(ii(j)).lt.0) phi(ii(j))=phi(ii(j))+twopi
                 enddo
              endif
              do j=1,4
               if(imask(j).ne.0) then
                 swt=swt+1./err2t(j)
                 swp=swp+1./err2p(j)
                 phav=phav+phi(j)/err2p(j)
                 thav=thav+theta(j)/err2t(j)
               endif
              enddo
              phav=phav/swp
              thav=thav/swt
              do j=1,4
               if(imask(j).ne.0) then
                chit=chit+(thav-theta(j))**2/err2t(j)
                chip=chip+(phav-phi(j))**2/err2p(j)
               endif
              enddo
              chisq=chit+chip
              if(offset.eq.1) then
                do j=2,4
                 if(imask(j).ne.0) phi(j)=phi(j)-resp
                enddo
                phav=phav-resp
              endif
c
              if(chisq.lt.Chi2Cut(2)) then
                 if(nhma.eq.maxhma) then
cc                   write(6,*) 'MFH123. too many matches (!).'
                   goto 9900
                 endif
                 nn=nn+1
                 nhma=nhma+1
                 xhma(1,nhma)=thav
                 xhma(2,nhma)=phav
                 xhma(3,nhma)=1./sqrt(swt)
                 xhma(4,nhma)=1./sqrt(swp)
                 xhma(5,nhma)=chit
                 xhma(6,nhma)=chip
                 j=6
                 do i=1,4
                  j=j+1
                  xhma(j,nhma)=theta(i)
                  j=j+1
                  xhma(j,nhma)=phi(i)
                 enddo
c--- kind of match
                 if(iopt.eq.1) then
                  khma(1,nhma)=11
                 elseif(iopt.eq.2) then
                  if(nctd(1).eq.0) then
                    khma(1,nhma)=14
                  else
                    khma(1,nhma)=12
                  endif
                 elseif(iopt.eq.3) then
                  if(nmip(1).eq.0) then
                    khma(1,nhma)=15
                  else
                    khma(1,nhma)=13
                  endif
                 endif
c-- keep track of pointers.
                 if(iopt.eq.1) then
                  khma(2,nhma)=i1
                  khma(3,nhma)=i2
                  khma(4,nhma)=i3
                 elseif(iopt.eq.2) then
                  khma(2,nhma)=i1
                  khma(3,nhma)=i2
                  khma(4,nhma)=0
                 elseif(iopt.eq.3) then
                  khma(2,nhma)=i1
                  khma(3,nhma)=0
                  khma(4,nhma)=i3
                 endif
              endif
c
  100       enddo
  200    enddo
  300 enddo

 9900 Continue
c
      End
c