diff --git a/subsystems/DecayVolume/README.md b/subsystems/DecayVolume/README.md index a190381..916670f 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,11 +36,17 @@ 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) ``` +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). @@ -47,6 +59,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 +95,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 +104,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/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..7c2edb0 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,13 @@ 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, 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; // ── Material names (resolved from SHiPMaterials) ──────────────────── std::string material_steel = "Iron"; @@ -54,6 +61,102 @@ 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). + /// + /// 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; + } + + // --- 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..3573a6e 100644 --- a/subsystems/DecayVolume/sbt.toml +++ b/subsystems/DecayVolume/sbt.toml @@ -44,10 +44,31 @@ 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 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 # ── 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..44a0d44 100644 --- a/subsystems/DecayVolume/src/DecayVolumeFactory.cpp +++ b/subsystems/DecayVolume/src/DecayVolumeFactory.cpp @@ -4,19 +4,16 @@ #include "DecayVolume/DecayVolumeFactory.h" #include "DecayVolume/SBTConfig.h" +#include "DecayVolume/SBTEnvelope.h" #include "DecayVolume/SBTSensorBuilder.h" #include "DecayVolume/SBTStructureBuilder.h" #include "SHiPGeometry/ConfigPath.h" #include "SHiPGeometry/SHiPMaterials.h" #include -#include #include #include -#include #include -#include -#include #include #include @@ -50,7 +47,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 +95,25 @@ 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 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 + // 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/SBTConfig.cpp b/subsystems/DecayVolume/src/SBTConfig.cpp index 092ae47..97de357 100644 --- a/subsystems/DecayVolume/src/SBTConfig.cpp +++ b/subsystems/DecayVolume/src/SBTConfig.cpp @@ -164,11 +164,24 @@ 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 feeds the helium inset = - // container_thickness_mm + helium_clearance_mm in DecayVolumeFactory). + // 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 + // 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) diff --git a/subsystems/DecayVolume/src/SBTEnvelope.cpp b/subsystems/DecayVolume/src/SBTEnvelope.cpp new file mode 100644 index 0000000..406fdfc --- /dev/null +++ b/subsystems/DecayVolume/src/SBTEnvelope.cpp @@ -0,0 +1,187 @@ +// 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.) + // + // 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); +} + +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..c2e1431 100644 --- a/subsystems/DecayVolume/src/SBTStructureBuilder.cpp +++ b/subsystems/DecayVolume/src/SBTStructureBuilder.cpp @@ -198,37 +198,26 @@ 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) { 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); @@ -264,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); @@ -286,15 +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 yBeamOffset = 0.5 * hBeamHw; // centred C-channel + const double xLo = xHalfAtZ(zLo_mm); + const double xHi = xHalfAtZ(zHi_mm); - 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)); + 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); @@ -387,13 +379,13 @@ 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); - 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..8e93739 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,396 @@ 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; +} + +// 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()); + 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. +// +// 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 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. +// > 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)); + } +}