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
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
std::unique_ptr< AlSymMatBasem_bigmatrix
 matrix to contain second derivative terms to be used for alignment
std::unique_ptr< AlVecm_bigvector
 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 56 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 59 of file MatrixTool.h.

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

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 ( )
virtualdefault

Virtual destructor.

Member Function Documentation

◆ accumulateFromBinaries()

bool MatrixTool::accumulateFromBinaries ( )

accumulates derivates from binary files

Definition at line 626 of file MatrixTool.cxx.

627 {
628
629
630 DataVector<AlignPar>* alignParList = m_alignModuleTool->alignParList1D();
631 int nDoF=alignParList->size();
632
633 std::map<int,unsigned long long> modIndexMap;
634 float dummyVersion(0.);
635 double totalscale=0.;
636 for (int ivec=0;ivec<(int)m_inputVectorFiles.size();ivec++) {
637
638 ATH_MSG_DEBUG("Reading vector "<<ivec<<" from file "<<m_inputVectorFiles[ivec]);
639
640 AlVec newVector(nDoF);
641 std::map<int,unsigned long long> newModIndexMap;
642 newVector.SetPathBin(m_pathbin.value()+m_prefixName.value());
643 newVector.SetPathTxt(m_pathtxt.value()+m_prefixName.value());
644 double scale=0;
645 StatusCode sc = newVector.ReadPartial(m_inputVectorFiles[ivec],scale,newModIndexMap,dummyVersion);
646 totalscale += scale;
647 if (sc==StatusCode::FAILURE) {
648 msg(MSG::FATAL)<<"Problem reading vector from "<<m_inputVectorFiles[ivec]<<endmsg;
649 return false;
650 }
651 if (newVector.size()!=m_bigvector->size()) {
652 msg(MSG::FATAL) <<"vector wrong size! newVector size "<<newVector.size()
653 <<", bigvector size "<<m_bigvector->size()<<endmsg;
654 return false;
655 }
656
657 // check modIndexMaps to make sure they are the same
658 if (ivec==0)
659 modIndexMap = std::move(newModIndexMap);
660 else if (modIndexMap!=newModIndexMap) {
661 msg(MSG::FATAL)<<"module index maps don't agree!"<<endmsg;
662 return false;
663 }
664 if (ivec>0)
665 *m_bigvector += newVector;
666 else
667 *m_bigvector = std::move(newVector);
668 }
669
670 m_scale = totalscale;
671
672 auto symBigMatrix=dynamic_cast<AlSymMat*>(m_bigmatrix.get());
673 auto spaBigMatrix=dynamic_cast<AlSpaMat*>(m_bigmatrix.get());
674
675
676 for (int imat=0;imat<(int)m_inputMatrixFiles.size();imat++) {
677 ATH_MSG_DEBUG("Reading matrix "<<imat<<" from file "<<m_inputMatrixFiles[imat]);
678
679 // create new matrix to read data from current file
680 int nDoF=modIndexMap.size();
681 bool triang;
683 if (symBigMatrix) {
684 AlSymMat newMatrix(nDoF);
685 sc = newMatrix.Read(m_inputMatrixFiles[imat],nDoF,triang,dummyVersion);
686 if (sc==StatusCode::SUCCESS)
687 *symBigMatrix += newMatrix;
688 }
689 else {
690 if (!spaBigMatrix) {
691 throw std::logic_error("Unhandled matrix type");
692 }
693
694 AlSpaMat newMatrix(nDoF);
695 sc = newMatrix.Read(m_inputMatrixFiles[imat],nDoF,triang,dummyVersion);
696
697 if (sc==StatusCode::SUCCESS) {
698 if (imat>0)
699 *spaBigMatrix += newMatrix;
700 else
701 *spaBigMatrix = newMatrix;
702 }
703 }
704
705 if (sc==StatusCode::FAILURE) {
706 msg(MSG::FATAL)<<"problem reading matrix from "<<m_inputMatrixFiles[imat]<<endmsg;
707 return false;
708 }
709
710 if (!m_useSparse && triang==m_wSqMatrix) {
711 ATH_MSG_WARNING("matrix not expected format! Changing m_wSqMatrix to "<<!triang);
712 m_wSqMatrix=!triang;
713 }
714
715 }
716
717 // accumulate hitmap from hitmap files
718 if(m_readHitmaps)
719 readHitmaps();
720
721 return true;
722 }
#define endmsg
#define ATH_MSG_DEBUG(x,...)
#define ATH_MSG_WARNING(x,...)
static Double_t sc
MsgStream & msg() const
size_type size() const noexcept
Returns the number of elements in the collection.
std::unique_ptr< AlSymMatBase > m_bigmatrix
matrix to contain second derivative terms to be used for alignment
Definition MatrixTool.h:147
Gaudi::Property< std::string > m_prefixName
Definition MatrixTool.h:223
Gaudi::Property< bool > m_readHitmaps
Definition MatrixTool.h:194
Gaudi::Property< bool > m_wSqMatrix
Definition MatrixTool.h:177
double m_scale
scale for big matrix and vector normalization
Definition MatrixTool.h:205
Gaudi::Property< std::vector< std::string > > m_inputVectorFiles
Definition MatrixTool.h:238
Gaudi::Property< std::vector< std::string > > m_inputMatrixFiles
Definition MatrixTool.h:235
PublicToolHandle< IAlignModuleTool > m_alignModuleTool
Definition MatrixTool.h:143
std::unique_ptr< AlVec > m_bigvector
vector to contain first derivative terms to be used for alignment
Definition MatrixTool.h:150
Gaudi::Property< std::string > m_pathbin
Definition MatrixTool.h:219
Gaudi::Property< bool > m_useSparse
flag to use AlSpaMat for the big matrix (default is AlSymMat)
Definition MatrixTool.h:153
Gaudi::Property< std::string > m_pathtxt
Definition MatrixTool.h:221
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 613 of file MatrixTool.cxx.

614 {
615
616 if(m_readTFiles){
617 ATH_MSG_INFO("Info to obtained from from TFiles");
618 return accumulateFromTFiles();
619 }else{
620 ATH_MSG_INFO("Info to obtained from from Binary files");
621 return accumulateFromBinaries();
622 }
623 }
#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:199

◆ accumulateFromTFiles()

bool MatrixTool::accumulateFromTFiles ( )

Store Files in a tfile.

Definition at line 804 of file MatrixTool.cxx.

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

1455 {
1456 (*m_bigvector)[irow] += firstderiv;
1457 }

◆ addFirstDerivatives() [1/2]

void MatrixTool::addFirstDerivatives ( AlVec * vector)
virtual

adds first derivative to vector

Implements Trk::IMatrixTool.

