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