C   20/10/81 303161051  MEMBER NAME  MUFFLS   (JADEMUS)     FORTRAN
C   28/03/81 106171601  MEMBER NAME  MUFFLS   (WORKS)       FORTRAN
C   24/03/80 004021728  MEMBER NAME  MUFFLS   (JADEMUS1)    FORTRAN
C   21/05/79 C9070301   MEMBER NAME  MUFFLS   (S)           FORTRAN
C LAST CHANGE 11.00 16/03/83 HUGH MCCANN  - BACK TO STANDARD M.SC. CONST
C      CHANGE 16.04 17/06/81 JOHN ALLISON - DPINA/DKNA CALCN. COMMENTED.
C      CHANGE 16.03 17/06/81 JOHN ALLISON - TO COPE WITH BACK-TRACKING.
C                              (ACTUALLY REQUIRED NO CHANGE. JA.)
C      CHANGE 08.00 10/04/80 HUGH MCCANN - CHANGE MULT. SCATT. CONSTANT
C      CHANGE 17.28 02/04/80 JOHN ALLISON.
C TAKEN FROM JADEMUS 20.15 24/03/80 TO ADD PIDK AND MULTSC.  JA.
C      CHANGE 01.28 03/07/79 JOHN ALLISON.
      SUBROUTINE MUFFLS(*)
C
      IMPLICIT INTEGER*2 (H)
C
C-----------------------------------------------------------------------
C
C STEPS THROUGH DETECTOR, CALCULATING ENERGY LOSS, MULTIPLE SCATTERING,
C    ETC.
C
C RETURN 1 IF PARTICLE STOPS.
C
C-----------------------------------------------------------------------
C
C COMMONS.
C
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   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
C-----------------------------------------------------------------------
C
C DATA INITILISATION, ETC,
C
C FOR CALCULATING PI AND K DECAY WE NEED A CONSTANT WHICH DEPENDS ON
C   MEANLIFE, ETC.  THE PROBABLITY OF DECAYING IN DISTANCE D IS
C   (BRANCHING RATIO)*D/(GAMMA*(MEANLIFE IN REST FRAME)*(VELOCITY))
C   = (BR)*(MASS)*D/((MOMENTUM)*(C-TAU)) WHERE C-TAU IS VELOCITY OF
C   LIGHT * (MEANLIFE IN REST FRAME) AND IS GIVEN, E.G., IN THE
C   PARTICLE DATA BOOKLET. GAMMA IS THE TIME-DILATION FACTOR.
C FOR THE PION, BR(PI->MU) = 1.,
C    C-TAU = 7804 MM, HENCE BR*M/(C-TAU) = 1.788E-05 GEV MM**-1.
C FOR THE KAON, BR(K->MU) = .64,
C    C-TAU = 3709 MM, HENCE BR*M/(C-TAU) = 8.518E-05 GEV MM**-1.
      DATA CPIDK,CKDK/1.788E-5,8.518E-5/
C
C FOR CALCULATING ABSORPTION PROBABILTY WE NED THE ABSORPTION CROSS-
C   SECTION.  THIS IS GIVEN FOR PROTONS ON VARIOUS NUCLEI IN THE
C   PARTICLE DATA BOOKLET.  WE HAVE USED THIS TO CALCULATE THE
C   ABSORPTION LENGTHS GIVEN IN THE BLOCK DATA STATEMENTS AFTER
C   SUBROUTINE MUFFLE (AND KEPT IN THE SAME MEMBER AS MUFFLE).
C   THEY FIGURES ARE A HIGH ENERGY (20GEV I BELIEVE) LIMIT. WE SHOULD
C   TAKE ACCOUNT OF THE ENERGY DEPENDENCE BUT FOR MOMENTA BELOW
C   ABOUT 1 GEV WE GET TOO MUCH BACKGROUND TO DISTIGUISH MUONS ANYWAY
C   AND SINCE 1 GEV IS 'HIGH ENERGY' IN THIS CONTEXT WE IGNORE THIS
C   EFFECT FOR THE TIME BEING.  HOWEVER WE DO INTRODUCE FACTORS
C   TO ACCOUNT FOR THE (SLIGHTLY) DIFFERENT CROSS-SECTIONS FOR PI,K,AND
C   P-BAR, AS FOLLOWS:
      DATA FACPB,FACPI,FACK/1.2,0.9,0.9/
C
C-----------------------------------------------------------------------
C
C LOOK AT ENERGIES AND DIVIDE STEP IF NECESSARY.
      NSTEP=1
      JSTEP=1
 1    CONTINUE
      IF(DESTEP.LT..05*E.OR.DESTEP.LT..010)GO TO 2
      NSTEP=2*NSTEP
      DSTEP=DSTEP/2.
      GMSTEP=GMSTEP/2.
      ABSTEP=ABSTEP/2.
      RDSTEP=RDSTEP/2.
      DESTEP=DESTEP/2.
      GO TO 1
