ATLAS Offline Software
Public Types | Public Member Functions | Static Public Member Functions | Protected Member Functions | Protected Attributes | Private Types | Private Member Functions | Static Private Member Functions | Private Attributes | List of all members
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. More...
 
virtual ~MatrixTool ()
 Virtual destructor. More...
 
StatusCode initialize ()
 initialize More...
 
StatusCode finalize ()
 initialize More...
 
StatusCode allocateMatrix (int nDoF=0)
 allocates memory for big matrix and big vector More...
 
void prepareBinaryFiles (int solveOption)
 reads/writes matrix entries from/to binary files as necessary More...
 
void addFirstDerivatives (AlVec *vector)
 adds first derivative to vector More...
 
void addFirstDerivatives (std::list< int, double > &derivatives)
 adds first derivative to vector for only some entries More...
 
void addFirstDerivative (int irow, double firstderiv)
 
void addSecondDerivatives (AlSymMatBase *matrix)
 adds second derivatives to matrix More...
 
void addSecondDerivatives (std::list< std::pair< int, int >, double > &derivatives)
 adds first derivative to vector for only some entries More...
 
void addSecondDerivative (int irow, int icol, double secondderiv)
 
bool accumulateFromFiles ()
 accumulates derivates from files. More...
 
bool accumulateFromBinaries ()
 accumulates derivates from binary files More...
 
int solve ()
 solves for alignment parameters More...
 
void storeInTFile (const TString &filename)
 Store Files in a tfile. More...
 
bool accumulateFromTFiles ()
 Store Files in a tfile. More...
 
void printModuleSolution (std::ostream &os, const AlignModule *module, const CLHEP::HepSymMatrix *cov) const
 namespace { class RestoreIOSFlags { public: RestoreIOSFlags (std::ostream &os) : m_os(&os), m_flags(m_os->flags()), m_precision(m_os->precision()) {} ~RestoreIOSFlags() { m_os->flags(m_flags); m_os->precision(m_precision); } private: std::ostream *m_os; std::ios_base::fmtflags m_flags; std::streamsize m_precision; }; } More...
 
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. More...
 
const ServiceHandle< StoreGateSvc > & evtStore () const
 The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc. More...
 
const ServiceHandle< StoreGateSvc > & detStore () const
 The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc. More...
 
virtual StatusCode sysInitialize () override
 Perform system initialization for an algorithm. More...
 
virtual StatusCode sysStart () override
 Handle START transition. More...
 
virtual std::vector< Gaudi::DataHandle * > inputHandles () const override
 Return this algorithm's input handles. More...
 
virtual std::vector< Gaudi::DataHandle * > outputHandles () const override
 Return this algorithm's output handles. More...
 
Gaudi::Details::PropertyBase & declareProperty (Gaudi::Property< T > &t)
 
Gaudi::Details::PropertyBase * declareProperty (const std::string &name, SG::VarHandleKey &hndl, const std::string &doc, const SG::VarHandleKeyType &)
 Declare a new Gaudi property. More...
 
Gaudi::Details::PropertyBase * declareProperty (const std::string &name, SG::VarHandleBase &hndl, const std::string &doc, const SG::VarHandleType &)
 Declare a new Gaudi property. More...
 
Gaudi::Details::PropertyBase * declareProperty (const std::string &name, SG::VarHandleKeyArray &hndArr, const std::string &doc, const SG::VarHandleKeyArrayType &)
 
Gaudi::Details::PropertyBase * declareProperty (const std::string &name, T &property, const std::string &doc, const SG::NotHandleType &)
 Declare a new Gaudi property. More...
 
Gaudi::Details::PropertyBase * declareProperty (const std::string &name, T &property, const std::string &doc="none")
 Declare a new Gaudi property. More...
 
void updateVHKA (Gaudi::Details::PropertyBase &)
 
MsgStream & msg () const
 
MsgStream & msg (const MSG::Level lvl) 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 More...
 
int nHits () const
 
void setNTracks (int n)
 set number of tracks More...
 
int nTracks () const
 
void setNMeasurements (int n)
 set number of measurements More...
 
int nMeasurements () const
 
virtual void setLogStream (std::ostream *os)
 sets the output stream for the logfile More...
 

Static Public Member Functions

static const InterfaceID & interfaceID ()
 Retrieve interface ID. More...
 

Protected Member Functions

void renounceArray (SG::VarHandleKeyArray &handlesArray)
 remove all handles from I/O resolution More...
 
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. More...
 

Protected Attributes

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

Private Types

typedef ServiceHandle< StoreGateSvc > StoreGateSvc_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 > &hndl, const SG::VarHandleKeyType &)
 specialization for handling Gaudi::Property<SG::VarHandleKey> More...
 
Gaudi::Details::PropertyBase & declareGaudiProperty (Gaudi::Property< T > &hndl, const SG::VarHandleKeyArrayType &)
 specialization for handling Gaudi::Property<SG::VarHandleKeyArray> More...
 
Gaudi::Details::PropertyBase & declareGaudiProperty (Gaudi::Property< T > &hndl, const SG::VarHandleType &)
 specialization for handling Gaudi::Property<SG::VarHandleBase> More...
 
Gaudi::Details::PropertyBase & declareGaudiProperty (Gaudi::Property< T > &t, const SG::NotHandleType &)
 specialization for handling everything that's not a Gaudi::Property<SG::VarHandleKey> or a <SG::VarHandleKeyArray> More...
 

Static Private Member Functions

static int fillVecMods ()
 

Private Attributes

ToolHandle< IAlignModuleTool > m_alignModuleTool
 
AlSymMatBase * m_bigmatrix
 matrix to contain second derivative terms to be used for alignment More...
 
AlVec * m_bigvector
 vector to contain first derivative terms to be used for alignment More...
 
bool m_useSparse
 flag to use AlSpaMat for the big matrix (default is AlSymMat) More...
 
bool m_diagonalize
 run diagonalization instead of inversion More...
 
double m_eigenvaluethreshold
 cut on the minimum eigenvalue More...
 
int m_solveOption
 solving option More...
 
int m_modcut
 cut on the weak modes which number is <par_modcut More...
 
int m_minNumHits
 cut on the minimum number of hits per module More...
 
int m_minNumTrks
 cut on the minimum number of tracks per module More...
 
float m_pullcut
 pull cut for the automatic weak mode removal method More...
 
float m_eigenvalueStep
 eigenvalue step for the second pass in the automatic weak mode removal method More...
 
float m_Align_db_step
 corr in the diagonal basis step for the third pass in the auto weak mode removal method More...
 
bool m_calDet
 Compute bigmatrix's determinant ? More...
 
bool m_wSqMatrix
 write a triangular matrix by default (true: square format) ? More...
 
bool m_writeMat
 write big matrix and vector into files ? More...
 
bool m_writeMatTxt
 also write big matrix and vector into txt files ? More...
 
bool m_writeEigenMat
 write eigenvalues and eigenvectors into files ? More...
 
bool m_writeEigenMatTxt
 also write eigenvalues and eigenvectors into txt files ? More...
 
bool m_writeModuleNames
 write module name instead of Identifier to vector file More...
 
bool m_writeHitmap
 write hitmap into file More...
 
bool m_writeHitmapTxt
 write hitmap into text file More...
 
bool m_readHitmaps
 accumulate hitymap from files More...
 
bool m_writeTFile
 write out files to a root file More...
 
bool m_readTFiles
 write out files to a root file More...
 
bool m_runLocal
 Run solving using Local method. More...
 
double m_scale
 scale for big matrix and vector normalization More...
 
bool m_scaleMatrix
 scale matrix by number of hits before solving More...
 
double m_softEigenmodeCut
 add constant to diagonal to effectively cut on weak eigenmodes More...
 
bool m_removeSpurious
 run spurious removal More...
 
bool m_calculateFullCovariance
 
std::string m_pathbin
 path binary files (in/out) More...
 
std::string m_pathtxt
 path ascii files (in/out) More...
 
std::string m_prefixName
 prefix string to filenames More...
 
std::string m_tfileName
 prefix string to filenames More...
 
std::string m_scalaMatName
 Scalapack matrix name. More...
 
std::string m_scalaVecName
 Scalapack vector name. More...
 
std::vector< std::string > m_inputMatrixFiles
 input binary files containing matrix terms More...
 
std::vector< std::string > m_inputVectorFiles
 input binary files containing vector terms More...
 
std::vector< std::string > m_inputHitmapFiles
 input binary files containing the hitmaps More...
 
std::vector< std::string > m_inputTFiles
 input binary files containing matrix terms More...
 
std::vector< int > m_activeIndices
 vector of indices which pass the min-hits cut More...
 
int m_aNDoF
 number of active DoF (size of m_activeIndices) More...
 
int m_maxReadErrors
 maximum number of reading TFile errors More...
 
bool m_AlignIBLbutNotPixel
 
bool m_AlignPixelbutNotIBL
 
bool m_DeactivateSCT_ECA_LastDisk
 
bool m_Remove_Pixel_Tx
 
bool m_Remove_Pixel_Ty
 
bool m_Remove_Pixel_Tz
 
bool m_Remove_Pixel_Rx
 
bool m_Remove_Pixel_Ry
 
bool m_Remove_Pixel_Rz
 
bool m_Remove_IBL_Tx
 
bool m_Remove_IBL_Ty
 
bool m_Remove_IBL_Tz
 
bool m_Remove_IBL_Rx
 
bool m_Remove_IBL_Ry
 
bool m_Remove_IBL_Rz
 
StoreGateSvc_t m_evtStore
 Pointer to StoreGate (event store by default) More...
 
StoreGateSvc_t m_detStore
 Pointer to StoreGate (detector store by default) More...
 
std::vector< SG::VarHandleKeyArray * > m_vhka
 
bool m_varHandleArraysDeclared
 
std::map< int, int > m_alignModuleMap
 
int m_nentries
 

Detailed Description

Definition at line 55 of file MatrixTool.h.

Member Typedef Documentation

◆ StoreGateSvc_t

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

Definition at line 388 of file AthCommonDataStore.h.

Member Enumeration Documentation

◆ SolveOption

enum Trk::MatrixTool::SolveOption
Enumerator
NONE 

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

SOLVE 

solving after data accumulation (LAPACK)

SOLVE_FAST 

Fast (Eigen method) solving after data accumulation.

DIRECT_SOLVE 

direct solving (LAPACK), already available matrix & vector

DIRECT_SOLVE_FAST 

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

DIRECT_SOLVE_CLUSTER 

computation of alignment parameters from SCALAPAK already solved matrix

SOLVE_ROOT 

computation using ROOT

SOLVE_CLHEP 

computation using CLHEP

Definition at line 58 of file MatrixTool.h.

58  {
59  NONE = 0,
60  SOLVE = 1,
61  SOLVE_FAST = 2,
62  DIRECT_SOLVE = 3,
63  DIRECT_SOLVE_FAST = 4,
64  DIRECT_SOLVE_CLUSTER = 5,
65  SOLVE_ROOT = 6,
66  SOLVE_CLHEP = 7
67  };

Constructor & Destructor Documentation

◆ MatrixTool()

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

Constructor.

Definition at line 58 of file MatrixTool.cxx.

60  : IMatrixTool()
61  , AthAlgTool(type,name,parent)
62  , m_alignModuleTool("Trk::AlignModuleTool/AlignModuleTool")
63  , m_bigmatrix(nullptr)
64  , m_bigvector(nullptr)
65  //m_useSparse
66  //m_diagonalize
67  //m_eigenvaluethreshold
68  //m_solveOption
69  //m_modcut
70  //m_minNumHits
71  //m_minNumTrks
72  //m_pullcut
73  //m_eigenvalueStep
74  //m_Align_db_step
75  //m_calDet
76  //m_wSqMatrix
77  //m_writeMat
78  //m_writeMatTxt
79  //m_writeEigenMat
80  //m_writeEigenMatTxt
81  //m_writeModuleNames
82  //m_writeHitmap
83  //m_writeHitmapTxt
84  //m_readHitmaps
85  //m_writeTFile
86  //m_readTFiles
87  //m_runLocal
88  , m_scale(-1.)
89  //m_scaleMatrix
90  //m_softEigenmodeCut
91  //m_removeSpurious
92  //m_calculateFullCovariance
93  //m_pathbin, m_pathtxt, m_prefixName, m_tfileName, m_scalaMatName, m_scalaVecName
94  //m_inputMatrixFiles, m_inputVectorFiles, m_inputHitmapFiles, m_inputTFiles
95  //m_activeIndices
96  , m_aNDoF(0)
97  //m_maxReadErrors
98  //m_AlignIBLbutNotPixel, m_AlignPixelbutNotIBL
99  //m_Remove_Pixel_Tx, m_Remove_Pixel_Ty, m_Remove_Pixel_Tz
100  //m_Remove_Pixel_Rx, m_Remove_Pixel_Ry, m_Remove_Pixel_Rz
101  //m_Remove_IBL_Tx, m_Remove_IBL_Ty, m_Remove_IBL_Tz
102  //m_Remove_IBL_Rx, m_Remove_IBL_Ry, m_Remove_IBL_Rz
103  {
104  declareInterface<IMatrixTool>(this);
105 
106  declareProperty("UseSparse", m_useSparse = false);
107  declareProperty("SolveOption", m_solveOption = NONE);
108 
109  declareProperty("Diagonalize", m_diagonalize = true);
110  declareProperty("EigenvalueThreshold", m_eigenvaluethreshold = 0.);
111 
112  declareProperty("ModCut", m_modcut = 0);
113  declareProperty("PullCut", m_pullcut = 1.0);
114  declareProperty("EigenvalueStep", m_eigenvalueStep = 1e3);
115  declareProperty("AlignCorrDBStep", m_Align_db_step = 10);
116  declareProperty("MinNumHitsPerModule", m_minNumHits = 0);
117  declareProperty("MinNumTrksPerModule", m_minNumTrks = 0);
118  declareProperty("RunLocalMethod", m_runLocal = true);
119 
120  declareProperty("MatrixDet", m_calDet = false);
121  declareProperty("WriteSquareMatrix", m_wSqMatrix = false);
122  declareProperty("WriteMat", m_writeMat = true);
123  declareProperty("WriteMatTxt", m_writeMatTxt = true);
124  declareProperty("WriteEigenMat", m_writeEigenMat = true);
125  declareProperty("WriteEigenMatTxt", m_writeEigenMatTxt = true);
126  declareProperty("WriteModuleNames", m_writeModuleNames = false);
127 
128  declareProperty("WriteTFile", m_writeTFile = false);
129  // if True then files will be read from TFiles instead of Binary files
130  declareProperty("ReadTFile", m_readTFiles = false);
131 
132 
133  std::vector<std::string> defaultMatInput,defaultVecInput,defaultTFile;
134  defaultMatInput.emplace_back("matrix.bin");
135  defaultVecInput.emplace_back("vector.bin");
136  defaultTFile.emplace_back("AlignmentTFile.root");
137 
138  declareProperty("InputMatrixFiles", m_inputMatrixFiles = defaultMatInput);
139  declareProperty("InputVectorFiles", m_inputVectorFiles = defaultVecInput);
140  declareProperty("InputTFiles", m_inputTFiles = defaultTFile);
141 
142  declareProperty("AlignModuleTool", m_alignModuleTool);
143 
144  declareProperty("PathBinName", m_pathbin = "./");
145  declareProperty("PathTxtName", m_pathtxt = "./");
146  declareProperty("PrefixName", m_prefixName = "");
147 
148  declareProperty("TFileName", m_tfileName = "AlignmentTFile.root");
149 
150  declareProperty("SoftEigenmodeCut", m_softEigenmodeCut = 0.);
151  declareProperty("RemoveSpurious", m_removeSpurious = false);
152  declareProperty("CalculateFullCovariance", m_calculateFullCovariance = true);
153 
154  // ScaLapack
155  declareProperty("ScalapackMatrixName", m_scalaMatName = "eigenvectors.bin");
156  declareProperty("ScalapackVectorName", m_scalaVecName = "eigenvalues.bin");
157 
158  declareProperty("ScaleMatrix", m_scaleMatrix = false);
159 
160  // Hitmap
161  declareProperty("WriteHitmap", m_writeHitmap = false);
162  declareProperty("WriteHitmapTxt", m_writeHitmapTxt = false);
163  declareProperty("ReadHitmaps", m_readHitmaps = false);
164  std::vector<std::string> defaultHitmapInput;
165  defaultMatInput.emplace_back("hitmap.bin");
166  declareProperty("InputHitmapFiles", m_inputHitmapFiles = defaultHitmapInput);
167 
168  declareProperty("MaxReadErrors", m_maxReadErrors = 10);
169  //To skip IBL or Pixel Alignment
170  declareProperty("AlignIBLbutNotPixel", m_AlignIBLbutNotPixel = false);
171  declareProperty("AlignPixelbutNotIBL", m_AlignPixelbutNotIBL = false);
172 
173  //To Skip Solving of SCT ECA Last Disk
174  declareProperty("DeactivateSCT_ECA_LastDisk", m_DeactivateSCT_ECA_LastDisk = false);
175 
176  //By Pixel DoF
177  declareProperty("Remove_Pixel_Tx", m_Remove_Pixel_Tx = false);
178  declareProperty("Remove_Pixel_Ty", m_Remove_Pixel_Ty = false);
179  declareProperty("Remove_Pixel_Tz", m_Remove_Pixel_Tz = false);
180  declareProperty("Remove_Pixel_Rx", m_Remove_Pixel_Rx = false);
181  declareProperty("Remove_Pixel_Ry", m_Remove_Pixel_Ry = false);
182  declareProperty("Remove_Pixel_Rz", m_Remove_Pixel_Rz = false);
183 
184  //By IBL DoF
185  declareProperty("Remove_IBL_Tx", m_Remove_IBL_Tx = false);
186  declareProperty("Remove_IBL_Ty", m_Remove_IBL_Ty = false);
187  declareProperty("Remove_IBL_Tz", m_Remove_IBL_Tz = false);
188  declareProperty("Remove_IBL_Rx", m_Remove_IBL_Rx = false);
189  declareProperty("Remove_IBL_Ry", m_Remove_IBL_Ry = false);
190  declareProperty("Remove_IBL_Rz", m_Remove_IBL_Rz = false);
191 
192  }

◆ ~MatrixTool()

MatrixTool::~MatrixTool ( )
virtual

Virtual destructor.

Definition at line 195 of file MatrixTool.cxx.

196  {
197  delete m_bigmatrix;
198  delete m_bigvector;
199  }

Member Function Documentation

◆ accumulateFromBinaries()

bool MatrixTool::accumulateFromBinaries ( )

accumulates derivates from binary files

Definition at line 763 of file MatrixTool.cxx.

