      SUBROUTINE UXINFO(KTRIG)
C
C UNPACK THE TRIGGER MASK INTO AN ARRAY
C  INPUT 
C    BANK XTEB
C OR      DTBP ( IF NO XTEB )
C
C OUTPUT - ARRAY KTRIG FILLED WITH TRIGGER PATTERN
C----------------------------------------------
      INCLUDE '/aleph/phy/qcde.inc'
      INCLUDE '/aleph/phy/qhac.inc'
C         
      LOGICAL FIRST
      INCLUDE '/aleph/phy/qmacro.inc'
          
      DATA    FIRST/.TRUE./,NWRIT/0/,ICALL/0/
          
      IF(FIRST)THEN
            NAXTEB= NAMIND('EVEH')
            NAXTEB= NAMIND('XTEB')
            NAXTEB= NAMIND('XTRB')
            NADTBP= NAMIND('DTBP')
            FIRST=.FALSE.
      ENDIF
          
      KEVEH = IW(NAEVEH)
      KXTEB = IW(NAXTEB)
      KDTBP = IW(NADTBP)
      KXTRB = IW(NAXTRB)
          

      IF(KDTBP.NE.0)THEN
                  KPOINT = KDTBP+3
      ELSEIF(KXTEB.NE.0)THEN
                  KPOINT = KXTEB+3
      ELSEIF(KEVEH.NE.0)THEN
C
C IF ALL ELSE FAILS TAKE IT OUT OF KEVEH
C
                  KPOINT = KEVEH+7
      ELSE
                  IF(NWRIT.LT.10)
     1            WRITE(6,*)'_UXINFO_ NO VALID TRIGGER BANK FOUND !!',
     1            KRUN,KEVT
                  NWRIT = NWRIT + 1
                  KTRIG = 0
                  GOTO 999
      ENDIF
      
          
      KTRIG =  IW(KPOINT+2)
      RETURN
          
999   CONTINUE
      RETURN
      END 
      SUBROUTINE ITCTRG(NITC,FAIL)
c
c find the number of tracks according to ITC trigger bits
C
      INCLUDE '/aleph/phy/qcde.inc'
      INCLUDE '/aleph/phy/qhac.inc'
      LOGICAL FIRST,FAIL
      DATA FIRST/.TRUE./
      INCLUDE '/aleph/phy/qmacro.inc'
   

      IF(FIRST)THEN
        NAX1RG = NAMIND('X1RG')
        NAXTRB = NAMIND('XTRB')
        FIRST = .FALSE.
      ENDIF
      FAIL = .TRUE.
      NITC=0
      IF(KRUN.GT.20000)THEN
              IF    (IW(NAX1RG).NE.0) THEN
                                FAIL = .FALSE.
                                RETURN
              ENDIF
              IAA = IW(NAX1RG)+LMHLEN
              IIBIT=IW(IAA+14)
      ELSE
              IF    (IW(NAXTRB).NE.0) THEN
                                FAIL = .FALSE.
                                RETURN
              ENDIF

      ENDIF
      IF(FAIL)RETURN
*
      IF (BTEST(IIBIT,8)) NITC = NITC + 1
      IF (BTEST(IIBIT,9)) NITC = NITC + 2
      IF (BTEST(IIBIT,10)) NITC = NITC + 4
      IF (BTEST(IIBIT,11)) NITC = NITC + 8
      IF (BTEST(IIBIT,12)) NITC = NITC + 16
      RETURN
      END

      REAL FUNCTION QUTEF(WVIS,PTMAX)

      INCLUDE '/aleph/phy/qcde.inc'
      INCLUDE '/aleph/phy/qhac.inc'

      CHARACTER*80 TEXT
      CHARACTER*4 NAME,CHAINT
      REAL INEFF
      LOGICAL FIRST,DOEF

      INCLUDE '/aleph/phy/qmacro.inc'

      DATA ICOUNT/0/
      DATA FIRST/.TRUE./,DOEF/.FALSE./
      IF(FIRST)THEN
            KDOEF = IW(NAMIND('DOEF',0))
            IF(KDOEF.NE.0)DOEF=.TRUE.
            IF(.NOT.DOEF)THEN
             WRITE(TEXT,*)
     1       '_QUTEF_ EFFICIENCY CALCULATION TURNED OFF '
             CALL  QUWR(TEXT)
            ELSE
             IF(IW(KDOEF).GT.0)THEN
              NAME = CHAINT(IW(KDOEF+1))
             ENDIF
             WRITE(TEXT,*)
     1       '_QUTEF_ EFFICIENCY CALCULATION TURNED ON '
             CALL  QUWR(TEXT)
              IF((NAME.NE.'WVIS').AND.(NAME.NE.'PTMX')
     1.AND.(NAME.NE.'PTCH'))THEN
                    NAME = 'WVIS'
               ENDIF
             WRITE(TEXT,*)' - USING ',NAME
             CALL  QUWR(TEXT)
            ENDIF
            FIRST = .FALSE.
      ENDIF


      QUTEF = 1.0
      IF(.NOT.DOEF)RETURN
      IF(NAME.EQ.'WVIS')THEN
C
C RETURNS TRIGGER EFFICIENCY FOR A GIVEN WVIS
C
C THESE ARE THE SNG_C_E2 VALUES
C
C      QUTEF = 0.98 - EXP(-1.57*(WVIS-0.28))
C
C EFFICIENCY IMPROVED IN 1991 ...
C
C         IF(KRUN.LT.10000)THEN
COLD            QUTEF = 1  - 1/((WVIS+0.77)**2.6)
CNEW CUTS 1/92
C                QUTEF = 1  - 1/((WVIS+0.3)**4.4)
C         ELSE
C            QUTEF = 1  - 1/((WVIS+0.46)**4.5)
C                QUTEF = 1
C         ENDIF
          IF(KRUN.LT.10000)THEN
C
C THIS IS 1990 DATA
C
                QUTEF=0.97-1/(WVIS+1.15)**1.93
          ELSEIF(KRUN.LT.20000)THEN
C
C THIS IS 1991/92 DATA
C
             IF(WVIS.GT.4.0) THEN
                QUTEF=1.0
             ELSE
                QUTEF=0.941-1/(WVIS+0.332)**4.43
             ENDIF
          ELSEIF(KRUN.LT.20839)THEN
