! *********************************************************
      SUBROUTINE AVAN(IP1,Z1,A1,XIA2P,T2,T4)
!
! Auxiliary routine for SAMCS
!

C      ------------------------------------------------------------------------------------
C       T2*XIC2         - < tet**2> due to MCS
C       T4*XIC2/PLAB**2 - < tet**4> due to MCS
C      ------------------------------------------------------------------------------------
C       IP1 - projectile type
C      -------------------------------------------------------------------------------------
C       Z1,A1 - charge and atomic mass of target  
C      -------------------------------------------------------------------------------------

      IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)
      COMMON/CSF/Z,A,IP,MOM
      DIMENSION S2(100,21),S4(100,21)
      DIMENSION IPS(21),ISP(8)
      DIMENSION ATBL(100)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      LOGICAL FIRSTCALL/.TRUE./
      SAVE FIRSTCALL,S2,S4      
      PARAMETER (SC=5.06765D0)
      EXTERNAL SEXT
      DATA IPS/1,0,2,2,3,3,4,4,0,4,4,1,0,5,6,7,8,4*0/
      DATA ISP/1,3,5,7,14,15,16,17/
      DATA ATBL/1.00797,4.0026,6.939,9.0122,10.811,12.01115,14.0067, 15
     *.9994,18.9984,20.183,22.9898,24.312,26.9815,28.088,30.9738, 32.064
     *,35.453,39.948,39.102,40.08,44.956,47.90,50.942,51.998, 54.9380,55
     *.847,58.9332,58.71,63.54,65.37,69.72,72.59,74.9216, 78.96,79.808,8
     *3.80,85.47,87.62,88.905,91.22,92.906,95.94,99.0, 101.07,102.905,10
     *6.4,107.87,112.4,114.82,118.69,121.75,127.60, 126.9044,131.30,132.
     *905,137.34,138.91, 140.12,140.907,144.24,147.,150.35,151.98,157.25
     *,158.924,162.50, 164.930,167.26,168.934,173.04,174.97,178.49,180.9
     *48,183.85, 186.2,190.2,192.2,195.08,196.987,200.59,204.37,207.19,2
     *08.980, 210.,210.,222.,223.,226.,227.,232.036,231.,238.03,237.,242
     *., 243.,247.,247.,248.,254.,253./
!
      IF ( FIRSTCALL ) THEN
         EPS=1.D-4  
         NZMAX=100
         IPMAX=8
         DO NZ=1,NZMAX 
            Z=1.D0*NZ
            A=ATBL(NZ)
            DO NIP=1,IPMAX
               IP=ISP(NIP)
               R2=MAX(R2P(IP),RAV2(Z,A)) ! fm**2
               Q1=SQRT(1.D-1/R2)/SC ! lower limit  GeV
               QMAX=SQRT(1.D3/R2)/SC ! upper limit  GeV
               MOM=2 
               CALL SGINT1(SEXT,Q1,QMAX,EPS,S2(NZ,NIP))
               MOM=4 
               CALL SGINT1(SEXT,Q1,QMAX,EPS,S4(NZ,NIP))
            ENDDO
         ENDDO
         FIRSTCALL = .FALSE.
      IP=IP1
      Z=Z1
      A=A1
      ENDIF
      R2=RAV2(Z1,A1)
      R2M=MAX(R2P(IP1),RAV2(Z1,A1))   ! fm**2
      Q12=1.D-1/R2M/SC**2 ! GeV**2
      IZ=Z1
      IPP=IPS(IP1)  
      R2T=R2+R2P(IP1)
      PAR=Q12*R2T*SC**2/3.D0
      T2=LOG(Q12/XIA2P)-1.D0-PAR+S2(IZ,IPP)
      T4=Q12*(1.D0-PAR/2.d0)+S4(IZ,IPP)  
      RETURN
      END
c  *********************************************************
      subroutine BETHE_BLOCH(iel,pavg,dens,dedxm)
!
! Get mean energy loss rate from Bethe-Bloch formula
!
!  iel: index for element arrays
!  pavg: momentum
!  dens: density effect parameter
!  dedxm: energy loss rate in GeV / (g/cm^2)
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
	real*8 ip,mel
!
      data bbcon/0.307075d-3/    ! [GeV/g cm^2]
c
      mel = mass(1)
      e = SQRT(pavg**2 + pmass**2)
      gam = e / pmass
      bet = pavg / pmass / gam
      bg = bet * gam
      bet2 = bet * bet
      zm = zelcomp(iel,krad)
      am = aelcomp(iel,krad)
      ip = ipelcomp(iel,krad)
      fac = bbcon * (ABS(pcharge)/ee/bet)**2 * zm/am
!
      tmax = 2d0*mel*bg*bg/(1d0+2d0*gam*mel/pmass+(mel/pmass)**2)
      arg = 2d0*mel*bg**2*tmax*1d18 / ip**2
      if( arg .le. 0. ) then
        iflag = -25
        go to 900
      end if
      dedxm = fac * ( 0.5d0*LOG(arg) - bet2 - dens/2d0 )
!
900   continue
      end
c  *********************************************************
      function CS_INEL(p,iptyp)
!
! Get total inelastic cross section [mb] for hadron interactions
!  in the range from 0 - 1 GeV/c  (0 - 2 GeV/c for protons)
!
! Cross section data from Rev. Part. Prop. (2000)  
!
      implicit double precision (a-h,o-z)
!
      if( iptyp .eq. 3 ) then        ! pi+
        if( p .le. 0.425 ) then
          cs_inel = 0.
        else if( p.gt.0.425 .and. p.lt.1.000 ) then
          cs_inel = -117.59 + 793.79*p  - 1958.95*p**2 + 2083.40*p**3
     &      - 789.68*p**4
        else
          cs_inel = 10
        end if
!
      else if( iptyp .eq. -3 ) then        ! pi-
        if( p .le. 0.100 ) then
          cs_inel = 7
        else if( p.gt.0.100 .and. p.lt.1.000 ) then
          cs_inel = 241.88 - 5790.80*p  + 51512.1*p**2 - 213391.*p**3
     &      + 466937.*p**4 - 557466.*p**5 + 343528.*p**6 - 85537.9*p**7
        else
          cs_inel = 33
        end if
!
      else if( iptyp .eq. 4 ) then      ! K+
        if( p .le. 0.265 ) then
          cs_inel = 0.
        else if( p.gt.0.265 .and. p.lt.0.550 ) then
          cs_inel = 1
        else if( p.gt.0.550 .and. p.lt.1.000 ) then
          cs_inel = 3.51 - 13.1985*p  + 13.2155*p**2
        else
          cs_inel = 5
        end if
!
      else if( iptyp .eq. -4 ) then      ! K-
        if( p .le. 0.245 ) then
          cs_inel = 60
        else if( p.gt.0.245 .and. p.lt.0.290 ) then
          cs_inel = 61 - 555.6*(p-0.245)
        else if( p.gt.0.290 .and. p.lt.0.365 ) then
          cs_inel = 36
        else if( p.gt.0.365 .and. p.lt.0.405 ) then
          cs_inel = 32 + 28*EXP(-2e4*(p-0.385)**2)
        else if( p.gt.0.405 .and. p.lt.1.000 ) then
          cs_inel = 4437.18 - 58983.5*p  + 325312*p**2 - 955527*p**3
     &      + 1.61593e6*p**4 - 1.5786e6*p**5 + 827978*p**6
     &      - 180518*p**7
        else
          cs_inel = 31
        end if
!
      else if( iptyp .eq. 5 ) then      ! proton
        if( p .le. 0.800 ) then
          cs_inel = 0.
        else if( p.gt.0.800 .and. p.lt.0.975 ) then
          cs_inel = 2
        else if( p.gt.0.975 .and. p.lt.2.000 ) then
          cs_inel = -75.5644 + 109.224*p  - 29.3942*p**2
        else
          cs_inel = 25
        end if
!
      else if( iptyp .eq. -5 ) then        ! pbar
        if( p .le. 0.240 ) then
          cs_inel = 172
        else if( p.gt.0.240 .and. p.lt.1.000 ) then
          cs_inel = 254.486 - 390.950*p  + 209.053*p**2
        else
          cs_inel = 72
        end if
      end if
!
      end
c  *********************************************************
      subroutine DECAY(pm)
c
c  Process particle decays
c
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      real*8 polp(3),jackson
      real*8 ctau(6),ppd(4),ppn(4)
      real*8 pcm(4,3),bet(4),bet0(4)
!
      data ctau/0d0, 658.654d0, 7.8045d0, 3.713d0, 0d0,0d0/  ! TRUE VALUES
!      data ctau/0.,1e12,1e-12,3.713,0.,0./  ! fast pi, slow mu decay  <<<<<<
c
c  Return immediately for e and p
      if( ABS(iptyp).eq.1 .or. ABS(iptyp).ge.5 ) go to 900
c
      if( fastdecay ) then    ! force immediate particle decays
        do i=2,4
        ctau(i) = 1.d-12
        end do 
      end if
!
c  Does this particle decay during this step?
      bg = pm / pmass
      beta = bg / SQRT( 1d0 + bg**2 )
      gamma = bg / beta
      arg = sdid / ( ctau(ABS(iptyp)) * bg )
!
      if( ABS(iptyp).eq.2 .and. declev.eq.5 ) then
! A patch for Bill to suppress muon decay
! Save accumulated normalized proper "time" as "SPINx"
        pol(1) = pol(1) + arg
        go to 900
      end if
!
      probdecay = 1d0 - EXP(-arg)
      rn = RAN1(irnarg)
      if( rn .gt. probdecay ) go to 900
!
! ---------------------------------------------------------
c  Particle decayed
      if( .not. dectrk ) then        ! stop tracking
        iflag = -52
        go to 900
      end if
!
      if( ABS(iptyp).eq.2 .and. (declev.eq.3 .or. declev.eq.4) ) then
! A patch for Bob to suppress muon decay
        go to 900
      end if
!          
! Scale momentum vector if PM doesn't match current PP vector
      call VMAG( pp,pmp )
      do i=1,3
      pp(i) = pp(i) * pm / pmp
      end do
!
! Save the parent-particle-in-LAB-frame information
      pmp = pm
      do i=1,3
      ppd(i) = pp(i)
      end do
      call PANG( ppd,cthp,sthp,cphp,sphp )        
!
! Follow the decay product( daughter particle)
!      ipnum = ipnum + 1
!
!  The equations in subr GLOREN require beta to be negative
      bet0(1) = 0d0
      bet0(2) = 0d0
      bet0(3) = -beta
      bet0(4) = gamma
!
! ---------------------------------------------------------
      if( ABS(iptyp) .eq. 2 ) then    ! Muon decay  mu -> e nu nu
!
!     Get the electron 4-vector in the muon rest frame
        call MUON_DECAY(pcm)
        if( iflag .lt. 0 ) go to 900
!        DO I=1,4
!        PPN(I) = PCM(I,1)
!        END DO
!        CALL OUTPUT_NU(1,PPN,1.d0)
!        IPTYP = 1
!        GO TO 900
!     Boost electron to frame with z axis aligned with muon in the LAB
        call GLOREN( bet0,pcm(1,1),ppd )  
!     Rotate electron momentum into frame with z axis aligned to
!      the problem geometry
        call GDROT( ppd,cthp,sthp,cphp,sphp )
        iptyp = SIGN(1,iptyp)
! n.b. the electron polarization is not implemented
!
! Neutrino production data
        if( neutrino.gt.19 .and. neutrino.lt.100 .and. 
     &     nnudk.gt.0 ) then
          wt = evtwt / nnudk
!
          do 50 in=1,nnudk
          call MUON_DECAY(pcm)
          if( iflag .lt. 0 ) go to 50
!
          call GLOREN( bet0,pcm(1,2),ppn )      ! boosted nu
!
          call GDROT( ppn,cthp,sthp,cphp,sphp )
          if( iptyp .gt. 0 ) then     ! mu+ decay
            id = 27
          else                        ! mu- decay
            id = 26
          end if
          call OUTPUT_NU(id,ppn,wt)
!
          call GLOREN( bet0,pcm(1,3),ppn )      ! boosted anti-nu
          call GDROT( ppn,cthp,sthp,cphp,sphp )
          if( iptyp .gt. 0 ) then     ! mu+ decay
            id = -26
          else                        ! mu- decay
            id = -27
          end if
          call OUTPUT_NU(id,ppn,wt)
50        continue
        end if
!
! ---------------------------------------------------------
      else if( ABS(iptyp) .eq. 3 ) then   ! Pion decay  pi -> mu nu 
!
        call GDECA2( mass(3),mass(2),0.d0,pcm )
!
        if( declev.eq.2 .or. declev.eq.4) then
! A patch for Bob to force production at 90 degrees in pion frame
          pmr = 0.02978d0
          phi = RANV(0.d0,twopi)
          pcm(1,1) = pmr * COS(phi)
          pcm(2,1) = pmr * SIN(phi)
          pcm(3,1) = 0d0
          pcm(4,1) = SQRT( pcm(1,1)**2 + pcm(2,1)**2 + mass(2)**2 )
        end if
!
        epilab = SQRT( pmp**2 + pmass**2 ) ! refers to parent particle
        do i=1,3
        bet(i) = -pp(i) / epilab
        end do
        bet(4) = bet0(4)
!
        call GLOREN( bet,pcm(1,1),ppd )
!
! Don't do the following rotation since we did the boost along
!  the pion direction in the lab and not along the z axis in the lab          
!        call GDROT( ppd,cthp,sthp,cphp,sphp )
        iptyp = SIGN(2,iptyp)
!
! Neutrino production data
        if( neutrino.gt.19 .and. neutrino.lt.100 .and. 
     &     nnudk.gt.0 ) then
          wt = evtwt / nnudk
          if( iptyp .gt. 0 ) then     ! pi+ decay
            id = 36
          else                        ! pi- decay
            id = -36
          end if
!
          do in=1,nnudk
          call GDECA2( mass(3),mass(2),0.d0,pcm )
          call GLOREN( bet0,pcm(1,2),ppn )
          call GDROT( ppn,cthp,sthp,cphp,sphp )
          call OUTPUT_NU(id,ppn,wt)
          end do
        end if
! .........................................................
! Compute muon polarization
!  Including many contributions from Y. Fukui, February 2000
!
        if( spin ) then
! Ref. S. Hayakawa, Phys. Rev. 108:1533, 1957, eqn. 11
          pmucm = SQRT( pcm(1,1)**2 + pcm(2,1)**2 + pcm(3,1)**2 )
          pmulab = SQRT( ppd(1)**2 + ppd(2)**2 + ppd(3)**2 )
          denom = pmulab * pmucm
          emucm = pcm(4,1)
          emulab = ppd(4)
          hayakawa = emulab*emucm/denom - gamma*mass(2)**2/denom
          hayakawa = -hayakawa * SIGN(1,iptyp) ! lab frame
!
! Use Wigner angle method (JCG 25 Oct 99)
          ppilab = pmp
!          cthcm = pcm(3,1) / pmucm
          cthcm = pp(1)*pcm(1,1) + pp(2)*pcm(2,1) + pp(3)*pcm(3,1)
          cthcm = cthcm / pmucm / ppilab
          coswig = epilab*pmucm + cthcm*emucm*ppilab
          coswig = coswig / mass(3) / pmulab
          wigner = -coswig * SIGN(1,iptyp)
!
          if( spintrk .eq. 0 ) then
            pol(1) = 0d0
            pol(2) = 0d0
            pol(3) = hayakawa
            if( jpar .le. nprnt ) then
              write(2,100) jpar,hayakawa,wigner
100           format(' ip,helicity : ',i3,1x,2f10.4)
            end if
            go to 800
!
          else
! If spin tracking is to be done, use instead a normalized
!  3-component spin vector. 
! 
! Consider first the PION rest frame where the helicity is -1 for mu+
! Let POLP be the unit vector along the momentum direction
            do i=1,3
!            polp(i) = -pcm(i,1) / pmucm * SIGN(1,iptyp)
            polp(i) = pcm(i,1) / pmucm 
            end do
!
! Boost the spin from the PION frame to the LAB
            do i=1,3
            bet(i) = pp(i) / epilab     ! pion beta in lab
            end do
            bds = polp(1)*bet(1)+polp(2)*bet(2)+polp(3)*bet(3)
            do i=1,3
            pol(i) = emucm*polp(i)
     &        + (gamma-1d0)*emucm/beta**2*bds*bet(i)
     &        + gamma*pmucm*bet(i)
            pol(i) = pol(i) / mass(2)       ! Hayakawa (8)   
            end do
            call VMAG(pol,smag)
            if( smag .lt. 1.d0 ) then
              iflag = -62
              go to 900
            end if
            s0 = SQRT( smag**2 - 1.d0 )       ! Hayakawa (4)
            if( emulab .lt. gamma*mass(2)**2/emucm ) s0 = -s0
            hlab = s0 * mass(2) / pmulab    ! Hayakawa (9)
            hlab = -hlab * SIGN(1,iptyp)
!
            if( jpar .le. nprnt ) then
              write(2,110) jpar,hayakawa,wigner,hlab
110           format(' ip,helicity : ',i3,1x,3f10.4)
              write(2,115) polp,pol
115           format(' Spi,Slab:',2(3f9.4,3x))
            end if
          end if      ! end of test on spin tracking
!
! Now let POLP be the spin vector in the muon rest frame
! Use Jackson, 3rd, ed, Eq. 11.158
          do i=1,3
          bet(i) = ppd(i) / emulab     ! muon beta in lab
          end do
          gmulab = emulab / mass(2)
          bds = pol(1)*bet(1)+pol(2)*bet(2)+pol(3)*bet(3)
          term = gmulab/(gmulab+1d0)*bds
          do i=1,3
          polp(i) = pol(i) - term*bet(i)
          end do
!
! Ensure proper normalization and load final spin into POL
          call VMAG(polp,smag)
          do i=1,3
          pol(i) = -polp(i) / smag * SIGN(1,iptyp)
!          pol(i) = polp(i) / smag 
          end do
!
! Get helicity from definition
          jackson = pol(1)*ppd(1) + pol(2)*ppd(2) + pol(3)*ppd(3)
!          jackson = -jackson / pmulab * SIGN(1,iptyp)
          jackson = jackson / pmulab
          if( jpar .le. nprnt ) then
            write(2,120) pol,jackson
120         format(' Smu,helicity:',3f9.4,5x,f9.4)
          end if
        end if     ! end of test on spin
! ---------------------------------------------------------
      else if( ABS(iptyp) .eq. 4 ) then   ! kaon decay 
!
! n.b. the muon (and electron) polarization in K decay is ignored        
c   Which K decay mode?
        rn = RAN1(irnarg)
!
        if( rn .le. 0.6351d0 ) then     ! K -> mu nu
          call GDECA2( mass(4),mass(2),0.d0,pcm )
          call GLOREN( bet0,pcm(1,1),ppd )  
          call GDROT( ppd,cthp,sthp,cphp,sphp )
          iptyp = SIGN(2,iptyp)
!
! Neutrino production data
          if( neutrino.gt.19 .and. neutrino.lt.100 .and. 
     &       nnudk.gt.0 ) then
            wt = evtwt / nnudk
            if( iptyp .gt. 0 ) then     ! K+ decay
              id = 46
            else                        ! K- decay
              id = -46
            end if
!
            do in=1,nnudk
            call GDECA2( mass(4),mass(2),0.d0,pcm )
            call GLOREN( bet0,pcm(1,2),ppn )
            call GDROT( ppn,cthp,sthp,cphp,sphp )
            call OUTPUT_NU(id,ppn,wt)
            end do
          end if
! .........................................................
        else if( rn .gt. 0.6351d0 .and. rn .le. 0.8467d0 ) then
!     K -> pi pi0 
          call GDECA2( mass(4),mass(2),0.134976d0,pcm )
          call GLOREN( bet0,pcm(1,1),ppd )  
          call GDROT( ppd,cthp,sthp,cphp,sphp )
          iptyp = SIGN(3,iptyp)
! .........................................................
        else if( rn .gt. 0.8467d0 .and. rn .le. 0.8949d0 ) then
!     K -> e pi0 nu
          call GDECA3( mass(4),mass(1),0.134976d0,0.d0,pcm )
          call GLOREN( bet0,pcm(1,1),ppd )  
          call GDROT( ppd,cthp,sthp,cphp,sphp )
          iptyp = SIGN(1,iptyp)
!
! Neutrino production data
          if( neutrino.gt.19 .and. neutrino.lt.100 .and. 
     &       nnudk.gt.0 ) then
            wt = evtwt / nnudk
            if( iptyp .gt. 0 ) then     ! K+ decay
              id = 47
            else                        ! K- decay
              id = -47
            end if
!
            do in=1,nnudk
            call GDECA3( mass(4),mass(1),0.134976d0,0.d0,pcm )
            call GLOREN( bet0,pcm(1,3),ppn )
            call GDROT( ppn,cthp,sthp,cphp,sphp )
            call OUTPUT_NU(id,ppn,wt)
            end do
          end if
! .........................................................
        else if( rn .gt. 0.8949d0 .and. rn .le. 0.9267d0 ) then
!     K -> mu pi0 nu
          call GDECA3( mass(4),mass(2),0.134976d0,0.d0,pcm )
          call GLOREN( bet0,pcm(1,1),ppd )  
          call GDROT( ppd,cthp,sthp,cphp,sphp )
          iptyp = SIGN(2,iptyp)
!
! Neutrino production data
          if( neutrino.gt.19 .and. neutrino.lt.100 .and. 
     &       nnudk.gt.0 ) then
            wt = evtwt / nnudk
            if( iptyp .gt. 0 ) then     ! K+ decay
              id = 46
            else                        ! K- decay
              id = -46
            end if
!
            do in=1,nnudk
            call GDECA3( mass(4),mass(2),0.134976d0,0.d0,pcm )
            call GLOREN( bet0,pcm(1,3),ppn )
            call GDROT( ppn,cthp,sthp,cphp,sphp )
            call OUTPUT_NU(id,ppn,wt)
            end do
          end if
! .........................................................
        else if( rn .gt. 0.9267d0 .and. rn .le. 0.9440d0 ) then
!     K -> pi pi0 pi0
          call GDECA3( mass(4),mass(3),0.134976d0,0.134976d0,pcm )
          call GLOREN( bet0,pcm(1,1),ppd )  
          call GDROT( ppd,cthp,sthp,cphp,sphp )
          iptyp = SIGN(3,iptyp)
! .........................................................
        else if( rn .gt. 0.9440d0 ) then
! K -> pi pi pi
! n.b. only the first pion is tracked; the other two are ignored
          call GDECA3( mass(4),mass(3),mass(3),mass(3),pcm )
          call GLOREN( bet0,pcm(1,1),ppd )  
          call GDROT( ppd,cthp,sthp,cphp,sphp )
          iptyp = SIGN(3,iptyp)
        end if
      end if      ! end of test on particle type
! ---------------------------------------------------------
!
800   continue
      do i=1,3
      pp(i) = ppd(i)
      end do
      call VMAG(pp,pm)        ! return new momentum value
!
900   continue
      return
      end
!  *********************************************************
      subroutine DEDX(hstp,pavg)
!
!  Ionization energy loss
!  This routine treats continuous energy loss and its fluctuations
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
c
      real*8 mel,kappa,landaur
c  The following function gives the maximum value for the
c   Landau lambda function;  cf. GEANT manual, PHYS332-4
      PARAMETER (P1=.60715,P2=.67794,P3=.52382E-1,P4=.94753,
     +                             P5=.74442,P6=1.1934)
      FLAND(X) = P1+P6*X+(P2+P3*X)*EXP(P4*X+P5)
!
      data c2/0.1571d-4/
      data c3/0.1534d-1/            ! [GeV] 
      data gam1/-0.422785/
!
! Find mean (or most probable) energy loss
      call DEDXMEAN(hstp,pavg,demean)
	dedxm = demean / hstp
!
! ---------------------------------------------------------
!
      de = 0.
      if( lstrag ) then
!  Determine fluctuation of energy loss around mean value
!
        mel = mass(1)
        e = SQRT(pavg**2 + pmass**2)
        gam = e / pmass
        bet = pavg / pmass / gam
        bg = bet * gam
        bet2 = bet * bet
        zm = mparm(1,krad)
        am = aelcomp(1,krad)
        dm = mparm(2,krad)
        densel = zm * densat            ! [# cm^-3]
!
        if( straglev .eq. 1 ) then    !  Gaussian distribution
! Ref. S. Ahlen, RMP 52:142, 1980
          arg = c2 * (ABS(pcharge)/ee)**2 * zm*dm/wtmol
     &      * gam**2 * (1.-0.5*bet2) * hstp
          sige = SQRT( arg )
          de = GAUSS( 0.d0,sige )
          go to 800
        end if
!
! Straggling chi parameter  Geant3 manual, p. 247
        xi = c3 * (ABS(pcharge)/ee/bet)**2 * zm*dm/wtmol * hstp     ! [GeV]
        eta = bg
        ratio = mel / pmass
        emax = (2.*mel*eta**2) / (1.+2.*ratio*gam+ratio**2)
        kappa = xi / emax
        if( kappa .lt. 0.d0 ) then
          iflag = -25
          go to 900
        else if( kappa .eq. 0.d0 ) then     ! allow zero pressure
          go to 800
        end if
        xmean = -bet2 - LOG(kappa) + gam1
        xlamx = FLAND( xmean )
!
        if( straglev .eq. 2 ) then     !  Landau distribution
100       continue
          call LANDAU(xlam)
          if( xlam .gt. xlamx ) go to 100
          de = xi * (xlam - xmean )
!
        else if( straglev .eq. 3 ) then
!  Alternative routine for the Landau distribution
110       continue
          rn = RAN1(irnarg)
! This test on rn exists in a commented-out section of GEANT routine GFLUCT
!          if( rn .gt. 0.980 ) go to 110       ! ???
! LANDAUR is a minimal repackaging of the GEANT routine GLANDR
! The comments at the head of that routine say it is a copy of the
! CERN library routine RANLAN and that it is "for restricted
! Landau distribution". However, there is no mention in the CERNLIB
! write up that RANLAN has any type of TCUT restrictions.
!
          xlam = LANDAUR(rn)
          if( xlam .gt. xlamx ) go to 110
          de = xi * (xlam - xmean )
!
        else if( straglev .eq. 4 ) then   !  Vavilov distribution
          if( kappa .lt. 0.01 ) then  ! Landau limit
120         continue
            call LANDAU(xlam)
            if( xlam .gt. xlamx ) go to 120
            de = xi * (xlam - xmean )
c
          else if( kappa .gt. 10. ) then   !  Gausssian limit
            arg = xi * emax * (1. - 0.5*bet2)
            sige = SQRT( arg )
            de = GAUSS( 0.d0,sige )
c
          else
            rn = RAN1(irnarg)
            xlam = VAVILOV(kappa,bet2,rn)
            de = xi * (xlam - xmean )
          end if      ! end of test on kappa
!
        else if( straglev .eq. 5 ) then   !  restricted fluctuations
! This takes into account the total energy loss
!  for continuous processes for energy transfers less than DCUTx
! RESTRICT_LOSS is different from the other fluctuation routines
! because it returns DE as the TOTAL loss and not the deviation
! from <dE/dx>
          call RESTRICT_LOSS(zm,hstp,pavg,e,dedxm,de)
          if( iflag .lt. 0 ) go to 900
          eloss = de
          go to 805
        end if        ! end of test on straglev
      end if          ! end of test on lstrag
!
! ---------------------------------------------------------
c  Correct the particle 4-momentum
c
800   continue
      eloss = demean + de   ! total kinetic energy loss
!
805   continue
!
      if( spinmatter.eq.1 .and. ABS(iptyp).eq.2 ) then
! Compute spin flip probability, cf. Rossmanith
        pflip = bet*bet * mel/mass(2) * eloss/e
        rn = RAN1(irnarg)
        if( rn .lt. pflip ) then
          do i=1,3
          pol(i) = -pol(i)
          end do
        end if
      end if
!
! Apply the computed energy loss to the actual particle 4-vector
      call VMAG(pp,pm)
      e = SQRT(pm**2 + pmass**2)
      ekin = e - pmass - eloss        ! remaining kinetic energy
      engmin = SQRT(pzmintrk**2 + pmass**2) - pmass
!
      if( ekin .gt. engmin ) then  ! apply loss and continue
        e = e - eloss            
!
      else                         ! bring particle to rest
        e = pmass    
        iflag = -43
        go to 900
      end if
!
      pmnew = SQRT(e**2 - pmass**2)
c
      do i=1,3
      pp(i) = pp(i) * pmnew / pm
      end do
c
      if( dedxm .gt. 0. ) then
        hnextde = MAX( epsf*ekin,engmin ) / dedxm
      else
        hnextde = zstep * scalestep
      end if
      pavg = pmnew
!
900   continue
      end
!  *********************************************************
      subroutine DEDXMEAN(hstp,pavg,demean)
!
! Find mean (or most probable) ionization energy loss
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
	real*8 mel,ip
      data bbcon/0.307075d-3/    ! [GeV/g cm^2]
!
      mel = mass(1)
      e = SQRT(pavg**2 + pmass**2)
      gam = e / pmass
      bet = pavg / pmass / gam
	bet2 = bet * bet
      bg = bet * gam
!
! Get density effect factor
      if( delev .gt. 1 ) then
	  call DENSITY(bg,dens)
	else
	  dens = 0d0
	end if
!
! Determine mean value of energy loss
!  Use Geant3 manual, p. 289, eq. 6
      dedxm = 0d0
!
      do k=1,nelcomp(krad)
      am = aelcomp(k,krad)
	wt = felcomp(k,krad)*am
!
      if( delev.eq.1 .or. delev.eq.2 )then
        call BETHE_BLOCH(k,pavg,dens,dedxk)
!
      else if( delev .eq. 3 ) then      ! Restricted Bethe-Bloch loss
! cf. Review of Particle Properties, 2006, eq. 27.7
! This takes into account the mean energy loss for continuous processes 
!  for energy transfers less than DCUTx
        zm = zelcomp(k,krad)
        ip = ipelcomp(k,krad)
        fac = bbcon * (ABS(pcharge)/ee/bet)**2 * zm/am
        if( ABS(iptyp) .eq. 1 ) then
          tcut = dcute
        else
          tcut = dcutm
        end if
        tmax = 2d0*mel*bg*bg/(1d0+2d0*gam*mel/pmass+(mel/pmass)**2)
        tupper = MIN(tcut,tmax)
!
        arg = 2.*mel*bg**2*tupper*1e18 / ip**2
        if( arg .le. 0. ) then
          iflag = -25
          go to 900
        end if
        dedxk = fac * (0.5*LOG(arg)-0.5*(1.+tupper/tmax)*bet2-dens/2.)
!
      else if( delev .eq. 4 ) then      ! most probable energy loss
! cf. Review of Particle Properties, 2010, eq. 27.10
        zm = zelcomp(k,krad)
        ip = ipelcomp(k,krad) * 1d-6        ! [MeV]
        fac = bbcon/2d0 * (ABS(pcharge)/ee/bet)**2 * zm/am  ! [GeV/g cm^2]
        arg = 2.*mel*1d3*bg**2 / ip
        psi = fac*1d3 * mparm(2,krad) * hstp*1d2    ! [MeV]
        dedxk = fac*(LOG(arg)+LOG(psi/ip)+0.2-bet2-dens)  ! [GeV/g cm^2]
!      
      else if( delev .eq. 5 ) then     ! test mode with dE = const * dz
! Use stepsize along z, not s
        call BETHE_BLOCH(k,pdelev5,dens,dedxk)
      end if
!
      dedxm = dedxm + wt*dedxk     
	end do            ! end of loop on nelcomp
!
      dedxm = dedxm / wtmol * mparm(2,krad) * 1d2       ! [GeV/m]
!
	if( delev .lt. 5 ) then
	  demean = dedxm * hstp           ! [GeV]
	else
	  demean = dedxm * hdid
	end if
!
900   continue
      end
!  *********************************************************
      subroutine DELTA_RAY(hstp,pm)
!
! Simulate elastic scattering from an atomic electron 
!  (production of delta ray)
!
! Adapted from GEANT 3.16/21 routines  GDRSGA and GDRAY
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      emass = mass(1)
      tkin = pmass * ( SQRT(1.+(pm/pmass)**2) - 1. )
      XE=TKIN+ pMASS
      GAM=XE/ pMASS
      GAM2=GAM*GAM
      BET2=1.-1./GAM2
      T=GAM-1.
      ndr = 0           ! number of delta rays
!
! ---------------------------------------------------------
      if( iptyp .eq. -1 ) then        ! electron
        IF(TKIN.GT.DCUTE)THEN
C
          X=DCUTE/(T*EMASS)
          X1=1.-X
C
C     Moller scattering
C
          IF(T.GT.2.*DCUTE/EMASS)THEN
            B1=T*T/GAM2
            B2=(2.*GAM-1.)/GAM2
            SIG=(B1*(0.5-X)+1./X-1./X1-B2*LOG(X1/X))/BET2
            SIG= drnorm *SIG/T
            IF(SIG.GT.0.) then
              xint = 1./SIG / densat / 100. ! [m]
              xm = hstp / xint
              ndr = POIDEV( xm,irnarg )
              if( ndr .le. 0 ) go to 900
            else 
              go to 900
            end if
!            rn = RAN1(idum)
!            prob = 1. - EXP(-hstp/xint)
!            if( rn .gt. prob ) go to 900
          ENDIF
        end if        ! end of test on tkin > dcute
!
        IF(X.GE.0.5) GO TO 900
        CCc=1.-2.*X
C
        do j=1,ndr
  10    rn = RAN1(irnarg)
        E=X/(1.-CCc * rn)
C
        B1=4./(9.*GAM2-10.*GAM+5.)
        B2=T*T*B1
        B3=(2.*GAM2+2.*GAM-1.)*B1
        E1=1.-E
C
        SCREJ=B2*E*E-B3*E/E1+B1*GAM2/(E1*E1)
C
        IF(RAN1(irnarg) .GT. SCREJ) GOTO 10
!
        call DRVEC( pm,xe,t,e )
        end do
!
! ---------------------------------------------------------
      else if( iptyp .eq. 1 ) then       ! positron
        X=DCUTE/(T*EMASS)
        X1=1.-X
        IF(TKIN.GT.DCUTE)THEN
C     Bhabha scattering
C
          IF(T.GT.DCUTE/EMASS)THEN
            Y=1./(GAM+1.)
            Y2=Y*Y
            C=1.-2.*Y
            B1=2.-Y2
            B2=C*(3.+Y2)
            C2=C*C
            B4=C2*C
            B3=C2+B4
            SIG=(1./X-1.)/BET2+B1*LOG(X)
     +         +B2*(1.-X)-B3*(1.-X*X)/2.+B4*(1.-X**3)/3.
            SIG= drnorm *SIG/T
            IF(SIG.GT.0.) then
              xint = 1./SIG / densat / 100. ! [m]
              xm = hstp / xint
              ndr = POIDEV( xm,irnarg )
              if( ndr .le. 0 ) go to 900
            else 
              go to 900
            end if
          ENDIF
        end if        ! end of test on tkin > dcute
!
        IF(X.GE.1.) GO TO 900
        X1=1.-X
!
        do j=1,ndr
20      rn = RAN1(irnarg)
        E=X/(1.-X1*rn)
C
        Y=1./(GAM+1.)
        Y2=Y*Y
        C=1.-2.*Y
        B1=2.-Y2
        B2=C*(3.+Y2)
        C2=C*C
        B4=C2*C
        B3=C2+B4
        B0=GAM2/(GAM2-1.)
C
        SCREJ=(((B4*E-B3)*E+B2)*E-B1)*E+B0
        SCREJ=SCREJ/((((B4*X-B3)*X+B2)*X-B1)*X+B0)
        IF( RAN1(irnarg) .GT.SCREJ) GOTO 20
!
        call DRVEC( pm,xe,t,e )
        end do
!
! ---------------------------------------------------------
      else if( ABS(iptyp).eq.2 .or. ABS(iptyp).eq.5 ) then
! spin 1/2 particles, here muon or proton
!
        IF(TKIN.GT.DCUTM)THEN
          emmu = pmass
          TMAX=2.*EMASS*(GAM2-1.)/
     +     (1.+2.*GAM*EMASS/EMMU+(EMASS/EMMU)**2)
!
          IF(TMAX.GT.DCUTM)THEN
            Y=DCUTM/TMAX
            SIG=(1.-Y+BET2*Y*LOG(Y))/DCUTM
            SIG=SIG+(TMAX-DCUTM)/(2.*XE**2)
            SIG=SIG* drnorm *EMASS/BET2
            IF(SIG.GT.0.) then
              xint = 1./SIG / densat / 100. ! [m]
              xm = hstp / xint
              ndr = POIDEV( xm,irnarg )
              if( ndr .le. 0 ) go to 900
            else 
              go to 900
            end if
          end if      ! end of test on tmax > dcutm
        end if        ! end of test on tkin > dcutm
!
        do j=1,ndr
40      rn = RAN1(irnarg)
        E=1./DCUTM+rn*(1./TMAX-1./DCUTM)
        E=1./E
        SCREJ=1.-BET2*(E/TMAX)
        SCREJ=SCREJ+0.5*(E/ xe)**2        ! extra term for spin 1/2
        rn = RAN1(irnarg)
        IF( rn .GT.SCREJ) GO TO 40
!
        E=E/(T*EMASS)
        call DRVEC( pm,xe,t,e )
        end do
!
! ---------------------------------------------------------
      else if( ABS(iptyp).eq.3 .or. ABS(iptyp).eq.4 ) then
! spin 0 particles, here pion or kaon
! cf. GEANT 3.16 manual, p.290
!
        IF(TKIN.GT.DCUTM)THEN
          emmu = pmass
          TMAX=2.*EMASS*(GAM2-1.)/
     +     (1.+2.*GAM*EMASS/EMMU+(EMASS/EMMU)**2)
!
          IF(TMAX.GT.DCUTM)THEN
            Y=DCUTM/TMAX
            SIG=(1.-Y+BET2*Y*LOG(Y))/DCUTM
            SIG=SIG* drnorm *EMASS/BET2
            IF(SIG.GT.0.) then
              xint = 1./SIG / densat / 100. ! [m]
              xm = hstp / xint
              ndr = POIDEV( xm,irnarg )
              if( ndr .le. 0 ) go to 900
            else 
              go to 900
            end if
          end if      ! end of test on tmax > dcutm
        end if        ! end of test on tkin > dcutm
!
        do j=1,ndr
50      rn = RAN1(irnarg)
        E=1./DCUTM+rn*(1./TMAX-1./DCUTM)
        E=1./E
        SCREJ=1.-BET2*(E/TMAX)
        IF(RAN1(irnarg) .GT.SCREJ) GO TO 50
!
        E=E/(T*EMASS)
        call DRVEC( pm,xe,t,e )
        end do
!
      end if          ! end of test on particle type
!
900   continue
      return
      end
! *********************************************************
      subroutine DENSITY(bg,dens)
!
! Get density effect parameter 
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      x0 = mparm(5,krad)
      x1 = mparm(6,krad)
      sa = mparm(7,krad)
      sc = mparm(8,krad)
      sm = mparm(9,krad)
      dens = 0.d0
!
      if( x0 .eq. 0d0 ) go to 900        ! undefined density parameters
!
      if( bg .le. 0. ) then
        iflag = -25                                       
        go to 900
      end if
!
      xs = LOG10( bg )    
      if( xs.ge.x0 .and. xs.le.x1 )then
        dens = 4.6052 * xs + sa * (x1 - xs)**sm + sc
      else if( xs .gt. x1 )then
        dens = 4.6052 * xs + sc
      end if
!
900   continue
	end
! *********************************************************
	REAL*8 FUNCTION DFL(AM2,AM3)
!
! Auxiliary routine for SAMCS
!
	IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)  
      common/rans/iseed
