Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 13:02:42

0001 /*
0002 //\class CSCGeometryValidate
0003 
0004  Description: CSC 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:  Thu, 05 March 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/CSCGeometry/interface/CSCGeometry.h"
0018 #include "Geometry/CSCGeometry/interface/CSCChamber.h"
0019 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0020 #include "Geometry/CSCGeometry/interface/CSCLayerGeometry.h"
0021 #include "Geometry/CSCGeometry/interface/CSCStripTopology.h"
0022 #include "Geometry/CSCGeometry/interface/CSCWireTopology.h"
0023 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0024 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0025 
0026 #include "Fireworks/Core/interface/FWGeometry.h"
0027 
0028 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0029 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0030 
0031 #include <TFile.h>
0032 #include <TH1.h>
0033 
0034 #include <limits>
0035 #include <string>
0036 #include <type_traits>
0037 #include <algorithm>
0038 
0039 using namespace std;
0040 
0041 template <class T>
0042 typename enable_if<!numeric_limits<T>::is_integer, bool>::type almost_equal(T x, T y, int ulp) {
0043   // the machine epsilon has to be scaled to the magnitude of the values used
0044   // and multiplied by the desired precision in ULPs (units in the last place)
0045   return abs(x - y) <= numeric_limits<T>::epsilon() * abs(x + y) * ulp
0046          // unless the result is subnormal
0047          || abs(x - y) < numeric_limits<T>::min();
0048 }
0049 
0050 using namespace edm;
0051 
0052 class CSCGeometryValidate : public one::EDAnalyzer<> {
0053 public:
0054   explicit CSCGeometryValidate(const ParameterSet&);
0055   ~CSCGeometryValidate() override {}
0056 
0057 private:
0058   void beginJob() override;
0059   void analyze(const edm::Event&, const edm::EventSetup&) override;
0060   void endJob() override;
0061 
0062   void validateCSCChamberGeometry();
0063   void validateCSCLayerGeometry();
0064 
0065   void compareTransform(const GlobalPoint&, const TGeoMatrix*);
0066   void compareShape(const GeomDet*, const float*);
0067 
0068   float getDistance(const GlobalPoint&, const GlobalPoint&);
0069   float getDiff(const float, const float);
0070 
0071   void makeHistograms(const char*);
0072   void makeHistograms2(const char*);
0073   void makeHistogram(const string&, vector<float>&);
0074 
0075   void clearData() {
0076     globalDistances_.clear();
0077     topWidths_.clear();
0078     bottomWidths_.clear();
0079     lengths_.clear();
0080     thicknesses_.clear();
0081   }
0082 
0083   void clearData2() {
0084     yAxisOrientation_.clear();
0085     sOffset_.clear();
0086     yCentreOfStripPlane_.clear();
0087     angularWidth_.clear();
0088     centreToIntersection_.clear();
0089     phiOfOneEdge_.clear();
0090     wireSpacing_.clear();
0091     wireAngle_.clear();
0092   }
0093 
0094   const edm::ESGetToken<CSCGeometry, MuonGeometryRecord> tokCSC_;
0095   const CSCGeometry* cscGeometry_;
0096   FWGeometry fwGeometry_;
0097   TFile* outFile_;
0098   //chambers
0099   vector<float> globalDistances_;
0100   vector<float> topWidths_;
0101   vector<float> bottomWidths_;
0102   vector<float> lengths_;
0103   vector<float> thicknesses_;
0104   // strips
0105   vector<float> yAxisOrientation_;
0106   vector<float> sOffset_;
0107   vector<float> yCentreOfStripPlane_;
0108   vector<float> angularWidth_;
0109   vector<float> centreToIntersection_;
0110   vector<float> phiOfOneEdge_;
0111   //wires
0112   vector<float> wireSpacing_;
0113   vector<float> wireAngle_;
0114   //files
0115   string infileName_;
0116   string outfileName_;
0117   int tolerance_;
0118 };
0119 
0120 CSCGeometryValidate::CSCGeometryValidate(const edm::ParameterSet& iConfig)
0121     : tokCSC_{esConsumes<CSCGeometry, MuonGeometryRecord>(edm::ESInputTag{})},
0122       infileName_(iConfig.getUntrackedParameter<string>("infileName", "cmsGeom10.root")),
0123       outfileName_(iConfig.getUntrackedParameter<string>("outfileName", "validateCSCGeometry.root")),
0124       tolerance_(iConfig.getUntrackedParameter<int>("tolerance", 6)) {
0125   fwGeometry_.loadMap(infileName_.c_str());
0126   outFile_ = TFile::Open(outfileName_.c_str(), "RECREATE");
0127 }
0128 
0129 void CSCGeometryValidate::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
0130   cscGeometry_ = &eventSetup.getData(tokCSC_);
0131   LogVerbatim("CSCGeometry") << "Validating CSC chamber geometry";
0132   validateCSCChamberGeometry();
0133   validateCSCLayerGeometry();
0134 }
0135 
0136 void CSCGeometryValidate::validateCSCChamberGeometry() {
0137   clearData();
0138 
0139   for (auto const& it : cscGeometry_->chambers()) {
0140     CSCDetId chId = it->id();
0141     GlobalPoint gp = it->surface().toGlobal(LocalPoint(0.0, 0.0, 0.0));
0142 
0143     const TGeoMatrix* matrix = fwGeometry_.getMatrix(chId.rawId());
0144 
0145     if (!matrix) {
0146       LogVerbatim("CSCGeometry") << "Failed to get matrix of CSC chamber with detid: " << chId.rawId();
0147       continue;
0148     }
0149     compareTransform(gp, matrix);
0150 
0151     auto const& shape = fwGeometry_.getShapePars(chId.rawId());
0152 
0153     if (!shape) {
0154       LogVerbatim("CSCGeometry") << "Failed to get shape of CSC chamber with detid: " << chId.rawId();
0155       continue;
0156     }
0157     compareShape(it, shape);
0158   }
0159   makeHistograms("CSC Chamber");
0160 }
0161 
0162 void CSCGeometryValidate::validateCSCLayerGeometry() {
0163   clearData2();
0164 
0165   for (auto const& it : cscGeometry_->layers()) {
0166     CSCDetId chId = it->id();
0167     const CSCLayerGeometry* laygeo = it->geometry();
0168     const int n_strips = laygeo->numberOfStrips();
0169     const int n_wire = laygeo->numberOfWires();
0170     const float strips_offset = laygeo->stripOffset();
0171     const CSCStripTopology* stopo = laygeo->topology();
0172     const float ycentre_of_strip_plane = stopo->yCentreOfStripPlane();
0173     const float angular_width = stopo->angularWidth();
0174     const float y_axis_orientation = stopo->yAxisOrientation();
0175     const float centre_to_intersection = stopo->centreToIntersection();
0176     const float phi_of_one_edge = stopo->phiOfOneEdge();
0177     const float* parameters = fwGeometry_.getParameters(chId.rawId());
0178     const CSCWireTopology* wiretopo = laygeo->wireTopology();
0179     const double wire_spacing = wiretopo->wireSpacing();
0180     const float wire_angle = wiretopo->wireAngle();
0181 
0182     if (n_strips) {
0183       for (int istrips = 1; istrips <= n_strips; istrips++) {
0184         yAxisOrientation_.push_back(fabs(y_axis_orientation - parameters[0]));
0185         sOffset_.push_back(fabs(strips_offset - parameters[4]));
0186         yCentreOfStripPlane_.push_back(fabs(ycentre_of_strip_plane - parameters[2]));
0187         angularWidth_.push_back(fabs(angular_width - parameters[5]));
0188         centreToIntersection_.push_back(fabs(centre_to_intersection - parameters[1]));
0189         phiOfOneEdge_.push_back(fabs(phi_of_one_edge - parameters[3]));
0190       }
0191     } else {
0192       LogVerbatim("CSCGeometry") << "ATTENTION! nStrips == 0";
0193     }
0194 
0195     if (n_wire) {
0196       for (int iwires = 1; iwires <= n_wire; iwires++) {
0197         wireSpacing_.push_back(fabs(wire_spacing - parameters[6]));
0198         wireAngle_.push_back(fabs(wire_angle - parameters[7]));
0199       }
0200     } else {
0201       LogVerbatim("CSCGeometry") << "ATTENTION! nWires == 0";
0202     }
0203   }
0204   makeHistograms2("CSC Layer");
0205 }
0206 
0207 void CSCGeometryValidate::compareTransform(const GlobalPoint& gp, const TGeoMatrix* matrix) {
0208   double local[3] = {0.0, 0.0, 0.0};
0209   double global[3];
0210 
0211   matrix->LocalToMaster(local, global);
0212 
0213   float distance = getDistance(GlobalPoint(global[0], global[1], global[2]), gp);
0214   if ((distance >= 0.0) && (distance < 1.0e-7))
0215     distance = 0.0;  // set a tollerance for the distance inside Histos
0216   globalDistances_.push_back(distance);
0217 }
0218 
0219 void CSCGeometryValidate::compareShape(const GeomDet* det, const float* shape) {
0220   float shapeTopWidth;
0221   float shapeBottomWidth;
0222   float shapeLength;
0223   float shapeThickness;
0224 
0225   if (shape[0] == 1) {
0226     shapeTopWidth = shape[2];
0227     shapeBottomWidth = shape[1];
0228     shapeLength = shape[4];
0229     shapeThickness = shape[3];
0230   } else if (shape[0] == 2) {
0231     shapeTopWidth = shape[1];
0232     shapeBottomWidth = shape[1];
0233     shapeLength = shape[2];
0234     shapeThickness = shape[3];
0235   } else {
0236     LogVerbatim("CSCGeometry") << "Failed to get box or trapezoid from shape";
0237 
0238     return;
0239   }
0240 
0241   float topWidth, bottomWidth;
0242   float length, thickness;
0243 
0244   const Bounds* bounds = &(det->surface().bounds());
0245   if (const TrapezoidalPlaneBounds* tpbs = dynamic_cast<const TrapezoidalPlaneBounds*>(bounds)) {
0246     array<const float, 4> const& ps = tpbs->parameters();
0247 
0248     assert(ps.size() == 4);
0249 
0250     bottomWidth = ps[0];
0251     topWidth = ps[1];
0252     thickness = ps[2];
0253     length = ps[3];
0254   } else if ((dynamic_cast<const RectangularPlaneBounds*>(bounds))) {
0255     length = det->surface().bounds().length() * 0.5;
0256     topWidth = det->surface().bounds().width() * 0.5;
0257     bottomWidth = topWidth;
0258     thickness = det->surface().bounds().thickness() * 0.5;
0259   } else {
0260     LogVerbatim("CSCGeometry") << "Failed to get bounds";
0261 
0262     return;
0263   }
0264   topWidths_.push_back(fabs(shapeTopWidth - topWidth));
0265   bottomWidths_.push_back(fabs(shapeBottomWidth - bottomWidth));
0266   lengths_.push_back(fabs(shapeLength - length));
0267   thicknesses_.push_back(fabs(shapeThickness - thickness));
0268 }
0269 
0270 float CSCGeometryValidate::getDistance(const GlobalPoint& p1, const GlobalPoint& p2) {
0271   return sqrt((p1.x() - p2.x()) * (p1.x() - p2.x()) + (p1.y() - p2.y()) * (p1.y() - p2.y()) +
0272               (p1.z() - p2.z()) * (p1.z() - p2.z()));
0273 }
0274 
0275 float CSCGeometryValidate::getDiff(const float val1, const float val2) {
0276   if (almost_equal(val1, val2, tolerance_))
0277     return 0.0f;
0278   else
0279     return (val1 - val2);
0280 }
0281 
0282 void CSCGeometryValidate::makeHistograms(const char* detector) {
0283   outFile_->cd();
0284 
0285   string d(detector);
0286 
0287   string gdn = d + ": distance between points in global coordinates";
0288   makeHistogram(gdn, globalDistances_);
0289 
0290   string twn = d + ": absolute difference between top widths (along X)";
0291   makeHistogram(twn, topWidths_);
0292 
0293   string bwn = d + ": absolute difference between bottom widths (along X)";
0294   makeHistogram(bwn, bottomWidths_);
0295 
0296   string ln = d + ": absolute difference between lengths (along Y)";
0297   makeHistogram(ln, lengths_);
0298 
0299   string tn = d + ": absolute difference between thicknesses (along Z)";
0300   makeHistogram(tn, thicknesses_);
0301 }
0302 
0303 void CSCGeometryValidate::makeHistograms2(const char* detector) {
0304   outFile_->cd();
0305 
0306   string d(detector);
0307 
0308   string ns = d + ": absolute difference between Y Axis Orientation of the Strips";
0309   makeHistogram(ns, yAxisOrientation_);
0310 
0311   string pi = d + ": absolute difference between Strips Offset";
0312   makeHistogram(pi, sOffset_);
0313 
0314   string pl = d + ": absolute difference between 'Y centre' of the Strips Planes";
0315   makeHistogram(pl, yCentreOfStripPlane_);
0316 
0317   string aw = d + ": absolute difference between 'angular width' of the Strips ";
0318   makeHistogram(aw, angularWidth_);
0319 
0320   string ci = d + ": absolute difference between 'centre to intersection' of the Strips ";
0321   makeHistogram(ci, centreToIntersection_);
0322 
0323   string po = d + ": absolute difference between 'phi of one edge' of the Strips ";
0324   makeHistogram(po, phiOfOneEdge_);
0325 
0326   string ws = d + ": absolute difference between 'wire spacing' of the Wires ";
0327   makeHistogram(ws, wireSpacing_);
0328 
0329   string wa = d + ": absolute difference between 'wire angle' of the Wires ";
0330   makeHistogram(wa, wireAngle_);
0331 }
0332 
0333 void CSCGeometryValidate::makeHistogram(const string& name, vector<float>& data) {
0334   if (data.empty())
0335     return;
0336 
0337   const auto [minE, maxE] = minmax_element(begin(data), end(data));
0338 
0339   TH1D hist(name.c_str(), name.c_str(), 100, *minE * (1 + 0.10), *maxE * (1 + 0.10));
0340 
0341   for (auto const& it : data)
0342     hist.Fill(it);
0343 
0344   hist.GetXaxis()->SetTitle("[cm]");
0345   hist.Write();
0346 }
0347 
0348 void CSCGeometryValidate::beginJob() { outFile_->cd(); }
0349 
0350 void CSCGeometryValidate::endJob() {
0351   LogVerbatim("CSCGeometry") << "Done.";
0352   LogVerbatim("CSCGeometry") << "Results written to " << outfileName_;
0353   outFile_->Close();
0354 }
0355 
0356 DEFINE_FWK_MODULE(CSCGeometryValidate);