/NextMassProduction

/TestMethodForNumu

/EventSpectra

/SensitivityStudies

/MarchStudies

Click on image to see full size figures.

1. 04b fluxes (Far/Near ratios) from JNUBEAM

These fluxes and Far/Near ratios were obtained using the official ntuples & macros from the JNUBEAM pages. Remember these are fluxes (no cross-section information yet).
Plot of (SK-C*2KM)/SK, where C is a constant to normalize the fluxes to the correct value. I chose C = (1800/295000)^2*1/(800*800) (1800m: distance of 2KM to the average neutrino production point ; distance to the source = 1839m). 1/(800*800) : convert the flux at 2KM to cm^-2. The same correction factor is used for the other F/N plots.

I produced the same plot using my own error bar computation (based on sum(weight^2) in each bin rather than the code in the official kumac with which I don't agree). The official kumac artificially underestimates the size of the error bars. (SK-2KM)/SK is compatible with zero above 1 GeV contrary to what the previous plot seems to indicate. We therefore cannot draw any conclusions about flatness of the F/N ratio above 1 GeV due to a lack of statitistics. The conclusion are unchanged between 0 and 1 GeV.

Same plot with different bin size to compensate for low statistics above 1 GeV.

Same plot, selecting only neutrinos with R<200cm (present 55.8t FV r = 175cm ; 100t FV r = 225cm). Note that selecting only "central" events reduces statistics in the first bin.

2. F/N with vector files

2.1. Explanation

The vector ntuples are the ntuples before simulation + reconstruction. They will therefore contain low charge events that don't make it through AFIT (and are not saved into the reconstructed ntuples).
They also contain Pauli blocked single pi events (modes 11-13, 31-33) ; after discussion with Hayato-san I use a cut he recommended two years ago.
Excerpt from his email :

cut 97 ((mode.le.10).or.((mode.ge.14).and.(mode.le.30)).or.(mode.ge.35))
cut 98 ((abspnu(2).eq.0).or.(abspnu(4).ge.0.217))
cut 99 (($97).or.($98))

 I checked the vector and there might be some problem
 I did not see any problematic events but potential problem
still exists. Therefore, please include those cuts in your
analysis for the safety.

explanation of the cut

keep events that are either :

• not single pi
• single pi on hydrogen (no Pauli blocking needed)

• single pi on oxygen with outgoing nucleon above the Fermi surface

Using the following criteria : vertex (=true vertex) inside the FV + Pauli blocked events removal, the following F/N ratio is obtained. I used the 55.8t FV and I now have the correct normalization obtained from Hayato-san. These plots were built using the same method (L_2km=1839m-39m) as the previous ones.

2.2. F/N and (SK-2KM)/SK at the vector level

Compare with F/N at the flux level.

COMMENTS : very low stats below 200 MeV or so. Since the 55.8t FV is in use, we select events close to the center of the tank, therefore we select an energy spectrum close to the one we expect at SK. We expect greater flatness.
The large peak at 1.5 GeV is under investigation. I suspect it comes from sampling errors : the flux is only imperfectly known (as explained in section 1). The vector ntuples are made using those events (or resampling according to the same histogram), but the statistical error at the flux level is not correctly propagated into the vector files. The errors are therefore underestimated at the vector level. See plots below.
The corresponding event spectra can be seen in the following picture :

2.3. Distribution of (SK-2KM)/SK

Using a very fine binning (200 bins from 0 to 1.5 GeV), it is possible to make the distribution of (SK-2KM)/SK (computed as above). There are therefore 200 entries in this histogram. The distribution should be gaussian. At 2KM I used the 56t FV (radius<175cm). The results are visible in the next plot :

2.4. Situation above 1 GeV + NEUT statistics

In the plots above it appears that there is a large fluctuation around 1-2 GeV, creating a large unexplained peak. The following plot shows equal bin histograms of the event spectra at 2KM & SK. It seems that :

• the error bars are underestimated above 1 GeV. I believe the output of JNUBEAM is piped into NEUT. Although I don't know the exact method, I can guess that having (very) limited statistics for the flux distribution introduces statistical errors that are propagated into NEUT. But they are not yet included in the error bars which appear to be too short.

• the fluctuations we see above 1 GeV seem to be statistical, see the plot in this section. I used equal bins even at higher energies. It appears that SK

has small fluctuations ( 225,000 events in the FV + much higher statistics fin NEUT). As expected the fluctuations are much higher at 2KM (66,000 events in FV56t and much higher errors from JNUBEAM). They could conspire to make this peak. The discrepancy between SK & 2KM seems to be smoothed out if we integrate all events between 1 and 5 GeV...

2.4.1. same events

Integrating above 1 GeV, and using the same events as above, the following plot is obtained :

The peak at 1.5 GeV is smoothed out. The difference between 2KM and SK is visible below in slightly different conditions that should not make any difference at high energies.

2.4.2. extra cut P,,mu,,>200 MeV/c and CC events only

An extra cut was added to make the plot below : I added mode < 30 (ie CC) and P_lepton > 200 MeV/c in order to remove low energy events that are lost before reconstruction. This is at best an approximation of the loss caused by the reconstruction, but it should not make any difference at high energies.

3. SKreweighted / 2KM : F/N with reconstructed, reweighted files

NOTE

All these plots are made using MC true neutrino energy. All the events generated inside the FV are used. There are no explicit reconstruction cuts but there could be implicit ones (that affect the low energy end of the spectrum)

The goal is to make the same plots as in the previous section, using the reconstructed ntuples. For SK we only have 03X reconstructed ntuples, that we reweight to 04b as explained below. At 2KM, the reconstructed ntuples are a subset of the ones in the previous section(vector-ntuples). Events for which no vertex can be found in AFIT (Need to check the conditions in details) are removed (any others ?).

3.1. Details of the reweighting procedure

WARNING : There is officially no SK04b. But 04a is the same as 04b as far as Super-Kamiokande is concerned.

3.1.1. Reweighting the SK ntuples to the new flux model

A word of explanation. This uses the MC true (ie vector-level) information in the reconstructed files.

in PAW pseudo-code :

1. nt/pro h2km //2KMntuples/1.pnu(1) $(generated in SK FV)
2. hsk is more complicated : nt/pro hsk //SKntuples/1.pnu(1) weight()*(generated in 2KM FV)

• where weight() = fluxSK04b(i)*xsec(i)/nevent(i)
  • with :
    1. xsec(i) : energy-dependent total nu_mu cross section
    2. i = int(pnu(1)/binwidth) (true neutrino energy bin) (+1 in fortran not in C)
    3. nevent(i) : number of events generated in FV with pnu(1) in bin i to get the correct normalization.

The point of this reweighting is to make the SK 03a files behave like 04b files. Our 2km events in the ntuples are sampled according to the 2KM flux04b*xsec distribution. The SK events are sampled according to SK03a*xsec, which we reweight to SK04b*xsec.

3.1.2. F/N & (SK-2KM)/2KM plots after reconstruction

In order to compare with the "theoretical" plots above, here's a plot of (SKreweighted-2KM)/2KM. It is possible he error bars are understimated.

The same trend (negative slope) as in the F/N plots obtained from the beam MC (flux level) is observed. however the first bins are way out of the Y axis range (statistics probably too limited to obtain anything meaningful).

Compare with plots at the flux level and at the vector level.

3.2. Event spectra

Using the same cut (ie numu events generated in the FV), the following event rates plots were obtained (see fig). SK was reweighted to 04b, assuming no oscillation takes place. 2KM was rescaled (not reweighted), with the following factor : 22500/55.8025*(1800.0/295000.0)²

It appears that the peak is shifted. Is there a difference in the cross-section model ? under study

same plot without the zoom :

3.3. Comparison between ''vector'' & ''reconstructed'' plots

Here reconstructed means made using MC true info from the reconstructed ntuples. Clearly at energies above roughly 200 MeV the plots made from 04b vector files and the plots made from reweighted SK03X and 2KM04b after reconstruction are very similar. This seems to show that the reweighted SK ntuples behave like 04b.
At the moment I have no explanation for the large differences at high energies between SK and 2KM.

4. Why are the spectra at 2KM & SK different

4.1. Top/Down effect

Effect of the finite size of the "detector" (defined as a 800cmx800cm square window in the direction of SK 1839 m away from the target). The beam points down in the Y direction (so neutrinos with y in (-400 cm; 0 cm) at 2KM are more on axis than neutrinos with y in (0 ; 400) ). The following plots illustrate this.

It's also interesting to bin the detector in the Y direction (vertical direction at the target site), and compute the average neutrino energy for each of these bins. The following plot illustrates this :

The difference is about 5.5% between the top and the bottom.

Same plot but using the peak energy. This is more difficult (in order to see the effect I need to reduce the bin size therefore the histogram becomes more sensitive to statistical fluctuations).

The effect is in fact 2-dimensional because the detector is also off axis in X. In JNUBEAM it is off axis in the direction of SK, so it is more off axis as X increases. In 'real-life' the 2KM will be off axis in the HK direction so we should see the exact opposite dependance in X.

Yet another way to look at this is to take 'iso-off-axis angle' curves. These should be circles (assuming the beam is a cone originating from a single point). In the following plot, the first calculation assumes that the beam source is at the target (0,0,0) ; in the right-hand side plot the origin is a point 40 m downstream along the beam. For 2.5 degrees, at a distance of ~2km, the radius is about 80m so the curves can be safely considered to be straight lines in the detector. The shapes are not changed when the origin is changed but the value of the off axis angle is (from 2.43 to 2.49 : the 2KM position is optimized to be 2.5 degrees off axis with the average production point). This plot suggests an other binning of all the previous plots to maximize the effect.

It is possible to compute the slope of those lines. To simplify let's assume that the beam source is at the target (and not at the average production point). This being said it is possible to compute the angle between the beam direction and any point, for example a point on the 2KM detector 'window', and therefore compute the iso off axis curve. It is an ellipse because this window is not orthogonal to the beam direction (it is vertical and the beam points 3.637 degrees downwards (gamma angle) according to JNUBEAM). The equation is : x² + cos²(gamma) *( y + tan(gamma)*Z)² = Z²*sin²(theta off axis). It is almost circular (would be if gamma were 0). Approximating the ellipse by its tangent, the slope of the iso off axis lines in the previous plots is -X/(Y+Z*tan(gamma)) = -0.33 (X,Y,Z are the coordinates of the 2KM center in JNUBEAM's master reference system).

