#include "Riostream.h"
#include "TROOT.h"
#include "TClass.h"
#include "TGeoManager.h"
#include "TGeoVolume.h"
#include "TGeoPolygon.h"
#include "TVirtualGeoPainter.h"
#include "TGeoXtru.h"
#include "TVirtualPad.h"
#include "TBuffer3D.h"
#include "TBuffer3DTypes.h"
ClassImp(TGeoXtru)
TGeoXtru::TGeoXtru()
{
SetShapeBit(TGeoShape::kGeoXtru);
fNvert = 0;
fNz = 0;
fZcurrent = 0.;
fPoly = 0;
fX = 0;
fY = 0;
fXc = 0;
fYc = 0;
fZ = 0;
fScale = 0;
fX0 = 0;
fY0 = 0;
fSeg = 0;
fIz = 0;
}
TGeoXtru::TGeoXtru(Int_t nz)
:TGeoBBox(0, 0, 0)
{
SetShapeBit(TGeoShape::kGeoXtru);
if (nz<2) {
Error("ctor", "Cannot create TGeoXtru %s with less than 2 Z planes", GetName());
return;
}
fNvert = 0;
fNz = nz;
fZcurrent = 0.;
fPoly = 0;
fX = 0;
fY = 0;
fXc = 0;
fYc = 0;
fZ = new Double_t[nz];
fScale = new Double_t[nz];
fX0 = new Double_t[nz];
fY0 = new Double_t[nz];
fSeg = 0;
fIz = 0;
}
TGeoXtru::TGeoXtru(Double_t *param)
:TGeoBBox(0, 0, 0)
{
SetShapeBit(TGeoShape::kGeoXtru);
fNvert = 0;
fNz = 0;
fZcurrent = 0.;
fPoly = 0;
fX = 0;
fY = 0;
fXc = 0;
fYc = 0;
fZ = 0;
fScale = 0;
fX0 = 0;
fY0 = 0;
fSeg = 0;
fIz = 0;
SetDimensions(param);
}
TGeoXtru::TGeoXtru(const TGeoXtru& xt) :
TGeoBBox(xt),
fNvert(xt.fNvert),
fNz(xt.fNz),
fZcurrent(xt.fZcurrent),
fPoly(xt.fPoly),
fX(xt.fX),
fY(xt.fY),
fXc(xt.fXc),
fYc(xt.fYc),
fZ(xt.fZ),
fScale(xt.fScale),
fX0(xt.fX0),
fY0(xt.fY0),
fSeg(xt.fSeg),
fIz(xt.fIz)
{
}
TGeoXtru& TGeoXtru::operator=(const TGeoXtru& xt)
{
if(this!=&xt) {
TGeoBBox::operator=(xt);
fNvert=xt.fNvert;
fNz=xt.fNz;
fZcurrent=xt.fZcurrent;
fPoly=xt.fPoly;
fX=xt.fX;
fY=xt.fY;
fXc=xt.fXc;
fYc=xt.fYc;
fZ=xt.fZ;
fScale=xt.fScale;
fX0=xt.fX0;
fY0=xt.fY0;
fSeg=xt.fSeg;
fIz=xt.fIz;
}
return *this;
}
TGeoXtru::~TGeoXtru()
{
if (fX) {delete[] fX; fX = 0;}
if (fY) {delete[] fY; fY = 0;}
if (fXc) {delete[] fXc; fXc = 0;}
if (fYc) {delete[] fYc; fYc = 0;}
if (fZ) {delete[] fZ; fZ = 0;}
if (fScale) {delete[] fScale; fScale = 0;}
if (fX0) {delete[] fX0; fX0 = 0;}
if (fY0) {delete[] fY0; fY0 = 0;}
}
Double_t TGeoXtru::Capacity() const
{
Int_t iz;
Double_t capacity = 0;
Double_t area, dz, sc1, sc2;
TGeoXtru *xtru = (TGeoXtru*)this;
xtru->SetCurrentVertices(0.,0.,1.);
area = fPoly->Area();
for (iz=0; iz<fNz-1; iz++) {
dz = fZ[iz+1]-fZ[iz];
if (dz==0) continue;
sc1 = fScale[iz];
sc2 = fScale[iz+1];
capacity += (area*dz/3.)*(sc1*sc1+sc1*sc2+sc2*sc2);
}
return capacity;
}
void TGeoXtru::ComputeBBox()
{
if (!fX || !fZ || !fNvert) {
Error("ComputeBBox", "In shape %s polygon not defined", GetName());
return;
}
Double_t zmin = fZ[0];
Double_t zmax = fZ[fNz-1];
Double_t xmin = TGeoShape::Big();
Double_t xmax = -TGeoShape::Big();
Double_t ymin = TGeoShape::Big();
Double_t ymax = -TGeoShape::Big();
for (Int_t i=0; i<fNz; i++) {
SetCurrentVertices(fX0[i], fY0[i], fScale[i]);
for (Int_t j=0; j<fNvert; j++) {
if (fXc[j]<xmin) xmin=fXc[j];
if (fXc[j]>xmax) xmax=fXc[j];
if (fYc[j]<ymin) ymin=fYc[j];
if (fYc[j]>ymax) ymax=fYc[j];
}
}
fOrigin[0] = 0.5*(xmin+xmax);
fOrigin[1] = 0.5*(ymin+ymax);
fOrigin[2] = 0.5*(zmin+zmax);
fDX = 0.5*(xmax-xmin);
fDY = 0.5*(ymax-ymin);
fDZ = 0.5*(zmax-zmin);
}
void TGeoXtru::ComputeNormal(Double_t * , Double_t *dir, Double_t *norm)
{
if (fIz<0) {
memset(norm,0,3*sizeof(Double_t));
norm[2] = (dir[2]>0)?1:-1;
return;
}
Double_t vert[12];
GetPlaneVertices(fIz, fSeg, vert);
GetPlaneNormal(vert, norm);
Double_t ndotd = norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2];
if (ndotd<0) {
norm[0] = -norm[0];
norm[1] = -norm[1];
norm[2] = -norm[2];
}
}
Bool_t TGeoXtru::Contains(Double_t *point) const
{
TGeoXtru *xtru = (TGeoXtru*)this;
if (point[2]<fZ[0]) return kFALSE;
if (point[2]>fZ[fNz-1]) return kFALSE;
Int_t iz = TMath::BinarySearch(fNz, fZ, point[2]);
if (point[2]==fZ[iz]) {
xtru->SetIz(-1);
xtru->SetCurrentVertices(fX0[iz],fY0[iz], fScale[iz]);
if (fPoly->Contains(point)) return kTRUE;
if (iz>1 && fZ[iz]==fZ[iz-1]) {
xtru->SetCurrentVertices(fX0[iz-1],fY0[iz-1], fScale[iz-1]);
return fPoly->Contains(point);
} else if (iz<fNz-2 && fZ[iz]==fZ[iz+1]) {
xtru->SetCurrentVertices(fX0[iz+1],fY0[iz+1], fScale[iz+1]);
return fPoly->Contains(point);
}
}
xtru->SetCurrentZ(point[2], iz);
if (TMath::Abs(point[2]-fZ[iz])<1.e-8 ||
TMath::Abs(fZ[iz+1]-point[2])<1.e-8) xtru->SetIz(-1);
return fPoly->Contains(point);
}
Int_t TGeoXtru::DistancetoPrimitive(Int_t px, Int_t py)
{
const Int_t numPoints = fNvert*fNz;
return ShapeDistancetoPrimitive(numPoints, px, py);
}
Double_t TGeoXtru::DistToPlane(Double_t *point, Double_t *dir, Int_t iz, Int_t ivert, Double_t stepmax, Bool_t in) const
{
Double_t snext;
Double_t vert[12];
Double_t norm[3];
Double_t znew;
Double_t pt[3];
Double_t safe;
if (fZ[iz]==fZ[iz+1] && !in) {
TGeoXtru *xtru = (TGeoXtru*)this;
snext = (fZ[iz]-point[2])/dir[2];
pt[0] = point[0]+snext*dir[0];
pt[1] = point[1]+snext*dir[1];
pt[2] = point[2]+snext*dir[2];
xtru->SetCurrentVertices(fX0[iz], fY0[iz], fScale[iz]);
if (!xtru->Contains(pt)) return TGeoShape::Big();
xtru->SetCurrentVertices(fX0[iz+1], fY0[iz+1], fScale[iz+1]);
if (!xtru->Contains(pt)) return TGeoShape::Big();
return snext;
}
GetPlaneVertices(iz, ivert, vert);
GetPlaneNormal(vert, norm);
Double_t ndotd = norm[0]*dir[0]+norm[1]*dir[1]+norm[2]*dir[2];
if (in) {
if (ndotd<=0) return TGeoShape::Big();
safe = (vert[0]-point[0])*norm[0]+
(vert[1]-point[1])*norm[1]+
(vert[2]-point[2])*norm[2];
if (safe<0) return TGeoShape::Big();
} else {
ndotd = -ndotd;
if (ndotd<=0) return TGeoShape::Big();
safe = (point[0]-vert[0])*norm[0]+
(point[1]-vert[1])*norm[1]+
(point[2]-vert[2])*norm[2];
if (safe<0) return TGeoShape::Big();
}
snext = safe/ndotd;
if (snext>stepmax) return TGeoShape::Big();
if (fZ[iz]<fZ[iz+1]) {
znew = point[2] + snext*dir[2];
if (znew<fZ[iz]) return TGeoShape::Big();
if (znew>fZ[iz+1]) return TGeoShape::Big();
}
pt[0] = point[0]+snext*dir[0];
pt[1] = point[1]+snext*dir[1];
pt[2] = point[2]+snext*dir[2];
if (!IsPointInsidePlane(pt, vert, norm)) return TGeoShape::Big();
return snext;
}
Double_t TGeoXtru::DistFromInside(Double_t *point, Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
{
if (iact<3 && safe) {
*safe = Safety(point, kTRUE);
if (iact==0) return TGeoShape::Big();
if (iact==1 && step<*safe) return TGeoShape::Big();
}
TGeoXtru *xtru = (TGeoXtru*)this;
Int_t iz = TMath::BinarySearch(fNz, fZ, point[2]);
if (iz==fNz-1) {
if (dir[2]>=0) {
xtru->SetIz(-1);
return 0.;
}
iz--;
} else {
if (iz>0) {
if (point[2] == fZ[iz]) {
if ((fZ[iz]==fZ[iz+1]) && (dir[2]<0)) iz++;
else if ((fZ[iz]==fZ[iz-1]) && (dir[2]>0)) iz--;
}
}
}
Bool_t convex = fPoly->IsConvex();
Double_t stepmax = step;
Double_t snext = TGeoShape::Big();
Double_t dist, sz;
Double_t pt[3];
Int_t iv, ipl, inext;
if (dir[2]==0) {
for (iv=0; iv<fNvert; iv++) {
xtru->SetIz(-1);
dist = DistToPlane(point,dir,iz,iv,stepmax,kTRUE);
if (dist<stepmax) {
stepmax = dist;
snext = dist;
xtru->SetSeg(iv);
if (convex) return snext;
}
}
return snext;
}
Int_t incseg = (dir[2]>0)?1:-1;
Int_t iznext = iz;
Bool_t zexit = kFALSE;
while (iz>=0 && iz<fNz-1) {
ipl = iz+((incseg+1)>>1);
inext = ipl+incseg;
sz = (fZ[ipl]-point[2])/dir[2];
if (sz<stepmax) {
iznext += incseg;
pt[0] = point[0]+sz*dir[0];
pt[1] = point[1]+sz*dir[1];
xtru->SetCurrentVertices(fX0[ipl],fY0[ipl],fScale[ipl]);
if (fPoly->Contains(pt)) {
if (ipl==0 || ipl==fNz-1) {
xtru->SetIz(-1);
if (convex) return sz;
zexit = kTRUE;
snext = sz;
stepmax = sz;
}
if (!zexit && fZ[ipl]==fZ[inext]) {
xtru->SetCurrentVertices(fX0[inext],fY0[inext],fScale[inext]);
if (!fPoly->Contains(pt)) {
xtru->SetIz(-1);
if (convex) return sz;
zexit = kTRUE;
snext = sz;
stepmax = sz;
} else {
iznext = inext;
}
}
}
} else {
iznext = fNz-1;
}
xtru->SetIz(iz);
for (iv=0; iv<fNvert; iv++) {
dist = DistToPlane(point,dir,iz,iv,stepmax,kTRUE);
if (dist<stepmax) {
xtru->SetSeg(iv);
snext = dist;
stepmax = dist;
if (convex) return snext;
zexit = kTRUE;
}
}
if (zexit) return snext;
iz = iznext;
}
return TGeoShape::Big();
}
Double_t TGeoXtru::DistFromOutside(Double_t *point, Double_t *dir, Int_t iact, Double_t step, Double_t *safe) const
{
if (iact<3 && safe) {
*safe = Safety(point, kTRUE);
if (iact==0) return TGeoShape::Big();
if (iact==1 && step<*safe) return TGeoShape::Big();
}
Double_t stepmax = step;
Double_t snext = 0.;
Double_t dist = TGeoShape::Big();
Int_t i, iv;
Double_t pt[3];
memcpy(pt,point,3*sizeof(Double_t));
TGeoXtru *xtru = (TGeoXtru*)this;
Int_t iz = TMath::BinarySearch(fNz, fZ, point[2]);
if (iz<0) {
if (dir[2]<=0) return TGeoShape::Big();
snext = (fZ[0] - point[2])/dir[2];
if (snext>stepmax) return TGeoShape::Big();
for (i=0; i<3; i++) pt[i] = point[i] + snext*dir[i];
xtru->SetCurrentVertices(fX0[0],fY0[0],fScale[0]);
if (fPoly->Contains(pt)) {
xtru->SetIz(-1);
return snext;
}
iz=0;
stepmax -= snext;
} else {
if (iz==fNz-1) {
if (dir[2]>=0) return TGeoShape::Big();
snext = (fZ[fNz-1] - point[2])/dir[2];
if (snext>stepmax) return TGeoShape::Big();
for (i=0; i<3; i++) pt[i] = point[i] + snext*dir[i];
xtru->SetCurrentVertices(fX0[fNz-1],fY0[fNz-1],fScale[fNz-1]);
if (fPoly->Contains(pt)) {
xtru->SetIz(-1);
return snext;
}
iz = fNz-2;
stepmax -= snext;
}
}
if (!TGeoBBox::Contains(pt)) {
dist = TGeoBBox::DistFromOutside(pt,dir,3);
if (dist>stepmax) return TGeoShape::Big();
if (dist>1E-6) dist-=1E-6;
for (i=0; i<3; i++) pt[i] += dist*dir[i];
iz = TMath::BinarySearch(fNz, fZ, pt[2]);
if (iz<0) iz=0;
else if (iz==fNz-1) iz = fNz-2;
snext += dist;
stepmax -= dist;
}
Bool_t convex = fPoly->IsConvex();
Bool_t hit = kFALSE;
if (dir[2]==0) {
xtru->SetIz(iz);
for (iv=0; iv<fNvert; iv++) {
dist = DistToPlane(pt,dir,iz,iv,stepmax,kFALSE);
if (dist<stepmax) {
xtru->SetSeg(iv);
if (convex) return (snext+dist);
stepmax = dist;
hit = kTRUE;
}
}
if (hit) return (snext+stepmax);
return TGeoShape::Big();
}
Int_t incseg = (dir[2]>0)?1:-1;
while (iz>=0 && iz<fNz-1) {
xtru->SetIz(iz);
if (fZ[iz]==fZ[iz+1]) xtru->SetIz(-1);
for (iv=0; iv<fNvert; iv++) {
dist = DistToPlane(pt,dir,iz,iv,stepmax,kFALSE);
if (dist<stepmax) {
xtru->SetSeg(iv);
if (convex) return (snext+dist);
stepmax = dist;
hit = kTRUE;
}
}
if (hit) return (snext+stepmax);
iz += incseg;
}
return TGeoShape::Big();
}
Bool_t TGeoXtru::DefinePolygon(Int_t nvert, const Double_t *xv, const Double_t *yv)
{
if (nvert<3) {
Error("DefinePolygon","In shape %s cannot create polygon with less than 3 vertices", GetName());
return kFALSE;
}
fNvert = nvert;
if (fX) delete [] fX;
fX = new Double_t[nvert];
if (fY) delete [] fY;
fY = new Double_t[nvert];
if (fXc) delete [] fXc;
fXc = new Double_t[nvert];
if (fYc) delete [] fYc;
fYc = new Double_t[nvert];
memcpy(fX,xv,nvert*sizeof(Double_t));
memcpy(fXc,xv,nvert*sizeof(Double_t));
memcpy(fY,yv,nvert*sizeof(Double_t));
memcpy(fYc,yv,nvert*sizeof(Double_t));
if (fPoly) delete fPoly;
fPoly = new TGeoPolygon(nvert);
fPoly->SetXY(fXc,fYc);
fPoly->FinishPolygon();
return kTRUE;
}
void TGeoXtru::DefineSection(Int_t snum, Double_t z, Double_t x0, Double_t y0, Double_t scale)
{
if ((snum<0) || (snum>=fNz)) return;
fZ[snum] = z;
fX0[snum] = x0;
fY0[snum] = y0;
fScale[snum] = scale;
if (snum) {
if (fZ[snum]<fZ[snum-1]) {
Warning("DefineSection", "In shape: %s, Z position of section "
"%i, z=%e, not in increasing order, %i, z=%e",
GetName(),snum,fZ[snum],snum-1,fZ[snum-1]);
return;
}
}
if (snum==(fNz-1)) ComputeBBox();
}
Double_t TGeoXtru::GetZ(Int_t ipl) const
{
if (ipl<0 || ipl>(fNz-1)) {
Error("GetZ","In shape %s, ipl=%i out of range (0,%i)",GetName(),ipl,0,fNz-1);
return 0.;
}
return fZ[ipl];
}
void TGeoXtru::GetPlaneNormal(const Double_t *vert, Double_t *norm) const
{
Double_t cross = 0.;
Double_t v1[3], v2[3];
v1[0] = vert[9]-vert[0];
v1[1] = vert[10]-vert[1];
v1[2] = vert[11]-vert[2];
v2[0] = vert[3]-vert[0];
v2[1] = vert[4]-vert[1];
v2[2] = vert[5]-vert[2];
norm[0] = v1[1]*v2[2]-v1[2]*v2[1];
cross += norm[0]*norm[0];
norm[1] = v1[2]*v2[0]-v1[0]*v2[2];
cross += norm[1]*norm[1];
norm[2] = v1[0]*v2[1]-v1[1]*v2[0];
cross += norm[2]*norm[2];
cross = 1./TMath::Sqrt(cross);
for (Int_t i=0; i<3; i++) norm[i] *= cross;
}
void TGeoXtru::GetPlaneVertices(Int_t iz, Int_t ivert, Double_t *vert) const
{
Double_t x,y,z1,z2;
Int_t iv1 = (ivert+1)%fNvert;
Int_t icrt = 0;
z1 = fZ[iz];
z2 = fZ[iz+1];
if (fPoly->IsClockwise()) {
x = fX[ivert]*fScale[iz]+fX0[iz];
y = fY[ivert]*fScale[iz]+fY0[iz];
vert[icrt++] = x;
vert[icrt++] = y;
vert[icrt++] = z1;
x = fX[iv1]*fScale[iz]+fX0[iz];
y = fY[iv1]*fScale[iz]+fY0[iz];
vert[icrt++] = x;
vert[icrt++] = y;
vert[icrt++] = z1;
x = fX[iv1]*fScale[iz+1]+fX0[iz+1];
y = fY[iv1]*fScale[iz+1]+fY0[iz+1];
vert[icrt++] = x;
vert[icrt++] = y;
vert[icrt++] = z2;
x = fX[ivert]*fScale[iz+1]+fX0[iz+1];
y = fY[ivert]*fScale[iz+1]+fY0[iz+1];
vert[icrt++] = x;
vert[icrt++] = y;
vert[icrt++] = z2;
} else {
x = fX[iv1]*fScale[iz]+fX0[iz];
y = fY[iv1]*fScale[iz]+fY0[iz];
vert[icrt++] = x;
vert[icrt++] = y;
vert[icrt++] = z1;
x = fX[ivert]*fScale[iz]+fX0[iz];
y = fY[ivert]*fScale[iz]+fY0[iz];
vert[icrt++] = x;
vert[icrt++] = y;
vert[icrt++] = z1;
x = fX[ivert]*fScale[iz+1]+fX0[iz+1];
y = fY[ivert]*fScale[iz+1]+fY0[iz+1];
vert[icrt++] = x;
vert[icrt++] = y;
vert[icrt++] = z2;
x = fX[iv1]*fScale[iz+1]+fX0[iz+1];
y = fY[iv1]*fScale[iz+1]+fY0[iz+1];
vert[icrt++] = x;
vert[icrt++] = y;
vert[icrt++] = z2;
}
}
Bool_t TGeoXtru::IsPointInsidePlane(Double_t *point, Double_t *vert, Double_t *norm) const
{
Double_t v1[3], v2[3];
Double_t cross;
Int_t j,k;
for (Int_t i=0; i<4; i++) {
j = 3*i;
k = 3*((i+1)%4);
v1[0] = point[0]-vert[j];
v1[1] = point[1]-vert[j+1];
v1[2] = point[2]-vert[j+2];
v2[0] = vert[k]-vert[j];
v2[1] = vert[k+1]-vert[j+1];
v2[2] = vert[k+2]-vert[j+2];
cross = (v1[1]*v2[2]-v1[2]*v2[1])*norm[0]+
(v1[2]*v2[0]-v1[0]*v2[2])*norm[1]+
(v1[0]*v2[1]-v1[1]*v2[0])*norm[2];
if (cross<0) return kFALSE;
}
return kTRUE;
}
void TGeoXtru::InspectShape() const
{
printf("*** Shape %s: TGeoXtru ***\n", GetName());
printf(" Nz = %i\n", fNz);
printf(" List of (x,y) of polygon vertices:\n");
for (Int_t ivert = 0; ivert<fNvert; ivert++)
printf(" x = %11.5f y = %11.5f\n", fX[ivert],fY[ivert]);
for (Int_t ipl=0; ipl<fNz; ipl++)
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]);
printf(" Bounding box:\n");
TGeoBBox::InspectShape();
}
TBuffer3D *TGeoXtru::MakeBuffer3D() const
{
Int_t nz = GetNz();
Int_t nvert = GetNvert();
Int_t nbPnts = nz*nvert;
Int_t nbSegs = nvert*(2*nz-1);
Int_t nbPols = nvert*(nz-1)+2;
TBuffer3D* buff = new TBuffer3D(TBuffer3DTypes::kGeneric,
nbPnts, 3*nbPnts, nbSegs, 3*nbSegs, nbPols, 6*(nbPols-2)+2*(2+nvert));
if (buff)
{
SetPoints(buff->fPnts);
SetSegsAndPols(*buff);
}
return buff;
}
void TGeoXtru::SetSegsAndPols(TBuffer3D &buff) const
{
Int_t nz = GetNz();
Int_t nvert = GetNvert();
Int_t c = GetBasicColor();
Int_t i,j;
Int_t indx, indx2, k;
indx = indx2 = 0;
for (i=0; i<nz; i++) {
indx2 = i*nvert;
for (j=0; j<nvert; j++) {
k = (j+1)%nvert;
buff.fSegs[indx++] = c;
buff.fSegs[indx++] = indx2+j;
buff.fSegs[indx++] = indx2+k;
}
}
for (i=0; i<nz-1; i++) {
indx2 = i*nvert;
for (j=0; j<nvert; j++) {
k = j + nvert;
buff.fSegs[indx++] = c;
buff.fSegs[indx++] = indx2+j;
buff.fSegs[indx++] = indx2+k;
}
}
indx = 0;
for (i=0; i<nz-1; i++) {
indx2 = i*nvert;
for (j=0; j<nvert; j++) {
k = (j+1)%nvert;
buff.fPols[indx++] = c+j%3;
buff.fPols[indx++] = 4;
buff.fPols[indx++] = indx2+j;
buff.fPols[indx++] = nz*nvert+indx2+k;
buff.fPols[indx++] = indx2+nvert+j;
buff.fPols[indx++] = nz*nvert+indx2+j;
}
}
buff.fPols[indx++] = c+2;
buff.fPols[indx++] = nvert;
indx2 = 0;
for (j = nvert - 1; j >= 0; --j) {
buff.fPols[indx++] = indx2+j;
}
buff.fPols[indx++] = c;
buff.fPols[indx++] = nvert;
indx2 = (nz-1)*nvert;
for (j=0; j<nvert; j++) {
buff.fPols[indx++] = indx2+j;
}
}
Double_t TGeoXtru::SafetyToSector(Double_t *point, Int_t iz, Double_t safmin)
{
Double_t safz = TGeoShape::Big();
Double_t saf1, saf2;
Bool_t in1, in2;
Int_t iseg;
Double_t safe = TGeoShape::Big();
if (fZ[iz] == fZ[iz+1]) {
safz = TMath::Abs(point[2]-fZ[iz]);
if (safz>safmin) return TGeoShape::Big();
SetCurrentVertices(fX0[iz], fY0[iz], fScale[iz]);
saf1 = fPoly->Safety(point, iseg);
in1 = fPoly->Contains(point);
if (!in1 && saf1>safmin) return TGeoShape::Big();
SetCurrentVertices(fX0[iz+1], fY0[iz+1], fScale[iz+1]);
saf2 = fPoly->Safety(point, iseg);
in2 = fPoly->Contains(point);
if ((in1&!in2)|(in2&!in1)) {
safe = safz;
} else {
safe = TMath::Min(saf1,saf2);
safe = TMath::Max(safe, safz);
}
if (safe>safmin) return TGeoShape::Big();
return safe;
}
safz = fZ[iz]-point[2];
if (safz>safmin) return TGeoShape::Big();
if (safz<0) {
saf1 = point[2]-fZ[iz+1];
if (saf1>safmin) return TGeoShape::Big();
if (saf1<0) {
safz = 0.;
} else {
safz = saf1;
}
}
SetCurrentZ(point[2],iz);
saf1 = fPoly->Safety(point, iseg);
Double_t vert[12];
Double_t norm[3];
GetPlaneVertices(iz,iseg,vert);
GetPlaneNormal(vert, norm);
saf1 = saf1*TMath::Sqrt(1.-norm[2]*norm[2]);
safe = TMath::Max(safz, saf1);
if (safe>safmin) return TGeoShape::Big();
return safe;
}
Double_t TGeoXtru::Safety(Double_t *point, Bool_t in) const
{
Double_t safmin = TGeoShape::Big();
Double_t safe;
Double_t safz = 0.;
TGeoXtru *xtru = (TGeoXtru*)this;
Int_t iz;
if (in) {
safmin = TMath::Min(point[2]-fZ[0], fZ[fNz-1]-point[2]);
for (iz=0; iz<fNz-1; iz++) {
safe = xtru->SafetyToSector(point, iz, safmin);
if (safe<safmin) safmin = safe;
}
return safmin;
}
iz = TMath::BinarySearch(fNz, fZ, point[2]);
if (iz<0) {
iz = 0;
safz = fZ[0] - point[2];
} else {
if (iz==fNz-1) {
iz = fNz-2;
safz = point[2] - fZ[fNz-1];
}
}
Int_t i;
for (i=iz; i<fNz-1; i++) {
safe = xtru->SafetyToSector(point,i,safmin);
if (safe<safmin) safmin=safe;
}
for (i=iz-1; i>0; i--) {
safe = xtru->SafetyToSector(point,i,safmin);
if (safe<safmin) safmin=safe;
}
safe = TMath::Max(safmin, safz);
return safe;
}
void TGeoXtru::SavePrimitive(ostream &out, Option_t * )
{
if (TObject::TestBit(kGeoSavePrimitive)) return;
out << " // Shape: " << GetName() << " type: " << ClassName() << endl;
out << " nz = " << fNz << ";" << endl;
out << " nvert = " << fNvert << ";" << endl;
out << " TGeoXtru *xtru = new TGeoXtru(nz);" << endl;
out << " xtru->SetName(\"" << GetName() << "\");" << endl;
Int_t i;
for (i=0; i<fNvert; i++) {
out << " xvert[" << i << "] = " << fX[i] << "; yvert[" << i << "] = " << fY[i] << ";" << endl;
}
out << " xtru->DefinePolygon(nvert,xvert,yvert);" << endl;
for (i=0; i<fNz; i++) {
out << " zsect = " << fZ[i] << ";" << endl;
out << " x0 = " << fX0[i] << ";" << endl;
out << " y0 = " << fY0[i] << ";" << endl;
out << " scale0 = " << fScale[i] << ";" << endl;
out << " xtru->DefineSection(" << i << ",zsect,x0,y0,scale0);" << endl;
}
out << " TGeoShape *" << GetPointerName() << " = xtru;" << endl;
TObject::SetBit(TGeoShape::kGeoSavePrimitive);
}
void TGeoXtru::SetCurrentZ(Double_t z, Int_t iz)
{
Double_t x0, y0, scale, a, b;
Int_t ind1, ind2;
ind1 = iz;
ind2 = iz+1;
Double_t invdz = 1./(fZ[ind2]-fZ[ind1]);
a = (fX0[ind1]*fZ[ind2]-fX0[ind2]*fZ[ind1])*invdz;
b = (fX0[ind2]-fX0[ind1])*invdz;
x0 = a+b*z;
a = (fY0[ind1]*fZ[ind2]-fY0[ind2]*fZ[ind1])*invdz;
b = (fY0[ind2]-fY0[ind1])*invdz;
y0 = a+b*z;
a = (fScale[ind1]*fZ[ind2]-fScale[ind2]*fZ[ind1])*invdz;
b = (fScale[ind2]-fScale[ind1])*invdz;
scale = a+b*z;
SetCurrentVertices(x0,y0,scale);
}
void TGeoXtru::SetCurrentVertices(Double_t x0, Double_t y0, Double_t scale)
{
for (Int_t i=0; i<fNvert; i++) {
fXc[i] = scale*fX[i] + x0;
fYc[i] = scale*fY[i] + y0;
}
}
void TGeoXtru::SetDimensions(Double_t *param)
{
fNz = (Int_t)param[0];
if (fNz<2) {
Error("SetDimensions","Cannot create TGeoXtru %s with less than 2 Z planes",GetName());
return;
}
if (fZ) delete [] fZ;
if (fScale) delete [] fScale;
if (fX0) delete [] fX0;
if (fY0) delete [] fY0;
fZ = new Double_t[fNz];
fScale = new Double_t[fNz];
fX0 = new Double_t[fNz];
fY0 = new Double_t[fNz];
for (Int_t i=0; i<fNz; i++)
DefineSection(i, param[1+4*i], param[2+4*i], param[3+4*i], param[4+4*i]);
}
void TGeoXtru::SetPoints(Double_t *points) const
{
Int_t i, j;
Int_t indx = 0;
TGeoXtru *xtru = (TGeoXtru*)this;
if (points) {
for (i = 0; i < fNz; i++) {
xtru->SetCurrentVertices(fX0[i], fY0[i], fScale[i]);
if (fPoly->IsClockwise()) {
for (j = 0; j < fNvert; j++) {
points[indx++] = fXc[j];
points[indx++] = fYc[j];
points[indx++] = fZ[i];
}
} else {
for (j = 0; j < fNvert; j++) {
points[indx++] = fXc[fNvert-1-j];
points[indx++] = fYc[fNvert-1-j];
points[indx++] = fZ[i];
}
}
}
}
}
void TGeoXtru::SetPoints(Float_t *points) const
{
Int_t i, j;
Int_t indx = 0;
TGeoXtru *xtru = (TGeoXtru*)this;
if (points) {
for (i = 0; i < fNz; i++) {
xtru->SetCurrentVertices(fX0[i], fY0[i], fScale[i]);
if (fPoly->IsClockwise()) {
for (j = 0; j < fNvert; j++) {
points[indx++] = fXc[j];
points[indx++] = fYc[j];
points[indx++] = fZ[i];
}
} else {
for (j = 0; j < fNvert; j++) {
points[indx++] = fXc[fNvert-1-j];
points[indx++] = fYc[fNvert-1-j];
points[indx++] = fZ[i];
}
}
}
}
}
Int_t TGeoXtru::GetNmeshVertices() const
{
Int_t numPoints = fNz*fNvert;
return numPoints;
}
void TGeoXtru::Sizeof3D() const
{
}
void TGeoXtru::Streamer(TBuffer &R__b)
{
if (R__b.IsReading()) {
TGeoXtru::Class()->ReadBuffer(R__b, this);
if (fPoly) fPoly->SetXY(fXc,fYc);
} else {
TGeoXtru::Class()->WriteBuffer(R__b, this);
}
}
const TBuffer3D & TGeoXtru::GetBuffer3D(Int_t reqSections, Bool_t localFrame) const
{
static TBuffer3D buffer(TBuffer3DTypes::kGeneric);
TGeoBBox::FillBuffer3D(buffer, reqSections, localFrame);
if (reqSections & TBuffer3D::kRawSizes) {
Int_t nz = GetNz();
Int_t nvert = GetNvert();
Int_t nbPnts = nz*nvert;
Int_t nbSegs = nvert*(2*nz-1);
Int_t nbPols = nvert*(nz-1)+2;
if (buffer.SetRawSizes(nbPnts, 3*nbPnts, nbSegs, 3*nbSegs, nbPols, 6*(nbPols-2)+2*(2+nvert))) {
buffer.SetSectionsValid(TBuffer3D::kRawSizes);
}
}
if ((reqSections & TBuffer3D::kRaw) && buffer.SectionsValid(TBuffer3D::kRawSizes)) {
SetPoints(buffer.fPnts);
if (!buffer.fLocalFrame) {
TransformPoints(buffer.fPnts, buffer.NbPnts());
}
SetSegsAndPols(buffer);
buffer.SetSectionsValid(TBuffer3D::kRaw);
}
return buffer;
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.