File indexing completed on 2023-01-21 00:19:08
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023 #include <memory>
0024 #include <unordered_map>
0025
0026
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
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
0080 const ScalerMap scalers_;
0081 const bool pullBadModulesToIdeal_;
0082 const double outlierPullToIdealCut_;
0083 bool firstEvent_{true};
0084 };
0085
0086
0087
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
0102
0103
0104
0105 void MCMisalignmentScaler::analyze(const edm::Event&, const edm::EventSetup& iSetup) {
0106 if (!firstEvent_)
0107 return;
0108 firstEvent_ = false;
0109
0110
0111 const SiPixelQuality* pixelModules = &iSetup.getData(pixelQualityToken_);
0112 const SiStripQuality* stripModules = &iSetup.getData(stripQualityToken_);
0113
0114
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
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
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;
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
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
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()) {
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
0257
0258 subDetMap[PixelSubdetector::PixelBarrel][std::abs(sides[0] - 2) + 1] *= factor;
0259 }
0260 } else if (name.find("PXB") != std::string::npos) {
0261
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
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
0320 DEFINE_FWK_MODULE(MCMisalignmentScaler);