764  {
765 
766 
767  DataVector<AlignPar>* alignParList = m_alignModuleTool->alignParList1D();
768  int nDoF=alignParList->size();
769 
770  std::map<int,unsigned long long> modIndexMap;
771  float dummyVersion(0.);
772  double totalscale=0.;
773  for (int ivec=0;ivec<(int)m_inputVectorFiles.size();ivec++) {
774 
775  ATH_MSG_DEBUG("Reading vector "<<ivec<<" from file "<<m_inputVectorFiles[ivec]);
776 
777  AlVec newVector(nDoF);
778  std::map<int,unsigned long long> newModIndexMap;
779  newVector.SetPathBin(m_pathbin+m_prefixName);
780  newVector.SetPathTxt(m_pathtxt+m_prefixName);
781  double scale=0;
782  StatusCode sc = newVector.ReadPartial(m_inputVectorFiles[ivec],scale,newModIndexMap,dummyVersion);
783  totalscale += scale;
784  if (sc==StatusCode::FAILURE) {
785  msg(MSG::FATAL)<<"Problem reading vector from "<<m_inputVectorFiles[ivec]<<endmsg;
786  return false;
787  }
788  if (newVector.size()!=m_bigvector->size()) {
789  msg(MSG::FATAL) <<"vector wrong size! newVector size "<<newVector.size()
790  <<", bigvector size "<<m_bigvector->size()<<endmsg;
791  return false;
792  }
793 
794  // check modIndexMaps to make sure they are the same
795  if (ivec==0)
796  modIndexMap = newModIndexMap;
797  else if (modIndexMap!=newModIndexMap) {
798  msg(MSG::FATAL)<<"module index maps don't agree!"<<endmsg;
799  return false;
800  }
801  if (ivec>0)
802  *m_bigvector += newVector;
803  else
804  *m_bigvector = newVector;
805  }
806 
807  m_scale = totalscale;
808 
809  AlSymMat * symBigMatrix=dynamic_cast<AlSymMat*>(m_bigmatrix);
810  AlSpaMat * spaBigMatrix=dynamic_cast<AlSpaMat*>(m_bigmatrix);
811 
812 
813  for (int imat=0;imat<(int)m_inputMatrixFiles.size();imat++) {
814  ATH_MSG_DEBUG("Reading matrix "<<imat<<" from file "<<m_inputMatrixFiles[imat]);
815 
816  // create new matrix to read data from current file
817  int nDoF=modIndexMap.size();
818  bool triang;
819  StatusCode sc;
820  if (symBigMatrix) {
821  AlSymMat newMatrix(nDoF);
822  sc = newMatrix.Read(m_inputMatrixFiles[imat],nDoF,triang,dummyVersion);
823  if (sc==StatusCode::SUCCESS)
824  *symBigMatrix += newMatrix;
825  }
826  else {
827  if (!spaBigMatrix) {
828  throw std::logic_error("Unhandled matrix type");
829  }
830 
831  AlSpaMat newMatrix(nDoF);
832  sc = newMatrix.Read(m_inputMatrixFiles[imat],nDoF,triang,dummyVersion);
833 
834  if (sc==StatusCode::SUCCESS) {
835  if (imat>0)
836  *spaBigMatrix += newMatrix;
837  else
838  *spaBigMatrix = newMatrix;
839  }
840  }
841 
842  if (sc==StatusCode::FAILURE) {
843  msg(MSG::FATAL)<<"problem reading matrix from "<<m_inputMatrixFiles[imat]<<endmsg;
844  return false;
845  }
846 
847  if (!m_useSparse && triang==m_wSqMatrix) {
848  ATH_MSG_WARNING("matrix not expected format! Changing m_wSqMatrix to "<<!triang);
849  m_wSqMatrix=!triang;
850  }
851 
852  }
853 
854  // accumulate hitmap from hitmap files
855  if(m_readHitmaps)
856  readHitmaps();
857 
858  return true;
859  }

◆ accumulateFromFiles()

bool MatrixTool::accumulateFromFiles ( )
virtual

accumulates derivates from files.

Flag decides if it is binary or TFiles

Implements Trk::IMatrixTool.

Definition at line 750 of file MatrixTool.cxx.

751  {
752 
753  if(m_readTFiles){
754  ATH_MSG_INFO("Info to obtained from from TFiles");
755  return accumulateFromTFiles();
756  }else{
757  ATH_MSG_INFO("Info to obtained from from Binary files");
758  return accumulateFromBinaries();
759  }
760  }

◆ accumulateFromTFiles()

bool MatrixTool::accumulateFromTFiles ( )

Store Files in a tfile.

Definition at line 957 of file MatrixTool.cxx.

958  {
959  DataVector<AlignPar>* alignParList = m_alignModuleTool->alignParList1D();
960  int nDoF=alignParList->size();
961  ATH_MSG_DEBUG("OPENING TFILES");
962 
963  std::map<int,unsigned long long> modIndexMap;
964  std::map<int,unsigned long long> DoFMap;
965  double totalscale=0.;
966 
967  AlSymMat * symBigMatrix=dynamic_cast<AlSymMat*>(m_bigmatrix);
968  AlSpaMat * spaBigMatrix=dynamic_cast<AlSpaMat*>(m_bigmatrix);
969  //TMatrixDSparse *accumMatrix(0);
970  AlSpaMat *accumMatrix = nullptr;
971 
972  const AlignModuleList * moduleList = m_alignModuleTool->alignModules1D();
973  int nModules = moduleList->size();
974 
975  TVectorD TotalHits(nModules);
976  TVectorD TotalTracks(nModules);
977 
978  int numberOfReadErrors = 0;
979 
980  struct rusage myusage{};
981  int itworked = getrusage(RUSAGE_SELF,&myusage);
982  if(itworked)//note: rusage returns zero if it succeeds!
983  ATH_MSG_DEBUG("ItWorked");
984 
985  long intialMemUse = myusage.ru_maxrss;
986 
987  for (int ifile = 0; ifile < (int)m_inputTFiles.size(); ifile++) {
988  if (numberOfReadErrors > m_maxReadErrors){
989  msg(MSG::FATAL) << " number of errors when reading the TFiles already exceed " << m_maxReadErrors << endmsg;
990  return false;
991  }
992 
993  ATH_MSG_DEBUG("Reading File number " << ifile << ", " << m_inputTFiles[ifile]);
994 
995  itworked = getrusage(RUSAGE_SELF,&myusage);
996  ATH_MSG_DEBUG("Memory usage [MB], total " << myusage.ru_maxrss/1024 << ", increase " << (myusage.ru_maxrss-intialMemUse)/1024);
997 
998  TFile* myFile = TFile::Open(m_inputTFiles[ifile].c_str());
999 
1000  if ( myFile->IsZombie() || !(myFile->IsOpen()) ) {
1001  ++numberOfReadErrors;
1002  ATH_MSG_ERROR( " Problem reading TFile " << m_inputTFiles[ifile] );
1003  continue;
1004  }
1005 
1006  std::map<int,unsigned long long> newModIndexMap;
1007 
1008  TVectorD* myModuleIDs;
1009  myModuleIDs = (TVectorD*)myFile->Get("ModuleID");
1010  if( !myModuleIDs ){
1011  ++numberOfReadErrors;
1012  ATH_MSG_ERROR("Modules ID not read!!!");
1013  continue;
1014  }
1015 
1016  for (int i(0); i<myModuleIDs->GetNrows(); ++i){
1017  //Coverting back from a double to a unvi signed long long
1018  double source = (*myModuleIDs)(i);
1019  uint64_t target;
1020  memcpy(&target, &source, sizeof(target));
1021  newModIndexMap[i]=target;
1022  //std::cout << i<< " " <<target <<std::endl;
1023  }
1024 
1025  delete myModuleIDs;
1026 
1027  std::map<int,unsigned long long> newDoFMap;
1028 
1029  TVectorD* myDoFs;
1030  myDoFs = (TVectorD*)myFile->Get("dof");
1031  if( !myDoFs ){
1032  ++numberOfReadErrors;
1033  ATH_MSG_ERROR("DoFs not read!!!");
1034  continue;
1035  }
1036 
1037  for (int i(0); i<myDoFs->GetNrows(); ++i){
1038  //Coverting back from a double to a unsigned long long
1039  double source = (*myDoFs)(i);
1040  uint64_t target;
1041  memcpy(&target, &source, sizeof(target));
1042  newDoFMap[i]=target;
1043  }
1044  delete myDoFs;
1045 
1046 
1047  TVectorD* Scale;
1048  Scale = (TVectorD*)myFile->Get("Scale");
1049  if( !Scale ){
1050  ++numberOfReadErrors;
1051  ATH_MSG_ERROR("Scale not read!!!");
1052  continue;
1053  }
1054 
1055  double scale=(*Scale)(0);
1056  totalscale += scale;
1057  delete Scale;
1058 
1059 
1060  ATH_MSG_DEBUG("Reading Vector");
1061  TVectorD* vector = (TVectorD*)myFile->Get("Vector");
1062  if( !vector ){
1063  ++numberOfReadErrors;
1064  ATH_MSG_ERROR("Vector not read!!!");
1065  continue;
1066  }
1067 
1068  AlVec* newVector = new AlVec(nDoF);
1069  newVector->SetPathBin(m_pathbin+m_prefixName);
1070  newVector->SetPathTxt(m_pathtxt+m_prefixName);
1071 
1072  if (newVector->size() != m_bigvector->size() ) {
1073  msg(MSG::FATAL) << "vector wrong size! newVector size " << newVector->size()
1074  << ", bigvector size " << m_bigvector->size()<<endmsg;
1075  delete newVector;
1076  delete vector;
1077  return false;
1078  }
1079 
1080  if (m_bigvector->size() != vector->GetNrows() ) {
1081  msg(MSG::FATAL) << "File vector wrong size! File Vector size " << vector->GetNrows()
1082  << ", bigvector size " << m_bigvector->size()<<endmsg;
1083  delete newVector;
1084  delete vector;
1085  return false;
1086  }
1087 
1088 
1089  for (int i=0;i<nDoF;i++) {
1090  (*newVector)[i] = (*vector)(i);
1091  }
1092  delete vector;
1093 
1094  // check modIndexMaps to make sure they are the same
1095  if (ifile == 0){
1096  DoFMap = newDoFMap;
1097  } else if (DoFMap!=newDoFMap) {
1098  msg(MSG::FATAL) << "module dofs don't agree!" << endmsg;
1099  return false;
1100  }
1101 
1102  if (ifile == 0){
1103  modIndexMap = newModIndexMap;
1104  } else if (modIndexMap!=newModIndexMap) {
1105  msg(MSG::FATAL) << "module index maps don't agree!" << endmsg;
1106  return false;
1107  }
1108 
1109  if (ifile>0){
1110  *m_bigvector += *newVector;
1111  delete newVector;
1112  } else {
1113  delete m_bigvector;
1114  m_bigvector = newVector;
1115  }
1116 
1117 
1118  ATH_MSG_DEBUG("Reading matrix ");
1119  TMatrixDSparse* matrix = (TMatrixDSparse*)myFile->Get("Matrix");
1120 
1121  if( !matrix ){
1122  ++numberOfReadErrors;
1123  ATH_MSG_ERROR("Matrix not read!!!");
1124  continue;
1125  }
1126 
1127 
1128  if (ifile == 0 ){
1129 
1130  accumMatrix = new AlSpaMat(nDoF);
1131  ATH_MSG_DEBUG("Matrix size b4 "<< accumMatrix->ptrMap()->size() );
1132 
1133  //This method is ok for large matrix files... really only access the non zero elements
1134  for (int ii=0;ii<nDoF;ii++) {
1135  const TMatrixTSparseRow_const<double> myRow = (*matrix)[ii];
1136  int i = myRow.GetRowIndex();
1137  for (int jj=0;jj<myRow.GetNindex();jj++) {
1138  int j = (myRow.GetColPtr())[jj];
1139  const double myElement= (myRow.GetDataPtr())[jj];
1140  if (i<j){
1141  ATH_MSG_DEBUG("i < j " );
1142  j = i;
1143  i = (myRow.GetColPtr())[jj];
1144  }
1145  (*accumMatrix)[i][j] = myElement;
1146  }
1147  }
1148  ATH_MSG_DEBUG("Matrix size AF "<< accumMatrix->ptrMap()->size() );
1149 
1150  } else if ( accumMatrix) {
1151  ATH_MSG_DEBUG("Matrix size b4 "<< accumMatrix->ptrMap()->size() );
1152 
1153  for (int ii=0;ii<nDoF;ii++) {
1154  const TMatrixTSparseRow_const<double> myRow = (*matrix)[ii];
1155  int i = myRow.GetRowIndex();
1156  for (int jj=0;jj<myRow.GetNindex();jj++) {
1157  int j = (myRow.GetColPtr())[jj];
1158  const double myElement= (myRow.GetDataPtr())[jj];
1159  if (i<j){
1160  ATH_MSG_DEBUG("i < j " );
1161  j = i;
1162  i = (myRow.GetColPtr())[jj];
1163  }
1164  (*accumMatrix)[i][j] += myElement;
1165  }
1166  }
1167  ATH_MSG_DEBUG("Matrix size AF "<< accumMatrix->ptrMap()->size() );
1168 
1169  } else {
1170  delete matrix;
1171  ++numberOfReadErrors;
1172  ATH_MSG_ERROR("Matrix allocation error!!!");
1173  continue;
1174  }
1175 
1176  delete matrix;
1177 
1178  TVectorD* hits;
1179  TVectorD* tracks;
1180 
1181  ATH_MSG_DEBUG("Reading hitmap ");
1182  hits = (TVectorD*)myFile->Get("Hits");
1183  if( !hits ){
1184  ++numberOfReadErrors;
1185  ATH_MSG_ERROR("Hitmap 1 not read!!!");
1186  continue;
1187  }
1188 
1189  tracks = (TVectorD*)myFile->Get("Tracks");
1190  if( !tracks ){
1191  delete hits;
1192  ++numberOfReadErrors;
1193  ATH_MSG_ERROR("Hitmap 2 not read!!!");
1194  continue;
1195  }
1196 
1197  if(hits->GetNrows() != TotalHits.GetNrows() ){
1198  delete hits;
1199  delete tracks;
1200  ++numberOfReadErrors;
1201  ATH_MSG_ERROR("Hitmap size incorrect!!!");
1202  continue;
1203  }
1204 
1205  TotalHits += (*hits);
1206  TotalTracks += (*tracks);
1207 
1208  delete hits;
1209  delete tracks;
1210 
1211  myFile->Close("R");
1212  delete myFile;
1213 
1214  itworked = getrusage(RUSAGE_SELF,&myusage);
1215  ATH_MSG_DEBUG("Memory usage [MB], total " << myusage.ru_maxrss/1024 << ", increase " << (myusage.ru_maxrss-intialMemUse)/1024);
1216 
1217  }
1218 
1219 
1220 
1221  // create new matrix to read data from current file
1222  if(accumMatrix){
1223  if (symBigMatrix) {
1224  AlSymMat newMatrix(nDoF);
1225  //This method is ok for small matrix files
1226  for (int i=0;i<nDoF;i++) {
1227  for (int j=0;j<=i;j++) {
1228  newMatrix[i][j] = (*accumMatrix)[i][j];
1229  }
1230  }
1231 
1232  *symBigMatrix += newMatrix;
1233  delete accumMatrix;
1234  } else if (spaBigMatrix) {
1235  ATH_MSG_DEBUG( "should reassign matrix "<< spaBigMatrix->ptrMap()->size() );
1236  *spaBigMatrix += *accumMatrix;
1237  ATH_MSG_DEBUG( "?????? "<< spaBigMatrix->ptrMap()->size() );
1238  delete accumMatrix;
1239  }
1240  }
1241 
1242  ATH_MSG_DEBUG( "?????? "<< m_bigmatrix->ptrMap()->size() );
1243 
1244  AlignModuleList::const_iterator imod = moduleList->begin();
1245  AlignModuleList::const_iterator imod_end = moduleList->end();
1246  int index = 0;
1247  int totalhits = 0;
1248  for(; imod != imod_end; ++imod, ++index ) {
1249  AlignModule * module = *imod;
1250  module->setNHits((int)TotalHits(index));
1251  module->setNTracks((int)TotalTracks(index));
1252  totalhits += (int)TotalHits(index);
1253  }
1254 
1255 
1256  m_nHits = totalhits;
1257  m_nTracks = 0;
1258  m_nMeasurements = 0;
1259  m_scale = totalscale;
1260 
1261  return true;
1262  }

◆ addFirstDerivative()

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

Implements Trk::IMatrixTool.

Definition at line 1609 of file MatrixTool.cxx.

1610  {
1611  (*m_bigvector)[irow] += firstderiv;
1612  }

◆ addFirstDerivatives() [1/2]

void MatrixTool::addFirstDerivatives ( AlVec *  vector)
virtual

adds first derivative to vector

Implements Trk::IMatrixTool.

Definition at line 1589 of file MatrixTool.cxx.

1590  {
1591  }

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

1600  {
1601  }

◆ 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; }

◆ addSecondDerivative()

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

Implements Trk::IMatrixTool.

Definition at line 1615 of file MatrixTool.cxx.

1616  {
1617  (*m_bigmatrix)[irow][icol] += secondderiv;
1618  }

◆ addSecondDerivatives() [1/2]

void MatrixTool::addSecondDerivatives ( AlSymMatBase *  matrix)
virtual

adds second derivatives to matrix

Implements Trk::IMatrixTool.

Definition at line 1594 of file MatrixTool.cxx.

1595  {
1596  }

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

1605  {
1606  }

◆ allocateMatrix()

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

allocates memory for big matrix and big vector

Implements Trk::IMatrixTool.

Definition at line 232 of file MatrixTool.cxx.

233  {
234  ATH_MSG_INFO("allocating matrix and vector with nDoF = "<<nDoF);
235 
236  if (nullptr!=m_bigmatrix || nullptr!=m_bigvector)
237  ATH_MSG_ERROR("big matrix already allocated!");
238 
239  // Decide upon the big matrix representation:
240  if( m_useSparse )
241  m_bigmatrix = new AlSpaMat(nDoF);
242  else
243  m_bigmatrix = new AlSymMat(nDoF);
244 
245  m_bigvector = new AlVec(nDoF);
246 
247  ATH_MSG_INFO(" After Matrix and Vector allocation");
248 
249  // set paths for matrix and vector output
250  m_bigmatrix->SetPathBin(m_pathbin+m_prefixName);
251  m_bigmatrix->SetPathTxt(m_pathtxt+m_prefixName);
252  m_bigvector->SetPathBin(m_pathbin+m_prefixName);
253  m_bigvector->SetPathTxt(m_pathtxt+m_prefixName);
254 
255  ATH_MSG_INFO("set path to "<<m_pathbin+m_prefixName);
256  return StatusCode::SUCCESS;
257  }

◆ declareGaudiProperty() [1/4]

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

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

Definition at line 170 of file AthCommonDataStore.h.

172  {
173  return *AthCommonDataStore<PBASE>::declareProperty(hndl.name(),
174  hndl.value(),
175  hndl.documentation());
176 
177  }

◆ declareGaudiProperty() [2/4]

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

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

Definition at line 156 of file AthCommonDataStore.h.

158  {
159  return *AthCommonDataStore<PBASE>::declareProperty(hndl.name(),
160  hndl.value(),
161  hndl.documentation());
162 
163  }

◆ declareGaudiProperty() [3/4]

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

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

Definition at line 184 of file AthCommonDataStore.h.

186  {
187  return *AthCommonDataStore<PBASE>::declareProperty(hndl.name(),
188  hndl.value(),
189  hndl.documentation());
190  }

◆ declareGaudiProperty() [4/4]

Gaudi::Details::PropertyBase& AthCommonDataStore< AthCommonMsg< AlgTool > >::declareGaudiProperty ( Gaudi::Property< T > &  t,
const SG::NotHandleType &   
)
inlineprivateinherited

specialization for handling everything that's not a Gaudi::Property<SG::VarHandleKey> or a <SG::VarHandleKeyArray>

