macro loadvecs
   ve/del *
   ve/read r,d
   ve/read e,x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x16,x17,x18 bnlhs/e889-xsec.vec
  
   ve/read e,nf00 bnlhs/numi-flux-00.vec
   ve/read e,nf10 bnlhs/numi-flux-10.vec
   ve/read e,bf00 bnlhs/e889-flux-00.vec
   ve/read e,bf15 bnlhs/e889-flux-15.vec

   ve/read e,a,b,c bnlhs/prob_0.8_1.0_0.05_5.0e-5_3.0e-3.vec
   sigma gev=e/1e9
 return
 
 macro plot-density
   sigma rkm=r/1e5
   exec vec#tohist 1 rkm d
   title 'PREM Earth Density'
   null 0 6346 0 14
   set dmod 1
   hi/pl 1 s
   atitle 'Earth radius (km)' 'Density (g/cm^3!)' '' 322
   hi/del 1
 return
 
 
 macro plot-flux
   exec vec#tohist 1 gev bf00
   exec vec#tohist 2 gev bf15
   exec vec#tohist 3 gev nf00
   exec vec#tohist 4 gev nf10
   
   ymax=5e-17
   ystart=0.9*[ymax]
   ystep=0.1*[ymax]

   title '1MW [n] flux at BNL-HS distance'
   null 0 10 0 [ymax]
   exec multi#plot-dmod 0 1 1 2 2 6 3 4 4
   exec multi#line lin 6 [ystart] [ystep] 1 1 'P889 0^o!' 2 2  'P889 1^o!' 6 3 'NuMI 0^o!' 4 4 'NuMI 1^o!'
   atitle 'E?[n]!(GeV)' '[n]/s/cm^2!/PoT/0.01GeV' '' 222

   hi/del 1
   hi/del 2
   hi/del 3
   hi/del 4
 return
 
 macro plot-prob xmax=10
   sigma a10=a*10
   exec vec#tohist 1 gev a10
   exec vec#tohist 2 gev b
   exec vec#tohist 3 gev c
   
   title 'BNL-HS 3 [n] osc, PREM'
   null 0 [xmax] 0 1
   exec multi#plot-dmod 0 1 1 2 2 4 3 
   exec multi#line lin 6 0.28 .05 1 1 '10*P([n]?[m]!"5#[n]?e!)' 2 2 'P([n]?[m]!"5#[n]?[m]!)' 4 3 'P([n]?[m]!"5#[n]?[t]!)'
   exec multi#itx 'lin' 5.0 0.65 0.06 1 'sin^2!2[q]?ij! = 0.8/1.0/0.4' 1 '[D]m?ij!^2! = 5e-5/3e-3 eV^2!'
   atitle 'E?[n]!(GeV)' 'Probability' '' 222
   ve/del a10
   hi/del 1
   hi/del 2
   hi/del 3
 return
 
 
 
 macro xsec-hists
   for n in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 16 17 18
     hid=100+[n]
     exec vec#tohist [hid] gev x[n]
   endfor
 return
 
 macro plot-all-xsec
   null 0 10 0 1.2
   for n in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 17 18
     hid=100+[n]
     hi/pl [hid] s
   endfor
 return

 macro plot-xsec
   title '[n]?l!N"5#l^-!P cross section'
   null 0 10 0 1.2
   exec vec#tohist 1 gev x2
   hi/pl 1 s
   atitle 'E?[n]!(GeV)' '10^-38!cm^2!' '' 222
   hi/del 1
 return
 
 macro plot-events flux desc xmax=10
* Na * gram/kton * 10 Neutrons per H2O * units of xsec
   conv=6.02e23*1e9*10.0*1e-38
   sigma mucc=[flux]*x2*[conv]
   sigma oscmucc=[flux]*b*x2*[conv]
   sigma oscecc=[flux]*a*x17*10*[conv]

   exec vec#tohist-centered 1 gev mucc
   exec vec#tohist-centered 2 gev oscmucc
   exec vec#tohist-centered 3 gev oscecc

   ymax=3e-21
   ystart=[ymax]*0.9
   ydiff=[ymax]*0.05

   title '[n] CC event rate, '//[desc]
   null 0 [xmax] 0 [ymax]
   exec multi#plot-dmod 0 1303 3 1 1 2 2 
   exec multi#line 'lin' 2 [ystart] [ydiff] 1 2 '[n]?[m] CC' 2 3 '[n]?[m]! CC osc.' 3 1 '10*[n]?e! CC osc.'
   atitle 'E?[n]!(GeV)' 'Events/KTon/PoT/0.01GeV' '' 222
   hi/del 1
   hi/del 2
   hi/del 3
 return
 
 macro plot-events-norm flux desc detmass=100 npotpspil=4e14 spilltime=1.6 secperyr=3e7
   norm=[detmass]*[npotpspil]*[secperyr]/[spilltime]
* Na * gram/kton * 10 Neutrons per H2O * units of xsec
   conv=6.02e23*1e9*10.0*1e-38*[norm]
   sigma mucc=[flux]*x2*[conv]
   sigma oscmucc=[flux]*b*x2*[conv]
   sigma oscecc=[flux]*a*x17*10*[conv]

   exec vec#tohist-centered 1 gev mucc
   exec vec#tohist-centered 2 gev oscmucc
   exec vec#tohist-centered 3 gev oscecc

   ymax=2000
   ytext=[ymax]*0.9
   ystart=[ymax]*0.5
   ydiff=[ymax]*0.05

   title '[n] CC event rate, '//[desc]
   null 0 10 0 [ymax]
   
   exec multi#itx 'lin' 3 [ytext] [ydiff] 1 [detmass]//' kton water target' 1 [npotpspil]//' PoT/spill' 1 [spilltime]//' sec spill' 1 [secperyr]//' sec/year'

   exec multi#plot-dmod 0 1303 3 1 1 2 2 
   exec multi#line 'lin' 6 [ystart] [ydiff] 1 1 '[n]?[m] CC' 2 2 '[n]?[m]! CC osc.' 3 3 '10*[n]?e! CC osc.'
   atitle 'E?[n]!(GeV)' 'Events/Year/10MeV' '' 222
   hi/del 1
   hi/del 2
   hi/del 3
 return
 
 macro plot-all
   fort/file 69 bnlhs-plots.ps
   meta 69 -114
   exec plot-density
   exec plot-prob
   exec plot-prob 4

   exec plot-flux
   exec plot-xsec
   exec plot-events bf00 '1MW P889 0^o!' 10
   exec plot-events bf00 '1MW P889 0^o!' 3
   exec plot-events bf15 '1MW P889 1.5^o!' 10
   exec plot-events bf15 '1MW P889 1.5^o!' 3
   exec plot-events nf00 '1MW NuMI 0^o!' 10
   exec plot-events nf00 '1MW NuMI 0^o!' 3
   exec plot-events nf10 '1MW NuMI 1^o!' 10
   exec plot-events nf10 '1MW NuMI 1^o!' 3
   
   exec plot-events-norm bf00 '1MW AGS P889 0^o!' 100 4.2e14 1.6
   exec plot-events-norm bf15 '1MW AGS P889 1.5^o!' 100 4.2e14 1.6
   exec plot-events-norm nf00 '1MW NuMI 0^o!' 100 9.7e13 1.9
   exec plot-events-norm nf10 '1MW NuMI 1^o!' 100 9.7e13 1.9
   
   close 69
 return