From 9c22b7f7677dd0528607cce4ac2a0dbc883b9d20 Mon Sep 17 00:00:00 2001 From: EricMEsch Date: Fri, 10 Jan 2025 14:10:21 +0100 Subject: [PATCH 1/5] Add distance to surface tests --- tests/CMakeLists.txt | 1 + tests/distances/CMakeLists.txt | 28 ++++ tests/distances/gdml/ge-array.gdml | 0 tests/distances/macros/B99000A.json | 36 +++++ tests/distances/macros/C99000A.json | 44 ++++++ tests/distances/macros/P99000A.json | 29 ++++ tests/distances/macros/V99000A.json | 44 ++++++ tests/distances/macros/test-ge-distance.mac | 17 +++ tests/distances/make_ge_gdml.py | 87 ++++++++++++ tests/distances/test_ge_distance.py | 141 ++++++++++++++++++++ 10 files changed, 427 insertions(+) create mode 100644 tests/distances/CMakeLists.txt create mode 100644 tests/distances/gdml/ge-array.gdml create mode 100644 tests/distances/macros/B99000A.json create mode 100644 tests/distances/macros/C99000A.json create mode 100644 tests/distances/macros/P99000A.json create mode 100644 tests/distances/macros/V99000A.json create mode 100644 tests/distances/macros/test-ge-distance.mac create mode 100644 tests/distances/make_ge_gdml.py create mode 100644 tests/distances/test_ge_distance.py diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 4b1ce4f0..67bfe4cf 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -9,6 +9,7 @@ get_target_property(PYTHONPATH python-virtualenv PYTHONPATH) add_subdirectory(basics) add_subdirectory(confinement) +add_subdirectory(distances) add_subdirectory(germanium) add_subdirectory(internals) add_subdirectory(output) diff --git a/tests/distances/CMakeLists.txt b/tests/distances/CMakeLists.txt new file mode 100644 index 00000000..5fe761f0 --- /dev/null +++ b/tests/distances/CMakeLists.txt @@ -0,0 +1,28 @@ +# collect auxiliary files +file( + GLOB _aux + RELATIVE ${PROJECT_SOURCE_DIR} + macros/*.mac macros/*.json gdml/*.gdml gdml/*.xml *.py) + +# copy them to the build area + +foreach(_file ${_aux}) + configure_file(${PROJECT_SOURCE_DIR}/${_file} ${PROJECT_BINARY_DIR}/${_file} COPYONLY) +endforeach() + +# generate the GDML file +add_test(NAME distances-ge/gen-gdml COMMAND ${PYTHONPATH} make_ge_gdml.py) +set_tests_properties(distances-ge/gen-gdml PROPERTIES LABELS extra FIXTURES_SETUP + distance-gdml-fixture) + +# test on HPGe containment +add_test(NAME distances-ge/gen-output COMMAND ${REMAGE_PYEXE} -g gdml/ge-array.gdml -w -o + test-distance.lh5 -- macros/test-ge-distance.mac) +set_tests_properties( + distances-ge/gen-output PROPERTIES LABELS extra FIXTURES_SETUP distance-output-fixture + FIXTURES_REQUIRED distance-gdml-fixture) + +add_test(NAME distances-ge/distance COMMAND ${PYTHONPATH} ./test_ge_distance.py) + +set_tests_properties(distances-ge/distance PROPERTIES LABELS extra FIXTURES_REQUIRED + distance-output-fixture) diff --git a/tests/distances/gdml/ge-array.gdml b/tests/distances/gdml/ge-array.gdml new file mode 100644 index 00000000..e69de29b diff --git a/tests/distances/macros/B99000A.json b/tests/distances/macros/B99000A.json new file mode 100644 index 00000000..d63577d3 --- /dev/null +++ b/tests/distances/macros/B99000A.json @@ -0,0 +1,36 @@ +{ + "name": "B99000A", + "type": "bege", + "production": { + "manufacturer": "Test", + "enrichment": { + "val": 0.75, + "unc": 0.05 + } + }, + "geometry": { + "height_in_mm": 40.0, + "radius_in_mm": 35.0, + "groove": { + "depth_in_mm": 2.0, + "radius_in_mm": { + "outer": 12.0, + "inner": 10.0 + } + }, + "pp_contact": { + "radius_in_mm": 7.5, + "depth_in_mm": 0 + }, + "taper": { + "top": { + "angle_in_deg": 0.0, + "height_in_mm": 0.0 + }, + "bottom": { + "angle_in_deg": 0.0, + "height_in_mm": 0.0 + } + } + } +} diff --git a/tests/distances/macros/C99000A.json b/tests/distances/macros/C99000A.json new file mode 100644 index 00000000..c6b58b33 --- /dev/null +++ b/tests/distances/macros/C99000A.json @@ -0,0 +1,44 @@ +{ + "name": "C99000A", + "type": "coax", + "production": { + "manufacturer": "Test", + "enrichment": { + "val": 0.75, + "unc": 0.05 + } + }, + "geometry": { + "height_in_mm": 80, + "radius_in_mm": 40, + "borehole": { + "radius_in_mm": 7, + "depth_in_mm": 70 + }, + "groove": { + "depth_in_mm": 2, + "radius_in_mm": { + "outer": 20, + "inner": 17 + } + }, + "pp_contact": { + "radius_in_mm": 17, + "depth_in_mm": 0 + }, + "taper": { + "top": { + "angle_in_deg": 45, + "height_in_mm": 5 + }, + "bottom": { + "angle_in_deg": 45, + "height_in_mm": 2 + }, + "borehole": { + "angle_in_deg": 0, + "height_in_mm": 0 + } + } + } +} diff --git a/tests/distances/macros/P99000A.json b/tests/distances/macros/P99000A.json new file mode 100644 index 00000000..ee684f2f --- /dev/null +++ b/tests/distances/macros/P99000A.json @@ -0,0 +1,29 @@ +{ + "name": "P99000A", + "type": "ppc", + "production": { + "manufacturer": "Test", + "enrichment": { + "val": 0.75, + "unc": 0.05 + } + }, + "geometry": { + "height_in_mm": 45.0, + "radius_in_mm": 35.0, + "pp_contact": { + "radius_in_mm": 2, + "depth_in_mm": 2 + }, + "taper": { + "top": { + "angle_in_deg": 0, + "height_in_mm": 0 + }, + "bottom": { + "angle_in_deg": 30, + "height_in_mm": 10 + } + } + } +} diff --git a/tests/distances/macros/V99000A.json b/tests/distances/macros/V99000A.json new file mode 100644 index 00000000..ab1b0906 --- /dev/null +++ b/tests/distances/macros/V99000A.json @@ -0,0 +1,44 @@ +{ + "name": "V99000A", + "type": "icpc", + "production": { + "manufacturer": "Test", + "enrichment": { + "val": 0.75, + "unc": 0.05 + } + }, + "geometry": { + "height_in_mm": 70, + "radius_in_mm": 35, + "borehole": { + "radius_in_mm": 5, + "depth_in_mm": 55 + }, + "groove": { + "depth_in_mm": 1, + "radius_in_mm": { + "outer": 10, + "inner": 9 + } + }, + "pp_contact": { + "radius_in_mm": 3, + "depth_in_mm": 2 + }, + "taper": { + "top": { + "angle_in_deg": 10, + "height_in_mm": 10 + }, + "bottom": { + "angle_in_deg": 0, + "height_in_mm": 0 + }, + "borehole": { + "angle_in_deg": 0, + "height_in_mm": 0 + } + } + } +} diff --git a/tests/distances/macros/test-ge-distance.mac b/tests/distances/macros/test-ge-distance.mac new file mode 100644 index 00000000..fefd0e98 --- /dev/null +++ b/tests/distances/macros/test-ge-distance.mac @@ -0,0 +1,17 @@ +/control/execute macros/detectors-fake.mac +/RMG/Output/NtuplePerDetector false + +/run/initialize + +/RMG/Generator/Confine Volume +/RMG/Generator/Confinement/Physical/AddVolume B.* +/RMG/Generator/Confinement/Physical/AddVolume C.* +/RMG/Generator/Confinement/Physical/AddVolume P.* +/RMG/Generator/Confinement/Physical/AddVolume V.* + + +/RMG/Generator/Select GPS +/gps/particle e- +/gps/energy 1 eV + +/run/beamOn 100000 diff --git a/tests/distances/make_ge_gdml.py b/tests/distances/make_ge_gdml.py new file mode 100644 index 00000000..6a266506 --- /dev/null +++ b/tests/distances/make_ge_gdml.py @@ -0,0 +1,87 @@ +from __future__ import annotations + +import json +from pathlib import Path + +import numpy as np +import pyg4ometry as pg4 +import pygeomtools as pytools +from legendhpges import make_hpge + +# read the configs +out_gdml = "gdml/ge-array.gdml" +det_macro = "macros/detectors-fake.mac" +config_dict = {} +for det_type in ["B", "P", "V", "C"]: + with Path.open(Path(f"macros/{det_type}99000A.json")) as file: + config_dict[det_type] = json.load(file) + + +def add_hpge(lar, reg, angle, radius, z, idx, dtype): + x = radius * np.sin(np.deg2rad(angle)) + y = radius * np.cos(np.deg2rad(angle)) + logical_detector = make_hpge(config_dict[dtype], name=f"{dtype}{idx}", registry=reg) + logical_detector.pygeom_color_rgba = (0, 1, 1, 0.2) + physical_detector = pg4.geant4.PhysicalVolume( + [0, 0, 0], [x, y, z], logical_detector, f"{dtype}{idx}", lar, reg + ) + + physical_detector.pygeom_active_dector = pytools.RemageDetectorInfo( + "germanium", + idx, + config_dict[dtype], + ) + return idx + 1 + + +# construct geometry +reg = pg4.geant4.Registry() +ws = pg4.geant4.solid.Box("ws", 500, 500, 500, reg, lunit="mm") +wl = pg4.geant4.LogicalVolume(ws, "G4_Galactic", "wl", reg) +wl.pygeom_color_rgba = (0.1, 1, 0.1, 0.5) + +reg.setWorld(wl) + + +# lar +lar_s = pg4.geant4.solid.Tubs( + "LAr_s", 0, 200, 250, 0, 2 * np.pi, registry=reg, lunit="mm" +) +lar_l = pg4.geant4.LogicalVolume(lar_s, "G4_lAr", "LAr_l", registry=reg) +lar_l.pygeom_color_rgba = (1, 0.1, 0, 0.2) +pg4.geant4.PhysicalVolume([0, 0, 0], [0, 0, 0], lar_l, "LAr", wl, registry=reg) + + +# hpge strings +string_radius = 85 +string_angles = [0, 90, 180, 270] +detectors = ["V", "P", "B", "C"] + + +n = 0 +lines = [] +for i, det in enumerate(detectors): + angle = string_angles[i] + n = add_hpge(lar_l, reg, angle, string_radius, -20, n, det) + + x = string_radius * np.sin(np.deg2rad(angle)) + y = string_radius * np.cos(np.deg2rad(angle)) + + lines.append("/RMG/Generator/Confinement/Geometrical/AddSolid Cylinder ") + lines.append(f"/RMG/Generator/Confinement/Geometrical/CenterPositionX {x} mm") + lines.append(f"/RMG/Generator/Confinement/Geometrical/CenterPositionY {y} mm ") + lines.append("/RMG/Generator/Confinement/Geometrical/CenterPositionZ 0 mm ") + lines.append("/RMG/Generator/Confinement/Geometrical/Cylinder/OuterRadius 44 mm ") + lines.append("/RMG/Generator/Confinement/Geometrical/Cylinder/Height 100 mm \n") + +lines_exclude = [line.replace("AddSolid", "AddExcludedSolid") for line in lines] + + +pytools.detectors.write_detector_auxvals(reg) +pytools.geometry.check_registry_sanity(reg, reg) + + +w = pg4.gdml.Writer() +w.addDetector(reg) +w.write(out_gdml) +pytools.detectors.generate_detector_macro(reg, det_macro) diff --git a/tests/distances/test_ge_distance.py b/tests/distances/test_ge_distance.py new file mode 100644 index 00000000..4b0d8f33 --- /dev/null +++ b/tests/distances/test_ge_distance.py @@ -0,0 +1,141 @@ +# test_ge_distances.py +# This test checks if the Geant4 calculated distances agree +# with the distances calculated by the HPGe class. +# Fails if the distances are not within a tolerance (default: 1 nm). + +from __future__ import annotations + +import copy + +import awkward as ak +import legendhpges as hpges +import numpy as np +import pyg4ometry as pg4 +from lgdo import lh5 +from matplotlib import pyplot as plt +from pygeomtools import get_sensvol_metadata + +plt.rcParams["lines.linewidth"] = 1 +plt.rcParams["font.size"] = 12 + +gdml = "gdml/ge-array.gdml" +outfile = "test-distance.lh5" + +# get the geometry +reg = pg4.gdml.Reader(gdml).getRegistry() +reg_tmp = pg4.geant4.Registry() +detectors = list(reg.physicalVolumeDict.keys()) + +detectors = [det for det in detectors if det[0] in ["V", "P", "B", "C"]] +# Map uid to hpge +det_map = { + det: { + "uint": int(det[1:]), + "pos": reg.physicalVolumeDict[det].position.eval(), + "hpge": hpges.make_hpge( + get_sensvol_metadata(reg, det), name=det, registry=reg_tmp + ), + } + for idx, det in enumerate(detectors) +} + +steps = lh5.read_as("stp/germanium", outfile, "ak") +# Transform the coordinates to the local detector frame +uint = np.array(np.full_like(steps.time, -1), dtype=int) +xlocal = np.array(1000 * steps.xloc) +ylocal = np.array(1000 * steps.yloc) +zlocal = np.array(1000 * steps.zloc) + +positions = np.array( + np.transpose(np.vstack([steps.xloc * 1000, steps.yloc * 1000, steps.zloc * 1000])) +) +for det in det_map: + local_positions = copy.copy(positions) + local_positions -= det_map[det]["pos"] + + is_inside = np.full(len(uint), False) + is_inside[uint == -1] = det_map[det]["hpge"].is_inside(local_positions[uint == -1]) + + uint[is_inside] = det_map[det]["uint"] + xlocal[is_inside] -= det_map[det]["pos"][0] + ylocal[is_inside] -= det_map[det]["pos"][1] + zlocal[is_inside] -= det_map[det]["pos"][2] + + +steps["xlocal"] = xlocal +steps["ylocal"] = ylocal +steps["zlocal"] = zlocal +steps = ak.unflatten(steps, ak.run_lengths(steps.evtid)) +uid = ak.fill_none(ak.firsts(steps.det_uid, axis=-1), -1) +dist_G4 = ak.fill_none(ak.firsts(steps.dist_to_surf * 1000, axis=-1), -1) +xlocal = ak.fill_none(ak.firsts(steps.xlocal, axis=-1), -1) +ylocal = ak.fill_none(ak.firsts(steps.ylocal, axis=-1), -1) +zlocal = ak.fill_none(ak.firsts(steps.zlocal, axis=-1), -1) +hits = ak.Array( + {"dist_G4": dist_G4, "uid": uid, "xloc": xlocal, "yloc": ylocal, "zloc": zlocal} +) + + +def make_plot(hit, tolerance=1e-6): + good_distance = True + for idx, det in enumerate(det_map): + temp = det_map[det]["uint"] + sel_hit = hit[hit.uid == temp] + + coords = np.column_stack((sel_hit.xloc, sel_hit.yloc, sel_hit.zloc)) + dist_py = det_map[det]["hpge"].distance_to_surface(coords) + # Check if all distances are within tolerance + if np.sum(np.abs(sel_hit.dist_G4 - dist_py) > tolerance): + good_distance = False + + # Now draw the plots + hpges.draw.plot_profile(det_map[det]["hpge"], axes=axs[idx]) + + rpos_loc = np.sqrt(sel_hit.xloc**2 + sel_hit.yloc**2) + rng = np.random.default_rng() + r = rng.choice([-1, 1], p=[0.5, 0.5], size=len(rpos_loc)) * rpos_loc + z = sel_hit.zloc + c = sel_hit.dist_G4 + cut = c < 5 + + s = axs[idx].scatter( + r[cut], + z[cut], + c=c[cut], + marker=".", + label="gen. points", + cmap="BuPu", + ) + + if idx == 0: + axs[idx].set_ylabel("Height [mm]") + c = plt.colorbar(s) + c.set_label("Distance [mm]") + + axs[idx].set_xlabel("Radius [mm]") + + return good_distance + + +num_dets = len(det_map) +cols = 2 +rows = (num_dets + (cols - 1)) // cols + +fig, axs = plt.subplots(rows, cols, figsize=(5 * cols, 5 * rows)) +axs = axs.flatten() + + +are_distances_good = make_plot(hits) +plt.suptitle("Distance check for HPGes") +caption = "The distance to the surface of the HPGe detectors, as calculated by Geant4, " +caption += "illustrated within the HPGe outline drawn with legendhpges.\n" +caption += ( + "The test will fail if the difference between the distances calculated by Geant4 " +) +caption += "compared to the distances calculated by legendhpges is not within a tolerance of (default) 1 nm." +plt.figtext(0.02, 0.04, caption, wrap=True, ha="left", fontsize=11) +plt.tight_layout(rect=[0, 0.12, 1, 1]) +# plt.tight_layout() +plt.savefig("distance-ge.output.pdf") + +assert are_distances_good From b65ee248970314f2fc4868a312f4cef9f58d8798 Mon Sep 17 00:00:00 2001 From: EricMEsch Date: Fri, 10 Jan 2025 14:36:52 +0100 Subject: [PATCH 2/5] fix imports --- tests/distances/test_ge_distance.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/distances/test_ge_distance.py b/tests/distances/test_ge_distance.py index 4b0d8f33..e251ef88 100644 --- a/tests/distances/test_ge_distance.py +++ b/tests/distances/test_ge_distance.py @@ -11,6 +11,7 @@ import legendhpges as hpges import numpy as np import pyg4ometry as pg4 +from legendhpges import draw from lgdo import lh5 from matplotlib import pyplot as plt from pygeomtools import get_sensvol_metadata @@ -89,7 +90,7 @@ def make_plot(hit, tolerance=1e-6): good_distance = False # Now draw the plots - hpges.draw.plot_profile(det_map[det]["hpge"], axes=axs[idx]) + draw.plot_profile(det_map[det]["hpge"], axes=axs[idx]) rpos_loc = np.sqrt(sel_hit.xloc**2 + sel_hit.yloc**2) rng = np.random.default_rng() From 1f1656d505c9393b257332aba6a9dad30b8726c5 Mon Sep 17 00:00:00 2001 From: EricMEsch Date: Fri, 10 Jan 2025 15:12:43 +0100 Subject: [PATCH 3/5] remove argon volume --- tests/distances/make_ge_gdml.py | 17 +++-------------- 1 file changed, 3 insertions(+), 14 deletions(-) diff --git a/tests/distances/make_ge_gdml.py b/tests/distances/make_ge_gdml.py index 6a266506..11c5993b 100644 --- a/tests/distances/make_ge_gdml.py +++ b/tests/distances/make_ge_gdml.py @@ -36,33 +36,22 @@ def add_hpge(lar, reg, angle, radius, z, idx, dtype): # construct geometry reg = pg4.geant4.Registry() -ws = pg4.geant4.solid.Box("ws", 500, 500, 500, reg, lunit="mm") +ws = pg4.geant4.solid.Box("ws", 300, 300, 300, reg, lunit="mm") wl = pg4.geant4.LogicalVolume(ws, "G4_Galactic", "wl", reg) wl.pygeom_color_rgba = (0.1, 1, 0.1, 0.5) reg.setWorld(wl) -# lar -lar_s = pg4.geant4.solid.Tubs( - "LAr_s", 0, 200, 250, 0, 2 * np.pi, registry=reg, lunit="mm" -) -lar_l = pg4.geant4.LogicalVolume(lar_s, "G4_lAr", "LAr_l", registry=reg) -lar_l.pygeom_color_rgba = (1, 0.1, 0, 0.2) -pg4.geant4.PhysicalVolume([0, 0, 0], [0, 0, 0], lar_l, "LAr", wl, registry=reg) - - # hpge strings string_radius = 85 string_angles = [0, 90, 180, 270] -detectors = ["V", "P", "B", "C"] - n = 0 lines = [] -for i, det in enumerate(detectors): +for i, det in enumerate(["V", "P", "B", "C"]): angle = string_angles[i] - n = add_hpge(lar_l, reg, angle, string_radius, -20, n, det) + n = add_hpge(wl, reg, angle, string_radius, -20, n, det) x = string_radius * np.sin(np.deg2rad(angle)) y = string_radius * np.cos(np.deg2rad(angle)) From 85e0f657d3b7859837923491489de4f240003af4 Mon Sep 17 00:00:00 2001 From: EricMEsch Date: Fri, 10 Jan 2025 16:37:23 +0100 Subject: [PATCH 4/5] Remove unnecessary lines --- tests/distances/make_ge_gdml.py | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/tests/distances/make_ge_gdml.py b/tests/distances/make_ge_gdml.py index 11c5993b..62a5464b 100644 --- a/tests/distances/make_ge_gdml.py +++ b/tests/distances/make_ge_gdml.py @@ -48,24 +48,10 @@ def add_hpge(lar, reg, angle, radius, z, idx, dtype): string_angles = [0, 90, 180, 270] n = 0 -lines = [] for i, det in enumerate(["V", "P", "B", "C"]): angle = string_angles[i] n = add_hpge(wl, reg, angle, string_radius, -20, n, det) - x = string_radius * np.sin(np.deg2rad(angle)) - y = string_radius * np.cos(np.deg2rad(angle)) - - lines.append("/RMG/Generator/Confinement/Geometrical/AddSolid Cylinder ") - lines.append(f"/RMG/Generator/Confinement/Geometrical/CenterPositionX {x} mm") - lines.append(f"/RMG/Generator/Confinement/Geometrical/CenterPositionY {y} mm ") - lines.append("/RMG/Generator/Confinement/Geometrical/CenterPositionZ 0 mm ") - lines.append("/RMG/Generator/Confinement/Geometrical/Cylinder/OuterRadius 44 mm ") - lines.append("/RMG/Generator/Confinement/Geometrical/Cylinder/Height 100 mm \n") - -lines_exclude = [line.replace("AddSolid", "AddExcludedSolid") for line in lines] - - pytools.detectors.write_detector_auxvals(reg) pytools.geometry.check_registry_sanity(reg, reg) From 10f7f51272de79dd89219a917068ee27c6aeab86 Mon Sep 17 00:00:00 2001 From: EricMEsch Date: Wed, 15 Jan 2025 17:23:16 +0100 Subject: [PATCH 5/5] Increase required accuracy --- tests/distances/test_ge_distance.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/distances/test_ge_distance.py b/tests/distances/test_ge_distance.py index e251ef88..4b7102e4 100644 --- a/tests/distances/test_ge_distance.py +++ b/tests/distances/test_ge_distance.py @@ -77,7 +77,7 @@ ) -def make_plot(hit, tolerance=1e-6): +def make_plot(hit, tolerance=1e-9): good_distance = True for idx, det in enumerate(det_map): temp = det_map[det]["uint"]