ATLAS Offline Software
Loading...
Searching...
No Matches
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),
51 m_parameters(nullptr),
52 m_perigee(nullptr),
54 m_perigeeWeight(nullptr),
55 m_rowsDM(0),
56 m_usePerigee(false){}
Amg::MatrixX m_perigeeDifference
std::vector< FitMeasurement * > * m_measurements
const Amg::MatrixX * m_perigeeWeight
bool m_constrainedAlignmentEffects
const Amg::VectorX * m_perigee
double m_largePhiWeight
FitParameters * m_parameters
fitMatrix m_fitMatrix
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.

◆ 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
639 if (m_parameters->fitEnergyDeposit() &&
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}
Amg::MatrixX m_weight

◆ 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 const 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}
#define endmsg
std::vector< int > m_firstRowForParameter
std::vector< int > m_lastRowForParameter

◆ 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 ( void )
inline

Definition at line 142 of file FitMatrices.h.

142 {
143 return m_derivativeMatrix;
144}
Amg::MatrixX m_derivativeMatrix

◆ finalCovariance()

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

Definition at line 146 of file FitMatrices.h.

146 {
147 return m_finalCovariance;
148}
Amg::MatrixX m_finalCovariance

◆ 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?
96 if (m_parameters->phiInstability()) {
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 const 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)
138 if (m_parameters->fitEnergyDeposit()) {
139 const 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}
Amg::MatrixX m_covariance
void avoidMomentumSingularity(void)
bool solveEquations(void)
bool hasPositiveDiagElems(const AmgSymMatrix(N) &mat)
Returns true if all diagonal elements of the covariance matrix are finite aka sane in the above defin...
row
Appending html table to final .html summary file.

◆ 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 {
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 {
163 m_perigeeDifference = m_parameters->parameterDifference(*m_perigee);
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;
176 std::vector<FitMeasurement*>::iterator m = m_measurements->begin();
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) {
276 int i = m_firstRowForParameter[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}
Amg::VectorX m_weightedDifference

◆ 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
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)
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()) {
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}
#define mxparam
Definition FitMatrix.h:9
#define mxmeas
Definition FitMatrix.h:8
std::vector< double > m_residuals
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorX
Dynamic Vector - dynamic allocation.

◆ 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
597 if (m_parameters->phiInstability())
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.

118{};

◆ m_derivativeMatrix

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

Definition at line 119 of file FitMatrices.h.

119{};

◆ m_finalCovariance

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

Definition at line 120 of file FitMatrices.h.

120{};

◆ 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: