Logo ROOT   6.18/05
Reference Guide
TGeoXtru.cxx
Go to the documentation of this file.
1// @(#)root/geom:$Id$
2// Author: Mihaela Gheata 24/01/04
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/** \class TGeoXtru
13\ingroup Geometry_classes
14 An extrusion with fixed outline shape in x-y and a sequence
15of z extents (segments). The overall scale of the outline scales
16linearly between z points and the center can have an x-y offset.
17
18Based on the initial implementation of R. Hatcher
19
20### Creation of TGeoXtru shape
21
22A TGeoXtru represents a polygonal extrusion. It is defined by the:
231. 'Blueprint' of the arbitrary polygon representing any Z section. This
24 is an arbitrary polygon (convex or not) defined by the X/Y positions of
25 its vertices.
261. A sequence of Z sections ordered on the Z axis. Each section defines the
27 'actual' parameters of the polygon at a given Z. The sections may be
28 translated with respect to the blueprint and/or scaled. The TGeoXtru
29 segment in between 2 Z sections is a solid represented by the linear
30 extrusion between the 2 polygons. Two consecutive sections may be defined
31 at same Z position.
32
331. `TGeoXtru *xtru = TGeoXtru(Int_t nz);`
34
35 where nz=number of Z planes
36
372. `Double_t x[nvertices]; // array of X positions of blueprint polygon vertices`
38
39 `Double_t y[nvertices]; // array of Y positions of blueprint polygon vertices`
40
413. `xtru->DefinePolygon(nvertices,x,y);`
42
434. `DefineSection(0, z0, x0, y0, scale0); // Z position, offset and scale for first section`
44
45 `DefineSection(1, z1, x1, y1, scale1); // -''- second section`
46 `....`
47 `DefineSection(nz-1, zn, xn, yn, scalen); // parameters for last section`
48
49#### NOTES
50Currently navigation functionality not fully implemented (only Contains()).
51Decomposition in concave polygons not implemented - drawing in solid mode
52within x3d produces incorrect end-faces
53*/
54
55#include "TGeoXtru.h"
56
57#include "Riostream.h"
58#include "TVirtualPad.h"
59#include "TBuffer3D.h"
60#include "TBuffer3DTypes.h"
61#include "TClass.h"
62#include "TMath.h"
63
64#include "TVirtualGeoPainter.h"
65#include "TGeoManager.h"
66#include "TGeoVolume.h"
67#include "TGeoPolygon.h"
68
70
71////////////////////////////////////////////////////////////////////////////////
72/// Constructor.
73
75 fSeg(0), fIz(0), fXc(0), fYc(0), fPoly(0)
76{
77}
78
79////////////////////////////////////////////////////////////////////////////////
80/// Destructor.
81
83{
84 delete [] fXc;
85 delete [] fYc;
86 delete fPoly;
87}
88
89////////////////////////////////////////////////////////////////////////////////
90
92{
93 if (!fThreadSize) ((TGeoXtru*)this)->CreateThreadData(1);
95 return *fThreadData[tid];
96}
97
98////////////////////////////////////////////////////////////////////////////////
99
101{
102 std::lock_guard<std::mutex> guard(fMutex);
103 std::vector<ThreadData_t*>::iterator i = fThreadData.begin();
104 while (i != fThreadData.end())
105 {
106 delete *i;
107 ++i;
108 }
109 fThreadData.clear();
110 fThreadSize = 0;
111}
112
113////////////////////////////////////////////////////////////////////////////////
114/// Create thread data for n threads max.
115
117{
118 std::lock_guard<std::mutex> guard(fMutex);
119 fThreadData.resize(nthreads);
120 fThreadSize = nthreads;
121 for (Int_t tid=0; tid<nthreads; tid++) {
122 if (fThreadData[tid] == 0) {
123 fThreadData[tid] = new ThreadData_t;
124 ThreadData_t &td = *fThreadData[tid];
125 td.fXc = new Double_t [fNvert];
126 td.fYc = new Double_t [fNvert];
127 memcpy(td.fXc, fX, fNvert*sizeof(Double_t));
128 memcpy(td.fYc, fY, fNvert*sizeof(Double_t));
129 td.fPoly = new TGeoPolygon(fNvert);
130 td.fPoly->SetXY(td.fXc, td.fYc); // initialize with current coordinates
131 td.fPoly->FinishPolygon();
132 if (tid == 0 && td.fPoly->IsIllegalCheck()) {
133 Error("DefinePolygon", "Shape %s of type XTRU has an illegal polygon.", GetName());
134 }
135 }
136 }
137}
138
139////////////////////////////////////////////////////////////////////////////////
140/// Set current z-plane.
141
143{
144 GetThreadData().fIz = iz;
145}
146////////////////////////////////////////////////////////////////////////////////
147/// Set current segment.
148
150{
151 GetThreadData().fSeg = iseg;
152}
153
154////////////////////////////////////////////////////////////////////////////////
155/// dummy ctor
156
158 :TGeoBBox(),
159 fNvert(0),
160 fNz(0),
161 fZcurrent(0.),
162 fX(0),
163 fY(0),
164 fZ(0),
165 fScale(0),
166 fX0(0),
167 fY0(0),
168 fThreadData(0),
169 fThreadSize(0)
170{
172}
173
174////////////////////////////////////////////////////////////////////////////////
175/// Default constructor
176
178 :TGeoBBox(0, 0, 0),
179 fNvert(0),
180 fNz(nz),
181 fZcurrent(0.),
182 fX(0),
183 fY(0),
184 fZ(new Double_t[nz]),
185 fScale(new Double_t[nz]),
186 fX0(new Double_t[nz]),
187 fY0(new Double_t[nz]),
188 fThreadData(0),
189 fThreadSize(0)
190{
192 if (nz<2) {
193 Error("ctor", "Cannot create TGeoXtru %s with less than 2 Z planes", GetName());
195 return;
196 }
197}
198
199////////////////////////////////////////////////////////////////////////////////
200/// Default constructor in GEANT3 style
201/// - param[0] = nz // number of z planes
202///
203/// - param[1] = z1 // Z position of first plane
204/// - param[2] = x1 // X position of first plane
205/// - param[3] = y1 // Y position of first plane
206/// - param[4] = scale1 // scale factor for first plane
207/// ...
208/// - param[4*(nz-1]+1] = zn
209/// - param[4*(nz-1)+2] = xn
210/// - param[4*(nz-1)+3] = yn
211/// - param[4*(nz-1)+4] = scalen
212
214 :TGeoBBox(0, 0, 0),
215 fNvert(0),
216 fNz(0),
217 fZcurrent(0.),
218 fX(0),
219 fY(0),
220 fZ(0),
221 fScale(0),
222 fX0(0),
223 fY0(0),
224 fThreadData(0),
225 fThreadSize(0)
226{
228 SetDimensions(param);
229}
230
231////////////////////////////////////////////////////////////////////////////////
232///copy constructor
233
235 TGeoBBox(xt),
236 fNvert(0),
237 fNz(0),
238 fZcurrent(0),
239 fX(0),
240 fY(0),
241 fZ(0),
242 fScale(0),
243 fX0(0),
244 fY0(0),
245 fThreadData(0),
246 fThreadSize(0)
247{
248}
249
250////////////////////////////////////////////////////////////////////////////////
251///assignment operator
252
254{
255 if(this!=&xt) {
257 fNvert=0;
258 fNz=0;
259 fZcurrent=0;
260 fX=0;
261 fY=0;
262 fZ=0;
263 fScale=0;
264 fX0=0;
265 fY0=0;
266 fThreadSize=0;
267 }
268 return *this;
269}
270
271////////////////////////////////////////////////////////////////////////////////
272/// destructor
273
275{
276 if (fX) {delete[] fX; fX = 0;}
277 if (fY) {delete[] fY; fY = 0;}
278 if (fZ) {delete[] fZ; fZ = 0;}
279 if (fScale) {delete[] fScale; fScale = 0;}
280 if (fX0) {delete[] fX0; fX0 = 0;}
281 if (fY0) {delete[] fY0; fY0 = 0;}
283}
284
285////////////////////////////////////////////////////////////////////////////////
286/// Compute capacity [length^3] of this shape.
287
289{
291 Int_t iz;
292 Double_t capacity = 0;
293 Double_t area, dz, sc1, sc2;
294 TGeoXtru *xtru = (TGeoXtru*)this;
295 xtru->SetCurrentVertices(0.,0.,1.);
296 area = td.fPoly->Area();
297 for (iz=0; iz<fNz-1; iz++) {
298 dz = fZ[iz+1]-fZ[iz];
299 if (TGeoShape::IsSameWithinTolerance(dz,0)) continue;
300 sc1 = fScale[iz];
301 sc2 = fScale[iz+1];
302 capacity += (area*dz/3.)*(sc1*sc1+sc1*sc2+sc2*sc2);
303 }
304 return capacity;
305}
306
307////////////////////////////////////////////////////////////////////////////////
308/// compute bounding box of the pcon
309
311{
313 if (!fX || !fZ || !fNvert) {
314 Error("ComputeBBox", "In shape %s polygon not defined", GetName());
316 return;
317 }
318 Double_t zmin = fZ[0];
319 Double_t zmax = fZ[fNz-1];
324 for (Int_t i=0; i<fNz; i++) {
325 SetCurrentVertices(fX0[i], fY0[i], fScale[i]);
326 for (Int_t j=0; j<fNvert; j++) {
327 if (td.fXc[j]<xmin) xmin=td.fXc[j];
328 if (td.fXc[j]>xmax) xmax=td.fXc[j];
329 if (td.fYc[j]<ymin) ymin=td.fYc[j];
330 if (td.fYc[j]>ymax) ymax=td.fYc[j];
331 }
332 }
333 fOrigin[0] = 0.5*(xmin+xmax);
334 fOrigin[1] = 0.5*(ymin+ymax);
335 fOrigin[2] = 0.5*(zmin+zmax);
336 fDX = 0.5*(xmax-xmin);
337 fDY = 0.5*(ymax-ymin);
338 fDZ = 0.5*(zmax-zmin);
339}
340
341////////////////////////////////////////////////////////////////////////////////
342/// Compute normal to closest surface from POINT.
343
344void TGeoXtru::ComputeNormal(const Double_t * /*point*/, const Double_t *dir, Double_t *norm)
345{
347 if (td.fIz<0) {
348 memset(norm,0,3*sizeof(Double_t));
349 norm[2] = (dir[2]>0)?1:-1;
350 return;
351 }
352 Double_t vert[12];
353 GetPlaneVertices(td.fIz, td.fSeg, vert);
354 GetPlaneNormal(vert, norm);
355 Double_t ndotd = norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2];
356 if (ndotd<0) {
357 norm[0] = -norm[0];
358 norm[1] = -norm[1];
359 norm[2] = -norm[2];
360 }
361}
362
363////////////////////////////////////////////////////////////////////////////////
364/// test if point is inside this shape
365
367{
369 // Check Z range
370 TGeoXtru *xtru = (TGeoXtru*)this;
371 if (point[2]<fZ[0]) return kFALSE;
372 if (point[2]>fZ[fNz-1]) return kFALSE;
373 Int_t iz = TMath::BinarySearch(fNz, fZ, point[2]);
374 if (iz<0 || iz==fNz-1) return kFALSE;
375 if (TGeoShape::IsSameWithinTolerance(point[2],fZ[iz])) {
376 xtru->SetIz(-1);
377 xtru->SetCurrentVertices(fX0[iz],fY0[iz], fScale[iz]);
378 if (td.fPoly->Contains(point)) return kTRUE;
379 if (iz>1 && TGeoShape::IsSameWithinTolerance(fZ[iz],fZ[iz-1])) {
380 xtru->SetCurrentVertices(fX0[iz-1],fY0[iz-1], fScale[iz-1]);
381 return td.fPoly->Contains(point);
382 } else if (iz<fNz-2 && TGeoShape::IsSameWithinTolerance(fZ[iz],fZ[iz+1])) {
383 xtru->SetCurrentVertices(fX0[iz+1],fY0[iz+1], fScale[iz+1]);
384 return td.fPoly->Contains(point);
385 }
386 }
387 xtru->SetCurrentZ(point[2], iz);
388 if (TMath::Abs(point[2]-fZ[iz])<TGeoShape::Tolerance() ||
389 TMath::Abs(fZ[iz+1]-point[2])<TGeoShape::Tolerance()) xtru->SetIz(-1);
390 // Now td.fXc,fYc represent the vertices of the section at point[2]
391 return td.fPoly->Contains(point);
392}
393
394////////////////////////////////////////////////////////////////////////////////
395/// compute closest distance from point px,py to each corner
396
398{
399 const Int_t numPoints = fNvert*fNz;
400 return ShapeDistancetoPrimitive(numPoints, px, py);
401}
402
403////////////////////////////////////////////////////////////////////////////////
404/// Draw the section polygon.
405
407 {
409 if (td.fPoly) td.fPoly->Draw(option);
410}
411
412////////////////////////////////////////////////////////////////////////////////
413/// Compute distance to a Xtru lateral surface.
414
415Double_t TGeoXtru::DistToPlane(const Double_t *point, const Double_t *dir, Int_t iz, Int_t ivert, Double_t stepmax, Bool_t in) const
416{
419 Double_t vert[12];
420 Double_t norm[3];
421 Double_t znew;
422 Double_t pt[3];
423 Double_t safe;
424 if (TGeoShape::IsSameWithinTolerance(fZ[iz],fZ[iz+1]) && !in) {
425 TGeoXtru *xtru = (TGeoXtru*)this;
426 snext = (fZ[iz]-point[2])/dir[2];
427 if (snext<0) return TGeoShape::Big();
428 pt[0] = point[0]+snext*dir[0];
429 pt[1] = point[1]+snext*dir[1];
430 pt[2] = point[2]+snext*dir[2];
431 if (dir[2] < 0.) xtru->SetCurrentVertices(fX0[iz], fY0[iz], fScale[iz]);
432 else xtru->SetCurrentVertices(fX0[iz+1], fY0[iz+1], fScale[iz+1]);
433 if (!td.fPoly->Contains(pt)) return TGeoShape::Big();
434 return snext;
435 }
436 GetPlaneVertices(iz, ivert, vert);
437 GetPlaneNormal(vert, norm);
438 Double_t ndotd = norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2];
439 if (in) {
440 if (ndotd<=0) return TGeoShape::Big();
441 safe = (vert[0]-point[0])*norm[0]+
442 (vert[1]-point[1])*norm[1]+
443 (vert[2]-point[2])*norm[2];
444 if (safe<-1.E-8) return TGeoShape::Big(); // direction outwards plane
445 } else {
446 ndotd = -ndotd;
447 if (ndotd<=0) return TGeoShape::Big();
448 safe = (point[0]-vert[0])*norm[0]+
449 (point[1]-vert[1])*norm[1]+
450 (point[2]-vert[2])*norm[2];
451 if (safe<-1.E-8) return TGeoShape::Big(); // direction outwards plane
452 }
453 snext = safe/ndotd;
454 if (snext>stepmax) return TGeoShape::Big();
455 if (fZ[iz]<fZ[iz+1]) {
456 znew = point[2] + snext*dir[2];
457 if (znew<fZ[iz]) return TGeoShape::Big();
458 if (znew>fZ[iz+1]) return TGeoShape::Big();
459 }
460 pt[0] = point[0]+snext*dir[0];
461 pt[1] = point[1]+snext*dir[1];
462 pt[2] = point[2]+snext*dir[2];
463 if (!IsPointInsidePlane(pt, vert, norm)) return TGeoShape::Big();
464 return TMath::Max(snext, 0.);
465}
466
467////////////////////////////////////////////////////////////////////////////////
468/// compute distance from inside point to surface of the polycone
469/// locate Z segment
470
471Double_t TGeoXtru::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
472{
474 if (iact<3 && safe) {
475 *safe = Safety(point, kTRUE);
476 if (iact==0) return TGeoShape::Big();
477 if (iact==1 && step<*safe) return TGeoShape::Big();
478 }
479 TGeoXtru *xtru = (TGeoXtru*)this;
480 Int_t iz = TMath::BinarySearch(fNz, fZ, point[2]);
481 if (iz < 0) {
482 if (dir[2]<=0) {
483 xtru->SetIz(-1);
484 return 0.;
485 }
486 iz = 0;
487 }
488 if (iz==fNz-1) {
489 if (dir[2]>=0) {
490 xtru->SetIz(-1);
491 return 0.;
492 }
493 iz--;
494 } else {
495 if (iz>0) {
496 if (TGeoShape::IsSameWithinTolerance(point[2],fZ[iz])) {
497 if (TGeoShape::IsSameWithinTolerance(fZ[iz],fZ[iz+1]) && dir[2]<0) iz++;
498 else if (TGeoShape::IsSameWithinTolerance(fZ[iz],fZ[iz-1]) && dir[2]>0) iz--;
499 }
500 }
501 }
502 Bool_t convex = td.fPoly->IsConvex();
503// Double_t stepmax = step;
504// if (stepmax>TGeoShape::Big()) stepmax = TGeoShape::Big();
506 Double_t dist, sz;
507 Double_t pt[3];
508 Int_t iv, ipl, inext;
509 // we treat the special case when dir[2]=0
510 if (TGeoShape::IsSameWithinTolerance(dir[2],0)) {
511 for (iv=0; iv<fNvert; iv++) {
512 xtru->SetIz(-1);
513 dist = DistToPlane(point,dir,iz,iv,TGeoShape::Big(),kTRUE);
514 if (dist<snext) {
515 snext = dist;
516 xtru->SetSeg(iv);
517 if (convex) return snext;
518 }
519 }
520 if (snext < 1.E10) return snext;
521 return TGeoShape::Tolerance();
522 }
523
524 // normal case
525 Int_t incseg = (dir[2]>0)?1:-1;
526 Int_t iznext = iz;
527 Bool_t zexit = kFALSE;
528 while (iz>=0 && iz<fNz-1) {
529 // find the distance to current segment end Z surface
530 ipl = iz+((incseg+1)>>1); // next plane
531 inext = ipl+incseg; // next next plane
532 sz = (fZ[ipl]-point[2])/dir[2];
533 if (sz<snext) {
534 iznext += incseg;
535 // we cross the next Z section before stepmax
536 pt[0] = point[0]+sz*dir[0];
537 pt[1] = point[1]+sz*dir[1];
538 xtru->SetCurrentVertices(fX0[ipl],fY0[ipl],fScale[ipl]);
539 if (td.fPoly->Contains(pt)) {
540 // ray gets through next polygon - is it the last one?
541 if (ipl==0 || ipl==fNz-1) {
542 xtru->SetIz(-1);
543 if (convex) return sz;
544 zexit = kTRUE;
545 snext = sz;
546 }
547 // maybe a Z discontinuity - check this
548 if (!zexit && TGeoShape::IsSameWithinTolerance(fZ[ipl],fZ[inext])) {
549 xtru->SetCurrentVertices(fX0[inext],fY0[inext],fScale[inext]);
550 // if we do not cross the next polygone, we are out
551 if (!td.fPoly->Contains(pt)) {
552 xtru->SetIz(-1);
553 if (convex) return sz;
554 zexit = kTRUE;
555 snext = sz;
556 } else {
557 iznext = inext;
558 }
559 }
560 }
561 } else {
562 iznext = fNz-1; // stop
563 }
564 // ray may cross the lateral surfaces of section iz
565 for (iv=0; iv<fNvert; iv++) {
566 dist = DistToPlane(point,dir,iz,iv,TGeoShape::Big(),kTRUE);
567 if (dist<snext) {
568 xtru->SetIz(iz);
569 xtru->SetSeg(iv);
570 snext = dist;
571 if (convex) return snext;
572 zexit = kTRUE;
573 }
574 }
575 if (zexit) return snext;
576 iz = iznext;
577 }
578 return TGeoShape::Tolerance();
579}
580
581////////////////////////////////////////////////////////////////////////////////
582/// compute distance from outside point to surface of the tube
583/// Warning("DistFromOutside", "not implemented");
584
585Double_t TGeoXtru::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
586{
588 if (iact<3 && safe) {
589 *safe = Safety(point, kTRUE);
590 if (iact==0) return TGeoShape::Big();
591 if (iact==1 && step<*safe) return TGeoShape::Big();
592 }
593// Check if the bounding box is crossed within the requested distance
594 Double_t sdist = TGeoBBox::DistFromOutside(point,dir, fDX, fDY, fDZ, fOrigin, step);
595 if (sdist>=step) return TGeoShape::Big();
596 Double_t stepmax = step;
597 if (stepmax>TGeoShape::Big()) stepmax = TGeoShape::Big();
598 Double_t snext = 0.;
600 Int_t i, iv;
601 Double_t pt[3];
602 memcpy(pt,point,3*sizeof(Double_t));
603 TGeoXtru *xtru = (TGeoXtru*)this;
604 // We might get out easy with Z checks
605 Int_t iz = TMath::BinarySearch(fNz, fZ, point[2]);
606 if (iz<0) {
607 if (dir[2]<=0) return TGeoShape::Big();
608 // propagate to first Z plane
609 snext = (fZ[0] - point[2])/dir[2];
610 if (snext>stepmax) return TGeoShape::Big();
611 for (i=0; i<3; i++) pt[i] = point[i] + snext*dir[i];
612 xtru->SetCurrentVertices(fX0[0],fY0[0],fScale[0]);
613 if (td.fPoly->Contains(pt)) {
614 xtru->SetIz(-1);
615 return snext;
616 }
617 iz=0; // valid starting value = first segment
618 stepmax -= snext;
619 } else {
620 if (iz==fNz-1) {
621 if (dir[2]>=0) return TGeoShape::Big();
622 // propagate to last Z plane
623 snext = (fZ[fNz-1] - point[2])/dir[2];
624 if (snext>stepmax) return TGeoShape::Big();
625 for (i=0; i<3; i++) pt[i] = point[i] + snext*dir[i];
626 xtru->SetCurrentVertices(fX0[fNz-1],fY0[fNz-1],fScale[fNz-1]);
627 if (td.fPoly->Contains(pt)) {
628 xtru->SetIz(-1);
629 return snext;
630 }
631 iz = fNz-2; // valid value = last segment
632 stepmax -= snext;
633 }
634 }
635 // Check if the bounding box is missed by the track
636 if (!TGeoBBox::Contains(pt)) {
638 if (dist>stepmax) return TGeoShape::Big();
639 if (dist>1E-6) dist-=1E-6; // decrease snext to make sure we do not cross the xtru
640 else dist = 0;
641 for (i=0; i<3; i++) pt[i] += dist*dir[i]; // we are now closer
642 iz = TMath::BinarySearch(fNz, fZ, pt[2]);
643 if (iz<0) iz=0;
644 else if (iz==fNz-1) iz = fNz-2;
645 snext += dist;
646 stepmax -= dist;
647 }
648 // not the case - we have to do some work...
649 // Start tracking from current iz
650 // - first solve particular case dir[2]=0
651 Bool_t convex = td.fPoly->IsConvex();
652 Bool_t hit = kFALSE;
653 if (TGeoShape::IsSameWithinTolerance(dir[2],0)) {
654 // loop lateral planes to see if we cross something
655 xtru->SetIz(iz);
656 for (iv=0; iv<fNvert; iv++) {
657 dist = DistToPlane(pt,dir,iz,iv,stepmax,kFALSE);
658 if (dist<stepmax) {
659 xtru->SetSeg(iv);
660 if (convex) return (snext+dist);
661 stepmax = dist;
662 hit = kTRUE;
663 }
664 }
665 if (hit) return (snext+stepmax);
666 return TGeoShape::Big();
667 }
668 // general case
669 Int_t incseg = (dir[2]>0)?1:-1;
670 while (iz>=0 && iz<fNz-1) {
671 // compute distance to lateral planes
672 xtru->SetIz(iz);
673 if (TGeoShape::IsSameWithinTolerance(fZ[iz],fZ[iz+1])) xtru->SetIz(-1);
674 for (iv=0; iv<fNvert; iv++) {
675 dist = DistToPlane(pt,dir,iz,iv,stepmax,kFALSE);
676 if (dist<stepmax) {
677 // HIT
678 xtru->SetSeg(iv);
679 if (convex) return (snext+dist);
680 stepmax = dist;
681 hit = kTRUE;
682 }
683 }
684 if (hit) return (snext+stepmax);
685 iz += incseg;
686 }
687 return TGeoShape::Big();
688}
689
690////////////////////////////////////////////////////////////////////////////////
691/// Creates the polygon representing the blueprint of any Xtru section.
692/// - nvert = number of vertices >2
693/// - xv[nvert] = array of X vertex positions
694/// - yv[nvert] = array of Y vertex positions
695///
696/// *NOTE* should be called before DefineSection or ctor with 'param'
697
699{
700 if (nvert<3) {
701 Error("DefinePolygon","In shape %s cannot create polygon with less than 3 vertices", GetName());
703 return kFALSE;
704 }
705 for (Int_t i=0; i<nvert-1; i++) {
706 for (Int_t j=i+1; j<nvert; j++) {
707 if (TMath::Abs(xv[i]-xv[j])<TGeoShape::Tolerance() &&
708 TMath::Abs(yv[i]-yv[j])<TGeoShape::Tolerance()) {
709 Error("DefinePolygon","In shape %s 2 vertices cannot be identical",GetName());
711// return kFALSE;
712 }
713 }
714 }
715 fNvert = nvert;
716 if (fX) delete [] fX;
717 fX = new Double_t[nvert];
718 if (fY) delete [] fY;
719 fY = new Double_t[nvert];
720 memcpy(fX,xv,nvert*sizeof(Double_t));
721 memcpy(fY,yv,nvert*sizeof(Double_t));
722
724
725 return kTRUE;
726}
727
728////////////////////////////////////////////////////////////////////////////////
729/// defines z position of a section plane, rmin and rmax at this z.
730
732{
733 if ((snum<0) || (snum>=fNz)) return;
734 fZ[snum] = z;
735 fX0[snum] = x0;
736 fY0[snum] = y0;
737 fScale[snum] = scale;
738 if (snum) {
739 if (fZ[snum]<fZ[snum-1]) {
740 Warning("DefineSection", "In shape: %s, Z position of section "
741 "%i, z=%e, not in increasing order, %i, z=%e",
742 GetName(),snum,fZ[snum],snum-1,fZ[snum-1]);
743 return;
744 }
745 }
746 if (snum==(fNz-1)) {
747 ComputeBBox();
749 }
750}
751
752////////////////////////////////////////////////////////////////////////////////
753/// Return the Z coordinate for segment ipl.
754
756{
757 if (ipl<0 || ipl>(fNz-1)) {
758 Error("GetZ","In shape %s, ipl=%i out of range (0,%i)",GetName(),ipl,fNz-1);
759 return 0.;
760 }
761 return fZ[ipl];
762}
763////////////////////////////////////////////////////////////////////////////////
764/// Returns normal vector to the planar quadrilateral defined by vector VERT.
765/// The normal points outwards the xtru.
766
767void TGeoXtru::GetPlaneNormal(const Double_t *vert, Double_t *norm) const
768{
769 Double_t cross = 0.;
770 Double_t v1[3], v2[3];
771 v1[0] = vert[9]-vert[0];
772 v1[1] = vert[10]-vert[1];
773 v1[2] = vert[11]-vert[2];
774 v2[0] = vert[3]-vert[0];
775 v2[1] = vert[4]-vert[1];
776 v2[2] = vert[5]-vert[2];
777 norm[0] = v1[1]*v2[2]-v1[2]*v2[1];
778 cross += norm[0]*norm[0];
779 norm[1] = v1[2]*v2[0]-v1[0]*v2[2];
780 cross += norm[1]*norm[1];
781 norm[2] = v1[0]*v2[1]-v1[1]*v2[0];
782 cross += norm[2]*norm[2];
783 if (cross < TGeoShape::Tolerance()) return;
784 cross = 1./TMath::Sqrt(cross);
785 for (Int_t i=0; i<3; i++) norm[i] *= cross;
786}
787
788////////////////////////////////////////////////////////////////////////////////
789/// Returns (x,y,z) of 3 vertices of the surface defined by Z sections (iz, iz+1)
790/// and polygon vertices (ivert, ivert+1). No range check.
791
792void TGeoXtru::GetPlaneVertices(Int_t iz, Int_t ivert, Double_t *vert) const
793{
795 Double_t x,y,z1,z2;
796 Int_t iv1 = (ivert+1)%fNvert;
797 Int_t icrt = 0;
798 z1 = fZ[iz];
799 z2 = fZ[iz+1];
800 if (td.fPoly->IsClockwise()) {
801 x = fX[ivert]*fScale[iz]+fX0[iz];
802 y = fY[ivert]*fScale[iz]+fY0[iz];
803 vert[icrt++] = x;
804 vert[icrt++] = y;
805 vert[icrt++] = z1;
806 x = fX[iv1]*fScale[iz]+fX0[iz];
807 y = fY[iv1]*fScale[iz]+fY0[iz];
808 vert[icrt++] = x;
809 vert[icrt++] = y;
810 vert[icrt++] = z1;
811 x = fX[iv1]*fScale[iz+1]+fX0[iz+1];
812 y = fY[iv1]*fScale[iz+1]+fY0[iz+1];
813 vert[icrt++] = x;
814 vert[icrt++] = y;
815 vert[icrt++] = z2;
816 x = fX[ivert]*fScale[iz+1]+fX0[iz+1];
817 y = fY[ivert]*fScale[iz+1]+fY0[iz+1];
818 vert[icrt++] = x;
819 vert[icrt++] = y;
820 vert[icrt++] = z2;
821 } else {
822 x = fX[iv1]*fScale[iz]+fX0[iz];
823 y = fY[iv1]*fScale[iz]+fY0[iz];
824 vert[icrt++] = x;
825 vert[icrt++] = y;
826 vert[icrt++] = z1;
827 x = fX[ivert]*fScale[iz]+fX0[iz];
828 y = fY[ivert]*fScale[iz]+fY0[iz];
829 vert[icrt++] = x;
830 vert[icrt++] = y;
831 vert[icrt++] = z1;
832 x = fX[ivert]*fScale[iz+1]+fX0[iz+1];
833 y = fY[ivert]*fScale[iz+1]+fY0[iz+1];
834 vert[icrt++] = x;
835 vert[icrt++] = y;
836 vert[icrt++] = z2;
837 x = fX[iv1]*fScale[iz+1]+fX0[iz+1];
838 y = fY[iv1]*fScale[iz+1]+fY0[iz+1];
839 vert[icrt++] = x;
840 vert[icrt++] = y;
841 vert[icrt++] = z2;
842 }
843}
844////////////////////////////////////////////////////////////////////////////////
845/// Check if the quadrilateral defined by VERT contains a coplanar POINT.
846
848{
849 Double_t v1[3], v2[3];
850 Double_t cross;
851 Int_t j,k;
852 for (Int_t i=0; i<4; i++) { // loop vertices
853 j = 3*i;
854 k = 3*((i+1)%4);
855 v1[0] = point[0]-vert[j];
856 v1[1] = point[1]-vert[j+1];
857 v1[2] = point[2]-vert[j+2];
858 v2[0] = vert[k]-vert[j];
859 v2[1] = vert[k+1]-vert[j+1];
860 v2[2] = vert[k+2]-vert[j+2];
861 cross = (v1[1]*v2[2]-v1[2]*v2[1])*norm[0]+
862 (v1[2]*v2[0]-v1[0]*v2[2])*norm[1]+
863 (v1[0]*v2[1]-v1[1]*v2[0])*norm[2];
864 if (cross<0) return kFALSE;
865 }
866 return kTRUE;
867}
868
869////////////////////////////////////////////////////////////////////////////////
870/// Print actual Xtru parameters.
871
873{
874 printf("*** Shape %s: TGeoXtru ***\n", GetName());
875 printf(" Nz = %i\n", fNz);
876 printf(" List of (x,y) of polygon vertices:\n");
877 for (Int_t ivert = 0; ivert<fNvert; ivert++)
878 printf(" x = %11.5f y = %11.5f\n", fX[ivert],fY[ivert]);
879 for (Int_t ipl=0; ipl<fNz; ipl++)
880 printf(" plane %i: z=%11.5f x0=%11.5f y0=%11.5f scale=%11.5f\n", ipl, fZ[ipl], fX0[ipl], fY0[ipl], fScale[ipl]);
881 printf(" Bounding box:\n");
883}
884
885////////////////////////////////////////////////////////////////////////////////
886/// Creates a TBuffer3D describing *this* shape.
887/// Coordinates are in local reference frame.
888
890{
891 Int_t nz = GetNz();
892 Int_t nvert = GetNvert();
893 Int_t nbPnts = nz*nvert;
894 Int_t nbSegs = nvert*(2*nz-1);
895 Int_t nbPols = nvert*(nz-1)+2;
896
898 nbPnts, 3*nbPnts, nbSegs, 3*nbSegs, nbPols, 6*(nbPols-2)+2*(2+nvert));
899 if (buff)
900 {
901 SetPoints(buff->fPnts);
902 SetSegsAndPols(*buff);
903 }
904
905 return buff;
906}
907
908////////////////////////////////////////////////////////////////////////////////
909/// Fill TBuffer3D structure for segments and polygons.
910
912{
913 Int_t nz = GetNz();
914 Int_t nvert = GetNvert();
916
917 Int_t i,j;
918 Int_t indx, indx2, k;
919 indx = indx2 = 0;
920 for (i=0; i<nz; i++) {
921 // loop Z planes
922 indx2 = i*nvert;
923 // loop polygon segments
924 for (j=0; j<nvert; j++) {
925 k = (j+1)%nvert;
926 buff.fSegs[indx++] = c;
927 buff.fSegs[indx++] = indx2+j;
928 buff.fSegs[indx++] = indx2+k;
929 }
930 } // total: nz*nvert polygon segments
931 for (i=0; i<nz-1; i++) {
932 // loop Z planes
933 indx2 = i*nvert;
934 // loop polygon segments
935 for (j=0; j<nvert; j++) {
936 k = j + nvert;
937 buff.fSegs[indx++] = c;
938 buff.fSegs[indx++] = indx2+j;
939 buff.fSegs[indx++] = indx2+k;
940 }
941 } // total (nz-1)*nvert lateral segments
942
943 indx = 0;
944
945 // fill lateral polygons
946 for (i=0; i<nz-1; i++) {
947 indx2 = i*nvert;
948 for (j=0; j<nvert; j++) {
949 k = (j+1)%nvert;
950 buff.fPols[indx++] = c+j%3;
951 buff.fPols[indx++] = 4;
952 buff.fPols[indx++] = indx2+j;
953 buff.fPols[indx++] = nz*nvert+indx2+k;
954 buff.fPols[indx++] = indx2+nvert+j;
955 buff.fPols[indx++] = nz*nvert+indx2+j;
956 }
957 } // total (nz-1)*nvert polys
958 buff.fPols[indx++] = c+2;
959 buff.fPols[indx++] = nvert;
960 indx2 = 0;
961 for (j = nvert - 1; j >= 0; --j) {
962 buff.fPols[indx++] = indx2+j;
963 }
964
965 buff.fPols[indx++] = c;
966 buff.fPols[indx++] = nvert;
967 indx2 = (nz-1)*nvert;
968
969 for (j=0; j<nvert; j++) {
970 buff.fPols[indx++] = indx2+j;
971 }
972}
973
974////////////////////////////////////////////////////////////////////////////////
975/// Compute safety to sector iz, returning also the closest segment index.
976
978{
980 Double_t safz = TGeoShape::Big();
981 Double_t saf1, saf2;
982 Bool_t in1, in2;
983 Int_t iseg;
984 Double_t safe = TGeoShape::Big();
985 // segment-break case
986 if (TGeoShape::IsSameWithinTolerance(fZ[iz],fZ[iz+1])) {
987 safz = TMath::Abs(point[2]-fZ[iz]);
988 if (safz>safmin) return TGeoShape::Big();
989 SetCurrentVertices(fX0[iz], fY0[iz], fScale[iz]);
990 saf1 = td.fPoly->Safety(point, iseg);
991 in1 = td.fPoly->Contains(point);
992// if (!in1 && saf1>safmin) return TGeoShape::Big();
993 SetCurrentVertices(fX0[iz+1], fY0[iz+1], fScale[iz+1]);
994 saf2 = td.fPoly->Safety(point, iseg);
995 in2 = td.fPoly->Contains(point);
996 if ((in1&!in2)|(in2&!in1)) {
997 safe = safz;
998 } else {
999 safe = TMath::Min(saf1,saf2);
1000 safe = TMath::Max(safe, safz);
1001 }
1002 if (safe>safmin) return TGeoShape::Big();
1003 return safe;
1004 }
1005 // normal case
1006 safz = fZ[iz]-point[2];
1007 if (safz>safmin) return TGeoShape::Big();
1008 if (safz<0) {
1009 saf1 = point[2]-fZ[iz+1];
1010 if (saf1>safmin) return TGeoShape::Big();
1011 if (saf1<0) {
1012 safz = TMath::Max(safz, saf1); // we are in between the 2 Z segments - we ignore safz
1013 } else {
1014 safz = saf1;
1015 }
1016 }
1017
1018 // loop segments
1019 Bool_t found = kFALSE;
1020 Double_t vert[12];
1021 Double_t norm[3];
1022// printf("plane %d: safz=%f in=%d\n", iz, safz, in);
1023 for (iseg=0; iseg<fNvert; iseg++) {
1024 GetPlaneVertices(iz,iseg,vert);
1025 GetPlaneNormal(vert, norm);
1026 saf1 = (point[0]-vert[0])*norm[0]+(point[1]-vert[1])*norm[1]+(point[2]-vert[2])*norm[2];
1027 if (in) saf1 = -saf1;
1028// printf("segment %d: (%f,%f)-(%f,%f) norm=(%f,%f,%f): saf1=%f\n", iseg, vert[0],vert[1],vert[3],vert[4],norm[0],norm[1],norm[2],saf1);
1029 if (saf1<-1.E-8) continue;
1030 safe = TMath::Max(safz, saf1);
1031 safe = TMath::Abs(safe);
1032 if (safe>safmin) continue;
1033 safmin = safe;
1034 found = kTRUE;
1035 }
1036 if (found) return safmin;
1037 return TGeoShape::Big();
1038}
1039
1040////////////////////////////////////////////////////////////////////////////////
1041/// computes the closest distance from given point to this shape, according
1042/// to option. The matching point on the shape is stored in spoint.
1043///> localize the Z segment
1044
1046{
1047 Double_t safmin = TGeoShape::Big();
1048 Double_t safe;
1049 Double_t safz = TGeoShape::Big();
1050 TGeoXtru *xtru = (TGeoXtru*)this;
1051 Int_t iz;
1052 if (in) {
1053 safmin = TMath::Min(point[2]-fZ[0], fZ[fNz-1]-point[2]);
1054 for (iz=0; iz<fNz-1; iz++) {
1055 safe = xtru->SafetyToSector(point, iz, safmin, in);
1056 if (safe<safmin) safmin = safe;
1057 }
1058 return safmin;
1059 }
1060 // Accurate safety is expensive, use the bounding box
1061 if (!TGeoBBox::Contains(point)) return TGeoBBox::Safety(point,in);
1062 iz = TMath::BinarySearch(fNz, fZ, point[2]);
1063 if (iz<0) {
1064 iz = 0;
1065 safz = fZ[0] - point[2];
1066 } else {
1067 if (iz==fNz-1) {
1068 iz = fNz-2;
1069 safz = point[2] - fZ[fNz-1];
1070 }
1071 }
1072 // loop segments from iz up
1073 Int_t i;
1074 for (i=iz; i<fNz-1; i++) {
1075 safe = xtru->SafetyToSector(point,i,safmin, in);
1076 if (safe<safmin) safmin=safe;
1077 }
1078 // loop segments from iz-1 down
1079 for (i=iz-1; i>=0; i--) {
1080 safe = xtru->SafetyToSector(point,i,safmin, in);
1081 if (safe<safmin) safmin=safe;
1082 }
1083 safe = TMath::Min(safmin, safz);
1084 return safe;
1085}
1086
1087////////////////////////////////////////////////////////////////////////////////
1088/// Save a primitive as a C++ statement(s) on output stream "out".
1089
1090void TGeoXtru::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
1091{
1093 out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
1094 out << " nz = " << fNz << ";" << std::endl;
1095 out << " nvert = " << fNvert << ";" << std::endl;
1096 out << " TGeoXtru *xtru = new TGeoXtru(nz);" << std::endl;
1097 out << " xtru->SetName(\"" << GetName() << "\");" << std::endl;
1098 Int_t i;
1099 for (i=0; i<fNvert; i++) {
1100 out << " xvert[" << i << "] = " << fX[i] << "; yvert[" << i << "] = " << fY[i] << ";" << std::endl;
1101 }
1102 out << " xtru->DefinePolygon(nvert,xvert,yvert);" << std::endl;
1103 for (i=0; i<fNz; i++) {
1104 out << " zsect = " << fZ[i] << ";" << std::endl;
1105 out << " x0 = " << fX0[i] << ";" << std::endl;
1106 out << " y0 = " << fY0[i] << ";" << std::endl;
1107 out << " scale0 = " << fScale[i] << ";" << std::endl;
1108 out << " xtru->DefineSection(" << i << ",zsect,x0,y0,scale0);" << std::endl;
1109 }
1110 out << " TGeoShape *" << GetPointerName() << " = xtru;" << std::endl;
1112}
1113
1114////////////////////////////////////////////////////////////////////////////////
1115/// Recompute current section vertices for a given Z position within range of section iz.
1116
1118{
1119 Double_t x0, y0, scale, a, b;
1120 Int_t ind1, ind2;
1121 ind1 = iz;
1122 ind2 = iz+1;
1123 Double_t invdz = 1./(fZ[ind2]-fZ[ind1]);
1124 a = (fX0[ind1]*fZ[ind2]-fX0[ind2]*fZ[ind1])*invdz;
1125 b = (fX0[ind2]-fX0[ind1])*invdz;
1126 x0 = a+b*z;
1127 a = (fY0[ind1]*fZ[ind2]-fY0[ind2]*fZ[ind1])*invdz;
1128 b = (fY0[ind2]-fY0[ind1])*invdz;
1129 y0 = a+b*z;
1130 a = (fScale[ind1]*fZ[ind2]-fScale[ind2]*fZ[ind1])*invdz;
1131 b = (fScale[ind2]-fScale[ind1])*invdz;
1132 scale = a+b*z;
1133 SetCurrentVertices(x0,y0,scale);
1134}
1135
1136////////////////////////////////////////////////////////////////////////////////
1137/// Set current vertex coordinates according X0, Y0 and SCALE.
1138
1140{
1142 for (Int_t i=0; i<fNvert; i++) {
1143 td.fXc[i] = scale*fX[i] + x0;
1144 td.fYc[i] = scale*fY[i] + y0;
1145 }
1146}
1147
1148////////////////////////////////////////////////////////////////////////////////
1149/// - param[0] = nz // number of z planes
1150///
1151/// - param[1] = z1 // Z position of first plane
1152/// - param[2] = x1 // X position of first plane
1153/// - param[3] = y1 // Y position of first plane
1154/// - param[4] = scale1 // scale factor for first plane
1155/// ...
1156/// - param[4*(nz-1]+1] = zn
1157/// - param[4*(nz-1)+2] = xn
1158/// - param[4*(nz-1)+3] = yn
1159/// - param[4*(nz-1)+4] = scalen
1160
1162{
1163 fNz = (Int_t)param[0];
1164 if (fNz<2) {
1165 Error("SetDimensions","Cannot create TGeoXtru %s with less than 2 Z planes",GetName());
1167 return;
1168 }
1169 if (fZ) delete [] fZ;
1170 if (fScale) delete [] fScale;
1171 if (fX0) delete [] fX0;
1172 if (fY0) delete [] fY0;
1173 fZ = new Double_t[fNz];
1174 fScale = new Double_t[fNz];
1175 fX0 = new Double_t[fNz];
1176 fY0 = new Double_t[fNz];
1177
1178 for (Int_t i=0; i<fNz; i++)
1179 DefineSection(i, param[1+4*i], param[2+4*i], param[3+4*i], param[4+4*i]);
1180}
1181
1182////////////////////////////////////////////////////////////////////////////////
1183/// create polycone mesh points
1184
1186{
1188 Int_t i, j;
1189 Int_t indx = 0;
1190 TGeoXtru *xtru = (TGeoXtru*)this;
1191 if (points) {
1192 for (i = 0; i < fNz; i++) {
1193 xtru->SetCurrentVertices(fX0[i], fY0[i], fScale[i]);
1194 if (td.fPoly->IsClockwise()) {
1195 for (j = 0; j < fNvert; j++) {
1196 points[indx++] = td.fXc[j];
1197 points[indx++] = td.fYc[j];
1198 points[indx++] = fZ[i];
1199 }
1200 } else {
1201 for (j = 0; j < fNvert; j++) {
1202 points[indx++] = td.fXc[fNvert-1-j];
1203 points[indx++] = td.fYc[fNvert-1-j];
1204 points[indx++] = fZ[i];
1205 }
1206 }
1207 }
1208 }
1209}
1210
1211////////////////////////////////////////////////////////////////////////////////
1212/// create polycone mesh points
1213
1215{
1217 Int_t i, j;
1218 Int_t indx = 0;
1219 TGeoXtru *xtru = (TGeoXtru*)this;
1220 if (points) {
1221 for (i = 0; i < fNz; i++) {
1222 xtru->SetCurrentVertices(fX0[i], fY0[i], fScale[i]);
1223 if (td.fPoly->IsClockwise()) {
1224 for (j = 0; j < fNvert; j++) {
1225 points[indx++] = td.fXc[j];
1226 points[indx++] = td.fYc[j];
1227 points[indx++] = fZ[i];
1228 }
1229 } else {
1230 for (j = 0; j < fNvert; j++) {
1231 points[indx++] = td.fXc[fNvert-1-j];
1232 points[indx++] = td.fYc[fNvert-1-j];
1233 points[indx++] = fZ[i];
1234 }
1235 }
1236 }
1237 }
1238}
1239
1240////////////////////////////////////////////////////////////////////////////////
1241/// Returns numbers of vertices, segments and polygons composing the shape mesh.
1242
1243void TGeoXtru::GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const
1244{
1245 Int_t nz = GetNz();
1246 Int_t nv = GetNvert();
1247 nvert = nz*nv;
1248 nsegs = nv*(2*nz-1);
1249 npols = nv*(nz-1)+2;
1250}
1251
1252////////////////////////////////////////////////////////////////////////////////
1253/// Return number of vertices of the mesh representation
1254
1256{
1257 Int_t numPoints = fNz*fNvert;
1258 return numPoints;
1259}
1260
1261////////////////////////////////////////////////////////////////////////////////
1262/// fill size of this 3-D object
1263
1265{
1266}
1267
1268////////////////////////////////////////////////////////////////////////////////
1269/// Fills a static 3D buffer and returns a reference.
1270
1271const TBuffer3D & TGeoXtru::GetBuffer3D(Int_t reqSections, Bool_t localFrame) const
1272{
1273 static TBuffer3D buffer(TBuffer3DTypes::kGeneric);
1274
1275 TGeoBBox::FillBuffer3D(buffer, reqSections, localFrame);
1276
1277 if (reqSections & TBuffer3D::kRawSizes) {
1278 Int_t nz = GetNz();
1279 Int_t nvert = GetNvert();
1280 Int_t nbPnts = nz*nvert;
1281 Int_t nbSegs = nvert*(2*nz-1);
1282 Int_t nbPols = nvert*(nz-1)+2;
1283 if (buffer.SetRawSizes(nbPnts, 3*nbPnts, nbSegs, 3*nbSegs, nbPols, 6*(nbPols-2)+2*(2+nvert))) {
1285 }
1286 }
1287 // TODO: Push down to TGeoShape?
1288 if ((reqSections & TBuffer3D::kRaw) && buffer.SectionsValid(TBuffer3D::kRawSizes)) {
1289 SetPoints(buffer.fPnts);
1290 if (!buffer.fLocalFrame) {
1291 TransformPoints(buffer.fPnts, buffer.NbPnts());
1292 }
1293
1294 SetSegsAndPols(buffer);
1296 }
1297
1298 return buffer;
1299}
1300
1301////////////////////////////////////////////////////////////////////////////////
1302/// Check the inside status for each of the points in the array.
1303/// Input: Array of point coordinates + vector size
1304/// Output: Array of Booleans for the inside of each point
1305
1306void TGeoXtru::Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
1307{
1308 for (Int_t i=0; i<vecsize; i++) inside[i] = Contains(&points[3*i]);
1309}
1310
1311////////////////////////////////////////////////////////////////////////////////
1312/// Compute the normal for an array o points so that norm.dot.dir is positive
1313/// Input: Arrays of point coordinates and directions + vector size
1314/// Output: Array of normal directions
1315
1316void TGeoXtru::ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
1317{
1318 for (Int_t i=0; i<vecsize; i++) ComputeNormal(&points[3*i], &dirs[3*i], &norms[3*i]);
1319}
1320
1321////////////////////////////////////////////////////////////////////////////////
1322/// Compute distance from array of input points having directions specified by dirs. Store output in dists
1323
1324void TGeoXtru::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
1325{
1326 for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromInside(&points[3*i], &dirs[3*i], 3, step[i]);
1327}
1328
1329////////////////////////////////////////////////////////////////////////////////
1330/// Compute distance from array of input points having directions specified by dirs. Store output in dists
1331
1332void TGeoXtru::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
1333{
1334 for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromOutside(&points[3*i], &dirs[3*i], 3, step[i]);
1335}
1336
1337////////////////////////////////////////////////////////////////////////////////
1338/// Compute safe distance from each of the points in the input array.
1339/// Input: Array of point coordinates, array of statuses for these points, size of the arrays
1340/// Output: Safety values
1341
1342void TGeoXtru::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
1343{
1344 for (Int_t i=0; i<vecsize; i++) safe[i] = Safety(&points[3*i], inside[i]);
1345}
#define b(i)
Definition: RSha256.hxx:100
#define c(i)
Definition: RSha256.hxx:101
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
float Float_t
Definition: RtypesCore.h:53
const Bool_t kTRUE
Definition: RtypesCore.h:87
const char Option_t
Definition: RtypesCore.h:62
#define ClassImp(name)
Definition: Rtypes.h:365
float xmin
Definition: THbookFile.cxx:93
float ymin
Definition: THbookFile.cxx:93
float xmax
Definition: THbookFile.cxx:93
float ymax
Definition: THbookFile.cxx:93
point * points
Definition: X3DBuffer.c:22
Generic 3D primitive description class.
Definition: TBuffer3D.h:18
Int_t * fPols
Definition: TBuffer3D.h:114
UInt_t NbPnts() const
Definition: TBuffer3D.h:80
Bool_t SectionsValid(UInt_t mask) const
Definition: TBuffer3D.h:67
@ kRawSizes
Definition: TBuffer3D.h:53
void SetSectionsValid(UInt_t mask)
Definition: TBuffer3D.h:65
Int_t * fSegs
Definition: TBuffer3D.h:113
Bool_t fLocalFrame
Definition: TBuffer3D.h:90
Bool_t SetRawSizes(UInt_t reqPnts, UInt_t reqPntsCapacity, UInt_t reqSegs, UInt_t reqSegsCapacity, UInt_t reqPols, UInt_t reqPolsCapacity)
Set kRaw tessellation section of buffer with supplied sizes.
Definition: TBuffer3D.cxx:359
Double_t * fPnts
Definition: TBuffer3D.h:112
Box class.
Definition: TGeoBBox.h:18
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:855
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:430
virtual void InspectShape() const
Prints shape parameters.
Definition: TGeoBBox.cxx:793
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:340
Double_t fDZ
Definition: TGeoBBox.h:23
virtual void FillBuffer3D(TBuffer3D &buffer, Int_t reqSections, Bool_t localFrame) const
Fills the supplied buffer, with sections in desired frame See TBuffer3D.h for explanation of sections...
Definition: TGeoBBox.cxx:1033
static Int_t ThreadId()
Translates the current thread id to an ordinal number.
An arbitrary polygon defined by vertices.
Definition: TGeoPolygon.h:20
Bool_t IsConvex() const
Definition: TGeoPolygon.h:59
Double_t Safety(const Double_t *point, Int_t &isegment) const
Compute minimum distance from POINT to any segment. Returns segment index.
void SetXY(Double_t *x, Double_t *y)
Set X/Y array pointer for the polygon and daughters.
Double_t Area() const
Computes area of the polygon in [length^2].
Bool_t Contains(const Double_t *point) const
Check if a point given by X = point[0], Y = point[1] is inside the polygon.
Bool_t IsIllegalCheck() const
Check for illegal crossings between non-consecutive segments.
virtual void Draw(Option_t *option="")
Draw the polygon.
Bool_t IsClockwise() const
Definition: TGeoPolygon.h:58
void FinishPolygon()
Decompose polygon in a convex outscribed part and a list of daughter polygons that have to be subtrac...
static Double_t Big()
Definition: TGeoShape.h:88
Int_t GetBasicColor() const
Get the basic color (0-7).
Definition: TGeoShape.cxx:673
void TransformPoints(Double_t *points, UInt_t NbPoints) const
Tranform a set of points (LocalToMaster)
Definition: TGeoShape.cxx:552
void SetShapeBit(UInt_t f, Bool_t set)
Equivalent of TObject::SetBit.
Definition: TGeoShape.cxx:524
static Bool_t IsSameWithinTolerance(Double_t a, Double_t b)
Check if two numbers differ with less than a tolerance.
Definition: TGeoShape.cxx:326
const char * GetPointerName() const
Provide a pointer name containing uid.
Definition: TGeoShape.cxx:699
Int_t ShapeDistancetoPrimitive(Int_t numpoints, Int_t px, Int_t py) const
Returns distance to shape primitive mesh.
Definition: TGeoShape.cxx:259
virtual const char * GetName() const
Get the shape name.
Definition: TGeoShape.cxx:248
@ kGeoSavePrimitive
Definition: TGeoShape.h:65
@ kGeoXtru
Definition: TGeoShape.h:61
static Double_t Tolerance()
Definition: TGeoShape.h:91
Bool_t TestShapeBit(UInt_t f) const
Definition: TGeoShape.h:163
An extrusion with fixed outline shape in x-y and a sequence of z extents (segments).
Definition: TGeoXtru.h:22
Double_t * fScale
Definition: TGeoXtru.h:47
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 tube Warning("DistFromOutside",...
Definition: TGeoXtru.cxx:585
Double_t * fX0
Definition: TGeoXtru.h:48
TGeoXtru()
dummy ctor
Definition: TGeoXtru.cxx:157
virtual void ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
Compute normal to closest surface from POINT.
Definition: TGeoXtru.cxx:344
virtual Bool_t Contains(const Double_t *point) const
test if point is inside this shape
Definition: TGeoXtru.cxx:366
ThreadData_t & GetThreadData() const
Definition: TGeoXtru.cxx:91
Int_t fNz
Definition: TGeoXtru.h:42
virtual void ComputeBBox()
compute bounding box of the pcon
Definition: TGeoXtru.cxx:310
Double_t * fZ
Definition: TGeoXtru.h:46
Int_t fThreadSize
Navigation data per thread.
Definition: TGeoXtru.h:52
Int_t fNvert
Definition: TGeoXtru.h:41
void SetCurrentVertices(Double_t x0, Double_t y0, Double_t scale)
Set current vertex coordinates according X0, Y0 and SCALE.
Definition: TGeoXtru.cxx:1139
virtual Double_t Capacity() const
Compute capacity [length^3] of this shape.
Definition: TGeoXtru.cxx:288
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: TGeoXtru.cxx:1316
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: TGeoXtru.cxx:1332
void GetPlaneVertices(Int_t iz, Int_t ivert, Double_t *vert) const
Returns (x,y,z) of 3 vertices of the surface defined by Z sections (iz, iz+1) and polygon vertices (i...
Definition: TGeoXtru.cxx:792
virtual void Sizeof3D() const
fill size of this 3-D object
Definition: TGeoXtru.cxx:1264
Double_t * GetZ() const
Definition: TGeoXtru.h:100
Bool_t IsPointInsidePlane(const Double_t *point, Double_t *vert, Double_t *norm) const
Check if the quadrilateral defined by VERT contains a coplanar POINT.
Definition: TGeoXtru.cxx:847
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save a primitive as a C++ statement(s) on output stream "out".
Definition: TGeoXtru.cxx:1090
TGeoXtru & operator=(const TGeoXtru &)
assignment operator
Definition: TGeoXtru.cxx:253
Int_t GetNvert() const
Definition: TGeoXtru.h:94
Double_t * fY0
Definition: TGeoXtru.h:49
virtual void SetPoints(Double_t *points) const
create polycone mesh points
Definition: TGeoXtru.cxx:1185
Bool_t DefinePolygon(Int_t nvert, const Double_t *xv, const Double_t *yv)
Creates the polygon representing the blueprint of any Xtru section.
Definition: TGeoXtru.cxx:698
virtual void CreateThreadData(Int_t nthreads)
Create thread data for n threads max.
Definition: TGeoXtru.cxx:116
void SetCurrentZ(Double_t z, Int_t iz)
Recompute current section vertices for a given Z position within range of section iz.
Definition: TGeoXtru.cxx:1117
virtual TBuffer3D * MakeBuffer3D() const
Creates a TBuffer3D describing this shape.
Definition: TGeoXtru.cxx:889
virtual void DefineSection(Int_t snum, Double_t z, Double_t x0=0., Double_t y0=0., Double_t scale=1.)
defines z position of a section plane, rmin and rmax at this z.
Definition: TGeoXtru.cxx:731
void DrawPolygon(Option_t *option="")
Draw the section polygon.
Definition: TGeoXtru.cxx:406
virtual const TBuffer3D & GetBuffer3D(Int_t reqSections, Bool_t localFrame) const
Fills a static 3D buffer and returns a reference.
Definition: TGeoXtru.cxx:1271
Double_t SafetyToSector(const Double_t *point, Int_t iz, Double_t safmin, Bool_t in)
Compute safety to sector iz, returning also the closest segment index.
Definition: TGeoXtru.cxx:977
Double_t * fY
Definition: TGeoXtru.h:45
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: TGeoXtru.cxx:1342
Double_t * fX
Definition: TGeoXtru.h:44
virtual void SetSegsAndPols(TBuffer3D &buff) const
Fill TBuffer3D structure for segments and polygons.
Definition: TGeoXtru.cxx:911
virtual void GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const
Returns numbers of vertices, segments and polygons composing the shape mesh.
Definition: TGeoXtru.cxx:1243
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: TGeoXtru.cxx:1324
virtual void SetDimensions(Double_t *param)
Definition: TGeoXtru.cxx:1161
virtual void ClearThreadData() const
Definition: TGeoXtru.cxx:100
virtual ~TGeoXtru()
destructor
Definition: TGeoXtru.cxx:274
Double_t DistToPlane(const Double_t *point, const Double_t *dir, Int_t iz, Int_t ivert, Double_t stepmax, Bool_t in) const
Compute distance to a Xtru lateral surface.
Definition: TGeoXtru.cxx:415
virtual Int_t DistancetoPrimitive(Int_t px, Int_t py)
compute closest distance from point px,py to each corner
Definition: TGeoXtru.cxx:397
void SetSeg(Int_t iseg)
Set current segment.
Definition: TGeoXtru.cxx:149
std::mutex fMutex
size of thread-specific array
Definition: TGeoXtru.h:53
virtual void InspectShape() const
Print actual Xtru parameters.
Definition: TGeoXtru.cxx:872
virtual Double_t Safety(const Double_t *point, Bool_t in=kTRUE) const
computes the closest distance from given point to this shape, according to option.
Definition: TGeoXtru.cxx:1045
void GetPlaneNormal(const Double_t *vert, Double_t *norm) const
Returns normal vector to the planar quadrilateral defined by vector VERT.
Definition: TGeoXtru.cxx:767
virtual Int_t GetNmeshVertices() const
Return number of vertices of the mesh representation.
Definition: TGeoXtru.cxx:1255
std::vector< ThreadData_t * > fThreadData
Definition: TGeoXtru.h:51
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: TGeoXtru.cxx:1306
void SetIz(Int_t iz)
Set current z-plane.
Definition: TGeoXtru.cxx:142
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 polycone locate Z segment
Definition: TGeoXtru.cxx:471
Int_t GetNz() const
Definition: TGeoXtru.h:93
Double_t fZcurrent
Definition: TGeoXtru.h:43
TNamed & operator=(const TNamed &rhs)
TNamed assignment operator.
Definition: TNamed.cxx:51
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:172
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:128
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:866
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:694
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
TPaveText * pt
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:212
constexpr Double_t E()
Base of natural log:
Definition: TMath.h:97
Double_t Sqrt(Double_t x)
Definition: TMath.h:679
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:180
Long64_t BinarySearch(Long64_t n, const T *array, T value)
Definition: TMathBase.h:278
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
~ThreadData_t()
Destructor.
Definition: TGeoXtru.cxx:82
TGeoPolygon * fPoly
Definition: TGeoXtru.h:30
ThreadData_t()
Constructor.
Definition: TGeoXtru.cxx:74
auto * a
Definition: textangle.C:12
#define snext(osub1, osub2)
Definition: triangle.c:1167