src/Digitization/SurfaceDrift.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 "ActsFatras/Digitization/SurfaceDrift.hpp"
0010
0011 #include "Acts/Definitions/Tolerance.hpp"
0012 #include "Acts/Utilities/Helpers.hpp"
0013 #include "ActsFatras/Digitization/DigitizationError.hpp"
0014
0015 #include <cmath>
0016
0017 namespace ActsFatras {
0018
0019 namespace {
0020
0021 /// Express the hit position and direction in the readout-local 3D frame.
0022 ///
0023 /// The returned position has its third component (the surface-normal
0024 /// coordinate) set to 0 — the readout reference plane — while the direction is
0025 /// a unit vector whose third component is the projection onto the surface
0026 /// normal. The first two components are the in-plane readout coordinates:
0027 ///   - Plane / Disc : Cartesian local (x, y)
0028 ///   - Cylinder     : unrolled (rPhi, z)
0029 ///
0030 /// Returns false if the local frame cannot be established.
0031 bool toLocalFrame(const Acts::GeometryContext& gctx,
0032                   const Acts::Surface& surface, const Acts::Vector3& pos,
0033                   const Acts::Vector3& dir, Acts::Vector3& pos3Local,
0034                   Acts::Vector3& dir3Local) {
0035   const auto invTransform = surface.localToGlobalTransform(gctx).inverse();
0036
0037   if (surface.type() == Acts::Surface::SurfaceType::Cylinder) {
0038     // Position/direction in the cylinder-axis frame (axis aligned with local z)
0039     const Acts::Vector3 posCyl = invTransform * pos;
0040     const Acts::Vector3 dirCyl = invTransform.linear() * dir.normalized();
0041
0042     const double phi = std::atan2(posCyl.y(), posCyl.x());
0043     const double R = std::hypot(posCyl.x(), posCyl.y());
0044     const double cphi = std::cos(phi);
0045     const double sphi = std::sin(phi);
0046
0047     // Orthonormal readout-local basis at the hit:
0048     //   tangential = (-sinφ, cosφ, 0), axial = (0, 0, 1), radial = (cosφ, sinφ,
0049     //   0)
0050     const double dsRPhi = -sphi * dirCyl.x() + cphi * dirCyl.y();
0051     const double dsZ = dirCyl.z();
0052     const double dsR = cphi * dirCyl.x() + sphi * dirCyl.y();
0053
0054     // (rPhi, z) readout position, normal coordinate = 0
0055     pos3Local = Acts::Vector3(R * phi, posCyl.z(), 0.);
0056     // (tangential, axial, radial); radial is the surface-normal component
0057     dir3Local = Acts::Vector3(dsRPhi, dsZ, dsR);
0058     return true;
0059   }
0060
0061   // Plane / Disc: the Cartesian local frame, surface normal = local z.
0062   pos3Local = invTransform * pos;
0063   dir3Local = invTransform.linear() * dir.normalized();
0064   return true;
0065 }
0066
0067 }  // namespace
0068
0069 Acts::Result<std::tuple<SurfaceDrift::Segment2D, SurfaceDrift::Segment3D>>
0070 SurfaceDrift::toReadout(const Acts::GeometryContext& gctx,
0071                         const Acts::Surface& surface, double thickness,
0072                         const Acts::Vector3& pos, const Acts::Vector3& dir,
0073                         const Acts::Vector3& driftDir) const {
0074   Acts::Vector3 pos3Local = Acts::Vector3::Zero();
0075   Acts::Vector3 seg3Local = Acts::Vector3::Zero();
0076   if (!toLocalFrame(gctx, surface, pos, dir, pos3Local, seg3Local)) {
0077     return DigitizationError::DriftError;
0078   }
0079
0080   // seg3Local is a unit direction; its z-component is the projection onto the
0081   // surface normal. If it (nearly) vanishes the track is parallel to the
0082   // surface and the drift is undefined.
0083   if (std::abs(seg3Local.z()) < Acts::s_epsilon) {
0084     return DigitizationError::DriftError;
0085   }
0086
0087   // Scale the unit vector to the thickness of the module
0088   const double scale = thickness / seg3Local.z();
0089   seg3Local *= scale;
0090   // The drift direction is in the local frame, so we need to transform it
0091   const Acts::Vector3 entry = pos3Local - 0.5 * seg3Local;
0092   const Acts::Vector3 exit = pos3Local + 0.5 * seg3Local;
0093   Acts::Vector3 driftedEntry = entry;
0094   Acts::Vector3 driftedExit = exit;
0095   // Apply a Lorentz drift if the drift direction has both an in-plane and a
0096   // perpendicular (normal) component. driftDir is expressed in the same
0097   // readout-local frame as seg3Local.
0098   if (Acts::VectorHelpers::perp(driftDir) > Acts::s_epsilon &&
0099       std::abs(driftDir.z()) > Acts::s_epsilon) {
0100     // Apply the drift to the entry and exit points
0101     const double driftScale = 0.5 * thickness / driftDir.z();
0102     driftedEntry += driftScale * driftDir;
0103     driftedExit -= driftScale * driftDir;
0104   }
0105
0106   // The drifted 2D readout segment and the undrifted 3D segment (whose length
0107   // is the true chord through the depletion zone, used downstream to rescale
0108   // the 2D activation path to 3D).
0109   return std::tuple{
0110       Segment2D{driftedEntry.segment<2>(0), driftedExit.segment<2>(0)},
0111       Segment3D{entry, exit}};
0112 }
0113
0114 }  // namespace ActsFatras