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