Loading [MathJax]/extensions/tex2jax.js
ATLAS Offline Software
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
Classes | Public Types | Public Member Functions | Static Public Member Functions | Private Types | Private Member Functions | Private Attributes | List of all members
MuonR4::MdtSegmentFitter Class Reference

#include <MdtSegmentFitter.h>

Inheritance diagram for MuonR4::MdtSegmentFitter:
Collaboration diagram for MuonR4::MdtSegmentFitter:

Classes

struct  Config
 
struct  LineWithPartials
 Store the partial derivative of the line w.r.t. More...
 
struct  ResidualAuxillaries
 Helper struct carrying the space for all auxillary variables needed to calculate the residual from wire measurements. More...
 
struct  ResidualWithPartials
 Helper struct carrying the residual with its derivatives. More...
 

Public Types

using ParamDefs = SegmentFit::ParamDefs
 
using Parameters = SegmentFit::Parameters
 
using HitType = std::unique_ptr< CalibratedSpacePoint >
 
using HitVec = std::vector< HitType >
 

Public Member Functions

 MdtSegmentFitter (const std::string &name, Config &&config)
 Standard constructor. More...
 
SegmentFitResult fitSegment (const EventContext &ctx, HitVec &&calibHits, const Parameters &startPars, const Amg::Transform3D &localToGlobal) const
 
void calculateWireResiduals (const LineWithPartials &line, const CalibratedSpacePoint &spacePoint, ResidualWithPartials &residual) const
 Calculates the residuals together with the corresponding derivatives for a drift-circle measurement. More...
 
void calculateStripResiduals (const LineWithPartials &line, const CalibratedSpacePoint &spacePoint, ResidualWithPartials &residual) const
 Calculates the residual together with hte correspdonding derivatives for strip measurements. More...
 
bool msgLvl (const MSG::Level lvl) const
 Test the output level. More...
 
MsgStream & msg () const
 The standard message stream. More...
 
MsgStream & msg (const MSG::Level lvl) const
 The standard message stream. More...
 
void setLevel (MSG::Level lvl)
 Change the current logging level. More...
 

Static Public Member Functions

static void updateLinePartials (const Parameters &fitPars, LineWithPartials &line)
 Updates the line parameters together with its first & second order partial derivatives. More...
 

Private Types

enum  UpdateStatus { UpdateStatus::allOkay = 0, UpdateStatus::outOfBounds = 1, UpdateStatus::noChange = 2 }
 
using MeasCov_t = CalibratedSpacePoint::Covariance_t
 Updates the chi2, its Gradient & Hessian from the measurement residual. More...
 

Private Member Functions

bool recalibrate (const EventContext &ctx, SegmentFitResult &fitResult) const
 Recalibrate the measurements participating in the fit & shift them into the centre-of gravity frame. More...
 
bool updateHitSummary (SegmentFitResult &fitResult) const
 Brief updates the hit summary from the contributing hits. More...
 
void updateDriftSigns (const Amg::Vector3D &segPos, const Amg::Vector3D &segDir, SegmentFitResult &fitRes) const
 Update the signs of the measurement. More...
 
void updateDerivatives (const ResidualWithPartials &fitMeas, const MeasCov_t &measCovariance, AmgVector(5)&gradient, AmgSymMatrix(5)&hessian, double &chi2, int startPar) const
 
template<unsigned int nDim>
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. More...
 
template<unsigned int nDim>
void blockCovariance (const AmgSymMatrix(5)&hessian, SegmentFit::Covariance &covariance) const
 
void initMessaging () const
 Initialize our message level and MessageSvc. More...
 

Private Attributes

Config m_cfg {}
 
std::string m_nm
 Message source name. More...
 
boost::thread_specific_ptr< MsgStream > m_msg_tls
 MsgStream instance (a std::cout like with print-out levels) More...
 
std::atomic< IMessageSvc * > m_imsg { nullptr }
 MessageSvc pointer. More...
 
std::atomic< MSG::Level > m_lvl { MSG::NIL }
 Current logging level. More...
 
std::atomic_flag m_initialized ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT
 Messaging initialized (initMessaging) More...
 

Detailed Description

Definition at line 19 of file MdtSegmentFitter.h.

Member Typedef Documentation

◆ HitType

Definition at line 23 of file MdtSegmentFitter.h.

◆ HitVec

Definition at line 24 of file MdtSegmentFitter.h.

