1 gaskelld 1.1 !!!SF is 40 Em bins of 5 MeV (0-200 MeV) by 40 Pm bins of 20 MeV (0-800 MeV/c)
  2              
  3              !SF is 80 Em bins of 5 MeV             by 40 Pm bins of 20 MeV.
  4              
  5              destroy *
  6              pi=3.14159265
  7              read\scalar benharsf_?1.dat pbins ebins
  8              read benharsf_?1.dat\(2) pmbin embin sfpbin sfnbin dp de
  9              
 10              stat\nomess sfpbin normfacp\sum
 11              stat\nomess sfnbin normfacn\sum
 12              sfpbin=sfpbin/normfacp
 13              sfnbin=sfnbin/normfacn
 14              
 15              sort embin pmbin sfpbin sfnbin
 16              em = embin[1:#:pbins]
 17              rebin sfpbin sfp_em pbins
 18              rebin sfnbin sfn_em pbins
 19              
 20              sort pmbin embin sfpbin sfnbin
 21              pm = pmbin[1:#:ebins]
 22 gaskelld 1.1 rebin sfpbin sfp_pm ebins
 23              rebin sfnbin sfn_pm ebins
 24              
 25              clear
 26              landscape
 27              window 5
 28              set ylog 10
 29              xlabel `Pm'
 30              ylabel `n(p) d<^>3<_>p'
 31              pch 2 .8 1
 32              auto on
 33              scales 0 1000 5 -4 0 4
 34              gr pm sfp_pm
 35              pch 12 .8 2
 36              cur pm sfn_pm
 37              
 38              copy sfp_pm sfpsum_pm
 39              copy sfn_pm sfnsum_pm
 40              do ind=[2:len(pm):1]
 41               stat\nomess sfp_pm[1:ind] sfpsum_pm[ind]\sum
 42               stat\nomess sfn_pm[1:ind] sfnsum_pm[ind]\sum
 43 gaskelld 1.1 enddo
 44              sfpsum_pm=sfpsum_pm/sfpsum_pm[len(pm)]
 45              sfnsum_pm=sfnsum_pm/sfnsum_pm[len(pm)]
 46              
 47              set %xloc 35
 48              set %yloc 85
 49              set cursor -7
 50              text `Raw Proton Norm. = '//rchar(nint(10000*normfacp)/10000)
 51              set %yloc 80
 52              text `Raw Neutron Norm. = '//rchar(nint(10000*normfacn)/10000)
 53              
 54              
 55              
 56              
 57              window 6
 58              k=pm
 59              nk_p=sfp_pm/pm**2
 60              nk_n=sfn_pm/pm**2
 61              
 62              
 63              !interpolate to get more bins.  Since each bin is a probability, we need to
 64 gaskelld 1.1 ! scale by the old number of bins vs. the new number of bins.
 65              
 66              xx=k
 67              yy=nk_p
 68              yy=yy*10**(xx/160)		!minimize the variation.
 69              yy=yy*(len(xx)/600)
 70              gen xxnew 0,,xx[#] 600
 71              yynew=interp(xx,yy,xxnew)
 72              yynew=yynew/10**(xxnew/160)
 73              k=xxnew
 74              nk_p=yynew
 75              
 76              yy=nk_n
 77              yy=yy*10**(xx/160)		!minimize the variation.
 78              yy=yy*(len(xx)/600)
 79              gen xxnew 0,,xx[#] 600
 80              yynew=interp(xx,yy,xxnew)
 81              yynew=yynew/10**(xxnew/160)
 82              k=xxnew
 83              nk_n=yynew
 84              
 85 gaskelld 1.1 
 86              xlabel `Pm'
 87              ylabel `n(p) d<^>3<_>p/p<^>2<_> = n(p)/<Delta>p'
 88              pch 2 .8 1
 89              scales 0 1000 5 0 1 1
 90              auto yaxis
 91              gr k nk_p
 92              pch 12 .8 2
 93              cur k nk_n
 94              
 95              !** proton normalization	!note: sfp_pm already normalized to PROBABILITY
 96              integrand = k**2*nk_p
 97              intp[1]=integrand[1]		!double count first k step
 98              do ind=[2:len(k)]
 99                intp[ind]=intp[ind-1]+integrand[ind]
100              enddo
101              pnorm1=nint(10000*intp[#])/10000
102              
103              intp=intp/intp[#]
104              
105              !** neutron normalization
106 gaskelld 1.1 integrand = k**2*nk_n
107              intn[1]=integrand[1]
108              do ind=[2:len(k)]
109                intn[ind]=intn[ind-1]+integrand[ind]
110              enddo
111              nnorm1=nint(10000*intn[#])/10000
112              intn=intn/intn[#]
113              
114              set %xloc 40
115              set %yloc 85
116              set cursor -7
117              text `Rebinned p Norm. = '//rchar(pnorm1)
118              set %yloc 80
119              text `Rebinned n Norm. = '//rchar(nnorm1)
120              
121              
122              window 7
123              ylabel `<Int> n(p) d<^>3<_>p'
124              pch 2 .8 1
125              scales 0 1000 5 0 1 1
126              set ylog 0
127 gaskelld 1.1 auto yaxis
128              gr k intp
129              pch 12 .8 2
130              cur k intn
131              
132              
133              window 8
134              ylabel `1 - <Int> n(p) d<^>3<_>p'
135              pch 2 .8 1
136              set ylog 10
137              scales 0 1000 5 -5 0 5
138              gr k 1-intp
139              pch 12 .8 2
140              cur k 1-intn
141              
142              
143              write benhar_pmdist_?1p.dat k intp
144              write benhar_pmdist_?1n.dat k intn
145              
146              terminal
147              
148 gaskelld 1.1 clear
149              window 0
150              znuc=?2
151              anuc=?3
152              nnuc=anuc-znuc
153              y=-k
154              sort y nk_p nk_n
155              fy_p=2*pi*integral(y,abs(y)*nk_p)/Znuc
156              fy_n=2*pi*integral(y,abs(y)*nk_n)/Nnuc
157              fy = (Znuc*fy_p + Nnuc*fy_n) / (Anuc)
158              xlabel `y'
159              ylabel `F(y)'
160              
161              scales -800 0 4 -8 -2 6
162              gr y fy
163              
164              write benharsf_?1.fy y fy