DTGeometryBuilder/plugins/DTGeometryValidate.cc

0001 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0002 #include "FWCore/Framework/interface/MakerMacros.h"
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007
0008 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0009 #include "Geometry/DTGeometry/interface/DTLayer.h"
0010 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0011
0012 #include "Fireworks/Core/interface/FWGeometry.h"
0013
0014 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0015 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0016
0017 #include <TFile.h>
0018 #include <TH1.h>
0019
0020 #include <cmath>
0021 #include <limits>
0022 #include <string>
0023 #include <type_traits>
0024 #include <algorithm>
0025
0026 using namespace std;
0027
0028 template <class T>
0029 typename enable_if<!numeric_limits<T>::is_integer, bool>::type almost_equal(T x, T y, int ulp) {
0030   // the machine epsilon has to be scaled to the magnitude of the values used
0031   // and multiplied by the desired precision in ULPs (units in the last place)
0032   return abs(x - y) <= numeric_limits<T>::epsilon() * abs(x + y) * ulp
0033          // unless the result is subnormal
0034          || abs(x - y) < numeric_limits<T>::min();
0035 }
0036
0037 using namespace edm;
0038
0039 class DTGeometryValidate : public one::EDAnalyzer<> {
0040 public:
0041   explicit DTGeometryValidate(const ParameterSet&);
0042   ~DTGeometryValidate() override {}
0043
0044 private:
0045   void beginJob() override;
0046   void analyze(const edm::Event&, const edm::EventSetup&) override;
0047   void endJob() override;
0048
0049   void validateDTChamberGeometry();
0050   void validateDTLayerGeometry();
0051
0052   void compareTransform(const GlobalPoint&, const TGeoMatrix*);
0053   void compareShape(const GeomDet*, const float*);
0054
0055   float getDistance(const GlobalPoint&, const GlobalPoint&);
0056   float getDiff(const float, const float);
0057
0058   void makeHistograms(const char*);
0059   void makeHistogram(const string&, vector<float>&);
0060
0061   void clearData() {
0062     globalDistances_.clear();
0063     topWidths_.clear();
0064     bottomWidths_.clear();
0065     lengths_.clear();
0066     thicknesses_.clear();
0067   }
0068
0069   edm::ESGetToken<DTGeometry, MuonGeometryRecord> dtGeometryToken_;
0070   edm::ESHandle<DTGeometry> dtGeometry_;
0071   FWGeometry fwGeometry_;
0072   TFile* outFile_;
0073   vector<float> globalDistances_;
0074   vector<float> topWidths_;
0075   vector<float> bottomWidths_;
0076   vector<float> lengths_;
0077   vector<float> thicknesses_;
0078   string infileName_;
0079   string outfileName_;
0080   int tolerance_;
0081 };
0082
0083 DTGeometryValidate::DTGeometryValidate(const edm::ParameterSet& iConfig)
0084     : dtGeometryToken_{esConsumes<DTGeometry, MuonGeometryRecord>(edm::ESInputTag{})},
0085       infileName_(
0086           iConfig.getUntrackedParameter<string>("infileName", "Geometry/DTGeometryBuilder/data/cmsRecoGeom-2021.root")),
0087       outfileName_(iConfig.getUntrackedParameter<string>("outfileName", "validateDTGeometry.root")),
0088       tolerance_(iConfig.getUntrackedParameter<int>("tolerance", 6)) {
0089   edm::FileInPath fp(infileName_.c_str());
0090   fwGeometry_.loadMap(fp.fullPath().c_str());
0091   outFile_ = TFile::Open(outfileName_.c_str(), "RECREATE");
0092 }
0093
0094 void DTGeometryValidate::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
0095   dtGeometry_ = eventSetup.getHandle(dtGeometryToken_);
0096
0097   if (dtGeometry_.isValid()) {
0098     LogVerbatim("DTGeometry") << "Validating DT chamber geometry";
0099     validateDTChamberGeometry();
0100     LogVerbatim("DTGeometry") << "Validating DT layer geometry";
0101     validateDTLayerGeometry();
0102   } else
0103     LogVerbatim("DTGeometry") << "Invalid DT geometry";
0104 }
0105
0106 void DTGeometryValidate::validateDTChamberGeometry() {
0107   clearData();
0108
0109   for (auto const& it : dtGeometry_->chambers()) {
0110     DTChamberId chId = it->id();
0111     GlobalPoint gp = it->surface().toGlobal(LocalPoint(0.0, 0.0, 0.0));
0112
0113     const TGeoMatrix* matrix = fwGeometry_.getMatrix(chId.rawId());
0114
0115     if (!matrix) {
0116       LogVerbatim("DTGeometry") << "Failed to get matrix of DT chamber with detid: " << chId.rawId();
0117       continue;
0118     }
0119
0120     compareTransform(gp, matrix);
0121
0122     auto const& shape = fwGeometry_.getShapePars(chId.rawId());
0123
0124     if (!shape) {
0125       LogVerbatim("DTGeometry") << "Failed to get shape of DT chamber with detid: " << chId.rawId();
0126       continue;
0127     }
0128
0129     compareShape(it, shape);
0130   }
0131
0132   makeHistograms("DT Chamber");
0133 }
0134
0135 void DTGeometryValidate::validateDTLayerGeometry() {
0136   clearData();
0137
0138   vector<float> wire_positions;
0139
0140   for (auto const& it : dtGeometry_->layers()) {
0141     DTLayerId layerId = it->id();
0142     GlobalPoint gp = it->surface().toGlobal(LocalPoint(0.0, 0.0, 0.0));
0143
0144     const TGeoMatrix* matrix = fwGeometry_.getMatrix(layerId.rawId());
0145
0146     if (!matrix) {
0147       LogVerbatim("DTGeometry") << "Failed to get matrix of DT layer with detid: " << layerId.rawId();
0148       continue;
0149     }
0150
0151     compareTransform(gp, matrix);
0152
0153     auto const& shape = fwGeometry_.getShapePars(layerId.rawId());
0154
0155     if (!shape) {
0156       LogVerbatim("DTGeometry") << "Failed to get shape of DT layer with detid: " << layerId.rawId();
0157       continue;
0158     }
0159
0160     compareShape(it, shape);
0161
0162     auto const& parameters = fwGeometry_.getParameters(layerId.rawId());
0163
0164     if (parameters == nullptr) {
0165       LogVerbatim("DTGeometry") << "Parameters empty for DT layer with detid: " << layerId.rawId();
0166       continue;
0167     }
0168
0169     float width = it->surface().bounds().width();
0170     assert(width == parameters[6]);
0171
0172     float thickness = it->surface().bounds().thickness();
0173     assert(thickness == parameters[7]);
0174
0175     float length = it->surface().bounds().length();
0176     assert(length == parameters[8]);
0177
0178     int firstChannel = it->specificTopology().firstChannel();
0179     assert(firstChannel == parameters[3]);
0180
0181     int lastChannel = it->specificTopology().lastChannel();
0182     int nChannels = parameters[5];
0183     assert(nChannels == (lastChannel - firstChannel) + 1);
0184
0185     for (int wireN = firstChannel; wireN - lastChannel <= 0; ++wireN) {
0186       float localX1 = it->specificTopology().wirePosition(wireN);
0187       float localX2 = (wireN - (firstChannel - 1) - 0.5f) * parameters[0] - nChannels / 2.0f * parameters[0];
0188       wire_positions.emplace_back(getDiff(localX1, localX2));
0189     }
0190   }
0191
0192   makeHistogram("DT Layer Wire localX", wire_positions);
0193   makeHistograms("DT Layer");
0194 }
0195
0196 void DTGeometryValidate::compareTransform(const GlobalPoint& gp, const TGeoMatrix* matrix) {
0197   double local[3] = {0.0, 0.0, 0.0};
0198   double global[3];
0199
0200   matrix->LocalToMaster(local, global);
0201
0202   float distance = getDistance(GlobalPoint(global[0], global[1], global[2]), gp);
0203   globalDistances_.push_back(distance);
0204 }
0205
0206 void DTGeometryValidate::compareShape(const GeomDet* det, const float* shape) {
0207   float shapeTopWidth;
0208   float shapeBottomWidth;
0209   float shapeLength;
0210   float shapeThickness;
0211
0212   if (shape[0] == 1) {
0213     shapeTopWidth = shape[2];
0214     shapeBottomWidth = shape[1];
0215     shapeLength = shape[4];
0216     shapeThickness = shape[3];
0217   } else if (shape[0] == 2) {
0218     shapeTopWidth = shape[1];
0219     shapeBottomWidth = shape[1];
0220     shapeLength = shape[2];
0221     shapeThickness = shape[3];
0222   } else {
0223     LogVerbatim("DTGeometry") << "Failed to get box or trapezoid from shape";
0224     return;
0225   }
0226
0227   float topWidth, bottomWidth;
0228   float length, thickness;
0229
0230   const Bounds* bounds = &(det->surface().bounds());
0231
0232   if (const TrapezoidalPlaneBounds* tpbs = dynamic_cast<const TrapezoidalPlaneBounds*>(bounds)) {
0233     array<const float, 4> const& ps = tpbs->parameters();
0234
0235     assert(ps.size() == 4);
0236
0237     bottomWidth = ps[0];
0238     topWidth = ps[1];
0239     thickness = ps[2];
0240     length = ps[3];
0241   } else if ((dynamic_cast<const RectangularPlaneBounds*>(bounds))) {
0242     length = det->surface().bounds().length() * 0.5;
0243     topWidth = det->surface().bounds().width() * 0.5;
0244     bottomWidth = topWidth;
0245     thickness = det->surface().bounds().thickness() * 0.5;
0246   } else {
0247     LogVerbatim("DTGeometry") << "Failed to get bounds";
0248     return;
0249   }
0250   topWidths_.push_back(fabs(shapeTopWidth - topWidth));
0251   bottomWidths_.push_back(fabs(shapeBottomWidth - bottomWidth));
0252   lengths_.push_back(fabs(shapeLength - length));
0253   thicknesses_.push_back(fabs(shapeThickness - thickness));
0254 }
0255
0256 float DTGeometryValidate::getDistance(const GlobalPoint& p1, const GlobalPoint& p2) {
0257   return sqrt((p1.x() - p2.x()) * (p1.x() - p2.x()) + (p1.y() - p2.y()) * (p1.y() - p2.y()) +
0258               (p1.z() - p2.z()) * (p1.z() - p2.z()));
0259 }
0260
0261 float DTGeometryValidate::getDiff(const float val1, const float val2) {
0262   if (almost_equal(val1, val2, tolerance_))
0263     return 0.0f;
0264   else
0265     return (val1 - val2);
0266 }
0267
0268 void DTGeometryValidate::makeHistograms(const char* detector) {
0269   outFile_->cd();
0270
0271   string d(detector);
0272
0273   string gdn = d + ": distance between points in global coordinates";
0274   makeHistogram(gdn, globalDistances_);
0275
0276   string twn = d + ": absolute difference between top widths (along X)";
0277   makeHistogram(twn, topWidths_);
0278
0279   string bwn = d + ": absolute difference between bottom widths (along X)";
0280   makeHistogram(bwn, bottomWidths_);
0281
0282   string ln = d + ": absolute difference between lengths (along Y)";
0283   makeHistogram(ln, lengths_);
0284
0285   string tn = d + ": absolute difference between thicknesses (along Z)";
0286   makeHistogram(tn, thicknesses_);
0287 }
0288
0289 void DTGeometryValidate::makeHistogram(const string& name, vector<float>& data) {
0290   if (data.empty())
0291     return;
0292
0293   const auto [minE, maxE] = minmax_element(begin(data), end(data));
0294
0295   TH1D hist(name.c_str(), name.c_str(), 100, *minE * (1 + 0.10), *maxE * (1 + 0.10));
0296
0297   for (auto const& it : data)
0298     hist.Fill(it);
0299
0300   hist.GetXaxis()->SetTitle("[cm]");
0301   hist.Write();
0302 }
0303
0304 void DTGeometryValidate::beginJob() { outFile_->cd(); }
0305
0306 void DTGeometryValidate::endJob() {
0307   LogVerbatim("DTGeometry") << "Done.";
0308   LogVerbatim("DTGeometry") << "Results written to " << outfileName_;
0309   outFile_->Close();
0310 }
0311
0312 DEFINE_FWK_MODULE(DTGeometryValidate);