No check for fit convergence during muon rings fitting

[jstvdk]

Please describe the use case that requires this feature.

During testing muon ring calibration I faced very large systematic uncertainties. After investigation, I found out that it is caused by the fitting procedure, which has not checked for convergence and boundary conditions. Therefore even if the fit is not converged or converged badly, its result is still used in further analysis, and we are getting optical throughput with high uncertainty.

I attach here just for example a few images of events that are considered as successfully fitted muon rings and optical throughput is extracted from them:

Describe the solution you'd like

Possible solution is the next:

Add several features to the code inside class MuonIntensityFitter from the ctapipe/image/muon/intensity_fitter.py, particularly:

• Check if fit is converged
• Put boundary conditions on the ring width as $mean \pm 2 \sigma$. Because in case of bad fit, this parameter always reach very large values.
• Use in further analysis only muon rings for which fit was converged and converged not into boundary conditions.

I tested this solution on my local ctapipe repository, and it showed very good improvement. For comparison I attach here figures for calculated optical efficiency together with relative uncertainty (in %) for same dataset before/after implementation above mentioned solution.

Before:

After:

[maxnoe]

Hi @jstvdk , thanks for the report.

Your analysis seems sound and I'd add those additional fields (e.g. is_valid for if the fit converged or is at boundary conditions and maybe a goodness-of-fit measure) to the output container and fill them during the fit:

https://github.com/cta-observatory/ctapipe/blob/b3ea362fb1c1d3a783bc82947a051639edfe1703/ctapipe/containers.py#L981

so that badly fit events can be ignored in the analysis.

Just to be sure: you do not improve the fitting itself? This is just about identifying events where the fitting was not successful?

[maxnoe]

An additional note: the example images you linked are clearly not nice muon events.

One is a muon ring overlayed with a shower image and one is just a shower image.

You should be able to get rid of most of these by choosing appropriate quality criteria for the MuonProcessor or later.

E.g. by comparing the light inside the muon ring to the dl1 image size or the concentration_cog parameter, which should be very low for muon rings.

[jstvdk]

Hi @maxnoe,

Just to be sure: you do not improve the fitting itself? This is just about identifying events where the fitting was not successful?

Yes, I was not improving the fit itself. Actually I just add 2 IF statement and 1 limitation to the fit parameters, particularly:

1. One if to check validity of fit
2. Set the fit parameter boundary: minuit.limits["ring_width"] = (0.001169, 0.001518) (this is $mean \pm 2 \sigma$)
3. One more if to check if fit was converged not into this boundary condition $\pm 1e-6$ (maybe there is some more efficient solution)
4. Finally if any of above mentioned conditions is not satisfied, then MuonEfficiencyContainer is filled by nan values

You should be able to get rid of most of these by choosing appropriate quality criteria for the MuonProcessor or later.

E.g. by comparing the light inside the muon ring to the dl1 image size or the concentration_cog parameter, which should be very low for muon rings.

Thanks for the idea, this can improve overall quality. But maybe some bad events still can survive and bring large uncertainty without proper fit control.

[maxnoe]

Finally if any of above mentioned conditions is not satisfied, then MuonEfficiencyContainer is filled by nan values

As said above, I'd add fields for these checks to the container, not replace the values with nan.

[jstvdk]

Hi @maxnoe I was investigating this problem more deeply, to find out some universal way of setting boundary condition on the ring width during the fit. And I found out in several publications (Vacanti et al. 1994; Humensky et al. 2005), that muon ring width should be less then $2 mrad \approx 0.12^{\degree}$. I checked this boundary condition value on my local repo, and it works really well.

[maxnoe]

that muon ring width should be less then 2 mrad ≈ 0.12°

That is a selection cut (in the VERITAS paper), not an upper limit in the fit. Also: this will depend on the PSF, if it is worse, the ring gets broader. Since we want to monitor the PSF, we shouldn't impose an upper limit to it.

[jstvdk]

That is a selection cut (in the VERITAS paper), not an upper limit in the fit. Also: this will depend on the PSF, if it is worse, the ring gets broader. Since we want to monitor the PSF, we shouldn't impose an upper limit to it.

Thanks for the explanation, I was just speculating on how can we choose a universal boundary condition for the muon ring width during the fit, because if this value has no boundary, then the fit quite often wrongly converge.

[maxnoe]

because if this value has no boundary, then the fit quite often wrongly converge.

What do you mean with "wrongly converge"? A nice muon ring gets a too large width? That shouldn't happen.

That the fit does not create a good result for an air shower or a mixed ring and an air shower is a different story and I think we should catch those situations with quality criteria on which rings we include in computatation, not by restricting the fit range.

