From c010a86ff3ce33ebd0b1dba0d3f01bb89188e60e Mon Sep 17 00:00:00 2001 From: Patrick Krause Date: Mon, 5 Jan 2026 15:38:44 -0800 Subject: [PATCH 01/16] Convert to Q frame and inject detector info into G frame --- src/mintanalysis/pmt/i3/daq_to_q_frame.py | 74 ++++++++ src/mintanalysis/pmt/i3/insert_aux.py | 205 ++++++++++++++++++++++ 2 files changed, 279 insertions(+) create mode 100644 src/mintanalysis/pmt/i3/daq_to_q_frame.py create mode 100644 src/mintanalysis/pmt/i3/insert_aux.py diff --git a/src/mintanalysis/pmt/i3/daq_to_q_frame.py b/src/mintanalysis/pmt/i3/daq_to_q_frame.py new file mode 100644 index 0000000..f017c93 --- /dev/null +++ b/src/mintanalysis/pmt/i3/daq_to_q_frame.py @@ -0,0 +1,74 @@ +import argparse +import logging +import os + +from icecube import icetray, pone_unfolding + +from mintanalysis.pmt.ana.utils import setup_logging + + +def main(): + parser = argparse.ArgumentParser(description="Build Q frame from DAQ data.") + parser.add_argument( + "-i", + "--f_i3", + default=None, + help="Path to i3 file (if omitted replaces all occurrences of daq in f_daq with i3)", + ) + parser.add_argument("-d", "--f_daq", help="Path to DAQ file", required=True) + parser.add_argument( + "-n", + "--n_events", + default=0, + type=int, + help="Number of events to process (omit or 0 for all)", + ) + parser.add_argument( + "-o", "--overwrite", action="store_true", help="Override existing output file" + ) + + args = parser.parse_args() + + daq_ext = args.f_daq.split(".")[-1] + f_i3 = ( + args.f_daq.replace("daq", "i3").replace(daq_ext, "i3") if args.f_i3 is None else args.f_i3 + ) + + # Create raw folders if not existing + dir = os.path.dirname(f_i3) + if dir: + os.makedirs(dir, exist_ok=True) + + f_log = f_i3.replace(f_i3.split(".")[-1], "log") + logger = setup_logging(log_file=f_log, level=logging.INFO) + + # Check if I3 file exists and if overwrite flag is set + # if so move it to back-up file (which will be deleted at the end) + if os.path.isfile(f_i3): + if args.overwrite: + os.rename(f_i3, f_i3 + ".bak") + else: + msg = f"I3 file {f_i3} exists and overwrite flag is not specified." + logger.error(msg) + return + + # Generate i3 file + msg = f"Start building i3 file from DAQ file {args.f_daq}" + logger.info(msg) + tray = icetray.I3Tray() + tray.AddModule(pone_unfolding.P1Reader, Input=args.f_daq, OM=1) + tray.AddModule("I3Writer", "writer", Filename=f_i3) + if args.n_events > 0: + tray.Execute(args.n_events) + else: + tray.Execute() + tray.Finish() + msg = f"Finished processing. i3 file {f_i3} generated" + logger.info(msg) + + if os.path.isfile(f_i3 + ".bak"): + os.remove(f_i3 + ".bak") + + +if __name__ == "__main__": + main() diff --git a/src/mintanalysis/pmt/i3/insert_aux.py b/src/mintanalysis/pmt/i3/insert_aux.py new file mode 100644 index 0000000..bae30b0 --- /dev/null +++ b/src/mintanalysis/pmt/i3/insert_aux.py @@ -0,0 +1,205 @@ +import argparse +import logging +import os +import re +from pathlib import Path + +import yaml +from icecube import dataclasses, icetray, p1_dataclasses +from pint import UnitRegistry + +from mintanalysis.pmt.ana.utils import get_physics_object, setup_logging + + +class injectDetectorInfo(icetray.I3Module): + def __init__(self, context): + icetray.I3Module.__init__(self, context) + self.AddParameter("config", "Aux dictionary of Module") + self.AddParameter("scale", "ADC scaling factor in mV (default: 0.25 mV)", 0.25) + self.ctr = True + + def Configure(self): + self.config = self.GetParameter("config") + self.scale = self.GetParameter("scale") + + def Process(self): + self.current_frame = self.PopFrame() + if self.current_frame.Stop == icetray.I3Frame.DAQ and self.ctr: + self.Inject() + self.ctr = False + self.PushFrame(self.current_frame) + + def Inject(self): + # Inject Aux data (geometry frame) + frame = icetray.I3Frame(icetray.I3Frame.Geometry) + if "channels" not in self.config: + return + channels = [icetray.OMKey(0, 1, i) for i in self.config["channels"]] + for k, v in self.config.items(): + if isinstance(v, list): + map = dataclasses.I3MapKeyDouble() + for i in range(len(channels)): + map[channels[i]] = v[i] + if k == "v10_in_V": + frame.Put("PMTVoltages", map) + else: + frame.Put(k, map) + else: + frame.Put(k, map) + + # Inject ChannelCalibrationMap (calibration frame) + cal = p1_dataclasses.P1ChannelCalibrationMap() + for ch in channels: + ccal = p1_dataclasses.P1ChannelCalibration() + ccal.SetBaseline(frame["channel_means"][ch]) + ccal.SetVoltageFactor(self.scale * icetray.I3Units.millivolt) + ccal.SetNoiseLevel(frame["noise_levels"][ch]) + cal[ch] = ccal + + self.PushFrame(frame) + frame = icetray.I3Frame(icetray.I3Frame.Calibration) + frame.Put("P1ChannelCalibration", cal) + self.PushFrame(frame) + + +def load_aux(aux_yaml: Path, key: str, ch_mask: int = 0xFFFF) -> dict: + # get list of channels from channel mask + channels = [] + for i in range(16): + bit_value = 1 << i + if ch_mask & bit_value: + channels.append(i) + + ureg = UnitRegistry() + if not aux_yaml.exists(): + msg = f"Aux file not found: {aux_yaml}" + raise FileNotFoundError(msg) + with open(aux_yaml) as f: + aux = yaml.safe_load(f) + + # convert to physics units + aux = get_physics_object(aux, ureg) + + aux = aux[key] + + ret = {"channels": [], "channel_thresholds": [], "channel_means": [], "noise_levels": []} + + for i in range(16): + if (i in aux) and (i in channels): + ret["channels"].append(i) + ret["channel_thresholds"].append( + aux.get("collector_config", {}).get("ChannelThresholds", [])[i] + ) + ret["channel_means"].append(aux.get("collector_config", {}).get("ChannelMeans", [])[i]) + ret["noise_levels"].append( + (ret["channel_thresholds"][-1] - ret["channel_means"][-1]) / 5 + ) # assumption: threshold is set at 5 sigma + + # make sure values are in desired units + for p in [("v10", "V"), ("vsup", "V"), ("di10", "mA"), ("isup", "mA")]: + if p[0] in aux[i]: + k = f"{p[0]}_in_{p[1]}" + if k not in ret: + ret[k] = [] + ret[k].append(aux[i][p[0]].to(p[1]).magnitude) + + if "runtime" in aux: + ret["runtime_in_s"] = aux["runtime"].to("s").magnitude + + for pth in ["pth_start", "pth_end"]: + if pth in aux: + for v in [ + ("humidity", "%", "percent"), + ("pressure", "Pa"), + ("temperature", "degree_Celsius"), + ]: + if v[0] in aux[pth]: + ret[f"{pth.split('_')[-1]}_{v[0]}_in_{v[1] if len(v) == 2 else v[2]}"] = ( + aux[pth][v[0]].to(v[1]).magnitude + ) + + return ret + + +def main(): + parser = argparse.ArgumentParser(description="Build G frame from aux data.") + parser.add_argument( + "-i", + "--f_i3_in", + help="Path to input i3 file", + ) + parser.add_argument( + "-o", + "--f_i3_out", + default=None, + help="Path to output i3 file (needed if append flag is not set)", + ) + parser.add_argument( + "-a", + "--f_aux", + help="Path to aux file", + ) + parser.add_argument( + "-k", + "--key", + default=None, + help="run key (if omitted key is inferred from i3 file name)", + ) + parser.add_argument("-e", "--expand", action="store_true", help="Expand existing i3 file") + parser.add_argument("-c", "--channel_mask", default=0xFFFF, type=int, help="Channel mask") + parser.add_argument( + "-s", + "--scale", + default=0.25, + type=float, + help="ADC scaling factor in mV (default 0.25 mV)", + ) + + args = parser.parse_args() + + if not args.expand: + if args.f_i3_out is None: + msg = "Output not set. Either use append flag or specify output path." + raise ValueError(msg) + if os.path.isfile(args.f_i3_out): + msg = "Output file exists." + raise ValueError(msg) + + i3_split = args.f_i3_in.split(".") + f_i3_out = i3_split[0] + "_aux." + "".join(i3_split[1:]) if args.expand else args.f_i3_out + f_log = f_i3_out.split(".")[0] + ".log" + + logger = setup_logging(log_file=f_log, level=logging.INFO) + + # if key is None see if we can find it in the path of the aux file + # we are looking for YYYY_MM_DD_HH_MM_SS + key = args.key + if key is None: + match = re.search(r"\d{4}_\d{2}_\d{2}_\d{2}_\d{2}_\d{2}", args.f_i3) + if match is not None: + key = match.group() + else: + msg = "Key not provided and can not infer from aux path" + raise ValueError(msg) + + # Generate i3 file + msg = f"Writing aux info to G frame {args.f_aux}" + logger.info(msg) + tray = icetray.I3Tray() + tray.AddModule("I3Reader", "reader", filename=args.f_i3_in) + tray.AddModule( + injectDetectorInfo, config=load_aux(Path(args.f_aux), key=key), scale=args.scale + ) + tray.AddModule("I3Writer", "writer", Filename=f_i3_out) + tray.Execute() + tray.Finish() + + if args.expand: + os.rename(f_i3_out, args.f_i3_in) + + msg = f"Finished processing. i3 file {f_i3_out if not args.expand else args.f_i3_in} generated" + logger.info(msg) + + +if __name__ == "__main__": + main() From 535c0291cd27c787453b2a8592221f6ddab8c918 Mon Sep 17 00:00:00 2001 From: Patrick Krause Date: Mon, 5 Jan 2026 18:27:39 -0800 Subject: [PATCH 02/16] add p frame injection (charge reconstruction) --- src/mintanalysis/pmt/i3/funny_heartbeat.py | 168 +++++++++++++++++++++ src/mintanalysis/pmt/i3/inject_p_frame.py | 104 +++++++++++++ src/mintanalysis/pmt/i3/insert_aux.py | 4 +- 3 files changed, 274 insertions(+), 2 deletions(-) create mode 100644 src/mintanalysis/pmt/i3/funny_heartbeat.py create mode 100644 src/mintanalysis/pmt/i3/inject_p_frame.py diff --git a/src/mintanalysis/pmt/i3/funny_heartbeat.py b/src/mintanalysis/pmt/i3/funny_heartbeat.py new file mode 100644 index 0000000..2de5e9a --- /dev/null +++ b/src/mintanalysis/pmt/i3/funny_heartbeat.py @@ -0,0 +1,168 @@ +#!/usr/bin/env python3 +""" +funny_heartbeat.py + +Displays absurd waveform-to-physics heartbeat lines in the terminal. +- Escalates absurdity over runtime +- Occasionally prints rare easter egg lines +""" + +import random +import threading +import time +from contextlib import contextmanager +import shutil +import sys + +# --- core word lists --- +VERBS = [ + "Gently coercing", + "Politely requesting", + "Forcing", + "Convincing", + "Bribing", + "Strongly encouraging", + "Shaming", + "Reassuring", + "Whispering to", + "Negotiating with", + "Taming", + "Pleasing", +] + +OBJECTS = [ + "the low-resolution PMT trace", + "these stubborn pulses", + "the mischievous waveform", + "the noisy PMT bins", + "this opinionated signal", + "the reluctant photon counts", + "the wiggly PMT output", + "the latent pulses", + "the under-sampled peaks", +] + +PHYSICS = [ + "for its hidden amplitudes", + "to recover true photon arrival times", + "for non-negative signal components", + "to unfold underlying pulses", + "for latent photon distribution", + "to reconstruct the high-resolution waveform", + "to infer pulse heights", + "for all scientifically plausible signal shapes", +] + +METHODS = [ + "using NNLS unfolding", + "via repeated non-negative least squares", + "by iterative pulse nudging", + "with gentle spline persuasion", + "through spectral telepathy", + "by brute-force amplitude whispering", + "using photon-level negotiations", + "via zero-error curve hugging", + "with mild FFT intimidation", +] + +ENDINGS = [ + "…", + "please cooperate…", + "results are approximate™…", + "the photons seem skeptical…", + "stand by for high-resolution enlightenment…", + "this might break reality…", + "science is optional…", + "the waveform resists all logic…", +] + +# Rare “easter egg” lines +EASTER_EGGS = [ + "🎉 The PMTs are throwing a photon party!", + "🛸 Aliens are tweaking your NNLS bins.", + "🧙‍♂️ Wizard applied magic to the unfolding.", + "💥 Chaos ensues in the photon distribution!", + "🐱 Cats are supervising the waveform reconstruction.", + "🌈 Rainbow pulses detected!", + "🪄 Unfolding complete: miracles observed.", +] + +SPINNER = ["⠋","⠙","⠹","⠸","⠼","⠴","⠦","⠧","⠇","⠏"] + +# --- generate a line with escalating absurdity --- +def funny_line(elapsed_seconds: float) -> str: + """Return an absurd waveform-to-physics line, escalation based on elapsed time.""" + # Escalate absurdity by expanding METHODS and ENDINGS over time + extra_methods = [ + "while loudly chanting to the PMTs", + "by negotiating with imaginary photons", + "through time-traveling pulse inference" + ] + extra_endings = [ + "prepare for interdimensional interference…", + "physics may be optional here…", + "please wait while reality stabilizes…" + ] + + methods_pool = METHODS + (extra_methods if elapsed_seconds > 60 else []) + endings_pool = ENDINGS + (extra_endings if elapsed_seconds > 120 else []) + + line = f"{random.choice(VERBS)} {random.choice(OBJECTS)} {random.choice(PHYSICS)} ({random.choice(methods_pool)}). {random.choice(endings_pool)}" + + # Very rare easter egg (1% chance) + if random.random() < 0.01: + line = random.choice(EASTER_EGGS) + + return line + +# --- heartbeat thread --- +def funny_heartbeat(stop_event, line_interval: float = 10.0, spinner_interval: float = 0.1): + start_time = time.time() + last_line_time = 0 + current_line = funny_line(0) + spinner_idx = 0 + + while not stop_event.is_set(): + elapsed = time.time() - start_time + + # Update the line every line_interval + if elapsed - last_line_time >= line_interval: + current_line = funny_line(elapsed) + last_line_time = elapsed + + # Update spinner + spinner_char = SPINNER[spinner_idx % len(SPINNER)] + spinner_idx += 1 + + # Get terminal width + width = shutil.get_terminal_size((120, 20)).columns + + # Prepare full text and pad to clear old content + full_text = f"{spinner_char} {current_line}" + padded_text = full_text[:width-1].ljust(width-1) # ensure fits terminal + + sys.stdout.write("\r" + padded_text) + sys.stdout.flush() + time.sleep(spinner_interval) + + # Clear line on exit + sys.stdout.write("\r" + " " * 120 + "\r") + sys.stdout.flush() + + +# --- context manager for easy use --- +@contextmanager +def heartbeat_ctx(line_interval: float = 10.0, spinner_interval: float = 0.1): + stop = threading.Event() + t = threading.Thread(target=funny_heartbeat, args=(stop, line_interval, spinner_interval), daemon=True) + t.start() + try: + yield + finally: + stop.set() + t.join() + +if __name__ == "__main__": + + with heartbeat_ctx(): + time.sleep(300) \ No newline at end of file diff --git a/src/mintanalysis/pmt/i3/inject_p_frame.py b/src/mintanalysis/pmt/i3/inject_p_frame.py new file mode 100644 index 0000000..5c91734 --- /dev/null +++ b/src/mintanalysis/pmt/i3/inject_p_frame.py @@ -0,0 +1,104 @@ +import os +import argparse +import copy +import logging +from pathlib import Path + +from mintanalysis.pmt.ana.utils import setup_logging +from mintanalysis.pmt.i3.funny_heartbeat import heartbeat_ctx + +from icecube import icetray, pone_unfolding + +class injectEmptyPhysics(icetray.I3Module): + def __init__(self, context): + icetray.I3Module.__init__(self, context) + self.ctr = 0 + + def Process(self): + self.current_frame = self.PopFrame() + if self.current_frame.Stop == icetray.I3Frame.DAQ: + if self.ctr >= 1: + self.Inject() + self.ctr += 1 + self.PushFrame(self.current_frame) + + def Inject(self): + # Inject empty physics frame + frame = icetray.I3Frame(icetray.I3Frame.Physics) + self.PushFrame(frame) + + +class WaveFormUnfold(pone_unfolding.WaveUnfold): + def DAQ(self, frame): + self.PushFrame(frame) + + + def Physics(self,frame): + ch_check_key= "P1ChannelCalibration" + channels = frame[ch_check_key].keys() + + # Tricky: Delete waveforms for channels not currently being used, to avoid needing calibration for them. + raw_data = copy.copy(frame["RawData"]) + for c in [c for c in raw_data.keys() if c not in channels]: + del raw_data[c] + frame.Delete("RawData") + frame.Put("RawData", raw_data, icetray.I3Frame.Stream('Q')) + + super().DAQ(frame) + + + + +def main(): + parser = argparse.ArgumentParser(description="Build P frame.") + parser.add_argument( + "-i", + "--f_i3_in", + help="Path to input i3 file", + ) + parser.add_argument( + "-o", + "--f_i3_out", + default=None, + help="Path to output i3 file (needed if append flag is not set)", + ) + parser.add_argument("-a", "--append", action="store_true", help="Append existing i3 file") + parser.add_argument("-u", "--upsample", type=int, default=100, help="Upsample factor (spe per bin)") + parser.add_argument("-t", "--tolerance", type=float, default=2, help="Chi^2 stopping tolerance") + parser.add_argument("-w", "--waveform", default="", help="Name of refolded waveform field (no waveforms stored if omitted or empty string)") + + args = parser.parse_args() + + if not args.append: + if args.f_i3_out is None: + msg = "Output not set. Either use append flag or specify output path." + raise ValueError(msg) + if os.path.isfile(args.f_i3_out): + msg = "Output file exists." + raise ValueError(msg) + + i3_split = args.f_i3_in.split(".") + f_i3_out = i3_split[0] + "_aux." + "".join(i3_split[1:]) if args.append else args.f_i3_out + f_log = f_i3_out.split(".")[0] + ".log" + + logger = setup_logging(log_file=f_log, level=logging.INFO) + + msg = f"Creating P frame. This may take a while" + logger.info(msg) + with heartbeat_ctx(): + tray = icetray.I3Tray() + tray.AddModule('I3Reader', 'reader',filename=args.f_i3_in) + tray.AddModule(injectEmptyPhysics) + tray.AddModule(WaveFormUnfold,SPEsPerBin=args.upsample, OutputWaveforms="", RefoldedWaveforms=args.waveform, Tolerance=args.tolerance) + tray.AddModule("I3Writer", 'writer', Filename=f_i3_out) + tray.Execute() + tray.Finish() + + if args.append: + os.rename(f_i3_out, args.f_i3_in) + + msg = f"Finished processing. i3 file {f_i3_out if not args.append else args.f_i3_in} generated" + logger.info(msg) + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/src/mintanalysis/pmt/i3/insert_aux.py b/src/mintanalysis/pmt/i3/insert_aux.py index bae30b0..4e03e73 100644 --- a/src/mintanalysis/pmt/i3/insert_aux.py +++ b/src/mintanalysis/pmt/i3/insert_aux.py @@ -5,7 +5,7 @@ from pathlib import Path import yaml -from icecube import dataclasses, icetray, p1_dataclasses +from icecube import dataclasses, icetray, p1_dataclasses, dataio from pint import UnitRegistry from mintanalysis.pmt.ana.utils import get_physics_object, setup_logging @@ -188,7 +188,7 @@ def main(): tray = icetray.I3Tray() tray.AddModule("I3Reader", "reader", filename=args.f_i3_in) tray.AddModule( - injectDetectorInfo, config=load_aux(Path(args.f_aux), key=key), scale=args.scale + injectDetectorInfo, config=load_aux(Path(args.f_aux), key=key, ch_mask = args.channel_mask), scale=args.scale, ) tray.AddModule("I3Writer", "writer", Filename=f_i3_out) tray.Execute() From 2f08fc5baeba43e69cce7497548cfd0be8729eb9 Mon Sep 17 00:00:00 2001 From: Patrick Krause Date: Mon, 5 Jan 2026 18:30:27 -0800 Subject: [PATCH 03/16] pre-commit fixes --- src/mintanalysis/pmt/i3/funny_heartbeat.py | 35 ++++++------ src/mintanalysis/pmt/i3/inject_p_frame.py | 63 +++++++++++++--------- src/mintanalysis/pmt/i3/insert_aux.py | 6 ++- 3 files changed, 58 insertions(+), 46 deletions(-) diff --git a/src/mintanalysis/pmt/i3/funny_heartbeat.py b/src/mintanalysis/pmt/i3/funny_heartbeat.py index 2de5e9a..d4b4989 100644 --- a/src/mintanalysis/pmt/i3/funny_heartbeat.py +++ b/src/mintanalysis/pmt/i3/funny_heartbeat.py @@ -1,18 +1,9 @@ -#!/usr/bin/env python3 -""" -funny_heartbeat.py - -Displays absurd waveform-to-physics heartbeat lines in the terminal. -- Escalates absurdity over runtime -- Occasionally prints rare easter egg lines -""" - import random +import shutil +import sys import threading import time from contextlib import contextmanager -import shutil -import sys # --- core word lists --- VERBS = [ @@ -87,7 +78,8 @@ "🪄 Unfolding complete: miracles observed.", ] -SPINNER = ["⠋","⠙","⠹","⠸","⠼","⠴","⠦","⠧","⠇","⠏"] +SPINNER = ["⠋", "⠙", "⠹", "⠸", "⠼", "⠴", "⠦", "⠧", "⠇", "⠏"] + # --- generate a line with escalating absurdity --- def funny_line(elapsed_seconds: float) -> str: @@ -96,18 +88,19 @@ def funny_line(elapsed_seconds: float) -> str: extra_methods = [ "while loudly chanting to the PMTs", "by negotiating with imaginary photons", - "through time-traveling pulse inference" + "through time-traveling pulse inference", ] extra_endings = [ "prepare for interdimensional interference…", "physics may be optional here…", - "please wait while reality stabilizes…" + "please wait while reality stabilizes…", ] methods_pool = METHODS + (extra_methods if elapsed_seconds > 60 else []) endings_pool = ENDINGS + (extra_endings if elapsed_seconds > 120 else []) - line = f"{random.choice(VERBS)} {random.choice(OBJECTS)} {random.choice(PHYSICS)} ({random.choice(methods_pool)}). {random.choice(endings_pool)}" + line = f"{random.choice(VERBS)} {random.choice(OBJECTS)} {random.choice(PHYSICS)}" + line += f" ({random.choice(methods_pool)}). {random.choice(endings_pool)}" # Very rare easter egg (1% chance) if random.random() < 0.01: @@ -115,6 +108,7 @@ def funny_line(elapsed_seconds: float) -> str: return line + # --- heartbeat thread --- def funny_heartbeat(stop_event, line_interval: float = 10.0, spinner_interval: float = 0.1): start_time = time.time() @@ -134,12 +128,12 @@ def funny_heartbeat(stop_event, line_interval: float = 10.0, spinner_interval: f spinner_char = SPINNER[spinner_idx % len(SPINNER)] spinner_idx += 1 - # Get terminal width + # Get terminal width width = shutil.get_terminal_size((120, 20)).columns # Prepare full text and pad to clear old content full_text = f"{spinner_char} {current_line}" - padded_text = full_text[:width-1].ljust(width-1) # ensure fits terminal + padded_text = full_text[: width - 1].ljust(width - 1) # ensure fits terminal sys.stdout.write("\r" + padded_text) sys.stdout.flush() @@ -154,7 +148,9 @@ def funny_heartbeat(stop_event, line_interval: float = 10.0, spinner_interval: f @contextmanager def heartbeat_ctx(line_interval: float = 10.0, spinner_interval: float = 0.1): stop = threading.Event() - t = threading.Thread(target=funny_heartbeat, args=(stop, line_interval, spinner_interval), daemon=True) + t = threading.Thread( + target=funny_heartbeat, args=(stop, line_interval, spinner_interval), daemon=True + ) t.start() try: yield @@ -162,7 +158,8 @@ def heartbeat_ctx(line_interval: float = 10.0, spinner_interval: float = 0.1): stop.set() t.join() + if __name__ == "__main__": with heartbeat_ctx(): - time.sleep(300) \ No newline at end of file + time.sleep(300) diff --git a/src/mintanalysis/pmt/i3/inject_p_frame.py b/src/mintanalysis/pmt/i3/inject_p_frame.py index 5c91734..e1b569b 100644 --- a/src/mintanalysis/pmt/i3/inject_p_frame.py +++ b/src/mintanalysis/pmt/i3/inject_p_frame.py @@ -1,19 +1,19 @@ -import os import argparse import copy import logging -from pathlib import Path +import os + +from icecube import icetray, pone_unfolding from mintanalysis.pmt.ana.utils import setup_logging from mintanalysis.pmt.i3.funny_heartbeat import heartbeat_ctx -from icecube import icetray, pone_unfolding class injectEmptyPhysics(icetray.I3Module): def __init__(self, context): icetray.I3Module.__init__(self, context) self.ctr = 0 - + def Process(self): self.current_frame = self.PopFrame() if self.current_frame.Stop == icetray.I3Frame.DAQ: @@ -21,34 +21,31 @@ def Process(self): self.Inject() self.ctr += 1 self.PushFrame(self.current_frame) - + def Inject(self): # Inject empty physics frame - frame = icetray.I3Frame(icetray.I3Frame.Physics) + frame = icetray.I3Frame(icetray.I3Frame.Physics) self.PushFrame(frame) - + class WaveFormUnfold(pone_unfolding.WaveUnfold): def DAQ(self, frame): self.PushFrame(frame) - - def Physics(self,frame): - ch_check_key= "P1ChannelCalibration" + def Physics(self, frame): + ch_check_key = "P1ChannelCalibration" channels = frame[ch_check_key].keys() - - # Tricky: Delete waveforms for channels not currently being used, to avoid needing calibration for them. + + # Delete waveforms for channels not currently being used raw_data = copy.copy(frame["RawData"]) - for c in [c for c in raw_data.keys() if c not in channels]: + for c in [c for c in raw_data if c not in channels]: del raw_data[c] frame.Delete("RawData") - frame.Put("RawData", raw_data, icetray.I3Frame.Stream('Q')) + frame.Put("RawData", raw_data, icetray.I3Frame.Stream("Q")) super().DAQ(frame) - - def main(): parser = argparse.ArgumentParser(description="Build P frame.") parser.add_argument( @@ -63,9 +60,18 @@ def main(): help="Path to output i3 file (needed if append flag is not set)", ) parser.add_argument("-a", "--append", action="store_true", help="Append existing i3 file") - parser.add_argument("-u", "--upsample", type=int, default=100, help="Upsample factor (spe per bin)") - parser.add_argument("-t", "--tolerance", type=float, default=2, help="Chi^2 stopping tolerance") - parser.add_argument("-w", "--waveform", default="", help="Name of refolded waveform field (no waveforms stored if omitted or empty string)") + parser.add_argument( + "-u", "--upsample", type=int, default=100, help="Upsample factor (spe per bin)" + ) + parser.add_argument( + "-t", "--tolerance", type=float, default=2, help="Chi^2 stopping tolerance" + ) + parser.add_argument( + "-w", + "--waveform", + default="", + help="Name of refolded waveform field (no waveforms stored if omitted or empty string)", + ) args = parser.parse_args() @@ -76,21 +82,27 @@ def main(): if os.path.isfile(args.f_i3_out): msg = "Output file exists." raise ValueError(msg) - + i3_split = args.f_i3_in.split(".") f_i3_out = i3_split[0] + "_aux." + "".join(i3_split[1:]) if args.append else args.f_i3_out f_log = f_i3_out.split(".")[0] + ".log" logger = setup_logging(log_file=f_log, level=logging.INFO) - msg = f"Creating P frame. This may take a while" + msg = "Creating P frame. This may take a while" logger.info(msg) with heartbeat_ctx(): tray = icetray.I3Tray() - tray.AddModule('I3Reader', 'reader',filename=args.f_i3_in) + tray.AddModule("I3Reader", "reader", filename=args.f_i3_in) tray.AddModule(injectEmptyPhysics) - tray.AddModule(WaveFormUnfold,SPEsPerBin=args.upsample, OutputWaveforms="", RefoldedWaveforms=args.waveform, Tolerance=args.tolerance) - tray.AddModule("I3Writer", 'writer', Filename=f_i3_out) + tray.AddModule( + WaveFormUnfold, + SPEsPerBin=args.upsample, + OutputWaveforms="", + RefoldedWaveforms=args.waveform, + Tolerance=args.tolerance, + ) + tray.AddModule("I3Writer", "writer", Filename=f_i3_out) tray.Execute() tray.Finish() @@ -100,5 +112,6 @@ def main(): msg = f"Finished processing. i3 file {f_i3_out if not args.append else args.f_i3_in} generated" logger.info(msg) + if __name__ == "__main__": - main() \ No newline at end of file + main() diff --git a/src/mintanalysis/pmt/i3/insert_aux.py b/src/mintanalysis/pmt/i3/insert_aux.py index 4e03e73..8cdce18 100644 --- a/src/mintanalysis/pmt/i3/insert_aux.py +++ b/src/mintanalysis/pmt/i3/insert_aux.py @@ -5,7 +5,7 @@ from pathlib import Path import yaml -from icecube import dataclasses, icetray, p1_dataclasses, dataio +from icecube import dataclasses, icetray, p1_dataclasses from pint import UnitRegistry from mintanalysis.pmt.ana.utils import get_physics_object, setup_logging @@ -188,7 +188,9 @@ def main(): tray = icetray.I3Tray() tray.AddModule("I3Reader", "reader", filename=args.f_i3_in) tray.AddModule( - injectDetectorInfo, config=load_aux(Path(args.f_aux), key=key, ch_mask = args.channel_mask), scale=args.scale, + injectDetectorInfo, + config=load_aux(Path(args.f_aux), key=key, ch_mask=args.channel_mask), + scale=args.scale, ) tray.AddModule("I3Writer", "writer", Filename=f_i3_out) tray.Execute() From 7b8a26f9ad01e675e55bca521ba3673da28d8da3 Mon Sep 17 00:00:00 2001 From: Patrick Krause Date: Wed, 7 Jan 2026 17:59:13 -0800 Subject: [PATCH 04/16] Split out histogram loading into own function --- .../pmt/ana/peSpectrumAnalyzer.py | 45 +++++++++++-------- 1 file changed, 26 insertions(+), 19 deletions(-) diff --git a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py index 79b042e..98e9ee4 100644 --- a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py +++ b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py @@ -243,24 +243,11 @@ def analyze_run(self, run_name: str, meta: dict[str, Any]) -> dict[int, dict[str self.logger.info("Wrote PDF for run %s to %s", run_name, pdf_path) return run_results - # ---------------------- - # Per-channel processing - # ---------------------- - def process_channel( + def _get_histo( self, - run_name: str, ch: str, - meta: dict[str, Any], - f_raw: Path, f_dsp: Path, - ) -> tuple[plt.Figure | None, dict[str, Any]]: - """Process channel. Returns (figure_or_None, channel_result_dict). - - Non-critical failures return a result dict with status 'skipped' - """ - result: dict[str, Any] = {} - ch_idx = int(ch[2:]) - result["voltage"] = meta[ch_idx]["v10"] + ) -> tuple[np.NDArray[Any], np.NDArray[Any]]: # load data try: @@ -277,13 +264,33 @@ def process_channel( msg = f"Failed to compute pe values for {ch}: {e}" self.logger.warning(msg) return None, {"status": "skipped", "reason": msg} + return np.histogram(pe_vals, bins=self.bins) + + # ---------------------- + # Per-channel processing + # ---------------------- + def process_channel( + self, + run_name: str, + ch: str, + meta: dict[str, Any], + f_raw: Path, + f_dsp: Path, + ) -> tuple[plt.Figure | None, dict[str, Any]]: + """Process channel. Returns (figure_or_None, channel_result_dict). + + Non-critical failures return a result dict with status 'skipped' + """ + ch_idx = int(ch[2:]) + result: dict[str, Any] = {} + result["voltage"] = meta[ch_idx]["v10"] # histogram fig, ax = plt.subplots(figsize=A4_LANDSCAPE) - n, bins, _ = ax.hist( - pe_vals, - bins=self.bins, - histtype="step", + n, bins = self._get_histo(ch, f_dsp) + n, bins, _ = ax.stairs( + values=n, + edges=bins, label=f"channel {ch} (PMT {ch_idx+1}) at {result['voltage'].magnitude:.2f}" f" {format(result['voltage'].units,'~')}", ) From 5d424b17aac520253dcdcea89aee51f7c3b55964 Mon Sep 17 00:00:00 2001 From: Patrick Krause Date: Wed, 7 Jan 2026 18:00:00 -0800 Subject: [PATCH 05/16] small bug fix --- src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py index 98e9ee4..7923690 100644 --- a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py +++ b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py @@ -288,7 +288,7 @@ def process_channel( # histogram fig, ax = plt.subplots(figsize=A4_LANDSCAPE) n, bins = self._get_histo(ch, f_dsp) - n, bins, _ = ax.stairs( + ax.stairs( values=n, edges=bins, label=f"channel {ch} (PMT {ch_idx+1}) at {result['voltage'].magnitude:.2f}" From 93e85d8b85f3531a20fe16093d8df4dbc997ad99 Mon Sep 17 00:00:00 2001 From: Patrick Krause Date: Wed, 7 Jan 2026 18:22:11 -0800 Subject: [PATCH 06/16] split out further --- .../pmt/ana/peSpectrumAnalyzer.py | 32 +++++++++---------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py index 7923690..0785778 100644 --- a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py +++ b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py @@ -228,7 +228,9 @@ def analyze_run(self, run_name: str, meta: dict[str, Any]) -> dict[int, dict[str self.logger.info("Run %s - channel %s (PMT %d)", run_name, ch, ch_idx + 1) try: - fig, chan_data = self.process_channel(run_name, ch, meta, f_raw, f_dsp) + n, bins = self._get_histo(ch, f_dsp) + raw_runtime = self._extract_runtime_if_present(f_raw, ch) + fig, chan_data = self.process_channel(run_name, ch, meta, n, bins, raw_runtime) # fig may be None if plotting skipped if fig is not None: pdf.savefig(fig) @@ -272,26 +274,25 @@ def _get_histo( def process_channel( self, run_name: str, - ch: str, + ch_idx: int, meta: dict[str, Any], - f_raw: Path, - f_dsp: Path, + n: np.NDArray[int], + bins: np.NDArray[float], + raw_runtime: float | None = None, ) -> tuple[plt.Figure | None, dict[str, Any]]: """Process channel. Returns (figure_or_None, channel_result_dict). Non-critical failures return a result dict with status 'skipped' """ - ch_idx = int(ch[2:]) result: dict[str, Any] = {} result["voltage"] = meta[ch_idx]["v10"] # histogram fig, ax = plt.subplots(figsize=A4_LANDSCAPE) - n, bins = self._get_histo(ch, f_dsp) ax.stairs( values=n, edges=bins, - label=f"channel {ch} (PMT {ch_idx+1}) at {result['voltage'].magnitude:.2f}" + label=f"channel {ch_idx} (PMT {ch_idx+1}) at {result['voltage'].magnitude:.2f}" f" {format(result['voltage'].units,'~')}", ) @@ -309,7 +310,6 @@ def process_channel( result["runtime"] = {} if "runtime" in meta: result["runtime"]["aux"] = meta["runtime"] - raw_runtime = self._extract_runtime_if_present(f_raw, ch) if raw_runtime is not None: result["runtime"]["raw"] = raw_runtime * self.ureg.seconds @@ -327,7 +327,7 @@ def process_channel( vi = valley_index_strict(n) if vi is None: msg = "Valley detection failed (no strict peak/valley)." - self.logger.warning("Run %s ch %s: %s", run_name, ch, msg) + self.logger.warning("Run %s ch %i: %s", run_name, ch_idx, msg) self._decorate_axis(ax) result["status"] = ("skipped",) result["reason"] = (msg,) @@ -340,7 +340,7 @@ def process_channel( pe_vi = np.argmax(sub) if pe_vi is None: msg = "1st-p.e. detection failed after noise peak." - self.logger.warning("Run %s ch %s: %s", run_name, ch, msg) + self.logger.warning("Run %s ch %i: %s", run_name, ch_idx, msg) self._decorate_axis(ax) return fig, {"status": "skipped", "reason": msg} @@ -359,13 +359,13 @@ def process_channel( if len(bin_centers_fit) < 3: msg = f"Insufficient bins for fitting (n_fit={len(bin_centers_fit)})." - self.logger.warning("Run %s ch %s: %s", run_name, ch, msg) + self.logger.warning("Run %s ch %i: %s", run_name, ch_idx, msg) self._decorate_axis(ax) return fig, {"status": "skipped", "reason": msg} amp0 = float(n[pe_peak_idx]) mu0 = float(bin_centers[pe_peak_idx]) - sigma0 = 100.0 + sigma0 = mu0 - bins[valley_idx] try: m = Minuit( @@ -377,14 +377,14 @@ def process_channel( m.errordef = Minuit.LIKELIHOOD m.migrad(iterate=10) except Exception as e: - msg = f"Minuit error for {ch}: {e}" - self.logger.warning("Run %s ch %s: %s", run_name, ch, msg) + msg = f"Minuit error for {ch_idx}: {e}" + self.logger.warning("Run %s ch %i: %s", run_name, ch_idx, msg) self._decorate_axis(ax) return fig, {"status": "skipped", "reason": msg} # Basic validity check if not getattr(m, "valid", True): - self.logger.warning("Minuit invalid result for run %s ch %s", run_name, ch) + self.logger.warning("Minuit invalid result for run %s ch %i", run_name, ch_idx) self._decorate_axis(ax) return fig, {"status": "skipped", "reason": "Minuit invalid"} @@ -442,7 +442,7 @@ def process_channel( } except Exception as e: self.logger.warning( - "Statistics computation failed for run %s ch %s: %s", run_name, ch, e + "Statistics computation failed for run %s ch %i: %s", run_name, ch_idx, e ) result.setdefault("statistics", {}) result["statistics"]["error"] = str(e) From 54a274dea85db46720723b62bf6420c90e8261b1 Mon Sep 17 00:00:00 2001 From: Patrick Krause Date: Fri, 9 Jan 2026 16:11:21 -0800 Subject: [PATCH 07/16] split up path setup into a helper function --- .../pmt/ana/peSpectrumAnalyzer.py | 26 ++++++++++++------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py index 0785778..b755a1b 100644 --- a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py +++ b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py @@ -82,15 +82,14 @@ def __init__( self.aux_yaml = aux_yaml self.keys = keys - self.plot_folder = self.aux_yaml.parent / "../ana/plots" - self.result_yaml = self.aux_yaml.parent / "../ana/results.yaml" + self._set_up_paths() self.bin_size = bin_size self.bins = np.arange(-100, 10000, bin_size) self.lim = lim self.override_results = override_results self.hemispheres = {} self.logger = logger or setup_logging() - self.plot_folder.mkdir(parents=True, exist_ok=True) + self.calibrator = calibrator if calibrator is None: self.calibrator = Calibration(24 / 240, 0.25e-3, 75.0, 4.8e-9, 1) @@ -105,10 +104,6 @@ def __init__( st = self.ureg.Quantity(self.calibrator.sampling_time, self.ureg.seconds) imp = self.ureg.Quantity(self.calibrator.adc_impedance, self.ureg.ohm) - nnls_coloumb_factor = (vadc * usr * st * rf) / imp - self.ureg.define(f"NNLS = {nnls_coloumb_factor.to('coulomb').magnitude} * coulomb") - self.ureg.define(f"ADC = {usr.magnitude}*NNLS") - cal_dict = { "vadc": vadc, "upsampling_ratio": usr, @@ -116,10 +111,21 @@ def __init__( "sampling_time": st, "adc_impedance": imp, } - self._save_results(cal_dict, "calibration_constants") + + nnls_coloumb_factor = (vadc * usr * st * rf) / imp + self.ureg.define(f"NNLS = {nnls_coloumb_factor.to('coulomb').magnitude} * coulomb") + self.ureg.define(f"ADC = {usr.magnitude}*NNLS") + + if self.calib != "None": + self._save_results(cal_dict, "calibration_constants") self._save_results(self.hemispheres, "hemispheres") + def _set_up_paths(self): + self.plot_folder = self.aux_yaml.parent / "../ana/plots" + self.result_yaml = self.aux_yaml.parent / "../ana/results.yaml" + self.plot_folder.mkdir(parents=True, exist_ok=True) + # ---------------------- # I/O helpers # ---------------------- @@ -164,9 +170,9 @@ def _save_results(self, results: dict[str, dict[int, dict[str, Any]]], key: str) with open(self.result_yaml) as f: existing = yaml.safe_load(f) or {} if key in existing and not self.override_results: - msg = key + " results already present and override flag is False." + msg = key + " results already present and override flag is False. Not saving" self.logger.error(msg) - raise RuntimeError(msg) + return existing[key] = quantity_to_dict(results) with open(self.result_yaml, "w") as f: From 79cfdf9980d746888ef81cbe72b129f4704707eb Mon Sep 17 00:00:00 2001 From: Patrick Krause Date: Fri, 9 Jan 2026 16:22:38 -0800 Subject: [PATCH 08/16] Cancle linear gain fit if not at least 2 datapoints are available (and throw warning on less than 3 points) --- src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py index b755a1b..fabb887 100644 --- a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py +++ b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py @@ -689,6 +689,15 @@ def compute_dark_count_rate(self, time_mode: str = "aux") -> None: def compute_linear_gain_fit(self) -> None: data = self._load_results() tmp_dic = {"used_keys": [], "temperatures": []} + + if len(data["pe_spectrum"]) < 2: + msg = f"Only found {len(data['pe_spectrum'])} datapoints <2 in {self.result_yaml}" + self.logger.error(msg) + return + if len(data["pe_spectrum"]) < 3: + msg = f"Only found 2 datapoints in {self.result_yaml}. Fit will have no uncertainties." + self.logger.warning(msg) + for key, run in data["pe_spectrum"].items(): tmp_dic["used_keys"].append(key) tmp_dic["temperatures"].append( From 4e3755aa70573043b0ef37987b5dc1331ac65c11 Mon Sep 17 00:00:00 2001 From: Patrick Krause Date: Fri, 9 Jan 2026 16:48:53 -0800 Subject: [PATCH 09/16] Add i3 datastream analyser class --- .../pmt/ana/peSpectrumAnalyzer.py | 176 +++++++++++++++++- 1 file changed, 166 insertions(+), 10 deletions(-) diff --git a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py index fabb887..95c10b4 100644 --- a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py +++ b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py @@ -11,6 +11,11 @@ - linear gain fit - SNR plot +IceTraySpectrumAnalyzer + +Child of PESpectrumAnalyzer to do analysis based on an icecube file +Needs icetray installed + Usage: Run via CLI with desired flags """ @@ -25,6 +30,7 @@ import matplotlib.pyplot as plt import numpy as np import yaml +from icecube import icetray from iminuit import Minuit from lgdo import lh5 from matplotlib.backends.backend_pdf import PdfPages @@ -807,6 +813,135 @@ def convert_charge_units(d): self.logger.info(msg) +class IceTraySpectrumAnalyzer(PESpectrumAnalyzer): + def __init__( + self, + aux_yaml, + f_i3: Path | None = None, + keys=None, + bin_size=0.1, + lim=0, + override_results=False, + logger=None, + calibrator=None, + calib="None", + ): + super().__init__( + aux_yaml, keys, bin_size, lim, override_results, logger, calibrator, calib + ) + self.bins = np.arange(0, 10, bin_size) + self.f_i3 = f_i3 + icetray.set_log_level(icetray.I3LogLevel.LOG_WARN) + + def _set_up_paths(self): + self.plot_folder = self.aux_yaml.parent / "../ana/i3_plots" + self.result_yaml = self.aux_yaml.parent / "../ana/i3_results.yaml" + self.plot_folder.mkdir(parents=True, exist_ok=True) + + def analyze_run(self, run_name, meta): + # build file paths + if self.f_i3 is None: + f_i3 = self.aux_yaml.parent / Path( + meta["daq"].replace("daq", "i3").replace("data", "i3") + ) + else: + f_i3 = self.f_i3 + if not f_i3.exists(): + msg = f"Raw file for run {run_name} not found: {f_i3}" + raise FileNotFoundError(msg) + + run_results: dict[int, dict[str, Any]] = {} + pdf_path = self.plot_folder / f"pe_spectra_{run_name}.pdf" + + charge_counts = {} + + def histogram_charges(frame, charge_counts, pulses_name="UnfoldedPulses"): + + if not frame.Has(pulses_name): + return + pulses = frame.Get(pulses_name) + + for channel, ch_pulses in pulses.items(): + if channel.pmt not in charge_counts: + charge_counts[channel.pmt] = [] + + charge = 0 + for pulse in ch_pulses: + if pulse.charge > self.lim: + charge += pulse.charge + charge_counts[channel.pmt].append(charge) + + tray = icetray.I3Tray() + tray.AddModule("I3Reader", "reader", filename=str(f_i3)) + tray.AddModule( + histogram_charges, charge_counts=charge_counts, Streams=[icetray.I3Frame.Physics] + ) + tray.Execute() + tray.Finish() + + # collect figures and write once + with PdfPages(pdf_path) as pdf: + for ch_idx, v in sorted(charge_counts.items()): + if ch_idx not in meta: + self.logger.warning( + "Run %s: channel %s (PMT %d) not in aux file, but data exists. Skipping", + run_name, + ch_idx, + ch_idx + 1, + ) + continue + + self.logger.info("Run %s - channel %s (PMT %d)", run_name, ch_idx, ch_idx + 1) + try: + n, bins = np.histogram(v, bins=self.bins) + raw_runtime = self._extract_runtime_if_present(str(f_i3)) + fig, chan_data = self.process_channel( + run_name, ch_idx, meta, n, bins, raw_runtime + ) + # fig may be None if plotting skipped + if fig is not None: + pdf.savefig(fig) + plt.close(fig) + run_results[ch_idx] = chan_data + except Exception as exc: + self.logger.exception( + "Channel-level error run=%s ch=%s: %s", run_name, ch_idx, exc + ) + run_results[ch_idx] = {"status": "error", "reason": str(exc)} + + self.logger.info("Wrote PDF for run %s to %s", run_name, pdf_path) + return run_results + + def _extract_runtime_if_present(self, f_i3): + times = [] + + def get_runtime(frame, times): + if not frame.Has("I3EventHeader"): + return + # Assume here all data is taken in the same year + times.append(frame.Get("I3EventHeader").start_time.utc_daq_time) + + tray = icetray.I3Tray() + tray.AddModule("I3Reader", "reader", filename=str(f_i3)) + tray.AddModule(get_runtime, times=times, Streams=[icetray.I3Frame.DAQ]) + tray.Execute() + tray.Finish() + return (max(times) - min(times)) * 1e-10 + + def _decorate_axis(self, ax: plt.Axes) -> None: + xmin = np.min(self.bins) + xmin = xmin - 0.05 * xmin + xmax = np.max(self.bins) + xmax = xmax + 0.05 * xmin + ax.set_xlim(xmin, xmax) + ax.set_ylim(0.5, None) + ax.set_yscale("log") + ax.set_ylabel(f"Counts/{self.bin_size} NNLS") + ax.set_xlabel("Charge (NNLS)") + ax.set_title(f"I3 reconstruction ({self.lim} units solution cut-off)") + ax.legend() + + # -------------------------- # CLI entrypoint # -------------------------- @@ -856,6 +991,12 @@ def main() -> None: parser.add_argument( "-o", "--override", action="store_true", help="Override results if existing" ) + parser.add_argument( + "-i", + "--icetray", + action="store_true", + help="Analysis of i3 datastream (Otherwise pygama datastream)", + ) parser.add_argument( "-k", "--keys", @@ -866,20 +1007,35 @@ def main() -> None: args = parser.parse_args() - f_log = Path(args.aux_file).parent / "../ana/analysis.log" + if args.icetray: + f_log = Path(args.aux_file).parent / "../ana/i3_analysis.log" + else: + f_log = Path(args.aux_file).parent / "../ana/analysis.log" + f_log.parent.mkdir(parents=True, exist_ok=True) logger = setup_logging(log_file=f_log, level=logging.INFO) try: - analyzer = PESpectrumAnalyzer( - logger=logger, - aux_yaml=Path(args.aux_file), - keys=args.keys, - bin_size=args.bin_size, - lim=args.nnls_limit, - override_results=args.override, - calib=args.calibrate, - ) + if args.icetray: + analyzer = IceTraySpectrumAnalyzer( + logger=logger, + aux_yaml=Path(args.aux_file), + keys=args.keys, + bin_size=args.bin_size, + lim=args.nnls_limit, + override_results=args.override, + calib=args.calibrate, + ) + else: + analyzer = PESpectrumAnalyzer( + logger=logger, + aux_yaml=Path(args.aux_file), + keys=args.keys, + bin_size=args.bin_size, + lim=args.nnls_limit, + override_results=args.override, + calib=args.calibrate, + ) if args.compute_pe: analyzer.run() if args.compute_dcr: From 7b3d4f432bad4c7b2506afd4157528f615cf1946 Mon Sep 17 00:00:00 2001 From: Patrick Krause Date: Wed, 18 Mar 2026 14:58:32 -0700 Subject: [PATCH 10/16] sync --- pyproject.toml | 4 ++++ src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py | 4 ++-- src/mintanalysis/pmt/i3/daq_to_q_frame.py | 2 +- src/mintanalysis/pmt/i3/insert_aux.py | 5 ++--- src/mintanalysis/pmt/raw/build_raw.py | 2 +- 5 files changed, 10 insertions(+), 7 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 7eec043..9953090 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -101,6 +101,10 @@ build_nnls_db = "mintanalysis.pmt.dsp.build_nnls_database:build_nnls_database_cl build_raw = "mintanalysis.pmt.raw:main" +build_i3_daq = "mintanalysis.pmt.i3.daq_to_q_frame:main" +build_i3_aux = "mintanalysis.pmt.i3.insert_aux:main" +build_i3_phy = "mintanalysis.pmt.i3.inject_p_frame:main" + #[[tool.uv.index]] # Optional name for the index. #name = "dspeed-pone" diff --git a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py index 95c10b4..9a22d40 100644 --- a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py +++ b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py @@ -242,7 +242,7 @@ def analyze_run(self, run_name: str, meta: dict[str, Any]) -> dict[int, dict[str try: n, bins = self._get_histo(ch, f_dsp) raw_runtime = self._extract_runtime_if_present(f_raw, ch) - fig, chan_data = self.process_channel(run_name, ch, meta, n, bins, raw_runtime) + fig, chan_data = self.process_channel(run_name, ch_idx, meta, n, bins, raw_runtime) # fig may be None if plotting skipped if fig is not None: pdf.savefig(fig) @@ -466,7 +466,7 @@ def process_channel( # Utilities # ---------------------- def _decorate_axis(self, ax: plt.Axes) -> None: - ax.set_xlim(-10, 2.5e3) + ax.set_xlim(-10, self.bins[-1]) ax.set_ylim(0.5, None) ax.set_yscale("log") ax.set_ylabel(f"Counts/{self.bin_size} NNLS") diff --git a/src/mintanalysis/pmt/i3/daq_to_q_frame.py b/src/mintanalysis/pmt/i3/daq_to_q_frame.py index f017c93..febd18c 100644 --- a/src/mintanalysis/pmt/i3/daq_to_q_frame.py +++ b/src/mintanalysis/pmt/i3/daq_to_q_frame.py @@ -2,7 +2,7 @@ import logging import os -from icecube import icetray, pone_unfolding +from icecube import icetray, pone_unfolding, dataio from mintanalysis.pmt.ana.utils import setup_logging diff --git a/src/mintanalysis/pmt/i3/insert_aux.py b/src/mintanalysis/pmt/i3/insert_aux.py index 8cdce18..14fc928 100644 --- a/src/mintanalysis/pmt/i3/insert_aux.py +++ b/src/mintanalysis/pmt/i3/insert_aux.py @@ -5,7 +5,7 @@ from pathlib import Path import yaml -from icecube import dataclasses, icetray, p1_dataclasses +from icecube import dataclasses, icetray, p1_dataclasses, dataio from pint import UnitRegistry from mintanalysis.pmt.ana.utils import get_physics_object, setup_logging @@ -34,7 +34,7 @@ def Inject(self): frame = icetray.I3Frame(icetray.I3Frame.Geometry) if "channels" not in self.config: return - channels = [icetray.OMKey(0, 1, i) for i in self.config["channels"]] + channels = [icetray.OMKey(1, 1, i) for i in self.config["channels"]] for k, v in self.config.items(): if isinstance(v, list): map = dataclasses.I3MapKeyDouble() @@ -79,7 +79,6 @@ def load_aux(aux_yaml: Path, key: str, ch_mask: int = 0xFFFF) -> dict: # convert to physics units aux = get_physics_object(aux, ureg) - aux = aux[key] ret = {"channels": [], "channel_thresholds": [], "channel_means": [], "noise_levels": []} diff --git a/src/mintanalysis/pmt/raw/build_raw.py b/src/mintanalysis/pmt/raw/build_raw.py index 2e2bda8..37f6d82 100644 --- a/src/mintanalysis/pmt/raw/build_raw.py +++ b/src/mintanalysis/pmt/raw/build_raw.py @@ -247,7 +247,7 @@ def main(): daq_ext = args.f_daq.split(".")[-1] f_raw = ( - args.f_daq.replace("daq", "raw").replace(daq_ext, "lh5") + args.f_daq.replace("daq", "raw").replace("."+daq_ext, ".lh5") if args.f_raw is None else args.f_raw ) From c3925b15354dc01ec56f4627054aa6bb020a6118 Mon Sep 17 00:00:00 2001 From: Patrick Krause Date: Fri, 20 Mar 2026 22:21:06 -0700 Subject: [PATCH 11/16] feinschliff --- pyproject.toml | 16 +- src/mintanalysis/pmt/ana/anaCli.py | 138 +++++++++ .../pmt/ana/peSpectrumAnalyzer.py | 262 +----------------- .../pmt/ana/peSpectrumAnalyzerIcetray.py | 141 ++++++++++ src/mintanalysis/pmt/dsp/build_dsp.py | 2 +- src/mintanalysis/pmt/i3/daq_to_q_frame.py | 8 +- src/mintanalysis/pmt/i3/inject_p_frame.py | 4 +- src/mintanalysis/pmt/i3/insert_aux.py | 6 +- src/mintanalysis/pmt/raw/build_raw.py | 4 +- 9 files changed, 309 insertions(+), 272 deletions(-) create mode 100644 src/mintanalysis/pmt/ana/anaCli.py create mode 100644 src/mintanalysis/pmt/ana/peSpectrumAnalyzerIcetray.py diff --git a/pyproject.toml b/pyproject.toml index 9953090..c79db60 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -94,16 +94,16 @@ docs = [ [project.scripts] # Format: = ":" -build_ana = "mintanalysis.pmt.ana.peSpectrumAnalyzer:main" -upload_ana = "mintanalysis.pmt.ana.uploadToDB:main" -build_dsp = "mintanalysis.pmt.dsp.build_dsp:build_dsp_cli" -build_nnls_db = "mintanalysis.pmt.dsp.build_nnls_database:build_nnls_database_cli" +mint_build_ana = "mintanalysis.pmt.ana.anaCli:main" +mint_upload_ana = "mintanalysis.pmt.ana.uploadToDB:main" +mint_build_dsp = "mintanalysis.pmt.dsp.build_dsp:build_dsp_cli" +mint_build_nnls_db = "mintanalysis.pmt.dsp.build_nnls_database:build_nnls_database_cli" -build_raw = "mintanalysis.pmt.raw:main" +mint_build_raw = "mintanalysis.pmt.raw:main" -build_i3_daq = "mintanalysis.pmt.i3.daq_to_q_frame:main" -build_i3_aux = "mintanalysis.pmt.i3.insert_aux:main" -build_i3_phy = "mintanalysis.pmt.i3.inject_p_frame:main" +mint_build_i3_daq = "mintanalysis.pmt.i3.daq_to_q_frame:main" +mint_build_i3_aux = "mintanalysis.pmt.i3.insert_aux:main" +mint_build_i3_phy = "mintanalysis.pmt.i3.inject_p_frame:main" #[[tool.uv.index]] # Optional name for the index. diff --git a/src/mintanalysis/pmt/ana/anaCli.py b/src/mintanalysis/pmt/ana/anaCli.py new file mode 100644 index 0000000..e701699 --- /dev/null +++ b/src/mintanalysis/pmt/ana/anaCli.py @@ -0,0 +1,138 @@ +from __future__ import annotations + +import argparse +import logging +from pathlib import Path + +from .utils import setup_logging + +# -------------------------- +# CLI entrypoint +# -------------------------- + + +def main() -> None: + parser = argparse.ArgumentParser( + description="PE Spectrum Analyzer with DCR and linear gain fit" + ) + parser.add_argument( + "-p", "--compute-pe", action="store_true", help="Do p.e. spectrum analysis" + ) + parser.add_argument( + "-d", "--compute-dcr", action="store_true", help="Compute DCR after analysis" + ) + parser.add_argument( + "-g", "--compute-gain", action="store_true", help="Compute linear gain fit after analysis" + ) + parser.add_argument( + "-s", "--compute-snr", action="store_true", help="Compute SNR after analysis" + ) + parser.add_argument( + "-c", + "--calibrate", + default="None", + help="Choose a charge calibration unit", + ) + parser.add_argument( + "-b", + "--bin_size", + type=int, + default=20, + help="Number of bins used for analysis", + ) + parser.add_argument( + "-l", + "--nnls_limit", + type=float, + default=20, + help="Lower limit for solutions in the NNLS solution vector to be accepted.", + ) + parser.add_argument( + "-a", + "--aux_file", + help="Path to auxiliary file", + ) + parser.add_argument( + "-o", "--override", action="store_true", help="Override results if existing" + ) + parser.add_argument( + "-i", + "--icetray", + action="store_true", + help="Analysis of i3 datastream (Otherwise pygama datastream)", + ) + parser.add_argument( + "-k", + "--keys", + nargs="+", + default=None, + help="Only analyse this keys, ignore all other keys in aux file", + ) + parser.add_argument( + "--ignore_keys", + nargs="+", + default=None, + help="Ignore these keys in aux file", + ) + + args = parser.parse_args() + + if args.icetray: + f_log = Path(args.aux_file).parent / "../ana/i3_analysis.log" + else: + f_log = Path(args.aux_file).parent / "../ana/analysis.log" + + f_log.parent.mkdir(parents=True, exist_ok=True) + + logger = setup_logging(log_file=f_log, level=logging.INFO) + try: + if args.icetray: + from mintanalysis.pmt.ana.peSpectrumAnalyzerIcetray import IceTraySpectrumAnalyzer + + analyzer = IceTraySpectrumAnalyzer( + logger=logger, + aux_yaml=Path(args.aux_file), + keys=args.keys, + ignore_keys=args.ignore_keys, + bin_size=args.bin_size, + lim=args.nnls_limit, + override_results=args.override, + calib=args.calibrate, + ) + else: + from mintanalysis.pmt.ana.peSpectrumAnalyzer import PESpectrumAnalyzer + + analyzer = PESpectrumAnalyzer( + logger=logger, + aux_yaml=Path(args.aux_file), + keys=args.keys, + ignore_keys=args.ignore_keys, + bin_size=args.bin_size, + lim=args.nnls_limit, + override_results=args.override, + calib=args.calibrate, + ) + if args.compute_pe: + analyzer.run() + if args.compute_dcr: + analyzer.compute_dark_count_rate(time_mode="aux") + if args.compute_gain: + analyzer.compute_linear_gain_fit() + if args.compute_snr: + analyzer.compute_snr() + if args.calibrate != "None": + analyzer.calibrate_nnls_values( + output_file=str(analyzer.result_yaml).replace(".yaml", "_calibrated.yaml") + ) + + logger.info("Processing complete.") + except Exception as e: + logger.exception("Fatal error: %s", e) + raise + + +# -------------------------- +# CLI entrypoint +# -------------------------- +if __name__ == "__main__": + main() diff --git a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py index 9a22d40..b2b79a5 100644 --- a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py +++ b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py @@ -21,7 +21,6 @@ from __future__ import annotations -import argparse import logging from dataclasses import dataclass from pathlib import Path @@ -30,7 +29,6 @@ import matplotlib.pyplot as plt import numpy as np import yaml -from icecube import icetray from iminuit import Minuit from lgdo import lh5 from matplotlib.backends.backend_pdf import PdfPages @@ -78,6 +76,7 @@ def __init__( self, aux_yaml: Path, keys: list | None = None, + ignore_keys: list | None = None, bin_size: int = 20, lim: float = 20, override_results: bool = False, @@ -87,7 +86,7 @@ def __init__( ) -> None: self.aux_yaml = aux_yaml self.keys = keys - + self.ignore_keys = ignore_keys self._set_up_paths() self.bin_size = bin_size self.bins = np.arange(-100, 10000, bin_size) @@ -157,6 +156,10 @@ def _load_aux(self) -> dict[str, Any]: msg = f"Key {k} not in aux file, skipping." self.logger.warning(msg) + for k in ret: + if k in self.ignore_keys(): + del ret[k] + return ret def _load_results(self) -> dict: @@ -242,7 +245,9 @@ def analyze_run(self, run_name: str, meta: dict[str, Any]) -> dict[int, dict[str try: n, bins = self._get_histo(ch, f_dsp) raw_runtime = self._extract_runtime_if_present(f_raw, ch) - fig, chan_data = self.process_channel(run_name, ch_idx, meta, n, bins, raw_runtime) + fig, chan_data = self.process_channel( + run_name, ch_idx, meta, n, bins, raw_runtime + ) # fig may be None if plotting skipped if fig is not None: pdf.savefig(fig) @@ -811,252 +816,3 @@ def convert_charge_units(d): msg = f"Calibrated results written to {output_file}" self.logger.info(msg) - - -class IceTraySpectrumAnalyzer(PESpectrumAnalyzer): - def __init__( - self, - aux_yaml, - f_i3: Path | None = None, - keys=None, - bin_size=0.1, - lim=0, - override_results=False, - logger=None, - calibrator=None, - calib="None", - ): - super().__init__( - aux_yaml, keys, bin_size, lim, override_results, logger, calibrator, calib - ) - self.bins = np.arange(0, 10, bin_size) - self.f_i3 = f_i3 - icetray.set_log_level(icetray.I3LogLevel.LOG_WARN) - - def _set_up_paths(self): - self.plot_folder = self.aux_yaml.parent / "../ana/i3_plots" - self.result_yaml = self.aux_yaml.parent / "../ana/i3_results.yaml" - self.plot_folder.mkdir(parents=True, exist_ok=True) - - def analyze_run(self, run_name, meta): - # build file paths - if self.f_i3 is None: - f_i3 = self.aux_yaml.parent / Path( - meta["daq"].replace("daq", "i3").replace("data", "i3") - ) - else: - f_i3 = self.f_i3 - if not f_i3.exists(): - msg = f"Raw file for run {run_name} not found: {f_i3}" - raise FileNotFoundError(msg) - - run_results: dict[int, dict[str, Any]] = {} - pdf_path = self.plot_folder / f"pe_spectra_{run_name}.pdf" - - charge_counts = {} - - def histogram_charges(frame, charge_counts, pulses_name="UnfoldedPulses"): - - if not frame.Has(pulses_name): - return - pulses = frame.Get(pulses_name) - - for channel, ch_pulses in pulses.items(): - if channel.pmt not in charge_counts: - charge_counts[channel.pmt] = [] - - charge = 0 - for pulse in ch_pulses: - if pulse.charge > self.lim: - charge += pulse.charge - charge_counts[channel.pmt].append(charge) - - tray = icetray.I3Tray() - tray.AddModule("I3Reader", "reader", filename=str(f_i3)) - tray.AddModule( - histogram_charges, charge_counts=charge_counts, Streams=[icetray.I3Frame.Physics] - ) - tray.Execute() - tray.Finish() - - # collect figures and write once - with PdfPages(pdf_path) as pdf: - for ch_idx, v in sorted(charge_counts.items()): - if ch_idx not in meta: - self.logger.warning( - "Run %s: channel %s (PMT %d) not in aux file, but data exists. Skipping", - run_name, - ch_idx, - ch_idx + 1, - ) - continue - - self.logger.info("Run %s - channel %s (PMT %d)", run_name, ch_idx, ch_idx + 1) - try: - n, bins = np.histogram(v, bins=self.bins) - raw_runtime = self._extract_runtime_if_present(str(f_i3)) - fig, chan_data = self.process_channel( - run_name, ch_idx, meta, n, bins, raw_runtime - ) - # fig may be None if plotting skipped - if fig is not None: - pdf.savefig(fig) - plt.close(fig) - run_results[ch_idx] = chan_data - except Exception as exc: - self.logger.exception( - "Channel-level error run=%s ch=%s: %s", run_name, ch_idx, exc - ) - run_results[ch_idx] = {"status": "error", "reason": str(exc)} - - self.logger.info("Wrote PDF for run %s to %s", run_name, pdf_path) - return run_results - - def _extract_runtime_if_present(self, f_i3): - times = [] - - def get_runtime(frame, times): - if not frame.Has("I3EventHeader"): - return - # Assume here all data is taken in the same year - times.append(frame.Get("I3EventHeader").start_time.utc_daq_time) - - tray = icetray.I3Tray() - tray.AddModule("I3Reader", "reader", filename=str(f_i3)) - tray.AddModule(get_runtime, times=times, Streams=[icetray.I3Frame.DAQ]) - tray.Execute() - tray.Finish() - return (max(times) - min(times)) * 1e-10 - - def _decorate_axis(self, ax: plt.Axes) -> None: - xmin = np.min(self.bins) - xmin = xmin - 0.05 * xmin - xmax = np.max(self.bins) - xmax = xmax + 0.05 * xmin - ax.set_xlim(xmin, xmax) - ax.set_ylim(0.5, None) - ax.set_yscale("log") - ax.set_ylabel(f"Counts/{self.bin_size} NNLS") - ax.set_xlabel("Charge (NNLS)") - ax.set_title(f"I3 reconstruction ({self.lim} units solution cut-off)") - ax.legend() - - -# -------------------------- -# CLI entrypoint -# -------------------------- - - -def main() -> None: - parser = argparse.ArgumentParser( - description="PE Spectrum Analyzer with DCR and linear gain fit" - ) - parser.add_argument( - "-p", "--compute-pe", action="store_true", help="Do p.e. spectrum analysis" - ) - parser.add_argument( - "-d", "--compute-dcr", action="store_true", help="Compute DCR after analysis" - ) - parser.add_argument( - "-g", "--compute-gain", action="store_true", help="Compute linear gain fit after analysis" - ) - parser.add_argument( - "-s", "--compute-snr", action="store_true", help="Compute SNR after analysis" - ) - parser.add_argument( - "-c", - "--calibrate", - default="None", - help="Choose a charge calibration unit", - ) - parser.add_argument( - "-b", - "--bin_size", - type=int, - default=20, - help="Number of bins used for analysis", - ) - parser.add_argument( - "-l", - "--nnls_limit", - type=float, - default=20, - help="Lower limit for solutions in the NNLS solution vector to be accepted.", - ) - parser.add_argument( - "-a", - "--aux_file", - help="Path to auxiliary file", - ) - parser.add_argument( - "-o", "--override", action="store_true", help="Override results if existing" - ) - parser.add_argument( - "-i", - "--icetray", - action="store_true", - help="Analysis of i3 datastream (Otherwise pygama datastream)", - ) - parser.add_argument( - "-k", - "--keys", - nargs="+", - default=None, - help="Only analyse this keys, ignore all other keys in aux file", - ) - - args = parser.parse_args() - - if args.icetray: - f_log = Path(args.aux_file).parent / "../ana/i3_analysis.log" - else: - f_log = Path(args.aux_file).parent / "../ana/analysis.log" - - f_log.parent.mkdir(parents=True, exist_ok=True) - - logger = setup_logging(log_file=f_log, level=logging.INFO) - try: - if args.icetray: - analyzer = IceTraySpectrumAnalyzer( - logger=logger, - aux_yaml=Path(args.aux_file), - keys=args.keys, - bin_size=args.bin_size, - lim=args.nnls_limit, - override_results=args.override, - calib=args.calibrate, - ) - else: - analyzer = PESpectrumAnalyzer( - logger=logger, - aux_yaml=Path(args.aux_file), - keys=args.keys, - bin_size=args.bin_size, - lim=args.nnls_limit, - override_results=args.override, - calib=args.calibrate, - ) - if args.compute_pe: - analyzer.run() - if args.compute_dcr: - analyzer.compute_dark_count_rate(time_mode="aux") - if args.compute_gain: - analyzer.compute_linear_gain_fit() - if args.compute_snr: - analyzer.compute_snr() - if args.calibrate != "None": - analyzer.calibrate_nnls_values( - output_file=str(analyzer.result_yaml).replace(".yaml", "_calibrated.yaml") - ) - - logger.info("Processing complete.") - except Exception as e: - logger.exception("Fatal error: %s", e) - raise - - -# -------------------------- -# CLI entrypoint -# -------------------------- -if __name__ == "__main__": - main() diff --git a/src/mintanalysis/pmt/ana/peSpectrumAnalyzerIcetray.py b/src/mintanalysis/pmt/ana/peSpectrumAnalyzerIcetray.py new file mode 100644 index 0000000..1b26b30 --- /dev/null +++ b/src/mintanalysis/pmt/ana/peSpectrumAnalyzerIcetray.py @@ -0,0 +1,141 @@ +from __future__ import annotations + +from pathlib import Path +from typing import Any + +import matplotlib.pyplot as plt +import numpy as np +from icecube import icetray +from matplotlib.backends.backend_pdf import PdfPages + +from mintanalysis.pmt.ana.peSpectrumAnalyzer import PESpectrumAnalyzer + + +class IceTraySpectrumAnalyzer(PESpectrumAnalyzer): + def __init__( + self, + aux_yaml, + f_i3: Path | None = None, + keys=None, + ignore_keys: list | None = None, + bin_size: int = 100, + lim=0, + override_results=False, + logger=None, + calibrator=None, + calib="None", + ): + super().__init__( + aux_yaml, keys, ignore_keys, bin_size, lim, override_results, logger, calibrator, calib + ) + self.bins = np.arange(0, 10, 10.0 / bin_size) + self.f_i3 = f_i3 + icetray.set_log_level(icetray.I3LogLevel.LOG_WARN) + + def _set_up_paths(self): + self.plot_folder = self.aux_yaml.parent / "../ana/i3_plots" + self.result_yaml = self.aux_yaml.parent / "../ana/i3_results.yaml" + self.plot_folder.mkdir(parents=True, exist_ok=True) + + def analyze_run(self, run_name, meta): + # build file paths + if self.f_i3 is None: + f_i3 = self.aux_yaml.parent / Path( + meta["daq"].replace("daq", "i3").replace("data", "i3") + ) + else: + f_i3 = self.f_i3 + if not f_i3.exists(): + msg = f"Raw file for run {run_name} not found: {f_i3}" + raise FileNotFoundError(msg) + + run_results: dict[int, dict[str, Any]] = {} + pdf_path = self.plot_folder / f"pe_spectra_{run_name}.pdf" + + charge_counts = {} + + def histogram_charges(frame, charge_counts, pulses_name="UnfoldedPulses"): + + if not frame.Has(pulses_name): + return + pulses = frame.Get(pulses_name) + + for channel, ch_pulses in pulses.items(): + if channel.pmt not in charge_counts: + charge_counts[channel.pmt] = [] + + charge = 0 + for pulse in ch_pulses: + if pulse.charge > self.lim: + charge += pulse.charge + charge_counts[channel.pmt].append(charge) + + tray = icetray.I3Tray() + tray.AddModule("I3Reader", "reader", filename=str(f_i3)) + tray.AddModule( + histogram_charges, charge_counts=charge_counts, Streams=[icetray.I3Frame.Physics] + ) + tray.Execute() + tray.Finish() + + # collect figures and write once + with PdfPages(pdf_path) as pdf: + for ch_idx, v in sorted(charge_counts.items()): + if ch_idx not in meta: + self.logger.warning( + "Run %s: channel %s (PMT %d) not in aux file, but data exists. Skipping", + run_name, + ch_idx, + ch_idx + 1, + ) + continue + + self.logger.info("Run %s - channel %s (PMT %d)", run_name, ch_idx, ch_idx + 1) + try: + n, bins = np.histogram(v, bins=self.bins) + raw_runtime = self._extract_runtime_if_present(str(f_i3)) + fig, chan_data = self.process_channel( + run_name, ch_idx, meta, n, bins, raw_runtime + ) + # fig may be None if plotting skipped + if fig is not None: + pdf.savefig(fig) + plt.close(fig) + run_results[ch_idx] = chan_data + except Exception as exc: + self.logger.exception( + "Channel-level error run=%s ch=%s: %s", run_name, ch_idx, exc + ) + run_results[ch_idx] = {"status": "error", "reason": str(exc)} + + self.logger.info("Wrote PDF for run %s to %s", run_name, pdf_path) + return run_results + + def _extract_runtime_if_present(self, f_i3): + times = [] + + def get_runtime(frame, times): + if not frame.Has("I3EventHeader"): + return + # Assume here all data is taken in the same year + times.append(frame.Get("I3EventHeader").start_time.utc_daq_time) + + tray = icetray.I3Tray() + tray.AddModule("I3Reader", "reader", filename=str(f_i3)) + tray.AddModule(get_runtime, times=times, Streams=[icetray.I3Frame.DAQ]) + tray.Execute() + tray.Finish() + return (max(times) - min(times)) * 1e-10 + + def _decorate_axis(self, ax: plt.Axes) -> None: + xmin = np.min(self.bins) + xmin = xmin - 0.05 * xmin + xmax = np.max(self.bins) + xmax = xmax + 0.05 * xmin + ax.set_xlim(xmin, xmax) + ax.set_ylim(0.5, None) + ax.set_yscale("log") + ax.set_ylabel(f"Counts/{self.bin_size} NNLS") + ax.set_xlabel("Charge (NNLS)") + ax.set_title(f"I3 reconstruction ({self.lim} units solution cut-off)") + ax.legend() diff --git a/src/mintanalysis/pmt/dsp/build_dsp.py b/src/mintanalysis/pmt/dsp/build_dsp.py index 608bed6..1a04a90 100644 --- a/src/mintanalysis/pmt/dsp/build_dsp.py +++ b/src/mintanalysis/pmt/dsp/build_dsp.py @@ -70,7 +70,7 @@ def build_dsp_cli(): sh.setFormatter(fmt) logger.addHandler(sh) - log_file = f_dsp.replace(f_dsp.split(".")[-1], "log") + log_file = f_dsp.replace("." + f_dsp.split(".")[-1], ".log") fh = logging.FileHandler(log_file, mode="w") fh.setLevel(log_level) fh.setFormatter(fmt) diff --git a/src/mintanalysis/pmt/i3/daq_to_q_frame.py b/src/mintanalysis/pmt/i3/daq_to_q_frame.py index febd18c..797e40c 100644 --- a/src/mintanalysis/pmt/i3/daq_to_q_frame.py +++ b/src/mintanalysis/pmt/i3/daq_to_q_frame.py @@ -2,7 +2,7 @@ import logging import os -from icecube import icetray, pone_unfolding, dataio +from icecube import icetray, pone_unfolding from mintanalysis.pmt.ana.utils import setup_logging @@ -31,7 +31,9 @@ def main(): daq_ext = args.f_daq.split(".")[-1] f_i3 = ( - args.f_daq.replace("daq", "i3").replace(daq_ext, "i3") if args.f_i3 is None else args.f_i3 + args.f_daq.replace("daq", "i3").replace("." + daq_ext, ".i3") + if args.f_i3 is None + else args.f_i3 ) # Create raw folders if not existing @@ -39,7 +41,7 @@ def main(): if dir: os.makedirs(dir, exist_ok=True) - f_log = f_i3.replace(f_i3.split(".")[-1], "log") + f_log = f_i3.replace("." + f_i3.split(".")[-1], ".log") logger = setup_logging(log_file=f_log, level=logging.INFO) # Check if I3 file exists and if overwrite flag is set diff --git a/src/mintanalysis/pmt/i3/inject_p_frame.py b/src/mintanalysis/pmt/i3/inject_p_frame.py index e1b569b..7488b99 100644 --- a/src/mintanalysis/pmt/i3/inject_p_frame.py +++ b/src/mintanalysis/pmt/i3/inject_p_frame.py @@ -61,7 +61,7 @@ def main(): ) parser.add_argument("-a", "--append", action="store_true", help="Append existing i3 file") parser.add_argument( - "-u", "--upsample", type=int, default=100, help="Upsample factor (spe per bin)" + "-u", "--upsample", type=int, default=10, help="Upsample factor (spe per bin)" ) parser.add_argument( "-t", "--tolerance", type=float, default=2, help="Chi^2 stopping tolerance" @@ -84,7 +84,7 @@ def main(): raise ValueError(msg) i3_split = args.f_i3_in.split(".") - f_i3_out = i3_split[0] + "_aux." + "".join(i3_split[1:]) if args.append else args.f_i3_out + f_i3_out = i3_split[0] + "_phy." + "".join(i3_split[1:]) if args.append else args.f_i3_out f_log = f_i3_out.split(".")[0] + ".log" logger = setup_logging(log_file=f_log, level=logging.INFO) diff --git a/src/mintanalysis/pmt/i3/insert_aux.py b/src/mintanalysis/pmt/i3/insert_aux.py index 14fc928..ac186e9 100644 --- a/src/mintanalysis/pmt/i3/insert_aux.py +++ b/src/mintanalysis/pmt/i3/insert_aux.py @@ -5,7 +5,7 @@ from pathlib import Path import yaml -from icecube import dataclasses, icetray, p1_dataclasses, dataio +from icecube import dataclasses, icetray, p1_dataclasses from pint import UnitRegistry from mintanalysis.pmt.ana.utils import get_physics_object, setup_logging @@ -34,7 +34,7 @@ def Inject(self): frame = icetray.I3Frame(icetray.I3Frame.Geometry) if "channels" not in self.config: return - channels = [icetray.OMKey(1, 1, i) for i in self.config["channels"]] + channels = [icetray.OMKey(0, 1, i) for i in self.config["channels"]] for k, v in self.config.items(): if isinstance(v, list): map = dataclasses.I3MapKeyDouble() @@ -174,7 +174,7 @@ def main(): # we are looking for YYYY_MM_DD_HH_MM_SS key = args.key if key is None: - match = re.search(r"\d{4}_\d{2}_\d{2}_\d{2}_\d{2}_\d{2}", args.f_i3) + match = re.search(r"\d{4}_\d{2}_\d{2}_\d{2}_\d{2}_\d{2}", args.f_i3_in) if match is not None: key = match.group() else: diff --git a/src/mintanalysis/pmt/raw/build_raw.py b/src/mintanalysis/pmt/raw/build_raw.py index 37f6d82..d2dd6d9 100644 --- a/src/mintanalysis/pmt/raw/build_raw.py +++ b/src/mintanalysis/pmt/raw/build_raw.py @@ -247,7 +247,7 @@ def main(): daq_ext = args.f_daq.split(".")[-1] f_raw = ( - args.f_daq.replace("daq", "raw").replace("."+daq_ext, ".lh5") + args.f_daq.replace("daq", "raw").replace("." + daq_ext, ".lh5") if args.f_raw is None else args.f_raw ) @@ -267,7 +267,7 @@ def main(): sh.setFormatter(fmt) logger.addHandler(sh) - log_file = f_raw.replace(f_raw.split(".")[-1], "log") + log_file = f_raw.replace("." + f_raw.split(".")[-1], ".log") fh = logging.FileHandler(log_file, mode="w") fh.setLevel(log_level) fh.setFormatter(fmt) From d0450af357fd698f28c47df6ce9c96918812f406 Mon Sep 17 00:00:00 2001 From: Patrick Krause Date: Sat, 21 Mar 2026 00:56:23 -0700 Subject: [PATCH 12/16] Adjust DCR calculation to Hamamatsu definition --- .../pmt/ana/peSpectrumAnalyzer.py | 32 ++++++++++++++----- .../pmt/ana/peSpectrumAnalyzerIcetray.py | 6 ++-- 2 files changed, 27 insertions(+), 11 deletions(-) diff --git a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py index b2b79a5..c7efc0c 100644 --- a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py +++ b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py @@ -426,13 +426,30 @@ def process_channel( # statistics try: + border = 0.75*fit_vals["mu"] + closest_index = np.abs(bin_centers - border).argmin() + + dcr_idx = np.abs(bin_centers-0.3*fit_vals["mu"]).argmin() + zero_idx = np.abs(bin_centers).argmin() + result["statistics"] = { "1st_pe_fit_integral_below_valley": self.ureg.Quantity( - ufloat(np.sum(y_fit[:valley_idx]), np.sum(y_fit[:valley_idx]) ** 0.5), + ufloat(np.sum(y_fit[zero_idx:valley_idx]), np.sum(y_fit[zero_idx:valley_idx]) ** 0.5), + "dimensionless", + ), + "1st_pe_fit_integral_below_dcr_low_lim": self.ureg.Quantity( + ufloat(np.sum(y_fit[zero_idx:dcr_idx]), np.sum(y_fit[zero_idx:dcr_idx]) ** 0.5), + "dimensionless", + ), + "1st_pe_fit_integral_below_border": self.ureg.Quantity( + ufloat(np.sum(y_fit[zero_idx:closest_index]), np.sum(y_fit[zero_idx:closest_index]) ** 0.5), "dimensionless", ), "cts_above_valley": self.ureg.Quantity( - ufloat(np.sum(n[:valley_idx]), np.sum(n[:valley_idx]) ** 0.5), "dimensionless" + ufloat(np.sum(n[valley_idx:]), np.sum(n[valley_idx:]) ** 0.5), "dimensionless" + ), + "cts_above_border": self.ureg.Quantity( + ufloat(np.sum(n[closest_index:]), np.sum(n[closest_index:]) ** 0.5), "dimensionless" ), "1st_pe_fit_integral": self.ureg.Quantity( ufloat(np.sum(y_fit), np.sum(y_fit) ** 0.5), "dimensionless" @@ -524,10 +541,11 @@ def _add_charge_axis(self, ax: plt.Axes, is_y) -> None: """ if self.calib == "gain": + echarge = 1.60217663E-19 label = "Gain (a.u.)" func = ( - lambda x: ((x * self.ureg.NNLS).to("C") / self.ureg.elementary_charge).m, - lambda y: ((y * self.ureg.elementary_charge).to("NNLS")).m, + lambda x: ((x * self.ureg.NNLS).to("elementary_charge").m), + lambda y: (y * self.ureg.elementary_charge).to("NNLS").m, ) else: if not self.ureg.NNLS.is_compatible_with(self.calib): @@ -660,9 +678,7 @@ def compute_dark_count_rate(self, time_mode: str = "aux") -> None: time_mode, ) continue - dcts = stats.get("1st_pe_fit_integral_below_valley") + stats.get( - "cts_above_valley" - ) + dcts = stats.get("1st_pe_fit_integral_below_border") + stats.get("cts_above_border") - stats.get("1st_pe_fit_integral_below_dcr_low_lim") runtime = runtime_info[time_mode] run_dcr[pmt] = dcts / runtime @@ -799,7 +815,7 @@ def convert_charge_units(d): convert_charge_units(v) # recurse elif isinstance(v, self.ureg.Quantity): if self.calib == "gain": - d[k] = self._unit_converter(v, "C", self.ureg.elementary_charge) + d[k] = self._unit_converter(v, self.ureg.elementary_charge)#, 1./(1.60217663E-19*self.ureg("C")) ) else: d[k] = self._unit_converter(v, self.calib) diff --git a/src/mintanalysis/pmt/ana/peSpectrumAnalyzerIcetray.py b/src/mintanalysis/pmt/ana/peSpectrumAnalyzerIcetray.py index 1b26b30..ec29752 100644 --- a/src/mintanalysis/pmt/ana/peSpectrumAnalyzerIcetray.py +++ b/src/mintanalysis/pmt/ana/peSpectrumAnalyzerIcetray.py @@ -5,7 +5,7 @@ import matplotlib.pyplot as plt import numpy as np -from icecube import icetray +from icecube import dataclasses, icetray, p1_dataclasses from matplotlib.backends.backend_pdf import PdfPages from mintanalysis.pmt.ana.peSpectrumAnalyzer import PESpectrumAnalyzer @@ -18,7 +18,7 @@ def __init__( f_i3: Path | None = None, keys=None, ignore_keys: list | None = None, - bin_size: int = 100, + bin_size: int = 200, lim=0, override_results=False, logger=None, @@ -28,7 +28,7 @@ def __init__( super().__init__( aux_yaml, keys, ignore_keys, bin_size, lim, override_results, logger, calibrator, calib ) - self.bins = np.arange(0, 10, 10.0 / bin_size) + self.bins = np.arange(0, 10, 20.0 / bin_size) self.f_i3 = f_i3 icetray.set_log_level(icetray.I3LogLevel.LOG_WARN) From 2945a594ea69eb9569fceb72e2e20b7e6c93e599 Mon Sep 17 00:00:00 2001 From: Patrick Krause Date: Sat, 21 Mar 2026 01:06:52 -0700 Subject: [PATCH 13/16] Change definition of signal to noise ratio --- analysis.log | 16 ---------------- src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py | 4 ++-- 2 files changed, 2 insertions(+), 18 deletions(-) delete mode 100644 analysis.log diff --git a/analysis.log b/analysis.log deleted file mode 100644 index 2041109..0000000 --- a/analysis.log +++ /dev/null @@ -1,16 +0,0 @@ -[2025-12-15 15:05:54,438] [PESpectrum - main] [INFO] collected measurement results -[2025-12-15 15:05:54,438] [PESpectrum - main] [INFO] collected PMT info -[2025-12-15 15:05:54,439] [PESpectrum - main] [INFO] collected environment info -[2025-12-15 15:05:54,439] [PESpectrum - main] [INFO] collected software info -[2025-12-15 15:05:58,512] [PESpectrum - _connect_tunnel] [INFO] [info] SSH tunnel established: 127.0.0.1:43307 -> 127.0.0.1:27017 via production.pacific-neutrino.org -[2025-12-15 15:05:58,620] [PESpectrum - _connect_client] [INFO] [info] Connected to MongoDB at 127.0.0.1:43307 with user p-one -[2025-12-15 15:05:58,753] [PESpectrum - main] [INFO] Build data object -[2025-12-15 15:05:58,753] [PESpectrum - main] [INFO] The following dict would be uploaded to the DB: {'measurement_type': 'dcr', 'measurement_location': 'MINT', 'devices_used': ['p-1-1-pmt-unit-KM56674_275', 'p-1-1-om-hs-01'], 'result': 26463.90079865678, 'result_unc': 30.15089782382669, 'units': 'Hz', 'pmt_info': {'v10_in_volt': 85.09310150146484, 'di10_in_mA': 2.9245901107788086, 'frac': 0.877996027469635}, 'env_info': {'temperature_in_celsius': 29.11357307688087, 'pressure_in_hpa': 1011.6844068085361, 'humidity_in_percent': 23.183785204174583, 'measurement_duration_in_s': 29.110747814178467}, 'software_info': {'framework': 'mint-xyz', 'pe_reconstruction': 'NNLS', 'sftp_path': '/mint/mint-data/p-1-1-om-hs-01/PMT', 'run_tags': '2025_12_13_03_31_52'}, 'mapping': None} -[2025-12-15 15:05:58,753] [PESpectrum - main] [INFO] collected measurement results -[2025-12-15 15:05:58,754] [PESpectrum - main] [INFO] collected PMT info -[2025-12-15 15:05:58,754] [PESpectrum - main] [INFO] collected environment info -[2025-12-15 15:05:58,754] [PESpectrum - main] [INFO] collected software info -[2025-12-15 15:05:59,167] [PESpectrum - _connect_tunnel] [INFO] [info] SSH tunnel established: 127.0.0.1:38145 -> 127.0.0.1:27017 via production.pacific-neutrino.org -[2025-12-15 15:05:59,242] [PESpectrum - _connect_client] [INFO] [info] Connected to MongoDB at 127.0.0.1:38145 with user p-one -[2025-12-15 15:05:59,379] [PESpectrum - main] [INFO] Build data object -[2025-12-15 15:05:59,380] [PESpectrum - main] [INFO] The following dict would be uploaded to the DB: {'measurement_type': 'dcr', 'measurement_location': 'MINT', 'devices_used': ['p-1-1-pmt-unit-KM54356_23', 'p-1-1-om-hs-01'], 'result': 26326.463667634333, 'result_unc': 30.07250335177778, 'units': 'Hz', 'pmt_info': {'v10_in_volt': 87.58090209960938, 'di10_in_mA': 1.8282999992370605, 'frac': 0.8174149990081787}, 'env_info': {'temperature_in_celsius': 29.11357307688087, 'pressure_in_hpa': 1011.6844068085361, 'humidity_in_percent': 23.183785204174583, 'measurement_duration_in_s': 29.110747814178467}, 'software_info': {'framework': 'mint-xyz', 'pe_reconstruction': 'NNLS', 'sftp_path': '/mint/mint-data/p-1-1-om-hs-01/PMT', 'run_tags': '2025_12_13_03_31_52'}, 'mapping': None} diff --git a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py index c7efc0c..f86072f 100644 --- a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py +++ b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py @@ -621,7 +621,7 @@ def compute_snr(self) -> None: signal = info.get("pe_peak_fit").get("amp") else: signal = info.get("statistics").get("1st_pe_guess").get("amp") - run_snr[pmt] = 1 - noise / signal + run_snr[pmt] = signal / noise except Exception as e: self.logger.warning( "Failed to compute SNR for run %s channel %s: %s", run_name, pmt, e @@ -644,7 +644,7 @@ def compute_snr(self) -> None: ax.errorbar(pmts, vals, errs, label=lbl, fmt="o") ax.set_xlabel("Channel") ax.set_ylabel("SNR (a.u.)") - ax.set_title("SNR per Channel (= 1 - valley/peak)") + ax.set_title("SNR per Channel (= peak / valley)") ax.legend() plt.tight_layout() plot_path = self.plot_folder / "snr_plot.png" From cc5e9265d4c177f1ec11950164756f14c547e714 Mon Sep 17 00:00:00 2001 From: Patrick Krause Date: Sat, 21 Mar 2026 01:09:15 -0700 Subject: [PATCH 14/16] pre-commit fixes --- .../pmt/ana/peSpectrumAnalyzer.py | 32 +++++++++++++------ .../pmt/ana/peSpectrumAnalyzerIcetray.py | 2 +- 2 files changed, 24 insertions(+), 10 deletions(-) diff --git a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py index f86072f..dbe150c 100644 --- a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py +++ b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py @@ -426,30 +426,39 @@ def process_channel( # statistics try: - border = 0.75*fit_vals["mu"] + border = 0.75 * fit_vals["mu"] closest_index = np.abs(bin_centers - border).argmin() - dcr_idx = np.abs(bin_centers-0.3*fit_vals["mu"]).argmin() + dcr_idx = np.abs(bin_centers - 0.3 * fit_vals["mu"]).argmin() zero_idx = np.abs(bin_centers).argmin() result["statistics"] = { "1st_pe_fit_integral_below_valley": self.ureg.Quantity( - ufloat(np.sum(y_fit[zero_idx:valley_idx]), np.sum(y_fit[zero_idx:valley_idx]) ** 0.5), + ufloat( + np.sum(y_fit[zero_idx:valley_idx]), + np.sum(y_fit[zero_idx:valley_idx]) ** 0.5, + ), "dimensionless", ), "1st_pe_fit_integral_below_dcr_low_lim": self.ureg.Quantity( - ufloat(np.sum(y_fit[zero_idx:dcr_idx]), np.sum(y_fit[zero_idx:dcr_idx]) ** 0.5), + ufloat( + np.sum(y_fit[zero_idx:dcr_idx]), np.sum(y_fit[zero_idx:dcr_idx]) ** 0.5 + ), "dimensionless", ), "1st_pe_fit_integral_below_border": self.ureg.Quantity( - ufloat(np.sum(y_fit[zero_idx:closest_index]), np.sum(y_fit[zero_idx:closest_index]) ** 0.5), + ufloat( + np.sum(y_fit[zero_idx:closest_index]), + np.sum(y_fit[zero_idx:closest_index]) ** 0.5, + ), "dimensionless", ), "cts_above_valley": self.ureg.Quantity( ufloat(np.sum(n[valley_idx:]), np.sum(n[valley_idx:]) ** 0.5), "dimensionless" ), "cts_above_border": self.ureg.Quantity( - ufloat(np.sum(n[closest_index:]), np.sum(n[closest_index:]) ** 0.5), "dimensionless" + ufloat(np.sum(n[closest_index:]), np.sum(n[closest_index:]) ** 0.5), + "dimensionless", ), "1st_pe_fit_integral": self.ureg.Quantity( ufloat(np.sum(y_fit), np.sum(y_fit) ** 0.5), "dimensionless" @@ -541,7 +550,6 @@ def _add_charge_axis(self, ax: plt.Axes, is_y) -> None: """ if self.calib == "gain": - echarge = 1.60217663E-19 label = "Gain (a.u.)" func = ( lambda x: ((x * self.ureg.NNLS).to("elementary_charge").m), @@ -678,7 +686,11 @@ def compute_dark_count_rate(self, time_mode: str = "aux") -> None: time_mode, ) continue - dcts = stats.get("1st_pe_fit_integral_below_border") + stats.get("cts_above_border") - stats.get("1st_pe_fit_integral_below_dcr_low_lim") + dcts = ( + stats.get("1st_pe_fit_integral_below_border") + + stats.get("cts_above_border") + - stats.get("1st_pe_fit_integral_below_dcr_low_lim") + ) runtime = runtime_info[time_mode] run_dcr[pmt] = dcts / runtime @@ -815,7 +827,9 @@ def convert_charge_units(d): convert_charge_units(v) # recurse elif isinstance(v, self.ureg.Quantity): if self.calib == "gain": - d[k] = self._unit_converter(v, self.ureg.elementary_charge)#, 1./(1.60217663E-19*self.ureg("C")) ) + d[k] = self._unit_converter( + v, self.ureg.elementary_charge + ) # , 1./(1.60217663E-19*self.ureg("C")) ) else: d[k] = self._unit_converter(v, self.calib) diff --git a/src/mintanalysis/pmt/ana/peSpectrumAnalyzerIcetray.py b/src/mintanalysis/pmt/ana/peSpectrumAnalyzerIcetray.py index ec29752..31e9d77 100644 --- a/src/mintanalysis/pmt/ana/peSpectrumAnalyzerIcetray.py +++ b/src/mintanalysis/pmt/ana/peSpectrumAnalyzerIcetray.py @@ -5,7 +5,7 @@ import matplotlib.pyplot as plt import numpy as np -from icecube import dataclasses, icetray, p1_dataclasses +from icecube import icetray from matplotlib.backends.backend_pdf import PdfPages from mintanalysis.pmt.ana.peSpectrumAnalyzer import PESpectrumAnalyzer From 1d34da1d932804ba7f05fa6772c3ba4739702312 Mon Sep 17 00:00:00 2001 From: Patrick Krause Date: Sat, 21 Mar 2026 17:03:52 -0700 Subject: [PATCH 15/16] leave DB uploader open during full runtime --- src/mintanalysis/pmt/ana/uploadToDB.py | 130 +++++++++++++------------ 1 file changed, 69 insertions(+), 61 deletions(-) diff --git a/src/mintanalysis/pmt/ana/uploadToDB.py b/src/mintanalysis/pmt/ana/uploadToDB.py index dac62b7..8b5a67d 100644 --- a/src/mintanalysis/pmt/ana/uploadToDB.py +++ b/src/mintanalysis/pmt/ana/uploadToDB.py @@ -361,77 +361,85 @@ def main(): msg = f"Tag {tag} not found in the result file" logger.error(msg) raise ValueError + try: + db = ProductionDatabase(**db_opts) + while ch_mask: + hs = aux["hemisphere_a"] + ch = (ch_mask & -ch_mask).bit_length() - 1 + pmt_no = ch + 1 + if pmt_no > 8: + pmt_no -= 8 + hs = aux["hemisphere_b"] + + # collect value(s) + vals, errs, units, mapping, _keys = get_values( + data[measurement][pmt_no] + if measurement == "linear_gain" + else data[measurement][tag][pmt_no] + ) + logger.info("collected measurement results") + + # collect pmt info + pmt_info = get_pmt_info(pmt_no, aux, tag) + logger.info("collected PMT info") + + # get environment info + env_info = get_env_info(aux, tag) + logger.info("collected environment info") + + # get software info + try: + framework = version("mint-analysis") + except PackageNotFoundError: + framework = "unknown" + + sw_info = SoftwareInfo( + framework="mint_analyis_" + framework, + pe_reconstruction=reco, + sftp_path=sftp_dir, + run_tags=tag, + ) + logger.info("collected software info") - while ch_mask: - ch = (ch_mask & -ch_mask).bit_length() - 1 - pmt_no = ch + 1 - - # collect value(s) - vals, errs, units, mapping, _keys = get_values( - data[measurement][pmt_no] - if measurement == "linear_gain" - else data[measurement][tag][pmt_no] - ) - logger.info("collected measurement results") - - # collect pmt info - pmt_info = get_pmt_info(pmt_no, aux, tag) - logger.info("collected PMT info") - - # get environment info - env_info = get_env_info(aux, tag) - logger.info("collected environment info") + # get used devices + # Get settings from DB + # TODO replace with module level logic - # get software info - try: - framework = version("mint-analysis") - except PackageNotFoundError: - framework = "unknown" - - sw_info = SoftwareInfo( - framework="mint_analyis_" + framework, - pe_reconstruction=reco, - sftp_path=sftp_dir, - run_tags=tag, - ) - logger.info("collected software info") + overview = db.get_hemisphere_overview(hs, "sfu") + pmt_units = {int(p.get("position")): p for p in overview.get("pmt-unit")} - # get used devices - # Get settings from DB - # TODO replace with module level logic - with ProductionDatabase(**db_opts) as db: - overview = db.get_hemisphere_overview(str(sftp_root_dir), "tum") pmt_obj = { "device_type": "pmt-unit", - "uid": overview.get("pmt-unit").get(str(pmt_no)).get("uid"), - "_id": overview.get("pmt-unit").get(str(pmt_no)).get("_id"), + "uid": pmt_units.get(pmt_no).get("uid"), + "_id": pmt_units.get(pmt_no).get("_id"), } - # Build measurement result - res = MeasurementResult( - measurement_type=measurement, - measurement_location="MINT", - devices_used=pmt_obj, - result=vals, - result_unc=errs, - units=units, - mapping=mapping, - pmt_info=pmt_info, - env_info=env_info, - software_info=sw_info, - ) - logger.info("Build data object") + # Build measurement result + res = MeasurementResult( + measurement_type=measurement, + measurement_location="MINT", + devices_used=pmt_obj, + result=vals, + result_unc=errs, + units=units, + mapping=mapping, + pmt_info=pmt_info, + env_info=env_info, + software_info=sw_info, + ) + logger.info("Build data object") - # Upload to database - if args.dry: - msg = f"The following dict would be uploaded to the DB: {asdict(res)}" - else: - with ProductionDatabase(**db_opts) as db: + # Upload to database + if args.dry: + msg = f"The following dict would be uploaded to the DB: {asdict(res)}" + else: db.client["mint"]["Measurements_Pmt"].insert_one(asdict(res)) - msg = f"Uploaded document {asdict(res)} to mint/Measurements_Pmt" + msg = f"Uploaded document {asdict(res)} to mint/Measurements_Pmt" - logger.info(msg) - ch_mask &= ch_mask - 1 # clear lowest set bit + logger.info(msg) + ch_mask &= ch_mask - 1 # clear lowest set bit + finally: + db.disconnect() # -------------------------- From f46a8e7b5a7c66223490f3e33b30f1cceddea74d Mon Sep 17 00:00:00 2001 From: Patrick Krause Date: Sat, 21 Mar 2026 17:52:11 -0700 Subject: [PATCH 16/16] encode used icetray version if used --- src/mintanalysis/pmt/ana/uploadToDB.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/mintanalysis/pmt/ana/uploadToDB.py b/src/mintanalysis/pmt/ana/uploadToDB.py index 8b5a67d..b1a7139 100644 --- a/src/mintanalysis/pmt/ana/uploadToDB.py +++ b/src/mintanalysis/pmt/ana/uploadToDB.py @@ -241,6 +241,12 @@ def main(): default=None, help="key to upload (will be ignored for linear_gain measurement)", ) + parser.add_argument( + "-i", + "--icetray", + default=None, + help="Path to icetray shell if data has been analyzed using the icetray framework.", + ) parser.add_argument( "-r", "--reco", default="NNLS", help="Reconstruction algorithm used in DSP" ) @@ -393,6 +399,20 @@ def main(): except PackageNotFoundError: framework = "unknown" + # get icetray version from shell file if used + if args.icetray: + ice = "icetray.unknown_version" + with open(args.icetray, encoding="utf-8") as file: + for line in file: + if line.strip().startswith('printctr "Version'): + ice = "icetray" + "_".join( + line.split('printctr "Version')[-1] + .replace('"', "") + .replace("icetray", "") + .split() + ) + framework += "_" + ice + sw_info = SoftwareInfo( framework="mint_analyis_" + framework, pe_reconstruction=reco,