Definition at line 199 of file AthCommonDataStore.h.

200  {
201  return PBASE::declareProperty(t);
202  }

◆ declareProperty() [1/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( const std::string &  name,
SG::VarHandleBase &  hndl,
const std::string &  doc,
const SG::VarHandleType &   
)
inlineinherited

Declare a new Gaudi property.

Parameters
nameName of the property.
hndlObject holding the property value.
docDocumentation string for the property.

This is the version for types that derive from SG::VarHandleBase. The property value object is put on the input and output lists as appropriate; then we forward to the base class.

Definition at line 245 of file AthCommonDataStore.h.

249  {
250  this->declare(hndl.vhKey());
251  hndl.vhKey().setOwner(this);
252 
253  return PBASE::declareProperty(name,hndl,doc);
254  }

◆ declareProperty() [2/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( const std::string &  name,
SG::VarHandleKey &  hndl,
const std::string &  doc,
const SG::VarHandleKeyType &   
)
inlineinherited

Declare a new Gaudi property.

Parameters
nameName of the property.
hndlObject holding the property value.
docDocumentation string for the property.

This is the version for types that derive from SG::VarHandleKey. The property value object is put on the input and output lists as appropriate; then we forward to the base class.

Definition at line 221 of file AthCommonDataStore.h.

225  {
226  this->declare(hndl);
227  hndl.setOwner(this);
228 
229  return PBASE::declareProperty(name,hndl,doc);
230  }

◆ declareProperty() [3/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( const std::string &  name,
SG::VarHandleKeyArray &  hndArr,
const std::string &  doc,
const SG::VarHandleKeyArrayType &   
)
inlineinherited

Definition at line 259 of file AthCommonDataStore.h.

263  {
264 
265  // std::ostringstream ost;
266  // ost << Algorithm::name() << " VHKA declareProp: " << name
267  // << " size: " << hndArr.keys().size()
268  // << " mode: " << hndArr.mode()
269  // << " vhka size: " << m_vhka.size()
270  // << "\n";
271  // debug() << ost.str() << endmsg;
272 
273  hndArr.setOwner(this);
274  m_vhka.push_back(&hndArr);
275 
276  Gaudi::Details::PropertyBase* p = PBASE::declareProperty(name, hndArr, doc);
277  if (p != 0) {
278  p->declareUpdateHandler(&AthCommonDataStore<PBASE>::updateVHKA, this);
279  } else {
280  ATH_MSG_ERROR("unable to call declareProperty on VarHandleKeyArray "
281  << name);
282  }
283 
284  return p;
285 
286  }

◆ declareProperty() [4/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( const std::string &  name,
T &  property,
const std::string &  doc,
const SG::NotHandleType &   
)
inlineinherited

Declare a new Gaudi property.

Parameters
nameName of the property.
propertyObject holding the property value.
docDocumentation string for the property.

This is the generic version, for types that do not derive from SG::VarHandleKey. It just forwards to the base class version of declareProperty.

Definition at line 333 of file AthCommonDataStore.h.

337  {
338  return PBASE::declareProperty(name, property, doc);
339  }

◆ declareProperty() [5/6]

Gaudi::Details::PropertyBase* AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( const std::string &  name,
T &  property,
const std::string &  doc = "none" 
)
inlineinherited

Declare a new Gaudi property.

Parameters
nameName of the property.
propertyObject holding the property value.
docDocumentation string for the property.

This dispatches to either the generic declareProperty or the one for VarHandle/Key/KeyArray.

Definition at line 352 of file AthCommonDataStore.h.

355  {
356  typedef typename SG::HandleClassifier<T>::type htype;
357  return declareProperty (name, property, doc, htype());
358  }

◆ declareProperty() [6/6]

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

Definition at line 145 of file AthCommonDataStore.h.

145  {
146  typedef typename SG::HandleClassifier<T>::type htype;
147  return AthCommonDataStore<PBASE>::declareGaudiProperty(t, htype());
148  }

◆ 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.

95 { return m_detStore; }

◆ 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() [1/2]

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.

85 { return m_evtStore; }

◆ evtStore() [2/2]

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

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

Definition at line 90 of file AthCommonDataStore.h.

90 { return m_evtStore; }

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

1760  {
1761  return 0;
1762  }

◆ finalize()

StatusCode MatrixTool::finalize ( )

initialize

Definition at line 223 of file MatrixTool.cxx.

224  {
225  ATH_MSG_DEBUG("finalize() of MatrixTool");
226 
227  return StatusCode::SUCCESS;
228  }

◆ initialize()

StatusCode MatrixTool::initialize ( )

initialize

Definition at line 202 of file MatrixTool.cxx.

203  {
204  ATH_MSG_DEBUG("initialize() of MatrixTool");
205 
206  // get AlignModuleTool
207  if (m_alignModuleTool.retrieve().isSuccess())
208  ATH_MSG_INFO("Retrieved " << m_alignModuleTool);
209  else{
210  msg(MSG::FATAL) << "Could not get " << m_alignModuleTool << endmsg;
211  return StatusCode::FAILURE;
212  }
213 
214  ATH_MSG_INFO("Retrieving data from the following files: ");
215  for (auto & inputVectorFile : m_inputVectorFiles) {
216  ATH_MSG_INFO(m_pathbin+inputVectorFile);
217  }
218 
219  return StatusCode::SUCCESS;
220  }

◆ 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  {
111  return IID_TRKALIGNINTERFACES_IMatrixTool;
112  }

◆ msg() [1/2]

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

Definition at line 24 of file AthCommonMsg.h.

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

◆ msg() [2/2]

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

Definition at line 27 of file AthCommonMsg.h.

27  {
28  return this->msgStream(lvl);
29  }

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

1959  {
1960  ATH_MSG_DEBUG("in postSolvinglapack()");
1961 
1962  if( z.ncol() != size) {
1963  msg(MSG::ERROR)<<"Eigenvector matrix has incorrect size : "<<z.ncol()<<" != "<<size<<endmsg;
1964  return;
1965  }
1966 
1967  if( (int)m_activeIndices.size() != size) {
1968  msg(MSG::ERROR)<<"Number of active parameters is incorrect : "<<m_activeIndices.size()<<" != "<<size<<endmsg;
1969  return;
1970  }
1971 
1972  // Compute bigvector in diagonal basis (Vb = Ut * bigvector)
1973  AlVec D(size);
1974  D = z*(*dChi2);
1975 
1976  if (m_writeEigenMat) {
1977 
1978  ATH_MSG_INFO("writing the eigenvectors in a matrix: "<< z.nrow() << "x" << z.ncol());
1979 
1980  // Set Path for the z matrix (eigenvector matrix)
1981  z.SetPathBin(m_pathbin+m_prefixName);
1982  z.SetPathTxt(m_pathtxt+m_prefixName);
1983 
1984  ATH_MSG_INFO("writing the eigenvector matrix: "<< m_scalaMatName);
1985  ATH_MSG_DEBUG("matrix will be in: "<< m_pathbin+m_prefixName+m_scalaMatName);
1986 
1987  StatusCode sc = z.Write("eigenvectors.bin",true); // write the eigenvector matrix
1988 
1989  if (sc!=StatusCode::SUCCESS)
1990  msg(MSG::ERROR)<<"Problem writing eigenvector matrix"<<endmsg;
1991 
1992  // Set Path for the w matrix (eigenvalues matrix - diagonal bigmatrix)
1993  w.SetPathBin(m_pathbin+m_prefixName);
1994  w.SetPathTxt(m_pathtxt+m_prefixName);
1995 
1996  ATH_MSG_INFO("writing the eigenvectors in a vector: "<< w.size());
1997  ATH_MSG_INFO("writing the eigenvalues vector (diagonal bigmatrix): "<< m_scalaVecName);
1998  ATH_MSG_DEBUG("vector will be in: "<< m_pathbin+m_prefixName+m_scalaVecName);
1999 
2000  sc = w.WriteEigenvalueVec("eigenvalues.bin",true); // write the eigenvalues vecor
2001 
2002  if (sc!=StatusCode::SUCCESS)
2003  msg(MSG::ERROR)<<"Problem writing eigenvector matrix"<<endmsg;
2004 
2005  if (m_writeEigenMatTxt) {
2006  sc = z.Write("eigenvectors.txt", false);
2007  if (sc!=StatusCode::SUCCESS)
2008  msg(MSG::ERROR)<<"Problem writing eigenvector matrix to text file"<<endmsg;
2009  sc = w.WriteEigenvalueVec("eigenvalues.txt", false);
2010  if (sc!=StatusCode::SUCCESS)
2011  msg(MSG::ERROR)<<"Problem writing eigenvalue vector to text file"<<endmsg;
2012  }
2013 
2014  }
2015 
2016  // Set eigenvalue thresholds
2017  const double eigenvalue_threshold = 1e-19;
2018 
2019  // weak mode removal
2020  if (m_modcut == -1) {
2021 
2022  ATH_MSG_INFO(" Starting the automatic Weak Mode Removal method");
2023 
2024  // create a pull vector for the alignment corrections in diagonal basis (db_pulls)
2025  //int nDoF=m_alignModuleTool->nAlignParameters();
2026  AlVec* Align_db = new AlVec(size);
2027  AlVec* Align_error_db = new AlVec(size);
2028  AlVec* AlignPull = new AlVec(size);
2029  ATH_MSG_DEBUG("AlignPull vector size is: "<< (*AlignPull).size());
2030 
2031  m_modcut = 0;
2032  bool wm_stop = false;
2033 
2034  // -----------------------------------------------------------------------
2035  // First pass: removing weak modes because of large pull
2036  // compute alignment pulls for corrections in diagonal basis (db)
2037  for(int i=0; i<size; i++) {
2038 
2039  (*Align_db)[i] = (-D[i]/w[i]);
2040  (*Align_error_db)[i] = sqrt(1.0/w[i]/m_scale);
2041 
2042  if (w[i]<eigenvalue_threshold) {
2043  ATH_MSG_INFO(" + EigenMode " << i
2044  << " removed as eigenvalue lower than the threshold " << eigenvalue_threshold
2045  << ": " << w[i]);
2046  (*AlignPull)[i] = 0.0;
2047  m_modcut++;
2048  }
2049  else
2050  (*AlignPull)[i] = (*Align_db)[i] / (*Align_error_db)[i];
2051 
2052  ATH_MSG_DEBUG(i << ". AlignPar: " << (*Align_db)[i] << " +- " << (*Align_error_db)[i]
2053  << " (pull: " << (*AlignPull)[i] << ") ; w[i]: " << w[i]);
2054  }
2055  ATH_MSG_INFO(" +++ Weak Mode removal ++ stop after mode "<< m_modcut << " (First pass)");
2056  // -----------------------------------------------------------------------
2057 
2058  // -----------------------------------------------------------------------
2059  // Second pass
2060  // if the error is greater than the correction -> cut this mode
2061  for(int i=m_modcut; (i<size && !wm_stop); i++) {
2062 
2063  // if the error is greater than the correction -> cut this mode
2064  if (fabs((*AlignPull)[i])<m_pullcut) {
2065  ATH_MSG_INFO(" + EigenMode " << i
2066  << " removed as pull is lower than " << m_pullcut << ": "
2067  << (*AlignPull)[i]);
2068  m_modcut++;
2069  }
2070  else
2071  wm_stop = true;
2072 
2073  }
2074  ATH_MSG_INFO(" +++ Weak Mode removal ++ stop after mode "<< m_modcut << " (Second pass)");
2075  // -----------------------------------------------------------------------
2076 
2077  wm_stop = false;
2078 
2079  // ----------------------------------------------------------------------
2080  // Third pass
2081  // Check if the next eigenvalues is far away. If it is two orders of
2082  // magnitude bigger remove also this mode and allow the search
2083  for(int i=m_modcut; (i<size && !wm_stop); i++) {
2084 
2085  // if the next eigenvalues is far away -> cut this mode
2086  if (m_eigenvalueStep*w[i]<w[i+1]) {
2087  ATH_MSG_INFO(" + EigenMode " << i
2088  << " removed as diff between eigenvalues, " << w[i] << " and " << w[i+1]
2089  << ", is greater than " << m_eigenvalueStep);
2090  m_modcut++;
2091  }
2092  else
2093  wm_stop = true;
2094 
2095  }
2096  ATH_MSG_INFO(" +++ Weak Mode removal ++ stop after mode "<< m_modcut << " (Third pass)");
2097  // -----------------------------------------------------------------------
2098 
2099  wm_stop = false;
2100 
2101  // -----------------------------------------------------------------------
2102  // Fourth pass
2103  // Check if the next eigenvalues is far away. If it is two orders of
2104  // magnitude bigger remove also this mode and allow the search
2105  for(int i=m_modcut; (i<size && !wm_stop); i++) {
2106 
2107  // if the next eigenvalues is far away -> cut this mode
2108  if ( fabs((*Align_db)[i]) > m_Align_db_step*fabs((*Align_db)[i+1]) ) {
2109  ATH_MSG_INFO(" + EigenMode " << i
2110  << " removed as diff between corrections, " << w[i] << " and " << w[i+1]
2111  << ", is greater than "
2112  << m_Align_db_step);
2113  m_modcut++;
2114  }
2115  else
2116  wm_stop = true;
2117 
2118  }
2119  ATH_MSG_INFO(" +++ Weak Mode removal ++ stop after mode "<< m_modcut << " (Fourth pass)");
2120  // -----------------------------------------------------------------------------------------------------
2121 
2122  // Free memory and clear the pointer to
2123  // prevent using invalid memory reference
2124  delete Align_db;
2125  delete Align_error_db;
2126  delete AlignPull;
2127  Align_db = nullptr;
2128  Align_error_db = nullptr;
2129  AlignPull = nullptr;
2130 
2131  } // end of if(m_modcut == -1)
2132 
2133  // Save some stuff to debug purposes
2134  /*
2135  if (m_storeDia) {
2136  std::string path = m_pathtxt+m_prefixName+"align_dia";
2137  std::fstream orthogon(path.c_str(), std::ios::out);
2138  orthogon.setf(std::ios::fixed);
2139  orthogon.setf(std::ios::showpoint);
2140  orthogon.precision(6);
2141 
2142  orthogon << std::setw(10)
2143  << "--------------------------------------------------------------------------------"
2144  << std::endl;
2145  orthogon << std::setw(10) << " ModeCut = " << m_modcut << std::endl;
2146  orthogon << std::setw(10)
2147  << "--------------------------------------------------------------------------------"
2148  << std::endl;
2149 
2150  orthogon << std::setw(10) << "mode"
2151  << std::setw(20) << "eigenmode"
2152  << std::setw(20) << "eigenmode error"
2153  << std::setw(25) << "eigenvalue"
2154  << std::endl;
2155 
2156  for( int m=0; m<size; m++) {
2157 
2158  // mode
2159  orthogon << std::setw(10) << m;
2160 
2161  // eigenmode (db)
2162  if( w[m]>1.0e-15) orthogon << std::setw(20) << -D[m]/w[m];
2163  else orthogon << std::setw(20) << 0.0;
2164 
2165  // error eigenmode (error_db)
2166  if( w[m]>1.0e-15) orthogon << std::setw(20) << sqrt(1.0/w[m]/m_scale);
2167  else orthogon << std::setw(20) << 0.0;
2168 
2169  // eigenvalues
2170  orthogon << std::setw(25) << w[m] << std::endl;
2171  }
2172  orthogon.close();
2173  } // end store align_dia.txt
2174  */
2175 
2176  AlVec delta(size);
2177  AlVec deltafull(size);
2178  AlVec errSq(size);
2179 
2180  // full covariance matrix
2181  CLHEP::HepSymMatrix * cov = nullptr;
2182  if(m_calculateFullCovariance)
2183  // Warning ! The matrix can be huge!
2184  // This can lead to memory problems
2185  cov = new CLHEP::HepSymMatrix(size,0);
2186 
2187  if(m_logStream)
2188  *m_logStream<<"/------ The Eigenvalue Spectrum -------"<<std::endl;
2189 
2190  for (int i=0;i<size;i++) {
2191  AlVec thisdelta(size);
2192  for(int j=0;j<size;j++)
2193  thisdelta[j] = z[i][j] * (-D[i]/w[i]);
2194  deltafull += thisdelta;
2195 
2196  ATH_MSG_DEBUG("eigenvalue "<<w[i]);
2197  if( i<m_modcut ) {
2198  ATH_MSG_INFO("skipping eigenvalue "<<w[i]<<" , modcut is "<<m_modcut);
2199  if(m_logStream)
2200  *m_logStream<<"| skipping eigenvalue "<<w[i]<<std::endl;
2201  }
2202  else if( w[i] < m_eigenvaluethreshold ) {
2203  ATH_MSG_INFO("skipping eigenvalue "<<w[i]<<" , cut is "<<m_eigenvaluethreshold);
2204  if(m_logStream)
2205  *m_logStream<<"| skipping eigenvalue "<<w[i]<<std::endl;
2206  }
2207  else {
2208  if(m_logStream)
2209  *m_logStream<<"| "<<w[i]<<std::endl;
2210 
2211  delta += thisdelta;
2212  for(int j=0;j<size;j++) {
2213  errSq[j] += z[i][j] * z[i][j] / w[i];
2214  if(m_calculateFullCovariance) {
2215  for(int k=0;k<=j;k++)
2216  (*cov)[j][k] += z[i][j] * z[i][k] / w[i];
2217  }
2218  }
2219  }
2220  }
2221 
2222  if(m_logStream)
2223  *m_logStream<<"\\----- End of Eigenvalue Spectrum -----"<<std::endl;
2224 
2225  ATH_MSG_DEBUG("Alignment constants:");
2226 
2227  DataVector<AlignPar> * alignParList = m_alignModuleTool->alignParList1D();
2228 
2229  // compute alignment corrections (translations in mm and rotations in rad) and their variances
2230  for(int i=0; i<size; i++) {
2231 
2232  double param = delta[i];
2233  double err = sqrt(2.*std::fabs(errSq[i]));
2234 
2235  int idof = m_activeIndices[i];
2236  AlignPar * alignPar=(*alignParList)[idof];
2237 
2238  // undo the sigma scaling
2239  double sigma = alignPar->sigma();
2240 
2241  param *= sigma;
2242  err *= sigma;
2243 
2244  // undo normalization scaling of error
2245  if(m_scaleMatrix && m_scale>0.)
2246  err /= sqrt(m_scale);
2247 
2248  ATH_MSG_DEBUG(i <<" : "<< param << " +/- "<< err);
2249  ATH_MSG_DEBUG("cov("<<i<<")="<<errSq[i]<<", sigma: "<<sigma);
2250  ATH_MSG_DEBUG("init par: "<<alignPar->initPar());
2251  alignPar->setPar(param, err);
2252  ATH_MSG_DEBUG("set param to "<<param<<" for alignPar "<<alignPar);
2253  ATH_MSG_DEBUG(*(*alignParList)[idof]);
2254  }
2255 
2256  if(m_logStream) {
2257  printGlobalSolution(*m_logStream, cov);
2258 
2259  // norm of first derivative
2260  double norm1st = dChi2->norm();
2261  if(m_scaleMatrix && m_scale>0.) // undo normalization scaling
2262  norm1st *= m_scale;
2263  *m_logStream<<"norm of first derivative : "<<norm1st<<std::endl;
2264 
2265  if(d2Chi2) {
2266  // distance to solution
2267  double dist = ( (*d2Chi2) * deltafull + (*dChi2) ).norm();
2268  if(m_scaleMatrix && m_scale>0.) // undo normalization scaling
2269  dist *= m_scale;
2270  *m_logStream<<"distance to solution : "<<dist<<std::endl;
2271 
2272  // calculate chi2 of the alignment change
2273  double chi2 = delta * (*d2Chi2) * delta * .5;
2274  if(m_scaleMatrix && m_scale>0.) // undo normalization scaling
2275  chi2 *= m_scale;
2276  *m_logStream<<"delta(chi2) of the alignment change : "<<chi2<<" / "<<size<<std::endl;
2277  }
2278  }
2279 
2280  delete cov;
2281  }

◆ prepareBinaryFiles()

void MatrixTool::prepareBinaryFiles ( int  solveOption)
virtual

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

Implements Trk::IMatrixTool.

Definition at line 260 of file MatrixTool.cxx.

261  {
262  }

◆ printGlobalSolution() [1/2]

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

Definition at line 1621 of file MatrixTool.cxx.

1622  {
1623  const AlignModuleList * alignModules = m_alignModuleTool->alignModules1D();
1624 
1625  AlignModuleList::const_iterator imod = alignModules->begin();
1626  AlignModuleList::const_iterator imod_end = alignModules->end();
1627  for( ; imod!=imod_end; ++imod) {
1628  AlignModule * module = *imod;
1629 
1630  DataVector<AlignPar> * alignPars = m_alignModuleTool->getAlignPars(module);
1631  int thisNDoF = alignPars->size();
1632 
1633  // fill local covariance matrix
1634  CLHEP::HepSymMatrix * covsub = nullptr;
1635  if(cov && module->nHits() >= m_minNumHits && module->nTracks() >= m_minNumTrks) {
1636  covsub = new CLHEP::HepSymMatrix(thisNDoF,0);
1637  for (int i=0;i<thisNDoF;++i) {
1638  int ipar = alignPars->at(i)->index();
1639  double sigma_i = alignPars->at(i)->sigma();
1640 
1641  std::vector<int>::iterator itActive = std::find(m_activeIndices.begin(),m_activeIndices.end(),ipar);
1642  if( itActive == m_activeIndices.end() )
1643  continue;
1644  int iActive = std::distance(m_activeIndices.begin(),itActive);
1645 
1646  for (int j=0;j<=i;++j) {
1647  int jpar = alignPars->at(j)->index();
1648  double sigma_j = alignPars->at(j)->sigma();
1649 
1650  std::vector<int>::iterator jtActive = std::find(m_activeIndices.begin(),m_activeIndices.end(),jpar);
1651  if( jtActive == m_activeIndices.end() )
1652  continue;
1653  int jActive = std::distance(m_activeIndices.begin(),jtActive);
1654 
1655  (*covsub)[i][j] = (*cov)[iActive][jActive] * sigma_i * sigma_j;
1656  }
1657  }
1658  }
1659 
1660  printModuleSolution(os,module,covsub);
1661 
1662  delete covsub;
1663  }
1664  os << "--------------------------------------------------------------------------------" << std::endl;
1665  }

◆ printGlobalSolution() [2/2]

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

Definition at line 1668 of file MatrixTool.cxx.

1669  {
1670  CLHEP::HepSymMatrix * cov = nullptr;
1671  if(cov0) {
1672  int nsize = cov0->GetNrows();
1673  cov = new CLHEP::HepSymMatrix(nsize,0);
1674 
1675  for(int i=0; i<nsize; i++)
1676  for(int j=0; j<=i; j++)
1677  (*cov)[i][j] = (*cov0)[i][j];
1678  }
1679 
1680  printGlobalSolution(os,cov);
1681 
1682  delete cov;
1683  }

◆ printModuleSolution()

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

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

Definition at line 1706 of file MatrixTool.cxx.

1707  {
1708  boost::io::ios_all_saver ias( os ); //save the stream state
1709 
1710  os << "--------------------------------------------------------------------------------" << std::endl;
1711  os << "Alignment parameters for module: " << module->name() << std::endl;
1712  os << "Number of tracks passing: " << module->nTracks() << std::endl;
1713  if(m_minNumHits>0 && module->nHits()<m_minNumHits) {
1714  os << "Number of hits too small: "<<module->nHits()<<" < "<<m_minNumHits<<" Skipping the module"<<std::endl;
1715  return;
1716  }
1717  if(m_minNumTrks>0 && module->nTracks()<m_minNumTrks) {
1718  os << "Number of tracks too small: "<<module->nTracks()<<" < "<<m_minNumTrks<<" Skipping the module"<<std::endl;
1719  return;
1720  }
1721  os << "Number of hits seen: " << module->nHits() << std::endl;
1722  os << "Number of tracks seen: " << module->nTracks() << std::endl;
1723 
1724  DataVector<AlignPar> * alignPars = m_alignModuleTool->getAlignPars(module);
1725  int thisNDoF = alignPars->size();
1726 
1727  if(alignPars->empty())
1728  os << "No active parameters" << std::endl;
1729  else
1730  {
1731  //RestoreIOSFlags restore_flags(os);
1732 
1733  os.unsetf(std::ios_base::floatfield);
1734  os << std::setiosflags(std::ios_base::left) << std::setprecision(5);
1735 
1736  // output alignment parameters and errors
1737  DataVector<AlignPar>::const_iterator ipar = alignPars->begin();
1738  DataVector<AlignPar>::const_iterator ipar_end = alignPars->end();
1739  for ( ; ipar != ipar_end; ++ipar) {
1740  const AlignPar * par = *ipar;
1741  os << std::setw(10) << par->dumpType()
1742  << std::setw(12) << par->par() << " +/- " << std::setw(12) << par->err()
1743  << std::endl;
1744  }
1745 
1746  if(cov) {
1747  // calculate local correlation matrix
1748  CLHEP::HepSymMatrix corrsub(thisNDoF,0);
1749  for(int irow=0; irow<thisNDoF; ++irow)
1750  for(int icol=0; icol<=irow; ++icol)
1751  corrsub[irow][icol] = (*cov)[irow][icol] / sqrt((*cov)[irow][irow] * (*cov)[icol][icol]);
1752  os << "Local correlation matrix: " << corrsub << std::flush;
1753  }
1754  }
1755  ias.restore(); //restore the stream state
1756  }

◆ readHitmaps()

void MatrixTool::readHitmaps ( )
private

Definition at line 2455 of file MatrixTool.cxx.

2456  {
2457  ATH_MSG_INFO("read hitmaps from files");
2458 
2459  const AlignModuleList * moduleList = m_alignModuleTool->alignModules1D();
2460  int nModules = moduleList->size();
2461 
2462  AlMat hitmap(nModules,2);
2463  int nFiles = (int)m_inputHitmapFiles.size();
2464  for(int imap=0;imap<nFiles; imap++) {
2465  AlMat nextmap(nModules,2);
2466 
2467  ATH_MSG_INFO("Reading hitmap "<<imap<<" from file "<<m_inputHitmapFiles[imap]);
2468 
2469  if(nextmap.ReadScalaPack(m_inputHitmapFiles[imap]).isFailure()) {
2470  ATH_MSG_WARNING("Problem reading hitmap from \'"<<m_inputHitmapFiles[imap]<<"\'. Skipping.");
2471  continue;
2472  }
2473 
2474  if(nextmap.nrow()!=nModules || nextmap.ncol()!=2) {
2475  ATH_MSG_WARNING("Matrix in file \'"<<m_inputHitmapFiles[imap]<<"\' has wrong size ("
2476  <<nextmap.nrow()<<" x "<<nextmap.ncol()<<"), should be ("<<nModules<<" x 2). Skipping.");
2477  continue;
2478  }
2479 
2480  hitmap += nextmap;
2481  }
2482 
2483  AlignModuleList::const_iterator imod = moduleList->begin();
2484  AlignModuleList::const_iterator imod_end = moduleList->end();
2485  int index = 0;
2486  int totalhits = 0;
2487  for(; imod != imod_end; ++imod) {
2488  AlignModule * module = *imod;
2489  //if ((int)hitmap[index][0]!=(int)module->identify().get_identifier32().get_compact()) ATH_MSG_ERROR("bad module identifier");
2490  //module->setIdentifier((Identifier)hitmap[index][0]);
2491  module->setNHits((int)hitmap[index][0]);
2492  module->setNTracks((int)hitmap[index][1]);
2493  totalhits += (int)hitmap[index][0];
2494  index++;
2495  }
2496 
2497  m_nHits = totalhits;
2498  m_nTracks = 0;
2499  m_nMeasurements = 0;
2500 
2501  ATH_MSG_INFO("Hitmap accumulated from "<<nFiles<<" files with total of "<<totalhits<<" hits.");
2502  }

◆ 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();
383  PBASE::renounce (h);
384  }

◆ 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  {
365  handlesArray.renounce();
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.

83 { m_nHits = n; }

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

1269  {
1270  // ============
1271  // solve
1272  // ============
1273  ATH_MSG_DEBUG("in MatrixTool::solve()");
1274 
1275  // set normalization scale to number of hits for now
1276  if(m_scale<0)
1277  m_scale = m_nHits;
1278 
1279  //-------------------------------------------------------
1280  // write matrix and vector to file
1281  if (m_writeMat) {
1282  // version has to be 2 to for reading matrices and vectors back in to work properly
1283  double dummyVersion(2.);
1284 
1285  // make map of matrix entry to module index (set by geometry manager tool)
1286  std::map<int,unsigned long long> modIndexMap;
1287  std::map<int,std::string> modNameMap;
1288  DataVector<AlignPar>* alignPars = m_alignModuleTool->alignParList1D();
1289  for (int i=0;i<(int)alignPars->size();i++) {
1290  modIndexMap[i]=(*alignPars)[i]->alignModule()->identify().get_compact();
1291  modNameMap [i]=(*alignPars)[i]->alignModule()->name();
1292  }
1293 
1294  // binary files
1295  ATH_MSG_DEBUG("writing binary files");
1296  StatusCode sc1 = m_bigmatrix->Write("matrix.bin",true,m_wSqMatrix,m_scale,dummyVersion);
1297  StatusCode sc2 = m_bigvector->WritePartial("vector.bin",true,m_scale,modIndexMap,dummyVersion);
1298  if (!sc1.isSuccess() || !sc2.isSuccess()) {
1299  msg(MSG::ERROR)<<"problem writing matrix or vector"<<endmsg;
1300  return -1;
1301  }
1302 
1303  if (m_writeMatTxt) {
1304 
1305  // text files
1306  ATH_MSG_DEBUG("writing text files");
1307  sc1 = m_bigmatrix->Write("matrix.txt",false,m_wSqMatrix,m_scale,dummyVersion);
1308  sc2 = m_writeModuleNames ?
1309  m_bigvector->WritePartial("vector.txt",false,m_scale,modNameMap,dummyVersion) :
1310  m_bigvector->WritePartial("vector.txt",false,m_scale,modIndexMap,dummyVersion);
1311 
1312  if (!sc1.isSuccess() || !sc2.isSuccess()) {
1313  msg(MSG::ERROR)<<"problem writing matrix or vector"<<endmsg;
1314  return -1;
1315  }
1316  }
1317 
1318  ATH_MSG_DEBUG("matrix and vector written to: "<<m_pathbin+m_prefixName<<"matrix.bin (.txt) and "<<m_pathbin+m_prefixName<<"vector.bin (.txt)");
1319  }
1320 
1321  //-------------------------------------------------------
1322  // write hitmap to file
1323  if (m_writeHitmap)
1324  writeHitmap();
1325 
1326  if(m_writeTFile)
1327  storeInTFile(m_pathbin+m_prefixName+m_tfileName);
1328 
1329 
1330  if(!m_runLocal && m_solveOption==0) {
1331  ATH_MSG_DEBUG("No solving requested.");
1332  return 1;
1333  }
1334 
1335  //-------------------------------------------------------
1336  // rescale the vector and the matrix according to sigmas
1337  // and apply soft mode cut
1338 
1339  ATH_MSG_DEBUG("rescaling the matrix/vector and applying the soft-mode-cut");
1340 
1341  DataVector<AlignPar>* alignParList = m_alignModuleTool->alignParList1D();
1342  int nDoF = alignParList->size();
1343 
1344 
1345  const AlSymMat * chkMatrix = dynamic_cast<const AlSymMat*>(m_bigmatrix);
1346  if(chkMatrix){
1347  // Method when using the dense matrix
1348  for (int i=0;i<nDoF;i++) {
1349  // scale the vector
1350  double sigma_i = (*alignParList)[i]->sigma();
1351  double softCut = 2 * pow( (*alignParList)[i]->softCut() , -2 );
1352  (*m_bigvector)[i] *= sigma_i;
1353 
1354  for (int j=0;j<=i;j++) {
1355  // scale the matrix
1356  if ((*chkMatrix)[i][j] != 0.) {
1357  double sigma_j = (*alignParList)[j]->sigma();
1358  (*m_bigmatrix)[i][j] *= sigma_i * sigma_j;
1359  }
1360  // apply soft-mode-cut
1361  if (i==j && m_softEigenmodeCut>0.){
1362  (*m_bigmatrix)[i][j] += m_softEigenmodeCut * softCut;
1363 
1364  }
1365 
1366  // set first and second derivatives on AlignPar
1367  if (i==j) {
1368  (*alignParList)[i]->setFirstDeriv((*m_bigvector)[i]/sigma_i);
1369  (*alignParList)[i]->setSecndDeriv((*m_bigmatrix)[i][i]/sigma_i/sigma_i);
1370  }
1371  }
1372  }
1373  } else {
1374  // Method when using the sparse matrix
1375  for (const datamap::value_type& p : *m_bigmatrix->ptrMap()) {
1376  int i = p.first.first;
1377  int j = p.first.second;
1378 
1379  // Scale matrix
1380  double sigma_i = (*alignParList)[i]->sigma();
1381  double sigma_j = (*alignParList)[j]->sigma();
1382 
1383  (*m_bigmatrix)[i][j] *= sigma_i * sigma_j;
1384 
1385  }
1386 
1387 
1388  for (int i=0;i<nDoF;i++) {
1389  // scale the vector
1390  double sigma_i = (*alignParList)[i]->sigma();
1391  (*m_bigvector)[i] *= sigma_i;
1392  if (m_softEigenmodeCut >0. ){
1393  double softCut = 2 * pow( (*alignParList)[i]->softCut() , -2 );
1394  (*m_bigmatrix)[i][i] += m_softEigenmodeCut * softCut;
1395  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]);
1396  }
1397  (*alignParList)[i]->setFirstDeriv((*m_bigvector)[i]/sigma_i);
1398  (*alignParList)[i]->setSecndDeriv((*m_bigmatrix)[i][i]/sigma_i/sigma_i);
1399  }
1400  }
1401 
1402  unsigned long long OldPixelIdentifier = 37769216; //Identifier for the Pixel Detector
1403  unsigned long long IBLIdentifier = 33574912; //Identifier for the Pixel Detector
1404 
1405  unsigned long long SCT_ECA_8_Identifier = 218116096; //Identifier for the SCT ECA Last Disk
1406  std::string SCT_ECA_8_Name = "SCT/EndcapA/Disk_8";
1407 
1408 
1409 
1410  ATH_MSG_INFO("rescaling done");
1411  ATH_MSG_INFO("Javi: Printing (*alignParList)[i]->alignModule()->identify32()");
1412 
1413  // select modules with non-zero tracks
1414  for(int i=0;i<nDoF;i++)
1415  {
1416  ATH_MSG_DEBUG(i);
1417  ATH_MSG_DEBUG((*alignParList)[i]->alignModule()->identify32());
1418  ATH_MSG_DEBUG((*alignParList)[i]->alignModule());
1419  ATH_MSG_DEBUG((*alignParList)[i]->alignModule()->name());
1420  ATH_MSG_DEBUG((*alignParList)[i]->paramType());
1421 
1422  //Skip solving for Pixel or IBL:
1423  const auto & theParameterList = *alignParList;
1424  const auto & thisIdentifier = theParameterList[i]->alignModule()->identify32();
1425  const auto & thisName = theParameterList[i]->alignModule()->name();
1426  const auto & thisParameterType = theParameterList[i]->paramType();
1427  const bool oldPixel = (thisIdentifier == OldPixelIdentifier);
1428  const bool ibl = (thisIdentifier == IBLIdentifier);
1429  const bool SCTECA8 = (thisIdentifier == SCT_ECA_8_Identifier);
1430  const bool SCTECA8_n = (thisName.find(SCT_ECA_8_Name)!= std::string::npos);
1431 
1432  if (m_DeactivateSCT_ECA_LastDisk){
1433  if (SCTECA8)
1434  {ATH_MSG_INFO( "SCT ECA Last Disk DoF have been skipped in the solving because DeactivateSCT_ECA_LastDisk is set to True");
1435  continue;}
1436  if (SCTECA8_n)
1437  {ATH_MSG_INFO( "SCT ECA Last Disk DoF have been skipped in the solving because DeactivateSCT_ECA_LastDisk is set to True");
1438  continue;}
1439  }
1440 
1441  if (m_AlignIBLbutNotPixel) //If m_AlignIBLbutNotPixel is set to True, Pixel will be skipped in the solving.
1442  if (oldPixel)
1443  {ATH_MSG_INFO( "Pixel DoF have been skipped in the solving because AlignIBLbutNotPixel is set to True");
1444  continue;}
1445 
1446  if (m_AlignPixelbutNotIBL) //If m_AlignPixelbutNotIBL is set to True, IBL will be skipped in the solving.
1447  if (ibl)
1448  {ATH_MSG_INFO( "IBL DoF have been skipped in the solving because AlignPixelbutNotIBL is set to True");
1449  continue;}
1450 
1451  //For specific DoF: (*alignParList)[i]->paramType() = 0,1,2,3,4,5,6 for Tx,Ty,Tz,Rx,Ry,Rz,Bx
1452  //Pixel Dofs:
1453  if (m_Remove_Pixel_Tx) //If m_Remove_Pixel_Tx is set to True, Pixel Tx will be skipped in the solving.
1454  if ( oldPixel and (thisParameterType == 0))
1455  {ATH_MSG_INFO( "Pixel Tx DoF has been skipped in the solving because Remove_Pixel_Tx is set to True");
1456  continue;}
1457 
1458  if (m_Remove_Pixel_Ty) //If m_Remove_Pixel_Ty is set to True, Pixel Ty will be skipped in the solving.
1459  if ( oldPixel and (thisParameterType == 1))
1460  {ATH_MSG_INFO( "Pixel Ty DoF has been skipped in the solving because Remove_Pixel_Ty is set to True");
1461  continue;}
1462 
1463  if (m_Remove_Pixel_Tz) //If m_Remove_Pixel_Tz is set to True, Pixel Tz will be skipped in the solving.
1464  if (oldPixel and (thisParameterType == 2))
1465  {ATH_MSG_INFO( "Pixel Tz DoF has been skipped in the solving because Remove_Pixel_Tz is set to True");
1466  continue;}
1467 
1468  if (m_Remove_Pixel_Rx) //If m_Remove_Pixel_Rx is set to True, Pixel Rx will be skipped in the solving.
1469  if (oldPixel and (thisParameterType == 3))
1470  {ATH_MSG_INFO( "Pixel Rx DoF has been skipped in the solving because Remove_Pixel_Rx is set to True");
1471  continue;}
1472 
1473  if (m_Remove_Pixel_Ry) //If m_Remove_Pixel_Ry is set to True, Pixel Ry will be skipped in the solving.
1474  if (oldPixel and (thisParameterType == 4))
1475  {ATH_MSG_INFO( "Pixel Ry DoF has been skipped in the solving because Remove_Pixel_Ry is set to True");
1476  continue;}
1477 
1478  if (m_Remove_Pixel_Rz) //If m_Remove_Pixel_Rz is set to True, Pixel Rz will be skipped in the solving.
1479  if (oldPixel and (thisParameterType == 5))
1480  {ATH_MSG_INFO( "Pixel Rz DoF has been skipped in the solving because Remove_Pixel_Rz is set to True");
1481  continue;}
1482 
1483  //IBL Dofs:
1484  if (m_Remove_IBL_Tx) //If m_Remove_IBL_Tx is set to True, IBL Tx will be skipped in the solving.
1485  if (ibl and (thisParameterType == 0))
1486  {ATH_MSG_INFO( "IBL Tx DoF has been skipped in the solving because Remove_IBL_Tx is set to True");
1487  continue;}
1488 
1489  if (m_Remove_IBL_Ty) //If m_Remove_IBL_Ty is set to True, IBL Ty will be skipped in the solving.
1490  if (ibl and (thisParameterType == 1))
1491  {ATH_MSG_INFO( "IBL Ty DoF has been skipped in the solving because Remove_IBL_Ty is set to True");
1492  continue;}
1493 
1494  if (m_Remove_IBL_Tz) //If m_Remove_IBL_Tz is set to True, IBL Tz will be skipped in the solving.
1495  if (ibl and (thisParameterType == 2))
1496  {ATH_MSG_INFO( "IBL Tz DoF has been skipped in the solving because Remove_IBL_Tz is set to True");
1497  continue;}
1498 
1499  if (m_Remove_IBL_Rx) //If m_Remove_IBL_Rx is set to True, IBL Rx will be skipped in the solving.
1500  if (ibl and (thisParameterType == 3))
1501  {ATH_MSG_INFO( "IBL Rx DoF has been skipped in the solving because Remove_IBL_Rx is set to True");
1502  continue;}
1503 
1504  if (m_Remove_IBL_Ry) //If m_Remove_IBL_Ry is set to True, IBL Ry will be skipped in the solving.
1505  if (ibl and (thisParameterType == 4))
1506  {ATH_MSG_INFO( "IBL Ry DoF has been skipped in the solving because Remove_IBL_Ry is set to True");
1507  continue;}
1508 
1509  if (m_Remove_IBL_Rz) //If m_Remove_IBL_Rz is set to True, IBL Rz will be skipped in the solving.
1510  if (ibl and (thisParameterType == 5))
1511  {ATH_MSG_INFO( "IBL Rz DoF has been skipped in the solving because Remove_IBL_Rz is set to True");
1512  continue;}
1513 
1514  if(theParameterList[i]->alignModule()->nHits() >= m_minNumHits && theParameterList[i]->alignModule()->nTracks() >= m_minNumTrks)
1515  m_activeIndices.push_back(i);
1516  }
1517  m_aNDoF = m_activeIndices.size();
1518  ATH_MSG_DEBUG("aNDoF/nDoF: "<<m_aNDoF<<"/"<<nDoF);
1519 
1520  // --------------------
1521  // now do the SOLVING
1522  // --------------------
1523 
1524  int info = 0;
1525 
1526  // first Local solving
1527  if (m_runLocal)
1528  info = solveLocal();
1529 
1530  // remove spurious modules and resize
1531  if (m_removeSpurious) {
1532 
1533  ATH_MSG_INFO("Spurious removal not implemented at the moment.");
1534 /* if (StatusCode::SUCCESS != spuriousRemoval()) {
1535  ATH_MSG_ERROR("Problem while trying to remove spurious. Stopping solving");
1536  return -1;
1537  }
1538 
1539  // if nDoF=0, bad job...
1540  int NDoF = m_alignModuleTool->nAlignParameters();
1541  if(NDoF==0) {
1542  ATH_MSG_WARNING("Removal removed everything: NDoF=" << NDoF << " !!!");
1543  return 1;
1544  }
1545  ATH_MSG_DEBUG("NDoF: " << NDoF);
1546 */
1547  }
1548 
1549  // --------------------
1550  // now Global solving
1551  switch(m_solveOption) {
1552 
1553  case NONE:
1554  ATH_MSG_DEBUG("No global solving requested.");
1555  break;
1556 
1557  case SOLVE_ROOT:
1558  info = solveROOT();
1559  break;
1560 
1561  case SOLVE_CLHEP:
1562  info = solveCLHEP();
1563  break;
1564 
1565  case SOLVE:
1566  case DIRECT_SOLVE:
1567  case DIRECT_SOLVE_CLUSTER:
1568  info = solveLapack();
1569  break;
1570 
1571  case SOLVE_FAST:
1572  case DIRECT_SOLVE_FAST:
1573  info = solveSparseEigen();
1574  break;
1575 
1576  default:
1577  ATH_MSG_INFO("Unknown solving option.");
1578  info = 0;
1579  break;
1580  }
1581 
1582  ATH_MSG_INFO("Return value from solving: "<<info);
1583 
1584  return info;
1585  }

◆ solveCLHEP()

int MatrixTool::solveCLHEP ( )
private

Definition at line 391 of file MatrixTool.cxx.

392  {
393  ATH_MSG_INFO("solving Global using CLHEP");
394  if(m_logStream) {
395  *m_logStream<<"*************************************************************"<<std::endl;
396  *m_logStream<<"************** solving using Global method ****************"<<std::endl;
397  *m_logStream<<"************** using CLHEP ****************"<<std::endl;
398  *m_logStream<<"*************************************************************"<<std::endl;
399  }
400 
401  // start measuring time
402  clock_t starttime = clock();
403 
404  DataVector<AlignPar> * alignParList = m_alignModuleTool->alignParList1D();
405  //const AlignModuleList* alignModules = m_alignModuleTool->alignModules1D();
406  for (int i=0;i<(int)alignParList->size();i++)
407  ATH_MSG_DEBUG("ap["<<i<<"]="<<(*alignParList)[i]);
408 
409  int nDoF = m_alignModuleTool->nAlignParameters();
410 
411  // some debugging output
412  if (msgLvl(MSG::DEBUG)) {
413  msg(MSG::DEBUG)<<"dumping matrix and vector to screen"<<endmsg;
414  for (int i=0;i<nDoF;i++)
415  for (int j=0;j<nDoF;j++)
416  //if (std::fabs((*m_bigmatrix)[i][j])>.0001)
417  msg(MSG::DEBUG)<<i<<", "<<j<<" : "<<(*m_bigmatrix)[i][j] <<endmsg;
418 
419  for (int i=0;i<nDoF;i++)
420  msg(MSG::DEBUG)<<i <<" : "<<(*m_bigvector)[i]<<endmsg;
421  }
422 
423  // get rescaled first and second derivatives
424  CLHEP::HepSymMatrix * d2Chi2 = new CLHEP::HepSymMatrix(m_aNDoF,0);
425  CLHEP::HepVector * dChi2 = new CLHEP::HepVector(m_aNDoF,0);
426  for (int iActive=0;iActive<m_aNDoF;iActive++) {
427  int i = m_activeIndices[iActive];
428  (*dChi2)[iActive] = (*m_bigvector)[i];
429  for (int jActive=0;jActive<m_aNDoF;jActive++) {
430  int j = m_activeIndices[jActive];
431  (*d2Chi2)[iActive][jActive] = (*m_bigmatrix)[i][j];
432  }
433  }
434 
435  ATH_MSG_DEBUG("First derivatives:" << (*dChi2));
436  ATH_MSG_DEBUG("Second derivatives:" << (*d2Chi2));
437 
438  CLHEP::HepSymMatrix cov(m_aNDoF,0);
439  CLHEP::HepVector delta(m_aNDoF,0);
440  CLHEP::HepVector deltafull(m_aNDoF,0);
441 
442  bool status=true;
443  int ierr(0);
444  if(!m_diagonalize) {
445  // ==========================================================
446  // Run Matrix Inversion
447  ATH_MSG_INFO("Running matrix inversion");
448  if(m_logStream)
449  *m_logStream<<"Running matrix inversion"<<std::endl;
450 
451  cov = *d2Chi2;
452  cov.invert(ierr);
453  if(ierr>0)
454  msg(MSG::ERROR)<<"CLHEP inversion status flag = "<<ierr<<endmsg;
455  else
456  ATH_MSG_INFO("CLHEP inversion OK");
457  if(m_logStream)
458  *m_logStream<<"CLHEP inversion status flag = "<<ierr<<std::endl;
459 
460  // calculate corrections
461  delta = cov * (*dChi2);
462 
463  // the covariance matrix is actually defined as 2 * d2Chi2^-1
464  ATH_MSG_DEBUG("Result: "<<delta);
465  ATH_MSG_DEBUG("cov: "<<cov*2);
466 
467  // -----------------------
468  // calculate also for matrix and vector multiplied by factor 0.5
469  // this should make no difference if everything is correct but
470  // it can go wrong if insensitive DoF is included
471  CLHEP::HepSymMatrix cov2 = *d2Chi2 * .5;
472  // invert the matrix
473  int ierr2 = 0;
474  cov2.invert(ierr2);
475  if(ierr2>0)
476  msg(MSG::WARNING)<<"Second CLHEP inversion status flag = "<<ierr2<<endmsg;
477 
478  CLHEP::HepVector delta2 = cov2 * (*dChi2) * .5;
479  for (int i=0;i<delta.num_row(); ++i)
480  if ( fabs((delta[i] - delta2[i])/delta[i]) > 1e-5 ) {
481  msg(MSG::WARNING)<<"Something's wrong with the matrix inversion: delta["<<i<<"] = "<<delta[i]<<" delta2["<<i<<"] = "<<delta2[i]<<endmsg;
482  status=false;
483  break;
484  }
485 
486  if(m_logStream && (ierr2>0 || !status)) {
487  *m_logStream<<"CLHEP inversion status flag for halfed matrix = "<<ierr2<<std::endl;
488  *m_logStream<<"Matrix inversion check failed"<<std::endl;
489  *m_logStream<<std::endl;
490  }
491  // -- end of check of matrix inversion
492  }
493  else {
494  // ==========================================================
495  // Run Diagonalization
496  ATH_MSG_INFO("Running diagonalization");
497  if(m_logStream)
498  *m_logStream<<"Running diagonalization"<<std::endl;
499 
500  CLHEP::HepSymMatrix D = *d2Chi2;
501  CLHEP::HepMatrix U = CLHEP::diagonalize( &D );
502 
503  ATH_MSG_INFO("Diagonalization done");
504  //sold = U*sdiag*U.T.
505 
506  // reorder eigenvalues ascending
507  // eigenvectors need to be reordered consistently
508  ATH_MSG_DEBUG(" Reordering eigenvalues ascending ");
509  for (int i=0; i<m_aNDoF-1; i++)
510  for (int j=i+1; j<m_aNDoF; j++)
511  if(D[j][j] < D[i][i]) {
512  // swap eigenvalues
513  double ei = D[i][i];
514  D[i][i] = D[j][j];
515  D[j][j] = ei;
516  // swap eigenvectors
517  for(int k=0;k<m_aNDoF; k++) {
518  double ev = U[k][i];
519  U[k][i] = U[k][j];
520  U[k][j] = ev;
521  }
522  }
523 
524  // how do I now get the eigenvalues? this cannot be the most
525  // efficient way ... CLHEP::HepSymMatrix D = d2Chi2->similarityT( U );
526  CLHEP::HepVector eigenvector(m_aNDoF);
527 
528  if(m_logStream)
529  *m_logStream<<"/------ The Eigenvalue Spectrum -------"<<std::endl;
530 
531  ATH_MSG_DEBUG("Calculating eigenvalues");
532  for(int imode=0; imode<m_aNDoF; ++imode) {
533 
534  // get the relevant eigenvector
535  for(int irow=0; irow<m_aNDoF; ++irow)
536  eigenvector[irow] = U[irow][imode];
537 
538  // calculate the eigenvalue
539  //double eigenvalue = d2Chi2->similarity( eigenvector );
540  double eigenvalue = D[imode][imode];
541  ATH_MSG_DEBUG("eigenvalue "<<eigenvalue);
542 
543  double evdotb = dot(*dChi2,eigenvector);
544  CLHEP::HepVector thisdelta = evdotb/eigenvalue * eigenvector;
545  deltafull += thisdelta;
546 
547  if(imode<m_modcut) {
548  ATH_MSG_INFO("skipping eigenvalue "<<imode<<" : "<<eigenvalue<<" , modcut is "<<m_modcut);
549  if(m_logStream)
550  *m_logStream<<"| skipping eigenvalue "<<eigenvalue<<std::endl;
551  }
552  else if( eigenvalue < m_eigenvaluethreshold ) {
553  ATH_MSG_INFO("skipping eigenvalue "<<eigenvalue<<" , cut is "<<m_eigenvaluethreshold);
554  if(m_logStream)
555  *m_logStream<<"| skipping eigenvalue "<<eigenvalue<<std::endl;
556  }
557  else {
558  if(m_logStream)
559  *m_logStream<<"| "<<eigenvalue<<std::endl;
560 
561  delta += thisdelta;
562 
563  // this is the time consuming part
564  for(int irow=0; irow<m_aNDoF; ++irow)
565  for(int icol=0; icol<=irow; ++icol)
566  cov[irow][icol] += eigenvector[irow] * eigenvector[icol] / eigenvalue;
567  }
568  }
569 
570  // the covariance matrix is actually defined as 2 * d2Chi2^-1
571 
572  ATH_MSG_DEBUG("Result: "<<delta);
573  ATH_MSG_DEBUG("cov: "<<cov);
574 
575  if(m_logStream)
576  *m_logStream<<"\\----- End of Eigenvalue Spectrum -----"<<std::endl;
577 
578  // end of diagonalization
579  // ==========================================================
580  }
581 
582  // stop measuring time
583  clock_t stoptime = clock();
584  double totaltime = (stoptime-starttime)/double(CLOCKS_PER_SEC);
585  ATH_MSG_INFO("Time spent in solveCLHEP: "<<totaltime<<" s");
586 
587  if(ierr==0 && status)
588  {
589  ATH_MSG_DEBUG("Alignment constants:");
590  for (int iAdof=0;iAdof<m_aNDoF;iAdof++) {
591 
592  int idof = m_activeIndices[iAdof];
593  AlignPar * alignPar=(*alignParList)[idof];
594 
595  double sigma = alignPar->sigma();
596  double param = -delta[iAdof] * sigma;
597  double err = std::sqrt(2.*std::fabs(cov[iAdof][iAdof])) * sigma;
598 
599  ATH_MSG_DEBUG(iAdof <<" : "<< param << " +/- "<< err);
600  ATH_MSG_DEBUG("cov("<<iAdof<<")="<<cov[iAdof][iAdof]<<", sigma: "<<sigma);
601  ATH_MSG_DEBUG("init par: "<<alignPar->initPar());
602  alignPar->setPar(param,err);
603  ATH_MSG_DEBUG("set param to "<<param<<" for alignPar "<<alignPar);
604  ATH_MSG_DEBUG(*(*alignParList)[idof]);
605  }
606 
607  if(m_logStream) {
608  printGlobalSolution(*m_logStream,&cov);
609 
610  // norm of first derivative
611  *m_logStream<<"norm of first derivative : "<<dChi2->norm()<<std::endl;
612 
613  // distance to solution
614  double dist = ( - (*d2Chi2) * deltafull + (*dChi2) ).norm();
615  *m_logStream<<"distance to solution : "<<dist<<std::endl;
616 
617  // calculate chi2 of the alignment change
618  double chi2 = d2Chi2->similarity(delta) * .5;
619  *m_logStream<<"delta(chi2) of the alignment change : "<<chi2<<" / "<<m_aNDoF<<std::endl;
620 
621  // time spent here
622  *m_logStream<<"time spent in solve : "<<totaltime<<" s"<<std::endl;
623  }
624  }
625 
626  delete d2Chi2;
627  delete dChi2;
628 
629  return 1;
630  }

◆ solveLapack()

int MatrixTool::solveLapack ( )
private

Definition at line 1765 of file MatrixTool.cxx.

1766  {
1767  ATH_MSG_INFO("solving Global using Lapack");
1768  if(m_logStream) {
1769  *m_logStream<<"*************************************************************"<<std::endl;
1770  *m_logStream<<"************** solving using Global method ****************"<<std::endl;
1771  *m_logStream<<"************** using LAPACK ****************"<<std::endl;
1772  *m_logStream<<"*************************************************************"<<std::endl;
1773  }
1774 
1775  // get rescaled first and second derivatives
1776  AlSymMat* aBetterMat = new AlSymMat(m_aNDoF);
1777  AlVec* aBetterVec = new AlVec(m_aNDoF);
1778  for (int iActive=0;iActive<m_aNDoF;iActive++) {
1779  int i = m_activeIndices[iActive];
1780  (*aBetterVec)[iActive] = (*m_bigvector)[i];
1781  for (int jActive=0;jActive<m_aNDoF;jActive++) {
1782  int j = m_activeIndices[jActive];
1783  (*aBetterMat)[iActive][jActive] = (*m_bigmatrix)[i][j];
1784  }
1785  }
1786 
1787  // normalize bigmatrix and bigvector
1788  if(m_scaleMatrix) {
1789  if(m_scale<=0.)
1790  ATH_MSG_WARNING("Scaling requested but scale not set. Not scaling matrix and vector.");
1791  else {
1792  (*aBetterVec) *= 1./m_scale;
1793  (*aBetterMat) *= 1./m_scale;
1794  }
1795  }
1796 
1797  ATH_MSG_DEBUG("Now Solving alignment using lapack diagonalization routine dspev...");
1798 
1799  if (m_calDet) {
1800  const double tol = 1.e-20;
1801  // compute final determinant
1802  double determ = (*aBetterMat).determinant();
1803  ATH_MSG_INFO("Determinant: " << determ);
1804  if (fabs(determ) < tol)
1805  ATH_MSG_WARNING("Matrix is singular!");
1806  }
1807 
1808  // store the original matrix for checks
1809  AlSymMat * d2Chi2 = nullptr;
1810  if (m_calculateFullCovariance)
1811  d2Chi2 = new AlSymMat(*aBetterMat);
1812 
1813  clock_t starttime = clock();
1814 
1815  // declare transition matrix + vector to store eigenvalues
1816  AlMat z(m_aNDoF,m_aNDoF);
1817  AlVec w(m_aNDoF); // vector to store the eigenvalues
1818  ATH_MSG_DEBUG("MatrixTool::after z/w allocation");
1819 
1820  char jobz = 'V';
1821  int info = (*aBetterMat).diagonalize(jobz,w,z);
1822  ATH_MSG_DEBUG(" info: " << info);
1823  ATH_MSG_INFO("MatrixTool::after diagonalization");
1824 
1825  // stop time calculation
1826  clock_t stoptime = clock();
1827  double time_diag = (stoptime-starttime)/double(CLOCKS_PER_SEC);
1828  ATH_MSG_INFO(" - time spent diagonalizing the matrix: "<<time_diag<<" s");
1829 
1830  double time_solve = 0.;
1831  if (info==0) {
1832  starttime = clock();
1833  postSolvingLapack(aBetterVec,d2Chi2,w,z,m_aNDoF);
1834  stoptime = clock();
1835  time_solve = (stoptime-starttime)/double(CLOCKS_PER_SEC);
1836  ATH_MSG_INFO(" - time spent solving the system: "<<time_solve<<" s");
1837  if(m_logStream) {
1838  *m_logStream<<"time spent for diagonalization: "<<time_diag<<" s"<<std::endl;
1839  *m_logStream<<"time spent for post-solving: "<<time_solve<<" s"<<std::endl;
1840  }
1841  }
1842  else {
1843  ATH_MSG_ERROR("Problem in diagonalization. Solving skipped.");
1844  if(m_logStream)
1845  *m_logStream<<"time spent for diagonalization: "<<time_diag<<" s"<<std::endl;
1846  }
1847 
1848  if(m_logStream) {
1849  *m_logStream<<"total time spent in solve: "<<time_diag+time_solve<<" s"<<std::endl;
1850  *m_logStream<<std::endl;
1851  }
1852 
1853  delete d2Chi2;
1854  delete aBetterMat;
1855  delete aBetterVec;
1856 
1857  // need to do this since success return value from Lapack is 0
1858  // and from solveLapack() it is 1
1859  if (info==0)
1860  info = 1;
1861 
1862  return info;
1863  }

◆ solveLocal()

int MatrixTool::solveLocal ( )
private

Definition at line 633 of file MatrixTool.cxx.

634  {
635  ATH_MSG_INFO("solving using Local method");
636  if(m_logStream) {
637  *m_logStream<<"*************************************************************"<<std::endl;
638  *m_logStream<<"************** solving using Local method *****************"<<std::endl;
639  *m_logStream<<"*************************************************************"<<std::endl;
640  }
641 
642  int totalNDoF(0);
643  double totalChi2(0.);
644 
645  const AlignModuleList * alignModules = m_alignModuleTool->alignModules1D();
646 
647  AlignModuleList::const_iterator imod = alignModules->begin();
648  AlignModuleList::const_iterator imod_end = alignModules->end();
649  for( ; imod!=imod_end; ++imod) {
650 
651  AlignModule * module = *imod;
652 
653  ATH_MSG_INFO("Solving for module: "<<module->name());
654 
655  DataVector<AlignPar> * alignPars = m_alignModuleTool->getAlignPars(module);
656 
657  int thisNDoF = alignPars->size();
658 
659  CLHEP::HepSymMatrix d2Chi2(thisNDoF,0);
660  CLHEP::HepVector dChi2(thisNDoF,0);
661  for (int i=0;i<thisNDoF;++i) {
662  int ipar = alignPars->at(i)->index();
663  dChi2[i] = (*m_bigvector)[ipar];
664  for (int j=0;j<thisNDoF;++j) {
665  int jpar = alignPars->at(j)->index();
666  d2Chi2[i][j] = (*m_bigmatrix)[ipar][jpar];
667  }
668  }
669 
670  ATH_MSG_DEBUG("First derivatives:" << dChi2);
671  ATH_MSG_DEBUG("Second derivatives:" << d2Chi2);
672 
673  if(module->nHits() < m_minNumHits || module->nTracks() < m_minNumTrks) {
674  ATH_MSG_INFO("Not enough hits in module \'"<<module->name()<<"\': "
675  <<module->nHits()<<" < "<<m_minNumHits<<" or "
676  <<module->nTracks()<<" < "<<m_minNumTrks
677  <<". Skipping...");
678 
679  // print module summary even when solving is not done
680  if(m_logStream)
681  printModuleSolution(*m_logStream,module,nullptr);
682 
683  continue;
684  }
685 
686  ATH_MSG_DEBUG("First derivatives:" << dChi2);
687  ATH_MSG_DEBUG("Second derivatives:" << d2Chi2);
688 
689 
690  totalNDoF += thisNDoF;
691 
692  CLHEP::HepSymMatrix cov(d2Chi2);
693 
694  // invert the matrix
695  int ierr = 0;
696  cov.invert(ierr);
697  if(ierr>0)
698  ATH_MSG_WARNING("CLHEP inversion status flag = "<<ierr);
699  else
700  ATH_MSG_DEBUG("CLHEP inversion status flag = "<<ierr);
701 
702  // calculate corrections
703  CLHEP::HepVector delta = cov * dChi2;
704  //ATH_MSG_DEBUG("d2Chi2: "<<d2Chi2);
705  //ATH_MSG_DEBUG("cov: "<<cov);
706  //ATH_MSG_DEBUG("d2Chi2*cov: "<<d2Chi2*cov);
707  //ATH_MSG_DEBUG("dChi2: "<<dChi2);
708  ATH_MSG_DEBUG("Result: "<<delta);
709 
710  ATH_MSG_DEBUG("Alignment constants:");
711  for (int idof=0;idof<thisNDoF;++idof) {
712  AlignPar * alignPar = alignPars->at(idof);
713 
714  double sigma = alignPar->sigma();
715  double param = -delta[idof] * sigma;
716  double err = std::sqrt(2.*std::fabs(cov[idof][idof])) * sigma;
717 
718  ATH_MSG_DEBUG(idof <<" : "<< param << " +/- "<< err);
719  ATH_MSG_DEBUG("cov("<<idof<<")="<<cov[idof][idof]<<", sigma: "<<sigma);
720  ATH_MSG_DEBUG("init par: "<<alignPar->initPar());
721 
722  ATH_MSG_DEBUG("Filling constants obtained using Local method");
723  alignPar->setPar(param,err);
724  ATH_MSG_DEBUG("set param to "<<param<<" for alignPar "<<alignPar);
725  ATH_MSG_DEBUG(*alignPar);
726  }
727 
728  if(m_logStream) {
729  printModuleSolution(*m_logStream,module,&cov);
730 
731  *m_logStream<<"CLHEP inversion status flag = "<<ierr<<std::endl;
732 
733  // calculate chi2 of the alignment change
734  double chi2 = d2Chi2.similarity(delta) * .5;
735  totalChi2 += chi2;
736  *m_logStream<<"delta(chi2) of the alignment change : "<<chi2<<" / "<<thisNDoF<<std::endl;
737  }
738  }
739 
740  if(m_logStream) {
741  *m_logStream<<"--------------------------------------------------------------------------------"<<std::endl;
742  *m_logStream<<"Total delta(chi2) of the alignment change from the local method : "<<totalChi2<<" / "<<totalNDoF<<std::endl;
743  *m_logStream<<std::endl;
744  }
745 
746  return 1;
747  }

◆ solveROOT()

int MatrixTool::solveROOT ( )
private

Definition at line 265 of file MatrixTool.cxx.

266  {
267  ATH_MSG_INFO("solving Global using ROOT");
268  if(m_logStream) {
269  *m_logStream<<"*************************************************************"<<std::endl;
270  *m_logStream<<"************** solving using Global method ****************"<<std::endl;
271  *m_logStream<<"************** using ROOT ****************"<<std::endl;
272  *m_logStream<<"*************************************************************"<<std::endl;
273  }
274 
275  // start measuring time
276  clock_t starttime = clock();
277 
278  DataVector<AlignPar>* alignParList = m_alignModuleTool->alignParList1D();
279  //const AlignModuleList* alignModules = m_alignModuleTool->alignModules1D();
280 
281  int nDoF=m_alignModuleTool->nAlignParameters();
282 
283  // some debugging output
284  if (msgLvl(MSG::VERBOSE)) {
285  msg(MSG::VERBOSE)<<"dumping matrix and vector to screen"<<endmsg;
286  for (int i=0;i<nDoF;i++)
287  for (int j=0;j<nDoF;j++)
288  //if (std::fabs((*m_bigmatrix)[i][j])>.0001)
289  msg(MSG::VERBOSE)<<i<<", "<<j<<" : "<<(*m_bigmatrix)[i][j] <<endmsg;
290 
291  for (int i=0;i<nDoF;i++)
292  msg(MSG::VERBOSE)<<i <<" : "<<(*m_bigvector)[i]<<endmsg;
293  }
294 
295  // get rescaled first and second derivatives
296  double * secderiv = new double[m_aNDoF*m_aNDoF];
297  double * firstderiv = new double[m_aNDoF];
298  for (int iActive=0;iActive<m_aNDoF;iActive++) {
299  int i = m_activeIndices[iActive];
300  firstderiv[iActive] = (*m_bigvector)[i];
301  for (int jActive=0;jActive<m_aNDoF;jActive++) {
302  int j = m_activeIndices[jActive];
303  secderiv[iActive*m_aNDoF+jActive] = (*m_bigmatrix)[i][j];
304  }
305  }
306 
307  // attention, the dimension of matrix a and b is m_aNDoF not nDoF,
308  // this means some alignment parameters have not been calculated
309  // if the corresponding modules did not satify the select cut
310 
311  TMatrixDSym a(m_aNDoF,secderiv);
312  TVectorD b(m_aNDoF,firstderiv);
313 
314  if(msgLvl(MSG::DEBUG)) {
315  msg(MSG::DEBUG)<<"First derivatives:"<<endmsg;
316  b.Print();
317  msg(MSG::DEBUG)<<"Second derivatives:"<<endmsg;
318  a.Print();
319  }
320 
321  TDecompBK c(a);
322  Bool_t status;
323  TMatrixDSym ainv(c.Invert(status));
324 
325  TVectorD r(b.GetNrows());
326  if(status)
327  r = c.Solve(b,status);
328 
329  // stop measuring time
330  clock_t stoptime = clock();
331  double totaltime = (stoptime-starttime)/double(CLOCKS_PER_SEC);
332  ATH_MSG_INFO("Time spent in solveROOT: "<<totaltime<<" s");
333 
334  if(!status) {
335  msg(MSG::ERROR)<<"ROOT inversion failed"<<endmsg;
336  if(m_logStream) {
337  *m_logStream<<"ROOT inversion failed"<<std::endl;
338  *m_logStream<<std::endl;
339  }
340  }
341  else {
342  ATH_MSG_INFO("ROOT inversion ok");
343 
344  ATH_MSG_DEBUG("Alignment constants:");
345  for (int iAdof=0;iAdof<m_aNDoF;iAdof++) {
346 
347  int idof = m_activeIndices[iAdof];
348  AlignPar * alignPar=(*alignParList)[idof];
349 
350  double sigma = alignPar->sigma();
351  double param = -r[iAdof] * sigma;
352  double err = std::sqrt(2.*std::fabs(ainv(iAdof,iAdof))) * sigma;
353 
354  ATH_MSG_DEBUG(iAdof <<" : "<< param << " +/- "<< err);
355  ATH_MSG_DEBUG("ainv("<<iAdof<<")="<<ainv(iAdof,iAdof)<<", sigma: "<<sigma);
356  ATH_MSG_DEBUG("init par: "<<alignPar->initPar());
357  alignPar->setPar(param,err);
358  ATH_MSG_DEBUG("set param to "<<param<<" for alignPar "<<alignPar);
359  ATH_MSG_DEBUG(*(*alignParList)[idof]);
360  }
361 
362  if(m_logStream)
363  {
364  *m_logStream<<"ROOT inversion ok"<<std::endl;
365 
366  printGlobalSolution(*m_logStream,&ainv);
367 
368  // norm of first derivative
369  *m_logStream<<"norm of first derivative : "<<sqrt(b.Norm2Sqr())<<std::endl;
370 
371  // distance to solution
372  double dist = sqrt( ( b - (a * r) ).Norm2Sqr() );
373  *m_logStream<<"distance to solution : "<<dist<<std::endl;
374 
375  // calculate chi2 of the alignment change
376  double chi2 = a.Similarity(r) * .5;
377  *m_logStream<<"delta(chi2) of the alignment change : "<<chi2<<" / "<<m_aNDoF<<std::endl;
378 
379  // time spent here
380  *m_logStream<<"time spent in solve : "<<totaltime<<" s"<<std::endl;
381  }
382  }
383 
384  delete [] secderiv;
385  delete [] firstderiv;
386 
387  return 1;
388  }

◆ solveSparseEigen()

int MatrixTool::solveSparseEigen ( )
private

Definition at line 2284 of file MatrixTool.cxx.

2285  {
2286  ATH_MSG_INFO("solving Global using SparseEigen");
2287  if(m_logStream) {
2288  *m_logStream<<"*************************************************************"<<std::endl;
2289  *m_logStream<<"************** solving using Global method ****************"<<std::endl;
2290  *m_logStream<<"************** using SparseEigen ****************"<<std::endl;
2291  *m_logStream<<"*************************************************************"<<std::endl;
2292  }
2293 
2294  // start measuring time
2295  clock_t starttime = clock();
2296 
2297  DataVector<AlignPar> * alignParList = m_alignModuleTool->alignParList1D();
2298 
2299  AlSpaMat * ABetterMat = nullptr;
2300  bool isCopy = false;
2301  if ( dynamic_cast<AlSymMat*>(m_bigmatrix) ) {
2302  ATH_MSG_INFO("Converting Matrix Format for fast solving");
2303  ABetterMat = new AlSpaMat(*(dynamic_cast<AlSymMat*>(m_bigmatrix)));
2304  isCopy = true;
2305  }
2306  else if ( dynamic_cast<AlSpaMat*>(m_bigmatrix) ) {
2307  ATH_MSG_INFO("Matrix format native to the fast solving");
2308  ABetterMat = (dynamic_cast<AlSpaMat*>(m_bigmatrix));
2309  }
2310  else {
2311  ATH_MSG_ERROR("Cannot cast to neither AlSymMat nor AlSpaMat");
2312  return 0;
2313  }
2314 
2315  ATH_MSG_DEBUG("checking active indices");
2316 
2317  // use const matrix when checking for non-zero elements to avoid
2318  // filling of the whole matrix
2319  const AlSpaMat * chkMatrix = ABetterMat;
2320 
2321  AlSpaMat * aBetterMat = new AlSpaMat(m_aNDoF);
2322  AlVec * aBetterVec = new AlVec(m_aNDoF);
2323 
2324 
2325  for (int iActive=0;iActive<m_aNDoF;iActive++) {
2326  int i = m_activeIndices[iActive];
2327  (*aBetterVec)[iActive] = (*m_bigvector)[i];
2328  for (int jActive=0;jActive<m_aNDoF;jActive++) {
2329  int j = m_activeIndices[jActive];
2330  // only fill if non-zero !!!
2331  if ( (*chkMatrix)[iActive][jActive] != 0. )
2332  (*aBetterMat)[iActive][jActive]=(*ABetterMat)[i][j];
2333  }
2334  }
2335 
2336  // store original vector for cross-checks
2337  AlVec origVec(*aBetterVec);
2338 
2339  ATH_MSG_DEBUG("running the solving");
2340 
2341  // solve
2342  int info = (*aBetterMat).SolveWithEigen(*aBetterVec);
2343 
2344  if(info == 0) {
2345  ATH_MSG_INFO("SolveWithEigen solving OK");
2346  if(m_logStream)
2347  *m_logStream<<"SolveWithEigen solving OK."<<std::endl;
2348  }
2349  else {
2350  ATH_MSG_ERROR( "SolveWithEigen error code (0 if OK) = "<<info );
2351  if(m_logStream)
2352  *m_logStream<<"SolveWithEigen error code (0 if OK) = "<<info<<std::endl;
2353  }
2354 
2355  if( isCopy )
2356  delete ABetterMat;
2357  ABetterMat = nullptr;
2358 
2359  // stop measuring time
2360  clock_t stoptime = clock();
2361  double totaltime = (stoptime-starttime)/double(CLOCKS_PER_SEC);
2362  ATH_MSG_INFO("Time spent in SolveWithEigen: "<<totaltime<<" s");
2363 
2364  ATH_MSG_DEBUG("Alignment constants:");
2365  // compute alignment corrections (translations in mm and rotations in rad)
2366  // for solveSparseEigen variances are not calculated
2367  for(int i=0; i<m_aNDoF; i++) {
2368 
2369  double param = -(*aBetterVec)[i];
2370  double err = 0.;
2371 
2372  int idof = m_activeIndices[i];
2373  AlignPar * alignPar=(*alignParList)[idof];
2374 
2375  // undo the sigma scaling
2376  double sigma = alignPar->sigma();
2377  param *= sigma;
2378 
2379  ATH_MSG_DEBUG(i <<" : "<< param << " +/- "<< err);
2380  ATH_MSG_DEBUG("sigma: "<<sigma);
2381  ATH_MSG_DEBUG("init par: "<<alignPar->initPar());
2382  alignPar->setPar(param, err);
2383  ATH_MSG_DEBUG("set param to "<<param<<" for alignPar "<<alignPar);
2384  ATH_MSG_DEBUG(*(*alignParList)[idof]);
2385  }
2386 
2387  if(m_logStream) {
2388  CLHEP::HepSymMatrix * cov = nullptr;
2389  printGlobalSolution(*m_logStream, cov);
2390 
2391  // norm of first derivative
2392  *m_logStream<<"norm of first derivative : "<<origVec.norm()<<std::endl;
2393 
2394  // distance to solution
2395  double dist = ( (*aBetterMat) * (*aBetterVec) - origVec ).norm();
2396  *m_logStream<<"distance to solution : "<<dist<<std::endl;
2397 
2398  // calculate chi2 of the alignment change
2399  double chi2 = (*aBetterVec) * (*aBetterMat) * (*aBetterVec) * .5;
2400  *m_logStream<<"delta(chi2) of the alignment change : "<<chi2<<" / "<<m_aNDoF<<std::endl;
2401 
2402  // time spent here
2403  *m_logStream<<"time spent in solve : "<<totaltime<<" s"<<std::endl;
2404  }
2405 
2406  delete aBetterMat;
2407  delete aBetterVec;
2408 
2409  // need to do this since success return value from Lapack is 0
2410  // and from SolveWithEigen() it is 1
2411  if (info==0)
2412  info = 1;
2413 
2414  return info;
2415  }

◆ spuriousRemoval()

StatusCode MatrixTool::spuriousRemoval ( )
private

Definition at line 1866 of file MatrixTool.cxx.

1867  {
1868 
1869  // copied from SiGlobalChi2Algs
1870  ATH_MSG_DEBUG("in spuriousRemoval");
1871 
1872  // compute determinant before resizing
1873  if (m_calDet) {
1874  const double tol = 1.e-20;
1875  double determ = m_bigmatrix->determinant();
1876  ATH_MSG_INFO("Determinant: " << determ);
1877  if (std::fabs(determ) < tol)
1878  ATH_MSG_WARNING("Matrix is singular!");
1879  }
1880 
1881  // fill vector with modules that need to be removed
1882  int fillvecmods=fillVecMods();
1883  if (fillvecmods==0) {
1884 
1885  //ATH_MSG_INFO(" No resize needed (NhitsCut = "
1886  // << m_hitscut << ")");
1887 
1888  if (msgLvl(MSG::DEBUG)) {
1889  //m_bigmatrix->Write("bigmatrix.txt", false, m_wSqMatrix, m_scale,
1890  // MatVersion);
1891  //m_bigvector->Write("bigvector.txt", false, m_scale, m_modcodemap,
1892  // VecVersion, m_alignProcessLevel, m_fitParam);
1893  }
1894 
1895  return StatusCode::SUCCESS;
1896  }
1897  else if (fillvecmods==2)
1898  return StatusCode::FAILURE;
1899 
1900  /* this is a bit difficult to implement for now....
1901 
1902  // remove matrix/vector elements
1903  int cont=0;
1904  ModuleIndexMap::iterator itcode;
1905  ATH_MSG_INFO("Eliminating module...");
1906 
1907 
1908 
1909  for (std::vector<int>::const_iterator it=m_dropmods.begin();
1910  it!= m_dropmods.end(); ++it) {
1911 
1912  itcode = m_modcodemap.find((*it));
1913 
1914  if (itcode == m_modcodemap.end()) {
1915  ATH_MSG_WARNING("Could not find module " << *it << " in map.");
1916  return StatusCode::FAILURE;
1917  }
1918 
1919  ATH_MSG_INFO(" - Removing mcode: " << (itcode->second)
1920  << " (old index: " << (*it) << " -> new index: " << (*it)-cont << ")");
1921 
1922  m_bigmatrix->RemoveModule((*it)-cont);
1923  m_bigvector->RemoveModule((*it)-cont);
1924  cont++;
1925  }
1926  ATH_MSG_INFO("Modules removed from the Matrix and from the Vector Successfully!");
1927  ATH_MSG_DEBUG("(DoF: " << nDoF << ")");
1928 
1929  // resizing...
1930  ATH_MSG_INFO("Resizing the bigvector in memory...");
1931  m_bigvector->reSize(nDoF-6*m_dropmods.size());
1932  ATH_MSG_INFO("Resizing the bigmatrix in memory...");
1933  m_bigmatrix->reSize(nDoF-6*m_dropmods.size());
1934 
1935  ATH_MSG_INFO(m_dropmods.size() << " modules eliminated from the matrix, i.e, "
1936  << 6*m_dropmods.size() << " DoFs");
1937  ATH_MSG_INFO(nDoF/6 - m_dropmods.size() << " modules to align (" << nDoF-6*m_dropmods.size() << " DoFs)");
1938  ATH_MSG_INFO("New bigmatrix size is: " << m_bigmatrix->size());
1939  ATH_MSG_INFO("New bigvector size is: " << (*m_bigvector).size());
1940 
1941  NDoF = nDoF-6*(m_dropmods.size());
1942 
1943  // resize (update) nDoF variable
1944  nDoF = NDoF;
1945 
1946  // Resizing vectors to store results
1947  m_alignPar->reSize(nDoF);
1948  m_alignSqErr->reSize(nDoF);
1949 
1950  // Fill the mapping module map and updating m_modcodemap
1951  UpdateModMap();
1952  */
1953 
1954  return StatusCode::SUCCESS;
1955  }

◆ storeInTFile()

void MatrixTool::storeInTFile ( const TString &  filename)

Store Files in a tfile.

Definition at line 862 of file MatrixTool.cxx.

863  {
864  //Store reults in a single TFile....
865  //Including Matrix Vector Hitmap.. Soluton EVs etc.
866 
867  ATH_MSG_DEBUG("Writing Results to a TFile");
868 
869  DataVector<AlignPar> * alignParList = m_alignModuleTool->alignParList1D();
870  int nDoF = alignParList->size();
871 
872  TMatrixDSparse* myTMatrix = m_bigmatrix->makeTMatrix();
873 
874  ATH_MSG_DEBUG( "Created TMatrixDSparse" );
875 
876 
877  double *val = new double[nDoF];
878  for (int i=0;i<nDoF;i++) {
879  val[i] = (*m_bigvector)[i];
880  }
881 
882  TVectorD myTVector(nDoF, val);
883  delete [] val;
884 
885  ATH_MSG_DEBUG( "Created TVectorD" );
886 
887 
888  const AlignModuleList * moduleList = m_alignModuleTool->alignModules1D();
889  int nModules = moduleList->size();
890 
891  double *hitmapA = new double[nModules];
892  double *hitmapB = new double[nModules];
893  AlignModuleList::const_iterator imod = moduleList->begin();
894  AlignModuleList::const_iterator imod_end = moduleList->end();
895  int index(0);
896  for(; imod != imod_end; ++imod) {
897  AlignModule * module = *imod;
898  hitmapA[index] = (double)module->nHits();
899  hitmapB[index] = (double)module->nTracks();
900  index++;
901  }
902 
903  TVectorD hitmapHits(nModules, hitmapA);
904  TVectorD hitmapTracks(nModules, hitmapB);
905 
906  delete [] hitmapA;
907  delete [] hitmapB;
908 
909 
910 
911  TFile myFile(filename,"recreate");
912  hitmapHits.Write("Hits");
913  hitmapTracks.Write("Tracks");
914  myTMatrix->Write("Matrix");
915  myTVector.Write("Vector");
916 
917  TVectorD scale(1, &m_scale) ;
918  scale.Write("Scale");
919 
920 
921  double *moduleInfoA = new double[nDoF];//unsigned long long
922  double *dofInfoA = new double[nDoF];//int
923 
924  if (sizeof(unsigned long long) != sizeof(double))
925  ATH_MSG_ERROR("Module Identifiers will not be saved. sizeof(double)!=sizeof(ulonglong)");
926  else{
927 
928  DataVector<AlignPar>* alignPars = m_alignModuleTool->alignParList1D();
929  for (int i=0;i<(int)alignPars->size();i++) {
930  //Do a direct memory copy to store unsigned long long in the momory of a double
931  double target;
932  uint64_t id = (*alignPars)[i]->alignModule()->identify().get_compact();
933  memcpy(&target, &id, sizeof(target));
934  moduleInfoA[i]=target;
935  //moduleInfoB[i]=(*alignPars)[i]->alignModule()->name();
936  uint64_t dof = (*alignPars)[i]->paramType();
937  memcpy(&target, &dof, sizeof(target));
938  dofInfoA[i]=target;
939  }
940 
941  TVectorD moduleIDs(nDoF, moduleInfoA) ;
942  TVectorD moduleDoFs(nDoF,dofInfoA);
943  delete [] moduleInfoA;
944  moduleIDs.Write("ModuleID");
945  moduleDoFs.Write("dof");
946  }
947 
948  myFile.Write();
949  myFile.Close();
950 
951  delete myTMatrix;
952  ATH_MSG_DEBUG("Finshed writing TFILE");
953 
954  }

◆ 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 DerivationFramework::CfAthAlgTool, AthCheckedComponent< AthAlgTool >, AthCheckedComponent<::AthAlgTool >, and asg::AsgMetadataTool.

◆ 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) {
312  std::vector<SG::VarHandleKey*> keys = a->keys();
313  for (auto k : keys) {
314  k->setOwner(this);
315  }
316  }
317  }

