Add muon analysis on calibrated data

magicctapipe/image/muons/muon_analysis.py

@@ -34,29 +34,30 @@ def perform_muon_analysis(muon_parameters, event, telescope_id, telescope_name,
 
     """
     if data_type == 'obs':
-        bad_pixels = event.mon.tel[telescope_id].calibration.unusable_pixels[0]
-        # Set to 0 unreliable pixels:
-        image = image * (~bad_pixels)
+        try:
+            bad_pixels = event.mon.tel[telescope_id].calibration.unusable_pixels[0]
+            # Set to 0 unreliable pixels:
+            image = image * (~bad_pixels)
+        except TypeError:
+            pass
     # process only promising events, in terms of # of pixels with large signals:
     thr_low = good_ring_config['thr_low'] if 'thr_low' in good_ring_config else 50
     if tag_pix_thr(image, thr_low=thr_low):
-        # re-calibrate r1 to obtain new dl1, using a more adequate pulse integrator for muon rings
-        numsamples = event.r1.tel[telescope_id].waveform.shape[1]  # not necessarily the same as in r0!
         if data_type == 'obs':
-            bad_pixels_hg = event.mon.tel[telescope_id].calibration.unusable_pixels[0]
-            bad_pixels_lg = event.mon.tel[telescope_id].calibration.unusable_pixels[1]
-            # Now set to 0 all samples in unreliable pixels. Important for global peak
-            # integrator in case of crazy pixels!  TBD: can this be done in a simpler
-            # way?
-            bad_pixels = bad_pixels_hg | bad_pixels_lg
-            bad_waveform = np.transpose(np.array(numsamples * [bad_pixels]))
-
-            # print('hg bad pixels:',np.where(bad_pixels_hg))
-            # print('lg bad pixels:',np.where(bad_pixels_lg))
-
-            event.r1.tel[telescope_id].waveform *= ~bad_waveform
-            image = image * (~bad_pixels)
-        r1_dl1_calibrator_for_muon_rings(event)
+            try:
+                bad_pixels_hg = event.mon.tel[telescope_id].calibration.unusable_pixels[0]
+                bad_pixels_lg = event.mon.tel[telescope_id].calibration.unusable_pixels[1]
+                bad_pixels = bad_pixels_hg | bad_pixels_lg
+                image = image * (~bad_pixels)
+            except TypeError:
+                pass
+        if r1_dl1_calibrator_for_muon_rings is not None:
+            # re-calibrate r1 to obtain new dl1, using a more adequate pulse integrator for muon rings
+            numsamples = event.r1.tel[telescope_id].waveform.shape[1]  # not necessarily the same as in r0!
+            if data_type == 'obs':
+                bad_waveform = np.transpose(np.array(numsamples * [bad_pixels]))
+                event.r1.tel[telescope_id].waveform *= ~bad_waveform
+            r1_dl1_calibrator_for_muon_rings(event)
 
         # Check again: with the extractor for muon rings (most likely GlobalPeakWindowSum)
         # perhaps the event is no longer promising (e.g. if it has a large time evolution)
@@ -72,26 +73,31 @@ def perform_muon_analysis(muon_parameters, event, telescope_id, telescope_name,
                 analyze_muon_event(subarray, telescope_id, event_id,
                                    image, good_ring_config,
                                    plot_rings=False, plots_path='')
-            #                      plot_rings=True, plots_path='../../../../data/'+telescope_name+'/')
+                                  #plot_rings=True, plots_path='../data/real'+telescope_name+'/')
             #           (test) plot muon rings as png files
 
-            # Now we want to obtain the waveform sample (in HG & LG) at which the ring light peaks:
-            bright_pixels = image > min_pe_for_muon_t_calc
-            selected_gain = event.r1.tel[telescope_id].selected_gain_channel
-            mask_hg = bright_pixels & (selected_gain == 0)
-            mask_lg = bright_pixels & (selected_gain == 1)
+            if r1_dl1_calibrator_for_muon_rings is not None:
+                # Now we want to obtain the waveform sample (in HG & LG) at which the ring light peaks:
+                bright_pixels = image > min_pe_for_muon_t_calc
+                selected_gain = event.r1.tel[telescope_id].selected_gain_channel
+                mask_hg = bright_pixels & (selected_gain == 0)
+                mask_lg = bright_pixels & (selected_gain == 1)
 
-            bright_pixels_waveforms_hg = event.r1.tel[telescope_id].waveform[mask_hg, :]
-            bright_pixels_waveforms_lg = event.r1.tel[telescope_id].waveform[mask_lg, :]
-            stacked_waveforms_hg = np.sum(bright_pixels_waveforms_hg, axis=0)
-            stacked_waveforms_lg = np.sum(bright_pixels_waveforms_lg, axis=0)
+                bright_pixels_waveforms_hg = event.r1.tel[telescope_id].waveform[mask_hg, :]
+                bright_pixels_waveforms_lg = event.r1.tel[telescope_id].waveform[mask_lg, :]
+                stacked_waveforms_hg = np.sum(bright_pixels_waveforms_hg, axis=0)
+                stacked_waveforms_lg = np.sum(bright_pixels_waveforms_lg, axis=0)
 
-            # stacked waveforms from all bright pixels; shape (ngains, nsamples)
-            hg_peak_sample = np.argmax(stacked_waveforms_hg, axis=-1)
-            lg_peak_sample = np.argmax(stacked_waveforms_lg, axis=-1)
+                # stacked waveforms from all bright pixels; shape (ngains, nsamples)
+                hg_peak_sample = np.argmax(stacked_waveforms_hg, axis=-1)
+                lg_peak_sample = np.argmax(stacked_waveforms_lg, axis=-1)
+            else:
+                hg_peak_sample, lg_peak_sample = -1, -1
+
+        mc_energy = event.simulation.shower.energy if data_type == 'mc' else -1
 
         if good_ring:
-            fill_muon_event(-1,
+            fill_muon_event(mc_energy,
                             muon_parameters,
                             good_ring,
                             event.index.event_id,
@@ -107,7 +113,6 @@ def perform_muon_analysis(muon_parameters, event, telescope_id, telescope_name,
                             hg_peak_sample, lg_peak_sample)
             muon_parameters['telescope_name'].append(telescope_name)
 
-
 # TODO replace this copy of dev version of lstchain (last release tag v0.9.4) by an import once v0.9.5+ is out
 import matplotlib.pyplot as plt
 import astropy.units as u

magicctapipe/scripts/lst1_magic/config.yaml

@@ -64,8 +64,8 @@ MAGIC:
 
     muon_ring:
         thr_low: 25
-        tailcut: [8, 4]
-        ring_completeness_threshold: 15
+        tailcut: [12, 8]
+        ring_completeness_threshold: 25
 
 
 event_coincidence:

magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py

@@ -196,7 +196,7 @@ def mc_dl0_to_dl1(input_file, output_dir, config, muons_analysis):
             extractor_muon_name_lst, subarray=subarray, config=config_extractor_lst
         )
         r1_dl1_calibrator_for_muon_rings[tel_id_lst1] = CameraCalibrator(subarray,
-                                                                        image_extractor=extractor_lst_muons)
+                                                                         image_extractor=extractor_lst_muons)
         extractor_muon_name_magic = 'GlobalPeakWindowSum'
         extractor_magic_muons = ImageExtractor.from_name(
             extractor_muon_name_magic, subarray=subarray, config=config_extractor_magic

magicctapipe/scripts/lst1_magic/muon_analysis_LST_or_MAGIC_data.py

@@ -0,0 +1,165 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+"""
+Author: Gabriel Emery
+
+This script processes MAGIC calibrated data (*_Y_*.root) to perform the muon ring analysis.
+Current naming allows LST but first implementation will only handle MAGIC as script for LST are available in lstchain.
+
+Usage:
+$ python muon_analysis_LST_or_MAGIC_data.py
+--input-file ./data/calibrated/20201216_M1_05093711.001_Y_CrabNebula-W0.40+035.root
+--output-dir ./data/muons
+--config-file ./config.yaml
+"""
+import argparse
+import yaml
+import logging
+import numpy as np
+from pathlib import Path
+
+from astropy.table import Table
+from ctapipe.io import EventSource
+from lstchain.image.muon import create_muon_table
+from magicctapipe.image import MAGICClean
+from magicctapipe.image.muons import perform_muon_analysis
+
+logger = logging.getLogger(__name__)
+logger.addHandler(logging.StreamHandler())
+logger.setLevel(logging.INFO)
+
+pedestal_types = [
+    'fundamental',
+    'from_extractor',
+    'from_extractor_rndm',
+]
+
+
+def magic_muons_from_cal(input_file, output_dir, config, process_run):
+    """
+    Process all event a single telescope MAGIC calibrated data run or subrun to perform the muon ring analysis.
+
+    Parameters
+    ----------
+    input_file:
+    output_dir:
+    config:
+    process_run:
+
+    """
+
+    event_source = EventSource(input_url=input_file)
+    subarray = event_source.subarray
+    obs_id = event_source.obs_ids[0]
+    tel_id = event_source.telescope
+
+    # Create the table which will contain the selected muon ring parameters
+    muon_parameters = create_muon_table()
+    muon_parameters['telescope_name'] = []
+
+    logger.info(f'\nProcess the following data (process_run = {process_run}):')
+
+    # Configure the MAGIC cleaning:
+    config_cleaning = config['MAGIC']['magic_clean']
+
+    if config_cleaning['find_hotpixels'] == 'auto':
+        config_cleaning.update({'find_hotpixels': True})
+
+    logger.info('\nConfiguration for the image cleaning:')
+    logger.info(config_cleaning)
+
+    ped_type = config_cleaning.pop('pedestal_type')
+    i_ped_type = np.where(np.array(pedestal_types) == ped_type)[0][0]
+
+    camera_geom = event_source.subarray.tel[tel_id].camera.geometry
+    magic_clean = MAGICClean(camera_geom, config_cleaning)
+
+    # Prepare for saving muons data to an output file:
+    Path(output_dir).mkdir(exist_ok=True, parents=True)
+
+    if process_run:
+        output_file = f'{output_dir}/muons_M{tel_id}.Run{obs_id:08}.fits'
+    else:
+        subrun_id = event_source.metadata['subrun_number']
+        output_file = f'{output_dir}/muons_M{tel_id}.Run{obs_id:08}.{subrun_id:03}.fits'
+
+    # Start processing events:
+    logger.info('\nProcessing the events:')
+
+    # load muon analysis config
+    muon_config = {}
+    if 'muon_ring' in config['MAGIC']:
+        muon_config = config['MAGIC']['muon_ring']
+    # Select the telescope name to be filed in muon_parameters['telescope_name']
+    if tel_id == 1:
+        tel_name = f'MAGIC-I'
+    elif tel_id == 2:
+        tel_name = f'MAGIC-II'
+    else:
+        tel_name = f'MAGIC_?'
+
+    for event in event_source:
+
+        # Apply the image cleaning:
+        dead_pixels = event.mon.tel[tel_id].pixel_status.hardware_failing_pixels[0]
+        try:
+            badrms_pixels = event.mon.tel[tel_id].pixel_status.pedestal_failing_pixels[i_ped_type]
+        except TypeError:
+            badrms_pixels = np.zeros(dead_pixels.shape, dtype=bool)
+        unsuitable_mask = np.logical_or(dead_pixels, badrms_pixels)
+
+        signal_pixels, image, peak_time = magic_clean.clean_image(
+            event.dl1.tel[tel_id].image, event.dl1.tel[tel_id].peak_time, unsuitable_mask,
+        )
+        perform_muon_analysis(muon_parameters,
+                              event=event,
+                              telescope_id=tel_id,
+                              telescope_name=tel_name,
+                              image=image,
+                              subarray=subarray,
+                              r1_dl1_calibrator_for_muon_rings=None,
+                              good_ring_config=muon_config,
+                              data_type='obs')
+
+    table = Table(muon_parameters)
+    table.write(output_file, format='fits', overwrite=True)
+    logger.info(f'\nOutput muons file: {output_file}')
+
+
+def main():
+
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument(
+        '--input-file', '-i', dest='input_file', type=str, required=True,
+        help='Path to an input MAGIC calibrated data file (*_Y_*.root).',
+    )
+
+    parser.add_argument(
+        '--output-dir', '-o', dest='output_dir', type=str, default='./data',
+        help='Path to a directory where to save an output muon file.',
+    )
+
+    parser.add_argument(
+        '--config-file', '-c', dest='config_file', type=str, default='./config.yaml',
+        help='Path to a yaml configuration file.',
+    )
+
+    parser.add_argument(
+        '--process-run', dest='process_run', action='store_true',
+        help='Processes all the sub-run files of the same observation ID at once.',
+    )
+
+    args = parser.parse_args()
+
+    with open(args.config_file, 'rb') as f:
+        config = yaml.safe_load(f)
+
+    # Process the input data:
+    magic_muons_from_cal(args.input_file, args.output_dir, config, args.process_run)
+
+    logger.info('\nDone.')
+
+if __name__ == '__main__':
+    main()