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