Definition at line 1434 of file MatrixTool.cxx.

1435 {
1436 }

◆ 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 1444 of file MatrixTool.cxx.

1445 {
1446 }

◆ 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 1460 of file MatrixTool.cxx.

1461 {
1462 (*m_bigmatrix)[irow][icol] += secondderiv;
1463 }

◆ addSecondDerivatives() [1/2]

void MatrixTool::addSecondDerivatives ( AlSymMatBase * matrix)
virtual

adds second derivatives to matrix

Implements Trk::IMatrixTool.

Definition at line 1439 of file MatrixTool.cxx.

1440 {
1441 }

◆ 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 1449 of file MatrixTool.cxx.

1450 {
1451 }

◆ allocateMatrix()

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

allocates memory for big matrix and big vector

Implements Trk::IMatrixTool.

Definition at line 97 of file MatrixTool.cxx.

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

◆ 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 1575 of file MatrixTool.cxx.

1576 {
1577 return 0;
1578 }

◆ finalize()

StatusCode MatrixTool::finalize ( )

initialize

Definition at line 88 of file MatrixTool.cxx.

89 {
90 ATH_MSG_DEBUG("finalize() of MatrixTool");
91
92 return StatusCode::SUCCESS;
93 }

◆ initialize()

StatusCode MatrixTool::initialize ( )

initialize

Definition at line 67 of file MatrixTool.cxx.

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

◆ 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 1773 of file MatrixTool.cxx.

