cylinder3D.hpp

// This file is part of the ACTS project.
//
// Copyright (C) 2016 CERN for the benefit of the ACTS project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
#pragma once
 
// Project include(s)
#include "detray/definitions/algebra.hpp"
#include "detray/definitions/containers.hpp"
#include "detray/definitions/detail/qualifiers.hpp"
#include "detray/definitions/indexing.hpp"
#include "detray/definitions/math.hpp"
#include "detray/definitions/units.hpp"
#include "detray/geometry/coordinates/cylindrical3D.hpp"
#include "detray/geometry/detail/shape_utils.hpp"
 
// System include(s)
#include <limits>
#include <ostream>
#include <stdexcept>
#include <string_view>
 
namespace detray {

class cylinder3D {
 public:
  static constexpr std::string_view name = "cylinder3D";

  enum boundaries : unsigned int {
    e_min_r = 0u,
    e_min_phi = 1u,
    e_min_z = 2u,
    e_max_r = 3u,
    e_max_phi = 4u,
    e_max_z = 5u,
    e_size = 6u,
  };

  template <concepts::scalar scalar_t>
  using bounds_type = darray<scalar_t, boundaries::e_size>;

  template <concepts::algebra algebra_t>
  using local_frame_type = cylindrical3D<algebra_t>;

  template <typename bool_t>
  using result_type = bool_t;

  static constexpr std::size_t dim{3u};

  template <concepts::scalar scalar_t, concepts::point point_t>
  DETRAY_HOST_DEVICE inline scalar_t min_dist_to_boundary(
      const bounds_type<scalar_t> &bounds, const point_t &loc_p) const {
    const scalar_t min_r_dist =
        math::min(math::fabs(loc_p[0] - bounds[e_min_r]),
                  math::fabs(bounds[e_max_r] - loc_p[0]));
    const scalar_t min_phi_dist =
        math::min(math::fabs(loc_p[1] - bounds[e_min_phi]),
                  math::fabs(bounds[e_max_phi] - loc_p[1]));
    const scalar_t min_z_dist =
        math::min(math::fabs(loc_p[2] - bounds[e_min_z]),
                  math::fabs(bounds[e_max_z] - loc_p[2]));
 
    // Use the chord for the phi distance
    return math::min(
        math::min(min_r_dist, 2.f * loc_p[0] * math::sin(0.5f * min_phi_dist)),
        min_z_dist);
  }

  template <concepts::algebra algebra_t>
  DETRAY_HOST_DEVICE constexpr result_type<dbool<algebra_t>> check_boundaries(
      const bounds_type<dscalar<algebra_t>> &bounds,
      const dtransform3D<algebra_t> &trf, const dpoint3D<algebra_t> &glob_p,
      const dscalar<algebra_t> tol =
          std::numeric_limits<dscalar<algebra_t>>::epsilon(),
      const dscalar<algebra_t> /*edge_tol*/ = 0.f) const {
    // Get the full local position
    const auto loc_p =
        local_frame_type<algebra_t>::global_to_local(trf, glob_p, {});
 
    return check_boundaries(bounds, loc_p, tol);
  }

  template <concepts::scalar scalar_t, concepts::point point_t>
  DETRAY_HOST_DEVICE constexpr auto check_boundaries(
      const bounds_type<scalar_t> &bounds, const point_t &loc_p,
      const scalar_t tol = std::numeric_limits<scalar_t>::epsilon(),
      const scalar_t /*edge_tol*/ = 0.f) const {
    const scalar_t phi_tol = detail::phi_tolerance(tol, loc_p[0]);
 
    return ((bounds[e_min_r] - tol) <= loc_p[0] &&
            (bounds[e_min_phi] - phi_tol) <= loc_p[1] &&
            (bounds[e_min_z] - tol) <= loc_p[2] &&
            loc_p[0] <= (bounds[e_max_r] + tol) &&
            loc_p[1] <= (bounds[e_max_phi] + phi_tol) &&
            loc_p[2] <= (bounds[e_max_z] + tol));
  }


