|
ATLAS Offline Software
|
#include <MatrixTool.h>
|
| MatrixTool (const std::string &type, const std::string &name, const IInterface *parent) |
| Constructor. More...
|
|
virtual | ~MatrixTool () |
| Virtual destructor. More...
|
|
StatusCode | initialize () |
| initialize More...
|
|
StatusCode | finalize () |
| initialize More...
|
|
StatusCode | allocateMatrix (int nDoF=0) |
| allocates memory for big matrix and big vector More...
|
|
void | prepareBinaryFiles (int solveOption) |
| reads/writes matrix entries from/to binary files as necessary More...
|
|
void | addFirstDerivatives (AlVec *vector) |
| adds first derivative to vector More...
|
|
void | addFirstDerivatives (std::list< int, double > &derivatives) |
| adds first derivative to vector for only some entries More...
|
|
void | addFirstDerivative (int irow, double firstderiv) |
|
void | addSecondDerivatives (AlSymMatBase *matrix) |
| adds second derivatives to matrix More...
|
|
void | addSecondDerivatives (std::list< std::pair< int, int >, double > &derivatives) |
| adds first derivative to vector for only some entries More...
|
|
void | addSecondDerivative (int irow, int icol, double secondderiv) |
|
bool | accumulateFromFiles () |
| accumulates derivates from files. More...
|
|
bool | accumulateFromBinaries () |
| accumulates derivates from binary files More...
|
|
int | solve () |
| solves for alignment parameters More...
|
|
void | storeInTFile (const TString &filename) |
| Store Files in a tfile. More...
|
|
bool | accumulateFromTFiles () |
| Store Files in a tfile. More...
|
|
void | printModuleSolution (std::ostream &os, const AlignModule *module, const CLHEP::HepSymMatrix *cov) const |
| namespace { class RestoreIOSFlags { public: RestoreIOSFlags (std::ostream &os) : m_os(&os), m_flags(m_os->flags()), m_precision(m_os->precision()) {} ~RestoreIOSFlags() { m_os->flags(m_flags); m_os->precision(m_precision); } private: std::ostream *m_os; std::ios_base::fmtflags m_flags; std::streamsize m_precision; }; } More...
|
|
void | printGlobalSolution (std::ostream &os, const CLHEP::HepSymMatrix *cov) |
|
void | printGlobalSolution (std::ostream &os, const TMatrixDSym *cov) |
|
ServiceHandle< StoreGateSvc > & | evtStore () |
| The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc . More...
|
|
const ServiceHandle< StoreGateSvc > & | evtStore () const |
| The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc . More...
|
|
const ServiceHandle< StoreGateSvc > & | detStore () const |
| The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc . More...
|
|
virtual StatusCode | sysInitialize () override |
| Perform system initialization for an algorithm. More...
|
|
virtual StatusCode | sysStart () override |
| Handle START transition. More...
|
|
virtual std::vector< Gaudi::DataHandle * > | inputHandles () const override |
| Return this algorithm's input handles. More...
|
|
virtual std::vector< Gaudi::DataHandle * > | outputHandles () const override |
| Return this algorithm's output handles. More...
|
|
Gaudi::Details::PropertyBase & | declareProperty (Gaudi::Property< T > &t) |
|
Gaudi::Details::PropertyBase * | declareProperty (const std::string &name, SG::VarHandleKey &hndl, const std::string &doc, const SG::VarHandleKeyType &) |
| Declare a new Gaudi property. More...
|
|
Gaudi::Details::PropertyBase * | declareProperty (const std::string &name, SG::VarHandleBase &hndl, const std::string &doc, const SG::VarHandleType &) |
| Declare a new Gaudi property. More...
|
|
Gaudi::Details::PropertyBase * | declareProperty (const std::string &name, SG::VarHandleKeyArray &hndArr, const std::string &doc, const SG::VarHandleKeyArrayType &) |
|
Gaudi::Details::PropertyBase * | declareProperty (const std::string &name, T &property, const std::string &doc, const SG::NotHandleType &) |
| Declare a new Gaudi property. More...
|
|
Gaudi::Details::PropertyBase * | declareProperty (const std::string &name, T &property, const std::string &doc="none") |
| Declare a new Gaudi property. More...
|
|
void | updateVHKA (Gaudi::Details::PropertyBase &) |
|
MsgStream & | msg () const |
|
MsgStream & | msg (const MSG::Level lvl) const |
|
bool | msgLvl (const MSG::Level lvl) const |
|
void | addModule (int alignModuleIndex, int nAlignParam) |
|
int | entryNumber (int alignModuleIndex) |
|
void | setNHits (int n) |
| set module identifier More...
|
|
int | nHits () const |
|
void | setNTracks (int n) |
| set number of tracks More...
|
|
int | nTracks () const |
|
void | setNMeasurements (int n) |
| set number of measurements More...
|
|
int | nMeasurements () const |
|
virtual void | setLogStream (std::ostream *os) |
| sets the output stream for the logfile More...
|
|
Definition at line 55 of file MatrixTool.h.
◆ StoreGateSvc_t
◆ SolveOption
Enumerator |
---|
NONE | not solve in any case (to be used when ipc)
|
SOLVE | solving after data accumulation (LAPACK)
|
SOLVE_FAST | Fast (Eigen method) solving after data accumulation.
|
DIRECT_SOLVE | direct solving (LAPACK), already available matrix & vector
|
DIRECT_SOLVE_FAST | direct Fast (Eigen method) solving, already available matrix & vector
|
DIRECT_SOLVE_CLUSTER | computation of alignment parameters from SCALAPAK already solved matrix
|
SOLVE_ROOT | computation using ROOT
|
SOLVE_CLHEP | computation using CLHEP
|
Definition at line 58 of file MatrixTool.h.
◆ MatrixTool()
MatrixTool::MatrixTool |
( |
const std::string & |
type, |
|
|
const std::string & |
name, |
|
|
const IInterface * |
parent |
|
) |
| |
Constructor.
Definition at line 58 of file MatrixTool.cxx.
104 declareInterface<IMatrixTool>(
this);
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");
164 std::vector<std::string> defaultHitmapInput;
165 defaultMatInput.emplace_back(
"hitmap.bin");
◆ ~MatrixTool()
MatrixTool::~MatrixTool |
( |
| ) |
|
|
virtual |
◆ accumulateFromBinaries()
bool MatrixTool::accumulateFromBinaries |
( |
| ) |
|
accumulates derivates from binary files
Definition at line 763 of file MatrixTool.cxx.
768 int nDoF=alignParList->
size();
770 std::map<int,unsigned long long> modIndexMap;
771 float dummyVersion(0.);
772 double totalscale=0.;
777 AlVec newVector(nDoF);
778 std::map<int,unsigned long long> newModIndexMap;
784 if (
sc==StatusCode::FAILURE) {
789 msg(
MSG::FATAL) <<
"vector wrong size! newVector size "<<newVector.size()
796 modIndexMap = newModIndexMap;
797 else if (modIndexMap!=newModIndexMap) {
809 AlSymMat * symBigMatrix=
dynamic_cast<AlSymMat*
>(
m_bigmatrix);
810 AlSpaMat * spaBigMatrix=
dynamic_cast<AlSpaMat*
>(
m_bigmatrix);
817 int nDoF=modIndexMap.size();
821 AlSymMat newMatrix(nDoF);
823 if (
sc==StatusCode::SUCCESS)
824 *symBigMatrix += newMatrix;
828 throw std::logic_error(
"Unhandled matrix type");
831 AlSpaMat newMatrix(nDoF);
834 if (
sc==StatusCode::SUCCESS) {
836 *spaBigMatrix += newMatrix;
838 *spaBigMatrix = newMatrix;
842 if (
sc==StatusCode::FAILURE) {
848 ATH_MSG_WARNING(
"matrix not expected format! Changing m_wSqMatrix to "<<!triang);
◆ accumulateFromFiles()
bool MatrixTool::accumulateFromFiles |
( |
| ) |
|
|
virtual |
◆ accumulateFromTFiles()
bool MatrixTool::accumulateFromTFiles |
( |
| ) |
|
Store Files in a tfile.
Definition at line 957 of file MatrixTool.cxx.
960 int nDoF=alignParList->
size();
963 std::map<int,unsigned long long> modIndexMap;
964 std::map<int,unsigned long long> DoFMap;
965 double totalscale=0.;
967 AlSymMat * symBigMatrix=
dynamic_cast<AlSymMat*
>(
m_bigmatrix);
968 AlSpaMat * spaBigMatrix=
dynamic_cast<AlSpaMat*
>(
m_bigmatrix);
970 AlSpaMat *accumMatrix =
nullptr;
978 int numberOfReadErrors = 0;
980 struct rusage myusage{};
981 int itworked = getrusage(RUSAGE_SELF,&myusage);
985 long intialMemUse = myusage.ru_maxrss;
995 itworked = getrusage(RUSAGE_SELF,&myusage);
996 ATH_MSG_DEBUG(
"Memory usage [MB], total " << myusage.ru_maxrss/1024 <<
", increase " << (myusage.ru_maxrss-intialMemUse)/1024);
1001 ++numberOfReadErrors;
1006 std::map<int,unsigned long long> newModIndexMap;
1008 TVectorD* myModuleIDs;
1009 myModuleIDs = (TVectorD*)
myFile->Get(
"ModuleID");
1011 ++numberOfReadErrors;
1016 for (
int i(0);
i<myModuleIDs->GetNrows(); ++
i){
1018 double source = (*myModuleIDs)(
i);
1027 std::map<int,unsigned long long> newDoFMap;
1030 myDoFs = (TVectorD*)
myFile->Get(
"dof");
1032 ++numberOfReadErrors;
1037 for (
int i(0);
i<myDoFs->GetNrows(); ++
i){
1050 ++numberOfReadErrors;
1055 double scale=(*Scale)(0);
1056 totalscale +=
scale;
1063 ++numberOfReadErrors;
1068 AlVec* newVector =
new AlVec(nDoF);
1073 msg(
MSG::FATAL) <<
"vector wrong size! newVector size " << newVector->size()
1089 for (
int i=0;
i<nDoF;
i++) {
1090 (*newVector)[
i] = (*vector)(
i);
1097 }
else if (DoFMap!=newDoFMap) {
1103 modIndexMap = newModIndexMap;
1104 }
else if (modIndexMap!=newModIndexMap) {
1122 ++numberOfReadErrors;
1130 accumMatrix =
new AlSpaMat(nDoF);
1131 ATH_MSG_DEBUG(
"Matrix size b4 "<< accumMatrix->ptrMap()->size() );
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];
1143 i = (myRow.GetColPtr())[jj];
1145 (*accumMatrix)[
i][j] = myElement;
1148 ATH_MSG_DEBUG(
"Matrix size AF "<< accumMatrix->ptrMap()->size() );
1150 }
else if ( accumMatrix) {
1151 ATH_MSG_DEBUG(
"Matrix size b4 "<< accumMatrix->ptrMap()->size() );
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];
1162 i = (myRow.GetColPtr())[jj];
1164 (*accumMatrix)[
i][j] += myElement;
1167 ATH_MSG_DEBUG(
"Matrix size AF "<< accumMatrix->ptrMap()->size() );
1171 ++numberOfReadErrors;
1184 ++numberOfReadErrors;
1189 tracks = (TVectorD*)
myFile->Get(
"Tracks");
1192 ++numberOfReadErrors;
1197 if(
hits->GetNrows() != TotalHits.GetNrows() ){
1200 ++numberOfReadErrors;
1205 TotalHits += (*hits);
1206 TotalTracks += (*tracks);
1214 itworked = getrusage(RUSAGE_SELF,&myusage);
1215 ATH_MSG_DEBUG(
"Memory usage [MB], total " << myusage.ru_maxrss/1024 <<
", increase " << (myusage.ru_maxrss-intialMemUse)/1024);
1224 AlSymMat newMatrix(nDoF);
1226 for (
int i=0;
i<nDoF;
i++) {
1227 for (
int j=0;j<=
i;j++) {
1228 newMatrix[
i][j] = (*accumMatrix)[
i][j];
1232 *symBigMatrix += newMatrix;
1234 }
else if (spaBigMatrix) {
1235 ATH_MSG_DEBUG(
"should reassign matrix "<< spaBigMatrix->ptrMap()->size() );
1236 *spaBigMatrix += *accumMatrix;
1237 ATH_MSG_DEBUG(
"?????? "<< spaBigMatrix->ptrMap()->size() );
1244 AlignModuleList::const_iterator imod =
moduleList->begin();
1245 AlignModuleList::const_iterator imod_end =
moduleList->end();
1248 for(; imod != imod_end; ++imod, ++
index ) {
1252 totalhits += (
int)TotalHits(
index);
◆ addFirstDerivative()
void MatrixTool::addFirstDerivative |
( |
int |
irow, |
|
|
double |
firstderiv |
|
) |
| |
|
virtual |
◆ addFirstDerivatives() [1/2]
void MatrixTool::addFirstDerivatives |
( |
AlVec * |
vector | ) |
|
|
virtual |
◆ addFirstDerivatives() [2/2]
void MatrixTool::addFirstDerivatives |
( |
std::list< int, double > & |
derivatives | ) |
|
|
virtual |
◆ addModule()
void IMatrixTool::addModule |
( |
int |
alignModuleIndex, |
|
|
int |
nAlignParam |
|
) |
| |
|
inlineinherited |
◆ addSecondDerivative()
void MatrixTool::addSecondDerivative |
( |
int |
irow, |
|
|
int |
icol, |
|
|
double |
secondderiv |
|
) |
| |
|
virtual |
◆ addSecondDerivatives() [1/2]
void MatrixTool::addSecondDerivatives |
( |
AlSymMatBase * |
matrix | ) |
|
|
virtual |
◆ addSecondDerivatives() [2/2]
void MatrixTool::addSecondDerivatives |
( |
std::list< std::pair< int, int >, double > & |
derivatives | ) |
|
|
virtual |
◆ allocateMatrix()
StatusCode MatrixTool::allocateMatrix |
( |
int |
nDoF = 0 | ) |
|
|
virtual |
◆ declareGaudiProperty() [1/4]
specialization for handling Gaudi::Property<SG::VarHandleKeyArray>
Definition at line 170 of file AthCommonDataStore.h.
175 hndl.documentation());
◆ declareGaudiProperty() [2/4]
specialization for handling Gaudi::Property<SG::VarHandleKey>
Definition at line 156 of file AthCommonDataStore.h.
161 hndl.documentation());
◆ declareGaudiProperty() [3/4]
specialization for handling Gaudi::Property<SG::VarHandleBase>
Definition at line 184 of file AthCommonDataStore.h.
189 hndl.documentation());
◆ declareGaudiProperty() [4/4]
◆ declareProperty() [1/6]
Declare a new Gaudi property.
- Parameters
-
name | Name of the property. |
hndl | Object holding the property value. |
doc | Documentation string for the property. |
This is the version for types that derive from SG::VarHandleBase
. The property value object is put on the input and output lists as appropriate; then we forward to the base class.
Definition at line 245 of file AthCommonDataStore.h.
250 this->declare(hndl.
vhKey());
251 hndl.
vhKey().setOwner(
this);
253 return PBASE::declareProperty(
name,hndl,
doc);
◆ declareProperty() [2/6]
Declare a new Gaudi property.
- Parameters
-
name | Name of the property. |
hndl | Object holding the property value. |
doc | Documentation string for the property. |
This is the version for types that derive from SG::VarHandleKey
. The property value object is put on the input and output lists as appropriate; then we forward to the base class.
Definition at line 221 of file AthCommonDataStore.h.
229 return PBASE::declareProperty(
name,hndl,
doc);
◆ declareProperty() [3/6]
◆ declareProperty() [4/6]
Declare a new Gaudi property.
- Parameters
-
name | Name of the property. |
property | Object holding the property value. |
doc | Documentation string for the property. |
This is the generic version, for types that do not derive from SG::VarHandleKey
. It just forwards to the base class version of declareProperty
.
Definition at line 333 of file AthCommonDataStore.h.
338 return PBASE::declareProperty(
name, property,
doc);
◆ declareProperty() [5/6]
Declare a new Gaudi property.
- Parameters
-
name | Name of the property. |
property | Object holding the property value. |
doc | Documentation string for the property. |
This dispatches to either the generic declareProperty
or the one for VarHandle/Key/KeyArray.
Definition at line 352 of file AthCommonDataStore.h.
◆ declareProperty() [6/6]
◆ detStore()
◆ entryNumber()
int IMatrixTool::entryNumber |
( |
int |
alignModuleIndex | ) |
|
|
inlineinherited |
◆ evtStore() [1/2]
◆ evtStore() [2/2]
◆ extraDeps_update_handler()
Add StoreName to extra input/output deps as needed.
use the logic of the VarHandleKey to parse the DataObjID keys supplied via the ExtraInputs and ExtraOuputs Properties to add the StoreName if it's not explicitly given
◆ fillVecMods()
int MatrixTool::fillVecMods |
( |
| ) |
|
|
staticprivate |
◆ finalize()
StatusCode MatrixTool::finalize |
( |
| ) |
|
initialize
Definition at line 223 of file MatrixTool.cxx.
227 return StatusCode::SUCCESS;
◆ initialize()
StatusCode MatrixTool::initialize |
( |
| ) |
|
initialize
Definition at line 202 of file MatrixTool.cxx.
211 return StatusCode::FAILURE;
214 ATH_MSG_INFO(
"Retrieving data from the following files: ");
219 return StatusCode::SUCCESS;
◆ inputHandles()
Return this algorithm's input handles.
We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.
◆ interfaceID()
const InterfaceID & IMatrixTool::interfaceID |
( |
| ) |
|
|
inlinestaticinherited |
Retrieve interface ID.
Definition at line 110 of file IMatrixTool.h.
111 return IID_TRKALIGNINTERFACES_IMatrixTool;
◆ msg() [1/2]
◆ msg() [2/2]
◆ msgLvl()
◆ nHits()
int Trk::IMatrixTool::nHits |
( |
| ) |
const |
|
inlineinherited |
◆ nMeasurements()
int Trk::IMatrixTool::nMeasurements |
( |
| ) |
const |
|
inlineinherited |
◆ nTracks()
int Trk::IMatrixTool::nTracks |
( |
| ) |
const |
|
inlineinherited |
◆ outputHandles()
Return this algorithm's output handles.
We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.
◆ postSolvingLapack()
Definition at line 1958 of file MatrixTool.cxx.
1962 if(
z.ncol() !=
size) {
1963 msg(MSG::ERROR)<<
"Eigenvector matrix has incorrect size : "<<
z.ncol()<<
" != "<<
size<<
endmsg;
1978 ATH_MSG_INFO(
"writing the eigenvectors in a matrix: "<<
z.nrow() <<
"x" <<
z.ncol());
1989 if (
sc!=StatusCode::SUCCESS)
1990 msg(MSG::ERROR)<<
"Problem writing eigenvector matrix"<<
endmsg;
1996 ATH_MSG_INFO(
"writing the eigenvectors in a vector: "<<
w.size());
2000 sc =
w.WriteEigenvalueVec(
"eigenvalues.bin",
true);
2002 if (
sc!=StatusCode::SUCCESS)
2003 msg(MSG::ERROR)<<
"Problem writing eigenvector matrix"<<
endmsg;
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;
2017 const double eigenvalue_threshold = 1
e-19;
2022 ATH_MSG_INFO(
" Starting the automatic Weak Mode Removal method");
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());
2032 bool wm_stop =
false;
2039 (*Align_db)[
i] = (-D[
i]/
w[
i]);
2040 (*Align_error_db)[
i] = sqrt(1.0/
w[
i]/
m_scale);
2042 if (
w[
i]<eigenvalue_threshold) {
2044 <<
" removed as eigenvalue lower than the threshold " << eigenvalue_threshold
2046 (*AlignPull)[
i] = 0.0;
2050 (*AlignPull)[
i] = (*Align_db)[
i] / (*Align_error_db)[
i];
2052 ATH_MSG_DEBUG(
i <<
". AlignPar: " << (*Align_db)[
i] <<
" +- " << (*Align_error_db)[
i]
2053 <<
" (pull: " << (*AlignPull)[
i] <<
") ; w[i]: " <<
w[
i]);
2066 <<
" removed as pull is lower than " <<
m_pullcut <<
": "
2067 << (*AlignPull)[
i]);
2088 <<
" removed as diff between eigenvalues, " <<
w[
i] <<
" and " <<
w[
i+1]
2110 <<
" removed as diff between corrections, " <<
w[
i] <<
" and " <<
w[
i+1]
2111 <<
", is greater than "
2125 delete Align_error_db;
2128 Align_error_db =
nullptr;
2129 AlignPull =
nullptr;
2177 AlVec deltafull(
size);
2181 CLHEP::HepSymMatrix *
cov =
nullptr;
2185 cov =
new CLHEP::HepSymMatrix(
size,0);
2188 *
m_logStream<<
"/------ The Eigenvalue Spectrum -------"<<std::endl;
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;
2212 for(
int j=0;j<
size;j++) {
2213 errSq[j] +=
z[
i][j] *
z[
i][j] /
w[
i];
2215 for(
int k=0;
k<=j;
k++)
2223 *
m_logStream<<
"\\----- End of Eigenvalue Spectrum -----"<<std::endl;
2232 double param = delta[
i];
2233 double err = sqrt(2.*std::fabs(errSq[
i]));
2236 AlignPar * alignPar=(*alignParList)[idof];
2239 double sigma = alignPar->sigma();
2251 alignPar->setPar(param,
err);
2252 ATH_MSG_DEBUG(
"set param to "<<param<<
" for alignPar "<<alignPar);
2260 double norm1st = dChi2->norm();
2263 *
m_logStream<<
"norm of first derivative : "<<norm1st<<std::endl;
2267 double dist = ( (*d2Chi2) * deltafull + (*dChi2) ).
norm();
2270 *
m_logStream<<
"distance to solution : "<<dist<<std::endl;
2273 double chi2 = delta * (*d2Chi2) * delta * .5;
◆ prepareBinaryFiles()
void MatrixTool::prepareBinaryFiles |
( |
int |
solveOption | ) |
|
|
virtual |
◆ printGlobalSolution() [1/2]
void MatrixTool::printGlobalSolution |
( |
std::ostream & |
os, |
|
|
const CLHEP::HepSymMatrix * |
cov |
|
) |
| |
Definition at line 1621 of file MatrixTool.cxx.
1625 AlignModuleList::const_iterator imod = alignModules->begin();
1626 AlignModuleList::const_iterator imod_end = alignModules->end();
1627 for( ; imod!=imod_end; ++imod) {
1631 int thisNDoF = alignPars->
size();
1634 CLHEP::HepSymMatrix * covsub =
nullptr;
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();
1646 for (
int j=0;j<=
i;++j) {
1647 int jpar = alignPars->
at(j)->index();
1648 double sigma_j = alignPars->
at(j)->sigma();
1655 (*covsub)[
i][j] = (*cov)[iActive][jActive] * sigma_i * sigma_j;
1664 os <<
"--------------------------------------------------------------------------------" << std::endl;
◆ printGlobalSolution() [2/2]
void MatrixTool::printGlobalSolution |
( |
std::ostream & |
os, |
|
|
const TMatrixDSym * |
cov |
|
) |
| |
Definition at line 1668 of file MatrixTool.cxx.
1670 CLHEP::HepSymMatrix *
cov =
nullptr;
1672 int nsize = cov0->GetNrows();
1673 cov =
new CLHEP::HepSymMatrix(nsize,0);
1675 for(
int i=0;
i<nsize;
i++)
1676 for(
int j=0; j<=
i; j++)
1677 (*
cov)[
i][j] = (*cov0)[
i][j];
◆ printModuleSolution()
void MatrixTool::printModuleSolution |
( |
std::ostream & |
os, |
|
|
const AlignModule * |
module, |
|
|
const CLHEP::HepSymMatrix * |
cov |
|
) |
| const |
namespace { class RestoreIOSFlags { public: RestoreIOSFlags (std::ostream &os) : m_os(&os), m_flags(m_os->flags()), m_precision(m_os->precision()) {} ~RestoreIOSFlags() { m_os->flags(m_flags); m_os->precision(m_precision); } private: std::ostream *m_os; std::ios_base::fmtflags m_flags; std::streamsize m_precision; }; }
Definition at line 1706 of file MatrixTool.cxx.
1708 boost::io::ios_all_saver ias(
os );
1710 os <<
"--------------------------------------------------------------------------------" << std::endl;
1711 os <<
"Alignment parameters for module: " <<
module->name() << std::endl;
1712 os <<
"Number of tracks passing: " <<
module->nTracks() << std::endl;
1714 os <<
"Number of hits too small: "<<
module->nHits()<<
" < "<<
m_minNumHits<<
" Skipping the module"<<std::endl;
1718 os <<
"Number of tracks too small: "<<
module->nTracks()<<
" < "<<
m_minNumTrks<<
" Skipping the module"<<std::endl;
1721 os <<
"Number of hits seen: " <<
module->nHits() << std::endl;
1722 os <<
"Number of tracks seen: " <<
module->nTracks() << std::endl;
1725 int thisNDoF = alignPars->
size();
1727 if(alignPars->
empty())
1728 os <<
"No active parameters" << std::endl;
1733 os.unsetf(std::ios_base::floatfield);
1734 os << std::setiosflags(std::ios_base::left) << std::setprecision(5);
1739 for ( ; ipar != ipar_end; ++ipar) {
1741 os << std::setw(10) <<
par->dumpType()
1742 << std::setw(12) <<
par->par() <<
" +/- " << std::setw(12) <<
par->err()
1748 CLHEP::HepSymMatrix corrsub(thisNDoF,0);
1749 for(
int irow=0; irow<thisNDoF; ++irow)
1752 os <<
"Local correlation matrix: " << corrsub <<
std::flush;
◆ readHitmaps()
void MatrixTool::readHitmaps |
( |
| ) |
|
|
private |
Definition at line 2455 of file MatrixTool.cxx.
2464 for(
int imap=0;imap<
nFiles; imap++) {
2474 if(nextmap.nrow()!=
nModules || nextmap.ncol()!=2) {
2476 <<nextmap.nrow()<<
" x "<<nextmap.ncol()<<
"), should be ("<<
nModules<<
" x 2). Skipping.");
2483 AlignModuleList::const_iterator imod =
moduleList->begin();
2484 AlignModuleList::const_iterator imod_end =
moduleList->end();
2487 for(; imod != imod_end; ++imod) {
2493 totalhits += (
int)hitmap[
index][0];
2501 ATH_MSG_INFO(
"Hitmap accumulated from "<<
nFiles<<
" files with total of "<<totalhits<<
" hits.");
◆ renounce()
◆ renounceArray()
◆ setLogStream()
virtual void Trk::IMatrixTool::setLogStream |
( |
std::ostream * |
os | ) |
|
|
inlinevirtualinherited |
sets the output stream for the logfile
Definition at line 95 of file IMatrixTool.h.
◆ setNHits()
void Trk::IMatrixTool::setNHits |
( |
int |
n | ) |
|
|
inlineinherited |
set module identifier
set number of hits
Definition at line 83 of file IMatrixTool.h.
◆ setNMeasurements()
void Trk::IMatrixTool::setNMeasurements |
( |
int |
n | ) |
|
|
inlineinherited |
◆ setNTracks()
void Trk::IMatrixTool::setNTracks |
( |
int |
n | ) |
|
|
inlineinherited |
◆ solve()
int MatrixTool::solve |
( |
| ) |
|
|
virtual |
solves for alignment parameters
Implements Trk::IMatrixTool.
Definition at line 1268 of file MatrixTool.cxx.
1283 double dummyVersion(2.);
1286 std::map<int,unsigned long long> modIndexMap;
1287 std::map<int,std::string> modNameMap;
1290 modIndexMap[
i]=(*alignPars)[
i]->alignModule()->identify().get_compact();
1291 modNameMap [
i]=(*alignPars)[
i]->alignModule()->name();
1298 if (!sc1.isSuccess() || !sc2.isSuccess()) {
1299 msg(MSG::ERROR)<<
"problem writing matrix or vector"<<
endmsg;
1312 if (!sc1.isSuccess() || !sc2.isSuccess()) {
1313 msg(MSG::ERROR)<<
"problem writing matrix or vector"<<
endmsg;
1339 ATH_MSG_DEBUG(
"rescaling the matrix/vector and applying the soft-mode-cut");
1342 int nDoF = alignParList->
size();
1345 const AlSymMat * chkMatrix =
dynamic_cast<const AlSymMat*
>(
m_bigmatrix);
1348 for (
int i=0;
i<nDoF;
i++) {
1350 double sigma_i = (*alignParList)[
i]->sigma();
1351 double softCut = 2 *
pow( (*alignParList)[
i]->softCut() , -2 );
1352 (*m_bigvector)[
i] *= sigma_i;
1354 for (
int j=0;j<=
i;j++) {
1356 if ((*chkMatrix)[
i][j] != 0.) {
1357 double sigma_j = (*alignParList)[j]->sigma();
1358 (*m_bigmatrix)[
i][j] *= sigma_i * sigma_j;
1368 (*alignParList)[
i]->setFirstDeriv((*
m_bigvector)[
i]/sigma_i);
1369 (*alignParList)[
i]->setSecndDeriv((*
m_bigmatrix)[
i][
i]/sigma_i/sigma_i);
1376 int i =
p.first.first;
1377 int j =
p.first.second;
1380 double sigma_i = (*alignParList)[
i]->sigma();
1381 double sigma_j = (*alignParList)[j]->sigma();
1383 (*m_bigmatrix)[
i][j] *= sigma_i * sigma_j;
1388 for (
int i=0;
i<nDoF;
i++) {
1390 double sigma_i = (*alignParList)[
i]->sigma();
1391 (*m_bigvector)[
i] *= sigma_i;
1393 double softCut = 2 *
pow( (*alignParList)[
i]->softCut() , -2 );
1397 (*alignParList)[
i]->setFirstDeriv((*
m_bigvector)[
i]/sigma_i);
1398 (*alignParList)[
i]->setSecndDeriv((*
m_bigmatrix)[
i][
i]/sigma_i/sigma_i);
1402 unsigned long long OldPixelIdentifier = 37769216;
1403 unsigned long long IBLIdentifier = 33574912;
1405 unsigned long long SCT_ECA_8_Identifier = 218116096;
1406 std::string SCT_ECA_8_Name =
"SCT/EndcapA/Disk_8";
1411 ATH_MSG_INFO(
"Javi: Printing (*alignParList)[i]->alignModule()->identify32()");
1414 for(
int i=0;
i<nDoF;
i++)
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);
1434 {
ATH_MSG_INFO(
"SCT ECA Last Disk DoF have been skipped in the solving because DeactivateSCT_ECA_LastDisk is set to True");
1437 {
ATH_MSG_INFO(
"SCT ECA Last Disk DoF have been skipped in the solving because DeactivateSCT_ECA_LastDisk is set to True");
1443 {
ATH_MSG_INFO(
"Pixel DoF have been skipped in the solving because AlignIBLbutNotPixel is set to True");
1448 {
ATH_MSG_INFO(
"IBL DoF have been skipped in the solving because AlignPixelbutNotIBL is set to True");
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");
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");
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");
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");
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");
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");
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");
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");
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");
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");
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");
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");
1533 ATH_MSG_INFO(
"Spurious removal not implemented at the moment.");
◆ solveCLHEP()
int MatrixTool::solveCLHEP |
( |
| ) |
|
|
private |
Definition at line 391 of file MatrixTool.cxx.
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;
402 clock_t starttime = clock();
414 for (
int i=0;
i<nDoF;
i++)
415 for (
int j=0;j<nDoF;j++)
419 for (
int i=0;
i<nDoF;
i++)
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++) {
428 (*dChi2)[iActive] = (*m_bigvector)[
i];
429 for (
int jActive=0;jActive<
m_aNDoF;jActive++) {
431 (*d2Chi2)[iActive][jActive] = (*m_bigmatrix)[
i][j];
439 CLHEP::HepVector delta(
m_aNDoF,0);
440 CLHEP::HepVector deltafull(
m_aNDoF,0);
449 *
m_logStream<<
"Running matrix inversion"<<std::endl;
454 msg(MSG::ERROR)<<
"CLHEP inversion status flag = "<<ierr<<
endmsg;
458 *
m_logStream<<
"CLHEP inversion status flag = "<<ierr<<std::endl;
461 delta =
cov * (*dChi2);
471 CLHEP::HepSymMatrix cov2 = *d2Chi2 * .5;
476 msg(MSG::WARNING)<<
"Second CLHEP inversion status flag = "<<ierr2<<
endmsg;
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]) > 1
e-5 ) {
481 msg(MSG::WARNING)<<
"Something's wrong with the matrix inversion: delta["<<
i<<
"] = "<<delta[
i]<<
" delta2["<<
i<<
"] = "<<delta2[
i]<<
endmsg;
487 *
m_logStream<<
"CLHEP inversion status flag for halfed matrix = "<<ierr2<<std::endl;
488 *
m_logStream<<
"Matrix inversion check failed"<<std::endl;
498 *
m_logStream<<
"Running diagonalization"<<std::endl;
500 CLHEP::HepSymMatrix D = *d2Chi2;
501 CLHEP::HepMatrix U = CLHEP::diagonalize( &D );
511 if(D[j][j] < D[
i][
i]) {
526 CLHEP::HepVector eigenvector(
m_aNDoF);
529 *
m_logStream<<
"/------ The Eigenvalue Spectrum -------"<<std::endl;
532 for(
int imode=0; imode<
m_aNDoF; ++imode) {
535 for(
int irow=0; irow<
m_aNDoF; ++irow)
536 eigenvector[irow] = U[irow][imode];
540 double eigenvalue = D[imode][imode];
543 double evdotb =
dot(*dChi2,eigenvector);
544 CLHEP::HepVector thisdelta = evdotb/eigenvalue * eigenvector;
545 deltafull += thisdelta;
550 *
m_logStream<<
"| skipping eigenvalue "<<eigenvalue<<std::endl;
555 *
m_logStream<<
"| skipping eigenvalue "<<eigenvalue<<std::endl;
564 for(
int irow=0; irow<
m_aNDoF; ++irow)
566 cov[irow][
icol] += eigenvector[irow] * eigenvector[
icol] / eigenvalue;
576 *
m_logStream<<
"\\----- End of Eigenvalue Spectrum -----"<<std::endl;
583 clock_t stoptime = clock();
584 double totaltime = (stoptime-starttime)/
double(CLOCKS_PER_SEC);
585 ATH_MSG_INFO(
"Time spent in solveCLHEP: "<<totaltime<<
" s");
590 for (
int iAdof=0;iAdof<
m_aNDoF;iAdof++) {
593 AlignPar * alignPar=(*alignParList)[idof];
595 double sigma = alignPar->sigma();
596 double param = -delta[iAdof] *
sigma;
597 double err = std::sqrt(2.*std::fabs(
cov[iAdof][iAdof])) *
sigma;
602 alignPar->setPar(param,
err);
603 ATH_MSG_DEBUG(
"set param to "<<param<<
" for alignPar "<<alignPar);
611 *
m_logStream<<
"norm of first derivative : "<<dChi2->norm()<<std::endl;
614 double dist = ( - (*d2Chi2) * deltafull + (*dChi2) ).
norm();
615 *
m_logStream<<
"distance to solution : "<<dist<<std::endl;
618 double chi2 = d2Chi2->similarity(delta) * .5;
622 *
m_logStream<<
"time spent in solve : "<<totaltime<<
" s"<<std::endl;
◆ solveLapack()
int MatrixTool::solveLapack |
( |
| ) |
|
|
private |
Definition at line 1765 of file MatrixTool.cxx.
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;
1776 AlSymMat* aBetterMat =
new AlSymMat(
m_aNDoF);
1777 AlVec* aBetterVec =
new AlVec(
m_aNDoF);
1778 for (
int iActive=0;iActive<
m_aNDoF;iActive++) {
1780 (*aBetterVec)[iActive] = (*m_bigvector)[
i];
1781 for (
int jActive=0;jActive<
m_aNDoF;jActive++) {
1783 (*aBetterMat)[iActive][jActive] = (*m_bigmatrix)[
i][j];
1790 ATH_MSG_WARNING(
"Scaling requested but scale not set. Not scaling matrix and vector.");
1797 ATH_MSG_DEBUG(
"Now Solving alignment using lapack diagonalization routine dspev...");
1800 const double tol = 1.e-20;
1802 double determ = (*aBetterMat).determinant();
1804 if (fabs(determ) < tol)
1809 AlSymMat * d2Chi2 =
nullptr;
1811 d2Chi2 =
new AlSymMat(*aBetterMat);
1813 clock_t starttime = clock();
1821 int info = (*aBetterMat).diagonalize(jobz,
w,
z);
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");
1830 double time_solve = 0.;
1832 starttime = clock();
1835 time_solve = (stoptime-starttime)/
double(CLOCKS_PER_SEC);
1836 ATH_MSG_INFO(
" - time spent solving the system: "<<time_solve<<
" s");
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;
1843 ATH_MSG_ERROR(
"Problem in diagonalization. Solving skipped.");
1845 *
m_logStream<<
"time spent for diagonalization: "<<time_diag<<
" s"<<std::endl;
1849 *
m_logStream<<
"total time spent in solve: "<<time_diag+time_solve<<
" s"<<std::endl;
◆ solveLocal()
int MatrixTool::solveLocal |
( |
| ) |
|
|
private |
Definition at line 633 of file MatrixTool.cxx.
637 *
m_logStream<<
"*************************************************************"<<std::endl;
638 *
m_logStream<<
"************** solving using Local method *****************"<<std::endl;
639 *
m_logStream<<
"*************************************************************"<<std::endl;
643 double totalChi2(0.);
647 AlignModuleList::const_iterator imod = alignModules->begin();
648 AlignModuleList::const_iterator imod_end = alignModules->end();
649 for( ; imod!=imod_end; ++imod) {
657 int thisNDoF = alignPars->
size();
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];
690 totalNDoF += thisNDoF;
692 CLHEP::HepSymMatrix
cov(d2Chi2);
703 CLHEP::HepVector delta =
cov * dChi2;
711 for (
int idof=0;idof<thisNDoF;++idof) {
714 double sigma = alignPar->sigma();
715 double param = -delta[idof] *
sigma;
716 double err = std::sqrt(2.*std::fabs(
cov[idof][idof])) *
sigma;
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);
731 *
m_logStream<<
"CLHEP inversion status flag = "<<ierr<<std::endl;
734 double chi2 = d2Chi2.similarity(delta) * .5;
736 *
m_logStream<<
"delta(chi2) of the alignment change : "<<
chi2<<
" / "<<thisNDoF<<std::endl;
741 *
m_logStream<<
"--------------------------------------------------------------------------------"<<std::endl;
742 *
m_logStream<<
"Total delta(chi2) of the alignment change from the local method : "<<totalChi2<<
" / "<<totalNDoF<<std::endl;
◆ solveROOT()
int MatrixTool::solveROOT |
( |
| ) |
|
|
private |
Definition at line 265 of file MatrixTool.cxx.
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;
276 clock_t starttime = clock();
286 for (
int i=0;
i<nDoF;
i++)
287 for (
int j=0;j<nDoF;j++)
291 for (
int i=0;
i<nDoF;
i++)
297 double * firstderiv =
new double[
m_aNDoF];
298 for (
int iActive=0;iActive<
m_aNDoF;iActive++) {
300 firstderiv[iActive] = (*m_bigvector)[
i];
301 for (
int jActive=0;jActive<
m_aNDoF;jActive++) {
303 secderiv[iActive*
m_aNDoF+jActive] = (*m_bigmatrix)[
i][j];
323 TMatrixDSym ainv(
c.Invert(
status));
325 TVectorD
r(
b.GetNrows());
330 clock_t stoptime = clock();
331 double totaltime = (stoptime-starttime)/
double(CLOCKS_PER_SEC);
332 ATH_MSG_INFO(
"Time spent in solveROOT: "<<totaltime<<
" s");
335 msg(MSG::ERROR)<<
"ROOT inversion failed"<<
endmsg;
345 for (
int iAdof=0;iAdof<
m_aNDoF;iAdof++) {
348 AlignPar * alignPar=(*alignParList)[idof];
350 double sigma = alignPar->sigma();
351 double param = -
r[iAdof] *
sigma;
352 double err = std::sqrt(2.*std::fabs(ainv(iAdof,iAdof))) *
sigma;
357 alignPar->setPar(param,
err);
358 ATH_MSG_DEBUG(
"set param to "<<param<<
" for alignPar "<<alignPar);
369 *
m_logStream<<
"norm of first derivative : "<<sqrt(
b.Norm2Sqr())<<std::endl;
372 double dist = sqrt( (
b - (
a *
r) ).Norm2Sqr() );
373 *
m_logStream<<
"distance to solution : "<<dist<<std::endl;
376 double chi2 =
a.Similarity(
r) * .5;
380 *
m_logStream<<
"time spent in solve : "<<totaltime<<
" s"<<std::endl;
385 delete [] firstderiv;
◆ solveSparseEigen()
int MatrixTool::solveSparseEigen |
( |
| ) |
|
|
private |
Definition at line 2284 of file MatrixTool.cxx.
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;
2295 clock_t starttime = clock();
2299 AlSpaMat * ABetterMat =
nullptr;
2300 bool isCopy =
false;
2302 ATH_MSG_INFO(
"Converting Matrix Format for fast solving");
2303 ABetterMat =
new AlSpaMat(*(
dynamic_cast<AlSymMat*
>(
m_bigmatrix)));
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));
2311 ATH_MSG_ERROR(
"Cannot cast to neither AlSymMat nor AlSpaMat");
2319 const AlSpaMat * chkMatrix = ABetterMat;
2321 AlSpaMat * aBetterMat =
new AlSpaMat(
m_aNDoF);
2322 AlVec * aBetterVec =
new AlVec(
m_aNDoF);
2325 for (
int iActive=0;iActive<
m_aNDoF;iActive++) {
2327 (*aBetterVec)[iActive] = (*m_bigvector)[
i];
2328 for (
int jActive=0;jActive<
m_aNDoF;jActive++) {
2331 if ( (*chkMatrix)[iActive][jActive] != 0. )
2332 (*aBetterMat)[iActive][jActive]=(*ABetterMat)[
i][j];
2337 AlVec origVec(*aBetterVec);
2342 int info = (*aBetterMat).SolveWithEigen(*aBetterVec);
2347 *
m_logStream<<
"SolveWithEigen solving OK."<<std::endl;
2352 *
m_logStream<<
"SolveWithEigen error code (0 if OK) = "<<
info<<std::endl;
2357 ABetterMat =
nullptr;
2360 clock_t stoptime = clock();
2361 double totaltime = (stoptime-starttime)/
double(CLOCKS_PER_SEC);
2362 ATH_MSG_INFO(
"Time spent in SolveWithEigen: "<<totaltime<<
" s");
2369 double param = -(*aBetterVec)[
i];
2373 AlignPar * alignPar=(*alignParList)[idof];
2376 double sigma = alignPar->sigma();
2382 alignPar->setPar(param,
err);
2383 ATH_MSG_DEBUG(
"set param to "<<param<<
" for alignPar "<<alignPar);
2388 CLHEP::HepSymMatrix *
cov =
nullptr;
2392 *
m_logStream<<
"norm of first derivative : "<<origVec.norm()<<std::endl;
2395 double dist = ( (*aBetterMat) * (*aBetterVec) - origVec ).
norm();
2396 *
m_logStream<<
"distance to solution : "<<dist<<std::endl;
2399 double chi2 = (*aBetterVec) * (*aBetterMat) * (*aBetterVec) * .5;
2403 *
m_logStream<<
"time spent in solve : "<<totaltime<<
" s"<<std::endl;
◆ spuriousRemoval()
StatusCode MatrixTool::spuriousRemoval |
( |
| ) |
|
|
private |
Definition at line 1866 of file MatrixTool.cxx.
1874 const double tol = 1.e-20;
1877 if (std::fabs(determ) < tol)
1883 if (fillvecmods==0) {
1895 return StatusCode::SUCCESS;
1897 else if (fillvecmods==2)
1898 return StatusCode::FAILURE;
1954 return StatusCode::SUCCESS;
◆ storeInTFile()
void MatrixTool::storeInTFile |
( |
const TString & |
filename | ) |
|
Store Files in a tfile.
Definition at line 862 of file MatrixTool.cxx.
870 int nDoF = alignParList->
size();
877 double *
val =
new double[nDoF];
878 for (
int i=0;
i<nDoF;
i++) {
879 val[
i] = (*m_bigvector)[
i];
882 TVectorD myTVector(nDoF,
val);
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();
896 for(; imod != imod_end; ++imod) {
903 TVectorD hitmapHits(
nModules, hitmapA);
904 TVectorD hitmapTracks(
nModules, hitmapB);
912 hitmapHits.Write(
"Hits");
913 hitmapTracks.Write(
"Tracks");
914 myTMatrix->Write(
"Matrix");
915 myTVector.Write(
"Vector");
918 scale.Write(
"Scale");
921 double *moduleInfoA =
new double[nDoF];
922 double *dofInfoA =
new double[nDoF];
924 if (
sizeof(
unsigned long long) !=
sizeof(
double))
925 ATH_MSG_ERROR(
"Module Identifiers will not be saved. sizeof(double)!=sizeof(ulonglong)");
932 uint64_t id = (*alignPars)[
i]->alignModule()->identify().get_compact();
936 uint64_t dof = (*alignPars)[
i]->paramType();
941 TVectorD moduleIDs(nDoF, moduleInfoA) ;
942 TVectorD moduleDoFs(nDoF,dofInfoA);
943 delete [] moduleInfoA;
944 moduleIDs.Write(
"ModuleID");
945 moduleDoFs.Write(
"dof");
◆ sysInitialize()
◆ sysStart()
Handle START transition.
We override this in order to make sure that conditions handle keys can cache a pointer to the conditions container.
◆ updateVHKA()
◆ writeHitmap()
void MatrixTool::writeHitmap |
( |
| ) |
|
|
private |
Definition at line 2418 of file MatrixTool.cxx.
2426 AlignModuleList::const_iterator imod =
moduleList->begin();
2427 AlignModuleList::const_iterator imod_end =
moduleList->end();
2429 for(; imod != imod_end; ++imod) {
2442 if (
sc!=StatusCode::SUCCESS)
2446 sc = hitmap.Write(
"hitmap.txt",
false, 0);
2447 if (
sc!=StatusCode::SUCCESS)
2448 ATH_MSG_ERROR(
"Problem writing hitmap matrix to text file");
◆ m_activeIndices
std::vector<int> Trk::MatrixTool::m_activeIndices |
|
private |
vector of indices which pass the min-hits cut
Definition at line 209 of file MatrixTool.h.
◆ m_Align_db_step
float Trk::MatrixTool::m_Align_db_step |
|
private |
corr in the diagonal basis step for the third pass in the auto weak mode removal method
Definition at line 162 of file MatrixTool.h.
◆ m_AlignIBLbutNotPixel
bool Trk::MatrixTool::m_AlignIBLbutNotPixel |
|
private |
◆ m_alignModuleMap
std::map<int,int> Trk::IMatrixTool::m_alignModuleMap |
|
privateinherited |
◆ m_alignModuleTool
◆ m_AlignPixelbutNotIBL
bool Trk::MatrixTool::m_AlignPixelbutNotIBL |
|
private |
◆ m_aNDoF
int Trk::MatrixTool::m_aNDoF |
|
private |
number of active DoF (size of m_activeIndices)
Definition at line 210 of file MatrixTool.h.
◆ m_bigmatrix
matrix to contain second derivative terms to be used for alignment
Definition at line 145 of file MatrixTool.h.
◆ m_bigvector
AlVec* Trk::MatrixTool::m_bigvector |
|
private |
vector to contain first derivative terms to be used for alignment
Definition at line 148 of file MatrixTool.h.
◆ m_calculateFullCovariance
bool Trk::MatrixTool::m_calculateFullCovariance |
|
private |
◆ m_calDet
bool Trk::MatrixTool::m_calDet |
|
private |
Compute bigmatrix's determinant ?
Definition at line 164 of file MatrixTool.h.
◆ m_DeactivateSCT_ECA_LastDisk
bool Trk::MatrixTool::m_DeactivateSCT_ECA_LastDisk |
|
private |
◆ m_detStore
◆ m_diagonalize
bool Trk::MatrixTool::m_diagonalize |
|
private |
run diagonalization instead of inversion
Definition at line 153 of file MatrixTool.h.
◆ m_eigenvalueStep
float Trk::MatrixTool::m_eigenvalueStep |
|
private |
eigenvalue step for the second pass in the automatic weak mode removal method
Definition at line 161 of file MatrixTool.h.
◆ m_eigenvaluethreshold
double Trk::MatrixTool::m_eigenvaluethreshold |
|
private |
◆ m_evtStore
◆ m_inputHitmapFiles
std::vector<std::string> Trk::MatrixTool::m_inputHitmapFiles |
|
private |
input binary files containing the hitmaps
Definition at line 205 of file MatrixTool.h.
◆ m_inputMatrixFiles
std::vector<std::string> Trk::MatrixTool::m_inputMatrixFiles |
|
private |
input binary files containing matrix terms
Definition at line 202 of file MatrixTool.h.
◆ m_inputTFiles
std::vector<std::string> Trk::MatrixTool::m_inputTFiles |
|
private |
input binary files containing matrix terms
Definition at line 207 of file MatrixTool.h.
◆ m_inputVectorFiles
std::vector<std::string> Trk::MatrixTool::m_inputVectorFiles |
|
private |
input binary files containing vector terms
Definition at line 203 of file MatrixTool.h.
◆ m_logStream
std::ostream* Trk::IMatrixTool::m_logStream |
|
protectedinherited |
◆ m_maxReadErrors
int Trk::MatrixTool::m_maxReadErrors |
|
private |
maximum number of reading TFile errors
Definition at line 212 of file MatrixTool.h.
◆ m_minNumHits
int Trk::MatrixTool::m_minNumHits |
|
private |
cut on the minimum number of hits per module
Definition at line 158 of file MatrixTool.h.
◆ m_minNumTrks
int Trk::MatrixTool::m_minNumTrks |
|
private |
cut on the minimum number of tracks per module
Definition at line 159 of file MatrixTool.h.
◆ m_modcut
int Trk::MatrixTool::m_modcut |
|
private |
cut on the weak modes which number is <par_modcut
Definition at line 157 of file MatrixTool.h.
◆ m_nentries
int Trk::IMatrixTool::m_nentries |
|
privateinherited |
◆ m_nHits
int Trk::IMatrixTool::m_nHits |
|
protectedinherited |
◆ m_nMeasurements
int Trk::IMatrixTool::m_nMeasurements |
|
protectedinherited |
◆ m_nTracks
int Trk::IMatrixTool::m_nTracks |
|
protectedinherited |
◆ m_pathbin
std::string Trk::MatrixTool::m_pathbin |
|
private |
◆ m_pathtxt
std::string Trk::MatrixTool::m_pathtxt |
|
private |
◆ m_prefixName
std::string Trk::MatrixTool::m_prefixName |
|
private |
◆ m_pullcut
float Trk::MatrixTool::m_pullcut |
|
private |
pull cut for the automatic weak mode removal method
Definition at line 160 of file MatrixTool.h.
◆ m_readHitmaps
bool Trk::MatrixTool::m_readHitmaps |
|
private |
◆ m_readTFiles
bool Trk::MatrixTool::m_readTFiles |
|
private |
write out files to a root file
Definition at line 177 of file MatrixTool.h.
◆ m_Remove_IBL_Rx
bool Trk::MatrixTool::m_Remove_IBL_Rx |
|
private |
◆ m_Remove_IBL_Ry
bool Trk::MatrixTool::m_Remove_IBL_Ry |
|
private |
◆ m_Remove_IBL_Rz
bool Trk::MatrixTool::m_Remove_IBL_Rz |
|
private |
◆ m_Remove_IBL_Tx
bool Trk::MatrixTool::m_Remove_IBL_Tx |
|
private |
◆ m_Remove_IBL_Ty
bool Trk::MatrixTool::m_Remove_IBL_Ty |
|
private |
◆ m_Remove_IBL_Tz
bool Trk::MatrixTool::m_Remove_IBL_Tz |
|
private |
◆ m_Remove_Pixel_Rx
bool Trk::MatrixTool::m_Remove_Pixel_Rx |
|
private |
◆ m_Remove_Pixel_Ry
bool Trk::MatrixTool::m_Remove_Pixel_Ry |
|
private |
◆ m_Remove_Pixel_Rz
bool Trk::MatrixTool::m_Remove_Pixel_Rz |
|
private |
◆ m_Remove_Pixel_Tx
bool Trk::MatrixTool::m_Remove_Pixel_Tx |
|
private |
◆ m_Remove_Pixel_Ty
bool Trk::MatrixTool::m_Remove_Pixel_Ty |
|
private |
◆ m_Remove_Pixel_Tz
bool Trk::MatrixTool::m_Remove_Pixel_Tz |
|
private |
◆ m_removeSpurious
bool Trk::MatrixTool::m_removeSpurious |
|
private |
◆ m_runLocal
bool Trk::MatrixTool::m_runLocal |
|
private |
Run solving using Local method.
Definition at line 179 of file MatrixTool.h.
◆ m_scalaMatName
std::string Trk::MatrixTool::m_scalaMatName |
|
private |
◆ m_scalaVecName
std::string Trk::MatrixTool::m_scalaVecName |
|
private |
◆ m_scale
double Trk::MatrixTool::m_scale |
|
private |
scale for big matrix and vector normalization
Definition at line 181 of file MatrixTool.h.
◆ m_scaleMatrix
bool Trk::MatrixTool::m_scaleMatrix |
|
private |
scale matrix by number of hits before solving
Definition at line 182 of file MatrixTool.h.
◆ m_softEigenmodeCut
double Trk::MatrixTool::m_softEigenmodeCut |
|
private |
add constant to diagonal to effectively cut on weak eigenmodes
Definition at line 184 of file MatrixTool.h.
◆ m_solveOption
int Trk::MatrixTool::m_solveOption |
|
private |
◆ m_tfileName
std::string Trk::MatrixTool::m_tfileName |
|
private |
◆ m_useSparse
bool Trk::MatrixTool::m_useSparse |
|
private |
◆ m_varHandleArraysDeclared
◆ m_vhka
◆ m_writeEigenMat
bool Trk::MatrixTool::m_writeEigenMat |
|
private |
write eigenvalues and eigenvectors into files ?
Definition at line 168 of file MatrixTool.h.
◆ m_writeEigenMatTxt
bool Trk::MatrixTool::m_writeEigenMatTxt |
|
private |
also write eigenvalues and eigenvectors into txt files ?
Definition at line 169 of file MatrixTool.h.
◆ m_writeHitmap
bool Trk::MatrixTool::m_writeHitmap |
|
private |
◆ m_writeHitmapTxt
bool Trk::MatrixTool::m_writeHitmapTxt |
|
private |
◆ m_writeMat
bool Trk::MatrixTool::m_writeMat |
|
private |
write big matrix and vector into files ?
Definition at line 166 of file MatrixTool.h.
◆ m_writeMatTxt
bool Trk::MatrixTool::m_writeMatTxt |
|
private |
also write big matrix and vector into txt files ?
Definition at line 167 of file MatrixTool.h.
◆ m_writeModuleNames
bool Trk::MatrixTool::m_writeModuleNames |
|
private |
◆ m_writeTFile
bool Trk::MatrixTool::m_writeTFile |
|
private |
write out files to a root file
Definition at line 176 of file MatrixTool.h.
◆ m_wSqMatrix
bool Trk::MatrixTool::m_wSqMatrix |
|
private |
write a triangular matrix by default (true: square format) ?
Definition at line 165 of file MatrixTool.h.
The documentation for this class was generated from the following files:
JetConstituentVector::iterator iterator
virtual double determinant()=0
@ z
global position (cartesian)
std::string find(const std::string &s)
return a remapped string
std::vector< AlignModule * > AlignModuleList
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
virtual void SetPathTxt(const std::string &)=0
StoreGateSvc_t m_evtStore
Pointer to StoreGate (event store by default)
std::vector< SG::VarHandleKeyArray * > m_vhka
bool msgLvl(const MSG::Level lvl) const
StatusCode WritePartial(const std::string &, bool, double, std::map< int, unsigned long long >, float)
virtual void setOwner(IDataHandleHolder *o)=0
const datamap * ptrMap() const
class TMatrixTSparse< double > TMatrixDSparse
AlignModule is a grouping of TrkDetElementBase objects, grouped according to the type of alignment,...
virtual TMatrixDSparse * makeTMatrix()=0
::StatusCode StatusCode
StatusCode definition for legacy code.
double chi2(TH1 *h0, TH1 *h1)
virtual StatusCode Write(const std::string &, bool, bool, double, float)=0
StoreGateSvc_t m_detStore
Pointer to StoreGate (detector store by default)
void Scale(TH1 *h, double d=1)
def dot(G, fn, nodesToHighlight=[])
virtual void renounce()=0
std::conditional< std::is_base_of< SG::VarHandleKeyArray, T >::value, VarHandleKeyArrayType, type2 >::type type
AlignPar contains all the information related to an alignment parameter of a particular align module ...
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
virtual void SetPathBin(const std::string &)=0
#define ATH_MSG_WARNING(x)
void SetPathBin(const std::string &)
SG::VarHandleKey & vhKey()
Return a non-const reference to the HandleKey.
void SetPathTxt(const std::string &)
const T * at(size_type n) const
Access an element, as an rvalue.
constexpr int pow(int base, int exp) noexcept
vec_fb< typename boost::int_t< sizeof(T) *8 >::exact, N > ivec
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
size_type size() const noexcept
Returns the number of elements in the collection.
Gaudi::Details::PropertyBase & declareGaudiProperty(Gaudi::Property< T > &hndl, const SG::VarHandleKeyType &)
specialization for handling Gaudi::Property<SG::VarHandleKey>
bool empty() const noexcept
Returns true if the collection is empty.
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.