!	REAL RNDM
	PARAMETER (PI=3.141592653589793227D+00)
	C=AM3/AM2**(1.5)
	U1=SQRT(1.D0+C**2/4.D0)
	ZU=C/2.D0-U1
	IF(ZU.GE.0.D0) THEN
	   U=(C/2.D0+U1)**(1.D0/3.D0)+ZU**(1.D0/3.D0)
	ELSE
	   U=(C/2.D0+U1)**(1.D0/3.D0)-(-ZU)**(1.D0/3.D0)
	ENDIF
	ALP=SQRT(LOG(1.D0+U**2))
!	R1=DBLE(RNDM(-1.))
!	R2=DBLE(RNDM(-1.))
	R1 = RAN1(iseed)
	R2 = RAN1(iseed)
	R12=SIN(2.D0*PI*R1)*SQRT(-2.D0*LOG(R2))
	DFL=SQRT(AM2)*(EXP(ALP*R12)/SQRT(1.D0+U**2)-1.D0)/U
	RETURN
	END
! *********************************************************
      real*8 function DOT(a,b)
!
! 4-vector product of A and B
! Adapted from J. Gallardo 10 April 2000
!
      implicit real*8(a-h,o-z)
      dimension a(4),b(4)
!
      prod=0
!
      do i=1,4
      if(i.eq.1)then
        prod=prod+a(i)*b(i)
      else
        prod=prod-a(i)*b(i)   ! METRIC (+,-,-,-)
      endif
      end do
!
      dot=prod
!
      end
! *********************************************************
      subroutine DRVEC(pm,xe,t,e)
!
! Update particle momentum vector for delta ray energy loss
!
! Adapted from GEANT 3.16/21 routine  GDRAY
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      real*8 p(3)
      real*8 pels(3),gkin(3),p4(3),me,mu
!
!      ndeltas = ndeltas + 1
!
      emass = mass(1)
      EEL=(T*E+1.)*EMASS        ! recoil electron total energy
      TEL=EEL-EMASS             ! kinetic energy
      PEL=SQRT(ABS((EEL+EMASS)*TEL))
      COSTH=(XE*EEL+EMASS*(TEL-XE))/(Pm*PEL)
      IF(COSTH.GE.1.) THEN
        COSTH=1.
        SINTH=0.
      ELSEIF(COSTH.LE.-1.) THEN
        COSTH=-1.
        SINTH=0.
      ELSE
        SINTH=SQRT((1.+COSTH)*(1.-COSTH))
      ENDIF
      PHI    = TWOPI* RAN1(irnarg)
      COSPHI = COS(PHI)
      SINPHI = SIN(PHI)
!
C           Polar co-ordinates to momentum components.
      GKIN(1)=PEL*SINTH*COSPHI
      GKIN(2)=PEL*SINTH*SINPHI
      GKIN(3)=PEL*COSTH
C
      PELS(1)=-GKIN(1)
      PELS(2)=-GKIN(2)
      PELS(3)=Pm -GKIN(3)
C
      do i=1,3
      p4(i) = pp(i)   ! at this point pp is still incident momentum
      end do
! Get angle of incident particle in LAB
      call PANG(p4,cthet,sthet,cphi,sphi)
!
! Rotate scattered particle to correct angle in LAB
      call GDROT(pels,cthet,sthet,cphi,sphi)
!
      if( spintrk .gt. 0 ) then 
        do i=1,3
        p(i) = pels(i)
        end do  
! Rotate the spin such that helicity is conserved
        call ROTATE_SPIN(p)      
!                   
        if( spinmatter.eq.2 .and. ABS(iptyp).eq.2 ) then
! Simulate spin flip in mu-e elastic scattering
          me = mass(1)
          mu = mass(2)
          e1 = SQRT( pm**2 + mu**2 )
          s = mu**2 + me**2 + 2.*e1*me
          ecm = SQRT( s )
! Get CM scattering angle from Jackson, 1st ed, Eq. 12.48
          pms = SQRT( pels(1)**2+pels(2)**2+pels(3)**2 )
          e3 = SQRT( pms**2 + mu**2 )
          cthcm = e3-0.5*(e1+me)*(1.+(mu**2-me**2)/ecm**2)
          arg = (1.-((mu+me)/ecm)**2)*(1.-((mu-me)/ecm)**2)
!          if( arg .ge. 0d0 ) then
!            write(*,*) 'Hello from DRVEC'
            denom = SQRT(arg)
!          else
!            stop 'Negative square root in DRVEC'
!          end if
          denom = 0.5*pm * denom
          cthcm = cthcm / denom
!
          if( cthcm .ge. 1.d0 ) then
            cthcm = 1.d0
          else if( cthcm .le. -1.d0 ) then
            cthcm = -1.d0
          end if
!
          call FLIP_ELASTIC( s,cthcm,pr )
!
          rn = RAN1(irnarg)
          if( rn .lt. pr ) then
            do i=1,3        ! flip spin direction
            pol(i) = -pol(i)
            end do
!            nflips = nflips + 1
          end if
        end if
      end if
!
      do i=1,3        ! load pp with new, rotated scattered momentum
      pp(i) = pels(i) 
      end do
      call VMAG(pp,pm)
!
      return
      end
c  *********************************************************
      subroutine ELMS(hstp,pavg)
c
c  Ionization energy loss
c
c  This routine treats continuous energy loss and its fluctuations
c
! Based on a routine by Simon Holmes, 22 NOV 2004
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
c
!      real demean
!      real*8 mel,ip,kappa,landaur

c elms variables

         real*8 dspran
!         integer jint      
c
         REAL*8 pigreek, c, density
         PARAMETER (pigreek = 3.1415926535d0)
         PARAMETER (c = 2.997924d8)
!         PARAMETER (elmsmass=105.658)
         PARAMETER (density=0.0708d0) !g cm2 
c         CHARACTER (100) trajectoryfolder, logname, stagenum, runnumber !!! NOW IN ICOMMON
!         REAL sran
!         REAL*8 dran
!         INTEGER*2 ihr,imin,isec,ihun
!         REAL*8 sigamin(8)
         REAL*8 meanin(8)
c        DB variables recorded in the directory (in icommon)
c
c        incoming vectors etc
!         INTEGER stage, rstage				!cooling stage and number of stages,  stage to be analysed
         INTEGER nstage				!cooling stage and number of stages,  stage to be analysed
         INTEGER maxtraj
         PARAMETER (maxtraj = 100001)		!max number of trajectories
         INTEGER ntraj				!actual number of trajectories
!         REAL*8 px1(maxtraj), py1(maxtraj), pz1(maxtraj)
         REAL*8 pzoffset
!         REAL*8 x1(maxtraj), y1(maxtraj), z1(maxtraj), t1(maxtraj)
c        input 6D phase space trajectories, and range (px, py, pz, x, y, z, t)
         REAL*8 slabthk, zend				!slab thickness in z (metres), end of this absorber
!         REAL rfboost, rfdelay1, rffreq						!pz peak boost per absorber at time rfdelay*stage. 
!         REAL rfdelay2
!         REAL betam
!         REAL drift									!drift distance between absorbers
!         REAL  forward	!output 6D phase space trajectories and range(z2)
!         REAL*8 px2(maxtraj), py2(maxtraj), pz2(maxtraj)
!         REAL*8 x2(maxtraj), y2(maxtraj), z2(maxtraj), t2(maxtraj)
c
!        working vars
!         INTEGER I, J, K, L, N, Ijunk, Kjunk, Itraj, LL, IL, Jmom,
!     &        KP, KPU, KPD
!         REAL*8 a1, a2, a3, a4, a5, a6, px, py, pz,
!     &        x, y, z, t, p, ke, pnew, dv 
!         REAL*8 deltapl, deltapt, deltae, steplen, vel
         REAL*8 UCpl1,UCpl2,UCpt1,UCpt2,UCang1,UCang2
         REAL*8 UCen1,UCen2     ! store two pt,pl for uncorrelated
         INTEGER UCINT
         REAL*8 llenwt, cphi, sphi, ctheta, stheta
         REAL*8 oldi, oldj, oldk, cphir, sphir, cthetas, sthetas
         REAL*8 v(7)
         EQUIVALENCE (v(1), px), (v(2), py), (v(3), pz), (v(4), x),
     &        (v(5), y), (v(6), z), (v(7), t)	!for trajectory concerned
!         REAL tthk, pmom
!         external ranlux

! ---------------------------------------------------------
          slabthk=hstp  !*** get from icool
          meanin(3)=pp(3)*1000 !*** get from icool
          pzoffset=meanin(3)          
          nstage = 1
          ntraj = 1
          drift = 0
          betam = meanin(3) / dsqrt (meanin(3)**2 + elmsmass**2)
          px=pp(1)*1000
          py=pp(2)*1000
          pz=meanin(3)          !*** get from icool
c          zend = stage * slabthk + (stage - 1) * drift
          zend=slabthk
          meanout = 0.0
          matrixout = 0.0
          meanboost = 0.0
          meanloss = 0.0
          Jmom = 2
          x=0.0
          y=0.0
          z=0.0
          t=0.0
 
!		K = Itraj / 10000       ! not used
c		IF (Itraj - 10000*K .eq. 0) write (6,*) Itraj
		x=0.0
                y=0.0
                z=0.0
                t=0.0
		en1 = dsqrt(px**2 + py**2 + pz**2 + elmsmass**2)
        DO 30000 LL=1,10000
	p = dSQRT(px**2 + py**2 + pz**2)
	ke = dSQRT(p**2 + elmsmass**2) - elmsmass
	IF( zend - z .LE. 1.0E-6 .OR. p .le. 5.0.OR.pz.le.5.0)
     &	GOTO 30010
	!unit vector in old direction
	oldi = px / p
        oldj = py / p
        oldk = pz / p
	!find the closest sample set in momentum
	IF (LL .EQ. 1 .OR. 2*P .GT. mtmdb(KP) + mtmdb(KPU)
     & 	 .OR. 2*P .LT. mtmdb(KP) + mtmdb(KPD)) THEN
	N = npdone
        Jmom = 1
	DO 30100 L = 1, 8
	N = min0((N+1)/2, npdone-Jmom-1)
	IF (P .GT. mtmdb(pdone(Jmom+N))) THEN
	Jmom = Jmom + N
	IF (Jmom .GE. npdone-1) GOTO 30110
	ELSEIF (N .EQ. 1) THEN
	GOTO 30110
	ENDIF
30100 CONTINUE
30110 a1 = 0.5 * (mtmdb(pdone(Jmom))+mtmdb(pdone(Jmom+1)))
	IF (P .GT. a1) THEN
	Jmom = Jmom + 1   !take nearest
	a2 = a1
	a1 = 0.5 * (mtmdb(pdone(Jmom))+mtmdb(pdone(Jmom+1)))
	ELSE
	a2 = 0.5 * (mtmdb(pdone(Jmom-1))+mtmdb(pdone(Jmom)))
	ENDIF
	KP = pdone(Jmom)
        KPU = pdone(Jmom+1)
        KPD = pdone(Jmom-1)
	IF (P .GT. a1 .OR. P .LT. a2) THEN
c		write(6,5721)P, mtmdb(KPD), mtmdb(KP), mtmdb(KPU) 
		write(2,5721)P, mtmdb(KPD), mtmdb(KP), mtmdb(KPU) 
5721	 Format(' Momentum interpol fail ',4F10.3)
		endif
		ENDIF
	IF (ELMSCOR.eq.1)THEN
	!now sample momentum Jmom
	!generate integer between 0 and norm(Jmom)
!	call ranlux(randvw,2)
      randvw(1) = RAN1(irnarg)
      randvw(2) = RAN1(irnarg)
c        dspran=(randvw(1)*1000000000+randvw(2))/1000000000  old technique (flawed)
 	dspran=(IDint(randvw(1)*32768.0D0)+randvw(2))/32768.0D0       ! combine two sp to get one dp
	a6 = normelms(KP) * dspran    
	!find it
	I = firstpoint(KP) 
	N = lastpoint(KP) - I + 1
	DO 30200 L = 1, 20
	N = min0((N+1)/2, lastpoint(KP) - I)
	IF (a6 .GT. probarray(I+N)) THEN
	I = I + N
	IF (I.eq.lastpoint(KP)) GOTO 30210
	ELSEIF (N .EQ. 1) THEN
	GOTO 30210
	ENDIF
	IF (I .GE. lastpoint(KP) .OR. N .EQ. 0) GOTO 30210  !N=0 shouldnt happen, misordered?
30200   CONTINUE
	!check algorithm
30210   IF (I .GT. lastpoint(KP)) THEN
c	write (6,3461) a6, probarray(lastpoint(KP)), I, N
	write (8,3461) a6, probarray(lastpoint(KP)), I, N
3461	 Format(' Select fail1', 2D15.8, 2I6)
	endif
	IF (L.GT.19.OR.N.EQ.0)THEN
c	write (6,3462) a6, L,N,I	 !L=20 would imply a choice of more than 10**6 entries
	write (8,3462) a6, L,N,I	
3462 	Format(' Select fail0', D15.8, 3I6)
	endif
	IF ((a6.LT.probarray(I-1).AND.I.GE.firstpoint(KP)) .OR.
     &	 (a6 .GT. probarray(I+1) .AND. I .NE. lastpoint(KP))) THEN 
c	write(6,3463)a6,probarray(I),probarray(I+1), L, N, I  
	write(2,3463)a6,probarray(I),probarray(I+1), L, N, I  
3463	Format( 'select fail3', 3E15.8, 3I6)
	ENDIF
	deltapt = dptelms(I) * (1.0 + pslope(KP) * (P - mtmdb(KP)))	  !interpolation
	deltae = elmsde(I)  * (1.0 + eslope(KP) * (P - mtmdb(KP)))
	deltapl = p - dsqrt(dabs(p**2 - deltapt**2 + deltae**2
     &        - 2.0 * deltae * (ke + elmsmass) ))
	!interpolation check
	a2 = abs(pslope(KP) * (P - mtmdb(KP)))
	a6 = abs(eslope(KP) * (P - mtmdb(KP)))
	IF (P.lt.100.)THEN
	a3 = 0.10	!check at 10% for low momentum
	ELSE
	a3 = 0.03	!otherwise 3%
	ENDIF
	IF(a2.gt.a3.or.a6.gt.a3)then
c	write(6,3464)P,a2,a6
	write(2,3464)P,a2,a6
3464 	Format(' interpolation too big: mtm, pt2, de', F8.1, 2F8.3)
	endif	
        ELSE !Uncorrolated version
30300   CONTINUE
        DO 30400 UCINT=1,2
        !now sample momentum Jmom
	!generate integer between 0 and norm(Jmom)
!	call ranlux(randvw,2)
      randvw(1) = RAN1(irnarg)
      randvw(2) = RAN1(irnarg)
 	dspran=(IDint(randvw(1)*32768.0D0)+randvw(2))/32768.0D0       ! combine two sp to get one dp
	a6 = normelms(KP) * dspran    
	!find it
	I = firstpoint(KP) 
	N = lastpoint(KP) - I + 1
	DO 30500 L = 1, 20
	N = min0((N+1)/2, lastpoint(KP) - I)
	IF (a6 .GT. probarray(I+N)) THEN
	I = I + N
	IF (I.eq.lastpoint(KP)) GOTO 30510
	ELSEIF (N .EQ. 1) THEN
	GOTO 30510
	ENDIF
	IF (I .GE. lastpoint(KP) .OR. N .EQ. 0) GOTO 30510  !N=0 shouldnt happen, misordered?
30500   CONTINUE
	!check algorithm
30510   IF (I .GT. lastpoint(KP)) THEN
	write (6,3461) a6, probarray(lastpoint(KP)), I, N
	write (8,3461) a6, probarray(lastpoint(KP)), I, N
	endif
	IF (L.GT.19.OR.N.EQ.0)THEN
	write (6,3462) a6, L,N,I	 !L=20 would imply a choice of more than 10**6 entries
	write (8,3462) a6, L,N,I	
	endif
	IF ((a6.LT.probarray(I-1).AND.I.GE.firstpoint(KP)) .OR.
     &	 (a6 .GT. probarray(I+1) .AND. I .NE. lastpoint(KP))) THEN 
	write(2,3463)a6,probarray(I),probarray(I+1), L, N, I  
!	write(8,3463)a6,probarray(I),probarray(I+1), L, N, I  
	ENDIF
	deltapt = dptelms(I) * (1.0 + pslope(KP) * (P - mtmdb(KP)))	  !interpolation
	deltae = elmsde(I)  * (1.0 + eslope(KP) * (P - mtmdb(KP)))
	deltapl = p - dsqrt(dabs(p**2 - deltapt**2 + deltae**2
     &        - 2.0 * deltae * (ke + elmsmass) ))
	!interpolation check
	a2 = abs(pslope(KP) * (P - mtmdb(KP)))
	a6 = abs(eslope(KP) * (P - mtmdb(KP)))
	IF (P.lt.100.)THEN
	a3 = 0.10	!check at 10% for low momentum
	ELSE
	a3 = 0.03	!otherwise 3%
	ENDIF
	IF(a2.gt.a3.or.a6.gt.a3)then
!	write(6,3464)P,a2,a6
	write(2,3464)P,a2,a6
	endif
        IF (UCINT.EQ.1) then
        UCPL1=deltapl
        UCPT1=deltapt
        UCen1=deltae
        UCang1=atan(deltapt/p)
        end if
        IF (UCINT.EQ.2) then
        UCPL2=deltapl
        UCPT2=deltapt
        UCen2=deltae
        UCang2=atan(deltapt/p)
        end if
30400   CONTINUE
c make uncorrelated - check!!!
        deltapt= P*tan(UCang2)
        deltapl= p - dsqrt(dabs(p**2 - deltapt**2 + UCen1**2
     &        - 2.0 * UCen1 * (ke + elmsmass) ))
        pnew = dsqrt((p - deltapl)**2 + deltapt**2)  
	IF(Pnew .lt.p*(1.0001))then
        GOTO 30310
        end if
        write(2,*)'increase p',pnew,p
        GOTO 30300
30310   CONTINUE
        ENDIF	
        !step shorten if it exits
	steplen = thkdb(KP) 
	llenwt = (zend - z) * p / (pz * steplen)
	IF(llenwt.le.0.0)then
	write(2,*)'negative lenwt',a2,llenwt,mtmdb(Jmom)
	goto 30010
	endif
	IF (llenwt .LT.1.0 .AND. pz .GT. 0.0) THEN				
	deltapl = deltapl * llenwt
        deltapt = deltapt * llenwt
  	  steplen = steplen * llenwt
		ENDIF  
			!combine this step, assuming that the energy loss and scattering occurs at the END of the step
		x = x + steplen * oldi
		y = y + steplen * oldj	    
		z = z + steplen * oldk
!		t = t + steplen * dsqrt(ivel2) / c        ! not used
		forward = sign(1d0, p - deltapl)
		pnew = dsqrt((p - deltapl)**2 + deltapt**2)  
		IF(Pnew .gt.p*(1.0001))then
		write(2,*)'increase p',pnew,p
			endif
		p = pnew
		IF	(p .LE. 5.0) THEN
		px = 0.0
                py = 0.0
                pz = 0.0
                GOTO 30010
		ELSE
			!unit vector in new direction
!                       call ranlux(randvw,1)
        randvw(1) = RAN1(irnarg)
		sran = randvw(1) * 2.0 * pigreek
		cphir = cos(sran)
                sphir = sin(sran)
		sthetas = min (1.0, deltapt / p)
		cthetas = dsqrt(1.0 - sthetas**2) * forward
		ctheta = oldk
		stheta = dsqrt(1.0 - ctheta**2)
		IF (stheta .GT. 1.0E-10) THEN												
		sphi = oldi / stheta
                cphi = oldj / stheta
		ELSE
		sphi = 0.0
                cphi = 1.0														  
		ENDIF
		px = p * (cphi * sthetas * sphir + sphi * ctheta
     &               * sthetas * cphir  + sphi * stheta * cthetas)
		py = p * (- sphi * sthetas * sphir + cphi * ctheta
     &               * sthetas * cphir  + cphi * stheta * cthetas)
		pz = p * (- stheta * sthetas * cphir +
     &		 ctheta * cthetas)
		ke = dsqrt(p**2+elmsmass**2)-elmsmass !jan 2004***********************************
		IF(dabs((px**2+py**2+pz**2)/p**2 - 1.0) .gt. 5.0E-9)
     &		THEN
		write(2,*) "fail 1", dabs((px**2+py**2+pz**2)/p**2 - 1.0) 
		ENDIF     !REcheck at E-10 level
		IF(abs((oldi*px+oldj*py+oldk*pz)/p - cthetas) .gt.
     &		 5.0E-9) stop "fail 2"   !Rechecks at E-10 level
		ENDIF
30000        CONTINUE
! ---------------------------------------------------------
30010        pp(1)= px/1000
                pp(2)= py/1000
                pp(3)= pz/1000
                call VMAG(pp,pm)
                pavg=pm
                e = SQRT(pm**2 + pmass**2)
                ekin=ke             
800   continue
805   continue
900   continue
      return
      end
! *********************************************************
      subroutine ELMS_INIT
!
! Set up arrays required by ELMS
!
! Based on a routine by Simon Holmes, 22 NOV 2004
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
      integer elb,plb,ne3,np3
      logical done
      character*100 runnumber
      character*50 elmsdir,elmsroot
!
c read in DB on first run only
c
      elmsmass = 105.658
!
!      myownifrun = 0
      OPEN(UNIT=12, file='elmscom.txt',status='unknown',iostat=ioc)
      if( ioc .ne. 0 ) then
        stop 'Couldnt open ELMSCOM.TXT'
      end if
      read(12,'(a50)') elmsdir
      read(12,'(a50)') elmsroot
!      read (12,1368) ELMSCOR !correlations on (1) or off (0)
!      read (12,1369) jint !random seed
      CLOSE (12)
!      jint=234561
!      call rluxgo(4,jint,0,0)
!         OPEN(UNIT=12, file='ELMSDB/RunDirectory.txt')
      open(unit=12,file=elmsdir(1:LASTNB(elmsdir)),
     &  status='unknown',iostat=ioc)
      if( ioc .ne. 0 ) then
        write(*,'(a,a)') 'ELMSDIR = ',elmsdir
        stop 'Couldnt open ELMS run directory'
      end if
!         elmsfirst = .true.
!
      DO 30600 nmtm = 1,80
      calculated(nmtm) = 0			
	iF (nmtm .EQ. 1)THEN
	  firstpoint(1) = 1
	ELSE
	  firstpoint(nmtm) = lastpoint(nmtm-1) + 1
	  lastpoint(nmtm) = lastpoint(nmtm-1)
	ENDIF
      READ (12, 1335) mrun, mtmdb(nmtm), thkdb(nmtm), 
     &       colprobperm(nmtm), meanex(nmtm), meanp2x(nmtm)  			     
	L = nmtm / 10
      K =  nmtm - 10 * L
      J = L/10
      L = L - 10 * J 		
!                runnumber =  'ELMSDB/ELMSFv3run' // char(48 + J) // 
!     &               char(48 + L) // char(48+K) // '.elm'
      runnumber = elmsroot(1:LASTNB(elmsroot)) // char(48 + J) // 
     &       char(48 + L) // char(48+K) // '.elm'
	INQUIRE(FILE = runnumber, EXIST = done )
	IF (.NOT. done) GOTO 30600
	calculated(nmtm) = 1
      OPEN (UNIT=13, ACCESS = 'SEQUENTIAL',file=runnumber,iostat=ioc)
      if( ioc .ne. 0 ) then
        write(*,'(a,a)') 'ELMS file = ',runnumber
        stop 'Couldnt open ELMS data file'
      end if
      LL = firstpoint(nmtm)
	READ(13,7832)IL, a1, a2, factorp, factore, pscale(2),
     &    escale(2),ELB, ne3, PLB, np3
      IF (IL.NE.nmtm) STOP 'ELMS db check fail'
      DO I = 2, maxe-1
	escale(I+1) = escale(I) * factore
	pscale(I+1) = pscale(I) * factorp
	ENDDO
      pscale(PLB) = 0.0
      escale(ELB) = 0.0
	energy = dsqrt(mtmdb(nmtm)**2 + elmsmass**2)
	a1 = 0.0
      a2 = 0.0
      DO J=ELB,ne3
!read actual map. Store nonzero values as *10**8 integers, integrated
!store mean pl (a1), mean de (a2), prob lostas being unphysical   
	READ(13,2335)(map3(I),I=PLB,np3)
      DO 30700 I=PLB,np3
	IF (map3(I) .LE. 1.0D-10) GOTO 30700
	IF (LL .EQ. firstpoint(nmtm)) THEN
	  probarray(LL) = map3(I)
	ELSE
	  probarray(LL) = map3(I) + probarray(LL-1)
	ENDIF
	elmsde(LL) = 0.5*(escale(J)+escale(J+1))
	dptelms(LL) = 0.5 * (pscale(I) + pscale(I+1))
	lastpoint(nmtm) = LL
!check physical
	a4 = (energy - elmsde(LL))**2 - dptelms(LL)**2 - elmsmass**2  !pl squared						   
	IF(a4.le.0.0)THEN	  !unphysical combination of pt and de. Note if sizeable	   
	  IF(map3(I).GT.1.0E-9) write(2,*) 'unphysical', elmsde(LL),
     &       dptelms(LL), a4, map3(I)								   
	  GOTO 30700
	ENDIF														   
!accumulate checks
	a2 = a2 + elmsde(LL) * map3(I)
	a1 = a1 + dptelms(LL)**2 * map3(I)
	LL = LL + 1					
30700 CONTINUE
      ENDDO
      CLOSE (13)
	normelms(nmtm) = probarray(lastpoint(nmtm))
	a1 = a1 / (meanp2x(nmtm) * thkdb(nmtm)) 
	a2 = a2 / (meanex(nmtm) * thkdb(nmtm))
	N = lastpoint(nmtm)-firstpoint(nmtm)+1
	write(2,7833) nmtm, mtmdb(nmtm), lastpoint(nmtm),
     &  normelms(nmtm), a2, a1
30600 CONTINUE
      nmtm = 138
      CLOSE(12)
1335  Format(I4,F13.3,E13.6,3E17.9)
7832  Format(I5,2E12.5,2F6.3,2E10.2,4I4)
7833  Format(I4, ' Mom', F8.0, I8, ' Norm', F10.7,
     &   ' Ratio E, Pt2', 2F8.5)    
2335  Format(10F13.10)
1334  Format(12E15.8)
1368  Format(I1)
1369  Format(I6)
!
!Find the value of "eslope" and "pslope" for each momentum for interpolation
      J = 0
      DO I = 1, nmtm
	IF (calculated(I).NE.1) CYCLE
      IF (J .NE. 0) THEN
	  eslope(I) = (meanex(I)/meanex(J) -1.0)/(mtmdb(I) - mtmdb(J))
	  pslope(I) = (dsqrt(meanp2x(I)/meanp2x(J)) -1.0)
     &       /(mtmdb(I) - mtmdb(J))
      ENDIF
      J = I
	ENDDO
!vector of momenta actually run
      DO I = 1, nmtm
	IF (calculated(I).NE.1)CYCLE
	npdone = npdone + 1
	pdone(npdone) = I
	ENDDO
!find max interpolations
	DO I = 1, npdone
	IF(I.EQ.1)CYCLE
      a1=0.5*eslope(pdone(I))*(mtmdb(pdone(I))-mtmdb(pdone(I-1)))
	a2=0.5*pslope(pdone(I))*(mtmdb(pdone(I))-mtmdb(pdone(I-1)))
	write(2,7145)mtmdb(pdone(I)),a1,a2
!	write(8,7145)mtmdb(pdone(I)),a1,a2
7145 	Format(' Max mtm interpol at ',F6.0,' MeV in de,dpt',2F6.3)
	ENDDO
!
      end
c  *********************************************************
      subroutine FANO(omega,chic2,pavg,th)