1774 {
1775 ATH_MSG_DEBUG("in postSolvinglapack()");
1776
1777 if( z.ncol() != size) {
1778 msg(MSG::ERROR)<<"Eigenvector matrix has incorrect size : "<<z.ncol()<<" != "<<size<<endmsg;
1779 return;
1780 }
1781
1782 if( (int)m_activeIndices.size() != size) {
1783 msg(MSG::ERROR)<<"Number of active parameters is incorrect : "<<m_activeIndices.size()<<" != "<<size<<endmsg;
1784 return;
1785 }
1786
1787 // Compute bigvector in diagonal basis (Vb = Ut * bigvector)
1788 AlVec D(size);
1789 D = z*(*dChi2);
1790
1791 if (m_writeEigenMat) {
1792
1793 ATH_MSG_INFO("writing the eigenvectors in a matrix: "<< z.nrow() << "x" << z.ncol());
1794
1795 // Set Path for the z matrix (eigenvector matrix)
1796 z.SetPathBin(m_pathbin.value()+m_prefixName.value());
1797 z.SetPathTxt(m_pathtxt.value()+m_prefixName.value());
1798
1799 ATH_MSG_INFO("writing the eigenvector matrix: "<< m_scalaMatName);
1800 ATH_MSG_DEBUG("matrix will be in: "<< m_pathbin.value()+m_prefixName.value()+m_scalaMatName.value());
1801
1802 StatusCode sc = z.Write(m_scalaMatName, true); // write the eigenvector matrix
1803
1804 if (sc!=StatusCode::SUCCESS)
1805 msg(MSG::ERROR)<<"Problem writing eigenvector matrix"<<endmsg;
1806
1807 // Set Path for the w matrix (eigenvalues matrix - diagonal bigmatrix)
1808 w.SetPathBin(m_pathbin.value()+m_prefixName.value());
1809 w.SetPathTxt(m_pathtxt.value()+m_prefixName.value());
1810
1811 ATH_MSG_INFO("writing the eigenvectors in a vector: "<< w.size());
1812 ATH_MSG_INFO("writing the eigenvalues vector (diagonal bigmatrix): "<< m_scalaVecName);
1813 ATH_MSG_DEBUG("vector will be in: "<< m_pathbin.value()+m_prefixName.value()+m_scalaVecName.value());
1814
1815 sc = w.WriteEigenvalueVec(m_scalaVecName, true); // write the eigenvalues vecor
1816
1817 if (sc!=StatusCode::SUCCESS)
1818 msg(MSG::ERROR)<<"Problem writing eigenvector matrix"<<endmsg;
1819
1820 if (m_writeEigenMatTxt) {
1821 sc = z.Write("eigenvectors.txt", false);
1822 if (sc!=StatusCode::SUCCESS)
1823 msg(MSG::ERROR)<<"Problem writing eigenvector matrix to text file"<<endmsg;
1824 sc = w.WriteEigenvalueVec("eigenvalues.txt", false);
1825 if (sc!=StatusCode::SUCCESS)
1826 msg(MSG::ERROR)<<"Problem writing eigenvalue vector to text file"<<endmsg;
1827 }
1828
1829 }
1830
1831 // Set eigenvalue thresholds
1832 const double eigenvalue_threshold = 1e-19;
1833
1834 // weak mode removal
1835 if (m_modcut == -1) {
1836
1837 ATH_MSG_INFO(" Starting the automatic Weak Mode Removal method");
1838
1839 // create a pull vector for the alignment corrections in diagonal basis (db_pulls)
1840 //int nDoF=m_alignModuleTool->nAlignParameters();
1841 AlVec* Align_db = new AlVec(size);
1842 AlVec* Align_error_db = new AlVec(size);
1843 AlVec* AlignPull = new AlVec(size);
1844 ATH_MSG_DEBUG("AlignPull vector size is: "<< (*AlignPull).size());
1845
1846 m_modcut = 0;
1847 bool wm_stop = false;
1848
1849 // -----------------------------------------------------------------------
1850 // First pass: removing weak modes because of large pull
1851 // compute alignment pulls for corrections in diagonal basis (db)
1852 for(int i=0; i<size; i++) {
1853
1854 (*Align_db)[i] = (-D[i]/w[i]);
1855 if(m_scale<=0.)[[unlikely]]{
1856 ATH_MSG_WARNING("postSolvingLapack: Scaling requested but scale not set. Not scaling matrix and vector.");
1857 } else {
1858 (*Align_error_db)[i] = sqrt(1.0/w[i]/m_scale);
1859 }
1860
1861 if (w[i]<eigenvalue_threshold) {
1862 ATH_MSG_INFO(" + EigenMode " << i
1863 << " removed as eigenvalue lower than the threshold " << eigenvalue_threshold
1864 << ": " << w[i]);
1865 (*AlignPull)[i] = 0.0;
1866 ++m_modcut;
1867 }
1868 else
1869 (*AlignPull)[i] = (*Align_db)[i] / (*Align_error_db)[i];
1870
1871 ATH_MSG_DEBUG(i << ". AlignPar: " << (*Align_db)[i] << " +- " << (*Align_error_db)[i]
1872 << " (pull: " << (*AlignPull)[i] << ") ; w[i]: " << w[i]);
1873 }
1874 ATH_MSG_INFO(" +++ Weak Mode removal ++ stop after mode "<< m_modcut << " (First pass)");
1875 // -----------------------------------------------------------------------
1876
1877 // -----------------------------------------------------------------------
1878 // Second pass
1879 // if the error is greater than the correction -> cut this mode
1880 for(int i=m_modcut; (i<size && !wm_stop); i++) {
1881
1882 // if the error is greater than the correction -> cut this mode
1883 if (fabs((*AlignPull)[i])<m_pullcut) {
1884 ATH_MSG_INFO(" + EigenMode " << i
1885 << " removed as pull is lower than " << m_pullcut << ": "
1886 << (*AlignPull)[i]);
1887 ++m_modcut;
1888 }
1889 else
1890 wm_stop = true;
1891
1892 }
1893 ATH_MSG_INFO(" +++ Weak Mode removal ++ stop after mode "<< m_modcut << " (Second pass)");
1894 // -----------------------------------------------------------------------
1895
1896 wm_stop = false;
1897
1898 // ----------------------------------------------------------------------
1899 // Third pass
1900 // Check if the next eigenvalues is far away. If it is two orders of
1901 // magnitude bigger remove also this mode and allow the search
1902 for(int i=m_modcut; (i<size && !wm_stop); i++) {
1903
1904 // if the next eigenvalues is far away -> cut this mode
1905 if (m_eigenvalueStep*w[i]<w[i+1]) {
1906 ATH_MSG_INFO(" + EigenMode " << i
1907 << " removed as diff between eigenvalues, " << w[i] << " and " << w[i+1]
1908 << ", is greater than " << m_eigenvalueStep);
1909 ++m_modcut;
1910 }
1911 else
1912 wm_stop = true;
1913
1914 }
1915 ATH_MSG_INFO(" +++ Weak Mode removal ++ stop after mode "<< m_modcut << " (Third pass)");
1916 // -----------------------------------------------------------------------
1917
1918 wm_stop = false;
1919
1920 // -----------------------------------------------------------------------
1921 // Fourth pass
1922 // Check if the next eigenvalues is far away. If it is two orders of
1923 // magnitude bigger remove also this mode and allow the search
1924 for(int i=m_modcut; (i<size && !wm_stop); i++) {
1925
1926 // if the next eigenvalues is far away -> cut this mode
1927 if ( fabs((*Align_db)[i]) > m_Align_db_step*fabs((*Align_db)[i+1]) ) {
1928 ATH_MSG_INFO(" + EigenMode " << i
1929 << " removed as diff between corrections, " << w[i] << " and " << w[i+1]
1930 << ", is greater than "
1931 << m_Align_db_step);
1932 ++m_modcut;
1933 }
1934 else
1935 wm_stop = true;
1936
1937 }
1938 ATH_MSG_INFO(" +++ Weak Mode removal ++ stop after mode "<< m_modcut << " (Fourth pass)");
1939 // -----------------------------------------------------------------------------------------------------
1940
1941 // Free memory and clear the pointer to
1942 // prevent using invalid memory reference
1943 delete Align_db;
1944 delete Align_error_db;
1945 delete AlignPull;
1946 Align_db = nullptr;
1947 Align_error_db = nullptr;
1948 AlignPull = nullptr;
1949
1950 } // end of if(m_modcut == -1)
1951
1952 // Save some stuff to debug purposes
1953 /*
1954 if (m_storeDia) {
1955 std::string path = m_pathtxt+m_prefixName+"align_dia";
1956 std::fstream orthogon(path.c_str(), std::ios::out);
1957 orthogon.setf(std::ios::fixed);
1958 orthogon.setf(std::ios::showpoint);
1959 orthogon.precision(6);
1960
1961 orthogon << std::setw(10)
1962 << "--------------------------------------------------------------------------------"
1963 << std::endl;
1964 orthogon << std::setw(10) << " ModeCut = " << m_modcut << std::endl;
1965 orthogon << std::setw(10)
1966 << "--------------------------------------------------------------------------------"
1967 << std::endl;
1968
1969 orthogon << std::setw(10) << "mode"
1970 << std::setw(20) << "eigenmode"
1971 << std::setw(20) << "eigenmode error"
1972 << std::setw(25) << "eigenvalue"
1973 << std::endl;
1974
1975 for( int m=0; m<size; m++) {
1976
1977 // mode
1978 orthogon << std::setw(10) << m;
1979
1980 // eigenmode (db)
1981 if( w[m]>1.0e-15) orthogon << std::setw(20) << -D[m]/w[m];
1982 else orthogon << std::setw(20) << 0.0;
1983
1984 // error eigenmode (error_db)
1985 if( w[m]>1.0e-15) orthogon << std::setw(20) << sqrt(1.0/w[m]/m_scale);
1986 else orthogon << std::setw(20) << 0.0;
1987
1988 // eigenvalues
1989 orthogon << std::setw(25) << w[m] << std::endl;
1990 }
1991 orthogon.close();
1992 } // end store align_dia.txt
1993 */
1994
1995 AlVec delta(size);
1996 AlVec deltafull(size);
1997 AlVec errSq(size);
1998
1999 // full covariance matrix
2000 std::unique_ptr<CLHEP::HepSymMatrix> cov;
2002 // Warning ! The matrix can be huge!
2003 // This can lead to memory problems
2004 cov = std::make_unique<CLHEP::HepSymMatrix>(size,0);
2005
2006 if(m_logStream)
2007 *m_logStream<<"/------ The Eigenvalue Spectrum -------\n";
2008
2009 for (int i=0;i<size;i++) {
2010 AlVec thisdelta(size);
2011 for(int j=0;j<size;j++)
2012 thisdelta[j] = z[i][j] * (-D[i]/w[i]);
2013 deltafull += thisdelta;
2014
2015 ATH_MSG_DEBUG("eigenvalue "<<w[i]);
2016 if( i<m_modcut ) {
2017 ATH_MSG_INFO("skipping eigenvalue "<<w[i]<<" , modcut is "<<m_modcut);
2018 if(m_logStream)
2019 *m_logStream<<"| skipping eigenvalue "<<w[i]<<"\n";
2020 }
2021 else if( w[i] < m_eigenvaluethreshold ) {
2022 ATH_MSG_INFO("skipping eigenvalue "<<w[i]<<" , cut is "<<m_eigenvaluethreshold);
2023 if(m_logStream)
2024 *m_logStream<<"| skipping eigenvalue "<<w[i]<<"\n";
2025 }
2026 else {
2027 if(m_logStream)
2028 *m_logStream<<"| "<<w[i]<<"\n";
2029
2030 delta += thisdelta;
2031 for(int j=0;j<size;j++) {
2032 errSq[j] += z[i][j] * z[i][j] / w[i];
2034 for(int k=0;k<=j;k++)
2035 (*cov)[j][k] += z[i][j] * z[i][k] / w[i];
2036 }
2037 }
2038 }
2039 }
2040
2041 if(m_logStream)
2042 *m_logStream<<"\\----- End of Eigenvalue Spectrum -----\n";
2043
2044 ATH_MSG_DEBUG("Alignment constants:");
2045
2046 DataVector<AlignPar> * alignParList = m_alignModuleTool->alignParList1D();
2047
2048 // compute alignment corrections (translations in mm and rotations in rad) and their variances
2049 for(int i=0; i<size; i++) {
2050
2051 double param = delta[i];
2052 double err = sqrt(2.*std::fabs(errSq[i]));
2053
2054 int idof = m_activeIndices[i];
2055 AlignPar * alignPar=(*alignParList)[idof];
2056
2057 // undo the sigma scaling
2058 double sigma = alignPar->sigma();
2059
2060 param *= sigma;
2061 err *= sigma;
2062
2063 // undo normalization scaling of error
2064 if(m_scaleMatrix && m_scale>0.)
2065 err /= sqrt(m_scale);
2066
2067 ATH_MSG_DEBUG(i <<" : "<< param << " +/- "<< err);
2068 ATH_MSG_DEBUG("cov("<<i<<")="<<errSq[i]<<", sigma: "<<sigma);
2069 ATH_MSG_DEBUG("init par: "<<alignPar->initPar());
2070 alignPar->setPar(param, err);
2071 ATH_MSG_DEBUG("set param to "<<param<<" for alignPar "<<alignPar);
2072 ATH_MSG_DEBUG(*(*alignParList)[idof]);
2073 }
2074
2075 if(m_logStream) {
2077
2078 // norm of first derivative
2079 double norm1st = dChi2->norm();
2080 if(m_scaleMatrix && m_scale>0.) // undo normalization scaling
2081 norm1st *= m_scale;
2082 *m_logStream<<"norm of first derivative : "<<norm1st<<"\n";
2083
2084 if(d2Chi2) {
2085 // distance to solution
2086 double dist = ( (*d2Chi2) * deltafull + (*dChi2) ).norm();
2087 if(m_scaleMatrix && m_scale>0.) // undo normalization scaling
2088 dist *= m_scale;
2089 *m_logStream<<"distance to solution : "<<dist<<"\n";
2090
2091 // calculate chi2 of the alignment change
2092 double chi2 = delta * (*d2Chi2) * delta * .5;
2093 if(m_scaleMatrix && m_scale>0.) // undo normalization scaling
2094 chi2 *= m_scale;
2095 *m_logStream<<"delta(chi2) of the alignment change : "<<chi2<<" / "<<size<<"\n";
2096 }
2097 }
2098
2099 }
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:185
std::vector< int > m_activeIndices
vector of indices which pass the min-hits cut
Definition MatrixTool.h:249
Gaudi::Property< double > m_calculateFullCovariance
Definition MatrixTool.h:216
Gaudi::Property< std::string > m_scalaVecName
Definition MatrixTool.h:232
Gaudi::Property< float > m_eigenvalueStep
Definition MatrixTool.h:170
Gaudi::Property< float > m_Align_db_step
Definition MatrixTool.h:172
Gaudi::Property< bool > m_scaleMatrix
Definition MatrixTool.h:206
Gaudi::Property< bool > m_writeEigenMat
Definition MatrixTool.h:183
Gaudi::Property< int > m_modcut
Definition MatrixTool.h:162
Gaudi::Property< double > m_eigenvaluethreshold
Definition MatrixTool.h:158
Gaudi::Property< std::string > m_scalaMatName
Definition MatrixTool.h:230
Gaudi::Property< float > m_pullcut
Definition MatrixTool.h:168
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 125 of file MatrixTool.cxx.

