add impactpoint_intensity_fitter.py by burmist-git · Pull Request #2814 · cta-observatory/ctapipe

burmist-git · 2025-08-20T16:14:18Z

An important error has been found and fixed (in the Chord calculation).
It shows an improvement in the resolution and provides more symmetric distributions.
With the current method, we achieve a 1.7 m resolution on the impact point (from simulation),
whereas ARTEMIS achieved 0.75 m.

We need to implement a new method to reconstruct the impact point.

I will separate the impact point reconstruction from the intensity fitter.
First provides the impact point with known precision.
Second provides the reconstructed and expected intensity and optical throughput.
The intensity fitter will use all available information about the muon to predict
expected intensity.

…st purposes only.

ctao-dpps-sonarqube · 2025-08-21T10:17:14Z

23.08% Duplicated Lines (%) on New Code (is greater than 3.00%)
73.60% Coverage on New Code (is less than 80.00%)

Analysis Details

2 Issues

0 Bugs
0 Vulnerabilities
2 Code Smells

Coverage and Duplications

73.60% Coverage (94.30% Estimated after merge)
23.08% Duplicated Code (0.80% Estimated after merge)

Project ID: cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

View in SonarQube

maxnoe · 2025-08-21T11:06:54Z

If there is an issue in the intensity fitter, why add a second implementation instead of fixing it?

burmist-git · 2025-08-21T11:20:20Z

I fixed the problems and bugs I found so far. It does help, but the resolution is still not enough.

The point is very simple :

I think the key point here is that we perform the fit with fewer parameters, and therefore achieve better resolution.

I will separate the impact point reconstruction from the intensity fitter.
First provides the impact point with known precision.
Second provides the reconstructed and expected intensity and optical throughput.

First, we need to achieve a better fit. If we succeed, we will also learn how to improve the existing one.

maxnoe · 2025-08-21T11:26:28Z

Note that the current intensity fit keeps the ring center fixed (opposed to just using it as an initial guess).

Maybe with the chord bug fixed, we can lift this restriction and the fit will perform better if it is allowed to adjust x/y?

burmist-git · 2025-08-21T11:32:58Z

I still think adding more free parameters will make it worse.

Moreover, it often fails to converge to the global maximum of the log-likelihood.

Another major issue (in current fit) is the camera shadowing, which is not well implemented.

…ng the impact point.

…om_muon_ring sceletons

…ry sensitive to outliers.

…ed_pe

mexanick

Overall, this PR needs significant improvements before it can be merged. At the moment it does not meet our project’s quality standards:
– There is extensive duplicated code instead of reuse, which will make maintenance harder.
– Documentation is largely missing, making it difficult to understand intent and usage.
– Tests are essentially absent (the current placeholder test is insufficient).
– Several PEP8/style conventions are not followed (e.g., single-letter/cryptic variable names, commented-out code, etc).

mexanick · 2025-09-19T08:25:36Z

