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