SBTConstants.h

Go to the documentation of this file.

// SPDX-License-Identifier: LGPL-3.0-or-later
// Copyright (C) CERN for the benefit of the SHiP Collaboration
 
#pragma once
 
#include <algorithm>
#include <string_view>
 
namespace SHiPGeometry::SBT {
 
struct SBTParams {
    // ── Frustum envelope ────────────────────────────────────────────────
    double xHalfEntrance = 1000.0;
    double yHalfEntrance = 1500.0;
    double xHalfExit = 2000.0;
    double yHalfExit = 3000.0;
    double totalLength = 50000.0;
    int nSubFrustum = 10;
    double yFloor = -3000.0;
    double zEntrance = -25000.0;
 
    // ── H-beam cross-section (HEA 260 approximation) ────────────────────
    double hbeamHeight = 250.0;
    double hbeamFlangeWidth = 260.0;
    double hbeamFlangeThickness = 12.5;
    double hbeamWebThickness = 7.5;
 
    // ── Sensor containers / cells ───────────────────────────────────────
    double containerThickness = 225.0;
    double cellWallThickness = 5.0;
    int nCells = 6;
    double sensorClearance = 1.0;
 
    // ── Helium decay region ─────────────────────────────────────────────
    // 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.
    double heliumClearance = 0.001;
};
 
inline constexpr SBTParams kSBT{};
 
// ── Air container enclosing the SBT structure + sensors and helium ──────
// The experiment's fixed envelope allocation for the decay region. Not an SBT
// dimension and deliberately not part of SBTParams: nothing is derived from
// it, and the three static_asserts at the bottom of this header are bounds the
// SBT must fit *inside*, not parameters it is built *from*.
inline constexpr double kEnvelopeHalfX = 2200.0;
inline constexpr double kEnvelopeHalfY = 3300.0;
inline constexpr double kEnvelopeHalfZ = 25200.0;
 
// ── Materials (names in SHiPMaterials) ──────────────────────────────────
inline constexpr std::string_view kAirMaterial = "Air";                 // envelope container
inline constexpr std::string_view kSteelMaterial = "Iron";              // H-beam structure
inline constexpr std::string_view kWallMaterial = "Aluminium";          // sensor cell walls
inline constexpr std::string_view kScintillatorMaterial = "LAB";        // sensor cells
inline constexpr std::string_view kHeliumMaterial = "PressurisedHe90";  // decay region
 
// ── Derived quantities ──────────────────────────────────────────────────
//
//  Convention for everything below: the parameters come last and default to
//  kSBT. Last rather than first because most of these already have a
//  meaningful leading argument (a Z, a half-width) and a defaulted parameter
//  cannot precede a non-defaulted one. The parameter is named `params`, not
//  `p`, because `p` is the conventional loop variable for a HeliumPiece in
//  SBTEnvelope and in the tests.
 
constexpr double subLength(const SBTParams& params = kSBT) {
    return params.totalLength / params.nSubFrustum;
}
constexpr double webHeight(const SBTParams& params = kSBT) {
    return params.hbeamHeight - 2.0 * params.hbeamFlangeThickness;
}
constexpr int nWalls(const SBTParams& params = kSBT) {
    return params.nCells + 1;
}
 
// ── Frustum profile ─────────────────────────────────────────────────────
constexpr double zExit(const SBTParams& params = kSBT) {
    return params.zEntrance + params.totalLength;
}
constexpr double xGrowth(const SBTParams& params = kSBT) {
    return (params.xHalfExit - params.xHalfEntrance) / params.totalLength;
}
constexpr double yGrowth(const SBTParams& params = kSBT) {
    return (params.yHalfExit - params.yHalfEntrance) / params.totalLength;
}
constexpr double xHalfAt(double z_mm, const SBTParams& params = kSBT) {
    return params.xHalfEntrance + (z_mm - params.zEntrance) * xGrowth(params);
}
constexpr double yHalfAt(double z_mm, const SBTParams& params = kSBT) {
    return params.yHalfEntrance + (z_mm - params.zEntrance) * yGrowth(params);
}
constexpr double zSubLo(int s, const SBTParams& params = kSBT) {
    return params.zEntrance + s * subLength(params);
}
 
// ── H-beam cross-section primitives ─────────────────────────────────────
constexpr double hbeamFlangeOffset(const SBTParams& params = kSBT) {
    return 0.5 * params.hbeamHeight - 0.5 * params.hbeamFlangeThickness;
}
constexpr double hbeamHalfHeight(const SBTParams& params = kSBT) {
    return 0.5 * params.hbeamHeight;
}
 
// ══════════════════════════════════════════════════════════════════════
//  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.
// ══════════════════════════════════════════════════════════════════════
 
constexpr double zSplitOffset(const SBTParams& params = kSBT) {
    return 0.5 * params.hbeamHeight + 0.5 * params.hbeamFlangeThickness;
}
 
// --- Side (±X) scintillator containers -------------------------------
constexpr double sideContainerHalfThickness(const SBTParams& params = kSBT) {
    return 0.5 * params.containerThickness;
}
constexpr double sideContainerCentreX(double x_half_mm, const SBTParams& params = kSBT) {
    return x_half_mm - 0.5 * params.hbeamFlangeWidth - sideContainerHalfThickness(params);
}
constexpr double sideSensorInnerX(double x_half_mm, const SBTParams& params = kSBT) {
    return sideContainerCentreX(x_half_mm, params) - sideContainerHalfThickness(params);
}
 
// --- Top/bottom (±Y) scintillator containers -------------------------
constexpr double topBottomContainerHalfThickness(const SBTParams& params = kSBT) {
    return 0.5 * params.containerThickness - params.sensorClearance;
}
constexpr double topBottomContainerCentreY(double y_half_mm, const SBTParams& params = kSBT) {
    return y_half_mm - 0.5 * params.containerThickness;
}
constexpr double topBottomSensorInnerY(double y_half_mm, const SBTParams& params = kSBT) {
    return topBottomContainerCentreY(y_half_mm, params) - topBottomContainerHalfThickness(params);
}
constexpr double topBottomAvailX(double x_half_mm, const SBTParams& params = kSBT) {
    return x_half_mm - 0.5 * params.hbeamFlangeWidth - 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.
constexpr double longBeamCentreY(double y_half_mm, const SBTParams& params = kSBT) {
    return y_half_mm - 0.5 * webHeight(params);
}
 
namespace detail {
 
constexpr double sqrtConst(double v) {
    if (v <= 0.0) {
        return 0.0;
    }
    double x = v;
    double prev = 0.0;
    for (int i = 0; i < 200 && x != prev; ++i) {
        prev = x;
        x = 0.5 * (x + v / x);
    }
    return x;
}
 
constexpr double longBeamTaper(const SBTParams& params = kSBT) {
    const double g = yGrowth(params);
    return sqrtConst(1.0 + g * g);
}
 
}  // namespace detail
 
constexpr double longBeamInnerY(double y_half_mm, const SBTParams& params = kSBT) {
    return longBeamCentreY(y_half_mm, params) -
           hbeamHalfHeight(params) * detail::longBeamTaper(params);
}
 
// ── Compile-time validation ─────────────────────────────────────────────
//
//  Whether an arbitrary SBTParams is buildable at all is SBTEnvelope's
//  leavesDecayRegion(), which the tests exercise on configurations that must
//  not build. What follows constrains the *shipped* parameters only.
 
constexpr bool isWellFormed(const SBTParams& params = kSBT) {
    return params.totalLength > 0.0 && params.xHalfEntrance > 0.0 && params.yHalfEntrance > 0.0 &&
           params.xHalfExit > 0.0 && params.yHalfExit > 0.0 && params.hbeamFlangeWidth > 0.0 &&
           params.hbeamFlangeThickness > 0.0 && params.hbeamWebThickness > 0.0 &&
           params.containerThickness > 0.0 && params.cellWallThickness > 0.0 &&
           params.nSubFrustum > 0 && params.nCells > 0 && params.sensorClearance > 0.0 &&
           params.sensorClearance < 0.5 * params.containerThickness &&
           params.hbeamHeight > 2.0 * params.hbeamFlangeThickness && params.heliumClearance >= 0.0;
}
 
static_assert(isWellFormed(kSBT),
              "the shipped SBT parameters are not self-consistent: check positivity, the "
              "sensor clearance against the container thickness, and the flange thickness "
              "against the beam height");
static_assert(zSplitOffset(kSBT) < subLength(kSBT),
              "the sensor containers' flat piece no longer fits inside a sub-frustum; "
              "reduce nSubFrustum or hbeamHeight");
 
// The SBT structure must fit the fixed decay-volume envelope. Bounds on the
// outermost structure reach; with the current values the reaches are
// ~2130 / ~3259 / ~25130 mm, inside the 2200 / 3300 / 25200 envelope.
// X: vertical columns sit at the wider of the two frustum faces, with a
//    half-flange overhang.
// Y: top/bottom cross-beams are shifted a full beam-height above yHalfExit
//    (plus a small frustum-growth term); +10 pads the assembly standoff.
// Z: the structure spans [zEntrance, zExit()]; cross-beam flanges extend
//    half a flange-width beyond the end rows.
//
// These bound kSBT against the envelope allocation and are deliberately not
// part of isWellFormed(): the parameter sweep in test_decayvolume varies the
// frustum well past this allocation on purpose, and builds into a generous
// throwaway container precisely because of it.
static_assert(std::max(kSBT.xHalfEntrance, kSBT.xHalfExit) + 0.5 * kSBT.hbeamFlangeWidth <=
                  kEnvelopeHalfX,
              "SBT structure pierces the decay-volume envelope in X");
static_assert(kSBT.yHalfExit + kSBT.hbeamHeight + 0.5 * kSBT.hbeamFlangeWidth * yGrowth(kSBT) +
                      10.0 <=
                  kEnvelopeHalfY,
              "SBT structure pierces the decay-volume envelope in Y");
static_assert(std::max(kSBT.zEntrance < 0.0 ? -kSBT.zEntrance : kSBT.zEntrance,
                       zExit(kSBT) < 0.0 ? -zExit(kSBT) : zExit(kSBT)) +
                      0.5 * kSBT.hbeamFlangeWidth <=
                  kEnvelopeHalfZ,
              "SBT structure pierces the decay-volume envelope in Z");
 
}  // namespace SHiPGeometry::SBT