Merge branch 'master' of https://gitlab.mpcdf.mpg.de/ievo/icrr-mpp-pipe

gti_generator.py → gti.py

@@ -153,129 +153,280 @@ def intersect_time_intervals(intervals1, intervals2):
  
     return joined_intervals
 
+class GTIGenerator:
+    def __init__(self, config=None, verbose=False):
+        self._config = config
+        self.verbose = verbose
+
+    @property
+    def config(self):
+        return self._config.copy()
+
+    @config.setter
+    def config(self, new_config):
+        if 'event_list' not in new_config:
+            raise ValueError('GTI generator error: the configuration dict is missing the "event_list" section.')
+        
+        self._config = new_config.copy()
 
-def gti_generator(config):
-    """
-    GTI list generator.
+    def _itentify_data_taking_time_edges(self, file_list, max_tdiff=1):
+        if not file_list:
+            raise ValueError("GTI generator: no files to process")
 
-    Parameters
-    ----------
-    config: dict
-        Configuration dictionary (e.g. read from a YAML file).
+        mjd = scipy.zeros(0)
+        tdiff = scipy.zeros(0)
 
-    Returns
-    -------
-    joint_intervals: list
-        A list of (TStart, TStop) lists, representing the identified GTIs.
+        if self.verbose:
+            info_message(f"identifying data taking time edges", "GTI generator")
 
-    """
+        for fnum, fname in enumerate(file_list):
+            with uproot.open(fname) as input_stream:
+                if self.verbose:
+                    info_message(f"processing file {fnum+1} / {len(file_list)}", "GTI generator")
 
-    if 'event_list' not in config:
-        raise ValueError('GTI generator error: the configuration file is missing the "event_list" section. Exiting.')
+                _tdiff = input_stream['Events']['MRawEvtHeader.fTimeDiff'].array()
+                
+                _mjd = input_stream['Events']['MTime.fMjd'].array()
+                _millisec = input_stream['Events']['MTime.fTime.fMilliSec'].array()
+                _nanosec = input_stream['Events']['MTime.fNanoSec'].array()
 
-    # Identifying the files to read
-    info_message("looking for the files", "GTI generator")
-    file_mask = config["data_files"]["data"]["test_sample"]["magic2"]["input_mask"]
+                _mjd = _mjd + (_millisec / 1e3 + _nanosec / 1e9) / 86400.0
+                
+                mjd = scipy.concatenate((mjd, _mjd))
+                tdiff = scipy.concatenate((tdiff, _tdiff))
 
-    file_list = glob.glob(file_mask)
+        sort_args = mjd.argsort()
 
-    if not file_list:
-        raise ValueError("No files to process")
+        not_edge = tdiff < max_tdiff
 
-    # Containers for the data points
-    dfs = {
-        'dc': pandas.DataFrame(),
-        'l3rate': pandas.DataFrame()
-    }
+        time_intervals = identify_time_edges(
+            mjd[sort_args], 
+            not_edge[sort_args], 
+            max_time_diff=max_tdiff / 86400.0
+        )
 
-    # Removing the containers, not specified in the configuration card
-    if "l3rate" not in config['event_list']['cuts']:
-        del dfs['l3rate']
+        return time_intervals
 
-    if "dc" not in config['event_list']['cuts']:
-        del dfs['dc']
+    def _itentify_dc_time_edges(self, file_list):
+        if not file_list:
+            raise ValueError("GTI generator: no files to process")
 
-    # Looping over the data files
-    for fnum, file_name in enumerate(file_list):
-        info_message(f"processing file {fnum+1} / {len(file_list)}", "GTI generator")
+        if "dc" not in self.config['event_list']['cuts']['quality']:
+            raise ValueError("GTI generator: no DC cuts given")
 
-        with uproot.open(file_name) as input_stream:
-            
-            # --- DC ---
-            if "dc" in dfs:
+        if self.verbose:
+            info_message(f"identifying DC time edges", "GTI generator")
+
+        df = pandas.DataFrame()
+
+        # Looping over the data files
+        for fnum, file_name in enumerate(file_list):
+            if self.verbose:
+                info_message(f"processing file {fnum+1} / {len(file_list)}", "GTI generator")
+
+            with uproot.open(file_name) as input_stream:
                 mjd = input_stream["Camera"]["MTimeCamera.fMjd"].array()
                 millisec = input_stream["Camera"]["MTimeCamera.fTime.fMilliSec"].array()
                 nanosec = input_stream["Camera"]["MTimeCamera.fNanoSec"].array()
 
                 df_ = pandas.DataFrame()
-                
+                    
                 df_['mjd'] = mjd + (millisec / 1e3 + nanosec / 1e9) / 86400
                 df_['value'] = input_stream["Camera"]["MReportCamera.fMedianDC"].array()
