C   01/06/83 306011450  MEMBER NAME  MUFFLB   (JADEMUS1)    FORTRAN
C
C  NEW S/R   13.45  06/05/83  CHRIS BOWDERY
C
C-----------------------------------------------------------------------
      SUBROUTINE MUFFLB
C-----------------------------------------------------------------------
C
C             PROCEDURE FOR BACK-TRACKING AND SETTING-UP 1ST PART
C             FOR MUFFLE STEPS 2 AND 3.
C
C-----------------------------------------------------------------------
C
C             COMMONS
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-------------------------------------
C   20/10/81 408081453  MEMBER NAME  CMUFFL   (S)           FORTRAN
C
C LAST CHANGE 15.00 08/08/84 - CKB PWA     - ADD NALLWM,PMXFLG
C      CHANGE 15.00 03/12/83 - HUGH MCCANN - NEW LEAD GLASS PARAMETERS.
C      CHANGE 08.00 12/10/81 - HUGH MCCANN - TO ADD OVLCUT(FOR OVERLAPS)
C      CHANGE 09.33 17/06/81 - JOHN ALLISON - TO ADD RESOLUTIONS.
C
C-------------START OF MACRO CMUFFL-------------------------------------
C
      COMMON / CMUFFL / XYSTEP,RCOILI,RCOILO,PGMID,WMU,WMUSQ,DUMF1(4),
     +                  NSPECI,XPBAR,XPROT,XPION,XKAON,DUMF2(5),
     +                  TBP,GMBP,ABBP,RDBP,DEBP,
     +                  TJETI,GMJETI,ABJETI,RDJETI,DEJETI,
     +                  TJET1,GMJET1,ABJET1,RDJET1,DEJET1,
     +                  TJET2,GMJET2,ABJET2,RDJET2,DEJET2,
     +                  TJET4,GMJET4,ABJET4,RDJET4,DEJET4,
     +                  TJETO,GMJETO,ABJETO,RDJETO,DEJETO,
     +                  TCOIL,GMCOIL,ABCOIL,RDCOIL,DECOIL,
     +                  ZJETE,TJETE,GMJETE,ABJETE,RDJETE,DEJETE,
     +                  ZEPLG,TEPLG,GMEPLG,ABEPLG,RDEPLG,DEEPLG,
     +                  RLG,ZLG,TLG,GMLG,ABLG,RDLG,DELG,
     +                  ZEP,TEP,GMEP,ABEP,RDEP,DEEP,
     +                  GMG(3),ABG(3),RDG(3),DEG(3),
     +                  VDRES,VLRES,
     +                  REST1,REST2,RESTMX,RESL,RESLMX,FACTOR,OVLCUT,
     +                  ZLGHAF,TLGHAF,ZLGC,GMLGC,ABLGC,RDLGC,DELGC,
     +                  NALLWM,PMXFLG
C
C FOR EXPLANATION SEE BLOCK DATA AFTER SUBROUTINE MUFFLE.
C
C--------------END OF MACRO CMUFFL--------------------------------------
C
C-----------------------------------------------------------------------
C
C             PRESERVE COORDS.
C
      CALL UCOPY(X,X0,NSPECI)
      X3 = X
      Y3 = Y
      Z3 = Z
C
C             CALCULATE DISTANCE  TRAVELLED  FROM  VERTEX  (INTERACTION
C             POINT), INCLUDING CORRECTION FOR CURVATURE.
C
      D3 = SQRT(X3**2 + Y3**2 + Z3**2)
      XYCHOR = SQRT(X3**2 + Y3**2)
      D3 = D3 + (CURV**2 * XYCHOR**3 / 24.) * COSEC
C
C             BACK-TRACK.
C
      X = X3 - DCX * D3
      Y = Y3 - DCY * D3
      Z = Z3 - DCZ * D3
C
C             MAKE A STEP TO MID-POINT OF TRACK. FOR  THIS  FIRST  STEP
C             THE QUANTITIES WHICH DETERMINE  MULTIPLE  SCATTERING  AND
C             ENERGY LOSS ARE SET TO  ZERO  BECAUSE  (A)  THE  MULTIPLE
C             SCATTERING IS DETERMINED BY THE TRACK FITTING  PARAMETERS
C             AT THE MID-POINT, AND (B) THE ENERGY IS DETERMINED BY THE
C             MOMENTUM EFFECTIVELY DETERMINED AT THE MID-POINT.
C
      DSTEP  = D3 / 2.0
      GMSTEP = (GMBP + GMJETI + GMJET1 + GMJET2) * COSEC
      ABSTEP = (ABBP + ABJETI + ABJET1 + ABJET2) * COSEC
      RDSTEP = 0.0
      DESTEP = 0.0
C
      RETURN
      END
