Merge pull request #1282 from cta-observatory/wf_noise_tuner

Waveform noise tuner

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
   },

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
   },

lstchain/image/tests/test_modifier.py

@@ -70,24 +70,26 @@ 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():
     import astropy.units as u
     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,23 +108,21 @@ 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)
     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]]
     )
-    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)

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)

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"))
     }
 

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")}"'