From 2d5bf2002b8e9ca773181d48a26700a903fa3a62 Mon Sep 17 00:00:00 2001
From: James Beacham <j.beacham@cern.ch>
Date: Wed, 13 Jun 2018 20:07:18 +0000
Subject: [PATCH] Merge branch 'Remove_MA27' into '21.0'

Remove MA27 and clean up private member variables

See merge request atlas/athena!11334

(cherry picked from commit cfffd84ad3a53503b3cd420ef5ab2d05d4e74e71 [formerly 974faa2b26cb119cf525a41e2fa7432d0bdd0a8c])

1139d41a Remove MA27 and clean up private member variables
9a68c581 Remove all references to MA27
fce0467e Remove reference to MA27 from iPatFitter
51c40550 Remove final reference to MA27
449fc553 Update CMakeLists file
014fa0a3 Merge remote-tracking branch 'upstream/21.0' into Remove_MA27
a4936b8d Change solver so something more suitable for the alignment like problems
2ee38b8d code cleanup -> NULL to nullptr
df76f86b start removing cout/cerr
7b7a6af6 Remove the remaining cout/cerr calls

Former-commit-id: eca54502bfc094d3a36e60302de5de3718e77294
---
 .../CalibrationLoopScripts/Solve_trf.py       |   4 +-
 .../python/IteratorGridClasses.py             |   4 +-
 .../share/AlignmentATNSimple.py               |   2 +-
 .../share/InDetAlignAlgSetup_DBM.py           |   4 +-
 .../share/NewInDetAlignAlgSetup.py            |   4 +-
 .../share/NewInDetAlignment_DBMRel17.py       |   2 +-
 .../share/NewInDetAlignment_Run2Rel17.py      |   2 +-
 .../share/NewInDetAlignment_Run2Rel19.py      |   2 +-
 .../TrkAlgebraUtils/CMakeLists.txt            |   8 +-
 .../TrkAlgebraUtils/TrkAlgebraUtils/AlMat.h   |  43 +-
 .../TrkAlgebraUtils/AlSpaMat.h                |  16 +-
 .../TrkAlgebraUtils/AlSymMat.h                |  12 +-
 .../TrkAlgebraUtils/AlSymMatBase.h            |  46 +-
 .../TrkAlgebraUtils/TrkAlgebraUtils/AlVec.h   |  30 +-
 .../TrkAlgebraUtils/TrkAlgebraUtils/IntVec.h  |   6 +-
 .../TrkAlgebraUtils/cmt/requirements          |   2 +-
 .../TrkAlgebraUtils/src/AlMat.cxx             | 252 +++++-----
 .../TrkAlgebraUtils/src/AlSpaMat.cxx          | 476 ++++++++----------
 .../TrkAlgebraUtils/src/AlSymMat.cxx          | 248 +++++----
 .../TrkAlgebraUtils/src/AlSymMatBase.cxx      |  22 +-
 .../TrkAlgebraUtils/src/AlVec.cxx             | 240 +++++----
 .../TrkAlgebraUtils/src/IntVec.cxx            | 118 ++---
 .../TrkAlignGenTools/MatrixTool.h             |   6 +-
 .../TrkAlignGenTools/src/MatrixTool.cxx       |  24 +-
 .../TrkAlignInterfaces/IDerivCalcTool.h       |   4 +-
 25 files changed, 737 insertions(+), 840 deletions(-)

diff --git a/InnerDetector/InDetExample/InDetAlignExample/CalibrationLoopScripts/Solve_trf.py b/InnerDetector/InDetExample/InDetAlignExample/CalibrationLoopScripts/Solve_trf.py
index bc2a3f94470..970f4702752 100644
--- a/InnerDetector/InDetExample/InDetAlignExample/CalibrationLoopScripts/Solve_trf.py
+++ b/InnerDetector/InDetExample/InDetAlignExample/CalibrationLoopScripts/Solve_trf.py
@@ -194,7 +194,7 @@ solvingOption =	1
 ## solving options
 ##   0 - No global solving
 ##   1 - Lapack
-##   2 - ***REMOVED***
+##   2 - Eigen
 ##   6 - ROOT
 ##   7 - CLHEP
 
@@ -342,7 +342,7 @@ solvingOption =	1
 ## solving options
 ##   0 - No global solving
 ##   1 - Lapack
-##   2 - ***REMOVED***
+##   2 - Eigen
 ##   6 - ROOT
 ##   7 - CLHEP
 
diff --git a/InnerDetector/InDetExample/InDetAlignExample/python/IteratorGridClasses.py b/InnerDetector/InDetExample/InDetAlignExample/python/IteratorGridClasses.py
index c7e37a97b98..c4a8c776d30 100644
--- a/InnerDetector/InDetExample/InDetAlignExample/python/IteratorGridClasses.py
+++ b/InnerDetector/InDetExample/InDetAlignExample/python/IteratorGridClasses.py
@@ -418,7 +418,7 @@ class writeScriptGridForTFile :
 				("pixelAlignmentLevel"        in self.AlignmentOptions and self.AlignmentOptions["pixelAlignmentLevel"]        == 3)  or 
 				("pixelAlignmentLevelBarrel"  in self.AlignmentOptions and self.AlignmentOptions["pixelAlignmentLevelBarrel"]  == 3)  or 
 				("pixelAlignmentLevelEndcaps" in self.AlignmentOptions and self.AlignmentOptions["pixelAlignmentLevelEndcaps"] == 3) ) :
-				print "hmn, you are going to run L3 alignment, ***REMOVED*** is using, so no eigen informations!!! "
+				print "hmn, you are going to run L3 alignment, Eigen is being used, so no eigen value information!!! "
 			 
 			else : 
 
@@ -1103,7 +1103,7 @@ class writeScriptGrid :
 				("pixelAlignmentLevel"        in self.AlignmentOptions and self.AlignmentOptions["pixelAlignmentLevel"]        == 3)  or 
 				("pixelAlignmentLevelBarrel"  in self.AlignmentOptions and self.AlignmentOptions["pixelAlignmentLevelBarrel"]  == 3)  or 
 				("pixelAlignmentLevelEndcaps" in self.AlignmentOptions and self.AlignmentOptions["pixelAlignmentLevelEndcaps"] == 3) ) : 
-				print "hmn, you are going to run L3 alignment, ***REMOVED*** is using, so no eigen informations!!! "
+				print "hmn, you are going to run L3 alignment, Eigen is going to be used, so no eigen value information!!! "
 
 			else : 
 				if ( "writeEigenMat"      in self.GridOptions  and self.GridOptions["writeEigenMat"]    == True  ) and ( self.AlignmentOptions["runLocal"] == False ):
diff --git a/InnerDetector/InDetExample/InDetAlignExample/share/AlignmentATNSimple.py b/InnerDetector/InDetExample/InDetAlignExample/share/AlignmentATNSimple.py
index c3bf2fb0c60..936b0d37d20 100644
--- a/InnerDetector/InDetExample/InDetAlignExample/share/AlignmentATNSimple.py
+++ b/InnerDetector/InDetExample/InDetAlignExample/share/AlignmentATNSimple.py
@@ -256,7 +256,7 @@ runLocal = False
 ## solving option
 ##   0 - No global solving
 ##   1 - Lapack
-##   2 - ***REMOVED***
+##   2 - Eigen
 ##   6 - ROOT
 ##   7 - CLHEP
 solvingOption = 1
diff --git a/InnerDetector/InDetExample/InDetAlignExample/share/InDetAlignAlgSetup_DBM.py b/InnerDetector/InDetExample/InDetAlignExample/share/InDetAlignAlgSetup_DBM.py
index 7a9ba8861b2..74cd5bc82eb 100644
--- a/InnerDetector/InDetExample/InDetAlignExample/share/InDetAlignAlgSetup_DBM.py
+++ b/InnerDetector/InDetExample/InDetAlignExample/share/InDetAlignAlgSetup_DBM.py
@@ -35,7 +35,7 @@ newInDetAlignAlg_Options = {
 	,"inputMatrixFiles" 	        : [ "matrix.bin" ] 		     #  list of matrix files when solving only
 	,"inputVectorFiles" 	        : [ "vector.bin" ]    		     #  list of vector files when solving only
 	,"inputTFiles"			: ["AlignmentTFile.root"]	     # list of the Alignment TFiles, used only if WriteTFile is True
-	,"solvingOption" 	        : 3              		     #  which global solver to run 0-none, 1-Lapack, 2-***REMOVED***, 6-ROOT, 7-CLHEP
+	,"solvingOption" 	        : 3              		     #  which global solver to run 0-none, 1-Lapack, 2-Eigen, 6-ROOT, 7-CLHEP
 	,"solveLocal" 		        : True            		     #  whether to run locaal solving
 	,"writeMatrixFile" 	        : True      			     #  whether to write matrix to file before solving
 	,"writeMatrixFileTxt"           : True     			     #  whether to write matrix to text file
@@ -264,7 +264,7 @@ if not newInDetAlignAlg_Options["runLocal"]:
 		):
 		matrixTool.WriteEigenMat = False
 		matrixTool.UseSparse = True
-		matrixTool.SolveOption = 2 # run ***REMOVED*** for L3 by default
+		matrixTool.SolveOption = 2 # run Eigen for L3 by default
 		matrixTool.WriteMatTxt = False
 		matrixTool.WriteHitmapTxt = False
 		matrixTool.CalculateFullCovariance = False
diff --git a/InnerDetector/InDetExample/InDetAlignExample/share/NewInDetAlignAlgSetup.py b/InnerDetector/InDetExample/InDetAlignExample/share/NewInDetAlignAlgSetup.py
index f021c34f418..64ef5afcffc 100644
--- a/InnerDetector/InDetExample/InDetAlignExample/share/NewInDetAlignAlgSetup.py
+++ b/InnerDetector/InDetExample/InDetAlignExample/share/NewInDetAlignAlgSetup.py
@@ -41,7 +41,7 @@ newInDetAlignAlg_Options = {
     ,"inputMatrixFiles"             : [ "matrix.bin" ]           #  list of matrix files when solving only
     ,"inputVectorFiles"             : [ "vector.bin" ]               #  list of vector files when solving only
     ,"inputTFiles"          : ["AlignmentTFile.root"]        # list of the Alignment TFiles, used only if WriteTFile is True
-    ,"solvingOption"            : 3                          #  which global solver to run 0-none, 1-Lapack, 2-***REMOVED***, 6-ROOT, 7-CLHEP
+    ,"solvingOption"            : 3                          #  which global solver to run 0-none, 1-Lapack, 2-Eigen, 6-ROOT, 7-CLHEP
     ,"solveLocal"               : True                       #  whether to run locaal solving
     ,"writeMatrixFile"          : True                       #  whether to write matrix to file before solving
     ,"writeMatrixFileTxt"           : True                   #  whether to write matrix to text file
@@ -307,7 +307,7 @@ if not newInDetAlignAlg_Options["runLocal"]:
         ):
         matrixTool.WriteEigenMat = False
         matrixTool.UseSparse = True
-        matrixTool.SolveOption = 2 # run ***REMOVED*** for L3 by default
+        matrixTool.SolveOption = 2 # run Eigen for L3 by default
         matrixTool.WriteMatTxt = False
         matrixTool.WriteHitmapTxt = False
         matrixTool.CalculateFullCovariance = False
diff --git a/InnerDetector/InDetExample/InDetAlignExample/share/NewInDetAlignment_DBMRel17.py b/InnerDetector/InDetExample/InDetAlignExample/share/NewInDetAlignment_DBMRel17.py
index ec431f65f8f..c5e33ca483e 100644
--- a/InnerDetector/InDetExample/InDetAlignExample/share/NewInDetAlignment_DBMRel17.py
+++ b/InnerDetector/InDetExample/InDetAlignExample/share/NewInDetAlignment_DBMRel17.py
@@ -156,7 +156,7 @@ runLocal = True
 ## solving option
 ##   0 - No global solving
 ##   1 - Lapack
-##   2 - ***REMOVED***
+##   2 - Eigen
 ##   6 - ROOT
 ##   7 - CLHEP
 solvingOption = 1
diff --git a/InnerDetector/InDetExample/InDetAlignExample/share/NewInDetAlignment_Run2Rel17.py b/InnerDetector/InDetExample/InDetAlignExample/share/NewInDetAlignment_Run2Rel17.py
index 99c1b14041e..c9574e5b9dd 100644
--- a/InnerDetector/InDetExample/InDetAlignExample/share/NewInDetAlignment_Run2Rel17.py
+++ b/InnerDetector/InDetExample/InDetAlignExample/share/NewInDetAlignment_Run2Rel17.py
@@ -326,7 +326,7 @@ runLocal = False
 ## solving option
 ##   0 - No global solving
 ##   1 - Lapack
-##   2 - ***REMOVED***
+##   2 - Eigen
 ##   6 - ROOT
 ##   7 - CLHEP
 solvingOption = 1
diff --git a/InnerDetector/InDetExample/InDetAlignExample/share/NewInDetAlignment_Run2Rel19.py b/InnerDetector/InDetExample/InDetAlignExample/share/NewInDetAlignment_Run2Rel19.py
index 1d6d95783e0..b1b7c210213 100644
--- a/InnerDetector/InDetExample/InDetAlignExample/share/NewInDetAlignment_Run2Rel19.py
+++ b/InnerDetector/InDetExample/InDetAlignExample/share/NewInDetAlignment_Run2Rel19.py
@@ -328,7 +328,7 @@ runLocal = True
 ## solving option
 ##   0 - No global solving
 ##   1 - Lapack
-##   2 - ***REMOVED***
+##   2 - Eigen
 ##   6 - ROOT
 ##   7 - CLHEP
 solvingOption = 1
diff --git a/Tracking/TrkAlignment/TrkAlgebraUtils/CMakeLists.txt b/Tracking/TrkAlignment/TrkAlgebraUtils/CMakeLists.txt
index e280a0e1f10..cb2f982db3f 100644
--- a/Tracking/TrkAlignment/TrkAlgebraUtils/CMakeLists.txt
+++ b/Tracking/TrkAlignment/TrkAlgebraUtils/CMakeLists.txt
@@ -11,14 +11,16 @@ atlas_depends_on_subdirs( PRIVATE
                           TestPolicy )
 
 # External dependencies:
+find_package( BLAS )
 find_package( LAPACK )
+find_package( Eigen )
 find_package( ROOT COMPONENTS Core Matrix Tree MathCore Hist RIO pthread )
 
 # Component(s) in the package:
 atlas_add_library( TrkAlgebraUtils
                    src/*.cxx
-                   ***REMOVED***/*.f
                    PUBLIC_HEADERS TrkAlgebraUtils
-                   PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${LAPACK_INCLUDE_DIRS}
-                   PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} ${LAPACK_LIBRARIES} GaudiKernel )
+                   PRIVATE_INCLUDE_DIRS ${ROOT_INCLUDE_DIRS} ${LAPACK_INCLUDE_DIRS} ${BLAS_INCLUDE_DIRS} ${EIGEN_INCLUDE_DIRS}
+                   PRIVATE_LINK_LIBRARIES ${ROOT_LIBRARIES} ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES} ${EIGEN_LIBRARIES} GaudiKernel ${CMAKE_Fortran_IMPLICIT_LINK_LIBRARIES})
+
 
diff --git a/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlMat.h b/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlMat.h
index 1aa61aa1139..3c495e3e75f 100644
--- a/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlMat.h
+++ b/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlMat.h
@@ -7,6 +7,7 @@
 
 #include <iostream>
 #include <string>
+#include <exception>
 
 class StatusCode;
 
