WORK FOR THE NEXT MASS PRODUCTION
What we want to do
We want to generate more neutrino interaction vectors in water with NEUT. Hayato-san kindly sent us the program he uses to take the output of JNUBEAM (neutrino 4-vectors, properly weighted) and turn it into ready-to-use interaction vectors.
* At SK : use the output of JNUBEAM to make energy histograms. then randomize the vertex position uniformly in the cylinder & randomize the energy according to the histogram (CERNLIB's RNHRAN). The neutrino direction is fixed (SK is approx. point-like). Use NEUT to make the interaction for this neutrino. Put into Ntuple. Simple & efficient.
* At 2KM : Loop through each 4-vector ntuple. Retrieve the information event by event. Get the max total cross section Mcr and the maximum weight Mw of the event in the ntuple. If ( random >= weight*total-xsec(of this event) / (Mw*Mcr) ) reject and go to the next event. Use the 4-momentum & vertex position from the 4-vector ntuple. Then make the vector using NEUT, and put into zbs/ntuple file.
CONCLUSION : At SK the 4-vector ntuples have high stats (because every parent hadron is forced to decay there). The distributions are known with a smaller statistical error. The energy spectrum used to sample nu energies is made with weighted events --> deal with the weights in a simple manner. We could apparently generate as many nus as we want, but we are limited by the stat error on the energy spectrum --> Don't forget !!
At 2KM, the situation is more complicated. If we want to increase stats, and keep using this method, then we must first generate more 4-vectors (using JNUBEAM). This method is the best way to correctly deal with position/energy correlations etc, but JNUBEAM is rather inefficient @ 2KM (because of the small solid angle subtended by the 2KM), so we have very few 4-vectors to sample from.
2 options:
* Oversample (probably already used to some extent in the official vectors ?? for nues if not numus ?). Use the same 4-vector several times to make different nu interactions. Effectively makes independant events, is probably effective for almost all studies except the top-down effect (studies that involve nu position vs energy). Then there is no gain.
* Get the PDF of (X,Y). Sample neutrinos from it. Use fixed directions for neutrinos (the spread in direction is smaller than our precision in the reconstruction). BUT we don't want to lose the correlation between (X,Y) and the energy. So based on (X,Y) chose the right energy spectrum to sample from and get the energy according to this spectrum.
I want to avoid (X,Y) binning as much as possible. So 1. Get the (X,Y) histogram with adapted binning (depends on the flavour because the stats are different) 1. Fit it with a plane (should be enough) 1. Bin the detector in (X,Y) and make one energy spectrum per bin 1. Sample using the beam profile PDF (need to normalize them of course). 1. Based on where (X,Y) was found, draw E according to the correct spectrum
FOR numus from pi decays (mode==11 in JNUBEAM) the last two steps could be adapted : make one energy spectrum in a (small) region around the detector center. Get (X,Y). Compute the off-axis value for (X,Y). Draw E according to the center spectrum and correct the energy for the correct off axis angle value. (avoids the binning in the top down plot...).
BEAM PROFILES AND ENERGY PROFILES
Plots with fits :
The beam profile is not flat for all flavours (geometrical effect). However the correlation between E & (x,y) is a by-product of the iff axis trick so it only works for numus from pions. So energy profiles of the beam are mostly flat except for numus.
Here are the fluxes for numu from K and muon decay, bin by bin (each surface bin (in X,Y) corresponds to one spectrum and is assigned a different color). 
BEAM PROFILES AFTER FOLDING WITH NEUT 4.5.1 CROSS SECTIONS