◆ writeHitmap()

void MatrixTool::writeHitmap ( )
private

Definition at line 2418 of file MatrixTool.cxx.

2419  {
2420  ATH_MSG_INFO("writing the hitmap to file");
2421 
2422  const AlignModuleList * moduleList = m_alignModuleTool->alignModules1D();
2423  int nModules = moduleList->size();
2424 
2425  AlMat hitmap(nModules,2);
2426  AlignModuleList::const_iterator imod = moduleList->begin();
2427  AlignModuleList::const_iterator imod_end = moduleList->end();
2428  int index(0);
2429  for(; imod != imod_end; ++imod) {
2430  AlignModule * module = *imod;
2431  hitmap[index][0] = module->nHits();
2432  hitmap[index][1] = module->nTracks();
2433  index++;
2434  }
2435 
2436  // Set Path for the hitmap matrix
2437  hitmap.SetPathBin(m_pathbin+m_prefixName);
2438  hitmap.SetPathTxt(m_pathtxt+m_prefixName);
2439 
2440  StatusCode sc = hitmap.Write("hitmap.bin",true); // write the hitmap matrix
2441 
2442  if (sc!=StatusCode::SUCCESS)
2443  ATH_MSG_ERROR("Problem writing hitmap matrix");
2444 
2445  if (m_writeHitmapTxt) {
2446  sc = hitmap.Write("hitmap.txt", false, 0);
2447  if (sc!=StatusCode::SUCCESS)
2448  ATH_MSG_ERROR("Problem writing hitmap matrix to text file");
2449  }
2450 
2451  ATH_MSG_DEBUG("hitmap written to: "<< m_pathbin+m_prefixName <<"hitmap.bin (.txt)");
2452  }

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