@@ -37,8 +38,8 @@ class AlMat {
     inline AlMat_row(AlMat&,int);
     double & operator[](int);
   private:
-    AlMat& _a;
-    int _r;
+    AlMat& m_a;
+    int m_r;
   };
 
   class AlMat_row_const {
@@ -46,8 +47,8 @@ class AlMat {
     inline AlMat_row_const(const AlMat&,int);
     const double & operator[](int) const;
   private:
-    const AlMat& _a;
-    int _r;
+    const AlMat& m_a;
+    int m_r;
   };
   // helper class to implement m[i][j]
 
@@ -104,12 +105,12 @@ class AlMat {
   void copy(const AlSymMatBase&  m);
 
  protected:
-  int _ncol, _nrow;
-  int _nele;
+  int m_ncol, m_nrow;
+  int m_nele;
   std::string m_pathbin;
   std::string m_pathtxt;
-  double* ptr_data;
-  bool transpose;
+  double* m_ptr_data;
+  bool m_transpose;
 
   virtual long int elem(long int, long int) const;
 
@@ -128,41 +129,41 @@ inline AlMat::AlMat_row_const AlMat::operator[] (int r) const{
 }
 
 inline double &AlMat::AlMat_row::operator[](int c){
-  if(_r<0||_r>=_a.nrow() || c<0||c>=_a.ncol())
-    std::cerr << "Range error in AlMat::operator[][]" << std::endl;
+  if(m_r<0||m_r>=m_a.nrow() || c<0||c>=m_a.ncol())
+    throw std::out_of_range( "Range error in AlMat::operator[][]" );
 
-  return *(_a.ptr_data+_r*_a.ncol()+c);
+  return *(m_a.m_ptr_data+m_r*m_a.ncol()+c);
 }
 
 inline const double & AlMat::AlMat_row_const::operator[](int c) const {
-  if(_r<0||_r>=_a.nrow() || c<0||c>=_a.ncol())
-    std::cerr << "Range error in AlMat::operator[][]" << std::endl;
+  if(m_r<0||m_r>=m_a.nrow() || c<0||c>=m_a.ncol())
+    throw std::out_of_range( "Range error in AlMat::operator[][]" );
 
-  return *(_a.ptr_data+_r*_a.ncol()+c);
+  return *(m_a.m_ptr_data+m_r*m_a.ncol()+c);
 }
 
 inline AlMat::AlMat_row::AlMat_row(AlMat &a,int r)
-	     : _a(a), _r(r)
+	     : m_a(a), m_r(r)
 {}
 
 inline AlMat::AlMat_row_const::AlMat_row_const
 (const AlMat&a,int r)
-	     : _a(a), _r(r)
+	     : m_a(a), m_r(r)
 {}
 
 
 inline long int AlMat::nrow() const {
-  if( transpose ) return _ncol;
-  return _nrow;
+  if( m_transpose ) return m_ncol;
+  return m_nrow;
 }
 
 inline long int AlMat::ncol() const {
-  if( transpose ) return _nrow;
-  return _ncol;
+  if( m_transpose ) return m_nrow;
+  return m_ncol;
 }
 
 inline double* AlMat::ptrData() const {
-  return ptr_data;
+  return m_ptr_data;
 }
 
 inline std::string AlMat::pathBin() const {
diff --git a/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlSpaMat.h b/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlSpaMat.h
index 4a771249ca1..6ae73cdfd77 100644
--- a/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlSpaMat.h
+++ b/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlSpaMat.h
@@ -7,7 +7,7 @@
 
 #include "TrkAlgebraUtils/AlSymMatBase.h"
 
-#include <iostream>
+#include <exception>
 
 class StatusCode;
 template<class Element> class TMatrixTSparse;
@@ -50,7 +50,7 @@ class AlSpaMat : public AlSymMatBase {
     virtual  AlSpaMat&  operator*=(const double&);
 
   // ADVANCED:
-    int ***REMOVED***Solve(AlVec& RHS);
+    int SolveWithEigen(AlVec& RHS);
     virtual  void RemoveModule( int);
     virtual  void RemoveAlignPar(int, int);
     virtual  void RemoveDoF( int, int nelem=1);
@@ -85,8 +85,8 @@ class AlSpaMat : public AlSymMatBase {
     std::string m_pathbin;
     std::string m_pathtxt;
 
-    int* ptr_row;
-    int* ptr_col;
+    int* m_ptr_row;
+    int* m_ptr_col;
 
 };
 
@@ -95,15 +95,15 @@ class AlSpaMat : public AlSymMatBase {
 inline void AlSpaMat::elem(const indices& key, long int& i, long int& j) const {
   i = key.first;
   j = key.second;
-  if(j>i) std::cerr << "AlSpaMat::elem: i<j ! " << std::endl;
+  if(j>i) throw std::out_of_range( "AlSpaMat::elem: i<j ! " );
   return;
 }
 
 inline long int AlSpaMat::nele() {
   // first check for intrinsic consistency:
-  if( int(ptr_map.size()) != _nele ) std::cerr << "AlSpaMat::_nele has been corrupted!" << std::endl;
-  _nele = ptr_map.size();
-  return _nele;
+  if( int(m_ptr_map.size()) != m_nele ) throw std::range_error( "AlSpaMat::m_nele has been corrupted!" );
+  m_nele = m_ptr_map.size();
+  return m_nele;
 }
 
 inline std::string AlSpaMat::pathBin() const {
diff --git a/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlSymMat.h b/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlSymMat.h
index f70a5c9d929..7c3e717c85e 100644
--- a/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlSymMat.h
+++ b/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlSymMat.h
@@ -83,7 +83,7 @@ class AlSymMat : public AlSymMatBase {
   virtual void copy(const AlSpaMat&  m);
   virtual void copy(const AlMat&  m);
 
-  double*     ptr_data;
+  double*     m_ptr_data;
 
   std::string m_pathbin;
   std::string m_pathtxt;
@@ -94,10 +94,10 @@ class AlSymMat : public AlSymMatBase {
 
 inline long int AlSymMat::elem(long int i,long int j) const {
 #ifdef _DEBUG
-  if( i<0 )     { std::cerr << "AlSymMat::elem: Index 1 < zero! " << i << std::endl;  return *(ptr_data); };
-  if( i>=size() ) { std::cerr << "AlSymMat::elem: Index 1 too large! " << i << std::endl;  return *(ptr_data); };
-  if( j<0 )     { std::cerr << "AlSymMat::elem: Index 2 < zero! " << j << std::endl;  return *(ptr_data);  };
-  if( j>=size() ) { std::cerr << "AlSymMat::elem: Index 2 too large! " << j << std::endl;  return *(ptr_data); };
+  if( i<0 )     { throw std::out_of_range( "AlSymMat::elem: Index 1 < zero! " ); };
+  if( i>=size() ) { throw std::out_of_range( "AlSymMat::elem: Index 1 too large! " ); };
+  if( j<0 )     { throw std::out_of_range( "AlSymMat::elem: Index 2 < zero! " );  };
+  if( j>=size() ) { throw std::out_of_range( "AlSymMat::elem: Index 2 too large! " ); };
 #endif
   if( j<=i ) {
     return  ((i+1)*i/2+j);
@@ -108,7 +108,7 @@ inline long int AlSymMat::elem(long int i,long int j) const {
 }
 
 inline double* AlSymMat::ptrData() const {
-  return ptr_data;
+  return m_ptr_data;
 }
 
 inline std::string AlSymMat::pathBin() const {
diff --git a/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlSymMatBase.h b/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlSymMatBase.h
index 93ee825e8c9..47353fe1a65 100644
--- a/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlSymMatBase.h
+++ b/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlSymMatBase.h
@@ -9,7 +9,7 @@
 // modified from SiGlobalChi2AlgebrUtils::AlSymMatBase, modified to make
 // pure virtual
 
-#include <iostream>
+#include <exception>
 #include <map>
 #include <vector>
 
@@ -46,8 +46,8 @@ class AlSymMatBase  {
     inline AlSymMatBase_row(AlSymMatBase&,long int);
     double & operator[](long int);
   private:
-    AlSymMatBase& _a;
-    long int  _r;
+    AlSymMatBase& m_a;
+    long int  m_r;
   };
 
   class AlSymMatBase_row_const {
@@ -55,8 +55,8 @@ class AlSymMatBase  {
     inline AlSymMatBase_row_const(const AlSymMatBase&,long int);
     double operator[](long int) const;
   private:
-    const AlSymMatBase& _a;
-    long int  _r;
+    const AlSymMatBase& m_a;
+    long int  m_r;
   };
   // helper class to implement m[i][j]
 
@@ -97,11 +97,11 @@ class AlSymMatBase  {
   AlSymMatBase(const AlSymMatBase&);
   AlSymMatBase& operator=(const AlSymMatBase&);
 
-  int      _matrix_type;
-  datamap  ptr_map;
+  int      m_matrix_type;
+  datamap  m_ptr_map;
 
-  long int _size;
-  long int _nele;
+  long int m_size;
+  long int m_nele;
 
 };
 
@@ -118,38 +118,36 @@ inline AlSymMatBase::AlSymMatBase_row_const AlSymMatBase::operator[] (long int r
 }
 
 inline double & AlSymMatBase::AlSymMatBase_row::operator[](long int c) {
-  if(_r<0||_r>=_a.nrow() || c<0||c>=_a.ncol()) {
-    std::cerr << "Range error in AlSymMatBase::operator[][]" << std::endl;
-    return _a.elemr(0,0);
+  if(m_r<0||m_r>=m_a.nrow() || c<0||c>=m_a.ncol()) {
+    throw std::out_of_range( "Range error in AlSymMatBase::operator[][]" );
   } else {
-    return _a.elemr(_r,c);
+    return m_a.elemr(m_r,c);
   }
 }
 
 inline double AlSymMatBase::AlSymMatBase_row_const::operator[](long int c) const {
-  if(_r<0||_r>=_a.nrow() || c<0||c>=_a.ncol()) {
-    std::cerr << "Range error in AlSymMatBase::operator[][]" << std::endl;
-    return _a.elemc(0,0);
+  if(m_r<0||m_r>=m_a.nrow() || c<0||c>=m_a.ncol()) {
+    throw std::out_of_range( "Range error in AlSymMatBase::operator[][]" );
   } else {
-    return _a.elemc(_r,c);
+    return m_a.elemc(m_r,c);
   }
 }
 
 inline AlSymMatBase::AlSymMatBase_row::AlSymMatBase_row(AlSymMatBase &a,long int r)
-  : _a(a), _r(r)
+  : m_a(a), m_r(r)
 {}
 
 inline AlSymMatBase::AlSymMatBase_row_const::AlSymMatBase_row_const
 (const AlSymMatBase&a,long int r)
-  : _a(a), _r(r)
+  : m_a(a), m_r(r)
 {}
 
 
-inline long int AlSymMatBase::size() const { return _size; }
-inline long int AlSymMatBase::nrow() const { return _size; }
-inline long int AlSymMatBase::ncol() const { return _size; }
-inline int AlSymMatBase::matrix_type() const  { return _matrix_type; }
-inline const datamap * AlSymMatBase::ptrMap() const { return &ptr_map; }
+inline long int AlSymMatBase::size() const { return m_size; }
+inline long int AlSymMatBase::nrow() const { return m_size; }
+inline long int AlSymMatBase::ncol() const { return m_size; }
+inline int AlSymMatBase::matrix_type() const  { return m_matrix_type; }
+inline const datamap * AlSymMatBase::ptrMap() const { return &m_ptr_map; }
 
 } // end namespace Trk
 
diff --git a/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlVec.h b/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlVec.h
index 18ba7c3f799..e32c8ca5c35 100644
--- a/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlVec.h
+++ b/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/AlVec.h
@@ -7,7 +7,7 @@
 
 // PBdR (17Apr2007)
 
-#include <iostream>
+#include <exception>
 #include <map>
 #include <vector>
 
@@ -66,8 +66,8 @@ class AlVec {
   inline double* ptrData() const;
 
  protected:
-  int _size;
-  double* ptr_data;
+  int m_size;
+  double* m_ptr_data;
   std::string m_pathbin;
   std::string m_pathtxt;
 
@@ -82,39 +82,35 @@ class AlVec {
 // inline methods:
 
 inline int AlVec::size() const {
-  return _size;
+  return m_size;
 }
 
 inline double* AlVec::ptrData() const {
-  return ptr_data;
+  return m_ptr_data;
 }
 
 inline double& AlVec::operator[](int i) {
   if( i < 0 ) {
-    std::cerr << "AlVec: Index < zero! " << std::endl;
-    return ptr_data[0];
+    throw std::out_of_range( "AlVec: Index < zero! " );
   }
 
-  if( i >= _size ) {
-    std::cerr << "AlVec: Index too large! " << std::endl;
-    return ptr_data[0];
+  if( i >= m_size ) {
+    throw std::out_of_range( "AlVec: Index too large! ");
   }
 
-  return *(ptr_data+i);
+  return *(m_ptr_data+i);
 }
 
 inline const double& AlVec::operator[](int i) const {
   if( i < 0 ) {
-    std::cerr << "AlVec: Index < zero! " << std::endl;
-    return ptr_data[0];
+    throw std::out_of_range( "AlVec: Index < zero! " );
   }
 
-  if( i >= _size ) {
-    std::cerr << "AlVec: Index too large! " << std::endl;
-    return ptr_data[0];
+  if( i >= m_size ) {
+    throw std::out_of_range( "AlVec: Index too large! " );
   }
 
-  return *(ptr_data+i);
+  return *(m_ptr_data+i);
 }
 
 } // end namespace Trk
diff --git a/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/IntVec.h b/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/IntVec.h
index 3d41b0cba46..ef2ebec8cae 100644
--- a/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/IntVec.h
+++ b/Tracking/TrkAlignment/TrkAlgebraUtils/TrkAlgebraUtils/IntVec.h
@@ -24,11 +24,11 @@ class IntVec {
   IntVec& operator-=(const IntVec&);
 
   void reSize(int);
-  int n_elem() const  { return Nele; }
+  int n_elem() const  { return m_Nele; }
 
  private:
-  int Nele;
-  int* ptr_data;
+  int m_Nele;
+  int* m_ptr_data;
 };
 
 } // end namespace Trk
diff --git a/Tracking/TrkAlignment/TrkAlgebraUtils/cmt/requirements b/Tracking/TrkAlignment/TrkAlgebraUtils/cmt/requirements
index c9abc1d7b0d..37ec1df5871 100644
--- a/Tracking/TrkAlignment/TrkAlgebraUtils/cmt/requirements
+++ b/Tracking/TrkAlignment/TrkAlgebraUtils/cmt/requirements
@@ -19,5 +19,5 @@ public
 
 macro_append ROOT_linkopts " -lMatrix"
 
-library TrkAlgebraUtils *.cxx -s=../***REMOVED*** *.f
+library TrkAlgebraUtils *.cxx
 apply_pattern installed_library
diff --git a/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlMat.cxx b/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlMat.cxx
index 3f0e87d02b0..27a4f310eee 100644
--- a/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlMat.cxx
+++ b/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlMat.cxx
@@ -14,49 +14,50 @@
 #include <fstream>
 #include <sstream>
 #include <stdint.h>
+#include <exception>
 
 namespace Trk {
 
 AlMat::AlMat() {
-  _ncol = 0; _nrow = 0; _nele = 0;
-  ptr_data = 0;  // set pointer to null
-  transpose = false;
+  m_ncol = 0; m_nrow = 0; m_nele = 0;
+  m_ptr_data = 0;  // set pointer to null
+  m_transpose = false;
   m_pathbin="./";
   m_pathtxt="./";
 }
 
 AlMat::AlMat(int N, int M) {
-  _nrow = N;
-  _ncol = M;
-  _nele = N*M;
-  ptr_data = new double[_nele];
-  transpose = false;
+  m_nrow = N;
+  m_ncol = M;
+  m_nele = N*M;
+  m_ptr_data = new double[m_nele];
+  m_transpose = false;
   m_pathbin="./";
   m_pathtxt="./";
 
-  double*  p = ptr_data + _nele;
-  while (p > ptr_data) *(--p) = 0.;
+  double*  p = m_ptr_data + m_nele;
+  while (p > m_ptr_data) *(--p) = 0.;
 
 }
 
 AlMat::AlMat(const AlMat& m) {
-  _nrow = m.nrow(); _ncol = m.ncol(); _nele = _nrow*_ncol;
-  ptr_data = new double[_nele];
-  transpose = false;
+  m_nrow = m.nrow(); m_ncol = m.ncol(); m_nele = m_nrow*m_ncol;
+  m_ptr_data = new double[m_nele];
+  m_transpose = false;
   m_pathbin = m.m_pathbin;
   m_pathtxt = m.m_pathtxt;
   copy(m);
 }
 
 AlMat::~AlMat()
-{if( ptr_data != NULL) delete [] ptr_data;}
+{if( m_ptr_data != nullptr ) delete [] m_ptr_data;}
 
 void AlMat::copy(const AlMat&  m) {
   int  nr = nrow();
   int  nc = ncol();
   if( nr != m.nrow() || nc != m.ncol() ) {
-    std::cerr << "AlMat::copy: sizes do not match!" << std::endl;
-    return; }
+    throw std::range_error( "AlMat::copy: sizes do not match!" );
+  }
 
   for( int i=0; i<nr; i++ ) {
     for( int j=0; j<nc; j++ ) {
@@ -68,57 +69,57 @@ void AlMat::copy(const AlMat&  m) {
 void AlMat::copy(const AlSymMatBase&  m) {
   long int  n = m.size();
   if( nrow() != n || ncol() != n ) {
-    std::cerr << "AlMat::copy: sizes do not match!" << std::endl;
-    return; }
+    throw std::range_error( "AlMat::copy: sizes do not match!" );
+  }
 
   for( int i=0; i<n; i++ ) {
     for( int j=0; j<n; j++ ) {
-      *(ptr_data+i*_ncol+j) = m.elemc(i, j);
+      *(m_ptr_data+i*m_ncol+j) = m.elemc(i, j);
     }
   }
 }
 
 double& AlMat::elemr(long int i,long int j) {
 #ifdef _DEBUG
-  if( i<0 )     { std::cerr << "AlMat::elemr: Index 1 < zero! " << i << std::endl;  return *(ptr_data); };
-  if( i>=size ) { std::cerr << "AlMat::elemr: Index 1 too large! " << i << std::endl;  return *(ptr_data); };
-  if( j<0 )     { std::cerr << "AlMat::elemr: Index 2 < zero! " << j << std::endl;  return *(ptr_data);  };
-  if( j>=size ) { std::cerr << "AlMat::elemr: Index 2 too large! " << j << std::endl;  return *(ptr_data); };
+  if( i<0 )     { throw std::out_of_range( "AlMat::elemr: Index 1 < zero! " ); };
+  if( i>=size ) { throw std::out_of_range( "AlMat::elemr: Index 1 too large!" ); };
+  if( j<0 )     { throw std::out_of_range( "AlMat::elemr: Index 2 < zero! " ); };
+  if( j>=size ) { throw std::out_of_range( "AlMat::elemr: Index 2 too large!" );  };
 #endif
 
-  if( transpose )   return  *(ptr_data+j*_ncol+i);
-  return  *(ptr_data+i*_ncol+j);
+  if( m_transpose )   return  *(m_ptr_data+j*m_ncol+i);
+  return  *(m_ptr_data+i*m_ncol+j);
 }
 
 double AlMat::elemc(long int i,long int j) const {
 #ifdef _DEBUG
-  if( i<0 )     { std::cerr << "AlMat::elemc: Index 1 < zero! " << i << std::endl;  return *(ptr_data); };
-  if( i>=size ) { std::cerr << "AlMat::elemc: Index 1 too large! " << i << std::endl;  return *(ptr_data); };
-  if( j<0 )     { std::cerr << "AlMat::elemc: Index 2 < zero! " << j << std::endl;  return *(ptr_data);  };
-  if( j>=size ) { std::cerr << "AlMat::elemc: Index 2 too large! " << j << std::endl;  return *(ptr_data); };
+  if( i<0 )     { throw std::out_of_range( "AlMat::elemr: Index 1 < zero! " ); };
+  if( i>=size ) { throw std::out_of_range( "AlMat::elemr: Index 1 too large!" ); };
+  if( j<0 )     { throw std::out_of_range( "AlMat::elemr: Index 2 < zero! " ); };
+  if( j>=size ) { throw std::out_of_range( "AlMat::elemr: Index 2 too large!" );  };
 #endif
 
-  if( transpose )   return  *(ptr_data+j*_ncol+i);
-  return  *(ptr_data+i*_ncol+j);
+  if( m_transpose )   return  *(m_ptr_data+j*m_ncol+i);
+  return  *(m_ptr_data+i*m_ncol+j);
 }
 
 long int AlMat::elem(long int i,long int j) const {
 #ifdef _DEBUG
-  if( i<0 )     { std::cerr << "AlMat::elem: Index 1 < zero! " << i << std::endl;  return *(ptr_data); };
-  if( i>=size ) { std::cerr << "AlMat::elem: Index 1 too large! " << i << std::endl;  return *(ptr_data); };
-  if( j<0 )     { std::cerr << "AlMat::elem: Index 2 < zero! " << j << std::endl;  return *(ptr_data);  };
-  if( j>=size ) { std::cerr << "AlMat::elem: Index 2 too large! " << j << std::endl;  return *(ptr_data); };
+  if( i<0 )     { throw std::out_of_range( "AlMat::elemr: Index 1 < zero! " ); };
+  if( i>=size ) { throw std::out_of_range( "AlMat::elemr: Index 1 too large!" ); };
+  if( j<0 )     { throw std::out_of_range( "AlMat::elemr: Index 2 < zero! " ); };
+  if( j>=size ) { throw std::out_of_range( "AlMat::elemr: Index 2 too large!" );  };
 #endif
 
-  if( transpose )   return  (j*_ncol+i);
-  return  (i*_ncol+j);
+  if( m_transpose )   return  (j*m_ncol+i);
+  return  (i*m_ncol+j);
 }
 
 
 AlMat&  AlMat::operator=(const double& d) {
 
-  double*  p = ptr_data + _nele;
-  while (p > ptr_data) *(--p) = d;
+  double*  p = m_ptr_data + m_nele;
+  while (p > m_ptr_data) *(--p) = d;
 
   return *this;
 }
@@ -126,10 +127,12 @@ AlMat&  AlMat::operator=(const double& d) {
 AlMat& AlMat::operator=(const AlMat& m) {
   if (this==&m) return *this;
 
-  if(( _nrow!=0 || _ncol!=0) && (_nrow != m.nrow() || _ncol != m.ncol() ))
-    { std::cerr << "AlMat=AlMat Assignment: size do not match!" << std::endl; return *this; }
+  if(( m_nrow!=0 || m_ncol!=0) && (m_nrow != m.nrow() || m_ncol != m.ncol() ))
+  { 
+    throw std::range_error( "AlMat=AlMat Assignment: size do not match!" );
+  }
 
-  if ( ptr_data != m.ptr_data ) {
+  if ( m_ptr_data != m.m_ptr_data ) {
     reSize(m.nrow(), m.ncol());
     copy(m);
   };
@@ -137,10 +140,12 @@ AlMat& AlMat::operator=(const AlMat& m) {
 }
 
 AlMat&  AlMat::operator=(const AlSymMat& m) {
-  if( ( _nrow!=0 || _ncol!=0) && (_nrow != m.size() || _ncol != m.size() ))
-    { std::cerr << "AlMat=AlSymMatBase Assignment: size do not match!" << std::endl; return *this; }
+  if( ( m_nrow!=0 || m_ncol!=0) && (m_nrow != m.size() || m_ncol != m.size() ))
+  {
+    throw std::range_error( "AlMat=AlSymMatBase Assignment: size do not match!" ); 
+  }
 
-  if ( ptr_data != m.ptrData() ) {
+  if ( m_ptr_data != m.ptrData() ) {
     reSize(m.size(), m.size());
     copy(m);
   }
@@ -148,59 +153,69 @@ AlMat&  AlMat::operator=(const AlSymMat& m) {
 }
 
 AlMat AlMat::operator+(const AlMat& m) const {
-  if( _nrow != m._nrow || _ncol != m._ncol )
-    { std::cerr << "AlMat: operator+: size do not match!" << std::endl; return *this; }
+  if( m_nrow != m.m_nrow || m_ncol != m.m_ncol )
+  { 
+    throw std::range_error( "AlMat: operator+: size do not match!" );
+  }
 
-  AlMat b( _nrow, _ncol );
+  AlMat b( m_nrow, m_ncol );
 
-  double*  p = ptr_data + _nele;
-  double*  q = m.ptr_data + _nele;
-  double*  r = b.ptr_data + _nele;
-  while (p > ptr_data) *(--r) = (*(--p))+(*(--q));
+  double*  p = m_ptr_data + m_nele;
+  double*  q = m.m_ptr_data + m_nele;
+  double*  r = b.m_ptr_data + m_nele;
+  while (p > m_ptr_data) *(--r) = (*(--p))+(*(--q));
 
   return b;
 }
 
 AlMat&  AlMat::operator+=(const AlMat& m) {
 
-  if( _nrow != m._nrow || _ncol != m._ncol )
-    { std::cerr << "AlMat: operator+=: size do not match!" << std::endl; return *this; }
+  if( m_nrow != m.m_nrow || m_ncol != m.m_ncol )
+  {  
+    throw std::range_error( "AlMat: operator+=: size do not match!" ); 
+  }
 
-  double*  p = ptr_data + _nele;
-  double*  q = m.ptr_data + _nele;
-  while (p > ptr_data) *(--p) += *(--q);
+  double*  p = m_ptr_data + m_nele;
+  double*  q = m.m_ptr_data + m_nele;
+  while (p > m_ptr_data) *(--p) += *(--q);
 
   return *this;
 }
 
 AlMat   AlMat::operator-(const AlMat& m) const {
-  if( _nrow != m._nrow || _ncol != m._ncol )
-    { std::cerr << "AlMat: operator-: size do not match!" << std::endl; return *this; }
+  if( m_nrow != m.m_nrow || m_ncol != m.m_ncol )
+  {  
+    throw std::range_error( "AlMat: operator-: size do not match!" ); 
+  }
 
-  AlMat b( _nrow, _ncol );
+  AlMat b( m_nrow, m_ncol );
 
-  double*  p = ptr_data + _nele;
-  double*  q = m.ptr_data + _nele;
-  double*  r = b.ptr_data + _nele;
-  while (p > ptr_data) *(--r) = (*(--p))-(*(--q));
+  double*  p = m_ptr_data + m_nele;
+  double*  q = m.m_ptr_data + m_nele;
+  double*  r = b.m_ptr_data + m_nele;
+  while (p > m_ptr_data) *(--r) = (*(--p))-(*(--q));
 
   return b;
 }
 
 AlMat&  AlMat::operator-=(const AlMat& m) {
-  if( _nrow != m._nrow || _ncol != m._ncol )
-    { std::cerr << "AlMat: operator-=: size do not match!" << std::endl; return *this; }
+  if( m_nrow != m.m_nrow || m_ncol != m.m_ncol )
+  { 
+    throw std::range_error( "AlMat: operator-=: size do not match!"); 
+  }
 
-  double*  p = ptr_data + _nele;
-  double*  q = m.ptr_data + _nele;
-  while (p > ptr_data) *(--p) -= *(--q);
+  double*  p = m_ptr_data + m_nele;
+  double*  q = m.m_ptr_data + m_nele;
+  while (p > m_ptr_data) *(--p) -= *(--q);
 
   return *this;
 }
 
 AlMat AlMat::operator*(const AlMat& m) const {
-  if( ncol() != m.nrow() ) {  std::cerr << "AlMat: operator*: size do not match!" << std::endl;
-  return m; }
+  if( ncol() != m.nrow() ) 
+  {  
+    throw std::range_error( "AlMat: operator*: size do not match!" );
+  }
 
   int k(nrow());
   int l(m.ncol());
@@ -216,8 +231,10 @@ AlMat AlMat::operator*(const AlMat& m) const {
 }
 
 AlMat AlMat::operator*(const AlSymMatBase& m) const {
-  if( ncol() != m.size()) {  std::cerr << "AlMat: operator*: size do not match!" << std::endl;
-  return *this; }
+  if( ncol() != m.size()) 
+  { 
+    throw std::range_error( "AlMat: operator*: size do not match!" ); 
+  }
 
   int k(nrow());
   int l(m.size());
@@ -234,7 +251,9 @@ AlMat AlMat::operator*(const AlSymMatBase& m) const {
 
 AlVec AlMat::operator*(const AlVec& v) const {
   if( ncol() != v.size() )
-    { std::cerr << "AlMat: operator*: size do not match! " << std::endl; return v; }
+  {  
+    throw std::range_error( "AlMat: operator*: size do not match! " );
+  }
 
   int k(nrow());
   int l(ncol());
@@ -244,7 +263,7 @@ AlVec AlMat::operator*(const AlVec& v) const {
   AlVec b(k);
   for( int i=0; i<k; i++ ) {
     p = b.ptrData()+i;
-    q = ptr_data+i*l;
+    q = m_ptr_data+i*l;
     for( int j=0; j<l; j++ )  *p += (*(q+j))*(*(v.ptrData()+j));
   }
 
@@ -252,8 +271,8 @@ AlVec AlMat::operator*(const AlVec& v) const {
 }
 
 AlMat&   AlMat::operator*=(const double& d) {
-  double*  p = ptr_data+_nele;
-  while (p > ptr_data)
+  double*  p = m_ptr_data+m_nele;
+  while (p > m_ptr_data)
     *(--p) *= d;
 
   return *this;
@@ -261,10 +280,10 @@ AlMat&   AlMat::operator*=(const double& d) {
 
 // transposition
 AlMat AlMat::T() const {
-  AlMat b(_ncol,_nrow);
-  for( int i=0; i<b._nrow; i++ ) {
-    for( int j=0; j<b._ncol; j++ ) {
-      b[i][j]= *(ptr_data+j*_ncol+i);
+  AlMat b(m_ncol,m_nrow);
+  for( int i=0; i<b.m_nrow; i++ ) {
+    for( int j=0; j<b.m_ncol; j++ ) {
+      b[i][j]= *(m_ptr_data+j*m_ncol+i);
     }
   }
   return b;
@@ -273,7 +292,7 @@ AlMat AlMat::T() const {
 // transposition
 AlMat& AlMat::Transpose() {
 
-  transpose = true;
+  m_transpose = true;
 
   return *this;
 }
@@ -281,7 +300,7 @@ AlMat& AlMat::Transpose() {
 // normal position
 AlMat& AlMat::Normal() {
 
-  transpose = false;
+  m_transpose = false;
 
   return *this;
 }
@@ -289,12 +308,11 @@ AlMat& AlMat::Normal() {
 
 //invert sym matrix declared as non-symetric for convenience
 void AlMat::invertS(int& ierr, double Norm = 1.) {
-  if(_nrow!=_ncol) {
-    std::cerr << "AlMat invertS: non-square matrix!" << std::endl;
-    return;
+  if(m_nrow!=m_ncol) {
+    throw std::range_error( "AlMat invertS: non-square matrix!" );
   }
 
-  AlSymMat b(_nrow);
+  AlSymMat b(m_nrow);
   for( int i=0; i<b.size(); i++ ) {
     for( int j=0; j<=i; j++ ) {
       b.elemr(i, j) = elemc(i, j);
@@ -311,23 +329,23 @@ void AlMat::invertS(int& ierr, double Norm = 1.) {
 
 // reSize
 void AlMat::reSize(int Nnew, int Mnew) {
-  if ( Nnew != _nrow || Mnew != _ncol ) {
+  if ( Nnew != m_nrow || Mnew != m_ncol ) {
 
-    double*  p = ptr_data;
-    _nele = Nnew*Mnew;
-    ptr_data = new double[_nele];
-    int _nrow_old = _nrow;
-    int _ncol_old = _ncol;
-    transpose = false;
+    double*  p = m_ptr_data;
+    m_nele = Nnew*Mnew;
+    m_ptr_data = new double[m_nele];
+    int m_nrow_old = m_nrow;
+    int m_ncol_old = m_ncol;
+    m_transpose = false;
 
-    _nrow = Nnew;
-    _ncol = Mnew;
-    int k = _nrow <= _nrow_old ? _nrow : _nrow_old;
-    int l = _ncol <= _ncol_old ? _ncol : _ncol_old;
+    m_nrow = Nnew;
+    m_ncol = Mnew;
+    int k = m_nrow <= m_nrow_old ? m_nrow : m_nrow_old;
+    int l = m_ncol <= m_ncol_old ? m_ncol : m_ncol_old;
 
     for( int i=0; i<k; i++ ) {
       for( int j=0; j<l; j++ ) {
-        *(ptr_data+i*_ncol+j) = *(p+i*_ncol_old+j);
+        *(m_ptr_data+i*m_ncol+j) = *(p+i*m_ncol_old+j);
       }
     }
 
@@ -341,25 +359,25 @@ StatusCode AlMat::ReadScalaPack(const std::string &filename){
   if(inmat.fail())
     return StatusCode::FAILURE;
 
-  int32_t mNrow = _nrow;
+  int32_t mNrow = m_nrow;
   inmat.read((char*)&mNrow, sizeof (mNrow));
-  int32_t mNcol = _ncol;
+  int32_t mNcol = m_ncol;
   inmat.read((char*)&mNcol, sizeof (mNcol));
 
-  _nrow=abs(mNrow);
-  _ncol=abs(mNcol);
-  _nele=_ncol*_nrow;
-  transpose = false;
+  m_nrow=abs(mNrow);
+  m_ncol=abs(mNcol);
+  m_nele=m_ncol*m_nrow;
+  m_transpose = false;
 
-  // printf("ALMat::nrow: %d \n",_nrow);
-  // printf("ALMat::ncol: %d \n",_ncol);
-  // printf("ALMat::nele: %d \n",_nele);
+  // printf("ALMat::nrow: %d \n",m_nrow);
+  // printf("ALMat::ncol: %d \n",m_ncol);
+  // printf("ALMat::nele: %d \n",m_nele);
 
   double melem=0;
-  for(int i=0; i<_nrow; i++) {
-    for(int j=0; j<_ncol; j++) {
+  for(int i=0; i<m_nrow; i++) {
+    for(int j=0; j<m_ncol; j++) {
       inmat.read((char*)&melem, sizeof (melem));
-      *(ptr_data+i*_ncol+j) = melem;
+      *(m_ptr_data+i*m_ncol+j) = melem;
       // printf("(%d,%d) = %.16lf \n",i,j,melem);
     }
   }
@@ -390,9 +408,9 @@ StatusCode AlMat::Write(const std::string &filename, bool binary, unsigned int p
     if(outmat.fail())
       return StatusCode::FAILURE;
 
-    int32_t mNrow = _nrow;
+    int32_t mNrow = m_nrow;
     outmat.write((char*)&mNrow, sizeof (mNrow));
-    int32_t mNcol = _ncol;
+    int32_t mNcol = m_ncol;
     outmat.write((char*)&mNcol, sizeof (mNcol));
   }
   else {
@@ -407,9 +425,9 @@ StatusCode AlMat::Write(const std::string &filename, bool binary, unsigned int p
 
   double melem=0;
 
-  for( int i=0; i<_nrow; i++) {
-    for( int j=0; j<_ncol; j++) {
-      melem =  *(ptr_data+i*_ncol+j);
+  for( int i=0; i<m_nrow; i++) {
+    for( int j=0; j<m_ncol; j++) {
+      melem =  *(m_ptr_data+i*m_ncol+j);
       if(binary)
         outmat.write((char*)&(melem), sizeof (melem));
       else
@@ -455,7 +473,7 @@ const std::string AlMat::Print(const int NColsPerSet)
     for (int row=0; row<nrow(); row++) {
       textmatrix << " |" << std::setw(6) << row << " | ";
       for (int col=FirstCol; col<LastCol; col++) {
-    double melem =  *(ptr_data+row*ncol()+col);
+    double melem =  *(m_ptr_data+row*ncol()+col);
     textmatrix << std::setprecision(5) << std:: setw(6) << melem << " | ";
       }
       textmatrix << std::endl;
diff --git a/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlSpaMat.cxx b/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlSpaMat.cxx
index 84b7bde47b1..528a0692668 100644
--- a/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlSpaMat.cxx
+++ b/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlSpaMat.cxx
@@ -19,88 +19,76 @@
 #include <math.h>
 #include <float.h>
 #include <stdint.h>
+#include <stdexcept>
 
 #include <TMatrixDSparse.h>
 
-extern"C" {
-  void ma27id_(int ICNTL[],double CNTL[]);
 
-  void ma27ad_(int *N,int *NZ,int IRN[],int ICN[],
-               int IW[],int *LIW,int IKEEP[],int IW1p[],
-               int *NSTEPS,int *IFLAG,int ICNTL[],double CNTL[],
-               int INFO[],double *OPS);
+#include <Eigen/Core>
+#include <Eigen/IterativeLinearSolvers>
+#include <Eigen/SparseCholesky>
 
-  void ma27bd_(int *N,int *NZ,int IRN[],int ICN[],
-               double A[],int *LA,int IW[],int *LIW,int IKEEP[],
-               int *NSTEPS,int *MAXFRT, int IW1[],int ICNTL[],
-               double CNTL[], int INFO[]);
-
-  void ma27cd_(int *N, double A[],int *LA,int IW[],int *LIW,
-               double W[],int *MAXFRT,double RHS[],int IW1[],int *NSTEPS,
-               int ICNTL[], int INFO[]);
-}
 
 namespace Trk {
 
 //______________________________________________________________________________
 AlSpaMat::AlSpaMat()
 {
-  _matrix_type = 2;
-  _size = 0;
-  _nele = 0;
-  ptr_row = NULL;
-  ptr_col = NULL;    // set pointer to null
+  m_matrix_type = 2;
+  m_size = 0;
+  m_nele = 0;
+  m_ptr_row = nullptr;
+  m_ptr_col = nullptr;    // set pointer to null
 }
 
 //______________________________________________________________________________
 AlSpaMat::AlSpaMat(long int N)
 {
-  _matrix_type = 2;
-  _size = N;
-  _nele = 0;
-  ptr_row = NULL;
-  ptr_col = NULL;    // set pointer to null
+  m_matrix_type = 2;
+  m_size = N;
+  m_nele = 0;
+  m_ptr_row = nullptr;
+  m_ptr_col = nullptr;    // set pointer to null
 }
 
 //______________________________________________________________________________
 AlSpaMat::AlSpaMat(const AlSpaMat& m)
  : AlSymMatBase(m)
 {
-  _matrix_type = 2;
-  _size = m.size();
-  ptr_row = NULL;
-  ptr_col = NULL;    // set pointer to null
+  m_matrix_type = 2;
+  m_size = m.size();
+  m_ptr_row = nullptr;
+  m_ptr_col = nullptr;    // set pointer to null
   copy(m);
 }
 
 //______________________________________________________________________________
 AlSpaMat::AlSpaMat(const AlSymMat& m)
 {
-  _matrix_type = 2;
-  _size = m.size();
-  ptr_row = NULL;
-  ptr_col = NULL;    // set pointer to null
+  m_matrix_type = 2;
+  m_size = m.size();
+  m_ptr_row = nullptr;
+  m_ptr_col = nullptr;    // set pointer to null
   copy(m);
 }
 
 //______________________________________________________________________________
 AlSpaMat::~AlSpaMat()
 {
-  ptr_map.clear();
-  if( ptr_row != NULL )  delete [] ptr_row;
-  if( ptr_col != NULL )  delete [] ptr_col;
+  m_ptr_map.clear();
+  if( m_ptr_row != nullptr )  delete [] m_ptr_row;
+  if( m_ptr_col != nullptr )  delete [] m_ptr_col;
 }
 
 //______________________________________________________________________________
 void AlSpaMat::copy(const AlSpaMat& m)
 {
   if( size() != m.size()) {
-    std::cerr << "AlSpaMat::copy: size do not match!" << std::endl;
-    return;
+    throw std::range_error( "AlSpaMat::copy: size do not match!" );
   }
-  ptr_map.clear();
-  _nele=m._nele;
-  ptr_map = m.ptr_map;
+  m_ptr_map.clear();
+  m_nele=m.m_nele;
+  m_ptr_map = m.m_ptr_map;
   return;
 }
 
@@ -108,18 +96,17 @@ void AlSpaMat::copy(const AlSpaMat& m)
 void AlSpaMat::copy(const AlSymMat& m)
 {
   if( size() != m.size()) {
-    std::cerr << "AlSpaMat::copy: size do not match!" << std::endl;
-    return;
+    throw std::range_error( "AlSpaMat::copy: size do not match!" );
   }
-  ptr_map.clear();
-  _nele=0;
+  m_ptr_map.clear();
+  m_nele=0;
 
   // Convert Storage System
-  for (int i = 0; i<_size; i++)
+  for (int i = 0; i<m_size; i++)
     for (int j = 0; j<=i; j++)
       if  (m.elemc(i,j) != 0.) {
-        ptr_map.insert(std::make_pair(std::make_pair(i,j), m.elemc(i,j)));
-        _nele++;
+        m_ptr_map.insert(std::make_pair(std::make_pair(i,j), m.elemc(i,j)));
+        m_nele++;
       }
 
   return;
@@ -129,20 +116,19 @@ void AlSpaMat::copy(const AlSymMat& m)
 void AlSpaMat::copy(const AlMat& m)
 {
   if( size() != m.nrow() || size() != m.ncol() ) {
-    std::cerr << "AlSpaMat::copy: size do not match!" << std::endl;
-    return;
+    throw std::range_error( "AlSpaMat::copy: size do not match!" );
   }
 
   // copy just the lower triangle:
-  ptr_map.clear();
-  _nele=0;
+  m_ptr_map.clear();
+  m_nele=0;
 
   //Convert Storage System
-  for (int i = 0; i<_size; i++)
+  for (int i = 0; i<m_size; i++)
     for (int j = 0; j<=i; j++)
       if  (m.elemc(i,j) != 0.) {
-        ptr_map.insert(std::make_pair(std::make_pair(i,j), m.elemc(i,j)));
-        _nele++;
+        m_ptr_map.insert(std::make_pair(std::make_pair(i,j), m.elemc(i,j)));
+        m_nele++;
       }
 
   return;
@@ -153,38 +139,34 @@ double& AlSpaMat::elemr(long int i,long int j)
 {
 #ifdef _DEBUG
   if( i<0 ) {
-    std::cerr << "AlSpaMat::elemr: Index 1 < zero! " << i << std::endl;
-    return ptr_map.begin()->second;
+    throw std::out_of_range( "AlSpaMat::elemr: Index 1 < zero! " );
   }
   if( i>=size() ) {
-    std::cerr << "AlSpaMat::elemr: Index 1 too large! " << i << std::endl;
-    return ptr_map.begin()->second;
+    throw std::out_of_range( "AlSpaMat::elemr: Index 1 too large! " );
   }
   if( j<0 ) {
-    std::cerr << "AlSpaMat::elemr: Index 2 < zero! " << j << std::endl;
-    return ptr_map.begin()->second;
+    throw std::out_of_range( "AlSpaMat::elemr: Index 2 < zero! " );
   }
   if( j>=size() ) {
-    std::cerr << "AlSpaMat::elemr: Index 2 too large! " << j << std::endl;
-    return ptr_map.begin()->second;
+    throw std::out_of_range("AlSpaMat::elemr: Index 2 too large! " );;
   }
 #endif
   // try fast referencing:
   mapiterator pos;
   indices key = j < i ? std::make_pair(i,j) : std::make_pair(j,i);
-  pos=ptr_map.find(key);
-  if( pos!=ptr_map.end() )    // already exists
+  pos=m_ptr_map.find(key);
+  if( pos!=m_ptr_map.end() )    // already exists
     return pos->second;
   else {                      // does not yet exist
     // create it with the null value:
-    _nele++;
-    return (ptr_map.insert(std::make_pair(key, 0.))).first->second;
+    m_nele++;
+    return (m_ptr_map.insert(std::make_pair(key, 0.))).first->second;
   }
   // simpler implementation of the above:
   /*
-    ptr_map[key]=0;
-    _nele = ptr_map.size();
-    return ptr_map[key];
+    m_ptr_map[key]=0;
+    m_nele = m_ptr_map.size();
+    return m_ptr_map[key];
   */
 }
 
@@ -193,28 +175,24 @@ double  AlSpaMat::elemc(long int i,long int j) const
 {
 #ifdef _DEBUG
   if( i<0 ) {
-    std::cerr << "AlSpaMat::elemc: Index 1 < zero! " << i << std::endl;
-    return 0.0;
+    throw std::out_of_range( "AlSpaMat::elemc: Index 1 < zero! " );
   }
   if( i>=size() ) {
-    std::cerr << "AlSpaMat::elemc: Index 1 too large! " << i << std::endl;
-    return 0.0;
+    throw std::out_of_range( "AlSpaMat::elemc: Index 1 too large! " );
   }
   if( j<0 ) {
-    std::cerr << "AlSpaMat::elemc: Index 2 < zero! " << j << std::endl;
-    return 0.0;
+    throw std::out_of_range( "AlSpaMat::elemc: Index 2 < zero! " );
   }
   if( j>=size() ) {
-    std::cerr << "AlSpaMat::elemc: Index 2 too large! " << j << std::endl;
-    return 0.0;
+    throw std::out_of_range( "AlSpaMat::elemc: Index 2 too large! " );
   }
 #endif
   // try fast referencing:
   const_mapiterator pos;
   indices key = j < i ? std::make_pair(i,j) : std::make_pair(j,i);
-  pos=ptr_map.find(key);
+  pos=m_ptr_map.find(key);
 
-  if( pos!=ptr_map.end() )    // already exists
+  if( pos!=m_ptr_map.end() )    // already exists
     return pos->second;
 
   // does not yet exist
@@ -228,20 +206,16 @@ indices  AlSpaMat::elem(long int i,long int j) const
   // ATTENTION! the key value is returned:
 #ifdef _DEBUG
   if( i<0 ) {
-    std::cerr << "AlSpaMat::elem: Index 1 < zero! " << i << std::endl;
-    return std::make_pair(0,0);
+    throw std::out_of_range( "AlSpaMat::elem: Index 1 < zero! " );
   }
   if( i>=size() ) {
-    std::cerr << "AlSpaMat::elem: Index 1 too large! " << i << std::endl;
-    return std::make_pair(0,0);
+    throw std::out_of_range( "AlSpaMat::elem: Index 1 too large! " );
   }
   if( j<0 ) {
-    std::cerr << "AlSpaMat::elem: Index 2 < zero! " << j << std::endl;
-    return std::make_pair(0,0);
+    throw std::out_of_range( "AlSpaMat::elem: Index 2 < zero! " );
   }
   if( j>=size() ) {
-    std::cerr << "AlSpaMat::elem: Index 2 too large! " << j << std::endl;
-    return std::make_pair(0,0);
+    throw std::out_of_range( "AlSpaMat::elem: Index 2 too large! " );
   }
 #endif
 
@@ -254,7 +228,7 @@ AlSpaMat&  AlSpaMat::operator=(const AlSpaMat& m)
 {
   
  if ( this != &m ) {
-    _size=m.size();
+    m_size=m.size();
     copy(m);
   }
   return *this;
@@ -263,7 +237,7 @@ AlSpaMat&  AlSpaMat::operator=(const AlSpaMat& m)
 //______________________________________________________________________________
 AlSpaMat&  AlSpaMat::operator=(const AlSymMat& m)
 {
-  _size=m.size();
+  m_size=m.size();
   copy(m);
   return *this;
 }
@@ -272,11 +246,10 @@ AlSpaMat&  AlSpaMat::operator=(const AlSymMat& m)
 AlSpaMat&  AlSpaMat::operator=(const AlMat& m)
 {
   if( m.nrow() != m.ncol() ) {
-    std::cerr << "AlSpaMat::=operator: allowed for square matrices only!" << std::endl;
-    return *this;
+    throw std::range_error( "AlSpaMat::=operator: allowed for square matrices only!" );
   }
 
-  _size=m.nrow();
+  m_size=m.nrow();
   copy(m);
   return *this;
 }
@@ -285,7 +258,7 @@ AlSpaMat&  AlSpaMat::operator=(const AlMat& m)
 AlSpaMat&  AlSpaMat::operator=(const double& d)
 {
   mapiterator pos;
-  for (pos = ptr_map.begin(); pos!=ptr_map.end(); pos++)
+  for (pos = m_ptr_map.begin(); pos!=m_ptr_map.end(); pos++)
     pos->second = d;
 
   return *this;
@@ -295,15 +268,14 @@ AlSpaMat&  AlSpaMat::operator=(const double& d)
 AlSpaMat AlSpaMat::operator+(const AlSpaMat& m) const
 {
   if( size() != m.size()) {
-    std::cerr << "AlSpaMat::operator+: size do not match!" << std::endl;
-    return *this;
+    throw std::range_error(  "AlSpaMat::operator+: size do not match!" );
   }
 
   AlSpaMat b(m);
   const_mapiterator pos;
-  for (pos = ptr_map.begin(); pos!=ptr_map.end(); pos++)
-    b.ptr_map[pos->first] += pos->second;
-  b._nele = b.ptr_map.size();
+  for (pos = m_ptr_map.begin(); pos!=m_ptr_map.end(); pos++)
+    b.m_ptr_map[pos->first] += pos->second;
+  b.m_nele = b.m_ptr_map.size();
 
   return b;
 }
@@ -312,14 +284,13 @@ AlSpaMat AlSpaMat::operator+(const AlSpaMat& m) const
 AlSpaMat&  AlSpaMat::operator+=(const AlSpaMat& m)
 {
   if( size() != m.size()) {
-    std::cerr << "AlSpaMat::operator+=: size do not match!" << std::endl;
-    return *this;
+    throw std::range_error( "AlSpaMat::operator+=: size do not match!" );
   }
 
   const_mapiterator pos;
-  for (pos = m.ptr_map.begin(); pos!=m.ptr_map.end(); pos++)
-    (*this).ptr_map[pos->first] += pos->second;
-  _nele = ptr_map.size();
+  for (pos = m.m_ptr_map.begin(); pos!=m.m_ptr_map.end(); pos++)
+    (*this).m_ptr_map[pos->first] += pos->second;
+  m_nele = m_ptr_map.size();
 
   return *this;
 }
@@ -328,15 +299,14 @@ AlSpaMat&  AlSpaMat::operator+=(const AlSpaMat& m)
 AlSpaMat AlSpaMat::operator-(const AlSpaMat& m) const
 {
   if( size() != m.size()) {
-    std::cerr << "AlSpaMat::operator-: size do not match!" << std::endl;
-    return *this;
+    throw std::range_error(  "AlSpaMat::operator-: size do not match!" );
   }
 
   AlSpaMat b(m);
   const_mapiterator pos;
-  for (pos = ptr_map.begin(); pos!=ptr_map.end(); pos++)
-    b.ptr_map[pos->first] -= pos->second;
-  b._nele = b.ptr_map.size();
+  for (pos = m_ptr_map.begin(); pos!=m_ptr_map.end(); pos++)
+    b.m_ptr_map[pos->first] -= pos->second;
+  b.m_nele = b.m_ptr_map.size();
 
   return b;
 }
@@ -345,14 +315,13 @@ AlSpaMat AlSpaMat::operator-(const AlSpaMat& m) const
 AlSpaMat&  AlSpaMat::operator-=(const AlSpaMat& m)
 {
   if( size() != m.size()) {
-    std::cerr << "AlSpaMat::operator-=: size do not match!" << std::endl;
-    return *this;
+    throw std::range_error(  "AlSpaMat::operator-=: size do not match!" );
   }
 
   const_mapiterator pos;
-  for (pos = m.ptr_map.begin(); pos!=m.ptr_map.end(); pos++)
-    (*this).ptr_map[pos->first] -= pos->second;
-  _nele = ptr_map.size();
+  for (pos = m.m_ptr_map.begin(); pos!=m.m_ptr_map.end(); pos++)
+    (*this).m_ptr_map[pos->first] -= pos->second;
+  m_nele = m_ptr_map.size();
 
   return *this;
 }
@@ -361,9 +330,7 @@ AlSpaMat&  AlSpaMat::operator-=(const AlSpaMat& m)
 AlMat AlSpaMat::operator*(const AlSymMatBase& m) const
 {
   if( size() != m.size() ) {
-    std::cerr << "AlSpaMat::operator*: size do not match!" << std::endl;
-    AlMat b( size(), m.size());
-    return b;
+    throw std::range_error(  "AlSpaMat::operator*: size do not match!" );
   }
 
   long int  isiz(size());
@@ -380,8 +347,7 @@ AlMat AlSpaMat::operator*(const AlSymMatBase& m) const
 //______________________________________________________________________________
 AlMat  AlSpaMat::operator*(const AlMat& m) const {
   if( size() != m.nrow() ) {
-    std::cerr << "AlSpaMat::operator*: size do not match!" << std::endl;
-    return m;
+    throw std::range_error(  "AlSpaMat::operator*: size do not match!" );
   }
 
   long int  isiz(size());
@@ -400,8 +366,7 @@ AlMat  AlSpaMat::operator*(const AlMat& m) const {
 AlVec AlSpaMat::operator*(const AlVec& v) const
 {
   if( size() != v.size() ) {
-    std::cerr << "AlSpaMat::operator*: size do not match! " << std::endl;
-    return v;
+    throw std::range_error(  "AlSpaMat::operator*: size do not match! " );
   }
 
   long int  isiz(size());
@@ -418,7 +383,7 @@ AlVec AlSpaMat::operator*(const AlVec& v) const
 AlSpaMat&  AlSpaMat::operator*=(const double& d)
 {
   mapiterator pos;
-  for (pos = ptr_map.begin(); pos!=ptr_map.end(); pos++)
+  for (pos = m_ptr_map.begin(); pos!=m_ptr_map.end(); pos++)
     pos->second *= d;
 
   return *this;
@@ -429,8 +394,8 @@ AlSpaMat  AlSpaMat::operator*(const double& d) const
 {
   AlSpaMat a(size());
   const_mapiterator pos;
-  for (pos = ptr_map.begin(); pos!=ptr_map.end(); pos++)
-    a.ptr_map.insert(std::make_pair(pos->first, (pos->second)*d));
+  for (pos = m_ptr_map.begin(); pos!=m_ptr_map.end(); pos++)
+    a.m_ptr_map.insert(std::make_pair(pos->first, (pos->second)*d));
 
   return a;
 }
@@ -438,118 +403,94 @@ AlSpaMat  AlSpaMat::operator*(const double& d) const
 //______________________________________________________________________________
 double AlSpaMat::determinant()
 {
-  double deter = 1.;
-
-  std::cerr << "AlSpaMat::determinant: not implemented!" << std::endl;
-
-  return deter;
+  throw std::invalid_argument( "AlSpaMat::determinant: not implemented!" );
 }
 
-//______________________________________________________________________________
-int AlSpaMat::***REMOVED***Solve(AlVec& RHS)
-{
-  int ICNTL[30];
-  double CNTL[5];
-  ma27id_(ICNTL,CNTL);
-
-  double* ptr_data = new double[_nele];
-  ptr_row = new int[_nele];
-  ptr_col = new int[_nele];
-
-  int Size(_size);
-  int Nele(_nele);
-
-  //Convert Storage System
+int AlSpaMat::SolveWithEigen(AlVec& RHS){
+ 
+  if(RHS.size() != size() ){
+    throw std::range_error(  "AlSpaMat::SolveWithEigen vector size is incorrect" );
+  }
+  
+  Eigen::VectorXd eigenBigVector( RHS.size() );
+  for(int i(0); i< RHS.size(); ++i ){
+    eigenBigVector[i] = RHS[i];
+  }
 
+  Eigen::SparseMatrix<double> eigenBigMatrix( m_size, m_size );
+  
+  typedef Eigen::Triplet<double> Triplet;
+  std::vector<Triplet> tripletList;
+  tripletList.reserve(m_nele);
   long int      i, j;
   long int counter(0);
   mapiterator pos;
-  for (pos = ptr_map.begin(); pos!=ptr_map.end(); pos++){
+  for (pos = m_ptr_map.begin(); pos!=m_ptr_map.end(); pos++){
     elem(pos->first, i, j);
-    *(ptr_data+counter)=pos->second;
-    *(ptr_row+counter)= i+1;
-    *(ptr_col+counter)= j+1;
+    tripletList.push_back(Triplet(i,j,pos->second));
+    if(i!=j) tripletList.push_back(Triplet(j,i,pos->second));
     counter++;
   }
+  eigenBigMatrix.setFromTriplets(tripletList.begin(), tripletList.end());
 
-  int LIW =3/2*(2*Nele+3*Size+1);
-  int* IW = new int[LIW];
-  int* IKEEP = new int[3*Size];
-  int* IW1 = new int[2*Size];
-  int NSTEPS;
-  int IFLAG = 0;
-  int INFO[20];
-  double OPS;
-
-  ma27ad_(&Size, &Nele,ptr_row,ptr_col,
-          IW,&LIW,IKEEP,IW1,
-          &NSTEPS,&IFLAG,ICNTL,CNTL,
-          INFO,&OPS);
-
-  int MAXFRT;
-  int LA =2*INFO[4];
+  // Eigen::CholmodSupernodalLLT  is much much quicker (x50) so it would be great to move to that in the future
+  // requires an external package SuiteSparse
+  // BiCGSTAB was the fastest iterative solver in Eigen from a quick test
+  // SimplicialLDLT was the fastest direct solver in Eigen (x2 slower than BiCGSTAB ) 
 
-  double* TempA = new double[Nele];
-  for (int i=0; i<Nele ;i++)
-    *(TempA+i) = *(ptr_data+i);
+  // Eigen::BiCGSTAB<Eigen::SparseMatrix<double> > solver;
+  Eigen::SimplicialLDLT<Eigen::SparseMatrix<double> > solver;
 
-  delete [] ptr_data;
-  ptr_data = new double[LA];
-
-  for (int i=0; i<LA ;i++)
-    *(ptr_data+i)=0;
-
-  for (int i=0; i<Nele ;i++)
-    *(ptr_data+i)= *(TempA+i);
-
-  delete [] TempA;
-
-  ma27bd_(&Size,&Nele,ptr_row,ptr_col,
-          ptr_data,&LA,IW,&LIW,IKEEP,
-          &NSTEPS,&MAXFRT,IW1,ICNTL,
-          CNTL,INFO);
+  solver.compute(eigenBigMatrix);
+  if(solver.info()!=Eigen::Success) {
+    // decomposition failed
+    throw std::domain_error("AlSpaMat::SolveWithEigen: failed to compute: -- is the input matrix singular?" );
+    return 1;
+  }
 
-  double* W =new double[MAXFRT];
+  Eigen::VectorXd x = solver.solve( eigenBigVector );
+  if(solver.info()!=Eigen::Success) {
+  // solving failed
+    throw std::domain_error("AlSpaMat::SolveWithEigen: Failed to solve: -- is your matrix singular? ");
+    return 2;
+  }
+  
+  //Copy results into vector
+  for(int i(0); i< RHS.size(); ++i ){
+    RHS[i] = x[i];
+  }
 
-  double* RHStemp =new double[Size];
-  for (int i=0; i<Size; i++)
-    RHStemp[i] = RHS[i];
+  //Check accuracy
+  Eigen::VectorXd residual = eigenBigMatrix * x - eigenBigVector;
+  double sumresidual = 0;
+  for( int i=0; i<residual.size(); ++i){
+    sumresidual += fabs(residual[i]);
+  }
+  sumresidual /= (double) residual.size();
+  
+  if( sumresidual > 1e-3 ){
+    throw std::overflow_error( "AlSpaMat::SolveWithEigen: your solution is no good! ");
+    return 3;
+  }
 
-  ma27cd_(&Size, ptr_data,&LA,IW,&LIW,
-          W,&MAXFRT,RHStemp,IW1,&NSTEPS,
-          ICNTL,INFO);
+  return 0;
 
-  for (int i=0; i<Size; i++)
-    RHS[i] = RHStemp[i];
 
-  delete [] RHStemp;
-  delete [] IW;
-  delete [] IKEEP;
-  delete [] IW1;
-  delete [] W;
+}
 
-  delete [] ptr_data;  ptr_data=NULL;
-  delete [] ptr_row;   ptr_row=NULL;
-  delete [] ptr_col;   ptr_col=NULL;
 
-  return INFO[0];
-}
 
 //______________________________________________________________________________
 //jlove int AlSpaMat::diagonalize(char jobz, AlVec& w, AlMat& z) {
 int AlSpaMat::diagonalize(char, AlVec&, AlMat&)
 {
-  std::cerr << "AlSpaMat::diagonalize: not implemented!" << std::endl;
-  int ierr = -1;
-  return ierr;
+  throw std::invalid_argument( "AlSpaMat::diagonalize: not implemented!" );
 }
 
 //______________________________________________________________________________
 int AlSpaMat::invert()
 {
-  std::cerr << "AlSpaMat::invert: not implemented!" << std::endl;
-  int ierr = -1;
-  return ierr;
+  throw std::invalid_argument( "AlSpaMat::invert: not implemented!" );
 }
 
 //______________________________________________________________________________
@@ -557,16 +498,16 @@ void AlSpaMat::RemoveDoF(int index, int nelem)
 {
   // index here is the first alignment parameter to remove, nelem number of consecutive ones
 
-  _size-=nelem;
+  m_size-=nelem;
   long int n=index+nelem-1;
   mapiterator pos;
-  mapiterator pos_obs=ptr_map.end();
-  for (pos = ptr_map.begin(); pos!=ptr_map.end(); pos++) {
+  mapiterator pos_obs=m_ptr_map.end();
+  for (pos = m_ptr_map.begin(); pos!=m_ptr_map.end(); pos++) {
 
-    if( pos_obs!=ptr_map.end() )
-      ptr_map.erase(pos_obs);
+    if( pos_obs!=m_ptr_map.end() )
+      m_ptr_map.erase(pos_obs);
 
-    pos_obs = ptr_map.end();
+    pos_obs = m_ptr_map.end();
 
     long int k, l;
     elem(pos->first, k, l);
@@ -585,9 +526,9 @@ void AlSpaMat::RemoveDoF(int index, int nelem)
     }
   }
 
-  if( pos_obs!=ptr_map.end() )
-    ptr_map.erase(pos_obs);
-  _nele = ptr_map.size();
+  if( pos_obs!=m_ptr_map.end() )
+    m_ptr_map.erase(pos_obs);
+  m_nele = m_ptr_map.size();
 }
 
 //______________________________________________________________________________
@@ -595,12 +536,11 @@ int AlSpaMat::RemoveCollsRows(std::vector<int> indices)
 {
   int n = indices.size();
   if (n==0) {
-    std::cerr<<"Vector of indices to remove is empty."<<std::endl;
-    return _size;
+    return m_size;
   }
-  if (n>_size) {
-    std::cerr<<"Vector of indices larger than matrix size."<<std::endl;
-    return _size;
+  if (n>m_size) {
+    throw std::invalid_argument( "AlSpaMat::RemoveCollsRows: Vector of indices larger than matrix size." );
+    return m_size;
   }
 
   // first sort the list of indices descending
@@ -615,14 +555,13 @@ int AlSpaMat::RemoveCollsRows(std::vector<int> indices)
 
   // remove rows and columns starting from largest indices
   for (int i=0;i<n;i++) {
-    if (indices[i] > _size-1) {
-      std::cerr<<"Index "<<indices[i]<<" goes beyond matrix (size "<<_size<<")."<<std::endl;
-      continue;
+    if (indices[i] > m_size-1) {
+      throw std::out_of_range( "AlSpaMat::RemoveCollsRows: Index goes beyond matrix." );
     }
     RemoveDoF(indices[i]);
   }
 
-  return _size;
+  return m_size;
 }
 
 //______________________________________________________________________________
@@ -652,15 +591,15 @@ void AlSpaMat::RemoveAlignPar(int index, int control)
         counterCol=0;
       }
 
-      pos = ptr_map.find(std::make_pair(row+shiftRow,col+shiftCol));
-      if( pos!=ptr_map.end() ) {    // exists
+      pos = m_ptr_map.find(std::make_pair(row+shiftRow,col+shiftCol));
+      if( pos!=m_ptr_map.end() ) {    // exists
         elemr(row,col) = pos->second;
-        ptr_map.erase(pos);
+        m_ptr_map.erase(pos);
       }
       else {
-        pos = ptr_map.find(std::make_pair(row,col));
-        if( pos!=ptr_map.end() )
-          ptr_map.erase(pos);
+        pos = m_ptr_map.find(std::make_pair(row,col));
+        if( pos!=m_ptr_map.end() )
+          m_ptr_map.erase(pos);
       }
       counterCol++;
       if (counterCol==5-control) {
@@ -676,7 +615,7 @@ void AlSpaMat::RemoveAlignPar(int index, int control)
     }
   }
 
-  _nele = ptr_map.size();
+  m_nele = m_ptr_map.size();
 }
 
 //______________________________________________________________________________
@@ -694,21 +633,21 @@ void AlSpaMat::RemoveModule(int index)
 void AlSpaMat::reSize(long int n)
 {
   // balanced binary tree!
-  if( n!=_size ) {
+  if( n!=m_size ) {
     // make a temporary copy:
     AlSpaMat m(*this);
-    ptr_map.clear();
-    _size = n;
+    m_ptr_map.clear();
+    m_size = n;
     long int i, j;
     mapiterator pos;
-    for (pos = m.ptr_map.begin(); pos!=m.ptr_map.end(); pos++) {
+    for (pos = m.m_ptr_map.begin(); pos!=m.m_ptr_map.end(); pos++) {
       m.elem(pos->first, i, j);
       if( i<n && j<n )
-        ptr_map.insert(*pos);
+        m_ptr_map.insert(*pos);
     }
   }
 
-  _nele = ptr_map.size();
+  m_nele = m_ptr_map.size();
   return;
 }
 
@@ -728,12 +667,12 @@ void AlSpaMat::SetPathTxt(const std::string &path)
 //jlove StatusCode AlSpaMat::Write(const std::string &filename, bool binary,
 //                           bool square, double scale, float version){
 StatusCode AlSpaMat::Write(const std::string &filename, bool binary,
-                           bool square, double, float version)
+                           bool /*square*/, double, float version)
 {
   std::ofstream outmat;
 
   int32_t msizz = 1000000+size();
-  int32_t nelem = ptr_map.size();
+  int32_t nelem = m_ptr_map.size();
 
   if(binary) {
     outmat.open((m_pathbin+filename).c_str(), std::ios::binary);
@@ -744,9 +683,6 @@ StatusCode AlSpaMat::Write(const std::string &filename, bool binary,
     outmat.write((char*)&msizz, sizeof (msizz));
     outmat.write((char*)&version, sizeof (version));
     outmat.write((char*)&nelem, sizeof (nelem));
-    // std::cout << "AlSpaMat::Write: msizz = " << msizz << std::endl;
-    // std::cout << "AlSpaMat::Write: version = " << version << std::endl;
-    // std::cout << "AlSpaMat::Write: nelem = " << nelem << std::endl;
   }
   else {
     outmat.open((m_pathtxt+filename).c_str());
@@ -765,7 +701,7 @@ StatusCode AlSpaMat::Write(const std::string &filename, bool binary,
   int32_t i, j;
 
   mapiterator pos;
-  for (pos = ptr_map.begin(); pos!=ptr_map.end(); pos++) {
+  for (pos = m_ptr_map.begin(); pos!=m_ptr_map.end(); pos++) {
     melem = pos->second;
     elem(pos->first, ii, jj);   i=ii;  j=jj;     // just a type conversion
     if(binary) {
@@ -777,8 +713,6 @@ StatusCode AlSpaMat::Write(const std::string &filename, bool binary,
       outmat << std::setw(6) << i  << std::setw(6) << j << std::setw(18) << melem << std::endl;
   }
   outmat.close();
-  if(square)
-    std::cout << "AlSpaMat::Write: square flag has been ignored!" << std::endl;
 
   return StatusCode::SUCCESS;
 }
@@ -790,7 +724,7 @@ StatusCode AlSpaMat::CheckMatVersion(const std::string filename, bool &StdUnits)
   if(inmat.fail())
     return StatusCode::FAILURE;
 
-  ptr_map.clear();
+  m_ptr_map.clear();
 
   int32_t msiz=0;
   inmat.read((char*)&msiz, sizeof (msiz));
@@ -805,8 +739,6 @@ StatusCode AlSpaMat::CheckMatVersion(const std::string filename, bool &StdUnits)
 
   inmat.close();
 
-  // std::cout << "AlSpaMat::StdUnits: " << StdUnits << std::endl;
-
   return StatusCode::SUCCESS;
 }
 
@@ -817,13 +749,11 @@ StatusCode AlSpaMat::Read(const std::string &filename, int &dofs, bool &triang,
   if (StatusCode::SUCCESS != CheckMatVersion(m_pathbin+filename, m_StdUnits))
     return StatusCode::FAILURE;
 
-  // std::cout << "AlSpaMat::StdUnits: " << m_StdUnits << std::endl;
-
   std::ifstream inmat((m_pathbin+filename).c_str(), std::ios::binary);
   if(inmat.fail())
     return StatusCode::FAILURE;
 
-  ptr_map.clear();
+  m_ptr_map.clear();
 
   int32_t msiz=0;
   int32_t nelem;
@@ -832,11 +762,10 @@ StatusCode AlSpaMat::Read(const std::string &filename, int &dofs, bool &triang,
     dofs = msiz-1000000;
   else
     dofs = abs(msiz);
-  _size = dofs;
+  m_size = dofs;
 
   if (m_StdUnits)
     inmat.read((char*)&version, sizeof (version));
-  // std::cout << "AlSpaMat::Write: version = " << version << std::endl;
 
   double melem=0;
   int32_t i, j;
@@ -844,13 +773,12 @@ StatusCode AlSpaMat::Read(const std::string &filename, int &dofs, bool &triang,
   if(msiz>999999) { // sparse format
     triang=false;
     inmat.read((char*)&nelem, sizeof (nelem));
-    _nele=nelem;
-    // std::cout << "AlSpaMat::Write: nelem = " << nelem << std::endl;
+    m_nele=nelem;
     for(int k=0; k<nelem; k++) {
       inmat.read((char*)&i, sizeof (i));
       inmat.read((char*)&j, sizeof (j));
       inmat.read((char*)&melem, sizeof (melem));
-      ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
+      m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
     }
   }
   else if(msiz>0) { // square format
@@ -858,26 +786,24 @@ StatusCode AlSpaMat::Read(const std::string &filename, int &dofs, bool &triang,
     for(int32_t i=0; i<msiz; i++) {
       for(int32_t j=0; j<msiz; j++) {
         inmat.read((char*)&melem, sizeof (melem));
-        // std::cout << "AlSpaMat::Write: nelem = " << nelem << std::endl;
         if( i>=j && melem!=0. )
-          ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
+          m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
       }
     }
   }
   else { // triangular format
     triang=true;
     msiz = (-1)*msiz;
-    // std::cout << "msiz="  << msiz << std::endl;
     for( int32_t i=0; i<msiz; i++) {
       for( int32_t j=0; j<=i; j++) {
         inmat.read((char*)&melem, sizeof (melem));
         if( melem!=0. )
-          ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
+          m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
       }
     }
   }
 
-  _nele = ptr_map.size();
+  m_nele = m_ptr_map.size();
   inmat.close();
   return StatusCode::SUCCESS;
 }
@@ -890,7 +816,7 @@ StatusCode AlSpaMat::ReadProjected(const std::string &filename, int &dofs,
   if(inmat.fail())
     return StatusCode::FAILURE;
 
-  ptr_map.clear();
+  m_ptr_map.clear();
 
   int32_t msiz=0;
   int32_t nelem;
@@ -899,7 +825,7 @@ StatusCode AlSpaMat::ReadProjected(const std::string &filename, int &dofs,
     dofs = msiz-1000000;
   else
     dofs = abs(msiz);
-  _size = dofs;
+  m_size = dofs;
 
   inmat.read((char*)&version, sizeof (version));
 
@@ -909,13 +835,12 @@ StatusCode AlSpaMat::ReadProjected(const std::string &filename, int &dofs,
   if(msiz>999999) { // sparse format
     triang=false;
     inmat.read((char*)&nelem, sizeof (nelem));
-    _nele=nelem;
-//    std::cout << "AlSpaMat::Write: nelem = " << nelem << std::endl;
+    m_nele=nelem;
     for(int k=0; k<nelem; k++) {
       inmat.read((char*)&i, sizeof (i));
       inmat.read((char*)&j, sizeof (j));
       inmat.read((char*)&melem, sizeof (melem));
-      ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
+      m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
     }
   }
   else if(msiz>0) { // square format
@@ -924,24 +849,23 @@ StatusCode AlSpaMat::ReadProjected(const std::string &filename, int &dofs,
       for(int32_t j=0; j<msiz; j++) {
         inmat.read((char*)&melem, sizeof (melem));
         if( i>=j && melem!=0. )
-          ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
+          m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
       }
     }
   }
   else { // triangular format
     triang=true;
     msiz = (-1)*msiz;
-//    std::cout << "msiz="  << msiz << std::endl;
     for( int32_t i=0; i<msiz; i++) {
       for( int32_t j=0; j<=i; j++) {
         inmat.read((char*)&melem, sizeof (melem));
         if( melem!=0. )
-          ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
+          m_ptr_map.insert(std::make_pair(std::make_pair(i,j), melem));
       }
     }
   }
 
-  _nele = ptr_map.size();
+  m_nele = m_ptr_map.size();
   inmat.close();
   return StatusCode::SUCCESS;
 }
@@ -949,15 +873,15 @@ StatusCode AlSpaMat::ReadProjected(const std::string &filename, int &dofs,
 TMatrixDSparse* AlSpaMat::makeTMatrix()
 {
 
-  long int nonZeroElements = ptr_map.size();
+  long int nonZeroElements = m_ptr_map.size();
   int    *irow = new int[nonZeroElements];
   int    *icol = new int[nonZeroElements];
   double *val  = new double[nonZeroElements];
 
   long int      i, j;
   long int counter(0);
-  const_mapiterator pos = ptr_map.begin();
-  for(pos=ptr_map.begin(); pos!=ptr_map.end(); pos++){
+  const_mapiterator pos = m_ptr_map.begin();
+  for(pos=m_ptr_map.begin(); pos!=m_ptr_map.end(); pos++){
     i = pos->first.first;
     j = pos->first.second;
     *(val+counter)=pos->second;
@@ -966,9 +890,7 @@ TMatrixDSparse* AlSpaMat::makeTMatrix()
      counter++;
   }
 
-  std::cout << counter << " " << nonZeroElements << " size " << _size  <<std::endl;
-
-  TMatrixDSparse* myTMatrix = new TMatrixDSparse(0,_size-1,0,_size-1);
+  TMatrixDSparse* myTMatrix = new TMatrixDSparse(0,m_size-1,0,m_size-1);
   myTMatrix->SetMatrixArray(nonZeroElements,irow,icol,val);
 
   delete [] irow;
diff --git a/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlSymMat.cxx b/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlSymMat.cxx
index c31c41e2ea5..b63210ce40c 100644
--- a/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlSymMat.cxx
+++ b/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlSymMat.cxx
@@ -16,6 +16,7 @@
 #include <math.h>
 #include <float.h> //for DBL_EPSILON
 #include <stdint.h>
+#include <stdexcept>
 
 #include <TMatrixDSparse.h>
 
@@ -38,10 +39,10 @@ namespace Trk {
 //______________________________________________________________________________
 AlSymMat::AlSymMat()
 {
-  _matrix_type = 1;
-  _size = 0;
-  _nele = 0;
-  ptr_data = NULL;  // set pointer to null
+  m_matrix_type = 1;
+  m_size = 0;
+  m_nele = 0;
+  m_ptr_data = nullptr;  // set pointer to null
   m_pathbin="./";
   m_pathtxt="./";
 }
@@ -50,15 +51,15 @@ AlSymMat::AlSymMat()
 //______________________________________________________________________________
 AlSymMat::AlSymMat(long int N)
 {
-  _matrix_type = 1;
-  _size = N;
-  _nele = N*(N+1)/2;
-  ptr_data = new double[_nele];
+  m_matrix_type = 1;
+  m_size = N;
+  m_nele = N*(N+1)/2;
+  m_ptr_data = new double[m_nele];
   m_pathbin="./";
   m_pathtxt="./";
 
-  double*  p = ptr_data + _nele;
-  while (p > ptr_data) *(--p) = 0.;
+  double*  p = m_ptr_data + m_nele;
+  while (p > m_ptr_data) *(--p) = 0.;
 
 }
 
@@ -66,12 +67,12 @@ AlSymMat::AlSymMat(long int N)
 AlSymMat::AlSymMat(const AlSymMat& m)
   : AlSymMatBase(m)
 {
-  _matrix_type = 1;
-  _size      = m.size();
-  _nele      = m._nele;
+  m_matrix_type = 1;
+  m_size      = m.size();
+  m_nele      = m.m_nele;
   m_pathbin = m.m_pathbin;
   m_pathtxt = m.m_pathtxt;
-  ptr_data = new double[_nele];
+  m_ptr_data = new double[m_nele];
   copy(m);
 }
 
@@ -79,25 +80,25 @@ AlSymMat::AlSymMat(const AlSymMat& m)
 //______________________________________________________________________________
 AlSymMat::AlSymMat(const AlSpaMat& m)
 {
-  _matrix_type = 1;
-  _size      = m.size();
-  _nele      = _size*(_size+1)/2;
+  m_matrix_type = 1;
+  m_size      = m.size();
+  m_nele      = m_size*(m_size+1)/2;
   m_pathbin = m.pathBin();
   m_pathtxt = m.pathTxt();
-  ptr_data = new double[_nele];
+  m_ptr_data = new double[m_nele];
   copy(m);
 }
 
 //______________________________________________________________________________
 AlSymMat& AlSymMat::operator=(const AlSpaMat& m)
 {
-  _matrix_type = 1;
-  _size = m.size();
-  _nele = _size*(_size+1)/2;
+  m_matrix_type = 1;
+  m_size = m.size();
+  m_nele = m_size*(m_size+1)/2;
   m_pathbin = m.pathBin();
   m_pathtxt = m.pathTxt();
 
-  ptr_data = new double[_nele];
+  m_ptr_data = new double[m_nele];
   copy(m);
   return *this;
 }
@@ -105,7 +106,7 @@ AlSymMat& AlSymMat::operator=(const AlSpaMat& m)
 //______________________________________________________________________________
 AlSymMat::~AlSymMat()
 {
-  if( ptr_data != NULL ) delete [] ptr_data;
+  if( m_ptr_data != nullptr ) delete [] m_ptr_data;
   //ptr_map.clear();
 }
 
@@ -114,13 +115,12 @@ void AlSymMat::copy(const AlSymMat& m)
 {
   // this one implements the fast copy alg.
   if( size() != m.size()) {
-    std::cerr << "AlSymMat::copy: size do not match!" << std::endl;
-    return;
+    throw std::range_error( "AlSymMat::copy: size do not match!" );
   }
 
-  double * p = ptr_data + _nele;
-  double*  r = m.ptr_data + _nele;
-  while (p > ptr_data) *(--p) = *(--r);
+  double * p = m_ptr_data + m_nele;
+  double*  r = m.m_ptr_data + m_nele;
+  while (p > m_ptr_data) *(--p) = *(--r);
 
   return;
 }
@@ -130,8 +130,7 @@ void AlSymMat::copy(const AlSpaMat& m)
 {
   //
   if( size() != m.size()) {
-    std::cerr << "AlSymMat::copy: size do not match!" << std::endl;
-    return;
+    throw std::range_error( "AlSymMat::copy: size do not match!" );
   }
 
   (*this) = 0.;
@@ -152,8 +151,7 @@ void AlSymMat::copy(const AlMat& m)
   // copy the lower triangle only!
   int  si = size();
   if( si != m.nrow() || si != m.ncol() ) {
-    std::cerr << "AlSymMat::copy: sizes do not match!" << std::endl;
-    return;
+    throw std::range_error("AlSymMat::copy: sizes do not match!" );
   }
 
   for( int i=0; i<si; i++ )
@@ -167,7 +165,7 @@ AlSymMat& AlSymMat::operator=(const AlSymMat& m)
   if (this==&m)
     return *this;
 
-  if ( ptr_data != m.ptr_data ) {
+  if ( m_ptr_data != m.m_ptr_data ) {
     reSize(m.size());
     copy(m);
   }
@@ -178,8 +176,7 @@ AlSymMat& AlSymMat::operator=(const AlSymMat& m)
 AlSymMat& AlSymMat::operator=(const AlMat& m)
 {
   if( m.nrow() != m.ncol() ) {
-    std::cerr << "AlSymMat::operator=: a squared matrix is required!" << std::endl;
-    return *this;
+    throw std::range_error( "AlSymMat::operator=: a squared matrix is required!" );
   }
 
   reSize(m.nrow());
@@ -191,8 +188,8 @@ AlSymMat& AlSymMat::operator=(const AlMat& m)
 //______________________________________________________________________________
 AlSymMat& AlSymMat::operator=(const double& d)
 {
-  double * p = ptr_data + _nele;
-  while (p > ptr_data) *(--p) = d;
+  double * p = m_ptr_data + m_nele;
+  while (p > m_ptr_data) *(--p) = d;
 
   return *this;
 }
@@ -201,15 +198,14 @@ AlSymMat& AlSymMat::operator=(const double& d)
 AlSymMat AlSymMat::operator+(const AlSymMat& m) const
 {
   if( size() != m.size()) {
-    std::cerr << "AlSymMat::operator+: size do not match!" << std::endl;
-    return *this;
+    throw std::range_error( "AlSymMat::operator+: size do not match!" );
   }
 
   AlSymMat b(size());
-  double * p = ptr_data + _nele;
-  double * q = m.ptr_data + _nele;
-  double * r = b.ptr_data + _nele;
-  while (p > ptr_data) *(--r) = (*(--p))+(*(--q));
+  double * p = m_ptr_data + m_nele;
+  double * q = m.m_ptr_data + m_nele;
+  double * r = b.m_ptr_data + m_nele;
+  while (p > m_ptr_data) *(--r) = (*(--p))+(*(--q));
 
   return b;
 }
@@ -218,13 +214,12 @@ AlSymMat AlSymMat::operator+(const AlSymMat& m) const
 AlSymMat&  AlSymMat::operator+=(const AlSymMat& m)
 {
   if( size() != m.size()){
-    std::cerr << "AlSymMat::operator+=: size do not match!" << std::endl;
-    return *this;
+    throw std::range_error( "AlSymMat::operator+=: size do not match!" );
   }
 
-  double * p = ptr_data + _nele;
-  double * q = m.ptr_data + _nele;
-  while (p > ptr_data) *(--p) += *(--q);
+  double * p = m_ptr_data + m_nele;
+  double * q = m.m_ptr_data + m_nele;
+  while (p > m_ptr_data) *(--p) += *(--q);
 
   return *this;
 }
@@ -233,15 +228,14 @@ AlSymMat&  AlSymMat::operator+=(const AlSymMat& m)
 AlSymMat AlSymMat::operator-(const AlSymMat& m) const
 {
   if( size() != m.size()) {
-    std::cerr << "AlSymMat::operator-: size do not match!" << std::endl;
-    return *this;
+    throw std::range_error( "AlSymMat::operator-: size do not match!" );
   }
 
   AlSymMat b(size());
-  double * p = ptr_data + _nele;
-  double * q = m.ptr_data + _nele;
-  double * r = b.ptr_data + _nele;
-  while (p > ptr_data) *(--r) = (*(--p))-(*(--q));
+  double * p = m_ptr_data + m_nele;
+  double * q = m.m_ptr_data + m_nele;
+  double * r = b.m_ptr_data + m_nele;
+  while (p > m_ptr_data) *(--r) = (*(--p))-(*(--q));
 
   return b;
 }
@@ -250,13 +244,12 @@ AlSymMat AlSymMat::operator-(const AlSymMat& m) const
 AlSymMat&  AlSymMat::operator-=(const AlSymMat& m)
 {
   if( size() != m.size()) {
-    std::cerr << "AlSymMat::operator-=: size do not match!" << std::endl;
-    return *this;
+    throw std::range_error( "AlSymMat::operator-=: size do not match!" );
   }
 
-  double * p = ptr_data + _nele;
-  double * q = m.ptr_data + _nele;
-  while (p > ptr_data) *(--p) -= *(--q);
+  double * p = m_ptr_data + m_nele;
+  double * q = m.m_ptr_data + m_nele;
+  while (p > m_ptr_data) *(--p) -= *(--q);
 
   return *this;
 }
@@ -265,9 +258,7 @@ AlSymMat&  AlSymMat::operator-=(const AlSymMat& m)
 AlMat AlSymMat::operator*(const AlSymMat& m) const
 {
   if( size() != m.size() ) {
-    std::cerr << "AlSymMat::operator*: size do not match!" << std::endl;
-    AlMat b( size(), m.size());
-    return b;
+    throw std::range_error( "AlSymMat::operator*: size do not match!" );
   }
 
   AlMat b( size(), size());
@@ -284,17 +275,17 @@ AlMat AlSymMat::operator*(const AlSymMat& m) const
       ii = (i+1)*i/2;
       jj = (j+1)*j/2;
       for( k=0; k<l; k++)
-        x += (*(ptr_data+ii+k))*(*(m.ptr_data+jj+k));
+        x += (*(m_ptr_data+ii+k))*(*(m.m_ptr_data+jj+k));
       if( i<j ) {
         for( k=l; k<n; k++)
-          x += (*(ptr_data+(k+1)*k/2+i))*(*(m.ptr_data+jj+k));
+          x += (*(m_ptr_data+(k+1)*k/2+i))*(*(m.m_ptr_data+jj+k));
       }
       else {
         for( k=l; k<n; k++)
-          x += (*(ptr_data+ii+k))*(*(m.ptr_data+(k+1)*k/2+j));
+          x += (*(m_ptr_data+ii+k))*(*(m.m_ptr_data+(k+1)*k/2+j));
       }
       for( k=n; k<siz; k++)
-        x += (*(ptr_data+(k+1)*k/2+i))*(*(m.ptr_data+(k+1)*k/2+j));
+        x += (*(m_ptr_data+(k+1)*k/2+i))*(*(m.m_ptr_data+(k+1)*k/2+j));
       b.elemr(i,j) = x;
     }
 
@@ -305,8 +296,7 @@ AlMat AlSymMat::operator*(const AlSymMat& m) const
 AlMat AlSymMat::operator*(const AlMat& m) const
 {
   if( size() != m.nrow() ) {
-    std::cerr << "AlSymMat::operator*: size do not match!" << std::endl;
-    return m;
+    throw std::range_error("AlSymMat::operator*: size do not match!" );
   }
 
   AlMat b( size(), m.ncol());
@@ -321,9 +311,9 @@ AlMat AlSymMat::operator*(const AlMat& m) const
       x = 0.;
       ii = (i+1)*i/2;
       for( k=0; k<i; k++)
-        x += (*(ptr_data+ii+k))*m.elemc(k,j);
+        x += (*(m_ptr_data+ii+k))*m.elemc(k,j);
       for( k=i; k<siz; k++)
-        x += (*(ptr_data+(k+1)*k/2+i))*m.elemc(k,j);
+        x += (*(m_ptr_data+(k+1)*k/2+i))*m.elemc(k,j);
       b.elemr(i,j) = x;
     }
 
@@ -334,8 +324,7 @@ AlMat AlSymMat::operator*(const AlMat& m) const
 AlVec AlSymMat::operator*(const AlVec& v) const
 {
   if( size() != v.size() ) {
-    std::cerr << "AlSymMat::operator*: size do not match! " << std::endl;
-    return v;
+    throw std::range_error( "AlSymMat::operator*: size do not match! " );
   }
 
   AlVec b(size());
@@ -345,9 +334,9 @@ AlVec AlSymMat::operator*(const AlVec& v) const
   for( long int i=0; i<size(); i++ ) {
     ii = (i+1)*i/2;
     for( long int j=0; j<i; j++ )
-      b[i] += (*(ptr_data+ii+j))*v[j];
+      b[i] += (*(m_ptr_data+ii+j))*v[j];
     for( long int j=i; j<size(); j++ )
-      b[i] += (*(ptr_data+(j+1)*j/2+i))*v[j];
+      b[i] += (*(m_ptr_data+(j+1)*j/2+i))*v[j];
   }
 
   return b;
@@ -356,8 +345,8 @@ AlVec AlSymMat::operator*(const AlVec& v) const
 //______________________________________________________________________________
 AlSymMat&  AlSymMat::operator*=(const double& d)
 {
-  double * p = ptr_data + _nele;
-  while (p > ptr_data) *(--p) *= d;
+  double * p = m_ptr_data + m_nele;
+  while (p > m_ptr_data) *(--p) *= d;
 
   return *this;
 }
@@ -366,9 +355,9 @@ AlSymMat&  AlSymMat::operator*=(const double& d)
 AlSymMat  AlSymMat::operator*(const double& d) const
 {
   AlSymMat a(size());
-  double * p = ptr_data + _nele;
-  double * q = a.ptr_data + _nele;
-  while (p > ptr_data) *(--q) = (*(--p))*d;
+  double * p = m_ptr_data + m_nele;
+  double * q = a.m_ptr_data + m_nele;
+  while (p > m_ptr_data) *(--q) = (*(--p))*d;
 
   return a;
 }
@@ -387,7 +376,7 @@ int AlSymMat::invert()
   int   N = size();
   int    * ipiv =  new int[N];
   double * work =  new double[N];
-  double * mx = ptr_data;
+  double * mx = m_ptr_data;
 
   dsptrf_(&uplo, &N, mx, ipiv, &ierr);
 
@@ -475,7 +464,7 @@ int AlSymMat::diagonalize(char jobz, AlVec &w, AlMat &z)
   int   N = size();
   double * work =  new double[3*N];
 
-  double * ap = ptr_data;
+  double * ap = m_ptr_data;
 
   dspev_(&jobz, &uplo, &N, ap, w.ptrData(), z.ptrData(), &N, work, &info);
 
@@ -494,12 +483,12 @@ void AlSymMat::RemoveModule(int index)
     for(int row=(shift*index); row<(size()-shift); row+=shift) {
       for(int k=0; k<shift; k++) {
         for(int col=0; col<row+k+1; col++) {
-          //  cout << "(" << row+k << "," << col << ") -> " << *(ptr_data+elem(row+k,col)) << endl;
+          //  cout << "(" << row+k << "," << col << ") -> " << *(m_ptr_data+elem(row+k,col)) << endl;
 
           if(col<shift*index)
-            *(ptr_data+elem(row+k,col))=*(ptr_data+elem(row+k+shift,col));
+            *(m_ptr_data+elem(row+k,col))=*(m_ptr_data+elem(row+k+shift,col));
           else
-            *(ptr_data+elem(row+k,col))=*(ptr_data+elem(row+k+shift,col+shift));
+            *(m_ptr_data+elem(row+k,col))=*(m_ptr_data+elem(row+k+shift,col+shift));
         }
       }
     }
@@ -512,12 +501,10 @@ int AlSymMat::RemoveCollsRows(std::vector<int> indices)
 {
   int n = indices.size();
   if (n==0) {
-    std::cerr<<"Vector of indices to remove is empty."<<std::endl;
-    return _size;
+    return m_size;
   }
-  if (n>_size) {
-    std::cerr<<"Vector of indices larger than matrix size."<<std::endl;
-    return _size;
+  if (n>m_size) {
+    throw std::range_error("Vector of indices larger than matrix size.");
   }
 
   // first sort the list of indices descending
@@ -533,22 +520,21 @@ int AlSymMat::RemoveCollsRows(std::vector<int> indices)
   // remove rows and columns starting from largest indices
   for (int i=0;i<n;i++) {
     int index = indices[i];
-    if (index > _size-1) {
-      std::cerr<<"Index "<<index<<" goes beyond matrix (size "<<_size<<")."<<std::endl;
-      continue;
+    if (index > m_size-1) {
+      throw std::out_of_range("AlSymMat::RemoveCollsRows: Index goes beyond matrix.");
     }
 
-    for (int j=index; j<_size-1; j++)
+    for (int j=index; j<m_size-1; j++)
       for (int col=0; col<=j; col++)
         if(col<index)
-          *(ptr_data+elem(j,col)) = *(ptr_data+elem(j+1,col));
+          *(m_ptr_data+elem(j,col)) = *(m_ptr_data+elem(j+1,col));
         else
-          *(ptr_data+elem(j,col)) = *(ptr_data+elem(j+1,col+1));
+          *(m_ptr_data+elem(j,col)) = *(m_ptr_data+elem(j+1,col+1));
 
-    _size--;
+    m_size--;
   }
 
-  return _size;
+  return m_size;
 }
 
 //______________________________________________________________________________
@@ -577,7 +563,7 @@ void AlSymMat::RemoveAlignPar(int index, int control)
         counterCol=0;
       }
 
-      *(ptr_data+elem(row,col))=*(ptr_data+elem(row+shiftRow,col+shiftCol));
+      *(m_ptr_data+elem(row,col))=*(m_ptr_data+elem(row+shiftRow,col+shiftCol));
 
       counterCol++;
       if (counterCol==5-control) {
@@ -639,7 +625,7 @@ StatusCode AlSymMat::Write(const std::string &filename, bool binary,
   if (!square) {
     for( int i=0; i<size(); i++) {
       for( int j=0; j<=i; j++) {
-        melem = *(ptr_data+(i+1)*i/2+j);
+        melem = *(m_ptr_data+(i+1)*i/2+j);
         if(binary)
           outmat.write((char*)&(melem), sizeof (melem));
         else
@@ -652,8 +638,8 @@ StatusCode AlSymMat::Write(const std::string &filename, bool binary,
   else {
     for( int i=0; i<size(); i++) {
       for( int j=0; j<size(); j++) {
-        if(i>=j) melem =  *(ptr_data+(i+1)*i/2+j);
-        else     melem =  *(ptr_data+(j+1)*j/2+i);
+        if(i>=j) melem =  *(m_ptr_data+(i+1)*i/2+j);
+        else     melem =  *(m_ptr_data+(j+1)*j/2+i);
         if(binary)
           outmat.write((char*)&(melem), sizeof (melem));
         else
@@ -709,7 +695,7 @@ StatusCode AlSymMat::Read(const std::string &filename, int &dofs,
   int32_t msiz=0;
   inmat.read((char*)&msiz, sizeof (msiz));
   dofs = abs(msiz);
-  _size = abs(msiz);
+  m_size = abs(msiz);
 
   if (m_StdUnits)
     inmat.read((char*)&version, sizeof (version));
@@ -722,7 +708,7 @@ StatusCode AlSymMat::Read(const std::string &filename, int &dofs,
       for(int j=0; j<msiz; j++) {
         inmat.read((char*)&melem, sizeof (melem));
         if( i>=j )
-          *(ptr_data+(i+1)*i/2+j) = melem;
+          *(m_ptr_data+(i+1)*i/2+j) = melem;
       }
     }
   }
@@ -733,7 +719,7 @@ StatusCode AlSymMat::Read(const std::string &filename, int &dofs,
     for( int i=0; i<msiz; i++) {
       for( int j=0; j<=i; j++) {
         inmat.read((char*)&melem, sizeof (melem));
-        *(ptr_data+(i+1)*i/2+j) = melem;
+        *(m_ptr_data+(i+1)*i/2+j) = melem;
       }
     }
   }
@@ -753,7 +739,7 @@ StatusCode AlSymMat::ReadProjected(const std::string &filename, int &dofs,
   int32_t msiz=0;
   inmat.read((char*)&msiz, sizeof (msiz));
   dofs = abs(msiz);
-  _size = abs(msiz);
+  m_size = abs(msiz);
 
   inmat.read((char*)&version, sizeof (version));
 
@@ -765,7 +751,7 @@ StatusCode AlSymMat::ReadProjected(const std::string &filename, int &dofs,
       for(int j=0; j<msiz; j++) {
         inmat.read((char*)&melem, sizeof (melem));
         if( i>=j )
-          *(ptr_data+(i+1)*i/2+j) = melem;
+          *(m_ptr_data+(i+1)*i/2+j) = melem;
       }
     }
   }
@@ -776,7 +762,7 @@ StatusCode AlSymMat::ReadProjected(const std::string &filename, int &dofs,
     for( int i=0; i<msiz; i++) {
       for( int j=0; j<=i; j++) {
         inmat.read((char*)&melem, sizeof (melem));
-        *(ptr_data+(i+1)*i/2+j) = melem;
+        *(m_ptr_data+(i+1)*i/2+j) = melem;
       }
     }
   }
@@ -789,18 +775,18 @@ StatusCode AlSymMat::ReadProjected(const std::string &filename, int &dofs,
 void AlSymMat::reSize(long int Nnew)
 {
   if ( Nnew != size() ) {
-    double * p = ptr_data;
+    double * p = m_ptr_data;
     long int ii;
 
-    _nele = Nnew*(Nnew+1)/2;
-    ptr_data = new double[_nele];
+    m_nele = Nnew*(Nnew+1)/2;
+    m_ptr_data = new double[m_nele];
     long int Size_old = size();
-    _size = Nnew;
+    m_size = Nnew;
     long int k = size() <= Size_old ? size() : Size_old;
     for( long int i=0; i<k; i++ ) {
       ii = (i+1)*i/2;
       for( long int j=0; j<=i; j++ )
-        *(ptr_data+ii+j) = *(p+ii+j);
+        *(m_ptr_data+ii+j) = *(p+ii+j);
     }
     delete [] p;
   }
@@ -811,26 +797,22 @@ double& AlSymMat::elemr(long int i,long int j)
 {
 #ifdef _DEBUG
   if( i<0 ) {
-    std::cerr << "AlSymMat::elemr: Index 1 < zero! " << i << std::endl;
-    return *(ptr_data);
+    throw std::underflow_error( "AlSymMat::elemr: Index 1 < zero!" );
   }
   if( i>=size() ) {
-    std::cerr << "AlSymMat::elemr: Index 1 too large! " << i << std::endl;
-    return *(ptr_data);
+    throw std::overflow_error( "AlSymMat::elemr: Index 1 too large!" );
   }
   if( j<0 ) {
-    std::cerr << "AlSymMat::elemr: Index 2 < zero! " << j << std::endl;
-    return *(ptr_data);
+    throw std::underflow_error( "AlSymMat::elemr: Index 2 < zero!" );
   }
   if( j>=size() ) {
-    std::cerr << "AlSymMat::elemr: Index 2 too large! " << j << std::endl;
-    return *(ptr_data);
+    throw std::overflow_error( "AlSymMat::elemr: Index 2 too large!" );
   }
 #endif
   if( j<=i )
-    return  *(ptr_data+(i+1)*i/2+j);
+    return  *(m_ptr_data+(i+1)*i/2+j);
 
-  return  *(ptr_data+(j+1)*j/2+i);
+  return  *(m_ptr_data+(j+1)*j/2+i);
 }
 
 //______________________________________________________________________________
@@ -838,26 +820,22 @@ double AlSymMat::elemc(long int i,long int j) const
 {
 #ifdef _DEBUG
   if( i<0 ) {
-    std::cerr << "AlSymMat::elemc: Index 1 < zero! " << i << std::endl;
-    return *(ptr_data);
+    throw std::underflow_error( "AlSymMat::elemc: Index 1 < zero!" );
   }
   if( i>=size() ) {
-    std::cerr << "AlSymMat::elemc: Index 1 too large! " << i << std::endl;
-    return *(ptr_data);
+    throw std::overflow_error( "AlSymMat::elemc: Index 1 too large!" );
   }
   if( j<0 ) {
-    std::cerr << "AlSymMat::elemc: Index 2 < zero! " << j << std::endl;
-    return *(ptr_data);
+    throw std::underflow_error( "AlSymMat::elemc: Index 2 < zero!" );
   }
   if( j>=size() ) {
-    std::cerr << "AlSymMat::elemc: Index 2 too large! " << j << std::endl;
-    return *(ptr_data);
+    throw std::overflow_error( "AlSymMat::elemc: Index 2 too large!" );
   }
 #endif
   if( j<=i )
-    return  *(ptr_data+(i+1)*i/2+j);
+    return  *(m_ptr_data+(i+1)*i/2+j);
 
-  return  *(ptr_data+(j+1)*j/2+i);
+  return  *(m_ptr_data+(j+1)*j/2+i);
 }
 
 
@@ -869,7 +847,7 @@ TMatrixDSparse* AlSymMat::makeTMatrix(){
 
   for( int i=0; i<si; i++ )
     for( int j=0; j<=i; j++ )
-			if( *(ptr_data+(i+1)*i/2+j)  != 0 ) 
+			if( *(m_ptr_data+(i+1)*i/2+j)  != 0 ) 
 				++nnz;
 
 	
@@ -880,10 +858,10 @@ TMatrixDSparse* AlSymMat::makeTMatrix(){
   
   for (int i=0;i<si;i++) {
     for (int j=0;j<=i;j++) {
-      if ( *(ptr_data+(i+1)*i/2+j) != 0.) {
+      if ( *(m_ptr_data+(i+1)*i/2+j) != 0.) {
         irow[n] = i;
         icol[n] = j;
-        val[n]  = *(ptr_data+(i+1)*i/2+j);
+        val[n]  = *(m_ptr_data+(i+1)*i/2+j);
         ++n;
       }
     }
diff --git a/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlSymMatBase.cxx b/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlSymMatBase.cxx
index 2da1a90f68a..f5c61b7a228 100644
--- a/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlSymMatBase.cxx
+++ b/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlSymMatBase.cxx
@@ -19,26 +19,26 @@ namespace Trk {
 //______________________________________________________________________________
 AlSymMatBase::AlSymMatBase()
 {
-  _matrix_type = 0;
-  _size = 0;
-  _nele = 0;
+  m_matrix_type = 0;
+  m_size = 0;
+  m_nele = 0;
 }
 
 //______________________________________________________________________________
 AlSymMatBase::AlSymMatBase(long int N)
 {
-  _matrix_type = 0;
-  _size = N;
-  _nele = N*(N+1)/2;
+  m_matrix_type = 0;
+  m_size = N;
+  m_nele = N*(N+1)/2;
 
 }
 
 //______________________________________________________________________________
 AlSymMatBase::AlSymMatBase(const AlSymMatBase & m)
 {
-  _matrix_type = 0;
-  _size      = m.size();
-  _nele      = m._nele;
+  m_matrix_type = 0;
+  m_size      = m.size();
+  m_nele      = m.m_nele;
 }
 
 //______________________________________________________________________________
@@ -46,8 +46,8 @@ AlSymMatBase & AlSymMatBase::operator=(const AlSymMatBase & m)
 {
   if (this==&m)
     return *this;
-  _matrix_type=m._matrix_type;
-  _nele=m._nele;
+  m_matrix_type=m.m_matrix_type;
+  m_nele=m.m_nele;
 
   return *this;
 }
diff --git a/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlVec.cxx b/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlVec.cxx
index 35b3893858a..ccdafd5be07 100644
--- a/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlVec.cxx
+++ b/Tracking/TrkAlignment/TrkAlgebraUtils/src/AlVec.cxx
@@ -12,34 +12,35 @@
 #include <fstream>
 #include <stdint.h>
 #include <math.h>
+#include <exception>
 
 namespace Trk {
 
 //______________________________________________________________________________
 AlVec::AlVec() {
-  _size = 0;
-  ptr_data = 0;  // set pointer to null
+  m_size = 0;
+  m_ptr_data = 0;  // set pointer to null
   m_pathbin="./";
   m_pathtxt="./";
 }
 
 //______________________________________________________________________________
 AlVec::AlVec(int N) {
-  _size = N;
-  ptr_data = new double[_size];
+  m_size = N;
+  m_ptr_data = new double[m_size];
   m_pathbin="./";
   m_pathtxt="./";
 
-  double*  p = ptr_data + _size;
-  while (p > ptr_data) *(--p) = 0.;
+  double*  p = m_ptr_data + m_size;
+  while (p > m_ptr_data) *(--p) = 0.;
 }
 
 //______________________________________________________________________________
 AlVec::AlVec(const AlVec& v) {
-  _size   = v._size;
+  m_size   = v.m_size;
   m_pathbin = v.m_pathbin;
   m_pathtxt = v.m_pathtxt;
-  ptr_data = new double[_size];
+  m_ptr_data = new double[m_size];
 
   copy(v);
 }
@@ -47,26 +48,25 @@ AlVec::AlVec(const AlVec& v) {
 //______________________________________________________________________________
 AlVec::~AlVec()
 {
-  delete [] ptr_data;
+  delete [] m_ptr_data;
 }
 
 //______________________________________________________________________________
 void AlVec::copy(const AlVec& v) {
-  if(_size!=v._size) {
-    std::cerr << "AlVec Assignment: size does not match!" << std::endl;
-    return;
+  if(m_size!=v.m_size) {
+    throw std::range_error( "AlVec Assignment: size does not match!" );
   }
 
-  double*  p = ptr_data + _size;
-  double*  q = v.ptr_data + _size;
-  while (p > ptr_data) *(--p) = (*(--q));
+  double*  p = m_ptr_data + m_size;
+  double*  q = v.m_ptr_data + m_size;
+  while (p > m_ptr_data) *(--p) = (*(--q));
 }
 
 //______________________________________________________________________________
 AlVec&  AlVec::operator=(const double& d) {
 
-  double*  p = ptr_data + _size;
-  while (p > ptr_data) *(--p) = d;
+  double*  p = m_ptr_data + m_size;
+  while (p > m_ptr_data) *(--p) = d;
 
   return *this;
 }
@@ -75,74 +75,69 @@ AlVec&  AlVec::operator=(const double& d) {
 AlVec& AlVec::operator=(const AlVec& v) {
   if (this==&v) return *this;
 
-  if(_size!=0 && _size!=v._size) {
-    std::cerr << "AlVec Assignment: size does not match!" << std::endl;
-    return *this;
+  if(m_size!=0 && m_size!=v.m_size) {
+    throw std::range_error( "AlVec Assignment: size does not match!" );
   }
 
-  if ( ptr_data != v.ptr_data )  copy(v);
+  if ( m_ptr_data != v.m_ptr_data )  copy(v);
 
   return *this;
 }
 
 //______________________________________________________________________________
 AlVec AlVec::operator+(const AlVec& v) const {
-  if( _size != v._size ) {
-    std::cerr << "operator+: vectors size does not match!" << std::endl;
-    return *this;
+  if( m_size != v.m_size ) {
+    throw std::range_error(  "operator+: vectors size does not match!" );
   }
 
-  AlVec b(_size);
+  AlVec b(m_size);
 
-  double*  p = ptr_data + _size;
-  double*  q = v.ptr_data + _size;
-  double*  r = b.ptr_data + _size;
-  while (p > ptr_data) *(--r) = (*(--p))+(*(--q));
+  double*  p = m_ptr_data + m_size;
+  double*  q = v.m_ptr_data + m_size;
+  double*  r = b.m_ptr_data + m_size;
+  while (p > m_ptr_data) *(--r) = (*(--p))+(*(--q));
 
   return b;
 }
 
 //______________________________________________________________________________
 AlVec& AlVec::operator+=(const AlVec& v) {
-  if( _size != v._size ) {
-    std::cerr << "operator+=: vectors size does not match!" << std::endl;
-    return *this;
+  if( m_size != v.m_size ) {
+    throw std::range_error(  "operator+=: vectors size does not match!" );
   }
 
-  double*  p = ptr_data + _size;
-  double*  q = v.ptr_data + _size;
-  while (p > ptr_data) *(--p) += (*(--q));
+  double*  p = m_ptr_data + m_size;
+  double*  q = v.m_ptr_data + m_size;
+  while (p > m_ptr_data) *(--p) += (*(--q));
 
   return *this;
 }
 
 //______________________________________________________________________________
 AlVec AlVec::operator-(const AlVec& v) const {
-  if( _size != v._size ) {
-    std::cerr << "operator-: vectors size does not match!" << std::endl;
-    return *this;
+  if( m_size != v.m_size ) {
+    throw std::range_error(  "operator-: vectors size does not match!" );
   }
 
-  AlVec b(_size);
+  AlVec b(m_size);
 
-  double*  p = ptr_data + _size;
-  double*  q = v.ptr_data + _size;
-  double*  r = b.ptr_data + _size;
-  while (p > ptr_data) *(--r) = (*(--p))-(*(--q));
+  double*  p = m_ptr_data + m_size;
+  double*  q = v.m_ptr_data + m_size;
+  double*  r = b.m_ptr_data + m_size;
+  while (p > m_ptr_data) *(--r) = (*(--p))-(*(--q));
 
   return b;
 }
 
 //______________________________________________________________________________
 AlVec& AlVec::operator-=(const AlVec& v) {
-  if( _size != v._size ) {
-    std::cerr << "operator+=: vectors size does not match!" << std::endl;
-    return *this;
+  if( m_size != v.m_size ) {
+    throw std::range_error(  "operator+=: vectors size does not match!" );
   }
 
-  double*  p = ptr_data + _size;
-  double*  q = v.ptr_data + _size;
-  while (p > ptr_data) *(--p) -= (*(--q));
+  double*  p = m_ptr_data + m_size;
+  double*  q = v.m_ptr_data + m_size;
+  while (p > m_ptr_data) *(--p) -= (*(--q));
 
   return *this;
 }
@@ -150,44 +145,41 @@ AlVec& AlVec::operator-=(const AlVec& v) {
 //______________________________________________________________________________
 double AlVec::operator*(const AlVec& v) const {
   double  b=0.;
-  if( _size != v._size ) {
-    std::cerr << "scalar product: vectors size does not match!" << std::endl;
-    return b;
+  if( m_size != v.m_size ) {
+    throw std::range_error(  "scalar product: vectors size does not match!" );
   }
 
-  double*  p = ptr_data + _size;
-  double*  q = v.ptr_data + _size;
-  while (p > ptr_data) b += *(--p)*(*(--q));
+  double*  p = m_ptr_data + m_size;
+  double*  q = v.m_ptr_data + m_size;
+  while (p > m_ptr_data) b += *(--p)*(*(--q));
 
   return b;
 }
 
 //______________________________________________________________________________
 AlVec AlVec::operator*(const AlMat& m) const {
-  if (_size != m.nrow()) {
-    std::cerr << "Left hand vector-matrix multiplication: size does not match!" << std::endl;
-    return *this;
+  if (m_size != m.nrow()) {
+    throw std::range_error(  "Left hand vector-matrix multiplication: size does not match!" );
   }
 
   AlVec b(m.ncol());
 
   for (int i=0;i<m.ncol();i++) {
-    for (int j=0;j<_size;j++)  *(b.ptr_data+i) += *(ptr_data+j)*m.elemc(j, i);
+    for (int j=0;j<m_size;j++)  *(b.m_ptr_data+i) += *(m_ptr_data+j)*m.elemc(j, i);
   }
   return b;
 }
 
 //______________________________________________________________________________
 AlVec  AlVec::operator*(const AlSymMatBase& m) const {
-  if (_size != m.size()) {
-    std::cerr << "Left hand vector-matrix multiplication: size does not match!" << std::endl;
-    return *this;
+  if (m_size != m.size()) {
+    throw std::range_error(  "Left hand vector-matrix multiplication: size does not match!" );
   }
 
-  AlVec b(_size);
+  AlVec b(m_size);
 
-  for (int i=0;i<_size;i++) {
-    for (int j=0;j<_size;j++)  *(b.ptr_data+i) += *(ptr_data+j)*m.elemc(j, i);
+  for (int i=0;i<m_size;i++) {
+    for (int j=0;j<m_size;j++)  *(b.m_ptr_data+i) += *(m_ptr_data+j)*m.elemc(j, i);
   }
   return b;
 }
@@ -195,8 +187,8 @@ AlVec  AlVec::operator*(const AlSymMatBase& m) const {
 //______________________________________________________________________________
 AlVec& AlVec::operator*=(const double& d) {
 
-  double*  p = ptr_data + _size;
-  while (p > ptr_data) *(--p) *= d;
+  double*  p = m_ptr_data + m_size;
+  while (p > m_ptr_data) *(--p) *= d;
 
   return *this;
 }
@@ -207,8 +199,8 @@ AlVec& AlVec::operator*=(const double& d) {
 double AlVec::norm() const {
 
   double result(0.);
-  double * p = ptr_data + _size;
-  while (p > ptr_data) {
+  double * p = m_ptr_data + m_size;
+  while (p > m_ptr_data) {
     --p;
     result += *p * *p ;
   }
@@ -218,16 +210,16 @@ double AlVec::norm() const {
 
 //______________________________________________________________________________
 void AlVec::reSize(int Nnew) {
-  if ( Nnew>=0 && Nnew != _size ) {
-    double*  p = ptr_data;
-    int _size_old = _size;
-    ptr_data = new double[Nnew];
-    _size = Nnew;
-    int k = _size <= _size_old ? _size : _size_old;
+  if ( Nnew>=0 && Nnew != m_size ) {
+    double*  p = m_ptr_data;
+    int m_size_old = m_size;
+    m_ptr_data = new double[Nnew];
+    m_size = Nnew;
+    int k = m_size <= m_size_old ? m_size : m_size_old;
 
     p += k;
-    double*  q = ptr_data + k;
-    while (q > ptr_data) *(--q) = *(--p);
+    double*  q = m_ptr_data + k;
+    while (q > m_ptr_data) *(--q) = *(--p);
 
     delete [] p;
   }
@@ -239,9 +231,9 @@ void AlVec::RemoveModule(int index) {
   const int shift=6;
   const int start = shift*index;
 
-  if(start<_size && index>=0){
-    for(int row=start; row<_size; row++) {
-      *(ptr_data+row)=*(ptr_data+row+shift);
+  if(start<m_size && index>=0){
+    for(int row=start; row<m_size; row++) {
+      *(m_ptr_data+row)=*(m_ptr_data+row+shift);
     }
   }
 }
@@ -256,8 +248,8 @@ void AlVec::RemoveAlignPar(int index, int control)
   index = index-control;
   //  std::cout << "index: " << index << " - control: " << control << std::endl;
 
-  for(int row=index; row<_size; row++) {
-    *(ptr_data+row) = *(ptr_data+row+shift);
+  for(int row=index; row<m_size; row++) {
+    *(m_ptr_data+row) = *(m_ptr_data+row+shift);
     counter++;
     if (counter==5-control) { counter=0; shift++; }
   }
@@ -268,12 +260,10 @@ int AlVec::RemoveElements(std::vector<int> indices)
 {
   int n = indices.size();
   if (n==0) {
-    std::cerr<<"Vector of indices to remove is empty."<<std::endl;
-    return _size;
+    return m_size;
   }
-  if (n>_size) {
-    std::cerr<<"Vector of indices larger than matrix size."<<std::endl;
-    return _size;
+  if (n>m_size) {
+    throw std::range_error( "Vector of indices larger than matrix size." );
   }
 
   // first sort the list of indices descending
@@ -289,18 +279,18 @@ int AlVec::RemoveElements(std::vector<int> indices)
   // remove elements starting from largest indices
   for (int i=0;i<n;i++) {
     int index = indices[i];
-    if (index > _size-1) {
-      std::cerr<<"Index "<<index<<" goes beyond matrix (size "<<_size<<")."<<std::endl;
+    if (index > m_size-1) {
+      throw std::out_of_range( "AlVec::RemoveElements: Index  goes beyond matrix " );
       continue;
     }
 
-    for (int j=index; j<_size-1; j++)
-      *(ptr_data+j) = *(ptr_data+j+1);
+    for (int j=index; j<m_size-1; j++)
+      *(m_ptr_data+j) = *(m_ptr_data+j+1);
 
-    _size--;
+    m_size--;
   }
 
-  return _size;
+  return m_size;
 }
 
 //______________________________________________________________________________
@@ -320,7 +310,7 @@ StatusCode AlVec::Write(const std::string &filename, bool binary, double scale,
                         std::map<int,unsigned long long> moduleIndexMap, float version)
 {
   std::ofstream outvec;
-  int32_t  io_size=_size;
+  int32_t  io_size=m_size;
 //  int32_t  io_scale=scale;
 
   if(binary) {
@@ -338,7 +328,7 @@ StatusCode AlVec::Write(const std::string &filename, bool binary, double scale,
     outvec.setf(std::ios::fixed);
     outvec.setf(std::ios::showpoint);
     outvec.precision(6);
-    outvec << "DoF: " << std::setw(6) << _size << std::endl;
+    outvec << "DoF: " << std::setw(6) << m_size << std::endl;
     outvec << "scale: " << std::setw(18) << scale << std::endl;
     outvec << "AlVec version: " << std::setw(6) << version << std::endl;
   }
@@ -347,13 +337,13 @@ StatusCode AlVec::Write(const std::string &filename, bool binary, double scale,
   double  velem=0;
   std::map<int,unsigned long long>::iterator itcod;
 
-  for( int i=0; i<_size; i++) {
+  for( int i=0; i<m_size; i++) {
     itcod=moduleIndexMap.find(i/6);
     if (itcod != moduleIndexMap.end()) {
       ielem = (itcod->second);
     }
 
-    velem = *(ptr_data+i);
+    velem = *(m_ptr_data+i);
 
     if(binary){
       outvec.write((char*)&(ielem), sizeof (ielem));
@@ -373,7 +363,7 @@ StatusCode AlVec::WritePartial(const std::string &filename, bool binary, double
   std::ofstream outvec;
 
   /*
-    int32_t  io_size=_size;
+    int32_t  io_size=m_size;
     //  int32_t  io_scale=scale;
     
     if(binary) {
@@ -392,7 +382,7 @@ StatusCode AlVec::WritePartial(const std::string &filename, bool binary, double
     outvec.setf(std::ios::fixed);
     outvec.setf(std::ios::showpoint);
     outvec.precision(6);
-    outvec << "DoF: " << std::setw(6) << _size << std::endl;
+    outvec << "DoF: " << std::setw(6) << m_size << std::endl;
     outvec << "scale: " << std::setw(18) << scale << std::endl;
     outvec << "AlVec version: " << std::setw(6) << version << std::endl;
     }
@@ -405,12 +395,12 @@ StatusCode AlVec::WritePartial(const std::string &filename, bool binary, double
   double velem=0;
   std::map<int,unsigned long long>::iterator itcod;
 
-  for( int i=0; i<_size; i++) {
+  for( int i=0; i<m_size; i++) {
     itcod=moduleIndexMap.find(i);
     if (itcod != moduleIndexMap.end())
       ielem = (itcod->second);
 
-    velem = *(ptr_data+i);
+    velem = *(m_ptr_data+i);
 
     if(binary){
       outvec.write((char*)&(ielem), sizeof (ielem));
@@ -436,13 +426,13 @@ StatusCode AlVec::WritePartial(const std::string &filename, bool binary, double
   double velem=0;
   std::map<int,std::string>::iterator itcod;
   
-  for( int i=0; i<_size; i++) {
+  for( int i=0; i<m_size; i++) {
     itcod=moduleNameMap.find(i);
     if (itcod != moduleNameMap.end()) {
       elem = (itcod->second);
     }
 
-    velem = *(ptr_data+i);
+    velem = *(m_ptr_data+i);
         
     if(binary){    
       std::cout<<"can't write module name in binary output!"<<std::endl;
@@ -459,7 +449,7 @@ StatusCode AlVec::WritePartial(const std::string &filename, bool binary, double
 StatusCode AlVec::InitializeOutputVector(const std::string& filename, bool binary, double scale, 
 					 float version, std::ofstream& outvec)
 {
-  int32_t  io_size=_size;
+  int32_t  io_size=m_size;
   
   if(binary) {
     outvec.open((m_pathbin+filename).c_str(), std::ios::binary);
@@ -477,7 +467,7 @@ StatusCode AlVec::InitializeOutputVector(const std::string& filename, bool binar
     outvec.setf(std::ios::fixed);
     outvec.setf(std::ios::showpoint);
     outvec.precision(6);
-    outvec << "DoF: " << std::setw(6) << _size << std::endl;
+    outvec << "DoF: " << std::setw(6) << m_size << std::endl;
     outvec << "scale: " << std::setw(18) << scale << std::endl;    
     outvec << "AlVec version: " << std::setw(6) << version << std::endl;    
   }
@@ -490,19 +480,19 @@ StatusCode AlVec::ReadPartial(const std::string &filename, double &scale,
 {
   bool m_StdUnits = true;
   if (StatusCode::SUCCESS != CheckVecVersion(m_pathbin+filename, m_StdUnits)) {
-    std::cout<<"CheckVecVersion failed"<<std::endl;
+    //std::cout<<"CheckVecVersion failed"<<std::endl;
     return StatusCode::FAILURE;
   }
 
   std::ifstream invec((m_pathbin+filename).c_str(), std::ios::binary);
   if(invec.fail()) {
-    std::cout<<"ifstream failed"<<std::endl;
+    //std::cout<<"ifstream failed"<<std::endl;
     return StatusCode::FAILURE;
   }
 
   int32_t vsiz=0;
   invec.read((char*)&vsiz, sizeof (vsiz));
-  _size  = vsiz;
+  m_size  = vsiz;
 
 //  int32_t io_scale;
   invec.read((char*)&scale, sizeof(scale));
@@ -517,12 +507,12 @@ StatusCode AlVec::ReadPartial(const std::string &filename, double &scale,
 
   int64_t ielem=0;
   double  velem=0.0;
-  for( int i=0; i<_size; i++) {
+  for( int i=0; i<m_size; i++) {
     invec.read((char*)&ielem, sizeof (ielem));
     modmap[i] = ielem;
 
     invec.read((char*)&velem, sizeof (velem));
-    *(ptr_data+i) = velem;
+    *(m_ptr_data+i) = velem;
 
     // std::cout << "AlVec (" << i << "):: " << ielem << ", " << velem << std::endl;
 
@@ -568,7 +558,7 @@ StatusCode AlVec::Read(const std::string &filename, double &scale,
 
   bool m_StdUnits = true;
   if (StatusCode::SUCCESS != CheckVecVersion(m_pathbin+filename, m_StdUnits)) {
-    std::cout<<"CheckVecVersion failed"<<std::endl;
+    //std::cout<<"CheckVecVersion failed"<<std::endl;
     return StatusCode::FAILURE;
   }
 
@@ -578,8 +568,8 @@ StatusCode AlVec::Read(const std::string &filename, double &scale,
 
   int32_t vsiz=0;
   invec.read((char*)&vsiz, sizeof (vsiz));
-  _size  = vsiz;
-//  std::cout<<"size="<<_size<<std::endl;
+  m_size  = vsiz;
+//  std::cout<<"size="<<m_size<<std::endl;
 
 //  int32_t io_scale;
   invec.read((char*)&scale, sizeof (scale));
@@ -595,12 +585,12 @@ StatusCode AlVec::Read(const std::string &filename, double &scale,
 
   int64_t ielem=0;
   double  velem=0.0;
-  for( int i=0; i<_size; i++) {
+  for( int i=0; i<m_size; i++) {
     invec.read((char*)&ielem, sizeof (ielem));
     modmap[i/6] = ielem;
 
     invec.read((char*)&velem, sizeof (velem));
-    *(ptr_data+i) = velem;
+    *(m_ptr_data+i) = velem;
 
     // std::cout << "AlVec (" << i << "):: " << ielem << ", " << velem << std::endl;
 
@@ -620,7 +610,7 @@ StatusCode AlVec::ReadProjected(const std::string &filename, double &scale,
 
   int32_t vsiz=0;
   invec.read((char*)&vsiz, sizeof (vsiz));
-  _size  = vsiz;
+  m_size  = vsiz;
 
 //  int32_t io_scale;
   invec.read((char*)&scale, sizeof (scale));
@@ -633,12 +623,12 @@ StatusCode AlVec::ReadProjected(const std::string &filename, double &scale,
 
   int64_t ielem=0;
   double  velem=0.0;
-  for( int i=0; i<_size; i++) {
+  for( int i=0; i<m_size; i++) {
     invec.read((char*)&ielem, sizeof (ielem));
     modmap[i/6] = ielem;
 
     invec.read((char*)&velem, sizeof (velem));
-    *(ptr_data+i) = velem;
+    *(m_ptr_data+i) = velem;
 
     // std::cout << "AlVec (" << i << "):: " << ielem << ", " << velem << std::endl;
 
@@ -656,13 +646,13 @@ StatusCode AlVec::ReadScalaPack(const std::string &filename){
 
   int32_t vsiz=0;
   eigenvec.read((char*)&vsiz, sizeof (vsiz));
-  _size=vsiz;
+  m_size=vsiz;
 
   double velem=0;
-  for( int i=0; i<_size; i++) {
+  for( int i=0; i<m_size; i++) {
     eigenvec.read((char*)&velem, sizeof (velem));
     // printf("v[%d] = %.16lf \n",i,velem);
-    *(ptr_data+i) = velem;
+    *(m_ptr_data+i) = velem;
   }
 
   eigenvec.close();
@@ -673,7 +663,7 @@ StatusCode AlVec::ReadScalaPack(const std::string &filename){
 StatusCode AlVec::WriteEigenvalueVec(const std::string &filename, bool binary){
 
   std::ofstream outvec;
-  int32_t  io_size=_size;
+  int32_t  io_size=m_size;
 
   if(binary) {
     outvec.open((m_pathbin+filename).c_str(), std::ios::binary);
@@ -692,16 +682,16 @@ StatusCode AlVec::WriteEigenvalueVec(const std::string &filename, bool binary){
     outvec.setf(std::ios::fixed);
     outvec.setf(std::ios::showpoint);
     outvec.precision(6);
-    outvec << "AlVec::DoF: " << std::setw(6) << _size << std::endl;
+    outvec << "AlVec::DoF: " << std::setw(6) << m_size << std::endl;
   }
 
   //jlove int32_t ielem=0;
   double  velem=0;
 //  std::map<int,int>::iterator itcod;
 
-  for( int i=0; i<_size; i++) {
+  for( int i=0; i<m_size; i++) {
 
-    velem = *(ptr_data+i);
+    velem = *(m_ptr_data+i);
 
     if(binary)
       outvec.write((char*)&(velem), sizeof (velem));
diff --git a/Tracking/TrkAlignment/TrkAlgebraUtils/src/IntVec.cxx b/Tracking/TrkAlignment/TrkAlgebraUtils/src/IntVec.cxx
index d52767e7b63..b9dbfde12b1 100644
--- a/Tracking/TrkAlignment/TrkAlgebraUtils/src/IntVec.cxx
+++ b/Tracking/TrkAlignment/TrkAlgebraUtils/src/IntVec.cxx
@@ -6,48 +6,48 @@
 #include "TrkAlgebraUtils/IntVec.h"
 #include <iostream>
 #include <stdint.h>
+#include <exception>
 
 namespace Trk {
 
 IntVec::IntVec()
-  : Nele(0), ptr_data(0)
+  : m_Nele(0), m_ptr_data(0)
 {}
 
 IntVec::IntVec(int N) {
-  Nele = N;
-  ptr_data = new int[Nele];
-  for(int i=0;i<Nele;i++)
-    *(ptr_data+i)=0;
+  m_Nele = N;
+  m_ptr_data = new int[m_Nele];
+  for(int i=0;i<m_Nele;i++)
+    *(m_ptr_data+i)=0;
 }
 
 IntVec::IntVec(int N, int init) {
-  Nele = N;
-  ptr_data = new int[Nele];
-  for(int i=0;i<Nele;i++)
-    *(ptr_data+i)=init;
+  m_Nele = N;
+  m_ptr_data = new int[m_Nele];
+  for(int i=0;i<m_Nele;i++)
+    *(m_ptr_data+i)=init;
 }
 
 IntVec::IntVec(const IntVec& v) {
-  Nele = v.Nele;
-  ptr_data = new int[Nele];
-  for(int i=0;i<Nele;i++)
-    *(ptr_data+i)=v[i];
+  m_Nele = v.m_Nele;
+  m_ptr_data = new int[m_Nele];
+  for(int i=0;i<m_Nele;i++)
+    *(m_ptr_data+i)=v[i];
 }
 
 IntVec::~IntVec(){
-  delete [] ptr_data;
+  delete [] m_ptr_data;
 }
 
 IntVec& IntVec::operator=(const IntVec& v) {
-  if(Nele!=0 && Nele!=v.Nele) {
-    std::cerr << "IntVec Assignment: size does not match!" << std::endl;
-    return *this;
+  if(m_Nele!=0 && m_Nele!=v.m_Nele) {
+    throw std::range_error( "IntVec Assignment: size does not match!" );
   }
 
-  if ( ptr_data != v.ptr_data ) {
-    reSize(v.Nele);
-    for(int i=0;i<Nele;i++)
-      *(ptr_data+i)=v[i];
+  if ( m_ptr_data != v.m_ptr_data ) {
+    reSize(v.m_Nele);
+    for(int i=0;i<m_Nele;i++)
+      *(m_ptr_data+i)=v[i];
   }
 
   return *this;
@@ -55,92 +55,84 @@ IntVec& IntVec::operator=(const IntVec& v) {
 
 int& IntVec::operator[](int i) {
   if(i<0) {
-    std::cerr << "IntVec: Index < zero! " << std::endl;
-    return ptr_data[0];
+    throw std::out_of_range( "IntVec: Index < zero! " );
   }
-  else if(i>=Nele) {
-    std::cerr << "IntVec: Index too large! " << std::endl;
-    return ptr_data[0];
+  else if(i>=m_Nele) {
+    throw std::out_of_range( "IntVec: Index too large! " );
   }
 
-  return *(ptr_data+i);
+  return *(m_ptr_data+i);
 }
 
 const int& IntVec::operator[](int i) const {
   if(i<0) {
-    std::cerr << "IntVec: Index < zero! " << std::endl;
-    return ptr_data[0];
+    throw std::out_of_range( "IntVec: Index < zero! " );
   }
-  else if(i>=Nele) {
-    std::cerr << "IntVec: Index too large! " << std::endl;
-    return ptr_data[0];
+  else if(i>=m_Nele) {
+    throw std::out_of_range( "IntVec: Index too large! " );
   }
 
-  return  *(ptr_data+i);
+  return  *(m_ptr_data+i);
 }
 
 IntVec IntVec::operator+(const IntVec& v) {
-  if( Nele != v.Nele ) {
-    std::cerr << "operator+: vectors size does not match!" << std::endl;
-    return *this;
+  if( m_Nele != v.m_Nele ) {
+    throw std::range_error( "operator+: vectors size does not match!" );
   }
 
-  IntVec b(Nele);
-  for (int i=0;i<Nele;i++)
-    b[i] = *(ptr_data+i) + v[i];
+  IntVec b(m_Nele);
+  for (int i=0;i<m_Nele;i++)
+    b[i] = *(m_ptr_data+i) + v[i];
 
   return b;
 }
 
 IntVec& IntVec::operator+=(const IntVec& v) {
-  if( Nele != v.Nele ) {
-    std::cerr << "operator+=: vectors size does not match!" << std::endl;
-    return *this;
+  if( m_Nele != v.m_Nele ) {
+    throw std::range_error(  "operator+=: vectors size does not match!" );
   }
 
-  for (int i=0;i<Nele;i++)
-    *(ptr_data+i)+=v[i];
+  for (int i=0;i<m_Nele;i++)
+    *(m_ptr_data+i)+=v[i];
 
   return *this;
 }
 
 IntVec  IntVec::operator-(const IntVec& v) {
-  if( Nele != v.Nele ) {
-    std::cerr << "operator+: vectors size does not match!" << std::endl;
-    return *this;
+  if( m_Nele != v.m_Nele ) {
+    throw std::range_error(  "operator+: vectors size does not match!" );
   }
 
-  IntVec b(Nele);
-  for (int i=0;i<Nele;i++)
-    b[i] = *(ptr_data+i) - v[i];
+  IntVec b(m_Nele);
+  for (int i=0;i<m_Nele;i++)
+    b[i] = *(m_ptr_data+i) - v[i];
 
   return b;
 }
 
 IntVec& IntVec::operator-=(const IntVec& v) {
-  if( Nele != v.Nele ) {
-    std::cerr << "operator+=: vectors size does not match!" << std::endl;
-    return *this;
+  if( m_Nele != v.m_Nele ) {
+    throw std::range_error(  "operator+=: vectors size does not match!" );
   }
 
-  for (int i=0;i<Nele;i++)
-    *(ptr_data+i)-=v[i];
+  for (int i=0;i<m_Nele;i++)
+    *(m_ptr_data+i)-=v[i];
 
   return *this;
 }
 
 
 void IntVec::reSize(int Nnew) {
-  if ( Nnew>=0 && Nnew != Nele ) {
-    int*  p = ptr_data;
-    int Nele_old = Nele;
-    ptr_data = new int[Nnew];
-    Nele = Nnew;
-    int k = Nele <= Nele_old ? Nele : Nele_old;
+  if ( Nnew>=0 && Nnew != m_Nele ) {
+    int*  p = m_ptr_data;
+    int Nele_old = m_Nele;
+    m_ptr_data = new int[Nnew];
+    m_Nele = Nnew;
+    int k = m_Nele <= Nele_old ? m_Nele : Nele_old;
 
     p += k;
-    int*  q = ptr_data + k;
-    while (q > ptr_data)
+    int*  q = m_ptr_data + k;
+    while (q > m_ptr_data)
       *(--q) = *(--p);
 
     delete [] p;
diff --git a/Tracking/TrkAlignment/TrkAlignGenTools/TrkAlignGenTools/MatrixTool.h b/Tracking/TrkAlignment/TrkAlignGenTools/TrkAlignGenTools/MatrixTool.h
index b3cf31e85c9..26a6bdb38b7 100644
--- a/Tracking/TrkAlignment/TrkAlignGenTools/TrkAlignGenTools/MatrixTool.h
+++ b/Tracking/TrkAlignment/TrkAlignGenTools/TrkAlignGenTools/MatrixTool.h
@@ -58,9 +58,9 @@ namespace Trk {
     enum SolveOption { 
       NONE                 = 0, //!< not solve in any case (to be used when ipc)
       SOLVE                = 1, //!< solving after data accumulation (LAPACK)
-      SOLVE_FAST           = 2, //!< Fast (***REMOVED*** method) solving after data accumulation
+      SOLVE_FAST           = 2, //!< Fast (Eigen method) solving after data accumulation
       DIRECT_SOLVE         = 3, //!< direct solving (LAPACK), already available matrix & vector
-      DIRECT_SOLVE_FAST    = 4, //!< direct Fast (***REMOVED*** method) solving, already available matrix & vector
+      DIRECT_SOLVE_FAST    = 4, //!< direct Fast (Eigen method) solving, already available matrix & vector
       DIRECT_SOLVE_CLUSTER = 5, //!< computation of alignment parameters from SCALAPAK already solved matrix
       SOLVE_ROOT           = 6, //!< computation using ROOT
       SOLVE_CLHEP          = 7  //!< computation using CLHEP
@@ -128,7 +128,7 @@ namespace Trk {
     int solveROOT();
     int solveCLHEP();
     int solveLapack();
-    int solve***REMOVED***();
+    int solveSparseEigen();
     int solveLocal();
 
     StatusCode spuriousRemoval();
diff --git a/Tracking/TrkAlignment/TrkAlignGenTools/src/MatrixTool.cxx b/Tracking/TrkAlignment/TrkAlignGenTools/src/MatrixTool.cxx
index a477f5f2540..b7c9846b05c 100644
--- a/Tracking/TrkAlignment/TrkAlignGenTools/src/MatrixTool.cxx
+++ b/Tracking/TrkAlignment/TrkAlignGenTools/src/MatrixTool.cxx
@@ -1426,7 +1426,7 @@ namespace Trk {
 
       case SOLVE_FAST:
       case DIRECT_SOLVE_FAST:
-        info = solve***REMOVED***();
+        info = solveSparseEigen();
         break;
 
       default:
@@ -2147,13 +2147,13 @@ namespace Trk {
   }
 
   //________________________________________________________________________
-  int MatrixTool::solve***REMOVED***()
+  int MatrixTool::solveSparseEigen()
   {
-    ATH_MSG_INFO("solving Global using ***REMOVED***");
+    ATH_MSG_INFO("solving Global using SparseEigen");
     if(m_logStream) {
       *m_logStream<<"*************************************************************"<<std::endl;
       *m_logStream<<"**************  solving using Global method  ****************"<<std::endl;
-      *m_logStream<<"**************          using ***REMOVED***           ****************"<<std::endl;
+      *m_logStream<<"**************      using SparseEigen        ****************"<<std::endl;
       *m_logStream<<"*************************************************************"<<std::endl;
     }
 
@@ -2205,17 +2205,17 @@ namespace Trk {
     ATH_MSG_DEBUG("running the solving");
 
     // solve
-    int info = (*aBetterMat).***REMOVED***Solve(*aBetterVec);
+    int info = (*aBetterMat).SolveWithEigen(*aBetterVec);
 
     if(info == 0) {
-      ATH_MSG_INFO("***REMOVED*** solving OK");
+      ATH_MSG_INFO("SolveWithEigen solving OK");
       if(m_logStream)
-        *m_logStream<<"***REMOVED*** solving OK."<<std::endl;
+        *m_logStream<<"SolveWithEigen solving OK."<<std::endl;
     }
     else {
-      msg(MSG::ERROR)<<"***REMOVED*** error code (0 if OK) = "<<info<<endreq;
+      msg(MSG::ERROR)<<"SolveWithEigen error code (0 if OK) = "<<info<<endreq;
       if(m_logStream)
-        *m_logStream<<"***REMOVED*** error code (0 if OK) = "<<info<<std::endl;
+        *m_logStream<<"SolveWithEigen error code (0 if OK) = "<<info<<std::endl;
     }
 
     if( isCopy )
@@ -2225,11 +2225,11 @@ namespace Trk {
     // stop measuring time
     clock_t stoptime = clock();
     double totaltime = (stoptime-starttime)/double(CLOCKS_PER_SEC);
-    ATH_MSG_INFO("Time spent in solve***REMOVED***: "<<totaltime<<" s");
+    ATH_MSG_INFO("Time spent in SolveWithEigen: "<<totaltime<<" s");
 
     ATH_MSG_DEBUG("Alignment constants:");
     // compute alignment corrections (translations in mm and rotations in rad)
-    // for ***REMOVED*** variances are not calculated
+    // for solveSparseEigen variances are not calculated
     for(int i=0; i<m_aNDoF; i++) {
 
       double param = -(*aBetterVec)[i];
@@ -2273,7 +2273,7 @@ namespace Trk {
     delete aBetterVec;
 
     // need to do this since success return value from Lapack is 0
-    // and from solve***REMOVED***() it is 1
+    // and from SolveWithEigen() it is 1
     if (info==0)
       info = 1;
 
diff --git a/Tracking/TrkAlignment/TrkAlignInterfaces/TrkAlignInterfaces/IDerivCalcTool.h b/Tracking/TrkAlignment/TrkAlignInterfaces/TrkAlignInterfaces/IDerivCalcTool.h
index d5bb42b7581..5d80d02c73f 100644
--- a/Tracking/TrkAlignment/TrkAlignInterfaces/TrkAlignInterfaces/IDerivCalcTool.h
+++ b/Tracking/TrkAlignment/TrkAlignInterfaces/TrkAlignInterfaces/IDerivCalcTool.h
@@ -41,9 +41,9 @@ namespace Trk {
     enum SolveOption { 
       NONE                 = 0, //!< not solve in any case (to be used when ipc)
       SOLVE                = 1, //!< solving after data accumulation (LAPACK)
-      SOLVE_FAST           = 2, //!< Fast (***REMOVED*** method) solving after data accumulation
+      SOLVE_FAST           = 2, //!< Fast (Eigen method) solving after data accumulation
       DIRECT_SOLVE         = 3, //!< direct solving (LAPACK), already available matrix & vector
-      DIRECT_SOLVE_FAST    = 4, //!< direct Fast (***REMOVED*** method) solving, already available matrix & vector
+      DIRECT_SOLVE_FAST    = 4, //!< direct Fast (Eigen method) solving, already available matrix & vector
       DIRECT_SOLVE_CLUSTER = 5  //!< computation of alignment parameters from SCALAPAK already solved matrix
     }; // this is also defined in TrkGlobAlign class
 
-- 
GitLab