CMS 3D CMS Logo

AlignableSiStripDet.cc

Go to the documentation of this file.
00001 /* 
00002  *  $Date: 2008/07/13 12:26:41 $
00003  *  $Revision: 1.5 $
00004  */
00005 
00006 #include "Alignment/TrackerAlignment/interface/AlignableSiStripDet.h"
00007  
00008 #include "Alignment/CommonAlignment/interface/AlignableSurface.h"
00009 
00010 #include "CondFormats/Alignment/interface/AlignmentErrors.h"
00011 #include "CondFormats/Alignment/interface/AlignTransformError.h"
00012 
00013 #include "DataFormats/GeometrySurface/interface/Bounds.h"
00014 #include "DataFormats/TrackingRecHit/interface/AlignmentPositionError.h"
00015 #include "DataFormats/GeometryCommonDetAlgo/interface/GlobalError.h"
00016 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
00017 #include "Geometry/TrackerGeometryBuilder/interface/PlaneBuilderForGluedDet.h"
00018 #include "Geometry/TrackerGeometryBuilder/interface/GluedGeomDet.h"
00019 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
00020 
00021 #include "FWCore/Utilities/interface/Exception.h"
00022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00023 
00024 #include <math.h>
00025 
00026 AlignableSiStripDet::AlignableSiStripDet(const GluedGeomDet *gluedDet) 
00027   : AlignableDet(gluedDet, true), // true: adding DetUnits
00028     theMonoBounds  (gluedDet->monoDet()  ->surface().bounds().clone()),
00029     theStereoBounds(gluedDet->stereoDet()->surface().bounds().clone()),
00030     theMonoType  (static_cast<const StripGeomDetUnit*>(gluedDet->monoDet())  ->specificType()),
00031     theStereoType(static_cast<const StripGeomDetUnit*>(gluedDet->stereoDet())->specificType())
00032 {
00033   // It is not allowed to store a pointer to GeomDet within objects with a life time
00034   // longer than an Event: 
00035   // GeomDet comes from TrackerGeometry that is created from GeometricDet that depends on
00036   // IdealGeometryRecord from EventSetup, so it could change in next event!
00037   //  ==> Need to store directly what I need from it. 
00038   // Unfortunately the current way with references for the type does not solve that,
00039   // either. But currently no way out, see header file.
00040 
00041   // check order mono/stereo
00042   const Alignables units(this->components());
00043   if (units.size() != 2
00044       || gluedDet->monoDet()->geographicalId() != units[0]->geomDetId()
00045       || gluedDet->stereoDet()->geographicalId() != units[1]->geomDetId()) {
00046     throw cms::Exception("LogicError")
00047       << "[AlignableSiStripDet] " << "Either != 2 components or "
00048       << "mono/stereo in wrong order for consistifyAlignments.";
00049   }
00050 }
00051 
00052 //__________________________________________________________________________________________________
00053 AlignableSiStripDet::~AlignableSiStripDet()
00054 {
00055   delete theMonoBounds;
00056   delete theStereoBounds;
00057 }
00058 
00059 //__________________________________________________________________________________________________
00060 Alignments* AlignableSiStripDet::alignments() const
00061 {
00062   const_cast<AlignableSiStripDet*>(this)->consistifyAlignments();
00063 
00064   return this->AlignableDet::alignments();
00065 }
00066 
00067 //__________________________________________________________________________________________________
00068 AlignmentErrors* AlignableSiStripDet::alignmentErrors() const
00069 {
00070   const_cast<AlignableSiStripDet*>(this)->consistifyAlignmentErrors();
00071 
00072   return this->AlignableDet::alignmentErrors();
00073 }
00074 
00075 //__________________________________________________________________________________________________
00076 void AlignableSiStripDet::consistifyAlignments()
00077 {
00078   // make alignments consistent with daughters, calling method from geometry
00079 
00080   // The aim of all this gymnastics is to have the alignments calculated by
00081   // PlaneBuilderForGluedDet::plane(const std::vector<const GeomDetUnit*> &detComps);
00082   // 
00083   // So we take the (new) position and orientation from the AligableDetUnits,
00084   // but bounds and GeomDetType from original GeomDetUnits to create new GeomDetUnits
00085   // that are passed to that routine.
00086 
00087   const Alignables aliUnits(this->components()); // order mono==0, stereo==1 checked in ctr.
00088 
00089   BoundPlane::BoundPlanePointer monoPlane
00090     = BoundPlane::build(aliUnits[0]->globalPosition(), aliUnits[0]->globalRotation(),
00091                         *theMonoBounds);
00092   // Fortunately we do not seem to need a GeometricDet pointer and can use 0:
00093   const StripGeomDetUnit monoDet(&(*monoPlane), &theMonoType, 0);
00094 
00095   BoundPlane::BoundPlanePointer stereoPlane
00096     = BoundPlane::build(aliUnits[1]->globalPosition(), aliUnits[1]->globalRotation(),
00097                         *theStereoBounds);
00098   // Fortunately we do not seem to need a GeometricDet pointer and can use 0:
00099   const StripGeomDetUnit stereoDet(&(*stereoPlane), &theStereoType, 0);
00100 
00101   std::vector<const GeomDetUnit*> detComps;
00102   detComps.push_back(&monoDet);   // order mono first, stereo second should be as in...
00103   detComps.push_back(&stereoDet); // ...TrackerGeomBuilderFromGeometricDet::buildGeomDet
00104 
00105   // Now we have all to calculate new position and rotation via PlaneBuilderForGluedDet.
00106   const PositionType oldPos(theSurface.position()); // From old surface for keeping...
00107   const RotationType oldRot(theSurface.rotation()); // ...track of changes.
00108 
00109   PlaneBuilderForGluedDet planeBuilder;
00110   theSurface = AlignableSurface(*planeBuilder.plane(detComps));
00111 
00112   // But do not forget to keep track of movements/rotations:
00113   const GlobalVector movement(theSurface.position().basicVector() - oldPos.basicVector());
00114   // Seems to be correct down to delta angles 4.*1e-8:
00115   const RotationType rotation(oldRot.multiplyInverse(theSurface.rotation()));
00116   this->addDisplacement(movement);
00117   this->addRotation(rotation);
00118 
00119 //   this->dumpCompareEuler(oldRot, theSurface.rotation());
00120 
00121 //   if (movement.mag2()) { // > 1.e-10) { 
00122 //     edm::LogWarning("Alignment") << "@SUB=consistifyAlignments" 
00123 //                               << "Delta: " << movement.x() << " " << movement.y() << " " << movement.z()
00124 //                               << "\nPos: " << oldPos.perp() << " " << oldPos.phi() << " " << oldPos.z();
00125 //   }
00126 
00127 }
00128 
00129 
00130 //__________________________________________________________________________________________________
00131 void AlignableSiStripDet::consistifyAlignmentErrors()
00132 {
00133   // make alignment errors consistent with daughters
00134 
00135   AlignmentErrors *oldErrs = this->AlignableDet::alignmentErrors();
00136 
00137   const Alignables units(this->components()); // order mono==0, stereo==1 does not matter here
00138 
00139   const AlignTransformError &gluedErr  = this->errorFromId(oldErrs->m_alignError,
00140                                                            this->geomDetId());
00141   const AlignTransformError &monoErr   = this->errorFromId(oldErrs->m_alignError,
00142                                                            units[0]->geomDetId());
00143   const AlignTransformError &stereoErr = this->errorFromId(oldErrs->m_alignError,
00144                                                            units[1]->geomDetId());
00145   const GlobalError errGlued (gluedErr.matrix());
00146   const GlobalError errMono  (monoErr.matrix());
00147   const GlobalError errStereo(stereoErr.matrix());
00148 
00149   //  const GlobalError newErrGlued((errMono + errStereo - errGlued).matrix_new() /= 4.);
00150   // The above line would be error propagation assuming:
00151   //   - Glued position is just the mean of its components.
00152   //   - Components APE is square sum of what has been set to glued and to components itself.
00153   // But this can be too small, e.g. 
00154   // - Simply by factor sqrt(4)=2 smaller than the old 'errGlued' in case APE (from position)
00155   //   was only set to glued (and 1-to-1 propagated to components).
00156   // - Factor sqrt(2) smaller than components APE in case APE (from position) was only
00157   //   directly applied to both components.
00158   //
00159   // So I choose the max of all three, ignoring correlations (which we do not have?), sigh!
00160   // And in this way it is safe against repetetive calls to this method!
00161   double maxX2 = (errMono.cxx() > errStereo.cxx() ? errMono.cxx() : errStereo.cxx());
00162   maxX2 = (maxX2 > errGlued.cxx() ? maxX2 : errGlued.cxx());
00163   double maxY2 = (errMono.cyy() > errStereo.cyy() ? errMono.cyy() : errStereo.cyy());
00164   maxY2 = (maxY2 > errGlued.cyy() ? maxY2 : errGlued.cyy());
00165   double maxZ2 = (errMono.czz() > errStereo.czz() ? errMono.czz() : errStereo.czz());
00166   maxZ2 = (maxZ2 > errGlued.czz() ? maxZ2 : errGlued.czz());
00167   const AlignmentPositionError newApeGlued(sqrt(maxX2), sqrt(maxY2), sqrt(maxZ2));
00168 
00169   // Now set new errors - and reset those of the components, since they get overwritten...
00170   this->setAlignmentPositionError(newApeGlued);
00171   units[0]->setAlignmentPositionError(AlignmentPositionError(errMono));
00172   units[1]->setAlignmentPositionError(AlignmentPositionError(errStereo));
00173 
00174 //   edm::LogWarning("Alignment") << "@SUB=consistifyAlignmentErrors" 
00175 //                             << "End Id " << this->geomDetId();
00176 //   AlignmentErrors *newErrs = this->AlignableDet::alignmentErrors();
00177 //   this->dumpCompareAPE(oldErrs->m_alignError, newErrs->m_alignError);
00178 //   delete newErrs;
00179 
00180   delete oldErrs;
00181 }
00182 
00183 //__________________________________________________________________________________________________
00184 const AlignTransformError& 
00185 AlignableSiStripDet::errorFromId(const std::vector<AlignTransformError> &trafoErrs,
00186                                  align::ID id) const
00187 {
00188   for (unsigned int i = 0; i < trafoErrs.size(); ++i) {
00189     if (trafoErrs[i].rawId() ==  id) return trafoErrs[i];
00190   }
00191 
00192   throw cms::Exception("Mismatch") << "[AlignableSiStripDet::indexFromId] "
00193                                    << id << " not found.";
00194 
00195   return trafoErrs.front(); // never reached due to exception (but pleasing the compiler)
00196 }
00197 
00198 // //__________________________________________________________________________________________________
00199 // void AlignableSiStripDet::dumpCompareAPE(const std::vector<AlignTransformError> &trafoErrs1,
00200 //                                       const std::vector<AlignTransformError> &trafoErrs2) const
00201 // {
00202 //   for (unsigned int i = 0; i < trafoErrs1.size() && i < trafoErrs2.size(); ++i) {
00203 //     if (trafoErrs1[i].rawId() != trafoErrs2[i].rawId()) {
00204 //       // complain
00205 //       break;
00206 //     }
00207 //     const GlobalError globErr1(trafoErrs1[i].matrix());
00208 //     const GlobalError globErr2(trafoErrs2[i].matrix());
00209 //     edm::LogVerbatim("Alignment") << trafoErrs1[i].rawId() << " | " 
00210 //                                << globErr1.cxx() << " " 
00211 //                                << globErr1.cyy() << " " 
00212 //                                << globErr1.czz() << " | "
00213 //                                << globErr2.cxx() << " " 
00214 //                                << globErr2.cyy() << " " 
00215 //                                << globErr2.czz();
00216 //   }
00217 // }
00218 //
00219 // #include "CLHEP/Vector/EulerAngles.h"
00220 // #include "CLHEP/Vector/Rotation.h"
00221 // //__________________________________________________________________________________________________
00222 // void AlignableSiStripDet::dumpCompareEuler(const RotationType &oldRot,
00223 //                                         const RotationType &newRot) const
00224 // {
00225 //   // 
00226 //   const HepRotation oldClhep(HepRep3x3(oldRot.xx(), oldRot.xy(), oldRot.xz(),
00227 //                                     oldRot.yx(), oldRot.yy(), oldRot.yz(),
00228 //                                     oldRot.zx(), oldRot.zy(), oldRot.zz()));
00229 //
00230 //   const HepRotation newClhep(HepRep3x3(newRot.xx(), newRot.xy(), newRot.xz(),
00231 //                                     newRot.yx(), newRot.yy(), newRot.yz(),
00232 //                                     newRot.zx(), newRot.zy(), newRot.zz()));
00233 //
00234 //   const RotationType rotationGlob(oldRot.multiplyInverse(newRot));
00235 //   const RotationType rotation(theSurface.toLocal(rotationGlob)); // not 100% correct: new global...
00236 //   const HepRotation diff(HepRep3x3(rotation.xx(), rotation.xy(), rotation.xz(),
00237 //                                     rotation.yx(), rotation.yy(), rotation.yz(),
00238 //                                     rotation.zx(), rotation.zy(), rotation.zz()));
00239 //
00240 //   edm::LogWarning("Alignment") << "@SUB=dumpCompareEuler" 
00241 //                             << "oldEuler " << oldClhep.eulerAngles()
00242 //                             << "\nnewEuler " << newClhep.eulerAngles()
00243 //                             << "\n diff_euler " << diff.eulerAngles()
00244 //                             << "\n diff_diag (" << diff.xx() << ", " << diff.yy() 
00245 //                             <<         ", " << diff.zz() << ")";
00246 // }

Generated on Tue Jun 9 17:25:02 2009 for CMSSW by  doxygen 1.5.4