Subroutine MFTK123(jzone,iopt,nn)
c    ===================================
c
      Implicit None
c
c
c  JZONE [I]  = 1/2/3    <---> F/B/R
c  IOPT  [I]  = 1/2/3/4  <---> mu-mip-ctd/mu-mip/mu-ctd/mip-ctd
c  NN    [O]  = # of matches
c
c
c
#include "partap.inc"
#include "fwmip.inc"
#include "fwtrk.inc"
#include "fwhmu.inc"
#include "mfsigm.inc"
#include "mfmtch.inc"
#include "mfcuts.inc"
#include "zdskey.inc"
#include "fmckin.inc"
#include "mfflg.inc"
c
      integer iopt,nn,jzone
      integer i1,i2,i3,j,i,ithz
c
      real theta(5),phi(5),err2t(5),err2p(5),x1,y1,z1,x2,y2,z2,x3,y3,
     .     z3,x12,y12,z12,x13,y13,z13,x23,y23,z23,r12,r13,r23,th,ph,
     .     cx12,cy12,cz12,cx13,cy13,cz13,cx23,cy23,cz23,aaa,s2x,s2y,
     .     r3sq,cst,snt,csp,snp,sig2xctd,sig2yctd,pmu2,cms,cms12,cms23,
     .     cms13,thav,phav,swt,swp,chip,chit,chisq,resp
c
      real pi,twopi,pi4
      save pi,twopi,pi4
c
      real prob
c
      integer nloop(3),imask(5),kk,ii(5)
c
      integer iucomp,kmu
c
      integer offset
      logical first /.true./
      save    first
c
      real delta_cut
      data delta_cut /0.2/
      save delta_cut
c
      real sigt(4,5,4),sigp(4,5,4)
      data sigt/
c iopt=1
     &          4*0.007,
     &          4*0.100,
     &          0.020,0.012,0.010,0.010,
     &          0.130,0.050,0.035,0.025,
     &          0.070,0.025,0.014,0.005,
c iopt=2
     &          0.006,0.006,0.006,0.006,
     &          0.012,0.012,0.012,0.012,
     &          4*0.,
     &          4*0.,
     &          4*0.,
c iopt=3
     &          0.006,0.006,0.006,0.006,
     &          4*0.,
     &          0.022,0.015,0.010,0.008,
     &          4*0.,
     &          0.090,0.040,0.016,0.004,
c iopt=4
     &          4*0.,4*0.,4*0.,
     &          0.130,0.050,0.035,0.025,
     &          0.070,0.025,0.014,0.005/
c
      data sigp/
c iopt=1
     &          0.020,0.020,0.020,0.020,
     &          0.075,0.075,0.075,0.075,
     &          0.090,0.085,0.070,0.020,
     &          0.130,0.130,0.120,0.050,
     &          0.120,0.090,0.075,0.025,
c iopt=2
     &          0.100,0.100,0.100,0.100,
     &          0.100,0.100,0.100,0.100,
     &          4*0.,
     &          4*0.,
     &          4*0.,
c iopt=3
     &          0.015,0.015,0.015,0.015,
     &          4*0.,
     &          0.060,0.035,0.020,0.015,
     &          4*0.,
     &          0.070,0.040,0.020,0.015,
c iopt=4
     &          4*0.,4*0.,4*0.,
     &          0.130,0.130,0.120,0.050,
     &          0.120,0.090,0.075,0.025/
c
c
      real scale
      data scale /2./
c
      real probt,probp,probtp
c
      save sigt,sigp,scale
c
      if(first) then
       first=.false.
       pi=acos(-1.)
       twopi=2.*pi
       pi4=pi/4.
      endif
c
      call vzero(theta,5)
      call vzero(phi,5)
c
      nn=0
c
c iopt = 1  mu_track -> mip -> ctd
c        2  mu_track -> mip
c        3  mu_track        -> ctd
c        4              mip -> ctd
c
      if(iopt.eq.1) then
        nloop(1)=ntmf
        nloop(2)=nmip(jzone)
        nloop(3)=nctd(jzone)
        kk=5
        do j=1,5
         imask(j)=1
         ii(j)=j
        enddo
      elseif(iopt.eq.2) then
        nloop(1)=ntmf
        nloop(2)=nmip(jzone)
        nloop(3)=1
        kk=2
        ii(1)=1
        ii(2)=2
        imask(1)=1
        imask(2)=1
        imask(3)=0
        imask(4)=0
        imask(5)=0
      elseif(iopt.eq.3) then
        nloop(1)=ntmf
        nloop(2)=1
        nloop(3)=nctd(jzone)
        kk=3
        ii(1)=1
        ii(2)=3
        ii(3)=5
        imask(1)=1
        imask(2)=0
        imask(3)=1
        imask(4)=0
        imask(5)=1
      elseif(iopt.eq.4) then
        nloop(1)=1
        nloop(2)=nmip(jzone)
        nloop(3)=nctd(jzone)
        kk=2
        ii(1)=4
        ii(2)=5
        imask(1)=0
        imask(2)=0
        imask(3)=0
        imask(4)=1
        imask(5)=1
      else
        write(6,*) 'MFTK123. IOPT error!'
        return
      endif
