diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 4b1ce4f..67bfe4c 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 0000000..5fe761f --- /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 0000000..e69de29 diff --git a/tests/distances/macros/B99000A.json b/tests/distances/macros/B99000A.json new file mode 100644 index 0000000..d63577d --- /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 0000000..c6b58b3 --- /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 0000000..ee684f2 --- /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 0000000..ab1b090 --- /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 0000000..fefd0e9 --- /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 0000000..62a5464 --- /dev/null +++ b/tests/distances/make_ge_gdml.py @@ -0,0 +1,62 @@ +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", 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) + + +# hpge strings +string_radius = 85 +string_angles = [0, 90, 180, 270] + +n = 0 +for i, det in enumerate(["V", "P", "B", "C"]): + angle = string_angles[i] + n = add_hpge(wl, reg, angle, string_radius, -20, n, det) + +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 0000000..4b7102e --- /dev/null +++ b/tests/distances/test_ge_distance.py @@ -0,0 +1,142 @@ +# 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 legendhpges import draw +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-9): + 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 + 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