4.2. Other possible effects

PRELIMINARY

There are three reasons (at least) that can account for differences in F & N spectrum (when N= 2KM).

• Finite size of the 2KM detector : even with perfectly focalized monochromatic pions decaying at the same point,

the spectrum at SK would only be identical to the spectrum at the center of the 2KM.

• The neutrino source is in fact not point-like, viewed from the 2KM. Different points of the decay tunnel correspond to different off-axis angles. Seen from SK the source is point-like, causing differences.
• The pion beam is not perfectly focalized (see below for an estimate of the angular divergence). Pions with different directions correspond to different off-axis angles.

4.3. Estimate of the pion beam divergence

4.3.1. A simple model

If components x&y (orthogonal to the beam direction) of the pion direction are gaussians with width sigma, then the pdf for the angle between the pion and the on-axis direction is theta*exp(-theta²/(2 sigma²)). Plots for these (pions only) can be seen below. If this equation were exact then the distribution of 1/2*(theta²) would be exponential and the plot on the right hand side would be a straight line. There is a deviation. UNDER INVESTIGATION.

There is a difference between the distributions found using the SK ntuples (in red) and the 2KM ntuples. This probably comes "acceptance" reasons (in JNUBEAM every pion that decays inside the tunnel is forced to decay to SK, ie a weight is assigned to the neutrino that would results if such a decay were found. For the NDs the situation is different : 1000 trial decays are generated and the resulting neutrinos are tested ; if they pass through the 2KM window they are saved ; obviously for a fraction of the pions no decays pass through the 2KM ; consequently less pions are found in the 2KM ntuples than the SK ntuples).

