Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TGeoHype.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Mihaela Gheata 20/11/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 #include "Riostream.h"
14 
15 #include "TGeoManager.h"
16 #include "TGeoVolume.h"
17 #include "TVirtualGeoPainter.h"
18 #include "TGeoHype.h"
19 #include "TVirtualPad.h"
20 #include "TBuffer3D.h"
21 #include "TBuffer3DTypes.h"
22 #include "TMath.h"
23 
24 /** \class TGeoHype
25 \ingroup Geometry_classes
26 
27 Hyperboloid class defined by 5 parameters. Bounded by:
28  - Two z planes at z=+/-dz
29  - Inner and outer lateral surfaces. These represent the surfaces
30  described by the revolution of 2 hyperbolas about the Z axis:
31  r^2 - (t*z)^2 = a^2
32 
33  - r = distance between hyperbola and Z axis at coordinate z
34  - t = tangent of the stereo angle (angle made by hyperbola
35  asymptotic lines and Z axis). t=0 means cylindrical surface.
36  - a = distance between hyperbola and Z axis at z=0
37 
38 The inner hyperbolic surface is described by:
39  r^2 - (tin*z)^2 = rin^2
40  - absence of the inner surface (filled hyperboloid can be forced
41  by rin=0 and sin=0
42 The outer hyperbolic surface is described by:
43  r^2 - (tout*z)^2 = rout^2
44 TGeoHype parameters: dz[cm], rin[cm], sin[deg], rout[cm], sout[deg].
45 MANDATORY conditions:
46 
47  - rin < rout
48  - rout > 0
49  - rin^2 + (tin*dz)^2 > rout^2 + (tout*dz)^2
50 
51 SUPPORTED CASES:
52 
53  - rin = 0, tin != 0 => inner surface conical
54  - tin=0 AND/OR tout=0 => corresponding surface(s) cylindrical
55  e.g. tin=0 AND tout=0 => shape becomes a tube with: rmin,rmax,dz
56 */
57 
58 ClassImp(TGeoHype);
59 
60 ////////////////////////////////////////////////////////////////////////////////
61 /// Default constructor
62 
63 TGeoHype::TGeoHype()
64 {
65  SetShapeBit(TGeoShape::kGeoHype);
66  fStIn = 0.;
67  fStOut = 0.;
68  fTin = 0.;
69  fTinsq = 0.;
70  fTout = 0.;
71  fToutsq = 0.;
72 }
73 
74 
75 ////////////////////////////////////////////////////////////////////////////////
76 /// Constructor specifying hyperboloid parameters.
77 
78 TGeoHype::TGeoHype(Double_t rin, Double_t stin, Double_t rout, Double_t stout, Double_t dz)
79  :TGeoTube(rin, rout, dz)
80 {
81  SetShapeBit(TGeoShape::kGeoHype);
82  SetHypeDimensions(rin, stin, rout, stout, dz);
83  // dz<0 can be used to force dz of hyperboloid fit the container volume
84  if (fDz<0) SetShapeBit(kGeoRunTimeShape);
85  ComputeBBox();
86 }
87 ////////////////////////////////////////////////////////////////////////////////
88 /// Constructor specifying parameters and name.
89 
90 TGeoHype::TGeoHype(const char *name,Double_t rin, Double_t stin, Double_t rout, Double_t stout, Double_t dz)
91  :TGeoTube(name, rin, rout, dz)
92 {
93  SetShapeBit(TGeoShape::kGeoHype);
94  SetHypeDimensions(rin, stin, rout, stout, dz);
95  // dz<0 can be used to force dz of hyperboloid fit the container volume
96  if (fDz<0) SetShapeBit(kGeoRunTimeShape);
97  ComputeBBox();
98 }
99 
100 ////////////////////////////////////////////////////////////////////////////////
101 /// Default constructor specifying a list of parameters
102 /// - param[0] = dz
103 /// - param[1] = rin
104 /// - param[2] = stin
105 /// - param[3] = rout
106 /// - param[4] = stout
107 
108 TGeoHype::TGeoHype(Double_t *param)
109  :TGeoTube(param[1],param[3],param[0])
110 {
111  SetShapeBit(TGeoShape::kGeoHype);
112  SetDimensions(param);
113  // dz<0 can be used to force dz of hyperboloid fit the container volume
114  if (fDz<0) SetShapeBit(kGeoRunTimeShape);
115  ComputeBBox();
116 }
117 
118 ////////////////////////////////////////////////////////////////////////////////
119 /// destructor
120 
121 TGeoHype::~TGeoHype()
122 {
123 }
124 
125 ////////////////////////////////////////////////////////////////////////////////
126 /// Computes capacity of the shape in [length^3]
127 
128 Double_t TGeoHype::Capacity() const
129 {
130  Double_t capacity = 2.*TMath::Pi()*fDz*(fRmax*fRmax-fRmin*fRmin) +
131  (2.*TMath::Pi()/3.)*fDz*fDz*fDz*(fToutsq-fTinsq);
132  return capacity;
133 }
134 
135 ////////////////////////////////////////////////////////////////////////////////
136 /// Compute bounding box of the hyperboloid
137 
138 void TGeoHype::ComputeBBox()
139 {
140  if (fRmin<0.) {
141  Warning("ComputeBBox", "Shape %s has invalid rmin=%g ! SET TO 0.", GetName(),fRmin);
142  fRmin = 0.;
143  }
144  if ((fRmin>fRmax) || (fRmin*fRmin+fTinsq*fDz*fDz > fRmax*fRmax+fToutsq*fDz*fDz)) {
145  SetShapeBit(kGeoInvalidShape);
146  Error("ComputeBBox", "Shape %s hyperbolic surfaces are malformed: rin=%g, stin=%g, rout=%g, stout=%g",
147  GetName(), fRmin, fStIn, fRmax, fStOut);
148  return;
149  }
150 
151  fDX = fDY = TMath::Sqrt(RadiusHypeSq(fDz, kFALSE));
152  fDZ = fDz;
153 }
154 
155 ////////////////////////////////////////////////////////////////////////////////
156 /// Compute normal to closest surface from POINT.
157 
158 void TGeoHype::ComputeNormal(const Double_t *point, const Double_t *dir, Double_t *norm)
159 {
160  Double_t saf[3];
161  Double_t rsq = point[0]*point[0]+point[1]*point[1];
162  Double_t r = TMath::Sqrt(rsq);
163  Double_t rin = (HasInner())?(TMath::Sqrt(RadiusHypeSq(point[2],kTRUE))):0.;
164  Double_t rout = TMath::Sqrt(RadiusHypeSq(point[2],kFALSE));
165  saf[0] = TMath::Abs(fDz-TMath::Abs(point[2]));
166  saf[1] = (HasInner())?TMath::Abs(rin-r):TGeoShape::Big();
167  saf[2] = TMath::Abs(rout-r);
168  Int_t i = TMath::LocMin(3,saf);
169  if (i==0 || r<1.E-10) {
170  norm[0] = norm[1] = 0.;
171  norm[2] = TMath::Sign(1.,dir[2]);
172  return;
173  }
174  Double_t t = (i==1)?fTinsq:fToutsq;;
175  t *= -point[2]/r;
176  Double_t ct = TMath::Sqrt(1./(1.+t*t));
177  Double_t st = t * ct;
178  Double_t phi = TMath::ATan2(point[1], point[0]);
179  Double_t cphi = TMath::Cos(phi);
180  Double_t sphi = TMath::Sin(phi);
181 
182  norm[0] = ct*cphi;
183  norm[1] = ct*sphi;
184  norm[2] = st;
185  if (norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2]<0) {
186  norm[0] = -norm[0];
187  norm[1] = -norm[1];
188  norm[2] = -norm[2];
189  }
190 }
191 
192 ////////////////////////////////////////////////////////////////////////////////
193 /// test if point is inside this tube
194 
195 Bool_t TGeoHype::Contains(const Double_t *point) const
196 {
197  if (TMath::Abs(point[2]) > fDz) return kFALSE;
198  Double_t r2 = point[0]*point[0]+point[1]*point[1];
199  Double_t routsq = RadiusHypeSq(point[2], kFALSE);
200  if (r2>routsq) return kFALSE;
201  if (!HasInner()) return kTRUE;
202  Double_t rinsq = RadiusHypeSq(point[2], kTRUE);
203  if (r2<rinsq) return kFALSE;
204  return kTRUE;
205 }
206 
207 ////////////////////////////////////////////////////////////////////////////////
208 /// compute closest distance from point px,py to each corner
209 
210 Int_t TGeoHype::DistancetoPrimitive(Int_t px, Int_t py)
211 {
212  Int_t numPoints = GetNmeshVertices();
213  return ShapeDistancetoPrimitive(numPoints, px, py);
214 }
215 
216 ////////////////////////////////////////////////////////////////////////////////
217 /// Compute distance from inside point to surface of the hyperboloid.
218 
219 Double_t TGeoHype::DistFromInside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
220 {
221  if (iact<3 && safe) {
222  *safe = Safety(point, kTRUE);
223  if (iact==0) return TGeoShape::Big();
224  if ((iact==1) && (*safe>step)) return TGeoShape::Big();
225  }
226  // compute distance to surface
227  // Do Z
228  Double_t sz = TGeoShape::Big();
229  if (dir[2]>0) {
230  sz = (fDz-point[2])/dir[2];
231  if (sz<=0.) return 0.;
232  } else {
233  if (dir[2]<0) {
234  sz = -(fDz+point[2])/dir[2];
235  if (sz<=0.) return 0.;
236  }
237  }
238 
239 
240  // Do R
241  Double_t srin = TGeoShape::Big();
242  Double_t srout = TGeoShape::Big();
243  Double_t sr;
244  // inner and outer surfaces
245  Double_t s[2];
246  Int_t npos;
247  npos = DistToHype(point, dir, s, kTRUE, kTRUE);
248  if (npos) srin = s[0];
249  npos = DistToHype(point, dir, s, kFALSE, kTRUE);
250  if (npos) srout = s[0];
251  sr = TMath::Min(srin, srout);
252  return TMath::Min(sz,sr);
253 }
254 
255 
256 ////////////////////////////////////////////////////////////////////////////////
257 /// compute distance from outside point to surface of the hyperboloid.
258 
259 Double_t TGeoHype::DistFromOutside(const Double_t *point, const Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
260 {
261  if (iact<3 && safe) {
262  *safe = Safety(point, kFALSE);
263  if (iact==0) return TGeoShape::Big();
264  if ((iact==1) && (step<=*safe)) return TGeoShape::Big();
265  }
266 // Check if the bounding box is crossed within the requested distance
267  Double_t sdist = TGeoBBox::DistFromOutside(point,dir, fDX, fDY, fDZ, fOrigin, step);
268  if (sdist>=step) return TGeoShape::Big();
269  // find distance to shape
270  // Do Z
271  Double_t xi, yi, zi;
272  Double_t sz = TGeoShape::Big();
273  if (TMath::Abs(point[2])>=fDz) {
274  // We might find Z plane crossing
275  if ((point[2]*dir[2]) < 0) {
276  // Compute distance to Z (always positive)
277  sz = (TMath::Abs(point[2])-fDz)/TMath::Abs(dir[2]);
278  // Extrapolate
279  xi = point[0]+sz*dir[0];
280  yi = point[1]+sz*dir[1];
281  Double_t r2 = xi*xi + yi*yi;
282  Double_t rmin2 = RadiusHypeSq(fDz, kTRUE);
283  if (r2 >= rmin2) {
284  Double_t rmax2 = RadiusHypeSq(fDz, kFALSE);
285  if (r2 <= rmax2) return sz;
286  }
287  }
288  }
289  // We do not cross Z planes.
290  Double_t sin = TGeoShape::Big();
291  Double_t sout = TGeoShape::Big();
292  Double_t s[2];
293  Int_t npos;
294  npos = DistToHype(point, dir, s, kTRUE, kFALSE);
295  if (npos) {
296  zi = point[2] + s[0]*dir[2];
297  if (TMath::Abs(zi) <= fDz) sin = s[0];
298  else if (npos==2) {
299  zi = point[2] + s[1]*dir[2];
300  if (TMath::Abs(zi) <= fDz) sin = s[1];
301  }
302  }
303  npos = DistToHype(point, dir, s, kFALSE, kFALSE);
304  if (npos) {
305  zi = point[2] + s[0]*dir[2];
306  if (TMath::Abs(zi) <= fDz) sout = s[0];
307  else if (npos==2) {
308  zi = point[2] + s[1]*dir[2];
309  if (TMath::Abs(zi) <= fDz) sout = s[1];
310  }
311  }
312  return TMath::Min(sin, sout);
313 }
314 
315 ////////////////////////////////////////////////////////////////////////////////
316 /// Compute distance from an arbitrary point to inner/outer surface of hyperboloid.
317 /// Returns number of positive solutions. S[2] contains the solutions.
318 
319 Int_t TGeoHype::DistToHype(const Double_t *point, const Double_t *dir, Double_t *s, Bool_t inner, Bool_t in) const
320 {
321  Double_t r0, t0, snext;
322  if (inner) {
323  if (!HasInner()) return 0;
324  r0 = fRmin;
325  t0 = fTinsq;
326  } else {
327  r0 = fRmax;
328  t0 = fToutsq;
329  }
330  Double_t a = dir[0]*dir[0] + dir[1]*dir[1] - t0*dir[2]*dir[2];
331  Double_t b = t0*point[2]*dir[2] - point[0]*dir[0] - point[1]*dir[1];
332  Double_t c = point[0]*point[0] + point[1]*point[1] - t0*point[2]*point[2] - r0*r0;
333 
334  if (TMath::Abs(a) < TGeoShape::Tolerance()) {
335  if (TMath::Abs(b) < TGeoShape::Tolerance()) return 0;
336  snext = 0.5*c/b;
337  if (snext < 0.) return 0;
338  s[0] = snext;
339  return 1;
340  }
341 
342  Double_t delta = b*b - a*c;
343  Double_t ainv = 1./a;
344  Int_t npos = 0;
345  if (delta < 0.) return 0;
346  delta = TMath::Sqrt(delta);
347  Double_t sone = TMath::Sign(1.,ainv);
348  Int_t i = -1;
349  while (i<2) {
350  snext = (b + i*sone*delta)*ainv;
351  i += 2;
352  if (snext<0) continue;
353  if (snext<1.E-8) {
354  Double_t r = TMath::Sqrt(point[0]*point[0]+point[1]*point[1]);
355  Double_t t = (inner)?fTinsq:fToutsq;
356  t *= -point[2]/r;
357  Double_t phi = TMath::ATan2(point[1], point[0]);
358  Double_t ndotd = TMath::Cos(phi)*dir[0]+TMath::Sin(phi)*dir[1]+t*dir[2];
359  if (inner) ndotd *= -1;
360  if (in) ndotd *= -1;
361  if (ndotd<0) s[npos++] = snext;
362  } else s[npos++] = snext;
363  }
364  return npos;
365 }
366 
367 ////////////////////////////////////////////////////////////////////////////////
368 /// Cannot divide hyperboloids.
369 
370 TGeoVolume *TGeoHype::Divide(TGeoVolume * /*voldiv*/, const char *divname, Int_t /*iaxis*/, Int_t /*ndiv*/,
371  Double_t /*start*/, Double_t /*step*/)
372 {
373  Error("Divide", "Hyperboloids cannot be divided. Division volume %s not created", divname);
374  return 0;
375 }
376 
377 ////////////////////////////////////////////////////////////////////////////////
378 /// Get range of shape for a given axis.
379 
380 Double_t TGeoHype::GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const
381 {
382  xlo = 0;
383  xhi = 0;
384  Double_t dx = 0;
385  switch (iaxis) {
386  case 1: // R
387  xlo = fRmin;
388  xhi = TMath::Sqrt(RadiusHypeSq(fDz, kFALSE));
389  dx = xhi-xlo;
390  return dx;
391  case 2: // Phi
392  xlo = 0;
393  xhi = 360;
394  dx = 360;
395  return dx;
396  case 3: // Z
397  xlo = -fDz;
398  xhi = fDz;
399  dx = xhi-xlo;
400  return dx;
401  }
402  return dx;
403 }
404 
405 ////////////////////////////////////////////////////////////////////////////////
406 /// Fill vector param[4] with the bounding cylinder parameters. The order
407 /// is the following : Rmin, Rmax, Phi1, Phi2, dZ
408 
409 void TGeoHype::GetBoundingCylinder(Double_t *param) const
410 {
411  param[0] = fRmin; // Rmin
412  param[0] *= param[0];
413  param[1] = TMath::Sqrt(RadiusHypeSq(fDz, kFALSE)); // Rmax
414  param[1] *= param[1];
415  param[2] = 0.; // Phi1
416  param[3] = 360.; // Phi2
417 }
418 
419 ////////////////////////////////////////////////////////////////////////////////
420 /// in case shape has some negative parameters, these has to be computed
421 /// in order to fit the mother
422 
423 TGeoShape *TGeoHype::GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix * /*mat*/) const
424 {
425  if (!TestShapeBit(kGeoRunTimeShape)) return 0;
426  Double_t dz;
427  Double_t zmin,zmax;
428  dz = fDz;
429  if (fDz<0) {
430  mother->GetAxisRange(3,zmin,zmax);
431  if (zmax<0) return 0;
432  dz=zmax;
433  } else {
434  Error("GetMakeRuntimeShape", "Shape %s does not have negative Z range", GetName());
435  return 0;
436  }
437  TGeoShape *hype = new TGeoHype(GetName(), dz, fRmax, fStOut, fRmin, fStIn);
438  return hype;
439 }
440 
441 ////////////////////////////////////////////////////////////////////////////////
442 /// print shape parameters
443 
444 void TGeoHype::InspectShape() const
445 {
446  printf("*** Shape %s: TGeoHype ***\n", GetName());
447  printf(" Rin = %11.5f\n", fRmin);
448  printf(" sin = %11.5f\n", fStIn);
449  printf(" Rout = %11.5f\n", fRmax);
450  printf(" sout = %11.5f\n", fStOut);
451  printf(" dz = %11.5f\n", fDz);
452 
453  printf(" Bounding box:\n");
454  TGeoBBox::InspectShape();
455 }
456 
457 ////////////////////////////////////////////////////////////////////////////////
458 /// Creates a TBuffer3D describing *this* shape.
459 /// Coordinates are in local reference frame.
460 
461 TBuffer3D *TGeoHype::MakeBuffer3D() const
462 {
463  Int_t n = gGeoManager->GetNsegments();
464  Bool_t hasRmin = HasInner();
465  Int_t nbPnts = (hasRmin)?(2*n*n):(n*n+2);
466  Int_t nbSegs = (hasRmin)?(4*n*n):(n*(2*n+1));
467  Int_t nbPols = (hasRmin)?(2*n*n):(n*(n+1));
468 
469  TBuffer3D* buff = new TBuffer3D(TBuffer3DTypes::kGeneric,
470  nbPnts, 3*nbPnts, nbSegs, 3*nbSegs, nbPols, 6*nbPols);
471  if (buff)
472  {
473  SetPoints(buff->fPnts);
474  SetSegsAndPols(*buff);
475  }
476 
477  return buff;
478 }
479 
480 ////////////////////////////////////////////////////////////////////////////////
481 /// Fill TBuffer3D structure for segments and polygons.
482 
483 void TGeoHype::SetSegsAndPols(TBuffer3D &buff) const
484 {
485  Int_t c = GetBasicColor();
486  Int_t i, j, n;
487  n = gGeoManager->GetNsegments();
488  Bool_t hasRmin = HasInner();
489  Int_t irin = 0;
490  Int_t irout = (hasRmin)?(n*n):2;
491  // Fill segments
492  // Case hasRmin:
493  // Inner circles: [isin = 0], n (per circle) * n ( circles)
494  // iseg = isin+n*i+j , i = 0, n-1 , j = 0, n-1
495  // seg(i=1,n; j=1,n) = [irin+n*i+j] and [irin+n*i+(j+1)%n]
496  // Inner generators: [isgenin = isin+n*n], n (per circle) *(n-1) (slices)
497  // iseg = isgenin + i*n + j, i=0,n-2, j=0,n-1
498  // seg(i,j) = [irin+n*i+j] and [irin+n*(i+1)+j]
499  // Outer circles: [isout = isgenin+n*(n-1)], n (per circle) * n ( circles)
500  // iseg = isout + i*n + j , iz = 0, n-1 , j = 0, n-1
501  // seg(i=1,n; j=1,n) = [irout+n*i+j] and [irout+n*i+(j+1)%n]
502  // Outer generators: [isgenout = isout+n*n], n (per circle) *(n-1) (slices)
503  // iseg = isgenout + i*n + j, i=0,n-2, j=0,n-1
504  // seg(i,j) = [irout+n*i+j] and [irout+n*(i+1)+j]
505  // Lower cap : [islow = isgenout + n*(n-1)], n radial segments
506  // iseg = islow + j, j=0,n-1
507  // seg(j) = [irin + j] and [irout+j]
508  // Upper cap: [ishi = islow + n], nradial segments
509  // iseg = ishi + j, j=0,n-1
510  // seg[j] = [irin + n*(n-1) + j] and [irout+n*(n-1) + j]
511  //
512  // Case !hasRmin:
513  // Outer circles: [isout=0], same outer circles (n*n)
514  // Outer generators: isgenout = isout + n*n
515  // Lower cap: [islow = isgenout+n*(n-1)], n seg.
516  // iseg = islow + j, j=0,n-1
517  // seg[j] = [irin] and [irout+j]
518  // Upper cap: [ishi = islow +n]
519  // iseg = ishi + j, j=0,n-1
520  // seg[j] = [irin+1] and [irout+n*(n-1) + j]
521 
522  Int_t isin = 0;
523  Int_t isgenin = (hasRmin)?(isin+n*n):0;
524  Int_t isout = (hasRmin)?(isgenin+n*(n-1)):0;
525  Int_t isgenout = isout+n*n;
526  Int_t islo = isgenout+n*(n-1);
527  Int_t ishi = islo + n;
528 
529  Int_t npt = 0;
530  // Fill inner circle segments (n*n)
531  if (hasRmin) {
532  for (i=0; i<n; i++) {
533  for (j=0; j<n; j++) {
534  npt = 3*(isin+n*i+j);
535  buff.fSegs[npt] = c;
536  buff.fSegs[npt+1] = irin+n*i+j;
537  buff.fSegs[npt+2] = irin+n*i+((j+1)%n);
538  }
539  }
540  // Fill inner generators (n*(n-1))
541  for (i=0; i<n-1; i++) {
542  for (j=0; j<n; j++) {
543  npt = 3*(isgenin+n*i+j);
544  buff.fSegs[npt] = c;
545  buff.fSegs[npt+1] = irin+n*i+j;
546  buff.fSegs[npt+2] = irin+n*(i+1)+j;
547  }
548  }
549  }
550  // Fill outer circle segments (n*n)
551  for (i=0; i<n; i++) {
552  for (j=0; j<n; j++) {
553  npt = 3*(isout + n*i+j);
554  buff.fSegs[npt] = c;
555  buff.fSegs[npt+1] = irout+n*i+j;
556  buff.fSegs[npt+2] = irout+n*i+((j+1)%n);
557  }
558  }
559  // Fill outer generators (n*(n-1))
560  for (i=0; i<n-1; i++) {
561  for (j=0; j<n; j++) {
562  npt = 3*(isgenout+n*i+j);
563  buff.fSegs[npt] = c;
564  buff.fSegs[npt+1] = irout+n*i+j;
565  buff.fSegs[npt+2] = irout+n*(i+1)+j;
566  }
567  }
568  // Fill lower cap (n)
569  for (j=0; j<n; j++) {
570  npt = 3*(islo+j);
571  buff.fSegs[npt] = c;
572  buff.fSegs[npt+1] = irin;
573  if (hasRmin) buff.fSegs[npt+1] += j;
574  buff.fSegs[npt+2] = irout + j;
575  }
576  // Fill upper cap (n)
577  for (j=0; j<n; j++) {
578  npt = 3*(ishi+j);
579  buff.fSegs[npt] = c;
580  buff.fSegs[npt+1] = irin+1;
581  if (hasRmin) buff.fSegs[npt+1] += n*(n-1)+j-1;
582  buff.fSegs[npt+2] = irout + n*(n-1)+j;
583  }
584 
585  // Fill polygons
586  // Inner polygons: [ipin = 0] (n-1) slices * n (edges)
587  // ipoly = ipin + n*i + j; i=0,n-2 j=0,n-1
588  // poly[i,j] = [isin+n*i+j] [isgenin+i*n+(j+1)%n] [isin+n*(i+1)+j] [isgenin+i*n+j]
589  // Outer polygons: [ipout = ipin+n*(n-1)] also (n-1)*n
590  // ipoly = ipout + n*i + j; i=0,n-2 j=0,n-1
591  // poly[i,j] = [isout+n*i+j] [isgenout+i*n+j] [isout+n*(i+1)+j] [isgenout+i*n+(j+1)%n]
592  // Lower cap: [iplow = ipout+n*(n-1): n polygons
593  // ipoly = iplow + j; j=0,n-1
594  // poly[i=0,j] = [isin+j] [islow+j] [isout+j] [islow+(j+1)%n]
595  // Upper cap: [ipup = iplow+n] : n polygons
596  // ipoly = ipup + j; j=0,n-1
597  // poly[i=n-1, j] = [isin+n*(n-1)+j] [ishi+(j+1)%n] [isout+n*(n-1)+j] [ishi+j]
598  //
599  // Case !hasRmin:
600  // ipin = 0 no inner polygons
601  // ipout = 0 same outer polygons
602  // Lower cap: iplow = ipout+n*(n-1): n polygons with 3 segments
603  // poly[i=0,j] = [isout+j] [islow+(j+1)%n] [islow+j]
604  // Upper cap: ipup = iplow+n;
605  // poly[i=n-1,j] = [isout+n*(n-1)+j] [ishi+j] [ishi+(j+1)%n]
606 
607  Int_t ipin = 0;
608  Int_t ipout = (hasRmin)?(ipin+n*(n-1)):0;
609  Int_t iplo = ipout+n*(n-1);
610  Int_t ipup = iplo+n;
611  // Inner polygons n*(n-1)
612  if (hasRmin) {
613  for (i=0; i<n-1; i++) {
614  for (j=0; j<n; j++) {
615  npt = 6*(ipin+n*i+j);
616  buff.fPols[npt] = c;
617  buff.fPols[npt+1] = 4;
618  buff.fPols[npt+2] = isin+n*i+j;
619  buff.fPols[npt+3] = isgenin+i*n+((j+1)%n);
620  buff.fPols[npt+4] = isin+n*(i+1)+j;
621  buff.fPols[npt+5] = isgenin+i*n+j;
622  }
623  }
624  }
625  // Outer polygons n*(n-1)
626  for (i=0; i<n-1; i++) {
627  for (j=0; j<n; j++) {
628  npt = 6*(ipout+n*i+j);
629  buff.fPols[npt] = c;
630  buff.fPols[npt+1] = 4;
631  buff.fPols[npt+2] = isout+n*i+j;
632  buff.fPols[npt+3] = isgenout+i*n+j;
633  buff.fPols[npt+4] = isout+n*(i+1)+j;
634  buff.fPols[npt+5] = isgenout+i*n+((j+1)%n);
635  }
636  }
637  // End caps
638  if (hasRmin) {
639  for (j=0; j<n; j++) {
640  npt = 6*(iplo+j);
641  buff.fPols[npt] = c+1;
642  buff.fPols[npt+1] = 4;
643  buff.fPols[npt+2] = isin+j;
644  buff.fPols[npt+3] = islo+j;
645  buff.fPols[npt+4] = isout+j;
646  buff.fPols[npt+5] = islo+((j+1)%n);
647  }
648  for (j=0; j<n; j++) {
649  npt = 6*(ipup+j);
650  buff.fPols[npt] = c+2;
651  buff.fPols[npt+1] = 4;
652  buff.fPols[npt+2] = isin+n*(n-1)+j;
653  buff.fPols[npt+3] = ishi+((j+1)%n);
654  buff.fPols[npt+4] = isout+n*(n-1)+j;
655  buff.fPols[npt+5] = ishi+j;
656  }
657  } else {
658  for (j=0; j<n; j++) {
659  npt = 6*iplo+5*j;
660  buff.fPols[npt] = c+1;
661  buff.fPols[npt+1] = 3;
662  buff.fPols[npt+2] = isout+j;
663  buff.fPols[npt+3] = islo+((j+1)%n);
664  buff.fPols[npt+4] = islo+j;
665  }
666  for (j=0; j<n; j++) {
667  npt = 6*iplo+5*(n+j);
668  buff.fPols[npt] = c+2;
669  buff.fPols[npt+1] = 3;
670  buff.fPols[npt+2] = isout+n*(n-1)+j;
671  buff.fPols[npt+3] = ishi+j;
672  buff.fPols[npt+4] = ishi+((j+1)%n);
673  }
674  }
675 }
676 
677 ////////////////////////////////////////////////////////////////////////////////
678 /// Compute r^2 = x^2 + y^2 at a given z coordinate, for either inner or outer hyperbolas.
679 
680 Double_t TGeoHype::RadiusHypeSq(Double_t z, Bool_t inner) const
681 {
682  Double_t r0, tsq;
683  if (inner) {
684  r0 = fRmin;
685  tsq = fTinsq;
686  } else {
687  r0 = fRmax;
688  tsq = fToutsq;
689  }
690  return (r0*r0+tsq*z*z);
691 }
692 
693 ////////////////////////////////////////////////////////////////////////////////
694 /// Compute z^2 at a given r^2, for either inner or outer hyperbolas.
695 
696 Double_t TGeoHype::ZHypeSq(Double_t r, Bool_t inner) const
697 {
698  Double_t r0, tsq;
699  if (inner) {
700  r0 = fRmin;
701  tsq = fTinsq;
702  } else {
703  r0 = fRmax;
704  tsq = fToutsq;
705  }
706  if (TMath::Abs(tsq) < TGeoShape::Tolerance()) return TGeoShape::Big();
707  return ((r*r-r0*r0)/tsq);
708 }
709 
710 ////////////////////////////////////////////////////////////////////////////////
711 /// computes the closest distance from given point to this shape, according
712 /// to option. The matching point on the shape is stored in spoint.
713 
714 Double_t TGeoHype::Safety(const Double_t *point, Bool_t in) const
715 {
716  Double_t safe, safrmin, safrmax;
717  if (in) {
718  safe = fDz-TMath::Abs(point[2]);
719  safrmin = SafetyToHype(point, kTRUE, in);
720  if (safrmin < safe) safe = safrmin;
721  safrmax = SafetyToHype(point, kFALSE,in);
722  if (safrmax < safe) safe = safrmax;
723  } else {
724  safe = -fDz+TMath::Abs(point[2]);
725  safrmin = SafetyToHype(point, kTRUE, in);
726  if (safrmin > safe) safe = safrmin;
727  safrmax = SafetyToHype(point, kFALSE,in);
728  if (safrmax > safe) safe = safrmax;
729  }
730  return safe;
731 }
732 
733 ////////////////////////////////////////////////////////////////////////////////
734 /// Compute an underestimate of the closest distance from a point to inner or
735 /// outer infinite hyperbolas.
736 
737 Double_t TGeoHype::SafetyToHype(const Double_t *point, Bool_t inner, Bool_t in) const
738 {
739  Double_t r, rsq, rhsq, rh, dr, tsq, saf;
740  if (inner && !HasInner()) return (in)?TGeoShape::Big():-TGeoShape::Big();
741  rsq = point[0]*point[0]+point[1]*point[1];
742  r = TMath::Sqrt(rsq);
743  rhsq = RadiusHypeSq(point[2], inner);
744  rh = TMath::Sqrt(rhsq);
745  dr = r - rh;
746  if (inner) {
747  if (!in && dr>0) return -TGeoShape::Big();
748  if (TMath::Abs(fStIn) < TGeoShape::Tolerance()) return TMath::Abs(dr);
749  if (fRmin<TGeoShape::Tolerance()) return TMath::Abs(dr/TMath::Sqrt(1.+ fTinsq));
750  tsq = fTinsq;
751  } else {
752  if (!in && dr<0) return -TGeoShape::Big();
753  if (TMath::Abs(fStOut) < TGeoShape::Tolerance()) return TMath::Abs(dr);
754  tsq = fToutsq;
755  }
756  if (TMath::Abs(dr)<TGeoShape::Tolerance()) return 0.;
757  // 1. dr<0 => approximate safety with distance to tangent to hyperbola in z = |point[2]|
758  Double_t m;
759  if (dr<0) {
760  m = rh/(tsq*TMath::Abs(point[2]));
761  saf = -m*dr/TMath::Sqrt(1.+m*m);
762  return saf;
763  }
764  // 2. dr>0 => approximate safety with distance from point to segment P1(r(z0),z0) and P2(r0, z(r0))
765  m = (TMath::Sqrt(ZHypeSq(r,inner)) - TMath::Abs(point[2]))/dr;
766  saf = m*dr/TMath::Sqrt(1.+m*m);
767  return saf;
768 }
769 
770 ////////////////////////////////////////////////////////////////////////////////
771 /// Save a primitive as a C++ statement(s) on output stream "out".
772 
773 void TGeoHype::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
774 {
775  if (TObject::TestBit(kGeoSavePrimitive)) return;
776  out << " // Shape: " << GetName() << " type: " << ClassName() << std::endl;
777  out << " rin = " << fRmin << ";" << std::endl;
778  out << " stin = " << fStIn << ";" << std::endl;
779  out << " rout = " << fRmax << ";" << std::endl;
780  out << " stout = " << fStOut << ";" << std::endl;
781  out << " dz = " << fDz << ";" << std::endl;
782  out << " TGeoShape *" << GetPointerName() << " = new TGeoHype(\"" << GetName() << "\",rin,stin,rout,stout,dz);" << std::endl;
783  TObject::SetBit(TGeoShape::kGeoSavePrimitive);
784 }
785 
786 ////////////////////////////////////////////////////////////////////////////////
787 /// Set dimensions of the hyperboloid.
788 
789 void TGeoHype::SetHypeDimensions(Double_t rin, Double_t stin, Double_t rout, Double_t stout, Double_t dz)
790 {
791  fRmin = rin;
792  fRmax = rout;
793  fDz = dz;
794  fStIn = stin;
795  fStOut = stout;
796  fTin = TMath::Tan(fStIn*TMath::DegToRad());
797  fTinsq = fTin*fTin;
798  fTout = TMath::Tan(fStOut*TMath::DegToRad());
799  fToutsq = fTout*fTout;
800  if ((fRmin==0) && (fStIn==0)) SetShapeBit(kGeoRSeg, kTRUE);
801  else SetShapeBit(kGeoRSeg, kFALSE);
802 }
803 
804 ////////////////////////////////////////////////////////////////////////////////
805 /// Set dimensions of the hyperboloid starting from an array.
806 /// - param[0] = dz
807 /// - param[1] = rin
808 /// - param[2] = stin
809 /// - param[3] = rout
810 /// - param[4] = stout
811 
812 void TGeoHype::SetDimensions(Double_t *param)
813 {
814  Double_t dz = param[0];
815  Double_t rin = param[1];
816  Double_t stin = param[2];
817  Double_t rout = param[3];
818  Double_t stout = param[4];
819  SetHypeDimensions(rin, stin, rout, stout, dz);
820 }
821 
822 ////////////////////////////////////////////////////////////////////////////////
823 /// create tube mesh points
824 
825 void TGeoHype::SetPoints(Double_t *points) const
826 {
827  Double_t z,dz,r;
828  Int_t i,j, n;
829  if (!points) return;
830  n = gGeoManager->GetNsegments();
831  Double_t dphi = 360./n;
832  Double_t phi = 0;
833  dz = 2.*fDz/(n-1);
834 
835  Int_t indx = 0;
836 
837  if (HasInner()) {
838  // Inner surface points
839  for (i=0; i<n; i++) {
840  z = -fDz+i*dz;
841  r = TMath::Sqrt(RadiusHypeSq(z, kTRUE));
842  for (j=0; j<n; j++) {
843  phi = j*dphi*TMath::DegToRad();
844  points[indx++] = r * TMath::Cos(phi);
845  points[indx++] = r * TMath::Sin(phi);
846  points[indx++] = z;
847  }
848  }
849  } else {
850  points[indx++] = 0.;
851  points[indx++] = 0.;
852  points[indx++] = -fDz;
853  points[indx++] = 0.;
854  points[indx++] = 0.;
855  points[indx++] = fDz;
856  }
857  // Outer surface points
858  for (i=0; i<n; i++) {
859  z = -fDz + i*dz;
860  r = TMath::Sqrt(RadiusHypeSq(z, kFALSE));
861  for (j=0; j<n; j++) {
862  phi = j*dphi*TMath::DegToRad();
863  points[indx++] = r * TMath::Cos(phi);
864  points[indx++] = r * TMath::Sin(phi);
865  points[indx++] = z;
866  }
867  }
868 }
869 
870 ////////////////////////////////////////////////////////////////////////////////
871 /// create tube mesh points
872 
873 void TGeoHype::SetPoints(Float_t *points) const
874 {
875  Double_t z,dz,r;
876  Int_t i,j, n;
877  if (!points) return;
878  n = gGeoManager->GetNsegments();
879  Double_t dphi = 360./n;
880  Double_t phi = 0;
881  dz = 2.*fDz/(n-1);
882 
883  Int_t indx = 0;
884 
885  if (HasInner()) {
886  // Inner surface points
887  for (i=0; i<n; i++) {
888  z = -fDz+i*dz;
889  r = TMath::Sqrt(RadiusHypeSq(z, kTRUE));
890  for (j=0; j<n; j++) {
891  phi = j*dphi*TMath::DegToRad();
892  points[indx++] = r * TMath::Cos(phi);
893  points[indx++] = r * TMath::Sin(phi);
894  points[indx++] = z;
895  }
896  }
897  } else {
898  points[indx++] = 0.;
899  points[indx++] = 0.;
900  points[indx++] = -fDz;
901  points[indx++] = 0.;
902  points[indx++] = 0.;
903  points[indx++] = fDz;
904  }
905  // Outer surface points
906  for (i=0; i<n; i++) {
907  z = -fDz + i*dz;
908  r = TMath::Sqrt(RadiusHypeSq(z, kFALSE));
909  for (j=0; j<n; j++) {
910  phi = j*dphi*TMath::DegToRad();
911  points[indx++] = r * TMath::Cos(phi);
912  points[indx++] = r * TMath::Sin(phi);
913  points[indx++] = z;
914  }
915  }
916 }
917 
918 ////////////////////////////////////////////////////////////////////////////////
919 /// Returns numbers of vertices, segments and polygons composing the shape mesh.
920 
921 void TGeoHype::GetMeshNumbers(Int_t &nvert, Int_t &nsegs, Int_t &npols) const
922 {
923  Int_t n = gGeoManager->GetNsegments();
924  Bool_t hasRmin = HasInner();
925  nvert = (hasRmin)?(2*n*n):(n*n+2);
926  nsegs = (hasRmin)?(4*n*n):(n*(2*n+1));
927  npols = (hasRmin)?(2*n*n):(n*(n+1));
928 }
929 
930 ////////////////////////////////////////////////////////////////////////////////
931 /// Return number of vertices of the mesh representation
932 
933 Int_t TGeoHype::GetNmeshVertices() const
934 {
935  Int_t n = gGeoManager->GetNsegments();
936  Int_t numPoints = (HasRmin())?(2*n*n):(n*n+2);
937  return numPoints;
938 }
939 
940 ////////////////////////////////////////////////////////////////////////////////
941 /// fill size of this 3-D object
942 
943 void TGeoHype::Sizeof3D() const
944 {
945 }
946 
947 ////////////////////////////////////////////////////////////////////////////////
948 /// Fills a static 3D buffer and returns a reference.
949 
950 const TBuffer3D & TGeoHype::GetBuffer3D(Int_t reqSections, Bool_t localFrame) const
951 {
952  static TBuffer3D buffer(TBuffer3DTypes::kGeneric);
953 
954  TGeoBBox::FillBuffer3D(buffer, reqSections, localFrame);
955 
956  if (reqSections & TBuffer3D::kRawSizes) {
957  Int_t n = gGeoManager->GetNsegments();
958  Bool_t hasRmin = HasInner();
959  Int_t nbPnts = (hasRmin)?(2*n*n):(n*n+2);
960  Int_t nbSegs = (hasRmin)?(4*n*n):(n*(2*n+1));
961  Int_t nbPols = (hasRmin)?(2*n*n):(n*(n+1));
962  if (buffer.SetRawSizes(nbPnts, 3*nbPnts, nbSegs, 3*nbSegs, nbPols, 6*nbPols)) {
963  buffer.SetSectionsValid(TBuffer3D::kRawSizes);
964  }
965  }
966  if ((reqSections & TBuffer3D::kRaw) && buffer.SectionsValid(TBuffer3D::kRawSizes)) {
967  SetPoints(buffer.fPnts);
968  if (!buffer.fLocalFrame) {
969  TransformPoints(buffer.fPnts, buffer.NbPnts());
970  }
971 
972  SetSegsAndPols(buffer);
973  buffer.SetSectionsValid(TBuffer3D::kRaw);
974  }
975 
976  return buffer;
977 }
978 
979 ////////////////////////////////////////////////////////////////////////////////
980 /// Check the inside status for each of the points in the array.
981 /// Input: Array of point coordinates + vector size
982 /// Output: Array of Booleans for the inside of each point
983 
984 void TGeoHype::Contains_v(const Double_t *points, Bool_t *inside, Int_t vecsize) const
985 {
986  for (Int_t i=0; i<vecsize; i++) inside[i] = Contains(&points[3*i]);
987 }
988 
989 ////////////////////////////////////////////////////////////////////////////////
990 /// Compute the normal for an array o points so that norm.dot.dir is positive
991 /// Input: Arrays of point coordinates and directions + vector size
992 /// Output: Array of normal directions
993 
994 void TGeoHype::ComputeNormal_v(const Double_t *points, const Double_t *dirs, Double_t *norms, Int_t vecsize)
995 {
996  for (Int_t i=0; i<vecsize; i++) ComputeNormal(&points[3*i], &dirs[3*i], &norms[3*i]);
997 }
998 
999 ////////////////////////////////////////////////////////////////////////////////
1000 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
1001 
1002 void TGeoHype::DistFromInside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
1003 {
1004  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromInside(&points[3*i], &dirs[3*i], 3, step[i]);
1005 }
1006 
1007 ////////////////////////////////////////////////////////////////////////////////
1008 /// Compute distance from array of input points having directions specified by dirs. Store output in dists
1009 
1010 void TGeoHype::DistFromOutside_v(const Double_t *points, const Double_t *dirs, Double_t *dists, Int_t vecsize, Double_t* step) const
1011 {
1012  for (Int_t i=0; i<vecsize; i++) dists[i] = DistFromOutside(&points[3*i], &dirs[3*i], 3, step[i]);
1013 }
1014 
1015 ////////////////////////////////////////////////////////////////////////////////
1016 /// Compute safe distance from each of the points in the input array.
1017 /// Input: Array of point coordinates, array of statuses for these points, size of the arrays
1018 /// Output: Safety values
1019 
1020 void TGeoHype::Safety_v(const Double_t *points, const Bool_t *inside, Double_t *safe, Int_t vecsize) const
1021 {
1022  for (Int_t i=0; i<vecsize; i++) safe[i] = Safety(&points[3*i], inside[i]);
1023 }