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