CommonTopologies/interface/GeometryAligner.h

0001 #ifndef Geometry_TrackingGeometryAligner_GeometryAligner_h
0002 #define Geometry_TrackingGeometryAligner_GeometryAligner_h
0003
0004 #include <vector>
0005 #include <algorithm>
0006 #include <iterator>
0007
0008 #include "FWCore/Utilities/interface/Exception.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include "Geometry/CommonTopologies/interface/DetPositioner.h"
0012
0013 #include "CondFormats/Alignment/interface/Alignments.h"
0014 #include "CondFormats/Alignment/interface/AlignmentErrors.h"
0015 #include "CondFormats/Alignment/interface/AlignmentErrorsExtended.h"
0016 #include "CondFormats/Alignment/interface/AlignTransform.h"
0017 #include "CondFormats/Alignment/interface/AlignTransformError.h"
0018 #include "CondFormats/Alignment/interface/AlignTransformErrorExtended.h"
0019 #include "CondFormats/Alignment/interface/AlignmentSurfaceDeformations.h"
0020
0021 #include "Geometry/CommonTopologies/interface/GeomDet.h"
0022 #include "DataFormats/GeometryCommonDetAlgo/interface/AlignmentPositionError.h"
0023 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
0024 #include "DataFormats/GeometrySurface/interface/Surface.h"
0025 #include "Geometry/CommonTopologies/interface/SurfaceDeformationFactory.h"
0026 #include "Geometry/CommonTopologies/interface/SurfaceDeformation.h"
0027
0028 class Alignments;
0029 class AlignmentSurfaceDeformations;
0030
0031 /// Class to update a given geometry with a set of alignments
0032
0033 class GeometryAligner : public DetPositioner {
0034 public:
0035   template <class C>
0036   void applyAlignments(const C* geometry,
0037                        const Alignments* alignments,
0038                        const AlignmentErrorsExtended* alignmentErrors,
0039                        const AlignTransform& globalCoordinates);
0040
0041   template <class C>
0042   void attachSurfaceDeformations(const C* geometry, const AlignmentSurfaceDeformations* surfaceDeformations);
0043
0044   inline void removeGlobalTransform(const Alignments* alignments,
0045                                     const AlignmentErrorsExtended* alignmentErrors,
0046                                     const AlignTransform& globalCoordinates,
0047                                     Alignments* newAlignments,
0048                                     AlignmentErrorsExtended* newAlignmentErrorsExtended);
0049 };
0050
0051 template <class C>
0052 void GeometryAligner::applyAlignments(const C* geometry,
0053                                       const Alignments* alignments,
0054                                       const AlignmentErrorsExtended* alignmentErrors,
0055                                       const AlignTransform& globalCoordinates) {
0056   edm::LogInfo("Alignment") << "@SUB=GeometryAligner::applyAlignments"
0057                             << "Starting to apply alignments.";
0058
0059   // Preliminary checks (or we can't loop!)
0060   if (alignments->m_align.size() != geometry->theMap.size())
0061     throw cms::Exception("GeometryMismatch") << "Size mismatch between geometry (size=" << geometry->theMap.size()
0062                                              << ") and alignments (size=" << alignments->m_align.size() << ")";
0063   if (alignments->m_align.size() != alignmentErrors->m_alignError.size())
0064     throw cms::Exception("GeometryMismatch")
0065         << "Size mismatch between geometry (size=" << geometry->theMap.size()
0066         << ") and alignment errors (size=" << alignmentErrors->m_alignError.size() << ")";
0067
0068   const AlignTransform::Translation& globalShift = globalCoordinates.translation();
0069   const AlignTransform::Rotation globalRotation = globalCoordinates.rotation();  // by value!
0070   const AlignTransform::Rotation inverseGlobalRotation = globalRotation.inverse();
0071
0072   // Parallel loop on alignments, alignment errors and geomdets
0073   std::vector<AlignTransform>::const_iterator iAlign = alignments->m_align.begin();
0074   std::vector<AlignTransformErrorExtended>::const_iterator iAlignError = alignmentErrors->m_alignError.begin();
0075   //copy  geometry->theMap to a real map to order it....
0076   std::map<unsigned int, GeomDet const*> theMap;
0077   std::copy(geometry->theMap.begin(), geometry->theMap.end(), std::inserter(theMap, theMap.begin()));
0078   unsigned int nAPE = 0;
0079   for (auto iPair = theMap.begin(); iPair != theMap.end(); ++iPair, ++iAlign, ++iAlignError) {
0080     // Check DetIds
0081     if ((*iPair).first != (*iAlign).rawId())
0082       throw cms::Exception("GeometryMismatch") << "DetId mismatch between geometry (rawId=" << (*iPair).first
0083                                                << ") and alignments (rawId=" << (*iAlign).rawId();
0084
0085     if ((*iPair).first != (*iAlignError).rawId())
0086       throw cms::Exception("GeometryMismatch") << "DetId mismatch between geometry (rawId=" << (*iPair).first
0087                                                << ") and alignment errors (rawId=" << (*iAlignError).rawId();
0088
0089     // Apply global correction
0090     CLHEP::Hep3Vector positionHep = globalRotation * CLHEP::Hep3Vector((*iAlign).translation()) + globalShift;
0091     CLHEP::HepRotation rotationHep = CLHEP::HepRotation((*iAlign).rotation()) * inverseGlobalRotation;
0092
0093     // Define new position/rotation objects and apply
0094     Surface::PositionType position(positionHep.x(), positionHep.y(), positionHep.z());
0095     Surface::RotationType rotation(rotationHep.xx(),
0096                                    rotationHep.xy(),
0097                                    rotationHep.xz(),
0098                                    rotationHep.yx(),
0099                                    rotationHep.yy(),
0100                                    rotationHep.yz(),
0101                                    rotationHep.zx(),
0102                                    rotationHep.zy(),
0103                                    rotationHep.zz());
0104     GeomDet* iGeomDet = const_cast<GeomDet*>((*iPair).second);
0105     this->setGeomDetPosition(*iGeomDet, position, rotation);
0106
0107     // Alignment Position Error only if non-zero to save memory
0108     GlobalErrorExtended error(asSMatrix<6>((*iAlignError).matrix()));
0109
0110     AlignmentPositionError ape(error);
0111     if (this->setAlignmentPositionError(*iGeomDet, ape))
0112       ++nAPE;
0113   }
0114
0115   edm::LogInfo("Alignment") << "@SUB=GeometryAligner::applyAlignments"
0116                             << "Finished to apply " << theMap.size() << " alignments with " << nAPE << " non-zero APE.";
0117 }
0118
0119 template <class C>
0120 void GeometryAligner::attachSurfaceDeformations(const C* geometry,
0121                                                 const AlignmentSurfaceDeformations* surfaceDeformations) {
0122   edm::LogInfo("Alignment") << "@SUB=GeometryAligner::attachSurfaceDeformations"
0123                             << "Starting to attach surface deformations.";
0124
0125   //copy geometry->theMapUnit to a real map to order it....
0126   std::map<unsigned int, GeomDetUnit const*> theMap;
0127   std::copy(geometry->theMapUnit.begin(), geometry->theMapUnit.end(), std::inserter(theMap, theMap.begin()));
0128
0129   unsigned int nSurfDef = 0;
0130   unsigned int itemIndex = 0;
0131   auto iPair = theMap.begin();
0132   for (std::vector<AlignmentSurfaceDeformations::Item>::const_iterator iItem = surfaceDeformations->items().begin();
0133        iItem != surfaceDeformations->items().end();
0134        ++iItem, ++iPair) {
0135     // Check DetIds
0136     // go forward in map of GeomDetUnits until DetId is found
0137     while ((*iPair).first != (*iItem).m_rawId) {
0138       // remove SurfaceDeformation from GeomDetUnit (i.e. set NULL pointer)
0139       GeomDetUnit* geomDetUnit = const_cast<GeomDetUnit*>((*iPair).second);
0140       this->setSurfaceDeformation(*geomDetUnit, nullptr);
0141
0142       ++iPair;
0143       if (iPair == theMap.end())
0144         throw cms::Exception("GeometryMismatch")
0145             << "GeomDetUnit with rawId=" << (*iItem).m_rawId << " not found in geometry";
0146     }
0147
0148     // get the parameters and put them into a vector
0149     AlignmentSurfaceDeformations::ParametersConstIteratorPair iteratorPair = surfaceDeformations->parameters(itemIndex);
0150     std::vector<double> parameters;
0151     std::copy(iteratorPair.first, iteratorPair.second, std::back_inserter(parameters));
0152
0153     // create SurfaceDeformation via factory
0154     SurfaceDeformation* surfDef = SurfaceDeformationFactory::create((*iItem).m_parametrizationType, parameters);
0155     GeomDetUnit* geomDetUnit = const_cast<GeomDetUnit*>((*iPair).second);
0156     this->setSurfaceDeformation(*geomDetUnit, surfDef);
0157     // delete is not needed since SurfaceDeformation is passed as a
0158     // DeepCopyPointerByClone which takes over ownership. Needs to be
0159     // cleaned up and checked once SurfaceDeformation are moved to
0160     // proxy topology classes
0161     //delete surfDef;
0162
0163     ++nSurfDef;
0164
0165     ++itemIndex;
0166   }
0167
0168   edm::LogInfo("Alignment") << "@SUB=GeometryAligner::attachSurfaceDeformations"
0169                             << "Finished to attach " << nSurfDef << " surface deformations.";
0170 }
0171
0172 void GeometryAligner::removeGlobalTransform(const Alignments* alignments,
0173                                             const AlignmentErrorsExtended* alignmentErrors,
0174                                             const AlignTransform& globalCoordinates,
0175                                             Alignments* newAlignments,
0176                                             AlignmentErrorsExtended* newAlignmentErrorsExtended) {
0177   edm::LogInfo("Alignment") << "@SUB=GeometryAligner::removeGlobalTransform"
0178                             << "Starting to remove global position from alignments and errors";
0179
0180   if (alignments->m_align.size() != alignmentErrors->m_alignError.size())
0181     throw cms::Exception("GeometryMismatch")
0182         << "Size mismatch between alignments (size=" << alignments->m_align.size()
0183         << ") and alignment errors (size=" << alignmentErrors->m_alignError.size() << ")";
0184
0185   const AlignTransform::Translation& globalShift = globalCoordinates.translation();
0186   const AlignTransform::Rotation globalRotation = globalCoordinates.rotation();  // by value!
0187   const AlignTransform::Rotation inverseGlobalRotation = globalRotation.inverse();
0188
0189   AlignTransform::Translation newPosition;
0190   AlignTransform::Rotation newRotation;
0191
0192   std::vector<AlignTransform>::const_iterator iAlign = alignments->m_align.begin();
0193   std::vector<AlignTransformErrorExtended>::const_iterator iAlignError = alignmentErrors->m_alignError.begin();
0194   unsigned int nAPE = 0;
0195   for (iAlign = alignments->m_align.begin(); iAlign != alignments->m_align.end(); ++iAlign, ++iAlignError) {
0196     // Remove global position transformation from alignment
0197     newPosition = inverseGlobalRotation * ((*iAlign).translation() - globalShift);
0198     newRotation = (*iAlign).rotation() * globalRotation;
0199
0200     newAlignments->m_align.emplace_back(AlignTransform(newPosition, newRotation, (*iAlign).rawId()));
0201
0202     // Don't remove global position transformation from APE
0203     // as it wasn't applied. Just fill vector with original
0204     // values
0205     GlobalErrorExtended error(asSMatrix<6>((*iAlignError).matrix()));
0206     newAlignmentErrorsExtended->m_alignError.emplace_back(
0207         AlignTransformErrorExtended((*iAlignError).matrix(), (*iAlignError).rawId()));
0208
0209     //if ( error.cxx() || error.cyy() || error.czz() ||
0210     //   error.cyx() || error.czx() || error.czy() ) {
0211     //      ++nAPE;
0212     //    }
0213
0214     // Code that removes the global postion transformation
0215     // from the APE.
0216     //
0217     //AlgebraicSymMatrix as(3,0);
0218     //as[0][0] = error.cxx();
0219     //as[1][0] = error.cyx(); as[1][1] = error.cyy();
0220     //as[2][0] = error.czx(); as[2][1] = error.czy(); as[2][2] = error.czz();
0221
0222     //AlgebraicMatrix am(3,3);
0223     //am[0][0] = inverseGlobalRotation.xx(); am[0][1] = inverseGlobalRotation.xy(); am[0][2] = inverseGlobalRotation.xz();
0224     //am[1][0] = inverseGlobalRotation.yx(); am[1][1] = inverseGlobalRotation.yy(); am[1][2] = inverseGlobalRotation.yz();
0225     //am[2][0] = inverseGlobalRotation.zx(); am[2][1] = inverseGlobalRotation.zy(); am[2][2] = inverseGlobalRotation.zz();
0226     //as = as.similarityT( am );
0227
0228     //GlobalErrorExtended newError( as );
0229     //newAlignmentErrorsExtended->m_alignError.emplace_back( AlignTransformErrorExtended( newError.matrix(),
0230     //                                                                 (*iAlignError).rawId() ) );
0231     //++nAPE;
0232   }
0233
0234   edm::LogInfo("Alignment") << "@SUB=GeometryAligner::removeGlobalTransform"
0235                             << "Finished to remove global transformation from " << alignments->m_align.size()
0236                             << " alignments with " << nAPE << " non-zero APE.";
0237 }
0238
0239 #endif