How to add peak arrival times?

[LukasNickel]

I was wondering how one could go ahead adding actual arrival times to the DL1 Container.
From my understanding ctapipe only offers peakpos as integer for the position of the extracted peak(?).

As this work is done in image/charge_extractors.py I should probably add a method to at least one of the classes there.

Now the Question becomes: Where would i put such a function? Put an abstract method in ChargeExtractor, implementation in the most specialized classes (-> Neighbourpeak, AverageWF) or do I follow the same scheme as peakpos (defining tons of intermediate functions for the classes between) or could I just add it to a more abstract class as it would be the same for each Integrator (although it probably wouldn't, I dont know)?

Thanks for all kind of input :)

[watsonjj]

I am refactoring the ChargeExtractor classes at the moment (PR coming soon). So this is a great time to have this discussion.

How is peak arrival time often defined? Is it calculated as part of the charge extraction routine (a byproduct like peakpos), or is it an independently calculated value?

[kosack]

The fact that peakpos is an integer is just an artefact of the simple charge integration method used so far (which is just an integer sum of waveform samples). If you want a float value (also for the integral), then we need a fancier charge integration method with interpolation, proper integration and peak detection. Alternately, we could apply some waveform filtering to up-sample and filter the waveforms, and your "integer" peakpos becomes a fraction (which is easy, but maybe not the fastest method).

[LukasNickel]

So we dont really two separate values but replace peakpos with an actual float value?

From my limited understanding (and from what @maxnoe told me a while ago) you could define an arrival time as weighted mean over the all samples where the weights would be the values for each sample.
This would not require interpolation of the waveforms as far as I understand it.

[LukasNickel]

So instead of doing sum_data.argmax(2).astype(np.int) to return the index of the sample with the highest value you could instead take the weighted mean over the sample indices and receive a float.
(sum_data.shape = n_channels, n_pixels, n_samples)

[thomasgas]

more here maybe (#907)

[maxnoe]

AFAIK, MAGIC uses a weighted average, so np.average(time, weights=charge)

FACT-Tools fits the rising edge of the pulse using a third order polynomial, and defines arrival time as its point of inflection. (Don't ask me how we ended here. I guess it's rather robust against noise).

[kpfrang]

To my knowledge the weighted average is also used in EventDisplay.

[watsonjj]

Great. So my intention is to replace the ChargeExtractor class with a WaveformReducer class, which returns an extracted charge and a peak time.

The inheritance approach of the ChargeExtractor will be abandoned for a composition design, where a WaveformReducer is built up from some base functions existing in the module.

The peak time returned will be a float calculated from a method such as weighted average. Some WaveformReducers can use this technique for extracting the time, other WaveformReducers may implement a different technique for extracting the peak time (and charge).

[watsonjj]

This would also remove the waveform cleaners module, as the waveform cleaning will be part of the charge/time extraction in the WaveformReducer

[watsonjj]

This was previously discussed in #907

[watsonjj]

When calculating the pulse time using the weighted average, should it be done for the samples inside the integration window defined by the peak finding, or should it be done using all samples in the waveform?

[kosack]

Probably depends on the peak-finding algorithm: for some methods a window that is larger than the integration window may be needed (especially if the integration only covers a sub-set of the full pulse).

[watsonjj]

And would one expect the pulse time from the NeighborPeakIntegrator to be calculated from the average waveform of the neighbours to the pixel? Or would it be calculated from the waveform of the pixel?

[watsonjj]

Addressed in #1033