ATLAS Offline Software
Loading...
Searching...
No Matches
Trk::MatrixTool Class Reference

#include <MatrixTool.h>

Inheritance diagram for Trk::MatrixTool:
Collaboration diagram for Trk::MatrixTool:

Public Types

enum  SolveOption {
  NONE = 0 , SOLVE = 1 , SOLVE_FAST = 2 , DIRECT_SOLVE = 3 ,
  DIRECT_SOLVE_FAST = 4 , DIRECT_SOLVE_CLUSTER = 5 , SOLVE_ROOT = 6 , SOLVE_CLHEP = 7
}

Public Member Functions

 MatrixTool (const std::string &type, const std::string &name, const IInterface *parent)
 Constructor.
virtual ~MatrixTool ()
 Virtual destructor.
StatusCode initialize ()
 initialize
StatusCode finalize ()
 initialize
StatusCode allocateMatrix (int nDoF=0)
 allocates memory for big matrix and big vector
void prepareBinaryFiles (int solveOption)
 reads/writes matrix entries from/to binary files as necessary
void addFirstDerivatives (AlVec *vector)
 adds first derivative to vector
void addFirstDerivatives (std::list< int, double > &derivatives)
 adds first derivative to vector for only some entries
void addFirstDerivative (int irow, double firstderiv)
void addSecondDerivatives (AlSymMatBase *matrix)
 adds second derivatives to matrix
void addSecondDerivatives (std::list< std::pair< int, int >, double > &derivatives)
 adds first derivative to vector for only some entries
void addSecondDerivative (int irow, int icol, double secondderiv)
bool accumulateFromFiles ()
 accumulates derivates from files.
bool accumulateFromBinaries ()
 accumulates derivates from binary files
int solve ()
 solves for alignment parameters
void storeInTFile (const TString &filename)
 Store Files in a tfile.
bool accumulateFromTFiles ()
 Store Files in a tfile.
void printModuleSolution (std::ostream &os, const AlignModule *module, const CLHEP::HepSymMatrix *cov) const
 namespace { class RestoreIOSFlags { public: RestoreIOSFlags (std::ostream &os) : m_os(&os), m_flags(m_os->flags()), m_precision(m_os->precision()) {} ~RestoreIOSFlags() { m_os->flags(m_flags); m_os->precision(m_precision); } private: std::ostream *m_os; std::ios_base::fmtflags m_flags; std::streamsize m_precision; }; }
void printGlobalSolution (std::ostream &os, const CLHEP::HepSymMatrix *cov)
void printGlobalSolution (std::ostream &os, const TMatrixDSym *cov)
ServiceHandle< StoreGateSvc > & evtStore ()
 The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
const ServiceHandle< StoreGateSvc > & detStore () const
 The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.
virtual StatusCode sysInitialize () override
 Perform system initialization for an algorithm.
virtual StatusCode sysStart () override
 Handle START transition.
virtual std::vector< Gaudi::DataHandle * > inputHandles () const override
 Return this algorithm's input handles.
virtual std::vector< Gaudi::DataHandle * > outputHandles () const override
 Return this algorithm's output handles.
Gaudi::Details::PropertyBase & declareProperty (Gaudi::Property< T, V, H > &t)
void updateVHKA (Gaudi::Details::PropertyBase &)
MsgStream & msg () const
bool msgLvl (const MSG::Level lvl) const
void addModule (int alignModuleIndex, int nAlignParam)
int entryNumber (int alignModuleIndex)
void setNHits (int n)
 set module identifier
int nHits () const
void setNTracks (int n)
 set number of tracks
int nTracks () const
void setNMeasurements (int n)
 set number of measurements
int nMeasurements () const
virtual void setLogStream (std::ostream *os)
 sets the output stream for the logfile

Static Public Member Functions

static const InterfaceID & interfaceID ()
 Retrieve interface ID.

Protected Member Functions

void renounceArray (SG::VarHandleKeyArray &handlesArray)
 remove all handles from I/O resolution
std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > renounce (T &h)
void extraDeps_update_handler (Gaudi::Details::PropertyBase &ExtraDeps)
 Add StoreName to extra input/output deps as needed.

Protected Attributes

std::ostream * m_logStream
 logfile output stream
int m_nHits
int m_nTracks
int m_nMeasurements

Private Types

typedef ServiceHandle< StoreGateSvcStoreGateSvc_t

Private Member Functions

