diff --git a/src/SHiPGeometry.cpp b/src/SHiPGeometry.cpp index 7c81059..4c594a0 100644 --- a/src/SHiPGeometry.cpp +++ b/src/SHiPGeometry.cpp @@ -48,22 +48,24 @@ GeoPhysVol* SHiPGeometryBuilder::build() { // Note: These are relative to the cave origin placeChild(world, target, "/SHiP/target", 1, GeoTrf::Translate3D(0.0, -14.45 * cm, 43.25 * cm)); - // Build and place MuonShieldArea - // GDML z range: 204–3148.66 cm → centre: 1676.33 cm = 16763.3 mm from world origin + // Build the muon shield, with the neutrino detector embedded inside it. + // + // The SND is an independent subsystem, but in volume terms it is a daughter + // of the muon-shield container: it sits within the shield region (WARM SND + // slot 26.40–31.50 m → centre 28.95 m). We build it first, then nest it in + // the shield container. The muon-shield block list must leave this slot + // free of iron. MuonShieldFactory muonShieldFactory(materials); - GeoPhysVol* muonShield = muonShieldFactory.build(); - placeChild(world, muonShield, "/SHiP/muon_shield", 2, - GeoTrf::Translate3D(0.0, 0.0, 16763.3 * mm)); - // Build and place the Scattering and Neutrino Detector (SND). - // Z: 26.40 to 31.50 m (WARM muon-shield configuration) → centre 28.95 m. - // The SND sits within the downstream end of the muon-shield region, so its - // envelope overlaps the muon-shield container by design (see test_consistency). NeutrinoDetectorFactory neutrinoDetectorFactory(materials); GeoPhysVol* neutrinoDetector = neutrinoDetectorFactory.build(); - placeChild(world, neutrinoDetector, "/SHiP/neutrino_detector", 9, - GeoTrf::Translate3D(0.0, 0.0, 28.95 * m)); + muonShieldFactory.embedDaughter(neutrinoDetector, 28.95 * 1000.0, "/SHiP/neutrino_detector"); + // The container is built centred on its own origin, so it is placed at the + // envelope centre reported by the factory after build(). + GeoPhysVol* muonShield = muonShieldFactory.build(); + placeChild(world, muonShield, "/SHiP/muon_shield", 2, + GeoTrf::Translate3D(0.0, 0.0, muonShieldFactory.centreZ_mm())); // Build and place UpstreamTagger (sensitive scintillator slab) // Z: 32.52 to 32.92 m → centre: 32.72 m SHiPUBTManager ubtManager; diff --git a/subsystems/MuonShield/CMakeLists.txt b/subsystems/MuonShield/CMakeLists.txt index 16222ae..8507094 100644 --- a/subsystems/MuonShield/CMakeLists.txt +++ b/subsystems/MuonShield/CMakeLists.txt @@ -1,4 +1,43 @@ # SPDX-License-Identifier: LGPL-3.0-or-later # Copyright (C) CERN for the benefit of the SHiP Collaboration +# Shared recipe: library (globs src/*.cpp), include dirs, kernel link, Catch test. ship_add_subsystem(MuonShield) + +# --- Extras: MuonShield parses MS.toml via MuonShieldConfig (toml++), same as +# the calorimeter. Everything below mirrors subsystems/Calorimeter/CMakeLists.txt. + +find_package(tomlplusplus REQUIRED) + +# toml++ is a build-time-only dependency (used only inside MuonShieldConfig.cpp, +# never in a public header). BUILD_INTERFACE keeps it out of the install export. +target_link_libraries( + MuonShield + PRIVATE $ +) + +# Stage MS.toml next to the tests and into the top-level build dir (for +# test_builder / OverlapCheck), so it is found when running from either. +configure_file( + ${CMAKE_CURRENT_SOURCE_DIR}/MS.toml + ${CMAKE_CURRENT_BINARY_DIR}/MS.toml + COPYONLY +) +configure_file( + ${CMAKE_CURRENT_SOURCE_DIR}/MS.toml + ${CMAKE_BINARY_DIR}/MS.toml + COPYONLY +) + +# Absolute fallbacks so the factory finds MS.toml when the CWD does not have it. +target_compile_definitions( + MuonShield + PRIVATE + MS_TOML_DEFAULT_PATH="${CMAKE_CURRENT_SOURCE_DIR}/MS.toml" + MS_TOML_INSTALL_PATH="${CMAKE_INSTALL_PREFIX}/${CMAKE_INSTALL_DATADIR}/SHiPGeometry/MS.toml" +) + +install( + FILES ${CMAKE_CURRENT_SOURCE_DIR}/MS.toml + DESTINATION ${CMAKE_INSTALL_DATADIR}/SHiPGeometry +) diff --git a/subsystems/MuonShield/MS.toml b/subsystems/MuonShield/MS.toml new file mode 100644 index 0000000..d853921 --- /dev/null +++ b/subsystems/MuonShield/MS.toml @@ -0,0 +1,58 @@ +# SPDX-License-Identifier: LGPL-3.0-or-later +# Copyright (C) CERN for the benefit of the SHiP Collaboration +# +# Muon shield geometry — solid-block approximation of the FairShip TRY_2026 +# magnet layout (FairShip PR #1334). Each magnet is one solid iron block from +# its outer envelope (dX*(1+ratio)+midGap+ceil(max(100/dY,gap)) x dY+dY_yoke); +# the aperture, yoke shape and field are dropped (field handled elsewhere). +# +# The muon shield is defined INDEPENDENTLY of the neutrino detector. The SND +# reservation is declared in NeutrinoDetector/SD.toml and subtracted from the +# iron at build time (A - B), so this file needs no knowledge of the SND. +# +# The FairShip absorber (row 0) is the magnetised hadron stopper and lives in +# the Target subsystem (2.14-4.44 m); the shield front is at 2.14 m, so magnet 1 +# begins at 4.59 m. + +block_material = "Iron" + +envelope_half_x_mm = 1760.0 +envelope_half_y_mm = 1320.0 +envelope_z_start_m = 4.5400 +envelope_z_end_m = 32.0800 + +# magnet 1 +[[block]] +start = [0.0, 0.0, 4590.0] +size = [2720.0, 1600.0, 3000.0] + +# magnet 2 +[[block]] +start = [0.0, 0.0, 7740.0] +size = [2980.0, 1780.0, 4500.0] + +# magnet 3 +[[block]] +start = [0.0, 0.0, 12400.0] +size = [3104.0, 1700.0, 4500.0] + +# magnet 4 +[[block]] +start = [0.0, 0.0, 17050.0] +size = [3184.0, 1200.0, 4500.0] + +# magnet 5 +[[block]] +start = [0.0, 0.0, 21710.0] +size = [2680.0, 1040.0, 2620.0] + +# magnet 6 +[[block]] +start = [0.0, 0.0, 24480.0] +size = [3000.0, 1260.0, 3360.0] + + +# magnet 7 +[[block]] +start = [0.0, 0.0, 28030.0] +size = [3420.0, 2540.0, 4000.0] diff --git a/subsystems/MuonShield/README.md b/subsystems/MuonShield/README.md index 2f5fbbf..2c4b096 100644 --- a/subsystems/MuonShield/README.md +++ b/subsystems/MuonShield/README.md @@ -1,47 +1,125 @@ # MuonShield -Active muon shield magnets for background suppression. +Passive iron muon shield for background suppression. ## Description -The MuonShield subsystem implements 6 magnet stations that deflect muons away from the detector acceptance. The current implementation uses bounding-box approximations of the GDML arb8 (trapezoidal) shapes — each station contains 8 iron pieces (left/right middle magnets, left/right return yokes, and 4 top/bottom corner pieces). +The MuonShield subsystem reserves an Air **envelope** and places an explicit +list of iron **blocks** inside it. Each block is a box, or a symmetric +trapezoid/frustum (`GeoTrd`) when tapered, anchored by its upstream face and +optionally rotated. The block list is the sole source of iron: holes and +apertures are represented by leaving a region free of blocks, or by arranging +several blocks around the gap. + +The default `MS.toml` is the **FairShip TRY_2026** magnet layout (FairShip +PR #1334) as a solid-block approximation: the 7 shield magnets, each reproduced +as one iron block from its outer envelope (aperture, yoke and field dropped). +It is generated by `fairship_to_mstoml.py`, which converts the FairShip +`ShipMuonShield` `params` table into blocks. The FairShip absorber (row 0) is +the magnetised hadron stopper and lives in the Target subsystem, so it is not +included here. ## Geometry Tree +The default `MS.toml` places the 7 TRY_2026 magnets as solid blocks. The SND +reservation (`NeutrinoDetector/SD.toml`) is subtracted from the iron at build +time, and the SND is embedded as a daughter that sits in the resulting cavity: + +```text +muon_shield (Air envelope, 3520 × 2640 × 27540 mm) centre z = 18.31 m + ├─ block_0 magnet 1 (Iron box) + ├─ … magnets 2–5 (Iron boxes; magnet 6 is a GeoTrd, widens in Y) + ├─ block_5 magnet 6 (Iron − SND box) ┐ carved by A − B + ├─ block_6 magnet 7 (Iron − SND box) ┘ (magnets the reservation intersects) + └─ neutrino_detector (Air box, 26.40 → 31.50 m) embedded, sits in the cavity ``` -MuonShieldArea (Air, 3620×3400×29448 mm) - ├─ MagnAbsorb_container (Air) z = -13568.3 mm - │ └─ 8 × Iron boxes (MiddleMagL/R, MagRetL/R, MagTop/BotLeft/Right) - ├─ Magn1_container (Air) z = -7263.3 mm - │ └─ 8 × Iron boxes - ├─ Magn2_container (Air) z = +591.5 mm - │ └─ 8 × Iron boxes - ├─ Magn3_container (Air) z = +5821.6 mm - │ └─ 8 × Iron boxes - ├─ Magn4_container (Air) z = +9096.9 mm - │ └─ 8 × Iron boxes - └─ Magn5_container (Air) z = +12385.1 mm - └─ 8 × Iron boxes + +The envelope values in `MS.toml` are generated to contain the magnets +(3520 × 2640 mm, 4.54–32.08 m); the container is placed in the world at its Z +centre (18.31 m) by `SHiPGeometryBuilder`. The SND cavity is a real Boolean +subtraction (like FairShip's `SetSNDSpace`): the reservation box is subtracted +from whichever magnets it intersects — so which magnets get carved follows from +the box position, not a hard-coded list. See the **Reserved space** section +below. + +## Configuration (`MS.toml`) + +| Key | Fallback | Meaning | +|----------------------|----------|------------------------------------------------| +| `block_material` | "Iron" | Absorber material (must exist in SHiPMaterials) | +| `envelope_half_x_mm` | 2100 | Container half-X (mm) | +| `envelope_half_y_mm` | 2300 | Container half-Y (mm) | +| `envelope_z_start_m` | 4.74 | Envelope start along Z (m, world coords) | +| `envelope_z_end_m` | 32.22 | Envelope end along Z (m, world coords) | +| `[[block]]` | — | One table per iron block (see below) | + +The **Fallback** column is the parser default used only when a key is *omitted*. +The shipped `MS.toml` sets the envelope explicitly (half-sizes 1760 × 1320 mm, +z = 4.54–32.08 m — the values in the geometry tree above), generated by +`fairship_to_mstoml.py` to contain the blocks. The **envelope** is the Air +container and the subsystem's world placement. The parser rejects an +inverted/zero envelope and any block with a non-positive size, a taper that +collapses the downstream face, or a face outside the envelope in Z. + +### Blocks + +Each `[[block]]` is a table with these fields (all lengths in mm, angles in +degrees, world/beamline coordinates): + +```toml +[[block]] +start = [0.0, 0.0, 4740.0] # centre of the UPSTREAM (-z) face +size = [3000.0, 2000.0, 21660.0] # near-face full x, y, and z length +rotation = [0.0, 2.0, 0.0] # optional: extrinsic X→Y→Z about `start` +taper = [5.0, 0.0] # optional: symmetric half-opening angles ``` -Position in world: z = 16763.3 mm (centre of MuonShieldArea). +- `start` is the centre of the block's upstream (−z) face; the block extends + downstream by `size[2]` before rotation. +- `rotation` (default `[0,0,0]`) is applied about `start`, extrinsically about + the world X then Y then Z axes. +- `taper` (default `[0,0]`) gives symmetric half-opening angles: the section + widens (or, if negative, narrows) toward the downstream face, so the far + half-x is `size[0]/2 + size[2]·tan(taper[0])` and similarly for Y. `[0,0]` + yields a `GeoBox`; otherwise a `GeoTrd`. + +### Reserved space (SND cavity) + +The shield iron is carved by Boolean subtraction (A − B) so external detectors +can define their space independently. `MuonShieldFactory::reserveSpace(centre, +size, rotation)` registers a box (world coords, mm); at build time every magnet +the box intersects has it subtracted (`GeoShapeSubtraction`), while other +magnets are untouched. The **neutrino detector** reservation is declared in +`NeutrinoDetector/SD.toml` and wired in by `SHiPGeometryBuilder`, which reads the +envelope, calls `reserveSpace(...)`, then `embedDaughter(...)` to place the SND +in the cavity. Because the carve follows from the box position, the muon shield +never names which magnets hold the SND — moving the box in `SD.toml` moves the +cavity. + +### Embedded daughters + +Other subsystems can be nested inside the shield container while remaining +independent subsystems. `MuonShieldFactory::embedDaughter(vol, worldCentreZ_mm, +name)` registers a pre-built volume, placed at build time at the given world-Z +inside the container. The SND keeps its own factory, config, and +`/SHiP/neutrino_detector` naming — only its position in the volume tree changes +(it becomes a daughter of `/SHiP/muon_shield` rather than a direct child of the +world), sitting in the cavity carved by its `reserveSpace` reservation. ## Materials -| Material | Density | Usage | -|----------|-----------|------------------| -| Air | 1.29 mg/cm³ | Container volumes | -| Iron | 7.87 g/cm³ | Magnet pieces | +| Material | Density | Usage | +|----------|-------------|-------------------| +| Air | 1.29 mg/cm³ | Container volume | +| Iron | 7.87 g/cm³ | Shield blocks | ## Status -- [x] C++ implementation (box approximations) -- [ ] Replace boxes with proper GeoTrap/arb8 shapes -- [ ] Verify field map integration points -- [ ] Verification against GDML +- [x] Explicit block list driven by `MS.toml` (position, rotation, taper) +- [ ] Magnetic field (handled elsewhere; blocks are passive iron) +- [ ] Populate the block list from the engineering reference layout ## TODO -- Replace bounding-box iron pieces with proper trapezoidal (GeoTrap or GeoGenericTrap) shapes matching the GDML arb8 vertices -- Add magnetic field regions (currently geometry only, no field) -- Verify station z-positions and piece dimensions against GDML reference +- Replace the placeholder default blocks with the real muon-shield magnet + layout supplied by the subsystem coordinator. diff --git a/subsystems/MuonShield/include/MuonShield/MuonShieldConfig.h b/subsystems/MuonShield/include/MuonShield/MuonShieldConfig.h new file mode 100644 index 0000000..62921b4 --- /dev/null +++ b/subsystems/MuonShield/include/MuonShield/MuonShieldConfig.h @@ -0,0 +1,90 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) CERN for the benefit of the SHiP Collaboration + +#pragma once + +#include +#include +#include + +namespace SHiPGeometry { + +/// Degrees-to-radians conversion, shared by the config parser and the factory. +inline constexpr double kDegToRad = 3.14159265358979323846 / 180.0; + +/** + * @brief A single muon-shield iron block. + * + * A block is a box (or, when tapered, a symmetric trapezoid/frustum) placed by + * its upstream (−z) face. All lengths are in millimetres and all angles in + * degrees, in world/beamline coordinates (target-front-face origin). + * + * Geometry conventions: + * - @c start is the centre of the block's upstream face; the block extends + * downstream (+z) by @c size[2] before rotation. + * - @c rotation is applied about @c start, extrinsically about the world X, Y + * then Z axes. + * - @c size = {sx, sy, sz}: the full X and Y of the *upstream* face and the Z + * length. + * - @c taper = {ax, ay}: symmetric half-opening angles. The cross-section + * grows (or, if negative, shrinks) toward the downstream face, so the far + * half-width is (sx/2 + sz·tan(ax)) and similarly for Y. {0, 0} → a box. + */ +struct MuonShieldBlock { + std::array start = {0.0, 0.0, 0.0}; ///< upstream-face centre (mm) + std::array size = {0.0, 0.0, 0.0}; ///< near-face sx, sy and length sz (mm) + std::array rotation_deg = {0.0, 0.0, 0.0}; ///< rotation about x, y, z (deg) + std::array taper_deg = {0.0, 0.0}; ///< half-opening angles ax, ay (deg) +}; + +/** + * @brief Configuration for the SHiP muon shield geometry. + * + * Populated by readMuonShieldConfig() from an MS.toml file. + * + * The shield is an explicit list of iron @c blocks placed inside an Air + * @c envelope container. The envelope sizes the container and fixes the + * subsystem's world placement (its Z centre); the blocks are the only iron. + * Holes/apertures are represented simply by the absence of a block, or by + * arranging several blocks around the gap. + */ +struct MuonShieldConfig { + // Absorber material for all blocks (must exist in SHiPMaterials). + std::string block_material = "Iron"; + + // Container (envelope) transverse half-sizes (mm). + double envelope_half_x_mm = 2100.0; + double envelope_half_y_mm = 2300.0; + + // Envelope Z extent in world coordinates (target-front-face origin, m). + // This is the Air container and the world placement of the subsystem. + double envelope_z_start_m = 4.74; + double envelope_z_end_m = 32.22; + + // The iron blocks. Empty = no iron. + std::vector blocks; + + // ── Derived helpers ───────────────────────────────────────────────── + /// Full Z length of the envelope (mm). + double envelopeLengthZ_mm() const { return (envelope_z_end_m - envelope_z_start_m) * 1000.0; } + /// World-Z centre of the envelope (mm). + double envelopeCentreZ_mm() const { + return 0.5 * (envelope_z_start_m + envelope_z_end_m) * 1000.0; + } +}; + +/** + * @brief Parse an MS.toml file and return a MuonShieldConfig. + * + * Uses toml++ for parsing. Unknown top-level keys are reported via stderr + * (helpful when a stale or mistyped key would otherwise be silently ignored), + * but do not cause the parse to fail. + * + * @throws std::runtime_error if the file cannot be opened, contains malformed + * TOML, has a non-positive or inverted Z envelope, or a block with a + * non-positive size, a taper that collapses the far face, or an + * upstream face outside the envelope in Z. + */ +MuonShieldConfig readMuonShieldConfig(const std::string& path); + +} // namespace SHiPGeometry diff --git a/subsystems/MuonShield/include/MuonShield/MuonShieldFactory.h b/subsystems/MuonShield/include/MuonShield/MuonShieldFactory.h index 366aad8..7e0f650 100644 --- a/subsystems/MuonShield/include/MuonShield/MuonShieldFactory.h +++ b/subsystems/MuonShield/include/MuonShield/MuonShieldFactory.h @@ -3,60 +3,110 @@ #pragma once -#include +#include +#include +#include class GeoPhysVol; namespace SHiPGeometry { class SHiPMaterials; +struct MuonShieldConfig; /** - * @brief Factory for the MuonShield (muon shield magnets) geometry + * @brief Factory for the MuonShield geometry. * - * Creates a MuonShieldArea (Air) container holding 6 magnet stations. - * Each station contains 8 Iron pieces — bounding-box approximations of the - * GDML arb8 shapes. Station z-positions (from GDML MuonShieldArea origin, cm): - * MagnAbsorb: 319.5, Magn1: 950, Magn2: 1735.48, - * Magn3: 2258.49, Magn4: 2586.02, Magn5: 2914.84 + * Builds an Air container spanning the muon-shield Z envelope and places an + * explicit list of iron blocks inside it (from MS.toml). Each block is a box, + * or a symmetric trapezoid/frustum (GeoTrd) when tapered, anchored by its + * upstream face and optionally rotated. + * + * Other subsystems can be embedded as daughter volumes of the shield container + * via embedDaughter() — used, for example, to nest the neutrino detector inside + * the muon shield while keeping it an independent subsystem. + * + * The container is built centred on its own origin; SHiPGeometryBuilder places + * it in the world at centreZ_mm() (the envelope centre from MS.toml). + * + * The config file is resolved at build() time: + * 1. "MS.toml" relative to the current working directory (works when running + * from the build directory, where CMake stages the file). + * 2. The absolute source-tree path baked in via MS_TOML_DEFAULT_PATH. + * 3. The installed data-dir path via MS_TOML_INSTALL_PATH. */ class MuonShieldFactory { public: - explicit MuonShieldFactory(SHiPMaterials& materials); + explicit MuonShieldFactory(SHiPMaterials& materials, std::string configPath = "MS.toml"); ~MuonShieldFactory() = default; - [[nodiscard]] GeoPhysVol* build(); + /** + * @brief Embed an external volume as a daughter of the shield container. + * + * Registers @p daughter to be placed inside the shield at world-Z + * @p worldCentreZ_mm. The daughter remains an independent subsystem; this + * only nests it in the volume tree. Since the iron is defined by an explicit + * block list, the caller is responsible for not listing blocks over the + * daughter's footprint. Call before build(). + * + * @param daughter Pre-built volume (centred on its own origin). + * @param worldCentreZ_mm - private: - SHiPMaterials& m_materials; - // Unit shorthand (GeoModel's native length unit is mm) - static constexpr double mm = GeoModelKernelUnits::mm; + centre where the daughter is placed (mm). + * @param name Name tag for the placement. + */ + void embedDaughter(GeoPhysVol* daughter, double worldCentreZ_mm, const std::string& name); - struct PieceData { - double halfX, halfY, halfZ; // bounding-box half-sizes - double centX, centY; // centre offset in station XY frame - const char* name; // piece name suffix - }; + /** + * @brief Reserve (carve) a box out of the shield iron via Boolean A - B. + * + * Registers a box to subtract from every magnet it intersects, so external + * detectors (e.g. the neutrino detector, from SD.toml) can define their + * space independently and have it removed from the iron at build() time. + * Positions are in world/beamline coordinates. Call before build(). + * + * @param worldCentre_mm Box centre in world coordinates (mm). + * @param size_mm Full box dimensions {x, y, z} (mm). + * @param rotation_deg Optional extrinsic X->Y->Z rotation about the centre. + */ + void reserveSpace(const std::array& worldCentre_mm, + const std::array& size_mm, + const std::array& rotation_deg = {0.0, 0.0, 0.0}); - struct StationData { - const char* name; // station name - double stationZ; // z in MuonShieldArea frame - double containerHalfX; - double containerHalfY; - double containerHalfZ; - PieceData pieces[8]; + /** Build and return the muon-shield container volume (centred on origin). */ + GeoPhysVol* build(); + + /** + * @brief World-Z centre (mm) at which the container should be placed. + * + * Valid only after build() has been called (it is read from MS.toml). + */ + double centreZ_mm() const { return m_centreZ_mm; } + + /** Return the config path that will actually be opened (after resolution). */ + std::string resolvedConfigPath() const; + + private: + struct EmbeddedDaughter { + GeoPhysVol* volume; + double worldCentreZ_mm; + std::string name; }; - // GDML-derived station data for all 6 stations × 8 pieces - static const StationData k_stations[6]; + struct ReservedBox { + std::array centre_mm; + std::array size_mm; + std::array rotation_deg; + }; - GeoPhysVol* buildStation(const StationData& station); + SHiPMaterials& m_materials; + std::string m_configPath; + std::vector m_daughters; + std::vector m_reservations; - // MuonShieldArea container dimensions - static constexpr double s_areaHalfX = 1810.0 * mm; - static constexpr double s_areaHalfY = 1700.0 * mm; - static constexpr double s_areaHalfZ = 14724.0 * mm; + // World-Z centre of the container, populated by build() from the config. + double m_centreZ_mm = 0.0; }; } // namespace SHiPGeometry diff --git a/subsystems/MuonShield/src/MuonShieldConfig.cpp b/subsystems/MuonShield/src/MuonShieldConfig.cpp new file mode 100644 index 0000000..67ee191 --- /dev/null +++ b/subsystems/MuonShield/src/MuonShieldConfig.cpp @@ -0,0 +1,201 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) CERN for the benefit of the SHiP Collaboration +// +// Muon shield configuration parser. Reads an MS.toml file and populates a +// MuonShieldConfig struct. Built on toml++ +// (https://github.com/marzer/tomlplusplus, single-header MIT, provided by the +// tomlplusplus package), so the file format is standard TOML. +// +// Behaviour notes +// =============== +// * Unknown top-level keys are NOT silently ignored. Each one is reported on +// stderr at parse time. This catches typos and stale keys. +// +// * Numeric fields accept both TOML integers and floats. + +#include "MuonShield/MuonShieldConfig.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace SHiPGeometry { + +namespace { + +using namespace std::string_view_literals; + +// Recognised top-level keys (sorted for binary search). Anything outside this +// set triggers a warning. "block" is the [[block]] array of tables. +static constexpr std::array kKnownKeys = { + "block"sv, + "block_material"sv, + "envelope_half_x_mm"sv, + "envelope_half_y_mm"sv, + "envelope_z_end_m"sv, + "envelope_z_start_m"sv, +}; + +// Mapping from TOML key name to MuonShieldConfig double member pointer. +struct NumericField { + const char* key; + double MuonShieldConfig::* member; +}; + +static constexpr NumericField kNumericFields[] = { + {"envelope_half_x_mm", &MuonShieldConfig::envelope_half_x_mm}, + {"envelope_half_y_mm", &MuonShieldConfig::envelope_half_y_mm}, + {"envelope_z_start_m", &MuonShieldConfig::envelope_z_start_m}, + {"envelope_z_end_m", &MuonShieldConfig::envelope_z_end_m}, +}; + +// Extract a TOML numeric (integer or float) as a double, if the node holds one. +std::optional asDouble(const toml::node* node) { + if (node) { + if (auto d = node->value()) + return *d; + if (auto i = node->value()) + return static_cast(*i); + } + return std::nullopt; +} + +// Read a double or integer as a double. +double readNumeric(const toml::node_view& node, const std::string& key) { + if (auto v = asDouble(node.node())) + return *v; + throw std::runtime_error("MuonShieldConfig: '" + key + "' must be a number"); +} + +// Read a single toml node as a double (int or float). +double nodeToNumber(const toml::node* node, const std::string& context, const std::string& path) { + if (auto v = asDouble(node)) + return *v; + throw std::runtime_error("MuonShieldConfig: " + context + " values must be numbers in " + path); +} + +// Read a fixed-length numeric array (e.g. size = [x, y, z]) from a block table. +template +std::array readFixedArray(const toml::table& block, const char* key, + const std::string& path, bool required, + const std::array& fallback) { + auto node = block[key]; + if (!node) { + if (required) + throw std::runtime_error("MuonShieldConfig: block is missing required '" + + std::string(key) + "' in " + path); + return fallback; + } + const toml::array* arr = node.as_array(); + if (!arr || arr->size() != N) + throw std::runtime_error("MuonShieldConfig: block '" + std::string(key) + "' must be an " + + std::to_string(N) + "-element array in " + path); + std::array out{}; + for (std::size_t i = 0; i < N; ++i) + out[i] = nodeToNumber(arr->get(i), "block '" + std::string(key) + "'", path); + return out; +} + +} // namespace + +MuonShieldConfig readMuonShieldConfig(const std::string& path) { + MuonShieldConfig cfg; + + toml::table table; + try { + table = toml::parse_file(path); + } catch (const toml::parse_error& e) { + throw std::runtime_error("MuonShieldConfig: failed to parse " + path + ": " + + std::string(e.description())); + } + + // First pass: warn about unknown keys. + for (const auto& [k, _] : table) { + if (!std::ranges::binary_search(kKnownKeys, std::string_view{k})) { + std::cerr << "MuonShieldConfig: warning: unknown key '" << k << "' in " << path + << " (typo? stale field? — value will be ignored)\n"; + } + } + + // Numeric (double) envelope fields. + for (const auto& [key, member] : kNumericFields) + if (auto n = table[key]; n) + cfg.*member = readNumeric(n, key); + + // String field. + if (auto n = table["block_material"]; n) { + if (auto s = n.value(); s) + cfg.block_material = *s; + else + throw std::runtime_error("MuonShieldConfig: 'block_material' must be a string"); + } + + // Blocks: an array of tables ([[block]]). + if (auto blocksNode = table["block"]; blocksNode) { + const toml::array* arr = blocksNode.as_array(); + if (!arr) + throw std::runtime_error("MuonShieldConfig: 'block' must be an array of tables in " + + path); + for (const auto& elem : *arr) { + const toml::table* bt = elem.as_table(); + if (!bt) + throw std::runtime_error( + "MuonShieldConfig: each [[block]] entry must be a table in " + path); + MuonShieldBlock block; + block.start = readFixedArray<3>(*bt, "start", path, true, {0.0, 0.0, 0.0}); + block.size = readFixedArray<3>(*bt, "size", path, true, {0.0, 0.0, 0.0}); + block.rotation_deg = readFixedArray<3>(*bt, "rotation", path, false, {0.0, 0.0, 0.0}); + block.taper_deg = readFixedArray<2>(*bt, "taper", path, false, {0.0, 0.0}); + cfg.blocks.push_back(block); + } + } + + // ── Validation ────────────────────────────────────────────────────── + if (cfg.envelope_z_end_m <= cfg.envelope_z_start_m) + throw std::runtime_error( + "MuonShieldConfig: envelope_z_end_m must be greater than envelope_z_start_m in " + + path); + if (cfg.envelope_half_x_mm <= 0.0 || cfg.envelope_half_y_mm <= 0.0) + throw std::runtime_error( + "MuonShieldConfig: envelope_half_x_mm and envelope_half_y_mm must be positive in " + + path); + + const double envStartMm = cfg.envelope_z_start_m * 1000.0; + const double envEndMm = cfg.envelope_z_end_m * 1000.0; + constexpr double kEps = 1e-6; // mm + + for (std::size_t i = 0; i < cfg.blocks.size(); ++i) { + const auto& b = cfg.blocks[i]; + const std::string where = "block " + std::to_string(i) + " in " + path; + if (b.size[0] <= 0.0 || b.size[1] <= 0.0 || b.size[2] <= 0.0) + throw std::runtime_error("MuonShieldConfig: " + where + " has a non-positive size"); + + // Taper must not collapse (or invert) the downstream face. + const double farHalfX = 0.5 * b.size[0] + b.size[2] * std::tan(b.taper_deg[0] * kDegToRad); + const double farHalfY = 0.5 * b.size[1] + b.size[2] * std::tan(b.taper_deg[1] * kDegToRad); + if (farHalfX <= 0.0 || farHalfY <= 0.0) + throw std::runtime_error("MuonShieldConfig: " + where + + " has a taper that collapses its downstream face"); + + // The block must sit within the envelope in Z (light sanity check; + // rotations are not accounted for here). + if (b.start[2] < envStartMm - kEps || b.start[2] > envEndMm + kEps) + throw std::runtime_error("MuonShieldConfig: " + where + " has its upstream face (z = " + + std::to_string(b.start[2]) + " mm) outside the envelope"); + if (b.start[2] + b.size[2] > envEndMm + kEps) + throw std::runtime_error( + "MuonShieldConfig: " + where + " has its downstream face (z = " + + std::to_string(b.start[2] + b.size[2]) + " mm) beyond the envelope end"); + } + + return cfg; +} + +} // namespace SHiPGeometry diff --git a/subsystems/MuonShield/src/MuonShieldFactory.cpp b/subsystems/MuonShield/src/MuonShieldFactory.cpp index f77bd59..d1cf46e 100644 --- a/subsystems/MuonShield/src/MuonShieldFactory.cpp +++ b/subsystems/MuonShield/src/MuonShieldFactory.cpp @@ -3,182 +3,159 @@ #include "MuonShield/MuonShieldFactory.h" +#include "MuonShield/MuonShieldConfig.h" +#include "SHiPGeometry/ConfigPath.h" #include "SHiPGeometry/SHiPMaterials.h" #include #include #include #include +#include #include #include +#include +#include +#include #include +#include +#include +#include #include +// Absolute fallback path baked in by CMake so out-of-source builds always find +// MS.toml even when the CWD doesn't contain a copy of it. +#ifndef MS_TOML_DEFAULT_PATH +#define MS_TOML_DEFAULT_PATH "MS.toml" +#endif +// Install-time data directory path, set by CMake during install configuration. +#ifndef MS_TOML_INSTALL_PATH +#define MS_TOML_INSTALL_PATH "" +#endif + namespace SHiPGeometry { -// --------------------------------------------------------------------------- -// GDML-derived station data (all dimensions in mm) -// Bounding-box centres and half-sizes are computed from the arb8 vertices in -// ship_geometry.gdml. Station z-positions are relative to the MuonShieldArea -// container centre (GDML station_z_cm - 1676.33 cm) × 10 mm/cm. -// --------------------------------------------------------------------------- -const MuonShieldFactory::StationData MuonShieldFactory::k_stations[6] = { - // ── MagnAbsorb (GDML z = 319.5 cm, dz = 115.5 cm) ────────────────── - {"magn_absorb", - -13568.3 * mm, - 1020.0 * mm, - 1691.0 * mm, - 1155.0 * mm, - { - {250.0 * mm, 1690.0 * mm, 1155.0 * mm, 250.0 * mm, -1.0 * mm, "middle_mag_l"}, - {250.0 * mm, 1690.0 * mm, 1155.0 * mm, -250.0 * mm, 1.0 * mm, "middle_mag_r"}, - {250.0 * mm, 1690.0 * mm, 1155.0 * mm, 770.0 * mm, 0.0, "mag_ret_l"}, - {250.0 * mm, 1690.0 * mm, 1155.0 * mm, -770.0 * mm, 0.0, "mag_ret_r"}, - {510.0 * mm, 250.0 * mm, 1155.0 * mm, 510.0 * mm, 1440.0 * mm, "mag_top_left"}, - {510.0 * mm, 250.0 * mm, 1155.0 * mm, -510.0 * mm, 1440.0 * mm, "mag_top_right"}, - {510.0 * mm, 250.0 * mm, 1155.0 * mm, 510.0 * mm, -1440.0 * mm, "mag_bot_left"}, - {510.0 * mm, 250.0 * mm, 1155.0 * mm, -510.0 * mm, -1440.0 * mm, "mag_bot_right"}, - }}, - - // ── Magn1 (GDML z = 950 cm, dz = 495 cm) ──────────────────────────── - {"magn_1", - -7263.3 * mm, - 1697.0 * mm, - 1230.0 * mm, - 4950.0 * mm, - { - {399.6 * mm, 1228.2 * mm, 4950.0 * mm, 399.6 * mm, 0.0, "middle_mag_l"}, - {399.6 * mm, 1228.2 * mm, 4950.0 * mm, -399.6 * mm, 0.0, "middle_mag_r"}, - {487.7 * mm, 1229.2 * mm, 4950.0 * mm, 1208.7 * mm, 0.0, "mag_ret_l"}, - {487.7 * mm, 1229.2 * mm, 4950.0 * mm, -1208.7 * mm, 0.0, "mag_ret_r"}, - {848.2 * mm, 479.6 * mm, 4950.0 * mm, 848.2 * mm, 749.6 * mm, "mag_top_left"}, - {848.2 * mm, 479.6 * mm, 4950.0 * mm, -848.2 * mm, 749.6 * mm, "mag_top_right"}, - {848.2 * mm, 479.6 * mm, 4950.0 * mm, 848.2 * mm, -749.6 * mm, "mag_bot_left"}, - {848.2 * mm, 479.6 * mm, 4950.0 * mm, -848.2 * mm, -749.6 * mm, "mag_bot_right"}, - }}, - - // ── Magn2 (GDML z = 1735.48 cm, dz = 280.48 cm) ───────────────────── - {"magn_2", - 591.5 * mm, - 1736.0 * mm, - 1056.0 * mm, - 2804.8 * mm, - { - {265.6 * mm, 1054.6 * mm, 2804.8 * mm, 265.6 * mm, 0.0, "middle_mag_l"}, - {265.6 * mm, 1054.6 * mm, 2804.8 * mm, -265.6 * mm, 0.0, "middle_mag_r"}, - {594.7 * mm, 1055.6 * mm, 2804.8 * mm, 1140.3 * mm, 0.0, "mag_ret_l"}, - {594.7 * mm, 1055.6 * mm, 2804.8 * mm, -1140.3 * mm, 0.0, "mag_ret_r"}, - {867.5 * mm, 312.8 * mm, 2804.8 * mm, 867.5 * mm, 742.8 * mm, "mag_top_left"}, - {867.5 * mm, 312.8 * mm, 2804.8 * mm, -867.5 * mm, 742.8 * mm, "mag_top_right"}, - {867.5 * mm, 312.8 * mm, 2804.8 * mm, 867.5 * mm, -742.8 * mm, "mag_bot_left"}, - {867.5 * mm, 312.8 * mm, 2804.8 * mm, -867.5 * mm, -742.8 * mm, "mag_bot_right"}, - }}, - - // ── Magn3 (GDML z = 2258.49 cm, dz = 232.53 cm) ───────────────────── - {"magn_3", - 5821.6 * mm, - 1781.0 * mm, - 597.0 * mm, - 2325.3 * mm, - { - {18.4 * mm, 595.8 * mm, 2325.3 * mm, 23.4 * mm, 0.0, "middle_mag_l"}, - {18.4 * mm, 595.8 * mm, 2325.3 * mm, -23.4 * mm, 0.0, "middle_mag_r"}, - {849.3 * mm, 596.8 * mm, 2325.3 * mm, 931.6 * mm, 0.0, "mag_ret_l"}, - {849.3 * mm, 596.8 * mm, 2325.3 * mm, -931.6 * mm, 0.0, "mag_ret_r"}, - {888.0 * mm, 18.4 * mm, 2325.3 * mm, 893.0 * mm, 578.4 * mm, "mag_top_left"}, - {888.0 * mm, 18.4 * mm, 2325.3 * mm, -893.0 * mm, 578.4 * mm, "mag_top_right"}, - {888.0 * mm, 18.4 * mm, 2325.3 * mm, 893.0 * mm, -578.4 * mm, "mag_bot_left"}, - {888.0 * mm, 18.4 * mm, 2325.3 * mm, -893.0 * mm, -578.4 * mm, "mag_bot_right"}, - }}, - - // ── Magn4 (GDML z = 2586.02 cm, dz = 85 cm) ───────────────────────── - {"magn_4", - 9096.9 * mm, - 1797.0 * mm, - 1332.0 * mm, - 850.0 * mm, - { - {535.6 * mm, 1330.2 * mm, 850.0 * mm, 535.6 * mm, 0.0, "middle_mag_l"}, - {535.6 * mm, 1330.2 * mm, 850.0 * mm, -535.6 * mm, 0.0, "middle_mag_r"}, - {713.0 * mm, 1331.2 * mm, 850.0 * mm, 1083.0 * mm, 0.0, "mag_ret_l"}, - {713.0 * mm, 1331.2 * mm, 850.0 * mm, -1083.0 * mm, 0.0, "mag_ret_r"}, - {898.0 * mm, 385.6 * mm, 850.0 * mm, 898.0 * mm, 945.6 * mm, "mag_top_left"}, - {898.0 * mm, 385.6 * mm, 850.0 * mm, -898.0 * mm, 945.6 * mm, "mag_top_right"}, - {898.0 * mm, 385.6 * mm, 850.0 * mm, 898.0 * mm, -945.6 * mm, "mag_bot_left"}, - {898.0 * mm, 385.6 * mm, 850.0 * mm, -898.0 * mm, -945.6 * mm, "mag_bot_right"}, - }}, - - // ── Magn5 (GDML z = 2914.84 cm, dz = 233.82 cm) ───────────────────── - {"magn_5", - 12385.1 * mm, - 1808.0 * mm, - 960.0 * mm, - 2338.2 * mm, - { - {200.0 * mm, 959.0 * mm, 2338.2 * mm, 200.0 * mm, 0.0, "middle_mag_l"}, - {200.0 * mm, 959.0 * mm, 2338.2 * mm, -200.0 * mm, 0.0, "middle_mag_r"}, - {728.9 * mm, 960.0 * mm, 2338.2 * mm, 1079.2 * mm, 0.0, "mag_ret_l"}, - {728.9 * mm, 960.0 * mm, 2338.2 * mm, -1079.2 * mm, 0.0, "mag_ret_r"}, - {904.0 * mm, 200.0 * mm, 2338.2 * mm, 904.0 * mm, 760.0 * mm, "mag_top_left"}, - {904.0 * mm, 200.0 * mm, 2338.2 * mm, -904.0 * mm, 760.0 * mm, "mag_top_right"}, - {904.0 * mm, 200.0 * mm, 2338.2 * mm, 904.0 * mm, -760.0 * mm, "mag_bot_left"}, - {904.0 * mm, 200.0 * mm, 2338.2 * mm, -904.0 * mm, -760.0 * mm, "mag_bot_right"}, - }}, -}; - -// --------------------------------------------------------------------------- - -MuonShieldFactory::MuonShieldFactory(SHiPMaterials& materials) : m_materials(materials) {} - -GeoPhysVol* MuonShieldFactory::buildStation(const StationData& station) { - const GeoMaterial* air = m_materials.requireMaterial("Air"); - const GeoMaterial* iron = m_materials.requireMaterial("Iron"); - - // Air container that spans all 8 pieces of this station - auto* stationBox = - new GeoBox(station.containerHalfX, station.containerHalfY, station.containerHalfZ); - std::string containerName = "/SHiP/muon_shield/" + std::string(station.name); - auto* stationLog = new GeoLogVol(containerName, stationBox, air); - auto* stationPhys = new GeoPhysVol(stationLog); - - // Place 8 Iron bounding-box approximations - for (const PieceData& piece : station.pieces) { - auto* pieceBox = new GeoBox(piece.halfX, piece.halfY, piece.halfZ); - std::string pieceName = "/SHiP/muon_shield/" + std::string(station.name) + "/" + piece.name; - auto* pieceLog = new GeoLogVol(pieceName, pieceBox, iron); - auto* piecePhys = new GeoPhysVol(pieceLog); - - GeoTrf::Transform3D trf = GeoTrf::Translate3D(piece.centX, piece.centY, 0.0); - stationPhys->add(new GeoNameTag(pieceName)); - stationPhys->add(new GeoIdentifierTag(static_cast(&piece - &station.pieces[0]))); - stationPhys->add(new GeoTransform(trf)); - stationPhys->add(piecePhys); - } +// ── constructor / accessors ────────────────────────────────────────────────── + +MuonShieldFactory::MuonShieldFactory(SHiPMaterials& materials, std::string configPath) + : m_materials(materials), m_configPath(std::move(configPath)) {} + +void MuonShieldFactory::embedDaughter(GeoPhysVol* daughter, double worldCentreZ_mm, + const std::string& name) { + m_daughters.push_back({daughter, worldCentreZ_mm, name}); +} - return stationPhys; +void MuonShieldFactory::reserveSpace(const std::array& worldCentre_mm, + const std::array& size_mm, + const std::array& rotation_deg) { + m_reservations.push_back({worldCentre_mm, size_mm, rotation_deg}); } +std::string MuonShieldFactory::resolvedConfigPath() const { + return resolveConfigPath(m_configPath, MS_TOML_DEFAULT_PATH, MS_TOML_INSTALL_PATH); +} + +// ── build ──────────────────────────────────────────────────────────────────── + GeoPhysVol* MuonShieldFactory::build() { - const GeoMaterial* air = m_materials.requireMaterial("Air"); - - // Overall MuonShieldArea container (Air) - auto* areaBox = new GeoBox(s_areaHalfX, s_areaHalfY, s_areaHalfZ); - auto* areaLog = new GeoLogVol("/SHiP/muon_shield", areaBox, air); - auto* areaPhys = new GeoPhysVol(areaLog); - - // Build and place 6 stations - for (const StationData& station : k_stations) { - GeoPhysVol* stationPhys = buildStation(station); - std::string stationName = "/SHiP/muon_shield/" + std::string(station.name); - GeoTrf::Transform3D trf = GeoTrf::Translate3D(0.0, 0.0, station.stationZ); - areaPhys->add(new GeoNameTag(stationName)); - areaPhys->add(new GeoIdentifierTag(static_cast(&station - &k_stations[0]))); - areaPhys->add(new GeoTransform(trf)); - areaPhys->add(stationPhys); + const MuonShieldConfig cfg = readMuonShieldConfig(resolvedConfigPath()); + + GeoMaterial* air = m_materials.requireMaterial("Air"); + GeoMaterial* absorber = m_materials.requireMaterial(cfg.block_material); + + // Envelope-derived dimensions (mm). + const double halfZ = 0.5 * cfg.envelopeLengthZ_mm(); + m_centreZ_mm = cfg.envelopeCentreZ_mm(); + + // Air container spanning the full muon-shield envelope, centred on its own + // origin (SHiPGeometryBuilder places it at centreZ_mm()). + auto* containerBox = new GeoBox(cfg.envelope_half_x_mm, cfg.envelope_half_y_mm, halfZ); + auto* containerLog = new GeoLogVol("/SHiP/muon_shield", containerBox, air); + auto* containerPhys = new GeoPhysVol(containerLog); + + // Place each iron block. A block is anchored by its upstream (−z) face at + // its `start`, extends downstream by size[2], is rotated about `start` + // (extrinsic X→Y→Z), and may taper symmetrically into a GeoTrd. + int childId = 0; + for (const auto& b : cfg.blocks) { + const double halfX = 0.5 * b.size[0]; + const double halfY = 0.5 * b.size[1]; + const double halfLenZ = 0.5 * b.size[2]; + + const double farHalfX = halfX + b.size[2] * std::tan(b.taper_deg[0] * kDegToRad); + const double farHalfY = halfY + b.size[2] * std::tan(b.taper_deg[1] * kDegToRad); + + // Placement transform: shift the solid (centred at its origin) so its + // upstream face sits at `start`, rotate about `start` (extrinsic + // X→Y→Z), then translate into place. Container is at m_centreZ_mm in Z. + const double anchorX = b.start[0]; + const double anchorY = b.start[1]; + const double anchorZ = b.start[2] - m_centreZ_mm; + const GeoTrf::Transform3D rotation = GeoTrf::RotateZ3D(b.rotation_deg[2] * kDegToRad) * + GeoTrf::RotateY3D(b.rotation_deg[1] * kDegToRad) * + GeoTrf::RotateX3D(b.rotation_deg[0] * kDegToRad); + const GeoTrf::Transform3D trf = GeoTrf::Translate3D(anchorX, anchorY, anchorZ) * rotation * + GeoTrf::Translate3D(0.0, 0.0, halfLenZ); + + GeoShape* shape = nullptr; + if (b.taper_deg[0] == 0.0 && b.taper_deg[1] == 0.0) + shape = new GeoBox(halfX, halfY, halfLenZ); + else + shape = new GeoTrd(halfX, farHalfX, halfY, farHalfY, halfLenZ); + + // Carve out any reserved box that intersects this magnet (A - B). The + // reserved box, given in world coords, is expressed in the block's own + // (shape) frame before subtracting. + const GeoTrf::Vector3D blockCentre = trf.translation(); + const double blockMaxHX = std::max(halfX, farHalfX); + const double blockMaxHY = std::max(halfY, farHalfY); + for (const auto& r : m_reservations) { + const double rcx = r.centre_mm[0]; + const double rcy = r.centre_mm[1]; + const double rcz = r.centre_mm[2] - m_centreZ_mm; + const double rhx = 0.5 * r.size_mm[0]; + const double rhy = 0.5 * r.size_mm[1]; + const double rhz = 0.5 * r.size_mm[2]; + const bool intersects = std::abs(blockCentre.x() - rcx) < blockMaxHX + rhx && + std::abs(blockCentre.y() - rcy) < blockMaxHY + rhy && + std::abs(blockCentre.z() - rcz) < halfLenZ + rhz; + if (!intersects) + continue; + const GeoTrf::Transform3D resRot = GeoTrf::RotateZ3D(r.rotation_deg[2] * kDegToRad) * + GeoTrf::RotateY3D(r.rotation_deg[1] * kDegToRad) * + GeoTrf::RotateX3D(r.rotation_deg[0] * kDegToRad); + const GeoTrf::Transform3D resInBlock = + trf.inverse() * GeoTrf::Translate3D(rcx, rcy, rcz) * resRot; + auto* voidBox = new GeoBox(rhx, rhy, rhz); + shape = new GeoShapeSubtraction(shape, new GeoShapeShift(voidBox, resInBlock)); + } + + const std::string blockName = "/SHiP/muon_shield/block_" + std::to_string(childId); + auto* blockLog = new GeoLogVol(blockName, shape, absorber); + auto* blockPhys = new GeoPhysVol(blockLog); + + containerPhys->add(new GeoNameTag(blockName)); + containerPhys->add(new GeoIdentifierTag(childId++)); + containerPhys->add(new GeoTransform(trf)); + containerPhys->add(blockPhys); + } + + // Place the embedded daughters inside the container, keeping them as + // independent, named subsystems. (Callers must not list iron blocks over a + // daughter's footprint — the block list is the sole source of iron.) + for (const auto& d : m_daughters) { + const double localZ = d.worldCentreZ_mm - m_centreZ_mm; + if (std::abs(localZ) > halfZ) + throw std::runtime_error("MuonShieldFactory: embedded daughter '" + d.name + + "' lies outside the shield envelope in Z"); + containerPhys->add(new GeoNameTag(d.name)); + containerPhys->add(new GeoIdentifierTag(childId++)); + containerPhys->add(new GeoTransform(GeoTrf::Translate3D(0.0, 0.0, localZ))); + containerPhys->add(d.volume); } - return areaPhys; + return containerPhys; } } // namespace SHiPGeometry diff --git a/subsystems/MuonShield/test_muonshield.cpp b/subsystems/MuonShield/test_muonshield.cpp index 71dbdee..215786a 100644 --- a/subsystems/MuonShield/test_muonshield.cpp +++ b/subsystems/MuonShield/test_muonshield.cpp @@ -1,26 +1,183 @@ // SPDX-License-Identifier: LGPL-3.0-or-later // Copyright (C) CERN for the benefit of the SHiP Collaboration +#include "MuonShield/MuonShieldConfig.h" #include "MuonShield/MuonShieldFactory.h" #include "SHiPGeometry/SHiPMaterials.h" #include #include #include +#include +#include #include +#include +#include +#include +#include +#include +using SHiPGeometry::MuonShieldConfig; +using SHiPGeometry::MuonShieldFactory; +using SHiPGeometry::readMuonShieldConfig; using SHiPGeometry::SHiPMaterials; -// MuonShieldArea container halfX ≤ 2100 (CSV WARM max half-width), -// halfY ≤ 2300 (CSV WARM max half-height) -TEST_CASE("MuonShieldWithinEnvelope", "[muonshield]") { +namespace { +// Write a temp MS.toml with the given body and return its path. +std::string writeTempToml(const std::string& name, const std::string& body) { + std::ofstream out(name); + out << body; + out.close(); + return name; +} +} // namespace + +// Default MS.toml: the 7 FairShip TRY_2026 magnets (solid-block approximation) +// inside an auto-sized envelope (4.54–32.08 m, 1760 × 1320 mm half-sizes). +TEST_CASE("MuonShieldBuilds", "[muonshield]") { SHiPMaterials materials; - SHiPGeometry::MuonShieldFactory factory(materials); + MuonShieldFactory factory(materials); GeoPhysVol* ms = factory.build(); REQUIRE(ms != nullptr); + auto* box = dynamic_cast(ms->getLogVol()->getShape()); REQUIRE(box != nullptr); - CHECK(box->getXHalfLength() <= 2100.0); - CHECK(box->getYHalfLength() <= 2300.0); + CHECK_THAT(box->getXHalfLength(), Catch::Matchers::WithinAbs(1760.0, 1e-6)); + CHECK_THAT(box->getYHalfLength(), Catch::Matchers::WithinAbs(1320.0, 1e-6)); + CHECK_THAT(box->getZHalfLength(), Catch::Matchers::WithinAbs(13770.0, 1e-6)); + + // Envelope centre = (4.54 + 32.08)/2 m = 18.31 m. + CHECK_THAT(factory.centreZ_mm(), Catch::Matchers::WithinAbs(18310.0, 1e-6)); +} + +TEST_CASE("MuonShieldDefaultLayout", "[muonshield]") { + SHiPMaterials materials; + MuonShieldFactory factory(materials); + GeoPhysVol* ms = factory.build(); + REQUIRE(ms != nullptr); + // 7 solid magnets (the SND cavity is carved by reserveSpace, not here). + REQUIRE(ms->getNChildVols() == 7u); // NOLINT(readability/check) + + // Magnet 1: straight box, upstream face at z = 4.59 m, 2720 × 1600 × 3000 mm. + auto* block0 = dynamic_cast(ms->getChildVol(0)->getLogVol()->getShape()); + REQUIRE(block0 != nullptr); + CHECK_THAT(block0->getXHalfLength(), Catch::Matchers::WithinAbs(1360.0, 1e-6)); + CHECK_THAT(block0->getYHalfLength(), Catch::Matchers::WithinAbs(800.0, 1e-6)); + CHECK_THAT(block0->getZHalfLength(), Catch::Matchers::WithinAbs(1500.0, 1e-6)); + + // Anchor (upstream face) at world 4590 → local (4590-18310); centre + halfLen. + const double blockCentreLocalZ = (4590.0 - 18310.0) + 1500.0; // = -12220 + CHECK_THAT(ms->getXToChildVol(0).translation().z(), + Catch::Matchers::WithinAbs(blockCentreLocalZ, 1e-3)); + + // Magnet 6 (index 5) widens in Y → GeoTrd. + CHECK(dynamic_cast(ms->getChildVol(5)->getLogVol()->getShape()) != nullptr); +} + +TEST_CASE("MuonShieldReservationCarvesIron", "[muonshield]") { + // A reserved box (the SND envelope: 800 × 800 × 5100 mm at z = 28.95 m) is + // subtracted (A - B) from every magnet it intersects, leaving upstream + // magnets untouched. + SHiPMaterials materials; + MuonShieldFactory factory(materials); + factory.reserveSpace({0.0, 0.0, 28950.0}, {800.0, 800.0, 5100.0}); + GeoPhysVol* ms = factory.build(); + REQUIRE(ms != nullptr); + REQUIRE(ms->getNChildVols() == 7u); + + // Magnet 1 (upstream, far from the SND) is untouched → still a plain box. + CHECK(dynamic_cast(ms->getChildVol(0)->getLogVol()->getShape()) != nullptr); + // The last magnet intersects the reservation → its shape is a subtraction. + CHECK(dynamic_cast(ms->getChildVol(6)->getLogVol()->getShape()) != + nullptr); +} + +TEST_CASE("MuonShieldTaperMakesTrd", "[muonshield]") { + // A block with an X taper becomes a GeoTrd that widens downstream. + const std::string path = writeTempToml( + "MS_taper.toml", + "[[block]]\nstart = [0,0,12000]\nsize = [3000,2000,2000]\ntaper = [5.0, 0.0]\n"); + SHiPMaterials materials; + MuonShieldFactory factory(materials, path); + GeoPhysVol* ms = factory.build(); + REQUIRE(ms != nullptr); + REQUIRE(ms->getNChildVols() == 1u); + + auto* trd = dynamic_cast(ms->getChildVol(0)->getLogVol()->getShape()); + REQUIRE(trd != nullptr); + CHECK_THAT(trd->getXHalfLength1(), Catch::Matchers::WithinAbs(1500.0, 1e-6)); + const double farHalfX = 1500.0 + 2000.0 * std::tan(5.0 * 3.14159265358979323846 / 180.0); + CHECK_THAT(trd->getXHalfLength2(), Catch::Matchers::WithinAbs(farHalfX, 1e-3)); + CHECK_THAT(trd->getYHalfLength1(), Catch::Matchers::WithinAbs(1000.0, 1e-6)); + CHECK_THAT(trd->getYHalfLength2(), Catch::Matchers::WithinAbs(1000.0, 1e-6)); +} + +TEST_CASE("MuonShieldParsesRotation", "[muonshield]") { + // The optional `rotation` field is parsed into rotation_deg (degrees, about + // x, y, z). The factory applies it as an extrinsic X→Y→Z rotation about the + // block's start; asserting the built orientation would require reading the + // child transform back, which is left to OverlapCheck / a future test. + const std::string path = writeTempToml( + "MS_rot.toml", + "[[block]]\nstart = [0,0,12000]\nsize = [200,200,400]\nrotation = [10.0, 20.0, 30.0]\n"); + const MuonShieldConfig cfg = readMuonShieldConfig(path); + REQUIRE(cfg.blocks.size() == 1u); + CHECK_THAT(cfg.blocks[0].rotation_deg[0], Catch::Matchers::WithinAbs(10.0, 1e-9)); + CHECK_THAT(cfg.blocks[0].rotation_deg[1], Catch::Matchers::WithinAbs(20.0, 1e-9)); + CHECK_THAT(cfg.blocks[0].rotation_deg[2], Catch::Matchers::WithinAbs(30.0, 1e-9)); + + // A rotated block still builds and yields a valid shape. + SHiPMaterials materials; + MuonShieldFactory factory(materials, path); + GeoPhysVol* ms = factory.build(); + REQUIRE(ms != nullptr); + REQUIRE(ms->getNChildVols() == 1u); + CHECK(dynamic_cast(ms->getChildVol(0)->getLogVol()->getShape()) != nullptr); +} + +TEST_CASE("MuonShieldEmbedsDaughter", "[muonshield]") { + SHiPMaterials materials; + auto* dBox = new GeoBox(100.0, 100.0, 500.0); + auto* dLog = new GeoLogVol("/SHiP/dummy", dBox, materials.requireMaterial("Air")); + auto* dPhys = new GeoPhysVol(dLog); + + MuonShieldFactory factory(materials); // default 7 solid magnets + factory.embedDaughter(dPhys, 28.95 * 1000.0, "/SHiP/dummy"); + GeoPhysVol* ms = factory.build(); + REQUIRE(ms != nullptr); + // 7 iron magnets + the embedded daughter. + CHECK(ms->getNChildVols() == 8u); // NOLINT(readability/check) + + bool found = false; + for (unsigned i = 0; i < ms->getNChildVols(); ++i) + if (ms->getChildVol(i)->getLogVol()->getName() == "/SHiP/dummy") + found = true; + CHECK(found); +} + +TEST_CASE("MuonShieldRejectsNonPositiveSize", "[muonshield]") { + const std::string path = writeTempToml( + "MS_badsize.toml", "[[block]]\nstart = [0,0,12000]\nsize = [-3000,2000,2000]\n"); + CHECK_THROWS_AS(readMuonShieldConfig(path), std::runtime_error); +} + +TEST_CASE("MuonShieldRejectsCollapsingTaper", "[muonshield]") { + const std::string path = writeTempToml( + "MS_badtaper.toml", + "[[block]]\nstart = [0,0,12000]\nsize = [3000,2000,2000]\ntaper = [-45.0, 0.0]\n"); + CHECK_THROWS_AS(readMuonShieldConfig(path), std::runtime_error); +} + +TEST_CASE("MuonShieldRejectsBlockOutsideEnvelope", "[muonshield]") { + const std::string path = writeTempToml( + "MS_outside.toml", "[[block]]\nstart = [0,0,40000]\nsize = [3000,2000,1000]\n"); + CHECK_THROWS_AS(readMuonShieldConfig(path), std::runtime_error); +} + +TEST_CASE("MuonShieldEmptyBlockList", "[muonshield]") { + // No blocks → container with no iron (valid; e.g. a placeholder). + const std::string path = writeTempToml("MS_empty.toml", "block_material = \"Iron\"\n"); + MuonShieldConfig cfg = readMuonShieldConfig(path); + CHECK(cfg.blocks.empty()); } diff --git a/subsystems/NeutrinoDetector/CMakeLists.txt b/subsystems/NeutrinoDetector/CMakeLists.txt index 95ce231..a1b4c8a 100644 --- a/subsystems/NeutrinoDetector/CMakeLists.txt +++ b/subsystems/NeutrinoDetector/CMakeLists.txt @@ -1,4 +1,45 @@ # SPDX-License-Identifier: LGPL-3.0-or-later # Copyright (C) CERN for the benefit of the SHiP Collaboration +# Shared recipe: library (globs src/*.cpp, incl. SNDEnvelope.cpp), include dirs, +# kernel link, Catch test. ship_add_subsystem(NeutrinoDetector) + +# --- Extras: NeutrinoDetector parses SD.toml (the SND reservation envelope) via +# SNDEnvelope (toml++), the same way MuonShield/Calorimeter parse their configs. + +find_package(tomlplusplus REQUIRED) + +# toml++ is a build-time-only dependency (used only inside SNDEnvelope.cpp, never +# in a public header). BUILD_INTERFACE keeps it out of the install export set. +target_link_libraries( + NeutrinoDetector + PRIVATE $ +) + +# Stage SD.toml next to the tests and into the top-level build dir (for +# test_builder / OverlapCheck), so it is found when running from either. +configure_file( + ${CMAKE_CURRENT_SOURCE_DIR}/SD.toml + ${CMAKE_CURRENT_BINARY_DIR}/SD.toml + COPYONLY +) +configure_file( + ${CMAKE_CURRENT_SOURCE_DIR}/SD.toml + ${CMAKE_BINARY_DIR}/SD.toml + COPYONLY +) + +install( + FILES ${CMAKE_CURRENT_SOURCE_DIR}/SD.toml + DESTINATION ${CMAKE_INSTALL_DATADIR}/SHiPGeometry +) + +# Absolute fallbacks for resolveConfigPath(), so the parser finds SD.toml when +# the CWD does not have it (mirrors the calorimeter's CALO_TOML_* defines). +target_compile_definitions( + NeutrinoDetector + PRIVATE + SD_TOML_DEFAULT_PATH="${CMAKE_CURRENT_SOURCE_DIR}/SD.toml" + SD_TOML_INSTALL_PATH="${CMAKE_INSTALL_PREFIX}/${CMAKE_INSTALL_DATADIR}/SHiPGeometry/SD.toml" +) diff --git a/subsystems/NeutrinoDetector/README.md b/subsystems/NeutrinoDetector/README.md index 99ca33c..1d1d0ef 100644 --- a/subsystems/NeutrinoDetector/README.md +++ b/subsystems/NeutrinoDetector/README.md @@ -35,7 +35,7 @@ bounded for navigation. ## Geometry tree ```text -/SHiP/neutrino_detector (Air, 640 x 640 x 4100 mm box) +/SHiP/neutrino_detector (Air, 800 x 800 x 5100 mm box) |- veto/P{0..2}_B{0..6} 21 PVT bars (3 planes x 7) |- target/L{0..119}/{W,Si_X,Si_Y} 120 W + 120 Si-X + 120 Si-Y \- hcal/S{0..2}_L{0..13}/ @@ -48,11 +48,15 @@ bounded for navigation. Direct container children: 21 + 360 + 4438 = 4819. Total volumes including the individual fibres: 340819 (~336k fibres). -Position in world: z = 28 950 mm (centre of the 26.40-31.50 m WARM SND slot). -The SND is part of the muon shield, so its box container sits within the -downstream end of the muon-shield region; the resulting envelope overlap is -declared in `tests/test_consistency.cpp`, analogous to the trackers/magnet -overlap. +Position: z = 28 950 mm (centre of the 26.40-31.50 m WARM SND slot). The SND is +part of the muon shield: in volume terms its container is a **daughter of the +muon-shield container**, nested by `SHiPGeometryBuilder` via +`MuonShieldFactory::embedDaughter`, which only places the pre-built daughter in +the shield's volume tree. The matching cavity in the iron is carved separately +by `MuonShieldFactory::reserveSpace` (a Boolean subtraction of the SND +reservation envelope from `SD.toml`), so the two do not overlap. It remains an +independent subsystem with its own factory, config, and +`/SHiP/neutrino_detector` naming — only its place in the volume tree changes. ## Materials diff --git a/subsystems/NeutrinoDetector/SD.toml b/subsystems/NeutrinoDetector/SD.toml new file mode 100644 index 0000000..ae6ba9a --- /dev/null +++ b/subsystems/NeutrinoDetector/SD.toml @@ -0,0 +1,16 @@ +# SPDX-License-Identifier: LGPL-3.0-or-later +# Copyright (C) CERN for the benefit of the SHiP Collaboration +# +# Scattering & Neutrino Detector (SND) reservation envelope. +# +# This declares the box that is subtracted from the muon-shield iron (A - B) to +# make room for the SND, so the two subsystems can be defined independently and +# merged at build time. For now this is only the reservation envelope; the SND's +# internal geometry will move into this file in a later, dedicated PR. +# +# Values are in world/beamline coordinates (target-front-face origin), mm. +# Defaults match the current SND container (800 x 800 x 5100 mm at z = 28.95 m). + +centre = [0.0, 0.0, 28950.0] # world centre of the reservation box +size = [800.0, 800.0, 5100.0] # full x, y, z +# rotation = [0.0, 0.0, 0.0] # optional, degrees, extrinsic X -> Y -> Z about centre diff --git a/subsystems/NeutrinoDetector/include/NeutrinoDetector/SNDEnvelope.h b/subsystems/NeutrinoDetector/include/NeutrinoDetector/SNDEnvelope.h new file mode 100644 index 0000000..cb1f26b --- /dev/null +++ b/subsystems/NeutrinoDetector/include/NeutrinoDetector/SNDEnvelope.h @@ -0,0 +1,39 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) CERN for the benefit of the SHiP Collaboration + +#pragma once + +#include +#include + +namespace SHiPGeometry { + +/** + * @brief The SND reservation envelope: a box subtracted from the muon-shield + * iron to make room for the neutrino detector. + * + * All values are in world/beamline coordinates (target-front-face origin). + * Populated by readSNDEnvelope() from SD.toml. + */ +struct SNDEnvelope { + std::array centre_mm = {0.0, 0.0, 28950.0}; ///< world centre (mm) + std::array size_mm = {800.0, 800.0, 5100.0}; ///< full x, y, z (mm) + std::array rotation_deg = {0.0, 0.0, 0.0}; ///< extrinsic X->Y->Z (deg) +}; + +/** + * @brief Parse an SD.toml file into an SNDEnvelope. + * + * The path is resolved with the shared resolveConfigPath() fallback chain + * (CWD, then the source tree, then the installed data dir), so callers can + * simply use the default. + * + * @throws std::runtime_error on malformed TOML, a wrong-length array, or a + * non-positive size. + */ +SNDEnvelope readSNDEnvelope(const std::string& path = "SD.toml"); + +/// Resolve the SD.toml path that readSNDEnvelope() would actually open. +std::string resolvedSNDConfigPath(const std::string& path = "SD.toml"); + +} // namespace SHiPGeometry diff --git a/subsystems/NeutrinoDetector/src/SNDEnvelope.cpp b/subsystems/NeutrinoDetector/src/SNDEnvelope.cpp new file mode 100644 index 0000000..f2c58d9 --- /dev/null +++ b/subsystems/NeutrinoDetector/src/SNDEnvelope.cpp @@ -0,0 +1,95 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) CERN for the benefit of the SHiP Collaboration + +#include "NeutrinoDetector/SNDEnvelope.h" + +#include "SHiPGeometry/ConfigPath.h" + +#include +#include +#include +#include +#include +#include +#include + +// Absolute fallback paths baked in by CMake so out-of-source builds always find +// SD.toml even when the CWD doesn't contain a copy of it. +#ifndef SD_TOML_DEFAULT_PATH +#define SD_TOML_DEFAULT_PATH "SD.toml" +#endif +#ifndef SD_TOML_INSTALL_PATH +#define SD_TOML_INSTALL_PATH "" +#endif + +namespace SHiPGeometry { + +namespace { + +double toNumber(const toml::node* node, const std::string& key, const std::string& path) { + if (node) { + if (auto d = node->value()) + return *d; + if (auto i = node->value()) + return static_cast(*i); + } + throw std::runtime_error("SNDEnvelope: '" + key + "' values must be numbers in " + path); +} + +std::array readVec3(const toml::table& table, const char* key, const std::string& path, + bool required, const std::array& fallback) { + auto node = table[key]; + if (!node) { + if (required) + throw std::runtime_error("SNDEnvelope: missing required '" + std::string(key) + + "' in " + path); + return fallback; + } + const toml::array* arr = node.as_array(); + if (!arr || arr->size() != 3) + throw std::runtime_error("SNDEnvelope: '" + std::string(key) + + "' must be a 3-element array in " + path); + std::array out{}; + for (std::size_t i = 0; i < 3; ++i) + out[i] = toNumber(arr->get(i), key, path); + return out; +} + +} // namespace + +std::string resolvedSNDConfigPath(const std::string& path) { + return resolveConfigPath(path, SD_TOML_DEFAULT_PATH, SD_TOML_INSTALL_PATH); +} + +SNDEnvelope readSNDEnvelope(const std::string& rawPath) { + const std::string path = resolvedSNDConfigPath(rawPath); + toml::table table; + try { + table = toml::parse_file(path); + } catch (const toml::parse_error& e) { + throw std::runtime_error("SNDEnvelope: failed to parse " + path + ": " + + std::string(e.description())); + } + + SNDEnvelope env; + + // Warn about unknown top-level keys (typos / stale fields), like + // readMuonShieldConfig. Only centre, size and rotation are recognised. + for (const auto& [k, _] : table) { + const std::string_view key{k}; + if (key != "centre" && key != "size" && key != "rotation") + std::cerr << "SNDEnvelope: warning: unknown key '" << k << "' in " << path + << " (ignored)\n"; + } + + env.centre_mm = readVec3(table, "centre", path, true, env.centre_mm); + env.size_mm = readVec3(table, "size", path, true, env.size_mm); + env.rotation_deg = readVec3(table, "rotation", path, false, env.rotation_deg); + + if (env.size_mm[0] <= 0.0 || env.size_mm[1] <= 0.0 || env.size_mm[2] <= 0.0) + throw std::runtime_error("SNDEnvelope: size must be positive in " + path); + + return env; +} + +} // namespace SHiPGeometry diff --git a/tests/test_consistency.cpp b/tests/test_consistency.cpp index ad8130c..4497ddd 100644 --- a/tests/test_consistency.cpp +++ b/tests/test_consistency.cpp @@ -76,9 +76,10 @@ TEST_CASE("ConsistencyTest.ExpectedSubsystemCount", "[consistency]") { REQUIRE(world != nullptr); auto subsystems = collectSubsystems(world); - // 9 subsystems: target, muon_shield, neutrino_detector, upstream_tagger, - // decay_volume, trackers, magnet, timing_detector, calorimeter - CHECK(subsystems.size() == 9u); // NOLINT(readability/check) + // 8 subsystems: target, muon_shield, upstream_tagger, decay_volume, + // trackers, magnet, timing_detector, calorimeter. The neutrino detector is + // nested inside the muon-shield container, so it is not a direct world child. + CHECK(subsystems.size() == 8u); // NOLINT(readability/check) } TEST_CASE("ConsistencyTest.SubsystemsGenerallyInZOrder", "[consistency]") { @@ -110,13 +111,9 @@ TEST_CASE("ConsistencyTest.NoUnexpectedZOverlaps", "[consistency]") { // The trackers container intentionally spans across the magnet // (stations 1-2 before, stations 3-4 after), so that pair is allowed to overlap. - // The SND sits inside the downstream end of the muon-shield region, so the - // muon_shield / neutrino_detector pair is also an intentional overlap. auto isAllowedOverlap = [](const std::string& a, const std::string& b) { return (a == "/SHiP/trackers" && b == "/SHiP/magnet") || - (a == "/SHiP/magnet" && b == "/SHiP/trackers") || - (a == "/SHiP/muon_shield" && b == "/SHiP/neutrino_detector") || - (a == "/SHiP/neutrino_detector" && b == "/SHiP/muon_shield"); + (a == "/SHiP/magnet" && b == "/SHiP/trackers"); }; // Sort by Z centre @@ -156,8 +153,10 @@ TEST_CASE("ConsistencyTest.PositionsSanity", "[consistency]") { // Centres as placed in SHiPGeometryBuilder::build() std::vector expected = { {"/SHiP/target", 432.5, 500.0}, - {"/SHiP/muon_shield", 16763.3, 500.0}, - {"/SHiP/neutrino_detector", 28950.0, 500.0}, + {"/SHiP/muon_shield", 18310.0, 500.0}, + // The neutrino detector is nested inside the muon shield (see + // MuonShieldFactory::embedDaughter), not a direct world child, so it is + // not among collectSubsystems(world). {"/SHiP/upstream_tagger", 32720.0, 500.0}, {"/SHiP/decay_volume", 58120.0, 500.0}, {"/SHiP/trackers", 89570.0, 500.0},