C
C THIS IS EARLY 1993 DATA
C
             IF(WVIS.GT.2.5) THEN
                QUTEF=1.0
             ELSE
                TEMP = 0.72*WVIS-0.49
                IF(TEMP.GT.1.0)TEMP=1.0
                QUTEF = TEMP
              ENDIF
          ELSE
C
C THIS IS MOST 1993 DATA
C
C                QUTEF=PAR(3)-1/(W+PAR(1))**PAR(2)
                 QUTEF=1.01-1/(WVIS+0.89429)**0.96101
C            QUTEF = 1  - 1/((WVIS+0.46)**4.5)
          ENDIF
      ELSEIF(NAME.EQ.'PTMX')THEN
         IF(PTMAX.GT.1.0)THEN
                QUTEF = 1.0
         ELSE
          IF(KRUN.LT.10000)THEN
                QUTEF = 1.07  - 1.0/((PTMAX+1.4)**3.0)
          ELSE
                QUTEF = 1.03  - 1/((PTMAX+1.7)**3.3)
C                QUTEF = 1
          ENDIF
        ENDIF
      ELSEIF(NAME.EQ.'PTCH')THEN
C
C CALCULATE INEFFICIENCY PRODUCT FOR EACH TRACK
C
        INEFF = 1.0      
        DO 10 ITK = KFCHT,KLCHT
           IF(XLOCK(ITK))GOTO 10
           IF(ABS(QCT(ITK)).GT.0.6)THEN
              P1 = 0.977
              P2 = 0.9963
              P3 = 3.026
           ELSE
C              IF(QPT(ITK).GT.1.0)THEN
C                 INEFF = 0.0
C                 GOTO 11
C              ELSE
                 P1 = 0.9948
                 P2 = 0.7988
                 P3 = 6.417
C                  FUNC = (0.8*QPT(ITK))**2
C                  INEFF = INEFF * (1-FUNC)
C                  GOTO 10
C              ENDIF
           ENDIF
           
           PT = QPT(ITK)
           FUNC = (PT + P2)**P3
           FUNC = 1/FUNC
           FUNC = P1 - FUNC
           IF(FUNC.LT.0.0)FUNC =0.0
           FUNC = 1.0 - FUNC
           INEFF = INEFF * FUNC
10      CONTINUE
11      CONTINUE
      QUTEF = 1.0 - INEFF
      IF(QUTEF.EQ.0.0)QUTEF = 0.0001
      TRGEFF = QUTEF
      ENDIF
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE QUXNAM
      INCLUDE '/aleph/phy/qcde.inc'
      INCLUDE '/aleph/phy/qhac.inc' 
      INTEGER   AGETDB
      CHARACTER*4 CHAINT,MNEMON(2),FULL(10)
      EXTERNAL    CHAINT,AGETDB
      INCLUDE '/aleph/phy/qmacro.inc'
