Coverage for polars_analysis / plotting / pedestal_plotting.py: 84%
438 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-04-20 15:59 -0400
1import logging
2import pathlib
3from itertools import product
4from pathlib import Path
5from typing import Any, Dict, List, Literal, Optional, Set, Union
7import matplotlib
8from typing_extensions import deprecated
10matplotlib.use("agg")
11import matplotlib.colors as mcolors
12import matplotlib.pyplot as plt
13import matplotlib.ticker as ticker
14import numpy as np
15import polars as pl
16import scipy
17from diptest import diptest # type: ignore
18from scipy.stats import gamma
20from polars_analysis.analysis import constants
21from polars_analysis.analysis.cut_thresholds import CutThresholds
22from polars_analysis.plotting import helper
23from polars_analysis.plotting.helper import Metadata
25# Instantiate logger
26log = logging.getLogger(__name__)
28"""
29Functions for plotting pedestal run data, which has already been processed.
30"""
33def plot_raw(
34 info: Metadata,
35 samples: np.ndarray,
36 plot_dir: pathlib.Path,
37) -> None:
38 """
39 Plot the raw pedestal samples for a given channel and gain.
41 :param run_number: The run number.
42 :type run_number: int
43 :param channel: The channel number.
44 :type channel: int
45 :param gain: The gain setting.
46 :type gain: Literal["hi", "lo"]
47 :param samples: The raw pedestal samples.
48 :type samples: np.ndarray
49 :param plot_dir: The directory to save the plot.
50 :type plot_dir: pathlib.Path
51 """
52 plt.figure(figsize=(10, 10))
53 plt.plot(
54 samples,
55 "k.",
56 label=rf"""Samples
57$\mu$ = {np.mean(samples):.02f}, $\sigma$ = {np.std(samples):.02f}""",
58 )
59 plt.legend(loc=1)
60 plt.title(
61 helper.plot_summary_string(
62 name="Pedestal",
63 board_id=info.board_id,
64 info=info,
65 )
66 )
67 plt.xlabel("Sample #")
68 plt.ylabel("ADC Counts")
69 plt.grid()
70 plt.tight_layout()
71 if info.channels is None:
72 info.channels = "999"
73 plt.savefig(Path(plot_dir) / f"rawPed_meas0_board_{info.board_id}_channel{info.channels.zfill(3)}_{info.gain}.png")
74 plt.cla()
75 plt.clf()
76 plt.close()
79@deprecated("outlier test is no longer used")
80def outlier_test(
81 run_number: int,
82 channel: int,
83 gain: Literal["hi", "lo"],
84 samples: np.ndarray,
85 plot_dir: pathlib.Path,
86 board_id: Optional[str] = None,
87 attenuation: Optional[str] = None,
88 pas_mode: Optional[Any] = None,
89 info: Optional[Metadata] = None,
90) -> None: # pragma: no cover
91 """
92 Perform an outlier test by creating a histogram of the deviation from mean,
93 divided by standard deviation.
95 :param run_number: The run number.
96 :type run_number: int
97 :param channel: The channel number.
98 :type channel: int
99 :param gain: The gain setting.
100 :type gain: Literal["hi", "lo"]
101 :param samples: The raw pedestal samples.
102 :type samples: np.ndarray
103 :param plot_dir: The directory to save the plot.
104 :type plot_dir: pathlib.Path
105 """
106 mean = np.mean(samples)
107 std = np.std(samples)
108 deviations = (samples - mean) / std # Calculate the deviation from the mean in units of standard deviation
110 plt.figure(figsize=(10, 10))
111 plt.hist(deviations, bins=100, log=True, color="blue", alpha=0.7) # Logarithmic y-axis
112 plt.title(
113 helper.plot_summary_string(
114 board_id=board_id,
115 run_numbers=run_number,
116 channels=helper.list_to_text_ranges(channel),
117 gain=gain,
118 name="Outlier Test",
119 attenuation=attenuation,
120 pas_mode=pas_mode,
121 info=info,
122 )
123 )
124 plt.xlabel("Deviation from Mean (σ)")
125 plt.ylabel("log(N samples)")
126 plt.grid(True)
127 plt.tight_layout()
128 plt.savefig(Path(plot_dir) / f"outlier_test_channel{channel:03}_{gain}.png")
129 plt.cla()
130 plt.clf()
131 plt.close()
134def plot_hist(
135 info: Metadata,
136 samples: np.ndarray,
137 plot_dir: pathlib.Path,
138) -> None:
139 """
140 Plot the histogram of the pedestal samples for a given channel and gain.
142 :param run_number: The run number.
143 :type run_number: int
144 :param channel: The channel number.
145 :type channel: int
146 :param gain: The gain setting.
147 :type gain: Literal["hi", "lo"]
148 :param samples: The pedestal samples.
149 :type samples: np.ndarray
150 :param plot_dir: The directory to save the plot.
151 :type plot_dir: pathlib.Path
152 """
153 min_samples = min(samples)
154 max_samples = max(samples)
155 stdev = np.std(samples)
156 skew = scipy.stats.skew(samples)
158 _, (ax1, ax2) = plt.subplots(2, gridspec_kw={"hspace": 0, "height_ratios": [3, 1]}, sharex=True)
159 ax1.label_outer()
160 ax2.label_outer()
162 bin_width = 1
163 bins = np.arange(min_samples, max_samples + bin_width, bin_width)
164 if len(bins) > 1:
165 bin_content, _, h = ax1.hist(samples, bins=bins, align="mid", edgecolor="black")
166 _, dip_pval = diptest(bin_content) # type: ignore
167 h[0].set_label(f"RMS={stdev:.01f}, γ={skew:.02f}, dip={dip_pval:.02f}")
168 else:
169 bin_content, _, h = ax1.hist(samples, align="mid", edgecolor="black", label=f"RMS={stdev:.01f}")
171 fit_pars = helper.calc_gaussian(samples, bins)
172 bin_centers = 0.5 * (bins[1:] + bins[:-1])
173 gauss_fit = helper.gauss(bin_centers, mu=fit_pars[0], sigma=fit_pars[2], N=fit_pars[4])
175 ax1.plot(
176 bin_centers,
177 gauss_fit,
178 color="r",
179 label=rf"""$\mu$ = {fit_pars[0]:.03f} $\pm$ {fit_pars[1]:.03f},
180$\sigma$ = {fit_pars[2]:.03f} $\pm$ {fit_pars[3]:.03f}""",
181 )
183 ax1.legend(loc="upper right")
184 ax1.set_xlabel("ADC Counts")
185 left, right = ax1.get_xlim()
186 ax1.set_xlim(left, right + (left - right) * -0.1)
187 bottom, top = ax1.get_ylim()
188 ax1.set_ylim(bottom, top + (bottom - top) * -0.1)
189 ax1.set_ylabel("Entries")
190 ax1.set_title(
191 helper.plot_summary_string(
192 name="Pedestal",
193 board_id=info.board_id,
194 info=info,
195 )
196 )
197 ax1.grid(True)
199 dof = len(bin_centers) - 3
200 a = 0.32 # 1 sigma
201 err_up = gamma.ppf(1 - a / 2, gauss_fit + 1, scale=1) - gauss_fit # approx up side of Poisson interval
202 err_dw = gauss_fit - gamma.ppf(a / 2, gauss_fit, scale=1) # approx down side of Poisson interval
203 residuals = bin_content - gauss_fit
204 err = np.where(residuals > 0, err_up, err_dw)
205 residual_sigma = residuals / err
206 chi2_dof = np.sum((bin_content - gauss_fit) ** 2 / err**2) / dof
208 ax2.plot(bin_centers, residual_sigma, color="r", label="Fit residual σ", marker=".", linestyle="None")
209 ax2.text(0.8, 0.8, f"Χ²/dof={chi2_dof:.02f}", transform=ax2.transAxes)
210 ax2.set_ylabel("σ")
211 ax2.set_xlabel("ADC Counts")
212 ax2.grid(True)
214 if info.channels is None:
215 info.channels = "999"
216 plt.savefig(Path(plot_dir) / f"channel{info.channels.zfill(3)}_{info.gain}_board_{info.board_id}_pedestal_hist.png")
217 plt.cla()
218 plt.clf()
219 plt.close()
222def plot_autocorrelation(
223 run_number: int,
224 channel: int,
225 gain: Literal["hi", "lo"],
226 autocorrelation: np.ndarray,
227 plot_dir: pathlib.Path,
228 board_id: Optional[str] = None,
229 attenuation: Optional[str] = None,
230 pas_mode: Optional[Any] = None,
231 info: Optional[Metadata] = None,
232) -> None:
233 """
234 Plot the autocorrelation of the pedestal samples for a given channel and gain.
236 :param run_number: The run number.
237 :type run_number: int
238 :param channel: The channel number.
239 :type channel: int
240 :param gain: The gain setting.
241 :type gain: Literal["hi", "lo"]
242 :param autocorrelation: The autocorrelation of the pedestal samples.
243 :type autocorrelation: np.ndarray
244 :param plot_dir: The directory to save the plot.
245 :type plot_dir: pathlib.Path
246 """
247 plt.plot(autocorrelation[0:50], "k.-", label="Autocorrelation")
248 plt.axhline(y=0, color="r", ls="--")
249 plt.legend(loc="upper right")
250 plt.title(
251 helper.plot_summary_string(
252 board_id=board_id,
253 run_numbers=run_number,
254 channels=helper.list_to_text_ranges(channel),
255 gain=gain,
256 name="Pedestal",
257 attenuation=attenuation,
258 pas_mode=pas_mode,
259 info=info,
260 )
261 )
262 plt.xlabel("Lag")
263 plt.ylabel("Autocorrelation")
264 plt.grid()
265 plt.tight_layout()
266 plt.savefig(Path(plot_dir) / f"autocorr_board_{board_id}_channel{channel:03}_{gain}.png")
267 plt.cla()
268 plt.clf()
269 plt.close()
272def plot_baseline_means_rms(
273 derived_df: pl.DataFrame,
274 plot_dir: pathlib.Path,
275 skip_channels_hi: Optional[List[int]] = None,
276 skip_channels_lo: Optional[List[int]] = None,
277 draw_bounds: Optional[bool] = False,
278 run_number: Optional[int] = None,
279 board_id: Optional[str] = None,
280 attenuation: Optional[str] = None,
281 pas_mode: Optional[Any] = None,
282 info: Optional[Metadata] = None,
283) -> None:
284 """
285 Plot the mean, RMS, and maxmin of the pedestal values for a given measurement and channels for high and low gain.
287 :param derived_df: DataFrame containing mean, std, and maxmin columns for hi and lo gains
288 :type derived_df: pl.DataFrame
289 :param plot_dir: The directory to save the plot.
290 :type plot_dir: pathlib.Path
291 """
292 n_channels = derived_df["channel"].n_unique()
294 color_dict = {"lo": "b", "hi": "r"}
295 title_dict = {"lo": "LG", "hi": "HG"}
297 fig, ax = plt.subplots()
298 plt.xticks(np.arange(0, n_channels, 4), rotation=70)
299 ax.xaxis.set_tick_params(pad=0.1)
300 fig2, ax2 = plt.subplots(1)
301 plt.xticks(np.arange(0, n_channels, 4), rotation=70)
302 ax2.xaxis.set_tick_params(pad=0.1)
303 fig3, ax3 = plt.subplots(1)
304 plt.xticks(np.arange(0, n_channels, 4), rotation=70)
305 ax3.xaxis.set_tick_params(pad=0.1)
306 max_means = 0
307 max_rms = 0
308 max_maxmins = 0
309 gains: List[Literal["hi", "lo"]] = ["hi", "lo"]
310 for gain in gains:
311 gain_df = derived_df.filter(pl.col("gain") == gain)
312 gain_df = gain_df.sort(by="measurement").unique("channel", keep="first").sort(by="channel")
313 if gain == "hi" and skip_channels_hi is not None:
314 gain_df = gain_df.filter(~pl.col("channel").is_in(skip_channels_hi))
315 elif gain == "lo" and skip_channels_lo is not None:
316 gain_df = gain_df.filter(~pl.col("channel").is_in(skip_channels_lo))
317 means = gain_df["mean"].to_numpy()
318 stds = gain_df["std"].to_numpy()
319 maxmins = gain_df["maxmin"].to_numpy()
320 channels = gain_df["channel"].to_numpy()
321 color = color_dict[gain]
322 title = title_dict[gain]
324 ax.grid(visible=True, zorder=0)
325 mean = np.mean(means)
326 std = np.std(means) / np.sqrt(len(means))
327 ax.bar(channels, means, fill=False, ec=color, label=f"{title} mean = {mean:.2f}±{std:.2f}", zorder=3)
328 max_means = max(max_means, max(means))
329 ax.margins(x=0)
331 ax2.grid(visible=True, zorder=0)
332 mean = np.mean(stds)
333 std = np.std(stds) / np.sqrt(len(stds))
334 ax2.bar(channels, stds, fill=False, ec=color, label=f"{title} mean = {mean:.2f}±{std:.2f}", zorder=3)
335 max_rms = max(max_rms, max(stds))
336 ax2.margins(x=0)
338 ax3.grid(visible=True, zorder=0)
339 mean = np.mean(maxmins)
340 std = np.std(maxmins) / np.sqrt(len(maxmins))
341 ax3.bar(channels, maxmins, fill=False, ec=color, label=f"{title} mean = {mean:.2f}±{std:.2f}", zorder=3)
342 max_maxmins = max(max_maxmins, max(maxmins))
343 ax3.margins(x=0)
344 # Draw acceptance cuts
345 if draw_bounds:
346 cuts = CutThresholds()
347 cuts.draw_on(ax, "mean", gain=gain)
348 cuts.draw_on(ax2, "std", gain=gain)
349 cuts.draw_on(ax3, "maxmin", gain=gain)
351 all_channels = derived_df["channel"].unique().sort().to_list()
352 ax.set_title(
353 helper.plot_summary_string(
354 board_id=board_id,
355 run_numbers=run_number,
356 channels=helper.list_to_text_ranges(all_channels),
357 name="Mean Pedestal Value",
358 attenuation=attenuation,
359 pas_mode=pas_mode,
360 info=info,
361 )
362 )
363 ax.set_xlabel("Channel", loc="right")
364 ax.set_ylabel("ADC Counts")
365 if max_means > 10 * np.median(means):
366 ax.set_yscale("log")
367 else:
368 ax.set_ylim(0, 1.33 * max_means)
369 ax.legend()
370 fig.tight_layout()
371 fig.savefig(f"{plot_dir}/mu_board_{board_id}_summary.png")
372 fig.clf()
374 ax2.set_title(
375 helper.plot_summary_string(
376 board_id=board_id,
377 run_numbers=run_number,
378 channels=helper.list_to_text_ranges(all_channels),
379 name="Pedestal RMS",
380 attenuation=attenuation,
381 pas_mode=pas_mode,
382 info=info,
383 )
384 )
385 ax2.set_xlabel("Channel", loc="right")
386 ax2.set_ylabel("ADC Counts")
387 if max_rms > 10 * np.median(stds):
388 ax2.set_yscale("log")
389 else:
390 ax2.set_ylim(0, 1.33 * max_rms)
391 ax2.legend()
392 fig2.tight_layout()
393 fig2.savefig(f"{plot_dir}/rms_board_{board_id}_summary.png")
394 fig2.clf()
396 ax3.set_title(
397 helper.plot_summary_string(
398 board_id=board_id,
399 run_numbers=run_number,
400 channels=helper.list_to_text_ranges(all_channels),
401 name="Pedestal Max-Min",
402 attenuation=attenuation,
403 pas_mode=pas_mode,
404 info=info,
405 )
406 )
407 ax3.set_xlabel("Channel", loc="right")
408 ax3.set_ylabel("ADC Counts")
409 if max_maxmins > 10 * np.median(maxmins):
410 ax3.set_yscale("log")
411 else:
412 ax3.set_ylim(0, 1.33 * max_maxmins)
413 ax3.legend()
414 fig3.tight_layout()
415 fig3.savefig(f"{plot_dir}/maxmin_board_{board_id}_summary.png")
416 fig3.clf()
418 plt.cla()
419 plt.clf()
420 plt.close()
423def plot_coherent_noise(
424 row: Dict[str, Any],
425 plot_dir: pathlib.Path,
426 meas_type: str = "pedestal",
427 board_id: Optional[str] = None,
428 attenuation: Optional[str] = None,
429 pas_mode: Optional[Any] = None,
430 use_log_scale: bool = False,
431) -> Union[tuple[List[str], Any], tuple[None, None]]:
432 """
433 Plot the coherent noise for a given measurement type and gain, given a DataFrame of coherent noise results.
435 :param row: A dictionary from pl.DataFrame.iter_rows(named=True) containing the coherent noise results.
436 Must contain the following columns:
438 * min_channel: The minimum channel number.
439 * n_channels: The number of channels.
440 * gain: The gain setting.
441 * data_sum_hist: Histogram of the per samples sum of the pedestal values across channels.
442 * data_sum_bins: Bin edges for data_sum_hist
443 * tot_noise: The total noise.
444 * ch_noise: The channel noise.
445 * avg_noise: The average noise.
446 * coh_noise: The coherent noise.
447 * pct_coh: The percentage of coherent noise.
448 * d_tot_noise: The total noise uncertainty.
449 * d_ch_noise: The channel noise uncertainty.
450 * d_avg: The average noise uncertainty.
451 * d_coh: The coherent noise uncertainty.
452 * d_pct: The percentage of coherent noise uncertainty.
453 :type row: dict[str, Any]
454 :param run_number: The run number.
455 :type run_number: int
456 :param plot_dir: The directory to save the plot.
457 :type plot_dir: pathlib.Path
458 :param meas_type: The measurement type, defaults to "pedestal".
459 :type meas_type: str, optional
460 """
461 run_number: int = row["run_number"]
462 if not board_id:
463 board_id = row["board_id"] if "board_id" in row else None
464 min_channel: int = row["min_channel"]
465 n_channels: int = row["n_channels"]
466 gain: Literal["hi", "lo"] = row["gain"]
467 data_sum_hist: np.ndarray = np.array(row["data_sum_hist"])
468 bins: np.ndarray = np.array(row["data_sum_bins"])
469 channel_list = row["channel_list"]
471 centers = 0.5 * (bins[1:] + bins[:-1])
472 mu = helper.hist_mean(centers, data_sum_hist)
473 skew = helper.hist_moment(centers, data_sum_hist, 3)
475 _, (ax1, ax2) = plt.subplots(2, gridspec_kw={"hspace": 0, "height_ratios": [3, 1]}, sharex=True)
476 ax1.label_outer()
477 ax2.label_outer()
479 _, dip_pval = diptest(data_sum_hist) # type: ignore
480 _ = ax1.bar(
481 x=bins[:-1],
482 height=data_sum_hist,
483 width=bins[1:] - bins[:-1],
484 align="edge",
485 color="steelblue",
486 fill=True,
487 edgecolor="black",
488 label=rf"""$\mu$ = {mu:.02f}, $\sigma$ = {row["tot_noise"]:.02f},
489γ={skew:.02f}, dip={dip_pval:.02f}""",
490 )
492 fit_pars = helper.calc_gaussian_from_bins(data_sum_hist, bins)
493 fit_mu, fit_sigma, fit_N = fit_pars[0], fit_pars[2], fit_pars[4]
494 e_fit_mu, e_fit_sigma = fit_pars[1], fit_pars[3]
495 bin_centers = 0.5 * (bins[1:] + bins[:-1])
496 gauss_fit = helper.gauss(bin_centers, fit_mu, fit_sigma, fit_N)
497 ax1.plot(
498 bin_centers,
499 gauss_fit,
500 color="r",
501 label=rf"""Gaussian Fit
502$\mu$ = {fit_mu:.02f}±{e_fit_mu:0.2f}, $\sigma$ = {fit_sigma:.02f}±{e_fit_sigma:.02f}""",
503 )
504 ax1.set_ylim(0, max(data_sum_hist) * 1.4)
505 ax1.set_xlabel("ADC Counts")
506 ax1.set_ylabel("Entries")
507 ax1.set_title(
508 helper.plot_summary_string(
509 board_id=board_id,
510 run_numbers=run_number,
511 channels=helper.list_to_text_ranges(channel_list),
512 gain=gain,
513 name=f"{meas_type} coherent noise",
514 attenuation=attenuation,
515 pas_mode=pas_mode,
516 )
517 )
518 ax1.grid()
519 ax1.text(
520 0.995,
521 0.84,
522 rf"""$\sqrt{{\sigma_i^2}}$ = {row["ch_noise"]:.02f} $\pm$ {row["d_ch_noise"]:.02f}
523 Avg. noise/ch = {row["avg_noise"]:.02f} $\pm$ {row["d_avg"]:.02f}
524 Coh. noise/ch = {row["coh_noise"]:.02f} $\pm$ {row["d_coh"]:.02f}
525 [%] Coh. noise = {row["pct_coh"]:.02f} $\pm$ {row["d_pct"]:.02f}""",
526 horizontalalignment="right",
527 verticalalignment="center",
528 transform=ax1.transAxes,
529 )
530 log_filename: Literal["", "_log"] = ""
531 if use_log_scale:
532 ax1.set_yscale("log")
533 ax1.set_ylim(bottom=0.5)
534 log_filename = "_log"
535 ax1.legend(loc="upper left", handlelength=1.0)
537 dof = len(bin_centers) - 3
538 a = 0.32 # 1 sigma
539 err_up = gamma.ppf(1 - a / 2, gauss_fit + 1, scale=1) - gauss_fit # approx up side of Poisson interval
540 err_dw = gauss_fit - gamma.ppf(a / 2, gauss_fit, scale=1) # approx down side of Poisson interval
541 residuals = data_sum_hist - gauss_fit
542 err = np.where(residuals > 0, err_up, err_dw)
543 residual_sigma = residuals / err
544 chi2_dof = np.sum((data_sum_hist - gauss_fit) ** 2 / err**2) / dof
546 ax2.plot(bin_centers, residual_sigma, color="r", label="Fit residual σ", marker=".", linestyle="None")
547 ax2.text(0.8, 0.8, f"Χ²/dof={chi2_dof:.02f}", transform=ax2.transAxes)
548 ax2.set_ylabel("σ")
549 ax2.set_xlabel("ADC Counts")
550 ax2.grid(True)
552 if n_channels % 128 == 0:
553 plt.savefig(plot_dir / f"coherence_all_{gain}_board_{board_id}_pedestal_hist{log_filename}.png")
554 else:
555 plt.savefig(
556 plot_dir / f"{gain}_board_{board_id}_cnoise_ch_{min_channel}_{min_channel + n_channels}{log_filename}.png"
557 )
558 plt.cla()
559 plt.clf()
560 plt.close()
562 if board_id is not None:
563 return board_id.split("_"), row["pct_coh"]
564 else:
565 return None, None
568def plot_coherent_noise_matrix(coh_matrix: Dict[str, Dict[str, Any]], plot_dir: pathlib.Path) -> None:
569 for gain in coh_matrix.keys():
570 # Parse keys to extract board combinations
571 unsorted_individual_boards: Set[str] = set()
572 board_combinations = {}
574 for key, value in coh_matrix[gain].items():
575 boards_in_key = key.split("_")
576 board_combinations[tuple(sorted(boards_in_key))] = value
577 unsorted_individual_boards.update(boards_in_key)
579 individual_boards: List[str] = sorted(unsorted_individual_boards)
580 n_boards = len(individual_boards)
582 if n_boards == 0:
583 log.error("Error: No boards found!")
584 return
586 # Find the combined noise value
587 all_boards_key = tuple(individual_boards)
588 combined_noise = board_combinations.get(all_boards_key, np.nan)
590 matrix = np.zeros((n_boards, n_boards))
592 for i, board_x in enumerate(individual_boards):
593 for j, board_y in enumerate(individual_boards):
594 if i == j:
595 # Diagonal: individual board coherent noise
596 single_board_key = (board_x,)
597 if single_board_key in board_combinations:
598 matrix[i, j] = board_combinations[single_board_key]
599 else:
600 log.warning(f"Warning: Individual board key {single_board_key} not found")
601 matrix[i, j] = np.nan
602 else:
603 # Off-diagonal: pairwise board coherent noise
604 pairwise_key = tuple(sorted([board_x, board_y]))
605 if pairwise_key in board_combinations:
606 matrix[i, j] = board_combinations[pairwise_key]
607 else:
608 log.warning(f"Warning: Pairwise key {pairwise_key} not found")
609 matrix[i, j] = np.nan
611 vmin = 0.1 # Small positive value
612 vmax = max(np.max(matrix).astype(float), 6)
613 norm = mcolors.TwoSlopeNorm(vmin=vmin, vcenter=5, vmax=vmax)
614 fig, ax = plt.subplots(figsize=(max(6, n_boards * 1.5), max(6, n_boards * 1.5)))
615 im = ax.imshow(matrix, cmap="Reds", norm=norm)
617 ax.set_xticks(np.arange(n_boards))
618 ax.set_yticks(np.arange(n_boards))
619 ax.set_xticklabels(individual_boards, rotation=45, fontsize=8)
620 ax.set_yticklabels(individual_boards, fontsize=8)
622 # Add text annotations
623 for i, j in product(range(n_boards), range(n_boards)):
624 if np.isnan(matrix[i, j]):
625 val_to_write = "N/A"
626 color = "gray"
627 else:
628 val_to_write = f"{matrix[i, j]:.1f}"
629 color = "k"
631 ax.text(j, i, val_to_write, ha="center", fontsize=16, va="center", color=color)
633 cbar = fig.colorbar(im, ax=ax)
634 cbar.set_label("Coherent Noise [%]", rotation=270, labelpad=15)
635 gain_str = "HI" if gain == "hi" else "LO"
636 if not np.isnan(combined_noise):
637 title = f"{gain_str} Gain Coherent Noise Per Channel Matrix [%]\n \
638 Combined value for all boards: {combined_noise:.1f}%"
639 else:
640 title = f"{gain_str} Gain Coherent Noise Per Channel Matrix [%]"
642 ax.set_title(title)
643 ax.set_xlabel("Board ID")
644 ax.set_ylabel("Board ID")
645 fig.tight_layout()
646 plt.savefig(plot_dir / f"coherent_noise_matrix_{gain}.png", dpi=150, bbox_inches="tight")
647 plt.close()
650def plot_correlation_matrix(
651 matrix: np.ndarray,
652 gain: Literal["hi", "lo"],
653 min_channel: int,
654 nchannels: int,
655 plot_dir: pathlib.Path,
656 run_number: Optional[int] = None,
657 board_id: Optional[str] = None,
658 attenuation: Optional[str] = None,
659 pas_mode: Optional[Any] = None,
660 info: Optional[Metadata] = None,
661 plot_numbers: bool = True,
662) -> None:
663 """
664 Plot the correlation matrix for a given gain, minimum channel, and number of channels.
666 :param matrix: The correlation matrix.
667 :type matrix: np.ndarray
668 :param gain: The gain setting.
669 :type gain: Literal["hi", "lo"]
670 :param min_channel: The minimum channel number.
671 :type min_channel: int
672 :param nchannels: The number of channels.
673 :type nchannels: int
674 :param plot_dir: The directory to save the plot.
675 :type plot_dir: pathlib.Path
676 """
677 log.debug("Plotting correlation matrix")
679 plot_all = nchannels % 128 == 0
680 max_channel = min_channel + nchannels
681 if min_channel > matrix.shape[0]:
682 # The given interval does not overlap with available channels
683 log.error(f"min_channel {min_channel} > matrix shape {matrix.shape[0]}!")
684 return
686 submatrix = matrix[min_channel:max_channel, min_channel:max_channel]
688 fig, ax = plt.subplots(figsize=(nchannels / 4, nchannels / 4))
689 _ = ax.imshow(submatrix, cmap="RdBu", vmin=-0.3, vmax=0.3)
690 ax.set_xticks(np.arange(0, nchannels, max(1, round(nchannels / 32))))
691 ax.set_yticks(np.arange(0, nchannels, max(1, round(nchannels / 32))))
692 ax.set_xticklabels(
693 np.arange(min_channel, min_channel + nchannels, max(1, round(nchannels / 32))),
694 rotation=90,
695 fontsize=7,
696 )
697 ax.set_yticklabels(np.arange(min_channel, min_channel + nchannels, max(1, round(nchannels / 32))), fontsize=7)
699 # Use grid lines at every 128 entries
700 for pos in np.arange(128, nchannels, 128):
701 ax.axvline(pos - 0.5, color="black", linewidth=3, alpha=0.7)
702 ax.axhline(pos - 0.5, color="black", linewidth=3, alpha=0.7)
704 if plot_numbers:
705 # for i, j in product(range(min_channel, min(128, max_channel)), range(min_channel, min(128, max_channel))):
706 for i, j in product(range(min_channel, max_channel), range(min_channel, max_channel)):
707 val_to_write = " " if np.isnan(matrix[i, j]) else f"{100 * matrix[i, j]:.0f}"
708 color = "w" if ((i == j) and (val_to_write != "0")) else "k"
709 _ = ax.text(
710 j - min_channel, i - min_channel, val_to_write, ha="center", fontsize=6, va="center", color=color
711 )
713 default_fontsize = plt.rcParams["font.size"]
714 ax.set_title(
715 helper.plot_summary_string(
716 board_id=board_id,
717 run_numbers=run_number,
718 channels=f"{min_channel}-{max_channel - 1}",
719 gain=gain,
720 name="pearson correlation [%]",
721 attenuation=attenuation,
722 pas_mode=pas_mode,
723 info=info,
724 ),
725 fontsize=default_fontsize / (1.2 if nchannels < 64 else 1),
726 )
727 fig.tight_layout()
728 if plot_all:
729 plt.savefig(plot_dir / f"{gain}_board_{board_id}_corr.png")
730 else:
731 plt.savefig(plot_dir / f"{gain}_board_{board_id}_corr_ch_{min_channel}_{max_channel}.png")
732 plt.cla()
733 plt.clf()
734 plt.close()
737def plot_fft(
738 channel: int,
739 gain: Literal["hi", "lo"],
740 freq: List[float],
741 psd: List[float],
742 peaks: List[int],
743 plot_dir: pathlib.Path,
744 run_number: Optional[int] = None,
745 board_id: Optional[str] = None,
746 attenuation: Optional[str] = None,
747 pas_mode: Optional[Any] = None,
748 info: Optional[Metadata] = None,
749) -> None:
750 """
751 Plot the FFT of the pedestal samples for a given channel and gain.
753 :param channel: The channel number.
754 :type channel: int
755 :param gain: The gain setting.
756 :type gain: Literal["hi", "lo"]
757 :param freq: The frequency array.
758 :type freq: np.ndarray
759 :param psd: The power spectral density array.
760 :type psd: np.ndarray
761 :param plot_dir: The directory to save the plot.
762 :type plot_dir: pathlib.Path
763 """
765 # # Convert to dBFS
766 psd_dbfs = 10 * np.log10(np.array(psd) / (2**constants.ADC_BITS) ** 2)
768 plt.plot(freq, psd_dbfs, color="k")
769 plt.plot([freq[peak] for peak in peaks], [psd_dbfs[peak] for peak in peaks], "xr")
770 for peak in peaks:
771 peak_psd_dbfs = float(psd_dbfs[peak])
772 plt.text(freq[peak] - 0.5, peak_psd_dbfs + 0.5, f"{freq[peak]:.02f}", rotation=45, color="k")
773 plt.ylim(top=max(psd_dbfs) + 2)
774 plt.xlabel("Frequency [MHz]")
775 plt.ylabel("PSD [dB Full Scale]")
776 plt.title(
777 helper.plot_summary_string(
778 board_id=board_id,
779 run_numbers=run_number,
780 channels=helper.list_to_text_ranges(channel),
781 gain=gain,
782 name="Power Spectral Density",
783 attenuation=attenuation,
784 pas_mode=pas_mode,
785 info=info,
786 )
787 )
788 plt.grid()
789 plt.tight_layout()
790 _ = plt.subplot(111)
791 plt.savefig(plot_dir / f"pedestal_FFT_board_{board_id}_channel{channel:03}_{gain}.png")
792 plt.cla()
793 plt.clf()
794 plt.close()
797def plot_fft2d(
798 freq: List[float],
799 ffts: Union[pl.Series, np.ndarray],
800 plot_dir: pathlib.Path,
801 channels: Union[pl.Series, np.ndarray, List[int]],
802 gain: Literal["hi", "lo"],
803 run_number: Optional[int] = None,
804 board_id: Optional[str] = None,
805 pas_mode: Optional[Any] = None,
806 unit: str = "MHz",
807 extra_filename: str = "",
808) -> None:
809 """
810 Plot the 2D FFT of the pedestal samples for a given set of channels and gain.
812 :param freq: The frequency of the fft.
813 :type freq: List[float]
814 :param ffts: The 2D list of ffts. ffts[i] is the fft for channel i
815 :type ffts: pl.Series
816 :param plot_dir: The directory to save the plot.
817 :type plot_dir: pathlib.Path
818 :param channels: The channels that should be plotted
819 :type channels: pl.Series
820 :param gain: The gain setting.
821 :type gain: Literal["hi", "lo"]
822 :param run_number: The run number.
823 :type run_number: int
824 :param board_id: The board id.
825 :type board_id: Optional[str]
826 :param pas_mode: The ALFE mode or HPS mode.
827 :type pas_mode: Optional[str]
828 """
829 plt.figure(figsize=(12, 4))
831 psd_dbfs = 10 * np.log10(np.stack(ffts) / (2**constants.ADC_BITS) ** 2) # type: ignore
833 arr = np.zeros((128, len(freq)))
835 for index, channel in enumerate(channels):
836 arr[channel] = psd_dbfs[index]
838 plt.pcolormesh(freq, list(range(128)), arr)
840 if isinstance(channels, pl.Series):
841 channels = channels.to_numpy()
842 plt.title(
843 helper.plot_summary_string(
844 board_id=board_id,
845 run_numbers=run_number,
846 channels=helper.list_to_text_ranges(channels),
847 gain=gain,
848 name="Power Spectral Density",
849 pas_mode=pas_mode,
850 )
851 )
853 plt.ylabel("Channel")
854 plt.xlabel(f"Frequency [{unit}]")
855 plt.gca().xaxis.set_major_locator(ticker.MultipleLocator(1))
857 cbar = plt.colorbar()
858 cbar.ax.get_yaxis().labelpad = 15
859 cbar.set_label("PSD [dB Full Scale]")
861 plt.tight_layout()
863 plt.savefig(plot_dir / f"pedestal_FFT_2D_board_{board_id}_{gain}{extra_filename}.png")
864 plt.cla()
865 plt.clf()
866 plt.close()
869def plot_coherence(
870 channel1: int,
871 channel2: int,
872 gain: Literal["hi", "lo"],
873 freq: np.ndarray,
874 coh: np.ndarray,
875 plot_dir: pathlib.Path,
876 run_number: Optional[int] = None,
877 board_id: Optional[str] = None,
878 attenuation: Optional[str] = None,
879 pas_mode: Optional[Any] = None,
880 info: Optional[Metadata] = None,
881) -> None:
882 """
883 Plot the coherence between two channels.
885 :param channel: The channel number.
886 :type channel: int
887 :param gain: The gain setting.
888 :type gain: Literal["hi", "lo"]
889 :param freq: The frequency array.
890 :type freq: np.ndarray
891 :param coh: The coherence array.
892 :type coh: np.ndarray
893 :param plot_dir: The directory to save the plot.
894 :type plot_dir: pathlib.Path
895 """
896 plt.semilogy(
897 freq,
898 coh,
899 color="k",
900 )
901 # plt.legend(loc=1)
902 plt.xlabel("Frequency [MHz]")
903 plt.ylabel("Coherence")
904 plt.title(
905 helper.plot_summary_string(
906 board_id=board_id,
907 run_numbers=run_number,
908 channels=f"{channel1} & {channel2}",
909 gain=gain,
910 name="Pedestal Coherence",
911 attenuation=attenuation,
912 pas_mode=pas_mode,
913 info=info,
914 )
915 )
916 yLow, _ = plt.gca().get_ylim()
917 plt.gca().set_ylim(yLow, 2)
918 plt.grid()
919 plt.tight_layout()
920 _ = plt.subplot(111)
921 plt.savefig(plot_dir / f"pedestal_coherence_channel{channel1}-{channel2}_{gain}.png")
922 plt.cla()
923 plt.clf()
924 plt.close()
926 log.info(f"Wrote {plot_dir}/pedestal_coherence_channel{channel1}-{channel2}_{gain}.png")
929def plot_chi2(
930 derived_df: pl.DataFrame,
931 gain: Literal["hi", "lo"],
932 plot_dir: pathlib.Path,
933 run_number: Optional[int] = None,
934 board_id: Optional[str] = None,
935 pas_mode: Optional[Any] = None,
936 info: Optional[Metadata] = None,
937) -> None:
938 """
939 Plot χ²/dof vs. channel for the specified gain using the derived DataFrame.
941 :param derived_df: Polars DataFrame containing "channel", "chi2_dof", and "gain" columns.
942 :param gain: Gain setting ("hi" or "lo").
943 :param plot_dir: Directory where the plot will be saved.
944 :param run_number: Optional run number metadata.
945 :param board_id: Optional board id metadata.
946 :param pas_mode: Optional ALFE or HPS mode metadata.
947 :param info: Optional additional metadata.
948 """
949 # Filter for the specified gain and sort by channel
950 df_gain = derived_df.filter(pl.col("gain") == gain).sort("channel")
951 channels = np.array(df_gain["channel"].to_numpy())
952 chi2_values = np.array(df_gain["chi2_dof"].to_numpy())
954 color = "red" if gain == "hi" else "blue"
956 fig, ax = plt.subplots(figsize=(10, 6))
957 ax.plot(channels, chi2_values, marker="o", linestyle="-", color=color, label=f"{gain.upper()} Gain")
958 ax.set_xlabel("Channel")
959 ax.set_ylabel("χ²/dof")
960 ax.set_yscale("log")
961 ax.set_title(
962 helper.plot_summary_string(
963 board_id=board_id,
964 run_numbers=run_number,
965 channels="All",
966 gain=gain,
967 name="χ²/dof vs. Channel",
968 pas_mode=pas_mode,
969 info=info,
970 )
971 )
972 ax.grid(True)
973 ax.legend()
975 fig.tight_layout()
976 plot_file = plot_dir / f"{gain}_board_{board_id}_chi2.png"
977 plt.savefig(plot_file)
978 plt.cla()
979 plt.clf()
980 plt.close()
983def calc_chi2_threshold(chi2_values, k: float = 6.0, j: float = -0.0005) -> float:
984 ### Calculate a threshold for chi² values using median, stdev and MAD ###
985 chi2_arr = np.asarray(chi2_values)
986 med = np.median(chi2_arr)
987 mad = np.median(np.abs(chi2_arr - med))
988 std = np.std(chi2_arr)
990 threshold = med + k * mad + j * std
991 return threshold.astype(float)
994def plot_chi2_hist(
995 derived_df: pl.DataFrame,
996 gain: Literal["hi", "lo"],
997 plot_dir: pathlib.Path,
998 chi2_threshold: Optional[float] = 0.0,
999 run_number: Optional[int] = None,
1000 board_id: Optional[str] = None,
1001 pas_mode: Optional[Any] = None,
1002 info: Optional[Metadata] = None,
1003) -> None:
1004 """
1005 Plot a histogram of the χ²/dof values for the specified gain.
1006 Adds a vertical line at χ²/dof = 1000 and a text box (occupying the right 25% of the plot)
1007 listing flagged channels (those with χ²/dof >= 1000) sorted in descending order.
1009 :param derived_df: Polars DataFrame containing "chi2_dof", "channel", and "gain" columns.
1010 :param gain: The gain setting ("hi" or "lo").
1011 :param plot_dir: The directory where the plot will be saved.
1012 :param run_number: Optional run number metadata.
1013 :param board_id: Optional board id metadata.
1014 :param pas_mode: Optional ALFE or HPS mode metadata.
1015 :param info: Optional additional metadata.
1016 """
1017 df_gain = derived_df.filter(pl.col("gain") == gain)
1018 chi2_values = np.array(df_gain["chi2_dof"].to_numpy())
1020 # Define threshold and identify flagged channels sorted by highest chi² to lowest.
1021 chi2_threshold = constants.CHI2_THRESHOLD
1022 if chi2_threshold == 0:
1023 chi2_threshold = calc_chi2_threshold(chi2_values)
1024 flagged_df = df_gain.filter(pl.col("chi2_dof") >= chi2_threshold).sort("chi2_dof", descending=True)
1025 flagged_list = flagged_df.select(["channel", "chi2_dof"]).to_dicts()
1027 color = "red" if gain == "hi" else "blue"
1029 fig, ax = plt.subplots(figsize=(10, 6))
1030 _ = ax.hist(chi2_values, bins="auto", edgecolor="black", color=color, alpha=0.7)
1032 if flagged_list:
1033 flagged_text = f"Flagged channels:\n(χ²/dof > {chi2_threshold:.4g})\n"
1034 for entry in flagged_list:
1035 flagged_text += f"Ch {entry['channel']}: {entry['chi2_dof']:.1f}\n"
1037 ax.axvline(x=chi2_threshold, color="black", linestyle="--", linewidth=2)
1038 ax.text(
1039 0.84,
1040 0.98,
1041 flagged_text,
1042 transform=ax.transAxes,
1043 fontsize=10,
1044 verticalalignment="top",
1045 horizontalalignment="left",
1046 bbox=dict(boxstyle="round", facecolor="white", alpha=0.8),
1047 )
1049 ax.set_xscale("log")
1050 ax.set_yscale("log")
1052 ax.set_xlabel("χ²/dof")
1053 ax.set_ylabel("Count")
1054 ax.set_title(
1055 helper.plot_summary_string(
1056 board_id=board_id,
1057 run_numbers=run_number,
1058 channels="All",
1059 gain=gain,
1060 name="χ²/dof Histogram",
1061 pas_mode=pas_mode,
1062 info=info,
1063 )
1064 )
1065 ax.grid(True)
1067 fig.tight_layout()
1068 plt.savefig(plot_dir / f"{gain}_board_{board_id}_chi2_hist.png")
1069 plt.cla()
1070 plt.clf()
1071 plt.close()