Skip to content
Snippets Groups Projects
Commit f0ae26d7 authored by Christos Anastopoulos's avatar Christos Anastopoulos Committed by Vakhtang Tsulaia
Browse files

FitMatrices: Avoid ptr to DynamicEigen/std::vector get rid of extraneous new/delete

FitMatrices: Avoid ptr to DynamicEigen/std::vector get rid of extraneous new/delete
parent 558f5b45
No related branches found
No related tags found
1 merge request!66690FitMatrices: Avoid ptr to DynamicEigen/std::vector get rid of extraneous new/delete
/* /*
Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
*/ */
////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////
...@@ -42,10 +42,14 @@ class FitMatrices { ...@@ -42,10 +42,14 @@ class FitMatrices {
public: public:
FitMatrices(bool constrainedAlignmentEffects); FitMatrices(bool constrainedAlignmentEffects);
~FitMatrices(void); // copy, assignment deleted
FitMatrices(const FitMatrices&) = delete;
// forbidden copy constructor FitMatrices& operator=(const FitMatrices&) = delete;
// forbidden assignment operator //move defaulted
FitMatrices(FitMatrices&&) = default;
FitMatrices& operator=(FitMatrices&&) = default;
//dtor default
~FitMatrices() = default ;
// debugging aid: check 'smart' pointers // debugging aid: check 'smart' pointers
void checkPointers(MsgStream& log) const; void checkPointers(MsgStream& log) const;
...@@ -54,13 +58,13 @@ class FitMatrices { ...@@ -54,13 +58,13 @@ class FitMatrices {
static double chiSquaredChange(void); static double chiSquaredChange(void);
// accessor to DerivativeMatrix (eigen) // accessor to DerivativeMatrix (eigen)
Amg::MatrixX* derivativeMatrix(void); Amg::MatrixX& derivativeMatrix();
// accessor to final covariance (5*5) // accessor to final covariance (5*5)
// includes leading material effects on construction, but // includes leading material effects on construction, but
// field gradient effects have to be added by MeasurementProcessor as // field gradient effects have to be added by MeasurementProcessor as
// 'extrapolation aware' // 'extrapolation aware'
Amg::MatrixX* finalCovariance(void); Amg::MatrixX& finalCovariance();
// produce full covariance (including scatterer parameters) i.e. invert weight // produce full covariance (including scatterer parameters) i.e. invert weight
// matrix NOTE: this does not contain the external errors which will be // matrix NOTE: this does not contain the external errors which will be
...@@ -103,11 +107,7 @@ class FitMatrices { ...@@ -103,11 +107,7 @@ class FitMatrices {
void usePerigee(const FitMeasurement& measurement); void usePerigee(const FitMeasurement& measurement);
private: private:
// copy, assignment: no semantics, no implementation // add perigee measurement
FitMatrices(const FitMatrices&) = delete;
FitMatrices& operator=(const FitMatrices&) = delete;
// add perigee measurement
void addPerigeeMeasurement(void); void addPerigeeMeasurement(void);
// fix for momentum singularity // fix for momentum singularity
void avoidMomentumSingularity(void); // using Eigen void avoidMomentumSingularity(void); // using Eigen
...@@ -115,35 +115,35 @@ class FitMatrices { ...@@ -115,35 +115,35 @@ class FitMatrices {
fitMatrix m_fitMatrix; fitMatrix m_fitMatrix;
int m_columnsDM; int m_columnsDM;
bool m_constrainedAlignmentEffects; bool m_constrainedAlignmentEffects;
Amg::MatrixX* m_covariance; Amg::MatrixX m_covariance{};
Amg::MatrixX* m_derivativeMatrix; Amg::MatrixX m_derivativeMatrix{};
Amg::MatrixX* m_finalCovariance; Amg::MatrixX m_finalCovariance{};
std::vector<int> m_firstRowForParameter; std::vector<int> m_firstRowForParameter;
double m_largePhiWeight; double m_largePhiWeight;
std::vector<int> m_lastRowForParameter; std::vector<int> m_lastRowForParameter;
bool m_matrixFromCLHEP; bool m_matrixFromCLHEP;
std::vector<FitMeasurement*>* m_measurements; std::vector<FitMeasurement*>* m_measurements; // not owning ptr
int m_numberDoF; int m_numberDoF;
int m_numberDriftCircles; int m_numberDriftCircles;
int m_numberPerigee; int m_numberPerigee;
FitParameters* m_parameters; FitParameters* m_parameters;
const Amg::VectorX* m_perigee; const Amg::VectorX* m_perigee; // not owning ptr
Amg::MatrixX m_perigeeDifference; Amg::MatrixX m_perigeeDifference;
const Amg::MatrixX* m_perigeeWeight; const Amg::MatrixX* m_perigeeWeight; // not owning ptr
std::vector<double>* m_residuals; std::vector<double> m_residuals;
int m_rowsDM; int m_rowsDM;
bool m_usePerigee; bool m_usePerigee;
Amg::MatrixX* m_weight; Amg::MatrixX m_weight;
Amg::VectorX* m_weightedDifference; Amg::VectorX m_weightedDifference;
}; };
//<<<<<< INLINE MEMBER FUNCTIONS >>>>>> //<<<<<< INLINE MEMBER FUNCTIONS >>>>>>
inline Amg::MatrixX* FitMatrices::derivativeMatrix(void) { inline Amg::MatrixX& FitMatrices::derivativeMatrix(void) {
return m_derivativeMatrix; return m_derivativeMatrix;
} }
inline Amg::MatrixX* FitMatrices::finalCovariance(void) { inline Amg::MatrixX& FitMatrices::finalCovariance(void) {
return m_finalCovariance; return m_finalCovariance;
} }
......
...@@ -42,9 +42,6 @@ FitMatrices::FitMatrices(bool constrainedAlignmentEffects) ...@@ -42,9 +42,6 @@ FitMatrices::FitMatrices(bool constrainedAlignmentEffects)
: m_fitMatrix{}, : m_fitMatrix{},
m_columnsDM(0), m_columnsDM(0),
m_constrainedAlignmentEffects(constrainedAlignmentEffects), m_constrainedAlignmentEffects(constrainedAlignmentEffects),
m_covariance(nullptr),
m_derivativeMatrix(nullptr),
m_finalCovariance(nullptr),
m_largePhiWeight(10000.), // arbitrary - equiv to 10um m_largePhiWeight(10000.), // arbitrary - equiv to 10um
m_matrixFromCLHEP(false), m_matrixFromCLHEP(false),
m_measurements(nullptr), m_measurements(nullptr),
...@@ -55,20 +52,8 @@ FitMatrices::FitMatrices(bool constrainedAlignmentEffects) ...@@ -55,20 +52,8 @@ FitMatrices::FitMatrices(bool constrainedAlignmentEffects)
m_perigee(nullptr), m_perigee(nullptr),
m_perigeeDifference(Amg::MatrixX(1, m_numberPerigee)), m_perigeeDifference(Amg::MatrixX(1, m_numberPerigee)),
m_perigeeWeight(nullptr), m_perigeeWeight(nullptr),
m_residuals(nullptr),
m_rowsDM(0), m_rowsDM(0),
m_usePerigee(false), m_usePerigee(false){}
m_weight(nullptr),
m_weightedDifference(nullptr) {}
FitMatrices::~FitMatrices(void) {
delete m_covariance;
delete m_derivativeMatrix;
delete m_finalCovariance;
delete m_residuals;
delete m_weight;
delete m_weightedDifference;
}
void FitMatrices::checkPointers(MsgStream& log) const { void FitMatrices::checkPointers(MsgStream& log) const {
// debugging: check smart pointers // debugging: check smart pointers
...@@ -102,10 +87,10 @@ double FitMatrices::chiSquaredChange(void) { ...@@ -102,10 +87,10 @@ double FitMatrices::chiSquaredChange(void) {
const Amg::MatrixX* FitMatrices::fullCovariance(void) { const Amg::MatrixX* FitMatrices::fullCovariance(void) {
// return result if matrix already inverted // return result if matrix already inverted
if (m_covariance) { if (m_covariance.size()!=0) {
return m_covariance; return &m_covariance;
} }
m_covariance = new Amg::MatrixX(m_columnsDM, m_columnsDM); m_covariance = Amg::MatrixX(m_columnsDM, m_columnsDM);
// fix weighting ???? shouldn't we just remove large phi weight? // fix weighting ???? shouldn't we just remove large phi weight?
if (m_parameters->phiInstability()) { if (m_parameters->phiInstability()) {
...@@ -113,7 +98,7 @@ const Amg::MatrixX* FitMatrices::fullCovariance(void) { ...@@ -113,7 +98,7 @@ const Amg::MatrixX* FitMatrices::fullCovariance(void) {
} }
// invert weight matrix // invert weight matrix
Amg::MatrixX& covariance = *m_covariance; Amg::MatrixX& covariance = m_covariance;
// avoid singularity through ill-defined momentum ???? again // avoid singularity through ill-defined momentum ???? again
avoidMomentumSingularity(); avoidMomentumSingularity();
...@@ -127,13 +112,12 @@ const Amg::MatrixX* FitMatrices::fullCovariance(void) { ...@@ -127,13 +112,12 @@ const Amg::MatrixX* FitMatrices::fullCovariance(void) {
weight.selfadjointView<0x2>(); weight.selfadjointView<0x2>();
// check if m_weights makes sense before inverting // check if m_weights makes sense before inverting
if (!Amg::saneCovarianceDiagonal(*m_weight)) { if (!Amg::saneCovarianceDiagonal(m_weight)) {
delete m_covariance; m_covariance.resize(0, 0);
m_covariance = nullptr;
return nullptr; return nullptr;
} }
weight = (*m_weight).inverse(); weight = (m_weight).inverse();
for (int row = 0; row != m_columnsDM; ++row) { for (int row = 0; row != m_columnsDM; ++row) {
for (int col = 0; col != m_columnsDM; ++col) for (int col = 0; col != m_columnsDM; ++col)
covariance(row, col) = weight(col, row); covariance(row, col) = weight(col, row);
...@@ -164,16 +148,15 @@ const Amg::MatrixX* FitMatrices::fullCovariance(void) { ...@@ -164,16 +148,15 @@ const Amg::MatrixX* FitMatrices::fullCovariance(void) {
// FIXME: errors underestimated on d0,z0 when large scatterer precedes precise // FIXME: errors underestimated on d0,z0 when large scatterer precedes precise
// measurements final covariance starts with 5*5 representing perigee from // measurements final covariance starts with 5*5 representing perigee from
// full covariance // full covariance
delete m_finalCovariance; m_finalCovariance = Amg::MatrixX(5, 5);
m_finalCovariance = new Amg::MatrixX(5, 5);
for (int i = 0; i != 5; ++i) { for (int i = 0; i != 5; ++i) {
for (int j = 0; j != 5; ++j) { for (int j = 0; j != 5; ++j) {
(*m_finalCovariance)(i, j) = covariance(i, j); (m_finalCovariance)(i, j) = covariance(i, j);
} }
} }
// return pointer to full covariance // return pointer to full covariance
return m_covariance; return &m_covariance;
} }
double FitMatrices::perigeeChiSquared(void) { double FitMatrices::perigeeChiSquared(void) {
...@@ -276,7 +259,7 @@ void FitMatrices::printWeightMatrix(void) { ...@@ -276,7 +259,7 @@ void FitMatrices::printWeightMatrix(void) {
<< row << " col 0 "; << row << " col 0 ";
for (int col = 0; col <= row; ++col) { for (int col = 0; col <= row; ++col) {
std::cout << std::setiosflags(std::ios::scientific) << std::setbase(10) std::cout << std::setiosflags(std::ios::scientific) << std::setbase(10)
<< std::setw(10) << (*m_weight)(row, col) << " "; << std::setw(10) << (m_weight)(row, col) << " ";
if ((col + 1) % 13 == 0 && col < row) if ((col + 1) % 13 == 0 && col < row)
std::cout << std::endl std::cout << std::endl
...@@ -306,12 +289,9 @@ void FitMatrices::refinePointers(void) { ...@@ -306,12 +289,9 @@ void FitMatrices::refinePointers(void) {
} }
void FitMatrices::releaseMemory(void) { void FitMatrices::releaseMemory(void) {
delete m_derivativeMatrix; m_derivativeMatrix.resize(0,0);
delete m_weight; m_weight.resize(0,0);
delete m_weightedDifference; m_weightedDifference.resize(0,0);
m_derivativeMatrix = nullptr;
m_weight = nullptr;
m_weightedDifference = nullptr;
} }
int FitMatrices::setDimensions(std::vector<FitMeasurement*>& measurements, int FitMatrices::setDimensions(std::vector<FitMeasurement*>& measurements,
...@@ -328,8 +308,7 @@ int FitMatrices::setDimensions(std::vector<FitMeasurement*>& measurements, ...@@ -328,8 +308,7 @@ int FitMatrices::setDimensions(std::vector<FitMeasurement*>& measurements,
m_lastRowForParameter.clear(); m_lastRowForParameter.clear();
m_lastRowForParameter.reserve(128); m_lastRowForParameter.reserve(128);
m_parameters = parameters; m_parameters = parameters;
delete m_residuals; m_residuals = std::vector<double>(2 * measurements.size(), 0.);
m_residuals = new std::vector<double>(2 * measurements.size(), 0.);
m_numberDriftCircles = 0; m_numberDriftCircles = 0;
bool haveMeasurement = false; bool haveMeasurement = false;
bool haveVertex = false; bool haveVertex = false;
...@@ -348,7 +327,7 @@ int FitMatrices::setDimensions(std::vector<FitMeasurement*>& measurements, ...@@ -348,7 +327,7 @@ int FitMatrices::setDimensions(std::vector<FitMeasurement*>& measurements,
// set pointers into big matrix // set pointers into big matrix
(**m).derivative(&m_fitMatrix.derivative[row][0]); (**m).derivative(&m_fitMatrix.derivative[row][0]);
(**m).residual(row + m_residuals->begin()); (**m).residual(row + m_residuals.begin());
++row; ++row;
if ((**m).is2Dimensional()) { if ((**m).is2Dimensional()) {
m_firstRowForParameter[row] = row; m_firstRowForParameter[row] = row;
...@@ -410,7 +389,7 @@ int FitMatrices::setDimensions(std::vector<FitMeasurement*>& measurements, ...@@ -410,7 +389,7 @@ int FitMatrices::setDimensions(std::vector<FitMeasurement*>& measurements,
// set pointers into big matrix // set pointers into big matrix
(**m).derivative(&m_fitMatrix.derivative[row][0]); (**m).derivative(&m_fitMatrix.derivative[row][0]);
(**m).residual(row + m_residuals->begin()); (**m).residual(row + m_residuals.begin());
++row; ++row;
if ((**m).is2Dimensional()) { if ((**m).is2Dimensional()) {
(**m).derivative2(&m_fitMatrix.derivative[row][0]); (**m).derivative2(&m_fitMatrix.derivative[row][0]);
...@@ -426,7 +405,7 @@ int FitMatrices::setDimensions(std::vector<FitMeasurement*>& measurements, ...@@ -426,7 +405,7 @@ int FitMatrices::setDimensions(std::vector<FitMeasurement*>& measurements,
// set pointers into big matrix // set pointers into big matrix
(**m).derivative(&m_fitMatrix.derivative[row][0]); (**m).derivative(&m_fitMatrix.derivative[row][0]);
(**m).residual(row + m_residuals->begin()); (**m).residual(row + m_residuals.begin());
++row; ++row;
if ((**m).is2Dimensional()) { if ((**m).is2Dimensional()) {
(**m).derivative2(&m_fitMatrix.derivative[row][0]); (**m).derivative2(&m_fitMatrix.derivative[row][0]);
...@@ -541,7 +520,7 @@ int FitMatrices::setDimensions(std::vector<FitMeasurement*>& measurements, ...@@ -541,7 +520,7 @@ int FitMatrices::setDimensions(std::vector<FitMeasurement*>& measurements,
for (int i = 0; i != numberParameters; ++i) for (int i = 0; i != numberParameters; ++i)
m_fitMatrix.derivative[row][i] = 0.; m_fitMatrix.derivative[row][i] = 0.;
(**m).derivative(&m_fitMatrix.derivative[row][0]); (**m).derivative(&m_fitMatrix.derivative[row][0]);
(**m).residual(row + m_residuals->begin()); (**m).residual(row + m_residuals.begin());
++row; ++row;
} }
for (int i = 0; i != numberParameters; ++i) for (int i = 0; i != numberParameters; ++i)
...@@ -562,30 +541,25 @@ int FitMatrices::setDimensions(std::vector<FitMeasurement*>& measurements, ...@@ -562,30 +541,25 @@ int FitMatrices::setDimensions(std::vector<FitMeasurement*>& measurements,
++m_numberDoF; ++m_numberDoF;
// we don't have any fit results yet // we don't have any fit results yet
delete m_covariance; m_covariance.resize(0,0);
m_covariance = nullptr; m_finalCovariance.resize(0,0);
delete m_finalCovariance;
m_finalCovariance = nullptr;
// reallocate to get correct matrix sizes // reallocate to get correct matrix sizes
if (!m_derivativeMatrix || !m_weight || numberParameters != m_columnsDM) { if (m_derivativeMatrix.size() == 0 || m_weight.size() == 0 ||
numberParameters != m_columnsDM) {
m_columnsDM = numberParameters; m_columnsDM = numberParameters;
delete m_derivativeMatrix; m_derivativeMatrix = Amg::MatrixX(m_rowsDM, m_columnsDM);
m_derivativeMatrix = new Amg::MatrixX(m_rowsDM, m_columnsDM); m_weight = Amg::MatrixX(m_columnsDM, m_columnsDM);
delete m_weight;
m_weight = new Amg::MatrixX(m_columnsDM, m_columnsDM);
// isn't this faster? indicating that we have a symmetric matrix <0x2> = // isn't this faster? indicating that we have a symmetric matrix <0x2> =
// <Upper> any gain seems to be negated by additional for loop copies to // <Upper> any gain seems to be negated by additional for loop copies to
// recover full cov matrix introduces some rounding differences - keep for // recover full cov matrix introduces some rounding differences - keep for
// release 21 to respect strict Tier0 policy // release 21 to respect strict Tier0 policy
(*m_weight).selfadjointView<0x2>(); (m_weight).selfadjointView<0x2>();
m_weightedDifference = Amg::VectorX(m_columnsDM);
delete m_weightedDifference;
m_weightedDifference = new Amg::VectorX(m_columnsDM);
} }
// this should never happen // this should never happen
if (!m_derivativeMatrix) if (m_derivativeMatrix.size()==0)
fitCode = 13; fitCode = 13;
return fitCode; return fitCode;
...@@ -604,12 +578,12 @@ ATH_FLATTEN ...@@ -604,12 +578,12 @@ ATH_FLATTEN
// storage // storage
// hence the loops are nested to optimise row-major and to form the // hence the loops are nested to optimise row-major and to form the
// Eigen transpose // Eigen transpose
Amg::MatrixX& weight = *m_weight; Amg::MatrixX& weight = m_weight;
Amg::VectorX& weightedDifference = *m_weightedDifference; Amg::VectorX& weightedDifference = m_weightedDifference;
Amg::MatrixX fitMatrixDerivativeT(m_columnsDM, m_rowsDM); Amg::MatrixX fitMatrixDerivativeT(m_columnsDM, m_rowsDM);
Amg::MatrixX residuals(m_rowsDM, 1); Amg::MatrixX residuals(m_rowsDM, 1);
for (int row = 0; row < m_rowsDM; ++row) { for (int row = 0; row < m_rowsDM; ++row) {
residuals(row, 0) = (*m_residuals)[row]; residuals(row, 0) = (m_residuals)[row];
for (int col = 0; col < m_columnsDM; ++col) { for (int col = 0; col < m_columnsDM; ++col) {
fitMatrixDerivativeT(col, row) = m_fitMatrix.derivative[row][col]; fitMatrixDerivativeT(col, row) = m_fitMatrix.derivative[row][col];
} }
...@@ -626,10 +600,10 @@ ATH_FLATTEN ...@@ -626,10 +600,10 @@ ATH_FLATTEN
// solve is faster than inverse: wait for explicit request for covariance // solve is faster than inverse: wait for explicit request for covariance
// before inversion // before inversion
*m_weightedDifference = m_weightedDifference =
weight.colPivHouseholderQr().solve(weightedDifference); weight.colPivHouseholderQr().solve(weightedDifference);
m_parameters->update(*m_weightedDifference); m_parameters->update(m_weightedDifference);
return true; return true;
} }
...@@ -659,7 +633,7 @@ void FitMatrices::addPerigeeMeasurement(void) { ...@@ -659,7 +633,7 @@ void FitMatrices::addPerigeeMeasurement(void) {
void FitMatrices::avoidMomentumSingularity(void) { void FitMatrices::avoidMomentumSingularity(void) {
// fix momentum if line-fit or fit attempted with negligible field integral // fix momentum if line-fit or fit attempted with negligible field integral
Amg::MatrixX& weight = *m_weight; Amg::MatrixX& weight = m_weight;
if (m_parameters->fitEnergyDeposit() && if (m_parameters->fitEnergyDeposit() &&
weight(5, 5) < 1. / Gaudi::Units::TeV) { weight(5, 5) < 1. / Gaudi::Units::TeV) {
for (int i = 0; i != m_columnsDM; ++i) { for (int i = 0; i != m_columnsDM; ++i) {
......
...@@ -370,7 +370,7 @@ const FitProcedureQuality& FitProcedure::execute( ...@@ -370,7 +370,7 @@ const FitProcedureQuality& FitProcedure::execute(
cache.numberParameters = parameters->numberParameters(); cache.numberParameters = parameters->numberParameters();
cache.worstMeasurement = 0; cache.worstMeasurement = 0;
MeasurementProcessor measurementProcessor( MeasurementProcessor measurementProcessor(
asymmetricCaloEnergy, *cache.fitMatrices->derivativeMatrix(), intersector, asymmetricCaloEnergy, cache.fitMatrices->derivativeMatrix(), intersector,
measurements, parameters, m_rungeKuttaIntersector, m_stepPropagator, measurements, parameters, m_rungeKuttaIntersector, m_stepPropagator,
m_useStepPropagator); m_useStepPropagator);
...@@ -593,14 +593,14 @@ const FitProcedureQuality& FitProcedure::execute( ...@@ -593,14 +593,14 @@ const FitProcedureQuality& FitProcedure::execute(
// field integral // field integral
const Amg::MatrixX* fullCovariance = cache.fitMatrices->fullCovariance(); const Amg::MatrixX* fullCovariance = cache.fitMatrices->fullCovariance();
if (fullCovariance) { if (fullCovariance) {
Amg::MatrixX* finalCovariance = cache.fitMatrices->finalCovariance(); Amg::MatrixX& finalCovariance = cache.fitMatrices->finalCovariance();
if (!for_iPatTrack) { if (!for_iPatTrack) {
if (!m_lineFit) if (!m_lineFit)
measurementProcessor.fieldIntegralUncertainty(*cache.log, measurementProcessor.fieldIntegralUncertainty(*cache.log,
*finalCovariance); finalCovariance);
measurementProcessor.propagationDerivatives(); measurementProcessor.propagationDerivatives();
} }
parameters->covariance(finalCovariance, fullCovariance); parameters->covariance(&finalCovariance, fullCovariance);
// fit quality // fit quality
if (perigeeQuality) { if (perigeeQuality) {
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment