C  PROGRAM mc866.for - Throw Drell-Yan & onia ->dimu in E866 configuration
c NOTE: the apparatus configuration is specified in the file E866.DAT
C SIGWT is weighted yield for one nominal month of running 
c       i.e. 1e16 protons incident; needs file WJS.DAT
      common/fourmu/pmu1l,pmu2l,smu(2),xf,pt,phimu,costh
	common /xyposn/xsp(28),ysp(28),xsn(28),ysn(28),phigam
      common/sm12/zb(28),ze(28),ptamp(28),yl(28),yh(28),xh(28),nm(28)
	common /vari/zstart,zwall,dedx(2),ztarg,targl,targrl,dmpcut
	common /retrace/pyzrt(2),xrt(2),yrt(2),thyrt(2),dert(2)
	common /restart/pyzs(2),xst(2),yst(2),thyst(2),ztrace
	common /seed/iseed,pi,amb,tau,x1,x2
	common /cross/ambmin,ambmax,xfmin,xfmax
C  COMMON FOR HBOOK4    
C SIZE OF HMEMOR MUST BE 128*GPAGES 
      COMMON/PAWC/H(409600)    
C NTUPLE (FOR PAW)  
	dimension pmu1l(4),pmu2l(4)
      INTEGER NTUPSIZ,izan(12)
      PARAMETER (NTUPSIZ = 96)  
      CHARACTER*8 NTnam(ntupsiz)
      REAL ANTUPLE(ntupsiz)  
