New 3D shower model

[StFroese]

This is just a draft for now to review some ideas.

This approach is based on this paper.

The shower is basically a 3D Gaussian distribution with a width $\sigma_t $ and length $\sigma_l$.
The orientation is given by the azimuthal and polar angle. Those are used to rotate the covariance matrix.
The barycentre of the shower is calculated using the height of the first interaction in the atmosphere and the length of the shower. The centre is shifted according to the intersection of the shower axis on ground.

The next step is to write an ImagePredictor class which will take a subarray, goes through each telescope and integrates the 3D distribution as it can be seen for each pixel in the camera. This could be done with scipy's tplquad but this is computationally very expensive. A better solution is to use the analytical solution provided by the paper (see Appendix).

[codecov]

Codecov Report

❗ No coverage uploaded for pull request base (main@f50ad38). Click here to learn what that means.
Patch has no changes to coverable lines.

❗ Current head 4866052 differs from pull request most recent head 4fcfdd9. Consider uploading reports for the commit 4fcfdd9 to get more accurate results

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #1996   +/-   ##
=======================================
  Coverage        ?   92.77%           
=======================================
  Files           ?      220           
  Lines           ?    18167           
  Branches        ?        0           
=======================================
  Hits            ?    16855           
  Misses          ?     1312           
  Partials        ?        0           

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

☔ View full report in Codecov by Sentry.
📢 Do you have feedback about the report comment? Let us know in this issue.

[vuillaut]

Appendix to be found where?

[StFroese]

@vuillaut Sry, I meant the appendix of the paper, p.29

[vuillaut]

@vuillaut Sry, I meant the appendix of the paper, p.29

Ha yes thanks !
Please include reference to the paper in the docstring of thhe function as well.

[vuillaut]

If reminds me of a forgotten package 😄
https://github.com/vuillaut/pschitt

[kosack]

If reminds me of a forgotten package 😄
https://github.com/vuillaut/pschitt

I remember remarking once when you presented Pschitt that it could be used for reconstruction! All you need to do is add a likelihood function and minimize it

[StFroese]

Predictor now generates telescope images for a given shower model and some input parameters for the telescopes/subarray. It's a little bit slow because of some loops (I'm sure this can be done faster :D).
Looks ok for me so far:

NominalFrame of 4 LSTs:

[kosack]

Nice work! This is progressing fast. I will try to do a more thorough review later, but for now here are just a few minor comments...

[StFroese]

I changed the parameter names to more recognizable ones as recommended by Karl. The _photons function of the Predictor is now vectorized and takes an array of pixels as input :)

[StFroese]

We also need to add an energy estimator to the model since we only get the number of photons. The paper recommends obtaining a relation from simulations based on number of photons, the telescope multiplicity, the telescope pointing and the direction of the shower axis. If I understand them correctly they interpolate their data to get a relation $E(N_{photons}, \vec{\theta})$ with $\vec{\theta}$ being the other parameters mentioned.

I think another approach would be to use the EnergyRegressor with the number of photons given as a feature, but I'm not sure. What do you think @kosack @maxnoe ?

[kosack]

another approach would be to use the EnergyRegressor with the number of photons given as a feature, but I'm not sure. What do you think @kosack @maxnoe ?

It think using an EnergyRegressor is a better idea than using a lookup table as they do. It makes it the same sort of training as for other reconstructions, and we avoid the usual problems with binning and edge cases.

[maxnoe]

If we'd really want to use a lookup table, my approach would be to implement one as a custom scikit learn regressor and then we can use the machine learning module structure the same way.

[kosack]

If we'd really want to use a lookup table, my approach would be to implement one as a custom scikit learn regressor and then we can use the machine learning module structure the same way.

This reminds me that we might want to also think about merging the ImPACT models and training with the ML module eventually - they are basically the same thing: a model that predicts an image, and should be trainable using the same mechanism. I'll open an issue to remind us

[StFroese]

The reconstructor seems to work now. The images show the reconstruction of a 4 LST event. Upper image is the reconstruction and the lower image is the true one.

[StFroese]

There are some open questions/tasks I think:

• Need a good seed for the total number of photons in the shower. I guess something like the average photon number from the HillasReconstructor times some scaling factor. This should be something like a ratio between telescope mirror area and area of the Cherenkov light cone?
• Faster code. This currently takes 6min for one event and 4 LSTs.
• Container for fit parameters. As suggested by @kosack I renamed the reconstructor Model3DReconstructor but did not include Geometry in the name since it should be able to also reconstruct the energy and also fitted parameters would be lost if we only write a ReconstructedGeometryContainer to the event structure. A structure that would make sense to me would be: Let Model3DReconstructor fit the event and save all reconstructed parameters to a new container Reconstructed3DContainer. Model3DGeometryReconstructor and Model3DEnergyReconstructor then can use the 3d container to fill the geometry and energy container

(additional: test_shower_processor_geometry failes, I will fix it with the next commit)

[maxnoe]

in the name since it should be able to also reconstruct the energy

I discussed this yesterday with Karl and we came to the conclusion that, no, this does not predict energy. Rather you'd use the number of photons this reconstructor produces as input feature for the machine learning EnergyRegressor.

[StFroese]

Rather you'd use the number of photons this reconstructor produces as input feature for the machine learning EnergyRegressor.

ok great, but this leaves the question where to store this number. The ReconstructedGeometryContainer only provides average_intensity.

[StFroese]

Test now takes ~8s instead of 6min

[maxnoe]

There are conflicts that need to be resolved, please rebase / merge the current main