Spline Validation Procedure

Full proccedure for how to check and validate splines.

Plot Event Rate Spectra

All Systematic Parameters

This is good to get an overall picture before delving into details.

Can either use all systematics or split into catagories (flux, int, modern int etc.).

1. Write an analysis/config file that contains the systematics you are interested in (for a specific detector).

2. Obtain tweaks file from VALOR:

valor-tweak_pdf_N \
  --config analysis/config_params.xml \
  --output tweaks_valor_params.root \
  --ntoys 1000 \
  --physics-hypothesis void

3. Obtain tweaks from universes

    valor_sbn-universe_variation_tweaks \
      --config analysis/config_params.xml \
      --input /path_to_dst_with_universes.root \
      --params params \
      --ntoys 500 \
      --output tweaks_uni_params.root \
      --sample numuCCincl

see Apps/universe_variation_tweaks.cxx for options for params.

I used only 500 toys here as I got an error with some of my configurations, will look into this later but realistically anything above a few hundred is sufficient.

4. (Optional) Plot overlay

python $VALOR_EXPT/macros/mock_data_challenge/draw_postfit_data_error_envelope.py \
  -c analysis/config_params.xml \
  -p tweaks_uni_params.root \
  -b tweaks_valor_params.root \
  -f root \
  -t tag \
  -e 1 \ // Experiment enum, SBND = 1 etc
  --plabel 'SBND universes #pm 1 #sigma' \
  --blabel "SBND VALOR #pm 1 #sigma"

Can remove -b/-blabel to only plot one.

sbnprism branch version looks nice (from Rhiannon), allows addition of -a and -alabel to include, for example, mock data.

5. Produce ratio plots

valor -b -l -q $VALOR_EXPT/macros/systematic_construction/plot_envelope_ratios.C'("tweaks_valor_params.root","tweaks_uni_params.root","output_tag")'

Output files will have form universe_valor_*_ratios_<output_tag>.root.

Looking at double ratios (using SBND as example):

• If using all params: ratio will probably up to the order of 5/6%
• If using flux: ratio will probably on the order of 1-2%
• If using proposal interaction: ratio will probably on the order of <1%
• If using modern interaction: ratio will probably on the order of 1-4%

Individual Systematic Parameters

Look at plots showing the event rate spectra from both valor and universes for 0, ±1, ±2, ±3 σ.

Best to use loops to quickly plot all.

Example for flux params in SBND:

params=("horncurrent_FluxUnisim" "expskin_FluxUnisim" "nucleoninexsec_FluxUnisim" "nucleontotxsec_FluxUnisim" "nucleonqexsec_FluxUnisim" "pioninexsec_FluxUnisim" "piontotxsec_FluxUnisim" "pionqexsec_FluxUnisim")

for j in ${!params[@]}; do
  p=${params[$j]}

  valor_sbn-tweak_pdf_nsigma_var \
  --config analysis/config.xml \
  --params $p \
  --physics-hypothesis void \
  --universe-spectra /path/to/spectra/file/valor_splines_numuCCincl_${p}_BNB_FHC_sbn_nd_spectra.root \
  --output nsigma_var_${p}_numuCCincl_all_sigma.root

  valor -l $VALOR_EXPT/macros/systematic_construction/spline_validation.C'("'${p}'","numuCCincl","BNB_FHC","sbn_nd","/output/files/location/","/spline/location/",true,1,false,"nsigma_var")'

done

It's a pain but you now have ~350 images to look at:

• Nominals: should have 0 variation universe from nominal and VALOR from nominal (0 variation universe from VALOR). We set the nominals explicitly to be equal so if these are different then theres a bug somewhere.
• ±Xσ: Total universe variation from VALOR should be 0 if perfect. May see some be 1 or 2 due to precision errors in subtraction. At least the order of the ΔRate/Rate(Average),% should only be a few percent.

Debugging

If either of the above sets of plots are not right, can use the following to debug

• Nominal: there are 2 types of nominal defined, true nominal at 0σ and a pseudo nominal which is the closest universe to 0σ.

  • Check that both use the true nominal in spline_validation.C.
  • Check number of events in standard and universe DSTs agree.

• Find a non-flat spline to investigate:

  • Plot all splines using view_splines.cc. Adapt to your setup accordingly:
    • Use a systematic that looked particularly dodgy
    • Loop limits for numu: true = 22, reco = 21 (bins numbered from 0)
    • Loop limits for nue: true = 33, reco = 12
  • Look at outputs and pick one that's not flat, make note of bin number (r, t)

• Get MCT bin event rate from VALOR:

  • Add --def-spline-debug, --reco-bin r, --true-bin t to arguments of tweak_pdf_nsigma_var.cxx (what we used to make the tweaked event rate spectras for all sigmas for individual parameters earlier)
  • 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.
  • Will have both pseudo nominal and true nominal (true is at the very end of output)

• Get MCT bin event rate from universes:

  • Add --def-spline-debug, --reco-bin r, --true-bin t to arguments of valor_sbn-make_flux_splines (or valor_sbn-make_xsec_splines if picked an interaction systematic)
  • 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

• Check that the Total event rate for reco bin r: X match (note that the universe ones are printed every 0.5 sigma and the VALOR ones are every whole sigma.)

• Manually construct Spline:

  • The spline maker will give us the variations from the universes in question and the corresponding event rate.
    • Divide each of the event rates (y-values) by the nominal value
    • Generate the TSpline3 (x-values = universe variations, y-values = scaled event rate)

• Load spline in root:

  • root valor_splines_numuCCincl_expskin_FluxUnisim_BNB_FHC_sbn_nd.root
  • TSpline3 *my_spline = (TSpline3*)spl_numuCCQE_expskin_FluxUnisim->At(22)
    • The 22 is reco-true bin (1,1) for numu
    • Calculate by doing: bin = (t_max*t) + r
    • For this example (21*1)+1 = 22 (spline bins numbered from 0, limits above are for a loop that is i < n_t_max so have to do n_t_max=1+t_max)

• Check spline:

  • Is my_spline->Draw() correct?
  • Does my_spline->Eval(0) = 1.0000000 ?

• Compare to manual spline using universe variation

  • True nominal MCT bin ER * my_spline->Eval(V) = MCT bin event rate from the universes for that variation
    • True nominal from end of output of valor_sbn-tweak_pdf_nsigma_var
    • ER from universes given under Looking at <param> with sigma variation = <V>, Printing MCT event rate breakdown for reco bin = r and true bin = t:

If all of the above works well, look at:

• 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.
  • Look at tick table: https://github.com/VALOR-fit/VALOR-SBN/wiki/Produce-systematic-parameter-tick-tables
    • Green tick: mode in modelist, param applied to beam, sample, detector, mode and nuclear target (Ar for SBND)
    • Yellow cross: mode not in modelist, param applied to ...
    • Red tick: mode in modelist, param not applied to ...
    • Grey: mode not in modelist, param not applied to ...
  • Can temporarily set reac_type to 'all' to check. This should use all the splines available, but I guess is bad practice.
• Are the binning definitions given in analysis/congig.xml and mc/config.xml the same?

After you've double checked all the above and it's still broken, sorry I can't help you anymore :(

Adapted from:

• https://github.com/VALOR-fit/VALOR-SBN/wiki/Produce-and-validate-event-rate-spectra-with-single-systematic-variations
• https://github.com/VALOR-fit/VALOR-SBN/wiki/Produce-error-envelopes-for-event-rate-spectra
• https://github.com/VALOR-fit/VALOR-SBN/wiki/Debugging-response-functions