src/Digitization/CartesianSegmentation.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2016-2018 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 http://mozilla.org/MPL/2.0/.
0008
0009 ///////////////////////////////////////////////////////////////////
0010 // CartesianSegmentation.cpp, Acts project
0011 ///////////////////////////////////////////////////////////////////
0012
0013 #include "Acts/Digitization/CartesianSegmentation.hpp"
0014
0015 #include "Acts/Surfaces/PlaneSurface.hpp"
0016 #include "Acts/Surfaces/RectangleBounds.hpp"
0017 #include "Acts/Surfaces/Surface.hpp"
0018 #include "Acts/Utilities/BinningType.hpp"
0019
0020 #include <algorithm>
0021 #include <cmath>
0022 #include <utility>
0023
0024 Acts::CartesianSegmentation::CartesianSegmentation(
0025     const std::shared_ptr<const PlanarBounds>& mBounds, std::size_t numCellsX,
0026     std::size_t numCellsY)
0027     : m_activeBounds(mBounds), m_binUtility(nullptr) {
0028   auto mutableBinUtility = std::make_shared<BinUtility>(
0029       numCellsX, -mBounds->boundingBox().halfLengthX(),
0030       mBounds->boundingBox().halfLengthX(), Acts::open, Acts::binX);
0031   (*mutableBinUtility) +=
0032       BinUtility(numCellsY, -mBounds->boundingBox().halfLengthY(),
0033                  mBounds->boundingBox().halfLengthY(), Acts::open, Acts::binY);
0034   m_binUtility = std::const_pointer_cast<const BinUtility>(mutableBinUtility);
0035 }
0036
0037 Acts::CartesianSegmentation::CartesianSegmentation(
0038     std::shared_ptr<const BinUtility> bUtility,
0039     std::shared_ptr<const PlanarBounds> mBounds)
0040     : m_activeBounds(std::move(mBounds)), m_binUtility(std::move(bUtility)) {
0041   if (!m_activeBounds) {
0042     m_activeBounds = std::make_shared<const RectangleBounds>(
0043         m_binUtility->max(0), m_binUtility->max(1));
0044   }
0045 }
0046
0047 Acts::CartesianSegmentation::~CartesianSegmentation() = default;
0048
0049 void Acts::CartesianSegmentation::createSegmentationSurfaces(
0050     SurfacePtrVector& boundarySurfaces, SurfacePtrVector& segmentationSurfacesX,
0051     SurfacePtrVector& segmentationSurfacesY, double halfThickness,
0052     int readoutDirection, double lorentzAngle) const {
0053   // may be needed throughout
0054   double lorentzAngleTan = tan(lorentzAngle);
0055   double lorentzPlaneShiftX = halfThickness * lorentzAngleTan;
0056
0057   // (A) --- top/bottom surfaces
0058   // -----------------------------------------------------------
0059   // let's create the top/botten surfaces first - we call them readout / counter
0060   // readout
0061   // there are some things to consider
0062   // - they share the RectangleBounds only if the lorentzAngle is 0
0063   // otherwise only the readout surface has full length bounds like the module
0064   std::shared_ptr<const PlanarBounds> moduleBounds(
0065       new RectangleBounds(m_activeBounds->boundingBox()));
0066   // - they are separated by half a thickness in z
0067   auto readoutPlaneTransform = Transform3::Identity();
0068   auto counterPlaneTransform = Transform3::Identity();
0069   // readout and counter readout bounds, the bounds of the readout plane are
0070   // like the active ones
0071   std::shared_ptr<const PlanarBounds> readoutPlaneBounds = moduleBounds;
0072   std::shared_ptr<const PlanarBounds> counterPlaneBounds(nullptr);
0073   // the transform of the readout plane is always centric
0074   readoutPlaneTransform.translation() =
0075       Vector3(0., 0., readoutDirection * halfThickness);
0076   // no lorentz angle and everything is straight-forward
0077   if (lorentzAngle == 0.) {
0078     counterPlaneBounds = moduleBounds;
0079     counterPlaneTransform.translation() =
0080         Vector3(0., 0., -readoutDirection * halfThickness);
0081   } else {
0082     // lorentz reduced Bounds
0083     double lorentzReducedHalfX =
0084         m_activeBounds->boundingBox().halfLengthX() - fabs(lorentzPlaneShiftX);
0085     std::shared_ptr<const PlanarBounds> lorentzReducedBounds(
0086         new RectangleBounds(lorentzReducedHalfX,
0087                             m_activeBounds->boundingBox().halfLengthY()));
0088     counterPlaneBounds = lorentzReducedBounds;
0089     // now we shift the counter plane in position - this depends on lorentz
0090     // angle
0091     double counterPlaneShift = -readoutDirection * lorentzPlaneShiftX;
0092     counterPlaneTransform.translation() =
0093         Vector3(counterPlaneShift, 0., -readoutDirection * halfThickness);
0094   }
0095   // - build the readout & counter readout surfaces
0096   boundarySurfaces.push_back(Surface::makeShared<PlaneSurface>(
0097       readoutPlaneTransform, readoutPlaneBounds));
0098   boundarySurfaces.push_back(Surface::makeShared<PlaneSurface>(
0099       counterPlaneTransform, counterPlaneBounds));
0100
0101   // (B) - bin X and lorentz surfaces
0102   // -----------------------------------------------------------
0103   // easy stuff first, constant pitch size and
0104   double pitchX =
0105       2. * m_activeBounds->boundingBox().halfLengthX() / m_binUtility->bins(0);
0106
0107   // now, let's create the shared bounds of all surfaces marking x bins - choice
0108   // fixes orientation of the matrix
0109   std::shared_ptr<const PlanarBounds> xBinBounds(new RectangleBounds(
0110       m_activeBounds->boundingBox().halfLengthY(), halfThickness));
0111   // now, let's create the shared bounds of all surfaces marking lorentz planes
0112   double lorentzPlaneHalfX = std::abs(halfThickness / cos(lorentzAngle));
0113   // the bounds of the lorentz plane
0114   std::shared_ptr<const PlanarBounds> lorentzPlaneBounds =
0115       (lorentzAngle == 0.)
0116           ? xBinBounds
0117           : std::shared_ptr<const PlanarBounds>(
0118                 new RectangleBounds(m_activeBounds->boundingBox().halfLengthY(),
0119                                     lorentzPlaneHalfX));
0120
0121   // now the rotation matrix for the xBins
0122   RotationMatrix3 xBinRotationMatrix;
0123   xBinRotationMatrix.col(0) = Vector3::UnitY();
0124   xBinRotationMatrix.col(1) = Vector3::UnitZ();
0125   xBinRotationMatrix.col(2) = Vector3::UnitX();
0126   // now the lorentz plane rotation should be the xBin rotation, rotated by the
0127   // lorentz angle around y
0128   RotationMatrix3 lorentzPlaneRotationMatrix =
0129       (lorentzAngle != 0.)
0130           ? xBinRotationMatrix * AngleAxis3(lorentzAngle, Vector3::UnitX())
0131           : xBinRotationMatrix;
0132
0133   // reserve, it's always (number of bins-1) as the boundaries are within the
0134   // boundarySurfaces
0135   segmentationSurfacesX.reserve(m_binUtility->bins(0));
0136   // create and fill them
0137   for (std::size_t ibinx = 0; ibinx <= m_binUtility->bins(0); ++ibinx) {
0138     // the current step x position
0139     double cPosX =
0140         -m_activeBounds->boundingBox().halfLengthX() + ibinx * pitchX;
0141     // (i) this is the low/high boundary --- ( ibin == 0/m_binUtility->bins(0) )
0142     if ((ibinx == 0u) || ibinx == m_binUtility->bins(0)) {
0143       // check if it is a straight boundary or not: always straight for no
0144       // lorentz angle, and either the first boundary or the last depending on
0145       // lorentz and readout
0146       bool boundaryStraight =
0147           (lorentzAngle == 0. ||
0148            ((ibinx == 0u) && readoutDirection * lorentzAngle > 0.) ||
0149            (ibinx == m_binUtility->bins(0) &&
0150             readoutDirection * lorentzAngle < 0));
0151       // set the low boundary parameters : position & rotation
0152       Vector3 boundaryXPosition =
0153           boundaryStraight
0154               ? Vector3(cPosX, 0., 0.)
0155               : Vector3(cPosX - readoutDirection * lorentzPlaneShiftX, 0., 0.);
0156       // rotation of the boundary: straight or lorentz
0157       const RotationMatrix3& boundaryXRotation =
0158           boundaryStraight ? xBinRotationMatrix : lorentzPlaneRotationMatrix;
0159       // build the rotation from it
0160       auto boundaryXTransform =
0161           Transform3(Translation3(boundaryXPosition) * boundaryXRotation);
0162       // the correct bounds for this
0163       std::shared_ptr<const PlanarBounds> boundaryXBounds =
0164           boundaryStraight ? xBinBounds : lorentzPlaneBounds;
0165       // boundary surfaces
0166       boundarySurfaces.push_back(Surface::makeShared<PlaneSurface>(
0167           boundaryXTransform, boundaryXBounds));
0168       // (ii) this is the in between bins  --- ( 1 <= ibin < m_mbnsX )
0169     } else {
0170       // shift by the lorentz angle
0171       Vector3 lorentzPlanePosition(
0172           cPosX - readoutDirection * lorentzPlaneShiftX, 0., 0.);
0173       auto lorentzPlaneTransform = Transform3(
0174           Translation3(lorentzPlanePosition) * lorentzPlaneRotationMatrix);
0175       // lorentz plane surfaces
0176       segmentationSurfacesX.push_back(Surface::makeShared<PlaneSurface>(
0177           lorentzPlaneTransform, lorentzPlaneBounds));
0178     }
0179   }
0180
0181   // (C) - bin Y surfaces - everything is defined
0182   // -----------------------------------------------------------
0183   // now the rotation matrix for the yBins - anticyclic
0184   RotationMatrix3 yBinRotationMatrix;
0185   yBinRotationMatrix.col(0) = Vector3::UnitX();
0186   yBinRotationMatrix.col(1) = Vector3::UnitZ();
0187   yBinRotationMatrix.col(2) = Vector3(0., -1., 0.);
0188   // easy stuff first, constant pitch in Y
0189   double pitchY =
0190       2. * m_activeBounds->boundingBox().halfLengthY() / m_binUtility->bins(1);
0191   // let's create the shared bounds of all surfaces marking y bins
0192   std::shared_ptr<const PlanarBounds> yBinBounds(new RectangleBounds(
0193       m_activeBounds->boundingBox().halfLengthX(), halfThickness));
0194   // reserve, it's always (number of bins-1) as the boundaries are within the
0195   // boundarySurfaces
0196   segmentationSurfacesY.reserve(m_binUtility->bins(1));
0197   for (std::size_t ibiny = 0; ibiny <= m_binUtility->bins(1); ++ibiny) {
0198     // the position of the bin surface
0199     double binPosY =
0200         -m_activeBounds->boundingBox().halfLengthY() + ibiny * pitchY;
0201     Vector3 binSurfaceCenter(0., binPosY, 0.);
0202     // the binning transform
0203     auto binTransform =
0204         Transform3(Translation3(binSurfaceCenter) * yBinRotationMatrix);
0205     // these are the boundaries
0206     if (ibiny == 0 || ibiny == m_binUtility->bins(1)) {
0207       boundarySurfaces.push_back(
0208           Surface::makeShared<PlaneSurface>(binTransform, yBinBounds));
0209     } else {  // these are the bin boundaries
0210       segmentationSurfacesY.push_back(
0211           Surface::makeShared<PlaneSurface>(binTransform, yBinBounds));
0212     }
0213   }
0214 }
0215
0216 Acts::Vector2 Acts::CartesianSegmentation::cellPosition(
0217     const DigitizationCell& dCell) const {
0218   double bX = m_binUtility->bins(0) > 1
0219                   ? m_binUtility->binningData()[0].center(dCell.channel0)
0220                   : 0.;
0221   double bY = m_binUtility->bins(1) > 1
0222                   ? m_binUtility->binningData()[1].center(dCell.channel1)
0223                   : 0.;
0224   return Vector2(bX, bY);
0225 }
0226
0227 /** Get the digitization cell from 3D position, it used the projection to the
0228  * readout surface to estimate the 2D position */
0229 Acts::DigitizationStep Acts::CartesianSegmentation::digitizationStep(
0230     const Vector3& startStep, const Vector3& endStep, double halfThickness,
0231     int readoutDirection, double lorentzAngle) const {
0232   Vector3 stepCenter = 0.5 * (startStep + endStep);
0233   // take the full drift length
0234   // this is the absolute drift in z
0235   double driftInZ = halfThickness - readoutDirection * stepCenter.z();
0236   // this is the absolute drift length
0237   double driftLength = driftInZ / cos(lorentzAngle);
0238   // project to parameter the readout surface
0239   double lorentzDeltaX = readoutDirection * driftInZ * tan(lorentzAngle);
0240   // the projected center, it has the lorentz shift applied
0241   Vector2 stepCenterProjected(stepCenter.x() + lorentzDeltaX, stepCenter.y());
0242   // the cell & its center
0243   Acts::DigitizationCell dCell = cell(stepCenterProjected);
0244   Vector2 cellCenter = cellPosition(dCell);
0245   // we are ready to return what we have
0246   return DigitizationStep((endStep - startStep).norm(), driftLength, dCell,
0247                           startStep, endStep, stepCenterProjected, cellCenter);
0248 }