C   08/06/83 308081851  MEMBER NAME  MUFFLX   (JADEMUS)     FORTRAN
C
C LAST CHANGE 03.11 14/07/83 HUGH MCCANN - VARIOUS CHANGES.
C      CHANGE 12.17 05/05/83 CHRIS BOWDERY- GENERAL TIDY-UP. NEW HAMB=2200000400
C      CHANGE 22.04 27/04/83 CHRIS BOWDERY- FIX BUG WITH IAPPRO+BAD HIT
C      CHANGE 19.10 16/03/82 CHRIS BOWDERY- FOR 'LET-OFF' CHECK (LOCH)
C      CHANGE 14.30 04/02/82 HUGH MCCANN  - CORRECT TREATMENT OF CASE
C                                           WHERE ALL HITS ARE BAD.
C      CHANGE 11.30 18/10/81 HUGH MCCANN  - CORRECT MUR2/3 COCK-UP.
C      CHANGE 10.05 15/10/81 HUGH MCCANN  - PROTECT LAYRAB (<= 10).
C      CHANGE 08.30 12/10/81 HUGH MCCANN  - FOR OVERLAP HITS.
C
C-----------------------------------------------------------------------
      SUBROUTINE MUFFLX(HAMB)
C-----------------------------------------------------------------------
C
C        THIS SUBROUTINE CALCULATES CHI-SQUARED AND SELECTS THE BEST
C        HIT COMBINATION.
C
C-----------------------------------------------------------------------
C
      IMPLICIT INTEGER*2 (H)
C
C             COMMONS AND ARRAY DECLARATIONS
C
C   20/10/81 202041521  MEMBER NAME  CMUFWORK (JADEMUS)     FORTRAN
C   12/10/81 110152138  MEMBER NAME  CMUFWORK (JADEMUS1)    FORTRAN
C-------------START OF MACRO CMUFWORK-----------------------------------
C LAST CHANGED 21.30 15/10/81 HUGH MCCANN  - EXTEND IOVLAP TO ALLOW
C                                            > 10 ASSOC HITS IN MUFFLY.
C LAST CHANGED 10.00 04/10/81 HUGH MCCANN  - FOR USE OF OVERLAP HITS.
C LAST CHANGED 14.46 17/06/81 JOHN ALLISON - TO ADD 'OVER' VARIABLES.
C LAST CHANGED 23.00 09/04/81 HUGH MCCANN - JADEMUS UPDATE.
C
      COMMON /CWORK/X,Y,Z,DCX,DCY,DCZ,P,PSQ,E,ESQ,D,GM,AB,RD,DE,PPIDK,
     *PKDK,PPBPEN,PPPEN,PPIPEN,PKPEN,DUMW1(29),X8,Y8,Z8,DCX8,DCY8,
     *DCZ8,P8,PSQ8,E8,ESQ8,D8,GM8,AB8,RD8,DE8,PPIDK8,PKDK8,PPBPE8,
     *PPPEN8,PPIPE8,PKPEN8,DUMW28(29),X0,Y0,Z0,DCX0,DCY0,DCZ0,P0,PSQ0,
     *E0,ESQ0,D0,GM0,AB0,RD0,DE0,PPIDK0,PKDK0,PPBPE0,PPPEN0,PPIPE0,
     *PKPEN0,DUMW2(29),X9,Y9,Z9,DCX9,DCY9,DCZ9,P9,PSQ9,E9,ESQ9,D9,GM9,
     *AB9,RD9,DE9,PPIDK9,PKDK9,PPBPE9,PPPEN9,PPIPE9,PKPEN9,DUMW29(29),
C-----200 WORDS UP TO HERE----------------------------------------------
     *X1,Y1,Z1,X2,Y2,Z2,XMID,YMID,ZMID,X3,Y3,Z3,D1,D2,D3,DTC,DUMW3(4),
     *DSTEP,GMSTEP,ABSTEP,RDSTEP,DESTEP,DUMW4(15),ADCZ,COSEC,STPINI,CURV
     *,DANG,DX,DY,DZ,R,RSQ,PT,DPINA,DKNA,DUMW5(7),VX,VY,VZ,RMSXY,RMSZ,
     *VTXYD,VTZD,VTXYAN,VTZANG,VMSANG,VMSD,CMS,DUMW6(8),SDS(10),DUMW7(1000002000
C
C    *SDX,SDY,SDZ,SDTXYD,SDTZD,SDTXYA,SDTZAN,SDMSAN,SDMSD,SDCMS,DM7(10),
C
C NOTE ABOVE CARD SHORTENED TO SDS(10),DUMW7(10), SO THAT STATEMENT DOES
C   NOT EXCEED 19 CONTINUATION CARDS.  I THINK SD--'S ARE NOT USED.
     *),CHISQT,NDFT,CHIPRT,CHISTL,NDFTL,CHIPTL,NBADCH,DUMW8(13),IDTRK,
     *NCEPTS,NTHIS,INEFF,NHLAYR,INFLAG,NELIPS,NGLAYR,DUMW9(72),IREG(200)
C-----600 WORDS UP TO HERE----------------------------------------------
     *,RDOTD(200),CINT(3,200),NFLAG(200),ITEM(200),EM(30,30),DALONG(10),
C-----2710 WORDS UP TO HERE---------------------------------------------
     *DEV(30),VXYDA(10),VZDA(10),CXYA(10),CZA(10),INCUT(30),BADA(30),G(200003100
     *0,20),DEVA(20),IPERM(10),MWORK(20),ITHISA(10),IJA(10),EM1(20,20),E
     *MG(20,20),ELIPSE(700),INCHAM(10),INREG(10),IAPPRO(30),HBADCH(1000)
     *,CTDASH(10),CLDASH(10),VTDASH(10),VLDASH(10),IHTREG(20),LAYRAB,IEF
C-----5421 WORDS UP TO HERE---------------------------------------------
     *FRG(10),DUMW10(69),AAA(10),BBB(10),CCC(10),DDD(10),DOVER,GMOVER,AB
     *OVER,RDOVER,DEOVER,DUMW11(55),X7(50),IHTP(30),IOVLAP(20,20)