@@ -0,0 +1,558 @@
+"""


I don't think creation of a new module is justified here. The new fitter class returns the same container as the MuonIntensityFitter in the ctapipe.image.muon.intensity_fitter, therefore it should be also implemented there.

It’s not a new module in my opinion — we need to completely remove the old fitter. The new one has a clearer logic. It is easy to follow and in addition fully vectorized that is why this optimization can be removed.

mexanick · 2025-09-19T08:28:20Z

+]
+
+
+def chord_length(radius, rho, phi0, phi):


This function is an almost exact copy-paste of the chord_length of ctapipe.image.muon.intensity_fitter. Please reuse the code instead of copying it, and do not strip vectorization and other performance optimizations. Besides, in case of the API update, the documentation shall be updated too. A new argument (phi0) was introduced, that breaks the order of arguments and lacks description.

It is far from being an exact copy. This function has two main modifications: first, it uses modulo
2pi, and therefore we remove the limits on the fitting parameter. Second, we explicitly include phi0 (defines the direction with respect to the mirror of the impact point) in the function instead of hiding it somewhere in the fitting procedure this makes it much easier to read.

If this is a general improvement, improve the function. Do not copy and leave a "suboptimal" slightly different version elsewhere.

I'm still not sure where this pull request is going and why a separate fitter is needed at all as opposed to freezing the efficiency and then re-fitting it with the impact frozen.

The current algorithm is inefficient, with many values stuck at the limits or default settings. An initial estimate produced 30% more fitted events, and the proposed method is straightforward and easy to follow. We also decouple the impact point fitting from the optical efficiency.

Then improve the current implementation, please.

I have already implemented a simple and easy-to-debug method, with decoupled impact point width and optical efficiency measurements.

I a parallel structure. If this new method doesn't have disadvantages, there is no point in this duplication. Just improve the existing one please

You can also enable one or the other approach by adding configuration options to MuonIntensityFitter if you think it's valuable to keep both.

Do you mean ‘replace’? Yes, that can be done very quickly. Or do you mean improving the current likelihood? Because the main improvement, in my view, comes from fitting only the phi distribution with only 18 bins in angle.

You can also enable one or the other approach by adding configuration options to MuonIntensityFitter if you think it's valuable to keep both.

In principle, no — we do not need it.

mexanick · 2025-09-19T08:46:04Z

+    return chord_length(radius, rho, 0, phi - phi0)
+
+
+def gauss_pedestal(x, A, mu, sigma, pedestal):


Is this function used to describe the ring width? Please elaborate more the fitting approach.

Previously, the ring width was fitted with a 2D Gaussian, now we are doing the same in 1D. This makes it simpler and easier to read and interpret.

No, the ring width is described using a 1d normal distribution:

ctapipe/src/ctapipe/image/muon/intensity_fitter.py

Line 224 in b2221d5

def gaussian_cdf(x, mu, sig):

@maxnoe It has radial symmetry, but the overall fit is done in 2D. https://github.com/cta-observatory/ctapipe/blob/b2221d5d12bb153c6e0bc0914c3d33ffdaf29bdf/src/ctapipe/image/muon/intensity_fitter.py#L527C1-L527C65

I'm not sure what your point is. The existing likelihood function transforms into polar coordinates as soon as possible and continues from there.

No 2d gaussian fit is performed as you wrote here:

Previously, the ring width was fitted with a 2D Gaussian,

The ring width could be fitted in 1D, whereas it is currently fitted effectively in 2D.

I don't understand your line of argument here. The ring width is a parameter in a likelihood which is described by a 1D probability distribution.

It's not 2d

The likelihood fit itself is in 2D.

And why do you think this is an issue? How else can you take the ring center location into account for the image prediction?

Just like in the 1D fit of the phi distribution.

mexanick · 2025-09-19T08:56:10Z

+    )
+
+    # Predicted total number of Cherenkov photons falling on the mirror
+    pred_total_Cher_phot = (


Why do you re-write this? There's a well-defined function image_prediction_no_units that provides you with this. Also, I believe the difference of normalization here is the reason why you have a different final number. As it is correctly noted in the documentation of image_prediction_no_units, since we are measuring the throughput relative to the simulation, the overall normalization doesn't matter as it cancels out.

This function has issues and errors:

-The expected measured number of signal events does not depend on the pixel FoV, but this function does. The new function fixes this issue.
-There is an analytical formula for the trivial case (rho=0 and no shadowing), which should be included in the tests. The current function gives the wrong result and new function fixes this issue.

mexanick · 2025-09-19T09:00:45Z

+    )
+    fit.errordef = Minuit.LEAST_SQUARES
+    #
+    fit.errors["A"] = 10


How these values are obtained and what units they represent?

This value defines the initial step of the fit. Usually, we set it much smaller and should be changed to an smaller number. Minuit works with scalars only, you cannot define values in units like meters, for example.

mexanick · 2025-09-19T09:03:27Z


-        self.intensity_fitter = MuonIntensityFitter(subarray=subarray, parent=self)
+        # self.intensity_fitter = MuonIntensityFitter(subarray=subarray, parent=self)
+        self.intensity_fitter = MuonImpactpointIntensityFitter(


If you're leaving both options, it should be possible to select one of the fitters.

This is clear form very fest comment: #2814 (comment)

And discussion and motivation of the issue.

mexanick · 2025-09-19T09:04:07Z

@@ -0,0 +1,34 @@
+import astropy.units as u


Please provide a proper testing of the code you've added.

Since I discovered new issues and problems in the code (the lat one today), we need to work through it in more detail. Yes, we definitely need to add more tests.

maxnoe · 2025-09-19T14:49:21Z

+
+    weights = (phi_y > 0).astype(int)
+
+    phi_x = ((np.roll(hist_phi[1], 1) + hist_phi[1]) / 2.0)[1:]


this a very complicated way to say 0.5 * (edges[:-1] + edges[1:])

…using a 1D phi distribution.

ctao-sonarqube · 2025-09-22T15:49:08Z

Quality Gate failed

Failed conditions
13 New issues
64.6% Coverage on New Code (required ≥ 80%)

See analysis details on SonarQube

Catch issues before they fail your Quality Gate with our IDE extension SonarQube for IDE

maxnoe · 2025-09-25T07:38:21Z

+        impact_x=rho * np.cos(Quantity(fit2.values["phi0"], u.rad)),
+        impact_y=rho * np.sin(Quantity(fit2.values["phi0"], u.rad)),
+        chord_fit_ampl=fit2.values["amplitude"],
+        phi_bin_0=phi_y[0],


You can just put the array into the container

Yes, I know. We decided to split this MR into parts, so I just pushed the code to switch the branches. Please ignore this draft MR, until it’s moved out of draft. Also can this MR can be completely deleted i will let you know. I just want to keep the history of changes for current development.

add impactpoint_intensity_fitter.py

47ddd19

This was referenced Aug 20, 2025

New method to reconstruct muon impact point #2815

Open

Enhancement of the muon analysis pipeline. #2760

Open

This comment has been minimized.

Sign in to view

burmist-git added 2 commits August 21, 2025 10:41

Remove the old description. New one to be added.

a3779e5

Removing unnecessary parameters

d050726

This comment has been minimized.

Sign in to view

burmist-git added 4 commits August 21, 2025 11:52

test: add preliminary dummy test

bca8069

Add MuonImpactpointIntensityFitter to the muon processor

b20b256

Change the return value to the default MuonEfficiencyContainer for te…

9bfb139

…st purposes only.

change the name LSTMono -> LST

cebf1c2

burmist-git added 13 commits August 21, 2025 17:16

Introduce a seed muon_ring_phi_distribution_fit function for estimati…

f4fb668

…ng the impact point.

Add staticmethod to count number of calls

f67a0c5

Add temporary printouts for development (to be removed)

bbd630b

add support function for dev. only. To be removed

1493c36

add support function for dev. only. To be removed (processor.py)

041a03a

Add event_id only for dev. TO be removed.

552badb

Add event_id only for dev. To be removed.

ca0bad5

Add smoothing phi dist

64f4b8d

modify chord_length with modulo

afb162f

Add the loss function for the phi-distribution of the Cherenkov light.

c19ef6e

Add working fit with Minuit.

2a771f7

Remove debugging printouts from processor.py

3bef25e

Remove debugging printouts from impactpoint_intensity_fitter.py

4dd8802

burmist-git added 13 commits September 1, 2025 16:23

Cleanup code, add units, and include MuonEfficiencyContainer output

796140b

Cleaning the code and adding calculations for the initial parameters.

2362269

add compute_muon_ring_width and comput_absolute_optical_efficiency_fr…

b742336

…om_muon_ring sceletons

add compute_muon_ring_width function

f8ee916

update save_histogram_to_csv for tests

cd78b24

Test non-fitting algorithm: simple standard deviation calculation. Ve…

d705e59

…ry sensitive to outliers.

add fit_muon_ring_width function

d356921

remove help/debug function

8868a02

Add compute_absolute_optical_efficiency_from_muon_ring and get_measur…

97184cd

…ed_pe

remove prints

71aca3b

save_histogram_to_csv debug function

e81f830

remove coment

13f3077

make the formula easier to read

c4d27df

burmist-git requested a review from mexanick September 18, 2025 13:11

mexanick requested changes Sep 19, 2025

View reviewed changes

maxnoe reviewed Sep 19, 2025

View reviewed changes

burmist-git added 3 commits September 22, 2025 17:11

Preparation to replace the current intensity fitter with a new one - …

28c86ae

…using a 1D phi distribution.

put back mask argument

fe167fa

Update the Class description

677f00b

burmist-git added 2 commits September 24, 2025 18:26

tmp. commit

8d0cf9b

tmp. commit

ad13e56

maxnoe reviewed Sep 25, 2025

View reviewed changes

		return chord_length(radius, rho, 0, phi - phi0)


		def gauss_pedestal(x, A, mu, sigma, pedestal):


		weights = (phi_y > 0).astype(int)

		phi_x = ((np.roll(hist_phi[1], 1) + hist_phi[1]) / 2.0)[1:]

Conversation

burmist-git commented Aug 20, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

This comment has been minimized.

This comment has been minimized.

This comment has been minimized.

ctao-dpps-sonarqube Bot commented Aug 21, 2025

Analysis Details

2 Issues

Coverage and Duplications

Uh oh!

maxnoe commented Aug 21, 2025

Uh oh!

burmist-git commented Aug 21, 2025

Uh oh!

maxnoe commented Aug 21, 2025

Uh oh!

burmist-git commented Aug 21, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mexanick left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

burmist-git Sep 19, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

burmist-git commented Aug 20, 2025 •

edited

Loading

burmist-git commented Aug 21, 2025 •

edited

Loading

burmist-git Sep 19, 2025 •

edited

Loading

burmist-git Sep 25, 2025 •

edited

Loading