Fix the normalization for muon analysis to account for different PixelShapes. by burmist-git · Pull Request #2845 · cta-observatory/ctapipe

burmist-git · 2025-09-25T12:19:08Z

issue : #2844 (comment)
Parent issue : #2815

The goal is to improve the existing intensity fitter and prepare it for adding another method to measure the impact point and optical efficiency.

This pull request contains two main tasks:

Include phi0 as a parameter in chord_length(radius, rho, phi) in image.muon.intensity_fitter:
Fix the normalization to account for different PixelShapes. Add a function to calculate the trivial solution for the number of photons (from a muon) incident on the telescope mirror, and include it in the tests.

…otons (from muon) incident on the telescope mirror.

…factor in this way.

…e to the mirror radius.

mexanick

@burmist-git thank you for this PR. I've posted a couple of quick comments above. In addition, please add the changelog statement and fix the issues identified by sonar.

mexanick · 2025-09-29T08:28:53Z

+@vectorize(
+    [double(double, double, double, double)], cache=not CTAPIPE_DISABLE_NUMBA_CACHE
+)
+def chord_length(radius, rho, phi0, phi):


What's the reason to break the order of mandatory arguments? Can't we get away with

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

?

Not just the order, also the meaning.

Rho is now absolute instead of a fraction of the radius. Why? That seems like completely needless breaking change, just moving one division or multiplication around.

Yes, the rho and phi0 were changed. And indeed, we can reduce some of the changes.
However, we are modifying this function for a good reason: to make it easier to read and to reuse the same function in the new method.

Old version arguments: (scalar) mirror radius, (scalar) impact parameter relative to the mirror, omitting phi0 and using only one vector phi.

New version arguments: (scalar) mirror radius, (scalar) impact parameter, (scalar) phi0, and one vector phi.

The latest implementation is more natural for three reasons:

It includes all the required parameters in a single function.

If Rmirror is zero, we need to check this before performing the division by Rmirror, since it enters the MINUIT minimizer. It is more natural to handle this inside the function itself, making it more self-contained. For example:

if Rmirror (e.g. for a hole) is 0, the function simply provides the correct answer (zero).

The argument ordering is clearer: (scalar), (scalar), (scalar), (vector).

This function is used only in this module, so I don’t really understand your concerns.

@maxnoe @mexanick please give feed back on this point. However, the current code already incorporates your suggestions, but my comments are still valid. I still believe we need to make this function more straightforward.

mexanick · 2025-09-29T08:30:22Z

+    phi = phi - phi0
+
+    phi_modulo = (phi + np.pi) % (2 * np.pi) - np.pi
+    if phi < 0:


use math.abs() or np.abs()

If you are commenting about this line:
phi_modulo *= -1
then it is correct. Since modulo returns sign less value and we need to return the minus sign back.

looking closer at your code, I don't understand why do you even need this. phi_modulo is used in

np.sin(phi_modulo) ** 2

np.cos(phi_modulo)

np.abs(phi_modulo)

The result of the operations above doesn't depend on the sign of phi_modulo.

Ok i see. Let me run the tests.

Isn't reflecting at the origin just wrong here? Shouldn't you add 2π instead?

