ATLAS Offline Software
Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes | List of all members
Trk::FitMatrices Class Reference

#include <FitMatrices.h>

Collaboration diagram for Trk::FitMatrices:

Public Member Functions

 FitMatrices (bool constrainedAlignmentEffects)
 
 FitMatrices (const FitMatrices &)=delete
 
FitMatricesoperator= (const FitMatrices &)=delete
 
 FitMatrices (FitMatrices &&)=default
 
FitMatricesoperator= (FitMatrices &&)=default
 
 ~FitMatrices ()=default
 
void checkPointers (MsgStream &log) const
 
Amg::MatrixXderivativeMatrix ()
 
Amg::MatrixXfinalCovariance ()
 
const Amg::MatrixXfullCovariance (void)
 
int numberDoF (void) const
 
int numberDriftCircles (void) const
 
double perigeeChiSquared (void)
 
void printDerivativeMatrix (void)
 
void printWeightMatrix (void)
 
void refinePointers (void)
 
void releaseMemory (void)
 
int setDimensions (std::vector< FitMeasurement * > &measurements, FitParameters *parameters)
 
bool solveEquations (void)
 
void usePerigee (const FitMeasurement &measurement)
 

Static Public Member Functions

static double chiSquaredChange (void)
 

Private Member Functions

void addPerigeeMeasurement (void)
 
void avoidMomentumSingularity (void)
 

Private Attributes

fitMatrix m_fitMatrix
 
int m_columnsDM
 
bool m_constrainedAlignmentEffects
 
Amg::MatrixX m_covariance {}
 
Amg::MatrixX m_derivativeMatrix {}
 
Amg::MatrixX m_finalCovariance {}
 
std::vector< int > m_firstRowForParameter
 
double m_largePhiWeight
 
std::vector< int > m_lastRowForParameter
 
bool m_matrixFromCLHEP
 
std::vector< FitMeasurement * > * m_measurements
 
int m_numberDoF
 
int m_numberDriftCircles
 
int m_numberPerigee
 
FitParametersm_parameters
 
const Amg::VectorXm_perigee
 
Amg::MatrixX m_perigeeDifference
 
const Amg::MatrixXm_perigeeWeight
 
std::vector< double > m_residuals
 
int m_rowsDM
 
bool m_usePerigee
 
Amg::MatrixX m_weight
 
Amg::VectorX m_weightedDifference
 

Detailed Description

Definition at line 41 of file FitMatrices.h.

Constructor & Destructor Documentation

◆ FitMatrices() [1/3]

Trk::FitMatrices::FitMatrices ( bool  constrainedAlignmentEffects)

Definition at line 41 of file FitMatrices.cxx.

42  : m_fitMatrix{},
43  m_columnsDM(0),
44  m_constrainedAlignmentEffects(constrainedAlignmentEffects),
45  m_largePhiWeight(10000.), // arbitrary - equiv to 10um
46  m_matrixFromCLHEP(false),
47  m_measurements(nullptr),
48  m_numberDoF(0),
50  m_numberPerigee(5),
51  m_parameters(nullptr),
52  m_perigee(nullptr),
54  m_perigeeWeight(nullptr),
55  m_rowsDM(0),
56  m_usePerigee(false){}

◆ FitMatrices() [2/3]

Trk::FitMatrices::FitMatrices ( const FitMatrices )
delete

◆ FitMatrices() [3/3]

Trk::FitMatrices::FitMatrices ( FitMatrices &&  )
default

◆ ~FitMatrices()

Trk::FitMatrices::~FitMatrices ( )
default

Member Function Documentation

◆ addPerigeeMeasurement()

void Trk::FitMatrices::addPerigeeMeasurement ( void  )
private

Definition at line 621 of file FitMatrices.cxx.

621  {
622  // TODO: needs eigen equiv !!
623  // const Amg::MatrixX& perigeeWeight = *m_perigeeWeight;
624  // Amg::MatrixX& weight = *m_weightCLHEP;
625  // AmgVectorX& weightedDifference = *m_weightedDifferenceCLHEP;
626  // AmgVectorX diff_vector =
627  // perigeeWeight*m_perigeeDifference.T(); for (int row = 0; row <
628  // m_numberPerigee; ++row)
629  // {
630  // weightedDifference[row] += diff_vector[row];
631  // for (int col = 0; col <= row; ++col) weight[row][col] +=
632  // perigeeWeight[row][col];
633  // }
634 }

◆ avoidMomentumSingularity()

void Trk::FitMatrices::avoidMomentumSingularity ( void  )
private

Definition at line 636 of file FitMatrices.cxx.