◆ MeasCov_t

Updates the chi2, its Gradient & Hessian from the measurement residual.

Depending on whether the time needs to be taken into account or not the template parameter nDim is either 3 or 2, respectively.

Parameters
fitMeasResidual of the prediction from the measurement together with its derivatives
measCovarianceCovariance matrix of the measurement
gradientReference to the gradient of the chi2 which is updated by the function
hessianReference to the Hessian matrix of the chi2 which is updated by the function
chi2Reference to the chi2 value itself which is updated by the function
startParNumber of the maximum fit parameter to consider

Definition at line 163 of file MdtSegmentFitter.h.

◆ ParamDefs

Definition at line 21 of file MdtSegmentFitter.h.

◆ Parameters

Definition at line 22 of file MdtSegmentFitter.h.

Member Enumeration Documentation

◆ UpdateStatus

Enumerator
allOkay 
outOfBounds 
noChange 

Definition at line 172 of file MdtSegmentFitter.h.

172  {
173  allOkay = 0,
174  outOfBounds = 1,
175  noChange = 2,
176  };

Constructor & Destructor Documentation

◆ MdtSegmentFitter()

MuonR4::MdtSegmentFitter::MdtSegmentFitter ( const std::string &  name,
Config &&  config 
)

Standard constructor.

Parameters
nameName to be printed in the messaging
configFit configuration parameters

Definition at line 51 of file MdtSegmentFitter.cxx.

51  :
53  m_cfg{std::move(config)}{}

Member Function Documentation

◆ blockCovariance()

template<unsigned int nDim>
void MuonR4::MdtSegmentFitter::blockCovariance ( const AmgSymMatrix(5)&  hessian,
SegmentFit::Covariance covariance 
) const
private

Definition at line 601 of file MdtSegmentFitter.cxx.

602  {
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  }

◆ calculateStripResiduals()

void MuonR4::MdtSegmentFitter::calculateStripResiduals ( const LineWithPartials line,
const CalibratedSpacePoint spacePoint,
ResidualWithPartials residual 
) const

Calculates the residual together with hte correspdonding derivatives for strip measurements.

The residual is calculated by a straight line propagation of the segment onto the plane & then taking the absolute disances in the [precision] & [non-precision direction]

Parameters
lineCurrent segment line together with the first & second order derivatives
spacePointReference to the measurement to calculate the residual to
residualOutput residual object to save the actual residual and the derivatives

Update the residual accordingly

Definition at line 292 of file MdtSegmentFitter.cxx.

294  {
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  }

◆ calculateWireResiduals()

void MuonR4::MdtSegmentFitter::calculateWireResiduals ( const LineWithPartials line,
const CalibratedSpacePoint spacePoint,
ResidualWithPartials residual 
) const

Calculates the residuals together with the corresponding derivatives for a drift-circle measurement.

The residual is calculated by projecting the hit into the plane perpendicular to the wire. If the hit is a twin-hit the distance along the wire is taken into account as well.

Parameters
lineCurrent segment line together with the first & second order derivatives
spacePointReference to the measurement to calculate the residual to
residualOutput residual object to save the actual residual and the derivatives

Fetch the hit position & direction

Cache the scalar product & the lengths as they're appearing in the formulas below more often

Project the segment onto the plane

Calculate the distance from the two reference points

Distance from the segment line to the tube wire

If the tube is a twin-tube, the hit position is no longer arbitrary along the wire. Calculate the distance along the wire towards the point of closest approach.

Calculate the first derivative of the residual

No variation of the track parameters w.r.t. the time

Loop to include the second order derivatvies

Pure angular derivatives of the residual

Angular & Spatial mixed terms

Definition at line 162 of file MdtSegmentFitter.cxx.

164  {
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),
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  }

◆ fitSegment()

SegmentFitResult MuonR4::MdtSegmentFitter::fitSegment ( const EventContext &  ctx,
HitVec &&  calibHits,
const Parameters startPars,
const Amg::Transform3D localToGlobal 
) const

Cache of the partial derivatives of the line parameters w.r.t the fit parameters

Partials of the residual w.r.t. the fit parameters

Update the partial derivatives of the direction vector

First step calibrate the hits

Reset chi2

Loop over the hits to calculate the partial derivatives

If the time is fitted the drift radii need to be updated. That's properly only possible with recalibration

Time residual calculation

Loop over hits is done. Symmetrise the Hessian

