GEMGeometryBuilder/plugins/GEMGeometryValidate.cc

0001 /*
0002 //\class GEMGeometryValidate
0003
0004  Description: GEM GeometryValidate from DD & DD4hep
0005
0006  //
0007 // Author:  Sergio Lo Meo (sergio.lo.meo@cern.ch) following what Ianna Osburne made for DTs (DD4hep migration)
0008 //          Created:  27 Jan 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/GEMGeometry.h"
0018 #include "Geometry/GEMGeometry/interface/GEMChamber.h"
0019 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0020 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0021
0022 #include "Fireworks/Core/interface/FWGeometry.h"
0023
0024 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0025 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0026
0027 #include <TFile.h>
0028 #include <TH1.h>
0029
0030 #include <limits>
0031 #include <string>
0032 #include <type_traits>
0033 #include <algorithm>
0034 #include <cmath>
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 GEMGeometryValidate : public one::EDAnalyzer<> {
0050 public:
0051   explicit GEMGeometryValidate(const ParameterSet&);
0052   ~GEMGeometryValidate() 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 validateGEMChamberGeometry();
0060   void validateGEMEtaPartitionGeometry();
0061
0062   void compareTransform(const GlobalPoint&, const TGeoMatrix*);
0063   void compareShape(const GeomDet*, const float*);
0064
0065   float getDistance(const GlobalPoint&, const GlobalPoint&);
0066   float getDiff(const float, const float);
0067
0068   void makeHistograms(const char*);
0069   void makeHistograms2(const char*);
0070   void makeHistogram(const string&, vector<float>&);
0071
0072   void clearData() {
0073     globalDistances_.clear();
0074     topWidths_.clear();
0075     bottomWidths_.clear();
0076     lengths_.clear();
0077     thicknesses_.clear();
0078   }
0079
0080   void clearData2() {
0081     nstrips_.clear();
0082     pitch_.clear();
0083     stripslen_.clear();
0084   }
0085
0086   const edm::ESGetToken<GEMGeometry, MuonGeometryRecord> tokGeom_;
0087   const GEMGeometry* gemGeometry_;
0088   FWGeometry fwGeometry_;
0089   TFile* outFile_;
0090   vector<float> globalDistances_;
0091   vector<float> topWidths_;
0092   vector<float> bottomWidths_;
0093   vector<float> lengths_;
0094   vector<float> thicknesses_;
0095   vector<float> nstrips_;
0096   vector<float> pitch_;
0097   vector<float> stripslen_;
0098   string infileName_;
0099   string outfileName_;
0100   int tolerance_;
0101 };
0102
0103 GEMGeometryValidate::GEMGeometryValidate(const edm::ParameterSet& iConfig)
0104     : tokGeom_{esConsumes<GEMGeometry, MuonGeometryRecord>(edm::ESInputTag{})},
0105       infileName_(iConfig.getUntrackedParameter<string>("infileName", "cmsRecoGeom-2021.root")),
0106       outfileName_(iConfig.getUntrackedParameter<string>("outfileName", "validateGEMGeometry.root")),
0107       tolerance_(iConfig.getUntrackedParameter<int>("tolerance", 6)) {
0108   fwGeometry_.loadMap(infileName_.c_str());
0109   outFile_ = TFile::Open(outfileName_.c_str(), "RECREATE");
0110 }
0111
0112 void GEMGeometryValidate::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
0113   gemGeometry_ = &eventSetup.getData(tokGeom_);
0114
0115   LogVerbatim("GEMGeometry") << "Validating GEM chamber geometry";
0116
0117   validateGEMChamberGeometry();
0118
0119   validateGEMEtaPartitionGeometry();
0120 }
0121
0122 void GEMGeometryValidate::validateGEMChamberGeometry() {
0123   clearData();
0124
0125   for (auto const& it : gemGeometry_->chambers()) {
0126     GEMDetId chId = it->id();
0127     GlobalPoint gp = it->surface().toGlobal(LocalPoint(0.0, 0.0, 0.0));
0128
0129     const TGeoMatrix* matrix = fwGeometry_.getMatrix(chId.rawId());
0130
0131     if (!matrix) {
0132       LogVerbatim("GEMGeometry") << "Failed to get matrix of GEM chamber with detid: " << chId.rawId();
0133       continue;
0134     }
0135     compareTransform(gp, matrix);
0136
0137     auto const& shape = fwGeometry_.getShapePars(chId.rawId());
0138
0139     if (!shape) {
0140       LogVerbatim("GEMGeometry") << "Failed to get shape of GEM chamber with detid: " << chId.rawId();
0141       continue;
0142     }
0143     compareShape(it, shape);
0144   }
0145   makeHistograms("GEM Chamber");
0146 }
0147
0148 void GEMGeometryValidate::validateGEMEtaPartitionGeometry() {
0149   clearData2();
0150
0151   for (auto const& it : gemGeometry_->etaPartitions()) {
0152     GEMDetId chId = it->id();
0153     const int n_strips = it->nstrips();
0154     const float n_pitch = it->pitch();
0155     const StripTopology& topo = it->specificTopology();
0156     const float stripLen = topo.stripLength();
0157     const float* parameters = fwGeometry_.getParameters(chId.rawId());
0158     nstrips_.push_back(abs(n_strips - parameters[0]));
0159
0160     if (n_strips) {
0161       for (int istrips = 1; istrips <= n_strips; istrips++) {
0162         pitch_.push_back(fabs(n_pitch - parameters[2]));
0163
0164         stripslen_.push_back(fabs(stripLen - parameters[1]));
0165       }
0166     } else {
0167       LogVerbatim("GEMGeometry") << "ATTENTION! nStrips == 0";
0168     }
0169   }
0170   makeHistograms2("GEM Eta Partition");
0171 }
0172
0173 void GEMGeometryValidate::compareTransform(const GlobalPoint& gp, const TGeoMatrix* matrix) {
0174   double local[3] = {0.0, 0.0, 0.0};
0175   double global[3];
0176
0177   matrix->LocalToMaster(local, global);
0178
0179   float distance = getDistance(GlobalPoint(global[0], global[1], global[2]), gp);
0180   if (abs(distance) < 1.0e-7)
0181     distance = 0.0;  // set a tollerance for the distance inside Histos
0182   globalDistances_.push_back(distance);
0183 }
0184
0185 void GEMGeometryValidate::compareShape(const GeomDet* det, const float* shape) {
0186   float shapeTopWidth;
0187   float shapeBottomWidth;
0188   float shapeLength;
0189   float shapeThickness;
0190
0191   if (shape[0] == 1) {
0192     shapeTopWidth = shape[2];
0193     shapeBottomWidth = shape[1];
0194     shapeLength = shape[4];
0195     shapeThickness = shape[3];
0196   } else if (shape[0] == 2) {
0197     shapeTopWidth = shape[1];
0198     shapeBottomWidth = shape[1];
0199     shapeLength = shape[2];
0200     shapeThickness = shape[3];
0201   } else {
0202     LogVerbatim("GEMGeometry") << "Failed to get box or trapezoid from shape";
0203
0204     return;
0205   }
0206
0207   float topWidth, bottomWidth;
0208   float length, thickness;
0209
0210   const Bounds* bounds = &(det->surface().bounds());
0211   if (const TrapezoidalPlaneBounds* tpbs = dynamic_cast<const TrapezoidalPlaneBounds*>(bounds)) {
0212     array<const float, 4> const& ps = tpbs->parameters();
0213
0214     assert(ps.size() == 4);
0215
0216     bottomWidth = ps[0];
0217     topWidth = ps[1];
0218     thickness = ps[2];
0219     length = ps[3];
0220   } else if ((dynamic_cast<const RectangularPlaneBounds*>(bounds))) {
0221     length = det->surface().bounds().length() * 0.5;
0222     topWidth = det->surface().bounds().width() * 0.5;
0223     bottomWidth = topWidth;
0224     thickness = det->surface().bounds().thickness() * 0.5;
0225   } else {
0226     LogVerbatim("GEMGeometry") << "Failed to get bounds";
0227
0228     return;
0229   }
0230   topWidths_.push_back(fabs(shapeTopWidth - topWidth));
0231   bottomWidths_.push_back(fabs(shapeBottomWidth - bottomWidth));
0232   lengths_.push_back(fabs(shapeLength - length));
0233   thicknesses_.push_back(fabs(shapeThickness - thickness));
0234 }
0235
0236 float GEMGeometryValidate::getDistance(const GlobalPoint& p1, const GlobalPoint& p2) {
0237   return sqrt((p1.x() - p2.x()) * (p1.x() - p2.x()) + (p1.y() - p2.y()) * (p1.y() - p2.y()) +
0238               (p1.z() - p2.z()) * (p1.z() - p2.z()));
0239 }
0240
0241 float GEMGeometryValidate::getDiff(const float val1, const float val2) {
0242   if (almost_equal(val1, val2, tolerance_))
0243     return 0.0f;
0244   else
0245     return (val1 - val2);
0246 }
0247
0248 void GEMGeometryValidate::makeHistograms(const char* detector) {
0249   outFile_->cd();
0250
0251   string d(detector);
0252
0253   string gdn = d + ": distance between points in global coordinates";
0254   makeHistogram(gdn, globalDistances_);
0255
0256   string twn = d + ": absolute difference between top widths (along X)";
0257   makeHistogram(twn, topWidths_);
0258
0259   string bwn = d + ": absolute difference between bottom widths (along X)";
0260   makeHistogram(bwn, bottomWidths_);
0261
0262   string ln = d + ": absolute difference between lengths (along Y)";
0263   makeHistogram(ln, lengths_);
0264
0265   string tn = d + ": absolute difference between thicknesses (along Z)";
0266   makeHistogram(tn, thicknesses_);
0267 }
0268
0269 void GEMGeometryValidate::makeHistograms2(const char* detector) {
0270   outFile_->cd();
0271
0272   string d(detector);
0273
0274   string ns = d + ": absolute difference between nStrips";
0275   makeHistogram(ns, nstrips_);
0276
0277   string pi = d + ": absolute difference between Strips Pitch";
0278   makeHistogram(pi, pitch_);
0279
0280   string pl = d + ": absolute difference between Strips Length";
0281   makeHistogram(pl, stripslen_);
0282 }
0283
0284 void GEMGeometryValidate::makeHistogram(const string& name, vector<float>& data) {
0285   if (data.empty())
0286     return;
0287
0288   const auto [minE, maxE] = minmax_element(begin(data), end(data));
0289
0290   TH1D hist(name.c_str(), name.c_str(), 100, *minE * (1 + 0.10), *maxE * (1 + 0.10));
0291
0292   for (auto const& it : data)
0293     hist.Fill(it);
0294
0295   hist.GetXaxis()->SetTitle("[cm]");
0296   hist.Write();
0297 }
0298
0299 void GEMGeometryValidate::beginJob() { outFile_->cd(); }
0300
0301 void GEMGeometryValidate::endJob() {
0302   LogVerbatim("GEMGeometry") << "Done.";
0303   LogVerbatim("GEMGeometry") << "Results written to " << outfileName_;
0304   outFile_->Close();
0305 }
0306
0307 DEFINE_FWK_MODULE(GEMGeometryValidate);