636  {
637  // fix momentum if line-fit or fit attempted with negligible field integral
640  weight(5, 5) < 1. / Gaudi::Units::TeV) {
641  for (int i = 0; i != m_columnsDM; ++i) {
642  weight(i, 5) = 0.;
643  weight(5, i) = 0.;
644  }
645  weight(5, 5) += 1. / Gaudi::Units::TeV;
646  }
647  if (!m_parameters->fitMomentum() || weight(4, 4) < 1. / Gaudi::Units::TeV) {
648  m_parameters->fitMomentum(false);
649  for (int i = 0; i != m_columnsDM; ++i) {
650  weight(i, 4) = 0.;
651  weight(4, i) = 0.;
652  }
653  weight(4, 4) += 1. / Gaudi::Units::TeV;
654  }
655 }

◆ checkPointers()

void Trk::FitMatrices::checkPointers ( MsgStream &  log) const

Definition at line 58 of file FitMatrices.cxx.

58  {
59  // debugging: check smart pointers
60  for (int col = 0; col < m_columnsDM; ++col) {
61  // firstRow
62  for (int i = 0; i < m_firstRowForParameter[col]; ++i) {
63  if (m_fitMatrix.derivative[i][col] != 0.)
64  log << " col " << col << " unexpected first nonzero DM element " << i
65  << " / " << m_firstRowForParameter[col] << endmsg;
66  }
67  int j = m_firstRowForParameter[col];
68  if (m_fitMatrix.derivative[j][col] == 0.)
69  log << " col " << col << " first nonzero DM element is zero! " << j
70  << " / " << m_firstRowForParameter[col] << endmsg;
71 
72  // lastRow
73  for (int i = m_lastRowForParameter[col]; i < m_rowsDM; ++i) {
74  if (m_fitMatrix.derivative[i][col] != 0.)
75  log << " col " << col << " unexpected last nonzero DM element " << i
76  << " / " << m_lastRowForParameter[col] << endmsg;
77  }
78  }
79 }

◆ chiSquaredChange()

double Trk::FitMatrices::chiSquaredChange ( void  )
static

Definition at line 81 of file FitMatrices.cxx.

81  {
83  std::cout << " unexpected :chiSquaredChange " << std::endl;
84  return 0.;
86 }

◆ derivativeMatrix()

Amg::MatrixX & Trk::FitMatrices::derivativeMatrix ( )
inline

Definition at line 142 of file FitMatrices.h.

142  {
143  return m_derivativeMatrix;
144 }

◆ finalCovariance()

Amg::MatrixX & Trk::FitMatrices::finalCovariance ( )
inline

Definition at line 146 of file FitMatrices.h.

146  {
147  return m_finalCovariance;
148 }

◆ fullCovariance()

const Amg::MatrixX * Trk::FitMatrices::fullCovariance ( void  )

Definition at line 88 of file FitMatrices.cxx.

88  {
89  // return result if matrix already inverted
90  if (m_covariance.size()!=0) {
91  return &m_covariance;
92  }
94 
95  // fix weighting ???? shouldn't we just remove large phi weight?
98  }
99 
100  // invert weight matrix
101  Amg::MatrixX& covariance = m_covariance;
102  // avoid singularity through ill-defined momentum ???? again
104 
105  // neater - but gives small rounding-like diffs wrt matrix copy version
106  // keep matrix copy for release 21 to avoid rounding changes at Tier0
107  // covariance = (*m_weight).inverse();
108 
109  // matrix copy version (legacy of older library which needed copy between
110  // matrix packages)
112  weight.selfadjointView<0x2>();
113 
114  // check if m_weights makes sense before inverting
116  m_covariance.resize(0, 0);
117  return nullptr;
118  }
119 
120  weight = (m_weight).inverse();
121  for (int row = 0; row != m_columnsDM; ++row) {
122  for (int col = 0; col != m_columnsDM; ++col)
123  covariance(row, col) = weight(col, row);
124  }
125 
126  // back convert curved fits to Tracking units (MeV)
127  if (m_parameters->fitMomentum()) {
128  // transform to MeV
129  double d4 = 1. / Gaudi::Units::TeV;
130  for (int row = 0; row < m_columnsDM; ++row) {
131  covariance(4, row) *= d4;
132  covariance(row, 4) = covariance(4, row);
133  }
134  covariance(4, 4) *= d4;
135 
136  // transform units for fitted energy deposit (for now fit qOverP at calo
137  // exit)
139  double d5 = 1. / Gaudi::Units::TeV;
140  for (int row = 0; row < m_columnsDM; ++row) {
141  covariance(5, row) *= d5;
142  covariance(row, 5) = covariance(5, row);
143  }
144  covariance(5, 5) *= d5;
145  }
146  }
147 
148  // FIXME: errors underestimated on d0,z0 when large scatterer precedes precise
149  // measurements final covariance starts with 5*5 representing perigee from
150  // full covariance
152  for (int i = 0; i != 5; ++i) {
153  for (int j = 0; j != 5; ++j) {
154  (m_finalCovariance)(i, j) = covariance(i, j);
155  }
156  }
157 
158  // return pointer to full covariance
159  return &m_covariance;
160 }