-                
-                dfs['dc'] = dfs['dc'].append(df_)
-            
-            # --- L3 rate ---
-            if "l3rate" in dfs:
+                    
+                df = df.append(df_)
+
+        df = df.sort_values(by=['mjd'])
+
+        cut = self.config['event_list']['cuts']['quality']['dc']
+           
+        selection = df.query(cut)
+
+        _, idx, _ = scipy.intersect1d(df["mjd"], selection["mjd"], return_indices=True)
+
+        criterion = scipy.repeat(False, len(df["mjd"]))
+        criterion[idx] = True
+
+        time_intervals = identify_time_edges(
+            df["mjd"], 
+            criterion, 
+            max_time_diff=self.config['event_list']['max_time_diff']
+        )
+
+        return time_intervals
+
+    def _itentify_l3rate_time_edges(self, file_list):
+        if not file_list:
+            raise ValueError("GTI generator: no files to process")
+
+        if "l3rate" not in self.config['event_list']['cuts']['quality']:
+            raise ValueError("GTI generator: no L3 rate cuts given")
+
+        if self.verbose:
+            info_message(f"identifying L3 rate time edges", "GTI generator")
+
+        df = pandas.DataFrame()
+
+        # Looping over the data files
+        for fnum, file_name in enumerate(file_list):
+            info_message(f"processing file {fnum+1} / {len(file_list)}", "GTI generator")
+
+            with uproot.open(file_name) as input_stream:
                 mjd = input_stream["Trigger"]["MTimeTrigger.fMjd"].array()
                 millisec = input_stream["Trigger"]["MTimeTrigger.fTime.fMilliSec"].array()
                 nanosec = input_stream["Trigger"]["MTimeTrigger.fNanoSec"].array()
 
                 df_ = pandas.DataFrame()
-                
+                    
                 df_['mjd'] = mjd + (millisec / 1e3 + nanosec / 1e9) / 86400
                 df_['value'] = input_stream["Trigger"]["MReportTrigger.fL3Rate"].array()
-                
-                
-                dfs['l3rate'] = dfs['l3rate'].append(df_)
+                    
+                    
+                df = df.append(df_)
+
+        df = df.sort_values(by=['mjd'])
 
-    # Sorting data points by date is needed for the time intervals identification
-    for key in dfs:
-        dfs[key] = dfs[key].sort_values(by=['mjd'])
+        cut = self.config['event_list']['cuts']['quality']['l3rate']
+           
+        selection = df.query(cut)
 
-    info_message("identifying GTIs", "GTI generator")
+        _, idx, _ = scipy.intersect1d(df["mjd"], selection["mjd"], return_indices=True)
 
-    time_intervals_list = []
+        criterion = scipy.repeat(False, len(df["mjd"]))
+        criterion[idx] = True
+
+        time_intervals = identify_time_edges(
+            df["mjd"], 
+            criterion, 
+            max_time_diff=self.config['event_list']['max_time_diff']
+        )
 
-    # Identifying DC-related GTIs
-    if "dc" in dfs:
-        cut = config['event_list']['cuts']['dc']
+        return time_intervals
         
-        selection = dfs['dc'].query(cut)
 
-        _, idx, _ = scipy.intersect1d(dfs['dc']["mjd"], selection["mjd"], return_indices=True)
+    def process_files(self, file_list):
+        """
+        GTI list generator.
 
-        criterion = scipy.repeat(False, len(dfs['dc']["mjd"]))
-        criterion[idx] = True
+        Parameters
+        ----------
+        config: dict
+            Configuration dictionary (e.g. read from a YAML file).
 
-        time_intervals = identify_time_edges(dfs['dc']["mjd"], criterion, max_time_diff=config['event_list']['max_time_diff'])
+        Returns
+        -------
+        joint_intervals: list
+            A list of (TStart, TStop) lists, representing the identified GTIs.
 
-        time_intervals_list.append(time_intervals)
+        """
 
-    # Identifying L3-related GTIs
-    if "l3rate" in dfs:
-        cut = config['event_list']['cuts']['l3rate']
-        
-        selection = dfs['l3rate'].query(cut)
+        if not self.config:
+            raise ValueError("GTIGenerator: configuration is not set")
 
-        _, idx, _ = scipy.intersect1d(dfs['l3rate']["mjd"], selection["mjd"], return_indices=True)
+        # # Identifying the files to read
+        # info_message("looking for the files", "GTI generator")
+        # file_list = glob.glob(file_mask)
 
