SBTEnvelope.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 "DecayVolume/SBTConstants.h"
 
#include <algorithm>
#include <array>
#include <cstddef>
#include <span>
 
class GeoMaterial;
class GeoPhysVol;
 
namespace SHiPGeometry::SBT {
 
struct HeliumPiece {
    double z_lo_mm = 0.0;
    double z_hi_mm = 0.0;
    double dx_lo_mm = 0.0;  
    double dx_hi_mm = 0.0;  
    double dy_lo_mm = 0.0;  
    double dy_hi_mm = 0.0;  
};
 
namespace detail {
 
// Index of the sub-frustum containing z (clamped at the exit face).
//
// The cast truncates towards zero rather than flooring, which lands on the
// same index for every f: the two agree for f >= 0, and a negative f clamps
// to 0 whichever way it was rounded. Truncating also keeps this header clear
// of C++23's constexpr <cmath>, which not every standard library ships yet.
constexpr int subFrustumAt(double z_mm, const SBTParams& params = kSBT) {
    const double f = (z_mm - params.zEntrance) / subLength(params);
    return std::clamp(static_cast<int>(f), 0, params.nSubFrustum - 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.
constexpr double sideTrackedXHalf(double z_mm, const SBTParams& params = kSBT) {
    const int s = subFrustumAt(z_mm, params);
    const double zLo = zSubLo(s, params);
    if (z_mm <= zLo + zSplitOffset(params)) {
        return xHalfAt(zLo, params);
    }
    return xHalfAt(z_mm, params);
}
 
}  // namespace detail
 
constexpr double innerFreeHalfX(double z_mm, const SBTParams& params = kSBT) {
    // 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 sideSensorInnerX(detail::sideTrackedXHalf(z_mm, params), params);
}
 
constexpr double innerFreeHalfY(double z_mm, const SBTParams& params = kSBT) {
    const double yHalf = yHalfAt(z_mm, params);
    // Top/bottom scintillator containers ...
    const double sensor = topBottomSensorInnerY(yHalf, params);
    // ... and the longitudinal beams' inner flange, which hangs below them.
    const double beam = longBeamInnerY(yHalf, params);
    // (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);
}
 
namespace detail {
 
constexpr double knotAt(int i, const SBTParams& params = kSBT) {
    if (i >= 2 * params.nSubFrustum) {
        return zExit(params);
    }
    const int s = i / 2;
    return (i % 2 == 0) ? zSubLo(s, params) : zSubLo(s, params) + zSplitOffset(params);
}
 
// 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.
constexpr double envelopeAtKnot(double z_mm, bool isX, const SBTParams& params = kSBT) {
    constexpr double kEps = 1e-6;
    const double lo = std::max(z_mm - kEps, params.zEntrance);
    const double hi = std::min(z_mm + kEps, zExit(params));
    return isX ? std::min(innerFreeHalfX(lo, params), innerFreeHalfX(hi, params))
               : std::min(innerFreeHalfY(lo, params), innerFreeHalfY(hi, params));
}
 
}  // namespace detail
 
constexpr HeliumPiece heliumPieceAt(int i, const SBTParams& params = kSBT) {
    const double zLo = detail::knotAt(i, params);
    const double zHi = detail::knotAt(i + 1, params);
    return HeliumPiece{
        .z_lo_mm = zLo,
        .z_hi_mm = zHi,
        .dx_lo_mm = detail::envelopeAtKnot(zLo, /*isX=*/true, params) - params.heliumClearance,
        .dx_hi_mm = detail::envelopeAtKnot(zHi, /*isX=*/true, params) - params.heliumClearance,
        .dy_lo_mm = detail::envelopeAtKnot(zLo, /*isX=*/false, params) - params.heliumClearance,
        .dy_hi_mm = detail::envelopeAtKnot(zHi, /*isX=*/false, params) - params.heliumClearance,
    };
}
 
constexpr bool leavesDecayRegion(const SBTParams& params = kSBT) {
    if (!isWellFormed(params)) {
        return false;
    }
    for (int i = 0; i < 2 * params.nSubFrustum; ++i) {
        const HeliumPiece piece = heliumPieceAt(i, params);
        if (piece.z_hi_mm <= piece.z_lo_mm) {
            return false;
        }
        if (piece.dx_lo_mm <= 0.0 || piece.dx_hi_mm <= 0.0 || piece.dy_lo_mm <= 0.0 ||
            piece.dy_hi_mm <= 0.0) {
            return false;
        }
    }
    return true;
}
 
inline constexpr auto kHeliumPieces = [] {
    std::array<HeliumPiece, 2 * kSBT.nSubFrustum> pieces{};
    for (int i = 0; i < 2 * kSBT.nSubFrustum; ++i) {
        pieces[static_cast<std::size_t>(i)] = heliumPieceAt(i, kSBT);
    }
    return pieces;
}();
 
// An SBT whose containers, beams or clearances swallow the decay region fails
// the build.
static_assert(leavesDecayRegion(kSBT),
              "the configured SBT leaves no decay region — check containerThickness, "
              "hbeamHeight, nSubFrustum and heliumClearance");
 
// ── Arbitrary parameters, for the tests ─────────────────────────────────
//
//  std::array's extent is part of its type, and test_decayvolume sweeps
//  nSubFrustum between 5 and 20, so a swept stack needs one type that covers
//  them all. This is a testing affordance, not a geometry parameter: the
//  shipped stack is kHeliumPieces above and is exactly sized.
 
inline constexpr int kMaxSubFrustum = 64;
static_assert(kSBT.nSubFrustum <= kMaxSubFrustum);
 
struct HeliumSlabs {
    std::array<HeliumPiece, 2u * kMaxSubFrustum> pieces{};
    std::size_t count = 0;
    [[nodiscard]] constexpr std::span<const HeliumPiece> view() const {
        return {pieces.data(), count};
    }
};
 
constexpr HeliumSlabs heliumSlabs(const SBTParams& params = kSBT) {
    HeliumSlabs slabs;
    if (params.nSubFrustum <= 0 || params.nSubFrustum > kMaxSubFrustum) {
        return slabs;
    }
    slabs.count = 2u * static_cast<std::size_t>(params.nSubFrustum);
    for (std::size_t i = 0; i < slabs.count; ++i) {
        slabs.pieces[i] = heliumPieceAt(static_cast<int>(i), params);
    }
    return slabs;
}
 
void buildHelium(GeoPhysVol* container, const GeoMaterial* helium,
                 std::span<const HeliumPiece> pieces = kHeliumPieces);
 
}  // namespace SHiPGeometry::SBT