diff --git a/magicctapipe/conftest.py b/magicctapipe/conftest.py index df4a3ebc0..6dba1792e 100644 --- a/magicctapipe/conftest.py +++ b/magicctapipe/conftest.py @@ -53,6 +53,11 @@ def temp_DL1_gamma(tmp_path_factory): return tmp_path_factory.mktemp("DL1_gammas") +@pytest.fixture(scope="session") +def temp_DL1_gamma_focal_exc(tmp_path_factory): + return tmp_path_factory.mktemp("DL1_gammas_focal_exc") + + @pytest.fixture(scope="session") def temp_DL1_gamma_train(tmp_path_factory): return tmp_path_factory.mktemp("DL1_gamma_train") @@ -188,6 +193,11 @@ def temp_coinc(tmp_path_factory): return tmp_path_factory.mktemp("coincidence") +@pytest.fixture(scope="session") +def temp_coinc_preoff(tmp_path_factory): + return tmp_path_factory.mktemp("coincidence_pre_offset") + + @pytest.fixture(scope="session") def temp_coinc_stereo(tmp_path_factory): return tmp_path_factory.mktemp("coincidence_stereo") @@ -374,6 +384,12 @@ def config(): return config_path +@pytest.fixture(scope="session") +def config_preoff(): + config_path = resource_file("test_config_preoffset.yaml") + return config_path + + @pytest.fixture(scope="session") def config_monly(): config_path = resource_file("test_config_monly.yaml") diff --git a/magicctapipe/io/__init__.py b/magicctapipe/io/__init__.py index 0b44e3700..986894f35 100644 --- a/magicctapipe/io/__init__.py +++ b/magicctapipe/io/__init__.py @@ -20,6 +20,7 @@ load_train_data_files_tel, save_pandas_data_in_table, telescope_combinations, + telescope_positions, ) from .gadf import ( create_event_hdu, @@ -50,4 +51,5 @@ "load_train_data_files_tel", "save_pandas_data_in_table", "telescope_combinations", + "telescope_positions", ] diff --git a/magicctapipe/io/io.py b/magicctapipe/io/io.py index c4f185d99..7d066ca83 100644 --- a/magicctapipe/io/io.py +++ b/magicctapipe/io/io.py @@ -36,6 +36,7 @@ "load_train_data_files_tel", "save_pandas_data_in_table", "telescope_combinations", + "telescope_positions", ] logger = logging.getLogger(__name__) @@ -108,6 +109,48 @@ def recursive_solution(current_tel, current_comb): return TEL_NAMES, TEL_COMBINATIONS +def telescope_positions(config): + """ + This function computes the telescope positions with respect to the array baricenter. + The array can have any configuration, e.g.: M1+M2+LST1+LST2, the full MAGIC+LST array, etc. + + Parameters + ---------- + config: dict + dictionary generated from an yaml file with information about the telescope IDs. + + Returns + ------- + TEL_POSITIONS: dict + Dictionary with telescope positions in the baricenter reference frame of the adopted array. + """ + + #Telescope positions in meters in a generic reference frame: + RELATIVE_POSITIONS = { + "LST-1" : [-70.930, -52.070, 43.00], + "LST-2" : [-35.270, 66.140, 32.00], + "LST-3" : [75.280 , 50.490, 28.70], + "LST-4" : [30.910 , -64.540, 32.00], + "MAGIC-I" : [-23.540, -191.750, 41.25], + "MAGIC-II" : [-94.05, -143.770, 42.42] + } + + telescopes_in_use = {} + tels=config["mc_tel_ids"] + tel_cp=tels.copy() + for k, v in tel_cp.copy().items(): + if v <= 0: + tel_cp.pop(k) + else: + telescopes_in_use[v] = RELATIVE_POSITIONS[k] + + average_xyz=np.array([RELATIVE_POSITIONS[k] for k in tel_cp.keys()]).mean(axis=0) + TEL_POSITIONS = {} + for k, v in telescopes_in_use.items(): + TEL_POSITIONS[k] = list(np.round(np.asarray(v)-average_xyz,2)) * u.m + return TEL_POSITIONS + + def format_object(input_object): """ Formats a object (dictionary or list) to show its elements. diff --git a/magicctapipe/io/tests/test_io.py b/magicctapipe/io/tests/test_io.py index c44cc8b31..cd31a56db 100644 --- a/magicctapipe/io/tests/test_io.py +++ b/magicctapipe/io/tests/test_io.py @@ -3,6 +3,7 @@ get_dl2_mean, get_stereo_events, load_train_data_files, + load_train_data_files_tel, load_mc_dl2_data_file, load_irf_files, save_pandas_data_in_table, @@ -10,11 +11,13 @@ load_lst_dl1_data_file, load_dl2_data_file, telescope_combinations, + telescope_positions, ) - +import glob import pytest import numpy as np import pandas as pd +import astropy.units as u def test_telescope_combinations(config_gen, config_gen_4lst): @@ -28,6 +31,24 @@ def test_telescope_combinations(config_gen, config_gen_4lst): assert LSTs == {1: 'LST-1', 3: 'LST-2', 2: 'LST-3', 5: 'LST-4'} assert LSTs_comb == {'LST-1_LST-2': [1, 3], 'LST-1_LST-2_LST-3': [1, 3, 2], 'LST-1_LST-2_LST-3_LST-4': [1, 3, 2, 5], 'LST-1_LST-2_LST-4': [1, 3, 5], 'LST-1_LST-3': [1, 2], 'LST-1_LST-3_LST-4': [1, 2, 5], 'LST-1_LST-4': [1, 5], 'LST-2_LST-3': [3, 2], 'LST-2_LST-3_LST-4': [3, 2, 5], 'LST-2_LST-4': [3, 5], 'LST-3_LST-4': [2, 5]} + +def test_telescope_positions(config_gen, config_gen_4lst): + """ + Simple check on telescope positions + """ + M_LST_position = telescope_positions(config_gen) + LSTs_position = telescope_positions(config_gen_4lst) + assert np.allclose(M_LST_position[1].value ,[-8.09 , 77.13 , 0.78 ] ) + assert np.allclose(M_LST_position[2].value ,[39.3 , -62.55 , -0.97 ] ) + assert np.allclose(M_LST_position[3].value ,[-31.21 , -14.57 , 0.2] ) + assert np.allclose(LSTs_position[1].value ,[-70.93, -52.08, 9.08] ) + assert np.allclose(LSTs_position[3].value ,[-35.27, 66.14, -1.92] ) + assert np.allclose(LSTs_position[2].value ,[75.28, 50.48, -5.22] ) + assert np.allclose(LSTs_position[5].value ,[ 30.91, -64.54, -1.92] ) + assert M_LST_position[2].unit==u.m + assert LSTs_position[5].unit==u.m + + def test_format_object(): """ Simple check on a string @@ -136,6 +157,54 @@ def test_load_train_data_files_exc(temp_train_exc): _ = load_train_data_files(str(temp_train_exc)) +def test_load_train_data_files_tel_p(p_stereo, config_gen): + """ + Check dictionary + """ + + events = load_train_data_files_tel(str(p_stereo[0]),config_gen) + assert list(events.keys()) == [1,2,3] + data = events[2] + assert "off_axis" in data.columns + assert "true_event_class" not in data.columns + + +def test_load_train_data_files_tel_g(gamma_stereo, config_gen): + """ + Check dictionary + """ + + events = load_train_data_files_tel(str(gamma_stereo[0]), config_gen) + assert list(events.keys()) == [1,2,3] + data = events[1] + assert "off_axis" in data.columns + assert "true_event_class" not in data.columns + + +def test_load_train_data_files_tel_off(gamma_stereo, config_gen): + """ + Check off-axis cut + """ + events = load_train_data_files_tel( + str(gamma_stereo[0]), config=config_gen, offaxis_min="0.2 deg", offaxis_max="0.5 deg" + ) + data = events[1] + assert np.all(data["off_axis"] >= 0.2) + assert np.all(data["off_axis"] <= 0.5) + + +def test_load_train_data_files_tel_exc(temp_train_exc, config_gen): + """ + Check on exceptions + """ + with pytest.raises( + FileNotFoundError, + match="Could not find any DL1-stereo data files in the input directory.", + ): + _ = load_train_data_files(str(temp_train_exc), config_gen) + + + def test_load_mc_dl2_data_file(p_dl2, gamma_dl2): """ Checks on default loading diff --git a/magicctapipe/io/tests/test_io_monly.py b/magicctapipe/io/tests/test_io_monly.py index 0823b3f54..8bf1bdd63 100644 --- a/magicctapipe/io/tests/test_io_monly.py +++ b/magicctapipe/io/tests/test_io_monly.py @@ -3,6 +3,7 @@ get_dl2_mean, get_stereo_events, load_train_data_files, + load_train_data_files_tel, load_mc_dl2_data_file, load_irf_files, save_pandas_data_in_table, @@ -124,6 +125,54 @@ def test_load_train_data_files_exc(temp_train_exc): _ = load_train_data_files(str(temp_train_exc)) +def test_load_train_data_files_tel_p(p_stereo_monly, config_gen): + """ + Check dictionary + """ + + events = load_train_data_files_tel(str(p_stereo_monly[0]),config_gen) + assert list(events.keys()) == [2,3] + data = events[2] + assert "off_axis" in data.columns + assert "true_event_class" not in data.columns + + +def test_load_train_data_files_tel_g(gamma_stereo_monly, config_gen): + """ + Check dictionary + """ + + events = load_train_data_files_tel(str(gamma_stereo_monly[0]), config_gen) + assert list(events.keys()) == [2,3] + data = events[3] + assert "off_axis" in data.columns + assert "true_event_class" not in data.columns + + +def test_load_train_data_files_tel_off(gamma_stereo_monly, config_gen): + """ + Check off-axis cut + """ + events = load_train_data_files_tel( + str(gamma_stereo_monly[0]), config=config_gen, offaxis_min="0.2 deg", offaxis_max="0.5 deg" + ) + data = events[2] + assert np.all(data["off_axis"] >= 0.2) + assert np.all(data["off_axis"] <= 0.5) + + +def test_load_train_data_files_tel_exc(temp_train_exc, config_gen): + """ + Check on exceptions + """ + with pytest.raises( + FileNotFoundError, + match="Could not find any DL1-stereo data files in the input directory.", + ): + _ = load_train_data_files(str(temp_train_exc), config_gen) + + + def test_load_mc_dl2_data_file(p_dl2_monly, gamma_dl2_monly): """ Checks on default loading diff --git a/magicctapipe/resources/test_config_calib.yaml b/magicctapipe/resources/test_config_calib.yaml index 22d8b3952..3a9444ca7 100644 --- a/magicctapipe/resources/test_config_calib.yaml +++ b/magicctapipe/resources/test_config_calib.yaml @@ -12,8 +12,8 @@ LST: extra_noise_in_bright_pixels: 2.08 increase_psf: - use: false - fraction: null + use: true + fraction: 1.1 tailcuts_clean: picture_thresh: 8 @@ -31,7 +31,7 @@ LST: threshold: 267 fraction: 0.03 - use_only_main_island: false + use_only_main_island: true MAGIC: diff --git a/magicctapipe/resources/test_config_preoffset.yaml b/magicctapipe/resources/test_config_preoffset.yaml new file mode 100644 index 000000000..0e2ddbbf3 --- /dev/null +++ b/magicctapipe/resources/test_config_preoffset.yaml @@ -0,0 +1,77 @@ +mc_tel_ids: + LST-1: 1 + LST-2: 0 + LST-3: 0 + LST-4: 0 + MAGIC-I: 2 + MAGIC-II: 3 + +LST: + image_extractor: + type: "LocalPeakWindowSum" + window_shift: 4 + window_width: 8 + + increase_nsb: + use: true + extra_noise_in_dim_pixels: 1.27 + extra_bias_in_dim_pixels: 0.665 + transition_charge: 8 + extra_noise_in_bright_pixels: 2.08 + + increase_psf: + use: false + fraction: null + + tailcuts_clean: + picture_thresh: 8 + boundary_thresh: 4 + keep_isolated_pixels: false + min_number_picture_neighbors: 2 + + time_delta_cleaning: + use: true + min_number_neighbors: 1 + time_limit: 2 + + dynamic_cleaning: + use: true + threshold: 267 + fraction: 0.03 + + use_only_main_island: false + + +MAGIC: + image_extractor: + type: "SlidingWindowMaxSum" + window_width: 5 + apply_integration_correction: false + + charge_correction: + use: true + factor: 1.143 + + magic_clean: + use_time: true + use_sum: true + picture_thresh: 6 + boundary_thresh: 3.5 + max_time_off: 4.5 + max_time_diff: 1.5 + find_hotpixels: true + pedestal_type: "from_extractor_rndm" + + muon_ring: + thr_low: 25 + tailcut: [12, 8] + ring_completeness_threshold: 25 + +event_coincidence: + timestamp_type_lst: "dragon_time" # select "dragon_time", "tib_time" or "ucts_time" + window_half_width: "300 ns" + pre_offset_search: true + n_pre_offset_search_events: 100 + time_offset: + start: "-10 us" + stop: "0 us" diff --git a/magicctapipe/scripts/lst1_magic/tests/test_coincidence.py b/magicctapipe/scripts/lst1_magic/tests/test_coincidence.py new file mode 100644 index 000000000..4ca01f7a3 --- /dev/null +++ b/magicctapipe/scripts/lst1_magic/tests/test_coincidence.py @@ -0,0 +1,44 @@ +import subprocess +import glob + + +def test_coincidence_preoffset(dl1_lst, merge_magic, temp_coinc_preoff, config_preoff): + """ + Coincidence pre-offset option + """ + + for file in dl1_lst: + _ = subprocess.run( + [ + "lst1_magic_event_coincidence", + f"-l{str(file)}", + f"-m{str(merge_magic)}", + f"-o{str(temp_coinc_preoff)}", + f"-c{str(config_preoff)}", + ] + ) + + assert len(glob.glob(f"{temp_coinc_preoff}/*.h5"))==1 + + +def test_coincidence(dl1_lst, merge_magic, temp_coinc, config): + """ + Coincidence + """ + + for file in dl1_lst: + _ = subprocess.run( + [ + "lst1_magic_event_coincidence", + f"-l{str(file)}", + f"-m{str(merge_magic)}", + f"-o{str(temp_coinc)}", + f"-c{str(config)}", + ] + ) + + assert len(glob.glob(f"{temp_coinc}/*.h5"))==1 + + + + diff --git a/magicctapipe/scripts/lst1_magic/tests/test_mc_dl0_dl1.py b/magicctapipe/scripts/lst1_magic/tests/test_mc_dl0_dl1.py new file mode 100644 index 000000000..3b3947142 --- /dev/null +++ b/magicctapipe/scripts/lst1_magic/tests/test_mc_dl0_dl1.py @@ -0,0 +1,41 @@ +import subprocess +import glob + + +def test_mc_dl0_dl1(temp_DL1_gamma, dl0_gamma, config): + """ + MC DL0 to DL1 + """ + + for file in dl0_gamma: + _ =subprocess.run( + [ + "lst1_magic_mc_dl0_to_dl1", + f"-i{str(file)}", + f"-o{str(temp_DL1_gamma)}", + f"-c{str(config)}", + ] + ) + + + assert len(glob.glob(f"{temp_DL1_gamma}/*.h5"))==4 + + +def test_mc_dl0_dl1_focal_exc(temp_DL1_gamma_focal_exc, dl0_gamma, config): + """ + MC DL0 to DL1 focal length exception + """ + + for file in dl0_gamma: + _ =subprocess.run( + [ + "lst1_magic_mc_dl0_to_dl1", + f"-i{str(file)}", + f"-o{str(temp_DL1_gamma_focal_exc)}", + f"-c{str(config)}", + "-f abc", + ] + ) + assert len(glob.glob(f"{temp_DL1_gamma_focal_exc}/*.h5"))==0 + +