macro all

exec init deg=1
exec makehist nbin=400 xmin=0 xmax=4 fact=1.0
exec savehist file=fluxkorea1_hist.hbk

exec init deg=15
exec makehist nbin=400 xmin=0 xmax=4 fact=1.0
exec savehist file=fluxkorea15_hist.hbk

exec init deg=2
exec makehist nbin=400 xmin=0 xmax=4 fact=1.0
exec savehist file=fluxkorea2_hist.hbk

exec init deg=25
exec makehist nbin=400 xmin=0 xmax=4 fact=1.0
exec savehist file=fluxkorea25_hist.hbk

return
**********************************************************************
macro savehist file=fluxkorea25_hist.hbk

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

return
**********************************************************************
macro init deg=25

nt/chain -nubeam 
do i=1,16
  file = fluxkorea[deg].[i].hbk
  nt/chain nubeam [file]
enddo
uwfunc //nubeam/2000 jnubeam.inc

global/creat pot 16.0E+21

return
**********************************************************************
macro makehist nbin=100 xmin=0 xmax=10 

global/import *
cd //pawc
h/del 0

* numu
h/cre/1d 10 'E[nm] from [p]+,K+'  [nbin] [xmin] [xmax]
h/cre/1d 11 'E[nm] from [p]+'     [nbin] [xmin] [xmax]
h/cre/1d 12 'E[nm] from K+'       [nbin] [xmin] [xmax]

h/cre/1d 210 '[nm] N(SK) (CC,NC) '   [nbin] [xmin] [xmax]
h/cre/1d 310 '[nm] N(SK) (CC)    '  [nbin] [xmin] [xmax]
h/cre/1d 215 '[nm] to [n]e N(SK) (CC,NC) '  [nbin] [xmin] [xmax]
h/cre/1d 315 '[nm] to [n]e N(SK) (CC)    '  [nbin] [xmin] [xmax]

* numu_bar
h/cre/1d 20 'E[nm]_bar from [p]-,K-'  [nbin] [xmin] [xmax]
h/cre/1d 21 'E[nm]_bar from [p]-'     [nbin] [xmin] [xmax]
h/cre/1d 22 'E[nm]_bar from K-'       [nbin] [xmin] [xmax]

h/cre/1d 220 '[nm]_bar N(SK) (CC,NC)'  [nbin] [xmin] [xmax]
h/cre/1d 320 '[nm]_bar N(SK) (CC)   '  [nbin] [xmin] [xmax]
h/cre/1d 225 '[nm]_bar to [n]e_bar N(SK) (CC,NC)'  [nbin] [xmin] [xmax]
h/cre/1d 325 '[nm]_bar to [n]e_bar N(SK) (CC)   '  [nbin] [xmin] [xmax]

* nue
h/cre/1d 30 'E[n]e from all'  [nbin] [xmin] [xmax]
h/cre/1d 31 'E[n]e from K+'   [nbin] [xmin] [xmax]
h/cre/1d 32 'E[n]e from K0L'  [nbin] [xmin] [xmax]
h/cre/1d 33 'E[n]e from [m]+' [nbin] [xmin] [xmax]

h/cre/1d 230 '[n]e N(SK) (CC,NC)'   [nbin] [xmin] [xmax]
h/cre/1d 330 '[n]e N(SK) (CC)'  [nbin] [xmin] [xmax]

* nue_bar
h/cre/1d 40 'E[n]e_bar from all'  [nbin] [xmin] [xmax]
h/cre/1d 41 'E[n]e_bar from K-'   [nbin] [xmin] [xmax]
h/cre/1d 42 'E[n]e_bar from K0L'  [nbin] [xmin] [xmax]
h/cre/1d 43 'E[n]e_bar from [m]-' [nbin] [xmin] [xmax]

h/cre/1d 240 '[n]e_bar N(SK) (CC,NC)'   [nbin] [xmin] [xmax]
h/cre/1d 340 '[n]e_bar N(SK) (CC)'  [nbin] [xmin] [xmax]

* Ppi from numu
h/cre/1d 110 '[p]+ spectrum contributed nu at SK'  [nbin] [xmin] [xmax]
h/cre/1d 111 '[p]+ spectrum contributed nu at SK'  [nbin] [xmin] [xmax]
h/cre/1d 112 '[p]+ spectrum contributed nu at SK'  [nbin] [xmin] [xmax]
* Ppi from numu_bar
h/cre/1d 120 '[p]- spectrum contributed nu at SK'  [nbin] [xmin] [xmax]
h/cre/1d 121 '[p]- spectrum contributed nu at SK'  [nbin] [xmin] [xmax]
h/cre/1d 122 '[p]- spectrum contributed nu at SK'  [nbin] [xmin] [xmax]

Appli comis quit
      function fill(fact)
      include 'jnubeam.inc'
      fact=fact*1.0E-6
      weight=fact*norm
      if ( mode.ge.10 .and. mode.le.12 ) then
         call hf1(10,enu,weight)
         call hf1(110,ppi,weight)
         if ( mode.eq.11 ) then
            call hf1(11,enu,weight)
            call hf1(111,ppi,weight)
         endif
         if ( mode.eq.12 ) then 
            call hf1(12,enu,weight)
            call hf1(112,ppi,weight)
         endif
      else if ( mode.ge.20 .and. mode.le.22 ) then
         call hf1(20,enu,weight)
         call hf1(120,ppi,weight)
         if ( mode.eq.21 ) then 
            call hf1(21,enu,weight)
            call hf1(121,ppi,weight)
         endif
         if ( mode.eq.22 ) then 
            call hf1(22,enu,weight)
            call hf1(22,ppi,weight)
         endif
      else if ( mode.ge.30 .and. mode.le.33 ) then
         call hf1(30,enu,weight)
         call hf1(130,ppi,weight)
         if ( mode.eq.31 ) call hf1(31,enu,weight)
         if ( mode.eq.32 ) call hf1(32,enu,weight)
         if ( mode.eq.33 ) call hf1(32,enu,weight)
      else if ( mode.ge.40 .and. mode.le.43 ) then
         call hf1(40,enu,weight)
         if ( mode.eq.41 ) call hf1(41,enu,weight)
         if ( mode.eq.42 ) call hf1(42,enu,weight)
         if ( mode.eq.43 ) call hf1(42,enu,weight)
      endif
      return
      end
quit

*
* flux unit is /10^21 POT/cm^2
*
fact = $rsigma(1.0E+21/[pot])
nt/loop //nubeam/2000 fill([fact]) 

return
**********************************************************************