2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
5 #ifndef ANALOGUECLUSTERINGTOOLIMPL_ICC
6 #define ANALOGUECLUSTERINGTOOLIMPL_ICC
10 #include <PixelReadoutGeometry/PixelModuleDesign.h>
11 #include "InDetReadoutGeometry/SiDetectorElement.h"
12 #include "InDetMeasurementUtilities/Helpers.h"
14 namespace ActsTrk::detail {
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)
24 template <typename calib_data_t, typename traj_t>
25 StatusCode AnalogueClusteringToolImpl<calib_data_t, traj_t>::initialize()
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");
34 ATH_MSG_DEBUG( m_thickness );
35 ATH_MSG_DEBUG( m_postCalibration );
36 ATH_MSG_DEBUG( m_correctCovariance );
37 ATH_MSG_DEBUG( m_calibratedCovarianceLowerBound );
39 return StatusCode::SUCCESS;
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
46 const InDetDD::SiDetectorElement* element=nullptr;
47 SG::ReadCondHandle<InDetDD::SiDetectorElementCollection> pixelDetEleHandle(
49 Gaudi::Hive::currentContext());
51 const InDetDD::SiDetectorElementCollection* detElements(*pixelDetEleHandle);
53 if (!pixelDetEleHandle.isValid() or detElements == nullptr) {
54 ATH_MSG_ERROR(m_pixelDetEleCollKey.fullKey() << " is not available.");
57 element = detElements->getDetectorElement(id);
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");
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
73 float projPhi = direction.dot(element.phiAxis());
74 float projEta = direction.dot(element.etaAxis());
75 float projNorm = direction.dot(element.normal());
77 float anglePhi = std::atan2(projPhi, projNorm);
78 float angleEta = std::atan2(projEta, projNorm);
80 // Map the angles of inward-going tracks onto [-PI/2, PI/2]
81 if (anglePhi > M_PI *0.5) {
84 if (anglePhi < -M_PI *0.5) {
88 // settle the sign/pi periodicity issues
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;
97 angleEta = -1 * log(tan(thetaloc * 0.5));
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);
104 return std::make_pair(anglePhi, angleEta);
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
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),
126 // if weighted position, we need to compute the centroid
127 const std::vector<Identifier>& rdos = cluster.rdoList();
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);
138 int col = m_pixelid->eta_index(rid);
139 colmin = std::min(col, colmin);
140 colmax = std::max(col, colmax);
143 const InDetDD::PixelModuleDesign& design =
144 dynamic_cast<const InDetDD::PixelModuleDesign&>(element.design());
146 InDetDD::SiLocalPosition pos1 =
147 design.positionFromColumnRow(colmin, rowmin);
149 InDetDD::SiLocalPosition pos2 =
150 design.positionFromColumnRow(colmax, rowmin);
152 InDetDD::SiLocalPosition pos3 =
153 design.positionFromColumnRow(colmin, rowmax);
155 InDetDD::SiLocalPosition pos4 =
156 design.positionFromColumnRow(colmax, rowmax);
158 InDetDD::SiLocalPosition centroid = 0.25 * (pos1 + pos2 + pos3 + pos4);
160 double shift = m_lorentzAngleTool->getLorentzShift(element.identifyHash(), Gaudi::Hive::currentContext());
162 return std::make_pair(centroid.xPhi() + shift, centroid.xEta());
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
169 SG::ReadCondHandle<calib_data_t> handle(
171 Gaudi::Hive::currentContext());
173 if (!handle.isValid()) {
174 ATH_MSG_ERROR(m_clusterErrorKey << " is not available.");
178 const error_data_t* data = handle->getClusterErrorData();
179 if (data == nullptr) {
180 ATH_MSG_ERROR("No cluster error data corresponding to " << m_clusterErrorKey);
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
195 auto& [anglePhi, angleEta] = angles;
197 std::optional<float> posX {std::nullopt};
198 std::optional<float> posY {std::nullopt};
200 int nrows = cluster.channelsInPhi();
201 int ncols = cluster.channelsInEta();
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,
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);
211 IdentifierHash idHash = element.identifyHash();
213 std::pair<double, double> delta =
214 errorData.getDelta(idHash,
221 posX = centroid.first + delta.first * (omX - 0.5);
224 posY = centroid.second + delta.second * (omY - 0.5);
225 } // check on omega values
226 } // check on cluster sizes
228 return std::make_pair(posX, posY);
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
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);
245 int nrows = cluster.channelsInPhi();
246 int ncols = cluster.channelsInEta();
248 auto& [anglePhi, angleEta] = angles;
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);
259 errY = *errY / std::sqrt(12);
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)
266 if (ncols > 1 and vY > 0)
270 return std::make_pair(errX, errY);
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
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));
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
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));
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
315 Pos pos = cluster.template localPosition<2>();
316 Cov cov = cluster.template localCovariance<2>();
318 assert(!cluster.rdoList().empty());
320 const error_data_t *errorData = getErrorData();
321 if (errorData == nullptr) {
322 throw std::runtime_error("PixelClusterErrorData is NULL");
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;
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);
337 if (errY.has_value()) {
338 double newErr = *errY;
339 cov(1, 1) = std::max(newErr * newErr, cov(1, 1) * m_calibratedCovarianceLowerBound);
342 return std::make_pair(pos, cov);
345 template <typename calib_data_t, typename traj_t>
346 void AnalogueClusteringToolImpl<calib_data_t, traj_t>::connect(OnTrackCalibrator<traj_t>& calibrator) const
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;
356 calibrator.pixelCalibrator. template connect<static_cast<CalibFuncPtr_t>(&AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate)>(this);
359 template <typename calib_data_t, typename traj_t>
360 void AnalogueClusteringToolImpl<calib_data_t, traj_t>::connectPixelCalibrator(IOnBoundStateCalibratorTool::PixelCalibrator& pixelCalibrator) const
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;
370 pixelCalibrator.template connect<static_cast<CalibFuncPtr_t>(&AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate)>(this);
373 template <typename calib_data_t, typename traj_t>
374 bool AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrateAfterMeasurementSelection() const
375 { return m_postCalibration.value(); }
377 } // namespace ActsTrk