* create input histogram file for nuance.

* This assumes vectors give Neutrinos/GeV/m^2/PoT at 1 km from the
* source.  Columns of: E(eV), Numuflux, Nueflux, Anumuflux.
 
* The resulting histograms will be scalled to give
* neutrinos/BIN_WIDTH/cm^2/PoT and reduced to the given baseline

* assume 10MeV/bin
 macro loadvec filename=se-00-all.vec
   close 0
   ve/del *
   hi/del *
   ve/read a,b,c,d [filename]
   exec vec#sparsify a 1 10 e
   exec vec#sparsify b 1 10 numu
   exec vec#sparsify c 1 10 nue
   exec vec#sparsify d 1 10 anumu

   sigma gev=e/1e9

 return
 
 macro genhist filename=nuance-input.hbk hid1=1 hid2=2 hid3=3
* bl^2*(100cm/m)^2*(100 10MeV/GeV)
   
   bin=gev(2)-gev(1)
   len=$vdim(e)
   e1=gev(1)
   e2=gev(2)
   elast=gev([len])
   ediff=[e2]-[e1]
   eend=[elast]+[ediff]

   message 'Bin size: '//[ediff]//' GeV'

   pot_per_yr=1.1e21		| with a 0.5 MW, 28 GeV beam and 10^7 sec/yr

   reduction=10000*(1.0/[ediff])/[pot_per_yr]
   sigma numu10=numu/[reduction]
   sigma nue10=nue/[reduction]
   sigma anumu10=anumu/[reduction]

   1d [hid1] 'Numu/bin/cm**2/yr at 1 km'     [len] [e1] [eend]
   1d [hid2] 'Nue/bin/cm**2/yr at 1 km'      [len] [e1] [eend]
   1d [hid3] 'Antinumu/bin/cm**2/yr at 1 km' [len] [e1] [eend]
   
   put/cont [hid1] numu10
   put/cont [hid2] nue10
   put/cont [hid3] anumu10
   
   hi/file 1 [filename] ! N
   hrout [hid1]
   hrout [hid2]
   hrout [hid3]
   close 1

 return
 
 macro plot 
   opt logy
   
   ystart=$hinfo(1,'max')
   ymax=[ystart]*1.2
   ymin=$hinfo(1,'min')
* ymin=[ymax]*1.0e-4
   xmax=$hinfo(1,'xmax')
   xmin=$hinfo(1,'xmin')
   nbins=$hinfo(1,'xbins')
   xstart=[xmax]*0.7
   binsiz=$format(([xmax]-[xmin])/[nbins],F3.1)

   title 'P889 0^o! beam components'
* null 0 [xmax] [ymin] [ymax]
* null 0 10 .9e7 5e11
   null 0 10 4e7 1e12
   exec multi#plot-dmod 0 1 1 2 2 4 3
   exec multi#line 'log' [xstart] [ystart] 2.0 1 1 '[n]?[m]!' 2 2 '[n]?e!'  4 3 '[n]&^__!?[m]!'
   atitle 'E?[n]! (GeV)' '[n]/'//[binsiz]//'GeV/cm^2!/10^7!sec/0.5MW' ' ' 321
 return