Creating a response matrix (effective aperture)

Introduction

For some programs or tasks a pre-computed response matrix including the effective aperture depending on stellar position is advantageous. Nyx (differentiable rendering of night sky background) makes heavy use of them. We implemented a convenience function to render these directly.

import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt

from iactrace import MCIntegrator, Telescope

Preparation

We first import the HESS_I telescope as an example (we only sample 1 ray per facet in this tutorial, for better accuracy increase the value but watch out for memory consumption):

telescope = Telescope.from_yaml('../../configs/HESS/CT3.yaml', MCIntegrator(1), key = jax.random.key(42))

Then we apply mirror roughness:

telescope = telescope.apply_roughness(24)

Generating response matrix for all pixels for one telescope configuration

We can then create a grid of stellar sources (with unit flux) within our FoV:

fov_deg = 5.5
fov_rad = fov_deg * jnp.pi / 180

n_side = 512

x1d = jnp.linspace(-fov_rad / 2, fov_rad / 2, n_side)
y1d = jnp.linspace(-fov_rad / 2, fov_rad / 2, n_side)

X, Y = jnp.meshgrid(x1d, y1d, indexing="xy")

x = X.ravel()
y = Y.ravel()
z = -jnp.ones(n_side**2)

test_sources = jnp.stack([x, y, z], axis=1)
test_sources = test_sources / jnp.linalg.norm(test_sources, axis=1, keepdims=True)

Then we import and use the convenience function to render a response matrix for the hexagonal sensor:

from iactrace.core import render_response_matrix

r_m = render_response_matrix(telescope, test_sources, jnp.ones(test_sources.shape[0]),
                             'parallel', sensor_idx=0)

We can visualize the response matrix from the view of the stars (getting the total effective aperture for each star position):

fig, ax = plt.subplots(figsize=(15, 12))

weights = r_m.sum(axis=-1).reshape(n_side, n_side)

im = ax.imshow(
    weights,
    origin="lower",
    extent=[x1d[0], x1d[-1], y1d[0], y1d[-1]],
)

cbar = fig.colorbar(im, ax=ax)
cbar.set_label("Effective aperture [m^2]")  # colorbar label

ax.set_xlabel("x [rad]")
ax.set_ylabel("y [rad]")

plt.show()

Alternatively we can look at the per-pixel effective aperture:

fig, ax = plt.subplots(figsize=(15, 12))

weights = r_m[:,300].reshape(n_side, n_side)

im = ax.imshow(
    weights,
    origin="lower",
    extent=[x1d[0], x1d[-1], y1d[0], y1d[-1]],
)

cbar = fig.colorbar(im, ax=ax)
cbar.set_label("Effective aperture [m^2]")  # colorbar label

ax.set_xlabel("x [rad]")
ax.set_ylabel("y [rad]")

plt.show()

Transforming to per-pixel effective aperture array

For further use, separating the pixel responses into individual arrays of the same shape is advantageous. For this we choose 64x64 windows, to make these better for acceleration structures. We search the center of the box where values > 0 and then choose the 64x64 windows center on this row and column:

weights_3d = r_m.reshape(n_side, n_side, 960).transpose(2, 0, 1)
mask = weights_3d > 0

rows = jnp.arange(n_side)[None, :, None]
cols = jnp.arange(n_side)[None, None, :]

row_min = jnp.where(mask, rows, n_side).min(axis=(1, 2))
row_max = jnp.where(mask, rows, -1).max(axis=(1, 2))
col_min = jnp.where(mask, cols, n_side).min(axis=(1, 2))
col_max = jnp.where(mask, cols, -1).max(axis=(1, 2))

center_row = (row_min + row_max) // 2
center_col = (col_min + col_max) // 2

tops = jnp.clip(center_row - 32, 0, n_side - 64)
lefts = jnp.clip(center_col - 32, 0, n_side - 64)

empty = ~mask.any(axis=(1, 2))
tops = jnp.where(empty, 0, tops)
lefts = jnp.where(empty, 0, lefts)

We can check how these individual pixel apertures look like:

fig, ax = plt.subplots(figsize=(15, 12))

id = 200

top = tops[id]
left = lefts[id]

im = ax.imshow(
    weights_3d[id][top:top+64,left:left+64],
    origin="lower",
    extent=[x1d[left:left+64][0], x1d[left:left+64][-1], y1d[top:top+64][0], y1d[top:top+64][-1]]
)

cbar = fig.colorbar(im, ax=ax)
cbar.set_label("Effective aperture [m^2]")

plt.show()

Switching to random mirror alignment + displacement:

In IACTs, we often also have to model random mirror misalignment and displacements. For this there are convenience functions implemented:

# Apply misalignment
new_tel = telescope.apply_misalignment_to_group(0, 15, 10, jax.random.key(0))
# Apply displacement
new_tel = new_tel.apply_displacement_to_group(0, 0.02, jax.random.key(1))
# Apply focal error
new_tel = new_tel.apply_focal_error_to_group(0, 0.07, jax.random.key(2))
# Render response matrix
r_m_new = render_response_matrix(new_tel, test_sources, jnp.ones(test_sources.shape[0]), 'parallel', sensor_idx=0).block_until_ready()

fig, ax = plt.subplots(figsize=(15, 12))

weights_new = r_m_new.sum(axis=-1).reshape(n_side, n_side)

im = ax.imshow(
    weights_new,
    origin="lower",
    extent=[x1d[0], x1d[-1], y1d[0], y1d[-1]]
)

cbar = fig.colorbar(im, ax=ax)
cbar.set_label("Effective aperture [m^2]")

ax.set_xlabel("x [rad]")
ax.set_ylabel("y [rad]")

plt.show()