Check whether the gradient is already sufficiently small

Pure eta segment fit

In the case that the time is fit & that there are no phi measurements -> compress matrix by swaping the time column with whaever the second column is... it's zero

Subtract 1 degree of freedom to take the time into account

Calculate the chi2 per measurement

Sort the measurements by ascending z

Update the covariance

Definition at line 366 of file MdtSegmentFitter.cxx.

369  {
370 
371  const Muon::IMuonIdHelperSvc* idHelperSvc{calibHits[0]->spacePoint()->msSector()->idHelperSvc()};
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()) {
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  }
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;
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);
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  }

◆ initMessaging()

void AthMessaging::initMessaging ( ) const
privateinherited

Initialize our message level and MessageSvc.

This method should only be called once.

Definition at line 39 of file AthMessaging.cxx.

40 {
42  m_lvl = m_imsg ?
43  static_cast<MSG::Level>( m_imsg.load()->outputLevel(m_nm) ) :
44  MSG::INFO;
45 }

◆ msg() [1/2]

MsgStream & AthMessaging::msg ( ) const
inlineinherited

The standard message stream.

Returns a reference to the default message stream May not be invoked before sysInitialize() has been invoked.

Definition at line 164 of file AthMessaging.h.

165 {
166  MsgStream* ms = m_msg_tls.get();
167  if (!ms) {
168  if (!m_initialized.test_and_set()) initMessaging();
169  ms = new MsgStream(m_imsg,m_nm);
170  m_msg_tls.reset( ms );
171  }
172 
173  ms->setLevel (m_lvl);
174  return *ms;
175 }

◆ msg() [2/2]

MsgStream & AthMessaging::msg ( const MSG::Level  lvl) const
inlineinherited

The standard message stream.

Returns a reference to the default message stream May not be invoked before sysInitialize() has been invoked.

Definition at line 179 of file AthMessaging.h.

180 { return msg() << lvl; }

◆ msgLvl()

bool AthMessaging::msgLvl ( const MSG::Level  lvl) const
inlineinherited

Test the output level.

Parameters
lvlThe message level to test against
Returns
boolean Indicating if messages at given level will be printed
Return values
trueMessages at level "lvl" will be printed

Definition at line 151 of file AthMessaging.h.

152 {
153  if (!m_initialized.test_and_set()) initMessaging();
154  if (m_lvl <= lvl) {
155  msg() << lvl;
156  return true;
157  } else {
158  return false;
159  }
160 }

◆ recalibrate()

bool MuonR4::MdtSegmentFitter::recalibrate ( const EventContext &  ctx,
SegmentFitResult fitResult 
) const
inlineprivate

Recalibrate the measurements participating in the fit & shift them into the centre-of gravity frame.

The parsed centerOfGravity position is updated to be the average measurement position weighted by the inverse covariance. Returns true if there're enough degrees of freedom left to fit

Parameters
ctxEventContext to fetch the calibration constants
fitResultCurrent fit state
centerOfGravityVector to save the center of gravity position later on.
invCovVector to cache the inverse covariances of the problem.

Switch off the time fit if too little degrees of freedom are left

Recalibrate the measurements

Definition at line 62 of file MdtSegmentFitter.cxx.

63  {
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  }

◆ setLevel()

void AthMessaging::setLevel ( MSG::Level  lvl)
inherited

Change the current logging level.

Use this rather than msg().setLevel() for proper operation with MT.

Definition at line 28 of file AthMessaging.cxx.

29 {
30  m_lvl = lvl;
31 }

◆ updateDerivatives()

void MuonR4::MdtSegmentFitter::updateDerivatives ( const ResidualWithPartials fitMeas,
const MeasCov_t measCovariance,
AmgVector(5)&  gradient,
AmgSymMatrix(5)&  hessian,
double &  chi2,
int  startPar 
) const
private

Definition at line 579 of file MdtSegmentFitter.cxx.

582  {
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  }

◆ updateDriftSigns()

void MuonR4::MdtSegmentFitter::updateDriftSigns ( const Amg::Vector3D segPos,
const Amg::Vector3D segDir,
SegmentFitResult fitRes 
) const
inlineprivate

Update the signs of the measurement.

Definition at line 55 of file MdtSegmentFitter.cxx.