c
c  Determine polar scattering angle from Fano distribution
!
!   Ref. U. Fano, Phys. Rev. 93:117-120, 1954
!   Derived from GEANT program  GMOLIE
c
C.    *                                                                *
C.    * THRED(NA)=reduced angles of Moliere theory                     *
C.    *                                                                *
C.    * F0I(NA),F1I(NA),F2I(NA)= integrals of Moliere functions        *
C.    *                                                                *
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      PARAMETER (ENEPER = 2.7182818d0) 
      DIMENSION TINT(40),ARG(4),VAL(4),THRED(40),F0I(40),F1I(40),F2I(40)
      DATA THRED/
     +   0.00, 0.10, 0.20, 0.30
     +,  0.40, 0.50, 0.60, 0.70
     +,  0.80, 0.90, 1.00, 1.10
     +,  1.20, 1.30, 1.40, 1.50
     +,  1.60, 1.70, 1.80, 1.90
     +,  2.00, 2.20, 2.40, 2.60
     +,  2.80, 3.00, 3.20, 3.40
     +,  3.60, 3.80, 4.00, 5.00
     +,  6.00, 7.00, 8.00, 9.00
     +, 10.00,11.00,12.00,13.00/
      DATA F0I/
     +  0.000000E+00 ,0.995016E-02 ,0.392106E-01 ,0.860688E-01
     + ,0.147856E+00 ,0.221199E+00 ,0.302324E+00 ,0.387374E+00
     + ,0.472708E+00 ,0.555142E+00 ,0.632121E+00 ,0.701803E+00
     + ,0.763072E+00 ,0.815480E+00 ,0.859142E+00 ,0.894601E+00
     + ,0.922695E+00 ,0.944424E+00 ,0.960836E+00 ,0.972948E+00
     + ,0.981684E+00 ,0.992093E+00 ,0.996849E+00 ,0.998841E+00
     + ,0.999606E+00 ,0.999877E+00 ,0.999964E+00 ,0.999990E+00
     + ,0.999998E+00 ,0.999999E+00 ,0.100000E+01 ,0.100000E+01
     + ,0.100000E+01 ,0.100000E+01 ,0.100000E+01 ,0.100000E+01
     + ,1.,1.,1.,1./
      DATA F1I/
     +  0.000000E+00,0.414985E-02,0.154894E-01,0.310312E-01
     + ,0.464438E-01,0.569008E-01,0.580763E-01,0.468264E-01
     + ,0.217924E-01,-0.163419E-01,-0.651205E-01,-0.120503E+00
     + ,-0.178272E+00,-0.233580E+00,-0.282442E+00,-0.321901E+00
     + ,-0.350115E+00,-0.366534E+00,-0.371831E+00,-0.367378E+00
     + ,-0.354994E+00,-0.314803E+00,-0.266539E+00,-0.220551E+00
     + ,-0.181546E+00,-0.150427E+00,-0.126404E+00,-0.107830E+00
     + ,-0.933106E-01,-0.817375E-01,-0.723389E-01,-0.436650E-01
     + ,-0.294700E-01,-0.212940E-01,-0.161406E-01,-0.126604E-01
     + ,-0.102042E-01,-0.840465E-02,-0.704261E-02,-0.598886E-02/
      DATA F2I/
     +  0.0,0.121500E-01,0.454999E-01,0.913000E-01
     + ,0.137300E+00,0.171400E+00,0.183900E+00,0.170300E+00
     + ,0.132200E+00,0.763000E-01,0.126500E-01,-0.473500E-01
     + ,-0.936000E-01,-0.119750E+00,-0.123450E+00,-0.106300E+00
     + ,-0.732800E-01,-0.312400E-01,0.128450E-01,0.528800E-01
     + ,0.844100E-01,0.114710E+00,0.106200E+00,0.765830E-01
     + ,0.435800E-01,0.173950E-01,0.695001E-03,-0.809500E-02
     + ,-0.117355E-01,-0.125449E-01,-0.120280E-01,-0.686530E-02
     + ,-0.385275E-02,-0.231115E-02,-0.147056E-02,-0.982480E-03
     + ,-0.682440E-03,-0.489715E-03,-0.361190E-03,-0.272582E-03/
*
*     ------------------------------------------------------------------
*
      IF(OMEGA.LE.ENEPER)GO TO 90
      CHIC  = SQRT(chic2)
      CNST=LOG(OMEGA)         ! Moliere b parameter
      bet = pavg/pmass / SQRT(1.+(pavg/pmass)**2)
      z = mparm(1,krad)
      if( NINT(z) .eq. 1 ) then
        uin = 3.6       ! Fano p.118
      else if( NINT(z) .eq. 2 ) then
        uin = 4.2       ! numerical integration of Hubbell S(x) data
      else if( NINT(z) .eq. 3 ) then
        uin = 4.6       ! Fano p.118
      else if( NINT(z) .eq. 4 ) then
        uin = 4.3       ! numerical integration of Hubbell S(x) data
      else              ! Z > 4
        uin = 5.   
      end if
      if( mt .eq. 'LIH' ) then
        uin = 4.1       ! use average of Li and H
      end if
!
      if( ABS(iptyp) .eq. 1 ) then        ! electron
        argm = 0.160 * z**(-0.67) * (1.+3.87e-3*z/bet)
        db = 1./(z+1.)*(LOG(argm)+uin)
      else                                ! heavy particle
        argm = 1130.*z**(-1.33)/(1./bet**2-1.)
        db = 1./z*(LOG(argm)+uin-0.5*bet*bet)
      end if
      cnst = cnst + db
! ---------------------------------------------------------
      B=5.
c
      DO 10 L=1,10
         IF(ABS(B).LT.1.E-10)THEN
            B=1.E-10
         ENDIF
         DB=-(B-LOG(ABS(B))-CNST)/(1.-1./B)
         B=B+DB
         IF(ABS(DB).LE.0.0001)GO TO 20
   10 CONTINUE
c
      GO TO 90
c
   20 CONTINUE
      IF(B.LE.0.)GO TO 90
      BINV = 1./B
      TINT(1) = 0.
c
      DO 30 JA=2,4
         TINT(JA)=F0I(JA)+(F1I(JA)+F2I(JA)*BINV)*BINV
   30 CONTINUE
c
      NMAX = 4
   40 CONTINUE
      XINT = RAN1(irnarg)
c
      DO 50 NA=3,40
         IF(NA.GT.NMAX) THEN
            TINT(NA)=F0I(NA)+(F1I(NA)+F2I(NA)*BINV)*BINV
            NMAX=NA
         ENDIF
         IF(XINT.LE.TINT(NA-1)) GO TO 60
   50 CONTINUE
c
      IF(XINT.LE.TINT(40)) THEN
         NA=40
         GOTO 60
      ELSE
         TMP=1.-(1.-B*(1.-XINT))**5
         IF(TMP.LE.0.)GO TO 40
         THRI=5./TMP
         GO TO 80
      ENDIF
c
   60 CONTINUE
      NA = MAX(NA-1,3)
      NA3 = NA-3
c
      DO 70 M=1,4
         NA3M=NA3+M
         ARG(M)=TINT(NA3M)
         VAL(M)=THRED(NA3M)**2
   70 CONTINUE
c
      call POLINT(arg,val,4,xint,thri,err)    ! Num Rec
c
   80 CONTINUE
      TH=CHIC*SQRT(ABS(B*THRI))
      IF(TH.GT.PI)GO TO 40
      SINTH=SIN(TH)
      TEST=TH*(RAN1(irnarg))**2
      IF(TEST.GT.SINTH)GO TO 40
      GOTO 900
!
   90 CONTINUE
      iflag = -24
c
900   continue
      END
! *********************************************************
      subroutine FLIP_COULOMB(p,h,pr)
!
! Compute probability of muon spin flip in mu-Z Coulomb scattering
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      real*8 me,mu
!
      data con/4.375e-8/,alp/7.299e-3/
!
      me = mass(1)            ! [GeV]
      mu = mass(2)
      e = SQRT( p*p + mu*mu )
      b = p / e
      z = mparm(1,krad)
      a = aelcomp(1,krad)
      r = mparm(2,krad)       ! [g cm^-3]
      thmin = alp*z**(0.333333)*me/p
      thmax = 2./alp/a**(0.333333)*me/p
      arg = (1.-COS(thmax)) / (1.-COS(thmin))
!
      pr = con*h*r/a*(z/b/p/e)**2 * LOG( arg )
      hel = HELICITY( pol )
      pr = pr * ABS( hel )
!
      end
! *********************************************************
      subroutine FLIP_ELASTIC(s,z,pr)
!
! Compute probability of muon spin flip in mu-e elastic scattering
! Adapted from J. Gallardo 10 April 2000
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      dimension s1(4),s2(4),q1(4),q2(4),p1(4),p2(4)
!
!      xme=0.51099907d-3  ! mass in GeV
!      xmu=105.658389d-3
      xme = mass(1)
      xmu = mass(2)
!
      if( ABS(z) .le. 1.d0 ) then
        sinz=sqrt((1.d0-z**2))
      else
        iflag = -64
        go to 900
      end if
!
c cm of momentum frame of reference
      emu1cm=0.5*(s+xmu**2-xme**2)/sqrt(s)
      qmu1cm=sqrt(emu1cm**2-xmu**2)
      t=-2.*qmu1cm**2*(1.d0-z)
!
      q1(1)=emu1cm
      q1(2)= 0.d0
      q1(3)= 0.d0
      q1(4)=qmu1cm
      p1(1)=0.5*(s+xme**2-xmu**2)/sqrt(s)
      p1(2)= 0.d0
      p1(3)= 0.d0
      p1(4)=-qmu1cm
c
      q2(1)=emu1cm
      q2(2)=sinz*qmu1cm
      q2(3)= 0.d0
      q2(4)=z*qmu1cm
      p2(1)=0.5*(s+xme**2-xmu**2)/sqrt(s)
      p2(2)=-sinz*qmu1cm
      p2(3)= 0.d0
      p2(4)=-z*qmu1cm
!
      s1(1)=qmu1cm/xmu
      s1(2)= 0.d0
      s1(3)= 0.d0
      s1(4)=q1(1)/xmu
      s2(1)=qmu1cm/xmu
      s2(2)=emu1cm*sinz/xmu
      s2(3)= 0.d0
      s2(4)=emu1cm*z/xmu
!
      denominator = 32*((s-(xme**2+xmu**2))**2+0.5*t**2+s*t)
c
      xnumerator = -32*dot(s1,s2)*((s-(xmu**2+xme**2))**2+s*t)+ 
     +  64*s*(dot(p1,s2)*dot(q2,s1)+dot(p2,s1)*dot(q1,s2)+
     +  dot(q1,s2)*dot(q2,s1))+
     +  64*t*(dot(p1,s2)*dot(p2,s1)+dot(p1,s2)*dot(q2,s1)+
     +  dot(p2,s1)*dot(q1,s2)+dot(q1,s2)*dot(q2,s1))-
     +  64*(dot(p1,s2)*dot(q2,s1)*(xmu**2+xme**2)+dot(p2,s1)*
     +  dot(q1,s2)*(xmu**2+xme**2)+
     +  (xmu**2+2*xme**2)*dot(q1,s2)*dot(q2,s1))
!
      polar=xnumerator/denominator  ! polarization
      pr =0.5*(1.0d0-polar) ! Prob(up --->down)
      f = 0.5d0*(1.0d0 - z)   ! spin flip vs helicity flip
      pr = f - pr
      hel = HELICITY( pol )
      pr = pr * ABS( hel )
!
900   continue
      end
! *********************************************************
      REAL*8 FUNCTION FNEAR(C)
!
! Auxiliary routine for SAMCS
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)
      common/rans/iseed
!      REAL RNDM
! 1    RNF=DBLE(RNDM(-1.))
 1    RNF = RAN1(iseed)
      F1=1.D0-C/(1.D0+2.D0*C)
      IF(RNF.LE.F1) THEN
!	 RNF=DBLE(RNDM(-1.))
	 RNF = RAN1(iseed)
	 X=-LOG(RNF)
      ELSE
!	 RNF1=DBLE(RNDM(-1.))
!	 RNF2=DBLE(RNDM(-1.))
!	 RNF3=DBLE(RNDM(-1.))
	 RNF1 = RAN1(iseed)
	 RNF2 = RAN1(iseed)
	 RNF3 = RAN1(iseed)
	 X=-LOG(RNF1*RNF2*RNF3)
      ENDIF
      F2=1.-2.D0*C*X/(1.D0+C+C*X**2/2.D0)
!      RNF=DBLE(RNDM(-1.))
      RNF = RAN1(iseed)
      IF(RNF-F2) 2,1,1
 2    FNEAR=X  
      RETURN
      END 
! *********************************************************
	DOUBLE PRECISION FUNCTION FORMN(QG,Z,A)
!
! Auxiliary routine for SAMCS
!
C      ------------------------------------------------------------------------------------
C       Electromagnetic form factor of target atom squared 
C      ------------------------------------------------------------------------------------
C       QG - momentum transfer in GeV/c
C      -------------------------------------------------------------------------------------
C       Z,A - charge and atomic mass of target  
C      -------------------------------------------------------------------------------------
C       Data from:
C        H. de Vries, C.W. de Jager and C. de Vries
C        Atomic Data and Nuclear Data Tables 36 (1987) 495-536
C
C       V.V. Burov et al.  Phys. Atom. Nucl. 61 (1998) 526-532
C  -----------------------------------------------------------------
        IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)
	DIMENSION AHO(8)
        PARAMETER (SC=5.06765D0)
	PARAMETER (PI=3.141592653589793227D+00)
	DATA AHO/0.85d0,0.,1.77d0,1.78d0,1.69d0,0.,1.745d0,1.833d0/
!
        Q=QG*SC

C  Gaussian formfactor - for comparizon with analitical calculation
c        if(z.le.300.d0) then
c        r2=rav2(z,a) 
c        formn=exp(-q**2*r2/3.d0)
c        return
c        endif  

	IF(Z.EQ.1.D0) THEN
        FORMN = 1.D0/(1.D0+(Q*AHO(1))**2/12.D0)**4 
        return
        endif
      
	IF(Z.eq.2.D0) then
        R=1.238
        B=0.38
	FORMN = ZFORM(Q,R,B)**2
        return
        endif

	IF(Z.eq.6.D0) then
        R=2.214
        B=0.488
	FORMN = ZFORM(Q,R,B)**2
        return
        endif

        IF(Z.LE.8.D0) then
        IZ=Z
        B2=AHO(IZ)**2
        A0=(B2-.4824)*A/(A-1.)
        Z1=1.-(Z-2.)*A0*Q**2/6./Z
	FORMN=(Z1*EXP(-B2*Q**2/4.))**2
        ELSE 
        A13=A**(1./3.)
        R=1.141*A13-1.047/A13
        R2=(0.898*A13+1.198/A13)**2
	B=R*SQRT(3.*(5.*R2/3./R**2-1.)/7.)/PI
	FORMN = ZFORM(Q,R,B)**2
        ENDIF
	RETURN
	END
! *********************************************************
	DOUBLE PRECISION FUNCTION FORMP(QG,IP)
!
! Auxiliary routine for SAMCS
!
C      ------------------------------------------------------------------------------------
C       Electromagnetic form factor of projectile squared 
C      ------------------------------------------------------------------------------------
        IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)

C ---------------  GeV to 1/fermi
        Q=QG*5.06765D0

        IF(IP.LE.13) THEN
C --------------- Dipole form factor for hadrons
           FORMP = 1.D0/(1.D0+Q**2*R2P(IP)/12.D0)**4 
        ELSE
C --------------- Gaussian form factor for light nuclei
           FORMP=EXP(-Q**2*R2P(IP)/3.D0)
        ENDIF
        RETURN
        END
! *********************************************************
      SUBROUTINE GDECA2(XM0,XM1,XM2,PCM)
C.
C.    ******************************************************************
C.    *                                                                *
C.    *  Simulates two body decay process with isotropic angular       *
C.    *  distribution in CMS.                                          *
C.    *                                                                *
C.    *       Author    G.Patrick  *********                           *
C.    *                                                                *
C.    ******************************************************************
C.
! Adapted from Geant 3.2.1
!
      implicit double precision (a-h,o-z)
      DIMENSION PCM(4,2)
      DIMENSION RNDM(2)
C.
      data twopi/6.2831853d0/
C.    ------------------------------------------------------------------
C.
C             Generate first decay product in CMS.
C
      E1=(XM0*XM0+XM1*XM1-XM2*XM2)/(2d0*XM0)
      P1=SQRT(ABS((E1-XM1)*(E1+XM1)))
C
C             Isotropic decay angular distribution.
C
      CALL GRNDM(RNDM,2)
      COSTH=2d0*RNDM(1)-1d0
      IF(ABS(COSTH).GE.1d0) THEN
         COSTH=SIGN(1.d0,COSTH)
         SINTH=0d0
      ELSE
         SINTH=SQRT((1d0-COSTH)*(1d0+COSTH))
      ENDIF
      PHI=TWOPI*RNDM(2)
!      sinth=0.173648      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!      costh=0.984808
!      phi=0.
C
C             Polar co-ordinates to momentum components.
C
      PCM(1,1)=P1*SINTH*COS(PHI)
      PCM(2,1)=P1*SINTH*SIN(PHI)
      PCM(3,1)=P1*COSTH
      PCM(4,1)=E1
C
C             Generate second decay product.
C
      PCM(1,2)=-PCM(1,1)
      PCM(2,2)=-PCM(2,1)
      PCM(3,2)=-PCM(3,1)
      PCM(4,2)=SQRT(PCM(1,2)**2+PCM(2,2)**2+PCM(3,2)**2+XM2*XM2)
C
      END
! *********************************************************
      SUBROUTINE GDECA3(XM0,XM1,XM2,XM3,PCM)
C.
C.    ******************************************************************
C.    *                                                                *
C.    *   This routine simulates three-body decays in center-of-mass.  *
C.    *   Written by Jurgen Schukraft and Vincent Hedberg for R807.    *
C.    *   Adapted for Geant3 by Sverker Johansson (Aug-1987).          *
C.    *                                                                *
C.    *   As it is, the decay is done according to pure phase-space.   *
C.    *   There is code in the routine to properly treat K decays (V-A)*
C.    *   but it is not used. (Requires an additional parameter        *
C.    *   in the call, to tell which type of decay.)                   *
C.    *                                                                *
C.    *       Author    S. Johansson  *********                        *
C.    *                                                                *
C.    ******************************************************************
C.
! Adapted from Geant 3.2.1
!
      implicit double precision (a-h,o-z)
      DIMENSION  P(20,20)
      DIMENSION  HM(4,7)
      DIMENSION  PCM(4,3)
      DIMENSION RNDM(3)
*
****************************************************
*   Translation from new to old input parameters :
*      ( S.J. )
****************************************************
      ID = 1
      HM(1,ID) = XM1
      HM(2,ID) = XM2
      HM(3,ID) = XM3
      HM(4,ID) = XM0
      K1 = 1
      K2 = 2
      K3 = 3
      PI2 = 2d0 * 3.141592654d0
***********************
*  Original code :
C--------------------------------------------------------
C     SIMULATION OF SECONDARY PARTICLE MOMENTA
C     IN THE REST FRAME OF PARENT PARTICLE
C--------------------------------------------------------
      T0=HM(4,ID)-HM(1,ID)-HM(2,ID)-HM(3,ID)
  10  CALL GRNDM(RNDM,2)
      G1=RNDM(1)
      G2=RNDM(2)
      GMI=MIN(G1,G2)
      GMA=MAX(G1,G2)
      TE=GMI*T0
      TM=(1d0-GMA)*T0
      TN=(GMA-GMI)*T0
      PA1=SQRT(TE**2+2d0*HM(1,ID)*TE)
      PA2=SQRT(TM**2+2d0*HM(2,ID)*TM)
      PA3=SQRT(TN**2+2d0*HM(3,ID)*TN)
C--------------------------------------------------------
C     CHECK OF THE MOMENTUM CONSERVATION
C--------------------------------------------------------
      PMMM=MAX(PA1,PA2)
      PMAX=MAX(PA3,PMMM)
      SP=PA1+PA2+PA3
      IF(PMAX.GE.SP-PMAX) GOTO 10
      P(4,K1)=TE+HM(1,ID)
      P(4,K2)=TM+HM(2,ID)
      P(4,K3)=TN+HM(3,ID)
C--------------------------------------------------------
C     ACCOUNTING FOR MATRIX ELEMENT (FOR ID=4 AND 5)
*     (** Calculates correct (V-A) decay for kaons.       **)
*     (** Now inactive.  Change ID to activate it. (S.J)  **)
C--------------------------------------------------------
      IF(ID.NE.4.AND.ID.NE.5) GOTO 20
      EEM=((HM(4,ID)-HM(2,ID))**2+HM(1,ID)**2-
     *      HM(3,ID)**2)/(2d0*(HM(4,ID)-HM(2,ID)))
      HMMAX=EEM*(HM(4,ID)-HM(2,ID)-EEM)+
     *      0.5d0*(EEM**2-HM(1,ID)**2)
      HMA=P(4,K1)*P(4,K3)+0.25d0*(PA1**2+PA3**2-PA2**2)
      CALL GRNDM(RNDM,1)
      GG=RNDM(1)
      IF(HMA.LT.GG*HMMAX) GOTO 10
C--------------------------------------------------------
C     CALCULATIONS OF MOMENTA
C--------------------------------------------------------
20    CONTINUE
      CALL GRNDM(RNDM,3)
      CTE=2d0*RNDM(1)-1d0
      STE=SQRT(ABS(1d0-CTE**2))
      FE =PI2*RNDM(2)
      CFE=COS(FE)
      SFE=SIN(FE)
      P(1,K1)=PA1*STE*CFE
      P(2,K1)=PA1*STE*SFE
      P(3,K1)=PA1*CTE
      CTEN=(PA2**2-PA1**2-PA3**2)/(2d0*PA1*PA3)
      STEN=SQRT(ABS(1d0-CTEN**2))
      FEN =PI2*RNDM(3)
      CFEN=COS(FEN)
      SFEN=SIN(FEN)
      P(1,K3)=PA3*(STEN*CFEN*CTE*CFE-STEN*SFEN*SFE+CTEN*STE*CFE)
      P(2,K3)=PA3*(STEN*CFEN*CTE*SFE+STEN*SFEN*CFE+CTEN*STE*SFE)
      P(3,K3)=PA3*(-STEN*CFEN*STE+CTEN*CTE)
      P(1,K2)=-P(1,K1)-P(1,K3)
      P(2,K2)=-P(2,K1)-P(2,K3)
      P(3,K2)=-P(3,K1)-P(3,K3)
C-------------------------------------------------------
***   Re-translation from old to new output arrays : ***
***     (S.J.)                                       ***
***--------------------------------------------------***
      DO 37 I=1,4
        DO 38 J=1,3
          PCM(I,J) = P(I,J)
 38     CONTINUE
 37   CONTINUE
      END
c  *********************************************************
      SUBROUTINE GDROT(P,COSTH,SINTH,COSPH,SINPH)
C.
C.    ******************************************************************
C.    *                                                                *
C.    *  Rotates vector from one reference system into another.        *
C.    *  THETA and PHI are anti-clockwise Eulerian angles between the  *
C.    *  two systems.                                                  *
C.    *                                                                *
C.    *       Author    M.Hansroul, G.Patrick  *********               *
C.    *                                                                *
C.    *                                                                *
C.    ******************************************************************
C.
      implicit double precision (a-h,o-z)
      DIMENSION P(3)
C.
C.    ------------------------------------------------------------------
C.
      P1=P(1)
      P2=P(2)
      P3=P(3)
      P(1)=P1*COSTH*COSPH - P2*SINPH + P3*SINTH*COSPH
      P(2)=P1*COSTH*SINPH + P2*COSPH + P3*SINTH*SINPH
      P(3)=-P1*SINTH                 + P3*COSTH
C
      END
! *********************************************************
      SUBROUTINE GLOREN(BETA,PA,PB)
C.
C.    ******************************************************************
C.    *                                                                *
C     *       Routine to transform momentum and energy from the        *
C     *       Lorentz frame A to the Lorentz frame B                   *
C     *                                                                *
C     *       PA(1)                                                    *
C     *       PA(2)     Momentum components in frame A                 *
C     *       PA(3)                                                    *
C     *       PA(4)     Energy                                         *
C     *       PB(..)   same quantities in frame B                      *
C     *                                                                *
C     *       BETA(1)    Components of velocity of frame B             *
C     *       BETA(2)        as seen from frame A                      *
C     *       BETA(3)                                                  *
C     *       BETA(4)    1./SQRT(1.-BETA**2)                           *
C.    *                                                                *
C.    *       Author    M.Hansroul  *********                          *
C.    *                                                                *
C.    ******************************************************************
C.
      implicit double precision (a-h,o-z)
      DIMENSION BETA(4),PA(4),PB(4)
C.
C.    ------------------------------------------------------------------
C.
!      gam = beta(4)
!      bg = beta(3)*beta(4)
!      pb(1) = pa(1)
!      pb(2) = pa(2)
!      pb(3) = gam*pa(3) - bg*pa(4)
!      pb(4) = gam*pa(4) - bg*pa(3)
!      return
!
      BETPA  = BETA(1)*PA(1) + BETA(2)*PA(2) + BETA(3)*PA(3)
      BPGAM  = (BETPA * BETA(4)/(BETA(4) + 1d0) - PA(4)) * BETA(4)
      PB(1) = PA(1) + BPGAM  * BETA(1)
      PB(2) = PA(2) + BPGAM  * BETA(2)
      PB(3) = PA(3) + BPGAM  * BETA(3)
      PB(4) =(PA(4) - BETPA) * BETA(4)
      END
c  *********************************************************
      SUBROUTINE GPOISS(AMVEC,NPVEC,LEN)
!
! Poisson random number interface for GEANT routines
!
C.
C.    ******************************************************************
C.    *                                                                *
C.    *  Generates  a vector NPVEC of LEN random numbers               *
C.    *  POISSON distribued with mean values AMVEC                     *
C.    *                                                                *
C.    *  If the mean value A greater than PLIM, N is calculated        *
C.    *  according to the Gaussian approximation of the Poisson        *
C.    *  distribution.                                                 *
C.    *                                                                *
C.    *                                                                *
C.    *  Author      : L.Urban                                         *
C.    *  Date        : 28.04.1988           Last update :  1.02.1990   *
C.    *                                                                *
C.    ******************************************************************
C.
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      REAL*8 AMVEC(*),RNDM(2), N
      INTEGER NPVEC(*)
      PARAMETER (PLIM=16.,HMXINT=2E+9)
!
!      do i=1,len
!      xm = amvec(i)
!      npvec(i) = POIDEV(xm,idum)
!      end do
!
      DO 30 I=1,LEN
*     Protection against negative mean values
         N=0.
         IF(AMVEC(I).GT.0.) THEN
            IF(AMVEC(I).LE.PLIM) THEN
               CALL GRNDM(RNDM,1)
               R=RNDM(1)
               P=EXP(-AMVEC(I))
               S=P
               IF(R.LE.S) GOTO 20
   10          N=N+1.
               P=P*AMVEC(I)/N
               S=S+P
               IF(S.LT.R.AND.P.GT.1.E-30) GOTO 10
            ELSE
               CALL GRNDM(RNDM,2)
               RR=SQRT(-2.*LOG(RNDM(1)))
               PHI=TWOPI*RNDM(2)
               X=RR*COS(PHI)
               N=MIN(MAX(AMVEC(I)+X*SQRT(AMVEC(I)),0.d0),HMXINT)
            ENDIF
         ENDIF
*
   20    NPVEC(I) = N
   30 CONTINUE
*
      END
c  *********************************************************
      subroutine GRNDM(rvec,len)
!
! Uniform random number interface for GEANT routines
!
      implicit double precision (a-h,o-z)
      real*8 rvec(*)
!
      do i=1,len
      rvec(i) = RANV(0.d0,1.d0)
      end do
!
      return
      end
c  *********************************************************
      real*8 function HELICITY(s)
!
! Compute current muon helicity in the LAB frame
!  s  is the current value of particle spin vector in muon REST frame
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      real*8 s(3),p1(3)
!
      call VMAG(pp,pm)
!
      do i=1,3
      p1(i) = pp(i) / pm       ! unit vector along p
      end do
!
      helicity = DOT3(s,p1)
!
      return
      end
! *********************************************************
      function ILVSDMI(a,n,inc)
!
! Find index of smallest element in array A(N)
! Replaces CERNLIB routine LVSDMI
! Used by SAMCS
!
      implicit double precision (a-h,o-z)
	dimension a(n)
!
      ilvsdmi = 1
	amax = 1d12
      do i=1,n,inc
	if( a(i) .lt. amax ) then
	  amax = a(i)
	  ilvsdmi = i
	end if
	end do
!
      end
! *********************************************************
      function ILVSDMX(a,n,inc)
!
! Find index of largest element in array A(N)
! Replaces CERNLIB routine LVSDMX
! Used by SAMCS
!
      implicit double precision (a-h,o-z)
	dimension a(n)
!
      ilvsdmx = 1
	amin = -1d12
      do i=1,n,inc
	if( a(i) .gt. amin ) then
	  amin = a(i)
	  ilvsdmx = i
	end if
	end do
!
      end
c  *********************************************************
      subroutine INTERACTION(hstp,pavg)
!
! Check for inelastic nuclear interactions of hadrons
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      if( ABS(iptyp) .lt. 3 ) go to 900   ! skip e and mu
      p4 = pavg 
      sig = CS_INEL(p4,iptyp) * 1e-27  ! get inelastic cross section
      if( sig .eq. 0. ) go to 900
      a = aelcomp(1,krad)
      xint = 0.01 / ( densat * sig * a**(0.71) )  ! [m]
      xm = hstp / xint
      nin = POIDEV( xm,irnarg )
      if( nin .ge. 1 ) then
        if( intlev.eq.2 .and. iptyp.eq.5 ) then
          call WANG
        else
          iflag = -73
        end if
      end if
!
900   continue
      end
c  *********************************************************
      logical function IN_ASPWEDGE(hstp)
!
!  Determine if particle step goes thru an
!   Azimuthally Symmetric Polynomial shaped wedge
!
!   Describe polynomial shape in 1st quadrant of coordinate system
!    centered on the wedge: r(z) = a0+a1*z+a2*z**2+a3*z**3
!
!  The wedge body has reflection symmetry thru x-y plane
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      logical in0,in1
!
      in_aspwedge = .false.
      if( mgt .ne. 'ASPW' ) go to 900
!
      in0 = .false.
      in1 = .false.
      dr0 = 0.
      dr1 = 0.
!
      zv = sr0 + gparm(1,krad)              ! z of wedge center
      zcut = gparm(2,krad)                  ! z offset from center to edge
      a0 = gparm(3,krad)
      a1 = gparm(4,krad)
      a2 = gparm(5,krad)
      a3 = gparm(6,krad)
!     
      rp = SQRT( xp(1)**2 + xp(2)**2 )
      zmin = zv - zcut
      zmax = zv + zcut
!
!  Check end point
!  n.b. xp contains coordinates at end of the time step 
      if( xp(3).gt.zmin .and. xp(3).lt.zmax ) then
        dz = ABS( xp(3) - zv )
        rsurf = a0 + a1*dz + a2*dz**2 + a3*dz**3
        dr1 = rp - rsurf
        if( dr1 .le. 0. ) in1 = .true.
      end if
!
!  Check starting point
      x1 = xp(1) - pp(1)/pp(3)*hdid
      y1 = xp(2) - pp(2)/pp(3)*hdid
      z1 = xp(3) - hdid
      rp = SQRT( x1*x1 + y1*y1 )
!
      if( z1.gt.zmin .and. z1.lt.zmax ) then
        dz = ABS( z1 - zv )
        rsurf = a0 + a1*dz + a2*dz**2 + a3*dz**3
        dr0 = rp - rsurf
        if( dr0 .le. 0. ) in0 = .true.
      end if
!
      if( in0 .and. in1 ) then
        in_aspwedge = .true.
!     
      else if( in0 .and. (.not. in1) ) then
        in_aspwedge = .true.
        hstp = hstp * ABS(dr0) / ( ABS(dr0)+ABS(dr1) )
!     
      else if( in1 .and. (.not. in0) ) then
        in_aspwedge = .true.
        hstp = hstp * ABS(dr1) / ( ABS(dr0)+ABS(dr1) )
      end if
!
900   continue
      return
      end
c  *********************************************************
      logical function IN_ASRWEDGE(hstp)
!
!  Determine if particle step goes thru an
!   Azimuthally Symmetric, Radial Wedge 
!   (defines thickness as a function of R)
!
!  Wedge is symmetric about xy plane at z=zv
!  Describe polynomial shape in terms of half-thickness 
!       (distance of edge from zv)
!    dz(r) = b0+b1*r+b2*r**2+b3*r**3
!    dz < 0 means no wedge at that radius
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      logical in0,in1
!
      in_asrwedge = .false.
      if( mgt .ne. 'ASRW' ) go to 900
!
      in0 = .false.
      in1 = .false.
      dr0 = 0.
      dr1 = 0.
!
      zv = sr0 + gparm(1,krad)              ! z of center of wedge
      zcut = gparm(2,krad)                  ! maximum allowed half-thickness
      b0 = gparm(3,krad)
      b1 = gparm(4,krad)
      b2 = gparm(5,krad)
      b3 = gparm(6,krad)
!     
!  Check end point
!  n.b. xp contains coordinates at end of the time step 
      rp = SQRT( xp(1)**2 + xp(2)**2 )
      dzsurf = b0 + b1*rp + b2*rp**2 + b3*rp**3
!  dzsurf <0 means no wedge at this radius
      if (dzsurf .lt. 0.) dzsurf = 0.
!  max thickness zcut
      if (dzsurf .gt. zcut) dzsurf = zcut
      zmin = zv - dzsurf
      zmax = zv + dzsurf
!
      if( xp(3).gt.zmin .and. xp(3).lt.zmax ) then
!  if inside, find out how far from left edge (in case crossing occurred)
        in1 = .true.
        dz1 = xp(3) - zmin 
      else
!  if outside, find out how far from right edge (in case crossing occurred)
        dz1 = zmax - xp(3) 
      end if
!
!  Check starting point
      x1 = xp(1) - pp(1)/pp(3)*hdid
      y1 = xp(2) - pp(2)/pp(3)*hdid
      z1 = xp(3) - hdid
      rp = SQRT( x1*x1 + y1*y1 )
      dzsurf = b0 + b1*rp + b2*rp**2 + b3*rp**3
!  dzsurf <0 means no wedge at this radius
      if (dzsurf .lt. 0.) dzsurf = 0.
!  max thickness zcut
      if (dzsurf .gt. zcut) dzsurf = zcut
      zmin = zv - dzsurf
      zmax = zv + dzsurf
!
      if( z1.gt.zmin .and. z1.lt.zmax ) then
!  if inside, find out how far from right edge (in case crossing occurred)
        in0 = .true.
        dz0 = zmax - z1 
      else
!  if outside, find out how far from left edge (in case crossing occurred)
        dz0 = z1 - zmin 
      end if