Restricting the fit range might have the opposite effect we want: produce a result looking like a nice muon ring for events that are indeed not nice rings.

[maxnoe]

I think event after good muon tagging, there are still can left some contaminated or bad muons, and fit can converge badly for for such events which will result in bad optical efficiency value.

I am not talking about criteria before the fit. Also selecting after the fit has been done.

[jstvdk]

What do you mean with "wrongly converge"? A nice muon ring gets a too large width? That shouldn't happen.

I think even after good muon tagging, there are still can leave some contaminated or bad muons, and fit can converge badly for such events which will result in a bad optical efficiency value. They can be very sparse, but their influence on the total uncertainty is huge. A good way to handle such events is to look after their width, because this should have some typical value for given physical parameters of incident muon.

That the fit does not create a good result for an air shower or a mixed ring and an air shower is a different story and I think we should catch those situations with quality criteria on which rings we include in computatation, not by restricting the fit range.

That's a good point of view, but with current quality criteria it is impossible.

Restricting the fit range might have the opposite effect we want: produce a result looking like a nice muon ring for events that are indeed not nice rings.

As far as I experimented with this, I never met a bad muon ring fitted as nice ring due to restriction on the ring width.

[jstvdk]

I am not talking about criteria before the fit. Also selecting after the fit has been done.

Didn't see this reply before my previous comment. That sounds good, I will think about it.

[maxnoe]

@jstvdk There is something weird going on... I wrote that comment as a reply to your previous comment, but now it appears as if it was posted before.

See that I quoted your reply?

[jstvdk]

@jstvdk There is something weird going on... I wrote that comment as a reply to your previous comment, but now it appears as if it was posted before.

See that I quoted your reply?

Yeah, that's my bad, I am sorry, I accidentally sent not finished comment, so I quickly deleted it in order to properly finish it.

I saw your quoted reply, thanks for the explanation.

[maxnoe]

Yeah, that's my bad, I am sorry, I accidentally sent not finished comment, so I quickly deleted it in order to properly finish it.

Ah ok, no problem. Maybe better use edit next time :)

[maxnoe]

That's a good point of view, but with current quality criteria it is impossible.

How come?

[jstvdk]

That's a good point of view, but with current quality criteria it is impossible.

How come?

I played a lot with configurable parameters, and couldn't make a set of criteria that will get rid of all bad muons. Also, I tried your proposed solution with concentration_cog, but there were still left some bad muons with low concentration_cog. Of course, I think it is possible to distinct between good and bad muons, but with some new algorithms or improvements to existing ones.

I will think about it, but earlier the way of controlling ring width seemed to me the easiest one.

[maxnoe]

Just to be sure: will you open a pull request with the changes to the container and the intensity fitter? Adding the fields for the validity, the parameters being inside limits and the likelihood in the minimum?

[jstvdk]

Just to be sure: will you open a pull request with the changes to the container and the intensity fitter? Adding the fields for the validity, the parameters being inside limits and the likelihood in the minimum?

I can do it, though I need more clarifications please as for me still is not clear how we agreed on fit boundary conditions.

Because the previously used mean and $\sigma$ values as boundary limits were calculated only for a certain dataset. And for other datasets, we need other values. So that's why I was thinking of some universal way of setting this.

[mdpunch]

Hi,

I'm wondering what you mean here by "bad muons"? Is this some identification which relies on a neural net (in someone's brain)?

It is indeed possible for muons to undergo scattering during their path, which would give a bad image not useful for calibration. Is it that kind of thing?

But mostly, the muons I see in the simulations of muons are pretty well identified, with a tail of atypical muons. For the calibration, though, surely we just need to choose only the typical muons? Doing this yields very nice muons in mono-telescope proton MCs, see this gif: https://gitlab.cta-observatory.org/punch/nectarcam-muon-candidates/-/raw/main/muons.gif

Is ctapipe intended to use the muons for the calibration? If so, my opinion is that there should be some kind of task force on it, including mainly the people who wrote the big CTA muon calibration paper...

Good Luck, Michael.

Some examples of parameter fits, versus mono-telescope proton images, for the NectarCAM (for muons >20 GeV and hitting the mirror, vs mono-telescope protons)... note the $\log$ y-scale. I could/should include muon width here...

[mdpunch]

Ah, you're probably referring to the fit of the muon ring azimuthally and radially, which requires to have the energy, direction, and impact parameter as input.

Looking at the images in the first post, I would recommend not to even try to fit images with $|I_{\rm ratio} - 1|>0.05$, containment 0.95 and with a limit on the MSE (which may be telescope-dependent).

