ATLAS Offline Software
MuonSegmentInOverlapResolvingTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "MuonSegmentInOverlapResolvingTool.h"
6 
7 #include "GaudiKernel/MsgStream.h"
8 #include "GeoPrimitives/GeoPrimitivesHelpers.h"
9 #include "GeoPrimitives/GeoPrimitivesToStringConverter.h"
10 #include "MuonRIO_OnTrack/MdtDriftCircleOnTrack.h"
11 #include "MuonReadoutGeometry/MdtReadoutElement.h"
12 #include "MuonSegment/MuonSegment.h"
13 #include "TrkDetDescrUtils/Intersection.h"
14 #include "TrkEventPrimitives/LocalDirection.h"
15 #include "TrkEventPrimitives/ResidualPull.h"
16 #include "TrkParameters/TrackParameters.h"
17 
18 namespace Muon {
19 
20 MuonSegmentInOverlapResolvingTool::MuonSegmentInOverlapResolvingTool(const std::string& ty, const std::string& na,
21  const IInterface* pa)
22  : AthAlgTool(ty, na, pa), m_magFieldProperties(Trk::NoField)
23 {
24  declareInterface<IMuonSegmentInOverlapResolvingTool>(this);
25 }
26 
27 StatusCode
28 MuonSegmentInOverlapResolvingTool::initialize()
29 {
30  ATH_CHECK(AthAlgTool::initialize());
31 
32  ATH_CHECK(m_edmHelperSvc.retrieve());
33  ATH_CHECK(m_printer.retrieve());
34  ATH_CHECK(m_idHelperSvc.retrieve());
35  ATH_CHECK(m_propagator.retrieve());
36  ATH_CHECK(m_pullCalculator.retrieve());
37 
38  return StatusCode::SUCCESS;
39 }
40 
41 Amg::Vector3D
42 MuonSegmentInOverlapResolvingTool::updateSegmentDirection(const MuonSegment& seg, double phi) const
43 {
44 
45  Amg::Transform3D gToStation = seg.associatedSurface().transform().inverse();
46  Amg::Vector3D phiDir(1., 0., 0.);
47  Amg::setThetaPhi(phiDir, seg.globalDirection().theta(), phi);
48  Amg::Vector3D lphiDir = gToStation.linear() * phiDir;
49  Amg::Vector3D lsegDir = gToStation.linear() * seg.globalDirection();
50  double road_dz = lphiDir.z();
51  double seg_dy = lsegDir.y();
52  double seg_dz = lsegDir.z();
53  if (road_dz * seg_dz < 0) {
54  seg_dy *= -1.;
55  seg_dz *= -1.;
56  }
57  if (std::abs(seg_dz) < 1e-6) {
58  seg_dz = 1e-6;
59  ATH_MSG_DEBUG(" Unexpected local direction of segment " << lsegDir);
60  }
61  double scale = road_dz / seg_dz;
62  seg_dy *= scale;
63  Amg::Vector3D locDir(lphiDir.x(), seg_dy, road_dz);
64  Amg::Vector3D gDir = seg.associatedSurface().transform().linear() * locDir;
65  return gDir;
66 }
67 
68 Amg::Vector3D
69 MuonSegmentInOverlapResolvingTool::estimateSegmentDirection(const MuonSegment& seg1, const MuonSegment& seg2,
70  double& phi, double& stereoangle) const
71 {
72  //
73  // From two segments the direction and position of the segments are determined
74  // assuming a straight line
75  //
76  // if a stereoangle is present one can determine: phi and x local for both segments
77  //
78 
79  // bool check = false;
80 
81  // Confusing use of objects:
82  // This is needed because otherwise the Transform3D DOES NOT know
83  // whether it transforms a position or a direction
84  //
85  const Amg::Transform3D& gToGlobal1 = seg1.associatedSurface().transform();
86  const Amg::Transform3D& gToGlobal2 = seg2.associatedSurface().transform();
87  Amg::Transform3D gToLocal1 = seg1.associatedSurface().transform().inverse();
88  Amg::Transform3D gToLocal2 = seg2.associatedSurface().transform().inverse();
89  //
90  const Amg::Vector3D& gDir1 = seg1.globalDirection();
91  const Amg::Vector3D& gDir2 = seg2.globalDirection();
92  // Confusing but Amg::Vector3D can hold a 3D local direction
93  Amg::Vector3D lDir1 = gToLocal1.linear() * gDir1;
94  // Amg::Vector3D lDir2 = gToLocal2.linear()*gDir2;
95  // segment 2 local direction in frame of segment 1
96  Amg::Vector3D lDir12 = gToLocal1.linear() * gDir2;
97  // Amg::Vector3D lDir21 = gToLocal2.linear()*gDir1;
98 
99  Amg::Vector3D segLocDiro(1., 0., 0.);
100  // orthogonal vector to segment 2 in local frame of segment 1
101  Amg::Vector3D lDiro12 = gToLocal1.linear() * (gToGlobal2.linear() * segLocDiro);
102  // orthogonal vector to segment 1in local frame of segment 2
103  Amg::Vector3D lDiro21 = gToLocal2.linear() * (gToGlobal1.linear() * segLocDiro);
104 
105  stereoangle = std::acos(lDiro12.x());
106 
107  // We have the following equations for segment 2 in local frame of segment 1
108  // dx1 = free
109  // dy1 = lDir1.y() = a * lDir12.y() + b * lDiro12.y()
110  // dz1 = lDir1.z() = a * lDir12.z() + b * lDiro12.z()
111 
112  double b = lDir1.y() * lDir12.z() - lDir1.z() * lDir12.y();
113  double a = lDir1.y() * lDiro12.z() - lDir1.z() * lDiro12.y();
114  double dxn = lDir12.x();
115  double dyn = lDir12.y();
116  double dzn = lDir12.z();
117  if (std::abs(a) > 1e-2) {
118  dxn = lDir12.x() - b * lDiro12.x() / a;
119  dyn = lDir12.y() - b * lDiro12.y() / a;
120  dzn = lDir12.z() - b * lDiro12.z() / a;
121  }
122  double norm = std::hypot(dxn, dyn, dzn);
123  if (norm < 1e-6) {
124  ATH_MSG_DEBUG(" Unexpected normalisation " << norm);
125  norm = 1e-6;
126  }
127  if (dxn * lDir1.x() + dyn * lDir1.y() + dzn * lDir1.z() < 0) norm = -norm;
128  dxn = dxn / norm;
129  dyn = dyn / norm;
130  dzn = dzn / norm;
131 
132  // The final result for the direction of the two segments
133 
134  // New Local direction of segment 1
135  Amg::Vector3D lDirn1(dxn, dyn, dzn);
136  // Global direction of both segments
137  Amg::Vector3D gDirn = gToGlobal1.linear() * lDirn1;
138  // New Local direction of segment 2
139  Amg::Vector3D lDirn2 = gToLocal2.linear() * gDirn;
140  phi = gDirn.phi();
141  double theta = gDirn.theta();
142 
143  const Amg::Vector3D& gPos1 = seg1.globalPosition();
144  // Confusing but Amg::Vector3D can hold a 3D local position
145  Amg::Vector3D lPos1 = gToLocal1 * gPos1;
146  Amg::Vector3D lPos21 = gToLocal2 * gPos1;
147  const Amg::Vector3D& gPos2 = seg2.globalPosition();
148  Amg::Vector3D lPos12 = gToLocal1 * gPos2;
149  Amg::Vector3D lPos2 = gToLocal2 * gPos2;
150  // In local frame of segment 2 shift segment 1 to obtain zero residual
151  double res21 = (lPos2.y() - lPos21.y()) * lDirn2.z() - (lPos2.z() - lPos21.z()) * lDirn2.y();
152  double localx1 = 0.;
153  double step = (lDiro21.y() * lDirn2.z() - lDiro21.z() * lDirn2.y());
154  if (std::abs(step) > 1e-5) {
155  localx1 = res21 / step;
156  }
157  ATH_MSG_DEBUG(" localx1 " << localx1 << " res21 " << res21 << " step " << step);
158 
159  // Result: New local and global position of segment 1
160  Amg::Vector3D lPosn1(lPos1.x() + localx1, lPos1.y(), lPos1.z());
161  Amg::Vector3D gPosn1 = gToGlobal1 * lPosn1;
162  // In local frame of segment 1 shift segment 2 to obtain zero residual
163  double res12 = (lPos1.y() - lPos12.y()) * lDirn1.z() - (lPos1.z() - lPos12.z()) * lDirn1.y();
164  step = (lDiro12.y() * lDirn1.z() - lDiro12.z() * lDirn1.y());
165  double localx2 = 0.;
166  if (std::abs(step) > 1e-5) {
167  localx2 = res12 / step;
168  }
169  ATH_MSG_DEBUG(" localx2 " << localx2 << " res12 " << res12 << " step " << step);
170 
171  // Result: New local and global position of segment 2
172  Amg::Vector3D lPosn2(lPos2.x() + localx2, lPos2.y(), lPos2.z());
173  Amg::Vector3D gPosn2 = gToGlobal2 * lPosn2;
174  ATH_MSG_DEBUG(" segment 1 local position "
175  << lPos1 << " new " << lPosn1 << std::endl
176  << " segment 1 global position " << gPos1 << " new " << gPosn1 << std::endl
177  << " segment 2 local position " << lPos2 << " new " << lPosn2 << std::endl
178  << " segment 2 global position " << gPos2 << " new " << gPosn2);
179 
180  // This gives the direction from the position of the segments
181  Amg::Vector3D gDirPos = gPosn2 - gPosn1;
182  if (gDir1.x() * gDirPos.x() + gDir1.y() * gDirPos.y() + gDir1.z() * gDirPos.z() < 0) {
183  gDirPos = -gPosn2 + gPosn1;
184  }
185  double dtheta = theta - gDirPos.theta();
186  double dphi = phi - gDirPos.phi();
187  ATH_MSG_DEBUG(" theta " << theta << " gDirPos theta " << gDirPos.theta() << " dtheta " << dtheta << " phi " << phi
188  << " gDirPos phi " << gDirPos.phi() << " dphi " << dphi);
189  return gDirn;
190 }
191 
192 
193 MuonSegmentInOverlapResolvingTool::SegmentPhiMatchResult
194 MuonSegmentInOverlapResolvingTool::bestPhiMatchAnalytic(const MuonSegment& seg1, const MuonSegment& seg2) const
195 {
196 
197  //
198  // From two segments the direction and position of the segments are determined
199  // assuming a straight line
200  //
201  // if a stereoangle is present one can determine: phi and x local for both segments
202  //
203 
204  // bool check = false;
205 
206  // Confusing use of objects:
207  // Trk::GlobalDirection can and will hold a 3D LOCAL direction
208  // Trk::GlobalPositionn can and will hold a 3D LOCAL position
209  // This is needed because otherwise the Amg::Transform3D DOES NOT know
210  // whether it transforms a position or a direction
211  //
212  const Amg::Transform3D& gToGlobal1 = seg1.associatedSurface().transform();
213  const Amg::Transform3D& gToGlobal2 = seg2.associatedSurface().transform();
214  Amg::Transform3D gToLocal1 = seg1.associatedSurface().transform().inverse();
215  // Amg::Transform3D gToLocal2 = seg2.associatedSurface().transform().inverse();
216  //
217  const Amg::Vector3D& gDir1 = seg1.globalDirection();
218  const Amg::Vector3D& gDir2 = seg2.globalDirection();
219  // Confusing but Amg::Vector3D can hold a 3D local direction
220  Amg::Vector3D lDir1 = gToLocal1.linear() * gDir1;
221  // Amg::Vector3D lDir2 = gToLocal2.linear()*gDir2;
222  // segment 2 local direction in frame of segment 1
223  Amg::Vector3D lDir12 = gToLocal1.linear() * gDir2;
224  // Amg::Vector3D lDir21 = gToLocal2.linear()*gDir1;
225 
226  Amg::Vector3D segLocDiro(1., 0., 0.);
227  // orthogonal vector to segment 2 in local frame of segment 1
228  Amg::Vector3D lDiro12 = gToLocal1.linear() * (gToGlobal2.linear() * segLocDiro);
229  // orthogonal vector to segment 1in local frame of segment 2
230  // Amg::Vector3D lDiro21 = gToLocal2.linear()*(gToGlobal1.linear()*segLocDiro);
231 
232  // stereoangle = acos(lDiro12.x());
233 
234  // We have the following equations for segment 2 in local frame of segment 1
235  // dx1 = free
236  // dy1 = lDir1.y() = a * lDir12.y() + b * lDiro12.y()
237  // dz1 = lDir1.z() = a * lDir12.z() + b * lDiro12.z()
238 
239  double b = lDir1.y() * lDir12.z() - lDir1.z() * lDir12.y();
240  double a = lDir1.y() * lDiro12.z() - lDir1.z() * lDiro12.y();
241  double dxn = lDir12.x();
242  double dyn = lDir12.y();
243  double dzn = lDir12.z();
244  if (std::abs(a) > 1e-2) {
245  dxn = lDir12.x() - b * lDiro12.x() / a;
246  dyn = lDir12.y() - b * lDiro12.y() / a;
247  dzn = lDir12.z() - b * lDiro12.z() / a;
248  }
249  double norm = std::hypot(dxn, dyn, dzn);
250  if (norm < 1e-6) {
251  ATH_MSG_DEBUG(" Unexpected normalisation " << norm);
252  norm = 1e-6;
253  }
254  if (dxn * lDir1.x() + dyn * lDir1.y() + dzn * lDir1.z() < 0) norm = -norm;
255  dxn = dxn / norm;
256  dyn = dyn / norm;
257  dzn = dzn / norm;
258 
259  // The final result for the direction of the two segments
260 
261  // New Local direction of segment 1
262  Amg::Vector3D lDirn1(dxn, dyn, dzn);
263  // Global direction of both segments
264  Amg::Vector3D segDir1Min = gToGlobal1.linear() * lDirn1;
265  const Amg::Vector3D& segDir2Min = segDir1Min; // updateSegmentDirection(seg2,segDir1Min.phi());
266  Trk::LocalDirection segLocDir1;
267  seg2.associatedSurface().globalToLocalDirection(segDir1Min, segLocDir1);
268  Trk::LocalDirection segLocDir2;
269  seg2.associatedSurface().globalToLocalDirection(segDir2Min, segLocDir2);
270  double dyz = std::abs(segLocDir1.angleYZ() - segLocDir2.angleYZ());
271  return {segDir1Min, segDir2Min, dyz};
272 }
273 
274 MuonSegmentInOverlapResolvingTool::SegmentPhiMatchResult
275 MuonSegmentInOverlapResolvingTool::bestPhiMatch(const MuonSegment& seg1, const MuonSegment& seg2) const
276 {
277 
278  unsigned int nbins = 11;
279  // unsigned int bestBin = 0;
280  double scanRange = 1.;
281  double scanStep = scanRange / (nbins - 1);
282  double phiStart = seg1.globalDirection().phi() - 0.5 * scanRange;
283 
284  double dthetaMin = 1e9;
285  Amg::Vector3D segDir1Min = seg1.globalDirection();
286  Trk::LocalDirection segLocDir2;
287  seg2.associatedSurface().globalToLocalDirection(seg2.globalDirection(), segLocDir2);
288  for (unsigned int i = 0; i < nbins; ++i) {
289  double phi = phiStart + scanStep * i;
290  Amg::Vector3D segDir1 = updateSegmentDirection(seg1, phi);
291  Trk::LocalDirection segLocDir12;
292  seg2.associatedSurface().globalToLocalDirection(segDir1, segLocDir12);
293  double dyz = std::abs(segLocDir12.angleYZ() - segLocDir2.angleYZ());
294  if (dyz < dthetaMin) {
295  dthetaMin = dyz;
296  segDir1Min = segDir1;
297  // bestBin = i;
298  }
299  }
300  return {segDir1Min, updateSegmentDirection(seg2, segDir1Min.phi()), dthetaMin};
301 }
302 
303 
304 MuonSegmentInOverlapResolvingTool::SegmentPositionMatchResult
305 MuonSegmentInOverlapResolvingTool::bestPositionAlongTubeMatch(const MuonSegment& seg1, const MuonSegment& seg2,
306  const Amg::Vector3D& segDir1Min) const
307 {
308  bool goodMatch = true;
309 
310  // get geomtry summary for segment
311  SegmentGeometrySummary segmentGeometry = segmentGeometrySummary(seg1);
312  SegmentGeometrySummary segmentGeometry2 = segmentGeometrySummary(seg2);
313 
314  int nbins = 2;
315  double tubeStep = segmentGeometry.hvPosInSegFrame - segmentGeometry.roPosInSegFrame;
316 
317  double posStep = tubeStep / (nbins - 1);
318 
319  double resfirst{1e9}, reslast{1e9}, posfirst{1e9}, poslast{1e9};
320  double locy = seg1.localParameters().contains(Trk::locY) ? seg1.localParameters()[Trk::locY] : 0.;
321 
322  for (int j = 0; j < nbins; ++j) {
323  double distToRO = segmentGeometry.roPosInSegFrame + posStep * j;
324 
325  // create local position
326  Amg::Vector2D lpos(distToRO, locy);
327 
328  // segment reference surface
329  const Trk::PlaneSurface& segmentSurface = seg2.associatedSurface();
330 
331  // get global position
332  Amg::Vector3D gpos;
333  seg1.associatedSurface().localToGlobal(lpos, segDir1Min, gpos);
334 
335  // calculate intersection
336  Trk::Intersection intersect = segmentSurface.straightLineIntersection(gpos, segDir1Min, false, false);
337  if (!intersect.valid || !segmentSurface.globalToLocal(intersect.position, segDir1Min, lpos)) {
338  // if( !intersect.valid || !segmentSurface.globalToLocal(intersect.intersection,segDir1Min,lpos) ){
339  ATH_MSG_WARNING(" Intersect with surface position " << Amg::toString(gpos) << " direction: phi "
340  << segDir1Min.phi() << " theta "
341  << segDir1Min.theta());
342  } else {
343 
344  // now we have the prediction of the track in the segment frame, lets look at residuals
345  double resy = lpos[Trk::locY] - seg2.localParameters()[Trk::locY];
346 
347  if (j == 0) {
348  resfirst = resy;
349  posfirst = distToRO;
350  }
351  if (j == nbins - 1) {
352  reslast = resy;
353  poslast = distToRO;
354  }
355  }
356  }
357 
358  double distPosMin2{1e9}, distPosInTube2{1e9}, distPosMin{1e9}, distPosInTube{1e9}, resyMin{1e9};
359  Amg::Vector3D segPos(0., 0., 0.);
360  double rangeCut = 1e5;
361  if (resfirst < rangeCut && reslast < rangeCut && posfirst < rangeCut && poslast < rangeCut) {
362  double resDif = reslast - resfirst;
363  double posDif = poslast - posfirst;
364  if (std::abs(resDif) < 1e-6) {
365  ATH_MSG_DEBUG(" Unexpected residual difference " << resDif);
366  resDif = resDif < 0. ? -1e-6 : 1e-6;
367  }
368  if (std::abs(posDif) < 1e-6) {
369  ATH_MSG_DEBUG(" Unexpected position difference " << posDif);
370  posDif = posDif < 0. ? -1e-6 : 1e-6;
371  }
372  distPosMin = posfirst - resfirst / (resDif) * (posDif);
373  distPosInTube = segmentGeometry.positionInsideTube(distPosMin);
374  resyMin = resfirst + (resDif) * (distPosInTube - posfirst) / (posDif);
375 
376  double locx = distPosInTube;
377  double locy = seg1.localParameters().contains(Trk::locY) ? seg1.localParameters()[Trk::locY] : 0.;
378 
379  // create local position
380  Amg::Vector2D lpos(locx, locy);
381 
382  // segment reference surface
383  const Trk::PlaneSurface& segmentSurface = seg2.associatedSurface();
384 
385  // get global position
386  Amg::Vector3D gpos;
387  seg1.associatedSurface().localToGlobal(lpos, segDir1Min, gpos);
388  segPos = gpos;
389 
390  // calculate intersection
391  Trk::Intersection intersect = segmentSurface.straightLineIntersection(gpos, segDir1Min, false, false);
392  if (!intersect.valid) {
393  ATH_MSG_WARNING(" Intersect with surface position " << Amg::toString(gpos) << " direction: phi "
394  << segDir1Min.phi() << " theta "
395  << segDir1Min.theta());
396  goodMatch = false;
397  } else {
398  Amg::Vector3D locExSeg2 = segmentGeometry2.globalToSeg * intersect.position;
399  distPosMin2 = locExSeg2.x();
400  distPosInTube2 = segmentGeometry2.positionInsideTube(distPosMin2);
401  }
402 
403  } else {
404  goodMatch = false;
405  }
406 
407  return {distPosMin, distPosInTube, segmentGeometry.shortestChannelLength, distPosMin2,
408  distPosInTube2, segmentGeometry2.shortestChannelLength, resyMin, goodMatch,
409  segPos};
410 }
411 
412 
413 MuonSegmentInOverlapResolvingTool::SegmentGeometrySummary
414 MuonSegmentInOverlapResolvingTool::segmentGeometrySummary(const MuonSegment& seg) const
415 {
416 
417  SegmentGeometrySummary summary;
418  // loop over hits
419  Identifier tubeId1;
420  const MuonGM::MdtReadoutElement* detEl = nullptr;
421  double shortestTubeLen = 1e9;
422  Amg::Vector3D roPos;
423  Amg::Vector3D tubeCenter;
424  bool hasMdt = false;
425  for (const Trk::MeasurementBase* meas :seg.containedMeasurements()) {
426 
427  const MdtDriftCircleOnTrack* mdt = dynamic_cast<const MdtDriftCircleOnTrack*>(meas);
428  if (mdt) {
429  hasMdt = true;
430  const Identifier& id = mdt->identify();
431  int layer = m_idHelperSvc->mdtIdHelper().tubeLayer(id);
432  int tube = m_idHelperSvc->mdtIdHelper().tube(id);
433  double tubelen = mdt->prepRawData()->detectorElement()->getActiveTubeLength(layer, tube);
434  if (tubelen < shortestTubeLen) {
435  shortestTubeLen = tubelen;
436  roPos = mdt->prepRawData()->detectorElement()->ROPos(id);
437  tubeCenter = mdt->prepRawData()->detectorElement()->surface(id).center();
438  detEl = mdt->prepRawData()->detectorElement();
439  tubeId1 = id;
440  }
441  }
442  }
443  summary.detEl = detEl;
444  summary.hasMdt = hasMdt;
445  summary.segToGlobal = seg.associatedSurface().transform();
446  summary.globalToSeg = summary.segToGlobal.inverse();
447  summary.roPosInSegFrame = (summary.globalToSeg * roPos).x();
448  double distTubeCenterFromRO = (summary.globalToSeg * tubeCenter).x() - summary.roPosInSegFrame;
449  summary.hvPosInSegFrame = summary.roPosInSegFrame + 2 * distTubeCenterFromRO;
450  summary.shortestChannelLength = std::abs(2 * distTubeCenterFromRO);
451 
452  return summary;
453 }
454 
455 MuonSegmentInOverlapResolvingTool::SegmentMatchResult
456 MuonSegmentInOverlapResolvingTool::matchResult(const EventContext& ctx, const MuonSegment& seg1, const MuonSegment& seg2) const
457 {
458 
459  ATH_MSG_DEBUG(" First segment " << m_printer->print(seg1) << std::endl
460  << " Second segment " << m_printer->print(seg2));
461  SegmentMatchResult result = SegmentMatchResult();
462 
463  // calculate the phi angle that matches the two local segment angles
464  result.phiResult = bestPhiMatchAnalytic(seg1, seg2);
465 
466  // calculate the position along the tube that minimizes the position residual for both segments
467  result.segmentResult1 = bestPositionAlongTubeMatch(seg1, seg2, result.phiResult.segmentDirection1);
468  result.segmentResult2 = bestPositionAlongTubeMatch(seg2, seg1, result.phiResult.segmentDirection2);
469 
470  // calculate the average pull of the phi hits on the segments with the new parameters
471  result.averagePhiHitPullSegment1 = checkPhiHitConsistency(ctx, seg1, result.phiResult, result.segmentResult1);
472  result.averagePhiHitPullSegment2 = checkPhiHitConsistency(ctx, seg2, result.phiResult, result.segmentResult2);
473 
474  if (result.segmentResult1.goodMatch && result.segmentResult2.goodMatch) {
475 
476  // calculate vector connecting the two new segment positions
477  Amg::Vector3D difPos = result.segmentResult2.segmentPosition - result.segmentResult1.segmentPosition;
478 
479  // check if vectors are pointing in the same direction, else flip
480  if (result.phiResult.segmentDirection1.y() * difPos.y() < 0.) {
481  difPos *= -1.;
482  }
483 
484  // calculate difference in angle between phi from phi match and the two new segment positions
485  result.angularDifferencePhi = difPos.deltaPhi(result.phiResult.segmentDirection1);
486  }
487  return result;
488 }
489 
490 
491 double
492 MuonSegmentInOverlapResolvingTool::checkPhiHitConsistency(const EventContext& ctx, const Muon::MuonSegment& segment,
493  SegmentPhiMatchResult& phiMatchResult,
494  SegmentPositionMatchResult& posMatchResult) const
495 {
496 
497  // calculate average pull of phi measurements
498  unsigned int nphiMeas = 0;
499  double averagePull = 0.;
500 
501  // calculate segment AtaPlane
502  // SegmentGeometrySummary geometry = segmentGeometrySummary( segment );
503  double locx = posMatchResult.positionInTube1;
504  double locy = segment.localParameters().contains(Trk::locY) ? segment.localParameters()[Trk::locY] : 0.;
505  Trk::AtaPlane segPars(locx, locy, phiMatchResult.segmentDirection1.phi(), phiMatchResult.segmentDirection1.theta(),
506  0., segment.associatedSurface());
507 
508  // loop over hits and calculate residuals for phi hits
509  for (const Trk::MeasurementBase* meas : segment.containedMeasurements()) {
510 
511  Identifier id = m_edmHelperSvc->getIdentifier(*meas);
512  if (!id.is_valid() || !m_idHelperSvc->measuresPhi(id)) continue;
513 
514  // predict onto segment surface
515  const Trk::Surface& measSurf = meas->associatedSurface();
516 
517  // propagate station parameters to segment
518  std::unique_ptr<const Trk::TrackParameters> exPars {
519  m_propagator->propagate(
520  ctx,
521  segPars,
522  measSurf,
523  Trk::anyDirection,
524  false,
525  m_magFieldProperties)};
526  if (!exPars) {
527  ATH_MSG_WARNING(" Failed to propagate parameter to segment surface" << std::endl
528  << " pars "
529  << m_printer->print(segPars));
530  continue;
531  }
532 
533  std::optional<Trk::ResidualPull> resPull { m_pullCalculator->residualPull(meas, exPars.get(), Trk::ResidualPull::Unbiased)};
534  if (!resPull) {
535  ATH_MSG_DEBUG(" calculation of residual/pull failed !!!!! ");
536  continue;
537  }
538 
539  // sanity check
540  if (resPull->pull().size() != 1) {
541  ATH_MSG_WARNING(" ResidualPull with empty pull vector for channel " << m_idHelperSvc->toString(id));
542  continue;
543  }
544  const double pull = resPull->pull().front();
545  averagePull += pull;
546  ++nphiMeas;
547 
548  }
549  if (nphiMeas != 0) averagePull /= nphiMeas;
550  return averagePull;
551 }
552 
553 } // namespace Muon
Trk::anyDirection
@ anyDirection
Definition: PropDirection.h:22
Trk::PlaneSurface::globalToLocal
virtual bool globalToLocal(const Amg::Vector3D &glob, const Amg::Vector3D &mom, Amg::Vector2D &loc) const override final
Specified for PlaneSurface: GlobalToLocal method without dynamic memory allocation - boolean checks i...
Definition: PlaneSurface.cxx:213
Muon::MuonSegmentInOverlapResolvingTool::estimateSegmentDirection
Amg::Vector3D estimateSegmentDirection(const MuonSegment &seg1, const MuonSegment &seg2, double &phi, double &stereoangle) const
Definition: MuonSegmentInOverlapResolvingTool.cxx:69
MdtReadoutElement.h
Muon::IMuonSegmentInOverlapResolvingTool::SegmentPhiMatchResult::segmentDirection1
Amg::Vector3D segmentDirection1
new segment direction for first segment with the update phi angle
Definition: IMuonSegmentInOverlapResolvingTool.h:42
PlotCalibFromCool.norm
norm
Definition: PlotCalibFromCool.py:100
Muon::MuonSegmentInOverlapResolvingTool::checkPhiHitConsistency
double checkPhiHitConsistency(const EventContext &ctx, const Muon::MuonSegment &segment, SegmentPhiMatchResult &phiMatchResult, SegmentPositionMatchResult &matchResult) const
compare phi hits with segment parameters, return average pull of the phi hits
Definition: MuonSegmentInOverlapResolvingTool.cxx:492
Muon::MuonSegmentInOverlapResolvingTool::m_printer
PublicToolHandle< MuonEDMPrinterTool > m_printer
EDM printer tool.
Definition: MuonSegmentInOverlapResolvingTool.h:95
Trk::Intersection
Definition: Intersection.h:24
get_generator_info.result
result
Definition: get_generator_info.py:21
TrackParameters.h
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
Muon::MdtDriftCircleOnTrack::prepRawData
virtual const MdtPrepData * prepRawData() const override final
Returns the PrepRawData used to create this corrected measurement.
Definition: MdtDriftCircleOnTrack.h:257
Trk::locY
@ locY
local cartesian
Definition: ParamDefs.h:44
Amg::Vector2D
Eigen::Matrix< double, 2, 1 > Vector2D
Definition: GeoPrimitives.h:48
Muon::MuonSegmentInOverlapResolvingTool::m_propagator
ToolHandle< Trk::IPropagator > m_propagator
Definition: MuonSegmentInOverlapResolvingTool.h:100
Trk::ParametersT
Dummy class used to allow special convertors to be called for surfaces owned by a detector element.
Definition: EMErrorDetail.h:25
Muon::IMuonSegmentInOverlapResolvingTool::SegmentPhiMatchResult
Definition: IMuonSegmentInOverlapResolvingTool.h:36
initialize
void initialize()
Definition: run_EoverP.cxx:894
Muon::IMuonSegmentInOverlapResolvingTool::SegmentGeometrySummary::hvPosInSegFrame
double hvPosInSegFrame
Definition: IMuonSegmentInOverlapResolvingTool.h:62
theta
Scalar theta() const
theta method
Definition: AmgMatrixBasePlugin.h:71
Muon::IMuonSegmentInOverlapResolvingTool::SegmentGeometrySummary::shortestChannelLength
double shortestChannelLength
Definition: IMuonSegmentInOverlapResolvingTool.h:56
Trk::ResidualPull::Unbiased
@ Unbiased
RP with track state that has measurement not included.
Definition: ResidualPull.h:57
Trk::LocalParameters::contains
bool contains(ParamDefs par) const
The simple check for the clients whether the parameter is contained.
Muon::MuonSegmentInOverlapResolvingTool::MuonSegmentInOverlapResolvingTool
MuonSegmentInOverlapResolvingTool(const std::string &, const std::string &, const IInterface *)
constructor
Definition: MuonSegmentInOverlapResolvingTool.cxx:20
Muon::MuonSegmentInOverlapResolvingTool::m_edmHelperSvc
ServiceHandle< IMuonEDMHelperSvc > m_edmHelperSvc
EDM Helper tool.
Definition: MuonSegmentInOverlapResolvingTool.h:88
MdtDriftCircleOnTrack.h
Trk::Surface::center
const Amg::Vector3D & center() const
Returns the center position of the Surface.
Muon
This class provides conversion from CSC RDO data to CSC Digits.
Definition: TrackSystemController.h:49
yodamerge_tmp.scale
scale
Definition: yodamerge_tmp.py:138
Muon::IMuonSegmentInOverlapResolvingTool::SegmentGeometrySummary::positionInsideTube
double positionInsideTube(double locPosX)
checks wether the position is inside the range, if not returns the position inside the range closest ...
Definition: IMuonSegmentInOverlapResolvingTool.h:73
x
#define x
Muon::MuonSegmentInOverlapResolvingTool::m_idHelperSvc
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc
Definition: MuonSegmentInOverlapResolvingTool.h:83
Trk::LocalDirection::angleYZ
double angleYZ() const
access method for angle of local YZ projection
Definition: LocalDirection.h:106
SCT_CalibAlgs::nbins
@ nbins
Definition: SCT_CalibNumbers.h:10
Muon::MuonSegmentInOverlapResolvingTool::matchResult
SegmentMatchResult matchResult(const EventContext &ctx, const MuonSegment &seg1, const MuonSegment &seg2) const
performance match and return result
Definition: MuonSegmentInOverlapResolvingTool.cxx:456
Muon::MuonSegmentInOverlapResolvingTool::initialize
StatusCode initialize()
AlgTool initilize.
Definition: MuonSegmentInOverlapResolvingTool.cxx:28
Amg::toString
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Definition: GeoPrimitivesToStringConverter.h:40
ResidualPull.h
Muon::MuonSegmentInOverlapResolvingTool::m_magFieldProperties
Trk::MagneticFieldProperties m_magFieldProperties
magnetic field properties
Definition: MuonSegmentInOverlapResolvingTool.h:111
MuonGM::MdtReadoutElement
Definition: MuonDetDescr/MuonReadoutGeometry/MuonReadoutGeometry/MdtReadoutElement.h:50
lumiFormat.i
int i
Definition: lumiFormat.py:92
Identifier
Definition: DetectorDescription/Identifier/Identifier/Identifier.h:32
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
Muon::MuonSegmentInOverlapResolvingTool::updateSegmentDirection
Amg::Vector3D updateSegmentDirection(const MuonSegment &seg, double phi) const
recalculate the segment direction give a new angle phi, keeps the angle in the precision plane fixed
Definition: MuonSegmentInOverlapResolvingTool.cxx:42
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
TRT::Hit::layer
@ layer
Definition: HitInfo.h:79
MuonGM::MdtReadoutElement::getActiveTubeLength
double getActiveTubeLength(const int tubeLayer, const int tube) const
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition: GeoPrimitives.h:46
Trk::Segment::containedMeasurements
const std::vector< const Trk::MeasurementBase * > & containedMeasurements() const
returns the vector of Trk::MeasurementBase objects
Definition: TrkEvent/TrkSegment/TrkSegment/Segment.h:166
Muon::MuonSegmentInOverlapResolvingTool::bestPositionAlongTubeMatch
SegmentPositionMatchResult bestPositionAlongTubeMatch(const MuonSegment &seg1, const MuonSegment &seg2, const Amg::Vector3D &segDir1Min) const
calculate the position along a tube direction of the first segment that results in the smallest posit...
Definition: MuonSegmentInOverlapResolvingTool.cxx:305
Muon::IMuonSegmentInOverlapResolvingTool::SegmentGeometrySummary
Definition: IMuonSegmentInOverlapResolvingTool.h:53
python.StandardJetMods.pull
pull
Definition: StandardJetMods.py:264
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Trk::PlaneSurface::globalToLocalDirection
void globalToLocalDirection(const Amg::Vector3D &glodir, Trk::LocalDirection &locdir) const
This method transforms the global direction to a local direction wrt the plane.
Definition: PlaneSurface.cxx:260
Muon::IMuonSegmentInOverlapResolvingTool::SegmentGeometrySummary::roPosInSegFrame
double roPosInSegFrame
Definition: IMuonSegmentInOverlapResolvingTool.h:61
Trk::LocalDirection
represents the three-dimensional global direction with respect to a planar surface frame.
Definition: LocalDirection.h:81
Muon::IMuonSegmentInOverlapResolvingTool::SegmentPositionMatchResult::positionInTube1
double positionInTube1
Definition: IMuonSegmentInOverlapResolvingTool.h:112
Trk::MeasurementBase
Definition: MeasurementBase.h:58
Trk::NoField
@ NoField
Field is set to 0., 0., 0.,.
Definition: MagneticFieldMode.h:18
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
Muon::IMuonSegmentInOverlapResolvingTool::SegmentPositionMatchResult
Definition: IMuonSegmentInOverlapResolvingTool.h:97
id
SG::auxid_t id
Definition: Control/AthContainers/Root/debug.cxx:194
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
Muon::MdtDriftCircleOnTrack
This class represents the corrected MDT measurements, where the corrections include the effects of wi...
Definition: MdtDriftCircleOnTrack.h:37
Muon::MuonSegmentInOverlapResolvingTool::bestPhiMatch
SegmentPhiMatchResult bestPhiMatch(const MuonSegment &seg1, const MuonSegment &seg2) const
calculate the angle phi for which the angular residual in the precision plane of the second segment i...
Definition: MuonSegmentInOverlapResolvingTool.cxx:275
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Trk::PlaneSurface::straightLineIntersection
virtual Intersection straightLineIntersection(const Amg::Vector3D &pos, const Amg::Vector3D &dir, bool forceDir, Trk::BoundaryCheck bchk) const override final
fast straight line intersection schema - standard: provides closest intersection and (signed) path le...
Definition: PlaneSurface.cxx:223
python.copyTCTOutput.locDir
locDir
Definition: copyTCTOutput.py:113
MuonGM::MdtReadoutElement::ROPos
Amg::Vector3D ROPos(const int tubelayer, const int tube) const
Definition: MuonDetDescr/MuonReadoutGeometry/src/MdtReadoutElement.cxx:288
Trk::MeasurementBase::localParameters
const LocalParameters & localParameters() const
Interface method to get the LocalParameters.
Definition: MeasurementBase.h:132
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
a
TList * a
Definition: liststreamerinfos.cxx:10
Muon::MuonSegmentInOverlapResolvingTool::m_pullCalculator
ToolHandle< Trk::IResidualPullCalculator > m_pullCalculator
Definition: MuonSegmentInOverlapResolvingTool.h:105
Amg::intersect
std::optional< double > intersect(const AmgVector(N)&posA, const AmgVector(N)&dirA, const AmgVector(N)&posB, const AmgVector(N)&dirB)
Calculates the closest approach of two lines.
Definition: GeoPrimitivesHelpers.h:302
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Trk::PlaneSurface
Definition: PlaneSurface.h:64
GeoPrimitivesHelpers.h
MuonSegmentInOverlapResolvingTool.h
Muon::IMuonSegmentInOverlapResolvingTool::SegmentMatchResult
Definition: IMuonSegmentInOverlapResolvingTool.h:133
MuonGM::MdtReadoutElement::surface
virtual const Trk::Surface & surface() const override final
Return surface associated with this detector element.
Definition: MuonDetDescr/MuonReadoutGeometry/src/MdtReadoutElement.cxx:875
Trk::RIO_OnTrack::identify
virtual Identifier identify() const final
return the identifier -extends MeasurementBase
Definition: RIO_OnTrack.h:155
GeoPrimitivesToStringConverter.h
LArCellBinning.step
step
Definition: LArCellBinning.py:158
MuonSegment.h
Muon::MuonSegment::globalPosition
virtual const Amg::Vector3D & globalPosition() const override final
global position
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:157
LocalDirection.h
Muon::MuonSegmentInOverlapResolvingTool::bestPhiMatchAnalytic
SegmentPhiMatchResult bestPhiMatchAnalytic(const MuonSegment &seg1, const MuonSegment &seg2) const
calculate the angle phi for which the angular residual in the precision plane of the second segment i...
Definition: MuonSegmentInOverlapResolvingTool.cxx:194
Muon::MuonSegment
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:45
Muon::IMuonSegmentInOverlapResolvingTool::SegmentGeometrySummary::globalToSeg
Amg::Transform3D globalToSeg
Definition: IMuonSegmentInOverlapResolvingTool.h:59
AthAlgTool
Definition: AthAlgTool.h:26
Trk::PlaneSurface::localToGlobal
virtual void localToGlobal(const Amg::Vector2D &locp, const Amg::Vector3D &mom, Amg::Vector3D &glob) const override final
Specified for PlaneSurface: LocalToGlobal method without dynamic memory allocation.
Definition: PlaneSurface.cxx:204
Muon::MuonSegmentInOverlapResolvingTool::segmentGeometrySummary
SegmentGeometrySummary segmentGeometrySummary(const MuonSegment &seg) const
calculate geometrical information for a segment
Definition: MuonSegmentInOverlapResolvingTool.cxx:414
Trk::Surface
Definition: Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:75
Trk::Surface::transform
const Amg::Transform3D & transform() const
Returns HepGeom::Transform3D by reference.
Intersection.h
Muon::MdtPrepData::detectorElement
virtual const MuonGM::MdtReadoutElement * detectorElement() const override
Returns the detector element corresponding to this PRD.
Definition: MdtPrepData.h:156
calibdata.tube
tube
Definition: calibdata.py:31
Muon::MuonSegment::associatedSurface
virtual const Trk::PlaneSurface & associatedSurface() const override final
returns the surface for the local to global transformation
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:175
Muon::MuonSegment::globalDirection
const Amg::Vector3D & globalDirection() const
global direction
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:163
Amg::setThetaPhi
void setThetaPhi(Amg::Vector3D &v, double theta, double phi)
sets the theta and phi angle of a vector without changing the magnitude
Definition: GeoPrimitivesHelpers.h:70
NSWL1::PadTriggerAdapter::segment
Muon::NSW_PadTriggerSegment segment(const NSWL1::PadTrigger &data)
Definition: PadTriggerAdapter.cxx:5
SCT_Monitoring::summary
@ summary
Definition: SCT_MonitoringNumbers.h:65
Generated on Thu Jul 4 2024 21:23:06 for ATLAS Offline Software by   doxygen 1.8.18