!
      if( in0 .and. in1 ) then
        in_asrwedge = .true.
!     
      else if( in0 .and. (.not. in1) ) then
        in_asrwedge = .true.
        hstp = hstp * ABS(dz0) / ( ABS(dz0)+ABS(dz1) )
!     
      else if( in1 .and. (.not. in0) ) then
        in_asrwedge = .true.
        hstp = hstp * ABS(dz1) / ( ABS(dz0)+ABS(dz1) )
      end if
!
900   continue
      return
      end
c  *********************************************************
      logical function IN_HEMIWIN()
!
!  Determine how a particle step proceeds thru a
!   hemispherical end region of an absorber
!
!  IHWIN returned as current ICOOL "radial" region
!  Returns IFLAG=1 if the step crosses a region boundary
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      real*8 dp(3),dz(3),x0(3),x(3)
!
      data eps/2.e-6/
!
      in_hemiwin = .false.
      if( mgt .ne. 'HWIN' ) go to 900
!
      i0 = 0      ! subregion of initial step position
      i1 = 0      ! subregion of final step position
! .........................................................
! Set geometry parameters
      rs1 = gparm(2,1)            ! inner spherical radius
      rs2 = rs1 + gparm(3,1)      ! outer spherical radius
      dzc = gparm(4,1)            ! offset of center of shell
      r3 = rhigh(1)               ! maximum cylindrical radius
      dz(1) = rs1 - dzc           ! absorber thickness
      dz(2) = gparm(3,1)          ! window thickness
      dz(3) = slen - rs2 + dzc    ! exterior thickness
!
      if( gparm(1,1) .lt. 0. ) then   ! ENTRANCE end
        zoff = slen
      else                            ! EXIT end
        zoff = 0.
      end if
!
! .........................................................     
! Check location of starting point
! Use coordinate system centered on the spherical shell
      x0(1) = xp0(1)
      x0(2) = xp0(2)
      x0(3) = zoff + dzc + gparm(1,1)*(xp0(3)-sr0)
      rsp0 = SQRT( x0(1)**2 + x0(2)**2 + x0(3)**2 )
!          
      if( rsp0 .lt. rs1 ) then             
        i0 = 1
      else if( rsp0.ge.rs1 .and. rsp0.le.rs2 ) then    
        i0 = 2
      else
        rp = SQRT( xp0(1)**2 + xp0(2)**2 ) 
        if( rp .lt. r3 ) then
          i0 = 3                             
        else
          go to 900
        end if
      end if   
! .........................................................
! Check location of end point
      x(1) = xp(1)
      x(2) = xp(2)
      x(3) = zoff + dzc + gparm(1,1)*(xp(3)-sr0)
      rsp = SQRT( x(1)**2 + x(2)**2 + x(3)**2 ) 
!         
      if( rsp .lt. rs1 ) then     
        i1 = 1
      else if( rsp.ge.rs1 .and. rsp.le.rs2 ) then 
        i1 = 2
      else
        rp = SQRT( xp(1)**2 + xp(2)**2 ) 
        if( rp .lt. r3 ) then
          i1 = 3       
        else
          go to 900
        end if
      end if   
! .........................................................
      if( iflag .eq. 1 ) then    ! this is the correction step
! Prepare stepping for the NEXT region
        ihwin = i0 + NINT( gparm(1,1) )
        hregn = dz(ihwin) / 10
        iflag = 0
        go to 800
      end if
!
      ihwin = i0              ! default values for current region
      hregn = dz(i0) / 10.
      iflag = 0
!
! ---------------------------------------------------------
! Step handling logic
! There are 18 cases to consider, including jumping over windows altogether
!
      if( gparm(1,1) .lt. 0. ) then   ! ENTRANCE end
        if( i0.eq.3 .and. i1.eq.3 ) then
          diff = ABS( rsp - rs2 )
          if( diff .le. eps ) then  ! step ends on 2nd sphere boundary
            ihwin = 2
            hregn = dz(2) / 10.     ! step size for NEXT step
          end if
          go to 800       ! good step
!
        else if( i0.eq.3 .and. i1.eq.2 ) then 
          diff0 = ABS( rsp0 - rs2 )
          diff = ABS( rsp - rs2 )
          if( diff0.le.eps .or. diff.le.eps ) then
! One end of the step ends on the 2nd sphere boundary
            ihwin = 2
            hregn = dz(2) / 10.
            if( diff0 .lt. eps ) i0 = 2 ! get right material
            go to 800       ! good step
          end if
          r2 = rs2**2       ! set radius for step adjustment logic
!
        else if( i0.eq.3 .and. i1.eq.1 ) then 
          diff0 = ABS( rsp0 - rs1 )
          if( diff0 .le. eps ) then   ! step starts on 2nd sphere boundary
            ihwin = 2
            hregn = dz(2) / 10.
            iflag = 1
            go to 900       ! redo step
          end if
          r2 = rs2**2       ! set radius for step adjustment logic
!
        else if( i0.eq.2 .and. i1.eq.2 ) then 
          diff = ABS( rsp - rs1 )
          if( diff .le. eps ) then  ! step ends on 1st sphere boundary
            ihwin = 1
            hregn = dz(1) / 10.     ! step size for NEXT step
          end if
          go to 800       ! good step
!
        else if( i0.eq.2 .and. i1.eq.1 ) then 
          diff0 = ABS( rsp0 - rs1 )
          diff = ABS( rsp - rs1 )
          if( diff0.le.eps .or. diff.le.eps ) then
! One end of the step ends on the 1st sphere boundary
            ihwin = 1
            hregn = dz(1) / 10.
            if( diff0 .lt. eps ) i0 = 1 ! get right material
            go to 800       ! good step
          end if
          r2 = rs1**2       ! set radius for step adjustment logic
!
        else if( i0.eq.1 .and. i1.eq.1 ) then 
          go to 800     ! good step
!
! It's also possible for trajectories to cross subregions backwards
        else if( i0.eq.1 .and. i1.eq.2 ) then 
          diff0 = ABS( rsp0 - rs1 )
          diff = ABS( rsp - rs1 )
          if( diff0.le.eps .or. diff.le.eps ) then
! One end of the step ends on the 1st sphere boundary
! Use the small stepping appropriate to region 2
            ihwin = 2
            hregn = dz(2) / 10.
            if( diff0 .lt. eps ) i0 = 2 ! use right material
            go to 800       ! good step
          end if
          r2 = rs1**2       ! set radius for step adjustment logic
!
        else if( i0.eq.2 .and. i1.eq.3 ) then 
          diff0 = ABS( rsp0 - rs2 )
          diff = ABS( rsp - rs2 )
          if( diff0.le.eps .or. diff.le.eps ) then
! One end of the step ends on the 2nd sphere boundary
            ihwin = 2
            hregn = dz(2) / 10.
            if( diff0 .lt. eps ) i0 = 3 ! use right material
            go to 800       ! good step
          end if
          r2 = rs2**2       ! set radius for step adjustment logic
!
        else if( i0.eq.1 .and. i1.eq.3 ) then 
          diff0 = ABS( rsp0 - rs1 )
          if( diff0 .le. eps ) then   ! step starts on 1st sphere boundary
            ihwin = 2
            hregn = dz(2) / 10.
            iflag = 1
            go to 900       ! redo step
          end if
          r2 = rs1**2       ! set radius for step adjustment logic
        end if        ! end of test on i0 and i1
!
!
      else                            ! EXIT end
        if( i0.eq.1 .and. i1.eq.1 ) then
          diff = ABS( rsp - rs1 )
          if( diff .le. eps ) then  ! step ends on 1st sphere boundary
            ihwin = 2
            hregn = dz(2) / 10.     ! step size for NEXT step
          end if
          go to 800       ! good step
!
        else if( i0.eq.1 .and. i1.eq.2 ) then 
          diff0 = ABS( rsp0 - rs1 )
          diff = ABS( rsp - rs1 )
          if( diff0.le.eps .or. diff.le.eps ) then
! One end of the step ends on the 1st sphere boundary
            ihwin = 2
            hregn = dz(2) / 10.
            if( diff0 .lt. eps ) i0 = 2 ! get right material
            go to 800       ! good step
          end if
          r2 = rs1**2       ! set radius for step adjustment logic
!
        else if( i0.eq.1 .and. i1.eq.3 ) then 
          diff0 = ABS( rsp0 - rs1 )
          if( diff0 .le. eps ) then   ! step starts on 1st sphere boundary
            ihwin = 2
            hregn = dz(2) / 10.
            iflag = 1
            go to 900       ! redo step
          end if
          r2 = rs1**2       ! set radius for step adjustment logic
!
        else if( i0.eq.2 .and. i1.eq.2 ) then 
          diff = ABS( rsp - rs2 )
          if( diff .le. eps ) then  ! step ends on 2nd sphere boundary
            ihwin = 3
            hregn = dz(3) / 10.     ! step size for NEXT step
          end if
          go to 800       ! good step
!
        else if( i0.eq.2 .and. i1.eq.3 ) then 
          diff0 = ABS( rsp0 - rs2 )
          diff = ABS( rsp - rs2 )
          if( diff0.le.eps .or. diff.le.eps ) then
! One end of the step ends on the 2nd sphere boundary
            ihwin = 3
            hregn = dz(3) / 10.
            if( diff0 .lt. eps ) i0 = 3 ! get right material
            go to 800       ! good step
          end if
          r2 = rs2**2       ! set radius for step adjustment logic
!
        else if( i0.eq.3 .and. i1.eq.3 ) then 
          go to 800     ! good step
!
! It's also possible for trajectories to cross subregions backwards
        else if( i0.eq.2 .and. i1.eq.1 ) then 
          diff0 = ABS( rsp0 - rs1 )
          diff = ABS( rsp - rs1 )
          if( diff0.le.eps .or. diff.le.eps ) then
! One end of the step ends on the 1st sphere boundary
! Leave the small stepping appropriate to region 2
            if( diff0 .lt. eps ) i0 = 1 ! but use right material
            go to 800       ! good step
          end if
          r2 = rs1**2       ! set radius for step adjustment logic
!
        else if( i0.eq.3 .and. i1.eq.2 ) then 
          diff0 = ABS( rsp0 - rs2 )
          diff = ABS( rsp - rs2 )
          if( diff0.le.eps .or. diff.le.eps ) then
! One end of the step ends on the 2nd sphere boundary
            ihwin = 2
            hregn = dz(2) / 10.
            if( diff0 .lt. eps ) i0 = 2 ! use right material
            go to 800       ! good step
          end if
          r2 = rs2**2       ! set radius for step adjustment logic
!
        else if( i0.eq.3 .and. i1.eq.1 ) then 
          diff0 = ABS( rsp0 - rs2 )
          if( diff0 .le. eps ) then   ! step starts on 2nd sphere boundary
            ihwin = 2
            hregn = dz(2) / 10.
            iflag = 1
            go to 900       ! redo step
          end if
          r2 = rs2**2       ! set radius for step adjustment logic
!
        end if    ! end of test on i0 and i1
      end if      ! end of test on gparm(1,1)
! ---------------------------------------------------------
! Step crosses significantly over the region boundary
! Readjust the step size to the remaining path length
!  in the incident region
      do i=1,3
      dp(i) = x(i) - x0(i)
      end do
!
      dp2 = dp(1)**2 + dp(2)**2 + dp(3)**2
!
      if( dp2 .eq. 0. ) then
        iflag = -58
        go to 900
      end if
!
      p12 = x0(1)**2 + x0(2)**2 + x0(3)**2
      aqf = dp2
      bqf = 2. * DOT3(x0,dp)
      cqf = p12 - r2
      u = -bqf + gparm(1,1)*SQRT( bqf**2 - 4.*aqf*cqf )
      u = u / (2.*aqf)
      if( u.lt.-eps .or. u.gt.1.+eps ) then
        iflag = -58
        go to 900
      end if
      hregn = u * (xp(3)-xp0(3))    ! fractional step size
      iflag = 1       ! signal to redo this step
      go to 900
! .........................................................
800   continue        ! good step
      mt = mtag(i0)
      if( mt .ne. 'VAC' ) then
        call LOAD_MATERIAL
        in_hemiwin = .true.
      end if
! .........................................................
900   continue
      return
      end
c  *********************************************************
      logical function IN_NIA(hstp)
!
!  Determine if particle step goes thru a Non-Isoceles Absorber
!
!  Ref. J.S. Berg, MUC-NOTE-COOL_THEORY-261, October 2002
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      real*8 p0(3),p1(3),u0(3),u1(3),x1(3),t(3),v(3)
      logical in0,in1
      data rad/0.01745329/     ! pi / 180
!
      dum = hstp
!
      in_nia = .false.
      if( mgt .ne. 'NIA' ) go to 900
!
      in0 = .false.
      in1 = .false.
!    
      zv = sr0 + gparm(1,krad)
      dz0 = gparm(2,krad)
      dz1 = gparm(3,krad)
      th0 = gparm(4,krad) * rad
      ph0 = gparm(5,krad) * rad
      th1 = gparm(6,krad) * rad
      ph1 = gparm(7,krad) * rad
!
      do i=1,2
      p0(i) = 0.d0
      p1(i) = 0.d0
      end do
      p0(3) = zv - dz0
      p1(3) = zv + dz1
!
      sth0 = SIN(th0)
      u0(1) = -sth0 * COS(ph0)
      u0(2) = -sth0 * SIN(ph0)
      u0(3) = -COS(th0)
!
      sth1 = SIN(th1)
      u1(1) = -sth1 * COS(ph1)
      u1(2) = -sth1 * SIN(ph1)
      u1(3) = COS(th1)
!
      call VMAG(pp,pm)
      t(1) = pp(1) / pm
      t(2) = pp(2) / pm
      t(3) = pp(3) / pm
!     
!  Check end point
      do i=1,3
      v(i) = xp(i) - p0(i)
      end do
      val0 = DOT3( v,u0 )
      do i=1,3
      v(i) = xp(i) - p1(i)
      end do
      val1 = DOT3( v,u1 )
      if( val0.lt.0. .and. val1.lt.0. ) in1 = .true.
!
!  Check starting point
      x1(1) = xp(1) - pp(1)/pp(3)*hdid
      x1(2) = xp(2) - pp(2)/pp(3)*hdid
      x1(3) = xp(3) - hdid
      do i=1,3
      v(i) = x1(i) - p0(i)
      end do
      val0 = DOT3( v,u0 )
      do i=1,3
      v(i) = x1(i) - p1(i)
      end do
      val1 = DOT3( v,u1 )
      if( val0.lt.0. .and. val1.lt.0. ) in0 = .true.
!
      if( in0 .and. in1 ) then
        in_nia = .true.
!
!      else if( in0 .and. (.not. in1) ) then
!        in_nia = .true.
!        hstp = DOT3( (p1-x1),t )
!     
!      else if( in1 .and. (.not. in0) ) then
!        in_nia = .true.
!        hstp = DOT3( (xp-p0),t )
      end if
!
900   continue
      return
      end
c  *********************************************************
      logical function IN_PWEDGE(hstp)
!
!  Determine if particle step goes thru an asymmetric, polynomial 
!   wedge oriented at an angle phi from the x-axis
!
!   Describe polynomial shape in upper half-plane of x-z coord system
!     dz(x) = a0+a1*x+a2*x**2+a3*x**3
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      real*8 x(2),p(2)
      logical in0,in1
      data rad/0.01745329/     ! pi / 180
!
      in_pwedge = .false.
      if( mgt .ne. 'PWEDGE' ) go to 900
!
      in0 = .false.
      in1 = .false.
      dz0 = 0.
      dz1 = 0.
!
      xv = gparm(2,krad)
      zv = sr0 + gparm(3,krad)
      phi = gparm(4,krad) * rad
      w = gparm(5,krad)              
      h = gparm(6,krad)
      a0 = gparm(7,krad)
      a1 = gparm(8,krad)
      a2 = gparm(9,krad)
      a3 = gparm(10,krad)
!
!  Get limits of problem box
      xmin = xv - w
      xmax = xv     
      ymin = -h / 2.
      ymax = h / 2.
      dz = a0 + a1*xmin + a2*xmin**2 + a3*xmin**3
      zmin = zv - dz
      zmax = zv + dz
!
!  Rotate particle coordinates into frame where wedge vertex lies
!   along the +ve x-axis
      ca = COS(phi)
      sa = SIN(phi)
      x(1) = xp(1)*ca + xp(2)*sa
           x(2) = -xp(1)*sa + xp(2)*ca
           p(1) = pp(1)*ca + pp(2)*sa
           p(2) = -pp(1)*sa + pp(2)*ca
!     
!  Check end point
      if( x(1).gt.xmin .and. x(1).lt.xmax .and.
     &  x(2).gt.ymin .and. x(2).lt.ymax .and. 
     &  xp(3).gt.zmin .and. xp(3).lt.zmax ) then
!  By construction zsurf is a positive quantity in a c.s. 
!   centered at zv
         zsurf = a0 + a1*x(1) + a2*x(1)**2 + a3*x(1)**3
         dz = xp(3) - zv
         if( dz .ge. 0. ) then
           dz1 = xp(3) - (zv + zsurf)
           if( dz1 .le. 0. ) in1 = .true.
         else     ! xp(3) < zv
           dz1 = xp(3) - (zv - zsurf)
           if( dz1 .ge. 0. ) in1 = .true.
         end if
      end if
!
!  Check starting point
      x1 = x(1) - p(1)/pp(3)*hdid
      y1 = x(2) - p(2)/pp(3)*hdid
      z1 = xp(3) - hdid
!
      if( x1.gt.xmin .and. x1.lt.xmax .and. y1.gt.ymin .and.
     &  y1.lt.ymax .and. z1.gt.zmin .and. z1.lt.zmax ) then
         zsurf = a0 + a1*x1 + a2*x1**2 + a3*x1**3
         dz = z1 - zv
         if( dz .ge. 0. ) then
           dz0 = z1 - (zv + zsurf)
           if( dz0 .le. 0. ) in0 = .true.
         else     ! z1 < zv
           dz0 = z1 - (zv - zsurf)
           if( dz1 .ge. 0. ) in0 = .true.
         end if
      end if
!
      if( in0 .and. in1 ) then
        in_pwedge = .true.
!
      else if( in0 .and. (.not. in1) ) then
        in_pwedge = .true.
        hstp = hstp * ABS(dz0) / ( ABS(dz0)+ABS(dz1) )
!     
      else if( in1 .and. (.not. in0) ) then
        in_pwedge = .true.
        hstp = hstp * ABS(dz1) / ( ABS(dz0)+ABS(dz1) )
      end if
!
900   continue
      return
      end
c  *********************************************************
      logical function IN_RING(hstp)
!
! Determine how much of particle step goes thru an annular ring of material
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      logical in1,in2
!
      in_ring = .false.
      stp = hstp
      if( mgt .ne. 'RING' ) go to 900
!
      in1 = .false.
      in2 = .false.
!
! Set geometrical parameters
      a1 = gparm(1,krad)      ! inner radius of ring
      a2 = gparm(2,krad)      ! outer radius of ring
!
! Check location of end of step
      r2 = SQRT( xp(1)**2 + xp(2)**2 )
      if( r2.ge.a1 .and. r2.le.a2 ) in2 = .true.
!
! Check location of start of step
      xpr = pp(1) / pp(3)
      ypr = pp(2) / pp(3)
      x1 = xp(1) - xpr*hdid
      y1 = xp(2) - ypr*hdid
      r1 = SQRT( x1*x1 + y1*y1 )
      if( r1.ge.a1 .and. r1.le.a2 ) in1 = .true.
!
! ---------------------------------------------------------
! Step handling logic (6 cases of step logic)
      if( in1 .and. in2 ) then   ! all of hstp is in material
        in_ring = .true.
        stp = hstp
        go to 900
      end if
!
      if( .not.in1 .and. .not.in2 ) then  ! none of hstp in the material
        stp = 0.
        go to 900
      end if
!
! Step crossed a radial boundary 
      z1 = xp(3) - hdid
      x2 = xp(1)
      y2 = xp(2)
      z2 = xp(3)
      denom = xpr*xpr + ypr*ypr
      f1 = x1*xpr + y1*ypr
!
      if( in1 .and. r2.lt.a1 ) then   ! case 3
        arg = f1*f1 - denom*(x1*x1+y1*y1-a1*a1)
        if( f1 .ge. 0. ) then
          zb = (-f1 + SQRT(arg)) / denom
        else
          zb = (-f1 - SQRT(arg)) / denom
        end if
        xb = xpr*zb
        yb = ypr*zb
        stp = SQRT( xb*xb + yb*yb + zb*zb )
!
      else if( in1 .and. r2.gt.a2 ) then  ! case 4
        arg = f1*f1 - denom*(x1*x1+y1*y1-a2*a2)
        if( f1 .ge. 0. ) then
          zb = (-f1 + SQRT(arg)) / denom
        else
          zb = (-f1 - SQRT(arg)) / denom
        end if
        xb = xpr*zb
        yb = ypr*zb
        stp = SQRT( xb*xb + yb*yb + zb*zb )
!
      else if( in2 .and. r1.lt.a1 ) then  ! case 5
        arg = f1*f1 - denom*(x1*x1+y1*y1-a1*a1)
        if( f1 .ge. 0. ) then
          zb = (-f1 + SQRT(arg)) / denom
        else
          zb = (-f1 - SQRT(arg)) / denom
        end if
        xb = xpr*zb
        yb = ypr*zb
        stp = hstp - SQRT( xb*xb + yb*yb + zb*zb )
!
      else if( in2 .and. r1.gt.a2 ) then  ! case 6
        arg = f1*f1 - denom*(x1*x1+y1*y1-a2*a2)
        if( f1 .ge. 0. ) then
          zb = (-f1 + SQRT(arg)) / denom
        else
          zb = (-f1 - SQRT(arg)) / denom
        end if
        xb = xpr*zb
        yb = ypr*zb
        stp = hstp - SQRT( xb*xb + yb*yb + zb*zb )
      end if
      in_ring = .true.
! ---------------------------------------------------------
900   continue
!      if( stp.lt.0. .or. stp.gt.hstp ) then
!        write(*,'(a,2i5,2e12.4)') 'RING ',j7rec,jpar,xp(3),stp
!        stop 
!      end if
      hstp = stp
      end
c  *********************************************************
      logical function IN_WEDGE(hstp)
!
!  Determine if particle step goes thru an asymmetric wedge
!   oriented at an angle phi from the x-axis
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      real*8 x(2),p(2)
      logical in1,in2
      data rad/0.01745329d0/     ! pi / 180
!
      in_wedge = .false.
      if( mgt .ne. 'WEDGE' ) go to 900
!
      in1 = .false.
      in2 = .false.
!
      alp = gparm(1,krad) * rad
      xv = gparm(2,krad)
      zv = sr0 + gparm(3,krad)
      phi = gparm(4,krad) * rad
      w = gparm(5,krad)              
      h = gparm(6,krad)
      ta = TAN(alp/2.d0)
      dz = w * ta
!
      xmin = xv - w
      xmax = xv     
      ymin = -h / 2.d0
      ymax = h / 2.d0
      zmin = zv - dz
      zmax = zv + dz
!
!  Rotate particle coordinates into frame where wedge vertex lies
!   along the +ve x-axis
      ca = COS(phi)
      sa = SIN(phi)
      x(1) = xp(1)*ca + xp(2)*sa
      x(2) = -xp(1)*sa + xp(2)*ca
      p(1) = pp(1)*ca + pp(2)*sa
      p(2) = -pp(1)*sa + pp(2)*ca
!     
!  Check end point
      if( x(1).gt.xmin .and. x(1).lt.xmax .and.
     &   x(2).gt.ymin .and. x(2).lt.ymax .and. 
     &   xp(3).gt.zmin .and. xp(3).lt.zmax ) then
        dx = ABS(x(1)-xv)
        dzs = dx * ta         ! dz for surface
        dzp = ABS(xp(3)-zv)   ! dz for particle
        if( dzp .le. dzs ) in2 = .true.
      end if
!
!  Check starting point
      x1 = x(1) - p(1)/pp(3)*hdid
      y1 = x(2) - p(2)/pp(3)*hdid
      z1 = xp(3) - hdid
!
      if( x1.gt.xmin .and. x1.lt.xmax .and. y1.gt.ymin .and.
     &  y1.lt.ymax .and. z1.gt.zmin .and. z1.lt.zmax ) then
        dx = ABS(x1-xv)
        dzs = dx * ta         ! dz for surface
        dzp = ABS(z1-zv)      ! dz for particle
        if( dzp .le. dzs ) in1 = .true.
      end if
!
      if( in1 .and. in2 ) then
        in_wedge = .true.
!
      else if( in1 .and. (.not. in2) ) then
        in_wedge = .true.
        dx1 = ABS(x1-xv)
        dx2 = ABS(xp(1)-xv)
        a = zv - z1 + dx1*ta
        b = xp(3) - dx2*ta
        hstp = hstp * ABS(a) / ( ABS(a)+ABS(b) )
!     
      else if( in2 .and. (.not. in1) ) then
        in_wedge = .true.
        dx1 = ABS(x1-xv)
        dx2 = ABS(xp(1)-xv)
        a = -dx1*ta - z1
        b = xp(3) - zv + dx2*ta
        hstp = hstp * ABS(b) / ( ABS(a)+ABS(b) )
      end if
!
900   continue
      end
c  *********************************************************
      subroutine LANDAU(yran)
c
c  Determine fluctuation in ionization energy loss
c   from Landau distribution
c
c  Derived from GEANT routine  GLANDG
C.
C.    ******************************************************************
C.    *                                                                *
C.    *  Copy of the CERN library routine GENLAN                       *
C.    *  Generation of LANDAU-distributed random numbers by 4-point    *
C.    *  interpolation in the previously-tabulated inverse cumulative  *
C.    *  distribution                                                  *
C.    *                                                                *
C.    ******************************************************************
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      PARAMETER (RANGE1=0.807069,RANGE2=0.994869)
      PARAMETER (X1BOT =0.0     ,X2BOT =0.791240)
      PARAMETER (GP1INV=118.836 ,GP2INV=486.178 )
C
      DIMENSION XCUM1(100) , XCUM2(100)
      DIMENSION RNDM(1)
      save xcum1,xcum2
!
C====>  1ST TABLE OF INVERSE CUMULATIVE LANDAU POINTS. ( 0.0<P<0.833 )
      DATA XCUM1(  1) / -2.57000 /   ,   XCUM1(  2) / -2.15412 /
      DATA XCUM1(  3) / -1.94167 /   ,   XCUM1(  4) / -1.79583 /
      DATA XCUM1(  5) / -1.67975 /   ,   XCUM1(  6) / -1.58090 /
      DATA XCUM1(  7) / -1.49341 /   ,   XCUM1(  8) / -1.41397 /
      DATA XCUM1(  9) / -1.34057 /   ,   XCUM1( 10) / -1.27185 /
      DATA XCUM1( 11) / -1.20686 /   ,   XCUM1( 12) / -1.14490 /
      DATA XCUM1( 13) / -1.08545 /   ,   XCUM1( 14) / -1.02809 /
      DATA XCUM1( 15) /  -.97249 /   ,   XCUM1( 16) /  -.91839 /
      DATA XCUM1( 17) /  -.86557 /   ,   XCUM1( 18) /  -.81386 /
      DATA XCUM1( 19) /  -.76308 /   ,   XCUM1( 20) /  -.71311 /
      DATA XCUM1( 21) /  -.66384 /   ,   XCUM1( 22) /  -.61515 /
      DATA XCUM1( 23) /  -.56696 /   ,   XCUM1( 24) /  -.51919 /
      DATA XCUM1( 25) /  -.47176 /   ,   XCUM1( 26) /  -.42462 /
      DATA XCUM1( 27) /  -.37769 /   ,   XCUM1( 28) /  -.33092 /
      DATA XCUM1( 29) /  -.28426 /   ,   XCUM1( 30) /  -.23767 /
      DATA XCUM1( 31) /  -.19109 /   ,   XCUM1( 32) /  -.14447 /
      DATA XCUM1( 33) /  -.09779 /   ,   XCUM1( 34) /  -.05100 /
      DATA XCUM1( 35) /  -.00405 /   ,   XCUM1( 36) /   .04308 /
      DATA XCUM1( 37) /   .09044 /   ,   XCUM1( 38) /   .13806 /
      DATA XCUM1( 39) /   .18598 /   ,   XCUM1( 40) /   .23423 /
      DATA XCUM1( 41) /   .28285 /   ,   XCUM1( 42) /   .33187 /
      DATA XCUM1( 43) /   .38134 /   ,   XCUM1( 44) /   .43128 /
      DATA XCUM1( 45) /   .48174 /   ,   XCUM1( 46) /   .53275 /
      DATA XCUM1( 47) /   .58435 /   ,   XCUM1( 48) /   .63658 /
      DATA XCUM1( 49) /   .68948 /   ,   XCUM1( 50) /   .74309 /
      DATA XCUM1( 51) /   .79746 /   ,   XCUM1( 52) /   .85263 /
      DATA XCUM1( 53) /   .90865 /   ,   XCUM1( 54) /   .96556 /
      DATA XCUM1( 55) /  1.02342 /   ,   XCUM1( 56) /  1.08228 /
      DATA XCUM1( 57) /  1.14219 /   ,   XCUM1( 58) /  1.20322 /
      DATA XCUM1( 59) /  1.26542 /   ,   XCUM1( 60) /  1.32887 /
      DATA XCUM1( 61) /  1.39362 /   ,   XCUM1( 62) /  1.45976 /
      DATA XCUM1( 63) /  1.52736 /   ,   XCUM1( 64) /  1.59650 /
      DATA XCUM1( 65) /  1.66727 /   ,   XCUM1( 66) /  1.73976 /
      DATA XCUM1( 67) /  1.81407 /   ,   XCUM1( 68) /  1.89032 /
      DATA XCUM1( 69) /  1.96860 /   ,   XCUM1( 70) /  2.04905 /
      DATA XCUM1( 71) /  2.13180 /   ,   XCUM1( 72) /  2.21699 /
      DATA XCUM1( 73) /  2.30477 /   ,   XCUM1( 74) /  2.39531 /
      DATA XCUM1( 75) /  2.48878 /   ,   XCUM1( 76) /  2.58540 /
      DATA XCUM1( 77) /  2.68535 /   ,   XCUM1( 78) /  2.78889 /
      DATA XCUM1( 79) /  2.89626 /   ,   XCUM1( 80) /  3.00775 /
      DATA XCUM1( 81) /  3.12364 /   ,   XCUM1( 82) /  3.24429 /
      DATA XCUM1( 83) /  3.37005 /   ,   XCUM1( 84) /  3.50136 /
      DATA XCUM1( 85) /  3.63866 /   ,   XCUM1( 86) /  3.78246 /
      DATA XCUM1( 87) /  3.93334 /   ,   XCUM1( 88) /  4.09194 /
      DATA XCUM1( 89) /  4.25900 /   ,   XCUM1( 90) /  4.43533 /
      DATA XCUM1( 91) /  4.62186 /   ,   XCUM1( 92) /  4.81960 /
      DATA XCUM1( 93) /  5.02974 /   ,   XCUM1( 94) /  5.25368 /
      DATA XCUM1( 95) /  5.49312 /   ,   XCUM1( 96) /  5.74987 /
      DATA XCUM1( 97) /  6.02605 /   ,   XCUM1( 98) /  6.32428 /
      DATA XCUM1( 99) /  6.64773 /   ,   XCUM1(100) /  7.00000 /
C
C====>  2ND TABLE OF INVERSE CUMULATIVE LANDAU POINTS. ( 0.791<P<0.995 )
      DATA XCUM2(  1) /   5.50000 /   ,   XCUM2(  2) /   5.56120 /
      DATA XCUM2(  3) /   5.62347 /   ,   XCUM2(  4) /   5.68684 /
      DATA XCUM2(  5) /   5.75133 /   ,   XCUM2(  6) /   5.81699 /
      DATA XCUM2(  7) /   5.88383 /   ,   XCUM2(  8) /   5.95191 /
      DATA XCUM2(  9) /   6.02125 /   ,   XCUM2( 10) /   6.09190 /
      DATA XCUM2( 11) /   6.16391 /   ,   XCUM2( 12) /   6.23732 /
      DATA XCUM2( 13) /   6.31219 /   ,   XCUM2( 14) /   6.38855 /
      DATA XCUM2( 15) /   6.46646 /   ,   XCUM2( 16) /   6.54597 /
      DATA XCUM2( 17) /   6.62712 /   ,   XCUM2( 18) /   6.70998 /
      DATA XCUM2( 19) /   6.79460 /   ,   XCUM2( 20) /   6.88103 /
      DATA XCUM2( 21) /   6.96935 /   ,   XCUM2( 22) /   7.05962 /
      DATA XCUM2( 23) /   7.15192 /   ,   XCUM2( 24) /   7.24633 /
      DATA XCUM2( 25) /   7.34294 /   ,   XCUM2( 26) /   7.44182 /
      DATA XCUM2( 27) /   7.54306 /   ,   XCUM2( 28) /   7.64676 /
      DATA XCUM2( 29) /   7.75300 /   ,   XCUM2( 30) /   7.86188 /
      DATA XCUM2( 31) /   7.97351 /   ,   XCUM2( 32) /   8.08800 /
      DATA XCUM2( 33) /   8.20548 /   ,   XCUM2( 34) /   8.32610 /
      DATA XCUM2( 35) /   8.44997 /   ,   XCUM2( 36) /   8.57725 /
      DATA XCUM2( 37) /   8.70808 /   ,   XCUM2( 38) /   8.84262 /
      DATA XCUM2( 39) /   8.98103 /   ,   XCUM2( 40) /   9.12349 /
      DATA XCUM2( 41) /   9.27021 /   ,   XCUM2( 42) /   9.42141 /
      DATA XCUM2( 43) /   9.57730 /   ,   XCUM2( 44) /   9.73812 /
      DATA XCUM2( 45) /   9.90410 /   ,   XCUM2( 46) /  10.07552 /
      DATA XCUM2( 47) /  10.25265 /   ,   XCUM2( 48) /  10.43584 /
      DATA XCUM2( 49) /  10.62540 /   ,   XCUM2( 50) /  10.82169 /
      DATA XCUM2( 51) /  11.02508 /   ,   XCUM2( 52) /  11.23598 /
      DATA XCUM2( 53) /  11.45487 /   ,   XCUM2( 54) /  11.68222 /
      DATA XCUM2( 55) /  11.91855 /   ,   XCUM2( 56) /  12.16441 /
      DATA XCUM2( 57) /  12.42045 /   ,   XCUM2( 58) /  12.68734 /
      DATA XCUM2( 59) /  12.96580 /   ,   XCUM2( 60) /  13.25663 /
      DATA XCUM2( 61) /  13.56075 /   ,   XCUM2( 62) /  13.87912 /
      DATA XCUM2( 63) /  14.21280 /   ,   XCUM2( 64) /  14.56296 /
      DATA XCUM2( 65) /  14.93088 /   ,   XCUM2( 66) /  15.31799 /
      DATA XCUM2( 67) /  15.72593 /   ,   XCUM2( 68) /  16.15650 /
      DATA XCUM2( 69) /  16.61170 /   ,   XCUM2( 70) /  17.09378 /
      DATA XCUM2( 71) /  17.60523 /   ,   XCUM2( 72) /  18.14884 /
      DATA XCUM2( 73) /  18.72792 /   ,   XCUM2( 74) /  19.34622 /
      DATA XCUM2( 75) /  20.00796 /   ,   XCUM2( 76) /  20.71792 /
      DATA XCUM2( 77) /  21.48176 /   ,   XCUM2( 78) /  22.30618 /
      DATA XCUM2( 79) /  23.19862 /   ,   XCUM2( 80) /  24.16809 /
      DATA XCUM2( 81) /  25.22547 /   ,   XCUM2( 82) /  26.38320 /
      DATA XCUM2( 83) /  27.65658 /   ,   XCUM2( 84) /  29.06490 /
      DATA XCUM2( 85) /  30.63048 /   ,   XCUM2( 86) /  32.38211 /
      DATA XCUM2( 87) /  34.35555 /   ,   XCUM2( 88) /  36.59607 /
      DATA XCUM2( 89) /  39.16212 /   ,   XCUM2( 90) /  42.13195 /
      DATA XCUM2( 91) /  45.61204 /   ,   XCUM2( 92) /  49.74758 /
      DATA XCUM2( 93) /  54.74189 /   ,   XCUM2( 94) /  60.89935 /
      DATA XCUM2( 95) /  68.67370 /   ,   XCUM2( 96) /  78.81480 /
      DATA XCUM2( 97) /  92.61047 /   ,   XCUM2( 98) / 112.50807 /
      DATA XCUM2( 99) / 143.78539 /   ,   XCUM2(100) / 200.00000 /
