```// @(#)root/geom:\$Name:  \$:\$Id: TGeoHype.cxx,v 1.4 2004/12/07 14:24:57 brun Exp \$
// Author: Mihaela Gheata   20/11/04

/*************************************************************************
*                                                                       *
* For the licensing terms see \$ROOTSYS/LICENSE.                         *
* For the list of contributors see \$ROOTSYS/README/CREDITS.             *
*************************************************************************/

#include "TROOT.h"

#include "TGeoManager.h"
#include "TGeoVolume.h"
#include "TVirtualGeoPainter.h"
#include "TGeoHype.h"
#include "TBuffer3D.h"

//_____________________________________________________________________________
// TGeoHype - Hyperboloid class defined by 5 parameters. Bounded by:
//            - Two z planes at z=+/-dz
//            - Inner and outer lateral surfaces. These represent the surfaces
//              described by the revolution of 2 hyperbolas about the Z axis:
//               r^2 - (t*z)^2 = a^2
//
//            r = distance between hyperbola and Z axis at coordinate z
//            t = tangent of the stereo angle (angle made by hyperbola
//                asimptotic lines and Z axis). t=0 means cylindrical surface.
//            a = distance between hyperbola and Z axis at z=0
//
//          The inner hyperbolic surface is described by:
//              r^2 - (tin*z)^2 = rin^2
//           - absence of the inner surface (filled hyperboloid can be forced
//             by rin=0 and sin=0
//          The outer hyperbolic surface is described by:
//              r^2 - (tout*z)^2 = rout^2
//  TGeoHype parameters: dz[cm], rin[cm], sin[deg], rout[cm], sout[deg].
//    MANDATORY conditions:
//           - rin < rout
//           - rout > 0
//           - rin^2 + (tin*dz)^2 > rout^2 + (tout*dz)^2
//    SUPPORTED CASES:
//           - rin = 0, tin != 0     => inner surface conical
//           - tin=0 AND/OR tout=0   => corresponding surface(s) cyllindrical
//             e.g. tin=0 AND tout=0 => shape becomes a tube with: rmin,rmax,dz
//
//_____________________________________________________________________________

ClassImp(TGeoHype)

//_____________________________________________________________________________
TGeoHype::TGeoHype()
{
// Default constructor
SetShapeBit(TGeoShape::kGeoHype);
fStIn = 0.;
fStOut = 0.;
fTin = 0.;
fTinsq = 0.;
fTout = 0.;
fToutsq = 0.;
}

//_____________________________________________________________________________
TGeoHype::TGeoHype(Double_t rin, Double_t stin, Double_t rout, Double_t stout, Double_t dz)
:TGeoTube(rin, rout, dz)
{
// Constructor specifying hyperboloid parameters.
SetShapeBit(TGeoShape::kGeoHype);
SetHypeDimensions(rin, stin, rout, stout, dz);
// dz<0 can be used to force dz of hyperboloid fit the container volume
if (fDz<0) SetShapeBit(kGeoRunTimeShape);
ComputeBBox();
}
//_____________________________________________________________________________
TGeoHype::TGeoHype(const char *name,Double_t rin, Double_t stin, Double_t rout, Double_t stout, Double_t dz)
:TGeoTube(name, rin, rout, dz)
{
// Constructor specifying parameters and name.
SetShapeBit(TGeoShape::kGeoHype);
SetHypeDimensions(rin, stin, rout, stout, dz);
// dz<0 can be used to force dz of hyperboloid fit the container volume
if (fDz<0) SetShapeBit(kGeoRunTimeShape);
ComputeBBox();
}

//_____________________________________________________________________________
TGeoHype::TGeoHype(Double_t *param)
:TGeoTube(param[1],param[3],param[0])
{
// Default constructor specifying a list of parameters
// param[0] = dz
// param[1] = rin
// param[2] = stin
// param[3] = rout
// param[4] = stout
SetShapeBit(TGeoShape::kGeoHype);
SetDimensions(param);
// dz<0 can be used to force dz of hyperboloid fit the container volume
if (fDz<0) SetShapeBit(kGeoRunTimeShape);
ComputeBBox();
}

//_____________________________________________________________________________
TGeoHype::~TGeoHype()
{
// destructor
}

//_____________________________________________________________________________
void TGeoHype::ComputeBBox()
{
// Compute bounding box of the hyperboloid
if (fRmin<0.) {
Warning("ComputeBBox", "Shape %s has invalid rmin=%g ! SET TO 0.", fRmin);
fRmin = 0.;
}
if ((fRmin>fRmax) || (fRmin*fRmin+fTinsq*fDz*fDz > fRmax*fRmax+fToutsq*fDz*fDz)) {
SetShapeBit(kGeoInvalidShape);
Error("ComputeBBox", "Shape %s hyperbolic surfaces are malformed: rin=%g, stin=%g, rout=%g, stout=%g",
GetName(), fRmin, fStIn, fRmax, fStOut);
return;
}

fDX = fDY = TMath::Sqrt(RadiusHypeSq(fDz, kFALSE));
fDZ = fDz;
}

//_____________________________________________________________________________
void TGeoHype::ComputeNormal(Double_t *point, Double_t *dir, Double_t *norm)
{
// Compute normal to closest surface from POINT.
Double_t saf[3];
Double_t rsq = point[0]*point[0]+point[1]*point[1];
Double_t r = TMath::Sqrt(rsq);
saf[0] = TMath::Abs(fDz-TMath::Abs(point[2]));
saf[1] = (HasInner())?TMath::Abs(rin-r):TGeoShape::Big();
saf[2] = TMath::Abs(rout-r);
Int_t i = TMath::LocMin(3,saf);
if (i==0 || r<1.E-10) {
norm[0] = norm[1] = 0.;
norm[2] = TMath::Sign(1.,dir[2]);
return;
}
Double_t t = (i==1)?fTinsq:fToutsq;;
t *= -point[2]/r;
Double_t ct = TMath::Sqrt(1./(1.+t*t));
Double_t st = t * ct;
Double_t phi = TMath::ATan2(point[1], point[0]);
Double_t cphi = TMath::Cos(phi);
Double_t sphi = TMath::Sin(phi);

norm[0] = ct*cphi;
norm[1] = ct*sphi;
norm[2] = st;
if (norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2]<0) {
norm[0] = -norm[0];
norm[1] = -norm[1];
norm[2] = -norm[2];
}
}

//_____________________________________________________________________________
Bool_t TGeoHype::Contains(Double_t *point) const
{
// test if point is inside this tube
if (TMath::Abs(point[2]) > fDz) return kFALSE;
Double_t r2 = point[0]*point[0]+point[1]*point[1];
if (r2>routsq) return kFALSE;
if (!HasInner()) return kTRUE;
if (r2<rinsq) return kFALSE;
return kTRUE;
}

//_____________________________________________________________________________
Int_t TGeoHype::DistancetoPrimitive(Int_t px, Int_t py)
{
// compute closest distance from point px,py to each corner
Int_t numPoints = GetNmeshVertices();
return ShapeDistancetoPrimitive(numPoints, px, py);
}

//_____________________________________________________________________________
Double_t TGeoHype::DistFromInside(Double_t *point, Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
{
// Compute distance from inside point to surface of the hyperboloid.
if (iact<3 && safe) {
*safe = Safety(point, kTRUE);
if (iact==0) return TGeoShape::Big();
if ((iact==1) && (*safe>step)) return TGeoShape::Big();
}
// compute distance to surface
// Do Z
Double_t sz = TGeoShape::Big();
if (dir[2]>0) {
sz = (fDz-point[2])/dir[2];
if (sz<=0.) return 0.;
} else {
if (dir[2]<0) {
sz = -(fDz+point[2])/dir[2];
if (sz<=0.) return 0.;
}
}

// Do R
Double_t srin = TGeoShape::Big();
Double_t srout = TGeoShape::Big();
Double_t sr;
// inner and outer surfaces
Double_t s[2];
Int_t npos;
npos = DistToHype(point, dir, s, kTRUE);
if (npos) srin = s[0];
npos = DistToHype(point, dir, s, kFALSE);
if (npos) srout = s[0];
sr = TMath::Min(srin, srout);
return TMath::Min(sz,sr);
}

//_____________________________________________________________________________
Double_t TGeoHype::DistFromOutside(Double_t *point, Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
{
// compute distance from outside point to surface of the hyperboloid.
if (iact<3 && safe) {
*safe = Safety(point, kFALSE);
if (iact==0) return TGeoShape::Big();
if ((iact==1) && (step<=*safe)) return TGeoShape::Big();
}
// find distance to shape
// Do Z
Double_t xi, yi, zi;
Double_t sz = TGeoShape::Big();
if (TMath::Abs(point[2])>=fDz) {
// We might find Z plane crossing
if ((point[2]*dir[2]) < 0) {
// Compute distance to Z (always positive)
sz = (TMath::Abs(point[2])-fDz)/TMath::Abs(dir[2]);
// Extrapolate
xi = point[0]+sz*dir[0];
yi = point[1]+sz*dir[1];
Double_t r2 = xi*xi + yi*yi;
if (r2 >= rmin2) {
if (r2 <= rmax2) return sz;
}
}
}
// We do not cross Z planes.
Double_t sin = TGeoShape::Big();
Double_t sout = TGeoShape::Big();
Double_t s[2];
Int_t npos;
npos = DistToHype(point, dir, s, kTRUE);
if (npos) {
zi = point[2] + s[0]*dir[2];
if (TMath::Abs(zi) <= fDz) sin = s[0];
else if (npos==2) {
zi = point[2] + s[1]*dir[2];
if (TMath::Abs(zi) <= fDz) sin = s[1];
}
}
npos = DistToHype(point, dir, s, kFALSE);
if (npos) {
zi = point[2] + s[0]*dir[2];
if (TMath::Abs(zi) <= fDz) sout = s[0];
else if (npos==2) {
zi = point[2] + s[1]*dir[2];
if (TMath::Abs(zi) <= fDz) sout = s[1];
}
}
return TMath::Min(sin, sout);
}

//_____________________________________________________________________________
Int_t TGeoHype::DistToHype(Double_t *point, Double_t *dir, Double_t *s, Bool_t inner) const
{
// Compute distance from an arbitrary point to inner/outer surface of hyperboloid.
// Returns number of positive solutions. S[2] contains the solutions.
Double_t r0, t0, snext;
if (inner) {
if (!HasInner()) return 0;
r0 = fRmin;
t0 = fTinsq;
} else {
r0 = fRmax;
t0 = fToutsq;
}
Double_t a = dir[0]*dir[0] + dir[1]*dir[1] - t0*dir[2]*dir[2];
Double_t b = t0*point[2]*dir[2] - point[0]*dir[0] - point[1]*dir[1];
Double_t c = point[0]*point[0] + point[1]*point[1] - t0*point[2]*point[2] - r0*r0;

if (a == 0.) {
if (b == 0.) return 0;
snext = 0.5*c/b;
if (snext < 0.) return 0;
s[0] = snext;
return 1;
}

Double_t delta = b*b - a*c;
Double_t ainv = 1./a;
Int_t npos = 0;
if (delta < 0.) return 0;
delta = TMath::Sqrt(delta);
Double_t sone = TMath::Sign(1.,ainv);
snext = (b - sone*delta)*ainv;
if (snext >= 0.) s[npos++] = snext;
snext = (b + sone*delta)*ainv;
if (snext >= 0.) s[npos++] = snext;
return npos;
}

//_____________________________________________________________________________
TGeoVolume *TGeoHype::Divide(TGeoVolume * /*voldiv*/, const char *divname, Int_t /*iaxis*/, Int_t /*ndiv*/,
Double_t /*start*/, Double_t /*step*/)
{
// Cannot divide hyperboloids.
Error("Divide", "Hyperboloids cannot be divided. Division volume %s not created", divname);
return 0;
}

//_____________________________________________________________________________
Double_t TGeoHype::GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const
{
// Get range of shape for a given axis.
xlo = 0;
xhi = 0;
Double_t dx = 0;
switch (iaxis) {
case 1: // R
xlo = fRmin;
dx = xhi-xlo;
return dx;
case 2: // Phi
xlo = 0;
xhi = 360;
dx = 360;
return dx;
case 3: // Z
xlo = -fDz;
xhi = fDz;
dx = xhi-xlo;
return dx;
}
return dx;
}

//_____________________________________________________________________________
void TGeoHype::GetBoundingCylinder(Double_t *param) const
{
//--- Fill vector param[4] with the bounding cylinder parameters. The order
// is the following : Rmin, Rmax, Phi1, Phi2, dZ
param[0] = fRmin; // Rmin
param[0] *= param[0];
param[1] = TMath::Sqrt(RadiusHypeSq(fDz, kFALSE)); // Rmax
param[1] *= param[1];
param[2] = 0.;    // Phi1
param[3] = 360.;  // Phi1
}

//_____________________________________________________________________________
TGeoShape *TGeoHype::GetMakeRuntimeShape(TGeoShape *mother, TGeoMatrix * /*mat*/) const
{
// in case shape has some negative parameters, these has to be computed
// in order to fit the mother
if (!TestShapeBit(kGeoRunTimeShape)) return 0;
Double_t rmin, rmax, dz;
Double_t zmin,zmax;
rmin = fRmin;
rmax = fRmax;
dz = fDz;
if (fDz<0) {
mother->GetAxisRange(3,zmin,zmax);
if (zmax<0) return 0;
dz=zmax;
} else {
Error("GetMakeRuntimeShape", "Shape %s does not have negative Z range", GetName());
return 0;
}
TGeoShape *hype = new TGeoHype(GetName(), dz, fRmax, fStOut, fRmin, fStIn);
return hype;
}

//_____________________________________________________________________________
void TGeoHype::InspectShape() const
{
// print shape parameters
printf("*** Shape %s: TGeoHype ***\n", GetName());
printf("    Rin  = %11.5f\n", fRmin);
printf("    sin  = %11.5f\n", fStIn);
printf("    Rout = %11.5f\n", fRmax);
printf("    sout = %11.5f\n", fStOut);
printf("    dz   = %11.5f\n", fDz);

printf(" Bounding box:\n");
TGeoBBox::InspectShape();
}

//_____________________________________________________________________________
TBuffer3D *TGeoHype::MakeBuffer3D() const
{
// Creates a TBuffer3D describing *this* shape.
// Coordinates are in local reference frame.

Int_t n = gGeoManager->GetNsegments();
Bool_t hasRmin = HasInner();
Int_t NbPnts = (hasRmin)?(2*n*n):(n*n+2);
Int_t NbSegs = (hasRmin)?(4*n*n):(n*(2*n+1));
Int_t NbPols = (hasRmin)?(2*n*n):(n*(n+1));

TBuffer3D* buff = new TBuffer3D(3*NbPnts, 3*NbSegs, 6*NbPols);

buff->fType = TBuffer3D::kANY; // should be kHYPE
buff->fNbPnts = NbPnts;
buff->fNbSegs = NbSegs;
buff->fNbPols = NbPols;

SetPoints(buff->fPnts);

SetSegsAndPols(buff);

return buff;
}

//_____________________________________________________________________________
void TGeoHype::Paint(Option_t *option)
{
// Paint this shape according to option

// Allocate the necessary spage in gPad->fBuffer3D to store this shape
Int_t n = gGeoManager->GetNsegments();
Bool_t hasRmin = HasInner();
Int_t NbPnts = (hasRmin)?(2*n*n):(n*n+2);
Int_t NbSegs = (hasRmin)?(4*n*n):(n*(2*n+1));
Int_t NbPols = (hasRmin)?(2*n*n):(n*(n+1));
TBuffer3D *buff = gPad->AllocateBuffer3D(3*NbPnts, 3*NbSegs, 6*NbPols);

buff->fType = TBuffer3D::kANY; // should be kHYPE
TGeoVolume *vol = gGeoManager->GetPaintVolume();
buff->fId   = vol;

// Fill gPad->fBuffer3D. Points coordinates are in Master space
buff->fNbPnts = NbPnts;
buff->fNbSegs = NbSegs;
buff->fNbPols = NbPols;
// In case of option "size" it is not necessary to fill the buffer
if (strstr(option,"size")) {
buff->Paint(option);
return;
}
// Fill points
// Case hasRmin:
//   irin = 0 , n (per circle) * n (circles) starting with z = -fDz
//            icin(j) = irin + j*n  (j=0,n-1) = index of first pt. on circle j
//   irout = n*n , n (per circle) * n (circles) starting with z = -fDz
//            icout(j) = irout + j*n  (j=0,n-1) = index of first pt. on circle j
// Case !hasRmin:
//   irin = 0, 2 points (lower and upper centers)
//   irout = 2
SetPoints(buff->fPnts);
TransformPoints(buff);

// Basic colors: 0, 1, ... 7
buff->fColor = vol->GetLineColor();
SetSegsAndPols(buff);
buff->Paint(option);
}

//_____________________________________________________________________________
void TGeoHype::SetSegsAndPols(TBuffer3D *buff) const
{
// Fill TBuffer3D structure for segments and polygons.
Int_t c = (((buff->fColor) %8) -1) * 4;
if (c < 0) c = 0;
Int_t i, j, n;
n = gGeoManager->GetNsegments();
Bool_t hasRmin = HasInner();
Int_t irin = 0;
Int_t irout = (hasRmin)?(n*n):2;
// Fill segments
// Case hasRmin:
//   Inner circles:  [isin = 0], n (per circle) * n ( circles)
//        iseg = isin+n*i+j , i = 0, n-1   , j = 0, n-1
//        seg(i=1,n; j=1,n) = [irin+n*i+j] and [irin+n*i+(j+1)%n]
//   Inner generators: [isgenin = isin+n*n], n (per circle) *(n-1) (slices)
//        iseg = isgenin + i*n + j, i=0,n-2,  j=0,n-1
//        seg(i,j) = [irin+n*i+j] and [irin+n*(i+1)+j]
//   Outer circles:  [isout = isgenin+n*(n-1)], n (per circle) * n ( circles)
//        iseg = isout + i*n + j , iz = 0, n-1   , j = 0, n-1
//        seg(i=1,n; j=1,n) = [irout+n*i+j] and [irout+n*i+(j+1)%n]
//   Outer generators: [isgenout = isout+n*n], n (per circle) *(n-1) (slices)
//        iseg = isgenout + i*n + j, i=0,n-2,  j=0,n-1
//        seg(i,j) = [irout+n*i+j] and [irout+n*(i+1)+j]
//   Lower cap : [islow = isgenout + n*(n-1)], n radial segments
//        iseg = islow + j,  j=0,n-1
//        seg(j) = [irin + j] and [irout+j]
//   Upper cap: [ishi = islow + n], nradial segments
//        iseg = ishi + j, j=0,n-1
//        seg[j] = [irin + n*(n-1) + j] and [irout+n*(n-1) + j]
//
// Case !hasRmin:
//   Outer circles: [isout=0], same outer circles (n*n)
// Outer generators: isgenout = isout + n*n
//   Lower cap: [islow = isgenout+n*(n-1)], n seg.
//        iseg = islow + j, j=0,n-1
//        seg[j] = [irin] and [irout+j]
//   Upper cap: [ishi = islow +n]
//        iseg = ishi + j, j=0,n-1
//        seg[j] = [irin+1] and [irout+n*(n-1) + j]

Int_t isin = 0;
Int_t isgenin = (hasRmin)?(isin+n*n):0;
Int_t isout = (hasRmin)?(isgenin+n*(n-1)):0;
Int_t isgenout  = isout+n*n;
Int_t islo = isgenout+n*(n-1);
Int_t ishi = islo + n;

Int_t npt = 0;
// Fill inner circle segments (n*n)
if (hasRmin) {
for (i=0; i<n; i++) {
for (j=0; j<n; j++) {
npt = 3*(isin+n*i+j);
buff->fSegs[npt]   = c;
buff->fSegs[npt+1] = irin+n*i+j;
buff->fSegs[npt+2] = irin+n*i+((j+1)%n);
}
}
// Fill inner generators (n*(n-1))
for (i=0; i<n-1; i++) {
for (j=0; j<n; j++) {
npt = 3*(isgenin+n*i+j);
buff->fSegs[npt]   = c;
buff->fSegs[npt+1] = irin+n*i+j;
buff->fSegs[npt+2] = irin+n*(i+1)+j;
}
}
}
// Fill outer circle segments (n*n)
for (i=0; i<n; i++) {
for (j=0; j<n; j++) {
npt = 3*(isout + n*i+j);
buff->fSegs[npt]   = c;
buff->fSegs[npt+1] = irout+n*i+j;
buff->fSegs[npt+2] = irout+n*i+((j+1)%n);
}
}
// Fill outer generators (n*(n-1))
for (i=0; i<n-1; i++) {
for (j=0; j<n; j++) {
npt = 3*(isgenout+n*i+j);
buff->fSegs[npt]   = c;
buff->fSegs[npt+1] = irout+n*i+j;
buff->fSegs[npt+2] = irout+n*(i+1)+j;
}
}
// Fill lower cap (n)
for (j=0; j<n; j++) {
npt = 3*(islo+j);
buff->fSegs[npt]   = c;
buff->fSegs[npt+1] = irin;
if (hasRmin) buff->fSegs[npt+1] += j;
buff->fSegs[npt+2] = irout + j;
}
// Fill upper cap (n)
for (j=0; j<n; j++) {
npt = 3*(ishi+j);
buff->fSegs[npt]   = c;
buff->fSegs[npt+1] = irin+1;
if (hasRmin) buff->fSegs[npt+1] += n*(n-1)+j-1;
buff->fSegs[npt+2] = irout + n*(n-1)+j;
}

// Fill polygons
// Inner polygons: [ipin = 0] (n-1) slices * n (edges)
//   ipoly = ipin + n*i + j;  i=0,n-2   j=0,n-1
//   poly[i,j] = [isin+n*i+j]  [isgenin+i*n+(j+1)%n]  [isin+n*(i+1)+j]  [isgenin+i*n+j]
// Outer polygons: [ipout = ipin+n*(n-1)]  also (n-1)*n
//   ipoly = ipout + n*i + j; i=0,n-2   j=0,n-1
//   poly[i,j] = [isout+n*i+j]  [isgenout+i*n+j]  [isout+n*(i+1)+j]  [isgenout+i*n+(j+1)%n]
// Lower cap: [iplow = ipout+n*(n-1):  n polygons
//   ipoly = iplow + j;  j=0,n-1
//   poly[i=0,j] = [isin+j] [islow+j] [isout+j] [islow+(j+1)%n]
// Upper cap: [ipup = iplow+n] : n polygons
//   ipoly = ipup + j;  j=0,n-1
//   poly[i=n-1, j] = [isin+n*(n-1)+j] [ishi+(j+1)%n] [isout+n*(n-1)+j] [ishi+j]
//
// Case !hasRmin:
// ipin = 0 no inner polygons
// ipout = 0 same outer polygons
// Lower cap: iplow = ipout+n*(n-1):  n polygons with 3 segments
//   poly[i=0,j] = [isout+j] [islow+(j+1)%n] [islow+j]
// Upper cap: ipup = iplow+n;
//   poly[i=n-1,j] = [isout+n*(n-1)+j] [ishi+j] [ishi+(j+1)%n]

Int_t ipin = 0;
Int_t ipout = (hasRmin)?(ipin+n*(n-1)):0;
Int_t iplo = ipout+n*(n-1);
Int_t ipup = iplo+n;
// Inner polygons n*(n-1)
if (hasRmin) {
for (i=0; i<n-1; i++) {
for (j=0; j<n; j++) {
npt = 6*(ipin+n*i+j);
buff->fPols[npt]   = c;
buff->fPols[npt+1] = 4;
buff->fPols[npt+2] = isin+n*i+j;
buff->fPols[npt+3] = isgenin+i*n+((j+1)%n);
buff->fPols[npt+4] = isin+n*(i+1)+j;
buff->fPols[npt+5] = isgenin+i*n+j;
}
}
}
// Outer polygons n*(n-1)
for (i=0; i<n-1; i++) {
for (j=0; j<n; j++) {
npt = 6*(ipout+n*i+j);
buff->fPols[npt]   = c;
buff->fPols[npt+1] = 4;
buff->fPols[npt+2] = isout+n*i+j;
buff->fPols[npt+3] = isgenout+i*n+j;
buff->fPols[npt+4] = isout+n*(i+1)+j;
buff->fPols[npt+5] = isgenout+i*n+((j+1)%n);
}
}
// End caps
if (hasRmin) {
for (j=0; j<n; j++) {
npt = 6*(iplo+j);
buff->fPols[npt]   = c+1;
buff->fPols[npt+1] = 4;
buff->fPols[npt+2] = isin+j;
buff->fPols[npt+3] = islo+j;
buff->fPols[npt+4] = isout+j;
buff->fPols[npt+5] = islo+((j+1)%n);
}
for (j=0; j<n; j++) {
npt = 6*(ipup+j);
buff->fPols[npt]   = c+2;
buff->fPols[npt+1] = 4;
buff->fPols[npt+2] = isin+n*(n-1)+j;
buff->fPols[npt+3] = ishi+((j+1)%n);
buff->fPols[npt+4] = isout+n*(n-1)+j;
buff->fPols[npt+5] = ishi+j;
}
} else {
for (j=0; j<n; j++) {
npt = 6*iplo+5*j;
buff->fPols[npt]   = c+1;
buff->fPols[npt+1] = 3;
buff->fPols[npt+2] = isout+j;
buff->fPols[npt+3] = islo+((j+1)%n);
buff->fPols[npt+4] = islo+j;
}
for (j=0; j<n; j++) {
npt = 6*iplo+5*(n+j);
buff->fPols[npt]   = c+2;
buff->fPols[npt+1] = 3;
buff->fPols[npt+2] = isout+n*(n-1)+j;
buff->fPols[npt+3] = ishi+j;
buff->fPols[npt+4] = ishi+((j+1)%n);
}
}
}

//_____________________________________________________________________________
Double_t TGeoHype::RadiusHypeSq(Double_t z, Bool_t inner) const
{
// Compute r^2 = x^2 + y^2 at a given z coordinate, for either inner or outer hyperbolas.
Double_t r0, tsq;
if (inner) {
r0 = fRmin;
tsq = fTinsq;
} else {
r0 = fRmax;
tsq = fToutsq;
}
return (r0*r0+tsq*z*z);
}

//_____________________________________________________________________________
Double_t TGeoHype::ZHypeSq(Double_t r, Bool_t inner) const
{
// Compute z^2 at a given  r^2, for either inner or outer hyperbolas.
Double_t r0, tsq;
if (inner) {
r0 = fRmin;
tsq = fTinsq;
} else {
r0 = fRmax;
tsq = fToutsq;
}
if (tsq==0) return TGeoShape::Big();
return ((r*r-r0*r0)/tsq);
}

//_____________________________________________________________________________
Double_t TGeoHype::Safety(Double_t *point, Bool_t in) const
{
// computes the closest distance from given point to this shape, according
// to option. The matching point on the shape is stored in spoint.
Double_t safe, safrmin, safrmax;
if (in) {
safe    = fDz-TMath::Abs(point[2]);
safrmin = SafetyToHype(point, kTRUE, in);
if (safrmin < safe) safe = safrmin;
safrmax = SafetyToHype(point, kFALSE,in);
if (safrmax < safe) safe = safrmax;
} else {
safe    = -fDz+TMath::Abs(point[2]);
safrmin = SafetyToHype(point, kTRUE, in);
if (safrmin > safe) safe = safrmin;
safrmax = SafetyToHype(point, kFALSE,in);
if (safrmax > safe) safe = safrmax;
}
return safe;
}

//_____________________________________________________________________________
Double_t TGeoHype::SafetyToHype(Double_t *point, Bool_t inner, Bool_t in) const
{
// Compute an underestimate of the closest distance from a point to inner or
// outer infinite hyperbolas.
Double_t r, rsq, rhsq, rh, dr, tsq, saf;
if (inner && !HasInner()) return (in)?TGeoShape::Big():-TGeoShape::Big();
rsq = point[0]*point[0]+point[1]*point[1];
r = TMath::Sqrt(rsq);
rh = TMath::Sqrt(rhsq);
dr = r - rh;
if (inner) {
if (!in && dr>0) return -TGeoShape::Big();
if (fStIn == 0) return TMath::Abs(dr);
if (fRmin==0) return TMath::Abs(dr/TMath::Sqrt(1.+ fTinsq));
tsq = fTinsq;
} else {
if (!in && dr<0) return -TGeoShape::Big();
if (fStOut == 0) return TMath::Abs(dr);
tsq = fToutsq;
}
if (dr==0) return 0.;
// 1. dr<0 => approximate safety with distance to tangent to hyperbola in z = |point[2]|
Double_t m;
if (dr<0) {
m = rh/(tsq*TMath::Abs(point[2]));
saf = -m*dr/TMath::Sqrt(1.+m*m);
return saf;
}
// 2. dr>0 => approximate safety with distance from point to segment P1(r(z0),z0) and P2(r0, z(r0))
m = (TMath::Sqrt(ZHypeSq(r,inner)) - TMath::Abs(point[2]))/dr;
saf = m*dr/TMath::Sqrt(1.+m*m);
return saf;
}

//_____________________________________________________________________________
void TGeoHype::SetHypeDimensions(Double_t rin, Double_t stin, Double_t rout, Double_t stout, Double_t dz)
{
fRmin = rin;
fRmax = rout;
fDz   = dz;
fStIn = stin;
fStOut = stout;
fTinsq = fTin*fTin;
fToutsq = fTout*fTout;
if ((fRmin==0) && (fStIn==0)) SetShapeBit(kGeoRSeg, kTRUE);
else                          SetShapeBit(kGeoRSeg, kFALSE);
}

//_____________________________________________________________________________
void TGeoHype::SetDimensions(Double_t *param)
{
// param[0] = dz
// param[1] = rin
// param[2] = stin
// param[3] = rout
// param[4] = stout
Double_t dz = param[0];
Double_t rin = param[1];
Double_t stin = param[2];
Double_t rout = param[3];
Double_t stout = param[4];
SetHypeDimensions(rin, stin, rout, stout, dz);
}

//_____________________________________________________________________________
void TGeoHype::SetPoints(Double_t *buff) const
{
// create tube mesh points
Double_t z,dz,r;
Int_t i,j, n;
if (!buff) return;
n = gGeoManager->GetNsegments();
Double_t dphi = 360./n;
Double_t phi = 0;
dz = 2.*fDz/(n-1);

Int_t indx = 0;

if (HasInner()) {
// Inner surface points
for (i=0; i<n; i++) {
z = -fDz+i*dz;
for (j=0; j<n; j++) {
buff[indx++] = r * TMath::Cos(phi);
buff[indx++] = r * TMath::Sin(phi);
buff[indx++] = z;
}
}
} else {
buff[indx++] = 0.;
buff[indx++] = 0.;
buff[indx++] = -fDz;
buff[indx++] = 0.;
buff[indx++] = 0.;
buff[indx++] = fDz;
}
// Outer surface points
for (i=0; i<n; i++) {
z = -fDz + i*dz;
for (j=0; j<n; j++) {
buff[indx++] = r * TMath::Cos(phi);
buff[indx++] = r * TMath::Sin(phi);
buff[indx++] = z;
}
}
}

//_____________________________________________________________________________
void TGeoHype::SetPoints(Float_t *buff) const
{
// create tube mesh points
Double_t z,dz,r;
Int_t i,j, n;
if (!buff) return;
n = gGeoManager->GetNsegments();
Double_t dphi = 360./n;
Double_t phi = 0;
dz = 2.*fDz/(n-1);

Int_t indx = 0;

if (HasInner()) {
// Inner surface points
for (i=0; i<n; i++) {
z = -fDz+i*dz;
for (j=0; j<n; j++) {
buff[indx++] = r * TMath::Cos(phi);
buff[indx++] = r * TMath::Sin(phi);
buff[indx++] = z;
}
}
} else {
buff[indx++] = 0.;
buff[indx++] = 0.;
buff[indx++] = -fDz;
buff[indx++] = 0.;
buff[indx++] = 0.;
buff[indx++] = fDz;
}
// Outer surface points
for (i=0; i<n; i++) {
z = -fDz + i*dz;
for (j=0; j<n; j++) {
buff[indx++] = r * TMath::Cos(phi);
buff[indx++] = r * TMath::Sin(phi);
buff[indx++] = z;
}
}
}

//_____________________________________________________________________________
Int_t TGeoHype::GetNmeshVertices() const
{
// Return number of vertices of the mesh representation
Int_t n = gGeoManager->GetNsegments();
Int_t numPoints = (HasRmin())?(2*n*n):(n*n+2);
return numPoints;
}

//_____________________________________________________________________________
void TGeoHype::Sizeof3D() const
{
///// fill size of this 3-D object
///    TVirtualGeoPainter *painter = gGeoManager->GetGeomPainter();
///    if (!painter) return;
///    Int_t n = gGeoManager->GetNsegments();
///    Int_t numPoints = n*4;
///    Int_t numSegs   = n*8;
///    Int_t numPolys  = n*4;