iactsim 0.1.0

Imaging Atmospheric Cherenkov Telescopes simulation

iactsim: Imaging Atmospheric Cherenkov Telescope Simulation

Overview

iactsim is a Python package designed for simulating the response of Imaging Atmospheric Cherenkov Telescopes (IACTs). It exploits the computational power of GPUs to accelerate computationally intensive tasks such as ray-tracing and SiPM response simulation. This project aims to provide a fast, flexible, and user-friendly simulation framework for IACT performance studies, instrument design and data analysis within the Python ecosystem.

• Software documentation
• Source code

Features

• GPU-accelerated ray-tracing
• GPU-accelerated camera response
• Runtime configuration
• And more...

Installation

Prerequisites

• Python: 3.10, 3.11, or 3.12.

• NVIDIA HPC SDK: iactsim requires the NVIDIA HPC Software Development Kit (SDK) to compile CUDA kernels.

• CuPy: iactsim utilizes CuPy for GPU offloading.

Installation Instructions (dev)

1. Clone the repository:

   git clone https://gitlab.com/davide.mollica/iactsim.git
   

2. It is recommended to use a virtual environment. For example, using mamba:

   mamba create -n simenv python=3.12
   mamba activate simenv
   

3. Configure the enviroment to use HPC SDK compilers (you can download HPC SDK from the NVIDIA website). We suggest to use Environment Modules to handle SDK configuration and then define the NVCC and CUDA_PATH enviromental variables:

   module load nvhpc
   export NVCC=$NVHPC_ROOT/compilers/bin/nvcc
   export CUDA_PATH=$NVHPC_ROOT/cuda
   

   With conda/mamba enviroments you can use the provided configuration script configure_conda_env.sh

   mamba activate simenv
   ./configure_conda_env.sh simenv
   mamba deactivate
   

   This adds an activation script and a deactivation script to the simenv enviroment that will automatically handle the configuration when it is activated or deactivated.

4. Install CuPy:

   pip install cupy-cuda<XXX>
   

   Replace with your CUDA version (e.g., cupy-cuda12x). For more detailed instructions on installing CuPy, refer to the CuPy documentation.

5. Move inside the cloned folder and install iactsim from source:

   cd iactsim
   python -m pip install . -v
   

   For the C++ part, by deafult zlib-ng (zlib data compression library for the next generation systems) will be used (cmake will clone the repo automatically). If you have to use zlib append -C cmake.args="-DUSE_ZLIBNG=OFF" to the install command.

Conda virtual environments and Jupyter

If you want to run Jupyter Notebook or JupyterLab from a different environment (usually the base environment) and access the new one, do the following:

1. Install nb_conda_kernels in the base environment:

   mamba deactivate
   mamba install nb_conda_kernels 
   

2. Make sure ipykernel is installed in the new virtual environment

   mamba install -n simenv ipykernel
   

Now you can find the new environment kernel in JupyterLab/Notebook kernels list (Kernel->Change kernel in the menu bar).

In order to let the tqdm progress bar work properly, install ipywidgets following the installation guide.

Usage

Optical system definition

import iactsim
import matplotlib.pyplot as plt

plt.style.use('iactsim.iactsim')

# Spherical mirror
mirror_curvature_radius = 20000
plate_scale = mirror_curvature_radius/57.296/2.
mirror = iactsim.optics.SphericalSurface(
    half_aperture=10000., 
    curvature=1./mirror_curvature_radius,
    position=(0,0,0),
    surface_type=iactsim.optics.SurfaceType.REFLECTIVE_OUT, # reflective in the pointing direction
    name = 'Mirror'
)

# Flat focal surface (5deg hexagon)
focal_plane = iactsim.optics.FlatSurface(
    half_aperture = 5*plate_scale, 
    position = (0,0,0.5*mirror_curvature_radius),
    aperture_shape = iactsim.optics.ApertureShape.HEXAGONAL,
    surface_type=iactsim.optics.SurfaceType.SENSITIVE_IN, # sensitive surface opposite to the pointing direction
    name = 'Focal Plane'
)

# Optical system
os = iactsim.optics.OpticalSystem(surfaces=[focal_plane, mirror], name='TEST-OS')

# Telescope position
pos = (0,0,0)

# Telescope pointing (alt,az)
poi = (0.,0.)

# IACT
telescope = iactsim.IACT(os, position=pos, pointing=poi)

