ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CartesianSegmentation.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CartesianSegmentation.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2016-2018 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 // CartesianSegmentation.cpp, Acts project
12 
13 #include <utility>
14 
19 
21  const std::shared_ptr<const PlanarBounds>& mBounds, size_t numCellsX,
22  size_t numCellsY)
23  : m_activeBounds(mBounds), m_binUtility(nullptr) {
24  auto mutableBinUtility = std::make_shared<BinUtility>(
25  numCellsX, -mBounds->boundingBox().halfLengthX(),
26  mBounds->boundingBox().halfLengthX(), Acts::open, Acts::binX);
27  (*mutableBinUtility) +=
28  BinUtility(numCellsY, -mBounds->boundingBox().halfLengthY(),
29  mBounds->boundingBox().halfLengthY(), Acts::open, Acts::binY);
30  m_binUtility = std::const_pointer_cast<const BinUtility>(mutableBinUtility);
31 }
32 
34  std::shared_ptr<const BinUtility> bUtility,
35  std::shared_ptr<const PlanarBounds> mBounds)
36  : m_activeBounds(std::move(mBounds)), m_binUtility(std::move(bUtility)) {
37  if (!m_activeBounds) {
38  m_activeBounds = std::make_shared<const RectangleBounds>(
39  m_binUtility->max(0), m_binUtility->max(1));
40  }
41 }
42 
44 
46  SurfacePtrVector& boundarySurfaces, SurfacePtrVector& segmentationSurfacesX,
47  SurfacePtrVector& segmentationSurfacesY, double halfThickness,
48  int readoutDirection, double lorentzAngle) const {
49  // may be needed throughout
50  double lorentzAngleTan = tan(lorentzAngle);
51  double lorentzPlaneShiftX = halfThickness * lorentzAngleTan;
52 
53  // (A) --- top/bottom surfaces
54  // -----------------------------------------------------------
55  // let's create the top/botten surfaces first - we call them readout / counter
56  // readout
57  // there are some things to consider
58  // - they share the RectangleBounds only if the lorentzAngle is 0
59  // otherwise only the readout surface has full length bounds like the module
60  std::shared_ptr<const PlanarBounds> moduleBounds(
61  new RectangleBounds(m_activeBounds->boundingBox()));
62  // - they are separated by half a thickness in z
63  auto mutableReadoutPlaneTransform =
64  std::make_shared<Transform3D>(Transform3D::Identity());
65  auto mutableCounterPlaneTransform =
66  std::make_shared<Transform3D>(Transform3D::Identity());
67  // readout and counter readout bounds, the bounds of the readout plane are
68  // like the active ones
69  std::shared_ptr<const PlanarBounds> readoutPlaneBounds = moduleBounds;
70  std::shared_ptr<const PlanarBounds> counterPlaneBounds(nullptr);
71  // the transform of the readout plane is always centric
72  (*mutableReadoutPlaneTransform).translation() =
73  Vector3D(0., 0., readoutDirection * halfThickness);
74  // no lorentz angle and everything is straight-forward
75  if (lorentzAngle == 0.) {
76  counterPlaneBounds = moduleBounds;
77  (*mutableCounterPlaneTransform).translation() =
78  Vector3D(0., 0., -readoutDirection * halfThickness);
79  } else {
80  // lorentz reduced Bounds
81  double lorentzReducedHalfX =
82  m_activeBounds->boundingBox().halfLengthX() - fabs(lorentzPlaneShiftX);
83  std::shared_ptr<const PlanarBounds> lorentzReducedBounds(
84  new RectangleBounds(lorentzReducedHalfX,
85  m_activeBounds->boundingBox().halfLengthY()));
86  counterPlaneBounds = lorentzReducedBounds;
87  // now we shift the counter plane in position - this depends on lorentz
88  // angle
89  double counterPlaneShift = -readoutDirection * lorentzPlaneShiftX;
90  (*mutableCounterPlaneTransform).translation() =
91  Vector3D(counterPlaneShift, 0., -readoutDirection * halfThickness);
92  }
93  // - finalize the transforms
94  auto readoutPlaneTransform =
95  std::const_pointer_cast<const Transform3D>(mutableReadoutPlaneTransform);
96  auto counterPlaneTransform =
97  std::const_pointer_cast<const Transform3D>(mutableCounterPlaneTransform);
98  // - build the readout & counter readout surfaces
99  boundarySurfaces.push_back(Surface::makeShared<PlaneSurface>(
100  readoutPlaneTransform, readoutPlaneBounds));
101  boundarySurfaces.push_back(Surface::makeShared<PlaneSurface>(
102  counterPlaneTransform, counterPlaneBounds));
103 
104  // (B) - bin X and lorentz surfaces
105  // -----------------------------------------------------------
106  // easy stuff first, constant pitch size and
107  double pitchX =
108  2. * m_activeBounds->boundingBox().halfLengthX() / m_binUtility->bins(0);
109 
110  // now, let's create the shared bounds of all surfaces marking x bins - choice
111  // fixes orientation of the matrix
112  std::shared_ptr<const PlanarBounds> xBinBounds(new RectangleBounds(
113  m_activeBounds->boundingBox().halfLengthY(), halfThickness));
114  // now, let's create the shared bounds of all surfaces marking lorentz planes
115  double lorentzPlaneHalfX = std::abs(halfThickness / cos(lorentzAngle));
116  // the bounds of the lorentz plane
117  std::shared_ptr<const PlanarBounds> lorentzPlaneBounds =
118  (lorentzAngle == 0.)
119  ? xBinBounds
120  : std::shared_ptr<const PlanarBounds>(
121  new RectangleBounds(m_activeBounds->boundingBox().halfLengthY(),
122  lorentzPlaneHalfX));
123 
124  // now the rotation matrix for the xBins
125  RotationMatrix3D xBinRotationMatrix;
126  xBinRotationMatrix.col(0) = Vector3D::UnitY();
127  xBinRotationMatrix.col(1) = Vector3D::UnitZ();
128  xBinRotationMatrix.col(2) = Vector3D::UnitX();
129  // now the lorentz plane rotation should be the xBin rotation, rotated by the
130  // lorentz angle around y
131  RotationMatrix3D lorentzPlaneRotationMatrix =
132  (lorentzAngle != 0.)
133  ? xBinRotationMatrix * AngleAxis3D(lorentzAngle, Vector3D::UnitX())
134  : xBinRotationMatrix;
135 
136  // reserve, it's always (number of bins-1) as the boundaries are within the
137  // boundarySurfaces
138  segmentationSurfacesX.reserve(m_binUtility->bins(0));
139  // create and fill them
140  for (size_t ibinx = 0; ibinx <= m_binUtility->bins(0); ++ibinx) {
141  // the current step x position
142  double cPosX =
143  -m_activeBounds->boundingBox().halfLengthX() + ibinx * pitchX;
144  // (i) this is the low/high boundary --- ( ibin == 0/m_binUtility->bins(0) )
145  if ((ibinx == 0u) || ibinx == m_binUtility->bins(0)) {
146  // check if it a straight boundary or not: always straight for no lorentz
147  // angle,
148  // and either the first boundary or the last dependening on lorentz &
149  // readout
150  bool boundaryStraight =
151  (lorentzAngle == 0. ||
152  ((ibinx == 0u) && readoutDirection * lorentzAngle > 0.) ||
153  (ibinx == m_binUtility->bins(0) &&
154  readoutDirection * lorentzAngle < 0));
155  // set the low boundary parameters : position & rotation
156  Vector3D boundaryXPosition =
157  boundaryStraight
158  ? Vector3D(cPosX, 0., 0.)
159  : Vector3D(cPosX - readoutDirection * lorentzPlaneShiftX, 0., 0.);
160  // rotation of the boundary: striaght or lorentz
161  const RotationMatrix3D& boundaryXRotation =
162  boundaryStraight ? xBinRotationMatrix : lorentzPlaneRotationMatrix;
163  // build the rotation from it
164  auto boundaryXTransform = std::make_shared<const Transform3D>(
165  Translation3D(boundaryXPosition) * boundaryXRotation);
166  // the correct bounds for this
167  std::shared_ptr<const PlanarBounds> boundaryXBounds =
168  boundaryStraight ? xBinBounds : lorentzPlaneBounds;
169  // boundary surfaces
170  boundarySurfaces.push_back(Surface::makeShared<PlaneSurface>(
171  boundaryXTransform, boundaryXBounds));
172  // (ii) this is the in between bins --- ( 1 <= ibin < m_mbnsX )
173  } else {
174  // shift by the lorentz angle
175  Vector3D lorentzPlanePosition(
176  cPosX - readoutDirection * lorentzPlaneShiftX, 0., 0.);
177  auto lorentzPlaneTransform = std::make_shared<const Transform3D>(
178  Translation3D(lorentzPlanePosition) * lorentzPlaneRotationMatrix);
179  // lorentz plane surfaces
180  segmentationSurfacesX.push_back(Surface::makeShared<PlaneSurface>(
181  lorentzPlaneTransform, lorentzPlaneBounds));
182  }
183  }
184 
185  // (C) - bin Y surfaces - everything is defined
186  // -----------------------------------------------------------
187  // now the rotation matrix for the yBins - anticyclic
188  RotationMatrix3D yBinRotationMatrix;
189  yBinRotationMatrix.col(0) = Vector3D::UnitX();
190  yBinRotationMatrix.col(1) = Vector3D::UnitZ();
191  yBinRotationMatrix.col(2) = Vector3D(0., -1., 0.);
192  // easy stuff first, constant pitch in Y
193  double pitchY =
194  2. * m_activeBounds->boundingBox().halfLengthY() / m_binUtility->bins(1);
195  // let's create the shared bounds of all surfaces marking y bins
196  std::shared_ptr<const PlanarBounds> yBinBounds(new RectangleBounds(
197  m_activeBounds->boundingBox().halfLengthX(), halfThickness));
198  // reserve, it's always (number of bins-1) as the boundaries are within the
199  // boundarySurfaces
200  segmentationSurfacesY.reserve(m_binUtility->bins(1));
201  for (size_t ibiny = 0; ibiny <= m_binUtility->bins(1); ++ibiny) {
202  // the position of the bin surface
203  double binPosY =
204  -m_activeBounds->boundingBox().halfLengthY() + ibiny * pitchY;
205  Vector3D binSurfaceCenter(0., binPosY, 0.);
206  // the binning transform
207  auto binTransform = std::make_shared<const Transform3D>(
208  Translation3D(binSurfaceCenter) * yBinRotationMatrix);
209  // these are the boundaries
210  if (ibiny == 0 || ibiny == m_binUtility->bins(1)) {
211  boundarySurfaces.push_back(
212  Surface::makeShared<PlaneSurface>(binTransform, yBinBounds));
213  } else { // these are the bin boundaries
214  segmentationSurfacesY.push_back(
215  Surface::makeShared<PlaneSurface>(binTransform, yBinBounds));
216  }
217  }
218 }
219 
221  const DigitizationCell& dCell) const {
222  double bX = m_binUtility->bins(0) > 1
223  ? m_binUtility->binningData()[0].center(dCell.channel0)
224  : 0.;
225  double bY = m_binUtility->bins(1) > 1
226  ? m_binUtility->binningData()[1].center(dCell.channel1)
227  : 0.;
228  return Vector2D(bX, bY);
229 }
230 
234  const Vector3D& startStep, const Vector3D& endStep, double halfThickness,
235  int readoutDirection, double lorentzAngle) const {
236  Vector3D stepCenter = 0.5 * (startStep + endStep);
237  // take the full drift length
238  // this is the absolute drift in z
239  double driftInZ = halfThickness - readoutDirection * stepCenter.z();
240  // this is the absolute drift length
241  double driftLength = driftInZ / cos(lorentzAngle);
242  // project to parameter the readout surface
243  double lorentzDeltaX = readoutDirection * driftInZ * tan(lorentzAngle);
244  // the projected center, it has the lorentz shift applied
245  Vector2D stepCenterProjected(stepCenter.x() + lorentzDeltaX, stepCenter.y());
246  // the cell & its center
247  Acts::DigitizationCell dCell = cell(stepCenterProjected);
248  Vector2D cellCenter = cellPosition(dCell);
249  // we are ready to return what we have
250  return DigitizationStep((endStep - startStep).norm(), driftLength, dCell,
251  startStep, endStep, stepCenterProjected, cellCenter);
252 }