Debugging Response Functions

Some things to consider if the results from Produce and validate event rate spectra with single systematic variations and/or Produce error envelopes for event rate spectra don't look good:

Generating the response functions 'correctly' (for nue)

• For nue, we need to generate independent sets of response functions for each sample.

  • Should do this for cosmic and dirt too? The cosmic sample exists in name only atm (after all the cuts we end up with no events despite having initially generated millions..). Dirt events seem to get mapped to some numu reac modes and not the dirt mode - something to look at in the future? For now have lumped cosmic and dirt samples with numu.
  • We'd expect the response functions from the different samples but for the same reac mode and bin to look different (e.g. nueCCQE, numuCCQE and numutonueCCQE would not be identical).

• When generating the response functions from the oscillated sample (for nue), we need to make sure we're actually reading some events. If the PhysHypothesis is set to kHypVoid we have no oscillations so no events to generate the splines.

  • Ideally want to take all events in the osc sample (with weight = 1, without applying oscillation weighting), but not straight forward to do with the existing apps. Easiest work around is just to set some osc params so we get some events (pick sinsq = 1 and dmsq = 100 to avoid dealing with very small or vanishing numbers in some Etrue bins). Oscillations and systematic effects should “commute”, so I think this approach should be okay.

Choices + things to tweak/change and be aware of

• In the spline construction app, there is a line similar to

  double event_ratio = (nevents0 > 0) ? nevents/double(nevents0) : 1; // Ratio of the variation rate to the nominal rate

• This requires more than N events (0 in this case) in a MCT bin before producing a non-flat spline. In previous iterations (for the numu analysis in particular), the check was set to nevents0 > 4. Not a 100% sure on the reasoning for choosing a non-zero number. Picking 0 should generate a spline which is a 'perfect' representation, but perhaps if splines look strange it might be caused by dividing by a tiny number.

• Generating the splines is done by comparing the nominal to nsigma variations (see line above). Previously (again for the numu analysis), the nominal was produced by using the universe with a variation closest to zero - have now changed to using the true nominal (i.e. no variations at all).

• Need to be careful as we can have a 'true' nominal spectrum (labelled 'nom') and a 'pseudo' nominal spectrum which is produced by using the universe variation closest to 0 (labelled 'p0'). For example the tweak_pdf_nsigma_var.cxx app generates both and when using the spline_validation.C macro to compare the spectra from the universes and from VALOR, we need to make sure we're using the same nominal spectrum.

• The universe_variation_tweaks.cxx app counts events differently to the other apps/macros. It doesn't care about any PhysHypothesis and does what is mentioned above in the part about generating response functions for the osc sample I believe. This means that the spectrum from the oscillated sample will be different and the nominal ratio isn't comparable to the others. The single and double ratios are of course still comparable. The intrinsic and numu samples are unaffected by this.

Deep dive into event rates

• Probably best to focus on an individual systematic and a single sample (if looking at the nue analysis) i.e. look at the plots from Produce and validate event rate spectra with single systematic variations.

• First check the nominal spectra agree. If not confirm that the standard and universe dsts have the same number of events and then check we're looking at the same nominal spectra (see above)

• If we've made it this far, it must be the case that either the splines are wrong and/or they are not being applied. I'm not aware of an intuitive way to know which splines we'd expect to be non-flat so lets just print them all out and look at them (See for example /VALOR-SBN/macros/contrib/tom/view_splines.cc and adapt it to your setup). Pick a non-flat spline - we now have a MCT bin to investigate further.

• Enable the debug option in the spline constructor (add --def-spline-debug and your choice of MCT bin --reco-bin X and --true-bin Y to the command line args). This will print out the MCT event rate for each mode and for each of the sigma variations for that bin + the total event rate in that reco bin (https://github.com/VALOR-fit/VALOR-SBN/blob/2a2e84fea4474cfef26f1499e8764f9efe29d9b0/src/sbn/Apps/make_xsec_splines.cxx#L425).

• For a given reco bin, the sum of all the associated MCT bins should give event rate in that reco bin seen in the 'standard' spectra.

• If we want to check the spline has been constructed correctly, we can manually construct the same spline and see how they compare. The spline maker will give us the variations from the universes in question and the corresponding event rate. Just need to divide each of the event rates (y-values) by the nominal value and then generate the TSpline3 (x-values = universe variations, y-values = scaled event rate). The main causes of any issues here I think would be due to;

  • Confusing the 'true' and 'pseudo' nominal as discussed above.
  • Requiring some sort of minimum event rate before constructing the splines.

• Print out the MCT bin event rate from VALOR using the splines (add --def-spline-debug and your choice of MCT bin --reco-bin X and --true-bin Y to the command line args for the tweak_pdf_nsigma_var.cxx app). Compare this to the output from the universes given by the spline maker. Note, the plots produced from this app are labeled as some integer sigma value, but are in fact using the actual variation value from the given universe so should match the output from the spline constructor app.

• Load up a spline of interest in root. e.g. Lets say I want the following spline spl_BNB_FHC_nuelikeCChigh_numutonueCCQE_horncurrent_FluxUnisim_32_7, then we can do something ala

  root valor_splines_nuelikeCChigh_horncurrent_FluxUnisim_BNB_FHC_sbn_nd.root
  TSpline3 *my_spline = (TSpline3*)spl_numutonueCCQE_horncurrent_FluxUnisim->At(391);
  my_spline->Draw();  \\ check this is indeed the spline we want
  my_spline->Eval(some variation);
  

• In this example we have 12 reco bins and 33 true bins. So to get the correct spline from the TObjArray, we do (32 * 12) + 7 = 391.

• As a check, my_spline->Eval(0) should always evaluate to 1.

• Now we can pick one of the universe variations, V. Then the true nominal MCT bin event rate * my_spline->Eval(V) should equal the MCT bin event rate from the universes for that variation.

• If splines exist and are correct, but valor and universes still don't agree, the first thing to check is the reac_type in the definition of the systematic. Probably what's happening is that some modes are having an impact, but we've said that they shouldn't so the corresponding splines are being ignored. Can temporarily set reac_type to 'all' to check. This should use all the splines available, but I guess is bad practice.

• Results from Produce and validate event rate spectra with single systematic variations can (and maybe should) agree perfectly down to the event level since the response function we generate are a perfect representation of the variations (subject to some caveats mentioned earlier). Plots from Produce error envelopes for event rate spectra are produced by taking the 1 sigma variation from N throws, so since we're taking an average, I don't think it's realistic to get a perfect agreement (but of course results should still be 'similar').