◆ m_Align_db_step

float Trk::MatrixTool::m_Align_db_step
private

corr in the diagonal basis step for the third pass in the auto weak mode removal method

Definition at line 162 of file MatrixTool.h.

◆ m_AlignIBLbutNotPixel

bool Trk::MatrixTool::m_AlignIBLbutNotPixel
private

Definition at line 214 of file MatrixTool.h.

◆ m_alignModuleMap

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

Definition at line 105 of file IMatrixTool.h.

◆ m_alignModuleTool

ToolHandle<IAlignModuleTool> Trk::MatrixTool::m_alignModuleTool
private

Definition at line 142 of file MatrixTool.h.

◆ m_AlignPixelbutNotIBL

bool Trk::MatrixTool::m_AlignPixelbutNotIBL
private

Definition at line 215 of file MatrixTool.h.

◆ m_aNDoF

int Trk::MatrixTool::m_aNDoF
private

number of active DoF (size of m_activeIndices)

Definition at line 210 of file MatrixTool.h.

◆ m_bigmatrix

AlSymMatBase* Trk::MatrixTool::m_bigmatrix
private

matrix to contain second derivative terms to be used for alignment

Definition at line 145 of file MatrixTool.h.

◆ m_bigvector

AlVec* Trk::MatrixTool::m_bigvector
private

