ATLAS Offline Software
MdtSegmentFitter.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2025 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 
13 #include <format>
14 /* * Residual strip hit: R= (P + <P-H|D>*D -H)
15  * dR
16  * --> -- dP = dP + <dP|D>*D
17  * dP
18  *
19  * dR
20  * -- dD = <P-H|D>dD + <P-H|dD>*D
21  * dD
22  *
23  * chi2 = <R|1./cov^{2} | R>
24  *
25  * dchi2 = <dR | 1./cov^{2} | R> + <R | 1./cov^{2} | dR>
26  */
27 namespace {
28  constexpr double c_inv = 1. / Gaudi::Units::c_light;
29  /* Cut off value for the determinant. Hessian matrices with a determinant smaller than this
30  are considered to be invalid */
31  constexpr double detCutOff = 1.e-8;
32 }
33 
34 namespace MuonR4{
35  using namespace SegmentFit;
36  using HitType = SegmentFitResult::HitType;
37  using HitVec = SegmentFitResult::HitVec;
38 
39  MdtSegmentFitter::Config::RangeArray
40  MdtSegmentFitter::Config::defaultRanges() {
41  RangeArray rng{};
42  constexpr double spatRang = 10.*Gaudi::Units::m;
43  constexpr double timeTange = 25 * Gaudi::Units::ns;
44  rng[toInt(ParamDefs::y0)] = std::array{-spatRang, spatRang};
45  rng[toInt(ParamDefs::x0)] = std::array{-spatRang, spatRang};
46  rng[toInt(ParamDefs::phi)] = std::array{-179.* Gaudi::Units::deg, 179. * Gaudi::Units::deg};
47  rng[toInt(ParamDefs::theta)] = std::array{-85. * Gaudi::Units::deg, 85. * Gaudi::Units::deg};
48  rng[toInt(ParamDefs::time)] = std::array{-timeTange, timeTange};
49  return rng;
50  }
51  MdtSegmentFitter::MdtSegmentFitter(const std::string& name, Config&& config):
52  AthMessaging{name},
53  m_cfg{std::move(config)}{}
54 
55  inline void MdtSegmentFitter::updateDriftSigns(const Amg::Vector3D& segPos, const Amg::Vector3D& segDir,
56  SegmentFitResult& fitResult) const {
57  for (std::unique_ptr<CalibratedSpacePoint>& meas : fitResult.calibMeasurements) {
58  meas->setDriftRadius(SegmentFitHelpers::driftSign(segPos,segDir, *meas, msg())*meas->driftRadius());
59  }
60  }
61 
62  inline bool MdtSegmentFitter::recalibrate(const EventContext& ctx,
63  SegmentFitResult& fitResult) const{
64 
65  const auto [segPos, segDir] = makeLine(fitResult.segmentPars);
66 
67  fitResult.calibMeasurements = m_cfg.calibrator->calibrate(ctx, std::move(fitResult.calibMeasurements), segPos, segDir,
68  fitResult.segmentPars[toInt(ParamDefs::time)]);
69  if (!updateHitSummary(fitResult)){
70  return false;
71  }
72  updateDriftSigns(segPos, segDir, fitResult);
74  if (fitResult.timeFit && fitResult.nDoF <= 1) {
75  fitResult.timeFit = false;
76  ATH_MSG_DEBUG("Switch of the time fit because nDoF: "<<fitResult.nDoF);
77  fitResult.segmentPars[toInt(ParamDefs::time)] = 0.;
79  fitResult.calibMeasurements = m_cfg.calibrator->calibrate(ctx, std::move(fitResult.calibMeasurements), segPos, segDir,
80  fitResult.segmentPars[toInt(ParamDefs::time)]);
82  } else if (!fitResult.timeFit && m_cfg.doTimeFit) {
83  ATH_MSG_DEBUG("Somehow a measurement is on the narrow ridge of validity. Let's try if the time can be fitted now ");
84  fitResult.timeFit = true;
85  }
86 
87  return true;
88  }
89  inline bool MdtSegmentFitter::updateHitSummary(SegmentFitResult& fitResult) const {
91  using State = CalibratedSpacePoint::State;
92  fitResult.nPhiMeas = fitResult.nDoF = fitResult.nTimeMeas = 0;
93  for (const HitType& hit : fitResult.calibMeasurements) {
94  if (hit->fitState() != State::Valid){
95  continue;
96  }
97 
98  fitResult.nPhiMeas+= hit->measuresPhi();
99  fitResult.nDoF+= hit->measuresPhi();
100  fitResult.nDoF+= hit->measuresEta();
101  fitResult.nPrecMeas+= (hit->type() == xAOD::UncalibMeasType::MdtDriftCircleType);
103  fitResult.nDoF += (m_cfg.doTimeFit && hit->type() != xAOD::UncalibMeasType::MdtDriftCircleType && hit->measuresTime());
104  fitResult.nTimeMeas+=hit->measuresTime();
105  }
106  if (!fitResult.nDoF) {
107  ATH_MSG_VERBOSE("Measurements rejected.");
108  return false;
109  }
110  if (!fitResult.nPhiMeas) {
111  ATH_MSG_VERBOSE("No phi measurements are left.");
112  fitResult.segmentPars[toInt(ParamDefs::phi)] = 90. * Gaudi::Units::deg;
113  const double nHits = fitResult.calibMeasurements.size();
114  double avgX{0.};
115  std::ranges::for_each(fitResult.calibMeasurements,[&avgX, nHits](const HitType& hit) {
116  return avgX += hit->positionInChamber().x() / nHits;
117  });
118  fitResult.segmentPars[toInt(ParamDefs::x0)] = avgX;
119  }
120 
121  fitResult.nDoF = fitResult.nDoF - 2 - (fitResult.nPhiMeas > 0 ? 2 : 0);
122 
123  return true;
124  }
125  void MdtSegmentFitter::updateLinePartials(const Parameters& params,
126  LineWithPartials& line) {
127 
128  line.pos[Amg::x] = params[toInt(ParamDefs::x0)];
129  line.pos[Amg::y] = params[toInt(ParamDefs::y0)];
130  line.dir = Amg::dirFromAngles(params[toInt(ParamDefs::phi)],
131  params[toInt(ParamDefs::theta)]);
132  line.gradient[toInt(ParamDefs::x0)] = Amg::Vector3D::UnitX();
133  line.gradient[toInt(ParamDefs::y0)] = Amg::Vector3D::UnitY();
144  const CxxUtils::sincos theta{params[toInt(ParamDefs::theta)]},
145  phi{params[toInt(ParamDefs::phi)]};
146  line.gradient[toInt(ParamDefs::theta)] = Amg::Vector3D{phi.cs*theta.cs, phi.sn*theta.cs, -theta.sn};
147  line.gradient[toInt(ParamDefs::phi)] = Amg::Vector3D{-theta.sn *phi.sn, theta.sn*phi.cs, 0};
148  /*********************************************************************************
149  * Non-vanishing second order derivatives
150  *
151  * d^{2} segDir d^{2} segDir cos phi
152  * ------------- = - segDir , ------------ = - sin(theta) sin phi
153  * d^{2} theta d^{2} phi 0
154  *
155  * d^{2} segDir -sin phi
156  * ------------- = cos(theta) cos phi
157  * d theta dPhi 0
158  ************************************************************************************/
159  constexpr unsigned nPars = LineWithPartials::nPars;
160  constexpr int idxThetaSq = vecIdxFromSymMat<nPars>(toInt(ParamDefs::theta), toInt(ParamDefs::theta));
161  constexpr int idxPhiSq = vecIdxFromSymMat<nPars>(toInt(ParamDefs::phi), toInt(ParamDefs::phi));
162  constexpr int idxPhiTheta = vecIdxFromSymMat<nPars>(toInt(ParamDefs::theta), toInt(ParamDefs::phi));
163  line.hessian[idxThetaSq] = - Amg::Vector3D{phi.cs*theta.sn,phi.sn*theta.sn, theta.cs};
164  line.hessian[idxPhiSq] = -theta.sn * Amg::Vector3D{phi.cs, phi.sn, 0.};
165  line.hessian[idxPhiTheta] = theta.cs * Amg::Vector3D{-phi.sn, phi.cs, 0.};
166  }
167  void MdtSegmentFitter::calculateWireResiduals(const LineWithPartials& line,
168  const CalibratedSpacePoint& hit,
169  ResidualWithPartials& resObj) const {
171  const Amg::Vector3D& hitPos{hit.positionInChamber()};
172  const Amg::Vector3D& hitDir{hit.directionInChamber()};
174  resObj.projIntoWirePlane = line.dir.dot(hitDir);
175  resObj.projDirLenSq = 1. - resObj.projIntoWirePlane*resObj.projIntoWirePlane;
176  if (resObj.projDirLenSq < std::numeric_limits<float>::epsilon()) {
177  resObj.residual.setZero();
178  for (Amg::Vector3D& grad : resObj.gradient){
179  grad.setZero();
180  }
181  if (m_cfg.useSecOrderDeriv) {
182  for (Amg::Vector3D& hess : resObj.hessian) {
183  hess.setZero();
184  }
185  }
186  ATH_MSG_WARNING("Segment is parallel along the wire");
187  return;
188  }
189  resObj.invProjLenSq = 1. / resObj.projDirLenSq;
190  resObj.invProjLen = std::sqrt(resObj.invProjLenSq);
192  resObj.projDir = (line.dir - resObj.projIntoWirePlane*hitDir) * resObj.invProjLen;
194  const Amg::Vector3D hitMinSeg = hitPos - line.pos ;
196  const double lineDist = resObj.projDir.cross(hitDir).dot(hitMinSeg);
197  const double resVal = (lineDist - hit.driftRadius());
198  resObj.residual = resVal * Amg::Vector3D::Unit(toInt(AxisDefs::eta));
199  ATH_MSG_VERBOSE("Mdt drift radius: "<<hit.driftRadius()<<" distance: "<<lineDist<<";"
200  <<Amg::signedDistance(hitPos, hitDir, line.pos, line.dir)<<", residual: "<<resVal);
203  if (hit.dimension() == 2) {
204  resObj.residual[toInt(AxisDefs::phi)] = (hitMinSeg.dot(line.dir)*resObj.projIntoWirePlane - hitMinSeg.dot(hitDir)) * resObj.invProjLenSq;
205  }
206  constexpr unsigned nLinePars = LineWithPartials::nPars;
208  for (const int param : {toInt(ParamDefs::y0), toInt(ParamDefs::x0),
209  toInt(ParamDefs::theta), toInt(ParamDefs::phi)}) {
210  if (!resObj.evalPhiPars && (param == toInt(ParamDefs::x0) || param == toInt(ParamDefs::phi))){
211  continue;
212  }
213  switch (param) {
214  case toInt(ParamDefs::theta):
215  case toInt(ParamDefs::phi): {
216  resObj.partWirePlaneProj[param] = line.gradient[param].dot(hitDir);
217  resObj.partProjDir[param] = (line.gradient[param] - resObj.partWirePlaneProj[param]*hitDir) * resObj.invProjLen
218  + resObj.partWirePlaneProj[param]*resObj.projIntoWirePlane* resObj.projDir * resObj.invProjLenSq;
219  const double partialDist = resObj.partProjDir[param].cross(hitDir).dot(hitMinSeg);
220 
221  resObj.gradient[param] = partialDist * Amg::Vector3D::Unit(toInt(AxisDefs::eta));
222  if (hit.dimension() == 2) {
223  resObj.gradient[param][toInt(AxisDefs::phi)] =
224  ( hitMinSeg.dot(line.gradient[param]) * resObj.projIntoWirePlane +
225  hitMinSeg.dot(line.dir)*resObj.partWirePlaneProj[param]) * resObj.invProjLenSq
226  +2.* resObj.residual[toInt(AxisDefs::phi)]*( resObj.projIntoWirePlane * resObj.partWirePlaneProj[param]) * resObj.invProjLenSq;
227  }
228  break;
229  } case toInt(ParamDefs::y0):
230  case toInt(ParamDefs::x0): {
231  const double partialDist = - resObj.projDir.cross(hitDir).dot(line.gradient[param]);
232  resObj.gradient[param] = partialDist * Amg::Vector3D::Unit(toInt(AxisDefs::eta));
233  if (hit.dimension() == 2) {
234  resObj.gradient[param][toInt(AxisDefs::phi)] = -(line.gradient[param].dot(line.dir) * resObj.projIntoWirePlane -
235  line.gradient[param].dot(hitDir)) * resObj.invProjLenSq;
236  }
237  break;
238  }
240  default:
241  break;
242  }
243  }
244  if (!m_cfg.useSecOrderDeriv) {
245  return;
246  }
248  for (int param = toInt(ParamDefs::phi); param >=0; --param){
249  if (!resObj.evalPhiPars && (param == toInt(ParamDefs::x0) || param == toInt(ParamDefs::phi))){
250  continue;
251  }
252  for (int param1 = param; param1>=0; --param1) {
253  if (!resObj.evalPhiPars && (param1 == toInt(ParamDefs::x0) || param1 == toInt(ParamDefs::phi))){
254  continue;
255  }
256  const int lineIdx = vecIdxFromSymMat<nLinePars>(param, param1);
257  const int resIdx = vecIdxFromSymMat<toInt(ParamDefs::nPars)>(param, param1);
259  if ( (param == toInt(ParamDefs::theta) || param == toInt(ParamDefs::phi)) &&
260  (param1 == toInt(ParamDefs::theta) || param1 == toInt(ParamDefs::phi))) {
261 
262  const double partSqLineProject = line.hessian[lineIdx].dot(hitDir);
263  const Amg::Vector3D projDirPartSq = (line.hessian[lineIdx] - partSqLineProject * hitDir) * resObj.invProjLen
264  + (resObj.partWirePlaneProj[param1] * resObj.projIntoWirePlane) * resObj.invProjLenSq * resObj.partProjDir[param]
265  + (resObj.partWirePlaneProj[param] * resObj.projIntoWirePlane) * resObj.invProjLenSq * resObj.partProjDir[param1]
266  + (partSqLineProject*resObj.projIntoWirePlane) * resObj.invProjLenSq * resObj.projDir
267  + (resObj.partWirePlaneProj[param1] * resObj.partWirePlaneProj[param]) * std::pow(resObj.invProjLenSq, 2) * resObj.projDir;
268 
269  const double partialSqDist = projDirPartSq.cross(hitDir).dot(hitMinSeg);
270  resObj.hessian[resIdx] = partialSqDist * Amg::Vector3D::Unit(toInt(AxisDefs::eta));
271  if (hit.dimension() == 2) {
272  const double partialSqAlongWire =2.*resObj.residual[toInt(AxisDefs::phi)]*resObj.projIntoWirePlane*partSqLineProject * resObj.invProjLenSq
273  +2.*resObj.residual[toInt(AxisDefs::phi)]*resObj.partWirePlaneProj[param]*resObj.partWirePlaneProj[param1]*resObj.invProjLenSq
274  +2.*resObj.gradient[param1][toInt(AxisDefs::phi)]*resObj.projIntoWirePlane*resObj.partWirePlaneProj[param]*resObj.invProjLenSq
275  +2.*resObj.gradient[param][toInt(AxisDefs::phi)]*resObj.projIntoWirePlane*resObj.partWirePlaneProj[param1]*resObj.invProjLenSq
276  + hitMinSeg.dot(line.hessian[lineIdx]) *resObj.projIntoWirePlane * resObj.invProjLenSq
277  + hitMinSeg.dot(line.dir)*partSqLineProject * resObj.invProjLenSq
278  + hitMinSeg.dot(line.gradient[param])*resObj.partWirePlaneProj[param1]*resObj.invProjLenSq
279  + hitMinSeg.dot(line.gradient[param1])*resObj.partWirePlaneProj[param]*resObj.invProjLenSq;
280  resObj.hessian[resIdx][toInt(AxisDefs::phi)] = partialSqAlongWire;
281  }
282  }
284  else if (param == toInt(ParamDefs::theta) || param == toInt(ParamDefs::phi)){
285  const double partialSqDist = - resObj.partProjDir[param].cross(hitDir).dot(line.gradient[param1]);
286  resObj.hessian[resIdx] = partialSqDist * Amg::Vector3D::Unit(toInt(AxisDefs::eta));
287  if (hit.dimension() == 2) {
288  const double partialSqAlongWire = -(line.gradient[param1].dot(line.gradient[param])*resObj.projIntoWirePlane +
289  line.gradient[param1].dot(line.dir)*resObj.partWirePlaneProj[param]) * resObj.invProjLenSq
290  + 2.* resObj.gradient[param1][toInt(AxisDefs::phi)]*( resObj.projIntoWirePlane * resObj.partWirePlaneProj[param]) * resObj.invProjLenSq;
291  resObj.hessian[resIdx][toInt(AxisDefs::phi)] = partialSqAlongWire;
292  }
293  }
294  }
295  }
296  }
297  void MdtSegmentFitter::calculateStripResiduals(const LineWithPartials& line,
298  const CalibratedSpacePoint& hit,
299  ResidualWithPartials& residual) const {
300 
301 
302  const Amg::Vector3D& hitPos{hit.positionInChamber()};
303  const Amg::Vector3D normal = hit.type() != xAOD::UncalibMeasType::Other
304  ? hit.spacePoint()->planeNormal() : Amg::Vector3D::UnitZ();
305  const double planeOffSet = normal.dot(hitPos);
306 
307  const double normDot = normal.dot(line.dir);
308  if (std::abs(normDot) < std::numeric_limits<double>::epsilon()){
309  residual.residual.setZero();
310  for (Amg::Vector3D& deriv: residual.gradient) {
311  deriv.setZero();
312  }
313  if (m_cfg.useSecOrderDeriv) {
314  for (Amg::Vector3D& hessian : residual.hessian){
315  hessian.setZero();
316  }
317  }
318  ATH_MSG_WARNING("The hit parallel with the segment line "<<Amg::toString(line.dir));
319  return;
320  }
321  const double travelledDist = (planeOffSet - line.pos.dot(normal)) / normDot;
322 
323  const Amg::Vector3D planeIsect = line.pos + travelledDist * line.dir;
325  residual.residual.block<2,1>(0,0) = (planeIsect - hitPos).block<2,1>(0,0);
326 
327  for (unsigned fitPar = 0; fitPar < line.gradient.size(); ++fitPar) {
328  switch (fitPar) {
329  case toInt(ParamDefs::phi):
330  case toInt(ParamDefs::theta):{
331  const double partialDist = - travelledDist / normDot * normal.dot(line.gradient[fitPar]);
332  residual.gradient[fitPar].block<2,1>(0,0) = (travelledDist * line.gradient[fitPar] + partialDist * line.dir).block<2,1>(0,0);
333  break;
334  } case toInt(ParamDefs::y0):
335  case toInt(ParamDefs::x0): {
336  residual.gradient[fitPar].block<2,1>(0,0) = (line.gradient[fitPar] - line.gradient[fitPar].dot(normal) / normDot * line.dir).block<2,1>(0,0);
337  break;
338  }
339  default: {
340  break;
341  }
342  }
343  }
344  if (!m_cfg.useSecOrderDeriv) {
345  return;
346  }
347  constexpr unsigned nLinePars = LineWithPartials::nPars;
348  for (int param = toInt(ParamDefs::phi); param >=0; --param){
349  for (int param1 = param; param1>=0; --param1) {
350  const int lineIdx = vecIdxFromSymMat<nLinePars>(param, param1);
351  const int resIdx = vecIdxFromSymMat<toInt(ParamDefs::nPars)>(param, param1);
352  if ( (param == toInt(ParamDefs::theta) || param == toInt(ParamDefs::phi)) &&
353  (param1 == toInt(ParamDefs::theta) || param1 == toInt(ParamDefs::phi))) {
354  residual.hessian[resIdx].block<2,1>(0,0) =(
355  travelledDist * line.hessian[lineIdx] - travelledDist *(normal.dot(line.hessian[lineIdx])) / normDot * line.dir
356  - normal.dot(line.gradient[param1]) / normDot * residual.gradient[param]
357  - normal.dot(line.gradient[param]) / normDot * residual.gradient[param1]).block<2,1>(0,0);
358  } else if (param == toInt(ParamDefs::theta) || param == toInt(ParamDefs::phi)) {
359  const double gradientDisplace = normal.dot(line.gradient[param1]);
360  if (gradientDisplace > std::numeric_limits<float>::epsilon()){
361  residual.hessian[resIdx].block<2,1>(0,0) = gradientDisplace*( normal.dot(line.gradient[param]) / std::pow(normDot,2)*line.dir
362  - line.gradient[param] / normDot ).block<2,1>(0,0);
363  } else {
364  residual.hessian[resIdx].setZero();
365  }
366  }
367  }
368  }
369  }
370 
371  SegmentFitResult MdtSegmentFitter::fitSegment(const EventContext& ctx,
372  HitVec&& calibHits,
373  const Parameters& startPars,
374  const Amg::Transform3D& localToGlobal) const {
375 
376  const Muon::IMuonIdHelperSvc* idHelperSvc{calibHits[0]->spacePoint()->msSector()->idHelperSvc()};
377  using State = CalibratedSpacePoint::State;
378 
379  if (msgLvl(MSG::VERBOSE)) {
380  std::stringstream hitStream{};
381  const auto [startPos, startDir] = makeLine(startPars);
382  for (const HitType& hit : calibHits) {
383  hitStream<<" **** "<<(hit->type() != xAOD::UncalibMeasType::Other ? idHelperSvc->toString(hit->spacePoint()->identify()): "beamspot" )
384  <<" position: "<<Amg::toString(hit->positionInChamber());
385  if (hit->type() == xAOD::UncalibMeasType::MdtDriftCircleType) {
386  hitStream<<", driftRadius: "<<SegmentFitHelpers::driftSign(startPos,startDir, *hit, msg())*hit->driftRadius() ;
387  }
388  hitStream<<", channel dir: "<<Amg::toString(hit->directionInChamber())<<std::endl;
389  }
390  ATH_MSG_VERBOSE("Start segment fit with parameters "<<toString(startPars)
391  <<", plane location: "<<Amg::toString(localToGlobal)<<std::endl<<hitStream.str());
392 
393  }
394 
395  SegmentFitResult fitResult{};
396  fitResult.segmentPars = startPars;
397  fitResult.timeFit = m_cfg.doTimeFit;
398  fitResult.calibMeasurements = std::move(calibHits);
399  if (!updateHitSummary(fitResult)) {
400  ATH_MSG_WARNING("No valid segment seed parsed from the beginning");
401  return fitResult;
402  }
403  if (!m_cfg.reCalibrate) {
404  const auto [segPos, segDir] = makeLine(startPars);
405  updateDriftSigns(segPos, segDir, fitResult);
406  }
407  Parameters gradient{AmgVector(5)::Zero()}, prevGrad{AmgVector(5)::Zero()}, prevPars{AmgVector(5)::Zero()};
408  AmgSymMatrix(5) hessian{AmgSymMatrix(5)::Zero()};
409 
411  LineWithPartials segmentLine{};
413  ResidualWithPartials residual{};
414 
415  unsigned int noChangeIter{0};
416  while (fitResult.nIter++ < m_cfg.nMaxCalls) {
417  ATH_MSG_DEBUG("Iteration: "<<fitResult.nIter<<" parameters: "<<toString(fitResult.segmentPars)<<", chi2: "<<fitResult.chi2);
419  updateLinePartials(fitResult.segmentPars, segmentLine);
420 
422  if (m_cfg.reCalibrate && !recalibrate(ctx, fitResult)) {
423  break;
424  }
426  fitResult.chi2 = 0;
427  hessian.setZero();
428  gradient.setZero();
429 
431  for (HitType& hit : fitResult.calibMeasurements) {
432  if (hit->fitState() != State::Valid) {
433  continue;
434  }
435  const Amg::Vector3D& hitPos{hit->positionInChamber()};
436  const Amg::Vector3D& hitDir{hit->directionInChamber()};
437  ATH_MSG_DEBUG("Update chi2 from measurement "<<(hit->spacePoint() ? idHelperSvc->toString(hit->spacePoint()->identify())
438  : "pseudo meas")<<" position: "<<Amg::toString(hitPos)<<" + "<<Amg::toString(hitDir));
439 
440  const int start = toInt(fitResult.nPhiMeas ? ParamDefs::phi : ParamDefs::theta);
441  switch (hit->type()) {
442  case xAOD::UncalibMeasType::MdtDriftCircleType: {
445  const double dSign = (hit->driftRadius() > 0 ? 1. : -1.);
446  if (fitResult.timeFit && !m_cfg.reCalibrate && fitResult.nIter > 1) {
447  hit = m_cfg.calibrator->calibrate(ctx, *hit, segmentLine.pos, segmentLine.dir,
448  fitResult.segmentPars[toInt(ParamDefs::time)]);
449  hit->setDriftRadius(dSign*hit->driftRadius());
450  }
451  calculateWireResiduals(segmentLine, *hit, residual);
452 
453  if (!fitResult.timeFit) {
454  break;
455  }
457  residual.gradient[toInt(ParamDefs::time)] = dSign*m_cfg.calibrator->driftVelocity(ctx, *hit) * Amg::Vector3D::Unit(toInt(AxisDefs::eta));
458  if (m_cfg.useSecOrderDeriv) {
459  constexpr unsigned hessIdx = vecIdxFromSymMat<ResidualWithPartials::nPars>(toInt(ParamDefs::time), toInt(ParamDefs::time));
460  residual.hessian[hessIdx] = dSign*m_cfg.calibrator->driftAcceleration(ctx, *hit) * Amg::Vector3D::Unit(toInt(AxisDefs::eta));
461  }
462  break;
463  }
464  case xAOD::UncalibMeasType::RpcStripType:
465  case xAOD::UncalibMeasType::TgcStripType:
466  case xAOD::UncalibMeasType::sTgcStripType:
467  case xAOD::UncalibMeasType::MMClusterType: {
468  calculateStripResiduals(segmentLine, *hit, residual);
469  if (!fitResult.timeFit || !hit->measuresTime()) {
470  break;
471  }
472  constexpr ParamDefs par = ParamDefs::time;
473  residual.gradient[toInt(par)] = -Amg::Vector3D::Unit(toInt(AxisDefs::t0));
474  const Amg::Vector3D planeIsect = residual.residual + hitPos;
475  const double totFlightDist = (localToGlobal * planeIsect).mag() * c_inv;
476  residual.residual[toInt(AxisDefs::t0)] = hit->time() - totFlightDist * c_inv - fitResult.segmentPars[toInt(ParamDefs::time)];
477 
478  for (int p = start; p >= 0; --p) {
479  residual.gradient[p][toInt(AxisDefs::t0)] = -residual.gradient[p].perp() * c_inv;
480  ATH_MSG_VERBOSE("Partial derivative of "<<idHelperSvc->toString(hit->spacePoint()->identify())
481  <<" residual "<<Amg::toString(residual.residual)<<" "<<Amg::toString(multiply(inverse(hit->covariance()), residual.residual))
482  <<" "<<std::endl<<toString(hit->covariance())<<std::endl<<" w.r.t "
483  <<toString(static_cast<ParamDefs>(p))<<"="<<Amg::toString(residual.gradient[p]));
484  }
485  break;
486  } case xAOD::UncalibMeasType::Other: {
487  calculateStripResiduals(segmentLine, *hit, residual);
488  break;
489  } default: {
490  const auto &measurement = *hit->spacePoint()->primaryMeasurement();
491  ATH_MSG_WARNING("MdtSegmentFitter() - Unsupported measurment type" <<typeid(measurement).name());
492  }
493  }
494 
495  ATH_MSG_VERBOSE("Update derivatives for hit "<< (hit->spacePoint() ? idHelperSvc->toString(hit->spacePoint()->identify()) : "beamspot"));
496  updateDerivatives(residual, hit->covariance(), gradient, hessian, fitResult.chi2,
497  fitResult.timeFit && hit->measuresTime() ? toInt(ParamDefs::time) : start);
498  }
499 
501  for (int k =1; k < 5 - (!fitResult.timeFit); ++k){
502  for (int l = 0; l< k; ++l){
503  hessian(l,k) = hessian(k,l);
504  }
505  }
507  if (gradient.mag() < m_cfg.tolerance) {
508  fitResult.converged = true;
509  ATH_MSG_VERBOSE("Fit converged after "<<fitResult.nIter<<" iterations with "<<fitResult.chi2);
510  break;
511  }
512  ATH_MSG_VERBOSE("Chi2: "<<fitResult.chi2<<", gradient: "<<Amg::toString(gradient)<<"hessian: "<<std::endl<<hessian);
514  UpdateStatus paramUpdate{UpdateStatus::outOfBounds};
515  if (!fitResult.nPhiMeas && !fitResult.timeFit) {
516  paramUpdate = updateParameters<2>(fitResult.segmentPars, prevPars, gradient, prevGrad, hessian);
517  } else if (!fitResult.nPhiMeas && fitResult.timeFit) {
520  hessian.col(2).swap(hessian.col(toInt(ParamDefs::time)));
521  hessian.row(2).swap(hessian.row(toInt(ParamDefs::time)));
522  std::swap(gradient[2], gradient[toInt(ParamDefs::time)]);
523  std::swap(fitResult.segmentPars[2], fitResult.segmentPars[toInt(ParamDefs::time)]);
524 
525  paramUpdate = updateParameters<3>(fitResult.segmentPars, prevPars, gradient, prevGrad, hessian);
526 
527  std::swap(fitResult.segmentPars[2], fitResult.segmentPars[toInt(ParamDefs::time)]);
528  hessian.col(2).swap(hessian.col(toInt(ParamDefs::time)));
529  hessian.row(2).swap(hessian.row(toInt(ParamDefs::time)));
530  std::swap(gradient[2], gradient[toInt(ParamDefs::time)]);
531  } else if (fitResult.nPhiMeas && !fitResult.timeFit) {
532  paramUpdate = updateParameters<4>(fitResult.segmentPars, prevPars, gradient, prevGrad, hessian);
533  } else if (fitResult.nPhiMeas && fitResult.timeFit) {
534  paramUpdate = updateParameters<5>(fitResult.segmentPars, prevPars, gradient, prevGrad, hessian);
535  }
536  if (paramUpdate == UpdateStatus::noChange) {
537  if ((++noChangeIter) >= m_cfg.noMoveIter){
538  fitResult.converged = true;
539  break;
540  }
541  } else if (paramUpdate == UpdateStatus::allOkay) {
542  noChangeIter = 0;
543  } else if (paramUpdate == UpdateStatus::outOfBounds){
544  return fitResult;
545  }
546  }
547 
549  fitResult.nDoF-=fitResult.timeFit;
550 
552  const auto [segPos, segDir] = fitResult.makeLine();
553  fitResult.chi2 =0.;
554 
555  std::optional<double> toF = fitResult.timeFit ? std::make_optional<double>((localToGlobal * segPos).mag() * c_inv) : std::nullopt;
557  std::ranges::stable_sort(fitResult.calibMeasurements, [](const HitType&a, const HitType& b){
558  return a->positionInChamber().z() < b->positionInChamber().z();
559  });
560 
561  /*** Remove the drift sign again */
562  fitResult.chi2 =0.;
563  for (const HitType& hit : fitResult.calibMeasurements) {
564  hit->setDriftRadius(std::abs(hit->driftRadius()));
565  fitResult.chi2 +=SegmentFitHelpers::chiSqTerm(segPos, segDir, fitResult.segmentPars[toInt(ParamDefs::time)], toF, *hit, msg());
566  }
568  if (!fitResult.nPhiMeas&& !fitResult.timeFit) {
569  blockCovariance<2>(std::move(hessian), /* std::move(gradient), std::move(prevGrad), */ fitResult.segmentParErrs);
570  }else if (!fitResult.nPhiMeas && fitResult.timeFit) {
571  hessian.col(2).swap(hessian.col(toInt(ParamDefs::time)));
572  hessian.row(2).swap(hessian.row(toInt(ParamDefs::time)));
573  blockCovariance<3>(std::move(hessian),/* std::move(gradient), std::move(prevGrad), */ fitResult.segmentParErrs);
574  fitResult.segmentParErrs.col(2).swap(fitResult.segmentParErrs.col(toInt(ParamDefs::time)));
575  fitResult.segmentParErrs.row(2).swap(fitResult.segmentParErrs.row(toInt(ParamDefs::time)));
576  } else if (fitResult.nPhiMeas) {
577  blockCovariance<4>(std::move(hessian), /* std::move(gradient), std::move(prevGrad), */ fitResult.segmentParErrs);
578  } else if (fitResult.nPhiMeas && fitResult.timeFit) {
579  blockCovariance<5>(std::move(hessian),/* std::move(gradient), std::move(prevGrad), */ fitResult.segmentParErrs);
580  }
581  return fitResult;
582  }
583 
584  void MdtSegmentFitter::updateDerivatives(const ResidualWithPartials& fitMeas,
585  const MeasCov_t& covariance,
586  AmgVector(5)& gradient, AmgSymMatrix(5)& hessian,
587  double& chi2, int startPar) const {
588  const MeasCov_t invCov{inverse(covariance)};
589  const Amg::Vector3D covRes = multiply(invCov, fitMeas.residual);
590  chi2 += covRes.dot(fitMeas.residual);
591  for (int p = startPar; p >=0 ; --p) {
592  gradient[p] +=2.*covRes.dot(fitMeas.gradient[p]);
593  for (int k=p; k>=0; --k) {
594  hessian(p,k)+= 2.*contract(invCov, fitMeas.gradient[p], fitMeas.gradient[k]);
595  if (m_cfg.useSecOrderDeriv) {
596  const int symMatIdx = vecIdxFromSymMat<ResidualWithPartials::nPars>(p,k);
597  hessian(p,k)+=contract(invCov, covRes, fitMeas.hessian[symMatIdx]);
598  }
599  }
600  }
601  ATH_MSG_VERBOSE("After derivative update --- chi2: "<<chi2<<"("<<covRes.dot(fitMeas.residual)<<"), gradient: "
602  <<toString(gradient)<<", Hessian:\n"<<hessian<<", measurement covariance\n"<<toString(invCov));
603  }
604 
605  template <unsigned int nDim>
606  void MdtSegmentFitter::blockCovariance(const AmgSymMatrix(5)& hessian,
607  SegmentFit::Covariance& covariance) const {
608 
609  covariance.setIdentity();
610  AmgSymMatrix(nDim) miniHessian = hessian.block<nDim, nDim>(0,0);
611  if (std::abs(miniHessian.determinant()) <= detCutOff) {
612  ATH_MSG_VERBOSE("Boeser mini hessian ("<<miniHessian.determinant()<<")\n"<<miniHessian<<"\n\n"<<hessian);
613  return;
614  }
615  ATH_MSG_VERBOSE("Hessian matrix: \n"<<hessian<<",\nblock Hessian:\n"<<miniHessian<<",\n determinant: "<<miniHessian.determinant());
616  covariance.block<nDim,nDim>(0,0) = miniHessian.inverse();
617  ATH_MSG_VERBOSE("covariance: \n"<<covariance);
618  }
619 
620  template <unsigned int nDim>
621  MdtSegmentFitter::UpdateStatus
622  MdtSegmentFitter::updateParameters(Parameters& currPars, Parameters& prevPars,
623  Parameters& currGrad, Parameters& prevGrad,
624  const AmgSymMatrix(5)& currentHessian) const {
625 
626  AmgSymMatrix(nDim) miniHessian = currentHessian.block<nDim, nDim>(0,0);
627  ATH_MSG_DEBUG("Parameter update -- \ncurrenPars: "<<toString(currPars)<<", \ngradient: "<<toString(currGrad)
628  <<", hessian ("<<miniHessian.determinant()<<")"<<std::endl<<miniHessian);
629 
630  double changeMag{0.};
631  if (std::abs(miniHessian.determinant()) > detCutOff) {
632  prevPars.block<nDim,1>(0,0) = currPars.block<nDim,1>(0,0);
633  // Update the parameters accrodingly to the hessian
634  const AmgVector(nDim) updateMe = miniHessian.inverse()* currGrad.block<nDim, 1>(0,0);
635  changeMag = std::sqrt(updateMe.dot(updateMe));
636  currPars.block<nDim,1>(0,0) -= updateMe;
637  prevGrad.block<nDim,1>(0,0) = currGrad.block<nDim,1>(0,0);
638  ATH_MSG_DEBUG("Hessian inverse:\n"<<miniHessian.inverse()<<"\n\nUpdate the parameters by -"
639  <<Amg::toString(updateMe));
640  } else {
641  const AmgVector(nDim) gradDiff = (currGrad - prevGrad).block<nDim,1>(0,0);
642  const double gradDiffMag = gradDiff.mag2();
643  double denom = (gradDiffMag > std::numeric_limits<float>::epsilon() ? gradDiffMag : 1.);
644  const double gamma = std::abs((currPars - prevPars).block<nDim,1>(0,0).dot(gradDiff))
645  / denom;
646  ATH_MSG_VERBOSE("Hessian determinant invalid. Try deepest descent - \nprev parameters: "
647  <<toString(prevPars)<<",\nprevious gradient: "<<toString(prevGrad)<<", gamma: "<<gamma);
648  prevPars.block<nDim, 1>(0,0) = currPars.block<nDim, 1>(0,0);
649  currPars.block<nDim, 1>(0,0) -= gamma* currGrad.block<nDim, 1>(0,0);
650  prevGrad.block<nDim, 1>(0,0) = currGrad.block<nDim,1>(0,0);
651  changeMag = std::abs(gamma) * currGrad.block<nDim, 1>(0,0).mag();
652  }
654  unsigned int nOutOfBound{0};
655  for (unsigned int p = 0; p< nDim; ++p) {
656  double& parValue{currPars[p]};
657  if constexpr(nDim == 3) {
659  if (p == 2) {
660  p = toInt(ParamDefs::time);
661  }
662  }
663  if (m_cfg.ranges[p][0] > parValue || m_cfg.ranges[p][1] < parValue) {
664  ATH_MSG_VERBOSE("The "<<p<<"-th parameter "<<toString(static_cast<ParamDefs>(p))<<" is out of range "<<parValue
665  <<" allowed ["<<m_cfg.ranges[p][0]<<"-"<<m_cfg.ranges[p][1]<<"]");
666  ++nOutOfBound;
667  parValue = std::clamp(parValue, m_cfg.ranges[p][0], m_cfg.ranges[p][1]);
668  }
669 
670  }
671  if (nOutOfBound > m_cfg.nParsOutOfBounds){
672  return UpdateStatus::outOfBounds;
673  }
674  if (changeMag <= m_cfg.tolerance) {
675  return UpdateStatus::noChange;
676  }
677  return UpdateStatus::allOkay;
678  }
679 }
MuonR4::MdtSegmentFitter::ResidualWithPartials::residual
Amg::Vector3D residual
Vector carrying the residual.
Definition: MdtSegmentFitter.h:91
MuonR4::MdtSegmentFitter::ResidualAuxillaries::invProjLen
double invProjLen
inverse of the unormalized dir porjection
Definition: MdtSegmentFitter.h:83
MuonR4::MdtSegmentFitter::updateDriftSigns
void updateDriftSigns(const Amg::Vector3D &segPos, const Amg::Vector3D &segDir, SegmentFitResult &fitRes) const
Update the signs of the measurement.
Definition: MdtSegmentFitter.cxx:55
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:606
MuonR4::MdtSegmentFitter::updateHitSummary
bool updateHitSummary(SegmentFitResult &fitResult) const
Brief updates the hit summary from the contributing hits.
Definition: MdtSegmentFitter.cxx:89
MuonR4::MdtSegmentFitter::calculateStripResiduals
void calculateStripResiduals(const LineWithPartials &line, const CalibratedSpacePoint &spacePoint, ResidualWithPartials &residual) const
Calculates the residual together with hte correspdonding derivatives for strip measurements.
Definition: MdtSegmentFitter.cxx:297
MuonR4::MdtSegmentFitter::LineWithPartials::nPars
static constexpr unsigned nPars
Free parameters of the line (x0,y0,theta,phi)
Definition: MdtSegmentFitter.h:56
MuonR4::MdtSegmentFitter::ResidualWithPartials
Helper struct carrying the residual with its derivatives.
Definition: MdtSegmentFitter.h:87
MuonR4::MdtSegmentFitter::UpdateStatus
UpdateStatus
Definition: MdtSegmentFitter.h:172
MuonR4::MdtSegmentFitter::ResidualAuxillaries::partProjDir
std::array< Amg::Vector3D, nLinePars > partProjDir
Partial derivatives of the dir projection w.r.t.
Definition: MdtSegmentFitter.h:73
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:22
MuonR4::SegmentFitResult::calibMeasurements
HitVec calibMeasurements
Calibrated measurements used in the fit.
Definition: SegmentFitterEventData.h:64
mergePhysValFiles.start
start
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:13
MuonR4::MdtSegmentFitter::UpdateStatus::outOfBounds
@ outOfBounds
MuonR4::SpacePoint::planeNormal
Amg::Vector3D planeNormal() const
Returns the vector pointing out of the measurement plane.
Definition: MuonSpectrometer/MuonPhaseII/Event/MuonSpacePoint/src/SpacePoint.cxx:68
Amg::y
@ y
Definition: GeoPrimitives.h:35
xAOD::UncalibMeasType::MMClusterType
@ MMClusterType
MuonR4::MdtSegmentFitter::ResidualAuxillaries::partWirePlaneProj
std::array< double, nLinePars > partWirePlaneProj
Partial derivatives of the dir projection lengths w.r.t line parameters.
Definition: MdtSegmentFitter.h:75
sincos.h
Helper to simultaneously calculate sin and cos of the same angle.
deg
#define deg
Definition: SbPolyhedron.cxx:17
ISpacePointCalibrator.h
UploadAMITag.l
list l
Definition: UploadAMITag.larcaf.py:157
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 drift-circle space point.
xAOD::UncalibMeasType::sTgcStripType
@ sTgcStripType
MuonR4::MdtSegmentFitter::Config::noMoveIter
unsigned int noMoveIter
How many iterations with changes below tolerance.
Definition: MdtSegmentFitter.h:42
D3PDTest::rng
uint32_t rng()
Definition: FillerAlg.cxx:40
MuonR4::SegmentFitHelpers::chiSqTerm
double chiSqTerm(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 measurement.
Definition: SegmentFitHelperFunctions.cxx:28
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:26
ParamDefs.h
MuonR4::MdtSegmentFitter::Config::tolerance
double tolerance
Gradient cut off.
Definition: MdtSegmentFitter.h:30
MuonR4::toString
std::string toString(const CalibratedSpacePoint::Covariance_t &mat)
Returns the matrix in string.
Definition: MuonSpectrometer/MuonPhaseII/Event/MuonSpacePoint/src/UtilFunctions.cxx:75
dq_defect_bulk_create_defects.line
line
Definition: dq_defect_bulk_create_defects.py:27
MuonR4::MdtSegmentFitter::fitSegment
SegmentFitResult fitSegment(const EventContext &ctx, HitVec &&calibHits, const Parameters &startPars, const Amg::Transform3D &localToGlobal) const
Definition: MdtSegmentFitter.cxx:371
config
Definition: PhysicsAnalysis/AnalysisCommon/AssociationUtils/python/config.py:1
MuonR4::SegmentFitResult::HitVec
std::vector< HitType > HitVec
Definition: SegmentFitterEventData.h:53
MuonR4::SegmentFit::ParamDefs::phi
@ phi
xAOD::UncalibMeasType::TgcStripType
@ TgcStripType
MuonR4::MdtSegmentFitter::Config::calibrator
const ISpacePointCalibrator * calibrator
Pointer to the calibrator tool.
Definition: MdtSegmentFitter.h:40
MuonR4::MdtSegmentFitter::ResidualAuxillaries::projDirLenSq
double projDirLenSq
Length squared of the projected direction.
Definition: MdtSegmentFitter.h:79
MuonR4::MdtSegmentFitter::ResidualWithPartials::gradient
std::array< Amg::Vector3D, nPars > gradient
First order derivatives.
Definition: MdtSegmentFitter.h:93
MuonR4::SegmentFitResult::nDoF
int nDoF
degrees of freedom
Definition: SegmentFitterEventData.h:68
MuonR4::MdtSegmentFitter::Config::ranges
RangeArray ranges
Definition: MdtSegmentFitter.h:48
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:239
MuonR4::SegmentFit::AxisDefs::t0
@ t0
MuonR4::MdtSegmentFitter::ResidualWithPartials::evalPhiPars
bool evalPhiPars
Flag whether the the residuals w.r.t phi shall be evaluated.
Definition: MdtSegmentFitter.h:97
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
MuonR4::MdtSegmentFitter::ResidualAuxillaries::invProjLenSq
double invProjLenSq
inverse squared of the unnormalized dir projection
Definition: MdtSegmentFitter.h:81
Amg::toString
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Definition: GeoPrimitivesToStringConverter.h:40
MuonR4::ISpacePointCalibrator::driftAcceleration
virtual double driftAcceleration(const EventContext &ctx, const CalibratedSpacePoint &spacePoint) const =0
Returns the drift acceleration for a given drift-circle space point.
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:27
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:163
SegmentFitHelperFunctions.h
AmgVector
AmgVector(4) T2BSTrackFilterTool
Definition: T2BSTrackFilterTool.cxx:114
Amg::x
@ x
Definition: GeoPrimitives.h:34
MuonR4::MdtSegmentFitter::HitType
std::unique_ptr< CalibratedSpacePoint > HitType
Definition: MdtSegmentFitter.h:23
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition: GeoPrimitives.h:46
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:525
MuonR4::MdtSegmentFitter::Config
Definition: MdtSegmentFitter.h:26
MuonR4::CalibratedSpacePoint::State
State
State flag to distinguish different space point states.
Definition: CalibratedSpacePoint.h:24
MuonR4::SegmentFitResult::nPrecMeas
unsigned nPrecMeas
Number of precision hits.
Definition: SegmentFitterEventData.h:70
MuonR4::MdtSegmentFitter::Config::useSecOrderDeriv
bool useSecOrderDeriv
Switch toggling whether the second order derivative shall be included.
Definition: MdtSegmentFitter.h:36
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:38
xAOD::Other
@ Other
MuonR4::SegmentFitResult::segmentPars
Parameters segmentPars
Final segment parameters.
Definition: SegmentFitterEventData.h:60
MuonR4::CalibratedSpacePoint::type
xAOD::UncalibMeasType type() const
Returns the space point type.
Definition: CalibratedSpacePoint.cxx:36
MuonR4::MdtSegmentFitter::ResidualAuxillaries::projIntoWirePlane
double projIntoWirePlane
projection of the segment direction along the wire
Definition: MdtSegmentFitter.h:77
MuonR4::SegmentFit::ParamDefs::x0
@ x0
MuonR4::MdtSegmentFitter::calculateWireResiduals
void calculateWireResiduals(const LineWithPartials &line, const CalibratedSpacePoint &spacePoint, ResidualWithPartials &residual) const
Calculates the residuals together with the corresponding derivatives for a drift-circle measurement.
Definition: MdtSegmentFitter.cxx:167
WriteCalibToCool.swap
swap
Definition: WriteCalibToCool.py:94
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
compute_lumi.denom
denom
Definition: compute_lumi.py:76
MuonR4::SegmentFit::ParamDefs::time
@ time
MuonR4::CalibratedSpacePoint::State::Valid
@ Valid
MuonR4::MdtSegmentFitter::updateDerivatives
void updateDerivatives(const ResidualWithPartials &fitMeas, const MeasCov_t &measCovariance, AmgVector(5)&gradient, AmgSymMatrix(5)&hessian, double &chi2, int startPar) const
Definition: MdtSegmentFitter.cxx:584
lumiFormat.array
array
Definition: lumiFormat.py:91
MuonR4::MdtSegmentFitter::m_cfg
Config m_cfg
Definition: MdtSegmentFitter.h:111
MuonR4::SegmentFit::ParamDefs::y0
@ y0
MuonR4::MdtSegmentFitter::MdtSegmentFitter
MdtSegmentFitter(const std::string &name, Config &&config)
Standard constructor.
Definition: MdtSegmentFitter.cxx:51
MuonR4::MdtSegmentFitter::Config::nMaxCalls
unsigned int nMaxCalls
How many calls shall be executed.
Definition: MdtSegmentFitter.h:28
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:28
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:76
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::Config::reCalibrate
bool reCalibrate
Switch toggling whether the calibrator shall be called at each iteration.
Definition: MdtSegmentFitter.h:34
MuonR4::MdtSegmentFitter::UpdateStatus::allOkay
@ allOkay
MuonR4::SegmentFitResult::HitType
std::unique_ptr< CalibratedSpacePoint > HitType
Definition: SegmentFitterEventData.h:52
python.PhysicalConstants.c_light
float c_light
Definition: PhysicalConstants.py:73
MuonR4::HitVec
SpacePointPerLayerSplitter::HitVec HitVec
Definition: SpacePointPerLayerSplitter.cxx:12
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
MuonR4::MdtSegmentFitter::LineWithPartials
Store the partial derivative of the line w.r.t.
Definition: MdtSegmentFitter.h:54
Amg::dirFromAngles
Amg::Vector3D dirFromAngles(const double phi, const double theta)
Constructs a direction vector from the azimuthal & polar angles.
Definition: GeoPrimitivesHelpers.h:299
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:44
MuonR4::SegmentFit::Covariance
AmgSymMatrix(toInt(ParamDefs::nPars)) Covariance
Definition: MuonHoughDefs.h:49
a
TList * a
Definition: liststreamerinfos.cxx:10
MdtSegmentFitter.h
MuonR4::CalibratedSpacePoint::dimension
unsigned dimension() const
Returns the local dimension of the measurement.
Definition: CalibratedSpacePoint.cxx:61
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:32
MuonR4::MdtSegmentFitter::HitVec
std::vector< HitType > HitVec
Definition: MdtSegmentFitter.h:24
MuonR4::SegmentFitResult::timeFit
bool timeFit
Was the time fitted.
Definition: SegmentFitterEventData.h:56
MuonR4::CalibratedSpacePoint
The calibrated Space point is created during the calibration process.
Definition: CalibratedSpacePoint.h:15
CxxUtils::sincos
Helper to simultaneously calculate sin and cos of the same angle.
Definition: sincos.h:39
Muon::IMuonIdHelperSvc
Interface for Helper service that creates muon Identifiers and can be used to print Identifiers.
Definition: IMuonIdHelperSvc.h:27
Amg::signedDistance
double signedDistance(const Amg::Vector3D &posA, const Amg::Vector3D &dirA, const Amg::Vector3D &posB, const Amg::Vector3D &dirB)
Calculates the signed distance between two lines in 3D space.
Definition: GeoPrimitivesHelpers.h:329
MuonR4::MdtSegmentFitter::ResidualAuxillaries::projDir
Amg::Vector3D projDir
projection of the segment direction onto the wire planes
Definition: MdtSegmentFitter.h:71
MuonR4::SegmentFitResult
Definition: SegmentFitterEventData.h:46
python.Constants.VERBOSE
int VERBOSE
Definition: Control/AthenaCommon/python/Constants.py:13
PowhegControl_ttFCNC_NLO.params
params
Definition: PowhegControl_ttFCNC_NLO.py:226
MuonR4::SegmentFitResult::nPhiMeas
unsigned nPhiMeas
How many phi measurements.
Definition: SegmentFitterEventData.h:72
MuonR4::CalibratedSpacePoint::spacePoint
const SpacePoint * spacePoint() const
The pointer to the space point out of which this space point has been built.
Definition: CalibratedSpacePoint.cxx:18
MuonR4::MdtSegmentFitter::Config::defaultRanges
static RangeArray defaultRanges()
Function that returns a set of predefined ranges for testing.
Definition: MdtSegmentFitter.cxx:40
MuonR4::AmgSymMatrix
const AmgSymMatrix(2) &SpacePoint
Definition: MuonSpectrometer/MuonPhaseII/Event/MuonSpacePoint/src/SpacePoint.cxx:90
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
MuonR4::MdtSegmentFitter::recalibrate
bool recalibrate(const EventContext &ctx, SegmentFitResult &fitResult) const
Recalibrate the measurements participating in the fit & shift them into the centre-of gravity frame.
Definition: MdtSegmentFitter.cxx:62
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:622
MuonR4::MdtSegmentFitter::ResidualWithPartials::hessian
std::array< Amg::Vector3D, sumUp(nPars)> hessian
Second order derivatives.
Definition: MdtSegmentFitter.h:95
MuonR4::CalibratedSpacePoint::directionInChamber
const Amg::Vector3D & directionInChamber() const
The direction of the calibrated space point inside the chamber.
Definition: CalibratedSpacePoint.cxx:24
MuonR4::SegmentFit::ParamDefs::nPars
@ nPars
xAOD::UncalibMeasType::RpcStripType
@ RpcStripType
MuonR4::HitType
SegmentFitResult::HitType HitType
Definition: MdtSegmentFitter.cxx:36
MuonR4::CalibratedSpacePoint::driftRadius
double driftRadius() const
The drift radius of the calibrated space point.
Definition: CalibratedSpacePoint.cxx:27
MuonR4::SegmentFitResult::nTimeMeas
unsigned nTimeMeas
How many measurements give time constaint.
Definition: SegmentFitterEventData.h:74
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
python.SystemOfUnits.ns
float ns
Definition: SystemOfUnits.py:146
python.SystemOfUnits.m
float m
Definition: SystemOfUnits.py:106
MuonR4::MdtSegmentFitter::updateLinePartials
static void updateLinePartials(const Parameters &fitPars, LineWithPartials &line)
Updates the line parameters together with its first & second order partial derivatives.
Definition: MdtSegmentFitter.cxx:125
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 Tue Sep 2 2025 21:16:14 for ATLAS Offline Software by   doxygen 1.8.18