126 {
127 }

◆ printGlobalSolution() [1/2]

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

Definition at line 1466 of file MatrixTool.cxx.

1467 {
1468 const AlignModuleList * alignModules = m_alignModuleTool->alignModules1D();
1469
1470 AlignModuleList::const_iterator imod = alignModules->begin();
1471 AlignModuleList::const_iterator imod_end = alignModules->end();
1472 for( ; imod!=imod_end; ++imod) {
1473 AlignModule * module = *imod;
1474
1475 DataVector<AlignPar> * alignPars = m_alignModuleTool->getAlignPars(module);
1476 int thisNDoF = alignPars->size();
1477
1478 // fill local covariance matrix
1479 std::unique_ptr<CLHEP::HepSymMatrix> covsub;;
1480 if(cov && module->nHits() >= m_minNumHits && module->nTracks() >= m_minNumTrks) {
1481 covsub = std::make_unique<CLHEP::HepSymMatrix>(thisNDoF,0);
1482 for (int i=0;i<thisNDoF;++i) {
1483 int ipar = alignPars->at(i)->index();
1484 double sigma_i = alignPars->at(i)->sigma();
1485
1486 std::vector<int>::iterator itActive = std::find(m_activeIndices.begin(),m_activeIndices.end(),ipar);
1487 if( itActive == m_activeIndices.end() )
1488 continue;
1489 int iActive = std::distance(m_activeIndices.begin(),itActive);
1490
1491 for (int j=0;j<=i;++j) {
1492 int jpar = alignPars->at(j)->index();
1493 double sigma_j = alignPars->at(j)->sigma();
1494
1495 std::vector<int>::iterator jtActive = std::find(m_activeIndices.begin(),m_activeIndices.end(),jpar);
1496 if( jtActive == m_activeIndices.end() )
1497 continue;
1498 int jActive = std::distance(m_activeIndices.begin(),jtActive);
1499
1500 (*covsub)[i][j] = (*cov)[iActive][jActive] * sigma_i * sigma_j;
1501 }
1502 }
1503 }
1504
1505 printModuleSolution(os,module,covsub.get());
1506
1507 }
1508 os << "--------------------------------------------------------------------------------" << std::endl;
1509 }
const T * at(size_type n) const
Access an element, as an rvalue.
Gaudi::Property< int > m_minNumHits
Definition MatrixTool.h:164
void printModuleSolution(std::ostream &os, const AlignModule *module, const CLHEP::HepSymMatrix *cov) const
Gaudi::Property< int > m_minNumTrks
Definition MatrixTool.h:166

◆ printGlobalSolution() [2/2]

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

Definition at line 1512 of file MatrixTool.cxx.