vector to contain first derivative terms to be used for alignment

Definition at line 148 of file MatrixTool.h.

◆ m_calculateFullCovariance

bool Trk::MatrixTool::m_calculateFullCovariance
private

Definition at line 188 of file MatrixTool.h.

◆ m_calDet

bool Trk::MatrixTool::m_calDet
private

Compute bigmatrix's determinant ?

Definition at line 164 of file MatrixTool.h.

◆ m_DeactivateSCT_ECA_LastDisk

bool Trk::MatrixTool::m_DeactivateSCT_ECA_LastDisk
private

Definition at line 216 of file MatrixTool.h.

◆ 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

bool Trk::MatrixTool::m_diagonalize
private

run diagonalization instead of inversion

Definition at line 153 of file MatrixTool.h.

◆ m_eigenvalueStep

float Trk::MatrixTool::m_eigenvalueStep
private

eigenvalue step for the second pass in the automatic weak mode removal method

Definition at line 161 of file MatrixTool.h.

◆ m_eigenvaluethreshold

double Trk::MatrixTool::m_eigenvaluethreshold
private

cut on the minimum eigenvalue

Definition at line 154 of file MatrixTool.h.

◆ 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

std::vector<std::string> Trk::MatrixTool::m_inputHitmapFiles
private

input binary files containing the hitmaps