int solveROOT ()
int solveCLHEP ()
int solveLapack ()
int solveSparseEigen ()
int solveLocal ()
StatusCode spuriousRemoval ()
void postSolvingLapack (AlVec *dChi2, AlSymMat *d2Chi2, AlVec &w, AlMat &z, int size)
void writeHitmap ()
void readHitmaps ()
Gaudi::Details::PropertyBase & declareGaudiProperty (Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
 specialization for handling Gaudi::Property<SG::VarHandleKey>

Static Private Member Functions

static int fillVecMods ()

Private Attributes

PublicToolHandle< IAlignModuleToolm_alignModuleTool
AlSymMatBasem_bigmatrix = nullptr
 matrix to contain second derivative terms to be used for alignment
AlVecm_bigvector = nullptr
 vector to contain first derivative terms to be used for alignment
Gaudi::Property< bool > m_useSparse {this, "UseSparse", false}
 flag to use AlSpaMat for the big matrix (default is AlSymMat)
Gaudi::Property< bool > m_diagonalize
Gaudi::Property< double > m_eigenvaluethreshold {this, "EigenvalueThreshold", 0., "cut on the minimum eigenvalue"}
Gaudi::Property< int > m_solveOption {this, "SolveOption", NONE, "solving option"}
Gaudi::Property< int > m_modcut
Gaudi::Property< int > m_minNumHits
Gaudi::Property< int > m_minNumTrks
Gaudi::Property< float > m_pullcut
Gaudi::Property< float > m_eigenvalueStep
Gaudi::Property< float > m_Align_db_step
Gaudi::Property< bool > m_calDet
Gaudi::Property< bool > m_wSqMatrix
Gaudi::Property< bool > m_writeMat
Gaudi::Property< bool > m_writeMatTxt
Gaudi::Property< bool > m_writeEigenMat
Gaudi::Property< bool > m_writeEigenMatTxt
Gaudi::Property< bool > m_writeModuleNames
Gaudi::Property< bool > m_writeHitmap
Gaudi::Property< bool > m_writeHitmapTxt
Gaudi::Property< bool > m_readHitmaps
Gaudi::Property< bool > m_writeTFile
Gaudi::Property< bool > m_readTFiles
Gaudi::Property< bool > m_runLocal
double m_scale = -1.
 scale for big matrix and vector normalization
Gaudi::Property< bool > m_scaleMatrix
Gaudi::Property< double > m_softEigenmodeCut
Gaudi::Property< double > m_removeSpurious
Gaudi::Property< double > m_calculateFullCovariance
Gaudi::Property< std::string > m_pathbin
Gaudi::Property< std::string > m_pathtxt
Gaudi::Property< std::string > m_prefixName
Gaudi::Property< std::string > m_tfileName {this, "TFileName", "AlignmentTFile.root", "prefix string to filenames"}
Gaudi::Property< std::string > m_scalaMatName {this, "ScalapackMatrixName", "eigenvectors.bin", "Scalapack matrix name"}
Gaudi::Property< std::string > m_scalaVecName {this, "ScalapackVectorName", "eigenvalues.bin", "Scalapack vector name"}
Gaudi::Property< std::vector< std::string > > m_inputMatrixFiles
Gaudi::Property< std::vector< std::string > > m_inputVectorFiles
Gaudi::Property< std::vector< std::string > > m_inputHitmapFiles
Gaudi::Property< std::vector< std::string > > m_inputTFiles
std::vector< int > m_activeIndices {}
 vector of indices which pass the min-hits cut
int m_aNDoF = 0
 number of active DoF (size of m_activeIndices)
Gaudi::Property< int > m_maxReadErrors
Gaudi::Property< bool > m_AlignIBLbutNotPixel {this, "AlignIBLbutNotPixel", false}
Gaudi::Property< bool > m_AlignPixelbutNotIBL {this, "AlignPixelbutNotIBL", false}
Gaudi::Property< bool > m_DeactivateSCT_ECA_LastDisk {this, "DeactivateSCT_ECA_LastDisk", false}
Gaudi::Property< bool > m_Remove_Pixel_Tx {this, "Remove_Pixel_Tx", false}
Gaudi::Property< bool > m_Remove_Pixel_Ty {this, "Remove_Pixel_Ty", false}
Gaudi::Property< bool > m_Remove_Pixel_Tz {this, "Remove_Pixel_Tz", false}
Gaudi::Property< bool > m_Remove_Pixel_Rx {this, "Remove_Pixel_Rx", false}
Gaudi::Property< bool > m_Remove_Pixel_Ry {this, "Remove_Pixel_Ry", false}
Gaudi::Property< bool > m_Remove_Pixel_Rz {this, "Remove_Pixel_Rz", false}
Gaudi::Property< bool > m_Remove_IBL_Tx {this, "Remove_IBL_Tx", false}
Gaudi::Property< bool > m_Remove_IBL_Ty {this, "Remove_IBL_Ty", false}
Gaudi::Property< bool > m_Remove_IBL_Tz {this, "Remove_IBL_Tz", false}
Gaudi::Property< bool > m_Remove_IBL_Rx {this, "Remove_IBL_Rx", false}
Gaudi::Property< bool > m_Remove_IBL_Ry {this, "Remove_IBL_Ry", false}
Gaudi::Property< bool > m_Remove_IBL_Rz {this, "Remove_IBL_Rz", false}
StoreGateSvc_t m_evtStore
 Pointer to StoreGate (event store by default).
StoreGateSvc_t m_detStore
 Pointer to StoreGate (detector store by default).
std::vector< SG::VarHandleKeyArray * > m_vhka
bool m_varHandleArraysDeclared
std::map< int, int > m_alignModuleMap
int m_nentries

Detailed Description

Definition at line 55 of file MatrixTool.h.

Member Typedef Documentation

◆ StoreGateSvc_t

typedef ServiceHandle<StoreGateSvc> AthCommonDataStore< AthCommonMsg< AlgTool > >::StoreGateSvc_t
privateinherited

Definition at line 388 of file AthCommonDataStore.h.

Member Enumeration Documentation

◆ SolveOption

Enumerator
NONE 

not solve in any case (to be used when ipc)

SOLVE 

solving after data accumulation (LAPACK)

SOLVE_FAST 

Fast (Eigen method) solving after data accumulation.

DIRECT_SOLVE 

direct solving (LAPACK), already available matrix & vector

DIRECT_SOLVE_FAST 

direct Fast (Eigen method) solving, already available matrix & vector

DIRECT_SOLVE_CLUSTER 

computation of alignment parameters from SCALAPAK already solved matrix

SOLVE_ROOT 

computation using ROOT

SOLVE_CLHEP 

computation using CLHEP

Definition at line 58 of file MatrixTool.h.

58 {
59 NONE = 0,
60 SOLVE = 1,
61 SOLVE_FAST = 2,
62 DIRECT_SOLVE = 3,
65 SOLVE_ROOT = 6,
66 SOLVE_CLHEP = 7
67 };
@ SOLVE_FAST
Fast (Eigen method) solving after data accumulation.
Definition MatrixTool.h:61
@ SOLVE_ROOT
computation using ROOT
Definition MatrixTool.h:65
@ SOLVE
solving after data accumulation (LAPACK)
Definition MatrixTool.h:60
@ DIRECT_SOLVE_FAST
direct Fast (Eigen method) solving, already available matrix & vector
Definition MatrixTool.h:63
@ DIRECT_SOLVE
direct solving (LAPACK), already available matrix & vector
Definition MatrixTool.h:62
@ DIRECT_SOLVE_CLUSTER
computation of alignment parameters from SCALAPAK already solved matrix
Definition MatrixTool.h:64
@ NONE
not solve in any case (to be used when ipc)
Definition MatrixTool.h:59
@ SOLVE_CLHEP
computation using CLHEP
Definition MatrixTool.h:66

Constructor & Destructor Documentation

◆ MatrixTool()

MatrixTool::MatrixTool ( const std::string & type,
const std::string & name,
const IInterface * parent )

Constructor.

Definition at line 54 of file MatrixTool.cxx.

56 : IMatrixTool()
57 , AthAlgTool(type,name,parent)
58 {
59 declareInterface<IMatrixTool>(this);
60 }
AthAlgTool()
Default constructor:
IMatrixTool()
constructor

◆ ~MatrixTool()

MatrixTool::~MatrixTool ( )
virtual

Virtual destructor.

Definition at line 63 of file MatrixTool.cxx.

64 {
65 delete m_bigmatrix;
66 delete m_bigvector;
67 }
AlVec * m_bigvector
vector to contain first derivative terms to be used for alignment
Definition MatrixTool.h:149
AlSymMatBase * m_bigmatrix
matrix to contain second derivative terms to be used for alignment
Definition MatrixTool.h:146

Member Function Documentation

◆ accumulateFromBinaries()

bool MatrixTool::accumulateFromBinaries ( )

accumulates derivates from binary files

Definition at line 631 of file MatrixTool.cxx.

632 {
633
634
635 DataVector<AlignPar>* alignParList = m_alignModuleTool->alignParList1D();
636 int nDoF=alignParList->size();
637
638 std::map<int,unsigned long long> modIndexMap;
639 float dummyVersion(0.);
640 double totalscale=0.;
641 for (int ivec=0;ivec<(int)m_inputVectorFiles.size();ivec++) {
642
643 ATH_MSG_DEBUG("Reading vector "<<ivec<<" from file "<<m_inputVectorFiles[ivec]);
644
645 AlVec newVector(nDoF);
646 std::map<int,unsigned long long> newModIndexMap;
647 newVector.SetPathBin(m_pathbin.value()+m_prefixName.value());
648 newVector.SetPathTxt(m_pathtxt.value()+m_prefixName.value());
649 double scale=0;
650 StatusCode sc = newVector.ReadPartial(m_inputVectorFiles[ivec],scale,newModIndexMap,dummyVersion);
651 totalscale += scale;
652 if (sc==StatusCode::FAILURE) {
653 msg(MSG::FATAL)<<"Problem reading vector from "<<m_inputVectorFiles[ivec]<<endmsg;
654 return false;
655 }
656 if (newVector.size()!=m_bigvector->size()) {
657 msg(MSG::FATAL) <<"vector wrong size! newVector size "<<newVector.size()
658 <<", bigvector size "<<m_bigvector->size()<<endmsg;
659 return false;
660 }
661
662 // check modIndexMaps to make sure they are the same
663 if (ivec==0)
664 modIndexMap = std::move(newModIndexMap);
665 else if (modIndexMap!=newModIndexMap) {
666 msg(MSG::FATAL)<<"module index maps don't agree!"<<endmsg;
667 return false;
668 }
669 if (ivec>0)
670 *m_bigvector += newVector;
671 else
672 *m_bigvector = std::move(newVector);
673 }
674
675 m_scale = totalscale;
676
677 AlSymMat * symBigMatrix=dynamic_cast<AlSymMat*>(m_bigmatrix);
678 AlSpaMat * spaBigMatrix=dynamic_cast<AlSpaMat*>(m_bigmatrix);
679
680
681 for (int imat=0;imat<(int)m_inputMatrixFiles.size();imat++) {
682 ATH_MSG_DEBUG("Reading matrix "<<imat<<" from file "<<m_inputMatrixFiles[imat]);
683
684 // create new matrix to read data from current file
685 int nDoF=modIndexMap.size();
686 bool triang;
688 if (symBigMatrix) {
689 AlSymMat newMatrix(nDoF);
690 sc = newMatrix.Read(m_inputMatrixFiles[imat],nDoF,triang,dummyVersion);
691 if (sc==StatusCode::SUCCESS)
692 *symBigMatrix += newMatrix;
693 }
694 else {
695 if (!spaBigMatrix) {
696 throw std::logic_error("Unhandled matrix type");
697 }
698
699 AlSpaMat newMatrix(nDoF);
700 sc = newMatrix.Read(m_inputMatrixFiles[imat],nDoF,triang,dummyVersion);
701
702 if (sc==StatusCode::SUCCESS) {
703 if (imat>0)
704 *spaBigMatrix += newMatrix;
705 else
706 *spaBigMatrix = newMatrix;
707 }
708 }
709
710 if (sc==StatusCode::FAILURE) {
711 msg(MSG::FATAL)<<"problem reading matrix from "<<m_inputMatrixFiles[imat]<<endmsg;
712 return false;
713 }
714
715 if (!m_useSparse && triang==m_wSqMatrix) {
716 ATH_MSG_WARNING("matrix not expected format! Changing m_wSqMatrix to "<<!triang);
717 m_wSqMatrix=!triang;
718 }
719
720 }
721
722 // accumulate hitmap from hitmap files
723 if(m_readHitmaps)
724 readHitmaps();
725
726 return true;
727 }
#define endmsg
#define ATH_MSG_WARNING(x)
#define ATH_MSG_DEBUG(x)
static Double_t sc
MsgStream & msg() const
size_type size() const noexcept
Returns the number of elements in the collection.
Gaudi::Property< std::string > m_prefixName
Definition MatrixTool.h:222
Gaudi::Property< bool > m_readHitmaps
Definition MatrixTool.h:193
Gaudi::Property< bool > m_wSqMatrix
Definition MatrixTool.h:176
double m_scale
scale for big matrix and vector normalization
Definition MatrixTool.h:204
Gaudi::Property< std::vector< std::string > > m_inputVectorFiles
Definition MatrixTool.h:237
Gaudi::Property< std::vector< std::string > > m_inputMatrixFiles
Definition MatrixTool.h:234
PublicToolHandle< IAlignModuleTool > m_alignModuleTool
Definition MatrixTool.h:142
Gaudi::Property< std::string > m_pathbin
Definition MatrixTool.h:218
Gaudi::Property< bool > m_useSparse
flag to use AlSpaMat for the big matrix (default is AlSymMat)
Definition MatrixTool.h:152
Gaudi::Property< std::string > m_pathtxt
Definition MatrixTool.h:220
vec_fb< CxxUtils::bit_int_t< sizeof(T) *8 >, N > ivec
Definition vec_fb.h:52
::StatusCode StatusCode
StatusCode definition for legacy code.

◆ accumulateFromFiles()

bool MatrixTool::accumulateFromFiles ( )
virtual

accumulates derivates from files.

Flag decides if it is binary or TFiles

Implements Trk::IMatrixTool.

Definition at line 618 of file MatrixTool.cxx.

619 {
620
621 if(m_readTFiles){
622 ATH_MSG_INFO("Info to obtained from from TFiles");
623 return accumulateFromTFiles();
624 }else{
625 ATH_MSG_INFO("Info to obtained from from Binary files");
626 return accumulateFromBinaries();
627 }
628 }
#define ATH_MSG_INFO(x)
bool accumulateFromTFiles()
Store Files in a tfile.
bool accumulateFromBinaries()
accumulates derivates from binary files
Gaudi::Property< bool > m_readTFiles
Definition MatrixTool.h:198

◆ accumulateFromTFiles()

bool MatrixTool::accumulateFromTFiles ( )

Store Files in a tfile.

Definition at line 825 of file MatrixTool.cxx.

826 {
827 DataVector<AlignPar>* alignParList = m_alignModuleTool->alignParList1D();
828 int nDoF=alignParList->size();
829 ATH_MSG_DEBUG("OPENING TFILES");
830
831 std::map<int,unsigned long long> modIndexMap;
832 std::map<int,unsigned long long> DoFMap;
833 double totalscale=0.;
834
835 AlSymMat * symBigMatrix=dynamic_cast<AlSymMat*>(m_bigmatrix);
836 AlSpaMat * spaBigMatrix=dynamic_cast<AlSpaMat*>(m_bigmatrix);
837 //TMatrixDSparse *accumMatrix(0);
838 AlSpaMat *accumMatrix = nullptr;
839
840 const AlignModuleList * moduleList = m_alignModuleTool->alignModules1D();
841 int nModules = moduleList->size();
842
843 TVectorD TotalHits(nModules);
844 TVectorD TotalTracks(nModules);
845
846 int numberOfReadErrors = 0;
847
848 struct rusage myusage{};
849 int itworked = getrusage(RUSAGE_SELF,&myusage);
850 if(itworked == 0)//note: rusage returns zero if it succeeds!
851 ATH_MSG_DEBUG("ItWorked");
852
853 long intialMemUse = myusage.ru_maxrss;
854
855 for (int ifile = 0; ifile < (int)m_inputTFiles.size(); ifile++) {
856 if (numberOfReadErrors > m_maxReadErrors){
857 msg(MSG::FATAL) << " number of errors when reading the TFiles already exceed " << m_maxReadErrors << endmsg;
858 delete accumMatrix;
859 return false;
860 }
861
862 ATH_MSG_DEBUG("Reading File number " << ifile << ", " << m_inputTFiles[ifile]);
863
864 itworked = getrusage(RUSAGE_SELF,&myusage);
865 if (itworked ==0){
866 ATH_MSG_DEBUG("Memory usage [MB], total " << myusage.ru_maxrss/1024 << ", increase " << (myusage.ru_maxrss-intialMemUse)/1024);
867 }
868 TFile* myFile = TFile::Open(m_inputTFiles[ifile].c_str());
869
870 if ( myFile->IsZombie() || !(myFile->IsOpen()) ) {
871 ++numberOfReadErrors;
872 ATH_MSG_ERROR( " Problem reading TFile " << m_inputTFiles[ifile] );
873 continue;
874 }
875
876 std::map<int,unsigned long long> newModIndexMap;
877
878 TVectorD* myModuleIDs;
879 myModuleIDs = (TVectorD*)myFile->Get("ModuleID");
880 if( !myModuleIDs ){
881 ++numberOfReadErrors;
882 ATH_MSG_ERROR("Modules ID not read!!!");
883 continue;
884 }
885
886 for (int i(0); i<myModuleIDs->GetNrows(); ++i){
887 //Coverting back from a double to a unvi signed long long
888 double source = (*myModuleIDs)(i);
890 memcpy(&target, &source, sizeof(target));
891 newModIndexMap[i]=target;
892 //std::cout << i<< " " <<target <<"\n";
893 }
894
895 delete myModuleIDs;
896
897 std::map<int,unsigned long long> newDoFMap;
898
899 TVectorD* myDoFs;
900 myDoFs = (TVectorD*)myFile->Get("dof");
901 if( !myDoFs ){
902 ++numberOfReadErrors;
903 ATH_MSG_ERROR("DoFs not read!!!");
904 continue;
905 }
906
907 for (int i(0); i<myDoFs->GetNrows(); ++i){
908 //Coverting back from a double to a unsigned long long
909 double source = (*myDoFs)(i);
911 memcpy(&target, &source, sizeof(target));
912 newDoFMap[i]=target;
913 }
914 delete myDoFs;
915
916
917 TVectorD* Scale;
918 Scale = (TVectorD*)myFile->Get("Scale");
919 if( !Scale ){
920 ++numberOfReadErrors;
921 ATH_MSG_ERROR("Scale not read!!!");
922 continue;
923 }
924
925 double scale=(*Scale)(0);
926 totalscale += scale;
927 delete Scale;
928
929
930 ATH_MSG_DEBUG("Reading Vector");
931 TVectorD* vector = (TVectorD*)myFile->Get("Vector");
932 if( !vector ){
933 ++numberOfReadErrors;
934 ATH_MSG_ERROR("Vector not read!!!");
935 continue;
936 }
937
938 AlVec* newVector = new AlVec(nDoF);
939 newVector->SetPathBin(m_pathbin.value()+m_prefixName.value());
940 newVector->SetPathTxt(m_pathtxt.value()+m_prefixName.value());
941
942 if (newVector->size() != m_bigvector->size() ) {
943 msg(MSG::FATAL) << "vector wrong size! newVector size " << newVector->size()
944 << ", bigvector size " << m_bigvector->size()<<endmsg;
945 delete newVector;
946 delete vector;
947 return false;
948 }
949
950 if (m_bigvector->size() != vector->GetNrows() ) {
951 msg(MSG::FATAL) << "File vector wrong size! File Vector size " << vector->GetNrows()
952 << ", bigvector size " << m_bigvector->size()<<endmsg;
953 delete newVector;
954 delete vector;
955 return false;
956 }
957
958
959 for (int i=0;i<nDoF;i++) {
960 (*newVector)[i] = (*vector)(i);
961 }
962 delete vector;
963
964 // check modIndexMaps to make sure they are the same
965 if (ifile == 0){
966 DoFMap = std::move(newDoFMap);
967 } else if (DoFMap!=newDoFMap) {
968 delete newVector;
969 delete vector;
970 msg(MSG::FATAL) << "module dofs don't agree!" << endmsg;
971 return false;
972 }
973
974 if (ifile == 0){
975 modIndexMap = std::move(newModIndexMap);
976 } else if (modIndexMap!=newModIndexMap) {
977 delete newVector;
978 delete vector;
979 msg(MSG::FATAL) << "module index maps don't agree!" << endmsg;
980 return false;
981 }
982
983 if (ifile>0){
984 *m_bigvector += *newVector;
985 delete newVector;
986 } else {
987 delete m_bigvector;
988 m_bigvector = newVector;
989 }
990
991
992 ATH_MSG_DEBUG("Reading matrix ");
993 TMatrixDSparse* matrix = (TMatrixDSparse*)myFile->Get("Matrix");
994
995 if( !matrix ){
996 ++numberOfReadErrors;
997 ATH_MSG_ERROR("Matrix not read!!!");
998 continue;
999 }
1000
1001
1002 if (ifile == 0 ){
1003
1004 accumMatrix = new AlSpaMat(nDoF);
1005 ATH_MSG_DEBUG("Matrix size b4 "<< accumMatrix->ptrMap()->size() );
1006
1007 //This method is ok for large matrix files... really only access the non zero elements
1008 for (int ii=0;ii<nDoF;ii++) {
1009 const TMatrixTSparseRow_const<double> myRow = (*matrix)[ii];
1010 int i = myRow.GetRowIndex();
1011 for (int jj=0;jj<myRow.GetNindex();jj++) {
1012 int j = (myRow.GetColPtr())[jj];
1013 const double myElement= (myRow.GetDataPtr())[jj];
1014 if (i<j){
1015 ATH_MSG_DEBUG("i < j " );
1016 j = i;
1017 i = (myRow.GetColPtr())[jj];
1018 }
1019 (*accumMatrix)[i][j] = myElement;
1020 }
1021 }
1022 ATH_MSG_DEBUG("Matrix size AF "<< accumMatrix->ptrMap()->size() );
1023
1024 } else if ( accumMatrix) {
1025 ATH_MSG_DEBUG("Matrix size b4 "<< accumMatrix->ptrMap()->size() );
1026
1027 for (int ii=0;ii<nDoF;ii++) {
1028 const TMatrixTSparseRow_const<double> myRow = (*matrix)[ii];
1029 int i = myRow.GetRowIndex();
1030 for (int jj=0;jj<myRow.GetNindex();jj++) {
1031 int j = (myRow.GetColPtr())[jj];
1032 const double myElement= (myRow.GetDataPtr())[jj];
1033 if (i<j){
1034 ATH_MSG_DEBUG("i < j " );
1035 j = i;
1036 i = (myRow.GetColPtr())[jj];
1037 }
1038 (*accumMatrix)[i][j] += myElement;
1039 }
1040 }
1041 ATH_MSG_DEBUG("Matrix size AF "<< accumMatrix->ptrMap()->size() );
1042
1043 } else {
1044 delete matrix;
1045 ++numberOfReadErrors;
1046 ATH_MSG_ERROR("Matrix allocation error!!!");
1047 continue;
1048 }
1049
1050 delete matrix;
1051
1052 TVectorD* hits;
1053 TVectorD* tracks;
1054
1055 ATH_MSG_DEBUG("Reading hitmap ");
1056 hits = (TVectorD*)myFile->Get("Hits");
1057 if( !hits ){
1058 ++numberOfReadErrors;
1059 ATH_MSG_ERROR("Hitmap 1 not read!!!");
1060 continue;
1061 }
1062
1063 tracks = (TVectorD*)myFile->Get("Tracks");
1064 if( !tracks ){
1065 delete hits;
1066 ++numberOfReadErrors;
1067 ATH_MSG_ERROR("Hitmap 2 not read!!!");
1068 continue;
1069 }
1070
1071 if(hits->GetNrows() != TotalHits.GetNrows() ){
1072 delete hits;
1073 delete tracks;
1074 ++numberOfReadErrors;
1075 ATH_MSG_ERROR("Hitmap size incorrect!!!");
1076 continue;
1077 }
1078
1079 TotalHits += (*hits);
1080 TotalTracks += (*tracks);
1081
1082 delete hits;
1083 delete tracks;
1084
1085 myFile->Close("R");
1086 delete myFile;
1087 //
1088 //setting 'itworked' here would be overwritten before use
1089 //
1090 ATH_MSG_DEBUG("Memory usage [MB], total " << myusage.ru_maxrss/1024 << ", increase " << (myusage.ru_maxrss-intialMemUse)/1024);
1091
1092 }
1093
1094
1095
1096 // create new matrix to read data from current file
1097 if(accumMatrix){
1098 if (symBigMatrix) {
1099 AlSymMat newMatrix(nDoF);
1100 //This method is ok for small matrix files
1101 for (int i=0;i<nDoF;i++) {
1102 for (int j=0;j<=i;j++) {
1103 newMatrix[i][j] = (*accumMatrix)[i][j];
1104 }
1105 }
1106
1107 *symBigMatrix += newMatrix;
1108 delete accumMatrix;
1109 } else if (spaBigMatrix) {
1110 ATH_MSG_DEBUG( "should reassign matrix "<< spaBigMatrix->ptrMap()->size() );
1111 *spaBigMatrix += *accumMatrix;
1112 ATH_MSG_DEBUG( "?????? "<< spaBigMatrix->ptrMap()->size() );
1113 delete accumMatrix;
1114 }
1115 }
1116
1117 ATH_MSG_DEBUG( "?????? "<< m_bigmatrix->ptrMap()->size() );
1118
1119 AlignModuleList::const_iterator imod = moduleList->begin();
1120 AlignModuleList::const_iterator imod_end = moduleList->end();
1121 int index = 0;
1122 int totalhits = 0;
1123 for(; imod != imod_end; ++imod, ++index ) {
1124 AlignModule * module = *imod;
1125 module->setNHits((int)TotalHits(index));
1126 module->setNTracks((int)TotalTracks(index));
1127 totalhits += (int)TotalHits(index);
1128 }
1129
1130
1131 m_nHits = totalhits;
1132 m_nTracks = 0;
1133 m_nMeasurements = 0;
1134 m_scale = totalscale;
1135
1136 return true;
1137 }
class TMatrixTSparse< double > TMatrixDSparse
Definition AlSpaMat.h:14
#define ATH_MSG_ERROR(x)
Gaudi::Property< int > m_maxReadErrors
Definition MatrixTool.h:251
Gaudi::Property< std::vector< std::string > > m_inputTFiles
Definition MatrixTool.h:245
void Scale(TH1 *h, double d=1)
str index
Definition DeMoScan.py:362
float j(const xAOD::IParticle &, const xAOD::TrackMeasurementValidation &hit, const Eigen::Matrix3d &jab_inv)
std::vector< AlignModule * > AlignModuleList

◆ addFirstDerivative()

void MatrixTool::addFirstDerivative ( int irow,
double firstderiv )
virtual

Implements Trk::IMatrixTool.

Definition at line 1484 of file MatrixTool.cxx.

1485 {
1486 (*m_bigvector)[irow] += firstderiv;
1487 }

◆ addFirstDerivatives() [1/2]

void MatrixTool::addFirstDerivatives ( AlVec * vector)
virtual

adds first derivative to vector

Implements Trk::IMatrixTool.

Definition at line 1464 of file MatrixTool.cxx.

1465 {
1466 }

◆ addFirstDerivatives() [2/2]

void MatrixTool::addFirstDerivatives ( std::list< int, double > & derivatives)
virtual

adds first derivative to vector for only some entries

Implements Trk::IMatrixTool.

Definition at line 1474 of file MatrixTool.cxx.

1475 {
1476 }

◆ addModule()

void IMatrixTool::addModule ( int alignModuleIndex,
int nAlignParam )
inlineinherited

Definition at line 116 of file IMatrixTool.h.

116{ m_alignModuleMap[alignModuleIndex]=m_nentries; m_nentries += nAlignParam; }
std::map< int, int > m_alignModuleMap

◆ addSecondDerivative()

void MatrixTool::addSecondDerivative ( int irow,
int icol,
double secondderiv )
virtual

Implements Trk::IMatrixTool.

Definition at line 1490 of file MatrixTool.cxx.

1491 {
1492 (*m_bigmatrix)[irow][icol] += secondderiv;
1493 }

◆ addSecondDerivatives() [1/2]

void MatrixTool::addSecondDerivatives ( AlSymMatBase * matrix)
virtual

adds second derivatives to matrix

Implements Trk::IMatrixTool.

Definition at line 1469 of file MatrixTool.cxx.

1470 {
1471 }

◆ addSecondDerivatives() [2/2]

void MatrixTool::addSecondDerivatives ( std::list< std::pair< int, int >, double > & derivatives)
virtual

adds first derivative to vector for only some entries

Implements Trk::IMatrixTool.

Definition at line 1479 of file MatrixTool.cxx.

1480 {
1481 }

◆ allocateMatrix()

StatusCode MatrixTool::allocateMatrix ( int nDoF = 0)
virtual

allocates memory for big matrix and big vector

Implements Trk::IMatrixTool.

Definition at line 100 of file MatrixTool.cxx.

101 {
102 ATH_MSG_INFO("allocating matrix and vector with nDoF = "<<nDoF);
103
104 if (nullptr!=m_bigmatrix || nullptr!=m_bigvector)
105 ATH_MSG_ERROR("big matrix already allocated!");
106
107 // Decide upon the big matrix representation:
108 if( m_useSparse )
109 m_bigmatrix = new AlSpaMat(nDoF);
110 else
111 m_bigmatrix = new AlSymMat(nDoF);
112
113 m_bigvector = new AlVec(nDoF);
114
115 ATH_MSG_INFO(" After Matrix and Vector allocation");
116
117 // set paths for matrix and vector output
118 m_bigmatrix->SetPathBin(m_pathbin.value()+m_prefixName);
119 m_bigmatrix->SetPathTxt(m_pathtxt.value()+m_prefixName);
120 m_bigvector->SetPathBin(m_pathbin.value()+m_prefixName);
121 m_bigvector->SetPathTxt(m_pathtxt.value()+m_prefixName);
122
123 ATH_MSG_INFO("set path to "<<m_pathbin.value()+m_prefixName.value());
124 return StatusCode::SUCCESS;
125 }

◆ declareGaudiProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< AlgTool > >::declareGaudiProperty ( Gaudi::Property< T, V, H > & hndl,
const SG::VarHandleKeyType &  )
inlineprivateinherited

specialization for handling Gaudi::Property<SG::VarHandleKey>

Definition at line 156 of file AthCommonDataStore.h.

158 {
160 hndl.value(),
161 hndl.documentation());
162
163 }
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)

◆ declareProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( Gaudi::Property< T, V, H > & t)
inlineinherited

Definition at line 145 of file AthCommonDataStore.h.

145 {
146 typedef typename SG::HandleClassifier<T>::type htype;
148 }
Gaudi::Details::PropertyBase & declareGaudiProperty(Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
specialization for handling Gaudi::Property<SG::VarHandleKey>

◆ detStore()

const ServiceHandle< StoreGateSvc > & AthCommonDataStore< AthCommonMsg< AlgTool > >::detStore ( ) const
inlineinherited

The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 95 of file AthCommonDataStore.h.

◆ entryNumber()

int IMatrixTool::entryNumber ( int alignModuleIndex)
inlineinherited

Definition at line 118 of file IMatrixTool.h.

118{ if ( m_alignModuleMap.find(alignModuleIndex) == m_alignModuleMap.end()) return -1; else return m_alignModuleMap[alignModuleIndex]; }

◆ evtStore()

ServiceHandle< StoreGateSvc > & AthCommonDataStore< AthCommonMsg< AlgTool > >::evtStore ( )
inlineinherited

The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 85 of file AthCommonDataStore.h.

◆ extraDeps_update_handler()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::extraDeps_update_handler ( Gaudi::Details::PropertyBase & ExtraDeps)
protectedinherited

Add StoreName to extra input/output deps as needed.

use the logic of the VarHandleKey to parse the DataObjID keys supplied via the ExtraInputs and ExtraOuputs Properties to add the StoreName if it's not explicitly given

◆ fillVecMods()

int MatrixTool::fillVecMods ( )
staticprivate

Definition at line 1627 of file MatrixTool.cxx.

1628 {
1629 return 0;
1630 }

◆ finalize()

StatusCode MatrixTool::finalize ( )

initialize

Definition at line 91 of file MatrixTool.cxx.

92 {
93 ATH_MSG_DEBUG("finalize() of MatrixTool");
94
95 return StatusCode::SUCCESS;
96 }

◆ initialize()

StatusCode MatrixTool::initialize ( )

initialize

Definition at line 70 of file MatrixTool.cxx.

71 {
72 ATH_MSG_DEBUG("initialize() of MatrixTool");
73
74 // get AlignModuleTool
75 if (m_alignModuleTool.retrieve().isSuccess())
76 ATH_MSG_INFO("Retrieved " << m_alignModuleTool);
77 else{
78 msg(MSG::FATAL) << "Could not get " << m_alignModuleTool << endmsg;
79 return StatusCode::FAILURE;
80 }
81
82 ATH_MSG_INFO("Retrieving data from the following files: ");
83 for (auto & inputVectorFile : m_inputVectorFiles) {
84 ATH_MSG_INFO(m_pathbin+inputVectorFile);
85 }
86
87 return StatusCode::SUCCESS;
88 }

◆ inputHandles()

virtual std::vector< Gaudi::DataHandle * > AthCommonDataStore< AthCommonMsg< AlgTool > >::inputHandles ( ) const
overridevirtualinherited

Return this algorithm's input handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

◆ interfaceID()

const InterfaceID & IMatrixTool::interfaceID ( )
inlinestaticinherited

Retrieve interface ID.

Definition at line 110 of file IMatrixTool.h.

110 {
112 }
static const InterfaceID IID_TRKALIGNINTERFACES_IMatrixTool("IMatrixTool", 1, 0)

◆ msg()

MsgStream & AthCommonMsg< AlgTool >::msg ( ) const
inlineinherited

Definition at line 24 of file AthCommonMsg.h.

24 {
25 return this->msgStream();
26 }

◆ msgLvl()

bool AthCommonMsg< AlgTool >::msgLvl ( const MSG::Level lvl) const
inlineinherited

Definition at line 30 of file AthCommonMsg.h.

30 {
31 return this->msgLevel(lvl);
32 }

◆ nHits()

int Trk::IMatrixTool::nHits ( ) const
inlineinherited

Definition at line 84 of file IMatrixTool.h.

84{ return m_nHits; }

◆ nMeasurements()

int Trk::IMatrixTool::nMeasurements ( ) const
inlineinherited

Definition at line 92 of file IMatrixTool.h.

92{ return m_nMeasurements; }

◆ nTracks()

int Trk::IMatrixTool::nTracks ( ) const
inlineinherited

Definition at line 88 of file IMatrixTool.h.

88{ return m_nTracks; }

◆ outputHandles()

virtual std::vector< Gaudi::DataHandle * > AthCommonDataStore< AthCommonMsg< AlgTool > >::outputHandles ( ) const
overridevirtualinherited

Return this algorithm's output handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

◆ postSolvingLapack()

void MatrixTool::postSolvingLapack ( AlVec * dChi2,
AlSymMat * d2Chi2,
AlVec & w,
AlMat & z,
int size )
private

Definition at line 1826 of file MatrixTool.cxx.

1827 {
1828 ATH_MSG_DEBUG("in postSolvinglapack()");
1829
1830 if( z.ncol() != size) {
1831 msg(MSG::ERROR)<<"Eigenvector matrix has incorrect size : "<<z.ncol()<<" != "<<size<<endmsg;
1832 return;
1833 }
1834
1835 if( (int)m_activeIndices.size() != size) {
1836 msg(MSG::ERROR)<<"Number of active parameters is incorrect : "<<m_activeIndices.size()<<" != "<<size<<endmsg;
1837 return;
1838 }
1839
1840 // Compute bigvector in diagonal basis (Vb = Ut * bigvector)
1841 AlVec D(size);
1842 D = z*(*dChi2);
1843
1844 if (m_writeEigenMat) {
1845
1846 ATH_MSG_INFO("writing the eigenvectors in a matrix: "<< z.nrow() << "x" << z.ncol());
1847
1848 // Set Path for the z matrix (eigenvector matrix)
1849 z.SetPathBin(m_pathbin.value()+m_prefixName.value());
1850 z.SetPathTxt(m_pathtxt.value()+m_prefixName.value());
1851
1852 ATH_MSG_INFO("writing the eigenvector matrix: "<< m_scalaMatName);
1853 ATH_MSG_DEBUG("matrix will be in: "<< m_pathbin.value()+m_prefixName.value()+m_scalaMatName.value());
1854
1855 StatusCode sc = z.Write(m_scalaMatName, true); // write the eigenvector matrix
1856
1857 if (sc!=StatusCode::SUCCESS)
1858 msg(MSG::ERROR)<<"Problem writing eigenvector matrix"<<endmsg;
1859
1860 // Set Path for the w matrix (eigenvalues matrix - diagonal bigmatrix)
1861 w.SetPathBin(m_pathbin.value()+m_prefixName.value());
1862 w.SetPathTxt(m_pathtxt.value()+m_prefixName.value());
1863
1864 ATH_MSG_INFO("writing the eigenvectors in a vector: "<< w.size());
1865 ATH_MSG_INFO("writing the eigenvalues vector (diagonal bigmatrix): "<< m_scalaVecName);
1866 ATH_MSG_DEBUG("vector will be in: "<< m_pathbin.value()+m_prefixName.value()+m_scalaVecName.value());
1867
1868 sc = w.WriteEigenvalueVec(m_scalaVecName, true); // write the eigenvalues vecor
1869
1870 if (sc!=StatusCode::SUCCESS)
1871 msg(MSG::ERROR)<<"Problem writing eigenvector matrix"<<endmsg;
1872
1873 if (m_writeEigenMatTxt) {
1874 sc = z.Write("eigenvectors.txt", false);
1875 if (sc!=StatusCode::SUCCESS)
1876 msg(MSG::ERROR)<<"Problem writing eigenvector matrix to text file"<<endmsg;
1877 sc = w.WriteEigenvalueVec("eigenvalues.txt", false);
1878 if (sc!=StatusCode::SUCCESS)
1879 msg(MSG::ERROR)<<"Problem writing eigenvalue vector to text file"<<endmsg;
1880 }
1881
1882 }
1883
1884 // Set eigenvalue thresholds
1885 const double eigenvalue_threshold = 1e-19;
1886
1887 // weak mode removal
1888 if (m_modcut == -1) {
1889
1890 ATH_MSG_INFO(" Starting the automatic Weak Mode Removal method");
1891
1892 // create a pull vector for the alignment corrections in diagonal basis (db_pulls)
1893 //int nDoF=m_alignModuleTool->nAlignParameters();
1894 AlVec* Align_db = new AlVec(size);
1895 AlVec* Align_error_db = new AlVec(size);
1896 AlVec* AlignPull = new AlVec(size);
1897 ATH_MSG_DEBUG("AlignPull vector size is: "<< (*AlignPull).size());
1898
1899 m_modcut = 0;
1900 bool wm_stop = false;
1901
1902 // -----------------------------------------------------------------------
1903 // First pass: removing weak modes because of large pull
1904 // compute alignment pulls for corrections in diagonal basis (db)
1905 for(int i=0; i<size; i++) {
1906
1907 (*Align_db)[i] = (-D[i]/w[i]);
1908 if(m_scale<=0.)[[unlikely]]{
1909 ATH_MSG_WARNING("postSolvingLapack: Scaling requested but scale not set. Not scaling matrix and vector.");
1910 } else {
1911 (*Align_error_db)[i] = sqrt(1.0/w[i]/m_scale);
1912 }
1913
1914 if (w[i]<eigenvalue_threshold) {
1915 ATH_MSG_INFO(" + EigenMode " << i
1916 << " removed as eigenvalue lower than the threshold " << eigenvalue_threshold
1917 << ": " << w[i]);
1918 (*AlignPull)[i] = 0.0;
1919 ++m_modcut;
1920 }
1921 else
1922 (*AlignPull)[i] = (*Align_db)[i] / (*Align_error_db)[i];
1923
1924 ATH_MSG_DEBUG(i << ". AlignPar: " << (*Align_db)[i] << " +- " << (*Align_error_db)[i]
1925 << " (pull: " << (*AlignPull)[i] << ") ; w[i]: " << w[i]);
1926 }
1927 ATH_MSG_INFO(" +++ Weak Mode removal ++ stop after mode "<< m_modcut << " (First pass)");
1928 // -----------------------------------------------------------------------
1929
1930 // -----------------------------------------------------------------------
1931 // Second pass
1932 // if the error is greater than the correction -> cut this mode
1933 for(int i=m_modcut; (i<size && !wm_stop); i++) {
1934
1935 // if the error is greater than the correction -> cut this mode
1936 if (fabs((*AlignPull)[i])<m_pullcut) {
1937 ATH_MSG_INFO(" + EigenMode " << i
1938 << " removed as pull is lower than " << m_pullcut << ": "
1939 << (*AlignPull)[i]);
1940 ++m_modcut;
1941 }
1942 else
1943 wm_stop = true;
1944
1945 }
1946 ATH_MSG_INFO(" +++ Weak Mode removal ++ stop after mode "<< m_modcut << " (Second pass)");
1947 // -----------------------------------------------------------------------
1948
1949 wm_stop = false;
1950
1951 // ----------------------------------------------------------------------
1952 // Third pass
1953 // Check if the next eigenvalues is far away. If it is two orders of
1954 // magnitude bigger remove also this mode and allow the search
1955 for(int i=m_modcut; (i<size && !wm_stop); i++) {
1956
1957 // if the next eigenvalues is far away -> cut this mode
1958 if (m_eigenvalueStep*w[i]<w[i+1]) {
1959 ATH_MSG_INFO(" + EigenMode " << i
1960 << " removed as diff between eigenvalues, " << w[i] << " and " << w[i+1]
1961 << ", is greater than " << m_eigenvalueStep);
1962 ++m_modcut;
1963 }
1964 else
1965 wm_stop = true;
1966
1967 }
1968 ATH_MSG_INFO(" +++ Weak Mode removal ++ stop after mode "<< m_modcut << " (Third pass)");
1969 // -----------------------------------------------------------------------
1970
1971 wm_stop = false;
1972
1973 // -----------------------------------------------------------------------
1974 // Fourth pass
1975 // Check if the next eigenvalues is far away. If it is two orders of
1976 // magnitude bigger remove also this mode and allow the search
1977 for(int i=m_modcut; (i<size && !wm_stop); i++) {
1978
1979 // if the next eigenvalues is far away -> cut this mode
1980 if ( fabs((*Align_db)[i]) > m_Align_db_step*fabs((*Align_db)[i+1]) ) {
1981 ATH_MSG_INFO(" + EigenMode " << i
1982 << " removed as diff between corrections, " << w[i] << " and " << w[i+1]
1983 << ", is greater than "
1984 << m_Align_db_step);
1985 ++m_modcut;
1986 }
1987 else
1988 wm_stop = true;
1989
1990 }
1991 ATH_MSG_INFO(" +++ Weak Mode removal ++ stop after mode "<< m_modcut << " (Fourth pass)");
1992 // -----------------------------------------------------------------------------------------------------
1993
1994 // Free memory and clear the pointer to
1995 // prevent using invalid memory reference
1996 delete Align_db;
1997 delete Align_error_db;
1998 delete AlignPull;
1999 Align_db = nullptr;
2000 Align_error_db = nullptr;
2001 AlignPull = nullptr;
2002
2003 } // end of if(m_modcut == -1)
2004
2005 // Save some stuff to debug purposes
2006 /*
2007 if (m_storeDia) {
2008 std::string path = m_pathtxt+m_prefixName+"align_dia";
2009 std::fstream orthogon(path.c_str(), std::ios::out);
2010 orthogon.setf(std::ios::fixed);
2011 orthogon.setf(std::ios::showpoint);
2012 orthogon.precision(6);
2013
2014 orthogon << std::setw(10)
2015 << "--------------------------------------------------------------------------------"
2016 << std::endl;
2017 orthogon << std::setw(10) << " ModeCut = " << m_modcut << std::endl;
2018 orthogon << std::setw(10)
2019 << "--------------------------------------------------------------------------------"
2020 << std::endl;
2021
2022 orthogon << std::setw(10) << "mode"
2023 << std::setw(20) << "eigenmode"
2024 << std::setw(20) << "eigenmode error"
2025 << std::setw(25) << "eigenvalue"
2026 << std::endl;
2027
2028 for( int m=0; m<size; m++) {
2029
2030 // mode
2031 orthogon << std::setw(10) << m;
2032
2033 // eigenmode (db)
2034 if( w[m]>1.0e-15) orthogon << std::setw(20) << -D[m]/w[m];
2035 else orthogon << std::setw(20) << 0.0;
2036
2037 // error eigenmode (error_db)
2038 if( w[m]>1.0e-15) orthogon << std::setw(20) << sqrt(1.0/w[m]/m_scale);
2039 else orthogon << std::setw(20) << 0.0;
2040
2041 // eigenvalues
2042 orthogon << std::setw(25) << w[m] << std::endl;
2043 }
2044 orthogon.close();
2045 } // end store align_dia.txt
2046 */
2047
2048 AlVec delta(size);
2049 AlVec deltafull(size);
2050 AlVec errSq(size);
2051
2052 // full covariance matrix
2053 CLHEP::HepSymMatrix * cov = nullptr;
2055 // Warning ! The matrix can be huge!
2056 // This can lead to memory problems
2057 cov = new CLHEP::HepSymMatrix(size,0);
2058
2059 if(m_logStream)
2060 *m_logStream<<"/------ The Eigenvalue Spectrum -------\n";
2061
2062 for (int i=0;i<size;i++) {
2063 AlVec thisdelta(size);
2064 for(int j=0;j<size;j++)
2065 thisdelta[j] = z[i][j] * (-D[i]/w[i]);
2066 deltafull += thisdelta;
2067
2068 ATH_MSG_DEBUG("eigenvalue "<<w[i]);
2069 if( i<m_modcut ) {
2070 ATH_MSG_INFO("skipping eigenvalue "<<w[i]<<" , modcut is "<<m_modcut);
2071 if(m_logStream)
2072 *m_logStream<<"| skipping eigenvalue "<<w[i]<<"\n";
2073 }
2074 else if( w[i] < m_eigenvaluethreshold ) {
2075 ATH_MSG_INFO("skipping eigenvalue "<<w[i]<<" , cut is "<<m_eigenvaluethreshold);
2076 if(m_logStream)
2077 *m_logStream<<"| skipping eigenvalue "<<w[i]<<"\n";
2078 }
2079 else {
2080 if(m_logStream)
2081 *m_logStream<<"| "<<w[i]<<"\n";
2082
2083 delta += thisdelta;
2084 for(int j=0;j<size;j++) {
2085 errSq[j] += z[i][j] * z[i][j] / w[i];
2087 for(int k=0;k<=j;k++)
2088 (*cov)[j][k] += z[i][j] * z[i][k] / w[i];
2089 }
2090 }
2091 }
2092 }
2093
2094 if(m_logStream)
2095 *m_logStream<<"\\----- End of Eigenvalue Spectrum -----\n";
2096
2097 ATH_MSG_DEBUG("Alignment constants:");
2098
2099 DataVector<AlignPar> * alignParList = m_alignModuleTool->alignParList1D();
2100
2101 // compute alignment corrections (translations in mm and rotations in rad) and their variances
2102 for(int i=0; i<size; i++) {
2103
2104 double param = delta[i];
2105 double err = sqrt(2.*std::fabs(errSq[i]));
2106
2107 int idof = m_activeIndices[i];
2108 AlignPar * alignPar=(*alignParList)[idof];
2109
2110 // undo the sigma scaling
2111 double sigma = alignPar->sigma();
2112
2113 param *= sigma;
2114 err *= sigma;
2115
2116 // undo normalization scaling of error
2117 if(m_scaleMatrix && m_scale>0.)
2118 err /= sqrt(m_scale);
2119
2120 ATH_MSG_DEBUG(i <<" : "<< param << " +/- "<< err);
2121 ATH_MSG_DEBUG("cov("<<i<<")="<<errSq[i]<<", sigma: "<<sigma);
2122 ATH_MSG_DEBUG("init par: "<<alignPar->initPar());
2123 alignPar->setPar(param, err);
2124 ATH_MSG_DEBUG("set param to "<<param<<" for alignPar "<<alignPar);
2125 ATH_MSG_DEBUG(*(*alignParList)[idof]);
2126 }
2127
2128 if(m_logStream) {
2130
2131 // norm of first derivative
2132 double norm1st = dChi2->norm();
2133 if(m_scaleMatrix && m_scale>0.) // undo normalization scaling
2134 norm1st *= m_scale;
2135 *m_logStream<<"norm of first derivative : "<<norm1st<<"\n";
2136
2137 if(d2Chi2) {
2138 // distance to solution
2139 double dist = ( (*d2Chi2) * deltafull + (*dChi2) ).norm();
2140 if(m_scaleMatrix && m_scale>0.) // undo normalization scaling
2141 dist *= m_scale;
2142 *m_logStream<<"distance to solution : "<<dist<<"\n";
2143
2144 // calculate chi2 of the alignment change
2145 double chi2 = delta * (*d2Chi2) * delta * .5;
2146 if(m_scaleMatrix && m_scale>0.) // undo normalization scaling
2147 chi2 *= m_scale;
2148 *m_logStream<<"delta(chi2) of the alignment change : "<<chi2<<" / "<<size<<"\n";
2149 }
2150 }
2151
2152 delete cov;
2153 }
size_t size() const
Number of registered mappings.
std::ostream * m_logStream
logfile output stream
Definition IMatrixTool.h:98
Gaudi::Property< bool > m_writeEigenMatTxt
Definition MatrixTool.h:184
std::vector< int > m_activeIndices
vector of indices which pass the min-hits cut
Definition MatrixTool.h:248
Gaudi::Property< double > m_calculateFullCovariance
Definition MatrixTool.h:215
Gaudi::Property< std::string > m_scalaVecName
Definition MatrixTool.h:231
Gaudi::Property< float > m_eigenvalueStep
Definition MatrixTool.h:169
Gaudi::Property< float > m_Align_db_step
Definition MatrixTool.h:171
Gaudi::Property< bool > m_scaleMatrix
Definition MatrixTool.h:205
Gaudi::Property< bool > m_writeEigenMat
Definition MatrixTool.h:182
Gaudi::Property< int > m_modcut
Definition MatrixTool.h:161
Gaudi::Property< double > m_eigenvaluethreshold
Definition MatrixTool.h:157
Gaudi::Property< std::string > m_scalaMatName
Definition MatrixTool.h:229
Gaudi::Property< float > m_pullcut
Definition MatrixTool.h:167
void printGlobalSolution(std::ostream &os, const CLHEP::HepSymMatrix *cov)
double chi2(TH1 *h0, TH1 *h1)
@ z
global position (cartesian)
Definition ParamDefs.h:57
#define unlikely(x)

