From 91cdd0bdb40ed2ae257b664153985f60372d6ced Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Wed, 31 Jul 2024 10:50:05 +0200
Subject: [PATCH 01/21] New version of the code to obtain the NSB tuning
 parameters

---
 lstchain/image/modifier.py | 247 +++++++++++++++++++++++++------------
 1 file changed, 166 insertions(+), 81 deletions(-)

diff --git a/lstchain/image/modifier.py b/lstchain/image/modifier.py
index 9f52390238..6676400b71 100644
--- a/lstchain/image/modifier.py
+++ b/lstchain/image/modifier.py
@@ -6,13 +6,17 @@
     EventSource,
     read_table,
 )
+from ctapipe.containers import EventType
 from numba import njit
 from scipy.interpolate import interp1d
 from traitlets.config import Config
 
-from lstchain.io import standard_config
+from lstchain.io import standard_config, read_configuration_file
 from lstchain.io.config import read_configuration_file
 from lstchain.reco.reconstructorCC import nsb_only_waveforms
+from lstchain.data import NormalizedPulseTemplate
+from lstchain.io.io import get_resource_path
+from scipy.optimize import curve_fit
 
 __all__ = [
     'add_noise_in_pixels',
@@ -407,7 +411,6 @@ def calculate_noise_parameters(simtel_filename, data_dl1_filename,
 
 
 def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config=None):
-    # TODO reduce duplicated code with 'calculate_noise_parameters'
     """
     Calculates the additional NSB needed in the MC waveforms
     to match a real data DL1 file
@@ -428,7 +431,8 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
 
     Returns
     -------
-    extra_nsb: Fraction of the additional NSB in data compared to MC.
+    extra_nsb: additional NSB rate in absolute units, p.e./ns (a.k.a. "GHz"),
+    that has to be added in the MC to match the real data
     data_ped_variance: Pedestal variance from data
     mc_ped_variance: Pedestal variance from MC
 
@@ -440,118 +444,200 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
         config = standard_config
 
     # Real data DL1 tables:
+    camera_images = read_table(data_dl1_filename,
+                               '/dl1/event/telescope/image/LST_LSTCam')
+    # parameters:
+    image_params = read_table(data_dl1_filename,
+                              '/dl1/event/telescope/parameters/LST_LSTCam')
+
+    # All events in a DL1 file should correspond to the same telescope,
+    # just take it from the first event:
+    tel_id = image_params['tel_id'][0]
+
     data_dl1_calibration = read_table(data_dl1_filename,
                                       '/dl1/event/telescope/monitoring/calibration')
-    data_dl1_pedestal = read_table(data_dl1_filename,
-                                   '/dl1/event/telescope/monitoring/pedestal')
     unusable = data_dl1_calibration['unusable_pixels']
+
     # Locate pixels with HG declared unusable either in original calibration or
     # in interleaved events:
     bad_pixels = unusable[0][0]  # original calibration
     for tf in unusable[1:][0]:  # calibrations with interleaved
         bad_pixels = np.logical_or(bad_pixels, tf)
     good_pixels = ~bad_pixels
-
     # First index:  1,2,... = values from interleaved (0 is for original
     # calibration run)
     # Second index: 0 = high gain
     # Third index: pixels
 
-    # HG adc to pe conversion factors from interleaved calibrations:
-    data_HG_dc_to_pe = data_dl1_calibration['dc_to_pe'][:, 0, :]
-    # Pixel-wise pedestal standard deviation (for an unbiased extractor),
-    # in adc counts:
-    data_HG_ped_std = data_dl1_pedestal['charge_std'][1:, 0, :]
-    # indices which connect each pedestal calculation to a given calibration:
-    calibration_id = data_dl1_pedestal['calibration_id'][1:]
-    # convert pedestal st deviations to p.e.
-    dummy = []
-    for i, x in enumerate(data_HG_ped_std[:, ]):
-        dummy.append(x * data_HG_dc_to_pe[calibration_id[i],])
-    dummy = np.array(dummy)
-
-    # Average for all interleaved calibrations (in case there are more than one)
-    data_HG_ped_std_pe = np.mean(dummy, axis=0)  # one value per pixel
-
-    # Identify noisy pixels, likely containing stars - we want to adjust MC to
-    # the average diffuse NSB across the camera
-    data_median_std_ped_pe = np.median(data_HG_ped_std_pe)
-    data_std_std_ped_pe = np.std(data_HG_ped_std_pe)
-    log.info(f'Real data: median across camera of good pixels\' pedestal std '
-             f'{data_median_std_ped_pe:.3f} p.e.')
-    brightness_limit = data_median_std_ped_pe + 3 * data_std_std_ped_pe
-    too_bright_pixels = (data_HG_ped_std_pe > brightness_limit)
-    log.info(f'Number of pixels beyond 3 std dev of median: '
-             f'{too_bright_pixels.sum()}, (above {brightness_limit:.2f} p.e.)')
-
-    # Exclude too bright pixels, besides those with unusable calibration:
-    good_pixels &= ~too_bright_pixels
-    # recalculate the median of the pixels' std dev, with good_pixels:
-    data_median_std_ped_pe = np.median(data_HG_ped_std_pe[good_pixels])
-
+    # Now compute the average pixel charge in real pedestal events, for good
+    # and "not too bright" pixels (the idea is to exclude stars, we want to
+    # estimate the overall "diffuse" NSB level)
+    interleaved_ped = image_params['event_type'] == EventType.SKY_PEDESTAL.value
+    ped_pixq = camera['image'][interleaved_ped]
+    ped_meanpixq = np.mean(ped_pixq, axis=1)
+    ped_stdpixq = np.std(ped_pixq, axis=1)
+    median_ped_meanpixq = np.median(ped_meanpixq[good_pixels])
+    # Exclude the brightest pixels, which may be affected by stars:
+    too_bright = (ped_meanpixq > median_ped_meanpixq + 3 * np.std(ped_meanpixq))
+    good_pixels &= ~too_bright
     log.info(f'Good and not too bright pixels: {good_pixels.sum()}')
+    # Recompyte median:
+    median_ped_meanpixq = np.median(ped_meanpixq[good_pixels])
+    median_ped_stdpixq = np.median(ped_stdpixq[good_pixels])
 
+    # Now we process the Monte Carlo:
     # Event reader for simtel file:
     mc_reader = EventSource(input_url=simtel_filename, config=Config(config))
 
-    # Obtain the configuration with which the pedestal calculations were
-    # performed:
-    ped_config = config['LSTCalibrationCalculator']['PedestalIntegrator']
-    tel_id = ped_config['tel_id']
-    # Obtain the (unbiased) extractor used for pedestal calculations:
-    pedestal_calibrator = CameraCalibrator(
-        image_extractor_type=ped_config['charge_product'],
-        config=Config(config['LSTCalibrationCalculator']),
-        subarray=mc_reader.subarray)
-
-    # Since these extractors are now for use on MC, we have to apply the pulse
-    # integration correction (in data that is currently, as of
-    # lstchain v0.7.5, replaced by an empirical (hard-coded) correction of the
+    subarray = mcreader.subarray
+    sampling_rate = subarray.tel[1].camera.readout.sampling_rate
+    dt = (1.0 / sampling_rate).to(u.ns)
+
+    # Get the single-pe response fluctuations:
+    spe_location = (config['waveform_nsb_tuning']['spe_location']
+                    if 'spe_location' in config['waveform_nsb_tuning']
+                       and config['waveform_nsb_tuning']['spe_location']
+                       is not None
+                    else get_resource_path(
+            "data/SinglePhE_ResponseInPhE_expo2Gaus.dat"))
+    spe = np.loadtxt(spe_location).T
+    spe_integral = np.cumsum(spe[1])
+    charge_spe_cumulative_pdf = interp1d(spe_integral, spe[0], kind='cubic',
+                                         bounds_error=False, fill_value=0.,
+                                         assume_sorted=True)
+    # Pulse template shape for a single p.e.:
+    pulse_templates = {tel_id: NormalizedPulseTemplate.load_from_eventsource(
+            subarray.tel[tel_id].camera.readout, resample=True) for tel_id in
+        config['source_config']['LSTEventSource']['allowed_tels']}
+
+    # Since the pulse integrator will now be used on MC, we have to apply the
+    # pulse integration correction (in data that is currently, as of
+    # lstchain v0.10, replaced by an empirical (hard-coded) correction of the
     # adc to pe conversion factors )
-    pedestal_calibrator.image_extractors[ped_config['charge_product']].apply_integration_correction = True
-
-    # MC pedestals integrated with the unbiased pedestal extractor
-    mc_ped_charges = []
+    config['LocalPeakWindowSum']['apply_integration_correction'] = True
+    r1_dl1_calibrator = CameraCalibrator(image_extractor_type=config[
+        'image_extractor'], config=Config(config), subarray=subarray)
+
+    numevents = 0
+    maxmcevents = 200 # Enough statistics to determine the right NSB level
+
+    # Simulated levels of additional NSB, rate in p.e./ns (a.k.a. GHz):
+    total_added_nsb = np.array([0, 0.125, 0.25, 0.5, 1, 2, 4])
+
+    # Just the product of added_nsb and original_nsb (first two arguments
+    # below) is relevant! Units "GHz" (p.e./ns)
+    nsb_tuner = [None]
+    # We now create the instances of WaveformNsbTunner to add the different
+    # levels of noise to the MC waveforms.
+    # NOTE: the waveform gets updated every time, the noise addition is
+    # cumulative (hence the np.diff):
+
+
+    for nsb_rate in np.diff(total_added_nsb):
+        # Create an array to put the NSB value for tel_id. It is an array for
+        # the future case in which there are more telescopes, but for now we
+        # just create it with a single value for tel_id
+        nsb = np.zeros(tel_id+1)
+        nsb[tel_id] = nsbrate
+        nsb_tuner.append(WaveformNsbTunner(nsb * u.GHz, pulse_templates,
+                                           charge_spe_cumulative_pdf, 10))
+    # last argument means it will precompute 10 * 1855 (pixels) * 2 (gains)
+    # noise waveforms to pick from during the actual introduction of the noise
+
+    # List of lists to keep the integrated pixel charges for different NSB
+    # levels:
+    modified_integrated_charge = [[] for i in range(len(nsb_tuner))]
 
     for event in mc_reader:
         if tel_id not in event.trigger.tels_with_trigger:
             continue
-        # Extract the signals as we do for pedestals (unbiased fixed window
-        # extractor):
-        pedestal_calibrator(event)
-        charges = event.dl1.tel[tel_id].image
+        numevents += 1
+        if numevents > maxmcevents:
+            break
 
-        # True number of pe's from Cherenkov photons (to identify noise-only pixels)
+        # Calibrate the event to get the integrated charges (DL1a):
+        r1_dl1_calibrator(event)
+
+        selected_gains = event.r1.tel[tel_id].selected_gain_channel
+        mask_high = (selected_gains == 0)
         true_image = event.simulation.tel[tel_id].true_image
-        mc_ped_charges.append(charges[true_image == 0])
+        # Use only pixels with no Cherenkov signal, just noise:
+        pedmask = mask_high & (true_image == 0)
 
-    # All pixels behave (for now) in the same way in MC, just put them together
-    mc_ped_charges = np.concatenate(mc_ped_charges)
-    mc_unbiased_std_ped_pe = np.std(mc_ped_charges)
+        # First store the charges with no added NSB:
+        modified_integrated_charge[0].extend(event.dl1.tel[tel_id].image[
+                                                 pedmask])
+
+        # Now add the differentlevels of NSB and recompute charges:
+        for ii, tuner in enumerate(nsb_tuner[1:]):
+            waveform = event.r1.tel[tel_id].waveform
+
+            # NOTE!! The line below modifies the waveform in event.r1
+            tuner.tune_nsb_on_waveform(waveform, tel_id, mask_high, subarray)
+            r1_dl1_calibrator(event)
+            modified_integrated_charge[ii + 1].extend(
+                    event.dl1.tel[1].image[pedmask])
 
-    # Find the additional noise (in data w.r.t. MC) for the unbiased extractor
-    # The idea is that pedestal variance scales with NSB
+        modified_integrated_charge = np.array(modified_integrated_charge)
+        # Fit the total added NSB rate vs. the average pixel charge:
+        params, _ = curve_fit(custom_function,
+                              np.mean(modified_integrated_charge, axis=1),
+                              total_added_nsb, p0=[-0.2, 0.2, 2])
+
+    # Obtain the right rate of NSB to be added to MC so that it matches the
+    # data:
+    extra_nsb = custom_function(median_ped_meanpixq,
+                                params[0], params[1], params[2])
+
+
+    # Now open the MC file again and test that the tuning was successful:
+    mc_reader = EventSource(input_url=simtel_filename, config=Config(config))
+    nsb = np.zeros(tel_id+1)
+    nsb[tel_id] = extra_nsb
+    tuner = WaveformNsbTunner(nsb * u.GHz , pulse_templates,
+                              charge_spe_cumulative_pdf, 10)
+    final_mc_qped = []
+    numevents = 0
+    for event in mcreader:
+        if tel_id not in event.trigger.tels_with_trigger:
+            continue
+        numevents += 1
+        if numevents > maxmcevents:
+            break
+        selected_gains = event.r1.tel[tel_id].selected_gain_channel
+        mask_high = (selected_gains == 0)
+        true_image = event.simulation.tel[tel_id].true_image
+        pedmask = mask_high & (true_image == 0)
+
+        waveform = event.r1.tel[tel_id].waveform
+        tuner.tune_nsb_on_waveform(waveform, tel_id, mask_high, subarray)
+        r1_dl1_calibrator(event)
+        final_mc_qped.extend(event.dl1.tel[tel_id].image[pedmask])
+
+    data_ped_variance = median_ped_stdpixq**2
+    mc_ped_variance = np.std(final_mc_qped)**2
 
-    data_ped_variance = data_median_std_ped_pe ** 2
-    mc_ped_variance = mc_unbiased_std_ped_pe ** 2
-    extra_nsb = ((data_ped_variance - mc_ped_variance)
-                 / mc_ped_variance)
     return extra_nsb, data_ped_variance, mc_ped_variance
 
 
+def custom_function(x, a, b, c):
+    """
+    Custom function to fit the "mean pix charge in pedestal events" vs. the
+    added level of NSB in waveforms
+    """
+    return a + b * x ** c
+
 class WaveformNsbTunner:
     """
     Handles the injection of additional NSB pulses in waveforms.
     """
-    def __init__(self, added_nsb_fraction, original_nsb, pulse_template, charge_spe_cumulative_pdf,
+    def __init__(self, added_nsb, pulse_template, charge_spe_cumulative_pdf,
                  pre_computed_multiplicity=10):
         """
         Parameters
         ----------
-        added_nsb_fraction: `float`
-            Fraction of `original_nsb` to be added during NSB injection. Can be more than 1.
-        original_nsb: `dict`[`int`,`astropy.units.Quantity`]
-            NSB frequency used in the original waveforms per telescope
+        added_nsb: array [`float`,`astropy.units.Quantity`]
+            NSB frequency (GHz) to be added to the waveforms (per telescope)
         pulse_template: `dict`[int,`lstchain.data.NormalizedPulseTemplate`]
             Single photo-electron pulse template per telescope
         charge_spe_cumulative_pdf: `scipy.interpolate.interp1d`
@@ -561,8 +647,7 @@ def __init__(self, added_nsb_fraction, original_nsb, pulse_template, charge_spe_
             waveforms. Later used during event modification. Set to 0 to always compute the correction on the fly.
 
         """
-        self.added_nsb_fraction = added_nsb_fraction
-        self.original_nsb = original_nsb
+        self.added_nsb = added_nsb
         self.pulse_template = pulse_template
         self.charge_spe_cumulative_pdf = charge_spe_cumulative_pdf
         self.multiplicity = pre_computed_multiplicity
@@ -596,11 +681,11 @@ def initialise_waveforms(self, waveform, dt, tel_id):
         """
         log.info(f"Pre-generating nsb waveforms for nsb tuning and telescope id {tel_id}.")
         n_pixels, n_samples = waveform.shape
-        baseline_correction = -(self.added_nsb_fraction * self.original_nsb[tel_id] * dt).to_value("")
+        baseline_correction = -(self.added_nsb[tel_id] * dt).to_value("")
         nsb_waveforms = np.full((self.multiplicity * n_pixels, 2, n_samples), baseline_correction)
         duration = (self.extra_samples + n_samples) * dt
         t = np.arange(-self.extra_samples, n_samples) * dt.to_value(u.ns)
-        mean_added_nsb = (self.added_nsb_fraction * self.original_nsb[tel_id] * duration).to_value("")
+        mean_added_nsb = (self.added_nsb[tel_id] * duration).to_value("")
         additional_nsb = self.rng.poisson(mean_added_nsb, self.multiplicity * n_pixels)
         added_nsb_time = self.rng.uniform(-self.extra_samples * dt.to_value(u.ns),
                                           -self.extra_samples * dt.to_value(u.ns) + duration.to_value(u.ns),
@@ -680,13 +765,13 @@ def _tune_nsb_on_waveform_direct(self, waveform, tel_id, is_high_gain, subarray)
         dt = (1.0 / sampling_rate).to(u.s)
         duration = (self.extra_samples + n_samples) * dt
         t = np.arange(-self.extra_samples, n_samples) * dt.to_value(u.ns)
-        mean_added_nsb = (self.added_nsb_fraction * self.original_nsb[tel_id] * duration).to_value("")
+        mean_added_nsb = (self.added_nsb[tel_id] * duration).to_value("")
         additional_nsb = self.rng.poisson(mean_added_nsb, n_pixels)
         added_nsb_time = self.rng.uniform(-self.extra_samples * dt.to_value(u.ns),
                                           -self.extra_samples * dt.to_value(u.ns) + duration.to_value(u.ns),
                                           (n_pixels, max(additional_nsb)))
         added_nsb_amp = self.charge_spe_cumulative_pdf(self.rng.uniform(size=(n_pixels, max(additional_nsb))))
-        baseline_correction = (self.added_nsb_fraction * self.original_nsb[tel_id] * dt).to_value("")
+        baseline_correction = (self.added_nsb[tel_id] * dt).to_value("")
         waveform -= baseline_correction
         waveform += nsb_only_waveforms(
             time=t[self.extra_samples:],

From 99a2d9dab2af1b28fb3587e4920e6d0bf2e6f72e Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Wed, 31 Jul 2024 11:26:53 +0200
Subject: [PATCH 02/21] Adapted all calls of WaveformNsbTunner to new arguments

---
 lstchain/data/lstchain_lhfit_config.json    |  2 +-
 lstchain/data/lstchain_standard_config.json |  2 +-
 lstchain/image/tests/test_modifier.py       |  8 ++---
 lstchain/reco/r0_to_dl1.py                  | 39 ++++++++++++---------
 4 files changed, 27 insertions(+), 24 deletions(-)

diff --git a/lstchain/data/lstchain_lhfit_config.json b/lstchain/data/lstchain_lhfit_config.json
index c761c19937..05adda2285 100644
--- a/lstchain/data/lstchain_lhfit_config.json
+++ b/lstchain/data/lstchain_lhfit_config.json
@@ -279,7 +279,7 @@
   },
   "waveform_nsb_tuning":{
     "nsb_tuning": false,
-    "nsb_tuning_ratio": 0.5,
+    "nsb_tuning_rate_GHz": 0.15,
     "spe_location": null,
     "pre_computed_multiplicity": 10
   },
diff --git a/lstchain/data/lstchain_standard_config.json b/lstchain/data/lstchain_standard_config.json
index 74312da664..3935fee3f0 100644
--- a/lstchain/data/lstchain_standard_config.json
+++ b/lstchain/data/lstchain_standard_config.json
@@ -274,7 +274,7 @@
   },
   "waveform_nsb_tuning":{
     "nsb_tuning": false,
-    "nsb_tuning_ratio": 0.5,
+    "nsb_tuning_rate_GHz": 0.15,
     "spe_location": null,
     "pre_computed_multiplicity": 10
   },
diff --git a/lstchain/image/tests/test_modifier.py b/lstchain/image/tests/test_modifier.py
index a23c25eff2..b9beb2c213 100644
--- a/lstchain/image/tests/test_modifier.py
+++ b/lstchain/image/tests/test_modifier.py
@@ -87,7 +87,7 @@ def test_tune_nsb_on_waveform():
         [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
          [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]]
     )
-    added_nsb_fraction, original_nsb = 1.0, {1: 0.1 * u.GHz}
+    added_nsb_rate = {1: 0.1 * u.GHz}
     subarray = LSTEventSource.create_subarray(1)
     amplitude_HG = np.zeros(40)
     amplitude_LG = np.zeros(40)
@@ -106,8 +106,7 @@ def test_tune_nsb_on_waveform():
     charge_spe_cumulative_pdf = interp1d(spe_integral, spe[0], kind='cubic',
                                          bounds_error=False, fill_value=0.,
                                          assume_sorted=True)
-    nsb_tunner = WaveformNsbTunner(added_nsb_fraction,
-                                   original_nsb,
+    nsb_tunner = WaveformNsbTunner(added_nsb_rate,
                                    pulse_templates,
                                    charge_spe_cumulative_pdf,
                                    pre_computed_multiplicity=0)
@@ -118,8 +117,7 @@ def test_tune_nsb_on_waveform():
         [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
          [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]]
     )
-    nsb_tunner = WaveformNsbTunner(added_nsb_fraction,
-                                   original_nsb,
+    nsb_tunner = WaveformNsbTunner(added_nsb_rate,
                                    pulse_templates,
                                    charge_spe_cumulative_pdf,
                                    pre_computed_multiplicity=10)
diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py
index ce08cbaa64..c54a37c9ff 100644
--- a/lstchain/reco/r0_to_dl1.py
+++ b/lstchain/reco/r0_to_dl1.py
@@ -403,7 +403,7 @@ def r0_to_dl1(
             metadata=metadata,
         )
     nsb_tuning = False
-    nsb_tunner = None
+    nsb_tuner = None
     if 'waveform_nsb_tuning' in config.keys():
         nsb_tuning = config['waveform_nsb_tuning']['nsb_tuning']
         if nsb_tuning:
@@ -412,16 +412,18 @@ def r0_to_dl1(
                 pulse_templates = {tel_id: NormalizedPulseTemplate.load_from_eventsource(
                     subarray.tel[tel_id].camera.readout, resample=True)
                     for tel_id in config['source_config']['LSTEventSource']['allowed_tels']}
-                if 'nsb_tuning_ratio' in config['waveform_nsb_tuning'].keys():
+                if 'nsb_tuning_rate' in config['waveform_nsb_tuning'].keys():
                     # get value from config to possibly extract it beforehand on multiple files for averaging purposes
                     # or gain time
-                    nsb_tuning_ratio = config['waveform_nsb_tuning']['nsb_tuning_ratio']
+                    nsb_tuning_rate = config['waveform_nsb_tuning'][
+                        'nsb_tuning_rate_GHz']
                 else:
                     # extract the pedestal variance difference between the current MC file and the target data
                     # FIXME? fails for multiple telescopes
-                    nsb_tuning_ratio = calculate_required_additional_nsb(input_filename,
-                                                                         config['waveform_nsb_tuning']['target_data'],
-                                                                         config=config)[0]
+                    nsb_tuning_rate, _, _ = calculate_required_additional_nsb(
+                            input_filename,
+                            config['waveform_nsb_tuning']['target_data'],
+                            config=config)
                 spe_location = (config['waveform_nsb_tuning']['spe_location']
                                 or get_resource_path("data/SinglePhE_ResponseInPhE_expo2Gaus.dat"))
                 spe = np.loadtxt(spe_location).T
@@ -430,15 +432,18 @@ def r0_to_dl1(
                                                      bounds_error=False, fill_value=0.,
                                                      assume_sorted=True)
                 pre_computed_multiplicity = config['waveform_nsb_tuning'].get('pre_computed_multiplicity', 10)
-                logger.info('Tuning NSB on MC waveform from '
-                            + str(np.asarray(nsb_original))
-                            + ' to {0:d}%'.format(int(nsb_tuning_ratio * 100 + 100.5))
-                            + ' for telescopes ids ' + str(config['source_config']['LSTEventSource']['allowed_tels']))
-                nsb_tunner = WaveformNsbTunner(nsb_tuning_ratio,
-                                               nsb_original,
-                                               pulse_templates,
-                                               charge_spe_cumulative_pdf,
-                                               pre_computed_multiplicity)
+
+                allowed_tels = config['source_config']['LSTEventSource'][
+                    'allowed_tels']
+                logger.info(f'Tuning NSB on MC waveform by adding ')
+                logger.info(f'{nsb_tuning_rate:.3f} GHz for telescope ids:')
+                logger.info(f'{allowed_tels}')
+
+                nsb_per_tel = np.ones(np.max(allowed_tels)+1) * nsb_tuning_rate
+                nsb_tuner = WaveformNsbTunner(nsb_per_tel * u.GHz,
+                                              pulse_templates,
+                                              charge_spe_cumulative_pdf,
+                                              pre_computed_multiplicity)
             else:
                 logger.warning('NSB tuning on waveform active in config but file is real data, option will be ignored')
                 nsb_tuning = False
@@ -451,7 +456,7 @@ def r0_to_dl1(
             tmp_source = EventSource(input_url=input_filename,
                                      config=Config(config["source_config"]))
             if is_simu:
-                lhfit_fitter.get_ped_from_true_signal_less(tmp_source, nsb_tunner)
+                lhfit_fitter.get_ped_from_true_signal_less(tmp_source, nsb_tuner)
             else:
                 lhfit_fitter.get_ped_from_interleaved(tmp_source)
             del tmp_source
@@ -570,7 +575,7 @@ def r0_to_dl1(
                     waveform = event.r1.tel[tel_id].waveform
                     selected_gains = event.r1.tel[tel_id].selected_gain_channel
                     mask_high = selected_gains == 0
-                    nsb_tunner.tune_nsb_on_waveform(waveform, tel_id, mask_high, subarray)
+                    nsb_tuner.tune_nsb_on_waveform(waveform, tel_id, mask_high, subarray)
 
             # create image for all events
             r1_dl1_calibrator(event)

From bf27816e7b336281fae7dc210a94b9989836f250 Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Wed, 31 Jul 2024 12:52:16 +0200
Subject: [PATCH 03/21] Solved a number of issues (wrong variable names etc)

---
 lstchain/image/modifier.py | 12 ++++--------
 lstchain/reco/r0_to_dl1.py |  1 -
 2 files changed, 4 insertions(+), 9 deletions(-)

diff --git a/lstchain/image/modifier.py b/lstchain/image/modifier.py
index 6676400b71..390ad15206 100644
--- a/lstchain/image/modifier.py
+++ b/lstchain/image/modifier.py
@@ -12,7 +12,6 @@
 from traitlets.config import Config
 
 from lstchain.io import standard_config, read_configuration_file
-from lstchain.io.config import read_configuration_file
 from lstchain.reco.reconstructorCC import nsb_only_waveforms
 from lstchain.data import NormalizedPulseTemplate
 from lstchain.io.io import get_resource_path
@@ -409,7 +408,6 @@ def calculate_noise_parameters(simtel_filename, data_dl1_filename,
     return extra_noise_in_dim_pixels, extra_bias_in_dim_pixels, \
         extra_noise_in_bright_pixels
 
-
 def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config=None):
     """
     Calculates the additional NSB needed in the MC waveforms
@@ -473,7 +471,7 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
     # and "not too bright" pixels (the idea is to exclude stars, we want to
     # estimate the overall "diffuse" NSB level)
     interleaved_ped = image_params['event_type'] == EventType.SKY_PEDESTAL.value
-    ped_pixq = camera['image'][interleaved_ped]
+    ped_pixq = camera_images['image'][interleaved_ped]
     ped_meanpixq = np.mean(ped_pixq, axis=1)
     ped_stdpixq = np.std(ped_pixq, axis=1)
     median_ped_meanpixq = np.median(ped_meanpixq[good_pixels])
@@ -489,9 +487,7 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
     # Event reader for simtel file:
     mc_reader = EventSource(input_url=simtel_filename, config=Config(config))
 
-    subarray = mcreader.subarray
-    sampling_rate = subarray.tel[1].camera.readout.sampling_rate
-    dt = (1.0 / sampling_rate).to(u.ns)
+    subarray = mc_reader.subarray
 
     # Get the single-pe response fluctuations:
     spe_location = (config['waveform_nsb_tuning']['spe_location']
@@ -538,7 +534,7 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
         # the future case in which there are more telescopes, but for now we
         # just create it with a single value for tel_id
         nsb = np.zeros(tel_id+1)
-        nsb[tel_id] = nsbrate
+        nsb[tel_id] = nsb_rate
         nsb_tuner.append(WaveformNsbTunner(nsb * u.GHz, pulse_templates,
                                            charge_spe_cumulative_pdf, 10))
     # last argument means it will precompute 10 * 1855 (pixels) * 2 (gains)
@@ -598,7 +594,7 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
                               charge_spe_cumulative_pdf, 10)
     final_mc_qped = []
     numevents = 0
-    for event in mcreader:
+    for event in mc_reader:
         if tel_id not in event.trigger.tels_with_trigger:
             continue
         numevents += 1
diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py
index c54a37c9ff..7f4153ac4b 100644
--- a/lstchain/reco/r0_to_dl1.py
+++ b/lstchain/reco/r0_to_dl1.py
@@ -408,7 +408,6 @@ def r0_to_dl1(
         nsb_tuning = config['waveform_nsb_tuning']['nsb_tuning']
         if nsb_tuning:
             if is_simu:
-                nsb_original = extract_simulation_nsb(input_filename)
                 pulse_templates = {tel_id: NormalizedPulseTemplate.load_from_eventsource(
                     subarray.tel[tel_id].camera.readout, resample=True)
                     for tel_id in config['source_config']['LSTEventSource']['allowed_tels']}

From b7237a6fc21a1fb3788ef11303014cdfba8b6fd1 Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Wed, 31 Jul 2024 13:44:54 +0200
Subject: [PATCH 04/21] Fix logger string. Remove unused import.

---
 lstchain/reco/r0_to_dl1.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py
index 7f4153ac4b..123daf17f4 100644
--- a/lstchain/reco/r0_to_dl1.py
+++ b/lstchain/reco/r0_to_dl1.py
@@ -59,7 +59,7 @@
     write_subarray_tables
 )
 
-from ..io.io import add_column_table, extract_simulation_nsb, dl1_params_lstcam_key, get_resource_path
+from ..io.io import add_column_table, dl1_params_lstcam_key, get_resource_path
 from ..io.lstcontainers import ExtraImageInfo, DL1MonitoringEventIndexContainer
 from ..paths import parse_r0_filename, run_to_dl1_filename, r0_to_dl1_filename
 from ..visualization.plot_reconstructor import plot_debug
@@ -434,7 +434,7 @@ def r0_to_dl1(
 
                 allowed_tels = config['source_config']['LSTEventSource'][
                     'allowed_tels']
-                logger.info(f'Tuning NSB on MC waveform by adding ')
+                logger.info('Tuning NSB on MC waveform by adding ')
                 logger.info(f'{nsb_tuning_rate:.3f} GHz for telescope ids:')
                 logger.info(f'{allowed_tels}')
 

From 7a2a68686cc3767d0a2c842a680858921e6d4e5b Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Wed, 31 Jul 2024 14:26:44 +0200
Subject: [PATCH 05/21] Fixed wrong averaging axis

---
 lstchain/image/modifier.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/lstchain/image/modifier.py b/lstchain/image/modifier.py
index 390ad15206..7ba6e6d2ab 100644
--- a/lstchain/image/modifier.py
+++ b/lstchain/image/modifier.py
@@ -472,8 +472,8 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
     # estimate the overall "diffuse" NSB level)
     interleaved_ped = image_params['event_type'] == EventType.SKY_PEDESTAL.value
     ped_pixq = camera_images['image'][interleaved_ped]
-    ped_meanpixq = np.mean(ped_pixq, axis=1)
-    ped_stdpixq = np.std(ped_pixq, axis=1)
+    ped_meanpixq = np.mean(ped_pixq, axis=0)
+    ped_stdpixq = np.std(ped_pixq, axis=0)
     median_ped_meanpixq = np.median(ped_meanpixq[good_pixels])
     # Exclude the brightest pixels, which may be affected by stars:
     too_bright = (ped_meanpixq > median_ped_meanpixq + 3 * np.std(ped_meanpixq))

From b3f07d19a77db16c43b55a8bbcf5d7b7295e5a62 Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Wed, 31 Jul 2024 14:47:04 +0200
Subject: [PATCH 06/21] Wrong indent!

---
 lstchain/image/modifier.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/lstchain/image/modifier.py b/lstchain/image/modifier.py
index 7ba6e6d2ab..5bb4922f59 100644
--- a/lstchain/image/modifier.py
+++ b/lstchain/image/modifier.py
@@ -574,11 +574,11 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
             modified_integrated_charge[ii + 1].extend(
                     event.dl1.tel[1].image[pedmask])
 
-        modified_integrated_charge = np.array(modified_integrated_charge)
-        # Fit the total added NSB rate vs. the average pixel charge:
-        params, _ = curve_fit(custom_function,
-                              np.mean(modified_integrated_charge, axis=1),
-                              total_added_nsb, p0=[-0.2, 0.2, 2])
+    modified_integrated_charge = np.array(modified_integrated_charge)
+    # Fit the total added NSB rate vs. the average pixel charge:
+    params, _ = curve_fit(custom_function,
+                          np.mean(modified_integrated_charge, axis=1),
+                          total_added_nsb, p0=[-0.2, 0.2, 2])
 
     # Obtain the right rate of NSB to be added to MC so that it matches the
     # data:

From cbbaaa59b3af90867c9a9007b720406cded33691 Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Wed, 31 Jul 2024 15:03:17 +0200
Subject: [PATCH 07/21] Replaced array by dict in arguments

---
 lstchain/image/modifier.py | 22 +++++++++++-----------
 1 file changed, 11 insertions(+), 11 deletions(-)

diff --git a/lstchain/image/modifier.py b/lstchain/image/modifier.py
index 5bb4922f59..fb49020017 100644
--- a/lstchain/image/modifier.py
+++ b/lstchain/image/modifier.py
@@ -530,12 +530,12 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
 
 
     for nsb_rate in np.diff(total_added_nsb):
-        # Create an array to put the NSB value for tel_id. It is an array for
-        # the future case in which there are more telescopes, but for now we
-        # just create it with a single value for tel_id
-        nsb = np.zeros(tel_id+1)
-        nsb[tel_id] = nsb_rate
-        nsb_tuner.append(WaveformNsbTunner(nsb * u.GHz, pulse_templates,
+        # Create a dict to put the NSB value for each tel_id. It is an array for
+        # the future case in which there are more telescopes, but for now it
+        # just creates it with a single value, for tel_id
+        nsb = {tel_id: nsb_rate * u.GHz for tel_id in
+               config['source_config']['LSTEventSource']['allowed_tels']}
+        nsb_tuner.append(WaveformNsbTunner(nsb, pulse_templates,
                                            charge_spe_cumulative_pdf, 10))
     # last argument means it will precompute 10 * 1855 (pixels) * 2 (gains)
     # noise waveforms to pick from during the actual introduction of the noise
@@ -588,9 +588,9 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
 
     # Now open the MC file again and test that the tuning was successful:
     mc_reader = EventSource(input_url=simtel_filename, config=Config(config))
-    nsb = np.zeros(tel_id+1)
-    nsb[tel_id] = extra_nsb
-    tuner = WaveformNsbTunner(nsb * u.GHz , pulse_templates,
+    nsb = {tel_id: extra_nsb * u.GHz for tel_id in
+           config['source_config']['LSTEventSource']['allowed_tels']}
+    tuner = WaveformNsbTunner(nsb, pulse_templates,
                               charge_spe_cumulative_pdf, 10)
     final_mc_qped = []
     numevents = 0
@@ -632,8 +632,8 @@ def __init__(self, added_nsb, pulse_template, charge_spe_cumulative_pdf,
         """
         Parameters
         ----------
-        added_nsb: array [`float`,`astropy.units.Quantity`]
-            NSB frequency (GHz) to be added to the waveforms (per telescope)
+        added_nsb: dict [`float`,`astropy.units.Quantity`]
+            NSB frequency (GHz) to be added to the waveforms (per tel_id)
         pulse_template: `dict`[int,`lstchain.data.NormalizedPulseTemplate`]
             Single photo-electron pulse template per telescope
         charge_spe_cumulative_pdf: `scipy.interpolate.interp1d`

From 0c96287fd7f812987f83860f0f2f1dc81ba541c0 Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Wed, 31 Jul 2024 15:18:27 +0200
Subject: [PATCH 08/21] Small refactoring

---
 lstchain/image/modifier.py | 77 ++++++++++++++++++++++----------------
 1 file changed, 44 insertions(+), 33 deletions(-)

diff --git a/lstchain/image/modifier.py b/lstchain/image/modifier.py
index fb49020017..dd3427c695 100644
--- a/lstchain/image/modifier.py
+++ b/lstchain/image/modifier.py
@@ -408,39 +408,10 @@ def calculate_noise_parameters(simtel_filename, data_dl1_filename,
     return extra_noise_in_dim_pixels, extra_bias_in_dim_pixels, \
         extra_noise_in_bright_pixels
 
-def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config=None):
-    """
-    Calculates the additional NSB needed in the MC waveforms
-    to match a real data DL1 file
-
-    Parameters
-    ----------
-    simtel_filename: a simtel file containing showers, from the production
-    (same NSB and telescope settings) as the one on which the correction will
-    be applied. It must contain pixel-wise info on true number of p.e.'s from
-    C-photons (will be used to identify pixels which only contain noise).
-    data_dl1_filename: a real data DL1 file (processed with calibration
-    settings corresponding to those with which the MC is to be processed).
-    It must contain calibrated images, i.e. "DL1a" data. This file has the
-    "target" NSB which we want to have in the MC files, for better
-    agreement of data and simulations.
-    config: configuration containing the calibration
-    settings used for processing both the data and the MC files above
-
-    Returns
-    -------
-    extra_nsb: additional NSB rate in absolute units, p.e./ns (a.k.a. "GHz"),
-    that has to be added in the MC to match the real data
-    data_ped_variance: Pedestal variance from data
-    mc_ped_variance: Pedestal variance from MC
 
+def get_pix_median_charges(data_dl1_filename, event_type):
+    """
     """
-
-    log.setLevel(logging.INFO)
-
-    if config is None:
-        config = standard_config
-
     # Real data DL1 tables:
     camera_images = read_table(data_dl1_filename,
                                '/dl1/event/telescope/image/LST_LSTCam')
@@ -470,7 +441,7 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
     # Now compute the average pixel charge in real pedestal events, for good
     # and "not too bright" pixels (the idea is to exclude stars, we want to
     # estimate the overall "diffuse" NSB level)
-    interleaved_ped = image_params['event_type'] == EventType.SKY_PEDESTAL.value
+    interleaved_ped = image_params['event_type'] == event_type.value
     ped_pixq = camera_images['image'][interleaved_ped]
     ped_meanpixq = np.mean(ped_pixq, axis=0)
     ped_stdpixq = np.std(ped_pixq, axis=0)
@@ -479,10 +450,50 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
     too_bright = (ped_meanpixq > median_ped_meanpixq + 3 * np.std(ped_meanpixq))
     good_pixels &= ~too_bright
     log.info(f'Good and not too bright pixels: {good_pixels.sum()}')
-    # Recompyte median:
+    # Recompute median:
     median_ped_meanpixq = np.median(ped_meanpixq[good_pixels])
     median_ped_stdpixq = np.median(ped_stdpixq[good_pixels])
 
+    return median_ped_meanpixq, median_ped_stdpixq
+
+
+
+def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config=None):
+    """
+    Calculates the additional NSB needed in the MC waveforms
+    to match a real data DL1 file
+
+    Parameters
+    ----------
+    simtel_filename: a simtel file containing showers, from the production
+    (same NSB and telescope settings) as the one on which the correction will
+    be applied. It must contain pixel-wise info on true number of p.e.'s from
+    C-photons (will be used to identify pixels which only contain noise).
+    data_dl1_filename: a real data DL1 file (processed with calibration
+    settings corresponding to those with which the MC is to be processed).
+    It must contain calibrated images, i.e. "DL1a" data. This file has the
+    "target" NSB which we want to have in the MC files, for better
+    agreement of data and simulations.
+    config: configuration containing the calibration
+    settings used for processing both the data and the MC files above
+
+    Returns
+    -------
+    extra_nsb: additional NSB rate in absolute units, p.e./ns (a.k.a. "GHz"),
+    that has to be added in the MC to match the real data
+    data_ped_variance: Pedestal variance from data
+    mc_ped_variance: Pedestal variance from MC
+
+    """
+
+    log.setLevel(logging.INFO)
+
+    if config is None:
+        config = standard_config
+
+    median_ped_meanpixq, median_ped_stdpixq = get_pix_median_charges(
+            data_dl1_filename, EventType.SKY_PEDESTAL)
+
     # Now we process the Monte Carlo:
     # Event reader for simtel file:
     mc_reader = EventSource(input_url=simtel_filename, config=Config(config))

From 481d90ce11711987c5cdbb270c51961d4c6abc31 Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Wed, 31 Jul 2024 15:20:15 +0200
Subject: [PATCH 09/21] Removed unused variable

---
 lstchain/image/modifier.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lstchain/image/modifier.py b/lstchain/image/modifier.py
index dd3427c695..f89dd3f30b 100644
--- a/lstchain/image/modifier.py
+++ b/lstchain/image/modifier.py
@@ -421,7 +421,7 @@ def get_pix_median_charges(data_dl1_filename, event_type):
 
     # All events in a DL1 file should correspond to the same telescope,
     # just take it from the first event:
-    tel_id = image_params['tel_id'][0]
+    # tel_id = image_params['tel_id'][0]
 
     data_dl1_calibration = read_table(data_dl1_filename,
                                       '/dl1/event/telescope/monitoring/calibration')

From 6b7d05c4531f1f162d0770c7caa37e9e173ccbc4 Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Wed, 31 Jul 2024 16:21:51 +0200
Subject: [PATCH 10/21] fix missing tel_id

---
 lstchain/image/modifier.py | 28 +++++++++++++++++++++++++---
 1 file changed, 25 insertions(+), 3 deletions(-)

diff --git a/lstchain/image/modifier.py b/lstchain/image/modifier.py
index f89dd3f30b..5f5490d3c4 100644
--- a/lstchain/image/modifier.py
+++ b/lstchain/image/modifier.py
@@ -411,7 +411,26 @@ def calculate_noise_parameters(simtel_filename, data_dl1_filename,
 
 def get_pix_median_charges(data_dl1_filename, event_type):
     """
+    Obtain from a DL1 real data file (containing camera images, i.e. DL1a)
+    the median (accross camera) of mean (and std dev) of pixel charge for
+    events of type event_type. The medians are obtained using only healthy
+    pixels and excluding too bright outliers (we want to get the typical
+    pixel values, unaffected by the few stars in the FoV)
+
+    Parameters
+    ----------
+    data_dl1_filename: DL1 file name
+    event_type (ctapipe.containers.EventType): type of events, e.g. pedestals
+
+    Returns
+    -------
+    median_ped_meanpixq (p.e.)
+    median_ped_stdpixq  (p.e.)  Median values for healthy and not too
+    noisy pixels
+    tel_id : telescope id to which the events in the file correspond (we assume
+    for now that a given DL1 file contains data from just one telescope)
     """
+
     # Real data DL1 tables:
     camera_images = read_table(data_dl1_filename,
                                '/dl1/event/telescope/image/LST_LSTCam')
@@ -421,7 +440,7 @@ def get_pix_median_charges(data_dl1_filename, event_type):
 
     # All events in a DL1 file should correspond to the same telescope,
     # just take it from the first event:
-    # tel_id = image_params['tel_id'][0]
+    tel_id = image_params['tel_id'][0]
 
     data_dl1_calibration = read_table(data_dl1_filename,
                                       '/dl1/event/telescope/monitoring/calibration')
@@ -454,7 +473,7 @@ def get_pix_median_charges(data_dl1_filename, event_type):
     median_ped_meanpixq = np.median(ped_meanpixq[good_pixels])
     median_ped_stdpixq = np.median(ped_stdpixq[good_pixels])
 
-    return median_ped_meanpixq, median_ped_stdpixq
+    return median_ped_meanpixq, median_ped_stdpixq, tel_id
 
 
 
@@ -491,7 +510,10 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
     if config is None:
         config = standard_config
 
-    median_ped_meanpixq, median_ped_stdpixq = get_pix_median_charges(
+    # Obtain the median (accross camera) of mean (and satd dev of) pixel charge
+    # taking into account only healthy pixels, and excluding outliers.
+    # We also get the tel_id to which the DL1 file corresponds:
+    median_ped_meanpixq, median_ped_stdpixq, tel_id = get_pix_median_charges(
             data_dl1_filename, EventType.SKY_PEDESTAL)
 
     # Now we process the Monte Carlo:

From 0f75e8be43f208834268a1022961eb0574d49abe Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Wed, 31 Jul 2024 16:57:00 +0200
Subject: [PATCH 11/21] Attempt to fix tests

---
 lstchain/scripts/lstchain_tune_nsb_waveform.py  |  4 ++--
 lstchain/scripts/tests/test_lstchain_scripts.py | 10 +++++-----
 2 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/lstchain/scripts/lstchain_tune_nsb_waveform.py b/lstchain/scripts/lstchain_tune_nsb_waveform.py
index 1a1f8ee0f4..615dbba3fc 100644
--- a/lstchain/scripts/lstchain_tune_nsb_waveform.py
+++ b/lstchain/scripts/lstchain_tune_nsb_waveform.py
@@ -82,13 +82,13 @@ def main():
 
     config = read_configuration_file(args.config)
 
-    nsb_correction_ratio, data_ped_variance, mc_ped_variance = \
+    extra_nsb_rate, data_ped_variance, mc_ped_variance = \
         calculate_required_additional_nsb(args.input_mc, args.input_data,
                                           config=Config(config))
 
     dict_nsb = {
         "nsb_tuning": True,
-        "nsb_tuning_ratio": np.round(nsb_correction_ratio, decimals=2),
+        "nsb_tuning_rate": np.round(extra_nsb_rate, decimals=3),
         "spe_location": str(get_resource_path("data/SinglePhE_ResponseInPhE_expo2Gaus.dat"))
     }
 
diff --git a/lstchain/scripts/tests/test_lstchain_scripts.py b/lstchain/scripts/tests/test_lstchain_scripts.py
index b7354ce2a5..0311225fc5 100644
--- a/lstchain/scripts/tests/test_lstchain_scripts.py
+++ b/lstchain/scripts/tests/test_lstchain_scripts.py
@@ -236,16 +236,16 @@ def test_validity_tune_nsb(tune_nsb):
 
 def test_validity_tune_nsb_waveform(tune_nsb_waveform):
     """
-    The resulting nsb_tuning_ratio value of -1 expected in this test is meaningless
-    because the input data do not allow a full test of the functionality.
-    This test is only a formal check that the script runs.
+    The resulting nsb_tuning_rate value of -1 expected in this test is
+    meaningless because the input data do not allow a full test of the
+    functionality. This test is only a formal check that the script runs.
     """
     output_lines = tune_nsb_waveform.stdout.splitlines()
     for line in output_lines:
         if '"nsb_tuning"' in line:
             assert line == '  "nsb_tuning": true,'
-        if '"nsb_tuning_ratio"' in line:
-            assert line == '  "nsb_tuning_ratio": -1.0,'
+        if '"nsb_tuning_rate"' in line:
+            assert line == '  "nsb_tuning_rate": -1.0,'
         if '"spe_location"' in line:
             assert line == f'  "spe_location": "{get_resource_path("data/SinglePhE_ResponseInPhE_expo2Gaus.dat")}"'
 

From 3cb7b2b0abde021813db54fd117d7e114c4b1a0d Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Wed, 31 Jul 2024 17:20:43 +0200
Subject: [PATCH 12/21] Correct config key

---
 lstchain/reco/r0_to_dl1.py                     | 9 ++++++---
 lstchain/scripts/lstchain_tune_nsb_waveform.py | 2 +-
 2 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py
index 123daf17f4..52ff541fc0 100644
--- a/lstchain/reco/r0_to_dl1.py
+++ b/lstchain/reco/r0_to_dl1.py
@@ -411,7 +411,8 @@ def r0_to_dl1(
                 pulse_templates = {tel_id: NormalizedPulseTemplate.load_from_eventsource(
                     subarray.tel[tel_id].camera.readout, resample=True)
                     for tel_id in config['source_config']['LSTEventSource']['allowed_tels']}
-                if 'nsb_tuning_rate' in config['waveform_nsb_tuning'].keys():
+                if 'nsb_tuning_rate_GHz' in config[
+                    'waveform_nsb_tuning'].keys():
                     # get value from config to possibly extract it beforehand on multiple files for averaging purposes
                     # or gain time
                     nsb_tuning_rate = config['waveform_nsb_tuning'][
@@ -438,8 +439,10 @@ def r0_to_dl1(
                 logger.info(f'{nsb_tuning_rate:.3f} GHz for telescope ids:')
                 logger.info(f'{allowed_tels}')
 
-                nsb_per_tel = np.ones(np.max(allowed_tels)+1) * nsb_tuning_rate
-                nsb_tuner = WaveformNsbTunner(nsb_per_tel * u.GHz,
+                nsb_per_tel = {tel_id: nsb_tuning_rate * u.GHz for tel_id in
+                               allowed_tels}
+
+                nsb_tuner = WaveformNsbTunner(nsb_per_tel,
                                               pulse_templates,
                                               charge_spe_cumulative_pdf,
                                               pre_computed_multiplicity)
diff --git a/lstchain/scripts/lstchain_tune_nsb_waveform.py b/lstchain/scripts/lstchain_tune_nsb_waveform.py
index 615dbba3fc..a91e940827 100644
--- a/lstchain/scripts/lstchain_tune_nsb_waveform.py
+++ b/lstchain/scripts/lstchain_tune_nsb_waveform.py
@@ -88,7 +88,7 @@ def main():
 
     dict_nsb = {
         "nsb_tuning": True,
-        "nsb_tuning_rate": np.round(extra_nsb_rate, decimals=3),
+        "nsb_tuning_rate_GHz": np.round(extra_nsb_rate, decimals=3),
         "spe_location": str(get_resource_path("data/SinglePhE_ResponseInPhE_expo2Gaus.dat"))
     }
 

From 21d725be5f1078b1812641c2983534f72406a184 Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Wed, 31 Jul 2024 17:44:19 +0200
Subject: [PATCH 13/21] Removed unnecesary ".keys()"

---
 lstchain/reco/r0_to_dl1.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py
index 52ff541fc0..492d462a7a 100644
--- a/lstchain/reco/r0_to_dl1.py
+++ b/lstchain/reco/r0_to_dl1.py
@@ -404,7 +404,7 @@ def r0_to_dl1(
         )
     nsb_tuning = False
     nsb_tuner = None
-    if 'waveform_nsb_tuning' in config.keys():
+    if 'waveform_nsb_tuning' in config:
         nsb_tuning = config['waveform_nsb_tuning']['nsb_tuning']
         if nsb_tuning:
             if is_simu:
@@ -412,7 +412,7 @@ def r0_to_dl1(
                     subarray.tel[tel_id].camera.readout, resample=True)
                     for tel_id in config['source_config']['LSTEventSource']['allowed_tels']}
                 if 'nsb_tuning_rate_GHz' in config[
-                    'waveform_nsb_tuning'].keys():
+                    'waveform_nsb_tuning']:
                     # get value from config to possibly extract it beforehand on multiple files for averaging purposes
                     # or gain time
                     nsb_tuning_rate = config['waveform_nsb_tuning'][
@@ -451,7 +451,7 @@ def r0_to_dl1(
                 nsb_tuning = False
 
     lhfit_fitter = None
-    if 'lh_fit_config' in config.keys():
+    if 'lh_fit_config' in config:
         lhfit_fitter_config = {'TimeWaveformFitter': config['lh_fit_config']}
         lhfit_fitter = TimeWaveformFitter(subarray=subarray, config=Config(lhfit_fitter_config))
         if lhfit_fitter_config['TimeWaveformFitter']['use_interleaved']:

From 52c8dd275f32a5a164768407e49dec4999aba48c Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Wed, 31 Jul 2024 18:42:27 +0200
Subject: [PATCH 14/21] Improved comment. Made test (with slightly random
 result) more tolerant

---
 lstchain/image/modifier.py            | 2 +-
 lstchain/image/tests/test_modifier.py | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/lstchain/image/modifier.py b/lstchain/image/modifier.py
index 5f5490d3c4..9f18cc1f11 100644
--- a/lstchain/image/modifier.py
+++ b/lstchain/image/modifier.py
@@ -501,7 +501,7 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
     extra_nsb: additional NSB rate in absolute units, p.e./ns (a.k.a. "GHz"),
     that has to be added in the MC to match the real data
     data_ped_variance: Pedestal variance from data
-    mc_ped_variance: Pedestal variance from MC
+    mc_ped_variance: Pedestal variance from MC, AFTER THE TUNING!
 
     """
 
diff --git a/lstchain/image/tests/test_modifier.py b/lstchain/image/tests/test_modifier.py
index b9beb2c213..a66bef1bf1 100644
--- a/lstchain/image/tests/test_modifier.py
+++ b/lstchain/image/tests/test_modifier.py
@@ -71,8 +71,8 @@ def test_calculate_required_additional_nsb(mc_gamma_testfile, observed_dl1_files
         observed_dl1_files["dl1_file1"]
     )
     assert np.isclose(data_ped_variance, 0.0, atol=0.1)
-    assert np.isclose(mc_ped_variance, 3.11, atol=0.01)
-    assert np.isclose(extra_nsb, -1.0)
+    assert np.isclose(mc_ped_variance, 3.2, atol=0.2)
+    assert np.isclose(extra_nsb, 0.22, atol=0.03)
 
 
 def test_tune_nsb_on_waveform():

From 38590fc5dae16ae5aec68dbe2cace923e418f394 Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Thu, 1 Aug 2024 09:47:14 +0200
Subject: [PATCH 15/21] Added comments. Don't allow extra_nsb to be <0

---
 lstchain/image/modifier.py            | 4 ++++
 lstchain/image/tests/test_modifier.py | 2 ++
 2 files changed, 6 insertions(+)

diff --git a/lstchain/image/modifier.py b/lstchain/image/modifier.py
index 9f18cc1f11..b5a1642ac8 100644
--- a/lstchain/image/modifier.py
+++ b/lstchain/image/modifier.py
@@ -618,6 +618,10 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
     extra_nsb = custom_function(median_ped_meanpixq,
                                 params[0], params[1], params[2])
 
+    # Since the level of noise in MC is low, it should hardly ever happen
+    # that extra_nsb is negative, but just in case:
+    if extra_nsb < 0:
+        extra_nsb = 0
 
     # Now open the MC file again and test that the tuning was successful:
     mc_reader = EventSource(input_url=simtel_filename, config=Config(config))
diff --git a/lstchain/image/tests/test_modifier.py b/lstchain/image/tests/test_modifier.py
index a66bef1bf1..c549c6d46c 100644
--- a/lstchain/image/tests/test_modifier.py
+++ b/lstchain/image/tests/test_modifier.py
@@ -70,6 +70,8 @@ def test_calculate_required_additional_nsb(mc_gamma_testfile, observed_dl1_files
         mc_gamma_testfile,
         observed_dl1_files["dl1_file1"]
     )
+    # Mean pixel charge variance in pedestals is 0 because there is only one
+    # pedestal event in the test file
     assert np.isclose(data_ped_variance, 0.0, atol=0.1)
     assert np.isclose(mc_ped_variance, 3.2, atol=0.2)
     assert np.isclose(extra_nsb, 0.22, atol=0.03)

From e5c2446bf4536f108d9225c34b38201a2290ade4 Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Thu, 1 Aug 2024 11:49:47 +0200
Subject: [PATCH 16/21] Fixed test

---
 lstchain/image/tests/test_modifier.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/lstchain/image/tests/test_modifier.py b/lstchain/image/tests/test_modifier.py
index c549c6d46c..5537fa0fe0 100644
--- a/lstchain/image/tests/test_modifier.py
+++ b/lstchain/image/tests/test_modifier.py
@@ -73,8 +73,8 @@ def test_calculate_required_additional_nsb(mc_gamma_testfile, observed_dl1_files
     # Mean pixel charge variance in pedestals is 0 because there is only one
     # pedestal event in the test file
     assert np.isclose(data_ped_variance, 0.0, atol=0.1)
-    assert np.isclose(mc_ped_variance, 3.2, atol=0.2)
-    assert np.isclose(extra_nsb, 0.22, atol=0.03)
+    assert np.isclose(mc_ped_variance, 2.3, atol=0.2)
+    assert np.isclose(extra_nsb, 0.08, atol=0.02)
 
 
 def test_tune_nsb_on_waveform():

From db0bdf72e818a19a71099123266cbada333c19ce Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo <moralejo@ifae.es>
Date: Thu, 1 Aug 2024 13:22:44 +0200
Subject: [PATCH 17/21] Update lstchain/image/modifier.py

Co-authored-by: Daniel Morcuende <dmorcuende@iaa.es>
---
 lstchain/image/modifier.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lstchain/image/modifier.py b/lstchain/image/modifier.py
index b5a1642ac8..03e7da4ab4 100644
--- a/lstchain/image/modifier.py
+++ b/lstchain/image/modifier.py
@@ -510,7 +510,7 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
     if config is None:
         config = standard_config
 
-    # Obtain the median (accross camera) of mean (and satd dev of) pixel charge
+    # Obtain the median (across camera) of mean (and std dev of) pixel charge
     # taking into account only healthy pixels, and excluding outliers.
     # We also get the tel_id to which the DL1 file corresponds:
     median_ped_meanpixq, median_ped_stdpixq, tel_id = get_pix_median_charges(

From 15f6f58369107726b31d2704dd891677569c837a Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo <moralejo@ifae.es>
Date: Thu, 1 Aug 2024 13:23:11 +0200
Subject: [PATCH 18/21] Update lstchain/image/modifier.py

Co-authored-by: Daniel Morcuende <dmorcuende@iaa.es>
---
 lstchain/image/modifier.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lstchain/image/modifier.py b/lstchain/image/modifier.py
index 03e7da4ab4..56d64644f4 100644
--- a/lstchain/image/modifier.py
+++ b/lstchain/image/modifier.py
@@ -412,7 +412,7 @@ def calculate_noise_parameters(simtel_filename, data_dl1_filename,
 def get_pix_median_charges(data_dl1_filename, event_type):
     """
     Obtain from a DL1 real data file (containing camera images, i.e. DL1a)
-    the median (accross camera) of mean (and std dev) of pixel charge for
+    the median (across camera) of mean (and std dev) of pixel charge for
     events of type event_type. The medians are obtained using only healthy
     pixels and excluding too bright outliers (we want to get the typical
     pixel values, unaffected by the few stars in the FoV)

From 933df16369e79b4f3241a0af779b973ad1e6992b Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo <moralejo@ifae.es>
Date: Thu, 1 Aug 2024 13:23:49 +0200
Subject: [PATCH 19/21] Use max() for more readable code

Co-authored-by: Daniel Morcuende <dmorcuende@iaa.es>
---
 lstchain/image/modifier.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/lstchain/image/modifier.py b/lstchain/image/modifier.py
index 56d64644f4..4c342e398b 100644
--- a/lstchain/image/modifier.py
+++ b/lstchain/image/modifier.py
@@ -620,8 +620,7 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
 
     # Since the level of noise in MC is low, it should hardly ever happen
     # that extra_nsb is negative, but just in case:
-    if extra_nsb < 0:
-        extra_nsb = 0
+    extra_nsb = max(extra_nsb, 0)
 
     # Now open the MC file again and test that the tuning was successful:
     mc_reader = EventSource(input_url=simtel_filename, config=Config(config))

From fbd2c1979b25c16d59fef1bdb5727b4771b240b1 Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo <moralejo@ifae.es>
Date: Thu, 1 Aug 2024 13:26:49 +0200
Subject: [PATCH 20/21] Typo

Co-authored-by: Daniel Morcuende <dmorcuende@iaa.es>
---
 lstchain/image/modifier.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lstchain/image/modifier.py b/lstchain/image/modifier.py
index 4c342e398b..89036257d8 100644
--- a/lstchain/image/modifier.py
+++ b/lstchain/image/modifier.py
@@ -597,7 +597,7 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
         modified_integrated_charge[0].extend(event.dl1.tel[tel_id].image[
                                                  pedmask])
 
-        # Now add the differentlevels of NSB and recompute charges:
+        # Now add the different levels of NSB and recompute charges:
         for ii, tuner in enumerate(nsb_tuner[1:]):
             waveform = event.r1.tel[tel_id].waveform
 

From 20d04bbbceff1d9487ced70856f447b56f5d811e Mon Sep 17 00:00:00 2001
From: Abelardo Moralejo Olaizola <moralejo@ifae.es>
Date: Thu, 1 Aug 2024 13:35:18 +0200
Subject: [PATCH 21/21] WaveformNsbTunner => WaveformNsbTuner

---
 lstchain/image/modifier.py            | 16 +++++++++-------
 lstchain/image/tests/test_modifier.py | 18 +++++++++---------
 lstchain/reco/r0_to_dl1.py            | 10 +++++-----
 3 files changed, 23 insertions(+), 21 deletions(-)

diff --git a/lstchain/image/modifier.py b/lstchain/image/modifier.py
index 89036257d8..046db78035 100644
--- a/lstchain/image/modifier.py
+++ b/lstchain/image/modifier.py
@@ -23,7 +23,7 @@
     'calculate_noise_parameters',
     'random_psf_smearer',
     'set_numba_seed',
-    'WaveformNsbTunner',
+    'WaveformNsbTuner',
 ]
 
 log = logging.getLogger(__name__)
@@ -556,7 +556,7 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
     # Just the product of added_nsb and original_nsb (first two arguments
     # below) is relevant! Units "GHz" (p.e./ns)
     nsb_tuner = [None]
-    # We now create the instances of WaveformNsbTunner to add the different
+    # We now create the instances of WaveformNsbTuner to add the different
     # levels of noise to the MC waveforms.
     # NOTE: the waveform gets updated every time, the noise addition is
     # cumulative (hence the np.diff):
@@ -568,8 +568,9 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
         # just creates it with a single value, for tel_id
         nsb = {tel_id: nsb_rate * u.GHz for tel_id in
                config['source_config']['LSTEventSource']['allowed_tels']}
-        nsb_tuner.append(WaveformNsbTunner(nsb, pulse_templates,
-                                           charge_spe_cumulative_pdf, 10))
+        nsb_tuner.append(WaveformNsbTuner(nsb, pulse_templates,
+                                          charge_spe_cumulative_pdf,
+                                          pre_computed_multiplicity=10))
     # last argument means it will precompute 10 * 1855 (pixels) * 2 (gains)
     # noise waveforms to pick from during the actual introduction of the noise
 
@@ -626,8 +627,9 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
     mc_reader = EventSource(input_url=simtel_filename, config=Config(config))
     nsb = {tel_id: extra_nsb * u.GHz for tel_id in
            config['source_config']['LSTEventSource']['allowed_tels']}
-    tuner = WaveformNsbTunner(nsb, pulse_templates,
-                              charge_spe_cumulative_pdf, 10)
+    tuner = WaveformNsbTuner(nsb, pulse_templates,
+                             charge_spe_cumulative_pdf,
+                             pre_computed_multiplicity=10)
     final_mc_qped = []
     numevents = 0
     for event in mc_reader:
@@ -659,7 +661,7 @@ def custom_function(x, a, b, c):
     """
     return a + b * x ** c
 
-class WaveformNsbTunner:
+class WaveformNsbTuner:
     """
     Handles the injection of additional NSB pulses in waveforms.
     """
diff --git a/lstchain/image/tests/test_modifier.py b/lstchain/image/tests/test_modifier.py
index 5537fa0fe0..5d4625b702 100644
--- a/lstchain/image/tests/test_modifier.py
+++ b/lstchain/image/tests/test_modifier.py
@@ -82,7 +82,7 @@ def test_tune_nsb_on_waveform():
     from ctapipe_io_lst import LSTEventSource
     from scipy.interpolate import interp1d
     from lstchain.io.io import get_resource_path
-    from lstchain.image.modifier import WaveformNsbTunner
+    from lstchain.image.modifier import WaveformNsbTuner
     from lstchain.data.normalised_pulse_template import NormalizedPulseTemplate
 
     waveform = np.array(
@@ -108,10 +108,10 @@ def test_tune_nsb_on_waveform():
     charge_spe_cumulative_pdf = interp1d(spe_integral, spe[0], kind='cubic',
                                          bounds_error=False, fill_value=0.,
                                          assume_sorted=True)
-    nsb_tunner = WaveformNsbTunner(added_nsb_rate,
-                                   pulse_templates,
-                                   charge_spe_cumulative_pdf,
-                                   pre_computed_multiplicity=0)
+    nsb_tunner = WaveformNsbTuner(added_nsb_rate,
+                                  pulse_templates,
+                                  charge_spe_cumulative_pdf,
+                                  pre_computed_multiplicity=0)
     nsb_tunner.tune_nsb_on_waveform(waveform, 1, gain, subarray)
     assert np.any(waveform != 0)
     assert np.isclose(np.mean(waveform), 0.0, atol=0.2)
@@ -119,10 +119,10 @@ def test_tune_nsb_on_waveform():
         [[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
          [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]]
     )
-    nsb_tunner = WaveformNsbTunner(added_nsb_rate,
-                                   pulse_templates,
-                                   charge_spe_cumulative_pdf,
-                                   pre_computed_multiplicity=10)
+    nsb_tunner = WaveformNsbTuner(added_nsb_rate,
+                                  pulse_templates,
+                                  charge_spe_cumulative_pdf,
+                                  pre_computed_multiplicity=10)
     nsb_tunner.tune_nsb_on_waveform(waveform, 1, gain, subarray)
     assert np.any(waveform != 0)
     assert np.isclose(np.mean(waveform), 0.0, atol=0.2)
diff --git a/lstchain/reco/r0_to_dl1.py b/lstchain/reco/r0_to_dl1.py
index 492d462a7a..41314154a3 100644
--- a/lstchain/reco/r0_to_dl1.py
+++ b/lstchain/reco/r0_to_dl1.py
@@ -37,7 +37,7 @@
 from ..calib.camera import load_calibrator_from_config
 from ..calib.camera.calibration_calculator import CalibrationCalculator
 from ..image.cleaning import apply_dynamic_cleaning
-from ..image.modifier import calculate_required_additional_nsb, WaveformNsbTunner
+from ..image.modifier import calculate_required_additional_nsb, WaveformNsbTuner
 from .reconstructor import TimeWaveformFitter
 from ..image.muon import analyze_muon_event, tag_pix_thr
 from ..image.muon import create_muon_table, fill_muon_event
@@ -442,10 +442,10 @@ def r0_to_dl1(
                 nsb_per_tel = {tel_id: nsb_tuning_rate * u.GHz for tel_id in
                                allowed_tels}
 
-                nsb_tuner = WaveformNsbTunner(nsb_per_tel,
-                                              pulse_templates,
-                                              charge_spe_cumulative_pdf,
-                                              pre_computed_multiplicity)
+                nsb_tuner = WaveformNsbTuner(nsb_per_tel,
+                                             pulse_templates,
+                                             charge_spe_cumulative_pdf,
+                                             pre_computed_multiplicity)
             else:
                 logger.warning('NSB tuning on waveform active in config but file is real data, option will be ignored')
                 nsb_tuning = False