diff --git a/src/legendoptics/cli.py b/src/legendoptics/cli.py index 11d4dba..44aa802 100644 --- a/src/legendoptics/cli.py +++ b/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") diff --git a/src/legendoptics/fibers.py b/src/legendoptics/fibers.py index 47436d7..7ba1394 100644 --- a/src/legendoptics/fibers.py +++ b/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" + ) diff --git a/tests/test_g4gps.py b/tests/test_g4gps.py index a955017..b3ccb75 100644 --- a/tests/test_g4gps.py +++ b/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) diff --git a/tests/test_import.py b/tests/test_import.py index 6a9b38f..87b6134 100644 --- a/tests/test_import.py +++ b/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() diff --git a/tests/test_pyg4.py b/tests/test_pyg4.py index d05ccaa..1a98fca 100644 --- a/tests/test_pyg4.py +++ b/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: