g4gps: add simple cli, fix energy ordering

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

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")

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"
+    )

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"
+    )

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")

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)