Using the mean of the right hand side plot @ 2KM the divergence (= sqrt(RMS) ) is estimated to be 20 mrad (still under investigation).

4.3.2. Effect on the off-axis angle

If the pion beam has a divergence of 20 mrad, what is the distribution of off-axis angles at 2KM ?

The following plot (obtained from a toy MC) is the answer (preliminary) :

This distribution is different at each point of the 2KM detector. Then it is interesting to see how large the effect is (is-it large enough to explain the variation of <E> or E_peak with y ?) for the top and bottom ( +/- 2 mrad up/down).

The results are (in degrees) :

• center point : mean = 2.60379 ; sigma = 1.04802
• bottom point : mean = 2.50799 ; sigma = 1.0373
• top point : mean = 2.70325 ; sigma = 1.05851

The peak energy is ~60MeV/(2*off-axis angle (rad)) ( 60 MeV is 2*E(nu) in the pion center of mass frame) leading to :

• centre : 660 MeV
• bottom : 675 MeV
• top : 636 MeV

Systematically ~20 MeV higher than in the E_peak vs y plot. Under study. But a variation pf the same order is observed.
It seems that the pion beam divergence is not enough to wash out this effect ?

5. Nu_mu analysis

Preliminary plots for numu analysis at 2KM. I 'removed' the extra events added to compensate for the lack of stats since they have different direction PDFs (see section Other Topics).

