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