56  {
57  for (std::unique_ptr<CalibratedSpacePoint>& meas : fitResult.calibMeasurements) {
58  meas->setDriftRadius(SegmentFitHelpers::driftSign(segPos,segDir, *meas, msg())*meas->driftRadius());
59  }
60  }

◆ updateHitSummary()

bool MuonR4::MdtSegmentFitter::updateHitSummary ( SegmentFitResult fitResult) const
inlineprivate

Brief updates the hit summary from the contributing hits.

Returns fals if there're too little valid hits.

Parameters
fitResultReference to the container containing all the measurements

Count the phi & time measurements measurements

Mdts are already counted in the measures eta category. Don't count them twice

Definition at line 89 of file MdtSegmentFitter.cxx.

89  {
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  }

◆ updateLinePartials()

void MuonR4::MdtSegmentFitter::updateLinePartials ( const Parameters fitPars,
LineWithPartials line 
)
static

Updates the line parameters together with its first & second order partial derivatives.

Parameters
fitParsSet of segment parameters in the current iteration
lineReference to the line object to be updated

x_{0} cos (phi) sin (theta) segPos = y_{0} , segDir = sin (phi) sin (theta) 0 cos theta

d segDir cos (phi) cos(theta) dSegDir -sin(phi) sin (theta) -------—= sin (phi) cos(theta) -----— = cos(phi) sin (theta) dTheta - sin (theta) dPhi 0

Definition at line 120 of file MdtSegmentFitter.cxx.

121  {
122 
127  line.gradient[toInt(ParamDefs::x0)] = Amg::Vector3D::UnitX();
128  line.gradient[toInt(ParamDefs::y0)] = Amg::Vector3D::UnitY();
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  }

◆ updateParameters()

template<unsigned int nDim>
MdtSegmentFitter::UpdateStatus MuonR4::MdtSegmentFitter::updateParameters ( Parameters currentPars,
Parameters previousPars,
Parameters currGrad,
Parameters prevGrad,
const AmgSymMatrix(5)&  hessian 
) const
private

Update step of the segment parameters using the Hessian and the gradient.

If the Hessian is definite, the currentParameters are updated according to x_{n+1} = x_{n} - H_{n}^{1} * grad(x_{n}) Otherwise, the method of steepest descent is attempted.

Parameters
currentParsBest segment estimator parameters
previousParsSegment estimator parameters from the last iteration
currGradGradient of the chi2 from this iteration
prevGradGradient of the chi2 from the previous iteration
hessianHessian estimator

Check that all parameters remain within the parameter boundary window

Recall that for 3x3 the dimensions are [y0, theta, time]

Definition at line 617 of file MdtSegmentFitter.cxx.

619  {
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) {
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){
668  }
669  if (changeMag <= m_cfg.tolerance) {
670  return UpdateStatus::noChange;
671  }
672  return UpdateStatus::allOkay;
673  }

Member Data Documentation

◆ ATLAS_THREAD_SAFE

std::atomic_flag m_initialized AthMessaging::ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT
mutableprivateinherited

Messaging initialized (initMessaging)

Definition at line 141 of file AthMessaging.h.

◆ m_cfg

Config MuonR4::MdtSegmentFitter::m_cfg {}
private

Definition at line 111 of file MdtSegmentFitter.h.

◆ m_imsg

std::atomic<IMessageSvc*> AthMessaging::m_imsg { nullptr }
mutableprivateinherited

MessageSvc pointer.

Definition at line 135 of file AthMessaging.h.

◆ m_lvl

std::atomic<MSG::Level> AthMessaging::m_lvl { MSG::NIL }
mutableprivateinherited

Current logging level.

Definition at line 138 of file AthMessaging.h.

◆ m_msg_tls

boost::thread_specific_ptr<MsgStream> AthMessaging::m_msg_tls
mutableprivateinherited

MsgStream instance (a std::cout like with print-out levels)

Definition at line 132 of file AthMessaging.h.

◆ m_nm

std::string AthMessaging::m_nm
privateinherited

Message source name.

Definition at line 129 of file AthMessaging.h.