Definition at line 205 of file MatrixTool.h.

◆ m_inputMatrixFiles

std::vector<std::string> Trk::MatrixTool::m_inputMatrixFiles
private

input binary files containing matrix terms

Definition at line 202 of file MatrixTool.h.

◆ m_inputTFiles

std::vector<std::string> Trk::MatrixTool::m_inputTFiles
private

input binary files containing matrix terms

Definition at line 207 of file MatrixTool.h.

◆ m_inputVectorFiles

std::vector<std::string> Trk::MatrixTool::m_inputVectorFiles
private

input binary files containing vector terms

Definition at line 203 of file MatrixTool.h.

◆ m_logStream

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

logfile output stream

Definition at line 98 of file IMatrixTool.h.

◆ m_maxReadErrors

int Trk::MatrixTool::m_maxReadErrors
private

maximum number of reading TFile errors

Definition at line 212 of file MatrixTool.h.

◆ m_minNumHits

int Trk::MatrixTool::m_minNumHits
private

cut on the minimum number of hits per module

Definition at line 158 of file MatrixTool.h.

◆ m_minNumTrks

int Trk::MatrixTool::m_minNumTrks
private

cut on the minimum number of tracks per module

Definition at line 159 of file MatrixTool.h.

◆ m_modcut

int Trk::MatrixTool::m_modcut
private

cut on the weak modes which number is <par_modcut

Definition at line 157 of file MatrixTool.h.

◆ 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

std::string Trk::MatrixTool::m_pathbin
private

path binary files (in/out)

Definition at line 190 of file MatrixTool.h.

◆ m_pathtxt

std::string Trk::MatrixTool::m_pathtxt
private

path ascii files (in/out)

Definition at line 191 of file MatrixTool.h.

◆ m_prefixName

std::string Trk::MatrixTool::m_prefixName
private

prefix string to filenames

Definition at line 192 of file MatrixTool.h.

◆ m_pullcut

float Trk::MatrixTool::m_pullcut
private

pull cut for the automatic weak mode removal method

Definition at line 160 of file MatrixTool.h.

◆ m_readHitmaps

bool Trk::MatrixTool::m_readHitmaps
private

accumulate hitymap from files

Definition at line 174 of file MatrixTool.h.

◆ m_readTFiles

bool Trk::MatrixTool::m_readTFiles
private

write out files to a root file

Definition at line 177 of file MatrixTool.h.

◆ m_Remove_IBL_Rx

bool Trk::MatrixTool::m_Remove_IBL_Rx
private

Definition at line 228 of file MatrixTool.h.

◆ m_Remove_IBL_Ry

bool Trk::MatrixTool::m_Remove_IBL_Ry
private

Definition at line 229 of file MatrixTool.h.

◆ m_Remove_IBL_Rz

bool Trk::MatrixTool::m_Remove_IBL_Rz
private

Definition at line 230 of file MatrixTool.h.

◆ m_Remove_IBL_Tx

bool Trk::MatrixTool::m_Remove_IBL_Tx
private

Definition at line 225 of file MatrixTool.h.

◆ m_Remove_IBL_Ty

bool Trk::MatrixTool::m_Remove_IBL_Ty
private

Definition at line 226 of file MatrixTool.h.

◆ m_Remove_IBL_Tz

bool Trk::MatrixTool::m_Remove_IBL_Tz
private

Definition at line 227 of file MatrixTool.h.

◆ m_Remove_Pixel_Rx

bool Trk::MatrixTool::m_Remove_Pixel_Rx
private

Definition at line 221 of file MatrixTool.h.

◆ m_Remove_Pixel_Ry

bool Trk::MatrixTool::m_Remove_Pixel_Ry
private

Definition at line 222 of file MatrixTool.h.

◆ m_Remove_Pixel_Rz

bool Trk::MatrixTool::m_Remove_Pixel_Rz
private

Definition at line 223 of file MatrixTool.h.

◆ m_Remove_Pixel_Tx

bool Trk::MatrixTool::m_Remove_Pixel_Tx
private

Definition at line 218 of file MatrixTool.h.

◆ m_Remove_Pixel_Ty

bool Trk::MatrixTool::m_Remove_Pixel_Ty
private

Definition at line 219 of file MatrixTool.h.

◆ m_Remove_Pixel_Tz

bool Trk::MatrixTool::m_Remove_Pixel_Tz
private

Definition at line 220 of file MatrixTool.h.

◆ m_removeSpurious

bool Trk::MatrixTool::m_removeSpurious
private

run spurious removal

Definition at line 186 of file MatrixTool.h.

◆ m_runLocal

bool Trk::MatrixTool::m_runLocal
private

Run solving using Local method.

Definition at line 179 of file MatrixTool.h.

◆ m_scalaMatName

std::string Trk::MatrixTool::m_scalaMatName
private

Scalapack matrix name.

Definition at line 197 of file MatrixTool.h.

◆ m_scalaVecName

std::string Trk::MatrixTool::m_scalaVecName
private

Scalapack vector name.

Definition at line 198 of file MatrixTool.h.

◆ m_scale

double Trk::MatrixTool::m_scale
private

scale for big matrix and vector normalization

Definition at line 181 of file MatrixTool.h.

◆ m_scaleMatrix

bool Trk::MatrixTool::m_scaleMatrix
private

scale matrix by number of hits before solving

Definition at line 182 of file MatrixTool.h.

◆ m_softEigenmodeCut

double Trk::MatrixTool::m_softEigenmodeCut
private

add constant to diagonal to effectively cut on weak eigenmodes

Definition at line 184 of file MatrixTool.h.

◆ m_solveOption

int Trk::MatrixTool::m_solveOption
private

solving option

Definition at line 156 of file MatrixTool.h.

◆ m_tfileName

std::string Trk::MatrixTool::m_tfileName
private

prefix string to filenames

Definition at line 194 of file MatrixTool.h.

◆ m_useSparse

bool Trk::MatrixTool::m_useSparse
private

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

Definition at line 151 of file MatrixTool.h.

◆ 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

bool Trk::MatrixTool::m_writeEigenMat
private

write eigenvalues and eigenvectors into files ?

Definition at line 168 of file MatrixTool.h.

◆ m_writeEigenMatTxt

bool Trk::MatrixTool::m_writeEigenMatTxt
private

also write eigenvalues and eigenvectors into txt files ?

Definition at line 169 of file MatrixTool.h.

◆ m_writeHitmap

bool Trk::MatrixTool::m_writeHitmap
private

write hitmap into file

Definition at line 172 of file MatrixTool.h.

◆ m_writeHitmapTxt

bool Trk::MatrixTool::m_writeHitmapTxt
private

write hitmap into text file

Definition at line 173 of file MatrixTool.h.

◆ m_writeMat

bool Trk::MatrixTool::m_writeMat
private

write big matrix and vector into files ?

Definition at line 166 of file MatrixTool.h.

◆ m_writeMatTxt

bool Trk::MatrixTool::m_writeMatTxt
private

also write big matrix and vector into txt files ?

Definition at line 167 of file MatrixTool.h.

◆ m_writeModuleNames

bool Trk::MatrixTool::m_writeModuleNames
private

write module name instead of Identifier to vector file

Definition at line 170 of file MatrixTool.h.

◆ m_writeTFile

bool Trk::MatrixTool::m_writeTFile
private

write out files to a root file

Definition at line 176 of file MatrixTool.h.

◆ m_wSqMatrix

bool Trk::MatrixTool::m_wSqMatrix
private

write a triangular matrix by default (true: square format) ?

Definition at line 165 of file MatrixTool.h.

