ATLAS Offline Software
ITkStripClusterOnTrackTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // Implementation file for class ITkStripClusterOnTrackTool
8 // (c) ATLAS Detector software
10 // AlgTool used for ITkStripClusterOnTrack object production
12 
14 
20 
21 #include <cmath>
22 
23 using CLHEP::micrometer;
24 using CLHEP::deg;
25 
27 // Constructor
29 
31  (const std::string &t, const std::string &n, const IInterface *p) :
32  AthAlgTool(t, n, p) {
33  declareInterface<IRIO_OnTrackCreator>(this);
34 }
35 
37 // Initialisation
39 
43 
44  ATH_MSG_DEBUG("Error strategy set to ");
45  switch (m_option_errorStrategy) {
46  case -1: ATH_MSG_DEBUG("keep the PRD errors");
47  break;
48 
49  case 0: ATH_MSG_DEBUG("apply width/sqrt(12) as errors");
50  break;
51 
52  default: ATH_MSG_ERROR(" -- NO, UNKNOWN. Pls check jobOptions!");
53  return StatusCode::FAILURE;
54  }
55  ATH_MSG_DEBUG(" will be applied during ITkStripClusterOnTrack making");
56 
57  ATH_MSG_DEBUG("Position correction strategy set to ");
58  switch (m_option_correctionStrategy) {
59  case -1: ATH_MSG_DEBUG("keep the global position as evaluated");
60  break;
61 
62  default: ATH_MSG_ERROR(" -- NO, UNKNOWN. Pls check jobOptions!");
63  return StatusCode::FAILURE;
64  }
65  ATH_MSG_DEBUG(" will be applied during ITkStripClusterOnTrack making");
66 
67 
68  // get the error scaling tool
69  if (!m_stripErrorScalingKey.key().empty()) {
70  ATH_CHECK(m_stripErrorScalingKey.initialize());
71  ATH_MSG_DEBUG("Detected need for scaling ITkStrip errors.");
72  }
73 
74  return sc;
75 }
76 
78 // Trk::ITkStripClusterOnTrack production
80 
83  (const Trk::PrepRawData &rio, const Trk::TrackParameters &trackPar, const EventContext& ctx ) const {
84  const InDet::SCT_Cluster *cluster = nullptr;
85 
86  if (!(cluster = dynamic_cast<const InDet::SCT_Cluster *> (&rio))) {
87  ATH_MSG_WARNING("Attempt to correct RIO which is not SCT_Cluster with ITk::StripClusterOnTrackTool: returning nullptr");
88  return nullptr;
89  }
90 
91  ATH_MSG_VERBOSE("STARTING CLUSTER ON TRACK CORRECTION... " << __func__ << " " << __LINE__);
92  ATH_MSG_VERBOSE(" DUMPING CLUSTER POSITION / COVARIANCE: ");
93  ATH_MSG_VERBOSE("STRIP CLUSTER POS --> " << cluster->localPosition()[0] << ", " << cluster->localPosition()[1]);
94  ATH_MSG_VERBOSE("STRIP CLUSTER COV --> " << cluster->localCovariance()(0, 0) << ", " << cluster->localCovariance()(0, 1));
95  ATH_MSG_VERBOSE("STRIP CLUSTER COV --> " << cluster->localCovariance()(1, 0) << ", " << cluster->localCovariance()(1, 1));
96  ATH_MSG_VERBOSE("STRIP CLUSTER GLOBAL POSITION = " << cluster->globalPosition().x() << ", " << cluster->globalPosition().y() << ", " << cluster->globalPosition().z());
97 
98  // Get pointer to detector element
99  //
100  const InDetDD::SiDetectorElement *detectorElement = cluster->detectorElement();
101  if (!detectorElement) {
102  return nullptr;
103  }
104 
105  // Get local position of track
106  //
107  Amg::Vector2D loct = trackPar.localPosition();
108  double sinAlpha = detectorElement->sinStereoLocal(cluster->localPosition());
109  double cosAlpha = std::sqrt(1 - sinAlpha * sinAlpha);
110  Amg::Vector3D localstripdir(-sinAlpha, cosAlpha, 0.);
111  ATH_MSG_VERBOSE("STRIP DIRECTION = " << localstripdir[0] << ", " << localstripdir[1]);
112  Amg::Vector3D globalstripdir = trackPar.associatedSurface().transform().linear() * localstripdir;
113 
114  // Evaluate distance between cluster and estimated track parameters
115  // used later on to estimate the corrected cluster position
116  double distance = (trackPar.position() - cluster->globalPosition()).mag();
117 
118  ATH_MSG_VERBOSE(" DUMPING TRACK PARAMETER POSITION / COVARIANCE: ");
119  ATH_MSG_VERBOSE("TRACK PAR LOCAL POS = " << loct[0] << ", " << loct[1]);
120  ATH_MSG_VERBOSE("TRACK PAR GLOBAL POSITION = " << trackPar.position().x() << ", " << trackPar.position().y() << ", " << trackPar.position().z());
121 
122  // phi pitch in radians, for endcap modules
123  double phiPitchInRad = 0.;
124 
125  // barrel or endcap treatment
126  if (detectorElement->isBarrel()) {
127  // barrel treatment:
128  // get the module half length from the associated surface bounds
129  const Trk::SurfaceBounds *bounds = &trackPar.associatedSurface().bounds();
130  double boundsy = (static_cast<const Trk::RectangleBounds *>(bounds))->halflengthY();
131  ATH_MSG_VERBOSE("BARREL ====>>>> DISTANCE*COSALPHA / HALF LENGTH --> " << distance*cosAlpha << " / " << boundsy);
132  // Check if distance between track parameter local position
133  // and cluster position is larger than surface bounds (including local stereo angle).
134  // If so, set distance to maximum (- tolerance)
135  if (distance*cosAlpha > boundsy){
136  ATH_MSG_VERBOSE("DISTANCE TO LARGE COMPARED TO BOUNDS, SETTING TO MAXIMUM");
137  distance = boundsy/cosAlpha - 1.; // use 1 mm as tolerance parameter
138  // if local position is negative, also the distance has to be negative
139  if (loct.y() < 0)
140  distance = -distance;
141  }
142  } else {
143  // endcap treatment:
144  // for annuli do something different, since we already have in-sensor
145  // stereo rotations which strip length accounts for
146  const InDetDD::StripStereoAnnulusDesign * design =
147  static_cast<const InDetDD::StripStereoAnnulusDesign *> (&detectorElement->design());
148  const InDetDD::SiCellId & siCellId = detectorElement->cellIdOfPosition(cluster->localPosition());
149  double striphalflength = design->stripLength(siCellId) / 2.0;
150  ATH_MSG_VERBOSE("ENDCAP ====>>>> DISTANCE / STRIP HALF LENGTH --> " << distance << " / " << striphalflength);
151  // caching phi pitch in radians
152  phiPitchInRad = design->phiWidth()/design->diodesInRow(0);
153  // Check if distance between track parameter local position
154  // and cluster position is larger than strip length.
155  // If so, set distance to maximum (- tolerance)
156  if (distance > striphalflength) {
157  ATH_MSG_VERBOSE("DISTANCE TO LARGE COMPARED TO BOUNDS, SETTING TO MAXIMUM");
158  distance = striphalflength - 1.; // use 1 mm as tolerance parameter
159  }
160  }
161 
162  Amg::MatrixX prevCov = cluster->localCovariance();
163  // Local position and error matrix production
164  //
165  constexpr double ONE_TWELFTH= 1.0/12.;
166  if (m_option_errorStrategy > -1) {
167  Amg::MatrixX mat(2, 2);
168  mat.setZero();
169  const InDet::SiWidth width = cluster->width();
170  switch (m_option_errorStrategy) {
171  case 0:
172  // apply width/sqrt(12) as errors
173  mat(0, 0) = std::pow(width.phiR(), 2) * ONE_TWELFTH;
174  mat(1, 1) = std::pow(width.z(), 2) * ONE_TWELFTH;
175  break;
176 
177  default:
178  // don't do anything....
179  break;
180  }
181 
182  ATH_MSG_VERBOSE("CLUSTER ON TRACK COVARIANCE = " << mat(0, 0) << ", " << mat(0, 1) );
183  ATH_MSG_VERBOSE(" " << mat(1, 0) << ", " << mat(1, 1) );
184 
185  // error matrix rotation for endcap clusters
186  if (not detectorElement->isBarrel()) {
187  // the cluster covariance is expressed in the module
188  // frame and needs to be rotated into the strip frame.
189  // This rotation is done using the jacobian:
190  //
191  // J = ( cos(alpha) sin(alpha)
192  // -sin(alpha) cos(alpha) )
193  //
194  // The rotated covariance becomes
195  // C = J C_{0} J^T
196  //
197  // where C_{0} is the cluster covariance,
198  // and alpha is the local stereo angle.
199  //
200  // C_{0} = ( var_{0}^2 0
201  // 0 var_{1}^2 )
202  //
203  // var_{0} also accounts for the pitch scaling.
204  //
205  double sinAlpha = detectorElement->sinStereoLocal(cluster->localPosition());
206  double sinAlpha2 = sinAlpha * sinAlpha;
207  double cosAlpha2 = (1. - sinAlpha) * (1. + sinAlpha);;
208  double weight = detectorElement->phiPitch(cluster->localPosition()) / detectorElement->phiPitch();
209  double v0 = mat(0, 0) * weight * weight;
210  double v1 = mat(1, 1);
211  mat(0, 0) = (cosAlpha2 * v0 + sinAlpha2 * v1);
212  mat(1, 0) = (sinAlpha * std::sqrt(cosAlpha2) * (v0 - v1));
213  mat(0, 1) = mat(1, 0);
214  mat(1, 1) = (sinAlpha2 * v0 + cosAlpha2 * v1);
215 
216  ATH_MSG_VERBOSE(" ROTATION OF ENDCAP STRIP CLUSTER COVARIANCE");
217  ATH_MSG_VERBOSE("sinAlpha / sinAlpha2 / cosAlpha2 / weight = " << sinAlpha << " / " << sinAlpha2 << " / " << cosAlpha2 << " / " << weight);
218  ATH_MSG_VERBOSE("v0 / v1 = " << v0 << " / " << v1);
219  ATH_MSG_VERBOSE("ROTATED CLUSTER COVARIANCE = " << mat(0, 0) << ", " << mat(0, 1) );
220  ATH_MSG_VERBOSE(" " << mat(1, 0) << ", " << mat(1, 1) );
221 
222  }
223  // updating covariance
224  prevCov = mat;
225  }
226 
227  Trk::LocalParameters localParameters;
228  Amg::MatrixX covariance(prevCov);
229  Amg::Vector3D globalPosition(cluster->globalPosition() + distance * globalstripdir);
230 
231  // construction of the local parameters to build cluster on track
232  if (detectorElement->isBarrel()) {
233  Trk::DefinedParameter lpos1dim(cluster->localPosition().x(), Trk::locX);
234  localParameters = Trk::LocalParameters(lpos1dim);
235  covariance = Amg::MatrixX(1, 1);
236  covariance(0, 0) = prevCov(0, 0);
237  // scaling errors if required
238  if (!m_stripErrorScalingKey.key().empty()) {
239  SG::ReadCondHandle<RIO_OnTrackErrorScaling> error_scaling( m_stripErrorScalingKey,ctx );
240  covariance = Trk::ErrorScalingCast<SCTRIO_OnTrackErrorScaling>(*error_scaling)
241  ->getScaledCovariance(std::move(covariance), false, 0.0);
242  }
243  } else {
244  localParameters = Trk::LocalParameters(cluster->localPosition());
245  // scaling errors if required
246  if (!m_stripErrorScalingKey.key().empty()) {
247  SG::ReadCondHandle<RIO_OnTrackErrorScaling> error_scaling(m_stripErrorScalingKey,ctx);
248  covariance = Trk::ErrorScalingCast<SCTRIO_OnTrackErrorScaling>(*error_scaling)
249  ->getScaledCovariance(std::move(covariance), true,
250  detectorElement->sinStereoLocal(cluster->localPosition()));
251  }
252 
253  // For endcap strip clusters, the error matrix needs to be scaled
254  // accordingly to the strip pitch (in mm) at the estimated track parameter.
255  // The scaling term is evaluated in the strip frame and
256  // transformed in the module frame.
257  double sinAlpha = detectorElement->sinStereoLocal(cluster->localPosition());
258  double sinAlpha2 = sinAlpha * sinAlpha;
259  double cosAlpha2 = (1. - sinAlpha) * (1. + sinAlpha);
260  double sinAlphaCosAlpha = sinAlpha * std::sqrt(cosAlpha2);
261  // Weight factor to express the strip pitch (in mm) from the module center to the
262  // estimated track parameter position.
263  double radiusAtLocPos = std::hypot(loct.x(), loct.y());
264  double phiPitchAtLocPos = phiPitchInRad*radiusAtLocPos;
265  double weight = phiPitchAtLocPos / detectorElement->phiPitch();
266  // Error matrix scaling term, this is expressend in the module frame.
267  // It is evaluated using the same jacobian used above.
268  double dV0 = (cosAlpha2 * covariance(0, 0) + sinAlpha2 * covariance(1, 1) +
269  2. * sinAlphaCosAlpha * covariance(1, 0)) * (weight * weight - 1.);
270 
271  // The final covariance matrix, which also accounts for the scaling:
272  // Scaling and rotation are done with the same transformation:
273  //
274  // C = J C_{0} J^T
275  //
276  // where C_{0} is also includes the scaling term dV0:
277  //
278  // C_{0} = ( var_{0}^2 + dV0 0
279  // 0 var_{1}^2 )
280  //
281  // The scaling terms are then added to the previously evaluated matrix:
282  covariance(0, 0) += (cosAlpha2 * dV0);
283  covariance(1, 0) += (sinAlphaCosAlpha * dV0);
284  covariance(0, 1) = covariance(1, 0);
285  covariance(1, 1) += (sinAlpha2 * dV0);
286 
287  ATH_MSG_VERBOSE(" SCALING OF ENDCAP STRIP CLUSTER COVARIANCE");
288  ATH_MSG_VERBOSE("sinAlpha / sinAlpha2 / cosAlpha2 / weight = " << sinAlpha << " / " << sinAlpha2 << " / " << cosAlpha2 << " / " << weight );
289  ATH_MSG_VERBOSE("dV0 = (" << cosAlpha2 * covariance(0, 0) << " + "
290  << sinAlpha2 * covariance(1, 1) << " + "
291  << 2. * sinAlphaCosAlpha * covariance(1, 0)
292  << " ) * " << (weight * weight - 1.) << " = "
293  << dV0);
294  ATH_MSG_VERBOSE("SCALED CLUSTER COVARIANCE = " << covariance(0, 0) << ", "
295  << covariance(0, 1));
296  ATH_MSG_VERBOSE(" " << covariance(1, 0) << ", " << covariance(1, 1) );
297  }
298  // final construction of clustr of track
299  bool isbroad = m_option_errorStrategy == 0;
300  return new InDet::SCT_ClusterOnTrack(cluster, std::move(localParameters), std::move(covariance),
301  detectorElement->identifyHash(), globalPosition, isbroad);
302 }
Trk::LocalParameters
Definition: LocalParameters.h:98
Trk::RectangleBounds
Definition: RectangleBounds.h:38
InDetDD::SolidStateDetectorElementBase::cellIdOfPosition
SiCellId cellIdOfPosition(const Amg::Vector2D &localPos) const
As in previous method but returns SiCellId.
Definition: SolidStateDetectorElementBase.cxx:224
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition: EventPrimitives.h:27
SG::ReadCondHandle
Definition: ReadCondHandle.h:44
RectangleBounds.h
Trk::locX
@ locX
Definition: ParamDefs.h:37
Trk::ParametersBase::position
const Amg::Vector3D & position() const
Access method for the position.
Amg::Vector2D
Eigen::Matrix< double, 2, 1 > Vector2D
Definition: GeoPrimitives.h:48
Trk::ParametersBase::associatedSurface
virtual const Surface & associatedSurface() const override=0
Access to the Surface associated to the Parameters.
Trk::SurfaceBounds
Definition: SurfaceBounds.h:47
mat
GeoMaterial * mat
Definition: LArDetectorConstructionTBEC.cxx:55
initialize
void initialize()
Definition: run_EoverP.cxx:894
Trk::PrepRawData::localCovariance
const Amg::MatrixX & localCovariance() const
return const ref to the error matrix
ITk::StripClusterOnTrackTool::initialize
virtual StatusCode initialize() override
AlgTool initialisation.
Definition: ITkStripClusterOnTrackTool.cxx:41
deg
#define deg
Definition: SbPolyhedron.cxx:17
InDetDD::StripStereoAnnulusDesign::phiWidth
double phiWidth() const
Definition: StripStereoAnnulusDesign.h:331
ITk::StripClusterOnTrackTool::StripClusterOnTrackTool
StripClusterOnTrackTool(const std::string &, const std::string &, const IInterface *)
AlgTool constructor.
Definition: ITkStripClusterOnTrackTool.cxx:31
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
InDetDD::SolidStateDetectorElementBase::identifyHash
virtual IdentifierHash identifyHash() const override final
identifier hash (inline)
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
Trk::DefinedParameter
std::pair< double, ParamDefs > DefinedParameter
Definition: DefinedParameter.h:27
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:189
InDetDD::StripStereoAnnulusDesign::diodesInRow
virtual int diodesInRow(const int row) const override
Definition: StripStereoAnnulusDesign.h:251
InDetDD::SiDetectorElement::phiPitch
double phiPitch() const
Pitch (inline methods)
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
AnnulusBounds.h
ITk::StripClusterOnTrackTool::correct
virtual InDet::SCT_ClusterOnTrack * correct(const Trk::PrepRawData &, const Trk::TrackParameters &, const EventContext &ctx=Gaudi::Hive::currentContext()) const override
produces an SCT_ClusterOnTrack for ITk strip clusters using the measured SCT_Cluster and the track pr...
Definition: ITkStripClusterOnTrackTool.cxx:83
parseMapping.v0
def v0
Definition: parseMapping.py:149
beamspotman.n
n
Definition: beamspotman.py:731
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
InDetDD::StripStereoAnnulusDesign
Definition: StripStereoAnnulusDesign.h:50
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Trk::ParametersBase
Definition: ParametersBase.h:55
python.SystemOfUnits.micrometer
int micrometer
Definition: SystemOfUnits.py:71
InDet::SCT_Cluster
Definition: InnerDetector/InDetRecEvent/InDetPrepRawData/InDetPrepRawData/SCT_Cluster.h:34
InDet::SiCluster::detectorElement
virtual const InDetDD::SiDetectorElement * detectorElement() const override final
return the detector element corresponding to this PRD The pointer will be zero if the det el is not d...
Trk::PrepRawData
Definition: PrepRawData.h:62
Trk::ParametersCommon::localPosition
Amg::Vector2D localPosition() const
Access method for the local coordinates, local parameter definitions differ for each surface type.
ITkStripClusterOnTrackTool.h
Trk::Surface::bounds
virtual const SurfaceBounds & bounds() const =0
Surface Bounds method.
Trk::PrepRawData::localPosition
const Amg::Vector2D & localPosition() const
return the local position reference
InDetDD::SiDetectorElement
Definition: SiDetectorElement.h:109
InDetDD::SiDetectorElement::isBarrel
bool isBarrel() const
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
ErrorScalingCast.h
InDetDD::SiCellId
Definition: SiCellId.h:29
StripStereoAnnulusDesign.h
InDet::SiCluster::globalPosition
const Amg::Vector3D & globalPosition() const
return global position reference
InDet::SiCluster::width
const InDet::SiWidth & width() const
return width class reference
Base_Fragment.width
width
Definition: Sherpa_i/share/common/Base_Fragment.py:59
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
InDetDD::SiDetectorElement::sinStereoLocal
double sinStereoLocal(const Amg::Vector2D &localPos) const
Angle of strip in local frame with respect to the etaAxis.
Definition: SiDetectorElement.cxx:288
InDet::SiWidth
Definition: SiWidth.h:25
ActsTrk::ONE_TWELFTH
constexpr double ONE_TWELFTH
Definition: StripClusteringTool.cxx:17
AthAlgTool
Definition: AthAlgTool.h:26
SiCellId.h
InDetDD::SiDetectorElement::design
virtual const SiDetectorDesign & design() const override final
access to the local description (inline):
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
Trk::Surface::transform
const Amg::Transform3D & transform() const
Returns HepGeom::Transform3D by reference.
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:26
InDet::SCT_ClusterOnTrack
Definition: SCT_ClusterOnTrack.h:44
InDetDD::StripStereoAnnulusDesign::stripLength
double stripLength(const SiCellId &cellId) const
Definition: StripStereoAnnulusDesign.cxx:618