# Copy data to the device
telescope.cuda_init()

# Photon source initialized on-axis 
source = iactsim.optics.sources.Source(telescope)
source.positions.radial_uniformity = False
source.positions.random = False

# Plot spot diagram at different off-axis angles
n_plots = 5
fig, axes = plt.subplots(1,n_plots,figsize=(3.5*n_plots,3.5))

source.directions.altitude -= 2.
for ax in axes:
    # Adjust photon position to match the mirror position
    source.set_target('Mirror')
    
    # Generate photons
    ps, vs, wls, ts = source.generate(10000)
    
    # Perform ray-tracing
    telescope.trace_photons(ps, vs, wls, ts)
    
    # Plot spot diagram
    iactsim.visualization.scatter(ps, s=0.2, ax=ax, color='black', alpha=0.5, scale=plate_scale)
    ax.set_xlabel('X (deg)')
    ax.set_ylabel('Y (deg)')
    ax.grid(ls='--')

    # Move the source
    source.directions.altitude += 1. # degree

plt.tight_layout()
plt.show()

Mirror segmentation

The following code provides an example of how to segment a surface (AsphericalSurface, SphericalSurface or FlatSurface) starting from a mother surface (in this case mirror). Note that each segment is an independent surface and does not need a mother surface, which is used here simply for convenience.

import numpy as np

# List of segments
segments = []

# Segment ID
k = 0

# Segments on a 10X10 grid, 80 total
n = 10

segment_distance = 2*mirror.half_aperture / (n+3)

for i in range(n+3):
    for j in range(n+3):
        offset = [-mirror.half_aperture+segment_distance*i, -mirror.half_aperture+segment_distance*j]
        r_segment = np.sqrt(offset[0]**2+offset[1]**2)
        
        # Do not create segments outside the original mirror aperture
        if r_segment > mirror.half_aperture-segment_distance*np.sqrt(2):
            continue

        # Ideal segment position
        segment_position = [
            offset[0],
            offset[1],
            mirror.sagitta(r_segment),
        ]
        
        # Create the surface
        segment = iactsim.optics.SphericalSurface(
            curvature=mirror.curvature,
            half_aperture=0.45*segment_distance, 
            position=segment_position,
            surface_type=mirror.type,
            name = f'Segment-{k}',
            aperture_shape=iactsim.optics.ApertureShape.SQUARE,
            tilt_angles=np.random.normal(0,1,3), # Big random dispersion
            scattering_dispersion=0.05
        )
        
        # Specify the segment offset
        # When a segment is created in this way:
        #  - it will be oriented with the same surface normal 
        #    of the mother surface at the specified offset 
        #  - `tilt_angles` attribute will define a deviation from this orientation.
        segment.offset = offset
        
        segments.append(segment)
        k += 1

# Optical system
segmented_os = iactsim.optics.OpticalSystem(surfaces=[focal_plane, *segments], name='SEGMENTED-TEST-OS')

# IACT
segmented_telescope = iactsim.IACT(segmented_os, position=pos, pointing=poi)
segmented_telescope.cuda_init()

# Plot spot diagram at different off-axis angles
n_plots = 5
fig, axes = plt.subplots(1,n_plots,figsize=(3.5*n_plots,3.5))

# Photon source initialized on-axis 
source = iactsim.optics.sources.Source(segmented_telescope)
source.positions.radial_uniformity = False
source.positions.random = False
source.positions.r_max = mirror.half_aperture*1.5

source.directions.altitude -= 2.
for ax in axes:
    # Adjust photon position to match the mirror position
    source.set_target()
    
    # Generate photons
    ps, vs, wls, ts = source.generate(10000)
    
    # Perform ray-tracing
    segmented_telescope.trace_photons(ps, vs, wls, ts)
    
    # Plot spot diagram
    iactsim.visualization.scatter(ps, s=0.2, ax=ax, color='black', alpha=0.5, scale=plate_scale)
    ax.set_xlabel('X (deg)')
    ax.set_ylabel('Y (deg)')
    ax.grid(ls='--')

    # Move the source
    source.directions.altitude += 1. # degree

plt.tight_layout()
plt.show()

Visualize your geometry

For optical systems with complex geometry, it is often useful to perform a visual check of the geometry. To do so, a VTK visualizer is provided:

renderer = iactsim.visualization.VTKOpticalSystem(segmented_telescope.optical_system)
renderer.start_render()