!!!SF is 40 Em bins of 5 MeV (0-200 MeV) by 40 Pm bins of 20 MeV (0-800 MeV/c)

!SF is 80 Em bins of 5 MeV             by 40 Pm bins of 20 MeV.

destroy *
read\scalar benharsf_?1.dat pbins ebins
read benharsf_?1.dat\(2) pmbin embin sfp sfn dp de

oneforproton_twoforneutron=?2
whichpart=?2
if (whichpart=1) then
 sfbin = sfp
endif
if (whichpart=2) then
 sfbin = sfn
endif

stat\nomess sfbin normfac\sum

sort embin pmbin sfbin
em = embin[1:#:pbins]
rebin sfbin sf_em pbins

sort pmbin embin sfbin
pm = pmbin[1:#:ebins]
rebin sfbin sf_pm ebins

clear
portrait
window 3
set ylog 10
xlabel `Pm'
pch 2 .5
auto on
gr pm sf_pm

copy sf_pm sfsum_pm
do ind=[2:len(pm):1]
 stat\nomess sf_pm[1:ind] sfsum_pm[ind]\sum
enddo
sfsum_pm=sfsum_pm/sfsum_pm[len(pm)]

set %xloc 50
set %yloc 80
set cursor -7
text `Normalization = '//rchar(nint(10000*normfac)/10000)

window 4
xlabel `Em'
set ylog 0
!gr em[2:#] sf_em[2:#]	!1st bin is zero
gr em sf_em

terminal

cl
set ylog 10
landscape
window 0
nbins=10		!must be mult of 40
scales 0 400 8 -8 0 8
graph\axesonly
do ind=[1:nbins:1]
 pmlo=nint(800/nbins)*(ind-1)+1
 pmhi=nint(800/nbins)*ind
 destroy x y
 copy embin sfbin xtmp ytmp IFF (pmbin>=pmlo & pmbin<=pmhi)
 sort xtmp ytmp
 x=xtmp[1:#:nint(40/nbins)]
 rebin ytmp y nint(40/nbins)
 pch ind 1 ind
 destroy x y IFF (y<1.d-30)

 stat\nomess y ysum\sum
 y=y/ysum		!fixed normalization

 cur x y
enddo
!replot