C
C-----------------------------------------------------------------------
C
C START OF STEPPING LOOP.
 2    CONTINUE
C
C-----------------------------------------------------------------------
C
C CALCULATE ACTUAL ENERGY LOSS.
      DEACT=DESTEP*ESQ/PSQ
C
C-----------------------------------------------------------------------
C
C CHECK IF STOPPED - RETURN 1 IF SO.
      IF(E-DEACT.LT.WMU+.010)RETURN1
C
C-----------------------------------------------------------------------
C
C PROPAGATE HADRON ABSORPTION AND DECAY PARAMETERS.
C   PDK1=PROBABILITY OF DECAYING THIS STEP.
C   PSN1=PROBABILITY OF SNEAKING, I.E. NOT BEING ABSORBEB THIS STEP.
C   (FOR DESCRIPTION OF OTHER VARIABLES SEE CMUFWORK.)
C
C FOR P-BAR.
      PSN1=EXP(-FACPB*ABSTEP)
      PPBPEN=PPBPEN*PSN1
C
C FOR PROTON.
      PSN1=EXP(-ABSTEP)
      PPPEN=PPPEN*PSN1
C
C FOR PION.
      PDK1=CPIDK*DSTEP/P
      PSN1=EXP(-FACPI*ABSTEP)
      PPIDK=PPIDK+PPIPEN*PDK1
      PPIPEN=PPIPEN*PSN1*(1.-PDK1)
C     DPINA=DPINA+DSTEP*PPIPEN
C
C FOR KAON.
      PDK1=CKDK*DSTEP/P
      PSN1=EXP(-FACK*ABSTEP)
      PKDK=PKDK+PKPEN*PDK1
      PKPEN=PKPEN*PSN1*(1.-PDK1)
C     DKNA=DKNA+DSTEP*PKPEN
C
C-----------------------------------------------------------------------
C
C PROPAGATE MULTIPLE SCATTERING.
C   THE FORMULAE FOR PROPAGATION OF MULTIPLE SCATTERING ANGLES, ETC.,
C   ARE JUSTIFIED IN JADE NOTE 47. WE ASSUME THAT, E.G.,
C   RMS SCATTERING ANGLE PROJECTED ONTO A PLANE =
C     ((15 MEV/C)/(P*BETA))*SQRT(NUMBER OF RADIATION LENGTHS).
C   (SEE PARTICLE DATA BOOKLET FOR REFINEMENTS AND REFERENCES.)
C   (SEE CMUFWORK FOR DEFINITION OF VARIABLES.)
      CON=RDSTEP*(0.015*E/P**2)**2
      VMSD=VMSD+(VMSANG+CON/3.)*DSTEP**2+2.*CMS*DSTEP
      CMS=CMS+(VMSANG+CON/2.)*DSTEP
      VMSANG=VMSANG+CON
C
C-----------------------------------------------------------------------
C
C PROPAGATE OTHER QUANTITIES.
C
      X=X+DSTEP*DCX
      Y=Y+DSTEP*DCY
      Z=Z+DSTEP*DCZ
      E=E-DEACT
      ESQ=E**2
      PSQ=ESQ-WMUSQ
      P=SQRT(PSQ)
      D=D+DSTEP
      GM=GM+GMSTEP
      AB=AB+ABSTEP
      RD=RD+RDSTEP
      DE=DE+DESTEP
C
C-----------------------------------------------------------------------
C
C END OF STEPPING LOOP.
      JSTEP=JSTEP+1
      IF(JSTEP.LE.NSTEP)GO TO 2
C
C-----------------------------------------------------------------------
C
 99   CONTINUE
C     IF(DSTEP.GE.0.)GO TO 999
C     WRITE(6,100)X,Y,Z,DCX,DCY,DCZ,E,P,D,GM,AB,RD,DE
C100  FORMAT(' X,Y,Z,DCX,DCY,DCZ,E,P,D,GM,AB,RD,DE',3F7.0,3F7.4,2F7.2,
C    * 2F7.0,3F7.2)
C     WRITE(6,101)PPIPEN,PPIDK,DPINA,PKPEN,PKDK,DKNA,PPPEN,PPBPEN
C101  FORMAT(' PPIPEN,PPIDK,DPINA,PKPEN,PKDK,DKNA,PPPEN,PPBPEN',8F10.4)
C     WRITE(6,102)CON,VMSANG,VMSD,CMS
C102  FORMAT(' CON,VMSANG,VMSD,CMS',4E15.4)
C999  CONTINUE
      RETURN
C     DEBUG INIT        (INSIST ON FORTRANG IF USE THIS.)
      END
