TrackerAlignment/plugins/TrackerSystematicMisalignments.cc

0001 /** \class TrackerSystematicMisalignments
0002  *
0003  *  Class to misaligned tracker from DB.
0004  *
0005  *  $Date: 2012/06/13 09:24:50 $
0006  *  $Revision: 1.5 $
0007  *  \author Chung Khim Lae
0008  */
0009
0010 // user include files
0011 #include "Alignment/CommonAlignment/interface/Alignable.h"
0012 #include "Alignment/CommonAlignment/interface/SurveyDet.h"
0013 #include "Alignment/TrackerAlignment/interface/AlignableTracker.h"
0014 #include "CLHEP/Random/RandGauss.h"
0015 #include "CondFormats/Alignment/interface/AlignmentErrorsExtended.h"
0016 #include "CondFormats/Alignment/interface/Alignments.h"
0017 #include "CondFormats/Alignment/interface/DetectorGlobalPosition.h"
0018 #include "CondFormats/AlignmentRecord/interface/GlobalPositionRcd.h"
0019 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentErrorExtendedRcd.h"
0020 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentRcd.h"
0021 #include "CondFormats/GeometryObjects/interface/PTrackerParameters.h"
0022 #include "DataFormats/DetId/interface/DetId.h"
0023 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"  // for enums TID/TIB/etc.
0024 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028 #include "Geometry/CommonTopologies/interface/GeometryAligner.h"
0029 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0030 #include "Geometry/Records/interface/PTrackerParametersRcd.h"
0031 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0032 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0033 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeomBuilderFromGeometricDet.h"
0034 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0035
0036 // Database
0037 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0038 #include "FWCore/ServiceRegistry/interface/Service.h"
0039
0040 class AlignableSurface;
0041 class Alignments;
0042
0043 namespace edm {
0044   class ParameterSet;
0045 }
0046
0047 class TrackerSystematicMisalignments : public edm::one::EDAnalyzer<> {
0048 public:
0049   TrackerSystematicMisalignments(const edm::ParameterSet&);
0050
0051   /// Read ideal tracker geometry from DB
0052   void beginJob() override;
0053
0054   void analyze(const edm::Event&, const edm::EventSetup&) override;
0055
0056 private:
0057   void applySystematicMisalignment(Alignable*);
0058   //align::GlobalVector findSystematicMis( align::PositionType );
0059   align::GlobalVector findSystematicMis(const align::PositionType&, const bool blindToZ, const bool blindToR);
0060
0061   const edm::ESGetToken<GeometricDet, IdealGeometryRecord> geomDetToken_;
0062   const edm::ESGetToken<PTrackerParameters, PTrackerParametersRcd> ptpToken_;
0063   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0064   const edm::ESGetToken<Alignments, TrackerAlignmentRcd> aliToken_;
0065   const edm::ESGetToken<AlignmentErrorsExtended, TrackerAlignmentErrorExtendedRcd> aliErrorToken_;
0066   const edm::ESGetToken<Alignments, GlobalPositionRcd> gprToken_;
0067   AlignableTracker* theAlignableTracker;
0068
0069   // configurables needed for the systematic misalignment
0070   bool m_fromDBGeom;
0071
0072   double m_radialEpsilon;
0073   double m_telescopeEpsilon;
0074   double m_layerRotEpsilon;
0075   double m_bowingEpsilon;
0076   double m_zExpEpsilon;
0077   double m_twistEpsilon;
0078   double m_ellipticalEpsilon;
0079   double m_skewEpsilon;
0080   double m_sagittaEpsilon;
0081
0082   //misalignment phases
0083   double m_ellipticalDelta;
0084   double m_skewDelta;
0085   double m_sagittaDelta;
0086
0087   // flag to steer suppression of blind movements
0088   bool suppressBlindMvmts;
0089
0090   // flag for old z behaviour, version <= 1.5
0091   bool oldMinusZconvention;
0092 };
0093
0094 // -----------------------------------------------------------------
0095 // 2010-05-20 Frank Meier
0096 // Changed sign of z-correction, i.e. z-expansion is now an expansion
0097 // made some variables constant, removed obviously dead code and comments
0098
0099 TrackerSystematicMisalignments::TrackerSystematicMisalignments(const edm::ParameterSet& cfg)
0100     : geomDetToken_(esConsumes()),
0101       ptpToken_(esConsumes()),
0102       topoToken_(esConsumes()),
0103       aliToken_(esConsumes()),
0104       aliErrorToken_(esConsumes()),
0105       gprToken_(esConsumes()),
0106       theAlignableTracker(nullptr) {
0107   // use existing geometry
0108   m_fromDBGeom = cfg.getUntrackedParameter<bool>("fromDBGeom");
0109
0110   // constants
0111   m_radialEpsilon = cfg.getUntrackedParameter<double>("radialEpsilon");
0112   m_telescopeEpsilon = cfg.getUntrackedParameter<double>("telescopeEpsilon");
0113   m_layerRotEpsilon = cfg.getUntrackedParameter<double>("layerRotEpsilon");
0114   m_bowingEpsilon = cfg.getUntrackedParameter<double>("bowingEpsilon");
0115   m_zExpEpsilon = cfg.getUntrackedParameter<double>("zExpEpsilon");
0116   m_twistEpsilon = cfg.getUntrackedParameter<double>("twistEpsilon");
0117   m_ellipticalEpsilon = cfg.getUntrackedParameter<double>("ellipticalEpsilon");
0118   m_skewEpsilon = cfg.getUntrackedParameter<double>("skewEpsilon");
0119   m_sagittaEpsilon = cfg.getUntrackedParameter<double>("sagittaEpsilon");
0120
0121   m_ellipticalDelta = cfg.getUntrackedParameter<double>("ellipticalDelta");
0122   m_skewDelta = cfg.getUntrackedParameter<double>("skewDelta");
0123   m_sagittaDelta = cfg.getUntrackedParameter<double>("sagittaDelta");
0124
0125   if (m_radialEpsilon > -990.0) {
0126     edm::LogWarning("MisalignedTracker") << "Applying radial ...";
0127   }
0128   if (m_telescopeEpsilon > -990.0) {
0129     edm::LogWarning("MisalignedTracker") << "Applying telescope ...";
0130   }
0131   if (m_layerRotEpsilon > -990.0) {
0132     edm::LogWarning("MisalignedTracker") << "Applying layer rotation ...";
0133   }
0134   if (m_bowingEpsilon > -990.0) {
0135     edm::LogWarning("MisalignedTracker") << "Applying bowing ...";
0136   }
0137   if (m_zExpEpsilon > -990.0) {
0138     edm::LogWarning("MisalignedTracker") << "Applying z-expansion ...";
0139   }
0140   if (m_twistEpsilon > -990.0) {
0141     edm::LogWarning("MisalignedTracker") << "Applying twist ...";
0142   }
0143   if (m_ellipticalEpsilon > -990.0) {
0144     edm::LogWarning("MisalignedTracker") << "Applying elliptical ...";
0145   }
0146   if (m_skewEpsilon > -990.0) {
0147     edm::LogWarning("MisalignedTracker") << "Applying skew ...";
0148   }
0149   if (m_sagittaEpsilon > -990.0) {
0150     edm::LogWarning("MisalignedTracker") << "Applying sagitta ...";
0151   }
0152
0153   // get flag for suppression of blind movements
0154   suppressBlindMvmts = cfg.getUntrackedParameter<bool>("suppressBlindMvmts");
0155   if (suppressBlindMvmts) {
0156     edm::LogWarning("MisalignedTracker") << "Blind movements suppressed (TIB/TOB in z, TID/TEC in r)";
0157   }
0158
0159   // compatibility with old (weird) z convention
0160   oldMinusZconvention = cfg.getUntrackedParameter<bool>("oldMinusZconvention");
0161   if (oldMinusZconvention) {
0162     edm::LogWarning("MisalignedTracker") << "Old z convention: dz --> -dz";
0163   } else {
0164     edm::LogWarning("MisalignedTracker") << "New z convention: dz --> dz";
0165   }
0166 }
0167
0168 void TrackerSystematicMisalignments::beginJob() {}
0169
0170 void TrackerSystematicMisalignments::analyze(const edm::Event& event, const edm::EventSetup& setup) {
0171   //Retrieve tracker topology from geometry
0172   const GeometricDet* geom = &setup.getData(geomDetToken_);
0173   const PTrackerParameters& ptp = setup.getData(ptpToken_);
0174   const TrackerTopology* tTopo = &setup.getData(topoToken_);
0175
0176   TrackerGeometry* tracker = TrackerGeomBuilderFromGeometricDet().build(geom, ptp, tTopo);
0177
0178   //take geometry from DB or randomly generate geometry
0179   if (m_fromDBGeom) {
0180     //build the tracker
0181     const Alignments* alignments = &setup.getData(aliToken_);
0182     const AlignmentErrorsExtended* alignmentErrors = &setup.getData(aliErrorToken_);
0183     const Alignments* globalPositionRcd = &setup.getData(gprToken_);
0184
0185     //apply the latest alignments
0186     GeometryAligner aligner;
0187     aligner.applyAlignments<TrackerGeometry>(&(*tracker),
0188                                              &(*alignments),
0189                                              &(*alignmentErrors),
0190                                              align::DetectorGlobalPosition(*globalPositionRcd, DetId(DetId::Tracker)));
0191   }
0192
0193   theAlignableTracker = new AlignableTracker(&(*tracker), tTopo);
0194
0195   applySystematicMisalignment(&(*theAlignableTracker));
0196
0197   // -------------- writing out to alignment record --------------
0198   Alignments myAlignments = *(theAlignableTracker->alignments());
0199   AlignmentErrorsExtended myAlignmentErrorsExtended = *(theAlignableTracker->alignmentErrors());
0200
0201   // Store alignment[Error]s to DB
0202   edm::Service<cond::service::PoolDBOutputService> poolDbService;
0203   std::string theAlignRecordName = "TrackerAlignmentRcd";
0204   std::string theErrorRecordName = "TrackerAlignmentErrorExtendedRcd";
0205
0206   // Call service
0207   if (!poolDbService.isAvailable())  // Die if not available
0208     throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
0209
0210   poolDbService->writeOneIOV<Alignments>(myAlignments, poolDbService->beginOfTime(), theAlignRecordName);
0211   poolDbService->writeOneIOV<AlignmentErrorsExtended>(
0212       myAlignmentErrorsExtended, poolDbService->beginOfTime(), theErrorRecordName);
0213 }
0214
0215 void TrackerSystematicMisalignments::applySystematicMisalignment(Alignable* ali) {
0216   const align::Alignables& comp = ali->components();
0217   unsigned int nComp = comp.size();
0218   //move then do for lower level object
0219   //for issue of det vs detunit
0220   bool usecomps = true;
0221   if ((ali->alignableObjectId() == 2) && (nComp >= 1))
0222     usecomps = false;
0223   for (unsigned int i = 0; i < nComp; ++i) {
0224     if (usecomps)
0225       applySystematicMisalignment(comp[i]);
0226   }
0227
0228   // if suppression of blind mvmts: check if subdet is blind to a certain mode
0229   bool blindToZ(false), blindToR(false);
0230   if (suppressBlindMvmts) {
0231     const int subdetid = ali->geomDetId().subdetId();
0232     switch (subdetid) {
0233       // TIB/TON blind to z
0234       case SiStripDetId::TIB:
0235       case SiStripDetId::TOB:
0236         blindToZ = true;
0237         break;
0238       // TID/TEC blind to R
0239       case SiStripDetId::TID:
0240       case SiStripDetId::TEC:
0241         blindToR = true;
0242         break;
0243       default:
0244         break;
0245     }
0246   }
0247
0248   const int level = ali->alignableObjectId();
0249   if ((level == 1) || (level == 2)) {
0250     const align::PositionType gP = ali->globalPosition();
0251     const align::GlobalVector gVec = findSystematicMis(gP, blindToZ, blindToR);
0252     ali->move(gVec);
0253   }
0254 }
0255
0256 align::GlobalVector TrackerSystematicMisalignments::findSystematicMis(const align::PositionType& globalPos,
0257                                                                       const bool blindToZ,
0258                                                                       const bool blindToR) {
0259   //align::GlobalVector TrackerSystematicMisalignments::findSystematicMis( align::PositionType globalPos ){
0260   // calculates shift for the current alignable
0261   // all corrections are calculated w.r.t. the original geometry
0262   double deltaX = 0.0;
0263   double deltaY = 0.0;
0264   double deltaZ = 0.0;
0265   const double oldX = globalPos.x();
0266   const double oldY = globalPos.y();
0267   const double oldZ = globalPos.z();
0268   const double oldPhi = globalPos.phi();
0269   const double oldR = sqrt(globalPos.x() * globalPos.x() + globalPos.y() * globalPos.y());
0270
0271   if (m_radialEpsilon > -990.0 && !blindToR) {
0272     deltaX += m_radialEpsilon * oldX;
0273     deltaY += m_radialEpsilon * oldY;
0274   }
0275   if (m_telescopeEpsilon > -990.0 && !blindToZ) {
0276     deltaZ += m_telescopeEpsilon * oldR;
0277   }
0278   if (m_layerRotEpsilon > -990.0) {
0279     // The following number was chosen such that the Layer Rotation systematic
0280     // misalignment would not cause an overall rotation of the tracker.
0281     const double Roffset = 57.0;
0282     const double xP = oldR * cos(oldPhi + m_layerRotEpsilon * (oldR - Roffset));
0283     const double yP = oldR * sin(oldPhi + m_layerRotEpsilon * (oldR - Roffset));
0284     deltaX += (xP - oldX);
0285     deltaY += (yP - oldY);
0286   }
0287   if (m_bowingEpsilon > -990.0 && !blindToR) {
0288     const double trackeredgePlusZ = 271.846;
0289     const double bowfactor = m_bowingEpsilon * (trackeredgePlusZ * trackeredgePlusZ - oldZ * oldZ);
0290     deltaX += oldX * bowfactor;
0291     deltaY += oldY * bowfactor;
0292   }
0293   if (m_zExpEpsilon > -990.0 && !blindToZ) {
0294     deltaZ += oldZ * m_zExpEpsilon;
0295   }
0296   if (m_twistEpsilon > -990.0) {
0297     const double xP = oldR * cos(oldPhi + m_twistEpsilon * oldZ);
0298     const double yP = oldR * sin(oldPhi + m_twistEpsilon * oldZ);
0299     deltaX += (xP - oldX);
0300     deltaY += (yP - oldY);
0301   }
0302   if (m_ellipticalEpsilon > -990.0 && !blindToR) {
0303     deltaX += oldX * m_ellipticalEpsilon * cos(2.0 * oldPhi + m_ellipticalDelta);
0304     deltaY += oldY * m_ellipticalEpsilon * cos(2.0 * oldPhi + m_ellipticalDelta);
0305   }
0306   if (m_skewEpsilon > -990.0 && !blindToZ) {
0307     deltaZ += m_skewEpsilon * cos(oldPhi + m_skewDelta);
0308   }
0309   if (m_sagittaEpsilon > -990.0) {
0310     // deltaX += oldX/fabs(oldX)*m_sagittaEpsilon; // old one...
0311     deltaX += oldR * m_sagittaEpsilon * sin(m_sagittaDelta);
0312     deltaY += oldR * m_sagittaEpsilon * cos(m_sagittaDelta);  //Delta y is cos so that delta=0 reflects the old behavior
0313   }
0314
0315   // Compatibility with old version <= 1.5
0316   if (oldMinusZconvention)
0317     deltaZ = -deltaZ;
0318
0319   align::GlobalVector gV(deltaX, deltaY, deltaZ);
0320   return gV;
0321 }
0322
0323 // Plug in to framework
0324
0325 #include "FWCore/Framework/interface/MakerMacros.h"
0326
0327 DEFINE_FWK_MODULE(TrackerSystematicMisalignments);