Consider e.g. a small negative value (you'd change -π/4 to +π/4 for example).

It is not needed there at all, the line above returns phi_modulo correctly normalized within [-π, π].

This is unnecessary
phi_modulo *= -1

Here is the test without it

mexanick · 2025-09-29T16:32:48Z

    from ctapipe.image.muon.intensity_fitter import chord_length

-    length = chord_length(radius, rho, phi.to_value(u.rad))
+    length = chord_length(radius, rho, phi0.to_value(u.rad), phi.to_value(u.rad))


Looks like the order of parameters is wrong now (phi <-> phi0)

mexanick · 2025-09-29T16:35:30Z



-def intersect_circle(mirror_radius, r, angle, hole_radius=0):
+def intersect_circle(mirror_radius, r, phi0, angle, hole_radius=0):


phi0 is not used in the function, so there's no need (at least at the moment) to add it.
Also, when adding, considering adding it as a last-order keyword argument to preserve the legacy calls.

This has been corrected. It was the reason this wrapper was written.

burmist-git · 2025-09-29T16:48:13Z

The code is ready for review. Sonar reports incorrect code coverage for functions using vectorization.

burmist-git · 2025-09-30T08:42:44Z

The wrapper did not change the time.

On my laptop, the log-likelihood method takes on average 7 minutes and 40 seconds to run for 9,809 muon fits.
With the new method, the runtime is measured to be about 3 minutes at most.

mexanick

LGTM

burmist-git · 2025-10-06T07:20:07Z

@maxnoe Any feed back on this MR ?

maxnoe · 2025-10-06T08:53:57Z

+
+    phi_modulo = (phi + np.pi) % (2 * np.pi) - np.pi
+
+    rho_r = np.abs(rho) / radius


What is the motivation for this redefinition? It seems like an arbitrary choice, so why change it now?

Please revert.

Also, calling abs internally on a quantity that you expect to be always positive might result in very hard debug issues as it just hides wrong input.

@maxnoe concerning this : #2845 (comment)
here is the reason why :

So, in other words, we can minimize the symmetric function with respect to rho. And of course, the same logic is applied to phi.

https://root.cern.ch/download/minuit.pdf

@maxnoe concerning this : #2845 (comment)

Here I put this as explanation :
#2845 (comment)

Then I'd expect this transformation as part of the likelihood definition, not in the low-level function.

I think here it is misleading.

What should we do if this value is negative? Raise an error?

If the value is negative, it still returns a result. This is non-physical and incorrect, so we cannot output it without handling it properly.

Yes, raising a ValueError(f"r must be > 0, got {r}") seems appropriate here.

maxnoe · 2025-10-06T08:57:18Z

@@ -317,7 +344,12 @@ def image_prediction_no_units(
    # diameter. In any case, since in the end we do a data-MC comparison of the muon
    # ring analysis outputs, it is not critical that this value is exact.


The comment here provided by @moralejo is no longer correct, since you introduce the correction for the pixel shape.

Please update it.

Please also have a look here @moralejo.

@moralejo Is this is ok to replace. In the last two sentences… :

# Actually, for the large rings (relative to pixel # size) we are concerned with, a good enough approximation is the ratio between a # circle's area and that of the square whose side is equal to the circle's # diameter or flat-to-flat distance for hexagonal pixel.

maxnoe

Please better explain why the change of API of these functions is needed or an improvement (both adding phi0 and changing the definition of rho).

If we decide that the absolute is better, then the parameter should also be renamed, since rho is used in the current definition by the linked references.

I.e. rho should only be used for the relative quantity.

In general, please separate bug fixes like the modulo check from improvements like the pixel type scaling from API changes.

burmist-git · 2025-10-06T13:39:36Z

Please better explain why the change of API of these functions is needed or an improvement (both adding phi0 and changing the definition of rho).

Adding phi0 and changing rho (to be renamed to impact_parameter) need to be done for at least three reasons:

The function becomes clearer — there is no confusion between the mirror radius and the impact parameter from the mirror’s center, as both are now expressed in the same units.
phi0 is explicitly defined, consistent with one of the reference publications we are using for the muon analysis code.
These changes are preparatory steps for integrating the new intensity fitter (the one you’ve already seen). This time, however, we’ll proceed step by step in separate MRs.

If we decide that the absolute is better, then the parameter should also be renamed, since rho is used in the current definition by the linked references.

Sure - we can rename it to impact_parameter.

I.e. rho should only be used for the relative quantity.

Ok.

In general, please separate bug fixes like the modulo check from improvements like the pixel type scaling from API changes.

Ok.

maxnoe · 2025-10-09T09:12:24Z

    radius: float or ndarray
        radius of circle
-    rho: float or ndarray
+    impact_point: float or ndarray


point implies 2d, but we pass radial coordinates r, phi0. So I'd call it impact_distance (as also the docsting already uses)

Dear @maxnoe, this is a friendly reminder about this MR — let’s merge it :-)

ctao-sonarqube · 2025-10-09T09:50:37Z

Quality Gate failed

Failed conditions
69.8% Coverage on New Code (required ≥ 80%)

See analysis details on SonarQube

maxnoe · 2025-10-20T09:09:38Z

    pixels_on_circle = int(circumference / pixel_diameter)

-    ang = phi + linspace_two_pi(pixels_on_circle * oversampling)
+    ang = linspace_two_pi(pixels_on_circle * oversampling) - phi


Please explain this change. Why do you change the order of this array?

https://github.com/cta-observatory/ctapipe/pull/2845/files/397903ad39453ed0227fb8928d124e4b95c26e92#r2444300605

These two changes are linked. The reason is that phi0 should be subtracted from phi. In this way, the fitted value of phi0 gives us the azimuth of the impact point, whereas before it was giving us the complementary angle of the impact point. This was leading to inconsistencies between the reconstructed and true xcore and ycore values. In addition to the coordinate frame transformation.

But phi0 is subtracted from phi. Now you again subtract phi from the array of all angles.

And: if you are changing the coordinate system of the muon fit, don't you think this is a major change that needs to be mentioned both in the motivation for the pull request and in the changelog?

Please, explain what you are doing so that people can actually review and understand.

Please note that phi0 is set to 0. As we agree do not change the current fitter.

Even more misleading.

Please... do not mix bug fixes with such changes.

I would propose to close this merge request.

Please re-open the many changes you do here in a way that they are separate, their impact can be properly explained, assessed and understood.

Well reuse this one for Fix the normalization for muon analysis to account for different PixelShapes.
I will rename the MR.

maxnoe · 2025-10-20T09:10:08Z

    ang = np.arctan2(dy, dx)
    # Add muon rotation angle
-    ang += phi_rad
+    ang -= phi_rad


Why do you change here + to -?

The answer is here : https://github.com/cta-observatory/ctapipe/pull/2845/files/397903ad39453ed0227fb8928d124e4b95c26e92#r2444298627

burmist-git added 4 commits September 25, 2025 10:57

Add a function to calculate the trivial solution for the number of ph…

dad7a92

…otons (from muon) incident on the telescope mirror.

Add the PixelShape. For image_prediction, we fixed the normalization …

02bbc77

…factor in this way.

add Add the PixelShape

9503e9b

Add test for normalisation factor

cb3a873

burmist-git marked this pull request as draft September 25, 2025 12:19

burmist-git added 2 commits September 26, 2025 10:25

Enhance the chord_length function by applying modulo two pi

81a1da3

Change the parameter rho to absolute values instead of values relativ…

2d92bc4

…e to the mirror radius.

burmist-git requested review from maxnoe and mexanick September 26, 2025 15:19

burmist-git marked this pull request as ready for review September 26, 2025 15:19

burmist-git added 2 commits September 26, 2025 17:21

Add phi0 and change phi signe.

bb86b80

Add phi0.

4290711

This was referenced Sep 26, 2025

Update the chord_length function to use the phi0, add add the PixelShape for image_prediction #2844

Open

Enhancement of the muon analysis pipeline. #2760

Open

mexanick requested changes Sep 29, 2025

View reviewed changes

burmist-git added 8 commits September 29, 2025 12:47

Documentation update. Add desctiption of phi0.

cd51979

Remove uncovered code.

1ad9d3b

Removing unnecessary sign verification and correction.

00a01d1

Add a wrapper to preserve the functional API.

5d30bbc

from camel case to snake case

3d5d805

from camel case to snake case

fa0e5a0

add the bugfix changelog statement

85659cc

Add the bugfix changelog statement - for chord calculation

b7dd762

mexanick reviewed Sep 29, 2025

View reviewed changes

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

mexanick reviewed Sep 29, 2025

View reviewed changes

burmist-git added 2 commits September 29, 2025 18:37

Get back to initail intersect_circle API

c2ca5fe

swap phi and phi0

037b8d4

mexanick self-requested a review September 30, 2025 08:49

mexanick previously approved these changes Sep 30, 2025

View reviewed changes

burmist-git requested a review from kosack October 6, 2025 07:26

maxnoe reviewed Oct 6, 2025

View reviewed changes

maxnoe requested changes Oct 6, 2025

View reviewed changes

rho renamed to impact parameter

b147ccb

burmist-git dismissed mexanick’s stale review via b147ccb October 9, 2025 09:09

maxnoe reviewed Oct 9, 2025

View reviewed changes

burmist-git added 2 commits October 9, 2025 11:19

Adopt the comment with respect to the current changes.

f85c09e

impact parameter renamed to impact_distance

397903a

burmist-git requested a review from maxnoe October 20, 2025 08:55

maxnoe reviewed Oct 20, 2025

View reviewed changes

burmist-git changed the title ~~Muon fit chord length with phi0~~ Fix the normalization for muon analysis to account for different PixelShapes. Oct 20, 2025

burmist-git closed this Oct 20, 2025

burmist-git deleted the muon_fit_chord_length_with_phi0 branch October 20, 2025 13:16

This was referenced Oct 20, 2025

Fix muon normalization pixel shapes #2865

Merged

Phi0 should be subtracted from phi not added. #2869

Merged

add phi0 to chord_length #2872

Merged



		def intersect_circle(mirror_radius, r, angle, hole_radius=0):
		def intersect_circle(mirror_radius, r, phi0, angle, hole_radius=0):


		phi_modulo = (phi + np.pi) % (2 * np.pi) - np.pi

		rho_r = np.abs(rho) / radius

		@@ -317,7 +344,12 @@ def image_prediction_no_units(
		# diameter. In any case, since in the end we do a data-MC comparison of the muon
		# ring analysis outputs, it is not critical that this value is exact.

Conversation

burmist-git commented Sep 25, 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!

burmist-git Sep 29, 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!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

burmist-git Sep 29, 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!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

burmist-git Sep 29, 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!

burmist-git commented Sep 29, 2025

Uh oh!

burmist-git commented Sep 30, 2025

Uh oh!

mexanick left a comment

Choose a reason for hiding this comment

Uh oh!

burmist-git commented Oct 6, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

burmist-git Oct 6, 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!

maxnoe Oct 20, 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!

maxnoe Oct 6, 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!

burmist-git Oct 6, 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!

burmist-git commented Sep 25, 2025 •

edited

Loading

burmist-git Sep 29, 2025 •

edited

Loading

burmist-git Sep 29, 2025 •

edited

Loading

burmist-git Sep 29, 2025 •

edited

Loading

burmist-git Oct 6, 2025 •

edited

Loading

maxnoe Oct 20, 2025 •

edited

Loading

maxnoe Oct 6, 2025 •

edited

Loading

burmist-git Oct 6, 2025 •

edited

Loading

ctao-sonarqube Bot commented Oct 9, 2025 •

edited

Loading

burmist-git Oct 20, 2025 •

edited

Loading