◆ numberDoF()

int Trk::FitMatrices::numberDoF ( void  ) const
inline

Definition at line 150 of file FitMatrices.h.

150  {
151  return m_numberDoF;
152 }

◆ numberDriftCircles()

int Trk::FitMatrices::numberDriftCircles ( void  ) const
inline

Definition at line 154 of file FitMatrices.h.

154  {
155  return m_numberDriftCircles;
156 }

◆ operator=() [1/2]

FitMatrices& Trk::FitMatrices::operator= ( const FitMatrices )
delete

◆ operator=() [2/2]

FitMatrices& Trk::FitMatrices::operator= ( FitMatrices &&  )
default

◆ perigeeChiSquared()

double Trk::FitMatrices::perigeeChiSquared ( void  )

Definition at line 162 of file FitMatrices.cxx.

162  {
164  return (m_perigeeDifference * (*m_perigeeWeight) *
165  m_perigeeDifference.transpose())(0, 0);
166 }

◆ printDerivativeMatrix()

void Trk::FitMatrices::printDerivativeMatrix ( void  )

Definition at line 168 of file FitMatrices.cxx.

168  {
169  std::cout << "DerivativeMatrix: rows * columns " << m_rowsDM << " * "
170  << m_columnsDM << " numberDoF " << m_numberDoF << std::endl;
171 
172  int firstCol = 0;
173  int lastCol = 0;
174  if (!m_measurements)
175  return;
177  std::vector<FitMeasurement*> alignmentfm;
178  FitMeasurement* fm = *m;
179  bool singleRow = true;
180 
181  for (int row = 0; row < m_rowsDM; ++row) {
182  // get corresponding FitMeasurement
183  if (singleRow) {
184  while (m != m_measurements->end() &&
185  (!(**m).numberDoF() || (**m).isAlignment())) {
186  if ((**m).isAlignment())
187  alignmentfm.push_back(*m);
188  ++m;
189  }
190  if (m != m_measurements->end()) {
191  fm = *m;
192  } else {
193  fm = alignmentfm.back();
194  alignmentfm.pop_back();
195  }
196 
197  firstCol = fm->firstParameter();
198  lastCol = fm->lastParameter() - 1;
199  std::cout << std::endl << std::setiosflags(std::ios::fixed);
200  if (fm->isPositionMeasurement()) {
201  std::cout << "measurement";
202  } else if (fm->isScatterer()) {
203  std::cout << "scatterer ";
204  } else if (fm->isEnergyDeposit()) {
205  std::cout << "energyDepos";
206  } else if (fm->isAlignment()) {
207  std::cout << "alignment ";
208  } else {
209  std::cout << " ?? ";
210  }
211  if (fm->is2Dimensional()) {
212  singleRow = false;
213  } else {
214  if (m != m_measurements->end())
215  ++m;
216  }
217  std::cout << " row " << std::setw(3) << row << " col 0 ";
218  } else {
219  if (m != m_measurements->end())
220  ++m;
221  singleRow = true;
222  std::cout << std::endl
223  << std::setiosflags(std::ios::fixed) << " row "
224  << std::setw(3) << row << " col 0 ";
225  }
226  for (int col = 0; col < m_columnsDM; ++col) {
227  if (col < firstCol || col > lastCol) // m_firstRowForParameter[row])
228  {
229  if (m_fitMatrix.derivative[row][col] == 0.) {
230  std::cout << " ";
231  } else {
232  // flag out-of-order
233  std::cout << std::setiosflags(std::ios::scientific)
234  << std::setbase(10) << "<" << std::setw(10)
235  << m_fitMatrix.derivative[row][col] << ">";
236  }
237  } else {
238  std::cout << std::setiosflags(std::ios::scientific) << std::setbase(10)
239  << std::setw(10) << m_fitMatrix.derivative[row][col] << " ";
240  }
241 
242  if ((col + 1) % 12 == 0 && col + 1 < m_columnsDM)
243  std::cout << std::endl
244  << std::setiosflags(std::ios::fixed)
245  << " col " << std::setw(3) << col + 1
246  << " ";
247  }
248  }
249  std::cout << std::endl;
250 }

