ATLAS Offline Software
MatrixTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // MatrixTool.cxx
6 
7 // AlgTool for creating big matrix and vector, manipulation of entries, and
8 // solving for alignment parameters
9 // Robert Harrington, started 1/5/08 based on SiGlobalChi2Align
10 
11 #include "GaudiKernel/StatusCode.h"
12 #include "GaudiKernel/MsgStream.h"
13 
14 #include "GaudiKernel/AlgTool.h"
15 
16 #include "TrkAlgebraUtils/AlVec.h"
17 #include "TrkAlgebraUtils/AlSymMatBase.h"
18 #include "TrkAlgebraUtils/AlMat.h"
19 #include "TrkAlgebraUtils/AlSymMat.h"
20 #include "TrkAlgebraUtils/AlSpaMat.h"
21 
22 #include "TrkAlignEvent/AlignModule.h"
23 #include "TrkAlignEvent/AlignModuleList.h"
24 #include "TrkAlignEvent/AlignPar.h"
25 #include "TrkAlignGenTools/MatrixTool.h"
26 #include "TrkAlignInterfaces/IAlignModuleTool.h"
27 
28 #include "CLHEP/Matrix/Matrix.h"
29 #include "CLHEP/Matrix/SymMatrix.h"
30 #include "CLHEP/Matrix/Vector.h"
31 
32 #include <TMatrixT.h>
33 #include <TMatrixDSym.h>
34 #include <TMatrixDSparse.h>
35 #include <TVectorD.h>
36 #include <TVectorT.h>
37 #include <TString.h>
38 #include <TFile.h>
39 
40 
41 #include <TDecompBK.h>
42 
43 #include <cmath>
44 #include <ctime>
45 #include <algorithm>
46 #include <iterator>
47 
48 //restore ostream
49 #include <boost/io/ios_state.hpp>
50 
51 
52 
53 #include <sys/resource.h>
54 
55 namespace Trk {
56 
57  //_______________________________________________________________________
58  MatrixTool::MatrixTool(const std::string& type, const std::string& name,
59  const IInterface* parent)
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  }
193 
194  //_______________________________________________________________________
195  MatrixTool::~MatrixTool()
196  {
197  delete m_bigmatrix;
198  delete m_bigvector;
199  }
200 
201  //_______________________________________________________________________
202  StatusCode MatrixTool::initialize()
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  }
221 
222  //_______________________________________________________________________
223  StatusCode MatrixTool::finalize()
224  {
225  ATH_MSG_DEBUG("finalize() of MatrixTool");
226 
227  return StatusCode::SUCCESS;
228  }
229 
230 
231  //_______________________________________________________________________
232  StatusCode MatrixTool::allocateMatrix(int nDoF)
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  }
258 
259  //_______________________________________________________________________
260  void MatrixTool::prepareBinaryFiles(int)
261  {
262  }
263 
264  //_______________________________________________________________________
265  int MatrixTool::solveROOT()
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  }
389 
390  //_______________________________________________________________________
391  int MatrixTool::solveCLHEP()
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  }
631 
632  //_______________________________________________________________________
633  int MatrixTool::solveLocal()
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  }
748 
749  //________________________________________________________________________
750  bool MatrixTool::accumulateFromFiles()
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  }
761 
762  //________________________________________________________________________
763  bool MatrixTool::accumulateFromBinaries()
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  }
860 
861  //_______________________________________________________________________
862  void MatrixTool::storeInTFile(const TString& filename)
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  }
955 
956 //________________________________________________________________________
957  bool MatrixTool::accumulateFromTFiles()
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  }
1263 
1264 
1265 
1266 
1267  //_______________________________________________________________________
1268  int MatrixTool::solve()
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  }
1586 
1587 
1588  //_______________________________________________________________________
1589  void MatrixTool::addFirstDerivatives(AlVec* )
1590  {
1591  }
1592 
1593  //_______________________________________________________________________
1594  void MatrixTool::addSecondDerivatives(AlSymMatBase* )
1595  {
1596  }
1597 
1598  //________________________________________________________________________
1599  void MatrixTool::addFirstDerivatives(std::list<int,double>& )
1600  {
1601  }
1602 
1603  //________________________________________________________________________
1604  void MatrixTool::addSecondDerivatives(std::list<std::pair<int,int>,double >&)
1605  {
1606  }
1607 
1608  //________________________________________________________________________
1609  void MatrixTool::addFirstDerivative(int irow, double firstderiv)
1610  {
1611  (*m_bigvector)[irow] += firstderiv;
1612  }
1613 
1614  //________________________________________________________________________
1615  void MatrixTool::addSecondDerivative(int irow, int icol, double secondderiv)
1616  {
1617  (*m_bigmatrix)[irow][icol] += secondderiv;
1618  }
1619 
1620  //________________________________________________________________________
1621  void MatrixTool::printGlobalSolution(std::ostream & os, const CLHEP::HepSymMatrix * cov)
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  }
1666 
1667  //________________________________________________________________________
1668  void MatrixTool::printGlobalSolution(std::ostream & os, const TMatrixDSym * cov0)
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  }
1705  //________________________________________________________________________
1706  void MatrixTool::printModuleSolution(std::ostream & os, const AlignModule * module, const CLHEP::HepSymMatrix * cov) const
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  }
1757 
1758  //________________________________________________________________________
1759  int MatrixTool::fillVecMods()
1760  {
1761  return 0;
1762  }
1763 
1764  //________________________________________________________________________
1765  int MatrixTool::solveLapack( )
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  }
1864 
1865  //________________________________________________________________________
1866  StatusCode MatrixTool::spuriousRemoval()
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  }
1956 
1957  //________________________________________________________________________
1958  void MatrixTool::postSolvingLapack(AlVec * dChi2, AlSymMat * d2Chi2, AlVec &w, AlMat &z, int size)
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  }
2282 
2283  //________________________________________________________________________
2284  int MatrixTool::solveSparseEigen()
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  }
2416 
2417  //________________________________________________________________________
2418  void MatrixTool::writeHitmap()
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  }
2453 
2454  //________________________________________________________________________
2455  void MatrixTool::readHitmaps()
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  }
2503 
2504 } // end of namespace
2505 
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::IMatrixTool
Definition: IMatrixTool.h:35
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
Trk::AlignPar::sigma
double sigma() const
returns sigma
Definition: AlignPar.h:65
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
Trk::MatrixTool::solve
int solve()
solves for alignment parameters
Definition: MatrixTool.cxx:1268
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::AlignPar::initPar
double initPar() const
returns initial parameter and error
Definition: AlignPar.h:53
Trk::MatrixTool::m_writeTFile
bool m_writeTFile
write out files to a root file
Definition: MatrixTool.h:176
Trk::AlVec
Definition: AlVec.h:23
Trk::MatrixTool::m_inputVectorFiles
std::vector< std::string > m_inputVectorFiles
input binary files containing vector terms
Definition: MatrixTool.h:203
AlSpaMat.h
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::AlSymMatBase
Definition: AlSymMatBase.h:38
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::accumulateFromFiles
bool accumulateFromFiles()
accumulates derivates from files.
Definition: MatrixTool.cxx:750
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
MatrixTool.h
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
Trk::AlMat::ReadScalaPack
StatusCode ReadScalaPack(const std::string &)
Definition: AlMat.cxx:364
PixelModuleFeMask_create_db.nModules
nModules
Definition: PixelModuleFeMask_create_db.py:47
Trk::AlignModule
Definition: AlignModule.h:45
AlVec.h
AthCommonMsg< AlgTool >::msgLvl
bool msgLvl(const MSG::Level lvl) const
Definition: AthCommonMsg.h:30
Trk::MatrixTool::m_removeSpurious
bool m_removeSpurious
run spurious removal
Definition: MatrixTool.h:186
Trk::AlMat::SetPathTxt
void SetPathTxt(const std::string &)
Definition: AlMat.cxx:403
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
AlSymMat.h
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::AlMat::ncol
long int ncol() const
Definition: AlMat.h:162
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
Trk::AlSymMat
Definition: AlSymMat.h:26
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
Trk::MatrixTool::m_Remove_Pixel_Ry
bool m_Remove_Pixel_Ry
Definition: MatrixTool.h:222
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
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
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::AlVec::norm
double norm() const
Definition: AlVec.cxx:206
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
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
Trk::AlignPar::setPar
void setPar(double par, double err)
sets final parameter and error
Definition: AlignPar.h:71
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::AlMat::SetPathBin
void SetPathBin(const std::string &)
Definition: AlMat.cxx:397
Trk::MatrixTool::solveLocal
int solveLocal()
Definition: MatrixTool.cxx:633
Trk::MatrixTool::m_AlignPixelbutNotIBL
bool m_AlignPixelbutNotIBL
Definition: MatrixTool.h:215
Trk::MatrixTool::spuriousRemoval
StatusCode spuriousRemoval()
Definition: MatrixTool.cxx:1866
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
Trk::AlSymMat::Read
virtual StatusCode Read(const std::string &, int &, bool &, float &) override final
Definition: AlSymMat.cxx:671
Trk::AlignPar
Definition: AlignPar.h:25
AlMat.h
test_pyathena.parent
parent
Definition: test_pyathena.py:15
Trk::MatrixTool::initialize
StatusCode initialize()
initialize
Definition: MatrixTool.cxx:202
Trk::AlMat
Definition: AlMat.h:27
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::MatrixTool::addSecondDerivative
void addSecondDerivative(int irow, int icol, double secondderiv)
Definition: MatrixTool.cxx:1615
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
histSizes.list
def list(name, path='/')
Definition: histSizes.py:38
Trk::MatrixTool::printGlobalSolution
void printGlobalSolution(std::ostream &os, const CLHEP::HepSymMatrix *cov)
Definition: MatrixTool.cxx:1621
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
Scale
void Scale(TH1 *h, double d=1)
Definition: comparitor.cxx:76
DataVector< AlignPar >
Trk::MatrixTool::addSecondDerivatives
void addSecondDerivatives(AlSymMatBase *matrix)
adds second derivatives to matrix
Definition: MatrixTool.cxx:1594
AlignModuleList.h
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::AlMat::nrow
long int nrow() const
Definition: AlMat.h:157
Trk::MatrixTool::m_wSqMatrix
bool m_wSqMatrix
write a triangular matrix by default (true: square format) ?
Definition: MatrixTool.h:165
Trk::MatrixTool::addFirstDerivatives
void addFirstDerivatives(AlVec *vector)
adds first derivative to vector
Definition: MatrixTool.cxx:1589
Trk::MatrixTool::~MatrixTool
virtual ~MatrixTool()
Virtual destructor.
Definition: MatrixTool.cxx:195
Trk::MatrixTool::m_activeIndices
std::vector< int > m_activeIndices
vector of indices which pass the min-hits cut
Definition: MatrixTool.h:209
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
Trk::MatrixTool::m_aNDoF
int m_aNDoF
number of active DoF (size of m_activeIndices)
Definition: MatrixTool.h:210
Trk::AlMat::Write
StatusCode Write(const std::string &, bool, unsigned int precision=6)
Definition: AlMat.cxx:409
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
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
AlSymMatBase.h
AlignModule.h
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::MatrixTool::allocateMatrix
StatusCode allocateMatrix(int nDoF=0)
allocates memory for big matrix and big vector
Definition: MatrixTool.cxx:232
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
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
Trk::AlVec::ReadPartial
StatusCode ReadPartial(const std::string &, double &, std::map< int, unsigned long long > &, float &)
Definition: AlVec.cxx:484
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
IAlignModuleTool.h
CaloCellTimeCorrFiller.filename
filename
Definition: CaloCellTimeCorrFiller.py:24
Trk::AlVec::SetPathBin
void SetPathBin(const std::string &)
Definition: AlVec.cxx:303
Trk::MatrixTool::m_Remove_IBL_Ty
bool m_Remove_IBL_Ty
Definition: MatrixTool.h:226
Trk::MatrixTool::prepareBinaryFiles
void prepareBinaryFiles(int solveOption)
reads/writes matrix entries from/to binary files as necessary
Definition: MatrixTool.cxx:260
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
AlignPar.h
Trk::MatrixTool::finalize
StatusCode finalize()
initialize
Definition: MatrixTool.cxx:223
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
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
AthAlgTool
Definition: AthAlgTool.h:26
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
Trk::MatrixTool::MatrixTool
MatrixTool(const std::string &type, const std::string &name, const IInterface *parent)
Constructor.
Definition: MatrixTool.cxx:58
CxxUtils::ivec
vec_fb< typename boost::int_t< sizeof(T) *8 >::exact, N > ivec
Definition: vec_fb.h:53
Trk::MatrixTool::addFirstDerivative
void addFirstDerivative(int irow, double firstderiv)
Definition: MatrixTool.cxx:1609
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.
DataVector::empty
bool empty() const noexcept
Returns true if the collection is empty.
Trk::AlSpaMat
Definition: AlSpaMat.h:27
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
Trk::AlSpaMat::Read
virtual StatusCode Read(const std::string &, int &, bool &, float &) override final
Definition: AlSpaMat.cxx:737
CscCalibQuery.nFiles
int nFiles
Definition: CscCalibQuery.py:332
Generated on Thu Nov 7 2024 21:20:29 for ATLAS Offline Software by   doxygen 1.8.18