Skip to content

Commit

Permalink
Add support for using toymodel in telescope frame
Browse files Browse the repository at this point in the history
  • Loading branch information
@maxnoe
maxnoe committed Jun 13, 2023
1 parent bb887b5 commit c520c9a
Show file tree
Hide file tree
Showing 2 changed files with 76 additions and 62 deletions.
6 changes: 3 additions & 3 deletions ctapipe/image/tests/test_toy.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@ def test_intensity(seed, monkeypatch, prod5_lst):
width = 0.05 * u.m
length = 0.15 * u.m
intensity = 200
psi = "30d"
psi = 30 * u.deg

# make sure we set a fixed seed for this test even when testing the
# API without giving the rng
Expand Down Expand Up @@ -69,7 +69,7 @@ def test_skewed(prod5_lst):
width = 0.05 * u.m
length = 0.15 * u.m
intensity = 50
psi = "30d"
psi = 30 * u.deg
skewness = 0.3

model = SkewedGaussian(
Expand All @@ -94,7 +94,7 @@ def test_compare(prod5_lst):
width = 0.05 * u.m
length = 0.15 * u.m
intensity = 50
psi = "30d"
psi = 30 * u.deg

skewed = SkewedGaussian(x=x, y=y, width=width, length=length, psi=psi, skewness=0)
normal = Gaussian(x=x, y=y, width=width, length=length, psi=psi)
Expand Down
132 changes: 73 additions & 59 deletions ctapipe/image/toymodel.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,15 +18,16 @@
>>> print(image.shape)
(400,)
"""
import numpy as np
from ctapipe.utils import linalg
from ctapipe.image.hillas import camera_to_shower_coordinates
import astropy.units as u
from astropy.coordinates import Angle
from scipy.stats import multivariate_normal, skewnorm, norm
from scipy.ndimage import convolve1d
from abc import ABCMeta, abstractmethod

import astropy.units as u
import numpy as np
from numpy.random import default_rng
from scipy.ndimage import convolve1d
from scipy.stats import multivariate_normal, norm, skewnorm

from ctapipe.image.hillas import camera_to_shower_coordinates
from ctapipe.utils import linalg

__all__ = [
"WaveformModel",
Expand All @@ -40,35 +41,33 @@
TOYMODEL_RNG = default_rng(0)


@u.quantity_input(
x=u.m,
y=u.m,
centroid_x=u.m,
centroid_y=u.m,
psi=u.deg,
time_gradient=u.ns / u.m,
time_intercept=u.ns,
)
def _obtain_image_no_units(
x, y, centroid_x, centroid_y, psi, time_gradient, time_intercept
):
longitudinal, _ = camera_to_shower_coordinates(x, y, centroid_x, centroid_y, psi)
return longitudinal * time_gradient + time_intercept


def obtain_time_image(x, y, centroid_x, centroid_y, psi, time_gradient, time_intercept):
"""Create a pulse time image for a toymodel shower. Assumes the time development
occurs only along the longitudinal (major) axis of the shower, and scales
linearly with distance along the axis.
Parameters
----------
x : u.Quantity[length]
x : u.Quantity[angle]
X camera coordinate to evaluate the time at.
Usually the array of pixel X positions
y : u.Quantity[length]
y : u.Quantity[angle]
Y camera coordinate to evaluate the time at.
Usually the array of pixel Y positions
centroid_x : u.Quantity[length]
centroid_x : u.Quantity[angle]
X camera coordinate for the centroid of the shower
centroid_y : u.Quantity[length]
centroid_y : u.Quantity[angle]
Y camera coordinate for the centroid of the shower
psi : convertible to `astropy.coordinates.Angle`
rotation angle about the centroid (0=x-axis)
time_gradient : u.Quantity[time/length]
time_gradient : u.Quantity[time/angle]
Rate at which the time changes with distance along the shower axis
time_intercept : u.Quantity[time]
Pulse time at the shower centroid
Expand All @@ -79,11 +78,27 @@ def obtain_time_image(x, y, centroid_x, centroid_y, psi, time_gradient, time_int
Pulse time in nanoseconds at (x, y)
"""
longitudinal, _ = camera_to_shower_coordinates(x, y, centroid_x, centroid_y, psi)
longitudinal_m = longitudinal.to_value(u.m)
time_gradient_ns_m = time_gradient.to_value(u.ns / u.m)
time_intercept_ns = time_intercept.to_value(u.ns)
return longitudinal_m * time_gradient_ns_m + time_intercept_ns
# camera frame
if x.unit == u.m:
return _obtain_image_no_units(
x=x.to_value(u.m),
y=y.to_value(u.m),
centroid_x=centroid_x.to_value(u.m),
centroid_y=centroid_y.to_value(u.m),
psi=psi.to_value(u.rad),
time_gradient=time_gradient.to_value(u.ns / u.m),
time_intercept=time_intercept.to_value(u.ns),
)

return _obtain_image_no_units(
x=x.to_value(u.deg),
y=y.to_value(u.deg),
centroid_x=centroid_x.to_value(u.dg),
centroid_y=centroid_y.to_value(u.deg),
psi=psi.to_value(u.rad),
time_gradient=time_gradient.to_value(u.ns / u.deg),
time_intercept=time_intercept.to_value(u.ns),
)


class WaveformModel:
Expand Down Expand Up @@ -183,11 +198,9 @@ def from_camera_readout(cls, readout, gain_channel=0):


class ImageModel(metaclass=ABCMeta):
@u.quantity_input(x=u.m, y=u.m)
@abstractmethod
def pdf(self, x, y):
"""Probability density function.
"""
"""Probability density function."""

def generate_image(self, camera, intensity=50, nsb_level_pe=20, rng=None):
"""Generate a randomized DL1 shower image.
Expand Down Expand Up @@ -244,7 +257,6 @@ def expected_signal(self, camera, intensity):


class Gaussian(ImageModel):
@u.quantity_input(x=u.m, y=u.m, length=u.m, width=u.m)
def __init__(self, x, y, length, width, psi):
"""Create 2D Gaussian model for a shower image in a camera.
Expand All @@ -264,32 +276,36 @@ def __init__(self, x, y, length, width, psi):
a `scipy.stats` object
"""
self.unit = x.unit
self.x = x
self.y = y
self.width = width
self.length = length
self.psi = psi

@u.quantity_input(x=u.m, y=u.m)
def pdf(self, x, y):
"""2d probability for photon electrons in the camera plane"""
aligned_covariance = np.array(
[[self.length.to_value(u.m) ** 2, 0], [0, self.width.to_value(u.m) ** 2]]
[
[self.length.to_value(self.unit) ** 2, 0],
[0, self.width.to_value(self.unit) ** 2],
]
)
# rotate by psi angle: C' = R C R+
rotation = linalg.rotation_matrix_2d(self.psi)
rotated_covariance = rotation @ aligned_covariance @ rotation.T
self.dist = multivariate_normal(
mean=u.Quantity([self.x, self.y]).to_value(self.unit),
cov=rotated_covariance,
)

return multivariate_normal(
mean=[self.x.to_value(u.m), self.y.to_value(u.m)], cov=rotated_covariance
).pdf(np.column_stack([x.to_value(u.m), y.to_value(u.m)]))
def pdf(self, x, y):
"""2d probability for photon electrons in the camera plane"""
X = np.column_stack([x.to_value(self.unit), y.to_value(self.unit)])
return self.dist.pdf(X)


class SkewedGaussian(ImageModel):
"""A shower image that has a skewness along the major axis.
"""
"""A shower image that has a skewness along the major axis."""

@u.quantity_input(x=u.m, y=u.m, length=u.m, width=u.m)
def __init__(self, x, y, length, width, psi, skewness):
"""Create 2D skewed Gaussian model for a shower image in a camera.
Skewness is only applied along the main shower axis.
Expand All @@ -313,40 +329,38 @@ def __init__(self, x, y, length, width, psi, skewness):
a `scipy.stats` object
"""
self.unit = x.unit
self.x = x
self.y = y
self.width = width
self.length = length
self.psi = psi
self.skewness = skewness

a, loc, scale = self._moments_to_parameters()
self.long_dist = skewnorm(a=a, loc=loc, scale=scale)
self.trans_dist = norm(loc=0, scale=self.width.to_value(self.unit))
self.rotation = linalg.rotation_matrix_2d(-self.psi)
self.mu = u.Quantity([self.x, self.y]).to_value(u.m)

def _moments_to_parameters(self):
"""Returns loc and scale from mean, std and skewnewss."""
# see https://en.wikipedia.org/wiki/Skew_normal_distribution#Estimation
skew23 = np.abs(self.skewness) ** (2 / 3)
delta = np.sign(self.skewness) * np.sqrt(
(np.pi / 2 * skew23) / (skew23 + (0.5 * (4 - np.pi)) ** (2 / 3))
)
a = delta / np.sqrt(1 - delta ** 2)
scale = self.length.to_value(u.m) / np.sqrt(1 - 2 * delta ** 2 / np.pi)
a = delta / np.sqrt(1 - delta**2)
scale = self.length.to_value(self.unit) / np.sqrt(1 - 2 * delta**2 / np.pi)
loc = -scale * delta * np.sqrt(2 / np.pi)

return a, loc, scale

@u.quantity_input(x=u.m, y=u.m)
def pdf(self, x, y):
"""2d probability for photon electrons in the camera plane."""
mu = u.Quantity([self.x, self.y]).to_value(u.m)

rotation = linalg.rotation_matrix_2d(-Angle(self.psi))
pos = np.column_stack([x.to_value(u.m), y.to_value(u.m)])
long, trans = rotation @ (pos - mu).T

trans_pdf = norm(loc=0, scale=self.width.to_value(u.m)).pdf(trans)

a, loc, scale = self._moments_to_parameters()

return trans_pdf * skewnorm(a=a, loc=loc, scale=scale).pdf(long)
pos = np.column_stack([x.to_value(self.unit), y.to_value(self.unit)])
long, trans = self.rotation @ (pos - self.mu).T
return self.trans_dist.pdf(trans) * self.long_dist.pdf(long)


class RingGaussian(ImageModel):
Expand All @@ -355,17 +369,17 @@ class RingGaussian(ImageModel):
"""

@u.quantity_input(x=u.m, y=u.m, radius=u.m, sigma=u.m)
def __init__(self, x, y, radius, sigma):
self.unit = x.unit
self.x = x
self.y = y
self.sigma = sigma
self.radius = radius
self.dist = norm(
self.radius.to_value(self.unit), self.sigma.to_value(self.unit)
)

@u.quantity_input(x=u.m, y=u.m)
def pdf(self, x, y):
"""2d probability for photon electrons in the camera plane."""

r = np.sqrt((x - self.x) ** 2 + (y - self.y) ** 2)

return norm(self.radius.to_value(u.m), self.sigma.to_value(u.m)).pdf(r)
return self.dist.pdf(r)

0 comments on commit c520c9a

Please sign in to comment.