C-----6080 WORDS UP TO HERE---------------------------------------------
      DIMENSION ILIPSE(700)
      INTEGER*2 HLIPSE(1400)
      DIMENSION IHTPER(10),IHTEMP(10),IPTEMP(10)
      EQUIVALENCE (ILIPSE(1),HLIPSE(1),ELIPSE(1))
      EQUIVALENCE (IHTP(1) ,IHTPER(1)),
     *            (IHTP(11),IHTEMP(1)),
     *            (IHTP(21),IPTEMP(1))
      LOGICAL INCUT,BADA
C
C-----------------------------------------------------------------------
C
C THE FIRST 50 LOCATIONS ARE RESERVED FOR SPECIAL WORKING VARIABLES.
C THE NEXT 150 LOCATIONS ARE RESERVED FOR COPIES OF THESE SPECIAL
C    VARIABLES, IF, E.G., ONE WISHES TO PRESERVE THEIR VALUES AT
C    THE BEGINNING OF A STEP.
C THE COPIES HAVE A SUFFIX 0 FOR VALUES AT LAST ASSOCIATED HIT.
C THE COPIES HAVE A SUFFIX 7 FOR VALUES AFTER LEAD-GLASS BARREL OR
C    YOKE END-PLUG (X7 IS LATER IN COMMON - SEE BELOW).
C THE COPIES HAVE A SUFFIX 8 FOR VALUES AT LAST BUT ONE GOOD
C    CHAMBER LAYER.
C THE COPIES HAVE A SUFFIX 9 FOR VALUES AT LAST GOOD CHAMBER LAYER.
C NSPECI (SET IN BLOCK DATA BEHIND MUFFLE) IS THE NUMBER OF SUCH
C    VARIABLES ACTUALLY USED.
C
C SPECIAL VARIABLES....
C
C X,Y,Z       = CURRENT COORDINATES.
C DCX,DCY,DCZ = CURRENT DIRECTION COSINES OF TRACK.
C P,PSQ       = CURRENT MOMENTUM (AND ITS SQUARE).
C E,ESQ       = CURRENT ENERGY (AND ITS SQUARE).
C D           = CURRENT DISTANCE FROM VERTEX (MM).
C GM          = MATERIAL TRAVERSED TO CURRENT POSITION (GM CM**-2).
C AB          = ABSORPTION LENGTHS SO FAR.
C RD          = RADIATION LENGTHS SO FAR.
C DE          = ENERGY LOSS (GEV).
C PPIDK       = PROBABILITY OF PION DECAYING TO MUON BEFORE INTERACTING.
C PKDK        = PROBABILITY OF KAON DECAYING TO MUON BEFORE INTERACTING.
C PPBPEN      = ANTI-PROTON ) ( PENETRATION PROBABILITY, I.E. PROBAB-
C PPPEN       = PROTON      ) ( ILITY OF NOT BEING ABSORBED BY A NUCLEUS
C PPIPEN      = PI          ) ( AND NOT DECAYING.  SEE SUBROUTINE
C PKPEN       = K           ) ( MUFFLS FOR FURTHER DETAILS.
C
C-----------------------------------------------------------------------
C
C THE NEXT 200 LOCATIONS ARE RESERVED FOR VARIABLES THAT DO NOT NEED
C    SPECIAL PRESERVATION.
C
C OTHER VARIABLES...
C
C X1,Y1,Z1    = COORDS OF 1ST POINT ON TRACK.
C X2,Y2,Z2    = COORDS OF LAST POINT ON TRACK.
C XMID,YMID,ZMID = COORDS OF MID POINT ON TRACK.
C X3,Y3,Z3    = COORDS OF POINT WHERE TRACK LEAVES MAGNETIC FIELD.
C D1          = DISTANCE TRAVELLED TO (X1,Y1,Z1).
C D2          = DISTANCE TRAVELLED TO (X2,Y2,Z2).
C D3          = DISTANCE TRAVELLED TO (X3,Y3,Z3).
C DTC         = DISTANCE TRAVELLED TO MID-POINT OF TRACK.
C
C DSTEP       = STEP LENGTH THIS STEP (FOR MUFFLU) (MM).
C GMSTEP      = MATERIAL THIS STEP (GM CM**-2).
C ABSTEP      = ABSORPTION LENGTHS THIS STEP.
C RDSTEP      = RADIATION LENGTHS THIS STEP.
C DESTEP      = ENERGY LOSS OF MINIMUM IONISING PARTICLE THIS STEP,
C
C ADCZ        = ABS(DCZ)
C COSEC       = ABS(COSEC(THETA)), WHERE THETA IS ANGLE TO BEAM.
C STPINI      = STEP LENGTH IN INITIAL TRACKING TO COIL OR END PLUG.
C CURV        = CURVATURE FROM 'PATR' BANK.
C DANG        = ANGLE OF TURN FOR EACH STEP.
C DX,DY,DZ    = CHANGES OF X,Y,Z.
C R,RSQ       = DISTANCE FROM BEAM (RADIUS) AND ITS SQUARE.
C PT          = MOMENTUM COMPONENT PERPENDICULAR TO BEAM.
C DPINA       = INTEGRAL OF (DECAY LENGTH)*(PROBABILTY OF NOT INTER-
C                 ACTING) FOR PION. (NO LONGER USED. CAN BE RESURECTED
C DKNA        = SIMILARLY FOR KAON. (BY UN-COMMENTING IN MUFFLS.
C
C VX,VY,VZ    = CURRENT VARIANCES ON X,Y,Z.
C RMSXY       = RMS DEV. NORMAL TO CHORD IN XY FIT.     )
C RMSZ        = RMS Z DEV. IN RZ FIT.                   ) FROM
C VTXYD       = VARIANCE ON DEVIATION NORMAL TO         ) TRACK
C                 TRACK IN XY PLANE AT TRACK CENTRE.    ) RECON-
C VTZD        = VARIANCE ON DEVIATION NORMAL TO         ) STRUCTION
C                 TRACK IN RZ PLANE AT TRACK CENTRE.    ) ERRORS.
C VTXYAN      = VARIANCE ON ANGLE IN XY PLANE (PHI)     )
C VTZANG      = VARIANCE ON ANGLE IN RZ PLANE (THETA)   )
C VMSANG      = VARIANCE OF MULTIPLE SCATTERING ANGLE.
C VMSD        = VARIANCE OF MULTIPLE SCATTERING DEVIATION NORMAL TO
C                 TRACK.
C CMS         = COVARIANCE OF ABOVE TWO QUANTITIES.
C SD--        = CORRESPONDING STANDARD DEVIATIONS (SQUARE ROOTS).
C
C CHISQT      = CHI-SQUARED OF DEVIATION OF HITS FROM   ) FOR
C                 PROJECTED TRACK.                      ) TRANSVERSE
C NDFT        = NUMBER OF DEGREES OF FREEDOM.           ) (DRIFT)
C CHIPRT      = CHI-SQUARED PROBABILITY.                ) COORDS.
C
C CHISTL      = CHI-SQUARED OF DEVIATION OF HITS FROM   ) FOR
C                 PROJECTED TRACK.                      ) TRANSVERSE
C NDFTL       = NUMBER OF DEGREES OF FREEDOM.           ) & LONG'L
C CHIPTL      = CHI-SQUARED PROBABILITY.                ) COORDS.
C
C NBADCH      = NO. OF BAD CHAMBERS (AS DEFINED BY HMCSTA).
C
C IDTRK       = INNER DETECTOR TRACK NUMBER.
C NCEPTS      = NUMBER OF INTERCEPTS FOUND BY MUREGY.
C NTHIS       = NUMBER OF MUON HITS ASSOCIATED WITH PROJECTED INNER
C                DETECTOR TRACK.
C INEFF       = NUMBER OF INEFFICIENCIES IN MUON SYSTEM FOUND THIS TRACK
C NHLAYR      = NUMBER WHICH REPRESENTS NUMBER OF MUON LAYERS WITH
C                ASSOCIATED MUON HITS. ADD 1 FOR THE INNER LAYER AND
C                2 FOR EACH SUBSEQUENT LAYER.
C INFLAG      = 1 IF LAST LAYER INTERCEPTED WAS INEFFICIENT.
C NELIPS      = NUMBER OF MULTIPLE SCATTERING ELIPSES ENTERED IN ELIPSE.
C NGLAYR      = AS NHLAYR, BUT FOR 'GOOD' LAYERS WHETHER THERE WAS AN
C                ASSOCIATED HIT OR NOT.
C
C-----------------------------------------------------------------------
C
C THE NEXT 1400 LOCATIONS USED BY MUREGY (CALCULATES REGION INTERCEPTS).
C
C IREG        = REGION NUMBER.                           )(IN PAIRS,I.E.
C RDOTD       = R.DOT.(DIR. COSINES) FOR EACH INTERCEPT. )( 2 INTERCEPTS
C CINT        = X,Y,Z COORDINATES OF EACH INTERCEPT.     )( PER REGION.
C NFLAG       = USED INTERNALLY BY MUREGY.
C ITEM        = USED INTERNALLY BY MUREGY.
C
C-----------------------------------------------------------------------
C
C THE NEXT 3600 LOCATIONS USED BY MUFFLE, ETC., TO ACCUMULATE
C   QUANTITIES NEEDED FOR CALCULATING CHI-SQUARED. IN THE DESCRIPTION
C   BELOW I,J GO AS FOLLOWS..
C     1   LEFT DRIFT COORDINATE    )  OF 1ST
C     2   RIGHT DRIFT COORDINATE   )  ASSOCIATED
C     3   LONGITUDINAL COORDINATE  )  HIT
C     4    ETC. FOR NEXT ASSOCIATED HIT.
C   K REFERS TO K'TH ASSOCIATED HIT.
C
C EM(I,J)     = 'ERROR MATRIX', I.E. VARIANCE (I=J) OR COVARIANCE (I.NE.
C                J) FOR COORDINATE I,J.
C DALONG(K)   = DISTANCE OF ASSOCIATED HIT K ALONG TRACK.
C DEV(I)      = DEVIATION OF COORDINATE I FROM EXTRAPOLATED INNER
C                DETECTOR TRACK.
C VXYDA(K)    = VARIANCE ON DEVIATION IN XY PLANE FOR HIT K.
C VZDA(K)     = VARIANCE ON DEVIATION IN RZ PLANE FOR HIT K.
C CXYA(K)     = COVARIANCE BETWEEN DEVIATION AND ANGLE IN XY PLANE.
C CZA(K)      = COVARIANCE BETWEEN DEVIATION AND ANGLE IN RZ PLANE.
C INCUT(I)    = .TRUE. IF COORDINATE I IS IN CUT (SEE MUFFLY).
C BADA(I)     = .TRUE. IF COORDINATE I IS BAD.
C G           = INVERSE ERROR MATRIX.
C DEVA        = CORRESPONDING DEVIATIONS.
C IPERM(K)    = USED TO PERMUTE LEFT/RIGHT AMBIGUITIES.
C MWORK       = USED BY MATRIX INVERTING ROUTINE AS WORKING SPACE.
C ITHISA(K)   = HIT NUMBER OF K'TH ASSOCIATED HIT.
C IJA(K)      = INDEX FOR K'TH ASSOCIATED HIT FOR ENTRY IN BANKS 2 & 3.
C EM1         = COPY OF RELEVANT PART OF ERROR MATRIX FOR TEST PURPOSES.
C EMG         = PRODUCT OF ERROR MATRIX AND ITS INVERSE. SHOULD BE UNIT
C                MATRIX. USED FOR TEST PURPOSES.
C ELIPSE      = TEMPORARY STORAGE FOR MULTIPLE SCATTERING ELIPSES.
C INCHAM      = LIST OF INEFFICIENT CHAMBERS ON THIS TRACK.
C INREG       = LIST OF INEFFICIENT REGIONS (PLANES) ON THIS TRACK.
C IAPPRO      = USED AS WORKING SPACE IN MUFFLX.
C HBADCH      = LIST OF BAD CHAMBERS (AS DEFINED BY HMCSTA).
C CTDASH(K)   = TRANSVERSE COORD. ESTIMATED FROM I.D. PROJECTION.
C CLDASH(K)   = LONGITUD'L COORD. ESTIMATED FROM I.D. PROJECTION.
C VTDASH(K)   = VARIANCE ON CTDASH.
C VLDASH(K)   = VARIANCE ON CLDASH (FOR K'TH ASSOCIATED HIT).
C
C-----------------------------------------------------------------------
C
C THE NEXT  380 LOCATIONS USED FOR VARIOUS IMPROVEMENTS TO MUFFLE, ETC.
C
C IHTREG      = FOR EACH ASSOCIATED HIT, THE REGION NO. IN WHICH
C               THE HIT OCCURS.
C LAYRAB      = ABSOLUTE NO. OF LAYERS WHICH HAVE ASSOCIATED HITS.
C IEFFRG      = LIST OF REGION NOS. OF THOSE REGIONS WHICH HAVE AN
C               ASSOCIATED HIT.
C   IHTPER = FLAG FOR EACH ASSOCIATED HIT TO TELL IF IT'S TO
C            GO FORWARD INTO THE CHI**2 CALCULATION . ONLY ONE
C            HIT PER LAYER IS TO BE USED ,APART FROM OVERLAP HITS.
C   IHTEMP = LIST OF HITS WHICH ARE TEMPORARILY IN THE HIT PERMUTATION
C            DUE TO OVERLAPS WITH THE CURRENT AMBIGUITYOF ONE OF THE
C            HITS IN IHTPER.
C   IPTEMP =  AS IHTEMP, BUT FOR L/R AMBIGUITIES.
C IOVLAP   =  MATRIX OF OVERLAP HITS. IOVLAP(I,J)=1 MEANS THAT THERE
C             IS AN OVERLAP BETWEEN ASSOCIATED HITS I & J ( I & J BOTH
C             IN THE RANGE 1 TO NTHIS ).
C AAA(K)      = TRANSFORMATION COEFF. A FOR K'TH ASSOCD. HIT.
C BBB(K)      = TRANSFORMATION COEFF. B FOR K'TH ASSOCD. HIT.
C CCC(K)      = TRANSFORMATION COEFF. C FOR K'TH ASSOCD. HIT.
C DDD(K)      = TRANSFORMATION COEFF. D FOR K'TH ASSOCD. HIT.
C
C DOVER, ETC. = STORES QUANTITIES LEFT OVER AFTER LEAVING LAST ABSORBER.
C
C X7          = VALUES OF SPECIAL VARIABLES AFTER LEAD-GLASS BARREL OR
C               YOKE END-PLUG - USED FOR RE-FITTING, ETC.
C
C-------------END OF MACRO CMUFWORK-------------------------------------
      COMMON /CMULOC/ ALLDEV(12,2),ALLABS(12),REGDEV(30),REGDCE(60),
     +                REGZE(60),LOCHRS(60),ABLOCH(60)
C
      DIMENSION HAMB(1)
      DIMENSION LHIT(10),LHC(10,10),LAMB(10),LPERM(10),LAYER(10)
      DIMENSION IHTMIN(10),ILRMIN(10),LOPERM(10),LOHPER(10)
C
C-----------------------------------------------------------------------
C             DEFINITIONS
C-----------------------------------------------------------------------
C
C   LHIT   = THE NO. OF ASSOCIATED HITS IN EACH OF THE 'LAYRAB'
C            LAYERS WHICH HAVE ASSOCIATED HITS.
C   LHC    = LAYER-HIT CORRELATION ARRAY :
C                  LHC(I,J) CONTAINS THE HIT NO. ( WITHIN 1 --> NTHISS )
C                           OF THE JTH. ASSOCIATED HIT IN THE ITH.
C                           LAYER WHICH HAS ASSOCIATED HITS.
C   NAMB   = THE NO. OF LAYERS WHICH HAVE > 1 ASSOCIATED HIT.
C   LAMB   = ARRAY CONTAINING THE INDEX ( WITHIN 1--> LAYRAB ) OF
C            THOSE LAYERS WHICH HAVE > 1 ASSOCIATED HIT.
C   LAYER  = FOR ASSD EACH HIT ( LABELLED 1 --> NTHISS ), ITS LAYER
C            NO. ( LABELLED 1 --> LAYRAB ) .
C
C   IHTPER = ARRAY (1 --> NTHISS ) WITH 1 OR 0 TO FLAG WHETHER A HIT IS
C            IN THE CURRENT HIT PERMUTATION
C   IHTEMP = ARRAY (1 --> NTHISS ) WITH 1 OR 0 TO FLAG WHETHER A HIT IS
C            AN ADDITIONAL, OVERLAP HIT IN THE CURRENT HIT PERMUTATION
C   LOHPER = IHTPER + IHTEMP OF CURRENT LOWEST CHI**2 /NDF FOR CURRENT
C            HIT PERMUTATION.
C   IPERM  = ARRAY (1 --> NTHISS ) WITH 1 OR -1 TO FLAG AMBIGUITY OF
C            HITS IN THE CURRENT PERMUTATION. 0 = HIT NOT IN THIS PERM.
C   IPTEMP = ARRAY (1 --> NTHISS ) WITH 1 OR -1 TO FLAG AMBIGUITY OF
C            ADDITIONAL OVERLAP HITS IN THE CURRENT PERMUTATION.
C   LOPERM = IPERM + IPTEMP OF CURRENT LOWEST CHI**2 /NDF FOR CURRENT
C            HIT PERMUTATION.
C   IHTMIN = ARRAY (1 --> NTHISS ) WITH IHTPER OF FINAL SOLUTION
C   ILRMIN = ARRAY (1 --> NTHISS ) WITH LOPERM OF FINAL SOLUTION
C
C            N.B. LOPERM * LOHPER ARE REQUIRED IN THIS VERSION, WHERE
C                 OVERLAP HITS ARE USED IN THE CHI**2 CALCULATION.
C                 SIMILARLY, FOR THE L/R AMBIGUITY PERMUTATION,
C                 CHI**2/NO. OF DEGREES OF FREEDOM IS NOW THE VARIABLE
C                 WHICH DEFINES THE BEST PERMUTATION.
C
C                 ALSO NOTE THAT NTHISS = NTHIS UNLESS NTHIS > 10
C
C--------- M E T H O D   O F   S O L U T I O N -------------------------
C
C  CHOVND = CHI**2 / NO. OF DEGREES OF FREEDOM (DRIFT COORDS ONLY)
C
C  THE SOLUTION IS THAT HIT * L/R PERMUTATION WITH LOWEST CHOVND.
C
C  DCHISQ / IDNDF = DRIFT CHI**2 / NDF OF SOLUTION
C  ZCHISQ / IZNDF = DRIFT + LONGL CHI**2 / NDF OF SOLUTION    (NOT DONE)
C
C-----------------------------------------------------------------------
C
C             INITIALISE THE ABOVE VARIABLES FOR THIS TRACK.
C
      CALL VZERO (REGDEV,30)
      CALL VZERO (IHTPER,10)
      CALL VZERO (IHTEMP,10)
      CALL VZERO (IPTEMP,10)
      CALL VZERO (LHIT  ,10)
      CALL VZERO (LHC  ,100)
      CALL VZERO (LAMB  ,10)
      CALL VZERO (LAYER ,10)
      CALL VZERO (IHTMIN,10)
      CALL VZERO (ILRMIN,10)
      CALL VZERO (LOPERM,10)
      CALL VZERO (LOHPER,10)
      NAMB   = 0
      CHOVND = 1.E19 ! 1.E50  PMF 06/11/98
      DCHISQ = 0.
      IDNDF  = 0
C
C++++
C     WRITE(6,2857)DALONG,DEV,VXYDA,VZDA,CXYA,CZA,INCUT,BADA,ITHISA,IJA,
C    *LAYRAB,IEFFRG,IHTREG
C
C2857 FORMAT('0DALONG',T20,10F10.1/
C    *' DEVIATIONS',T20,10F10.1/2(T20,10F10.1/)
C    *' VXYDA',T20,10F10.0/
C    *' VZDA',T20,10F10.0/
C    *' CXYA',T20,10F10.0/
C    *' CZA',T20,10F10.0/
C    *' INCUT',T20,30L2/
C    *' BADA',T20,30L2/
C    *' ITHISA',T20,10I5/
C    *' IJA',T20,10I5/
C    *' LAYRAB',I3,'  IEFFRG',10I5/
C    *' IHTREG',10I5)
C
C     DO  2849  I = 1,10
C       WRITE(6,2858) (IOVLAP(I,J),J=1,10)
C2858   FORMAT(10I2)
C2849 CONTINUE
C++++
C
C
      IF(LAYRAB.GT.10) LAYRAB = 10
      NTHISS = NTHIS
      IF(NTHIS.GT.10)  NTHISS = 10
      N2 = 2 * NTHISS
      N3 = 3 * NTHISS
C
C             FILL NAMB,LHIT * LHC.
C
      DO 10  I = 1,LAYRAB
C
        DO  5  J = 1,NTHISS
          IF(IHTREG(J).NE.IEFFRG(I)) GO TO 5
C+++
          IF(IHTREG(J).LT.0)GO TO 5
C+++
          LHIT(I) = LHIT(I) + 1
          LHC(I,LHIT(I)) = J
          LAYER(J) = I
    5   CONTINUE
C
CCCC    IF(LHIT(I).EQ.1) GO TO 8
C+++
        IF(LHIT(I).LE.1) GO TO 8
C+++
          NAMB       = NAMB + 1
          LAMB(NAMB) = I
          GO TO 10
    8   CONTINUE
C
C              FILL IHTPER FOR THOSE HITS IN
C              LAYERS WITH ONLY 1 ASSOCIATED HIT.
C
CCCC    IHTPER( LHC(I,1) ) = 1
C+++
        IF(LHIT(I).EQ.1)IHTPER( LHC(I,1) ) = 1
C+++
   10 CONTINUE
C
C
C     DO  2010  I = 1,10
C       WRITE(6,2009)(LHC(I,J),J=1,10),LHIT(I)
C2009   FORMAT('0 ',10I3,'  LHIT=',I3)
C2010 CONTINUE
C     WRITE(6,2015)NAMB,LAMB,IHTPER
C2015 FORMAT('0 NAMB=',I3,'  LAMB=',10I3/'0 IHTPER=',10I3)
C
C              FILL OUT REST OF ERROR MATRIX.
C
      N21 = N2 - 1
      DO  1  J = 1,N21
        J1  = J + 1
        DO  3  I = J1,N2
          EM(I,J) = EM(J,I)
 3      CONTINUE
 1    CONTINUE
C
C     WRITE(6,2856)((EM(I,J),I=1,12),J=1,12)
C2856 FORMAT('0FIRST 12 ENTRIES IN ERROR MATRIX'/(12F10.0))
C
C
C-----------------------------------------------------------------------
C
C             BEFORE GETTING FIRST PERMUTATION OF HITS IN AMBIGUOUS
C             LAYERS, INITIALISE LPERM TO ZERO.
C
      CALL VZERO(LPERM,10)
C
C             IF THERE ARE NO LAYERS WITH MORE THAN 1 ASSOCIATED HIT,
C             JUST PERMUTE THE LEFT/RIGHT AMBIGUITIES.
C
      IF(NAMB.EQ.0) GO TO 27
C
C-----------------------------------------------------------------------
C
C                  LOOP OVER PERMUTATIONS OF HITS
C                  ------------------------------
C
C             ALWAYS GO THRO 13 TO GET NEXT PERMUTATION OF HITS.
C
C LPERM(I) ( I=1,NAMB ) COUNTS HOW MANY TIME AMBIGUOUS LAYER 'I' HAS HAD
C             ITS HIT REPRESENTATIVE CHANGED IN THE PERMUTING ALGORITHM
C             SINCE THE (I-1)TH. LAYER LAST CHANGED ITS REPRESENTATIVE.
C             ALL PERMUTATIONS ARE EXHAUSTED WHEN :
C
C             LPERM(I)=LHIT( LAMB(I) )        FOR ALL I = 1,NAMB.
C
C-----------------------------------------------------------------------
C
  13  CONTINUE
C
C             GET NEXT PERMUTATION OF HITS.
C
      DO  25  I = 1,NAMB
        ILAYR = LAMB(I)
        NHL   = LHIT(ILAYR)
C
C             IF THIS IS THE LAST AMBIG. LAYER, CYCLE THRO' ALL ITS HITS
C             WITH THE CURRENT REPS OF AMBIG. LAYERS 1--> NAMB-1  FIXED.
C
        IF(I.EQ.NAMB) GO TO 17
C
C             THIS LAYER ONLY CHANGES ITS REP IF, FOR ITS CURRENT REP,
C             THE AMBIGUOUS LAYERS WITH HIGHER I HAVE ALL GONE THRO' ALL
C             POSSIBLE PERMUTATIONS I.E. LPERM(J)=LHIT(LAMB(J))
C                                                ALL J > I  ( COND A )
          I1 = I + 1
          DO  15  J = I1,NAMB
            IF(LPERM(I).NE.0.AND.LPERM(J).LT.LHIT(LAMB(J))) GO TO 25
 15       CONTINUE
C
C             CHANGE THE REP OF THIS AMBIG LAYER.
C
 17     DO  21  J = 1,NHL
          KHIT = LHC(ILAYR,J)
          IF(J.LT.LPERM(I)) GO TO 21
C******* THE STATEMENT BELOW MAY BE REDUNDANT ********
            IF(IHTPER(KHIT).EQ.1)GO TO 19
              IHTPER(KHIT)=1
              J1=J
   18         J1=J1+1
              IF(J1.GT.NHL)GO TO 23
                KHIT = LHC(ILAYR,J1)
                IHTPER(KHIT) = 0
              GO TO 18
   19       IHTPER(KHIT)=0
   21   CONTINUE
C
   23   LPERM(I)=LPERM(I)+1
        IF(I.EQ.NAMB)GO TO 25
C
C             COND A MUST HAVE BEEN SATISFIED FOR THIS LAYER. THEREFORE,
C             START ALL OVER AGAIN WITH THIS REP OF LAYER 'I'. SO , FOR
C             ALL J > I  ( UP TO NAMB ), MUST PUT LPERM(J)=0.
C
        DO  24  J = I1,NAMB
          LPERM(J) = 0
   24   CONTINUE
C
   25 CONTINUE
C
   27 CONTINUE
C
C++++
C     WRITE(6,2025)LPERM,IHTPER
C2025 FORMAT('0 LPERM=',10I3/'0 IHTPER=',10I3)
C++++
C
C-----------------------------------------------------------------------
C
C             CALCULATE CHI-SQUARED FOR DRIFT COORDINATES ALONE.
C
      CHINLO = 1.E19 ! 1.E50 PMF 06/11/98
      NDFTMI = 0
      NPERM  = 1
C
C             SET UP LEFT/RIGHT PERMUTATION VECTOR IPERM
C             IPERM = 0 IF NO AMBIGUITY.  INITIALLY SET TO -1 TO SELECT
C             LEFT AMBIGUITY IF AMBIGUITY EXISTS. THEN EACH NON-ZERO
C             VALUE OF IPERM IS SET TO +1 UNTIL ALL PERMUTATIONS HAVE
C             BEEN TRIED.
C
      DO  2  I = 1,NTHISS
        I1 = 3 * I - 2
        I2 = I1 + 1
        IT = 2 * I - 1
        IPERM(I) = 0
        IF(BADA(IT)) GO TO 2
C+++
        IF(IHTREG(I).LT.0)GO TO 2
C+++
        IF(IHTPER(I).NE.1) GO TO 2
        IF(INCUT(I1).AND.INCUT(I2)) IPERM(I) = -1
 2    CONTINUE
C
C             START OF LEFT/RIGHT AMBIGUITY PERMUTATION LOOP.
C
 30   CONTINUE
C
C             CLEAR IHTEMP , IPTEMP.
C
      CALL VZERO(IHTEMP,10)
      CALL VZERO(IPTEMP,10)
C
C             FOR THIS PERMUTATION OF HITS * AMBIGUITIES,
C             SELECT APPROPRIATE BITS OF DEVIATION VECTOR.
C
      CALL VZERO(IAPPRO,N3)
      DO  35  I = 1,NTHISS
        I1 = 3 * I - 2
        IT = 2 * I - 1
C
C             IS THIS HIT IN THE PRESENT HIT PERMUTATION ?
C
        IF(IHTPER(I).NE.1) GO TO 35
        IF(BADA(IT))       GO TO 35
        I2 = I1 + 1
        IF(IPERM(I)) 33,31,32
 31       IF(INCUT(I1))      GO TO 33
 32       II = I2
          JSIGN = 1
        GO TO 34
 33       II = I1
          JSIGN = -1
 34     IAPPRO(II) = 1
C
C             IS THERE AN OVERLAP HIT FOR THIS AMBIGUITY OF THIS HIT ?
C             IF THERE ARE NO AMBIGUOUS LAYERS , NO ;
C             IF THIS HIT IS NOT IN AN AMBIGUOUS LAYER, NO ;
C
        IF(NAMB.EQ.0) GO TO 35
        IF(LHIT(LAYER(I)).LT.2) GO TO 35
C
C             MAYBE, SO LOOK AT THE OVERLAP MATRIX :
C
        DO  345  K = 1,NTHISS
C
C             OF COURSE , ONLY CONSIDER OTHER HITS IN THE SAME LAYER.
C
          IF(LAYER(K).NE.LAYER(I)) GO TO 345
          IF(JSIGN.GT.0) GO TO 343
            IF(IOVLAP(I,K).EQ.0) GO TO 345
            GO TO 344
 343      IF(IOVLAP(K,I).EQ.0) GO TO 345
 344      IHTEMP(K) = 1
          IPTEMP(K) = -JSIGN
C
C                  DONT MAKE ENTRY IN IAPPRO IF BAD HIT
C
          IF(BADA(2*K-1)) GO TO 35
          KK = 3 * K - 1
          IF(JSIGN.EQ.1) KK = KK-1
          IAPPRO(KK) = 1
          GO TO 35
 345    CONTINUE
 35   CONTINUE
C
C             CALL MUFMAT TO GET ERROR MATRIX FOR THIS HIT PERMUTATION.
C             NOTE THAT THIS IS DONE FOR EVERY AMB PERM OF EVERY HIT
C             PERM. THIS IS NECESSARY DUE TO THE FACT THAT CHANGING THE
C             AMB PERM MAY ADD HITS TO THE HIT PERM DUE TO OVERLAPS. IN
C             PRINCIPLE , ONE COULD DO THIS ONLY WHEN AN OVERLAP DOES
C             EXIST. BUT THEN HAVE TO CALCULATE THE ERROR MATRIX AGAIN
C             FOR THE BASIC HIT PERM BEFORE NEXT AMB PERM, SO GAIN IS
C                                                               SMALL.
      CALL MUFMAT(NTHISS,MERR)
      IF(MERR.EQ.1)GO TO 53
C
C             FILL DEVA WITH APPROPRIATE BITS OF EM.
C
      K = 0
      DO  37  I = 1,N3
        IF(IAPPRO(I).EQ.0) GO TO 37
        K = K + 1
        DEVA(K) = DEV(I)
 37   CONTINUE
C
C++++
C     WRITE(6,2859)(DEVA(I),I=1,NDFT)
C2859 FORMAT(' DEVIATIONS'/(12F10.0))
C++++
C
C             K SHOULD EQUAL THE NUMBER OF DEGREES OF FREEDOM
C             IF THE PROGRAM IS WORKING CORRECTLY
C
      IF(K.EQ.NDFT) GO TO 38
        CALL MUERRY('MUFFLX',K,
     +  'SERIOUS--> NO. OF MARKED HITS .NE. NO. OF DEGREES OF FREEDOM^')
C
C             FORM CHI-SQUARED.
C
 38   CHISQT = 0.
      DO  50  I = 1,NDFT
        CH = 0.
C
        DO  51  J = 1,NDFT
          CH = CH + G(I,J)*DEVA(J)
 51     CONTINUE
C
        CHISQT = CHISQT + DEVA(I)*CH
 50   CONTINUE
C
      CHINDF = CHISQT / NDFT
      GO TO 55
 53   CHINDF = 1.1E19! 1.1E50 PMF 06/11/98
      NDFT=0
 55   CONTINUE
C
C++++
C     WRITE(6,2834)CHISQT,NDFT
C2834 FORMAT(' CHI-SQUARED',F10.4,'  NDFT',I4)
C++++
C
C             RECORD MINIMUM DRIFT CHI**2 PER DEGREE OF FREEDOM
C             FOR THIS HIT PERMUTATION.
C
      IF(CHINDF.GE.CHINLO) GO TO 60
      CHISMI = CHISQT
      NDFTMI = NDFT
      CHINLO = CHINDF
      MPERM  = NPERM
      DO  58  I = 1,NTHISS
        LOHPER(I) = IHTPER(I)
        IF(IHTEMP(I).NE.0) LOHPER(I) = IHTEMP(I)
        LOPERM(I) = IPERM(I)
        IF(IPTEMP(I).NE.0) LOPERM(I) = IPTEMP(I)
 58   CONTINUE
C
C             PERMUTE L/R AMBIGUITIES. GO BACK TO 30 IF ANOTHER
C             L/R AMBIGUITY IS STILL TO BE TRIED.
C
 60   NPERM=NPERM+1
      DO  63  I = 1,NTHISS
        IPERM(I) = - IPERM(I)
        IF(IPERM(I).EQ. 1) GO TO 30
 63   CONTINUE
C
      QUOT = 0.
      IF(NDFTMI.GT.0) QUOT = CHINLO
C
C++++
C     WRITE(6,1735)CHISMI,NDFTMI,QUOT
C1735 FORMAT(' CHISMI,NDFTMI',F10.4,I5,'  QUOT=',F10.3)
C++++
C
C-----------------------------------------------------------------------
C
C             CHECK HIT PERMUTATION.
C
      IF(NDFTMI.EQ.0) GO TO 170
C
C             IS THIS A NEW MINIMUM OF CHOVND ?  IF NO THEN JUMP.
C
      IF(QUOT.GT.CHOVND) GO TO 170
C
C             **  N.B. COULD HAVE MORE CONDITIONS ON THE ABOVE SELECTION
C             **       E.G. IF NDFTMI > IDNDF + 2 , ACCEPT PRESENT PERM
C             **          OVER CURRRENT MINIMUM IF QUOT < CHOVND + 0.05.
C
C             YES. SO RECORD CURRENT BEST HIT AND AMBIGUITY PERMUTATION
C                  AND NEW MIN CHOVND * CORRESPONDING QUANTITIES.
C
      DO  165  I = 1,NTHISS
        IHTMIN(I)=LOHPER(I)
        ILRMIN(I)=LOPERM(I)
  165 CONTINUE
C
      CHOVND = QUOT
      DCHISQ = CHISMI
      IDNDF  = NDFTMI
C
C             ARE THERE ANY FURTHER HIT PERMUTATIONS ? IF SO, GO TO 13.
C
 170  IF(NAMB.EQ.0) GO TO 180
      DO  175  I = 1,NAMB
        IF(LPERM(I) .LT. LHIT(LAMB(I)) ) GO TO 13
  175 CONTINUE
C
C             ALL PERMUTATIONS OF HITS ARE EXHAUSTED. CALCULATE THE
C             CHI**2 PROBABILITY.
C
 180  IF(IDNDF.GT.0) CHIPRT = PROB(DCHISQ,IDNDF)
C
C++++
C     WRITE(6,2185)IHTMIN,ILRMIN,CHOVND,DCHISQ,IDNDF,CHIPRT
C2185 FORMAT('0 IHTMIN=',10I3/'0 ILRMIN=',10I3/
C    *'0 CHOVND=',F10.3,' DCHISQ=',F10.3,' IDNDF=',I5,' CHIPRT=',G15.5)
C++++
C
C-----------------------------------------------------------------------
C
C ====>  AT PRESENT THE LONGITUDINAL CHI-SQUARED CALCULATION IS NOT DONE
C
      CHIPTL = 0.
C
C CALCULATE CHI-SQUARED FOR DRIFT AND LONGITUDINAL COORDINATES.
C
C SELECT APPROPRIATE BITS OF DEVIATION VECTOR.
C     CALL VZERO(IAPPRO,N3)
C     NDFTL=0
C     DO 136 I=1,NTHISS
C     I1=3*I-2
C     IT=2*I-1
C IS THIS HIT IN THE PRESENT HIT PERMUTATION ?
C     IF(IHTPER(I).NE.1)GO TO 136
C     IF(BADA(IT))GO TO 135
C     I2=I1+1
C     IF(IPERM(I))133,131,132
C131  IF(INCUT(I1))GO TO 133
C132  II=I2
C     GO TO 134
C133  II=I1
C134  IAPPRO(II)=1
C     NDFTL=NDFTL+1
C135  II=I1+2
C     IL=IT+1
C     IF(BADA(IL))GO TO 136
C     IAPPRO(II)=1
C     NDFTL=NDFTL+1
C136  CONTINUE
C     IF(NDFTL.LE.0)GO TO 152
C
C
C FILL G AND DEVA WITH APPROPRIATE INFORMATION.
C     K=1
C     DO 138 I=1,N3
C     IF(IAPPRO(I).EQ.0)GO TO 138
C     IE=2*((I-1)/3)+1
C     IF(MOD(I,3).EQ.0)IE=IE+1
C     L=1
C     DO 137 J=1,N3
C     IF(IAPPRO(J).EQ.0)GO TO 137
C     JE=2*((J-1)/3)+1
C     IF(MOD(J,3).EQ.0)JE=JE+1
C     G(K,L)=EM(IE,JE)
C     L=L+1
C137  CONTINUE
C     DEVA(K)=DEV(I)
C     K=K+1
C138  CONTINUE
C
C++++
C     WRITE(6,3848)NDFTL
C3848 FORMAT(I5,' COORDS IN DRIFT AND LONG. CALCULATION.',
C    *' SPECIFIC ERROR MATRIX IS')
C     DO 3849 I=1,NDFTL
C     WRITE(6,2858)(G(I,J),J=1,NDFTL)
C3849 CONTINUE
C     WRITE(6,2859)(DEVA(I),I=1,NDFTL)
C++++
C
C INVERT MATRIX.
C     CALL SYMIN4(G,20,NDFTL,MWORK,IER)
C     IF(IER.EQ.0)GO TO 98
C       CALL MUERRY('MUFFLX',IER,'ERROR RETURN FROM SYMIN4.^')
C       RETURN
C98   CONTINUE
C++++
C     WRITE(6,2861)IER
C     WRITE(6,2838)
C     DO 3839 I=1,NDFTL
C     WRITE(6,2840)(G(I,J),J=1,NDFTL)
C3839 CONTINUE
C++++
C
C
C FORM CHI-SQUARED.
C     CHISTL=0.
C     DO 150 I=1,NDFTL
C     CH=0.
C     DO 151 J=1,NDFTL
C151  CH=CH+G(I,J)*DEVA(J)
C150  CHISTL=CHISTL+DEVA(I)*CH
C     GO TO 153
C
C152  CONTINUE
C     CHISTL=0.
C
C153  CONTINUE
C     CHIPTL=PROB(CHISTL,NDFTL)
C++++
C     WRITE(6,3735)CHISTL,NDFTL,CHIPTL
C3735 FORMAT(' CHISTL,NDFTL,CHIPTL',F10.4,I5,G15.5)
C++++
C
C-----------------------------------------------------------------------
C
C             RECORD HIT AND L/R ASSIGNMENTS OF SOLUTION.
C             (IDEVP IS THE POINTER FOR REGDEV NEEDED BY MULOCH)
C
  155 IDEVP = 1
      DO  185  I = 1,NTHISS
        IF(IHTMIN(I).NE.1) GO TO 185
        I1 = 3*I-2
        I2 = I1+1
        IAMB = ILRMIN(I)
        IF(IAMB.NE.0) GO TO 186
C
C             PUT IAMB TO CORRECT VALUE FOR THOSE HITS  WHICH  HAD  ONLY
C             ONE AMB WITHIN  FACTOR*SIGMA ( * THEREFORE HAD IPERM = 0 *
C             ILRMIN = 0 ) .
C
        IAMB = -1
        IF(INCUT(I2)) IAMB =1
C
C             FOR LOCH ('LET-OFF' CHECK), COPY THE DEVIATION OF THE
C             SELECTED HIT AMBIGUITY INTO REGDEV AFTER SAVING THE
C             REGION NUMBER. (WE MAY HAVE 2 HITS PER REGION)
C
  186   REGDEV(IDEVP)   = IHTREG(I)
        REGDEV(IDEVP+1) = ALLABS(I)
        REGDEV(IDEVP+2) = ALLDEV(I,(IAMB+3)/2 )
        IDEVP           = IDEVP + 3
C
C             PRESENT SYSTEM :  IF ONLY LEFT/RIGHT AMB IS WITHIN 3 SIGMA
C             HAMB =  -2 / +2 . IF LEFT/RIGHT AMB IS IN THE BEST HIT/AMB
C                               PERM AND THE OTHER AMB IS WITHIN 3 SIGMA
C             HAMB = -12 / +12 .
C
        J = IJA(I)
        IF(J.LT.0) GO TO 185
        HAMB(J) = 2*IAMB
        IF( INCUT(I1).AND.INCUT(I2) ) HAMB(J) = HAMB(J) + IAMB*10
        IT = 2*I - 1
        IF(.NOT.BADA(IT)) GO TO 185
C
C             THIS IS A BAD HIT WHICH  FINDS  ITSELF  IN  THE  BEST  HIT
C             PERMUTATION. ARE BOTH AMBIGUITIES WITHIN 3 SIGMA OF  TRACK
C             OR WHICH ?   IF BOTH, MAKE HAMB = 22.
C
        IF(INCUT(I1)) HAMB(J) = -2
        IF(INCUT(I2)) HAMB(J) =  2
        IF( INCUT(I1).AND.INCUT(I2) ) HAMB(J) = 22
  185 CONTINUE
C
C++++
C     CALL BPRS('MUR2',2)
C     CALL BPRS('MUR2',3)
C++++
C
      RETURN
      END
