GEMGeometryBuilder/plugins/ME0GeometryValidate.cc

0001 /*
0002 //\class ME0GeometryValidate
0003
0004  Description: ME0 GeometryValidate for DD4hep
0005
0006  //
0007 // Author:  Sergio Lo Meo (sergio.lo.meo@cern.ch) following what Ianna Osborne made for DTs (DD4hep migration)
0008 //          Created:  29 Apr 2020
0009 */
0010
0011 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016
0017 #include "Geometry/GEMGeometry/interface/ME0Geometry.h"
0018 #include "Geometry/GEMGeometry/interface/ME0Layer.h"
0019 #include "Geometry/GEMGeometry/interface/ME0Chamber.h"
0020 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0021 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0022
0023 #include "Fireworks/Core/interface/FWGeometry.h"
0024
0025 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0026 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0027
0028 #include <TFile.h>
0029 #include <TH1.h>
0030
0031 #include <limits>
0032 #include <string>
0033 #include <type_traits>
0034 #include <algorithm>
0035
0036 using namespace std;
0037
0038 template <class T>
0039 typename enable_if<!numeric_limits<T>::is_integer, bool>::type almost_equal(T x, T y, int ulp) {
0040   // the machine epsilon has to be scaled to the magnitude of the values used
0041   // and multiplied by the desired precision in ULPs (units in the last place)
0042   return abs(x - y) <= numeric_limits<T>::epsilon() * abs(x + y) * ulp
0043          // unless the result is subnormal
0044          || abs(x - y) < numeric_limits<T>::min();
0045 }
0046
0047 using namespace edm;
0048
0049 class ME0GeometryValidate : public one::EDAnalyzer<> {
0050 public:
0051   explicit ME0GeometryValidate(const ParameterSet&);
0052   ~ME0GeometryValidate() override {}
0053
0054 private:
0055   void beginJob() override;
0056   void analyze(const edm::Event&, const edm::EventSetup&) override;
0057   void endJob() override;
0058
0059   void validateME0ChamberGeometry();
0060   void validateME0EtaPartitionGeometry();
0061   void validateME0EtaPartitionGeometry2();
0062
0063   void compareTransform(const GlobalPoint&, const TGeoMatrix*);
0064   void compareShape(const GeomDet*, const float*);
0065
0066   float getDistance(const GlobalPoint&, const GlobalPoint&);
0067   float getDiff(const float, const float);
0068
0069   void makeHistograms(const char*);
0070   void makeHistograms2(const char*);
0071   void makeHistogram(const string&, vector<float>&);
0072
0073   void clearData() {
0074     globalDistances_.clear();
0075     topWidths_.clear();
0076     bottomWidths_.clear();
0077     lengths_.clear();
0078     thicknesses_.clear();
0079   }
0080
0081   void clearData2() {
0082     nstrips_.clear();
0083     pitch_.clear();
0084     stripslen_.clear();
0085   }
0086
0087   const edm::ESGetToken<ME0Geometry, MuonGeometryRecord> tokGeom_;
0088   const ME0Geometry* me0Geometry_;
0089   FWGeometry fwGeometry_;
0090   TFile* outFile_;
0091   vector<float> globalDistances_;
0092   vector<float> topWidths_;
0093   vector<float> bottomWidths_;
0094   vector<float> lengths_;
0095   vector<float> thicknesses_;
0096   vector<float> nstrips_;
0097   vector<float> pitch_;
0098   vector<float> stripslen_;
0099   string infileName_;
0100   string outfileName_;
0101   int tolerance_;
0102 };
0103
0104 ME0GeometryValidate::ME0GeometryValidate(const edm::ParameterSet& iConfig)
0105     : tokGeom_{esConsumes<ME0Geometry, MuonGeometryRecord>(edm::ESInputTag{})},
0106       infileName_(iConfig.getUntrackedParameter<string>("infileName", "cmsRecoGeom-2026.root")),
0107       outfileName_(iConfig.getUntrackedParameter<string>("outfileName", "validateME0Geometry.root")),
0108       tolerance_(iConfig.getUntrackedParameter<int>("tolerance", 6)) {
0109   fwGeometry_.loadMap(infileName_.c_str());
0110   outFile_ = TFile::Open(outfileName_.c_str(), "RECREATE");
0111 }
0112
0113 void ME0GeometryValidate::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
0114   me0Geometry_ = &eventSetup.getData(tokGeom_);
0115   LogTrace("ME0Geometry") << "Validating ME0 chamber geometry";
0116   validateME0ChamberGeometry();
0117   validateME0EtaPartitionGeometry2();
0118   validateME0EtaPartitionGeometry();
0119 }
0120
0121 void ME0GeometryValidate::validateME0ChamberGeometry() {
0122   clearData();
0123
0124   for (auto const& it : me0Geometry_->chambers()) {
0125     ME0DetId chId = it->id();
0126     GlobalPoint gp = it->surface().toGlobal(LocalPoint(0.0, 0.0, 0.0));
0127     const TGeoMatrix* matrix = fwGeometry_.getMatrix(chId.rawId());
0128
0129     if (!matrix) {
0130       LogTrace("ME0Geometry") << "Failed to get matrix of ME0 chamber with detid: " << chId.rawId();
0131       continue;
0132     }
0133     compareTransform(gp, matrix);
0134
0135     auto const& shape = fwGeometry_.getShapePars(chId.rawId());
0136
0137     if (!shape) {
0138       LogTrace("ME0Geometry") << "Failed to get shape of ME0 chamber with detid: " << chId.rawId();
0139       continue;
0140     }
0141     compareShape(it, shape);
0142   }
0143   makeHistograms("ME0 Chamber");
0144 }
0145
0146 void ME0GeometryValidate::validateME0EtaPartitionGeometry2() {
0147   clearData();
0148
0149   for (auto const& it : me0Geometry_->etaPartitions()) {
0150     ME0DetId chId = it->id();
0151     GlobalPoint gp = it->surface().toGlobal(LocalPoint(0.0, 0.0, 0.0));
0152     const TGeoMatrix* matrix = fwGeometry_.getMatrix(chId.rawId());
0153
0154     if (!matrix) {
0155       LogTrace("ME0Geometry") << "Failed to get matrix of ME0 eta partition 2 with detid: " << chId.rawId();
0156       continue;
0157     }
0158     compareTransform(gp, matrix);
0159
0160     auto const& shape = fwGeometry_.getShapePars(chId.rawId());
0161
0162     if (!shape) {
0163       LogTrace("ME0Geometry") << "Failed to get shape of ME0 eta partition 2 with detid: " << chId.rawId();
0164       continue;
0165     }
0166     compareShape(it, shape);
0167   }
0168   makeHistograms("ME0 Eta Partition");
0169 }
0170
0171 void ME0GeometryValidate::validateME0EtaPartitionGeometry() {
0172   clearData2();
0173
0174   for (auto const& it : me0Geometry_->etaPartitions()) {
0175     ME0DetId chId = it->id();
0176     const int n_strips = it->nstrips();
0177     const float n_pitch = it->pitch();
0178     const StripTopology& topo = it->specificTopology();
0179     const float stripLen = topo.stripLength();
0180     const float* parameters = fwGeometry_.getParameters(chId.rawId());
0181     nstrips_.push_back(std::abs(n_strips - parameters[0]));
0182
0183     if (n_strips) {
0184       for (int istrips = 1; istrips <= n_strips; istrips++) {
0185         if (std::abs(n_pitch - parameters[2]) < 1.0e-05)
0186           pitch_.push_back(0.0f);
0187         else
0188           pitch_.push_back(std::abs(n_pitch - parameters[2]));
0189         if (std::abs(stripLen - parameters[1]) < 1.0e-05)
0190           pitch_.push_back(0.0f);
0191         else
0192           stripslen_.push_back(std::abs(stripLen - parameters[1]));
0193       }
0194     } else {
0195       LogWarning("ME0Geometry") << "ATTENTION! nStrips == 0";
0196     }
0197   }
0198   makeHistograms2("ME0 Eta Partition");
0199 }
0200
0201 void ME0GeometryValidate::compareTransform(const GlobalPoint& gp, const TGeoMatrix* matrix) {
0202   double local[3] = {0.0, 0.0, 0.0};
0203   double global[3];
0204
0205   matrix->LocalToMaster(local, global);
0206
0207   float distance = getDistance(GlobalPoint(global[0], global[1], global[2]), gp);
0208   if (abs(distance) < 1.0e-7)
0209     distance = 0.0;  // set a tollerance for the distance inside Histos
0210   globalDistances_.push_back(distance);
0211 }
0212
0213 void ME0GeometryValidate::compareShape(const GeomDet* det, const float* shape) {
0214   float shapeTopWidth;
0215   float shapeBottomWidth;
0216   float shapeLength;
0217   float shapeThickness;
0218
0219   if (shape[0] == 1) {
0220     shapeTopWidth = shape[2];
0221     shapeBottomWidth = shape[1];
0222     shapeLength = shape[4];
0223     shapeThickness = shape[3];
0224   } else if (shape[0] == 2) {
0225     shapeTopWidth = shape[1];
0226     shapeBottomWidth = shape[1];
0227     shapeLength = shape[2];
0228     shapeThickness = shape[3];
0229   } else {
0230     LogTrace("ME0Geometry") << "Failed to get box or trapezoid from shape";
0231
0232     return;
0233   }
0234
0235   float topWidth, bottomWidth;
0236   float length, thickness;
0237
0238   const Bounds* bounds = &(det->surface().bounds());
0239   if (const TrapezoidalPlaneBounds* tpbs = dynamic_cast<const TrapezoidalPlaneBounds*>(bounds)) {
0240     array<const float, 4> const& ps = tpbs->parameters();
0241
0242     assert(ps.size() == 4);
0243
0244     bottomWidth = ps[0];
0245     topWidth = ps[1];
0246     thickness = ps[2];
0247     length = ps[3];
0248
0249   } else if ((dynamic_cast<const RectangularPlaneBounds*>(bounds))) {
0250     length = det->surface().bounds().length() * 0.5;
0251     topWidth = det->surface().bounds().width() * 0.5;
0252     bottomWidth = topWidth;
0253     thickness = det->surface().bounds().thickness() * 0.5;
0254
0255   } else {
0256     LogTrace("ME0Geometry") << "Failed to get bounds";
0257
0258     return;
0259   }
0260   topWidths_.push_back(std::abs(shapeTopWidth - topWidth));
0261   bottomWidths_.push_back(std::abs(shapeBottomWidth - bottomWidth));
0262   lengths_.push_back(std::abs(shapeLength - length));
0263   thicknesses_.push_back(std::abs(shapeThickness - thickness));
0264 }
0265
0266 float ME0GeometryValidate::getDistance(const GlobalPoint& p1, const GlobalPoint& p2) {
0267   return sqrt((p1.x() - p2.x()) * (p1.x() - p2.x()) + (p1.y() - p2.y()) * (p1.y() - p2.y()) +
0268               (p1.z() - p2.z()) * (p1.z() - p2.z()));
0269 }
0270
0271 float ME0GeometryValidate::getDiff(const float val1, const float val2) {
0272   if (almost_equal(val1, val2, tolerance_))
0273     return 0.0f;
0274   else
0275     return (val1 - val2);
0276 }
0277
0278 void ME0GeometryValidate::makeHistograms(const char* detector) {
0279   outFile_->cd();
0280
0281   string d(detector);
0282
0283   string gdn = d + ": distance between points in global coordinates";
0284   makeHistogram(gdn, globalDistances_);
0285
0286   string twn = d + ": absolute difference between top widths (along X)";
0287   makeHistogram(twn, topWidths_);
0288
0289   string bwn = d + ": absolute difference between bottom widths (along X)";
0290   makeHistogram(bwn, bottomWidths_);
0291
0292   string ln = d + ": absolute difference between lengths (along Y)";
0293   makeHistogram(ln, lengths_);
0294
0295   string tn = d + ": absolute difference between thicknesses (along Z)";
0296   makeHistogram(tn, thicknesses_);
0297 }
0298
0299 void ME0GeometryValidate::makeHistograms2(const char* detector) {
0300   outFile_->cd();
0301
0302   string d(detector);
0303
0304   string ns = d + ": absolute difference between nStrips";
0305   makeHistogram(ns, nstrips_);
0306
0307   string pi = d + ": absolute difference between Strips Pitch";
0308   makeHistogram(pi, pitch_);
0309
0310   string pl = d + ": absolute difference between Strips Length";
0311   makeHistogram(pl, stripslen_);
0312 }
0313
0314 void ME0GeometryValidate::makeHistogram(const string& name, vector<float>& data) {
0315   if (data.empty())
0316     return;
0317
0318   const auto [minE, maxE] = minmax_element(begin(data), end(data));
0319
0320   TH1D hist(name.c_str(), name.c_str(), 100, *minE * (1 + 0.10), *maxE * (1 + 0.10));
0321
0322   for (auto const& it : data)
0323     hist.Fill(it);
0324
0325   hist.GetXaxis()->SetTitle("[cm]");
0326   hist.Write();
0327 }
0328
0329 void ME0GeometryValidate::beginJob() { outFile_->cd(); }
0330
0331 void ME0GeometryValidate::endJob() {
0332   LogTrace("ME0Geometry") << "Done.";
0333   LogTrace("ME0Geometry") << "Results written to " << outfileName_;
0334   outFile_->Close();
0335 }
0336
0337 DEFINE_FWK_MODULE(ME0GeometryValidate);