  template <concepts::scalar scalar_t>
  DETRAY_HOST_DEVICE constexpr scalar_t measure(
      const bounds_type<scalar_t> &bounds) const {
    return volume(bounds);
  }

  template <concepts::scalar scalar_t>
  DETRAY_HOST_DEVICE constexpr scalar_t volume(
      const bounds_type<scalar_t> &bounds) const {
    return constant<scalar_t>::pi * (bounds[e_max_z] - bounds[e_min_z]) *
           (bounds[e_max_r] * bounds[e_max_r] -
            bounds[e_min_r] * bounds[e_min_r]);
  }

  template <concepts::scalar scalar_t>
  DETRAY_HOST_DEVICE constexpr bounds_type<scalar_t> merge(
      const bounds_type<scalar_t> &bounds,
      const bounds_type<scalar_t> &o_bounds) const {
    bounds_type<scalar_t> new_bounds{};
 
    new_bounds[e_min_r] = math::min(bounds[e_min_r], o_bounds[e_min_r]);
    new_bounds[e_min_phi] = math::min(bounds[e_min_phi], o_bounds[e_min_phi]);
    new_bounds[e_min_z] = math::min(bounds[e_min_z], o_bounds[e_min_z]);
    new_bounds[e_max_r] = math::max(bounds[e_max_r], o_bounds[e_max_r]);
    new_bounds[e_max_phi] = math::max(bounds[e_max_phi], o_bounds[e_max_phi]);
    new_bounds[e_max_z] = math::max(bounds[e_max_z], o_bounds[e_max_z]);
 
    return new_bounds;
  }

  // @todo: Look at phi - range for a better fit
  template <concepts::algebra algebra_t>
  DETRAY_HOST_DEVICE inline darray<dscalar<algebra_t>, 6> local_min_bounds(
      const bounds_type<dscalar<algebra_t>> &bounds,
      const dscalar<algebra_t> env =
          std::numeric_limits<dscalar<algebra_t>>::epsilon()) const {
    assert(env > 0.f);
    const dscalar<algebra_t> r_bound{bounds[e_max_r] + env};
    return {-r_bound, -r_bound, bounds[e_min_z] - env,
            r_bound,  r_bound,  bounds[e_max_z] + env};
  }

  template <concepts::algebra algebra_t>
  DETRAY_HOST_DEVICE dpoint3D<algebra_t> centroid(
      const bounds_type<dscalar<algebra_t>> &bounds) const {
    return 0.5f * dpoint3D<algebra_t>{static_cast<dscalar<algebra_t>>(0.f),
                                      (bounds[e_min_phi] + bounds[e_max_phi]),
                                      (bounds[e_min_z] + bounds[e_max_z])};
  }

  template <concepts::algebra algebra_t>
  DETRAY_HOST dvector<dpoint3D<algebra_t>> vertices(
      const bounds_type<dscalar<algebra_t>> & /*unused*/,
      dindex /*unused*/) const {
    throw std::runtime_error(
        "Vertex generation for 3D cylinders is not implemented");
    return {};
  }

  template <concepts::scalar scalar_t>
  DETRAY_HOST constexpr bool check_consistency(
      const bounds_type<scalar_t> &bounds, std::ostream &os) const {
    constexpr auto tol{10.f * std::numeric_limits<scalar_t>::epsilon()};
 
    if (bounds[e_min_r] < tol) {
      os << "DETRAY ERROR (HOST): Radii must be in the range (0, "
            "numeric_max)"
         << std::endl;
      return false;
    }
    if (bounds[e_min_r] >= bounds[e_max_r] ||
        math::fabs(bounds[e_min_r] - bounds[e_max_r]) < tol) {
      os << "DETRAY ERROR (HOST): Min Radius must be smaller than max "
            "Radius.";
      return false;
    }
    if (bounds[e_min_z] >= bounds[e_max_z] ||
        math::fabs(bounds[e_min_z] - bounds[e_max_z]) < tol) {
      os << "DETRAY ERROR (HOST): Min z must be smaller than max z.";
      return false;
    }
 
    return true;
  }
};
 
}  // namespace detray