The documentation for this class was generated from the following files:
AthMessaging::m_lvl
std::atomic< MSG::Level > m_lvl
Current logging level.
Definition: AthMessaging.h:138
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::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
MuonR4::MdtSegmentFitter::LineWithPartials::nPars
static constexpr unsigned nPars
Free parameters of the line (x0,y0,theta,phi)
Definition: MdtSegmentFitter.h:56
MuonR4::MdtSegmentFitter::UpdateStatus
UpdateStatus
Definition: MdtSegmentFitter.h:172
collListGuids.line
string line
Definition: collListGuids.py:77
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
mergePhysValFiles.start
start
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:14
MuonR4::MdtSegmentFitter::UpdateStatus::outOfBounds
@ outOfBounds
Amg::y
@ y
Definition: GeoPrimitives.h:35
xAOD::UncalibMeasType::MMClusterType
@ MMClusterType
deg
#define deg
Definition: SbPolyhedron.cxx:17
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
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
AthMessaging::m_imsg
std::atomic< IMessageSvc * > m_imsg
MessageSvc pointer.
Definition: AthMessaging.h:135
python.SystemOfUnits.ms
int ms
Definition: SystemOfUnits.py:132
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
Athena::getMessageSvc
IMessageSvc * getMessageSvc(bool quiet=false)
Definition: getMessageSvc.cxx:20
config
Definition: PhysicsAnalysis/AnalysisCommon/AssociationUtils/python/config.py:1
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
xAOD::phi
setEt phi
Definition: TrigEMCluster_v1.cxx:29
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
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
Amg::toString
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Definition: GeoPrimitivesToStringConverter.h:40
TrigConf::MSGTC::Level
Level
Definition: Trigger/TrigConfiguration/TrigConfBase/TrigConfBase/MsgStream.h:21
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
MuonR4::MdtSegmentFitter::MeasCov_t
CalibratedSpacePoint::Covariance_t MeasCov_t
Updates the chi2, its Gradient & Hessian from the measurement residual.
Definition: MdtSegmentFitter.h:163
AmgVector
AmgVector(4) T2BSTrackFilterTool
Definition: T2BSTrackFilterTool.cxx:114
Amg::x
@ x
Definition: GeoPrimitives.h:34
MuonR4::State
CalibratedSpacePoint::State State
Definition: SpacePointCalibrator.cxx:24
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:525
MuonR4::CalibratedSpacePoint::State
State
State flag to distinguish different space point states.
Definition: CalibratedSpacePoint.h:24
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::SegmentFit::ParamDefs::x0
@ x0
jobOption.theta
theta
Definition: jobOption.ParticleGun_fwd_sequence.py:13
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
HitType
HitType
Definition: FPGATrackSimHit.h:38
MuonR4::SegmentFit::AxisDefs::eta
@ eta
dot.dot
def dot(G, fn, nodesToHighlight=[])
Definition: dot.py:5
AthMessaging::msg
MsgStream & msg() const
The standard message stream.
Definition: AthMessaging.h:164
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
MuonR4::MdtSegmentFitter::m_cfg
Config m_cfg
Definition: MdtSegmentFitter.h:111
MuonR4::SegmentFit::ParamDefs::y0
@ y0
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::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
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Amg::dirFromAngles
Amg::Vector3D dirFromAngles(const double phi, const double theta)
Constructs a direction vector from the azimuthal & polar angles.
Definition: GeoPrimitivesHelpers.h:299
a
TList * a
Definition: liststreamerinfos.cxx:10
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
AthMessaging::m_nm
std::string m_nm
Message source name.
Definition: AthMessaging.h:129
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
python.Constants.VERBOSE
int VERBOSE
Definition: Control/AthenaCommon/python/Constants.py:13
AthMessaging::initMessaging
void initMessaging() const
Initialize our message level and MessageSvc.
Definition: AthMessaging.cxx:39
PowhegControl_ttFCNC_NLO.params
params
Definition: PowhegControl_ttFCNC_NLO.py:226
AthMessaging::m_msg_tls
boost::thread_specific_ptr< MsgStream > m_msg_tls
MsgStream instance (a std::cout like with print-out levels)
Definition: AthMessaging.h:132
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::SegmentFit::ParamDefs::nPars
@ nPars
xAOD::UncalibMeasType::RpcStripType
@ RpcStripType
xAOD::UncalibMeasType::MdtDriftCircleType
@ MdtDriftCircleType
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:26
MuonR4::MdtSegmentFitter::ParamDefs
SegmentFit::ParamDefs ParamDefs
Definition: MdtSegmentFitter.h:21
MuonR4::SegmentFit::ParamDefs::theta
@ theta
fitman.k
k
Definition: fitman.py:528
generate::Zero
void Zero(TH1D *hin)
Definition: generate.cxx:32
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