2 Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
5 #ifndef ANALOGUECLUSTERINGTOOLIMPL_ICC
6 #define ANALOGUECLUSTERINGTOOLIMPL_ICC
10 #include <PixelReadoutGeometry/PixelModuleDesign.h>
14 template <typename calib_data_t, typename traj_t>
15 AnalogueClusteringToolImpl<calib_data_t, traj_t>::AnalogueClusteringToolImpl(
16 const std::string& type,
17 const std::string& name,
18 const IInterface* parent)
19 : base_class(type, name, parent)
22 template <typename calib_data_t, typename traj_t>
23 StatusCode AnalogueClusteringToolImpl<calib_data_t, traj_t>::initialize()
25 ATH_MSG_DEBUG("Initializing " << AthAlgTool::name() << " ...");
26 ATH_CHECK(m_pixelDetEleCollKey.initialize());
27 ATH_CHECK(m_clusterErrorKey.initialize());
28 ATH_CHECK(m_lorentzAngleTool.retrieve());
29 ATH_CHECK(AthAlgTool::detStore()->retrieve(m_pixelid, "PixelID"));
30 ATH_MSG_DEBUG(AthAlgTool::name() << " successfully initialized");
31 return StatusCode::SUCCESS;
34 template <typename calib_data_t, typename traj_t>
35 const InDetDD::SiDetectorElement&
36 AnalogueClusteringToolImpl<calib_data_t, traj_t>::getDetectorElement(xAOD::DetectorIDHashType id) const
38 const InDetDD::SiDetectorElement* element=nullptr;
39 SG::ReadCondHandle<InDetDD::SiDetectorElementCollection> pixelDetEleHandle(
41 Gaudi::Hive::currentContext());
43 const InDetDD::SiDetectorElementCollection* detElements(*pixelDetEleHandle);
45 if (!pixelDetEleHandle.isValid() or detElements == nullptr) {
46 ATH_MSG_ERROR(m_pixelDetEleCollKey.fullKey() << " is not available.");
49 element = detElements->getDetectorElement(id);
51 if (element == nullptr) {
52 ATH_MSG_ERROR("No element corresponding to hash " << id << " for " << m_pixelDetEleCollKey.fullKey());
53 throw std::runtime_error("SiDetectorElement is NULL");
59 template <typename calib_data_t, typename traj_t>
60 std::pair<float, float>
61 AnalogueClusteringToolImpl<calib_data_t, traj_t>::anglesOfIncidence(const InDetDD::SiDetectorElement& element,
62 const Acts::Vector3& direction) const
65 float projPhi = direction.dot(element.phiAxis());
66 float projEta = direction.dot(element.etaAxis());
67 float projNorm = direction.dot(element.normal());
69 float anglePhi = std::atan2(projPhi, projNorm);
70 float angleEta = std::atan2(projEta, projNorm);
72 // Map the angles of inward-going tracks onto [-PI/2, PI/2]
73 if (anglePhi > M_PI *0.5) {
76 if (anglePhi < -M_PI *0.5) {
80 // settle the sign/pi periodicity issues
82 if (angleEta > -0.5 * M_PI && angleEta < M_PI / 2.) {
83 thetaloc = M_PI_2 - angleEta;
84 } else if (angleEta > M_PI_2 && angleEta < M_PI) {
85 thetaloc = 1.5 * M_PI - angleEta;
86 } else { // 3rd quadrant
87 thetaloc = -M_PI_2 - angleEta;
89 angleEta = -1 * log(tan(thetaloc * 0.5));
92 // Subtract the Lorentz angle effect
93 float angleShift = m_lorentzAngleTool->getTanLorentzAngle(element.identifyHash());
94 anglePhi = std::atan(std::tan(anglePhi) - element.design().readoutSide() * angleShift);
96 return std::make_pair(anglePhi, angleEta);
100 template <typename calib_data_t, typename traj_t>
101 std::pair<float, float>
102 AnalogueClusteringToolImpl<calib_data_t, traj_t>::getCentroid(
103 const std::vector<Identifier>& rdos,
104 const InDetDD::SiDetectorElement& element) const
107 int rowmin = std::numeric_limits<int>::max();
108 int rowmax = std::numeric_limits<int>::min();
109 int colmin = std::numeric_limits<int>::max();
110 int colmax = std::numeric_limits<int>::min();
111 for (const Identifier& rid : rdos) {
112 int row = m_pixelid->phi_index(rid);
113 rowmin = std::min(row, rowmin);
114 rowmax = std::max(row, rowmax);
116 int col = m_pixelid->eta_index(rid);
117 colmin = std::min(col, colmin);
118 colmax = std::max(col, colmax);
121 const InDetDD::PixelModuleDesign& design =
122 dynamic_cast<const InDetDD::PixelModuleDesign&>(element.design());
124 InDetDD::SiLocalPosition pos1 =
125 design.positionFromColumnRow(colmin, rowmin);
127 InDetDD::SiLocalPosition pos2 =
128 design.positionFromColumnRow(colmax, rowmin);
130 InDetDD::SiLocalPosition pos3 =
131 design.positionFromColumnRow(colmin, rowmax);
133 InDetDD::SiLocalPosition pos4 =
134 design.positionFromColumnRow(colmax, rowmax);
136 InDetDD::SiLocalPosition centroid = 0.25 * (pos1 + pos2 + pos3 + pos4);
138 double shift = m_lorentzAngleTool->getLorentzShift(element.identifyHash());
140 return std::make_pair(centroid.xPhi() + shift, centroid.xEta());
143 template <typename calib_data_t, typename traj_t>
144 const typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::error_data_t*
145 AnalogueClusteringToolImpl<calib_data_t, traj_t>::getErrorData() const
147 SG::ReadCondHandle<calib_data_t> handle(
149 Gaudi::Hive::currentContext());
151 if (!handle.isValid()) {
152 ATH_MSG_ERROR(m_clusterErrorKey << " is not available.");
156 const error_data_t* data = handle->getClusterErrorData();
157 if (data == nullptr) {
158 ATH_MSG_ERROR("No cluster error data corresponding to " << m_clusterErrorKey);
165 template <typename calib_data_t, typename traj_t>
166 std::pair<std::optional<float>, std::optional<float>>
167 AnalogueClusteringToolImpl<calib_data_t, traj_t>::getCorrectedPosition(
168 const std::vector<Identifier>& rdos,
169 const error_data_t& errorData,
170 const InDetDD::SiDetectorElement& element,
171 const std::pair<float, float>& angles,
172 const xAOD::PixelCluster& cluster) const
174 // TODO validate these angles
175 auto& [anglePhi, angleEta] = angles;
177 std::optional<float> posX;
178 std::optional<float> posY;
180 float omX = cluster.omegaX();
181 float omY = cluster.omegaY();
182 int nrows = cluster.channelsInPhi();
183 int ncols = cluster.channelsInEta();
185 if (omX > -0.5 && omY > -0.5 && (nrows > 1 || ncols > 1)) {
187 std::pair<float, float> centroid = getCentroid(rdos, element);
189 Identifier id = element.identify();
190 std::pair<double, double> delta =
199 posX = centroid.first + delta.first * (omX - 0.5);
202 posY = centroid.second + delta.second * (omY - 0.5);
205 return std::make_pair(posX, posY);
208 template <typename calib_data_t, typename traj_t>
209 std::pair<std::optional<float>, std::optional<float>>
210 AnalogueClusteringToolImpl<calib_data_t, traj_t>::getCorrectedError(
211 const error_data_t& errorData,
212 const InDetDD::SiDetectorElement& element,
213 const std::pair<float, float>& angles,
214 const xAOD::PixelCluster& cluster) const
216 std::optional<float> errX;
217 std::optional<float> errY;
219 int nrows = cluster.channelsInPhi();
220 int ncols = cluster.channelsInEta();
222 auto& [anglePhi, angleEta] = angles;
224 // Special case for very shallow tracks
225 // Error estimated from geometrical projection of
226 // the track path in silicon onto the module surface
227 if (std::abs(anglePhi) > 1) {
228 errX = m_thickness * Acts::UnitConstants::um * std::tan(std::abs(anglePhi)) / std::sqrt(12);
229 errY = m_thickness * Acts::UnitConstants::um * std::tan(std::abs(angleEta));
230 if (cluster.widthInEta() > errY) {
231 errY = cluster.widthInEta() / std::sqrt(12);
233 errY = *errY / std::sqrt(12);
235 } else if (nrows > 1 && ncols > 1) {
236 Identifier id = element.identify();
237 auto [vX, vY] = errorData.getDeltaError(&id);
244 return std::make_pair(errX, errY);
247 template <typename calib_data_t, typename traj_t>
248 std::pair<typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Pos,
249 typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Cov>
250 AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate(
251 const Acts::GeometryContext& gctx,
252 const Acts::CalibrationContext& cctx,
253 const xAOD::PixelCluster& cluster,
254 const TrackStateProxy& state) const
256 const InDetDD::SiDetectorElement &detElement = getDetectorElement(cluster.identifierHash());
257 Acts::Vector3 direction = Acts::makeDirectionFromPhiTheta(
258 state.parameters()[Acts::eBoundPhi],
259 state.parameters()[Acts::eBoundTheta]);
260 return calibrate(gctx, cctx, cluster, detElement, anglesOfIncidence(detElement, direction));
263 template <typename calib_data_t, typename traj_t>
264 std::pair<typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Pos,
265 typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Cov>
266 AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate(
267 const Acts::GeometryContext& gctx,
268 const Acts::CalibrationContext& cctx,
269 const xAOD::PixelCluster& cluster,
270 const Acts::BoundTrackParameters& bound_parameters) const
272 const InDetDD::SiDetectorElement &detElement = getDetectorElement(cluster.identifierHash());
273 Acts::Vector3 direction = Acts::makeDirectionFromPhiTheta(
274 bound_parameters.parameters()[Acts::eBoundPhi],
275 bound_parameters.parameters()[Acts::eBoundTheta]);
276 return calibrate(gctx, cctx, cluster, detElement, anglesOfIncidence(detElement, direction));
279 template <typename calib_data_t, typename traj_t>
280 std::pair<typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Pos,
281 typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Cov>
282 AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate(
283 const Acts::GeometryContext& /*gctx*/,
284 const Acts::CalibrationContext& /*cctx*/,
285 const xAOD::PixelCluster& cluster,
286 const InDetDD::SiDetectorElement& detElement,
287 const std::pair<float, float>& angles) const
289 Pos pos = cluster.template localPosition<2>();
290 Cov cov = cluster.template localCovariance<2>();
292 assert(!cluster.rdoList().empty());
294 const error_data_t *errorData = getErrorData();
295 if (errorData == nullptr) {
296 throw std::runtime_error("PixelClusterErrorData is NULL");
300 auto [posX, posY] = getCorrectedPosition(cluster.rdoList(), *errorData, detElement, angles, cluster);
301 if (posX.has_value())
302 pos[Acts::eBoundLoc0] = *posX;
303 if (posY.has_value())
304 pos[Acts::eBoundLoc1] = *posY;
306 auto [errX, errY] = getCorrectedError(*errorData, detElement, angles, cluster);
307 if (errX.has_value())
308 cov(0, 0) = (*errX) * (*errX);
309 if (errY.has_value())
310 cov(1, 1) = (*errY) * (*errY);
312 return std::make_pair(pos, cov);
315 template <typename calib_data_t, typename traj_t>
316 void AnalogueClusteringToolImpl<calib_data_t, traj_t>::connect(OnTrackCalibrator<traj_t>& calibrator) const
318 using CalibFuncPtr_t =
319 std::pair<typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Pos,
320 typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Cov>
321 (AnalogueClusteringToolImpl<calib_data_t, traj_t>:: *)(const Acts::GeometryContext&,
322 const Acts::CalibrationContext&,
323 const xAOD::PixelCluster&,
324 const TrackStateProxy&) const;
326 calibrator.pixel_calibrator. template connect<static_cast<CalibFuncPtr_t>(&AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate)>(this);
329 template <typename calib_data_t, typename traj_t>
330 void AnalogueClusteringToolImpl<calib_data_t, traj_t>::connectPixelCalibrator(IOnBoundStateCalibratorTool::PixelCalibrator& pixel_calibrator) const
332 using CalibFuncPtr_t =
333 std::pair<typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Pos,
334 typename AnalogueClusteringToolImpl<calib_data_t, traj_t>::Cov>
335 (AnalogueClusteringToolImpl<calib_data_t, traj_t>:: *)(const Acts::GeometryContext&,
336 const Acts::CalibrationContext&,
337 const xAOD::PixelCluster&,
338 const Acts::BoundTrackParameters&) const;
340 pixel_calibrator. template connect<static_cast<CalibFuncPtr_t>(&AnalogueClusteringToolImpl<calib_data_t, traj_t>::calibrate)>(this);
343 } // namespace ActsTrk