The selected events are : FV (55.8t), FC, 1 ring, mulike, evis > 30.0 MeV , Pmu > 200.0 MeV

Note that :

• FC = pomax < 100.0 @ 2KM

• FC = nhitac< 10 @ SK

I neglected mis-ID nue events in this study since this effect is known to be extremely small (ie only numu ntuples were used).

2KM was scaled to SK using the same factor as nues : M_sk/M_2km*(L_2km/L_SK)² = 22.5kt/55.8t*(1800.0/295000.0)² (I chose L2km=1800.0 m meaning the average production point is estimated to be 39 m after the target).

I used the 'standard' numu energy reconstruction equation, ie I did not correct for the nuclear potential (-27 MeV).

• Fully reconstructed spectra and ratio

• same events as above, but using the true neutrino energy

Same plots for CCQE events only. This time the analysis cuts are : FV (55.8t), FC, evis > 30.0, Pmu>200.0, 'mode==1' (ie CCQE)

• CCQE, reconstructed neutrino energy

• CCQE, true neutrino energy

6. Muon momentum reconstruction study

scatter plots : reconstructed mu momentum (amomm[0]) vs true momentum (ideally this should be pnu(3) but G4 does not sort particles in that order so it is necessary to search for the outgoing muon by hand).

All the muon analysis cuts were used (FC,FV,1r,mu,evis>30) except Pmu>200.0 . 2KM is the left panel (green) and SK is the right panel (red).

Same plot for CCQE events only : the reconstruction performances seem to be similar in both detectors for CCQE (but the acceptances are different).

7. Preparing the next 2KM mass production

7.1. New digitizer in G4

New features :

• qisk *= G4RandGaus::shoot(1.0,0.05) --> new charge smearing

• add a trigger t0 : smears the starting point of the simulation (same as dshigh.F) ; now tisk = firsthit + t0, not firstHit + constant --> broadens timing dsitribution

• add integration gates (300ns) and an event gate (-t0+offset+950ns). qisk = sum(hits with timing < min(firsthit+300ns,-t0+offset+950ns)).

This will remove delayed Cherenkov light coming from decay electrons and hadronic interactions. Note that it doesn't affect through-going muons very much because hardly any hits fall outside the gate for such events (no delayed light). Still, it appeared that it was necessary to smear the charge by an extra 5% to match data.

The main flaws are still present :

• reflection peak (small second peak) still too large
• still too much indirect light (charge ratio)

Note for people working on the digitizer in the future :

• The efficiency coefficient was found by comparing the charge scale with data (previous plot). It is adjusted by trial and error...
  

