Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-08 22:26:19

0001 // -*- C++ -*-
0002 //
0003 //
0004 
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 
0015 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0016 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0017 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0018 #include "Geometry/DTGeometry/interface/DTLayer.h"
0019 
0020 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0021 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0022 
0023 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0024 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0025 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0026 
0027 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0028 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0029 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0030 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
0031 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetType.h"
0032 
0033 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0034 
0035 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0036 #include "Geometry/CommonTopologies/interface/RectangularStripTopology.h"
0037 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
0038 
0039 #include "Fireworks/Core/interface/FWGeometry.h"
0040 #include "Fireworks/Core/interface/fwLog.h"
0041 #include "Fireworks/Tracks/interface/TrackUtils.h"
0042 
0043 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0044 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0045 #include "DataFormats/MuonDetId/interface/DTWireId.h"
0046 #include "DataFormats/HcalDetId/interface/HcalCastorDetId.h"
0047 #include "DataFormats/HcalDetId/interface/HcalZDCDetId.h"
0048 
0049 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0050 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0051 
0052 #include <string>
0053 #include <iostream>
0054 
0055 #include <TEveGeoNode.h>
0056 #include <TGeoVolume.h>
0057 #include <TGeoBBox.h>
0058 #include <TGeoArb8.h>
0059 #include <TFile.h>
0060 #include <TH1.h>
0061 
0062 class ValidateGeometry : public edm::one::EDAnalyzer<> {
0063 public:
0064   explicit ValidateGeometry(const edm::ParameterSet&);
0065   ~ValidateGeometry() override;
0066 
0067 private:
0068   void beginJob() override;
0069   void analyze(const edm::Event&, const edm::EventSetup&) override;
0070   void endJob() override;
0071 
0072   void validateRPCGeometry(const int regionNumber, const char* regionName);
0073 
0074   void validateDTChamberGeometry();
0075   void validateDTLayerGeometry();
0076 
0077   void validateCSChamberGeometry(const int endcap, const char* detname);
0078 
0079   void validateCSCLayerGeometry(const int endcap, const char* detname);
0080 
0081   void validateCaloGeometry(DetId::Detector detector, int subdetector, const char* detname);
0082 
0083   void validateTrackerGeometry(const TrackerGeometry::DetContainer& dets, const char* detname);
0084 
0085   void validatePixelTopology(const TrackerGeometry::DetContainer& dets, const char* detname);
0086 
0087   void validateStripTopology(const TrackerGeometry::DetContainer& dets, const char* detname);
0088 
0089   void compareTransform(const GlobalPoint& point, const TGeoMatrix* matrix);
0090 
0091   void compareShape(const GeomDet* det, const float* shape);
0092 
0093   double getDistance(const GlobalPoint& point1, const GlobalPoint& point2);
0094 
0095   void makeHistograms(const char* detector);
0096   void makeHistogram(const std::string& name, std::vector<double>& data);
0097 
0098   std::string infileName_;
0099   std::string outfileName_;
0100 
0101   edm::ESGetToken<RPCGeometry, MuonGeometryRecord> rpcGeometryToken_;
0102   edm::ESGetToken<DTGeometry, MuonGeometryRecord> dtGeometryToken_;
0103   edm::ESGetToken<CSCGeometry, MuonGeometryRecord> cscGeometryToken_;
0104   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0105   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeometryToken_;
0106 
0107   edm::ESHandle<RPCGeometry> rpcGeometry_;
0108   edm::ESHandle<DTGeometry> dtGeometry_;
0109   edm::ESHandle<CSCGeometry> cscGeometry_;
0110   edm::ESHandle<CaloGeometry> caloGeometry_;
0111   edm::ESHandle<TrackerGeometry> trackerGeometry_;
0112 
0113   FWGeometry fwGeometry_;
0114 
0115   TFile* outFile_;
0116 
0117   std::vector<double> globalDistances_;
0118   std::vector<double> topWidths_;
0119   std::vector<double> bottomWidths_;
0120   std::vector<double> lengths_;
0121   std::vector<double> thicknesses_;
0122 
0123   void clearData() {
0124     globalDistances_.clear();
0125     topWidths_.clear();
0126     bottomWidths_.clear();
0127     lengths_.clear();
0128     thicknesses_.clear();
0129   }
0130 
0131   bool dataEmpty() {
0132     return (globalDistances_.empty() && topWidths_.empty() && bottomWidths_.empty() && lengths_.empty() &&
0133             thicknesses_.empty());
0134   }
0135 
0136   bool doTracker_;
0137   bool doMuon_;
0138   bool doCalo_;
0139 };
0140 
0141 ValidateGeometry::ValidateGeometry(const edm::ParameterSet& iConfig)
0142     : infileName_(iConfig.getUntrackedParameter<std::string>("infileName")),
0143       outfileName_(iConfig.getUntrackedParameter<std::string>("outfileName")) {
0144   doTracker_ = iConfig.getUntrackedParameter<bool>("Tracker", true);
0145   doMuon_ = iConfig.getUntrackedParameter<bool>("Muon", true);
0146   doCalo_ = iConfig.getUntrackedParameter<bool>("Calo", true);
0147 
0148   if (doMuon_) {
0149     rpcGeometryToken_ = esConsumes();
0150     dtGeometryToken_ = esConsumes();
0151     cscGeometryToken_ = esConsumes();
0152   }
0153   if (doCalo_) {
0154     caloGeometryToken_ = esConsumes();
0155   }
0156   if (doTracker_) {
0157     trackerGeometryToken_ = esConsumes();
0158   }
0159 
0160   fwGeometry_.loadMap(infileName_.c_str());
0161 
0162   outFile_ = new TFile(outfileName_.c_str(), "RECREATE");
0163 }
0164 
0165 ValidateGeometry::~ValidateGeometry() {}
0166 
0167 void ValidateGeometry::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
0168   if (doMuon_) {
0169     rpcGeometry_ = eventSetup.getHandle(rpcGeometryToken_);
0170 
0171     if (rpcGeometry_.isValid()) {
0172       std::cout << "Validating RPC -z endcap geometry" << std::endl;
0173       validateRPCGeometry(-1, "RPC -z endcap");
0174 
0175       std::cout << "Validating RPC +z endcap geometry" << std::endl;
0176       validateRPCGeometry(+1, "RPC +z endcap");
0177 
0178       std::cout << "Validating RPC barrel geometry" << std::endl;
0179       validateRPCGeometry(0, "RPC barrel");
0180     } else
0181       fwLog(fwlog::kWarning) << "Invalid RPC geometry" << std::endl;
0182 
0183     dtGeometry_ = eventSetup.getHandle(dtGeometryToken_);
0184 
0185     if (dtGeometry_.isValid()) {
0186       std::cout << "Validating DT chamber geometry" << std::endl;
0187       validateDTChamberGeometry();
0188 
0189       std::cout << "Validating DT layer geometry" << std::endl;
0190       validateDTLayerGeometry();
0191     } else
0192       fwLog(fwlog::kWarning) << "Invalid DT geometry" << std::endl;
0193 
0194     cscGeometry_ = eventSetup.getHandle(cscGeometryToken_);
0195 
0196     if (cscGeometry_.isValid()) {
0197       std::cout << "Validating CSC -z geometry" << std::endl;
0198       validateCSChamberGeometry(-1, "CSC chamber -z endcap");
0199 
0200       std::cout << "Validating CSC +z geometry" << std::endl;
0201       validateCSChamberGeometry(+1, "CSC chamber +z endcap");
0202 
0203       std::cout << "Validating CSC layer -z geometry" << std::endl;
0204       validateCSCLayerGeometry(-1, "CSC layer -z endcap");
0205 
0206       std::cout << "Validating CSC layer +z geometry" << std::endl;
0207       validateCSCLayerGeometry(+1, "CSC layer +z endcap");
0208     } else
0209       fwLog(fwlog::kWarning) << "Invalid CSC geometry" << std::endl;
0210   }
0211 
0212   if (doTracker_) {
0213     trackerGeometry_ = eventSetup.getHandle(trackerGeometryToken_);
0214 
0215     if (trackerGeometry_.isValid()) {
0216       std::cout << "Validating TIB geometry and topology" << std::endl;
0217       validateTrackerGeometry(trackerGeometry_->detsTIB(), "TIB");
0218       validateStripTopology(trackerGeometry_->detsTIB(), "TIB");
0219 
0220       std::cout << "Validating TOB geometry and topology" << std::endl;
0221       validateTrackerGeometry(trackerGeometry_->detsTOB(), "TOB");
0222       validateStripTopology(trackerGeometry_->detsTOB(), "TOB");
0223 
0224       std::cout << "Validating TEC geometry and topology" << std::endl;
0225       validateTrackerGeometry(trackerGeometry_->detsTEC(), "TEC");
0226       validateStripTopology(trackerGeometry_->detsTEC(), "TEC");
0227 
0228       std::cout << "Validating TID geometry and topology" << std::endl;
0229       validateTrackerGeometry(trackerGeometry_->detsTID(), "TID");
0230       validateStripTopology(trackerGeometry_->detsTID(), "TID");
0231 
0232       std::cout << "Validating PXB geometry and topology" << std::endl;
0233       validateTrackerGeometry(trackerGeometry_->detsPXB(), "PXB");
0234       validatePixelTopology(trackerGeometry_->detsPXB(), "PXB");
0235 
0236       std::cout << "Validating PXF geometry and topology" << std::endl;
0237       validateTrackerGeometry(trackerGeometry_->detsPXF(), "PXF");
0238       validatePixelTopology(trackerGeometry_->detsPXF(), "PXF");
0239     } else
0240       fwLog(fwlog::kWarning) << "Invalid Tracker geometry" << std::endl;
0241   }
0242 
0243   if (doCalo_) {
0244     caloGeometry_ = eventSetup.getHandle(caloGeometryToken_);
0245 
0246     if (caloGeometry_.isValid()) {
0247       std::cout << "Validating EB geometry" << std::endl;
0248       validateCaloGeometry(DetId::Ecal, EcalBarrel, "EB");
0249 
0250       std::cout << "Validating EE geometry" << std::endl;
0251       validateCaloGeometry(DetId::Ecal, EcalEndcap, "EE");
0252 
0253       std::cout << "Validating ES geometry" << std::endl;
0254       validateCaloGeometry(DetId::Ecal, EcalPreshower, "ES");
0255 
0256       std::cout << "Validating HB geometry" << std::endl;
0257       validateCaloGeometry(DetId::Hcal, HcalBarrel, "HB");
0258 
0259       std::cout << "Validating HE geometry" << std::endl;
0260       validateCaloGeometry(DetId::Hcal, HcalEndcap, "HE");
0261 
0262       std::cout << "Validating HO geometry" << std::endl;
0263       validateCaloGeometry(DetId::Hcal, HcalOuter, "HO");
0264 
0265       std::cout << "Validating HF geometry" << std::endl;
0266       validateCaloGeometry(DetId::Hcal, HcalForward, "HF");
0267 
0268       std::cout << "Validating Castor geometry" << std::endl;
0269       validateCaloGeometry(DetId::Calo, HcalCastorDetId::SubdetectorId, "Castor");
0270 
0271       std::cout << "Validating ZDC geometry" << std::endl;
0272       validateCaloGeometry(DetId::Calo, HcalZDCDetId::SubdetectorId, "ZDC");
0273     } else
0274       fwLog(fwlog::kWarning) << "Invalid Calo geometry" << std::endl;
0275   }
0276 }
0277 
0278 void ValidateGeometry::validateRPCGeometry(const int regionNumber, const char* regionName) {
0279   clearData();
0280 
0281   std::vector<double> centers;
0282 
0283   auto const& rolls = rpcGeometry_->rolls();
0284 
0285   for (auto it = rolls.begin(), itEnd = rolls.end(); it != itEnd; ++it) {
0286     const RPCRoll* roll = *it;
0287 
0288     if (roll) {
0289       RPCDetId rpcDetId = roll->id();
0290 
0291       if (rpcDetId.region() == regionNumber) {
0292         const GeomDetUnit* det = rpcGeometry_->idToDetUnit(rpcDetId);
0293         GlobalPoint gp = det->surface().toGlobal(LocalPoint(0.0, 0.0, 0.0));
0294 
0295         const TGeoMatrix* matrix = fwGeometry_.getMatrix(rpcDetId.rawId());
0296 
0297         if (!matrix) {
0298           std::cout << "Failed to get matrix of RPC with detid: " << rpcDetId.rawId() << std::endl;
0299           continue;
0300         }
0301 
0302         compareTransform(gp, matrix);
0303 
0304         const float* shape = fwGeometry_.getShapePars(rpcDetId.rawId());
0305 
0306         if (!shape) {
0307           std::cout << "Failed to get shape of RPC with detid: " << rpcDetId.rawId() << std::endl;
0308           continue;
0309         }
0310 
0311         compareShape(det, shape);
0312 
0313         const float* parameters = fwGeometry_.getParameters(rpcDetId.rawId());
0314 
0315         if (parameters == nullptr) {
0316           std::cout << "Parameters empty for RPC with detid: " << rpcDetId.rawId() << std::endl;
0317           continue;
0318         }
0319 
0320         // Yes, I know that below I'm comparing the equivalence
0321         // of floating point numbers
0322 
0323         int nStrips = roll->nstrips();
0324         assert(nStrips == parameters[0]);
0325 
0326         float stripLength = roll->specificTopology().stripLength();
0327         assert(stripLength == parameters[1]);
0328 
0329         float pitch = roll->specificTopology().pitch();
0330         assert(pitch == parameters[2]);
0331 
0332         float offset = -0.5 * nStrips * pitch;
0333 
0334         for (int strip = 1; strip <= roll->nstrips(); ++strip) {
0335           LocalPoint centreOfStrip1 = roll->centreOfStrip(strip);
0336           LocalPoint centreOfStrip2 = LocalPoint((strip - 0.5) * pitch + offset, 0.0);
0337 
0338           centers.push_back(centreOfStrip1.x() - centreOfStrip2.x());
0339         }
0340       }
0341     }
0342   }
0343 
0344   std::string hn(regionName);
0345   makeHistogram(hn + ": centreOfStrip", centers);
0346 
0347   makeHistograms(regionName);
0348 }
0349 
0350 void ValidateGeometry::validateDTChamberGeometry() {
0351   clearData();
0352 
0353   auto const& chambers = dtGeometry_->chambers();
0354 
0355   for (auto it = chambers.begin(), itEnd = chambers.end(); it != itEnd; ++it) {
0356     const DTChamber* chamber = *it;
0357 
0358     if (chamber) {
0359       DTChamberId chId = chamber->id();
0360       GlobalPoint gp = chamber->surface().toGlobal(LocalPoint(0.0, 0.0, 0.0));
0361 
0362       const TGeoMatrix* matrix = fwGeometry_.getMatrix(chId.rawId());
0363 
0364       if (!matrix) {
0365         std::cout << "Failed to get matrix of DT chamber with detid: " << chId.rawId() << std::endl;
0366         continue;
0367       }
0368 
0369       compareTransform(gp, matrix);
0370 
0371       const float* shape = fwGeometry_.getShapePars(chId.rawId());
0372 
0373       if (!shape) {
0374         std::cout << "Failed to get shape of DT chamber with detid: " << chId.rawId() << std::endl;
0375         continue;
0376       }
0377 
0378       compareShape(chamber, shape);
0379     }
0380   }
0381 
0382   makeHistograms("DT chamber");
0383 }
0384 
0385 void ValidateGeometry::validateDTLayerGeometry() {
0386   clearData();
0387 
0388   std::vector<double> wire_positions;
0389 
0390   auto const& layers = dtGeometry_->layers();
0391 
0392   for (auto it = layers.begin(), itEnd = layers.end(); it != itEnd; ++it) {
0393     const DTLayer* layer = *it;
0394 
0395     if (layer) {
0396       DTLayerId layerId = layer->id();
0397       GlobalPoint gp = layer->surface().toGlobal(LocalPoint(0.0, 0.0, 0.0));
0398 
0399       const TGeoMatrix* matrix = fwGeometry_.getMatrix(layerId.rawId());
0400 
0401       if (!matrix) {
0402         std::cout << "Failed to get matrix of DT layer with detid: " << layerId.rawId() << std::endl;
0403         continue;
0404       }
0405 
0406       compareTransform(gp, matrix);
0407 
0408       const float* shape = fwGeometry_.getShapePars(layerId.rawId());
0409 
0410       if (!shape) {
0411         std::cout << "Failed to get shape of DT layer with detid: " << layerId.rawId() << std::endl;
0412         continue;
0413       }
0414 
0415       compareShape(layer, shape);
0416 
0417       const float* parameters = fwGeometry_.getParameters(layerId.rawId());
0418 
0419       if (parameters == nullptr) {
0420         std::cout << "Parameters empty for DT layer with detid: " << layerId.rawId() << std::endl;
0421         continue;
0422       }
0423 
0424       float width = layer->surface().bounds().width();
0425       assert(width == parameters[6]);
0426 
0427       float thickness = layer->surface().bounds().thickness();
0428       assert(thickness == parameters[7]);
0429 
0430       float length = layer->surface().bounds().length();
0431       assert(length == parameters[8]);
0432 
0433       int firstChannel = layer->specificTopology().firstChannel();
0434       assert(firstChannel == parameters[3]);
0435 
0436       int lastChannel = layer->specificTopology().lastChannel();
0437       int nChannels = parameters[5];
0438       assert(nChannels == (lastChannel - firstChannel) + 1);
0439 
0440       for (int wireN = firstChannel; wireN - lastChannel <= 0; ++wireN) {
0441         double localX1 = layer->specificTopology().wirePosition(wireN);
0442         double localX2 = (wireN - (firstChannel - 1) - 0.5) * parameters[0] - nChannels / 2.0 * parameters[0];
0443 
0444         wire_positions.push_back(localX1 - localX2);
0445 
0446         //std::cout<<"wireN, localXpos: "<< wireN <<" "<< localX1 <<" "<< localX2 <<std::endl;
0447       }
0448     }
0449   }
0450 
0451   makeHistogram("DT layer wire localX", wire_positions);
0452 
0453   makeHistograms("DT layer");
0454 }
0455 
0456 void ValidateGeometry::validateCSChamberGeometry(const int endcap, const char* detname) {
0457   clearData();
0458 
0459   auto const& chambers = cscGeometry_->chambers();
0460 
0461   for (auto it = chambers.begin(), itEnd = chambers.end(); it != itEnd; ++it) {
0462     const CSCChamber* chamber = *it;
0463 
0464     if (chamber && chamber->id().endcap() == endcap) {
0465       DetId detId = chamber->geographicalId();
0466       GlobalPoint gp = chamber->surface().toGlobal(LocalPoint(0.0, 0.0, 0.0));
0467 
0468       const TGeoMatrix* matrix = fwGeometry_.getMatrix(detId.rawId());
0469 
0470       if (!matrix) {
0471         std::cout << "Failed to get matrix of CSC chamber with detid: " << detId.rawId() << std::endl;
0472         continue;
0473       }
0474 
0475       compareTransform(gp, matrix);
0476 
0477       const float* shape = fwGeometry_.getShapePars(detId.rawId());
0478 
0479       if (!shape) {
0480         std::cout << "Failed to get shape of CSC chamber with detid: " << detId.rawId() << std::endl;
0481         continue;
0482       }
0483 
0484       compareShape(chamber, shape);
0485     }
0486   }
0487 
0488   makeHistograms(detname);
0489 }
0490 
0491 void ValidateGeometry::validateCSCLayerGeometry(const int endcap, const char* detname) {
0492   clearData();
0493   std::vector<double> strip_positions;
0494   std::vector<double> wire_positions;
0495 
0496   std::vector<double> me11_wiresLocal;
0497   std::vector<double> me12_wiresLocal;
0498   std::vector<double> me13_wiresLocal;
0499   std::vector<double> me14_wiresLocal;
0500   std::vector<double> me21_wiresLocal;
0501   std::vector<double> me22_wiresLocal;
0502   std::vector<double> me31_wiresLocal;
0503   std::vector<double> me32_wiresLocal;
0504   std::vector<double> me41_wiresLocal;
0505   std::vector<double> me42_wiresLocal;
0506 
0507   auto const& layers = cscGeometry_->layers();
0508 
0509   for (auto it = layers.begin(), itEnd = layers.end(); it != itEnd; ++it) {
0510     const CSCLayer* layer = *it;
0511 
0512     if (layer && layer->id().endcap() == endcap) {
0513       DetId detId = layer->geographicalId();
0514       GlobalPoint gp = layer->surface().toGlobal(LocalPoint(0.0, 0.0, 0.0));
0515 
0516       const TGeoMatrix* matrix = fwGeometry_.getMatrix(detId.rawId());
0517 
0518       if (!matrix) {
0519         std::cout << "Failed to get matrix of CSC layer with detid: " << detId.rawId() << std::endl;
0520         continue;
0521       }
0522 
0523       compareTransform(gp, matrix);
0524 
0525       const float* shape = fwGeometry_.getShapePars(detId.rawId());
0526 
0527       if (!shape) {
0528         std::cout << "Failed to get shape of CSC layer with detid: " << detId.rawId() << std::endl;
0529         continue;
0530       }
0531 
0532       compareShape(layer, shape);
0533 
0534       double length;
0535 
0536       if (shape[0] == 1) {
0537         length = shape[4];
0538       }
0539 
0540       else {
0541         std::cout << "Failed to get trapezoid from shape for CSC layer with detid: " << detId.rawId() << std::endl;
0542         continue;
0543       }
0544 
0545       const float* parameters = fwGeometry_.getParameters(detId.rawId());
0546 
0547       if (parameters == nullptr) {
0548         std::cout << "Parameters empty for CSC layer with detid: " << detId.rawId() << std::endl;
0549         continue;
0550       }
0551 
0552       int yAxisOrientation = layer->geometry()->topology()->yAxisOrientation();
0553       assert(yAxisOrientation == parameters[0]);
0554 
0555       float centreToIntersection = layer->geometry()->topology()->centreToIntersection();
0556       assert(centreToIntersection == parameters[1]);
0557 
0558       float yCentre = layer->geometry()->topology()->yCentreOfStripPlane();
0559       assert(yCentre == parameters[2]);
0560 
0561       float phiOfOneEdge = layer->geometry()->topology()->phiOfOneEdge();
0562       assert(phiOfOneEdge == parameters[3]);
0563 
0564       float stripOffset = layer->geometry()->topology()->stripOffset();
0565       assert(stripOffset == parameters[4]);
0566 
0567       float angularWidth = layer->geometry()->topology()->angularWidth();
0568       assert(angularWidth == parameters[5]);
0569 
0570       for (int nStrip = 1; nStrip <= layer->geometry()->numberOfStrips(); ++nStrip) {
0571         float xOfStrip1 = layer->geometry()->xOfStrip(nStrip);
0572 
0573         double stripAngle = phiOfOneEdge + yAxisOrientation * (nStrip - (0.5 - stripOffset)) * angularWidth;
0574         double xOfStrip2 = yAxisOrientation * (centreToIntersection - yCentre) * tan(stripAngle);
0575 
0576         strip_positions.push_back(xOfStrip1 - xOfStrip2);
0577       }
0578 
0579       int station = layer->id().station();
0580       int ring = layer->id().ring();
0581 
0582       double wireSpacingInGroup = layer->geometry()->wireTopology()->wireSpacing();
0583       assert(wireSpacingInGroup == parameters[6]);
0584 
0585       double wireSpacing = 0.0;
0586       // we calculate an average wire group
0587       // spacing from radialExtentOfTheWirePlane / numOfWireGroups
0588 
0589       double extentOfWirePlane = 0.0;
0590 
0591       if (ring == 2) {
0592         if (station == 1)
0593           extentOfWirePlane = 174.81;  //wireSpacing = 174.81/64;
0594         else
0595           extentOfWirePlane = 323.38;  //wireSpacing = 323.38/64;
0596       } else if (station == 1 && (ring == 1 || ring == 4))
0597         extentOfWirePlane = 150.5;  //wireSpacing = 150.5/48;
0598       else if (station == 1 && ring == 3)
0599         extentOfWirePlane = 164.47;  //wireSpacing = 164.47/32;
0600       else if (station == 2 && ring == 1)
0601         extentOfWirePlane = 189.97;  //wireSpacing = 189.97/112;
0602       else if (station == 3 && ring == 1)
0603         extentOfWirePlane = 170.01;  //wireSpacing = 170.01/96;
0604       else if (station == 4 && ring == 1)
0605         extentOfWirePlane = 149.73;  //wireSpacing = 149.73/96;
0606 
0607       float wireAngle = layer->geometry()->wireTopology()->wireAngle();
0608       assert(wireAngle == parameters[7]);
0609 
0610       //float cosWireAngle = cos(wireAngle);
0611 
0612       /* NOTE
0613          Some parameters don't seem available in a public interface
0614          so have to perhaps hard-code. This may not be too bad as there
0615          seems to be a lot of degeneracy. 
0616       */
0617 
0618       double alignmentPinToFirstWire;
0619       double yAlignmentFrame = 3.49;
0620 
0621       if (station == 1) {
0622         if (ring == 1 || ring == 4) {
0623           alignmentPinToFirstWire = 1.065;
0624           yAlignmentFrame = 0.0;
0625         }
0626 
0627         else  // ME12, ME 13
0628           alignmentPinToFirstWire = 2.85;
0629       }
0630 
0631       else if (station == 4 && ring == 1)
0632         alignmentPinToFirstWire = 3.04;
0633 
0634       else if (station == 3 && ring == 1)
0635         alignmentPinToFirstWire = 2.84;
0636 
0637       else  // ME21, ME22, ME32, ME42
0638         alignmentPinToFirstWire = 2.87;
0639 
0640       double yOfFirstWire = (yAlignmentFrame - length) + alignmentPinToFirstWire;
0641 
0642       int nWireGroups = layer->geometry()->numberOfWireGroups();
0643       double E = extentOfWirePlane / nWireGroups;
0644 
0645       for (int nWireGroup = 1; nWireGroup <= nWireGroups; ++nWireGroup) {
0646         LocalPoint centerOfWireGroup = layer->geometry()->localCenterOfWireGroup(nWireGroup);
0647         double yOfWire1 = centerOfWireGroup.y();
0648 
0649         //double yOfWire2 = (-0.5 - (nWireGroups*0.5 - 1) + (nWireGroup-1))*E;
0650         //yOfWire2 += 0.5*E;
0651         double yOfWire2 = yOfFirstWire + ((nWireGroup - 1) * E);
0652         yOfWire2 += wireSpacing * 0.5;
0653 
0654         double ydiff_local = yOfWire1 - yOfWire2;
0655         wire_positions.push_back(ydiff_local);
0656 
0657         //GlobalPoint globalPoint = layer->surface().toGlobal(LocalPoint(0.0,yOfWire1,0.0));
0658 
0659         /*
0660         float fwLocalPoint[3] = 
0661         {
0662           0.0, yOfWire2, 0.0
0663         };
0664         
0665         float fwGlobalPoint[3]; 
0666         fwGeometry_.localToGlobal(detId.rawId(), fwLocalPoint, fwGlobalPoint);
0667         double ydiff_global = globalPoint.y() - fwGlobalPoint[1]; 
0668         */
0669 
0670         if (station == 1) {
0671           if (ring == 1) {
0672             me11_wiresLocal.push_back(ydiff_local);
0673           } else if (ring == 2) {
0674             me12_wiresLocal.push_back(ydiff_local);
0675           } else if (ring == 3) {
0676             me13_wiresLocal.push_back(ydiff_local);
0677           } else if (ring == 4) {
0678             me14_wiresLocal.push_back(ydiff_local);
0679           }
0680         } else if (station == 2) {
0681           if (ring == 1) {
0682             me21_wiresLocal.push_back(ydiff_local);
0683           } else if (ring == 2) {
0684             me22_wiresLocal.push_back(ydiff_local);
0685           }
0686         } else if (station == 3) {
0687           if (ring == 1) {
0688             me31_wiresLocal.push_back(ydiff_local);
0689           } else if (ring == 2) {
0690             me32_wiresLocal.push_back(ydiff_local);
0691           }
0692         } else if (station == 4) {
0693           if (ring == 1) {
0694             me41_wiresLocal.push_back(ydiff_local);
0695           } else if (ring == 2) {
0696             me42_wiresLocal.push_back(ydiff_local);
0697           }
0698         }
0699       }
0700     }
0701   }
0702 
0703   std::string hn(detname);
0704   makeHistogram(hn + ": xOfStrip", strip_positions);
0705 
0706   makeHistogram(hn + ": local yOfWire", wire_positions);
0707 
0708   makeHistogram("ME11: local yOfWire", me11_wiresLocal);
0709   makeHistogram("ME12: local yOfWire", me12_wiresLocal);
0710   makeHistogram("ME13: local yOfWire", me13_wiresLocal);
0711   makeHistogram("ME14: local yOfWire", me14_wiresLocal);
0712   makeHistogram("ME21: local yOfWire", me21_wiresLocal);
0713   makeHistogram("ME22: local yOfWire", me22_wiresLocal);
0714   makeHistogram("ME31: local yOfWire", me31_wiresLocal);
0715   makeHistogram("ME32: local yOfWire", me32_wiresLocal);
0716   makeHistogram("ME41: local yOfWire", me41_wiresLocal);
0717   makeHistogram("ME42: local yOfWire", me42_wiresLocal);
0718 
0719   makeHistograms(detname);
0720 }
0721 
0722 void ValidateGeometry::validateCaloGeometry(DetId::Detector detector, int subdetector, const char* detname) {
0723   clearData();
0724 
0725   const CaloSubdetectorGeometry* geometry = caloGeometry_->getSubdetectorGeometry(detector, subdetector);
0726 
0727   const std::vector<DetId>& ids = geometry->getValidDetIds(detector, subdetector);
0728 
0729   for (auto it = ids.begin(), iEnd = ids.end(); it != iEnd; ++it) {
0730     unsigned int rawId = (*it).rawId();
0731 
0732     const float* points = fwGeometry_.getCorners(rawId);
0733 
0734     if (points == nullptr) {
0735       std::cout << "Failed to get points of " << detname << " element with detid: " << rawId << std::endl;
0736       continue;
0737     }
0738 
0739     auto cellGeometry = geometry->getGeometry(*it);
0740     const CaloCellGeometry::CornersVec& corners = cellGeometry->getCorners();
0741 
0742     assert(corners.size() == 8);
0743 
0744     for (unsigned int i = 0, offset = 0; i < 8; ++i) {
0745       offset = 2 * i;
0746 
0747       double distance = getDistance(GlobalPoint(points[i + offset], points[i + 1 + offset], points[i + 2 + offset]),
0748                                     GlobalPoint(corners[i].x(), corners[i].y(), corners[i].z()));
0749 
0750       globalDistances_.push_back(distance);
0751     }
0752   }
0753 
0754   makeHistograms(detname);
0755 }
0756 
0757 void ValidateGeometry::validateTrackerGeometry(const TrackerGeometry::DetContainer& dets, const char* detname) {
0758   clearData();
0759 
0760   for (TrackerGeometry::DetContainer::const_iterator it = dets.begin(), itEnd = dets.end(); it != itEnd; ++it) {
0761     GlobalPoint gp =
0762         (trackerGeometry_->idToDet((*it)->geographicalId()))->surface().toGlobal(LocalPoint(0.0, 0.0, 0.0));
0763     unsigned int rawId = (*it)->geographicalId().rawId();
0764 
0765     const TGeoMatrix* matrix = fwGeometry_.getMatrix(rawId);
0766 
0767     if (!matrix) {
0768       std::cout << "Failed to get matrix of " << detname << " element with detid: " << rawId << std::endl;
0769       continue;
0770     }
0771 
0772     compareTransform(gp, matrix);
0773 
0774     const float* shape = fwGeometry_.getShapePars(rawId);
0775 
0776     if (!shape) {
0777       std::cout << "Failed to get shape of " << detname << " element with detid: " << rawId << std::endl;
0778       continue;
0779     }
0780 
0781     compareShape(*it, shape);
0782   }
0783 
0784   makeHistograms(detname);
0785 }
0786 
0787 void ValidateGeometry::validatePixelTopology(const TrackerGeometry::DetContainer& dets, const char* detname) {
0788   std::vector<double> pixelLocalXs;
0789   std::vector<double> pixelLocalYs;
0790 
0791   for (TrackerGeometry::DetContainer::const_iterator it = dets.begin(), itEnd = dets.end(); it != itEnd; ++it) {
0792     unsigned int rawId = (*it)->geographicalId().rawId();
0793 
0794     const float* parameters = fwGeometry_.getParameters(rawId);
0795 
0796     if (parameters == nullptr) {
0797       std::cout << "Parameters empty for " << detname << " element with detid: " << rawId << std::endl;
0798       continue;
0799     }
0800 
0801     if (const PixelGeomDetUnit* det =
0802             dynamic_cast<const PixelGeomDetUnit*>(trackerGeometry_->idToDetUnit((*it)->geographicalId()))) {
0803       if (const PixelTopology* rpt = &det->specificTopology()) {
0804         int nrows = rpt->nrows();
0805         int ncolumns = rpt->ncolumns();
0806 
0807         for (int row = 1; row <= nrows; ++row) {
0808           for (int column = 1; column <= ncolumns; ++column) {
0809             LocalPoint localPoint = rpt->localPosition(MeasurementPoint(row, column));
0810 
0811             pixelLocalXs.push_back(localPoint.x() - fireworks::pixelLocalX(row, parameters));
0812             pixelLocalYs.push_back(localPoint.y() - fireworks::pixelLocalY(column, parameters));
0813           }
0814         }
0815       }
0816 
0817       else
0818         std::cout << "No topology for " << detname << " " << rawId << std::endl;
0819     }
0820 
0821     else
0822       std::cout << "No geomDetUnit for " << detname << " " << rawId << std::endl;
0823   }
0824 
0825   std::string hn(detname);
0826   makeHistogram(hn + " pixelLocalX", pixelLocalXs);
0827   makeHistogram(hn + " pixelLocalY", pixelLocalYs);
0828 }
0829 
0830 void ValidateGeometry::validateStripTopology(const TrackerGeometry::DetContainer& dets, const char* detname) {
0831   std::vector<double> radialStripLocalXs;
0832   std::vector<double> rectangularStripLocalXs;
0833 
0834   for (TrackerGeometry::DetContainer::const_iterator it = dets.begin(), itEnd = dets.end(); it != itEnd; ++it) {
0835     unsigned int rawId = (*it)->geographicalId().rawId();
0836 
0837     const float* parameters = fwGeometry_.getParameters(rawId);
0838 
0839     if (parameters == nullptr) {
0840       std::cout << "Parameters empty for " << detname << " element with detid: " << rawId << std::endl;
0841       continue;
0842     }
0843 
0844     if (const StripGeomDetUnit* det =
0845             dynamic_cast<const StripGeomDetUnit*>(trackerGeometry_->idToDet((*it)->geographicalId()))) {
0846       // NOTE: why the difference in dets vs. units between these and pixels? The dynamic cast above
0847       // fails for many of the detids...
0848 
0849       const StripTopology* st = dynamic_cast<const StripTopology*>(&det->specificTopology());
0850 
0851       if (st) {
0852         //assert(parameters[0] == 0);
0853         int nstrips = st->nstrips();
0854         assert(parameters[1] == nstrips);
0855         assert(parameters[2] == st->stripLength());
0856 
0857         if (const RadialStripTopology* rst =
0858                 dynamic_cast<const RadialStripTopology*>(&(det->specificType().specificTopology()))) {
0859           assert(parameters[0] == 1);
0860           assert(parameters[3] == rst->yAxisOrientation());
0861           assert(parameters[4] == rst->originToIntersection());
0862           assert(parameters[5] == rst->phiOfOneEdge());
0863           assert(parameters[6] == rst->angularWidth());
0864 
0865           for (uint16_t strip = 1; strip <= nstrips; ++strip) {
0866             float stripAngle1 = rst->stripAngle(strip);
0867             float stripAngle2 = parameters[3] * (parameters[5] + strip * parameters[6]);
0868 
0869             assert((stripAngle1 - stripAngle2) == 0);
0870 
0871             LocalPoint stripPosition = st->localPosition(strip);
0872 
0873             float stripX = parameters[4] * tan(stripAngle2);
0874             radialStripLocalXs.push_back(stripPosition.x() - stripX);
0875           }
0876         }
0877 
0878         else if (dynamic_cast<const RectangularStripTopology*>(&(det->specificType().specificTopology()))) {
0879           assert(parameters[0] == 2);
0880           assert(parameters[3] == st->pitch());
0881 
0882           for (uint16_t strip = 1; strip <= nstrips; ++strip) {
0883             LocalPoint stripPosition = st->localPosition(strip);
0884             float stripX = -parameters[1] * 0.5 * parameters[3];
0885             stripX += strip * parameters[3];
0886             rectangularStripLocalXs.push_back(stripPosition.x() - stripX);
0887           }
0888         }
0889 
0890         else if (dynamic_cast<const TrapezoidalStripTopology*>(&(det->specificType().specificTopology()))) {
0891           assert(parameters[0] == 3);
0892           assert(parameters[3] == st->pitch());
0893         }
0894 
0895         else
0896           std::cout << "Failed to get pitch for " << detname << " " << rawId << std::endl;
0897       }
0898 
0899       else
0900         std::cout << "Failed cast to StripTopology for " << detname << " " << rawId << std::endl;
0901     }
0902 
0903     //else
0904     //  std::cout<<"Failed cast to StripGeomDetUnit for "<< detname <<" "<< rawId <<std::endl;
0905   }
0906 
0907   std::string hn(detname);
0908   makeHistogram(hn + " radial strip localX", radialStripLocalXs);
0909   makeHistogram(hn + " rectangular strip localX", rectangularStripLocalXs);
0910 }
0911 
0912 void ValidateGeometry::compareTransform(const GlobalPoint& gp, const TGeoMatrix* matrix) {
0913   double local[3] = {0.0, 0.0, 0.0};
0914 
0915   double global[3];
0916 
0917   matrix->LocalToMaster(local, global);
0918 
0919   double distance = getDistance(GlobalPoint(global[0], global[1], global[2]), gp);
0920   globalDistances_.push_back(distance);
0921 }
0922 
0923 void ValidateGeometry::compareShape(const GeomDet* det, const float* shape) {
0924   double shape_topWidth;
0925   double shape_bottomWidth;
0926   double shape_length;
0927   double shape_thickness;
0928 
0929   if (shape[0] == 1) {
0930     shape_topWidth = shape[2];
0931     shape_bottomWidth = shape[1];
0932     shape_length = shape[4];
0933     shape_thickness = shape[3];
0934   }
0935 
0936   else if (shape[0] == 2) {
0937     shape_topWidth = shape[1];
0938     shape_bottomWidth = shape[1];
0939     shape_length = shape[2];
0940     shape_thickness = shape[3];
0941   }
0942 
0943   else {
0944     std::cout << "Failed to get box or trapezoid from shape" << std::endl;
0945     return;
0946   }
0947 
0948   double topWidth, bottomWidth;
0949   double length, thickness;
0950 
0951   const Bounds* bounds = &(det->surface().bounds());
0952 
0953   if (const TrapezoidalPlaneBounds* tpbs = dynamic_cast<const TrapezoidalPlaneBounds*>(bounds)) {
0954     std::array<const float, 4> const& ps = tpbs->parameters();
0955 
0956     assert(ps.size() == 4);
0957 
0958     bottomWidth = ps[0];
0959     topWidth = ps[1];
0960     thickness = ps[2];
0961     length = ps[3];
0962   }
0963 
0964   else if ((dynamic_cast<const RectangularPlaneBounds*>(bounds))) {
0965     length = det->surface().bounds().length() * 0.5;
0966     topWidth = det->surface().bounds().width() * 0.5;
0967     bottomWidth = topWidth;
0968     thickness = det->surface().bounds().thickness() * 0.5;
0969   }
0970 
0971   else {
0972     std::cout << "Failed to get bounds" << std::endl;
0973     return;
0974   }
0975 
0976   //assert((tgeotrap && trapezoid) || (tgeobbox && rectangle));
0977 
0978   /*
0979   std::cout<<"topWidth: "<< shape_topWidth <<" "<< topWidth <<std::endl;
0980   std::cout<<"bottomWidth: "<< shape_bottomWidth <<" "<< bottomWidth <<std::endl;
0981   std::cout<<"length: "<< shape_length <<" "<< length <<std::endl;
0982   std::cout<<"thickness: "<< shape_thickness <<" "<< thickness <<std::endl;
0983   */
0984 
0985   topWidths_.push_back(fabs(shape_topWidth - topWidth));
0986   bottomWidths_.push_back(fabs(shape_bottomWidth - bottomWidth));
0987   lengths_.push_back(fabs(shape_length - length));
0988   thicknesses_.push_back(fabs(shape_thickness - thickness));
0989 
0990   return;
0991 }
0992 
0993 double ValidateGeometry::getDistance(const GlobalPoint& p1, const GlobalPoint& p2) {
0994   /*
0995   std::cout<<"X: "<< p1.x() <<" "<< p2.x() <<std::endl;
0996   std::cout<<"Y: "<< p1.y() <<" "<< p2.y() <<std::endl;
0997   std::cout<<"Z: "<< p1.z() <<" "<< p2.z() <<std::endl;
0998   */
0999 
1000   return sqrt((p1.x() - p2.x()) * (p1.x() - p2.x()) + (p1.y() - p2.y()) * (p1.y() - p2.y()) +
1001               (p1.z() - p2.z()) * (p1.z() - p2.z()));
1002 }
1003 
1004 void ValidateGeometry::makeHistograms(const char* detector) {
1005   outFile_->cd();
1006 
1007   std::string d(detector);
1008 
1009   std::string gdn = d + ": distance between points in global coordinates";
1010   makeHistogram(gdn, globalDistances_);
1011 
1012   std::string twn = d + ": absolute difference between top widths (along X)";
1013   makeHistogram(twn, topWidths_);
1014 
1015   std::string bwn = d + ": absolute difference between bottom widths (along X)";
1016   makeHistogram(bwn, bottomWidths_);
1017 
1018   std::string ln = d + ": absolute difference between lengths (along Y)";
1019   makeHistogram(ln, lengths_);
1020 
1021   std::string tn = d + ": absolute difference between thicknesses (along Z)";
1022   makeHistogram(tn, thicknesses_);
1023 
1024   return;
1025 }
1026 
1027 void ValidateGeometry::makeHistogram(const std::string& name, std::vector<double>& data) {
1028   if (data.empty())
1029     return;
1030 
1031   std::vector<double>::iterator it;
1032 
1033   it = std::min_element(data.begin(), data.end());
1034   double minE = *it;
1035 
1036   it = std::max_element(data.begin(), data.end());
1037   double maxE = *it;
1038 
1039   std::vector<double>::iterator itEnd = data.end();
1040 
1041   TH1D hist(name.c_str(), name.c_str(), 100, minE * (1 + 0.10), maxE * (1 + 0.10));
1042 
1043   for (it = data.begin(); it != itEnd; ++it)
1044     hist.Fill(*it);
1045 
1046   hist.GetXaxis()->SetTitle("[cm]");
1047   hist.Write();
1048 }
1049 
1050 void ValidateGeometry::beginJob() { outFile_->cd(); }
1051 
1052 void ValidateGeometry::endJob() {
1053   std::cout << "Done. " << std::endl;
1054   std::cout << "Results written to " << outfileName_ << std::endl;
1055   outFile_->Close();
1056 }
1057 
1058 DEFINE_FWK_MODULE(ValidateGeometry);