C
C WRITE OUT A SUMMARY OF ALL XTBN INFORMATION
C ON THE DAF (unfortunately - its now kept on
C the data so this is less usefull!
C

C - Get bank with trigger bit names (XTBN)                              
C      ISTAT = AGETDB('XTBN',KRUN)                                       
C      IF(ISTAT.EQ.0)THEN                                                
C         CALL ALTELL('+++QUXNAM+++ Bank XTBN not found',0,'RETURN')     
C         GOTO 999                                                       
C      ENDIF                                                             
      LUN = JUNIDB(0) 
      WRITE(6,*)' DATABASE IS ON ',LUN
      NR = 999999
1     NR = NDANR(LUN,'XTBN','LE',NR-1)
      IF(NR.EQ.0)RETURN
      WRITE(6,*)'    '
      WRITE(6,*)' XTBN ',NR
      WRITE(6,*)'    '
      JXTBN = MDARD(IW,LUN,'XTBN',NR)
C      JXTBN = IW(NAMIND('XTBN'))                                        
      IF(JXTBN.NE.0)THEN                        
        WRITE(6,*)' ID   Validity Range Defined  Bit   Mnemonic ',
     1'    Name '
        WRITE(6,*)'                      Flag   Number '
                       
        DO 50 K=1,LROWS(JXTBN)                                          

         MNEMON(1) = CHAINT(ITABL(JXTBN,K,JXTBBM))
         MNEMON(2) = CHAINT(ITABL(JXTBN,K,JXTBBM+1))
         DO 49 I = 0,9
49       FULL(I+1) = CHAINT(ITABL(JXTBN,K,JXTBBN+I))

         WRITE(6,1000)ITABL(JXTBN,K,JXTBID),
     1ITABL(JXTBN,K,JXTBVR),ITABL(JXTBN,K,JXTBVR+1),
     1ITABL(JXTBN,K,JXTBTN),ITABL(JXTBN,K,JXTBTB),
     1MNEMON,FULL
1000    FORMAT(1X,I3,1X,I8,1X,I8,1X,I5,1X,I5,1X,2A4,2X,10A4)
   50   CONTINUE                                                        
      ENDIF                                                             
   55 CONTINUE                                                          
      GOTO 1
C  999 RETURN                                                           
      END                  

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C----------------------------------
      SUBROUTINE QVTX(XYZV,IER)
C         
C FIND A VERTEX...
C         
C-------------------- /QCDE/ --- COMMONs, DIMENSIONs, etc. ------------
      INCLUDE '/aleph/phy/qcde.inc'
          
C         
C Working arrays for vertex fitter
C*CA UVPWRK
      PARAMETER (MXTRA=80)
      COMMON/UVPWRK/ XYI(2,MXTRA),X0(MXTRA),Y0(MXTRA),Z0(MXTRA),
     1               RADI(MXTRA),Q(MXTRA),CPHC(MXTRA),SPHC(MXTRA),
     2               TDIP(MXTRA),SIGX(MXTRA),SIGZ(MXTRA),
     3               G(3),GG(3,3),GGINV(3,3),XS(3),XV(3),
     4               CHISQ8,CHISQO,CHISQN
C         
      DOUBLE PRECISION XYI,X0,Y0,Z0,RADI,Q,CPHC,SPHC,TDIP,SIGX,SIGZ,
     1                 G,GG,GGINV,XS,XV,CHISQ8,CHISQO,CHISQN
C         
      DIMENSION TPAR(6,40),XYZV(3),DXYZ(3),ERMXV(6)
      DIMENSION ERRD(40),ERRZ(40)
C Cut on chisq/deg of freedom for accepting vertex fit
C AJF 10. -> 100.
      DATA CHCUT/100./
      INCLUDE '/aleph/phy/qmacro.inc'
C         
          
C  Do vertex fit
C         
C         
      IER = 0
      NT = 0

      XV(1) = 0.0
      XV(2) = 0.0
      XV(3) = 0.0
      SUM = 0.0

C
C 22/3/93 - CHANGED TO USE EFLOW CHARGED TRACKS
C
      DO 80 ITK= KFEFT,KLEFT
         IF(QCH(ITK).EQ.0)GOTO 80
         IF(XLOCK(ITK))GOTO 80
         IF(QSQT(QDBS2(ITK)).EQ.0)GOTO 80
         IF(QSQT(QZBS2(ITK)).EQ.0)GOTO 80

c         IF(QP(ITK).LT.0.1)GOTO 80
         NT = NT + 1
         TPAR(1,NT) = -1./QFRFIR(ITK)
         TPAR(2,NT) = QFRFD0(ITK)*SIGN(1.,-TPAR(1,NT))
         TPAR(3,NT) = QFRFP0(ITK)
         TPAR(4,NT) = QFRFZ0(ITK)
         TPAR(5,NT) = QFRFTL(ITK)
         TPAR(6,NT) = 0.
         ERRD(NT) = QSQT(QDBS2(ITK))
         ERRZ(NT) = QSQT(QZBS2(ITK))
C
C CALC STARTING  XYZ
C
         XV(1) = XV(1) + QFRFD0(ITK)*COS(QFRFP0(ITK))*QP(ITK)
         XV(2) = XV(2) + QFRFD0(ITK)*SIN(QFRFP0(ITK))*QP(ITK)
         XV(3) = XV(3) + QFRFZ0(ITK)*QP(ITK)
         SUM = SUM + QP(ITK)
   80 CONTINUE
      IF(NT.LE.0)THEN
C            WRITE(6,*)' _QVTX : NO TRACKS TO FIT '
            IER = 1
            GOTO 999
      ENDIF

      XV(1) = XV(1) /SUM
      XV(2) = XV(2) /SUM
      XV(3) = XV(3) /SUM
C         
      NPAR = 3
      CALL UVPFIT(NT,TPAR,ERRD,ERRZ,NPAR,  NIT,NDF,CHISQ,XYZV,ERMXV)
      IF(NIT.LE.0)THEN
            IER = 2
            GOTO 999
      ENDIF
C         
C      DXYZ(1) = SQRT(1./ERMXV(1))
C      DXYZ(2) = SQRT(1./ERMXV(3))
C      DXYZ(3) = SQRT(1./ERMXV(6))
      CHPD = CHISQ/NDF
C         
      IF (CHPD.GT.CHCUT) THEN
         IER= 3
         GO TO 999
      ENDIF
          
999   RETURN
      END 
C----------------------------------
      SUBROUTINE QVTX2(XYZV,XYZE,IER)
C         
C FIND A VERTEX...
C         
C-------------------- /QCDE/ --- COMMONs, DIMENSIONs, etc. ------------
      INCLUDE '/aleph/phy/qcde.inc'
          
C         
C Working arrays for vertex fitter
C*CA UVPWRK
      PARAMETER (MXTRA=80)
      COMMON/UVPWRK/ XYI(2,MXTRA),X0(MXTRA),Y0(MXTRA),Z0(MXTRA),
     1               RADI(MXTRA),Q(MXTRA),CPHC(MXTRA),SPHC(MXTRA),
     2               TDIP(MXTRA),SIGX(MXTRA),SIGZ(MXTRA),
     3               G(3),GG(3,3),GGINV(3,3),XS(3),XV(3),
     4               CHISQ8,CHISQO,CHISQN
C         
      DOUBLE PRECISION XYI,X0,Y0,Z0,RADI,Q,CPHC,SPHC,TDIP,SIGX,SIGZ,
     1                 G,GG,GGINV,XS,XV,CHISQ8,CHISQO,CHISQN
C         
      DIMENSION TPAR(6,40),XYZV(3),DXYZ(3),ERMXV(6)
      DIMENSION ERRD(40),ERRZ(40)
      DIMENSION XYZE(3,3)
C Cut on chisq/deg of freedom for accepting vertex fit
C AJF 10. -> 100.
      DATA CHCUT/100./
      INCLUDE '/aleph/phy/qmacro.inc'
C         
          
C  Do vertex fit
C         
C         
      IER = 0
      NT = 0

      XV(1) = 0.0
      XV(2) = 0.0
      XV(3) = 0.0
      SUM = 0.0

C
C 22/3/93 - CHANGED TO USE EFLOW CHARGED TRACKS
C
      DO 80 ITK= KFEFT,KLEFT
         IF(QCH(ITK).EQ.0)GOTO 80
         IF(XLOCK(ITK))GOTO 80
         IF(QSQT(QDBS2(ITK)).EQ.0)GOTO 80
         IF(QSQT(QZBS2(ITK)).EQ.0)GOTO 80

c         IF(QP(ITK).LT.0.1)GOTO 80
         NT = NT + 1
         TPAR(1,NT) = -1./QFRFIR(ITK)
         TPAR(2,NT) = QFRFD0(ITK)*SIGN(1.,-TPAR(1,NT))
         TPAR(3,NT) = QFRFP0(ITK)
         TPAR(4,NT) = QFRFZ0(ITK)
         TPAR(5,NT) = QFRFTL(ITK)
         TPAR(6,NT) = 0.
         ERRD(NT) = QSQT(QDBS2(ITK))
         ERRZ(NT) = QSQT(QZBS2(ITK))
C
C CALC STARTING  XYZ
C
         XV(1) = XV(1) + QFRFD0(ITK)*COS(QFRFP0(ITK))*QP(ITK)
         XV(2) = XV(2) + QFRFD0(ITK)*SIN(QFRFP0(ITK))*QP(ITK)
         XV(3) = XV(3) + QFRFZ0(ITK)*QP(ITK)
         SUM = SUM + QP(ITK)
   80 CONTINUE
      IF(NT.LE.0)THEN
C            WRITE(6,*)' _QVTX2 : NO TRACKS TO FIT '
            IER = 1
            GOTO 999
      ENDIF

      XV(1) = XV(1) /SUM
      XV(2) = XV(2) /SUM
      XV(3) = XV(3) /SUM
C         
      NPAR = 3
      CALL UVPFIT(NT,TPAR,ERRD,ERRZ,NPAR,  NIT,NDF,CHISQ,XYZV,ERMXV)
      IF(NIT.LE.0)THEN
            IER = 2
            GOTO 999
      ENDIF
C         
C      DXYZ(1) = SQRT(1./ERMXV(1))
C      DXYZ(2) = SQRT(1./ERMXV(3))
C      DXYZ(3) = SQRT(1./ERMXV(6))
      CHPD = CHISQ/NDF
C         
      IF (CHPD.GT.CHCUT) THEN
         IER= 3
         GO TO 999
      ENDIF

      XYZE(1,1) = ERMXV(1)
      XYZE(1,2) = ERMXV(2)
      XYZE(2,1) = ERMXV(2)
      XYZE(2,2) = ERMXV(3)
      XYZE(1,3) = ERMXV(4)
      XYZE(3,1) = ERMXV(4)
      XYZE(2,3) = ERMXV(5)
      XYZE(3,2) = ERMXV(5)
      XYZE(3,3) = ERMXV(6)
          
999   RETURN
      END 
      SUBROUTINE UVPFIT(NTRA,TRPAR,ERRX,ERRZ,NPAR,
     &                         NIT,NDF,CHISQ,VXP,ERRVX)
C--------------------------------------------------------------------*
C! 3-D vertex fitter                                                 *
C  Author: M.E.Mermikides    27 May 1987                             *
C                                                                    *
C   INPUT:                                                           *
C           NTRA = Number of tracks to be fitted                     *
C           TRPAR(6*NTRA) = Array of parameters of each track        *
C                            (rho,d0,phi0,zo,tanl,Theta)             *
C           ERRX(NTRA) = Errors on track deviation from vertex in    *
C                        (x,y) plane. (For tracks coming from near   *
C                        the origin this corresponds closely to the  *
C                        error on d0)                                *
C           ERRZ(NTRA) = Errors on track deviation from vertex in    *
C                        z plane. (For tracks coming from near       *
C                        the origin this corresponds closely to the  *
C                        error on z0)                                *
C           NPAR = No of parameters to be fitted (2 -> xy, 3 -> xyz) *
C                                                                    *
C   OUTPUT:                                                          *
C          NIT  = Number of iterations (<0 if fit failed)            *
C          CHISQ  = chisquared of fit                                *
C          NDF = number of degrees of freedom                        *
C          VXP(NPAR) = fitted vertex coordinates (x,y(,z))           *
C          ERRVX(NDIM) = Error matrix (=Inverse covariance matrix)   *
C                 of fitted vertex position in triangular form:      *
C                (xx,xy,yy(,xz,yz.zz));    (NDIM = NPAR*(NPAR+1)/2)  *
C                                                                    *
C      CHISQ IS CALCULATED, 1ST DERIVATIVES IN G(3), 2ND IN GG(3,3)  *
C                                                                    *
C   Description:                                                     *
C   ===========                                                      *
C                                                                    *
C       The method is based on the formalism by D.H. Saxon           *
C   (NIM A234 (1985) 258-266), generalised to 3-dimensions for       *
C   use with the ALEPH tracking detectors.   The procedure is to     *
C   minimise the deviations of the tracks from the common vertex     *
C   at the coordinates (A,B,C).                                      *
C                                                                    *
C       In the x-y plane and in z, the deviations of a track         *
C   from the vertex point (A,B,C) are given by:                      *
C                                                                    *
C       d(xy) = s - rho                                              *
C                                                                    *
C       d(z) = z - C - rho*tdip*delta                                *
C                                                                    *
C       where                                                        *
C          s =sqrt[(A-xc)**2 - (B-yc)**2]                            *
C          rho = radius of curvature of track circle                 *
C          tdip = tangent of track helix dip angle                   *
C          delta = 2*arcsin(sqrt[(xt-x0)**2 + (yt-y0)**2]/(2*rho))   *
C          xt = xc +/- rho*(A-xc)/s                                  *
C          yt = yc +/- rho*(B-xc)/s                                  *
C          x0,y0 = coordinates at point of closest approach of track *
C                  to origin (d0)                                    *
C                                                                    *
C    Fitting is done using a matrix inversion method on the vector   *
C    (A,B,C)                                                         *
C                                                                    *
C   Chisq = d(xy)**2/SIGX**2 + d(z)**2/SIGZ**2                       *
C                                                                    *
C--------------------------------------------------------------------*
C*IF -DEF,DOC ---------------------------------------------
C         
C*CA UVPWRK
      PARAMETER (MXTRA=80)
      COMMON/UVPWRK/ XYI(2,MXTRA),X0(MXTRA),Y0(MXTRA),Z0(MXTRA),
     1               RADI(MXTRA),Q(MXTRA),CPHC(MXTRA),SPHC(MXTRA),
     2               TDIP(MXTRA),SIGX(MXTRA),SIGZ(MXTRA),
     3               G(3),GG(3,3),GGINV(3,3),XS(3),XV(3),
     4               CHISQ8,CHISQO,CHISQN
C         
      DOUBLE PRECISION XYI,X0,Y0,Z0,RADI,Q,CPHC,SPHC,TDIP,SIGX,SIGZ,
     1                 G,GG,GGINV,XS,XV,CHISQ8,CHISQO,CHISQN
C         
C*CC UVPWRK
      DIMENSION TRPAR(6*NTRA),ERRX(NTRA),ERRZ(NTRA),VXP(3),ERRVX(6)
C         
      DOUBLE PRECISION DRES1,DRES2,CUT(4)
      DOUBLE PRECISION D0,PHI0,PHIC, STPMAX,STEP,ST
C         
C         
      PARAMETER (PI=3.141592653589, PIBY2=PI/2.)
C         
      DATA CUT/0.03,0.01,0.,-0.03/
      DATA NITMAX/10/
      DATA STPMAX/2./
C         
      VXP(1) = 0.0001
      VXP(2) = 0.0001
      VXP(3) = 0.0001
      NDIM = NPAR*(NPAR+1)/2
      DO 50 I=1,NDIM
         ERRVX(I) = 0.
   50 CONTINUE
C         
C Starting values for vertex position (A,B,C)
C         
C      XV(1) = 0.
C      XV(2) = 0.
C      XV(3) = 0.
      NTK = MIN(NTRA,MXTRA)
C No of degrees of freedom
      NDF = (NPAR-1)*NTK - NPAR
C         
C Loop over tracks to get centre, radius and errors
C         
      DO 100 I = 1,NTK
         IND = (I-1)*6
         D0 = TRPAR(IND+2)
         PHI0 = TRPAR(IND+3)
         Q(I) = SIGN(1.,TRPAR(IND+1))
         RADI(I) = ABS(TRPAR(IND+1))
         PHIC = PHI0 - Q(I)*PIBY2
         CPHC(I) = DCOS(PHIC)
         SPHC(I) = DSIN(PHIC)
C Coordinates of circle centre (xc(i),yc(i))
         XYI(1,I) = (RADI(I) - D0)*CPHC(I)
         XYI(2,I) = (RADI(I) - D0)*SPHC(I)
C Coordinates at D0
         X0(I) = - D0*CPHC(I)
         Y0(I) = - D0*SPHC(I)
         Z0(I) = TRPAR(IND+4)
         TDIP(I) = TRPAR(IND + 5)
C Errors  
         SIGX(I) = ERRX(I)
         SIGZ(I) = ERRZ(I)
  100 CONTINUE
C         
      CHISQO = 1.D10
      NIT = 0
C         
C  Start of iteration       <-------------------------------------
  200 NIT = NIT + 1
C         
C Compute derivatives and matrices
C         
      CALL UVPMAT(NTRA,NPAR)
C         
      DO 11 J=1,3
      DO 11 K=1,3
         GGINV(K,J)= GG(K,J)
   11 CONTINUE
      CALL UMXINV(GGINV(1,1),3,3,NPAR,KFLG)
      IF (KFLG.NE.0) THEN
         NIT = -1
         GO TO 999
      ENDIF
C         
C  CHISQO = Chisq of previous iteration
C  CHISQ8 = Current value of chi squared
C  CHISQN = Updated chisq using current corrections
C         
      CHISQN = CHISQ8
      DO 410 K=1,NPAR
         XS(K) = 0.
         DO 400 L=1,NPAR
            XS(K) = XS(K) - GGINV(K,L)*G(L)
  400    CONTINUE
         CHISQN = CHISQN + 0.5*G(K)*XS(K)
  410 CONTINUE
C         
C  Test for convergence and adjust step size
C         
      DRES1 = CHISQO - CHISQ8
      DRES2 = CHISQ8 - CHISQN
      STEP = 1.
      IF(DRES2.LT.CUT(4)) STEP=0.5
C Step cut
      IF(DRES2.GE.CHISQ8) STEP = CHISQ8/DRES2
      DO 414 K = 1,NPAR
         ST = XS(K)*STEP
         IF(K.EQ.3) ST = ST * 0.25
         IF(ST.LE.0.0) ST = -ST
         IF(ST.GE.STPMAX) STEP = STEP*STPMAX/ST
  414 CONTINUE
      IF(STEP.NE.1.0) THEN
         DO 412 K = 1,3
            XS(K) = XS(K)*STEP
  412    CONTINUE
      ENDIF
C         
C Convergence test 1:  Change in chisq within given limits
C         
      IF(DRES1.LT.CUT(1).AND.DRES1.GT.CUT(4)) THEN
          GOTO 500
      ENDIF
C         
C Update parameter vector
C         
      DO 450 K = 1,NPAR
         XV(K) = XV(K) + XS(K)
  450 CONTINUE
C         
C Convergence test 2:  Chisq correction for full step below given limit
C         
      IF(DRES2.LT.CUT(2).AND.STEP.EQ.1.0) THEN
         CHISQ8 = CHISQN
         GO TO 500
      ENDIF
C         
      CHISQO = CHISQ8
C         
C Do another iteration   ------------------------------->
C         
      IF(NIT.LT.NITMAX) GOTO 200
C         
C Stop after NITMAX iterations, keep best chisq
C         
      IF(CHISQ8.LT.CUT(3).OR.CUT(3).LE.0.0) CHISQ8 = CHISQN
C         
C   Save final parameters and error matrix
C         
  500 VXP(1) = XV(1)
      VXP(2) = XV(2)
      VXP(3) = XV(3)
      CHISQ = CHISQ8
      ERRVX(1) = 0.5*GG(1,1)
      ERRVX(2) = 0.5*GG(1,2)
      ERRVX(3) = 0.5*GG(2,2)
      IF (NPAR.EQ.3) THEN
         ERRVX(4) = 0.5*GG(3,1)
         ERRVX(5) = 0.5*GG(3,2)
         ERRVX(6) = 0.5*GG(3,3)
      ENDIF
C         
  999 RETURN
      END 
C*EI      
C*HE 09/13/89 19:53:22 C
C*DK UVPMAT
      SUBROUTINE UVPMAT(NTK,NPAR)
C--------------------------------------------------------------------*
C! Auxiliary routine to UVPFIT to compute chisq and derivative       *
C  matrices                                                          *
C  Author: M.E.Mermikides    27 May 1987                             *
C--------------------------------------------------------------------*
C*IF -DEF,DOC ---------------------------------------------
C*CA UVPWRK
      PARAMETER (MXTRA=80)
      COMMON/UVPWRK/ XYI(2,MXTRA),X0(MXTRA),Y0(MXTRA),Z0(MXTRA),
     1               RADI(MXTRA),Q(MXTRA),CPHC(MXTRA),SPHC(MXTRA),
     2               TDIP(MXTRA),SIGX(MXTRA),SIGZ(MXTRA),
     3               G(3),GG(3,3),GGINV(3,3),XS(3),XV(3),
     4               CHISQ8,CHISQO,CHISQN
C         
      DOUBLE PRECISION XYI,X0,Y0,Z0,RADI,Q,CPHC,SPHC,TDIP,SIGX,SIGZ,
     1                 G,GG,GGINV,XS,XV,CHISQ8,CHISQO,CHISQN
C         
C*CC UVPWRK
C         
      DOUBLE PRECISION DDDV(3),DZDV(3),D2DDV2(3,3),D2ZDV2(3,3)
      DOUBLE PRECISION DELTA,S,SS,SSS,T,TT,DXY,DZ,
     &                 XTA,YTA,D2A,XTB,YTB,D2B,XT,YT,U0,U1,V1,
     &                 SIGX2,SIGZ2,RV2,DIR,DX,DY, DXTDA,DYTDA,
     &                 DXTDB,DYTDB,DTDA,DTDB
C         
      PARAMETER (PI=3.141592653589, PIBY2=PI/2.)
C         
      CHISQ8 = 0.
C         
      DO 10 I=1,3
        G(I)=0.0
        DO 12 J=1,3
           GG(J,I) = 0.0
   12   CONTINUE
   10 CONTINUE
C         
      RV2 = XV(1)*XV(1) + XV(2)*XV(2)
      DO 100 I = 1,NTK
         DX = XV(1) - XYI(1,I)
         DY = XV(2) - XYI(2,I)
         SS = DX*DX + DY*DY
         S = DSQRT(SS)
         SSS = SS*S
C Deviation in (x,y)
         DXY = S - RADI(I)
C Derivatives of d(xy)
         DDDV(1) = DX/S
         DDDV(2) = DY/S
         DDDV(3) = 0.
         DO 21 J=1,3
         DO 21 K=1,3
            D2DDV2(K,J) = 0.0
   21    CONTINUE
         D2DDV2(1,1) = 1./S - DX*DX/SSS
         D2DDV2(1,2) = - DX*DY/SSS
         D2DDV2(2,2) = 1./S - DY*DY/SSS
         D2DDV2(2,1) = D2DDV2(1,2)
C 2D fit  
         IF (NPAR.EQ.2) THEN
            DZ = 0.
            DO 31 J=1,3
               DZDV(J)=0.0
               DO 32 K=1,3
                  D2ZDV2(K,J) = 0.0
   32          CONTINUE
   31       CONTINUE
C 3D fit  
         ELSE
C Compute deviations in Z and derivatives.  First we find the
C coordinates at the point on the track closest to (A,B)
            XTA = XYI(1,I) + DX*RADI(I)/S
            YTA = XYI(2,I) + DY*RADI(I)/S
            D2A = ABS(XTA*XTA + YTA*YTA - RV2)
            XTB = XYI(1,I) - DX*RADI(I)/S
            YTB = XYI(2,I) - DY*RADI(I)/S
            D2B = ABS(XTB*XTB + YTB*YTB - RV2)
            IF (D2A.LT.D2B) THEN
               SOL = +1.
               XT = XTA
               YT = YTA
            ELSE
               SOL = -1.
               XT = XTB
               YT = YTB
            ENDIF
C Find z at (xt,yt); we have to take care of sign of turning
C angle from start of track to this point.
            U0 = XYI(1,I)*CPHC(I)+ XYI(2,I)*SPHC(I) - RADI(I)
            U1 = XT*CPHC(I) + YT*SPHC(I)
            V1 = YT*CPHC(I) - XT*SPHC(I)
            DIR = DSIGN(1.D0,Q(I)*V1)
            T = DSQRT((XT-X0(I))**2 + (YT-Y0(I))**2)
          
            IF(T.EQ.00000)THEN
                  WRITE(6,*)' T = 0 IN UPVMAT '
                  T = 0.000000001
            ENDIF
          
            DELTA = 2.0*DIR*DASIN( DMIN1( T/(2.*RADI(I)),1.D0 ) )
C         
            DZ = XV(3) - (Z0(I) + RADI(I)*DELTA*TDIP(I))
C         
C Derivatives of d(z), calculated in approximation of arc length by
C chord:  
C DZ = XV(3) - Z0 - DIR*TDIP(I)*SQRT((XT-X0)**2 + (YT-Y0)**2)
C         
            DXTDA = SOL*RADI(I)*(1. - DX*DX/SS) /S
            DYTDA =-SOL*RADI(I)*DX*DY/SSS
            DXTDB = DYTDA
            DYTDB = SOL*RADI(I)*(1. - DY*DY/SS) /S
            DTDA = ((XT-X0(I))*DXTDA + (YT-Y0(I))*DYTDA)/T
            DTDB = ((XT-X0(I))*DXTDB + (YT-Y0(I))*DYTDB)/T
            DZDV(1) = -DIR*TDIP(I)*DTDA
            DZDV(2) = -DIR*TDIP(I)*DTDB
            DZDV(3) = 1.
            DO 41 J=1,3
               DO 41 K=1,3
                  D2ZDV2(K,J) = 0.0
   41       CONTINUE
C Second derivatives of dz are set to zero
         ENDIF
C         
C  Sum terms of chisq.
C         
         SIGX2 = SIGX(I)**2
         SIGZ2 = SIGZ(I)**2
         CHISQ8 = CHISQ8 + DXY*DXY/SIGX2 + DZ*DZ/SIGZ2
C  Derivatives of chisq.
         DO 320 K = 1,NPAR
            G(K) = G(K) + 2.*DXY*DDDV(K)/SIGX2 + 2.*DZ*DZDV(K)/SIGZ2
            DO 325 L = 1,K
               GG(K,L) = GG(K,L) +
     &          2.*(DXY*D2DDV2(K,L)+DDDV(K)*DDDV(L))/SIGX2 +
     &          2.*( DZ*D2ZDV2(K,L)+DZDV(K)*DZDV(L))/SIGZ2
                IF (L.LT.K) GG(L,K) = GG(K,L)
  325       CONTINUE
  320    CONTINUE
C         
  100 CONTINUE
      RETURN
      END 
C*EI      
C*DK UMXINV
      SUBROUTINE UMXINV(A,L,M,N,IFAIL)
C MODIFIED BY JOE FROM CERNLIB ROUTINE SYMINV TO BE DOUBLE PRECISION
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION A(L,M),P(100),Q(100)
C      REAL*8 B(30,30),C(30,30)
      INTEGER R(100)
C         
      IFAIL=0
C      EPSILN=1.E-6
      EPSILN=1.E-14
      IF(L.LT.N) GO TO 95
      IF(M.LT.N) GO TO 96
C         
C                   CONSTRUCT TRUTH TABLE
C         
      DO 10 I=1,N
   10 R(I)=1
C         
C                   BEGIN PROGRAMME
C         
      DO 65 I=1,N
      K=0 
C         
C                   SEARCH FOR PIVOT
C         
      BIG=0.
      DO 37 J=1,N
      TEST=ABS(A(J,J))
      IF(TEST-BIG)37,37,31
   31 CONTINUE
      IBOMB=1
      IF(R(J))100,37,32
   32 BIG=TEST
      K=J 
   37 CONTINUE
C         
C                   TEST FOR ZERO MATRIX
C         
      IBOMB=2
      IF(K)100,100,38
C         
C                   TEST FOR LINEARITY
C         
   38 IF(I.EQ.1) PIVOT1=A(K,K)
      IBOMB=3
      IF(ABS(A(K,K)/PIVOT1)-EPSILN) 100,39,39
C         
C                   PREPARATION FOR ELIMINATION STEP1
C         
   39 R(K)=0
      Q(K)=1./A(K,K)
      P(K)=1.
      A(K,K)=0.0
      KP1=K+1
      KM1=K-1
      IBOMB=4
      IF(KM1)100,50,40
   40 DO 49 J=1,KM1
      P(J)=A(J,K)
      Q(J)=A(J,K)*Q(K)
      IF(R(J))100,49,42
   42 Q(J)=-Q(J)
   49 A(J,K)=0.
   50 IBOMB=5
      IF(K-N)51,60,100
   51 DO 59 J=KP1,N
      P(J)=A(K,J)
      Q(J)=-A(K,J)*Q(K)
      IBOMB=6
      IF(R(J))100,52,59
   52 P(J)=-P(J)
   59 A(K,J)=0.0
C         
C                   ELIMINATION PROPER
C         
   60 DO 65 J=1,N
      DO 65 K=J,N
   65 A(J,K)=A(J,K)+P(J)*Q(K)
C         
C                   ELEMENTS OF LEFT DIAGONAL
C         
      DO 70 J=2,N
      JM1=J-1
      DO 70 K=1,JM1
   70 A(J,K)=A(K,J)
C         
      RETURN
C         
C                   FAILURE RETURN
C         
   95 CONTINUE
   96 CONTINUE
  100 IFAIL=1
      RETURN
      END 
          
C-----------------------------------------------------------------------
      SUBROUTINE QUWR(TEXT)
C-----------------------------------------------------------------------
      CHARACTER*(*) TEXT
      WRITE(6,*)TEXT
      WRITE(2,*)TEXT
      RETURN
      END
                  
C-----------------------------------------------------------------------
      SUBROUTINE QUTHR(THRUST,THRPT,ANGLE)
C--------------------------------------------------------------------
C
C ROUTINE TO ANALYSE THRUST, RETURNS THRUST VALUE , AND PT
C
C--------------------------------------------------------------------
      INCLUDE '/aleph/phy/qcde.inc'
      INTEGER  IARRAY(200)
      REAL PVECT(3,2)
      DATA NCALL/0/
      DATA KROLD/0/,KEOLD/0/,THROLD/0.0/,TPTOLD/0.0/
      INCLUDE '/aleph/phy/qmacro.inc'
C-----------------------------------------------------------------------
C
C BOOST ALL UNLOCKED TRACKS INTO THEIR REST FRAME
C
C-----------------------------------------------------------------------
      IF(KRUN.EQ.KROLD.AND.KEVT.EQ.KEOLD)THEN
                THRUST = THROLD
                THRPT =  TPTOLD
                RETURN
      ENDIF
      NCALL = NCALL + 1
      CALL QUCMS(ICMS,IARRAY,NTRACK,0)
      IF(NTRACK.LE.0)THEN
            THRUST = 0.0
            ANGLE  = 0.0
            GOTO 999
      ENDIF
C-----------------------------------------------------------------------
C
C CALCULATE THRUST OF THESE TRACKS
C
C-----------------------------------------------------------------------
C
C
      XJTHRU = .TRUE.
      CALL QJTHRU(THRUST,'THRUST',ICMS)
      IF(THRUST.EQ.0.0)THEN
            THRUST = 0.0
            ANGLE  = 0.0
            GOTO 999
      ENDIF
      ANGLE = -1.0
C-----------------------------------------------------------------------
C
C LOOP OVER ALL THE TRACKS AND ADD UP MOMENTA IN EACH HEMISPHERE
C
C-----------------------------------------------------------------------
      ITHR = KPDIR('THRUST',ICMS)
       ANGLE = QCT(ITHR)
       THRPT = QPT(ITHR)/2.0
200   CONTINUE
 
999   CONTINUE
      CALL QVDROP(' ',ICMS)
      CALL QJOPTR('RE',' ')
      KROLD  = KRUN
      KEOLD  =  KEVT
      THROLD = THRUST 
      TPTOLD = THRPT

      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE QUTHC(THRUST,THRPT,ANGLE)
C--------------------------------------------------------------------
C
C ROUTINE TO ANALYSE THRUST of charged tracks
C  RETURNS THRUST VALUE , AND PT
C
C--------------------------------------------------------------------
      INCLUDE '/aleph/phy/qcde.inc'
      INTEGER  IARRAY(200)
      REAL PVECT(3,2)
      DATA NCALL/0/
      DATA KROLD/0/,KEOLD/0/,THROLD/0.0/,TPTOLD/0.0/
      INCLUDE '/aleph/phy/qmacro.inc'
C-----------------------------------------------------------------------
C
C BOOST ALL UNLOCKED CHarged TRACKS INTO THEIR REST FRAME
C
C-----------------------------------------------------------------------
      IF(KRUN.EQ.KROLD.AND.KEVT.EQ.KEOLD)THEN
                THRUST = THROLD
                THRPT =  TPTOLD
                RETURN
      ENDIF
      NCALL = NCALL + 1
      CALL QUCMS(ICMS,IARRAY,NTRACK,1)
      IF(NTRACK.LE.0)THEN
            THRUST = 0.0
            ANGLE  = 0.0
            GOTO 999
      ENDIF
C-----------------------------------------------------------------------
C
C CALCULATE THRUST OF THESE TRACKS
C
C-----------------------------------------------------------------------
C
C NOTE CHANGE - I AM USING KRECO NOT ICMS
C OH YES I AM ! AMY,PLUTO,CELLO USE CMS FRAM
C
C
      XJTHRU = .TRUE.
      CALL QJTHRU(THRUST,'THRUST',ICMS)
      IF(THRUST.EQ.0.0)THEN
            THRUST = 0.0
            ANGLE  = 0.0
            GOTO 999
      ENDIF
      ANGLE = -1.0
C-----------------------------------------------------------------------
C
C LOOP OVER ALL THE TRACKS AND ADD UP MOMENTA IN EACH HEMISPHERE
C
C-----------------------------------------------------------------------
      ITHR = KPDIR('THRUST',ICMS)
       ANGLE = QCT(ITHR)
       THRPT = QPT(ITHR)/2.0
200   CONTINUE
 
999   CONTINUE
      CALL QVDROP(' ',ICMS)
      CALL QJOPTR('RE',' ')
      KROLD  = KRUN
      KEOLD  =  KEVT
      THROLD = THRUST 
      TPTOLD = THRPT

      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE QUTHN(THRUST,THRPT,ANGLE)
C--------------------------------------------------------------------
C
C ROUTINE TO ANALYSE THRUST, RETURNS THRUST VALUE , AND PT
C
C--------------------------------------------------------------------
      INCLUDE '/aleph/phy/qcde.inc'
      INTEGER  IARRAY(200)
      REAL PVECT(3,2)
      DATA NCALL/0/
      DATA KROLD/0/,KEOLD/0/,THROLD/0.0/,TPTOLD/0.0/
      INCLUDE '/aleph/phy/qmacro.inc'
C-----------------------------------------------------------------------
C
C BOOST ALL UNLOCKED  GAT PHOTONS INTO THEIR REST FRAME
C
C-----------------------------------------------------------------------
      IF(KRUN.EQ.KROLD.AND.KEVT.EQ.KEOLD)THEN
                THRUST = THROLD
                THRPT =  TPTOLD
                ANGLE = ANGOLD
                RETURN
      ENDIF
      NCALL = NCALL + 1
      CALL QUCMS(ICMS,IARRAY,NTRACK,2)
      IF(NTRACK.LE.0)THEN
            THRUST = 0.0
            ANGLE  = 0.0
            GOTO 999
      ENDIF
C-----------------------------------------------------------------------
C
C CALCULATE THRUST OF THESE TRACKS
C
C-----------------------------------------------------------------------
C
C NOTE CHANGE - I AM USING KRECO NOT ICMS
C OH YES I AM ! AMY,PLUTO,CELLO USE CMS FRAM
C
C
      XJTHRU = .TRUE.
      CALL QJTHRU(THRUST,'THRUST',ICMS)
      ANGLE = -1.1
      IF(THRUST.EQ.0.0)THEN
            ANGLE  = -1.1
            GOTO 999
      ENDIF
      ITHR = KPDIR('THRUST',ICMS)
       ANGLE = QCT(ITHR)
       THRPT = QPT(ITHR)/2.0
200   CONTINUE

 
999   CONTINUE
      CALL QVDROP(' ',ICMS)
      CALL QJOPTR('RE',' ')

      KROLD  = KRUN
      KEOLD  =  KEVT
      THROLD = THRUST 
      TPTOLD = THRPT
      ANGOLD = ANGLE 

      RETURN
      END
C
      SUBROUTINE QUCMS(ICMS,IARRAY,III,IFLAG)
C
      INCLUDE '/aleph/phy/qcde.inc'
C
      INTEGER IARRAY(*)
      CHARACTER*6 NAME(3)
 
      LOGICAL FIRST
      LOGICAL TALL
      DATA    TALL/.FALSE./
      DATA    FIRST/.TRUE./
      DATA    NAME/'CMSALL','CMSCHT','CMSGAT'/
      INCLUDE '/aleph/phy/qmacro.inc'
C
C----------------- CODE STARTS HERE ------------------------------------
C
C CREATE A 'CMS' TRACK
C
       IF(IFLAG.EQ.0)THEN
             KFIRST = KFEFT
             KLAST  = KLEFT
       ELSEIF(IFLAG.EQ.1)THEN
             KFIRST = KFCHT
             KLAST  = KLCHT
       ELSEIF(IFLAG.EQ.2)THEN
             KFIRST = KFGAT
             KLAST  = KLGAT
       ENDIF

       ICMS = KVNEW(DUMMY)
       IUSED = 0

       DO  ITK = KFIRST,KLAST
         IF (.NOT.XLOCK(ITK))THEN
          IUSED = IUSED + 1
          CALL QVADDN(ICMS,ITK)
        ENDIF
       ENDDO
      IF(IUSED.LE.1)THEN
        III = 0
        IARRAY(III+1) = 0
        RETURN
      ENDIF
C
C NOW BOOST TRACKS INTO THIS FRAME
C AND MAKE A NAMED COPY
C
      III = 0
C
      DO 80 ITK = KFIRST,KLAST
        IF (.NOT.(XLOCK(ITK)))THEN
            III = III + 1
            ITEMP =  KTLOR(ITK,ICMS)
            IARRAY(III) = KVSAVE(ITEMP,NAME(IFLAG))
        ENDIF
80    ENDDO
C
      CALL QJOPTR('NO',NAME(IFLAG))
C
      IARRAY(III+1) = 0
      RETURN
      END
 
      SUBROUTINE MCANAL(WTRUE)
      INCLUDE '/aleph/phy/qcde.inc'
      CHARACTER*2 NAME
      DATA IFAIL/0/,IDEB/0/,ICALL/0/
      INCLUDE '/aleph/phy/qmacro.inc'
C
C
C CALC W TRUE FROM SCATTERED ELECTRONS
C
C SCATTERED ELECTRONS IDENTIFIED AS FIRST TWO ELECTRONS LEAVING PRIMARY VERTEX
C ( only makes sense for a gamma gamma generator !)
C
 
      IDEB = IDEB + 1
      WTRUE = 0
      IEL1 = 0
      EE = 0.0
      IEL2 = 0
      NAME = 'E+'
1     ITK = KPDIR(NAME,KMONTE)
90    IF(ITK.EQ.0)GOTO 100
C
C CHECK THE VERTEX, AND ENERGY
C
        IF(KORIV(ITK).NE.KFMCV)GOTO 95
        IF(QE(ITK).LT.EE)GOTO 95
              IF(NAME.EQ.'E+')IEL1 = ITK
              IF(NAME.EQ.'E-')IEL2 = ITK
95    ITK = KFOLLO(ITK)
      GOTO 90
100   CONTINUE
      IF(NAME.EQ.'E+')THEN
           NAME = 'E-'
           EE = 0.0
           GOTO 1
      ENDIF
 
 
      IF(IEL1.EQ.0.AND.IEL2.EQ.0)THEN
           IFAIL = IFAIL + 1
           IF(IFAIL.LT.10)
     1     WRITE(6,*)' SCATTERED ELECTRONS NOT FOUND '
      ELSE
            PTOTS = (QX(IEL1)+QX(IEL2))**2+
     1              (QY(IEL1)+QY(IEL2))**2+
     1              (QZ(IEL1)+QZ(IEL2))**2
            ETOTM  = QE(IEL1)+QE(IEL2)
            WTRUE = SQRT((QELEP-ETOTM)**2 - PTOTS )
c            IF(IDEB.LT.10)THEN
c                 WRITE(6,*)'_MCANAL_ ELECTRONS ... '
c                 CALL QWITK(IEL1,'  ')
c                 CALL QWITK(IEL2,'  ')
c            ENDIF
 
      ENDIF
      END
