TrackerAlignment/plugins/MCMisalignmentScaler.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    Alignment/TrackerAlignment
0004 // Class:      MCMisalignmentScaler
0005 //
0006 /**\class MCMisalignmentScaler MCMisalignmentScaler.cc Alignment/TrackerAlignment/plugins/MCMisalignmentScaler.cc
0007
0008    Description: Plugin to rescale misalignment wrt. ideal geometry
0009
0010    Implementation:
0011
0012    The plugin takes the ideal geometry and the alignment object and rescales
0013    the position difference by the scaling factor provided by the user.
0014
0015 */
0016 //
0017 // Original Author:  Gregor Mittag
0018 //         Created:  Tue, 10 Oct 2017 15:49:17 GMT
0019 //
0020 //
0021
0022 // system include files
0023 #include <memory>
0024 #include <unordered_map>
0025
0026 // user include files
0027 #include "FWCore/Framework/interface/Frameworkfwd.h"
0028 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/EventSetup.h"
0031 #include "FWCore/Framework/interface/ESHandle.h"
0032 #include "FWCore/Framework/interface/MakerMacros.h"
0033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0034 #include "FWCore/ServiceRegistry/interface/Service.h"
0035
0036 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
0037 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
0038
0039 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0040
0041 #include "CondFormats/Alignment/interface/Alignments.h"
0042 #include "CondFormats/AlignmentRecord/interface/TrackerAlignmentRcd.h"
0043 #include "CondFormats/DataRecord/interface/SiPixelQualityRcd.h"
0044 #include "CondFormats/GeometryObjects/interface/PTrackerParameters.h"
0045 #include "CondFormats/SiPixelObjects/interface/SiPixelQuality.h"
0046
0047 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0048 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0049
0050 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0051 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeomBuilderFromGeometricDet.h"
0052 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0053 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0054 #include "Geometry/Records/interface/PTrackerParametersRcd.h"
0055
0056 //
0057 // class declaration
0058 //
0059
0060 class MCMisalignmentScaler : public edm::one::EDAnalyzer<> {
0061 public:
0062   explicit MCMisalignmentScaler(const edm::ParameterSet&);
0063   ~MCMisalignmentScaler() override = default;
0064
0065   static void fillDescriptions(edm::ConfigurationDescriptions&);
0066
0067 private:
0068   const edm::ESGetToken<SiPixelQuality, SiPixelQualityRcd> pixelQualityToken_;
0069   const edm::ESGetToken<SiStripQuality, SiStripQualityRcd> stripQualityToken_;
0070   const edm::ESGetToken<GeometricDet, IdealGeometryRecord> geomDetToken_;
0071   const edm::ESGetToken<PTrackerParameters, PTrackerParametersRcd> ptpToken_;
0072   const edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0073   const edm::ESGetToken<Alignments, TrackerAlignmentRcd> aliToken_;
0074   using ScalerMap = std::unordered_map<unsigned int, std::unordered_map<int, double> >;
0075
0076   void analyze(const edm::Event&, const edm::EventSetup&) override;
0077   ScalerMap decodeSubDetectors(const edm::VParameterSet&);
0078
0079   // ----------member data ---------------------------
0080   const ScalerMap scalers_;
0081   const bool pullBadModulesToIdeal_;
0082   const double outlierPullToIdealCut_;
0083   bool firstEvent_{true};
0084 };
0085
0086 //
0087 // constructors and destructor
0088 //
0089 MCMisalignmentScaler::MCMisalignmentScaler(const edm::ParameterSet& iConfig)
0090     : pixelQualityToken_(esConsumes()),
0091       stripQualityToken_(esConsumes()),
0092       geomDetToken_(esConsumes()),
0093       ptpToken_(esConsumes()),
0094       topoToken_(esConsumes()),
0095       aliToken_(esConsumes()),
0096       scalers_{decodeSubDetectors(iConfig.getParameter<edm::VParameterSet>("scalers"))},
0097       pullBadModulesToIdeal_{iConfig.getUntrackedParameter<bool>("pullBadModulesToIdeal")},
0098       outlierPullToIdealCut_{iConfig.getUntrackedParameter<double>("outlierPullToIdealCut")} {}
0099
0100 //
0101 // member functions
0102 //
0103
0104 // ------------ method called for each event  ------------
0105 void MCMisalignmentScaler::analyze(const edm::Event&, const edm::EventSetup& iSetup) {
0106   if (!firstEvent_)
0107     return;
0108   firstEvent_ = false;
0109
0110   // get handle on bad modules
0111   const SiPixelQuality* pixelModules = &iSetup.getData(pixelQualityToken_);
0112   const SiStripQuality* stripModules = &iSetup.getData(stripQualityToken_);
0113
0114   // get the tracker geometry
0115   const GeometricDet* geometricDet = &iSetup.getData(geomDetToken_);
0116   const PTrackerParameters& ptp = iSetup.getData(ptpToken_);
0117   const TrackerTopology* topology = &iSetup.getData(topoToken_);
0118
0119   TrackerGeomBuilderFromGeometricDet trackerBuilder;
0120   auto tracker = std::unique_ptr<TrackerGeometry>{trackerBuilder.build(geometricDet, ptp, topology)};
0121
0122   auto dets = tracker->dets();
0123   std::sort(dets.begin(), dets.end(), [](const auto& a, const auto& b) {
0124     return a->geographicalId().rawId() < b->geographicalId().rawId();
0125   });
0126
0127   // get the input alignment
0128   const Alignments* alignments = &iSetup.getData(aliToken_);
0129
0130   if (dets.size() != alignments->m_align.size()) {
0131     throw cms::Exception("GeometryMismatch") << "Size mismatch between alignments (size=" << alignments->m_align.size()
0132                                              << ") and ideal geometry (size=" << dets.size() << ")";
0133   }
0134
0135   Alignments rescaledAlignments{};
0136   {
0137     auto outlierCounter{0};
0138     auto ideal = dets.cbegin();
0139     const auto& ideal_end = dets.cend();
0140     auto misaligned = alignments->m_align.cbegin();
0141     for (; ideal != ideal_end; ++ideal, ++misaligned) {
0142       if ((*ideal)->geographicalId().rawId() != misaligned->rawId()) {
0143         throw cms::Exception("GeometryMismatch") << "Order differs between Dets in alignments ideal geometry.";
0144       }
0145
0146       // determine scale factor
0147       const auto& subDetId = (*ideal)->geographicalId().subdetId();
0148       auto side = topology->side((*ideal)->geographicalId());
0149       if (side == 0) {
0150         switch (subDetId) {
0151           case PixelSubdetector::PixelBarrel:
0152             side = 1;  // both sides are treated identical -> pick one of them
0153             break;
0154           case StripSubdetector::TIB:
0155             side = topology->tibSide((*ideal)->geographicalId());
0156             break;
0157           case StripSubdetector::TOB:
0158             side = topology->tobSide((*ideal)->geographicalId());
0159             break;
0160           default:
0161             break;
0162         }
0163       }
0164       auto scaleFactor = scalers_.find(subDetId)->second.find(side)->second;
0165
0166       if (pullBadModulesToIdeal_ &&
0167           (pixelModules->IsModuleBad(misaligned->rawId()) || stripModules->IsModuleBad(misaligned->rawId()))) {
0168         scaleFactor = 0.0;
0169       }
0170
0171       auto x_diff = misaligned->translation().x() - (*ideal)->position().x();
0172       auto y_diff = misaligned->translation().y() - (*ideal)->position().y();
0173       auto z_diff = misaligned->translation().z() - (*ideal)->position().z();
0174
0175       auto xx_diff = misaligned->rotation().xx() - (*ideal)->rotation().xx();
0176       auto xy_diff = misaligned->rotation().xy() - (*ideal)->rotation().xy();
0177       auto xz_diff = misaligned->rotation().xz() - (*ideal)->rotation().xz();
0178       auto yx_diff = misaligned->rotation().yx() - (*ideal)->rotation().yx();
0179       auto yy_diff = misaligned->rotation().yy() - (*ideal)->rotation().yy();
0180       auto yz_diff = misaligned->rotation().yz() - (*ideal)->rotation().yz();
0181       auto zx_diff = misaligned->rotation().zx() - (*ideal)->rotation().zx();
0182       auto zy_diff = misaligned->rotation().zy() - (*ideal)->rotation().zy();
0183       auto zz_diff = misaligned->rotation().zz() - (*ideal)->rotation().zz();
0184
0185       if (outlierPullToIdealCut_ > 0.0 &&
0186           (x_diff * x_diff + y_diff * y_diff + z_diff * z_diff) > outlierPullToIdealCut_ * outlierPullToIdealCut_) {
0187         ++outlierCounter;
0188         edm::LogInfo("Alignment") << outlierCounter << ") Outlier found in subdetector " << subDetId
0189                                   << ":  delta x: " << x_diff << ",  delta y: " << y_diff << ",  delta z: " << z_diff
0190                                   << ",  delta xx: " << xx_diff << ",  delta xy: " << xy_diff
0191                                   << ",  delta xz: " << xz_diff << ",  delta yx: " << yx_diff
0192                                   << ",  delta yx: " << yy_diff << ",  delta yy: " << yz_diff
0193                                   << ",  delta zz: " << zx_diff << ",  delta zy: " << zy_diff
0194                                   << ",  delta zz: " << zz_diff << "\n";
0195         scaleFactor = 0.0;
0196       }
0197
0198       const AlignTransform::Translation rescaledTranslation{(*ideal)->position().x() + scaleFactor * x_diff,
0199                                                             (*ideal)->position().y() + scaleFactor * y_diff,
0200                                                             (*ideal)->position().z() + scaleFactor * z_diff};
0201
0202       const AlignTransform::Rotation rescaledRotation{
0203           CLHEP::HepRep3x3{(*ideal)->rotation().xx() + scaleFactor * xx_diff,
0204                            (*ideal)->rotation().xy() + scaleFactor * xy_diff,
0205                            (*ideal)->rotation().xz() + scaleFactor * xz_diff,
0206                            (*ideal)->rotation().yx() + scaleFactor * yx_diff,
0207                            (*ideal)->rotation().yy() + scaleFactor * yy_diff,
0208                            (*ideal)->rotation().yz() + scaleFactor * yz_diff,
0209                            (*ideal)->rotation().zx() + scaleFactor * zx_diff,
0210                            (*ideal)->rotation().zy() + scaleFactor * zy_diff,
0211                            (*ideal)->rotation().zz() + scaleFactor * zz_diff}};
0212
0213       const AlignTransform rescaledTransform{rescaledTranslation, rescaledRotation, misaligned->rawId()};
0214       rescaledAlignments.m_align.emplace_back(rescaledTransform);
0215     }
0216   }
0217
0218   edm::Service<cond::service::PoolDBOutputService> poolDb;
0219   if (!poolDb.isAvailable()) {
0220     throw cms::Exception("NotAvailable") << "PoolDBOutputService not available";
0221   }
0222   edm::LogInfo("Alignment") << "Writing rescaled tracker-alignment record.";
0223   const auto& since = cond::timeTypeSpecs[cond::runnumber].beginValue;
0224   poolDb->writeOneIOV(rescaledAlignments, since, "TrackerAlignmentRcd");
0225 }
0226
0227 MCMisalignmentScaler::ScalerMap MCMisalignmentScaler::decodeSubDetectors(const edm::VParameterSet& psets) {
0228   // initialize scaler map
0229   ScalerMap subDetMap;
0230   for (unsigned int subDetId = 1; subDetId <= 6; ++subDetId) {
0231     subDetMap[subDetId][1] = 1.0;
0232     subDetMap[subDetId][2] = 1.0;
0233   }
0234
0235   // apply scale factors from configuration
0236   for (const auto& pset : psets) {
0237     const auto& name = pset.getUntrackedParameter<std::string>("subDetector");
0238     const auto& factor = pset.getUntrackedParameter<double>("factor");
0239
0240     std::vector<int> sides;
0241     if (name.find('-') != std::string::npos)
0242       sides.push_back(1);
0243     if (name.find('+') != std::string::npos)
0244       sides.push_back(2);
0245     if (sides.empty()) {  // -> use both sides
0246       sides.push_back(1);
0247       sides.push_back(2);
0248     }
0249
0250     if (name.find("Tracker") != std::string::npos) {
0251       for (unsigned int subDetId = 1; subDetId <= 6; ++subDetId) {
0252         for (const auto& side : sides)
0253           subDetMap[subDetId][side] *= factor;
0254       }
0255       if (sides.size() == 1) {
0256         // if only one side to be scaled
0257         // -> scale also the other side for PXB (subdetid = 1)
0258         subDetMap[PixelSubdetector::PixelBarrel][std::abs(sides[0] - 2) + 1] *= factor;
0259       }
0260     } else if (name.find("PXB") != std::string::npos) {
0261       // ignore sides for PXB
0262       subDetMap[PixelSubdetector::PixelBarrel][1] *= factor;
0263       subDetMap[PixelSubdetector::PixelBarrel][2] *= factor;
0264     } else if (name.find("PXF") != std::string::npos) {
0265       for (const auto& side : sides)
0266         subDetMap[PixelSubdetector::PixelEndcap][side] *= factor;
0267     } else if (name.find("TIB") != std::string::npos) {
0268       for (const auto& side : sides)
0269         subDetMap[StripSubdetector::TIB][side] *= factor;
0270     } else if (name.find("TOB") != std::string::npos) {
0271       for (const auto& side : sides)
0272         subDetMap[StripSubdetector::TOB][side] *= factor;
0273     } else if (name.find("TID") != std::string::npos) {
0274       for (const auto& side : sides)
0275         subDetMap[StripSubdetector::TID][side] *= factor;
0276     } else if (name.find("TEC") != std::string::npos) {
0277       for (const auto& side : sides)
0278         subDetMap[StripSubdetector::TEC][side] *= factor;
0279     } else {
0280       throw cms::Exception("BadConfig") << "@SUB=MCMisalignmentScaler::decodeSubDetectors\n"
0281                                         << "Unknown tracker subdetector: " << name
0282                                         << "\nSupported options: Tracker, PXB, PXF, TIB, TOB, TID, TEC "
0283                                         << "(possibly decorated with '+' or '-')";
0284     }
0285   }
0286
0287   std::stringstream logInfo;
0288   logInfo << "MC misalignment scale factors:\n";
0289   for (const auto& subdet : subDetMap) {
0290     logInfo << "  Subdet " << subdet.first << "\n";
0291     for (const auto& side : subdet.second) {
0292       logInfo << "    side " << side.first << ": " << side.second << "\n";
0293     }
0294     logInfo << "\n";
0295   }
0296   edm::LogInfo("Alignment") << logInfo.str();
0297
0298   return subDetMap;
0299 }
0300
0301 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0302 void MCMisalignmentScaler::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0303   edm::ParameterSetDescription desc;
0304   desc.setComment(
0305       "Creates rescaled MC misalignment scenario. "
0306       "PoolDBOutputService must be set up for 'TrackerAlignmentRcd'.");
0307   edm::ParameterSetDescription descScaler;
0308   descScaler.setComment(
0309       "ParameterSet specifying the tracker part to be scaled "
0310       "by a given factor.");
0311   descScaler.addUntracked<std::string>("subDetector", "Tracker");
0312   descScaler.addUntracked<double>("factor", 1.0);
0313   desc.addVPSet("scalers", descScaler, std::vector<edm::ParameterSet>(1));
0314   desc.addUntracked<bool>("pullBadModulesToIdeal", false);
0315   desc.addUntracked<double>("outlierPullToIdealCut", -1.0);
0316   descriptions.add("mcMisalignmentScaler", desc);
0317 }
0318
0319 //define this as a plug-in
0320 DEFINE_FWK_MODULE(MCMisalignmentScaler);