From 801ae547a8437286b9713aee228d4e15a49e3394 Mon Sep 17 00:00:00 2001 From: Jarred Green Date: Thu, 24 Oct 2024 16:12:03 +0200 Subject: [PATCH 1/3] Make `followup` module more robust (Addresses #37) --- gravitational_wave_toy/followup.py | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/gravitational_wave_toy/followup.py b/gravitational_wave_toy/followup.py index 97bdeff..ab6ad78 100644 --- a/gravitational_wave_toy/followup.py +++ b/gravitational_wave_toy/followup.py @@ -66,14 +66,22 @@ def extrapolate_obs_time( if delay < min(event_dict.keys()): res["error_message"] = f"Minimum delay is {min(event_dict.keys())} seconds for this simulation" res["obs_time"] = -1 + raise ValueError(f"Minimum delay is {min(event_dict.keys())} seconds for this simulation [{delay}s requested]") + elif delay > max(event_dict.keys()): + print(f"Warning: delay is greater than maximum delay of {max(event_dict.keys())}s for this simulation [{delay}s requested], value will be extrapolated.") # remove negative values pos_event_dict = {k: v for k, v in event_dict.items() if v > 0} + + if not pos_event_dict: + res["error_message"] = f"Event is never detectable under the observation conditions {filters}" + res["obs_time"] = -1 + return res # perform log interpolation log_event_dict = {log10(k): log10(v) for k, v in pos_event_dict.items()} - interp = interp1d(list(log_event_dict.keys()), list(log_event_dict.values()), kind="linear") + interp = interp1d(list(log_event_dict.keys()), list(log_event_dict.values()), kind="linear", fill_value="extrapolate") try: res["obs_time"] = 10**interp(log10(delay)) @@ -187,12 +195,22 @@ def get_exposure( obs_time = obs_info["obs_time"] if obs_time > 0: - obs_time = round(obs_time / target_precision.value) * target_precision + if obs_time > max_time.value: + obs_info["error_message"] = f"Exposure time of {int(obs_time)} s exceeds maximum time" + obs_time = -1 + else: + obs_time = round(obs_time / target_precision.value) * target_precision # rename key - obs_info["angle"] = obs_info.pop("theta_view") + obs_info["angle"] = obs_info.pop("theta_view") * u.deg obs_info["ebl_model"] = obs_info.pop("irf_ebl_model") + # add other units + obs_info["long"] = obs_info["long"] * u.deg + obs_info["lat"] = obs_info["lat"] * u.deg + obs_info["eiso"] = obs_info["eiso"] * u.erg + obs_info["dist"] = obs_info["dist"] * u.kpc + other_info = { "min_energy": min_energy, "max_energy": max_energy, From 90892d5397edbb039d92e585082b2cbd45eaa9fb Mon Sep 17 00:00:00 2001 From: Jarred Green Date: Thu, 24 Oct 2024 16:18:05 +0200 Subject: [PATCH 2/3] Set units of lat and long to radians --- gravitational_wave_toy/followup.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gravitational_wave_toy/followup.py b/gravitational_wave_toy/followup.py index ab6ad78..c7037f6 100644 --- a/gravitational_wave_toy/followup.py +++ b/gravitational_wave_toy/followup.py @@ -206,8 +206,8 @@ def get_exposure( obs_info["ebl_model"] = obs_info.pop("irf_ebl_model") # add other units - obs_info["long"] = obs_info["long"] * u.deg - obs_info["lat"] = obs_info["lat"] * u.deg + obs_info["long"] = obs_info["long"] * u.rad + obs_info["lat"] = obs_info["lat"] * u.rad obs_info["eiso"] = obs_info["eiso"] * u.erg obs_info["dist"] = obs_info["dist"] * u.kpc From 235c8e2efecf90ed8e8155841aa8e0b9088029d0 Mon Sep 17 00:00:00 2001 From: Jarred Green Date: Thu, 24 Oct 2024 16:24:34 +0200 Subject: [PATCH 3/3] Add followup example to the readme --- README.md | 52 +++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 51 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 948037a..4718db6 100644 --- a/README.md +++ b/README.md @@ -27,6 +27,8 @@ - [Authors](#️authors) - [Instructions](#instructions) - [GW Observations](#instructions-gw-obs) + - [Example Simulation](#example-simulation) + - [Example followup calculation](#example-followup) - [Reading Results](#instructions-reading) - [Plotting Heatmaps](#instructions-plotting) @@ -91,7 +93,7 @@ This code simulates observations of simulated gravitational wave events to deter - a dictionary object containing detection information and parameters of the event itself (extracted from the model) -#### Example +#### Example Simulation ```python from astropy import units as u @@ -146,6 +148,54 @@ print(f"Observation time at delay={delay_time} is {res_ebl['obs_time']} with EBL # Obs time at delay=1800.0 s is 1292.0 s with EBL=franceschini ``` +### Example followup calculation + +```python +import astropy.units as u +import pandas as pd +from gravitational_wave_toy import followup + +lookup_talbe = "./O5_gammapy_observations_v4.parquet" + +# optional, but it's recommended to load the DataFrame first save time +# otherwise you can directly pass the filepath to the get_exposure method +lookup_df = pd.read_parquet(lookup_talbe) + +event_id = 1 +delay = 10 * u.s +site = "north" +zenith = 60 +ebl = "franceschini" + + +followup.get_exposure( + event_id=event_id, + delay=delay, + site=site, + zenith=zenith, + extrapolation_df=lookup_df, + ebl=ebl, +) + +# returns, e.g. +# { +# 'long': , +# 'lat': , +# 'eiso': , +# 'dist': , +# 'obs_time': , +# 'error_message': '', +# 'angle': , +# 'ebl_model': 'franceschini', +# 'min_energy': , +# 'max_energy': , +# 'seen': True, +# 'start_time': , +# 'end_time': , +# 'id': 4 +# } +``` + ### Reading Results Note: The most recent simulations for CTA for the O5 observing run are [available on the CTA XWiki](https://cta.cloud.xwiki.com/xwiki/wiki/phys/view/Transients%20WG/Chasing%20the%20counterpart%20of%20GW%20alerts%20with%20CTA/O5%20Observation%20times%20with%20gw-toy%20package/) for CTA members.