C.
      RNDM(1)=RAN1(irnarg)
      X = RNDM(1)
      if( x .le. 1e-9 ) x = 1e-9      ! P.L. correction for x=0
c
      IF(X .LT. 0.004) THEN
c  Extreme left-hand tail
         YRAN = -SQRT(ABS(LOG(X)))
*
      ELSE IF(X .LE. RANGE1) THEN
c  4-point interpolation in the first cumulative table
         TABPO1 = (X-X1BOT)*GP1INV
         J = TABPO1 + 1
         J = MAX(J,2)
         J = MIN(J,98)
         P = TABPO1 - (J-1)
         A = (P+1.0) * XCUM1(J+2) - (P-2.0) * XCUM1(J-1)
         B = (P-1.0) * XCUM1(J)   - P       * XCUM1(J+1)
         YRAN = A*P*(P-1.0)*0.16666667+B*(P+1.0)*(P-2.0)*0.5
*
      ELSE IF(X .LE. RANGE2) THEN
c  4-point interpolation in the second cumulative table
         TABPO2 = (X-X2BOT)*GP2INV
         J = TABPO2 + 1
         J = MAX(J,2)
         J = MIN(J,98)
         P = TABPO2 - (J-1)
         A = (P+1.0) * XCUM2(J+2) - (P-2.0) * XCUM2(J-1)
         B = (P-1.0) * XCUM2(J)   - P       * XCUM2(J+1)
         YRAN = A*P*(P-1.0)*0.16666667+B*(P+1.0)*(P-2.0)*0.5
*
      ELSE
c  1/x**2 sampling for extreme Landau tail
              RNDM(1)=RAN1(irnarg)
              YRAN=200./RNDM(1)
      ENDIF
c
      return
      END

c  *********************************************************
      FUNCTION LANDAUR(X)
C.
C.    ******************************************************************
C.    *                                                                *
C.    *  Copy of the CERN library routine RANLAN for restricted        *
C.    *  Landau distribution.                                          *
C.    *                                                                *
!     *  Modified from GEANT routine  GLANDR                           *
C.    ******************************************************************
C.
      implicit double precision (a-h,o-z)
      double precision landaur
      DIMENSION F(-1:1002)
      DATA (F(I),I= -1,100)
     A/     7*0.,                                        -2.244733,
     B -2.204365,-2.168163,-2.135219,-2.104898,-2.076740,-2.050397,
     C -2.025605,-2.002150,-1.979866,-1.958612,-1.938275,-1.918760,
     D -1.899984,-1.881879,-1.864385,-1.847451,-1.831030,-1.815083,
     E -1.799574,-1.784473,-1.769751,-1.755383,-1.741346,-1.727620,
     F -1.714187,-1.701029,-1.688130,-1.675477,-1.663057,-1.650858,
     G -1.638868,-1.627078,-1.615477,-1.604058,-1.592811,-1.581729,
     H -1.570806,-1.560034,-1.549407,-1.538919,-1.528565,-1.518339,
     I -1.508237,-1.498254,-1.488386,-1.478628,-1.468976,-1.459428,
     J -1.449979,-1.440626,-1.431365,-1.422195,-1.413111,-1.404112,
     K -1.395194,-1.386356,-1.377594,-1.368906,-1.360291,-1.351746,
     L -1.343269,-1.334859,-1.326512,-1.318229,-1.310006,-1.301843,
     M -1.293737,-1.285688,-1.277693,-1.269752,-1.261863,-1.254024,
     N -1.246235,-1.238494,-1.230800,-1.223153,-1.215550,-1.207990,
     O -1.200474,-1.192999,-1.185566,-1.178172,-1.170817,-1.163500,
     P -1.156220,-1.148977,-1.141770,-1.134598,-1.127459,-1.120354,
     Q -1.113282,-1.106242,-1.099233,-1.092255/
C
      DATA (F(I),I=101,200)
     A/-1.085306,-1.078388,-1.071498,-1.064636,-1.057802,-1.050996,
     B -1.044215,-1.037461,-1.030733,-1.024029,-1.017350,-1.010695,
     C -1.004064, -.997456, -.990871, -.984308, -.977767, -.971247,
     D  -.964749, -.958271, -.951813, -.945375, -.938957, -.932558,
     E  -.926178, -.919816, -.913472, -.907146, -.900838, -.894547,
     F  -.888272, -.882014, -.875773, -.869547, -.863337, -.857142,
     G  -.850963, -.844798, -.838648, -.832512, -.826390, -.820282,
     H  -.814187, -.808106, -.802038, -.795982, -.789940, -.783909,
     I  -.777891, -.771884, -.765889, -.759906, -.753934, -.747973,
     J  -.742023, -.736084, -.730155, -.724237, -.718328, -.712429,
     K  -.706541, -.700661, -.694791, -.688931, -.683079, -.677236,
     L  -.671402, -.665576, -.659759, -.653950, -.648149, -.642356,
     M  -.636570, -.630793, -.625022, -.619259, -.613503, -.607754,
     N  -.602012, -.596276, -.590548, -.584825, -.579109, -.573399,
     O  -.567695, -.561997, -.556305, -.550618, -.544937, -.539262,
     P  -.533592, -.527926, -.522266, -.516611, -.510961, -.505315,
     Q  -.499674, -.494037, -.488405, -.482777/
C
      DATA (F(I),I=201,300)
     A/ -.477153, -.471533, -.465917, -.460305, -.454697, -.449092,
     B  -.443491, -.437893, -.432299, -.426707, -.421119, -.415534,
     C  -.409951, -.404372, -.398795, -.393221, -.387649, -.382080,
     D  -.376513, -.370949, -.365387, -.359826, -.354268, -.348712,
     E  -.343157, -.337604, -.332053, -.326503, -.320955, -.315408,
     F  -.309863, -.304318, -.298775, -.293233, -.287692, -.282152,
     G  -.276613, -.271074, -.265536, -.259999, -.254462, -.248926,
     H  -.243389, -.237854, -.232318, -.226783, -.221247, -.215712,
     I  -.210176, -.204641, -.199105, -.193568, -.188032, -.182495,
     J  -.176957, -.171419, -.165880, -.160341, -.154800, -.149259,
     K  -.143717, -.138173, -.132629, -.127083, -.121537, -.115989,
     L  -.110439, -.104889, -.099336, -.093782, -.088227, -.082670,
     M  -.077111, -.071550, -.065987, -.060423, -.054856, -.049288,
     N  -.043717, -.038144, -.032569, -.026991, -.021411, -.015828,
     O  -.010243, -.004656,  .000934,  .006527,  .012123,  .017722,
     P   .023323,  .028928,  .034535,  .040146,  .045759,  .051376,
     Q   .056997,  .062620,  .068247,  .073877/
C
      DATA (F(I),I=301,400)
     A/  .079511,  .085149,  .090790,  .096435,  .102083,  .107736,
     B   .113392,  .119052,  .124716,  .130385,  .136057,  .141734,
     C   .147414,  .153100,  .158789,  .164483,  .170181,  .175884,
     D   .181592,  .187304,  .193021,  .198743,  .204469,  .210201,
     E   .215937,  .221678,  .227425,  .233177,  .238933,  .244696,
     F   .250463,  .256236,  .262014,  .267798,  .273587,  .279382,
     G   .285183,  .290989,  .296801,  .302619,  .308443,  .314273,
     H   .320109,  .325951,  .331799,  .337654,  .343515,  .349382,
     I   .355255,  .361135,  .367022,  .372915,  .378815,  .384721,
     J   .390634,  .396554,  .402481,  .408415,  .414356,  .420304,
     K   .426260,  .432222,  .438192,  .444169,  .450153,  .456145,
     L   .462144,  .468151,  .474166,  .480188,  .486218,  .492256,
     M   .498302,  .504356,  .510418,  .516488,  .522566,  .528653,
     N   .534747,  .540850,  .546962,  .553082,  .559210,  .565347,
     O   .571493,  .577648,  .583811,  .589983,  .596164,  .602355,
     P   .608554,  .614762,  .620980,  .627207,  .633444,  .639689,
     Q   .645945,  .652210,  .658484,  .664768/
C
      DATA (F(I),I=401,500)
     A/  .671062,  .677366,  .683680,  .690004,  .696338,  .702682,
     B   .709036,  .715400,  .721775,  .728160,  .734556,  .740963,
     C   .747379,  .753807,  .760246,  .766695,  .773155,  .779627,
     D   .786109,  .792603,  .799107,  .805624,  .812151,  .818690,
     E   .825241,  .831803,  .838377,  .844962,  .851560,  .858170,
     F   .864791,  .871425,  .878071,  .884729,  .891399,  .898082,
     G   .904778,  .911486,  .918206,  .924940,  .931686,  .938446,
     H   .945218,  .952003,  .958802,  .965614,  .972439,  .979278,
     I   .986130,  .992996,  .999875, 1.006769, 1.013676, 1.020597,
     J  1.027533, 1.034482, 1.041446, 1.048424, 1.055417, 1.062424,
     K  1.069446, 1.076482, 1.083534, 1.090600, 1.097681, 1.104778,
     L  1.111889, 1.119016, 1.126159, 1.133316, 1.140490, 1.147679,
     M  1.154884, 1.162105, 1.169342, 1.176595, 1.183864, 1.191149,
     N  1.198451, 1.205770, 1.213105, 1.220457, 1.227826, 1.235211,
     O  1.242614, 1.250034, 1.257471, 1.264926, 1.272398, 1.279888,
     P  1.287395, 1.294921, 1.302464, 1.310026, 1.317605, 1.325203,
     Q  1.332819, 1.340454, 1.348108, 1.355780/
C
      DATA (F(I),I=501,600)
     A/ 1.363472, 1.371182, 1.378912, 1.386660, 1.394429, 1.402216,
     B  1.410024, 1.417851, 1.425698, 1.433565, 1.441453, 1.449360,
     C  1.457288, 1.465237, 1.473206, 1.481196, 1.489208, 1.497240,
     D  1.505293, 1.513368, 1.521465, 1.529583, 1.537723, 1.545885,
     E  1.554068, 1.562275, 1.570503, 1.578754, 1.587028, 1.595325,
     F  1.603644, 1.611987, 1.620353, 1.628743, 1.637156, 1.645593,
     G  1.654053, 1.662538, 1.671047, 1.679581, 1.688139, 1.696721,
     H  1.705329, 1.713961, 1.722619, 1.731303, 1.740011, 1.748746,
     I  1.757506, 1.766293, 1.775106, 1.783945, 1.792810, 1.801703,
     J  1.810623, 1.819569, 1.828543, 1.837545, 1.846574, 1.855631,
     K  1.864717, 1.873830, 1.882972, 1.892143, 1.901343, 1.910572,
     L  1.919830, 1.929117, 1.938434, 1.947781, 1.957158, 1.966566,
     M  1.976004, 1.985473, 1.994972, 2.004503, 2.014065, 2.023659,
     N  2.033285, 2.042943, 2.052633, 2.062355, 2.072110, 2.081899,
     O  2.091720, 2.101575, 2.111464, 2.121386, 2.131343, 2.141334,
     P  2.151360, 2.161421, 2.171517, 2.181648, 2.191815, 2.202018,
     Q  2.212257, 2.222533, 2.232845, 2.243195/
C
      DATA (F(I),I=601,700)
     A/ 2.253582, 2.264006, 2.274468, 2.284968, 2.295507, 2.306084,
     B  2.316701, 2.327356, 2.338051, 2.348786, 2.359562, 2.370377,
     C  2.381234, 2.392131, 2.403070, 2.414051, 2.425073, 2.436138,
     D  2.447246, 2.458397, 2.469591, 2.480828, 2.492110, 2.503436,
     E  2.514807, 2.526222, 2.537684, 2.549190, 2.560743, 2.572343,
     F  2.583989, 2.595682, 2.607423, 2.619212, 2.631050, 2.642936,
     G  2.654871, 2.666855, 2.678890, 2.690975, 2.703110, 2.715297,
     H  2.727535, 2.739825, 2.752168, 2.764563, 2.777012, 2.789514,
     I  2.802070, 2.814681, 2.827347, 2.840069, 2.852846, 2.865680,
     J  2.878570, 2.891518, 2.904524, 2.917588, 2.930712, 2.943894,
     K  2.957136, 2.970439, 2.983802, 2.997227, 3.010714, 3.024263,
     L  3.037875, 3.051551, 3.065290, 3.079095, 3.092965, 3.106900,
     M  3.120902, 3.134971, 3.149107, 3.163312, 3.177585, 3.191928,
     N  3.206340, 3.220824, 3.235378, 3.250005, 3.264704, 3.279477,
     O  3.294323, 3.309244, 3.324240, 3.339312, 3.354461, 3.369687,
     P  3.384992, 3.400375, 3.415838, 3.431381, 3.447005, 3.462711,
     Q  3.478500, 3.494372, 3.510328, 3.526370/
C
      DATA (F(I),I=701,800)
     A/ 3.542497, 3.558711, 3.575012, 3.591402, 3.607881, 3.624450,
     B  3.641111, 3.657863, 3.674708, 3.691646, 3.708680, 3.725809,
     C  3.743034, 3.760357, 3.777779, 3.795300, 3.812921, 3.830645,
     D  3.848470, 3.866400, 3.884434, 3.902574, 3.920821, 3.939176,
     E  3.957640, 3.976215, 3.994901, 4.013699, 4.032612, 4.051639,
     F  4.070783, 4.090045, 4.109425, 4.128925, 4.148547, 4.168292,
     G  4.188160, 4.208154, 4.228275, 4.248524, 4.268903, 4.289413,
     H  4.310056, 4.330832, 4.351745, 4.372794, 4.393982, 4.415310,
     I  4.436781, 4.458395, 4.480154, 4.502060, 4.524114, 4.546319,
     J  4.568676, 4.591187, 4.613854, 4.636678, 4.659662, 4.682807,
     K  4.706116, 4.729590, 4.753231, 4.777041, 4.801024, 4.825179,
     L  4.849511, 4.874020, 4.898710, 4.923582, 4.948639, 4.973883,
     M  4.999316, 5.024942, 5.050761, 5.076778, 5.102993, 5.129411,
     N  5.156034, 5.182864, 5.209903, 5.237156, 5.264625, 5.292312,
     O  5.320220, 5.348354, 5.376714, 5.405306, 5.434131, 5.463193,
     P  5.492496, 5.522042, 5.551836, 5.581880, 5.612178, 5.642734,
     Q  5.673552, 5.704634, 5.735986, 5.767610/
C
      DATA (F(I),I=801,900)
     A/ 5.799512, 5.831694, 5.864161, 5.896918, 5.929968, 5.963316,
     B  5.996967, 6.030925, 6.065194, 6.099780, 6.134687, 6.169921,
     C  6.205486, 6.241387, 6.277630, 6.314220, 6.351163, 6.388465,
     D  6.426130, 6.464166, 6.502578, 6.541371, 6.580553, 6.620130,
     E  6.660109, 6.700495, 6.741297, 6.782520, 6.824173, 6.866262,
     F  6.908795, 6.951780, 6.995225, 7.039137, 7.083525, 7.128398,
     G  7.173764, 7.219632, 7.266011, 7.312910, 7.360339, 7.408308,
     H  7.456827, 7.505905, 7.555554, 7.605785, 7.656608, 7.708035,
     I  7.760077, 7.812747, 7.866057, 7.920019, 7.974647, 8.029953,
     J  8.085952, 8.142657, 8.200083, 8.258245, 8.317158, 8.376837,
     K  8.437300, 8.498562, 8.560641, 8.623554, 8.687319, 8.751955,
     L  8.817481, 8.883916, 8.951282, 9.019600, 9.088889, 9.159174,
     M  9.230477, 9.302822, 9.376233, 9.450735, 9.526355, 9.603118,
     N  9.681054, 9.760191, 9.840558, 9.922186,10.005107,10.089353,
     O 10.174959,10.261958,10.350389,10.440287,10.531693,10.624646,
     P 10.719188,10.815362,10.913214,11.012789,11.114137,11.217307,
     Q 11.322352,11.429325,11.538283,11.649285/
C
      DATA (F(I),I=901,1002)
     A/11.762390,11.877664,11.995170,12.114979,12.237161,12.361791,
     B 12.488946,12.618708,12.751161,12.886394,13.024498,13.165570,
     C 13.309711,13.457026,13.607625,13.761625,13.919145,14.080314,
     D 14.245263,14.414134,14.587072,14.764233,14.945778,15.131877,
     E 15.322712,15.518470,15.719353,15.925570,16.137345,16.354912,
     F 16.578520,16.808433,17.044929,17.288305,17.538873,17.796967,
     G 18.062943,18.337176,18.620068,18.912049,19.213574,19.525133,
     H 19.847249,20.180480,20.525429,20.882738,21.253102,21.637266,
     I 22.036036,22.450278,22.880933,23.329017,23.795634,24.281981,
     J 24.789364,25.319207,25.873062,26.452634,27.059789,27.696581,
     K 28.365274,29.068370,29.808638,30.589157,31.413354,32.285060,
     L 33.208568,34.188705,35.230920,36.341388,37.527131,38.796172,
     M 40.157721,41.622399,43.202525,44.912465,46.769077,48.792279,
     N 51.005773,53.437996,56.123356,59.103894,              20*0./
C.
C.    ------------------------------------------------------------------
C.
      U=1000.0*X
      I=U
      U=U-FLOAT(I)
      LANDauR=F(I)+U*(F(I+1)-F(I))
      IF(I .GE. 70 .AND. I .LE. 800) RETURN
      LANDauR=LANDauR-0.25*U*(1.0-U)*(F(I+2)-F(I+1)-F(I)+F(I-1))
      IF(I .GE.  7 .AND. I .LE. 980) RETURN
      IF(I. LT. 7) THEN
         IF(X.LT.1.E-10)THEN
            V=-30.
         ELSE
            V=LOG(X)
         ENDIF
         U=1.0/V
         LANDauR=((0.99858950+(3.45213058E1+1.70854528E1*U)*U)/
     1           (1.0       +(3.41760202E1+4.01244582  *U)*U))*
     2           (-LOG(-0.91893853-V)-1.0)
      ELSE
         U=1.0-X
         V=U**2
         IF(X .LE. 0.999) THEN
            LANDauR=(1.00060006+2.63991156E2*U+4.37320068E3*V)/
     1             ((1.0       +2.57368075E2*U+3.41448018E3*V)*U)
         ELSE
            LANDauR=(1.00001538+6.07514119E3*U+7.34266409E5*V)/
     1             ((1.0       +6.06511919E3*U+6.94021044E5*V)*U)
         ENDIF
      ENDIF
      END

! *********************************************************
	SUBROUTINE MCSDIS(PL,TMY2,P0,A,Z,I,XIA2,FTOT)
!
! Auxiliary routine for SAMCS
!
	IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)  	
!	REAL RNDM
      common/rans/iseed
	PARAMETER (PI=3.141592653589793227D+00)
!	CALL RNPSSN(SNGL(PL),NS,IERR)
      ns = POIDEV(pl,iseed)
	FTOT=0.D0
	IF(NS.GT.0) THEN
	   DO NN=1,NS
! 5            GT1=DBLE(RNDM(-1.))
 5            GT1 = RAN1(iseed)
	      T12=(TMY2+XIA2)/GT1-XIA2
              QMIN=P0*SQRT(TMY2)
              Q=P0*SQRT(T12)
	      F2MIN=FORMN(QMIN,Z,A)*FORMP(QMIN,I)
              F2Q=FORMN(Q,Z,A)*FORMP(Q,I) 
!	      GT2=DBLE(RNDM(-1.))
	      GT2 = RAN1(iseed)
	      TEST=F2Q/GT2/F2MIN
	      IF(TEST.LT.1.D0) GO TO 5
	      TS=SQRT(T12)
!	      G1=2.D0*PI*DBLE(RNDM(-1.))
	      G1 = 2.D0*PI*RAN1(iseed)
	      FTOT=FTOT+T12-2.D0*SQRT(FTOT)*TS*SIN(G1)
	END DO
	END IF
      RETURN
      END
c  *********************************************************
      subroutine MOLIERE(omega,chic2,th)
c
c  Determine scattering angle from Moliere distribution
!
!   Ref. H. Bethe, Phys. Rev. 89:1256-1266, 1953
!   Derived from GEANT program  GMOLIE
c
c     omega : mean number of scatters
c     chic2 : characteristic scattering angle (squared)
C.    *                                                                *
C.    * THRED(NA)=reduced angles of Moliere theory                     *
C.    *                                                                *
C.    * F0I(NA),F1I(NA),F2I(NA)= integrale of Moliere functions        *
C.    *                                                                *
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      PARAMETER (ENEPER = 2.7182818d0) !,pi=3.141592654d0)
      DIMENSION TINT(40),ARG(4),VAL(4),THRED(40),F0I(40),F1I(40),F2I(40)
      DATA THRED/
     +   0.00, 0.10, 0.20, 0.30
     +,  0.40, 0.50, 0.60, 0.70
     +,  0.80, 0.90, 1.00, 1.10
     +,  1.20, 1.30, 1.40, 1.50
     +,  1.60, 1.70, 1.80, 1.90
     +,  2.00, 2.20, 2.40, 2.60
     +,  2.80, 3.00, 3.20, 3.40
     +,  3.60, 3.80, 4.00, 5.00
     +,  6.00, 7.00, 8.00, 9.00
     +, 10.00,11.00,12.00,13.00/
      DATA F0I/
     +  0.000000E+00 ,0.995016E-02 ,0.392106E-01 ,0.860688E-01
     + ,0.147856E+00 ,0.221199E+00 ,0.302324E+00 ,0.387374E+00
     + ,0.472708E+00 ,0.555142E+00 ,0.632121E+00 ,0.701803E+00
     + ,0.763072E+00 ,0.815480E+00 ,0.859142E+00 ,0.894601E+00
     + ,0.922695E+00 ,0.944424E+00 ,0.960836E+00 ,0.972948E+00
     + ,0.981684E+00 ,0.992093E+00 ,0.996849E+00 ,0.998841E+00
     + ,0.999606E+00 ,0.999877E+00 ,0.999964E+00 ,0.999990E+00
     + ,0.999998E+00 ,0.999999E+00 ,0.100000E+01 ,0.100000E+01
     + ,0.100000E+01 ,0.100000E+01 ,0.100000E+01 ,0.100000E+01
     + ,1.,1.,1.,1./
      DATA F1I/
     +  0.000000E+00,0.414985E-02,0.154894E-01,0.310312E-01
     + ,0.464438E-01,0.569008E-01,0.580763E-01,0.468264E-01
     + ,0.217924E-01,-0.163419E-01,-0.651205E-01,-0.120503E+00
     + ,-0.178272E+00,-0.233580E+00,-0.282442E+00,-0.321901E+00
     + ,-0.350115E+00,-0.366534E+00,-0.371831E+00,-0.367378E+00
     + ,-0.354994E+00,-0.314803E+00,-0.266539E+00,-0.220551E+00
     + ,-0.181546E+00,-0.150427E+00,-0.126404E+00,-0.107830E+00
     + ,-0.933106E-01,-0.817375E-01,-0.723389E-01,-0.436650E-01
     + ,-0.294700E-01,-0.212940E-01,-0.161406E-01,-0.126604E-01
     + ,-0.102042E-01,-0.840465E-02,-0.704261E-02,-0.598886E-02/
      DATA F2I/
     +  0.0,0.121500E-01,0.454999E-01,0.913000E-01
     + ,0.137300E+00,0.171400E+00,0.183900E+00,0.170300E+00
     + ,0.132200E+00,0.763000E-01,0.126500E-01,-0.473500E-01
     + ,-0.936000E-01,-0.119750E+00,-0.123450E+00,-0.106300E+00
     + ,-0.732800E-01,-0.312400E-01,0.128450E-01,0.528800E-01
     + ,0.844100E-01,0.114710E+00,0.106200E+00,0.765830E-01
     + ,0.435800E-01,0.173950E-01,0.695001E-03,-0.809500E-02
     + ,-0.117355E-01,-0.125449E-01,-0.120280E-01,-0.686530E-02
     + ,-0.385275E-02,-0.231115E-02,-0.147056E-02,-0.982480E-03
     + ,-0.682440E-03,-0.489715E-03,-0.361190E-03,-0.272582E-03/
*
*     ------------------------------------------------------------------
*
*
* *** Compute Theta angle from Moliere distribution
      IF(OMEGA.LE.ENEPER)GO TO 90
      CHIC  = sqrt(chic2)
      CNST=LOG(OMEGA)
      B=5.
c
      DO 10 L=1,10
         IF(ABS(B).LT.1.E-10)THEN
            B=1.E-10
         ENDIF
         DB=-(B-LOG(ABS(B))-CNST)/(1.-1./B)
         B=B+DB
         IF(ABS(DB).LE.0.0001)GO TO 20
   10 CONTINUE
c
      GO TO 90
c
   20 CONTINUE
      IF(B.LE.0.)GO TO 90
      BINV = 1./B
      TINT(1) = 0.
c
      DO 30 JA=2,4
         TINT(JA)=F0I(JA)+(F1I(JA)+F2I(JA)*BINV)*BINV
   30 CONTINUE
c
      NMAX = 4
   40 CONTINUE
      XINT = RAN1(irnarg)
c
      DO 50 NA=3,40
         IF(NA.GT.NMAX) THEN
            TINT(NA)=F0I(NA)+(F1I(NA)+F2I(NA)*BINV)*BINV
            NMAX=NA
         ENDIF
         IF(XINT.LE.TINT(NA-1)) GO TO 60
   50 CONTINUE
c
      IF(XINT.LE.TINT(40)) THEN
         NA=40
         GOTO 60
      ELSE
         TMP=1.-(1.-B*(1.-XINT))**5
         IF(TMP.LE.0.)GO TO 40
         THRI=5./TMP
         GO TO 80
      ENDIF
c
   60 CONTINUE
      NA = MAX(NA-1,3)
      NA3 = NA-3
c
      DO 70 M=1,4
         NA3M=NA3+M
         ARG(M)=TINT(NA3M)
         VAL(M)=THRED(NA3M)**2
   70 CONTINUE
c
      call POLINT(arg,val,4,xint,thri,err)    ! Num Rec
c
   80 CONTINUE
      TH=CHIC*SQRT(ABS(B*THRI))
      IF(TH.GT.PI)GO TO 40
      SINTH=SIN(TH)
      TEST=TH*(RAN1(irnarg))**2
      IF(TEST.GT.SINTH)GO TO 40
      GOTO 900
!
   90 CONTINUE
      iflag = -85
c
900   continue
      return
      END
c  *********************************************************
      subroutine MUON_DECAY(pcm)
!
! Find muon decay particle 4-vectors in muon rest frame
!  Ignore muon polarization
!
! 4-vectors stored in array PCM in order 
!   1: e    2: neutrino     3:anti-neutrino
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
      real*8 pcm(4,3)
!
! Sample energy from appropriate  V-A   x  distribution
!  Ignore correlations among the 3 decay products
      if( iptyp .gt. 0 ) then       ! mu+ decay
! e+ distribution
100   continue
        f1 = RANV(0.d0,1.333d0)
        x = (0.75d0*f1)**0.6667d0
        yhi = 2d0*SQRT(x)            ! comparison function
        y = RANV(0.d0,yhi)
        p = 2d0*x*x*(3d0-2d0*x)        ! probability distribution
        if( y .gt. p ) go to 100          ! reject it
        pcm(4,1) = mass(2) /2d0 * x
        arg = pcm(4,1)**2 - mass(1)**2
        if( arg .gt. 0d0 ) then
          pm = SQRT( arg )
        else
          iflag = -30
          go to 900
        end if
        cth = RANV(-1.d0,1.d0)
        sth = SQRT(1.d0-cth*cth)
        phi = RANV(0.d0,twopi)
        pcm(1,1) = pm * sth * COS(phi)
        pcm(2,1) = pm * sth * SIN(phi)
        pcm(3,1) = pm * cth
!
! anti nu_mu distribution
120   continue
        f1 = RANV(0.d0,1.333d0)
        x = (0.75d0*f1)**0.6667d0
        yhi = 2d0*SQRT(x)            ! comparison function
        y = RANV(0.d0,yhi)
        p = 2d0*x*x*(3d0-2d0*x)
        if( y .gt. p ) go to 120          ! reject it
        pcm(4,3) = mass(2) /2d0 * x
        pm = pcm(4,3)
        cth = RANV(-1.d0,1.d0)
        sth = SQRT(1.d0-cth*cth)
        phi = RANV(0.d0,twopi)
        pcm(1,3) = pm * sth * COS(phi)
        pcm(2,3) = pm * sth * SIN(phi)
        pcm(3,3) = pm * cth
!
! nu_e distribution
140   continue
        g1 = RANV(0.d0,1.448d0)
        x = (g1/1.45d0)**0.6667d0
        yhi = 2.18d0*SQRT(x)            ! comparison function
        y = RANV(0.d0,yhi)
        p = 12d0*x*x*(1d0-x)
        if( y .gt. p ) go to 140          ! reject it
        pcm(4,2) = mass(2) /2d0 * x
        pm = pcm(4,2)
        cth = RANV(-1.d0,1.d0)
        sth = SQRT(1.d0-cth*cth)
        phi = RANV(0.d0,twopi)
        pcm(1,2) = pm * sth * COS(phi)
        pcm(2,2) = pm * sth * SIN(phi)
        pcm(3,2) = pm * cth
!
      else                          ! mu- decay
! e- distribution
200   continue
        f1 = RANV(0.d0,1.333d0)
        x = (0.75d0*f1)**0.6667d0
        yhi = 2d0*SQRT(x)            ! comparison function
        y = RANV(0.d0,yhi)
        p = 2d0*x*x*(3d0-2d0*x)
        if( y .gt. p ) go to 200          ! reject it
        pcm(4,1) = mass(2) /2d0 * x
        arg = pcm(4,1)**2 - mass(1)**2
        if( arg .gt. 0d0 ) then
          pm = SQRT( arg )
        else
          iflag = -30
          go to 900
        end if
        cth = RANV(-1.d0,1.d0)
        sth = SQRT(1.d0-cth*cth)
        phi = RANV(0.d0,twopi)
        pcm(1,1) = pm * sth * COS(phi)
        pcm(2,1) = pm * sth * SIN(phi)
        pcm(3,1) = pm * cth
!
! nu_mu distribution
220   continue
        f1 = RANV(0.d0,1.333d0)
        x = (0.75d0*f1)**0.6667d0
        yhi = 2d0*SQRT(x)            ! comparison function
        y = RANV(0.d0,yhi)
        p = 2d0*x*x*(3d0-2d0*x)
        if( y .gt. p ) go to 220          ! reject it
        pcm(4,2) = mass(2) /2d0 * x
        pm = pcm(4,2)
        cth = RANV(-1.d0,1.d0)
        sth = SQRT(1.d0-cth*cth)
        phi = RANV(0.d0,twopi)
        pcm(1,2) = pm * sth * COS(phi)
        pcm(2,2) = pm * sth * SIN(phi)
        pcm(3,2) = pm * cth
!
! anti nu_e distribution
240   continue
        g1 = RANV(0.d0,1.448d0)
        x = (g1/1.45d0)**0.6667d0
        yhi = 2.18d0*SQRT(x)            ! comparison function
        y = RANV(0.d0,yhi)
        p = 12d0*x*x*(1d0-x)
        if( y .gt. p ) go to 240          ! reject it
        pcm(4,3) = mass(2) /2d0 * x
        pm = pcm(4,3)
        cth = RANV(-1.d0,1.d0)
        sth = SQRT(1.d0-cth*cth)
        phi = RANV(0.d0,twopi)
        pcm(1,3) = pm * sth * COS(phi)
        pcm(2,3) = pm * sth * SIN(phi)
        pcm(3,3) = pm * cth
      end if
!
900   continue
      end
c  *********************************************************
      subroutine OUTPUT_NU(ipt,ppn,wt)