1513 {
1514 std::unique_ptr<CLHEP::HepSymMatrix> cov;
1515 if(cov0) {
1516 int nsize = cov0->GetNrows();
1517 cov = std::make_unique<CLHEP::HepSymMatrix>(nsize,0);
1518
1519 for(int i=0; i<nsize; i++)
1520 for(int j=0; j<=i; j++)
1521 (*cov)[i][j] = (*cov0)[i][j];
1522 }
1523
1524 printGlobalSolution(os,cov.get());
1525
1526 }

◆ printModuleSolution()

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

Definition at line 1529 of file MatrixTool.cxx.

1530 {
1531 os << "--------------------------------------------------------------------------------" << std::endl;
1532 os << "Alignment parameters for module: " << module->name() << std::endl;
1533 os << "Number of tracks passing: " << module->nTracks() << std::endl;
1534 if(m_minNumHits>0 && module->nHits()<m_minNumHits) {
1535 os << "Number of hits too small: "<<module->nHits()<<" < "<<m_minNumHits<<" Skipping the module\n";
1536 return;
1537 }
1538 if(m_minNumTrks>0 && module->nTracks()<m_minNumTrks) {
1539 os << "Number of tracks too small: "<<module->nTracks()<<" < "<<m_minNumTrks<<" Skipping the module\n";
1540 return;
1541 }
1542 os << "Number of hits seen: " << module->nHits() << std::endl;
1543 os << "Number of tracks seen: " << module->nTracks() << std::endl;
1544
1545 DataVector<AlignPar> * alignPars = m_alignModuleTool->getAlignPars(module);
1546 int thisNDoF = alignPars->size();
1547
1548 if(alignPars->empty())
1549 os << "No active parameters" << std::endl;
1550 else
1551 {
1552 // output alignment parameters and errors
1553 DataVector<AlignPar>::const_iterator ipar = alignPars->begin();
1554 DataVector<AlignPar>::const_iterator ipar_end = alignPars->end();
1555 for ( ; ipar != ipar_end; ++ipar) {
1556 const AlignPar * par = *ipar;
1557 os << std::format("{:<10}{:<12.5g} +/- {:<12.5g}\n",
1558 par->dumpType(),
1559 par->par(),
1560 par->err());
1561 }
1562
1563 if(cov) {
1564 // calculate local correlation matrix
1565 CLHEP::HepSymMatrix corrsub(thisNDoF,0);
1566 for(int irow=0; irow<thisNDoF; ++irow)
1567 for(int icol=0; icol<=irow; ++icol)
1568 corrsub[irow][icol] = (*cov)[irow][icol] / sqrt((*cov)[irow][irow] * (*cov)[icol][icol]);
1569 os << "Local correlation matrix: " << corrsub << std::flush;
1570 }
1571 }
1572 }
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 2273 of file MatrixTool.cxx.

2274 {
2275 ATH_MSG_INFO("read hitmaps from files");
2276
2277 const AlignModuleList * moduleList = m_alignModuleTool->alignModules1D();
2278 int nModules = moduleList->size();
2279
2280 AlMat hitmap(nModules,2);
2281 int nFiles = (int)m_inputHitmapFiles.size();
2282 for(int imap=0;imap<nFiles; imap++) {
2283 AlMat nextmap(nModules,2);
2284
2285 ATH_MSG_INFO("Reading hitmap "<<imap<<" from file "<<m_inputHitmapFiles[imap]);
2286
2287 if(nextmap.ReadScalaPack(m_inputHitmapFiles[imap]).isFailure()) {
2288 ATH_MSG_WARNING("Problem reading hitmap from \'"<<m_inputHitmapFiles[imap]<<"\'. Skipping.");
2289 continue;
2290 }
2291
2292 if(nextmap.nrow()!=nModules || nextmap.ncol()!=2) {
2293 ATH_MSG_WARNING("Matrix in file \'"<<m_inputHitmapFiles[imap]<<"\' has wrong size ("
2294 <<nextmap.nrow()<<" x "<<nextmap.ncol()<<"), should be ("<<nModules<<" x 2). Skipping.");
2295 continue;
2296 }
2297
2298 hitmap += nextmap;
2299 }
2300
2301 AlignModuleList::const_iterator imod = moduleList->begin();
2302 AlignModuleList::const_iterator imod_end = moduleList->end();
2303 int index = 0;
2304 int totalhits = 0;
2305 for(; imod != imod_end; ++imod) {
2306 AlignModule * module = *imod;
2307 //if ((int)hitmap[index][0]!=(int)module->identify().get_identifier32().get_compact()) ATH_MSG_ERROR("bad module identifier");
2308 //module->setIdentifier((Identifier)hitmap[index][0]);
2309 module->setNHits((int)hitmap[index][0]);
2310 module->setNTracks((int)hitmap[index][1]);
2311 totalhits += (int)hitmap[index][0];
2312 index++;
2313 }
2314
2315 m_nHits = totalhits;
2316 m_nTracks = 0;
2317 m_nMeasurements = 0;
2318
2319 ATH_MSG_INFO("Hitmap accumulated from "<<nFiles<<" files with total of "<<totalhits<<" hits.");
2320 }
Gaudi::Property< std::vector< std::string > > m_inputHitmapFiles
Definition MatrixTool.h:242

◆ 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 1113 of file MatrixTool.cxx.

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

◆ solveCLHEP()

int MatrixTool::solveCLHEP ( )
private

Definition at line 255 of file MatrixTool.cxx.

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

