File indexing completed on 2021-11-09 23:32:01
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include "Geometry/RPCGeometryBuilder/src/RPCGeometryBuilder.h"
0012 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0013 #include "Geometry/RPCGeometry/interface/RPCRollSpecs.h"
0014 #include <DetectorDescription/Core/interface/DDFilter.h>
0015 #include <DetectorDescription/Core/interface/DDFilteredView.h>
0016 #include <DetectorDescription/DDCMS/interface/DDFilteredView.h>
0017 #include <DetectorDescription/DDCMS/interface/DDCompactView.h>
0018 #include <DetectorDescription/Core/interface/DDSolid.h>
0019 #include "Geometry/MuonNumbering/interface/MuonGeometryNumbering.h"
0020 #include "Geometry/MuonNumbering/interface/MuonBaseNumber.h"
0021 #include "Geometry/MuonNumbering/interface/RPCNumberingScheme.h"
0022 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
0023 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0024 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0025 #include "DataFormats/GeometryVector/interface/Basic3DVector.h"
0026 #include <iostream>
0027 #include <algorithm>
0028 #include "DetectorDescription/DDCMS/interface/DDSpecParRegistry.h"
0029 #include "DataFormats/Math/interface/CMSUnits.h"
0030 #include "DataFormats/Math/interface/GeantUnits.h"
0031 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0032
0033 using namespace cms_units::operators;
0034
0035 RPCGeometryBuilder::RPCGeometryBuilder() {}
0036
0037
0038
0039 std::unique_ptr<RPCGeometry> RPCGeometryBuilder::build(const DDCompactView* cview,
0040 const MuonGeometryConstants& muonConstants) {
0041 const std::string attribute = "ReadOutName";
0042 const std::string value = "MuonRPCHits";
0043 DDSpecificsMatchesValueFilter filter{DDValue(attribute, value, 0.0)};
0044 DDFilteredView fview(*cview, filter);
0045 return this->buildGeometry(fview, muonConstants);
0046 }
0047
0048
0049 std::unique_ptr<RPCGeometry> RPCGeometryBuilder::build(const cms::DDCompactView* cview,
0050 const MuonGeometryConstants& muonConstants) {
0051 const std::string attribute = "ReadOutName";
0052 const std::string value = "MuonRPCHits";
0053 const cms::DDFilter filter(attribute, value);
0054 cms::DDFilteredView fview(*cview, filter);
0055 return this->buildGeometry(fview, muonConstants);
0056 }
0057
0058
0059 std::unique_ptr<RPCGeometry> RPCGeometryBuilder::buildGeometry(DDFilteredView& fview,
0060 const MuonGeometryConstants& muonConstants) {
0061 LogDebug("RPCGeometryBuilder") << "Building the geometry service";
0062 std::unique_ptr<RPCGeometry> geometry = std::make_unique<RPCGeometry>();
0063 LogDebug("RPCGeometryBuilder") << "About to run through the RPC structure\n"
0064 << " First logical part " << fview.logicalPart().name().name();
0065 bool doSubDets = fview.firstChild();
0066 LogDebug("RPCGeometryBuilder") << "doSubDets = " << doSubDets;
0067 while (doSubDets) {
0068 LogDebug("RPCGeometryBuilder") << "start the loop";
0069 MuonGeometryNumbering mdddnum(muonConstants);
0070 LogDebug("RPCGeometryBuilder") << "Getting the Muon base Number";
0071 MuonBaseNumber mbn = mdddnum.geoHistoryToBaseNumber(fview.geoHistory());
0072 LogDebug("RPCGeometryBuilder") << "Start the Rpc Numbering Schema";
0073 RPCNumberingScheme rpcnum(muonConstants);
0074 LogDebug("RPCGeometryBuilder") << "Getting the Unit Number";
0075 const int detid = rpcnum.baseNumberToUnitNumber(mbn);
0076 LogDebug("RPCGeometryBuilder") << "Getting the RPC det Id " << detid;
0077 RPCDetId rpcid(detid);
0078 RPCDetId chid(rpcid.region(), rpcid.ring(), rpcid.station(), rpcid.sector(), rpcid.layer(), rpcid.subsector(), 0);
0079 LogDebug("RPCGeometryBuilder") << "The RPCDetid is " << rpcid;
0080
0081 DDValue numbOfStrips("nStrips");
0082
0083 std::vector<const DDsvalues_type*> specs(fview.specifics());
0084 int nStrips = 0;
0085 for (auto& spec : specs) {
0086 if (DDfetch(spec, numbOfStrips)) {
0087 nStrips = int(numbOfStrips.doubles()[0]);
0088 }
0089 }
0090
0091 LogDebug("RPCGeometryBuilder") << ((nStrips == 0) ? ("No strip found!!") : (""));
0092
0093 std::vector<double> dpar = fview.logicalPart().solid().parameters();
0094 std::string name = fview.logicalPart().name().name();
0095
0096 edm::LogVerbatim("RPCGeometryBuilder")
0097 << "(1) "
0098 << "detid: " << detid << " name: " << name << " number of Strips: " << nStrips;
0099
0100 DDTranslation tran = fview.translation();
0101 DDRotationMatrix rota = fview.rotation();
0102 Surface::PositionType pos(geant_units::operators::convertMmToCm(tran.x()),
0103 geant_units::operators::convertMmToCm(tran.y()),
0104 geant_units::operators::convertMmToCm(tran.z()));
0105 edm::LogVerbatim("RPCGeometryBuilder") << "(2), tran.x() " << geant_units::operators::convertMmToCm(tran.x())
0106 << " tran.y(): " << geant_units::operators::convertMmToCm(tran.y())
0107 << " tran.z(): " << geant_units::operators::convertMmToCm(tran.z());
0108
0109 DD3Vector x, y, z;
0110 rota.GetComponents(x, y, z);
0111 Surface::RotationType rot(float(x.X()),
0112 float(x.Y()),
0113 float(x.Z()),
0114 float(y.X()),
0115 float(y.Y()),
0116 float(y.Z()),
0117 float(z.X()),
0118 float(z.Y()),
0119 float(z.Z()));
0120
0121 RPCRollSpecs* rollspecs = nullptr;
0122 Bounds* bounds = nullptr;
0123
0124 if (dpar.size() == 3) {
0125 const float width = geant_units::operators::convertMmToCm(dpar[0]);
0126 const float length = geant_units::operators::convertMmToCm(dpar[1]);
0127 const float thickness = geant_units::operators::convertMmToCm(dpar[2]);
0128
0129
0130 edm::LogVerbatim("RPCGeometryBuilder")
0131 << "(3) dpar.size() == 3, width: " << width << " length: " << length << " thickness: " << thickness;
0132
0133 bounds = new RectangularPlaneBounds(width, length, thickness);
0134 const std::vector<float> pars = {width, length, float(numbOfStrips.doubles()[0])};
0135
0136 rollspecs = new RPCRollSpecs(GeomDetEnumerators::RPCBarrel, name, pars);
0137
0138 } else {
0139 const float be = geant_units::operators::convertMmToCm(dpar[4]);
0140 const float te = geant_units::operators::convertMmToCm(dpar[8]);
0141 const float ap = geant_units::operators::convertMmToCm(dpar[0]);
0142 const float ti = 0.4;
0143
0144 bounds = new TrapezoidalPlaneBounds(be, te, ap, ti);
0145
0146 const std::vector<float> pars = {float(geant_units::operators::convertMmToCm(dpar[4])),
0147 float(geant_units::operators::convertMmToCm(dpar[8])),
0148 float(geant_units::operators::convertMmToCm(dpar[0])),
0149 float(numbOfStrips.doubles()[0])};
0150
0151 edm::LogVerbatim("RPCGeometryBuilder")
0152 << "(4), else, dpar[4]: " << be << " dpar[8]: " << te << " dpar[0]: " << ap << " ti: " << ti;
0153
0154 rollspecs = new RPCRollSpecs(GeomDetEnumerators::RPCEndcap, name, pars);
0155
0156 Basic3DVector<float> newX(1., 0., 0.);
0157 Basic3DVector<float> newY(0., 0., 1.);
0158 newY *= -1;
0159 Basic3DVector<float> newZ(0., 1., 0.);
0160 rot.rotateAxes(newX, newY, newZ);
0161 }
0162 LogDebug("RPCGeometryBuilder") << " Number of strips " << nStrips;
0163
0164 BoundPlane* bp = new BoundPlane(pos, rot, bounds);
0165 ReferenceCountingPointer<BoundPlane> surf(bp);
0166 RPCRoll* r = new RPCRoll(rpcid, surf, rollspecs);
0167 geometry->add(r);
0168
0169 auto rls = chids.find(chid);
0170 if (rls == chids.end())
0171 rls = chids.insert(std::make_pair(chid, std::list<RPCRoll*>())).first;
0172 rls->second.emplace_back(r);
0173
0174 doSubDets = fview.nextSibling();
0175 }
0176 for (auto& ich : chids) {
0177 const RPCDetId& chid = ich.first;
0178 const auto& rls = ich.second;
0179
0180 BoundPlane* bp = nullptr;
0181 if (!rls.empty()) {
0182 const auto& refSurf = (*rls.begin())->surface();
0183 if (chid.region() == 0) {
0184 float corners[6] = {0, 0, 0, 0, 0, 0};
0185 for (auto rl : rls) {
0186 const double h2 = rl->surface().bounds().length() / 2;
0187 const double w2 = rl->surface().bounds().width() / 2;
0188 const auto x1y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-w2, -h2, 0)));
0189 const auto x2y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+w2, +h2, 0)));
0190 corners[0] = std::min(corners[0], x1y1AtRef.x());
0191 corners[1] = std::min(corners[1], x1y1AtRef.y());
0192 corners[2] = std::max(corners[2], x2y2AtRef.x());
0193 corners[3] = std::max(corners[3], x2y2AtRef.y());
0194 corners[4] = std::min(corners[4], x1y1AtRef.z());
0195 corners[5] = std::max(corners[5], x1y1AtRef.z());
0196 }
0197 const LocalPoint lpOfCentre((corners[0] + corners[2]) / 2, (corners[1] + corners[3]) / 2, 0);
0198 const auto gpOfCentre = refSurf.toGlobal(lpOfCentre);
0199 auto bounds = new RectangularPlaneBounds(
0200 (corners[2] - corners[0]) / 2, (corners[3] - corners[1]) / 2, (corners[5] - corners[4]) + 0.5);
0201 bp = new BoundPlane(gpOfCentre, refSurf.rotation(), bounds);
0202 } else {
0203 float cornersLo[3] = {0, 0, 0}, cornersHi[3] = {0, 0, 0};
0204 float cornersZ[2] = {0, 0};
0205 for (auto rl : rls) {
0206 const double h2 = rl->surface().bounds().length() / 2;
0207 const double w2 = rl->surface().bounds().width() / 2;
0208 const auto& topo = dynamic_cast<const TrapezoidalStripTopology&>(rl->specificTopology());
0209 const double r = topo.radius();
0210 const double wAtLo = w2 / r * (r - h2);
0211 const double wAtHi = w2 / r * (r + h2);
0212
0213 const auto x1y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-wAtLo, -h2, 0)));
0214 const auto x2y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+wAtLo, -h2, 0)));
0215 const auto x1y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-wAtHi, +h2, 0)));
0216 const auto x2y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+wAtHi, +h2, 0)));
0217
0218 cornersLo[0] = std::min(cornersLo[0], x1y1AtRef.x());
0219 cornersLo[1] = std::max(cornersLo[1], x2y1AtRef.x());
0220 cornersLo[2] = std::min(cornersLo[2], x1y1AtRef.y());
0221
0222 cornersHi[0] = std::min(cornersHi[0], x1y2AtRef.x());
0223 cornersHi[1] = std::max(cornersHi[1], x2y2AtRef.x());
0224 cornersHi[2] = std::max(cornersHi[2], x1y2AtRef.y());
0225
0226 cornersZ[0] = std::min(cornersZ[0], x1y1AtRef.z());
0227 cornersZ[1] = std::max(cornersZ[1], x1y1AtRef.z());
0228 }
0229 const LocalPoint lpOfCentre((cornersHi[0] + cornersHi[1]) / 2, (cornersLo[2] + cornersHi[2]) / 2, 0);
0230 const auto gpOfCentre = refSurf.toGlobal(lpOfCentre);
0231 auto bounds = new TrapezoidalPlaneBounds((cornersLo[1] - cornersLo[0]) / 2,
0232 (cornersHi[1] - cornersHi[0]) / 2,
0233 (cornersHi[2] - cornersLo[2]) / 2,
0234 (cornersZ[1] - cornersZ[0]) + 0.5);
0235 bp = new BoundPlane(gpOfCentre, refSurf.rotation(), bounds);
0236 }
0237 }
0238
0239 ReferenceCountingPointer<BoundPlane> surf(bp);
0240 RPCChamber* ch = new RPCChamber(chid, surf);
0241 for (auto rl : rls)
0242 ch->add(rl);
0243 geometry->add(ch);
0244 }
0245
0246 return geometry;
0247 }
0248
0249
0250
0251 std::unique_ptr<RPCGeometry> RPCGeometryBuilder::buildGeometry(cms::DDFilteredView& fview,
0252 const MuonGeometryConstants& muonConstants) {
0253 std::unique_ptr<RPCGeometry> geometry = std::make_unique<RPCGeometry>();
0254
0255 while (fview.firstChild()) {
0256 MuonGeometryNumbering mdddnum(muonConstants);
0257 RPCNumberingScheme rpcnum(muonConstants);
0258 int rawidCh = rpcnum.baseNumberToUnitNumber(mdddnum.geoHistoryToBaseNumber(fview.history()));
0259 RPCDetId rpcid = RPCDetId(rawidCh);
0260
0261 RPCDetId chid(rpcid.region(), rpcid.ring(), rpcid.station(), rpcid.sector(), rpcid.layer(), rpcid.subsector(), 0);
0262
0263 auto nStrips = fview.get<double>("nStrips");
0264
0265 std::vector<double> dpar = fview.parameters();
0266
0267 std::string_view name = fview.name();
0268
0269 edm::LogVerbatim("RPCGeometryBuilder")
0270 << "(1), detid: " << rawidCh << " name: " << std::string(name) << " number of Strips: " << nStrips;
0271
0272 const Double_t* tran = fview.trans();
0273
0274 DDRotationMatrix rota;
0275 fview.rot(rota);
0276
0277 Surface::PositionType pos(tran[0] / dd4hep::cm, tran[1] / dd4hep::cm, tran[2] / dd4hep::cm);
0278 edm::LogVerbatim("RPCGeometryBuilder")
0279 << "(2), tran.x(): " << tran[0] / dd4hep::cm << " tran.y(): " << tran[1] / dd4hep::cm
0280 << " tran.z(): " << tran[2] / dd4hep::cm;
0281 DD3Vector x, y, z;
0282 rota.GetComponents(x, y, z);
0283 Surface::RotationType rot(float(x.X()),
0284 float(x.Y()),
0285 float(x.Z()),
0286 float(y.X()),
0287 float(y.Y()),
0288 float(y.Z()),
0289 float(z.X()),
0290 float(z.Y()),
0291 float(z.Z()));
0292
0293 RPCRollSpecs* rollspecs = nullptr;
0294 Bounds* bounds = nullptr;
0295
0296 if (dd4hep::isA<dd4hep::Box>(fview.solid())) {
0297 const float width = dpar[0] / dd4hep::cm;
0298 const float length = dpar[1] / dd4hep::cm;
0299 const float thickness = dpar[2] / dd4hep::cm;
0300 edm::LogVerbatim("RPCGeometryBuilder")
0301 << "(3), dd4hep::Box, width: " << dpar[0] / dd4hep::cm << " length: " << dpar[1] / dd4hep::cm
0302 << " thickness: " << dpar[2] / dd4hep::cm;
0303 bounds = new RectangularPlaneBounds(width, length, thickness);
0304
0305 const std::vector<float> pars = {width, length, float(nStrips)};
0306
0307 rollspecs = new RPCRollSpecs(GeomDetEnumerators::RPCBarrel, std::string(name), pars);
0308
0309 } else {
0310 const float be = dpar[0] / dd4hep::cm;
0311 const float te = dpar[1] / dd4hep::cm;
0312 const float ap = dpar[3] / dd4hep::cm;
0313 const float ti = 0.4 / dd4hep::cm;
0314
0315 bounds = new TrapezoidalPlaneBounds(be, te, ap, ti);
0316 const std::vector<float> pars = {
0317 float(dpar[0] / dd4hep::cm), float(dpar[1] / dd4hep::cm), float(dpar[3] / dd4hep::cm), float(nStrips)};
0318 edm::LogVerbatim("RPCGeometryBuilder")
0319 << "(4), else, dpar[0] (i.e. dpar[4] for DD): " << dpar[0] / dd4hep::cm
0320 << " dpar[1] (i.e. dpar[8] for DD): " << dpar[1] / dd4hep::cm
0321 << " dpar[3] (i.e. dpar[0] for DD): " << dpar[3] / dd4hep::cm << " ti: " << ti / dd4hep::cm;
0322 rollspecs = new RPCRollSpecs(GeomDetEnumerators::RPCEndcap, std::string(name), pars);
0323
0324 Basic3DVector<float> newX(1., 0., 0.);
0325 Basic3DVector<float> newY(0., 0., 1.);
0326 newY *= -1;
0327 Basic3DVector<float> newZ(0., 1., 0.);
0328 rot.rotateAxes(newX, newY, newZ);
0329 }
0330
0331 BoundPlane* bp = new BoundPlane(pos, rot, bounds);
0332 ReferenceCountingPointer<BoundPlane> surf(bp);
0333 RPCRoll* r = new RPCRoll(rpcid, surf, rollspecs);
0334 geometry->add(r);
0335
0336 auto rls = chids.find(chid);
0337 if (rls == chids.end())
0338 rls = chids.insert(std::make_pair(chid, std::list<RPCRoll*>())).first;
0339 rls->second.emplace_back(r);
0340 }
0341
0342 for (auto& ich : chids) {
0343 const RPCDetId& chid = ich.first;
0344 const auto& rls = ich.second;
0345
0346 BoundPlane* bp = nullptr;
0347 if (!rls.empty()) {
0348 const auto& refSurf = (*rls.begin())->surface();
0349 if (chid.region() == 0) {
0350 float corners[6] = {0, 0, 0, 0, 0, 0};
0351 for (auto rl : rls) {
0352 const double h2 = rl->surface().bounds().length() / 2;
0353 const double w2 = rl->surface().bounds().width() / 2;
0354 const auto x1y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-w2, -h2, 0)));
0355 const auto x2y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+w2, +h2, 0)));
0356 corners[0] = std::min(corners[0], x1y1AtRef.x());
0357 corners[1] = std::min(corners[1], x1y1AtRef.y());
0358 corners[2] = std::max(corners[2], x2y2AtRef.x());
0359 corners[3] = std::max(corners[3], x2y2AtRef.y());
0360
0361 corners[4] = std::min(corners[4], x1y1AtRef.z());
0362 corners[5] = std::max(corners[5], x1y1AtRef.z());
0363 }
0364 const LocalPoint lpOfCentre((corners[0] + corners[2]) / 2, (corners[1] + corners[3]) / 2, 0);
0365 const auto gpOfCentre = refSurf.toGlobal(lpOfCentre);
0366 auto bounds = new RectangularPlaneBounds(
0367 (corners[2] - corners[0]) / 2, (corners[3] - corners[1]) / 2, (corners[5] - corners[4]) + 0.5);
0368 bp = new BoundPlane(gpOfCentre, refSurf.rotation(), bounds);
0369
0370 } else {
0371 float cornersLo[3] = {0, 0, 0}, cornersHi[3] = {0, 0, 0};
0372 float cornersZ[2] = {0, 0};
0373 for (auto rl : rls) {
0374 const double h2 = rl->surface().bounds().length() / 2;
0375 const double w2 = rl->surface().bounds().width() / 2;
0376 const auto& topo = dynamic_cast<const TrapezoidalStripTopology&>(rl->specificTopology());
0377 const double r = topo.radius();
0378 const double wAtLo = w2 / r * (r - h2);
0379 const double wAtHi = w2 / r * (r + h2);
0380 const auto x1y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-wAtLo, -h2, 0)));
0381 const auto x2y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+wAtLo, -h2, 0)));
0382 const auto x1y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-wAtHi, +h2, 0)));
0383 const auto x2y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+wAtHi, +h2, 0)));
0384
0385 cornersLo[0] = std::min(cornersLo[0], x1y1AtRef.x());
0386 cornersLo[1] = std::max(cornersLo[1], x2y1AtRef.x());
0387 cornersLo[2] = std::min(cornersLo[2], x1y1AtRef.y());
0388
0389 cornersHi[0] = std::min(cornersHi[0], x1y2AtRef.x());
0390 cornersHi[1] = std::max(cornersHi[1], x2y2AtRef.x());
0391 cornersHi[2] = std::max(cornersHi[2], x1y2AtRef.y());
0392
0393 cornersZ[0] = std::min(cornersZ[0], x1y1AtRef.z());
0394 cornersZ[1] = std::max(cornersZ[1], x1y1AtRef.z());
0395 }
0396 const LocalPoint lpOfCentre((cornersHi[0] + cornersHi[1]) / 2, (cornersLo[2] + cornersHi[2]) / 2, 0);
0397 const auto gpOfCentre = refSurf.toGlobal(lpOfCentre);
0398 auto bounds = new TrapezoidalPlaneBounds((cornersLo[1] - cornersLo[0]) / 2,
0399 (cornersHi[1] - cornersHi[0]) / 2,
0400 (cornersHi[2] - cornersLo[2]) / 2,
0401 (cornersZ[1] - cornersZ[0]) + 0.5);
0402 bp = new BoundPlane(gpOfCentre, refSurf.rotation(), bounds);
0403 }
0404 }
0405
0406 ReferenceCountingPointer<BoundPlane> surf(bp);
0407
0408 RPCChamber* ch = new RPCChamber(chid, surf);
0409
0410 for (auto rl : rls)
0411 ch->add(rl);
0412
0413 geometry->add(ch);
0414 }
0415 return geometry;
0416 }