ATLAS Offline Software
MdtSegmentFitter.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 #include <MuonPatternHelpers/MdtSegmentFitter.h>
5 #include <MuonPatternHelpers/SegmentFitHelperFunctions.h>
6 #include <MuonRecToolInterfacesR4/ISpacePointCalibrator.h>
7 #include <TrkEventPrimitives/ParamDefs.h>
8 #include <EventPrimitives/EventPrimitivesCovarianceHelpers.h>
9 #include <GaudiKernel/PhysicalConstants.h>
10 #include <MuonSpacePoint/UtilFunctions.h>
11 #include <CxxUtils/sincos.h>
12 /* * Residual strip hit: R= (P + <P-H|D>*D -H)
13  * dR
14  * --> -- dP = dP + <dP|D>*D
15  * dP
16  *
17  * dR
18  * -- dD = <P-H|D>dD + <P-H|dD>*D
19  * dD
20  *
21  * chi2 = <R|1./cov^{2} | R>
22  *
23  * dchi2 = <dR | 1./cov^{2} | R> + <R | 1./cov^{2} | dR>
24  */
25 namespace {
26  constexpr double c_inv = 1. / Gaudi::Units::c_light;
27  /* Cut off value for the determinant. Hessian matrices with a determinant smaller than this
28  are considered to be invalid */
29  constexpr double detCutOff = 1.e-8;
30 }
31 
32 namespace MuonR4{
33  using namespace SegmentFit;
34  using HitType = SegmentFitResult::HitType;
35  using HitVec = SegmentFitResult::HitVec;
36 
37 
38  MdtSegmentFitter::Config::RangeArray
39  MdtSegmentFitter::Config::defaultRanges() {
40  RangeArray rng{};
41  constexpr double spatRang = 10.*Gaudi::Units::m;
42  constexpr double timeTange = 25 * Gaudi::Units::ns;
43  rng[toInt(ParamDefs::y0)] = std::array{-spatRang, spatRang};
44  rng[toInt(ParamDefs::x0)] = std::array{-spatRang, spatRang};
45  rng[toInt(ParamDefs::phi)] = std::array{-180.* Gaudi::Units::deg, 180. * Gaudi::Units::deg};
46  rng[toInt(ParamDefs::theta)] = std::array{-90. * Gaudi::Units::deg, 90. * Gaudi::Units::deg};
47  rng[toInt(ParamDefs::time)] = std::array{-timeTange, timeTange};
48  return rng;
49  }
50  MdtSegmentFitter::MdtSegmentFitter(const std::string& name, Config&& config):
51  AthMessaging{name},
52  m_cfg{std::move(config)}{}
53 
54 
55  inline Amg::Vector3D MdtSegmentFitter::partialPlaneIntersect(const Amg::Vector3D& normal, const double offset,
56  const Amg::Vector3D& segPos, const Amg::Vector3D& segDir,
57  const LinePartialArray& linePartials, const ParamDefs fitPar) {
58  const double normDot = normal.dot(segDir);
59  switch (fitPar) {
60  case ParamDefs::phi:
61  case ParamDefs::theta:{
62  const double travelledDist = (offset - segPos.dot(normal)) / normDot;
63  const double partialDist = - travelledDist / normDot * normal.dot(linePartials[toInt(fitPar)]);
64  return travelledDist * linePartials[toInt(fitPar)] + partialDist * segDir;
65  break;
66  } case ParamDefs::y0:
67  case ParamDefs::x0:
68  return linePartials[toInt(fitPar)] - linePartials[toInt(fitPar)].dot(normal) / normDot * segDir;
69  break;
70  default:
71  break;
72  }
73  return Amg::Vector3D::Zero();
74  }
75 
76 
77  inline Amg::Vector3D MdtSegmentFitter::partialClosestApproach(const MuonR4::CalibratedSpacePoint& sp,
78  const Amg::Vector3D& segPos, const Amg::Vector3D& segDir,
79  const LinePartialArray& linePartials, const ParamDefs fitPar) {
80 
81  const Amg::Vector3D& hitDir{sp.directionInChamber()};
82  const Amg::Vector3D& hitPos{sp.positionInChamber()};
83 
84  const double dirDots = hitDir.dot(segDir);
85  const double divisor = (1. - dirDots * dirDots);
86 
87  switch (fitPar) {
88  case ParamDefs::phi:
89  case ParamDefs::theta: {
90  const Amg::Vector3D& segDirPartial{linePartials[toInt(fitPar)]};
91  const Amg::Vector3D AminusB = hitPos - segPos;
92 
93  const double AminusBdotHit = AminusB.dot(hitDir);
94  const double numerator = (AminusB.dot(segDir) - AminusBdotHit * dirDots);
95  const double travelledDist = numerator / divisor;
96  const double partialDirDots = hitDir.dot(segDirPartial);
97 
98 
99  const double derivativeDist = (divisor * (AminusB.dot(segDirPartial) - AminusBdotHit* partialDirDots)
100  + 2. *numerator * dirDots *partialDirDots) / (divisor*divisor);
101 
102  return travelledDist *segDirPartial + derivativeDist * segDir;
103  break;
104  }
105  case ParamDefs::x0:
106  case ParamDefs::y0: {
107  const Amg::Vector3D AminusB = -linePartials[toInt(fitPar)];
108  const double numerator = (AminusB.dot(segDir) - AminusB.dot(hitDir) * dirDots);
109  const double travelledDist = numerator / divisor;
110  return linePartials[toInt(fitPar)] + travelledDist * segDir;
111  }
112  default:
113  break;
114  }
115  return Amg::Vector3D::Zero();
116  }
117  inline void MdtSegmentFitter::updateLinePartials(const Parameters& fitPars, LinePartialArray& linePartials) const{
128  const CxxUtils::sincos theta{fitPars[toInt(ParamDefs::theta)]},
129  phi{fitPars[toInt(ParamDefs::phi)]};
130  linePartials[toInt(ParamDefs::theta)] = Amg::Vector3D{phi.cs*theta.cs, phi.sn*theta.cs, -theta.sn};
131  linePartials[toInt(ParamDefs::phi)] = Amg::Vector3D{-theta.sn *phi.sn, theta.sn*phi.cs, 0};
132  ATH_MSG_VERBOSE("Directional derivatives: dDir/dTheta = "<<Amg::toString(linePartials[toInt(ParamDefs::theta)])
133  <<", dDir/dPhi ="<<Amg::toString(linePartials[toInt(ParamDefs::phi)]));
134  }
135 
136 
137  SegmentFitResult MdtSegmentFitter::fitSegment(const EventContext& ctx,
138  HitVec&& calibHits,
139  const Parameters& startPars,
140  const Amg::Transform3D& localToGlobal) const {
141 
142  const Muon::IMuonIdHelperSvc* idHelperSvc{calibHits[0]->spacePoint()->msSector()->idHelperSvc()};
143  using State = CalibratedSpacePoint::State;
144 
145  if (msgLvl(MSG::VERBOSE)) {
146  std::stringstream hitStream{};
147  const auto [startPos, startDir] = makeLine(startPars);
148  for (const HitType& hit : calibHits) {
149  hitStream<<" **** "<<(hit->type() != xAOD::UncalibMeasType::Other ? idHelperSvc->toString(hit->spacePoint()->identify()): "beamspot" )
150  <<" position: "<<Amg::toString(hit->positionInChamber());
151  if (hit->type() == xAOD::UncalibMeasType::MdtDriftCircleType) {
152  hitStream<<", driftRadius: "<<SegmentFitHelpers::driftSign(startPos,startDir, *hit, msg())*hit->driftRadius() ;
153  }
154  hitStream<<", channel dir: "<<Amg::toString(hit->directionInChamber())<<std::endl;
155  }
156  ATH_MSG_VERBOSE("Start segment fit with parameters "<<toString(startPars)
157  <<", plane location: "<<Amg::toString(localToGlobal)<<std::endl<<hitStream.str());
158 
159  }
160 
161  SegmentFitResult fitResult{};
162  fitResult.segmentPars = startPars;
163  fitResult.timeFit = m_cfg.doTimeFit;
164  fitResult.calibMeasurements = std::move(calibHits);
165 
166  Parameters gradient{AmgVector(5)::Zero()}, prevGrad{AmgVector(5)::Zero()}, prevPars{AmgVector(5)::Zero()};
167  AmgSymMatrix(5) hessian{AmgSymMatrix(5)::Zero()};
168 
170  LinePartialArray linePartials{make_array<Amg::Vector3D, toInt(ParamDefs::nPars)>(Amg::Vector3D::Zero())};
171  linePartials[toInt(ParamDefs::x0)] = Amg::Vector3D::UnitX();
172  linePartials[toInt(ParamDefs::y0)] = Amg::Vector3D::UnitY();
173 
175  LinePartialArray partialsResidual{make_array<Amg::Vector3D,toInt(ParamDefs::nPars)>(Amg::Vector3D::Zero())};
176  Amg::Vector3D residual{Amg::Vector3D::Zero()};
177 
178  unsigned int noChangeIter{0};
179  while (fitResult.nIter++ < m_cfg.nMaxCalls) {
180  ATH_MSG_VERBOSE("Iteration: "<<fitResult.nIter<<" parameters: "<<toString(fitResult.segmentPars)<<"chi2: "<<fitResult.chi2);
182  const auto [segPos, segDir] = fitResult.makeLine();
184  fitResult.calibMeasurements = m_cfg.calibrator->calibrate(ctx, std::move(fitResult.calibMeasurements), segPos, segDir,
185  fitResult.segmentPars[toInt(ParamDefs::time)]);
187  fitResult.nPhiMeas = fitResult.nDoF = fitResult.nTimeMeas = 0;
188 
189  for (const HitType& hit : fitResult.calibMeasurements) {
190  if (hit->fitState() != State::Valid){
191  continue;
192  }
193  fitResult.nPhiMeas+= hit->measuresPhi();
194  fitResult.nDoF+= hit->measuresPhi();
195  fitResult.nDoF+= hit->measuresEta();
197  fitResult.nDoF += (m_cfg.doTimeFit && hit->type() != xAOD::UncalibMeasType::MdtDriftCircleType && hit->measuresTime());
198  fitResult.nTimeMeas+=hit->measuresTime();
199  }
200  if (!fitResult.nDoF) {
201  ATH_MSG_WARNING("TSCHUUUUUUUUSS Measurements...");
202  break;
203  }
204  if (!fitResult.nPhiMeas) {
205  ATH_MSG_VERBOSE("No phi measurements are left.");
206  fitResult.segmentPars[toInt(ParamDefs::phi)] = 90. * Gaudi::Units::deg;
207  fitResult.segmentPars[toInt(ParamDefs::x0)] = 0;
208  }
209 
210  fitResult.nDoF = fitResult.nDoF - 2 - (fitResult.nPhiMeas > 0 ? 2 : 0);
211 
213  if (fitResult.timeFit && fitResult.nDoF <= 1) {
214  fitResult.timeFit = false;
215  ATH_MSG_DEBUG("Switch of the time fit because nDoF: "<<fitResult.nDoF);
216  fitResult.segmentPars[toInt(ParamDefs::time)] = 0.;
218  fitResult.calibMeasurements = m_cfg.calibrator->calibrate(ctx, std::move(fitResult.calibMeasurements), segPos, segDir,
219  fitResult.segmentPars[toInt(ParamDefs::time)]);
221  } else if (!fitResult.timeFit && m_cfg.doTimeFit) {
222  ATH_MSG_DEBUG("Somehow a measurement is on the narrow ridge of validity. Let's try if the time can be fitted now ");
223  fitResult.timeFit = true;
224  }
226  fitResult.chi2 = 0;
227  hessian.setZero();
228  gradient.setZero();
230  updateLinePartials(fitResult.segmentPars, linePartials);
231 
233  for (const HitType& hit : fitResult.calibMeasurements) {
234  if (hit->fitState() != State::Valid) {
235  continue;
236  }
237  const Amg::Vector3D& hitPos{hit->positionInChamber()};
238  const Amg::Vector3D& hitDir{hit->directionInChamber()};
239  const unsigned int dim = hit->type() == xAOD::UncalibMeasType::Other ? 2 : hit->spacePoint()->dimension();
241  const int start = toInt(fitResult.nPhiMeas ? ParamDefs::phi : ParamDefs::theta);
242  switch (hit->type()) {
243  case xAOD::UncalibMeasType::MdtDriftCircleType: {
245  const double travelledDist = Amg::intersect(hitPos, hitDir, segPos, segDir).value_or(0);
246  const Amg::Vector3D closePointSeg = segPos + travelledDist* segDir;
248  const Amg::Vector3D closePointWire = hitPos + hitDir.dot(closePointSeg - hitPos) * hitDir;
249 
250  const Amg::Vector3D lineConnect = (closePointWire - closePointSeg);
251  const double lineDist = lineConnect.mag();
253  residual[toInt(AxisDefs::t0)] = 0.;
254  residual[toInt(AxisDefs::eta)] = lineDist - hit->driftRadius();
255  const Amg::Vector3D globApproach = (localToGlobal*closePointSeg).unit();
256  // For twin tubes, the hit position is not updated during the calibration. Make use of
257  // additional constraint to fit phi
258  residual[toInt(AxisDefs::phi)] = dim == 2 ? (hitPos - closePointSeg).x() : 0.;
260  const double driftV{m_cfg.calibrator->driftVelocity(ctx, *hit)};
261  for (int p = start; p >= 0; --p) {
262  const ParamDefs par{static_cast<ParamDefs>(p)};
264  const Amg::Vector3D partialClosePointOnSeg = partialClosestApproach(*hit, segPos, segDir, linePartials, par);
266  const Amg::Vector3D partialClosePointW = hitDir.dot(partialClosePointOnSeg) * hitDir;
267 
268  const double dR = -driftV * c_inv * globApproach.dot(linePartials[p]);
269 
270  partialsResidual[p][toInt(AxisDefs::phi)] = dim ==2 ? -partialClosePointOnSeg.x() : 0.;
271  partialsResidual[p][toInt(AxisDefs::eta)] = lineConnect.dot(partialClosePointW - partialClosePointOnSeg) / lineDist;
272  partialsResidual[p][toInt(AxisDefs::t0)] = 0.;
273  ATH_MSG_VERBOSE("Partial derivative of "<<idHelperSvc->toString(hit->spacePoint()->identify())
274  <<" residual "<<Amg::toString(residual)<<" w.r.t "<<toString(par)<<"="
275  <<Amg::toString(partialsResidual[p])<<", drift velocity: "<<driftV
276  <<" in terms of beta: "<<(driftV * c_inv)<<" -> dR: "<<dR);
277 
278  }
280  if (fitResult.timeFit) {
281  constexpr int par = toInt(ParamDefs::time);
282  partialsResidual[par] = driftV * Amg::Vector3D::Unit(toInt(AxisDefs::eta));
283  ATH_MSG_VERBOSE("Partial derivative of "<<idHelperSvc->toString(hit->spacePoint()->identify())
284  <<", driftRadius: "<<hit->driftRadius()<<", residual "<<Amg::toString(residual)
285  <<" w.r.t "<<toString(ParamDefs::time)<<"="<<Amg::toString(partialsResidual[par]));
286  }
287  break;
288  }
289  case xAOD::UncalibMeasType::RpcStripType:
290  case xAOD::UncalibMeasType::TgcStripType:{
291  const Amg::Vector3D normal = hit->spacePoint()->planeNormal();
292  const double planeOffSet = normal.dot(hit->positionInChamber());
293  const Amg::Vector3D planeIsect = segPos + Amg::intersect<3>(segPos, segDir, normal, planeOffSet).value_or(0)* segDir;
294 
296  residual.block<2,1>(0,0) = (hitPos - planeIsect).block<2,1>(0,0);
297 
298  if (fitResult.timeFit && hit->measuresTime()) {
300  const double totFlightDist = (localToGlobal * planeIsect).mag();
301  residual[toInt(AxisDefs::t0)] = hit->time() - totFlightDist * c_inv - fitResult.segmentPars[toInt(ParamDefs::time)];
302  }
303 
304  for (int p = start; p >= 0; --p) {
305  const ParamDefs par{static_cast<ParamDefs>(p)};
306  partialsResidual[toInt(par)].block<2,1>(0, 0) = - partialPlaneIntersect(normal, planeOffSet,
307  segPos, segDir,linePartials, par).block<2,1>(0,0);
308 
309  if (fitResult.timeFit && hit->measuresTime()) {
310  partialsResidual[toInt(par)][toInt(AxisDefs::t0)] = -partialsResidual[toInt(par)].perp() * c_inv;
311  }
312  ATH_MSG_VERBOSE("Partial derivative of "<<idHelperSvc->toString(hit->spacePoint()->identify())
313  <<" residual "<<Amg::toString(residual)<<" "<<
314  Amg::toString(multiply(inverse(hit->covariance()), residual))
315  <<" "<<std::endl<<toString(hit->covariance())<<std::endl<<" w.r.t "<<toString(par)<<"="
316  <<Amg::toString(partialsResidual[toInt(par)]));
317 
318  }
319  if (fitResult.timeFit && hit->measuresTime()) {
320  constexpr ParamDefs par = ParamDefs::time;
321  partialsResidual[toInt(par)] = -Amg::Vector3D::Unit(toInt(AxisDefs::t0));
322  ATH_MSG_VERBOSE("Partial derivative of "<<idHelperSvc->toString(hit->spacePoint()->identify())
323  <<" residual "<<Amg::toString(residual)<<" w.r.t "<<toString(par)<<"="
324  <<Amg::toString(partialsResidual[toInt(par)]));
325  }
326  break;
327  } case xAOD::UncalibMeasType::Other:{
328  static const Amg::Vector3D normal = Amg::Vector3D::UnitZ();
329  const double planeOffSet = normal.dot(hit->positionInChamber());
330  const Amg::Vector3D planeIsect = segPos + Amg::intersect<3>(segPos, segDir, normal, planeOffSet).value_or(0)* segDir;
331 
332  residual.block<2,1>(0,0) = (hitPos - planeIsect).block<2,1>(0,0);
333  residual[toInt(AxisDefs::t0)] =0.;
334  for (int p = start; p >= 0; --p) {
335  const ParamDefs par{static_cast<ParamDefs>(p)};
336  partialsResidual[toInt(par)].block<2,1>(0, 0) = - partialPlaneIntersect(normal, planeOffSet,
337  segPos, segDir,linePartials, par).block<2,1>(0,0);
338 
339  partialsResidual[toInt(par)][toInt(AxisDefs::t0)] = 0;
340  ATH_MSG_VERBOSE("Partial derivative of "<<idHelperSvc->toString(hit->spacePoint()->identify())
341  <<" residual "<<Amg::toString(residual)<<" w.r.t "<<toString(par)<<"="
342  <<Amg::toString(partialsResidual[toInt(par)]));
343 
344  }
345  break;
346  } default:
347  ATH_MSG_WARNING("MdtSegmentFitter() - Unsupported measurment type" <<typeid(*hit->spacePoint()).name());
348  }
349 
350  ATH_MSG_VERBOSE("Update derivatives for hit "<< (hit->spacePoint() ? idHelperSvc->toString(hit->spacePoint()->identify()) : "beamspot"));
351  updateDerivatives(residual, partialsResidual, hit->covariance(), gradient, hessian, fitResult.chi2,
352  fitResult.timeFit && hit->measuresTime() ? toInt(ParamDefs::time) : start);
353  }
354 
356  for (int k =1; k < 5 - (!fitResult.timeFit); ++k){
357  for (int l = 0; l< k; ++l){
358  hessian(l,k) = hessian(k,l);
359  }
360  }
362  if (gradient.mag() < m_cfg.tolerance) {
363  fitResult.converged = true;
364  ATH_MSG_VERBOSE("Fit converged after "<<fitResult.nIter<<" iterations with "<<fitResult.chi2);
365  break;
366  }
367  ATH_MSG_VERBOSE("Chi2: "<<fitResult.chi2<<", gradient: "<<Amg::toString(gradient)<<"hessian: "<<std::endl<<hessian);
369  UpdateStatus paramUpdate{UpdateStatus::outOfBounds};
370  if (!fitResult.nPhiMeas && !fitResult.timeFit) {
371  paramUpdate = updateParameters<2>(fitResult.segmentPars, prevPars, gradient, prevGrad, hessian);
372  } else if (!fitResult.nPhiMeas && fitResult.timeFit) {
375  hessian.col(2).swap(hessian.col(toInt(ParamDefs::time)));
376  hessian.row(2).swap(hessian.row(toInt(ParamDefs::time)));
377  std::swap(gradient[2], gradient[toInt(ParamDefs::time)]);
378  std::swap(fitResult.segmentPars[2], fitResult.segmentPars[toInt(ParamDefs::time)]);
379 
380  paramUpdate = updateParameters<3>(fitResult.segmentPars, prevPars, gradient, prevGrad, hessian);
381 
382  std::swap(fitResult.segmentPars[2], fitResult.segmentPars[toInt(ParamDefs::time)]);
383  hessian.col(2).swap(hessian.col(toInt(ParamDefs::time)));
384  hessian.row(2).swap(hessian.row(toInt(ParamDefs::time)));
385  std::swap(gradient[2], gradient[toInt(ParamDefs::time)]);
386  } else if (fitResult.nPhiMeas && !fitResult.timeFit) {
387  paramUpdate = updateParameters<4>(fitResult.segmentPars, prevPars, gradient, prevGrad, hessian);
388  } else if (fitResult.nPhiMeas && fitResult.timeFit) {
389  paramUpdate = updateParameters<5>(fitResult.segmentPars, prevPars, gradient, prevGrad, hessian);
390  }
391  if (paramUpdate == UpdateStatus::noChange) {
392  if ((++noChangeIter) >= m_cfg.noMoveIter){
393  fitResult.converged = true;
394  break;
395  }
396  } else if (paramUpdate == UpdateStatus::allOkay) {
397  noChangeIter = 0;
398  } else if (paramUpdate == UpdateStatus::outOfBounds){
399  fitResult.chi2PerMeasurement.resize(fitResult.calibMeasurements.size(), -1.);
400  return fitResult;
401  }
402  }
403 
405  fitResult.nDoF-=fitResult.timeFit;
406 
408  const auto [segPos, segDir] = fitResult.makeLine();
409  fitResult.chi2 =0.;
410 
411  std::optional<double> toF = fitResult.timeFit ? std::make_optional<double>((localToGlobal * segPos).mag() * c_inv) : std::nullopt;
413  std::ranges::stable_sort(fitResult.calibMeasurements, [](const HitType&a, const HitType& b){
414  return a->positionInChamber().z() > b->positionInChamber().z();
415  });
416 
417  auto [chi2Term, chi2]= SegmentFitHelpers::postFitChi2PerMas(fitResult.segmentPars, toF, fitResult.calibMeasurements, msg());
418  fitResult.chi2PerMeasurement = std::move(chi2Term);
419  fitResult.chi2 = chi2;
421  if (!fitResult.nPhiMeas&& !fitResult.timeFit) {
422  blockCovariance<2>(std::move(hessian), /* std::move(gradient), std::move(prevGrad), */ fitResult.segmentParErrs);
423  }else if (!fitResult.nPhiMeas && fitResult.timeFit) {
424  hessian.col(2).swap(hessian.col(toInt(ParamDefs::time)));
425  hessian.row(2).swap(hessian.row(toInt(ParamDefs::time)));
426  blockCovariance<3>(std::move(hessian),/* std::move(gradient), std::move(prevGrad), */ fitResult.segmentParErrs);
427  fitResult.segmentParErrs.col(2).swap(fitResult.segmentParErrs.col(toInt(ParamDefs::time)));
428  fitResult.segmentParErrs.row(2).swap(fitResult.segmentParErrs.row(toInt(ParamDefs::time)));
429  } else if (fitResult.nPhiMeas) {
430  blockCovariance<4>(std::move(hessian), /* std::move(gradient), std::move(prevGrad), */ fitResult.segmentParErrs);
431  } else if (fitResult.nPhiMeas && fitResult.timeFit) {
432  blockCovariance<5>(std::move(hessian),/* std::move(gradient), std::move(prevGrad), */ fitResult.segmentParErrs);
433  }
434  return fitResult;
435  }
436 
437  void MdtSegmentFitter::updateDerivatives(const Amg::Vector3D& residual,
438  const LinePartialArray& residualPartials,
439  const MeasCov_t& measCovariance,
440  AmgVector(5)& gradient,
441  AmgSymMatrix(5)& hessian,
442  double& chi2, int startPar) const {
443 
444  const MeasCov_t invCov{inverse(measCovariance)};
445 
446  const Amg::Vector3D covRes = multiply(invCov, residual);
447  chi2 += covRes.dot(residual);
448  for (int p = startPar; p >=0 ; --p) {
449  gradient[p] +=2.*covRes.dot(residualPartials[p]);
450  for (int k=p; k>=0; --k) {
451  hessian(p,k)+= 2.*contract(invCov, residualPartials[p], residualPartials[k]);
452  }
453  }
454  ATH_MSG_VERBOSE("After derivative update --- chi2: "<<chi2<<"("<<covRes.dot(residual)<<"), gradient: "
455  <<toString(gradient)<<", Hessian:\n"<<hessian<<", measurement covariance\n"<<toString(invCov));
456  }
457 
458  template <unsigned int nDim>
459  void MdtSegmentFitter::blockCovariance(const AmgSymMatrix(5)& hessian,
460  SegmentFit::Covariance& covariance) const {
461 
462  covariance.setIdentity();
463  AmgSymMatrix(nDim) miniHessian = hessian.block<nDim, nDim>(0,0);
464  if (std::abs(miniHessian.determinant()) <= detCutOff) {
465  ATH_MSG_VERBOSE("Boeser mini hessian ("<<miniHessian.determinant()<<")\n"<<miniHessian<<"\n\n"<<hessian);
466  return;
467  }
468  ATH_MSG_VERBOSE("Hessian matrix: \n"<<hessian<<",\nblock Hessian:\n"<<miniHessian<<",\n determinant: "<<miniHessian.determinant());
469  covariance.block<nDim,nDim>(0,0) = miniHessian.inverse();
470  ATH_MSG_VERBOSE("covariance: \n"<<covariance);
471  }
472 
473  template <unsigned int nDim>
474  MdtSegmentFitter::UpdateStatus
475  MdtSegmentFitter::updateParameters(Parameters& currPars, Parameters& prevPars,
476  Parameters& currGrad, Parameters& prevGrad,
477  const AmgSymMatrix(5)& currentHessian) const {
478 
479  AmgSymMatrix(nDim) miniHessian = currentHessian.block<nDim, nDim>(0,0);
480  ATH_MSG_VERBOSE("Parameter update -- \ncurrenPars: "<<toString(currPars)<<", \ngradient: "<<toString(currGrad)
481  <<", hessian ("<<miniHessian.determinant()<<")"<<std::endl<<miniHessian);
482 
483  double changeMag{0.};
484  if (std::abs(miniHessian.determinant()) > detCutOff) {
485  prevPars.block<nDim,1>(0,0) = currPars.block<nDim,1>(0,0);
486  // Update the parameters accrodingly to the hessian
487  const AmgVector(nDim) updateMe = miniHessian.inverse()* currGrad.block<nDim, 1>(0,0);
488  changeMag = std::sqrt(updateMe.dot(updateMe));
489  currPars.block<nDim,1>(0,0) -= updateMe;
490  prevGrad.block<nDim,1>(0,0) = currGrad.block<nDim,1>(0,0);
491  ATH_MSG_VERBOSE("Hessian inverse:\n"<<miniHessian.inverse()<<"\n\nUpdate the parameters by -"
492  <<Amg::toString(updateMe));
493  } else {
494  const AmgVector(nDim) gradDiff = (currGrad - prevGrad).block<nDim,1>(0,0);
495  const double gradDiffMag = gradDiff.mag2();
496  const double gamma = std::abs((currPars - prevPars).block<nDim,1>(0,0).dot(gradDiff))
497  / gradDiffMag > std::numeric_limits<float>::epsilon() ? gradDiffMag : 1.;
498  ATH_MSG_VERBOSE("Hessian determinant invalid. Try deepest descent - \nprev parameters: "
499  <<toString(prevPars)<<",\nprevious gradient: "<<toString(prevGrad)<<", gamma: "<<gamma);
500  prevPars.block<nDim, 1>(0,0) = currPars.block<nDim, 1>(0,0);
501  currPars.block<nDim, 1>(0,0) -= gamma* currGrad.block<nDim, 1>(0,0);
502  prevGrad.block<nDim, 1>(0,0) = currGrad.block<nDim,1>(0,0);
503  changeMag = std::abs(gamma) * currGrad.block<nDim, 1>(0,0).mag();
504  }
506  unsigned int nOutOfBound{0};
507  for (unsigned int p =0; p< nDim; ++p) {
508  if (m_cfg.ranges[p][0] > currPars[p] || m_cfg.ranges[p][1]< currPars[p]) {
509  ATH_MSG_WARNING("The "<<p<<"-th parameter "<<toString(static_cast<ParamDefs>(p))<<" is out of range "<<currPars[p]
510  <<"["<<m_cfg.ranges[p][0]<<"-"<<m_cfg.ranges[p][1]<<"]");
511  ++nOutOfBound;
512  }
513  currPars[p] = std::clamp(currPars[p], m_cfg.ranges[p][0], m_cfg.ranges[p][1]);
514  }
515  if (nOutOfBound > m_cfg.nParsOutOfBounds){
516  return UpdateStatus::outOfBounds;
517  }
518  if (changeMag <= m_cfg.tolerance) {
519  return UpdateStatus::noChange;
520  }
521  return UpdateStatus::allOkay;
522  }
523 }
MuonR4::SegmentFit::ParamDefs
ParamDefs
This file defines the parameter enums in the Trk namespace.
Definition: MuonHoughDefs.h:29
MuonR4::MdtSegmentFitter::blockCovariance
void blockCovariance(const AmgSymMatrix(5)&hessian, SegmentFit::Covariance &covariance) const
Definition: MdtSegmentFitter.cxx:459
yodamerge_tmp.dim
dim
Definition: yodamerge_tmp.py:239
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
MuonR4::MdtSegmentFitter::UpdateStatus
UpdateStatus
Definition: MdtSegmentFitter.h:118
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
ClusterSeg::residual
@ residual
Definition: ClusterNtuple.h:20
MuonR4::multiply
Amg::Vector3D multiply(const CalibratedSpacePoint::Covariance_t &mat, const Amg::Vector3D &v)
Multiplies a 3D vector with the covariance matrix which can be either 2x2 or 3x3 matrix.
Definition: MuonSpectrometer/MuonPhaseII/Event/MuonSpacePoint/src/UtilFunctions.cxx:40
MuonR4::MdtSegmentFitter::Parameters
SegmentFit::Parameters Parameters
Definition: MdtSegmentFitter.h:21
MuonR4::MdtSegmentFitter::partialPlaneIntersect
static Amg::Vector3D partialPlaneIntersect(const Amg::Vector3D &normal, const double offset, const Amg::Vector3D &segPos, const Amg::Vector3D &segDir, const LinePartialArray &linePartials, const ParamDefs fitPar)
Calculates the partial derivative of the intersection point between the segment line and the measurem...
Definition: MdtSegmentFitter.cxx:55
mergePhysValFiles.start
start
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:14
MuonR4::MdtSegmentFitter::UpdateStatus::outOfBounds
@ outOfBounds
compute_lumi.divisor
divisor
Definition: compute_lumi.py:60
sincos.h
Helper to simultaneously calculate sin and cos of the same angle.
deg
#define deg
Definition: SbPolyhedron.cxx:17
MuonR4::HitVec
SpacePointPerLayerSorter::HitVec HitVec
Definition: SpacePointPerLayerSorter.cxx:9
ISpacePointCalibrator.h
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:158
MuonR4::ISpacePointCalibrator::calibrate
virtual CalibSpacePointPtr calibrate(const EventContext &ctx, const SpacePoint *spacePoint, const Amg::Vector3D &seedPosInChamb, const Amg::Vector3D &seedDirInChamb, const double timeDelay) const =0
Calibrates a single space point.
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
MuonR4::ISpacePointCalibrator::driftVelocity
virtual double driftVelocity(const EventContext &ctx, const CalibratedSpacePoint &spacePoint) const =0
Returns the drift velocity for a given space point.
MuonR4::MdtSegmentFitter::Config::noMoveIter
unsigned int noMoveIter
How many iterations with changes below tolerance.
Definition: MdtSegmentFitter.h:37
D3PDTest::rng
uint32_t rng()
Definition: FillerAlg.cxx:40
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
x
#define x
ParamDefs.h
MuonR4::MdtSegmentFitter::Config::tolerance
double tolerance
Gradient cut off.
Definition: MdtSegmentFitter.h:29
MuonR4::toString
std::string toString(const CalibratedSpacePoint::Covariance_t &mat)
Returns the matrix in string.
Definition: MuonSpectrometer/MuonPhaseII/Event/MuonSpacePoint/src/UtilFunctions.cxx:75
MuonR4::MdtSegmentFitter::fitSegment
SegmentFitResult fitSegment(const EventContext &ctx, HitVec &&calibHits, const Parameters &startPars, const Amg::Transform3D &localToGlobal) const
Definition: MdtSegmentFitter.cxx:137
config
Definition: PhysicsAnalysis/AnalysisCommon/AssociationUtils/python/config.py:1
MuonR4::SegmentFitResult::HitVec
std::vector< HitType > HitVec
Definition: SegmentFitterEventData.h:57
MuonR4::SegmentFit::ParamDefs::phi
@ phi
xAOD::UncalibMeasType::TgcStripType
@ TgcStripType
MuonR4::MdtSegmentFitter::Config::calibrator
const ISpacePointCalibrator * calibrator
Pointer to the calibrator tool.
Definition: MdtSegmentFitter.h:35
MuonR4::MdtSegmentFitter::partialClosestApproach
static Amg::Vector3D partialClosestApproach(const MuonR4::CalibratedSpacePoint &sp, const Amg::Vector3D &segPos, const Amg::Vector3D &segDir, const LinePartialArray &linePartials, const ParamDefs fitPar)
Calculates the partial derivative of the point of closest approach to the meausrement's wire w....
Definition: MdtSegmentFitter.cxx:77
MuonR4::MdtSegmentFitter::Config::ranges
RangeArray ranges
Definition: MdtSegmentFitter.h:43
MuonR4::SegmentFit::AxisDefs::phi
@ phi
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
MuonR4::SegmentFit::AxisDefs::t0
@ t0
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
Amg::toString
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Definition: GeoPrimitivesToStringConverter.h:40
MuonR4::MdtSegmentFitter::UpdateStatus::noChange
@ noChange
AthMessaging::msgLvl
bool msgLvl(const MSG::Level lvl) const
Test the output level.
Definition: AthMessaging.h:151
TrigVtx::gamma
@ gamma
Definition: TrigParticleTable.h:26
MuonR4::inverse
CalibratedSpacePoint::Covariance_t inverse(const CalibratedSpacePoint::Covariance_t &mat)
Inverts the parsed matrix.
Definition: MuonSpectrometer/MuonPhaseII/Event/MuonSpacePoint/src/UtilFunctions.cxx:65
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
UtilFunctions.h
MuonR4::MdtSegmentFitter::MeasCov_t
CalibratedSpacePoint::Covariance_t MeasCov_t
Updates the chi2, its Gradient & Hessian from the measurement residual.
Definition: MdtSegmentFitter.h:108
SegmentFitHelperFunctions.h
AmgVector
AmgVector(4) T2BSTrackFilterTool
Definition: T2BSTrackFilterTool.cxx:114
MuonR4::MdtSegmentFitter::HitType
std::unique_ptr< CalibratedSpacePoint > HitType
Definition: MdtSegmentFitter.h:22
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition: GeoPrimitives.h:46
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:522
MuonR4::MdtSegmentFitter::Config
Definition: MdtSegmentFitter.h:25
MuonR4::CalibratedSpacePoint::State
State
State flag to distinguish different space point states.
Definition: CalibratedSpacePoint.h:24
AthMessaging
Class to provide easy MsgStream access and capabilities.
Definition: AthMessaging.h:55
MuonR4::MdtSegmentFitter::Config::nParsOutOfBounds
unsigned int nParsOutOfBounds
Abort the fit as soon as more than n parameters leave the fit range.
Definition: MdtSegmentFitter.h:33
xAOD::Other
@ Other
MuonR4::SegmentFitResult::segmentPars
Parameters segmentPars
Final segment parameters.
Definition: SegmentFitterEventData.h:64
MuonR4::MdtSegmentFitter::updateDerivatives
void updateDerivatives(const Amg::Vector3D &residual, const LinePartialArray &residualPartials, const MeasCov_t &measCovariance, AmgVector(5)&gradient, AmgSymMatrix(5)&hessian, double &chi2, int startPar) const
Definition: MdtSegmentFitter.cxx:437
MuonR4::SegmentFit::ParamDefs::x0
@ x0
WriteCalibToCool.swap
swap
Definition: WriteCalibToCool.py:94
HitType
HitType
Definition: FPGATrackSimHit.h:38
MuonR4::SegmentFit::AxisDefs::eta
@ eta
dot.dot
def dot(G, fn, nodesToHighlight=[])
Definition: dot.py:5
AthMessaging::msg
MsgStream & msg() const
The standard message stream.
Definition: AthMessaging.h:164
MuonR4::SegmentFit::ParamDefs::time
@ time
MuonR4::CalibratedSpacePoint::State::Valid
@ Valid
lumiFormat.array
array
Definition: lumiFormat.py:91
MuonR4::MdtSegmentFitter::m_cfg
Config m_cfg
Definition: MdtSegmentFitter.h:56
MuonR4::MdtSegmentFitter::updateLinePartials
void updateLinePartials(const Parameters &fitPars, LinePartialArray &linePartials) const
Updates the partial derivaitves of the line w.r.t the fit parameters.
Definition: MdtSegmentFitter.cxx:117
MuonR4::SegmentFit::ParamDefs::y0
@ y0
MuonR4::MdtSegmentFitter::MdtSegmentFitter
MdtSegmentFitter(const std::string &name, Config &&config)
Definition: MdtSegmentFitter.cxx:50
MuonR4::MdtSegmentFitter::Config::nMaxCalls
unsigned int nMaxCalls
How many calls shall be executed.
Definition: MdtSegmentFitter.h:27
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:221
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:29
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
MuonR4::SegmentFit::toInt
constexpr int toInt(const ParamDefs p)
Definition: MuonHoughDefs.h:42
MuonR4::CalibratedSpacePoint::positionInChamber
const Amg::Vector3D & positionInChamber() const
The position of the calibrated space point inside the chamber.
Definition: CalibratedSpacePoint.cxx:21
MuonR4::MdtSegmentFitter::UpdateStatus::allOkay
@ allOkay
MuonR4::SegmentFitResult::HitType
std::unique_ptr< CalibratedSpacePoint > HitType
Definition: SegmentFitterEventData.h:56
python.PhysicalConstants.c_light
float c_light
Definition: PhysicalConstants.py:63
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
MuonR4
This header ties the generic definitions in this package.
Definition: HoughEventData.h:16
EventPrimitivesCovarianceHelpers.h
MuonR4::MdtSegmentFitter::Config::RangeArray
std::array< std::array< double, 2 >, SegmentFit::toInt(ParamDefs::nPars)> RangeArray
Allowed parameter ranges.
Definition: MdtSegmentFitter.h:39
MuonR4::MdtSegmentFitter::LinePartialArray
std::array< Amg::Vector3D, toInt(ParamDefs::nPars)> LinePartialArray
Store the partial derivative of the line w.r.t.
Definition: MdtSegmentFitter.h:61
MuonR4::SegmentFit::Covariance
AmgSymMatrix(toInt(ParamDefs::nPars)) Covariance
Definition: MuonHoughDefs.h:49
a
TList * a
Definition: liststreamerinfos.cxx:10
MdtSegmentFitter.h
Amg::intersect
std::optional< double > intersect(const AmgVector(N)&posA, const AmgVector(N)&dirA, const AmgVector(N)&posB, const AmgVector(N)&dirB)
Calculates the point of closest approach of two lines.
Definition: GeoPrimitivesHelpers.h:325
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
MuonR4::MdtSegmentFitter::Config::doTimeFit
bool doTimeFit
Switch toggling whether the T0 shall be fitted or not.
Definition: MdtSegmentFitter.h:31
MuonR4::MdtSegmentFitter::HitVec
std::vector< HitType > HitVec
Definition: MdtSegmentFitter.h:23
unit
const PlainObject unit() const
This is a plugin that makes Eigen look like CLHEP & defines some convenience methods.
Definition: AmgMatrixBasePlugin.h:21
MuonR4::CalibratedSpacePoint
The calibrated Space point is created during the calibration process.
Definition: CalibratedSpacePoint.h:15
convertTimingResiduals.offset
offset
Definition: convertTimingResiduals.py:71
CxxUtils::sincos
Helper to simultaneously calculate sin and cos of the same angle.
Definition: sincos.h:76
Muon::IMuonIdHelperSvc
Interface for Helper service that creates muon Identifiers and can be used to print Identifiers.
Definition: IMuonIdHelperSvc.h:26
python.SystemOfUnits.ns
int ns
Definition: SystemOfUnits.py:130
MuonR4::SegmentFitResult
Definition: SegmentFitterEventData.h:50
python.Constants.VERBOSE
int VERBOSE
Definition: Control/AthenaCommon/python/Constants.py:14
MuonR4::MdtSegmentFitter::Config::defaultRanges
static RangeArray defaultRanges()
Function that returns a set of predefined ranges for testing.
Definition: MdtSegmentFitter.cxx:39
MuonR4::AmgSymMatrix
const AmgSymMatrix(2) &SpacePoint
Definition: MuonSpectrometer/MuonPhaseII/Event/MuonSpacePoint/src/SpacePoint.cxx:150
MuonR4::MdtSegmentFitter::updateParameters
UpdateStatus updateParameters(Parameters &currentPars, Parameters &previousPars, Parameters &currGrad, Parameters &prevGrad, const AmgSymMatrix(5)&hessian) const
Update step of the segment parameters using the Hessian and the gradient.
Definition: MdtSegmentFitter.cxx:475
MuonR4::CalibratedSpacePoint::directionInChamber
const Amg::Vector3D & directionInChamber() const
The direction of the calibrated space point inside the chamber.
Definition: CalibratedSpacePoint.cxx:24
xAOD::UncalibMeasType::RpcStripType
@ RpcStripType
xAOD::UncalibMeasType::MdtDriftCircleType
@ MdtDriftCircleType
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:26
MuonR4::SegmentFit::ParamDefs::theta
@ theta
fitman.k
k
Definition: fitman.py:528
generate::Zero
void Zero(TH1D *hin)
Definition: generate.cxx:32
MuonR4::contract
double contract(const CalibratedSpacePoint::Covariance_t &mat, const Amg::Vector3D &a, const Amg::Vector3D &b)
Definition: MuonSpectrometer/MuonPhaseII/Event/MuonSpacePoint/src/UtilFunctions.cxx:13
Generated on Thu Nov 7 2024 21:20:46 for ATLAS Offline Software by   doxygen 1.8.18