ATLAS Offline Software
AnalogueClusteringToolImpl.icc
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #ifndef ANALOGUECLUSTERINGTOOLIMPL_ICC
6 #define ANALOGUECLUSTERINGTOOLIMPL_ICC
7 
8 #include <limits>
9 
10 #include <PixelReadoutGeometry/PixelModuleDesign.h>
11 #include "InDetReadoutGeometry/SiDetectorElement.h"
12 #include "InDetMeasurementUtilities/Helpers.h"
13 
14 namespace ActsTrk::detail {
15 
16 template <typename calib_data_t, typename traj_t>
17 AnalogueClusteringToolImpl<calib_data_t, traj_t>::AnalogueClusteringToolImpl(
18  const std::string& type,
19  const std::string& name,
20  const IInterface* parent)
21  : base_class(type, name, parent)
22 {}
23 
24 template <typename calib_data_t, typename traj_t>
25 StatusCode AnalogueClusteringToolImpl<calib_data_t, traj_t>::initialize()
26 {
27  ATH_MSG_DEBUG("Initializing " << AthAlgTool::name() << " ...");
28  ATH_CHECK(m_pixelDetEleCollKey.initialize());
29  ATH_CHECK(m_clusterErrorKey.initialize());
30  ATH_CHECK(m_lorentzAngleTool.retrieve());
31  ATH_CHECK(AthAlgTool::detStore()->retrieve(m_pixelid, "PixelID"));
32  ATH_MSG_DEBUG(AthAlgTool::name() << " successfully initialized");
33 
34  ATH_MSG_DEBUG( m_thickness );
35  ATH_MSG_DEBUG( m_postCalibration );
36  ATH_MSG_DEBUG( m_correctCovariance );
37  ATH_MSG_DEBUG( m_calibratedCovarianceLowerBound );
38 
39  return StatusCode::SUCCESS;
40 }
41 
42 template <typename calib_data_t, typename traj_t>
43 const InDetDD::SiDetectorElement&
44 AnalogueClusteringToolImpl<calib_data_t, traj_t>::getDetectorElement(xAOD::DetectorIDHashType id) const
45 {
46  const InDetDD::SiDetectorElement* element=nullptr;
47  SG::ReadCondHandle<InDetDD::SiDetectorElementCollection> pixelDetEleHandle(
48  m_pixelDetEleCollKey,
49  Gaudi::Hive::currentContext());
50 
51  const InDetDD::SiDetectorElementCollection* detElements(*pixelDetEleHandle);
52 
53  if (!pixelDetEleHandle.isValid() or detElements == nullptr) {
54  ATH_MSG_ERROR(m_pixelDetEleCollKey.fullKey() << " is not available.");
55  }
56  else {
57  element = detElements->getDetectorElement(id);
58  }
59  if (element == nullptr) {
60  ATH_MSG_ERROR("No element corresponding to hash " << id << " for " << m_pixelDetEleCollKey.fullKey());
61  throw std::runtime_error("SiDetectorElement is NULL");
62  }
63 
64  return *element;
65 }
66 
67 template <typename calib_data_t, typename traj_t>
68 std::pair<float, float>
69 AnalogueClusteringToolImpl<calib_data_t, traj_t>::anglesOfIncidence(const InDetDD::SiDetectorElement& element,
70  const Acts::Vector3& direction) const
71 {
72 
73  float projPhi = direction.dot(element.phiAxis());
74  float projEta = direction.dot(element.etaAxis());
75  float projNorm = direction.dot(element.normal());
76 
77  float anglePhi = std::atan2(projPhi, projNorm);
78  float angleEta = std::atan2(projEta, projNorm);
79 
80  // Map the angles of inward-going tracks onto [-PI/2, PI/2]
81  if (anglePhi > M_PI *0.5) {
82  anglePhi -= M_PI;
83  }
84  if (anglePhi < -M_PI *0.5) {
85  anglePhi += M_PI;
86  }
87 
88  // settle the sign/pi periodicity issues
89  float thetaloc;
90  if (angleEta > -0.5 * M_PI && angleEta < M_PI / 2.) {
91  thetaloc = M_PI_2 - angleEta;
92  } else if (angleEta > M_PI_2 && angleEta < M_PI) {
93  thetaloc = 1.5 * M_PI - angleEta;
94  } else { // 3rd quadrant
95  thetaloc = -M_PI_2 - angleEta;
96  }
97  angleEta = -1 * log(tan(thetaloc * 0.5));
98 
99 
100  // Subtract the Lorentz angle effect
101  float angleShift = m_lorentzAngleTool->getTanLorentzAngle(element.identifyHash(), Gaudi::Hive::currentContext());
102  anglePhi = std::atan(std::tan(anglePhi) - element.design().readoutSide() * angleShift);
103 
104  return std::make_pair(anglePhi, angleEta);
105 }
106 
107 
108 template <typename calib_data_t, typename traj_t>
109 std::pair<float, float>
110 AnalogueClusteringToolImpl<calib_data_t, traj_t>::getCentroid(
111  const xAOD::PixelCluster& cluster,
112  const InDetDD::SiDetectorElement& element) const
113 {
114  // Different behaviours if pixel uncalibrated cluster use digital or
115  // charge-weighted local position
116  if (not m_useWeightedPos) {
117  // We return the cluster local position here
118  // This is actually the digital position, which only minimally
119  // differs from the centroid
120  // Eventually we will have a NN for the on-track calibration instead this instead
121  const auto& localPos = cluster.template localPosition<2>();
122  return std::make_pair(localPos(0, 0),
123  localPos(1, 0));
124  }
125 
126  // if weighted position, we need to compute the centroid
127  const std::vector<Identifier>& rdos = cluster.rdoList();
128 
129  int rowmin = std::numeric_limits<int>::max();
130  int rowmax = std::numeric_limits<int>::min();
131  int colmin = std::numeric_limits<int>::max();
132  int colmax = std::numeric_limits<int>::min();
133  for (const Identifier& rid : rdos) {
134  int row = m_pixelid->phi_index(rid);
135  rowmin = std::min(row, rowmin);
136  rowmax = std::max(row, rowmax);
137 
138  int col = m_pixelid->eta_index(rid);
139  colmin = std::min(col, colmin);
140  colmax = std::max(col, colmax);
141  }
142 
143  const InDetDD::PixelModuleDesign& design =
144  dynamic_cast<const InDetDD::PixelModuleDesign&>(element.design());
145 
146  InDetDD::SiLocalPosition pos1 =
147  design.positionFromColumnRow(colmin, rowmin);
148 
149  InDetDD::SiLocalPosition pos2 =
150  design.positionFromColumnRow(colmax, rowmin);
151 
152  InDetDD::SiLocalPosition pos3 =
153  design.positionFromColumnRow(colmin, rowmax);
154 
155  InDetDD::SiLocalPosition pos4 =
156  design.positionFromColumnRow(colmax, rowmax);
157 
158  InDetDD::SiLocalPosition centroid = 0.25 * (pos1 + pos2 + pos3 + pos4);
159 
160  double shift = m_lorentzAngleTool->getLorentzShift(element.identifyHash(), Gaudi::Hive::currentContext());
161 
162  return std::make_pair(centroid.xPhi() + shift, centroid.xEta());
163 }
164 
165 template <typename calib_data_t, typename traj_t>
166 const typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::error_data_t*
167 AnalogueClusteringToolImpl<calib_data_t, traj_t>::getErrorData() const
168 {
169  SG::ReadCondHandle<calib_data_t> handle(
170  m_clusterErrorKey,
171  Gaudi::Hive::currentContext());
172 
173  if (!handle.isValid()) {
174  ATH_MSG_ERROR(m_clusterErrorKey << " is not available.");
175  return nullptr;
176  }
177 
178  const error_data_t* data = handle->getClusterErrorData();
179  if (data == nullptr) {
180  ATH_MSG_ERROR("No cluster error data corresponding to " << m_clusterErrorKey);
181  return nullptr;
182  }
183 
184  return data;
185 }
186 
187 template <typename calib_data_t, typename traj_t>
188 std::pair<std::optional<float>, std::optional<float>>
189 AnalogueClusteringToolImpl<calib_data_t, traj_t>::getCorrectedPosition(
190  const xAOD::PixelCluster& cluster,
191  const error_data_t& errorData,
192  const InDetDD::SiDetectorElement& element,
193  const std::pair<float, float>& angles) const
194 {
195  auto& [anglePhi, angleEta] = angles;
196 
197  std::optional<float> posX {std::nullopt};
198  std::optional<float> posY {std::nullopt};
199 
200  int nrows = cluster.channelsInPhi();
201  int ncols = cluster.channelsInEta();
202 
203  // If size of cluster (in any direction) is higher than 1, we may need to apply a correction
204  if (nrows > 1 or ncols > 1) {
205  const auto& [omX, omY] = TrackingUtilities::computeOmegas(cluster,
206  *m_pixelid);
207  // if the omegas are valid, we can apply the correction
208  if (omX > -0.5 and omY > -0.5) {
209  std::pair<float, float> centroid = getCentroid(cluster, element);
210 
211  IdentifierHash idHash = element.identifyHash();
212 
213  std::pair<double, double> delta =
214  errorData.getDelta(idHash,
215  nrows,
216  anglePhi,
217  ncols,
218  angleEta);
219 
220  if (nrows > 1)
221  posX = centroid.first + delta.first * (omX - 0.5);
222 
223  if (ncols > 1)
224  posY = centroid.second + delta.second * (omY - 0.5);
225  } // check on omega values
226  } // check on cluster sizes
227 
228  return std::make_pair(posX, posY);
229 }
230 
231 template <typename calib_data_t, typename traj_t>
232 std::pair<std::optional<float>, std::optional<float>>
233 AnalogueClusteringToolImpl<calib_data_t, traj_t>::getCorrectedError(
234  const error_data_t& errorData,
235  const InDetDD::SiDetectorElement& element,
236  const std::pair<float, float>& angles,
237  const xAOD::PixelCluster& cluster) const
238 {
239  std::optional<float> errX {std::nullopt};
240  std::optional<float> errY {std::nullopt};
241  if (not m_correctCovariance) {
242  return std::make_pair(errX, errY);
243  }
244 
245  int nrows = cluster.channelsInPhi();
246  int ncols = cluster.channelsInEta();
247 
248  auto& [anglePhi, angleEta] = angles;
249 
250  // Special case for very shallow tracks
251  // Error estimated from geometrical projection of
252  // the track path in silicon onto the module surface
253  if (std::abs(anglePhi) > 1) {
254  errX = m_thickness * Acts::UnitConstants::um * std::tan(std::abs(anglePhi)) / std::sqrt(12);
255  errY = m_thickness * Acts::UnitConstants::um * std::tan(std::abs(angleEta));
256  if (cluster.widthInEta() > errY) {
257  errY = cluster.widthInEta() / std::sqrt(12);
258  } else {
259  errY = *errY / std::sqrt(12);
260  }
261  } else if (nrows > 1 or ncols > 1) {
262  IdentifierHash idHash = element.identifyHash();
263  auto [vX, vY] = errorData.getDeltaError(idHash);
264  if (nrows > 1 and vX > 0)
265  errX = vX;
266  if (ncols > 1 and vY > 0)
267  errY = vY;
268  }
269 
270  return std::make_pair(errX, errY);
271 }
272 
273 template <typename calib_data_t, typename traj_t>
274 std::pair<typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Pos,
275  typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Cov>
276 AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate(
277  const Acts::GeometryContext& gctx,
278  const Acts::CalibrationContext& cctx,
279  const xAOD::PixelCluster& cluster,
280  const TrackStateProxy& state) const
281 {
282  const InDetDD::SiDetectorElement &detElement = getDetectorElement(cluster.identifierHash());
283  Acts::Vector3 direction = Acts::makeDirectionFromPhiTheta(
284  state.parameters()[Acts::eBoundPhi],
285  state.parameters()[Acts::eBoundTheta]);
286  return calibrate(gctx, cctx, cluster, detElement, anglesOfIncidence(detElement, direction));
287 }
288 
289 template <typename calib_data_t, typename traj_t>
290 std::pair<typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Pos,
291  typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Cov>
292 AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate(
293  const Acts::GeometryContext& gctx,
294  const Acts::CalibrationContext& cctx,
295  const xAOD::PixelCluster& cluster,
296  const Acts::BoundTrackParameters& bound_parameters) const
297 {
298  const InDetDD::SiDetectorElement &detElement = getDetectorElement(cluster.identifierHash());
299  Acts::Vector3 direction = Acts::makeDirectionFromPhiTheta(
300  bound_parameters.parameters()[Acts::eBoundPhi],
301  bound_parameters.parameters()[Acts::eBoundTheta]);
302  return calibrate(gctx, cctx, cluster, detElement, anglesOfIncidence(detElement, direction));
303 }
304 
305 template <typename calib_data_t, typename traj_t>
306 std::pair<typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Pos,
307  typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Cov>
308 AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate(
309  const Acts::GeometryContext& /*gctx*/,
310  const Acts::CalibrationContext& /*cctx*/,
311  const xAOD::PixelCluster& cluster,
312  const InDetDD::SiDetectorElement& detElement,
313  const std::pair<float, float>& angles) const
314 {
315  Pos pos = cluster.template localPosition<2>();
316  Cov cov = cluster.template localCovariance<2>();
317 
318  assert(!cluster.rdoList().empty());
319 
320  const error_data_t *errorData = getErrorData();
321  if (errorData == nullptr) {
322  throw std::runtime_error("PixelClusterErrorData is NULL");
323  }
324 
325 
326  auto [posX, posY] = getCorrectedPosition(cluster, *errorData, detElement, angles);
327  if (posX.has_value())
328  pos[Acts::eBoundLoc0] = *posX;
329  if (posY.has_value())
330  pos[Acts::eBoundLoc1] = *posY;
331 
332  auto [errX, errY] = getCorrectedError(*errorData, detElement, angles, cluster);
333  if (errX.has_value()) {
334  double newErr = *errX;
335  cov(0, 0) = std::max(newErr * newErr, cov(0, 0) * m_calibratedCovarianceLowerBound);
336  }
337  if (errY.has_value()) {
338  double newErr = *errY;
339  cov(1, 1) = std::max(newErr * newErr, cov(1, 1) * m_calibratedCovarianceLowerBound);
340  }
341 
342  return std::make_pair(pos, cov);
343 }
344 
345 template <typename calib_data_t, typename traj_t>
346 void AnalogueClusteringToolImpl<calib_data_t, traj_t>::connect(OnTrackCalibrator<traj_t>& calibrator) const
347 {
348  using CalibFuncPtr_t =
349  std::pair<typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Pos,
350  typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Cov>
351  (AnalogueClusteringToolImpl<calib_data_t, traj_t>:: *)(const Acts::GeometryContext&,
352  const Acts::CalibrationContext&,
353  const xAOD::PixelCluster&,
354  const TrackStateProxy&) const;
355 
356  calibrator.pixelCalibrator. template connect<static_cast<CalibFuncPtr_t>(&AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate)>(this);
357 }
358 
359 template <typename calib_data_t, typename traj_t>
360 void AnalogueClusteringToolImpl<calib_data_t, traj_t>::connectPixelCalibrator(IOnBoundStateCalibratorTool::PixelCalibrator& pixelCalibrator) const
361 {
362  using CalibFuncPtr_t =
363  std::pair<typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Pos,
364  typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Cov>
365  (AnalogueClusteringToolImpl<calib_data_t, traj_t>:: *)(const Acts::GeometryContext&,
366  const Acts::CalibrationContext&,
367  const xAOD::PixelCluster&,
368  const Acts::BoundTrackParameters&) const;
369 
370  pixelCalibrator.template connect<static_cast<CalibFuncPtr_t>(&AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate)>(this);
371 }
372 
373 template <typename calib_data_t, typename traj_t>
374 bool AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrateAfterMeasurementSelection() const
375 { return m_postCalibration.value(); }
376 
377 } // namespace ActsTrk
378 
379 #endif