◆ printWeightMatrix()

void Trk::FitMatrices::printWeightMatrix ( void  )

Definition at line 252 of file FitMatrices.cxx.

252  {
253  std::cout << std::endl
254  << "WeightMatrix: symmetric with rank " << m_columnsDM;
255 
256  for (int row = 0; row < m_columnsDM; ++row) {
257  std::cout << std::endl
258  << std::setiosflags(std::ios::fixed) << " row " << std::setw(3)
259  << row << " col 0 ";
260  for (int col = 0; col <= row; ++col) {
261  std::cout << std::setiosflags(std::ios::scientific) << std::setbase(10)
262  << std::setw(10) << (m_weight)(row, col) << " ";
263 
264  if ((col + 1) % 13 == 0 && col < row)
265  std::cout << std::endl
266  << std::setiosflags(std::ios::fixed) << " col "
267  << std::setw(3) << col + 1 << " ";
268  }
269  }
270  std::cout << std::endl;
271 }

◆ refinePointers()

void Trk::FitMatrices::refinePointers ( void  )

Definition at line 273 of file FitMatrices.cxx.

273  {
274  // remove leading and trailing zeroes from smart pointers
275  for (int col = 0; col < m_columnsDM; ++col) {
277  int j = m_lastRowForParameter[col];
278  if (i < --j && m_fitMatrix.derivative[i][col] == 0. ) {
279  while (i != j && m_fitMatrix.derivative[i][col] == 0.)
280  ++i;
282  }
283  if (m_fitMatrix.derivative[j][col] == 0. && j > i) {
284  while (j != i && m_fitMatrix.derivative[j][col] == 0.)
285  --j;
286  m_lastRowForParameter[col] = ++j;
287  }
288  }
289 }

◆ releaseMemory()

void Trk::FitMatrices::releaseMemory ( void  )

Definition at line 291 of file FitMatrices.cxx.

291  {
292  m_derivativeMatrix.resize(0,0);
293  m_weight.resize(0,0);
294  m_weightedDifference.resize(0,0);
295 }

◆ setDimensions()

int Trk::FitMatrices::setDimensions ( std::vector< FitMeasurement * > &  measurements,
FitParameters parameters 
)

Definition at line 297 of file FitMatrices.cxx.

