diff --git a/Alignment/AlignKernel/AlignKernel/AlMat.h b/Alignment/AlignKernel/AlignKernel/AlMat.h
index ad79b6a32d10349444cf772a5d1a5706f57961ca..fe29b9db4de105276d3245d0b3ca1386715ac23b 100644
--- a/Alignment/AlignKernel/AlignKernel/AlMat.h
+++ b/Alignment/AlignKernel/AlignKernel/AlMat.h
@@ -1,5 +1,5 @@
/*****************************************************************************\
-* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration *
+* (c) Copyright 2000-2019 CERN for the benefit of the LHCb Collaboration *
* *
* This software is distributed under the terms of the GNU General Public *
* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
@@ -8,131 +8,21 @@
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
-#ifndef SOLVKERNEL_ALMAT_H
-#define SOLVKERNEL_ALMAT_H
-#include <iomanip>
-#include <iostream>
+#pragma once
-class AlVec;
-class AlSymMat;
-
-class AlMat
-
-{
+#pragma GCC diagnostic push
+#if !defined( __clang__ )
+# pragma GCC diagnostic ignored "-Wdeprecated-copy"
+#endif
+#include <Eigen/Dense>
+#pragma GCC diagnostic pop
+/// Basically just a wrapper around the Eigen class, but with zero initialization
+class AlMat : public Eigen::MatrixXd {
public:
- AlMat( size_t N, size_t M );
- AlMat();
- AlMat( const AlMat& m );
- ~AlMat();
-
- class AlMat_row {
- public:
- inline AlMat_row( AlMat&, size_t );
- double& operator[]( size_t );
-
- private:
- AlMat& _a;
- size_t _r;
- };
- class AlMat_row_const {
- public:
- inline AlMat_row_const( const AlMat&, size_t );
- const double& operator[]( size_t ) const;
-
- private:
- const AlMat& _a;
- size_t _r;
- };
- // helper class to implement m[i][j]
-
- inline AlMat_row operator[]( size_t );
- inline AlMat_row_const operator[]( size_t ) const;
- double operator()( size_t irow, size_t icol ) const { return ( *this )[irow][icol]; }
- double& operator()( size_t irow, size_t icol ) { return ( *this )[irow][icol]; }
-
- AlMat& operator=( const AlMat& m );
- AlMat& operator=( const AlSymMat& m );
- AlMat operator+( const AlMat& m );
- AlMat operator+( const AlMat& m ) const;
- AlMat& operator+=( const AlMat& m );
- AlMat operator-( const AlMat& m );
- AlMat operator-( const AlMat& m ) const;
- AlMat& operator-=( const AlMat& m );
- AlMat operator*( const AlMat& m );
- AlMat operator*( const AlMat& m ) const;
- AlMat operator*( const AlSymMat& m );
- AlMat operator*( const AlSymMat& m ) const;
- AlVec operator*( const AlVec& v );
- AlVec operator*( const AlVec& v ) const;
- AlMat& operator*=( const double& d );
-
- // transposed matrix
- AlMat T() const;
- //
- // invert sym matrix declared as non-symetric for convenience
-
- void invertS( int& ierr, double Norm );
-
- size_t nrow() const;
- size_t ncol() const;
- // reSize method
- void reSize( size_t Nnew, size_t Mnew );
- void removeRow( size_t Nr );
- void removeCol( size_t Nr );
-
- void copyS( const AlSymMat& m );
-
- // void invert(int ierr); // only if Ncol=Nrow !!
-
-private:
- friend class AlVec;
- friend class AlSymMat;
- friend class AlMat_row;
- friend class AlMat_row_const;
-
- void copy( const AlMat& m );
- // void copyS(const AlSymMat& m);
-
- size_t Ncol, Nrow;
- size_t Nele;
- double* ptr_data;
-
- // double* mat_p;
- // double** mat;
+ using Base = Eigen::MatrixXd;
+ AlMat( Index n, Index m ) : Base( Base::Zero( n, m ) ) {}
+ AlMat( size_t n, size_t m ) : Base( Base::Zero( n, m ) ) {}
+ using Base::Base;
};
-
-// inline operators
-
-inline AlMat::AlMat_row AlMat::operator[]( size_t r ) {
- AlMat_row b( *this, r );
- return b;
-}
-
-inline AlMat::AlMat_row_const AlMat::operator[]( size_t r ) const
-
-{
- const AlMat_row_const b( *this, r );
-
- return b;
-}
-
-inline double& AlMat::AlMat_row::operator[]( size_t c ) {
- if ( _r >= _a.Nrow || c >= _a.Ncol ) std::cerr << "Range error in AlMat::operator[][]" << std::endl;
-
- return *( _a.ptr_data + _r * _a.Ncol + c );
-}
-
-inline const double& AlMat::AlMat_row_const::operator[]( size_t c ) const {
-
- if ( _r >= _a.Nrow || c >= _a.Ncol ) std::cerr << "Range error in AlMat::operator[][]" << std::endl;
-
- return *( _a.ptr_data + _r * _a.Ncol + c );
-}
-
-inline AlMat::AlMat_row::AlMat_row( AlMat& a, size_t r ) : _a( a ), _r( r ) {}
-
-inline AlMat::AlMat_row_const::AlMat_row_const( const AlMat& a, size_t r ) : _a( a ), _r( r ) {}
-
-#endif // SOLVKERNEL_ALMAT_H
diff --git a/Alignment/AlignKernel/AlignKernel/AlSymMat.h b/Alignment/AlignKernel/AlignKernel/AlSymMat.h
index af050b29b75efb717f7c51d9c87524e82aef958b..28560205ed558f4e61bdbd2cd4f99bbcdac0b4ac 100644
--- a/Alignment/AlignKernel/AlignKernel/AlSymMat.h
+++ b/Alignment/AlignKernel/AlignKernel/AlSymMat.h
@@ -8,155 +8,93 @@
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
-#ifndef SOLVKERNEL_ALSYMMAT_H
-#define SOLVKERNEL_ALSYMMAT_H
+
+#pragma once
// from STD
#include <iomanip>
#include <iostream>
-#include <ostream>
#include <vector>
-class AlVec;
-class AlMat;
+#include "AlignKernel/AlMat.h"
+#include "AlignKernel/AlVec.h"
class AlSymMat {
public:
- AlSymMat( size_t N );
- AlSymMat();
- AlSymMat( const AlSymMat& m );
- ~AlSymMat();
-
- class AlSymMat_row {
- public:
- inline AlSymMat_row( AlSymMat&, size_t );
- double& operator[]( size_t );
-
- private:
- AlSymMat& _a;
- size_t _r;
- };
- class AlSymMat_row_const {
- public:
- inline AlSymMat_row_const( const AlSymMat&, size_t );
- const double& operator[]( size_t ) const;
-
- private:
- const AlSymMat& _a;
- size_t _r;
- };
- // helper class to implement m[i][j]
-
- inline AlSymMat_row operator[]( size_t );
- inline AlSymMat_row_const operator[]( size_t ) const;
+ AlSymMat( size_t N ) { reSize( N ); }
+ AlSymMat() = default;
+ AlSymMat( const AlSymMat& m ) = default;
+ AlSymMat( AlSymMat&& m ) = default;
+ AlSymMat& operator=( const AlSymMat& m ) = default;
+ AlSymMat& operator=( AlSymMat&& m ) = default;
// access method that requires irow>=icol without checking
- double fast( size_t irow, size_t icol ) const { return *( ptr_data + ( irow + 1 ) * irow / 2 + icol ); }
- double& fast( size_t irow, size_t icol ) { return *( ptr_data + ( irow + 1 ) * irow / 2 + icol ); }
+ double fast( size_t irow, size_t icol ) const { return m_data[( irow + 1 ) * irow / 2 + icol]; }
+ double& fast( size_t irow, size_t icol ) { return m_data[( irow + 1 ) * irow / 2 + icol]; }
double operator()( size_t irow, size_t icol ) const { return irow >= icol ? fast( irow, icol ) : fast( icol, irow ); }
double& operator()( size_t irow, size_t icol ) { return irow >= icol ? fast( irow, icol ) : fast( icol, irow ); }
- AlSymMat& operator=( const AlSymMat& m );
- AlSymMat operator+( const AlSymMat& m );
- AlSymMat operator+( const AlSymMat& m ) const;
AlSymMat& operator+=( const AlSymMat& m );
- AlSymMat operator-( const AlSymMat& m );
- AlSymMat operator-( const AlSymMat& m ) const;
- AlMat operator*( const AlSymMat& m );
- AlMat operator*( const AlSymMat& m ) const;
- AlVec operator*( const AlVec& v );
- AlVec operator*( const AlVec& v ) const;
- AlMat operator*( const AlMat& m );
- AlMat operator*( const AlMat& m ) const;
+ AlSymMat& operator-=( const AlSymMat& m );
AlSymMat& operator*=( const double& d );
+ AlSymMat operator-( const AlSymMat& m ) const;
+ AlSymMat operator+( const AlSymMat& m ) const;
// insert block from SMatrix types
template <typename T>
void insert( unsigned col, unsigned row, const T& rhs ) {
for ( unsigned i = 0; i < T::kCols; ++i )
- for ( unsigned j = 0; j < T::kRows; ++j ) ( *this )[col + i][row + j] = rhs( i, j );
+ for ( unsigned j = 0; j < T::kRows; ++j ) ( *this )( col + i, row + j ) = rhs( i, j );
}
// reSize method
- void reSize( size_t Nnew );
- size_t size() const { return Size; }
- size_t nrow() const { return size(); }
- size_t ncol() const { return size(); }
-
- void removeElement( size_t Nr );
- void ShrinkDofs( std::vector<bool>& mask );
+ void reSize( size_t N ) {
+ m_N = N;
+ m_data.resize( N * ( N + 1 ) / 2, 0.0 );
+ }
+ size_t size() const { return m_N; }
+ size_t rows() const { return m_N; }
+ size_t cols() const { return m_N; }
- void invert( int& ierr );
- double determinant();
- int diagonalize( char jobz, AlVec& w, AlMat& z );
- void diagonalize_GSL( AlVec& egval, AlMat& egvec );
+ int diagonalize_GSL( AlVec& egval, AlMat& egvec ) const;
+ double determinant() const;
+ inline AlMat toMatrix() const;
+ // explicit operator AlMat() const { return toMatrix() ; }
private:
- friend class AlVec;
- friend class AlMat;
- friend class AlSymMat_row;
- friend class AlSymMat_row_const;
-
- void copy( const AlSymMat& m );
- double* ptr_data;
- size_t Size;
- size_t Nele;
+ std::vector<double> m_data;
+ size_t m_N{0};
};
// inline operators
inline std::ostream& operator<<( std::ostream& lhs, const AlSymMat& rhs ) {
lhs << "[ ";
for ( size_t i = 0u; i < rhs.size(); ++i ) {
- for ( size_t j = 0u; j < rhs.size(); ++j ) { lhs << " " << rhs[i][j] << " "; }
+ for ( size_t j = 0u; j < rhs.size(); ++j ) { lhs << " " << rhs( i, j ) << " "; }
lhs << std::endl;
}
lhs << " ]" << std::endl;
return lhs;
}
-inline AlSymMat::AlSymMat_row AlSymMat::operator[]( size_t r ) {
- AlSymMat_row b( *this, r );
- return b;
-}
-
-inline AlSymMat::AlSymMat_row_const AlSymMat::operator[]( size_t r ) const
-
-{
- const AlSymMat_row_const b( *this, r );
-
- return b;
-}
-
-inline double& AlSymMat::AlSymMat_row::operator[]( size_t c ) {
- if ( _r >= _a.Size || c >= _a.Size ) std::cerr << "Range error in AlSymMat::operator[][]" << std::endl;
-
- if ( _r >= c ) {
- return *( _a.ptr_data + ( _r + 1 ) * _r / 2 + c );
- } else {
- return *( _a.ptr_data + ( c + 1 ) * c / 2 + _r );
- }
+inline AlMat AlSymMat::toMatrix() const {
+ const auto& matrix = *this;
+ const auto N = m_N;
+ Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> rc{N, N};
+ for ( size_t irow = 0; irow < N; ++irow )
+ for ( size_t icol = 0; icol <= irow; ++icol ) rc( irow, icol ) = rc( icol, irow ) = matrix.fast( irow, icol );
+ return rc;
}
-inline const double& AlSymMat::AlSymMat_row_const::operator[]( size_t c ) const {
+inline AlMat operator*( const AlSymMat& lhs, const AlSymMat& rhs ) { return lhs.toMatrix() * rhs.toMatrix(); }
- if ( _r >= _a.Size || c >= _a.Size ) std::cerr << "Range error in AlSymMat::operator[][]" << std::endl;
+inline AlMat operator*( const AlMat& lhs, const AlSymMat& rhs ) { return lhs * rhs.toMatrix(); }
- if ( _r >= c ) {
- return *( _a.ptr_data + ( _r + 1 ) * _r / 2 + c );
- } else {
- return *( _a.ptr_data + ( c + 1 ) * c / 2 + _r );
- }
-}
-
-inline AlSymMat::AlSymMat_row::AlSymMat_row( AlSymMat& a, size_t r ) : _a( a ), _r( r ) {}
-
-inline AlSymMat::AlSymMat_row_const::AlSymMat_row_const( const AlSymMat& a, size_t r ) : _a( a ), _r( r ) {}
+inline AlMat operator*( const AlSymMat& lhs, const AlMat& rhs ) { return lhs.toMatrix() * rhs; }
inline AlSymMat operator*( double d, const AlSymMat& m ) {
AlSymMat rc = m;
rc *= d;
return rc;
}
-
-#endif // SOLVKERNEL_ALSYMMAT_H
diff --git a/Alignment/AlignKernel/AlignKernel/AlVec.h b/Alignment/AlignKernel/AlignKernel/AlVec.h
index 2a15ae601e5a88beacce415d1780f90e83e1c1ca..a87a912e638e08930793f7dea104c21aa687e022 100644
--- a/Alignment/AlignKernel/AlignKernel/AlVec.h
+++ b/Alignment/AlignKernel/AlignKernel/AlVec.h
@@ -1,5 +1,5 @@
/*****************************************************************************\
-* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration *
+* (c) Copyright 2000-2019 CERN for the benefit of the LHCb Collaboration *
* *
* This software is distributed under the terms of the GNU General Public *
* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
@@ -8,78 +8,21 @@
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
-#ifndef SOLVKERNEL_ALVEC_H
-#define SOLVKERNEL_ALVEC_H
-// from STD
-#include <iostream>
-#include <ostream>
-#include <vector>
+#pragma once
-class AlSymMat;
-class AlMat;
-
-class AlVec {
+#pragma GCC diagnostic push
+#if !defined( __clang__ )
+# pragma GCC diagnostic ignored "-Wdeprecated-copy"
+#endif
+#include <Eigen/Dense>
+#pragma GCC diagnostic pop
+/// Basically just a wrapper around the Eigen class, but with zero initialization
+class AlVec : public Eigen::VectorXd {
public:
- AlVec( size_t N );
- AlVec();
- AlVec( const AlVec& v );
- ~AlVec();
-
- double& operator[]( size_t i );
- const double& operator[]( size_t i ) const;
- double& operator()( size_t i ) { return *( ptr_data + i ); }
- const double& operator()( size_t i ) const { return *( ptr_data + i ); }
-
- AlVec& operator=( const AlVec& v );
- // AlVec& operator=(const double& v);
- AlVec operator+( const AlVec& v );
- AlVec operator+( const AlVec& v ) const;
- AlVec& operator+=( const AlVec& v );
- AlVec operator-( const AlVec& v );
- AlVec operator-( const AlVec& v ) const;
- AlVec& operator-=( const AlVec& v );
- double operator*( const AlVec& v );
- double operator*( const AlVec& v ) const;
- AlVec operator*( const AlMat& m );
- AlVec operator*( const AlMat& m ) const;
- AlVec operator*( const AlSymMat& m );
- AlVec operator*( const AlSymMat& m ) const;
- AlVec operator*( const double& d );
- AlVec operator*( const double& d ) const;
- AlVec& operator*=( const double& d );
-
- // insert block from SMatrix types
- template <typename T>
- void insert( unsigned row, const T& rhs ) {
- for ( unsigned i = 0; i < T::kSize; ++i ) ( *this )[row + i] = rhs( i );
- }
-
- void reSize( size_t Nnew );
- void removeElement( size_t Nr );
- void ShrinkDofs( std::vector<bool>& mask );
- size_t size() const { return Nele; }
-
-private:
- friend class AlSymMat;
- friend class AlMat;
-
- size_t Nele;
- double* ptr_data;
+ using Base = Eigen::VectorXd;
+ AlVec( Index n ) : Base( Base::Zero( n ) ) {}
+ AlVec( size_t n ) : Base( Base::Zero( n ) ) {}
+ using Base::Base;
};
-
-inline std::ostream& operator<<( std::ostream& lhs, const AlVec& rhs ) {
- lhs << "[ ";
- for ( size_t i = 0u; i < rhs.size(); ++i ) lhs << " " << rhs[i] << " ";
- lhs << " ]" << std::endl;
- return lhs;
-}
-
-inline AlVec operator*( double d, const AlVec& vec ) {
- AlVec rc = vec;
- rc *= d;
- return rc;
-}
-
-#endif // SOLVKERNEL_ALVEC_H
diff --git a/Alignment/AlignKernel/CMakeLists.txt b/Alignment/AlignKernel/CMakeLists.txt
index d5c21e4b1f7b2beef0ab888012ee3bd6376eedb4..2a5ba3621e1ff79fe76c8a8df8d24549b14cd4d7 100644
--- a/Alignment/AlignKernel/CMakeLists.txt
+++ b/Alignment/AlignKernel/CMakeLists.txt
@@ -14,42 +14,22 @@ gaudi_depends_on_subdirs(GaudiAlg
GaudiGSL
GaudiKernel)
-# helper function used as fall-back solution to find BLAS/LAPACK
-macro(simple_find_lib _name)
- find_library(${_name}_LIB ${_name})
- if(${_name}_LIB)
- message(STATUS "Found ${_name}: ${${_name}_LIB}")
- set(${_name}_LIBRARIES ${${_name}_LIB})
- set(${_name}_FOUND True CACHE INTERNAL "${_name} library found")
- else()
- message(FATAL_ERROR "Cannot find ${_name}")
- endif()
-endmacro()
-
find_package(ROOT)
find_package(GSL)
include_directories(SYSTEM ${ROOT_INCLUDE_DIRS})
-# Try regular LAPACK finder
-find_package(LAPACK QUIET)
-
-# Fall-back solution for HEPTools versions of LAPACK and BLAS
-if(NOT LAPACK_FOUND)
- simple_find_lib(BLAS)
- simple_find_lib(LAPACK)
- set(LAPACK_LIBRARIES ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES})
-endif()
-
+find_package(Eigen)
+include_directories(SYSTEM ${EIGEN_INCLUDE_DIRS})
gaudi_add_library(AlignKernel
src/*.cpp
PUBLIC_HEADERS AlignKernel
- INCLUDE_DIRS GSL
- LINK_LIBRARIES GSL LAPACK gfortran GaudiAlgLib GaudiGSLLib GaudiKernel)
+ INCLUDE_DIRS GSL Eigen
+ LINK_LIBRARIES GSL Eigen GaudiAlgLib GaudiGSLLib GaudiKernel)
gaudi_add_dictionary(AlignKernel
dict/AlignKernelDict.h
dict/AlignKernelDict.xml
- INCLUDE_DIRS GSL
- LINK_LIBRARIES GSL GaudiAlgLib GaudiGSLLib GaudiKernel AlignKernel)
+ INCLUDE_DIRS GSL Eigen
+ LINK_LIBRARIES GSL Eigen GaudiAlgLib GaudiGSLLib GaudiKernel AlignKernel)
diff --git a/Alignment/AlignKernel/src/AlMat.cpp b/Alignment/AlignKernel/src/AlMat.cpp
deleted file mode 100644
index 70c554f0ea01d1e650a551aa3d47ba797dc106d7..0000000000000000000000000000000000000000
--- a/Alignment/AlignKernel/src/AlMat.cpp
+++ /dev/null
@@ -1,405 +0,0 @@
-/*****************************************************************************\
-* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration *
-* *
-* This software is distributed under the terms of the GNU General Public *
-* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
-* *
-* In applying this licence, CERN does not waive the privileges and immunities *
-* granted to it by virtue of its status as an Intergovernmental Organization *
-* or submit itself to any jurisdiction. *
-\*****************************************************************************/
-#include "AlignKernel/AlMat.h"
-#include "AlignKernel/AlSymMat.h"
-#include "AlignKernel/AlVec.h"
-#include <iomanip>
-#include <iostream>
-
-AlMat::AlMat() {
-
- Ncol = 0;
- Nrow = 0;
- Nele = 0;
- ptr_data = 0; // set pointer to null
-}
-
-AlMat::AlMat( size_t N, size_t M ) {
-
- Nrow = N;
- Ncol = M;
- Nele = N * M;
- ptr_data = new double[Nele];
-
- for ( size_t i = 0; i < Nrow; i++ ) {
- for ( size_t j = 0; j < Ncol; j++ ) { *( ptr_data + i * Ncol + j ) = 0.; }
- }
-}
-
-AlMat::AlMat( const AlMat& m ) {
- Nrow = m.Nrow;
- Ncol = m.Ncol;
- Nele = Nrow * Ncol;
- ptr_data = new double[Nele];
- copy( m );
-}
-
-AlMat::~AlMat() { delete[] ptr_data; }
-
-void AlMat::copy( const AlMat& m ) {
-
- size_t n = Nrow <= m.Nrow ? Nrow : m.Nrow;
- size_t l = Ncol <= m.Ncol ? Ncol : m.Ncol;
-
- for ( size_t i = 0; i < n; i++ ) {
- for ( size_t j = 0; j < l; j++ ) { *( ptr_data + i * Ncol + j ) = *( m.ptr_data + i * m.Ncol + j ); }
- }
-}
-
-void AlMat::copyS( const AlSymMat& m ) {
-
- size_t n = Nrow <= m.Size ? Nrow : m.Size;
- size_t l = Ncol <= m.Size ? Ncol : m.Size;
-
- for ( size_t i = 0; i < n; i++ ) {
- for ( size_t j = 0; j < l; j++ ) {
-
- if ( j <= i ) {
- *( ptr_data + i * Ncol + j ) = *( m.ptr_data + ( i + 1 ) * i / 2 + j );
- } else {
- *( ptr_data + i * Ncol + j ) = *( m.ptr_data + ( j + 1 ) * j / 2 + i );
- }
- }
- }
-}
-
-AlMat& AlMat::operator=( const AlMat& m ) {
-
- 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 ( ptr_data != m.ptr_data ) {
- reSize( m.Nrow, m.Ncol );
- copy( m );
- };
- return ( *this );
-}
-
-AlMat& AlMat::operator=( const AlSymMat& m ) {
-
- if ( ( Nrow != 0 || Ncol != 0 ) && ( Nrow != m.Size || Ncol != m.Size ) ) {
- std::cerr << "AlMat AlSymMat Assignment: size do not match!" << std::endl;
- return *this;
- }
-
- if ( ptr_data != m.ptr_data ) {
- reSize( m.Size, m.Size );
- copyS( m );
- }
- return ( *this );
-}
-
-AlMat AlMat::operator+( const AlMat& m ) {
-
- if ( Nrow != m.Nrow || Ncol != m.Ncol ) {
- std::cerr << "operator+: size do not match!" << std::endl;
- return *this;
- }
-
- AlMat b( Nrow, Ncol );
-
- for ( size_t i = 0; i < Nrow; i++ ) {
- for ( size_t j = 0; j < Ncol; j++ ) { *( b.ptr_data + i * b.Ncol + j ) = *( ptr_data + i * Ncol + j ) + m[i][j]; }
- }
-
- return b;
-}
-
-AlMat AlMat::operator+( const AlMat& m ) const {
-
- if ( Nrow != m.Nrow || Ncol != m.Ncol ) {
- std::cerr << "operator+: size do not match!" << std::endl;
- return *this;
- }
-
- AlMat b( Nrow, Ncol );
-
- for ( size_t i = 0; i < Nrow; i++ ) {
- for ( size_t j = 0; j < Ncol; j++ ) { *( b.ptr_data + i * b.Ncol + j ) = *( ptr_data + i * Ncol + j ) + m[i][j]; }
- }
-
- return b;
-}
-
-AlMat& AlMat::operator+=( const AlMat& m ) {
-
- if ( Nrow != m.Nrow || Ncol != m.Ncol ) {
- std::cerr << "operator+=: size do not match!" << std::endl;
- return *this;
- }
-
- for ( size_t i = 0; i < Nrow; i++ ) {
- for ( size_t j = 0; j < Ncol; j++ ) { *( ptr_data + i * Ncol + j ) += *( m.ptr_data + i * m.Ncol + j ); };
- };
-
- return *this;
-}
-
-AlMat AlMat::operator-( const AlMat& m ) {
- if ( Nrow != m.Nrow || Ncol != m.Ncol ) {
- std::cerr << "operator-: size do not match!" << std::endl;
- return *this;
- }
-
- AlMat b( Nrow, Ncol );
-
- for ( size_t i = 0; i < Nrow; i++ ) {
- for ( size_t j = 0; j < Ncol; j++ ) { *( b.ptr_data + i * b.Ncol + j ) = *( ptr_data + i * Ncol + j ) - m[i][j]; }
- }
-
- return b;
-}
-
-AlMat AlMat::operator-( const AlMat& m ) const {
- if ( Nrow != m.Nrow || Ncol != m.Ncol ) {
- std::cerr << "operator-: size do not match!" << std::endl;
- return *this;
- }
-
- AlMat b( Nrow, Ncol );
-
- for ( size_t i = 0; i < Nrow; i++ ) {
- for ( size_t j = 0; j < Ncol; j++ ) { *( b.ptr_data + i * b.Ncol + j ) = *( ptr_data + i * Ncol + j ) - m[i][j]; }
- }
-
- return b;
-}
-
-AlMat& AlMat::operator-=( const AlMat& m ) {
-
- if ( Nrow != m.Nrow || Ncol != m.Ncol ) {
- std::cerr << "operator-=: size do not match!" << std::endl;
- return *this;
- }
-
- for ( size_t i = 0; i < Nrow; i++ ) {
- for ( size_t j = 0; j < Ncol; j++ ) { *( ptr_data + i * Ncol + j ) -= *( m.ptr_data + i * m.Ncol + j ); };
- };
-
- return *this;
-}
-
-AlMat AlMat::operator*( const AlMat& m ) {
-
- if ( Ncol != m.Nrow ) {
- std::cerr << "operator*: size do not match!" << std::endl;
- return m;
- }
-
- AlMat b( Nrow, m.Ncol );
-
- for ( size_t i = 0; i < Nrow; i++ ) {
- for ( size_t j = 0; j < m.Ncol; j++ ) {
- for ( size_t k = 0; k < m.Nrow; k++ ) b[i][j] += *( ptr_data + i * Ncol + k ) * ( m[k][j] );
- }
- }
-
- return b;
-}
-
-AlMat AlMat::operator*( const AlMat& m ) const {
-
- if ( Ncol != m.Nrow ) {
- std::cerr << "operator*: size do not match!" << std::endl;
- return m;
- }
-
- AlMat b( Nrow, m.Ncol );
-
- for ( size_t i = 0; i < Nrow; i++ ) {
- for ( size_t j = 0; j < m.Ncol; j++ ) {
- for ( size_t k = 0; k < m.Nrow; k++ ) b[i][j] += *( ptr_data + i * Ncol + k ) * ( m[k][j] );
- }
- }
-
- return b;
-}
-
-AlMat AlMat::operator*( const AlSymMat& m ) {
- if ( Ncol != m.Size ) {
- std::cerr << "operator*: size do not match!" << std::endl;
- return *this;
- }
-
- AlMat b( Nrow, m.Size );
-
- for ( size_t i = 0; i < Nrow; i++ ) {
- for ( size_t j = 0; j < m.Size; j++ ) {
- for ( size_t k = 0; k < m.Size; k++ ) *( b.ptr_data + i * b.Ncol + j ) += *( ptr_data + i * Ncol + k ) * m[k][j];
- }
- }
- return b;
-}
-
-AlMat AlMat::operator*( const AlSymMat& m ) const {
- if ( Ncol != m.Size ) {
- std::cerr << "operator*: size do not match!" << std::endl;
- return *this;
- }
-
- AlMat b( Nrow, m.Size );
-
- for ( size_t i = 0; i < Nrow; i++ ) {
- for ( size_t j = 0; j < m.Size; j++ ) {
- for ( size_t k = 0; k < m.Size; k++ ) *( b.ptr_data + i * b.Ncol + j ) += *( ptr_data + i * Ncol + k ) * m[k][j];
- }
- }
- return b;
-}
-
-AlVec AlMat::operator*( const AlVec& v ) {
-
- if ( Ncol != v.Nele ) {
- std::cerr << "operator*: size do not match! " << std::endl;
- return v;
- }
-
- AlVec b( Nrow );
-
- for ( size_t i = 0; i < Nrow; i++ ) {
-
- for ( size_t j = 0; j < Ncol; j++ ) b[i] += ( *( ptr_data + i * Ncol + j ) ) * ( v[j] );
- }
-
- return b;
-}
-
-AlVec AlMat::operator*( const AlVec& v ) const {
-
- if ( Ncol != v.Nele ) {
- std::cerr << "operator*: size do not match! " << std::endl;
- return v;
- }
-
- AlVec b( Nrow );
-
- for ( size_t i = 0; i < Nrow; i++ ) {
-
- for ( size_t j = 0; j < Ncol; j++ ) b[i] += ( *( ptr_data + i * Ncol + j ) ) * ( v[j] );
- }
-
- return b;
-}
-
-AlMat& AlMat::operator*=( const double& d ) {
-
- double* p = ptr_data + Nele;
-
- while ( p > ptr_data ) *( --p ) *= d;
-
- return *this;
-}
-
-// transposition
-
-AlMat AlMat::T() const {
- AlMat b( Ncol, Nrow );
- for ( size_t i = 0; i < b.Nrow; i++ ) {
- for ( size_t j = 0; j < b.Ncol; j++ ) { b[i][j] = *( ptr_data + j * Ncol + i ); }
- }
- return b;
-}
-
-// 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;
- }
-
- AlSymMat b( Nrow );
- // double minelem=99999999.;
- // double maxelem=0.;
-
- for ( size_t i = 0; i < b.Size; i++ ) {
- for ( size_t j = 0; j <= i; j++ ) {
-
- b[i][j] = ( *this )[i][j];
- // if (fabs(b[i][j])<fabs(minelem)) minelem=b[i][j];
- // if (fabs(b[i][j])>fabs(maxelem)) maxelem=b[i][j];
- }
- }
-
- // proceed to regularization before and after inversion, if needed
-
- // std::cout<<" >>>>> MAXELEM: "<<maxelem<<std::endl;
- // Norm = 1./maxelem;
- // std::cout<<" >>>>> NORM: "<<Norm<<std::endl;
- /* if ((log(fabs(maxelem)))<0.) {
- size_t factor = size_t(log(fabs(maxelem)))-1;
- Norm=1./1e(factor);
- }*/
- b *= Norm;
-
- b.invert( ierr );
- b *= Norm;
-
- if ( ierr == 0 ) ( *this ).copyS( b );
-}
-
-// reSize
-
-void AlMat::reSize( size_t Nnew, size_t Mnew ) {
-
- if ( Nnew != Nrow || Mnew != Ncol ) {
-
- double* p = ptr_data;
- Nele = Nnew * Mnew;
- ptr_data = new double[Nele];
- size_t Nrow_old = Nrow;
- size_t Ncol_old = Ncol;
-
- Nrow = Nnew;
- Ncol = Mnew;
- size_t k = Nrow <= Nrow_old ? Nrow : Nrow_old;
- size_t l = Ncol <= Ncol_old ? Ncol : Ncol_old;
-
- for ( size_t i = 0; i < k; i++ ) {
- for ( size_t j = 0; j < l; j++ ) { *( ptr_data + i * Ncol + j ) = *( p + i * Ncol_old + j ); }
- }
-
- delete[] p;
- }
-}
-
-void AlMat::removeRow( size_t Nr ) {
-
- // reorganize matrix by filling holes, etc...
- for ( size_t i = Nr; i < Nrow - 1; i++ ) {
- for ( size_t j = 0; j < Ncol; j++ ) { ( *this )[i][j] = ( *this )[i + 1][j]; }
- }
- // delete obsolete space and shrink matrix row size by one unit
- size_t max = Nrow - 1;
- for ( size_t col = 0; col < Ncol; col++ ) delete ( ptr_data + max * Ncol + col );
-
- Nrow--;
-}
-
-void AlMat::removeCol( size_t Nr ) {
-
- // reorganize matrix by filling holes, etc...
- for ( size_t i = 0; i < Nrow; i++ ) {
- for ( size_t j = Nr; j < Ncol - 1; j++ ) { ( *this )[i][j] = ( *this )[i][j + 1]; }
- }
- // delete obsolete space and shrink matrix column size by one unit
- size_t max = Ncol - 1;
- for ( size_t row = 0; row < Nrow; row++ ) delete ( ptr_data + row * Ncol + max );
-
- Ncol--;
-}
-
-size_t AlMat::nrow() const { return Nrow; }
-
-size_t AlMat::ncol() const { return Ncol; }
diff --git a/Alignment/AlignKernel/src/AlSymMat.cpp b/Alignment/AlignKernel/src/AlSymMat.cpp
index 6abfa51b06a8c969d723f254a63e4f1d43d1a3de..8e60c2393999c438a688038425dc5b9e63ace918 100644
--- a/Alignment/AlignKernel/src/AlSymMat.cpp
+++ b/Alignment/AlignKernel/src/AlSymMat.cpp
@@ -11,389 +11,69 @@
#include "gsl/gsl_eigen.h"
#include "gsl/gsl_math.h"
#include "gsl/gsl_matrix.h"
-#include <float.h> //for DBL_EPSILON
-#include <iomanip>
-#include <iostream>
-#include <math.h>
-#include <string.h>
-//#include <Nag/nag.h>
-//#include <Nag/nagf03.h>
-#include "AlignKernel/AlMat.h"
-#include "AlignKernel/AlSymMat.h"
-#include "AlignKernel/AlVec.h"
-
-extern "C" {
-void dsptrf_( char*, const int*, double*, int*, int* );
-void dsptri_( char*, const int*, double*, int*, double*, int* );
-void dspev_( char*, char*, const int*, double*, double*, double*, const int*, double*, int* );
-}
-
-AlSymMat::AlSymMat() {
- Size = 0;
- Nele = 0;
- ptr_data = 0; // set pointer to null
-}
-
-AlSymMat::AlSymMat( size_t N ) {
-
- Size = N;
- Nele = N * ( N + 1 ) / 2;
- ptr_data = new double[Nele];
-
- for ( size_t i = 0; i < Size; i++ ) {
- for ( size_t j = 0; j <= i; j++ ) { *( ptr_data + ( i + 1 ) * i / 2 + j ) = 0.; }
- }
-}
-
-AlSymMat::AlSymMat( const AlSymMat& m ) {
- Size = m.Size;
- Nele = m.Nele;
- ptr_data = new double[Nele];
- copy( m );
-}
-
-AlSymMat::~AlSymMat() { delete[] ptr_data; }
-
-void AlSymMat::copy( const AlSymMat& m ) {
-
- size_t n = Size <= m.Size ? Size : m.Size;
-
- for ( size_t i = 0; i < n; i++ ) {
- for ( size_t j = 0; j <= i; j++ ) { *( ptr_data + ( i + 1 ) * i / 2 + j ) = m[i][j]; }
- }
-}
-
-AlSymMat& AlSymMat::operator=( const AlSymMat& m ) {
- if ( size() != m.size() ) {
- delete[] ptr_data;
- Nele = m.Nele;
- Size = m.Size;
- ptr_data = new double[Nele];
- }
- memcpy( ptr_data, m.ptr_data, Nele * sizeof( double ) );
- return *this;
-}
-
-AlSymMat AlSymMat::operator+( const AlSymMat& m ) {
-
- if ( Size != m.Size ) {
- std::cerr << "operator+: size do not match!" << std::endl;
- return *this;
- }
-
- AlSymMat b( Size );
-
- for ( size_t i = 0; i < Size; i++ ) {
- for ( size_t j = 0; j <= i; j++ ) { b[i][j] = *( ptr_data + ( i + 1 ) * i / 2 + j ) + m[i][j]; }
- }
-
- return b;
-}
-
-AlSymMat AlSymMat::operator+( const AlSymMat& m ) const {
-
- if ( Size != m.Size ) {
- std::cerr << "operator+: size do not match!" << std::endl;
- return *this;
- }
-
- AlSymMat b( Size );
-
- for ( size_t i = 0; i < Size; i++ ) {
- for ( size_t j = 0; j <= i; j++ ) { b[i][j] = *( ptr_data + ( i + 1 ) * i / 2 + j ) + m[i][j]; }
- }
+#include "AlignKernel/AlSymMat.h"
+AlSymMat AlSymMat::operator+( const AlSymMat& rhs ) const {
+ auto b = *this;
+ b += rhs;
return b;
}
-AlSymMat& AlSymMat::operator+=( const AlSymMat& m ) {
+AlSymMat& AlSymMat::operator+=( const AlSymMat& rhs ) {
- if ( Size != m.Size ) {
+ if ( size() != rhs.size() ) {
std::cerr << "operator+=: size do not match!" << std::endl;
return *this;
}
- for ( size_t i = 0; i < Size; i++ ) {
- for ( size_t j = 0; j <= i; j++ ) { *( ptr_data + ( i + 1 ) * i / 2 + j ) += m[i][j]; };
- };
-
+ for ( size_t i = 0; i < m_data.size(); ++i ) m_data[i] += rhs.m_data[i];
return *this;
}
-AlSymMat AlSymMat::operator-( const AlSymMat& m ) {
-
- if ( Size != m.Size ) {
- std::cerr << "operator-: size do not match!" << std::endl;
- return *this;
- }
-
- AlSymMat b( Size );
-
- for ( size_t i = 0; i < Size; i++ ) {
- for ( size_t j = 0; j <= i; j++ ) { b[i][j] = *( ptr_data + ( i + 1 ) * i / 2 + j ) - m[i][j]; }
- }
-
+AlSymMat AlSymMat::operator-( const AlSymMat& rhs ) const {
+ auto b = *this;
+ b -= rhs;
return b;
}
-AlSymMat AlSymMat::operator-( const AlSymMat& m ) const {
+AlSymMat& AlSymMat::operator-=( const AlSymMat& rhs ) {
- if ( Size != m.Size ) {
- std::cerr << "operator-: size do not match!" << std::endl;
+ if ( size() != rhs.size() ) {
+ std::cerr << "operator+=: size do not match!" << std::endl;
return *this;
}
- AlSymMat b( Size );
-
- for ( size_t i = 0; i < Size; i++ ) {
- for ( size_t j = 0; j <= i; j++ ) { b[i][j] = *( ptr_data + ( i + 1 ) * i / 2 + j ) - m[i][j]; }
- }
-
- return b;
-}
-
-AlMat AlSymMat::operator*( const AlSymMat& m ) {
-
- if ( Size != m.Size ) {
- std::cerr << "operator*: size do not match!" << std::endl;
-
- AlMat b( Size, m.Size );
-
- return b;
- }
-
- AlMat b( Size, Size );
-
- for ( size_t i = 0; i < Size; i++ ) {
- for ( size_t j = 0; j < Size; j++ ) {
- for ( size_t k = 0; k < Size; k++ ) {
- if ( k <= i ) b[i][j] += *( ptr_data + ( i + 1 ) * i / 2 + k ) * ( m[k][j] );
- if ( k > i ) b[i][j] += *( ptr_data + ( k + 1 ) * k / 2 + i ) * ( m[k][j] );
- }
- }
- }
-
- return b;
-}
-
-AlMat AlSymMat::operator*( const AlSymMat& m ) const {
-
- if ( Size != m.Size ) {
- std::cerr << "operator*: size do not match!" << std::endl;
-
- AlMat b( Size, m.Size );
-
- return b;
- }
-
- AlMat b( Size, Size );
-
- for ( size_t i = 0; i < Size; i++ ) {
- for ( size_t j = 0; j < Size; j++ ) {
- for ( size_t k = 0; k < Size; k++ ) {
- if ( k <= i ) b[i][j] += *( ptr_data + ( i + 1 ) * i / 2 + k ) * ( m[k][j] );
- if ( k > i ) b[i][j] += *( ptr_data + ( k + 1 ) * k / 2 + i ) * ( m[k][j] );
- }
- }
- }
-
- return b;
-}
-
-AlMat AlSymMat::operator*( const AlMat& m ) {
-
- if ( Size != m.Nrow ) {
- std::cerr << "operator*: size do not match!" << std::endl;
- return m;
- }
-
- AlMat b( Size, m.Ncol );
-
- for ( size_t i = 0; i < Size; i++ ) {
- for ( size_t j = 0; j < m.Ncol; j++ ) {
- for ( size_t k = 0; k < Size; k++ ) {
- if ( k <= i ) b[i][j] += *( ptr_data + ( i + 1 ) * i / 2 + k ) * ( m[k][j] );
- if ( k > i ) b[i][j] += *( ptr_data + ( k + 1 ) * k / 2 + i ) * ( m[k][j] );
- }
- }
- }
-
- return b;
-}
-
-AlMat AlSymMat::operator*( const AlMat& m ) const {
-
- if ( Size != m.Nrow ) {
- std::cerr << "operator*: size do not match!" << std::endl;
- return m;
- }
-
- AlMat b( Size, m.Ncol );
-
- for ( size_t i = 0; i < Size; i++ ) {
- for ( size_t j = 0; j < m.Ncol; j++ ) {
- for ( size_t k = 0; k < Size; k++ ) {
- if ( k <= i ) b[i][j] += *( ptr_data + ( i + 1 ) * i / 2 + k ) * ( m[k][j] );
- if ( k > i ) b[i][j] += *( ptr_data + ( k + 1 ) * k / 2 + i ) * ( m[k][j] );
- }
- }
- }
-
- return b;
-}
-
-AlVec AlSymMat::operator*( const AlVec& v ) {
-
- if ( Size != v.Nele ) {
- std::cerr << "operator*: size do not match! " << std::endl;
- return v;
- }
-
- AlVec b( Size );
-
- for ( size_t i = 0; i < Size; i++ ) {
- for ( size_t j = 0; j < Size; j++ ) {
- if ( j <= i ) b[i] += *( ptr_data + ( i + 1 ) * i / 2 + j ) * v[j];
- if ( j > i ) b[i] += *( ptr_data + ( j + 1 ) * j / 2 + i ) * v[j];
- }
- }
-
- return b;
-}
-
-AlVec AlSymMat::operator*( const AlVec& v ) const {
-
- if ( Size != v.Nele ) {
- std::cerr << "operator*: size do not match! " << std::endl;
- return v;
- }
-
- AlVec b( Size );
-
- for ( size_t i = 0; i < Size; i++ ) {
- for ( size_t j = 0; j < Size; j++ ) {
- if ( j <= i ) b[i] += *( ptr_data + ( i + 1 ) * i / 2 + j ) * v[j];
- if ( j > i ) b[i] += *( ptr_data + ( j + 1 ) * j / 2 + i ) * v[j];
- }
- }
-
- return b;
+ for ( size_t i = 0; i < m_data.size(); ++i ) m_data[i] -= rhs.m_data[i];
+ return *this;
}
AlSymMat& AlSymMat::operator*=( const double& d ) {
-
- double* p = ptr_data + Nele;
-
- while ( p > ptr_data ) *( --p ) *= d;
-
+ for ( auto& e : m_data ) e *= d;
return *this;
}
-void AlSymMat::reSize( size_t Nnew ) {
-
- if ( Nnew != Size && Size != 0 ) {
-
- double* p = ptr_data;
- Nele = Nnew * ( Nnew + 1 ) / 2;
- ptr_data = new double[Nele];
- size_t Size_old = Size;
- Size = Nnew;
- size_t k = Size <= Size_old ? Size : Size_old;
- for ( size_t i = 0; i < k; i++ ) {
- for ( size_t j = 0; j <= i; j++ ) { *( ptr_data + ( i + 1 ) * i / 2 + j ) = *( p + ( i + 1 ) * i / 2 + j ); }
- }
-
- delete[] p;
- }
-
- if ( Size == 0 ) *this = AlSymMat( Nnew );
-}
-
-void AlSymMat::removeElement( size_t Nr ) {
-
- // reorganize matrix by filling holes, etc...
- for ( size_t i = Nr; i < Size - 1; i++ ) {
- for ( size_t j = 0; j <= i; j++ ) {
- if ( j < Nr )
- ( *this )[i][j] = ( *this )[i + 1][j];
- else
- ( *this )[i][j] = ( *this )[i + 1][j + 1];
- }
- }
- // delete obsolete space and shrink matrix size by one unit
- size_t max = Size - 1;
- for ( size_t col = 0; col < Size; col++ ) delete ( ptr_data + ( max + 1 ) * max / 2 + col );
- Size--;
-}
-
-void AlSymMat::ShrinkDofs( std::vector<bool>& mask ) {
-
- if ( mask.size() != Size ) {
- std::cerr << "Beware, mask size different from matrix size, skip..." << std::endl;
- return;
- }
-
- size_t offset = 0;
-
- for ( size_t iR = 0; iR < mask.size(); iR++ ) {
- if ( !mask[iR] ) {
- removeElement( iR - offset );
- offset++;
- }
- }
-}
-
-void AlSymMat::invert( int& ierr ) {
-
- // add a protection to detect singular matrices:
-
- double det = ( *this ).determinant();
- double detcut = 1.E-20;
-
- if ( fabs( det ) > detcut ) {
-
- char uplo = 'U';
- int N = Size;
- int* ipiv = new int[N];
- double* work = new double[N];
- double* mx = ptr_data;
-
- dsptrf_( &uplo, &N, mx, ipiv, &ierr );
-
- dsptri_( &uplo, &N, mx, ipiv, work, &ierr );
-
- delete[] ipiv;
- delete[] work;
-
- } else {
-
- ierr = 999;
- }
-}
-
-double AlSymMat::determinant() {
+double AlSymMat::determinant() const {
double deter = 1.;
// get determinant using LU decomposition + Crout algorithm
- AlMat A( Size, Size );
- A.copyS( *this );
+ AlMat A = this->toMatrix();
- double AMAX, DUM, SUM;
- int CODE;
- double TINY = 1.E-20;
- double D;
- size_t N = Size;
+ double AMAX, DUM, SUM;
+ int CODE;
+ double TINY = 1.E-20;
+ double D;
+ const auto N = rows();
D = 1.;
CODE = 0;
- for ( size_t I = 0; I < N; I++ ) {
+ for ( size_t I = 0; I < N; ++I ) {
AMAX = 0.;
- for ( size_t J = 0; J < N; J++ ) {
- if ( fabs( A[I][J] ) > AMAX ) AMAX = fabs( A[I][J] );
+ for ( size_t J = 0; J < N; ++J ) {
+ if ( std::abs( A( I, J ) ) > AMAX ) AMAX = std::abs( A( I, J ) );
}
if ( AMAX < TINY ) {
CODE = 1; // matrix is singular
@@ -405,38 +85,38 @@ double AlSymMat::determinant() {
for ( size_t I = 0; I < N; I++ ) {
- SUM = A[I][J];
+ SUM = A( I, J );
size_t KMAX = ( I < J ) ? I : J;
- for ( size_t K = 0; K < KMAX; K++ ) SUM = SUM - A[I][K] * A[K][J];
- A[I][J] = SUM;
+ for ( size_t K = 0; K < KMAX; K++ ) SUM = SUM - A( I, K ) * A( K, J );
+ A( I, J ) = SUM;
}
// find pivot and exchange if necessary
size_t IMAX = J;
for ( size_t I = J + 1; I < N; I++ ) {
- if ( fabs( A[I][J] ) > fabs( A[IMAX][J] ) ) IMAX = I;
+ if ( std::abs( A( I, J ) ) > std::abs( A( IMAX, J ) ) ) IMAX = I;
}
if ( IMAX != J ) {
for ( size_t K = 0; K < N; K++ ) {
- DUM = A[IMAX][K];
- A[IMAX][K] = A[J][K];
- A[J][K] = DUM;
+ DUM = A( IMAX, K );
+ A( IMAX, K ) = A( J, K );
+ A( J, K ) = DUM;
}
D = -D;
}
- if ( J < N && ( fabs( A[J][J] ) > TINY ) ) {
+ if ( J < N && ( std::abs( A( J, J ) ) > TINY ) ) {
- DUM = 1. / A[J][J];
- for ( size_t I = J + 1; I < N; I++ ) A[I][J] = A[I][J] * DUM;
+ DUM = 1. / A( J, J );
+ for ( size_t I = J + 1; I < N; I++ ) A( I, J ) = A( I, J ) * DUM;
}
}
if ( CODE == 0 ) {
deter = deter * D;
- for ( size_t I = 0; I < N; I++ ) { deter = deter * A[I][I]; }
+ for ( size_t I = 0; I < N; I++ ) { deter = deter * A( I, I ); }
return deter;
} else {
@@ -444,73 +124,47 @@ double AlSymMat::determinant() {
}
}
-int AlSymMat::diagonalize( char jobz, AlVec& w, AlMat& z ) {
-
- int info = 0;
- char uplo = 'U';
- int N = Size;
- double* work = new double[3 * N];
+int AlSymMat::diagonalize_GSL( AlVec& egval, AlMat& egvec ) const {
- double* ap = ptr_data;
+ const auto N = rows();
- dspev_( &jobz, &uplo, &N, ap, w.ptr_data, z.ptr_data, &N, work, &info );
-
- delete[] work;
-
- return info;
-}
-
-void AlSymMat::diagonalize_GSL( AlVec& egval, AlMat& egvec ) {
-
- size_t N = Size;
-
- std::cout << " Debug...before gsl_matrix_alloc" << std::endl;
+ // std::cout << " Debug...before gsl_matrix_alloc" << std::endl;
gsl_matrix* m = gsl_matrix_alloc( N, N );
- std::cout << " Debug...after gsl_matrix_alloc" << std::endl;
+ // std::cout << " Debug...after gsl_matrix_alloc" << std::endl;
// if (1) {
// AlMat V(N,N);
// V.copyS((*this));
// m->data = V.ptr_data;
for ( size_t ig = 0; ig < N; ig++ ) {
- for ( size_t jg = 0; jg < N; jg++ ) { gsl_matrix_set( m, ig, jg, ( *this )[ig][jg] ); }
+ for ( size_t jg = 0; jg < N; jg++ ) { gsl_matrix_set( m, ig, jg, ( *this )( ig, jg ) ); }
}
// gsl_matrix_set(m,0,0,temp00);
// std::cout<<"gsl00: "<<gsl_matrix_get(m,0,0)<<std::endl;
// }
- std::cout << "Debut....printing gsl matrix elements" << std::endl;
-
- for ( size_t ig = 0; ig < N; ig++ ) {
- for ( size_t jg = 0; jg < N; jg++ ) { std::cout << gsl_matrix_get( m, ig, jg ) << " "; }
- std::cout << std::endl;
- }
-
- std::cout << " Debug...after copy ptr data" << std::endl;
-
gsl_vector* eval = gsl_vector_alloc( N );
gsl_matrix* evec = gsl_matrix_alloc( N, N );
- gsl_eigen_symmv_workspace* w = gsl_eigen_symmv_alloc( N );
- gsl_eigen_symmv( m, eval, evec, w );
+ gsl_eigen_symmv_workspace* w = gsl_eigen_symmv_alloc( N );
+ int rc = gsl_eigen_symmv( m, eval, evec, w );
gsl_eigen_symmv_free( w );
- std::cout << " Debug...after diagonalization" << std::endl;
+ // std::cout << " Debug...after diagonalization" << std::endl;
// sort eigen values/eigen vectors in ascending order
gsl_eigen_symmv_sort( eval, evec, GSL_EIGEN_SORT_VAL_ASC );
- egval.ptr_data = eval->data;
- egvec.ptr_data = evec->data;
+ for ( size_t ig = 0; ig < N; ++ig ) {
+ egval( ig ) = gsl_vector_get( eval, ig );
+ for ( size_t jg = 0; jg < N; ++jg ) { egvec( ig, jg ) = gsl_matrix_get( evec, ig, jg ); }
+ }
+ gsl_matrix_free( evec );
+ gsl_vector_free( eval );
- // for(int ig=0;ig<N;ig++) {
- // egval[ig]=gsl_vector_get(eval,ig);
- // for(int jg=0;jg<N;jg++) {
- // egvec[ig][jg]=gsl_matrix_get(evec,ig,jg);
- // }
- //}
+ return rc;
}
diff --git a/Alignment/AlignKernel/src/AlVec.cpp b/Alignment/AlignKernel/src/AlVec.cpp
deleted file mode 100644
index c1eebc41b15eac7bdf1d4fc8f72551895b239a46..0000000000000000000000000000000000000000
--- a/Alignment/AlignKernel/src/AlVec.cpp
+++ /dev/null
@@ -1,289 +0,0 @@
-/*****************************************************************************\
-* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration *
-* *
-* This software is distributed under the terms of the GNU General Public *
-* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
-* *
-* In applying this licence, CERN does not waive the privileges and immunities *
-* granted to it by virtue of its status as an Intergovernmental Organization *
-* or submit itself to any jurisdiction. *
-\*****************************************************************************/
-#include "AlignKernel/AlVec.h"
-#include "AlignKernel/AlMat.h"
-#include "AlignKernel/AlSymMat.h"
-#include <iomanip>
-#include <iostream>
-#include <string.h>
-
-AlVec::AlVec() {
-
- Nele = 0;
- ptr_data = 0; // set pointer to null
-}
-
-AlVec::AlVec( size_t N ) {
-
- Nele = N;
- ptr_data = new double[Nele];
- for ( size_t i = 0; i < Nele; i++ ) *( ptr_data + i ) = 0.;
-}
-
-AlVec::AlVec( const AlVec& v ) {
-
- Nele = v.Nele;
- ptr_data = new double[Nele];
- memcpy( ptr_data, v.ptr_data, Nele * sizeof( double ) );
-}
-
-AlVec::~AlVec() { delete[] ptr_data; }
-
-double& AlVec::operator[]( size_t i ) {
-
- if ( i >= Nele ) {
- std::cerr << "AlVec: Index too large! " << std::endl;
- return ptr_data[0];
- };
-
- return *( ptr_data + i );
-}
-
-const double& AlVec::operator[]( size_t i ) const {
-
- if ( i >= Nele ) {
- std::cerr << "AlVec: Index too large! " << std::endl;
- return ptr_data[0];
- };
-
- return *( ptr_data + i );
-}
-
-AlVec& AlVec::operator=( const AlVec& v ) {
- if ( size() != v.size() ) {
- delete[] ptr_data;
- Nele = v.Nele;
- ptr_data = new double[Nele];
- }
- memcpy( ptr_data, v.ptr_data, Nele * sizeof( double ) );
- return *this;
-}
-
-AlVec AlVec::operator+( const AlVec& v ) {
-
- if ( Nele != v.Nele ) {
- std::cerr << "operator+: vectors size does not match!" << std::endl;
- return *this;
- };
- AlVec b( Nele );
-
- for ( size_t i = 0; i < Nele; i++ ) b[i] = *( ptr_data + i ) + v[i];
- return b;
-}
-
-AlVec AlVec::operator+( const AlVec& v ) const {
-
- if ( Nele != v.Nele ) {
- std::cerr << "operator+: vectors size does not match!" << std::endl;
- return *this;
- };
- AlVec b( Nele );
-
- for ( size_t i = 0; i < Nele; i++ ) b[i] = *( ptr_data + i ) + v[i];
- return b;
-}
-
-AlVec& AlVec::operator+=( const AlVec& v ) {
-
- if ( Nele != v.Nele ) {
- std::cerr << "operator+=: vectors size does not match!" << std::endl;
- return *this;
- };
-
- for ( size_t i = 0; i < Nele; i++ ) *( ptr_data + i ) += v[i];
- return *this;
-}
-
-AlVec AlVec::operator-( const AlVec& v ) {
-
- if ( Nele != v.Nele ) {
- std::cerr << "operator+: vectors size does not match!" << std::endl;
- return *this;
- }
- AlVec b( Nele );
-
- for ( size_t i = 0; i < Nele; i++ ) b[i] = *( ptr_data + i ) - v[i];
- return b;
-}
-
-AlVec AlVec::operator-( const AlVec& v ) const {
-
- if ( Nele != v.Nele ) {
- std::cerr << "operator+: vectors size does not match!" << std::endl;
- return *this;
- }
- AlVec b( Nele );
-
- for ( size_t i = 0; i < Nele; i++ ) b[i] = *( ptr_data + i ) - v[i];
- return b;
-}
-
-AlVec& AlVec::operator-=( const AlVec& v ) {
-
- if ( Nele != v.Nele ) {
- std::cerr << "operator+=: vectors size does not match!" << std::endl;
- return *this;
- }
-
- for ( size_t i = 0; i < Nele; i++ ) *( ptr_data + i ) -= v[i];
- return *this;
-}
-
-double AlVec::operator*( const AlVec& v ) {
- double b = 0.;
- if ( Nele != v.Nele ) {
- std::cerr << "scalar product: vectors size does not match!" << std::endl;
- return b;
- }
-
- for ( size_t i = 0; i < Nele; i++ ) b += *( ptr_data + i ) * v[i];
- return b;
-}
-
-double AlVec::operator*( const AlVec& v ) const {
- double b = 0.;
- if ( Nele != v.Nele ) {
- std::cerr << "scalar product: vectors size does not match!" << std::endl;
- return b;
- }
-
- for ( size_t i = 0; i < Nele; i++ ) b += *( ptr_data + i ) * v[i];
- return b;
-}
-
-AlVec AlVec::operator*( const AlMat& m ) {
-
- if ( Nele != m.Nrow ) {
- std::cerr << "Left hand vector-matrix multiplication: size does not match!" << std::endl;
- return *this;
- }
-
- AlVec b( m.Ncol );
-
- for ( size_t i = 0; i < m.Ncol; i++ ) {
- for ( size_t j = 0; j < Nele; j++ ) b[i] += *( ptr_data + j ) * m[j][i];
- }
-
- return b;
-}
-
-AlVec AlVec::operator*( const AlMat& m ) const {
-
- if ( Nele != m.Nrow ) {
- std::cerr << "Left hand vector-matrix multiplication: size does not match!" << std::endl;
- return *this;
- }
-
- AlVec b( m.Ncol );
-
- for ( size_t i = 0; i < m.Ncol; i++ ) {
- for ( size_t j = 0; j < Nele; j++ ) b[i] += *( ptr_data + j ) * m[j][i];
- }
-
- return b;
-}
-
-AlVec AlVec::operator*( const AlSymMat& m ) {
-
- if ( Nele != m.Size ) {
- std::cerr << "Left hand vector-matrix multiplication: size does not match!" << std::endl;
- return *this;
- }
-
- AlVec b( m.Size );
-
- for ( size_t i = 0; i < m.Size; i++ ) {
- for ( size_t j = 0; j < Nele; j++ ) b[i] += *( ptr_data + j ) * m[j][i];
- }
-
- return b;
-}
-
-AlVec AlVec::operator*( const AlSymMat& m ) const {
-
- if ( Nele != m.Size ) {
- std::cerr << "Left hand vector-matrix multiplication: size does not match!" << std::endl;
- return *this;
- }
-
- AlVec b( m.Size );
-
- for ( size_t i = 0; i < m.Size; i++ ) {
- for ( size_t j = 0; j < Nele; j++ ) b[i] += *( ptr_data + j ) * m[j][i];
- }
-
- return b;
-}
-
-AlVec& AlVec::operator*=( const double& d ) {
-
- for ( size_t i = 0; i < Nele; i++ ) *( ptr_data + i ) *= d;
-
- return *this;
-}
-
-AlVec AlVec::operator*( const double& d ) {
- AlVec b( Nele );
- for ( size_t i = 0; i < Nele; i++ ) b[i] = *( ptr_data + i ) * d;
- return b;
-}
-
-AlVec AlVec::operator*( const double& d ) const {
- AlVec b( Nele );
- for ( size_t i = 0; i < Nele; i++ ) b[i] = *( ptr_data + i ) * d;
- return b;
-}
-
-void AlVec::reSize( size_t Nnew ) {
-
- if ( Nnew != Nele && Nele != 0 ) {
- double* p = ptr_data;
- size_t Nele_old = Nele;
- ptr_data = new double[Nnew];
- Nele = Nnew;
- size_t k = Nele <= Nele_old ? Nele : Nele_old;
-
- p += k;
- double* q = ptr_data + k;
- while ( q > ptr_data ) *( --q ) = *( --p );
-
- delete[] p;
- }
-
- if ( Nele == 0 ) *this = AlVec( Nnew );
-}
-
-void AlVec::removeElement( size_t Nr ) {
-
- // reorganize vector by filling holes, etc...
- for ( size_t i = Nr; i < Nele - 1; i++ ) ( *this )[i] = ( *this )[i + 1];
-
- // delete obsolete space and shrink vector size by one unit
- delete ( ptr_data + Nele - 1 );
- Nele--;
-}
-
-void AlVec::ShrinkDofs( std::vector<bool>& mask ) {
-
- if ( mask.size() != Nele ) {
- std::cerr << "Beware, mask size different from vector size, skip..." << std::endl;
- return;
- }
-
- size_t offset = 0;
-
- for ( size_t iR = 0; iR < mask.size(); iR++ ) {
- if ( !mask[iR] ) {
- removeElement( iR - offset );
- offset++;
- }
- }
-}
diff --git a/Alignment/AlignSolvTools/src/CLHEPSolver.cpp b/Alignment/AlignSolvTools/src/CLHEPSolver.cpp
index 5241345bb607f1583a9702500298bb2059917dfa..f1042f25dff9f4272a8e5cc588c273e0a3e0e9c1 100755
--- a/Alignment/AlignSolvTools/src/CLHEPSolver.cpp
+++ b/Alignment/AlignSolvTools/src/CLHEPSolver.cpp
@@ -27,15 +27,15 @@ class CLHEPSolver : public extends<GaudiTool, IAlignSolvTool> {
public:
/// Standard constructor
- CLHEPSolver( const std::string& type, const std::string& name, const IInterface* parent );
+ using extends::extends;
using IAlignSolvTool::compute; ///< avoids hiding the original function definitions
/// Solves Ax = b using gsl_linalg_SV_decomp (gsl_matrix * A, gsl_matrix * V, gsl_vector * S, gsl_vector * work)
- bool compute( AlSymMat& symMatrix, AlVec& vector ) const override;
+ StatusCode compute( AlSymMat& symMatrix, AlVec& vector ) const override;
private:
- bool m_useSVD;
- size_t m_numberOfPrintedEigenvalues;
+ Gaudi::Property<size_t> m_numberOfPrintedEigenvalues{this, "NumberOfPrintedEigenvalues", 20};
+ Gaudi::Property<bool> m_useSVD{this, "UseSVD", true};
};
//-----------------------------------------------------------------------------
@@ -51,22 +51,12 @@ DECLARE_COMPONENT( CLHEPSolver )
#include "CLHEP/Matrix/SymMatrix.h"
#include "CLHEP/Matrix/Vector.h"
-//=============================================================================
-// Standard constructor, initializes variables
-//=============================================================================
-CLHEPSolver::CLHEPSolver( const std::string& type, const std::string& name, const IInterface* parent )
- : base_class( type, name, parent ) {
- declareInterface<IAlignSolvTool>( this );
- declareProperty( "NumberOfPrintedEigenvalues", m_numberOfPrintedEigenvalues = 20 );
- declareProperty( "UseSVD", m_useSVD = true );
-}
-
// fix CLHEP bug exposed in gcc 4.3
namespace CLHEP {
double dot( const HepVector& v1, const HepVector& v2 );
}
-bool CLHEPSolver::compute( AlSymMat& symMatrix, AlVec& vector ) const {
+StatusCode CLHEPSolver::compute( AlSymMat& symMatrix, AlVec& vector ) const {
// copy system into CLHEP matrices
size_t size = symMatrix.size();
@@ -74,7 +64,7 @@ bool CLHEPSolver::compute( AlSymMat& symMatrix, AlVec& vector ) const {
CLHEP::HepVector b( size );
for ( size_t irow = 0; irow < size; ++irow ) {
b( irow + 1 ) = vector[irow];
- for ( size_t icol = 0; icol <= irow; ++icol ) A.fast( irow + 1, icol + 1 ) = symMatrix[irow][icol];
+ for ( size_t icol = 0; icol <= irow; ++icol ) A.fast( irow + 1, icol + 1 ) = symMatrix.fast( irow, icol );
}
clock_t starttime = clock();
@@ -129,8 +119,8 @@ bool CLHEPSolver::compute( AlSymMat& symMatrix, AlVec& vector ) const {
// now copy everything back
for ( size_t irow = 0; irow < size; ++irow ) {
vector[irow] = delta( irow + 1 );
- for ( size_t icol = 0; icol <= irow; ++icol ) symMatrix[irow][icol] = cov.fast( irow + 1, icol + 1 );
+ for ( size_t icol = 0; icol <= irow; ++icol ) symMatrix.fast( irow, icol ) = cov.fast( irow + 1, icol + 1 );
}
- return true;
+ return StatusCode::SUCCESS;
}
diff --git a/Alignment/AlignSolvTools/src/DiagSolvTool.cpp b/Alignment/AlignSolvTools/src/DiagSolvTool.cpp
deleted file mode 100755
index b99fbb87c00ad28b90a3a158be95b30078ac197d..0000000000000000000000000000000000000000
--- a/Alignment/AlignSolvTools/src/DiagSolvTool.cpp
+++ /dev/null
@@ -1,237 +0,0 @@
-/*****************************************************************************\
-* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration *
-* *
-* This software is distributed under the terms of the GNU General Public *
-* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
-* *
-* In applying this licence, CERN does not waive the privileges and immunities *
-* granted to it by virtue of its status as an Intergovernmental Organization *
-* or submit itself to any jurisdiction. *
-\*****************************************************************************/
-// Include files
-
-#include <fstream>
-#include <iomanip>
-#include <iostream>
-#include <math.h>
-#include <set>
-#include <stdio.h>
-#include <string>
-
-#include "AlignKernel/AlMat.h"
-#include "AlignKernel/AlSymMat.h"
-#include "AlignKernel/AlVec.h"
-
-// local
-#include "DiagSolvTool.h"
-
-//-----------------------------------------------------------------------------
-// Implementation file for class : DiagSolvTool
-//
-// 2007-06 : Adlene Hicheur
-//-----------------------------------------------------------------------------
-
-// Declaration of the Tool Factory
-DECLARE_COMPONENT( DiagSolvTool )
-
-//=============================================================================
-// Standard constructor, initializes variables
-//=============================================================================
-DiagSolvTool::DiagSolvTool( const std::string& type, const std::string& name, const IInterface* parent )
- : GaudiTool( type, name, parent ) {
- declareInterface<IAlignSolvTool>( this );
- declareProperty( "NumberOfPrintedEigenvalues", m_numberOfPrintedEigenvalues = 20 );
- declareProperty( "EigenValueThreshold", m_eigenValueThreshold = -1 );
- declareProperty( "ApplyScaling", m_applyScaling = false );
- declareProperty( "MinEigenModeChisquare", m_minEigenModeChisquare = -1 );
-}
-
-//=============================================================================
-// Destructor
-//=============================================================================
-DiagSolvTool::~DiagSolvTool() {}
-
-//=============================================================================
-// Initialize
-//=============================================================================
-StatusCode DiagSolvTool::initialize() {
- StatusCode sc = GaudiTool::initialize();
- info() << "EigenValueThreshold = " << m_eigenValueThreshold << endmsg;
- info() << "MinEigenModeChisquare = " << m_minEigenModeChisquare << endmsg;
- return sc;
-}
-
-//=============================================================================
-
-bool DiagSolvTool::compute( AlSymMat& m, AlVec& b ) const {
- IAlignSolvTool::SolutionInfo info;
- StatusCode sc = compute( m, b, info );
- return sc.isSuccess();
-}
-
-namespace {
- struct SortByAbsSecond {
- template <class T>
- bool operator()( const T& lhs, const T& rhs ) const {
- return std::abs( lhs.second ) < std::abs( rhs.second );
- }
- };
-} // namespace
-
-StatusCode DiagSolvTool::compute( AlSymMat& m_bigmatrix, AlVec& m_bigvector,
- IAlignSolvTool::SolutionInfo& solinfo ) const {
- StatusCode sc = StatusCode::SUCCESS;
-
- info() << "Solving with diagonalization method " << endmsg;
-
- char jobz = 'V';
- size_t N = m_bigmatrix.size();
-
- // declare transition matrix + vector to store eigenvalues
- AlMat z( N, N );
- AlVec w( N );
- AlVec rweight( N ), cweight( N );
-
- info() << "After z/w allocation" << endmsg;
-
- // do solving, default diagonalization via LAPACK
-
- int infjob = m_bigmatrix.diagonalize( jobz, w, z );
-
- // diagonalize via GSL
- // m_bigmatrix.diagonalize_GSL(w,z);
-
- info() << "After diagonalization" << endmsg;
-
- // Dump the 10 smallest eigenvalues. We'll sort these really in abs,
- // ignoring the sign.
- {
- // count the number of negative eigenvalues and number smaller than one
- solinfo.numNegativeEigenvalues = solinfo.numSmallEigenvalues = 0;
- std::vector<std::pair<size_t, double>> sortedev( N );
- for ( size_t ipar = 0; ipar < N; ++ipar ) {
- sortedev[ipar] = std::make_pair( ipar, w[ipar] );
- if ( std::abs( w[ipar] ) < 1.0 ) ++solinfo.numSmallEigenvalues;
- if ( w[ipar] < 0.0 ) ++solinfo.numNegativeEigenvalues;
- }
-
- std::sort( sortedev.begin(), sortedev.end(), SortByAbsSecond() );
- std::ostringstream logmessage;
- logmessage << "Smallest eigen values: [ " << std::setprecision( 4 );
- size_t maxpar = std::min( N, m_numberOfPrintedEigenvalues );
- for ( size_t ipar = 0; ipar < maxpar; ++ipar )
- logmessage << sortedev[ipar].second << ( ( ipar < maxpar - 1 ) ? ", " : " ]" );
-
- // consider the smallest eigenvalue
- if ( N > 0 ) {
- // dump the corresponding eigenvector
- if ( N <= m_numberOfPrintedEigenvalues ) {
- logmessage << std::endl;
- logmessage << "Eigenvector for smallest eigenvalue: [ ";
- for ( size_t ipar = 0; ipar < N; ++ipar ) {
- logmessage << z[sortedev[0].first][ipar] << ( ( ipar < N - 1 ) ? ", " : " ]" );
- }
- }
- // figure out which parameter contributes most to this eigenvector
- size_t thepar( 0 );
- for ( size_t ipar = 0; ipar < N; ++ipar ) {
- if ( std::abs( z[sortedev[0].first][thepar] ) < std::abs( z[sortedev[0].first][ipar] ) ) thepar = ipar;
- }
- solinfo.minEigenValue = sortedev[0].second;
- solinfo.weakestModeMaxPar = thepar;
- solinfo.weakestModeMaxParCoef = z[sortedev[0].first][thepar];
- }
- info() << logmessage.str() << endmsg;
- }
-
- if ( infjob == 0 ) {
-
- // Compute bigvector in diagonal basis
- AlVec D( N );
-
- D = z * m_bigvector;
- // Warning: with GSL, the definition of z is transposed:
- // D = z.T()*m_bigvector;
-
- // Issue a warning for each negative eigenvalue
- for ( size_t ipar = 0; ipar < N; ++ipar )
- if ( !( w[ipar] > 0 ) )
- warning() << "Negative eigenvalue (i,val,chi2)= " << ipar << " " << w[ipar] << " "
- << D[ipar] * D[ipar] / w[ipar] << endmsg;
-
- // compute alignment corrections and their variance. first flag
- // modes that we cut away by eigenvalue. (constraints have large
- // negative value, so the 'abs' should do)
- std::vector<bool> keepEigenValue( N, true );
- if ( m_eigenValueThreshold > 0 )
- for ( size_t i = 0; i < N; i++ )
- keepEigenValue[i] = std::abs( w[i] ) > m_eigenValueThreshold ||
- ( m_minEigenModeChisquare > 0 && D[i] * D[i] / w[i] > m_minEigenModeChisquare );
-
- double sumchisqaccepted( 0 ), sumchisqrejected( 0 ), maxchisq( 0 );
- size_t numrejected( 0 );
- std::multiset<double> chisqvalues;
- for ( size_t i = 0; i < N; i++ ) {
- double thischisq = D[i] * D[i] / w[i];
- if ( !keepEigenValue[i] ) {
- info() << "Rejecting eigenvalue: val = " << w[i] << " chisq = " << thischisq << endmsg;
- ++numrejected;
- sumchisqrejected += thischisq;
- } else {
- chisqvalues.insert( thischisq );
- if ( maxchisq < thischisq ) maxchisq = thischisq;
- sumchisqaccepted += thischisq;
- if ( std::abs( w[i] ) < m_eigenValueThreshold )
- info() << "Accepting eigenvalue: val = " << w[i] << " chisq = " << thischisq << endmsg;
- }
- }
-
- solinfo.totalChisq = sumchisqaccepted;
-
- info() << "Number / total chi2 of rejected eigenvalues: " << numrejected << " " << sumchisqrejected << endmsg;
- info() << "Total chi2 of accepted eigenvalues: " << sumchisqaccepted << endmsg;
- info() << "Maximum chi2 of individual mode: " << maxchisq << endmsg;
-
- // because the lagrange constraints lead to negative ev, they also
- // give negative chi2 contributions. that's all fine for the total
- // chi2, but it leads to problems if you want to identify the
- // 'largest' chi2 contribution: such contributions may be
- // compensated by large negative contributions. Let's solve this
- // pragmatically:
- // - make a multiset of all chisq contributions
- // - while the smallest one is negative, remove the smallest and
- // largest from the set, and reinsert the sum
- // - the largest remaining eigenvalue is the one we need;
- while ( !chisqvalues.empty() && *( chisqvalues.begin() ) < 0 && *( chisqvalues.rbegin() ) > 0 ) {
- double a = *( chisqvalues.begin() );
- chisqvalues.erase( chisqvalues.begin() );
- double b = *( chisqvalues.rbegin() );
- chisqvalues.erase( --( chisqvalues.rbegin().base() ) );
- chisqvalues.insert( a + b );
- }
- if ( !chisqvalues.empty() ) solinfo.maxChisqEigenMode = *( chisqvalues.rbegin() );
-
- info() << "Maximum chi2 of individual mode after correcting for lagrange constraints: " << solinfo.maxChisqEigenMode
- << endmsg;
-
- // reset the input
- for ( size_t i = 0; i < N; i++ ) {
- m_bigvector[i] = 0;
- for ( size_t j = 0; j <= i; j++ ) m_bigmatrix[i][j] = 0.;
- }
-
- // now fill it
- for ( size_t k = 0; k < N; ++k )
- if ( keepEigenValue[k] )
- for ( size_t i = 0; i < N; ++i ) {
- m_bigvector[i] += z[k][i] * D[k] / w[k];
- for ( size_t j = 0; j <= i; ++j ) m_bigmatrix[i][j] += ( z[k][i] * z[k][j] / ( w[k] ) );
- }
-
- } else {
- error() << "inversion (diagonalization) of big matrix failed" << endmsg;
- sc = StatusCode::FAILURE;
- }
-
- return sc;
-}
diff --git a/Alignment/AlignSolvTools/src/DiagSolvToolBase.cpp b/Alignment/AlignSolvTools/src/DiagSolvToolBase.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c1d33c2708f3d48cb6533689c8ceedc0b5060346
--- /dev/null
+++ b/Alignment/AlignSolvTools/src/DiagSolvToolBase.cpp
@@ -0,0 +1,199 @@
+/*****************************************************************************\
+* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration *
+* *
+* This software is distributed under the terms of the GNU General Public *
+* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
+* *
+* In applying this licence, CERN does not waive the privileges and immunities *
+* granted to it by virtue of its status as an Intergovernmental Organization *
+* or submit itself to any jurisdiction. *
+\*****************************************************************************/
+// Include files
+
+#include <iomanip>
+#include <iostream>
+#include <set>
+
+// local
+#include "DiagSolvToolBase.h"
+
+//-----------------------------------------------------------------------------
+// Implementation file for class : DiagSolvTool
+//
+// 2007-06 : Adlene Hicheur
+//-----------------------------------------------------------------------------
+
+// Declaration of the Tool Factory
+// DECLARE_COMPONENT( DiagSolvTool )
+
+//=============================================================================
+// Initialize
+//=============================================================================
+StatusCode DiagSolvToolBase::initialize() {
+ StatusCode sc = GaudiTool::initialize();
+ info() << "EigenValueThreshold = " << m_eigenValueThreshold << endmsg;
+ info() << "MinEigenModeChisquare = " << m_minEigenModeChisquare << endmsg;
+ return sc;
+}
+
+namespace {
+ struct SortByAbsSecond {
+ template <class T>
+ bool operator()( const T& lhs, const T& rhs ) const {
+ return std::abs( lhs.second ) < std::abs( rhs.second );
+ }
+ };
+} // namespace
+
+StatusCode DiagSolvToolBase::compute( AlSymMat& bigmatrix, AlVec& bigvector,
+ IAlignSolvTool::SolutionInfo& solinfo ) const {
+ info() << "Solving with diagonalization method" << endmsg;
+
+ // declare transition matrix + vector to store eigenvalues
+ const size_t N = bigmatrix.rows();
+ AlMat U( N, N );
+ AlMat V( N, N );
+ AlVec S( N );
+
+ debug() << "Before diagonalization" << endmsg;
+
+ // defer the diagonalization to whatever tool derives from this.
+ // bigmatrix = U S V^T
+ StatusCode sc = diagonalize( bigmatrix, U, V, S );
+ if ( !sc.isSuccess() ) {
+ error() << "Diagonalization failed" << endmsg;
+ return sc;
+ }
+
+ debug() << "After diagonalization" << endmsg;
+
+ // Compute bigvector in diagonal basis. Here we need the transpose of U, just like it would be in GSL.
+ const Eigen::VectorXd D = U.transpose() * bigvector;
+
+ // Perform some analysis of the eigenvalues.
+ solinfo.numNegativeEigenvalues = solinfo.numSmallEigenvalues = 0;
+ {
+ // count the number of negative eigenvalues and number smaller than one
+ std::vector<std::pair<size_t, double>> sortedev( N );
+ for ( size_t ipar = 0; ipar < N; ++ipar ) {
+ sortedev[ipar] = std::make_pair( ipar, S[ipar] );
+ if ( std::abs( S[ipar] ) < 1.0 ) ++solinfo.numSmallEigenvalues;
+
+ // Eigen and Gsl put the 'sign' of negative eigenvalues in a
+ // difference between U and V. This explains why we do not see
+ // negative Eigenvalues and need to use both U and V in the
+ // solution.
+ const double sign = U.col( ipar ).dot( V.col( ipar ) );
+ if ( !( sign * S[ipar] > 0.0 ) ) {
+ ++solinfo.numNegativeEigenvalues;
+ // Issue a warning for each negative eigenvalue.
+ warning() << "Negative eigenvalue (i,val,chi2)= " << ipar << " " << sign * S( ipar ) << " "
+ << D( ipar ) * D( ipar ) / S( ipar ) << endmsg;
+ }
+ }
+
+ // Dump the 10 smallest eigenvalues. We'll sort these really in
+ // abs, ignoring the sign.
+ std::sort( sortedev.begin(), sortedev.end(), SortByAbsSecond() );
+ std::ostringstream logmessage;
+ logmessage << "Smallest eigen values: [ " << std::setprecision( 4 );
+ size_t maxpar = std::min( N, m_numberOfPrintedEigenvalues.value() );
+ for ( size_t ipar = 0; ipar < maxpar; ++ipar )
+ logmessage << sortedev[ipar].second << ( ( ipar < maxpar - 1 ) ? ", " : " ]" );
+
+ // consider the smallest eigenvalue
+ if ( N > 0 ) {
+ // dump the corresponding eigenvector
+ if ( N <= m_numberOfPrintedEigenvalues ) {
+ logmessage << std::endl;
+ logmessage << "Eigenvector for smallest eigenvalue: [ ";
+ for ( size_t ipar = 0; ipar < N; ++ipar ) {
+ logmessage << U( sortedev[0].first, ipar ) << ( ( ipar < N - 1 ) ? ", " : " ]" );
+ }
+ }
+ // figure out which parameter contributes most to this eigenvector
+ size_t thepar( 0 );
+ for ( size_t ipar = 0; ipar < N; ++ipar ) {
+ if ( std::abs( U( sortedev[0].first, thepar ) ) < std::abs( U( sortedev[0].first, ipar ) ) ) thepar = ipar;
+ }
+ solinfo.minEigenValue = sortedev[0].second;
+ solinfo.weakestModeMaxPar = thepar;
+ solinfo.weakestModeMaxParCoef = U( sortedev[0].first, thepar );
+ }
+ info() << logmessage.str() << endmsg;
+ }
+
+ // compute alignment corrections and their variance. first flag
+ // modes that we cut away by eigenvalue. (constraints have large
+ // negative value, so the 'abs' should do)
+ std::vector<bool> keepEigenValue( N, true );
+ if ( m_eigenValueThreshold > 0 )
+ for ( size_t i = 0; i < N; i++ )
+ keepEigenValue[i] = std::abs( S( i ) ) > m_eigenValueThreshold ||
+ ( m_minEigenModeChisquare > 0 && D( i ) * D( i ) / S( i ) > m_minEigenModeChisquare );
+
+ double sumchisqaccepted( 0 ), sumchisqrejected( 0 ), maxchisq( 0 );
+ size_t numrejected( 0 );
+ std::multiset<double> chisqvalues;
+ for ( size_t i = 0; i < N; i++ ) {
+ double thischisq = D( i ) * D( i ) / S( i );
+ if ( !keepEigenValue[i] ) {
+ info() << "Rejecting eigenvalue: val = " << S( i ) << " chisq = " << thischisq << endmsg;
+ ++numrejected;
+ sumchisqrejected += thischisq;
+ } else {
+ chisqvalues.insert( thischisq );
+ if ( maxchisq < thischisq ) maxchisq = thischisq;
+ sumchisqaccepted += thischisq;
+ if ( std::abs( S( i ) ) < m_eigenValueThreshold )
+ info() << "Accepting low eigenvalue: val = " << S( i ) << " chisq = " << thischisq << endmsg;
+ }
+ }
+
+ solinfo.totalChisq = sumchisqaccepted;
+
+ info() << "Number / total chi2 of rejected eigenvalues: " << numrejected << " " << sumchisqrejected << endmsg;
+ info() << "Total chi2 of accepted eigenvalues: " << sumchisqaccepted << endmsg;
+ info() << "Maximum chi2 of individual mode: " << maxchisq << endmsg;
+
+ // because the lagrange constraints lead to negative ev, they also
+ // give negative chi2 contributions. that's all fine for the total
+ // chi2, but it leads to problems if you want to identify the
+ // 'largest' chi2 contribution: such contributions may be
+ // compensated by large negative contributions. Let's solve this
+ // pragmatically:
+ // - make a multiset of all chisq contributions
+ // - while the smallest one is negative, remove the smallest and
+ // largest from the set, and reinsert the sum
+ // - the largest remaining eigenvalue is the one we need;
+ while ( !chisqvalues.empty() && *( chisqvalues.begin() ) < 0 && *( chisqvalues.rbegin() ) > 0 ) {
+ double a = *( chisqvalues.begin() );
+ chisqvalues.erase( chisqvalues.begin() );
+ double b = *( chisqvalues.rbegin() );
+ chisqvalues.erase( --( chisqvalues.rbegin().base() ) );
+ chisqvalues.insert( a + b );
+ }
+ if ( !chisqvalues.empty() ) solinfo.maxChisqEigenMode = *( chisqvalues.rbegin() );
+
+ info() << "Maximum chi2 of individual mode after correcting for lagrange constraints: " << solinfo.maxChisqEigenMode
+ << endmsg;
+
+ // reset the input
+ for ( size_t i = 0; i < N; i++ ) {
+ bigvector( i ) = 0.;
+ for ( size_t j = 0; j <= i; j++ ) bigmatrix( i, j ) = 0.;
+ }
+
+ // now fill it
+ for ( size_t k = 0; k < N; ++k )
+ if ( keepEigenValue[k] ) {
+ for ( size_t i = 0; i < N; ++i ) {
+ bigvector( i ) += V( i, k ) * D( k ) / S( k );
+ for ( size_t j = 0; j <= i; ++j ) bigmatrix( i, j ) += V( i, k ) * U( j, k ) / S( k );
+ }
+ } else {
+ warning() << "Rejecting eigenvalue: " << k << " " << S( k ) << endmsg;
+ }
+
+ return sc;
+}
diff --git a/Alignment/AlignSolvTools/src/DiagSolvTool.h b/Alignment/AlignSolvTools/src/DiagSolvToolBase.h
old mode 100755
new mode 100644
similarity index 63%
rename from Alignment/AlignSolvTools/src/DiagSolvTool.h
rename to Alignment/AlignSolvTools/src/DiagSolvToolBase.h
index b168cc43b15c85d11ea10d98acd8998944c176ef..58d19f73467f4e940b769d5b0abb82625f87bc7d
--- a/Alignment/AlignSolvTools/src/DiagSolvTool.h
+++ b/Alignment/AlignSolvTools/src/DiagSolvToolBase.h
@@ -8,37 +8,34 @@
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
-// $Id: DiagSolvTool.h,v 1.11 2010-03-02 15:55:26 wouter Exp $
-#ifndef ALIGNSOLVTOOLS_DIAGSOLVTOOL_H
-#define ALIGNSOLVTOOLS_DIAGSOLVTOOL_H 1
+#pragma once
// Include files
// from Gaudi
#include "AlignmentInterfaces/IAlignSolvTool.h" // Interface
#include "GaudiAlg/GaudiTool.h"
-class AlSymMat;
-class AlMat;
-class AlVec;
-
/** @class DiagSolvTool DiagSolvTool.h AlignSolvTools/DiagSolvTool.h
*
*
* @author Adlene Hicheur
* @date 2007-06
*/
-class DiagSolvTool : public GaudiTool, virtual public IAlignSolvTool {
+
+class DiagSolvToolBase : public extends<GaudiTool, IAlignSolvTool> {
public:
/// Standard constructor
- DiagSolvTool( const std::string& type, const std::string& name, const IInterface* parent );
-
- virtual ~DiagSolvTool(); ///< Destructor
+ using extends::extends;
using IAlignSolvTool::compute; ///< avoids hiding the original function definitions
+
///< Call method to compute the solution, get symetric
/// matrix and vector in input and replaces them by
/// inverse matrix and solution vector in output.
- bool compute( AlSymMat&, AlVec& ) const override;
+ StatusCode compute( AlSymMat& A, AlVec& b ) const override {
+ IAlignSolvTool::SolutionInfo info;
+ return compute( A, b, info );
+ }
// Solves the system Ax=b for x. Return A=A^{-1} and
// b=x. Furthermore, returns the total chi2 (b^T*A^{-1}*b) and the
@@ -46,19 +43,16 @@ public:
StatusCode compute( AlSymMat& A, AlVec& b, IAlignSolvTool::SolutionInfo& ) const override;
StatusCode initialize() override;
- // enum Solver {SPMINV,
- // DIAG,
- // MINRES};
-protected:
-private:
- // parameters
- size_t m_numberOfPrintedEigenvalues;
- bool m_applyScaling;
- double m_eigenValueThreshold;
- double m_minEigenModeChisquare;
+ // Diagonalizes A symmetric matrix A such that A = U D V^T
+ // * S is the vector with eigenvalues, the diagonal elements of D
+ // * the columns of U are the eigen vectors of A
+ // * for symmetric matrices a solution exists for which V=U, but ... both Eigen and Gsl do not obey that!
+ virtual StatusCode diagonalize( const AlSymMat& A, AlMat& U, AlMat& V, AlVec& S ) const = 0;
- void MonitorDiag( AlMat& z, AlVec& w, AlVec& D, double scale );
- int SolvMinRes( AlSymMat& m_bigmatrix, AlVec& m_bigvector, AlVec& result );
+protected:
+ Gaudi::Property<size_t> m_numberOfPrintedEigenvalues{this, "NumberOfPrintedEigenvalues", 20};
+ Gaudi::Property<double> m_eigenValueThreshold{this, "EigenValueThreshold", -1};
+ Gaudi::Property<bool> m_applyScaling{this, "ApplyScaling", false};
+ Gaudi::Property<double> m_minEigenModeChisquare{this, "MinEigenModeChisquare", -1};
};
-#endif // ALIGNSOLVTOOLS_ALIGNSOLVTOOL_H
diff --git a/Alignment/AlignSolvTools/src/EigenDiagSolvTool.cpp b/Alignment/AlignSolvTools/src/EigenDiagSolvTool.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..55a162cc52200be2e5d270bca05b8894ef5154ee
--- /dev/null
+++ b/Alignment/AlignSolvTools/src/EigenDiagSolvTool.cpp
@@ -0,0 +1,35 @@
+/*****************************************************************************\
+* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration *
+* *
+* This software is distributed under the terms of the GNU General Public *
+* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
+* *
+* In applying this licence, CERN does not waive the privileges and immunities *
+* granted to it by virtue of its status as an Intergovernmental Organization *
+* or submit itself to any jurisdiction. *
+\*****************************************************************************/
+// Include files
+
+#include "DiagSolvToolBase.h"
+
+class EigenDiagSolvTool : public DiagSolvToolBase {
+public:
+ using DiagSolvToolBase::DiagSolvToolBase;
+
+ using IAlignSolvTool::compute; ///< avoids hiding the original function definitions
+
+ StatusCode diagonalize( const AlSymMat& A, AlMat& U, AlMat& V, AlVec& S ) const override;
+};
+
+// Declaration of the Tool Factory
+DECLARE_COMPONENT( EigenDiagSolvTool )
+
+StatusCode EigenDiagSolvTool::diagonalize( const AlSymMat& A, AlMat& U, AlMat& V, AlVec& S ) const {
+ // There are several options in Eigen. We just start with one.
+ Eigen::MatrixXd eigenA = A.toMatrix();
+ Eigen::BDCSVD<Eigen::MatrixXd> svdcalculator{eigenA, Eigen::ComputeFullU | Eigen::ComputeFullV};
+ S = svdcalculator.singularValues();
+ U = svdcalculator.matrixU();
+ V = svdcalculator.matrixV();
+ return StatusCode::SUCCESS;
+}
diff --git a/Alignment/AlignSolvTools/src/GslDiagSolvTool.cpp b/Alignment/AlignSolvTools/src/GslDiagSolvTool.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c693e46031ad9663c62c225e83f88764edbdd948
--- /dev/null
+++ b/Alignment/AlignSolvTools/src/GslDiagSolvTool.cpp
@@ -0,0 +1,117 @@
+/*****************************************************************************\
+* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration *
+* *
+* This software is distributed under the terms of the GNU General Public *
+* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
+* *
+* In applying this licence, CERN does not waive the privileges and immunities *
+* granted to it by virtue of its status as an Intergovernmental Organization *
+* or submit itself to any jurisdiction. *
+\*****************************************************************************/
+// Include files
+
+#include <cmath>
+
+// from GSL
+#include "gsl/gsl_linalg.h"
+
+// Interface
+#include "DiagSolvToolBase.h"
+
+class GslDiagSolvTool : public DiagSolvToolBase {
+public:
+ using DiagSolvToolBase::DiagSolvToolBase;
+
+ using IAlignSolvTool::compute; ///< avoids hiding the original function definitions
+
+ StatusCode diagonalize( const AlSymMat& A, AlMat& U, AlMat& V, AlVec& S ) const override;
+
+private:
+ Gaudi::Property<bool> m_svdJacobi{this, "SVDJacobi", false}; ///< Use Jacobi SVD?
+};
+
+// Declaration of the Tool Factory
+DECLARE_COMPONENT( GslDiagSolvTool )
+
+/*
+namespace {
+
+ std::ostream& operator<<( std::ostream& lhs, const gsl_vector& rhs ) {
+ lhs << std::endl;
+ lhs << "[ " << std::endl;
+ for ( unsigned i = 0u; i < rhs.size; ++i ) {
+ lhs << " " << ( *gsl_vector_const_ptr( &rhs, i ) ) << " ";
+ lhs << std::endl;
+ }
+ lhs << " ]" << std::endl;
+ return lhs;
+ }
+
+ std::ostream& operator<<( std::ostream& lhs, const gsl_matrix& rhs ) {
+ lhs << std::endl;
+ lhs << "[ " << std::endl;
+ for ( unsigned i = 0u; i < rhs.size1; ++i ) {
+ for ( unsigned j = 0u; j < rhs.size2; ++j ) { lhs << " " << ( *gsl_matrix_const_ptr( &rhs, i, j ) ) << " "; }
+ lhs << std::endl;
+ }
+ lhs << " ]" << std::endl;
+ return lhs;
+ }
+} // namespace
+*/
+
+StatusCode GslDiagSolvTool::diagonalize( const AlSymMat& A, AlMat& U, AlMat& V, AlVec& S ) const {
+
+ const auto N = A.rows();
+
+ debug() << "Size of AlSymMat A = " << N << endmsg;
+ /// Allocate matrix A
+ gsl_matrix* matrixA = gsl_matrix_alloc( N, N );
+ debug() << "Size of gsl_matrix A = " << matrixA->size1 << endmsg;
+ /// Fill matrix A
+ for ( unsigned i = 0u; i < N; ++i ) {
+ for ( unsigned j = 0u; j < N; ++j ) {
+ debug() << "Element (i,j) of AlSymMat A = " << A( i, j ) << endmsg;
+ gsl_matrix_set( matrixA, i, j, A( i, j ) );
+ debug() << "Element (i,j) of gsl_matrix A = " << gsl_matrix_get( matrixA, i, j ) << endmsg;
+ }
+ }
+
+ /// Allocate required matrix, vector and workspace
+ gsl_matrix* matrixV = gsl_matrix_alloc( N, N );
+ gsl_vector* vectorS = gsl_vector_alloc( N );
+ gsl_vector* vectorW = gsl_vector_alloc( N );
+ // gsl_matrix* matrixU = matrixA ;
+
+ debug() << "Factorising matrix A" << endmsg;
+ /// Factorise A into the SVD A = USV^T. Note matrix A is replaced with matrix U.
+ /// GSL returns 0 if successful
+ int factorised = 1;
+ if ( !m_svdJacobi.value() ) {
+ info() << "==> SVD standard" << endmsg;
+ factorised = gsl_linalg_SV_decomp( matrixA, matrixV, vectorS, vectorW );
+ } else {
+ info() << "==> SVD Jacobi" << endmsg;
+ factorised = gsl_linalg_SV_decomp_jacobi( matrixA, matrixV, vectorS );
+ }
+ if ( factorised != 0 ) {
+ error() << "==> Couldn't factorise matrix" << endmsg;
+ return StatusCode::FAILURE;
+ }
+ debug() << "Done factorising matrix A" << endmsg;
+
+ //
+ for ( size_t irow = 0; irow < N; ++irow ) S( irow ) = gsl_vector_get( vectorS, irow );
+ for ( size_t irow = 0; irow < N; ++irow )
+ for ( size_t icol = 0; icol < N; ++icol ) {
+ U( irow, icol ) = gsl_matrix_get( matrixA, irow, icol );
+ V( irow, icol ) = gsl_matrix_get( matrixV, irow, icol );
+ }
+
+ gsl_matrix_free( matrixA );
+ gsl_matrix_free( matrixV );
+ gsl_vector_free( vectorS );
+ gsl_vector_free( vectorW );
+
+ return StatusCode::SUCCESS;
+}
diff --git a/Alignment/AlignSolvTools/src/SolvExample.cpp b/Alignment/AlignSolvTools/src/SolvExample.cpp
index 1ecc435d220435678af36b9129a652f469198b9c..c20e7174075756b1902283e08c5c846fa9e80b3d 100755
--- a/Alignment/AlignSolvTools/src/SolvExample.cpp
+++ b/Alignment/AlignSolvTools/src/SolvExample.cpp
@@ -8,15 +8,6 @@
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
-// Include files
-
-// local
-#include "SolvExample.h"
-
-// Matrix, vector classes
-#include "AlignKernel/AlMat.h"
-#include "AlignKernel/AlSymMat.h"
-#include "AlignKernel/AlVec.h"
//-----------------------------------------------------------------------------
// Implementation file for class : SolvExample
@@ -24,35 +15,32 @@
// 2007-03-07 : Adlene Hicheur
//-----------------------------------------------------------------------------
-// Declaration of the Algorithm Factory
-DECLARE_COMPONENT( SolvExample )
+// Include files
+// from Gaudi
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "GaudiKernel/ToolHandle.h"
-//=============================================================================
-// Standard constructor, initializes variables
-//=============================================================================
-SolvExample::SolvExample( const std::string& name, ISvcLocator* pSvcLocator ) : GaudiAlgorithm( name, pSvcLocator ) {
- declareProperty( "SolvMethod", m_solvtoolname = "SpmInvTool" );
-}
-//=============================================================================
-// Destructor
-//=============================================================================
-SolvExample::~SolvExample() {}
+// Solver
+#include "AlignmentInterfaces/IAlignSolvTool.h"
-//=============================================================================
-// Initialization
-//=============================================================================
-StatusCode SolvExample::initialize() {
- StatusCode sc = GaudiAlgorithm::initialize(); // must be executed first
- if ( sc.isFailure() ) return sc; // error printed already by GaudiAlgorithm
+class IAlignSolvTool;
- // solver tool
+class SolvExample : public GaudiAlgorithm {
+public:
+ /// Standard constructor
+ using GaudiAlgorithm::GaudiAlgorithm;
+ StatusCode execute() override; ///< Algorithm execution
+private:
+ ToolHandle<IAlignSolvTool> m_solver{this, "SolvMethod", "SpmInvTool"};
+};
- m_solver = tool<IAlignSolvTool>( m_solvtoolname );
-
- debug() << "==> Initialize" << endmsg;
+// Matrix, vector classes
+#include "AlignKernel/AlMat.h"
+#include "AlignKernel/AlSymMat.h"
+#include "AlignKernel/AlVec.h"
- return StatusCode::SUCCESS;
-}
+// Declaration of the Algorithm Factory
+DECLARE_COMPONENT( SolvExample )
//=============================================================================
// Main execution
@@ -75,7 +63,7 @@ StatusCode SolvExample::execute() {
}
b -= c;
- debug() << "scalar product c*d..." << ( b * c ) << endmsg;
+ debug() << "scalar product c*d..." << b.dot( c ) << endmsg;
AlSymMat md( Nvec );
AlSymMat ma( Nvec );
@@ -83,9 +71,9 @@ StatusCode SolvExample::execute() {
for ( unsigned i = 0; i < md.size(); ++i ) {
debug() << "AlSymMat md values line " << i << " -- ";
for ( unsigned j = 0; j <= i; ++j ) {
- md[i][j] = i - j;
- if ( i == j ) ma[i][j] = 3;
- debug() << md[i][j] << " - ";
+ md( i, j ) = i - j;
+ if ( i == j ) ma( i, j ) = 3;
+ debug() << md( i, j ) << " - ";
}
debug() << "finish line" << endmsg;
}
@@ -98,12 +86,11 @@ StatusCode SolvExample::execute() {
for ( unsigned i = 0; i < mdT.size(); ++i ) {
debug() << "AlSymMat mdT values line " << i << " -- ";
- for ( unsigned j = 0; j < mdT.size(); ++j ) { debug() << mdT[i][j] << " - "; }
+ for ( unsigned j = 0; j < mdT.size(); ++j ) { debug() << mdT( i, j ) << " - "; }
debug() << endmsg;
}
- AlVec e;
- e = md * c;
+ AlVec e = md * c;
debug() << "md*c Nelem..." << e.size() << endmsg;
debug() << "md*c values..." << endmsg;
for ( unsigned i = 0; i < e.size(); ++i ) debug() << e[i] << " - ";
@@ -113,32 +100,32 @@ StatusCode SolvExample::execute() {
AlMat mk = ma * me;
AlMat mj = me * ma;
- for ( size_t i = 0; i < me.nrow(); i++ ) {
+ for ( auto i = 0; i < me.rows(); i++ ) {
debug() << "AlMat me values line " << i << " -- ";
- for ( size_t j = 0; j < me.ncol(); j++ ) { debug() << me[i][j] << " - "; }
+ for ( auto j = 0; j < me.cols(); j++ ) { debug() << me( i, j ) << " - "; }
debug() << endmsg;
}
- for ( size_t i = 0; i < mk.nrow(); i++ ) {
+ for ( auto i = 0; i < mk.rows(); i++ ) {
debug() << "AlSymMat*AlMat mk=ma*me values line " << i << " -- ";
- for ( size_t j = 0; j < mk.ncol(); j++ ) { debug() << mk[i][j] << " - "; }
+ for ( auto j = 0; j < mk.cols(); j++ ) { debug() << mk( i, j ) << " - "; }
debug() << endmsg;
}
- for ( size_t i = 0; i < mj.nrow(); i++ ) {
+ for ( auto i = 0; i < mj.rows(); i++ ) {
debug() << "AlMat*AlSymMat mj=me*ma values line " << i << " -- ";
- for ( size_t j = 0; j < mj.ncol(); j++ ) { debug() << mj[i][j] << " - "; }
+ for ( auto j = 0; j < mj.cols(); j++ ) { debug() << mj( i, j ) << " - "; }
debug() << endmsg;
}
me *= 5.;
- AlSymMat mf( me.nrow() );
+ AlSymMat mf( me.rows() );
- for ( size_t i = 0; i < me.nrow(); i++ ) {
+ for ( auto i = 0; i < me.rows(); i++ ) {
debug() << "AlMat me*5 values line " << i << " -- ";
- for ( size_t j = 0; j < me.ncol(); j++ ) {
- debug() << me[i][j] << " - ";
- mf[i][j] = me[i][j];
+ for ( auto j = 0; j < me.cols(); j++ ) {
+ debug() << me( i, j ) << " - ";
+ mf( i, j ) = me( i, j );
}
debug() << endmsg;
}
@@ -150,7 +137,7 @@ StatusCode SolvExample::execute() {
int ierr = 0;
- mf.invert( ierr );
+ // mf.invert(ierr);
if ( ierr == 0 ) {
@@ -158,32 +145,30 @@ StatusCode SolvExample::execute() {
for ( unsigned i = 0; i < mf.size(); ++i ) {
debug() << "AlSymMat mf values line " << i << " -- ";
- for ( unsigned j = 0; j < mf.size(); ++j ) { debug() << mf[i][j] << " - "; }
+ for ( unsigned j = 0; j < mf.size(); ++j ) { debug() << mf( i, j ) << " - "; }
debug() << endmsg;
}
debug() << "AlSymMat me*me^-1 values line " << endmsg;
- for ( size_t i = 0; i < ( mf * mfi ).nrow(); i++ ) {
+ for ( auto i = 0; i < ( mf * mfi ).rows(); i++ ) {
debug() << "AlSymMat me*me^-1 values line " << i << " -- ";
- for ( size_t j = 0; j < ( mf * mfi ).ncol(); j++ ) { debug() << ( mf * mfi )[i][j] << " - "; }
+ for ( auto j = 0; j < ( mf * mfi ).cols(); j++ ) { debug() << ( mf * mfi )( i, j ) << " - "; }
debug() << endmsg;
}
}
debug() << "Now using indirect inversion via diagonalization..." << endmsg;
- int info = 999;
- char jobz = 'V';
// declare transition matrix + vector to store eigenvalues
- size_t dim = mf.size();
+ const size_t dim = mf.size();
AlMat z( dim, dim );
AlVec w( dim );
AlSymMat mfi_bis( mfi );
- info = mfi.diagonalize( jobz, w, z );
+ const auto info = mfi.diagonalize_GSL( w, z );
if ( info == 0 ) {
debug() << "*** successful diagonalization ***" << endmsg;
@@ -192,7 +177,7 @@ StatusCode SolvExample::execute() {
for ( unsigned i = 0; i < invmat.size(); ++i ) {
debug() << w[i] << " " << endmsg;
for ( unsigned j = 0; j <= i; ++j ) {
- for ( size_t k = 0; k < dim; k++ ) invmat[i][j] = invmat[i][j] + ( z[k][i] * z[k][j] / w[k] );
+ for ( size_t k = 0; k < dim; k++ ) invmat( i, j ) = invmat( i, j ) + ( z( k, i ) * z( k, j ) / w[k] );
}
}
@@ -200,30 +185,17 @@ StatusCode SolvExample::execute() {
for ( unsigned i = 0; i < invmat.size(); ++i ) {
// debug()<<"AlSymMat me*me^-1 values line "<<i<<" -- ";
- for ( unsigned j = 0; j < invmat.size(); ++j ) { debug() << ( me * invmat )[i][j] << " "; }
+ for ( unsigned j = 0; j < invmat.size(); ++j ) { debug() << ( me * invmat )( i, j ) << " "; }
debug() << endmsg;
}
}
// Calling the solver
- m_solver->compute( mf, e );
+ m_solver->compute( mf, e ).ignore();
debug() << "Printing solution of mf*X = e..." << endmsg;
for ( unsigned i = 0; i < e.size(); ++i ) debug() << e[i] << " - ";
debug() << endmsg;
return StatusCode::SUCCESS;
}
-
-//=============================================================================
-// Finalize
-//=============================================================================
-StatusCode SolvExample::finalize() {
-
- debug() << "==> Finalize" << endmsg;
- // finalize the tool to close it properly
- m_solver->finalize().ignore();
- return GaudiAlgorithm::finalize(); // must be called after all other actions
-}
-
-//=============================================================================
diff --git a/Alignment/AlignSolvTools/src/SparseSolver.cpp b/Alignment/AlignSolvTools/src/SparseSolver.cpp
index 2ec4b5aad7ea9b1bd6f386a4f5075fe649e3416b..55edeb58120d94170e596dee392866db8eb25948 100644
--- a/Alignment/AlignSolvTools/src/SparseSolver.cpp
+++ b/Alignment/AlignSolvTools/src/SparseSolver.cpp
@@ -23,17 +23,15 @@
* @date 2007-07-24
*/
-class SparseSolver : public GaudiTool, virtual public IAlignSolvTool {
+class SparseSolver : public extends<GaudiTool, IAlignSolvTool> {
public:
/// Standard constructor
- SparseSolver( const std::string& type, const std::string& name, const IInterface* parent );
-
- virtual ~SparseSolver(); ///< Destructor
+ using extends::extends;
using IAlignSolvTool::compute; ///< avoids hiding the original function definitions
- /// Solves Ax = b using gsl_linalg_SV_decomp (gsl_matrix * A, gsl_matrix * V, gsl_vector * S, gsl_vector * work)
- bool compute( AlSymMat& symMatrix, AlVec& vector ) const override;
+ /// Solves Ax = b
+ StatusCode compute( AlSymMat& symMatrix, AlVec& vector ) const override;
private:
};
@@ -52,19 +50,7 @@ DECLARE_COMPONENT( SparseSolver )
#include "TMatrixDSym.h"
#include "TVectorD.h"
-//=============================================================================
-// Standard constructor, initializes variables
-//=============================================================================
-SparseSolver::SparseSolver( const std::string& type, const std::string& name, const IInterface* parent )
- : GaudiTool( type, name, parent ) {
- declareInterface<IAlignSolvTool>( this );
-}
-//=============================================================================
-// Destructor
-//=============================================================================
-SparseSolver::~SparseSolver() {}
-
-bool SparseSolver::compute( AlSymMat& symMatrix, AlVec& vector ) const {
+StatusCode SparseSolver::compute( AlSymMat& symMatrix, AlVec& vector ) const {
// create a sparce matrix. I think that this is most efficient via
// TMatrix, because the structure of the TSparseMatrix isn't
// entirely self evident.
@@ -98,5 +84,5 @@ bool SparseSolver::compute( AlSymMat& symMatrix, AlVec& vector ) const {
tmp.fast( irow, irow ) = c != 0 ? 1 / c : 0;
}
symMatrix = tmp;
- return ierr ? false : true;
+ return ierr ? StatusCode::FAILURE : StatusCode::SUCCESS;
}
diff --git a/Alignment/AlignSolvTools/src/SpmInvTool.cpp b/Alignment/AlignSolvTools/src/SpmInvTool.cpp
index 99bdc222da1ac74da5cb3cc31e7655fa78c338cb..5aa56e8603f0e32c71c5fdc53d4f638621bedd17 100755
--- a/Alignment/AlignSolvTools/src/SpmInvTool.cpp
+++ b/Alignment/AlignSolvTools/src/SpmInvTool.cpp
@@ -8,87 +8,78 @@
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
-// Include files
+/** @class SpmInvTool SpmInvTool.cpp AlignSolvTools/SpmInvTool.cpp
+ *
+ *
+ * @author Adlene Hicheur
+ * @date 2007-06
+ */
-#include <fstream>
-#include <iomanip>
-#include <iostream>
-#include <math.h>
-#include <stdio.h>
-#include <string>
+// Include files
+#include <vector>
// from Gaudi
-#include "GaudiKernel/INTupleSvc.h"
-#include "GaudiKernel/NTuple.h"
-#include "GaudiKernel/SmartDataPtr.h"
+#include "GaudiAlg/GaudiTool.h"
+
+// Interface
+#include "AlignmentInterfaces/IAlignSolvTool.h"
+
+class SpmInvTool : public extends<GaudiTool, IAlignSolvTool> {
+public:
+ /// Standard constructor
+ using extends::extends;
+
+ using IAlignSolvTool::compute; ///< avoids hiding the original function definitions
+ // Call method to compute the solution, get symetric matrix and
+ // vector in input and replaces them by inverse matrix and solution vector in output.
+ StatusCode compute( AlSymMat&, AlVec& ) const override;
-// local
-#include "SpmInvTool.h"
+private:
+ // parameters
-//-----------------------------------------------------------------------------
-// Implementation file for class : SpmInvTool
-//
-// 2007-06 : Adlene Hicheur
-//-----------------------------------------------------------------------------
+ int SolvSpmInv( AlSymMat& M, AlVec& B );
+
+ void Precond( AlSymMat& M, AlVec& r, AlVec& c, bool equIn );
+};
// Declaration of the Tool Factory
DECLARE_COMPONENT( SpmInvTool )
-//=============================================================================
-// Standard constructor, initializes variables
-//=============================================================================
-SpmInvTool::SpmInvTool( const std::string& type, const std::string& name, const IInterface* parent )
- : GaudiTool( type, name, parent ) {
- declareInterface<IAlignSolvTool>( this );
-}
-//=============================================================================
-// Destructor
-//=============================================================================
-SpmInvTool::~SpmInvTool() {}
-
//=============================================================================
-bool SpmInvTool::compute( AlSymMat& m, AlVec& b ) const {
+StatusCode SpmInvTool::compute( AlSymMat& m, AlVec& b ) const {
if ( m.size() > 0 ) {
const_cast<SpmInvTool*>( this )->SolvSpmInv( m, b );
- return true;
+ return StatusCode::SUCCESS;
} else {
error() << "Error, null matrix size, don't call the solving!" << endmsg;
- return false;
+ return StatusCode::FAILURE;
}
}
int SpmInvTool::SolvSpmInv( AlSymMat& M, AlVec& B ) {
- AlMat V( M.size(), M.size() );
- V.copyS( M );
+ AlMat V = M.toMatrix();
int N = M.size();
- AlVec diag( N );
- bool* flag = new bool[N];
+ AlVec diag( N );
+ std::vector<bool> flag( N, true );
int i, j, jj, k, nrank;
double vkk;
double eps = 1e-16;
- AlVec r( N );
- AlVec c( N );
-
- AlVec temp( N );
-
- for ( i = 0; i < N; i++ ) {
- flag[i] = true;
-
- // for (j=0; j<=i; j++) {if (V[j][i] == 0) {V[j][i] = V[i][j];}}
- }
+ AlVec r{N};
+ AlVec c{N};
+ AlVec temp{N};
// Small loop for matrix equilibration (gives a better conditioning)
for ( i = 0; i < N; i++ ) {
for ( j = 0; j < N; j++ ) {
- if ( fabs( V[i][j] ) >= r[i] ) r[i] = fabs( V[i][j] ); // Max elemt of row i
- if ( fabs( V[j][i] ) >= c[i] ) c[i] = fabs( V[j][i] ); // Max elemt of column i
+ if ( fabs( V( i, j ) ) >= r[i] ) r[i] = fabs( V( i, j ) ); // Max elemt of row i
+ if ( fabs( V( j, i ) ) >= c[i] ) c[i] = fabs( V( j, i ) ); // Max elemt of column i
}
}
for ( i = 0; i < N; i++ ) {
@@ -98,13 +89,13 @@ int SpmInvTool::SolvSpmInv( AlSymMat& M, AlVec& B ) {
for ( i = 0; i < N; i++ ) // Equilibrate the V matrix
{
- for ( j = 0; j < N; j++ ) { V[i][j] = sqrt( r[i] ) * V[i][j] * sqrt( c[j] ); }
+ for ( j = 0; j < N; j++ ) { V( i, j ) = sqrt( r[i] ) * V( i, j ) * sqrt( c[j] ); }
}
nrank = 0;
// save diagonal elem absolute values
- for ( i = 0; i < N; i++ ) { diag[i] = fabs( V[i][i] ); }
+ for ( i = 0; i < N; i++ ) { diag[i] = fabs( V( i, i ) ); }
for ( i = 0; i < N; i++ ) {
vkk = 0.0;
@@ -112,8 +103,8 @@ int SpmInvTool::SolvSpmInv( AlSymMat& M, AlVec& B ) {
for ( j = 0; j < N; j++ ) // First look for the pivot, ie max unused diagonal element
{
- if ( flag[j] && ( fabs( V[j][j] ) > std::max( fabs( vkk ), eps * diag[j] ) ) ) {
- vkk = V[j][j];
+ if ( flag[j] && ( fabs( V( j, j ) ) > std::max( fabs( vkk ), eps * diag[j] ) ) ) {
+ vkk = V( j, j );
k = j;
}
}
@@ -122,16 +113,16 @@ int SpmInvTool::SolvSpmInv( AlSymMat& M, AlVec& B ) {
{
debug() << "Pivot value :" << vkk << endmsg;
nrank++;
- flag[k] = false; // This value is used
- vkk = 1.0 / vkk;
- V[k][k] = -vkk; // Replace pivot by its inverse
+ flag[k] = false; // This value is used
+ vkk = 1.0 / vkk;
+ V( k, k ) = -vkk; // Replace pivot by its inverse
for ( j = 0; j < N; j++ ) {
for ( jj = 0; jj < N; jj++ ) {
if ( j != k && jj != k ) // Other elements (!!! do them first as you use old Vk][j]'s !!!)
{
- V[j][jj] = V[j][jj] - vkk * V[j][k] * V[k][jj];
- // V[j][jj] = V[j][jj] + vkk*V[j][k]*V[k][jj];
+ V( j, jj ) = V( j, jj ) - vkk * V( j, k ) * V( k, jj );
+ // V(j,jj) = V(j,jj) + vkk*V(j,k)*V(k,jj);
}
}
}
@@ -139,8 +130,8 @@ int SpmInvTool::SolvSpmInv( AlSymMat& M, AlVec& B ) {
for ( j = 0; j < N; j++ ) {
if ( j != k ) // Pivot row or column elements
{
- V[j][k] = ( V[j][k] ) * vkk; // Column
- V[k][j] = ( V[k][j] ) * vkk; // Line
+ V( j, k ) = ( V( j, k ) ) * vkk; // Column
+ V( k, j ) = ( V( k, j ) ) * vkk; // Line
}
}
} else // No more pivot value (clear those elements)
@@ -150,7 +141,7 @@ int SpmInvTool::SolvSpmInv( AlSymMat& M, AlVec& B ) {
B[j] = 0.0;
for ( k = 0; k <= j; k++ ) {
- V[j][k] = 0.0;
+ V( j, k ) = 0.0;
// V[k][j] = 0.0;
}
}
@@ -161,7 +152,7 @@ int SpmInvTool::SolvSpmInv( AlSymMat& M, AlVec& B ) {
}
for ( i = 0; i < N; i++ ) // Correct matrix V
{
- for ( j = 0; j < N; j++ ) { V[i][j] = sqrt( c[i] ) * V[i][j] * sqrt( r[j] ); }
+ for ( j = 0; j < N; j++ ) { V( i, j ) = sqrt( c[i] ) * V( i, j ) * sqrt( r[j] ); }
}
for ( j = 0; j < N; j++ ) {
@@ -169,23 +160,20 @@ int SpmInvTool::SolvSpmInv( AlSymMat& M, AlVec& B ) {
for ( jj = 0; jj < N; jj++ ) // Reverse matrix elements
{
- V[j][jj] = -V[j][jj];
- M[j][jj] = V[j][jj];
- temp[j] += V[j][jj] * B[jj];
+ V( j, jj ) = -V( j, jj );
+ M( j, jj ) = V( j, jj );
+ temp[j] += V( j, jj ) * B[jj];
}
}
for ( j = 0; j < N; j++ ) { B[j] = temp[j]; } // The final result
- delete[] flag;
return nrank;
}
// Routine to equilibrate the matrix for better conditioning
void SpmInvTool::Precond( AlSymMat& M, AlVec& r, AlVec& c, bool equIn = true ) {
-
- AlMat V( M.size(), M.size() );
- V.copyS( M );
+ AlMat V = M.toMatrix();
int N = M.size();
@@ -199,8 +187,8 @@ void SpmInvTool::Precond( AlSymMat& M, AlVec& r, AlVec& c, bool equIn = true ) {
for ( i = 0; i < N; i++ ) {
for ( j = 0; j < N; j++ ) {
- if ( fabs( V[i][j] ) >= r[i] ) r[i] = fabs( V[i][j] ); // Max elemt of row i
- if ( fabs( V[j][i] ) >= c[i] ) c[i] = fabs( V[j][i] ); // Max elemt of column i
+ if ( fabs( V( i, j ) ) >= r[i] ) r[i] = fabs( V( i, j ) ); // Max elemt of row i
+ if ( fabs( V( j, i ) ) >= c[i] ) c[i] = fabs( V( j, i ) ); // Max elemt of column i
}
}
@@ -213,8 +201,8 @@ void SpmInvTool::Precond( AlSymMat& M, AlVec& r, AlVec& c, bool equIn = true ) {
{
for ( j = 0; j < N; j++ ) {
- V[i][j] = r[i] * V[i][j] * c[j];
- M[i][j] = V[i][j];
+ V( i, j ) = r[i] * V( i, j ) * c[j];
+ M( i, j ) = V( i, j );
}
}
@@ -224,8 +212,8 @@ void SpmInvTool::Precond( AlSymMat& M, AlVec& r, AlVec& c, bool equIn = true ) {
{
for ( j = 0; j < N; j++ ) {
- V[i][j] = c[i] * V[i][j] * r[j];
- M[i][j] = V[i][j];
+ V( i, j ) = c[i] * V( i, j ) * r[j];
+ M( i, j ) = V( i, j );
}
}
}
diff --git a/Alignment/AlignSolvTools/src/SpmInvTool.h b/Alignment/AlignSolvTools/src/SpmInvTool.h
deleted file mode 100755
index b94aee5ab440bdaf0af8107ebc524c1ab9bd9355..0000000000000000000000000000000000000000
--- a/Alignment/AlignSolvTools/src/SpmInvTool.h
+++ /dev/null
@@ -1,49 +0,0 @@
-/*****************************************************************************\
-* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration *
-* *
-* This software is distributed under the terms of the GNU General Public *
-* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
-* *
-* In applying this licence, CERN does not waive the privileges and immunities *
-* granted to it by virtue of its status as an Intergovernmental Organization *
-* or submit itself to any jurisdiction. *
-\*****************************************************************************/
-// $Id: SpmInvTool.h,v 1.4 2009-08-16 14:16:24 wouter Exp $
-#ifndef ALIGNSOLVTOOLS_SPMINVTOOL_H
-#define ALIGNSOLVTOOLS_SPMINVTOOL_H 1
-
-// Include files
-// from Gaudi
-#include "GaudiAlg/GaudiTool.h"
-
-// Interface
-#include "AlignmentInterfaces/IAlignSolvTool.h"
-
-/** @class SpmInvTool SpmInvTool.h AlignSolvTools/SpmInvTool.h
- *
- *
- * @author Adlene Hicheur
- * @date 2007-06
- */
-
-class SpmInvTool : public GaudiTool, virtual public IAlignSolvTool {
-public:
- /// Standard constructor
- SpmInvTool( const std::string& type, const std::string& name, const IInterface* parent );
-
- virtual ~SpmInvTool(); ///< Destructor
-
- using IAlignSolvTool::compute; ///< avoids hiding the original function definitions
- // Call method to compute the solution, get symetric matrix and
- // vector in input and replaces them by inverse matrix and solution vector in output.
- bool compute( AlSymMat&, AlVec& ) const override;
-
-protected:
-private:
- // parameters
-
- int SolvSpmInv( AlSymMat& M, AlVec& B );
-
- void Precond( AlSymMat& M, AlVec& r, AlVec& c, bool equIn );
-};
-#endif // ALIGNSOLVTOOLS_ALIGNSOLVTOOL_H
diff --git a/Alignment/AlignSolvTools/src/gslSVDsolver.cpp b/Alignment/AlignSolvTools/src/gslSVDsolver.cpp
index 297c747bee14274becef082cf0219178e225f1ff..b15eed9aa69f75fc4f3dfe8cfbb350ec00812193 100755
--- a/Alignment/AlignSolvTools/src/gslSVDsolver.cpp
+++ b/Alignment/AlignSolvTools/src/gslSVDsolver.cpp
@@ -1,5 +1,5 @@
/*****************************************************************************\
-* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration *
+* (c) Copyright 2000-2019 CERN for the benefit of the LHCb Collaboration *
* *
* This software is distributed under the terms of the GNU General Public *
* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
@@ -12,12 +12,55 @@
#include <cmath>
-// local
-#include "gslSVDsolver.h"
+// Include files
+// from Gaudi
+#include "GaudiAlg/GaudiTool.h"
+
+// from GSL
+#include "gsl/gsl_linalg.h"
+
+// Interface
+#include "AlignmentInterfaces/IAlignSolvTool.h"
+
+/** @class gslSVDsolver gslSVDsolver.h
+ *
+ *
+ * @author Jan Amoraal
+ * @date 2007-07-24
+ */
+
+class gslSVDsolver : public GaudiTool, virtual public IAlignSolvTool {
+
+public:
+ /// Standard constructor
+ gslSVDsolver( const std::string& type, const std::string& name, const IInterface* parent );
+
+ virtual ~gslSVDsolver(); ///< Destructor
+
+ using IAlignSolvTool::compute; ///< avoids hiding the original function definitions
+ /// Solves Ax = b using gsl_linalg_SV_decomp (gsl_matrix * A, gsl_matrix * V, gsl_vector * S, gsl_vector * work)
+ StatusCode compute( AlSymMat& symMatrix, AlVec& vector, AlVec& evvector ) const override;
+
+ StatusCode compute( AlSymMat& symMatrix, AlVec& vector ) const override {
+ AlVec dummy( vector.size() );
+ return compute( symMatrix, vector, dummy );
+ }
+
+protected:
+private:
+ bool m_svdJacobi; ///< Use Jacobi SVD?
+ double m_svdEpsilon; ///< SVD threshold. The singular values are stored in gsl_vector S.
+ unsigned int m_nZero; ///< remove the n smallest eigenvalues
+ double m_eigenValueThreshold;
+ size_t m_numberOfPrintedEigenvalues;
+};
//-----------------------------------------------------------------------------
// Implementation file for class : gslSVDsolver
//
+// See also
+// http://home.thep.lu.se/~jari/documents/gsl-ref.html/Singular-Value-Decomposition.html
+//
// 2007-07-24 : Jan Amoraal
//-----------------------------------------------------------------------------
@@ -41,7 +84,7 @@ gslSVDsolver::gslSVDsolver( const std::string& type, const std::string& name, co
//=============================================================================
gslSVDsolver::~gslSVDsolver() {}
-bool gslSVDsolver::compute( AlSymMat& symMatrix, AlVec& vector, AlVec& evvector ) const {
+StatusCode gslSVDsolver::compute( AlSymMat& symMatrix, AlVec& vector, AlVec& evvector ) const {
size_t size = symMatrix.size();
@@ -52,8 +95,8 @@ bool gslSVDsolver::compute( AlSymMat& symMatrix, AlVec& vector, AlVec& evvector
/// Fill matrix A
for ( unsigned i = 0u; i < size; ++i ) {
for ( unsigned j = 0u; j < size; ++j ) {
- debug() << "Element (i,j) of AlSymMat A = " << symMatrix[i][j] << endmsg;
- gsl_matrix_set( matrixA, i, j, symMatrix[i][j] );
+ debug() << "Element (i,j) of AlSymMat A = " << symMatrix( i, j ) << endmsg;
+ gsl_matrix_set( matrixA, i, j, symMatrix( i, j ) );
debug() << "Element (i,j) of gsl_matrix A = " << gsl_matrix_get( matrixA, i, j ) << endmsg;
}
}
@@ -64,8 +107,6 @@ bool gslSVDsolver::compute( AlSymMat& symMatrix, AlVec& vector, AlVec& evvector
gsl_vector* vectorW = gsl_vector_alloc( size );
gsl_matrix* matrixU = matrixA;
- debug() << "==> Matrix A = " << ( *matrixA ) << endmsg;
-
debug() << "Factorising matrix A" << endmsg;
/// Factorise A into the SVD A = USV^T. Note matrix A is replaced with matrix U.
/// GSL returns 0 if successful
@@ -79,14 +120,10 @@ bool gslSVDsolver::compute( AlSymMat& symMatrix, AlVec& vector, AlVec& evvector
}
if ( factorised != 0 ) {
error() << "==> Couldn't factorise matrix" << endmsg;
- return false;
+ return StatusCode::FAILURE;
}
debug() << "Done factorising matrix A" << endmsg;
- debug() << "==> Matrix U = " << ( *matrixA ) << endmsg;
- debug() << "==> Vector S = " << ( *vectorS ) << endmsg;
- debug() << "==> Matrix V = " << ( *matrixV ) << endmsg;
-
{
// print the eigenvalues, but make sure to insert minus signs where needed
std::vector<double> eigenvalues;
@@ -139,11 +176,10 @@ bool gslSVDsolver::compute( AlSymMat& symMatrix, AlVec& vector, AlVec& evvector
}
info() << "Number of removed eigenvalues = " << numremoved << endmsg;
- debug() << "==> Regularised Vector S = " << ( *vectorS ) << endmsg;
// Replace symMatrix with its inverse (the covariance matrix)
for ( unsigned irow = 0; irow < size; ++irow )
- for ( unsigned int icol = 0; icol <= irow; ++icol ) symMatrix[irow][icol] = 0;
+ for ( unsigned int icol = 0; icol <= irow; ++icol ) symMatrix( irow, icol ) = 0;
for ( unsigned int k = 0; k < size; ++k ) {
double s = gsl_vector_get( vectorS, k );
if ( s > 0 ) {
@@ -152,7 +188,7 @@ bool gslSVDsolver::compute( AlSymMat& symMatrix, AlVec& vector, AlVec& evvector
// honestly: don't have a clue about the order of U and V
// here, but they should be identical for our symmetric
// systems
- symMatrix[irow][icol] += gsl_matrix_get( matrixV, irow, k ) * gsl_matrix_get( matrixU, icol, k ) / s;
+ symMatrix( irow, icol ) += gsl_matrix_get( matrixV, irow, k ) * gsl_matrix_get( matrixU, icol, k ) / s;
}
}
}
@@ -165,14 +201,12 @@ bool gslSVDsolver::compute( AlSymMat& symMatrix, AlVec& vector, AlVec& evvector
int solved = gsl_linalg_SV_solve( matrixA, matrixV, vectorS, vectorB, vectorX );
if ( solved != 0 ) {
error() << "==> Couldn't solve system Ax=b" << endmsg;
- return false;
+ return StatusCode::FAILURE;
}
/// Fill AlVec
for ( unsigned i = 0; i < size; ++i ) vector[i] = ( *gsl_vector_const_ptr( vectorX, i ) );
- debug() << "==> Vector x = " << ( *vectorX ) << endmsg;
-
/// free gsl vector and matrices
gsl_matrix_free( matrixA );
gsl_matrix_free( matrixV );
@@ -181,7 +215,7 @@ bool gslSVDsolver::compute( AlSymMat& symMatrix, AlVec& vector, AlVec& evvector
gsl_vector_free( vectorB );
gsl_vector_free( vectorX );
- return true;
+ return StatusCode::SUCCESS;
}
//=============================================================================
diff --git a/Alignment/AlignSolvTools/src/gslSVDsolver.h b/Alignment/AlignSolvTools/src/gslSVDsolver.h
deleted file mode 100755
index 469f562fe039c5761403907163514bb17920ba4a..0000000000000000000000000000000000000000
--- a/Alignment/AlignSolvTools/src/gslSVDsolver.h
+++ /dev/null
@@ -1,80 +0,0 @@
-/*****************************************************************************\
-* (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration *
-* *
-* This software is distributed under the terms of the GNU General Public *
-* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
-* *
-* In applying this licence, CERN does not waive the privileges and immunities *
-* granted to it by virtue of its status as an Intergovernmental Organization *
-* or submit itself to any jurisdiction. *
-\*****************************************************************************/
-// $Id: gslSVDsolver.h,v 1.9 2009-08-16 14:31:51 wouter Exp $
-#ifndef GSLSVDSOLVER_H
-#define GSLSVDSOLVER_H 1
-
-// Include files
-// from Gaudi
-#include "GaudiAlg/GaudiTool.h"
-
-// from GSL
-#include "gsl/gsl_linalg.h"
-
-// Interface
-#include "AlignmentInterfaces/IAlignSolvTool.h"
-
-/** @class gslSVDsolver gslSVDsolver.h
- *
- *
- * @author Jan Amoraal
- * @date 2007-07-24
- */
-
-class gslSVDsolver : public GaudiTool, virtual public IAlignSolvTool {
-
-public:
- /// Standard constructor
- gslSVDsolver( const std::string& type, const std::string& name, const IInterface* parent );
-
- virtual ~gslSVDsolver(); ///< Destructor
-
- using IAlignSolvTool::compute; ///< avoids hiding the original function definitions
- /// Solves Ax = b using gsl_linalg_SV_decomp (gsl_matrix * A, gsl_matrix * V, gsl_vector * S, gsl_vector * work)
- bool compute( AlSymMat& symMatrix, AlVec& vector, AlVec& evvector ) const override;
-
- bool compute( AlSymMat& symMatrix, AlVec& vector ) const override {
- AlVec dummy( vector.size() );
- return compute( symMatrix, vector, dummy );
- }
-
-protected:
-private:
- bool m_svdJacobi; ///< Use Jacobi SVD?
- double m_svdEpsilon; ///< SVD threshold. The singular values are stored in gsl_vector S.
- unsigned int m_nZero; ///< remove the n smallest eigenvalues
- double m_eigenValueThreshold;
- size_t m_numberOfPrintedEigenvalues;
-};
-
-inline std::ostream& operator<<( std::ostream& lhs, const gsl_vector& rhs ) {
- lhs << std::endl;
- lhs << "[ " << std::endl;
- for ( unsigned i = 0u; i < rhs.size; ++i ) {
- lhs << " " << ( *gsl_vector_const_ptr( &rhs, i ) ) << " ";
- lhs << std::endl;
- }
- lhs << " ]" << std::endl;
- return lhs;
-}
-
-inline std::ostream& operator<<( std::ostream& lhs, const gsl_matrix& rhs ) {
- lhs << std::endl;
- lhs << "[ " << std::endl;
- for ( unsigned i = 0u; i < rhs.size1; ++i ) {
- for ( unsigned j = 0u; j < rhs.size2; ++j ) { lhs << " " << ( *gsl_matrix_const_ptr( &rhs, i, j ) ) << " "; }
- lhs << std::endl;
- }
- lhs << " ]" << std::endl;
- return lhs;
-}
-
-#endif // GSLSVDSOLVER_H
diff --git a/Alignment/AlignTrTools/src/ATrackSelector.h b/Alignment/AlignTrTools/src/ATrackSelector.h
index 254969f4da9faef1646a4e7c791b2cb635ce1d5b..b5d3695063df9dbb3befe977596f0ecec56e4579 100755
--- a/Alignment/AlignTrTools/src/ATrackSelector.h
+++ b/Alignment/AlignTrTools/src/ATrackSelector.h
@@ -72,9 +72,8 @@ private:
private:
// Interfaces:
- IATrackSelectorTool* m_trackselector;
- ITrackExtrapolator* m_extrapolator;
- ITrackCaloMatch* m_trackenergy;
+ ITrackExtrapolator* m_extrapolator;
+ ITrackCaloMatch* m_trackenergy;
bool Unify( const LHCb::Track& );
bool uniformCut( int& );
diff --git a/Alignment/AlignTrTools/src/AlignSaveTuple.h b/Alignment/AlignTrTools/src/AlignSaveTuple.h
index 1bbfa7204e46050b8ede5e7e786fa24cb3d2b45a..6d63cc70759b17ab91567fc1e42877367fbf8531 100755
--- a/Alignment/AlignTrTools/src/AlignSaveTuple.h
+++ b/Alignment/AlignTrTools/src/AlignSaveTuple.h
@@ -44,9 +44,9 @@
//===========================================================================
// Forward declarations
//===========================================================================
-class ITrackExtrapolator;
+struct ITrackExtrapolator;
class IMagneticFieldSvc;
-class ITrackGhostClassification;
+struct ITrackGhostClassification;
class GhostTrackInfo;
//===========================================================================
diff --git a/Alignment/AlignTrTools/src/AlignTrackMonitor.h b/Alignment/AlignTrTools/src/AlignTrackMonitor.h
index 5994ba6e45cb4b24df56cac08ae2357910da3db3..a13b8d08927f28854ab609e433fdf2efe85927c2 100755
--- a/Alignment/AlignTrTools/src/AlignTrackMonitor.h
+++ b/Alignment/AlignTrTools/src/AlignTrackMonitor.h
@@ -43,7 +43,7 @@
//===========================================================================
// Forward declarations
//===========================================================================
-class ITrackExtrapolator;
+struct ITrackExtrapolator;
class IMagneticFieldSvc;
//===========================================================================
diff --git a/Alignment/AlignTrTools/src/DToKPiTwoProng.h b/Alignment/AlignTrTools/src/DToKPiTwoProng.h
index 9fdd8697fd8c86f018fdef5aafef4366ddf9a754..3116198aeb7b12b11b86edb16344b13ec56b0e53 100644
--- a/Alignment/AlignTrTools/src/DToKPiTwoProng.h
+++ b/Alignment/AlignTrTools/src/DToKPiTwoProng.h
@@ -24,10 +24,9 @@ namespace LHCb {
class ILHCbMagnetSvc;
class IHitExpectation;
-class ITrackVertexer;
-class ITrackExtrapolator;
-
-class ITrackFitter;
+struct ITrackVertexer;
+struct ITrackExtrapolator;
+struct ITrackFitter;
class DToKPiTwoProng final : public GaudiTupleAlg {
diff --git a/Alignment/AlignTrTools/src/MakeMuonTracks.h b/Alignment/AlignTrTools/src/MakeMuonTracks.h
index 956c671f6a73b3759ddfdc972b09604e190d7657..ccf3e90b1edf4004c73cca8c69da78975a4f5957 100644
--- a/Alignment/AlignTrTools/src/MakeMuonTracks.h
+++ b/Alignment/AlignTrTools/src/MakeMuonTracks.h
@@ -21,7 +21,7 @@
#include "MuonDet/DeMuonDetector.h"
#include "MuonInterfaces/IMuonTrackRec.h"
-class ITrackMomentumEstimate;
+struct ITrackMomentumEstimate;
/** @class MakeMuonTracks MakeMuonTracks.h
* Simple algorithm to make "LHCb" tracks from MuonNNet tracks
diff --git a/Alignment/AlignTrTools/src/MuonTrackSelector.cpp b/Alignment/AlignTrTools/src/MuonTrackSelector.cpp
index f0691e1600ec3ee49568f0f7f4ffbb5e649b7a96..35c563b78c5b2f86902e2e7008dda4c490cac9d8 100644
--- a/Alignment/AlignTrTools/src/MuonTrackSelector.cpp
+++ b/Alignment/AlignTrTools/src/MuonTrackSelector.cpp
@@ -31,8 +31,7 @@ MuonTrackSelector::MuonTrackSelector( const std::string& name, ISvcLocator* pSvc
, m_tracksInputContainer( "" )
, m_tracksOutputContainer( "" )
, m_trackType( "" )
- , m_trackSelectorName( "" )
- , m_trackSelector( 0 ) {
+ , m_trackSelectorName( "" ) {
declareProperty( "TracksInputContainer", m_tracksInputContainer = TrackLocation::Default );
declareProperty( "TracksOutputContainer", m_tracksOutputContainer = "Alignment/FilteredTracks" );
declareProperty( "TrackType", m_trackType = "Long" );
diff --git a/Alignment/AlignTrTools/src/MuonTrackSelector.h b/Alignment/AlignTrTools/src/MuonTrackSelector.h
index cc103f8e88de5dde657431af6f003395b5d39037..5d375fdcfcdd23e72cbe3801608312b514423df3 100644
--- a/Alignment/AlignTrTools/src/MuonTrackSelector.h
+++ b/Alignment/AlignTrTools/src/MuonTrackSelector.h
@@ -26,9 +26,6 @@
// from TrackEvent
#include "Event/Track.h"
-// from AlignmentInterfaces
-#include "TrackInterfaces/ITrackSelector.h"
-
/** @class TrackFilterAlg TrackFilterAlg.h
* A small algorithm that filters tracks of a certain type.
*
@@ -55,17 +52,16 @@ private:
void filterMuonTrack( LHCb::Track* track, LHCb::Tracks* outputContainer );
- std::string m_tracksInputContainer; ///< Tracks input container
- std::string m_tracksOutputContainer; ///< Tracks output container
- std::string m_trackType; ///< Type of tracks to filter
- std::string m_trackSelectorName; ///< Name of track selector for alignment
- ITrackSelector* m_trackSelector; ///< Pointer to track selector tool for alignment
- bool m_calo;
- bool m_noOverlap;
- double m_pcut;
- double m_muonChisquareCut;
- int m_nStation;
- int m_theR;
+ std::string m_tracksInputContainer; ///< Tracks input container
+ std::string m_tracksOutputContainer; ///< Tracks output container
+ std::string m_trackType; ///< Type of tracks to filter
+ std::string m_trackSelectorName; ///< Name of track selector for alignment
+ bool m_calo;
+ bool m_noOverlap;
+ double m_pcut;
+ double m_muonChisquareCut;
+ int m_nStation;
+ int m_theR;
};
/*
diff --git a/Alignment/AlignTrTools/src/OTMuonCosmicsMatching.h b/Alignment/AlignTrTools/src/OTMuonCosmicsMatching.h
index 4f5025e76212a45530ca49f9d12b804fc65a13d4..5e3525b4733736109f078104f794c1805d975d68 100755
--- a/Alignment/AlignTrTools/src/OTMuonCosmicsMatching.h
+++ b/Alignment/AlignTrTools/src/OTMuonCosmicsMatching.h
@@ -34,7 +34,7 @@
#include "TrackInterfaces/ITrackChi2Calculator.h"
namespace MatchingHelpers {
- class MatchedTrack;
+ struct MatchedTrack;
}
/** @class OTMuonCosmicsMatching OTMuonCosmicsMatching.h
diff --git a/Alignment/AlignTrTools/src/ParticleToTrackContainer.cpp b/Alignment/AlignTrTools/src/ParticleToTrackContainer.cpp
index b6e78656941bc5e6e6248582067119ad170f4c22..ff602300b5023449aba27f2f852828147bd9bdc3 100644
--- a/Alignment/AlignTrTools/src/ParticleToTrackContainer.cpp
+++ b/Alignment/AlignTrTools/src/ParticleToTrackContainer.cpp
@@ -40,8 +40,6 @@ public:
private:
std::string m_inputLocation;
std::string m_outputLocation;
- double m_minMass;
- double m_maxMass;
ToolHandle<ITrackSelector> m_selector;
};
diff --git a/Alignment/TAlignment/src/TStation.cpp b/Alignment/AlignTrTools/src/TStation.cpp
similarity index 100%
rename from Alignment/TAlignment/src/TStation.cpp
rename to Alignment/AlignTrTools/src/TStation.cpp
diff --git a/Alignment/TAlignment/src/TStation.h b/Alignment/AlignTrTools/src/TStation.h
similarity index 100%
rename from Alignment/TAlignment/src/TStation.h
rename to Alignment/AlignTrTools/src/TStation.h
diff --git a/Alignment/AlignTrTools/src/TrackMuonMatching.h b/Alignment/AlignTrTools/src/TrackMuonMatching.h
index 81f2aae2868ccdb6017f336b9439f05fc8f77f44..dacbbdece8c0f87147aae2f0b992272fa70635af 100644
--- a/Alignment/AlignTrTools/src/TrackMuonMatching.h
+++ b/Alignment/AlignTrTools/src/TrackMuonMatching.h
@@ -66,7 +66,6 @@ private:
std::string m_tracksOutputLocation;
double m_matchAtZ;
bool m_matchAtFirstMuonHit;
- bool m_matchAtOT;
double m_matchChi2Cut;
bool m_allCombinations;
DeMuonDetector* m_muonDet;
diff --git a/Alignment/AlignTrTools/src/TrackPVRefitter.cpp b/Alignment/AlignTrTools/src/TrackPVRefitter.cpp
index 883cbe17310da7bf76b1a10cf0dff34c95cc4f81..032f26b5f56edc8371a86ca20664cf4148678c8c 100644
--- a/Alignment/AlignTrTools/src/TrackPVRefitter.cpp
+++ b/Alignment/AlignTrTools/src/TrackPVRefitter.cpp
@@ -35,7 +35,6 @@ public:
private:
std::string m_trackContainerName;
std::string m_pvContainerName;
- double m_maxLongTrackChisqPerDof;
ToolHandle<ITrackVertexer> m_vertexer;
ToolHandle<ITrackFitter> m_trackfitter;
diff --git a/Alignment/AlignTrTools/src/TrackParticleRefitter.cpp b/Alignment/AlignTrTools/src/TrackParticleRefitter.cpp
index 3824887ca58ae818cf6055f2d1c99ee2c450d59f..fd3e7536602a49ccbb1e5896069d94a21e9a604f 100644
--- a/Alignment/AlignTrTools/src/TrackParticleRefitter.cpp
+++ b/Alignment/AlignTrTools/src/TrackParticleRefitter.cpp
@@ -19,143 +19,81 @@
#include "Event/KalmanFitResult.h"
#include "Event/Particle.h"
#include "Event/Track.h"
-#include "Event/Vertex.h"
#include "GaudiAlg/GaudiAlgorithm.h"
-#include "GaudiKernel/SharedObjectsContainer.h"
#include "GaudiKernel/ToolHandle.h"
-#include "Kernel/IParticlePropertySvc.h"
-#include "Kernel/ParticleID.h"
-#include "Kernel/ParticleProperty.h"
+#include "Kernel/IVertexFit.h"
#include "LHCbMath/LorentzVectorWithError.h"
#include "TrackInterfaces/ITrackFitter.h"
#include "TrackInterfaces/ITrackStateProvider.h"
-#include "TrackKernel/TrackStateVertex.h"
-#include <boost/foreach.hpp>
+
#include <string>
class TrackParticleRefitter : public GaudiAlgorithm {
public:
// Constructors and destructor
- TrackParticleRefitter( const std::string& name, ISvcLocator* pSvcLocator );
- virtual ~TrackParticleRefitter();
-
- StatusCode initialize() override;
- StatusCode finalize() override;
+ using GaudiAlgorithm::GaudiAlgorithm;
StatusCode execute() override;
private:
- std::string m_inputLocation;
- std::string m_particleName;
- int m_pid;
- ToolHandle<ITrackFitter> m_trackfitter;
- ToolHandle<ITrackStateProvider> m_stateprovider;
+ StatusCode refittree( const LHCb::Particle& p, std::set<const LHCb::Particle*>& refittedparticles );
+
+private:
+ Gaudi::Property<std::string> m_inputLocation{this, "ParticleLocation"};
+ Gaudi::Property<bool> m_forceTrackRefit{this, "ForceTrackRefit", false};
+ ToolHandle<IVertexFit> m_vertexfitter{this, "VeretxFitter", "ParticleVertexFitter"};
+ ToolHandle<ITrackFitter> m_trackfitter{this, "TrackFitter", "TrackMasterFitter"};
+ PublicToolHandle<ITrackStateProvider> m_stateprovider{this, "TrackStateProvider", "TrackStateProvider"};
};
DECLARE_COMPONENT( TrackParticleRefitter )
-TrackParticleRefitter::TrackParticleRefitter( const std::string& name, ISvcLocator* pSvcLocator )
- : GaudiAlgorithm( name, pSvcLocator )
- , m_pid( 0 )
- , m_trackfitter( "TrackMasterFitter", this )
- , m_stateprovider( "TrackStateProvider" ) {
- // constructor
- declareProperty( "ParticleName", m_particleName = "" );
- declareProperty( "ParticleLocation", m_inputLocation );
- declareProperty( "TrackFitter", m_trackfitter );
-}
-
-StatusCode TrackParticleRefitter::initialize() {
- StatusCode sc = GaudiAlgorithm::initialize();
- if ( sc.isSuccess() ) sc = m_trackfitter.retrieve();
- if ( sc.isSuccess() ) sc = m_stateprovider.retrieve();
- if ( sc.isSuccess() && !m_particleName.empty() ) {
- LHCb::IParticlePropertySvc* propertysvc = svc<LHCb::IParticlePropertySvc>( "LHCb::ParticlePropertySvc", true );
- const LHCb::ParticleProperty* prop = propertysvc->find( m_particleName );
- if ( prop != 0 )
- m_pid = prop->particleID().pid();
- else {
- error() << "Cannot find particle property for \'" << m_particleName << "\'" << endmsg;
- sc = StatusCode::FAILURE;
- }
- }
- return sc;
-}
-
-StatusCode TrackParticleRefitter::finalize() {
- m_trackfitter.release().ignore();
- m_stateprovider.release().ignore();
- return GaudiAlgorithm::finalize();
-}
-
-TrackParticleRefitter::~TrackParticleRefitter() {
- // destructor
-}
-
-namespace {
- void addTracks( const LHCb::Particle& p, std::vector<const LHCb::Track*>& tracks, std::vector<double>& masshypos ) {
+StatusCode TrackParticleRefitter::refittree( const LHCb::Particle& p,
+ std::set<const LHCb::Particle*>& refittedparticles ) {
+ StatusCode sc = StatusCode::SUCCESS;
+ if ( refittedparticles.find( &p ) == refittedparticles.end() ) {
+ // if we are a simple track, then refit the track
if ( p.proto() && p.proto()->track() ) {
- tracks.push_back( p.proto()->track() );
- masshypos.push_back( p.momentum().M() );
+ auto tr = p.proto()->track();
+ if ( m_forceTrackRefit || !dynamic_cast<const LHCb::KalmanFitResult*>( tr->fitResult() ) ) {
+ StatusCode thissc = m_trackfitter->operator()( const_cast<LHCb::Track&>( *tr ) );
+ if ( !sc.isSuccess() ) {
+ warning() << "problem fitting track" << endmsg;
+ sc = thissc;
+ }
+ m_stateprovider->clearCache( *tr ); // better safe than sorry
+ }
} else {
- BOOST_FOREACH ( const LHCb::Particle* dau, p.daughters() )
- addTracks( *dau, tracks, masshypos );
+ // make sure to refit all daughters
+ for ( const LHCb::Particle* dau : p.daughters() ) {
+ StatusCode thissc = refittree( *dau, refittedparticles );
+ if ( !thissc.isSuccess() ) sc = thissc;
+ }
+ // if we have an endvertex, do a vertex fit
+ if ( p.endVertex() ) {
+ LHCb::Particle& ncp = const_cast<LHCb::Particle&>( p );
+ StatusCode thissc = m_vertexfitter->reFit( ncp );
+ // never quite understood why we need to set this separately
+ ncp.setMeasuredMass( ncp.momentum().M() );
+ ncp.setMeasuredMassErr( Gaudi::Math::LorentzVectorWithError( ncp.momentum(), ncp.momCovMatrix() ).m().error() );
+ if ( !thissc.isSuccess() ) sc = thissc;
+ }
}
+ // add this particle to the list of particles that we already treated
+ refittedparticles.insert( &p );
}
-} // namespace
+ return sc;
+}
StatusCode TrackParticleRefitter::execute() {
- LHCb::Particle::Range inputparticles = get<LHCb::Particle::Range>( m_inputLocation );
- BOOST_FOREACH ( const LHCb::Particle* p, inputparticles ) {
- double z = p->referencePoint().z();
-
- std::vector<const LHCb::Track*> tracks;
- std::vector<double> masshypos;
- addTracks( *p, tracks, masshypos );
- if ( tracks.size() < 2 ) {
- warning() << "Not enough tracks! .. skipping particle" << p->daughters().size() << " " << tracks.size() << endmsg;
-
- } else {
- std::vector<const LHCb::State*> states;
- BOOST_FOREACH ( const LHCb::Track* tr, tracks ) {
- if ( !dynamic_cast<const LHCb::KalmanFitResult*>( tr->fitResult() ) ) {
- StatusCode sc = m_trackfitter->operator()( const_cast<LHCb::Track&>( *tr ) );
- if ( !sc.isSuccess() ) warning() << "problem fitting track" << endmsg;
- }
- LHCb::State* state = new LHCb::State();
- StatusCode sc = m_stateprovider->stateFromTrajectory( *state, *tr, z );
- if ( !sc.isSuccess() ) warning() << "problem getting state from stateprovider" << endmsg;
- states.push_back( state );
- }
-
- LHCb::TrackStateVertex vertex( states );
- vertex.fit();
- // now copy everything
- LHCb::Particle* ncp = const_cast<LHCb::Particle*>( p );
- if ( m_pid ) ncp->setParticleID( LHCb::ParticleID( m_pid ) );
- ncp->setMomentum( vertex.p4( masshypos ) );
- ncp->setReferencePoint( vertex.position() );
- Gaudi::SymMatrix7x7 cov7x7 = vertex.covMatrix7x7( masshypos );
- ncp->setMomCovMatrix( cov7x7.Sub<Gaudi::SymMatrix4x4>( 3, 3 ) );
- ncp->setPosCovMatrix( cov7x7.Sub<Gaudi::SymMatrix3x3>( 0, 0 ) );
- ncp->setPosMomCovMatrix( cov7x7.Sub<Gaudi::Matrix4x3>( 3, 0 ) );
- ncp->setMeasuredMass( ncp->momentum().M() );
- // ncp->setMeasuredMassErr( vertex.massErr(masshypos) ) ;
- ncp->setMeasuredMassErr(
- Gaudi::Math::LorentzVectorWithError( ncp->momentum(), ncp->momCovMatrix() ).m().error() );
- LHCb::Vertex* endvertex = ncp->endVertex();
- if ( endvertex ) {
- endvertex->setPosition( vertex.position() );
- endvertex->setChi2AndDoF( vertex.chi2(), vertex.nDoF() );
- endvertex->setCovMatrix( vertex.covMatrix() );
- }
-
- // don't forget to delete all states !
- BOOST_FOREACH ( const LHCb::State* s, states )
- delete s;
- }
+ LHCb::Particle::Range inputparticles = get<LHCb::Particle::Range>( m_inputLocation );
+ // we keep track of which particles were treated, so we don't do them more than once.
+ std::set<const LHCb::Particle*> refittedparticles;
+ // now do the vertex fits
+ for ( const LHCb::Particle* p : inputparticles ) {
+ StatusCode sc = refittree( *p, refittedparticles );
+ if ( !sc.isSuccess() ) warning() << "Vertex fit failed" << endmsg;
}
-
return StatusCode::SUCCESS;
}
diff --git a/Alignment/AlignmentInterfaces/AlignmentInterfaces/IAlignSolvTool.h b/Alignment/AlignmentInterfaces/AlignmentInterfaces/IAlignSolvTool.h
index f2dab86c3d357542cec58926eed9eae54aa685bd..c8ddf5a390582b3df5e6729feda4d6ee64f2c219 100755
--- a/Alignment/AlignmentInterfaces/AlignmentInterfaces/IAlignSolvTool.h
+++ b/Alignment/AlignmentInterfaces/AlignmentInterfaces/IAlignSolvTool.h
@@ -8,13 +8,11 @@
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
-// $Id: IAlignSolvTool.h,v 1.7 2009-07-06 14:06:44 wouter Exp $
-#ifndef ALIGNSOLVTOOLS_IALIGNSOLVTOOL_H
-#define ALIGNSOLVTOOLS_IALIGNSOLVTOOL_H 1
+#pragma once
// Include files
// from STL
-#include <string>
+//#include <string>
// from Gaudi
#include "GaudiKernel/IAlgTool.h"
@@ -47,11 +45,11 @@ public:
};
// Solves the system Ax=b for x. Returns A=A^{-1} and b=x.
- virtual bool compute( AlSymMat& A, AlVec& b ) const = 0;
+ virtual StatusCode compute( AlSymMat& A, AlVec& b ) const = 0;
// Solves the system Ax=b for x. Returns A=A^{-1}, b=x and S =
// eigenvalues of A.
- virtual bool compute( AlSymMat& A, AlVec& b, AlVec& /*S*/ ) const { return compute( A, b ); }
+ virtual StatusCode compute( AlSymMat& A, AlVec& b, AlVec& /*S*/ ) const { return compute( A, b ); }
// Solves the system Ax=b for x. Return A=A^{-1} and
// b=x. Furthermore, returns an object with some characteristics of
@@ -60,4 +58,3 @@ public:
return compute( A, b ) ? StatusCode::SUCCESS : StatusCode::FAILURE;
}
};
-#endif // ALIGNSOLVTOOLS_IALIGNSOLVTOOL_H
diff --git a/Alignment/Escher/scripts/solvertest.py b/Alignment/Escher/scripts/solvertest.py
new file mode 100644
index 0000000000000000000000000000000000000000..79901021320e0a9e035bad248c793c8e6c0e7deb
--- /dev/null
+++ b/Alignment/Escher/scripts/solvertest.py
@@ -0,0 +1,171 @@
+#!/usr/bin/env python
+###############################################################################
+# (c) Copyright 2000-2018 CERN for the benefit of the LHCb Collaboration #
+# #
+# This software is distributed under the terms of the GNU General Public #
+# Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". #
+# #
+# In applying this licence, CERN does not waive the privileges and immunities #
+# granted to it by virtue of its status as an Intergovernmental Organization #
+# or submit itself to any jurisdiction. #
+###############################################################################
+
+# first parse all options
+from optparse import OptionParser
+parser = OptionParser(usage="%prog [options] <opts_file> ...")
+parser.add_option(
+ "-n",
+ "--numiter",
+ type="int",
+ dest="numiter",
+ help="number of iterations",
+ default=1)
+parser.add_option(
+ "-e",
+ "--numevents",
+ type="int",
+ dest="numevents",
+ help="number of events",
+ default=-1)
+parser.add_option(
+ "-d",
+ "--aligndb",
+ action='append',
+ dest="aligndb",
+ help="path to file with CALIBOFF/LHCBCOND/SIMCOND database layer")
+parser.add_option(
+ "--dddb",
+ action='append',
+ dest="dddb",
+ help="path to file with DDDB database layer")
+parser.add_option(
+ "--skipprintsequence",
+ action="store_true",
+ help="skip printing the sequence")
+parser.add_option(
+ "-l",
+ "--lhcbcond",
+ action="store_true",
+ dest="lhcbcondtag",
+ help="activate if constants in lhcbcondDB",
+ default=False)
+(opts, args) = parser.parse_args()
+
+# Prepare the "configuration script" to parse (like this it is easier than
+# having a list with files and python commands, with an if statements that
+# decides to do importOptions or exec)
+options = ["importOptions(%r)" % f for f in args]
+
+# The option lines are inserted into the list of commands using their
+# position on the command line
+#optlines = list(opts.options)
+#optlines.reverse() # this allows to avoid to have to care about corrections of the positions
+#for pos, l in optlines:
+# options.insert(pos,l)
+
+# "execute" the configuration script generated (if any)
+from Gaudi.Configuration import logging
+if options:
+ g = {}
+ l = {}
+ exec "from Gaudi.Configuration import *" in g, l
+ for o in options:
+ logging.debug(o)
+ exec o in g, l
+
+#options = [ "importOptions(%r)" % f for f in args ]
+
+# make sure that the algorithms know how many iterations are coming
+from Configurables import TAlignment
+TAlignment().UpdateInFinalize = False
+
+# set the database layer
+if opts.aligndb:
+ counter = 1
+ for db in opts.aligndb:
+ from Configurables import CondDB
+ if ".db" in db:
+ # this is an sqlite file
+ from Configurables import CondDBAccessSvc
+ alignCond = CondDBAccessSvc('AlignCond' + str(counter))
+ alignCond.ConnectionString = 'sqlite_file:' + db + condtag
+ CondDB().addLayer(alignCond)
+ counter += 1
+ else:
+ # this is a git layer
+ CondDB().addLayer(db)
+ print "adding layer: ", db
+
+if opts.dddb:
+ counter = 1
+ for db in opts.dddb:
+ from Configurables import (CondDB, CondDBAccessSvc)
+ alignCond = CondDBAccessSvc('AlignDDDB' + str(counter))
+ alignCond.ConnectionString = 'sqlite_file:' + db + '/DDDB'
+ CondDB().addLayer(alignCond)
+ counter += 1
+
+## configure all solver tools (before creating the appmgr)
+solvers = [
+ "GslDiagSolvTool", "EigenDiagSolvTool", "SpmInvTool", "gslSVDsolver",
+ "CLHEPSolver", "SparseSolver"
+]
+updatetools = []
+for solvername in solvers:
+ # first configure it, then retrieve it via the toolsvc. does that make sense?
+ from Configurables import Al__AlignUpdateTool as AlignUpdateTool
+ updatetoolname = "UpdateTool" + solvername
+ dir(AlignUpdateTool())
+ AlignUpdateTool().clone(updatetoolname)
+
+ updatetool = AlignUpdateTool(
+ updatetoolname,
+ MatrixSolverTool=solvername,
+ SkipGeometryUpdate=True,
+ LogFile="alignlog_" + solvername + ".txt",
+ OutputLevel=3)
+ # for now we do it like this
+ from Configurables import SurveyConstraints, TAlignment
+ SurveyConstraints().configureTool(updatetool.SurveyConstraintTool)
+ updatetool.LagrangeConstraintTool.Constraints = TAlignment().getProp(
+ "Constraints")
+ updatetool.LagrangeConstraintTool.UseWeightedAverage = TAlignment(
+ ).getProp("UseWeightedAverageConstraint")
+ updatetools.append(updatetoolname)
+
+## Instantiate application manager
+from GaudiPython.Bindings import AppMgr
+appMgr = AppMgr()
+
+## print the sequence
+if not opts.skipprintsequence:
+ from Escher.Utils import printsequence
+ printsequence(appMgr)
+
+evtSel = appMgr.evtSel()
+evtSel.OutputLevel = 1
+mainSeq = appMgr.algorithm("EscherSequencer")
+
+print evtSel.Input
+
+# event loop
+appMgr.run(opts.numevents)
+
+# get the derivatives
+import copy
+det = appMgr.detsvc()
+alignderivatives = det['AlignDerivativeData']
+print "alignderivatives:", alignderivatives
+for updatetoolname in updatetools:
+ print ["*" * 40]
+ print "UpdateTool: ", updatetoolname
+ updatetool = appMgr.toolsvc().create(
+ typ="Al::AlignUpdateTool",
+ name=updatetoolname,
+ interface="Al::IAlignUpdateTool")
+
+ equations = copy.deepcopy(alignderivatives.equations())
+ updatetool.process(equations, 0, 1)
+
+#exit the appmgr for finalize
+#appMgr.exit()
diff --git a/Alignment/TAlignment/python/TAlignment/Configuration.py b/Alignment/TAlignment/python/TAlignment/Configuration.py
index ff0fc60cc8dcaf4c1651418892058d0fd52fb7ab..ac6995c24bc8d9db91a05840cb8684be43a8078e 100644
--- a/Alignment/TAlignment/python/TAlignment/Configuration.py
+++ b/Alignment/TAlignment/python/TAlignment/Configuration.py
@@ -84,7 +84,7 @@ class TAlignment(LHCbConfigurableUser):
True # Pre-conditioning
,
"SolvTool":
- "DiagSolvTool" # Solver to use
+ "EigenDiagSolvTool" # Solver to use
,
"EigenValueThreshold":
-1 # Eigenvalue threshold for cutting out weak modes
@@ -396,7 +396,7 @@ class TAlignment(LHCbConfigurableUser):
from Configurables import (
AlignAlgorithm, GetElementsToBeAligned, gslSVDsolver,
- CLHEPSolver, SparseSolver, DiagSolvTool,
+ CLHEPSolver, SparseSolver, EigenDiagSolvTool, GslDiagSolvTool,
Al__AlignConstraintTool, Al__AlignUpdateTool,
Al__TrackResidualTool, WriteMultiAlignmentConditionsTool)
@@ -457,7 +457,9 @@ class TAlignment(LHCbConfigurableUser):
# and these too
gslSVDsolver().EigenValueThreshold = self.getProp(
"EigenValueThreshold")
- DiagSolvTool().EigenValueThreshold = self.getProp(
+ EigenDiagSolvTool().EigenValueThreshold = self.getProp(
+ "EigenValueThreshold")
+ GslDiagSolvTool().EigenValueThreshold = self.getProp(
"EigenValueThreshold")
alignSequencer.Members.append(alignAlg)
diff --git a/Alignment/TAlignment/src/AlElementHistos.h b/Alignment/TAlignment/src/AlElementHistos.h
index d4555cb12ce99d0f498fb69d0ce3428f5018c335..a9ef61a0f3b1c260c8a58387b77713ab1a565205 100755
--- a/Alignment/TAlignment/src/AlElementHistos.h
+++ b/Alignment/TAlignment/src/AlElementHistos.h
@@ -18,7 +18,7 @@ namespace AIDA {
class IHistogram2D;
} // namespace AIDA
class AlignmentElement;
-class GaudiHistoAlg;
+struct GaudiHistoAlg;
struct AlElementHistos {
AlElementHistos( GaudiHistoAlg& parent, const AlignmentElement& elem, std::size_t niter );
diff --git a/Alignment/TAlignment/src/AlParameters.cpp b/Alignment/TAlignment/src/AlParameters.cpp
index 1c7b2bb33e2e0e2a2a98217a24d8b5be7cbdc109..316af10b77b88f3654c802d80fdfa4d59ed3c3ef 100755
--- a/Alignment/TAlignment/src/AlParameters.cpp
+++ b/Alignment/TAlignment/src/AlParameters.cpp
@@ -28,8 +28,8 @@ AlParameters::AlParameters( DofMask mask ) : m_mask( mask ), m_parameters( dim()
AlParameters::AlParameters( const Vector& parameters, const Covariance& covariance, DofMask mask, size_t offset )
: m_mask( mask ), m_parameters( dim() ), m_covariance( dim() ), m_weightmatrix( dim() ) {
for ( unsigned int i = 0u; i < dim(); ++i ) {
- m_parameters[i] = parameters[i + offset];
- for ( unsigned int j = 0u; j <= i; ++j ) m_covariance[i][j] = covariance[offset + i][offset + j];
+ m_parameters( i ) = parameters( i + offset );
+ for ( unsigned int j = 0u; j <= i; ++j ) m_covariance( i, j ) = covariance( offset + i, offset + j );
}
}
@@ -37,10 +37,10 @@ AlParameters::AlParameters( const Vector& parameters, const Covariance& covarian
DofMask mask, size_t offset )
: m_mask( mask ), m_parameters( dim() ), m_covariance( dim() ), m_weightmatrix( dim() ) {
for ( unsigned int i = 0u; i < dim(); ++i ) {
- m_parameters[i] = parameters[i + offset];
+ m_parameters( i ) = parameters( i + offset );
for ( unsigned int j = 0u; j <= i; ++j ) {
- m_covariance[i][j] = covariance[offset + i][offset + j];
- m_weightmatrix[i][j] = weightmatrix[offset + i][offset + j];
+ m_covariance( i, j ) = covariance( offset + i, offset + j );
+ m_weightmatrix( i, j ) = weightmatrix( offset + i, offset + j );
}
}
}
@@ -48,29 +48,29 @@ AlParameters::AlParameters( const Vector& parameters, const Covariance& covarian
AlParameters::AlParameters( const TransformParameters& parameters, const TransformCovariance& covariance )
: m_mask( NumPars ), m_parameters( 6 ), m_covariance( 6 ), m_weightmatrix( 6 ) {
for ( unsigned int i = 0u; i < dim(); ++i ) {
- m_parameters[i] = parameters[i];
- for ( unsigned int j = 0u; j <= i; ++j ) m_covariance[i][j] = covariance[i][j];
+ m_parameters( i ) = parameters( i );
+ for ( unsigned int j = 0u; j <= i; ++j ) m_covariance( i, j ) = covariance( i, j );
}
}
AlParameters::AlParameters( const double parameters[6], DofMask mask )
: m_mask( mask ), m_parameters( dim() ), m_covariance( dim() ) {
- for ( unsigned int i = 0u; i < dim(); ++i ) m_parameters[i] = parameters[mask.parIndex( i )];
+ for ( unsigned int i = 0u; i < dim(); ++i ) m_parameters( i ) = parameters[mask.parIndex( i )];
}
AlParameters::AlParameters( const TransformParameters& parameters, DofMask mask )
: m_mask( mask ), m_parameters( dim() ), m_covariance( dim() ) {
- for ( unsigned int i = 0u; i < dim(); ++i ) m_parameters[i] = parameters[mask.parIndex( i )];
+ for ( unsigned int i = 0u; i < dim(); ++i ) m_parameters( i ) = parameters( mask.parIndex( i ) );
}
AlParameters::AlParameters( const ROOT::Math::Transform3D& transform, DofMask mask )
: m_mask( mask ), m_parameters( dim() ), m_covariance( dim() ) {
TransformParameters parameters = transformParameters( transform );
- for ( unsigned int i = 0u; i < dim(); ++i ) m_parameters[i] = parameters[mask.parIndex( i )];
+ for ( unsigned int i = 0u; i < dim(); ++i ) m_parameters( i ) = parameters( mask.parIndex( i ) );
}
void AlParameters::setParameters( const TransformParameters& parameters ) {
- for ( unsigned int i = 0u; i < dim(); ++i ) m_parameters[i] = parameters[m_mask.parIndex( i )];
+ for ( unsigned int i = 0u; i < dim(); ++i ) m_parameters( i ) = parameters( m_mask.parIndex( i ) );
}
std::string AlParameters::parName( int parindex ) {
@@ -106,30 +106,30 @@ double AlParameters::signedSqrt( double root ) {
std::ostream& AlParameters::fillStream( std::ostream& os ) const {
for ( unsigned int iactive = 0u; iactive < dim(); ++iactive ) {
- os << std::setw( 6 ) << parName( m_mask.parIndex( iactive ) ) << ": " << std::setw( 12 ) << m_parameters[iactive]
- << " +/- " << std::setw( 12 ) << signedSqrt( m_covariance[iactive][iactive] ) << std::endl;
+ os << std::setw( 6 ) << parName( m_mask.parIndex( iactive ) ) << ": " << std::setw( 12 ) << m_parameters( iactive )
+ << " +/- " << std::setw( 12 ) << signedSqrt( m_covariance( iactive, iactive ) ) << std::endl;
}
return os;
}
std::vector<double> AlParameters::translation() const {
std::vector<double> rc( 3, 0.0 );
- rc[Tx] = m_mask.isActive( Tx ) ? m_parameters[m_mask.activeParIndex( Tx )] : 0;
- rc[Ty] = m_mask.isActive( Ty ) ? m_parameters[m_mask.activeParIndex( Ty )] : 0;
- rc[Tz] = m_mask.isActive( Tz ) ? m_parameters[m_mask.activeParIndex( Tz )] : 0;
+ rc[Tx] = m_mask.isActive( Tx ) ? m_parameters( m_mask.activeParIndex( Tx ) ) : 0;
+ rc[Ty] = m_mask.isActive( Ty ) ? m_parameters( m_mask.activeParIndex( Ty ) ) : 0;
+ rc[Tz] = m_mask.isActive( Tz ) ? m_parameters( m_mask.activeParIndex( Tz ) ) : 0;
return rc;
}
std::vector<double> AlParameters::errTranslation() const {
std::vector<double> errT( 3, 0.0 );
errT[Tx] = m_mask.isActive( Tx )
- ? signedSqrt( m_covariance[m_mask.activeParIndex( Tx )][m_mask.activeParIndex( Tx )] )
+ ? signedSqrt( m_covariance( m_mask.activeParIndex( Tx ), m_mask.activeParIndex( Tx ) ) )
: 0.0;
errT[Ty] = m_mask.isActive( Ty )
- ? signedSqrt( m_covariance[m_mask.activeParIndex( Ty )][m_mask.activeParIndex( Ty )] )
+ ? signedSqrt( m_covariance( m_mask.activeParIndex( Ty ), m_mask.activeParIndex( Ty ) ) )
: 0.0;
errT[Tz] = m_mask.isActive( Tz )
- ? signedSqrt( m_covariance[m_mask.activeParIndex( Tz )][m_mask.activeParIndex( Tz )] )
+ ? signedSqrt( m_covariance( m_mask.activeParIndex( Tz ), m_mask.activeParIndex( Tz ) ) )
: 0.0;
return errT;
}
@@ -145,13 +145,13 @@ std::vector<double> AlParameters::rotation() const {
std::vector<double> AlParameters::errRotation() const {
std::vector<double> errR( 3, 0.0 );
errR[Rx - 3] = m_mask.isActive( Rx )
- ? signedSqrt( m_covariance[m_mask.activeParIndex( Rx )][m_mask.activeParIndex( Rx )] )
+ ? signedSqrt( m_covariance( m_mask.activeParIndex( Rx ), m_mask.activeParIndex( Rx ) ) )
: 0.0;
errR[Ry - 3] = m_mask.isActive( Ry )
- ? signedSqrt( m_covariance[m_mask.activeParIndex( Ry )][m_mask.activeParIndex( Ry )] )
+ ? signedSqrt( m_covariance( m_mask.activeParIndex( Ry ), m_mask.activeParIndex( Ry ) ) )
: 0.0;
errR[Rz - 3] = m_mask.isActive( Rz )
- ? signedSqrt( m_covariance[m_mask.activeParIndex( Rz )][m_mask.activeParIndex( Rz )] )
+ ? signedSqrt( m_covariance( m_mask.activeParIndex( Rz ), m_mask.activeParIndex( Rz ) ) )
: 0.0;
return errR;
}
@@ -322,7 +322,7 @@ AlParameters::Matrix6x6 AlParameters::jacobianNumeric( const ROOT::Math::Transfo
}
double AlParameters::globalCorrelationCoefficient( int iactive ) const {
- double varweight = m_covariance[iactive][iactive] * m_weightmatrix[iactive][iactive];
+ double varweight = m_covariance( iactive, iactive ) * m_weightmatrix( iactive, iactive );
return varweight > 1 ? std::sqrt( 1 - 1 / varweight ) : 0;
}
@@ -330,6 +330,6 @@ double AlParameters::measurementCoordinateSigma( double weightR ) const {
double sum( 0 );
for ( size_t iactive = 0; iactive < dim(); ++iactive )
for ( size_t jactive = 0; jactive < dim(); ++jactive )
- sum += m_covariance[iactive][jactive] * m_weightmatrix[iactive][jactive];
+ sum += m_covariance( iactive, jactive ) * m_weightmatrix( iactive, jactive );
return ( sum / weightR ) > 0 ? std::sqrt( sum / weightR ) : 0;
}
diff --git a/Alignment/TAlignment/src/AlignAlgorithm.cpp b/Alignment/TAlignment/src/AlignAlgorithm.cpp
index 98256ce32fd7ebfcf5f426c4af391c0611afd8df..e38e6e76ab7a3442a010f9aee304d8caefeb9ebc 100755
--- a/Alignment/TAlignment/src/AlignAlgorithm.cpp
+++ b/Alignment/TAlignment/src/AlignAlgorithm.cpp
@@ -96,53 +96,13 @@ public:
};
AlignAlgorithm::AlignAlgorithm( const std::string& name, ISvcLocator* pSvcLocator )
- : GaudiHistoAlg( name, pSvcLocator )
- , m_iteration( 0u )
- , m_nIterations( 0u )
- , m_nTracks( 0u )
- , m_covFailure( 0u )
- , m_align( 0 )
- , m_projSelector( 0 )
- , m_trackresidualtool( "Al::TrackResidualTool" )
- , m_vertexresidualtool( "Al::VertexResidualTool" )
- , m_updatetool( "Al::AlignUpdateTool" )
- , m_vertextrackselector( "TrackSelector", this )
- , m_xmlWritersTool( "WriteMultiAlignmentConditionsTool" )
- , m_fireRunChange( false )
- , m_equations( 0 )
- , m_resetHistos( false ) {
+ : GaudiHistoAlg( name, pSvcLocator ) {
m_HistoUpdater = new HistoUpdater();
- declareProperty( "NumberOfIterations", m_nIterations = 10u );
- declareProperty( "TracksLocation", m_tracksLocation = TrackLocation::Default );
- declareProperty( "VertexLocation", m_vertexLocation = "" );
- declareProperty( "DimuonLocation", m_dimuonLocation = "" );
- declareProperty( "ParticleLocation", m_particleLocation = "" );
- declareProperty( "ProjectorSelector", m_projSelectorName = "TrackProjectorSelector" );
- declareProperty( "UseCorrelations", m_correlation = true );
- declareProperty( "UpdateInFinalize", m_updateInFinalize = false );
- declareProperty( "OutputDataFile", m_outputDataFileName = "" );
- declareProperty( "InputDataFiles", m_inputDataFileNames );
- declareProperty( "Chi2Outlier", m_chi2Outlier = 10000 );
- declareProperty( "UpdateTool", m_updatetool );
- declareProperty( "TrackResidualTool", m_trackresidualtool );
- declareProperty( "VertexResidualTool", m_vertexresidualtool );
- declareProperty( "XmlWritersTool", m_xmlWritersTool );
- declareProperty( "MaxTracksPerVertex", m_maxTracksPerVertex = 8 );
- declareProperty( "MinTracksPerVertex", m_minTracksPerVertex = 2 );
-
- declareProperty( "AlignSummaryDataSvc", m_alignSummaryDataSvc = "DetectorDataSvc" );
- //"HistogramDataSvc" ) ;
- declareProperty( "AlignSummaryLocation", m_alignSummaryLocation = "AlignDerivativeData" );
- declareProperty( "FillHistos", m_fillHistos = false );
- declareProperty( "ForcedInitialTime", m_forcedInitTime = 0 );
- declareProperty( "OnlineMode", m_Online = false );
+
// declareProperty("XmlWriters",m_xmlWriterNames) ;
// For Online Running
declareProperty( "PartitionName", m_HistoUpdater->m_PartitionName = "LHCbA" );
declareProperty( "ReferenceFile", m_HistoUpdater->m_RefFileName = "" );
- declareProperty( "RunList", m_RunList );
-
- m_runnr = 0;
}
AlignAlgorithm::~AlignAlgorithm() { delete m_HistoUpdater; }
@@ -162,23 +122,11 @@ StatusCode AlignAlgorithm::initialize() {
incSvc->addListener( this, "DAQ_PAUSE" );
/// Get tool to align detector
- m_align = tool<IGetElementsToBeAligned>( "GetElementsToBeAligned" );
- if ( !m_align ) return Error( "==> Failed to retrieve detector selector tool!", StatusCode::FAILURE );
-
- sc = m_trackresidualtool.retrieve();
- if ( sc.isFailure() ) return sc;
-
- sc = m_vertexresidualtool.retrieve();
- if ( sc.isFailure() ) return sc;
-
- sc = m_updatetool.retrieve();
- if ( sc.isFailure() ) return sc;
-
- sc = m_vertextrackselector.retrieve();
- if ( sc.isFailure() ) return sc;
+ // m_elementtool = tool<IGetElementsToBeAligned>( "GetElementsToBeAligned" );
+ // if ( !m_elementtool ) return Error( "==> Failed to retrieve detector selector tool!", StatusCode::FAILURE );
/// Get range detector elements
- const Elements& elements = m_align->elements();
+ const Elements& elements = m_elementtool->elements();
if ( printDebug() ) {
debug() << "==> Got " << elements.size() << " elements to align!" << endmsg;
@@ -205,10 +153,6 @@ StatusCode AlignAlgorithm::initialize() {
if ( !sc.isSuccess() ) { error() << "cannot register object. statuscode = " << sc << endmsg; }
m_equations = &( m_summaryData->equations() );
- /// Get projector selector tool
- m_projSelector = tool<ITrackProjectorSelector>( m_projSelectorName, "Projector", this );
- if ( !m_projSelector ) return Error( "==> Failed to retrieve projector selector tool", StatusCode::FAILURE );
-
/// Monitoring
/// Book residual histograms
/// Residuals
@@ -235,9 +179,6 @@ StatusCode AlignAlgorithm::finalize() {
for ( std::vector<AlElementHistos*>::iterator ielem = m_elemHistos.begin(); ielem != m_elemHistos.end(); ++ielem )
delete *ielem;
- m_trackresidualtool.release().ignore();
- m_vertexresidualtool.release().ignore();
- m_updatetool.release().ignore();
if ( m_Online ) { m_HistoUpdater->m_MonSvc = 0; }
return GaudiHistoAlg::finalize();
}
@@ -247,7 +188,7 @@ StatusCode AlignAlgorithm::stop() {
if ( m_RunList.size() > 0 ) {
if ( m_RunList[0] != "*" ) { m_HistoUpdater->Update( ::atoi( m_RunList[0].c_str() ) ); }
}
- if ( !m_outputDataFileName.empty() ) m_equations->writeToFile( m_outputDataFileName.c_str() );
+ if ( !m_outputDataFileName.empty() ) m_equations->writeToFile( m_outputDataFileName.value().c_str() );
m_fireRunChange = true;
}
return StatusCode::SUCCESS;
@@ -281,8 +222,8 @@ StatusCode AlignAlgorithm::execute() {
if ( m_resetHistos ) resetHistos();
// Make sure that the alignment frames are properly initialized on the first event
- if ( m_equations->initTime() != 0 && m_equations->initTime() != m_align->initTime() ) {
- error() << "Mismatch in init time: " << m_equations->initTime() << " " << m_align->initTime() << ". Aborting."
+ if ( m_equations->initTime() != 0 && m_equations->initTime() != m_elementtool->initTime() ) {
+ error() << "Mismatch in init time: " << m_equations->initTime() << " " << m_elementtool->initTime() << ". Aborting."
<< endmsg;
return StatusCode::FAILURE;
}
@@ -296,10 +237,10 @@ StatusCode AlignAlgorithm::execute() {
if ( m_equations->initTime() == 0 ) {
if ( m_forcedInitTime )
- m_align->initAlignmentFrame( m_forcedInitTime ).ignore();
+ m_elementtool->initAlignmentFrame( m_forcedInitTime.value() ).ignore();
else
- m_align->initAlignmentFrame( eventtime ).ignore();
- m_align->initEquations( *m_equations );
+ m_elementtool->initAlignmentFrame( eventtime ).ignore();
+ m_elementtool->initEquations( *m_equations );
}
// Get tracks. Copy them into a vector, because that's what we need for dealing with vertices.
@@ -375,7 +316,7 @@ StatusCode AlignAlgorithm::execute() {
if ( !m_particleLocation.empty() ) {
LHCb::Particle::Range particles = get<LHCb::Particle::Range>( m_particleLocation );
for ( const LHCb::Particle* p : particles ) {
- const Al::MultiTrackResiduals* res = m_vertexresidualtool->get( *p );
+ const auto res = m_vertexresidualtool->get( *p );
if ( res && accumulate( *res ) ) {
m_equations->addVertexChi2Summary( res->vertexChi2(), res->vertexNDoF() );
++numuseddimuons;
@@ -393,7 +334,7 @@ StatusCode AlignAlgorithm::execute() {
// rebuild this vertex from tracks in the selected track list
const LHCb::RecVertex* clone = cloneVertex( **ivertex, selectedtracks );
if ( clone ) {
- const Al::MultiTrackResiduals* res = m_vertexresidualtool->get( *clone );
+ const auto res = m_vertexresidualtool->get( *clone );
if ( res && accumulate( *res ) ) {
m_equations->addVertexChi2Summary( res->vertexChi2(), res->vertexNDoF() );
++numuseddimuons;
@@ -417,7 +358,7 @@ StatusCode AlignAlgorithm::execute() {
VertexContainer splitvertices;
splitVertex( **ivertex, selectedtracks, splitvertices );
for ( VertexContainer::iterator isubver = splitvertices.begin(); isubver != splitvertices.end(); ++isubver ) {
- const Al::MultiTrackResiduals* res = m_vertexresidualtool->get( *isubver );
+ const auto res = m_vertexresidualtool->get( *isubver );
if ( res && accumulate( *res ) ) {
m_equations->addVertexChi2Summary( res->vertexChi2(), res->vertexNDoF() );
++numusedvertices;
@@ -611,14 +552,14 @@ void AlignAlgorithm::update() {
warning() << "Adding derivatives from input file: " << *ifile << " " << tmp.numHits() << " "
<< tmp.totalChiSquare() << " " << m_equations->totalChiSquare() << endmsg;
}
- if ( m_align->initTime() == 0 ) m_align->initAlignmentFrame( m_equations->initTime() ).ignore();
+ if ( m_elementtool->initTime() == 0 ) m_elementtool->initAlignmentFrame( m_equations->initTime() ).ignore();
}
update( *m_equations );
}
}
void AlignAlgorithm::update( const Al::Equations& equations ) {
- if ( !m_outputDataFileName.empty() ) m_equations->writeToFile( m_outputDataFileName.c_str() );
+ if ( !m_outputDataFileName.empty() ) m_equations->writeToFile( m_outputDataFileName.value().c_str() );
m_updatetool->process( equations, m_iteration, m_nIterations ).ignore();
// write the xml
m_xmlWritersTool->write().ignore();
diff --git a/Alignment/TAlignment/src/AlignAlgorithm.h b/Alignment/TAlignment/src/AlignAlgorithm.h
index 9b3c0ffc14400b6b98f48ae051af8a747918e6ab..272b9f98bc97bc59df81db4cb337bf47c3557386 100755
--- a/Alignment/TAlignment/src/AlignAlgorithm.h
+++ b/Alignment/TAlignment/src/AlignAlgorithm.h
@@ -8,13 +8,10 @@
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
-// $Id: AlignAlgorithm.h,v 1.39 2010-04-17 16:21:39 wouter Exp $
-#ifndef TALIGNMENT_ALIGNALGORITHM_H
-#define TALIGNMENT_ALIGNALGORITHM_H 1
+#pragma once
// Include files
// from std
-#include <map>
#include <string>
#include <utility>
#include <vector>
@@ -75,7 +72,6 @@ namespace Al {
* @date 2007-03-05
*/
-class IAlignWork;
class HistoUpdater;
class AlignAlgorithm : public GaudiHistoAlg, virtual public IIncidentListener {
@@ -109,14 +105,6 @@ public:
void reset();
- /** Method to get alignment constants, posXYZ and rotXYZ for a given set
- * of detector elements
- * @param elements flat vector of detector elements, i.e. std::vector<DetectorElements>
- * @param alignConstants reference to a flat vector of alignment constants, i.e. std::vector<double>
- */
- void getAlignmentConstants( const Elements& elements, AlignConstants& alignConstants ) const;
- std::vector<std::string> m_RunList;
-
protected:
bool printDebug() const { return msgLevel( MSG::DEBUG ); };
bool printVerbose() const { return msgLevel( MSG::VERBOSE ); };
@@ -133,50 +121,54 @@ protected:
LHCb::RecVertex* cloneVertex( const LHCb::RecVertex& vertex, const TrackContainer& selectedtracks ) const;
bool testNodes( const LHCb::Track& track ) const;
- HistoUpdater* m_HistoUpdater;
+private:
+ // properties
+ Gaudi::Property<size_t> m_nIterations{this, "NumberOfIterations", 10u, "Number of iterations used on histogram axes"};
+ PublicToolHandle<IGetElementsToBeAligned> m_elementtool{
+ this, "AlignmentElements", "GetElementsToBeAligned"}; ///< Pointer to tool to align detector
+ Gaudi::Property<std::string> m_tracksLocation{this, "TracksLocation", LHCb::TrackLocation::Default};
+ Gaudi::Property<std::string> m_vertexLocation{this, "VertexLocation"}; ///< Vertices for alignment
+ Gaudi::Property<std::string> m_dimuonLocation{this, "DimuonLocation"}; ///< J/psi vertcies for alignment
+ Gaudi::Property<std::string> m_particleLocation{this,
+ "ParticleLocation"}; ///< particles with mass constraint for alignment
+ ToolHandle<ITrackProjectorSelector> m_projSelector{this, "ProjectorSelector",
+ "TrackProjectorSelector"}; ///< Pointer to projector selector tool
+ PublicToolHandle<Al::ITrackResidualTool> m_trackresidualtool{this, "TrackResidualTool", "Al::TrackResidualTool"};
+ PublicToolHandle<Al::IVertexResidualTool> m_vertexresidualtool{this, "VertexResidualTool", "Al::VertexResidualTool"};
+ PublicToolHandle<Al::IAlignUpdateTool> m_updatetool{this, "UpdateTool", "Al::AlignUpdateTool"};
+ ToolHandle<ITrackSelector> m_vertextrackselector{this, "VertexTrackSelector", "TrackSelector"};
+ PublicToolHandle<IWriteAlignmentConditionsTool> m_xmlWritersTool{this, "XmlWritersTool",
+ "WriteMultiAlignmentConditionsTool"};
+ Gaudi::Property<std::string> m_alignSummaryDataSvc{this, "AlignSummaryDataSvc", "DetectorDataSvc"};
+ Gaudi::Property<std::string> m_alignSummaryLocation{this, "AlignSummaryLocation", "AlignDerivativeData"};
+ Gaudi::Property<bool> m_fillHistos{this, "FillHistos", false};
+
+ Gaudi::Property<bool> m_correlation{this, "UseCorrelations",
+ true}; ///< Do we take into account correlations between residuals?
+ Gaudi::Property<bool> m_updateInFinalize{this, "UpdateInFinalize", false}; ///< Call update from finalize
+ Gaudi::Property<double> m_chi2Outlier{this, "Chi2Outlier", 1e4}; ///< Chi2 threshold for outlier rejection
+ Gaudi::Property<size_t> m_minTracksPerVertex{this, "MinTracksPerVertex", 2};
+ Gaudi::Property<size_t> m_maxTracksPerVertex{this, "MaxTracksPerVertex", 8};
+ Gaudi::Property<std::string> m_outputDataFileName{this, "OutputDataFile"};
+ Gaudi::Property<std::vector<std::string>> m_inputDataFileNames{this, "InputDataFiles"};
+ Gaudi::Property<long long> m_forcedInitTime{
+ this, "ForcedInitialTime", 0,
+ "force the alignment geometry to initialize with this time (rather than first event)"};
+ Gaudi::Property<bool> m_Online{this, "OnlineMode", false};
+ Gaudi::Property<std::vector<std::string>> m_RunList{this, "RunList"};
private:
- size_t m_iteration; ///< Iteration counter
- size_t m_nIterations; ///< Number of iterations
- size_t m_nTracks; ///< Number of tracks used
- size_t m_covFailure; ///< Number of covariance calculation failures
- IGetElementsToBeAligned* m_align; ///< Pointer to tool to align detector
- std::string m_tracksLocation; ///< Tracks for alignment
- std::string m_vertexLocation; ///< Vertices for alignment
- std::string m_dimuonLocation; ///< J/psi vertcies for alignment
- std::string m_particleLocation; ///< particles with mass constraint for alignment
- std::string m_projSelectorName; ///< Name of projector selector tool
- ITrackProjectorSelector* m_projSelector; ///< Pointer to projector selector tool
- ToolHandle<Al::ITrackResidualTool> m_trackresidualtool;
- ToolHandle<Al::IVertexResidualTool> m_vertexresidualtool;
- ToolHandle<Al::IAlignUpdateTool> m_updatetool;
- ToolHandle<ITrackSelector> m_vertextrackselector;
- ToolHandle<IWriteAlignmentConditionsTool> m_xmlWritersTool;
- // std::vector<std::string> m_xmlWriterNames ;
- // std::vector<IWriteAlignmentConditionsTool*> m_xmlWriters ;
- std::string m_alignSummaryDataSvc;
- std::string m_alignSummaryLocation;
- bool m_fillHistos;
- bool m_fireRunChange;
-
- Al::Equations* m_equations; ///< Equations to solve
- bool m_correlation; ///< Do we take into account correlations between residuals?
- bool m_updateInFinalize; ///< Call update from finalize
- double m_chi2Outlier; ///< Chi2 threshold for outlier rejection
- bool m_usePreconditioning; ///< Precondition the system of equations before calling solver
- size_t m_minTracksPerVertex;
- size_t m_maxTracksPerVertex;
- std::string m_outputDataFileName;
- std::vector<std::string> m_inputDataFileNames;
- unsigned int m_runnr;
+ // state (some of these can be changed in to counters)
+ size_t m_iteration{0u}; ///< Iteration counter
+ size_t m_nTracks{0u}; ///< Number of tracks used
+ size_t m_covFailure{0u}; ///< Number of covariance calculation failures
+ Al::Equations* m_equations{nullptr}; ///< Equations to solve
+ bool m_fireRunChange{false};
+ unsigned int m_runnr{0};
+ bool m_resetHistos{false}; // reset histos on next event processing
+ HistoUpdater* m_HistoUpdater{nullptr};
/// Monitoring
// @todo: Move this to a monitoring tool
std::vector<AlElementHistos*> m_elemHistos;
- bool m_resetHistos; // reset histos on next event processing
- long long m_forcedInitTime; // force the alignment geometry to initialize with this time (rather than first event)
- bool m_Online;
- IAlignWork* m_IAlwork;
};
-
-#endif // TALIGNMENT_ALIGNALGORITHM_H
diff --git a/Alignment/TAlignment/src/AlignChisqConstraintTool.cpp b/Alignment/TAlignment/src/AlignChisqConstraintTool.cpp
index 60b00943f5859b2dc6dacb74db75ca5e420cd2bf..f65033066b96c57687cc4e3a0ba067ebd37622a9 100644
--- a/Alignment/TAlignment/src/AlignChisqConstraintTool.cpp
+++ b/Alignment/TAlignment/src/AlignChisqConstraintTool.cpp
@@ -40,10 +40,9 @@ namespace Al {
double pivot[3];
};
- class AlignChisqConstraintTool : public GaudiTool, virtual public IAlignChisqConstraintTool {
+ class AlignChisqConstraintTool : public extends<GaudiTool, IAlignChisqConstraintTool> {
public:
- AlignChisqConstraintTool( const std::string& type, const std::string& name, const IInterface* parent );
- ~AlignChisqConstraintTool();
+ using extends::extends;
StatusCode initialize() override;
LHCb::ChiSquare addConstraints( const Elements& elements, AlVec& halfDChi2DAlpha, AlSymMat& halfD2Chi2DAlpha2,
std::ostream& stream ) const override;
@@ -59,9 +58,9 @@ namespace Al {
const Gaudi::Vector6* findXmlUncertainty( const std::string& name ) const;
private:
- std::vector<std::string> m_constraintNames;
- std::vector<std::string> m_xmlFiles;
- std::vector<std::string> m_xmlUncertainties;
+ Gaudi::Property<std::vector<std::string>> m_constraintNames{this, "Constraints"};
+ Gaudi::Property<std::vector<std::string>> m_xmlFiles{this, "XmlFiles"};
+ Gaudi::Property<std::vector<std::string>> m_xmlUncertainties{this, "XmlUncertainties"};
typedef std::map<std::string, XmlSurveyData> XmlData;
XmlData m_xmldata;
typedef std::vector<ConfiguredSurveyData> ConstraintData;
@@ -74,18 +73,6 @@ namespace Al {
DECLARE_COMPONENT( AlignChisqConstraintTool )
- AlignChisqConstraintTool::~AlignChisqConstraintTool() {}
-
- AlignChisqConstraintTool::AlignChisqConstraintTool( const std::string& type, const std::string& name,
- const IInterface* parent )
- : GaudiTool( type, name, parent ) {
- // interfaces
- declareInterface<IAlignChisqConstraintTool>( this );
- declareProperty( "Constraints", m_constraintNames );
- declareProperty( "XmlFiles", m_xmlFiles );
- declareProperty( "XmlUncertainties", m_xmlUncertainties );
- }
-
StatusCode AlignChisqConstraintTool::initialize() {
StatusCode sc = GaudiTool::initialize();
for ( std::vector<std::string>::const_iterator ifile = m_constraintNames.begin();
diff --git a/Alignment/TAlignment/src/AlignConstraintTool.cpp b/Alignment/TAlignment/src/AlignConstraintTool.cpp
index 841d76ac6166886ebb70724e50a79335a3b53dc1..2a522ce3cb64f3cb77f4aaf6da5450ae65192986 100755
--- a/Alignment/TAlignment/src/AlignConstraintTool.cpp
+++ b/Alignment/TAlignment/src/AlignConstraintTool.cpp
@@ -91,8 +91,8 @@ namespace Al {
ConstraintDerivatives::ConstraintDerivatives( size_t dim, const std::vector<std::string>& activeconstraints,
const std::string& nameprefix, int offset )
: m_dofmask( size_t( NumConstraints ) )
- , m_derivatives( NumConstraints, dim )
- , m_residuals( NumConstraints )
+ , m_derivatives{size_t( NumConstraints ), dim}
+ , m_residuals{size_t( NumConstraints )}
, m_nameprefix( nameprefix )
, m_activeParOffset( offset ) {
std::vector<bool> active( NumConstraints, false );
@@ -107,11 +107,12 @@ namespace Al {
boost::assign::list_of( "Tx" )( "Ty" )( "Tz" )( "Rx" )( "Ry" )( "Rz" )( "Szx" )( "Szy" )( "Szz" )( "Sxx" )(
"SRz" )( "Trx" )( "Try" )( "Trtx" )( "Trty" )( "Trcur" )( "PVx" )( "PVy" )( "PVz" )( "Sz2x" )( "Sz2y" );
- class AlignConstraintTool : public GaudiTool, virtual public IAlignConstraintTool {
+ class AlignConstraintTool : public extends<GaudiTool, IAlignConstraintTool> {
public:
typedef IGetElementsToBeAligned::Elements Elements;
typedef std::vector<Al::ConstraintDerivatives*> DerivativeContainer;
- AlignConstraintTool( const std::string& type, const std::string& name, const IInterface* parent );
+
+ using extends::extends;
~AlignConstraintTool();
StatusCode initialize() override;
@@ -153,26 +154,16 @@ namespace Al {
};
private:
- std::vector<std::string> m_constraintNames;
- std::vector<std::string> m_chiSquareConstraintNames;
- std::vector<ConstraintDefinition> m_definitions;
- mutable DerivativeContainer m_derivatives;
- bool m_useWeightedAverage;
- bool m_fixRedundantMotherDofs;
- const IGetElementsToBeAligned* m_elementtool;
+ Gaudi::Property<std::vector<std::string>> m_constraintNames{this, "Constraints"};
+ Gaudi::Property<bool> m_useWeightedAverage{this, "UseWeightedAverage", false};
+ Gaudi::Property<bool> m_fixRedundantMotherDofs{this, "FixRedundantMotherDofs", false};
+ PublicToolHandle<IGetElementsToBeAligned> m_elementtool{this, "ElementTool", "GetElementsToBeAligned"};
+ std::vector<ConstraintDefinition> m_definitions;
+ mutable DerivativeContainer m_derivatives;
};
DECLARE_COMPONENT( AlignConstraintTool )
- AlignConstraintTool::AlignConstraintTool( const std::string& type, const std::string& name, const IInterface* parent )
- : GaudiTool( type, name, parent ), m_derivatives( 0 ) {
- // interfaces
- declareInterface<IAlignConstraintTool>( this );
- declareProperty( "Constraints", m_constraintNames );
- declareProperty( "UseWeightedAverage", m_useWeightedAverage = false );
- declareProperty( "FixRedundantMotherDofs", m_fixRedundantMotherDofs = false );
- }
-
AlignConstraintTool::~AlignConstraintTool() { clearConstraintDerivatives(); }
#include <string>
@@ -203,8 +194,6 @@ namespace Al {
// check that constraints exist
StatusCode sc = GaudiTool::initialize();
- m_elementtool = tool<IGetElementsToBeAligned>( "GetElementsToBeAligned" );
-
if ( sc.isSuccess() ) {
// now we need to decode the string with constraints
ConstraintDefinition common( "" );
@@ -480,6 +469,8 @@ namespace Al {
AlVec halfDChi2DAlphaNew( dim );
AlSymMat halfD2Chi2DAlpha2New( dim );
+ debug() << "In AlignConstraintTool: before vec=" << halfDChi2DAlpha << endmsg;
+
// copy the old matrices in there
for ( size_t irow = 0; irow < numpars; ++irow ) {
halfDChi2DAlphaNew( irow ) = halfDChi2DAlpha( irow );
@@ -490,6 +481,7 @@ namespace Al {
// now add the constraints, only to the 2nd derivative (residuals are zero for now)
for ( DerivativeContainer::const_iterator ic = m_derivatives.begin(); ic != m_derivatives.end(); ++ic )
if ( ( *ic )->nActive() > 0 ) {
+ debug() << "residuals: " << ( *ic )->residuals() << endmsg;
for ( size_t irow = 0; irow < Al::ConstraintDerivatives::NumConstraints; ++irow ) {
int iactive = ( *ic )->activeParIndex( irow );
if ( 0 <= iactive ) {
@@ -501,14 +493,17 @@ namespace Al {
for ( size_t irow = 0; irow < Al::ConstraintDerivatives::NumConstraints; ++irow ) {
int iactive = ( *ic )->activeParIndex( irow );
- if ( 0 <= iactive )
+ if ( 0 <= iactive ) {
halfDChi2DAlphaNew( ( *ic )->activeParOffset() + iactive ) = ( *ic )->residuals()( irow );
+ debug() << "constraint residual: " << ( *ic )->residuals()( irow ) << endmsg;
+ }
}
}
// copy the result back
halfDChi2DAlpha = halfDChi2DAlphaNew;
halfD2Chi2DAlpha2 = halfD2Chi2DAlpha2New;
+ debug() << "In AlignConstraintTool: after vec=" << halfDChi2DAlpha << endmsg;
}
return numactive;
@@ -527,11 +522,12 @@ namespace Al {
totalmask.push_back( ( *ic )->dofMask() );
// first extract the part concerning the lagrange multipliers
- AlVec lambda( numConstraintDerivatives );
- AlSymMat covlambda( numConstraintDerivatives );
+ AlVec lambda{numConstraintDerivatives};
+ AlMat covlambda{numConstraintDerivatives, numConstraintDerivatives};
for ( size_t i = 0u; i < numConstraintDerivatives; ++i ) {
- lambda[i] = parameters[numPars + i];
- for ( size_t j = 0u; j <= i; ++j ) covlambda[i][j] = covariance[numPars + i][numPars + j];
+ lambda( i ) = parameters[numPars + i];
+ for ( size_t j = 0u; j < numConstraintDerivatives; ++j )
+ covlambda( i, j ) = covariance( numPars + i, numPars + j );
}
// For a lagrage constraint defined as delta-chisqyare = 2 * lambda * constraint, the solution for lambda is
@@ -543,11 +539,10 @@ namespace Al {
// chisquare = - lambda * constraint
// Note that this all needs to work for vectors.
- AlSymMat V = -1 * covlambda; // Al$%$%^&$%&&*&^
- int ierr;
- V.invert( ierr );
- AlVec x = ( V * lambda ); // This might need a minus sign ...
- double chisquare = ( lambda * V * lambda );
+ const auto minusV = covlambda.inverse();
+ const auto V = -minusV;
+ const auto x = V * lambda; // This might need a minus sign ...
+ const double chisquare = ( lambda.transpose() * V * lambda );
logmessage << "Canonical constraint values: " << std::endl;
for ( size_t iactive = 0u; iactive < numConstraintDerivatives; ++iactive ) {
size_t parindex = totalmask.parIndex( iactive );
@@ -555,23 +550,23 @@ namespace Al {
<< m_derivatives[parindex / Al::ConstraintDerivatives::NumConstraints]->name(
parindex % Al::ConstraintDerivatives::NumConstraints )
<< std::setw( 12 ) << x( iactive ) << " "
- << " +/- " << std::setw( 12 ) << AlParameters::signedSqrt( V.fast( iactive, iactive ) )
- << " gcc^2: " << std::setw( 12 )
- << 1 + 1 / ( V.fast( iactive, iactive ) * covlambda.fast( iactive, iactive ) ) << std::endl;
+ << " +/- " << std::setw( 12 ) << AlParameters::signedSqrt( V( iactive, iactive ) )
+ << " gcc^2: " << std::setw( 12 ) << 1 + 1 / ( V( iactive, iactive ) * covlambda( iactive, iactive ) )
+ << std::endl;
}
logmessage << "Canonical constraint chisquare: " << chisquare << std::endl;
// Now the inactive constraints. Shrink the parameter and covariance matrix to the parameters only part
- AlVec subparameters( numPars );
- AlSymMat subcovariance( numPars );
+ AlVec subparameters{numPars};
+ AlMat subcovariance{numPars, numPars};
for ( size_t irow = 0; irow < numPars; ++irow ) {
subparameters[irow] = parameters[irow];
- for ( size_t icol = 0; icol <= irow; ++icol ) subcovariance.fast( irow, icol ) = covariance.fast( irow, icol );
+ for ( size_t icol = 0; icol < numPars; ++icol ) subcovariance( irow, icol ) = covariance( irow, icol );
}
// Calculate the delta
- AlVec constraintdelta = m_derivatives.front()->derivatives() * subparameters;
- AlMat constraintcov =
- m_derivatives.front()->derivatives() * subcovariance * m_derivatives.front()->derivatives().T();
+ const auto derivatives = m_derivatives.front()->derivatives();
+ const auto constraintdelta = derivatives * subparameters;
+ const auto constraintcov = derivatives * subcovariance * derivatives.transpose();
logmessage << "Values of constraint equations after solution (X=active,O=inactive): " << std::endl;
for ( size_t i = 0; i < Al::ConstraintDerivatives::NumConstraints; ++i ) {
diff --git a/Alignment/TAlignment/src/AlignUpdateTool.cpp b/Alignment/TAlignment/src/AlignUpdateTool.cpp
index 911a55e7f1222fbaac2486b38ca660ae404a21ea..c30fb7a3572cfe1aa08802ab61a31c90a6a24dde 100755
--- a/Alignment/TAlignment/src/AlignUpdateTool.cpp
+++ b/Alignment/TAlignment/src/AlignUpdateTool.cpp
@@ -31,14 +31,16 @@
namespace Al {
- class AlignUpdateTool : public GaudiHistoTool, virtual public IAlignUpdateTool {
+ class AlignUpdateTool : public extends<GaudiHistoTool, IAlignUpdateTool> {
public:
- typedef IGetElementsToBeAligned::Elements Elements;
- typedef std::vector<double> AlignConstants;
- AlignUpdateTool( const std::string& type, const std::string& name, const IInterface* parent );
- ~AlignUpdateTool();
+ /// Typedefs etc.
+ using Elements = IGetElementsToBeAligned::Elements;
+ using AlignConstants = std::vector<double>;
+
+ /// Standard constructor
+ using extends::extends;
+
StatusCode initialize() override;
- StatusCode finalize() override;
StatusCode process( const Al::Equations& equations, size_t iter, size_t maxiter ) const override;
StatusCode process( const Al::Equations& equations, ConvergenceStatus& convergencestatus ) const override;
StatusCode process( const Al::Equations& equations, ConvergenceStatus& convergencestatus, size_t iter,
@@ -64,81 +66,45 @@ namespace Al {
static double computeLookElsewhereChi2Cut( double chi2, double n );
private:
- std::string m_matrixSolverToolName; ///< Name of linear algebra solver tool
- IAlignSolvTool* m_matrixSolverTool; ///< Pointer to linear algebra solver tool
- IGetElementsToBeAligned* m_elementProvider;
- ToolHandle<IAlignConstraintTool> m_lagrangeconstrainttool;
- ToolHandle<IAlignChisqConstraintTool> m_chisqconstrainttool;
- size_t m_minNumberOfHits; ///< Minimum number of hits for an Alignable to be aligned
- bool m_usePreconditioning; ///< Precondition the system of equations before calling solver
- double m_maxDeltaChi2PDofForConvergence;
- double m_maxModeDeltaChi2ForConvergence;
- std::string m_logFileName;
- mutable std::ostringstream m_logMessage;
- mutable unsigned int m_iteration;
+ PublicToolHandle<IAlignSolvTool> m_matrixSolverTool{this, "MatrixSolverTool", "EigenDiagSolvTool",
+ "Linear algebra solver tool"};
+ PublicToolHandle<IGetElementsToBeAligned> m_elementProvider{this, "ElementProvider", "GetElementsToBeAligned"};
+ ToolHandle<IAlignConstraintTool> m_lagrangeconstrainttool{this, "LagrangeConstraintTool",
+ "Al::AlignConstraintTool"};
+ ToolHandle<IAlignChisqConstraintTool> m_chisqconstrainttool{this, "SurveyConstraintTool",
+ "Al::AlignChisqConstraintTool"};
+ Gaudi::Property<size_t> m_minNumberOfHits{this, "MinNumberOfHits", 100u,
+ "Minimum number of hits for an Alignable to be aligned"};
+ Gaudi::Property<bool> m_usePreconditioning{this, "UsePreconditioning", true,
+ "Precondition the system of equations before calling solver"};
+ Gaudi::Property<double> m_maxDeltaChi2PDofForConvergence{this, "MaxDeltaChi2PDofForConvergence", 4};
+ Gaudi::Property<double> m_maxModeDeltaChi2ForConvergence{this, "MaxModeDeltaChi2ForConvergence", 25};
+ Gaudi::Property<std::string> m_logFileName{this, "LogFile", "alignlog.txt"};
+ Gaudi::Property<bool> m_skipGeometryUpdate{this, "SkipGeometryUpdate", false, "Skip the update of the geometry"};
+ mutable std::ostringstream m_logMessage;
+ mutable unsigned int m_iteration{0};
mutable std::map<std::string, double> m_initConsts;
mutable std::map<std::string, double> m_finalConsts;
};
DECLARE_COMPONENT( AlignUpdateTool )
- AlignUpdateTool::AlignUpdateTool( const std::string& type, const std::string& name, const IInterface* parent )
- : GaudiHistoTool( type, name, parent )
- , m_lagrangeconstrainttool( "Al::AlignConstraintTool", this )
- , m_chisqconstrainttool( "Al::AlignChisqConstraintTool", this )
- , m_iteration( 0 ) {
- // interfaces
- declareInterface<IAlignUpdateTool>( this );
- declareProperty( "MatrixSolverTool", m_matrixSolverToolName = "SpmInvTool" );
- declareProperty( "MinNumberOfHits", m_minNumberOfHits = 100u );
- declareProperty( "UsePreconditioning", m_usePreconditioning = true );
- declareProperty( "LogFile", m_logFileName = "alignlog.txt" );
- declareProperty( "SurveyConstraintTool", m_chisqconstrainttool );
- declareProperty( "LagrangeConstraintTool", m_lagrangeconstrainttool );
- declareProperty( "MaxDeltaChi2PDofForConvergence", m_maxDeltaChi2PDofForConvergence = 4 );
- declareProperty( "MaxModeDeltaChi2ForConvergence", m_maxModeDeltaChi2ForConvergence = 25 );
- }
-
- AlignUpdateTool::~AlignUpdateTool() {}
-
StatusCode AlignUpdateTool::initialize() {
// check that constraints exist
StatusCode sc = GaudiHistoTool::initialize();
if ( sc.isFailure() ) return sc; // error printed already by GaudiAlgorithm
- setHistoTopDir( "" );
- setHistoDir( "Alignment" );
-
- // Get tool to align detector, from the tool service
- m_elementProvider = tool<IGetElementsToBeAligned>( "GetElementsToBeAligned" );
- if ( !m_elementProvider ) return Error( "==> Failed to retrieve detector selector tool!", StatusCode::FAILURE );
-
- // Get the constraint tool, also from the toolsvc
- sc = m_lagrangeconstrainttool.retrieve();
- if ( !sc.isSuccess() ) return Error( "==> Failed to retrieve lagrange constraint tool!", StatusCode::FAILURE );
-
- // Get the constraint tool, also from the toolsvc
- sc = m_chisqconstrainttool.retrieve();
- if ( !sc.isSuccess() ) return Error( "==> Failed to retrieve chisq constraint tool!", StatusCode::FAILURE );
-
- // Get matrix solver tool
- // m_matrixSolverTool = tool<IAlignSolvTool>(m_matrixSolverToolName, "MatrixSolver", this);
- m_matrixSolverTool = tool<IAlignSolvTool>( m_matrixSolverToolName );
- if ( !m_matrixSolverTool ) return Error( "==> Failed to retrieve matrix solver tool", StatusCode::FAILURE );
-
+ if ( name() == std::string( "ToolSvc.AlignUpdateTool" ) ) {
+ setHistoTopDir( "" );
+ setHistoDir( "Alignment" );
+ }
return sc;
}
- StatusCode AlignUpdateTool::finalize() {
- m_lagrangeconstrainttool.release().ignore();
- m_chisqconstrainttool.release().ignore();
- return GaudiHistoTool::finalize();
- }
-
void AlignUpdateTool::preCondition( const Elements& elements, const Al::Equations& equations, AlVec& halfDChi2DAlpha,
AlSymMat& halfD2Chi2DAlpha2, AlVec& scale ) const {
// This is just not sufficiently fool proof!
- size_t size = halfDChi2DAlpha.size();
- scale.reSize( size );
+ size_t size = halfDChi2DAlpha.rows();
+ scale = AlVec{size};
size_t ipar( 0 );
// initialize all scales with 1
for ( ipar = 0u; ipar < size; ++ipar ) scale[ipar] = 1;
@@ -151,7 +117,7 @@ namespace Al {
size_t N = std::max( equations.element( iElem ).numHits(), size_t( 1 ) );
for ( ipar = offset; ipar < offset + ndof; ++ipar ) {
assert( ipar < size );
- if ( halfD2Chi2DAlpha2[ipar][ipar] > 0 ) scale[ipar] = std::sqrt( N / halfD2Chi2DAlpha2[ipar][ipar] );
+ if ( halfD2Chi2DAlpha2( ipar, ipar ) > 0 ) scale[ipar] = std::sqrt( N / halfD2Chi2DAlpha2( ipar, ipar ) );
}
}
}
@@ -177,25 +143,25 @@ namespace Al {
// into account the scales of the parameters)
double norm2( 0 );
for ( size_t jpar = 0; jpar < numParameters; ++jpar ) {
- double derivative = halfD2Chi2DAlpha2[ipar][jpar] * scale[jpar];
+ double derivative = halfD2Chi2DAlpha2( ipar, jpar ) * scale( jpar );
norm2 += derivative * derivative;
}
// now set the scale
- if ( norm2 > 0 ) scale[ipar] = defaultnorm / std::sqrt( norm2 ) * ( ipar - numParameters + 1 );
+ if ( norm2 > 0 ) scale( ipar ) = defaultnorm / std::sqrt( norm2 ) * ( ipar - numParameters + 1 );
}
// now apply the scales
for ( size_t i = 0u, iEnd = size; i < iEnd; ++i ) {
- halfDChi2DAlpha[i] *= scale[i];
- for ( size_t j = 0u; j <= i; ++j ) halfD2Chi2DAlpha2[i][j] *= scale[i] * scale[j];
+ halfDChi2DAlpha( i ) *= scale( i );
+ for ( size_t j = 0u; j <= i; ++j ) halfD2Chi2DAlpha2( i, j ) *= scale( i ) * scale( j );
}
}
void AlignUpdateTool::postCondition( AlVec& halfDChi2DAlpha, AlSymMat& halfD2Chi2DAlpha2, const AlVec& scale ) const {
size_t size = halfDChi2DAlpha.size();
for ( size_t i = 0u, iEnd = size; i < iEnd; ++i ) {
- halfDChi2DAlpha[i] *= scale[i];
- for ( size_t j = 0u; j <= i; ++j ) halfD2Chi2DAlpha2[i][j] *= scale[i] * scale[j];
+ halfDChi2DAlpha( i ) *= scale( i );
+ for ( size_t j = 0u; j <= i; ++j ) halfD2Chi2DAlpha2( i, j ) *= scale( i ) * scale( j );
}
}
@@ -359,9 +325,15 @@ namespace Al {
bool AlignUpdateTool::checkGeometry( const Al::Equations& equations ) const {
// test that the alignment constants in AlEquations actually correspond to what is in the geometry
- const double maxmag2 = 1e-10; // could be finetuned a bit
- bool checkpassed = true;
- const Elements& elements = m_elementProvider->elements();
+ const Elements& elements = m_elementProvider->elements();
+ if ( elements.size() != equations.nElem() ) {
+ error() << "Number of elements in elementprovider and equations do not match:" << elements.size() << " "
+ << equations.nElem() << endmsg;
+ return false;
+ }
+
+ const double maxmag2 = 1e-10; // could be finetuned a bit
+ bool checkpassed = true;
for ( auto ielem : elements ) {
bool thisBad = false;
// compare the current delta
@@ -403,7 +375,8 @@ namespace Al {
if ( m_elementProvider->initTime() != equations.initTime() || !checkGeometry( constequations ) ) {
warning() << "Geometry does not match AlEquations: " << m_elementProvider->initTime().ns() << " "
<< equations.initTime().ns() << ". Will reinitalize alignment frame." << endmsg;
- sc = m_elementProvider->initAlignmentFrame( equations.initTime() );
+ auto nonconstelementprovider = const_cast<IGetElementsToBeAligned*>( &( *m_elementProvider ) );
+ sc = nonconstelementprovider->initAlignmentFrame( equations.initTime() );
if ( !sc.isSuccess() ) return sc;
}
@@ -418,7 +391,6 @@ namespace Al {
return sc;
}
- info() << "\n";
debug() << "==> iteration " << iteration << " : Initial alignment conditions : [";
const Elements& elements = m_elementProvider->elements();
// Store constants before the update
@@ -552,6 +524,17 @@ namespace Al {
info() << "AlSymMat Matrix = " << halfD2Chi2dX2 << endmsg;
} else {
info() << "Matrices too big to dump to std" << endmsg;
+ // check that second derivative looks okay:
+ // 1. check diagonal elements:
+ for ( size_t i = 0; i < numParameters; ++i )
+ if ( halfD2Chi2dX2( i, i ) <= 0 )
+ info() << "Negative diagonal element: " << i << " " << halfD2Chi2dX2( i, i ) << endmsg;
+ for ( size_t i = 0; i < numParameters; ++i )
+ for ( size_t j = 0; j < i; ++j ) {
+ const auto eta = halfD2Chi2dX2( i, j ) / std::sqrt( halfD2Chi2dX2( i, i ) * halfD2Chi2dX2( j, j ) );
+ if ( eta < -1 || eta > 1 )
+ info() << "Off-diagonal element out of range: i, j, eta:" << i << "," << j << "," << eta << endmsg;
+ }
}
// Tool returns H^-1 and alignment constants. Need to copy the 2nd derivative because we still need it.
@@ -564,17 +547,26 @@ namespace Al {
}
info() << "Calling solver tool" << endmsg;
IAlignSolvTool::SolutionInfo solinfo;
- sc = m_matrixSolverTool->compute( covmatrix, solution, solinfo );
+ using clock = std::chrono::high_resolution_clock;
+ const auto t1 = clock::now();
+ sc = m_matrixSolverTool->compute( covmatrix, solution, solinfo );
+ const auto t2 = clock::now();
+ logmessage << "Time to compute solution [ms]: "
+ << std::chrono::duration<double, std::nano>( t2 - t1 ).count() * 1.0e-6 << std::endl;
if ( sc.isSuccess() ) {
-
// undo the scaling
if ( m_usePreconditioning ) {
info() << "Applying post-conditioning" << endmsg;
postCondition( solution, covmatrix, scale );
}
+ // test the numerical accuracy of the solution
+ const AlVec errorvec = halfD2Chi2dX2 * solution - halfDChi2dX;
+ const double error2 = errorvec.dot( errorvec );
+ logmessage << "Error2: " << error2 << std::endl;
+
info() << "computing delta-chi2" << endmsg;
- double deltaChi2 = solution * halfDChi2dX; // somewhere we have been a bit sloppy with two minus signs!
+ double deltaChi2 = solution.dot( halfDChi2dX ); // somewhere we have been a bit sloppy with two minus signs!
if ( deltaChi2 < 0 ) {
error() << "Serious problem in update: delta-chi2 is negative: " << deltaChi2 << endmsg;
//<< "Will recopute solution without scaling to see what that does! " << endmsg ;
@@ -656,8 +648,9 @@ namespace Al {
<< " curtot= " << std::setw( 12 ) << reftotaldelta.parameters()[iactive]
<< " cur= " << std::setw( 12 ) << refdelta.parameters()[iactive]
<< " delta= " << std::setw( 12 ) << delta.parameters()[iactive] << " +/- " << std::setw( 12 )
- << AlParameters::signedSqrt( delta.covariance()[iactive][iactive] ) << " [ " << std::setw( 12 )
- << AlParameters::signedSqrt( 1 / delta.weightMatrix()[iactive][iactive] ) << " ] "
+ << AlParameters::signedSqrt( delta.covariance()( iactive, iactive ) ) << " [ "
+ << std::setw( 12 ) << AlParameters::signedSqrt( 1 / delta.weightMatrix()( iactive, iactive ) )
+ << " ] "
<< " gcc= " << delta.globalCorrelationCoefficient( iactive ) << std::endl;
double contributionToCoordinateError = delta.measurementCoordinateSigma( elemdata.weightR() );
double coordinateError = std::sqrt( elemdata.numHits() / elemdata.weightV() );
@@ -676,8 +669,10 @@ namespace Al {
totalLocalDeltaChi2 += thisLocalDeltaChi2;
// need const_cast because loki range givess access only to const values
- StatusCode sc = ( const_cast<AlignmentElement*>( *it ) )->updateGeometry( delta );
- if ( !sc.isSuccess() ) error() << "Failed to set alignment condition for " << ( *it )->name() << endmsg;
+ if ( !m_skipGeometryUpdate ) {
+ StatusCode sc = ( const_cast<AlignmentElement*>( *it ) )->updateGeometry( delta );
+ if ( !sc.isSuccess() ) error() << "Failed to set alignment condition for " << ( *it )->name() << endmsg;
+ }
std::string name = ( *it )->name();
std::string dirname = "element" + boost::lexical_cast<std::string>( ( *it )->index() ) + "/"; // name + "/"
@@ -768,7 +763,7 @@ namespace Al {
// FIXME: we don't want to do this on every update! Yet, we cannot wait for finalize in online.
if ( !m_logFileName.empty() ) {
- std::ofstream logfile( m_logFileName.c_str() );
+ std::ofstream logfile( m_logFileName.value().c_str() );
logfile << m_logMessage.str() << std::endl;
}
return sc;
diff --git a/Alignment/TAlignment/src/CountingPrescaler.cpp b/Alignment/TAlignment/src/CountingPrescaler.cpp
index 4dc4b00142c936b130e4a1a924944e127bd24384..9803d7f13a0c0e382bf5ae8ecf2284d74ae1eb9a 100644
--- a/Alignment/TAlignment/src/CountingPrescaler.cpp
+++ b/Alignment/TAlignment/src/CountingPrescaler.cpp
@@ -16,27 +16,19 @@
class CountingPrescaler : public GaudiAlgorithm {
public:
- // Constructors and destructor
- CountingPrescaler( const std::string& name, ISvcLocator* pSvcLocator );
+ using GaudiAlgorithm::GaudiAlgorithm;
StatusCode initialize() override;
StatusCode execute() override;
private:
- size_t m_interval;
- size_t m_offset;
- size_t m_printfreq;
- size_t m_counter;
+ Gaudi::Property<size_t> m_interval{this, "Interval", 1};
+ Gaudi::Property<size_t> m_offset{this, "Offset", 0};
+ Gaudi::Property<size_t> m_printfreq{this, "PrintFreq", 0};
+ size_t m_counter{0};
};
DECLARE_COMPONENT( CountingPrescaler )
-CountingPrescaler::CountingPrescaler( const std::string& name, ISvcLocator* pSvcLocator )
- : GaudiAlgorithm( name, pSvcLocator ), m_counter( 0 ) {
- declareProperty( "Interval", m_interval = 1 );
- declareProperty( "Offset", m_offset = 0 );
- declareProperty( "PrintFreq", m_printfreq = 0 );
-}
-
StatusCode CountingPrescaler::initialize() {
info() << "Interval, offset = " << m_interval << ", " << m_offset << endmsg;
m_counter = 0;
diff --git a/Alignment/TAlignment/src/GetElementsToBeAligned.cpp b/Alignment/TAlignment/src/GetElementsToBeAligned.cpp
index e96584e94b260e64fe3c5df5953c2d061d4e959e..23091e87c6eca006eaaf122f163cb67ffa740c8e 100755
--- a/Alignment/TAlignment/src/GetElementsToBeAligned.cpp
+++ b/Alignment/TAlignment/src/GetElementsToBeAligned.cpp
@@ -46,16 +46,6 @@
// Declaration of the Tool Factory
DECLARE_COMPONENT( GetElementsToBeAligned )
-GetElementsToBeAligned::GetElementsToBeAligned( const std::string& type, const std::string& name,
- const IInterface* parent )
- : GaudiTool( type, name, parent ), m_useLocalFrame( true ), m_elemsToBeAligned(), m_elementMap(), m_initTime( 0 ) {
- declareInterface<IGetElementsToBeAligned>( this );
- declareProperty( "Elements", m_elemsToBeAligned );
- declareProperty( "UseLocalFrame", m_useLocalFrame );
-}
-
-GetElementsToBeAligned::~GetElementsToBeAligned() {}
-
StatusCode GetElementsToBeAligned::initialize() {
StatusCode sc = GaudiTool::initialize();
diff --git a/Alignment/TAlignment/src/GetElementsToBeAligned.h b/Alignment/TAlignment/src/GetElementsToBeAligned.h
index d868fb3b21d9d13a0c6958498b4db2bb102072dc..c1dadd1f8cfa026b2dbb1c761e688bbc3fb9d84c 100755
--- a/Alignment/TAlignment/src/GetElementsToBeAligned.h
+++ b/Alignment/TAlignment/src/GetElementsToBeAligned.h
@@ -8,9 +8,7 @@
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
-// $Id: GetElementsToBeAligned.h,v 1.17 2010-02-19 08:54:36 rlambert Exp $
-#ifndef GETELEMENTSTOBEALIGNED_H
-#define GETELEMENTSTOBEALIGNED_H 1
+#pragma once
// Include files
// from STL
@@ -37,8 +35,7 @@
* @date 2007-10-08
*/
-class GetElementsToBeAligned : public GaudiTool, virtual public IGetElementsToBeAligned {
-
+class GetElementsToBeAligned : public extends<GaudiTool, IGetElementsToBeAligned> {
public:
typedef IDetectorElement::IDEContainer::const_iterator IDetIter; ///< IDetector element iterator
typedef boost::char_separator<char> Separator; ///< Char separator
@@ -46,9 +43,7 @@ public:
typedef std::list<boost::regex> RegExs; ///< List of regular expressions
/// Standard constructor
- GetElementsToBeAligned( const std::string& type, const std::string& name, const IInterface* parent );
-
- ~GetElementsToBeAligned(); ///< Destructor
+ using extends::extends;
StatusCode initialize() override;
StatusCode finalize() override;
@@ -134,12 +129,10 @@ private:
};
private:
- bool m_useLocalFrame; ///< Use local frame as alignmentframe
- std::vector<std::string> m_elemsToBeAligned; ///< Elemenst : Path to elements
- mutable IGetElementsToBeAligned::Elements m_elements; ///< Flat vector of alignment elements
- typedef std::map<const DetectorElement*, const AlignmentElement*> ElementMap;
+ Gaudi::Property<bool> m_useLocalFrame{this, "UseLocalFrame", true}; ///< Use local frame as alignmentframe
+ Gaudi::Property<std::vector<std::string>> m_elemsToBeAligned{this, "Elements"}; ///< Elemenst : Path to elements
+ mutable IGetElementsToBeAligned::Elements m_elements; ///< Flat vector of alignment elements
+ using ElementMap = std::map<const DetectorElement*, const AlignmentElement*>;
ElementMap m_elementMap; ///< Map of detector elements to alignment element
- Gaudi::Time m_initTime;
+ Gaudi::Time m_initTime{0};
};
-
-#endif // GETELEMENTSTOBEALIGNED_H
diff --git a/Alignment/TAlignment/src/IAlignChisqConstraintTool.h b/Alignment/TAlignment/src/IAlignChisqConstraintTool.h
index 6433ef2dc11eaaba7592d9ad895b5618a1c8cc52..ace5c906de6e6562a649447106f04492f7935ee7 100644
--- a/Alignment/TAlignment/src/IAlignChisqConstraintTool.h
+++ b/Alignment/TAlignment/src/IAlignChisqConstraintTool.h
@@ -8,36 +8,27 @@
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
-#ifndef TALIGNMENT_IALIGNCHISQCONSTRAINTTOOL_H
-# define TALIGNMENT_IALIGNCHISQCONSTRAINTTOOL_H
+#pragma once
-# include "AlignKernel/AlEquations.h"
-# include "AlignKernel/AlSymMat.h"
-# include "AlignKernel/AlVec.h"
-# include "Event/ChiSquare.h"
-# include "GaudiKernel/IAlgTool.h"
-# include "IGetElementsToBeAligned.h"
+#include "AlignKernel/AlEquations.h"
+#include "AlignKernel/AlSymMat.h"
+#include "AlignKernel/AlVec.h"
+#include "Event/ChiSquare.h"
+#include "GaudiKernel/IAlgTool.h"
+#include "IGetElementsToBeAligned.h"
namespace Al {
- /// Static ID object
- static const InterfaceID IID_IAlignChisqConstraintTool( "Al::IAlignChisqConstraintTool", 0, 0 );
+ struct IAlignChisqConstraintTool : public extend_interfaces<IAlgTool> {
+ DeclareInterfaceID( IAlignChisqConstraintTool, 2, 0 );
- class IAlignChisqConstraintTool : virtual public IAlgTool {
- public:
typedef IGetElementsToBeAligned::Elements Elements;
- // Retrieve interface ID
- static const InterfaceID& interfaceID() { return IID_IAlignChisqConstraintTool; }
- virtual LHCb::ChiSquare addConstraints( const Elements& elements, AlVec& halfDChi2DAlpha,
- AlSymMat& halfD2Chi2DAlpha2, std::ostream& logmessage ) const = 0;
- virtual LHCb::ChiSquare addConstraints( const Elements& inputelements, Al::Equations& equations,
- std::ostream& logmessage ) const = 0;
- virtual LHCb::ChiSquare chiSquare( const AlignmentElement& element, bool activeonly = true ) const = 0;
+ virtual LHCb::ChiSquare addConstraints( const Elements& elements, AlVec& halfDChi2DAlpha,
+ AlSymMat& halfD2Chi2DAlpha2, std::ostream& logmessage ) const = 0;
+ virtual LHCb::ChiSquare addConstraints( const Elements& inputelements, Al::Equations& equations,
+ std::ostream& logmessage ) const = 0;
+ virtual LHCb::ChiSquare chiSquare( const AlignmentElement& element, bool activeonly = true ) const = 0;
// this returns the survey in the AlignmentFrame of the element, taking correlations into account.
virtual const AlParameters* surveyParameters( const AlignmentElement& element ) const = 0;
};
} // namespace Al
-
-#endif
-
-// LocalWords: TALIGNMENT
diff --git a/Alignment/TAlignment/src/IAlignConstraintTool.h b/Alignment/TAlignment/src/IAlignConstraintTool.h
index e8342a046e9ddfc76091143cc527d17f3058453d..5cd9949f41bbda9e002223a2f066cb9d8b46576a 100755
--- a/Alignment/TAlignment/src/IAlignConstraintTool.h
+++ b/Alignment/TAlignment/src/IAlignConstraintTool.h
@@ -8,8 +8,7 @@
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
-#ifndef TALIGNMENT_IALIGNCONSTRAINTTOOL_H
-#define TALIGNMENT_IALIGNCONSTRAINTTOOL_H
+#pragma once
#include "AlignKernel/AlEquations.h"
#include "AlignKernel/AlSymMat.h"
@@ -18,16 +17,11 @@
#include "IGetElementsToBeAligned.h"
namespace Al {
- /// Static ID object
- static const InterfaceID IID_IAlignConstraintTool( "Al::IAlignConstraintTool", 0, 0 );
+ struct IAlignConstraintTool : public extend_interfaces<IAlgTool> {
+ DeclareInterfaceID( IAlignConstraintTool, 2, 0 );
- class IAlignConstraintTool : virtual public IAlgTool {
- public:
typedef IGetElementsToBeAligned::Elements Elements;
- // Retrieve interface ID
- static const InterfaceID& interfaceID() { return IID_IAlignConstraintTool; }
-
virtual size_t addConstraints( const Elements& elements, const Al::Equations& equations, AlVec& halfDChi2DAlpha,
AlSymMat& halfD2Chi2DAlpha2 ) const = 0;
@@ -35,5 +29,3 @@ namespace Al {
std::ostream& logmessage ) const = 0;
};
} // namespace Al
-
-#endif
diff --git a/Alignment/TAlignment/src/IGetElementsToBeAligned.h b/Alignment/TAlignment/src/IGetElementsToBeAligned.h
index 2d09e3771dc49b42bf8d7c5c838362175b0b2348..26e608b5609051cd15c34469055d75ba8be21223 100755
--- a/Alignment/TAlignment/src/IGetElementsToBeAligned.h
+++ b/Alignment/TAlignment/src/IGetElementsToBeAligned.h
@@ -30,8 +30,6 @@
/// Local
#include "AlignmentElement.h"
-static const InterfaceID IID_IGetElementsToBeAligned( "IGetElementsToBeAligned", 1, 0 );
-
/** @class IGetElementsToBeAligned IGetElementsToBeAligned.h
*
*
@@ -52,13 +50,13 @@ namespace Al {
class Equations;
}
-class IGetElementsToBeAligned : virtual public IAlgTool {
+class IGetElementsToBeAligned : public extend_interfaces<IAlgTool> {
public:
typedef std::vector<const AlignmentElement*> Elements;
// Return the interface ID
- static const InterfaceID& interfaceID() { return IID_IGetElementsToBeAligned; };
+ DeclareInterfaceID( IGetElementsToBeAligned, 2, 0 );
// Return pair of forward begin iter and forward end iter
virtual const Elements& elements() const = 0;
diff --git a/Alignment/TAlignment/src/ITrackResidualTool.h b/Alignment/TAlignment/src/ITrackResidualTool.h
index 98b62452b0c173c1877a209789e66e5b03de3d51..51d8b4df3ac8e6e9cc0799eeae44470da61627f1 100755
--- a/Alignment/TAlignment/src/ITrackResidualTool.h
+++ b/Alignment/TAlignment/src/ITrackResidualTool.h
@@ -17,16 +17,12 @@
namespace LHCb {}
namespace Al {
- /// Static ID object
- static const InterfaceID IID_ITrackResidualTool( "Al::ITrackResidualTool", 0, 0 );
-
- class ITrackResidualTool : virtual public IAlgTool {
- public:
- // Retrieve interface ID
- static const InterfaceID& interfaceID() { return IID_ITrackResidualTool; }
-
+ struct ITrackResidualTool : public extend_interfaces<IAlgTool> {
+ DeclareInterfaceID( ITrackResidualTool, 2, 0 );
// the only method
virtual const Al::TrackResiduals* get( const LHCb::Track& track ) const = 0;
+ // and just another one, temporarily
+ virtual void clearCache() const = 0;
};
} // namespace Al
diff --git a/Alignment/TAlignment/src/IVertexResidualTool.h b/Alignment/TAlignment/src/IVertexResidualTool.h
index 72cc16b3ab90495dffcefe82aa26bc9fec6fa536..1a9895258226bef85fca45cac35135e5e321a9a8 100755
--- a/Alignment/TAlignment/src/IVertexResidualTool.h
+++ b/Alignment/TAlignment/src/IVertexResidualTool.h
@@ -8,8 +8,7 @@
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
\*****************************************************************************/
-#ifndef TALIGNMENT_IVERTEXRESIDUALTOOL_H
-#define TALIGNMENT_IVERTEXRESIDUALTOOL_H
+#pragma once
#include "AlResiduals.h"
#include "Event/Track.h"
@@ -20,21 +19,15 @@ namespace LHCb {
} // namespace LHCb
namespace Al {
- /// Static ID object
- static const InterfaceID IID_IVertexResidualTool( "Al::IVertexResidualTool", 0, 0 );
- class IVertexResidualTool : virtual public IAlgTool {
- public:
- // Retrieve interface ID
- static const InterfaceID& interfaceID() { return IID_IVertexResidualTool; }
+ struct IVertexResidualTool : public extend_interfaces<IAlgTool> {
+ DeclareInterfaceID( IVertexResidualTool, 2, 0 );
// the only method
- typedef std::vector<const LHCb::Track*> TrackContainer;
+ using ReturnType = std::unique_ptr<const Al::MultiTrackResiduals>;
// process a vertex (primary vertex of twoprongvertex)
- virtual const Al::MultiTrackResiduals* get( const LHCb::RecVertex& vertex ) const = 0;
+ virtual ReturnType get( const LHCb::RecVertex& vertex ) const = 0;
// process a particle. this applies a mass constraint.
- virtual const Al::MultiTrackResiduals* get( const LHCb::Particle& p ) const = 0;
+ virtual ReturnType get( const LHCb::Particle& p ) const = 0;
};
} // namespace Al
-
-#endif
diff --git a/Alignment/TAlignment/src/TrackResidualTool.cpp b/Alignment/TAlignment/src/TrackResidualTool.cpp
index 1ccbafb5c5709e4b22ddc87d966672fc90848260..81939639498f7193ed2c73dadc4b054a6576d2d9 100755
--- a/Alignment/TAlignment/src/TrackResidualTool.cpp
+++ b/Alignment/TAlignment/src/TrackResidualTool.cpp
@@ -22,31 +22,33 @@
namespace Al {
- class TrackResidualTool : public GaudiTool, virtual public ITrackResidualTool, virtual public IIncidentListener {
+ class TrackResidualTool : public extends<GaudiTool, ITrackResidualTool>, virtual public IIncidentListener {
public:
// constructor
- TrackResidualTool( const std::string& type, const std::string& name, const IInterface* parent );
- // destructor
- ~TrackResidualTool(){};
+ using extends::extends;
// initialisation
StatusCode initialize() override;
- // finalize
- StatusCode finalize() override;
// incident service handle
void handle( const Incident& incident ) override;
+
+ // clear the cache
+ void clearCache() const override { m_residuals.clear(); }
+
// get the residuals for this track
const Al::TrackResiduals* get( const LHCb::Track& track ) const override;
+ std::unique_ptr<Al::TrackResiduals> compute( const LHCb::Track& track ) const;
+
private:
- const Al::TrackResiduals* compute( const LHCb::Track& track ) const;
- Al::TrackResiduals* create( const LHCb::Track& track, const std::vector<const LHCb::FitNode*>& nodes,
- std::vector<size_t>& residualnodeindices, size_t& refnodeindex ) const;
+ std::unique_ptr<Al::TrackResiduals> create( const LHCb::Track& track,
+ const std::vector<const LHCb::FitNode*>& nodes,
+ std::vector<size_t>& residualnodeindices, size_t& refnodeindex ) const;
private:
- typedef std::map<const LHCb::Track*, const Al::TrackResiduals*> ResidualMap;
- mutable ResidualMap m_residuals;
- bool m_testPosDef;
- ToolHandle<IGetElementsToBeAligned> m_elementTool;
+ using ResidualMap = std::map<const LHCb::Track*, std::unique_ptr<Al::TrackResiduals>>;
+ mutable ResidualMap m_residuals;
+ Gaudi::Property<bool> m_testPosDef{this, "TestPosDef", false};
+ PublicToolHandle<IGetElementsToBeAligned> m_elementTool{this, "ElementTool", "GetElementsToBeAligned"};
};
} // namespace Al
@@ -61,82 +63,62 @@ namespace Al {
DECLARE_COMPONENT( TrackResidualTool )
- TrackResidualTool::TrackResidualTool( const std::string& type, const std::string& name, const IInterface* parent )
- : GaudiTool( type, name, parent ), m_elementTool( "GetElementsToBeAligned" ) {
- // interfaces
- declareInterface<ITrackResidualTool>( this );
- declareProperty( "TestPosDef", m_testPosDef = false );
- }
-
StatusCode TrackResidualTool::initialize() {
StatusCode sc = GaudiTool::initialize();
- if ( sc.isFailure() ) return Error( "Failed to initialize", sc );
- sc = m_elementTool.retrieve();
incSvc()->addListener( this, IncidentType::EndEvent );
return sc;
}
- StatusCode TrackResidualTool::finalize() {
- m_elementTool.release().ignore();
- return GaudiTool::finalize();
- }
-
void TrackResidualTool::handle( const Incident& incident ) {
- if ( IncidentType::EndEvent == incident.type() ) {
- for ( ResidualMap::iterator it = m_residuals.begin(); it != m_residuals.end(); ++it ) delete it->second;
- m_residuals.clear();
- }
+ if ( IncidentType::EndEvent == incident.type() ) clearCache();
}
const Al::TrackResiduals* TrackResidualTool::get( const LHCb::Track& track ) const {
- const Al::TrackResiduals* rc( 0 );
- ResidualMap::const_iterator it = m_residuals.find( &track );
- if ( it != m_residuals.end() )
- rc = it->second;
- else {
- rc = compute( track );
- m_residuals[&track] = rc;
- }
- return rc;
+ // is this the most efficient way to create a cache?
+ auto entry = m_residuals.insert( {&track, std::unique_ptr<Al::TrackResiduals>{}} );
+ if ( entry.second ) entry.first->second = compute( track );
+ return entry.first->second.get();
}
- template <class Matrix>
- Matrix correlationMatrix( const Matrix& M ) {
- Matrix N;
- for ( int i = 0; i < Matrix::kCols; ++i ) {
- N( i, i ) = 1;
- for ( int j = 0; j < i; ++j ) N( i, j ) = M( i, j ) / std::sqrt( M( i, i ) * M( j, j ) );
+ namespace {
+ template <class Matrix>
+ Matrix correlationMatrix( const Matrix& M ) {
+ Matrix N;
+ for ( int i = 0; i < Matrix::kCols; ++i ) {
+ N( i, i ) = 1;
+ for ( int j = 0; j < i; ++j ) N( i, j ) = M( i, j ) / std::sqrt( M( i, i ) * M( j, j ) );
+ }
+ return N;
}
- return N;
- }
- template <class Matrix>
- bool testPosDef( const Matrix& M ) {
- Matrix N = correlationMatrix( M );
- double det;
- /*bool success =*/N.Det( det );
- const double tolerance = 1e-9;
- return /*success &&*/ det + tolerance > 0;
- }
+ template <class Matrix>
+ bool testPosDef( const Matrix& M ) {
+ Matrix N = correlationMatrix( M );
+ double det;
+ /*bool success =*/N.Det( det );
+ const double tolerance = 1e-9;
+ return /*success &&*/ det + tolerance > 0;
+ }
- template <int N>
- ROOT::Math::SMatrix<double, 2 * N, 2 * N, ROOT::Math::MatRepSym<double, 2 * N>>
- constructNxN( const Gaudi::SymMatrix5x5& MAA, const Gaudi::SymMatrix5x5& MBB, const Gaudi::Matrix5x5& MAB ) {
- ROOT::Math::SMatrix<double, 2 * N, 2 * N, ROOT::Math::MatRepSym<double, 2 * N>> M;
- for ( int i = 0; i < N; ++i ) {
- for ( int j = 0; j <= i; ++j ) {
- M( i, j ) = MAA( i, j );
- M( i + N, j + N ) = MBB( i, j );
+ template <int N>
+ ROOT::Math::SMatrix<double, 2 * N, 2 * N, ROOT::Math::MatRepSym<double, 2 * N>>
+ constructNxN( const Gaudi::SymMatrix5x5& MAA, const Gaudi::SymMatrix5x5& MBB, const Gaudi::Matrix5x5& MAB ) {
+ ROOT::Math::SMatrix<double, 2 * N, 2 * N, ROOT::Math::MatRepSym<double, 2 * N>> M;
+ for ( int i = 0; i < N; ++i ) {
+ for ( int j = 0; j <= i; ++j ) {
+ M( i, j ) = MAA( i, j );
+ M( i + N, j + N ) = MBB( i, j );
+ }
+ for ( int j = 0; j < N; ++j ) { M( i, j + N ) = MAB( i, j ); }
}
- for ( int j = 0; j < N; ++j ) { M( i, j + N ) = MAB( i, j ); }
+ return M;
}
- return M;
- }
+ } // namespace
- Al::TrackResiduals* TrackResidualTool::create( const LHCb::Track& track,
- const std::vector<const LHCb::FitNode*>& nodes,
- std::vector<size_t>& residualnodeindices,
- size_t& refnodeindex ) const {
+ std::unique_ptr<Al::TrackResiduals> TrackResidualTool::create( const LHCb::Track& track,
+ const std::vector<const LHCb::FitNode*>& nodes,
+ std::vector<size_t>& residualnodeindices,
+ size_t& refnodeindex ) const {
// first select nodes with a measurement and equip them with an index in the original list
typedef std::pair<const LHCb::FitNode*, size_t> NodeWithIndex;
std::vector<NodeWithIndex> nodeswithindex;
@@ -158,7 +140,7 @@ namespace Al {
// now select the subset of nodes associated to an alignable object
residualnodeindices.clear();
residualnodeindices.reserve( nodeswithindex.size() );
- Al::TrackResiduals* rc = new Al::TrackResiduals( track );
+ auto rc = std::unique_ptr<Al::TrackResiduals>{new Al::TrackResiduals{track}};
rc->m_residuals.reserve( nodeswithindex.size() );
rc->m_state = nodeswithindex.front().first->classicalSmoothedState();
for ( std::vector<NodeWithIndex>::const_iterator it = nodeswithindex.begin(); it != nodeswithindex.end(); ++it ) {
@@ -187,7 +169,7 @@ namespace Al {
return *offdiagcov[l][k];
}
- const Al::TrackResiduals* TrackResidualTool::compute( const LHCb::Track& track ) const {
+ std::unique_ptr<Al::TrackResiduals> TrackResidualTool::compute( const LHCb::Track& track ) const {
debug() << "In TrackResidualTool::Compute" << endmsg;
// this is not stricktly necessary, but let's do it anyway
@@ -276,7 +258,7 @@ namespace Al {
// this is also where things become really slow.
if ( !error && ( error = fr->inError() ) ) warning() << "Error in KalmanFitResult: " << fr->getError() << endmsg;
- Al::TrackResiduals* residuals( 0 );
+ std::unique_ptr<Al::TrackResiduals> residuals;
if ( !error ) {
@@ -338,17 +320,14 @@ namespace Al {
warning() << "Track type = " << track.type() << " " << track.history() << " " << track.flag() << endmsg;
}
}
- if ( error ) {
- delete residuals;
- residuals = 0;
- }
+ if ( error ) residuals.reset();
}
// let's clean up:
for ( size_t irow = 0; irow < numnodes; ++irow )
for ( size_t icol = 0; icol < irow; ++icol ) delete offdiagcov[irow][icol];
- debug() << "End of TrackResidualTool::Compute " << residuals << endmsg;
+ debug() << "End of TrackResidualTool::Compute " << residuals.get() << endmsg;
counter( "Fraction of failures" ) += bool( residuals == 0 );
return residuals;
}
diff --git a/Alignment/TAlignment/src/VertexResidualTool.cpp b/Alignment/TAlignment/src/VertexResidualTool.cpp
index 25cd2e49241810ff5f5e56941b441391fa5e2c0b..6b65a740c71bba983fdbd8c04bfc9207124f3218 100755
--- a/Alignment/TAlignment/src/VertexResidualTool.cpp
+++ b/Alignment/TAlignment/src/VertexResidualTool.cpp
@@ -10,7 +10,7 @@
\*****************************************************************************/
#include "Event/FitNode.h"
#include "GaudiAlg/GaudiTool.h"
-#include "GaudiKernel/IIncidentListener.h"
+#include "GaudiKernel/ServiceHandle.h"
#include "GaudiKernel/ToolHandle.h"
#include "IVertexResidualTool.h"
#include "Kernel/IParticlePropertySvc.h"
@@ -21,7 +21,7 @@
// TAlignment
#include "ITrackResidualTool.h"
-class ITrackExtrapolator;
+struct ITrackExtrapolator;
namespace LHCb {
class RecVertex;
}
@@ -36,53 +36,41 @@ namespace {
} // namespace
namespace Al {
- class ITrackResidualTool;
+ struct ITrackResidualTool;
struct TrackContribution;
- class VertexResidualTool : public GaudiTool, virtual public IVertexResidualTool, virtual public IIncidentListener {
+ class VertexResidualTool : public extends<GaudiTool, IVertexResidualTool> {
public:
//
typedef std::vector<const Al::TrackResiduals*> TrackResidualContainer;
// constructor
- VertexResidualTool( const std::string& type, const std::string& name, const IInterface* parent );
- // destructor
- virtual ~VertexResidualTool(){};
- // initialize
- StatusCode initialize() override;
- // finalize
- StatusCode finalize() override;
- // incident service handle
- void handle( const Incident& incident ) override;
+ using extends::extends;
// process a vertex (primary vertex of twoprongvertex)
- const Al::MultiTrackResiduals* get( const LHCb::RecVertex& vertex ) const override;
+ ReturnType get( const LHCb::RecVertex& vertex ) const override;
// process a particle. this applies a mass constraint.
- const Al::MultiTrackResiduals* get( const LHCb::Particle& particle ) const override;
+ ReturnType get( const LHCb::Particle& particle ) const override;
private:
// create a new MultiTrackResiduals
- const Al::MultiTrackResiduals* compute( const TrackResidualContainer& tracks, const Gaudi::XYZPoint& vertexestimate,
- const MassConstraint* massconstraint ) const;
+ ReturnType compute( const TrackResidualContainer& tracks, const Gaudi::XYZPoint& vertexestimate,
+ const MassConstraint* massconstraint = nullptr ) const;
// extrapolate the state in trackresiduals to position z
StatusCode extrapolate( const Al::TrackResiduals& trackin, double z, TrackContribution& trackout ) const;
private:
- ToolHandle<ITrackResidualTool> m_trackresidualtool;
- ToolHandle<ITrackExtrapolator> m_extrapolator;
- LHCb::IParticlePropertySvc* m_propertysvc;
- double m_chiSquarePerDofCut;
- double m_maxMassConstraintChi2;
- bool m_computeCorrelations;
- size_t m_maxHitsPerTrackForCorrelations;
- typedef std::vector<const Al::MultiTrackResiduals*> ResidualContainer;
- mutable ResidualContainer m_residuals;
- bool m_debugLevel;
+ PublicToolHandle<ITrackResidualTool> m_trackresidualtool{this, "TrackResidualTool", "Al::TrackResidualTool"};
+ PublicToolHandle<ITrackExtrapolator> m_extrapolator{this, "Extrapolator", "TrackMasterExtrapolator"};
+ ServiceHandle<LHCb::IParticlePropertySvc> m_propertysvc{this, "ParticlePropertyService",
+ "LHCb::ParticlePropertySvc"};
+ Gaudi::Property<double> m_chiSquarePerDofCut{this, "MaxVertexChi2PerDoF", 10};
+ Gaudi::Property<double> m_maxMassConstraintChi2{this, "MaxMassConstraintChi2", 16};
+ Gaudi::Property<int> m_maxHitsPerTrackForCorrelations{this, "MaxHitsPerTrackForCorrelations", -1};
};
} // namespace Al
/////////////////////////////////////////////////////////////////////////////////////////////////////////
-#include "GaudiKernel/IIncidentSvc.h"
#include <algorithm>
#include <iostream>
@@ -100,54 +88,11 @@ namespace Al {
DECLARE_COMPONENT( VertexResidualTool )
- VertexResidualTool::VertexResidualTool( const std::string& type, const std::string& name, const IInterface* parent )
- : GaudiTool( type, name, parent )
- , m_trackresidualtool( "Al::TrackResidualTool" )
- , // important: use the toolsvc, because of caching!
- m_extrapolator( "TrackMasterExtrapolator" )
- , m_propertysvc( 0 )
- , m_maxHitsPerTrackForCorrelations( 9999 ) {
- // interfaces
- declareInterface<IVertexResidualTool>( this );
- declareProperty( "TrackResidualTool", m_trackresidualtool );
- declareProperty( "Extrapolator", m_extrapolator );
- declareProperty( "MaxHitsPerTrackForCorrelations", m_maxHitsPerTrackForCorrelations );
- declareProperty( "MaxVertexChi2PerDoF", m_chiSquarePerDofCut = 10 );
- declareProperty( "MaxMassConstraintChi2", m_maxMassConstraintChi2 = 16 );
- }
-
- StatusCode VertexResidualTool::initialize() {
- StatusCode sc = GaudiTool::initialize();
- if ( sc.isFailure() ) { return Error( "Failed to initialize", sc ); }
- m_trackresidualtool.retrieve().ignore();
- m_extrapolator.retrieve().ignore();
- incSvc()->addListener( this, IncidentType::EndEvent );
- m_propertysvc = svc<LHCb::IParticlePropertySvc>( "LHCb::ParticlePropertySvc", true );
- m_debugLevel = msgLevel( MSG::DEBUG );
- return sc;
- }
-
- StatusCode VertexResidualTool::finalize() {
- bool success = m_trackresidualtool.release().isSuccess() && m_extrapolator.release().isSuccess() &&
- GaudiTool::finalize().isSuccess();
-
- return success ? StatusCode::SUCCESS : StatusCode::FAILURE;
- }
-
- void VertexResidualTool::handle( const Incident& incident ) {
- if ( IncidentType::EndEvent == incident.type() ) {
- for ( ResidualContainer::iterator it = m_residuals.begin(); it != m_residuals.end(); ++it ) delete *it;
- m_residuals.clear();
- }
- }
-
- const Al::MultiTrackResiduals* VertexResidualTool::get( const LHCb::RecVertex& vertex ) const {
+ VertexResidualTool::ReturnType VertexResidualTool::get( const LHCb::RecVertex& vertex ) const {
// loop over the list of vertices, collect tracks in the vertex
- const Al::MultiTrackResiduals* rc( 0 );
- TrackResidualContainer trackresiduals;
- for ( SmartRefVector<LHCb::Track>::const_iterator itrack = vertex.tracks().begin(); itrack != vertex.tracks().end();
- ++itrack ) {
- const Al::TrackResiduals* trackres = m_trackresidualtool->get( **itrack );
+ TrackResidualContainer trackresiduals;
+ for ( const auto& track : vertex.tracks() ) {
+ const Al::TrackResiduals* trackres = m_trackresidualtool->get( *track );
if ( trackres ) {
trackresiduals.push_back( trackres );
} else {
@@ -155,20 +100,14 @@ namespace Al {
assert( 0 );
}
}
-
- const MassConstraint* mconstraint( 0 );
- rc = compute( trackresiduals, vertex.position(), mconstraint );
- if ( rc ) m_residuals.push_back( rc );
-
- delete mconstraint;
- return rc;
+ return std::unique_ptr<const Al::MultiTrackResiduals>( compute( trackresiduals, vertex.position() ) );
}
- const Al::MultiTrackResiduals* VertexResidualTool::get( const LHCb::Particle& p ) const {
+ VertexResidualTool::ReturnType VertexResidualTool::get( const LHCb::Particle& p ) const {
// loop over the list of vertices, collect tracks in the vertex
- const Al::MultiTrackResiduals* rc( 0 );
- TrackResidualContainer trackresiduals;
- MassConstraint mconstraint;
+ ReturnType rc;
+ TrackResidualContainer trackresiduals;
+ MassConstraint mconstraint;
mconstraint.mothermass = m_propertysvc->find( p.particleID() )->mass();
mconstraint.motherwidth = m_propertysvc->find( p.particleID() )->width();
// info() << "Mothermass: " << mconstraint.mothermass << "\t width: " << mconstraint.motherwidth << endmsg ;
@@ -200,10 +139,7 @@ namespace Al {
}
}
- if ( success ) {
- rc = compute( trackresiduals, p.referencePoint(), &mconstraint );
- if ( rc ) m_residuals.push_back( rc );
- }
+ if ( success ) rc = compute( trackresiduals, p.referencePoint(), &mconstraint );
return rc;
}
@@ -217,11 +153,12 @@ namespace Al {
size_t offset;
};
- const Al::MultiTrackResiduals* VertexResidualTool::compute( const TrackResidualContainer& tracks,
+ VertexResidualTool::ReturnType VertexResidualTool::compute( const TrackResidualContainer& tracks,
const Gaudi::XYZPoint& vertexestimate,
const MassConstraint* mconstraint ) const {
- Al::MultiTrackResiduals* rc( 0 );
- bool success = true;
+ Al::MultiTrackResiduals* rc{nullptr};
+ bool success = true;
+ bool m_debugLevel = msgLevel( MSG::DEBUG );
// first we need to propagate all tracks to the vertex and keep
// track of the correlation between the state at the vertex and
@@ -397,11 +334,12 @@ namespace Al {
if ( m_maxHitsPerTrackForCorrelations > 0 && !computationerror )
for ( size_t j = 0; j < i && !computationerror; ++j ) {
size_t joffset = states[j].offset;
- size_t maxirow = std::min( states[i].trackresiduals->size(), m_maxHitsPerTrackForCorrelations );
+ size_t maxirow = std::min( states[i].trackresiduals->size(), size_t( m_maxHitsPerTrackForCorrelations ) );
for ( size_t irow = 0; irow < maxirow; ++irow ) {
// store this intermediate matrix too save some time
Gaudi::TrackProjectionMatrix A = states[i].dResidualdState[irow] * vertex.stateCovariance( i, j );
- size_t maxjcol = std::min( states[j].trackresiduals->size(), m_maxHitsPerTrackForCorrelations );
+ size_t maxjcol =
+ std::min( states[j].trackresiduals->size(), size_t( m_maxHitsPerTrackForCorrelations ) );
for ( size_t jcol = 0; jcol < maxjcol; ++jcol ) {
double deltaHCH = ( A * ROOT::Math::Transpose( states[j].dResidualdState[jcol] ) )( 0, 0 );
rc->m_HCHElements.push_back( Al::CovarianceElement( irow + ioffset, jcol + joffset, deltaHCH ) );
@@ -419,7 +357,7 @@ namespace Al {
if ( computationerror ) {
delete rc;
- rc = 0;
+ rc = nullptr;
warning() << "VertexResidualTool::compute failed" << endmsg;
}
} else {
@@ -429,7 +367,7 @@ namespace Al {
} else {
warning() << "not enough tracks for vertex anymore" << endmsg;
}
- return rc;
+ return ReturnType{rc};
}
StatusCode VertexResidualTool::extrapolate( const Al::TrackResiduals& track, double z, TrackContribution& tc ) const {
diff --git a/Alignment/VeloAlignment/src/PVHisto.cpp b/Alignment/VeloAlignment/src/PVHisto.cpp
index 800e2a46ae490cc5536072acece67a5fd2da372b..f2b55e95d5f173565c973df5dbdf77c283f5e627 100755
--- a/Alignment/VeloAlignment/src/PVHisto.cpp
+++ b/Alignment/VeloAlignment/src/PVHisto.cpp
@@ -29,7 +29,6 @@ PVHisto::PVHisto() {}
PVHisto::PVHisto( double min, double max, double width )
: m_Min( min )
- , m_Max( max )
, m_BinWidth( width )
, m_Bins( int( 0.99 + ( max - min ) / width ) )
, m_Histo( m_Bins, 0 )
@@ -37,12 +36,7 @@ PVHisto::PVHisto( double min, double max, double width )
m_MaxBin = m_Histo.begin();
}
PVHisto::PVHisto( double min, double max, int bins )
- : m_Min( min )
- , m_Max( max )
- , m_BinWidth( ( max - min ) / bins )
- , m_Bins( bins )
- , m_Histo( bins, 0 )
- , m_MaxVal( 0 ) {
+ : m_Min( min ), m_BinWidth( ( max - min ) / bins ), m_Bins( bins ), m_Histo( bins, 0 ), m_MaxVal( 0 ) {
m_MaxBin = m_Histo.begin();
}
//=============================================================================
diff --git a/Alignment/VeloAlignment/src/PVHisto.h b/Alignment/VeloAlignment/src/PVHisto.h
index 46a0d7f5a3d0621bf56d405f6b46503f3a808287..fb6619fc252e3654a5ae01be1e0d0053faa2b534 100755
--- a/Alignment/VeloAlignment/src/PVHisto.h
+++ b/Alignment/VeloAlignment/src/PVHisto.h
@@ -48,7 +48,7 @@ public:
protected:
private:
- double m_Min, m_Max;
+ double m_Min;
double m_BinWidth;
int m_Bins;
std::vector<int> m_Histo;
diff --git a/Alignment/VeloAlignment/src/VAlign.h b/Alignment/VeloAlignment/src/VAlign.h
index 4514d5aad75b09c88a5f55fbe584ed473222f34c..23101af01e54010536ab61eeac1cff22228fb5b6 100755
--- a/Alignment/VeloAlignment/src/VAlign.h
+++ b/Alignment/VeloAlignment/src/VAlign.h
@@ -104,8 +104,6 @@ public:
protected:
private:
- const LHCb::RecVertices* pvcontainer;
-
// TrackStore Variables
MilleConfig* my_Config;
@@ -173,15 +171,10 @@ private:
std::vector<double> m_sigmab;
std::vector<bool> m_constrainb;
std::vector<double> m_residual_cutb;
- double m_z_range;
unsigned int m_PV_trackmin;
- double m_z_sigma;
- double m_IPmax;
- double m_TrIPmax;
double m_starfactr;
unsigned int m_maxtrack;
- double m_maxChi2;
bool m_moni_constants;
bool m_moni_PV;
bool m_moni_overlaps;
@@ -198,12 +191,8 @@ private:
std::vector<double> mis_pull;
// VELO geometry info
-
DeVelo* m_velo;
- const DeVeloRType* rDet;
- const DeVeloPhiType* phiDet;
-
double VELOmap[42];
double zmoy_R, zmoy_L, szmoy_R, szmoy_L;
diff --git a/Calibration/Pi0Calibration/Pi0Calibration/MMapVector.h b/Calibration/Pi0Calibration/Pi0Calibration/MMapVector.h
index 90112d16719768ecb9471406946bd2b287037f5e..98abfd714d3176a4da304899c6166503e29fdd96 100644
--- a/Calibration/Pi0Calibration/Pi0Calibration/MMapVector.h
+++ b/Calibration/Pi0Calibration/Pi0Calibration/MMapVector.h
@@ -553,7 +553,7 @@ public:
/// clear the vector
void clear() noexcept( noexcept( destroy_range( begin(), end() ) ) ) {
destroy_range( begin(), end() );
- m_size -= m_size;
+ m_size = 0;
}
/// resize the vector (default-construct new elements if needed)
diff --git a/Calibration/Pi0Calibration/src/Pi0MassFitter.cpp b/Calibration/Pi0Calibration/src/Pi0MassFitter.cpp
index b58edc11aeb1e3295e92cd7d3ce7134fb89fd0c9..1e53dc5a0f19f7ccc8767d7b0bf76dac3ceee1af 100644
--- a/Calibration/Pi0Calibration/src/Pi0MassFitter.cpp
+++ b/Calibration/Pi0Calibration/src/Pi0MassFitter.cpp
@@ -62,13 +62,13 @@ using namespace boost::filesystem;
using namespace RooFit;
using namespace Calibration::Pi0Calibration;
-const unsigned int BitsCol = 6;
-const unsigned int BitsRow = 6;
-const unsigned int BitsArea = 2;
-const unsigned int ShiftCol = 0;
-const unsigned int ShiftRow = ShiftCol + BitsCol;
-const unsigned int ShiftArea = ShiftRow + BitsRow;
-const unsigned int MaskArea = ( ( ( (unsigned int)1 ) << BitsArea ) - 1 ) << ShiftArea;
+// const unsigned int BitsCol = 6;
+// const unsigned int BitsRow = 6;
+// const unsigned int BitsArea = 2;
+// const unsigned int ShiftCol = 0;
+// const unsigned int ShiftRow = ShiftCol + BitsCol;
+// const unsigned int ShiftArea = ShiftRow + BitsRow;
+// const unsigned int MaskArea = ( ( ( (unsigned int)1 ) << BitsArea ) - 1 ) << ShiftArea;
Double_t Pi0MassFitter::polyn( Double_t* x, Double_t* par ) { return par[0] + par[1] * x[0] + par[2] * x[0] * x[0]; }
Double_t Pi0MassFitter::polynSB( Double_t* x, Double_t* par ) {