Coverage for src/km3modules/plot.py: 14%
183 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-04-23 03:15 +0000
1# Filename: plot.py
2# -*- coding: utf-8 -*-
3# pylint: disable=locally-disabled
4"""
5A collection of plotting functions and modules.
7"""
8from collections import Counter
9from datetime import datetime
10import gc
11import os
12import multiprocessing as mp
13import time
15import numpy as np
16import matplotlib.pyplot as plt # noqa
17import matplotlib.ticker as ticker
18from matplotlib import pylab # noqa
19import km3db
21import km3pipe as kp # noqa
22import km3pipe.style # noqa
24from km3modules.hits import count_multiplicities
27def plot_dom_parameters(
28 data,
29 detector,
30 filename,
31 label,
32 title,
33 vmin=0.0,
34 vmax=10.0,
35 cmap="cividis",
36 under="deepskyblue",
37 over="deeppink",
38 underfactor=1.0,
39 overfactor=1.0,
40 missing="lightgray",
41 hide_limits=False,
42):
43 """Creates a plot in the classical monitoring.km3net.de style.
45 Parameters
46 ----------
47 data: dict((du, floor) -> value)
48 detector: km3pipe.hardware.Detector() instance
49 filename: filename or filepath
50 label: str
51 title: str
52 underfactor: a scale factor for the points used for underflow values
53 overfactor: a scale factor for the points used for overflow values
54 hide_limits: do not show under/overflows in the plot
56 """
57 x, y, _ = zip(*detector.doms.values())
58 fig, ax = plt.subplots(figsize=(10, 6))
59 cmap = plt.get_cmap(cmap).copy()
60 cmap.set_over(over, 1.0)
61 cmap.set_under(under, 1.0)
63 m_size = 100
64 scatter_args = {
65 "edgecolors": "None",
66 "vmin": vmin,
67 "vmax": vmax,
68 }
69 sc_inactive = ax.scatter(
70 x, y, c=missing, label="missing", s=m_size * 0.9, **scatter_args
71 )
73 xa, ya = map(np.array, zip(*data.keys()))
74 zs = np.array(list(data.values()))
75 in_range_idx = np.logical_and(zs >= vmin, zs <= vmax)
76 sc = ax.scatter(
77 xa[in_range_idx],
78 ya[in_range_idx],
79 c=zs[in_range_idx],
80 cmap=cmap,
81 s=m_size,
82 **scatter_args
83 )
84 if not hide_limits:
85 under_idx = zs < vmin
86 ax.scatter(
87 xa[under_idx],
88 ya[under_idx],
89 c=under,
90 label="< {0}".format(vmin),
91 s=m_size * underfactor,
92 **scatter_args
93 )
94 over_idx = zs > vmax
95 ax.scatter(
96 xa[over_idx],
97 ya[over_idx],
98 c=over,
99 label="> {0}".format(vmax),
100 s=m_size * overfactor,
101 **scatter_args
102 )
104 cb = plt.colorbar(sc)
105 cb.set_label(label)
107 ax.set_title("{0}\n{1} UTC".format(title, datetime.utcnow().strftime("%c")))
108 ax.set_xlabel("DU")
109 ax.set_ylabel("DOM")
110 ax.set_ylim(-2)
111 ax.set_yticks(range(1, 18 + 1))
112 major_locator = pylab.MaxNLocator(integer=True)
113 sc_inactive.axes.xaxis.set_major_locator(major_locator)
115 ax.legend(
116 bbox_to_anchor=(0.0, -0.16, 1.0, 0.102),
117 loc=1,
118 ncol=2,
119 mode="expand",
120 borderaxespad=0.0,
121 )
123 fig.tight_layout()
125 plt.savefig(filename, dpi=120, bbox_inches="tight")
126 plt.close("all")
127 gc.collect()
130def make_dom_map(pmt_directions, values, nside=512, d=0.2, smoothing=0.1):
131 """Create a mollweide projection of a DOM with given PMTs.
133 The output can be used to call the `healpy.mollview` function.
134 """
135 import healpy as hp
137 discs = [hp.query_disc(nside, dir, 0.2) for dir in pmt_directions]
138 npix = hp.nside2npix(nside)
139 pixels = np.zeros(npix)
140 for disc, value in zip(discs, values):
141 for d in disc:
142 pixels[d] = value
143 if smoothing > 0:
144 return hp.sphtfunc.smoothing(pixels, fwhm=smoothing, iter=1)
145 return pixels
148class IntraDOMCalibrationPlotter(kp.Module):
149 def configure(self):
150 self.plots_path = self.get("plots_path", default=os.getcwd())
151 self.data_path = self.get("data_path", default=os.getcwd())
152 self.det_oid = self.require("det_oid")
153 self.clbmap = km3db.CLBMap(self.det_oid)
155 def process(self, blob):
156 calibration = blob["IntraDOMCalibration"]
157 for process in (self.create_plot, self.save_hdf5):
158 proc = mp.Process(target=process, args=(calibration,))
159 proc.daemon = True
160 proc.start()
161 proc.join()
162 return blob
164 def create_plot(self, calibration):
165 print("Creating plot...")
166 fig, axes = plt.subplots(6, 3, figsize=(16, 20), sharex=True, sharey=True)
167 sorted_dom_ids = sorted(
168 calibration.keys(),
169 key=lambda d: (self.clbmap.dom_ids[d].du, self.clbmap.dom_ids[d].floor),
170 ) # by DU and FLOOR, note that DET OID is needed!
171 for ax, dom_id in zip(axes.flatten(), sorted_dom_ids):
172 calib = calibration[dom_id]
173 ax.plot(np.cos(calib["angles"]), calib["means"], ".")
174 ax.plot(np.cos(calib["angles"]), calib["corrected_means"], ".")
175 du = self.clbmap.dom_ids[dom_id].du
176 floor = self.clbmap.dom_ids[dom_id].floor
177 ax.set_title("{0} - {1}".format("DU{}-DOM{}".format(du, floor), dom_id))
178 ax.set_ylim((-10, 10))
179 plt.suptitle("{0} UTC".format(datetime.utcnow().strftime("%c")))
180 plt.savefig(os.path.join(self.plots_path, "intradom.png"), bbox_inches="tight")
181 plt.close("all")
183 fig, axes = plt.subplots(6, 3, figsize=(16, 20), sharex=True, sharey=True)
184 for ax, dom_id in zip(axes.flatten(), sorted_dom_ids):
185 calib = calibration[dom_id]
186 ax.plot(np.cos(calib["angles"]), calib["rates"], ".")
187 ax.plot(np.cos(calib["angles"]), calib["corrected_rates"], ".")
188 du = self.clbmap.dom_ids[dom_id].du
189 floor = self.clbmap.dom_ids[dom_id].floor
190 ax.set_title("{0} - {1}".format("DU{}-DOM{}".format(du, floor), dom_id))
191 ax.set_ylim((0, 10))
192 plt.suptitle("{0} UTC".format(datetime.utcnow().strftime("%c")))
193 plt.savefig(
194 os.path.join(self.plots_path, "angular_k40rate_distribution.png"),
195 bbox_inches="tight",
196 )
197 plt.close("all")
198 gc.collect()
200 def save_hdf5(self, calibration):
201 print("Saving calibration information...")
202 pd = kp.extras.pandas()
203 store = pd.HDFStore(os.path.join(self.data_path, "k40calib.h5"), mode="a")
204 now = int(time.time())
205 timestamps = (now,) * 31
206 for dom_id, calib in calibration.items():
207 tdc_channels = range(31)
208 t0s = calib["opt_t0s"].x
209 dom_ids = (dom_id,) * 31
210 df = pd.DataFrame(
211 {
212 "timestamp": timestamps,
213 "dom_id": dom_ids,
214 "tdc_channel": tdc_channels,
215 "t0s": t0s,
216 }
217 )
218 store.append("t0s", df, format="table", data_columns=True)
219 store.close()
222def ztplot(
223 hits,
224 filename=None,
225 title=None,
226 max_z=None,
227 figsize=(16, 8),
228 n_dus=4,
229 ytick_distance=200,
230 max_multiplicity_entries=10,
231 grid_lines=[],
232):
233 """Creates a ztplot like shown in the online monitoring"""
234 fontsize = 16
236 dus = set(hits.du)
238 if n_dus is not None:
239 dus = [c[0] for c in Counter(hits.du).most_common(n_dus)]
240 mask = [du in dus for du in hits.du]
241 hits = hits[mask]
243 dus = sorted(dus)
244 doms = set(hits.dom_id)
246 hits = hits.append_columns("multiplicity", np.ones(len(hits))).sorted(by="time")
248 if max_z is None:
249 max_z = int(np.ceil(np.max(hits.pos_z) / 100.0)) * 100 * 1.05
251 for dom in doms:
252 dom_hits = hits[hits.dom_id == dom]
253 mltps, m_ids = count_multiplicities(dom_hits.time)
254 hits["multiplicity"][hits.dom_id == dom] = mltps
256 if np.any(hits.triggered):
257 time_offset = np.min(hits[hits.triggered > 0].time)
258 hits.time -= time_offset
260 n_plots = len(dus)
261 n_cols = int(np.ceil(np.sqrt(n_plots)))
262 n_rows = int(n_plots / n_cols) + (n_plots % n_cols > 0)
263 marker_fig, marker_axes = plt.subplots()
264 # for the marker size hack...
265 fig, axes = plt.subplots(
266 ncols=n_cols,
267 nrows=n_rows,
268 sharex=True,
269 sharey=True,
270 figsize=figsize,
271 constrained_layout=True,
272 )
274 axes = [axes] if n_plots == 1 else trim_axes(axes, n_plots)
276 for ax, du in zip(axes, dus):
277 du_hits = hits[hits.du == du]
278 for grid_line in grid_lines:
279 ax.axhline(grid_line, lw=1, color="b", ls="--", alpha=0.15)
280 trig_hits = du_hits[du_hits.triggered > 0]
282 ax.scatter(
283 du_hits.time,
284 du_hits.pos_z,
285 s=du_hits.multiplicity * 30,
286 c="#09A9DE",
287 label="hit",
288 alpha=0.5,
289 )
290 ax.scatter(
291 trig_hits.time,
292 trig_hits.pos_z,
293 s=trig_hits.multiplicity * 30,
294 alpha=0.8,
295 marker="+",
296 c="#FF6363",
297 label="triggered hit",
298 )
299 ax.set_title("DU{0}".format(int(du)), fontsize=fontsize, fontweight="bold")
301 # The only way I could create a legend with matching marker sizes
302 max_multiplicity = int(np.max(du_hits.multiplicity))
303 markers = list(
304 range(
305 0,
306 max_multiplicity,
307 np.ceil(max_multiplicity / max_multiplicity_entries).astype(int),
308 )
309 )
310 custom_markers = [
311 marker_axes.scatter([], [], s=mult * 30, color="#09A9DE", lw=0, alpha=0.5)
312 for mult in markers
313 ] + [marker_axes.scatter([], [], s=30, marker="+", c="#FF6363")]
314 ax.legend(
315 custom_markers,
316 ["multiplicity"] + [" %d" % m for m in markers[1:]] + ["triggered"],
317 scatterpoints=1,
318 markerscale=1,
319 loc="lower right",
320 frameon=True,
321 framealpha=0.7,
322 )
324 for idx, ax in enumerate(axes):
325 ax.set_ylim(0, max_z)
326 ax.tick_params(labelsize=fontsize)
327 ax.yaxis.set_major_locator(ticker.MultipleLocator(ytick_distance))
328 xlabels = ax.get_xticklabels()
329 for label in xlabels:
330 label.set_rotation(45)
332 if idx % n_cols == 0:
333 ax.set_ylabel("z [m]", fontsize=fontsize)
334 if idx >= len(axes) - n_cols:
335 ax.set_xlabel("time [ns]", fontsize=fontsize)
337 if title is not None:
338 plt.suptitle(title, fontsize=fontsize, y=1.05)
339 if filename is not None:
340 plt.savefig(filename, dpi=120, bbox_inches="tight")
342 gc.collect()
343 return fig
346def trim_axes(axes, n):
347 """little helper to massage the axes list to have correct length..."""
348 axes = axes.flat
349 for ax in axes[n:]:
350 ax.remove()
351 return axes[:n]
354def cumulative_run_livetime(qtable, kind="runs"):
355 """Create a figure which plots the cumulative livetime of runs
357 Parameters
358 ----------
359 qtable: pandas.DataFrame
360 A table which has the run number as index and columns for
361 'livetime_s', 'timestamp' and 'datetime' (pandas datetime).
362 kind: str
363 'runs' to plot for each run or 'timeline' to plot based
364 on the actual run time.
366 Returns
367 -------
368 matplotlib.Figure
369 """
370 fig, ax = plt.subplots()
372 options = {
373 "runs": {"xlabel": "run", "xs": qtable.index},
374 "timeline": {"xlabel": None, "xs": qtable.datetime},
375 }
377 actual_livetime = np.cumsum(qtable["livetime_s"])
378 optimal_livetime = np.cumsum(qtable.timestamp.diff())
380 ax.plot(options[kind]["xs"], actual_livetime, label="actual livetime")
381 ax.plot(options[kind]["xs"], optimal_livetime, label="100% livetime")
382 ax.set_xlabel(options[kind]["xlabel"])
383 ax.set_ylabel("time / s")
384 ax.legend()
386 return fig