File indexing completed on 2024-04-06 12:11:47
0001
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
0321
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
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
0587
0588
0589 double extentOfWirePlane = 0.0;
0590
0591 if (ring == 2) {
0592 if (station == 1)
0593 extentOfWirePlane = 174.81;
0594 else
0595 extentOfWirePlane = 323.38;
0596 } else if (station == 1 && (ring == 1 || ring == 4))
0597 extentOfWirePlane = 150.5;
0598 else if (station == 1 && ring == 3)
0599 extentOfWirePlane = 164.47;
0600 else if (station == 2 && ring == 1)
0601 extentOfWirePlane = 189.97;
0602 else if (station == 3 && ring == 1)
0603 extentOfWirePlane = 170.01;
0604 else if (station == 4 && ring == 1)
0605 extentOfWirePlane = 149.73;
0606
0607 float wireAngle = layer->geometry()->wireTopology()->wireAngle();
0608 assert(wireAngle == parameters[7]);
0609
0610
0611
0612
0613
0614
0615
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
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
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
0650
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
0658
0659
0660
0661
0662
0663
0664
0665
0666
0667
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
0847
0848
0849 const StripTopology* st = dynamic_cast<const StripTopology*>(&det->specificTopology());
0850
0851 if (st) {
0852
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
0904
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
0977
0978
0979
0980
0981
0982
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
0996
0997
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);
