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 <iostream>
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)) {
316 left->PaintComposite(option);
317 } else if (fLeft) {
318 const TBuffer3D & leftBuffer = fLeft->GetBuffer3D(TBuffer3D::kAll, localFrame);
319 viewer->AddObject(leftBuffer);
320 }
321
322 // Setup matrix and fetch/add the right component buffer
323 *glmat = &mat;
324 glmat->Multiply(fRightMat);
325 //fRight->Paint(option);
326 if (TGeoCompositeShape *right = dynamic_cast<TGeoCompositeShape *>(fRight)) {
327 right->PaintComposite(option);
328 } else if (fRight) {
329 const TBuffer3D & rightBuffer = fRight->GetBuffer3D(TBuffer3D::kAll, localFrame);
330 viewer->AddObject(rightBuffer);
331 }
332
333 *glmat = &mat;
334}
335
336////////////////////////////////////////////////////////////////////////////////
337/// Register all matrices of the boolean node and descendents.
338
340{
343 if (fLeft->IsComposite()) ((TGeoCompositeShape*)fLeft)->GetBoolNode()->RegisterMatrices();
344 if (fRight->IsComposite()) ((TGeoCompositeShape*)fRight)->GetBoolNode()->RegisterMatrices();
345}
346
347////////////////////////////////////////////////////////////////////////////////
348/// Replace one of the matrices. Does not work with TGeoIdentity. Returns true
349/// if replacement was successful.
350
352{
353 if (mat==gGeoIdentity || newmat==gGeoIdentity) {
354 Error("ReplaceMatrix", "Matrices should not be gGeoIdentity. Use default matrix constructor to represent identities.");
355 return kFALSE;
356 }
357 if (!mat || !newmat) {
358 Error("ReplaceMatrix", "Matrices should not be null pointers.");
359 return kFALSE;
360 }
361 Bool_t replaced = kFALSE;
362 if (fLeftMat == mat) {
363 fLeftMat = newmat;
364 replaced = kTRUE;
365 }
366 if (fRightMat == mat) {
367 fRightMat = newmat;
368 replaced = kTRUE;
369 }
370 return replaced;
371}
372
373////////////////////////////////////////////////////////////////////////////////
374/// Save a primitive as a C++ statement(s) on output stream "out".
375
376void TGeoBoolNode::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
377{
378 fLeft->SavePrimitive(out,option);
379 fRight->SavePrimitive(out,option);
380 if (!fLeftMat->IsIdentity()) {
382 fLeftMat->SavePrimitive(out,option);
383 }
384 if (!fRightMat->IsIdentity()) {
386 fRightMat->SavePrimitive(out,option);
387 }
388}
389
390////////////////////////////////////////////////////////////////////////////////
391/// Fill buffer with shape vertices.
392
394{
395 TGeoBoolNode *bn = (TGeoBoolNode*)this;
396 Int_t npoints = bn->GetNpoints();
397 memcpy(points, fPoints, 3*npoints*sizeof(Double_t));
398}
399
400////////////////////////////////////////////////////////////////////////////////
401/// Fill buffer with shape vertices.
402
404{
405 TGeoBoolNode *bn = (TGeoBoolNode*)this;
406 Int_t npoints = bn->GetNpoints();
407 for (Int_t i=0; i<3*npoints; i++) points[i] = fPoints[i];
408}
409
410////////////////////////////////////////////////////////////////////////////////
411/// Register size of this 3D object
412
414{
415 fLeft->Sizeof3D();
416 fRight->Sizeof3D();
417}
418
420
421////////////////////////////////////////////////////////////////////////////////
422/// Make a clone of this. Pointers are preserved.
423
425{
426 return new TGeoUnion(fLeft, fRight, fLeftMat, fRightMat);
427}
428
429////////////////////////////////////////////////////////////////////////////////
430/// Paint method.
431
433{
434 TVirtualViewer3D *viewer = gPad->GetViewer3D();
435
436 if (!viewer) {
437 Error("Paint", "gPad->GetViewer3D() returned 0, cannot work with composite!\n");
438 return;
439 }
440
442
443 TGeoBoolNode::Paint(option);
444}
445
446////////////////////////////////////////////////////////////////////////////////
447/// Default constructor
448
450{
451}
452
453////////////////////////////////////////////////////////////////////////////////
454/// Constructor
455
456TGeoUnion::TGeoUnion(const char *expr1, const char *expr2)
457 :TGeoBoolNode(expr1, expr2)
458{
459}
460
461////////////////////////////////////////////////////////////////////////////////
462/// Constructor providing pointers to components
463
465 :TGeoBoolNode(left,right,lmat,rmat)
466{
468 Fatal("TGeoUnion", "Unions with a half-space (%s + %s) not allowed", left->GetName(), right->GetName());
469 }
470}
471
472////////////////////////////////////////////////////////////////////////////////
473/// Destructor
474/// --- deletion of components handled by TGeoManager class.
475
477{
478}
479
480////////////////////////////////////////////////////////////////////////////////
481/// Compute bounding box corresponding to a union of two shapes.
482
484{
485 if (((TGeoBBox*)fLeft)->IsNullBox()) fLeft->ComputeBBox();
486 if (((TGeoBBox*)fRight)->IsNullBox()) fRight->ComputeBBox();
487 Double_t vert[48];
488 Double_t pt[3];
489 Int_t i;
490 Double_t xmin, xmax, ymin, ymax, zmin, zmax;
491 xmin = ymin = zmin = TGeoShape::Big();
492 xmax = ymax = zmax = -TGeoShape::Big();
493 ((TGeoBBox*)fLeft)->SetBoxPoints(&vert[0]);
494 ((TGeoBBox*)fRight)->SetBoxPoints(&vert[24]);
495 for (i=0; i<8; i++) {
496 fLeftMat->LocalToMaster(&vert[3*i], &pt[0]);
497 if (pt[0]<xmin) xmin=pt[0];
498 if (pt[0]>xmax) xmax=pt[0];
499 if (pt[1]<ymin) ymin=pt[1];
500 if (pt[1]>ymax) ymax=pt[1];
501 if (pt[2]<zmin) zmin=pt[2];
502 if (pt[2]>zmax) zmax=pt[2];
503 }
504 for (i=8; i<16; i++) {
505 fRightMat->LocalToMaster(&vert[3*i], &pt[0]);
506 if (pt[0]<xmin) xmin=pt[0];
507 if (pt[0]>xmax) xmax=pt[0];
508 if (pt[1]<ymin) ymin=pt[1];
509 if (pt[1]>ymax) ymax=pt[1];
510 if (pt[2]<zmin) zmin=pt[2];
511 if (pt[2]>zmax) zmax=pt[2];
512 }
513 dx = 0.5*(xmax-xmin);
514 origin[0] = 0.5*(xmin+xmax);
515 dy = 0.5*(ymax-ymin);
516 origin[1] = 0.5*(ymin+ymax);
517 dz = 0.5*(zmax-zmin);
518 origin[2] = 0.5*(zmin+zmax);
519}
520
521////////////////////////////////////////////////////////////////////////////////
522/// Find if a union of two shapes contains a given point
523
525{
526 Double_t local[3];
527 fLeftMat->MasterToLocal(point, &local[0]);
528 Bool_t inside = fLeft->Contains(&local[0]);
529 if (inside) return kTRUE;
530 fRightMat->MasterToLocal(point, &local[0]);
531 inside = fRight->Contains(&local[0]);
532 return inside;
533}
534
535////////////////////////////////////////////////////////////////////////////////
536/// Normal computation in POINT. The orientation is chosen so that DIR.dot.NORM>0.
537
538void TGeoUnion::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
539{
541 norm[0] = norm[1] = 0.;
542 norm[2] = 1.;
543 Double_t local[3];
544 Double_t ldir[3], lnorm[3];
545 if (td.fSelected == 1) {
546 fLeftMat->MasterToLocal(point, local);
547 fLeftMat->MasterToLocalVect(dir, ldir);
548 fLeft->ComputeNormal(local,ldir,lnorm);
549 fLeftMat->LocalToMasterVect(lnorm, norm);
550 return;
551 }
552 if (td.fSelected == 2) {
553 fRightMat->MasterToLocal(point, local);
554 fRightMat->MasterToLocalVect(dir, ldir);
555 fRight->ComputeNormal(local,ldir,lnorm);
556 fRightMat->LocalToMasterVect(lnorm, norm);
557 return;
558 }
559 fLeftMat->MasterToLocal(point, local);
560 if (fLeft->Contains(local)) {
561 fLeftMat->MasterToLocalVect(dir, ldir);
562 fLeft->ComputeNormal(local,ldir,lnorm);
563 fLeftMat->LocalToMasterVect(lnorm, norm);
564 return;
565 }
566 fRightMat->MasterToLocal(point, local);
567 if (fRight->Contains(local)) {
568 fRightMat->MasterToLocalVect(dir, ldir);
569 fRight->ComputeNormal(local,ldir,lnorm);
570 fRightMat->LocalToMasterVect(lnorm, norm);
571 return;
572 }
573 // Propagate forward/backward to see which of the components is intersected first
574 local[0] = point[0] + 1E-5*dir[0];
575 local[1] = point[1] + 1E-5*dir[1];
576 local[2] = point[2] + 1E-5*dir[2];
577
578 if (!Contains(local)) {
579 local[0] = point[0] - 1E-5*dir[0];
580 local[1] = point[1] - 1E-5*dir[1];
581 local[2] = point[2] - 1E-5*dir[2];
582 if (!Contains(local)) return;
583 }
584 ComputeNormal(local,dir,norm);
585}
586
587////////////////////////////////////////////////////////////////////////////////
588/// Compute minimum distance to shape vertices.
589
591{
592 return 9999;
593}
594
595////////////////////////////////////////////////////////////////////////////////
596/// Computes distance from a given point inside the shape to its boundary.
597
599 Double_t step, Double_t *safe) const
600{
601 if (iact<3 && safe) {
602 // compute safe distance
603 *safe = Safety(point,kTRUE);
604 if (iact==0) return TGeoShape::Big();
605 if (iact==1 && step<*safe) return TGeoShape::Big();
606 }
607
608 Double_t local[3], local1[3], master[3], ldir[3], rdir[3], pushed[3];
609 memcpy(master, point, 3*sizeof(Double_t));
610 Int_t i;
611 TGeoBoolNode *node = (TGeoBoolNode*)this;
612 Double_t d1=0., d2=0., snxt=0., eps=0.;
613 fLeftMat->MasterToLocalVect(dir, ldir);
614 fRightMat->MasterToLocalVect(dir, rdir);
615 fLeftMat->MasterToLocal(point, local);
616 Bool_t inside1 = fLeft->Contains(local);
617 if (inside1) d1 = fLeft->DistFromInside(local, ldir, 3);
618 else memcpy(local1, local, 3*sizeof(Double_t));
619 fRightMat->MasterToLocal(point, local);
620 Bool_t inside2 = fRight->Contains(local);
621 if (inside2) d2 = fRight->DistFromInside(local, rdir, 3);
622 if (!(inside1 | inside2)) {
623 // This is a pathological case when the point is on the boundary
624 d1 = fLeft->DistFromOutside(local1, ldir, 3);
625 if (d1<1.E-3) {
626 eps = d1+TGeoShape::Tolerance();
627 for (i=0; i<3; i++) local1[i] += eps*ldir[i];
628 inside1 = kTRUE;
629 d1 = fLeft->DistFromInside(local1, ldir, 3);
630 d1 += eps;
631 } else {
632 d2 = fRight->DistFromOutside(local, rdir, 3);
633 if (d2<1.E-3) {
634 eps = d2+TGeoShape::Tolerance();
635 for (i=0; i<3; i++) local[i] += eps*rdir[i];
636 inside2 = kTRUE;
637 d2 = fRight->DistFromInside(local, rdir, 3);
638 d2 += eps;
639 }
640 }
641 }
642 while (inside1 || inside2) {
643 if (inside1 && inside2) {
644 if (d1<d2) {
645 snxt += d1;
646 node->SetSelected(1);
647 // propagate to exit of left shape
648 inside1 = kFALSE;
649 for (i=0; i<3; i++) master[i] += d1*dir[i];
650 // check if propagated point is in right shape
651 fRightMat->MasterToLocal(master, local);
652 inside2 = fRight->Contains(local);
653 if (!inside2) return snxt;
654 d2 = fRight->DistFromInside(local, rdir, 3);
655 if (d2 < TGeoShape::Tolerance()) return snxt;
656 } else {
657 snxt += d2;
658 node->SetSelected(2);
659 // propagate to exit of right shape
660 inside2 = kFALSE;
661 for (i=0; i<3; i++) master[i] += d2*dir[i];
662 // check if propagated point is in left shape
663 fLeftMat->MasterToLocal(master, local);
664 inside1 = fLeft->Contains(local);
665 if (!inside1) return snxt;
666 d1 = fLeft->DistFromInside(local, ldir, 3);
667 if (d1 < TGeoShape::Tolerance()) return snxt;
668 }
669 }
670 if (inside1) {
671 snxt += d1;
672 node->SetSelected(1);
673 // propagate to exit of left shape
674 inside1 = kFALSE;
675 for (i=0; i<3; i++) {
676 master[i] += d1*dir[i];
677 pushed[i] = master[i]+(1.+d1)*TGeoShape::Tolerance()*dir[i];
678 }
679 // check if propagated point is in right shape
680 fRightMat->MasterToLocal(pushed, local);
681 inside2 = fRight->Contains(local);
682 if (!inside2) return snxt;
683 d2 = fRight->DistFromInside(local, rdir, 3);
684 if (d2 < TGeoShape::Tolerance()) return snxt;
685 d2 += (1.+d1)*TGeoShape::Tolerance();
686 }
687 if (inside2) {
688 snxt += d2;
689 node->SetSelected(2);
690 // propagate to exit of right shape
691 inside2 = kFALSE;
692 for (i=0; i<3; i++) {
693 master[i] += d2*dir[i];
694 pushed[i] = master[i]+(1.+d2)*TGeoShape::Tolerance()*dir[i];
695 }
696 // check if propagated point is in left shape
697 fLeftMat->MasterToLocal(pushed, local);
698 inside1 = fLeft->Contains(local);
699 if (!inside1) return snxt;
700 d1 = fLeft->DistFromInside(local, ldir, 3);
701 if (d1 < TGeoShape::Tolerance()) return snxt;
702 d1 += (1.+d2)*TGeoShape::Tolerance();
703 }
704 }
705 return snxt;
706}
707
708////////////////////////////////////////////////////////////////////////////////
709/// Compute distance from a given outside point to the shape.
710
712 Double_t step, Double_t *safe) const
713{
714 if (iact<3 && safe) {
715 // compute safe distance
716 *safe = Safety(point,kFALSE);
717 if (iact==0) return TGeoShape::Big();
718 if (iact==1 && step<*safe) return TGeoShape::Big();
719 }
720 TGeoBoolNode *node = (TGeoBoolNode*)this;
721 Double_t local[3], ldir[3], rdir[3];
722 Double_t d1, d2, snxt;
723 fLeftMat->MasterToLocal(point, &local[0]);
724 fLeftMat->MasterToLocalVect(dir, &ldir[0]);
725 fRightMat->MasterToLocalVect(dir, &rdir[0]);
726 d1 = fLeft->DistFromOutside(&local[0], &ldir[0], iact, step, safe);
727 fRightMat->MasterToLocal(point, &local[0]);
728 d2 = fRight->DistFromOutside(&local[0], &rdir[0], iact, step, safe);
729 if (d1<d2) {
730 snxt = d1;
731 node->SetSelected(1);
732 } else {
733 snxt = d2;
734 node->SetSelected(2);
735 }
736 return snxt;
737}
738
739////////////////////////////////////////////////////////////////////////////////
740/// Returns number of vertices for the composite shape described by this union.
741
743{
744 Int_t itot=0;
745 Double_t point[3];
746 Double_t tolerance = TGeoShape::Tolerance();
747 if (fNpoints) return fNpoints;
748 // Local points for the left shape
749 Int_t nleft = fLeft->GetNmeshVertices();
750 Double_t *points1 = new Double_t[3*nleft];
751 fLeft->SetPoints(points1);
752 // Local points for the right shape
753 Int_t nright = fRight->GetNmeshVertices();
754 Double_t *points2 = new Double_t[3*nright];
755 fRight->SetPoints(points2);
756 Double_t *points = new Double_t[3*(nleft+nright)];
757 for (Int_t i=0; i<nleft; i++) {
758 if (TMath::Abs(points1[3*i])<tolerance && TMath::Abs(points1[3*i+1])<tolerance) continue;
759 fLeftMat->LocalToMaster(&points1[3*i], &points[3*itot]);
760 fRightMat->MasterToLocal(&points[3*itot], point);
761 if (!fRight->Contains(point)) itot++;
762 }
763 for (Int_t i=0; i<nright; i++) {
764 if (TMath::Abs(points2[3*i])<tolerance && TMath::Abs(points2[3*i+1])<tolerance) continue;
765 fRightMat->LocalToMaster(&points2[3*i], &points[3*itot]);
766 fLeftMat->MasterToLocal(&points[3*itot], point);
767 if (!fLeft->Contains(point)) itot++;
768 }
769
770 AssignPoints(itot, points);
771
772 delete [] points1;
773 delete [] points2;
774 delete [] points;
775 return fNpoints;
776}
777
778////////////////////////////////////////////////////////////////////////////////
779/// Compute safety distance for a union node;
780
782{
783 Double_t local1[3], local2[3];
784 fLeftMat->MasterToLocal(point,local1);
785 Bool_t in1 = fLeft->Contains(local1);
786 fRightMat->MasterToLocal(point,local2);
787 Bool_t in2 = fRight->Contains(local2);
788 Bool_t intrue = in1 | in2;
789 if (intrue^in) return 0.0;
790 Double_t saf1 = fLeft->Safety(local1, in1);
791 Double_t saf2 = fRight->Safety(local2, in2);
792 if (in1 && in2) return TMath::Min(saf1, saf2);
793 if (in1) return saf1;
794 if (in2) return saf2;
795 return TMath::Min(saf1,saf2);
796}
797
798////////////////////////////////////////////////////////////////////////////////
799/// Save a primitive as a C++ statement(s) on output stream "out".
800
801void TGeoUnion::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
802{
803 TGeoBoolNode::SavePrimitive(out,option);
804 out << " pBoolNode = new TGeoUnion(";
805 out << fLeft->GetPointerName() << ",";
806 out << fRight->GetPointerName() << ",";
807 if (!fLeftMat->IsIdentity()) out << fLeftMat->GetPointerName() << ",";
808 else out << "0,";
809 if (!fRightMat->IsIdentity()) out << fRightMat->GetPointerName() << ");" << std::endl;
810 else out << "0);" << std::endl;
811}
812
813////////////////////////////////////////////////////////////////////////////////
814/// Register 3D size of this shape.
815
817{
819}
820
821
823
824////////////////////////////////////////////////////////////////////////////////
825/// Make a clone of this. Pointers are preserved.
826
828{
830}
831
832////////////////////////////////////////////////////////////////////////////////
833/// Paint method.
834
836{
837 TVirtualViewer3D *viewer = gPad->GetViewer3D();
838
839 if (!viewer) {
840 Error("Paint", "gPad->GetViewer3D() returned 0, cannot work with composite!\n");
841 return;
842 }
843
845
846 TGeoBoolNode::Paint(option);
847}
848
849////////////////////////////////////////////////////////////////////////////////
850/// Default constructor
851
853{
854}
855
856////////////////////////////////////////////////////////////////////////////////
857/// Constructor
858
859TGeoSubtraction::TGeoSubtraction(const char *expr1, const char *expr2)
860 :TGeoBoolNode(expr1, expr2)
861{
862}
863
864////////////////////////////////////////////////////////////////////////////////
865/// Constructor providing pointers to components
866
868 :TGeoBoolNode(left,right,lmat,rmat)
869{
871 Fatal("TGeoSubstraction", "Subtractions from a half-space (%s) not allowed", left->GetName());
872 }
873}
874
875////////////////////////////////////////////////////////////////////////////////
876/// Destructor
877/// --- deletion of components handled by TGeoManager class.
878
880{
881}
882
883////////////////////////////////////////////////////////////////////////////////
884/// Compute bounding box corresponding to a subtraction of two shapes.
885
887{
889 if (box->IsNullBox()) fLeft->ComputeBBox();
890 Double_t vert[24];
891 Double_t pt[3];
892 Int_t i;
893 Double_t xmin, xmax, ymin, ymax, zmin, zmax;
894 xmin = ymin = zmin = TGeoShape::Big();
895 xmax = ymax = zmax = -TGeoShape::Big();
896 box->SetBoxPoints(&vert[0]);
897 for (i=0; i<8; i++) {
898 fLeftMat->LocalToMaster(&vert[3*i], &pt[0]);
899 if (pt[0]<xmin) xmin=pt[0];
900 if (pt[0]>xmax) xmax=pt[0];
901 if (pt[1]<ymin) ymin=pt[1];
902 if (pt[1]>ymax) ymax=pt[1];
903 if (pt[2]<zmin) zmin=pt[2];
904 if (pt[2]>zmax) zmax=pt[2];
905 }
906 dx = 0.5*(xmax-xmin);
907 origin[0] = 0.5*(xmin+xmax);
908 dy = 0.5*(ymax-ymin);
909 origin[1] = 0.5*(ymin+ymax);
910 dz = 0.5*(zmax-zmin);
911 origin[2] = 0.5*(zmin+zmax);
912}
913
914////////////////////////////////////////////////////////////////////////////////
915/// Normal computation in POINT. The orientation is chosen so that DIR.dot.NORM>0.
916
917void TGeoSubtraction::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
918{
920 norm[0] = norm[1] = 0.;
921 norm[2] = 1.;
922 Double_t local[3], ldir[3], lnorm[3];
923 if (td.fSelected == 1) {
924 fLeftMat->MasterToLocal(point, local);
925 fLeftMat->MasterToLocalVect(dir, ldir);
926 fLeft->ComputeNormal(local,ldir,lnorm);
927 fLeftMat->LocalToMasterVect(lnorm, norm);
928 return;
929 }
930 if (td.fSelected == 2) {
931 fRightMat->MasterToLocal(point, local);
932 fRightMat->MasterToLocalVect(dir, ldir);
933 fRight->ComputeNormal(local,ldir,lnorm);
934 fRightMat->LocalToMasterVect(lnorm, norm);
935 return;
936 }
937 fRightMat->MasterToLocal(point,local);
938 if (fRight->Contains(local)) {
939 fRightMat->MasterToLocalVect(dir,ldir);
940 fRight->ComputeNormal(local,ldir, lnorm);
941 fRightMat->LocalToMasterVect(lnorm,norm);
942 return;
943 }
944 fLeftMat->MasterToLocal(point,local);
945 if (!fLeft->Contains(local)) {
946 fLeftMat->MasterToLocalVect(dir,ldir);
947 fLeft->ComputeNormal(local,ldir, lnorm);
948 fLeftMat->LocalToMasterVect(lnorm,norm);
949 return;
950 }
951 // point is inside left shape, but not inside the right
952 local[0] = point[0]+1E-5*dir[0];
953 local[1] = point[1]+1E-5*dir[1];
954 local[2] = point[2]+1E-5*dir[2];
955 if (Contains(local)) {
956 local[0] = point[0]-1E-5*dir[0];
957 local[1] = point[1]-1E-5*dir[1];
958 local[2] = point[2]-1E-5*dir[2];
959 if (Contains(local)) return;
960 }
961 ComputeNormal(local,dir,norm);
962}
963
964////////////////////////////////////////////////////////////////////////////////
965/// Find if a subtraction of two shapes contains a given point
966
968{
969 Double_t local[3];
970 fLeftMat->MasterToLocal(point, &local[0]);
971 Bool_t inside = fLeft->Contains(&local[0]);
972 if (!inside) return kFALSE;
973 fRightMat->MasterToLocal(point, &local[0]);
974 inside = !fRight->Contains(&local[0]);
975 return inside;
976}
977
978////////////////////////////////////////////////////////////////////////////////
979/// Compute minimum distance to shape vertices
980
982{
983 return 9999;
984}
985
986////////////////////////////////////////////////////////////////////////////////
987/// Compute distance from a given point inside to the shape boundary.
988
990 Double_t step, Double_t *safe) const
991{
992 if (iact<3 && safe) {
993 // compute safe distance
994 *safe = Safety(point,kTRUE);
995 if (iact==0) return TGeoShape::Big();
996 if (iact==1 && step<*safe) return TGeoShape::Big();
997 }
998 TGeoBoolNode *node = (TGeoBoolNode*)this;
999 Double_t local[3], ldir[3], rdir[3];
1000 Double_t d1, d2, snxt=0.;
1001 fLeftMat->MasterToLocal(point, &local[0]);
1002 fLeftMat->MasterToLocalVect(dir, &ldir[0]);
1003 fRightMat->MasterToLocalVect(dir, &rdir[0]);
1004 d1 = fLeft->DistFromInside(&local[0], &ldir[0], iact, step, safe);
1005 fRightMat->MasterToLocal(point, &local[0]);
1006 d2 = fRight->DistFromOutside(&local[0], &rdir[0], iact, step, safe);
1007 if (d1<d2) {
1008 snxt = d1;
1009 node->SetSelected(1);
1010 } else {
1011 snxt = d2;
1012 node->SetSelected(2);
1013 }
1014 return snxt;
1015}
1016
1017////////////////////////////////////////////////////////////////////////////////
1018/// Compute distance from a given point outside to the shape.
1019
1021 Double_t step, Double_t *safe) const
1022{
1023 if (iact<3 && safe) {
1024 // compute safe distance
1025 *safe = Safety(point,kFALSE);
1026 if (iact==0) return TGeoShape::Big();
1027 if (iact==1 && step<*safe) return TGeoShape::Big();
1028 }
1029 TGeoBoolNode *node = (TGeoBoolNode*)this;
1030 Double_t local[3], master[3], ldir[3], rdir[3];
1031 memcpy(&master[0], point, 3*sizeof(Double_t));
1032 Int_t i;
1033 Double_t d1, d2, snxt=0.;
1034 fRightMat->MasterToLocal(point, &local[0]);
1035 fLeftMat->MasterToLocalVect(dir, &ldir[0]);
1036 fRightMat->MasterToLocalVect(dir, &rdir[0]);
1037 // check if inside '-'
1038 Bool_t inside = fRight->Contains(&local[0]);
1039 Double_t epsil = 0.;
1040 while (1) {
1041 if (inside) {
1042 // propagate to outside of '-'
1043 node->SetSelected(2);
1044 d1 = fRight->DistFromInside(&local[0], &rdir[0], iact, step, safe);
1045 snxt += d1+epsil;
1046 for (i=0; i<3; i++) master[i] += (d1+1E-8)*dir[i];
1047 epsil = 1.E-8;
1048 // now master outside '-'; check if inside '+'
1049 fLeftMat->MasterToLocal(&master[0], &local[0]);
1050 if (fLeft->Contains(&local[0])) return snxt;
1051 }
1052 // master outside '-' and outside '+' ; find distances to both
1053 node->SetSelected(1);
1054 fLeftMat->MasterToLocal(&master[0], &local[0]);
1055 d2 = fLeft->DistFromOutside(&local[0], &ldir[0], iact, step, safe);
1056 if (d2>1E20) return TGeoShape::Big();
1057
1058 fRightMat->MasterToLocal(&master[0], &local[0]);
1059 d1 = fRight->DistFromOutside(&local[0], &rdir[0], iact, step, safe);
1060 if (d2<d1-TGeoShape::Tolerance()) {
1061 snxt += d2+epsil;
1062 return snxt;
1063 }
1064 // propagate to '-'
1065 snxt += d1+epsil;
1066 for (i=0; i<3; i++) master[i] += (d1+1E-8)*dir[i];
1067 epsil = 1.E-8;
1068 // now inside '-' and not inside '+'
1069 fRightMat->MasterToLocal(&master[0], &local[0]);
1070 inside = kTRUE;
1071 }
1072}
1073
1074////////////////////////////////////////////////////////////////////////////////
1075/// Returns number of vertices for the composite shape described by this subtraction.
1076
1078{
1079 Int_t itot=0;
1080 Double_t point[3];
1081 Double_t tolerance = TGeoShape::Tolerance();
1082 if (fNpoints) return fNpoints;
1083 Int_t nleft = fLeft->GetNmeshVertices();
1084 Int_t nright = fRight->GetNmeshVertices();
1085 Double_t *points = new Double_t[3*(nleft+nright)];
1086 Double_t *points1 = new Double_t[3*nleft];
1087 fLeft->SetPoints(points1);
1088 for (Int_t i=0; i<nleft; i++) {
1089 if (TMath::Abs(points1[3*i])<tolerance && TMath::Abs(points1[3*i+1])<tolerance) continue;
1090 fLeftMat->LocalToMaster(&points1[3*i], &points[3*itot]);
1091 fRightMat->MasterToLocal(&points[3*itot], point);
1092 if (!fRight->Contains(point)) itot++;
1093 }
1094 Double_t *points2 = new Double_t[3*nright];
1095 fRight->SetPoints(points2);
1096 for (Int_t i=0; i<nright; i++) {
1097 if (TMath::Abs(points2[3*i])<tolerance && TMath::Abs(points2[3*i+1])<tolerance) continue;
1098 fRightMat->LocalToMaster(&points2[3*i], &points[3*itot]);
1099 fLeftMat->MasterToLocal(&points[3*itot], point);
1100 if (fLeft->Contains(point)) itot++;
1101 }
1102
1103 AssignPoints(itot, points);
1104
1105 delete [] points1;
1106 delete [] points2;
1107 delete [] points;
1108 return fNpoints;
1109}
1110
1111////////////////////////////////////////////////////////////////////////////////
1112/// Compute safety distance for a union node;
1113
1115{
1116 Double_t local1[3], local2[3];
1117 fLeftMat->MasterToLocal(point,local1);
1118 Bool_t in1 = fLeft->Contains(local1);
1119 fRightMat->MasterToLocal(point,local2);
1120 Bool_t in2 = fRight->Contains(local2);
1121 Bool_t intrue = in1 && (!in2);
1122 if (in^intrue) return 0.0;
1123 Double_t saf1 = fLeft->Safety(local1, in1);
1124 Double_t saf2 = fRight->Safety(local2, in2);
1125 if (in1 && in2) return saf2;
1126 if (in1) return TMath::Min(saf1,saf2);
1127 if (in2) return TMath::Max(saf1,saf2);
1128 return saf1;
1129}
1130
1131////////////////////////////////////////////////////////////////////////////////
1132/// Save a primitive as a C++ statement(s) on output stream "out".
1133
1134void TGeoSubtraction::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
1135{
1136 TGeoBoolNode::SavePrimitive(out,option);
1137 out << " pBoolNode = new TGeoSubtraction(";
1138 out << fLeft->GetPointerName() << ",";
1139 out << fRight->GetPointerName() << ",";
1140 if (!fLeftMat->IsIdentity()) out << fLeftMat->GetPointerName() << ",";
1141 else out << "0,";
1142 if (!fRightMat->IsIdentity()) out << fRightMat->GetPointerName() << ");" << std::endl;
1143 else out << "0);" << std::endl;
1144}
1145
1146////////////////////////////////////////////////////////////////////////////////
1147/// Register 3D size of this shape.
1148
1150{
1152}
1153
1155
1156////////////////////////////////////////////////////////////////////////////////
1157/// Make a clone of this. Pointers are preserved.
1158
1160{
1162}
1163
1164////////////////////////////////////////////////////////////////////////////////
1165/// Paint method.
1166
1168{
1169 TVirtualViewer3D *viewer = gPad->GetViewer3D();
1170
1171 if (!viewer) {
1172 Error("Paint", "gPad->GetViewer3D() returned 0, cannot work with composite!\n");
1173 return;
1174 }
1175
1177
1178 TGeoBoolNode::Paint(option);
1179}
1180
1181////////////////////////////////////////////////////////////////////////////////
1182/// Default constructor
1183
1185{
1186}
1187
1188////////////////////////////////////////////////////////////////////////////////
1189/// Constructor
1190
1191TGeoIntersection::TGeoIntersection(const char *expr1, const char *expr2)
1192 :TGeoBoolNode(expr1, expr2)
1193{
1194}
1195
1196////////////////////////////////////////////////////////////////////////////////
1197/// Constructor providing pointers to components
1198
1200 :TGeoBoolNode(left,right,lmat,rmat)
1201{
1204 if (hs1 && hs2) Fatal("ctor", "cannot intersect two half-spaces: %s * %s", left->GetName(), right->GetName());
1205}
1206
1207////////////////////////////////////////////////////////////////////////////////
1208/// Destructor
1209/// --- deletion of components handled by TGeoManager class.
1210
1212{
1213}
1214
1215////////////////////////////////////////////////////////////////////////////////
1216/// Compute bounding box corresponding to a intersection of two shapes.
1217
1219{
1222 Double_t vert[48];
1223 Double_t pt[3];
1224 Int_t i;
1225 Double_t xmin1, xmax1, ymin1, ymax1, zmin1, zmax1;
1226 Double_t xmin2, xmax2, ymin2, ymax2, zmin2, zmax2;
1227 xmin1 = ymin1 = zmin1 = xmin2 = ymin2 = zmin2 = TGeoShape::Big();
1228 xmax1 = ymax1 = zmax1 = xmax2 = ymax2 = zmax2 = -TGeoShape::Big();
1229 if (!hs1) {
1230 if (((TGeoBBox*)fLeft)->IsNullBox()) fLeft->ComputeBBox();
1231 ((TGeoBBox*)fLeft)->SetBoxPoints(&vert[0]);
1232 for (i=0; i<8; i++) {
1233 fLeftMat->LocalToMaster(&vert[3*i], &pt[0]);
1234 if (pt[0]<xmin1) xmin1=pt[0];
1235 if (pt[0]>xmax1) xmax1=pt[0];
1236 if (pt[1]<ymin1) ymin1=pt[1];
1237 if (pt[1]>ymax1) ymax1=pt[1];
1238 if (pt[2]<zmin1) zmin1=pt[2];
1239 if (pt[2]>zmax1) zmax1=pt[2];
1240 }
1241 }
1242 if (!hs2) {
1243 if (((TGeoBBox*)fRight)->IsNullBox()) fRight->ComputeBBox();
1244 ((TGeoBBox*)fRight)->SetBoxPoints(&vert[24]);
1245 for (i=8; i<16; i++) {
1246 fRightMat->LocalToMaster(&vert[3*i], &pt[0]);
1247 if (pt[0]<xmin2) xmin2=pt[0];
1248 if (pt[0]>xmax2) xmax2=pt[0];
1249 if (pt[1]<ymin2) ymin2=pt[1];
1250 if (pt[1]>ymax2) ymax2=pt[1];
1251 if (pt[2]<zmin2) zmin2=pt[2];
1252 if (pt[2]>zmax2) zmax2=pt[2];
1253 }
1254 }
1255 if (hs1) {
1256 dx = 0.5*(xmax2-xmin2);
1257 origin[0] = 0.5*(xmax2+xmin2);
1258 dy = 0.5*(ymax2-ymin2);
1259 origin[1] = 0.5*(ymax2+ymin2);
1260 dz = 0.5*(zmax2-zmin2);
1261 origin[2] = 0.5*(zmax2+zmin2);
1262 return;
1263 }
1264 if (hs2) {
1265 dx = 0.5*(xmax1-xmin1);
1266 origin[0] = 0.5*(xmax1+xmin1);
1267 dy = 0.5*(ymax1-ymin1);
1268 origin[1] = 0.5*(ymax1+ymin1);
1269 dz = 0.5*(zmax1-zmin1);
1270 origin[2] = 0.5*(zmax1+zmin1);
1271 return;
1272 }
1273 Double_t sort[4];
1274 Int_t isort[4];
1275 sort[0] = xmin1;
1276 sort[1] = xmax1;
1277 sort[2] = xmin2;
1278 sort[3] = xmax2;
1279 TMath::Sort(4, &sort[0], &isort[0], kFALSE);
1280 if (isort[1]%2) {
1281 Warning("ComputeBBox", "shapes %s and %s do not intersect", fLeft->GetName(), fRight->GetName());
1282 dx = dy = dz = 0;
1283 memset(origin, 0, 3*sizeof(Double_t));
1284 return;
1285 }
1286 dx = 0.5*(sort[isort[2]]-sort[isort[1]]);
1287 origin[0] = 0.5*(sort[isort[1]]+sort[isort[2]]);
1288 sort[0] = ymin1;
1289 sort[1] = ymax1;
1290 sort[2] = ymin2;
1291 sort[3] = ymax2;
1292 TMath::Sort(4, &sort[0], &isort[0], kFALSE);
1293 if (isort[1]%2) {
1294 Warning("ComputeBBox", "shapes %s and %s do not intersect", fLeft->GetName(), fRight->GetName());
1295 dx = dy = dz = 0;
1296 memset(origin, 0, 3*sizeof(Double_t));
1297 return;
1298 }
1299 dy = 0.5*(sort[isort[2]]-sort[isort[1]]);
1300 origin[1] = 0.5*(sort[isort[1]]+sort[isort[2]]);
1301 sort[0] = zmin1;
1302 sort[1] = zmax1;
1303 sort[2] = zmin2;
1304 sort[3] = zmax2;
1305 TMath::Sort(4, &sort[0], &isort[0], kFALSE);
1306 if (isort[1]%2) {
1307 Warning("ComputeBBox", "shapes %s and %s do not intersect", fLeft->GetName(), fRight->GetName());
1308 dx = dy = dz = 0;
1309 memset(origin, 0, 3*sizeof(Double_t));
1310 return;
1311 }
1312 dz = 0.5*(sort[isort[2]]-sort[isort[1]]);
1313 origin[2] = 0.5*(sort[isort[1]]+sort[isort[2]]);
1314}
1315
1316////////////////////////////////////////////////////////////////////////////////
1317/// Normal computation in POINT. The orientation is chosen so that DIR.dot.NORM>0.
1318
1319void TGeoIntersection::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
1320{
1322 Double_t local[3], ldir[3], lnorm[3];
1323 norm[0] = norm[1] = 0.;
1324 norm[2] = 1.;
1325 if (td.fSelected == 1) {
1326 fLeftMat->MasterToLocal(point, local);
1327 fLeftMat->MasterToLocalVect(dir, ldir);
1328 fLeft->ComputeNormal(local,ldir,lnorm);
1329 fLeftMat->LocalToMasterVect(lnorm, norm);
1330 return;
1331 }
1332 if (td.fSelected == 2) {
1333 fRightMat->MasterToLocal(point, local);
1334 fRightMat->MasterToLocalVect(dir, ldir);
1335 fRight->ComputeNormal(local,ldir,lnorm);
1336 fRightMat->LocalToMasterVect(lnorm, norm);
1337 return;
1338 }
1339 fLeftMat->MasterToLocal(point,local);
1340 if (!fLeft->Contains(local)) {
1341 fLeftMat->MasterToLocalVect(dir,ldir);
1342 fLeft->ComputeNormal(local,ldir,lnorm);
1343 fLeftMat->LocalToMasterVect(lnorm,norm);
1344 return;
1345 }
1346 fRightMat->MasterToLocal(point,local);
1347 if (!fRight->Contains(local)) {
1348 fRightMat->MasterToLocalVect(dir,ldir);
1349 fRight->ComputeNormal(local,ldir,lnorm);
1350 fRightMat->LocalToMasterVect(lnorm,norm);
1351 return;
1352 }
1353 // point is inside intersection.
1354 local[0] = point[0] + 1E-5*dir[0];
1355 local[1] = point[1] + 1E-5*dir[1];
1356 local[2] = point[2] + 1E-5*dir[2];
1357 if (Contains(local)) {
1358 local[0] = point[0] - 1E-5*dir[0];
1359 local[1] = point[1] - 1E-5*dir[1];
1360 local[2] = point[2] - 1E-5*dir[2];
1361 if (Contains(local)) return;
1362 }
1363 ComputeNormal(local,dir,norm);
1364}
1365
1366////////////////////////////////////////////////////////////////////////////////
1367/// Find if a intersection of two shapes contains a given point
1368
1370{
1371 Double_t local[3];
1372 fLeftMat->MasterToLocal(point, &local[0]);
1373 Bool_t inside = fLeft->Contains(&local[0]);
1374 if (!inside) return kFALSE;
1375 fRightMat->MasterToLocal(point, &local[0]);
1376 inside = fRight->Contains(&local[0]);
1377 return inside;
1378}
1379
1380////////////////////////////////////////////////////////////////////////////////
1381/// Compute minimum distance to shape vertices
1382
1384{
1385 return 9999;
1386}
1387
1388////////////////////////////////////////////////////////////////////////////////
1389/// Compute distance from a given point inside to the shape boundary.
1390
1392 Double_t step, Double_t *safe) const
1393{
1394 if (iact<3 && safe) {
1395 // compute safe distance
1396 *safe = Safety(point,kTRUE);
1397 if (iact==0) return TGeoShape::Big();
1398 if (iact==1 && step<*safe) return TGeoShape::Big();
1399 }
1400 TGeoBoolNode *node = (TGeoBoolNode*)this;
1401 Double_t local[3], ldir[3], rdir[3];
1402 Double_t d1, d2, snxt=0.;
1403 fLeftMat->MasterToLocal(point, &local[0]);
1404 fLeftMat->MasterToLocalVect(dir, &ldir[0]);
1405 fRightMat->MasterToLocalVect(dir, &rdir[0]);
1406 d1 = fLeft->DistFromInside(&local[0], &ldir[0], iact, step, safe);
1407 fRightMat->MasterToLocal(point, &local[0]);
1408 d2 = fRight->DistFromInside(&local[0], &rdir[0], iact, step, safe);
1409 if (d1<d2) {
1410 snxt = d1;
1411 node->SetSelected(1);
1412 } else {
1413 snxt = d2;
1414 node->SetSelected(2);
1415 }
1416 return snxt;
1417}
1418
1419////////////////////////////////////////////////////////////////////////////////
1420/// Compute distance from a given point outside to the shape.
1421
1423 Double_t step, Double_t *safe) const
1424{
1426 if (iact<3 && safe) {
1427 // compute safe distance
1428 *safe = Safety(point,kFALSE);
1429 if (iact==0) return TGeoShape::Big();
1430 if (iact==1 && step<*safe) return TGeoShape::Big();
1431 }
1432 TGeoBoolNode *node = (TGeoBoolNode*)this;
1433 Double_t lpt[3], rpt[3], master[3], ldir[3], rdir[3];
1434 memcpy(master, point, 3*sizeof(Double_t));
1435 Int_t i;
1436 Double_t d1 = 0.;
1437 Double_t d2 = 0.;
1438 fLeftMat->MasterToLocal(point, lpt);
1439 fRightMat->MasterToLocal(point, rpt);
1440 fLeftMat->MasterToLocalVect(dir, ldir);
1441 fRightMat->MasterToLocalVect(dir, rdir);
1442 Bool_t inleft = fLeft->Contains(lpt);
1443 Bool_t inright = fRight->Contains(rpt);
1444 node->SetSelected(0);
1445 Double_t snext = 0.0;
1446 if (inleft && inright) {
1447 // It is vey likely to have a numerical issue and the point should
1448 // be logically outside one of the shapes
1449 d1 = fLeft->DistFromInside(lpt,ldir,3);
1450 d2 = fRight->DistFromInside(rpt,rdir,3);
1451 if (d1<1.E-3) inleft = kFALSE;
1452 if (d2<1.E-3) inright = kFALSE;
1453 if (inleft && inright) return snext;
1454 }
1455
1456 while (1) {
1457 d1 = d2 = 0;
1458 if (!inleft) {
1459 d1 = fLeft->DistFromOutside(lpt,ldir,3);
1460 d1 = TMath::Max(d1,tol);
1461 if (d1>1E20) return TGeoShape::Big();
1462 }
1463 if (!inright) {
1464 d2 = fRight->DistFromOutside(rpt,rdir,3);
1465 d2 = TMath::Max(d2,tol);
1466 if (d2>1E20) return TGeoShape::Big();
1467 }
1468
1469 if (d1>d2) {
1470 // propagate to left shape
1471 snext += d1;
1472 node->SetSelected(1);
1473 inleft = kTRUE;
1474 for (i=0; i<3; i++) master[i] += d1*dir[i];
1475 fRightMat->MasterToLocal(master,rpt);
1476 // Push rpt to avoid a bad boundary condition
1477 for (i=0; i<3; i++) rpt[i] += tol*rdir[i];
1478 // check if propagated point is inside right shape
1479 inright = fRight->Contains(rpt);
1480 if (inright) return snext;
1481 // here inleft=true, inright=false
1482 } else {
1483 // propagate to right shape
1484 snext += d2;
1485 node->SetSelected(2);
1486 inright = kTRUE;
1487 for (i=0; i<3; i++) master[i] += d2*dir[i];
1488 fLeftMat->MasterToLocal(master,lpt);
1489 // Push lpt to avoid a bad boundary condition
1490 for (i=0; i<3; i++) lpt[i] += tol*ldir[i];
1491 // check if propagated point is inside left shape
1492 inleft = fLeft->Contains(lpt);
1493 if (inleft) return snext;
1494 // here inleft=false, inright=true
1495 }
1496 }
1497 return snext;
1498}
1499
1500////////////////////////////////////////////////////////////////////////////////
1501/// Returns number of vertices for the composite shape described by this intersection.
1502
1504{
1505 Int_t itot=0;
1506 Double_t point[3];
1507 Double_t tolerance = TGeoShape::Tolerance();
1508 if (fNpoints) return fNpoints;
1509 Int_t nleft = fLeft->GetNmeshVertices();
1510 Int_t nright = fRight->GetNmeshVertices();
1511 Double_t *points = new Double_t[3*(nleft+nright)];
1512 Double_t *points1 = new Double_t[3*nleft];
1513 fLeft->SetPoints(points1);
1514 for (Int_t i=0; i<nleft; i++) {
1515 if (TMath::Abs(points1[3*i])<tolerance && TMath::Abs(points1[3*i+1])<tolerance) continue;
1516 fLeftMat->LocalToMaster(&points1[3*i], &points[3*itot]);
1517 fRightMat->MasterToLocal(&points[3*itot], point);
1518 if (fRight->Contains(point)) itot++;
1519 }
1520 Double_t *points2 = new Double_t[3*nright];
1521 fRight->SetPoints(points2);
1522 for (Int_t i=0; i<nright; i++) {
1523 if (TMath::Abs(points2[3*i])<tolerance && TMath::Abs(points2[3*i+1])<tolerance) continue;
1524 fRightMat->LocalToMaster(&points2[3*i], &points[3*itot]);
1525 fLeftMat->MasterToLocal(&points[3*itot], point);
1526 if (fLeft->Contains(point)) itot++;
1527 }
1528
1529 AssignPoints(itot, points);
1530
1531 delete [] points1;
1532 delete [] points2;
1533 delete [] points;
1534 return fNpoints;
1535}
1536
1537////////////////////////////////////////////////////////////////////////////////
1538/// Compute safety distance for a union node;
1539
1541{
1542 Double_t local1[3], local2[3];
1543 fLeftMat->MasterToLocal(point,local1);
1544 Bool_t in1 = fLeft->Contains(local1);
1545 fRightMat->MasterToLocal(point,local2);
1546 Bool_t in2 = fRight->Contains(local2);
1547 Bool_t intrue = in1 & in2;
1548 if (in^intrue) return 0.0;
1549 Double_t saf1 = fLeft->Safety(local1, in1);
1550 Double_t saf2 = fRight->Safety(local2, in2);
1551 if (in1 && in2) return TMath::Min(saf1, saf2);
1552 if (in1) return saf2;
1553 if (in2) return saf1;
1554 return TMath::Max(saf1,saf2);
1555}
1556
1557////////////////////////////////////////////////////////////////////////////////
1558/// Save a primitive as a C++ statement(s) on output stream "out".
1559
1560void TGeoIntersection::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
1561{
1562 TGeoBoolNode::SavePrimitive(out,option);
1563 out << " pBoolNode = new TGeoIntersection(";
1564 out << fLeft->GetPointerName() << ",";
1565 out << fRight->GetPointerName() << ",";
1566 if (!fLeftMat->IsIdentity()) out << fLeftMat->GetPointerName() << ",";
1567 else out << "0,";
1568 if (!fRightMat->IsIdentity()) out << fRightMat->GetPointerName() << ");" << std::endl;
1569 else out << "0);" << std::endl;
1570}
1571
1572////////////////////////////////////////////////////////////////////////////////
1573/// Register 3D size of this shape.
1574
1576{
1578}
int Int_t
Definition: RtypesCore.h:45
const Bool_t kFALSE
Definition: RtypesCore.h:101
bool Bool_t
Definition: RtypesCore.h:63
double Double_t
Definition: RtypesCore.h:59
float Float_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:100
const char Option_t
Definition: RtypesCore.h:66
#define ClassImp(name)
Definition: Rtypes.h:364
R__EXTERN TGeoManager * gGeoManager
Definition: TGeoManager.h:602
R__EXTERN TGeoIdentity * gGeoIdentity
Definition: TGeoMatrix.h:478
float xmin
Definition: THbookFile.cxx:95
float ymin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
float ymax
Definition: THbookFile.cxx:95
#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:26
virtual void Sizeof3D() const
Register size of this 3D object.
Int_t fNpoints
Definition: TGeoBoolNode.h:52
Bool_t MakeBranch(const char *expr, Bool_t left)
Mutex for thread data access.
TGeoMatrix * fLeftMat
Definition: TGeoBoolNode.h:50
void ClearThreadData() const
TGeoShape * fLeft
Definition: TGeoBoolNode.h:48
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:56
std::mutex fMutex
Size for the navigation data array.
Definition: TGeoBoolNode.h:57
std::vector< ThreadData_t * > fThreadData
array of mesh points
Definition: TGeoBoolNode.h:55
void RegisterMatrices()
Register all matrices of the boolean node and descendents.
TGeoShape * fRight
Definition: TGeoBoolNode.h:49
virtual ~TGeoBoolNode()
Destructor.
virtual Int_t GetNpoints()=0
Double_t * fPoints
number of points on the mesh
Definition: TGeoBoolNode.h:53
TGeoBoolNode()
Default constructor.
TGeoMatrix * fRightMat
Definition: TGeoBoolNode.h:51
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:490
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:495
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:879
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save a primitive as a C++ statement(s) on output stream "out".
Definition: TObject.cxx:666
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:893
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:921
Basic string class.
Definition: TString.h:136
Ssiz_t Length() const
Definition: TString.h:410
const char * Data() const
Definition: TString.h:369
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:96
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:1168