298  {
299  // keep pointer for debug purposes
300  m_measurements = &measurements;
301 
302  // only use perigee on request (from special fit types)
303  m_usePerigee = false;
304 
305  // count rows, misalignments and scatterers from loop over FitMeasurements
306  m_firstRowForParameter.clear();
307  m_firstRowForParameter.reserve(128);
308  m_lastRowForParameter.clear();
309  m_lastRowForParameter.reserve(128);
311  m_residuals = std::vector<double>(2 * measurements.size(), 0.);
313  bool haveMeasurement = false;
314  bool haveVertex = false;
315  int numberAlignments = 0;
316  int numberEnergyDeposits = 0;
317  int numberParameters = 5;
318  int numberScatterers = 0;
319  int row = 0;
320 
321  // keep first row with measurements up to each number of parameters
322  m_firstRowForParameter = std::vector<int>(numberParameters, -1);
323  std::vector<FitMeasurement*>::iterator m = measurements.begin();
324  if ((**m).isVertex()) {
325  haveVertex = true;
327 
328  // set pointers into big matrix
329  (**m).derivative(&m_fitMatrix.derivative[row][0]);
330  (**m).residual(row + m_residuals.begin());
331  ++row;
332  if ((**m).is2Dimensional()) {
334  (**m).derivative2(&m_fitMatrix.derivative[row][0]);
335  ++row;
336  }
337  ++m;
338  }
339 
340  // allocate rows to fitted measurements (DoF > 0)
341  for (; m != measurements.end(); ++m) {
342  if (!(**m).numberDoF())
343  continue;
344 
345  // alignment rows come after scattering
346  if ((**m).isAlignment())
347  continue;
348 
349  // identify leading material
350  if (!haveMeasurement) {
351  if ((**m).isPositionMeasurement()) {
352  haveMeasurement = true;
353  for (int i = 0; i < numberParameters; ++i)
354  if (m_firstRowForParameter[i] < 0)
356  } else if (!haveVertex && (**m).isScatterer()) {
357  (**m).numberDoF(0);
358  continue;
359  } else {
360  // row += (**m).numberDoF();
361  }
362  }
363 
364  // only allocate rows to fitted measurements (DoF > 0)
365  // if (! (**m).numberDoF()) continue;
366  if ((**m).isDrift())
368 
369  // fit energyDeposit unless momentum fixed or near infinite
370  if ((**m).isEnergyDeposit()) {
371  // if (! m_parameters->fitMomentum() ||
372  // m_parameters->extremeMomentum())
373  if (!m_parameters->fitMomentum()) {
374  (**m).numberDoF(0);
375  continue;
376  } else {
377  if (m_firstRowForParameter.back() < 0)
378  m_firstRowForParameter.back() = row;
379 
380  // m_firstRowForParameter[numberParameters-1] = row;
381  m_firstRowForParameter.push_back(row);
382  m_parameters->fitEnergyDeposit((**m).minEnergyDeposit());
383  ++numberEnergyDeposits;
384  ++numberParameters;
385  }
386  }
387  if ((**m).isScatterer())
388  ++numberScatterers;
389 
390  // set pointers into big matrix
391  (**m).derivative(&m_fitMatrix.derivative[row][0]);
392  (**m).residual(row + m_residuals.begin());
393  ++row;
394  if ((**m).is2Dimensional()) {
395  (**m).derivative2(&m_fitMatrix.derivative[row][0]);
396  ++row;
397  }
398  }
399 
400  // second loop puts alignment rows at bottom of matrix
401  for (m = measurements.begin(); m != measurements.end(); ++m) {
402  if (!(**m).isAlignment())
403  continue;
404  ++numberAlignments;
405 
406  // set pointers into big matrix
407  (**m).derivative(&m_fitMatrix.derivative[row][0]);
408  (**m).residual(row + m_residuals.begin());
409  ++row;
410  if ((**m).is2Dimensional()) {
411  (**m).derivative2(&m_fitMatrix.derivative[row][0]);
412  ++row;
413  }
414  }
415 
416  // keep first row with measurements up to each number of parameters
417  parameters->numberAlignments(numberAlignments);
418  parameters->numberScatterers(numberScatterers);
419  bool afterCalo = false;
420  int lastRow = 0;
421  m_rowsDM = 0;
422  for (m = measurements.begin(); m != measurements.end(); ++m) {
423  if ((**m).numberDoF()) {
424  if ((**m).isEnergyDeposit()) {
425  afterCalo = true;
426  // m_lastRowForParameter[4] = m_firstRowForParameter[5] + 1;
427  } else if ((**m).isScatterer()) {
428  m_firstRowForParameter.push_back(m_rowsDM);
429  m_firstRowForParameter.push_back(++m_rowsDM);
430  (**m).lastParameter(m_firstRowForParameter.size(), afterCalo);
431  parameters->addScatterer((**m).scattererPhi(), (**m).scattererTheta());
432  lastRow = ++m_rowsDM;
433  continue;
434  } else if ((**m).isAlignment()) {
435  // defer alignment
436  continue;
437  }
438  m_rowsDM += (**m).numberDoF();
439  lastRow = m_rowsDM;
440  }
441  (**m).lastParameter(m_firstRowForParameter.size(), afterCalo);
442  }
443 
444  numberParameters = m_firstRowForParameter.size();
445  m_lastRowForParameter = std::vector<int>(numberParameters, lastRow);
446  if (afterCalo) {
447  // cppcheck-suppress containerOutOfBounds; false positive
449  }
450 
451  // following loop puts any alignment info into the final rows
452  if (numberAlignments) {
453  row = 0;
454  unsigned alignmentParameter = 0;
455  int firstAlignmentRow = 0;
456  int firstAlignment2Row = 0;
457  for (m = measurements.begin(); m != measurements.end(); ++m) {
458  if ((**m).alignmentParameter()) {
459  if ((**m).alignmentParameter() > alignmentParameter) {
460  if (!firstAlignmentRow)
461  firstAlignmentRow = row;
462  }
463  if ((**m).alignmentParameter2()) {
464  if (!firstAlignment2Row)
465  firstAlignment2Row = row;
466  }
467  }
468  if (!(**m).numberDoF())
469  continue;
470  if (!(**m).isAlignment()) {
471  row += (**m).numberDoF();
472  continue;
473  }
474 
476  (**m).alignmentAngle(), (**m).alignmentOffset());
477  m_firstRowForParameter.push_back(firstAlignmentRow);
478  m_firstRowForParameter.push_back(firstAlignmentRow);
479  m_lastRowForParameter.push_back(++m_rowsDM);
480  m_lastRowForParameter.push_back(++m_rowsDM);
481  (**m).firstParameter(numberParameters);
482  numberParameters += 2;
483  alignmentParameter = parameters->numberAlignments();
484  (**m).alignmentParameter(alignmentParameter);
485  (**m).lastParameter(numberParameters, afterCalo);
486  // m_rowsDM += (**m).numberDoF();
487 
488  // some bug to fix here...
489  // firstAlignmentRow = firstAlignment2Row;
490  // firstAlignment2Row = 0;
491  }
492  }
493 
494  // initialize number of parameters (including alignments and scatterers)
495  numberParameters = m_firstRowForParameter.size();
496  parameters->numberParameters(numberParameters);
497 
498  // and degrees of freedom
499  m_numberDoF = m_rowsDM - numberParameters;
500 
501  // make some checks: return fitCode in case of problem
502  int fitCode = 0;
503  // if (row > mxmeas) fitCode = 2; // too many measurements for
504  // fit matrix size
505  if (2 * measurements.size() > mxmeas)
506  fitCode = 2; // too many measurements for fit matrix size
507  if (numberParameters > mxparam)
508  fitCode = 3; // too many parameters for fit matrix size
509  if (m_numberDoF < 0)
510  fitCode = 4; // unconstrained fit: negative numberDoF
511  if (numberEnergyDeposits > 1)
512  fitCode = 12; // too many EnergyDeposit parameters
513  if (fitCode)
514  return fitCode;
515 
516  // reserve derivatives for jacobian propagation to 'non-measurements'
517  row = m_rowsDM;
518  for (m = measurements.begin(); m != measurements.end(); ++m) {
519  if (!(**m).isPositionMeasurement() || (**m).numberDoF() > 1)
520  continue;
521  if (!(**m).numberDoF()) {
522  for (int i = 0; i != numberParameters; ++i)
523  m_fitMatrix.derivative[row][i] = 0.;
524  (**m).derivative(&m_fitMatrix.derivative[row][0]);
525  (**m).residual(row + m_residuals.begin());
526  ++row;
527  }
528  for (int i = 0; i != numberParameters; ++i)
529  m_fitMatrix.derivative[row][i] = 0.;
530  (**m).derivative2(&m_fitMatrix.derivative[row][0]);
531  ++row;
532  }
533 
534  // update partitioning of fit matrices
535  for (int row = 0; row < m_rowsDM; ++row) {
536  for (int param = 0; param < numberParameters; ++param) {
537  m_fitMatrix.derivative[row][param] = 0.;
538  }
539  }
540 
541  // fix degrees of freedom for external customers
542  if (!m_parameters->fitMomentum())
543  ++m_numberDoF;
544 
545  // we don't have any fit results yet
546  m_covariance.resize(0,0);
547  m_finalCovariance.resize(0,0);
548 
549  // reallocate to get correct matrix sizes
550  if (m_derivativeMatrix.size() == 0 || m_weight.size() == 0 ||
551  numberParameters != m_columnsDM) {
552  m_columnsDM = numberParameters;
555  // isn't this faster? indicating that we have a symmetric matrix <0x2> =
556  // <Upper> any gain seems to be negated by additional for loop copies to
557  // recover full cov matrix introduces some rounding differences - keep for
558  // release 21 to respect strict Tier0 policy
559  (m_weight).selfadjointView<0x2>();
561  }
562 
563  // this should never happen
564  if (m_derivativeMatrix.size()==0)
565  fitCode = 13;
566 
567  return fitCode;
568 }

