src/Geometry/Extent.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Geometry/Extent.hpp"
0010
0011 #include "Acts/Utilities/VectorHelpers.hpp"
0012
0013 #include <algorithm>
0014 #include <cmath>
0015 #include <cstddef>
0016 #include <iomanip>
0017 #include <limits>
0018 #include <numbers>
0019
0020 Acts::Extent::Extent(const ExtentEnvelope& envelope)
0021     : m_constrains(0), m_envelope(envelope) {
0022   m_range[toUnderlying(AxisDirection::AxisR)] =
0023       Range1D<double>(0., std::numeric_limits<double>::max());
0024   m_range[toUnderlying(AxisDirection::AxisPhi)] =
0025       Range1D<double>(-std::numbers::pi, std::numbers::pi);
0026   m_range[toUnderlying(AxisDirection::AxisRPhi)] =
0027       Range1D<double>(0., std::numeric_limits<double>::max());
0028   m_range[toUnderlying(AxisDirection::AxisMag)] =
0029       Range1D<double>(0., std::numeric_limits<double>::max());
0030 }
0031
0032 void Acts::Extent::extend(const Vector3& vtx,
0033                           const std::vector<AxisDirection>& aDirs,
0034                           bool applyEnv, bool fillHistograms) {
0035   for (auto aDir : aDirs) {
0036     // Get the casted value given the binning value description
0037     double cValue = VectorHelpers::cast(vtx, aDir);
0038     if (fillHistograms) {
0039       m_valueHistograms[toUnderlying(aDir)].push_back(cValue);
0040     }
0041     // Apply envelope as suggested
0042     double lEnv = applyEnv ? m_envelope[aDir][0] : 0.;
0043     double hEnv = applyEnv ? m_envelope[aDir][1] : 0.;
0044     double mValue = cValue - lEnv;
0045     // Special protection for radial value
0046     if (aDir == AxisDirection::AxisR && mValue < 0.) {
0047       mValue = std::max(mValue, 0.);
0048     }
0049     if (constrains(aDir)) {
0050       m_range[toUnderlying(aDir)].expand(mValue, cValue + hEnv);
0051     } else {
0052       m_range[toUnderlying(aDir)].shrink(mValue, cValue + hEnv);
0053     }
0054     m_constrains.set(toUnderlying(aDir));
0055   }
0056 }
0057
0058 void Acts::Extent::extend(const Extent& rhs,
0059                           const std::vector<AxisDirection>& aDirs,
0060                           bool applyEnv) {
0061   for (auto aDir : aDirs) {
0062     // The value is constraint, envelope can be optional
0063     if (rhs.constrains(aDir)) {
0064       double lEnv = applyEnv ? m_envelope[aDir][0] : 0.;
0065       double hEnv = applyEnv ? m_envelope[aDir][1] : 0.;
0066       if (constrains(aDir)) {
0067         m_range[toUnderlying(aDir)].expand(
0068             rhs.range()[toUnderlying(aDir)].min() - lEnv,
0069             rhs.range()[toUnderlying(aDir)].max() + hEnv);
0070       } else {
0071         m_range[toUnderlying(aDir)].shrink(
0072             rhs.range()[toUnderlying(aDir)].min() - lEnv,
0073             rhs.range()[toUnderlying(aDir)].max() + hEnv);
0074       }
0075       m_constrains.set(toUnderlying(aDir));
0076     } else if (rhs.envelope()[aDir] != zeroEnvelope) {
0077       // Only an envelope given, but value is not constraint -> apply envelope
0078       m_range[toUnderlying(aDir)].expand(
0079           m_range[toUnderlying(aDir)].min() - rhs.envelope()[aDir][0],
0080           m_range[toUnderlying(aDir)].max() + rhs.envelope()[aDir][1]);
0081       m_constrains.set(toUnderlying(aDir));
0082     }
0083   }
0084 }
0085
0086 void Acts::Extent::addConstrain(const Acts::Extent& rhs,
0087                                 const ExtentEnvelope& envelope) {
0088   for (const auto& aDir : allAxisDirections()) {
0089     if (rhs.constrains(aDir) && !constrains(aDir)) {
0090       const auto& cRange = rhs.range(aDir);
0091       m_range[toUnderlying(aDir)].setMin(cRange.min() - envelope[aDir][0u]);
0092       m_range[toUnderlying(aDir)].setMax(cRange.max() + envelope[aDir][1u]);
0093       m_constrains.set(toUnderlying(aDir));
0094     }
0095   }
0096 }
0097
0098 void Acts::Extent::set(AxisDirection aDir, double min, double max) {
0099   double minval = min;
0100   if (aDir == AxisDirection::AxisR && minval < 0.) {
0101     minval = 0.;
0102   }
0103   m_range[toUnderlying(aDir)] = Range1D<double>{minval, max};
0104   m_constrains.set(toUnderlying(aDir));
0105 }
0106
0107 void Acts::Extent::setMin(AxisDirection aDir, double min) {
0108   double minval = min;
0109   if (aDir == AxisDirection::AxisR && minval < 0.) {
0110     minval = 0.;
0111   }
0112   m_range[toUnderlying(aDir)].setMin(0u, minval);
0113   m_constrains.set(toUnderlying(aDir));
0114 }
0115
0116 void Acts::Extent::setMax(AxisDirection aDir, double max) {
0117   m_range[toUnderlying(aDir)].setMax(0u, max);
0118   m_constrains.set(toUnderlying(aDir));
0119 }
0120
0121 void Acts::Extent::setEnvelope(const ExtentEnvelope& envelope) {
0122   m_envelope = envelope;
0123 }
0124
0125 bool Acts::Extent::contains(const Vector3& vtx) const {
0126   Extent checkExtent;
0127   for (const auto& bv : allAxisDirections()) {
0128     if (constrains(bv)) {
0129       double vtxVal = VectorHelpers::cast(vtx, bv);
0130       checkExtent.set(bv, vtxVal, vtxVal);
0131     }
0132   }
0133   return contains(checkExtent);
0134 }
0135
0136 bool Acts::Extent::contains(const Extent& rhs,
0137                             std::optional<AxisDirection> aDir) const {
0138   // Helper to check including a constraint bit set check
0139   auto checkContainment = [&](AxisDirection bvc) -> bool {
0140     if (!constrains(bvc)) {
0141       return true;
0142     }
0143     return (rhs.range()[toUnderlying(bvc)] <= m_range[toUnderlying(bvc)]);
0144   };
0145
0146   // Check all
0147   if (!aDir.has_value()) {
0148     return std::ranges::all_of(allAxisDirections(), checkContainment);
0149   }
0150   // Check specific
0151   return checkContainment(aDir.value());
0152 }
0153
0154 bool Acts::Extent::intersects(const Extent& rhs,
0155                               std::optional<AxisDirection> aDir) const {
0156   // Helper to check including a constraint bit set check
0157   auto checkIntersect = [&](AxisDirection bvc) -> bool {
0158     if (!constrains(bvc) || !rhs.constrains(bvc)) {
0159       return false;
0160     }
0161     return (m_range[toUnderlying(bvc)] && rhs.range()[toUnderlying(bvc)]);
0162   };
0163
0164   // Check all
0165   if (!aDir.has_value()) {
0166     return std::ranges::any_of(allAxisDirections(), checkIntersect);
0167   }
0168   // Check specific
0169   return checkIntersect(aDir.value());
0170 }
0171
0172 bool Acts::Extent::constrains(AxisDirection aDir) const {
0173   return m_constrains.test(static_cast<std::size_t>(aDir));
0174 }
0175
0176 bool Acts::Extent::constrains() const {
0177   return m_constrains.count() > 0;
0178 }
0179
0180 bool Acts::Extent::operator==(const Extent& e) const {
0181   if (m_constrains != e.m_constrains) {
0182     return false;
0183   }
0184   if (m_envelope != e.m_envelope) {
0185     return false;
0186   }
0187   if (!(m_range == e.m_range)) {
0188     return false;
0189   }
0190   if (m_valueHistograms != e.m_valueHistograms) {
0191     return false;
0192   }
0193   return true;
0194 }
0195
0196 std::string Acts::Extent::toString(const std::string& indent) const {
0197   std::stringstream sl;
0198   sl << indent << "Extent in space :" << std::endl;
0199   for (const auto& bv : allAxisDirections()) {
0200     if (constrains(bv)) {
0201       sl << indent << "  - value :" << std::setw(10) << axisDirectionName(bv)
0202          << " | range = [" << m_range[toUnderlying(bv)].min() << ", "
0203          << m_range[toUnderlying(bv)].max() << "]" << std::endl;
0204     }
0205   }
0206   return sl.str();
0207 }
0208
0209 // Overload of << operator for std::ostream for debug output
0210 std::ostream& Acts::operator<<(std::ostream& sl, const Extent& rhs) {
0211   sl << rhs.toString();
0212   return sl;
0213 }