Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TGeoBoolNode.cxx
Go to the documentation of this file.
1// @(#):$Id$
2// Author: Andrei Gheata 30/05/02
3// TGeoBoolNode::Contains and parser implemented by Mihaela Gheata
4
5/*************************************************************************
6 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
7 * All rights reserved. *
8 * *
9 * For the licensing terms see $ROOTSYS/LICENSE. *
10 * For the list of contributors see $ROOTSYS/README/CREDITS. *
11 *************************************************************************/
12#include "TGeoBoolNode.h"
13
14#include <iostream>
15
16#include "TVirtualPad.h"
17#include "TVirtualViewer3D.h"
18#include "TBuffer3D.h"
19#include "TBuffer3DTypes.h"
20#include "TMath.h"
21#include "TGeoCompositeShape.h"
22#include "TGeoMatrix.h"
23#include "TGeoManager.h"
24
25/** \class TGeoBoolNode
26\ingroup Geometry_classes
27
28Base class for Boolean operations between two shapes.
29
30A Boolean node describes a Boolean operation between 'left' and 'right'
31shapes positioned with respect to an ARBITRARY reference frame. The boolean
32node is referenced by a mother composite shape and its shape components may
33be primitive but also composite shapes. The later situation leads to a binary
34tree hierarchy. When the parent composite shape is used to create a volume,
35the reference frame of the volume is chosen to match the frame in which
36node shape components were defined.
37
38The positioned shape components may or may not be disjoint. The specific
39implementations for Boolean nodes are:
40
41 - TGeoUnion - representing the Boolean union of two positioned shapes
42 - TGeoSubtraction - representing the Boolean subtraction of two positioned shapes
43 - TGeoIntersection - representing the Boolean intersection of two positioned shapes
44*/
45
47
48////////////////////////////////////////////////////////////////////////////////
49/// Constructor.
50
52
53////////////////////////////////////////////////////////////////////////////////
54/// Destructor.
55
57
58////////////////////////////////////////////////////////////////////////////////
59
61{
63 /*
64 std::lock_guard<std::mutex> guard(fMutex);
65 if (tid >= fThreadSize) {
66 Error("GetThreadData", "Thread id=%d bigger than maximum declared thread number %d. \nUse
67 TGeoManager::SetMaxThreads properly !!!", tid, fThreadSize);
68 }
69 if (tid >= fThreadSize)
70 {
71 fThreadData.resize(tid + 1);
72 fThreadSize = tid + 1;
73 }
74 if (fThreadData[tid] == 0)
75 {
76 if (fThreadData[tid] == 0)
77 fThreadData[tid] = new ThreadData_t;
78 }
79 */
80 return *fThreadData[tid];
81}
82
83////////////////////////////////////////////////////////////////////////////////
84
86{
87 std::lock_guard<std::mutex> guard(fMutex);
88 std::vector<ThreadData_t *>::iterator i = fThreadData.begin();
89 while (i != fThreadData.end()) {
90 delete *i;
91 ++i;
92 }
93 fThreadData.clear();
94 fThreadSize = 0;
95}
96
97////////////////////////////////////////////////////////////////////////////////
98/// Create thread data for n threads max.
99
101{
102 std::lock_guard<std::mutex> guard(fMutex);
103 fThreadData.resize(nthreads);
104 fThreadSize = nthreads;
105 for (Int_t tid = 0; tid < nthreads; tid++) {
106 if (fThreadData[tid] == nullptr) {
107 fThreadData[tid] = new ThreadData_t;
108 }
109 }
110 // Propagate to components
111 if (fLeft)
112 fLeft->CreateThreadData(nthreads);
113 if (fRight)
114 fRight->CreateThreadData(nthreads);
115}
116
117////////////////////////////////////////////////////////////////////////////////
118/// Set the selected branch.
119
121{
123}
124
125////////////////////////////////////////////////////////////////////////////////
126/// Default constructor
127
129{
130 fLeft = nullptr;
131 fRight = nullptr;
132 fLeftMat = nullptr;
133 fRightMat = nullptr;
134 fNpoints = 0;
135 fPoints = nullptr;
136 fThreadSize = 0;
138}
139
140////////////////////////////////////////////////////////////////////////////////
141/// Constructor called by TGeoCompositeShape providing 2 subexpressions for the 2 branches.
142
143TGeoBoolNode::TGeoBoolNode(const char *expr1, const char *expr2)
144{
145 fLeft = nullptr;
146 fRight = nullptr;
147 fLeftMat = nullptr;
148 fRightMat = nullptr;
149 fNpoints = 0;
150 fPoints = nullptr;
151 fThreadSize = 0;
153 if (!MakeBranch(expr1, kTRUE)) {
154 return;
155 }
156 if (!MakeBranch(expr2, kFALSE)) {
157 return;
158 }
159}
160
161////////////////////////////////////////////////////////////////////////////////
162/// Constructor providing left and right shapes and matrices (in the Boolean operation).
163
165{
166 fLeft = left;
167 fRight = right;
168 fLeftMat = lmat;
169 fNpoints = 0;
170 fPoints = nullptr;
171 fThreadSize = 0;
173 if (!fLeftMat)
175 else
177 fRightMat = rmat;
178 if (!fRightMat)
180 else
182 if (!fLeft) {
183 Error("ctor", "left shape is NULL");
184 return;
185 }
186 if (!fRight) {
187 Error("ctor", "right shape is NULL");
188 return;
189 }
190}
191
192////////////////////////////////////////////////////////////////////////////////
193/// Destructor.
194/// --- deletion of components handled by TGeoManager class.
195
197{
198 if (fPoints)
199 delete[] fPoints;
201}
202
203////////////////////////////////////////////////////////////////////////////////
204/// Set fPoints array
205
207{
208 if (fPoints) {
209 delete[] fPoints;
210 fPoints = nullptr;
211 fNpoints = 0;
212 }
213 if (points) {
214 fNpoints = npoints;
215 fPoints = new Double_t[3 * fNpoints];
216 memcpy(fPoints, points, 3 * fNpoints * sizeof(Double_t));
217 }
218}
219
220////////////////////////////////////////////////////////////////////////////////
221/// Expands the boolean expression either on left or right branch, creating
222/// component elements (composite shapes and boolean nodes). Returns true on success.
223
225{
226 TString sleft, sright, stransf;
227 Int_t boolop = TGeoManager::Parse(expr, sleft, sright, stransf);
228 if (boolop < 0) {
229 Error("MakeBranch", "invalid expression");
230 return kFALSE;
231 }
232 TGeoShape *shape = nullptr;
233 TGeoMatrix *mat;
234 TString newshape;
235
236 if (stransf.Length() == 0) {
237 mat = gGeoIdentity;
238 } else {
240 }
241 if (!mat) {
242 Error("MakeBranch", "transformation %s not found", stransf.Data());
243 return kFALSE;
244 }
245 switch (boolop) {
246 case 0:
247 // elementary shape
248 shape = (TGeoShape *)gGeoManager->GetListOfShapes()->FindObject(sleft.Data());
249 if (!shape) {
250 Error("MakeBranch", "shape %s not found", sleft.Data());
251 return kFALSE;
252 }
253 break;
254 case 1:
255 // composite shape - union
256 newshape = sleft;
257 newshape += "+";
258 newshape += sright;
259 shape = new TGeoCompositeShape(newshape.Data());
260 break;
261 case 2:
262 // composite shape - difference
263 newshape = sleft;
264 newshape += "-";
265 newshape += sright;
266 shape = new TGeoCompositeShape(newshape.Data());
267 break;
268 case 3:
269 // composite shape - intersection
270 newshape = sleft;
271 newshape += "*";
272 newshape += sright;
273 shape = new TGeoCompositeShape(newshape.Data());
274 break;
275 }
276 if (boolop && (!shape || !shape->IsValid())) {
277 Error("MakeBranch", "Shape %s not valid", newshape.Data());
278 if (shape)
279 delete shape;
280 return kFALSE;
281 }
282 if (left) {
283 fLeft = shape;
284 fLeftMat = mat;
285 } else {
286 fRight = shape;
287 fRightMat = mat;
288 }
289 return kTRUE;
290}
291
292////////////////////////////////////////////////////////////////////////////////
293/// Special schema for feeding the 3D buffers to the painter client.
294
296{
297 TVirtualViewer3D *viewer = gPad->GetViewer3D();
298 if (!viewer)
299 return;
300
301 // Components of composite shape hierarchies for local frame viewers are painted
302 // in coordinate frame of the top level composite shape. So we force
303 // conversion to this. See TGeoPainter::PaintNode for loading of GLMatrix.
304 Bool_t localFrame = kFALSE; // viewer->PreferLocalFrame();
305
307 TGeoHMatrix mat;
308 mat = glmat; // keep a copy
309
310 // Now perform fetch and add of the two components buffers.
311 // Note we assume that composite shapes are always completely added
312 // so don't bother to get addDaughters flag from viewer->AddObject()
313
314 // Setup matrix and fetch/add the left component buffer
315 glmat->Multiply(fLeftMat);
316 // fLeft->Paint(option);
317 if (TGeoCompositeShape *left = dynamic_cast<TGeoCompositeShape *>(fLeft)) {
318 left->PaintComposite(option);
319 } else if (fLeft) {
320 const TBuffer3D &leftBuffer = fLeft->GetBuffer3D(TBuffer3D::kAll, localFrame);
321 viewer->AddObject(leftBuffer);
322 }
323
324 // Setup matrix and fetch/add the right component buffer
325 *glmat = &mat;
326 glmat->Multiply(fRightMat);
327 // fRight->Paint(option);
328 if (TGeoCompositeShape *right = dynamic_cast<TGeoCompositeShape *>(fRight)) {
329 right->PaintComposite(option);
330 } else if (fRight) {
331 const TBuffer3D &rightBuffer = fRight->GetBuffer3D(TBuffer3D::kAll, localFrame);
332 viewer->AddObject(rightBuffer);
333 }
334
335 *glmat = &mat;
336}
337
338////////////////////////////////////////////////////////////////////////////////
339/// Register all matrices of the boolean node and descendents.
340
342{
343 if (!fLeftMat->IsIdentity())
345 if (!fRightMat->IsIdentity())
347 if (fLeft->IsComposite())
348 ((TGeoCompositeShape *)fLeft)->GetBoolNode()->RegisterMatrices();
349 if (fRight->IsComposite())
350 ((TGeoCompositeShape *)fRight)->GetBoolNode()->RegisterMatrices();
351}
352
353////////////////////////////////////////////////////////////////////////////////
354/// Replace one of the matrices. Does not work with TGeoIdentity. Returns true
355/// if replacement was successful.
356
358{
359 if (mat == gGeoIdentity || newmat == gGeoIdentity) {
360 Error("ReplaceMatrix",
361 "Matrices should not be gGeoIdentity. Use default matrix constructor to represent identities.");
362 return kFALSE;
363 }
364 if (!mat || !newmat) {
365 Error("ReplaceMatrix", "Matrices should not be null pointers.");
366 return kFALSE;
367 }
368 Bool_t replaced = kFALSE;
369 if (fLeftMat == mat) {
370 fLeftMat = newmat;
371 replaced = kTRUE;
372 }
373 if (fRightMat == mat) {
374 fRightMat = newmat;
375 replaced = kTRUE;
376 }
377 return replaced;
378}
379
380////////////////////////////////////////////////////////////////////////////////
381/// Save a primitive as a C++ statement(s) on output stream "out".
382
383void TGeoBoolNode::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
384{
387 if (!fLeftMat->IsIdentity()) {
390 }
391 if (!fRightMat->IsIdentity()) {
394 }
395}
396
397////////////////////////////////////////////////////////////////////////////////
398/// Fill buffer with shape vertices.
399
401{
402 TGeoBoolNode *bn = (TGeoBoolNode *)this;
403 Int_t npoints = bn->GetNpoints();
404 memcpy(points, fPoints, 3 * npoints * sizeof(Double_t));
405}
406
407////////////////////////////////////////////////////////////////////////////////
408/// Fill buffer with shape vertices.
409
411{
412 TGeoBoolNode *bn = (TGeoBoolNode *)this;
413 Int_t npoints = bn->GetNpoints();
414 for (Int_t i = 0; i < 3 * npoints; i++)
415 points[i] = fPoints[i];
416}
417
418////////////////////////////////////////////////////////////////////////////////
419/// Register size of this 3D object
420
422{
423 fLeft->Sizeof3D();
424 fRight->Sizeof3D();
425}
426
428
429////////////////////////////////////////////////////////////////////////////////
430/// Make a clone of this. Pointers are preserved.
431
433{
434 return new TGeoUnion(fLeft, fRight, fLeftMat, fRightMat);
435}
436
437////////////////////////////////////////////////////////////////////////////////
438/// Paint method.
439
441{
442 TVirtualViewer3D *viewer = gPad->GetViewer3D();
443
444 if (!viewer) {
445 Error("Paint", "gPad->GetViewer3D() returned 0, cannot work with composite!\n");
446 return;
447 }
448
450
452}
453
454////////////////////////////////////////////////////////////////////////////////
455/// Default constructor
456
458
459////////////////////////////////////////////////////////////////////////////////
460/// Constructor
461
462TGeoUnion::TGeoUnion(const char *expr1, const char *expr2) : TGeoBoolNode(expr1, expr2) {}
463
464////////////////////////////////////////////////////////////////////////////////
465/// Constructor providing pointers to components
466
468 : TGeoBoolNode(left, right, lmat, rmat)
469{
471 Fatal("TGeoUnion", "Unions with a half-space (%s + %s) not allowed", left->GetName(), right->GetName());
472 }
473}
474
475////////////////////////////////////////////////////////////////////////////////
476/// Destructor
477/// --- deletion of components handled by TGeoManager class.
478
480
481////////////////////////////////////////////////////////////////////////////////
482/// Compute bounding box corresponding to a union of two shapes.
483
485{
486 if (((TGeoBBox *)fLeft)->IsNullBox())
488 if (((TGeoBBox *)fRight)->IsNullBox())
490 Double_t vert[48];
491 Double_t pt[3];
492 Int_t i;
493 Double_t xmin, xmax, ymin, ymax, zmin, zmax;
494 xmin = ymin = zmin = TGeoShape::Big();
495 xmax = ymax = zmax = -TGeoShape::Big();
496 ((TGeoBBox *)fLeft)->SetBoxPoints(&vert[0]);
497 ((TGeoBBox *)fRight)->SetBoxPoints(&vert[24]);
498 for (i = 0; i < 8; i++) {
499 fLeftMat->LocalToMaster(&vert[3 * i], &pt[0]);
500 if (pt[0] < xmin)
501 xmin = pt[0];
502 if (pt[0] > xmax)
503 xmax = pt[0];
504 if (pt[1] < ymin)
505 ymin = pt[1];
506 if (pt[1] > ymax)
507 ymax = pt[1];
508 if (pt[2] < zmin)
509 zmin = pt[2];
510 if (pt[2] > zmax)
511 zmax = pt[2];
512 }
513 for (i = 8; i < 16; i++) {
514 fRightMat->LocalToMaster(&vert[3 * i], &pt[0]);
515 if (pt[0] < xmin)
516 xmin = pt[0];
517 if (pt[0] > xmax)
518 xmax = pt[0];
519 if (pt[1] < ymin)
520 ymin = pt[1];
521 if (pt[1] > ymax)
522 ymax = pt[1];
523 if (pt[2] < zmin)
524 zmin = pt[2];
525 if (pt[2] > zmax)
526 zmax = pt[2];
527 }
528 dx = 0.5 * (xmax - xmin);
529 origin[0] = 0.5 * (xmin + xmax);
530 dy = 0.5 * (ymax - ymin);
531 origin[1] = 0.5 * (ymin + ymax);
532 dz = 0.5 * (zmax - zmin);
533 origin[2] = 0.5 * (zmin + zmax);
534}
535
536////////////////////////////////////////////////////////////////////////////////
537/// Find if a union of two shapes contains a given point
538
540{
541 Double_t local[3];
542 fLeftMat->MasterToLocal(point, &local[0]);
543 Bool_t inside = fLeft->Contains(&local[0]);
544 if (inside)
545 return kTRUE;
546 fRightMat->MasterToLocal(point, &local[0]);
547 inside = fRight->Contains(&local[0]);
548 return inside;
549}
550
551////////////////////////////////////////////////////////////////////////////////
552/// Normal computation in POINT. The orientation is chosen so that DIR.dot.NORM>0.
553
554void TGeoUnion::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
555{
557 norm[0] = norm[1] = 0.;
558 norm[2] = 1.;
559 Double_t local[3];
560 Double_t ldir[3], lnorm[3];
561 if (td.fSelected == 1) {
562 fLeftMat->MasterToLocal(point, local);
563 fLeftMat->MasterToLocalVect(dir, ldir);
564 fLeft->ComputeNormal(local, ldir, lnorm);
565 fLeftMat->LocalToMasterVect(lnorm, norm);
566 return;
567 }
568 if (td.fSelected == 2) {
569 fRightMat->MasterToLocal(point, local);
570 fRightMat->MasterToLocalVect(dir, ldir);
571 fRight->ComputeNormal(local, ldir, lnorm);
572 fRightMat->LocalToMasterVect(lnorm, norm);
573 return;
574 }
575 fLeftMat->MasterToLocal(point, local);
576 if (fLeft->Contains(local)) {
577 fLeftMat->MasterToLocalVect(dir, ldir);
578 fLeft->ComputeNormal(local, ldir, lnorm);
579 fLeftMat->LocalToMasterVect(lnorm, norm);
580 return;
581 }
582 fRightMat->MasterToLocal(point, local);
583 if (fRight->Contains(local)) {
584 fRightMat->MasterToLocalVect(dir, ldir);
585 fRight->ComputeNormal(local, ldir, lnorm);
586 fRightMat->LocalToMasterVect(lnorm, norm);
587 return;
588 }
589 // Propagate forward/backward to see which of the components is intersected first
590 local[0] = point[0] + 1E-5 * dir[0];
591 local[1] = point[1] + 1E-5 * dir[1];
592 local[2] = point[2] + 1E-5 * dir[2];
593
594 if (!Contains(local)) {
595 local[0] = point[0] - 1E-5 * dir[0];
596 local[1] = point[1] - 1E-5 * dir[1];
597 local[2] = point[2] - 1E-5 * dir[2];
598 if (!Contains(local))
599 return;
600 }
601 ComputeNormal(local, dir, norm);
602}
603
604////////////////////////////////////////////////////////////////////////////////
605/// Compute minimum distance to shape vertices.
606
608{
609 return 9999;
610}
611
612////////////////////////////////////////////////////////////////////////////////
613/// Computes distance from a given point inside the shape to its boundary.
614
616TGeoUnion::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
617{
618 if (iact < 3 && safe) {
619 // compute safe distance
620 *safe = Safety(point, kTRUE);
621 if (iact == 0)
622 return TGeoShape::Big();
623 if (iact == 1 && step < *safe)
624 return TGeoShape::Big();
625 }
626
627 Double_t local[3], local1[3], master[3], ldir[3], rdir[3], pushed[3];
628 memcpy(master, point, 3 * sizeof(Double_t));
629 Int_t i;
630 TGeoBoolNode *node = (TGeoBoolNode *)this;
631 Double_t d1 = 0., d2 = 0., snxt = 0., eps = 0.;
632 fLeftMat->MasterToLocalVect(dir, ldir);
633 fRightMat->MasterToLocalVect(dir, rdir);
634 fLeftMat->MasterToLocal(point, local);
635 Bool_t inside1 = fLeft->Contains(local);
636 if (inside1)
637 d1 = fLeft->DistFromInside(local, ldir, 3);
638 else
639 memcpy(local1, local, 3 * sizeof(Double_t));
640 fRightMat->MasterToLocal(point, local);
641 Bool_t inside2 = fRight->Contains(local);
642 if (inside2)
643 d2 = fRight->DistFromInside(local, rdir, 3);
644 if (!(inside1 | inside2)) {
645 // This is a pathological case when the point is on the boundary
646 d1 = fLeft->DistFromOutside(local1, ldir, 3);
647 if (d1 < 1.E-3) {
648 eps = d1 + TGeoShape::Tolerance();
649 for (i = 0; i < 3; i++)
650 local1[i] += eps * ldir[i];
651 inside1 = kTRUE;
652 d1 = fLeft->DistFromInside(local1, ldir, 3);
653 d1 += eps;
654 } else {
655 d2 = fRight->DistFromOutside(local, rdir, 3);
656 if (d2 < 1.E-3) {
657 eps = d2 + TGeoShape::Tolerance();
658 for (i = 0; i < 3; i++)
659 local[i] += eps * rdir[i];
660 inside2 = kTRUE;
661 d2 = fRight->DistFromInside(local, rdir, 3);
662 d2 += eps;
663 }
664 }
665 }
666 while (inside1 || inside2) {
667 if (inside1 && inside2) {
668 if (d1 < d2) {
669 snxt += d1;
670 node->SetSelected(1);
671 // propagate to exit of left shape
672 inside1 = kFALSE;
673 for (i = 0; i < 3; i++)
674 master[i] += d1 * dir[i];
675 // check if propagated point is in right shape
676 fRightMat->MasterToLocal(master, local);
677 inside2 = fRight->Contains(local);
678 if (!inside2)
679 return snxt;
680 d2 = fRight->DistFromInside(local, rdir, 3);
681 if (d2 < TGeoShape::Tolerance())
682 return snxt;
683 } else {
684 snxt += d2;
685 node->SetSelected(2);
686 // propagate to exit of right shape
687 inside2 = kFALSE;
688 for (i = 0; i < 3; i++)
689 master[i] += d2 * dir[i];
690 // check if propagated point is in left shape
691 fLeftMat->MasterToLocal(master, local);
692 inside1 = fLeft->Contains(local);
693 if (!inside1)
694 return snxt;
695 d1 = fLeft->DistFromInside(local, ldir, 3);
696 if (d1 < TGeoShape::Tolerance())
697 return snxt;
698 }
699 }
700 if (inside1) {
701 snxt += d1;
702 node->SetSelected(1);
703 // propagate to exit of left shape
704 inside1 = kFALSE;
705 for (i = 0; i < 3; i++) {
706 master[i] += d1 * dir[i];
707 pushed[i] = master[i] + (1. + d1) * TGeoShape::Tolerance() * dir[i];
708 }
709 // check if propagated point is in right shape
710 fRightMat->MasterToLocal(pushed, local);
711 inside2 = fRight->Contains(local);
712 if (!inside2)
713 return snxt;
714 d2 = fRight->DistFromInside(local, rdir, 3);
715 if (d2 < TGeoShape::Tolerance())
716 return snxt;
717 d2 += (1. + d1) * TGeoShape::Tolerance();
718 }
719 if (inside2) {
720 snxt += d2;
721 node->SetSelected(2);
722 // propagate to exit of right shape
723 inside2 = kFALSE;
724 for (i = 0; i < 3; i++) {
725 master[i] += d2 * dir[i];
726 pushed[i] = master[i] + (1. + d2) * TGeoShape::Tolerance() * dir[i];
727 }
728 // check if propagated point is in left shape
729 fLeftMat->MasterToLocal(pushed, local);
730 inside1 = fLeft->Contains(local);
731 if (!inside1)
732 return snxt;
733 d1 = fLeft->DistFromInside(local, ldir, 3);
734 if (d1 < TGeoShape::Tolerance())
735 return snxt;
736 d1 += (1. + d2) * TGeoShape::Tolerance();
737 }
738 }
739 return snxt;
740}
741
742////////////////////////////////////////////////////////////////////////////////
743/// Compute distance from a given outside point to the shape.
744
746TGeoUnion::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
747{
748 if (iact < 3 && safe) {
749 // compute safe distance
750 *safe = Safety(point, kFALSE);
751 if (iact == 0)
752 return TGeoShape::Big();
753 if (iact == 1 && step < *safe)
754 return TGeoShape::Big();
755 }
756 TGeoBoolNode *node = (TGeoBoolNode *)this;
757 Double_t local[3], ldir[3], rdir[3];
758 Double_t d1, d2, snxt;
759 fLeftMat->MasterToLocal(point, &local[0]);
760 fLeftMat->MasterToLocalVect(dir, &ldir[0]);
761 fRightMat->MasterToLocalVect(dir, &rdir[0]);
762 d1 = fLeft->DistFromOutside(&local[0], &ldir[0], iact, step, safe);
763 fRightMat->MasterToLocal(point, &local[0]);
764 d2 = fRight->DistFromOutside(&local[0], &rdir[0], iact, step, safe);
765 if (d1 < d2) {
766 snxt = d1;
767 node->SetSelected(1);
768 } else {
769 snxt = d2;
770 node->SetSelected(2);
771 }
772 return snxt;
773}
774
775////////////////////////////////////////////////////////////////////////////////
776/// Returns number of vertices for the composite shape described by this union.
777
779{
780 Int_t itot = 0;
781 Double_t point[3];
782 Double_t tolerance = TGeoShape::Tolerance();
783 if (fNpoints)
784 return fNpoints;
785 // Local points for the left shape
786 Int_t nleft = fLeft->GetNmeshVertices();
787 Double_t *points1 = new Double_t[3 * nleft];
788 fLeft->SetPoints(points1);
789 // Local points for the right shape
790 Int_t nright = fRight->GetNmeshVertices();
791 Double_t *points2 = new Double_t[3 * nright];
792 fRight->SetPoints(points2);
793 Double_t *points = new Double_t[3 * (nleft + nright)];
794 for (Int_t i = 0; i < nleft; i++) {
795 if (TMath::Abs(points1[3 * i]) < tolerance && TMath::Abs(points1[3 * i + 1]) < tolerance)
796 continue;
797 fLeftMat->LocalToMaster(&points1[3 * i], &points[3 * itot]);
798 fRightMat->MasterToLocal(&points[3 * itot], point);
799 if (!fRight->Contains(point))
800 itot++;
801 }
802 for (Int_t i = 0; i < nright; i++) {
803 if (TMath::Abs(points2[3 * i]) < tolerance && TMath::Abs(points2[3 * i + 1]) < tolerance)
804 continue;
805 fRightMat->LocalToMaster(&points2[3 * i], &points[3 * itot]);
806 fLeftMat->MasterToLocal(&points[3 * itot], point);
807 if (!fLeft->Contains(point))
808 itot++;
809 }
810
811 AssignPoints(itot, points);
812
813 delete[] points1;
814 delete[] points2;
815 delete[] points;
816 return fNpoints;
817}
818
819////////////////////////////////////////////////////////////////////////////////
820/// Compute safety distance for a union node;
821
823{
824 Double_t local1[3], local2[3];
825 fLeftMat->MasterToLocal(point, local1);
826 Bool_t in1 = fLeft->Contains(local1);
827 fRightMat->MasterToLocal(point, local2);
828 Bool_t in2 = fRight->Contains(local2);
829 Bool_t intrue = in1 | in2;
830 if (intrue ^ in)
831 return 0.0;
832 Double_t saf1 = fLeft->Safety(local1, in1);
833 Double_t saf2 = fRight->Safety(local2, in2);
834 if (in1 && in2)
835 return TMath::Min(saf1, saf2);
836 if (in1)
837 return saf1;
838 if (in2)
839 return saf2;
840 return TMath::Min(saf1, saf2);
841}
842
843////////////////////////////////////////////////////////////////////////////////
844/// Save a primitive as a C++ statement(s) on output stream "out".
845
846void TGeoUnion::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
847{
849 out << " pBoolNode = new TGeoUnion(";
850 out << fLeft->GetPointerName() << ",";
851 out << fRight->GetPointerName() << ",";
852 if (!fLeftMat->IsIdentity())
853 out << fLeftMat->GetPointerName() << ",";
854 else
855 out << "0,";
856 if (!fRightMat->IsIdentity())
857 out << fRightMat->GetPointerName() << ");" << std::endl;
858 else
859 out << "0);" << std::endl;
860}
861
862////////////////////////////////////////////////////////////////////////////////
863/// Register 3D size of this shape.
864
866{
868}
869
871
872////////////////////////////////////////////////////////////////////////////////
873/// Make a clone of this. Pointers are preserved.
874
876{
878}
879
880////////////////////////////////////////////////////////////////////////////////
881/// Paint method.
882
884{
885 TVirtualViewer3D *viewer = gPad->GetViewer3D();
886
887 if (!viewer) {
888 Error("Paint", "gPad->GetViewer3D() returned 0, cannot work with composite!\n");
889 return;
890 }
891
893
895}
896
897////////////////////////////////////////////////////////////////////////////////
898/// Default constructor
899
901
902////////////////////////////////////////////////////////////////////////////////
903/// Constructor
904
905TGeoSubtraction::TGeoSubtraction(const char *expr1, const char *expr2) : TGeoBoolNode(expr1, expr2) {}
906
907////////////////////////////////////////////////////////////////////////////////
908/// Constructor providing pointers to components
909
911 : TGeoBoolNode(left, right, lmat, rmat)
912{
914 Fatal("TGeoSubstraction", "Subtractions from a half-space (%s) not allowed", left->GetName());
915 }
916}
917
918////////////////////////////////////////////////////////////////////////////////
919/// Destructor
920/// --- deletion of components handled by TGeoManager class.
921
923
924////////////////////////////////////////////////////////////////////////////////
925/// Compute bounding box corresponding to a subtraction of two shapes.
926
928{
930 if (box->IsNullBox())
932 Double_t vert[24];
933 Double_t pt[3];
934 Int_t i;
935 Double_t xmin, xmax, ymin, ymax, zmin, zmax;
936 xmin = ymin = zmin = TGeoShape::Big();
937 xmax = ymax = zmax = -TGeoShape::Big();
938 box->SetBoxPoints(&vert[0]);
939 for (i = 0; i < 8; i++) {
940 fLeftMat->LocalToMaster(&vert[3 * i], &pt[0]);
941 if (pt[0] < xmin)
942 xmin = pt[0];
943 if (pt[0] > xmax)
944 xmax = pt[0];
945 if (pt[1] < ymin)
946 ymin = pt[1];
947 if (pt[1] > ymax)
948 ymax = pt[1];
949 if (pt[2] < zmin)
950 zmin = pt[2];
951 if (pt[2] > zmax)
952 zmax = pt[2];
953 }
954 dx = 0.5 * (xmax - xmin);
955 origin[0] = 0.5 * (xmin + xmax);
956 dy = 0.5 * (ymax - ymin);
957 origin[1] = 0.5 * (ymin + ymax);
958 dz = 0.5 * (zmax - zmin);
959 origin[2] = 0.5 * (zmin + zmax);
960}
961
962////////////////////////////////////////////////////////////////////////////////
963/// Normal computation in POINT. The orientation is chosen so that DIR.dot.NORM>0.
964
965void TGeoSubtraction::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
966{
968 norm[0] = norm[1] = 0.;
969 norm[2] = 1.;
970 Double_t local[3], ldir[3], lnorm[3];
971 if (td.fSelected == 1) {
972 fLeftMat->MasterToLocal(point, local);
973 fLeftMat->MasterToLocalVect(dir, ldir);
974 fLeft->ComputeNormal(local, ldir, lnorm);
975 fLeftMat->LocalToMasterVect(lnorm, norm);
976 return;
977 }
978 if (td.fSelected == 2) {
979 fRightMat->MasterToLocal(point, local);
980 fRightMat->MasterToLocalVect(dir, ldir);
981 fRight->ComputeNormal(local, ldir, lnorm);
982 fRightMat->LocalToMasterVect(lnorm, norm);
983 return;
984 }
985 fRightMat->MasterToLocal(point, local);
986 if (fRight->Contains(local)) {
987 fRightMat->MasterToLocalVect(dir, ldir);
988 fRight->ComputeNormal(local, ldir, lnorm);
989 fRightMat->LocalToMasterVect(lnorm, norm);
990 return;
991 }
992 fLeftMat->MasterToLocal(point, local);
993 if (!fLeft->Contains(local)) {
994 fLeftMat->MasterToLocalVect(dir, ldir);
995 fLeft->ComputeNormal(local, ldir, lnorm);
996 fLeftMat->LocalToMasterVect(lnorm, norm);
997 return;
998 }
999 // point is inside left shape, but not inside the right
1000 local[0] = point[0] + 1E-5 * dir[0];
1001 local[1] = point[1] + 1E-5 * dir[1];
1002 local[2] = point[2] + 1E-5 * dir[2];
1003 if (Contains(local)) {
1004 local[0] = point[0] - 1E-5 * dir[0];
1005 local[1] = point[1] - 1E-5 * dir[1];
1006 local[2] = point[2] - 1E-5 * dir[2];
1007 if (Contains(local))
1008 return;
1009 }
1010 ComputeNormal(local, dir, norm);
1011}
1012
1013////////////////////////////////////////////////////////////////////////////////
1014/// Find if a subtraction of two shapes contains a given point
1015
1017{
1018 Double_t local[3];
1019 fLeftMat->MasterToLocal(point, &local[0]);
1020 Bool_t inside = fLeft->Contains(&local[0]);
1021 if (!inside)
1022 return kFALSE;
1023 fRightMat->MasterToLocal(point, &local[0]);
1024 inside = !fRight->Contains(&local[0]);
1025 return inside;
1026}
1027
1028////////////////////////////////////////////////////////////////////////////////
1029/// Compute minimum distance to shape vertices
1030
1032{
1033 return 9999;
1034}
1035
1036////////////////////////////////////////////////////////////////////////////////
1037/// Compute distance from a given point inside to the shape boundary.
1038
1040 Double_t *safe) const
1041{
1042 if (iact < 3 && safe) {
1043 // compute safe distance
1044 *safe = Safety(point, kTRUE);
1045 if (iact == 0)
1046 return TGeoShape::Big();
1047 if (iact == 1 && step < *safe)
1048 return TGeoShape::Big();
1049 }
1050 TGeoBoolNode *node = (TGeoBoolNode *)this;
1051 Double_t local[3], ldir[3], rdir[3];
1052 Double_t d1, d2, snxt = 0.;
1053 fLeftMat->MasterToLocal(point, &local[0]);
1054 fLeftMat->MasterToLocalVect(dir, &ldir[0]);
1055 fRightMat->MasterToLocalVect(dir, &rdir[0]);
1056 d1 = fLeft->DistFromInside(&local[0], &ldir[0], iact, step, safe);
1057 fRightMat->MasterToLocal(point, &local[0]);
1058 d2 = fRight->DistFromOutside(&local[0], &rdir[0], iact, step, safe);
1059 if (d1 < d2) {
1060 snxt = d1;
1061 node->SetSelected(1);
1062 } else {
1063 snxt = d2;
1064 node->SetSelected(2);
1065 }
1066 return snxt;
1067}
1068
1069////////////////////////////////////////////////////////////////////////////////
1070/// Compute distance from a given point outside to the shape.
1071
1073 Double_t *safe) const
1074{
1075 if (iact < 3 && safe) {
1076 // compute safe distance
1077 *safe = Safety(point, kFALSE);
1078 if (iact == 0)
1079 return TGeoShape::Big();
1080 if (iact == 1 && step < *safe)
1081 return TGeoShape::Big();
1082 }
1083 TGeoBoolNode *node = (TGeoBoolNode *)this;
1084 Double_t local[3], master[3], ldir[3], rdir[3];
1085 memcpy(&master[0], point, 3 * sizeof(Double_t));
1086 Int_t i;
1087 Double_t d1, d2, snxt = 0.;
1088 fRightMat->MasterToLocal(point, &local[0]);
1089 fLeftMat->MasterToLocalVect(dir, &ldir[0]);
1090 fRightMat->MasterToLocalVect(dir, &rdir[0]);
1091 // check if inside '-'
1092 Bool_t inside = fRight->Contains(&local[0]);
1093 Double_t epsil = 0.;
1094 while (true) {
1095 if (inside) {
1096 // propagate to outside of '-'
1097 node->SetSelected(2);
1098 d1 = fRight->DistFromInside(&local[0], &rdir[0], iact, step, safe);
1099 snxt += d1 + epsil;
1100 for (i = 0; i < 3; i++)
1101 master[i] += (d1 + 1E-8) * dir[i];
1102 epsil = 1.E-8;
1103 // now master outside '-'; check if inside '+'
1104 fLeftMat->MasterToLocal(&master[0], &local[0]);
1105 if (fLeft->Contains(&local[0]))
1106 return snxt;
1107 }
1108 // master outside '-' and outside '+' ; find distances to both
1109 node->SetSelected(1);
1110 fLeftMat->MasterToLocal(&master[0], &local[0]);
1111 d2 = fLeft->DistFromOutside(&local[0], &ldir[0], iact, step, safe);
1112 if (d2 > 1E20)
1113 return TGeoShape::Big();
1114
1115 fRightMat->MasterToLocal(&master[0], &local[0]);
1116 d1 = fRight->DistFromOutside(&local[0], &rdir[0], iact, step, safe);
1117 if (d2 < d1 - TGeoShape::Tolerance()) {
1118 snxt += d2 + epsil;
1119 return snxt;
1120 }
1121 // propagate to '-'
1122 snxt += d1 + epsil;
1123 for (i = 0; i < 3; i++)
1124 master[i] += (d1 + 1E-8) * dir[i];
1125 epsil = 1.E-8;
1126 // now inside '-' and not inside '+'
1127 fRightMat->MasterToLocal(&master[0], &local[0]);
1128 inside = kTRUE;
1129 }
1130}
1131
1132////////////////////////////////////////////////////////////////////////////////
1133/// Returns number of vertices for the composite shape described by this subtraction.
1134
1136{
1137 Int_t itot = 0;
1138 Double_t point[3];
1139 Double_t tolerance = TGeoShape::Tolerance();
1140 if (fNpoints)
1141 return fNpoints;
1142 Int_t nleft = fLeft->GetNmeshVertices();
1143 Int_t nright = fRight->GetNmeshVertices();
1144 Double_t *points = new Double_t[3 * (nleft + nright)];
1145 Double_t *points1 = new Double_t[3 * nleft];
1146 fLeft->SetPoints(points1);
1147 for (Int_t i = 0; i < nleft; i++) {
1148 if (TMath::Abs(points1[3 * i]) < tolerance && TMath::Abs(points1[3 * i + 1]) < tolerance)
1149 continue;
1150 fLeftMat->LocalToMaster(&points1[3 * i], &points[3 * itot]);
1151 fRightMat->MasterToLocal(&points[3 * itot], point);
1152 if (!fRight->Contains(point))
1153 itot++;
1154 }
1155 Double_t *points2 = new Double_t[3 * nright];
1156 fRight->SetPoints(points2);
1157 for (Int_t i = 0; i < nright; i++) {
1158 if (TMath::Abs(points2[3 * i]) < tolerance && TMath::Abs(points2[3 * i + 1]) < tolerance)
1159 continue;
1160 fRightMat->LocalToMaster(&points2[3 * i], &points[3 * itot]);
1161 fLeftMat->MasterToLocal(&points[3 * itot], point);
1162 if (fLeft->Contains(point))
1163 itot++;
1164 }
1165
1166 AssignPoints(itot, points);
1167
1168 delete[] points1;
1169 delete[] points2;
1170 delete[] points;
1171 return fNpoints;
1172}
1173
1174////////////////////////////////////////////////////////////////////////////////
1175/// Compute safety distance for a union node;
1176
1178{
1179 Double_t local1[3], local2[3];
1180 fLeftMat->MasterToLocal(point, local1);
1181 Bool_t in1 = fLeft->Contains(local1);
1182 fRightMat->MasterToLocal(point, local2);
1183 Bool_t in2 = fRight->Contains(local2);
1184 Bool_t intrue = in1 && (!in2);
1185 if (in ^ intrue)
1186 return 0.0;
1187 Double_t saf1 = fLeft->Safety(local1, in1);
1188 Double_t saf2 = fRight->Safety(local2, in2);
1189 if (in1 && in2)
1190 return saf2;
1191 if (in1)
1192 return TMath::Min(saf1, saf2);
1193 if (in2)
1194 return TMath::Max(saf1, saf2);
1195 return saf1;
1196}
1197
1198////////////////////////////////////////////////////////////////////////////////
1199/// Save a primitive as a C++ statement(s) on output stream "out".
1200
1201void TGeoSubtraction::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
1202{
1204 out << " pBoolNode = new TGeoSubtraction(";
1205 out << fLeft->GetPointerName() << ",";
1206 out << fRight->GetPointerName() << ",";
1207 if (!fLeftMat->IsIdentity())
1208 out << fLeftMat->GetPointerName() << ",";
1209 else
1210 out << "0,";
1211 if (!fRightMat->IsIdentity())
1212 out << fRightMat->GetPointerName() << ");" << std::endl;
1213 else
1214 out << "0);" << std::endl;
1215}
1216
1217////////////////////////////////////////////////////////////////////////////////
1218/// Register 3D size of this shape.
1219
1221{
1223}
1224
1226
1227////////////////////////////////////////////////////////////////////////////////
1228/// Make a clone of this. Pointers are preserved.
1229
1231{
1233}
1234
1235////////////////////////////////////////////////////////////////////////////////
1236/// Paint method.
1237
1239{
1240 TVirtualViewer3D *viewer = gPad->GetViewer3D();
1241
1242 if (!viewer) {
1243 Error("Paint", "gPad->GetViewer3D() returned 0, cannot work with composite!\n");
1244 return;
1245 }
1246
1248
1250}
1251
1252////////////////////////////////////////////////////////////////////////////////
1253/// Default constructor
1254
1256
1257////////////////////////////////////////////////////////////////////////////////
1258/// Constructor
1259
1260TGeoIntersection::TGeoIntersection(const char *expr1, const char *expr2) : TGeoBoolNode(expr1, expr2) {}
1261
1262////////////////////////////////////////////////////////////////////////////////
1263/// Constructor providing pointers to components
1264
1266 : TGeoBoolNode(left, right, lmat, rmat)
1267{
1270 if (hs1 && hs2)
1271 Fatal("ctor", "cannot intersect two half-spaces: %s * %s", left->GetName(), right->GetName());
1272}
1273
1274////////////////////////////////////////////////////////////////////////////////
1275/// Destructor
1276/// --- deletion of components handled by TGeoManager class.
1277
1279
1280////////////////////////////////////////////////////////////////////////////////
1281/// Compute bounding box corresponding to a intersection of two shapes.
1282
1284{
1287 Double_t vert[48];
1288 Double_t pt[3];
1289 Int_t i;
1290 Double_t xmin1, xmax1, ymin1, ymax1, zmin1, zmax1;
1291 Double_t xmin2, xmax2, ymin2, ymax2, zmin2, zmax2;
1292 xmin1 = ymin1 = zmin1 = xmin2 = ymin2 = zmin2 = TGeoShape::Big();
1293 xmax1 = ymax1 = zmax1 = xmax2 = ymax2 = zmax2 = -TGeoShape::Big();
1294 if (!hs1) {
1295 if (((TGeoBBox *)fLeft)->IsNullBox())
1296 fLeft->ComputeBBox();
1297 ((TGeoBBox *)fLeft)->SetBoxPoints(&vert[0]);
1298 for (i = 0; i < 8; i++) {
1299 fLeftMat->LocalToMaster(&vert[3 * i], &pt[0]);
1300 if (pt[0] < xmin1)
1301 xmin1 = pt[0];
1302 if (pt[0] > xmax1)
1303 xmax1 = pt[0];
1304 if (pt[1] < ymin1)
1305 ymin1 = pt[1];
1306 if (pt[1] > ymax1)
1307 ymax1 = pt[1];
1308 if (pt[2] < zmin1)
1309 zmin1 = pt[2];
1310 if (pt[2] > zmax1)
1311 zmax1 = pt[2];
1312 }
1313 }
1314 if (!hs2) {
1315 if (((TGeoBBox *)fRight)->IsNullBox())
1317 ((TGeoBBox *)fRight)->SetBoxPoints(&vert[24]);
1318 for (i = 8; i < 16; i++) {
1319 fRightMat->LocalToMaster(&vert[3 * i], &pt[0]);
1320 if (pt[0] < xmin2)
1321 xmin2 = pt[0];
1322 if (pt[0] > xmax2)
1323 xmax2 = pt[0];
1324 if (pt[1] < ymin2)
1325 ymin2 = pt[1];
1326 if (pt[1] > ymax2)
1327 ymax2 = pt[1];
1328 if (pt[2] < zmin2)
1329 zmin2 = pt[2];
1330 if (pt[2] > zmax2)
1331 zmax2 = pt[2];
1332 }
1333 }
1334 if (hs1) {
1335 dx = 0.5 * (xmax2 - xmin2);
1336 origin[0] = 0.5 * (xmax2 + xmin2);
1337 dy = 0.5 * (ymax2 - ymin2);
1338 origin[1] = 0.5 * (ymax2 + ymin2);
1339 dz = 0.5 * (zmax2 - zmin2);
1340 origin[2] = 0.5 * (zmax2 + zmin2);
1341 return;
1342 }
1343 if (hs2) {
1344 dx = 0.5 * (xmax1 - xmin1);
1345 origin[0] = 0.5 * (xmax1 + xmin1);
1346 dy = 0.5 * (ymax1 - ymin1);
1347 origin[1] = 0.5 * (ymax1 + ymin1);
1348 dz = 0.5 * (zmax1 - zmin1);
1349 origin[2] = 0.5 * (zmax1 + zmin1);
1350 return;
1351 }
1352 Double_t sort[4];
1353 Int_t isort[4];
1354 sort[0] = xmin1;
1355 sort[1] = xmax1;
1356 sort[2] = xmin2;
1357 sort[3] = xmax2;
1358 TMath::Sort(4, &sort[0], &isort[0], kFALSE);
1359 if (isort[1] % 2) {
1360 Warning("ComputeBBox", "shapes %s and %s do not intersect", fLeft->GetName(), fRight->GetName());
1361 dx = dy = dz = 0;
1362 memset(origin, 0, 3 * sizeof(Double_t));
1363 return;
1364 }
1365 dx = 0.5 * (sort[isort[2]] - sort[isort[1]]);
1366 origin[0] = 0.5 * (sort[isort[1]] + sort[isort[2]]);
1367 sort[0] = ymin1;
1368 sort[1] = ymax1;
1369 sort[2] = ymin2;
1370 sort[3] = ymax2;
1371 TMath::Sort(4, &sort[0], &isort[0], kFALSE);
1372 if (isort[1] % 2) {
1373 Warning("ComputeBBox", "shapes %s and %s do not intersect", fLeft->GetName(), fRight->GetName());
1374 dx = dy = dz = 0;
1375 memset(origin, 0, 3 * sizeof(Double_t));
1376 return;
1377 }
1378 dy = 0.5 * (sort[isort[2]] - sort[isort[1]]);
1379 origin[1] = 0.5 * (sort[isort[1]] + sort[isort[2]]);
1380 sort[0] = zmin1;
1381 sort[1] = zmax1;
1382 sort[2] = zmin2;
1383 sort[3] = zmax2;
1384 TMath::Sort(4, &sort[0], &isort[0], kFALSE);
1385 if (isort[1] % 2) {
1386 Warning("ComputeBBox", "shapes %s and %s do not intersect", fLeft->GetName(), fRight->GetName());
1387 dx = dy = dz = 0;
1388 memset(origin, 0, 3 * sizeof(Double_t));
1389 return;
1390 }
1391 dz = 0.5 * (sort[isort[2]] - sort[isort[1]]);
1392 origin[2] = 0.5 * (sort[isort[1]] + sort[isort[2]]);
1393}
1394
1395////////////////////////////////////////////////////////////////////////////////
1396/// Normal computation in POINT. The orientation is chosen so that DIR.dot.NORM>0.
1397
1398void TGeoIntersection::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
1399{
1401 Double_t local[3], ldir[3], lnorm[3];
1402 norm[0] = norm[1] = 0.;
1403 norm[2] = 1.;
1404 if (td.fSelected == 1) {
1405 fLeftMat->MasterToLocal(point, local);
1406 fLeftMat->MasterToLocalVect(dir, ldir);
1407 fLeft->ComputeNormal(local, ldir, lnorm);
1408 fLeftMat->LocalToMasterVect(lnorm, norm);
1409 return;
1410 }
1411 if (td.fSelected == 2) {
1412 fRightMat->MasterToLocal(point, local);
1413 fRightMat->MasterToLocalVect(dir, ldir);
1414 fRight->ComputeNormal(local, ldir, lnorm);
1415 fRightMat->LocalToMasterVect(lnorm, norm);
1416 return;
1417 }
1418 fLeftMat->MasterToLocal(point, local);
1419 if (!fLeft->Contains(local)) {
1420 fLeftMat->MasterToLocalVect(dir, ldir);
1421 fLeft->ComputeNormal(local, ldir, lnorm);
1422 fLeftMat->LocalToMasterVect(lnorm, norm);
1423 return;
1424 }
1425 fRightMat->MasterToLocal(point, local);
1426 if (!fRight->Contains(local)) {
1427 fRightMat->MasterToLocalVect(dir, ldir);
1428 fRight->ComputeNormal(local, ldir, lnorm);
1429 fRightMat->LocalToMasterVect(lnorm, norm);
1430 return;
1431 }
1432 // point is inside intersection.
1433 local[0] = point[0] + 1E-5 * dir[0];
1434 local[1] = point[1] + 1E-5 * dir[1];
1435 local[2] = point[2] + 1E-5 * dir[2];
1436 if (Contains(local)) {
1437 local[0] = point[0] - 1E-5 * dir[0];
1438 local[1] = point[1] - 1E-5 * dir[1];
1439 local[2] = point[2] - 1E-5 * dir[2];
1440 if (Contains(local))
1441 return;
1442 }
1443 ComputeNormal(local, dir, norm);
1444}
1445
1446////////////////////////////////////////////////////////////////////////////////
1447/// Find if a intersection of two shapes contains a given point
1448
1450{
1451 Double_t local[3];
1452 fLeftMat->MasterToLocal(point, &local[0]);
1453 Bool_t inside = fLeft->Contains(&local[0]);
1454 if (!inside)
1455 return kFALSE;
1456 fRightMat->MasterToLocal(point, &local[0]);
1457 inside = fRight->Contains(&local[0]);
1458 return inside;
1459}
1460
1461////////////////////////////////////////////////////////////////////////////////
1462/// Compute minimum distance to shape vertices
1463
1465{
1466 return 9999;
1467}
1468
1469////////////////////////////////////////////////////////////////////////////////
1470/// Compute distance from a given point inside to the shape boundary.
1471
1473 Double_t *safe) const
1474{
1475 if (iact < 3 && safe) {
1476 // compute safe distance
1477 *safe = Safety(point, kTRUE);
1478 if (iact == 0)
1479 return TGeoShape::Big();
1480 if (iact == 1 && step < *safe)
1481 return TGeoShape::Big();
1482 }
1483 TGeoBoolNode *node = (TGeoBoolNode *)this;
1484 Double_t local[3], ldir[3], rdir[3];
1485 Double_t d1, d2, snxt = 0.;
1486 fLeftMat->MasterToLocal(point, &local[0]);
1487 fLeftMat->MasterToLocalVect(dir, &ldir[0]);
1488 fRightMat->MasterToLocalVect(dir, &rdir[0]);
1489 d1 = fLeft->DistFromInside(&local[0], &ldir[0], iact, step, safe);
1490 fRightMat->MasterToLocal(point, &local[0]);
1491 d2 = fRight->DistFromInside(&local[0], &rdir[0], iact, step, safe);
1492 if (d1 < d2) {
1493 snxt = d1;
1494 node->SetSelected(1);
1495 } else {
1496 snxt = d2;
1497 node->SetSelected(2);
1498 }
1499 return snxt;
1500}
1501
1502////////////////////////////////////////////////////////////////////////////////
1503/// Compute distance from a given point outside to the shape.
1504
1506 Double_t *safe) const
1507{
1509 if (iact < 3 && safe) {
1510 // compute safe distance
1511 *safe = Safety(point, kFALSE);
1512 if (iact == 0)
1513 return TGeoShape::Big();
1514 if (iact == 1 && step < *safe)
1515 return TGeoShape::Big();
1516 }
1517 TGeoBoolNode *node = (TGeoBoolNode *)this;
1518 Double_t lpt[3], rpt[3], master[3], ldir[3], rdir[3];
1519 memcpy(master, point, 3 * sizeof(Double_t));
1520 Int_t i;
1521 Double_t d1 = 0.;
1522 Double_t d2 = 0.;
1523 fLeftMat->MasterToLocal(point, lpt);
1524 fRightMat->MasterToLocal(point, rpt);
1525 fLeftMat->MasterToLocalVect(dir, ldir);
1526 fRightMat->MasterToLocalVect(dir, rdir);
1527 Bool_t inleft = fLeft->Contains(lpt);
1528 Bool_t inright = fRight->Contains(rpt);
1529 node->SetSelected(0);
1530 Double_t snext = 0.0;
1531 if (inleft && inright) {
1532 // It is vey likely to have a numerical issue and the point should
1533 // be logically outside one of the shapes
1534 d1 = fLeft->DistFromInside(lpt, ldir, 3);
1535 d2 = fRight->DistFromInside(rpt, rdir, 3);
1536 if (d1 < 1.E-3)
1537 inleft = kFALSE;
1538 if (d2 < 1.E-3)
1539 inright = kFALSE;
1540 if (inleft && inright)
1541 return snext;
1542 }
1543
1544 while (true) {
1545 d1 = d2 = 0;
1546 if (!inleft) {
1547 d1 = fLeft->DistFromOutside(lpt, ldir, 3);
1548 d1 = TMath::Max(d1, tol);
1549 if (d1 > 1E20)
1550 return TGeoShape::Big();
1551 }
1552 if (!inright) {
1553 d2 = fRight->DistFromOutside(rpt, rdir, 3);
1554 d2 = TMath::Max(d2, tol);
1555 if (d2 > 1E20)
1556 return TGeoShape::Big();
1557 }
1558
1559 if (d1 > d2) {
1560 // propagate to left shape
1561 snext += d1;
1562 node->SetSelected(1);
1563 inleft = kTRUE;
1564 for (i = 0; i < 3; i++)
1565 master[i] += d1 * dir[i];
1566 fRightMat->MasterToLocal(master, rpt);
1567 // Push rpt to avoid a bad boundary condition
1568 for (i = 0; i < 3; i++)
1569 rpt[i] += tol * rdir[i];
1570 // check if propagated point is inside right shape
1571 inright = fRight->Contains(rpt);
1572 if (inright)
1573 return snext;
1574 // here inleft=true, inright=false
1575 } else {
1576 // propagate to right shape
1577 snext += d2;
1578 node->SetSelected(2);
1579 inright = kTRUE;
1580 for (i = 0; i < 3; i++)
1581 master[i] += d2 * dir[i];
1582 fLeftMat->MasterToLocal(master, lpt);
1583 // Push lpt to avoid a bad boundary condition
1584 for (i = 0; i < 3; i++)
1585 lpt[i] += tol * ldir[i];
1586 // check if propagated point is inside left shape
1587 inleft = fLeft->Contains(lpt);
1588 if (inleft)
1589 return snext;
1590 // here inleft=false, inright=true
1591 }
1592 }
1593 return snext;
1594}
1595
1596////////////////////////////////////////////////////////////////////////////////
1597/// Returns number of vertices for the composite shape described by this intersection.
1598
1600{
1601 Int_t itot = 0;
1602 Double_t point[3];
1603 Double_t tolerance = TGeoShape::Tolerance();
1604 if (fNpoints)
1605 return fNpoints;
1606 Int_t nleft = fLeft->GetNmeshVertices();
1607 Int_t nright = fRight->GetNmeshVertices();
1608 Double_t *points = new Double_t[3 * (nleft + nright)];
1609 Double_t *points1 = new Double_t[3 * nleft];
1610 fLeft->SetPoints(points1);
1611 for (Int_t i = 0; i < nleft; i++) {
1612 if (TMath::Abs(points1[3 * i]) < tolerance && TMath::Abs(points1[3 * i + 1]) < tolerance)
1613 continue;
1614 fLeftMat->LocalToMaster(&points1[3 * i], &points[3 * itot]);
1615 fRightMat->MasterToLocal(&points[3 * itot], point);
1616 if (fRight->Contains(point))
1617 itot++;
1618 }
1619 Double_t *points2 = new Double_t[3 * nright];
1620 fRight->SetPoints(points2);
1621 for (Int_t i = 0; i < nright; i++) {
1622 if (TMath::Abs(points2[3 * i]) < tolerance && TMath::Abs(points2[3 * i + 1]) < tolerance)
1623 continue;
1624 fRightMat->LocalToMaster(&points2[3 * i], &points[3 * itot]);
1625 fLeftMat->MasterToLocal(&points[3 * itot], point);
1626 if (fLeft->Contains(point))
1627 itot++;
1628 }
1629
1630 AssignPoints(itot, points);
1631
1632 delete[] points1;
1633 delete[] points2;
1634 delete[] points;
1635 return fNpoints;
1636}
1637
1638////////////////////////////////////////////////////////////////////////////////
1639/// Compute safety distance for a union node;
1640
1642{
1643 Double_t local1[3], local2[3];
1644 fLeftMat->MasterToLocal(point, local1);
1645 Bool_t in1 = fLeft->Contains(local1);
1646 fRightMat->MasterToLocal(point, local2);
1647 Bool_t in2 = fRight->Contains(local2);
1648 Bool_t intrue = in1 & in2;
1649 if (in ^ intrue)
1650 return 0.0;
1651 Double_t saf1 = fLeft->Safety(local1, in1);
1652 Double_t saf2 = fRight->Safety(local2, in2);
1653 if (in1 && in2)
1654 return TMath::Min(saf1, saf2);
1655 if (in1)
1656 return saf2;
1657 if (in2)
1658 return saf1;
1659 return TMath::Max(saf1, saf2);
1660}
1661
1662////////////////////////////////////////////////////////////////////////////////
1663/// Save a primitive as a C++ statement(s) on output stream "out".
1664
1665void TGeoIntersection::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
1666{
1668 out << " pBoolNode = new TGeoIntersection(";
1669 out << fLeft->GetPointerName() << ",";
1670 out << fRight->GetPointerName() << ",";
1671 if (!fLeftMat->IsIdentity())
1672 out << fLeftMat->GetPointerName() << ",";
1673 else
1674 out << "0,";
1675 if (!fRightMat->IsIdentity())
1676 out << fRightMat->GetPointerName() << ");" << std::endl;
1677 else
1678 out << "0);" << std::endl;
1679}
1680
1681////////////////////////////////////////////////////////////////////////////////
1682/// Register 3D size of this shape.
1683
1685{
1687}
float Float_t
Definition RtypesCore.h:57
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
double Double_t
Definition RtypesCore.h:59
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:377
Option_t Option_t option
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t sel
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t points
R__EXTERN TGeoManager * gGeoManager
R__EXTERN TGeoIdentity * gGeoIdentity
Definition TGeoMatrix.h:537
float xmin
float ymin
float xmax
float ymax
#define gPad
Generic 3D primitive description class.
Definition TBuffer3D.h:18
@ kCSDifference
Definition TBuffer3D.h:43
@ kCSIntersection
Definition TBuffer3D.h:43
Box class.
Definition TGeoBBox.h:17
Base class for Boolean operations between two shapes.
virtual void Sizeof3D() const
Register size of this 3D object.
Bool_t MakeBranch(const char *expr, Bool_t left)
Mutex for thread data access.
TGeoMatrix * fLeftMat
void ClearThreadData() const
std::vector< ThreadData_t * > fThreadData
array of mesh points
~TGeoBoolNode() override
Destructor.
TGeoShape * fLeft
Bool_t ReplaceMatrix(TGeoMatrix *mat, TGeoMatrix *newmat)
Replace one of the matrices.
void AssignPoints(Int_t npoints, Double_t *points)
Set fPoints array.
void SavePrimitive(std::ostream &out, Option_t *option="") override
Save a primitive as a C++ statement(s) on output stream "out".
void CreateThreadData(Int_t nthreads)
Create thread data for n threads max.
void Paint(Option_t *option) override
Special schema for feeding the 3D buffers to the painter client.
virtual void SetPoints(Double_t *points) const
Fill buffer with shape vertices.
Int_t fThreadSize
Navigation data per thread.
std::mutex fMutex
Size for the navigation data array.
void RegisterMatrices()
Register all matrices of the boolean node and descendents.
TGeoShape * fRight
virtual Int_t GetNpoints()=0
Double_t * fPoints
number of points on the mesh
TGeoBoolNode()
Default constructor.
TGeoMatrix * fRightMat
ThreadData_t & GetThreadData() const
void SetSelected(Int_t sel)
Set the selected branch.
Composite shapes are Boolean combinations of two or more shape components.
Matrix class used for computing global transformations Should NOT be used for node definition.
Definition TGeoMatrix.h:458
void Multiply(const TGeoMatrix *right)
multiply to the right with an other transformation if right is identity matrix, just return
Int_t DistanceToPrimitive(Int_t px, Int_t py) override
Compute minimum distance to shape vertices.
TGeoBoolNode * MakeClone() const override
Make a clone of this. Pointers are preserved.
TGeoIntersection()
Default constructor.
void Sizeof3D() const override
Register 3D size of this shape.
void ComputeBBox(Double_t &dx, Double_t &dy, Double_t &dz, Double_t *origin) override
Compute bounding box corresponding to a intersection of two shapes.
Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=0, Double_t *safe=nullptr) const override
Compute distance from a given point outside to the shape.
Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const override
Compute safety distance for a union node;.
void SavePrimitive(std::ostream &out, Option_t *option="") override
Save a primitive as a C++ statement(s) on output stream "out".
Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=0, Double_t *safe=nullptr) const override
Compute distance from a given point inside to the shape boundary.
void Paint(Option_t *option) override
Paint method.
~TGeoIntersection() override
Destructor — deletion of components handled by TGeoManager class.
void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm) override
Normal computation in POINT. The orientation is chosen so that DIR.dot.NORM>0.
Bool_t Contains(const Double_t *point) const override
Find if a intersection of two shapes contains a given point.
Int_t GetNpoints() override
Returns number of vertices for the composite shape described by this intersection.
TObjArray * GetListOfMatrices() const
static Int_t Parse(const char *expr, TString &expr1, TString &expr2, TString &expr3)
Parse a string boolean expression and do a syntax check.
TObjArray * GetListOfShapes() const
static Int_t ThreadId()
Translates the current thread id to an ordinal number.
Geometrical transformation package.
Definition TGeoMatrix.h:38
virtual void LocalToMasterVect(const Double_t *local, Double_t *master) const
convert a vector by multiplying its column vector (x, y, z, 1) to matrix inverse
virtual void MasterToLocal(const Double_t *master, Double_t *local) const
convert a point by multiplying its column vector (x, y, z, 1) to matrix
virtual void MasterToLocalVect(const Double_t *master, Double_t *local) const
convert a point by multiplying its column vector (x, y, z, 1) to matrix
virtual void RegisterYourself()
Register the matrix in the current manager, which will become the owner.
virtual void LocalToMaster(const Double_t *local, Double_t *master) const
convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse
Bool_t IsIdentity() const
Definition TGeoMatrix.h:63
const char * GetPointerName() const
Provide a pointer name containing uid.
Base abstract class for all shapes.
Definition TGeoShape.h:25
virtual void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)=0
virtual const TBuffer3D & GetBuffer3D(Int_t reqSections, Bool_t localFrame) const
Stub implementation to avoid forcing implementation at this stage.
static Double_t Big()
Definition TGeoShape.h:87
virtual Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const =0
virtual void CreateThreadData(Int_t)
Definition TGeoShape.h:67
Bool_t IsValid() const
Definition TGeoShape.h:143
virtual Int_t GetNmeshVertices() const
Definition TGeoShape.h:127
virtual void Sizeof3D() const =0
virtual Bool_t IsComposite() const
Definition TGeoShape.h:130
virtual Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const =0
const char * GetPointerName() const
Provide a pointer name containing uid.
virtual Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=nullptr) const =0
const char * GetName() const override
Get the shape name.
virtual void ComputeBBox()=0
static TGeoMatrix * GetTransform()
Returns current transformation matrix that applies to shape.
virtual Bool_t Contains(const Double_t *point) const =0
@ kGeoHalfSpace
Definition TGeoShape.h:62
static Double_t Tolerance()
Definition TGeoShape.h:90
virtual void SetPoints(Double_t *points) const =0
Bool_t TestShapeBit(UInt_t f) const
Definition TGeoShape.h:167
TGeoSubtraction()
Default constructor.
void ComputeBBox(Double_t &dx, Double_t &dy, Double_t &dz, Double_t *origin) override
Compute bounding box corresponding to a subtraction of two shapes.
~TGeoSubtraction() override
Destructor — deletion of components handled by TGeoManager class.
Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const override
Compute safety distance for a union node;.
void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm) override
Normal computation in POINT. The orientation is chosen so that DIR.dot.NORM>0.
TGeoBoolNode * MakeClone() const override
Make a clone of this. Pointers are preserved.
Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=0, Double_t *safe=nullptr) const override
Compute distance from a given point outside to the shape.
void SavePrimitive(std::ostream &out, Option_t *option="") override
Save a primitive as a C++ statement(s) on output stream "out".
Bool_t Contains(const Double_t *point) const override
Find if a subtraction of two shapes contains a given point.
void Sizeof3D() const override
Register 3D size of this shape.
Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=0, Double_t *safe=nullptr) const override
Compute distance from a given point inside to the shape boundary.
void Paint(Option_t *option) override
Paint method.
Int_t DistanceToPrimitive(Int_t px, Int_t py) override
Compute minimum distance to shape vertices.
Int_t GetNpoints() override
Returns number of vertices for the composite shape described by this subtraction.
Int_t DistanceToPrimitive(Int_t px, Int_t py) override
Compute minimum distance to shape vertices.
TGeoBoolNode * MakeClone() const override
Make a clone of this. Pointers are preserved.
void ComputeBBox(Double_t &dx, Double_t &dy, Double_t &dz, Double_t *origin) override
Compute bounding box corresponding to a union of two shapes.
Int_t GetNpoints() override
Returns number of vertices for the composite shape described by this union.
Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=0, Double_t *safe=nullptr) const override
Compute distance from a given outside point to the shape.
Bool_t Contains(const Double_t *point) const override
Find if a union of two shapes contains a given point.
Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const override
Compute safety distance for a union node;.
Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=0, Double_t *safe=nullptr) const override
Computes distance from a given point inside the shape to its boundary.
~TGeoUnion() override
Destructor — deletion of components handled by TGeoManager class.
TGeoUnion()
Default constructor.
void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm) override
Normal computation in POINT. The orientation is chosen so that DIR.dot.NORM>0.
void Paint(Option_t *option) override
Paint method.
void SavePrimitive(std::ostream &out, Option_t *option="") override
Save a primitive as a C++ statement(s) on output stream "out".
void Sizeof3D() const override
Register 3D size of this shape.
TObject * FindObject(const char *name) const override
Find an object in this collection using its name.
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:962
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save a primitive as a C++ statement(s) on output stream "out".
Definition TObject.cxx:751
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:976
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition TObject.cxx:1004
Basic string class.
Definition TString.h:139
Ssiz_t Length() const
Definition TString.h:417
const char * Data() const
Definition TString.h:376
Abstract 3D shapes viewer.
virtual void AddCompositeOp(UInt_t operation)=0
virtual Int_t AddObject(const TBuffer3D &buffer, Bool_t *addChildren=nullptr)=0
TPaveText * pt
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Definition fillpatterns.C:1
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:250
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:198
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Sort the n elements of the array a of generic templated type Element.
Definition TMathBase.h:431
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123