!
!  Write out neutrino production data
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
      real*8 ppn(3)
!
      pm = SQRT( ppn(1)**2 + ppn(2)**2 + ppn(3)**2 )
      thet = ACOS( ppn(3)/pm )
!
      if( thet.ge.nuthmin .and. thet.le.nuthmax ) then
        write(neutrino,100)ievt,ipnum,ipt,ipflg,jsrg,tp,xp,ppn,wt
100     format(i7,2i4,i5,i6,e16.8,8e13.5)
      end if
!
      return
      end
! *********************************************************
      subroutine PANG(p,cthet,sthet,cphi,sphi)
!
! Get sine & cosine of spherical angles theta & phi for momentum p
!
      implicit double precision (a-h,o-z)
      real*8 p(3)
!
      pm = SQRT( p(1)*p(1) + p(2)*p(2) + p(3)*p(3) )
      cthet = p(3) / pm
      thet = ACOS( cthet )
      sthet = SIN( thet )
! Don't use  SQRT  or you dont get negative SIN values
!      sthet = SQRT( 1. - cthet*cthet )
      pxy = SQRT( p(1)*p(1) + p(2)*p(2) )
      if( pxy .gt. 0d0 ) then
        phi = ATAN2( p(2),p(1) )
        cphi = COS( phi )
        sphi = SIN( phi )
      else
        cphi = 1d0
        sphi = 0d0
      end if
!
      return
      end
! ********************************************************* 
	REAL*8 FUNCTION POT(Z)
!
! Auxiliary routine for SAMCS
!
C     Ionization potential (MeV)
	IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)  
C     FOR H_2 GAS I=19.2 EV, FOR H_2 LIQUID I=21.8 EV	
	DIMENSION POTELM(100)
	DATA POTELM/
c     $     21.8000d-06, 41.8000d-06, 40.0000d-06, 63.70d-06,  76.0d-06, 
     $     19.2000d-06, 41.8000d-06, 40.0000d-06, 63.70d-06,  76.0d-06, 
     $     78.0000d-06, 82.0000d-06, 95.0000d-06, 115.0d-06, 137.0d-06,
     $     149.000d-06, 156.000d-06, 166.000d-06, 173.0d-06, 173.0d-06, 
     $     180.000d-06, 174.000d-06, 188.000d-06, 190.0d-06, 191.0d-06,
     $     216.000d-06, 233.000d-06, 245.000d-06, 257.0d-06, 272.0d-06, 
     $     286.000d-06, 297.000d-06, 311.000d-06, 322.0d-06, 330.0d-06,
     $     334.000d-06, 350.000d-06, 347.000d-06, 348.0d-06, 357.0d-06, 
     $     352.000d-06, 363.000d-06, 366.000d-06, 379.0d-06, 393.0d-06,
     $     417.000d-06, 424.000d-06, 428.000d-06, 441.0d-06, 449.0d-06, 
     $     470.000d-06, 470.000d-06, 469.000d-06, 488.0d-06, 488.0d-06,
     $     487.000d-06, 485.000d-06, 491.000d-06, 482.0d-06, 488.0d-06, 
     $     491.000d-06, 501.000d-06, 523.000d-06, 535.0d-06, 546.0d-06,
     $     560.000d-06, 574.000d-06, 580.000d-06, 591.0d-06, 614.0d-06, 
     $     628.000d-06, 650.000d-06, 658.000d-06, 674.0d-06, 684.0d-06,
     $     694.000d-06, 705.000d-06, 718.000d-06, 727.0d-06, 736.0d-06, 
     $     746.000d-06, 757.000d-06, 790.000d-06, 790.0d-06, 800.0d-06,
     $     810.000d-06, 823.000d-06, 823.000d-06, 830.0d-06, 825.0d-06, 
     $     794.000d-06, 827.000d-06, 826.000d-06, 841.0d-06, 847.0d-06,
     $     878.000d-06, 890.000d-06, 902.000d-06, 921.0d-06, 934.0d-06, 
     $     939.000d-06, 952.000d-06, 966.000d-06, 980.0d-06, 994.0d-06/

	I=NINT(Z)
	POT=POTELM(I)
	RETURN
	END
! *********************************************************
	DOUBLE PRECISION FUNCTION RAV2(Z,A)
!
! Auxiliary routine for SAMCS
!
c    mean squared radius in fm**2 
        IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)
	DIMENSION AHO(8)
	PARAMETER (PI=3.141592653589793227D+00)
	DATA AHO/0.,0.,1.77d0,1.78d0,1.69d0,0.,1.745d0,1.833d0/
!
	IF(Z.EQ.1.D0) THEN
	RAV2=0.7225D0
        RETURN
        ENDIF
        
	IF(Z.EQ.2.D0) then
        R=1.238
        B=0.38
	Z1=PI*B/R
	Z2=Z1*Z1
	RAV2=3.*R**2*(1.+7.*Z2/3.)/5.
        return  
        endif

	IF(Z.eq.6.D0) then
        R=2.214
        B=0.488
	Z1=PI*B/R
	Z2=Z1*Z1
	RAV2=3.*R**2*(1.+7.*Z2/3.)/5.
        return
        endif

        IF(Z.LE.8.D0) then
        IZ=Z
        B2=AHO(IZ)**2
        A02=(B2-.4824)*A/(A-1.)
	RAV2=A02*(Z-2)/Z+1.5*B2
        ELSE 
        A13=A**(1./3.)
        RAV2=(0.898*A13+1.198/A13)**2
        ENDIF
	RETURN
	END
! *********************************************************
	DOUBLE PRECISION FUNCTION RAV4(Z,A)
!
! Auxiliary routine for SAMCS
!
        IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)
	PARAMETER (PI=3.141592653589793227D+00)
        IF(Z.EQ.1.D0) THEN
        RAV4=5.D0*(RAV2(Z,A))**2/2.D0
        RETURN
        ENDIF  
        IF(Z.EQ.2.D0) THEN
        R=1.238
        B=0.38
	Z1=PI*B/R
	Z2=Z1*Z1
	RAV4=3.*R**4*(1.+6.*Z2+31.*Z2**2/3.)/7.	
        RETURN
        ENDIF
        IF(Z.EQ.6.D0) THEN
        R=2.214
        B=0.488
	Z1=PI*B/R
	Z2=Z1*Z1
	RAV4=3.*R**4*(1.+6.*Z2+31.*Z2**2/3.)/7.	
        RETURN
        ENDIF
        IF(Z.LE.8.D0) THEN
        RAV4=5.D0*(RAV2(Z,A))**2/3.D0  
        ELSE 
        A13=A**(1./3.)
        R=1.141*A13-1.047/A13
        R2=(0.898*A13+1.198/A13)**2
	B=R*SQRT(3.*(5.*R2/3./R**2-1.)/7.)/PI
	Z1=PI*B/R
	Z2=Z1*Z1
	RAV4=3.*R**4*(1.+6.*Z2+31.*Z2**2/3.)/7.	
        ENDIF
        RETURN
        END
c  *********************************************************
      subroutine RESTRICT_LOSS(z,step,p,e,dedx,de)
!
! Mean energy loss and fluctuations in the restricted energy range
!  below DCUTx
!
! Modified from GEANT routine  GLANDZ
!
C.    ******************************************************************
C.    *                                                                *
C.    *  Energy straggling using a Monte Carlo model.                  *
C.    *  It can be used with or without delta ray generation.          *
C.    *                                                                *
C.    *   It is a NEW VERSION of the model , which reproduces          *
C.    *   the experimental data rather well.                           *
C.    *                                                                *
! n.b we assume STEP is in [m] in the following
C.    *  Input : STEP  = current step-length (cm)                      *
C.    *                                                                *
C.    *  Output: DE    = actual energy loss (Gev)                      *
C.    *                 ( NOT the fluctuation DE/DX-<DE/DX> !)         *
C.    *                                                                *
C.    *  Author      : L.Urban                                         *
C.    *  Date        : 28.04.1988       Last update :  1.02.90         *
C.    *                                                                *
C.    ******************************************************************
C.
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      PARAMETER (MAXRND=100)
      DIMENSION RNDM(MAXRND),APOIS(3),NPOIS(3),FPOIS(3)
      PARAMETER (HMXINT=2**30)
!      PARAMETER (HMXINT=2**45)
      PARAMETER (ONE=1,HALF=ONE/2,ZERO=0)
      PARAMETER (RCD=0.40, RCD1=1.-RCD, PROBLM=0.01)
      PARAMETER( C1=4, C2=16)
*
*     ------------------------------------------------------------------
*
      amass = pmass
      emass = mass(1)
!
      POTI=1.6E-8*Z**0.9
      POTIL=LOG(POTI)
*
      IF(Z.GT.2.) THEN
         F2=2./Z
      ELSE
         F2=0.
      ENDIF
      F1=1.-F2
*
      E2 = 1.E-8*Z*Z
      E2L= LOG(E2)
      E1L= (POTIL-F2*E2L)/F1
      E1 = EXP(E1L)
*
*
      P2=P*P
      B2=P2/(E*E)
      BG2=P2/(AMASS*AMASS)
!      IF(ITRTYP.EQ.2) THEN
      if( ABS(iptyp) .eq. 1 ) then
         TM=P2/(E+EMASS)
!         IF(CHARGE.LT.0.) TM=TM/2.
         IF(pCHARGE.LT.0.) TM=TM/2.
         TM=TM-POTI
         IF(TM.GT.DCUTE) TM=DCUTE
      ELSE
         TM=EMASS*P2/(0.5*AMASS*AMASS+EMASS*E)
         TM=TM-POTI
         IF(TM.GT.DCUTM) TM=DCUTM
      ENDIF
*
*
* *** Protection against negative TM    ---------------------
*     TM can be negative only for heavy particles with  a very
*     low kinetic energy (e.g. for proton with T  100-300 kev)
      TM=MAX(TM,ZERO)
*
      W  = TM+POTI
      WW = W/POTI
      WWW= 2.*EMASS*BG2
      WL = LOG(WWW)
      CSB=STEP*RCD1*DEDX/(WL-POTIL-B2)
      APOIS(1)=CSB*F1*(WL-E1L-B2)/E1
      APOIS(2)=CSB*F2*(WL-E2L-B2)/E2
*
      IF(TM.GT.0.) THEN
         APOIS(3)=RCD*DEDX*STEP*TM/(POTI*W*LOG(WW))
      ELSE
         APOIS(1)=APOIS(1)/RCD1
         APOIS(2)=APOIS(2)/RCD1
         APOIS(3)=0.
      ENDIF
*
*    calculate the probability of the zero energy loss
*
      APSUM=APOIS(1)+APOIS(2)+APOIS(3)
      if( apsum .lt. 0d0 ) then
        iflag = -29
        go to 900
      end if
      IF(APSUM.LT.50. ) THEN
         PROB=EXP(-APSUM)
      ELSE
         PROB=0.
      ENDIF
*
*
*      do it differently if prob > problm  <====================
      IF(PROB.GT.PROBLM) THEN
         E0=1.E-8
         EMEAN=DEDX*STEP
         IF(TM.LE.0.) THEN
*      excitation only ....
            APOIS(1)=EMEAN/E0
*
            CALL GPOISS(APOIS,NPOIS,1)
            FPOIS(1)=NPOIS(1)
            DE=FPOIS(1)*E0
*
         ELSE
*         ionization only ....
            EM=TM+E0
            APOIS(1)=EMEAN*(EM-E0)/(EM*E0*LOG(EM/E0))
            CALL GPOISS(APOIS,NPOIS,1)
            NN=NPOIS(1)
            DE=0.
*
            IF(NN.GT.0) THEN
               RCORR=1.
               IF(NN.GT.MAXRND) THEN
                  RCORR=FLOAT(NN)/MAXRND
                  NN=MAXRND
*
               ENDIF
               W=(EM-E0)/EM
               CALL GRNDM(RNDM,NN)
               DO 10 I=1,NN
                  DE=DE+E0/(1.-W*RNDM(I))
   10          CONTINUE
               DE=RCORR*DE
*
            ENDIF
         ENDIF
         GOTO 999
      ENDIF
*
      IF(MAX(APOIS(1),APOIS(2),APOIS(3)).LT.HMXINT) THEN
         CALL GPOISS(APOIS,NPOIS,3)
         FPOIS(1)=NPOIS(1)
         FPOIS(2)=NPOIS(2)
         FPOIS(3)=NPOIS(3)
      ELSE
         DO 20 JPOIS=1, 3
            IF(APOIS(JPOIS).LT.HMXINT) THEN
               CALL GPOISS(APOIS(JPOIS),NPOIS(JPOIS),1)
               FPOIS(JPOIS)=NPOIS(JPOIS)
            ELSE
               CALL GRNDM(RNDM,2)
               FPOIS(JPOIS)=ABS(SQRT(-2.*LOG(RNDM(1)*ONE)
     +         *APOIS(JPOIS))*SIN(TWOPI*RNDM(2)*ONE)+APOIS(JPOIS))
            ENDIF
   20    CONTINUE
      ENDIF

*
*          Now we have all three numbers in REAL/DOUBLE
*          variables. REALK is actually an INTEGER that now may
*          exceed the machine representation limit for integers.
*
      DE=FPOIS(1)*E1+FPOIS(2)*E2
*
*     smear to avoid peaks in the energy loss (note: E1<<E2)
*
      IF(DE.GT.0.) THEN
         CALL GRNDM(RNDM,1)
         DE=DE+E1*(1.-2.*RNDM(1))
      ENDIF
*
      ALFA=1.
      REALK=0
      DEC=0.
*
      ANC=FPOIS(3)
      IF(ANC.GE.C2) THEN
         R=ANC/(C2+ANC)
         AN=ANC*R
         SN=C1*R
         CALL GRNDM(RNDM,2)
         RR=SQRT(-2.*LOG(RNDM(1)))
         PHI=TWOPI*RNDM(2)
         X=RR*COS(PHI)
         AK=AN+SN*X
         ALFA=WW*(C2+ANC)/(C2*WW+ANC)
         EA=AK*POTI*ALFA*LOG(ALFA)/(ALFA-1.)
         SA=SQRT(ABS(AK*ALFA*POTI*POTI-EA*EA/AK))
         AKL=(EA-C1*SA)/POTI
         IF(AK.LE.AKL) THEN
            X=RR*SIN(PHI)
            DEC=EA+SA*X
            REALK=AK+HALF-MOD(AK+HALF,ONE)
         ELSE
            ALFA=1.
         ENDIF
      ENDIF
      NN=NINT(FPOIS(3)-REALK)
      IF(NN.GT.MAXRND) THEN
         W=1.-ALFA/WW
         WW=POTI*ALFA
*
*     Here we take a gaussian distribution to avoid loosing
*     too much time in computing
*
         AVERAG=-LOG(1.-W)/W
         SIGMA =SQRT(NN*(1./(1.-W)-AVERAG**2))
         CALL GRNDM(RNDM,2)
         DEC   = DEC+WW*(NN*AVERAG+SIGMA*SQRT(-2.*LOG(RNDM(1)))*
     +           SIN(TWOPI*RNDM(2)))
         DE=DE+DEC
      ELSEIF(NN.GT.0) THEN
         W=1.-ALFA/WW
         WW=POTI*ALFA
         CALL GRNDM(RNDM,NN)
         DO 30 I=1,NN
            DEC=DEC+WW/(1.-W*RNDM(I))
   30    CONTINUE
         DE=DE+DEC
      ENDIF
*
900   continue
  999 END
c  *********************************************************
      subroutine ROTATE_SPIN(p)
!
! Rotate spin vector such that helicity is conserved
!  p is momentum after scattering
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
!
      real*8 p1(3),cr2(3),cr3(3),p(3)
!
      call VMAG(pp,pm)
      h = HELICITY(pol)        ! initial helicity
      hp = SQRT(1. - h*h)
      q = SQRT( p(1)**2 + p(2)**2 + p(3)**2 )
      do i=1,3
      p1(i) = p(i) / q
      end do
      call DCROSS(p1,pol,cr2)
      call DCROSS(cr2,p1,cr3)
      call VMAG(cr3,cr3m)
!
      if( cr3m .gt. 0. ) then
        do i=1,3
        pol(i) = h*p1(i) + hp*cr3(i)/cr3m
        end do
      end if
!
      end
c  *********************************************************
      subroutine RUTHSCAT(omega,chia2,thet,phi)
c
c  Determine scattering angle from Rutherford distribution
c
c  Derived form GEANT program  GMCOUL
c
c     omega : number of scatters
c     chia2 : Coulomb screening angle (squared)
c
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
c
c  Get actual # of scatters from Poisson distribution
      omega0 = 1.167d0 * omega
      nsc = POIDEV(omega0,irnarg)    ! Num Rec
!
      if( nsc .gt. 0 ) then
        sumx = 0.d0
        sumy = 0.d0
!
        do i=1,nsc
        rn = RAN1(irnarg)
        thet = SQRT( chia2 * ((1./rn) - 1.) )
        phi = twopi * RAN1(irnarg)
! Projected angles have scaler sums
        sumx = sumx + thet*COS(phi)
        sumy = sumy + thet*SIN(phi)
        end do
!
        thet = SQRT( sumx**2 + sumy**2 )
        if( sumx .ne. 0. ) then
          phi = ATAN2(sumy,sumx)
        else 
          phi = 0.
        end if
      else
        thet = 0.
        phi = 0.
      end if
c
      return
      end
! *********************************************************
	DOUBLE PRECISION FUNCTION R2P(IP)
!
! Auxiliary routine for SAMCS
!
C      ------------------------------------------------------------------------------------
C       Average squared radius of projectile in fermi 
C      ------------------------------------------------------------------------------------
c     1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18   19  20   21
c     p   n  pi+ pi-  K+  K- mu+ mu-  g  e-  e+  ap  pi0  d   t  He3 He4 num nuam nue nuae
C      -------------------------------------------------------------------------------------
C       Proton and pion data - from Review of Particle Physics (Phys. Rev. D66 (2002) 010001)
C       Kaon data - from Phys. Rev. Lett 45 (1980) 232 and Phys. Lett. B178 (1986) 435 
C       Light nuclei data  - from Atomic Data and Nuclear Data Tables 36 (1987) 495-536
C      -------------------------------------------------------------------------------------
        IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)
        DIMENSION R2(21)
        DATA R2/0.7225D0,0.D0,2*0.452D0,2*0.31D0,5*0.D0,0.7225D0,0.D0,
     !   4.433D0,2.9584D0,3.6062D0,2.8258D0,4*0.D0/
        R2P=R2(IP)
        RETURN
        END
! *********************************************************
      SUBROUTINE SAMCS(EKIN,STEP,IP,NELT,AD,ZD,XD,irns,
     &     TETM,DELION)
!
! Correlated scattering and straggling using model of Sergei Striganov
!
	IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)   
c   input EKIN  - kinetic energy in GeV
c   input STEP  - step length in g/cm**2
c   input IP    - index of projectile   
c   input NELT  - number of elements in mixture
c   input AD(NELT) - array of atomic masses
c   input ZD(NELT) - array of atomic numbers
c   input XD(NELT) - array of relative weights
c-------------------------------------------------
c   output TETM - angle in rad (due to Coulomb scattering on nuclear field and atomic electrons)
c   output DELION - energy loss in GeV (due to scattering on atomic electrons)
c-------------------------------------------------
c   SAMCS calls CERNLIB routine/function RNDM, RNPSSN  
!!!!! this stuff has single precision in CERNLIB !!!!!!!!

!	REAL RNDM
      common/rans/iseed
	PARAMETER (NMIX=100)
        DIMENSION AD(NMIX),ZD(NMIX),XD(NMIX),ASING(NMIX)
	DIMENSION XIC2(NMIX),XIA2(NMIX),AMY2(NMIX)
        DIMENSION XIC2I(NMIX),XII2(NMIX),SIGMA(NMIX)
        DIMENSION RI1(NMIX),RI2(NMIX),ELEMIN(NMIX)
        DIMENSION ARA(NMIX),ALI(NMIX),FDI(NMIX)
c    projectile properties 
c ------------------------------------------------------------------------------------------    
c IP  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18   19  20   21
c     p   n  pi+ pi-  K+  K- mu+ mu-  g  e-  e+  ap  pi0  d   t  He3 He4 num nuam nue nuae
c ------------------------------------------------------------------------------------------
c    PMASS     - mass in GeV
c    ZPR       - charge
c    R2P(IP)   - average squared radius in fm**2. See function R2P(IP)
c    SPII      - spin 
C    -------------------------------------------------------------------------------------

	PARAMETER (NPROJ=21)
        DIMENSION PMASS(NPROJ),ZPR(NPROJ),SPII(NPROJ)

	PARAMETER (SME=0.51099906D-03)
	PARAMETER (PI=3.141592653589793227D+00)
	PARAMETER (CEUL=5.772156649D-1)
	PARAMETER (DELT=3.D-2)
	PARAMETER (DELT2=1.D-1)
	PARAMETER (DELTI=1.D-1)
	PARAMETER (DELTI2=1.D-1)
        DATA PMASS/0.938272310D+00,0.939565630D+00,2*0.139569950D+00,    
     &    2*0.493677000D+00,2*0.105658389D+00,0.000000000D+00,    
     &    2*0.510999070D-03,0.938272310D+00,0.134976400D+00,
     &    0.187783794D+01,  0.281740357D+01,0.281611025D+01, 
     &    0.375567588D+01,  4*0.0D+00/   
        DATA ZPR/1.D0,0.D0,1.D0,-1.D0,1.D0,-1.D0,1.D0,-1.D0,0.D0,
     &        -1.D0,1.D0,-1.D0,0.D0,1.D0,1.D0,2.D0,2.D0,4*0.D0/
        DATA SPII/2*5.D-1,4*0.D0,2*5.D-1,0.D0,3*5.D-1,9*0.D0/
!
      iseed = irns
        SPIN=SPII(IP)
	AMASS=PMASS(IP)
	ETOT=EKIN+AMASS
	PLAB=SQRT(ETOT**2-AMASS**2)
	BET2=(PLAB/ETOT)**2	
	EMAX=2.*SME*PLAB**2/(SME**2+2.*SME*ETOT+AMASS**2)
	C1=2.D0*SME/PLAB**2
	C2=ETOT*EMAX/PLAB**2
        C3=SPIN*(EMAX/ETOT)**2
        XM=0.D0
	YM=0.D0
	ZM=0.D0
	WM1=0.D0
	T2=0.D0
	T4=0.D0
        ZMI=0.D0 
        T2I=0.D0
	T4I=0.D0
	XMIXT=0.D0
        YMIXT=0.D0
	ZMIXT=0.D0
        DO N=1,NELT
	   TX=STEP*XD(N)
	   XIC2(N)=1.5692D-7*ZD(N)*ZD(N)*TX*ZPR(IP)**2
     &             /(AD(N)*PLAB**2*BET2)
	   XIA2(N)=2.00694D-11*(1.D0+(ZD(N)*ZPR(IP)/75.12D0)
     &              **2/BET2)*ZD(N)**(2.D0/3.D0)/PLAB**2
c    r^2_nuc from formfactor of nuclear	
c	   RA2=(RAV2(ZD(N),AD(N))+R2P(IP))*5.06765D0**2
c	   AMY2(N)=3.D0/PLAB**2/RA2
c           ASING(N)=SNGL(AMY2(N))
	   XM=XM+XIC2(N)
	   YM=YM+XIC2(N)*LOG(1.D0/XIA2(N))
	   ZM=ZM+XIC2(N)/XIA2(N)
           XIA2P=XIA2(N)*PLAB**2
	   CALL AVAN(IP,ZD(N),AD(N),XIA2P,T2C,T4C) 
	   T2=T2+XIC2(N)*T2C
	   T4=T4+XIC2(N)*T4C/PLAB**2
c    r^2_nuc from average angle squared
           AMY2(N)=EXP(T2C+1.+CEUL)*XIA2(N)
 	   ASING(N)=AMY2(N)
           WM1=WM1+XIC2(N)/AMY2(N)
	   ELEMIN(N)=POT(ZD(N))*1.D-2
           ALI(N)=ELEMIN(N)/EMAX
           IF(ALI(N).LE.1.D0) THEN
              ARA(N)=2.D0*ELEMIN(N)*(SIPAR(ZD(N)))**2/SME
              XII2(N)=(SME/SIPAR(ZD(N))/PLAB)**2
              XIC2I(N)=XIC2(N)/ZD(N)
              ZMI=ZMI+XIC2I(N)/XII2(N)
              RI1(N)=1.D0-LOG((ARA(N)+1.D0)/ARA(N))-1.D0/(1.D0+ARA(N))
              RI2(N)=2.D0*LOG(ARA(N)+1.D0)+1.D0/(1.D0+ARA(N))-1.D0
              XMIXT=XMIXT+XIC2I(N)
              YMIXT=YMIXT+XIC2I(N)*(LOG(XII2(N))+RI1(N))
              YMAX=1.D0
              CALL T2ION(ALI(N),YMAX,PLAB,AMASS,SPIN,RES2)
              CALL T4ION(ALI(N),YMAX,PLAB,AMASS,SPIN,RES4)
              T2I=T2I+XIC2I(N)*(LOG(YMAX/ALI(N))-RI1(N)
     &            +LOG(ARA(N))+RES2)
              T4I=T4I+XIC2I(N)*(C1*EMAX*(RES4+ALI(N))-XII2(N)*RI2(N))
           ENDIF
	ENDDO

C Nuclear field  Z**2

c        AMYMAX1=MINVAL(ASING(1:NELT))
c  if there are no MAXVAL as intrinistic function 
c  use CERNLIB LVSDMX routine 
!        LMIN=LVSDMI(ASING,NELT,1)
        LMIN = ILVSDMI(ASING,NELT,1)
        AMYMAX=ASING(LMIN)

	CMY=T4/2.D0/T2**2
	IF(CMY.LE.1.D-2) THEN
	   THTOT=T2*FNEAR(CMY)
	ELSE   
	   IF(ZM.LT.3.D1) THEN
	      TMY2=0.D0
	      FG2=0.D0
	   ELSE
	      BMY=(LOG(8.D0*DELT)+LOG(XM)-1.D0+YM/XM)/2.D0
	      IF(BMY.GT.1.D0) THEN
		 XX=BMY*(1.D0+LOG(BMY)/(BMY-1.D0))*
     &	         (1.D0+SQRT(1.D0-1.D0/BMY))/2.D0
		 YMY=8.D0*DELT*XX*XX
		 TMY21=YMY*XM
		 TMY22=AMYMAX*DELT2
		 TMY2=MIN(TMY21,TMY22)
	      ELSE
		 TMY2=EXP(3.D0-YM/XM)
	      ENDIF
	      DMY=XM*LOG(TMY2)+YM-XM-TMY2*WM1
	      C=XM*TMY2*(1.D0-TMY2*WM1/2.D0/XM)/2./DMY**2
	      FG2=DMY*FNEAR(C)
	   ENDIF
	   THTOT=FG2
           DO ND=1,NELT
	      IF(TMY2.GT.0.D0) THEN
                 QMIN=PLAB*SQRT(TMY2)   
		 PL=XIC2(ND)*SIGD(IP,ZD(ND),AD(ND),QMIN)*PLAB**2
   	         DTMY=TMY2/AMY2(ND)
	      ELSE
		 AXIA=XIA2(ND)/AMY2(ND)
		 PL=XIC2(ND)/XIA2(ND)*(1.D0-AXIA*(CEUL-LOG(AXIA)))
	      ENDIF
	      CALL MCSDIS(PL,TMY2,PLAB,AD(ND),ZD(ND),IP,
     &                    XIA2(ND),TET2)
!	      G1=2.D0*PI*DBLE(RNDM(-1.))
	      G1=2.D0*PI* RAN1(iseed)
	      THTOT=THTOT+TET2-2.D0*SQRT(THTOT)*SQRT(TET2)*SIN(G1)
	   ENDDO
	ENDIF

C  Electrons  Z
        DELION=0.D0
        TEL2  =0.D0
c        YMIN=MAXVAL(ALI(1:NELT))
c  if there are no MAXVAL as intrinistic function 
c  use CERNLIB LVSDMX routine 
!        LMAX=LVSDMX(ALI,NELT,1)
        LMAX = ILVSDMX(ALI,NELT,1)
        YMIN=ALI(LMAX)
        IF(YMIN.LE.1.D0) THEN 
        CMYI=T4I/2.D0/T2I**2
	IF(CMYI.LE.1.D-3) THEN
	   TEL2=T2I*FNEAR(CMYI)
	   PSI=XMIXT/C1
	   CAPPAI=PSI/EMAX
	   IF(CAPPAI.GE.1.D0) THEN
	      AM2=PSI*EMAX*(1.D0-BET2/2.D0+C3/3.D0)
	      AM3=PSI*EMAX**2*(5.D-1-BET2/3.D0+C3/4.D0)
	      DELION=DFL(AM2,AM3)
	   ELSE
	      EPOR=PSI
	      XEPOR=EPOR/EMAX
	      XEE2=(EPOR/ETOT)**2
	      AM2=PSI*EPOR*(1.D0-BET2*XEPOR/2.D0+SPIN*XEE2/3.D0)
	      AM3=PSI*EPOR**2*(5.D-1-BET2*XEPOR/3.D0+SPIN*XEE2/4.D0)
	      DLAR=PSI*(-LOG(XEPOR)-BET2*(1.-XEPOR)
     &             +C3*(1.D0-XEPOR**2)/2.D0)
	      DELION=DFL(AM2,AM3)-DLAR
	      PLI=PSI/EPOR*(1.D0-XEPOR+BET2*XEPOR*LOG(XEPOR)+
     &	      C3*(1.D0-XEPOR)*XEPOR)
!	      CALL RNPSSN(SNGL(PLI),NION,IERR)
            nion = POIDEV(pli,iseed)
	      IF(NION.GT.0) THEN
		 DO NI=1,NION
		    DEL1=SELEI(SPIN,EPOR,XEPOR,BET2,EMAX,ETOT)
		    DELION=DELION+DEL1
		 ENDDO
	      ELSE
	      ENDIF
	   ENDIF
        ELSE
           IF(ZMI.LT.3.D1) THEN
	      PLI=0.D0
	      DO N=1,NELT
		 FDI(N)=1.D0-ALI(N)+BET2*ALI(N)*LOG(ALI(N))+
     &                  C3*ALI(N)*(1.D0-ALI(N))
		 SIGMA(N)=XIC2I(N)*(ARA(N)/(1.D0+ARA(N))+
     &                         FDI(N)/ARA(N))/XII2(N)
		 PLI=PLI+SIGMA(N)
	      ENDDO
!	      CALL RNPSSN(SNGL(PLI),NION,IERR)
            nion = POIDEV(pli,iseed)
	      TEL2=0.D0
	      DELION=0.D0
	      DO NI=1,NION
                 SK=0.D0  
!                 GT0=DBLE(RNDM(-1.))
                 GT0 = RAN1(iseed)
		 DO IK=1,NELT
		    SK=SK+SIGMA(IK)/PLI
		    IF(GT0.LE.SK) THEN
		       NM=IK
                    GO TO 1
		    ELSE
		    ENDIF
		 ENDDO
 1               FATOM1=ARA(NM)/(1.D0+ARA(NM))
		 FATOM2=FATOM1/(FATOM1+FDI(NM)/ARA(NM))
!		 GT1=DBLE(RNDM(-1.))
		 GT1 = RAN1(iseed)
		 IF(GT1.LE.FATOM2) THEN
!		    GT2=DBLE(RNDM(-1.))
		    GT2 = RAN1(iseed)
		    A1=XII2(NM)
		    A2=A1*(1.D0+ARA(NM))
		    T2HI=A1/(1.D0-GT2*(1.D0-A1/A2))-A1
!		    GN=2.D0*PI*DBLE(RNDM(-1.))
		    GN = 2.D0*PI*RAN1(iseed)
		    TEL2=T2HI+TEL2-2.D0*SQRT(T2HI)*SQRT(TEL2)*SIN(GN)
	         ELSE
		    DEL1=SELEI(SPIN,ELEMIN(NM),ALI(NM),BET2,EMAX,ETOT)
		    DELION=DELION+DEL1
		    RATI=DEL1/EMAX
		    T2HI=C1*DEL1*(1.D0-RATI)/(1.D0-C2*RATI)**2
!		    GN=2.D0*PI*DBLE(RNDM(-1.))
		    GN = 2.D0*PI*RAN1(iseed)
		    TEL2=T2HI+TEL2-2.D0*SQRT(T2HI)*SQRT(TEL2)*SIN(GN)
                if(t2hi.lt.0.d0) write(8,*) 'd/emax,t2hi',rati,t2hi
		 ENDIF
	      ENDDO
	   ELSE
	      BMYI=(LOG(8.D0*DELTI)+LOG(XMIXT)-YMIXT/XMIXT)/2.D0
	      IF(BMYI.GT.1.D0) THEN
		 XXI=BMYI*(1.D0+LOG(BMYI)/(BMYI-1.D0))*
     &	         (1.D0+SQRT(1.D0-1.D0/BMYI))/2.D0
		 YMYI=8.D0*DELTI*XXI*XXI
		 TMY2A=YMYI*XMIXT             ! t4/2/t2=delti     
		 TMY2B=C1*EMAX*DELTI2         ! D_MIN/E_MAX << 1   
		 TMY2C=XMIXT                  ! CAPPA=1 
		 TMY2I=MIN(TMY2A,TMY2B,TMY2C)
	      ELSE
		 TMY2I=EXP(2.D0+YMIXT/XMIXT)
	      ENDIF
	      YMAX=TMY2I/C1/EMAX
	      DMYI2=0.D0
	      DMYI4=0.D0
              DO N=1,NELT
	      YI=TMY2I/XII2(N)
	      CALL T2ION(ALI(N),YMAX,PLAB,AMASS,SPIN,RES2)
	      CALL T4ION(ALI(N),YMAX,PLAB,AMASS,SPIN,RES4)
	      DMYI2=DMYI2+XIC2I(N)*(LOG(YI)-RI1(N)+RES2)
              DMYI4=DMYI4+XIC2I(N)*(C1*EMAX*(RES4+ALI(N))
     &              -XII2(N)*RI2(N))
	      ENDDO
	      CI=DMYI4/2.D0/DMYI2**2
	      TEL2=DMYI2*FNEAR(CI)
              EPOR=TMY2I/C1
	      XEPOR=EPOR/EMAX
	      XEE2=(EPOR/ETOT)**2
              PSI=XMIXT/C1
	      AM2=PSI*EPOR*(1.D0-BET2*XEPOR/2.D0+SPIN*XEE2/3.D0)
	      AM3=PSI*EPOR**2*(5.D-1-BET2*XEPOR/3.D0+SPIN*XEE2/4.D0)
	      DLAR=PSI*(-LOG(XEPOR)-BET2*(1.-XEPOR)
     &             +C3*(1.D0-XEPOR**2)/2.D0)
	      DELION=DFL(AM2,AM3)-DLAR
	      PLI=XMIXT/TMY2I*(1.D0-XEPOR+BET2*XEPOR*LOG(XEPOR)+
     &	      C3*(1.D0-XEPOR)*XEPOR)
