Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
47 ClassImp(TGeoBoolNode);
48 
49 ////////////////////////////////////////////////////////////////////////////////
50 /// Constructor.
51 
52 TGeoBoolNode::ThreadData_t::ThreadData_t() :
53  fSelected(0)
54 {
55 }
56 
57 ////////////////////////////////////////////////////////////////////////////////
58 /// Destructor.
59 
60 TGeoBoolNode::ThreadData_t::~ThreadData_t()
61 {
62 }
63 
64 ////////////////////////////////////////////////////////////////////////////////
65 
66 TGeoBoolNode::ThreadData_t& TGeoBoolNode::GetThreadData() const
67 {
68  Int_t tid = TGeoManager::ThreadId();
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 
91 void TGeoBoolNode::ClearThreadData() const
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 
107 void TGeoBoolNode::CreateThreadData(Int_t nthreads)
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 
125 void TGeoBoolNode::SetSelected(Int_t sel)
126 {
127  GetThreadData().fSelected = sel;
128 }
129 
130 ////////////////////////////////////////////////////////////////////////////////
131 /// Default constructor
132 
133 TGeoBoolNode::TGeoBoolNode()
134 {
135  fLeft = 0;
136  fRight = 0;
137  fLeftMat = 0;
138  fRightMat = 0;
139  fNpoints = 0;
140  fPoints = 0;
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 = 0;
151  fRight = 0;
152  fLeftMat = 0;
153  fRightMat = 0;
154  fNpoints = 0;
155  fPoints = 0;
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 
169 TGeoBoolNode::TGeoBoolNode(TGeoShape *left, TGeoShape *right, TGeoMatrix *lmat, TGeoMatrix *rmat)
170 {
171  fLeft = left;
172  fRight = right;
173  fLeftMat = lmat;
174  fNpoints = 0;
175  fPoints = 0;
176  fThreadSize = 0;
177  CreateThreadData(1);
178  if (!fLeftMat) fLeftMat = gGeoIdentity;
179  else fLeftMat->RegisterYourself();
180  fRightMat = rmat;
181  if (!fRightMat) fRightMat = gGeoIdentity;
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 
197 TGeoBoolNode::~TGeoBoolNode()
198 {
199  if (fPoints) delete [] fPoints;
200  ClearThreadData();
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 
207 Bool_t TGeoBoolNode::MakeBranch(const char *expr, Bool_t left)
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 {
222  mat = (TGeoMatrix*)gGeoManager->GetListOfMatrices()->FindObject(stransf.Data());
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 
277 void TGeoBoolNode::Paint(Option_t * option)
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 
287  TGeoHMatrix *glmat = (TGeoHMatrix*)TGeoShape::GetTransform();
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 
320 void TGeoBoolNode::RegisterMatrices()
321 {
322  if (!fLeftMat->IsIdentity()) fLeftMat->RegisterYourself();
323  if (!fRightMat->IsIdentity()) fRightMat->RegisterYourself();
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 
332 Bool_t TGeoBoolNode::ReplaceMatrix(TGeoMatrix *mat, TGeoMatrix *newmat)
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 
357 void TGeoBoolNode::SavePrimitive(std::ostream &out, Option_t *option /*= ""*/)
358 {
359  fLeft->SavePrimitive(out,option);
360  fRight->SavePrimitive(out,option);
361  if (!fLeftMat->IsIdentity()) {
362  fLeftMat->RegisterYourself();
363  fLeftMat->SavePrimitive(out,option);
364  }
365  if (!fRightMat->IsIdentity()) {
366  fRightMat->RegisterYourself();
367  fRightMat->SavePrimitive(out,option);
368  }
369 }
370 
371 ////////////////////////////////////////////////////////////////////////////////
372 /// Fill buffer with shape vertices.
373 
374 void TGeoBoolNode::SetPoints(Double_t *points) const
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 
384 void TGeoBoolNode::SetPoints(Float_t *points) const
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 
394 void TGeoBoolNode::Sizeof3D() const
395 {
396  fLeft->Sizeof3D();
397  fRight->Sizeof3D();
398 }
399 
400 ClassImp(TGeoUnion);
401 
402 ////////////////////////////////////////////////////////////////////////////////
403 /// Make a clone of this. Pointers are preserved.
404 
405 TGeoBoolNode *TGeoUnion::MakeClone() const
406 {
407  return new TGeoUnion(fLeft, fRight, fLeftMat, fRightMat);
408 }
409 
410 ////////////////////////////////////////////////////////////////////////////////
411 /// Paint method.
412 
413 void TGeoUnion::Paint(Option_t *option)
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 
422  viewer->AddCompositeOp(TBuffer3D::kCSUnion);
423 
424  TGeoBoolNode::Paint(option);
425 }
426 
427 ////////////////////////////////////////////////////////////////////////////////
428 /// Default constructor
429 
430 TGeoUnion::TGeoUnion()
431 {
432 }
433 
434 ////////////////////////////////////////////////////////////////////////////////
435 /// Constructor
436 
437 TGeoUnion::TGeoUnion(const char *expr1, const char *expr2)
438  :TGeoBoolNode(expr1, expr2)
439 {
440 }
441 
442 ////////////////////////////////////////////////////////////////////////////////
443 /// Constructor providing pointers to components
444 
445 TGeoUnion::TGeoUnion(TGeoShape *left, TGeoShape *right, TGeoMatrix *lmat, TGeoMatrix *rmat)
446  :TGeoBoolNode(left,right,lmat,rmat)
447 {
448  if (left->TestShapeBit(TGeoShape::kGeoHalfSpace) || right->TestShapeBit(TGeoShape::kGeoHalfSpace)) {
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 
457 TGeoUnion::~TGeoUnion()
458 {
459 }
460 
461 ////////////////////////////////////////////////////////////////////////////////
462 /// Compute bounding box corresponding to a union of two shapes.
463 
464 void TGeoUnion::ComputeBBox(Double_t &dx, Double_t &dy, Double_t &dz, Double_t *origin)
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 
505 Bool_t TGeoUnion::Contains(const Double_t *point) const
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 
519 void TGeoUnion::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
520 {
521  ThreadData_t& td = GetThreadData();
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 
571 Int_t TGeoUnion::DistanceToPrimitive(Int_t /*px*/, Int_t /*py*/)
572 {
573  return 9999;
574 }
575 
576 ////////////////////////////////////////////////////////////////////////////////
577 /// Computes distance from a given point inside the shape to its boundary.
578 
579 Double_t TGeoUnion::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact,
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 
692 Double_t TGeoUnion::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact,
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 
723 Int_t TGeoUnion::GetNpoints()
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 
762 Double_t TGeoUnion::Safety(const Double_t *point, Bool_t in) const
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 
782 void 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 
797 void TGeoUnion::Sizeof3D() const
798 {
799  TGeoBoolNode::Sizeof3D();
800 }
801 
802 
803 ClassImp(TGeoSubtraction);
804 
805 ////////////////////////////////////////////////////////////////////////////////
806 /// Make a clone of this. Pointers are preserved.
807 
808 TGeoBoolNode *TGeoSubtraction::MakeClone() const
809 {
810  return new TGeoSubtraction(fLeft, fRight, fLeftMat, fRightMat);
811 }
812 
813 ////////////////////////////////////////////////////////////////////////////////
814 /// Paint method.
815 
816 void TGeoSubtraction::Paint(Option_t *option)
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 
825  viewer->AddCompositeOp(TBuffer3D::kCSDifference);
826 
827  TGeoBoolNode::Paint(option);
828 }
829 
830 ////////////////////////////////////////////////////////////////////////////////
831 /// Default constructor
832 
833 TGeoSubtraction::TGeoSubtraction()
834 {
835 }
836 
837 ////////////////////////////////////////////////////////////////////////////////
838 /// Constructor
839 
840 TGeoSubtraction::TGeoSubtraction(const char *expr1, const char *expr2)
841  :TGeoBoolNode(expr1, expr2)
842 {
843 }
844 
845 ////////////////////////////////////////////////////////////////////////////////
846 /// Constructor providing pointers to components
847 
848 TGeoSubtraction::TGeoSubtraction(TGeoShape *left, TGeoShape *right, TGeoMatrix *lmat, TGeoMatrix *rmat)
849  :TGeoBoolNode(left,right,lmat,rmat)
850 {
851  if (left->TestShapeBit(TGeoShape::kGeoHalfSpace)) {
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 
860 TGeoSubtraction::~TGeoSubtraction()
861 {
862 }
863 
864 ////////////////////////////////////////////////////////////////////////////////
865 /// Compute bounding box corresponding to a subtraction of two shapes.
866 
867 void TGeoSubtraction::ComputeBBox(Double_t &dx, Double_t &dy, Double_t &dz, Double_t *origin)
868 {
869  TGeoBBox *box = (TGeoBBox*)fLeft;
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 
898 void TGeoSubtraction::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
899 {
900  ThreadData_t& td = GetThreadData();
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 
948 Bool_t TGeoSubtraction::Contains(const Double_t *point) const
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 
962 Int_t TGeoSubtraction::DistanceToPrimitive(Int_t /*px*/, Int_t /*py*/)
963 {
964  return 9999;
965 }
966 
967 ////////////////////////////////////////////////////////////////////////////////
968 /// Compute distance from a given point inside to the shape boundary.
969 
970 Double_t TGeoSubtraction::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact,
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 
1001 Double_t TGeoSubtraction::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact,
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 
1058 Int_t TGeoSubtraction::GetNpoints()
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 
1095 Double_t TGeoSubtraction::Safety(const Double_t *point, Bool_t in) const
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 
1115 void 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 
1130 void TGeoSubtraction::Sizeof3D() const
1131 {
1132  TGeoBoolNode::Sizeof3D();
1133 }
1134 
1135 ClassImp(TGeoIntersection);
1136 
1137 ////////////////////////////////////////////////////////////////////////////////
1138 /// Make a clone of this. Pointers are preserved.
1139 
1140 TGeoBoolNode *TGeoIntersection::MakeClone() const
1141 {
1142  return new TGeoIntersection(fLeft, fRight, fLeftMat, fRightMat);
1143 }
1144 
1145 ////////////////////////////////////////////////////////////////////////////////
1146 /// Paint method.
1147 
1148 void TGeoIntersection::Paint(Option_t *option)
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 
1157  viewer->AddCompositeOp(TBuffer3D::kCSIntersection);
1158 
1159  TGeoBoolNode::Paint(option);
1160 }
1161 
1162 ////////////////////////////////////////////////////////////////////////////////
1163 /// Default constructor
1164 
1165 TGeoIntersection::TGeoIntersection()
1166 {
1167 }
1168 
1169 ////////////////////////////////////////////////////////////////////////////////
1170 /// Constructor
1171 
1172 TGeoIntersection::TGeoIntersection(const char *expr1, const char *expr2)
1173  :TGeoBoolNode(expr1, expr2)
1174 {
1175 }
1176 
1177 ////////////////////////////////////////////////////////////////////////////////
1178 /// Constructor providing pointers to components
1179 
1180 TGeoIntersection::TGeoIntersection(TGeoShape *left, TGeoShape *right, TGeoMatrix *lmat, TGeoMatrix *rmat)
1181  :TGeoBoolNode(left,right,lmat,rmat)
1182 {
1183  Bool_t hs1 = (fLeft->TestShapeBit(TGeoShape::kGeoHalfSpace))?kTRUE:kFALSE;
1184  Bool_t hs2 = (fRight->TestShapeBit(TGeoShape::kGeoHalfSpace))?kTRUE:kFALSE;
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 
1192 TGeoIntersection::~TGeoIntersection()
1193 {
1194 }
1195 
1196 ////////////////////////////////////////////////////////////////////////////////
1197 /// Compute bounding box corresponding to a intersection of two shapes.
1198 
1199 void TGeoIntersection::ComputeBBox(Double_t &dx, Double_t &dy, Double_t &dz, Double_t *origin)
1200 {
1201  Bool_t hs1 = (fLeft->TestShapeBit(TGeoShape::kGeoHalfSpace))?kTRUE:kFALSE;
1202  Bool_t hs2 = (fRight->TestShapeBit(TGeoShape::kGeoHalfSpace))?kTRUE:kFALSE;
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 
1300 void TGeoIntersection::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
1301 {
1302  ThreadData_t& td = GetThreadData();
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 
1350 Bool_t TGeoIntersection::Contains(const Double_t *point) const
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 
1364 Int_t TGeoIntersection::DistanceToPrimitive(Int_t /*px*/, Int_t /*py*/)
1365 {
1366  return 9999;
1367 }
1368 
1369 ////////////////////////////////////////////////////////////////////////////////
1370 /// Compute distance from a given point inside to the shape boundary.
1371 
1372 Double_t TGeoIntersection::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact,
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 
1403 Double_t TGeoIntersection::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact,
1404  Double_t step, Double_t *safe) const
1405 {
1406  Double_t tol = TGeoShape::Tolerance();
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 
1484 Int_t TGeoIntersection::GetNpoints()
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 
1521 Double_t TGeoIntersection::Safety(const Double_t *point, Bool_t in) const
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 
1541 void 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 
1556 void TGeoIntersection::Sizeof3D() const
1557 {
1558  TGeoBoolNode::Sizeof3D();
1559 }