show what was thrown

plenoirf/summary/scripts/0020_what_was_thrown.py → plenoirf/summary/scripts/0020_thrown_primary_directions.py

@@ -21,16 +21,6 @@
     os.path.join(paths["analysis_dir"], "0005_common_binning", "energy.json")
 )["trigger_acceptance"]
 
-passing_trigger = json_utils.tree.read(
-    os.path.join(paths["analysis_dir"], "0055_passing_trigger")
-)
-
-nat_bin = binning_utils.Binning(
-    bin_edges=np.geomspace(1, 100_000, energy_bin["num_bins"])
-)
-
-POINTNIG_ZENITH_BIN = res.ZenithBinning("once")
-
 hh = spherical_histogram.HemisphereHistogram(
     num_vertices=12_000,
     max_zenith_distance_rad=np.deg2rad(90),

plenoirf/summary/scripts/0021_thrown_primary_energy.py

@@ -0,0 +1,111 @@
+#!/usr/bin/python
+import sys
+import numpy as np
+import plenoirf as irf
+import sparse_numeric_table as snt
+import os
+import cosmic_fluxes
+import json_utils
+import spherical_coordinates
+import solid_angle_utils
+import binning_utils
+import spherical_histogram
+import sebastians_matplotlib_addons as seb
+
+
+paths = irf.summary.paths_from_argv(sys.argv)
+res = irf.summary.Resources.from_argv(sys.argv)
+os.makedirs(paths["out_dir"], exist_ok=True)
+seb.matplotlib.rcParams.update(res.analysis["plot"]["matplotlib"])
+
+energy_bin = json_utils.read(
+    os.path.join(paths["analysis_dir"], "0005_common_binning", "energy.json")
+)["trigger_acceptance"]
+
+POINTNIG_ZENITH_BIN = res.ZenithBinning("once")
+
+eee = {}
+for pk in res.PARTICLES:
+    with res.open_event_table(particle_key=pk) as arc:
+        table = arc.read_table(
+            levels_and_columns={
+                "primary": [snt.IDX, "energy_GeV", "azimuth_rad", "zenith_rad"]
+            }
+        )
+
+    eee[pk] = {}
+    eee[pk]["bin_counts"] = np.zeros(
+        shape=(POINTNIG_ZENITH_BIN.num, energy_bin["num_bins"])
+    )
+
+    for zdbin in range(POINTNIG_ZENITH_BIN.num):
+        zenith_start_rad = POINTNIG_ZENITH_BIN.edges[zdbin]
+        zenith_stop_rad = POINTNIG_ZENITH_BIN.edges[zdbin + 1]
+
+        mask = np.logical_and(
+            table["primary"]["zenith_rad"] >= zenith_start_rad,
+            table["primary"]["zenith_rad"] < zenith_stop_rad,
+        )
+
+        eee[pk]["bin_counts"][zdbin, :] = np.histogram(
+            table["primary"]["energy_GeV"][mask], bins=energy_bin["edges"]
+        )[0]
+
+
+# differential SED style histogram
+# --------------------------------
+linestyles = ["-", "--", ":"]
+fig = seb.figure(irf.summary.figure.FIGURE_STYLE)
+ax = seb.add_axes(fig=fig, span=irf.summary.figure.AX_SPAN)
+for pk in res.PARTICLES:
+    for zdbin in range(POINTNIG_ZENITH_BIN.num):
+        ax.plot(
+            energy_bin["centers"],
+            (eee[pk]["bin_counts"][zdbin, :] / energy_bin["widths"])
+            * energy_bin["centers"] ** (1.5),
+            linestyle=linestyles[zdbin],
+            color=res.PARTICLE_COLORS[pk],
+        )
+res.ax_add_site_marker(ax)
+ax.loglog()
+ax.set_xlim(energy_bin["limits"])
+ax.set_xlabel("energy / GeV")
+ax.set_ylabel(
+    r"(energy)$^{1.5}$ differential intensity /"
+    + "\n"
+    + r"(GeV)$^{1.5}$ (GeV)$^{-1}$"
+)
+fig.savefig(
+    os.path.join(
+        paths["out_dir"],
+        "thrown_primary_energy_differential_spectral_energy_distribution.jpg",
+    )
+)
+seb.close(fig)
+
+# simple histogram
+# ----------------
+fig = seb.figure(irf.summary.figure.FIGURE_STYLE)
+ax = seb.add_axes(fig=fig, span=irf.summary.figure.AX_SPAN)
+for pk in res.PARTICLES:
+    for zdbin in range(POINTNIG_ZENITH_BIN.num):
+        seb.ax_add_histogram(
+            ax=ax,
+            bin_edges=energy_bin["edges"],
+            bincounts=eee[pk]["bin_counts"][zdbin, :],
+            linestyle=linestyles[zdbin],
+            linecolor=res.PARTICLE_COLORS[pk],
+            linealpha=1.0,
+            bincounts_upper=None,
+            bincounts_lower=None,
+            face_color=res.PARTICLE_COLORS[pk],
+            face_alpha=0.3,
+            draw_bin_walls=True,
+        )
+res.ax_add_site_marker(ax)
+ax.loglog()
+ax.set_xlim(energy_bin["limits"])
+ax.set_xlabel("energy / GeV")
+ax.set_ylabel(r"intensity / 1")
+fig.savefig(os.path.join(paths["out_dir"], "thrown_primary_energy.jpg"))
+seb.close(fig)

plenoirf/summary/scripts/0025_drawn_instrument_pointing.py

@@ -0,0 +1,124 @@
+#!/usr/bin/python
+import sys
+import numpy as np
+import plenoirf as irf
+import sparse_numeric_table as snt
+import os
+import json_utils
+import spherical_coordinates
+import solid_angle_utils
+import binning_utils
+import spherical_histogram
+import sebastians_matplotlib_addons as seb
+
+
+paths = irf.summary.paths_from_argv(sys.argv)
+res = irf.summary.Resources.from_argv(sys.argv)
+os.makedirs(paths["out_dir"], exist_ok=True)
+seb.matplotlib.rcParams.update(res.analysis["plot"]["matplotlib"])
+
+energy_bin = json_utils.read(
+    os.path.join(paths["analysis_dir"], "0005_common_binning", "energy.json")
+)["trigger_acceptance"]
+
+hh = spherical_histogram.HemisphereHistogram(
+    num_vertices=12_000,
+    max_zenith_distance_rad=np.deg2rad(90),
+)
+
+cmap = seb.plt.colormaps["inferno"].resampled(256)
+
+for pk in res.PARTICLES:
+    with res.open_event_table(particle_key=pk) as arc:
+        table = arc.read_table(
+            levels_and_columns={"instrument_pointing": "__all__"}
+        )
+
+    hh.reset()
+    hh.assign_azimuth_zenith(
+        azimuth_rad=table["instrument_pointing"]["azimuth_rad"],
+        zenith_rad=table["instrument_pointing"]["zenith_rad"],
+    )
+
+    pointing_range_color = "yellowgreen"
+    fig = seb.figure(irf.summary.figure.FIGURE_STYLE)
+    ax = seb.add_axes(fig=fig, span=[0.0, 0.0, 1, 1], style=seb.AXES_BLANK)
+    ax_legend = seb.add_axes(
+        fig=fig, span=[0.8, 0.05, 0.2, 0.9], style=seb.AXES_BLANK
+    )
+    ax_legend.set_xlim([0, 1])
+    ax_legend.set_ylim([0, 1])
+    ax_colorbar = seb.add_axes(
+        fig=fig,
+        span=[0.8, 0.4, 0.02, 0.5],
+    )
+
+    vx, vy, vz = hh.bin_geometry.vertices.T
+    vaz, vzd = spherical_coordinates.cx_cy_cz_to_az_zd(vx, vy, vz)
+    faces_counts_per_solid_angle = (
+        hh.bin_counts / hh.bin_geometry.faces_solid_angles
+    )
+    faces_counts_per_solid_angle_p99 = np.percentile(
+        faces_counts_per_solid_angle, q=99
+    )
+    faces_colors = cmap(
+        faces_counts_per_solid_angle / faces_counts_per_solid_angle_p99
+    )
+    seb.hemisphere.ax_add_faces(
+        ax=ax,
+        azimuths_rad=vaz,
+        zeniths_rad=vzd,
+        faces=hh.bin_geometry.faces,
+        faces_colors=faces_colors,
+    )
+    ax.set_aspect("equal")
+    ptg_circ_az = np.linspace(0, np.pi * 2, 360)
+    ptg_circ_zd = res.config["pointing"]["range"][
+        "max_zenith_distance_rad"
+    ] * np.ones(shape=ptg_circ_az.shape)
+    seb.hemisphere.ax_add_plot(
+        ax=ax,
+        azimuths_rad=ptg_circ_az,
+        zeniths_rad=ptg_circ_zd,
+        linewidth=0.66,
+        color=pointing_range_color,
+    )
+    seb.hemisphere.ax_add_grid_stellarium_style(
+        ax=ax, color="grey", linewidth=0.33
+    )
+    seb.hemisphere.ax_add_ticklabel_text(ax=ax, color="grey")
+    ax_legend.plot(
+        [0.0, 0.1],
+        [0.1, 0.1],
+        color=pointing_range_color,
+        linewidth=0.66,
+    )
+    ax_legend.text(
+        x=0.15,
+        y=0.1,
+        s="pointing\nrange",
+        fontsize=9,
+        verticalalignment="center",
+    )
+    ax_legend.text(
+        x=0.0,
+        y=0.2,
+        s=f"site: {res.SITE['name']:s}",
+        fontsize=9,
+        verticalalignment="center",
+    )
+    ax_legend.text(
+        x=0.0,
+        y=0.3,
+        s=f"particle: {pk:s}",
+        fontsize=9,
+        verticalalignment="center",
+    )
+    _norm = seb.plt_colors.Normalize(
+        vmin=0.0, vmax=faces_counts_per_solid_angle_p99, clip=False
+    )
+    _mappable = seb.plt.cm.ScalarMappable(norm=_norm, cmap=cmap)
+    seb.plt.colorbar(mappable=_mappable, cax=ax_colorbar)
+    ax_colorbar.set_ylabel(r"intensity / sr$^{-1}$")
+    fig.savefig(os.path.join(paths["out_dir"], f"{pk:s}_pointings_drawn.jpg"))
+    seb.close(fig)