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