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