EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ConeSurface.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ConeSurface.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2016-2020 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
10 
16 
17 #include <cassert>
18 #include <cmath>
19 
22 
24  : GeometryObject(), Surface(other), m_bounds(other.m_bounds) {}
25 
27  const ConeSurface& other,
28  const Transform3D& shift)
29  : GeometryObject(), Surface(gctx, other, shift), m_bounds(other.m_bounds) {}
30 
32  bool symmetric)
33  : GeometryObject(),
34  Surface(transform),
35  m_bounds(std::make_shared<const ConeBounds>(alpha, symmetric)) {}
36 
38  double zmin, double zmax, double halfPhi)
39  : GeometryObject(),
40  Surface(transform),
41  m_bounds(std::make_shared<const ConeBounds>(alpha, zmin, zmax, halfPhi)) {
42 }
43 
45  const std::shared_ptr<const ConeBounds>& cbounds)
46  : GeometryObject(), Surface(transform), m_bounds(cbounds) {
47  throw_assert(cbounds, "ConeBounds must not be nullptr");
48 }
49 
51  const GeometryContext& gctx, Acts::BinningValue bValue) const {
52  const Vector3D& sfCenter = center(gctx);
53 
54  // special binning type for R-type methods
55  if (bValue == Acts::binR || bValue == Acts::binRPhi) {
56  return Vector3D(sfCenter.x() + bounds().r(sfCenter.z()), sfCenter.y(),
57  sfCenter.z());
58  }
59  // give the center as default for all of these binning types
60  // binX, binY, binZ, binR, binPhi, binRPhi, binH, binEta
61  return sfCenter;
62 }
63 
65  return Surface::Cone;
66 }
67 
69  if (this != &other) {
70  Surface::operator=(other);
71  m_bounds = other.m_bounds;
72  }
73  return *this;
74 }
75 
77  const GeometryContext& gctx) const {
78  return std::move(transform(gctx).matrix().block<3, 1>(0, 2));
79 }
80 
82  const GeometryContext& gctx, const Vector3D& position,
83  const Vector3D& /*unused*/) const {
84  RotationMatrix3D mFrame;
85  // construct the measurement frame
86  // measured Y is the local z axis
87  Vector3D measY = rotSymmetryAxis(gctx);
88  // measured z is the position transverse normalized
89  Vector3D measDepth = Vector3D(position.x(), position.y(), 0.).normalized();
90  // measured X is what comoes out of it
91  Acts::Vector3D measX(measY.cross(measDepth).normalized());
92  // the columnes
93  mFrame.col(0) = measX;
94  mFrame.col(1) = measY;
95  mFrame.col(2) = measDepth;
96  // return the rotation matrix
98  // return it
99  return mFrame;
100 }
101 
103  const GeometryContext& gctx, const Vector2D& lposition,
104  const Vector3D& /*unused*/) const {
105  // create the position in the local 3d frame
106  double r = lposition[Acts::eBoundLoc1] * bounds().tanAlpha();
107  double phi = lposition[Acts::eBoundLoc0] / r;
108  Vector3D loc3Dframe(r * cos(phi), r * sin(phi), lposition[Acts::eBoundLoc1]);
109  return transform(gctx) * loc3Dframe;
110 }
111 
113  const GeometryContext& gctx, const Vector3D& position,
114  const Vector3D& /*unused*/) const {
115  Vector3D loc3Dframe = transform(gctx).inverse() * position;
116  double r = loc3Dframe.z() * bounds().tanAlpha();
117  if (std::abs(perp(loc3Dframe) - r) > s_onSurfaceTolerance) {
118  return Result<Vector2D>::failure(SurfaceError::GlobalPositionNotOnSurface);
119  }
121  {r * atan2(loc3Dframe.y(), loc3Dframe.x()), loc3Dframe.z()});
122 }
123 
125  const Vector3D& position,
126  const Vector3D& direction) const {
127  // (cos phi cos alpha, sin phi cos alpha, sgn z sin alpha)
128  Vector3D posLocal = transform(gctx).inverse() * position;
129  double phi = VectorHelpers::phi(posLocal);
130  double sgn = posLocal.z() > 0. ? -1. : +1.;
131  double cosAlpha = std::cos(bounds().get(ConeBounds::eAlpha));
132  double sinAlpha = std::sin(bounds().get(ConeBounds::eAlpha));
133  Vector3D normalC(cos(phi) * cosAlpha, sin(phi) * cosAlpha, sgn * sinAlpha);
134  normalC = transform(gctx) * normalC;
135  // Back to the global frame
136  double cAlpha = normalC.dot(direction);
137  return std::abs(1. / cAlpha);
138 }
139 
140 std::string Acts::ConeSurface::name() const {
141  return "Acts::ConeSurface";
142 }
143 
145  const GeometryContext& gctx, const Acts::Vector2D& lposition) const {
146  // (cos phi cos alpha, sin phi cos alpha, sgn z sin alpha)
147  double phi = lposition[Acts::eBoundLoc0] /
148  (bounds().r(lposition[Acts::eBoundLoc1])),
149  sgn = lposition[Acts::eBoundLoc1] > 0 ? -1. : +1.;
150  double cosAlpha = std::cos(bounds().get(ConeBounds::eAlpha));
151  double sinAlpha = std::sin(bounds().get(ConeBounds::eAlpha));
152  Vector3D localNormal(cos(phi) * cosAlpha, sin(phi) * cosAlpha,
153  sgn * sinAlpha);
154  return Vector3D(transform(gctx).linear() * localNormal);
155 }
156 
158  const Acts::Vector3D& position) const {
159  // get it into the cylinder frame if needed
160  // @todo respect opening angle
161  Vector3D pos3D = transform(gctx).inverse() * position;
162  pos3D.z() = 0;
163  return pos3D.normalized();
164 }
165 
167  // is safe because no constructor w/o bounds exists
168  return (*m_bounds.get());
169 }
170 
172  const GeometryContext& gctx, size_t lseg) const {
173  // Prepare vertices and faces
174  std::vector<Vector3D> vertices;
175  std::vector<Polyhedron::FaceType> faces;
176  std::vector<Polyhedron::FaceType> triangularMesh;
177 
178  double minZ = bounds().get(ConeBounds::eMinZ);
179  double maxZ = bounds().get(ConeBounds::eMaxZ);
180 
181  if (minZ == -std::numeric_limits<double>::infinity() or
182  maxZ == std::numeric_limits<double>::infinity()) {
183  throw std::domain_error(
184  "Polyhedron repr of boundless surface not possible");
185  }
186 
187  auto ctransform = transform(gctx);
188 
189  // The tip - created only once and only, if the it's not a cut-off cone
190  bool tipExists = false;
191  if (minZ * maxZ <= s_onSurfaceTolerance) {
192  vertices.push_back(ctransform * Vector3D(0., 0., 0.));
193  tipExists = true;
194  }
195 
196  // Cone parameters
197  double hPhiSec = bounds().get(ConeBounds::eHalfPhiSector);
198  double avgPhi = bounds().get(ConeBounds::eAveragePhi);
199  bool fullCone = (hPhiSec == M_PI);
200 
201  // Get the phi segments from the helper
202  auto phiSegs = fullCone ? detail::VerticesHelper::phiSegments()
204  avgPhi - hPhiSec, avgPhi + hPhiSec, {avgPhi});
205 
206  // Negative cone if exists
207  std::vector<double> coneSides;
208  if (std::abs(minZ) > s_onSurfaceTolerance) {
209  coneSides.push_back(minZ);
210  }
211  if (std::abs(maxZ) > s_onSurfaceTolerance) {
212  coneSides.push_back(maxZ);
213  }
214  for (auto& z : coneSides) {
215  // Remember the first vertex
216  size_t firstIv = vertices.size();
217  // Radius and z offset
218  double r = std::abs(z) * bounds().tanAlpha();
219  Vector3D zoffset(0., 0., z);
220  for (unsigned int iseg = 0; iseg < phiSegs.size() - 1; ++iseg) {
221  int addon = (iseg == phiSegs.size() - 2 and not fullCone) ? 1 : 0;
222  detail::VerticesHelper::createSegment(vertices, {r, r}, phiSegs[iseg],
223  phiSegs[iseg + 1], lseg, addon,
224  zoffset, ctransform);
225  }
226  // Create the faces
227  if (tipExists) {
228  for (size_t iv = firstIv + 2; iv < vertices.size() + 1; ++iv) {
229  size_t one = 0, two = iv - 1, three = iv - 2;
230  if (z < 0.) {
231  std::swap(two, three);
232  }
233  faces.push_back({one, two, three});
234  }
235  // Complete cone if necessary
236  if (fullCone) {
237  if (z > 0.) {
238  faces.push_back({0, firstIv, vertices.size() - 1});
239  } else {
240  faces.push_back({0, vertices.size() - 1, firstIv});
241  }
242  }
243  }
244  }
245  // if no tip exists, connect the two bows
246  if (tipExists) {
247  triangularMesh = faces;
248  } else {
249  auto facesMesh =
250  detail::FacesHelper::cylindricalFaceMesh(vertices, fullCone);
251  faces = facesMesh.first;
252  triangularMesh = facesMesh.second;
253  }
254  return Polyhedron(vertices, faces, triangularMesh, false);
255 }