1582 {
1583 ATH_MSG_INFO("solving Global using Lapack");
1584 if(m_logStream) {
1585 *m_logStream<<"*************************************************************\n";
1586 *m_logStream<<"************** solving using Global method ****************\n";
1587 *m_logStream<<"************** using LAPACK ****************\n";
1588 *m_logStream<<"*************************************************************\n";
1589 }
1590
1591 // get rescaled first and second derivatives
1592 AlSymMat* aBetterMat = new AlSymMat(m_aNDoF);
1593 AlVec* aBetterVec = new AlVec(m_aNDoF);
1594 for (int iActive=0;iActive<m_aNDoF;iActive++) {
1595 int i = m_activeIndices[iActive];
1596 (*aBetterVec)[iActive] = (*m_bigvector)[i];
1597 for (int jActive=0;jActive<m_aNDoF;jActive++) {
1598 int j = m_activeIndices[jActive];
1599 (*aBetterMat)[iActive][jActive] = (*m_bigmatrix)[i][j];
1600 }
1601 }
1602
1603 // normalize bigmatrix and bigvector
1604 if(m_scaleMatrix) {
1605 if(m_scale<=0.)
1606 ATH_MSG_WARNING("Scaling requested but scale not set. Not scaling matrix and vector.");
1607 else {
1608 (*aBetterVec) *= 1./m_scale;
1609 (*aBetterMat) *= 1./m_scale;
1610 }
1611 }
1612
1613 ATH_MSG_DEBUG("Now Solving alignment using lapack diagonalization routine dspev...");
1614
1615 if (m_calDet) {
1616 const double tol = 1.e-20;
1617 // compute final determinant
1618 double determ = (*aBetterMat).determinant();
1619 ATH_MSG_INFO("Determinant: " << determ);
1620 if (fabs(determ) < tol)
1621 ATH_MSG_WARNING("Matrix is singular!");
1622 }
1623
1624 // store the original matrix for checks
1625 std::unique_ptr<AlSymMat> d2Chi2 ;
1627 d2Chi2 = std::make_unique<AlSymMat>(*aBetterMat);
1628
1629 clock_t starttime = clock();
1630
1631 // declare transition matrix + vector to store eigenvalues
1632 AlMat z(m_aNDoF,m_aNDoF);
1633 AlVec w(m_aNDoF); // vector to store the eigenvalues
1634 ATH_MSG_DEBUG("MatrixTool::after z/w allocation");
1635
1636 char jobz = 'V';
1637 int info = (*aBetterMat).diagonalize(jobz,w,z);
1638 ATH_MSG_DEBUG(" info: " << info);
1639 ATH_MSG_INFO("MatrixTool::after diagonalization");
1640
1641 // stop time calculation
1642 clock_t stoptime = clock();
1643 double time_diag = (stoptime-starttime)/double(CLOCKS_PER_SEC);
1644 ATH_MSG_INFO(" - time spent diagonalizing the matrix: "<<time_diag<<" s");
1645
1646 double time_solve = 0.;
1647 if (info==0) {
1648 starttime = clock();
1649 postSolvingLapack(aBetterVec,d2Chi2.get(),w,z,m_aNDoF);
1650 stoptime = clock();
1651 time_solve = (stoptime-starttime)/double(CLOCKS_PER_SEC);
1652 ATH_MSG_INFO(" - time spent solving the system: "<<time_solve<<" s");
1653 if(m_logStream) {
1654 *m_logStream<<"time spent for diagonalization: "<<time_diag<<" s\n";
1655 *m_logStream<<"time spent for post-solving: "<<time_solve<<" s\n";
1656 }
1657 }
1658 else {
1659 ATH_MSG_ERROR("Problem in diagonalization. Solving skipped.");
1660 if(m_logStream)
1661 *m_logStream<<"time spent for diagonalization: "<<time_diag<<" s\n";
1662 }
1663
1664 if(m_logStream) {
1665 *m_logStream<<"total time spent in solve: "<<time_diag+time_solve<<" s\n";
1666 *m_logStream<<"\n";
1667 }
1668
1669 delete aBetterMat;
1670 delete aBetterVec;
1671
1672 // need to do this since success return value from Lapack is 0
1673 // and from solveLapack() it is 1
1674 if (info==0)
1675 info = 1;
1676
1677 return info;
1678 }
void postSolvingLapack(AlVec *dChi2, AlSymMat *d2Chi2, AlVec &w, AlMat &z, int size)
Gaudi::Property< bool > m_calDet
Definition MatrixTool.h:175

◆ solveLocal()

int MatrixTool::solveLocal ( )
private

Definition at line 496 of file MatrixTool.cxx.

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

◆ solveROOT()

int MatrixTool::solveROOT ( )
private

Definition at line 130 of file MatrixTool.cxx.

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

◆ solveSparseEigen()

int MatrixTool::solveSparseEigen ( )
private

Definition at line 2102 of file MatrixTool.cxx.

2103 {
2104 ATH_MSG_INFO("solving Global using SparseEigen");
2105 if(m_logStream) {
2106 *m_logStream<<"*************************************************************\n";
2107 *m_logStream<<"************** solving using Global method ****************\n";
2108 *m_logStream<<"************** using SparseEigen ****************\n";
2109 *m_logStream<<"*************************************************************\n";
2110 }
2111
2112 // start measuring time
2113 clock_t starttime = clock();
2114
2115 DataVector<AlignPar> * alignParList = m_alignModuleTool->alignParList1D();
2116
2117 AlSpaMat * ABetterMat = nullptr;
2118 bool isCopy = false;
2119 if ( dynamic_cast<AlSymMat*>(m_bigmatrix.get()) ) {
2120 ATH_MSG_INFO("Converting Matrix Format for fast solving");
2121 ABetterMat = new AlSpaMat(*(dynamic_cast<AlSymMat*>(m_bigmatrix.get())));
2122 isCopy = true;
2123 }
2124 else if ( dynamic_cast<AlSpaMat*>(m_bigmatrix.get()) ) {
2125 ATH_MSG_INFO("Matrix format native to the fast solving");
2126 ABetterMat = (dynamic_cast<AlSpaMat*>(m_bigmatrix.get()));
2127 }
2128 else {
2129 ATH_MSG_ERROR("Cannot cast to neither AlSymMat nor AlSpaMat");
2130 return 0;
2131 }
2132
2133 ATH_MSG_DEBUG("checking active indices");
2134
2135 // use const matrix when checking for non-zero elements to avoid
2136 // filling of the whole matrix
2137 const AlSpaMat * chkMatrix = ABetterMat;
2138
2139 AlSpaMat * aBetterMat = new AlSpaMat(m_aNDoF);
2140 AlVec * aBetterVec = new AlVec(m_aNDoF);
2141
2142
2143 for (int iActive=0;iActive<m_aNDoF;iActive++) {
2144 int i = m_activeIndices[iActive];
2145 (*aBetterVec)[iActive] = (*m_bigvector)[i];
2146 for (int jActive=0;jActive<m_aNDoF;jActive++) {
2147 int j = m_activeIndices[jActive];
2148 // only fill if non-zero !!!
2149 if ( (*chkMatrix)[iActive][jActive] != 0. )
2150 (*aBetterMat)[iActive][jActive]=(*ABetterMat)[i][j];
2151 }
2152 }
2153
2154 // store original vector for cross-checks
2155 AlVec origVec(*aBetterVec);
2156
2157 ATH_MSG_DEBUG("running the solving");
2158
2159 // solve
2160 int info = (*aBetterMat).SolveWithEigen(*aBetterVec);
2161
2162 if(info == 0) {
2163 ATH_MSG_INFO("SolveWithEigen solving OK");
2164 if(m_logStream)
2165 *m_logStream<<"SolveWithEigen solving OK.\n";
2166 }
2167 else {
2168 ATH_MSG_ERROR( "SolveWithEigen error code (0 if OK) = "<<info );
2169 if(m_logStream)
2170 *m_logStream<<"SolveWithEigen error code (0 if OK) = "<<info<<"\n";
2171 }
2172
2173 if( isCopy )
2174 delete ABetterMat;
2175 ABetterMat = nullptr;
2176
2177 // stop measuring time
2178 clock_t stoptime = clock();
2179 double totaltime = (stoptime-starttime)/double(CLOCKS_PER_SEC);
2180 ATH_MSG_INFO("Time spent in SolveWithEigen: "<<totaltime<<" s");
2181
2182 ATH_MSG_DEBUG("Alignment constants:");
2183 // compute alignment corrections (translations in mm and rotations in rad)
2184 // for solveSparseEigen variances are not calculated
2185 for(int i=0; i<m_aNDoF; i++) {
2186
2187 double param = -(*aBetterVec)[i];
2188 double err = 0.;
2189
2190 int idof = m_activeIndices[i];
2191 AlignPar * alignPar=(*alignParList)[idof];
2192
2193 // undo the sigma scaling
2194 double sigma = alignPar->sigma();
2195 param *= sigma;
2196
2197 ATH_MSG_DEBUG(i <<" : "<< param << " +/- "<< err);
2198 ATH_MSG_DEBUG("sigma: "<<sigma);
2199 ATH_MSG_DEBUG("init par: "<<alignPar->initPar());
2200 alignPar->setPar(param, err);
2201 ATH_MSG_DEBUG("set param to "<<param<<" for alignPar "<<alignPar);
2202 ATH_MSG_DEBUG(*(*alignParList)[idof]);
2203 }
2204
2205 if(m_logStream) {
2206 std::unique_ptr<CLHEP::HepSymMatrix> cov;
2208
2209 // norm of first derivative
2210 *m_logStream<<"norm of first derivative : "<<origVec.norm()<<"\n";
2211
2212 // distance to solution
2213 double dist = ( (*aBetterMat) * (*aBetterVec) - origVec ).norm();
2214 *m_logStream<<"distance to solution : "<<dist<<"\n";
2215
2216 // calculate chi2 of the alignment change
2217 double chi2 = (*aBetterVec) * (*aBetterMat) * (*aBetterVec) * .5;
2218 *m_logStream<<"delta(chi2) of the alignment change : "<<chi2<<" / "<<m_aNDoF<<"\n";
2219
2220 // time spent here
2221 *m_logStream<<"time spent in solve : "<<totaltime<<" s\n";
2222 }
2223
2224 delete aBetterMat;
2225 delete aBetterVec;
2226
2227 // need to do this since success return value from Lapack is 0
2228 // and from SolveWithEigen() it is 1
2229 if (info==0)
2230 info = 1;
2231
2232 return info;
2233 }

