ATLAS Offline Software
SegmentFittingAlg.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "SegmentFittingAlg.h"
6 
7 #include <TrkEventPrimitives/ParticleHypothesis.h>
8 #include <GeoPrimitives/GeoPrimitivesHelpers.h>
9 
10 #include <MuonPatternHelpers/SegmentFitHelperFunctions.h>
11 #include <MuonPatternHelpers/MdtSegmentSeedGenerator.h>
12 #include <MuonPatternHelpers/CalibSegmentChi2Minimizer.h>
13 #include <MuonPatternHelpers/MdtSegmentFitter.h>
14 
15 #include <MuonSpacePoint/SpacePointPerLayerSorter.h>
16 #include <MuonSpacePoint/UtilFunctions.h>
17 
18 #include <xAODMuonPrepData/RpcMeasurement.h>
19 #include <xAODMuonPrepData/TgcStrip.h>
20 #include "xAODMuonSimHit/MuonSimHitContainer.h"
21 
22 #include <GaudiKernel/PhysicalConstants.h>
23 #include <Minuit2/Minuit2Minimizer.h>
24 #include <Math/Minimizer.h>
25 
26 #include <MuonVisualizationHelpersR4/VisualizationHelpers.h>
27 
28 namespace MuonR4 {
29  using namespace SegmentFit;
30  using namespace MuonValR4;
31 
32  constexpr double inv_c = 1./ Gaudi::Units::c_light;
33 
34 
35  MuonR4::SegmentFitResult::HitVec copy(const MuonR4::SegmentFitResult::HitVec& hits) {
36  MuonR4::SegmentFitResult::HitVec copied{};
37  copied.reserve(hits.size());
38  std::ranges::transform(hits, std::back_inserter(copied),
39  [](const auto& hit) { return std::make_unique<MuonR4::CalibratedSpacePoint>(*hit);});
40  return copied;
41  }
42  MuonR4::SegmentFitResult copy(const MuonR4::SegmentFitResult& toCopy) {
43  MuonR4::SegmentFitResult toRet{};
44  toRet.segmentPars = toCopy.segmentPars;
45  toRet.calibMeasurements = copy(toCopy.calibMeasurements);
46  toRet.chi2 = toCopy.chi2;
47  toRet.nDoF = toCopy.nDoF;
48  toRet.timeFit = toCopy.timeFit;
49  return toRet;
50  }
51  bool removeBeamSpot(MuonR4::SegmentFitResult::HitVec& hits){
52  std::size_t before = hits.size();
53  hits.erase(std::remove_if(hits.begin(), hits.end(),
54  [](const auto& a){
55  return a->type() == xAOD::UncalibMeasType::Other;
56  }), hits.end());
57  return before != hits.size();
58  }
59 
60  using PrimitiveVec = MuonValR4::IPatternVisualizationTool::PrimitiveVec;
61 
62  SegmentFittingAlg::~SegmentFittingAlg() = default;
63  StatusCode SegmentFittingAlg::initialize() {
64  ATH_CHECK(m_geoCtxKey.initialize());
65  ATH_CHECK(m_seedKey.initialize());
66  ATH_CHECK(m_outSegments.initialize());
67  ATH_CHECK(m_calibTool.retrieve());
68  ATH_CHECK(m_idHelperSvc.retrieve());
69  ATH_CHECK(m_visionTool.retrieve(EnableTool{!m_visionTool.empty()}));
70  SegmentAmbiSolver::Config cfg{};
71  m_ambiSolver = std::make_unique<SegmentAmbiSolver>(name(), std::move(cfg));
72  return StatusCode::SUCCESS;
73  }
74  StatusCode SegmentFittingAlg::execute(const EventContext& ctx) const {
75 
76  const ActsGeometryContext* gctx{nullptr};
77  ATH_CHECK(retrieveContainer(ctx, m_geoCtxKey, gctx));
78  const SegmentSeedContainer* segmentSeeds=nullptr;
79  ATH_CHECK(retrieveContainer(ctx, m_seedKey, segmentSeeds));
80 
81  SG::WriteHandle writeSegments{m_outSegments, ctx};
82  ATH_CHECK(writeSegments.record(std::make_unique<SegmentContainer>()));
83  std::vector<std::unique_ptr<Segment>> allSegments{};
84  for (const SegmentSeed* seed : *segmentSeeds) {
85  std::vector<std::unique_ptr<Segment>> segments = fitSegmentSeed(ctx, *gctx, seed);
86  if (m_visionTool.isEnabled() && segments.size() > 1) {
87  auto drawFinalReco = [this, &segments, &gctx, &ctx,&seed](const std::string& nameTag) {
88  PrimitiveVec segmentLines{};
89  double yLegend{0.85};
90  segmentLines.push_back(drawLabel(std::format("# segments: {:d}", segments.size()), 0.12, yLegend, 14));
91  yLegend-=0.04;
92  for (const std::unique_ptr<Segment>& seg : segments) {
93  const Parameters pars = localSegmentPars(*gctx, *seg);
94  const auto [locPos, locDir] = makeLine(pars);
95 
96  segmentLines.emplace_back(drawLine(pars, -Gaudi::Units::m, Gaudi::Units::m, kRed));
97  std::stringstream signStream{};
98  signStream<<"#chi^{2}/nDoF: "<<std::format("{:.2f}", seg->chi2() / seg->nDoF())<<"("<<seg->nDoF()<<"), ";
99  signStream<<"y_{0}="<<std::format("{:.2f}",pars[toInt(ParamDefs::y0)])<<", ";
100  signStream<<std::format("#theta={:.2f}", pars[toInt(ParamDefs::theta)]/ Gaudi::Units::deg )<<", ";
101  for (const Segment::MeasType& m : seg->measurements()) {
102  if (m->type() == xAOD::UncalibMeasType::MdtDriftCircleType && m->fitState() == CalibratedSpacePoint::State::Valid) {
103  signStream<<(SegmentFitHelpers::driftSign(locPos, locDir, *m, msgStream()) == -1 ? "L" : "R");
104  }
105  }
106  segmentLines.push_back(drawLabel(signStream.str(), 0.12, yLegend, 12));
107  yLegend-=0.03;
108  }
109 
110  m_visionTool->visualizeBucket(ctx, *seed->parentBucket(), nameTag, std::move(segmentLines));
111  };
112  drawFinalReco("all segments");
113  const unsigned int nBeforeAmbi = segments.size();
114  segments = m_ambiSolver->resolveAmbiguity(*gctx, std::move(segments));
115  if (nBeforeAmbi != segments.size()) {
116  drawFinalReco("post ambiguity");
117  }
118  }
119  allSegments.insert(allSegments.end(), std::make_move_iterator(segments.begin()),
120  std::make_move_iterator(segments.end()));
121  }
122  resolveAmbiguities(*gctx, allSegments);
123  writeSegments->insert(writeSegments->end(),
124  std::make_move_iterator(allSegments.begin()),
125  std::make_move_iterator(allSegments.end()));
126  ATH_MSG_VERBOSE("Found in total "<<writeSegments->size()<<" segments. ");
127  return StatusCode::SUCCESS;
128  }
129 
130  template <class ContainerType>
131  StatusCode SegmentFittingAlg::retrieveContainer(const EventContext& ctx,
132  const SG::ReadHandleKey<ContainerType>& key,
133  const ContainerType*& contToPush) const {
134  contToPush = nullptr;
135  if (key.empty()) {
136  ATH_MSG_VERBOSE("No key has been parsed for object "<< typeid(ContainerType).name());
137  return StatusCode::SUCCESS;
138  }
139  SG::ReadHandle readHandle{key, ctx};
140  ATH_CHECK(readHandle.isPresent());
141  contToPush = readHandle.cptr();
142  return StatusCode::SUCCESS;
143  }
144 
145  SegmentFitResult SegmentFittingAlg::fitSegmentHits(const EventContext& ctx,
146  const ActsGeometryContext& gctx,
147  const Parameters& startPars,
148  SegmentFitResult::HitVec&& calibHits) const {
149  SegmentFitResult data{};
150  if (calibHits.empty()) {
151  ATH_MSG_WARNING("No hits, no segment");
152  return data;
153  }
154 
156  if (m_doBeamspotConstraint) {
157  unsigned int numPhi = std::accumulate(calibHits.begin(), calibHits.end(),0,
158  [](unsigned int n, const HitVec::value_type& hit){
159  return n + (hit->fitState() == CalibratedSpacePoint::State::Valid &&
160  hit->measuresPhi());
161  });
162  if (numPhi) {
163  const Amg::Transform3D globToLoc{calibHits[0]->spacePoint()->msSector()->globalToLocalTrans(gctx)};
164  Amg::Vector3D beamSpot{globToLoc.translation()};
165  AmgSymMatrix(3) covariance{AmgSymMatrix(3)::Identity()};
167  covariance(0,0) = std::pow(m_beamSpotR, 2);
168  covariance(1,1) = std::pow(m_beamSpotR, 2);
169  covariance(2,2) = std::pow(m_beamSpotL, 2);
170  AmgSymMatrix(3) jacobian = globToLoc.linear();
171  covariance = jacobian * covariance * jacobian.transpose();
172  AmgSymMatrix(2) beamSpotCov{covariance.block<2,2>(0,0)};
173  auto beamSpotSP = std::make_unique<CalibratedSpacePoint>(nullptr, std::move(beamSpot), Amg::Vector3D::Zero());
174  beamSpotSP->setCovariance<2>(std::move(beamSpotCov));
175  ATH_MSG_VERBOSE("Beam spot constraint "<<Amg::toString(beamSpotSP->positionInChamber())<<", "<<toString(beamSpotSP->covariance()));
176  calibHits.push_back(std::move(beamSpotSP));
177  }
178  }
179 
180  const Amg::Transform3D& locToGlob{calibHits[0]->spacePoint()->msSector()->localToGlobalTrans(gctx)};
181 
182  if (!m_useMinuit) {
183  MdtSegmentFitter::Config fitCfg{};
184  fitCfg.calibrator = m_calibTool.get();
185  fitCfg.doTimeFit = m_doT0Fit;
186 
187  MdtSegmentFitter fitter{name(), std::move(fitCfg)};
188  return fitter.fitSegment(ctx, std::move(calibHits), startPars, locToGlob);
189  }
190 
191  data.segmentPars = startPars;
192 
193  CalibSegmentChi2Minimizer c2f{name(), ctx, locToGlob, copy(calibHits), m_calibTool.get(), m_doT0Fit};
194  data.hasPhi = c2f.hasPhiMeas();
195  data.timeFit = c2f.doTimeFit();
196  data.nDoF = c2f.nDoF();
197  if (data.nDoF <= 0) {
198  ATH_MSG_DEBUG("Reject fit due to 0 degrees of freedom");
199  return data;
200  }
201  ROOT::Minuit2::Minuit2Minimizer minimizer((name() + std::to_string(ctx.eventID().event_number())).c_str());
203  minimizer.SetMaxFunctionCalls(100000);
204  minimizer.SetTolerance(0.0001);
205  minimizer.SetPrintLevel(-1);
206  minimizer.SetStrategy(1);
207 
208  minimizer.SetVariable(toInt(ParamDefs::y0), "y0", startPars[toInt(ParamDefs::y0)], 1.e-5);
209  minimizer.SetVariable(toInt(ParamDefs::theta), "theta", startPars[toInt(ParamDefs::theta)], 1.e-5);
210  minimizer.SetVariableLimits(toInt(ParamDefs::y0),
211  startPars[toInt(ParamDefs::y0)] - 60. *Gaudi::Units::cm,
212  startPars[toInt(ParamDefs::y0)] + 60. *Gaudi::Units::cm);
213  minimizer.SetVariableLimits(toInt(ParamDefs::theta),
214  startPars[toInt(ParamDefs::theta)] - 0.6,
215  startPars[toInt(ParamDefs::theta)] + 0.6);
216 
217  if (data.hasPhi) {
218  minimizer.SetVariable(toInt(ParamDefs::x0), "x0", startPars[toInt(ParamDefs::x0)], 1.e-5);
219  minimizer.SetVariable(toInt(ParamDefs::phi), "phi", startPars[toInt(ParamDefs::phi)], 1.e-5);
220  minimizer.SetVariableLimits(toInt(ParamDefs::x0),
221  startPars[toInt(ParamDefs::x0)] - 600,
222  startPars[toInt(ParamDefs::x0)] + 600);
223  minimizer.SetVariableLimits(toInt(ParamDefs::phi),
224  startPars[toInt(ParamDefs::phi)] - 0.6,
225  startPars[toInt(ParamDefs::phi)] + 0.6);
226  } else {
227  minimizer.SetFixedVariable(toInt(ParamDefs::x0), "x0", 0.);
228  minimizer.SetFixedVariable(toInt(ParamDefs::phi), "phi", 90.*Gaudi::Units::deg);
229  }
231  if (data.timeFit) {
232  minimizer.SetVariable(toInt(ParamDefs::time), "t0", startPars[toInt(ParamDefs::time)] , 1.);
233  minimizer.SetVariableLimits(toInt(ParamDefs::time),
234  startPars[toInt(ParamDefs::time)] - 50,
235  startPars[toInt(ParamDefs::time)] + 50);
236  } else{
237  minimizer.SetFixedVariable(toInt(ParamDefs::time), "t0", 0.);
238  }
239  minimizer.SetFunction(c2f);
241  if (!minimizer.Minimize() || !minimizer.Hesse()) {
242  data.calibMeasurements = std::move(calibHits);
243  } else {
244  const double* xs = minimizer.X();
245  const double* errs = minimizer.Errors();
246 
247  for (unsigned int p = 0; p < toInt(ParamDefs::nPars); ++p) {
248  data.segmentPars[p] = xs[p];
249  data.segmentParErrs(p,p) = errs[p];
250  }
251  std::optional<double> timeOfArrival{std::nullopt};
252  const auto [locPos, locDir] = data.makeLine();
253  if (data.timeFit) {
254  timeOfArrival = std::make_optional<double>((locToGlob*locPos).mag() * inv_c);
255  }
256  data.nIter = minimizer.NCalls();
257  data.calibMeasurements = c2f.release(xs);
258 
259  auto [chiPerMeas, finalChi2] = SegmentFitHelpers::postFitChi2PerMas(data.segmentPars, timeOfArrival,
260  data.calibMeasurements, msgStream());
261  data.chi2PerMeasurement = std::move(chiPerMeas);
262  data.chi2 = finalChi2;
263  data.converged = true;
264  }
265  return data;
266  }
267  std::vector<std::unique_ptr<Segment>>
268  SegmentFittingAlg::fitSegmentSeed(const EventContext& ctx,
269  const ActsGeometryContext& gctx,
270  const SegmentSeed* patternSeed) const {
271 
272  const Amg::Transform3D& locToGlob{patternSeed->msSector()->localToGlobalTrans(gctx)};
273  std::vector<std::unique_ptr<Segment>> segments{};
274 
275  MdtSegmentSeedGenerator::Config genCfg{};
276  genCfg.hitPullCut = m_seedHitChi2;
277  genCfg.recalibSeedCircles = m_recalibSeed;
278  genCfg.fastSeedFit = m_refineSeed;
279  genCfg.calibrator = m_calibTool.get();
280 
281 
283  genCfg.busyLayerLimit = 2 + 2*(patternSeed->parameters()[toInt(ParamDefs::theta)] > 50 * Gaudi::Units::deg);
285  if (m_visionTool.isEnabled()) {
286  PrimitiveVec seedLines{};
287  MdtSegmentSeedGenerator drawMe{name(), patternSeed, genCfg};
288  while(auto s = drawMe.nextSeed(ctx)) {
289  seedLines.push_back(drawLine(s->parameters, -Gaudi::Units::m, Gaudi::Units::m, kViolet));
290  }
291  seedLines.push_back(drawLabel(std::format("possible seeds: {:d}", drawMe.numGenerated()), 0.15, 0.85, 14));
292  m_visionTool->visualizeSeed(ctx, *patternSeed, "pattern", std::move(seedLines));
293  }
294 
295  MdtSegmentSeedGenerator seedGen{name(), patternSeed, std::move(genCfg)};
296  while (auto seed = seedGen.nextSeed(ctx)) {
297  SegmentFitResult data{};
298  data.segmentPars = seed->parameters;
299  data.calibMeasurements = std::move(seed->measurements);
301  if (m_visionTool.isEnabled()) {
302  auto seedCopy = convertToSegment(locToGlob, patternSeed, copy(data));
303  m_visionTool->visualizeSegment(ctx, *seedCopy, std::format("Pre fit {:d}", seedGen.numGenerated()));
304  }
305  data = fitSegmentHits(ctx, gctx, seed->parameters, std::move(data.calibMeasurements));
306  data.nIter += seed->nIter;
307  if (m_visionTool.isEnabled() && data.converged) {
308  auto seedCopy = convertToSegment(locToGlob, patternSeed, copy(data));
309  m_visionTool->visualizeSegment(ctx, *seedCopy, std::format("Intermediate fit {:d}", seedGen.numGenerated()));
310  }
311  if (!removeOutliers(ctx, gctx, *patternSeed, data)) {
312  continue;
313  }
314  if (!plugHoles(ctx, gctx, *patternSeed, data)) {
315  continue;
316  }
317  if (m_visionTool.isEnabled()) {
318  auto seedCopy = convertToSegment(locToGlob, patternSeed, copy(data));
319  m_visionTool->visualizeSegment(ctx, *seedCopy, std::format("Final fit {:d}", seedGen.numGenerated()));
320  }
321  segments.push_back(convertToSegment(locToGlob, patternSeed, std::move(data)));
322  }
323  return segments;
324  }
325  std::unique_ptr<Segment> SegmentFittingAlg::convertToSegment(const Amg::Transform3D& locToGlob,
326  const SegmentSeed* patternSeed,
327  SegmentFitResult&& data) {
328  const auto [locPos, locDir] = data.makeLine();
329  Amg::Vector3D globPos = locToGlob * locPos;
330  Amg::Vector3D globDir = locToGlob.linear()* locDir;
331 
332 
333  auto finalSeg = std::make_unique<Segment>(std::move(globPos), std::move(globDir),
334  patternSeed, std::move(data.calibMeasurements),
335  data.chi2, data.nDoF);
336  finalSeg->setCallsToConverge(data.nIter);
337  finalSeg->setChi2PerMeasurement(std::move(data.chi2PerMeasurement));
338  finalSeg->setParUncertainties(std::move(data.segmentParErrs));
339  if (data.timeFit) {
340  finalSeg->setSegmentT0(data.segmentPars[toInt(ParamDefs::time)]);
341 
342  }
343 
344 
345  return finalSeg;
346  }
347 
348  bool SegmentFittingAlg::removeOutliers(const EventContext& ctx,
349  const ActsGeometryContext& gctx,
350  const SegmentSeed& seed,
351  SegmentFitResult& data) const {
352 
354  if (data.nDoF<=0 || data.calibMeasurements.empty()) {
355  ATH_MSG_VERBOSE("No degree of freedom available. What shall be removed?!. nDoF: "
356  <<data.nDoF<<", n-meas: "<<data.calibMeasurements);
357  return false;
358  }
359 
360  const auto [segPos, segDir] = data.makeLine();
361 
362  if (data.converged && data.chi2 / data.nDoF < m_outlierRemovalCut) {
363  ATH_MSG_VERBOSE("The segment "<<Amg::toString(segPos)<<" + "<<Amg::toString(segDir)
364  <<" is already of good quality "<<data.chi2 / std::max(data.nDoF, 1)
365  <<". Don't remove outliers");
366  return true;
367  }
368  ATH_MSG_VERBOSE("Segment "<<toString(data.segmentPars)<<" is of badish quality.");
371  if (m_doBeamspotConstraint && removeBeamSpot(data.calibMeasurements)) {
372  data.nDoF-=2;
373  data.nPhiMeas-=1;
374  }
375 
377  std::sort(data.calibMeasurements.begin(), data.calibMeasurements.end(),
378  [&, this](const HitVec::value_type& a, const HitVec::value_type& b){
379  return SegmentFitHelpers::chiSqTerm(segPos, segDir, data.segmentPars[toInt(ParamDefs::time)], std::nullopt, *a, msgStream()) <
380  SegmentFitHelpers::chiSqTerm(segPos, segDir, data.segmentPars[toInt(ParamDefs::time)], std::nullopt, *b, msgStream());
381  });
382 
384  data.calibMeasurements.back()->setFitState(CalibratedSpacePoint::State::Outlier);
385  data.nDoF -= data.calibMeasurements.back()->measuresEta();
386  data.nDoF -= data.calibMeasurements.back()->measuresPhi();
387  if (m_doT0Fit && data.calibMeasurements.back()->measuresTime() &&
388  data.calibMeasurements.back()->type() != xAOD::UncalibMeasType::MdtDriftCircleType) {
389  --data.nDoF;
390  }
392  std::vector<HoughHitType> uncalib{};
393  for (const HitVec::value_type& calib : data.calibMeasurements) {
394  uncalib.push_back(calib->spacePoint());
395  }
396  SegmentFitResult newAttempt = fitSegmentHits(ctx, gctx, data.segmentPars, std::move(data.calibMeasurements));
397  if (newAttempt.converged) {
398  newAttempt.nIter+=data.nIter;
399  data = std::move(newAttempt);
400  if (m_visionTool.isEnabled()) {
401  auto seedCopy = convertToSegment(seed.msSector()->localToGlobalTrans(gctx), &seed, copy(data));
402  m_visionTool->visualizeSegment(ctx, *seedCopy, "Bad fit recovery");
403  }
404  } else {
405  data.calibMeasurements = std::move(newAttempt.calibMeasurements);
406  }
407  return removeOutliers(ctx, gctx, seed, data);
408  }
409 
410  bool SegmentFittingAlg::plugHoles(const EventContext& ctx,
411  const ActsGeometryContext& gctx,
412  const SegmentSeed& seed,
413  SegmentFitResult& beforeRecov) const {
415  ATH_MSG_VERBOSE("plugHoles() -- segment "<<toString(beforeRecov.segmentPars)
416  <<", chi2: "<<beforeRecov.chi2<<", nDoF: "<<beforeRecov.nDoF<<", "
417  <<beforeRecov.chi2 /std::max(beforeRecov.nDoF, 1) );
419  std::unordered_set<const SpacePoint*> usedSpacePoint{};
420  for (const HitVec::value_type& hit : beforeRecov.calibMeasurements) {
421  usedSpacePoint.insert(hit->spacePoint());
422  }
423 
424  HitVec candidateHits{};
425  SpacePointPerLayerSorter hitLayers{*seed.parentBucket()};
426  bool hasCandidate{false};
427  const auto [locPos, locDir] = beforeRecov.makeLine();
428  for (const std::vector<HoughHitType>& mdtLayer : hitLayers.mdtHits()) {
429  for (const SpacePoint* mdtHit: mdtLayer) {
431  if (usedSpacePoint.count(mdtHit)) {
432  continue;
433  }
434  const double dist = Amg::lineDistance(locPos, locDir, mdtHit->positionInChamber(), mdtHit->directionInChamber());
435  const auto* dc = static_cast<const xAOD::MdtDriftCircle*>(mdtHit->primaryMeasurement());
436  if (dist >= dc->readoutElement()->innerTubeRadius()) {
437  continue;
438  }
439  HitVec::value_type calibHit{m_calibTool->calibrate(ctx, mdtHit, locPos, locDir, beforeRecov.segmentPars[toInt(ParamDefs::time)])};
440  const double pull = std::sqrt(SegmentFitHelpers::chiSqTermMdt(locPos, locDir, *calibHit, msgStream()));
441  ATH_MSG_VERBOSE(__func__<<"() :"<<__LINE__<<" Candidate hit for recovery "<<m_idHelperSvc->toString(mdtHit->identify())<<", chi2: "<<pull);
442  if (pull <= m_recoveryPull) {
443  hasCandidate |= calibHit->fitState() == CalibratedSpacePoint::State::Valid;
444  candidateHits.push_back(std::move(calibHit));
445  } else {
446  calibHit->setFitState(CalibratedSpacePoint::State::Outlier);
447  candidateHits.push_back(std::move(calibHit));
448  }
449  }
450  }
451  if (!hasCandidate) {
452  ATH_MSG_VERBOSE("No space point candidates for recovery were found");
453  beforeRecov.calibMeasurements.insert(beforeRecov.calibMeasurements.end(),
454  std::make_move_iterator(candidateHits.begin()),
455  std::make_move_iterator(candidateHits.end()));
456  eraseWrongHits(gctx, beforeRecov);
457  return beforeRecov.nDoF > 0;
458  }
459 
460  HitVec copied = copy(beforeRecov.calibMeasurements), copiedCandidates = copy(candidateHits);
462  if (m_doBeamspotConstraint) {
463  removeBeamSpot(copied);
464  }
465 
466  candidateHits.insert(candidateHits.end(), std::make_move_iterator(copied.begin()), std::make_move_iterator(copied.end()));
467 
468  SegmentFitResult recovered = fitSegmentHits(ctx, gctx, beforeRecov.segmentPars, std::move(candidateHits));
469  if (!recovered.converged) {
470  return false;
471  }
473  if (recovered.nDoF + recovered.timeFit <= beforeRecov.nDoF + beforeRecov.timeFit) {
474  for (HitVec::value_type& hit : copiedCandidates) {
475  hit->setFitState(CalibratedSpacePoint::State::Outlier);
476  beforeRecov.calibMeasurements.push_back(std::move(hit));
477  }
478  eraseWrongHits(gctx, beforeRecov);
479  return true;
480  }
481  ATH_MSG_VERBOSE("Chi2, nDOF before:"<<beforeRecov.chi2<<", "<<beforeRecov.nDoF<<" after recovery: "<<recovered.chi2<<", "<<recovered.nDoF);
482  double redChi2 = recovered.chi2 / std::max(recovered.nDoF,1);
484  if (redChi2 < m_outlierRemovalCut || (beforeRecov.nDoF == 0) || redChi2 < beforeRecov.chi2 / beforeRecov.nDoF) {
485  ATH_MSG_VERBOSE("Accept segment with recovered "<<(recovered.nDoF + recovered.timeFit) - (beforeRecov.nDoF + beforeRecov.timeFit)<<" hits.");
486  recovered.nIter += beforeRecov.nIter;
487  beforeRecov = std::move(recovered);
489  while (true) {
490  bool runAnotherTrial = false;
491  copied = copy(beforeRecov.calibMeasurements);
492  if (m_doBeamspotConstraint) {
493  removeBeamSpot(copied);
494  }
495  for (unsigned int m = 0; m < beforeRecov.calibMeasurements.size(); ++m) {
496  if (beforeRecov.calibMeasurements[m]->fitState() == CalibratedSpacePoint::State::Outlier &&
497  std::sqrt(beforeRecov.chi2PerMeasurement[m]) < m_recoveryPull) {
498  copied[m]->setFitState(CalibratedSpacePoint::State::Valid);
499  runAnotherTrial = true;
500  }
501  }
502  if (!runAnotherTrial) {
503  break;
504  }
505  recovered = fitSegmentHits(ctx, gctx, beforeRecov.segmentPars, std::move(copied));
506  if (!recovered.converged) {
507  break;
508  }
509  if (recovered.nDoF + recovered.timeFit <= beforeRecov.nDoF + beforeRecov.timeFit) {
510  break;
511  }
512  redChi2 = recovered.chi2 / std::max(recovered.nDoF, 1);
513  if (redChi2 < m_outlierRemovalCut || redChi2 < beforeRecov.chi2 / beforeRecov.nDoF) {
514  recovered.nIter += beforeRecov.nIter;
515  beforeRecov = std::move(recovered);
516  } else {
517  break;
518  }
519  }
521  eraseWrongHits(gctx, beforeRecov);
522  } else{
523  for (HitVec::value_type& hit : copiedCandidates) {
524  hit->setFitState(CalibratedSpacePoint::State::Outlier);
525  beforeRecov.calibMeasurements.push_back(std::move(hit));
526  }
527  }
528  return true;
529  }
530  void SegmentFittingAlg::eraseWrongHits(const ActsGeometryContext& gctx, SegmentFitResult& candidate) const {
531  auto [segPos, segDir] = makeLine(candidate.segmentPars);
532  candidate.calibMeasurements.erase(std::remove_if(candidate.calibMeasurements.begin(), candidate.calibMeasurements.end(),
533  [&segPos, &segDir](const HitVec::value_type& hit){
534  if (hit->fitState() != CalibratedSpacePoint::State::Outlier) {
535  return false;
536  }
538  if (hit->type() == xAOD::UncalibMeasType::MdtDriftCircleType) {
539  const double dist = Amg::lineDistance(segPos, segDir, hit->positionInChamber(), hit->directionInChamber());
540  const auto* dc = static_cast<const xAOD::MdtDriftCircle*>(hit->spacePoint()->primaryMeasurement());
541  return dist >= dc->readoutElement()->innerTubeRadius();
542  }
543  return true;
544  }), candidate.calibMeasurements.end());
545 
546  std::ranges::sort(candidate.calibMeasurements, [](const Segment::MeasType& a, const Segment::MeasType& b){
547  return a->positionInChamber().z() < b->positionInChamber().z();
548  });
549  const MuonGMR4::SpectrometerSector* chamber{nullptr};
550  for (const auto& hit : candidate.calibMeasurements) {
551  if (hit->type() != xAOD::UncalibMeasType::Other){
552  chamber = hit->spacePoint()->msSector();
553  break;
554  }
555  }
556  std::optional<double> timeOfFlight = candidate.timeFit ? std::make_optional<double>((chamber->localToGlobalTrans(gctx)*segPos).mag() * inv_c)
557  : std::nullopt;
558  auto [updatedMeasChi2, updateChi2] = SegmentFitHelpers::postFitChi2PerMas(candidate.segmentPars, timeOfFlight,
559  candidate.calibMeasurements, msgStream());
560  ATH_MSG_VERBOSE("The measurements before "<<candidate.chi2PerMeasurement<<", after: "<<updatedMeasChi2<<", chi2: "<<updateChi2);
561  candidate.chi2PerMeasurement = std::move(updatedMeasChi2);
562  }
563  void SegmentFittingAlg::resolveAmbiguities(const ActsGeometryContext& gctx,
564  std::vector<std::unique_ptr<Segment>>& segmentCandidates) const {
565  using SegmentVec = std::vector<std::unique_ptr<Segment>>;
566  ATH_MSG_VERBOSE("Resolve ambiguities amongst "<<segmentCandidates.size()<<" segment candidates. ");
567  std::unordered_map<const MuonGMR4::SpectrometerSector*, SegmentVec> candidatesPerChamber{};
568 
569  for (std::unique_ptr<Segment>& sortMe : segmentCandidates) {
570  const MuonGMR4::SpectrometerSector* chamb = sortMe->msSector();
571  candidatesPerChamber[chamb].push_back(std::move(sortMe));
572  }
573  segmentCandidates.clear();
574  for (auto& [chamber, resolveMe] : candidatesPerChamber) {
575  SegmentVec resolvedSegments = m_ambiSolver->resolveAmbiguity(gctx, std::move(resolveMe));
576  segmentCandidates.insert(segmentCandidates.end(),
577  std::make_move_iterator(resolvedSegments.begin()),
578  std::make_move_iterator(resolvedSegments.end()));
579  }
580  }
581 }
MuonR4::PrimitiveVec
MuonValR4::IPatternVisualizationTool::PrimitiveVec PrimitiveVec
Definition: SegmentFittingAlg.cxx:60
MuonR4::inv_c
constexpr double inv_c
Definition: SegmentFittingAlg.cxx:32
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
make_hlt_rep.pars
pars
Definition: make_hlt_rep.py:90
LArSamples::FitterData::fitter
const ShapeFitter * fitter
Definition: ShapeFitter.cxx:23
MuonR4::MdtSegmentFitter
Definition: MdtSegmentFitter.h:18
data
char data[hepevt_bytes_allocation_ATLAS]
Definition: HepEvt.cxx:11
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
MuonR4::SegmentFittingAlg::convertToSegment
static std::unique_ptr< Segment > convertToSegment(const Amg::Transform3D &locToGlobTrf, const SegmentSeed *parentSeed, SegmentFitResult &&toConvert)
Converts the fit result into a segment object.
Definition: SegmentFittingAlg.cxx:325
TRTCalib_Extractor.hits
hits
Definition: TRTCalib_Extractor.py:35
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
max
#define max(a, b)
Definition: cfImp.cxx:41
MuonGMR4::SpectrometerSector
A spectrometer sector forms the envelope of all chambers that are placed in the same MS sector & laye...
Definition: SpectrometerSector.h:39
vtune_athena.format
format
Definition: vtune_athena.py:14
MuonR4::SegmentFitHelpers::postFitChi2PerMas
std::pair< std::vector< double >, double > postFitChi2PerMas(const SegmentFit::Parameters &segPars, std::optional< double > arrivalTime, std::vector< std::unique_ptr< CalibratedSpacePoint >> &hits, MsgStream &msg)
Calculates the chi2 per measurement and the chi2 itself after the fit is finished.
Definition: SegmentFitHelperFunctions.cxx:250
MuonR4::SpacePointPerLayerSorter
The SpacePointPerLayerSorter groups the space points by their layer Identifier.
Definition: SpacePointPerLayerSorter.h:14
MuonR4::SegmentFittingAlg::eraseWrongHits
void eraseWrongHits(const ActsGeometryContext &gctx, SegmentFitResult &candidate) const
Removes all hits from the segment which are obvious outliers.
Definition: SegmentFittingAlg.cxx:530
calibdata.chamber
chamber
Definition: calibdata.py:32
VisualizationHelpers.h
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:70
MuonR4::SegmentFittingAlg::fitSegmentHits
SegmentFitResult fitSegmentHits(const EventContext &ctx, const ActsGeometryContext &gctx, const Parameters &startPars, SegmentFitResult::HitVec &&calibHits) const
Executes the segment fit with start parameters.
Definition: SegmentFittingAlg.cxx:145
MuonR4::SegmentFitResult::calibMeasurements
HitVec calibMeasurements
Calibrated measurements used in the fit.
Definition: SegmentFitterEventData.h:68
accumulate
bool accumulate(AccumulateMap &map, std::vector< module_t > const &modules, FPGATrackSimMatrixAccumulator const &acc)
Accumulates an accumulator (e.g.
Definition: FPGATrackSimMatrixAccumulator.cxx:22
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
SegmentFittingAlg.h
deg
#define deg
Definition: SbPolyhedron.cxx:17
MuonR4::SegmentFittingAlg::initialize
virtual StatusCode initialize() override
Definition: SegmentFittingAlg.cxx:63
MuonR4::SegmentFittingAlg::HitVec
SegmentFitResult::HitVec HitVec
Definition: SegmentFittingAlg.h:38
MuonR4::CalibratedSpacePoint::State::Outlier
@ Outlier
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
MuonR4::SegmentAmbiSolver::Config
Definition: SegmentAmbiSolver.h:18
SG::ReadHandleKey< ContainerType >
MuonR4::SegmentFittingAlg::Parameters
SegmentFit::Parameters Parameters
Definition: SegmentFittingAlg.h:41
MuonR4::SegmentFittingAlg::resolveAmbiguities
void resolveAmbiguities(const ActsGeometryContext &gctx, std::vector< std::unique_ptr< Segment >> &segmentCandidates) const
Definition: SegmentFittingAlg.cxx:563
MuonR4::SegmentFit::makeLine
std::pair< Amg::Vector3D, Amg::Vector3D > makeLine(const Parameters &pars)
Returns the parsed parameters into an Eigen line parametrization.
Definition: SegmentFitterEventData.cxx:30
MuonValR4::drawLabel
std::unique_ptr< TLatex > drawLabel(const std::string &text, const double xPos, const double yPos, const unsigned int fontSize=18)
Create a TLatex label,.
Definition: VisualizationHelpers.cxx:22
MuonR4::SegmentFitResult::HitVec
std::vector< HitType > HitVec
Definition: SegmentFitterEventData.h:57
MuonR4::SegmentFit::ParamDefs::phi
@ phi
MuonR4::MdtSegmentSeedGenerator
Helper class to generate valid seeds for the segment fit.
Definition: MdtSegmentSeedGenerator.h:24
MuonR4::MdtSegmentFitter::Config::calibrator
const ISpacePointCalibrator * calibrator
Pointer to the calibrator tool.
Definition: MdtSegmentFitter.h:35
MuonR4::SegmentFitResult::nDoF
int nDoF
degrees of freedom
Definition: SegmentFitterEventData.h:74
MuonR4::SegmentFitHelpers::driftSign
int driftSign(const Amg::Vector3D &posInChamber, const Amg::Vector3D &dirInChamber, const SpacePoint &uncalibHit, MsgStream &msg)
Calculates whether a segement line travereses the tube measurement on the left (-1) or right (1) side...
Definition: SegmentFitHelperFunctions.cxx:221
cm
const double cm
Definition: Simulation/ISF/ISF_FastCaloSim/ISF_FastCaloSimParametrization/tools/FCAL_ChannelMap.cxx:25
MuonR4::Segment::MeasType
std::unique_ptr< CalibratedSpacePoint > MeasType
Definition: MuonSpectrometer/MuonPhaseII/Event/MuonPatternEvent/MuonPatternEvent/Segment.h:22
MuonR4::copy
MuonR4::SegmentFitResult::HitVec copy(const MuonR4::SegmentFitResult::HitVec &hits)
Definition: SegmentFittingAlg.cxx:35
MuonR4::MdtSegmentSeedGenerator::Config
Configuration switches of the module
Definition: MdtSegmentSeedGenerator.h:29
MdtSegmentSeedGenerator.h
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
MuonR4::SegmentFittingAlg::fitSegmentSeed
std::vector< std::unique_ptr< Segment > > fitSegmentSeed(const EventContext &ctx, const ActsGeometryContext &gctx, const SegmentSeed *seed) const
Definition: SegmentFittingAlg.cxx:268
Amg::toString
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Definition: GeoPrimitivesToStringConverter.h:40
MuonR4::SegmentFittingAlg::plugHoles
bool plugHoles(const EventContext &ctx, const ActsGeometryContext &gctx, const SegmentSeed &seed, SegmentFitResult &toRecover) const
Recovery of missed hits.
Definition: SegmentFittingAlg.cxx:410
MuonR4::SegmentFitHelpers::chiSqTermMdt
double chiSqTermMdt(const Amg::Vector3D &posInChamber, const Amg::Vector3D &dirInChamber, const SpacePoint &measurement, MsgStream &msg)
Calculates the chi2 contribuation to a linear segment line from an uncalibrated Mdt measurement.
Definition: SegmentFitHelperFunctions.cxx:49
MuonR4::SegmentFittingAlg::removeOutliers
bool removeOutliers(const EventContext &ctx, const ActsGeometryContext &gctx, const SegmentSeed &seed, SegmentFitResult &fitResult) const
Spot hits with large discrepancy from the estimated parameters and remove them from the list.
Definition: SegmentFittingAlg.cxx:348
MuonR4::SegmentFitResult::converged
bool converged
Is the fit converged.
Definition: SegmentFitterEventData.h:80
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
UtilFunctions.h
SegmentFitHelperFunctions.h
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition: GeoPrimitives.h:46
Amg::transform
Amg::Vector3D transform(Amg::Vector3D &v, Amg::Transform3D &tr)
Transform a point from a Trasformation3D.
Definition: GeoPrimitivesHelpers.h:156
MuonR4::SegmentFitResult::chi2PerMeasurement
std::vector< double > chi2PerMeasurement
Chis per measurement to identify outliers.
Definition: SegmentFitterEventData.h:70
MuonR4::MdtSegmentFitter::Config
Definition: MdtSegmentFitter.h:25
xAOD::Other
@ Other
python.StandardJetMods.pull
pull
Definition: StandardJetMods.py:282
MuonSimHitContainer.h
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
MuonR4::SegmentFitResult::segmentPars
Parameters segmentPars
Final segment parameters.
Definition: SegmentFitterEventData.h:64
MuonR4::SegmentFit::toString
std::string toString(const Parameters &pars)
Definition: SegmentFitterEventData.cxx:63
MuonR4::removeBeamSpot
bool removeBeamSpot(MuonR4::SegmentFitResult::HitVec &hits)
Definition: SegmentFittingAlg.cxx:51
MuonR4::MdtSegmentSeedGenerator::Config::hitPullCut
double hitPullCut
Upper cut on the hit chi2 w.r.t.
Definition: MdtSegmentSeedGenerator.h:35
MuonR4::SegmentFit::ParamDefs::x0
@ x0
DataVector
Derived DataVector<T>.
Definition: DataVector.h:581
ParticleHypothesis.h
ActsGeometryContext
Include the GeoPrimitives which need to be put first.
Definition: ActsGeometryContext.h:27
MuonR4::SegmentFit::ParamDefs::time
@ time
MuonR4::CalibratedSpacePoint::State::Valid
@ Valid
MuonR4::SpacePoint
The muon space point is the combination of two uncalibrated measurements one of them measures the eta...
Definition: MuonSpectrometer/MuonPhaseII/Event/MuonSpacePoint/MuonSpacePoint/SpacePoint.h:18
MuonValR4
Lightweight algorithm to read xAOD MDT sim hits and (fast-digitised) drift circles from SG and fill a...
Definition: IPatternVisualizationTool.h:23
TgcStrip.h
MuonR4::SegmentFit::ParamDefs::y0
@ y0
MuonValR4::IPatternVisualizationTool::PrimitiveVec
std::vector< PrimitivePtr > PrimitiveVec
Definition: IPatternVisualizationTool.h:31
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:221
ActsTrk::to_string
std::string to_string(const DetectorType &type)
Definition: GeometryDefs.h:34
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
MuonR4::SegmentFit::toInt
constexpr int toInt(const ParamDefs p)
Definition: MuonHoughDefs.h:42
RpcMeasurement.h
WriteCaloSwCorrections.cfg
cfg
Definition: WriteCaloSwCorrections.py:23
MuonR4::SegmentFittingAlg::retrieveContainer
StatusCode retrieveContainer(const EventContext &ctx, const SG::ReadHandleKey< ContainerType > &key, const ContainerType *&contToPush) const
Helper method to fetch data from StoreGate.
Definition: SegmentFittingAlg.cxx:131
python.PhysicalConstants.c_light
float c_light
Definition: PhysicalConstants.py:63
Amg::lineDistance
double lineDistance(const AmgVector(N)&posA, const AmgVector(N)&dirA, const AmgVector(N)&posB, const AmgVector(N)&dirB)
: Calculates the shortest distance between two lines
Definition: GeoPrimitivesHelpers.h:308
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
python.copyTCTOutput.locDir
locDir
Definition: copyTCTOutput.py:113
MuonR4
This header ties the generic definitions in this package.
Definition: HoughEventData.h:16
PlotSFuncertainty.calib
calib
Definition: PlotSFuncertainty.py:110
SG::WriteHandle
Definition: StoreGate/StoreGate/WriteHandle.h:76
MuonR4::SegmentVec
SegmentAmbiSolver::SegmentVec SegmentVec
Definition: SegmentAmbiSolver.cxx:10
a
TList * a
Definition: liststreamerinfos.cxx:10
MdtSegmentFitter.h
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
MuonGMR4::SpectrometerSector::localToGlobalTrans
const Amg::Transform3D & localToGlobalTrans(const ActsGeometryContext &gctx) const
Returns the local -> global tarnsformation from the sector.
Definition: SpectrometerSector.cxx:51
GeoPrimitivesHelpers.h
MuonR4::SegmentFitResult::nIter
unsigned int nIter
Number of iterations called to reach the minimum.
Definition: SegmentFitterEventData.h:82
MuonR4::SegmentFitResult::timeFit
bool timeFit
Was the time fitted.
Definition: SegmentFitterEventData.h:60
MuonR4::SegmentSeed
Representation of a segment seed (a fully processed hough maximum) produced by the hough transform.
Definition: SegmentSeed.h:14
if
if(febId1==febId2)
Definition: LArRodBlockPhysicsV0.cxx:567
MuonR4::SegmentSeed::parameters
const Parameters & parameters() const
Returns the parameter array.
Definition: SegmentSeed.cxx:35
MuonR4::SegmentFit::localSegmentPars
Parameters localSegmentPars(const xAOD::MuonSegment &seg)
Returns the localSegPars decoration from a xAODMuon::Segment.
Definition: SegmentFitterEventData.cxx:36
python.BuildSignatureFlags.beamSpot
AthConfigFlags beamSpot(AthConfigFlags flags, str instanceName, str recoMode)
Definition: BuildSignatureFlags.py:454
MuonR4::SegmentFitResult
Definition: SegmentFitterEventData.h:50
MuonR4::SegmentFitResult::chi2
double chi2
chi2 of the fit
Definition: SegmentFitterEventData.h:72
drawLine
void drawLine(std::ostream &os)
Definition: PrintMC.cxx:16
MuonR4::CalibSegmentChi2Minimizer
Definition: CalibSegmentChi2Minimizer.h:21
xAOD::MdtDriftCircle_v1
https://gitlab.cern.ch/atlas/athena/-/blob/master/MuonSpectrometer/MuonReconstruction/MuonRecEvent/Mu...
Definition: MdtDriftCircle_v1.h:21
MuonR4::AmgSymMatrix
const AmgSymMatrix(2) &SpacePoint
Definition: MuonSpectrometer/MuonPhaseII/Event/MuonSpacePoint/src/SpacePoint.cxx:150
MuonR4::SegmentFittingAlg::execute
virtual StatusCode execute(const EventContext &ctx) const override
Definition: SegmentFittingAlg.cxx:74
MuonR4::SegmentFit::ParamDefs::nPars
@ nPars
xAOD::UncalibMeasType::MdtDriftCircleType
@ MdtDriftCircleType
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:26
SpacePointPerLayerSorter.h
MuonR4::SegmentFittingAlg::~SegmentFittingAlg
virtual ~SegmentFittingAlg()
CalibSegmentChi2Minimizer.h
MuonR4::SegmentFit::ParamDefs::theta
@ theta
MuonR4::SegmentSeed::msSector
const MuonGMR4::SpectrometerSector * msSector() const
Returns the associated chamber.
Definition: SegmentSeed.cxx:39
generate::Zero
void Zero(TH1D *hin)
Definition: generate.cxx:32
MuonR4::SegmentFitResult::makeLine
std::pair< Amg::Vector3D, Amg::Vector3D > makeLine() const
Returns the defining parameters as a pair of Amg::Vector3D The first part is the position expressed a...
Definition: SegmentFitterEventData.h:87
mapkey::key
key
Definition: TElectronEfficiencyCorrectionTool.cxx:37
Generated on Mon Nov 18 2024 21:25:46 for ATLAS Offline Software by   doxygen 1.8.18