• All Cherenkov light raw hits are present in the event ( 'time' member in JHF2kmWCHit) (However raw Ch. hits are not saved in the ROOT files to save disk space). With the new integration gate they are simply no longer digitized. It is possible to check whether there are delayed hits by using the 'triggerhisto' buffer that is now filled up in JHF2kmWCDigitizer::Digitize(). Then it's up to you to find a smart way of saving them that reproduces data...

7.2. Tuning to through-going muons

Warning

• This is a work in progress. It uses unphysical scattering parameters that probably contradict laser data. Note that the prevoious situation used imperfect tuning that contradicted reflections in cosmic ray data...

• In Red : K2K 1kton data
• In Black : G4

The Rayleigh scattering lengths were divided by 2 (ie Rayleigh scattering is increased). Reflections off PMTs are now described differently (PMT glass index is 1.55 and the 'cathode' is polished). The idea is to reduce the amount of reflected light. We get better agreement in terms of timing. The shape of the charge ratio plots is still different but we have less indirect light this time. This could probably be refined but we need to think about the strategy :

• we probably need to implement scattering off large centers (Mie scattering)

• we may need to change the angular distributions of scattering off small centers (Rayleigh or Einstein-Smoluchowski), and the power law for the corresponing scattering lengths.

See the Handbook of Optics for interesting details.

At the moment we use unphysical scattering parameters...

8. Other topics

8.1. Extra events added to increase statistics in the present analysis

• In black : events with vertices inside the cut
• In red : events that were moved inside the vertex area. Disrupts the beam direction.

The explanation goes as follows : these extra events are taken from the outer part of the detector window. Therefore they point in slightly different directions (depending on which quadrant of the detector they are in). The distributions of their direction components have the same shapes, but the means are shifted. The red curves correspond to a sum of these shitfed distributions.
Notice the neutrino energy is also shifted.

9. PID comparisons

9.1. tests with spexppe

Following the talk I gave at the 2KM videomeeting, here are some pattern comparisons between G4 and spexppe (used to compute expected light for PID purposes).

I used about 150 stable muons, near the center of the tank (25t cylinder like the standard FV but centered). The left-hand side plot show the standard result, the right-hand side plot show what happens when the CONV factor in 'spperini' is increased by a factor of 10. Please send me your comments. Note that only PMTs in a 70 degree cone are used so the agreement at very large angles is not that important.

With the "PMT use flag" turned on :

• 

9.2. tests with sppemse + comparison with kton 02a muon MC

Using a 1yr-old muon file from the kton 02a MC (with muon decay turned off), and an up-to-date G4 'stable' muon file :

• red and dashed red : computation from sppemse (includes scattering, reflections, knock-on e-), using the events from the kton file and the G4 file respectively. Should be and actually are similar.

• black : K2K 1kton MC 02a.

• blue : geant4

Interesting features : in both MC's the inside of the ring seems to be more filled than the computation. Am I forgetting anything in the initialization ?. The ring also seems to be less sharp in G4 than in ktdetsim...

Same plots for e- : left-hand side is 1kton, right hand side is G4. The red histos are expected profiles computed with sppemse.

9.3. comparison between G4 and 'old' kton 02a -- RESCALED PLOTS

Left hand plot::

• e- profile (500 MeV/c) ; blue = geant4, direct light only ; black = G3 (all light, 02a, but ARAS=0.6 instead of 1.5)

Right hand plot::

• mu- profiles (500 Mev/c) ; see the legend.

9.4. SUMMARY PLOTS-- Using SK2 routines

Using Okumura-san's SK2 routines, and without any rescaling of the plots. Electrons on the left, muons on the right.

9.5. Attempt

I also took profiles using the tabulated charge pattern engine for muons called 'pfmuonth' , which is not used in the standard PID/MSFIT code I replaced all the calculations in '*mkexp' and 'sppemse' with the appropriate calls to pfelecth and pfmuonth (mimicking the call to pfelecth that is already present in pfmkexp.F), and re-ran on the latest K2K 1kton numu MC files I have (with increased scattering). As expected we get PID peaks away further away from zero than the data. STILL UNDER STUDY. At the moment I am retuning G4 using the SK 1pe probablility function. Hopefully the patterns will agree with data better, especially in the parts with low pe hits.

||