◆ spuriousRemoval()

StatusCode MatrixTool::spuriousRemoval ( )
private

Definition at line 1681 of file MatrixTool.cxx.

1682 {
1683
1684 // copied from SiGlobalChi2Algs
1685 ATH_MSG_DEBUG("in spuriousRemoval");
1686
1687 // compute determinant before resizing
1688 if (m_calDet) {
1689 const double tol = 1.e-20;
1690 double determ = m_bigmatrix->determinant();
1691 ATH_MSG_INFO("Determinant: " << determ);
1692 if (std::fabs(determ) < tol)
1693 ATH_MSG_WARNING("Matrix is singular!");
1694 }
1695
1696 // fill vector with modules that need to be removed
1697 int fillvecmods=fillVecMods();
1698 if (fillvecmods==0) {
1699
1700 //ATH_MSG_INFO(" No resize needed (NhitsCut = "
1701 // << m_hitscut << ")");
1702
1703 if (msgLvl(MSG::DEBUG)) {
1704 //m_bigmatrix->Write("bigmatrix.txt", false, m_wSqMatrix, m_scale,
1705 // MatVersion);
1706 //m_bigvector->Write("bigvector.txt", false, m_scale, m_modcodemap,
1707 // VecVersion, m_alignProcessLevel, m_fitParam);
1708 }
1709
1710 return StatusCode::SUCCESS;
1711 }
1712 else if (fillvecmods==2)
1713 return StatusCode::FAILURE;
1714
1715 /* this is a bit difficult to implement for now....
1716
1717 // remove matrix/vector elements
1718 int cont=0;
1719 ModuleIndexMap::iterator itcode;
1720 ATH_MSG_INFO("Eliminating module...");
1721
1722
1723
1724 for (std::vector<int>::const_iterator it=m_dropmods.begin();
1725 it!= m_dropmods.end(); ++it) {
1726
1727 itcode = m_modcodemap.find((*it));
1728
1729 if (itcode == m_modcodemap.end()) {
1730 ATH_MSG_WARNING("Could not find module " << *it << " in map.");
1731 return StatusCode::FAILURE;
1732 }
1733
1734 ATH_MSG_INFO(" - Removing mcode: " << (itcode->second)
1735 << " (old index: " << (*it) << " -> new index: " << (*it)-cont << ")");
1736
1737 m_bigmatrix->RemoveModule((*it)-cont);
1738 m_bigvector->RemoveModule((*it)-cont);
1739 cont++;
1740 }
1741 ATH_MSG_INFO("Modules removed from the Matrix and from the Vector Successfully!");
1742 ATH_MSG_DEBUG("(DoF: " << nDoF << ")");
1743
1744 // resizing...
1745 ATH_MSG_INFO("Resizing the bigvector in memory...");
1746 m_bigvector->reSize(nDoF-6*m_dropmods.size());
1747 ATH_MSG_INFO("Resizing the bigmatrix in memory...");
1748 m_bigmatrix->reSize(nDoF-6*m_dropmods.size());
1749
1750 ATH_MSG_INFO(m_dropmods.size() << " modules eliminated from the matrix, i.e, "
1751 << 6*m_dropmods.size() << " DoFs");
1752 ATH_MSG_INFO(nDoF/6 - m_dropmods.size() << " modules to align (" << nDoF-6*m_dropmods.size() << " DoFs)");
1753 ATH_MSG_INFO("New bigmatrix size is: " << m_bigmatrix->size());
1754 ATH_MSG_INFO("New bigvector size is: " << (*m_bigvector).size());
1755
1756 NDoF = nDoF-6*(m_dropmods.size());
1757
1758 // resize (update) nDoF variable
1759 nDoF = NDoF;
1760
1761 // Resizing vectors to store results
1762 m_alignPar->reSize(nDoF);
1763 m_alignSqErr->reSize(nDoF);
1764
1765 // Fill the mapping module map and updating m_modcodemap
1766 UpdateModMap();
1767 */
1768
1769 return StatusCode::SUCCESS;
1770 }
static int fillVecMods()

◆ storeInTFile()

