Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-11-09 23:32:01 

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   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       // Barrel
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       //Forward
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 // for DD4hep
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!