VerticesHelper.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file VerticesHelper.hpp

// This file is part of the Acts project.
//
// Copyright (C) 2020 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 http://mozilla.org/MPL/2.0/.

#pragma once

#include "Acts/Definitions/Algebra.hpp"
#include "Acts/Definitions/Tolerance.hpp"

#include <algorithm>
#include <cmath>
#include <utility>
#include <vector>

namespace Acts {
namespace detail {
namespace VerticesHelper {

std::vector<ActsScalar> phiSegments(ActsScalar phiMin = -M_PI,
                                    ActsScalar phiMax = M_PI,
                                    const std::vector<ActsScalar>& phiRefs = {},
                                    ActsScalar phiTolerance = 1e-6);

template <typename vertex_t, typename transform_t>
void createSegment(std::vector<vertex_t>& vertices,
                   std::pair<ActsScalar, ActsScalar> rxy, ActsScalar phi1,
                   ActsScalar phi2, unsigned int lseg, int addon = 0,
                   const vertex_t& offset = vertex_t::Zero(),
                   const transform_t& transform = transform_t::Identity()) {
  // Calculate the number of segments - 1 is the minimum
  unsigned int segs =
      static_cast<unsigned int>(std::abs(phi2 - phi1) / (2 * M_PI) * lseg);
  segs = segs > 0 ? segs : 1;
  ActsScalar phistep = (phi2 - phi1) / segs;
  // Create the segments
  for (unsigned int iphi = 0; iphi < segs + addon; ++iphi) {
    ActsScalar phi = phi1 + iphi * phistep;
    vertex_t vertex = vertex_t::Zero();
    vertex(0) = rxy.first * std::cos(phi);
    vertex(1) = rxy.second * std::sin(phi);

    vertex = vertex + offset;
    vertices.push_back(transform * vertex);
  }
}

std::vector<Vector2> ellipsoidVertices(ActsScalar innerRx, ActsScalar innerRy,
                                       ActsScalar outerRx, ActsScalar outerRy,
                                       ActsScalar avgPhi = 0.,
                                       ActsScalar halfPhi = M_PI,
                                       unsigned int lseg = 1);

std::vector<Vector2> circularVertices(ActsScalar innerR, ActsScalar outerR,
                                      ActsScalar avgPhi = 0.,
                                      ActsScalar halfPhi = M_PI,
                                      unsigned int lseg = 1);
template <typename vertex_t, typename vertex_container_t>
bool isInsidePolygon(const vertex_t& point,
                     const vertex_container_t& vertices) {
  // when we move along the edges of a convex polygon, a point on the inside of
  // the polygon will always appear on the same side of each edge.
  // a point on the outside will switch sides at least once.

  // returns which side of the connecting line between `ll0` and `ll1` the point
  // `p` is on. computes the sign of the z-component of the cross-product
  // between the line normal vector and the vector from `ll0` to `p`.
  auto lineSide = [&](auto&& ll0, auto&& ll1) {
    auto normal = ll1 - ll0;
    auto delta = point - ll0;
    return std::signbit((normal[0] * delta[1]) - (normal[1] * delta[0]));
  };

  auto iv = std::begin(vertices);
  auto l0 = *iv;
  auto l1 = *(++iv);
  // use vertex0 to vertex1 to define reference sign and compare w/ all edges
  auto reference = lineSide(l0, l1);
  for (++iv; iv != std::end(vertices); ++iv) {
    l0 = l1;
    l1 = *iv;
    if (lineSide(l0, l1) != reference) {
      return false;
    }
  }
  // manual check for last edge from last vertex back to the first vertex
  if (lineSide(l1, *std::begin(vertices)) != reference) {
    return false;
  }
  // point was always on the same side. point must be inside.
  return true;
}

template <typename vertex_t>
bool isInsideRectangle(const vertex_t& point, const vertex_t& lowerLeft,
                       const vertex_t& upperRight) {
  return (lowerLeft[0] <= point[0]) && (point[0] < upperRight[0]) &&
         (lowerLeft[1] <= point[1]) && (point[1] < upperRight[1]);
}

bool onHyperPlane(const std::vector<Vector3>& vertices,
                  ActsScalar tolerance = s_onSurfaceTolerance);

}  // namespace VerticesHelper
}  // namespace detail
}  // namespace Acts