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/modifier.py b/lstchain/image/modifier.py
index 9f52390238..046db78035 100644
--- a/lstchain/image/modifier.py
+++ b/lstchain/image/modifier.py
@@ -6,13 +6,16 @@
     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.config import read_configuration_file
+from lstchain.io import standard_config, 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',
@@ -20,7 +23,7 @@
     'calculate_noise_parameters',
     'random_psf_smearer',
     'set_numba_seed',
-    'WaveformNsbTunner',
+    'WaveformNsbTuner',
 ]
 
 log = logging.getLogger(__name__)
@@ -406,8 +409,75 @@ def calculate_noise_parameters(simtel_filename, data_dl1_filename,
         extra_noise_in_bright_pixels
 
 
+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 (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)
+
+    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')
+    # 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')
+    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
+
+    # 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'] == 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)
+    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()}')
+    # 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, tel_id
+
+
+
 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,9 +498,10 @@ 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
+    mc_ped_variance: Pedestal variance from MC, AFTER THE TUNING!
 
     """
 
@@ -439,119 +510,168 @@ def calculate_required_additional_nsb(simtel_filename, data_dl1_filename, config
     if config is None:
         config = standard_config
 
-    # Real data DL1 tables:
-    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])
-
-    log.info(f'Good and not too bright pixels: {good_pixels.sum()}')
+    # 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(
+            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))
 
-    # 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 = mc_reader.subarray
+
+    # 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 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):
+
+
+    for nsb_rate in np.diff(total_added_nsb):
+        # 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(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
+
+    # 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)
+
+        # First store the charges with no added NSB:
+        modified_integrated_charge[0].extend(event.dl1.tel[tel_id].image[
+                                                 pedmask])
+
+        # 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
+
+            # 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])
+
+    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])
+
+    # Since the level of noise in MC is low, it should hardly ever happen
+    # that extra_nsb is negative, but just in case:
+    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))
+    nsb = {tel_id: extra_nsb * u.GHz for tel_id in
+           config['source_config']['LSTEventSource']['allowed_tels']}
+    tuner = WaveformNsbTuner(nsb, pulse_templates,
+                             charge_spe_cumulative_pdf,
+                             pre_computed_multiplicity=10)
+    final_mc_qped = []
+    numevents = 0
+    for event in mc_reader:
+        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)
 
-    # 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)
+        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])
 
-    # Find the additional noise (in data w.r.t. MC) for the unbiased extractor
-    # The idea is that pedestal variance scales with NSB
+    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
 
 
-class WaveformNsbTunner:
+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 WaveformNsbTuner:
     """
     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: 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`
@@ -561,8 +681,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 +715,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 +799,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:],
diff --git a/lstchain/image/tests/test_modifier.py b/lstchain/image/tests/test_modifier.py
index a23c25eff2..5d4625b702 100644
--- a/lstchain/image/tests/test_modifier.py
+++ b/lstchain/image/tests/test_modifier.py
@@ -70,9 +70,11 @@ 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.11, atol=0.01)
-    assert np.isclose(extra_nsb, -1.0)
+    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():
@@ -80,14 +82,14 @@ 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(
         [[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,11 +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_fraction,
-                                   original_nsb,
-                                   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)
@@ -118,11 +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_fraction,
-                                   original_nsb,
-                                   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 ce08cbaa64..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
@@ -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
@@ -403,25 +403,27 @@ def r0_to_dl1(
             metadata=metadata,
         )
     nsb_tuning = False
-    nsb_tunner = None
-    if 'waveform_nsb_tuning' in config.keys():
+    nsb_tuner = None
+    if 'waveform_nsb_tuning' in config:
         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']}
-                if 'nsb_tuning_ratio' in config['waveform_nsb_tuning'].keys():
+                if 'nsb_tuning_rate_GHz' in config[
+                    'waveform_nsb_tuning']:
                     # 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,28 +432,33 @@ 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('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 = {tel_id: nsb_tuning_rate * u.GHz for tel_id in
+                               allowed_tels}
+
+                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
 
     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']:
             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 +577,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)
diff --git a/lstchain/scripts/lstchain_tune_nsb_waveform.py b/lstchain/scripts/lstchain_tune_nsb_waveform.py
index 1a1f8ee0f4..a91e940827 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_GHz": 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")}"'