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  fitResult.segmentPars[toInt(ParamDefs::x0)] = 0;
114  }
115 
116  fitResult.nDoF = fitResult.nDoF - 2 - (fitResult.nPhiMeas > 0 ? 2 : 0);
117 
118  return true;
119  }
120  void MdtSegmentFitter::updateLinePartials(const Parameters& params,
121  LineWithPartials& line) {
122 
123  line.pos[Amg::x] = params[toInt(ParamDefs::x0)];
124  line.pos[Amg::y] = params[toInt(ParamDefs::y0)];
125  line.dir = Amg::dirFromAngles(params[toInt(ParamDefs::phi)],
126  params[toInt(ParamDefs::theta)]);
127  line.gradient[toInt(ParamDefs::x0)] = Amg::Vector3D::UnitX();
128  line.gradient[toInt(ParamDefs::y0)] = Amg::Vector3D::UnitY();
139  const CxxUtils::sincos theta{params[toInt(ParamDefs::theta)]},
140  phi{params[toInt(ParamDefs::phi)]};
141  line.gradient[toInt(ParamDefs::theta)] = Amg::Vector3D{phi.cs*theta.cs, phi.sn*theta.cs, -theta.sn};
142  line.gradient[toInt(ParamDefs::phi)] = Amg::Vector3D{-theta.sn *phi.sn, theta.sn*phi.cs, 0};
143  /*********************************************************************************
144  * Non-vanishing second order derivatives
145  *
146  * d^{2} segDir d^{2} segDir cos phi
147  * ------------- = - segDir , ------------ = - sin(theta) sin phi
148  * d^{2} theta d^{2} phi 0
149  *
150  * d^{2} segDir -sin phi
151  * ------------- = cos(theta) cos phi
152  * d theta dPhi 0
153  ************************************************************************************/
154  constexpr unsigned nPars = LineWithPartials::nPars;
155  constexpr int idxThetaSq = vecIdxFromSymMat<nPars>(toInt(ParamDefs::theta), toInt(ParamDefs::theta));
156  constexpr int idxPhiSq = vecIdxFromSymMat<nPars>(toInt(ParamDefs::phi), toInt(ParamDefs::phi));
157  constexpr int idxPhiTheta = vecIdxFromSymMat<nPars>(toInt(ParamDefs::theta), toInt(ParamDefs::phi));
158  line.hessian[idxThetaSq] = - Amg::Vector3D{phi.cs*theta.sn,phi.sn*theta.sn, theta.cs};
159  line.hessian[idxPhiSq] = -theta.sn * Amg::Vector3D{phi.cs, phi.sn, 0.};
160  line.hessian[idxPhiTheta] = theta.cs * Amg::Vector3D{-phi.sn, phi.cs, 0.};
161  }
162  void MdtSegmentFitter::calculateWireResiduals(const LineWithPartials& line,
163  const CalibratedSpacePoint& hit,
164  ResidualWithPartials& resObj) const {
166  const Amg::Vector3D& hitPos{hit.positionInChamber()};
167  const Amg::Vector3D& hitDir{hit.directionInChamber()};
169  resObj.projIntoWirePlane = line.dir.dot(hitDir);
170  resObj.projDirLenSq = 1. - resObj.projIntoWirePlane*resObj.projIntoWirePlane;
171  if (resObj.projDirLenSq < std::numeric_limits<float>::epsilon()) {
172  resObj.residual.setZero();
173  for (Amg::Vector3D& grad : resObj.gradient){
174  grad.setZero();
175  }
176  if (m_cfg.useSecOrderDeriv) {
177  for (Amg::Vector3D& hess : resObj.hessian) {
178  hess.setZero();
179  }
180  }
181  ATH_MSG_WARNING("Segment is parallel along the wire");
182  return;
183  }
184  resObj.invProjLenSq = 1. / resObj.projDirLenSq;
185  resObj.invProjLen = std::sqrt(resObj.invProjLenSq);
187  resObj.projDir = (line.dir - resObj.projIntoWirePlane*hitDir) * resObj.invProjLen;
189  const Amg::Vector3D hitMinSeg = hitPos - line.pos ;
191  const double lineDist = resObj.projDir.cross(hitDir).dot(hitMinSeg);
192  const double resVal = (lineDist - hit.driftRadius());
193  resObj.residual = resVal * Amg::Vector3D::Unit(toInt(AxisDefs::eta));
194  ATH_MSG_VERBOSE("Mdt drift radius: "<<hit.driftRadius()<<" distance: "<<lineDist<<";"
195  <<Amg::signedDistance(hitPos, hitDir, line.pos, line.dir)<<", residual: "<<resVal);
198  if (hit.dimension() == 2) {
199  resObj.residual[toInt(AxisDefs::phi)] = (hitMinSeg.dot(line.dir)*resObj.projIntoWirePlane - hitMinSeg.dot(hitDir)) * resObj.invProjLenSq;
200  }
201  constexpr unsigned nLinePars = LineWithPartials::nPars;
203  for (const int param : {toInt(ParamDefs::y0), toInt(ParamDefs::x0),
204  toInt(ParamDefs::theta), toInt(ParamDefs::phi)}) {
205  if (!resObj.evalPhiPars && (param == toInt(ParamDefs::x0) || param == toInt(ParamDefs::phi))){
206  continue;
207  }
208  switch (param) {
209  case toInt(ParamDefs::theta):
210  case toInt(ParamDefs::phi): {
211  resObj.partWirePlaneProj[param] = line.gradient[param].dot(hitDir);
212  resObj.partProjDir[param] = (line.gradient[param] - resObj.partWirePlaneProj[param]*hitDir) * resObj.invProjLen
213  + resObj.partWirePlaneProj[param]*resObj.projIntoWirePlane* resObj.projDir * resObj.invProjLenSq;
214  const double partialDist = resObj.partProjDir[param].cross(hitDir).dot(hitMinSeg);
215 
216  resObj.gradient[param] = partialDist * Amg::Vector3D::Unit(toInt(AxisDefs::eta));
217  if (hit.dimension() == 2) {
218  resObj.gradient[param][toInt(AxisDefs::phi)] =
219  ( hitMinSeg.dot(line.gradient[param]) * resObj.projIntoWirePlane +
220  hitMinSeg.dot(line.dir)*resObj.partWirePlaneProj[param]) * resObj.invProjLenSq
221  +2.* resObj.residual[toInt(AxisDefs::phi)]*( resObj.projIntoWirePlane * resObj.partWirePlaneProj[param]) * resObj.invProjLenSq;
222  }
223  break;
224  } case toInt(ParamDefs::y0):
225  case toInt(ParamDefs::x0): {
226  const double partialDist = - resObj.projDir.cross(hitDir).dot(line.gradient[param]);
227  resObj.gradient[param] = partialDist * Amg::Vector3D::Unit(toInt(AxisDefs::eta));
228  if (hit.dimension() == 2) {
229  resObj.gradient[param][toInt(AxisDefs::phi)] = -(line.gradient[param].dot(line.dir) * resObj.projIntoWirePlane -
230  line.gradient[param].dot(hitDir)) * resObj.invProjLenSq;
231  }
232  break;
233  }
235  default:
236  break;
237  }
238  }
239  if (!m_cfg.useSecOrderDeriv) {
240  return;
241  }
243  for (int param = toInt(ParamDefs::phi); param >=0; --param){
244  if (!resObj.evalPhiPars && (param == toInt(ParamDefs::x0) || param == toInt(ParamDefs::phi))){
245  continue;
246  }
247  for (int param1 = param; param1>=0; --param1) {
248  if (!resObj.evalPhiPars && (param1 == toInt(ParamDefs::x0) || param1 == toInt(ParamDefs::phi))){
249  continue;
250  }
251  const int lineIdx = vecIdxFromSymMat<nLinePars>(param, param1);
252  const int resIdx = vecIdxFromSymMat<toInt(ParamDefs::nPars)>(param, param1);
254  if ( (param == toInt(ParamDefs::theta) || param == toInt(ParamDefs::phi)) &&
255  (param1 == toInt(ParamDefs::theta) || param1 == toInt(ParamDefs::phi))) {
256 
257  const double partSqLineProject = line.hessian[lineIdx].dot(hitDir);
258  const Amg::Vector3D projDirPartSq = (line.hessian[lineIdx] - partSqLineProject * hitDir) * resObj.invProjLen
259  + (resObj.partWirePlaneProj[param1] * resObj.projIntoWirePlane) * resObj.invProjLenSq * resObj.partProjDir[param]
260  + (resObj.partWirePlaneProj[param] * resObj.projIntoWirePlane) * resObj.invProjLenSq * resObj.partProjDir[param1]
261  + (partSqLineProject*resObj.projIntoWirePlane) * resObj.invProjLenSq * resObj.projDir
262  + (resObj.partWirePlaneProj[param1] * resObj.partWirePlaneProj[param]) * std::pow(resObj.invProjLenSq, 2) * resObj.projDir;
263 
264  const double partialSqDist = projDirPartSq.cross(hitDir).dot(hitMinSeg);
265  resObj.hessian[resIdx] = partialSqDist * Amg::Vector3D::Unit(toInt(AxisDefs::eta));
266  if (hit.dimension() == 2) {
267  const double partialSqAlongWire =2.*resObj.residual[toInt(AxisDefs::phi)]*resObj.projIntoWirePlane*partSqLineProject * resObj.invProjLenSq
268  +2.*resObj.residual[toInt(AxisDefs::phi)]*resObj.partWirePlaneProj[param]*resObj.partWirePlaneProj[param1]*resObj.invProjLenSq
269  +2.*resObj.gradient[param1][toInt(AxisDefs::phi)]*resObj.projIntoWirePlane*resObj.partWirePlaneProj[param]*resObj.invProjLenSq
270  +2.*resObj.gradient[param][toInt(AxisDefs::phi)]*resObj.projIntoWirePlane*resObj.partWirePlaneProj[param1]*resObj.invProjLenSq
271  + hitMinSeg.dot(line.hessian[lineIdx]) *resObj.projIntoWirePlane * resObj.invProjLenSq
272  + hitMinSeg.dot(line.dir)*partSqLineProject * resObj.invProjLenSq
273  + hitMinSeg.dot(line.gradient[param])*resObj.partWirePlaneProj[param1]*resObj.invProjLenSq
274  + hitMinSeg.dot(line.gradient[param1])*resObj.partWirePlaneProj[param]*resObj.invProjLenSq;
275  resObj.hessian[resIdx][toInt(AxisDefs::phi)] = partialSqAlongWire;
276  }
277  }
279  else if (param == toInt(ParamDefs::theta) || param == toInt(ParamDefs::phi)){
280  const double partialSqDist = - resObj.partProjDir[param].cross(hitDir).dot(line.gradient[param1]);
281  resObj.hessian[resIdx] = partialSqDist * Amg::Vector3D::Unit(toInt(AxisDefs::eta));
282  if (hit.dimension() == 2) {
283  const double partialSqAlongWire = -(line.gradient[param1].dot(line.gradient[param])*resObj.projIntoWirePlane +
284  line.gradient[param1].dot(line.dir)*resObj.partWirePlaneProj[param]) * resObj.invProjLenSq
285  + 2.* resObj.gradient[param1][toInt(AxisDefs::phi)]*( resObj.projIntoWirePlane * resObj.partWirePlaneProj[param]) * resObj.invProjLenSq;
286  resObj.hessian[resIdx][toInt(AxisDefs::phi)] = partialSqAlongWire;
287  }
288  }
289  }
290  }
291  }
292  void MdtSegmentFitter::calculateStripResiduals(const LineWithPartials& line,
293  const CalibratedSpacePoint& hit,
294  ResidualWithPartials& residual) const {
295 
296 
297  const Amg::Vector3D& hitPos{hit.positionInChamber()};
298  const Amg::Vector3D normal = hit.type() != xAOD::UncalibMeasType::Other
299  ? hit.spacePoint()->planeNormal() : Amg::Vector3D::UnitZ();
300  const double planeOffSet = normal.dot(hitPos);
301 
302  const double normDot = normal.dot(line.dir);
303  if (std::abs(normDot) < std::numeric_limits<double>::epsilon()){
304  residual.residual.setZero();
305  for (Amg::Vector3D& deriv: residual.gradient) {
306  deriv.setZero();
307  }
308  if (m_cfg.useSecOrderDeriv) {
309  for (Amg::Vector3D& hessian : residual.hessian){
310  hessian.setZero();
311  }
312  }
313  ATH_MSG_WARNING("The hit parallel with the segment line "<<Amg::toString(line.dir));
314  return;
315  }
316  const double travelledDist = (planeOffSet - line.pos.dot(normal)) / normDot;
317 
318  const Amg::Vector3D planeIsect = line.pos + travelledDist * line.dir;
320  residual.residual.block<2,1>(0,0) = (planeIsect - hitPos).block<2,1>(0,0);
321 
322  for (unsigned fitPar = 0; fitPar < line.gradient.size(); ++fitPar) {
323  switch (fitPar) {
324  case toInt(ParamDefs::phi):
325  case toInt(ParamDefs::theta):{
326  const double partialDist = - travelledDist / normDot * normal.dot(line.gradient[fitPar]);
327  residual.gradient[fitPar].block<2,1>(0,0) = (travelledDist * line.gradient[fitPar] + partialDist * line.dir).block<2,1>(0,0);
328  break;
329  } case toInt(ParamDefs::y0):
330  case toInt(ParamDefs::x0): {
331  residual.gradient[fitPar].block<2,1>(0,0) = (line.gradient[fitPar] - line.gradient[fitPar].dot(normal) / normDot * line.dir).block<2,1>(0,0);
332  break;
333  }
334  default: {
335  break;
336  }
337  }
338  }
339  if (!m_cfg.useSecOrderDeriv) {
340  return;
341  }
342  constexpr unsigned nLinePars = LineWithPartials::nPars;
343  for (int param = toInt(ParamDefs::phi); param >=0; --param){
344  for (int param1 = param; param1>=0; --param1) {
345  const int lineIdx = vecIdxFromSymMat<nLinePars>(param, param1);
346  const int resIdx = vecIdxFromSymMat<toInt(ParamDefs::nPars)>(param, param1);
347  if ( (param == toInt(ParamDefs::theta) || param == toInt(ParamDefs::phi)) &&
348  (param1 == toInt(ParamDefs::theta) || param1 == toInt(ParamDefs::phi))) {
349  residual.hessian[resIdx].block<2,1>(0,0) =(
350  travelledDist * line.hessian[lineIdx] - travelledDist *(normal.dot(line.hessian[lineIdx])) / normDot * line.dir
351  - normal.dot(line.gradient[param1]) / normDot * residual.gradient[param]
352  - normal.dot(line.gradient[param]) / normDot * residual.gradient[param1]).block<2,1>(0,0);
353  } else if (param == toInt(ParamDefs::theta) || param == toInt(ParamDefs::phi)) {
354  const double gradientDisplace = normal.dot(line.gradient[param1]);
355  if (gradientDisplace > std::numeric_limits<float>::epsilon()){
356  residual.hessian[resIdx].block<2,1>(0,0) = gradientDisplace*( normal.dot(line.gradient[param]) / std::pow(normDot,2)*line.dir
357  - line.gradient[param] / normDot ).block<2,1>(0,0);
358  } else {
359  residual.hessian[resIdx].setZero();
360  }
361  }
362  }
363  }
364  }
365 
366  SegmentFitResult MdtSegmentFitter::fitSegment(const EventContext& ctx,
367  HitVec&& calibHits,
368  const Parameters& startPars,
369  const Amg::Transform3D& localToGlobal) const {
370 
371  const Muon::IMuonIdHelperSvc* idHelperSvc{calibHits[0]->spacePoint()->msSector()->idHelperSvc()};
372  using State = CalibratedSpacePoint::State;
373 
374  if (msgLvl(MSG::VERBOSE)) {
375  std::stringstream hitStream{};
376  const auto [startPos, startDir] = makeLine(startPars);
377  for (const HitType& hit : calibHits) {
378  hitStream<<" **** "<<(hit->type() != xAOD::UncalibMeasType::Other ? idHelperSvc->toString(hit->spacePoint()->identify()): "beamspot" )
379  <<" position: "<<Amg::toString(hit->positionInChamber());
380  if (hit->type() == xAOD::UncalibMeasType::MdtDriftCircleType) {
381  hitStream<<", driftRadius: "<<SegmentFitHelpers::driftSign(startPos,startDir, *hit, msg())*hit->driftRadius() ;
382  }
383  hitStream<<", channel dir: "<<Amg::toString(hit->directionInChamber())<<std::endl;
384  }
385  ATH_MSG_VERBOSE("Start segment fit with parameters "<<toString(startPars)
386  <<", plane location: "<<Amg::toString(localToGlobal)<<std::endl<<hitStream.str());
387 
388  }
389 
390  SegmentFitResult fitResult{};
391  fitResult.segmentPars = startPars;
392  fitResult.timeFit = m_cfg.doTimeFit;
393  fitResult.calibMeasurements = std::move(calibHits);
394  if (!updateHitSummary(fitResult)) {
395  ATH_MSG_WARNING("No valid segment seed parsed from the beginning");
396  return fitResult;
397  }
398  if (!m_cfg.reCalibrate) {
399  const auto [segPos, segDir] = makeLine(startPars);
400  updateDriftSigns(segPos, segDir, fitResult);
401  }
402  Parameters gradient{AmgVector(5)::Zero()}, prevGrad{AmgVector(5)::Zero()}, prevPars{AmgVector(5)::Zero()};
403  AmgSymMatrix(5) hessian{AmgSymMatrix(5)::Zero()};
404 
406  LineWithPartials segmentLine{};
408  ResidualWithPartials residual{};
409 
410  unsigned int noChangeIter{0};
411  while (fitResult.nIter++ < m_cfg.nMaxCalls) {
412  ATH_MSG_DEBUG("Iteration: "<<fitResult.nIter<<" parameters: "<<toString(fitResult.segmentPars)<<", chi2: "<<fitResult.chi2);
414  updateLinePartials(fitResult.segmentPars, segmentLine);
415 
417  if (m_cfg.reCalibrate && !recalibrate(ctx, fitResult)) {
418  break;
419  }
421  fitResult.chi2 = 0;
422  hessian.setZero();
423  gradient.setZero();
424 
426  for (HitType& hit : fitResult.calibMeasurements) {
427  if (hit->fitState() != State::Valid) {
428  continue;
429  }
430  const Amg::Vector3D& hitPos{hit->positionInChamber()};
431  const Amg::Vector3D& hitDir{hit->directionInChamber()};
432  ATH_MSG_DEBUG("Update chi2 from measurement "<<(hit->spacePoint() ? idHelperSvc->toString(hit->spacePoint()->identify())
433  : "pseudo meas")<<" position: "<<Amg::toString(hitPos)<<" + "<<Amg::toString(hitDir));
434 
435  const int start = toInt(fitResult.nPhiMeas ? ParamDefs::phi : ParamDefs::theta);
436  switch (hit->type()) {
437  case xAOD::UncalibMeasType::MdtDriftCircleType: {
440  const double dSign = (hit->driftRadius() > 0 ? 1. : -1.);
441  if (fitResult.timeFit && !m_cfg.reCalibrate && fitResult.nIter > 1) {
442  hit = m_cfg.calibrator->calibrate(ctx, *hit, segmentLine.pos, segmentLine.dir,
443  fitResult.segmentPars[toInt(ParamDefs::time)]);
444  hit->setDriftRadius(dSign*hit->driftRadius());
445  }
446  calculateWireResiduals(segmentLine, *hit, residual);
447 
448  if (!fitResult.timeFit) {
449  break;
450  }
452  residual.gradient[toInt(ParamDefs::time)] = dSign*m_cfg.calibrator->driftVelocity(ctx, *hit) * Amg::Vector3D::Unit(toInt(AxisDefs::eta));
453  if (m_cfg.useSecOrderDeriv) {
454  constexpr unsigned hessIdx = vecIdxFromSymMat<ResidualWithPartials::nPars>(toInt(ParamDefs::time), toInt(ParamDefs::time));
455  residual.hessian[hessIdx] = dSign*m_cfg.calibrator->driftAcceleration(ctx, *hit) * Amg::Vector3D::Unit(toInt(AxisDefs::eta));
456  }
457  break;
458  }
459  case xAOD::UncalibMeasType::RpcStripType:
460  case xAOD::UncalibMeasType::TgcStripType:
461  case xAOD::UncalibMeasType::sTgcStripType:
462  case xAOD::UncalibMeasType::MMClusterType: {
463  calculateStripResiduals(segmentLine, *hit, residual);
464  if (!fitResult.timeFit || !hit->measuresTime()) {
465  break;
466  }
467  constexpr ParamDefs par = ParamDefs::time;
468  residual.gradient[toInt(par)] = -Amg::Vector3D::Unit(toInt(AxisDefs::t0));
469  const Amg::Vector3D planeIsect = residual.residual + hitPos;
470  const double totFlightDist = (localToGlobal * planeIsect).mag() * c_inv;
471  residual.residual[toInt(AxisDefs::t0)] = hit->time() - totFlightDist * c_inv - fitResult.segmentPars[toInt(ParamDefs::time)];
472 
473  for (int p = start; p >= 0; --p) {
474  residual.gradient[p][toInt(AxisDefs::t0)] = -residual.gradient[p].perp() * c_inv;
475  ATH_MSG_VERBOSE("Partial derivative of "<<idHelperSvc->toString(hit->spacePoint()->identify())
476  <<" residual "<<Amg::toString(residual.residual)<<" "<<Amg::toString(multiply(inverse(hit->covariance()), residual.residual))
477  <<" "<<std::endl<<toString(hit->covariance())<<std::endl<<" w.r.t "
478  <<toString(static_cast<ParamDefs>(p))<<"="<<Amg::toString(residual.gradient[p]));
479  }
480  break;
481  } case xAOD::UncalibMeasType::Other: {
482  calculateStripResiduals(segmentLine, *hit, residual);
483  break;
484  } default: {
485  const auto &measurement = *hit->spacePoint()->primaryMeasurement();
486  ATH_MSG_WARNING("MdtSegmentFitter() - Unsupported measurment type" <<typeid(measurement).name());
487  }
488  }
489 
490  ATH_MSG_VERBOSE("Update derivatives for hit "<< (hit->spacePoint() ? idHelperSvc->toString(hit->spacePoint()->identify()) : "beamspot"));
491  updateDerivatives(residual, hit->covariance(), gradient, hessian, fitResult.chi2,
492  fitResult.timeFit && hit->measuresTime() ? toInt(ParamDefs::time) : start);
493  }
494 
496  for (int k =1; k < 5 - (!fitResult.timeFit); ++k){
497  for (int l = 0; l< k; ++l){
498  hessian(l,k) = hessian(k,l);
499  }
500  }
502  if (gradient.mag() < m_cfg.tolerance) {
503  fitResult.converged = true;
504  ATH_MSG_VERBOSE("Fit converged after "<<fitResult.nIter<<" iterations with "<<fitResult.chi2);
505  break;
506  }
507  ATH_MSG_VERBOSE("Chi2: "<<fitResult.chi2<<", gradient: "<<Amg::toString(gradient)<<"hessian: "<<std::endl<<hessian);
509  UpdateStatus paramUpdate{UpdateStatus::outOfBounds};
510  if (!fitResult.nPhiMeas && !fitResult.timeFit) {
511  paramUpdate = updateParameters<2>(fitResult.segmentPars, prevPars, gradient, prevGrad, hessian);
512  } else if (!fitResult.nPhiMeas && fitResult.timeFit) {
515  hessian.col(2).swap(hessian.col(toInt(ParamDefs::time)));
516  hessian.row(2).swap(hessian.row(toInt(ParamDefs::time)));
517  std::swap(gradient[2], gradient[toInt(ParamDefs::time)]);
518  std::swap(fitResult.segmentPars[2], fitResult.segmentPars[toInt(ParamDefs::time)]);
519 
520  paramUpdate = updateParameters<3>(fitResult.segmentPars, prevPars, gradient, prevGrad, hessian);
521 
522  std::swap(fitResult.segmentPars[2], fitResult.segmentPars[toInt(ParamDefs::time)]);
523  hessian.col(2).swap(hessian.col(toInt(ParamDefs::time)));
524  hessian.row(2).swap(hessian.row(toInt(ParamDefs::time)));
525  std::swap(gradient[2], gradient[toInt(ParamDefs::time)]);
526  } else if (fitResult.nPhiMeas && !fitResult.timeFit) {
527  paramUpdate = updateParameters<4>(fitResult.segmentPars, prevPars, gradient, prevGrad, hessian);
528  } else if (fitResult.nPhiMeas && fitResult.timeFit) {
529  paramUpdate = updateParameters<5>(fitResult.segmentPars, prevPars, gradient, prevGrad, hessian);
530  }
531  if (paramUpdate == UpdateStatus::noChange) {
532  if ((++noChangeIter) >= m_cfg.noMoveIter){
533  fitResult.converged = true;
534  break;
535  }
536  } else if (paramUpdate == UpdateStatus::allOkay) {
537  noChangeIter = 0;
538  } else if (paramUpdate == UpdateStatus::outOfBounds){
539  return fitResult;
540  }
541  }
542 
544  fitResult.nDoF-=fitResult.timeFit;
545 
547  const auto [segPos, segDir] = fitResult.makeLine();
548  fitResult.chi2 =0.;
549 
550  std::optional<double> toF = fitResult.timeFit ? std::make_optional<double>((localToGlobal * segPos).mag() * c_inv) : std::nullopt;
552  std::ranges::stable_sort(fitResult.calibMeasurements, [](const HitType&a, const HitType& b){
553  return a->positionInChamber().z() < b->positionInChamber().z();
554  });
555 
556  /*** Remove the drift sign again */
557  fitResult.chi2 =0.;
558  for (const HitType& hit : fitResult.calibMeasurements) {
559  hit->setDriftRadius(std::abs(hit->driftRadius()));
560  fitResult.chi2 +=SegmentFitHelpers::chiSqTerm(segPos, segDir, fitResult.segmentPars[toInt(ParamDefs::time)], toF, *hit, msg());
561  }
563  if (!fitResult.nPhiMeas&& !fitResult.timeFit) {
564  blockCovariance<2>(std::move(hessian), /* std::move(gradient), std::move(prevGrad), */ fitResult.segmentParErrs);
565  }else if (!fitResult.nPhiMeas && fitResult.timeFit) {
566  hessian.col(2).swap(hessian.col(toInt(ParamDefs::time)));
567  hessian.row(2).swap(hessian.row(toInt(ParamDefs::time)));
568  blockCovariance<3>(std::move(hessian),/* std::move(gradient), std::move(prevGrad), */ fitResult.segmentParErrs);
569  fitResult.segmentParErrs.col(2).swap(fitResult.segmentParErrs.col(toInt(ParamDefs::time)));
570  fitResult.segmentParErrs.row(2).swap(fitResult.segmentParErrs.row(toInt(ParamDefs::time)));
571  } else if (fitResult.nPhiMeas) {
572  blockCovariance<4>(std::move(hessian), /* std::move(gradient), std::move(prevGrad), */ fitResult.segmentParErrs);
573  } else if (fitResult.nPhiMeas && fitResult.timeFit) {
574  blockCovariance<5>(std::move(hessian),/* std::move(gradient), std::move(prevGrad), */ fitResult.segmentParErrs);
575  }
576  return fitResult;
577  }
578 
579  void MdtSegmentFitter::updateDerivatives(const ResidualWithPartials& fitMeas,
580  const MeasCov_t& covariance,
581  AmgVector(5)& gradient, AmgSymMatrix(5)& hessian,
582  double& chi2, int startPar) const {
583  const MeasCov_t invCov{inverse(covariance)};
584  const Amg::Vector3D covRes = multiply(invCov, fitMeas.residual);
585  chi2 += covRes.dot(fitMeas.residual);
586  for (int p = startPar; p >=0 ; --p) {
587  gradient[p] +=2.*covRes.dot(fitMeas.gradient[p]);
588  for (int k=p; k>=0; --k) {
589  hessian(p,k)+= 2.*contract(invCov, fitMeas.gradient[p], fitMeas.gradient[k]);
590  if (m_cfg.useSecOrderDeriv) {
591  const int symMatIdx = vecIdxFromSymMat<ResidualWithPartials::nPars>(p,k);
592  hessian(p,k)+=contract(invCov, covRes, fitMeas.hessian[symMatIdx]);
593  }
594  }
595  }
596  ATH_MSG_VERBOSE("After derivative update --- chi2: "<<chi2<<"("<<covRes.dot(fitMeas.residual)<<"), gradient: "
597  <<toString(gradient)<<", Hessian:\n"<<hessian<<", measurement covariance\n"<<toString(invCov));
598  }
599 
600  template <unsigned int nDim>
601  void MdtSegmentFitter::blockCovariance(const AmgSymMatrix(5)& hessian,
602  SegmentFit::Covariance& covariance) const {
603 
604  covariance.setIdentity();
605  AmgSymMatrix(nDim) miniHessian = hessian.block<nDim, nDim>(0,0);
606  if (std::abs(miniHessian.determinant()) <= detCutOff) {
607  ATH_MSG_VERBOSE("Boeser mini hessian ("<<miniHessian.determinant()<<")\n"<<miniHessian<<"\n\n"<<hessian);
608  return;
609  }
610  ATH_MSG_VERBOSE("Hessian matrix: \n"<<hessian<<",\nblock Hessian:\n"<<miniHessian<<",\n determinant: "<<miniHessian.determinant());
611  covariance.block<nDim,nDim>(0,0) = miniHessian.inverse();
612  ATH_MSG_VERBOSE("covariance: \n"<<covariance);
613  }
614 
615  template <unsigned int nDim>
616  MdtSegmentFitter::UpdateStatus
617  MdtSegmentFitter::updateParameters(Parameters& currPars, Parameters& prevPars,
618  Parameters& currGrad, Parameters& prevGrad,
619  const AmgSymMatrix(5)& currentHessian) const {
620 
621  AmgSymMatrix(nDim) miniHessian = currentHessian.block<nDim, nDim>(0,0);
622  ATH_MSG_DEBUG("Parameter update -- \ncurrenPars: "<<toString(currPars)<<", \ngradient: "<<toString(currGrad)
623  <<", hessian ("<<miniHessian.determinant()<<")"<<std::endl<<miniHessian);
624 
625  double changeMag{0.};
626  if (std::abs(miniHessian.determinant()) > detCutOff) {
627  prevPars.block<nDim,1>(0,0) = currPars.block<nDim,1>(0,0);
628  // Update the parameters accrodingly to the hessian
629  const AmgVector(nDim) updateMe = miniHessian.inverse()* currGrad.block<nDim, 1>(0,0);
630  changeMag = std::sqrt(updateMe.dot(updateMe));
631  currPars.block<nDim,1>(0,0) -= updateMe;
632  prevGrad.block<nDim,1>(0,0) = currGrad.block<nDim,1>(0,0);
633  ATH_MSG_DEBUG("Hessian inverse:\n"<<miniHessian.inverse()<<"\n\nUpdate the parameters by -"
634  <<Amg::toString(updateMe));
635  } else {
636  const AmgVector(nDim) gradDiff = (currGrad - prevGrad).block<nDim,1>(0,0);
637  const double gradDiffMag = gradDiff.mag2();
638  double denom = (gradDiffMag > std::numeric_limits<float>::epsilon() ? gradDiffMag : 1.);
639  const double gamma = std::abs((currPars - prevPars).block<nDim,1>(0,0).dot(gradDiff))
640  / denom;
641  ATH_MSG_VERBOSE("Hessian determinant invalid. Try deepest descent - \nprev parameters: "
642  <<toString(prevPars)<<",\nprevious gradient: "<<toString(prevGrad)<<", gamma: "<<gamma);
643  prevPars.block<nDim, 1>(0,0) = currPars.block<nDim, 1>(0,0);
644  currPars.block<nDim, 1>(0,0) -= gamma* currGrad.block<nDim, 1>(0,0);
645  prevGrad.block<nDim, 1>(0,0) = currGrad.block<nDim,1>(0,0);
646  changeMag = std::abs(gamma) * currGrad.block<nDim, 1>(0,0).mag();
647  }
649  unsigned int nOutOfBound{0};
650  for (unsigned int p = 0; p< nDim; ++p) {
651  double& parValue{currPars[p]};
652  if constexpr(nDim == 3) {
654  if (p == 2) {
655  p = toInt(ParamDefs::time);
656  }
657  }
658  if (m_cfg.ranges[p][0] > parValue || m_cfg.ranges[p][1] < parValue) {
659  ATH_MSG_VERBOSE("The "<<p<<"-th parameter "<<toString(static_cast<ParamDefs>(p))<<" is out of range "<<parValue
660  <<" allowed ["<<m_cfg.ranges[p][0]<<"-"<<m_cfg.ranges[p][1]<<"]");
661  ++nOutOfBound;
662  parValue = std::clamp(parValue, m_cfg.ranges[p][0], m_cfg.ranges[p][1]);
663  }
664 
665  }
666  if (nOutOfBound > m_cfg.nParsOutOfBounds){
667  return UpdateStatus::outOfBounds;
668  }
669  if (changeMag <= m_cfg.tolerance) {
670  return UpdateStatus::noChange;
671  }
672  return UpdateStatus::allOkay;
673  }
674 }
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:601
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:292
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
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
collListGuids.line
string line
Definition: collListGuids.py:77
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:14
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:128
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: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 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
MuonR4::MdtSegmentFitter::fitSegment
SegmentFitResult fitSegment(const EventContext &ctx, HitVec &&calibHits, const Parameters &startPars, const Amg::Transform3D &localToGlobal) const
Definition: MdtSegmentFitter.cxx:366
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:235
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:210
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: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: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:162
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:579
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: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::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:63
MuonR4::HitVec
SpacePointPerLayerSplitter::HitVec HitVec
Definition: SpacePointPerLayerSplitter.cxx:11
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
python.SystemOfUnits.ns
int ns
Definition: SystemOfUnits.py:130
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:150
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:617
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
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:120
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 Sun Apr 20 2025 21:14:46 for ATLAS Offline Software by   doxygen 1.8.18