[maxnoe]

@mdpunch the "task force" is the CalibPipe IKC, of which @jstvdk is a member, and yes, the plan is to use ctapipe.

[mdpunch]

@mdpunch the "task force" is the CalibPipe IKC, of which @jstvdk is a member, and yes, the plan is to use ctapipe.

That's good! Is there a task force page? Are there members from the big muon paper also?

Somehow, though working on muons (for the NectarCAM trigger question) I missed this effort.

[jstvdk]

I'm wondering what you mean here by "bad muons"? Is this some identification which relies on a neural net (in someone's brain)?

It is indeed possible for muons to undergo scattering during their path, which would give a bad image not useful for calibration. Is it that kind of thing?

Hi Michael, thanks for the detailed answer,

I used the term "Bad muons" here as some entity that is not clear muon ring, but can be met in the dataset with only muons, and which is fitted as a muon ring (you can see them on the pictures at the very beginning of this issue, and there are a lot more of them). Due to the successful fit of such events, we get large uncertainty on optical efficiency. We can get rid of such events with a proper muon tagging algorithm, but with the current situation, we don't have such one.

Therefore my idea was to get rid of them by controlling fit and especially ring width because I recognized that strongly contaminated rings and "rings" which are not rings, always characterized by large values of ring width.

[maxnoe]

CalibPipe is part of DPPS, its about using the muons taht are tagged by your work on the trigger (and also the other forms of array calibration).

The team gave a couple of presentations in the computing / dpps sessions of the last meetings.

As far as I know, the authors of the paper are bot directly involved, but I also know that they are quite busy with other things ;)

[mdpunch]

CalibPipe is part of DPPS, its about using the muons taht are tagged by your work on the trigger (and also the other forms of array calibration).

The team gave a couple of presentations in the computing / dpps sessions of the last meetings.

As far as I know, the authors of the paper are bot directly involved, but I also know that they are quite busy with other things ;)

Okay, thanks Max!

I found the talk at https://indico.cta-observatory.org/event/4497/contributions/39216/attachments/23598/34009/dpps_calibration_granada.pdf, though for now it seems to be just a UniGe effort.

I will see if I can drum up interest on the NectarCAM side, unless it is intended to be only UniGe...

[maxnoe]

Looking at the images in the first post, I would recommend not to even try to fit images with $|I_{\rm ratio} - 1|>0.05$, containment 0.95 and with a limit on the MSE (which may be telescope-dependent).

That would be something like this in the config for the MuonProcessor:

MuonProcessor:
  RingQuery:
    quality_criteria:
      - ["np.abs(parameters.intensity_ratio - 1) < 0.05"]
      - ...

[jstvdk]

That would be something like this in the config for the MuonProcessor

Thanks for the detailed explanation, I will try it.

[maxnoe]

I can do it, though I need more clarifications please as for me still is not clear how we agreed on fit boundary conditions.

We can do this in small steps. We agreed we should store the _is_valid flag and maybe the loss function in the minimum. that should be easy and straight forward to add.

I can also open the PR so you what changes that would be and you can work on the rest then.

[jstvdk]

We can do this in small steps. We agreed we should store the _is_valid flag and maybe the loss function in the minimum. that should be easy and straight forward to add.

Yeah, I agree. I just created pull request with starting changes.

[mexanick]

Hi guys, thank you for a nice discussion and useful insights! @maxnoe , considering the boundary conditions and fit convergency: in my opinion, setting reasonable boundary conditions and checking whether fit converged within or at them is always a good idea, given that the boundary conditions are correct. For each telescope, we know (or will know) approximately the size of its PSF and can convolute it with the expected muon ring width quoted by @jstvdk and thus set reasonable boundary conditions, thus optimizing the fit procedure and not altering the analysis of PSF with muons in the future. This will spare us from the successful fit results converging to the muon ring width of 1° or more. @mdpunch I'm very happy to see the NectarCam muon tagging results! I will contact you to organize a chat on your algorithm and simulations, so we can estimate whether we are going to get sufficient amount and sufficient quality of tagged muons for subsequent analysis. From what you've shown above, looks like there's no issues with the quality of muon rings, but the question of statistics is still open to me and we need to check it.

[maxnoe]

This will spare us from the successful fit results converging to the muon ring width of 1° or more.

ok, agreed. I guess that a width of something like 0.2 degrees is a reasonable upper limit then for all telescopes. Of course that could (should?) also be a configurable TelescopeParameter.

[maxnoe]

The fit-convergence status was added to the output container in #2381