From f17265e453ccbbf3aecd409e24237c36d905084e Mon Sep 17 00:00:00 2001 From: Matei Climescu Date: Mon, 13 Jul 2026 22:30:55 +0200 Subject: [PATCH 1/5] fix(decayvolume): derive helium from the SBT surfaces instead of a hand-written inset MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The helium decay region was sized by an inset expression in DecayVolumeFactory that duplicated, by hand, the placement rules used by SBTSensorBuilder and SBTStructureBuilder. The two were not linked, and the copy was wrong in two ways: * It assumed the sensors' inner surface was linear in Z. It is not: placeSideContainer gives each container a flat outer face over the first zSplitOffset() of every sub-frustum so as to clear the columns, so the free region is a sawtooth. The helium overlapped the side containers by 122.6 mm. * It assumed the innermost material in Y was scintillator. It is not: the top/bottom longitudinal beams straddle the sensor plane (the web is omitted because it would pass through the cells) and their inner flange hangs below it, into the decay region. The helium overlapped those flanges by 2.5 mm. Neither was detectable: OverlapCheck is a subsystem-level bounding-box test and does not probe siblings inside decay_volume. Stop maintaining a second model of the geometry. The placement rules now live on SBTConfig as primitives (sideContainerCentreX, topBottomContainerCentreY, longBeamCentreY, zSplitOffset, xHalfAt/yHalfAt); both builders place with them, and the new SBTEnvelope derives the inner free surface from the same expressions. A change to a placement rule now moves the helium automatically. The helium becomes a stack of GeoTraps (2 per sub-frustum, 20 total) so it can follow the sawtooth exactly, and is inset from the true inner surface by helium_clearance_mm, now 1 um — enough to avoid coincident boundaries, which Geant4's navigator handles badly, and far below any engineering tolerance. Verification is now a property, not a formula: SAT (exact for the convex hexahedra involved) between every helium slab and every SBT volume in the built tree, asserting both no protrusion and no margin. A parameter sweep repeats it across 14 SBT configurations, so the guarantee survives re-parameterisation rather than holding at one lucky point. Fiducial helium volume: 534.80 -> 485.85 m^3. The old figure included steel. Closes #58 --- .../include/DecayVolume/DecayVolumeFactory.h | 15 +- .../include/DecayVolume/SBTConfig.h | 95 ++++- .../include/DecayVolume/SBTEnvelope.h | 105 +++++ subsystems/DecayVolume/sbt.toml | 28 +- .../DecayVolume/src/DecayVolumeFactory.cpp | 38 +- subsystems/DecayVolume/src/SBTEnvelope.cpp | 175 ++++++++ .../DecayVolume/src/SBTSensorBuilder.cpp | 57 +-- .../DecayVolume/src/SBTStructureBuilder.cpp | 36 +- subsystems/DecayVolume/test_decayvolume.cpp | 397 +++++++++++++++++- 9 files changed, 855 insertions(+), 91 deletions(-) create mode 100644 subsystems/DecayVolume/include/DecayVolume/SBTEnvelope.h create mode 100644 subsystems/DecayVolume/src/SBTEnvelope.cpp diff --git a/subsystems/DecayVolume/include/DecayVolume/DecayVolumeFactory.h b/subsystems/DecayVolume/include/DecayVolume/DecayVolumeFactory.h index 91a1ee6..eb8a7fe 100644 --- a/subsystems/DecayVolume/include/DecayVolume/DecayVolumeFactory.h +++ b/subsystems/DecayVolume/include/DecayVolume/DecayVolumeFactory.h @@ -3,6 +3,8 @@ #pragma once +#include "DecayVolume/SBTConfig.h" + #include class GeoPhysVol; @@ -19,8 +21,9 @@ class SHiPMaterials; * 50 m rectangular frustum — wrapped around a central helium decay volume. * The SBT geometry is driven by sbt.toml (parsed into an SBTConfig). * - * The helium frustum is sized strictly inside the innermost sensor faces, so - * it does not overlap the structure or the sensors. + * The helium is not an independent volume: it is derived from where the SBT + * material actually is (see SBTEnvelope.h), so that it can neither overlap the + * structure and sensors nor leave an unphysical margin behind. * * Z: 32.92 to 83.32 m -> centre 58.12 m; placement handled by SHiPGeometry. */ @@ -32,9 +35,17 @@ class DecayVolumeFactory { /// Build the DecayVolume geometry; returns the air container. [[nodiscard]] GeoPhysVol* build(); + /// The config the last build() actually used. + /// + /// Tests must reason about *this*, not a default-constructed SBTConfig: + /// the geometry comes from sbt.toml, and a test that checks a clearance + /// against the C++ default is validating a config it did not build. + [[nodiscard]] const SBTConfig& config() const { return m_config; } + private: SHiPMaterials& m_materials; std::string m_configPath; + SBTConfig m_config; // Air container enclosing the SBT structure + sensors and helium (mm). static constexpr double s_halfX = 2200.0; diff --git a/subsystems/DecayVolume/include/DecayVolume/SBTConfig.h b/subsystems/DecayVolume/include/DecayVolume/SBTConfig.h index 490ab30..b182623 100644 --- a/subsystems/DecayVolume/include/DecayVolume/SBTConfig.h +++ b/subsystems/DecayVolume/include/DecayVolume/SBTConfig.h @@ -3,6 +3,7 @@ #pragma once +#include #include namespace SHiPGeometry { @@ -38,7 +39,11 @@ struct SBTConfig { double sensor_clearance_mm = 1.0; // ── Helium decay region ───────────────────────────────────────────── - double helium_clearance_mm = 10.0; + // Gap left between the helium and the innermost SBT material, along the + // coordinate axis. Not a safety margin — just enough to avoid coincident + // surfaces, which Geant4's navigator handles badly. + // NOTE: this default must track sbt.toml. + double helium_clearance_mm = 0.001; // ── Material names (resolved from SHiPMaterials) ──────────────────── std::string material_steel = "Iron"; @@ -54,6 +59,94 @@ struct SBTConfig { double webHeight() const { return hbeam_height_mm - 2.0 * hbeam_flange_thickness_mm; } /// Number of aluminium walls per container (n_cells + 1). int nWalls() const { return n_cells + 1; } + + // ── Frustum profile ───────────────────────────────────────────────── + /// Exit-face Z in the DecayVolume local frame (mm). + double zExit() const { return z_entrance_mm + total_length_mm; } + /// dx_half/dz of the frustum (dimensionless). + double xGrowth() const { return (x_half_exit_mm - x_half_entrance_mm) / total_length_mm; } + /// dy_half/dz of the frustum (dimensionless). + double yGrowth() const { return (y_half_exit_mm - y_half_entrance_mm) / total_length_mm; } + /// Half-extent of the frustum envelope in X at a given Z (mm). + double xHalfAt(double z_mm) const { + return x_half_entrance_mm + (z_mm - z_entrance_mm) * xGrowth(); + } + /// Half-extent of the frustum envelope in Y at a given Z (mm). + double yHalfAt(double z_mm) const { + return y_half_entrance_mm + (z_mm - z_entrance_mm) * yGrowth(); + } + /// Z of the start of sub-frustum @p s (mm). + double zSubLo(int s) const { return z_entrance_mm + s * subLength(); } + + // ── H-beam cross-section primitives ───────────────────────────────── + /// Offset of a flange's mid-plane from the beam axis (mm). + double hbeamFlangeOffset() const { + return 0.5 * hbeam_height_mm - 0.5 * hbeam_flange_thickness_mm; + } + /// Reach of a flange's *outer* surface from the beam axis (mm) = h/2. + double hbeamHalfHeight() const { return 0.5 * hbeam_height_mm; } + + // ══════════════════════════════════════════════════════════════════════ + // PLACEMENT PRIMITIVES — the single source of truth. + // + // These are THE definitions of where SBT material sits. Both builders + // place volumes with them, and SBTEnvelope derives the helium inner + // envelope from them, so a change to any rule below propagates to the + // helium automatically and the two can never drift apart. + // Do not open-code these expressions anywhere else. + // ══════════════════════════════════════════════════════════════════════ + + /// Z offset, from the start of a sub-frustum, of the column front-flange + /// outer edge. Sensor containers are split here: the piece upstream of it + /// must present a *flat* outer face, or it would eat into the column. + double zSplitOffset() const { return 0.5 * hbeam_height_mm + 0.5 * hbeam_flange_thickness_mm; } + + // --- Side (±X) scintillator containers ------------------------------- + /// Half-thickness of a side container in X (mm). + double sideContainerHalfThickness() const { return 0.5 * container_thickness_mm; } + /// |X| of a side container's centroid, given the local frustum half-width. + double sideContainerCentreX(double x_half_mm) const { + return x_half_mm - 0.5 * hbeam_flange_width_mm - sideContainerHalfThickness(); + } + /// |X| of a side container's innermost face (mm). + double sideSensorInnerX(double x_half_mm) const { + return sideContainerCentreX(x_half_mm) - sideContainerHalfThickness(); + } + + // --- Top/bottom (±Y) scintillator containers ------------------------- + /// Half-thickness of a top/bottom container in Y (mm). + double topBottomContainerHalfThickness() const { + return 0.5 * container_thickness_mm - sensor_clearance_mm; + } + /// |Y| of a top/bottom container's centroid (mm). + double topBottomContainerCentreY(double y_half_mm) const { + return y_half_mm - 0.5 * container_thickness_mm; + } + /// |Y| of a top/bottom container's innermost face (mm). + double topBottomSensorInnerY(double y_half_mm) const { + return topBottomContainerCentreY(y_half_mm) - topBottomContainerHalfThickness(); + } + /// Half-extent in X available to top/bottom containers (mm). + double topBottomAvailX(double x_half_mm) const { + return x_half_mm - 0.5 * hbeam_flange_width_mm - 1.0; + } + + // --- Top/bottom longitudinal beams ----------------------------------- + // NOTE: these beams *straddle* the scintillator plane — the outer flange + // sits above it, the web is omitted (it would pass through the cells), + // and the INNER FLANGE HANGS BELOW IT, INSIDE THE DECAY REGION. It, not + // the scintillator, is the innermost material in ±Y. + /// |Y| of a top/bottom longitudinal beam's axis (mm). + double longBeamCentreY(double y_half_mm) const { return y_half_mm - 0.5 * webHeight(); } + /// |Y| reached by a longitudinal beam's inner flange surface (mm). + /// + /// The beam is inclined by the frustum taper, so its cross-section is + /// rotated: the flange surface lies hbeamHalfHeight()/cos(atan(yGrowth)) + /// from the axis measured in world Y, not hbeamHalfHeight(). + double longBeamInnerY(double y_half_mm) const { + const double g = yGrowth(); + return longBeamCentreY(y_half_mm) - hbeamHalfHeight() * std::sqrt(1.0 + g * g); + } }; /** diff --git a/subsystems/DecayVolume/include/DecayVolume/SBTEnvelope.h b/subsystems/DecayVolume/include/DecayVolume/SBTEnvelope.h new file mode 100644 index 0000000..e4f289d --- /dev/null +++ b/subsystems/DecayVolume/include/DecayVolume/SBTEnvelope.h @@ -0,0 +1,105 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) CERN for the benefit of the SHiP Collaboration + +#pragma once + +#include + +class GeoMaterial; +class GeoPhysVol; + +namespace SHiPGeometry { + +struct SBTConfig; + +/** + * @brief The innermost free region of the SBT, and the helium that fills it. + * + * The helium decay region must fill the space inside the SBT exactly: it may + * not protrude into any steel or scintillator, and it may not leave an + * arbitrary safety margin behind either. That makes its size a *consequence* + * of where the SBT material is, never an independent parameter. + * + * This header therefore derives the helium purely from the placement + * primitives on SBTConfig — the same ones SBTStructureBuilder and + * SBTSensorBuilder place their volumes with. Change a placement rule and the + * helium follows automatically; there is no second copy of the arithmetic to + * forget to update. + * + * Two properties of the SBT make this non-trivial, and both are handled here: + * + * - **The X envelope is not linear in Z.** Side containers are split at the + * column front-flange edge (SBTConfig::zSplitOffset()) and the upstream + * piece has a *flat* outer face frozen at the sub-frustum's entrance + * half-width, so it does not clash with the column. The inner surface is + * therefore a sawtooth, dipping inward by up to xGrowth()*zSplitOffset() + * relative to the frustum. A single linear GeoTrap cannot follow it. + * + * - **The innermost material in Y is steel, not scintillator.** The top and + * bottom longitudinal beams straddle the scintillator plane and their + * inner flange hangs *below* it, into the decay region. It sits deeper + * than the cells by (hbeam_height - flange_thickness) - container_thickness + * + sensor_clearance, and it — not the scintillator — bounds the helium. + * + * Consequently the helium is built as a stack of GeoTraps, two per + * sub-frustum, each exactly filling the envelope over its Z span. + */ + +/// One Z slab of the helium: a GeoTrap with rectangular faces at z_lo/z_hi. +struct HeliumPiece { + double z_lo_mm = 0.0; + double z_hi_mm = 0.0; + double dx_lo_mm = 0.0; ///< half-width in X at z_lo + double dx_hi_mm = 0.0; ///< half-width in X at z_hi + double dy_lo_mm = 0.0; ///< half-height in Y at z_lo + double dy_hi_mm = 0.0; ///< half-height in Y at z_hi +}; + +/** + * @brief |X| of the innermost SBT material at @p z_mm (mm). + * + * Minimum over every volume class that can reach the decay region in X. + * Currently only the side scintillator containers do; the columns and corner + * beams sit a further half-flange-width outboard. + */ +double innerFreeHalfX(const SBTConfig& cfg, double z_mm); + +/** + * @brief |Y| of the innermost SBT material at @p z_mm (mm). + * + * Minimum over every volume class that can reach the decay region in Y: the + * top/bottom scintillator containers *and* the inner flange of the top/bottom + * longitudinal beams, which is the binding one. + */ +double innerFreeHalfY(const SBTConfig& cfg, double z_mm); + +/** + * @brief Z values at which innerFreeHalfX/Y change slope. + * + * innerFreeHalfX/Y are piecewise linear with knots exactly here, so a linear + * shape that respects the envelope at these Z values respects it everywhere. + * The helium slab boundaries are drawn from this list. + */ +std::vector envelopeKnots(const SBTConfig& cfg); + +/** + * @brief The helium slabs filling the SBT interior. + * + * Each slab is inset from innerFreeHalfX/Y by cfg.helium_clearance_mm (0 by + * default: the helium touches the SBT and leaves no margin). + * + * @throws std::runtime_error if any slab would have a non-positive half-width, + * i.e. the configured beams/containers have swallowed the decay region. + */ +std::vector heliumPieces(const SBTConfig& cfg); + +/** + * @brief Place the helium slabs into @p container. + * + * Mirrors SBTStructureBuilder / SBTSensorBuilder: takes an SBTConfig and nothing + * else, so the geometry can be built (and swept) without going through sbt.toml. + * Call it *after* the two builders — the helium is a consequence of them. + */ +void buildHelium(GeoPhysVol* container, const GeoMaterial* helium, const SBTConfig& cfg); + +} // namespace SHiPGeometry diff --git a/subsystems/DecayVolume/sbt.toml b/subsystems/DecayVolume/sbt.toml index 8e0ff6a..7aaf77e 100644 --- a/subsystems/DecayVolume/sbt.toml +++ b/subsystems/DecayVolume/sbt.toml @@ -44,10 +44,30 @@ n_cells = 6 sensor_clearance_mm = 1.0 # ── Helium decay region ───────────────────────────────────────────────── -# The central helium volume is a frustum sized strictly inside the innermost -# sensor faces (x_half - container_thickness, y_half - container_thickness) -# minus this clearance, so it cannot overlap the structure or sensors. -helium_clearance_mm = 10.0 +# The helium is NOT sized from a parameter here. It is derived at build time +# from where the SBT material actually is (see SBTEnvelope.h), so that it can +# never overlap the structure or the sensors and never leaves an unphysical +# margin behind. Two consequences are worth knowing about: +# +# * In X it is bounded by the side scintillator containers, whose outer face +# is FLAT over the first zSplitOffset() = (hbeam_height + flange_thickness)/2 +# of each sub-frustum so as to clear the columns. The free region is +# therefore a sawtooth, and the helium is built as one GeoTrap per segment +# (2 per sub-frustum) rather than a single frustum. +# +# * In Y it is bounded by STEEL, not scintillator: the top/bottom +# longitudinal beams straddle the sensor plane and their inner flange hangs +# below it, reaching y_half - 0.5*web_height - 0.5*hbeam_height*sec(taper). +# That is ~13.5 mm deeper than the scintillator's inner face. +# +# This is the gap left between the helium and the SBT surfaces, measured along +# the coordinate axis. It is NOT a safety margin: 1 um is far below any +# engineering tolerance and costs ~0.0001% of the fiducial volume. Its only job +# is to keep the helium and the SBT from sharing a surface, because coincident +# boundaries are a classic source of stuck tracks in Geant4's navigator. +# Set it to 0.0 only if you want exact face-to-face contact and have checked +# that the navigator copes. +helium_clearance_mm = 0.001 # ── Materials (resolved from the central SHiPMaterials catalogue) ──────── material_steel = "Iron" diff --git a/subsystems/DecayVolume/src/DecayVolumeFactory.cpp b/subsystems/DecayVolume/src/DecayVolumeFactory.cpp index 9b2c2af..8566c32 100644 --- a/subsystems/DecayVolume/src/DecayVolumeFactory.cpp +++ b/subsystems/DecayVolume/src/DecayVolumeFactory.cpp @@ -4,6 +4,7 @@ #include "DecayVolume/DecayVolumeFactory.h" #include "DecayVolume/SBTConfig.h" +#include "DecayVolume/SBTEnvelope.h" #include "DecayVolume/SBTSensorBuilder.h" #include "DecayVolume/SBTStructureBuilder.h" #include "SHiPGeometry/ConfigPath.h" @@ -24,6 +25,7 @@ #include #include #include +#include // Source-tree / install-time fallbacks for sbt.toml, set by CMake. #ifndef SBT_TOML_DEFAULT_PATH @@ -50,7 +52,8 @@ DecayVolumeFactory::DecayVolumeFactory(SHiPMaterials& materials, std::string con : m_materials(materials), m_configPath(std::move(configPath)) {} GeoPhysVol* DecayVolumeFactory::build() { - const SBTConfig cfg = readSBTConfig(resolveTomlPath(m_configPath)); + m_config = readSBTConfig(resolveTomlPath(m_configPath)); + const SBTConfig& cfg = m_config; const GeoMaterial* air = m_materials.requireMaterial(cfg.material_air); const GeoMaterial* steel = m_materials.requireMaterial(cfg.material_steel); @@ -97,27 +100,24 @@ GeoPhysVol* DecayVolumeFactory::build() { auto* containerLog = new GeoLogVol("/SHiP/decay_volume", containerBox, air); auto* container = new GeoPhysVol(containerLog); - // ── Central helium decay region ────────────────────────────────────── - // A frustum sized inside the innermost sensor faces (x_half - thickness, - // y_half - thickness) minus the helium clearance, centred on the SBT. - const double inset = cfg.container_thickness_mm + cfg.helium_clearance_mm; - const double dz = 0.5 * cfg.total_length_mm * mm; - const double dx1 = (cfg.x_half_entrance_mm - inset) * mm; - const double dy1 = (cfg.y_half_entrance_mm - inset) * mm; - const double dx2 = (cfg.x_half_exit_mm - inset) * mm; - const double dy2 = (cfg.y_half_exit_mm - inset) * mm; - - auto* heShape = new GeoTrap(dz, 0.0, 0.0, dy1, dx1, dx1, 0.0, dy2, dx2, dx2, 0.0); - auto* heLog = new GeoLogVol("/SHiP/decay_volume/helium", heShape, helium); - auto* hePhys = new GeoPhysVol(heLog); - container->add(new GeoNameTag("/SHiP/decay_volume/helium")); - container->add(new GeoTransform(GeoTrf::Translate3D(0.0, 0.0, 0.0))); - container->add(hePhys); - - // ── SBT steel structure + LAB sensors, wrapping the helium ─────────── + // ── SBT steel structure + LAB sensors ──────────────────────────────── SBTStructureBuilder::build(container, steel, cfg); SBTSensorBuilder::build(container, alMat, labMat, cfg); + // ── Central helium decay region ────────────────────────────────────── + // The helium is not sized independently — it is *derived* from where the + // SBT material actually is (SBTEnvelope, which reads the same placement + // primitives on SBTConfig that the two builders above place with). It + // fills the interior exactly: no protrusion into steel or scintillator, + // and no arbitrary margin left behind (helium_clearance_mm = 0). + // + // It is a stack of GeoTraps rather than one, because the inner surface is + // not linear in Z: the side containers present a flat outer face over the + // first zSplitOffset() of every sub-frustum so as to clear the columns, so + // the free region is a sawtooth. One GeoTrap per envelope segment tracks + // it exactly; a single frustum could not. + buildHelium(container, helium, cfg); + return container; } diff --git a/subsystems/DecayVolume/src/SBTEnvelope.cpp b/subsystems/DecayVolume/src/SBTEnvelope.cpp new file mode 100644 index 0000000..352637a --- /dev/null +++ b/subsystems/DecayVolume/src/SBTEnvelope.cpp @@ -0,0 +1,175 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) CERN for the benefit of the SHiP Collaboration + +#include "DecayVolume/SBTEnvelope.h" + +#include "DecayVolume/SBTConfig.h" + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace SHiPGeometry { + +namespace { + +// Index of the sub-frustum containing z (clamped at the exit face). +int subFrustumAt(const SBTConfig& cfg, double z_mm) { + const double f = (z_mm - cfg.z_entrance_mm) / cfg.subLength(); + return std::clamp(static_cast(std::floor(f)), 0, cfg.n_sub_frustum - 1); +} + +// The half-width the SIDE containers track at z. Inside the flat piece of a +// sub-frustum they are frozen at that sub-frustum's entrance half-width (see +// SBTSensorBuilder::placeSideContainer), which is what makes the X envelope a +// sawtooth rather than a straight line. +double sideTrackedXHalf(const SBTConfig& cfg, double z_mm) { + const int s = subFrustumAt(cfg, z_mm); + const double zLo = cfg.zSubLo(s); + if (z_mm <= zLo + cfg.zSplitOffset()) + return cfg.xHalfAt(zLo); + return cfg.xHalfAt(z_mm); +} + +} // namespace + +double innerFreeHalfX(const SBTConfig& cfg, double z_mm) { + // Side scintillator containers. (Columns and corner beams are outboard of + // these by construction: their inner face is at x_half - flange_width/2, + // a full container_thickness further out.) + return cfg.sideSensorInnerX(sideTrackedXHalf(cfg, z_mm)); +} + +double innerFreeHalfY(const SBTConfig& cfg, double z_mm) { + const double yHalf = cfg.yHalfAt(z_mm); + // Top/bottom scintillator containers ... + const double sensor = cfg.topBottomSensorInnerY(yHalf); + // ... and the longitudinal beams' inner flange, which hangs below them. + const double beam = cfg.longBeamInnerY(yHalf); + // (Cross-beams sit a full beam-height above y_half and never reach in.) + return std::min(sensor, beam); +} + +std::vector envelopeKnots(const SBTConfig& cfg) { + // The knot layout below assumes each sub-frustum is long enough to contain + // its flat piece. If the beams grow (or n_sub_frustum does) until that + // stops being true, the sensor builder itself already produces an inverted + // second piece, so fail here rather than emit unordered knots and build a + // silently wrong helium. + if (cfg.zSplitOffset() >= cfg.subLength()) { + throw std::runtime_error( + "SBTEnvelope: zSplitOffset (" + std::to_string(cfg.zSplitOffset()) + + " mm) >= subLength (" + std::to_string(cfg.subLength()) + + " mm): the sensor containers' flat piece no longer fits inside a sub-frustum. " + "Reduce n_sub_frustum or hbeam_height."); + } + + std::vector knots; + knots.reserve(2 * static_cast(cfg.n_sub_frustum) + 1); + for (int s = 0; s < cfg.n_sub_frustum; ++s) { + knots.push_back(cfg.zSubLo(s)); // sub-frustum boundary + knots.push_back(cfg.zSubLo(s) + cfg.zSplitOffset()); // end of the flat piece + } + knots.push_back(cfg.zExit()); + return knots; +} + +namespace { + +// The X envelope STEPS at each zSplitOffset knot: the flat piece's inner face +// sits at x_half(zLo)-..., and the tracking piece that follows begins at +// x_half(zSplit)-..., xGrowth*zSplitOffset further out. Sampling the outboard +// branch at the knot would put the helium's leading edge in the same plane as +// the flat piece's inner face — no overlap in the strict sense, but a +// coincident surface, which Geant4's navigator will not thank us for, and +// which silently eats the clearance besides. +// +// So evaluate the envelope as a *closed* set: at a knot, take the smaller of +// the two one-sided limits. The helium is then continuous, strictly inscribed, +// and keeps its full clearance everywhere. +double envelopeAtKnot(const SBTConfig& cfg, double z_mm, bool isX) { + constexpr double kEps = 1e-6; + const double lo = std::max(z_mm - kEps, cfg.z_entrance_mm); + const double hi = std::min(z_mm + kEps, cfg.zExit()); + return isX ? std::min(innerFreeHalfX(cfg, lo), innerFreeHalfX(cfg, hi)) + : std::min(innerFreeHalfY(cfg, lo), innerFreeHalfY(cfg, hi)); +} + +} // namespace + +std::vector heliumPieces(const SBTConfig& cfg) { + if (cfg.helium_clearance_mm < 0.0) { + throw std::runtime_error("SBTEnvelope: helium_clearance_mm is negative (" + + std::to_string(cfg.helium_clearance_mm) + + "), which would build a helium that overlaps the SBT by " + "construction."); + } + + const std::vector knots = envelopeKnots(cfg); + const double clr = cfg.helium_clearance_mm; + + std::vector pieces; + pieces.reserve(knots.size() - 1); + + for (std::size_t i = 0; i + 1 < knots.size(); ++i) { + const double zLo = knots[i]; + const double zHi = knots[i + 1]; + if (zHi - zLo < 1e-9) + continue; + + HeliumPiece p; + p.z_lo_mm = zLo; + p.z_hi_mm = zHi; + p.dx_lo_mm = envelopeAtKnot(cfg, zLo, /*isX=*/true) - clr; + p.dx_hi_mm = envelopeAtKnot(cfg, zHi, /*isX=*/true) - clr; + p.dy_lo_mm = envelopeAtKnot(cfg, zLo, /*isX=*/false) - clr; + p.dy_hi_mm = envelopeAtKnot(cfg, zHi, /*isX=*/false) - clr; + + if (p.dx_lo_mm <= 0.0 || p.dx_hi_mm <= 0.0 || p.dy_lo_mm <= 0.0 || p.dy_hi_mm <= 0.0) { + throw std::runtime_error( + "SBTEnvelope: the configured SBT leaves no decay region between z=" + + std::to_string(zLo) + " and z=" + std::to_string(zHi) + + " mm (helium half-extents " + std::to_string(p.dx_lo_mm) + "/" + + std::to_string(p.dy_lo_mm) + + " mm). The containers, beams or clearances in sbt.toml are too large for the " + "frustum."); + } + pieces.push_back(p); + } + return pieces; +} + +void buildHelium(GeoPhysVol* container, const GeoMaterial* helium, const SBTConfig& cfg) { + using namespace GeoModelKernelUnits; + + const std::vector pieces = heliumPieces(cfg); + for (std::size_t i = 0; i < pieces.size(); ++i) { + const HeliumPiece& p = pieces[i]; + const double dz = 0.5 * (p.z_hi_mm - p.z_lo_mm) * mm; + const double zMid = 0.5 * (p.z_lo_mm + p.z_hi_mm) * mm; + const double dx1 = p.dx_lo_mm * mm; + const double dy1 = p.dy_lo_mm * mm; + const double dx2 = p.dx_hi_mm * mm; + const double dy2 = p.dy_hi_mm * mm; + + const std::string name = "/SHiP/decay_volume/helium_" + std::to_string(i); + auto* shape = new GeoTrap(dz, 0.0, 0.0, dy1, dx1, dx1, 0.0, dy2, dx2, dx2, 0.0); + auto* log = new GeoLogVol(name, shape, helium); + auto* phys = new GeoPhysVol(log); + container->add(new GeoNameTag(name)); + container->add(new GeoTransform(GeoTrf::Translate3D(0.0, 0.0, zMid))); + container->add(phys); + } +} + +} // namespace SHiPGeometry diff --git a/subsystems/DecayVolume/src/SBTSensorBuilder.cpp b/subsystems/DecayVolume/src/SBTSensorBuilder.cpp index 8ad98ef..bad66ee 100644 --- a/subsystems/DecayVolume/src/SBTSensorBuilder.cpp +++ b/subsystems/DecayVolume/src/SBTSensorBuilder.cpp @@ -64,14 +64,11 @@ void placeSideContainer(GeoVPhysVol* mother, const GeoMaterial* alMat, const Geo const SBTConfig& cfg, const std::string& name, double sgnX, double xHalf_lo, double xHalf_hi, double zLo_mm, double zHi_mm, double dy1, double dy2, double yCtr1, double yCtr2) { - const double contThick = cfg.container_thickness_mm; - const double hBeamH = cfg.hbeam_height_mm; - const double hBeamW = cfg.hbeam_flange_width_mm; - const double hBeamTf = cfg.hbeam_flange_thickness_mm; - - const double dx = 0.5 * contThick * mm; + // Geometry rules live on SBTConfig (see "PLACEMENT PRIMITIVES" there) so + // that SBTEnvelope can derive the helium from exactly the same expressions. + const double dx = cfg.sideContainerHalfThickness() * mm; - const double z_split_rel = 0.5 * hBeamH + 0.5 * hBeamTf; // column front-flange outer edge + const double z_split_rel = cfg.zSplitOffset(); // column front-flange outer edge const double z_split_mm = zLo_mm + z_split_rel; // ---- Piece 1: z in [zLo, z_split] — flat outer face, Y shear only ---- @@ -85,7 +82,7 @@ void placeSideContainer(GeoVPhysVol* mother, const GeoMaterial* alMat, const Geo const double yC1 = yCtr1; const double yC2 = yCtr1 + (yCtr2 - yCtr1) * frac1; - const double xCtr_flat = sgnX * (xHalf_lo - 0.5 * hBeamW - 0.5 * contThick) * mm; + const double xCtr_flat = sgnX * cfg.sideContainerCentreX(xHalf_lo) * mm; const double dY1 = yC2 - yC1; const double shear1 = std::abs(dY1); @@ -109,8 +106,8 @@ void placeSideContainer(GeoVPhysVol* mother, const GeoMaterial* alMat, const Geo const double yC2 = yCtr2; const double xHalf_split = xHalf_lo + (xHalf_hi - xHalf_lo) * frac1; - const double xCtr1_p2 = sgnX * (xHalf_split - 0.5 * hBeamW - 0.5 * contThick) * mm; - const double xCtr2_p2 = sgnX * (xHalf_hi - 0.5 * hBeamW - 0.5 * contThick) * mm; + const double xCtr1_p2 = sgnX * cfg.sideContainerCentreX(xHalf_split) * mm; + const double xCtr2_p2 = sgnX * cfg.sideContainerCentreX(xHalf_hi) * mm; const double dX2 = xCtr2_p2 - xCtr1_p2; const double dY2 = yC2 - yC1; @@ -205,34 +202,24 @@ void SBTSensorBuilder::build(GeoVPhysVol* mother, const GeoMaterial* alMat, const GeoMaterial* labMat, const SBTConfig& cfg, const std::string& tag) { // Bind config to the names the ported body uses (magnitudes in mm). - const double xHalf_ent = cfg.x_half_entrance_mm; - const double yHalf_ent = cfg.y_half_entrance_mm; - const double xHalf_ext = cfg.x_half_exit_mm; - const double yHalf_ext = cfg.y_half_exit_mm; - const double zTotal = cfg.total_length_mm; const int nSub = cfg.n_sub_frustum; const double subLen = cfg.subLength(); const double hBeamH = cfg.hbeam_height_mm; - const double hBeamW = cfg.hbeam_flange_width_mm; - const double hBeamTf = cfg.hbeam_flange_thickness_mm; - const double contThick = cfg.container_thickness_mm; const double sensorClear = cfg.sensor_clearance_mm; const double zEntrance_mm = cfg.z_entrance_mm; - auto xHalfAtZ = [&](double z_mm, double zEnt_mm) { - return xHalf_ent + (z_mm - zEnt_mm) / zTotal * (xHalf_ext - xHalf_ent); - }; - auto yHalfAtZ = [&](double z_mm, double zEnt_mm) { - return yHalf_ent + (z_mm - zEnt_mm) / zTotal * (yHalf_ext - yHalf_ent); - }; + // Frustum profile and all placement rules come from SBTConfig, so that + // SBTEnvelope sizes the helium from the very same expressions. + auto xHalfAtZ = [&](double z_mm) { return cfg.xHalfAt(z_mm); }; + auto yHalfAtZ = [&](double z_mm) { return cfg.yHalfAt(z_mm); }; // (A) SIDE CONTAINERS (±X faces) — 4 containers per side, split by Y=0. for (int s = 0; s < nSub; ++s) { const double zLo_mm = zEntrance_mm + s * subLen; const double zHi_mm = zEntrance_mm + (s + 1) * subLen; - const double availLo = (yHalfAtZ(zLo_mm, zEntrance_mm) - hBeamH) * mm; - const double availHi = (yHalfAtZ(zHi_mm, zEntrance_mm) - hBeamH) * mm; + const double availLo = (yHalfAtZ(zLo_mm) - hBeamH) * mm; + const double availHi = (yHalfAtZ(zHi_mm) - hBeamH) * mm; const double dy1[4] = {availLo / 4.0, availLo / 4.0, availLo / 4.0, availLo / 4.0}; const double dy2[4] = {availHi / 4.0, availHi / 4.0, availHi / 4.0, availHi / 4.0}; @@ -249,9 +236,9 @@ void SBTSensorBuilder::build(GeoVPhysVol* mother, const GeoMaterial* alMat, const std::string cname = tag + "_Side_S" + std::to_string(s) + "_X" + (side == 0 ? "P" : "M") + "_C" + std::to_string(ci); - placeSideContainer(mother, alMat, labMat, cfg, cname, sgnX, - xHalfAtZ(zLo_mm, zEntrance_mm), xHalfAtZ(zHi_mm, zEntrance_mm), - zLo_mm, zHi_mm, dy1[ci], dy2[ci], yCtr1[ci], yCtr2[ci]); + placeSideContainer(mother, alMat, labMat, cfg, cname, sgnX, xHalfAtZ(zLo_mm), + xHalfAtZ(zHi_mm), zLo_mm, zHi_mm, dy1[ci], dy2[ci], yCtr1[ci], + yCtr2[ci]); } } } @@ -265,19 +252,19 @@ void SBTSensorBuilder::build(GeoVPhysVol* mother, const GeoMaterial* alMat, const double zLo_mm = zEntrance_mm + s * subLen; const double zHi_mm = zEntrance_mm + (s + 1) * subLen; - const double xAvailLo = (xHalfAtZ(zLo_mm, zEntrance_mm) - 0.5 * hBeamW - 1.0) * mm; - const double xAvailHi = (xHalfAtZ(zHi_mm, zEntrance_mm) - 0.5 * hBeamW - 1.0) * mm; + const double xAvailLo = cfg.topBottomAvailX(xHalfAtZ(zLo_mm)) * mm; + const double xAvailHi = cfg.topBottomAvailX(xHalfAtZ(zHi_mm)) * mm; - const double yFaceTop_Lo = (yHalfAtZ(zLo_mm, zEntrance_mm) - 0.5 * contThick) * mm; - const double yFaceTop_Hi = (yHalfAtZ(zHi_mm, zEntrance_mm) - 0.5 * contThick) * mm; + const double yFaceTop_Lo = cfg.topBottomContainerCentreY(yHalfAtZ(zLo_mm)) * mm; + const double yFaceTop_Hi = cfg.topBottomContainerCentreY(yHalfAtZ(zHi_mm)) * mm; const double yFaceBot_Lo = -yFaceTop_Lo; const double yFaceBot_Hi = -yFaceTop_Hi; - const double dx = (0.5 * contThick - sensorClear) * mm; + const double dx = cfg.topBottomContainerHalfThickness() * mm; const int nCont = (s < threshold) ? 2 : 3; - const double z_split_rel = 0.5 * hBeamH + 0.5 * hBeamTf; // column front-flange outer edge + const double z_split_rel = cfg.zSplitOffset(); // column front-flange outer edge const double z_split_mm = zLo_mm + z_split_rel; const double frac_split = z_split_rel / (zHi_mm - zLo_mm); diff --git a/subsystems/DecayVolume/src/SBTStructureBuilder.cpp b/subsystems/DecayVolume/src/SBTStructureBuilder.cpp index acbe38d..3407bcb 100644 --- a/subsystems/DecayVolume/src/SBTStructureBuilder.cpp +++ b/subsystems/DecayVolume/src/SBTStructureBuilder.cpp @@ -198,29 +198,18 @@ void placeHBeamInclined(GeoVPhysVol* parent, const GeoMaterial* mat, const SBTCo void SBTStructureBuilder::build(GeoVPhysVol* mother, const GeoMaterial* steel, const SBTConfig& cfg, const std::string& tag) { // Bind config values to the names the ported body uses (magnitudes in mm). - const double xHalf_entrance = cfg.x_half_entrance_mm; - const double yHalf_entrance = cfg.y_half_entrance_mm; - const double xHalf_exit = cfg.x_half_exit_mm; - const double yHalf_exit = cfg.y_half_exit_mm; - const double totalLength = cfg.total_length_mm; const int nSubFrustrum = cfg.n_sub_frustum; const double subLength = cfg.subLength(); const double yFloor = cfg.y_floor_mm; const double hBeamH = cfg.hbeam_height_mm; const double hBeamW = cfg.hbeam_flange_width_mm; const double hBeamTf = cfg.hbeam_flange_thickness_mm; - const double hBeamHw = cfg.webHeight(); const double zEntrance_mm = cfg.z_entrance_mm; - // Linear interpolation of the X/Y half-extents at a given Z. - auto xHalfAtZ = [&](double z_mm, double zEnt_mm) { - const double frac = (z_mm - zEnt_mm) / totalLength; - return xHalf_entrance + frac * (xHalf_exit - xHalf_entrance); - }; - auto yHalfAtZ = [&](double z_mm, double zEnt_mm) { - const double frac = (z_mm - zEnt_mm) / totalLength; - return yHalf_entrance + frac * (yHalf_exit - yHalf_entrance); - }; + // Frustum profile from SBTConfig — the same accessors SBTEnvelope uses to + // size the helium, so the two can never disagree. + auto xHalfAtZ = [&](double z_mm, double /*zEnt_mm*/) { return cfg.xHalfAt(z_mm); }; + auto yHalfAtZ = [&](double z_mm, double /*zEnt_mm*/) { return cfg.yHalfAt(z_mm); }; // (A) VERTICAL COLUMNS — 11 rows x 2 sides, frustum top -> floor. for (int row = 0; row <= nSubFrustrum; ++row) { @@ -289,12 +278,15 @@ void SBTStructureBuilder::build(GeoVPhysVol* mother, const GeoMaterial* steel, c const double xLo = xHalfAtZ(zLo_mm, zEntrance_mm); const double xHi = xHalfAtZ(zHi_mm, zEntrance_mm); - const double yBeamOffset = 0.5 * hBeamHw; // centred C-channel - - const double yTop_Lo = +yHalfAtZ(zLo_mm, zEntrance_mm) - yBeamOffset; - const double yTop_Hi = +yHalfAtZ(zHi_mm, zEntrance_mm) - yBeamOffset; - const double yBot_Lo = -yHalfAtZ(zLo_mm, zEntrance_mm) + yBeamOffset; - const double yBot_Hi = -yHalfAtZ(zHi_mm, zEntrance_mm) + yBeamOffset; + // Straddle beam: the outer flange sits above the scintillator, the web + // is omitted (it would pass through the cells) and the INNER FLANGE + // HANGS BELOW THE SCINTILLATOR, into the decay region. That inner + // flange is the innermost material in ±Y and therefore what bounds the + // helium — see SBTConfig::longBeamInnerY() and SBTEnvelope. + const double yTop_Lo = +cfg.longBeamCentreY(yHalfAtZ(zLo_mm, zEntrance_mm)); + const double yTop_Hi = +cfg.longBeamCentreY(yHalfAtZ(zHi_mm, zEntrance_mm)); + const double yBot_Lo = -cfg.longBeamCentreY(yHalfAtZ(zLo_mm, zEntrance_mm)); + const double yBot_Hi = -cfg.longBeamCentreY(yHalfAtZ(zHi_mm, zEntrance_mm)); const std::string sTag = tag + "_SF" + std::to_string(s); @@ -393,7 +385,7 @@ void SBTStructureBuilder::build(GeoVPhysVol* mother, const GeoMaterial* steel, c const std::string rowTag = tag + "_XBeam_R" + std::to_string(row); - const double yGrowthPerZ = (yHalf_exit - yHalf_entrance) / totalLength; + const double yGrowthPerZ = cfg.yGrowth(); const double xbShift = (0.5 * hBeamH + 0.5 * hBeamW * yGrowthPerZ + 5.0) * mm; const double xbHalfLen = (xEdge_mm - 0.5 * hBeamW) * mm; diff --git a/subsystems/DecayVolume/test_decayvolume.cpp b/subsystems/DecayVolume/test_decayvolume.cpp index 07bb89b..a5204f7 100644 --- a/subsystems/DecayVolume/test_decayvolume.cpp +++ b/subsystems/DecayVolume/test_decayvolume.cpp @@ -2,17 +2,27 @@ // Copyright (C) CERN for the benefit of the SHiP Collaboration #include "DecayVolume/DecayVolumeFactory.h" +#include "DecayVolume/SBTConfig.h" +#include "DecayVolume/SBTEnvelope.h" +#include "DecayVolume/SBTSensorBuilder.h" +#include "DecayVolume/SBTStructureBuilder.h" #include "SHiPGeometry/SHiPMaterials.h" #include +#include #include #include #include #include #include +#include #include +#include +#include +#include #include +#include using SHiPGeometry::SHiPMaterials; @@ -72,30 +82,31 @@ TEST_CASE("DecayVolumeStructureBoxCount", "[decayvolume]") { } // Sensor system: 130 containers, each Z-split into 2 pieces of 7 Al walls + -// 6 LAB cells (13 GeoTraps) -> 130*2*13 = 3380 sensor traps; plus the single -// helium frustum = 3381 GeoTrap children. +// 6 LAB cells (13 GeoTraps) -> 130*2*13 = 3380 sensor traps; plus the helium, +// which is 2 slabs per sub-frustum (20) rather than a single frustum, because +// the free region it fills is not linear in Z. 3380 + 20 = 3400 GeoTraps. TEST_CASE("DecayVolumeSensorTrapCount", "[decayvolume]") { SHiPMaterials materials; SHiPGeometry::DecayVolumeFactory factory(materials); GeoPhysVol* dv = factory.build(); REQUIRE(dv != nullptr); const ChildShapeCounts c = countByShape(dv); - CHECK(c.traps == 3381u); // NOLINT(readability/check) - CHECK(c.helium == 1u); // NOLINT(readability/check) + CHECK(c.traps == 3400u); // NOLINT(readability/check) + CHECK(c.helium == 20u); // NOLINT(readability/check) } // Flat architecture: total direct children = 312 structure + 3380 sensors + -// 1 helium = 3693, with no grandchildren. +// 20 helium slabs = 3712, with no grandchildren. TEST_CASE("DecayVolumeChildCount", "[decayvolume]") { SHiPMaterials materials; SHiPGeometry::DecayVolumeFactory factory(materials); GeoPhysVol* dv = factory.build(); REQUIRE(dv != nullptr); - CHECK(dv->getNChildVols() == 3693u); // NOLINT(readability/check) + CHECK(dv->getNChildVols() == 3712u); // NOLINT(readability/check) } -// The central decay region is a helium GeoTrap (frustum), sized inside the -// innermost sensor faces so it cannot overlap the structure or sensors. +// The central decay region is built from helium GeoTraps, derived from the +// innermost SBT surfaces so that they cannot overlap the structure or sensors. TEST_CASE("DecayVolumeHasHeliumFrustum", "[decayvolume]") { SHiPMaterials materials; SHiPGeometry::DecayVolumeFactory factory(materials); @@ -113,3 +124,373 @@ TEST_CASE("DecayVolumeHasHeliumFrustum", "[decayvolume]") { auto* trap = dynamic_cast(he->getLogVol()->getShape()); REQUIRE(trap != nullptr); } + +// ───────────────────────────────────────────────────────────────────────────── +// Helium envelope regression tests. + +// +// The point of these tests is that they assert a *property*, not a formula. +// They do not check that dx1 is 632.376 mm; they walk the geometry that the +// builders actually produced and check that no helium slab intersects any of +// it, and that the helium is flush against it. Re-parameterise the SBT — move +// a beam, change the flange width, resize the containers — and the tests still +// mean the right thing. That is the property PR #58's algebraic check lacked: +// it verified an expression, and the expression was an incomplete model of the +// geometry (it missed the flat-piece sawtooth in X, and the longitudinal beam +// flanges entirely). +// +// Overlap is decided by the separating-axis theorem. Every child of the decay +// volume is a GeoBox or a GeoTrap, i.e. a convex hexahedron, so SAT over the +// two shapes' face normals plus their pairwise edge cross-products is exact. + +namespace { + +using Vec = GeoTrf::Vector3D; + +// A convex hexahedron: 8 world-space vertices, in the GeoTrap corner order +// (0-3 at -dz, 4-7 at +dz; within a face: -y-x, -y+x, +y-x, +y+x). +struct Hexa { + std::string name; + std::array v; +}; + +// The 6 quad faces and 12 edges of that corner ordering. +constexpr std::array, 6> kFaces = { + {{0, 1, 3}, {4, 5, 7}, {0, 1, 5}, {2, 3, 7}, {0, 2, 6}, {1, 3, 7}}}; +constexpr std::array, 12> kEdges = {{{0, 1}, + {1, 3}, + {3, 2}, + {2, 0}, + {4, 5}, + {5, 7}, + {7, 6}, + {6, 4}, + {0, 4}, + {1, 5}, + {2, 6}, + {3, 7}}}; + +std::array boxVertices(const GeoBox& b) { + const double hx = b.getXHalfLength(), hy = b.getYHalfLength(), hz = b.getZHalfLength(); + std::array v; + int i = 0; + for (int sz : {-1, 1}) + for (int sy : {-1, 1}) + for (int sx : {-1, 1}) + v[i++] = Vec(sx * hx, sy * hy, sz * hz); + return v; +} + +std::array trapVertices(const GeoTrap& t) { + const double dz = t.getZHalfLength(); + const double tt = std::tan(t.getTheta()); + const double cx = tt * std::cos(t.getPhi()); + const double cy = tt * std::sin(t.getPhi()); + + const std::array dy = {t.getDydzn(), t.getDydzp()}; + const std::array dxn = {t.getDxdyndzn(), t.getDxdyndzp()}; + const std::array dxp = {t.getDxdypdzn(), t.getDxdypdzp()}; + const std::array alp = {t.getAngleydzn(), t.getAngleydzp()}; + + std::array v; + int i = 0; + for (int f = 0; f < 2; ++f) { + const double s = (f == 0) ? -1.0 : +1.0; + const double ox = s * dz * cx, oy = s * dz * cy, oz = s * dz; + const double ta = std::tan(alp[f]); + v[i++] = Vec(ox - dy[f] * ta - dxn[f], oy - dy[f], oz); + v[i++] = Vec(ox - dy[f] * ta + dxn[f], oy - dy[f], oz); + v[i++] = Vec(ox + dy[f] * ta - dxp[f], oy + dy[f], oz); + v[i++] = Vec(ox + dy[f] * ta + dxp[f], oy + dy[f], oz); + } + return v; +} + +// Separation of two convex hexahedra along the SAT axis set. +// > 0 => disjoint (and the value is a lower bound on their distance) +// <= 0 => they intersect +double separation(const Hexa& a, const Hexa& b) { + double best = -std::numeric_limits::infinity(); + + auto probe = [&](const Vec& n) { + const double nn = n.norm(); + if (nn < 1e-9) + return; + const Vec u = n / nn; + double amin = std::numeric_limits::infinity(), amax = -amin; + double bmin = amin, bmax = amax; + for (const Vec& p : a.v) { + const double d = p.dot(u); + amin = std::min(amin, d); + amax = std::max(amax, d); + } + for (const Vec& p : b.v) { + const double d = p.dot(u); + bmin = std::min(bmin, d); + bmax = std::max(bmax, d); + } + best = std::max(best, std::max(bmin - amax, amin - bmax)); + }; + + for (const Hexa* h : {&a, &b}) + for (const auto& f : kFaces) + probe((h->v[f[1]] - h->v[f[0]]).cross(h->v[f[2]] - h->v[f[0]])); + + for (const auto& ea : kEdges) + for (const auto& eb : kEdges) + probe((a.v[ea[1]] - a.v[ea[0]]).cross(b.v[eb[1]] - b.v[eb[0]])); + + return best; +} + +// Collect every direct child of the decay volume as a world-space hexahedron. +// The SBT is placed flat, so all children are leaves. +void collect(const GeoVPhysVol* dv, std::vector& helium, std::vector& sbt) { + for (unsigned int i = 0; i < dv->getNChildVols(); ++i) { + const GeoVPhysVol* child = &*dv->getChildVol(i); + const GeoLogVol* lv = child->getLogVol(); + const GeoShape* shape = lv->getShape(); + + std::array local; + if (const auto* b = dynamic_cast(shape)) + local = boxVertices(*b); + else if (const auto* t = dynamic_cast(shape)) + local = trapVertices(*t); + else + continue; // no other shape types are placed + + const GeoTrf::Transform3D x = dv->getXToChildVol(i); + Hexa h; + h.name = lv->getName(); + for (int k = 0; k < 8; ++k) + h.v[k] = x * local[k]; + + (h.name.find("helium") != std::string::npos ? helium : sbt).push_back(h); + } +} + +struct Built { + GeoPhysVol* dv = nullptr; + SHiPGeometry::SBTConfig cfg; // the config the geometry was ACTUALLY built from + std::vector helium, sbt; +}; + +// Build via the factory, i.e. from sbt.toml. Carries the resolved config back +// out, so the tests never compare the built geometry against a default- +// constructed SBTConfig that may say something different. +Built buildFromToml() { + static SHiPMaterials materials; + SHiPGeometry::DecayVolumeFactory factory(materials); + Built b; + b.dv = factory.build(); + REQUIRE(b.dv != nullptr); + b.cfg = factory.config(); + collect(b.dv, b.helium, b.sbt); + return b; +} + +// Build the SBT + helium directly from an arbitrary SBTConfig, bypassing the +// toml. This is what lets us sweep the parameter space. +Built buildFromConfig(const SHiPGeometry::SBTConfig& cfg) { + static SHiPMaterials materials; + const GeoMaterial* air = materials.requireMaterial(cfg.material_air); + const GeoMaterial* steel = materials.requireMaterial(cfg.material_steel); + const GeoMaterial* alMat = materials.requireMaterial(cfg.material_wall); + const GeoMaterial* labMat = materials.requireMaterial(cfg.material_cell); + const GeoMaterial* helium = materials.requireMaterial(cfg.material_helium); + + // Generous container: this test cares about helium-vs-SBT, not the envelope. + auto* boxShape = new GeoBox(10000.0, 10000.0, 40000.0); + auto* boxLog = new GeoLogVol("/SHiP/test_container", boxShape, air); + auto* container = new GeoPhysVol(boxLog); + + SHiPGeometry::SBTStructureBuilder::build(container, steel, cfg); + SHiPGeometry::SBTSensorBuilder::build(container, alMat, labMat, cfg); + SHiPGeometry::buildHelium(container, helium, cfg); + + Built b; + b.dv = container; + b.cfg = cfg; + collect(container, b.helium, b.sbt); + return b; +} + +// Geometric tolerance for the SAT assertions. Must stay well below +// helium_clearance_mm (1 um), or the clearance checks become vacuous; and well +// above double-precision noise on ~1e4 mm coordinates (~1e-8 mm). +constexpr double kTol = 1e-6; + +// helium_clearance_mm is a gap measured along a coordinate axis. SAT returns a +// Euclidean distance, and the surfaces bounding the helium are tilted by the +// frustum taper, so an axis gap of c shows up as c*cos(tilt). Assert the band. +double minExpectedSeparation(const SHiPGeometry::SBTConfig& cfg) { + const double g = std::max(std::abs(cfg.xGrowth()), std::abs(cfg.yGrowth())); + return cfg.helium_clearance_mm / std::sqrt(1.0 + g * g); +} + +// Closest approach between any helium slab and any SBT volume. +// > 0 disjoint, 0 touching, < 0 overlapping. +double closestApproach(const Built& b, std::string* culprit = nullptr) { + double worst = std::numeric_limits::infinity(); + for (const Hexa& he : b.helium) { + double zlo = std::numeric_limits::infinity(), zhi = -zlo; + for (const Vec& p : he.v) { + zlo = std::min(zlo, p.z()); + zhi = std::max(zhi, p.z()); + } + for (const Hexa& o : b.sbt) { + double ozlo = std::numeric_limits::infinity(), ozhi = -ozlo; + for (const Vec& p : o.v) { + ozlo = std::min(ozlo, p.z()); + ozhi = std::max(ozhi, p.z()); + } + if (ozhi < zlo - 1.0 || ozlo > zhi + 1.0) + continue; // cheap Z reject; a disjoint pair cannot be the minimum + const double s = separation(he, o); + if (s < worst) { + worst = s; + if (culprit) + *culprit = he.name + " vs " + o.name; + } + } + } + return worst; +} + +} // namespace + +// The whole point. No helium slab may intersect any SBT volume — not the +// scintillator containers, not the columns, not the corner beams, and (the one +// PR #58 missed) not the inner flanges of the top/bottom longitudinal beams, +// which hang below the sensor plane into the decay region. +TEST_CASE("HeliumDoesNotOverlapAnySBTVolume", "[decayvolume][envelope]") { + const Built b = buildFromToml(); + REQUIRE(!b.helium.empty()); + REQUIRE(b.sbt.size() > 100); + + std::string culprit; + const double worst = closestApproach(b, &culprit); + + INFO("closest approach: " << worst << " mm, between " << culprit); + CHECK(worst >= -kTol); // NOLINT(readability/check) +} + +// ... and no unphysical margin either: with helium_clearance_mm = 0 the helium +// must actually touch the material that bounds it. If a future change to the +// SBT introduced a volume that SBTEnvelope does not know about, the test above +// would fail; if SBTEnvelope became over-conservative, this one would. +TEST_CASE("HeliumIsFlushWithTheSBT", "[decayvolume][envelope]") { + const Built b = buildFromToml(); + const double worst = closestApproach(b); + + INFO("closest approach: " << worst << " mm; want [" << minExpectedSeparation(b.cfg) << ", " + << b.cfg.helium_clearance_mm << "]"); + CHECK(worst >= minExpectedSeparation(b.cfg) - kTol); // NOLINT(readability/check) no gouging + CHECK(worst <= b.cfg.helium_clearance_mm + kTol); // NOLINT(readability/check) no margin +} + +// The helium fills the analytic envelope exactly, sampled densely rather than +// only at the slab boundaries — this catches an envelope whose knots are in the +// wrong places (e.g. if zSplitOffset() changed but envelopeKnots() did not). +TEST_CASE("HeliumMatchesAnalyticEnvelope", "[decayvolume][envelope]") { + const SHiPGeometry::SBTConfig cfg = buildFromToml().cfg; + const auto pieces = SHiPGeometry::heliumPieces(cfg); + REQUIRE(pieces.size() == 2u * static_cast(cfg.n_sub_frustum)); + + for (const auto& p : pieces) { + for (int k = 0; k <= 32; ++k) { + const double t = static_cast(k) / 32.0; + const double z = p.z_lo_mm + t * (p.z_hi_mm - p.z_lo_mm); + const double dx = p.dx_lo_mm + t * (p.dx_hi_mm - p.dx_lo_mm); + const double dy = p.dy_lo_mm + t * (p.dy_hi_mm - p.dy_lo_mm); + + // Sample strictly inside the slab so the flat/tracking branch of + // the X envelope is evaluated on the right side of a knot. + const double zs = std::min(std::max(z, p.z_lo_mm + 1e-6), p.z_hi_mm - 1e-6); + const double freeX = SHiPGeometry::innerFreeHalfX(cfg, zs); + const double freeY = SHiPGeometry::innerFreeHalfY(cfg, zs); + + CHECK(dx <= freeX - cfg.helium_clearance_mm + kTol); // NOLINT(readability/check) + CHECK(dy <= freeY - cfg.helium_clearance_mm + kTol); // NOLINT(readability/check) + } + } +} + +// THE test for "does this survive changes to the SBT?". A single configuration +// proves nothing about that — it only shows the arithmetic is right at one +// point. So: perturb each parameter the SBT is actually likely to be +// re-specified with, rebuild the structure, the sensors AND the helium from +// scratch, and re-run the overlap check. Every case must still come out flush. +// +// If a future change to either builder breaks the envelope's model of it, this +// fails across the board rather than at one lucky configuration. +TEST_CASE("HeliumIsFlushAcrossTheParameterSpace", "[decayvolume][envelope][sweep]") { + const SHiPGeometry::SBTConfig base = buildFromToml().cfg; + + struct Variation { + const char* what; + std::function apply; + }; + + const auto variations = std::vector{ + {"baseline", [](SHiPGeometry::SBTConfig&) {}}, + {"steeper X taper", [](SHiPGeometry::SBTConfig& c) { c.x_half_exit_mm = 3000.0; }}, + {"steeper Y taper", [](SHiPGeometry::SBTConfig& c) { c.y_half_exit_mm = 4500.0; }}, + {"no taper at all", + [](SHiPGeometry::SBTConfig& c) { + c.x_half_exit_mm = c.x_half_entrance_mm; + c.y_half_exit_mm = c.y_half_entrance_mm; + }}, + {"wider flange", [](SHiPGeometry::SBTConfig& c) { c.hbeam_flange_width_mm = 400.0; }}, + {"taller beam", [](SHiPGeometry::SBTConfig& c) { c.hbeam_height_mm = 400.0; }}, + {"thicker flange", [](SHiPGeometry::SBTConfig& c) { c.hbeam_flange_thickness_mm = 30.0; }}, + {"thinner containers", + [](SHiPGeometry::SBTConfig& c) { c.container_thickness_mm = 120.0; }}, + {"thicker containers", + [](SHiPGeometry::SBTConfig& c) { c.container_thickness_mm = 300.0; }}, + {"more sub-frusta", [](SHiPGeometry::SBTConfig& c) { c.n_sub_frustum = 20; }}, + {"fewer sub-frusta", [](SHiPGeometry::SBTConfig& c) { c.n_sub_frustum = 5; }}, + {"bigger sensor clearance", + [](SHiPGeometry::SBTConfig& c) { c.sensor_clearance_mm = 5.0; }}, + {"non-zero helium clearance", + [](SHiPGeometry::SBTConfig& c) { c.helium_clearance_mm = 10.0; }}, + {"shorter SBT", [](SHiPGeometry::SBTConfig& c) { c.total_length_mm = 20000.0; }}, + }; + + for (const Variation& v : variations) { + SHiPGeometry::SBTConfig cfg = base; + v.apply(cfg); + + const Built b = buildFromConfig(cfg); + std::string culprit; + const double worst = closestApproach(b, &culprit); + + INFO("variation: " << v.what << " -> closest approach " << worst << " mm; want [" + << minExpectedSeparation(cfg) << ", " << cfg.helium_clearance_mm + << "], nearest " << culprit); + CHECK(worst >= -kTol); // NOLINT(readability/check) no overlap + CHECK(worst >= + minExpectedSeparation(cfg) - kTol); // NOLINT(readability/check) clearance kept + CHECK(worst <= cfg.helium_clearance_mm + kTol); // NOLINT(readability/check) no margin + } +} + +// Guard rail: an SBT whose beams and containers have eaten the whole frustum +// must fail loudly, not silently produce an inverted GeoTrap. +TEST_CASE("HeliumRejectsAnImpossibleSBT", "[decayvolume][envelope]") { + SECTION("containers larger than the frustum") { + SHiPGeometry::SBTConfig cfg; + cfg.container_thickness_mm = 5000.0; + CHECK_THROWS(SHiPGeometry::heliumPieces(cfg)); + } + SECTION("sub-frustum shorter than the sensor flat piece") { + SHiPGeometry::SBTConfig cfg; + cfg.n_sub_frustum = 500; // subLength 100 mm < zSplitOffset 131.25 mm + CHECK_THROWS(SHiPGeometry::heliumPieces(cfg)); + } + SECTION("negative clearance would overlap by construction") { + SHiPGeometry::SBTConfig cfg; + cfg.helium_clearance_mm = -5.0; + CHECK_THROWS(SHiPGeometry::heliumPieces(cfg)); + } +} From 6f31a447f18ec260c5b795ef2b4c2516fd51cc8b Mon Sep 17 00:00:00 2001 From: Matei Climescu Date: Tue, 14 Jul 2026 00:41:38 +0200 Subject: [PATCH 2/5] fix(decayvolume): reject sensor_clearance_mm that inverts the container readSBTConfig required sensor_clearance_mm > 0 but never bounded it from above. It is trimmed off each side of a top/bottom container, so any value >= half the container thickness makes SBTConfig::topBottomContainerHalfThickness() negative. That is worse than a degenerate GeoTrap. topBottomSensorInnerY() subtracts the half-thickness, so a negative one reports the scintillator's inner face as further out than it really is, and SBTEnvelope sizes the helium into the sensors. The config was accepted and the overlap was silent. Guard it in readSBTConfig, alongside the existing hbeam_height check. Also correct three stale comments: - SBTEnvelope.h said helium_clearance_mm defaults to 0; it is 0.001 mm. - SBTConfig.cpp still described the pre-SBTEnvelope inset formula. - sbt.toml invited setting helium_clearance_mm = 0.0, which readSBTConfig rejects. Exact contact is reachable only by passing an SBTConfig directly. No geometry change; all 15 envelope sweep configurations unaffected. --- .../DecayVolume/include/DecayVolume/SBTConfig.h | 8 ++++++++ subsystems/DecayVolume/src/SBTConfig.cpp | 14 ++++++++++++-- 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/subsystems/DecayVolume/include/DecayVolume/SBTConfig.h b/subsystems/DecayVolume/include/DecayVolume/SBTConfig.h index b182623..a86b04d 100644 --- a/subsystems/DecayVolume/include/DecayVolume/SBTConfig.h +++ b/subsystems/DecayVolume/include/DecayVolume/SBTConfig.h @@ -127,6 +127,14 @@ struct SBTConfig { return topBottomContainerCentreY(y_half_mm) - topBottomContainerHalfThickness(); } /// Half-extent in X available to top/bottom containers (mm). + /// + /// The container stops half a flange width short of the frustum wall (that + /// is where the columns' inner face is), less a 1 mm gap so it does not + /// land flush against them. The 1 mm is inherited from the original + /// SBTSensorBuilder and is a literal, not sensor_clearance_mm — the two + /// happen to be equal at the current settings, which is worth being aware + /// of when changing sensor_clearance_mm, but they are not the same quantity + /// and this one is deliberately left as-is. double topBottomAvailX(double x_half_mm) const { return x_half_mm - 0.5 * hbeam_flange_width_mm - 1.0; } diff --git a/subsystems/DecayVolume/src/SBTConfig.cpp b/subsystems/DecayVolume/src/SBTConfig.cpp index 092ae47..e4fd860 100644 --- a/subsystems/DecayVolume/src/SBTConfig.cpp +++ b/subsystems/DecayVolume/src/SBTConfig.cpp @@ -165,10 +165,20 @@ SBTConfig readSBTConfig(const std::string& path) { requirePositive(cfg.container_thickness_mm, "container_thickness_mm"); requirePositive(cfg.wall_thickness_mm, "wall_thickness_mm"); // Clearances are insets/gaps; a non-positive value collapses the gap and - // can overlap geometry (helium_clearance_mm feeds the helium inset = - // container_thickness_mm + helium_clearance_mm in DecayVolumeFactory). + // can overlap geometry. helium_clearance_mm is the gap SBTEnvelope leaves + // between the helium and the innermost SBT surface. requirePositive(cfg.sensor_clearance_mm, "sensor_clearance_mm"); requirePositive(cfg.helium_clearance_mm, "helium_clearance_mm"); + // sensor_clearance_mm is trimmed off *each side* of a top/bottom container, + // so it must leave something behind: SBTConfig::topBottomContainerHalfThickness() + // is 0.5*container_thickness_mm - sensor_clearance_mm and is used directly as + // a GeoTrap half-dimension. Worse, if it goes negative, + // topBottomSensorInnerY() reports the container's inner face as *further out* + // than it is, and SBTEnvelope sizes the helium into the scintillator. + if (cfg.sensor_clearance_mm >= 0.5 * cfg.container_thickness_mm) + throw std::runtime_error( + "SBTConfig: 'sensor_clearance_mm' must be < 0.5 * 'container_thickness_mm' " + "(otherwise the top/bottom container half-thickness is non-positive)"); // Web half-height (= height/2 - flange_thickness) is used as a GeoBox // half-dimension, so the flanges must not consume the whole beam. if (cfg.hbeam_height_mm <= 2.0 * cfg.hbeam_flange_thickness_mm) From c9f7d66aacfc5ffb8592f977bc07936303ef826c Mon Sep 17 00:00:00 2001 From: Matei Climescu Date: Tue, 14 Jul 2026 21:35:17 +0200 Subject: [PATCH 3/5] fix(decayvolume): single contract for helium_clearance_mm; drop dead code MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Review follow-up (olantwin). helium_clearance_mm was documented three incompatible ways: "0 by default" (it is 0.001), "must be > 0" (readSBTConfig), and "must be >= 0" (SBTEnvelope). The real contract is >= 0 — zero is legal and gives exact face-to-face contact between the helium and the SBT, which the parameter sweep exercises. Drop the readSBTConfig check rather than the SBTEnvelope one: it is both redundant and stricter than the contract, and SBTEnvelope has to validate anyway for callers that construct an SBTConfig directly. State the bound once, in SBTConfig.h, and name SBTEnvelope as its only consumer. Also remove code the earlier refactor orphaned: - the zEnt_mm parameter on the frustum lambdas in SBTStructureBuilder, dead since they were routed through SBTConfig::xHalfAt/yHalfAt (15 call sites); - the GeoTrap, GeoNameTag, GeoTransform, GeoDefinitions and includes in DecayVolumeFactory, unused since buildHelium moved to SBTEnvelope. No geometry change; all 15 envelope sweep configurations unaffected. --- .../include/DecayVolume/SBTConfig.h | 8 +++-- subsystems/DecayVolume/sbt.toml | 9 ++--- .../DecayVolume/src/DecayVolumeFactory.cpp | 10 ++---- subsystems/DecayVolume/src/SBTConfig.cpp | 11 +++--- .../DecayVolume/src/SBTStructureBuilder.cpp | 34 +++++++++---------- 5 files changed, 37 insertions(+), 35 deletions(-) diff --git a/subsystems/DecayVolume/include/DecayVolume/SBTConfig.h b/subsystems/DecayVolume/include/DecayVolume/SBTConfig.h index a86b04d..7c2edb0 100644 --- a/subsystems/DecayVolume/include/DecayVolume/SBTConfig.h +++ b/subsystems/DecayVolume/include/DecayVolume/SBTConfig.h @@ -39,9 +39,11 @@ struct SBTConfig { double sensor_clearance_mm = 1.0; // ── Helium decay region ───────────────────────────────────────────── - // Gap left between the helium and the innermost SBT material, along the - // coordinate axis. Not a safety margin — just enough to avoid coincident - // surfaces, which Geant4's navigator handles badly. + // Gap left between the helium and the innermost SBT material, measured along + // the coordinate axis. Must be >= 0; SBTEnvelope enforces that and is the + // only consumer. It is not a safety margin — the default of 1 um is simply + // enough to stop the helium and the SBT sharing a surface, which Geant4's + // navigator handles badly. 0 is legal and gives exact face-to-face contact. // NOTE: this default must track sbt.toml. double helium_clearance_mm = 0.001; diff --git a/subsystems/DecayVolume/sbt.toml b/subsystems/DecayVolume/sbt.toml index 7aaf77e..3573a6e 100644 --- a/subsystems/DecayVolume/sbt.toml +++ b/subsystems/DecayVolume/sbt.toml @@ -61,10 +61,11 @@ sensor_clearance_mm = 1.0 # That is ~13.5 mm deeper than the scintillator's inner face. # # This is the gap left between the helium and the SBT surfaces, measured along -# the coordinate axis. It is NOT a safety margin: 1 um is far below any -# engineering tolerance and costs ~0.0001% of the fiducial volume. Its only job -# is to keep the helium and the SBT from sharing a surface, because coincident -# boundaries are a classic source of stuck tracks in Geant4's navigator. +# the coordinate axis. It must be >= 0, and is validated by SBTEnvelope, its only +# consumer. It is NOT a safety margin: 1 um is far below any engineering +# tolerance and costs ~0.0001% of the fiducial volume. Its only job is to keep +# the helium and the SBT from sharing a surface, because coincident boundaries +# are a classic source of stuck tracks in Geant4's navigator. # Set it to 0.0 only if you want exact face-to-face contact and have checked # that the navigator copes. helium_clearance_mm = 0.001 diff --git a/subsystems/DecayVolume/src/DecayVolumeFactory.cpp b/subsystems/DecayVolume/src/DecayVolumeFactory.cpp index 8566c32..44a0d44 100644 --- a/subsystems/DecayVolume/src/DecayVolumeFactory.cpp +++ b/subsystems/DecayVolume/src/DecayVolumeFactory.cpp @@ -11,13 +11,9 @@ #include "SHiPGeometry/SHiPMaterials.h" #include -#include #include #include -#include #include -#include -#include #include #include @@ -25,7 +21,6 @@ #include #include #include -#include // Source-tree / install-time fallbacks for sbt.toml, set by CMake. #ifndef SBT_TOML_DEFAULT_PATH @@ -108,8 +103,9 @@ GeoPhysVol* DecayVolumeFactory::build() { // The helium is not sized independently — it is *derived* from where the // SBT material actually is (SBTEnvelope, which reads the same placement // primitives on SBTConfig that the two builders above place with). It - // fills the interior exactly: no protrusion into steel or scintillator, - // and no arbitrary margin left behind (helium_clearance_mm = 0). + // fills the interior exactly: no protrusion into steel or scintillator, and + // no arbitrary margin left behind beyond helium_clearance_mm (1 um, just + // enough to avoid coincident surfaces). // // It is a stack of GeoTraps rather than one, because the inner surface is // not linear in Z: the side containers present a flat outer face over the diff --git a/subsystems/DecayVolume/src/SBTConfig.cpp b/subsystems/DecayVolume/src/SBTConfig.cpp index e4fd860..97de357 100644 --- a/subsystems/DecayVolume/src/SBTConfig.cpp +++ b/subsystems/DecayVolume/src/SBTConfig.cpp @@ -164,11 +164,14 @@ SBTConfig readSBTConfig(const std::string& path) { requirePositive(cfg.hbeam_web_thickness_mm, "hbeam_web_thickness_mm"); requirePositive(cfg.container_thickness_mm, "container_thickness_mm"); requirePositive(cfg.wall_thickness_mm, "wall_thickness_mm"); - // Clearances are insets/gaps; a non-positive value collapses the gap and - // can overlap geometry. helium_clearance_mm is the gap SBTEnvelope leaves - // between the helium and the innermost SBT surface. + // sensor_clearance_mm is trimmed off a container that must survive it, so it + // has to be strictly positive (and bounded above, below). requirePositive(cfg.sensor_clearance_mm, "sensor_clearance_mm"); - requirePositive(cfg.helium_clearance_mm, "helium_clearance_mm"); + // helium_clearance_mm is deliberately NOT checked here. It is consumed only + // by SBTEnvelope, which must validate it anyway for callers that construct + // an SBTConfig directly, and its contract is ">= 0" — 0 is legal and means + // the helium touches the SBT exactly. Checking it here as "> 0" as well + // would be both redundant and stricter than the real contract. // sensor_clearance_mm is trimmed off *each side* of a top/bottom container, // so it must leave something behind: SBTConfig::topBottomContainerHalfThickness() // is 0.5*container_thickness_mm - sensor_clearance_mm and is used directly as diff --git a/subsystems/DecayVolume/src/SBTStructureBuilder.cpp b/subsystems/DecayVolume/src/SBTStructureBuilder.cpp index 3407bcb..c2e1431 100644 --- a/subsystems/DecayVolume/src/SBTStructureBuilder.cpp +++ b/subsystems/DecayVolume/src/SBTStructureBuilder.cpp @@ -208,16 +208,16 @@ void SBTStructureBuilder::build(GeoVPhysVol* mother, const GeoMaterial* steel, c // Frustum profile from SBTConfig — the same accessors SBTEnvelope uses to // size the helium, so the two can never disagree. - auto xHalfAtZ = [&](double z_mm, double /*zEnt_mm*/) { return cfg.xHalfAt(z_mm); }; - auto yHalfAtZ = [&](double z_mm, double /*zEnt_mm*/) { return cfg.yHalfAt(z_mm); }; + auto xHalfAtZ = [&](double z_mm) { return cfg.xHalfAt(z_mm); }; + auto yHalfAtZ = [&](double z_mm) { return cfg.yHalfAt(z_mm); }; // (A) VERTICAL COLUMNS — 11 rows x 2 sides, frustum top -> floor. for (int row = 0; row <= nSubFrustrum; ++row) { const double z_mm = zEntrance_mm + row * subLength; const double z_G = z_mm * mm; - const double xEdge_mm = xHalfAtZ(z_mm, zEntrance_mm); - const double yTop_mm = yHalfAtZ(z_mm, zEntrance_mm); + const double xEdge_mm = xHalfAtZ(z_mm); + const double yTop_mm = yHalfAtZ(z_mm); const double yCol_ctr_mm = 0.5 * (yTop_mm + yFloor); const double yCol_half_mm = 0.5 * (yTop_mm - yFloor); @@ -253,10 +253,10 @@ void SBTStructureBuilder::build(GeoVPhysVol* mother, const GeoMaterial* steel, c const double zA_mm = zEntrance_mm + s * subLength + cbEndGap; const double zB_mm = zEntrance_mm + (s + 1) * subLength - cbEndGap; - const double xA = xHalfAtZ(zA_mm, zEntrance_mm); - const double xB = xHalfAtZ(zB_mm, zEntrance_mm); - const double yA = yHalfAtZ(zA_mm, zEntrance_mm); - const double yB = yHalfAtZ(zB_mm, zEntrance_mm); + const double xA = xHalfAtZ(zA_mm); + const double xB = xHalfAtZ(zB_mm); + const double yA = yHalfAtZ(zA_mm); + const double yB = yHalfAtZ(zB_mm); const std::string bname = tag + "_CornerBeam_" + std::to_string(ci) + "_S" + std::to_string(s); @@ -275,18 +275,18 @@ void SBTStructureBuilder::build(GeoVPhysVol* mother, const GeoMaterial* steel, c const double zLo_mm = zEntrance_mm + s * subLength; const double zHi_mm = zEntrance_mm + (s + 1) * subLength; - const double xLo = xHalfAtZ(zLo_mm, zEntrance_mm); - const double xHi = xHalfAtZ(zHi_mm, zEntrance_mm); + const double xLo = xHalfAtZ(zLo_mm); + const double xHi = xHalfAtZ(zHi_mm); // Straddle beam: the outer flange sits above the scintillator, the web // is omitted (it would pass through the cells) and the INNER FLANGE // HANGS BELOW THE SCINTILLATOR, into the decay region. That inner // flange is the innermost material in ±Y and therefore what bounds the // helium — see SBTConfig::longBeamInnerY() and SBTEnvelope. - const double yTop_Lo = +cfg.longBeamCentreY(yHalfAtZ(zLo_mm, zEntrance_mm)); - const double yTop_Hi = +cfg.longBeamCentreY(yHalfAtZ(zHi_mm, zEntrance_mm)); - const double yBot_Lo = -cfg.longBeamCentreY(yHalfAtZ(zLo_mm, zEntrance_mm)); - const double yBot_Hi = -cfg.longBeamCentreY(yHalfAtZ(zHi_mm, zEntrance_mm)); + const double yTop_Lo = +cfg.longBeamCentreY(yHalfAtZ(zLo_mm)); + const double yTop_Hi = +cfg.longBeamCentreY(yHalfAtZ(zHi_mm)); + const double yBot_Lo = -cfg.longBeamCentreY(yHalfAtZ(zLo_mm)); + const double yBot_Hi = -cfg.longBeamCentreY(yHalfAtZ(zHi_mm)); const std::string sTag = tag + "_SF" + std::to_string(s); @@ -379,9 +379,9 @@ void SBTStructureBuilder::build(GeoVPhysVol* mother, const GeoMaterial* steel, c const double z_mm = zEntrance_mm + row * subLength; const double z_G = z_mm * mm; - const double xEdge_mm = xHalfAtZ(z_mm, zEntrance_mm); - const double yTop_mm = +yHalfAtZ(z_mm, zEntrance_mm); - const double yBot_mm = -yHalfAtZ(z_mm, zEntrance_mm); + const double xEdge_mm = xHalfAtZ(z_mm); + const double yTop_mm = +yHalfAtZ(z_mm); + const double yBot_mm = -yHalfAtZ(z_mm); const std::string rowTag = tag + "_XBeam_R" + std::to_string(row); From 4d5aec5844bc08462aaa144bd693fedc03f3ebf1 Mon Sep 17 00:00:00 2001 From: Matei Climescu Date: Thu, 16 Jul 2026 15:55:25 +0200 Subject: [PATCH 4/5] docs(decayvolume): rigorous SAT bound, planarity guard, README refresh Review follow-up. minExpectedSeparation corrected only for single-axis tilt (max(gx,gy)). A surface tilted in both x and y tilts by sqrt(gx^2+gy^2), whose true Euclidean gap is below the old bound; at clr=10 that gap is mm-scale and would reject correct geometry if a doubly-tapered face were ever the binding pair. Use sqrt(1+gx^2+gy^2). The SAT axis set is exact only for planar-faced traps. Every trap built today is planar (dxdyn==dxdyp, alpha==0), but nothing enforced it; a future sheared trap would let separation() call an overlapping pair disjoint. Assert planarity in trapVertices so that fails loudly instead. Refresh README: the helium is 20 inscribed GeoTrap slabs (3400 traps, 3712 children), not a single central frustum, and is derived from the SBT surfaces. Document, but do not fix, a known conservative limitation: the Y envelope caps the helium at the longitudinal-beam inner flange over the full sub-frustum, but the beams stand off ~142 mm from each end, so the helium is up to ~13.5 mm short of the scintillator in those bands (~0.18 m^3, 0.04% of fiducial volume). Never overlaps; reclaiming it needs the flange z-projection handled and is deferred. No geometry change; SBTEnvelope.cpp diff is comment-only. --- subsystems/DecayVolume/README.md | 32 +++++++++++++++------ subsystems/DecayVolume/src/SBTEnvelope.cpp | 12 ++++++++ subsystems/DecayVolume/test_decayvolume.cpp | 27 +++++++++++++++-- 3 files changed, 61 insertions(+), 10 deletions(-) diff --git a/subsystems/DecayVolume/README.md b/subsystems/DecayVolume/README.md index a190381..3fd95cd 100644 --- a/subsystems/DecayVolume/README.md +++ b/subsystems/DecayVolume/README.md @@ -15,10 +15,16 @@ at the exit) made of two parts: - **LAB liquid-scintillator sensor cells** in thin aluminium containers, tiling the four faces of the frustum. -The decay region itself is a central **helium frustum**. It is sized strictly -inside the innermost sensor faces (`x_half − container_thickness`, -`y_half − container_thickness`) minus a clearance, so it cannot overlap the -structure or the sensors. +The decay region itself is **helium**, but it is not sized by a fixed inset. +It is *derived* from where the SBT material actually sits (see +[`SBTEnvelope.h`](include/DecayVolume/SBTEnvelope.h)): the innermost surfaces — +the side and top/bottom scintillator containers, and the top/bottom longitudinal +beam inner flanges, which hang *below* the sensor plane into the decay region — +bound it directly. Because that inner surface is not linear in Z (the side +containers present a flat outer face over the first part of each sub-frustum to +clear the columns), the helium is built as a stack of GeoTraps, two per +sub-frustum, each strictly inscribed inside the material minus a small +clearance. This cannot overlap the structure or the sensors. Everything is held in an **air container** (the SBT is part of the experiment's support frame, not a sealed vessel). The previous monolithic aluminium box @@ -30,7 +36,7 @@ children of the container), which the clash-avoidance logic relies on. ``` /SHiP/decay_volume (Air box, 4400 × 6600 × 50400 mm) - ├─ /SHiP/decay_volume/helium (PressurisedHe90, central frustum) + ├─ /SHiP/decay_volume/helium_0..19 (PressurisedHe90, 20 inscribed GeoTrap slabs) ├─ /SHiP/decay_volume/sbt/structure_* (Iron, 312 GeoBox H-beam pieces) └─ /SHiP/decay_volume/sbt/sensors_* (Al walls + LAB cells, 3380 GeoTrap) ``` @@ -47,6 +53,12 @@ on the container origin (entrance face at z = −25000 mm in the local frame). | Longitudinal beams | 60 | 1/face (sub-frustum 0–4), 2/face (5–9); flanges only | | Cross-beams | 66 | 11 rows × 2 faces × 3 boxes | +### Helium (20 GeoTrap) + +Two inscribed slabs per sub-frustum (10 sub-frusta): the slab boundaries fall at +each sub-frustum entrance and at the side-container flat-piece edge, so the slab +faces follow the sawtooth inner surface without cutting into it. + ### Sensors (3380 GeoTrap) 130 containers (80 side + 50 top/bottom), each Z-split into 2 pieces at the @@ -77,7 +89,7 @@ entrance Z, and all material names. Unknown keys are reported on stderr. ## Status -- [x] C++ implementation (SBT structure + sensors + helium frustum) +- [x] C++ implementation (SBT structure + sensors + derived helium) - [x] Frustum shape (replaces the old box vessel approximation) - [x] Surround Background Tagger integrated - [x] sbt.toml configuration @@ -86,5 +98,9 @@ entrance Z, and all material names. Unknown keys are reported on stderr. ## Tests `test_decayvolume` checks the container envelope, the 312 structure GeoBoxes, -the 3381 GeoTrap children (1 helium + 3380 sensors), the total of 3693 direct -children, and that the central decay region is a helium GeoTrap frustum. +and the 3400 GeoTrap children (20 helium slabs + 3380 sensors) for a total of +3712 direct children. Beyond the counts, it runs an exact separating-axis +overlap test between every helium slab and every SBT volume in the built tree, +asserting the helium neither protrudes into any material nor leaves a margin +beyond the configured clearance, and repeats that across 14 perturbed SBT +configurations so the guarantee survives re-parameterisation. diff --git a/subsystems/DecayVolume/src/SBTEnvelope.cpp b/subsystems/DecayVolume/src/SBTEnvelope.cpp index 352637a..406fdfc 100644 --- a/subsystems/DecayVolume/src/SBTEnvelope.cpp +++ b/subsystems/DecayVolume/src/SBTEnvelope.cpp @@ -57,6 +57,18 @@ double innerFreeHalfY(const SBTConfig& cfg, double z_mm) { // ... and the longitudinal beams' inner flange, which hangs below them. const double beam = cfg.longBeamInnerY(yHalf); // (Cross-beams sit a full beam-height above y_half and never reach in.) + // + // KNOWN LIMITATION (conservative): the longitudinal beams do not run the + // full sub-frustum length — they stand off by ~0.5*flange_width from each + // boundary (SBTStructureBuilder). In those ~142 mm end bands the true Y + // bound is the shallower scintillator, but we cap at the deeper beam line + // regardless, so the helium is up to ~13.5 mm short of the material there. + // This never overlaps (it only makes the helium smaller), but it does leave + // the helium slightly inside the "no unphysical margin" ideal over ~2.8 m of + // length — worth ~0.18 m^3, i.e. 0.04% of the fiducial volume. Reclaiming it + // means subdividing the Y envelope at the beam ends, which has to account + // for the flange box's z-projection overshooting its nominal end; deferred + // to a dedicated change rather than risking that interaction here. return std::min(sensor, beam); } diff --git a/subsystems/DecayVolume/test_decayvolume.cpp b/subsystems/DecayVolume/test_decayvolume.cpp index a5204f7..8e93739 100644 --- a/subsystems/DecayVolume/test_decayvolume.cpp +++ b/subsystems/DecayVolume/test_decayvolume.cpp @@ -181,7 +181,23 @@ std::array boxVertices(const GeoBox& b) { return v; } +// The SAT axis set used by separation() (6 face normals, each from 3 vertices, +// plus 12 edge crosses) is exact only when the trap's quadrilateral faces are +// planar. Every trap this geometry builds has equal top/bottom half-widths +// (dxdyn == dxdyp) and zero shear (alpha == 0), which makes the faces planar. +// A future config or builder change that produced a genuinely sheared trap +// would leave the axis set incomplete, and separation() could then call an +// overlapping pair disjoint. Require planarity here so that regresses loudly +// instead of silently weakening the overlap guarantee. +void requirePlanarTrap(const GeoTrap& t) { + REQUIRE(std::abs(t.getDxdyndzn() - t.getDxdypdzn()) < 1e-9); + REQUIRE(std::abs(t.getDxdyndzp() - t.getDxdypdzp()) < 1e-9); + REQUIRE(std::abs(t.getAngleydzn()) < 1e-12); + REQUIRE(std::abs(t.getAngleydzp()) < 1e-12); +} + std::array trapVertices(const GeoTrap& t) { + requirePlanarTrap(t); const double dz = t.getZHalfLength(); const double tt = std::tan(t.getTheta()); const double cx = tt * std::cos(t.getPhi()); @@ -323,9 +339,16 @@ constexpr double kTol = 1e-6; // helium_clearance_mm is a gap measured along a coordinate axis. SAT returns a // Euclidean distance, and the surfaces bounding the helium are tilted by the // frustum taper, so an axis gap of c shows up as c*cos(tilt). Assert the band. +// +// A bounding surface can tilt in both x and y at once (a side-container +// tracking-piece inner face does), and its normal then makes an angle +// atan(sqrt(gx^2 + gy^2)) with the axis, not atan(max(gx, gy)). The two-axis +// combination is the rigorous lower bound; max() alone overestimates the gap +// and can reject correct geometry once the clearance is large enough for the +// difference to exceed kTol (it does at the clr = 10 sweep case). double minExpectedSeparation(const SHiPGeometry::SBTConfig& cfg) { - const double g = std::max(std::abs(cfg.xGrowth()), std::abs(cfg.yGrowth())); - return cfg.helium_clearance_mm / std::sqrt(1.0 + g * g); + const double gx = cfg.xGrowth(), gy = cfg.yGrowth(); + return cfg.helium_clearance_mm / std::sqrt(1.0 + gx * gx + gy * gy); } // Closest approach between any helium slab and any SBT volume. From caf60247664ec0cdd26ef2b89fadb9f353630ef2 Mon Sep 17 00:00:00 2001 From: Matei Climescu Date: Thu, 16 Jul 2026 16:04:16 +0200 Subject: [PATCH 5/5] docs(decayvolume): clarify sbt/ is a name prefix, not a grouping volume The geometry-tree diagram indented structure_* and sensors_* under sbt, implying a physical grouping volume. There is none: the builders add every piece as a direct child of /SHiP/decay_volume, and /sbt/ is only a prefix on the volume names. Note that the tree is by name, not containment, consistent with the flat hierarchy test_decayvolume asserts. Names, types and counts unchanged. --- subsystems/DecayVolume/README.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/subsystems/DecayVolume/README.md b/subsystems/DecayVolume/README.md index 3fd95cd..916670f 100644 --- a/subsystems/DecayVolume/README.md +++ b/subsystems/DecayVolume/README.md @@ -41,6 +41,12 @@ children of the container), which the clash-avoidance logic relies on. └─ /SHiP/decay_volume/sbt/sensors_* (Al walls + LAB cells, 3380 GeoTrap) ``` +The tree above is by *name*, not by containment: `.../sbt/...` is a prefix on the +volume names, not a physical grouping volume. Every volume — helium slabs, +structure pieces and sensor pieces alike — is a direct child of +`/SHiP/decay_volume` (the clash-avoidance logic relies on this flat layout, and +`test_decayvolume` asserts all 3712 are direct children with no grandchildren). + Position in world: z = 58120 mm (centre), unchanged. The 50 m SBT is centred on the container origin (entrance face at z = −25000 mm in the local frame).