The documentation for this class was generated from the following files:
grepfile.info
info
Definition: grepfile.py:38
Trk::IMatrixTool::m_nHits
int m_nHits
Definition: IMatrixTool.h:100
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
Trk::MatrixTool::SOLVE_CLHEP
@ SOLVE_CLHEP
computation using CLHEP
Definition: MatrixTool.h:66
beamspotman.r
def r
Definition: beamspotman.py:676
Trk::MatrixTool::writeHitmap
void writeHitmap()
Definition: MatrixTool.cxx:2418
Trk::MatrixTool::m_alignModuleTool
ToolHandle< IAlignModuleTool > m_alignModuleTool
Definition: MatrixTool.h:142
Trk::MatrixTool::DIRECT_SOLVE_FAST
@ DIRECT_SOLVE_FAST
direct Fast (Eigen method) solving, already available matrix & vector
Definition: MatrixTool.h:63
PlotCalibFromCool.norm
norm
Definition: PlotCalibFromCool.py:100
Trk::MatrixTool::SOLVE_ROOT
@ SOLVE_ROOT
computation using ROOT
Definition: MatrixTool.h:65
pdg_comparison.sigma
sigma
Definition: pdg_comparison.py:324
Trk::AlSymMatBase::determinant
virtual double determinant()=0
TRTCalib_Extractor.hits
hits
Definition: TRTCalib_Extractor.py:35
Trk::MatrixTool::m_Remove_Pixel_Tz
bool m_Remove_Pixel_Tz
Definition: MatrixTool.h:220
Trk::z
@ z
global position (cartesian)
Definition: ParamDefs.h:57
Trk::MatrixTool::accumulateFromTFiles
bool accumulateFromTFiles()
Store Files in a tfile.
Definition: MatrixTool.cxx:957
Trk::MatrixTool::m_Remove_IBL_Tz
bool m_Remove_IBL_Tz
Definition: MatrixTool.h:227
python.Constants.FATAL
int FATAL
Definition: Control/AthenaCommon/python/Constants.py:19
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
FullCPAlgorithmsTest_eljob.flush
flush
Definition: FullCPAlgorithmsTest_eljob.py:186
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
Trk::MatrixTool::DIRECT_SOLVE_CLUSTER
@ DIRECT_SOLVE_CLUSTER
computation of alignment parameters from SCALAPAK already solved matrix
Definition: MatrixTool.h:64
Trk::AlignModuleList
std::vector< AlignModule * > AlignModuleList
Definition: AlignModuleList.h:37
Trk::MatrixTool::m_writeTFile
bool m_writeTFile
write out files to a root file
Definition: MatrixTool.h:176
Trk::MatrixTool::m_inputVectorFiles
std::vector< std::string > m_inputVectorFiles
input binary files containing vector terms
Definition: MatrixTool.h:203
Trk::MatrixTool::m_calDet
bool m_calDet
Compute bigmatrix's determinant ?
Definition: MatrixTool.h:164
Trk::IMatrixTool::nTracks
int nTracks() const
Definition: IMatrixTool.h:88
index
Definition: index.py:1
Trk::MatrixTool::solveCLHEP
int solveCLHEP()
Definition: MatrixTool.cxx:391
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
Trk::MatrixTool::DIRECT_SOLVE
@ DIRECT_SOLVE
direct solving (LAPACK), already available matrix & vector
Definition: MatrixTool.h:62
Trk::MatrixTool::m_eigenvaluethreshold
double m_eigenvaluethreshold
cut on the minimum eigenvalue
Definition: MatrixTool.h:154
Trk::MatrixTool::SOLVE_FAST
@ SOLVE_FAST
Fast (Eigen method) solving after data accumulation.
Definition: MatrixTool.h:61
Trk::MatrixTool::m_scale
double m_scale
scale for big matrix and vector normalization
Definition: MatrixTool.h:181
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
Trk::MatrixTool::m_tfileName
std::string m_tfileName
prefix string to filenames
Definition: MatrixTool.h:194
Trk::MatrixTool::printModuleSolution
void printModuleSolution(std::ostream &os, const AlignModule *module, const CLHEP::HepSymMatrix *cov) const
namespace { class RestoreIOSFlags { public: RestoreIOSFlags (std::ostream &os) : m_os(&os),...
Definition: MatrixTool.cxx:1706
Trk::MatrixTool::m_Remove_Pixel_Rx
bool m_Remove_Pixel_Rx
Definition: MatrixTool.h:221
plotBeamSpotVxVal.cov
cov
Definition: plotBeamSpotVxVal.py:201
Trk::AlSymMatBase::SetPathTxt
virtual void SetPathTxt(const std::string &)=0
Trk::MatrixTool::m_writeMatTxt
bool m_writeMatTxt
also write big matrix and vector into txt files ?
Definition: MatrixTool.h:167
Trk::MatrixTool::m_maxReadErrors
int m_maxReadErrors
maximum number of reading TFile errors
Definition: MatrixTool.h:212
AthCommonDataStore< AthCommonMsg< AlgTool > >::m_evtStore
StoreGateSvc_t m_evtStore
Pointer to StoreGate (event store by default)
Definition: AthCommonDataStore.h:390
AthCommonDataStore< AthCommonMsg< AlgTool > >::m_vhka
std::vector< SG::VarHandleKeyArray * > m_vhka
Definition: AthCommonDataStore.h:398
PixelModuleFeMask_create_db.nModules
nModules
Definition: PixelModuleFeMask_create_db.py:47
AthCommonMsg< AlgTool >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition: AthCommonMsg.h:30
Trk::IMatrixTool::IMatrixTool
IMatrixTool()
constructor
Definition: IMatrixTool.h:114
Trk::MatrixTool::m_removeSpurious
bool m_removeSpurious
run spurious removal
Definition: MatrixTool.h:186
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
Trk::MatrixTool::m_Remove_IBL_Tx
bool m_Remove_IBL_Tx
Definition: MatrixTool.h:225
Trk::MatrixTool::m_useSparse
bool m_useSparse
flag to use AlSpaMat for the big matrix (default is AlSymMat)
Definition: MatrixTool.h:151
yodamerge_tmp.scale
scale
Definition: yodamerge_tmp.py:138
Trk::MatrixTool::m_readHitmaps
bool m_readHitmaps
accumulate hitymap from files
Definition: MatrixTool.h:174
Trk::MatrixTool::solveSparseEigen
int solveSparseEigen()
Definition: MatrixTool.cxx:2284
atlasStyleMacro.icol
int icol
Definition: atlasStyleMacro.py:13
Trk::AlVec::WritePartial
StatusCode WritePartial(const std::string &, bool, double, std::map< int, unsigned long long >, float)
Definition: AlVec.cxx:366
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
Trk::MatrixTool::m_Remove_Pixel_Ry
bool m_Remove_Pixel_Ry
Definition: MatrixTool.h:222
SG::VarHandleKeyArray::setOwner
virtual void setOwner(IDataHandleHolder *o)=0
Trk::MatrixTool::m_scalaVecName
std::string m_scalaVecName
Scalapack vector name.
Definition: MatrixTool.h:198
Trk::AlSymMatBase::ptrMap
const datamap * ptrMap() const
Definition: AlSymMatBase.h:149
TMatrixDSparse
class TMatrixTSparse< double > TMatrixDSparse
Definition: AlSpaMat.h:14
IDTPMcnv.htype
htype
Definition: IDTPMcnv.py:27
Trk::MatrixTool::m_calculateFullCovariance
bool m_calculateFullCovariance
Definition: MatrixTool.h:188
Trk::MatrixTool::m_softEigenmodeCut
double m_softEigenmodeCut
add constant to diagonal to effectively cut on weak eigenmodes
Definition: MatrixTool.h:184
Trk::MatrixTool::m_Remove_Pixel_Tx
bool m_Remove_Pixel_Tx
Definition: MatrixTool.h:218
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
AlignModule
AlignModule is a grouping of TrkDetElementBase objects, grouped according to the type of alignment,...
python.PyAthena.module
module
Definition: PyAthena.py:131
Trk::MatrixTool::m_writeEigenMatTxt
bool m_writeEigenMatTxt
also write eigenvalues and eigenvectors into txt files ?
Definition: MatrixTool.h:169
Trk::MatrixTool::m_writeMat
bool m_writeMat
write big matrix and vector into files ?
Definition: MatrixTool.h:166
Trk::AlSymMatBase::makeTMatrix
virtual TMatrixDSparse * makeTMatrix()=0
Trk::MatrixTool::m_modcut
int m_modcut
cut on the weak modes which number is <par_modcut
Definition: MatrixTool.h:157
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
Trk::MatrixTool::m_minNumHits
int m_minNumHits
cut on the minimum number of hits per module
Definition: MatrixTool.h:158
ev
int ev
Definition: globals.cxx:25
AthCommonDataStore
Definition: AthCommonDataStore.h:52
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
CheckAppliedSFs.e3
e3
Definition: CheckAppliedSFs.py:264
Trk::MatrixTool::accumulateFromBinaries
bool accumulateFromBinaries()
accumulates derivates from binary files
Definition: MatrixTool.cxx:763
dqt_zlumi_pandas.err
err
Definition: dqt_zlumi_pandas.py:182
Trk::MatrixTool::m_writeEigenMat
bool m_writeEigenMat
write eigenvalues and eigenvectors into files ?
Definition: MatrixTool.h:168
lumiFormat.i
int i
Definition: lumiFormat.py:85
Trk::AlVec::size
int size() const
Definition: AlVec.h:109
beamspotman.n
n
Definition: beamspotman.py:731
Trk::MatrixTool::m_readTFiles
bool m_readTFiles
write out files to a root file
Definition: MatrixTool.h:177
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
Trk::MatrixTool::m_bigvector
AlVec * m_bigvector
vector to contain first derivative terms to be used for alignment
Definition: MatrixTool.h:148
vector
Definition: MultiHisto.h:13
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
Trk::MatrixTool::m_inputHitmapFiles
std::vector< std::string > m_inputHitmapFiles
input binary files containing the hitmaps
Definition: MatrixTool.h:205
Trk::MatrixTool::SOLVE
@ SOLVE
solving after data accumulation (LAPACK)
Definition: MatrixTool.h:60
Trk::MatrixTool::postSolvingLapack
void postSolvingLapack(AlVec *dChi2, AlSymMat *d2Chi2, AlVec &w, AlMat &z, int size)
Definition: MatrixTool.cxx:1958
Trk::MatrixTool::m_bigmatrix
AlSymMatBase * m_bigmatrix
matrix to contain second derivative terms to be used for alignment
Definition: MatrixTool.h:145
Trk::MatrixTool::solveLocal
int solveLocal()
Definition: MatrixTool.cxx:633
Trk::MatrixTool::m_AlignPixelbutNotIBL
bool m_AlignPixelbutNotIBL
Definition: MatrixTool.h:215
Trk::MatrixTool::storeInTFile
void storeInTFile(const TString &filename)
Store Files in a tfile.
Definition: MatrixTool.cxx:862
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:522
test_pyathena.parent
parent
Definition: test_pyathena.py:15
Trk::MatrixTool::m_writeHitmapTxt
bool m_writeHitmapTxt
write hitmap into text file
Definition: MatrixTool.h:173
xAOD::uint64_t
uint64_t
Definition: EventInfo_v1.cxx:123
Trk::AlSymMatBase::Write
virtual StatusCode Write(const std::string &, bool, bool, double, float)=0
Trk::MatrixTool::m_scaleMatrix
bool m_scaleMatrix
scale matrix by number of hits before solving
Definition: MatrixTool.h:182
Trk::MatrixTool::printGlobalSolution
void printGlobalSolution(std::ostream &os, const CLHEP::HepSymMatrix *cov)
Definition: MatrixTool.cxx:1621
AthCommonDataStore< AthCommonMsg< AlgTool > >::m_detStore
StoreGateSvc_t m_detStore
Pointer to StoreGate (detector store by default)
Definition: AthCommonDataStore.h:393
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
Scale
void Scale(TH1 *h, double d=1)
Definition: comparitor.cxx:76
DataVector< AlignPar >
Trk::MatrixTool::m_Align_db_step
float m_Align_db_step
corr in the diagonal basis step for the third pass in the auto weak mode removal method
Definition: MatrixTool.h:162
Trk::MatrixTool::solveROOT
int solveROOT()
Definition: MatrixTool.cxx:265
dot.dot
def dot(G, fn, nodesToHighlight=[])
Definition: dot.py:5
ReadFromCoolCompare.os
os
Definition: ReadFromCoolCompare.py:231
Trk::MatrixTool::m_wSqMatrix
bool m_wSqMatrix
write a triangular matrix by default (true: square format) ?
Definition: MatrixTool.h:165
AthAlgTool::AthAlgTool
AthAlgTool()
Default constructor:
Trk::MatrixTool::m_activeIndices
std::vector< int > m_activeIndices
vector of indices which pass the min-hits cut
Definition: MatrixTool.h:209
SG::VarHandleKeyArray::renounce
virtual void renounce()=0
SG::HandleClassifier::type
std::conditional< std::is_base_of< SG::VarHandleKeyArray, T >::value, VarHandleKeyArrayType, type2 >::type type
Definition: HandleClassifier.h:54
Trk::IMatrixTool::m_nentries
int m_nentries
Definition: IMatrixTool.h:106
Trk::MatrixTool::m_aNDoF
int m_aNDoF
number of active DoF (size of m_activeIndices)
Definition: MatrixTool.h:210
merge_scale_histograms.doc
string doc
Definition: merge_scale_histograms.py:9
Trk::MatrixTool::m_Remove_IBL_Ry
bool m_Remove_IBL_Ry
Definition: MatrixTool.h:229
Trk::IMatrixTool::m_logStream
std::ostream * m_logStream
logfile output stream
Definition: IMatrixTool.h:98
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:221
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:29
AlignPar
AlignPar contains all the information related to an alignment parameter of a particular align module ...
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
Trk::IMatrixTool::m_alignModuleMap
std::map< int, int > m_alignModuleMap
Definition: IMatrixTool.h:105
Trk::MatrixTool::fillVecMods
static int fillVecMods()
Definition: MatrixTool.cxx:1759
Trk::MatrixTool::m_Remove_IBL_Rz
bool m_Remove_IBL_Rz
Definition: MatrixTool.h:230
Trk::IMatrixTool::m_nMeasurements
int m_nMeasurements
Definition: IMatrixTool.h:102
Trk::MatrixTool::readHitmaps
void readHitmaps()
Definition: MatrixTool.cxx:2455
Trk::MatrixTool::m_runLocal
bool m_runLocal
Run solving using Local method.
Definition: MatrixTool.h:179
Trk::IMatrixTool::m_nTracks
int m_nTracks
Definition: IMatrixTool.h:101
Trk::MatrixTool::m_AlignIBLbutNotPixel
bool m_AlignIBLbutNotPixel
Definition: MatrixTool.h:214
DataVector::end
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
Trk::IMatrixTool::nHits
int nHits() const
Definition: IMatrixTool.h:84
Trk::MatrixTool::NONE
@ NONE
not solve in any case (to be used when ipc)
Definition: MatrixTool.h:59
DeMoScan.index
string index
Definition: DeMoScan.py:364
python.testIfMatch.matrix
matrix
Definition: testIfMatch.py:66
Trk::MatrixTool::m_prefixName
std::string m_prefixName
prefix string to filenames
Definition: MatrixTool.h:192
Trk::MatrixTool::m_writeHitmap
bool m_writeHitmap
write hitmap into file
Definition: MatrixTool.h:172
a
TList * a
Definition: liststreamerinfos.cxx:10
Trk::AlSymMatBase::SetPathBin
virtual void SetPathBin(const std::string &)=0
h
Trk::MatrixTool::m_solveOption
int m_solveOption
solving option
Definition: MatrixTool.h:156
Trk::MatrixTool::m_DeactivateSCT_ECA_LastDisk
bool m_DeactivateSCT_ECA_LastDisk
Definition: MatrixTool.h:216
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
copySelective.target
string target
Definition: copySelective.py:37
Pythia8_RapidityOrderMPI.val
val
Definition: Pythia8_RapidityOrderMPI.py:14
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
DEBUG
#define DEBUG
Definition: page_access.h:11
Trk::MatrixTool::m_Remove_Pixel_Ty
bool m_Remove_Pixel_Ty
Definition: MatrixTool.h:219
AthCommonMsg< AlgTool >::msg
MsgStream & msg() const
Definition: AthCommonMsg.h:24
Trk::MatrixTool::m_pathbin
std::string m_pathbin
path binary files (in/out)
Definition: MatrixTool.h:190
Trk::MatrixTool::m_pathtxt
std::string m_pathtxt
path ascii files (in/out)
Definition: MatrixTool.h:191
CaloCellTimeCorrFiller.filename
filename
Definition: CaloCellTimeCorrFiller.py:24
Trk::AlVec::SetPathBin
void SetPathBin(const std::string &)
Definition: AlVec.cxx:303
SG::VarHandleBase::vhKey
SG::VarHandleKey & vhKey()
Return a non-const reference to the HandleKey.
Definition: StoreGate/src/VarHandleBase.cxx:623
Trk::MatrixTool::m_Remove_IBL_Ty
bool m_Remove_IBL_Ty
Definition: MatrixTool.h:226
Trk::MatrixTool::m_minNumTrks
int m_minNumTrks
cut on the minimum number of tracks per module
Definition: MatrixTool.h:159
Trk::AlVec::SetPathTxt
void SetPathTxt(const std::string &)
Definition: AlVec.cxx:309
Trk::MatrixTool::m_scalaMatName
std::string m_scalaMatName
Scalapack matrix name.
Definition: MatrixTool.h:197
Trk::MatrixTool::m_writeModuleNames
bool m_writeModuleNames
write module name instead of Identifier to vector file
Definition: MatrixTool.h:170
copySelective.source
string source
Definition: copySelective.py:32
python.Bindings.keys
keys
Definition: Control/AthenaPython/python/Bindings.py:798
Trk::MatrixTool::m_Remove_IBL_Rx
bool m_Remove_IBL_Rx
Definition: MatrixTool.h:228
merge.status
status
Definition: merge.py:17
Trk::MatrixTool::m_inputMatrixFiles
std::vector< std::string > m_inputMatrixFiles
input binary files containing matrix terms
Definition: MatrixTool.h:202
Trk::MatrixTool::m_eigenvalueStep
float m_eigenvalueStep
eigenvalue step for the second pass in the automatic weak mode removal method
Definition: MatrixTool.h:161
LArCellNtuple.ifile
string ifile
Definition: LArCellNtuple.py:133
Trk::MatrixTool::solveLapack
int solveLapack()
Definition: MatrixTool.cxx:1765
python.Constants.VERBOSE
int VERBOSE
Definition: Control/AthenaCommon/python/Constants.py:14
DataVector::at
const T * at(size_type n) const
Access an element, as an rvalue.
Trk::MatrixTool::m_diagonalize
bool m_diagonalize
run diagonalization instead of inversion
Definition: MatrixTool.h:153
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:200
Trk::MatrixTool::m_Remove_Pixel_Rz
bool m_Remove_Pixel_Rz
Definition: MatrixTool.h:223
Trk::MatrixTool::m_pullcut
float m_pullcut
pull cut for the automatic weak mode removal method
Definition: MatrixTool.h:160
CxxUtils::ivec
vec_fb< typename boost::int_t< sizeof(T) *8 >::exact, N > ivec
Definition: vec_fb.h:53
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
python.compressB64.c
def c
Definition: compressB64.py:93
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
AthCommonDataStore::declareGaudiProperty
Gaudi::Details::PropertyBase & declareGaudiProperty(Gaudi::Property< T > &hndl, const SG::VarHandleKeyType &)
specialization for handling Gaudi::Property<SG::VarHandleKey>
Definition: AthCommonDataStore.h:156
DataVector::empty
bool empty() const noexcept
Returns true if the collection is empty.
extractSporadic.myFile
myFile
Definition: extractSporadic.py:87
Trk::MatrixTool::m_inputTFiles
std::vector< std::string > m_inputTFiles
input binary files containing matrix terms
Definition: MatrixTool.h:207
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
WriteBchToCool.moduleList
moduleList
Definition: WriteBchToCool.py:72
fitman.k
k
Definition: fitman.py:528
CscCalibQuery.nFiles
int nFiles
Definition: CscCalibQuery.py:332
Generated on Thu Nov 7 2024 22:28:45 for ATLAS Offline Software by   doxygen 1.8.18