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