00001 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentExtendedCorrelationsStore.h"
00002
00003 #include "Alignment/CommonAlignment/interface/Alignable.h"
00004 #include "Alignment/CommonAlignment/interface/AlignmentParameters.h"
00005 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentExtendedCorrelationsEntry.h"
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00008 #include "FWCore/Utilities/interface/Exception.h"
00009
00010
00011 AlignmentExtendedCorrelationsStore::AlignmentExtendedCorrelationsStore( const edm::ParameterSet& config )
00012 {
00013 theMaxUpdates = config.getParameter<int>( "MaxUpdates" );
00014 theCut = config.getParameter<double>( "CutValue" );
00015 theWeight = config.getParameter<double>( "Weight" );
00016
00017 edm::LogInfo("Alignment") << "@SUB=AlignmentExtendedCorrelationsStore::AlignmentExtendedCorrelationsStore"
00018 << "Created.";
00019 }
00020
00021
00022 void AlignmentExtendedCorrelationsStore::correlations( Alignable* ap1, Alignable* ap2,
00023 AlgebraicSymMatrix& cov, int row, int col ) const
00024 {
00025 static Alignable* previousAlignable = 0;
00026 static ExtendedCorrelationsTable* previousCorrelations;
00027
00028
00029 if ( ap1 == 0 ) { previousAlignable = 0; return; }
00030
00031 bool transpose = ( ap2 > ap1 );
00032 if ( transpose ) std::swap( ap1, ap2 );
00033
00034 if ( ap1 == previousAlignable )
00035 {
00036 ExtendedCorrelationsTable::const_iterator itC2 = previousCorrelations->find( ap2 );
00037 if ( itC2 != previousCorrelations->end() )
00038 {
00039 transpose ?
00040 fillCovarianceT( ap1, ap2, (*itC2).second, cov, row, col ) :
00041 fillCovariance( ap1, ap2, (*itC2).second, cov, row, col );
00042 }
00043 }
00044 else
00045 {
00046 ExtendedCorrelations::const_iterator itC1 = theCorrelations.find( ap1 );
00047 if ( itC1 != theCorrelations.end() )
00048 {
00049 previousAlignable = ap1;
00050 previousCorrelations = (*itC1).second;
00051
00052 ExtendedCorrelationsTable::const_iterator itC2 = (*itC1).second->find( ap2 );
00053 if ( itC2 != (*itC1).second->end() )
00054 {
00055 transpose ?
00056 fillCovarianceT( ap1, ap2, (*itC2).second, cov, row, col ) :
00057 fillCovariance( ap1, ap2, (*itC2).second, cov, row, col );
00058 }
00059 }
00060 }
00061
00062
00063 return;
00064 }
00065
00066
00067 void AlignmentExtendedCorrelationsStore::setCorrelations( Alignable* ap1, Alignable* ap2,
00068 const AlgebraicSymMatrix& cov, int row, int col )
00069 {
00070 static Alignable* previousAlignable = 0;
00071 static ExtendedCorrelationsTable* previousCorrelations;
00072
00073
00074 if ( ap1 == 0 ) { previousAlignable = 0; return; }
00075
00076 bool transpose = ( ap2 > ap1 );
00077 if ( transpose ) std::swap( ap1, ap2 );
00078
00079 if ( ap1 == previousAlignable )
00080 {
00081 fillCorrelationsTable( ap1, ap2, previousCorrelations, cov, row, col, transpose );
00082 }
00083 else
00084 {
00085 ExtendedCorrelations::iterator itC = theCorrelations.find( ap1 );
00086 if ( itC != theCorrelations.end() )
00087 {
00088 fillCorrelationsTable( ap1, ap2, itC->second, cov, row, col, transpose );
00089 previousAlignable = ap1;
00090 previousCorrelations = itC->second;
00091 }
00092 else
00093 {
00094
00095 ExtendedCorrelationsTable* newTable = new ExtendedCorrelationsTable;
00096 fillCorrelationsTable( ap1, ap2, newTable, cov, row, col, transpose );
00097
00098 theCorrelations[ap1] = newTable;
00099
00100 previousAlignable = ap1;
00101 previousCorrelations = newTable;
00102 }
00103 }
00104 }
00105
00106
00107 void AlignmentExtendedCorrelationsStore::setCorrelations( Alignable* ap1, Alignable* ap2, AlgebraicMatrix& mat )
00108 {
00109 bool transpose = ( ap2 > ap1 );
00110 if ( transpose ) std::swap( ap1, ap2 );
00111
00112 ExtendedCorrelations::iterator itC1 = theCorrelations.find( ap1 );
00113 if ( itC1 != theCorrelations.end() )
00114 {
00115 ExtendedCorrelationsTable::iterator itC2 = itC1->second->find( ap1 );
00116 if ( itC2 != itC1->second->end() )
00117 {
00118 itC2->second = transpose ? ExtendedCorrelationsEntry( mat.T() ) : ExtendedCorrelationsEntry( mat );
00119 }
00120 else
00121 {
00122 (*itC1->second)[ap2] = transpose ? ExtendedCorrelationsEntry( mat.T() ) : ExtendedCorrelationsEntry( mat );
00123 }
00124 }
00125 else
00126 {
00127 ExtendedCorrelationsTable* newTable = new ExtendedCorrelationsTable;
00128 (*newTable)[ap2] = transpose ? ExtendedCorrelationsEntry( mat.T() ) : ExtendedCorrelationsEntry( mat );
00129 theCorrelations[ap1] = newTable;
00130 }
00131 }
00132
00133
00134 void AlignmentExtendedCorrelationsStore::getCorrelations( Alignable* ap1, Alignable* ap2, AlgebraicMatrix& mat ) const
00135 {
00136 bool transpose = ( ap2 > ap1 );
00137 if ( transpose ) std::swap( ap1, ap2 );
00138
00139 ExtendedCorrelations::const_iterator itC1 = theCorrelations.find( ap1 );
00140 if ( itC1 != theCorrelations.end() )
00141 {
00142 ExtendedCorrelationsTable::const_iterator itC2 = itC1->second->find( ap2 );
00143 if ( itC2 != itC1->second->end() )
00144 {
00145 mat = transpose ? itC2->second.matrix().T() : itC2->second.matrix();
00146 return;
00147 }
00148 }
00149
00150 mat = AlgebraicMatrix();
00151 }
00152
00153
00154 bool AlignmentExtendedCorrelationsStore::correlationsAvailable( Alignable* ap1, Alignable* ap2 ) const
00155 {
00156 bool transpose = ( ap2 > ap1 );
00157 if ( transpose ) std::swap( ap1, ap2 );
00158
00159 ExtendedCorrelations::const_iterator itC1 = theCorrelations.find( ap1 );
00160 if ( itC1 != theCorrelations.end() )
00161 {
00162 ExtendedCorrelationsTable::const_iterator itC2 = itC1->second->find( ap2 );
00163 if ( itC2 != itC1->second->end() ) return true;
00164 }
00165 return false;
00166 }
00167
00168
00169 void AlignmentExtendedCorrelationsStore::resetCorrelations( void )
00170 {
00171 ExtendedCorrelations::iterator itC;
00172 for ( itC = theCorrelations.begin(); itC != theCorrelations.end(); ++itC ) delete (*itC).second;
00173 theCorrelations.erase( theCorrelations.begin(), theCorrelations.end() );
00174
00175
00176 AlgebraicSymMatrix dummy;
00177 correlations( 0, 0, dummy, 0, 0 );
00178 setCorrelations( 0, 0, dummy, 0, 0 );
00179 }
00180
00181
00182 unsigned int AlignmentExtendedCorrelationsStore::size( void ) const
00183 {
00184 unsigned int size = 0;
00185 ExtendedCorrelations::const_iterator itC;
00186 for ( itC = theCorrelations.begin(); itC != theCorrelations.end(); ++itC )
00187 size += itC->second->size();
00188
00189 return size;
00190 }
00191
00192
00193 void
00194 AlignmentExtendedCorrelationsStore::fillCorrelationsTable( Alignable* ap1, Alignable* ap2,
00195 ExtendedCorrelationsTable* table,
00196 const AlgebraicSymMatrix& cov,
00197 int row, int col, bool transpose )
00198 {
00199 ExtendedCorrelationsTable::iterator itC = table->find( ap2 );
00200
00201 if ( itC != table->end() )
00202 {
00203
00204
00205 transpose ?
00206 readFromCovarianceT( ap1, ap2, itC->second, cov, row, col ) :
00207 readFromCovariance( ap1, ap2, itC->second, cov, row, col );
00208
00209
00210 }
00211 else
00212 {
00213 int nRow = ap1->alignmentParameters()->numSelected();
00214 int nCol = ap2->alignmentParameters()->numSelected();
00215 ExtendedCorrelationsEntry newEntry( nRow, nCol );
00216
00217 transpose ?
00218 readFromCovarianceT( ap1, ap2, newEntry, cov, row, col ) :
00219 readFromCovariance( ap1, ap2, newEntry, cov, row, col );
00220
00221 (*table)[ap2] = newEntry;
00222 }
00223 }
00224
00225
00226 void
00227 AlignmentExtendedCorrelationsStore::fillCovariance( Alignable* ap1, Alignable* ap2, const ExtendedCorrelationsEntry& entry,
00228 AlgebraicSymMatrix& cov, int row, int col ) const
00229 {
00230 int nRow = entry.numRow();
00231 int nCol = entry.numCol();
00232
00233 for ( int iRow = 0; iRow < nRow; ++iRow )
00234 {
00235 double factor = sqrt(cov[row+iRow][row+iRow]);
00236 if ( std::isnan(factor) ) throw cms::Exception("LogicError") << "[AlignmentExtendedCorrelationsStore::fillCovariance] "
00237 << "NaN-factor: sqrt(" << cov[row+iRow][row+iRow] << ")";
00238
00239 for ( int jCol = 0; jCol < nCol; ++jCol )
00240 cov[row+iRow][col+jCol] = entry( iRow, jCol )*factor;
00241 }
00242
00243 for ( int jCol = 0; jCol < nCol; ++jCol )
00244 {
00245 double factor = sqrt(cov[col+jCol][col+jCol]);
00246 if ( std::isnan(factor) ) throw cms::Exception("LogicError") << "[AlignmentExtendedCorrelationsStore::fillCovariance] "
00247 << "NaN-factor: sqrt(" << cov[col+jCol][col+jCol] << ")";
00248
00249 for ( int iRow = 0; iRow < nRow; ++iRow )
00250 cov[row+iRow][col+jCol] *= factor;
00251 }
00252 }
00253
00254
00255 void
00256 AlignmentExtendedCorrelationsStore::fillCovarianceT( Alignable* ap1, Alignable* ap2, const ExtendedCorrelationsEntry& entry,
00257 AlgebraicSymMatrix& cov, int row, int col ) const
00258 {
00259 int nRow = entry.numRow();
00260 int nCol = entry.numCol();
00261
00262 for ( int iRow = 0; iRow < nRow; ++iRow )
00263 {
00264 double factor = sqrt(cov[col+iRow][col+iRow]);
00265 if ( std::isnan(factor) ) throw cms::Exception("LogicError") << "[AlignmentExtendedCorrelationsStore::fillCovarianceT] "
00266 << "NaN-factor: sqrt(" << cov[col+iRow][col+iRow] << ")";
00267 for ( int jCol = 0; jCol < nCol; ++jCol )
00268 cov[row+jCol][col+iRow] = entry( iRow, jCol )*factor;
00269 }
00270
00271 for ( int jCol = 0; jCol < nCol; ++jCol )
00272 {
00273 double factor = sqrt(cov[row+jCol][row+jCol]);
00274 if ( std::isnan(factor) ) throw cms::Exception("LogicError") << "[AlignmentExtendedCorrelationsStore::fillCovarianceT] "
00275 << "NaN-factor: sqrt(" << cov[row+jCol][row+jCol] << ")";
00276 for ( int iRow = 0; iRow < nRow; ++iRow )
00277 cov[row+jCol][col+iRow] *= factor;
00278 }
00279
00280 }
00281
00282
00283 void
00284 AlignmentExtendedCorrelationsStore::readFromCovariance( Alignable* ap1, Alignable* ap2, ExtendedCorrelationsEntry& entry,
00285 const AlgebraicSymMatrix& cov, int row, int col )
00286 {
00287 int nRow = entry.numRow();
00288 int nCol = entry.numCol();
00289
00290 for ( int iRow = 0; iRow < nRow; ++iRow )
00291 {
00292 double factor = sqrt(cov[row+iRow][row+iRow]);
00293 for ( int jCol = 0; jCol < nCol; ++jCol )
00294 entry( iRow, jCol ) = cov[row+iRow][col+jCol]/factor;
00295 }
00296
00297 double maxCorr = 0;
00298
00299 for ( int jCol = 0; jCol < nCol; ++jCol )
00300 {
00301 double factor = sqrt(cov[col+jCol][col+jCol]);
00302 for ( int iRow = 0; iRow < nRow; ++iRow )
00303 {
00304 entry( iRow, jCol ) /= factor;
00305 if ( fabs( entry( iRow, jCol ) ) > maxCorr ) maxCorr = fabs( entry( iRow, jCol ) );
00306 }
00307 }
00308
00309 resizeCorruptCorrelations( entry, maxCorr );
00310 }
00311
00312
00313 void
00314 AlignmentExtendedCorrelationsStore::readFromCovarianceT( Alignable* ap1, Alignable* ap2, ExtendedCorrelationsEntry& entry,
00315 const AlgebraicSymMatrix& cov, int row, int col )
00316 {
00317 int nRow = entry.numRow();
00318 int nCol = entry.numCol();
00319
00320 for ( int iRow = 0; iRow < nRow; ++iRow )
00321 {
00322 double factor = sqrt(cov[col+iRow][col+iRow]);
00323 for ( int jCol = 0; jCol < nCol; ++jCol )
00324 entry( iRow, jCol ) = cov[row+jCol][col+iRow]/factor;
00325 }
00326
00327 double maxCorr = 0;
00328
00329 for ( int jCol = 0; jCol < nCol; ++jCol )
00330 {
00331 double factor = sqrt(cov[row+jCol][row+jCol]);
00332 for ( int iRow = 0; iRow < nRow; ++iRow )
00333 {
00334 entry( iRow, jCol ) /= factor;
00335 if ( fabs( entry( iRow, jCol ) ) > maxCorr ) maxCorr = fabs( entry( iRow, jCol ) );
00336 }
00337 }
00338
00339 resizeCorruptCorrelations( entry, maxCorr );
00340 }
00341
00342
00343 void
00344 AlignmentExtendedCorrelationsStore::resizeCorruptCorrelations( ExtendedCorrelationsEntry& entry,
00345 double maxCorr )
00346 {
00347 if ( maxCorr > 1. )
00348 {
00349 entry *= theWeight/maxCorr;
00350 }
00351 else if ( maxCorr > theCut )
00352 {
00353 entry *= 1. - ( maxCorr - theCut )/( 1. - theCut )*( 1. - theWeight );
00354 }
00355 }