diff --git a/.github/workflows/CI-linter.yml b/.github/workflows/CI-linter.yml index 4684132..fc89812 100644 --- a/.github/workflows/CI-linter.yml +++ b/.github/workflows/CI-linter.yml @@ -54,6 +54,7 @@ jobs: run: | pip install pylint pip install pre-commit + pip install gammasimtools - name: Pre-commit run: | diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index cbccd60..96b4747 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -22,6 +22,7 @@ repos: [ "-rn", # Only display messages "-sn", # Don't display the score + "--disable=W1203", # Disable Use lazy % formatting in logging functions ] # https://github.com/pre-commit/pre-commit-hooks - repo: https://github.com/pre-commit/pre-commit-hooks diff --git a/README.md b/README.md index 0d0b691..2ebc4d1 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,12 @@ # simtools-extra This repository contains additional tools and scripts that are useful for working with the [simtools](https://github.com/gammasim/simtools/tree/main). + +## Tools to generate or test simulation model files + +### Tests of reference configuration files in simtools + +Allow to test reference configuration files for prod5/prod6 against the original sim_telarray +configuration files and to compare all data tables. +See docstring in [test_reference_configs.py](src/simulation_model/test_reference_configs.py) and +[test_data_tables.py](src/simulation_model/test_data_tables.py) for more information. diff --git a/src/.gitkeep b/src/.gitkeep deleted file mode 100644 index e69de29..0000000 diff --git a/src/simulation_model/cfg/CTA-PROD5-Paranal-Alpha.cfg b/src/simulation_model/cfg/CTA-PROD5-Paranal-Alpha.cfg new file mode 100644 index 0000000..e9de6dd --- /dev/null +++ b/src/simulation_model/cfg/CTA-PROD5-Paranal-Alpha.cfg @@ -0,0 +1,122 @@ +#ifndef TELESCOPE +# define TELESCOPE 0 +#endif + +#ifndef NO_GSL_RNG + random_generator = mt19937 % Faster than ranlux. Will fail if not compiled with -DWITH_GSL_RNG. +#endif + +#define MST_TYPE unknown +#if TELESCOPE == 0 + echo + echo ************************************************ +#ifndef ONLY_SSTS + echo Paranal Alpha layout for prod-5 setup. +# if NECTARCAM != 0 + echo (With all MSTs assumed to be a MST-NectarCam.) +# define MST_TYPE MST-NectarCam +# elif SCT != 0 + echo (With all MSTs assumed to be a SCT.) +# define MST_TYPE SCT +# else + echo (With all MSTs assumed to be a MST-FlashCam.) +# define MST_TYPE MST-FlashCam +# endif + echo (With all SSTs of the unified type.) +#else + echo Paranal SST-only layout for prod-5 setup. +#endif + echo ************************************************ + echo +#else +# if NECTARCAM != 0 +# define MST_TYPE MST-NectarCam +# elif SCT != 0 +# define MST_TYPE SCT +# else +# define MST_TYPE MST-FlashCam +# endif +#endif + +array_config_name = Paranal-Alpha-prod5 +#ifndef ONLY_SSTS +array_config_variant = LST/$(MST_TYPE)/SST at CTA South with Alpha layout (plus extra positions: 4/25/50) +#else +array_config_variant = SSTs-only at CTA South for all positions +#endif +array_config_version = 2022-01-20 + +% What transmission option to use (see CTA-PROD4-site.cfg): +#define CTA_SOUTH 1 +#define ATMOSPHERE_PARANAL 1 +#define LOW_EXTINCTION 1 + +#ifdef ONLY_SSTS + +# include + +#else + +# if TELESCOPE == 0 + +% Global and default configuration for things missing in telescope-specific config. +# include + +# elif TELESCOPE < 5 + +# include + +# elif TELESCOPE < 30 + +% Default MST type here is FlashCam. + +# if NECTARCAM != 0 +# include +# elif SCT != 0 +# include +# else +# include +# endif + +# elif TELESCOPE < 80 + +# include + +# else + +Error Invalid telescope for CTA-PROD5 Paranal Alpha layout configuration. + +# endif + +#endif + +% nsb_autoscale_airmass = 0.7,0.15 +% nsb_scaling_factor = 1.0 + +#ifdef NO_STEREO_TRIGGER + + trigger_telescopes = 1 + array_trigger = none + +#else + + trigger_telescopes = 2 % We apply loose stereo trigger immediately + +# ifndef ONLY_SSTS + % To be more strict we need a matching array trigger definition. + array_trigger = array_trigger_prod5_paranal_alpha.dat +# else + array_trigger = none +# endif + +#endif + +#if defined(EXTPRIM_FILE) +echo Making use of externally produced primary particles from $(EXTPRIM_FILE) +IACT EXTPRIM $(EXTPRIM_FILE) +#endif + +#if defined(TELSAMPLE_PARAM) +echo Non-uniformly sampled core offsets according to $(TELSAMPLE_PARAM) +IACT TELSAMPLE $(TELSAMPLE_PARAM) +#endif diff --git a/src/simulation_model/cfg/array_trigger_prod5_paranal_alpha.dat b/src/simulation_model/cfg/array_trigger_prod5_paranal_alpha.dat new file mode 100644 index 0000000..f6257e8 --- /dev/null +++ b/src/simulation_model/cfg/array_trigger_prod5_paranal_alpha.dat @@ -0,0 +1,10 @@ +# The LST and MST+SST stereo triggers are basically independent but any +# extra triggered MST or SST is included in the readout after a LST stereo +# trigger. Single triggered LSTs are not to be read out after a MST+SST trigger. +# Without 'hardstereo', such events would have to be cleaned of the LST data in +# the analysis stage to emulate the hardware stereo readout. +# The array layout supported here is the 4+25+50 Alpha layout plus optional positions +# or the strict prod-5 baseline layout (4+25+70) without alternate MST camera type. + +Trigger 2 of 1, 2, 3, 4 width 120 hardstereo # LST array (hardware) stereo trigger +Trigger 2 of 5 to 29, 30 to 99 width 400 # Alpha/Baseline MST+SST array (software) stereo trigger diff --git a/src/simulation_model/download_configuration_from_zenodo.py b/src/simulation_model/download_configuration_from_zenodo.py new file mode 100644 index 0000000..ec70636 --- /dev/null +++ b/src/simulation_model/download_configuration_from_zenodo.py @@ -0,0 +1,58 @@ +#!/usr/bin/python3 +r""" + Download reference configurations for prod5/prod6 from Zenodo. + +""" + +import logging +import os +import tarfile +import warnings +from pathlib import Path + +import requests +from dotenv import load_dotenv + +logger = logging.getLogger(__name__) +logging.basicConfig(level=logging.INFO) + + +def download_and_unpack(url, target_directory, filename): + """Download and unpack tarball with sim_telarray configuration.""" + + tar_path = Path(target_directory) / filename + logger.info(f"Downloading sim_telarray configuration from {url}") + response = requests.get(url, stream=True, timeout=10) + response.raise_for_status() + + with open(tar_path, "wb") as file: + for chunk in response.iter_content(chunk_size=8192): + file.write(chunk) + + logger.info(f"Extracting {tar_path} to {target_directory}") + with warnings.catch_warnings(): # warnings; should go away with python 3.12 + warnings.simplefilter("ignore", category=RuntimeWarning) + with tarfile.open(tar_path, "r:*") as tar: + tar.extractall(path=target_directory) + tar_path.unlink() + + +def download_configuration_from_zenodo(): + """Download reference configurations for prod5/prod6 from Zenodo.""" + + load_dotenv(".env") + urls = [ + "https://zenodo.org/records/6218687/files/sim_telarray_config_prod5b.tar.gz?download=1", + "https://zenodo.org/records/14198379/files/sim_telarray_config_prod6.tar.gz?download=1", + ] + target_directory = os.getenv("SIMTOOLS_SIMTEL_PATH") or "/workdir/sim_telarray" + for url in urls: + download_and_unpack( + url=url, + target_directory=target_directory, + filename="download", + ) + + +if __name__ == "__main__": + download_configuration_from_zenodo() diff --git a/src/simulation_model/test_data_tables.py b/src/simulation_model/test_data_tables.py new file mode 100644 index 0000000..7bac18a --- /dev/null +++ b/src/simulation_model/test_data_tables.py @@ -0,0 +1,117 @@ +#!/usr/bin/python3 +r""" + Test data tables for prod5/prod6 from 'Files' collection. + + Requires first do download the reference configurations for prod5/prod6 from Zenodo + using the script 'download_configuration_from_zenodo.py'. + + List of files ignored fine-tuned and might need updates in future. + + To test all files run e.g.,: + + ```bash + python ./test_data_tables.py ../../../simulation-models/simulation-models/model_parameters/Files + ``` + +""" + +import argparse +import logging +import os +from pathlib import Path + +from dotenv import load_dotenv + + +logger = logging.getLogger() +logging.basicConfig(level=logging.INFO) + + +def compare_file(file: Path, simtel_path: Path): + """Compare two files.""" + + with open(file, "r", encoding="utf-8") as file1, open( + simtel_path, "r", encoding="utf-8" + ) as file2: + lines1 = file1.readlines() + lines2 = file2.readlines() + # remove empty lines + lines1 = [line for line in lines1 if line.strip()] + lines2 = [line for line in lines2 if line.strip()] + # remove trailing spaces + lines1 = [line.strip() for line in lines1] + lines2 = [line.strip() for line in lines2] + + if len(lines1) != len(lines2): + logger.error(f"Files {file} and {simtel_path} have different number of lines") + return + + for line1, line2 in zip(lines1, lines2): + if line1 != line2: + logger.error(f"Files {file} and {simtel_path} differ") + logger.error(f"Line1: {line1}") + logger.error(f"Line2: {line2}") + + logger.info(f"Files {file} and {simtel_path} are identical") + + +def test_simtel_files(file_directory: str, simtel_path: str): + """Compare all data tables (files) found in 'file_directory'""" + + exclude_list = [ + "ray-tracing-North-LST-1-d10.0-za20.0_validate_optics.ecsv", # simtools-derived + "ray-tracing-North-MST-NectarCam-D-d10.0-za20.0_validate_optics.ecsv", # simtools-derived + "ray-tracing-South-MST-FlashCam-D-d10.0-za20.0_validate_optics.ecsv", # simtools-derived + "ray-tracing-South-SST-D-d10.0-za20.0_validate_optics.ecsv", # simtools-derived + "ray-tracing-South-LST-D-d10.0-za20.0_validate_optics.ecsv", # simtools-derived + "array_coordinates_LaPalma_alpha.dat", + "array_coordinates_Paranal_alpha.dat", + "LaPalma_coords.lis", + "Paranal_coords.lis", + "sct_photon_incidence_angle_focal_surface.ecsv", # simtools-derived + "sst_photon_incidence_angle_secondary_mirror.ecsv", # simtools-derived + "sst_photon_incidence_angle_camera_window.ecsv", # simtools-derived + "sst_photon_incidence_angle_primary_mirror.ecsv", # simtools-derived + "Benn_LaPalma_sky_converted.lis", # needed for testeff and not for nominal simulations + "atm_trans_2200_1_3_0_0_0.dat", # needed for testeff and not for nominal simulations + ] + common_list = [ + "atm_trans_2158_1_3_2_0_0_0.1_0.1.dat", + "atm_trans_2156_1_3_2_0_0_0.1_0.1.dat", + "atm_trans_2147_1_10_2_0_2147.dat", + "funnel_perfect.dat", + "ref_AlSiO2HfO2.dat", + ] + + files = list(Path(file_directory).rglob("*")) + for file in files: + if file.name in exclude_list: + logger.info(f"Skipping {file}") + continue + if file.name in common_list: + simtel_file = simtel_path / "common" / file.name + else: + simtel_file = simtel_path / "CTA" / file.name + + logger.info(f"Comparing {file} with {simtel_file}") + compare_file(file, simtel_file) + + +def main(): + """Main function.""" + load_dotenv(".env") + + parser = argparse.ArgumentParser(description="Directory with files to be tested.") + parser.add_argument("file_directory", type=str, help="Directory with files to be tested.") + args = parser.parse_args() + + simtel_path = os.getenv("SIMTOOLS_SIMTEL_PATH") or "/workdir/sim_telarray" + + test_simtel_files( + file_directory=args.file_directory, + simtel_path=Path(simtel_path) / "sim_telarray/cfg/", + ) + + +if __name__ == "__main__": + main() diff --git a/src/simulation_model/test_reference_simtel_configuration.py b/src/simulation_model/test_reference_simtel_configuration.py new file mode 100644 index 0000000..637fba6 --- /dev/null +++ b/src/simulation_model/test_reference_simtel_configuration.py @@ -0,0 +1,299 @@ +#!/usr/bin/python3 +r""" + Test reference sim_telarray configurations against configurations generated sim_telarray. + + Reference configuration files are in `tests/resources/sim_telarray_configurations/` and are + used in the simtools configuration tests to compare simtools-derived model parameters / + configuration files using the simulation model database with those obtained using sim_telarray. + + These tests compare files with model version 5.0.0 (prod5) and 6.0.0 (prod6) configuration + generated by sim_telarray (no later models are part of the sim_telarray distribution). + + This tools should be executed in the simtools-dev container environment, as it modifies the + 'cfg' directory in the sim_telarray installation: + + 1. start a 'simtools-dev' container (as outlined in simtools/docker/README.md) + 2. download the reference configuration file from Zenodo using the script + `download_configuration_from_zenodo.py` + 2. run this script: + + ``` + python ../simtools-extra/src/simulation_model/test_reference_simtel_configuration.py prod5-south + ``` + + Note: tool gives incorrect answers for prod(5,6)-north CT3 and CT4. + + Allowed options are: 'prod5-north', 'prod5-south', 'prod6-north', 'prod6-south'. + + Original sim_telarray configurations are downloaded from their respective Zenodo repositories + (with the addition of two missing prod5 files provided as part of this repository). + +""" + +import argparse +import logging +import os +import shutil +import subprocess + +from pathlib import Path + +from dotenv import load_dotenv + +logger = logging.getLogger() +logging.basicConfig(level=logging.INFO) + + +configuration = { + "prod6-north": { + "site": "North", + "model_version": "6.0.0", + "cfg_name": "CTA-PROD6-LaPalma.cfg", + "mst_option": "-DNECTARCAM", + "n_telescopes": 30, + "cfg_replace": { + "CTA-PROD6-LST-prototype.cfg": "CTA-North-LSTN-01-6.0.0_test.cfg", + "CTA-PROD6-LST.cfg": "CTA-North-LSTN-02-6.0.0_test.cfg", + "CTA-PROD6-MST-NectarCam.cfg": "CTA-North-MSTN-01-6.0.0_test.cfg", + }, + "telescopes_to_test": ["CT1\t", "CT2\t", "CT3\t", "CT4\t", "CT5\t"], + }, + "prod6-south": { + "site": "South", + "model_version": "6.0.0", + "cfg_name": "CTA-PROD6-Paranal.cfg", + "mst_option": "-DFLASHCAM", + "n_telescopes": 87, + "cfg_replace": { + "CTA-PROD6-LST-prototype.cfg": "CTA-South-LSTS-01-6.0.0_test.cfg", + "CTA-PROD6-MST-FlashCam.cfg": "CTA-South-MSTS-01-6.0.0_test.cfg", + "CTA-PROD6-MST-NectarCam.cfg": "CTA-South-MSTS-01-6.0.0_test.cfg", + "CTA-PROD6-SST.cfg": "CTA-South-SSTS-01-6.0.0_test.cfg", + "CTA-PROD6-SCT.cfg": "CTA-South-SCTS-01-6.0.0_test.cfg", + }, + "telescopes_to_test": ["CT1\t", "CT2\t", "CT3\t", "CT4\t", "CT5\t", "CT30\t"], + }, + "prod5-north": { + "site": "North", + "model_version": "5.0.0", + "cfg_name": "CTA-PROD5-LaPalma-baseline.cfg", + "mst_option": "-DNECTARCAM", + "n_telescopes": 30, + "cfg_replace": { + "CTA-PROD4-LST-prototype.cfg": "CTA-North-LSTN-01-5.0.0_test.cfg", + "CTA-PROD4-LST.cfg": "CTA-North-LSTN-02-5.0.0_test.cfg", + "CTA-PROD4-MST-NectarCam.cfg": "CTA-North-MSTN-01-5.0.0_test.cfg", + }, + "telescopes_to_test": ["CT1\t", "CT2\t", "CT3\t", "CT4\t", "CT5\t"], + }, + "prod5-south": { + "site": "South", + "model_version": "5.0.0", + "cfg_name": "CTA-PROD5-Paranal-Alpha.cfg", + "mst_option": "-DFLASHCAM", + "n_telescopes": 79, + "cfg_replace": { + "CTA-PROD4-LST.cfg": "CTA-South-LSTS-01-5.0.0_test.cfg", + "CTA-PROD4-MST-FlashCam.cfg": "CTA-South-MSTS-01-5.0.0_test.cfg", + "CTA-PROD4-MST-NectarCam.cfg": "CTA-South-MSTS-01-5.0.0_test.cfg", + "CTA-PROD5-SST.cfg": "CTA-South-SSTS-01-5.0.0_test.cfg", + }, + "telescopes_to_test": ["CT1\t", "CT2\t", "CT3\t", "CT4\t", "CT5\t", "CT79\t"], + }, +} + + +def _generate_config(production, new_cfg_name=None): + """ + Generate a sim_telarray configuration file. + + Parameters + ---------- + cfg_name: str + Configuration file name. + + Returns + ------- + list + List of lines from the generated configuration file. + + """ + cfg_name = configuration[production]["cfg_name"] if new_cfg_name is None else new_cfg_name + command = [ + "./sim_telarray/bin/sim_telarray", + "-c", + f"sim_telarray/cfg/CTA/{cfg_name}", + "-C", + "limits=no-internal", + "-C", + "initlist=no-internal", + "-C", + "list=no-internal", + "-C", + "typelist=no-internal", + "-C", + f"maximum_telescopes={configuration[production]['n_telescopes']}", + "-DNSB_AUTOSCALE", + f"{configuration[production]['mst_option']}", + "-DHYPER_LAYOUT", + f"-DNUM_TELESCOPES={configuration[production]['n_telescopes']}", + "/dev/null", + ] + logger.info(f"Command to generate sim_telarray configuration: {(' ').join(command)}") + with subprocess.Popen( + command, + cwd=os.getenv("SIMTOOLS_SIMTEL_PATH"), + stderr=subprocess.DEVNULL, + stdout=subprocess.PIPE, + universal_newlines=True, + ) as process: + stdout, _ = process.communicate() + + # make sure that all telescopes are found + for telescope in configuration[production]["telescopes_to_test"]: + if telescope not in stdout: + raise ValueError(f"Telescope {telescope} not found in configuration") + return [ + line + for line in stdout.splitlines() + if any(telescope in line for telescope in configuration[production]["telescopes_to_test"]) + ] + + +def _modify_cfg(production): + """ + Modify original sim_telarray array configuration file and to use simtools files. + + - copy original sim_telarray array configuration file to a new file + - replace telescope configuration files within this files + - copy simtools-generated configuration files to sim_telarray cfg directory + (from 'tests/resources/sim_telarray_configurations') + + Returns + ------- + str + Name of the new (modified) configuration file. + + """ + simtel_cfg_path = Path(os.getenv("SIMTOOLS_SIMTEL_PATH")) / "sim_telarray/cfg/CTA/" + + # modify main cfg file + cfg_name = configuration[production]["cfg_name"] + replacements = configuration[production]["cfg_replace"] + new_cfg_name = f"new_{cfg_name}" + with open(simtel_cfg_path / cfg_name, encoding="utf-8") as file: + lines = file.readlines() + new_cfg_name = f"new_{cfg_name}" + with open(simtel_cfg_path / new_cfg_name, "w", encoding="utf-8") as file: + for line in lines: + for old, new in replacements.items(): + if old in line: + line = line.replace(old, new) + file.write(line) + + # copy simtools-generated cfg files + files_to_copy = Path("./tests/resources/sim_telarray_configurations").glob("*.cfg") + for file in files_to_copy: + logger.info(f"Copying {file} to {simtel_cfg_path}") + shutil.copy(file, simtel_cfg_path) + + return new_cfg_name + + +def test_reference_simtel_configuration(production): + """Test reference sim_telarray configuration files.""" + + logger.info( + "Testing reference sim_telarray configuration files for " + f"{configuration[production]['cfg_name']} with " + f" {configuration[production]['n_telescopes']} telescopes and " + f"model version {configuration[production]['model_version']} " + f"(MST option: {configuration[production]['mst_option']})." + ) + + simtel_cfg = _generate_config(production) + simtools_cfg = _generate_config(production, _modify_cfg(production)) + + _compare_configuration_files(simtel_cfg, simtools_cfg) + + +def _compare_configuration_files(simtel_cfg, simtools_cfg): + """ + Compare sim_telarray and simtools configurations. + + Includes a list of lines to be ignored as they are process dependent. + """ + + if len(simtel_cfg) == 0 or len(simtools_cfg) == 0: + raise ValueError( + "No configuration file found " + f"(simtel: {len(simtel_cfg)}; simtools: {len(simtools_cfg)})." + ) + + if len(simtel_cfg) != len(simtools_cfg): + raise ValueError( + "Number of telescopes in sim_telarray and simtools configuration files do not match." + ) + + ignore_keywords = { + "CONFIG_RELEASE", + "CONFIG_VERSION", + "OPTICS_CONFIG_NAME", + "OPTICS_CONFIG_VARIANT", + "CAMERA_CONFIG_NAME", + "CAMERA_CONFIG_VARIANT", + "TAILCUT_SCALE", + "PULSE_ANALYSIS", + "SUM_BEFORE_PEAK", + "SUM_AFTER_PEAK", + } + + n_par_differ = 0 + for simtel_line, simtools_line in zip(simtel_cfg, simtools_cfg): + if any(keyword in simtel_line for keyword in ignore_keywords): + continue + if simtel_line != simtools_line: + logger.error("Lines differ:") + logger.error(f"\t SIMTEL : {simtel_line}") + logger.error(f"\t SIMTOOLS:, {simtools_line}") + n_par_differ += 1 + + # Test site parameter (one-by-one; as cfg files are different in structure for sites) + site_parameters = { + "atmospheric_transmission", + } + for par in site_parameters: + simtel_par = [line for line in simtel_cfg if par.upper() in line.upper()] + simtools_par = [line for line in simtools_cfg if par.upper() in line.upper()] + if simtel_par != simtools_par: + logger.error("Site parameter differ:") + logger.error(f"\t SIMTEL : {simtel_par}") + logger.error(f"\t SIMTOOLS: {simtools_par}") + n_par_differ += 1 + + logger.info(f"Number of different parameters: {n_par_differ}") + + +def main(): + """Main function.""" + load_dotenv(".env") + + parser = argparse.ArgumentParser(description="Test reference sim_telarray configuration files.") + parser.add_argument( + "production", choices=configuration.keys(), help="Production configuration to test." + ) + args = parser.parse_args() + production = args.production + + # Prod5 files missing in Zenodo archive + files_to_copy = {"cfg/CTA-PROD5-Paranal-Alpha.cfg", "cfg/array_trigger_prod5_paranal_alpha.dat"} + cfg_dir = Path(os.getenv("SIMTOOLS_SIMTEL_PATH")) / "sim_telarray/cfg/CTA/" + for file in files_to_copy: + logger.info(f"Copying {file} to {cfg_dir}") + shutil.copy(Path(__file__).parent / file, cfg_dir) + + test_reference_simtel_configuration(production) + + +if __name__ == "__main__": + main()