◆ solveEquations()

ATH_FLATTEN bool Trk::FitMatrices::solveEquations ( void  )

Definition at line 577 of file FitMatrices.cxx.

577  {
578  // use Eigen Matrix multiplication ATR-15723
579  // note: fitMatrix (struct) is row-major whereas Eigen prefers column-major
580  // storage
581  // hence the loops are nested to optimise row-major and to form the
582  // Eigen transpose
584  Amg::VectorX& weightedDifference = m_weightedDifference;
585  Amg::MatrixX fitMatrixDerivativeT(m_columnsDM, m_rowsDM);
586  Amg::MatrixX residuals(m_rowsDM, 1);
587  for (int row = 0; row < m_rowsDM; ++row) {
588  residuals(row, 0) = (m_residuals)[row];
589  for (int col = 0; col < m_columnsDM; ++col) {
590  fitMatrixDerivativeT(col, row) = m_fitMatrix.derivative[row][col];
591  }
592  }
593  weight = fitMatrixDerivativeT * fitMatrixDerivativeT.transpose();
594  weightedDifference = fitMatrixDerivativeT * residuals;
595 
596  // stabilize fits with badly determined phi
598  weight(0, 0) += m_largePhiWeight;
599 
600  // avoid some possible singularities in matrix inversion
602 
603  // solve is faster than inverse: wait for explicit request for covariance
604  // before inversion
606  weight.colPivHouseholderQr().solve(weightedDifference);
607 
609  return true;
610 }

