From 293654500972cfe83fa192c0df7a44193af116a6 Mon Sep 17 00:00:00 2001 From: Manuel Huber Date: Wed, 2 Oct 2024 00:03:12 +0200 Subject: [PATCH] g4gps: add simple cli, fix energy ordering --- pyproject.toml | 3 +++ src/legendoptics/cli.py | 34 ++++++++++++++++++++++++++++++++++ src/legendoptics/lar.py | 6 ++++-- src/legendoptics/pen.py | 6 ++++-- src/legendoptics/utils.py | 25 +++++++++++++++++++++++-- tests/test_g4gps.py | 23 +++++++++++++++++++++-- 6 files changed, 89 insertions(+), 8 deletions(-) create mode 100644 src/legendoptics/cli.py diff --git a/pyproject.toml b/pyproject.toml index 03711f8..9e27ef0 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -69,6 +69,9 @@ test = [ "pytest-cov", ] +[project.scripts] +legend-pygeom-optics = "legendoptics.cli:optics_cli" + [tool.setuptools] include-package-data = true zip-safe = false diff --git a/src/legendoptics/cli.py b/src/legendoptics/cli.py new file mode 100644 index 0000000..11d4dba --- /dev/null +++ b/src/legendoptics/cli.py @@ -0,0 +1,34 @@ +from __future__ import annotations + +import argparse + + +def optics_cli() -> None: + parser = argparse.ArgumentParser( + prog="legend-pygeom-optics", + description="%(prog)s command line interface", + ) + + subparsers = parser.add_subparsers(dest="command", required=True) + + g4gps_parser = subparsers.add_parser( + "g4gps", help="build evt file from remage hit file" + ) + g4gps_parser.add_argument( + "--type", + choices=("arb_file", "macro"), + help="file that contains a list of detector ids that are part of the input file. default: %(default)e", + default="macro", + ) + g4gps_parser.add_argument("spectrum", choices=("lar_emission", "pen_emission")) + g4gps_parser.add_argument("output", help="output file") + + args = parser.parse_args() + + if args.command == "g4gps": + if args.spectrum == "lar_emission": + from legendoptics.lar import g4gps_lar_emissions_spectrum as g4gps_spectrum + elif args.spectrum == "pen_emission": + from legendoptics.pen import g4gps_pen_emissions_spectrum as g4gps_spectrum + + g4gps_spectrum(args.output, args.type == "macro") diff --git a/src/legendoptics/lar.py b/src/legendoptics/lar.py index 38afb50..6af2cc8 100644 --- a/src/legendoptics/lar.py +++ b/src/legendoptics/lar.py @@ -490,7 +490,7 @@ def pyg4_lar_attach_scintillation( pyg4_def_scint_by_particle_type(lar_mat, scint_params) -def g4gps_lar_emissions_spectrum(filename: str) -> None: +def g4gps_lar_emissions_spectrum(filename: str, output_macro: bool) -> None: """Write a LAr emission energy spectrum for G4GeneralParticleSource. See Also @@ -508,4 +508,6 @@ def g4gps_lar_emissions_spectrum(filename: str) -> None: scint_em[0] = 0 scint_em[-1] = 0 - g4gps_write_emission_spectrum(filename, λ_peak, scint_em) + g4gps_write_emission_spectrum( + filename, output_macro, λ_peak, scint_em, "lar_emissions_spectrum" + ) diff --git a/src/legendoptics/pen.py b/src/legendoptics/pen.py index 87bf9b4..d0fcffe 100644 --- a/src/legendoptics/pen.py +++ b/src/legendoptics/pen.py @@ -222,7 +222,7 @@ def pyg4_pen_attach_scintillation(mat, reg) -> None: pyg4_def_scint_by_particle_type(mat, pen_scintillation_params()) -def g4gps_pen_emissions_spectrum(filename: str) -> None: +def g4gps_pen_emissions_spectrum(filename: str, output_macro: bool) -> None: """Write a PEN emission energy spectrum for G4GeneralParticleSource. See Also @@ -239,4 +239,6 @@ def g4gps_pen_emissions_spectrum(filename: str) -> None: scint_em[0] = 0 scint_em[-1] = 0 - g4gps_write_emission_spectrum(filename, λ_scint, scint_em) + g4gps_write_emission_spectrum( + filename, output_macro, λ_scint, scint_em, "pen_emissions_spectrum" + ) diff --git a/src/legendoptics/utils.py b/src/legendoptics/utils.py index 21c1510..6ade357 100644 --- a/src/legendoptics/utils.py +++ b/src/legendoptics/utils.py @@ -1,6 +1,7 @@ from __future__ import annotations import logging +from pathlib import Path import numpy as np import pint @@ -108,7 +109,11 @@ def __call__(self, pts: Quantity) -> Quantity: def g4gps_write_emission_spectrum( - filename: str, λ_peak: Quantity, scint_em: Quantity + filename: str, + output_macro: bool, + λ_peak: Quantity, + scint_em: Quantity, + quantity_name: str, ) -> None: """Write a energy spectrum for use with G4GeneralParticleSource @@ -125,7 +130,23 @@ def g4gps_write_emission_spectrum( with u.context("sp"): pointwise = np.array([λ_peak.to("MeV").m, scint_em.m]).T + # reorder the values to be in ascending energy order. + sort = np.argsort(pointwise[:, 0]) + pointwise[:, 0] = pointwise[:, 0][sort] + pointwise[:, 1] = pointwise[:, 1][sort] + if pointwise.shape[0] > 1024: log.warning("G4GeneralParticleSource spectrum can only have 1024 bins.") - np.savetxt(filename, pointwise) + if not output_macro: + np.savetxt(filename, pointwise) + return + + with Path.open(filename, "wt") as f: + f.write(f"# {quantity_name} | legendoptics\n\n") + f.write("/gps/ene/type Arb\n") + f.write("/gps/ene/diffspec true\n") + f.write("/gps/hist/type arb\n") + f.write("/gps/hist/inter Lin\n") + for point in pointwise: + f.write(f"/gps/hist/point {point[0]} {point[1]}\n") diff --git a/tests/test_g4gps.py b/tests/test_g4gps.py index 7788764..a955017 100644 --- a/tests/test_g4gps.py +++ b/tests/test_g4gps.py @@ -2,18 +2,37 @@ from __future__ import annotations +import numpy as np import pint u = pint.get_application_registry() +def test_g4gps(tmp_path) -> None: + from legendoptics.utils import g4gps_write_emission_spectrum + + λ_peak = np.array([200, 300, 400]) * u.nm + scint_em = np.array([10, 11, 12]) * u.dimensionless + tmp_file = tmp_path / "test.csv" + g4gps_write_emission_spectrum(tmp_file, False, λ_peak, scint_em, "test") + + saved = np.loadtxt(tmp_file) + # ensure monotonicity (increasing _energy_ order). + assert saved[0, 0] < saved[1, 0] + assert saved[1, 0] < saved[2, 0] + assert saved[0, 1] > saved[1, 1] + assert saved[1, 1] > saved[2, 1] + + def test_g4gps_lar(tmp_path) -> None: import legendoptics.lar - legendoptics.lar.g4gps_lar_emissions_spectrum(tmp_path / "lar.csv") + legendoptics.lar.g4gps_lar_emissions_spectrum(tmp_path / "lar.csv", False) + legendoptics.lar.g4gps_lar_emissions_spectrum(tmp_path / "lar.mac", True) def test_g4gps_pen(tmp_path) -> None: import legendoptics.pen - legendoptics.pen.g4gps_pen_emissions_spectrum(tmp_path / "pen.csv") + legendoptics.pen.g4gps_pen_emissions_spectrum(tmp_path / "pen.csv", False) + legendoptics.pen.g4gps_pen_emissions_spectrum(tmp_path / "pen.mac", True)