macro all

exec makehist input=fluxkorea1_hist.hbk  output=prob2event_korea1_4gev.hbk
exec makehist input=fluxkorea15_hist.hbk output=prob2event_korea15_4gev.hbk
exec makehist input=fluxkorea2_hist.hbk  output=prob2event_korea2_4gev.hbk
exec makehist input=fluxkorea25_hist.hbk output=prob2event_korea25_4gev.hbk

return
**********************************************************************
macro makehist input=fluxkorea25_hist.hbk output=prob2event_korea25_4gev.hbk

*
*   kumac file to make event spectrum of (numu flux) X (cross section)
*   at kamioka and korea site
*

*
* import numu flux 
*   flux unit should be 10^6/cm^2/10^21POT
*   10^21 POT corresponds to 0.75MWx1yr
*
h/file 1 [input]
hrin 10  ! 10000
close 1

*
* import nue and anti-nue cross section
* this histogram is created by 'xsec' program
*  
*  histogram:
*    hid  103 : sigma(nue+n->(e-)+p)*0.444
*    hid 1103 : sigma(nueb+p->(e+)+p)*(~0.566)
*      0.444 : fraction of neutron in H2O
*     ~0.566 : fraction of proton in H2O (include free proton effect)
*  unit is 10^(-38)xcm^2
*
h/file 1 xsec_0-4gev_400bin.hbk
hrin 103  ! 20000 | nue CCqe
hrin 1103 ! 20000  | nue-bar CCqe
close 1

*
*  target mass
*
mass = 1.0E+12                   | 1Mton=10^6ton=10^12g
target = [mass]*6.02E+23         | number of nucleon in mass
fact_mw = 1.0/0.75               | 0.75MW -> 1MW
fact = [target]*[fact_mw]*1.0E-38*1.0E+6
mess fact is [fact]

*multi 10010 20101 991      | nue 
*multi 10010 21101 992     | nue-bar
multi 10010 20103 991      | nue 
multi 10010 21103 992     | nue-bar

* 
*   this normalization method is obsolute (26-Sep-2006)
* *
* *  normalize to 1Mtonx1yr event spectrum (from nunokawa-table)
* *    kam   nue   37148.4 ev/1Mton/1yr   (0-1.5GeV)
* *    kam   nueb  10874.8 ev/1Mton/1yr   (0-1.5GeV)
* *    korea nue    2932.3 ev/1Mton/1yr   (0-1.5GeV)
* *    korea nueb    858.4 ev/1Mton/1yr   (0-1.5GeV)
* *
* *
* *    sum in hid 991 in 0-1.5 GeV   7.723114    (OA2.5deg,fluxkorea25_hist.hbk)
* *    sum in hid 992 in 0-1.5 GeV   2.708750    (OA2.5deg,fluxkorea25_hist.hbk)
* *
* add 991 991 101 $eval(2932.3/7.723114) 0 
* add 992 992 102 $eval( 858.4/2.708750) 0 
* 
add 991 991 101 [fact] 0 
add 992 992 102 [fact] 0 

* convert x-axis from GeV to MeV
nbin = $hinfo(101,'xbins')
v/cre tmp3([nbin])
v/cre tmp4([nbin])
get/con 101 tmp3 
get/con 102 tmp4 
h/del 101,102
h/cre/1d 101 'event spec. at korea nue' 400 0 4000
h/cre/1d 102 'event spec. at korea nue' 400 0 4000
put/con 101 tmp3 
put/con 102 tmp4 



h/file 2  [output] ! 'N'
hrout 0 
close 2


return