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/pyproject.toml b/pyproject.toml index 7eec043..c79db60 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -94,12 +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" + +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 79b042e..dbe150c 100644 --- a/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py +++ b/src/mintanalysis/pmt/ana/peSpectrumAnalyzer.py @@ -11,12 +11,16 @@ - 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 """ from __future__ import annotations -import argparse import logging from dataclasses import dataclass from pathlib import Path @@ -72,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, @@ -81,16 +86,15 @@ def __init__( ) -> None: 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.ignore_keys = ignore_keys + 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 +109,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 +116,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 # ---------------------- @@ -145,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: @@ -164,9 +179,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: @@ -228,7 +243,11 @@ 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_idx, meta, n, bins, raw_runtime + ) # fig may be None if plotting skipped if fig is not None: pdf.savefig(fig) @@ -243,24 +262,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,14 +283,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_idx: int, + meta: dict[str, Any], + 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' + """ + 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", - label=f"channel {ch} (PMT {ch_idx+1}) at {result['voltage'].magnitude:.2f}" + ax.stairs( + values=n, + edges=bins, + label=f"channel {ch_idx} (PMT {ch_idx+1}) at {result['voltage'].magnitude:.2f}" f" {format(result['voltage'].units,'~')}", ) @@ -302,7 +327,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 @@ -320,7 +344,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,) @@ -333,7 +357,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} @@ -352,13 +376,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( @@ -370,14 +394,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"} @@ -402,13 +426,39 @@ 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" @@ -435,7 +485,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) @@ -447,7 +497,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") @@ -502,8 +552,8 @@ def _add_charge_axis(self, ax: plt.Axes, is_y) -> None: if self.calib == "gain": 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): @@ -579,7 +629,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 @@ -602,7 +652,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" @@ -636,8 +686,10 @@ 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 @@ -676,6 +728,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( @@ -766,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, "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) @@ -783,102 +846,3 @@ def convert_charge_units(d): msg = f"Calibrated results written to {output_file}" self.logger.info(msg) - - -# -------------------------- -# 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( - "-k", - "--keys", - nargs="+", - default=None, - help="Only analyse this keys, ignore all other keys in aux file", - ) - - args = parser.parse_args() - - 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.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..31e9d77 --- /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 = 200, + 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, 20.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/ana/uploadToDB.py b/src/mintanalysis/pmt/ana/uploadToDB.py index dac62b7..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" ) @@ -361,77 +367,99 @@ 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" + + # 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, + 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() # -------------------------- 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 new file mode 100644 index 0000000..797e40c --- /dev/null +++ b/src/mintanalysis/pmt/i3/daq_to_q_frame.py @@ -0,0 +1,76 @@ +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/funny_heartbeat.py b/src/mintanalysis/pmt/i3/funny_heartbeat.py new file mode 100644 index 0000000..d4b4989 --- /dev/null +++ b/src/mintanalysis/pmt/i3/funny_heartbeat.py @@ -0,0 +1,165 @@ +import random +import shutil +import sys +import threading +import time +from contextlib import contextmanager + +# --- 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)}" + line += f" ({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) 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..7488b99 --- /dev/null +++ b/src/mintanalysis/pmt/i3/inject_p_frame.py @@ -0,0 +1,117 @@ +import argparse +import copy +import logging +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 + + +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() + + # Delete waveforms for channels not currently being used + raw_data = copy.copy(frame["RawData"]) + 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")) + + 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=10, 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] + "_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) + + 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(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() diff --git a/src/mintanalysis/pmt/i3/insert_aux.py b/src/mintanalysis/pmt/i3/insert_aux.py new file mode 100644 index 0000000..ac186e9 --- /dev/null +++ b/src/mintanalysis/pmt/i3/insert_aux.py @@ -0,0 +1,206 @@ +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_in) + 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, ch_mask=args.channel_mask), + 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() diff --git a/src/mintanalysis/pmt/raw/build_raw.py b/src/mintanalysis/pmt/raw/build_raw.py index 2e2bda8..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)