ATLAS Offline Software
GMTreeBrowser.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // GMTreeBrowser.cxx, (c) ATLAS Detector software
8 
10 // GeoModel includes
11 #include <cmath>
12 #include <iomanip>
13 #include <iostream>
14 
15 #include "GeoModelKernel/GeoBox.h"
16 #include "GeoModelKernel/GeoCons.h"
17 #include "GeoModelKernel/GeoPara.h"
18 #include "GeoModelKernel/GeoPcon.h"
19 #include "GeoModelKernel/GeoPgon.h"
20 #include "GeoModelKernel/GeoShape.h"
21 #include "GeoModelKernel/GeoShapeIntersection.h"
22 #include "GeoModelKernel/GeoShapeShift.h"
23 #include "GeoModelKernel/GeoShapeSubtraction.h"
24 #include "GeoModelKernel/GeoShapeUnion.h"
25 #include "GeoModelKernel/GeoSimplePolygonBrep.h"
26 #include "GeoModelKernel/GeoTrap.h"
27 #include "GeoModelKernel/GeoTrd.h"
28 #include "GeoModelKernel/GeoTube.h"
29 #include "GeoModelKernel/GeoTubs.h"
30 #include "GeoModelKernel/GeoVPhysVol.h"
32 #include "CxxUtils/inline_hints.h"
33 
34 #if defined(FLATTEN)
35  // We compile this package with optimization, even in debug builds; otherwise,
36  // the heavy use of Eigen makes it too slow. However, from here we may call
37  // to out-of-line Eigen code that is linked from other DSOs; in that case,
38  // it would not be optimized. Avoid this by forcing all Eigen code
39  // to be inlined here if possible.
41 #endif
42 int Trk::GMTreeBrowser::compareGeoVolumes(const GeoVPhysVol* gv1,
43  const GeoVPhysVol* gv2,
44  double tolerance, bool dumpInfo,
45  int level) const {
46 
47  int diff = 0;
48 
49  // CASE 1: naming difference
50  if (gv1->getLogVol()->getName() != gv2->getLogVol()->getName()) {
51  diff = 1000 * level + 1;
52  if (dumpInfo) {
53  std::cout << "CASE 1: names differ at level:" << level << ":"
54  << gv1->getLogVol()->getName() << ":"
55  << gv2->getLogVol()->getName() << std::endl;
56  }
57  // else return m_diff; // the naming change is harmless and can mask
58  // more serious problems, continue the loop
59  }
60  // CASE 2: material type difference
61  if (gv1->getLogVol()->getMaterial()->getName() !=
62  gv2->getLogVol()->getMaterial()->getName()) {
63  diff = 1000 * level + 2;
64  if (dumpInfo) {
65  std::cout << "CASE 2: material types differ for volume:"
66  << gv1->getLogVol()->getName() << ":at level:" << level
67  << ":" << gv1->getLogVol()->getMaterial()->getName()
68  << ":differs from:"
69  << gv2->getLogVol()->getMaterial()->getName()
70  << std::endl;
71  } else
72  return diff;
73  }
74  // CASE 3: shape type difference
75  if (gv1->getLogVol()->getShape()->typeID() !=
76  gv2->getLogVol()->getShape()->typeID()) {
77  diff = 1000 * level + 3;
78  if (dumpInfo) {
79  std::cout << "CASE 3: shape types differ at level:" << level
80  << ":volume name:" << gv1->getLogVol()->getName()
81  << ":shape:" << gv1->getLogVol()->getShape()->type()
82  << ":shape ref:" << gv2->getLogVol()->getShape()->type()
83  << std::endl;
84  } else
85  return diff;
86  }
87  // CASE 4: difference in shape definition
88  if (!compareShapes(gv1->getLogVol()->getShape(),
89  gv2->getLogVol()->getShape(), tolerance)) {
90  diff = 1000 * level + 4;
91  if (dumpInfo) {
92  std::cout << "CASE 4: shape definition differ at level:" << level
93  << std::endl;
94  } else
95  return diff;
96  }
97  unsigned int nChild1 = gv1->getNChildVols();
98  // CASE 5: difference in the number of child volumes
99  if (nChild1 != gv2->getNChildVols()) {
100  diff = 1000 * level + 5;
101  if (dumpInfo) {
102  std::cout << "CASE 5: number of child vols differ at level:"
103  << level << ":volume name:" << gv1->getLogVol()->getName()
104  << ":nChildVols:" << gv1->getNChildVols()
105  << ":nChildVols ref:" << gv2->getNChildVols()
106  << std::endl;
107  } else
108  return diff;
109  }
110 
111  // CASE 6 & 7: transform to child difference
112  // We used to do this with something like
113  // for (unsigned int ic = 0; ic < gv1->getNChildVols(); ic++) {
114  // GeoTrf::Transform3D transf1 = gv1->getXToChildVol(ic);
115  // const GeoVPhysVol* cv1 = &(*(gv1->getChildVol(ic)));
116  //
117  // But getXToChildVol and getChildVol need to walk all the children
118  // until they reach the given index. So this would be N^2,
119  // and each time we repeat it for the transform and the volume.
120  // Better to use geoGetVolumes so that we only need do the walk once.
121  // And examination of profiling data shows that almost never
122  // fail a comparison from here on, so we'll almost always be examining
123  // all children anyway.
124  // (It would be even better if GeoVPhysVol has some sort of iterator
125  // interface. Maybe we can use a generator with C++23...)
126  GeoVolumeVec_t children1 = geoGetVolumes (gv1, 1, nChild1);
127  GeoVolumeVec_t children2 = geoGetVolumes (gv2, 1, nChild1);
128  assert (children1.size() == nChild1 && children2.size() == nChild1);
129  for (unsigned int ic = 0; ic < nChild1; ic++) {
130  GeoTrf::Transform3D& transf1 = children1.at(ic).second;
131  GeoTrf::Transform3D& transf2 = children2.at(ic).second;
132 
133  const GeoVPhysVol* cv1 = children1.at(ic).first;
134  const GeoVPhysVol* cv2 = children2.at(ic).first;
135 
136  if ((transf1.translation() - transf2.translation()).norm() >
137  tolerance) {
138  diff = 1000 * level + 10 * ic + 6;
139  if (dumpInfo) {
140  std::cout << "CASE 6: translation differs at level:" << level
141  << ": between mother and child:"
142  << gv1->getLogVol()->getName() << ":"
143  << cv1->getLogVol()->getName() << std::endl;
145  } else
146  return diff;
147  }
148  // For rotation matrices, transpose is the same as inverse.
150  transf1.rotation() * transf2.rotation().transpose();
151  if (std::abs(rot(0, 1)) > tolerance ||
152  std::abs(rot(0, 2)) > tolerance ||
153  std::abs(rot(1, 2)) > tolerance) {
154  diff = 1000 * level + 10 * ic + 7;
155  if (dumpInfo) {
156  std::cout << "CASE 7: rotation differs at level:" << level
157  << ":between mother and child:"
158  << gv1->getLogVol()->getName() << ":"
159  << cv1->getLogVol()->getName() << std::endl;
161  } else
162  return diff;
163  }
164 
165  int child_comp =
166  compareGeoVolumes(cv1, cv2, tolerance, dumpInfo, level + 1);
167  if (child_comp != 0)
168  return child_comp;
169  }
170 
171  return diff;
172 }
173 
174 #if defined(FLATTEN)
175  // We compile this package with optimization, even in debug builds; otherwise,
176  // the heavy use of Eigen makes it too slow. However, from here we may call
177  // to out-of-line Eigen code that is linked from other DSOs; in that case,
178  // it would not be optimized. Avoid this by forcing all Eigen code
179  // to be inlined here if possible.
181 #endif
182 bool Trk::GMTreeBrowser::compareShapes(const GeoShape* sh1, const GeoShape* sh2,
183  double tol) const {
184 
185  if (sh1->typeID() != sh2->typeID())
186  return false;
187 
188  if (sh1->type() == "Pgon") {
189  const GeoPgon* pgon1 = static_cast<const GeoPgon*>(sh1);
190  const GeoPgon* pgon2 = static_cast<const GeoPgon*>(sh2);
191  if (!pgon1 || !pgon2)
192  return false;
193 
194  if (pgon1->getNPlanes() != pgon2->getNPlanes())
195  return false;
196  if (pgon1->getNSides() != pgon2->getNSides())
197  return false;
198  if (std::abs(pgon1->getSPhi() - pgon2->getSPhi()) > tol)
199  return false;
200  if (std::abs(pgon1->getDPhi() - pgon2->getDPhi()) > tol)
201  return false;
202 
203  return true;
204 
205  } else if (sh1->type() == "Trd") {
206  const GeoTrd* trd1 = static_cast<const GeoTrd*>(sh1);
207  const GeoTrd* trd2 = static_cast<const GeoTrd*>(sh2);
208 
209  if (std::abs(trd1->getXHalfLength1() - trd2->getXHalfLength1()) > tol)
210  return false;
211  if (std::abs(trd1->getXHalfLength2() - trd2->getXHalfLength2()) > tol)
212  return false;
213  if (std::abs(trd1->getYHalfLength1() - trd2->getYHalfLength1()) > tol)
214  return false;
215  if (std::abs(trd1->getYHalfLength2() - trd2->getYHalfLength2()) > tol)
216  return false;
217  if (std::abs(trd1->getZHalfLength() - trd2->getZHalfLength()) > tol)
218  return false;
219 
220  return true;
221 
222  } else if (sh1->type() == "Box") {
223  const GeoBox* box1 = static_cast<const GeoBox*>(sh1);
224  const GeoBox* box2 = static_cast<const GeoBox*>(sh2);
225 
226  if (std::abs(box1->getXHalfLength() - box2->getXHalfLength()) > tol)
227  return false;
228  if (std::abs(box1->getYHalfLength() - box2->getYHalfLength()) > tol)
229  return false;
230  if (std::abs(box1->getZHalfLength() - box2->getZHalfLength()) > tol)
231  return false;
232 
233  return true;
234 
235  } else if (sh1->type() == "Tube") {
236  const GeoTube* tube1 = static_cast<const GeoTube*>(sh1);
237  const GeoTube* tube2 = static_cast<const GeoTube*>(sh2);
238 
239  if (std::abs(tube1->getRMin() - tube2->getRMin()) > tol)
240  return false;
241  if (std::abs(tube1->getRMax() - tube2->getRMax()) > tol)
242  return false;
243  if (std::abs(tube1->getZHalfLength() - tube2->getZHalfLength()) > tol)
244  return false;
245 
246  return true;
247 
248  } else if (sh1->type() == "Tubs") {
249  const GeoTubs* tubs1 = static_cast<const GeoTubs*>(sh1);
250  const GeoTubs* tubs2 = static_cast<const GeoTubs*>(sh2);
251 
252  if (std::abs(tubs1->getRMin() - tubs2->getRMin()) > tol)
253  return false;
254  if (std::abs(tubs1->getRMax() - tubs2->getRMax()) > tol)
255  return false;
256  if (std::abs(tubs1->getZHalfLength() - tubs2->getZHalfLength()) > tol)
257  return false;
258  if (std::abs(tubs1->getSPhi() - tubs2->getSPhi()) > tol)
259  return false;
260  if (std::abs(tubs1->getDPhi() - tubs2->getDPhi()) > tol)
261  return false;
262 
263  return true;
264 
265  } else if (sh1->type() == "Cons") {
266  const GeoCons* cons1 = static_cast<const GeoCons*>(sh1);
267  const GeoCons* cons2 = static_cast<const GeoCons*>(sh2);
268 
269  if (std::abs(cons1->getRMin1() - cons2->getRMin1()) > tol)
270  return false;
271  if (std::abs(cons1->getRMin2() - cons2->getRMin2()) > tol)
272  return false;
273  if (std::abs(cons1->getRMax1() - cons2->getRMax1()) > tol)
274  return false;
275  if (std::abs(cons1->getRMax2() - cons2->getRMax2()) > tol)
276  return false;
277  if (std::abs(cons1->getDZ() - cons2->getDZ()) > tol)
278  return false;
279  if (std::abs(cons1->getSPhi() - cons2->getSPhi()) > tol)
280  return false;
281  if (std::abs(cons1->getDPhi() - cons2->getDPhi()) > tol)
282  return false;
283 
284  return true;
285 
286  } else if (sh1->type() == "SimplePolygonBrep") {
287  const GeoSimplePolygonBrep* spb1 =
288  static_cast<const GeoSimplePolygonBrep*>(sh1);
289  const GeoSimplePolygonBrep* spb2 =
290  static_cast<const GeoSimplePolygonBrep*>(sh2);
291  if (!spb1 || !spb2)
292  return false;
293 
294  unsigned int nv1 = spb1->getNVertices();
295  unsigned int nv2 = spb2->getNVertices();
296  if (nv1 != nv2)
297  return false;
298  if (std::abs(spb1->getDZ() - spb2->getDZ()) > tol)
299  return false;
300 
301  for (unsigned int iv = 0; iv < nv1; iv++) {
302 
303  if (std::abs(spb1->getXVertex(iv) - spb2->getXVertex(iv)) > tol)
304  return false;
305  if (std::abs(spb1->getYVertex(iv) - spb2->getYVertex(iv)) > tol)
306  return false;
307  }
308 
309  return true;
310 
311  } else if (sh1->type() == "Pcon") {
312  const GeoPcon* pc1 = static_cast<const GeoPcon*>(sh1);
313  const GeoPcon* pc2 = static_cast<const GeoPcon*>(sh2);
314  if (!pc1 || !pc2)
315  return false;
316 
317  if (std::abs(pc1->getSPhi() - pc2->getSPhi()) > tol)
318  return false;
319  if (std::abs(pc1->getDPhi() - pc2->getDPhi()) > tol)
320  return false;
321 
322  unsigned int nv1 = pc1->getNPlanes();
323  unsigned int nv2 = pc2->getNPlanes();
324  if (nv1 != nv2)
325  return false;
326 
327  for (unsigned int iv = 0; iv < nv1; iv++) {
328 
329  if (std::abs(pc1->getZPlane(iv) - pc2->getZPlane(iv)) > tol)
330  return false;
331  if (std::abs(pc1->getRMinPlane(iv) - pc2->getRMinPlane(iv)) > tol)
332  return false;
333  if (std::abs(pc1->getRMaxPlane(iv) - pc2->getRMaxPlane(iv)) > tol)
334  return false;
335  }
336 
337  return true;
338 
339  } else if (sh1->type() == "Subtraction") {
340  const GeoShapeSubtraction* sub1 =
341  static_cast<const GeoShapeSubtraction*>(sh1);
342  const GeoShapeSubtraction* sub2 =
343  static_cast<const GeoShapeSubtraction*>(sh2);
344 
345  if (!sub1 || !sub2)
346  return false;
347 
348  if (!compareShapes(sub1->getOpA(), sub2->getOpA(), tol))
349  return false;
350  if (!compareShapes(sub1->getOpB(), sub2->getOpB(), tol))
351  return false;
352 
353  return true;
354 
355  } else if (sh1->type() == "Union") {
356  const GeoShapeUnion* sub1 = static_cast<const GeoShapeUnion*>(sh1);
357  const GeoShapeUnion* sub2 = static_cast<const GeoShapeUnion*>(sh2);
358 
359  if (!sub1 || !sub2)
360  return false;
361 
362  if (!compareShapes(sub1->getOpA(), sub2->getOpA(), tol))
363  return false;
364  if (!compareShapes(sub1->getOpB(), sub2->getOpB(), tol))
365  return false;
366 
367  return true;
368 
369  } else if (sh1->type() == "Shift") {
370  const GeoShapeShift* shift1 = static_cast<const GeoShapeShift*>(sh1);
371  const GeoShapeShift* shift2 = static_cast<const GeoShapeShift*>(sh2);
372 
373  if (!shift1 || !shift2)
374  return false;
375 
376  if (!compareShapes(shift1->getOp(), shift2->getOp(), tol))
377  return false;
378 
379  const GeoTrf::Transform3D& transf1 = shift1->getX();
380  const GeoTrf::Transform3D& transf2 = shift2->getX();
381 
382  if ((transf1.translation() - transf2.translation()).norm() > tol)
383  return false;
384 
385  // For rotation matrices, transpose is the same as inverse.
386  if (!identity_check(transf1.rotation() * transf2.rotation().transpose(),
387  tol))
388  return false;
389 
390  return true;
391  }
392 
393  std::cout << "unknown shape to compare:" << sh1->type() << std::endl;
394 
395  return false;
396 }
397 
398 bool Trk::GMTreeBrowser::findNamePattern(const GeoVPhysVol* gv,
399  std::string_view name) const {
400 
401  if (gv->getLogVol()->getName().find(name) != std::string::npos)
402  return true;
403 
404  for (unsigned int ic = 0; ic < gv->getNChildVols(); ic++) {
405 
406  const GeoVPhysVol* cv = &(*(gv->getChildVol(ic)));
407  if (this->findNamePattern(cv, name))
408  return true;
409  }
410 
411  return false;
412 }
413 
414 namespace {
415 
416 class GeoFindTopName : public GeoNodeAction
417 {
418 public:
419  explicit GeoFindTopName (const std::string_view name) : m_name (name)
420  {
421  }
422  virtual void handlePhysVol (const GeoPhysVol* v) override
423  { handleVol (v); }
424  virtual void handleFullPhysVol (const GeoFullPhysVol* v) override
425  { handleVol (v); }
426  const GeoVPhysVol* topName() const { return m_topName; }
427 private:
428  void handleVol (const GeoVPhysVol* v)
429  {
430  const GeoLogVol* clv = v->getLogVol();
431  if (clv && clv->getName().find(m_name) != std::string::npos) {
432  m_topName = v->getParent();
433  this->terminate();
434  }
435  }
436  const std::string_view m_name;
437  const GeoVPhysVol* m_topName = nullptr;
438 };
439 
440 }
441 
443  const GeoVPhysVol* gv, std::string_view name) {
444 
445  GeoFindTopName topName (name);
446  gv->exec (&topName);
447  return topName.topName();
448 }
449 
451  double tol) {
452 
453  if (std::abs(rotation(0, 1)) > tol)
454  return false;
455  if (std::abs(rotation(0, 2)) > tol)
456  return false;
457  if (std::abs(rotation(1, 2)) > tol)
458  return false;
459 
460  return true;
461 }
462 
464  GeoTrf::Transform3D tr_ref,
465  double tolerance) {
466  std::ios oldState(nullptr);
467  oldState.copyfmt(std::cout);
468  //
469  std::cout << std::fixed << std::setprecision(4);
470  std::cout << "test translation:x:y:z:" << tr_test.translation().x() << ":"
471  << tr_test.translation().y() << ":" << tr_test.translation().z()
472  << std::endl;
473  std::cout << " ref translation:x:y:z:" << tr_ref.translation().x() << ":"
474  << tr_ref.translation().y() << ":" << tr_ref.translation().z()
475  << std::endl;
476  std::cout << " absolute shift :"
477  << (tr_test.translation() - tr_ref.translation()).norm()
478  << ": to be compared with the tolerance limit:" << tolerance
479  << std::endl;
480  std::cout.copyfmt(oldState); // restore ostream state
481 }
482 
484  const GeoTrf::Transform3D& tr_ref,
485  double tolerance) {
486 
487  GeoTrf::RotationMatrix3D rotest = tr_test.rotation();
488  GeoTrf::RotationMatrix3D rotref = tr_ref.rotation();
489  GeoTrf::RotationMatrix3D rotid = rotest * rotref.inverse();
490  std::ios oldState(nullptr);
491  oldState.copyfmt(std::cout);
492 
493  std::cout << std::fixed << std::setprecision(4);
494  std::cout << "test rotation:" << rotest(0, 0) << ":" << rotest(0, 1) << ":"
495  << rotest(0, 2) << std::endl;
496  std::cout << " " << rotest(1, 0) << ":" << rotest(1, 1)
497  << ":" << rotest(1, 2) << std::endl;
498  std::cout << " " << rotest(2, 0) << ":" << rotest(2, 1)
499  << ":" << rotest(2, 2) << std::endl;
500  std::cout << " ref rotation:" << rotref(0, 0) << ":" << rotref(0, 1) << ":"
501  << rotref(0, 2) << std::endl;
502  std::cout << " " << rotref(1, 0) << ":" << rotref(1, 1)
503  << ":" << rotref(1, 2) << std::endl;
504  std::cout << " " << rotref(2, 0) << ":" << rotref(2, 1)
505  << ":" << rotref(2, 2) << std::endl;
506  std::cout << "test*inv(ref):" << rotid(0, 0) << ":" << rotid(0, 1) << ":"
507  << rotid(0, 2) << std::endl;
508  std::cout << " " << rotid(1, 0) << ":" << rotid(1, 1)
509  << ":" << rotid(1, 2) << std::endl;
510  std::cout << " " << rotid(2, 0) << ":" << rotid(2, 1)
511  << ":" << rotid(2, 2) << std::endl;
512  std::cout << " identity check fails within the tolerance limit:"
513  << tolerance << std::endl;
514  std::cout.copyfmt(oldState); // restore ostream state
515 }
geoGetVolumes
GeoVolumeVec_t geoGetVolumes(const GeoGraphNode *node, int depthLimit=1, int sizeHint=20)
Return the child volumes and associated transforms.
Definition: GeoVisitVolumes.cxx:211
GMTreeBrowser.h
PlotCalibFromCool.norm
norm
Definition: PlotCalibFromCool.py:100
inline_hints.h
GeoVolumeVec_t
std::vector< std::pair< const GeoVPhysVol *, GeoTrf::Transform3D > > GeoVolumeVec_t
Return the child volumes and associated transforms.
Definition: GeoVisitVolumes.h:219
Trk::GMTreeBrowser::printRotationDiff
static void printRotationDiff(const GeoTrf::Transform3D &trtest, const GeoTrf::Transform3D &trref, double tolerance)
Definition: GMTreeBrowser.cxx:483
mc.diff
diff
Definition: mc.SFGenPy8_MuMu_DD.py:14
Trk::GMTreeBrowser::findNamePattern
bool findNamePattern(const GeoVPhysVol *gv, std::string_view name) const
search of matching name patterns
Definition: GMTreeBrowser.cxx:398
python.iconfTool.models.loaders.level
level
Definition: loaders.py:20
Trk::GMTreeBrowser::identity_check
static bool identity_check(GeoTrf::RotationMatrix3D rotation, double tol)
check of rotation invariance
Definition: GMTreeBrowser.cxx:450
Trk::GMTreeBrowser::compareGeoVolumes
int compareGeoVolumes(const GeoVPhysVol *gv1, const GeoVPhysVol *gv2, double tolerance, bool printFullInfo=false, int level=0) const
Recursive comparison of trees/branches/volumes : in quiet mode (printFullInfo=False) ,...
Definition: GMTreeBrowser.cxx:42
GeoVisitVolumes.h
Visitor to process all volumes under a GeoModel node.
Trk::GMTreeBrowser::printTranslationDiff
static void printTranslationDiff(GeoTrf::Transform3D trtest, GeoTrf::Transform3D trref, double tolerance)
printout diff - unify output
Definition: GMTreeBrowser.cxx:463
ATH_FLATTEN
#define ATH_FLATTEN
Definition: inline_hints.h:52
xAOD::rotation
rotation
Definition: TrackSurface_v1.cxx:15
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition: GeoPrimitives.h:46
Trk::GMTreeBrowser::findTopBranch
static const GeoVPhysVol * findTopBranch(const GeoVPhysVol *gv, std::string_view name)
search of top branch : returns mother volume for children matching name
Definition: GMTreeBrowser.cxx:442
grepfile.ic
int ic
Definition: grepfile.py:33
tolerance
Definition: suep_shower.h:17
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:221
Amg::RotationMatrix3D
Eigen::Matrix< double, 3, 3 > RotationMatrix3D
Definition: GeoPrimitives.h:49
Trk::GMTreeBrowser::compareShapes
bool compareShapes(const GeoShape *gs1, const GeoShape *gv2, double tolerance) const
shape comparison
Definition: GMTreeBrowser.cxx:182
Trk::v
@ v
Definition: ParamDefs.h:78