CommonAlignmentParametrization/src/CompositeAlignmentParameters.cc

0001 ///  $Date: 2007/10/08 15:56:00 $
0002 ///  $Revision: 1.12 $
0003 /// (last update by $Author: cklae $)
0004
0005 #include "Alignment/CommonAlignmentParametrization/interface/CompositeAlignmentParameters.h"
0006
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 #include "Alignment/CommonAlignment/interface/Alignable.h"
0010 #include "Alignment/CommonAlignment/interface/AlignableDetOrUnitPtr.h"
0011 #include "Alignment/CommonAlignmentParametrization/interface/CompositeAlignmentDerivativesExtractor.h"
0012 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0013
0014 #include "Alignment/CommonAlignment/interface/AlignmentParameters.h"
0015
0016 //__________________________________________________________________________________________________
0017 CompositeAlignmentParameters::CompositeAlignmentParameters(const AlgebraicVector &par,
0018                                                            const AlgebraicSymMatrix &cov,
0019                                                            const Components &comp)
0020     : theData(DataContainer(new AlignmentParametersData(par, cov))), theComponents(comp) {}
0021
0022 //__________________________________________________________________________________________________
0023 CompositeAlignmentParameters::CompositeAlignmentParameters(const AlgebraicVector &par,
0024                                                            const AlgebraicSymMatrix &cov,
0025                                                            const Components &comp,
0026                                                            const AlignableDetToAlignableMap &alimap,
0027                                                            const Aliposmap &aliposmap,
0028                                                            const Alilenmap &alilenmap)
0029     : theData(DataContainer(new AlignmentParametersData(par, cov))),
0030       theComponents(comp),
0031       theAlignableDetToAlignableMap(alimap),
0032       theAliposmap(aliposmap),
0033       theAlilenmap(alilenmap) {}
0034
0035 //__________________________________________________________________________________________________
0036 CompositeAlignmentParameters::CompositeAlignmentParameters(const DataContainer &data,
0037                                                            const Components &comp,
0038                                                            const AlignableDetToAlignableMap &alimap,
0039                                                            const Aliposmap &aliposmap,
0040                                                            const Alilenmap &alilenmap)
0041     : theData(data),
0042       theComponents(comp),
0043       theAlignableDetToAlignableMap(alimap),
0044       theAliposmap(aliposmap),
0045       theAlilenmap(alilenmap) {}
0046
0047 //__________________________________________________________________________________________________
0048 CompositeAlignmentParameters::~CompositeAlignmentParameters() {}
0049
0050 //__________________________________________________________________________________________________
0051 CompositeAlignmentParameters *CompositeAlignmentParameters::clone(const AlgebraicVector &par,
0052                                                                   const AlgebraicSymMatrix &cov) const {
0053   CompositeAlignmentParameters *cap = new CompositeAlignmentParameters(par, cov, components());
0054
0055   return cap;
0056 }
0057
0058 //__________________________________________________________________________________________________
0059 CompositeAlignmentParameters *CompositeAlignmentParameters::clone(const AlgebraicVector &par,
0060                                                                   const AlgebraicSymMatrix &cov,
0061                                                                   const AlignableDetToAlignableMap &alimap,
0062                                                                   const Aliposmap &aliposmap,
0063                                                                   const Alilenmap &alilenmap) const {
0064   CompositeAlignmentParameters *cap =
0065       new CompositeAlignmentParameters(par, cov, components(), alimap, aliposmap, alilenmap);
0066
0067   return cap;
0068 }
0069
0070 //__________________________________________________________________________________________________
0071 CompositeAlignmentParameters::Components CompositeAlignmentParameters::components() const { return theComponents; }
0072
0073 //__________________________________________________________________________________________________
0074 // full derivatives for a composed object
0075 AlgebraicMatrix CompositeAlignmentParameters::derivatives(const std::vector<TrajectoryStateOnSurface> &tsosvec,
0076                                                           const std::vector<AlignableDet *> &alidetvec) const {
0077   std::vector<AlignableDetOrUnitPtr> detOrUnits;
0078   this->convert(alidetvec, detOrUnits);
0079
0080   return this->derivatives(tsosvec, detOrUnits);
0081 }
0082
0083 AlgebraicMatrix CompositeAlignmentParameters::derivatives(const std::vector<TrajectoryStateOnSurface> &tsosvec,
0084                                                           const std::vector<AlignableDetOrUnitPtr> &alidetvec) const {
0085   align::Alignables alivec;
0086   for (std::vector<AlignableDetOrUnitPtr>::const_iterator it = alidetvec.begin(); it != alidetvec.end(); ++it)
0087     alivec.push_back(alignableFromAlignableDet(*it));
0088
0089   CompositeAlignmentDerivativesExtractor extractor(alivec, alidetvec, tsosvec);
0090   return extractor.derivatives();
0091 }
0092
0093 //__________________________________________________________________________________________________
0094 AlgebraicVector CompositeAlignmentParameters::correctionTerm(const std::vector<TrajectoryStateOnSurface> &tsosvec,
0095                                                              const std::vector<AlignableDet *> &alidetvec) const {
0096   std::vector<AlignableDetOrUnitPtr> detOrUnits;
0097   this->convert(alidetvec, detOrUnits);
0098
0099   return this->correctionTerm(tsosvec, detOrUnits);
0100 }
0101
0102 //__________________________________________________________________________________________________
0103 AlgebraicVector CompositeAlignmentParameters::correctionTerm(
0104     const std::vector<TrajectoryStateOnSurface> &tsosvec, const std::vector<AlignableDetOrUnitPtr> &alidetvec) const {
0105   align::Alignables alivec;
0106   for (std::vector<AlignableDetOrUnitPtr>::const_iterator it = alidetvec.begin(); it != alidetvec.end(); ++it)
0107     alivec.push_back(alignableFromAlignableDet(*it));
0108
0109   CompositeAlignmentDerivativesExtractor extractor(alivec, alidetvec, tsosvec);
0110   return extractor.correctionTerm();
0111 }
0112
0113 //__________________________________________________________________________________________________
0114 // assume all are selected
0115 AlgebraicMatrix CompositeAlignmentParameters::selectedDerivatives(const std::vector<TrajectoryStateOnSurface> &tsosvec,
0116                                                                   const std::vector<AlignableDet *> &alidetvec) const {
0117   return derivatives(tsosvec, alidetvec);
0118 }
0119 //__________________________________________________________________________________________________
0120 // assume all are selected
0121 AlgebraicMatrix CompositeAlignmentParameters::selectedDerivatives(
0122     const std::vector<TrajectoryStateOnSurface> &tsosvec, const std::vector<AlignableDetOrUnitPtr> &alidetvec) const {
0123   return derivatives(tsosvec, alidetvec);
0124 }
0125
0126 //__________________________________________________________________________________________________
0127 // only one (tsos,AlignableDet) as argument [for compatibility with base class]
0128 AlgebraicMatrix CompositeAlignmentParameters::derivatives(const TrajectoryStateOnSurface &tsos,
0129                                                           const AlignableDetOrUnitPtr &alidet) const {
0130   std::vector<TrajectoryStateOnSurface> tsosvec;
0131   std::vector<AlignableDetOrUnitPtr> alidetvec;
0132   tsosvec.push_back(tsos);
0133   alidetvec.push_back(alidet);
0134   return derivatives(tsosvec, alidetvec);
0135 }
0136
0137 //__________________________________________________________________________________________________
0138 // assume all are selected
0139 AlgebraicMatrix CompositeAlignmentParameters::selectedDerivatives(const TrajectoryStateOnSurface &tsos,
0140                                                                   const AlignableDetOrUnitPtr &alidet) const {
0141   return derivatives(tsos, alidet);
0142 }
0143
0144 // Derivatives ----------------------------------------------------------------
0145 // legacy methods
0146 // full derivatives for a composed object
0147 AlgebraicMatrix CompositeAlignmentParameters::derivativesLegacy(const std::vector<TrajectoryStateOnSurface> &tsosvec,
0148                                                                 const std::vector<AlignableDet *> &alidetvec) const {
0149   // sanity check: length of parameter argument vectors must be equal
0150   if (alidetvec.size() != tsosvec.size()) {
0151     edm::LogError("BadArgument") << " Inconsistent length of argument vectors! ";
0152     AlgebraicMatrix selderiv(1, 0);
0153     return selderiv;
0154   }
0155
0156   std::vector<AlgebraicMatrix> vecderiv;
0157   int nparam = 0;
0158
0159   std::vector<TrajectoryStateOnSurface>::const_iterator itsos = tsosvec.begin();
0160   for (std::vector<AlignableDet *>::const_iterator it = alidetvec.begin(); it != alidetvec.end(); ++it, ++itsos) {
0161     AlignableDet *ad = (*it);
0162     Alignable *ali = alignableFromAlignableDet(ad);
0163     AlignmentParameters *ap = ali->alignmentParameters();
0164     AlgebraicMatrix thisselderiv = ap->selectedDerivatives(*itsos, ad);
0165     vecderiv.push_back(thisselderiv);
0166     nparam += thisselderiv.num_row();
0167   }
0168
0169   int ipos = 1;
0170   AlgebraicMatrix selderiv(nparam, 2);
0171   for (std::vector<AlgebraicMatrix>::const_iterator imat = vecderiv.begin(); imat != vecderiv.end(); ++imat) {
0172     const AlgebraicMatrix &thisselderiv = (*imat);
0173     int npar = thisselderiv.num_row();
0174     selderiv.sub(ipos, 1, thisselderiv);
0175     ipos += npar;
0176   }
0177
0178   return selderiv;
0179 }
0180
0181 //__________________________________________________________________________________________________
0182 // assume all are selected
0183 AlgebraicMatrix CompositeAlignmentParameters::selectedDerivativesLegacy(
0184     const std::vector<TrajectoryStateOnSurface> &tsosvec, const std::vector<AlignableDet *> &alidetvec) const {
0185   return derivativesLegacy(tsosvec, alidetvec);
0186 }
0187
0188 //__________________________________________________________________________________________________
0189 // only one (tsos,AlignableDet) as argument [for compatibility with base class]
0190 AlgebraicMatrix CompositeAlignmentParameters::derivativesLegacy(const TrajectoryStateOnSurface &tsos,
0191                                                                 AlignableDet *alidet) const {
0192   std::vector<TrajectoryStateOnSurface> tsosvec;
0193   std::vector<AlignableDet *> alidetvec;
0194   tsosvec.push_back(tsos);
0195   alidetvec.push_back(alidet);
0196   return derivativesLegacy(tsosvec, alidetvec);
0197 }
0198
0199 //__________________________________________________________________________________________________
0200 // assume all are selected
0201 AlgebraicMatrix CompositeAlignmentParameters::selectedDerivativesLegacy(const TrajectoryStateOnSurface &tsos,
0202                                                                         AlignableDet *alidet) const {
0203   return derivativesLegacy(tsos, alidet);
0204 }
0205
0206 //__________________________________________________________________________________________________
0207 // finds Alignable corresponding to AlignableDet
0208 Alignable *CompositeAlignmentParameters::alignableFromAlignableDet(const AlignableDetOrUnitPtr &adet) const {
0209   AlignableDetToAlignableMap::const_iterator iali = theAlignableDetToAlignableMap.find(adet);
0210   if (iali != theAlignableDetToAlignableMap.end())
0211     return (*iali).second;
0212   else
0213     return nullptr;
0214 }
0215
0216 //__________________________________________________________________________________________________
0217 AlgebraicVector CompositeAlignmentParameters::parameterSubset(const align::Alignables &vec) const {
0218   const auto &sel = extractAlignables(vec);
0219
0220   const unsigned int nali = sel.size();
0221   int ndim = 0;
0222
0223   std::vector<int> posvec;
0224   std::vector<int> lenvec;
0225
0226   posvec.reserve(nali);
0227   lenvec.reserve(nali);
0228
0229   // iterate over input vector of alignables to determine size of result vector
0230   if (!extractPositionAndLength(sel, posvec, lenvec, ndim))
0231     return AlgebraicVector();
0232
0233   // OK, let's do the real work now
0234   AlgebraicVector result(ndim);
0235
0236   int resi = 0;
0237   for (unsigned int iali = 0; iali < nali; ++iali) {
0238     int posi = posvec[iali];
0239     int leni = lenvec[iali];
0240
0241     for (int ir = 0; ir < leni; ++ir)
0242       result[resi + ir] = theData->parameters()[posi - 1 + ir];
0243
0244     resi += leni;
0245   }
0246
0247   return result;
0248 }
0249
0250 //__________________________________________________________________________________________________
0251 // extract covariance matrix for a subset of alignables
0252 AlgebraicSymMatrix CompositeAlignmentParameters::covarianceSubset(const align::Alignables &vec) const {
0253   const auto &sel = extractAlignables(vec);
0254
0255   const unsigned int nali = sel.size();
0256   int ndim = 0;
0257
0258   std::vector<int> posvec;
0259   std::vector<int> lenvec;
0260
0261   posvec.reserve(nali);
0262   lenvec.reserve(nali);
0263
0264   // iterate over input vectors of alignables
0265   // to determine dimensions of result matrix
0266   if (!extractPositionAndLength(sel, posvec, lenvec, ndim))
0267     return AlgebraicSymMatrix();
0268
0269   // OK, let's do the real work now
0270   AlgebraicSymMatrix result(ndim);
0271
0272   int resi = 0;
0273   for (unsigned int iali = 0; iali < nali; ++iali) {
0274     int posi = posvec[iali];
0275     int leni = lenvec[iali];
0276
0277     int resj = 0;
0278     for (unsigned int jali = 0; jali <= iali; ++jali) {
0279       int posj = posvec[jali];
0280       int lenj = lenvec[jali];
0281
0282       for (int ir = 0; ir < leni; ++ir)
0283         for (int ic = 0; ic < lenj; ++ic)
0284           result[resi + ir][resj + ic] = theData->covariance()[posi - 1 + ir][posj - 1 + ic];
0285
0286       resj += lenj;
0287     }
0288     resi += leni;
0289   }
0290
0291   return result;
0292 }
0293
0294 //__________________________________________________________________________________________________
0295 // extract covariance matrix elements between two subsets of alignables
0296 AlgebraicMatrix CompositeAlignmentParameters::covarianceSubset(const align::Alignables &veci,
0297                                                                const align::Alignables &vecj) const {
0298   const auto &seli = extractAlignables(veci);
0299   const auto &selj = extractAlignables(vecj);
0300
0301   int ndimi = 0;
0302   int ndimj = 0;
0303
0304   std::vector<int> posveci;
0305   std::vector<int> lenveci;
0306   std::vector<int> posvecj;
0307   std::vector<int> lenvecj;
0308
0309   posveci.reserve(seli.size());
0310   lenveci.reserve(seli.size());
0311   posvecj.reserve(selj.size());
0312   lenvecj.reserve(selj.size());
0313
0314   // iterate over input vectors of alignables
0315   // to determine dimensions of result matrix
0316   if (!extractPositionAndLength(seli, posveci, lenveci, ndimi))
0317     return AlgebraicSymMatrix();
0318   // vector vecj
0319   if (!extractPositionAndLength(selj, posvecj, lenvecj, ndimj))
0320     return AlgebraicSymMatrix();
0321
0322   // OK, let's do the real work now
0323   AlgebraicMatrix result(ndimi, ndimj);
0324
0325   int resi = 0;
0326   for (unsigned int iali = 0; iali < seli.size(); ++iali) {
0327     int posi = posveci[iali];
0328     int leni = lenveci[iali];
0329
0330     int resj = 0;
0331     for (unsigned int jali = 0; jali < selj.size(); ++jali) {
0332       int posj = posvecj[jali];
0333       int lenj = lenvecj[jali];
0334
0335       for (int ir = 0; ir < leni; ++ir)
0336         for (int ic = 0; ic < lenj; ++ic)
0337           result[resi + ir][resj + ic] = theData->covariance()[posi - 1 + ir][posj - 1 + ic];
0338
0339       resj += lenj;
0340     }
0341     resi += leni;
0342   }
0343
0344   return result;
0345 }
0346
0347 //__________________________________________________________________________________________________
0348 // Extract position and length of parameters for a subset of Alignables.
0349 bool CompositeAlignmentParameters::extractPositionAndLength(const align::Alignables &alignables,
0350                                                             std::vector<int> &posvec,
0351                                                             std::vector<int> &lenvec,
0352                                                             int &length) const {
0353   length = 0;
0354
0355   for (const auto &it : alignables) {
0356     // check if in components
0357     if (std::find(theComponents.begin(), theComponents.end(), it) == theComponents.end()) {
0358       edm::LogError("NotFound") << "@SUB=CompositeAlignmentParameters::extractPositionAndLength"
0359                                 << "Alignable not found in components!";
0360       return false;
0361     }
0362
0363     // get pos/length
0364     Aliposmap::const_iterator iposmap = theAliposmap.find(it);
0365     Alilenmap::const_iterator ilenmap = theAlilenmap.find(it);
0366     if (iposmap == theAliposmap.end() || ilenmap == theAlilenmap.end()) {
0367       edm::LogError("NotFound") << "@SUB=CompositeAlignmentParameters::extractPositionAndLength"
0368                                 << "position/length not found for Alignable in maps!";
0369       return false;
0370     }
0371     posvec.push_back(iposmap->second);
0372     lenvec.push_back(ilenmap->second);
0373     length += ilenmap->second;
0374   }
0375
0376   return true;
0377 }
0378
0379 //__________________________________________________________________________________________________
0380 align::Alignables CompositeAlignmentParameters::extractAlignables(const align::Alignables &alignables) const {
0381   align::Alignables result;
0382
0383   for (const auto &itA : alignables) {
0384     if (std::find(result.begin(), result.end(), itA) == result.end())
0385       result.push_back(itA);
0386   }
0387
0388   return result;
0389 }
0390
0391 //__________________________________________________________________________________________________
0392 void CompositeAlignmentParameters::convert(const std::vector<AlignableDet *> &input,
0393                                            std::vector<AlignableDetOrUnitPtr> &output) const {
0394   output.clear();
0395   output.reserve(input.size());
0396
0397   std::vector<AlignableDet *>::const_iterator it, itEnd;
0398   for (it = input.begin(), itEnd = input.end(); it != itEnd; ++it)
0399     output.push_back(AlignableDetOrUnitPtr(*it));
0400 }