      Subroutine MFFASE(Ntr,Ntrp,Recp,Angr,Match)
***********************************************************
*                                                         *
*            FMUON Fast Selection routine                 *
*          --------------------------------               *
*                                                         *
* Input: None                                             *
*                                                         *
* Output:                                                 *
*         Ntr(INTEGER):  number of octants with tracks    *
*         Ntrp(INTEGER): number of octants with           *
*                        pointing tracks                  *
*         Recp(REAL):    maximum reconstructed momentum   *
*                        of pointing tracks               *
*         Angr(INTEGER): theta region flag                *
*         Match(REAL):   momentum.flag                    *
*                        momentum = max(FMU,CTD)          *
*                        flag = 1 --> match FMU-CTD-MIP   *
*                        flag = 2 --> match FMU-MIP       *
*                                                         *
***********************************************************
      Implicit None
*
#include "mfrtz.inc"
#include "partap.inc"
*
      Integer Ntr,Ntrp,Angr
      Integer Ntrq
      Real Recp,Match
      Real Angl1,Angl2
      Parameter (Angl1=.13963,Angl2=.27925)
      Integer Nmf,i,coct(8),coctp(8),Thr(3),Lt1,Lt2,Lt3,Lt4,Lt5,Lt6
      Integer Dc1,Dc2,Dc3,Dc4,Nlt,Ndc,Nhit
      Real rho,theta,X0,Y0,rhov,mfzextr
*
* match MAMMA
      call FMUFLAG(Match)
*
      Nmf=Coutab(Mfrtz)
      If(Nmf.eq.0)Return
      Ntr=0
      Ntrp=0
      Recp=0.
      Angr=0
      Call vzero(coct,8)
      Call vzero(coctp,8)
      Call vzero(Thr,3)
* Track loop
      Do i=1,Nmf
       Call fettab(Mfrtz,id,i)
* octant
       coct(mfrtz_octant)=1
* Track lenght
       Lt1=mod(mfrtz_hitlt,2)
       Lt2=mod(mfrtz_hitlt/10,2)
       Lt3=mod(mfrtz_hitlt/100,2)
       Lt4=mod(mfrtz_hitlt/1000,2)
       Lt5=mod(mfrtz_hitlt/10000,2)
       Lt6=mod(mfrtz_hitlt/100000,2)
       Nlt=Lt1+Lt2+Lt3+Lt4+Lt5+Lt6
       Dc1=mod(mfrtz_hitdc,2)
       Dc2=mod(mfrtz_hitdc/10,2)
       Dc3=mod(mfrtz_hitdc/100,2)
       Dc4=mod(mfrtz_hitdc/1000,2)
       Ndc=Dc1+Dc2+Dc3+Dc4
* Rhov
       mfzextr=mfrtz_X(6)-(mfrtz_X(1)*mfrtz_X(3)+mfrtz_X(2)*mfrtz_X(4))/
     .                    (mfrtz_X(3)**2+mfrtz_X(4)**2)
       Rhov=sqrt((mfrtz_X(1)+mfrtz_X(3)*(mfzextr-mfrtz_X(6)))**2+
     .           (mfrtz_X(2)+mfrtz_X(4)*(mfzextr-mfrtz_X(6)))**2)
C       X0=mfrtz_x(1)-(mfrtz_x(6)*mfrtz_x(3))
C       Y0=mfrtz_x(2)-(mfrtz_x(6)*mfrtz_x(4))
C       Rhov=sqrt(X0*X0+Y0*Y0)
* Define pointing track
       Ntrq=Lt5+Lt6+Dc3+Dc4
       if(Rhov.lt.50.and.Ntrq.gt.0)Then
        coctp(mfrtz_octant)=1
        If(abs(1./mfrtz_x(5)).gt.Recp)Recp=abs(1./mfrtz_x(5))
       Endif
* theta range
       rho=sqrt(mfrtz_x(1)*mfrtz_x(1)+mfrtz_x(2)*mfrtz_x(2))
       theta=atan2(rho,mfrtz_x(6))
       if(theta.lt.Angl1)then
        thr(1)=1
       elseif(theta.ge.Angl1.and.theta.lt.Angl2)then
        thr(2)=1
       elseif(theta.ge.Angl2)then
        thr(3)=1
       endif
      Enddo
* Octant with track(s)
      Ntr=Coct(1)+Coct(2)+Coct(3)+Coct(4)+Coct(5)+Coct(6)
     >+Coct(7)+Coct(8)
* Octant with pointing track(s)
      Ntrp=Coctp(1)+Coctp(2)+Coctp(3)+Coctp(4)+Coctp(5)+Coctp(6)
     >+Coctp(7)+Coctp(8)     
      Angr=Thr(1)+10*Thr(2)+100*Thr(3)
      Return
      End 
C-----------------------------------------------------------------------

      Subroutine FMUFLAG(result)
c
      implicit none
c
#include "mfmtch.inc"
#include "fwhmu.inc"
#include "fwtrk.inc"
c
      real result
      integer iok
      integer n1,n2
      real p1(maxhma),p2(maxhma),pmax
      integer i,i1,i2,i3,jz,ind(maxhma)
c
      result=0.
      n1=0
      n2=0
c
      Call MAMMA(2,-1.,-1.,iok)
      if(iok.lt.1) return
c
      do i=1,nhma
       if(khma(1,i).eq.1) then
        if(n1.lt.maxhma) then
         n1=n1+1
         i1=khma(2,i)
         jz=mod(khma(4,i),10)
         i3=khma(4,i)/1000
         p1(n1)=max(abs(xtmf(6,i1)),pctdh(i3,jz))
        endif
       elseif(khma(1,i).eq.2.or.khma(1,i).eq.4) then
        if(n2.lt.maxhma) then
         n2=n2+1
         i1=khma(2,i)
         i2=khma(3,i)/1000
         jz=mod(khma(3,i),10)
         p2(n2)=abs(xtmf(6,i1))
        endif
       endif
      enddo
c
      if(n1+n2.lt.1) goto 9900
c
      if(n1.gt.0) then
       if(n1.gt.1) then
        call sortzv(p1,ind,n1,1,1,0)
        i1=ind(1)
       else
        i1=1
       endif
       pmax=p1(i1)
       result=float(int(pmax+0.5))+0.1
      elseif(n2.gt.0) then
       if(n2.gt.1) then
        call sortzv(p2,ind,n2,1,1,0)
        i1=ind(1)
       else
        i1=1
       endif
       pmax=p2(i1)
       result=float(int(pmax+0.5))+0.2
      endif
c
 9900 continue
c
      end
C-----------------------------------------------------------------------