c
c---> loop over mu elements
      do i1=1,nloop(1)
       if(ntmf.gt.0.and.imask(1).gt.0) then
        theta(1)=XTMF(4,i1)
        phi(1)=XTMF(5,i1)
        x1=XTMF(1,i1)
        y1=XTMF(2,i1)
        z1=XTMF(3,i1)
c
       endif
c-----> loop over mips elements
        do i2=1,nloop(2)
         if(nmip(jzone).gt.0.and.iopt.ne.3) then
          x2=xmip(i2,jzone)
          y2=ymip(i2,jzone)
          z2=zmip(i2,jzone)
          if(imask(2).gt.0) then
           x12=x1-x2
           y12=y1-y2
           z12=z1-z2
           Call mfCaPo(x12,y12,z12,r12,th,ph)
           theta(2)=th
           phi(2)=ph
          endif
         endif
c-------> loop over ctd elements
         do i3=1,nloop(3)
          if(nctd(jzone).gt.0.and.iopt.ne.2) then
           x3=xctd(i3,jzone) 
           y3=yctd(i3,jzone)
           z3=zctd(i3,jzone)
           theta(5)=thctd(i3,jzone)
           phi(5)=psctd(i3,jzone)
           if(imask(3).gt.0) then
            x13=x1-x3
            y13=y1-y3
            z13=z1-z3
            Call mfCaPo(x13,y13,z13,r13,th,ph)
            theta(3)=th
            phi(3)=ph
           endif
           if(imask(4).gt.0) then
            x23=x2-x3
            y23=y2-y3
            z23=z2-z3
            Call mfCaPo(x23,y23,z23,r23,th,ph)
            theta(4)=th
            phi(4)=ph
           endif
c
          endif
c
          if(iopt.eq.1.or.iopt.eq.3.or.iopt.eq.4) then
             if(ndfctd(i3,jzone).lt.10) then
                ithz=1
             elseif(ndfctd(i3,jzone).lt.20) then
                ithz=2
             elseif(ndfctd(i3,jzone).lt.30) then
                ithz=3
             else
                ithz=4
             endif
          else
             ithz=1
          endif
c
          do j=1,kk
             err2t(ii(j))=sigt(ithz,ii(j),iopt)**2
             err2p(ii(j))=sigp(ithz,ii(j),iopt)**2
             if(ii(j).eq.3.or.ii(j).eq.4.or.ii(j).eq.5) then
                err2t(ii(j))=err2t(ii(j))*scale
                err2p(ii(j))=err2p(ii(j))*scale
             endif
          enddo
c
c--------> chisquares
           phav=0.
           thav=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,5
            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,5
            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.gt.0) then
            do j=2,5
             if(imask(j).ne.0) phi(j)=phi(j)-resp
            enddo
            phav=phav-resp
           endif
           if(chisq.lt.Chi2Cut(1)) then
             if(nhma.eq.maxhma) then
cc               write(6,*) 'MFTK123. 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,5
              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)=1
              khma(2,nhma)=i1
              khma(3,nhma)=i2*1000+jzone
              khma(4,nhma)=i3*1000+jzone
              ictdused(i3,jzone)=1
              imipused(i2,jzone)=1
             elseif(iopt.eq.2) then
              khma(2,nhma)=i1
              khma(3,nhma)=i2*1000+jzone
              khma(4,nhma)=0
              if(nctd(jzone).eq.0) then
               khma(1,nhma)=4
              else
               khma(1,nhma)=2
              endif
              imipused(i2,jzone)=1
             elseif(iopt.eq.3) then
              khma(2,nhma)=i1
              khma(3,nhma)=0
              khma(4,nhma)=i3*1000+jzone
              if(nmip(jzone).eq.0) then
               khma(1,nhma)=5
              else
               khma(1,nhma)=3
              endif
              ictdused(i3,jzone)=1
             elseif(iopt.eq.4) then
              khma(2,nhma)=0
              khma(3,nhma)=i2*1000+jzone
              khma(4,nhma)=i3*1000+jzone
              khma(1,nhma)=6
              ictdused(i3,jzone)=1
              imipused(i2,jzone)=1
             endif
           endif
  100     enddo
  200   enddo
  300 enddo
c
 9900 continue
      end
c