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