PSF model

[ctoennis]

I am adding PSF models to ctapipe.image.psf_model. I made a parent class called PSFModel and a PSF model based on pure coma abberation called ComaModel.

This will be used in the pointing calculation using star tracking and PSF fitting using stars.

[kosack]

Quick question: are these models applicable for both gamma-ray and optical PSF? If not, maybe call it OpticalPSFModel to be more clear where it should be used.

[maxnoe]

I think this is purely optical, not IRF. Also, I think the appropriate place would be in ctapipe.instrument.optics

[ctoennis]

Quick question: are these models applicable for both gamma-ray and optical PSF? If not, maybe call it OpticalPSFModel to be more clear where it should be used.

The model accounts for coma abbarations and the default parameters i have would be for photons in the optical range. With different parameters this should work with gamma photons, but i am not sure. I could call the ComeModel OpticalComaModel and keep the parent class as is. Then later models specifically for gamma photons can be added.

[ctoennis]

I think this is purely optical, not IRF. Also, I think the appropriate place would be in ctapipe.instrument.optics

Makes sense. I will move it to optics.

[kosack]

Side note: we should check with what is used in SimPipe's instrument model, to ensure what is output here is compatible with what the simulations expect.

[maxnoe]

Side note: we should check with what is used in SimPipe's instrument model, to ensure what is output here is compatible with what the simulations expect.

I don't think there is a model like this in SimPipe, but @orelgueta or @GernotMaier can confirm. Rather coma is naturally resulting from the definition of the mirror facets.

[kosack]

I think the appropriate place would be in ctapipe.instrument.optics

The instrument module is where we define the InstrumentDescription model (which currently is collected under SubarrayDescription class here in ctapipe). Would be useful then to integrate it with the rest of the model, i.e. add a place for it in OpticsDescription, and add serialization so that it can be written and read. A collection of these could be serialized to a new table (one row per telescope), as we do with other OpticsDescription parts (SubarrayDescription.to_table(kind="optical_psf"), or else just in OpticsDescription... ) . That could be in another PR if you like, but in that case we should open an issue to make sure this new part of the model is properly integrated

Edit: to be more clear, I would expect something like this to work:

psf : Optional[PSFModel]  = subarray.tel[tel_id].optics.psf

if psf: 
     plot_psf(psf)  # some general function that plots any PSFModel

[ctoennis]

I think the appropriate place would be in ctapipe.instrument.optics

The instrument module is where we define the InstrumentDescription model (which currently is collected under SubarrayDescription class here in ctapipe). Would be useful then to integrate it with the rest of the model, i.e. add a place for it in OpticsDescription, and add serialization so that it can be written and read. A collection of these could be serialized to a new table (one row per telescope), as we do with other OpticsDescription parts (SubarrayDescription.to_table(kind="optical_psf"), or else just in OpticsDescription... ) . That could be in another PR if you like, but in that case we should open an issue to make sure this new part of the model is properly integrated

Edit: to be more clear, I would expect something like this to work:

psf : Optional[PSFModel]  = subarray.tel[tel_id].optics.psf

if psf: 
     plot_psf(psf)  # some general function that plots any PSFModel

if it is okay then i will open an issue for this

[ctoennis]

I think the appropriate place would be in ctapipe.instrument.optics

The instrument module is where we define the InstrumentDescription model (which currently is collected under SubarrayDescription class here in ctapipe). Would be useful then to integrate it with the rest of the model, i.e. add a place for it in OpticsDescription, and add serialization so that it can be written and read. A collection of these could be serialized to a new table (one row per telescope), as we do with other OpticsDescription parts (SubarrayDescription.to_table(kind="optical_psf"), or else just in OpticsDescription... ) . That could be in another PR if you like, but in that case we should open an issue to make sure this new part of the model is properly integrated
Edit: to be more clear, I would expect something like this to work:

psf : Optional[PSFModel]  = subarray.tel[tel_id].optics.psf

if psf: 
     plot_psf(psf)  # some general function that plots any PSFModel

if it is okay then i will open an issue for this

Here it is: #2647

[orelgueta]

Side note: we should check with what is used in SimPipe's instrument model, to ensure what is output here is compatible with what the simulations expect.

I don't think there is a model like this in SimPipe, but @orelgueta or @GernotMaier can confirm. Rather coma is naturally resulting from the definition of the mirror facets.

