File indexing completed on 2024-04-06 12:15:20
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 edm::LogVerbatim("RPCGeometryBuilder") << "Building the geometry service";
0062 std::unique_ptr<RPCGeometry> geometry = std::make_unique<RPCGeometry>();
0063 edm::LogVerbatim("RPCGeometryBuilder") << "About to run through the RPC structure\n"
0064 << " First logical part " << fview.logicalPart().name().name();
0065 bool doSubDets = fview.firstChild();
0066 edm::LogVerbatim("RPCGeometryBuilder") << "doSubDets = " << doSubDets;
0067 while (doSubDets) {
0068 edm::LogVerbatim("RPCGeometryBuilder") << "start the loop";
0069 MuonGeometryNumbering mdddnum(muonConstants);
0070 edm::LogVerbatim("RPCGeometryBuilder") << "Getting the Muon base Number";
0071 MuonBaseNumber mbn = mdddnum.geoHistoryToBaseNumber(fview.geoHistory());
0072
0073 RPCNumberingScheme rpcnum(muonConstants);
0074 edm::LogVerbatim("RPCGeometryBuilder") << "Getting the Unit Number";
0075 const int detid = rpcnum.baseNumberToUnitNumber(mbn);
0076 edm::LogVerbatim("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 edm::LogVerbatim("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 if (nStrips == 0)
0092 edm::LogVerbatim("RPCGeometryBuilder") << "No strip found!!";
0093
0094 std::vector<double> dpar = fview.logicalPart().solid().parameters();
0095 std::string name = fview.logicalPart().name().name();
0096
0097 edm::LogVerbatim("RPCGeometryBuilder")
0098 << "(1) "
0099 << "detid: " << detid << " name: " << name << " number of Strips: " << nStrips;
0100
0101 DDTranslation tran = fview.translation();
0102 DDRotationMatrix rota = fview.rotation();
0103 Surface::PositionType pos(geant_units::operators::convertMmToCm(tran.x()),
0104 geant_units::operators::convertMmToCm(tran.y()),
0105 geant_units::operators::convertMmToCm(tran.z()));
0106 edm::LogVerbatim("RPCGeometryBuilder") << "(2) tran.x(): " << geant_units::operators::convertMmToCm(tran.x())
0107 << " tran.y(): " << geant_units::operators::convertMmToCm(tran.y())
0108 << " tran.z(): " << geant_units::operators::convertMmToCm(tran.z());
0109
0110 DD3Vector x, y, z;
0111 rota.GetComponents(x, y, z);
0112 Surface::RotationType rot(float(x.X()),
0113 float(x.Y()),
0114 float(x.Z()),
0115 float(y.X()),
0116 float(y.Y()),
0117 float(y.Z()),
0118 float(z.X()),
0119 float(z.Y()),
0120 float(z.Z()));
0121
0122 RPCRollSpecs* rollspecs = nullptr;
0123 Bounds* bounds = nullptr;
0124
0125 if (dpar.size() == 3) {
0126 const float width = geant_units::operators::convertMmToCm(dpar[0]);
0127 const float length = geant_units::operators::convertMmToCm(dpar[1]);
0128 const float thickness = geant_units::operators::convertMmToCm(dpar[2]);
0129
0130
0131 edm::LogVerbatim("RPCGeometryBuilder")
0132 << "(3) Box, width: " << width << " length: " << length << " thickness: " << thickness;
0133
0134 bounds = new RectangularPlaneBounds(width, length, thickness);
0135 const std::vector<float> pars = {width, length, float(numbOfStrips.doubles()[0])};
0136
0137 rollspecs = new RPCRollSpecs(GeomDetEnumerators::RPCBarrel, name, pars);
0138
0139 } else {
0140 const float be = geant_units::operators::convertMmToCm(dpar[4]);
0141 const float te = geant_units::operators::convertMmToCm(dpar[8]);
0142 const float ap = geant_units::operators::convertMmToCm(dpar[0]);
0143 const float ti = 0.4;
0144
0145 bounds = new TrapezoidalPlaneBounds(be, te, ap, ti);
0146
0147 const std::vector<float> pars = {be, te, ap, float(numbOfStrips.doubles()[0])};
0148
0149 edm::LogVerbatim("RPCGeometryBuilder") << "(4) be: " << be << " te: " << te << " ap: " << ap << " ti: " << ti
0150 << " strips " << numbOfStrips.doubles()[0];
0151
0152 rollspecs = new RPCRollSpecs(GeomDetEnumerators::RPCEndcap, name, pars);
0153
0154 Basic3DVector<float> newX(1., 0., 0.);
0155 Basic3DVector<float> newY(0., 0., 1.);
0156 newY *= -1;
0157 Basic3DVector<float> newZ(0., 1., 0.);
0158 rot.rotateAxes(newX, newY, newZ);
0159 }
0160 edm::LogVerbatim("RPCGeometryBuilder") << " Number of strips " << nStrips;
0161
0162 BoundPlane* bp = new BoundPlane(pos, rot, bounds);
0163 ReferenceCountingPointer<BoundPlane> surf(bp);
0164 RPCRoll* r = new RPCRoll(rpcid, surf, rollspecs);
0165 geometry->add(r);
0166
0167 auto rls = chids.find(chid);
0168 if (rls == chids.end())
0169 rls = chids.insert(std::make_pair(chid, std::list<RPCRoll*>())).first;
0170 rls->second.emplace_back(r);
0171
0172 doSubDets = fview.nextSibling();
0173 }
0174 for (auto& ich : chids) {
0175 const RPCDetId& chid = ich.first;
0176 const auto& rls = ich.second;
0177
0178 BoundPlane* bp = nullptr;
0179 if (!rls.empty()) {
0180 const auto& refSurf = (*rls.begin())->surface();
0181 if (chid.region() == 0) {
0182 float corners[6] = {0, 0, 0, 0, 0, 0};
0183 for (auto rl : rls) {
0184 const double h2 = rl->surface().bounds().length() / 2;
0185 const double w2 = rl->surface().bounds().width() / 2;
0186 const auto x1y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-w2, -h2, 0)));
0187 const auto x2y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+w2, +h2, 0)));
0188 corners[0] = std::min(corners[0], x1y1AtRef.x());
0189 corners[1] = std::min(corners[1], x1y1AtRef.y());
0190 corners[2] = std::max(corners[2], x2y2AtRef.x());
0191 corners[3] = std::max(corners[3], x2y2AtRef.y());
0192 corners[4] = std::min(corners[4], x1y1AtRef.z());
0193 corners[5] = std::max(corners[5], x1y1AtRef.z());
0194 }
0195 const LocalPoint lpOfCentre((corners[0] + corners[2]) / 2, (corners[1] + corners[3]) / 2, 0);
0196 const auto gpOfCentre = refSurf.toGlobal(lpOfCentre);
0197 auto bounds = new RectangularPlaneBounds(
0198 (corners[2] - corners[0]) / 2, (corners[3] - corners[1]) / 2, (corners[5] - corners[4]) + 0.5);
0199 bp = new BoundPlane(gpOfCentre, refSurf.rotation(), bounds);
0200 } else {
0201 float cornersLo[3] = {0, 0, 0}, cornersHi[3] = {0, 0, 0};
0202 float cornersZ[2] = {0, 0};
0203 for (auto rl : rls) {
0204 const double h2 = rl->surface().bounds().length() / 2;
0205 const double w2 = rl->surface().bounds().width() / 2;
0206 const auto& topo = dynamic_cast<const TrapezoidalStripTopology&>(rl->specificTopology());
0207 const double r = topo.radius();
0208 const double wAtLo = w2 / r * (r - h2);
0209 const double wAtHi = w2 / r * (r + h2);
0210
0211 const auto x1y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-wAtLo, -h2, 0)));
0212 const auto x2y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+wAtLo, -h2, 0)));
0213 const auto x1y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-wAtHi, +h2, 0)));
0214 const auto x2y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+wAtHi, +h2, 0)));
0215
0216 cornersLo[0] = std::min(cornersLo[0], x1y1AtRef.x());
0217 cornersLo[1] = std::max(cornersLo[1], x2y1AtRef.x());
0218 cornersLo[2] = std::min(cornersLo[2], x1y1AtRef.y());
0219
0220 cornersHi[0] = std::min(cornersHi[0], x1y2AtRef.x());
0221 cornersHi[1] = std::max(cornersHi[1], x2y2AtRef.x());
0222 cornersHi[2] = std::max(cornersHi[2], x1y2AtRef.y());
0223
0224 cornersZ[0] = std::min(cornersZ[0], x1y1AtRef.z());
0225 cornersZ[1] = std::max(cornersZ[1], x1y1AtRef.z());
0226 }
0227 const LocalPoint lpOfCentre((cornersHi[0] + cornersHi[1]) / 2, (cornersLo[2] + cornersHi[2]) / 2, 0);
0228 const auto gpOfCentre = refSurf.toGlobal(lpOfCentre);
0229 auto bounds = new TrapezoidalPlaneBounds((cornersLo[1] - cornersLo[0]) / 2,
0230 (cornersHi[1] - cornersHi[0]) / 2,
0231 (cornersHi[2] - cornersLo[2]) / 2,
0232 (cornersZ[1] - cornersZ[0]) + 0.5);
0233 bp = new BoundPlane(gpOfCentre, refSurf.rotation(), bounds);
0234 }
0235 }
0236
0237 ReferenceCountingPointer<BoundPlane> surf(bp);
0238 RPCChamber* ch = new RPCChamber(chid, surf);
0239 for (auto rl : rls)
0240 ch->add(rl);
0241 geometry->add(ch);
0242 }
0243
0244 return geometry;
0245 }
0246
0247
0248
0249 std::unique_ptr<RPCGeometry> RPCGeometryBuilder::buildGeometry(cms::DDFilteredView& fview,
0250 const MuonGeometryConstants& muonConstants) {
0251 edm::LogVerbatim("RPCGeometryBuilder") << "Building the geometry service";
0252 std::unique_ptr<RPCGeometry> geometry = std::make_unique<RPCGeometry>();
0253
0254 while (fview.firstChild()) {
0255 edm::LogVerbatim("RPCGeometryBuilder") << "start the loop";
0256 MuonGeometryNumbering mdddnum(muonConstants);
0257 edm::LogVerbatim("RPCGeometryBuilder") << "Getting the Muon base Number";
0258 RPCNumberingScheme rpcnum(muonConstants);
0259 edm::LogVerbatim("RPCGeometryBuilder") << "Getting the Unit Number";
0260 int rawidCh = rpcnum.baseNumberToUnitNumber(mdddnum.geoHistoryToBaseNumber(fview.history()));
0261 edm::LogVerbatim("RPCGeometryBuilder") << "Getting the RPC det Id " << rawidCh;
0262 RPCDetId rpcid = RPCDetId(rawidCh);
0263 edm::LogVerbatim("RPCGeometryBuilder") << "The RPCDetid is " << rpcid;
0264
0265 RPCDetId chid(rpcid.region(), rpcid.ring(), rpcid.station(), rpcid.sector(), rpcid.layer(), rpcid.subsector(), 0);
0266
0267 auto nStrips = fview.get<double>("nStrips");
0268
0269 std::vector<double> dpar = fview.parameters();
0270
0271 std::string_view name = fview.name();
0272
0273 edm::LogVerbatim("RPCGeometryBuilder")
0274 << "(1) detid: " << rawidCh << " name: " << std::string(name) << " number of Strips: " << nStrips;
0275
0276 const Double_t* tran = fview.trans();
0277
0278 DDRotationMatrix rota;
0279 fview.rot(rota);
0280
0281 Surface::PositionType pos(tran[0] / dd4hep::cm, tran[1] / dd4hep::cm, tran[2] / dd4hep::cm);
0282 edm::LogVerbatim("RPCGeometryBuilder")
0283 << "(2) tran.x(): " << tran[0] / dd4hep::cm << " tran.y(): " << tran[1] / dd4hep::cm
0284 << " tran.z(): " << tran[2] / dd4hep::cm;
0285 DD3Vector x, y, z;
0286 rota.GetComponents(x, y, z);
0287 Surface::RotationType rot(float(x.X()),
0288 float(x.Y()),
0289 float(x.Z()),
0290 float(y.X()),
0291 float(y.Y()),
0292 float(y.Z()),
0293 float(z.X()),
0294 float(z.Y()),
0295 float(z.Z()));
0296
0297 RPCRollSpecs* rollspecs = nullptr;
0298 Bounds* bounds = nullptr;
0299
0300 if (dd4hep::isA<dd4hep::Box>(fview.solid())) {
0301 const float width = dpar[0] / dd4hep::cm;
0302 const float length = dpar[1] / dd4hep::cm;
0303 const float thickness = dpar[2] / dd4hep::cm;
0304 edm::LogVerbatim("RPCGeometryBuilder")
0305 << "(3) Box, width: " << dpar[0] / dd4hep::cm << " length: " << dpar[1] / dd4hep::cm
0306 << " thickness: " << dpar[2] / dd4hep::cm;
0307 bounds = new RectangularPlaneBounds(width, length, thickness);
0308
0309 const std::vector<float> pars = {width, length, float(nStrips)};
0310
0311 rollspecs = new RPCRollSpecs(GeomDetEnumerators::RPCBarrel, std::string(name), pars);
0312
0313 } else {
0314 const float be = dpar[0] / dd4hep::cm;
0315 const float te = dpar[1] / dd4hep::cm;
0316 const float ap = dpar[3] / dd4hep::cm;
0317 const float ti = 0.4;
0318
0319 bounds = new TrapezoidalPlaneBounds(be, te, ap, ti);
0320 const std::vector<float> pars = {be, te, ap, float(nStrips)};
0321 edm::LogVerbatim("RPCGeometryBuilder")
0322 << "(4) be: " << be << " te: " << te << " ap: " << ap << " ti: " << ti << " strips " << nStrips;
0323 rollspecs = new RPCRollSpecs(GeomDetEnumerators::RPCEndcap, std::string(name), pars);
0324
0325 Basic3DVector<float> newX(1., 0., 0.);
0326 Basic3DVector<float> newY(0., 0., 1.);
0327 newY *= -1;
0328 Basic3DVector<float> newZ(0., 1., 0.);
0329 rot.rotateAxes(newX, newY, newZ);
0330 }
0331 edm::LogVerbatim("RPCGeometryBuilder") << " Number of strips " << nStrips;
0332
0333 BoundPlane* bp = new BoundPlane(pos, rot, bounds);
0334 ReferenceCountingPointer<BoundPlane> surf(bp);
0335 RPCRoll* r = new RPCRoll(rpcid, surf, rollspecs);
0336 geometry->add(r);
0337
0338 auto rls = chids.find(chid);
0339 if (rls == chids.end())
0340 rls = chids.insert(std::make_pair(chid, std::list<RPCRoll*>())).first;
0341 rls->second.emplace_back(r);
0342 }
0343
0344 for (auto& ich : chids) {
0345 const RPCDetId& chid = ich.first;
0346 const auto& rls = ich.second;
0347
0348 BoundPlane* bp = nullptr;
0349 if (!rls.empty()) {
0350 const auto& refSurf = (*rls.begin())->surface();
0351 if (chid.region() == 0) {
0352 float corners[6] = {0, 0, 0, 0, 0, 0};
0353 for (auto rl : rls) {
0354 const double h2 = rl->surface().bounds().length() / 2;
0355 const double w2 = rl->surface().bounds().width() / 2;
0356 const auto x1y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-w2, -h2, 0)));
0357 const auto x2y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+w2, +h2, 0)));
0358 corners[0] = std::min(corners[0], x1y1AtRef.x());
0359 corners[1] = std::min(corners[1], x1y1AtRef.y());
0360 corners[2] = std::max(corners[2], x2y2AtRef.x());
0361 corners[3] = std::max(corners[3], x2y2AtRef.y());
0362
0363 corners[4] = std::min(corners[4], x1y1AtRef.z());
0364 corners[5] = std::max(corners[5], x1y1AtRef.z());
0365 }
0366 const LocalPoint lpOfCentre((corners[0] + corners[2]) / 2, (corners[1] + corners[3]) / 2, 0);
0367 const auto gpOfCentre = refSurf.toGlobal(lpOfCentre);
0368 auto bounds = new RectangularPlaneBounds(
0369 (corners[2] - corners[0]) / 2, (corners[3] - corners[1]) / 2, (corners[5] - corners[4]) + 0.5);
0370 bp = new BoundPlane(gpOfCentre, refSurf.rotation(), bounds);
0371
0372 } else {
0373 float cornersLo[3] = {0, 0, 0}, cornersHi[3] = {0, 0, 0};
0374 float cornersZ[2] = {0, 0};
0375 for (auto rl : rls) {
0376 const double h2 = rl->surface().bounds().length() / 2;
0377 const double w2 = rl->surface().bounds().width() / 2;
0378 const auto& topo = dynamic_cast<const TrapezoidalStripTopology&>(rl->specificTopology());
0379 const double r = topo.radius();
0380 const double wAtLo = w2 / r * (r - h2);
0381 const double wAtHi = w2 / r * (r + h2);
0382 const auto x1y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-wAtLo, -h2, 0)));
0383 const auto x2y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+wAtLo, -h2, 0)));
0384 const auto x1y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-wAtHi, +h2, 0)));
0385 const auto x2y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+wAtHi, +h2, 0)));
0386
0387 cornersLo[0] = std::min(cornersLo[0], x1y1AtRef.x());
0388 cornersLo[1] = std::max(cornersLo[1], x2y1AtRef.x());
0389 cornersLo[2] = std::min(cornersLo[2], x1y1AtRef.y());
0390
0391 cornersHi[0] = std::min(cornersHi[0], x1y2AtRef.x());
0392 cornersHi[1] = std::max(cornersHi[1], x2y2AtRef.x());
0393 cornersHi[2] = std::max(cornersHi[2], x1y2AtRef.y());
0394
0395 cornersZ[0] = std::min(cornersZ[0], x1y1AtRef.z());
0396 cornersZ[1] = std::max(cornersZ[1], x1y1AtRef.z());
0397 }
0398 const LocalPoint lpOfCentre((cornersHi[0] + cornersHi[1]) / 2, (cornersLo[2] + cornersHi[2]) / 2, 0);
0399 const auto gpOfCentre = refSurf.toGlobal(lpOfCentre);
0400 auto bounds = new TrapezoidalPlaneBounds((cornersLo[1] - cornersLo[0]) / 2,
0401 (cornersHi[1] - cornersHi[0]) / 2,
0402 (cornersHi[2] - cornersLo[2]) / 2,
0403 (cornersZ[1] - cornersZ[0]) + 0.5);
0404 bp = new BoundPlane(gpOfCentre, refSurf.rotation(), bounds);
0405 }
0406 }
0407
0408 ReferenceCountingPointer<BoundPlane> surf(bp);
0409
0410 RPCChamber* ch = new RPCChamber(chid, surf);
0411
0412 for (auto rl : rls)
0413 ch->add(rl);
0414
0415 geometry->add(ch);
0416 }
0417 return geometry;
0418 }