◆ prepareBinaryFiles()

void MatrixTool::prepareBinaryFiles ( int solveOption)
virtual

reads/writes matrix entries from/to binary files as necessary

Implements Trk::IMatrixTool.

Definition at line 128 of file MatrixTool.cxx.

129 {
130 }

◆ printGlobalSolution() [1/2]

void MatrixTool::printGlobalSolution ( std::ostream & os,
const CLHEP::HepSymMatrix * cov )

Definition at line 1496 of file MatrixTool.cxx.

1497 {
1498 const AlignModuleList * alignModules = m_alignModuleTool->alignModules1D();
1499
1500 AlignModuleList::const_iterator imod = alignModules->begin();
1501 AlignModuleList::const_iterator imod_end = alignModules->end();
1502 for( ; imod!=imod_end; ++imod) {
1503 AlignModule * module = *imod;
1504
1505 DataVector<AlignPar> * alignPars = m_alignModuleTool->getAlignPars(module);
1506 int thisNDoF = alignPars->size();
1507
1508 // fill local covariance matrix
1509 CLHEP::HepSymMatrix * covsub = nullptr;
1510 if(cov && module->nHits() >= m_minNumHits && module->nTracks() >= m_minNumTrks) {
1511 covsub = new CLHEP::HepSymMatrix(thisNDoF,0);
1512 for (int i=0;i<thisNDoF;++i) {
1513 int ipar = alignPars->at(i)->index();
1514 double sigma_i = alignPars->at(i)->sigma();
1515
1516 std::vector<int>::iterator itActive = std::find(m_activeIndices.begin(),m_activeIndices.end(),ipar);
1517 if( itActive == m_activeIndices.end() )
1518 continue;
1519 int iActive = std::distance(m_activeIndices.begin(),itActive);
1520
1521 for (int j=0;j<=i;++j) {
1522 int jpar = alignPars->at(j)->index();
1523 double sigma_j = alignPars->at(j)->sigma();
1524
1525 std::vector<int>::iterator jtActive = std::find(m_activeIndices.begin(),m_activeIndices.end(),jpar);
1526 if( jtActive == m_activeIndices.end() )
1527 continue;
1528 int jActive = std::distance(m_activeIndices.begin(),jtActive);
1529
1530 (*covsub)[i][j] = (*cov)[iActive][jActive] * sigma_i * sigma_j;
1531 }
1532 }
1533 }
1534
1535 printModuleSolution(os,module,covsub);
1536
1537 delete covsub;
1538 }
1539 os << "--------------------------------------------------------------------------------" << std::endl;
1540 }
const T * at(size_type n) const
Access an element, as an rvalue.
Gaudi::Property< int > m_minNumHits
Definition MatrixTool.h:163
void printModuleSolution(std::ostream &os, const AlignModule *module, const CLHEP::HepSymMatrix *cov) const
namespace { class RestoreIOSFlags { public: RestoreIOSFlags (std::ostream &os) : m_os(&os),...
Gaudi::Property< int > m_minNumTrks
Definition MatrixTool.h:165

◆ printGlobalSolution() [2/2]

void MatrixTool::printGlobalSolution ( std::ostream & os,
const TMatrixDSym * cov )

Definition at line 1543 of file MatrixTool.cxx.

1544 {
1545 CLHEP::HepSymMatrix * cov = nullptr;
1546 if(cov0) {
1547 int nsize = cov0->GetNrows();
1548 cov = new CLHEP::HepSymMatrix(nsize,0);
1549
1550 for(int i=0; i<nsize; i++)
1551 for(int j=0; j<=i; j++)
1552 (*cov)[i][j] = (*cov0)[i][j];
1553 }
1554
1555 printGlobalSolution(os,cov);
1556
1557 delete cov;
1558 }

◆ printModuleSolution()

void MatrixTool::printModuleSolution ( std::ostream & os,
const AlignModule * module,
const CLHEP::HepSymMatrix * cov ) const

namespace { class RestoreIOSFlags { public: RestoreIOSFlags (std::ostream &os) : m_os(&os), m_flags(m_os->flags()), m_precision(m_os->precision()) {} ~RestoreIOSFlags() { m_os->flags(m_flags); m_os->precision(m_precision); } private: std::ostream *m_os; std::ios_base::fmtflags m_flags; std::streamsize m_precision; }; }

Definition at line 1581 of file MatrixTool.cxx.

1582 {
1583 os << "--------------------------------------------------------------------------------" << std::endl;
1584 os << "Alignment parameters for module: " << module->name() << std::endl;
1585 os << "Number of tracks passing: " << module->nTracks() << std::endl;
1586 if(m_minNumHits>0 && module->nHits()<m_minNumHits) {
1587 os << "Number of hits too small: "<<module->nHits()<<" < "<<m_minNumHits<<" Skipping the module\n";
1588 return;
1589 }
1590 if(m_minNumTrks>0 && module->nTracks()<m_minNumTrks) {
1591 os << "Number of tracks too small: "<<module->nTracks()<<" < "<<m_minNumTrks<<" Skipping the module\n";
1592 return;
1593 }
1594 os << "Number of hits seen: " << module->nHits() << std::endl;
1595 os << "Number of tracks seen: " << module->nTracks() << std::endl;
1596
1597 DataVector<AlignPar> * alignPars = m_alignModuleTool->getAlignPars(module);
1598 int thisNDoF = alignPars->size();
1599
1600 if(alignPars->empty())
1601 os << "No active parameters" << std::endl;
1602 else
1603 {
1604 // output alignment parameters and errors
1605 DataVector<AlignPar>::const_iterator ipar = alignPars->begin();
1606 DataVector<AlignPar>::const_iterator ipar_end = alignPars->end();
1607 for ( ; ipar != ipar_end; ++ipar) {
1608 const AlignPar * par = *ipar;
1609 os << std::format("{:<10}{:<12.5g} +/- {:<12.5g}\n",
1610 par->dumpType(),
1611 par->par(),
1612 par->err());
1613 }
1614
1615 if(cov) {
1616 // calculate local correlation matrix
1617 CLHEP::HepSymMatrix corrsub(thisNDoF,0);
1618 for(int irow=0; irow<thisNDoF; ++irow)
1619 for(int icol=0; icol<=irow; ++icol)
1620 corrsub[irow][icol] = (*cov)[irow][icol] / sqrt((*cov)[irow][irow] * (*cov)[icol][icol]);
1621 os << "Local correlation matrix: " << corrsub << std::flush;
1622 }
1623 }
1624 }
DataModel_detail::const_iterator< DataVector > const_iterator
Standard const_iterator.
Definition DataVector.h:838
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
bool empty() const noexcept
Returns true if the collection is empty.

◆ readHitmaps()

void MatrixTool::readHitmaps ( )
private

Definition at line 2327 of file MatrixTool.cxx.

2328 {
2329 ATH_MSG_INFO("read hitmaps from files");
2330
2331 const AlignModuleList * moduleList = m_alignModuleTool->alignModules1D();
2332 int nModules = moduleList->size();
2333
2334 AlMat hitmap(nModules,2);
2335 int nFiles = (int)m_inputHitmapFiles.size();
2336 for(int imap=0;imap<nFiles; imap++) {
2337 AlMat nextmap(nModules,2);
2338
2339 ATH_MSG_INFO("Reading hitmap "<<imap<<" from file "<<m_inputHitmapFiles[imap]);
2340
2341 if(nextmap.ReadScalaPack(m_inputHitmapFiles[imap]).isFailure()) {
2342 ATH_MSG_WARNING("Problem reading hitmap from \'"<<m_inputHitmapFiles[imap]<<"\'. Skipping.");
2343 continue;
2344 }
2345
2346 if(nextmap.nrow()!=nModules || nextmap.ncol()!=2) {
2347 ATH_MSG_WARNING("Matrix in file \'"<<m_inputHitmapFiles[imap]<<"\' has wrong size ("
2348 <<nextmap.nrow()<<" x "<<nextmap.ncol()<<"), should be ("<<nModules<<" x 2). Skipping.");
2349 continue;
2350 }
2351
2352 hitmap += nextmap;
2353 }
2354
2355 AlignModuleList::const_iterator imod = moduleList->begin();
2356 AlignModuleList::const_iterator imod_end = moduleList->end();
2357 int index = 0;
2358 int totalhits = 0;
2359 for(; imod != imod_end; ++imod) {
2360 AlignModule * module = *imod;
2361 //if ((int)hitmap[index][0]!=(int)module->identify().get_identifier32().get_compact()) ATH_MSG_ERROR("bad module identifier");
2362 //module->setIdentifier((Identifier)hitmap[index][0]);
2363 module->setNHits((int)hitmap[index][0]);
2364 module->setNTracks((int)hitmap[index][1]);
2365 totalhits += (int)hitmap[index][0];
2366 index++;
2367 }
2368
2369 m_nHits = totalhits;
2370 m_nTracks = 0;
2371 m_nMeasurements = 0;
2372
2373 ATH_MSG_INFO("Hitmap accumulated from "<<nFiles<<" files with total of "<<totalhits<<" hits.");
2374 }
Gaudi::Property< std::vector< std::string > > m_inputHitmapFiles
Definition MatrixTool.h:241

◆ renounce()

std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > AthCommonDataStore< AthCommonMsg< AlgTool > >::renounce ( T & h)
inlineprotectedinherited

Definition at line 380 of file AthCommonDataStore.h.

381 {
382 h.renounce();
384 }
std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > renounce(T &h)

◆ renounceArray()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::renounceArray ( SG::VarHandleKeyArray & handlesArray)
inlineprotectedinherited

remove all handles from I/O resolution

Definition at line 364 of file AthCommonDataStore.h.

364 {
366 }

◆ setLogStream()

virtual void Trk::IMatrixTool::setLogStream ( std::ostream * os)
inlinevirtualinherited

sets the output stream for the logfile

Definition at line 95 of file IMatrixTool.h.

95{ m_logStream = os; }

◆ setNHits()

void Trk::IMatrixTool::setNHits ( int n)
inlineinherited

set module identifier

set number of hits

Definition at line 83 of file IMatrixTool.h.

◆ setNMeasurements()

void Trk::IMatrixTool::setNMeasurements ( int n)
inlineinherited

set number of measurements

Definition at line 91 of file IMatrixTool.h.

91{ m_nMeasurements = n; }

◆ setNTracks()

void Trk::IMatrixTool::setNTracks ( int n)
inlineinherited

set number of tracks

Definition at line 87 of file IMatrixTool.h.

87{ m_nTracks = n; }

◆ solve()

int MatrixTool::solve ( )
virtual

solves for alignment parameters

Implements Trk::IMatrixTool.

Definition at line 1143 of file MatrixTool.cxx.

1144 {
1145 // ============
1146 // solve
1147 // ============
1148 ATH_MSG_DEBUG("in MatrixTool::solve()");
1149
1150 // set normalization scale to number of hits for now
1151 if(m_scale<0)
1152 m_scale = m_nHits;
1153
1154 //-------------------------------------------------------
1155 // write matrix and vector to file
1156 if (m_writeMat) {
1157 // version has to be 2 to for reading matrices and vectors back in to work properly
1158 double dummyVersion(2.);
1159
1160 // make map of matrix entry to module index (set by geometry manager tool)
1161 std::map<int,unsigned long long> modIndexMap;
1162 std::map<int,std::string> modNameMap;
1163 DataVector<AlignPar>* alignPars = m_alignModuleTool->alignParList1D();
1164 for (int i=0;i<(int)alignPars->size();i++) {
1165 modIndexMap[i]=(*alignPars)[i]->alignModule()->identify().get_compact();
1166 modNameMap [i]=(*alignPars)[i]->alignModule()->name();
1167 }
1168
1169 // binary files
1170 ATH_MSG_DEBUG("writing binary files");
1171 StatusCode sc1 = m_bigmatrix->Write("matrix.bin",true,m_wSqMatrix,m_scale,dummyVersion);
1172 StatusCode sc2 = m_bigvector->WritePartial("vector.bin",true,m_scale,modIndexMap,dummyVersion);
1173 if (!sc1.isSuccess() || !sc2.isSuccess()) {
1174 msg(MSG::ERROR)<<"problem writing matrix or vector"<<endmsg;
1175 return -1;
1176 }
1177
1178 if (m_writeMatTxt) {
1179
1180 // text files
1181 ATH_MSG_DEBUG("writing text files");
1182 sc1 = m_bigmatrix->Write("matrix.txt",false,m_wSqMatrix,m_scale,dummyVersion);
1183 sc2 = m_writeModuleNames ?
1184 m_bigvector->WritePartial("vector.txt",false,m_scale,modNameMap,dummyVersion) :
1185 m_bigvector->WritePartial("vector.txt",false,m_scale,modIndexMap,dummyVersion);
1186
1187 if (!sc1.isSuccess() || !sc2.isSuccess()) {
1188 msg(MSG::ERROR)<<"problem writing matrix or vector"<<endmsg;
1189 return -1;
1190 }
1191 }
1192
1193 ATH_MSG_DEBUG("matrix and vector written to: "<<m_pathbin.value()+m_prefixName.value()<<"matrix.bin (.txt) and "<<m_pathbin.value()+m_prefixName.value()<<"vector.bin (.txt)");
1194 }
1195
1196 //-------------------------------------------------------
1197 // write hitmap to file
1198 if (m_writeHitmap)
1199 writeHitmap();
1200
1201 if(m_writeTFile)
1202 storeInTFile(m_pathbin.value()+m_prefixName.value()+m_tfileName.value());
1203
1204
1205 if(!m_runLocal && m_solveOption==0) {
1206 ATH_MSG_DEBUG("No solving requested.");
1207 return 1;
1208 }
1209
1210 //-------------------------------------------------------
1211 // rescale the vector and the matrix according to sigmas
1212 // and apply soft mode cut
1213
1214 ATH_MSG_DEBUG("rescaling the matrix/vector and applying the soft-mode-cut");
1215
1216 DataVector<AlignPar>* alignParList = m_alignModuleTool->alignParList1D();
1217 int nDoF = alignParList->size();
1218
1219
1220 const AlSymMat * chkMatrix = dynamic_cast<const AlSymMat*>(m_bigmatrix);
1221 if(chkMatrix){
1222 // Method when using the dense matrix
1223 for (int i=0;i<nDoF;i++) {
1224 // scale the vector
1225 double sigma_i = (*alignParList)[i]->sigma();
1226 double softCut = 2 * pow( (*alignParList)[i]->softCut() , -2 );
1227 (*m_bigvector)[i] *= sigma_i;
1228
1229 for (int j=0;j<=i;j++) {
1230 // scale the matrix
1231 if ((*chkMatrix)[i][j] != 0.) {
1232 double sigma_j = (*alignParList)[j]->sigma();
1233 (*m_bigmatrix)[i][j] *= sigma_i * sigma_j;
1234 }
1235 // apply soft-mode-cut
1236 if (i==j && m_softEigenmodeCut>0.){
1237 (*m_bigmatrix)[i][j] += m_softEigenmodeCut * softCut;
1238
1239 }
1240
1241 // set first and second derivatives on AlignPar
1242 if (i==j) {
1243 (*alignParList)[i]->setFirstDeriv((*m_bigvector)[i]/sigma_i);
1244 (*alignParList)[i]->setSecndDeriv((*m_bigmatrix)[i][i]/sigma_i/sigma_i);
1245 }
1246 }
1247 }
1248 } else {
1249 // Method when using the sparse matrix
1250 for (const datamap::value_type& p : *m_bigmatrix->ptrMap()) {
1251 int i = p.first.first;
1252 int j = p.first.second;
1253
1254 // Scale matrix
1255 double sigma_i = (*alignParList)[i]->sigma();
1256 double sigma_j = (*alignParList)[j]->sigma();
1257
1258 (*m_bigmatrix)[i][j] *= sigma_i * sigma_j;
1259
1260 }
1261
1262
1263 for (int i=0;i<nDoF;i++) {
1264 // scale the vector
1265 double sigma_i = (*alignParList)[i]->sigma();
1266 (*m_bigvector)[i] *= sigma_i;
1267 if (m_softEigenmodeCut >0. ){
1268 double softCut = 2 * pow( (*alignParList)[i]->softCut() , -2 );
1269 (*m_bigmatrix)[i][i] += m_softEigenmodeCut * softCut;
1270 ATH_MSG_DEBUG( "DOF "<< i <<" Nhits "<< (*alignParList)[i]->alignModule()->nHits() << " soft-mode-cut "<< (*alignParList)[i]->softCut() <<" -> " << m_softEigenmodeCut * softCut << " sigma_i "<< sigma_i << " Matrix: " << (*m_bigmatrix)[i][i] << " Vector: " << (*m_bigvector)[i]);
1271 }
1272 (*alignParList)[i]->setFirstDeriv((*m_bigvector)[i]/sigma_i);
1273 (*alignParList)[i]->setSecndDeriv((*m_bigmatrix)[i][i]/sigma_i/sigma_i);
1274 }
1275 }
1276
1277 unsigned long long OldPixelIdentifier = 37769216; //Identifier for the Pixel Detector
1278 unsigned long long IBLIdentifier = 33574912; //Identifier for the Pixel Detector
1279
1280 unsigned long long SCT_ECA_8_Identifier = 218116096; //Identifier for the SCT ECA Last Disk
1281 std::string SCT_ECA_8_Name = "SCT/EndcapA/Disk_8";
1282
1283
1284
1285 ATH_MSG_INFO("rescaling done");
1286 ATH_MSG_INFO("Javi: Printing (*alignParList)[i]->alignModule()->identify32()");
1287
1288 // select modules with non-zero tracks
1289 for(int i=0;i<nDoF;i++)
1290 {
1291 ATH_MSG_DEBUG(i);
1292 ATH_MSG_DEBUG((*alignParList)[i]->alignModule()->identify32());
1293 ATH_MSG_DEBUG((*alignParList)[i]->alignModule());
1294 ATH_MSG_DEBUG((*alignParList)[i]->alignModule()->name());
1295 ATH_MSG_DEBUG((*alignParList)[i]->paramType());
1296
1297 //Skip solving for Pixel or IBL:
1298 const auto & theParameterList = *alignParList;
1299 const auto & thisIdentifier = theParameterList[i]->alignModule()->identify32();
1300 const auto & thisName = theParameterList[i]->alignModule()->name();
1301 const auto & thisParameterType = theParameterList[i]->paramType();
1302 const bool oldPixel = (thisIdentifier == OldPixelIdentifier);
1303 const bool ibl = (thisIdentifier == IBLIdentifier);
1304 const bool SCTECA8 = (thisIdentifier == SCT_ECA_8_Identifier);
1305 const bool SCTECA8_n = (thisName.find(SCT_ECA_8_Name)!= std::string::npos);
1306
1308 if (SCTECA8)
1309 {ATH_MSG_INFO( "SCT ECA Last Disk DoF have been skipped in the solving because DeactivateSCT_ECA_LastDisk is set to True");
1310 continue;}
1311 if (SCTECA8_n)
1312 {ATH_MSG_INFO( "SCT ECA Last Disk DoF have been skipped in the solving because DeactivateSCT_ECA_LastDisk is set to True");
1313 continue;}
1314 }
1315
1316 if (m_AlignIBLbutNotPixel) //If m_AlignIBLbutNotPixel is set to True, Pixel will be skipped in the solving.
1317 if (oldPixel)
1318 {ATH_MSG_INFO( "Pixel DoF have been skipped in the solving because AlignIBLbutNotPixel is set to True");
1319 continue;}
1320
1321 if (m_AlignPixelbutNotIBL) //If m_AlignPixelbutNotIBL is set to True, IBL will be skipped in the solving.
1322 if (ibl)
1323 {ATH_MSG_INFO( "IBL DoF have been skipped in the solving because AlignPixelbutNotIBL is set to True");
1324 continue;}
1325
1326 //For specific DoF: (*alignParList)[i]->paramType() = 0,1,2,3,4,5,6 for Tx,Ty,Tz,Rx,Ry,Rz,Bx
1327 //Pixel Dofs:
1328 if (m_Remove_Pixel_Tx) //If m_Remove_Pixel_Tx is set to True, Pixel Tx will be skipped in the solving.
1329 if ( oldPixel and (thisParameterType == 0))
1330 {ATH_MSG_INFO( "Pixel Tx DoF has been skipped in the solving because Remove_Pixel_Tx is set to True");
1331 continue;}
1332
1333 if (m_Remove_Pixel_Ty) //If m_Remove_Pixel_Ty is set to True, Pixel Ty will be skipped in the solving.
1334 if ( oldPixel and (thisParameterType == 1))
1335 {ATH_MSG_INFO( "Pixel Ty DoF has been skipped in the solving because Remove_Pixel_Ty is set to True");
1336 continue;}
1337
1338 if (m_Remove_Pixel_Tz) //If m_Remove_Pixel_Tz is set to True, Pixel Tz will be skipped in the solving.
1339 if (oldPixel and (thisParameterType == 2))
1340 {ATH_MSG_INFO( "Pixel Tz DoF has been skipped in the solving because Remove_Pixel_Tz is set to True");
1341 continue;}
1342
1343 if (m_Remove_Pixel_Rx) //If m_Remove_Pixel_Rx is set to True, Pixel Rx will be skipped in the solving.
1344 if (oldPixel and (thisParameterType == 3))
1345 {ATH_MSG_INFO( "Pixel Rx DoF has been skipped in the solving because Remove_Pixel_Rx is set to True");
1346 continue;}
1347
1348 if (m_Remove_Pixel_Ry) //If m_Remove_Pixel_Ry is set to True, Pixel Ry will be skipped in the solving.
1349 if (oldPixel and (thisParameterType == 4))
1350 {ATH_MSG_INFO( "Pixel Ry DoF has been skipped in the solving because Remove_Pixel_Ry is set to True");
1351 continue;}
1352
1353 if (m_Remove_Pixel_Rz) //If m_Remove_Pixel_Rz is set to True, Pixel Rz will be skipped in the solving.
1354 if (oldPixel and (thisParameterType == 5))
1355 {ATH_MSG_INFO( "Pixel Rz DoF has been skipped in the solving because Remove_Pixel_Rz is set to True");
1356 continue;}
1357
1358 //IBL Dofs:
1359 if (m_Remove_IBL_Tx) //If m_Remove_IBL_Tx is set to True, IBL Tx will be skipped in the solving.
1360 if (ibl and (thisParameterType == 0))
1361 {ATH_MSG_INFO( "IBL Tx DoF has been skipped in the solving because Remove_IBL_Tx is set to True");
1362 continue;}
1363
1364 if (m_Remove_IBL_Ty) //If m_Remove_IBL_Ty is set to True, IBL Ty will be skipped in the solving.
1365 if (ibl and (thisParameterType == 1))
1366 {ATH_MSG_INFO( "IBL Ty DoF has been skipped in the solving because Remove_IBL_Ty is set to True");
1367 continue;}
1368
1369 if (m_Remove_IBL_Tz) //If m_Remove_IBL_Tz is set to True, IBL Tz will be skipped in the solving.
1370 if (ibl and (thisParameterType == 2))
1371 {ATH_MSG_INFO( "IBL Tz DoF has been skipped in the solving because Remove_IBL_Tz is set to True");
1372 continue;}
1373
1374 if (m_Remove_IBL_Rx) //If m_Remove_IBL_Rx is set to True, IBL Rx will be skipped in the solving.
1375 if (ibl and (thisParameterType == 3))
1376 {ATH_MSG_INFO( "IBL Rx DoF has been skipped in the solving because Remove_IBL_Rx is set to True");
1377 continue;}
1378
1379 if (m_Remove_IBL_Ry) //If m_Remove_IBL_Ry is set to True, IBL Ry will be skipped in the solving.
1380 if (ibl and (thisParameterType == 4))
1381 {ATH_MSG_INFO( "IBL Ry DoF has been skipped in the solving because Remove_IBL_Ry is set to True");
1382 continue;}
1383
1384 if (m_Remove_IBL_Rz) //If m_Remove_IBL_Rz is set to True, IBL Rz will be skipped in the solving.
1385 if (ibl and (thisParameterType == 5))
1386 {ATH_MSG_INFO( "IBL Rz DoF has been skipped in the solving because Remove_IBL_Rz is set to True");
1387 continue;}
1388
1389 if(theParameterList[i]->alignModule()->nHits() >= m_minNumHits && theParameterList[i]->alignModule()->nTracks() >= m_minNumTrks)
1390 m_activeIndices.push_back(i);
1391 }
1392 m_aNDoF = m_activeIndices.size();
1393 ATH_MSG_DEBUG("aNDoF/nDoF: "<<m_aNDoF<<"/"<<nDoF);
1394
1395 // --------------------
1396 // now do the SOLVING
1397 // --------------------
1398
1399 int info = 0;
1400
1401 // first Local solving
1402 if (m_runLocal)
1403 info = solveLocal();
1404
1405 // remove spurious modules and resize
1406 if (m_removeSpurious) {
1407
1408 ATH_MSG_INFO("Spurious removal not implemented at the moment.");
1409/* if (StatusCode::SUCCESS != spuriousRemoval()) {
1410 ATH_MSG_ERROR("Problem while trying to remove spurious. Stopping solving");
1411 return -1;
1412 }
1413
1414 // if nDoF=0, bad job...
1415 int NDoF = m_alignModuleTool->nAlignParameters();
1416 if(NDoF==0) {
1417 ATH_MSG_WARNING("Removal removed everything: NDoF=" << NDoF << " !!!");
1418 return 1;
1419 }
1420 ATH_MSG_DEBUG("NDoF: " << NDoF);
1421*/
1422 }
1423
1424 // --------------------
1425 // now Global solving
1426 switch(m_solveOption) {
1427
1428 case NONE:
1429 ATH_MSG_DEBUG("No global solving requested.");
1430 break;
1431
1432 case SOLVE_ROOT:
1433 info = solveROOT();
1434 break;
1435
1436 case SOLVE_CLHEP:
1437 info = solveCLHEP();
1438 break;
1439
1440 case SOLVE:
1441 case DIRECT_SOLVE:
1443 info = solveLapack();
1444 break;
1445
1446 case SOLVE_FAST:
1447 case DIRECT_SOLVE_FAST:
1449 break;
1450
1451 default:
1452 ATH_MSG_INFO("Unknown solving option.");
1453 info = 0;
1454 break;
1455 }
1456
1457 ATH_MSG_INFO("Return value from solving: "<<info);
1458
1459 return info;
1460 }
int nTracks() const
Definition IMatrixTool.h:88
int nHits() const
Definition IMatrixTool.h:84
Gaudi::Property< std::string > m_tfileName
Definition MatrixTool.h:226
Gaudi::Property< bool > m_AlignIBLbutNotPixel
Definition MatrixTool.h:256
Gaudi::Property< bool > m_Remove_Pixel_Ry
Definition MatrixTool.h:268
Gaudi::Property< bool > m_writeHitmap
Definition MatrixTool.h:189
Gaudi::Property< bool > m_Remove_Pixel_Rz
Definition MatrixTool.h:269
Gaudi::Property< bool > m_Remove_IBL_Rx
Definition MatrixTool.h:275
Gaudi::Property< bool > m_runLocal
Definition MatrixTool.h:201
Gaudi::Property< double > m_removeSpurious
Definition MatrixTool.h:211
Gaudi::Property< bool > m_Remove_Pixel_Tz
Definition MatrixTool.h:266
Gaudi::Property< bool > m_Remove_IBL_Rz
Definition MatrixTool.h:277
Gaudi::Property< bool > m_Remove_Pixel_Ty
Definition MatrixTool.h:265
Gaudi::Property< bool > m_Remove_IBL_Tx
Definition MatrixTool.h:272
int m_aNDoF
number of active DoF (size of m_activeIndices)
Definition MatrixTool.h:249
Gaudi::Property< bool > m_Remove_Pixel_Tx
Definition MatrixTool.h:264
Gaudi::Property< bool > m_Remove_IBL_Tz
Definition MatrixTool.h:274
Gaudi::Property< bool > m_writeMatTxt
Definition MatrixTool.h:180
Gaudi::Property< int > m_solveOption
Definition MatrixTool.h:160
Gaudi::Property< bool > m_DeactivateSCT_ECA_LastDisk
Definition MatrixTool.h:261
Gaudi::Property< bool > m_writeMat
Definition MatrixTool.h:178
void storeInTFile(const TString &filename)
Store Files in a tfile.
Gaudi::Property< double > m_softEigenmodeCut
Definition MatrixTool.h:208
Gaudi::Property< bool > m_writeTFile
Definition MatrixTool.h:196
Gaudi::Property< bool > m_Remove_Pixel_Rx
Definition MatrixTool.h:267
Gaudi::Property< bool > m_AlignPixelbutNotIBL
Definition MatrixTool.h:258
Gaudi::Property< bool > m_Remove_IBL_Ry
Definition MatrixTool.h:276
Gaudi::Property< bool > m_Remove_IBL_Ty
Definition MatrixTool.h:273
Gaudi::Property< bool > m_writeModuleNames
Definition MatrixTool.h:186
constexpr int pow(int x)
Definition conifer.h:27

◆ solveCLHEP()

int MatrixTool::solveCLHEP ( )
private

Definition at line 259 of file MatrixTool.cxx.

260 {
261 ATH_MSG_INFO("solving Global using CLHEP");
262 if(m_logStream) {
263 *m_logStream<<"*************************************************************\n";
264 *m_logStream<<"************** solving using Global method ****************\n";
265 *m_logStream<<"************** using CLHEP ****************\n";
266 *m_logStream<<"*************************************************************\n";
267 }
268
269 // start measuring time
270 clock_t starttime = clock();
271
272 DataVector<AlignPar> * alignParList = m_alignModuleTool->alignParList1D();
273 //const AlignModuleList* alignModules = m_alignModuleTool->alignModules1D();
274 for (int i=0;i<(int)alignParList->size();i++)
275 ATH_MSG_DEBUG("ap["<<i<<"]="<<(*alignParList)[i]);
276
277 int nDoF = m_alignModuleTool->nAlignParameters();
278
279 // some debugging output
280 if (msgLvl(MSG::DEBUG)) {
281 msg(MSG::DEBUG)<<"dumping matrix and vector to screen"<<endmsg;
282 for (int i=0;i<nDoF;i++)
283 for (int j=0;j<nDoF;j++)
284 //if (std::fabs((*m_bigmatrix)[i][j])>.0001)
285 msg(MSG::DEBUG)<<i<<", "<<j<<" : "<<(*m_bigmatrix)[i][j] <<endmsg;
286
287 for (int i=0;i<nDoF;i++)
288 msg(MSG::DEBUG)<<i <<" : "<<(*m_bigvector)[i]<<endmsg;
289 }
290
291 // get rescaled first and second derivatives
292 CLHEP::HepSymMatrix * d2Chi2 = new CLHEP::HepSymMatrix(m_aNDoF,0);
293 CLHEP::HepVector * dChi2 = new CLHEP::HepVector(m_aNDoF,0);
294 for (int iActive=0;iActive<m_aNDoF;iActive++) {
295 int i = m_activeIndices[iActive];
296 (*dChi2)[iActive] = (*m_bigvector)[i];
297 for (int jActive=0;jActive<m_aNDoF;jActive++) {
298 int j = m_activeIndices[jActive];
299 (*d2Chi2)[iActive][jActive] = (*m_bigmatrix)[i][j];
300 }
301 }
302
303 ATH_MSG_DEBUG("First derivatives:" << (*dChi2));
304 ATH_MSG_DEBUG("Second derivatives:" << (*d2Chi2));
305
306 CLHEP::HepSymMatrix cov(m_aNDoF,0);
307 CLHEP::HepVector delta(m_aNDoF,0);
308 CLHEP::HepVector deltafull(m_aNDoF,0);
309
310 bool status=true;
311 int ierr(0);
312 if(!m_diagonalize) {
313 // ==========================================================
314 // Run Matrix Inversion
315 ATH_MSG_INFO("Running matrix inversion");
316 if(m_logStream)
317 *m_logStream<<"Running matrix inversion\n";
318
319 cov = *d2Chi2;
320 cov.invert(ierr);
321 if(ierr>0)
322 msg(MSG::ERROR)<<"CLHEP inversion status flag = "<<ierr<<endmsg;
323 else
324 ATH_MSG_INFO("CLHEP inversion OK");
325 if(m_logStream)
326 *m_logStream<<"CLHEP inversion status flag = "<<ierr<<"\n";
327
328 // calculate corrections
329 delta = cov * (*dChi2);
330
331 // the covariance matrix is actually defined as 2 * d2Chi2^-1
332 ATH_MSG_DEBUG("Result: "<<delta);
333 ATH_MSG_DEBUG("cov: "<<cov*2);
334
335 // -----------------------
336 // calculate also for matrix and vector multiplied by factor 0.5
337 // this should make no difference if everything is correct but
338 // it can go wrong if insensitive DoF is included
339 CLHEP::HepSymMatrix cov2 = *d2Chi2 * .5;
340 // invert the matrix
341 int ierr2 = 0;
342 cov2.invert(ierr2);
343 if(ierr2>0)
344 msg(MSG::WARNING)<<"Second CLHEP inversion status flag = "<<ierr2<<endmsg;
345
346 CLHEP::HepVector delta2 = cov2 * (*dChi2) * .5;
347 for (int i=0;i<delta.num_row(); ++i)
348 if ( fabs((delta[i] - delta2[i])/delta[i]) > 1e-5 ) {
349 msg(MSG::WARNING)<<"Something's wrong with the matrix inversion: delta["<<i<<"] = "<<delta[i]<<" delta2["<<i<<"] = "<<delta2[i]<<endmsg;
350 status=false;
351 break;
352 }
353
354 if(m_logStream && (ierr2>0 || !status)) {
355 *m_logStream<<"CLHEP inversion status flag for halfed matrix = "<<ierr2<<"\n";
356 *m_logStream<<"Matrix inversion check failed\n";
357 *m_logStream<<"\n";
358 }
359 // -- end of check of matrix inversion
360 }
361 else {
362 // ==========================================================
363 // Run Diagonalization
364 ATH_MSG_INFO("Running diagonalization");
365 if(m_logStream)
366 *m_logStream<<"Running diagonalization\n";
367
368 CLHEP::HepSymMatrix D = *d2Chi2;
369 CLHEP::HepMatrix U = CLHEP::diagonalize( &D );
370
371 ATH_MSG_INFO("Diagonalization done");
372 //sold = U*sdiag*U.T.
373
374 // reorder eigenvalues ascending
375 // eigenvectors need to be reordered consistently
376 ATH_MSG_DEBUG(" Reordering eigenvalues ascending ");
377 for (int i=0; i<m_aNDoF-1; i++)
378 for (int j=i+1; j<m_aNDoF; j++)
379 if(D[j][j] < D[i][i]) {
380 // swap eigenvalues
381 double ei = D[i][i];
382 D[i][i] = D[j][j];
383 D[j][j] = ei;
384 // swap eigenvectors
385 for(int k=0;k<m_aNDoF; k++) {
386 double ev = U[k][i];
387 U[k][i] = U[k][j];
388 U[k][j] = ev;
389 }
390 }
391
392 // how do I now get the eigenvalues? this cannot be the most
393 // efficient way ... CLHEP::HepSymMatrix D = d2Chi2->similarityT( U );
394 CLHEP::HepVector eigenvector(m_aNDoF);
395
396 if(m_logStream)
397 *m_logStream<<"/------ The Eigenvalue Spectrum -------\n";
398
399 ATH_MSG_DEBUG("Calculating eigenvalues");
400 for(int imode=0; imode<m_aNDoF; ++imode) {
401
402 // get the relevant eigenvector
403 for(int irow=0; irow<m_aNDoF; ++irow)
404 eigenvector[irow] = U[irow][imode];
405
406 // calculate the eigenvalue
407 //double eigenvalue = d2Chi2->similarity( eigenvector );
408 double eigenvalue = D[imode][imode];
409 ATH_MSG_DEBUG("eigenvalue "<<eigenvalue);
410
411 double evdotb = dot(*dChi2,eigenvector);
412 CLHEP::HepVector thisdelta = evdotb/eigenvalue * eigenvector;
413 deltafull += thisdelta;
414
415 if(imode<m_modcut) {
416 ATH_MSG_INFO("skipping eigenvalue "<<imode<<" : "<<eigenvalue<<" , modcut is "<<m_modcut);
417 if(m_logStream)
418 *m_logStream<<"| skipping eigenvalue "<<eigenvalue<<"\n";
419 }
420 else if( eigenvalue < m_eigenvaluethreshold ) {
421 ATH_MSG_INFO("skipping eigenvalue "<<eigenvalue<<" , cut is "<<m_eigenvaluethreshold);
422 if(m_logStream)
423 *m_logStream<<"| skipping eigenvalue "<<eigenvalue<<"\n";
424 }
425 else {
426 if(m_logStream)
427 *m_logStream<<"| "<<eigenvalue<<"\n";
428
429 delta += thisdelta;
430
431 // this is the time consuming part
432 for(int irow=0; irow<m_aNDoF; ++irow)
433 for(int icol=0; icol<=irow; ++icol)
434 cov[irow][icol] += eigenvector[irow] * eigenvector[icol] / eigenvalue;
435 }
436 }
437
438 // the covariance matrix is actually defined as 2 * d2Chi2^-1
439
440 ATH_MSG_DEBUG("Result: "<<delta);
441 ATH_MSG_DEBUG("cov: "<<cov);
442
443 if(m_logStream)
444 *m_logStream<<"\\----- End of Eigenvalue Spectrum -----\n";
445
446 // end of diagonalization
447 // ==========================================================
448 }
449
450 // stop measuring time
451 clock_t stoptime = clock();
452 double totaltime = (stoptime-starttime)/double(CLOCKS_PER_SEC);
453 ATH_MSG_INFO("Time spent in solveCLHEP: "<<totaltime<<" s");
454
455 if(ierr==0 && status)
456 {
457 ATH_MSG_DEBUG("Alignment constants:");
458 for (int iAdof=0;iAdof<m_aNDoF;iAdof++) {
459
460 int idof = m_activeIndices[iAdof];
461 AlignPar * alignPar=(*alignParList)[idof];
462
463 double sigma = alignPar->sigma();
464 double param = -delta[iAdof] * sigma;
465 double err = std::sqrt(2.*std::fabs(cov[iAdof][iAdof])) * sigma;
466
467 ATH_MSG_DEBUG(iAdof <<" : "<< param << " +/- "<< err);
468 ATH_MSG_DEBUG("cov("<<iAdof<<")="<<cov[iAdof][iAdof]<<", sigma: "<<sigma);
469 ATH_MSG_DEBUG("init par: "<<alignPar->initPar());
470 alignPar->setPar(param,err);
471 ATH_MSG_DEBUG("set param to "<<param<<" for alignPar "<<alignPar);
472 ATH_MSG_DEBUG(*(*alignParList)[idof]);
473 }
474
475 if(m_logStream) {
477
478 // norm of first derivative
479 *m_logStream<<"norm of first derivative : "<<dChi2->norm()<<"\n";
480
481 // distance to solution
482 double dist = ( - (*d2Chi2) * deltafull + (*dChi2) ).norm();
483 *m_logStream<<"distance to solution : "<<dist<<"\n";
484
485 // calculate chi2 of the alignment change
486 double chi2 = d2Chi2->similarity(delta) * .5;
487 *m_logStream<<"delta(chi2) of the alignment change : "<<chi2<<" / "<<m_aNDoF<<"\n";
488
489 // time spent here
490 *m_logStream<<"time spent in solve : "<<totaltime<<" s\n";
491 }
492 }
493
494 delete d2Chi2;
495 delete dChi2;
496
497 return 1;
498 }
bool msgLvl(const MSG::Level lvl) const
Gaudi::Property< bool > m_diagonalize
Definition MatrixTool.h:154
int ev
Definition globals.cxx:25
dot(G, fn, nodesToHighlight=[])
Definition dot.py:5
status
Definition merge.py:16

◆ solveLapack()

int MatrixTool::solveLapack ( )
private

Definition at line 1633 of file MatrixTool.cxx.

1634 {
1635 ATH_MSG_INFO("solving Global using Lapack");
1636 if(m_logStream) {
1637 *m_logStream<<"*************************************************************\n";
1638 *m_logStream<<"************** solving using Global method ****************\n";
1639 *m_logStream<<"************** using LAPACK ****************\n";
1640 *m_logStream<<"*************************************************************\n";
1641 }
1642
1643 // get rescaled first and second derivatives
1644 AlSymMat* aBetterMat = new AlSymMat(m_aNDoF);
1645 AlVec* aBetterVec = new AlVec(m_aNDoF);
1646 for (int iActive=0;iActive<m_aNDoF;iActive++) {
1647 int i = m_activeIndices[iActive];
1648 (*aBetterVec)[iActive] = (*m_bigvector)[i];
1649 for (int jActive=0;jActive<m_aNDoF;jActive++) {
1650 int j = m_activeIndices[jActive];
1651 (*aBetterMat)[iActive][jActive] = (*m_bigmatrix)[i][j];
1652 }
1653 }
1654
1655 // normalize bigmatrix and bigvector
1656 if(m_scaleMatrix) {
1657 if(m_scale<=0.)
1658 ATH_MSG_WARNING("Scaling requested but scale not set. Not scaling matrix and vector.");
1659 else {
1660 (*aBetterVec) *= 1./m_scale;
1661 (*aBetterMat) *= 1./m_scale;
1662 }
1663 }
1664
1665 ATH_MSG_DEBUG("Now Solving alignment using lapack diagonalization routine dspev...");
1666
1667 if (m_calDet) {
1668 const double tol = 1.e-20;
1669 // compute final determinant
1670 double determ = (*aBetterMat).determinant();
1671 ATH_MSG_INFO("Determinant: " << determ);
1672 if (fabs(determ) < tol)
1673 ATH_MSG_WARNING("Matrix is singular!");
1674 }
1675
1676 // store the original matrix for checks
1677 AlSymMat * d2Chi2 = nullptr;
1679 d2Chi2 = new AlSymMat(*aBetterMat);
1680
1681 clock_t starttime = clock();
1682
1683 // declare transition matrix + vector to store eigenvalues
1684 AlMat z(m_aNDoF,m_aNDoF);
1685 AlVec w(m_aNDoF); // vector to store the eigenvalues
1686 ATH_MSG_DEBUG("MatrixTool::after z/w allocation");
1687
1688 char jobz = 'V';
1689 int info = (*aBetterMat).diagonalize(jobz,w,z);
1690 ATH_MSG_DEBUG(" info: " << info);
1691 ATH_MSG_INFO("MatrixTool::after diagonalization");
1692
1693 // stop time calculation
1694 clock_t stoptime = clock();
1695 double time_diag = (stoptime-starttime)/double(CLOCKS_PER_SEC);
1696 ATH_MSG_INFO(" - time spent diagonalizing the matrix: "<<time_diag<<" s");
1697
1698 double time_solve = 0.;
1699 if (info==0) {
1700 starttime = clock();
1701 postSolvingLapack(aBetterVec,d2Chi2,w,z,m_aNDoF);
1702 stoptime = clock();
1703 time_solve = (stoptime-starttime)/double(CLOCKS_PER_SEC);
1704 ATH_MSG_INFO(" - time spent solving the system: "<<time_solve<<" s");
1705 if(m_logStream) {
1706 *m_logStream<<"time spent for diagonalization: "<<time_diag<<" s\n";
1707 *m_logStream<<"time spent for post-solving: "<<time_solve<<" s\n";
1708 }
1709 }
1710 else {
1711 ATH_MSG_ERROR("Problem in diagonalization. Solving skipped.");
1712 if(m_logStream)
1713 *m_logStream<<"time spent for diagonalization: "<<time_diag<<" s\n";
1714 }
1715
1716 if(m_logStream) {
1717 *m_logStream<<"total time spent in solve: "<<time_diag+time_solve<<" s\n";
1718 *m_logStream<<"\n";
1719 }
1720
1721 delete d2Chi2;
1722 delete aBetterMat;
1723 delete aBetterVec;
1724
1725 // need to do this since success return value from Lapack is 0
1726 // and from solveLapack() it is 1
1727 if (info==0)
1728 info = 1;
1729
1730 return info;
1731 }
void postSolvingLapack(AlVec *dChi2, AlSymMat *d2Chi2, AlVec &w, AlMat &z, int size)
Gaudi::Property< bool > m_calDet
Definition MatrixTool.h:174

◆ solveLocal()

int MatrixTool::solveLocal ( )
private

Definition at line 501 of file MatrixTool.cxx.

502 {
503 ATH_MSG_INFO("solving using Local method");
504 if(m_logStream) {
505 *m_logStream<<"*************************************************************\n";
506 *m_logStream<<"************** solving using Local method *****************\n";
507 *m_logStream<<"*************************************************************\n";
508 }
509
510 int totalNDoF(0);
511 double totalChi2(0.);
512
513 const AlignModuleList * alignModules = m_alignModuleTool->alignModules1D();
514
515 AlignModuleList::const_iterator imod = alignModules->begin();
516 AlignModuleList::const_iterator imod_end = alignModules->end();
517 for( ; imod!=imod_end; ++imod) {
518
519 AlignModule * module = *imod;
520
521 ATH_MSG_INFO("Solving for module: "<<module->name());
522
523 DataVector<AlignPar> * alignPars = m_alignModuleTool->getAlignPars(module);
524
525 int thisNDoF = alignPars->size();
526
527 CLHEP::HepSymMatrix d2Chi2(thisNDoF,0);
528 CLHEP::HepVector dChi2(thisNDoF,0);
529 for (int i=0;i<thisNDoF;++i) {
530 int ipar = alignPars->at(i)->index();
531 dChi2[i] = (*m_bigvector)[ipar];
532 for (int j=0;j<thisNDoF;++j) {
533 int jpar = alignPars->at(j)->index();
534 d2Chi2[i][j] = (*m_bigmatrix)[ipar][jpar];
535 }
536 }
537
538 ATH_MSG_DEBUG("First derivatives:" << dChi2);
539 ATH_MSG_DEBUG("Second derivatives:" << d2Chi2);
540
541 if(module->nHits() < m_minNumHits || module->nTracks() < m_minNumTrks) {
542 ATH_MSG_INFO("Not enough hits in module \'"<<module->name()<<"\': "
543 <<module->nHits()<<" < "<<m_minNumHits<<" or "
544 <<module->nTracks()<<" < "<<m_minNumTrks
545 <<". Skipping...");
546
547 // print module summary even when solving is not done
548 if(m_logStream)
549 printModuleSolution(*m_logStream,module,nullptr);
550
551 continue;
552 }
553
554 ATH_MSG_DEBUG("First derivatives:" << dChi2);
555 ATH_MSG_DEBUG("Second derivatives:" << d2Chi2);
556
557
558 totalNDoF += thisNDoF;
559
560 CLHEP::HepSymMatrix cov(d2Chi2);
561
562 // invert the matrix
563 int ierr = 0;
564 cov.invert(ierr);
565 if(ierr>0)
566 ATH_MSG_WARNING("CLHEP inversion status flag = "<<ierr);
567 else
568 ATH_MSG_DEBUG("CLHEP inversion status flag = "<<ierr);
569
570 // calculate corrections
571 CLHEP::HepVector delta = cov * dChi2;
572 //ATH_MSG_DEBUG("d2Chi2: "<<d2Chi2);
573 //ATH_MSG_DEBUG("cov: "<<cov);
574 //ATH_MSG_DEBUG("d2Chi2*cov: "<<d2Chi2*cov);
575 //ATH_MSG_DEBUG("dChi2: "<<dChi2);
576 ATH_MSG_DEBUG("Result: "<<delta);
577
578 ATH_MSG_DEBUG("Alignment constants:");
579 for (int idof=0;idof<thisNDoF;++idof) {
580 AlignPar * alignPar = alignPars->at(idof);
581
582 double sigma = alignPar->sigma();
583 double param = -delta[idof] * sigma;
584 double err = std::sqrt(2.*std::fabs(cov[idof][idof])) * sigma;
585
586 ATH_MSG_DEBUG(idof <<" : "<< param << " +/- "<< err);
587 ATH_MSG_DEBUG("cov("<<idof<<")="<<cov[idof][idof]<<", sigma: "<<sigma);
588 ATH_MSG_DEBUG("init par: "<<alignPar->initPar());
589
590 ATH_MSG_DEBUG("Filling constants obtained using Local method");
591 alignPar->setPar(param,err);
592 ATH_MSG_DEBUG("set param to "<<param<<" for alignPar "<<alignPar);
593 ATH_MSG_DEBUG(*alignPar);
594 }
595
596 if(m_logStream) {
597 printModuleSolution(*m_logStream,module,&cov);
598
599 *m_logStream<<"CLHEP inversion status flag = "<<ierr<<"\n";
600
601 // calculate chi2 of the alignment change
602 double chi2 = d2Chi2.similarity(delta) * .5;
603 totalChi2 += chi2;
604 *m_logStream<<"delta(chi2) of the alignment change : "<<chi2<<" / "<<thisNDoF<<"\n";
605 }
606 }
607
608 if(m_logStream) {
609 *m_logStream<<"--------------------------------------------------------------------------------\n";
610 *m_logStream<<"Total delta(chi2) of the alignment change from the local method : "<<totalChi2<<" / "<<totalNDoF<<"\n";
611 *m_logStream<<"\n";
612 }
613
614 return 1;
615 }

◆ solveROOT()

int MatrixTool::solveROOT ( )
private

Definition at line 133 of file MatrixTool.cxx.

134 {
135 ATH_MSG_INFO("solving Global using ROOT");
136 if(m_logStream) {
137 *m_logStream<<"*************************************************************\n";
138 *m_logStream<<"************** solving using Global method ****************\n";
139 *m_logStream<<"************** using ROOT ****************\n";
140 *m_logStream<<"*************************************************************\n";
141 }
142
143 // start measuring time
144 clock_t starttime = clock();
145
146 DataVector<AlignPar>* alignParList = m_alignModuleTool->alignParList1D();
147 //const AlignModuleList* alignModules = m_alignModuleTool->alignModules1D();
148
149 int nDoF=m_alignModuleTool->nAlignParameters();
150
151 // some debugging output
152 if (msgLvl(MSG::VERBOSE)) {
153 msg(MSG::VERBOSE)<<"dumping matrix and vector to screen"<<endmsg;
154 for (int i=0;i<nDoF;i++)
155 for (int j=0;j<nDoF;j++)
156 //if (std::fabs((*m_bigmatrix)[i][j])>.0001)
157 msg(MSG::VERBOSE)<<i<<", "<<j<<" : "<<(*m_bigmatrix)[i][j] <<endmsg;
158
159 for (int i=0;i<nDoF;i++)
160 msg(MSG::VERBOSE)<<i <<" : "<<(*m_bigvector)[i]<<endmsg;
161 }
162
163 // get rescaled first and second derivatives
164 double * secderiv = new double[m_aNDoF*m_aNDoF];
165 double * firstderiv = new double[m_aNDoF];
166 for (int iActive=0;iActive<m_aNDoF;iActive++) {
167 int i = m_activeIndices[iActive];
168 firstderiv[iActive] = (*m_bigvector)[i];
169 for (int jActive=0;jActive<m_aNDoF;jActive++) {
170 int j = m_activeIndices[jActive];
171 secderiv[iActive*m_aNDoF+jActive] = (*m_bigmatrix)[i][j];
172 }
173 }
174
175 // attention, the dimension of matrix a and b is m_aNDoF not nDoF,
176 // this means some alignment parameters have not been calculated
177 // if the corresponding modules did not satify the select cut
178
179 TMatrixDSym a(m_aNDoF,secderiv);
180 TVectorD b(m_aNDoF,firstderiv);
181
182 if(msgLvl(MSG::DEBUG)) {
183 msg(MSG::DEBUG)<<"First derivatives:"<<endmsg;
184 b.Print();
185 msg(MSG::DEBUG)<<"Second derivatives:"<<endmsg;
186 a.Print();
187 }
188
189 TDecompBK c(a);
190 Bool_t status;
191 TMatrixDSym ainv(c.Invert(status));
192
193 TVectorD r(b.GetNrows());
194 if(status)
195 r = c.Solve(b,status);
196
197 // stop measuring time
198 clock_t stoptime = clock();
199 double totaltime = (stoptime-starttime)/double(CLOCKS_PER_SEC);
200 ATH_MSG_INFO("Time spent in solveROOT: "<<totaltime<<" s");
201
202 if(!status) {
203 msg(MSG::ERROR)<<"ROOT inversion failed"<<endmsg;
204 if(m_logStream) {
205 *m_logStream<<"ROOT inversion failed\n";
206 *m_logStream<<"\n";
207 }
208 }
209 else {
210 ATH_MSG_INFO("ROOT inversion ok");
211
212 ATH_MSG_DEBUG("Alignment constants:");
213 for (int iAdof=0;iAdof<m_aNDoF;iAdof++) {
214
215 int idof = m_activeIndices[iAdof];
216 AlignPar * alignPar=(*alignParList)[idof];
217
218 double sigma = alignPar->sigma();
219 double param = -r[iAdof] * sigma;
220 double err = std::sqrt(2.*std::fabs(ainv(iAdof,iAdof))) * sigma;
221
222 ATH_MSG_DEBUG(iAdof <<" : "<< param << " +/- "<< err);
223 ATH_MSG_DEBUG("ainv("<<iAdof<<")="<<ainv(iAdof,iAdof)<<", sigma: "<<sigma);
224 ATH_MSG_DEBUG("init par: "<<alignPar->initPar());
225 alignPar->setPar(param,err);
226 ATH_MSG_DEBUG("set param to "<<param<<" for alignPar "<<alignPar);
227 ATH_MSG_DEBUG(*(*alignParList)[idof]);
228 }
229
230 if(m_logStream)
231 {
232 *m_logStream<<"ROOT inversion ok\n";
233
235
236 // norm of first derivative
237 *m_logStream<<"norm of first derivative : "<<sqrt(b.Norm2Sqr())<<"\n";
238
239 // distance to solution
240 double dist = sqrt( ( b - (a * r) ).Norm2Sqr() );
241 *m_logStream<<"distance to solution : "<<dist<<"\n";
242
243 // calculate chi2 of the alignment change
244 double chi2 = a.Similarity(r) * .5;
245 *m_logStream<<"delta(chi2) of the alignment change : "<<chi2<<" / "<<m_aNDoF<<"\n";
246
247 // time spent here
248 *m_logStream<<"time spent in solve : "<<totaltime<<" s\n";
249 }
250 }
251
252 delete [] secderiv;
253 delete [] firstderiv;
254
255 return 1;
256 }
static Double_t a
int r
Definition globals.cxx:22

◆ solveSparseEigen()

int MatrixTool::solveSparseEigen ( )
private

Definition at line 2156 of file MatrixTool.cxx.

2157 {
2158 ATH_MSG_INFO("solving Global using SparseEigen");
2159 if(m_logStream) {
2160 *m_logStream<<"*************************************************************\n";
2161 *m_logStream<<"************** solving using Global method ****************\n";
2162 *m_logStream<<"************** using SparseEigen ****************\n";
2163 *m_logStream<<"*************************************************************\n";
2164 }
2165
2166 // start measuring time
2167 clock_t starttime = clock();
2168
2169 DataVector<AlignPar> * alignParList = m_alignModuleTool->alignParList1D();
2170
2171 AlSpaMat * ABetterMat = nullptr;
2172 bool isCopy = false;
2173 if ( dynamic_cast<AlSymMat*>(m_bigmatrix) ) {
2174 ATH_MSG_INFO("Converting Matrix Format for fast solving");
2175 ABetterMat = new AlSpaMat(*(dynamic_cast<AlSymMat*>(m_bigmatrix)));
2176 isCopy = true;
2177 }
2178 else if ( dynamic_cast<AlSpaMat*>(m_bigmatrix) ) {
2179 ATH_MSG_INFO("Matrix format native to the fast solving");
2180 ABetterMat = (dynamic_cast<AlSpaMat*>(m_bigmatrix));
2181 }
2182 else {
2183 ATH_MSG_ERROR("Cannot cast to neither AlSymMat nor AlSpaMat");
2184 return 0;
2185 }
2186
2187 ATH_MSG_DEBUG("checking active indices");
2188
2189 // use const matrix when checking for non-zero elements to avoid
2190 // filling of the whole matrix
2191 const AlSpaMat * chkMatrix = ABetterMat;
2192
2193 AlSpaMat * aBetterMat = new AlSpaMat(m_aNDoF);
2194 AlVec * aBetterVec = new AlVec(m_aNDoF);
2195
2196
2197 for (int iActive=0;iActive<m_aNDoF;iActive++) {
2198 int i = m_activeIndices[iActive];
2199 (*aBetterVec)[iActive] = (*m_bigvector)[i];
2200 for (int jActive=0;jActive<m_aNDoF;jActive++) {
2201 int j = m_activeIndices[jActive];
2202 // only fill if non-zero !!!
2203 if ( (*chkMatrix)[iActive][jActive] != 0. )
2204 (*aBetterMat)[iActive][jActive]=(*ABetterMat)[i][j];
2205 }
2206 }
2207
2208 // store original vector for cross-checks
2209 AlVec origVec(*aBetterVec);
2210
2211 ATH_MSG_DEBUG("running the solving");
2212
2213 // solve
2214 int info = (*aBetterMat).SolveWithEigen(*aBetterVec);
2215
2216 if(info == 0) {
2217 ATH_MSG_INFO("SolveWithEigen solving OK");
2218 if(m_logStream)
2219 *m_logStream<<"SolveWithEigen solving OK.\n";
2220 }
2221 else {
2222 ATH_MSG_ERROR( "SolveWithEigen error code (0 if OK) = "<<info );
2223 if(m_logStream)
2224 *m_logStream<<"SolveWithEigen error code (0 if OK) = "<<info<<"\n";
2225 }
2226
2227 if( isCopy )
2228 delete ABetterMat;
2229 ABetterMat = nullptr;
2230
2231 // stop measuring time
2232 clock_t stoptime = clock();
2233 double totaltime = (stoptime-starttime)/double(CLOCKS_PER_SEC);
2234 ATH_MSG_INFO("Time spent in SolveWithEigen: "<<totaltime<<" s");
2235
2236 ATH_MSG_DEBUG("Alignment constants:");
2237 // compute alignment corrections (translations in mm and rotations in rad)
2238 // for solveSparseEigen variances are not calculated
2239 for(int i=0; i<m_aNDoF; i++) {
2240
2241 double param = -(*aBetterVec)[i];
2242 double err = 0.;
2243
2244 int idof = m_activeIndices[i];
2245 AlignPar * alignPar=(*alignParList)[idof];
2246
2247 // undo the sigma scaling
2248 double sigma = alignPar->sigma();
2249 param *= sigma;
2250
2251 ATH_MSG_DEBUG(i <<" : "<< param << " +/- "<< err);
2252 ATH_MSG_DEBUG("sigma: "<<sigma);
2253 ATH_MSG_DEBUG("init par: "<<alignPar->initPar());
2254 alignPar->setPar(param, err);
2255 ATH_MSG_DEBUG("set param to "<<param<<" for alignPar "<<alignPar);
2256 ATH_MSG_DEBUG(*(*alignParList)[idof]);
2257 }
2258
2259 if(m_logStream) {
2260 CLHEP::HepSymMatrix * cov = nullptr;
2262
2263 // norm of first derivative
2264 *m_logStream<<"norm of first derivative : "<<origVec.norm()<<"\n";
2265
2266 // distance to solution
2267 double dist = ( (*aBetterMat) * (*aBetterVec) - origVec ).norm();
2268 *m_logStream<<"distance to solution : "<<dist<<"\n";
2269
2270 // calculate chi2 of the alignment change
2271 double chi2 = (*aBetterVec) * (*aBetterMat) * (*aBetterVec) * .5;
2272 *m_logStream<<"delta(chi2) of the alignment change : "<<chi2<<" / "<<m_aNDoF<<"\n";
2273
2274 // time spent here
2275 *m_logStream<<"time spent in solve : "<<totaltime<<" s\n";
2276 }
2277
2278 delete aBetterMat;
2279 delete aBetterVec;
2280
2281 // need to do this since success return value from Lapack is 0
2282 // and from SolveWithEigen() it is 1
2283 if (info==0)
2284 info = 1;
2285
2286 return info;
2287 }

◆ spuriousRemoval()

StatusCode MatrixTool::spuriousRemoval ( )
private

Definition at line 1734 of file MatrixTool.cxx.

1735 {
1736
1737 // copied from SiGlobalChi2Algs
1738 ATH_MSG_DEBUG("in spuriousRemoval");
1739
1740 // compute determinant before resizing
1741 if (m_calDet) {
1742 const double tol = 1.e-20;
1743 double determ = m_bigmatrix->determinant();
1744 ATH_MSG_INFO("Determinant: " << determ);
1745 if (std::fabs(determ) < tol)
1746 ATH_MSG_WARNING("Matrix is singular!");
1747 }
1748
1749 // fill vector with modules that need to be removed
1750 int fillvecmods=fillVecMods();
1751 if (fillvecmods==0) {
1752
1753 //ATH_MSG_INFO(" No resize needed (NhitsCut = "
1754 // << m_hitscut << ")");
1755
1756 if (msgLvl(MSG::DEBUG)) {
1757 //m_bigmatrix->Write("bigmatrix.txt", false, m_wSqMatrix, m_scale,
1758 // MatVersion);
1759 //m_bigvector->Write("bigvector.txt", false, m_scale, m_modcodemap,
1760 // VecVersion, m_alignProcessLevel, m_fitParam);
1761 }
1762
1763 return StatusCode::SUCCESS;
1764 }
1765 else if (fillvecmods==2)
1766 return StatusCode::FAILURE;
1767
1768 /* this is a bit difficult to implement for now....
1769
1770 // remove matrix/vector elements
1771 int cont=0;
1772 ModuleIndexMap::iterator itcode;
1773 ATH_MSG_INFO("Eliminating module...");
1774
1775
1776
1777 for (std::vector<int>::const_iterator it=m_dropmods.begin();
1778 it!= m_dropmods.end(); ++it) {
1779
1780 itcode = m_modcodemap.find((*it));
1781
1782 if (itcode == m_modcodemap.end()) {
1783 ATH_MSG_WARNING("Could not find module " << *it << " in map.");
1784 return StatusCode::FAILURE;
1785 }
1786
1787 ATH_MSG_INFO(" - Removing mcode: " << (itcode->second)
1788 << " (old index: " << (*it) << " -> new index: " << (*it)-cont << ")");
1789
1790 m_bigmatrix->RemoveModule((*it)-cont);
1791 m_bigvector->RemoveModule((*it)-cont);
1792 cont++;
1793 }
1794 ATH_MSG_INFO("Modules removed from the Matrix and from the Vector Successfully!");
1795 ATH_MSG_DEBUG("(DoF: " << nDoF << ")");
1796
1797 // resizing...
1798 ATH_MSG_INFO("Resizing the bigvector in memory...");
1799 m_bigvector->reSize(nDoF-6*m_dropmods.size());
1800 ATH_MSG_INFO("Resizing the bigmatrix in memory...");
1801 m_bigmatrix->reSize(nDoF-6*m_dropmods.size());
1802
1803 ATH_MSG_INFO(m_dropmods.size() << " modules eliminated from the matrix, i.e, "
1804 << 6*m_dropmods.size() << " DoFs");
1805 ATH_MSG_INFO(nDoF/6 - m_dropmods.size() << " modules to align (" << nDoF-6*m_dropmods.size() << " DoFs)");
1806 ATH_MSG_INFO("New bigmatrix size is: " << m_bigmatrix->size());
1807 ATH_MSG_INFO("New bigvector size is: " << (*m_bigvector).size());
1808
1809 NDoF = nDoF-6*(m_dropmods.size());
1810
1811 // resize (update) nDoF variable
1812 nDoF = NDoF;
1813
1814 // Resizing vectors to store results
1815 m_alignPar->reSize(nDoF);
1816 m_alignSqErr->reSize(nDoF);
1817
1818 // Fill the mapping module map and updating m_modcodemap
1819 UpdateModMap();
1820 */
1821
1822 return StatusCode::SUCCESS;
1823 }
static int fillVecMods()

◆ storeInTFile()

void MatrixTool::storeInTFile ( const TString & filename)

Store Files in a tfile.

Definition at line 730 of file MatrixTool.cxx.

731 {
732 //Store reults in a single TFile....
733 //Including Matrix Vector Hitmap.. Soluton EVs etc.
734
735 ATH_MSG_DEBUG("Writing Results to a TFile");
736
737 DataVector<AlignPar> * alignParList = m_alignModuleTool->alignParList1D();
738 int nDoF = alignParList->size();
739
740 TMatrixDSparse* myTMatrix = m_bigmatrix->makeTMatrix();
741
742 ATH_MSG_DEBUG( "Created TMatrixDSparse" );
743
744
745 double *val = new double[nDoF];
746 for (int i=0;i<nDoF;i++) {
747 val[i] = (*m_bigvector)[i];
748 }
749
750 TVectorD myTVector(nDoF, val);
751 delete [] val;
752
753 ATH_MSG_DEBUG( "Created TVectorD" );
754
755
756 const AlignModuleList * moduleList = m_alignModuleTool->alignModules1D();
757 int nModules = moduleList->size();
758
759 double *hitmapA = new double[nModules];
760 double *hitmapB = new double[nModules];
761 AlignModuleList::const_iterator imod = moduleList->begin();
762 AlignModuleList::const_iterator imod_end = moduleList->end();
763 int index(0);
764 for(; imod != imod_end; ++imod) {
765 AlignModule * module = *imod;
766 hitmapA[index] = (double)module->nHits();
767 hitmapB[index] = (double)module->nTracks();
768 index++;
769 }
770
771 TVectorD hitmapHits(nModules, hitmapA);
772 TVectorD hitmapTracks(nModules, hitmapB);
773
774 delete [] hitmapA;
775 delete [] hitmapB;
776
777
778
779 TFile myFile(filename,"recreate");
780 hitmapHits.Write("Hits");
781 hitmapTracks.Write("Tracks");
782 myTMatrix->Write("Matrix");
783 myTVector.Write("Vector");
784
785 TVectorD scale(1, &m_scale) ;
786 scale.Write("Scale");
787
788
789 double *moduleInfoA = new double[nDoF];//unsigned long long
790 double *dofInfoA = new double[nDoF];//int
791
792 if (sizeof(unsigned long long) != sizeof(double))
793 ATH_MSG_ERROR("Module Identifiers will not be saved. sizeof(double)!=sizeof(ulonglong)");
794 else{
795
796 DataVector<AlignPar>* alignPars = m_alignModuleTool->alignParList1D();
797 for (int i=0;i<(int)alignPars->size();i++) {
798 //Do a direct memory copy to store unsigned long long in the momory of a double
799 double target;
800 uint64_t id = (*alignPars)[i]->alignModule()->identify().get_compact();
801 memcpy(&target, &id, sizeof(target));
802 moduleInfoA[i]=target;
803 //moduleInfoB[i]=(*alignPars)[i]->alignModule()->name();
804 uint64_t dof = (*alignPars)[i]->paramType();
805 memcpy(&target, &dof, sizeof(target));
806 dofInfoA[i]=target;
807 }
808
809 TVectorD moduleIDs(nDoF, moduleInfoA) ;
810 TVectorD moduleDoFs(nDoF,dofInfoA);
811 delete [] moduleInfoA;
812 moduleIDs.Write("ModuleID");
813 moduleDoFs.Write("dof");
814 }
815
816 myFile.Write();
817 myFile.Close();
818
819 delete myTMatrix;
820 ATH_MSG_DEBUG("Finshed writing TFILE");
821
822 }

◆ sysInitialize()

virtual StatusCode AthCommonDataStore< AthCommonMsg< AlgTool > >::sysInitialize ( )
overridevirtualinherited

Perform system initialization for an algorithm.

We override this to declare all the elements of handle key arrays at the end of initialization. See comments on updateVHKA.

Reimplemented in asg::AsgMetadataTool, AthCheckedComponent< AthAlgTool >, and AthCheckedComponent<::AthAlgTool >.

◆ sysStart()

virtual StatusCode AthCommonDataStore< AthCommonMsg< AlgTool > >::sysStart ( )
overridevirtualinherited

Handle START transition.

We override this in order to make sure that conditions handle keys can cache a pointer to the conditions container.

◆ updateVHKA()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::updateVHKA ( Gaudi::Details::PropertyBase & )
inlineinherited

Definition at line 308 of file AthCommonDataStore.h.

308 {
309 // debug() << "updateVHKA for property " << p.name() << " " << p.toString()
310 // << " size: " << m_vhka.size() << endmsg;
311 for (auto &a : m_vhka) {
313 for (auto k : keys) {
314 k->setOwner(this);
315 }
316 }
317 }
std::vector< SG::VarHandleKeyArray * > m_vhka

◆ writeHitmap()

void MatrixTool::writeHitmap ( )
private

Definition at line 2290 of file MatrixTool.cxx.

2291 {
2292 ATH_MSG_INFO("writing the hitmap to file");
2293
2294 const AlignModuleList * moduleList = m_alignModuleTool->alignModules1D();
2295 int nModules = moduleList->size();
2296
2297 AlMat hitmap(nModules,2);
2298 AlignModuleList::const_iterator imod = moduleList->begin();
2299 AlignModuleList::const_iterator imod_end = moduleList->end();
2300 int index(0);
2301 for(; imod != imod_end; ++imod) {
2302 AlignModule * module = *imod;
2303 hitmap[index][0] = module->nHits();
2304 hitmap[index][1] = module->nTracks();
2305 index++;
2306 }
2307
2308 // Set Path for the hitmap matrix
2309 hitmap.SetPathBin(m_pathbin.value()+m_prefixName.value());
2310 hitmap.SetPathTxt(m_pathtxt.value()+m_prefixName.value());
2311
2312 StatusCode sc = hitmap.Write("hitmap.bin",true); // write the hitmap matrix
2313
2314 if (sc!=StatusCode::SUCCESS)
2315 ATH_MSG_ERROR("Problem writing hitmap matrix");
2316
2317 if (m_writeHitmapTxt) {
2318 sc = hitmap.Write("hitmap.txt", false, 0);
2319 if (sc!=StatusCode::SUCCESS)
2320 ATH_MSG_ERROR("Problem writing hitmap matrix to text file");
2321 }
2322
2323 ATH_MSG_DEBUG("hitmap written to: "<< m_pathbin.value()+m_prefixName.value() <<"hitmap.bin (.txt)");
2324 }
Gaudi::Property< bool > m_writeHitmapTxt
Definition MatrixTool.h:191

Member Data Documentation

◆ m_activeIndices

std::vector<int> Trk::MatrixTool::m_activeIndices {}
private

vector of indices which pass the min-hits cut

Definition at line 248 of file MatrixTool.h.

248{};

◆ m_Align_db_step

Gaudi::Property<float> Trk::MatrixTool::m_Align_db_step
private
Initial value:
{this, "AlignCorrDBStep", 10.,
"corr in the diagonal basis step for the third pass in the auto weak mode removal method"}

Definition at line 171 of file MatrixTool.h.

171 {this, "AlignCorrDBStep", 10.,
172 "corr in the diagonal basis step for the third pass in the auto weak mode removal method"};

◆ m_AlignIBLbutNotPixel

Gaudi::Property<bool> Trk::MatrixTool::m_AlignIBLbutNotPixel {this, "AlignIBLbutNotPixel", false}
private

Definition at line 255 of file MatrixTool.h.

256{this, "AlignIBLbutNotPixel", false};

◆ m_alignModuleMap

std::map<int,int> Trk::IMatrixTool::m_alignModuleMap
privateinherited

Definition at line 105 of file IMatrixTool.h.

◆ m_alignModuleTool

PublicToolHandle<IAlignModuleTool> Trk::MatrixTool::m_alignModuleTool
private
Initial value:
{
this, "AlignModuleTool", "Trk::AlignModuleTool/AlignModuleTool"}

Definition at line 142 of file MatrixTool.h.

142 {
143 this, "AlignModuleTool", "Trk::AlignModuleTool/AlignModuleTool"};

◆ m_AlignPixelbutNotIBL

Gaudi::Property<bool> Trk::MatrixTool::m_AlignPixelbutNotIBL {this, "AlignPixelbutNotIBL", false}
private

Definition at line 257 of file MatrixTool.h.

258{this, "AlignPixelbutNotIBL", false};

◆ m_aNDoF

int Trk::MatrixTool::m_aNDoF = 0
private

number of active DoF (size of m_activeIndices)

Definition at line 249 of file MatrixTool.h.

◆ m_bigmatrix

AlSymMatBase* Trk::MatrixTool::m_bigmatrix = nullptr
private

matrix to contain second derivative terms to be used for alignment

Definition at line 146 of file MatrixTool.h.

◆ m_bigvector

AlVec* Trk::MatrixTool::m_bigvector = nullptr
private

vector to contain first derivative terms to be used for alignment

Definition at line 149 of file MatrixTool.h.

◆ m_calculateFullCovariance

Gaudi::Property<double> Trk::MatrixTool::m_calculateFullCovariance
private
Initial value:
{this, "CalculateFullCovariance", true,
"calculate full covariance matrix for Lapack"}

Definition at line 214 of file MatrixTool.h.

215 {this, "CalculateFullCovariance", true,
216 "calculate full covariance matrix for Lapack"};

◆ m_calDet

Gaudi::Property<bool> Trk::MatrixTool::m_calDet
private
Initial value:
{this, "MatrixDet", false,
"compute bigmatrix's determinant ?"}

Definition at line 174 of file MatrixTool.h.

174 {this, "MatrixDet", false,
175 "compute bigmatrix's determinant ?"};

◆ m_DeactivateSCT_ECA_LastDisk

Gaudi::Property<bool> Trk::MatrixTool::m_DeactivateSCT_ECA_LastDisk {this, "DeactivateSCT_ECA_LastDisk", false}
private

Definition at line 260 of file MatrixTool.h.

261{this, "DeactivateSCT_ECA_LastDisk", false};

◆ m_detStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< AlgTool > >::m_detStore
privateinherited

Pointer to StoreGate (detector store by default).

Definition at line 393 of file AthCommonDataStore.h.

◆ m_diagonalize

Gaudi::Property<bool> Trk::MatrixTool::m_diagonalize
private
Initial value:
{this, "Diagonalize", true,
"run diagonalization instead of inversion"}

Definition at line 154 of file MatrixTool.h.

154 {this, "Diagonalize", true,
155 "run diagonalization instead of inversion"};

◆ m_eigenvalueStep

Gaudi::Property<float> Trk::MatrixTool::m_eigenvalueStep
private
Initial value:
{this, "EigenvalueStep", 1e3,
"eigenvalue step for the second pass in the automatic weak mode removal method"}

Definition at line 169 of file MatrixTool.h.

169 {this, "EigenvalueStep", 1e3,
170 "eigenvalue step for the second pass in the automatic weak mode removal method"};

◆ m_eigenvaluethreshold

Gaudi::Property<double> Trk::MatrixTool::m_eigenvaluethreshold {this, "EigenvalueThreshold", 0., "cut on the minimum eigenvalue"}
private

Definition at line 156 of file MatrixTool.h.

157{this, "EigenvalueThreshold", 0., "cut on the minimum eigenvalue"};

◆ m_evtStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< AlgTool > >::m_evtStore
privateinherited

Pointer to StoreGate (event store by default).

Definition at line 390 of file AthCommonDataStore.h.

◆ m_inputHitmapFiles

Gaudi::Property<std::vector<std::string> > Trk::MatrixTool::m_inputHitmapFiles
private
Initial value:
{this, "InputHitmapFiles", {"hitmap.bin"},
"input binary files containing the hitmaps"}

Definition at line 240 of file MatrixTool.h.

241 {this, "InputHitmapFiles", {"hitmap.bin"},
242 "input binary files containing the hitmaps"};

◆ m_inputMatrixFiles

Gaudi::Property<std::vector<std::string> > Trk::MatrixTool::m_inputMatrixFiles
private
Initial value:
{this, "InputMatrixFiles", {"matrix.bin"},
"input binary files containing matrix terms"}

Definition at line 233 of file MatrixTool.h.

234 {this, "InputMatrixFiles", {"matrix.bin"},
235 "input binary files containing matrix terms"};

◆ m_inputTFiles

Gaudi::Property<std::vector<std::string> > Trk::MatrixTool::m_inputTFiles
private
Initial value:
{this, "InputTFiles", {"AlignmentTFile.root"},
"input binary files containing matrix terms"}

Definition at line 244 of file MatrixTool.h.

245 {this, "InputTFiles", {"AlignmentTFile.root"},
246 "input binary files containing matrix terms"};

◆ m_inputVectorFiles

Gaudi::Property<std::vector<std::string> > Trk::MatrixTool::m_inputVectorFiles
private
Initial value:
{this, "InputVectorFiles", {"vector.bin"},
"input binary files containing vector terms"}

Definition at line 236 of file MatrixTool.h.

237 {this, "InputVectorFiles", {"vector.bin"},
238 "input binary files containing vector terms"};

◆ m_logStream

std::ostream* Trk::IMatrixTool::m_logStream
protectedinherited

logfile output stream

Definition at line 98 of file IMatrixTool.h.

◆ m_maxReadErrors

Gaudi::Property<int> Trk::MatrixTool::m_maxReadErrors
private
Initial value:
{this, "MaxReadErrors", 10,
"maximum number of reading TFile errors"}

Definition at line 251 of file MatrixTool.h.

251 {this, "MaxReadErrors", 10,
252 "maximum number of reading TFile errors"};

◆ m_minNumHits

Gaudi::Property<int> Trk::MatrixTool::m_minNumHits
private
Initial value:
{this, "MinNumHitsPerModule", 0,
"cut on the minimum number of hits per module"}

Definition at line 163 of file MatrixTool.h.

163 {this, "MinNumHitsPerModule", 0,
164 "cut on the minimum number of hits per module"};

◆ m_minNumTrks

Gaudi::Property<int> Trk::MatrixTool::m_minNumTrks
private
Initial value:
{this, "MinNumTrksPerModule", 0,
"cut on the minimum number of tracks per module"}

Definition at line 165 of file MatrixTool.h.

165 {this, "MinNumTrksPerModule", 0,
166 "cut on the minimum number of tracks per module"};

◆ m_modcut

Gaudi::Property<int> Trk::MatrixTool::m_modcut
private
Initial value:
{this, "ModCut", 0,
"cut on the weak modes which number is <par_modcut"}

Definition at line 161 of file MatrixTool.h.

161 {this, "ModCut", 0,
162 "cut on the weak modes which number is <par_modcut"};

◆ m_nentries

int Trk::IMatrixTool::m_nentries
privateinherited

Definition at line 106 of file IMatrixTool.h.

◆ m_nHits

int Trk::IMatrixTool::m_nHits
protectedinherited

Definition at line 100 of file IMatrixTool.h.

◆ m_nMeasurements

int Trk::IMatrixTool::m_nMeasurements
protectedinherited

Definition at line 102 of file IMatrixTool.h.

◆ m_nTracks

int Trk::IMatrixTool::m_nTracks
protectedinherited

Definition at line 101 of file IMatrixTool.h.

◆ m_pathbin

Gaudi::Property<std::string> Trk::MatrixTool::m_pathbin
private
Initial value:
{this, "PathBinName", "./",
"path binary files (in/out)"}

Definition at line 218 of file MatrixTool.h.

218 {this, "PathBinName", "./",
219 "path binary files (in/out)"};

◆ m_pathtxt

Gaudi::Property<std::string> Trk::MatrixTool::m_pathtxt
private
Initial value:
{this, "PathTxtName", "./",
"path ascii files (in/out)"}

Definition at line 220 of file MatrixTool.h.

220 {this, "PathTxtName", "./",
221 "path ascii files (in/out)"};

◆ m_prefixName

Gaudi::Property<std::string> Trk::MatrixTool::m_prefixName
private
Initial value:
{this, "PrefixName", "",
"prefix string to filenames"}

Definition at line 222 of file MatrixTool.h.

222 {this, "PrefixName", "",
223 "prefix string to filenames"};

◆ m_pullcut

Gaudi::Property<float> Trk::MatrixTool::m_pullcut
private
Initial value:
{this, "PullCut", 1.0,
"pull cut for the automatic weak mode removal method"}

Definition at line 167 of file MatrixTool.h.

167 {this, "PullCut", 1.0,
168 "pull cut for the automatic weak mode removal method"};

◆ m_readHitmaps

Gaudi::Property<bool> Trk::MatrixTool::m_readHitmaps
private
Initial value:
{this, "ReadHitmaps", false,
"accumulate hitymap from files"}

Definition at line 193 of file MatrixTool.h.

193 {this, "ReadHitmaps", false,
194 "accumulate hitymap from files"};

◆ m_readTFiles

Gaudi::Property<bool> Trk::MatrixTool::m_readTFiles
private
Initial value:
{this, "ReadTFile", false,
"if True then files will be read from TFiles instead of Binary files"}

Definition at line 198 of file MatrixTool.h.

198 {this, "ReadTFile", false,
199 "if True then files will be read from TFiles instead of Binary files"};

◆ m_Remove_IBL_Rx

Gaudi::Property<bool> Trk::MatrixTool::m_Remove_IBL_Rx {this, "Remove_IBL_Rx", false}
private

Definition at line 275 of file MatrixTool.h.

275{this, "Remove_IBL_Rx", false};

◆ m_Remove_IBL_Ry

Gaudi::Property<bool> Trk::MatrixTool::m_Remove_IBL_Ry {this, "Remove_IBL_Ry", false}
private

Definition at line 276 of file MatrixTool.h.

276{this, "Remove_IBL_Ry", false};

◆ m_Remove_IBL_Rz

Gaudi::Property<bool> Trk::MatrixTool::m_Remove_IBL_Rz {this, "Remove_IBL_Rz", false}
private

Definition at line 277 of file MatrixTool.h.

277{this, "Remove_IBL_Rz", false};

◆ m_Remove_IBL_Tx

Gaudi::Property<bool> Trk::MatrixTool::m_Remove_IBL_Tx {this, "Remove_IBL_Tx", false}
private

Definition at line 272 of file MatrixTool.h.

272{this, "Remove_IBL_Tx", false};

◆ m_Remove_IBL_Ty

Gaudi::Property<bool> Trk::MatrixTool::m_Remove_IBL_Ty {this, "Remove_IBL_Ty", false}
private

Definition at line 273 of file MatrixTool.h.

273{this, "Remove_IBL_Ty", false};

◆ m_Remove_IBL_Tz

Gaudi::Property<bool> Trk::MatrixTool::m_Remove_IBL_Tz {this, "Remove_IBL_Tz", false}
private

Definition at line 274 of file MatrixTool.h.

274{this, "Remove_IBL_Tz", false};

◆ m_Remove_Pixel_Rx

Gaudi::Property<bool> Trk::MatrixTool::m_Remove_Pixel_Rx {this, "Remove_Pixel_Rx", false}
private

Definition at line 267 of file MatrixTool.h.

267{this, "Remove_Pixel_Rx", false};

◆ m_Remove_Pixel_Ry

Gaudi::Property<bool> Trk::MatrixTool::m_Remove_Pixel_Ry {this, "Remove_Pixel_Ry", false}
private

Definition at line 268 of file MatrixTool.h.

268{this, "Remove_Pixel_Ry", false};

◆ m_Remove_Pixel_Rz

Gaudi::Property<bool> Trk::MatrixTool::m_Remove_Pixel_Rz {this, "Remove_Pixel_Rz", false}
private

Definition at line 269 of file MatrixTool.h.

269{this, "Remove_Pixel_Rz", false};

◆ m_Remove_Pixel_Tx

Gaudi::Property<bool> Trk::MatrixTool::m_Remove_Pixel_Tx {this, "Remove_Pixel_Tx", false}
private

Definition at line 264 of file MatrixTool.h.

264{this, "Remove_Pixel_Tx", false};

◆ m_Remove_Pixel_Ty

Gaudi::Property<bool> Trk::MatrixTool::m_Remove_Pixel_Ty {this, "Remove_Pixel_Ty", false}
private

Definition at line 265 of file MatrixTool.h.

265{this, "Remove_Pixel_Ty", false};

◆ m_Remove_Pixel_Tz

Gaudi::Property<bool> Trk::MatrixTool::m_Remove_Pixel_Tz {this, "Remove_Pixel_Tz", false}
private

Definition at line 266 of file MatrixTool.h.

266{this, "Remove_Pixel_Tz", false};

◆ m_removeSpurious

Gaudi::Property<double> Trk::MatrixTool::m_removeSpurious
private
Initial value:
{this, "RemoveSpurious", false,
"run spurious removal"}

Definition at line 211 of file MatrixTool.h.

211 {this, "RemoveSpurious", false,
212 "run spurious removal"};

◆ m_runLocal

Gaudi::Property<bool> Trk::MatrixTool::m_runLocal
private
Initial value:
{this, "RunLocalMethod", true,
"Run solving using Local method"}

Definition at line 201 of file MatrixTool.h.

201 {this, "RunLocalMethod", true,
202 "Run solving using Local method"};

◆ m_scalaMatName

Gaudi::Property<std::string> Trk::MatrixTool::m_scalaMatName {this, "ScalapackMatrixName", "eigenvectors.bin", "Scalapack matrix name"}
private

Definition at line 228 of file MatrixTool.h.

229{this, "ScalapackMatrixName", "eigenvectors.bin", "Scalapack matrix name"};

◆ m_scalaVecName

Gaudi::Property<std::string> Trk::MatrixTool::m_scalaVecName {this, "ScalapackVectorName", "eigenvalues.bin", "Scalapack vector name"}
private

Definition at line 230 of file MatrixTool.h.

231{this, "ScalapackVectorName", "eigenvalues.bin", "Scalapack vector name"};

◆ m_scale

double Trk::MatrixTool::m_scale = -1.
private

scale for big matrix and vector normalization

Definition at line 204 of file MatrixTool.h.

◆ m_scaleMatrix

Gaudi::Property<bool> Trk::MatrixTool::m_scaleMatrix
private
Initial value:
{this, "ScaleMatrix", false,
"scale matrix by number of hits before solving"}

Definition at line 205 of file MatrixTool.h.

205 {this, "ScaleMatrix", false,
206 "scale matrix by number of hits before solving"};

◆ m_softEigenmodeCut

Gaudi::Property<double> Trk::MatrixTool::m_softEigenmodeCut
private
Initial value:
{this, "SoftEigenmodeCut", 0.,
"add constant to diagonal to effectively cut on weak eigenmodes"}

Definition at line 208 of file MatrixTool.h.

208 {this, "SoftEigenmodeCut", 0.,
209 "add constant to diagonal to effectively cut on weak eigenmodes"};

◆ m_solveOption

Gaudi::Property<int> Trk::MatrixTool::m_solveOption {this, "SolveOption", NONE, "solving option"}
private

Definition at line 159 of file MatrixTool.h.

160{this, "SolveOption", NONE, "solving option"};

◆ m_tfileName

Gaudi::Property<std::string> Trk::MatrixTool::m_tfileName {this, "TFileName", "AlignmentTFile.root", "prefix string to filenames"}
private

Definition at line 225 of file MatrixTool.h.

226{this, "TFileName", "AlignmentTFile.root", "prefix string to filenames"};

◆ m_useSparse

Gaudi::Property<bool> Trk::MatrixTool::m_useSparse {this, "UseSparse", false}
private

flag to use AlSpaMat for the big matrix (default is AlSymMat)

Definition at line 152 of file MatrixTool.h.

152{this, "UseSparse", false};

◆ m_varHandleArraysDeclared

bool AthCommonDataStore< AthCommonMsg< AlgTool > >::m_varHandleArraysDeclared
privateinherited

Definition at line 399 of file AthCommonDataStore.h.

◆ m_vhka

std::vector<SG::VarHandleKeyArray*> AthCommonDataStore< AthCommonMsg< AlgTool > >::m_vhka
privateinherited

Definition at line 398 of file AthCommonDataStore.h.

◆ m_writeEigenMat

Gaudi::Property<bool> Trk::MatrixTool::m_writeEigenMat
private
Initial value:
{this, "WriteEigenMat", true,
"write eigenvalues and eigenvectors into files ?"}

Definition at line 182 of file MatrixTool.h.

182 {this, "WriteEigenMat", true,
183 "write eigenvalues and eigenvectors into files ?"};

◆ m_writeEigenMatTxt

Gaudi::Property<bool> Trk::MatrixTool::m_writeEigenMatTxt
private
Initial value:
{this, "WriteEigenMatTxt", true,
"also write eigenvalues and eigenvectors into txt files ?"}

Definition at line 184 of file MatrixTool.h.

184 {this, "WriteEigenMatTxt", true,
185 "also write eigenvalues and eigenvectors into txt files ?"};

◆ m_writeHitmap

Gaudi::Property<bool> Trk::MatrixTool::m_writeHitmap
private
Initial value:
{this, "WriteHitmap", false,
"write hitmap into file"}

Definition at line 189 of file MatrixTool.h.

189 {this, "WriteHitmap", false,
190 "write hitmap into file"};

◆ m_writeHitmapTxt

Gaudi::Property<bool> Trk::MatrixTool::m_writeHitmapTxt
private
Initial value:
{this, "WriteHitmapTxt", false,
"write hitmap into text file"}

Definition at line 191 of file MatrixTool.h.

191 {this, "WriteHitmapTxt", false,
192 "write hitmap into text file"};

◆ m_writeMat

Gaudi::Property<bool> Trk::MatrixTool::m_writeMat
private
Initial value:
{this, "WriteMat", true,
"write big matrix and vector into files ?"}

Definition at line 178 of file MatrixTool.h.

178 {this, "WriteMat", true,
179 "write big matrix and vector into files ?"};

◆ m_writeMatTxt

Gaudi::Property<bool> Trk::MatrixTool::m_writeMatTxt
private
Initial value:
{this, "WriteMatTxt", true,
"also write big matrix and vector into txt files ?"}

Definition at line 180 of file MatrixTool.h.

180 {this, "WriteMatTxt", true,
181 "also write big matrix and vector into txt files ?"};

◆ m_writeModuleNames

Gaudi::Property<bool> Trk::MatrixTool::m_writeModuleNames
private
Initial value:
{this, "WriteModuleNames", false,
"write module name instead of Identifier to vector file"}

Definition at line 186 of file MatrixTool.h.

186 {this, "WriteModuleNames", false,
187 "write module name instead of Identifier to vector file"};

◆ m_writeTFile

Gaudi::Property<bool> Trk::MatrixTool::m_writeTFile
private
Initial value:
{this, "WriteTFile", false,
"write out files to a root file"}

Definition at line 196 of file MatrixTool.h.

196 {this, "WriteTFile", false,
197 "write out files to a root file"};

◆ m_wSqMatrix

Gaudi::Property<bool> Trk::MatrixTool::m_wSqMatrix
private
Initial value:
{this, "WriteSquareMatrix", false,
"write a triangular matrix by default (true: square format) ?"}

Definition at line 176 of file MatrixTool.h.

176 {this, "WriteSquareMatrix", false,
177 "write a triangular matrix by default (true: square format) ?"};

The documentation for this class was generated from the following files: