Implement simulated histogram shower-distribution event weighting

[kosack]

Please describe the use case that requires this feature.

For upcoming CTAO simulation campaigns, they will start to test importance-sampling, i.e, generating non-powerlaw shower energy distributions and also non-flat spatial distributions. Our current event weighting code (and pyIRF itself) only support a simple power-law weighting function read from the SimulationConfig. We do however, have the full histogram and propegate it (/simulation/service/shower_distribution), and this can be used as the correct weighting function.

Importance sampled simulations are expected in ≈4-6 months.

Describe the solution you'd like

• Allow 2D weighting functions in the IRF/etc codes: energy (already there for the 1D case), offset from array center (needs to be added)
• Create a HistogramShowerDistribution function that is a 2D interpolator over these parameters, using the shower_distribution.

Some preliminary implementation is in #2927 in the spectrum_from_simulation_config function, but only a stub where we would put this code. Some modification is needed in pyirf if we want to have the weighting function there, or we could keep pyirf simple and create a local weighting function in ctapipe.

[maxnoe]

@orelgueta what importance sampling are you going to use?

Only the one resulting in non unity area weight? (#1577)

Or something else?

[kosack]

For context, here is what is in the simulated histogram for a merged DL2 file:

Which if you plot area-corrected is just flat

But this supports at least non-uniform radial distributions. Is there an importance sampling version that does it in 2D area? If so, we can either make this a 3D hist, or just store event weights as in the issue @maxnoe linked above.

Code

from ctapipe.io import read_table
from astropy.visualization import quantity_support
from matplotlib import pyplot as plt
from matplotlib.colors import PowerNorm
quantity_support()

table = read_table(
    "gamma_20deg_180deg_run002012___cta-prod6-2156m-LaPalma-dark+magic_cone10.01_alpha_test_applied.DL2.h5",
    "/simulation/service/shower_distribution",
)

E = table["bins_energy"].quantity[0]
D = table["bins_core_dist"].quantity[0]
hist = table["histogram"].sum(axis=0)

E_cen = 0.5 * (E[1:] + E[:-1])
D_cen = 0.5 * (D[1:] + D[:-1])

plt.pcolormesh(E_cen, D_cen, hist.T, norm="log")
plt.xscale("log")
plt.title("Simulated showers")
plt.colorbar()

[orelgueta]

At the moment the tests are only for the non-unity area weights. @mariakherlakian has only used power law distributions so far. She looked into simulating with a harder spectrum, which leads to better trigger efficiency, so we will most likely do that, but I guess the event weights mentioned in #1577 should work with that as well.

If in the future we consider other energy distributions, we will revisit this point.

[maxnoe]

Ok thanks. In this case we only need the area weights.