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