◆ usePerigee()

void Trk::FitMatrices::usePerigee ( const FitMeasurement measurement)

Definition at line 612 of file FitMatrices.cxx.

612  {
613  m_perigee = &(measurement.perigee());
614  m_perigeeWeight = &(measurement.perigeeWeight());
615  // TODO: needs eigen equiv !!
616  if (m_matrixFromCLHEP) {
617  m_usePerigee = true;
618  }
619 }

Member Data Documentation

◆ m_columnsDM

int Trk::FitMatrices::m_columnsDM
private

Definition at line 116 of file FitMatrices.h.

◆ m_constrainedAlignmentEffects

bool Trk::FitMatrices::m_constrainedAlignmentEffects
private

Definition at line 117 of file FitMatrices.h.

◆ m_covariance

Amg::MatrixX Trk::FitMatrices::m_covariance {}
private

Definition at line 118 of file FitMatrices.h.

◆ m_derivativeMatrix

Amg::MatrixX Trk::FitMatrices::m_derivativeMatrix {}
private

Definition at line 119 of file FitMatrices.h.

◆ m_finalCovariance

Amg::MatrixX Trk::FitMatrices::m_finalCovariance {}
private

Definition at line 120 of file FitMatrices.h.

◆ m_firstRowForParameter

std::vector<int> Trk::FitMatrices::m_firstRowForParameter
private

Definition at line 121 of file FitMatrices.h.

◆ m_fitMatrix

fitMatrix Trk::FitMatrices::m_fitMatrix
private

Definition at line 115 of file FitMatrices.h.

◆ m_largePhiWeight

double Trk::FitMatrices::m_largePhiWeight
private

Definition at line 122 of file FitMatrices.h.

◆ m_lastRowForParameter

std::vector<int> Trk::FitMatrices::m_lastRowForParameter
private

Definition at line 123 of file FitMatrices.h.

◆ m_matrixFromCLHEP

bool Trk::FitMatrices::m_matrixFromCLHEP
private

Definition at line 124 of file FitMatrices.h.

◆ m_measurements

std::vector<FitMeasurement*>* Trk::FitMatrices::m_measurements
private

Definition at line 125 of file FitMatrices.h.

◆ m_numberDoF

int Trk::FitMatrices::m_numberDoF
private

Definition at line 126 of file FitMatrices.h.

◆ m_numberDriftCircles

int Trk::FitMatrices::m_numberDriftCircles
private

Definition at line 127 of file FitMatrices.h.

◆ m_numberPerigee

int Trk::FitMatrices::m_numberPerigee
private

Definition at line 128 of file FitMatrices.h.

◆ m_parameters

FitParameters* Trk::FitMatrices::m_parameters
private

Definition at line 129 of file FitMatrices.h.

◆ m_perigee

const Amg::VectorX* Trk::FitMatrices::m_perigee
private

Definition at line 130 of file FitMatrices.h.

◆ m_perigeeDifference

Amg::MatrixX Trk::FitMatrices::m_perigeeDifference
private

Definition at line 131 of file FitMatrices.h.

◆ m_perigeeWeight

const Amg::MatrixX* Trk::FitMatrices::m_perigeeWeight
private

Definition at line 132 of file FitMatrices.h.

◆ m_residuals

std::vector<double> Trk::FitMatrices::m_residuals
private

Definition at line 133 of file FitMatrices.h.

◆ m_rowsDM

int Trk::FitMatrices::m_rowsDM
private

Definition at line 134 of file FitMatrices.h.

◆ m_usePerigee

bool Trk::FitMatrices::m_usePerigee
private

Definition at line 135 of file FitMatrices.h.

◆ m_weight

Amg::MatrixX Trk::FitMatrices::m_weight
private

Definition at line 136 of file FitMatrices.h.

◆ m_weightedDifference

Amg::VectorX Trk::FitMatrices::m_weightedDifference
private

Definition at line 137 of file FitMatrices.h.


