fiber: add scintillation

src/legendoptics/cli.py

@@ -20,15 +20,19 @@ def optics_cli() -> None:
         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(
+        "spectrum", choices=("lar_emission", "pen_emission", "fiber_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
+            from legendoptics.lar import g4gps_lar_emissions_spectrum as g4gps_spec
         elif args.spectrum == "pen_emission":
-            from legendoptics.pen import g4gps_pen_emissions_spectrum as g4gps_spectrum
+            from legendoptics.pen import g4gps_pen_emissions_spectrum as g4gps_spec
+        elif args.spectrum == "fiber_emission":
+            from legendoptics.fibers import g4gps_fiber_emissions_spectrum as g4gps_spec
 
-        g4gps_spectrum(args.output, args.type == "macro")
+        g4gps_spec(args.output, args.type == "macro")

src/legendoptics/fibers.py

@@ -13,7 +13,12 @@
 from pint import Quantity
 
 from legendoptics import store
-from legendoptics.utils import InterpolatingGraph, readdatafile
+from legendoptics.scintillate import ScintConfig, ScintParticle
+from legendoptics.utils import (
+    InterpolatingGraph,
+    g4gps_write_emission_spectrum,
+    readdatafile,
+)
 
 log = logging.getLogger(__name__)
 u = pint.get_application_registry()
@@ -129,6 +134,33 @@ def fiber_absorption_path_length() -> Quantity:
     return fiber_absorption_length() * 1.21
 
 
+@store.register_pluggable
+def fiber_core_scint_light_yield() -> Quantity:
+    """fiber scintillation yield for electrons, from TODO
+
+    .. optics-const::
+    """
+    return 8000 / u.MeV
+
+
+@store.register_pluggable
+def fiber_core_scintillation_params() -> ScintConfig:
+    """Get a :class:`ScintConfig` object for fibers.
+
+    See Also
+    --------
+    .fiber_core_scint_light_yield
+    """
+    # setup scintillation response just for electrons.
+    return ScintConfig(
+        flat_top=fiber_core_scint_light_yield(),
+        fano_factor=None,
+        particles=[
+            ScintParticle("electron", yield_factor=1, exc_ratio=None),
+        ],
+    )
+
+
 def pyg4_fiber_cladding1_attach_rindex(mat, reg) -> None:
     """Attach the refractive index to the given fiber cladding 1 material instance.
 
@@ -223,3 +255,58 @@ def pyg4_fiber_core_attach_absorption(
     absorption = np.array([length.m] * λ_full.shape[0]) * length.u
     with u.context("sp"):
         mat.addVecPropertyPint("ABSLENGTH", λ_full.to("eV"), absorption)
+
+
+def pyg4_fiber_core_attach_scintillation(mat, reg) -> None:
+    """Attach Geant4 properties for fiber scintillation response to the given material instance.
+
+    .. note:: This currently only adds scintillation for energy deposited by electrons.
+
+    See Also
+    --------
+    .fiber_core_scint_light_yield
+    .fiber_wls_emission
+    .fiber_wls_timeconstant
+    """
+    from legendoptics.pyg4utils import pyg4_def_scint_by_particle_type, pyg4_sample_λ
+
+    # sample the measured emission spectrum.
+    λ_scint = pyg4_sample_λ(350 * u.nm, 650 * u.nm, 200)
+    scint_em = InterpolatingGraph(*fiber_wls_emission(), min_idx=350 * u.nm)(λ_scint)
+    # make sure that the scintillation spectrum is zero at the boundaries.
+    scint_em[0] = 0
+    scint_em[-1] = 0
+
+    with u.context("sp"):
+        mat.addVecPropertyPint("SCINTILLATIONCOMPONENT1", λ_scint.to("eV"), scint_em)
+
+    # time constant is unknown.
+    mat.addConstPropertyPint("SCINTILLATIONTIMECONSTANT1", fiber_wls_timeconstant())
+
+    # We do not know the fiber fano factor. Geant4 calculates σ = RESOLUTIONSCALE × √mean,
+    # so we set it to 1; otherwise Geant4 will crash.
+    mat.addConstPropertyPint("RESOLUTIONSCALE", 1)
+
+    pyg4_def_scint_by_particle_type(mat, fiber_core_scintillation_params())
+
+
+def g4gps_fiber_emissions_spectrum(filename: str, output_macro: bool) -> None:
+    """Write a fiber emission energy spectrum for G4GeneralParticleSource.
+
+    See Also
+    --------
+    .fiber_wls_emission
+    utils.g4gps_write_emission_spectrum
+    """
+    from legendoptics.pyg4utils import pyg4_sample_λ
+
+    # sample the measured emission spectrum.
+    λ_scint = pyg4_sample_λ(350 * u.nm, 650 * u.nm, 200)
+    scint_em = InterpolatingGraph(*fiber_wls_emission(), min_idx=350 * u.nm)(λ_scint)
+    # make sure that the scintillation spectrum is zero at the boundaries.
+    scint_em[0] = 0
+    scint_em[-1] = 0
+
+    g4gps_write_emission_spectrum(
+        filename, output_macro, λ_scint, scint_em, "fiber_emissions_spectrum"
+    )

tests/test_g4gps.py

@@ -25,14 +25,21 @@ def test_g4gps(tmp_path) -> None:
 
 
 def test_g4gps_lar(tmp_path) -> None:
-    import legendoptics.lar
+    from legendoptics import lar
 
-    legendoptics.lar.g4gps_lar_emissions_spectrum(tmp_path / "lar.csv", False)
-    legendoptics.lar.g4gps_lar_emissions_spectrum(tmp_path / "lar.mac", True)
+    lar.g4gps_lar_emissions_spectrum(tmp_path / "lar.csv", False)
+    lar.g4gps_lar_emissions_spectrum(tmp_path / "lar.mac", True)
 
 
 def test_g4gps_pen(tmp_path) -> None:
-    import legendoptics.pen
+    from legendoptics import pen
 
-    legendoptics.pen.g4gps_pen_emissions_spectrum(tmp_path / "pen.csv", False)
-    legendoptics.pen.g4gps_pen_emissions_spectrum(tmp_path / "pen.mac", True)
+    pen.g4gps_pen_emissions_spectrum(tmp_path / "pen.csv", False)
+    pen.g4gps_pen_emissions_spectrum(tmp_path / "pen.mac", True)
+
+
+def test_g4gps_fibers(tmp_path) -> None:
+    from legendoptics import fibers
+
+    fibers.g4gps_fiber_emissions_spectrum(tmp_path / "fibers.csv", False)
+    fibers.g4gps_fiber_emissions_spectrum(tmp_path / "fibers.mac", True)

tests/test_import.py

@@ -46,6 +46,8 @@ def test_import():
     legendoptics.fibers.fiber_wls_timeconstant()
     legendoptics.fibers.fiber_absorption_length()
     legendoptics.fibers.fiber_absorption_path_length()
+    legendoptics.fibers.fiber_core_scint_light_yield()
+    legendoptics.fibers.fiber_core_scintillation_params()
 
     legendoptics.tetratex.tetratex_reflectivity()
 
@@ -67,3 +69,4 @@ def test_import():
     legendoptics.pen.pen_wls_absorption()
     legendoptics.pen.pen_wls_absorption()
     legendoptics.pen.pen_wls_absorption()
+    legendoptics.pen.pen_scintillation_params()

tests/test_pyg4.py

@@ -43,6 +43,7 @@ def test_pyg4_attach_fibers() -> None:
     legendoptics.fibers.pyg4_fiber_core_attach_rindex(mat, reg)
     legendoptics.fibers.pyg4_fiber_core_attach_wls(mat, reg)
     legendoptics.fibers.pyg4_fiber_core_attach_absorption(mat, reg)
+    legendoptics.fibers.pyg4_fiber_core_attach_scintillation(mat, reg)
 
 
 def test_pyg4_attach_tetratex() -> None: