Logo ROOT  
Reference Guide
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/*************************************************************************
7 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
8 * All rights reserved. *
9 * *
10 * For the licensing terms see $ROOTSYS/LICENSE. *
11 * For the list of contributors see $ROOTSYS/README/CREDITS. *
12 *************************************************************************/
13#include "TGeoBoolNode.h"
14
15#include "Riostream.h"
16
17#include "TVirtualPad.h"
18#include "TVirtualViewer3D.h"
19#include "TBuffer3D.h"
20#include "TBuffer3DTypes.h"
21#include "TMath.h"
22#include "TGeoCompositeShape.h"
23#include "TGeoMatrix.h"
24#include "TGeoManager.h"
25
26/** \class TGeoBoolNode
27\ingroup Geometry_classes
28
29Base class for Boolean operations between two shapes.
30
31A Boolean node describes a Boolean operation between 'left' and 'right'
32shapes positioned with respect to an ARBITRARY reference frame. The boolean
33node is referenced by a mother composite shape and its shape components may
34be primitive but also composite shapes. The later situation leads to a binary
35tree hierarchy. When the parent composite shape is used to create a volume,
36the reference frame of the volume is chosen to match the frame in which
37node shape components were defined.
38
39The positioned shape components may or may not be disjoint. The specific
40implementations for Boolean nodes are:
41
42 - TGeoUnion - representing the Boolean union of two positioned shapes
43 - TGeoSubtraction - representing the Boolean subtraction of two positioned shapes
44 - TGeoIntersection - representing the Boolean intersection of two positioned shapes
45*/
46
48
49////////////////////////////////////////////////////////////////////////////////
50/// Constructor.
51
53 fSelected(0)
54{
55}
56
57////////////////////////////////////////////////////////////////////////////////
58/// Destructor.
59
61{
62}
63
64////////////////////////////////////////////////////////////////////////////////
65
67{
69/*
70 std::lock_guard<std::mutex> guard(fMutex);
71 if (tid >= fThreadSize) {
72 Error("GetThreadData", "Thread id=%d bigger than maximum declared thread number %d. \nUse TGeoManager::SetMaxThreads properly !!!",
73 tid, fThreadSize);
74 }
75 if (tid >= fThreadSize)
76 {
77 fThreadData.resize(tid + 1);
78 fThreadSize = tid + 1;
79 }
80 if (fThreadData[tid] == 0)
81 {
82 if (fThreadData[tid] == 0)
83 fThreadData[tid] = new ThreadData_t;
84 }
85*/
86 return *fThreadData[tid];
87}
88
89////////////////////////////////////////////////////////////////////////////////
90
92{
93 std::lock_guard<std::mutex> guard(fMutex);
94 std::vector<ThreadData_t*>::iterator i = fThreadData.begin();
95 while (i != fThreadData.end())
96 {
97 delete *i;
98 ++i;
99 }
100 fThreadData.clear();
101 fThreadSize = 0;
102}
103
104////////////////////////////////////////////////////////////////////////////////
105/// Create thread data for n threads max.
106
108{
109 std::lock_guard<std::mutex> guard(fMutex);
110 fThreadData.resize(nthreads);
111 fThreadSize = nthreads;
112 for (Int_t tid=0; tid<nthreads; tid++) {
113 if (fThreadData[tid] == 0) {
114 fThreadData[tid] = new ThreadData_t;
115 }
116 }
117 // Propagate to components
118 if (fLeft) fLeft->CreateThreadData(nthreads);
119 if (fRight) fRight->CreateThreadData(nthreads);
120}
121
122////////////////////////////////////////////////////////////////////////////////
123/// Set the selected branch.
124
126{
127 GetThreadData().fSelected = sel;
128}
129
130////////////////////////////////////////////////////////////////////////////////
131/// Default constructor
132
134{
135 fLeft = nullptr;
136 fRight = nullptr;
137 fLeftMat = nullptr;
138 fRightMat = nullptr;
139 fNpoints = 0;
140 fPoints = nullptr;
141 fThreadSize = 0;
143}
144
145////////////////////////////////////////////////////////////////////////////////
146/// Constructor called by TGeoCompositeShape providing 2 subexpressions for the 2 branches.
147
148TGeoBoolNode::TGeoBoolNode(const char *expr1, const char *expr2)
149{
150 fLeft = nullptr;
151 fRight = nullptr;
152 fLeftMat = nullptr;
153 fRightMat = nullptr;
154 fNpoints = 0;
155 fPoints = nullptr;
156 fThreadSize = 0;
158 if (!MakeBranch(expr1, kTRUE)) {
159 return;
160 }
161 if (!MakeBranch(expr2, kFALSE)) {
162 return;
163 }
164}
165
166////////////////////////////////////////////////////////////////////////////////
167/// Constructor providing left and right shapes and matrices (in the Boolean operation).
168
170{
171 fLeft = left;
172 fRight = right;
173 fLeftMat = lmat;
174 fNpoints = 0;
175 fPoints = nullptr;
176 fThreadSize = 0;
180 fRightMat = rmat;
183 if (!fLeft) {
184 Error("ctor", "left shape is NULL");
185 return;
186 }
187 if (!fRight) {
188 Error("ctor", "right shape is NULL");
189 return;
190 }
191}
192
193////////////////////////////////////////////////////////////////////////////////
194/// Destructor.
195/// --- deletion of components handled by TGeoManager class.
196
198{
199 if (fPoints) 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) delete shape;
279 return kFALSE;
280 }
281 if (left) {
282 fLeft = shape;
283 fLeftMat = mat;
284 } else {
285 fRight = shape;
286 fRightMat = mat;
287 }
288 return kTRUE;
289}
290
291////////////////////////////////////////////////////////////////////////////////
292/// Special schema for feeding the 3D buffers to the painter client.
293
295{
296 TVirtualViewer3D * viewer = gPad->GetViewer3D();
297 if (!viewer) return;
298
299 // Components of composite shape hierarchies for local frame viewers are painted
300 // in coordinate frame of the top level composite shape. So we force
301 // conversion to this. See TGeoPainter::PaintNode for loading of GLMatrix.
302 Bool_t localFrame = kFALSE; //viewer->PreferLocalFrame();
303
305 TGeoHMatrix mat;
306 mat = glmat; // keep a copy
307
308 // Now perform fetch and add of the two components buffers.
309 // Note we assume that composite shapes are always completely added
310 // so don't bother to get addDaughters flag from viewer->AddObject()
311
312 // Setup matrix and fetch/add the left component buffer
313 glmat->Multiply(fLeftMat);
314 //fLeft->Paint(option);
315 if (TGeoCompositeShape *left = dynamic_cast<TGeoCompositeShape *>(fLeft)) left->PaintComposite(option);
316 else {
317 const TBuffer3D & leftBuffer = fLeft->GetBuffer3D(TBuffer3D::kAll, localFrame);
318 viewer->AddObject(leftBuffer);
319 }
320
321 // Setup matrix and fetch/add the right component buffer
322 *glmat = &mat;
323 glmat->Multiply(fRightMat);
324 //fRight->Paint(option);
325 if (TGeoCompositeShape *right = dynamic_cast<TGeoCompositeShape *>(fRight)) right->PaintComposite(option);
326 else {
327 const TBuffer3D & rightBuffer = fRight->GetBuffer3D(TBuffer3D::kAll, localFrame);
328 viewer->AddObject(rightBuffer);
329 }
330
331 *glmat = &mat;
332}
333
334////////////////////////////////////////////////////////////////////////////////
335/// Register all matrices of the boolean node and descendents.
336
338{
341 if (fLeft->IsComposite()) ((TGeoCompositeShape*)fLeft)->GetBoolNode()->RegisterMatrices();
342 if (fRight->IsComposite()) ((TGeoCompositeShape*)fRight)->GetBoolNode()->RegisterMatrices();
343}
344
345////////////////////////////////////////////////////////////////////////////////
346/// Replace one of the matrices. Does not work with TGeoIdentity. Returns true
347/// if replacement was successful.
348
350{
351 if (mat==gGeoIdentity || newmat==gGeoIdentity) {
352 Error("ReplaceMatrix", "Matrices should not be gGeoIdentity. Use default matrix constructor to represent identities.");
353 return kFALSE;
354 }
355 if (!mat || !newmat) {
356 Error("ReplaceMatrix", "Matrices should not be null pointers.");
357 return kFALSE;
358 }
359 Bool_t replaced = kFALSE;
360 if (fLeftMat == mat) {
361 fLeftMat = newmat;
362 replaced = kTRUE;
363 }
364 if (fRightMat == mat) {
365 fRightMat = newmat;
366 replaced = kTRUE;
367 }
368 return replaced;
369}
370
371////////////////////////////////////////////////////////////////////////////////
372/// Save a primitive as a C++ statement(s) on output stream "out".
373
374void TGeoBoolNode::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
375{
376 fLeft->SavePrimitive(out,option);
377 fRight->SavePrimitive(out,option);
378 if (!fLeftMat->IsIdentity()) {
380 fLeftMat->SavePrimitive(out,option);
381 }
382 if (!fRightMat->IsIdentity()) {
384 fRightMat->SavePrimitive(out,option);
385 }
386}
387
388////////////////////////////////////////////////////////////////////////////////
389/// Fill buffer with shape vertices.
390
392{
393 TGeoBoolNode *bn = (TGeoBoolNode*)this;
394 Int_t npoints = bn->GetNpoints();
395 memcpy(points, fPoints, 3*npoints*sizeof(Double_t));
396}
397
398////////////////////////////////////////////////////////////////////////////////
399/// Fill buffer with shape vertices.
400
402{
403 TGeoBoolNode *bn = (TGeoBoolNode*)this;
404 Int_t npoints = bn->GetNpoints();
405 for (Int_t i=0; i<3*npoints; i++) points[i] = fPoints[i];
406}
407
408////////////////////////////////////////////////////////////////////////////////
409/// Register size of this 3D object
410
412{
413 fLeft->Sizeof3D();
414 fRight->Sizeof3D();
415}
416
418
419////////////////////////////////////////////////////////////////////////////////
420/// Make a clone of this. Pointers are preserved.
421
423{
424 return new TGeoUnion(fLeft, fRight, fLeftMat, fRightMat);
425}
426
427////////////////////////////////////////////////////////////////////////////////
428/// Paint method.
429
431{
432 TVirtualViewer3D *viewer = gPad->GetViewer3D();
433
434 if (!viewer) {
435 Error("Paint", "gPad->GetViewer3D() returned 0, cannot work with composite!\n");
436 return;
437 }
438
440
441 TGeoBoolNode::Paint(option);
442}
443
444////////////////////////////////////////////////////////////////////////////////
445/// Default constructor
446
448{
449}
450
451////////////////////////////////////////////////////////////////////////////////
452/// Constructor
453
454TGeoUnion::TGeoUnion(const char *expr1, const char *expr2)
455 :TGeoBoolNode(expr1, expr2)
456{
457}
458
459////////////////////////////////////////////////////////////////////////////////
460/// Constructor providing pointers to components
461
463 :TGeoBoolNode(left,right,lmat,rmat)
464{
466 Fatal("TGeoUnion", "Unions with a half-space (%s + %s) not allowed", left->GetName(), right->GetName());
467 }
468}
469
470////////////////////////////////////////////////////////////////////////////////
471/// Destructor
472/// --- deletion of components handled by TGeoManager class.
473
475{
476}
477
478////////////////////////////////////////////////////////////////////////////////
479/// Compute bounding box corresponding to a union of two shapes.
480
482{
483 if (((TGeoBBox*)fLeft)->IsNullBox()) fLeft->ComputeBBox();
484 if (((TGeoBBox*)fRight)->IsNullBox()) fRight->ComputeBBox();
485 Double_t vert[48];
486 Double_t pt[3];
487 Int_t i;
488 Double_t xmin, xmax, ymin, ymax, zmin, zmax;
489 xmin = ymin = zmin = TGeoShape::Big();
490 xmax = ymax = zmax = -TGeoShape::Big();
491 ((TGeoBBox*)fLeft)->SetBoxPoints(&vert[0]);
492 ((TGeoBBox*)fRight)->SetBoxPoints(&vert[24]);
493 for (i=0; i<8; i++) {
494 fLeftMat->LocalToMaster(&vert[3*i], &pt[0]);
495 if (pt[0]<xmin) xmin=pt[0];
496 if (pt[0]>xmax) xmax=pt[0];
497 if (pt[1]<ymin) ymin=pt[1];
498 if (pt[1]>ymax) ymax=pt[1];
499 if (pt[2]<zmin) zmin=pt[2];
500 if (pt[2]>zmax) zmax=pt[2];
501 }
502 for (i=8; i<16; i++) {
503 fRightMat->LocalToMaster(&vert[3*i], &pt[0]);
504 if (pt[0]<xmin) xmin=pt[0];
505 if (pt[0]>xmax) xmax=pt[0];
506 if (pt[1]<ymin) ymin=pt[1];
507 if (pt[1]>ymax) ymax=pt[1];
508 if (pt[2]<zmin) zmin=pt[2];
509 if (pt[2]>zmax) zmax=pt[2];
510 }
511 dx = 0.5*(xmax-xmin);
512 origin[0] = 0.5*(xmin+xmax);
513 dy = 0.5*(ymax-ymin);
514 origin[1] = 0.5*(ymin+ymax);
515 dz = 0.5*(zmax-zmin);
516 origin[2] = 0.5*(zmin+zmax);
517}
518
519////////////////////////////////////////////////////////////////////////////////
520/// Find if a union of two shapes contains a given point
521
523{
524 Double_t local[3];
525 fLeftMat->MasterToLocal(point, &local[0]);
526 Bool_t inside = fLeft->Contains(&local[0]);
527 if (inside) return kTRUE;
528 fRightMat->MasterToLocal(point, &local[0]);
529 inside = fRight->Contains(&local[0]);
530 return inside;
531}
532
533////////////////////////////////////////////////////////////////////////////////
534/// Normal computation in POINT. The orientation is chosen so that DIR.dot.NORM>0.
535
536void TGeoUnion::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
537{
539 norm[0] = norm[1] = 0.;
540 norm[2] = 1.;
541 Double_t local[3];
542 Double_t ldir[3], lnorm[3];
543 if (td.fSelected == 1) {
544 fLeftMat->MasterToLocal(point, local);
545 fLeftMat->MasterToLocalVect(dir, ldir);
546 fLeft->ComputeNormal(local,ldir,lnorm);
547 fLeftMat->LocalToMasterVect(lnorm, norm);
548 return;
549 }
550 if (td.fSelected == 2) {
551 fRightMat->MasterToLocal(point, local);
552 fRightMat->MasterToLocalVect(dir, ldir);
553 fRight->ComputeNormal(local,ldir,lnorm);
554 fRightMat->LocalToMasterVect(lnorm, norm);
555 return;
556 }
557 fLeftMat->MasterToLocal(point, local);
558 if (fLeft->Contains(local)) {
559 fLeftMat->MasterToLocalVect(dir, ldir);
560 fLeft->ComputeNormal(local,ldir,lnorm);
561 fLeftMat->LocalToMasterVect(lnorm, norm);
562 return;
563 }
564 fRightMat->MasterToLocal(point, local);
565 if (fRight->Contains(local)) {
566 fRightMat->MasterToLocalVect(dir, ldir);
567 fRight->ComputeNormal(local,ldir,lnorm);
568 fRightMat->LocalToMasterVect(lnorm, norm);
569 return;
570 }
571 // Propagate forward/backward to see which of the components is intersected first
572 local[0] = point[0] + 1E-5*dir[0];
573 local[1] = point[1] + 1E-5*dir[1];
574 local[2] = point[2] + 1E-5*dir[2];
575
576 if (!Contains(local)) {
577 local[0] = point[0] - 1E-5*dir[0];
578 local[1] = point[1] - 1E-5*dir[1];
579 local[2] = point[2] - 1E-5*dir[2];
580 if (!Contains(local)) return;
581 }
582 ComputeNormal(local,dir,norm);
583}
584
585////////////////////////////////////////////////////////////////////////////////
586/// Compute minimum distance to shape vertices.
587
589{
590 return 9999;
591}
592
593////////////////////////////////////////////////////////////////////////////////
594/// Computes distance from a given point inside the shape to its boundary.
595
597 Double_t step, Double_t *safe) const
598{
599 if (iact<3 && safe) {
600 // compute safe distance
601 *safe = Safety(point,kTRUE);
602 if (iact==0) return TGeoShape::Big();
603 if (iact==1 && step<*safe) return TGeoShape::Big();
604 }
605
606 Double_t local[3], local1[3], master[3], ldir[3], rdir[3], pushed[3];
607 memcpy(master, point, 3*sizeof(Double_t));
608 Int_t i;
609 TGeoBoolNode *node = (TGeoBoolNode*)this;
610 Double_t d1=0., d2=0., snxt=0., eps=0.;
611 fLeftMat->MasterToLocalVect(dir, ldir);
612 fRightMat->MasterToLocalVect(dir, rdir);
613 fLeftMat->MasterToLocal(point, local);
614 Bool_t inside1 = fLeft->Contains(local);
615 if (inside1) d1 = fLeft->DistFromInside(local, ldir, 3);
616 else memcpy(local1, local, 3*sizeof(Double_t));
617 fRightMat->MasterToLocal(point, local);
618 Bool_t inside2 = fRight->Contains(local);
619 if (inside2) d2 = fRight->DistFromInside(local, rdir, 3);
620 if (!(inside1 | inside2)) {
621 // This is a pathological case when the point is on the boundary
622 d1 = fLeft->DistFromOutside(local1, ldir, 3);
623 if (d1<1.E-3) {
624 eps = d1+TGeoShape::Tolerance();
625 for (i=0; i<3; i++) local1[i] += eps*ldir[i];
626 inside1 = kTRUE;
627 d1 = fLeft->DistFromInside(local1, ldir, 3);
628 d1 += eps;
629 } else {
630 d2 = fRight->DistFromOutside(local, rdir, 3);
631 if (d2<1.E-3) {
632 eps = d2+TGeoShape::Tolerance();
633 for (i=0; i<3; i++) local[i] += eps*rdir[i];
634 inside2 = kTRUE;
635 d2 = fRight->DistFromInside(local, rdir, 3);
636 d2 += eps;
637 }
638 }
639 }
640 while (inside1 || inside2) {
641 if (inside1 && inside2) {
642 if (d1<d2) {
643 snxt += d1;
644 node->SetSelected(1);
645 // propagate to exit of left shape
646 inside1 = kFALSE;
647 for (i=0; i<3; i++) master[i] += d1*dir[i];
648 // check if propagated point is in right shape
649 fRightMat->MasterToLocal(master, local);
650 inside2 = fRight->Contains(local);
651 if (!inside2) return snxt;
652 d2 = fRight->DistFromInside(local, rdir, 3);
653 if (d2 < TGeoShape::Tolerance()) return snxt;
654 } else {
655 snxt += d2;
656 node->SetSelected(2);
657 // propagate to exit of right shape
658 inside2 = kFALSE;
659 for (i=0; i<3; i++) master[i] += d2*dir[i];
660 // check if propagated point is in left shape
661 fLeftMat->MasterToLocal(master, local);
662 inside1 = fLeft->Contains(local);
663 if (!inside1) return snxt;
664 d1 = fLeft->DistFromInside(local, ldir, 3);
665 if (d1 < TGeoShape::Tolerance()) return snxt;
666 }
667 }
668 if (inside1) {
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 pushed[i] = master[i]+(1.+d1)*TGeoShape::Tolerance()*dir[i];
676 }
677 // check if propagated point is in right shape
678 fRightMat->MasterToLocal(pushed, local);
679 inside2 = fRight->Contains(local);
680 if (!inside2) return snxt;
681 d2 = fRight->DistFromInside(local, rdir, 3);
682 if (d2 < TGeoShape::Tolerance()) return snxt;
683 d2 += (1.+d1)*TGeoShape::Tolerance();
684 }
685 if (inside2) {
686 snxt += d2;
687 node->SetSelected(2);
688 // propagate to exit of right shape
689 inside2 = kFALSE;
690 for (i=0; i<3; i++) {
691 master[i] += d2*dir[i];
692 pushed[i] = master[i]+(1.+d2)*TGeoShape::Tolerance()*dir[i];
693 }
694 // check if propagated point is in left shape
695 fLeftMat->MasterToLocal(pushed, local);
696 inside1 = fLeft->Contains(local);
697 if (!inside1) return snxt;
698 d1 = fLeft->DistFromInside(local, ldir, 3);
699 if (d1 < TGeoShape::Tolerance()) return snxt;
700 d1 += (1.+d2)*TGeoShape::Tolerance();
701 }
702 }
703 return snxt;
704}
705
706////////////////////////////////////////////////////////////////////////////////
707/// Compute distance from a given outside point to the shape.
708
710 Double_t step, Double_t *safe) const
711{
712 if (iact<3 && safe) {
713 // compute safe distance
714 *safe = Safety(point,kFALSE);
715 if (iact==0) return TGeoShape::Big();
716 if (iact==1 && step<*safe) return TGeoShape::Big();
717 }
718 TGeoBoolNode *node = (TGeoBoolNode*)this;
719 Double_t local[3], ldir[3], rdir[3];
720 Double_t d1, d2, snxt;
721 fLeftMat->MasterToLocal(point, &local[0]);
722 fLeftMat->MasterToLocalVect(dir, &ldir[0]);
723 fRightMat->MasterToLocalVect(dir, &rdir[0]);
724 d1 = fLeft->DistFromOutside(&local[0], &ldir[0], iact, step, safe);
725 fRightMat->MasterToLocal(point, &local[0]);
726 d2 = fRight->DistFromOutside(&local[0], &rdir[0], iact, step, safe);
727 if (d1<d2) {
728 snxt = d1;
729 node->SetSelected(1);
730 } else {
731 snxt = d2;
732 node->SetSelected(2);
733 }
734 return snxt;
735}
736
737////////////////////////////////////////////////////////////////////////////////
738/// Returns number of vertices for the composite shape described by this union.
739
741{
742 Int_t itot=0;
743 Double_t point[3];
744 Double_t tolerance = TGeoShape::Tolerance();
745 if (fNpoints) return fNpoints;
746 // Local points for the left shape
747 Int_t nleft = fLeft->GetNmeshVertices();
748 Double_t *points1 = new Double_t[3*nleft];
749 fLeft->SetPoints(points1);
750 // Local points for the right shape
751 Int_t nright = fRight->GetNmeshVertices();
752 Double_t *points2 = new Double_t[3*nright];
753 fRight->SetPoints(points2);
754 Double_t *points = new Double_t[3*(nleft+nright)];
755 for (Int_t i=0; i<nleft; i++) {
756 if (TMath::Abs(points1[3*i])<tolerance && TMath::Abs(points1[3*i+1])<tolerance) continue;
757 fLeftMat->LocalToMaster(&points1[3*i], &points[3*itot]);
758 fRightMat->MasterToLocal(&points[3*itot], point);
759 if (!fRight->Contains(point)) itot++;
760 }
761 for (Int_t i=0; i<nright; i++) {
762 if (TMath::Abs(points2[3*i])<tolerance && TMath::Abs(points2[3*i+1])<tolerance) continue;
763 fRightMat->LocalToMaster(&points2[3*i], &points[3*itot]);
764 fLeftMat->MasterToLocal(&points[3*itot], point);
765 if (!fLeft->Contains(point)) itot++;
766 }
767
768 AssignPoints(itot, points);
769
770 delete [] points1;
771 delete [] points2;
772 delete [] points;
773 return fNpoints;
774}
775
776////////////////////////////////////////////////////////////////////////////////
777/// Compute safety distance for a union node;
778
780{
781 Double_t local1[3], local2[3];
782 fLeftMat->MasterToLocal(point,local1);
783 Bool_t in1 = fLeft->Contains(local1);
784 fRightMat->MasterToLocal(point,local2);
785 Bool_t in2 = fRight->Contains(local2);
786 Bool_t intrue = in1 | in2;
787 if (intrue^in) return 0.0;
788 Double_t saf1 = fLeft->Safety(local1, in1);
789 Double_t saf2 = fRight->Safety(local2, in2);
790 if (in1 && in2) return TMath::Min(saf1, saf2);
791 if (in1) return saf1;
792 if (in2) return saf2;
793 return TMath::Min(saf1,saf2);
794}
795
796////////////////////////////////////////////////////////////////////////////////
797/// Save a primitive as a C++ statement(s) on output stream "out".
798
799void TGeoUnion::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
800{
801 TGeoBoolNode::SavePrimitive(out,option);
802 out << " pBoolNode = new TGeoUnion(";
803 out << fLeft->GetPointerName() << ",";
804 out << fRight->GetPointerName() << ",";
805 if (!fLeftMat->IsIdentity()) out << fLeftMat->GetPointerName() << ",";
806 else out << "0,";
807 if (!fRightMat->IsIdentity()) out << fRightMat->GetPointerName() << ");" << std::endl;
808 else out << "0);" << std::endl;
809}
810
811////////////////////////////////////////////////////////////////////////////////
812/// Register 3D size of this shape.
813
815{
817}
818
819
821
822////////////////////////////////////////////////////////////////////////////////
823/// Make a clone of this. Pointers are preserved.
824
826{
828}
829
830////////////////////////////////////////////////////////////////////////////////
831/// Paint method.
832
834{
835 TVirtualViewer3D *viewer = gPad->GetViewer3D();
836
837 if (!viewer) {
838 Error("Paint", "gPad->GetViewer3D() returned 0, cannot work with composite!\n");
839 return;
840 }
841
843
844 TGeoBoolNode::Paint(option);
845}
846
847////////////////////////////////////////////////////////////////////////////////
848/// Default constructor
849
851{
852}
853
854////////////////////////////////////////////////////////////////////////////////
855/// Constructor
856
857TGeoSubtraction::TGeoSubtraction(const char *expr1, const char *expr2)
858 :TGeoBoolNode(expr1, expr2)
859{
860}
861
862////////////////////////////////////////////////////////////////////////////////
863/// Constructor providing pointers to components
864
866 :TGeoBoolNode(left,right,lmat,rmat)
867{
869 Fatal("TGeoSubstraction", "Subtractions from a half-space (%s) not allowed", left->GetName());
870 }
871}
872
873////////////////////////////////////////////////////////////////////////////////
874/// Destructor
875/// --- deletion of components handled by TGeoManager class.
876
878{
879}
880
881////////////////////////////////////////////////////////////////////////////////
882/// Compute bounding box corresponding to a subtraction of two shapes.
883
885{
887 if (box->IsNullBox()) fLeft->ComputeBBox();
888 Double_t vert[24];
889 Double_t pt[3];
890 Int_t i;
891 Double_t xmin, xmax, ymin, ymax, zmin, zmax;
892 xmin = ymin = zmin = TGeoShape::Big();
893 xmax = ymax = zmax = -TGeoShape::Big();
894 box->SetBoxPoints(&vert[0]);
895 for (i=0; i<8; i++) {
896 fLeftMat->LocalToMaster(&vert[3*i], &pt[0]);
897 if (pt[0]<xmin) xmin=pt[0];
898 if (pt[0]>xmax) xmax=pt[0];
899 if (pt[1]<ymin) ymin=pt[1];
900 if (pt[1]>ymax) ymax=pt[1];
901 if (pt[2]<zmin) zmin=pt[2];
902 if (pt[2]>zmax) zmax=pt[2];
903 }
904 dx = 0.5*(xmax-xmin);
905 origin[0] = 0.5*(xmin+xmax);
906 dy = 0.5*(ymax-ymin);
907 origin[1] = 0.5*(ymin+ymax);
908 dz = 0.5*(zmax-zmin);
909 origin[2] = 0.5*(zmin+zmax);
910}
911
912////////////////////////////////////////////////////////////////////////////////
913/// Normal computation in POINT. The orientation is chosen so that DIR.dot.NORM>0.
914
915void TGeoSubtraction::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
916{
918 norm[0] = norm[1] = 0.;
919 norm[2] = 1.;
920 Double_t local[3], ldir[3], lnorm[3];
921 if (td.fSelected == 1) {
922 fLeftMat->MasterToLocal(point, local);
923 fLeftMat->MasterToLocalVect(dir, ldir);
924 fLeft->ComputeNormal(local,ldir,lnorm);
925 fLeftMat->LocalToMasterVect(lnorm, norm);
926 return;
927 }
928 if (td.fSelected == 2) {
929 fRightMat->MasterToLocal(point, local);
930 fRightMat->MasterToLocalVect(dir, ldir);
931 fRight->ComputeNormal(local,ldir,lnorm);
932 fRightMat->LocalToMasterVect(lnorm, norm);
933 return;
934 }
935 fRightMat->MasterToLocal(point,local);
936 if (fRight->Contains(local)) {
937 fRightMat->MasterToLocalVect(dir,ldir);
938 fRight->ComputeNormal(local,ldir, lnorm);
939 fRightMat->LocalToMasterVect(lnorm,norm);
940 return;
941 }
942 fLeftMat->MasterToLocal(point,local);
943 if (!fLeft->Contains(local)) {
944 fLeftMat->MasterToLocalVect(dir,ldir);
945 fLeft->ComputeNormal(local,ldir, lnorm);
946 fLeftMat->LocalToMasterVect(lnorm,norm);
947 return;
948 }
949 // point is inside left shape, but not inside the right
950 local[0] = point[0]+1E-5*dir[0];
951 local[1] = point[1]+1E-5*dir[1];
952 local[2] = point[2]+1E-5*dir[2];
953 if (Contains(local)) {
954 local[0] = point[0]-1E-5*dir[0];
955 local[1] = point[1]-1E-5*dir[1];
956 local[2] = point[2]-1E-5*dir[2];
957 if (Contains(local)) return;
958 }
959 ComputeNormal(local,dir,norm);
960}
961
962////////////////////////////////////////////////////////////////////////////////
963/// Find if a subtraction of two shapes contains a given point
964
966{
967 Double_t local[3];
968 fLeftMat->MasterToLocal(point, &local[0]);
969 Bool_t inside = fLeft->Contains(&local[0]);
970 if (!inside) return kFALSE;
971 fRightMat->MasterToLocal(point, &local[0]);
972 inside = !fRight->Contains(&local[0]);
973 return inside;
974}
975
976////////////////////////////////////////////////////////////////////////////////
977/// Compute minimum distance to shape vertices
978
980{
981 return 9999;
982}
983
984////////////////////////////////////////////////////////////////////////////////
985/// Compute distance from a given point inside to the shape boundary.
986
988 Double_t step, Double_t *safe) const
989{
990 if (iact<3 && safe) {
991 // compute safe distance
992 *safe = Safety(point,kTRUE);
993 if (iact==0) return TGeoShape::Big();
994 if (iact==1 && step<*safe) return TGeoShape::Big();
995 }
996 TGeoBoolNode *node = (TGeoBoolNode*)this;
997 Double_t local[3], ldir[3], rdir[3];
998 Double_t d1, d2, snxt=0.;
999 fLeftMat->MasterToLocal(point, &local[0]);
1000 fLeftMat->MasterToLocalVect(dir, &ldir[0]);
1001 fRightMat->MasterToLocalVect(dir, &rdir[0]);
1002 d1 = fLeft->DistFromInside(&local[0], &ldir[0], iact, step, safe);
1003 fRightMat->MasterToLocal(point, &local[0]);
1004 d2 = fRight->DistFromOutside(&local[0], &rdir[0], iact, step, safe);
1005 if (d1<d2) {
1006 snxt = d1;
1007 node->SetSelected(1);
1008 } else {
1009 snxt = d2;
1010 node->SetSelected(2);
1011 }
1012 return snxt;
1013}
1014
1015////////////////////////////////////////////////////////////////////////////////
1016/// Compute distance from a given point outside to the shape.
1017
1019 Double_t step, Double_t *safe) const
1020{
1021 if (iact<3 && safe) {
1022 // compute safe distance
1023 *safe = Safety(point,kFALSE);
1024 if (iact==0) return TGeoShape::Big();
1025 if (iact==1 && step<*safe) return TGeoShape::Big();
1026 }
1027 TGeoBoolNode *node = (TGeoBoolNode*)this;
1028 Double_t local[3], master[3], ldir[3], rdir[3];
1029 memcpy(&master[0], point, 3*sizeof(Double_t));
1030 Int_t i;
1031 Double_t d1, d2, snxt=0.;
1032 fRightMat->MasterToLocal(point, &local[0]);
1033 fLeftMat->MasterToLocalVect(dir, &ldir[0]);
1034 fRightMat->MasterToLocalVect(dir, &rdir[0]);
1035 // check if inside '-'
1036 Bool_t inside = fRight->Contains(&local[0]);
1037 Double_t epsil = 0.;
1038 while (1) {
1039 if (inside) {
1040 // propagate to outside of '-'
1041 node->SetSelected(2);
1042 d1 = fRight->DistFromInside(&local[0], &rdir[0], iact, step, safe);
1043 snxt += d1+epsil;
1044 for (i=0; i<3; i++) master[i] += (d1+1E-8)*dir[i];
1045 epsil = 1.E-8;
1046 // now master outside '-'; check if inside '+'
1047 fLeftMat->MasterToLocal(&master[0], &local[0]);
1048 if (fLeft->Contains(&local[0])) return snxt;
1049 }
1050 // master outside '-' and outside '+' ; find distances to both
1051 node->SetSelected(1);
1052 fLeftMat->MasterToLocal(&master[0], &local[0]);
1053 d2 = fLeft->DistFromOutside(&local[0], &ldir[0], iact, step, safe);
1054 if (d2>1E20) return TGeoShape::Big();
1055
1056 fRightMat->MasterToLocal(&master[0], &local[0]);
1057 d1 = fRight->DistFromOutside(&local[0], &rdir[0], iact, step, safe);
1058 if (d2<d1-TGeoShape::Tolerance()) {
1059 snxt += d2+epsil;
1060 return snxt;
1061 }
1062 // propagate to '-'
1063 snxt += d1+epsil;
1064 for (i=0; i<3; i++) master[i] += (d1+1E-8)*dir[i];
1065 epsil = 1.E-8;
1066 // now inside '-' and not inside '+'
1067 fRightMat->MasterToLocal(&master[0], &local[0]);
1068 inside = kTRUE;
1069 }
1070}
1071
1072////////////////////////////////////////////////////////////////////////////////
1073/// Returns number of vertices for the composite shape described by this subtraction.
1074
1076{
1077 Int_t itot=0;
1078 Double_t point[3];
1079 Double_t tolerance = TGeoShape::Tolerance();
1080 if (fNpoints) return fNpoints;
1081 Int_t nleft = fLeft->GetNmeshVertices();
1082 Int_t nright = fRight->GetNmeshVertices();
1083 Double_t *points = new Double_t[3*(nleft+nright)];
1084 Double_t *points1 = new Double_t[3*nleft];
1085 fLeft->SetPoints(points1);
1086 for (Int_t i=0; i<nleft; i++) {
1087 if (TMath::Abs(points1[3*i])<tolerance && TMath::Abs(points1[3*i+1])<tolerance) continue;
1088 fLeftMat->LocalToMaster(&points1[3*i], &points[3*itot]);
1089 fRightMat->MasterToLocal(&points[3*itot], point);
1090 if (!fRight->Contains(point)) itot++;
1091 }
1092 Double_t *points2 = new Double_t[3*nright];
1093 fRight->SetPoints(points2);
1094 for (Int_t i=0; i<nright; i++) {
1095 if (TMath::Abs(points2[3*i])<tolerance && TMath::Abs(points2[3*i+1])<tolerance) continue;
1096 fRightMat->LocalToMaster(&points2[3*i], &points[3*itot]);
1097 fLeftMat->MasterToLocal(&points[3*itot], point);
1098 if (fLeft->Contains(point)) itot++;
1099 }
1100
1101 AssignPoints(itot, points);
1102
1103 delete [] points1;
1104 delete [] points2;
1105 delete [] points;
1106 return fNpoints;
1107}
1108
1109////////////////////////////////////////////////////////////////////////////////
1110/// Compute safety distance for a union node;
1111
1113{
1114 Double_t local1[3], local2[3];
1115 fLeftMat->MasterToLocal(point,local1);
1116 Bool_t in1 = fLeft->Contains(local1);
1117 fRightMat->MasterToLocal(point,local2);
1118 Bool_t in2 = fRight->Contains(local2);
1119 Bool_t intrue = in1 && (!in2);
1120 if (in^intrue) return 0.0;
1121 Double_t saf1 = fLeft->Safety(local1, in1);
1122 Double_t saf2 = fRight->Safety(local2, in2);
1123 if (in1 && in2) return saf2;
1124 if (in1) return TMath::Min(saf1,saf2);
1125 if (in2) return TMath::Max(saf1,saf2);
1126 return saf1;
1127}
1128
1129////////////////////////////////////////////////////////////////////////////////
1130/// Save a primitive as a C++ statement(s) on output stream "out".
1131
1132void TGeoSubtraction::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
1133{
1134 TGeoBoolNode::SavePrimitive(out,option);
1135 out << " pBoolNode = new TGeoSubtraction(";
1136 out << fLeft->GetPointerName() << ",";
1137 out << fRight->GetPointerName() << ",";
1138 if (!fLeftMat->IsIdentity()) out << fLeftMat->GetPointerName() << ",";
1139 else out << "0,";
1140 if (!fRightMat->IsIdentity()) out << fRightMat->GetPointerName() << ");" << std::endl;
1141 else out << "0);" << std::endl;
1142}
1143
1144////////////////////////////////////////////////////////////////////////////////
1145/// Register 3D size of this shape.
1146
1148{
1150}
1151
1153
1154////////////////////////////////////////////////////////////////////////////////
1155/// Make a clone of this. Pointers are preserved.
1156
1158{
1160}
1161
1162////////////////////////////////////////////////////////////////////////////////
1163/// Paint method.
1164
1166{
1167 TVirtualViewer3D *viewer = gPad->GetViewer3D();
1168
1169 if (!viewer) {
1170 Error("Paint", "gPad->GetViewer3D() returned 0, cannot work with composite!\n");
1171 return;
1172 }
1173
1175
1176 TGeoBoolNode::Paint(option);
1177}
1178
1179////////////////////////////////////////////////////////////////////////////////
1180/// Default constructor
1181
1183{
1184}
1185
1186////////////////////////////////////////////////////////////////////////////////
1187/// Constructor
1188
1189TGeoIntersection::TGeoIntersection(const char *expr1, const char *expr2)
1190 :TGeoBoolNode(expr1, expr2)
1191{
1192}
1193
1194////////////////////////////////////////////////////////////////////////////////
1195/// Constructor providing pointers to components
1196
1198 :TGeoBoolNode(left,right,lmat,rmat)
1199{
1202 if (hs1 && hs2) Fatal("ctor", "cannot intersect two half-spaces: %s * %s", left->GetName(), right->GetName());
1203}
1204
1205////////////////////////////////////////////////////////////////////////////////
1206/// Destructor
1207/// --- deletion of components handled by TGeoManager class.
1208
1210{
1211}
1212
1213////////////////////////////////////////////////////////////////////////////////
1214/// Compute bounding box corresponding to a intersection of two shapes.
1215
1217{
1220 Double_t vert[48];
1221 Double_t pt[3];
1222 Int_t i;
1223 Double_t xmin1, xmax1, ymin1, ymax1, zmin1, zmax1;
1224 Double_t xmin2, xmax2, ymin2, ymax2, zmin2, zmax2;
1225 xmin1 = ymin1 = zmin1 = xmin2 = ymin2 = zmin2 = TGeoShape::Big();
1226 xmax1 = ymax1 = zmax1 = xmax2 = ymax2 = zmax2 = -TGeoShape::Big();
1227 if (!hs1) {
1228 if (((TGeoBBox*)fLeft)->IsNullBox()) fLeft->ComputeBBox();
1229 ((TGeoBBox*)fLeft)->SetBoxPoints(&vert[0]);
1230 for (i=0; i<8; i++) {
1231 fLeftMat->LocalToMaster(&vert[3*i], &pt[0]);
1232 if (pt[0]<xmin1) xmin1=pt[0];
1233 if (pt[0]>xmax1) xmax1=pt[0];
1234 if (pt[1]<ymin1) ymin1=pt[1];
1235 if (pt[1]>ymax1) ymax1=pt[1];
1236 if (pt[2]<zmin1) zmin1=pt[2];
1237 if (pt[2]>zmax1) zmax1=pt[2];
1238 }
1239 }
1240 if (!hs2) {
1241 if (((TGeoBBox*)fRight)->IsNullBox()) fRight->ComputeBBox();
1242 ((TGeoBBox*)fRight)->SetBoxPoints(&vert[24]);
1243 for (i=8; i<16; i++) {
1244 fRightMat->LocalToMaster(&vert[3*i], &pt[0]);
1245 if (pt[0]<xmin2) xmin2=pt[0];
1246 if (pt[0]>xmax2) xmax2=pt[0];
1247 if (pt[1]<ymin2) ymin2=pt[1];
1248 if (pt[1]>ymax2) ymax2=pt[1];
1249 if (pt[2]<zmin2) zmin2=pt[2];
1250 if (pt[2]>zmax2) zmax2=pt[2];
1251 }
1252 }
1253 if (hs1) {
1254 dx = 0.5*(xmax2-xmin2);
1255 origin[0] = 0.5*(xmax2+xmin2);
1256 dy = 0.5*(ymax2-ymin2);
1257 origin[1] = 0.5*(ymax2+ymin2);
1258 dz = 0.5*(zmax2-zmin2);
1259 origin[2] = 0.5*(zmax2+zmin2);
1260 return;
1261 }
1262 if (hs2) {
1263 dx = 0.5*(xmax1-xmin1);
1264 origin[0] = 0.5*(xmax1+xmin1);
1265 dy = 0.5*(ymax1-ymin1);
1266 origin[1] = 0.5*(ymax1+ymin1);
1267 dz = 0.5*(zmax1-zmin1);
1268 origin[2] = 0.5*(zmax1+zmin1);
1269 return;
1270 }
1271 Double_t sort[4];
1272 Int_t isort[4];
1273 sort[0] = xmin1;
1274 sort[1] = xmax1;
1275 sort[2] = xmin2;
1276 sort[3] = xmax2;
1277 TMath::Sort(4, &sort[0], &isort[0], kFALSE);
1278 if (isort[1]%2) {
1279 Warning("ComputeBBox", "shapes %s and %s do not intersect", fLeft->GetName(), fRight->GetName());
1280 dx = dy = dz = 0;
1281 memset(origin, 0, 3*sizeof(Double_t));
1282 return;
1283 }
1284 dx = 0.5*(sort[isort[2]]-sort[isort[1]]);
1285 origin[0] = 0.5*(sort[isort[1]]+sort[isort[2]]);
1286 sort[0] = ymin1;
1287 sort[1] = ymax1;
1288 sort[2] = ymin2;
1289 sort[3] = ymax2;
1290 TMath::Sort(4, &sort[0], &isort[0], kFALSE);
1291 if (isort[1]%2) {
1292 Warning("ComputeBBox", "shapes %s and %s do not intersect", fLeft->GetName(), fRight->GetName());
1293 dx = dy = dz = 0;
1294 memset(origin, 0, 3*sizeof(Double_t));
1295 return;
1296 }
1297 dy = 0.5*(sort[isort[2]]-sort[isort[1]]);
1298 origin[1] = 0.5*(sort[isort[1]]+sort[isort[2]]);
1299 sort[0] = zmin1;
1300 sort[1] = zmax1;
1301 sort[2] = zmin2;
1302 sort[3] = zmax2;
1303 TMath::Sort(4, &sort[0], &isort[0], kFALSE);
1304 if (isort[1]%2) {
1305 Warning("ComputeBBox", "shapes %s and %s do not intersect", fLeft->GetName(), fRight->GetName());
1306 dx = dy = dz = 0;
1307 memset(origin, 0, 3*sizeof(Double_t));
1308 return;
1309 }
1310 dz = 0.5*(sort[isort[2]]-sort[isort[1]]);
1311 origin[2] = 0.5*(sort[isort[1]]+sort[isort[2]]);
1312}
1313
1314////////////////////////////////////////////////////////////////////////////////
1315/// Normal computation in POINT. The orientation is chosen so that DIR.dot.NORM>0.
1316
1317void TGeoIntersection::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
1318{
1320 Double_t local[3], ldir[3], lnorm[3];
1321 norm[0] = norm[1] = 0.;
1322 norm[2] = 1.;
1323 if (td.fSelected == 1) {
1324 fLeftMat->MasterToLocal(point, local);
1325 fLeftMat->MasterToLocalVect(dir, ldir);
1326 fLeft->ComputeNormal(local,ldir,lnorm);
1327 fLeftMat->LocalToMasterVect(lnorm, norm);
1328 return;
1329 }
1330 if (td.fSelected == 2) {
1331 fRightMat->MasterToLocal(point, local);
1332 fRightMat->MasterToLocalVect(dir, ldir);
1333 fRight->ComputeNormal(local,ldir,lnorm);
1334 fRightMat->LocalToMasterVect(lnorm, norm);
1335 return;
1336 }
1337 fLeftMat->MasterToLocal(point,local);
1338 if (!fLeft->Contains(local)) {
1339 fLeftMat->MasterToLocalVect(dir,ldir);
1340 fLeft->ComputeNormal(local,ldir,lnorm);
1341 fLeftMat->LocalToMasterVect(lnorm,norm);
1342 return;
1343 }
1344 fRightMat->MasterToLocal(point,local);
1345 if (!fRight->Contains(local)) {
1346 fRightMat->MasterToLocalVect(dir,ldir);
1347 fRight->ComputeNormal(local,ldir,lnorm);
1348 fRightMat->LocalToMasterVect(lnorm,norm);
1349 return;
1350 }
1351 // point is inside intersection.
1352 local[0] = point[0] + 1E-5*dir[0];
1353 local[1] = point[1] + 1E-5*dir[1];
1354 local[2] = point[2] + 1E-5*dir[2];
1355 if (Contains(local)) {
1356 local[0] = point[0] - 1E-5*dir[0];
1357 local[1] = point[1] - 1E-5*dir[1];
1358 local[2] = point[2] - 1E-5*dir[2];
1359 if (Contains(local)) return;
1360 }
1361 ComputeNormal(local,dir,norm);
1362}
1363
1364////////////////////////////////////////////////////////////////////////////////
1365/// Find if a intersection of two shapes contains a given point
1366
1368{
1369 Double_t local[3];
1370 fLeftMat->MasterToLocal(point, &local[0]);
1371 Bool_t inside = fLeft->Contains(&local[0]);
1372 if (!inside) return kFALSE;
1373 fRightMat->MasterToLocal(point, &local[0]);
1374 inside = fRight->Contains(&local[0]);
1375 return inside;
1376}
1377
1378////////////////////////////////////////////////////////////////////////////////
1379/// Compute minimum distance to shape vertices
1380
1382{
1383 return 9999;
1384}
1385
1386////////////////////////////////////////////////////////////////////////////////
1387/// Compute distance from a given point inside to the shape boundary.
1388
1390 Double_t step, Double_t *safe) const
1391{
1392 if (iact<3 && safe) {
1393 // compute safe distance
1394 *safe = Safety(point,kTRUE);
1395 if (iact==0) return TGeoShape::Big();
1396 if (iact==1 && step<*safe) return TGeoShape::Big();
1397 }
1398 TGeoBoolNode *node = (TGeoBoolNode*)this;
1399 Double_t local[3], ldir[3], rdir[3];
1400 Double_t d1, d2, snxt=0.;
1401 fLeftMat->MasterToLocal(point, &local[0]);
1402 fLeftMat->MasterToLocalVect(dir, &ldir[0]);
1403 fRightMat->MasterToLocalVect(dir, &rdir[0]);
1404 d1 = fLeft->DistFromInside(&local[0], &ldir[0], iact, step, safe);
1405 fRightMat->MasterToLocal(point, &local[0]);
1406 d2 = fRight->DistFromInside(&local[0], &rdir[0], iact, step, safe);
1407 if (d1<d2) {
1408 snxt = d1;
1409 node->SetSelected(1);
1410 } else {
1411 snxt = d2;
1412 node->SetSelected(2);
1413 }
1414 return snxt;
1415}
1416
1417////////////////////////////////////////////////////////////////////////////////
1418/// Compute distance from a given point outside to the shape.
1419
1421 Double_t step, Double_t *safe) const
1422{
1424 if (iact<3 && safe) {
1425 // compute safe distance
1426 *safe = Safety(point,kFALSE);
1427 if (iact==0) return TGeoShape::Big();
1428 if (iact==1 && step<*safe) return TGeoShape::Big();
1429 }
1430 TGeoBoolNode *node = (TGeoBoolNode*)this;
1431 Double_t lpt[3], rpt[3], master[3], ldir[3], rdir[3];
1432 memcpy(master, point, 3*sizeof(Double_t));
1433 Int_t i;
1434 Double_t d1 = 0.;
1435 Double_t d2 = 0.;
1436 fLeftMat->MasterToLocal(point, lpt);
1437 fRightMat->MasterToLocal(point, rpt);
1438 fLeftMat->MasterToLocalVect(dir, ldir);
1439 fRightMat->MasterToLocalVect(dir, rdir);
1440 Bool_t inleft = fLeft->Contains(lpt);
1441 Bool_t inright = fRight->Contains(rpt);
1442 node->SetSelected(0);
1443 Double_t snext = 0.0;
1444 if (inleft && inright) {
1445 // It is vey likely to have a numerical issue and the point should
1446 // be logically outside one of the shapes
1447 d1 = fLeft->DistFromInside(lpt,ldir,3);
1448 d2 = fRight->DistFromInside(rpt,rdir,3);
1449 if (d1<1.E-3) inleft = kFALSE;
1450 if (d2<1.E-3) inright = kFALSE;
1451 if (inleft && inright) return snext;
1452 }
1453
1454 while (1) {
1455 d1 = d2 = 0;
1456 if (!inleft) {
1457 d1 = fLeft->DistFromOutside(lpt,ldir,3);
1458 d1 = TMath::Max(d1,tol);
1459 if (d1>1E20) return TGeoShape::Big();
1460 }
1461 if (!inright) {
1462 d2 = fRight->DistFromOutside(rpt,rdir,3);
1463 d2 = TMath::Max(d2,tol);
1464 if (d2>1E20) return TGeoShape::Big();
1465 }
1466
1467 if (d1>d2) {
1468 // propagate to left shape
1469 snext += d1;
1470 node->SetSelected(1);
1471 inleft = kTRUE;
1472 for (i=0; i<3; i++) master[i] += d1*dir[i];
1473 fRightMat->MasterToLocal(master,rpt);
1474 // Push rpt to avoid a bad boundary condition
1475 for (i=0; i<3; i++) rpt[i] += tol*rdir[i];
1476 // check if propagated point is inside right shape
1477 inright = fRight->Contains(rpt);
1478 if (inright) return snext;
1479 // here inleft=true, inright=false
1480 } else {
1481 // propagate to right shape
1482 snext += d2;
1483 node->SetSelected(2);
1484 inright = kTRUE;
1485 for (i=0; i<3; i++) master[i] += d2*dir[i];
1486 fLeftMat->MasterToLocal(master,lpt);
1487 // Push lpt to avoid a bad boundary condition
1488 for (i=0; i<3; i++) lpt[i] += tol*ldir[i];
1489 // check if propagated point is inside left shape
1490 inleft = fLeft->Contains(lpt);
1491 if (inleft) return snext;
1492 // here inleft=false, inright=true
1493 }
1494 }
1495 return snext;
1496}
1497
1498////////////////////////////////////////////////////////////////////////////////
1499/// Returns number of vertices for the composite shape described by this intersection.
1500
1502{
1503 Int_t itot=0;
1504 Double_t point[3];
1505 Double_t tolerance = TGeoShape::Tolerance();
1506 if (fNpoints) return fNpoints;
1507 Int_t nleft = fLeft->GetNmeshVertices();
1508 Int_t nright = fRight->GetNmeshVertices();
1509 Double_t *points = new Double_t[3*(nleft+nright)];
1510 Double_t *points1 = new Double_t[3*nleft];
1511 fLeft->SetPoints(points1);
1512 for (Int_t i=0; i<nleft; i++) {
1513 if (TMath::Abs(points1[3*i])<tolerance && TMath::Abs(points1[3*i+1])<tolerance) continue;
1514 fLeftMat->LocalToMaster(&points1[3*i], &points[3*itot]);
1515 fRightMat->MasterToLocal(&points[3*itot], point);
1516 if (fRight->Contains(point)) itot++;
1517 }
1518 Double_t *points2 = new Double_t[3*nright];
1519 fRight->SetPoints(points2);
1520 for (Int_t i=0; i<nright; i++) {
1521 if (TMath::Abs(points2[3*i])<tolerance && TMath::Abs(points2[3*i+1])<tolerance) continue;
1522 fRightMat->LocalToMaster(&points2[3*i], &points[3*itot]);
1523 fLeftMat->MasterToLocal(&points[3*itot], point);
1524 if (fLeft->Contains(point)) itot++;
1525 }
1526
1527 AssignPoints(itot, points);
1528
1529 delete [] points1;
1530 delete [] points2;
1531 delete [] points;
1532 return fNpoints;
1533}
1534
1535////////////////////////////////////////////////////////////////////////////////
1536/// Compute safety distance for a union node;
1537
1539{
1540 Double_t local1[3], local2[3];
1541 fLeftMat->MasterToLocal(point,local1);
1542 Bool_t in1 = fLeft->Contains(local1);
1543 fRightMat->MasterToLocal(point,local2);
1544 Bool_t in2 = fRight->Contains(local2);
1545 Bool_t intrue = in1 & in2;
1546 if (in^intrue) return 0.0;
1547 Double_t saf1 = fLeft->Safety(local1, in1);
1548 Double_t saf2 = fRight->Safety(local2, in2);
1549 if (in1 && in2) return TMath::Min(saf1, saf2);
1550 if (in1) return saf2;
1551 if (in2) return saf1;
1552 return TMath::Max(saf1,saf2);
1553}
1554
1555////////////////////////////////////////////////////////////////////////////////
1556/// Save a primitive as a C++ statement(s) on output stream "out".
1557
1558void TGeoIntersection::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
1559{
1560 TGeoBoolNode::SavePrimitive(out,option);
1561 out << " pBoolNode = new TGeoIntersection(";
1562 out << fLeft->GetPointerName() << ",";
1563 out << fRight->GetPointerName() << ",";
1564 if (!fLeftMat->IsIdentity()) out << fLeftMat->GetPointerName() << ",";
1565 else out << "0,";
1566 if (!fRightMat->IsIdentity()) out << fRightMat->GetPointerName() << ");" << std::endl;
1567 else out << "0);" << std::endl;
1568}
1569
1570////////////////////////////////////////////////////////////////////////////////
1571/// Register 3D size of this shape.
1572
1574{
1576}
int Int_t
Definition: RtypesCore.h:43
const Bool_t kFALSE
Definition: RtypesCore.h:90
bool Bool_t
Definition: RtypesCore.h:61
double Double_t
Definition: RtypesCore.h:57
float Float_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:89
const char Option_t
Definition: RtypesCore.h:64
#define ClassImp(name)
Definition: Rtypes.h:361
R__EXTERN TGeoManager * gGeoManager
Definition: TGeoManager.h:600
R__EXTERN TGeoIdentity * gGeoIdentity
Definition: TGeoMatrix.h:478
float xmin
Definition: THbookFile.cxx:93
float ymin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
float ymax
Definition: THbookFile.cxx:93
#define gPad
Definition: TVirtualPad.h:287
point * points
Definition: X3DBuffer.c:22
Generic 3D primitive description class.
Definition: TBuffer3D.h:18
@ kCSUnion
Definition: TBuffer3D.h:43
@ kCSDifference
Definition: TBuffer3D.h:43
@ kCSIntersection
Definition: TBuffer3D.h:43
Box class.
Definition: TGeoBBox.h:18
Base class for Boolean operations between two shapes.
Definition: TGeoBoolNode.h:25
virtual void Sizeof3D() const
Register size of this 3D object.
Int_t fNpoints
Definition: TGeoBoolNode.h:51
Bool_t MakeBranch(const char *expr, Bool_t left)
Mutex for thread data access.
TGeoMatrix * fLeftMat
Definition: TGeoBoolNode.h:49
void ClearThreadData() const
TGeoShape * fLeft
Definition: TGeoBoolNode.h:47
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.
Definition: TGeoBoolNode.h:55
std::mutex fMutex
Size for the navigation data array.
Definition: TGeoBoolNode.h:56
std::vector< ThreadData_t * > fThreadData
array of mesh points
Definition: TGeoBoolNode.h:54
void RegisterMatrices()
Register all matrices of the boolean node and descendents.
TGeoShape * fRight
Definition: TGeoBoolNode.h:48
virtual ~TGeoBoolNode()
Destructor.
virtual Int_t GetNpoints()=0
Double_t * fPoints
number of points on the mesh
Definition: TGeoBoolNode.h:52
TGeoBoolNode()
Default constructor.
TGeoMatrix * fRightMat
Definition: TGeoBoolNode.h:50
ThreadData_t & GetThreadData() const
void SetSelected(Int_t sel)
Set the selected branch.
Class handling Boolean composition of shapes.
Matrix class used for computing global transformations Should NOT be used for node definition.
Definition: TGeoMatrix.h:421
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.
virtual ~TGeoIntersection()
Destructor — deletion of components handled by TGeoManager class.
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.
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
Definition: TGeoManager.h:488
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
Definition: TGeoManager.h:493
static Int_t ThreadId()
Translates the current thread id to an ordinal number.
Geometrical transformation package.
Definition: TGeoMatrix.h:41
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
Definition: TGeoMatrix.cxx:363
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
Definition: TGeoMatrix.cxx:406
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
Definition: TGeoMatrix.cxx:431
virtual void RegisterYourself()
Register the matrix in the current manager, which will become the owner.
Definition: TGeoMatrix.cxx:526
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
Definition: TGeoMatrix.cxx:339
Bool_t IsIdentity() const
Definition: TGeoMatrix.h:66
char * GetPointerName() const
Provide a pointer name containing uid.
Definition: TGeoMatrix.cxx:294
Base abstract class for all shapes.
Definition: TGeoShape.h:26
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.
Definition: TGeoShape.cxx:689
static Double_t Big()
Definition: TGeoShape.h:88
virtual void CreateThreadData(Int_t)
Definition: TGeoShape.h:68
Bool_t IsValid() const
Definition: TGeoShape.h:140
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.
Definition: TGeoShape.cxx:699
virtual const char * GetName() const
Get the shape name.
Definition: TGeoShape.cxx:248
virtual Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=0) const =0
virtual void ComputeBBox()=0
virtual Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=0) const =0
static TGeoMatrix * GetTransform()
Returns current transformation matrix that applies to shape.
Definition: TGeoShape.cxx:536
virtual Bool_t Contains(const Double_t *point) const =0
@ kGeoHalfSpace
Definition: TGeoShape.h:63
static Double_t Tolerance()
Definition: TGeoShape.h:91
virtual void SetPoints(Double_t *points) const =0
Bool_t TestShapeBit(UInt_t f) const
Definition: TGeoShape.h:163
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.
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.
virtual ~TGeoSubtraction()
Destructor — deletion of components handled by TGeoManager class.
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.
virtual ~TGeoUnion()
Destructor — deletion of components handled by TGeoManager class.
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()
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.
virtual TObject * FindObject(const char *name) const
Find an object in this collection using its name.
Definition: TObjArray.cxx:415
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:877
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save a primitive as a C++ statement(s) on output stream "out".
Definition: TObject.cxx:664
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:891
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:919
Basic string class.
Definition: TString.h:131
Ssiz_t Length() const
Definition: TString.h:405
const char * Data() const
Definition: TString.h:364
Abstract 3D shapes viewer.
virtual Int_t AddObject(const TBuffer3D &buffer, Bool_t *addChildren=0)=0
virtual void AddCompositeOp(UInt_t operation)=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)
Definition: TMathBase.h:212
constexpr Double_t E()
Base of natural log:
Definition: TMath.h:97
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
void Sort(Index n, const Element *a, Index *index, Bool_t down=kTRUE)
Definition: TMathBase.h:362
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
#define snext(osub1, osub2)
Definition: triangle.c:1167