ROOT   6.14/05 Reference Guide
TGeoArb8.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Andrei Gheata 31/01/02
3
4 /*************************************************************************
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. * 9 * For the list of contributors see$ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11
12 #include "TGeoArb8.h"
13
14 #include "Riostream.h"
15 #include "TBuffer.h"
16 #include "TGeoManager.h"
17 #include "TGeoVolume.h"
18 #include "TGeoMatrix.h"
19 #include "TMath.h"
20
22
23 /** \class TGeoArb8
24 \ingroup Geometry_classes
25
26 An arbitrary trapezoid with less than 8 vertices standing on
27 two parallel planes perpendicular to Z axis. Parameters :
28  - dz - half length in Z;
29  - xy[8][2] - vector of (x,y) coordinates of vertices
30  - first four points (xy[i][j], i<4, j<2) are the (x,y)
31  coordinates of the vertices sitting on the -dz plane;
32  - last four points (xy[i][j], i>=4, j<2) are the (x,y)
33  coordinates of the vertices sitting on the +dz plane;
34
35 The order of defining the vertices of an arb8 is the following :
36  - point 0 is connected with points 1,3,4
37  - point 1 is connected with points 0,2,5
38  - point 2 is connected with points 1,3,6
39  - point 3 is connected with points 0,2,7
40  - point 4 is connected with points 0,5,7
41  - point 5 is connected with points 1,4,6
42  - point 6 is connected with points 2,5,7
43  - point 7 is connected with points 3,4,6
44
45 Points can be identical in order to create shapes with less than
46 8 vertices.
47
48 Begin_Macro(source)
49 {
50  TCanvas *c = new TCanvas("c", "c",0,0,600,600);
51  new TGeoManager("arb8", "poza12");
52  TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
53  TGeoMedium *med = new TGeoMedium("MED",1,mat);
54  TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
55  gGeoManager->SetTopVolume(top);
56  TGeoArb8 *arb = new TGeoArb8(20);
57  arb->SetVertex(0,-30,-25);
58  arb->SetVertex(1,-25,25);
59  arb->SetVertex(2,5,25);
60  arb->SetVertex(3,25,-25);
61  arb->SetVertex(4,-28,-23);
62  arb->SetVertex(5,-23,27);
63  arb->SetVertex(6,-23,27);
64  arb->SetVertex(7,13,-27);
65  TGeoVolume *vol = new TGeoVolume("ARB8",arb,med);
66  vol->SetLineWidth(2);
68  gGeoManager->CloseGeometry();
69  gGeoManager->SetNsegments(80);
70  top->Draw();
72  view->ShowAxis();
73 }
74 End_Macro
75 */
76
77 /** \class TGeoGtra
78 \ingroup Geometry_classes
79
80 Gtra is a twisted trapezoid.
81 i.e. one for which the faces perpendicular
82 to z are trapezia and their centres are not the same x, y. It has 12
83 parameters: the half length in z, the polar angles from the centre of
84 the face at low z to that at high z, twist, H1 the half length in y at low z,
85 LB1 the half length in x at low z and y low edge, LB2 the half length
86 in x at low z and y high edge, TH1 the angle w.r.t. the y axis from the
87 centre of low y edge to the centre of the high y edge, and H2, LB2,
88 LH2, TH2, the corresponding quantities at high z.
89
90 Begin_Macro(source)
91 {
92  TCanvas *c = new TCanvas("c", "c",0,0,600,600);
93  new TGeoManager("gtra", "poza11");
94  TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
95  TGeoMedium *med = new TGeoMedium("MED",1,mat);
96  TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
97  gGeoManager->SetTopVolume(top);
98  TGeoVolume *vol = gGeoManager->MakeGtra("Gtra",med, 30,15,30,30,20,10,15,0,20,10,15,0);
99  vol->SetLineWidth(2);
101  gGeoManager->CloseGeometry();
102  gGeoManager->SetNsegments(80);
103  top->Draw();
105  view->ShowAxis();
106 }
107 End_Macro
108 */
109
110 /** \class TGeoTrap
111 \ingroup Geometry_classes
112
113 TRAP is a general trapezoid, i.e. one for which the faces perpendicular
114 to z are trapezia and their centres are not the same x, y. It has 11
115 parameters: the half length in z, the polar angles from the centre of
116 the face at low z to that at high z, H1 the half length in y at low z,
117 LB1 the half length in x at low z and y low edge, LB2 the half length
118 in x at low z and y high edge, TH1 the angle w.r.t. the y axis from the
119 centre of low y edge to the centre of the high y edge, and H2, LB2,
120 LH2, TH2, the corresponding quantities at high z.
121
122 Begin_Macro(source)
123 {
124  TCanvas *c = new TCanvas("c", "c",0,0,600,600);
125  new TGeoManager("trap", "poza10");
126  TGeoMaterial *mat = new TGeoMaterial("Al", 26.98,13,2.7);
127  TGeoMedium *med = new TGeoMedium("MED",1,mat);
128  TGeoVolume *top = gGeoManager->MakeBox("TOP",med,100,100,100);
129  gGeoManager->SetTopVolume(top);
130  TGeoVolume *vol = gGeoManager->MakeTrap("Trap",med, 30,15,30,20,10,15,0,20,10,15,0);
131  vol->SetLineWidth(2);
133  gGeoManager->CloseGeometry();
134  gGeoManager->SetNsegments(80);
135  top->Draw();
137  view->ShowAxis();
138 }
139 End_Macro
140 */
141
142 ////////////////////////////////////////////////////////////////////////////////
143 /// Default constructor.
144
146 {
147  fDz = 0;
148  fTwist = 0;
149  for (Int_t i=0; i<8; i++) {
150  fXY[i][0] = 0.0;
151  fXY[i][1] = 0.0;
152  }
154 }
155
156 ////////////////////////////////////////////////////////////////////////////////
157 /// Constructor. If the array of vertices is not null, this should be
158 /// in the format : (x0, y0, x1, y1, ... , x7, y7)
159
161  :TGeoBBox(0,0,0)
162 {
163  fDz = dz;
164  fTwist = 0;
166  if (vertices) {
167  for (Int_t i=0; i<8; i++) {
168  fXY[i][0] = vertices[2*i];
169  fXY[i][1] = vertices[2*i+1];
170  }
171  ComputeTwist();
172  ComputeBBox();
173  } else {
174  for (Int_t i=0; i<8; i++) {
175  fXY[i][0] = 0.0;
176  fXY[i][1] = 0.0;
177  }
178  }
179 }
180
181 ////////////////////////////////////////////////////////////////////////////////
182 /// Named constructor. If the array of vertices is not null, this should be
183 /// in the format : (x0, y0, x1, y1, ... , x7, y7)
184
185 TGeoArb8::TGeoArb8(const char *name, Double_t dz, Double_t *vertices)
186  :TGeoBBox(name, 0,0,0)
187 {
188  fDz = dz;
189  fTwist = 0;
191  if (vertices) {
192  for (Int_t i=0; i<8; i++) {
193  fXY[i][0] = vertices[2*i];
194  fXY[i][1] = vertices[2*i+1];
195  }
196  ComputeTwist();
197  ComputeBBox();
198  } else {
199  for (Int_t i=0; i<8; i++) {
200  fXY[i][0] = 0.0;
201  fXY[i][1] = 0.0;
202  }
203  }
204 }
205
206 ////////////////////////////////////////////////////////////////////////////////
207 /// Copy constructor.
208
210  TGeoBBox(ga8),
211  fDz(ga8.fDz),
212  fTwist(ga8.fTwist)
213 {
214  for(Int_t i=0; i<8; i++) {
215  fXY[i][0]=ga8.fXY[i][0];
216  fXY[i][1]=ga8.fXY[i][1];
217  }
218 }
219
220 ////////////////////////////////////////////////////////////////////////////////
221 /// Assignment operator.
222
224 {
225  if(this!=&ga8) {
226  TGeoBBox::operator=(ga8);
227  fDz=ga8.fDz;
228  fTwist=ga8.fTwist;
229  for(Int_t i=0; i<8; i++) {
230  fXY[i][0]=ga8.fXY[i][0];
231  fXY[i][1]=ga8.fXY[i][1];
232  }
233  }
234  return *this;
235 }
236
237 ////////////////////////////////////////////////////////////////////////////////
238 /// Destructor.
239
241 {
242  if (fTwist) delete [] fTwist;
243 }
244
245 ////////////////////////////////////////////////////////////////////////////////
246 /// Computes capacity of the shape in [length^3].
247
249 {
250  Int_t i,j;
251  Double_t capacity = 0;
252  for (i=0; i<4; i++) {
253  j = (i+1)%4;
254  capacity += 0.25*fDz*((fXY[i][0]+fXY[i+4][0])*(fXY[j][1]+fXY[j+4][1]) -
255  (fXY[j][0]+fXY[j+4][0])*(fXY[i][1]+fXY[i+4][1]) +
256  (1./3)*((fXY[i+4][0]-fXY[i][0])*(fXY[j+4][1]-fXY[j][1]) -
257  (fXY[j][0]-fXY[j+4][0])*(fXY[i][1]-fXY[i+4][1])));
258  }
259  return TMath::Abs(capacity);
260 }
261
262 ////////////////////////////////////////////////////////////////////////////////
263 /// Computes bounding box for an Arb8 shape.
264
266 {
268  xmin = xmax = fXY[0][0];
269  ymin = ymax = fXY[0][1];
270
271  for (Int_t i=1; i<8; i++) {
272  if (xmin>fXY[i][0]) xmin=fXY[i][0];
273  if (xmax<fXY[i][0]) xmax=fXY[i][0];
274  if (ymin>fXY[i][1]) ymin=fXY[i][1];
275  if (ymax<fXY[i][1]) ymax=fXY[i][1];
276  }
277  fDX = 0.5*(xmax-xmin);
278  fDY = 0.5*(ymax-ymin);
279  fDZ = fDz;
280  fOrigin[0] = 0.5*(xmax+xmin);
281  fOrigin[1] = 0.5*(ymax+ymin);
282  fOrigin[2] = 0;
284 }
285
286 ////////////////////////////////////////////////////////////////////////////////
287 /// Computes tangents of twist angles (angles between projections on XY plane
288 /// of corresponding -dz +dz edges). Computes also if the vertices are defined
289 /// clockwise or anti-clockwise.
290
292 {
293  Double_t twist[4];
294  Bool_t twisted = kFALSE;
295  Double_t dx1, dy1, dx2, dy2;
296  Bool_t singleBottom = kTRUE;
297  Bool_t singleTop = kTRUE;
298  Int_t i;
299  for (i=0; i<4; i++) {
300  dx1 = fXY[(i+1)%4][0]-fXY[i][0];
301  dy1 = fXY[(i+1)%4][1]-fXY[i][1];
303  twist[i] = 0;
304  continue;
305  }
306  singleBottom = kFALSE;
307  dx2 = fXY[4+(i+1)%4][0]-fXY[4+i][0];
308  dy2 = fXY[4+(i+1)%4][1]-fXY[4+i][1];
310  twist[i] = 0;
311  continue;
312  }
313  singleTop = kFALSE;
314  twist[i] = dy1*dx2 - dx1*dy2;
315  if (TMath::Abs(twist[i])<TGeoShape::Tolerance()) {
316  twist[i] = 0;
317  continue;
318  }
319  twist[i] = TMath::Sign(1.,twist[i]);
320  twisted = kTRUE;
321  }
322  if (twisted) {
323  if (fTwist) delete [] fTwist;
324  fTwist = new Double_t[4];
325  memcpy(fTwist, &twist[0], 4*sizeof(Double_t));
326  }
327  if (singleBottom) {
328  for (i=0; i<4; i++) {
329  fXY[i][0] += 1.E-8*fXY[i+4][0];
330  fXY[i][1] += 1.E-8*fXY[i+4][1];
331  }
332  }
333  if (singleTop) {
334  for (i=0; i<4; i++) {
335  fXY[i+4][0] += 1.E-8*fXY[i][0];
336  fXY[i+4][1] += 1.E-8*fXY[i][1];
337  }
338  }
339  Double_t sum1 = 0.;
340  Double_t sum2 = 0.;
341  Int_t j;
342  for (i=0; i<4; i++) {
343  j = (i+1)%4;
344  sum1 += fXY[i][0]*fXY[j][1]-fXY[j][0]*fXY[i][1];
345  sum2 += fXY[i+4][0]*fXY[j+4][1]-fXY[j+4][0]*fXY[i+4][1];
346  }
347  if (sum1*sum2 < -TGeoShape::Tolerance()) {
348  Fatal("ComputeTwist", "Shape %s type Arb8: Lower/upper faces defined with opposite clockwise", GetName());
349  return;
350  }
351  if (sum1>TGeoShape::Tolerance()) {
352  Error("ComputeTwist", "Shape %s type Arb8: Vertices must be defined clockwise in XY planes. Re-ordering...", GetName());
353  Double_t xtemp, ytemp;
354  xtemp = fXY[1][0];
355  ytemp = fXY[1][1];
356  fXY[1][0] = fXY[3][0];
357  fXY[1][1] = fXY[3][1];
358  fXY[3][0] = xtemp;
359  fXY[3][1] = ytemp;
360  xtemp = fXY[5][0];
361  ytemp = fXY[5][1];
362  fXY[5][0] = fXY[7][0];
363  fXY[5][1] = fXY[7][1];
364  fXY[7][0] = xtemp;
365  fXY[7][1] = ytemp;
366  }
367  // Check for illegal crossings.
368  Bool_t illegal_cross = kFALSE;
369  illegal_cross = TGeoShape::IsSegCrossing(fXY[0][0],fXY[0][1],fXY[1][0],fXY[1][1],
370  fXY[2][0],fXY[2][1],fXY[3][0],fXY[3][1]);
371  if (!illegal_cross)
372  illegal_cross = TGeoShape::IsSegCrossing(fXY[4][0],fXY[4][1],fXY[5][0],fXY[5][1],
373  fXY[6][0],fXY[6][1],fXY[7][0],fXY[7][1]);
374  if (illegal_cross) {
375  Error("ComputeTwist", "Shape %s type Arb8: Malformed polygon with crossing opposite segments", GetName());
376  InspectShape();
377  }
378 }
379
380 ////////////////////////////////////////////////////////////////////////////////
381 /// Get twist for segment I in range [0,3]
382
384 {
385  if (!fTwist) return 0.;
386  if (iseg<0 || iseg>3) return 0.;
387  return fTwist[iseg];
388 }
389
390 ////////////////////////////////////////////////////////////////////////////////
391 /// Get index of the edge of the quadrilater represented by vert closest to point.
392 /// If [P1,P2] is the closest segment and P is the point, the function returns the fraction of the
393 /// projection of (P1P) over (P1P2). If projection of P is not in range [P1,P2] return -1.
394
395 Double_t TGeoArb8::GetClosestEdge(const Double_t *point, Double_t *vert, Int_t &isegment) const
396 {
397  isegment = 0;
398  Int_t isegmin = 0;
399  Int_t i1, i2;
400  Double_t p1[2], p2[2];
401  Double_t lsq, ssq, dx, dy, dpx, dpy, u;
402  Double_t umin = -1.;
403  Double_t safe=1E30;
404  for (i1=0; i1<4; i1++) {
405  if (TGeoShape::IsSameWithinTolerance(safe,0)) {
406  isegment = isegmin;
407  return umin;
408  }
409  i2 = (i1+1)%4;
410  p1[0] = vert[2*i1];
411  p1[1] = vert[2*i1+1];
412  p2[0] = vert[2*i2];
413  p2[1] = vert[2*i2+1];
414  dx = p2[0] - p1[0];
415  dy = p2[1] - p1[1];
416  dpx = point[0] - p1[0];
417  dpy = point[1] - p1[1];
418  lsq = dx*dx + dy*dy;
419  // Current segment collapsed to a point
421  ssq = dpx*dpx + dpy*dpy;
422  if (ssq < safe) {
423  safe = ssq;
424  isegmin = i1;
425  umin = -1;
426  }
427  continue;
428  }
429  // Projection fraction
430  u = (dpx*dx + dpy*dy)/lsq;
431  // If projection of P is outside P1P2 limits, take the distance from P to the closest of P1 and P2
432  if (u>1) {
433  // Outside, closer to P2
434  dpx = point[0]-p2[0];
435  dpy = point[1]-p2[1];
436  u = -1.;
437  } else {
438  if (u>=0) {
439  // Projection inside
440  dpx -= u*dx;
441  dpy -= u*dy;
442  } else {
443  // Outside, closer to P1
444  u = -1.;
445  }
446  }
447  ssq = dpx*dpx + dpy*dpy;
448  if (ssq < safe) {
449  safe = ssq;
450  isegmin = i1;
451  umin = u;
452  }
453  }
454  isegment = isegmin;
455  // safe = TMath::Sqrt(safe);
456  return umin;
457 }
458
459 ////////////////////////////////////////////////////////////////////////////////
460 /// Compute normal to closest surface from POINT.
461
462 void TGeoArb8::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
463 {
464  Double_t safc;
465  Double_t x0, y0, z0, x1, y1, z1, x2, y2, z2;
466  Double_t ax, ay, az, bx, by, bz;
467  Double_t fn;
468  safc = fDz-TMath::Abs(point[2]);
469  if (safc<10.*TGeoShape::Tolerance()) {
470  memset(norm,0,3*sizeof(Double_t));
471  norm[2] = (dir[2]>0)?1:(-1);
472  return;
473  }
474  Double_t vert[8];
475  SetPlaneVertices(point[2], vert);
476  // Get the closest edge (point should be on this edge within tolerance)
477  Int_t iseg;
478  Double_t frac = GetClosestEdge(point, vert, iseg);
479  if (frac<0) frac = 0.;
480  Int_t jseg = (iseg+1)%4;
481  x0 = vert[2*iseg];
482  y0 = vert[2*iseg+1];
483  z0 = point[2];
484  x2 = vert[2*jseg];
485  y2 = vert[2*jseg+1];
486  z2 = point[2];
487  x0 += frac*(x2-x0);
488  y0 += frac*(y2-y0);
489  x1 = fXY[iseg+4][0];
490  y1 = fXY[iseg+4][1];
491  z1 = fDz;
492  x1 += frac*(fXY[jseg+4][0]-x1);
493  y1 += frac*(fXY[jseg+4][1]-y1);
494  ax = x1-x0;
495  ay = y1-y0;
496  az = z1-z0;
497  bx = x2-x0;
498  by = y2-y0;
499  bz = z2-z0;
500  // Cross product of the vector given by the section segment (that contains the point) at z=point[2]
501  // and the vector connecting the point projection to its correspondent on the top edge (hard to swallow, isn't it?)
502  norm[0] = ay*bz-az*by;
503  norm[1] = az*bx-ax*bz;
504  norm[2] = ax*by-ay*bx;
505  fn = TMath::Sqrt(norm[0]*norm[0]+norm[1]*norm[1]+norm[2]*norm[2]);
506  // If point is on one edge, fn may be very small and the normal does not make sense - avoid divzero
507  if (fn<1E-10) {
508  norm[0] = 1.;
509  norm[1] = 0.;
510  norm[2] = 0.;
511  } else {
512  norm[0] /= fn;
513  norm[1] /= fn;
514  norm[2] /= fn;
515  }
516  // Make sure the dot product of the normal and the direction is positive.
517  if (dir[0]>-2. && dir[0]*norm[0]+dir[1]*norm[1]+dir[2]*norm[2] < 0) {
518  norm[0] = -norm[0];
519  norm[1] = -norm[1];
520  norm[2] = -norm[2];
521  }
522 }
523
524 ////////////////////////////////////////////////////////////////////////////////
525 /// Test if point is inside this shape.
526
527 Bool_t TGeoArb8::Contains(const Double_t *point) const
528 {
529  // first check Z range
530  if (TMath::Abs(point[2]) > fDz) return kFALSE;
531  // compute intersection between Z plane containing point and the arb8
532  Double_t poly[8];
533  // memset(&poly[0], 0, 8*sizeof(Double_t));
534  // SetPlaneVertices(point[2], &poly[0]);
535  Double_t cf = 0.5*(fDz-point[2])/fDz;
536  Int_t i;
537  for (i=0; i<4; i++) {
538  poly[2*i] = fXY[i+4][0]+cf*(fXY[i][0]-fXY[i+4][0]);
539  poly[2*i+1] = fXY[i+4][1]+cf*(fXY[i][1]-fXY[i+4][1]);
540  }
541  return InsidePolygon(point[0],point[1],poly);
542 }
543
544 ////////////////////////////////////////////////////////////////////////////////
545 /// Computes distance to plane ipl :
546 /// - ipl=0 : points 0,4,1,5
547 /// - ipl=1 : points 1,5,2,6
548 /// - ipl=2 : points 2,6,3,7
549 /// - ipl=3 : points 3,7,0,4
550
551 Double_t TGeoArb8::DistToPlane(const Double_t *point, const Double_t *dir, Int_t ipl, Bool_t in) const
552 {
553  Double_t xa,xb,xc,xd;
554  Double_t ya,yb,yc,yd;
555  Double_t eps = 10.*TGeoShape::Tolerance();
556  Double_t norm[3];
557  Double_t dirp[3];
558  Double_t ndotd = 0;
559  Int_t j = (ipl+1)%4;
560  xa=fXY[ipl][0];
561  ya=fXY[ipl][1];
562  xb=fXY[ipl+4][0];
563  yb=fXY[ipl+4][1];
564  xc=fXY[j][0];
565  yc=fXY[j][1];
566  xd=fXY[4+j][0];
567  yd=fXY[4+j][1];
568  Double_t dz2 =0.5/fDz;
569  Double_t tx1 =dz2*(xb-xa);
570  Double_t ty1 =dz2*(yb-ya);
571  Double_t tx2 =dz2*(xd-xc);
572  Double_t ty2 =dz2*(yd-yc);
573  Double_t dzp =fDz+point[2];
574  Double_t xs1 =xa+tx1*dzp;
575  Double_t ys1 =ya+ty1*dzp;
576  Double_t xs2 =xc+tx2*dzp;
577  Double_t ys2 =yc+ty2*dzp;
578  Double_t dxs =xs2-xs1;
579  Double_t dys =ys2-ys1;
580  Double_t dtx =tx2-tx1;
581  Double_t dty =ty2-ty1;
582  Double_t a=(dtx*dir[1]-dty*dir[0]+(tx1*ty2-tx2*ty1)*dir[2])*dir[2];
583  Double_t signa = TMath::Sign(1.,a);
584  Double_t b=dxs*dir[1]-dys*dir[0]+(dtx*point[1]-dty*point[0]+ty2*xs1-ty1*xs2
585  +tx1*ys2-tx2*ys1)*dir[2];
586  Double_t c=dxs*point[1]-dys*point[0]+xs1*ys2-xs2*ys1;
588  Double_t x1,x2,y1,y2,xp,yp,zi;
589  if (TMath::Abs(a)<eps) {
590  // Surface is planar
591  if (TMath::Abs(b)<eps) return TGeoShape::Big(); // Track parallel to surface
592  s=-c/b;
593  if (TMath::Abs(s)<1.E-6 && TMath::Abs(TMath::Abs(point[2])-fDz)>eps) {
594  memcpy(dirp,dir,3*sizeof(Double_t));
595  dirp[0] = -3;
596  // Compute normal pointing outside
597  ((TGeoArb8*)this)->ComputeNormal(point,dirp,norm);
598  ndotd = dir[0]*norm[0]+dir[1]*norm[1]+dir[2]*norm[2];
599  if (!in) ndotd*=-1.;
600  if (ndotd>0) {
601  s = TMath::Max(0.,s);
602  zi = (point[0]-xs1)*(point[0]-xs2)+(point[1]-ys1)*(point[1]-ys2);
603  if (zi<=0) return s;
604  }
605  return TGeoShape::Big();
606  }
607  if (s<0) return TGeoShape::Big();
608  } else {
609  Double_t d=b*b-4*a*c;
610  if (d<0) return TGeoShape::Big();
611  Double_t smin=0.5*(-b-signa*TMath::Sqrt(d))/a;
612  Double_t smax=0.5*(-b+signa*TMath::Sqrt(d))/a;
613  s = smin;
614  if (TMath::Abs(s)<1.E-6 && TMath::Abs(TMath::Abs(point[2])-fDz)>eps) {
615  memcpy(dirp,dir,3*sizeof(Double_t));
616  dirp[0] = -3;
617  // Compute normal pointing outside
618  ((TGeoArb8*)this)->ComputeNormal(point,dirp,norm);
619  ndotd = dir[0]*norm[0]+dir[1]*norm[1]+dir[2]*norm[2];
620  if (!in) ndotd*=-1.;
621  if (ndotd>0) return TMath::Max(0.,s);
622  s = 0.; // ignore
623  }
624  if (s>eps) {
625  // Check smin
626  zi=point[2]+s*dir[2];
627  if (TMath::Abs(zi)<fDz) {
628  x1=xs1+tx1*dir[2]*s;
629  x2=xs2+tx2*dir[2]*s;
630  xp=point[0]+s*dir[0];
631  y1=ys1+ty1*dir[2]*s;
632  y2=ys2+ty2*dir[2]*s;
633  yp=point[1]+s*dir[1];
634  zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
635  if (zi<=0) return s;
636  }
637  }
638  // Smin failed, try smax
639  s=smax;
640  if (TMath::Abs(s)<1.E-6 && TMath::Abs(TMath::Abs(point[2])-fDz)>eps) {
641  memcpy(dirp,dir,3*sizeof(Double_t));
642  dirp[0] = -3;
643  // Compute normal pointing outside
644  ((TGeoArb8*)this)->ComputeNormal(point,dirp,norm);
645  ndotd = dir[0]*norm[0]+dir[1]*norm[1]+dir[2]*norm[2];
646  if (!in) ndotd*=-1.;
647  if (ndotd>0) s = TMath::Max(0.,s);
648  else s = TGeoShape::Big();
649  return s;
650  }
651  }
652  if (s>eps) {
653  // Check smin
654  zi=point[2]+s*dir[2];
655  if (TMath::Abs(zi)<fDz) {
656  x1=xs1+tx1*dir[2]*s;
657  x2=xs2+tx2*dir[2]*s;
658  xp=point[0]+s*dir[0];
659  y1=ys1+ty1*dir[2]*s;
660  y2=ys2+ty2*dir[2]*s;
661  yp=point[1]+s*dir[1];
662  zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
663  if (zi<=0) return s;
664  }
665  }
666  return TGeoShape::Big();
667 }
668
669 ////////////////////////////////////////////////////////////////////////////////
670 /// Computes distance from outside point to surface of the shape.
671
672 Double_t TGeoArb8::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t /*iact*/, Double_t step, Double_t * /*safe*/) const
673 {
674  Double_t sdist = TGeoBBox::DistFromOutside(point,dir, fDX, fDY, fDZ, fOrigin, step);
675  if (sdist>=step) return TGeoShape::Big();
676  Double_t snext;
677  // check Z planes
678  if (TMath::Abs(point[2])>fDz-1.E-8) {
679  Double_t pt[3];
680  if (point[2]*dir[2]<0) {
681  pt[2]=fDz*TMath::Sign(1.,point[2]);
682  snext = TMath::Max((pt[2]-point[2])/dir[2],0.);
683  for (Int_t j=0; j<2; j++) pt[j]=point[j]+snext*dir[j];
684  if (Contains(&pt[0])) return snext;
685  }
686  }
687  // check lateral faces
688  Double_t dist;
689  snext = TGeoShape::Big();
690  for (Int_t i=0; i<4; i++) {
691  dist = DistToPlane(point, dir, i, kFALSE);
692  if (dist<snext) snext = dist;
693  }
694  return snext;
695 }
696
697 ////////////////////////////////////////////////////////////////////////////////
698 /// Compute distance from inside point to surface of the shape.
699
700 Double_t TGeoArb8::DistFromInside(const Double_t *point, const Double_t *dir, Int_t /*iact*/, Double_t /*step*/, Double_t * /*safe*/) const
701 {
702  Int_t i;
703  Double_t distz = TGeoShape::Big();
704  Double_t distl = TGeoShape::Big();
705  Double_t dist;
706  Double_t pt[3] = {0.,0.,0.};
707  if (dir[2]<0) {
708  distz=(-fDz-point[2])/dir[2];
709  pt[2] = -fDz;
710  } else {
711  if (dir[2]>0) distz=(fDz-point[2])/dir[2];
712  pt[2] = fDz;
713  }
714  for (i=0; i<4; i++) {
715  dist=DistToPlane(point, dir, i, kTRUE);
716  if (dist<distl) distl = dist;
717  }
718  if (distz<distl) {
719  pt[0] = point[0]+distz*dir[0];
720  pt[1] = point[1]+distz*dir[1];
721  if (!Contains(pt)) distz = TGeoShape::Big();
722  }
723  dist = TMath::Min(distz, distl);
724  if (dist<0 || dist>1.E10) return 0.;
725  return dist;
726 #ifdef OLDALGORITHM
727 //#else
728 // compute distance to plane ipl :
729 // ipl=0 : points 0,4,1,5
730 // ipl=1 : points 1,5,2,6
731 // ipl=2 : points 2,6,3,7
732 // ipl=3 : points 3,7,0,4
733  Double_t distmin;
734  Bool_t lateral_cross = kFALSE;
735  if (dir[2]<0) {
736  distmin=(-fDz-point[2])/dir[2];
737  } else {
738  if (dir[2]>0) distmin =(fDz-point[2])/dir[2];
739  else distmin = TGeoShape::Big();
740  }
741  Double_t dz2 =0.5/fDz;
742  Double_t xa,xb,xc,xd;
743  Double_t ya,yb,yc,yd;
744  Double_t eps = 100.*TGeoShape::Tolerance();
745  for (Int_t ipl=0;ipl<4;ipl++) {
746  Int_t j = (ipl+1)%4;
747  xa=fXY[ipl][0];
748  ya=fXY[ipl][1];
749  xb=fXY[ipl+4][0];
750  yb=fXY[ipl+4][1];
751  xc=fXY[j][0];
752  yc=fXY[j][1];
753  xd=fXY[4+j][0];
754  yd=fXY[4+j][1];
755
756  Double_t tx1 =dz2*(xb-xa);
757  Double_t ty1 =dz2*(yb-ya);
758  Double_t tx2 =dz2*(xd-xc);
759  Double_t ty2 =dz2*(yd-yc);
760  Double_t dzp =fDz+point[2];
761  Double_t xs1 =xa+tx1*dzp;
762  Double_t ys1 =ya+ty1*dzp;
763  Double_t xs2 =xc+tx2*dzp;
764  Double_t ys2 =yc+ty2*dzp;
765  Double_t dxs =xs2-xs1;
766  Double_t dys =ys2-ys1;
767  Double_t dtx =tx2-tx1;
768  Double_t dty =ty2-ty1;
769  Double_t a=(dtx*dir[1]-dty*dir[0]+(tx1*ty2-tx2*ty1)*dir[2])*dir[2];
770  Double_t b=dxs*dir[1]-dys*dir[0]+(dtx*point[1]-dty*point[0]+ty2*xs1-ty1*xs2
771  +tx1*ys2-tx2*ys1)*dir[2];
772  Double_t c=dxs*point[1]-dys*point[0]+xs1*ys2-xs2*ys1;
774  if (TMath::Abs(a)<eps) {
775  if (TMath::Abs(b)<eps) continue;
776  s=-c/b;
777  if (s>eps && s < distmin) {
778  distmin =s;
779  lateral_cross=kTRUE;
780  }
781  continue;
782  }
783  Double_t d=b*b-4*a*c;
784  if (d>=0.) {
785  if (a>0) s=0.5*(-b-TMath::Sqrt(d))/a;
786  else s=0.5*(-b+TMath::Sqrt(d))/a;
787  if (s>eps) {
788  if (s < distmin) {
789  distmin = s;
790  lateral_cross = kTRUE;
791  }
792  } else {
793  if (a>0) s=0.5*(-b+TMath::Sqrt(d))/a;
794  else s=0.5*(-b-TMath::Sqrt(d))/a;
795  if (s>eps && s < distmin) {
796  distmin =s;
797  lateral_cross = kTRUE;
798  }
799  }
800  }
801  }
802
803  if (!lateral_cross) {
804  // We have to make sure that track crosses the top or bottom.
805  if (distmin > 1.E10) return TGeoShape::Tolerance();
806  Double_t pt[2];
807  pt[0] = point[0]+distmin*dir[0];
808  pt[1] = point[1]+distmin*dir[1];
809  // Check if propagated point is in the polygon
810  Double_t poly[8];
811  Int_t i = 0;
812  if (dir[2]>0.) i=4;
813  for (Int_t j=0; j<4; j++) {
814  poly[2*j] = fXY[j+i][0];
815  poly[2*j+1] = fXY[j+i][1];
816  }
817  if (!InsidePolygon(pt[0],pt[1],poly)) return TGeoShape::Tolerance();
818  }
819  return distmin;
820 #endif
821 }
822
823 ////////////////////////////////////////////////////////////////////////////////
824 /// Divide this shape along one axis.
825
826 TGeoVolume *TGeoArb8::Divide(TGeoVolume *voldiv, const char * /*divname*/, Int_t /*iaxis*/, Int_t /*ndiv*/,
827  Double_t /*start*/, Double_t /*step*/)
828 {
829  Error("Divide", "Division of an arbitrary trapezoid not implemented");
830  return voldiv;
831 }
832
833 ////////////////////////////////////////////////////////////////////////////////
834 /// Get shape range on a given axis.
835
837 {
838  xlo = 0;
839  xhi = 0;
840  Double_t dx = 0;
841  if (iaxis==3) {
842  xlo = -fDz;
843  xhi = fDz;
844  dx = xhi-xlo;
845  return dx;
846  }
847  return dx;
848 }
849
850 ////////////////////////////////////////////////////////////////////////////////
851 /// Fill vector param[4] with the bounding cylinder parameters. The order
852 /// is the following : Rmin, Rmax, Phi1, Phi2
853
855 {
856  // first compute rmin/rmax
857  Double_t rmaxsq = 0;
858  Double_t rsq;
859  Int_t i;
860  for (i=0; i<8; i++) {
861  rsq = fXY[i][0]*fXY[i][0] + fXY[i][1]*fXY[i][1];
862  rmaxsq = TMath::Max(rsq, rmaxsq);
863  }
864  param[0] = 0.; // Rmin
865  param[1] = rmaxsq; // Rmax
866  param[2] = 0.; // Phi1
867  param[3] = 360.; // Phi2
868 }
869
870 ////////////////////////////////////////////////////////////////////////////////
871 /// Fills real parameters of a positioned box inside this arb8. Returns 0 if successful.
872
873 Int_t TGeoArb8::GetFittingBox(const TGeoBBox *parambox, TGeoMatrix *mat, Double_t &dx, Double_t &dy, Double_t &dz) const
874 {
875  dx=dy=dz=0;
876  if (mat->IsRotation()) {
877  Error("GetFittingBox", "cannot handle parametrized rotated volumes");
878  return 1; // ### rotation not accepted ###
879  }
880  //--> translate the origin of the parametrized box to the frame of this box.
881  Double_t origin[3];
882  mat->LocalToMaster(parambox->GetOrigin(), origin);
883  if (!Contains(origin)) {
884  Error("GetFittingBox", "wrong matrix - parametrized box is outside this");
885  return 1; // ### wrong matrix ###
886  }
887  //--> now we have to get the valid range for all parametrized axis
888  Double_t dd[3];
889  dd[0] = parambox->GetDX();
890  dd[1] = parambox->GetDY();
891  dd[2] = parambox->GetDZ();
892  //-> check if Z range is fixed
893  if (dd[2]<0) {
894  dd[2] = TMath::Min(origin[2]+fDz, fDz-origin[2]);
895  if (dd[2]<0) {
896  Error("GetFittingBox", "wrong matrix");
897  return 1;
898  }
899  }
900  if (dd[0]>=0 && dd[1]>=0) {
901  dx = dd[0];
902  dy = dd[1];
903  dz = dd[2];
904  return 0;
905  }
906  //-> check now vertices at Z = origin[2] +/- dd[2]
907  Double_t upper[8];
908  Double_t lower[8];
909  SetPlaneVertices(origin[2]-dd[2], lower);
910  SetPlaneVertices(origin[2]+dd[2], upper);
911  Double_t ddmin=TGeoShape::Big();
912  for (Int_t iaxis=0; iaxis<2; iaxis++) {
913  if (dd[iaxis]>=0) continue;
914  ddmin=TGeoShape::Big();
915  for (Int_t ivert=0; ivert<4; ivert++) {
916  ddmin = TMath::Min(ddmin, TMath::Abs(origin[iaxis]-lower[2*ivert+iaxis]));
917  ddmin = TMath::Min(ddmin, TMath::Abs(origin[iaxis]-upper[2*ivert+iaxis]));
918  }
919  dd[iaxis] = ddmin;
920  }
921  dx = dd[0];
922  dy = dd[1];
923  dz = dd[2];
924  return 0;
925 }
926
927 ////////////////////////////////////////////////////////////////////////////////
928 /// Computes normal to plane defined by P1, P2 and P3
929
931 {
932  Double_t cross = 0.;
933  Double_t v1[3], v2[3];
934  Int_t i;
935  for (i=0; i<3; i++) {
936  v1[i] = p2[i] - p1[i];
937  v2[i] = p3[i] - p1[i];
938  }
939  norm[0] = v1[1]*v2[2]-v1[2]*v2[1];
940  cross += norm[0]*norm[0];
941  norm[1] = v1[2]*v2[0]-v1[0]*v2[2];
942  cross += norm[1]*norm[1];
943  norm[2] = v1[0]*v2[1]-v1[1]*v2[0];
944  cross += norm[2]*norm[2];
945  if (TMath::Abs(cross) < TGeoShape::Tolerance()) return;
946  cross = 1./TMath::Sqrt(cross);
947  for (i=0; i<3; i++) norm[i] *= cross;
948 }
949
950 ////////////////////////////////////////////////////////////////////////////////
951 /// Fills array with n random points located on the surface of indexed facet.
952 /// The output array must be provided with a length of minimum 3*npoints. Returns
953 /// true if operation succeeded.
954 /// Possible index values:
955 /// - 0 - all facets together
956 /// - 1 to 6 - facet index from bottom to top Z
957
958 Bool_t TGeoArb8::GetPointsOnFacet(Int_t /*index*/, Int_t /*npoints*/, Double_t * /* array */) const
959 {
960  return kFALSE;
961 /*
962  if (index<0 || index>6) return kFALSE;
963  if (index==0) {
964  // Just generate same number of points on each facet
965  Int_t npts = npoints/6.;
966  Int_t count = 0;
967  for (Int_t ifacet=0; ifacet<6; ifacet++) {
968  if (GetPointsOnFacet(ifacet+1, npts, &array[3*count])) count += npts;
969  if (ifacet<5) npts = (npoints-count)/(5.-ifacet);
970  }
971  if (count>0) return kTRUE;
972  return kFALSE;
973  }
974  Double_t z, cf;
975  Double_t xmin=TGeoShape::Big();
976  Double_t xmax=-xmin;
977  Double_t ymin=TGeoShape::Big();
978  Double_t ymax=-ymin;
979  Double_t dy=0.;
980  Double_t poly[8];
981  Double_t point[2];
982  Int_t i;
983  if (index==1 || index==6) {
984  z = (index==1)?-fDz:fDz;
985  cf = 0.5*(fDz-z)/fDz;
986  for (i=0; i<4; i++) {
987  poly[2*i] = fXY[i+4][0]+cf*(fXY[i][0]-fXY[i+4][0]);
988  poly[2*i+1] = fXY[i+4][1]+cf*(fXY[i][1]-fXY[i+4][1]);
989  xmin = TMath::Min(xmin, poly[2*i]);
990  xmax = TMath::Max(xmax, poly[2*i]);
991  ymin = TMath::Min(ymin, poly[2*i]);
992  ymax = TMath::Max(ymax, poly[2*i]);
993  }
994  }
995  Int_t nshoot = 0;
996  Int_t nmiss = 0;
997  for (i=0; i<npoints; i++) {
998  Double_t *point = &array[3*i];
999  switch (surfindex) {
1000  case 1:
1001  case 6:
1002  while (nmiss<1000) {
1003  point[0] = xmin + (xmax-xmin)*gRandom->Rndm();
1004  point[1] = ymin + (ymax-ymin)*gRandom->Rndm();
1005  }
1006
1007  return InsidePolygon(point[0],point[1],poly);
1008 */
1009 }
1010
1011 ////////////////////////////////////////////////////////////////////////////////
1012 /// Finds if a point in XY plane is inside the polygon defines by PTS.
1013
1015 {
1016  Int_t i,j;
1017  Double_t x1,y1,x2,y2;
1018  Double_t cross;
1019  for (i=0; i<4; i++) {
1020  j = (i+1)%4;
1021  x1 = pts[i<<1];
1022  y1 = pts[(i<<1)+1];
1023  x2 = pts[j<<1];
1024  y2 = pts[(j<<1)+1];
1025  cross = (x-x1)*(y2-y1)-(y-y1)*(x2-x1);
1026  if (cross<0) return kFALSE;
1027  }
1028  return kTRUE;
1029 }
1030
1031 ////////////////////////////////////////////////////////////////////////////////
1032 /// Prints shape parameters
1033
1035 {
1036  printf("*** Shape %s: TGeoArb8 ***\n", GetName());
1037  if (IsTwisted()) printf(" = TWISTED\n");
1038  for (Int_t ip=0; ip<8; ip++) {
1039  printf(" point #%i : x=%11.5f y=%11.5f z=%11.5f\n",
1040  ip, fXY[ip][0], fXY[ip][1], fDz*((ip<4)?-1:1));
1041  }
1042  printf(" Bounding box:\n");
1044 }
1045
1046 ////////////////////////////////////////////////////////////////////////////////
1047 /// Computes the closest distance from given point to this shape.
1048
1049 Double_t TGeoArb8::Safety(const Double_t *point, Bool_t in) const
1050 {
1051  Double_t safz = fDz-TMath::Abs(point[2]);
1052  if (!in) safz = -safz;
1053  Int_t iseg;
1054  Double_t safe = TGeoShape::Big();
1055  Double_t lsq, ssq, dx, dy, dpx, dpy, u;
1056  if (IsTwisted()) {
1057  if (!in) {
1058  if (!TGeoBBox::Contains(point)) return TGeoBBox::Safety(point,kFALSE);
1059  }
1060  // Point is also in the bounding box ;-(
1061  // Compute closest distance to any segment
1062  Double_t vert[8];
1063  Double_t *p1, *p2;
1064  Int_t isegmin=0;
1065  Double_t umin = 0.;
1066  SetPlaneVertices (point[2], vert);
1067  for (iseg=0; iseg<4; iseg++) {
1068  if (safe<TGeoShape::Tolerance()) return 0.;
1069  p1 = &vert[2*iseg];
1070  p2 = &vert[2*((iseg+1)%4)];
1071  dx = p2[0] - p1[0];
1072  dy = p2[1] - p1[1];
1073  dpx = point[0] - p1[0];
1074  dpy = point[1] - p1[1];
1075
1076  lsq = dx*dx + dy*dy;
1077  u = (dpx*dx + dpy*dy)/lsq;
1078  if (u>1) {
1079  dpx = point[0]-p2[0];
1080  dpy = point[1]-p2[1];
1081  } else {
1082  if (u>=0) {
1083  dpx -= u*dx;
1084  dpy -= u*dy;
1085  }
1086  }
1087  ssq = dpx*dpx + dpy*dpy;
1088  if (ssq < safe) {
1089  isegmin = iseg;
1090  umin = u;
1091  safe = ssq;
1092  }
1093  }
1094  if (umin<0) umin = 0.;
1095  if (umin>1) {
1096  isegmin = (isegmin+1)%4;
1097  umin = 0.;
1098  }
1099  Int_t i1 = isegmin;
1100  Int_t i2 = (isegmin+1)%4;
1101  Double_t dx1 = fXY[i2][0]-fXY[i1][0];
1102  Double_t dx2 = fXY[i2+4][0]-fXY[i1+4][0];
1103  Double_t dy1 = fXY[i2][1]-fXY[i1][1];
1104  Double_t dy2 = fXY[i2+4][1]-fXY[i1+4][1];
1105  dx = dx1 + umin*(dx2-dx1);
1106  dy = dy1 + umin*(dy2-dy1);
1107  safe *= 1.- 4.*fDz*fDz/(dx*dx+dy*dy+4.*fDz*fDz);
1108  safe = TMath::Sqrt(safe);
1109  if (in) return TMath::Min(safz,safe);
1110  return TMath::Max(safz,safe);
1111  }
1112
1113  Double_t saf[5];
1114  saf[0] = safz;
1115
1116  for (iseg=0; iseg<4; iseg++) saf[iseg+1] = SafetyToFace(point,iseg,in);
1117  if (in) safe = saf[TMath::LocMin(5, saf)];
1118  else safe = saf[TMath::LocMax(5, saf)];
1119  if (safe<0) return 0.;
1120  return safe;
1121 }
1122
1123 ////////////////////////////////////////////////////////////////////////////////
1124 /// Estimate safety to lateral plane defined by segment iseg in range [0,3]
1125 /// Might be negative: plane seen only from inside.
1126
1127 Double_t TGeoArb8::SafetyToFace(const Double_t *point, Int_t iseg, Bool_t in) const
1128 {
1129  Double_t vertices[12];
1130  Int_t ipln = (iseg+1)%4;
1131  // point 1
1132  vertices[0] = fXY[iseg][0];
1133  vertices[1] = fXY[iseg][1];
1134  vertices[2] = -fDz;
1135  // point 2
1136  vertices[3] = fXY[ipln][0];
1137  vertices[4] = fXY[ipln][1];
1138  vertices[5] = -fDz;
1139  // point 3
1140  vertices[6] = fXY[ipln+4][0];
1141  vertices[7] = fXY[ipln+4][1];
1142  vertices[8] = fDz;
1143  // point 4
1144  vertices[9] = fXY[iseg+4][0];
1145  vertices[10] = fXY[iseg+4][1];
1146  vertices[11] = fDz;
1147  Double_t safe;
1148  Double_t norm[3];
1149  Double_t *p1, *p2, *p3;
1150  p1 = &vertices[0];
1151  p2 = &vertices[9];
1152  p3 = &vertices[6];
1153  if (IsSamePoint(p2,p3)) {
1154  p3 = &vertices[3];
1155  if (IsSamePoint(p1,p3)) return -TGeoShape::Big(); // skip single segment
1156  }
1157  GetPlaneNormal(p1,p2,p3,norm);
1158  safe = (point[0]-p1[0])*norm[0]+(point[1]-p1[1])*norm[1]+(point[2]-p1[2])*norm[2];
1159  if (in) return (-safe);
1160  return safe;
1161 }
1162
1163 ////////////////////////////////////////////////////////////////////////////////
1164 /// Save a primitive as a C++ statement(s) on output stream "out".
1165
1166 void TGeoArb8::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
1167 {
1168  if (TObject::TestBit(kGeoSavePrimitive)) return;
1169  out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
1170  out << " dz = " << fDz << ";" << std::endl;
1171  out << " vert[0] = " << fXY[0][0] << ";" << std::endl;
1172  out << " vert[1] = " << fXY[0][1] << ";" << std::endl;
1173  out << " vert[2] = " << fXY[1][0] << ";" << std::endl;
1174  out << " vert[3] = " << fXY[1][1] << ";" << std::endl;
1175  out << " vert[4] = " << fXY[2][0] << ";" << std::endl;
1176  out << " vert[5] = " << fXY[2][1] << ";" << std::endl;
1177  out << " vert[6] = " << fXY[3][0] << ";" << std::endl;
1178  out << " vert[7] = " << fXY[3][1] << ";" << std::endl;
1179  out << " vert[8] = " << fXY[4][0] << ";" << std::endl;
1180  out << " vert[9] = " << fXY[4][1] << ";" << std::endl;
1181  out << " vert[10] = " << fXY[5][0] << ";" << std::endl;
1182  out << " vert[11] = " << fXY[5][1] << ";" << std::endl;
1183  out << " vert[12] = " << fXY[6][0] << ";" << std::endl;
1184  out << " vert[13] = " << fXY[6][1] << ";" << std::endl;
1185  out << " vert[14] = " << fXY[7][0] << ";" << std::endl;
1186  out << " vert[15] = " << fXY[7][1] << ";" << std::endl;
1187  out << " TGeoShape *" << GetPointerName() << " = new TGeoArb8(\"" << GetName() << "\", dz,vert);" << std::endl;
1189 }
1190
1191 ////////////////////////////////////////////////////////////////////////////////
1192 /// Computes intersection points between plane at zpl and non-horizontal edges.
1193
1195 {
1196  Double_t cf = 0.5*(fDz-zpl)/fDz;
1197  for (Int_t i=0; i<4; i++) {
1198  vertices[2*i] = fXY[i+4][0]+cf*(fXY[i][0]-fXY[i+4][0]);
1199  vertices[2*i+1] = fXY[i+4][1]+cf*(fXY[i][1]-fXY[i+4][1]);
1200  }
1201 }
1202
1203 ////////////////////////////////////////////////////////////////////////////////
1204 /// Set all arb8 params in one step.
1205 /// param[0] = dz
1206 /// param[1] = x0
1207 /// param[2] = y0
1208 /// ...
1209
1211 {
1212  fDz = param[0];
1213  for (Int_t i=0; i<8; i++) {
1214  fXY[i][0] = param[2*i+1];
1215  fXY[i][1] = param[2*i+2];
1216  }
1217  ComputeTwist();
1218  ComputeBBox();
1219 }
1220
1221 ////////////////////////////////////////////////////////////////////////////////
1222 /// Creates arb8 mesh points
1223
1225 {
1226  for (Int_t i=0; i<8; i++) {
1227  points[3*i] = fXY[i][0];
1228  points[3*i+1] = fXY[i][1];
1229  points[3*i+2] = (i<4)?-fDz:fDz;
1230  }
1231 }
1232
1233 ////////////////////////////////////////////////////////////////////////////////
1234 /// Creates arb8 mesh points
1235
1237 {
1238  for (Int_t i=0; i<8; i++) {
1239  points[3*i] = fXY[i][0];
1240  points[3*i+1] = fXY[i][1];
1241  points[3*i+2] = (i<4)?-fDz:fDz;
1242  }
1243 }
1244
1245 ////////////////////////////////////////////////////////////////////////////////
1246 /// Set values for a given vertex.
1247
1249 {
1250  if (vnum<0 || vnum >7) {
1251  Error("SetVertex", "Invalid vertex number");
1252  return;
1253  }
1254  fXY[vnum][0] = x;
1255  fXY[vnum][1] = y;
1256  if (vnum == 7) {
1257  ComputeTwist();
1258  ComputeBBox();
1259  }
1260 }
1261
1262 ////////////////////////////////////////////////////////////////////////////////
1263 /// Fill size of this 3-D object
1264
1266 {
1268 }
1269
1270 ////////////////////////////////////////////////////////////////////////////////
1271 /// Stream an object of class TGeoManager.
1272
1273 void TGeoArb8::Streamer(TBuffer &R__b)
1274 {
1277  ComputeTwist();
1278  } else {
1279  R__b.WriteClassBuffer(TGeoArb8::Class(), this);
1280  }
1281 }
1282
1283 ////////////////////////////////////////////////////////////////////////////////
1284 /// Check the inside status for each of the points in the array.
1285 /// Input: Array of point coordinates + vector size
1286 /// Output: Array of Booleans for the inside of each point
1287
1288 void TGeoArb8::Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
1289 {
1290  for (Int_t i=0; i<vecsize; i++) inside[i] = Contains(&points[3*i]);
1291 }
1292
1293 ////////////////////////////////////////////////////////////////////////////////
1294 /// Compute the normal for an array o points so that norm.dot.dir is positive
1295 /// Input: Arrays of point coordinates and directions + vector size
1296 /// Output: Array of normal directions
1297
1298 void TGeoArb8::ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
1299 {
1300  for (Int_t i=0; i<vecsize; i++) ComputeNormal(&points[3*i], &dirs[3*i], &norms[3*i]);
1301 }
1302
1303 ////////////////////////////////////////////////////////////////////////////////
1304 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
1305
1306 void TGeoArb8::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
1307 {
1308  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromInside(&points[3*i], &dirs[3*i], 3, step[i]);
1309 }
1310
1311 ////////////////////////////////////////////////////////////////////////////////
1312 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
1313
1314 void TGeoArb8::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
1315 {
1316  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromOutside(&points[3*i], &dirs[3*i], 3, step[i]);
1317 }
1318
1319 ////////////////////////////////////////////////////////////////////////////////
1320 /// Compute safe distance from each of the points in the input array.
1321 /// Input: Array of point coordinates, array of statuses for these points, size of the arrays
1322 /// Output: Safety values
1323
1324 void TGeoArb8::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
1325 {
1326  for (Int_t i=0; i<vecsize; i++) safe[i] = Safety(&points[3*i], inside[i]);
1327 }
1328
1330
1331 ////////////////////////////////////////////////////////////////////////////////
1332 /// Default ctor
1333
1335 {
1336  fDz = 0;
1337  fTheta = 0;
1338  fPhi = 0;
1339  fH1 = fH2 = fBl1 = fBl2 = fTl1 = fTl2 = fAlpha1 = fAlpha2 = 0;
1340 }
1341
1342 ////////////////////////////////////////////////////////////////////////////////
1343 /// Constructor providing just a range in Z, theta and phi.
1344
1346  :TGeoArb8("", 0, 0)
1347 {
1348  fDz = dz;
1349  fTheta = theta;
1350  fPhi = phi;
1351  fH1 = fH2 = fBl1 = fBl2 = fTl1 = fTl2 = fAlpha1 = fAlpha2 = 0;
1352 }
1353
1354 ////////////////////////////////////////////////////////////////////////////////
1355 /// Normal constructor.
1356
1358  Double_t bl1, Double_t tl1, Double_t alpha1, Double_t h2, Double_t bl2,
1359  Double_t tl2, Double_t alpha2)
1360  :TGeoArb8("", 0, 0)
1361 {
1362  fDz = dz;
1363  fTheta = theta;
1364  fPhi = phi;
1365  fH1 = h1;
1366  fH2 = h2;
1367  fBl1 = bl1;
1368  fBl2 = bl2;
1369  fTl1 = tl1;
1370  fTl2 = tl2;
1371  fAlpha1 = alpha1;
1372  fAlpha2 = alpha2;
1377  fXY[0][0] = -dz*tx-h1*ta1-bl1; fXY[0][1] = -dz*ty-h1;
1378  fXY[1][0] = -dz*tx+h1*ta1-tl1; fXY[1][1] = -dz*ty+h1;
1379  fXY[2][0] = -dz*tx+h1*ta1+tl1; fXY[2][1] = -dz*ty+h1;
1380  fXY[3][0] = -dz*tx-h1*ta1+bl1; fXY[3][1] = -dz*ty-h1;
1381  fXY[4][0] = dz*tx-h2*ta2-bl2; fXY[4][1] = dz*ty-h2;
1382  fXY[5][0] = dz*tx+h2*ta2-tl2; fXY[5][1] = dz*ty+h2;
1383  fXY[6][0] = dz*tx+h2*ta2+tl2; fXY[6][1] = dz*ty+h2;
1384  fXY[7][0] = dz*tx-h2*ta2+bl2; fXY[7][1] = dz*ty-h2;
1385  ComputeTwist();
1386  if ((dz<0) || (h1<0) || (bl1<0) || (tl1<0) ||
1387  (h2<0) || (bl2<0) || (tl2<0)) {
1389  }
1390  else TGeoArb8::ComputeBBox();
1391 }
1392
1393 ////////////////////////////////////////////////////////////////////////////////
1394 /// Constructor with name.
1395
1397  Double_t bl1, Double_t tl1, Double_t alpha1, Double_t h2, Double_t bl2,
1398  Double_t tl2, Double_t alpha2)
1399  :TGeoArb8(name, 0, 0)
1400 {
1401  SetName(name);
1402  fDz = dz;
1403  fTheta = theta;
1404  fPhi = phi;
1405  fH1 = h1;
1406  fH2 = h2;
1407  fBl1 = bl1;
1408  fBl2 = bl2;
1409  fTl1 = tl1;
1410  fTl2 = tl2;
1411  fAlpha1 = alpha1;
1412  fAlpha2 = alpha2;
1413  for (Int_t i=0; i<8; i++) {
1414  fXY[i][0] = 0.0;
1415  fXY[i][1] = 0.0;
1416  }
1421  fXY[0][0] = -dz*tx-h1*ta1-bl1; fXY[0][1] = -dz*ty-h1;
1422  fXY[1][0] = -dz*tx+h1*ta1-tl1; fXY[1][1] = -dz*ty+h1;
1423  fXY[2][0] = -dz*tx+h1*ta1+tl1; fXY[2][1] = -dz*ty+h1;
1424  fXY[3][0] = -dz*tx-h1*ta1+bl1; fXY[3][1] = -dz*ty-h1;
1425  fXY[4][0] = dz*tx-h2*ta2-bl2; fXY[4][1] = dz*ty-h2;
1426  fXY[5][0] = dz*tx+h2*ta2-tl2; fXY[5][1] = dz*ty+h2;
1427  fXY[6][0] = dz*tx+h2*ta2+tl2; fXY[6][1] = dz*ty+h2;
1428  fXY[7][0] = dz*tx-h2*ta2+bl2; fXY[7][1] = dz*ty-h2;
1429  ComputeTwist();
1430  if ((dz<0) || (h1<0) || (bl1<0) || (tl1<0) ||
1431  (h2<0) || (bl2<0) || (tl2<0)) {
1433  }
1434  else TGeoArb8::ComputeBBox();
1435 }
1436
1437 ////////////////////////////////////////////////////////////////////////////////
1438 /// Destructor.
1439
1441 {
1442 }
1443
1444 ////////////////////////////////////////////////////////////////////////////////
1445 /// Compute distance from inside point to surface of the trapezoid
1446
1447 Double_t TGeoTrap::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
1448 {
1449  if (iact<3 && safe) {
1450  // compute safe distance
1451  *safe = Safety(point, kTRUE);
1452  if (iact==0) return TGeoShape::Big();
1453  if (iact==1 && step<*safe) return TGeoShape::Big();
1454  }
1455  // compute distance to get outside this shape
1456  // return TGeoArb8::DistFromInside(point, dir, iact, step, safe);
1457
1458 // compute distance to plane ipl :
1459 // ipl=0 : points 0,4,1,5
1460 // ipl=1 : points 1,5,2,6
1461 // ipl=2 : points 2,6,3,7
1462 // ipl=3 : points 3,7,0,4
1463  Double_t distmin;
1464  if (dir[2]<0) {
1465  distmin=(-fDz-point[2])/dir[2];
1466  } else {
1467  if (dir[2]>0) distmin =(fDz-point[2])/dir[2];
1468  else distmin = TGeoShape::Big();
1469  }
1470  Double_t xa,xb,xc;
1471  Double_t ya,yb,yc;
1472  for (Int_t ipl=0;ipl<4;ipl++) {
1473  Int_t j = (ipl+1)%4;
1474  xa=fXY[ipl][0];
1475  ya=fXY[ipl][1];
1476  xb=fXY[ipl+4][0];
1477  yb=fXY[ipl+4][1];
1478  xc=fXY[j][0];
1479  yc=fXY[j][1];
1480  Double_t ax,ay,az;
1481  ax = xb-xa;
1482  ay = yb-ya;
1483  az = 2.*fDz;
1484  Double_t bx,by;
1485  bx = xc-xa;
1486  by = yc-ya;
1487  Double_t ddotn = -dir[0]*az*by + dir[1]*az*bx+dir[2]*(ax*by-ay*bx);
1488  if (ddotn<=0) continue; // entering
1489  Double_t saf = -(point[0]-xa)*az*by + (point[1]-ya)*az*bx + (point[2]+fDz)*(ax*by-ay*bx);
1490  if (saf>=0.0) return 0.0;
1491  Double_t s = -saf/ddotn;
1492  if (s<distmin) distmin=s;
1493  }
1494  return distmin;
1495 }
1496
1497 ////////////////////////////////////////////////////////////////////////////////
1498 /// Compute distance from outside point to surface of the trapezoid
1499
1500 Double_t TGeoTrap::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
1501 {
1502  if (iact<3 && safe) {
1503  // compute safe distance
1504  *safe = Safety(point, kFALSE);
1505  if (iact==0) return TGeoShape::Big();
1506  if (iact==1 && step<*safe) return TGeoShape::Big();
1507  }
1508  // Check if the bounding box is crossed within the requested distance
1509  Double_t sdist = TGeoBBox::DistFromOutside(point,dir, fDX, fDY, fDZ, fOrigin, step);
1510  if (sdist>=step) return TGeoShape::Big();
1511  // compute distance to get outside this shape
1512  Bool_t in = kTRUE;
1513  Double_t pts[8];
1514  Double_t snxt;
1515  Double_t xnew, ynew, znew;
1516  Int_t i,j;
1517  if (point[2]<-fDz+TGeoShape::Tolerance()) {
1518  if (dir[2]<=0) return TGeoShape::Big();
1519  in = kFALSE;
1520  snxt = -(fDz+point[2])/dir[2];
1521  xnew = point[0] + snxt*dir[0];
1522  ynew = point[1] + snxt*dir[1];
1523  for (i=0;i<4;i++) {
1524  j = i<<1;
1525  pts[j] = fXY[i][0];
1526  pts[j+1] = fXY[i][1];
1527  }
1528  if (InsidePolygon(xnew,ynew,pts)) return snxt;
1529  } else if (point[2]>fDz-TGeoShape::Tolerance()) {
1530  if (dir[2]>=0) return TGeoShape::Big();
1531  in = kFALSE;
1532  snxt = (fDz-point[2])/dir[2];
1533  xnew = point[0] + snxt*dir[0];
1534  ynew = point[1] + snxt*dir[1];
1535  for (i=0;i<4;i++) {
1536  j = i<<1;
1537  pts[j] = fXY[i+4][0];
1538  pts[j+1] = fXY[i+4][1];
1539  }
1540  if (InsidePolygon(xnew,ynew,pts)) return snxt;
1541  }
1542  snxt = TGeoShape::Big();
1543
1544
1545  // check lateral faces
1546  Double_t dz2 =0.5/fDz;
1547  Double_t xa,xb,xc,xd;
1548  Double_t ya,yb,yc,yd;
1549  Double_t ax,ay,az;
1550  Double_t bx,by;
1551  Double_t ddotn, saf;
1552  Double_t safmin = TGeoShape::Big();
1553  Bool_t exiting = kFALSE;
1554  for (i=0; i<4; i++) {
1555  j = (i+1)%4;
1556  xa=fXY[i][0];
1557  ya=fXY[i][1];
1558  xb=fXY[i+4][0];
1559  yb=fXY[i+4][1];
1560  xc=fXY[j][0];
1561  yc=fXY[j][1];
1562  xd=fXY[4+j][0];
1563  yd=fXY[4+j][1];
1564  ax = xb-xa;
1565  ay = yb-ya;
1566  az = 2.*fDz;
1567  bx = xc-xa;
1568  by = yc-ya;
1569  ddotn = -dir[0]*az*by + dir[1]*az*bx+dir[2]*(ax*by-ay*bx);
1570  saf = (point[0]-xa)*az*by - (point[1]-ya)*az*bx - (point[2]+fDz)*(ax*by-ay*bx);
1571
1572  if (saf<=0) {
1573  // face visible from point outside
1574  in = kFALSE;
1575  if (ddotn>=0) return TGeoShape::Big();
1576  snxt = saf/ddotn;
1577  znew = point[2]+snxt*dir[2];
1578  if (TMath::Abs(znew)<=fDz) {
1579  xnew = point[0]+snxt*dir[0];
1580  ynew = point[1]+snxt*dir[1];
1581  Double_t tx1 =dz2*(xb-xa);
1582  Double_t ty1 =dz2*(yb-ya);
1583  Double_t tx2 =dz2*(xd-xc);
1584  Double_t ty2 =dz2*(yd-yc);
1585  Double_t dzp =fDz+znew;
1586  Double_t xs1 =xa+tx1*dzp;
1587  Double_t ys1 =ya+ty1*dzp;
1588  Double_t xs2 =xc+tx2*dzp;
1589  Double_t ys2 =yc+ty2*dzp;
1590  if (TMath::Abs(xs1-xs2)>TMath::Abs(ys1-ys2)) {
1591  if ((xnew-xs1)*(xs2-xnew)>=0) return snxt;
1592  } else {
1593  if ((ynew-ys1)*(ys2-ynew)>=0) return snxt;
1594  }
1595  }
1596  } else {
1597  if (saf<safmin) {
1598  safmin = saf;
1599  if (ddotn>=0) exiting = kTRUE;
1600  else exiting = kFALSE;
1601  }
1602  }
1603  }
1604  // Check also Z boundaries (point may be inside and close to Z) - Christian Hammann
1605  saf = fDz - TMath::Abs(point[2]);
1606  if (saf>0 && saf<safmin) exiting = (point[2]*dir[2] > 0)?kTRUE:kFALSE;
1607  if (!in) return TGeoShape::Big();
1608  if (exiting) return TGeoShape::Big();
1609  return 0.0;
1610 }
1611
1612 ////////////////////////////////////////////////////////////////////////////////
1613 /// Divide this trapezoid shape belonging to volume "voldiv" into ndiv volumes
1614 /// called divname, from start position with the given step. Only Z divisions
1615 /// are supported. For Z divisions just return the pointer to the volume to be
1616 /// divided. In case a wrong division axis is supplied, returns pointer to
1617 /// volume that was divided.
1618
1619 TGeoVolume *TGeoTrap::Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv,
1620  Double_t start, Double_t step)
1621 {
1622  TGeoShape *shape; //--- shape to be created
1623  TGeoVolume *vol; //--- division volume to be created
1624  TGeoVolumeMulti *vmulti; //--- generic divided volume
1625  TGeoPatternFinder *finder; //--- finder to be attached
1626  TString opt = ""; //--- option to be attached
1627  if (iaxis!=3) {
1628  Error("Divide", "cannot divide trapezoids on other axis than Z");
1629  return 0;
1630  }
1631  Double_t end = start+ndiv*step;
1632  Double_t points_lo[8];
1633  Double_t points_hi[8];
1634  finder = new TGeoPatternTrapZ(voldiv, ndiv, start, end);
1635  voldiv->SetFinder(finder);
1636  finder->SetDivIndex(voldiv->GetNdaughters());
1637  opt = "Z";
1638  vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
1639  Double_t txz = ((TGeoPatternTrapZ*)finder)->GetTxz();
1640  Double_t tyz = ((TGeoPatternTrapZ*)finder)->GetTyz();
1641  Double_t zmin, zmax, ox,oy,oz;
1642  for (Int_t idiv=0; idiv<ndiv; idiv++) {
1643  zmin = start+idiv*step;
1644  zmax = start+(idiv+1)*step;
1645  oz = start+idiv*step+step/2;
1646  ox = oz*txz;
1647  oy = oz*tyz;
1648  SetPlaneVertices(zmin, &points_lo[0]);
1649  SetPlaneVertices(zmax, &points_hi[0]);
1650  shape = new TGeoTrap(step/2, fTheta, fPhi);
1651  for (Int_t vert1=0; vert1<4; vert1++)
1652  ((TGeoArb8*)shape)->SetVertex(vert1, points_lo[2*vert1]-ox, points_lo[2*vert1+1]-oy);
1653  for (Int_t vert2=0; vert2<4; vert2++)
1654  ((TGeoArb8*)shape)->SetVertex(vert2+4, points_hi[2*vert2]-ox, points_hi[2*vert2+1]-oy);
1655  vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
1658  ((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder);
1659  }
1660  return vmulti;
1661 }
1662
1663 ////////////////////////////////////////////////////////////////////////////////
1664 /// In case shape has some negative parameters, these have to be computed
1665 /// in order to fit the mother.
1666
1668 {
1669  if (!TestShapeBit(kGeoRunTimeShape)) return 0;
1670  if (mother->IsRunTimeShape()) {
1671  Error("GetMakeRuntimeShape", "invalid mother");
1672  return 0;
1673  }
1674  Double_t dz, h1, bl1, tl1, h2, bl2, tl2;
1675  if (fDz<0) dz=((TGeoTrap*)mother)->GetDz();
1676  else dz=fDz;
1677
1678  if (fH1<0) h1 = ((TGeoTrap*)mother)->GetH1();
1679  else h1 = fH1;
1680
1681  if (fH2<0) h2 = ((TGeoTrap*)mother)->GetH2();
1682  else h2 = fH2;
1683
1684  if (fBl1<0) bl1 = ((TGeoTrap*)mother)->GetBl1();
1685  else bl1 = fBl1;
1686
1687  if (fBl2<0) bl2 = ((TGeoTrap*)mother)->GetBl2();
1688  else bl2 = fBl2;
1689
1690  if (fTl1<0) tl1 = ((TGeoTrap*)mother)->GetTl1();
1691  else tl1 = fTl1;
1692
1693  if (fTl2<0) tl2 = ((TGeoTrap*)mother)->GetTl2();
1694  else tl2 = fTl2;
1695
1696  return (new TGeoTrap(dz, fTheta, fPhi, h1, bl1, tl1, fAlpha1, h2, bl2, tl2, fAlpha2));
1697 }
1698
1699 ////////////////////////////////////////////////////////////////////////////////
1700 /// Computes the closest distance from given point to this shape.
1701
1702 Double_t TGeoTrap::Safety(const Double_t *point, Bool_t in) const
1703 {
1704  Double_t safe = TGeoShape::Big();
1705  Double_t saf[5];
1706  Double_t norm[3]; // normal to current facette
1707  Int_t i, j; // current facette index
1708  Double_t x0, y0, z0=-fDz, x1, y1, z1=fDz, x2, y2;
1709  Double_t ax, ay, az=z1-z0, bx, by;
1710  Double_t fn;
1711  //---> compute safety for lateral planes
1712  for (i=0; i<4; i++) {
1713  if (in) saf[i] = TGeoShape::Big();
1714  else saf[i] = 0.;
1715  x0 = fXY[i][0];
1716  y0 = fXY[i][1];
1717  x1 = fXY[i+4][0];
1718  y1 = fXY[i+4][1];
1719  ax = x1-x0;
1720  ay = y1-y0;
1721  az = z1-z0;
1722  j = (i+1)%4;
1723  x2 = fXY[j][0];
1724  y2 = fXY[j][1];
1725  bx = x2-x0;
1726  by = y2-y0;
1728  x2 = fXY[4+j][0];
1729  y2 = fXY[4+j][1];
1730  bx = x2-x1;
1731  by = y2-y1;
1732  if (TMath::Abs(bx)<TGeoShape::Tolerance() && TMath::Abs(by)<TGeoShape::Tolerance()) continue;
1733  }
1734  norm[0] = -az*by;
1735  norm[1] = az*bx;
1736  norm[2] = ax*by-ay*bx;
1737  fn = TMath::Sqrt(norm[0]*norm[0]+norm[1]*norm[1]+norm[2]*norm[2]);
1738  if (fn<1E-10) continue;
1739  saf[i] = (x0-point[0])*norm[0]+(y0-point[1])*norm[1]+(-fDz-point[2])*norm[2];
1740  if (in) {
1741  saf[i]=TMath::Abs(saf[i])/fn; // they should be all positive anyway
1742  } else {
1743  saf[i] = -saf[i]/fn; // only negative values are interesting
1744  }
1745  }
1746  saf[4] = fDz-TMath::Abs(point[2]);
1747  if (in) {
1748  safe = saf[0];
1749  for (j=1;j<5;j++) if (saf[j] <safe) safe = saf[j];
1750  } else {
1751  saf[4]=-saf[4];
1752  safe = saf[0];
1753  for (j=1;j<5;j++) if (saf[j] >safe) safe = saf[j];
1754  }
1755  return safe;
1756 }
1757
1758 ////////////////////////////////////////////////////////////////////////////////
1759 /// Save a primitive as a C++ statement(s) on output stream "out".
1760
1761 void TGeoTrap::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
1762 {
1763  if (TObject::TestBit(kGeoSavePrimitive)) return;
1764  out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
1765  out << " dz = " << fDz << ";" << std::endl;
1766  out << " theta = " << fTheta << ";" << std::endl;
1767  out << " phi = " << fPhi << ";" << std::endl;
1768  out << " h1 = " << fH1<< ";" << std::endl;
1769  out << " bl1 = " << fBl1<< ";" << std::endl;
1770  out << " tl1 = " << fTl1<< ";" << std::endl;
1771  out << " alpha1 = " << fAlpha1 << ";" << std::endl;
1772  out << " h2 = " << fH2 << ";" << std::endl;
1773  out << " bl2 = " << fBl2<< ";" << std::endl;
1774  out << " tl2 = " << fTl2<< ";" << std::endl;
1775  out << " alpha2 = " << fAlpha2 << ";" << std::endl;
1776  out << " TGeoShape *" << GetPointerName() << " = new TGeoTrap(\"" << GetName() << "\", dz,theta,phi,h1,bl1,tl1,alpha1,h2,bl2,tl2,alpha2);" << std::endl;
1778 }
1779
1780 ////////////////////////////////////////////////////////////////////////////////
1781 /// Set all arb8 params in one step.
1782 /// - param[0] = dz
1783 /// - param[1] = theta
1784 /// - param[2] = phi
1785 /// - param[3] = h1
1786 /// - param[4] = bl1
1787 /// - param[5] = tl1
1788 /// - param[6] = alpha1
1789 /// - param[7] = h2
1790 /// - param[8] = bl2
1791 /// - param[9] = tl2
1792 /// - param[10] = alpha2
1793
1795 {
1796  fDz = param[0];
1797  fTheta = param[1];
1798  fPhi = param[2];
1799  fH1 = param[3];
1800  fH2 = param[7];
1801  fBl1 = param[4];
1802  fBl2 = param[8];
1803  fTl1 = param[5];
1804  fTl2 = param[9];
1805  fAlpha1 = param[6];
1806  fAlpha2 = param[10];
1811  fXY[0][0] = -fDz*tx-fH1*ta1-fBl1; fXY[0][1] = -fDz*ty-fH1;
1812  fXY[1][0] = -fDz*tx+fH1*ta1-fTl1; fXY[1][1] = -fDz*ty+fH1;
1813  fXY[2][0] = -fDz*tx+fH1*ta1+fTl1; fXY[2][1] = -fDz*ty+fH1;
1814  fXY[3][0] = -fDz*tx-fH1*ta1+fBl1; fXY[3][1] = -fDz*ty-fH1;
1815  fXY[4][0] = fDz*tx-fH2*ta2-fBl2; fXY[4][1] = fDz*ty-fH2;
1816  fXY[5][0] = fDz*tx+fH2*ta2-fTl2; fXY[5][1] = fDz*ty+fH2;
1817  fXY[6][0] = fDz*tx+fH2*ta2+fTl2; fXY[6][1] = fDz*ty+fH2;
1818  fXY[7][0] = fDz*tx-fH2*ta2+fBl2; fXY[7][1] = fDz*ty-fH2;
1819  ComputeTwist();
1820  if ((fDz<0) || (fH1<0) || (fBl1<0) || (fTl1<0) ||
1821  (fH2<0) || (fBl2<0) || (fTl2<0)) {
1823  }
1824  else TGeoArb8::ComputeBBox();
1825 }
1826
1827 ////////////////////////////////////////////////////////////////////////////////
1828 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
1829
1830 void TGeoTrap::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
1831 {
1832  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromInside(&points[3*i], &dirs[3*i], 3, step[i]);
1833 }
1834
1835 ////////////////////////////////////////////////////////////////////////////////
1836 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
1837
1838 void TGeoTrap::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
1839 {
1840  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromOutside(&points[3*i], &dirs[3*i], 3, step[i]);
1841 }
1842
1843 ////////////////////////////////////////////////////////////////////////////////
1844 /// Compute safe distance from each of the points in the input array.
1845 /// Input: Array of point coordinates, array of statuses for these points, size of the arrays
1846 /// Output: Safety values
1847
1848 void TGeoTrap::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
1849 {
1850  for (Int_t i=0; i<vecsize; i++) safe[i] = Safety(&points[3*i], inside[i]);
1851 }
1852
1854
1855 ////////////////////////////////////////////////////////////////////////////////
1856 /// Default ctor
1857
1859 {
1860  fTwistAngle = 0;
1861
1862 }
1863
1864 ////////////////////////////////////////////////////////////////////////////////
1865 /// Constructor.
1866
1868  Double_t bl1, Double_t tl1, Double_t alpha1, Double_t h2, Double_t bl2,
1869  Double_t tl2, Double_t alpha2)
1870  :TGeoTrap(dz, theta, phi, h1, bl1, tl1, alpha1, h2, bl2, tl2, alpha2)
1871 {
1872  fTwistAngle = twist;
1873  Double_t x,y;
1876  // Coordinates of the center of the bottom face
1877  Double_t xc = -dz*TMath::Sin(th)*TMath::Cos(ph);
1878  Double_t yc = -dz*TMath::Sin(th)*TMath::Sin(ph);
1879
1880  Int_t i;
1881
1882  for (i=0; i<4; i++) {
1883  x = fXY[i][0] - xc;
1884  y = fXY[i][1] - yc;
1887  }
1888  // Coordinates of the center of the top face
1889  xc = -xc;
1890  yc = -yc;
1891  for (i=4; i<8; i++) {
1892  x = fXY[i][0] - xc;
1893  y = fXY[i][1] - yc;
1896  }
1897  ComputeTwist();
1898  if ((dz<0) || (h1<0) || (bl1<0) || (tl1<0) ||
1899  (h2<0) || (bl2<0) || (tl2<0)) SetShapeBit(kGeoRunTimeShape);
1900  else TGeoArb8::ComputeBBox();
1901 }
1902
1903 ////////////////////////////////////////////////////////////////////////////////
1904 /// Constructor providing the name of the shape.
1905
1906 TGeoGtra::TGeoGtra(const char *name, Double_t dz, Double_t theta, Double_t phi, Double_t twist, Double_t h1,
1907  Double_t bl1, Double_t tl1, Double_t alpha1, Double_t h2, Double_t bl2,
1908  Double_t tl2, Double_t alpha2)
1909  :TGeoTrap(name, dz, theta, phi, h1, bl1, tl1, alpha1, h2, bl2, tl2, alpha2)
1910 {
1911  fTwistAngle = twist;
1912  Double_t x,y;
1915  // Coordinates of the center of the bottom face
1916  Double_t xc = -dz*TMath::Sin(th)*TMath::Cos(ph);
1917  Double_t yc = -dz*TMath::Sin(th)*TMath::Sin(ph);
1918
1919  Int_t i;
1920
1921  for (i=0; i<4; i++) {
1922  x = fXY[i][0] - xc;
1923  y = fXY[i][1] - yc;
1926  }
1927  // Coordinates of the center of the top face
1928  xc = -xc;
1929  yc = -yc;
1930  for (i=4; i<8; i++) {
1931  x = fXY[i][0] - xc;
1932  y = fXY[i][1] - yc;
1935  }
1936  ComputeTwist();
1937  if ((dz<0) || (h1<0) || (bl1<0) || (tl1<0) ||
1938  (h2<0) || (bl2<0) || (tl2<0)) SetShapeBit(kGeoRunTimeShape);
1939  else TGeoArb8::ComputeBBox();
1940 }
1941
1942 ////////////////////////////////////////////////////////////////////////////////
1943 /// Destructor.
1944
1946 {
1947 }
1948
1949 ////////////////////////////////////////////////////////////////////////////////
1950 /// Compute distance from inside point to surface of the shape.
1951
1952 Double_t TGeoGtra::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
1953 {
1954  if (iact<3 && safe) {
1955  // compute safe distance
1956  *safe = Safety(point, kTRUE);
1957  if (iact==0) return TGeoShape::Big();
1958  if (iact==1 && step<*safe) return TGeoShape::Big();
1959  }
1960  // compute distance to get outside this shape
1961  return TGeoArb8::DistFromInside(point, dir, iact, step, safe);
1962 }
1963
1964 ////////////////////////////////////////////////////////////////////////////////
1965 /// Compute distance from inside point to surface of the shape.
1966
1967 Double_t TGeoGtra::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
1968 {
1969  if (iact<3 && safe) {
1970  // compute safe distance
1971  *safe = Safety(point, kTRUE);
1972  if (iact==0) return TGeoShape::Big();
1973  if (iact==1 && step<*safe) return TGeoShape::Big();
1974  }
1975  // compute distance to get outside this shape
1976  return TGeoArb8::DistFromOutside(point, dir, iact, step, safe);
1977 }
1978
1979 ////////////////////////////////////////////////////////////////////////////////
1980 /// In case shape has some negative parameters, these has to be computed
1981 /// in order to fit the mother
1982
1984 {
1985  if (!TestShapeBit(kGeoRunTimeShape)) return 0;
1986  if (mother->IsRunTimeShape()) {
1987  Error("GetMakeRuntimeShape", "invalid mother");
1988  return 0;
1989  }
1990  Double_t dz, h1, bl1, tl1, h2, bl2, tl2;
1991  if (fDz<0) dz=((TGeoTrap*)mother)->GetDz();
1992  else dz=fDz;
1993  if (fH1<0)
1994  h1 = ((TGeoTrap*)mother)->GetH1();
1995  else
1996  h1 = fH1;
1997  if (fH2<0)
1998  h2 = ((TGeoTrap*)mother)->GetH2();
1999  else
2000  h2 = fH2;
2001  if (fBl1<0)
2002  bl1 = ((TGeoTrap*)mother)->GetBl1();
2003  else
2004  bl1 = fBl1;
2005  if (fBl2<0)
2006  bl2 = ((TGeoTrap*)mother)->GetBl2();
2007  else
2008  bl2 = fBl2;
2009  if (fTl1<0)
2010  tl1 = ((TGeoTrap*)mother)->GetTl1();
2011  else
2012  tl1 = fTl1;
2013  if (fTl2<0)
2014  tl2 = ((TGeoTrap*)mother)->GetTl2();
2015  else
2016  tl2 = fTl2;
2017  return (new TGeoGtra(dz, fTheta, fPhi, fTwistAngle ,h1, bl1, tl1, fAlpha1, h2, bl2, tl2, fAlpha2));
2018 }
2019
2020 ////////////////////////////////////////////////////////////////////////////////
2021 /// Computes the closest distance from given point to this shape.
2022
2023 Double_t TGeoGtra::Safety(const Double_t *point, Bool_t in) const
2024 {
2025  return TGeoArb8::Safety(point,in);
2026 }
2027
2028 ////////////////////////////////////////////////////////////////////////////////
2029 /// Save a primitive as a C++ statement(s) on output stream "out".
2030
2031 void TGeoGtra::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
2032 {
2033  if (TObject::TestBit(kGeoSavePrimitive)) return;
2034  out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
2035  out << " dz = " << fDz << ";" << std::endl;
2036  out << " theta = " << fTheta << ";" << std::endl;
2037  out << " phi = " << fPhi << ";" << std::endl;
2038  out << " twist = " << fTwistAngle << ";" << std::endl;
2039  out << " h1 = " << fH1<< ";" << std::endl;
2040  out << " bl1 = " << fBl1<< ";" << std::endl;
2041  out << " tl1 = " << fTl1<< ";" << std::endl;
2042  out << " alpha1 = " << fAlpha1 << ";" << std::endl;
2043  out << " h2 = " << fH2 << ";" << std::endl;
2044  out << " bl2 = " << fBl2<< ";" << std::endl;
2045  out << " tl2 = " << fTl2<< ";" << std::endl;
2046  out << " alpha2 = " << fAlpha2 << ";" << std::endl;
2047  out << " TGeoShape *" << GetPointerName() << " = new TGeoGtra(\"" << GetName() << "\", dz,theta,phi,twist,h1,bl1,tl1,alpha1,h2,bl2,tl2,alpha2);" << std::endl;
2049 }
2050
2051 ////////////////////////////////////////////////////////////////////////////////
2052 /// Set all arb8 params in one step.
2053 /// - param[0] = dz
2054 /// - param[1] = theta
2055 /// - param[2] = phi
2056 /// - param[3] = h1
2057 /// - param[4] = bl1
2058 /// - param[5] = tl1
2059 /// - param[6] = alpha1
2060 /// - param[7] = h2
2061 /// - param[8] = bl2
2062 /// - param[9] = tl2
2063 /// - param[10] = alpha2
2064 /// - param[11] = twist
2065
2067 {
2068  TGeoTrap::SetDimensions(param);
2069  fTwistAngle = param[11];
2070  Double_t x,y;
2071  Double_t twist = fTwistAngle;
2074  // Coordinates of the center of the bottom face
2075  Double_t xc = -fDz*TMath::Sin(th)*TMath::Cos(ph);
2076  Double_t yc = -fDz*TMath::Sin(th)*TMath::Sin(ph);
2077
2078  Int_t i;
2079
2080  for (i=0; i<4; i++) {
2081  x = fXY[i][0] - xc;
2082  y = fXY[i][1] - yc;
2085  }
2086  // Coordinates of the center of the top face
2087  xc = -xc;
2088  yc = -yc;
2089  for (i=4; i<8; i++) {
2090  x = fXY[i][0] - xc;
2091  y = fXY[i][1] - yc;
2094  }
2095  ComputeTwist();
2096  if ((fDz<0) || (fH1<0) || (fBl1<0) || (fTl1<0) ||
2097  (fH2<0) || (fBl2<0) || (fTl2<0)) SetShapeBit(kGeoRunTimeShape);
2098  else TGeoArb8::ComputeBBox();
2099 }
2100
2101 ////////////////////////////////////////////////////////////////////////////////
2102 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
2103
2104 void TGeoGtra::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
2105 {
2106  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromInside(&points[3*i], &dirs[3*i], 3, step[i]);
2107 }
2108
2109 ////////////////////////////////////////////////////////////////////////////////
2110 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
2111
2112 void TGeoGtra::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
2113 {
2114  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromOutside(&points[3*i], &dirs[3*i], 3, step[i]);
2115 }
2116
2117 ////////////////////////////////////////////////////////////////////////////////
2118 /// Compute safe distance from each of the points in the input array.
2119 /// Input: Array of point coordinates, array of statuses for these points, size of the arrays
2120 /// Output: Safety values
2121
2122 void TGeoGtra::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
2123 {
2124  for (Int_t i=0; i<vecsize; i++) safe[i] = Safety(&points[3*i], inside[i]);
2125 }
TGeoVolumeMulti * MakeVolumeMulti(const char *name, TGeoMedium *medium)
Make a TGeoVolumeMulti handling a list of volumes.
virtual Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=0) const
Computes distance from outside point to surface of the shape.
Definition: TGeoArb8.cxx:672
Double_t GetClosestEdge(const Double_t *point, Double_t *vert, Int_t &isegment) const
Get index of the edge of the quadrilater represented by vert closest to point.
Definition: TGeoArb8.cxx:395
Definition: TBuffer.h:83
Volume families.
Definition: TGeoVolume.h:256
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
float xmin
Definition: THbookFile.cxx:93
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
#define snext(osub1, osub2)
Definition: triangle.c:1167
Box class.
Definition: TGeoBBox.h:17
Long64_t LocMax(Long64_t n, const T *a)
Return index of array with the maximum element.
Definition: TMath.h:999
static double p3(double t, double a, double b, double c, double d)
virtual ~TGeoGtra()
Destructor.
Definition: TGeoArb8.cxx:1945
T1 Sign(T1 a, T2 b)
Definition: TMathBase.h:153
void SetFinder(TGeoPatternFinder *finder)
Definition: TGeoVolume.h:234
virtual TGeoVolume * Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv, Double_t start, Double_t step)
Divide this shape along one axis.
Definition: TGeoArb8.cxx:826
virtual Double_t GetDX() const
Definition: TGeoBBox.h:70
float Float_t
Definition: RtypesCore.h:53
void AddNodeOffset(TGeoVolume *vol, Int_t copy_no, Double_t offset=0, Option_t *option="")
Add a division node to the list of nodes.
virtual void GetBoundingCylinder(Double_t *param) const
Fill vector param[4] with the bounding cylinder parameters.
Definition: TGeoArb8.cxx:854
Gtra is a twisted trapezoid.
Definition: TGeoArb8.h:143
const char Option_t
Definition: RtypesCore.h:62
virtual Double_t Capacity() const
Computes capacity of the shape in [length^3].
Definition: TGeoArb8.cxx:248
Geometrical transformation package.
Definition: TGeoMatrix.h:40
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save a primitive as a C++ statement(s) on output stream "out".
Definition: TGeoArb8.cxx:1761
float ymin
Definition: THbookFile.cxx:93
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition: TNamed.cxx:140
virtual void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
Compute safe distance from each of the points in the input array.
Definition: TGeoArb8.cxx:2122
TGeoVolume, TGeoVolumeMulti, TGeoVolumeAssembly are the volume classes.
Definition: TGeoVolume.h:48
void SetPlaneVertices(Double_t zpl, Double_t *vertices) const
Computes intersection points between plane at zpl and non-horizontal edges.
Definition: TGeoArb8.cxx:1194
virtual Bool_t Contains(const Double_t *point) const
Test if point is inside this shape.
Definition: TGeoArb8.cxx:527
Buffer base class used for serializing objects.
Definition: TBuffer.h:40
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:172
Double_t fOrigin[3]
Definition: TGeoBBox.h:24
Double_t fBl1
Definition: TGeoArb8.h:96
Basic string class.
Definition: TString.h:131
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:168
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
virtual Double_t GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const
Get shape range on a given axis.
Definition: TGeoArb8.cxx:836
virtual Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=0) const
Compute distance from inside point to surface of the trapezoid.
Definition: TGeoArb8.cxx:1447
virtual void InspectShape() const
Prints shape parameters.
Definition: TGeoBBox.cxx:793
TObject * At(Int_t idx) const
Definition: TObjArray.h:165
virtual void DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const
Compute distance from array of input points having directions specified by dirs. Store output in dist...
Definition: TGeoArb8.cxx:1306
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
virtual void InspectShape() const
Prints shape parameters.
Definition: TGeoArb8.cxx:1034
static Bool_t IsSameWithinTolerance(Double_t a, Double_t b)
Check if two numbers differ with less than a tolerance.
Definition: TGeoShape.cxx:326
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
TObjArray * GetNodes()
Definition: TGeoVolume.h:170
static Double_t Tolerance()
Definition: TGeoShape.h:91
virtual void Sizeof3D() const
Fill size of this 3-D object.
Definition: TGeoArb8.cxx:1265
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save a primitive as a C++ statement(s) on output stream "out".
Definition: TGeoArb8.cxx:2031
Int_t GetNdaughters() const
Definition: TGeoVolume.h:350
Double_t SafetyToFace(const Double_t *point, Int_t iseg, Bool_t in) const
Estimate safety to lateral plane defined by segment iseg in range [0,3] Might be negative: plane seen...
Definition: TGeoArb8.cxx:1127
Double_t GetTwist(Int_t iseg) const
Get twist for segment I in range [0,3].
Definition: TGeoArb8.cxx:383
static const double x2[5]
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:128
Double_t x[n]
Definition: legend1.C:17
Double_t fDZ
Definition: TGeoBBox.h:23
Bool_t IsTwisted() const
Definition: TGeoArb8.h:73
void Class()
Definition: Class.C:29
Add a volume with valid shape to the list of volumes.
static double p2(double t, double a, double b, double c)
Conversion from degree to radian: .
Definition: TMath.h:82
TGeoMedium * GetMedium() const
Definition: TGeoVolume.h:176
Double_t fXY[8][2]
[4] tangents of twist angles
Definition: TGeoArb8.h:29
virtual void SetVertex(Int_t vnum, Double_t x, Double_t y)
Set values for a given vertex.
Definition: TGeoArb8.cxx:1248
Double_t fTl1
Definition: TGeoArb8.h:97
Double_t fAlpha1
Definition: TGeoArb8.h:98
TRAP is a general trapezoid, i.e.
Definition: TGeoArb8.h:89
virtual void DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const
Compute distance from array of input points having directions specified by dirs. Store output in dist...
Definition: TGeoArb8.cxx:1830
virtual void SetPoints(Double_t *points) const
Creates arb8 mesh points.
Definition: TGeoArb8.cxx:1224
virtual void ComputeBBox()
Computes bounding box for an Arb8 shape.
Definition: TGeoArb8.cxx:265
TH1F * h1
Definition: legend1.C:5
virtual Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=0) const
Compute distance from inside point to surface of the shape.
Definition: TGeoArb8.cxx:1952
virtual void Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
Check the inside status for each of the points in the array.
Definition: TGeoArb8.cxx:1288
virtual void ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
Compute the normal for an array o points so that norm.dot.dir is positive Input: Arrays of point coor...
Definition: TGeoArb8.cxx:1298
Double_t fDz
Definition: TGeoArb8.h:27
virtual void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
Compute normal to closest surface from POINT.
Definition: TGeoArb8.cxx:462
const char * GetPointerName() const
Provide a pointer name containing uid.
Definition: TGeoShape.cxx:699
Base finder class for patterns.
point * points
Definition: X3DBuffer.c:20
static Bool_t IsSegCrossing(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Double_t x3, Double_t y3, Double_t x4, Double_t y4)
Check if segments (A,B) and (C,D) are crossing, where: A(x1,y1), B(x2,y2), C(x3,y3), D(x4,y4)
Definition: TGeoShape.cxx:336
virtual const char * GetName() const
Get the shape name.
Definition: TGeoShape.cxx:248
virtual Int_t GetFittingBox(const TGeoBBox *parambox, TGeoMatrix *mat, Double_t &dx, Double_t &dy, Double_t &dz) const
Fills real parameters of a positioned box inside this arb8. Returns 0 if successful.
Definition: TGeoArb8.cxx:873
TGeoArb8()
Default constructor.
Definition: TGeoArb8.cxx:145
virtual Bool_t GetPointsOnFacet(Int_t, Int_t, Double_t *) const
Fills array with n random points located on the surface of indexed facet.
Definition: TGeoArb8.cxx:958
float ymax
Definition: THbookFile.cxx:93
TNamed & operator=(const TNamed &rhs)
TNamed assignment operator.
Definition: TNamed.cxx:51
Base abstract class for all shapes.
Definition: TGeoShape.h:25
TPaveText * pt
virtual ~TGeoTrap()
Destructor.
Definition: TGeoArb8.cxx:1440
virtual Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const
Computes the closest distance from given point to this shape.
Definition: TGeoArb8.cxx:2023
auto * a
Definition: textangle.C:12
static Bool_t InsidePolygon(Double_t x, Double_t y, Double_t *pts)
Finds if a point in XY plane is inside the polygon defines by PTS.
Definition: TGeoArb8.cxx:1014
Double_t fPhi
Definition: TGeoArb8.h:94
Double_t fAlpha2
Definition: TGeoArb8.h:102
virtual void Sizeof3D() const
Definition: TGeoBBox.cxx:996
void ComputeTwist()
Computes tangents of twist angles (angles between projections on XY plane of corresponding -dz +dz ed...
Definition: TGeoArb8.cxx:291
Double_t fTheta
Definition: TGeoArb8.h:93
Double_t fH2
Definition: TGeoArb8.h:99
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
virtual const Double_t * GetOrigin() const
Definition: TGeoBBox.h:73
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save a primitive as a C++ statement(s) on output stream "out".
Definition: TGeoArb8.cxx:1166
static Bool_t IsSamePoint(const Double_t *p1, const Double_t *p2)
Definition: TGeoArb8.h:70
static double p1(double t, double a, double b)
float xmax
Definition: THbookFile.cxx:93
virtual Double_t GetDY() const
Definition: TGeoBBox.h:71
TGeoGtra()
Default ctor.
Definition: TGeoArb8.cxx:1858
void SetDivIndex(Int_t index)
virtual Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=0) const
Compute distance from outside point to surface of the trapezoid.
Definition: TGeoArb8.cxx:1500
constexpr Double_t E()
Base of natural log: .
Definition: TMath.h:97
virtual void LocalToMaster(const Double_t *local, Double_t *master) const
convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse
Definition: TGeoMatrix.cxx:339
Double_t Cos(Double_t)
Definition: TMath.h:640
const Bool_t kFALSE
Definition: RtypesCore.h:88
virtual Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const
Computes the closest distance from given point to this shape.
Definition: TGeoBBox.cxx:855
virtual Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const
Computes the closest distance from given point to this shape.
Definition: TGeoArb8.cxx:1049
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
Bool_t IsRotation() const
Definition: TGeoMatrix.h:68
#define d(i)
Definition: RSha256.hxx:102
virtual void DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const
Compute distance from array of input points having directions specified by dirs. Store output in dist...
Definition: TGeoArb8.cxx:1314
virtual TGeoShape * GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix *mat) const
In case shape has some negative parameters, these have to be computed in order to fit the mother...
Definition: TGeoArb8.cxx:1667
static const double x1[5]
#define ClassImp(name)
Definition: Rtypes.h:359
R__EXTERN TGeoManager * gGeoManager
Definition: TGeoManager.h:562
virtual void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
Compute safe distance from each of the points in the input array.
Definition: TGeoArb8.cxx:1848
Double_t fTwistAngle
Definition: TGeoArb8.h:147
double Double_t
Definition: RtypesCore.h:55
virtual Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=0) const
Compute distance from outside point to surface of the box.
Definition: TGeoBBox.cxx:430
Bool_t TestShapeBit(UInt_t f) const
Definition: TGeoShape.h:163
virtual void Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
Compute safe distance from each of the points in the input array.
Definition: TGeoArb8.cxx:1324
Double_t y[n]
Definition: legend1.C:17
Node containing an offset.
Definition: TGeoNode.h:181
static constexpr double s
An arbitrary trapezoid with less than 8 vertices standing on two parallel planes perpendicular to Z a...
Definition: TGeoArb8.h:17
static Double_t Big()
Definition: TGeoShape.h:88
static void GetPlaneNormal(Double_t *p1, Double_t *p2, Double_t *p3, Double_t *norm)
Computes normal to plane defined by P1, P2 and P3.
Definition: TGeoArb8.cxx:930
Double_t fDY
Definition: TGeoBBox.h:22
void SetShapeBit(UInt_t f, Bool_t set)
Equivalent of TObject::SetBit.
Definition: TGeoShape.cxx:524
virtual ~TGeoArb8()
Destructor.
Definition: TGeoArb8.cxx:240
Double_t DistToPlane(const Double_t *point, const Double_t *dir, Int_t ipl, Bool_t in) const
Computes distance to plane ipl :
Definition: TGeoArb8.cxx:551
Double_t fH1
Definition: TGeoArb8.h:95
virtual Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const
Computes the closest distance from given point to this shape.
Definition: TGeoArb8.cxx:1702
virtual void DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const
Compute distance from array of input points having directions specified by dirs. Store output in dist...
Definition: TGeoArb8.cxx:2104
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:200
Double_t Sin(Double_t)
Definition: TMath.h:636
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
virtual Bool_t Contains(const Double_t *point) const
Test if point is inside this shape.
Definition: TGeoBBox.cxx:340
Double_t fDX
Definition: TGeoBBox.h:21
#define c(i)
Definition: RSha256.hxx:101
Double_t fTl2
Definition: TGeoArb8.h:101
TGeoArb8 & operator=(const TGeoArb8 &)
Assignment operator.
Definition: TGeoArb8.cxx:223
Double_t fBl2
Definition: TGeoArb8.h:100
TGeoTrap()
Default ctor.
Definition: TGeoArb8.cxx:1334
Double_t Sqrt(Double_t x)
Definition: TMath.h:690
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:908
virtual void SetDimensions(Double_t *param)
Set all arb8 params in one step.
Definition: TGeoArb8.cxx:1210
virtual void SetDimensions(Double_t *param)
Set all arb8 params in one step.
Definition: TGeoArb8.cxx:2066
virtual TGeoVolume * Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv, Double_t start, Double_t step)
Divide this trapezoid shape belonging to volume "voldiv" into ndiv volumes called divname...
Definition: TGeoArb8.cxx:1619
Long64_t LocMin(Long64_t n, const T *a)
Return index of array with the minimum element.
Definition: TMath.h:971
virtual TGeoShape * GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix *mat) const
In case shape has some negative parameters, these has to be computed in order to fit the mother...
Definition: TGeoArb8.cxx:1983
const Bool_t kTRUE
Definition: RtypesCore.h:87
virtual void DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const
Compute distance from array of input points having directions specified by dirs. Store output in dist...
Definition: TGeoArb8.cxx:1838
Double_t Tan(Double_t)
Definition: TMath.h:644
virtual Double_t DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=0) const
Compute distance from inside point to surface of the shape.
Definition: TGeoArb8.cxx:1967
char name[80]
Definition: TGX11.cxx:109
virtual void SetDimensions(Double_t *param)
Set all arb8 params in one step.
Definition: TGeoArb8.cxx:1794
Double_t * fTwist
Definition: TGeoArb8.h:28
virtual Double_t DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact=1, Double_t step=TGeoShape::Big(), Double_t *safe=0) const
Compute distance from inside point to surface of the shape.
Definition: TGeoArb8.cxx:700
virtual Double_t GetDZ() const
Definition: TGeoBBox.h:72
Bool_t IsRunTimeShape() const
Definition: TGeoShape.h:139
const char * Data() const
Definition: TString.h:364
virtual void DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t *step) const
Compute distance from array of input points having directions specified by dirs. Store output in dist...
Definition: TGeoArb8.cxx:2112