-        criterion = scipy.repeat(False, len(dfs['l3rate']["mjd"]))
-        criterion[idx] = True
+        # if not file_list:
+        #     raise ValueError("No files to process")
+
+        # # Containers for the data points
+        # dfs = {
+        #     'dc': pandas.DataFrame(),
+        #     'l3rate': pandas.DataFrame()
+        # }
+
+        # # Removing the containers, not specified in the configuration card
+        # if "l3rate" not in config['event_list']['cuts']:
+        #     del dfs['l3rate']
+
+        # if "dc" not in config['event_list']['cuts']:
+        #     del dfs['dc']
+
+        # # Looping over the data files
+        # for fnum, file_name in enumerate(file_list):
+        #     info_message(f"processing file {fnum+1} / {len(file_list)}", "GTI generator")
+
+        #     with uproot.open(file_name) as input_stream:
+                
+        #         # --- DC ---
+        #         if "dc" in dfs:
+        #             mjd = input_stream["Camera"]["MTimeCamera.fMjd"].array()
+        #             millisec = input_stream["Camera"]["MTimeCamera.fTime.fMilliSec"].array()
+        #             nanosec = input_stream["Camera"]["MTimeCamera.fNanoSec"].array()
+
+        #             df_ = pandas.DataFrame()
+                    
+        #             df_['mjd'] = mjd + (millisec / 1e3 + nanosec / 1e9) / 86400
+        #             df_['value'] = input_stream["Camera"]["MReportCamera.fMedianDC"].array()
+                    
+        #             dfs['dc'] = dfs['dc'].append(df_)
+                
+        #         # --- L3 rate ---
+        #         if "l3rate" in dfs:
+        #             mjd = input_stream["Trigger"]["MTimeTrigger.fMjd"].array()
+        #             millisec = input_stream["Trigger"]["MTimeTrigger.fTime.fMilliSec"].array()
+        #             nanosec = input_stream["Trigger"]["MTimeTrigger.fNanoSec"].array()
+
+        #             df_ = pandas.DataFrame()
+                    
+        #             df_['mjd'] = mjd + (millisec / 1e3 + nanosec / 1e9) / 86400
+        #             df_['value'] = input_stream["Trigger"]["MReportTrigger.fL3Rate"].array()
+                    
+                    
+        #             dfs['l3rate'] = dfs['l3rate'].append(df_)
+
+        # # Sorting data points by date is needed for the time intervals identification
+        # for key in dfs:
+        #     dfs[key] = dfs[key].sort_values(by=['mjd'])
+
+        # info_message("identifying GTIs", "GTI generator")
+
+        # time_intervals_list = []
+
+        # # Identifying DC-related GTIs
+        # if "dc" in dfs:
+        #     cut = config['event_list']['cuts']['dc']
+            
+        #     selection = dfs['dc'].query(cut)
+
+        #     _, idx, _ = scipy.intersect1d(dfs['dc']["mjd"], selection["mjd"], return_indices=True)
+
+        #     criterion = scipy.repeat(False, len(dfs['dc']["mjd"]))
+        #     criterion[idx] = True
+
+        #     time_intervals = identify_time_edges(dfs['dc']["mjd"], criterion, max_time_diff=config['event_list']['max_time_diff'])
+
+        #     time_intervals_list.append(time_intervals)
+
+        # # Identifying L3-related GTIs
+        # if "l3rate" in dfs:
+        #     cut = config['event_list']['cuts']['l3rate']
+            
+        #     selection = dfs['l3rate'].query(cut)
+
+        #     _, idx, _ = scipy.intersect1d(dfs['l3rate']["mjd"], selection["mjd"], return_indices=True)
+
+        #     criterion = scipy.repeat(False, len(dfs['l3rate']["mjd"]))
+        #     criterion[idx] = True
 
-        time_intervals = identify_time_edges(dfs['l3rate']["mjd"], criterion, max_time_diff=config['event_list']['max_time_diff'])
+        #     time_intervals = identify_time_edges(dfs['l3rate']["mjd"], criterion, max_time_diff=config['event_list']['max_time_diff'])
 
-        time_intervals_list.append(time_intervals)
+        #     time_intervals_list.append(time_intervals)
 
-    if not len(time_intervals_list):
-        raise ValueError("No valid selection cuts provided in the configuration file.")
+        time_intervals_list = [
+            self._itentify_data_taking_time_edges(file_list),
+            self._itentify_dc_time_edges(file_list),
+            self._itentify_l3rate_time_edges(file_list)
+        ]
 
-    joint_intervals = time_intervals_list[0]
+        joint_intervals = time_intervals_list[0]
 
-    # Joining all found GTIs
-    for i in range(1, len(time_intervals_list)):
-        joint_intervals = intersect_time_intervals(joint_intervals, time_intervals_list[i])
+        # Joining all found GTIs
+        for i in range(1, len(time_intervals_list)):
+            joint_intervals = intersect_time_intervals(joint_intervals, time_intervals_list[i])
 
-    return joint_intervals
+        return joint_intervals
     
 
 # =================

make_irf.py

@@ -656,7 +656,7 @@ irf_generator.set_energy_binning(min_energy=0.1, max_energy=30, n_energy_bins=10
 irf_generator.set_theta_binning(min_theta=0.0, max_theta=1.5, n_theta_bins=5)
 irf_generator.set_migra_binning(min_migra=0.2, max_migra=5.0, n_migra_bins=5)
 
-irf_generator.set_cuts('(multiplicity > 1) & (abs(pos_angle_shift_reco - 0.5) > 0.4) & (event_class_0 > 0.8)')
+irf_generator.set_cuts(config['event_list']['cuts']['selection'])
 
 irf_generator.generate_irf('crab_irf.fits')