The documentation for this class was generated from the following files:
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
query_example.row
row
Definition: query_example.py:24
Trk::FitMatrices::m_weightedDifference
Amg::VectorX m_weightedDifference
Definition: FitMatrices.h:137
Amg::VectorX
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorX
Dynamic Vector - dynamic allocation.
Definition: EventPrimitives.h:30
Amg::hasPositiveDiagElems
bool hasPositiveDiagElems(const AmgSymMatrix(N) &mat)
Returns true if all diagonal elements of the covariance matrix are finite aka sane in the above defin...
Definition: EventPrimitivesCovarianceHelpers.h:96
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition: EventPrimitives.h:27
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
Trk::FitParameters::fitMomentum
bool fitMomentum(void) const
Definition: FitParameters.h:225
Trk::FitMatrices::m_perigeeDifference
Amg::MatrixX m_perigeeDifference
Definition: FitMatrices.h:131
Trk::FitMatrices::avoidMomentumSingularity
void avoidMomentumSingularity(void)
Definition: FitMatrices.cxx:636
Trk::FitMatrices::m_weight
Amg::MatrixX m_weight
Definition: FitMatrices.h:136
Trk::FitParameters::fitEnergyDeposit
bool fitEnergyDeposit(void) const
Definition: FitParameters.h:221
Trk::FitMatrices::m_derivativeMatrix
Amg::MatrixX m_derivativeMatrix
Definition: FitMatrices.h:119
Trk::FitMatrices::m_rowsDM
int m_rowsDM
Definition: FitMatrices.h:134
python.SystemOfUnits.TeV
int TeV
Definition: SystemOfUnits.py:158
Trk::FitMatrices::m_perigee
const Amg::VectorX * m_perigee
Definition: FitMatrices.h:130
Trk::FitMatrices::m_fitMatrix
fitMatrix m_fitMatrix
Definition: FitMatrices.h:115
Trk::FitMatrices::m_numberDriftCircles
int m_numberDriftCircles
Definition: FitMatrices.h:127
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:189
Trk::fitMatrix::derivative
double derivative[mxmeas][mxparam]
Definition: FitMatrix.h:14
Trk::FitMatrices::m_lastRowForParameter
std::vector< int > m_lastRowForParameter
Definition: FitMatrices.h:123
mxparam
#define mxparam
Definition: FitMatrix.h:9
lumiFormat.i
int i
Definition: lumiFormat.py:85
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
MuonR4::inverse
CalibratedSpacePoint::Covariance_t inverse(const CalibratedSpacePoint::Covariance_t &mat)
Inverts the parsed matrix.
Definition: MuonSpectrometer/MuonPhaseII/Event/MuonSpacePoint/src/UtilFunctions.cxx:65
Trk::FitMatrices::m_usePerigee
bool m_usePerigee
Definition: FitMatrices.h:135
Trk::FitParameters::update
void update(const Amg::VectorX &differences)
Definition: FitParameters.cxx:603
Trk::FitMatrices::m_finalCovariance
Amg::MatrixX m_finalCovariance
Definition: FitMatrices.h:120
Trk::FitMatrices::m_matrixFromCLHEP
bool m_matrixFromCLHEP
Definition: FitMatrices.h:124
Trk::FitMatrices::m_numberPerigee
int m_numberPerigee
Definition: FitMatrices.h:128
Trk::FitMatrices::m_residuals
std::vector< double > m_residuals
Definition: FitMatrices.h:133
Trk::FitParameters::parameterDifference
const Amg::MatrixX parameterDifference(const Amg::VectorX &parameters) const
Definition: FitParameters.cxx:186
Trk::FitMatrices::m_largePhiWeight
double m_largePhiWeight
Definition: FitMatrices.h:122
Trk::FitMatrices::m_covariance
Amg::MatrixX m_covariance
Definition: FitMatrices.h:118
query_example.col
col
Definition: query_example.py:7
Trk::FitMatrices::m_parameters
FitParameters * m_parameters
Definition: FitMatrices.h:129
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
physics_parameters.parameters
parameters
Definition: physics_parameters.py:144
Trk::FitMatrices::m_columnsDM
int m_columnsDM
Definition: FitMatrices.h:116
Trk::FitMatrices::m_constrainedAlignmentEffects
bool m_constrainedAlignmentEffects
Definition: FitMatrices.h:117
Trk::FitMatrices::m_measurements
std::vector< FitMeasurement * > * m_measurements
Definition: FitMatrices.h:125
Trk::FitMatrices::m_numberDoF
int m_numberDoF
Definition: FitMatrices.h:126
Trk::FitParameters::phiInstability
bool phiInstability(void) const
Definition: FitParameters.cxx:244
Trk::FitMatrices::m_firstRowForParameter
std::vector< int > m_firstRowForParameter
Definition: FitMatrices.h:121
Trk::FitMatrices::solveEquations
bool solveEquations(void)
Definition: FitMatrices.cxx:577
mxmeas
#define mxmeas
Definition: FitMatrix.h:8
Trk::FitMatrices::m_perigeeWeight
const Amg::MatrixX * m_perigeeWeight
Definition: FitMatrices.h:132