Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 09:50:04

0001 /*
0002 //\class RPCGeometryBuilder
0003 
0004  Description: RPC Geometry builder from DD & DD4hep
0005               DD4hep part added to the original old file (DD version) made by M. Maggi (INFN Bari)
0006 //
0007 // Author:  Sergio Lo Meo (sergio.lo.meo@cern.ch) following what Ianna Osborne made for DTs (DD4hep migration)
0008 //          Created:  Fri, 20 Sep 2019 
0009 //          Modified: Fri, 29 May 2020, following what Sunanda Banerjee made in PR #29842 PR #29943 and Ianna Osborne in PR #29954    
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 // for DDD
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 // for DD4hep
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 // for DDD
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     //  edm::LogVerbatim("RPCGeometryBuilder") << "Start the Rpc Numbering Schema";
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       // Barrel
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       //Forward
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 // for DD4hep
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 }