Subroutine MFTRF
c   ====================
c
c This is an interface routine to fill /FWHMU/ starting from the
c /MFRTZ/ tracks. /FWHMU/ is the reference common for the tracks
c to be matched.
c
      Implicit None
c
#include "partap.inc"
#include "mfrtz.inc"
#include "fwhmu.inc"
#include "mfsigm.inc"
#include "mfflg.inc"
c
#include "mfct.inc"
c
      integer itr
      integer i
      real cz,x2,y2,xy,s2,xf(3),pf(3),angf(2),pmu,cms,cv(15),
     .     dz,dz2,errp,rd
c
      real zxtfm
      data zxtfm /450./
      save zxtfm
c
      real twopi
      save twopi
c
      integer ierr
c
      logical first/.true./
      save    first
c
      if(first) then
       first=.false.
       twopi=2.*acos(-1.)
       rd=360./twopi
      endif
c
      ntmf=0
c
      do 10 i=1,coutab(mfrtz)
       call fettab(mfrtz,id,i)
       if(mfrtz_status.lt.0) goto 10
       if(mfrtz_x(6).gt.550.) then
        Call MFXT(ZXTFM,XF,PF,ANGF,IERR)
        if(ierr.gt.0) goto 10
        ntmf=ntmf+1
        if(ntmf.gt.mxtm) then
          ntmf=mxtm
cc          write(6,*) '... too many fmu tracks. array cut... '
          goto 9900
        endif
        XTMF(1,ntmf)=xf(1)
        XTMF(2,ntmf)=xf(2)
        XTMF(3,ntmf)=xf(3)
        XTMF(4,ntmf)=angf(1)
        XTMF(5,ntmf)=angf(2)
        errp=mfrtz_cov(15)/mfrtz_x(5)/mfrtz_x(5)
        pmu=sqrt(pf(1)**2+pf(2)**2+pf(3)**2)
        XTMF(6,ntmf)=sign(1.,mfrtz_x(5))*pmu
        XTMF(7,ntmf)=mfrtz_octant
        XTMF(8,ntmf)=mfrtz_id
        XTMF(9,ntmf)=mfrtz_status
c---> approximate: consider linear error prop and then scale.
        DZ=mfrtz_x(6)-zxtfm
        dz2=dz*dz
        call ucopy(mfrtz_cov,cv,15)
        cv(1)=cv(1)+2.*cv(4)*dz+cv(6)*dz2
        cv(2)=cv(2)+(cv(5)+cv(7))*dz+cv(9)*dz2
        cv(3)=cv(3)+2.*cv(8)*dz+cv(10)*dz2
        cv(4)=cv(4)+cv(6)*dz
        cv(5)=cv(5)+cv(9)*dz
        cv(7)=cv(7)+cv(9)*dz
        cv(8)=cv(8)+cv(10)*dz
        cv(11)=cv(11)+cv(13)*dz
        cv(12)=cv(12)+cv(14)*dz
c--> mulsca contribution: assume 1 m of iron crossed.
        cms=(0.013*sqrt(pmu**2+0.011)/pmu/pmu)**2
        cz=pf(3)/pmu
        x2=(pf(1)/pf(3))**2
        y2=(pf(2)/pf(3))**2
        s2=x2+y2
        xy=pf(1)*pf(2)/pf(3)/pf(3)
        ETMF(1,ntmf)=cv(1)*scalefm
        ETMF(2,ntmf)=cv(3)*scalefm
        ETMF(3,ntmf)=(4.*cz**2*(x2*cv(6)+y2*cv(10)+
     .               2.*xy*cv(9))/s2)*scalefm + cms
        ETMF(4,ntmf)=((y2*cv(6)+x2*cv(10)-2.*xy*cv(9))/s2/s2)*scalefm +
     .               cms/(1.-cz**2)
        ETMF(5,ntmf)=errp*XTMF(6,ntmf)**2
       else
        ntmf=ntmf+1
        if(ntmf.gt.mxtm) then
          ntmf=mxtm
cc          write(6,*) '... too many fmu tracks. array cut... '
          goto 9900
        endif
        XTMF(1,ntmf)=mfrtz_x(1)
        XTMF(2,ntmf)=mfrtz_x(2)
        XTMF(3,ntmf)=mfrtz_x(6)
        x2=mfrtz_x(3)**2
        y2=mfrtz_x(4)**2
        xy=mfrtz_x(3)*mfrtz_x(4)
        s2=x2+y2
        cz=1./sqrt(1.+s2)
        XTMF(4,ntmf)=acos(cz)
        XTMF(5,ntmf)=atan2(mfrtz_x(4),mfrtz_x(3))
        if(XTMF(5,ntmf).lt.0.) XTMF(5,ntmf)=XTMF(5,ntmf)+twopi
        XTMF(6,ntmf)=1./mfrtz_x(5)
        XTMF(7,ntmf)=mfrtz_octant
        XTMF(8,ntmf)=mfrtz_id
        XTMF(9,ntmf)=mfrtz_status
        ETMF(1,ntmf)=mfrtz_cov(1)
        ETMF(2,ntmf)=mfrtz_cov(3)
        ETMF(3,ntmf)=4.*cz**2*(x2*mfrtz_cov(6)+y2*mfrtz_cov(10)+
     .               2.*xy*mfrtz_cov(9))/s2
        ETMF(4,ntmf)=(y2*mfrtz_cov(6)+x2*mfrtz_cov(10)-
     .               2.*xy*mfrtz_cov(9))/s2/s2
        ETMF(5,ntmf)=mfrtz_cov(15)
       endif
   10 enddo
c
      if(Mdebug.gt.0) then
        write(6,*)
        write(6,*) ' Number of FMU tracks: ',ntmf
        write(6,*) ' ===================== '
        do i=1,ntmf
         write(6,*) 'Track ',i,' Theta/Phi ',XTMF(4,i)*rd,XTMF(5,i)*rd
        enddo
        write(6,*)
      endif
c
 9900 continue
      end
c