!	      CALL RNPSSN(SNGL(PLI),NION,IERR)
            nion = POIDEV(pli,iseed)
              IF(NION.GT.0) THEN
	      DO NI=1,NION
		 DEL1=SELEI(SPIN,EPOR,XEPOR,BET2,EMAX,ETOT)
		 DELION=DELION+DEL1
		 RATI=DEL1/EMAX
		 T2HI=C1*DEL1*(1.D0-RATI)/(1.D0-C2*RATI)**2
!		 GN=2.D0*PI*DBLE(RNDM(-1.))
		 GN = 2.D0*PI*RAN1(iseed)
		 TEL2=T2HI+TEL2-2.D0*SQRT(T2HI)*SQRT(TEL2)*SIN(GN)
	      ENDDO
	      ELSE
	      ENDIF
	   ENDIF
	ENDIF
        ENDIF
C  TOTAL
! 	G1=2.D0*PI*DBLE(RNDM(-1.))
 	G1 = 2.D0*PI*RAN1(iseed)
	THTOT=THTOT+TEL2-2.D0*SQRT(THTOT)*SQRT(TEL2)*SIN(G1)
	TETM=SQRT(THTOT)
	irns = iseed
	RETURN
	END
c  *********************************************************
	subroutine SAMCS_DRIVER(hstp,pm)
!
! Interface ICOOL with SAMCS for scattering and straggling
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
	dimension pms(3),ad(100),zd(100),xd(100)
	integer ndx(11)
	data ndx/12,6,4,8,10,0,11,7,3,5,1/
!
      pm0 = pm          ! save incoming momentum
! Find mean energy loss
      call DEDXMEAN(hstp,pm,demean)
	dedxm = demean / hstp
!
      e0 = SQRT(pm**2 + pmass**2)   ! initial total energy
      gam = e0 / pmass
      bet = pm / pmass / gam
      bg = bet * gam
	dm = mparm(2,krad)
!
! Get correlated scattering and straggling
      ekin = e0 - pmass
	step = hstp * dm * 1d2        
	ip = ndx(iptyp+6)
	nelt = nelcomp(krad)
	do i=1,nelt
	ad(i) = aelcomp(i,krad)
	zd(i) = zelcomp(i,krad)
	xd(i) = felcomp(i,krad)
	end do
      irns = irnarg
!
      call SAMCS(EKIN,STEP,IP,NELT,AD,ZD,XD,irns,TETM,DELION)
!
      irnarg = irns
! Apply the computed energy loss to the actual particle 4-vector
      eloss = demean + delion    ! actual kinetic energy loss
!
      ekin = e0 - pmass - eloss        ! remaining kinetic energy
      engmin = SQRT(pzmintrk**2 + pmass**2) - pmass   ! minimum allowed
!
      if( ekin .gt. engmin ) then  ! apply loss and continue
        e = e0 - eloss            
!
      else                         ! bring particle to rest
        iflag = -43
        go to 900
      end if
!
      if( dedxm .gt. 0. ) then
        hnextde = MAX( epsf*ekin,engmin ) / dedxm
      else
        hnextde = zstep * scalestep
      end if
!
! Scale down the momentum vector
      pm = SQRT(e**2 - pmass**2)
      do i=1,3
      pp(i) = pp(i) * pm / pm0
      end do
c
! Apply the computed scattering to the energy-corrected 4-vector
      sthms = SIN(tetm)
      phims = RANV(0.d0,twopi)
      pms(1) = pm * sthms * COS(phims)   
      pms(2) = pm * sthms * SIN(phims)
      pms(3) = pm * COS(tetm)
!
      cth = pp(3) / pm        ! initial direction in LAB
      pt = SQRT(pp(1)**2 + pp(2)**2)
      sth = pt / pm
      if( pt .gt. 0. ) then
        cph = pp(1) / pt
        sph = pp(2) / pt
      else
        cph = 1.
        sph = 0.
      end if
! New direction in LAB frame
      pp(1) = pms(1)*cth*cph - pms(2)*sph + pms(3)*sth*cph
      pp(2) = pms(1)*cth*sph + pms(2)*cph + pms(3)*sth*sph
      pp(3) = -pms(1)*sth + pms(3)*cth
!
      call VMAG(pp,pm)
!
900   continue
      end
c  *********************************************************
      subroutine SCATTER(hstp,pavg)
c
c  Correct the momentum vector for multiple scattering
!
!  Use pavg to obtain new multiple scattering angle
! 
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
c
      real*8 pscat(3),pms(3)
!      real*8 delp(3),dpm
!      LOGICAL FIRST
!      DATA FIRST/.TRUE./
!
!      IF(FIRST)THEN
!        OPEN(UNIT=63,FILE='SCATTER.DAT',STATUS='UNKNOWN')
!        FIRST = .FALSE.
!      END IF
      bet = pavg/pmass / SQRT(1.+(pavg/pmass)**2)
      bet2 = bet * bet
      emev = 1e3 * SQRT(pavg**2 + pmass**2 )
c
c  Get multiple scattering angles
      if( scatlev .eq. 1 ) then      ! Rossi-Greisen
        rl = mparm(4,krad) / 100.    ! [m]
        arg = hstp / rl
        thmsw = 21.2e-3 / (bet*pavg) * SQRT(arg) * ABS(pcharge)/ee
        hnextms = (epsf*pavg*bet/21.2e-3)**2 * rl
!
        thms = GAUSS(0.d0,thmsw)
        phims = RANV(0.d0,twopi)
c
      else if( scatlev .eq. 2 ) then    ! Highland
        rl = mparm(4,krad) / 100.       ! [m]
        arg = hstp / rl
        if( arg .le. 0. ) then
          iflag = -27
          go to 900
        end if
        thmsw = 19.2e-3/(bet*pavg) * ABS(pcharge)/ee
     &    * SQRT(arg) *(1.+0.038*LOG(arg))
        hnextms = (epsf*pavg*bet/19.2e-3)**2 * rl
!
        thms = GAUSS(0.d0,thmsw)
        phims = RANV(0.d0,twopi)
c ---------------------------------------------------------
      else if( scatlev .gt. 2 ) then
! Geant3 manual p. 233, eq. 3
!  chiccmcs has units [MeV^2 / cm]
        chic2 = chiccmcs / (pavg*1e3)**2 * hstp*100. 
     &    * (ABS(pcharge)/ee/bet)**2
!
        zxmcs = 0.
! Ref. Geant3 manual, p. 234
        if( scatlev.eq.6 .and. ABS(iptyp).ne.1 ) then  ! Fano, heavy
          do i=1,nelcomp(krad)
          zxmcs = zxmcs + felcomp(i,krad) * 
     &     zelcomp(i,krad)**2 * 
     &     LOG( 1. + 3.34 * (zelcomp(i,krad)/137.036/bet)**2 ) 
          end do
        else
          do i=1,nelcomp(krad)
          zxmcs = zxmcs + felcomp(i,krad) * 
     &     zelcomp(i,krad) * (zelcomp(i,krad)+1.) * 
     &     LOG( 1. + 3.34 * (zelcomp(i,krad)/137.036/bet)**2 ) 
          end do
        end if
!
        chia2 = 2.016e-5 * facmms /
     &    (pavg*1e3)**2 * EXP( (zxmcs-zemcs)/zsmcs )
!
        omega = chic2 / (1.167 * chia2)
!
        bc = omega * bet2 / (100.*hstp)      ! [cm^-1]
        tbethe = (emev * bet2)**2 /
     &    ( chiccmcs * LOG(bc * emev**2 * bet2 / chiccmcs) )
        hnextms = epsf * tbethe/100.        ! [m]
c .........................................................
        if( scatlev .eq. 3 ) then     ! Lynch-Dahl
          v = 25. * omega    ! v = 0.5*omega/(1-f) ; f=0.98
          arg = v + 1.
          if( arg .le. 0. ) then
            iflag = -27
            go to 900
          end if
          thmsw = 0.510*chic2 * ( (1.+v)/v * LOG(arg)-1. )
          thmsw = SQRT(thmsw)
!
          thms = GAUSS(0.d0,thmsw)
          phims = RANV(0.d0,twopi)
c .........................................................
        else if( scatlev .eq. 4 ) then    ! Moliere
          if( omega .lt. 20. ) then    ! Moliere not valid if om < 20
            call RUTHSCAT(omega,chia2,thms,phims)
          else
c
            call MOLIERE(omega,chic2,thms)
            phims = RANV(0.d0,twopi)
c
          end if
c .........................................................
        else if( scatlev .eq. 5 ) then    ! Rutherford
c
            call RUTHSCAT(omega,chia2,thms,phims)
c .........................................................
        else if( scatlev .eq. 6 ) then    ! Fano-modified Moliere
          if( omega .lt. 20. ) then    ! Moliere not valid if om < 20
            call RUTHSCAT(omega,chia2,thms,phims)
          else
c
            call FANO(omega,chic2,pavg,thms)
            phims = RANV(0.d0,twopi)
c
          end if
c .........................................................
        else if( scatlev .eq. 7 ) then    ! Tollestrup-modified Moliere
          if( omega .lt. 20. ) then    ! Moliere not valid if om < 20
            call RUTHSCAT(omega,chia2,thms,phims)
          else
c
            call TOLLESTRUP(omega,chic2,chia2,hstp,pavg,thms)
            phims = RANV(0.d0,twopi)
c
          end if
        end if
      end if
! ---------------------------------------------------------
! Update actual particle 4-vector
      call VMAG(pp,pm)
      sthms = SIN(thms)
!   Particle frame
      pms(1) = pm * sthms * COS(phims)   
      pms(2) = pm * sthms * SIN(phims)
      pms(3) = pm * COS(thms)
!
      cth = pp(3) / pm
      pt = SQRT(pp(1)**2 + pp(2)**2)
      sth = pt / pm
      if( pt .gt. 0. ) then
        cph = pp(1) / pt
        sph = pp(2) / pt
      else
        cph = 1.
        sph = 0.
      end if
!   LAB frame
      pscat(1) = pms(1)*cth*cph - pms(2)*sph + pms(3)*sth*cph
      pscat(2) = pms(1)*cth*sph + pms(2)*cph + pms(3)*sth*sph
      pscat(3) = -pms(1)*sth + pms(3)*cth
!
      do i=1,3
      pp(i) = pscat(i)
      end do
!
      if( spintrk .gt. 0 ) then   
! Rotate the spin such that helicity is conserved
        call ROTATE_SPIN(pp)
!
        if( spinmatter.eq.2 .and. ABS(iptyp).eq.2 ) then
! Simulate spin flip in mu-Z Coulomb scattering
!
          call FLIP_COULOMB( pavg,hstp,pr )
!
          rn = RAN1(irnarg)
          if( rn .lt. pr ) then
            do i=1,3        ! flip spin direction
            pol(i) = -pol(i)
            end do
          end if
        end if
      end if      ! end of test on spintrk
!
900   continue
      return
      end
! *********************************************************
      REAL*8 FUNCTION SELEI(SPIN,EPOR,R,B2,EMAX,ET)
!
! Auxiliary routine for SAMCS
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)
      common/rans/iseed
!      REAL RNDM
      DE1=1.D0-R
! 1    GRP=DBLE(RNDM(-1.))
 1    GRP = RAN1(iseed)
      SELEI=EPOR/(1.D0-GRP*DE1)
!      G2=DBLE(RNDM(-1.))
      G2 = RAN1(iseed)
      FOT=1.D0-B2*SELEI/EMAX+SELEI*SELEI*SPIN/ET**2
      IF(G2.GT.FOT) GO TO 1
      RETURN
      END
! *********************************************************
      SUBROUTINE SEXT(X,EXT)
!
! Auxiliary routine for SAMCS
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)
      COMMON/CSF/Z,A,IP,MOM
      PARAMETER (SC=5.06765D0)
      EXT=FORMP(X,IP)*FORMN(X,Z,A)*2.D0/X**3*X**MOM 
      RETURN
      END 

! *********************************************************
      SUBROUTINE SGINT1(SEXT,A,B,EPS,GINT1)
!
! Auxiliary routine for SAMCS
!
C....................................................
C     GAUSSIAN INTEGRATION
C     REVISION: 01-DEC-1998 BY NVM
C....................................................
      IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)

      EXTERNAL SEXT

      DIMENSION U(10),W(10)

      DATA M/4/,U/.861136311594D0,.339981043585D0,8*0.D0/
      DATA      W/.347854845137D0,.652145154863D0,8*0.D0/
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      N=1
      OTB=0.D0
 1    Y=OTB
      OTB=0.D0
      D=(B-A)*.5D0/N
      DO I=1,N
         DO K=1,M
            K1=M+1-K
            UK=U(K)
            WK=W(K)
            IF(K.GE.K1) THEN
               UK=-U(K1)
               WK=W(K1)
            END IF
            XD=A+D*(I+I-1-UK)
            CALL SEXT(XD,EXT)
            OTB=OTB+WK*EXT
         END DO
      END DO

      OTB=OTB*D
      N=2*N
      GINT1=OTB
      IF(EPS.LE.0.D0) RETURN
      IF(ABS(OTB-Y).GT.ABS(EPS*OTB)) GO TO 1

      RETURN
      END
! *********************************************************
      DOUBLE PRECISION FUNCTION SIGD(IP1,Z1,A1,QMIN)

!
! Auxiliary routine for SAMCS
!
C      ------------------------------------------------------------------------------------
C       SIGD*XIC2*PLAB**2 - Average number of Couomb interaction for momentum transfer
C       greater than QMIN 
C       QMIN - threshold of discrete interactions in (GeV/c)
C      ------------------------------------------------------------------------------------
C       Z1,A1 - charge and atomic mass of target  
C      -------------------------------------------------------------------------------------

      IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)
      COMMON/CSF/Z,A,IP,MOM
      DIMENSION S(100,8)
      DIMENSION IPS(21),ISP(8)
      DIMENSION ATBL(100)
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      LOGICAL FIRSTCALL/.TRUE./
      SAVE FIRSTCALL,S      
      PARAMETER (SC=5.06765D0)
      EXTERNAL SEXT
      DATA IPS/1,0,2,2,3,3,4,4,0,4,4,1,0,5,6,7,8,4*0/
      DATA ISP/1,3,5,7,14,15,16,17/
      DATA ATBL/1.00797,4.0026,6.939,9.0122,10.811,12.01115,14.0067, 15
     *.9994,18.9984,20.183,22.9898,24.312,26.9815,28.088,30.9738, 32.064
     *,35.453,39.948,39.102,40.08,44.956,47.90,50.942,51.998, 54.9380,55
     *.847,58.9332,58.71,63.54,65.37,69.72,72.59,74.9216, 78.96,79.808,8
     *3.80,85.47,87.62,88.905,91.22,92.906,95.94,99.0, 101.07,102.905,10
     *6.4,107.87,112.4,114.82,118.69,121.75,127.60, 126.9044,131.30,132.
     *905,137.34,138.91, 140.12,140.907,144.24,147.,150.35,151.98,157.25
     *,158.924,162.50, 164.930,167.26,168.934,173.04,174.97,178.49,180.9
     *48,183.85, 186.2,190.2,192.2,195.08,196.987,200.59,204.37,207.19,2
     *08.980, 210.,210.,222.,223.,226.,227.,232.036,231.,238.03,237.,242
     *., 243.,247.,247.,248.,254.,253./
!
      IF ( FIRSTCALL ) THEN
         EPS=1.D-4  
         NZMAX=100
         IPMAX=8
         MOM=0 
         DO NZ=1,NZMAX 
            Z=1.D0*NZ
            A=ATBL(NZ)
            DO NIP=1,IPMAX
               IP=ISP(NIP)
               R2=MAX(R2P(IP),RAV2(Z,A)) ! fm**2
               Q1=SQRT(3.D-1/R2)/SC ! lower limit  GeV
               QMAX=SQRT(1.D3/R2)/SC ! upper limit  GeV
               CALL SGINT1(SEXT,Q1,QMAX,EPS,S(NZ,NIP))
            ENDDO
         ENDDO
         FIRSTCALL = .FALSE.
      IP=IP1
      Z=Z1
      A=A1
      ENDIF
      R2=RAV2(Z1,A1)
      R2M=MAX(R2P(IP1),R2) ! fm**2
      Q1=SQRT(3.D-1/R2M)/SC ! GeV
      Y=(QMIN/Q1)**2
      IZ=Z1
      IPP=IPS(IP1)  
      A4=R2**2*(1.D0+3.D0*RAV4(Z1,A1)/R2**2/5.D0)/36.D0+
     !   R2*R2P(IP1)/9.D0+5.D0*R2P(IP1)**2/72.D0
      QG2=(QMIN*SC)**2 
      SIGD=(1.D0-Y+(R2+R2P(IP1))*QG2*LOG(Y)/3.D0
     ! +A4*QG2**2*(1./Y-1.D0))/QMIN**2+S(IZ,IPP)
      RETURN
      END
! *********************************************************
	REAL*8 FUNCTION SIPAR(Z)
!
! Auxiliary routine for SAMCS
!
	IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N) 
	DIMENSION ARI(5)
	DATA ARI/282.4,265.8,418.6,571.4,724.2/
	IF(Z.LE.4.D0) THEN
	   NLOW=Z
	   SIP=ARI(NLOW)
	ELSE
	   SIP=ARI(5)
	ENDIF
        SIPAR=SIP/Z**(2.D0/3.D0)
	RETURN
	END
c  *********************************************************
      subroutine TOLLESTRUP(omega,chic2,chia2,hstp,pavg,th)
c
c  Determine polar scattering angle from Tollestrup distribution
!
!   Ref. A. Tollestrup & J. Monroe, MC note 176, September 2000
!   Derived from GEANT program  GMOLIE
c
C.    *                                                                *
C.    * THRED(NA)=reduced angles of Moliere theory                     *
C.    *                                                                *
C.    * F0I(NA),F1I(NA),F2I(NA)= integrals of Moliere functions        *
C.    *                                                                *
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
      dimension tab1el(4),tab1in(4),tab3(4)
!
      PARAMETER (ENEPER = 2.7182818d0) 
      DIMENSION TINT(40),ARG(4),VAL(4),THRED(40),F0I(40),F1I(40),F2I(40)
      DATA THRED/
     +   0.00, 0.10, 0.20, 0.30
     +,  0.40, 0.50, 0.60, 0.70
     +,  0.80, 0.90, 1.00, 1.10
     +,  1.20, 1.30, 1.40, 1.50
     +,  1.60, 1.70, 1.80, 1.90
     +,  2.00, 2.20, 2.40, 2.60
     +,  2.80, 3.00, 3.20, 3.40
     +,  3.60, 3.80, 4.00, 5.00
     +,  6.00, 7.00, 8.00, 9.00
     +, 10.00,11.00,12.00,13.00/
      DATA F0I/
     +  0.000000E+00 ,0.995016E-02 ,0.392106E-01 ,0.860688E-01
     + ,0.147856E+00 ,0.221199E+00 ,0.302324E+00 ,0.387374E+00
     + ,0.472708E+00 ,0.555142E+00 ,0.632121E+00 ,0.701803E+00
     + ,0.763072E+00 ,0.815480E+00 ,0.859142E+00 ,0.894601E+00
     + ,0.922695E+00 ,0.944424E+00 ,0.960836E+00 ,0.972948E+00
     + ,0.981684E+00 ,0.992093E+00 ,0.996849E+00 ,0.998841E+00
     + ,0.999606E+00 ,0.999877E+00 ,0.999964E+00 ,0.999990E+00
     + ,0.999998E+00 ,0.999999E+00 ,0.100000E+01 ,0.100000E+01
     + ,0.100000E+01 ,0.100000E+01 ,0.100000E+01 ,0.100000E+01
     + ,1.,1.,1.,1./
      DATA F1I/
     +  0.000000E+00,0.414985E-02,0.154894E-01,0.310312E-01
     + ,0.464438E-01,0.569008E-01,0.580763E-01,0.468264E-01
     + ,0.217924E-01,-0.163419E-01,-0.651205E-01,-0.120503E+00
     + ,-0.178272E+00,-0.233580E+00,-0.282442E+00,-0.321901E+00
     + ,-0.350115E+00,-0.366534E+00,-0.371831E+00,-0.367378E+00
     + ,-0.354994E+00,-0.314803E+00,-0.266539E+00,-0.220551E+00
     + ,-0.181546E+00,-0.150427E+00,-0.126404E+00,-0.107830E+00
     + ,-0.933106E-01,-0.817375E-01,-0.723389E-01,-0.436650E-01
     + ,-0.294700E-01,-0.212940E-01,-0.161406E-01,-0.126604E-01
     + ,-0.102042E-01,-0.840465E-02,-0.704261E-02,-0.598886E-02/
      DATA F2I/
     +  0.0,0.121500E-01,0.454999E-01,0.913000E-01
     + ,0.137300E+00,0.171400E+00,0.183900E+00,0.170300E+00
     + ,0.132200E+00,0.763000E-01,0.126500E-01,-0.473500E-01
     + ,-0.936000E-01,-0.119750E+00,-0.123450E+00,-0.106300E+00
     + ,-0.732800E-01,-0.312400E-01,0.128450E-01,0.528800E-01
     + ,0.844100E-01,0.114710E+00,0.106200E+00,0.765830E-01
     + ,0.435800E-01,0.173950E-01,0.695001E-03,-0.809500E-02
     + ,-0.117355E-01,-0.125449E-01,-0.120280E-01,-0.686530E-02
     + ,-0.385275E-02,-0.231115E-02,-0.147056E-02,-0.982480E-03
     + ,-0.682440E-03,-0.489715E-03,-0.361190E-03,-0.272582E-03/
!
! Following is LOG of elastic form factor data in Tollestrup Table 1
      data tab1el/8.3337,8.8505,8.9049,8.9418/
! Following is LOG of inelastic form factor data in Tollestrup Table 1
      data tab1in/7.5004,8.0393,7.8623,7.9883/
! Following is inelastic integral data in Tollestrup Table 3
      data tab3/3.51257,2.97368,3.15061,3.0246/
      data th0/4.84e-3/       ! maximum scattering angle for mu-e [rad]
*     ------------------------------------------------------------------
*
      IF(OMEGA.LE.ENEPER)GO TO 90
! Use Lynch-Dahl algorithm to estimate width of scattering distribution
!  Ref. NIM B58:6-10, 1991
      v = 25. * omega    ! v = 0.5*omega/(1-f) ; f=0.98
      argm = v + 1.
      if( argm .le. 0. ) then
        iflag = -26
        go to 900
      end if
      thmsw = 0.510*chic2 * ( argm/v*LOG(argm) - 1. )
      thmsw = SQRT(thmsw)
      z = mparm(1,krad)
      a = aelcomp(1,krad)
      dens = mparm(2,krad)
      bet = pavg/pmass / SQRT(1.+(pavg/pmass)**2)
      n = NINT(z)
! ---------------------------------------------------------
      if( thmsw.le.th0 .or. ABS(iptyp).eq.1 ) then   ! use Tollestrup eqs. 14-18
        pc2 = 0.1571d12*dens/a*(hstp*100.)*zsmcs/(bet*bet)  ! [eV^2]
        chicsq = pc2 / (pavg*1d9)**2
        if( n .le. 4 ) then
          rlnpa = (tab1el(n)+tab1in(n)/z)/(1.+1./z)
          cnst = LOG(pc2) - 2.*rlnpa - 0.154       ! Moliere b parameter
        else            ! use standard Moliere-Bethe method with Z(Z+1)
          cnst = LOG(chicsq/(1.167*chia2))
        end if
        if( mt .eq. 'LIH' ) then          ! use eq. 17
          pc2H = 0.1571d12*0.0708/1.008*(hstp*100.)*2./(bet*bet)  ! [eV^2]
          pc2Li = 0.1571d12*0.534/6.941*(hstp*100.)*12./(bet*bet)  ! [eV^2]
          const = LOG(pc2H) - (tab1el(1)+tab1in(1))
     &      +  LOG(pc2Li) - 2.*(tab1el(3)+tab1in(3)/3.)/1.33 - 0.154         
        end if
! .........................................................      
      else         ! heavy particle with th > th0, use Tollestrup eqs. 19-23
        zsmcs = 0.
        do i=1,nelcomp(krad)
        zsmcs = zsmcs + felcomp(i,krad) * zelcomp(i,krad)**2
        end do
        pc2 = 0.1571d12*dens/a*(hstp*100.)*zsmcs/(bet*bet)
        chicsq = pc2 / (pavg*1d9)**2
        if( n .le. 4 ) then
          bel = tab3(n) + LOG(th0*(pavg*1d9)/1d5)
          cnst = LOG(pc2) - 2.*tab1el(n) - 0.154 + 2.*bel/z
        else            ! use standard Moliere method with Z*Z
! Recompute compound sums with Z*Z
          zemcs = 0.
          zxmcs = 0.
          do i=1,nelcomp(krad)
          zemcs = zemcs + felcomp(i,krad) * 
     &     zelcomp(i,krad)**2 * 
     &     LOG( zelcomp(i,krad)**(-0.6667) ) 
          zxmcs = zxmcs + felcomp(i,krad) * 
     &     zelcomp(i,krad)**2 * 
     &     LOG( 1. + 3.34 * (zelcomp(i,krad)/137.036/bet)**2 ) 
          end do
          chia2 = 2.016e-5 * facmms /
     &     (pavg*1e3)**2 * EXP( (zxmcs-zemcs)/zsmcs )
          cnst = LOG(chicsq/(1.167*chia2))
        end if          ! end of test on n (=Z)
        if( mt .eq. 'LIH' ) then          ! use eq. 21
          pc2H = 0.1571d12*0.0708/1.008*(hstp*100.)*1./(bet*bet)  ! [eV^2]
          pc2Li = 0.1571d12*0.534/6.941*(hstp*100.)*9./(bet*bet)  ! [eV^2]
          const = LOG(pc2H)-2.*tab1el(1) + LOG(pc2Li)-2.*tab1el(3)
     &      - 0.154 + 2.*tab3(1) + 2.*tab3(3)/3. 
     &     + LOG(th0*(pavg*1d9)/1d5)         
        end if
      end if            ! end of test on thmsw
! ---------------------------------------------------------
      chic = SQRT(chicsq)
!
! Only CHIC and CONST affect the rest of this routine
      B=5.
c
      DO 10 L=1,10
         IF(ABS(B).LT.1.E-10)THEN
            B=1.E-10
         ENDIF
         DB=-(B-LOG(ABS(B))-CNST)/(1.-1./B)
         B=B+DB
         IF(ABS(DB).LE.0.0001)GO TO 20
   10 CONTINUE
c
      GO TO 90
c
   20 CONTINUE
      IF(B.LE.0.)GO TO 90
      BINV = 1./B
      TINT(1) = 0.
c
      DO 30 JA=2,4
         TINT(JA)=F0I(JA)+(F1I(JA)+F2I(JA)*BINV)*BINV
   30 CONTINUE
c
      NMAX = 4
   40 CONTINUE
      XINT = RAN1(irnarg)
c
      DO 50 NA=3,40
         IF(NA.GT.NMAX) THEN
            TINT(NA)=F0I(NA)+(F1I(NA)+F2I(NA)*BINV)*BINV
            NMAX=NA
         ENDIF
         IF(XINT.LE.TINT(NA-1)) GO TO 60
   50 CONTINUE
c
      IF(XINT.LE.TINT(40)) THEN
         NA=40
         GOTO 60
      ELSE
         TMP=1.-(1.-B*(1.-XINT))**5
         IF(TMP.LE.0.)GO TO 40
         THRI=5./TMP
         GO TO 80
      ENDIF
c
   60 CONTINUE
      NA = MAX(NA-1,3)
      NA3 = NA-3
c
      DO 70 M=1,4
         NA3M=NA3+M
         ARG(M)=TINT(NA3M)
         VAL(M)=THRED(NA3M)**2
   70 CONTINUE
c
      call POLINT(arg,val,4,xint,thri,err)    ! Num Rec
c
   80 CONTINUE
      TH=CHIC*SQRT(ABS(B*THRI))
      IF(TH.GT.PI)GO TO 40
      SINTH=SIN(TH)
      TEST=TH*(RAN1(irnarg))**2
      IF(TEST.GT.SINTH)GO TO 40
!
! Check if this angle is in the tail for thmsw < th0 case
      if( thmsw .le. th0  ) then
        if( th .gt. th0 ) then 
          test = z / (z+1.)
          if( RAN1(irnarg) .gt. test ) then
            go to 40     ! deplete tail, find some other angle for this step
          end if
        end if   
      end if
!
      GOTO 900
!
   90 CONTINUE
      iflag = -26
c
900   continue
      END
! *********************************************************
      SUBROUTINE T2ION(XMIN,XMAX,P,PM,SPIN,RES)
!
! Auxiliary routine for SAMCS
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)
      PARAMETER (PME=0.51099906D-03)
      E=SQRT(P**2+PM**2)
      ARA=2.*PME*E/(PME**2+PM**2)
      C2=ARA/(1.D0+ARA)
      BET2=(P/E)**2
      a1aa=-1.+bet2*(1.+ara)**2/ara**2
     !     -bet2**2*(2.+ara)/ara*spin
      a2new=(-1.+bet2/c2-bet2**2*spin)*(1.-c2)
      a2t=(a2new/(1.-c2*xmin)/(1.-c2*xmax)-bet2**2*spin)
     !    *(xmax-xmin)
      res=a1aa*log((1.-c2*xmax)/(1.-c2*xmin))+a2t
      return
      end
! *********************************************************
      SUBROUTINE T4ION(YMIN,YMAX,P,PM,SPIN,RES)
!
! Auxiliary routine for SAMCS
!
      IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)
      PARAMETER (PME=0.51099906D-03)
      E=SQRT(P**2+PM**2)
      R=2.*PME*E/(PME**2+PM**2)
      C2=R/(1.D0+R)
      BET2=(P/E)**2
      C3=BET2**2*C2**2*SPIN
      d1=-2.-bet2
      d2=1.+2.*bet2+c3
      d3=-bet2-2.*c3
      zmax=1.-c2*Ymax
      zmin=1.-c2*Ymin
      f1=1.+d1/c2+   d2/c2**2+   d3/c2**3+   c3/c2**4
      f2=  -d1/c2-2.*d2/c2**2-3.*d3/c2**3-4.*c3/c2**4
      f3=            d2/c2**2+3.*d3/c2**3+6.*c3/c2**4
      f4=                       -d3/c2**3-4.*c3/c2**4
      f5=                                    c3/c2**4
      s1=-f1*(-1./zmax**3+1./zmin**3)/3./c2
      s2=-f2*(-1./zmax**2+1./zmin**2)/2./c2 
      s3=-f3*(-1./zmax+1./zmin)/c2
      s4=-f4*log(zmax/zmin)/c2
      s5=-f5*(zmax-zmin)/c2 

      RES=s1+s2+s3+s4+s5
     
      RETURN
      END
c  *********************************************************
      function VAVILOV(rkappa,beta2,ran)
