From 6b76747f8467c8b48645bb7d62c4bb8367324c85 Mon Sep 17 00:00:00 2001 From: Gabriel Emery Date: Wed, 27 Apr 2022 10:04:44 +0200 Subject: [PATCH 1/2] Add script to perfom the muon analysis on calibrated data and modified the muon analysis function to also work for this case. --- magicctapipe/image/muons/muon_analysis.py | 73 ++++---- magicctapipe/scripts/lst1_magic/config.yaml | 4 +- .../lst1_magic/lst1_magic_mc_dl0_to_dl1.py | 2 +- .../muon_analysis_LST_or_MAGIC_data.py | 169 ++++++++++++++++++ 4 files changed, 211 insertions(+), 37 deletions(-) create mode 100644 magicctapipe/scripts/lst1_magic/muon_analysis_LST_or_MAGIC_data.py diff --git a/magicctapipe/image/muons/muon_analysis.py b/magicctapipe/image/muons/muon_analysis.py index 63055990a..c9a2f7ccd 100644 --- a/magicctapipe/image/muons/muon_analysis.py +++ b/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 diff --git a/magicctapipe/scripts/lst1_magic/config.yaml b/magicctapipe/scripts/lst1_magic/config.yaml index 96492cb9c..bf8fb5f15 100644 --- a/magicctapipe/scripts/lst1_magic/config.yaml +++ b/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: diff --git a/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py b/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py index 1eaeb0658..69564564e 100644 --- a/magicctapipe/scripts/lst1_magic/lst1_magic_mc_dl0_to_dl1.py +++ b/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 diff --git a/magicctapipe/scripts/lst1_magic/muon_analysis_LST_or_MAGIC_data.py b/magicctapipe/scripts/lst1_magic/muon_analysis_LST_or_MAGIC_data.py new file mode 100644 index 000000000..67316593d --- /dev/null +++ b/magicctapipe/scripts/lst1_magic/muon_analysis_LST_or_MAGIC_data.py @@ -0,0 +1,169 @@ +#!/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_magic import MAGICEventSource +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 = MAGICEventSource(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:') + + # We give None as calibrator for the muons analysis as we work with extracted charges + r1_dl1_calibrator_for_muon_rings_magic = None + + # 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= + r1_dl1_calibrator_for_muon_rings_magic, + 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() From f8723966760ac370db94163f2f0f867fa67fde69 Mon Sep 17 00:00:00 2001 From: Gabriel Emery Date: Fri, 29 Apr 2022 13:13:24 +0200 Subject: [PATCH 2/2] Generalise EventSource call, remove intermediate variable. --- .../lst1_magic/muon_analysis_LST_or_MAGIC_data.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/magicctapipe/scripts/lst1_magic/muon_analysis_LST_or_MAGIC_data.py b/magicctapipe/scripts/lst1_magic/muon_analysis_LST_or_MAGIC_data.py index 67316593d..c5ed19c83 100644 --- a/magicctapipe/scripts/lst1_magic/muon_analysis_LST_or_MAGIC_data.py +++ b/magicctapipe/scripts/lst1_magic/muon_analysis_LST_or_MAGIC_data.py @@ -20,7 +20,7 @@ from pathlib import Path from astropy.table import Table -from ctapipe_io_magic import MAGICEventSource +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 @@ -49,7 +49,7 @@ def magic_muons_from_cal(input_file, output_dir, config, process_run): """ - event_source = MAGICEventSource(input_url=input_file) + event_source = EventSource(input_url=input_file) subarray = event_source.subarray obs_id = event_source.obs_ids[0] tel_id = event_source.telescope @@ -87,9 +87,6 @@ def magic_muons_from_cal(input_file, output_dir, config, process_run): # Start processing events: logger.info('\nProcessing the events:') - # We give None as calibrator for the muons analysis as we work with extracted charges - r1_dl1_calibrator_for_muon_rings_magic = None - # load muon analysis config muon_config = {} if 'muon_ring' in config['MAGIC']: @@ -121,8 +118,7 @@ def magic_muons_from_cal(input_file, output_dir, config, process_run): telescope_name=tel_name, image=image, subarray=subarray, - r1_dl1_calibrator_for_muon_rings= - r1_dl1_calibrator_for_muon_rings_magic, + r1_dl1_calibrator_for_muon_rings=None, good_ring_config=muon_config, data_type='obs')