void MatrixTool::storeInTFile ( const TString & filename)

Store Files in a tfile.

Definition at line 725 of file MatrixTool.cxx.

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

◆ 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 2236 of file MatrixTool.cxx.

2237 {
2238 ATH_MSG_INFO("writing the hitmap to file");
2239
2240 const AlignModuleList * moduleList = m_alignModuleTool->alignModules1D();
2241 int nModules = moduleList->size();
2242
2243 AlMat hitmap(nModules,2);
2244 AlignModuleList::const_iterator imod = moduleList->begin();
2245 AlignModuleList::const_iterator imod_end = moduleList->end();
2246 int index(0);
2247 for(; imod != imod_end; ++imod) {
2248 AlignModule * module = *imod;
2249 hitmap[index][0] = module->nHits();
2250 hitmap[index][1] = module->nTracks();
2251 index++;
2252 }
2253
2254 // Set Path for the hitmap matrix
2255 hitmap.SetPathBin(m_pathbin.value()+m_prefixName.value());
2256 hitmap.SetPathTxt(m_pathtxt.value()+m_prefixName.value());
2257
2258 StatusCode sc = hitmap.Write("hitmap.bin",true); // write the hitmap matrix
2259
2260 if (sc!=StatusCode::SUCCESS)
2261 ATH_MSG_ERROR("Problem writing hitmap matrix");
2262
2263 if (m_writeHitmapTxt) {
2264 sc = hitmap.Write("hitmap.txt", false, 0);
2265 if (sc!=StatusCode::SUCCESS)
2266 ATH_MSG_ERROR("Problem writing hitmap matrix to text file");
2267 }
2268
2269 ATH_MSG_DEBUG("hitmap written to: "<< m_pathbin.value()+m_prefixName.value() <<"hitmap.bin (.txt)");
2270 }
Gaudi::Property< bool > m_writeHitmapTxt
Definition MatrixTool.h:192

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 249 of file MatrixTool.h.

249{};

◆ 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 172 of file MatrixTool.h.

172 {this, "AlignCorrDBStep", 10.,
173 "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 256 of file MatrixTool.h.

257{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 143 of file MatrixTool.h.

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

◆ m_AlignPixelbutNotIBL

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

Definition at line 258 of file MatrixTool.h.

259{this, "AlignPixelbutNotIBL", false};

◆ m_aNDoF

int Trk::MatrixTool::m_aNDoF = 0
private

number of active DoF (size of m_activeIndices)

Definition at line 250 of file MatrixTool.h.

◆ m_bigmatrix

std::unique_ptr<AlSymMatBase> Trk::MatrixTool::m_bigmatrix
private

matrix to contain second derivative terms to be used for alignment

Definition at line 147 of file MatrixTool.h.

◆ m_bigvector

std::unique_ptr<AlVec> Trk::MatrixTool::m_bigvector
private

vector to contain first derivative terms to be used for alignment

Definition at line 150 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 215 of file MatrixTool.h.

216 {this, "CalculateFullCovariance", true,
217 "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 175 of file MatrixTool.h.

175 {this, "MatrixDet", false,
176 "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 261 of file MatrixTool.h.

262{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 155 of file MatrixTool.h.

155 {this, "Diagonalize", true,
156 "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 170 of file MatrixTool.h.

170 {this, "EigenvalueStep", 1e3,
171 "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 157 of file MatrixTool.h.

158{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 241 of file MatrixTool.h.

242 {this, "InputHitmapFiles", {"hitmap.bin"},
243 "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 234 of file MatrixTool.h.

235 {this, "InputMatrixFiles", {"matrix.bin"},
236 "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 245 of file MatrixTool.h.

246 {this, "InputTFiles", {"AlignmentTFile.root"},
247 "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 237 of file MatrixTool.h.

238 {this, "InputVectorFiles", {"vector.bin"},
239 "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 252 of file MatrixTool.h.

252 {this, "MaxReadErrors", 10,
253 "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 164 of file MatrixTool.h.

164 {this, "MinNumHitsPerModule", 0,
165 "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 166 of file MatrixTool.h.

166 {this, "MinNumTrksPerModule", 0,
167 "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 162 of file MatrixTool.h.

162 {this, "ModCut", 0,
163 "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 219 of file MatrixTool.h.

219 {this, "PathBinName", "./",
220 "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 221 of file MatrixTool.h.

221 {this, "PathTxtName", "./",
222 "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 223 of file MatrixTool.h.

223 {this, "PrefixName", "",
224 "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 168 of file MatrixTool.h.

168 {this, "PullCut", 1.0,
169 "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 194 of file MatrixTool.h.

194 {this, "ReadHitmaps", false,
195 "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 199 of file MatrixTool.h.

199 {this, "ReadTFile", false,
200 "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 276 of file MatrixTool.h.

276{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 277 of file MatrixTool.h.

277{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 278 of file MatrixTool.h.

278{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 273 of file MatrixTool.h.

273{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 274 of file MatrixTool.h.

274{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 275 of file MatrixTool.h.

275{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 268 of file MatrixTool.h.

268{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 269 of file MatrixTool.h.

269{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 270 of file MatrixTool.h.

270{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 265 of file MatrixTool.h.

265{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 266 of file MatrixTool.h.

266{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 267 of file MatrixTool.h.

267{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 212 of file MatrixTool.h.

212 {this, "RemoveSpurious", false,
213 "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 202 of file MatrixTool.h.

202 {this, "RunLocalMethod", true,
203 "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 229 of file MatrixTool.h.

230{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 231 of file MatrixTool.h.

232{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 205 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 206 of file MatrixTool.h.

206 {this, "ScaleMatrix", false,
207 "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 209 of file MatrixTool.h.

209 {this, "SoftEigenmodeCut", 0.,
210 "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 160 of file MatrixTool.h.

161{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 226 of file MatrixTool.h.

227{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 153 of file MatrixTool.h.

153{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 183 of file MatrixTool.h.

183 {this, "WriteEigenMat", true,
184 "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 185 of file MatrixTool.h.

185 {this, "WriteEigenMatTxt", true,
186 "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 190 of file MatrixTool.h.

190 {this, "WriteHitmap", false,
191 "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 192 of file MatrixTool.h.

192 {this, "WriteHitmapTxt", false,
193 "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 179 of file MatrixTool.h.

179 {this, "WriteMat", true,
180 "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 181 of file MatrixTool.h.

181 {this, "WriteMatTxt", true,
182 "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 187 of file MatrixTool.h.

187 {this, "WriteModuleNames", false,
188 "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 197 of file MatrixTool.h.

197 {this, "WriteTFile", false,
198 "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 177 of file MatrixTool.h.

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

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