c
c  Determine ionization energy loss fluctuation
c   from Vavilov distribution
c
c  Derived from GEANT routine  GVAVIV
C.
C.    ******************************************************************
C.    *                                                                *
C.    *  Initializes the parameters for a Vavilov distribution         *
C.    *                                                                *
C.    *  BETA2 is the particle velocity squared                        *
C.    *  RKAPPA the usual straggling parameter                         *
C.    *                                                                *
C.    *  This routine has been extracted from a proposed CERNLIB       *
C.    *  set of routines (VAVCOE,VAVFSM).                              *
C.    *                                                                *
C.    *    Alberto Rotondi, Paolo Montagna                             *
C.    *    Fast calculation of Vavilov distribution NIM B47:215, 1990  *
C.    *                                                                *
C.    *       Author    K.S.Koelbig     ********                       *
C.    *                                                                *
C.    ******************************************************************
C.
      implicit double precision (a-h,o-z)
      DIMENSION AC(0:12),HC(0:8),H(9)
      PARAMETER
     1(BKMNX1 = 0.02, BKMNY1 = 0.05, BKMNX2 = 0.12, BKMNY2 = 0.05,
     2 BKMNX3 = 0.22, BKMNY3 = 0.05, BKMXX1 = 0.1 , BKMXY1 = 1   ,
     3 BKMXX2 = 0.2 , BKMXY2 = 1   , BKMXX3 = 0.3 , BKMXY3 = 1   )
      PARAMETER
     1(FBKX1 = 2/(BKMXX1-BKMNX1), FBKX2 = 2/(BKMXX2-BKMNX2),
     2 FBKX3 = 2/(BKMXX3-BKMNX3), FBKY1 = 2/(BKMXY1-BKMNY1),
     3 FBKY2 = 2/(BKMXY2-BKMNY2), FBKY3 = 2/(BKMXY3-BKMNY3))

      DIMENSION EDGEC(2:7),FNINV(5),DRK(5),DSIGM(5),ALFA(2:5)
      DIMENSION U1(13),U2(13),U3(13),U4(12),U5(13),U6(13),U7( 8),U8(13)
      DIMENSION V1(12),V2(12),V3(13),V4(12),V5(13),V6(13),V7(11),V8(11)
      DIMENSION W1(13),W2(11),W3(13),W4(13),W5(13),W6(13),       W8( 8)

      DATA FNINV /1, 0.5, 0.33333333, 0.25, 0.2/

      DATA (EDGEC(J),J=2,7)
     1/ 0.16666667E+0, 0.41666667E-1, 0.83333333E-2,
     2  0.13888889E-1, 0.69444444E-2, 0.77160493E-3/

      DATA (U1(K),K=1,13)
     1/ 0.25850868E+0,  0.32477982E-1, -0.59020496E-2,
     2  0.            , 0.24880692E-1,  0.47404356E-2,
     3 -0.74445130E-3,  0.73225731E-2,  0.           ,
     4  0.11668284E-2,  0.           , -0.15727318E-2,-0.11210142E-2/

      DATA (U2(K),K=1,13)
     1/ 0.43142611E+0,  0.40797543E-1, -0.91490215E-2,
     2  0.           ,  0.42127077E-1,  0.73167928E-2,
     3 -0.14026047E-2,  0.16195241E-1,  0.24714789E-2,
     4  0.20751278E-2,  0.           , -0.25141668E-2,-0.14064022E-2/

      DATA (U3(K),K=1,13)
     1/ 0.25225955E+0,  0.64820468E-1, -0.23615759E-1,
     2  0.           ,  0.23834176E-1,  0.21624675E-2,
     3 -0.26865597E-2, -0.54891384E-2,  0.39800522E-2,
     4  0.48447456E-2, -0.89439554E-2, -0.62756944E-2,-0.24655436E-2/

      DATA (U4(K),K=1,12)
     1/ 0.12593231E+1, -0.20374501E+0,  0.95055662E-1,
     2 -0.20771531E-1, -0.46865180E-1, -0.77222986E-2,
     3  0.32241039E-2,  0.89882920E-2, -0.67167236E-2,
     4 -0.13049241E-1,  0.18786468E-1,  0.14484097E-1/

      DATA (U5(K),K=1,13)
     1/-0.24864376E-1, -0.10368495E-2,  0.14330117E-2,
     2  0.20052730E-3,  0.18751903E-2,  0.12668869E-2,
     3  0.48736023E-3,  0.34850854E-2,  0.           ,
     4 -0.36597173E-3,  0.19372124E-2,  0.70761825E-3, 0.46898375E-3/

      DATA (U6(K),K=1,13)
     1/ 0.35855696E-1, -0.27542114E-1,  0.12631023E-1,
     2 -0.30188807E-2, -0.84479939E-3,  0.           ,
     3  0.45675843E-3, -0.69836141E-2,  0.39876546E-2,
     4 -0.36055679E-2,  0.           ,  0.15298434E-2, 0.19247256E-2/

      DATA (U7(K),K=1,8)
     1/ 0.10234691E+2, -0.35619655E+1,  0.69387764E+0,
     2 -0.14047599E+0, -0.19952390E+1, -0.45679694E+0,
     3  0.           ,  0.50505298E+0/

      DATA (U8(K),K=1,13)
     1/ 0.21487518E+2, -0.11825253E+2,  0.43133087E+1,
     2 -0.14500543E+1, -0.34343169E+1, -0.11063164E+1,
     3 -0.21000819E+0,  0.17891643E+1, -0.89601916E+0,
     4  0.39120793E+0,  0.73410606E+0,  0.           ,-0.32454506E+0/

      DATA (V1(K),K=1,12)
     1/ 0.27827257E+0, -0.14227603E-2,  0.24848327E-2,
     2  0.           ,  0.45091424E-1,  0.80559636E-2,
     3 -0.38974523E-2,  0.           , -0.30634124E-2,
     4  0.75633702E-3,  0.54730726E-2,  0.19792507E-2/

      DATA (V2(K),K=1,12)
     1/ 0.41421789E+0, -0.30061649E-1,  0.52249697E-2,
     2  0.           ,  0.12693873E+0,  0.22999801E-1,
     3 -0.86792801E-2,  0.31875584E-1, -0.61757928E-2,
     4  0.           ,  0.19716857E-1,  0.32596742E-2/

      DATA (V3(K),K=1,13)
     1/ 0.20191056E+0, -0.46831422E-1,  0.96777473E-2,
     2 -0.17995317E-2,  0.53921588E-1,  0.35068740E-2,
     3 -0.12621494E-1, -0.54996531E-2, -0.90029985E-2,
     4  0.34958743E-2,  0.18513506E-1,  0.68332334E-2,-0.12940502E-2/

      DATA (V4(K),K=1,12)
     1/ 0.13206081E+1,  0.10036618E+0, -0.22015201E-1,
     2  0.61667091E-2, -0.14986093E+0, -0.12720568E-1,
     3  0.24972042E-1, -0.97751962E-2,  0.26087455E-1,
     4 -0.11399062E-1, -0.48282515E-1, -0.98552378E-2/

      DATA (V5(K),K=1,13)
     1/ 0.16435243E-1,  0.36051400E-1,  0.23036520E-2,
     2 -0.61666343E-3, -0.10775802E-1,  0.51476061E-2,
     3  0.56856517E-2, -0.13438433E-1,  0.           ,
     4  0.           , -0.25421507E-2,  0.20169108E-2,-0.15144931E-2/

      DATA (V6(K),K=1,13)
     1/ 0.33432405E-1,  0.60583916E-2, -0.23381379E-2,
     2  0.83846081E-3, -0.13346861E-1, -0.17402116E-2,
     3  0.21052496E-2,  0.15528195E-2,  0.21900670E-2,
     4 -0.13202847E-2, -0.45124157E-2, -0.15629454E-2, 0.22499176E-3/

      DATA (V7(K),K=1,11)
     1/ 0.54529572E+1, -0.90906096E+0,  0.86122438E-1,
     2  0.           , -0.12218009E+1, -0.32324120E+0,
     3 -0.27373591E-1,  0.12173464E+0,  0.           ,
     4  0.           ,  0.40917471E-1/

      DATA (V8(K),K=1,11)
     1/ 0.93841352E+1, -0.16276904E+1,  0.16571423E+0,
     2  0.           , -0.18160479E+1, -0.50919193E+0,
     3 -0.51384654E-1,  0.21413992E+0,  0.           ,
     4  0.           ,  0.66596366E-1/

      DATA (W1(K),K=1,13)
     1/ 0.29712951E+0,  0.97572934E-2,  0.           ,
     2 -0.15291686E-2,  0.35707399E-1,  0.96221631E-2,
     3 -0.18402821E-2, -0.49821585E-2,  0.18831112E-2,
     4  0.43541673E-2,  0.20301312E-2, -0.18723311E-2,-0.73403108E-3/

      DATA (W2(K),K=1,11)
     1/ 0.40882635E+0,  0.14474912E-1,  0.25023704E-2,
     2 -0.37707379E-2,  0.18719727E+0,  0.56954987E-1,
     3  0.           ,  0.23020158E-1,  0.50574313E-2,
     4  0.94550140E-2,  0.19300232E-1/

      DATA (W3(K),K=1,13)
     1/ 0.16861629E+0,  0.           ,  0.36317285E-2,
     2 -0.43657818E-2,  0.30144338E-1,  0.13891826E-1,
     3 -0.58030495E-2, -0.38717547E-2,  0.85359607E-2,
     4  0.14507659E-1,  0.82387775E-2, -0.10116105E-1,-0.55135670E-2/

      DATA (W4(K),K=1,13)
     1/ 0.13493891E+1, -0.26863185E-2, -0.35216040E-2,
     2  0.24434909E-1, -0.83447911E-1, -0.48061360E-1,
     3  0.76473951E-2,  0.24494430E-1, -0.16209200E-1,
     4 -0.37768479E-1, -0.47890063E-1,  0.17778596E-1, 0.13179324E-1/

      DATA (W5(K),K=1,13)
     1/ 0.10264945E+0,  0.32738857E-1,  0.           ,
     2  0.43608779E-2, -0.43097757E-1, -0.22647176E-2,
     3  0.94531290E-2, -0.12442571E-1, -0.32283517E-2,
     4 -0.75640352E-2, -0.88293329E-2,  0.52537299E-2, 0.13340546E-2/

      DATA (W6(K),K=1,13)
     1/ 0.29568177E-1, -0.16300060E-2, -0.21119745E-3,
     2  0.23599053E-2, -0.48515387E-2, -0.40797531E-2,
     3  0.40403265E-3,  0.18200105E-2, -0.14346306E-2,
     4 -0.39165276E-2, -0.37432073E-2,  0.19950380E-2, 0.12222675E-2/

      DATA (W8(K),K=1,8)
     1/ 0.66184645E+1, -0.73866379E+0,  0.44693973E-1,
     2  0.           , -0.14540925E+1, -0.39529833E+0,
     3 -0.44293243E-1,  0.88741049E-1/

      vavilov=0
      IF(RKAPPA .LT. 0.01 .OR. RKAPPA .GT. 12) RETURN
      IF(RKAPPA .GE. 0.29) THEN
       ITYPE=1
       NPT=100
       WK=1/SQRT(RKAPPA)
       AC(0)=(-0.032227*BETA2-0.074275)*RKAPPA+
     1   (0.24533*BETA2+0.070152)*WK+(-0.55610*BETA2-3.1579)
       AC(8)=(-0.013483*BETA2-0.048801)*RKAPPA+
     1   (-1.6921*BETA2+8.3656)*WK+(-0.73275*BETA2-3.5226)
       DRK(1)=WK**2
       DSIGM(1)=SQRT(RKAPPA/(1-0.5*BETA2))
       DO 1 J = 1,4
       DRK(J+1)=DRK(1)*DRK(J)
       DSIGM(J+1)=DSIGM(1)*DSIGM(J)
    1  ALFA(J+1)=(FNINV(J)-BETA2*FNINV(J+1))*DRK(J)

       HC(0)=LOG(RKAPPA)+BETA2+0.42278434
       HC(1)=DSIGM(1)
       HC(2)=ALFA(3)*DSIGM(3)
       HC(3)=(3*ALFA(2)**2+ALFA(4))*DSIGM(4)-3
       HC(4)=(10*ALFA(2)*ALFA(3)+ALFA(5))*DSIGM(5)-10*HC(2)
       HC(5)=HC(2)**2
       HC(6)=HC(2)*HC(3)
       HC(7)=HC(2)*HC(5)
       DO 2 J = 2,7
    2  HC(J)=EDGEC(J)*HC(J)
       HC(8)=0.39894228*HC(1)
      ELSEIF(RKAPPA .GE. 0.22) THEN
       ITYPE=2
       NPT=150
       X=1+(RKAPPA-BKMXX3)*FBKX3
       Y=1+(SQRT(BETA2)-BKMXY3)*FBKY3
       XX=2*X
       YY=2*Y
       X2=XX*X-1
       X3=XX*X2-X
       Y2=YY*Y-1
       Y3=YY*Y2-Y
       XY=X*Y
       P2=X2*Y
       P3=X3*Y
       Q2=Y2*X
       Q3=Y3*X
       PQ=X2*Y2
       AC(1)=W1(1)+W1(2)*X+W1(4)*X3+W1(5)*Y+W1(6)*Y2+W1(7)*Y3+
     1  W1(8)*XY+W1(9)*P2+W1(10)*P3+W1(11)*Q2+W1(12)*Q3+W1(13)*PQ
       AC(2)=W2(1)+W2(2)*X+W2(3)*X2+W2(4)*X3+W2(5)*Y+W2(6)*Y2+
     1  W2(8)*XY+W2(9)*P2+W2(10)*P3+W2(11)*Q2
       AC(3)=W3(1)+W3(3)*X2+W3(4)*X3+W3(5)*Y+W3(6)*Y2+W3(7)*Y3+
     1  W3(8)*XY+W3(9)*P2+W3(10)*P3+W3(11)*Q2+W3(12)*Q3+W3(13)*PQ
       AC(4)=W4(1)+W4(2)*X+W4(3)*X2+W4(4)*X3+W4(5)*Y+W4(6)*Y2+W4(7)*Y3+
     1  W4(8)*XY+W4(9)*P2+W4(10)*P3+W4(11)*Q2+W4(12)*Q3+W4(13)*PQ
       AC(5)=W5(1)+W5(2)*X+W5(4)*X3+W5(5)*Y+W5(6)*Y2+W5(7)*Y3+
     1  W5(8)*XY+W5(9)*P2+W5(10)*P3+W5(11)*Q2+W5(12)*Q3+W5(13)*PQ
       AC(6)=W6(1)+W6(2)*X+W6(3)*X2+W6(4)*X3+W6(5)*Y+W6(6)*Y2+W6(7)*Y3+
     1  W6(8)*XY+W6(9)*P2+W6(10)*P3+W6(11)*Q2+W6(12)*Q3+W6(13)*PQ
       AC(8)=W8(1)+W8(2)*X+W8(3)*X2+W8(5)*Y+W8(6)*Y2+W8(7)*Y3+W8(8)*XY
       AC(0)=-3.05
      ELSEIF(RKAPPA .GE. 0.12) THEN
       ITYPE=3
       NPT=200
       X=1+(RKAPPA-BKMXX2)*FBKX2
       Y=1+(SQRT(BETA2)-BKMXY2)*FBKY2
       XX=2*X
       YY=2*Y
       X2=XX*X-1
       X3=XX*X2-X
       Y2=YY*Y-1
       Y3=YY*Y2-Y
       XY=X*Y
       P2=X2*Y
       P3=X3*Y
       Q2=Y2*X
       Q3=Y3*X
       PQ=X2*Y2
       AC(1)=V1(1)+V1(2)*X+V1(3)*X2+V1(5)*Y+V1(6)*Y2+V1(7)*Y3+
     1  V1(9)*P2+V1(10)*P3+V1(11)*Q2+V1(12)*Q3
       AC(2)=V2(1)+V2(2)*X+V2(3)*X2+V2(5)*Y+V2(6)*Y2+V2(7)*Y3+
     1  V2(8)*XY+V2(9)*P2+V2(11)*Q2+V2(12)*Q3
       AC(3)=V3(1)+V3(2)*X+V3(3)*X2+V3(4)*X3+V3(5)*Y+V3(6)*Y2+V3(7)*Y3+
     1  V3(8)*XY+V3(9)*P2+V3(10)*P3+V3(11)*Q2+V3(12)*Q3+V3(13)*PQ
       AC(4)=V4(1)+V4(2)*X+V4(3)*X2+V4(4)*X3+V4(5)*Y+V4(6)*Y2+V4(7)*Y3+
     1  V4(8)*XY+V4(9)*P2+V4(10)*P3+V4(11)*Q2+V4(12)*Q3
       AC(5)=V5(1)+V5(2)*X+V5(3)*X2+V5(4)*X3+V5(5)*Y+V5(6)*Y2+V5(7)*Y3+
     1  V5(8)*XY+V5(11)*Q2+V5(12)*Q3+V5(13)*PQ
       AC(6)=V6(1)+V6(2)*X+V6(3)*X2+V6(4)*X3+V6(5)*Y+V6(6)*Y2+V6(7)*Y3+
     1  V6(8)*XY+V6(9)*P2+V6(10)*P3+V6(11)*Q2+V6(12)*Q3+V6(13)*PQ
       AC(7)=V7(1)+V7(2)*X+V7(3)*X2+V7(5)*Y+V7(6)*Y2+V7(7)*Y3+
     1  V7(8)*XY+V7(11)*Q2
       AC(8)=V8(1)+V8(2)*X+V8(3)*X2+V8(5)*Y+V8(6)*Y2+V8(7)*Y3+
     1  V8(8)*XY+V8(11)*Q2
       AC(0)=-3.04
      ELSE
       ITYPE=4
       IF(RKAPPA .GE. 0.02) ITYPE=3
       NPT=200
       X=1+(RKAPPA-BKMXX1)*FBKX1
       Y=1+(SQRT(BETA2)-BKMXY1)*FBKY1
       XX=2*X
       YY=2*Y
       X2=XX*X-1
       X3=XX*X2-X
       Y2=YY*Y-1
       Y3=YY*Y2-Y
       XY=X*Y
       P2=X2*Y
       P3=X3*Y
       Q2=Y2*X
       Q3=Y3*X
       PQ=X2*Y2
       IF(ITYPE .EQ. 4) GO TO 4
       AC(1)=U1(1)+U1(2)*X+U1(3)*X2+U1(5)*Y+U1(6)*Y2+U1(7)*Y3+
     1  U1(8)*XY+U1(10)*P3+U1(12)*Q3+U1(13)*PQ
       AC(2)=U2(1)+U2(2)*X+U2(3)*X2+U2(5)*Y+U2(6)*Y2+U2(7)*Y3+
     1  U2(8)*XY+U2(9)*P2+U2(10)*P3+U2(12)*Q3+U2(13)*PQ
       AC(3)=U3(1)+U3(2)*X+U3(3)*X2+U3(5)*Y+U3(6)*Y2+U3(7)*Y3+
     1  U3(8)*XY+U3(9)*P2+U3(10)*P3+U3(11)*Q2+U3(12)*Q3+U3(13)*PQ
       AC(4)=U4(1)+U4(2)*X+U4(3)*X2+U4(4)*X3+U4(5)*Y+U4(6)*Y2+U4(7)*Y3+
     1  U4(8)*XY+U4(9)*P2+U4(10)*P3+U4(11)*Q2+U4(12)*Q3
       AC(5)=U5(1)+U5(2)*X+U5(3)*X2+U5(4)*X3+U5(5)*Y+U5(6)*Y2+U5(7)*Y3+
     1  U5(8)*XY+U5(10)*P3+U5(11)*Q2+U5(12)*Q3+U5(13)*PQ
       AC(6)=U6(1)+U6(2)*X+U6(3)*X2+U6(4)*X3+U6(5)*Y+U6(7)*Y3+
     1  U6(8)*XY+U6(9)*P2+U6(10)*P3+U6(12)*Q3+U6(13)*PQ
    4  AC(7)=U7(1)+U7(2)*X+U7(3)*X2+U7(4)*X3+U7(5)*Y+U7(6)*Y2+U7(8)*XY
       AC(8)=U8(1)+U8(2)*X+U8(3)*X2+U8(4)*X3+U8(5)*Y+U8(6)*Y2+U8(7)*Y3+
     1  U8(8)*XY+U8(9)*P2+U8(10)*P3+U8(11)*Q2+U8(13)*PQ
       AC(0)=-3.03
      ENDIF
      AC(9)=(AC(8)-AC(0))/NPT
      IF(ITYPE .EQ. 3) THEN
       X=(AC(7)-AC(8))/(AC(7)*AC(8))
       Y=1/LOG(AC(8)/AC(7))
       P2=AC(7)**2
       AC(11)=P2*(AC(1)*EXP(-AC(2)*(AC(7)+AC(5)*P2)-
     1            AC(3)*EXP(-AC(4)*(AC(7)+AC(6)*P2)))-0.045*Y/AC(7))/
     2            (1+X*Y*AC(7))
       AC(12)=(0.045+X*AC(11))*Y
      ENDIF
      IF(ITYPE .EQ. 4) AC(10)=0.995/vavilov_s(AC(8))

      T = 2 * RAN / AC(9)
      RLAM=AC(0)
      FL=0
      S=0
      DO 21 N = 1,NPT
      RLAM=RLAM+AC(9)
      IF(ITYPE .EQ. 1) THEN
       FN=1
       X=(RLAM+HC(0))*HC(1)
       H(1)=X
       H(2)=X**2-1
       DO 31 K = 2,8
       FN=FN+1
   31  H(K+1)=X*H(K)-FN*H(K-1)
       Y=1+HC(7)*H(9)
       DO 32 K = 2,6
   32  Y=Y+HC(K)*H(K+1)
       FU=HC(8)*EXP(-0.5*X**2)*MAX(Y,0.d0)
      ELSEIF(ITYPE .EQ. 2) THEN
       X=RLAM**2
       FU=AC(1)*EXP(-AC(2)*(RLAM+AC(5)*X)-
     1    AC(3)*EXP(-AC(4)*(RLAM+AC(6)*X)))
      ELSEIF(ITYPE .EQ. 3) THEN
       IF(RLAM .LT. AC(7)) THEN
        X=RLAM**2
        FU=AC(1)*EXP(-AC(2)*(RLAM+AC(5)*X)-
     1     AC(3)*EXP(-AC(4)*(RLAM+AC(6)*X)))
       ELSE
        X=1/RLAM
        FU=(AC(11)*X+AC(12))*X
       ENDIF
      ELSE
       FU=AC(10)*vavilov_e(RLAM)
      ENDIF
      S=S+FL+FU
      IF(S .GT. T) GO TO 22
   21 FL=FU

   22 S0=S-FL-FU
      vavilov=RLAM-AC(9)
      IF(S .GT. S0) vavilov=vavilov+AC(9)*(T-S0)/(S-S0)
      RETURN
      END

c  *********************************************************
       function VAVILOV_E(x)
c
c  Utility routine for calculating Vavilov distribution
c
c  Derived from GEANT routine  GLANDE
c
C.    ******************************************************************
C.    *                                                                *
C.    *  Copy of the CERN library routine DENLAN (G110)                *
C.    *                                                                *
C.    ******************************************************************
C.
      implicit double precision (a-h,o-z)
      DIMENSION P1(0:4),P2(0:4),P3(0:4),P4(0:4),P5(0:4),P6(0:4)
      DIMENSION Q1(0:4),Q2(0:4),Q3(0:4),Q4(0:4),Q5(0:4),Q6(0:4)
      DIMENSION A1(1:3),A2(1:2)
C
      DATA (P1(I),I=0,4),(Q1(J),J=0,4)
     1/ 0.42598 94875E+0,-0.12497 62550E+0, 0.39842 43700E-1,
     2 -0.62982 87635E-2, 0.15111 62253E-2,
     3  1.0             ,-0.33882 60629E+0, 0.95943 93323E-1,
     4 -0.16080 42283E-1, 0.37789 42063E-2/
C
      DATA (P2(I),I=0,4),(Q2(J),J=0,4)
     1/ 0.17885 41609E+0, 0.11739 57403E+0, 0.14888 50518E-1,
     2 -0.13949 89411E-2, 0.12836 17211E-3,
     3  1.0             , 0.74287 95082E+0, 0.31539 32961E+0,
     4  0.66942 19548E-1, 0.87906 09714E-2/
C
      DATA (P3(I),I=0,4),(Q3(J),J=0,4)
     1/ 0.17885 44503E+0, 0.93591 61662E-1, 0.63253 87654E-2,
     2  0.66116 67319E-4,-0.20310 49101E-5,
     3  1.0             , 0.60978 09921E+0, 0.25606 16665E+0,
     4  0.47467 22384E-1, 0.69573 01675E-2/
C
      DATA (P4(I),I=0,4),(Q4(J),J=0,4)
     1/ 0.98740 54407E+0, 0.11867 23273E+3, 0.84927 94360E+3,
     2 -0.74377 92444E+3, 0.42702 62186E+3,
     3  1.0             , 0.10686 15961E+3, 0.33764 96214E+3,
     4  0.20167 12389E+4, 0.15970 63511E+4/
C
      DATA (P5(I),I=0,4),(Q5(J),J=0,4)
     1/ 0.10036 75074E+1, 0.16757 02434E+3, 0.47897 11289E+4,
     2  0.21217 86767E+5,-0.22324 94910E+5,
     3  1.0             , 0.15694 24537E+3, 0.37453 10488E+4,
     4  0.98346 98876E+4, 0.66924 28357E+5/
C
      DATA (P6(I),I=0,4),(Q6(J),J=0,4)
     1/ 0.10008 27619E+1, 0.66491 43136E+3, 0.62972 92665E+5,
     2  0.47555 46998E+6,-0.57436 09109E+7,
     3  1.0             , 0.65141 01098E+3, 0.56974 73333E+5,
     4  0.16591 74725E+6,-0.28157 59939E+7/
C
      DATA (A1(I),I=1,3)
     1/ 0.41666 66667E-1,-0.19965 27778E-1, 0.27095 38966E-1/
C
      DATA (A2(I),I=1,2)
     1/-0.18455 68670E+1,-0.42846 40743E+1/
C
      V=X
      IF(V .LT. -5.5) THEN
       U=EXP(V+1.0)
       vavilov_e=0.3989422803*(EXP(-1.0/U)/SQRT(U))*
     1        (1.0+(A1(1)+(A1(2)+A1(3)*U)*U)*U)
      ELSE IF(V .LT. -1.0) THEN
       U=EXP(-V-1.0)
       vavilov_e=EXP(-U)*SQRT(U)*
     1        (P1(0)+(P1(1)+(P1(2)+(P1(3)+P1(4)*V)*V)*V)*V)/
     2        (Q1(0)+(Q1(1)+(Q1(2)+(Q1(3)+Q1(4)*V)*V)*V)*V)
      ELSE IF(V .LT. 1.0) THEN
       vavilov_e=(P2(0)+(P2(1)+(P2(2)+(P2(3)+P2(4)*V)*V)*V)*V)/
     1        (Q2(0)+(Q2(1)+(Q2(2)+(Q2(3)+Q2(4)*V)*V)*V)*V)
      ELSE IF(V .LT. 5.0) THEN
       vavilov_e=(P3(0)+(P3(1)+(P3(2)+(P3(3)+P3(4)*V)*V)*V)*V)/
     1        (Q3(0)+(Q3(1)+(Q3(2)+(Q3(3)+Q3(4)*V)*V)*V)*V)
      ELSE IF(V .LT. 12.0) THEN
       U=1.0/V
       vavilov_e=U**2*(P4(0)+(P4(1)+(P4(2)+(P4(3)+P4(4)*U)*U)*U)*U)/
     1             (Q4(0)+(Q4(1)+(Q4(2)+(Q4(3)+Q4(4)*U)*U)*U)*U)
      ELSE IF(V .LT. 50.0) THEN
       U=1.0/V
       vavilov_e=U**2*(P5(0)+(P5(1)+(P5(2)+(P5(3)+P5(4)*U)*U)*U)*U)/
     1             (Q5(0)+(Q5(1)+(Q5(2)+(Q5(3)+Q5(4)*U)*U)*U)*U)
      ELSE IF(V .LT. 300.0) THEN
       U=1.0/V
       vavilov_e=U**2*(P6(0)+(P6(1)+(P6(2)+(P6(3)+P6(4)*U)*U)*U)*U)/
     1             (Q6(0)+(Q6(1)+(Q6(2)+(Q6(3)+Q6(4)*U)*U)*U)*U)
      ELSE
       U=1.0/(V-V*LOG(V)/(V+1.0))
       vavilov_e=U**2*(1.0+(A2(1)+A2(2)*U)*U)
      END IF
      RETURN
      END

c  *********************************************************
       function VAVILOV_S(x)
c
c  Utility routine for calculating Vavilov distribution
c
c  Derived from GEANT routine  GLANDS
c
C.    ******************************************************************
C.    *                                                                *
C.    *  Copy of the CERN library routine DSTLAN (G110)                *
C.    *                                                                *
C.    ******************************************************************
C.
      implicit double precision (a-h,o-z)
      DIMENSION P1(0:4),P2(0:3),P3(0:3),P4(0:3),P5(0:3),P6(0:3)
      DIMENSION Q1(0:4),Q2(0:3),Q3(0:3),Q4(0:3),Q5(0:3),Q6(0:3)
      DIMENSION A1(1:3),A2(1:3)
C
      DATA (P1(I),I=0,4),(Q1(J),J=0,4)
     1/ 0.25140 91491E+0,-0.62505 80444E-1, 0.14583 81230E-1,
     2 -0.21088 17737E-2, 0.74112 47290E-3,
     3  1.0             ,-0.55711 75625E-2, 0.62253 10236E-1,
     4 -0.31373 78427E-2, 0.19314 96439E-2/
C
      DATA (P2(I),I=0,3),(Q2(J),J=0,3)
     1/ 0.28683 28584E+0, 0.35643 63231E+0, 0.15235 18695E+0,
     2  0.22513 04883E-1,
     3  1.0             , 0.61911 36137E+0, 0.17207 21448E+0,
     4  0.22785 94771E-1/
C
      DATA (P3(I),I=0,3),(Q3(J),J=0,3)
     1/ 0.28683 29066E+0, 0.30038 28436E+0, 0.99509 51941E-1,
     2  0.87338 27185E-2,
     3  1.0             , 0.42371 90502E+0, 0.10956 31512E+0,
     4  0.86938 51567E-2/
C
      DATA (P4(I),I=0,3),(Q4(J),J=0,3)
     1/ 0.10003 51630E+1, 0.45035 92498E+1, 0.10858 83880E+2,
     2  0.75360 52269E+1,
     3  1.0             , 0.55399 69678E+1, 0.19335 81111E+2,
     4  0.27213 21508E+2/
C
      DATA (P5(I),I=0,3),(Q5(J),J=0,3)
     1/ 0.10000 06517E+1, 0.49094 14111E+2, 0.85055 44753E+2,
     2  0.15321 53455E+3,
     3  1.0             , 0.50099 28881E+2, 0.13998 19104E+3,
     4  0.42000 02909E+3/
C
      DATA (P6(I),I=0,3),(Q6(J),J=0,3)
     1/ 0.10000 00983E+1, 0.13298 68456E+3, 0.91621 49244E+3,
     2 -0.96050 54274E+3,
     3  1.0             , 0.13398 87843E+3, 0.10559 90413E+4,
     4  0.55322 24619E+3/
C
      DATA (A1(I),I=1,3)
     1/-0.45833 33333E+0, 0.66753 47222E+0,-0.16417 41416E+1/
C
      DATA (A2(I),I=1,3)
     1/ 1.0             ,-0.42278 43351E+0,-0.20434 03138E+1/
C
      V=X
      IF(V .LT. -5.5) THEN
       U=EXP(V+1.0)
       vavilov_s=0.3989422803*EXP(-1.0/U)*SQRT(U)*
     1        (1.0+(A1(1)+(A1(2)+A1(3)*U)*U)*U)
      ELSE IF(V .LT. -1.0) THEN
       U=EXP(-V-1.0)
       vavilov_s=(EXP(-U)/SQRT(U))*
     1        (P1(0)+(P1(1)+(P1(2)+(P1(3)+P1(4)*V)*V)*V)*V)/
     2        (Q1(0)+(Q1(1)+(Q1(2)+(Q1(3)+Q1(4)*V)*V)*V)*V)
      ELSE IF(V .LT. 1.0) THEN
       vavilov_s=(P2(0)+(P2(1)+(P2(2)+P2(3)*V)*V)*V)/
     1        (Q2(0)+(Q2(1)+(Q2(2)+Q2(3)*V)*V)*V)
      ELSE IF(V .LT. 4.0) THEN
       vavilov_s=(P3(0)+(P3(1)+(P3(2)+P3(3)*V)*V)*V)/
     1        (Q3(0)+(Q3(1)+(Q3(2)+Q3(3)*V)*V)*V)
      ELSE IF(V .LT. 12.0) THEN
       U=1.0/V
       vavilov_s=(P4(0)+(P4(1)+(P4(2)+P4(3)*U)*U)*U)/
     1        (Q4(0)+(Q4(1)+(Q4(2)+Q4(3)*U)*U)*U)
      ELSE IF(V .LT. 50.0) THEN
       U=1.0/V
       vavilov_s=(P5(0)+(P5(1)+(P5(2)+P5(3)*U)*U)*U)/
     1        (Q5(0)+(Q5(1)+(Q5(2)+Q5(3)*U)*U)*U)
      ELSE IF(V .LT. 300.0) THEN
       U=1.0/V
       vavilov_s=(P6(0)+(P6(1)+(P6(2)+P6(3)*U)*U)*U)/
     1        (Q6(0)+(Q6(1)+(Q6(2)+Q6(3)*U)*U)*U)
      ELSE
       U=1.0/(V-V*LOG(V)/(V+1.0))
       vavilov_s=1.0-(A2(1)+(A2(2)+A2(3)*U)*U)*U
      END IF
      RETURN
      END
! *********************************************************
      subroutine WANG
!
! Create a pion from WANG distribution
!
      implicit double precision (a-h,o-z)
      include 'icommon.inc'
      dimension pms(3)
!
      knt = 0
      pmax = ABS( wangpmx )
!
100   continue
      knt = knt + 1
      if( knt .gt. 1000 ) then
        write(*,*) 'WANG sampling does not converge'
        write(2,*) 'WANG sampling does not converge'
        stop
      end if
      x = RANV(0d0,1d0)
      pt = RANV(0d0,pmax)
      arg = -wangb*x**wangc - wangd*pt
      cs = wanga * pmax * x * (1d0-x) * EXP(arg)
      v = RANV(0d0,wangfmx)
      if( v .gt. cs ) go to 100
!
! Convert this proton track into a pion
!      ipnum = ipnum + 1
      if( wangpmx .gt. 0d0 ) then
        iptyp = 3
      else
        iptyp = -3
      end if
      phi = RANV(0d0,twopi)
! Pion in fixed coordinate system
      pms(1) = pt * COS(phi)
      pms(2) = pt * SIN(phi)
      pms(3) = x * pmax
! Angles of proton in LAB 
      call VMAG(pp,pm)
      cth = pp(3) / pm
      pt = SQRT(pp(1)**2 + pp(2)**2)
      sth = pt / pm
      if( pt .gt. 0. ) then
        cph = pp(1) / pt
        sph = pp(2) / pt
      else
        cph = 1.
        sph = 0.
      end if
!  Pion in the LAB frame
      pp(1) = pms(1)*cth*cph - pms(2)*sph + pms(3)*sth*cph
      pp(2) = pms(1)*cth*sph + pms(2)*cph + pms(3)*sth*sph
      pp(3) = -pms(1)*sth + pms(3)*cth
!
      end
! *********************************************************
	DOUBLE PRECISION FUNCTION ZFORM(Q,R,B)
!
! Auxiliary routine for SAMCS
!
        IMPLICIT DOUBLE PRECISION (A-H,O-Z), INTEGER (I-N)
	PARAMETER (PI=3.141592653589793227D+00)
	Z=PI*B/R
	Z2=Z*Z
	R2=3.*R**2*(1.+7.*Z2/3.)/5.
	R4=3.*R**4*(1.+6.*Z2+31.*Z2**2/3.)/7.
	Q2=Q**2
	DD=0.02*R2/R4
	IF(Q2.GT.DD) THEN
	EX=EXP(-PI*B*Q)
	EX2=EX*EX
	CT=(EX2+1.)/(1.-EX2)
	ST=1.-EX2
	SK=-COS(Q*R) + Z*SIN(Q*R)*CT
	ZFORM=6.*PI*B*EX*SK/Q/ST/R/R/(1.+Z2)
        ELSE 
        ZFORM=1.-Q2*R2/6.+Q2**2*R4/120.
        ENDIF
	RETURN
	END
