Simulation of increased NSB at the waveforms level

[moralejo]

Related to #1325

We intend to use the reconstructed pixel charges in interleaved pedestal events to quickly determine (via a parametrization) the NSB rate (p.e./ns) that has to be added to "dark MC" waveforms in order to match the observations. Note that lstchain.image.modifier.calculate_required_additional_nsb already does that by actually adding the noise to the MC and comparing the variances of the charges. That is rather slow; what we aim for with #1325 is a faster version of that, which can be used in the analysis pipeline to determine the NSB level and also the cleaning settings for a given run. The idea is simply to parametrize the needed additional NSB rate as a function of some measure of the typical distribution of pixel charges for the interleaved events.

We first tried to use the median reconstructed pixel charge (for LocalPeak pulse integrator) in interleaved pedestals as a proxy for the NSB level. We compute that for real data runs (details in the code proposed in #1325) and also for MC produced with different levels of NSB rate added to the waveforms. Below is what we get for a sample of 100 subruns from different runs. Each plot shows a number of subruns, compared to the two MCs that bracket them in terms of median pedestal charge (and hence in the "added NSB rate" which we obtain from a parabolic parametrization of the behaviour for the MC):

Comments:

• Ideally, the real data curves should be contained between the two bracketing MC curves
• In the plots above we used each single subrun that is plotted to find its matching MC. In real analysis we will find matching values using a subsample of the run's subruns (10, uniformly distributed), because we want to keep the same analysis settings in the whole run (this of course means that if the run is unstable the matching will be worse than shown above).
• At low NSB, there is a disagreement at large charge values. This comes from the fact that the simulated single-p.e-response spectrum ("gain fluctuations") is not perfect, and in particular the MC has a larger tail than the one in real data (from afterpulses).
• The NSB added on the waveforms in lstchain uses (I believe) a single p.e. response which is only for the prompt signal, and has no afterpulsing component (whereas in sim_telarray afterpulses are simulated). This would explain that as the NSB increases the agreement between data and simulation becomes better (since in real data the absolute afterpulsing rate would increase, and we are not reproducing that in our NSB-adding approach)
• For the highest NSB levels the tail in the data actually becomes longer than in the simulation. It seems that the median charge is no longer a good parameter to find the matching MC. For any given value of the median (the same for data and MC), the data charge distribution is wider than the MC one. Note that what matters for the later analysis is to have a good match on the right-side tail of the distribution, i.e. in the frequency of large pulses that have a chance to survive the image cleaning, so having the same median is not enough for a realistic simulation of the noise.
• We tried to use instead the pixel-wise pedestal charge variance (the median of that over the camera), which, being more sensitive to tails, improves the data-MC match, especially for high NSB. The variance is also what lstchain.image.modifier.calculate_required_additional_nsb uses.
• An even better match is achieved by using the 0.95-quantile pedestal pixel charge (the median of it over the camera). By construction, this achieve a better data-MC match on the right-side tail. The result is this one:

With the new approach (95% quantile as parameter to determine the NSB level) the agreement between data and MC on the right-side tail is much better. Only at the highest simulated NSB levels, and for pulse frequencies below ~1% of the peak frequency we see some of the the real data curves slightly exceeding the limits of the MC bracketing curves.

[moralejo]

Note for future improvement:

Use a better single-pe-response in sim_telarray, and also for adding NSB at the waveforms level (including the afterpulsing). See:

cta-observatory/lst-sim-config#40

I do not think however the effect of this improvement can be large, because the discrepancy from the afterpulsing starts at pulses with frequencies of less than 0.1% of the peak frequency. With our standard cleaning approach such pulses will very rarely survive cleaning.