C   
c possible target positions, zstart = -139, -12, 84, for LH2, wheel or dump
c   note: targlength=targl in inches, targden in gm/cm^3, 
c   target radiation length = targrl in cm
c   sigwt is events per 10^16 protons per 1 GeV bin per .1 xf bin
c   dmpcut is the required separation of each track in inches from the 
c   dump point (z=84") for target positions upstream of the dump point
c   ZWALL must be set to the middle of the most dense part of the absorber
c   mass range ambmin<amb<ambmax and xf range thrown
	data ambmin,ambmax,xfmin,xfmax/1.,13.,-0.99,.99/
c set for absorber wall at z=236", ztrace must be downstream of absorber
	data ztrace/500./
cDump production	data targl,zstart,zwall,nmonte/6.,84.,236.,100000/
cLH2/LD2 target
	data targl,zstart,zwall,nmonte/20.,-139.,236.,300000/
cSolid wheel target	data targl,zstart,zwall,nmonte/3.,-12.,236.,200000/
cCarbon	data targden,targrl,dmpcut/2.1,20.,2./
cBeryllium	data targden,targrl,dmpcut/1.848,35.3,2./
cAluminum	data targden,targrl,dmpcut/2.70,8.9,2./
cTitanium 	data targden,targrl,dmpcut/4.54,3.56,2./
cCopper	data targden,targrl,dmpcut/8.96,1.43,2./
cTungsten	data targden,targrl,dmpcut/19.3,.35,2./
cLH2	
	data targden,targrl,dmpcut/0.071,865.,2./
cLD2		data targden,targrl,dmpcut/0.162,757.,2./
c z positions for loading track positions into ntuple
	data izan/2,3,5,6,7,9,11,15,25,26,27,28/
c setup multiple scattering of 4" thick absorbers;  1=copper, 2=carbon
        z=29.
        av=63.55
        rha=8.96
        tm=4.*2.54
        en=100000.
	pmass=105.66
        zee=1.
        call mcscinit(1,Z,AV,RHA,TM,EN,PMASS,ZEE)
        z=6.
        av=12.01
        rha=2.2
        call mcscinit(2,Z,AV,RHA,TM,EN,PMASS,ZEE)
	call hlimit(70000)
        data ntnam/'mass','xf','x1',
     * 'x2','pt','pxlab','pylab','pzlab','elab',
     * 'p1x','p1y','p1z','p2x','p2y','p2z','xsm0p',
     * 'xboxp','xnosep','xdmpp','x104p','x176p','x236p',
     * 'x330p','xs1p','xs2p','xs3p','xs4p','ysm0p','yboxp',
     * 'ynosep','ydmpp','y104p','y176p','y236p','y330p',
     * 'ys1p','ys2p','ys3p','ys4p','xsm0n','xboxn','xnosen',
     * 'xdmpn','x104n','x176n','x236n','x330n','xs1n','xs2n',
     * 'xs3n','xs4n','ysm0n','yboxn','ynosen','ydmpn',   
     * 'y104n','y176n','y236n','y330n','ys1n','ys2n',
     * 'ys3n','ys4n','y1p','y2p','y3p','y4p','y1n','y2n',
     * 'y3n','y4n','x4p','x4N','dedxp','dedxn','mrtr','zrtr',
     * 'ztarg','yrtp','yrtn','sigwt','costh','phimu','phigam',
     * 'xrtp','xrtn','pyzrtp','pyzrtn','thyrtp','thyrtn',
     * 'pyzsp','pyzsn','thystp','thystn','dertp','dertn'/
       	lrecl=1024
	npr=1000
	call hropen(2,'aptuple','accept.nt','N',lrecl,istat)
	call hbookn(1,'sta123',ntupsiz,'aptuple',npr,ntnam)
C   
	ntry=nmonte
	xtry=nmonte
c factor for weight if 1e16 proptons incident = 1.e16*6e23*1.e-27 = 6e12
	tgm=targl*2.54*targden
	if(zstart.gt.80.)tgm=100.
	sigfac=6.e12*tgm*(ambmax-ambmin)*(xfmax-xfmin)/xtry
c use trick of choosing muon sign if target upstream of dump (see MUPAIR)
	if(zstart.lt.80.)ntry=nmonte/2
	n1=0
	n2=0
	n12=0
c acceptance for low mass = 1000., medium mass = 2000., high mass = 2750.
	amp=2000.
	call sm12fld(amp)
	do 5555 itry=1,ntry
c	type *,' itry',itry
	call mupair(good)
c	type *,' pmu1l',pmu1l,good
	if(good.lt.0.5)go to 5555
	call swim(1,iflag1)
c	type *,' flag1',iflag1
c iflag is the number of the zslice where the track failed
	if(iflag1.lt.28)go to 5555
	n1=n1+1
c	type *,' pmu2l',pmu2l
	call swim(2,iflag2)
c	type *,' flag2',iflag2
	if(iflag2.lt.28)go to 5555
	n2=n2+1
	 n12=n12+1
C   Good Dimuon Pair!  --  FILL NTUPLE WITH PHYSICS VARIABLES
         ANTUPLE(1)=amb   
         ANTUPLE(2)=xf
         ANTUPLE(3)=x1
         ANTUPLE(4)=x2
         ANTUPLE(5)=pt
         ANTUPLE(6)=pmu1l(1)+pmu2l(1)   
         ANTUPLE(7)=pmu1l(2)+pmu2l(2)
         ANTUPLE(8)=pmu1l(3)+pmu2l(3)
         ANTUPLE(9)=pmu1l(4)+pmu2l(4)
         ANTUPLE(10)=pmu1l(1)   
         ANTUPLE(11)=pmu1l(2)   
         ANTUPLE(12)=pmu1l(3)   
         ANTUPLE(13)=pmu2l(1)   
         ANTUPLE(14)=pmu2l(2) 
         ANTUPLE(15)=pmu2l(3)
c fill in track positions at various z slices
	 do 533 ib=1,12
	 ibz=izan(ib)
         ANTUPLE(ib+15)=xsp(ibz)   
         ANTUPLE(ib+27)=ysp(ibz)   
         ANTUPLE(ib+39)=xsn(ibz)   
         ANTUPLE(ib+51)=ysn(ibz)   
 533	 continue
c fill in hodoscope counter number
	 is1p=(ysp(25)+20.-.3)/2.5+1.
	 is2p=(ysp(26)+24.-.05)/3.+1.
	 is3p=(ysp(27)+45.5-.65)/7.+1.
	 ich1p=(ysp(25)+20.-.3)/.127+1.
	 ich2p=(ysp(26)+24.-.05)/.254+1.
	 ich3p=(ysp(27)+45.5-.65)/.508+1.
c	write(26,3842)ich1p,ich2p,ich3p
 3842	format(3i10)
	 is4p=(ysp(28)+49.-.496)/7.+1.
	 is1n=(ysn(25)+20.-.3)/2.5+1.
	 is2n=(ysn(26)+24.-.05)/3.+1.
	 is3n=(ysn(27)+45.5-.65)/7.+1.
	 is4n=(ysn(28)+49.-.496)/7.+1.
	 isx4p=(xsp(28)+60.)/8.+1.
	 isx4n=(xsn(28)+60.)/8.+1.
         ANTUPLE(64)=is1p   
         ANTUPLE(65)=is2p   
         ANTUPLE(66)=is3p   
         ANTUPLE(67)=is4p   
         ANTUPLE(68)=is1n   
         ANTUPLE(69)=is2n   
         ANTUPLE(70)=is3n   
         ANTUPLE(71)=is4n   
         ANTUPLE(72)=isx4p   
         ANTUPLE(73)=isx4n   
         ANTUPLE(74)=dedx(1)
         ANTUPLE(75)=dedx(2)
c fill yield and retrace variables
	dty=yrt(1)-yrt(2)
	zrtr=zstart-dty/(thyrt(1)-thyrt(2))   
	zbest=zstart
	ty=thyrt(1)-thyrt(2)+dty/(zwall-zbest)
	tx=(xrt(1)-xrt(2))/(zwall-zbest)
	amrtr=sqrt(pyzrt(1)*pyzrt(2)*(tx**2+ty**2))
         ANTUPLE(76)=amrtr
         ANTUPLE(77)=zrtr
         ANTUPLE(78)=ztarg
         ANTUPLE(79)=yrt(1)
	antuple(80)=yrt(2)
c   sigwt is events per 10^16 protons per 1 GeV bin per .1 xf bin
	call dysig(amb,x1,x2,sigwt)
c must unweight the gaussian thrown in xf
	unwt=sqrt(2.*pi*.27**2)*exp(+(xf-.05)**2/2./.27**2)
	antuple(81)=sigfac*unwt*sigwt/(amb**3)
        antuple(82)=costh
        antuple(83)=phimu
        antuple(84)=phigam
        antuple(85)=xrt(1)
        antuple(86)=xrt(2)
        antuple(87)=pyzrt(1)
        antuple(88)=pyzrt(2)
        antuple(89)=thyrt(1)
        antuple(90)=thyrt(2)
        antuple(91)=pyzs(1)
        antuple(92)=pyzs(2)
        antuple(93)=thyst(1)
        antuple(94)=thyst(2)
        antuple(95)=dert(1)
        antuple(96)=dert(2)
	 CALL HFN(1,ANTUPLE)   
 5555	continue
 2736	format(' n,n1,n2,n12',4i10)
	write(3,2736)nmonte,n1,n2,n12
	type *,' n,n1,n2,n12',nmonte,n1,n2,n12
c	type *,'       '
 2	continue
C-----CLOSE THE NTUPLE FILE--------------------------------------   
      CALL HROUT(0,I,' ')   
      CALL HREND('aptuple')   
	stop
	end
c
	subroutine mupair(good)
	common /seed/iseed,pi,amb,tau,x1,x2
      common/fourmu/pmu1l,pmu2l,smu(2),xf,pt,phimu,costh
	common /xyposn/xsp(28),ysp(28),xsn(28),ysn(28),phigam
	common /vari/zstart,zwall,dedx(2),ztarg,targl,targrl,dmpcut
	common /retrace/pyzrt(2),xrt(2),yrt(2),thyrt(2),dert(2)
	common /restart/pyzs(2),xst(2),yst(2),thyst(2),ztrace
	common /cross/ambmin,ambmax,xfmin,xfmax
      REAL PSIN2TH
      EXTERNAL PSIN2TH
	dimension pmu1l(4),pmu2l(4)
        logical first
	data pi,iseed,amu/3.1415927,1001001,.10566/
	data amb,pp,pm/5.273,800.0,0.93827/
        data first/.true./
c dimu production distribution a la Moreno, McGaughey, etc. et al.
        if(first)then
        CALL HBFUN1(194,'PSIN2TH',200,0.,3.14159265,PSIN2TH)
            EE=SQRT(PP**2+PM**2)
            SS=2.*PM**2+2.*PM*EE
            RS=SQRT(SS)
            BETL=PP/(EE+PM)
            GAML=(EE+PM)/RS
	first=.false.
        end if
 	good=0.
 7272	continue
c now pick the virtual gamma four momentum, add in onia as appropriate
	amb=ambmin+(ambmax-ambmin)*ran(iseed)
	if(amb.lt.2.1.and.amb.gt.1.3)amb=9.456
	if(amb.lt.1.3.and.amb.gt.1.1)amb=10.023
	if(amb.lt.1.1.and.amb.gt.1.0)amb=10.355
	if(amb.gt.12.5)amb=3.097
	tau=amb*amb/ss
c	type *,' tau,amb,s',tau,amb,ss	
c	xf=xfmin+(xfmax-xfmin)*ran(iseed)
c for DY, xf is approximately gaussian with slight offset of +.05 !!
	call gauss(xfa,xfb)
	xf=.27*xfa+.05
c	type *,' xf',xf
	x1=.5*xf+sqrt(.25*xf*xf+tau)
	if(x1.gt.0.99)return
	x2=-.5*xf+sqrt(.25*xf*xf+tau)
c	type *,' x1,2',x1,x2
	if(x2.lt.0.0005)return
constant in the Kaplan formula is 2.8 for JPsi, 3 for DY, 3.7 for Upsilon
	pt=3.*sqrt(1./((ran(iseed))**.2)-1)
        if(pt.gt.6.)go to 7272
	phigam=2.*pi*ran(iseed)
C      MUON MOMENTUM IN DIMUON C.M.S.
C      MU IS THE MUON MASS
      EMU=AMb/2.
      PMU=SQRT(EMU**2-aMU**2)
c if flat in cos(thet)                                                          
c      costh=-1.+2.*ran(iseed)
c if 1+cos(thet)**2 in Collins-Soper frame
        costh=cos(HRNDM1(194))
        PMUZ=PMU*costh
        sinth=sqrt(1.-costh*costh)
      PMUXY=PMU*sinth
C      BOOST ALONG PT
	ET=SQRT(PT**2+amb**2)
	BET=PT/ET
	GAM=ET/amb
C      THE MOM. ALONG THE DIRECTION OF PT
        pphimu=2.*pi*ran(iseed)
c phimu is the phi of the positive mu with respect to the pt of the gamma
	phimu=pphimu-phigam
	if(phimu.lt.0.0)phimu=phimu+2.*pi
      PPRL=PMUXY*COS(PHImu)
C      THE MOM. PERPENDICULAR TO THE PT
      PPRP=PMUXY*SIN(PHImu)
C      ENERGY OF MOUN+
      pmu1l(4)=GAM*(EMU+BET*PPRL)
C      ENERGY OF MUON-
      pmu2l(4)=GAM*(EMU-BET*PPRL)
C      LOG. MOM. OF MUON+                                                         
      pmu1l(3)=PMUZ
C      LOG. MOM. OF MUON-
      pmu2l(3)=-PMUZ
C      TAN. MOM. OF MUON+
      PPRL1=GAM*(PPRL+BET*EMU)
C      TAN. MOM. OF MUON-
      PPRL2=GAM*(-PPRL+BET*EMU)                                                   
      C=COS(PPHImu)
      S=SIN(PPHImu)
      pmu1l(1)=C*PPRL1-S*PPRP
      Pmu1l(2)=C*PPRP+S*PPRL1
      PMu2l(1)=C*PPRL2+S*PPRP                                                   
      PMu2l(2)=-C*PPRP+S*PPRL2
c	type *,' pmu1,2l',pmu1l,pmu2l
C      BOOST TO CM
      PL=XF*RS/2.
      ECM=SQRT(ET**2+PL**2)
      BET=PL/ECM
      GAM=ECM/ET
      CALL LORENZ(BET,GAM,pmu1l(4),pmu1l(3))
      CALL LORENZ(BET,GAM,pmu2l(4),pmu2l(3))
C      BOOST TO LAB, bomb out if muon z-momentum is below 10 GeV
      CALL LORENZ(BETL,GAML,pmu1l(4),pmu1l(3))
	if(pmu1l(3).lt.10.)return
      CALL LORENZ(BETL,GAML,pmu2l(4),pmu2l(3))                                    
	if(pmu2l(3).lt.10.)return
	smu(1)=1.
	smu(2)=-1.
c reflect in x-z plane if pmu1ly<0, factor of two trick -> ntry=nmonte/2
	if(pmu1l(2).lt.0..and.zstart.lt.80.)then
	 pmu1l(2)=-pmu1l(2)
	 pmu2l(2)=-pmu2l(2)
	end if
c	type 4179,pmu1l,pmu2l
 4179	format(' pmu1,2 in lab',4f10.2,/14x,4f10.2)
	good=1.
 999	continue
	return
	end
c
      SUBROUTINE LORENZ(B,G,E,P)
C LORENTZ TRANSFORMATION ALONG P
      ET=G*(E+B*P)
      P=G*(P+B*E)
      E=ET
      RETURN
      END
c
      FUNCTION PSIN2TH(R)                                                          
      DATA A /1.0/
c      DATA A /0.0/
      PSIN2TH=SIN(R)*(1.+A*COS(R)**2)
      RETURN
      END
c
	subroutine gauss(r1,r2)
c generates two independent gaussian variables, each with unit sigma
	common /seed/iseed,pi,amb,tau,x1,x2
 1	rr=ran(iseed)
	if(rr.le.0.)go to 1
	r=sqrt(-2.*alog(rr))
	thet=6.2831853*ran(iseed)
	r1=r*cos(thet)
	r2=r*sin(thet)
	return
	end
c
	SUBROUTINE SWIM(i,iFLAGG) !TRACE MUONS, GENERATE HITS
	common /seed/iseed,pi,amb,tau,x1,x2
      common/fourmu/pmu1l,pmu2l,smu(2),xf,pt,phimu,costh
	common /xyposn/xsp(28),ysp(28),xsn(28),ysn(28),phigam
	common /vari/zstart,zwall,dedx(2),ztarg,targl,targrl,dmpcut
	common /retrace/pyzrt(2),xrt(2),yrt(2),thyrt(2),dert(2)
	common /restart/pyzs(2),xst(2),yst(2),thyst(2),ztrace
      common/sm12/zb(28),ze(28),ptamp(28),yl(28),yh(28),xh(28),nm(28)
	common/xymax/imx,xmax(28),ymax(28),ymin(28),xtemp(28),ytemp(28)
	dimension pmu1l(4),pmu2l(4)
	data imx,xmax,ymax,ymin/0,56*0.,28*240./
Commence: ...ALL UNITS ARE INCHES!
	iap=0
	iflagg=0
	if(i.eq.1)then
	do 2394 ib=1,28
	 xsp(ib)=0.
	 ysp(ib)=0.
	 xsn(ib)=0.
 2394	 ysn(ib)=0.
c add multiple scattering inside finite length target
	 zback=targl*(ran(iseed))
	 ztarg=zstart+targl/2.-zback
	 dpms=.015*sqrt(zback*2.54/targrl)
	 call gauss(dpx,dpy)
	 thy=(pmu1l(2)+dpy*dpms)/pmu1l(3)
	 thx=(pmu1l(1)+dpx*dpms)/pmu1l(3)
	 pyz=sqrt(pmu1l(2)**2+pmu1l(3)**2)
	else
	 call gauss(dpx,dpy)
	 thy=(pmu2l(2)+dpy*dpms)/pmu2l(3)
	 thx=(pmu2l(1)+dpx*dpms)/pmu2l(3)
	 pyz=sqrt(pmu2l(2)**2+pmu2l(3)**2)
	end if
	pyzs(i)=pyz
	ptstart=pyz*sqrt(thx**2+thy**2)
	zold=ztarg
	deav=0.
	xst(i)=0.
	yst(i)=0.
	thyst(i)=thy
	xe=0.
	ye=0.
	izst=0
c - main swimming loop !!!!!!!
	do 50 ib=1,28
c	type *,' ib,zold,zb(ib)', ib,zold,zb(ib)
c skip to target point
	if(zold.gt.zb(ib))go to 50
	if(izst.eq.0)izst=ib
	dz=zb(ib)-zold
	yb=ye+thy*dz/2.
	xb=xe+thx*dz/2.
c is track inside absorber material?
	 if(abs(yb).lt.yl(ib)) then
	 if(nm(ib).eq.1)call decu24(pyz,dpyz,deav)
	 if(nm(ib).eq.2)call decar24(pyz,dpyz,deav)
	 call mcscat(nm(ib),pyz,dz,xb,yb,thx,thy)
	 pyz=pyz-dpyz*dz/24.
	 end if
	yb=yb+thy*dz/2.
	xb=xb+thx*dz/2.
c apply magnetic kick to track
	thy=thy-smu(i)*ptamp(ib)/pyz
	dz=ze(ib)-zb(ib)
	ye=yb+thy*dz/2.
        xe=xb+thx*dz/2.
c is track inside absorber material?
	 if(abs(yb).lt.yl(ib))then
	 if(nm(ib).eq.1)call decu24(pyz,dpyz,deav)
	 if(nm(ib).eq.2)call decar24(pyz,dpyz,deav)
	 call mcscat(nm(ib),pyz,dz,xb,yb,thx,thy)
	 pyz=pyz-dpyz*dz/24.
         end if	
	ye=ye+thy*dz/2.
        xe=xe+thx*dz/2.
	zold=ze(ib)
c       if(ib.lt.11)type *,' ye,ty,zold',ye,thy,zold,pyz
c -  retrace the scattered track at this point
	if(zold.ne.ztrace)go to 223
	ibr=ib
c	call gauss(ddx,ddy)
	yer=ye+.0000001*ddy
	thyr=thy-.0000001*ddy/(958.-ztrace)
	zoldr=zold
	pyzr=pyz-pyz**2*.0000001*ddy/(958.-ztrace)/.91
	pyz0=pyzr
c        type *,' yr,tyr,zr',yer,thyr,zoldr,pyzr
	do 222 itr=ibr,izst,-1
	ybr=yer+thyr*(zb(itr)-zoldr)
	if(zb(itr).lt.zwall.and.pyzr.eq.pyz0)pyzr=pyz+deav
	thyr=thyr+smu(i)*ptamp(itr)/pyzr
	if(itr.eq.izst)zoldr=zstart
	if(itr.gt.izst)zoldr=ze(itr-1)
	 if(zoldr.lt.zwall.and.pyzr.eq.pyz0)pyzr=pyz+deav
	yer=ybr+thyr*(zoldr-zb(itr))
c        type *,' yr,tyr,zr,pr',yer,thyr,zoldr,pyzr
	if(zoldr.eq.zstart)go to 2222
 222	continue
c  - now store the 'retraced' variables 
 2222	continue
c         type *,' yr,tyr,zr,pr',yer,thyr,zstart,pyzr
	 yrt(i)=yer
	 thyrt(i)=thyr
	 pyzrt(i)=pyzr
	 xrt(i)=xe-thx*(ztrace-zwall)
	 dedx(i)=pyzs(i)-pyz
	 dert(i)=deav
c 	 type *,' story,ty,p,x,tx',yer,thyr,pyzr,xrt(i),thx
 223	continue
c check dump cuts, i.e. minimum radius from dump point for upstream target
c        or minimum pt for dump production
	if(ib.ne.6)go to 3651
	 if(zstart.gt.80.)then
	  if(ptstart.lt.1.0)go to 999
	 else
	  dr=sqrt(ye*ye+xe*xe)
	  if(dr.lt.(75./pyz).or.dr.lt.dmpcut)go to 999
	 end if
 3651	continue
c make aperture cuts
	if(abs(ye).gt.yh(ib))go to 999
c Note: the yl(ib) cut is superceded by multuple scattering above
c	if(abs(ye).lt.yl(ib))go to 999
	if(abs(xe).gt.xh(ib))go to 999
c Store the traced track
	if(smu(i).gt.0.)then
	 xsp(ib)=xe
	 ysp(ib)=ye
	else
	 xsn(ib)=xe
	 ysn(ib)=ye
	end if
	iflagg=ib
 50	continue
 999	continue
	return
	END
c
	subroutine sm12fld(amp)
	common /seed/iseed,pi,amb,tau,x1,x2
      common/sm12/zb(28),ze(28),ptamp(28),yl(28),yh(28),xh(28),nm(28)
	dimension pt20(28),pt40(28)
c read magnetic field and aperture file
	open(1,file='e866.dat',status='old')
	do 11 i=1,28
	read(1,10)zb(i),ze(i),pt20(i),pt40(i),yl(i),yh(i),xh(i),nm(i)
 10	format(4f12.4,3f5.1,i2)
 11	continue
	close(unit=1)
	ptta=0.
	pttd=0.
	za=0.
	zd=0.
c convert the 2000 amp magnet map to requested SM12 current
	do 111 i=1,28
	if(i.gt.2.and.i.lt.25)ptamp(i)=amp*pt20(i)/2000.
	if(i.le.2.or.i.ge.25)ptamp(i)=pt20(i)
	ptta=ptta+ptamp(i)
	za=za+zb(i)*ptamp(i)
	if(zb(i).gt.84.)pttd=pttd+ptamp(i)
	if(zb(i).gt.84.)zd=zd+zb(i)*ptamp(i)
 111	continue
	za=za/ptta
	zd=zd/pttd
	do 159 i=1,28
	write(3,13),zb(i),ze(i),ptamp(i),yl(i),yh(i),xh(i),nm(i)
	type 13,zb(i),ze(i),ptamp(i),yl(i),yh(i),xh(i),nm(i)
 13	format(3f12.3,3f6.1,i2)
 159	continue
	write(3,12)amp,ptta,za,pttd,zd
c	type 12,amp,ptta,za,pttd,zd
 12	format(' amps, ptkick,zbend,ptdump,zd',f8.0,4f8.3)
	first=.false.
	return
	end
c
       subroutine DECU24(EMU,demu,deav)  
C      ENERGY LOSS IN 24" COPPER WALL.  FROM DETAILED CONVOLUTION OF ELEMENTARY 
C      CROSS SECTIONS USING VOSS.  SEE CONVLT.0XX ON FNAL07::   
C      THESE NUMBERS TAKEN FROM 
	common /seed/iseed,pi,amb,tau,x1,x2
	DIMENSION DINV(100,9),EL(9)   
       LOGICAL FIRST 
       DIMENSION EMUS(9),NENR(9),ELOS(100,9),DIST(100,9),DMAX(9) 
       DIMENSION ELOS1(100),ELOS2(100),ELOS3(100),ELOS4(100),    
     1      ELOS5(100),ELOS6(100),ELOS7(100),ELOS8(100),ELOS9(100)  
      DIMENSION DIST1(100),DIST2(100),DIST3(100),DIST4(100),    
     1      DIST5(100),DIST6(100),DIST7(100),DIST8(100),DIST9(100)  
       EQUIVALENCE (ELOS1,ELOS(1,1)),(ELOS2,ELOS(1,2)),  
     1      (ELOS3,ELOS(1,3)),(ELOS4,ELOS(1,4)),(ELOS5,ELOS(1,5)),  
     2      (ELOS6,ELOS(1,6)),(ELOS7,ELOS(1,7)),(ELOS8,ELOS(1,8)),  
     3      (ELOS9,ELOS(1,9))   
       EQUIVALENCE (DIST1,DIST(1,1)),(DIST2,DIST(1,2)),  
     1      (DIST3,DIST(1,3)),(DIST4,DIST(1,4)),(DIST5,DIST(1,5)),  
     2      (DIST6,DIST(1,6)),(DIST7,DIST(1,7)),(DIST8,DIST(1,8)),  
     3      (DIST9,DIST(1,9))   
       DATA FIRST/.TRUE./    
c       DATA A1CU/.95977/,A2CU/0.00004086/
       DATA A1CU/.90654/,A2CU/.0008997/,a3cu/.0000018939/
      DATA EMUS/10.,50.,100.,150.,200.,250.,300.,400.,500./ 
      DATA NENR/30,48,54,57,58,59,61,63,65/ 
        DATA ELOS1/0.590E+03,0.644E+03,0.702E+03,0.765E+03, 
     .  0.835E+03,0.910E+03,0.993E+03,0.108E+04,0.118E+04,0.129E+04,    
     .  0.140E+04,0.153E+04,0.167E+04,0.182E+04,0.199E+04,0.217E+04,    
     .  0.236E+04,0.257E+04,0.281E+04,0.306E+04,0.334E+04,0.364E+04,    
     .  0.397E+04,0.433E+04,0.472E+04,0.515E+04,0.562E+04,0.612E+04,    
     .  0.668E+04,0.728E+04,0.794E+04,0.866E+04,0.944E+04,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
        DATA ELOS2/0.633E+03,0.691E+03,0.753E+03,0.822E+03, 
     .  0.896E+03,0.977E+03,0.107E+04,0.116E+04,0.127E+04,0.138E+04,    
     .  0.151E+04,0.164E+04,0.179E+04,0.195E+04,0.213E+04,0.232E+04,    
     .  0.253E+04,0.276E+04,0.301E+04,0.329E+04,0.358E+04,0.391E+04,    
     .  0.426E+04,0.465E+04,0.507E+04,0.553E+04,0.603E+04,0.657E+04,    
     .  0.717E+04,0.782E+04,0.852E+04,0.929E+04,0.101E+05,0.111E+05,    
     .  0.121E+05,0.131E+05,0.143E+05,0.156E+05,0.170E+05,0.186E+05,    
     .  0.203E+05,0.221E+05,0.241E+05,0.263E+05,0.287E+05,0.313E+05,    
     .  0.341E+05,0.372E+05,0.405E+05,0.442E+05,0.482E+05,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
        DATA ELOS3/0.615E+03,0.671E+03,0.732E+03,0.798E+03, 
     .  0.870E+03,0.949E+03,0.104E+04,0.113E+04,0.123E+04,0.134E+04,    
     .  0.146E+04,0.160E+04,0.174E+04,0.190E+04,0.207E+04,0.226E+04,    
     .  0.246E+04,0.268E+04,0.293E+04,0.319E+04,0.348E+04,0.380E+04,    
     .  0.414E+04,0.451E+04,0.492E+04,0.537E+04,0.586E+04,0.639E+04,    
     .  0.696E+04,0.759E+04,0.828E+04,0.903E+04,0.985E+04,0.107E+05,    
     .  0.117E+05,0.128E+05,0.139E+05,0.152E+05,0.166E+05,0.181E+05,    
     .  0.197E+05,0.215E+05,0.234E+05,0.255E+05,0.279E+05,0.304E+05,    
     .  0.331E+05,0.361E+05,0.394E+05,0.430E+05,0.468E+05,0.511E+05,    
     .  0.557E+05,0.607E+05,0.662E+05,0.722E+05,0.788E+05,0.859E+05,    
     .  0.937E+05,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
         DATA ELOS4/0.613E+03,0.668E+03,0.729E+03,0.795E+03, 
     .  0.867E+03,0.945E+03,0.103E+04,0.112E+04,0.123E+04,0.134E+04,    
     .  0.146E+04,0.159E+04,0.173E+04,0.189E+04,0.206E+04,0.225E+04,    
     .  0.245E+04,0.267E+04,0.292E+04,0.318E+04,0.347E+04,0.378E+04,    
     .  0.412E+04,0.450E+04,0.490E+04,0.535E+04,0.583E+04,0.636E+04,    
     .  0.693E+04,0.756E+04,0.825E+04,0.899E+04,0.981E+04,0.107E+05,    
     .  0.117E+05,0.127E+05,0.139E+05,0.151E+05,0.165E+05,0.180E+05,    
     .  0.196E+05,0.214E+05,0.233E+05,0.254E+05,0.277E+05,0.302E+05,    
     .  0.330E+05,0.360E+05,0.392E+05,0.428E+05,0.466E+05,0.509E+05,    
     .  0.555E+05,0.605E+05,0.660E+05,0.719E+05,0.784E+05,0.855E+05,    
     .  0.933E+05,0.102E+06,0.111E+06,0.121E+06,0.132E+06,0.144E+06,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
        DATA ELOS5/0.613E+03,0.668E+03,0.728E+03,0.794E+03, 
     .  0.866E+03,0.945E+03,0.103E+04,0.112E+04,0.123E+04,0.134E+04,    
     .  0.146E+04,0.159E+04,0.173E+04,0.189E+04,0.206E+04,0.225E+04,    
     .  0.245E+04,0.267E+04,0.291E+04,0.318E+04,0.347E+04,0.378E+04,    
     .  0.412E+04,0.449E+04,0.490E+04,0.534E+04,0.583E+04,0.636E+04,    
     .  0.693E+04,0.756E+04,0.824E+04,0.899E+04,0.980E+04,0.107E+05,    
     .  0.117E+05,0.127E+05,0.139E+05,0.151E+05,0.165E+05,0.180E+05,    
     .  0.196E+05,0.214E+05,0.233E+05,0.254E+05,0.277E+05,0.302E+05,    
     .  0.330E+05,0.360E+05,0.392E+05,0.428E+05,0.466E+05,0.508E+05,    
     .  0.554E+05,0.605E+05,0.659E+05,0.719E+05,0.784E+05,0.855E+05,    
     .  0.932E+05,0.102E+06,0.111E+06,0.121E+06,0.132E+06,0.144E+06,    
     .  0.157E+06,0.171E+06,0.186E+06,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
        DATA ELOS6/0.613E+03,0.668E+03,0.728E+03,0.794E+03, 
     .  0.866E+03,0.945E+03,0.103E+04,0.112E+04,0.123E+04,0.134E+04,    
     .  0.146E+04,0.159E+04,0.173E+04,0.189E+04,0.206E+04,0.225E+04,    
     .  0.245E+04,0.267E+04,0.291E+04,0.318E+04,0.346E+04,0.378E+04,    
     .  0.412E+04,0.449E+04,0.490E+04,0.534E+04,0.583E+04,0.635E+04,    
     .  0.693E+04,0.756E+04,0.824E+04,0.899E+04,0.980E+04,0.107E+05,    
     .  0.117E+05,0.127E+05,0.139E+05,0.151E+05,0.165E+05,0.180E+05,    
     .  0.196E+05,0.214E+05,0.233E+05,0.254E+05,0.277E+05,0.302E+05,    
     .  0.330E+05,0.359E+05,0.392E+05,0.427E+05,0.466E+05,0.508E+05,    
     .  0.554E+05,0.605E+05,0.659E+05,0.719E+05,0.784E+05,0.855E+05,    
     .  0.932E+05,0.102E+06,0.111E+06,0.121E+06,0.132E+06,0.144E+06,    
     .  0.157E+06,0.171E+06,0.186E+06,0.203E+06,0.222E+06,0.242E+06,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
         DATA ELOS7/0.612E+03,0.668E+03,0.728E+03,0.794E+03, 
     .  0.866E+03,0.945E+03,0.103E+04,0.112E+04,0.122E+04,0.134E+04,    
     .  0.146E+04,0.159E+04,0.173E+04,0.189E+04,0.206E+04,0.225E+04,    
     .  0.245E+04,0.267E+04,0.291E+04,0.318E+04,0.346E+04,0.378E+04,    
     .  0.412E+04,0.449E+04,0.490E+04,0.534E+04,0.583E+04,0.635E+04,    
     .  0.693E+04,0.756E+04,0.824E+04,0.899E+04,0.980E+04,0.107E+05,    
     .  0.117E+05,0.127E+05,0.139E+05,0.151E+05,0.165E+05,0.180E+05,    
     .  0.196E+05,0.214E+05,0.233E+05,0.254E+05,0.277E+05,0.302E+05,    
     .  0.330E+05,0.359E+05,0.392E+05,0.427E+05,0.466E+05,0.508E+05,    
     .  0.554E+05,0.605E+05,0.659E+05,0.719E+05,0.784E+05,0.855E+05,    
     .  0.932E+05,0.102E+06,0.111E+06,0.121E+06,0.132E+06,0.144E+06,    
     .  0.157E+06,0.171E+06,0.186E+06,0.203E+06,0.222E+06,0.242E+06,    
     .  0.264E+06,0.288E+06,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
        DATA ELOS8/0.612E+03,0.668E+03,0.728E+03,0.794E+03, 
     .  0.866E+03,0.944E+03,0.103E+04,0.112E+04,0.122E+04,0.134E+04,    
     .  0.146E+04,0.159E+04,0.173E+04,0.189E+04,0.206E+04,0.225E+04,    
     .  0.245E+04,0.267E+04,0.291E+04,0.318E+04,0.346E+04,0.378E+04,    
     .  0.412E+04,0.449E+04,0.490E+04,0.534E+04,0.583E+04,0.635E+04,    
     .  0.693E+04,0.756E+04,0.824E+04,0.899E+04,0.980E+04,0.107E+05,    
     .  0.117E+05,0.127E+05,0.139E+05,0.151E+05,0.165E+05,0.180E+05,    
     .  0.196E+05,0.214E+05,0.233E+05,0.254E+05,0.277E+05,0.302E+05,    
     .  0.330E+05,0.359E+05,0.392E+05,0.427E+05,0.466E+05,0.508E+05,    
     .  0.554E+05,0.604E+05,0.659E+05,0.719E+05,0.784E+05,0.855E+05,    
     .  0.932E+05,0.102E+06,0.111E+06,0.121E+06,0.132E+06,0.144E+06,    
     .  0.157E+06,0.171E+06,0.186E+06,0.203E+06,0.222E+06,0.242E+06,    
     .  0.264E+06,0.288E+06,0.314E+06,0.342E+06,0.373E+06,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
        DATA ELOS9/0.612E+03,0.668E+03,0.728E+03,0.794E+03, 
     .  0.866E+03,0.944E+03,0.103E+04,0.112E+04,0.122E+04,0.134E+04,    
     .  0.146E+04,0.159E+04,0.173E+04,0.189E+04,0.206E+04,0.225E+04,    
     .  0.245E+04,0.267E+04,0.291E+04,0.318E+04,0.346E+04,0.378E+04,    
     .  0.412E+04,0.449E+04,0.490E+04,0.534E+04,0.583E+04,0.635E+04,    
     .  0.693E+04,0.756E+04,0.824E+04,0.898E+04,0.980E+04,0.107E+05,    
     .  0.117E+05,0.127E+05,0.139E+05,0.151E+05,0.165E+05,0.180E+05,    
     .  0.196E+05,0.214E+05,0.233E+05,0.254E+05,0.277E+05,0.302E+05,    
     .  0.330E+05,0.359E+05,0.392E+05,0.427E+05,0.466E+05,0.508E+05,    
     .  0.554E+05,0.604E+05,0.659E+05,0.719E+05,0.784E+05,0.855E+05,    
     .  0.932E+05,0.102E+06,0.111E+06,0.121E+06,0.132E+06,0.144E+06,    
     .  0.157E+06,0.171E+06,0.186E+06,0.203E+06,0.222E+06,0.242E+06,    
     .  0.264E+06,0.288E+06,0.314E+06,0.342E+06,0.373E+06,0.407E+06,    
     .  0.443E+06,0.484E+06,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
         DATA DIST1/0.578E-30,0.755E-14,0.943E-07,0.104E-03, 
     .  0.182E-02,0.395E-02,0.276E-02,0.133E-02,0.646E-03,0.343E-03,    
     .  0.193E-03,0.121E-03,0.776E-04,0.485E-04,0.330E-04,0.250E-04,    
     .  0.178E-04,0.117E-04,0.807E-05,0.602E-05,0.466E-05,0.340E-05,    
     .  0.232E-05,0.155E-05,0.105E-05,0.711E-06,0.457E-06,0.240E-06,    
     .  0.849E-07,-.180E-07,-.879E-07,-.137E-06,-.181E-06,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
        DATA DIST2/0.258E-18,0.902E-09,0.117E-04,0.632E-03, 
     .  0.267E-02,0.287E-02,0.178E-02,0.100E-02,0.590E-03,0.356E-03,    
     .  0.227E-03,0.153E-03,0.992E-04,0.662E-04,0.491E-04,0.370E-04,    
     .  0.257E-04,0.180E-04,0.135E-04,0.108E-04,0.866E-05,0.660E-05,    
     .  0.494E-05,0.380E-05,0.305E-05,0.255E-05,0.208E-05,0.163E-05,    
     .  0.129E-05,0.104E-05,0.849E-06,0.708E-06,0.613E-06,0.524E-06,    
     .  0.420E-06,0.336E-06,0.273E-06,0.224E-06,0.185E-06,0.153E-06,    
     .  0.127E-06,0.105E-06,0.873E-07,0.725E-07,0.602E-07,0.499E-07,    
     .  0.414E-07,0.341E-07,0.278E-07,0.218E-07,0.138E-07,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
        DATA DIST3/0.375E-24,0.217E-11,0.764E-06,0.175E-03, 
     .  0.156E-02,0.260E-02,0.198E-02,0.125E-02,0.791E-03,0.517E-03,    
     .  0.344E-03,0.244E-03,0.169E-03,0.116E-03,0.846E-04,0.657E-04,    
     .  0.476E-04,0.334E-04,0.244E-04,0.190E-04,0.151E-04,0.116E-04,    
     .  0.862E-05,0.647E-05,0.502E-05,0.407E-05,0.332E-05,0.260E-05,    
     .  0.203E-05,0.161E-05,0.130E-05,0.106E-05,0.898E-06,0.774E-06,    
     .  0.635E-06,0.510E-06,0.416E-06,0.343E-06,0.286E-06,0.239E-06,    
     .  0.201E-06,0.170E-06,0.144E-06,0.122E-06,0.103E-06,0.877E-07,    
     .  0.746E-07,0.634E-07,0.540E-07,0.460E-07,0.392E-07,0.333E-07,    
     .  0.284E-07,0.242E-07,0.207E-07,0.176E-07,0.149E-07,0.123E-07,    
     .  0.919E-08,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
         DATA DIST4/0.316E-25,0.653E-12,0.408E-06,0.121E-03, 
     .  0.124E-02,0.226E-02,0.186E-02,0.125E-02,0.830E-03,0.570E-03,    
     .  0.396E-03,0.291E-03,0.211E-03,0.149E-03,0.112E-03,0.888E-04,    
     .  0.664E-04,0.477E-04,0.354E-04,0.276E-04,0.221E-04,0.171E-04,    
     .  0.128E-04,0.962E-05,0.741E-05,0.593E-05,0.479E-05,0.373E-05,    
     .  0.288E-05,0.226E-05,0.180E-05,0.145E-05,0.120E-05,0.102E-05,    
     .  0.831E-06,0.663E-06,0.535E-06,0.439E-06,0.363E-06,0.302E-06,    
     .  0.253E-06,0.213E-06,0.180E-06,0.152E-06,0.129E-06,0.110E-06,    
     .  0.941E-07,0.805E-07,0.689E-07,0.591E-07,0.507E-07,0.436E-07,    
     .  0.375E-07,0.322E-07,0.277E-07,0.238E-07,0.205E-07,0.177E-07,    
     .  0.152E-07,0.131E-07,0.113E-07,0.963E-08,0.795E-08,0.558E-08,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
        DATA DIST5/0.209E-25,0.497E-12,0.333E-06,0.102E-03, 
     .  0.107E-02,0.202E-02,0.172E-02,0.120E-02,0.829E-03,0.589E-03,    
     .  0.422E-03,0.319E-03,0.237E-03,0.173E-03,0.132E-03,0.107E-03,    
     .  0.818E-04,0.601E-04,0.453E-04,0.358E-04,0.290E-04,0.227E-04,    
     .  0.172E-04,0.130E-04,0.101E-04,0.806E-05,0.649E-05,0.505E-05,    
     .  0.390E-05,0.304E-05,0.240E-05,0.192E-05,0.157E-05,0.132E-05,    
     .  0.106E-05,0.840E-06,0.672E-06,0.545E-06,0.446E-06,0.368E-06,    
     .  0.306E-06,0.255E-06,0.214E-06,0.180E-06,0.152E-06,0.129E-06,    
     .  0.110E-06,0.941E-07,0.806E-07,0.692E-07,0.595E-07,0.512E-07,    
     .  0.442E-07,0.381E-07,0.329E-07,0.285E-07,0.246E-07,0.213E-07,    
     .  0.184E-07,0.160E-07,0.138E-07,0.120E-07,0.104E-07,0.903E-08,    
     .  0.777E-08,0.657E-08,0.512E-08,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
         DATA DIST6/0.175E-25,0.427E-12,0.290E-06,0.899E-04, 
     .  0.955E-03,0.184E-02,0.161E-02,0.115E-02,0.817E-03,0.595E-03,    
     .  0.436E-03,0.337E-03,0.256E-03,0.190E-03,0.148E-03,0.121E-03,    
     .  0.945E-04,0.707E-04,0.541E-04,0.433E-04,0.354E-04,0.281E-04,    
     .  0.215E-04,0.165E-04,0.128E-04,0.103E-04,0.833E-05,0.651E-05,    
     .  0.503E-05,0.392E-05,0.309E-05,0.246E-05,0.200E-05,0.166E-05,    
     .  0.133E-05,0.105E-05,0.832E-06,0.669E-06,0.543E-06,0.444E-06,    
     .  0.365E-06,0.302E-06,0.251E-06,0.210E-06,0.176E-06,0.149E-06,    
     .  0.126E-06,0.107E-06,0.914E-07,0.783E-07,0.672E-07,0.578E-07,    
     .  0.498E-07,0.430E-07,0.372E-07,0.323E-07,0.280E-07,0.243E-07,    
     .  0.211E-07,0.183E-07,0.159E-07,0.138E-07,0.120E-07,0.104E-07,    
     .  0.910E-08,0.792E-08,0.689E-08,0.594E-08,0.496E-08,0.342E-08,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
        DATA DIST7/0.150E-25,0.375E-12,0.257E-06,0.805E-04, 
     .  0.864E-03,0.169E-02,0.151E-02,0.111E-02,0.801E-03,0.594E-03,    
     .  0.443E-03,0.348E-03,0.269E-03,0.203E-03,0.160E-03,0.133E-03,    
     .  0.105E-03,0.797E-04,0.618E-04,0.500E-04,0.413E-04,0.331E-04,    
     .  0.257E-04,0.198E-04,0.156E-04,0.126E-04,0.102E-04,0.803E-05,    
     .  0.624E-05,0.488E-05,0.384E-05,0.306E-05,0.249E-05,0.205E-05,    
     .  0.164E-05,0.128E-05,0.102E-05,0.812E-06,0.654E-06,0.531E-06,    
     .  0.433E-06,0.356E-06,0.293E-06,0.243E-06,0.203E-06,0.170E-06,    
     .  0.143E-06,0.121E-06,0.102E-06,0.873E-07,0.746E-07,0.640E-07,    
     .  0.551E-07,0.475E-07,0.411E-07,0.356E-07,0.309E-07,0.268E-07,    
     .  0.233E-07,0.203E-07,0.176E-07,0.154E-07,0.134E-07,0.117E-07,    
     .  0.102E-07,0.889E-08,0.776E-08,0.678E-08,0.591E-08,0.512E-08,    
     .  0.432E-08,0.317E-08,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
        DATA DIST8/0.117E-25,0.302E-12,0.211E-06,0.670E-04, 
     .  0.730E-03,0.146E-02,0.135E-02,0.102E-02,0.764E-03,0.584E-03,     
     .  0.447E-03,0.359E-03,0.284E-03,0.219E-03,0.176E-03,0.150E-03,    
     .  0.121E-03,0.939E-04,0.743E-04,0.610E-04,0.513E-04,0.418E-04,    
     .  0.331E-04,0.260E-04,0.208E-04,0.170E-04,0.139E-04,0.111E-04,    
     .  0.872E-05,0.688E-05,0.546E-05,0.436E-05,0.355E-05,0.292E-05,    
     .  0.233E-05,0.182E-05,0.144E-05,0.114E-05,0.913E-06,0.734E-06,    
     .  0.593E-06,0.481E-06,0.392E-06,0.321E-06,0.264E-06,0.218E-06,    
     .  0.181E-06,0.151E-06,0.127E-06,0.107E-06,0.904E-07,0.768E-07,     
     .  0.656E-07,0.562E-07,0.483E-07,0.417E-07,0.360E-07,0.312E-07,    
     .  0.271E-07,0.236E-07,0.206E-07,0.180E-07,0.157E-07,0.137E-07,    
     .  0.120E-07,0.105E-07,0.917E-08,0.803E-08,0.704E-08,0.617E-08,    
     .  0.542E-08,0.475E-08,0.415E-08,0.357E-08,0.286E-08,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
        DATA DIST9/0.960E-26,0.254E-12,0.179E-06,0.576E-04, 
     .  0.635E-03,0.130E-02,0.123E-02,0.953E-03,0.728E-03,0.567E-03,    
     .  0.443E-03,0.362E-03,0.291E-03,0.229E-03,0.187E-03,0.161E-03,    
     .  0.132E-03,0.104E-03,0.836E-04,0.696E-04,0.592E-04,0.490E-04,    
     .  0.394E-04,0.314E-04,0.254E-04,0.210E-04,0.174E-04,0.140E-04,    
     .  0.112E-04,0.890E-05,0.712E-05,0.573E-05,0.468E-05,0.386E-05,    
     .  0.309E-05,0.243E-05,0.192E-05,0.152E-05,0.121E-05,0.972E-06,    
     .  0.781E-06,0.629E-06,0.509E-06,0.413E-06,0.337E-06,0.275E-06,    
     .  0.226E-06,0.187E-06,0.155E-06,0.129E-06,0.108E-06,0.910E-07,    
     .  0.769E-07,0.653E-07,0.557E-07,0.477E-07,0.410E-07,0.354E-07,    
     .  0.306E-07,0.266E-07,0.231E-07,0.202E-07,0.176E-07,0.154E-07,       
     .  0.135E-07,0.118E-07,0.103E-07,0.906E-08,0.795E-08,0.697E-08,    
     .  0.612E-08,0.539E-08,0.474E-08,0.417E-08,0.367E-08,0.321E-08,    
     .  0.273E-08,0.198E-08,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,    
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/    
         DATA DMAX/.395E-02,.287E-02,.26E-02,.226E-02,.202E-02,  
     .  .184E-02,.169E-02,.146E-02,.13E-02/ 
C      NEXT LINE is average energy loss
       deav=deav+A1CU+A2CU*EMU+a3cu*emu*emu  
C      now calculate FLUCTUATED ENERGY LOSS    
             DO 100 J=2,9    
             IF(EMUS(J).GT.EMU) GO TO 200    
  100            CONTINUE    
             IF(J.GT.9) J=9  
  200            I=J-1   
             F=(EMU-EMUS(I))/(EMUS(J)-EMUS(I))   
C   
C     TO GENERATE THE INVERSE FUNCTIONS 
c   
       IF (FIRST) THEN   
       FIRST=.FALSE. 
       DO 250 L=1,9  
       DINV(1,L)=0.  
       SUM=0.    
       DO 260 LL=2,NENR(L)   
       SUM=SUM+(DIST(LL-1,L)+DIST(LL,L))*0.5*(ELOS(LL,L)-ELOS(LL-1,L))   
       DINV(LL,L)=SUM    
 260   CONTINUE  
       DO 270 LI=2,NENR(L)   
       DINV(LI,L)=DINV(LI,L)/SUM 
 270   CONTINUE  
 250   CONTINUE  
       ENDIF 
C   
C     TO FIND THE ENERGY LOSS FOR THE TWO DISTRIBUTIONS 
C   
       RAND=Ran(iseed)   
       DO 280 IJ=I,J 
       DO 290 IK=2,NENR(IJ)  
       IF (RAND.LT.DINV(IK,IJ)) GO TO 291    
 290   CONTINUE  
 291   CONTINUE  
       I1=IK-1   
       I2=IK 
       EL(IJ)=ELOS(I1,IJ)+(ELOS(I2,IJ)-ELOS(I1,IJ))* 
     1 (RAND-DINV(I1,IJ))/(DINV(I2,IJ)-DINV(I1,IJ)) 
 280   CONTINUE  
C   
C     INTERPOLATE THE ENERGY LOSS   
C   
      ELOSS=EL(I)+F*(EL(J)-EL(I))   
      demu=ELOSS/1000.  
C   
       RETURN    
       END   
c
	subroutine dysig(zmb,z1,z2,zigwt)
c	program to calculate the cross section for the production of a Drell-
c	Yan dimuon of mass=am, xf, and pt  
c	-- uses the Structure Functions of MRS
      IMPLICIT REAL*8(A-H,O-Y)
	s=2.*.938*800.
	am=zmb
	fudge=5.
cross section M^3d2sig/dMdxf is: K*(8*pi*alpha**2/9)*(x1*x2/(x1+x2))*
c	sigma{ei**2*(q(x1)*qbar(x2)+qbar(x1)*q(x2))}   (where K=.39 GeV^2mbarn)
	fac1=.39*8*3.14159/9/137.**2
c
	ak=.39
	scale=am
	mode=1
	x1=z1
	x2=z2
	if(x1.gt..99)go to 999
	if(x2.lt.0.0005)go to 999
	fac2=x1*x2/(x1+x2)
        call MRSEB(X1,SCALE,MODE,UPV1,DNV1,SEA1,STR,CHM,BOT,GLU)
c	type *,' x1,uv,dv,s,g',x1,upv1/x1,dnv1/x1,sea1/x1,glu/x1
        call MRSEB(X2,SCALE,MODE,UPV2,DNV2,SEA2,STR,CHM,BOT,GLU)
c	type *,' x1,uv,dv,s,g',x2,upv2/x2,dnv2/x2,sea2/x2,glu/x2
	up1=(upv1+sea1)/x1
	dn1=(dnv1+sea1)/x1
	dy1=(4.*up1+dn1)/9.*sea2/x2
	up2=(upv2+sea2)/x2
	dn2=(dnv2+sea2)/x2
	dy2=(4.*up2+dn2)/9.*sea1/x1
	dy=dy1+dy2
c	write(33,22)x1,up1,dn1,sea1,dy1
c	write(33,22)x2,up2,dn2,sea2,dy2
 22	format(3f12.5)
c	type *,' u1,d1,s1,dy1',up1,dn1,sea1,dy1
c	type *,' u2,d2,s2,dy2',up2,dn2,sea2,dy2
c	format(' x,up,dn,val,dy',6f12.5)
c	type *,' ak,f1,f2,dy',ak,fac1,fac2,dy
	zigwt=fudge*ak*fac1*fac2*dy
c	write(33,36)xf,cross
 36	format(f8.3,e14.4)
c	type 35,am,xf,x1,x2,cross
 35	format(' m,xf,x1,x2',f8.1,3f9.4,e12.3)
c	type *,' cross = ',cross
	return
 999	continue
	zigwt=0.
	return
	end
c
c From:	UK%"WJS@HEP.DURHAM.AC.UK" 14-APR-1989 16:02
c To:	B92493
c Subj:	
c Received: From UKACRL(MAILER) by FNALF with RSCS id 6332
c          for B92493@FNAL; Fri, 14 Apr 89 16:02 CDT
c Received: from RL.IB by UKACRL.BITNET (Mailer X1.25) with BSMTP id 8388; Fri,
c  14 Apr 89 22:02:33 BST
c Received:
c Via:           UK.AC.DUR.HEP; 14 APR 89 22:02:28 BST
c Date:          14-APR-1989 22:01:17
c From:          WJS@HEP.DURHAM.AC.UK
c To:            B92493@FNAL
c Sender:        JANET "WJS@UK.AC.DURHAM.HEP" <WJS@HEP.DURHAM.AC.UK>
c Message-Id:    <206036A2_0006417C.00923489BAC51D80$4_1@UK.AC.DURHAM.HEP>
c Originally-to: EARN%FNAL::B92493
c Originally-from: WJS
c Mailer:        Janet_Mailshr V3.2 (27-Oct-1988)
      SUBROUTINE MRSEB(X,SCALE,MODE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
C***************************************************************C
C                                                               C
C                                                               C
C     NEW VERSIONS !!!! JANUARY 1989  (AS DESCRIBED IN          C
C     "IMPROVED PARTON DISTRIBUTIONS ... " A.D. MARTIN,         C
C     R.G. ROBERTS AND W.J. STIRLING PREPRINT RAL-88-113 )      C
C                                                               C
C  MODE 1 CORRESPONDS TO                                        C
C  MARTIN, ROBERTS, STIRLING (EMC FIT)    WITH LAMBDA= 100 MEV  C
C                                                               C
C  MODE 2  CORRESPONDS TO                                       C
C  MARTIN, ROBERTS, STIRLING (BCDMS FIT)  WITH LAMBDA= 200 MEV  C
C                                                               C
C  (  SOFT GLUE :  X G(X,Q0) = A (1-X)**4.4 )                   C
C                                                               C
C                         -*-                                   C
C                                                               C
C    (NOTE THAT X TIMES THE PARTON DISTRIBUTION FUNCTION        C
C    IS RETURNED I.E. G(X) = GLU/X ETC, AND THAT "SEA"          C
C    IS THE LIGHT QUARK SEA I.E. UBAR(X)=DBAR(X)=  ...          C
C    = SEA/X FOR A PROTON.  IF IN DOUBT, CHECK THE              C
C    MOMENTUM SUM RULE! NOTE ALSO THAT SCALE=Q IN GEV)          C
C                                                               C
C                         -*-                                   C
C                                                               C
C     (THE RANGE OF APPLICABILITY IS CURRENTLY:                 C
C     10**-4 < X < 1  AND  5 < Q**2 < 1.31 * 10**6              C
C     HIGHER Q**2 VALUES CAN BE SUPPLIED ON REQUEST             C
C     - PROBLEMS, COMMENTS ETC TO SRG$T3@GEN                    C
C                                                               C
C                                                               C
C***************************************************************C
      IMPLICIT REAL*8(A-H,O-Z)
      IF(MODE.EQ.1) CALL STRUC1(X,SCALE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
      IF(MODE.EQ.2) CALL STRUC2(X,SCALE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
      RETURN
      END
C
      SUBROUTINE STRUC1(X,SCALE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
C :::::::::::: MARTIN ROBERTS STIRLING :::::SOFT GLUE:::: 91 MEV:::
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION F(6,42,19),G(6),XX(42),N0(6)
      DATA XX/1.D-4,2.D-4,4.D-4,6.D-4,8.D-4,
     .       1.D-3,2.D-3,4.D-3,6.D-3,8.D-3,
     .       1.D-2,2.D-2,4.D-2,6.D-2,8.D-2,
     .     .1D0,.125D0,.15D0,.175D0,.2D0,.225D0,.25D0,.275D0,
     .     .3D0,.325D0,.35D0,.375D0,.4D0,.425D0,.45D0,.475D0,
     .     .5D0,.525D0,.55D0,.575D0,.6D0,.65D0,.7D0,.75D0,
     .     .8D0,.9D0,1.D0/
      DATA XMIN,XMAX,QSQMIN,QSQMAX/1.D-4,1.D0,5.D0,1310720.D0/
      DATA N0/2,5,4,5,0,0/
      DATA INIT/0/
      IF(INIT.NE.0) GOTO 10
      INIT=1
	open(11,file='wjs.dat',status='old')
      DO 20 N=1,41
      DO 20 M=1,19
      READ(11,50)F(1,N,M),F(2,N,M),F(3,N,M),F(4,N,M),F(5,N,M),F(6,N,M)
         DO 25 I=1,6
  25     F(I,N,M)=F(I,N,M)/(1.D0-XX(N))**N0(I)
  20  CONTINUE
      DO 30 J=1,15
      XX(J)=DLOG10(XX(J))+1.1D0
      DO 30 I=1,5
      DO 30 K=1,19
  30  F(I,J,K)=DLOG(F(I,J,K))*F(I,16,K)/DLOG(F(I,16,K))
  50  FORMAT(6F10.5)
      DO 40 I=1,6
      DO 40 M=1,19
  40  F(I,42,M)=0.D0
  10  CONTINUE
      IF(X.LT.XMIN) X=XMIN
      IF(X.GT.XMAX) X=XMAX
      QSQ=SCALE**2
      IF(QSQ.LT.QSQMIN) QSQ=QSQMIN
      IF(QSQ.GT.QSQMAX) QSQ=QSQMAX
      XXX=X
      IF(X.LT.1.D-1) XXX=DLOG10(X)+1.1D0
      N=0
  70  N=N+1
      IF(XXX.GT.XX(N+1)) GOTO 70
      A=(XXX-XX(N))/(XX(N+1)-XX(N))
      RM=DLOG(QSQ/QSQMIN)/DLOG(2.D0)
      B=RM-DINT(RM)
      M=1+IDINT(RM)
      DO 60 I=1,6
      G(I)= (1.D0-A)*(1.D0-B)*F(I,N,M)+(1.D0-A)*B*F(I,N,M+1)
     .    + A*(1.D0-B)*F(I,N+1,M)  + A*B*F(I,N+1,M+1)
      IF(N.GE.16) GOTO 65
      IF(I.EQ.6) GOTO 65
          FAC=(1.D0-B)*F(I,16,M)+B*F(I,16,M+1)
          G(I)=FAC**(G(I)/FAC)
  65  CONTINUE
      G(I)=G(I)*(1.D0-X)**N0(I)
  60  CONTINUE
      UPV=G(1)
      DNV=G(2)
      SEA=G(4)
      STR=G(4)
      CHM=G(5)
      GLU=G(3)
      BOT=G(6)
      RETURN
      END
      SUBROUTINE STRUC2(X,SCALE,UPV,DNV,SEA,STR,CHM,BOT,GLU)
C :::::::::::: MARTIN ROBERTS STIRLING :::::SOFT GLUE::::228 MEV:::
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION F(6,42,19),G(6),XX(42),N0(6)
      DATA XX/1.D-4,2.D-4,4.D-4,6.D-4,8.D-4,
     .       1.D-3,2.D-3,4.D-3,6.D-3,8.D-3,
     .       1.D-2,2.D-2,4.D-2,6.D-2,8.D-2,
     .     .1D0,.125D0,.15D0,.175D0,.2D0,.225D0,.25D0,.275D0,
     .     .3D0,.325D0,.35D0,.375D0,.4D0,.425D0,.45D0,.475D0,
     .     .5D0,.525D0,.55D0,.575D0,.6D0,.65D0,.7D0,.75D0,
     .     .8D0,.9D0,1.D0/
      DATA XMIN,XMAX,QSQMIN,QSQMAX/1.D-4,1.D0,5.D0,1310720.D0/
      DATA N0/2,5,4,5,0,0/
      DATA INIT/0/
      IF(INIT.NE.0) GOTO 10
      INIT=1
      DO 20 N=1,41
      DO 20 M=1,19
      READ(12,50)F(1,N,M),F(2,N,M),F(3,N,M),F(4,N,M),F(5,N,M),F(6,N,M)
         DO 25 I=1,6
  25     F(I,N,M)=F(I,N,M)/(1.D0-XX(N))**N0(I)
  20  CONTINUE
      DO 30 J=1,15
      XX(J)=DLOG10(XX(J))+1.1D0
      DO 30 I=1,5
      DO 30 K=1,19
  30  F(I,J,K)=DLOG(F(I,J,K))*F(I,16,K)/DLOG(F(I,16,K))
  50  FORMAT(6F10.5)
      DO 40 I=1,6
      DO 40 M=1,19
  40  F(I,42,M)=0.D0
  10  CONTINUE
      IF(X.LT.XMIN) X=XMIN
      IF(X.GT.XMAX) X=XMAX
      QSQ=SCALE**2
      IF(QSQ.LT.QSQMIN) QSQ=QSQMIN
      IF(QSQ.GT.QSQMAX) QSQ=QSQMAX
      XXX=X
      IF(X.LT.1.D-1) XXX=DLOG10(X)+1.1D0
      N=0
  70  N=N+1
      IF(XXX.GT.XX(N+1)) GOTO 70
      A=(XXX-XX(N))/(XX(N+1)-XX(N))
      RM=DLOG(QSQ/QSQMIN)/DLOG(2.D0)
      B=RM-DINT(RM)
      M=1+IDINT(RM)
      DO 60 I=1,6
      G(I)= (1.D0-A)*(1.D0-B)*F(I,N,M)+(1.D0-A)*B*F(I,N,M+1)
     .    + A*(1.D0-B)*F(I,N+1,M)  + A*B*F(I,N+1,M+1)
      IF(N.GE.16) GOTO 65
      IF(I.EQ.6) GOTO 65
          FAC=(1.D0-B)*F(I,16,M)+B*F(I,16,M+1)
          G(I)=FAC**(G(I)/FAC)
  65  CONTINUE
      G(I)=G(I)*(1.D0-X)**N0(I)
  60  CONTINUE
      UPV=G(1)
      DNV=G(2)
      SEA=G(4)
      STR=G(4)
      CHM=G(5)
      GLU=G(3)
      BOT=G(6)
      RETURN
      END
c
      SUBROUTINE mcscinit(mat,Z,AV,RHA,TM,EN,PMASS,ZE)
c
c Initialize fast multiple coulomb scattering routine MCSCAT
c
c 	Generates inverse probability table for Moliere scattering
c	according to specified inputs.
c
c	Assumes that particles are all ultrarelativistic, so that differences
c	in scattering among different types of particles can be neglected
c	and scattering angle distribution depends on momentum only in
c	inverse proportionality.  (I haven't checked whether this is valid
c	for electrons, but it's OK from pions to protons.)
c
c	Note that EN is reference value of momentum, which should be
c	a typical value for subsequent calls to mcscAT.
c	Each call to mcscAT scales the distribution by 1/p.
c
c	Dan Kaplan	12/27/88	adapted from CERNLIB routine MLR
C
c	Allow more than one material	DMK	12/27/92
c
c	mat = material number (1 or 2)
c
C     GENERATION OF MULTIPLE SCATTERING (MOLIERE FORMULA)
C     ***************************************************
C     Z= ATOMIC NUMBER OF SCATTERER
C     AV= ATOMIC WEIGHT OF THE SCATTERER
C     RHA= DENSITY OF THE SCATTERER
C     TM= SCATTERER THICKNESS (CM)
C     EN= INCIDENT PARTICLE MOMENTUM (MEV/C)
C     PMASS= MASS OF INCIDENT PARTICLE (MEV)
C     ZE= CHARGE OF INCIDENT PARTICLE (IN UNITS OF THE ELECTRON CHARGE)
C     FI= GENERATED SCATTERING ANGLE
C     ATA(NA),FK1(NA),FK2(NA),FK3(NA)= TABLES USED IN MOLIERE FORMULA
C     FOR MULTIPLE SCATTERING
C
C--
C--      MODIFICATIONS FOR PAM VERSION GEN 1.14 OF 841016   HRR
C--      CHANGE ALOG(ABS(BM)) AT LABEL 10 TO LOG(BM) AND ADD TEST ON
C--      BM LT 1.0 TO BRANCH TO EXIT WITH FI= 0.0 SINCE ALGORITHM
C--      BECOMES INVALID. CONVERT INTERNAL WORKINGS (BUT NOT ARGUMENTS) TO
C--      DOUBLE PRECISION AND GENERIC FUNCTIONS HENCE FORTRAN-77 ONLY.
C--
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	common /mscat/ bm(2),ylim(2),fi1(2),table(1000,2),matmin,matmax
      REAL Z, AV, RHA, TM, EN, PMASS, ZE, FI, YFL, X
      DIMENSION A(60),ARG(60),VAL(4),ATA(55),FK1(55),FK2(55),FK3(55)
      DATA ATA/0.0,0.09999990,0.19999999,0.29999983,0.39999998,
     *0.49999976,0.59999990,0.69999975,0.79999989,0.89999968,0.99999988,
     *1.19999981,1.39999962,1.59999943,1.79999924,1.99999905,2.19999981,
     *2.39999962,2.59999943,2.79999924,2.99999905,3.19999981,3.39999962,
     *3.59999943,3.79999924,3.99999905,4.19999886,4.39999962,4.59999943,
     *4.79999924,4.99999905,5.19999886,5.39999962,5.59999943,5.79999924,
     *5.99999905,6.19999886,6.39999962,6.59999943,6.79999924,6.99999905,
     *7.19999886,7.39999771,7.59999943,7.79999924,7.99999905,8.19999886,
     *8.39999771,8.59999943,8.79999924,8.99999905,9.19999886,9.39999771,
     *9.59999943,9.79999924/
      DATA FK1/0.0,0.00995016,0.03921056,0.08606875,0.14785618,
     *0.22119910,0.30232364,0.38737333,0.47270751,0.55514169,0.63212043,
     *0.76307184,0.85914123,0.92269504,0.96083593,0.98168433,0.99209297,
     *0.99684882,0.99884081,0.99960631,0.99987662,0.99996430,0.99999040,
     *0.99999768,0.99999946,0.99999988,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,
     *1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1./
      DATA FK2/0.0,0.00415007,0.01550778,0.03102755,0.04642320,
     *0.05693358,0.05814465,0.04697935,0.02206861,-0.01593698,
     *-0.06454855,-0.17762768,-0.28186089,-0.34970075,-0.37165517,
     *-0.35487890,-0.31477380,-0.26649439,-0.22037679,-0.18129718,
     *-0.15022457,-0.12617821,-0.10760975,-0.09308898,-0.08151746,
     *-0.07211959,-0.06436849,-0.05788440,-0.05234956,-0.04762324,
     *-0.04352657,-0.03994689,-0.03680661,-0.03402448,-0.03155859,
     *-0.02935604,-0.02738180,-0.02561169,-0.02401571,-0.02256661,
     *-0.02124067,-0.02002795,-0.01892212,-0.01790888,-0.01697539,
     *-0.01611052,-0.01530945,-0.01456896,-0.01388225,-0.01324315,
     *-0.01264622,-0.01208691,-0.01156165,-0.01106801,-0.01060481/
      DATA FK3/0.0,0.01223041,0.04564741,0.09133768,0.13733602,
     *0.17154741,0.18427747,0.17116070,0.13350528,0.07785457,0.01460415,
     *-0.09175974,-0.12234902,-0.07293552,0.01300466,0.08469558,
     *0.11535949,0.10710764,0.07755226,0.04448753,0.01822135,0.00146502,
     *-0.00734811,-0.01098190,-0.01178313,-0.01125726,-0.01019054,
     *-0.00896781,-0.00786496,-0.00683298,-0.00592750,-0.00514682,
     *-0.00446538,-0.00387442,-0.00336000,-0.00291120,-0.00251946,
     *-0.00217772,-0.00187836,-0.00161443,-0.00137980,-0.00117159,
     *-0.00098778,-0.00082460,-0.00067865,-0.00054704,-0.00042855,
     *-0.00032243,-0.00022696,-0.00014061,-0.00006204,0.00000983,
     *0.00007585,0.00013652,0.00019206/
	data matmin,matmax/1000000,-1000000/
c
	print *
	print *,'Initialize Silicon Scattering, material',mat,
     1	'.  Parameters:'
	print *
	print *,' Z=',z,', A=',av,', rho=',rha,' T=',tm,' cm,'
	print *,' Pinc=',en,' MeV, mass=',pmass,' MeV, z=',ze
	print *
	if (mat.lt.1.or.mat.gt.2) then
	  type *,'MCSCINIT - bad material=',mat
	  stop
	else
	  matmin=min0(matmin,mat)
	  matmax=max0(matmax,mat)
	endif
      BETA=EN/SQRT(PMASS**2+EN**2)
      IF(PMASS.LT.1.0) Q=SQRT(Z*(Z+1.))
      IF(PMASS.GE.1.0) Q=Z
      RTN=SQRT(RHA*TM/AV)
      IF(PMASS.LT.1.0) QSML=(Z+1.)*Z**(1./3.)
      IF(PMASS.GE.1.0) QSML=Z**(4./3.)
      EPSN=1.13+3.76*(Z*ZE/(137.*BETA))**2
      GAMMA=(8.831E3*QSML*RHA*TM)/(BETA**2*AV*EPSN)
      CNST= LOG(GAMMA)-0.154
      bm(mat)=5.
      L=0
   10 DBM=-(bm(mat)- LOG(bm(mat))-CNST)/(1.-1./bm(mat))
      bm(mat)=bm(mat)+DBM
      IF (bm(mat).LT.1.0) GO TO 190
      L=L+1
      IF(L-10) 20,20,190
   20 IF(ABS(DBM)-0.0001) 30,30,10
   30 fi1(mat)=3.965E-1*Q*ZE/(EN*BETA)*RTN
c
c generate inverse probability table
c
	do 1000 ifl=2,1000
	  yfl=(ifl-1)*(1.+fk2(55)/bm(mat))/999.
c
      A(4)=FK1(4)+FK2(4)/bm(mat)+FK3(4)/bm(mat)**2
      NA=0
   40 NA=NA+1
      A(NA)=FK1(NA)+FK2(NA)/bm(mat)+FK3(NA)/bm(mat)**2
      S=YFL-A(NA)
      IF(S) 60,60,50
   50 CONTINUE
      GO TO 40
   60 CONTINUE
      IF(NA.EQ.55) NA=54
      A(NA+1)=FK1(NA+1)+FK2(NA+1)/bm(mat)+FK3(NA+1)/bm(mat)**2
C
C     INTERPOLATION ACCORDING TO AITKEN"S METHOD
C     ******************************************
C     YFL= ARGUMENT VALUE
C     ARG= ARGUMENT VALUES OF THE TABLE (ARG,VAL)
C     VAL= FUNCTION VALUES OF THE TABLE (ARG,VAL)
C     NDIM= NUMBER OF POINTS IN TABLE (ARG,VAL)
C     ATI= THE RESULTING INTERPOLATED FUNCTION VALUE
C
      NDIM=4
      IF(NA.LT.3) NA=3
      DO 70 M=1,NDIM
      ARG(M)=A(NA-NDIM+1+M)
      VAL(M)=ATA(NA-NDIM+1+M)
   70 CONTINUE
      F=ATA(NA)/100.
      DELT2=0.
   80 DO 130 J=2,NDIM
      DELT1=DELT2
      IEND=J-1
      DO 90 I=1,IEND
      H=ARG(I)-ARG(J)
      IF(H) 90,140,90
   90 VAL(J)=(VAL(I)*(YFL-ARG(J))-VAL(J)*(YFL-ARG(I)))/H
      DELT2=ABS(VAL(J)-VAL(IEND))
      IF(J-2) 130,130,100
  100 IF(DELT2-F) 160,160,110
  110 IF(J-5) 130,120,120
  120 IF(DELT2-DELT1) 130,140,140
  130 CONTINUE
      GO TO 150
  140 J=IEND
      GO TO 160
  150 J=NDIM
  160 ATI=VAL(J)
  170 FI=ATI*fi1(mat)*SQRT(bm(mat))
c
	  table(ifl,mat)=fi*en
c
1000	continue
	table(1,mat)=0.
	ylim(mat)=1.+fk2(55)/bm(mat)
	fi1(mat)=fi1(mat)*en
c
	return
c
190	print *,'MCSCINIT - bad input (label 190).'
	return
	end
c
	subroutine mcscat(mat,p,dz,x,y,thx,thy)
c Fast generation of silicon-detector multiple scattering distribution
c	Dan Kaplan	12/27/88	adapted from CERNLIB routine MLR
c
c	Allow more than one material	DMK	12/27/92
c
c	Set up for dz inches of material in 4" segments - CNB 2/12/95
c
	common /mscat/ bm(2),ylim(2),fi1(2),table(1000,2),matmin,matmax
	double precision bm,ylim,fi1,table
	common /seed/iseed,pi,amb,tau,x1,x2
c
	if (mat.lt.matmin.or.mat.gt.matmax) then
	  type *,'MCSCAT - bad material=',mat
	  stop
	endif
c do scattering in hunks of 4" each in z
	dthx=0.
	dthy=0.
	izhunk=(dz+2.)/4
	dzhunk=dz/2./izhunk
	do 456 iz=1,izhunk
 4836      YFL=Ran(iseed)
      IF(YFL.Ge.(ylim(mat))) GO TO 180
c interpolate linearly within inverse probability table
	r=1000.*yfl+1.
	index=r
	if(index.gt.998)go to 4836
	if(index.gt.999.or.mat.gt.2)then
	type *,' error,mat,index',mat,index
	stop
	endif
	fi=(1.-(r-index))*table(index,mat) + 
     1		(r-index)*table(index+1,mat)
	fi=fi/(1000.*p)	!Change to p in GeV 
	go to 200
  180 ATI=SQRT(5./(1.-(1.-bm(mat)*(1.-YFL))**5))
      FI=ATI*fi1(mat)*SQRT(bm(mat))/(1000.*p)
  200   PHI=2.*PI*ran(iseed)
c add on the deflections in position
	x=x+dzhunk*dthx
	y=y+dzhunk*dthy
        dTHX=dTHX+FI*COS(PHI)
        dTHY=dTHY+FI*sin(phi)
	x=x+dzhunk*dthx
	y=y+dzhunk*dthy
 456	continue
c now return the new track angles
	thx=thx+dthx
	thy=thy+dthy
	RETURN
      END
c
	subroutine DEcar24(EMU,demu,deav)
C	Energy loss in 24" GRAPHITE wall.  From detailed convolution of elementary
C	cross sections using VOSS.  See CONVLT.0XX on FNAL07::
	common /seed/iseed,pi,amb,tau,x1,x2
	COMMON /VERS/NAMVER(10),DATVER(10)
	CHARACTER*10 NAMVER,DATVER
	DATA NAMVER(6)/'DEDXLOSB1'/
	DATA DATVER(6)/'08/03/88'/
	COMMON /LOSNAM/LOSNAME
	CHARACTER*40 LOSNAME
	DATA LOSNAME/'DEDXLOSB-dE/dx loss via VOSS'/
	COMMON /FLAGS/FLG(20)
	LOGICAL*1 FLG
	dimension dinv(100,9),el(9)
	logical first
	data first/.true./
	DATA A1C/.2600/,A2C/0.0/	!These numbers taken from
	DIMENSION EMUS(9),NENR(9),ELOS(100,9),DIST(100,9),DMAX(9)
	DIMENSION ELOS1(100),ELOS2(100),ELOS3(100),ELOS4(100),
	1	ELOS5(100),ELOS6(100),ELOS7(100),ELOS8(100),ELOS9(100)
	DIMENSION DIST1(100),DIST2(100),DIST3(100),DIST4(100),
	1	DIST5(100),DIST6(100),DIST7(100),DIST8(100),DIST9(100)
	EQUIVALENCE (ELOS1,ELOS(1,1)),(ELOS2,ELOS(1,2)),
	1	(ELOS3,ELOS(1,3)),(ELOS4,ELOS(1,4)),(ELOS5,ELOS(1,5)),
	2	(ELOS6,ELOS(1,6)),(ELOS7,ELOS(1,7)),(ELOS8,ELOS(1,8)),
	3	(ELOS9,ELOS(1,9))
	EQUIVALENCE (DIST1,DIST(1,1)),(DIST2,DIST(1,2)),
	1	(DIST3,DIST(1,3)),(DIST4,DIST(1,4)),(DIST5,DIST(1,5)),
	2	(DIST6,DIST(1,6)),(DIST7,DIST(1,7)),(DIST8,DIST(1,8)),
	3	(DIST9,DIST(1,9))
	DATA EMUS/10.,50.,100.,150.,200.,250.,300.,400.,500./
	DATA NENR/47,66,74,79,82,85,87,90,93/
        DATA ELOS1/0.172E+03,0.187E+03,0.204E+03,0.223E+03,
     .  0.243E+03,0.265E+03,0.289E+03,0.315E+03,0.344E+03,0.375E+03,
     .  0.409E+03,0.446E+03,0.486E+03,0.530E+03,0.578E+03,0.630E+03,
     .  0.687E+03,0.749E+03,0.817E+03,0.891E+03,0.972E+03,0.106E+04,
     .  0.116E+04,0.126E+04,0.137E+04,0.150E+04,0.163E+04,0.178E+04,
     .  0.194E+04,0.212E+04,0.231E+04,0.252E+04,0.275E+04,0.300E+04,
     .  0.327E+04,0.356E+04,0.389E+04,0.424E+04,0.462E+04,0.504E+04,
     .  0.550E+04,0.600E+04,0.654E+04,0.713E+04,0.778E+04,0.848E+04,
     .  0.925E+04,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA ELOS2/0.169E+03,0.185E+03,0.201E+03,0.219E+03,
     .  0.239E+03,0.261E+03,0.285E+03,0.310E+03,0.338E+03,0.369E+03,
     .  0.402E+03,0.439E+03,0.479E+03,0.522E+03,0.569E+03,0.621E+03,
     .  0.677E+03,0.738E+03,0.805E+03,0.878E+03,0.957E+03,0.104E+04,
     .  0.114E+04,0.124E+04,0.135E+04,0.148E+04,0.161E+04,0.176E+04,
     .  0.191E+04,0.209E+04,0.228E+04,0.248E+04,0.271E+04,0.295E+04,
     .  0.322E+04,0.351E+04,0.383E+04,0.418E+04,0.455E+04,0.497E+04,
     .  0.541E+04,0.590E+04,0.644E+04,0.702E+04,0.766E+04,0.835E+04,
     .  0.911E+04,0.993E+04,0.108E+05,0.118E+05,0.129E+05,0.140E+05,
     .  0.153E+05,0.167E+05,0.182E+05,0.199E+05,0.217E+05,0.236E+05,
     .  0.258E+05,0.281E+05,0.306E+05,0.334E+05,0.364E+05,0.397E+05,
     .  0.433E+05,0.472E+05,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA ELOS3/0.169E+03,0.184E+03,0.201E+03,0.219E+03,
     .  0.239E+03,0.261E+03,0.284E+03,0.310E+03,0.338E+03,0.369E+03,
     .  0.402E+03,0.439E+03,0.478E+03,0.522E+03,0.569E+03,0.620E+03,
     .  0.677E+03,0.738E+03,0.805E+03,0.877E+03,0.957E+03,0.104E+04,
     .  0.114E+04,0.124E+04,0.135E+04,0.148E+04,0.161E+04,0.175E+04,
     .  0.191E+04,0.209E+04,0.228E+04,0.248E+04,0.271E+04,0.295E+04,
     .  0.322E+04,0.351E+04,0.383E+04,0.417E+04,0.455E+04,0.496E+04,
     .  0.541E+04,0.590E+04,0.644E+04,0.702E+04,0.765E+04,0.835E+04,
     .  0.910E+04,0.993E+04,0.108E+05,0.118E+05,0.129E+05,0.140E+05,
     .  0.153E+05,0.167E+05,0.182E+05,0.199E+05,0.217E+05,0.236E+05,
     .  0.257E+05,0.281E+05,0.306E+05,0.334E+05,0.364E+05,0.397E+05,
     .  0.433E+05,0.472E+05,0.515E+05,0.562E+05,0.612E+05,0.668E+05,
     .  0.728E+05,0.794E+05,0.866E+05,0.944E+05,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA ELOS4/0.169E+03,0.184E+03,0.201E+03,0.219E+03,
     .  0.239E+03,0.261E+03,0.284E+03,0.310E+03,0.338E+03,0.369E+03,
     .  0.402E+03,0.439E+03,0.478E+03,0.522E+03,0.569E+03,0.620E+03,
     .  0.676E+03,0.738E+03,0.804E+03,0.877E+03,0.957E+03,0.104E+04,
     .  0.114E+04,0.124E+04,0.135E+04,0.148E+04,0.161E+04,0.175E+04,
     .  0.191E+04,0.209E+04,0.228E+04,0.248E+04,0.271E+04,0.295E+04,
     .  0.322E+04,0.351E+04,0.383E+04,0.417E+04,0.455E+04,0.496E+04,
     .  0.541E+04,0.590E+04,0.644E+04,0.702E+04,0.765E+04,0.835E+04,
     .  0.910E+04,0.993E+04,0.108E+05,0.118E+05,0.129E+05,0.140E+05,
     .  0.153E+05,0.167E+05,0.182E+05,0.199E+05,0.216E+05,0.236E+05,
     .  0.257E+05,0.281E+05,0.306E+05,0.334E+05,0.364E+05,0.397E+05,
     .  0.433E+05,0.472E+05,0.515E+05,0.561E+05,0.612E+05,0.668E+05,
     .  0.728E+05,0.794E+05,0.866E+05,0.944E+05,0.103E+06,0.112E+06,
     .  0.122E+06,0.134E+06,0.146E+06,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA ELOS5/0.169E+03,0.184E+03,0.201E+03,0.219E+03,
     .  0.239E+03,0.261E+03,0.284E+03,0.310E+03,0.338E+03,0.369E+03,
     .  0.402E+03,0.439E+03,0.478E+03,0.522E+03,0.569E+03,0.620E+03,
     .  0.676E+03,0.738E+03,0.804E+03,0.877E+03,0.957E+03,0.104E+04,
     .  0.114E+04,0.124E+04,0.135E+04,0.148E+04,0.161E+04,0.175E+04,
     .  0.191E+04,0.209E+04,0.228E+04,0.248E+04,0.271E+04,0.295E+04,
     .  0.322E+04,0.351E+04,0.383E+04,0.417E+04,0.455E+04,0.496E+04,
     .  0.541E+04,0.590E+04,0.644E+04,0.702E+04,0.765E+04,0.835E+04,
     .  0.910E+04,0.992E+04,0.108E+05,0.118E+05,0.129E+05,0.140E+05,
     .  0.153E+05,0.167E+05,0.182E+05,0.198E+05,0.216E+05,0.236E+05,
     .  0.257E+05,0.281E+05,0.306E+05,0.334E+05,0.364E+05,0.397E+05,
     .  0.433E+05,0.472E+05,0.515E+05,0.561E+05,0.612E+05,0.668E+05,
     .  0.728E+05,0.794E+05,0.866E+05,0.944E+05,0.103E+06,0.112E+06,
     .  0.122E+06,0.134E+06,0.146E+06,0.159E+06,0.173E+06,0.189E+06,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA ELOS6/0.169E+03,0.184E+03,0.201E+03,0.219E+03,
     .  0.239E+03,0.261E+03,0.284E+03,0.310E+03,0.338E+03,0.369E+03,
     .  0.402E+03,0.439E+03,0.478E+03,0.522E+03,0.569E+03,0.620E+03,
     .  0.676E+03,0.738E+03,0.804E+03,0.877E+03,0.957E+03,0.104E+04,
     .  0.114E+04,0.124E+04,0.135E+04,0.148E+04,0.161E+04,0.175E+04,
     .  0.191E+04,0.209E+04,0.227E+04,0.248E+04,0.271E+04,0.295E+04,
     .  0.322E+04,0.351E+04,0.383E+04,0.417E+04,0.455E+04,0.496E+04,
     .  0.541E+04,0.590E+04,0.643E+04,0.702E+04,0.765E+04,0.834E+04,
     .  0.910E+04,0.992E+04,0.108E+05,0.118E+05,0.129E+05,0.140E+05,
     .  0.153E+05,0.167E+05,0.182E+05,0.198E+05,0.216E+05,0.236E+05,
     .  0.257E+05,0.281E+05,0.306E+05,0.334E+05,0.364E+05,0.397E+05,
     .  0.433E+05,0.472E+05,0.515E+05,0.561E+05,0.612E+05,0.668E+05,
     .  0.728E+05,0.794E+05,0.866E+05,0.944E+05,0.103E+06,0.112E+06,
     .  0.122E+06,0.134E+06,0.146E+06,0.159E+06,0.173E+06,0.189E+06,
     .  0.206E+06,0.225E+06,0.245E+06,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA ELOS7/0.169E+03,0.184E+03,0.201E+03,0.219E+03,
     .  0.239E+03,0.261E+03,0.284E+03,0.310E+03,0.338E+03,0.369E+03,
     .  0.402E+03,0.439E+03,0.478E+03,0.521E+03,0.569E+03,0.620E+03,
     .  0.676E+03,0.738E+03,0.804E+03,0.877E+03,0.956E+03,0.104E+04,
     .  0.114E+04,0.124E+04,0.135E+04,0.148E+04,0.161E+04,0.175E+04,
     .  0.191E+04,0.209E+04,0.227E+04,0.248E+04,0.271E+04,0.295E+04,
     .  0.322E+04,0.351E+04,0.383E+04,0.417E+04,0.455E+04,0.496E+04,
     .  0.541E+04,0.590E+04,0.643E+04,0.702E+04,0.765E+04,0.834E+04,
     .  0.910E+04,0.992E+04,0.108E+05,0.118E+05,0.129E+05,0.140E+05,
     .  0.153E+05,0.167E+05,0.182E+05,0.198E+05,0.216E+05,0.236E+05,
     .  0.257E+05,0.281E+05,0.306E+05,0.334E+05,0.364E+05,0.397E+05,
     .  0.433E+05,0.472E+05,0.515E+05,0.561E+05,0.612E+05,0.668E+05,
     .  0.728E+05,0.794E+05,0.866E+05,0.944E+05,0.103E+06,0.112E+06,
     .  0.122E+06,0.134E+06,0.146E+06,0.159E+06,0.173E+06,0.189E+06,
     .  0.206E+06,0.225E+06,0.245E+06,0.267E+06,0.291E+06,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA ELOS8/0.169E+03,0.184E+03,0.201E+03,0.219E+03,
     .  0.239E+03,0.261E+03,0.284E+03,0.310E+03,0.338E+03,0.369E+03,
     .  0.402E+03,0.438E+03,0.478E+03,0.521E+03,0.569E+03,0.620E+03,
     .  0.676E+03,0.737E+03,0.804E+03,0.877E+03,0.956E+03,0.104E+04,
     .  0.114E+04,0.124E+04,0.135E+04,0.147E+04,0.161E+04,0.175E+04,
     .  0.191E+04,0.209E+04,0.227E+04,0.248E+04,0.270E+04,0.295E+04,
     .  0.322E+04,0.351E+04,0.383E+04,0.417E+04,0.455E+04,0.496E+04,
     .  0.541E+04,0.590E+04,0.643E+04,0.702E+04,0.765E+04,0.834E+04,
     .  0.910E+04,0.992E+04,0.108E+05,0.118E+05,0.129E+05,0.140E+05,
     .  0.153E+05,0.167E+05,0.182E+05,0.198E+05,0.216E+05,0.236E+05,
     .  0.257E+05,0.281E+05,0.306E+05,0.334E+05,0.364E+05,0.397E+05,
     .  0.433E+05,0.472E+05,0.515E+05,0.561E+05,0.612E+05,0.667E+05,
     .  0.728E+05,0.794E+05,0.866E+05,0.944E+05,0.103E+06,0.112E+06,
     .  0.122E+06,0.133E+06,0.146E+06,0.159E+06,0.173E+06,0.189E+06,
     .  0.206E+06,0.224E+06,0.245E+06,0.267E+06,0.291E+06,0.317E+06,
     .  0.346E+06,0.378E+06,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA ELOS9/0.169E+03,0.184E+03,0.201E+03,0.219E+03,
     .  0.239E+03,0.261E+03,0.284E+03,0.310E+03,0.338E+03,0.369E+03,
     .  0.402E+03,0.438E+03,0.478E+03,0.521E+03,0.569E+03,0.620E+03,
     .  0.676E+03,0.737E+03,0.804E+03,0.877E+03,0.956E+03,0.104E+04,
     .  0.114E+04,0.124E+04,0.135E+04,0.147E+04,0.161E+04,0.175E+04,
     .  0.191E+04,0.209E+04,0.227E+04,0.248E+04,0.270E+04,0.295E+04,
     .  0.322E+04,0.351E+04,0.382E+04,0.417E+04,0.455E+04,0.496E+04,
     .  0.541E+04,0.590E+04,0.643E+04,0.701E+04,0.765E+04,0.834E+04,
     .  0.910E+04,0.992E+04,0.108E+05,0.118E+05,0.129E+05,0.140E+05,
     .  0.153E+05,0.167E+05,0.182E+05,0.198E+05,0.216E+05,0.236E+05,
     .  0.257E+05,0.281E+05,0.306E+05,0.334E+05,0.364E+05,0.397E+05,
     .  0.433E+05,0.472E+05,0.515E+05,0.561E+05,0.612E+05,0.667E+05,
     .  0.728E+05,0.794E+05,0.865E+05,0.944E+05,0.103E+06,0.112E+06,
     .  0.122E+06,0.133E+06,0.146E+06,0.159E+06,0.173E+06,0.189E+06,
     .  0.206E+06,0.224E+06,0.245E+06,0.267E+06,0.291E+06,0.317E+06,
     .  0.346E+06,0.378E+06,0.412E+06,0.449E+06,0.490E+06,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA DIST1/0.282E-21,0.136E-09,0.150E-04,0.192E-02,
     .  0.116E-01,0.131E-01,0.688E-02,0.319E-02,0.163E-02,0.905E-03,
     .  0.543E-03,0.361E-03,0.231E-03,0.152E-03,0.113E-03,0.886E-04,
     .  0.616E-04,0.425E-04,0.317E-04,0.258E-04,0.210E-04,0.159E-04,
     .  0.117E-04,0.877E-05,0.688E-05,0.568E-05,0.457E-05,0.344E-05,
     .  0.260E-05,0.200E-05,0.155E-05,0.121E-05,0.976E-06,0.776E-06,
     .  0.557E-06,0.385E-06,0.266E-06,0.179E-06,0.113E-06,0.632E-07,
     .  0.256E-07,-.257E-08,-.232E-07,-.379E-07,-.480E-07,-.544E-07,
     .  -.583E-07,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA DIST2/0.106E-24,0.420E-11,0.334E-05,0.101E-02,
     .  0.932E-02,0.135E-01,0.777E-02,0.368E-02,0.188E-02,0.106E-02,
     .  0.642E-03,0.433E-03,0.286E-03,0.188E-03,0.139E-03,0.111E-03,
     .  0.803E-04,0.556E-04,0.414E-04,0.336E-04,0.280E-04,0.219E-04,
     .  0.164E-04,0.124E-04,0.985E-05,0.825E-05,0.690E-05,0.541E-05,
     .  0.422E-05,0.337E-05,0.273E-05,0.225E-05,0.192E-05,0.168E-05,
     .  0.135E-05,0.106E-05,0.847E-06,0.689E-06,0.566E-06,0.467E-06,
     .  0.386E-06,0.319E-06,0.264E-06,0.219E-06,0.181E-06,0.150E-06,
     .  0.124E-06,0.103E-06,0.852E-07,0.704E-07,0.581E-07,0.479E-07,
     .  0.394E-07,0.324E-07,0.266E-07,0.218E-07,0.178E-07,0.145E-07,
     .  0.119E-07,0.965E-08,0.785E-08,0.638E-08,0.518E-08,0.420E-08,
     .  0.338E-08,0.258E-08,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA DIST3/0.867E-25,0.381E-11,0.317E-05,0.979E-03,
     .  0.909E-02,0.133E-01,0.770E-02,0.369E-02,0.191E-02,0.109E-02,
     .  0.668E-03,0.457E-03,0.306E-03,0.205E-03,0.153E-03,0.124E-03,
     .  0.909E-04,0.637E-04,0.479E-04,0.392E-04,0.329E-04,0.259E-04,
     .  0.195E-04,0.149E-04,0.118E-04,0.993E-05,0.832E-05,0.652E-05,
     .  0.509E-05,0.406E-05,0.328E-05,0.269E-05,0.230E-05,0.200E-05,
     .  0.161E-05,0.126E-05,0.100E-05,0.816E-06,0.669E-06,0.552E-06,
     .  0.457E-06,0.378E-06,0.314E-06,0.261E-06,0.217E-06,0.181E-06,
     .  0.151E-06,0.126E-06,0.105E-06,0.880E-07,0.736E-07,0.615E-07,
     .  0.515E-07,0.431E-07,0.361E-07,0.302E-07,0.253E-07,0.212E-07,
     .  0.177E-07,0.149E-07,0.124E-07,0.104E-07,0.870E-08,0.727E-08,
     .  0.608E-08,0.508E-08,0.425E-08,0.356E-08,0.298E-08,0.251E-08,
     .  0.211E-08,0.177E-08,0.148E-08,0.116E-08,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA DIST4/0.776E-25,0.360E-11,0.306E-05,0.958E-03,
     .  0.895E-02,0.131E-01,0.764E-02,0.368E-02,0.192E-02,0.111E-02,
     .  0.683E-03,0.472E-03,0.319E-03,0.215E-03,0.162E-03,0.133E-03,
     .  0.984E-04,0.697E-04,0.529E-04,0.436E-04,0.369E-04,0.293E-04,
     .  0.223E-04,0.171E-04,0.137E-04,0.115E-04,0.971E-05,0.765E-05,
     .  0.599E-05,0.479E-05,0.388E-05,0.318E-05,0.271E-05,0.236E-05,
     .  0.189E-05,0.148E-05,0.118E-05,0.953E-06,0.780E-06,0.641E-06,
     .  0.528E-06,0.436E-06,0.361E-06,0.299E-06,0.248E-06,0.206E-06,
     .  0.172E-06,0.143E-06,0.119E-06,0.996E-07,0.832E-07,0.697E-07,
     .  0.584E-07,0.490E-07,0.411E-07,0.346E-07,0.291E-07,0.245E-07,
     .  0.206E-07,0.174E-07,0.147E-07,0.124E-07,0.105E-07,0.883E-08,
     .  0.746E-08,0.630E-08,0.532E-08,0.450E-08,0.380E-08,0.322E-08,
     .  0.272E-08,0.231E-08,0.196E-08,0.166E-08,0.142E-08,0.121E-08,
     .  0.103E-08,0.863E-09,0.640E-09,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA DIST5/0.715E-25,0.345E-11,0.299E-05,0.942E-03,
     .  0.884E-02,0.130E-01,0.760E-02,0.367E-02,0.193E-02,0.112E-02,
     .  0.693E-03,0.482E-03,0.328E-03,0.223E-03,0.169E-03,0.140E-03,
     .  0.104E-03,0.743E-04,0.567E-04,0.472E-04,0.402E-04,0.321E-04,
     .  0.246E-04,0.190E-04,0.153E-04,0.130E-04,0.110E-04,0.870E-05,
     .  0.685E-05,0.550E-05,0.447E-05,0.368E-05,0.314E-05,0.273E-05,
     .  0.220E-05,0.171E-05,0.137E-05,0.111E-05,0.903E-06,0.741E-06,
     .  0.609E-06,0.501E-06,0.413E-06,0.341E-06,0.282E-06,0.234E-06,
     .  0.194E-06,0.161E-06,0.134E-06,0.111E-06,0.928E-07,0.775E-07,
     .  0.648E-07,0.543E-07,0.455E-07,0.382E-07,0.321E-07,0.271E-07,
     .  0.228E-07,0.193E-07,0.163E-07,0.138E-07,0.117E-07,0.990E-08,
     .  0.840E-08,0.713E-08,0.606E-08,0.515E-08,0.438E-08,0.372E-08,
     .  0.317E-08,0.270E-08,0.230E-08,0.196E-08,0.167E-08,0.143E-08,
     .  0.123E-08,0.105E-08,0.906E-09,0.781E-09,0.665E-09,0.534E-09,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA DIST6/0.679E-25,0.335E-11,0.294E-05,0.931E-03,
     .  0.875E-02,0.129E-01,0.755E-02,0.367E-02,0.193E-02,0.112E-02,
     .  0.700E-03,0.489E-03,0.335E-03,0.229E-03,0.175E-03,0.145E-03,
     .  0.109E-03,0.780E-04,0.599E-04,0.501E-04,0.429E-04,0.345E-04,
     .  0.266E-04,0.206E-04,0.167E-04,0.142E-04,0.121E-04,0.966E-05,
     .  0.764E-05,0.616E-05,0.503E-05,0.416E-05,0.356E-05,0.311E-05,
     .  0.251E-05,0.196E-05,0.156E-05,0.127E-05,0.103E-05,0.848E-06,
     .  0.697E-06,0.573E-06,0.471E-06,0.388E-06,0.320E-06,0.264E-06,
     .  0.218E-06,0.181E-06,0.150E-06,0.124E-06,0.103E-06,0.859E-07,
     .  0.716E-07,0.598E-07,0.500E-07,0.419E-07,0.351E-07,0.295E-07,
     .  0.249E-07,0.210E-07,0.177E-07,0.150E-07,0.127E-07,0.108E-07,
     .  0.916E-08,0.779E-08,0.663E-08,0.565E-08,0.482E-08,0.412E-08,
     .  0.352E-08,0.301E-08,0.257E-08,0.220E-08,0.189E-08,0.162E-08,
     .  0.139E-08,0.119E-08,0.103E-08,0.885E-09,0.765E-09,0.663E-09,
     .  0.572E-09,0.484E-09,0.345E-09,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA DIST7/0.643E-25,0.326E-11,0.289E-05,0.920E-03,
     .  0.867E-02,0.128E-01,0.752E-02,0.366E-02,0.193E-02,0.113E-02,
     .  0.706E-03,0.495E-03,0.340E-03,0.233E-03,0.179E-03,0.150E-03,
     .  0.112E-03,0.810E-04,0.625E-04,0.525E-04,0.452E-04,0.366E-04,
     .  0.283E-04,0.221E-04,0.180E-04,0.154E-04,0.132E-04,0.105E-04,
     .  0.837E-05,0.678E-05,0.555E-05,0.461E-05,0.397E-05,0.347E-05,
     .  0.281E-05,0.220E-05,0.176E-05,0.143E-05,0.117E-05,0.960E-06,
     .  0.789E-06,0.648E-06,0.533E-06,0.439E-06,0.361E-06,0.298E-06,
     .  0.245E-06,0.203E-06,0.167E-06,0.138E-06,0.115E-06,0.950E-07,
     .  0.790E-07,0.657E-07,0.548E-07,0.457E-07,0.383E-07,0.321E-07,
     .  0.269E-07,0.227E-07,0.191E-07,0.161E-07,0.137E-07,0.116E-07,
     .  0.984E-08,0.837E-08,0.713E-08,0.609E-08,0.520E-08,0.445E-08,
     .  0.381E-08,0.326E-08,0.280E-08,0.240E-08,0.206E-08,0.177E-08,
     .  0.152E-08,0.131E-08,0.113E-08,0.975E-09,0.843E-09,0.730E-09,
     .  0.634E-09,0.552E-09,0.479E-09,0.408E-09,0.308E-09,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA DIST8/0.593E-25,0.312E-11,0.282E-05,0.903E-03,
     .  0.855E-02,0.126E-01,0.745E-02,0.364E-02,0.193E-02,0.114E-02,
     .  0.714E-03,0.504E-03,0.349E-03,0.241E-03,0.186E-03,0.156E-03,
     .  0.118E-03,0.859E-04,0.668E-04,0.564E-04,0.490E-04,0.399E-04,
     .  0.311E-04,0.245E-04,0.201E-04,0.173E-04,0.149E-04,0.120E-04,
     .  0.965E-05,0.787E-05,0.650E-05,0.543E-05,0.471E-05,0.415E-05,
     .  0.338E-05,0.267E-05,0.215E-05,0.175E-05,0.144E-05,0.119E-05,
     .  0.977E-06,0.805E-06,0.663E-06,0.546E-06,0.450E-06,0.370E-06,
     .  0.304E-06,0.251E-06,0.206E-06,0.170E-06,0.140E-06,0.116E-06,
     .  0.955E-07,0.790E-07,0.655E-07,0.543E-07,0.452E-07,0.376E-07,
     .  0.314E-07,0.263E-07,0.221E-07,0.185E-07,0.156E-07,0.132E-07,
     .  0.112E-07,0.947E-08,0.806E-08,0.686E-08,0.586E-08,0.501E-08,
     .  0.430E-08,0.369E-08,0.317E-08,0.273E-08,0.235E-08,0.203E-08,
     .  0.175E-08,0.151E-08,0.131E-08,0.113E-08,0.978E-09,0.847E-09,
     .  0.735E-09,0.639E-09,0.557E-09,0.486E-09,0.425E-09,0.372E-09,
     .  0.322E-09,0.263E-09,0.000E+00,0.000E+00,0.000E+00,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA DIST9/0.558E-25,0.302E-11,0.276E-05,0.890E-03,
     .  0.845E-02,0.125E-01,0.740E-02,0.363E-02,0.193E-02,0.114E-02,
     .  0.720E-03,0.511E-03,0.355E-03,0.246E-03,0.191E-03,0.161E-03,
     .  0.123E-03,0.897E-04,0.700E-04,0.595E-04,0.519E-04,0.426E-04,
     .  0.334E-04,0.264E-04,0.218E-04,0.189E-04,0.164E-04,0.133E-04,
     .  0.107E-04,0.882E-05,0.732E-05,0.616E-05,0.536E-05,0.475E-05,
     .  0.390E-05,0.310E-05,0.251E-05,0.206E-05,0.170E-05,0.141E-05,
     .  0.116E-05,0.963E-06,0.796E-06,0.657E-06,0.542E-06,0.447E-06,
     .  0.368E-06,0.303E-06,0.249E-06,0.205E-06,0.168E-06,0.139E-06,
     .  0.114E-06,0.941E-07,0.776E-07,0.641E-07,0.530E-07,0.440E-07,
     .  0.365E-07,0.304E-07,0.253E-07,0.212E-07,0.177E-07,0.149E-07,
     .  0.126E-07,0.106E-07,0.898E-08,0.762E-08,0.649E-08,0.554E-08,
     .  0.474E-08,0.406E-08,0.349E-08,0.300E-08,0.259E-08,0.223E-08,
     .  0.193E-08,0.167E-08,0.145E-08,0.126E-08,0.109E-08,0.946E-09,
     .  0.822E-09,0.715E-09,0.622E-09,0.542E-09,0.474E-09,0.415E-09,
     .  0.364E-09,0.320E-09,0.280E-09,0.241E-09,0.175E-09,0.000E+00,
     .  0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00,0.000E+00/
        DATA DMAX/.131E-01,.135E-01,.133E-01,.131,.13E-01,.129E-01,
     .  .128E-01,.126E-01,.125E-01/
C	NEXT LINE HAS TO BE CHANGED FOR GRAPHITE
	DEav=deav+A1C+A2C*EMU
c	!Fluctuate energy loss 
		DO 100 J=2,9
		IF(EMUS(J).GT.EMU) GO TO 200
 100		CONTINUE
		IF(J.GT.9) J=9
 200		I=J-1
		F=(EMU-EMUS(I))/(EMUS(J)-EMUS(I))
C   
C     TO GENERATE THE INVERSE FUNCTIONS 
c   
       IF (FIRST) THEN   
       FIRST=.FALSE. 
       DO 250 L=1,9  
       DINV(1,L)=0.  
       SUM=0.    
       DO 260 LL=2,NENR(L)   
       SUM=SUM+(DIST(LL-1,L)+DIST(LL,L))*0.5*(ELOS(LL,L)-ELOS(LL-1,L))   
       DINV(LL,L)=SUM    
 260   CONTINUE  
       DO 270 LI=2,NENR(L)   
       DINV(LI,L)=DINV(LI,L)/SUM 
 270   CONTINUE  
 250   CONTINUE  
       ENDIF 
C   
C     TO FIND THE ENERGY LOSS FOR THE TWO DISTRIBUTIONS 
C   
       RAND=Ran(iseed)   
       DO 280 IJ=I,J 
       DO 290 IK=2,NENR(IJ)  
       IF (RAND.LT.DINV(IK,IJ)) GO TO 291    
 290   CONTINUE  
 291   CONTINUE  
       I1=IK-1   
       I2=IK 
       EL(IJ)=ELOS(I1,IJ)+(ELOS(I2,IJ)-ELOS(I1,IJ))* 
     1 (RAND-DINV(I1,IJ))/(DINV(I2,IJ)-DINV(I1,IJ)) 
 280   CONTINUE  
C   
C     INTERPOLATE THE ENERGY LOSS   
C   
      ELOSS=EL(I)+F*(EL(J)-EL(I))   
             DEmu=ELOSS/1000.  
C   
	return
	end
