Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TGeoBBox.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$// Author: Andrei Gheata 24/10/01
2 
3 // Contains() and DistFromOutside/Out() implemented by Mihaela Gheata
4 
5 /*************************************************************************
6  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
7  * All rights reserved. *
8  * *
9  * For the licensing terms see $ROOTSYS/LICENSE. *
10  * For the list of contributors see $ROOTSYS/README/CREDITS. *
11  *************************************************************************/
12 
13 /** \class TGeoBBox
14 \ingroup Geometry_classes
15 
16 Box class. All shape primitives inherit from this, their
17  constructor filling automatically the parameters of the box that bounds
18  the given shape. Defined by 6 parameters :
19  - fDX, fDY, fDZ - half lengths on X, Y and Z axis
20  - fOrigin[3] - position of box origin
21 
22 
23 ### Building boxes
24 
25  Normally a box has to be build only with 3 parameters : dx, dy, dz
26 representing the half lengths on X, Y and Z axis. In this case, the origin
27 of the box will match the one of its reference frame. The translation of the
28 origin is used only by the constructors of all other shapes in order to
29 define their own bounding boxes. Users should be aware that building a
30 translated box that will represent a physical shape by itself will affect any
31 further positioning of other shapes inside. Therefore in order to build a
32 positioned box one should follow the recipe described in class TGeoNode.
33 
34 #### Creation of boxes
35 
36  - TGeoBBox *box = new TGeoBBox("BOX", 20, 30, 40);
37 
38 Begin_Macro(source)
39 {
40  TCanvas *c = new TCanvas("c", "c",0,0,600,600);
41  new TGeoManager("box", "poza1");
42  TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
43  TGeoMedium *med = new TGeoMedium("MED",1,mat);
44  TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
45  gGeoManager->SetTopVolume(top);
46  TGeoVolume *vol = gGeoManager->MakeBox("BOX",med, 20,30,40);
47  vol->SetLineWidth(2);
48  top->AddNode(vol,1);
49  gGeoManager->CloseGeometry();
50  gGeoManager->SetNsegments(80);
51  top->Draw();
52  TView *view = gPad->GetView();
53  view->ShowAxis();
54 }
55 End_Macro
56 
57  - A volume having a box shape can be built in one step:
58  `TGeoVolume *vbox = gGeoManager->MakeBox("vbox", ptrMed, 20,30,40);`
59 
60 #### Divisions of boxes.
61 
62  Volumes having box shape can be divided with equal-length slices on
63 X, Y or Z axis. The following options are supported:
64 
65  - Dividing the full range of one axis in N slices
66  `TGeoVolume *divx = vbox->Divide("SLICEX", 1, N);`
67  - here 1 stands for the division axis (1-X, 2-Y, 3-Z)
68 
69 Begin_Macro(source)
70 {
71  TCanvas *c = new TCanvas("c", "c",0,0,600,600);
72  new TGeoManager("box", "poza1");
73  TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
74  TGeoMedium *med = new TGeoMedium("MED",1,mat);
75  TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
76  gGeoManager->SetTopVolume(top);
77  TGeoVolume *vol = gGeoManager->MakeBox("BOX",med, 20,30,40);
78  vol->SetLineWidth(2);
79  top->AddNode(vol,1);
80  TGeoVolume *divx = vol->Divide("SLICE",1,8,0,0);
81  gGeoManager->CloseGeometry();
82  gGeoManager->SetNsegments(80);
83  top->Draw();
84  TView *view = gPad->GetView();
85  view->ShowAxis();
86 }
87 End_Macro
88 
89  - Dividing in a limited range - general case.
90  `TGeoVolume *divy = vbox->Divide("SLICEY",2,N,start,step);`
91  - start = starting offset within (-fDY, fDY)
92  - step = slicing step
93 
94 Begin_Macro(source)
95 {
96  TCanvas *c = new TCanvas("c", "c",0,0,600,600);
97  new TGeoManager("box", "poza1");
98  TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
99  TGeoMedium *med = new TGeoMedium("MED",1,mat);
100  TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
101  gGeoManager->SetTopVolume(top);
102  TGeoVolume *vol = gGeoManager->MakeBox("BOX",med, 20,30,40);
103  vol->SetLineWidth(2);
104  top->AddNode(vol,1);
105  TGeoVolume *divx = vol->Divide("SLICE",2,8,2,3);
106  gGeoManager->CloseGeometry();
107  gGeoManager->SetNsegments(80);
108  top->Draw();
109  TView *view = gPad->GetView();
110  view->ShowAxis();
111 }
112 End_Macro
113 
114 Both cases are supported by all shapes.
115 See also class TGeoShape for utility methods provided by any particular shape.
116 */
117 
118 #include "Riostream.h"
119 
120 #include "TGeoManager.h"
121 #include "TGeoMatrix.h"
122 #include "TGeoVolume.h"
123 #include "TVirtualGeoPainter.h"
124 #include "TGeoBBox.h"
125 #include "TVirtualPad.h"
126 #include "TBuffer3D.h"
127 #include "TBuffer3DTypes.h"
128 #include "TMath.h"
129 #include "TRandom.h"
130 
131 ClassImp(TGeoBBox);
132 
133 ////////////////////////////////////////////////////////////////////////////////
134 /// Default constructor
135 
136 TGeoBBox::TGeoBBox()
137 {
138  SetShapeBit(TGeoShape::kGeoBox);
139  fDX = fDY = fDZ = 0;
140  fOrigin[0] = fOrigin[1] = fOrigin[2] = 0.0;
141 }
142 
143 ////////////////////////////////////////////////////////////////////////////////
144 /// Constructor where half-lengths are provided.
145 
146 TGeoBBox::TGeoBBox(Double_t dx, Double_t dy, Double_t dz, Double_t *origin)
147  :TGeoShape("")
148 {
149  SetShapeBit(TGeoShape::kGeoBox);
150  fOrigin[0] = fOrigin[1] = fOrigin[2] = 0.0;
151  SetBoxDimensions(dx, dy, dz, origin);
152 }
153 
154 ////////////////////////////////////////////////////////////////////////////////
155 /// Constructor with shape name.
156 
157 TGeoBBox::TGeoBBox(const char *name, Double_t dx, Double_t dy, Double_t dz, Double_t *origin)
158  :TGeoShape(name)
159 {
160  SetShapeBit(TGeoShape::kGeoBox);
161  fOrigin[0] = fOrigin[1] = fOrigin[2] = 0.0;
162  SetBoxDimensions(dx, dy, dz, origin);
163 }
164 
165 ////////////////////////////////////////////////////////////////////////////////
166 /// Constructor based on the array of parameters.
167 /// - param[0] - half-length in x
168 /// - param[1] - half-length in y
169 /// - param[2] - half-length in z
170 
171 TGeoBBox::TGeoBBox(Double_t *param)
172  :TGeoShape("")
173 {
174  SetShapeBit(TGeoShape::kGeoBox);
175  fOrigin[0] = fOrigin[1] = fOrigin[2] = 0.0;
176  SetDimensions(param);
177 }
178 
179 ////////////////////////////////////////////////////////////////////////////////
180 /// Destructor
181 
182 TGeoBBox::~TGeoBBox()
183 {
184 }
185 
186 ////////////////////////////////////////////////////////////////////////////////
187 /// Check if 2 positioned boxes overlap.
188 
189 Bool_t TGeoBBox::AreOverlapping(const TGeoBBox *box1, const TGeoMatrix *mat1, const TGeoBBox *box2, const TGeoMatrix *mat2)
190 {
191  Double_t master[3];
192  Double_t local[3];
193  Double_t ldir1[3], ldir2[3];
194  const Double_t *o1 = box1->GetOrigin();
195  const Double_t *o2 = box2->GetOrigin();
196  // Convert center of first box to the local frame of second
197  mat1->LocalToMaster(o1, master);
198  mat2->MasterToLocal(master, local);
199  if (TGeoBBox::Contains(local,box2->GetDX(),box2->GetDY(),box2->GetDZ(),o2)) return kTRUE;
200  Double_t distsq = (local[0]-o2[0])*(local[0]-o2[0]) +
201  (local[1]-o2[1])*(local[1]-o2[1]) +
202  (local[2]-o2[2])*(local[2]-o2[2]);
203  // Compute distance between box centers and compare with max value
204  Double_t rmaxsq = (box1->GetDX()+box2->GetDX())*(box1->GetDX()+box2->GetDX()) +
205  (box1->GetDY()+box2->GetDY())*(box1->GetDY()+box2->GetDY()) +
206  (box1->GetDZ()+box2->GetDZ())*(box1->GetDZ()+box2->GetDZ());
207  if (distsq > rmaxsq + TGeoShape::Tolerance()) return kFALSE;
208  // We are still not sure: shoot a ray from the center of "1" towards the
209  // center of 2.
210  Double_t dir[3];
211  mat1->LocalToMaster(o1, ldir1);
212  mat2->LocalToMaster(o2, ldir2);
213  distsq = 1./TMath::Sqrt(distsq);
214  dir[0] = (ldir2[0]-ldir1[0])*distsq;
215  dir[1] = (ldir2[1]-ldir1[1])*distsq;
216  dir[2] = (ldir2[2]-ldir1[2])*distsq;
217  mat1->MasterToLocalVect(dir, ldir1);
218  mat2->MasterToLocalVect(dir, ldir2);
219  // Distance to exit from o1
220  Double_t dist1 = TGeoBBox::DistFromInside(o1,ldir1,box1->GetDX(),box1->GetDY(),box1->GetDZ(),o1);
221  // Distance to enter from o2
222  Double_t dist2 = TGeoBBox::DistFromOutside(local,ldir2,box2->GetDX(),box2->GetDY(),box2->GetDZ(),o2);
223  if (dist1 > dist2) return kTRUE;
224  return kFALSE;
225 }
226 
227 ////////////////////////////////////////////////////////////////////////////////
228 /// Computes capacity of the shape in [length^3].
229 
230 Double_t TGeoBBox::Capacity() const
231 {
232  return (8.*fDX*fDY*fDZ);
233 }
234 
235 ////////////////////////////////////////////////////////////////////////////////
236 /// Computes normal to closest surface from POINT.
237 
238 void TGeoBBox::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
239 {
240  memset(norm,0,3*sizeof(Double_t));
241  Double_t saf[3];
242  Int_t i;
243  saf[0]=TMath::Abs(TMath::Abs(point[0]-fOrigin[0])-fDX);
244  saf[1]=TMath::Abs(TMath::Abs(point[1]-fOrigin[1])-fDY);
245  saf[2]=TMath::Abs(TMath::Abs(point[2]-fOrigin[2])-fDZ);
246  i = TMath::LocMin(3,saf);
247  norm[i] = (dir[i]>0)?1:(-1);
248 }
249 
250 ////////////////////////////////////////////////////////////////////////////////
251 /// Decides fast if the bounding box could be crossed by a vector.
252 
253 Bool_t TGeoBBox::CouldBeCrossed(const Double_t *point, const Double_t *dir) const
254 {
255  Double_t mind = fDX;
256  if (fDY<mind) mind=fDY;
257  if (fDZ<mind) mind=fDZ;
258  Double_t dx = fOrigin[0]-point[0];
259  Double_t dy = fOrigin[1]-point[1];
260  Double_t dz = fOrigin[2]-point[2];
261  Double_t do2 = dx*dx+dy*dy+dz*dz;
262  if (do2<=(mind*mind)) return kTRUE;
263  Double_t rmax2 = fDX*fDX+fDY*fDY+fDZ*fDZ;
264  if (do2<=rmax2) return kTRUE;
265  // inside bounding sphere
266  Double_t doct = dx*dir[0]+dy*dir[1]+dz*dir[2];
267  // leaving ray
268  if (doct<=0) return kFALSE;
269  Double_t dirnorm=dir[0]*dir[0]+dir[1]*dir[1]+dir[2]*dir[2];
270  if ((doct*doct)>=(do2-rmax2)*dirnorm) return kTRUE;
271  return kFALSE;
272 }
273 
274 ////////////////////////////////////////////////////////////////////////////////
275 /// Compute closest distance from point px,py to each corner.
276 
277 Int_t TGeoBBox::DistancetoPrimitive(Int_t px, Int_t py)
278 {
279  const Int_t numPoints = 8;
280  return ShapeDistancetoPrimitive(numPoints, px, py);
281 }
282 
283 ////////////////////////////////////////////////////////////////////////////////
284 /// Divide this box shape belonging to volume "voldiv" into ndiv equal volumes
285 /// called divname, from start position with the given step. Returns pointer
286 /// to created division cell volume. In case a wrong division axis is supplied,
287 /// returns pointer to volume to be divided.
288 
289 TGeoVolume *TGeoBBox::Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv,
290  Double_t start, Double_t step)
291 {
292  TGeoShape *shape; //--- shape to be created
293  TGeoVolume *vol; //--- division volume to be created
294  TGeoVolumeMulti *vmulti; //--- generic divided volume
295  TGeoPatternFinder *finder; //--- finder to be attached
296  TString opt = ""; //--- option to be attached
297  Double_t end = start+ndiv*step;
298  switch (iaxis) {
299  case 1: //--- divide on X
300  shape = new TGeoBBox(step/2., fDY, fDZ);
301  finder = new TGeoPatternX(voldiv, ndiv, start, end);
302  opt = "X";
303  break;
304  case 2: //--- divide on Y
305  shape = new TGeoBBox(fDX, step/2., fDZ);
306  finder = new TGeoPatternY(voldiv, ndiv, start, end);
307  opt = "Y";
308  break;
309  case 3: //--- divide on Z
310  shape = new TGeoBBox(fDX, fDY, step/2.);
311  finder = new TGeoPatternZ(voldiv, ndiv, start, end);
312  opt = "Z";
313  break;
314  default:
315  Error("Divide", "Wrong axis type for division");
316  return 0;
317  }
318  vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
319  vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
320  vmulti->AddVolume(vol);
321  voldiv->SetFinder(finder);
322  finder->SetDivIndex(voldiv->GetNdaughters());
323  for (Int_t ic=0; ic<ndiv; ic++) {
324  voldiv->AddNodeOffset(vol, ic, start+step/2.+ic*step, opt.Data());
325  ((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder);
326  }
327  return vmulti;
328 }
329 
330 ////////////////////////////////////////////////////////////////////////////////
331 /// Compute bounding box - nothing to do in this case.
332 
333 void TGeoBBox::ComputeBBox()
334 {
335 }
336 
337 ////////////////////////////////////////////////////////////////////////////////
338 /// Test if point is inside this shape.
339 
340 Bool_t TGeoBBox::Contains(const Double_t *point) const
341 {
342  if (TMath::Abs(point[2]-fOrigin[2]) > fDZ) return kFALSE;
343  if (TMath::Abs(point[0]-fOrigin[0]) > fDX) return kFALSE;
344  if (TMath::Abs(point[1]-fOrigin[1]) > fDY) return kFALSE;
345  return kTRUE;
346 }
347 
348 ////////////////////////////////////////////////////////////////////////////////
349 /// Static method to check if point[3] is located inside a box of having dx, dy, dz
350 /// as half-lengths.
351 
352 Bool_t TGeoBBox::Contains(const Double_t *point, Double_t dx, Double_t dy, Double_t dz, const Double_t *origin)
353 {
354  if (TMath::Abs(point[2]-origin[2]) > dz) return kFALSE;
355  if (TMath::Abs(point[0]-origin[0]) > dx) return kFALSE;
356  if (TMath::Abs(point[1]-origin[1]) > dy) return kFALSE;
357  return kTRUE;
358 }
359 
360 ////////////////////////////////////////////////////////////////////////////////
361 /// Compute distance from inside point to surface of the box.
362 /// Boundary safe algorithm.
363 
364 Double_t TGeoBBox::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
365 {
366  Double_t s,smin,saf[6];
367  Double_t newpt[3];
368  Int_t i;
369  for (i=0; i<3; i++) newpt[i] = point[i] - fOrigin[i];
370  saf[0] = fDX+newpt[0];
371  saf[1] = fDX-newpt[0];
372  saf[2] = fDY+newpt[1];
373  saf[3] = fDY-newpt[1];
374  saf[4] = fDZ+newpt[2];
375  saf[5] = fDZ-newpt[2];
376  if (iact<3 && safe) {
377  smin = saf[0];
378  // compute safe distance
379  for (i=1;i<6;i++) if (saf[i] < smin) smin = saf[i];
380  *safe = smin;
381  if (smin<0) *safe = 0.0;
382  if (iact==0) return TGeoShape::Big();
383  if (iact==1 && step<*safe) return TGeoShape::Big();
384  }
385  // compute distance to surface
386  smin=TGeoShape::Big();
387  for (i=0; i<3; i++) {
388  if (dir[i]!=0) {
389  s = (dir[i]>0)?(saf[(i<<1)+1]/dir[i]):(-saf[i<<1]/dir[i]);
390  if (s < 0) return 0.0;
391  if (s < smin) smin = s;
392  }
393  }
394  return smin;
395 }
396 
397 ////////////////////////////////////////////////////////////////////////////////
398 /// Compute distance from inside point to surface of the box.
399 /// Boundary safe algorithm.
400 
401 Double_t TGeoBBox::DistFromInside(const Double_t *point,const Double_t *dir,
402  Double_t dx, Double_t dy, Double_t dz, const Double_t *origin, Double_t /*stepmax*/)
403 {
404  Double_t s,smin,saf[6];
405  Double_t newpt[3];
406  Int_t i;
407  for (i=0; i<3; i++) newpt[i] = point[i] - origin[i];
408  saf[0] = dx+newpt[0];
409  saf[1] = dx-newpt[0];
410  saf[2] = dy+newpt[1];
411  saf[3] = dy-newpt[1];
412  saf[4] = dz+newpt[2];
413  saf[5] = dz-newpt[2];
414  // compute distance to surface
415  smin=TGeoShape::Big();
416  for (i=0; i<3; i++) {
417  if (dir[i]!=0) {
418  s = (dir[i]>0)?(saf[(i<<1)+1]/dir[i]):(-saf[i<<1]/dir[i]);
419  if (s < 0) return 0.0;
420  if (s < smin) smin = s;
421  }
422  }
423  return smin;
424 }
425 
426 ////////////////////////////////////////////////////////////////////////////////
427 /// Compute distance from outside point to surface of the box.
428 /// Boundary safe algorithm.
429 
430 Double_t TGeoBBox::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
431 {
432  Bool_t in = kTRUE;
433  Double_t saf[3];
434  Double_t par[3];
435  Double_t newpt[3];
436  Int_t i,j;
437  for (i=0; i<3; i++) newpt[i] = point[i] - fOrigin[i];
438  par[0] = fDX;
439  par[1] = fDY;
440  par[2] = fDZ;
441  for (i=0; i<3; i++) {
442  saf[i] = TMath::Abs(newpt[i])-par[i];
443  if (saf[i]>=step) return TGeoShape::Big();
444  if (in && saf[i]>0) in=kFALSE;
445  }
446  if (iact<3 && safe) {
447  // compute safe distance
448  if (in) {
449  *safe = 0.0;
450  } else {
451  *safe = saf[0];
452  if (saf[1] > *safe) *safe = saf[1];
453  if (saf[2] > *safe) *safe = saf[2];
454  }
455  if (iact==0) return TGeoShape::Big();
456  if (iact==1 && step<*safe) return TGeoShape::Big();
457  }
458  // compute distance from point to box
459  Double_t coord, snxt=TGeoShape::Big();
460  Int_t ibreak=0;
461  // protection in case point is actually inside box
462  if (in) {
463  j = 0;
464  Double_t ss = saf[0];
465  if (saf[1]>ss) {
466  ss = saf[1];
467  j = 1;
468  }
469  if (saf[2]>ss) j = 2;
470  if (newpt[j]*dir[j]>0) return TGeoShape::Big(); // in fact exiting
471  return 0.0;
472  }
473  for (i=0; i<3; i++) {
474  if (saf[i]<0) continue;
475  if (newpt[i]*dir[i] >= 0) continue;
476  snxt = saf[i]/TMath::Abs(dir[i]);
477  ibreak = 0;
478  for (j=0; j<3; j++) {
479  if (j==i) continue;
480  coord=newpt[j]+snxt*dir[j];
481  if (TMath::Abs(coord)>par[j]) {
482  ibreak=1;
483  break;
484  }
485  }
486  if (!ibreak) return snxt;
487  }
488  return TGeoShape::Big();
489 }
490 
491 ////////////////////////////////////////////////////////////////////////////////
492 /// Compute distance from outside point to surface of the box.
493 /// Boundary safe algorithm.
494 
495 Double_t TGeoBBox::DistFromOutside(const Double_t *point,const Double_t *dir,
496  Double_t dx, Double_t dy, Double_t dz, const Double_t *origin, Double_t stepmax)
497 {
498  Bool_t in = kTRUE;
499  Double_t saf[3];
500  Double_t par[3];
501  Double_t newpt[3];
502  Int_t i,j;
503  for (i=0; i<3; i++) newpt[i] = point[i] - origin[i];
504  par[0] = dx;
505  par[1] = dy;
506  par[2] = dz;
507  for (i=0; i<3; i++) {
508  saf[i] = TMath::Abs(newpt[i])-par[i];
509  if (saf[i]>=stepmax) return TGeoShape::Big();
510  if (in && saf[i]>0) in=kFALSE;
511  }
512  // In case point is inside return ZERO
513  if (in) return 0.0;
514  Double_t coord, snxt=TGeoShape::Big();
515  Int_t ibreak=0;
516  for (i=0; i<3; i++) {
517  if (saf[i]<0) continue;
518  if (newpt[i]*dir[i] >= 0) continue;
519  snxt = saf[i]/TMath::Abs(dir[i]);
520  ibreak = 0;
521  for (j=0; j<3; j++) {
522  if (j==i) continue;
523  coord=newpt[j]+snxt*dir[j];
524  if (TMath::Abs(coord)>par[j]) {
525  ibreak=1;
526  break;
527  }
528  }
529  if (!ibreak) return snxt;
530  }
531  return TGeoShape::Big();
532 }
533 
534 ////////////////////////////////////////////////////////////////////////////////
535 /// Returns name of axis IAXIS.
536 
537 const char *TGeoBBox::GetAxisName(Int_t iaxis) const
538 {
539  switch (iaxis) {
540  case 1:
541  return "X";
542  case 2:
543  return "Y";
544  case 3:
545  return "Z";
546  default:
547  return "UNDEFINED";
548  }
549 }
550 
551 ////////////////////////////////////////////////////////////////////////////////
552 /// Get range of shape for a given axis.
553 
554 Double_t TGeoBBox::GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const
555 {
556  xlo = 0;
557  xhi = 0;
558  Double_t dx = 0;
559  switch (iaxis) {
560  case 1:
561  xlo = fOrigin[0]-fDX;
562  xhi = fOrigin[0]+fDX;
563  dx = 2*fDX;
564  return dx;
565  case 2:
566  xlo = fOrigin[1]-fDY;
567  xhi = fOrigin[1]+fDY;
568  dx = 2*fDY;
569  return dx;
570  case 3:
571  xlo = fOrigin[2]-fDZ;
572  xhi = fOrigin[2]+fDZ;
573  dx = 2*fDZ;
574  return dx;
575  }
576  return dx;
577 }
578 
579 ////////////////////////////////////////////////////////////////////////////////
580 /// Fill vector param[4] with the bounding cylinder parameters. The order
581 /// is the following : Rmin, Rmax, Phi1, Phi2
582 
583 void TGeoBBox::GetBoundingCylinder(Double_t *param) const
584 {
585  param[0] = 0.; // Rmin
586  param[1] = fDX*fDX+fDY*fDY; // Rmax
587  param[2] = 0.; // Phi1
588  param[3] = 360.; // Phi2
589 }
590 
591 ////////////////////////////////////////////////////////////////////////////////
592 /// Get area in internal units of the facet with a given index.
593 /// Possible index values:
594 /// - 0 - all facets together
595 /// - 1 to 6 - facet index from bottom to top Z
596 
597 Double_t TGeoBBox::GetFacetArea(Int_t index) const
598 {
599  Double_t area = 0.;
600  switch (index) {
601  case 0:
602  area = 8.*(fDX*fDY + fDX*fDZ + fDY*fDZ);
603  return area;
604  case 1:
605  case 6:
606  area = 4.*fDX*fDY;
607  return area;
608  case 2:
609  case 4:
610  area = 4.*fDX*fDZ;
611  return area;
612  case 3:
613  case 5:
614  area = 4.*fDY*fDZ;
615  return area;
616  }
617  return area;
618 }
619 
620 ////////////////////////////////////////////////////////////////////////////////
621 /// Fills array with n random points located on the surface of indexed facet.
622 /// The output array must be provided with a length of minimum 3*npoints. Returns
623 /// true if operation succeeded.
624 /// Possible index values:
625 /// - 0 - all facets together
626 /// - 1 to 6 - facet index from bottom to top Z
627 
628 Bool_t TGeoBBox::GetPointsOnFacet(Int_t index, Int_t npoints, Double_t *array) const
629 {
630  if (index<0 || index>6) return kFALSE;
631  Double_t surf[6];
632  Double_t area = 0.;
633  if (index==0) {
634  for (Int_t isurf=0; isurf<6; isurf++) {
635  surf[isurf] = TGeoBBox::GetFacetArea(isurf+1);
636  if (isurf>0) surf[isurf] += surf[isurf-1];
637  }
638  area = surf[5];
639  }
640 
641  for (Int_t i=0; i<npoints; i++) {
642  // Generate randomly a surface index if needed.
643  Double_t *point = &array[3*i];
644  Int_t surfindex = index;
645  if (surfindex==0) {
646  Double_t val = area*gRandom->Rndm();
647  surfindex = 2+TMath::BinarySearch(6, surf, val);
648  if (surfindex>6) surfindex=6;
649  }
650  switch (surfindex) {
651  case 1:
652  point[0] = -fDX + 2*fDX*gRandom->Rndm();
653  point[1] = -fDY + 2*fDY*gRandom->Rndm();
654  point[2] = -fDZ;
655  break;
656  case 2:
657  point[0] = -fDX + 2*fDX*gRandom->Rndm();
658  point[1] = -fDY;
659  point[2] = -fDZ + 2*fDZ*gRandom->Rndm();
660  break;
661  case 3:
662  point[0] = -fDX;
663  point[1] = -fDY + 2*fDY*gRandom->Rndm();
664  point[2] = -fDZ + 2*fDZ*gRandom->Rndm();
665  break;
666  case 4:
667  point[0] = -fDX + 2*fDX*gRandom->Rndm();
668  point[1] = fDY;
669  point[2] = -fDZ + 2*fDZ*gRandom->Rndm();
670  break;
671  case 5:
672  point[0] = fDX;
673  point[1] = -fDY + 2*fDY*gRandom->Rndm();
674  point[2] = -fDZ + 2*fDZ*gRandom->Rndm();
675  break;
676  case 6:
677  point[0] = -fDX + 2*fDX*gRandom->Rndm();
678  point[1] = -fDY + 2*fDY*gRandom->Rndm();
679  point[2] = fDZ;
680  break;
681  }
682  }
683  return kTRUE;
684 }
685 
686 ////////////////////////////////////////////////////////////////////////////////
687 /// Fills array with n random points located on the line segments of the shape mesh.
688 /// The output array must be provided with a length of minimum 3*npoints. Returns
689 /// true if operation is implemented.
690 
691 Bool_t TGeoBBox::GetPointsOnSegments(Int_t npoints, Double_t *array) const
692 {
693  if (npoints<GetNmeshVertices()) {
694  Error("GetPointsOnSegments", "You should require at least %d points", GetNmeshVertices());
695  return kFALSE;
696  }
697  TBuffer3D &buff = (TBuffer3D &)GetBuffer3D(TBuffer3D::kRawSizes | TBuffer3D::kRaw, kTRUE);
698  Int_t npnts = buff.NbPnts();
699  Int_t nsegs = buff.NbSegs();
700  // Copy buffered points in the array
701  memcpy(array, buff.fPnts, 3*npnts*sizeof(Double_t));
702  Int_t ipoints = npoints - npnts;
703  Int_t icrt = 3*npnts;
704  Int_t nperseg = (Int_t)(Double_t(ipoints)/nsegs);
705  Double_t *p0, *p1;
706  Double_t x,y,z, dx,dy,dz;
707  for (Int_t i=0; i<nsegs; i++) {
708  p0 = &array[3*buff.fSegs[3*i+1]];
709  p1 = &array[3*buff.fSegs[3*i+2]];
710  if (i==(nsegs-1)) nperseg = ipoints;
711  dx = (p1[0]-p0[0])/(nperseg+1);
712  dy = (p1[1]-p0[1])/(nperseg+1);
713  dz = (p1[2]-p0[2])/(nperseg+1);
714  for (Int_t j=0; j<nperseg; j++) {
715  x = p0[0] + (j+1)*dx;
716  y = p0[1] + (j+1)*dy;
717  z = p0[2] + (j+1)*dz;
718  array[icrt++] = x; array[icrt++] = y; array[icrt++] = z;
719  ipoints--;
720  }
721  }
722  return kTRUE;
723 }
724 
725 ////////////////////////////////////////////////////////////////////////////////
726 /// Fills real parameters of a positioned box inside this one. Returns 0 if successful.
727 
728 Int_t TGeoBBox::GetFittingBox(const TGeoBBox *parambox, TGeoMatrix *mat, Double_t &dx, Double_t &dy, Double_t &dz) const
729 {
730  dx=dy=dz=0;
731  if (mat->IsRotation()) {
732  Error("GetFittingBox", "cannot handle parametrized rotated volumes");
733  return 1; // ### rotation not accepted ###
734  }
735  //--> translate the origin of the parametrized box to the frame of this box.
736  Double_t origin[3];
737  mat->LocalToMaster(parambox->GetOrigin(), origin);
738  if (!Contains(origin)) {
739  Error("GetFittingBox", "wrong matrix - parametrized box is outside this");
740  return 1; // ### wrong matrix ###
741  }
742  //--> now we have to get the valid range for all parametrized axis
743  Double_t xlo=0, xhi=0;
744  Double_t dd[3];
745  dd[0] = parambox->GetDX();
746  dd[1] = parambox->GetDY();
747  dd[2] = parambox->GetDZ();
748  for (Int_t iaxis=0; iaxis<3; iaxis++) {
749  if (dd[iaxis]>=0) continue;
750  TGeoBBox::GetAxisRange(iaxis+1, xlo, xhi);
751  //-> compute best fitting parameter
752  dd[iaxis] = TMath::Min(origin[iaxis]-xlo, xhi-origin[iaxis]);
753  if (dd[iaxis]<0) {
754  Error("GetFittingBox", "wrong matrix");
755  return 1;
756  }
757  }
758  dx = dd[0];
759  dy = dd[1];
760  dz = dd[2];
761  return 0;
762 }
763 
764 ////////////////////////////////////////////////////////////////////////////////
765 /// In case shape has some negative parameters, these has to be computed
766 /// in order to fit the mother
767 
768 TGeoShape *TGeoBBox::GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix *mat) const
769 {
770  if (!TestShapeBit(kGeoRunTimeShape)) return 0;
771  Double_t dx, dy, dz;
772  Int_t ierr = mother->GetFittingBox(this, mat, dx, dy, dz);
773  if (ierr) {
774  Error("GetMakeRuntimeShape", "cannot fit this to mother");
775  return 0;
776  }
777  return (new TGeoBBox(dx, dy, dz));
778 }
779 
780 ////////////////////////////////////////////////////////////////////////////////
781 /// Returns numbers of vertices, segments and polygons composing the shape mesh.
782 
783 void TGeoBBox::GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const
784 {
785  nvert = 8;
786  nsegs = 12;
787  npols = 6;
788 }
789 
790 ////////////////////////////////////////////////////////////////////////////////
791 /// Prints shape parameters
792 
793 void TGeoBBox::InspectShape() const
794 {
795  printf("*** Shape %s: TGeoBBox ***\n", GetName());
796  printf(" dX = %11.5f\n", fDX);
797  printf(" dY = %11.5f\n", fDY);
798  printf(" dZ = %11.5f\n", fDZ);
799  printf(" origin: x=%11.5f y=%11.5f z=%11.5f\n", fOrigin[0], fOrigin[1], fOrigin[2]);
800 }
801 
802 ////////////////////////////////////////////////////////////////////////////////
803 /// Creates a TBuffer3D describing *this* shape.
804 /// Coordinates are in local reference frame.
805 
806 TBuffer3D *TGeoBBox::MakeBuffer3D() const
807 {
808  TBuffer3D* buff = new TBuffer3D(TBuffer3DTypes::kGeneric, 8, 24, 12, 36, 6, 36);
809  if (buff)
810  {
811  SetPoints(buff->fPnts);
812  SetSegsAndPols(*buff);
813  }
814 
815  return buff;
816 }
817 
818 ////////////////////////////////////////////////////////////////////////////////
819 /// Fills TBuffer3D structure for segments and polygons.
820 
821 void TGeoBBox::SetSegsAndPols(TBuffer3D &buff) const
822 {
823  Int_t c = GetBasicColor();
824 
825  buff.fSegs[ 0] = c ; buff.fSegs[ 1] = 0 ; buff.fSegs[ 2] = 1 ;
826  buff.fSegs[ 3] = c+1 ; buff.fSegs[ 4] = 1 ; buff.fSegs[ 5] = 2 ;
827  buff.fSegs[ 6] = c+1 ; buff.fSegs[ 7] = 2 ; buff.fSegs[ 8] = 3 ;
828  buff.fSegs[ 9] = c ; buff.fSegs[10] = 3 ; buff.fSegs[11] = 0 ;
829  buff.fSegs[12] = c+2 ; buff.fSegs[13] = 4 ; buff.fSegs[14] = 5 ;
830  buff.fSegs[15] = c+2 ; buff.fSegs[16] = 5 ; buff.fSegs[17] = 6 ;
831  buff.fSegs[18] = c+3 ; buff.fSegs[19] = 6 ; buff.fSegs[20] = 7 ;
832  buff.fSegs[21] = c+3 ; buff.fSegs[22] = 7 ; buff.fSegs[23] = 4 ;
833  buff.fSegs[24] = c ; buff.fSegs[25] = 0 ; buff.fSegs[26] = 4 ;
834  buff.fSegs[27] = c+2 ; buff.fSegs[28] = 1 ; buff.fSegs[29] = 5 ;
835  buff.fSegs[30] = c+1 ; buff.fSegs[31] = 2 ; buff.fSegs[32] = 6 ;
836  buff.fSegs[33] = c+3 ; buff.fSegs[34] = 3 ; buff.fSegs[35] = 7 ;
837 
838  buff.fPols[ 0] = c ; buff.fPols[ 1] = 4 ; buff.fPols[ 2] = 0 ;
839  buff.fPols[ 3] = 9 ; buff.fPols[ 4] = 4 ; buff.fPols[ 5] = 8 ;
840  buff.fPols[ 6] = c+1 ; buff.fPols[ 7] = 4 ; buff.fPols[ 8] = 1 ;
841  buff.fPols[ 9] = 10 ; buff.fPols[10] = 5 ; buff.fPols[11] = 9 ;
842  buff.fPols[12] = c ; buff.fPols[13] = 4 ; buff.fPols[14] = 2 ;
843  buff.fPols[15] = 11 ; buff.fPols[16] = 6 ; buff.fPols[17] = 10 ;
844  buff.fPols[18] = c+1 ; buff.fPols[19] = 4 ; buff.fPols[20] = 3 ;
845  buff.fPols[21] = 8 ; buff.fPols[22] = 7 ; buff.fPols[23] = 11 ;
846  buff.fPols[24] = c+2 ; buff.fPols[25] = 4 ; buff.fPols[26] = 0 ;
847  buff.fPols[27] = 3 ; buff.fPols[28] = 2 ; buff.fPols[29] = 1 ;
848  buff.fPols[30] = c+3 ; buff.fPols[31] = 4 ; buff.fPols[32] = 4 ;
849  buff.fPols[33] = 5 ; buff.fPols[34] = 6 ; buff.fPols[35] = 7 ;
850 }
851 
852 ////////////////////////////////////////////////////////////////////////////////
853 /// Computes the closest distance from given point to this shape.
854 
855 Double_t TGeoBBox::Safety(const Double_t *point, Bool_t in) const
856 {
857  Double_t safe, safy, safz;
858  if (in) {
859  safe = fDX - TMath::Abs(point[0]-fOrigin[0]);
860  safy = fDY - TMath::Abs(point[1]-fOrigin[1]);
861  safz = fDZ - TMath::Abs(point[2]-fOrigin[2]);
862  if (safy < safe) safe = safy;
863  if (safz < safe) safe = safz;
864  } else {
865  safe = -fDX + TMath::Abs(point[0]-fOrigin[0]);
866  safy = -fDY + TMath::Abs(point[1]-fOrigin[1]);
867  safz = -fDZ + TMath::Abs(point[2]-fOrigin[2]);
868  if (safy > safe) safe = safy;
869  if (safz > safe) safe = safz;
870  }
871  return safe;
872 }
873 
874 ////////////////////////////////////////////////////////////////////////////////
875 /// Save a primitive as a C++ statement(s) on output stream "out".
876 
877 void TGeoBBox::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
878 {
879  if (TObject::TestBit(kGeoSavePrimitive)) return;
880  out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
881  out << " dx = " << fDX << ";" << std::endl;
882  out << " dy = " << fDY << ";" << std::endl;
883  out << " dz = " << fDZ << ";" << std::endl;
884  if (!TGeoShape::IsSameWithinTolerance(fOrigin[0],0) ||
885  !TGeoShape::IsSameWithinTolerance(fOrigin[1],0) ||
886  !TGeoShape::IsSameWithinTolerance(fOrigin[2],0)) {
887  out << " origin[0] = " << fOrigin[0] << ";" << std::endl;
888  out << " origin[1] = " << fOrigin[1] << ";" << std::endl;
889  out << " origin[2] = " << fOrigin[2] << ";" << std::endl;
890  out << " TGeoShape *" << GetPointerName() << " = new TGeoBBox(\"" << GetName() << "\", dx,dy,dz,origin);" << std::endl;
891  } else {
892  out << " TGeoShape *" << GetPointerName() << " = new TGeoBBox(\"" << GetName() << "\", dx,dy,dz);" << std::endl;
893  }
894  TObject::SetBit(TGeoShape::kGeoSavePrimitive);
895 }
896 
897 ////////////////////////////////////////////////////////////////////////////////
898 /// Set parameters of the box.
899 
900 void TGeoBBox::SetBoxDimensions(Double_t dx, Double_t dy, Double_t dz, Double_t *origin)
901 {
902  fDX = dx;
903  fDY = dy;
904  fDZ = dz;
905  if (origin) {
906  fOrigin[0] = origin[0];
907  fOrigin[1] = origin[1];
908  fOrigin[2] = origin[2];
909  }
910  if (TMath::Abs(fDX)<TGeoShape::Tolerance() &&
911  TMath::Abs(fDY)<TGeoShape::Tolerance() &&
912  TMath::Abs(fDZ)<TGeoShape::Tolerance()) return;
913  if ((fDX<0) || (fDY<0) || (fDZ<0)) SetShapeBit(kGeoRunTimeShape);
914 }
915 
916 ////////////////////////////////////////////////////////////////////////////////
917 /// Set dimensions based on the array of parameters
918 /// param[0] - half-length in x
919 /// param[1] - half-length in y
920 /// param[2] - half-length in z
921 
922 void TGeoBBox::SetDimensions(Double_t *param)
923 {
924  if (!param) {
925  Error("SetDimensions", "null parameters");
926  return;
927  }
928  fDX = param[0];
929  fDY = param[1];
930  fDZ = param[2];
931  if (TMath::Abs(fDX)<TGeoShape::Tolerance() &&
932  TMath::Abs(fDY)<TGeoShape::Tolerance() &&
933  TMath::Abs(fDZ)<TGeoShape::Tolerance()) return;
934  if ((fDX<0) || (fDY<0) || (fDZ<0)) SetShapeBit(kGeoRunTimeShape);
935 }
936 
937 ////////////////////////////////////////////////////////////////////////////////
938 /// Fill box vertices to an array.
939 
940 void TGeoBBox::SetBoxPoints(Double_t *points) const
941 {
942  TGeoBBox::SetPoints(points);
943 }
944 
945 ////////////////////////////////////////////////////////////////////////////////
946 /// Fill box points.
947 
948 void TGeoBBox::SetPoints(Double_t *points) const
949 {
950  if (!points) return;
951  Double_t xmin,xmax,ymin,ymax,zmin,zmax;
952  xmin = -fDX+fOrigin[0];
953  xmax = fDX+fOrigin[0];
954  ymin = -fDY+fOrigin[1];
955  ymax = fDY+fOrigin[1];
956  zmin = -fDZ+fOrigin[2];
957  zmax = fDZ+fOrigin[2];
958  points[ 0] = xmin; points[ 1] = ymin; points[ 2] = zmin;
959  points[ 3] = xmin; points[ 4] = ymax; points[ 5] = zmin;
960  points[ 6] = xmax; points[ 7] = ymax; points[ 8] = zmin;
961  points[ 9] = xmax; points[10] = ymin; points[11] = zmin;
962  points[12] = xmin; points[13] = ymin; points[14] = zmax;
963  points[15] = xmin; points[16] = ymax; points[17] = zmax;
964  points[18] = xmax; points[19] = ymax; points[20] = zmax;
965  points[21] = xmax; points[22] = ymin; points[23] = zmax;
966 }
967 
968 ////////////////////////////////////////////////////////////////////////////////
969 /// Fill box points.
970 
971 void TGeoBBox::SetPoints(Float_t *points) const
972 {
973  if (!points) return;
974  Double_t xmin,xmax,ymin,ymax,zmin,zmax;
975  xmin = -fDX+fOrigin[0];
976  xmax = fDX+fOrigin[0];
977  ymin = -fDY+fOrigin[1];
978  ymax = fDY+fOrigin[1];
979  zmin = -fDZ+fOrigin[2];
980  zmax = fDZ+fOrigin[2];
981  points[ 0] = xmin; points[ 1] = ymin; points[ 2] = zmin;
982  points[ 3] = xmin; points[ 4] = ymax; points[ 5] = zmin;
983  points[ 6] = xmax; points[ 7] = ymax; points[ 8] = zmin;
984  points[ 9] = xmax; points[10] = ymin; points[11] = zmin;
985  points[12] = xmin; points[13] = ymin; points[14] = zmax;
986  points[15] = xmin; points[16] = ymax; points[17] = zmax;
987  points[18] = xmax; points[19] = ymax; points[20] = zmax;
988  points[21] = xmax; points[22] = ymin; points[23] = zmax;
989 }
990 
991 ////////////////////////////////////////////////////////////////////////////////
992 ////// fill size of this 3-D object
993 //// TVirtualGeoPainter *painter = gGeoManager->GetGeomPainter();
994 //// if (painter) painter->AddSize3D(8, 12, 6);
995 
996 void TGeoBBox::Sizeof3D() const
997 {
998 }
999 
1000 ////////////////////////////////////////////////////////////////////////////////
1001 /// Fills a static 3D buffer and returns a reference.
1002 
1003 const TBuffer3D & TGeoBBox::GetBuffer3D(Int_t reqSections, Bool_t localFrame) const
1004 {
1005  static TBuffer3D buffer(TBuffer3DTypes::kGeneric);
1006 
1007  FillBuffer3D(buffer, reqSections, localFrame);
1008 
1009  // TODO: A box itself has has nothing more as already described
1010  // by bounding box. How will viewer interpret?
1011  if (reqSections & TBuffer3D::kRawSizes) {
1012  if (buffer.SetRawSizes(8, 3*8, 12, 3*12, 6, 6*6)) {
1013  buffer.SetSectionsValid(TBuffer3D::kRawSizes);
1014  }
1015  }
1016  if ((reqSections & TBuffer3D::kRaw) && buffer.SectionsValid(TBuffer3D::kRawSizes)) {
1017  SetPoints(buffer.fPnts);
1018  if (!buffer.fLocalFrame) {
1019  TransformPoints(buffer.fPnts, buffer.NbPnts());
1020  }
1021 
1022  SetSegsAndPols(buffer);
1023  buffer.SetSectionsValid(TBuffer3D::kRaw);
1024  }
1025 
1026  return buffer;
1027 }
1028 
1029 ////////////////////////////////////////////////////////////////////////////////
1030 /// Fills the supplied buffer, with sections in desired frame
1031 /// See TBuffer3D.h for explanation of sections, frame etc.
1032 
1033 void TGeoBBox::FillBuffer3D(TBuffer3D & buffer, Int_t reqSections, Bool_t localFrame) const
1034 {
1035  TGeoShape::FillBuffer3D(buffer, reqSections, localFrame);
1036 
1037  if (reqSections & TBuffer3D::kBoundingBox) {
1038  Double_t halfLengths[3] = { fDX, fDY, fDZ };
1039  buffer.SetAABoundingBox(fOrigin, halfLengths);
1040 
1041  if (!buffer.fLocalFrame) {
1042  TransformPoints(buffer.fBBVertex[0], 8);
1043  }
1044  buffer.SetSectionsValid(TBuffer3D::kBoundingBox);
1045  }
1046 }
1047 
1048 ////////////////////////////////////////////////////////////////////////////////
1049 /// Check the inside status for each of the points in the array.
1050 /// Input: Array of point coordinates + vector size
1051 /// Output: Array of Booleans for the inside of each point
1052 
1053 void TGeoBBox::Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
1054 {
1055  for (Int_t i=0; i<vecsize; i++) inside[i] = Contains(&points[3*i]);
1056 }
1057 
1058 ////////////////////////////////////////////////////////////////////////////////
1059 /// Compute the normal for an array o points so that norm.dot.dir is positive
1060 /// Input: Arrays of point coordinates and directions + vector size
1061 /// Output: Array of normal directions
1062 
1063 void TGeoBBox::ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
1064 {
1065  for (Int_t i=0; i<vecsize; i++) ComputeNormal(&points[3*i], &dirs[3*i], &norms[3*i]);
1066 }
1067 
1068 ////////////////////////////////////////////////////////////////////////////////
1069 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
1070 
1071 void TGeoBBox::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
1072 {
1073  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromInside(&points[3*i], &dirs[3*i], 3, step[i]);
1074 }
1075 
1076 ////////////////////////////////////////////////////////////////////////////////
1077 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
1078 
1079 void TGeoBBox::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
1080 {
1081  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromOutside(&points[3*i], &dirs[3*i], 3, step[i]);
1082 }
1083 
1084 ////////////////////////////////////////////////////////////////////////////////
1085 /// Compute safe distance from each of the points in the input array.
1086 /// Input: Array of point coordinates, array of statuses for these points, size of the arrays
1087 /// Output: Safety values
1088 
1089 void TGeoBBox::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
1090 {
1091  for (Int_t i=0; i<vecsize; i++) safe[i] = Safety(&points[3*i], inside[i]);
1092 }