2-Pass waveform integration à-la-CTA/MARS

[HealthyPear]

Fixes #602

Prerequisites:

• Merge PR Cleaning for biggest cluster #1131 - continuous integration will fail until this is done.

Description:

• ctapipe-based version of the calibration method used in CTA/MARS and described also here.
• triggered from development of protopipe and comparison against CTA/MARS (in particular DL1, see
  Comparison with CTA-MARS: Calibration protopipe#31 and Comparison with CTA-MARS: Image cleaning & parametrization protopipe#37 for more information)

Caveats:

1. There is at least one cut inside this image extractor, between the 1st and 2nd pass (n_pixels > 3); in ctapipe we shouldn't do this, aka any image extractor should always return charge and pulse times; what I am testing in protopipe and I propose to do also here, is to add a 3rd return argument (at this point to any image extractor) named status, which acts like a "quality label", so later the user can select events on the basis of image extraction quality for more complex image extractors that allow internal selection (like this one),
2. This triggers another problem: in the current protopipe/CTA-MARS what is compared and enters in the benchmark of the charge resolution is the 1st pass, whereas what enters in the final image cleaning is the 2nd pass; this image extractor should then export both passes, but right now in ctapipe only 1 charge and 1 pulse time per image is exported by the call to the calibration method
3. To be as much as close to CTA/MARS as possible, the ctapipe.image.timing_parameters should be modified in order to use a robust fitting algorithm (like scipy.stats.siegelslopes), though I couldn't find an easily importable solution which treated also the fit errors.

[kosack]

I'd suggest removing your name here, and instead adding it to the CODEOWERS file at the top level (we decided against putting names in the comments). You also should also update AUTHORS and .mailmap with your information

[kosack]

However, if there is a reference for the MARS analysis, it would be nice to add it also to docs/references.rst and add a citation link. Note that it should not be in the Returns section (sphinx won't like that), but in either the description or See Also section (see https://numpydoc.readthedocs.io/)

[HealthyPear]

I'll try to see what could be the right citation

[GernotMaier]

If I understood it correctly is this the two-pass integration method described in section 4 of https://ui.adsabs.harvard.edu/abs/2006APh....25..391H/abstract ?

If this is correct, then it maybe shouldn't be called CTA/mars analysis :-)

[HealthyPear]

Hi Gernot, I apologize for the late reply.

I don't know if the version of the algorithm I coded here comes directly from the reference that you provide (in any case, from what you say it seems a necessary reference to be added to the documentation so I will do it).
I just tried with as much fidelity as possible to translate in ctapipe-language the version that is used inside CTA/MARS.
That is the reason for the name of the PR (if this is what you refer to).
I also explicited in the class documentation that "[...] This is in particular the CTA-MARS version.".

For me there is no problem in changing names, perhaps we can ask @moralejo for a confirmation on this.

[moralejo]

Hi, the Eventdisplay reference is indeed to be included: I implemented this in MARS for the CTA analysis (it is not used in MAGIC) to get closer to Eventdisplay. Still, in the spirit of having fully independent analysis codes, I just took the idea, "make a time fit with significant pixels, then do a second round of pulse integration on the rest of the pixels using the times 'predicted' by the fit". The implementation was done in a completely independent manner, and it is bound to differ in the details. I do not think it makes much sense to copy every detail of the Mars version for protopipe - as you @HealthyPear realized, some things (like ROOT's "robust" fit details) would take quite some effort and are probably not essential.

[HealthyPear]

Hi, the Eventdisplay reference is indeed to be included: I implemented this in MARS for the CTA analysis (it is not used in MAGIC) to get closer to Eventdisplay. Still, in the spirit of having fully independent analysis codes, I just took the idea, "make a time fit with significant pixels, then do a second round of pulse integration on the rest of the pixels using the times 'predicted' by the fit". The implementation was done in a completely independent manner, and it is bound to differ in the details.

I will then add the Eventdisplay one; do you have any suggestion for the reference on the CTA/MARS side? Meanwhile, I took the liberty to add the URL of the wiki page as a citation.

[HealthyPear]

I'd suggest removing your name here,

Done

and instead adding it to the CODEOWERS file at the top level (we decided against putting names in the comments).

Done - I hope is the correct syntax

You also should also update AUTHORS and .mailmap with your information

I have already done this sometime ago, do I have to change or update something?

[moralejo]

I will then add the Eventdisplay one; do you have any suggestion for the reference on the CTA/MARS side? Meanwhile, I took the liberty to add the URL of the wiki page as a citation.

Unfortunately we have no public reference for that, since as I said it is not a method we use in the MAGIC analysis. Wiki page is ok for internal reference.

[kosack]

possibly the ICRC proceeding by Maxim Shayduk - I think that's the origin of the method.

[maxnoe]

use relative imports here.

[maxnoe]

@HealthyPear #1267 is merged, if you merge or rebase vs master, this should work now.

[HealthyPear]

ok lots of discussion here! Thanks for you input

[watsonjj]

There are some codacy complaints you should fix

[kosack]

I think if you fix a few Codacy complaints + implement @watsonjj's suggestion to split the __call__() function into two sub-methods, it will be ready.

For call, something like:

def __call__(self, waveforms, telid, selected_gain_channel):
     charge1, pulse_time1 = self.apply_pass_1(waveforms, telid, selected_gain_channel)
     self.apply_pass_2(waveforms, telid, selected_gain_channel, charge1, pulse_time1)

That way, you can just call apply_pass_1 for doing the benchmarks.

[HealthyPear]

I think if you fix a few Codacy complaints + implement @watsonjj's suggestion to split the __call__() function into two sub-methods, it will be ready.

For call, something like:

def __call__(self, waveforms, telid, selected_gain_channel):
     charge1, pulse_time1 = self.apply_pass_1(waveforms, telid, selected_gain_channel)
     self.apply_pass_2(waveforms, telid, selected_gain_channel, charge1, pulse_time1)

That way, you can just call apply_pass_1 for doing the benchmarks.

Yes, it looks like a neat solution....

Also in between the two, there is the preliminary image cleaning from which the actual return values are obtained by checking the result.