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