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
#include "cmuffl.for"
#include "cmufwork.for"
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
