Logo ROOT  

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