ATLAS Offline Software
Loading...
Searching...
No Matches
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
9#include "TrkDetDescrGeoModelCnv/GMTreeBrowser.h"
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"
31#include "GeoModelUtilities/GeoVisitVolumes.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.
40 ATH_FLATTEN
41#endif
42int 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;
144 Trk::GMTreeBrowser::printTranslationDiff(transf1, transf2, tolerance);
145 } else
146 return diff;
147 }
148 // For rotation matrices, transpose is the same as inverse.
149 GeoTrf::RotationMatrix3D rot =
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;
160 Trk::GMTreeBrowser::printRotationDiff(transf1, transf2, tolerance);
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.
180 ATH_FLATTEN
181#endif
182bool 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
398bool 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
414namespace {
415
416class GeoFindTopName : public GeoNodeAction
417{
418public:
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; }
427private:
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
442const GeoVPhysVol* Trk::GMTreeBrowser::findTopBranch(
443 const GeoVPhysVol* gv, std::string_view name) {
444
445 GeoFindTopName topName (name);
446 gv->exec (&topName);
447 return topName.topName();
448}
449
450bool Trk::GMTreeBrowser::identity_check(GeoTrf::RotationMatrix3D rotation,
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
463void Trk::GMTreeBrowser::printTranslationDiff(GeoTrf::Transform3D tr_test,
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
483void Trk::GMTreeBrowser::printRotationDiff(const GeoTrf::Transform3D& tr_test,
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}
GeoVisitVolumes.h
Visitor to process all volumes under a GeoModel node.
geoGetVolumes
GeoVolumeVec_t geoGetVolumes(const GeoGraphNode *node, int depthLimit=1, int sizeHint=20)
Return the child volumes and associated transforms.
Definition GeoVisitVolumes.cxx:211
GeoVolumeVec_t
std::vector< std::pair< const GeoVPhysVol *, GeoTrf::Transform3D > > GeoVolumeVec_t
Return the child volumes and associated transforms.
Definition GeoVisitVolumes.h:219
diff
void diff(const Jet &rJet1, const Jet &rJet2, std::map< std::string, double > varDiff)
Difference between jets - Non-Class function required by trigger.
Definition Jet.cxx:631
Trk::GMTreeBrowser::findNamePattern
bool findNamePattern(const GeoVPhysVol *gv, std::string_view name) const
search of matching name patterns
Definition GMTreeBrowser.cxx:398
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
Trk::GMTreeBrowser::printTranslationDiff
static void printTranslationDiff(GeoTrf::Transform3D trtest, GeoTrf::Transform3D trref, double tolerance)
printout diff - unify output
Definition GMTreeBrowser.cxx:463
Trk::GMTreeBrowser::identity_check
static bool identity_check(GeoTrf::RotationMatrix3D rotation, double tol)
check of rotation invariance
Definition GMTreeBrowser.cxx:450
Trk::GMTreeBrowser::compareShapes
bool compareShapes(const GeoShape *gs1, const GeoShape *gv2, double tolerance) const
shape comparison
Definition GMTreeBrowser.cxx:182
Trk::GMTreeBrowser::printRotationDiff
static void printRotationDiff(const GeoTrf::Transform3D &trtest, const GeoTrf::Transform3D &trref, double tolerance)
Definition GMTreeBrowser.cxx:483
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
inline_hints.h
ATH_FLATTEN
#define ATH_FLATTEN
Definition inline_hints.h:52
Generated on for ATLAS Offline Software by doxygen 1.14.0