Skip to content

Plot injs mod #5243

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
sebastiangomezlopez wants to merge 2 commits into
base: master
Choose a base branch
from
Draft

Plot injs mod #5243

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9dd1e65
Update pycbc_pygrb_plot_injs_results
sebastiangomezlopez Dec 15, 2025
2c3554d
Update pycbc_pygrb_plot_injs_results
sebastiangomezlopez Dec 15, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 26 additions & 9 deletions bin/pygrb/pycbc_pygrb_plot_injs_results
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,9 @@ __program__ = "pycbc_pygrb_plot_injs_results"
# =============================================================================
# Functions
# =============================================================================
def ra_to_ra_mollweide(ra):
return (ra + np.pi) % (2*np.pi) - np.pi

def process_var_strings(qty):
"""Add underscores to match HDF column name conventions"""

Expand Down Expand Up @@ -169,18 +172,19 @@ def complete_inj_data(injs, keys, tag, ifos=[]):
data_dict[key] = complete_incl_data(injs, key, tag)
elif key == 'sky_error':
data_dict[key] = complete_sky_error_data(injs, tag)
elif key == 'eff_dist':
eff_dists = np.empty((injs[tag + '/tc'].size, len(ifos)))
elif key.startswith('eff_dist'):
eff_dist = 0
for i, ifo in enumerate(ifos):
eff_dists[:, i] = Detector(ifo).effective_distance(
antenna = Detector(ifo)
data_dict['eff_dist_'+ifo] = antenna.effective_distance(
injs[tag + '/distance'],
injs[tag + '/ra'],
injs[tag + '/dec'],
injs[tag + '/polarization'],
injs[tag + '/tc'],
injs[tag + '/inclination'],
)
data_dict[key] = 1 / (1 / eff_dists).sum(axis=1)
injs[tag + '/inclination'])
eff_dist += 1.0 / data_dict['eff_dist_'+ifo]
data_dict['eff_dist'] = 1.0 / eff_dist
else:
data_dict[key] = injs[tag+'/'+key]
except KeyError:
Expand Down Expand Up @@ -212,7 +216,8 @@ admitted_vars = easy_keys + ['mtotal', 'q', 'mchirp',
'cosincl', 'absincl', 'cosabsincl',
'sky_error', 'skyerror', 'end_time', 'endtime',
'dec', 'ra', 'coaphase', 'coa_phase',
'eff_site_dist', 'eff_dist',
'eff_site_dist', 'eff_dist', 'eff_dist_H1',
'eff_dist_L1', 'eff_dist_V1', 'eff_dist_sqr_quadsum',
Comment on lines +219 to +220
Copy link
Contributor

@titodalcanton titodalcanton Jan 8, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that we already have "eff_site_dist" to plot the effective distance for a specific site. What is the motivation for adding the explicit per-site ones?

'effsitedist', 'effdist', 'chip', 'chi_p']
admitted_vars = sorted(set(admitted_vars))

Expand Down Expand Up @@ -427,6 +432,10 @@ axis_labels_dict = {'mchirp': "Chirp Mass (solar masses)",
'distance': "Distance (Mpc)",
'eff_site_dist': f"{opts.ifo} effective distance (Mpc)",
'eff_dist': "Inverse sum of inverse effective distances (Mpc)",
'eff_dist_H1': "Eff. dist. Hanford (Mpc)",
'eff_dist_L1': "Eff. dist. Livingston (Mpc)",
'eff_dist_V1': "Eff. dist. Virgo (Mpc)",
'eff_dist_sqr_quadsum': "sqrt of quadrature sum of effective distances (Mpc)",
'end_time': "GPS time (s)",
'sky_error': "Rec. sky error (rad)",
'coa_phase': "Phase of complex SNR (rad)",
Expand All @@ -448,9 +457,13 @@ axis_labels_dict = {'mchirp': "Chirp Mass (solar masses)",
'spin2y': "Spin y-component of 2nd binary component",
'spin2z': "Spin z-component of 2nd binary component",
'chi_p': "Effective precession spin"}

fig = plt.figure()
ax = fig.gca()
if 'dec' in y_qty and 'ra' in x_qty:
ax = fig.add_subplot(111, projection="mollweide")
missed_inj['ra'] = ra_to_ra_mollweide(missed_inj['ra'])
found_inj['ra'] = ra_to_ra_mollweide(found_inj['ra'])
Comment on lines +462 to +464
Copy link
Contributor

@titodalcanton titodalcanton Jan 8, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There are some cases where the all-sky projection plot is detrimental, typically when the input skymap is very precise, as all points are bunched together into a single blob. Perhaps we could have different variable names to control the projection, e.g. "ra"/"dec" and "ra_mollweide"/"dec_mollweide". Another way (probably my preferred one actually) would be to add a command line argument like --sky-projection mollweide that is only used when "ra" and "dec" are being plotted.

else:
ax = fig.add_subplot(111)
ax.set_xscale("log" if opts.x_log else "linear")
ax.set_yscale("log" if opts.y_log else "linear")
ax.set_xlabel(axis_labels_dict[x_qty])
Expand Down Expand Up @@ -607,6 +620,10 @@ if plot_title is None:
'distance': "distance (Mpc)",
'eff_dist': "inverse sum of inverse effective distances",
'eff_site_dist': "site specific effective distance",
'eff_dist_H1': "Eff. dist. Hanford (Mpc)",
'eff_dist_L1': "Eff. dist. Livingston (Mpc)",
'eff_dist_V1': "Eff. dist. Virgo (Mpc)",
'eff_dist_sqr_quadsum': "sqrt of quadrature sum of effective distances (Mpc)",
'end_time': "time",
'coa_phase': "phase of complex SNR",
'dec': "declination",
Expand Down
Loading