Merge pull request #1272 from cta-observatory/Fix_mc_nsb_reading

Fix original MC nsb reading

lstchain/io/io.py

@@ -1,6 +1,5 @@
 import logging
 import os
-import re
 import warnings
 from multiprocessing import Pool
 from contextlib import ExitStack
@@ -22,9 +21,8 @@
 from ctapipe.instrument import SubarrayDescription
 from ctapipe.io import HDF5TableReader, HDF5TableWriter
 
-from eventio import Histograms, EventIOFile
-from eventio.search_utils import yield_toplevel_of_type, yield_all_subobjects
-from eventio.simtel.objects import History, HistoryConfig
+from eventio import Histograms, SimTelFile
+from eventio.search_utils import yield_toplevel_of_type
 
 from pyirf.simulations import SimulatedEventsInfo
 
@@ -59,7 +57,6 @@
     'global_metadata',
     'merge_dl2_runs',
     'merging_check',
-    'parse_cfg_bytestring',
     'read_data_dl2_to_QTable',
     'read_dl2_params',
     'read_mc_dl2_to_QTable',
@@ -1264,37 +1261,34 @@ def remove_duplicated_events(data):
     data.remove_rows(remove_row_list)
 
 
-def parse_cfg_bytestring(bytestring):
-    """
-    Parse configuration as read by eventio
-    :param bytes bytestring: A ``Bytes`` object with configuration data for one parameter
-    :return: Tuple in form ``('parameter_name', 'value')``
+def extract_simulation_nsb(filename):
     """
-    line_decoded = bytestring.decode('utf-8').rstrip()
-    if 'ECHO' in line_decoded or '#' in line_decoded:
-        return None
-    line_list = line_decoded.split('%', 1)[0]  # drop comment
-    res = re.sub(' +', ' ', line_list).strip().split(' ', 1)  # remove extra whitespaces and split
-    return res[0].upper(), res[1]
+    Get current run NSB from configuration in simtel file.
 
+    WARNING : In current MC, correct NSB are logged after 'STORE_PHOTOELECTRONS' entries
+    In any new production, behaviour needs to be verified.
+    New version of simtel will allow to use better metadata.
 
-def extract_simulation_nsb(filename):
-    """
-    Get current run NSB from configuration in simtel file
     :param str filename: Input file name
-    :return array of `float` by tel_id: NSB rate
-    """
-    nsb = []
-    with EventIOFile(filename) as f:
-        for o in yield_all_subobjects(f, [History, HistoryConfig]):
-            if hasattr(o, 'parse'):
-                try:
-                    cfg_element = parse_cfg_bytestring(o.parse()[1])
-                    if cfg_element is not None:
-                        if cfg_element[0] == 'NIGHTSKY_BACKGROUND':
-                            nsb.append(float(cfg_element[1].strip('all:')))
-                except Exception as e:
-                    print('Unexpected end of %s,\n caught exception %s', filename, e)
+    :return dict of `float` by tel_id: NSB rate
+    """
+    nsb = {}
+    next_nsb = False
+    tel_id = 1
+    with SimTelFile(filename) as f:
+        for _, line in f.history:
+            line = line.decode('utf-8').strip().split(' ')
+            if next_nsb and line[0] == 'NIGHTSKY_BACKGROUND':
+                nsb[tel_id] = float(line[1].strip('all:'))
+                tel_id = tel_id+1
+            if line[0] == 'STORE_PHOTOELECTRONS':
+                next_nsb = True
+            else:
+                next_nsb = False
+    log.warning('Original MC night sky background extracted from the config history in the simtel file.\n'
+                'This is done for existing LST MC such as the one created using: '
+                'https://github.com/cta-observatory/lst-sim-config/tree/sim-tel_LSTProd2_MAGICST0316'
+                '\nExtracted values are: ' + str(np.asarray(nsb)) + 'GHz. Check that it corresponds to expectations.')
     return nsb
 
 

lstchain/io/tests/test_io.py

@@ -156,8 +156,7 @@ def test_extract_simulation_nsb(mc_gamma_testfile):
     from lstchain.io.io import extract_simulation_nsb
 
     nsb = extract_simulation_nsb(mc_gamma_testfile)
-    assert np.isclose(nsb[0], 0.246, rtol=0.1)
-    assert np.isclose(nsb[1], 0.217, rtol=0.1)
+    assert np.isclose(nsb[1], 0.246, rtol=0.1)
 
 
 def test_remove_duplicated_events():

lstchain/reco/r0_to_dl1.py

@@ -433,10 +433,8 @@ def r0_to_dl1(
                 charge_spe_cumulative_pdf = interp1d(spe_integral, spe[0], kind='cubic',
                                                      bounds_error=False, fill_value=0.,
                                                      assume_sorted=True)
-                allowed_tel = np.zeros(len(nsb_original), dtype=bool)
-                allowed_tel[np.array(config['source_config']['LSTEventSource']['allowed_tels'])] = True
                 logger.info('Tuning NSB on MC waveform from '
-                            + str(np.asarray(nsb_original)[allowed_tel])
+                            + str(np.asarray(nsb_original))
                             + 'GHz to {0:d}%'.format(int(nsb_tuning_ratio * 100 + 100.5))
                             + ' for telescopes ids ' + str(config['source_config']['LSTEventSource']['allowed_tels']))
                 nsb_tuning_args = [nsb_tuning_ratio, nsb_original, pulse_template, charge_spe_cumulative_pdf]