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