File indexing completed on 2023-03-17 13:03:52
0001
0002
0003
0004
0005
0006
0007
0008
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/RPCGeometry/interface/RPCGeometry.h"
0018 #include "Geometry/RPCGeometry/interface/RPCChamber.h"
0019 #include "Geometry/RPCGeometry/interface/RPCRoll.h"
0020 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0021 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0022
0023 #include "Fireworks/Core/interface/FWGeometry.h"
0024
0025 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0026 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0027
0028 #include <TFile.h>
0029 #include <TH1.h>
0030
0031 #include <limits>
0032 #include <string>
0033 #include <type_traits>
0034 #include <algorithm>
0035
0036 using namespace std;
0037
0038 template <class T>
0039 typename enable_if<!numeric_limits<T>::is_integer, bool>::type almost_equal(T x, T y, int ulp) {
0040
0041
0042 return abs(x - y) <= numeric_limits<T>::epsilon() * abs(x + y) * ulp
0043
0044 || abs(x - y) < numeric_limits<T>::min();
0045 }
0046
0047 using namespace edm;
0048
0049 class RPCGeometryValidate : public one::EDAnalyzer<> {
0050 public:
0051 explicit RPCGeometryValidate(const ParameterSet&);
0052 ~RPCGeometryValidate() 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 validateRPCChamberGeometry();
0060 void validateRPCStripsGeometry();
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<RPCGeometry, MuonGeometryRecord> tokRPC_;
0087 const RPCGeometry* rpcGeometry_;
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 RPCGeometryValidate::RPCGeometryValidate(const edm::ParameterSet& iConfig)
0104 : tokRPC_{esConsumes<RPCGeometry, MuonGeometryRecord>(edm::ESInputTag{})},
0105 infileName_(iConfig.getUntrackedParameter<string>("infileName", "cmsGeom10.root")),
0106 outfileName_(iConfig.getUntrackedParameter<string>("outfileName", "validateRPCGeometry.root")),
0107 tolerance_(iConfig.getUntrackedParameter<int>("tolerance", 6)) {
0108 fwGeometry_.loadMap(infileName_.c_str());
0109 outFile_ = new TFile(outfileName_.c_str(), "RECREATE");
0110 }
0111
0112 void RPCGeometryValidate::analyze(const edm::Event& event, const edm::EventSetup& eventSetup) {
0113 rpcGeometry_ = &eventSetup.getData(tokRPC_);
0114 LogVerbatim("RPCGeometry") << "Validating RPC chamber geometry";
0115 validateRPCChamberGeometry();
0116 validateRPCStripsGeometry();
0117 }
0118
0119 void RPCGeometryValidate::validateRPCChamberGeometry() {
0120 clearData();
0121
0122 for (auto const& it : rpcGeometry_->rolls()) {
0123 RPCDetId chId = it->id();
0124 GlobalPoint gp = it->surface().toGlobal(LocalPoint(0.0, 0.0, 0.0));
0125
0126 const TGeoMatrix* matrix = fwGeometry_.getMatrix(chId.rawId());
0127
0128 if (!matrix) {
0129 LogVerbatim("RPCGeometry") << "Failed to get matrix of RPC chamber with detid: " << chId.rawId();
0130 continue;
0131 }
0132 compareTransform(gp, matrix);
0133
0134 auto const& shape = fwGeometry_.getShapePars(chId.rawId());
0135
0136 if (!shape) {
0137 LogVerbatim("RPCGeometry") << "Failed to get shape of RPC chamber with detid: " << chId.rawId();
0138 continue;
0139 }
0140 compareShape(it, shape);
0141 }
0142 makeHistograms("RPC Chamber");
0143 }
0144
0145 void RPCGeometryValidate::validateRPCStripsGeometry() {
0146 clearData2();
0147
0148 for (auto const& it : rpcGeometry_->rolls()) {
0149 RPCDetId chId = it->id();
0150 const int n_strips = it->nstrips();
0151 const float n_pitch = it->pitch();
0152 const StripTopology& topo = it->specificTopology();
0153 const float stripLen = topo.stripLength();
0154 const float* parameters = fwGeometry_.getParameters(chId.rawId());
0155
0156 if (n_strips) {
0157 for (int istrips = 1; istrips <= n_strips; istrips++) {
0158 nstrips_.push_back(fabs(n_strips - parameters[0]));
0159 pitch_.push_back(fabs(n_pitch - parameters[2]));
0160 stripslen_.push_back(fabs(stripLen - parameters[1]));
0161 }
0162 } else {
0163 LogVerbatim("RPCGeometry") << "ATTENTION! nStrips == 0";
0164 }
0165 }
0166 makeHistograms2("RPC Strips");
0167 }
0168
0169 void RPCGeometryValidate::compareTransform(const GlobalPoint& gp, const TGeoMatrix* matrix) {
0170 double local[3] = {0.0, 0.0, 0.0};
0171 double global[3];
0172
0173 matrix->LocalToMaster(local, global);
0174
0175 float distance = getDistance(GlobalPoint(global[0], global[1], global[2]), gp);
0176 if ((distance >= 0.0) && (distance < 1.0e-7))
0177 distance = 0.0;
0178 globalDistances_.push_back(distance);
0179 }
0180
0181 void RPCGeometryValidate::compareShape(const GeomDet* det, const float* shape) {
0182 float shapeTopWidth;
0183 float shapeBottomWidth;
0184 float shapeLength;
0185 float shapeThickness;
0186
0187 if (shape[0] == 1) {
0188 shapeTopWidth = shape[2];
0189 shapeBottomWidth = shape[1];
0190 shapeLength = shape[4];
0191 shapeThickness = shape[3];
0192 } else if (shape[0] == 2) {
0193 shapeTopWidth = shape[1];
0194 shapeBottomWidth = shape[1];
0195 shapeLength = shape[2];
0196 shapeThickness = shape[3];
0197 } else {
0198 LogVerbatim("RPCGeometry") << "Failed to get box or trapezoid from shape";
0199
0200 return;
0201 }
0202
0203 float topWidth, bottomWidth;
0204 float length, thickness;
0205
0206 const Bounds* bounds = &(det->surface().bounds());
0207 if (const TrapezoidalPlaneBounds* tpbs = dynamic_cast<const TrapezoidalPlaneBounds*>(bounds)) {
0208 array<const float, 4> const& ps = tpbs->parameters();
0209
0210 assert(ps.size() == 4);
0211
0212 bottomWidth = ps[0];
0213 topWidth = ps[1];
0214 thickness = ps[2];
0215 length = ps[3];
0216 } else if ((dynamic_cast<const RectangularPlaneBounds*>(bounds))) {
0217 length = det->surface().bounds().length() * 0.5;
0218 topWidth = det->surface().bounds().width() * 0.5;
0219 bottomWidth = topWidth;
0220 thickness = det->surface().bounds().thickness() * 0.5;
0221 } else {
0222 LogVerbatim("RPCGeometry") << "Failed to get bounds";
0223
0224 return;
0225 }
0226 topWidths_.push_back(fabs(shapeTopWidth - topWidth));
0227 bottomWidths_.push_back(fabs(shapeBottomWidth - bottomWidth));
0228 lengths_.push_back(fabs(shapeLength - length));
0229 thicknesses_.push_back(fabs(shapeThickness - thickness));
0230 }
0231
0232 float RPCGeometryValidate::getDistance(const GlobalPoint& p1, const GlobalPoint& p2) {
0233 return sqrt((p1.x() - p2.x()) * (p1.x() - p2.x()) + (p1.y() - p2.y()) * (p1.y() - p2.y()) +
0234 (p1.z() - p2.z()) * (p1.z() - p2.z()));
0235 }
0236
0237 float RPCGeometryValidate::getDiff(const float val1, const float val2) {
0238 if (almost_equal(val1, val2, tolerance_))
0239 return 0.0f;
0240 else
0241 return (val1 - val2);
0242 }
0243
0244 void RPCGeometryValidate::makeHistograms(const char* detector) {
0245 outFile_->cd();
0246
0247 string d(detector);
0248
0249 string gdn = d + ": distance between points in global coordinates";
0250 makeHistogram(gdn, globalDistances_);
0251
0252 string twn = d + ": absolute difference between top widths (along X)";
0253 makeHistogram(twn, topWidths_);
0254
0255 string bwn = d + ": absolute difference between bottom widths (along X)";
0256 makeHistogram(bwn, bottomWidths_);
0257
0258 string ln = d + ": absolute difference between lengths (along Y)";
0259 makeHistogram(ln, lengths_);
0260
0261 string tn = d + ": absolute difference between thicknesses (along Z)";
0262 makeHistogram(tn, thicknesses_);
0263 }
0264
0265 void RPCGeometryValidate::makeHistograms2(const char* detector) {
0266 outFile_->cd();
0267
0268 string d(detector);
0269
0270 string ns = d + ": absolute difference between nStrips";
0271 makeHistogram(ns, nstrips_);
0272
0273 string pi = d + ": absolute difference between Strips Pitch";
0274 makeHistogram(pi, pitch_);
0275
0276 string pl = d + ": absolute difference between Strips Length";
0277 makeHistogram(pl, stripslen_);
0278 }
0279
0280 void RPCGeometryValidate::makeHistogram(const string& name, vector<float>& data) {
0281 if (data.empty())
0282 return;
0283
0284 const auto [minE, maxE] = minmax_element(begin(data), end(data));
0285
0286 TH1D hist(name.c_str(), name.c_str(), 100, *minE * (1 + 0.10), *maxE * (1 + 0.10));
0287
0288 for (auto const& it : data)
0289 hist.Fill(it);
0290
0291 hist.GetXaxis()->SetTitle("[cm]");
0292 hist.Write();
0293 }
0294
0295 void RPCGeometryValidate::beginJob() { outFile_->cd(); }
0296
0297 void RPCGeometryValidate::endJob() {
0298 LogVerbatim("RPCGeometry") << "Done.";
0299 LogVerbatim("RPCGeometry") << "Results written to " << outfileName_;
0300 outFile_->Close();
0301 }
0302
0303 DEFINE_FWK_MODULE(RPCGeometryValidate);