You are right, it comes out naturally through ray tracing with mirror roughness and misalignment.

[ctoennis]

@kosack @maxnoe I am preparing a PR for the tools to process the variance imaes for the StarTracker that will need the code here. Could you give it another look and see what needs to be changed?

[ctao-dpps-sonarqube]

Analysis Details

0 Issues

• 0 Bugs
• 0 Vulnerabilities
• 0 Code Smells

Coverage and Duplications

• 98.00% Coverage (93.90% Estimated after merge)
• 0.00% Duplicated Code (0.70% Estimated after merge)

Project ID: cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

View in SonarQube

[maxnoe]

Fixtures are not really meant for simple data like this.

You could just enter them in the test itself.

Also: why these highly specific values just for a unit test? Wouldn't nice round values in the right order of magnitude be easier on the eyes?

[maxnoe]

These docstrings are confusing to me. I don't know how to understand the "from where / at where" distinction.

One is the position of the point source and one the position where the PSF is evaluated, right?

Also: We don't really have anything that is using radial coordinates in the camera plane yet.

For the API, I think I'd prefer taking CameraFrame coordinate instances or x / y values. The conversions to polar coordinates can then be done for methods where this is needed for evaluation.

[ctoennis]

i changed the API to use carthesian coordinates and fixed the docustrings to be more clear. I also changed the parameters to round numbers that give reasonable values for some LST-sized camera.

[maxnoe]

Here it still has r,f,r0,f0

[maxnoe]

And if meters are required, then this should probably use u.quantity_input(..) with astropy units

[maxnoe]

This should use the traitlets.validate decorator, one for each traitlet instead of being called manually:

See https://traitlets.readthedocs.io/en/stable/api.html#validating-proposed-changes

[ctoennis]

Sounds good. I will get to it in about an hour.

[maxnoe]

Here it still has r,f,r0,f0

[mexanick]

Looks ok for the moment. I've posted one more suggestion of change, but this can be addressed in the future PR.

[maxnoe]

I find calling this azimuth confusing. This is not the azimuth (of the AltAz coordinate system), it's a polar angle in the camera plane

[maxnoe]

Call this ComaPSFModel?

[maxnoe]

Need an empty line here after the description.

[maxnoe]

The two new classes are currently not documented as we don't go down to optics. Best for now would be to expose the two classes also in ctapipe.instrument. So import them in src/ctapipe/instrument/__init__.py and add them to the __all__ there.

[maxnoe]

I wanted to make a plot of the psf, but this is giving an error:

import astropy.units as u
import numpy as np
from ctapipe.instrument.optics import ComaModel
from ctapipe.instrument import SubarrayDescription
import matplotlib.pyplot as plt

subarray = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst")

lst1 = subarray.select_subarray([1])
psf_model = ComaModel(
    subarray=lst1,
    asymmetry_params=[0.5, 10, 0.15],
    radial_scale_params=[0.015, -0.1, 0.06, 0.03],
    az_scale_params=[0.25, 7.5, 0.02],
)
edges_x = np.linspace(-0.1, 0.1, 250) * u.m
edges_y = np.linspace(-0.1, 0.1, 250) * u.m
centers_x = 0.5 * (edges_x[:-1] + edges_x[1:])
centers_y = 0.5 * (edges_y[:-1] + edges_y[1:])
x, y = np.meshgrid(centers_x, centers_y)

psf = psf_model.pdf(x, y, 0.0 * u.m, 0.0 * u.m)

fig, ax = plt.subplots()

ax.pcolormesh(psf, edges_x.to_value(u.m), edges_y.to_value(u.m), cmap='inferno')

plt.show()

results in:

Traceback (most recent call last):
  File "/home/mnoethe/Uni/CTA/ctapipe/psf.py", line 23, in <module>
    psf = psf_model.pdf(x, y, 0.0 * u.m, 0.0 * u.m)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/mnoethe/Uni/CTA/ctapipe/src/ctapipe/instrument/optics.py", line 411, in pdf
    f, _, r = cartesian_to_spherical((x, y, 0.0))
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (3,) + inhomogeneous part.

It seems the ComaModel.pdf does not support broadcasting, which is I think should be supported.

[maxnoe]

Calling the pdf method should not assign parameters to self

[maxnoe]

needs a [::-1] on the parameters as polyval goes from highest power to lowest power (Horner schema)

Max found a few issues, have to revisit the code.