CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
RPCGeometryBuilder.cc
Go to the documentation of this file.
1 /*
2 //\class RPCGeometryBuilder
3 
4  Description: RPC Geometry builder from DD & DD4hep
5  DD4hep part added to the original old file (DD version) made by M. Maggi (INFN Bari)
6 //
7 // Author: Sergio Lo Meo (sergio.lo.meo@cern.ch) following what Ianna Osborne made for DTs (DD4hep migration)
8 // Created: Fri, 20 Sep 2019
9 // Modified: Fri, 29 May 2020, following what Sunanda Banerjee made in PR #29842 PR #29943 and Ianna Osborne in PR #29954
10 */
26 #include <iostream>
27 #include <algorithm>
32 
33 using namespace cms_units::operators;
34 
36 
37 // for DDD
38 
39 std::unique_ptr<RPCGeometry> RPCGeometryBuilder::build(const DDCompactView* cview,
40  const MuonGeometryConstants& muonConstants) {
41  const std::string attribute = "ReadOutName";
42  const std::string value = "MuonRPCHits";
43  DDSpecificsMatchesValueFilter filter{DDValue(attribute, value, 0.0)};
44  DDFilteredView fview(*cview, filter);
45  return this->buildGeometry(fview, muonConstants);
46 }
47 // for DD4hep
48 
49 std::unique_ptr<RPCGeometry> RPCGeometryBuilder::build(const cms::DDCompactView* cview,
50  const MuonGeometryConstants& muonConstants) {
51  const std::string attribute = "ReadOutName";
52  const std::string value = "MuonRPCHits";
53  const cms::DDFilter filter(attribute, value);
54  cms::DDFilteredView fview(*cview, filter);
55  return this->buildGeometry(fview, muonConstants);
56 }
57 // for DDD
58 
59 std::unique_ptr<RPCGeometry> RPCGeometryBuilder::buildGeometry(DDFilteredView& fview,
60  const MuonGeometryConstants& muonConstants) {
61  LogDebug("RPCGeometryBuilder") << "Building the geometry service";
62  std::unique_ptr<RPCGeometry> geometry = std::make_unique<RPCGeometry>();
63  LogDebug("RPCGeometryBuilder") << "About to run through the RPC structure\n"
64  << " First logical part " << fview.logicalPart().name().name();
65  bool doSubDets = fview.firstChild();
66  LogDebug("RPCGeometryBuilder") << "doSubDets = " << doSubDets;
67  while (doSubDets) {
68  LogDebug("RPCGeometryBuilder") << "start the loop";
69  MuonGeometryNumbering mdddnum(muonConstants);
70  LogDebug("RPCGeometryBuilder") << "Getting the Muon base Number";
71  MuonBaseNumber mbn = mdddnum.geoHistoryToBaseNumber(fview.geoHistory());
72  LogDebug("RPCGeometryBuilder") << "Start the Rpc Numbering Schema";
73  RPCNumberingScheme rpcnum(muonConstants);
74  LogDebug("RPCGeometryBuilder") << "Getting the Unit Number";
75  const int detid = rpcnum.baseNumberToUnitNumber(mbn);
76  LogDebug("RPCGeometryBuilder") << "Getting the RPC det Id " << detid;
77  RPCDetId rpcid(detid);
78  RPCDetId chid(rpcid.region(), rpcid.ring(), rpcid.station(), rpcid.sector(), rpcid.layer(), rpcid.subsector(), 0);
79  LogDebug("RPCGeometryBuilder") << "The RPCDetid is " << rpcid;
80 
81  DDValue numbOfStrips("nStrips");
82 
83  std::vector<const DDsvalues_type*> specs(fview.specifics());
84  int nStrips = 0;
85  for (auto& spec : specs) {
86  if (DDfetch(spec, numbOfStrips)) {
87  nStrips = int(numbOfStrips.doubles()[0]);
88  }
89  }
90 
91  LogDebug("RPCGeometryBuilder") << ((nStrips == 0) ? ("No strip found!!") : (""));
92 
93  std::vector<double> dpar = fview.logicalPart().solid().parameters();
94  std::string name = fview.logicalPart().name().name();
95 
96  edm::LogVerbatim("RPCGeometryBuilder")
97  << "(1) "
98  << "detid: " << detid << " name: " << name << " number of Strips: " << nStrips;
99 
100  DDTranslation tran = fview.translation();
101  DDRotationMatrix rota = fview.rotation();
105  edm::LogVerbatim("RPCGeometryBuilder") << "(2), tran.x() " << geant_units::operators::convertMmToCm(tran.x())
106  << " tran.y(): " << geant_units::operators::convertMmToCm(tran.y())
107  << " tran.z(): " << geant_units::operators::convertMmToCm(tran.z());
108 
109  DD3Vector x, y, z;
110  rota.GetComponents(x, y, z);
111  Surface::RotationType rot(float(x.X()),
112  float(x.Y()),
113  float(x.Z()),
114  float(y.X()),
115  float(y.Y()),
116  float(y.Z()),
117  float(z.X()),
118  float(z.Y()),
119  float(z.Z()));
120 
121  RPCRollSpecs* rollspecs = nullptr;
122  Bounds* bounds = nullptr;
123 
124  if (dpar.size() == 3) {
125  const float width = geant_units::operators::convertMmToCm(dpar[0]);
126  const float length = geant_units::operators::convertMmToCm(dpar[1]);
127  const float thickness = geant_units::operators::convertMmToCm(dpar[2]);
128 
129  // Barrel
130  edm::LogVerbatim("RPCGeometryBuilder")
131  << "(3) dpar.size() == 3, width: " << width << " length: " << length << " thickness: " << thickness;
132 
133  bounds = new RectangularPlaneBounds(width, length, thickness);
134  const std::vector<float> pars = {width, length, float(numbOfStrips.doubles()[0])};
135 
136  rollspecs = new RPCRollSpecs(GeomDetEnumerators::RPCBarrel, name, pars);
137 
138  } else {
139  const float be = geant_units::operators::convertMmToCm(dpar[4]);
140  const float te = geant_units::operators::convertMmToCm(dpar[8]);
141  const float ap = geant_units::operators::convertMmToCm(dpar[0]);
142  const float ti = 0.4;
143 
144  bounds = new TrapezoidalPlaneBounds(be, te, ap, ti);
145 
146  const std::vector<float> pars = {float(geant_units::operators::convertMmToCm(dpar[4])),
149  float(numbOfStrips.doubles()[0])};
150  //Forward
151  edm::LogVerbatim("RPCGeometryBuilder")
152  << "(4), else, dpar[4]: " << be << " dpar[8]: " << te << " dpar[0]: " << ap << " ti: " << ti;
153 
154  rollspecs = new RPCRollSpecs(GeomDetEnumerators::RPCEndcap, name, pars);
155 
156  Basic3DVector<float> newX(1., 0., 0.);
157  Basic3DVector<float> newY(0., 0., 1.);
158  newY *= -1;
159  Basic3DVector<float> newZ(0., 1., 0.);
160  rot.rotateAxes(newX, newY, newZ);
161  }
162  LogDebug("RPCGeometryBuilder") << " Number of strips " << nStrips;
163 
164  BoundPlane* bp = new BoundPlane(pos, rot, bounds);
166  RPCRoll* r = new RPCRoll(rpcid, surf, rollspecs);
167  geometry->add(r);
168 
169  auto rls = chids.find(chid);
170  if (rls == chids.end())
171  rls = chids.insert(std::make_pair(chid, std::list<RPCRoll*>())).first;
172  rls->second.emplace_back(r);
173 
174  doSubDets = fview.nextSibling();
175  }
176  for (auto& ich : chids) {
177  const RPCDetId& chid = ich.first;
178  const auto& rls = ich.second;
179 
180  BoundPlane* bp = nullptr;
181  if (!rls.empty()) {
182  const auto& refSurf = (*rls.begin())->surface();
183  if (chid.region() == 0) {
184  float corners[6] = {0, 0, 0, 0, 0, 0};
185  for (auto rl : rls) {
186  const double h2 = rl->surface().bounds().length() / 2;
187  const double w2 = rl->surface().bounds().width() / 2;
188  const auto x1y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-w2, -h2, 0)));
189  const auto x2y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+w2, +h2, 0)));
190  corners[0] = std::min(corners[0], x1y1AtRef.x());
191  corners[1] = std::min(corners[1], x1y1AtRef.y());
192  corners[2] = std::max(corners[2], x2y2AtRef.x());
193  corners[3] = std::max(corners[3], x2y2AtRef.y());
194  corners[4] = std::min(corners[4], x1y1AtRef.z());
195  corners[5] = std::max(corners[5], x1y1AtRef.z());
196  }
197  const LocalPoint lpOfCentre((corners[0] + corners[2]) / 2, (corners[1] + corners[3]) / 2, 0);
198  const auto gpOfCentre = refSurf.toGlobal(lpOfCentre);
199  auto bounds = new RectangularPlaneBounds(
200  (corners[2] - corners[0]) / 2, (corners[3] - corners[1]) / 2, (corners[5] - corners[4]) + 0.5);
201  bp = new BoundPlane(gpOfCentre, refSurf.rotation(), bounds);
202  } else {
203  float cornersLo[3] = {0, 0, 0}, cornersHi[3] = {0, 0, 0};
204  float cornersZ[2] = {0, 0};
205  for (auto rl : rls) {
206  const double h2 = rl->surface().bounds().length() / 2;
207  const double w2 = rl->surface().bounds().width() / 2;
208  const auto& topo = dynamic_cast<const TrapezoidalStripTopology&>(rl->specificTopology());
209  const double r = topo.radius();
210  const double wAtLo = w2 / r * (r - h2);
211  const double wAtHi = w2 / r * (r + h2);
212 
213  const auto x1y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-wAtLo, -h2, 0)));
214  const auto x2y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+wAtLo, -h2, 0)));
215  const auto x1y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-wAtHi, +h2, 0)));
216  const auto x2y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+wAtHi, +h2, 0)));
217 
218  cornersLo[0] = std::min(cornersLo[0], x1y1AtRef.x());
219  cornersLo[1] = std::max(cornersLo[1], x2y1AtRef.x());
220  cornersLo[2] = std::min(cornersLo[2], x1y1AtRef.y());
221 
222  cornersHi[0] = std::min(cornersHi[0], x1y2AtRef.x());
223  cornersHi[1] = std::max(cornersHi[1], x2y2AtRef.x());
224  cornersHi[2] = std::max(cornersHi[2], x1y2AtRef.y());
225 
226  cornersZ[0] = std::min(cornersZ[0], x1y1AtRef.z());
227  cornersZ[1] = std::max(cornersZ[1], x1y1AtRef.z());
228  }
229  const LocalPoint lpOfCentre((cornersHi[0] + cornersHi[1]) / 2, (cornersLo[2] + cornersHi[2]) / 2, 0);
230  const auto gpOfCentre = refSurf.toGlobal(lpOfCentre);
231  auto bounds = new TrapezoidalPlaneBounds((cornersLo[1] - cornersLo[0]) / 2,
232  (cornersHi[1] - cornersHi[0]) / 2,
233  (cornersHi[2] - cornersLo[2]) / 2,
234  (cornersZ[1] - cornersZ[0]) + 0.5);
235  bp = new BoundPlane(gpOfCentre, refSurf.rotation(), bounds);
236  }
237  }
238 
240  RPCChamber* ch = new RPCChamber(chid, surf);
241  for (auto rl : rls)
242  ch->add(rl);
243  geometry->add(ch);
244  }
245 
246  return geometry;
247 }
248 
249 // for DD4hep
250 
251 std::unique_ptr<RPCGeometry> RPCGeometryBuilder::buildGeometry(cms::DDFilteredView& fview,
252  const MuonGeometryConstants& muonConstants) {
253  std::unique_ptr<RPCGeometry> geometry = std::make_unique<RPCGeometry>();
254 
255  while (fview.firstChild()) {
256  MuonGeometryNumbering mdddnum(muonConstants);
257  RPCNumberingScheme rpcnum(muonConstants);
258  int rawidCh = rpcnum.baseNumberToUnitNumber(mdddnum.geoHistoryToBaseNumber(fview.history()));
259  RPCDetId rpcid = RPCDetId(rawidCh);
260 
261  RPCDetId chid(rpcid.region(), rpcid.ring(), rpcid.station(), rpcid.sector(), rpcid.layer(), rpcid.subsector(), 0);
262 
263  auto nStrips = fview.get<double>("nStrips");
264 
265  std::vector<double> dpar = fview.parameters();
266 
267  std::string_view name = fview.name();
268 
269  edm::LogVerbatim("RPCGeometryBuilder")
270  << "(1), detid: " << rawidCh << " name: " << std::string(name) << " number of Strips: " << nStrips;
271 
272  const Double_t* tran = fview.trans();
273 
274  DDRotationMatrix rota;
275  fview.rot(rota);
276 
277  Surface::PositionType pos(tran[0] / dd4hep::cm, tran[1] / dd4hep::cm, tran[2] / dd4hep::cm);
278  edm::LogVerbatim("RPCGeometryBuilder")
279  << "(2), tran.x(): " << tran[0] / dd4hep::cm << " tran.y(): " << tran[1] / dd4hep::cm
280  << " tran.z(): " << tran[2] / dd4hep::cm;
281  DD3Vector x, y, z;
282  rota.GetComponents(x, y, z);
283  Surface::RotationType rot(float(x.X()),
284  float(x.Y()),
285  float(x.Z()),
286  float(y.X()),
287  float(y.Y()),
288  float(y.Z()),
289  float(z.X()),
290  float(z.Y()),
291  float(z.Z()));
292 
293  RPCRollSpecs* rollspecs = nullptr;
294  Bounds* bounds = nullptr;
295 
296  if (dd4hep::isA<dd4hep::Box>(fview.solid())) {
297  const float width = dpar[0] / dd4hep::cm;
298  const float length = dpar[1] / dd4hep::cm;
299  const float thickness = dpar[2] / dd4hep::cm;
300  edm::LogVerbatim("RPCGeometryBuilder")
301  << "(3), dd4hep::Box, width: " << dpar[0] / dd4hep::cm << " length: " << dpar[1] / dd4hep::cm
302  << " thickness: " << dpar[2] / dd4hep::cm;
303  bounds = new RectangularPlaneBounds(width, length, thickness);
304 
305  const std::vector<float> pars = {width, length, float(nStrips)};
306 
307  rollspecs = new RPCRollSpecs(GeomDetEnumerators::RPCBarrel, std::string(name), pars);
308 
309  } else {
310  const float be = dpar[0] / dd4hep::cm;
311  const float te = dpar[1] / dd4hep::cm;
312  const float ap = dpar[3] / dd4hep::cm;
313  const float ti = 0.4 / dd4hep::cm;
314 
315  bounds = new TrapezoidalPlaneBounds(be, te, ap, ti);
316  const std::vector<float> pars = {
317  float(dpar[0] / dd4hep::cm), float(dpar[1] / dd4hep::cm), float(dpar[3] / dd4hep::cm), float(nStrips)};
318  edm::LogVerbatim("RPCGeometryBuilder")
319  << "(4), else, dpar[0] (i.e. dpar[4] for DD): " << dpar[0] / dd4hep::cm
320  << " dpar[1] (i.e. dpar[8] for DD): " << dpar[1] / dd4hep::cm
321  << " dpar[3] (i.e. dpar[0] for DD): " << dpar[3] / dd4hep::cm << " ti: " << ti / dd4hep::cm;
322  rollspecs = new RPCRollSpecs(GeomDetEnumerators::RPCEndcap, std::string(name), pars);
323 
324  Basic3DVector<float> newX(1., 0., 0.);
325  Basic3DVector<float> newY(0., 0., 1.);
326  newY *= -1;
327  Basic3DVector<float> newZ(0., 1., 0.);
328  rot.rotateAxes(newX, newY, newZ);
329  }
330 
331  BoundPlane* bp = new BoundPlane(pos, rot, bounds);
333  RPCRoll* r = new RPCRoll(rpcid, surf, rollspecs);
334  geometry->add(r);
335 
336  auto rls = chids.find(chid);
337  if (rls == chids.end())
338  rls = chids.insert(std::make_pair(chid, std::list<RPCRoll*>())).first;
339  rls->second.emplace_back(r);
340  }
341 
342  for (auto& ich : chids) {
343  const RPCDetId& chid = ich.first;
344  const auto& rls = ich.second;
345 
346  BoundPlane* bp = nullptr;
347  if (!rls.empty()) {
348  const auto& refSurf = (*rls.begin())->surface();
349  if (chid.region() == 0) {
350  float corners[6] = {0, 0, 0, 0, 0, 0};
351  for (auto rl : rls) {
352  const double h2 = rl->surface().bounds().length() / 2;
353  const double w2 = rl->surface().bounds().width() / 2;
354  const auto x1y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-w2, -h2, 0)));
355  const auto x2y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+w2, +h2, 0)));
356  corners[0] = std::min(corners[0], x1y1AtRef.x());
357  corners[1] = std::min(corners[1], x1y1AtRef.y());
358  corners[2] = std::max(corners[2], x2y2AtRef.x());
359  corners[3] = std::max(corners[3], x2y2AtRef.y());
360 
361  corners[4] = std::min(corners[4], x1y1AtRef.z());
362  corners[5] = std::max(corners[5], x1y1AtRef.z());
363  }
364  const LocalPoint lpOfCentre((corners[0] + corners[2]) / 2, (corners[1] + corners[3]) / 2, 0);
365  const auto gpOfCentre = refSurf.toGlobal(lpOfCentre);
366  auto bounds = new RectangularPlaneBounds(
367  (corners[2] - corners[0]) / 2, (corners[3] - corners[1]) / 2, (corners[5] - corners[4]) + 0.5);
368  bp = new BoundPlane(gpOfCentre, refSurf.rotation(), bounds);
369 
370  } else {
371  float cornersLo[3] = {0, 0, 0}, cornersHi[3] = {0, 0, 0};
372  float cornersZ[2] = {0, 0};
373  for (auto rl : rls) {
374  const double h2 = rl->surface().bounds().length() / 2;
375  const double w2 = rl->surface().bounds().width() / 2;
376  const auto& topo = dynamic_cast<const TrapezoidalStripTopology&>(rl->specificTopology());
377  const double r = topo.radius();
378  const double wAtLo = w2 / r * (r - h2);
379  const double wAtHi = w2 / r * (r + h2);
380  const auto x1y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-wAtLo, -h2, 0)));
381  const auto x2y1AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+wAtLo, -h2, 0)));
382  const auto x1y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(-wAtHi, +h2, 0)));
383  const auto x2y2AtRef = refSurf.toLocal(rl->toGlobal(LocalPoint(+wAtHi, +h2, 0)));
384 
385  cornersLo[0] = std::min(cornersLo[0], x1y1AtRef.x());
386  cornersLo[1] = std::max(cornersLo[1], x2y1AtRef.x());
387  cornersLo[2] = std::min(cornersLo[2], x1y1AtRef.y());
388 
389  cornersHi[0] = std::min(cornersHi[0], x1y2AtRef.x());
390  cornersHi[1] = std::max(cornersHi[1], x2y2AtRef.x());
391  cornersHi[2] = std::max(cornersHi[2], x1y2AtRef.y());
392 
393  cornersZ[0] = std::min(cornersZ[0], x1y1AtRef.z());
394  cornersZ[1] = std::max(cornersZ[1], x1y1AtRef.z());
395  }
396  const LocalPoint lpOfCentre((cornersHi[0] + cornersHi[1]) / 2, (cornersLo[2] + cornersHi[2]) / 2, 0);
397  const auto gpOfCentre = refSurf.toGlobal(lpOfCentre);
398  auto bounds = new TrapezoidalPlaneBounds((cornersLo[1] - cornersLo[0]) / 2,
399  (cornersHi[1] - cornersHi[0]) / 2,
400  (cornersHi[2] - cornersLo[2]) / 2,
401  (cornersZ[1] - cornersZ[0]) + 0.5);
402  bp = new BoundPlane(gpOfCentre, refSurf.rotation(), bounds);
403  }
404  }
405 
407 
408  RPCChamber* ch = new RPCChamber(chid, surf);
409 
410  for (auto rl : rls)
411  ch->add(rl);
412 
413  geometry->add(ch);
414  }
415  return geometry;
416 }
Log< level::Info, true > LogVerbatim
void add(RPCRoll *rl)
Add Roll to the chamber which takes ownership.
Definition: RPCChamber.cc:32
const std::vector< double > & parameters(void) const
Give the parameters of the solid.
Definition: DDSolid.cc:125
MuonBaseNumber geoHistoryToBaseNumber(const DDGeoHistory &history) const
const DDLogicalPart & logicalPart() const
The logical-part of the current node in the filtered-view.
const std::vector< double > & doubles() const
a reference to the double-valued values stored in the given instance of DDValue
Definition: DDValue.cc:111
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double >> DD3Vector
const N & name() const
Definition: DDBase.h:59
virtual float length() const =0
common ppss p3p6s2 common epss epspn46 common const1 w2
Definition: inclppp.h:1
bool nextSibling()
set the current node to the next sibling ...
int baseNumberToUnitNumber(const MuonBaseNumber &) const override
const DDRotationMatrix & rotation() const
The absolute rotation of the current node.
Plane BoundPlane
Definition: Plane.h:94
const Bounds & bounds() const
Definition: Surface.h:87
const DDSolid & solid(void) const
Returns a reference object of the solid being the shape of this LogicalPart.
Compact representation of the geometrical detector hierarchy.
Definition: DDCompactView.h:81
const DDGeoHistory & geoHistory() const
The list of ancestors up to the root-node of the current node.
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
Definition: DDsvalues.cc:79
tuple nStrips
1.2 is to make the matching window safely the two nearest strips 0.35 is the size of an ME0 chamber i...
TkRotation & rotateAxes(const Basic3DVector< T > &newX, const Basic3DVector< T > &newY, const Basic3DVector< T > &newZ)
std::unique_ptr< RPCGeometry > buildGeometry(DDFilteredView &fview, const MuonGeometryConstants &muonConstants)
ROOT::Math::Rotation3D DDRotationMatrix
A DDRotationMatrix is currently implemented with a ROOT Rotation3D.
const ExpandedNodes & history()
The numbering history of the current node.
int ring() const
Definition: RPCDetId.h:59
const Double_t * rot() const
The absolute rotation of the current node.
std::string_view name() const
bool firstChild()
set the current node to the first child
int layer() const
Definition: RPCDetId.h:85
T get(const std::string &)
extract attribute value
constexpr NumType convertMmToCm(NumType millimeters)
Definition: angle_units.h:44
const Double_t * trans() const
The absolute translation of the current node.
std::unique_ptr< RPCGeometry > build(const DDCompactView *cview, const MuonGeometryConstants &muonConstants)
int sector() const
Sector id: the group of chambers at same phi (and increasing r)
Definition: RPCDetId.h:81
int subsector() const
SubSector id : some sectors are divided along the phi direction in subsectors (from 1 to 4 in Barrel...
Definition: RPCDetId.h:88
bool firstChild()
set the current node to the first child ...
Definition: Bounds.h:18
const std::vector< double > parameters() const
extract shape parameters
const DDTranslation & translation() const
The absolute translation of the current node.
std::vector< const DDsvalues_type * > specifics() const
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double > > DDTranslation
Definition: DDTranslation.h:7
const std::string & name() const
Returns the name.
Definition: DDName.cc:41
dd4hep::Solid solid() const
#define LogDebug(id)
int region() const
Region id: 0 for Barrel, +/-1 For +/- Endcap.
Definition: RPCDetId.h:53
int station() const
Definition: RPCDetId.h:78