// @(#)root/geom:$Name: $:$Id: TGeoArb8.cxx,v 1.44 2005/09/04 15:12:08 brun Exp $
// Author: Andrei Gheata 31/01/02
/*************************************************************************
* Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
* All rights reserved. *
* *
* For the licensing terms see $ROOTSYS/LICENSE. *
* For the list of contributors see $ROOTSYS/README/CREDITS. *
*************************************************************************/
#include "Riostream.h"
#include "TROOT.h"
#include "TGeoManager.h"
#include "TGeoVolume.h"
#include "TGeoArb8.h"
ClassImp(TGeoArb8)
//________________________________________________________________________
// TGeoArb8 - a arbitrary trapezoid with less than 8 vertices standing on
// two paralel planes perpendicular to Z axis. Parameters :
// - dz - half length in Z;
// - xy[8][2] - vector of (x,y) coordinates of vertices
// - first four points (xy[i][j], i<4, j<2) are the (x,y)
// coordinates of the vertices sitting on the -dz plane;
// - last four points (xy[i][j], i>=4, j<2) are the (x,y)
// coordinates of the vertices sitting on the +dz plane;
// The order of defining the vertices of an arb8 is the following :
// - point 0 is connected with points 1,3,4
// - point 1 is connected with points 0,2,5
// - point 2 is connected with points 1,3,6
// - point 3 is connected with points 0,2,7
// - point 4 is connected with points 0,5,7
// - point 5 is connected with points 1,4,6
// - point 6 is connected with points 2,5,7
// - point 7 is connected with points 3,4,6
// Points can be identical in order to create shapes with less than
// 8 vertices.
//
//
/*
*/
//
////////////////////////////////////////////////////////////////////////////
// //
// TGeoTrap //
// //
// TRAP is a general trapezoid, i.e. one for which the faces perpendicular//
// to z are trapezia and their centres are not the same x, y. It has 11 //
// parameters: the half length in z, the polar angles from the centre of //
// the face at low z to that at high z, H1 the half length in y at low z, //
// LB1 the half length in x at low z and y low edge, LB2 the half length //
// in x at low z and y high edge, TH1 the angle w.r.t. the y axis from the//
// centre of low y edge to the centre of the high y edge, and H2, LB2, //
// LH2, TH2, the corresponding quantities at high z. //
// //
////////////////////////////////////////////////////////////////////////////
//
/*
*/
//
//
//
/*
*/
//
////////////////////////////////////////////////////////////////////////////
// //
// TGeoGtra //
// //
// Gtra is a twisted trapezoid, i.e. one for which the faces perpendicular//
// to z are trapezia and their centres are not the same x, y. It has 12 //
// parameters: the half length in z, the polar angles from the centre of //
// the face at low z to that at high z, twist, H1 the half length in y at low z, //
// LB1 the half length in x at low z and y low edge, LB2 the half length //
// in x at low z and y high edge, TH1 the angle w.r.t. the y axis from the//
// centre of low y edge to the centre of the high y edge, and H2, LB2, //
// LH2, TH2, the corresponding quantities at high z. //
// //
////////////////////////////////////////////////////////////////////////////
//
/*
*/
//
//
//
/*
*/
//
//_____________________________________________________________________________
TGeoArb8::TGeoArb8()
{
// dummy ctor
fDz = 0;
fTwist = 0;
for (Int_t i=0; i<8; i++) {
fXY[i][0] = 0.0;
fXY[i][1] = 0.0;
}
SetShapeBit(kGeoArb8);
}
//_____________________________________________________________________________
TGeoArb8::TGeoArb8(Double_t dz, Double_t *vertices)
:TGeoBBox(0,0,0)
{
// constructor. If the array of vertices is not null, this should be
// in the format : (x0, y0, x1, y1, ... , x7, y7)
fDz = dz;
fTwist = 0;
SetShapeBit(kGeoArb8);
if (vertices) {
for (Int_t i=0; i<8; i++) {
fXY[i][0] = vertices[2*i];
fXY[i][1] = vertices[2*i+1];
}
ComputeTwist();
ComputeBBox();
} else {
for (Int_t i=0; i<8; i++) {
fXY[i][0] = 0.0;
fXY[i][1] = 0.0;
}
}
}
//_____________________________________________________________________________
TGeoArb8::TGeoArb8(const char *name, Double_t dz, Double_t *vertices)
:TGeoBBox(name, 0,0,0)
{
// constructor. If the array of vertices is not null, this should be
// in the format : (x0, y0, x1, y1, ... , x7, y7)
fDz = dz;
fTwist = 0;
SetShapeBit(kGeoArb8);
if (vertices) {
for (Int_t i=0; i<8; i++) {
fXY[i][0] = vertices[2*i];
fXY[i][1] = vertices[2*i+1];
}
ComputeTwist();
ComputeBBox();
} else {
for (Int_t i=0; i<8; i++) {
fXY[i][0] = 0.0;
fXY[i][1] = 0.0;
}
}
}
//_____________________________________________________________________________
TGeoArb8::~TGeoArb8()
{
// destructor
if (fTwist) delete [] fTwist;
}
//_____________________________________________________________________________
void TGeoArb8::ComputeBBox()
{
// compute bounding box for a Arb8
Double_t xmin, xmax, ymin, ymax;
xmin = xmax = fXY[0][0];
ymin = ymax = fXY[0][1];
for (Int_t i=1; i<8; i++) {
if (xmin>fXY[i][0]) xmin=fXY[i][0];
if (xmax<fXY[i][0]) xmax=fXY[i][0];
if (ymin>fXY[i][1]) ymin=fXY[i][1];
if (ymax<fXY[i][1]) ymax=fXY[i][1];
}
fDX = 0.5*(xmax-xmin);
fDY = 0.5*(ymax-ymin);
fDZ = fDz;
fOrigin[0] = 0.5*(xmax+xmin);
fOrigin[1] = 0.5*(ymax+ymin);
fOrigin[2] = 0;
SetShapeBit(kGeoClosedShape);
}
//_____________________________________________________________________________
void TGeoArb8::ComputeTwist()
{
// compute tangents of twist angles (angles between projections on XY plane
// of corresponding -dz +dz edges). Called after last point [7] was set.
Double_t twist[4];
Bool_t twisted = kFALSE;
Double_t dx1, dy1, dx2, dy2;
for (Int_t i=0; i<4; i++) {
dx1 = fXY[(i+1)%4][0]-fXY[i][0];
dy1 = fXY[(i+1)%4][1]-fXY[i][1];
if (dx1==0 && dy1==0) {
twist[i] = 0;
continue;
}
dx2 = fXY[4+(i+1)%4][0]-fXY[4+i][0];
dy2 = fXY[4+(i+1)%4][1]-fXY[4+i][1];
if (dx2==0 && dy2==0) {
twist[i] = 0;
continue;
}
twist[i] = dy1*dx2 - dx1*dy2;
if (TMath::Abs(twist[i])<1E-3) {
twist[i] = 0;
continue;
}
twist[i] = TMath::Sign(1.,twist[i]);
twisted = kTRUE;
}
if (!twisted) return;
if (fTwist) delete [] fTwist;
fTwist = new Double_t[4];
memcpy(fTwist, &twist[0], 4*sizeof(Double_t));
}
//_____________________________________________________________________________
Double_t TGeoArb8::GetTwist(Int_t iseg) const
{
// Get twist for segment I in range [0,3]
if (!fTwist) return 0.;
if (iseg<0 || iseg>3) return 0.;
return fTwist[iseg];
}
//_____________________________________________________________________________
void TGeoArb8::ComputeNormal(Double_t *point, Double_t *dir, Double_t *norm)
{
// Compute normal to closest surface from POINT.
Double_t safe = TGeoShape::Big();
Double_t safc;
Int_t i; // current facette index
Double_t x0, y0, z0, x1, y1, z1, x2, y2;
Double_t ax, ay, az, bx, by;
Double_t fn;
safc = fDz-TMath::Abs(point[2]);
if (safc<1E-4) {
memset(norm,0,3*sizeof(Double_t));
norm[2] = (dir[2]>0)?1:(-1);
return;
}
Double_t vert[8], lnorm[3];
SetPlaneVertices(point[2], vert);
//---> compute safety for lateral planes
for (i=0; i<4; i++) {
x0 = vert[2*i];
y0 = vert[2*i+1];
z0 = point[2];
x1 = fXY[i+4][0];
y1 = fXY[i+4][1];
z1 = fDz;
ax = x1-x0;
ay = y1-y0;
az = z1-z0;
x2 = vert[2*((i+1)%4)];
y2 = vert[2*((i+1)%4)+1];
bx = x2-x0;
by = y2-y0;
lnorm[0] = -az*by;
lnorm[1] = az*bx;
lnorm[2] = ax*by-ay*bx;
fn = TMath::Sqrt(lnorm[0]*lnorm[0]+lnorm[1]*lnorm[1]+lnorm[2]*lnorm[2]);
if (fn<1E-10) continue;
lnorm[0] /= fn;
lnorm[1] /= fn;
lnorm[2] /= fn;
safc = (x0-point[0])*lnorm[0]+(y0-point[1])*lnorm[1]+(z0-point[2])*lnorm[2];
safc = TMath::Abs(safc);
// printf("plane %i : (%g, %g, %g) safe=%g\n", i, lnorm[0],lnorm[1],lnorm[2],safc);
if (safc<safe) {
safe = safc;
memcpy(norm,lnorm,3*sizeof(Double_t));
}
}
if (dir[0]*norm[0]+dir[1]*norm[1]+dir[2]*norm[2] < 0) {
norm[0] = -norm[0];
norm[1] = -norm[1];
norm[2] = -norm[2];
}
}
//_____________________________________________________________________________
Bool_t TGeoArb8::Contains(Double_t *point) const
{
// test if point is inside this sphere
// first check Z range
if (TMath::Abs(point[2]) > fDz) return kFALSE;
// compute intersection between Z plane containing point and the arb8
Double_t poly[8];
// memset(&poly[0], 0, 8*sizeof(Double_t));
//SetPlaneVertices(point[2], &poly[0]);
Double_t cf = 0.5*(fDz-point[2])/fDz;
Int_t i;
for (i=0; i<4; i++) {
poly[2*i] = fXY[i+4][0]+cf*(fXY[i][0]-fXY[i+4][0]);
poly[2*i+1] = fXY[i+4][1]+cf*(fXY[i][1]-fXY[i+4][1]);
}
return InsidePolygon(point[0],point[1],poly);
}
//_____________________________________________________________________________
Double_t TGeoArb8::DistToPlane(Double_t *point, Double_t *dir, Int_t ipl, Bool_t in) const
{
// compute distance to plane ipl :
// ipl=0 : points 0,4,1,5
// ipl=1 : points 1,5,2,6
// ipl=2 : points 2,6,3,7
// ipl=3 : points 3,7,0,4
Double_t xa,xb,xc,xd;
Double_t ya,yb,yc,yd;
Int_t j = (ipl+1)%4;
xa=fXY[ipl][0];
ya=fXY[ipl][1];
xb=fXY[ipl+4][0];
yb=fXY[ipl+4][1];
xc=fXY[j][0];
yc=fXY[j][1];
xd=fXY[4+j][0];
yd=fXY[4+j][1];
Double_t dz2 =0.5/fDz;
Double_t tx1 =dz2*(xb-xa);
Double_t ty1 =dz2*(yb-ya);
Double_t tx2 =dz2*(xd-xc);
Double_t ty2 =dz2*(yd-yc);
Double_t dzp =fDz+point[2];
Double_t xs1 =xa+tx1*dzp;
Double_t ys1 =ya+ty1*dzp;
Double_t xs2 =xc+tx2*dzp;
Double_t ys2 =yc+ty2*dzp;
Double_t dxs =xs2-xs1;
Double_t dys =ys2-ys1;
Double_t dtx =tx2-tx1;
Double_t dty =ty2-ty1;
Double_t a=(dtx*dir[1]-dty*dir[0]+(tx1*ty2-tx2*ty1)*dir[2])*dir[2];
Double_t b=dxs*dir[1]-dys*dir[0]+(dtx*point[1]-dty*point[0]+ty2*xs1-ty1*xs2
+tx1*ys2-tx2*ys1)*dir[2];
Double_t c=dxs*point[1]-dys*point[0]+xs1*ys2-xs2*ys1;
Double_t s=TGeoShape::Big();
Double_t x1,x2,y1,y2,xp,yp,zi;
if (TMath::Abs(a)<1E-10) {
if (b==0) return TGeoShape::Big();
s=-c/b;
if (s>0) {
if (in) return s;
zi=point[2]+s*dir[2];
if (TMath::Abs(zi)<fDz) {
x1=xs1+tx1*dir[2]*s;
x2=xs2+tx2*dir[2]*s;
xp=point[0]+s*dir[0];
y1=ys1+ty1*dir[2]*s;
y2=ys2+ty2*dir[2]*s;
yp=point[1]+s*dir[1];
zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
if (zi<=0) return s;
}
}
return TGeoShape::Big();
}
b=0.5*b/a;
c=c/a;
Double_t d=b*b-c;
if (d>=0) {
Double_t sqd = TMath::Sqrt(d);
s=-b-sqd;
if (s>0) {
if (in) return s;
zi=point[2]+s*dir[2];
if (TMath::Abs(zi)<fDz) {
x1=xs1+tx1*dir[2]*s;
x2=xs2+tx2*dir[2]*s;
xp=point[0]+s*dir[0];
y1=ys1+ty1*dir[2]*s;
y2=ys2+ty2*dir[2]*s;
yp=point[1]+s*dir[1];
zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
if (zi<=0) return s;
}
}
s=-b+sqd;
if (s>0) {
if (in) return s;
zi=point[2]+s*dir[2];
if (TMath::Abs(zi)<fDz) {
x1=xs1+tx1*dir[2]*s;
x2=xs2+tx2*dir[2]*s;
xp=point[0]+s*dir[0];
y1=ys1+ty1*dir[2]*s;
y2=ys2+ty2*dir[2]*s;
yp=point[1]+s*dir[1];
zi = (xp-x1)*(xp-x2)+(yp-y1)*(yp-y2);
if (zi<=0) return s;
}
}
}
return TGeoShape::Big();
}
//_____________________________________________________________________________
Double_t TGeoArb8::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 arb8
Double_t snxt=TGeoShape::Big();
if (!TGeoBBox::Contains(point)) {
snxt=TGeoBBox::DistFromOutside(point,dir,3);
if (snxt>1E20) return snxt;
}
Double_t dist[5];
// check lateral faces
Int_t i;
for (i=0; i<4; i++) {
dist[i]=DistToPlane(point, dir, i, kFALSE);
}
// check Z planes
dist[4]=TGeoShape::Big();
if (TMath::Abs(point[2])>fDz) {
if (dir[2]!=0) {
Double_t pt[3];
if (point[2]>0) {
dist[4] = (fDz-point[2])/dir[2];
pt[2]=fDz;
} else {
dist[4] = (-fDz-point[2])/dir[2];
pt[2]=-fDz;
}
if (dist[4]<0) {
dist[4]=TGeoShape::Big();
} else {
for (Int_t j=0; j<2; j++) pt[j]=point[j]+dist[4]*dir[j];
if (!Contains(&pt[0])) dist[4]=TGeoShape::Big();
}
}
}
Double_t distmin = dist[0];
for (i=1;i<5;i++) if (dist[i] < distmin) distmin = dist[i];
return distmin;
}
//_____________________________________________________________________________
Double_t TGeoArb8::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 arb8
#ifdef OLDALGORITHM
Int_t i;
Double_t dist[6];
dist[0]=dist[1]=TGeoShape::Big();
if (dir[2]<0) {
dist[0]=(-fDz-point[2])/dir[2];
} else {
if (dir[2]>0) dist[1]=(fDz-point[2])/dir[2];
}
for (i=0; i<4; i++) {
dist[i+2]=DistToPlane(point, dir, i, kTRUE);
}
Double_t distmin = dist[0];
for (i=1;i<6;i++) if (dist[i] < distmin) distmin = dist[i];
return distmin;
#else
// compute distance to plane ipl :
// ipl=0 : points 0,4,1,5
// ipl=1 : points 1,5,2,6
// ipl=2 : points 2,6,3,7
// ipl=3 : points 3,7,0,4
Double_t distmin;
if (dir[2]<0) {
distmin=(-fDz-point[2])/dir[2];
} else {
if (dir[2]>0) distmin =(fDz-point[2])/dir[2];
else distmin = TGeoShape::Big();
}
Double_t dz2 =0.5/fDz;
Double_t xa,xb,xc,xd;
Double_t ya,yb,yc,yd;
for (Int_t ipl=0;ipl<4;ipl++) {
Int_t j = (ipl+1)%4;
xa=fXY[ipl][0];
ya=fXY[ipl][1];
xb=fXY[ipl+4][0];
yb=fXY[ipl+4][1];
xc=fXY[j][0];
yc=fXY[j][1];
xd=fXY[4+j][0];
yd=fXY[4+j][1];
Double_t tx1 =dz2*(xb-xa);
Double_t ty1 =dz2*(yb-ya);
Double_t tx2 =dz2*(xd-xc);
Double_t ty2 =dz2*(yd-yc);
Double_t dzp =fDz+point[2];
Double_t xs1 =xa+tx1*dzp;
Double_t ys1 =ya+ty1*dzp;
Double_t xs2 =xc+tx2*dzp;
Double_t ys2 =yc+ty2*dzp;
Double_t dxs =xs2-xs1;
Double_t dys =ys2-ys1;
Double_t dtx =tx2-tx1;
Double_t dty =ty2-ty1;
Double_t a=(dtx*dir[1]-dty*dir[0]+(tx1*ty2-tx2*ty1)*dir[2])*dir[2];
Double_t b=dxs*dir[1]-dys*dir[0]+(dtx*point[1]-dty*point[0]+ty2*xs1-ty1*xs2
+tx1*ys2-tx2*ys1)*dir[2];
Double_t c=dxs*point[1]-dys*point[0]+xs1*ys2-xs2*ys1;
Double_t s=TGeoShape::Big();
if (TMath::Abs(a)<1E-10) {
if (b==0) continue;
s=-c/b;
if (s>0 && s < distmin) distmin =s;
continue;
}
b=0.5*b/a;
c=c/a;
Double_t d=b*b-c;
if (d>=0) {
Double_t sqd = TMath::Sqrt(d);
s=-b-sqd;
if (s>0) {
if (s < distmin) distmin = s;
} else {
s=-b+sqd;
if (s>0 && s < distmin) distmin =s;
}
}
}
return distmin;
#endif
}
//_____________________________________________________________________________
TGeoVolume *TGeoArb8::Divide(TGeoVolume *voldiv, const char * /*divname*/, Int_t /*iaxis*/, Int_t /*ndiv*/,
Double_t /*start*/, Double_t /*step*/)
{
Error("Divide", "Division of an arbitrary trapezoid not implemented");
return voldiv;
}
//_____________________________________________________________________________
Double_t TGeoArb8::GetAxisRange(Int_t iaxis, Double_t &xlo, Double_t &xhi) const
{
xlo = 0;
xhi = 0;
Double_t dx = 0;
if (iaxis==3) {
xlo = -fDz;
xhi = fDz;
dx = xhi-xlo;
return dx;
}
return dx;
}
//_____________________________________________________________________________
void TGeoArb8::GetBoundingCylinder(Double_t *param) const
{
//--- Fill vector param[4] with the bounding cylinder parameters. The order
// is the following : Rmin, Rmax, Phi1, Phi2
//--- first compute rmin/rmax
Double_t rmaxsq = 0;
Double_t rsq;
Int_t i;
for (i=0; i<8; i++) {
rsq = fXY[i][0]*fXY[i][0] + fXY[i][1]*fXY[i][1];
rmaxsq = TMath::Max(rsq, rmaxsq);
}
param[0] = 0.; // Rmin
param[1] = rmaxsq; // Rmax
param[2] = 0.; // Phi1
param[3] = 360.; // Phi2
}
//_____________________________________________________________________________
Int_t TGeoArb8::GetFittingBox(const TGeoBBox *parambox, TGeoMatrix *mat, Double_t &dx, Double_t &dy, Double_t &dz) const
{
// Fills real parameters of a positioned box inside this arb8. Returns 0 if successfull.
dx=dy=dz=0;
if (mat->IsRotation()) {
Error("GetFittingBox", "cannot handle parametrized rotated volumes");
return 1; // ### rotation not accepted ###
}
//--> translate the origin of the parametrized box to the frame of this box.
Double_t origin[3];
mat->LocalToMaster(parambox->GetOrigin(), origin);
if (!Contains(origin)) {
Error("GetFittingBox", "wrong matrix - parametrized box is outside this");
return 1; // ### wrong matrix ###
}
//--> now we have to get the valid range for all parametrized axis
Double_t dd[3];
dd[0] = parambox->GetDX();
dd[1] = parambox->GetDY();
dd[2] = parambox->GetDZ();
//-> check if Z range is fixed
if (dd[2]<0) {
dd[2] = TMath::Min(origin[2]+fDz, fDz-origin[2]);
if (dd[2]<0) {
Error("GetFittingBox", "wrong matrix");
return 1;
}
}
if (dd[0]>=0 && dd[1]>=0) {
dx = dd[0];
dy = dd[1];
dz = dd[2];
return 0;
}
//-> check now vertices at Z = origin[2] +/- dd[2]
Double_t upper[8];
Double_t lower[8];
SetPlaneVertices(origin[2]-dd[2], lower);
SetPlaneVertices(origin[2]+dd[2], upper);
Double_t ddmin=TGeoShape::Big();
for (Int_t iaxis=0; iaxis<2; iaxis++) {
if (dd[iaxis]>=0) continue;
ddmin=TGeoShape::Big();
for (Int_t ivert=0; ivert<4; ivert++) {
ddmin = TMath::Min(ddmin, TMath::Abs(origin[iaxis]-lower[2*ivert+iaxis]));
ddmin = TMath::Min(ddmin, TMath::Abs(origin[iaxis]-upper[2*ivert+iaxis]));
}
dd[iaxis] = ddmin;
}
dx = dd[0];
dy = dd[1];
dz = dd[2];
return 0;
}
//_____________________________________________________________________________
void TGeoArb8::GetPlaneNormal(Double_t *p1, Double_t *p2, Double_t *p3, Double_t *norm)
{
// Compute normal to plane defined by P1, P2 and P3
Double_t cross = 0.;
Double_t v1[3], v2[3];
Int_t i;
for (i=0; i<3; i++) {
v1[i] = p2[i] - p1[i];
v2[i] = p3[i] - p1[i];
}
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];
if (cross == 0.) return;
cross = 1./TMath::Sqrt(cross);
for (i=0; i<3; i++) norm[i] *= cross;
}
//_____________________________________________________________________________
Bool_t TGeoArb8::InsidePolygon(Double_t x, Double_t y, Double_t *pts)
{
// Find if a point in XY plane is inside the polygon defines by PTS.
Int_t i,j;
Double_t x1,y1,x2,y2;
Double_t cross;
for (i=0; i<4; i++) {
j = (i+1)%4;
x1 = pts[i<<1];
y1 = pts[(i<<1)+1];
x2 = pts[j<<1];
y2 = pts[(j<<1)+1];
cross = (x-x1)*(y2-y1)-(y-y1)*(x2-x1);
if (cross<0) return kFALSE;
}
return kTRUE;
}
//_____________________________________________________________________________
void TGeoArb8::InspectShape() const
{
// print shape parameters
printf("*** Shape %s: TGeoArb8 ***\n", GetName());
if (IsTwisted()) printf(" = TWISTED\n");
for (Int_t ip=0; ip<8; ip++) {
printf(" point #%i : x=%11.5f y=%11.5f z=%11.5f\n",
ip, fXY[ip][0], fXY[ip][1], fDz*((ip<4)?-1:1));
}
printf(" Bounding box:\n");
TGeoBBox::InspectShape();
}
//_____________________________________________________________________________
Double_t TGeoArb8::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 safz = fDz-TMath::Abs(point[2]);
if (!in) safz = -safz;
Int_t iseg;
Double_t safmin = TGeoShape::Big();
Double_t safe = TGeoShape::Big();
Double_t lsq, ssq, dx, dy, dpx, dpy, u;
if (IsTwisted()) {
if (!in) {
if (!TGeoBBox::Contains(point)) return TGeoBBox::Safety(point,kFALSE);
}
// Point is also in the bounding box ;-(
// Compute closest distance to any segment
Double_t vert[8];
Double_t *p1, *p2;
Int_t isegmin=0;
Double_t umin = 0.;
SetPlaneVertices (point[2], vert);
for (iseg=0; iseg<4; iseg++) {
if (safe==0.) return 0.;
p1 = &vert[2*iseg];
p2 = &vert[2*((iseg+1)%4)];
dx = p2[0] - p1[0];
dy = p2[1] - p1[1];
dpx = point[0] - p1[0];
dpy = point[1] - p1[1];
lsq = dx*dx + dy*dy;
u = (dpx*dx + dpy*dy)/lsq;
if (u>1) {
dpx = point[0]-p2[0];
dpy = point[1]-p2[1];
} else {
if (u>=0) {
dpx -= u*dx;
dpy -= u*dy;
}
}
ssq = dpx*dpx + dpy*dpy;
if (ssq < safe) {
isegmin = iseg;
umin = u;
safe = ssq;
}
}
if (umin<0) umin = 0.;
if (umin>1) {
isegmin = (isegmin+1)%4;
umin = 0.;
}
Int_t i1 = isegmin;
Int_t i2 = (isegmin+1)%4;
Double_t dx1 = fXY[i2][0]-fXY[i1][0];
Double_t dx2 = fXY[i2+4][0]-fXY[i1+4][0];
Double_t dy1 = fXY[i2][1]-fXY[i1][1];
Double_t dy2 = fXY[i2+4][1]-fXY[i1+4][1];
dx = dx1 + umin*(dx2-dx1);
dy = dy1 + umin*(dy2-dy1);
safe *= 1.- 4.*fDz*fDz/(dx*dx+dy*dy+4.*fDz*fDz);
safe = TMath::Sqrt(safe);
return safe;
}
for (iseg=0; iseg<4; iseg++) {
safe = SafetyToFace(point,iseg,in);
if (safe>0) {
if (in && safe<safmin) {
safmin = safe;
continue;
}
if (!in && safe<1E10) {
if (safmin<1E10) safe = TMath::Max(safe,safmin);
else safmin=safe;
}
}
}
if (in) return TMath::Min(safmin, safz);
return TMath::Max(safmin, safz);
}
//_____________________________________________________________________________
Double_t TGeoArb8::SafetyToFace(Double_t *point, Int_t iseg, Bool_t in) const
{
// Estimate safety to lateral plane defined by segment iseg in range [0,3]
// might be negative: plane seen only from inside
Double_t vertices[12];
Int_t ipln = (iseg+1)%4;
// point 1
vertices[0] = fXY[iseg][0];
vertices[1] = fXY[iseg][1];
vertices[2] = -fDz;
// point 2
vertices[3] = fXY[ipln][0];
vertices[4] = fXY[ipln][1];
vertices[5] = -fDz;
// point 3
vertices[6] = fXY[ipln+4][0];
vertices[7] = fXY[ipln+4][1];
vertices[8] = fDz;
// point 4
vertices[9] = fXY[iseg+4][0];
vertices[10] = fXY[iseg+4][1];
vertices[11] = fDz;
Double_t twist = GetTwist(iseg);
Double_t safe;
Double_t norm[3];
Double_t *p1, *p2, *p3;
if (twist ==0) {
p1 = &vertices[0];
p2 = &vertices[9];
p3 = &vertices[6];
if (IsSamePoint(p2,p3)) {
p3 = &vertices[3];
if (IsSamePoint(p1,p3)) return TGeoShape::Big(); // skip single segment
}
GetPlaneNormal(p1,p2,p3,norm);
safe = (point[0]-p1[0])*norm[0]+(point[1]-p1[1])*norm[1]+(point[2]-p1[2])*norm[2];
if (in) return (-safe);
return safe;
}
// The face is twisted
return TGeoShape::Big();
}
//_____________________________________________________________________________
void TGeoArb8::SavePrimitive(ofstream &out, Option_t * /*option*/)
{
// Save a primitive as a C++ statement(s) on output stream "out".
if (TObject::TestBit(kGeoSavePrimitive)) return;
out << " // Shape: " << GetName() << " type: " << ClassName() << endl;
out << " dz = " << fDz << ";" << endl;
out << " vert[0] = " << fXY[0][0] << ";" << endl;
out << " vert[1] = " << fXY[0][1] << ";" << endl;
out << " vert[2] = " << fXY[1][0] << ";" << endl;
out << " vert[3] = " << fXY[1][1] << ";" << endl;
out << " vert[4] = " << fXY[2][0] << ";" << endl;
out << " vert[5] = " << fXY[2][1] << ";" << endl;
out << " vert[6] = " << fXY[3][0] << ";" << endl;
out << " vert[7] = " << fXY[3][1] << ";" << endl;
out << " vert[8] = " << fXY[4][0] << ";" << endl;
out << " vert[9] = " << fXY[4][1] << ";" << endl;
out << " vert[10] = " << fXY[5][0] << ";" << endl;
out << " vert[11] = " << fXY[5][1] << ";" << endl;
out << " vert[12] = " << fXY[6][0] << ";" << endl;
out << " vert[13] = " << fXY[6][1] << ";" << endl;
out << " vert[14] = " << fXY[7][0] << ";" << endl;
out << " vert[15] = " << fXY[7][1] << ";" << endl;
out << " TGeoShape *" << GetPointerName() << " = new TGeoArb8(\"" << GetName() << "\", dz,vert);" << endl;
TObject::SetBit(TGeoShape::kGeoSavePrimitive);
}
//_____________________________________________________________________________
void TGeoArb8::SetPlaneVertices(Double_t zpl, Double_t *vertices) const
{
// compute intersection points between plane at zpl and non-horizontal edges
Double_t cf = 0.5*(fDz-zpl)/fDz;
for (Int_t i=0; i<4; i++) {
vertices[2*i] = fXY[i+4][0]+cf*(fXY[i][0]-fXY[i+4][0]);
vertices[2*i+1] = fXY[i+4][1]+cf*(fXY[i][1]-fXY[i+4][1]);
}
}
//_____________________________________________________________________________
void TGeoArb8::SetDimensions(Double_t *param)
{
// set arb8 params in one step :
fDz = param[0];
for (Int_t i=0; i<8; i++) {
fXY[i][0] = param[2*i];
fXY[i][1] = param[2*i+1];
}
ComputeTwist();
ComputeBBox();
}
//_____________________________________________________________________________
void TGeoArb8::SetPoints(Double_t *points) const
{
// create arb8 mesh points
for (Int_t i=0; i<8; i++) {
points[3*i] = fXY[i][0];
points[3*i+1] = fXY[i][1];
points[3*i+2] = (i<4)?-fDz:fDz;
}
}
//_____________________________________________________________________________
void TGeoArb8::SetPoints(Float_t *points) const
{
// create arb8 mesh points
for (Int_t i=0; i<8; i++) {
points[3*i] = fXY[i][0];
points[3*i+1] = fXY[i][1];
points[3*i+2] = (i<4)?-fDz:fDz;
}
}
//_____________________________________________________________________________
void TGeoArb8::SetVertex(Int_t vnum, Double_t x, Double_t y)
{
// set values for a given vertex
if (vnum<0 || vnum >7) {
Error("SetVertex", "Invalid vertex number");
return;
}
fXY[vnum][0] = x;
fXY[vnum][1] = y;
if (vnum == 7) {
ComputeTwist();
ComputeBBox();
}
}
//_____________________________________________________________________________
void TGeoArb8::Sizeof3D() const
{
// fill size of this 3-D object
TGeoBBox::Sizeof3D();
}
ClassImp(TGeoTrap)
//_____________________________________________________________________________
TGeoTrap::TGeoTrap()
{
// dummy ctor
fDz = 0;
fTheta = 0;
fPhi = 0;
fH1 = fH2 = fBl1 = fBl2 = fTl1 = fTl2 = fAlpha1 = fAlpha2 = 0;
}
//_____________________________________________________________________________
TGeoTrap::TGeoTrap(Double_t dz, Double_t theta, Double_t phi)
:TGeoArb8("", 0, 0)
{
fDz = dz;
fTheta = theta;
fPhi = phi;
fH1 = fH2 = fBl1 = fBl2 = fTl1 = fTl2 = fAlpha1 = fAlpha2 = 0;
}
//_____________________________________________________________________________
TGeoTrap::TGeoTrap(Double_t dz, Double_t theta, Double_t phi, Double_t h1,
Double_t bl1, Double_t tl1, Double_t alpha1, Double_t h2, Double_t bl2,
Double_t tl2, Double_t alpha2)
:TGeoArb8("", 0, 0)
{
// constructor.
fDz = dz;
fTheta = theta;
fPhi = phi;
fH1 = h1;
fH2 = h2;
fBl1 = bl1;
fBl2 = bl2;
fTl1 = tl1;
fTl2 = tl2;
fAlpha1 = alpha1;
fAlpha2 = alpha2;
Double_t tx = TMath::Tan(theta*TMath::DegToRad())*TMath::Cos(phi*TMath::DegToRad());
Double_t ty = TMath::Tan(theta*TMath::DegToRad())*TMath::Sin(phi*TMath::DegToRad());
Double_t ta1 = TMath::Tan(alpha1*TMath::DegToRad());
Double_t ta2 = TMath::Tan(alpha2*TMath::DegToRad());
fXY[0][0] = -dz*tx-h1*ta1-bl1; fXY[0][1] = -dz*ty-h1;
fXY[1][0] = -dz*tx+h1*ta1-tl1; fXY[1][1] = -dz*ty+h1;
fXY[2][0] = -dz*tx+h1*ta1+tl1; fXY[2][1] = -dz*ty+h1;
fXY[3][0] = -dz*tx-h1*ta1+bl1; fXY[3][1] = -dz*ty-h1;
fXY[4][0] = dz*tx-h2*ta2-bl2; fXY[4][1] = dz*ty-h2;
fXY[5][0] = dz*tx+h2*ta2-tl2; fXY[5][1] = dz*ty+h2;
fXY[6][0] = dz*tx+h2*ta2+tl2; fXY[6][1] = dz*ty+h2;
fXY[7][0] = dz*tx-h2*ta2+bl2; fXY[7][1] = dz*ty-h2;
ComputeTwist();
if ((dz<0) || (h1<0) || (bl1<0) || (tl1<0) ||
(h2<0) || (bl2<0) || (tl2<0)) {
SetShapeBit(kGeoRunTimeShape);
}
else TGeoArb8::ComputeBBox();
}
//_____________________________________________________________________________
TGeoTrap::TGeoTrap(const char *name, Double_t dz, Double_t theta, Double_t phi, Double_t h1,
Double_t bl1, Double_t tl1, Double_t alpha1, Double_t h2, Double_t bl2,
Double_t tl2, Double_t alpha2)
:TGeoArb8(name, 0, 0)
{
// constructor with name
SetName(name);
fDz = dz;
fTheta = theta;
fPhi = phi;
fH1 = h1;
fH2 = h2;
fBl1 = bl1;
fBl2 = bl2;
fTl1 = tl1;
fTl2 = tl2;
fAlpha1 = alpha1;
fAlpha2 = alpha2;
for (Int_t i=0; i<8; i++) {
fXY[i][0] = 0.0;
fXY[i][1] = 0.0;
}
Double_t tx = TMath::Tan(theta*TMath::DegToRad())*TMath::Cos(phi*TMath::DegToRad());
Double_t ty = TMath::Tan(theta*TMath::DegToRad())*TMath::Sin(phi*TMath::DegToRad());
Double_t ta1 = TMath::Tan(alpha1*TMath::DegToRad());
Double_t ta2 = TMath::Tan(alpha2*TMath::DegToRad());
fXY[0][0] = -dz*tx-h1*ta1-bl1; fXY[0][1] = -dz*ty-h1;
fXY[1][0] = -dz*tx+h1*ta1-tl1; fXY[1][1] = -dz*ty+h1;
fXY[2][0] = -dz*tx+h1*ta1+tl1; fXY[2][1] = -dz*ty+h1;
fXY[3][0] = -dz*tx-h1*ta1+bl1; fXY[3][1] = -dz*ty-h1;
fXY[4][0] = dz*tx-h2*ta2-bl2; fXY[4][1] = dz*ty-h2;
fXY[5][0] = dz*tx+h2*ta2-tl2; fXY[5][1] = dz*ty+h2;
fXY[6][0] = dz*tx+h2*ta2+tl2; fXY[6][1] = dz*ty+h2;
fXY[7][0] = dz*tx-h2*ta2+bl2; fXY[7][1] = dz*ty-h2;
ComputeTwist();
if ((dz<0) || (h1<0) || (bl1<0) || (tl1<0) ||
(h2<0) || (bl2<0) || (tl2<0)) {
SetShapeBit(kGeoRunTimeShape);
}
else TGeoArb8::ComputeBBox();
}
//_____________________________________________________________________________
TGeoTrap::~TGeoTrap()
{
// destructor
}
//_____________________________________________________________________________
Double_t TGeoTrap::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 arb8
if (iact<3 && safe) {
// compute safe distance
*safe = Safety(point, kTRUE);
if (iact==0) return TGeoShape::Big();
if (iact==1 && step<*safe) return TGeoShape::Big();
}
// compute distance to get ouside this shape
// return TGeoArb8::DistFromInside(point, dir, iact, step, safe);
// compute distance to plane ipl :
// ipl=0 : points 0,4,1,5
// ipl=1 : points 1,5,2,6
// ipl=2 : points 2,6,3,7
// ipl=3 : points 3,7,0,4
Double_t distmin;
if (dir[2]<0) {
distmin=(-fDz-point[2])/dir[2];
} else {
if (dir[2]>0) distmin =(fDz-point[2])/dir[2];
else distmin = TGeoShape::Big();
}
Double_t xa,xb,xc;
Double_t ya,yb,yc;
for (Int_t ipl=0;ipl<4;ipl++) {
Int_t j = (ipl+1)%4;
xa=fXY[ipl][0];
ya=fXY[ipl][1];
xb=fXY[ipl+4][0];
yb=fXY[ipl+4][1];
xc=fXY[j][0];
yc=fXY[j][1];
Double_t ax,ay,az;
ax = xb-xa;
ay = yb-ya;
az = 2.*fDz;
Double_t bx,by;
bx = xc-xa;
by = yc-ya;
Double_t ddotn = -dir[0]*az*by + dir[1]*az*bx+dir[2]*(ax*by-ay*bx);
if (ddotn<=0) continue; // entering
Double_t saf = -(point[0]-xa)*az*by + (point[1]-ya)*az*bx + (point[2]+fDz)*(ax*by-ay*bx);
if (saf>=0.0) return 0.0;
Double_t s = -saf/ddotn;
if (s<distmin) distmin=s;
}
return distmin;
}
//_____________________________________________________________________________
Double_t TGeoTrap::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 arb8
if (iact<3 && safe) {
// compute safe distance
*safe = Safety(point, kFALSE);
if (iact==0) return TGeoShape::Big();
if (iact==1 && step<*safe) return TGeoShape::Big();
}
// compute distance to get ouside this shape
Bool_t in = kTRUE;
Double_t pts[8];
Double_t snxt;
Double_t xnew, ynew, znew;
Int_t i,j;
if (point[2]<-fDz+TGeoShape::Tolerance()) {
if (dir[2]<=0) return TGeoShape::Big();
in = kFALSE;
snxt = -(fDz+point[2])/dir[2];
xnew = point[0] + snxt*dir[0];
ynew = point[1] + snxt*dir[1];
for (i=0;i<4;i++) {
j = i<<1;
pts[j] = fXY[i][0];
pts[j+1] = fXY[i][1];
}
if (InsidePolygon(xnew,ynew,pts)) return snxt;
} else if (point[2]>fDz-TGeoShape::Tolerance()) {
if (dir[2]>=0) return TGeoShape::Big();
in = kFALSE;
snxt = (fDz-point[2])/dir[2];
xnew = point[0] + snxt*dir[0];
ynew = point[1] + snxt*dir[1];
for (i=0;i<4;i++) {
j = i<<1;
pts[j] = fXY[i+4][0];
pts[j+1] = fXY[i+4][1];
}
if (InsidePolygon(xnew,ynew,pts)) return snxt;
}
snxt = TGeoShape::Big();
// check lateral faces
Double_t dz2 =0.5/fDz;
Double_t xa,xb,xc,xd;
Double_t ya,yb,yc,yd;
Double_t ax,ay,az;
Double_t bx,by;
Double_t ddotn, saf;
Double_t safmin = TGeoShape::Big();
Bool_t exiting = kFALSE;
for (i=0; i<4; i++) {
j = (i+1)%4;
xa=fXY[i][0];
ya=fXY[i][1];
xb=fXY[i+4][0];
yb=fXY[i+4][1];
xc=fXY[j][0];
yc=fXY[j][1];
xd=fXY[4+j][0];
yd=fXY[4+j][1];
ax = xb-xa;
ay = yb-ya;
az = 2.*fDz;
bx = xc-xa;
by = yc-ya;
ddotn = -dir[0]*az*by + dir[1]*az*bx+dir[2]*(ax*by-ay*bx);
saf = (point[0]-xa)*az*by - (point[1]-ya)*az*bx - (point[2]+fDz)*(ax*by-ay*bx);
if (saf<=0) {
// face visible from point outside
in = kFALSE;
if (ddotn>=0) return TGeoShape::Big();
snxt = saf/ddotn;
znew = point[2]+snxt*dir[2];
if (TMath::Abs(znew)<=fDz) {
xnew = point[0]+snxt*dir[0];
ynew = point[1]+snxt*dir[1];
Double_t tx1 =dz2*(xb-xa);
Double_t ty1 =dz2*(yb-ya);
Double_t tx2 =dz2*(xd-xc);
Double_t ty2 =dz2*(yd-yc);
Double_t dzp =fDz+znew;
Double_t xs1 =xa+tx1*dzp;
Double_t ys1 =ya+ty1*dzp;
Double_t xs2 =xc+tx2*dzp;
Double_t ys2 =yc+ty2*dzp;
if (TMath::Abs(xs1-xs2)>TMath::Abs(ys1-ys2)) {
if ((xnew-xs1)*(xs2-xnew)>=0) return snxt;
} else {
if ((ynew-ys1)*(ys2-ynew)>=0) return snxt;
}
}
} else {
if (saf<safmin) {
safmin = saf;
if (ddotn>=0) exiting = kTRUE;
else exiting = kFALSE;
}
}
}
if (!in) return TGeoShape::Big();
if (exiting) return TGeoShape::Big();
return 0.0;
}
//_____________________________________________________________________________
TGeoVolume *TGeoTrap::Divide(TGeoVolume *voldiv, const char *divname, Int_t iaxis, Int_t ndiv,
Double_t start, Double_t step)
{
//--- Divide this trapezoid shape belonging to volume "voldiv" into ndiv volumes
// called divname, from start position with the given step. Only Z divisions
// are supported. For Z divisions just return the pointer to the volume to be
// divided. In case a wrong division axis is supplied, returns pointer to
// volume that was divided.
TGeoShape *shape; //--- shape to be created
TGeoVolume *vol; //--- division volume to be created
TGeoVolumeMulti *vmulti; //--- generic divided volume
TGeoPatternFinder *finder; //--- finder to be attached
TString opt = ""; //--- option to be attached
if (iaxis!=3) {
Error("Divide", "cannot divide trapezoids on other axis than Z");
return 0;
}
Double_t end = start+ndiv*step;
Double_t points_lo[8];
Double_t points_hi[8];
finder = new TGeoPatternTrapZ(voldiv, ndiv, start, end);
voldiv->SetFinder(finder);
finder->SetDivIndex(voldiv->GetNdaughters());
opt = "Z";
vmulti = gGeoManager->MakeVolumeMulti(divname, voldiv->GetMedium());
Double_t txz = ((TGeoPatternTrapZ*)finder)->GetTxz();
Double_t tyz = ((TGeoPatternTrapZ*)finder)->GetTyz();
Double_t zmin, zmax, ox,oy,oz;
for (Int_t idiv=0; idiv<ndiv; idiv++) {
zmin = start+idiv*step;
zmax = start+(idiv+1)*step;
oz = start+idiv*step+step/2;
ox = oz*txz;
oy = oz*tyz;
SetPlaneVertices(zmin, &points_lo[0]);
SetPlaneVertices(zmax, &points_hi[0]);
shape = new TGeoTrap(step/2, fTheta, fPhi);
for (Int_t vert1=0; vert1<4; vert1++)
((TGeoArb8*)shape)->SetVertex(vert1, points_lo[2*vert1]-ox, points_lo[2*vert1+1]-oy);
for (Int_t vert2=0; vert2<4; vert2++)
((TGeoArb8*)shape)->SetVertex(vert2+4, points_hi[2*vert2]-ox, points_hi[2*vert2+1]-oy);
vol = new TGeoVolume(divname, shape, voldiv->GetMedium());
vmulti->AddVolume(vol);
voldiv->AddNodeOffset(vol, idiv, oz, opt.Data());
((TGeoNodeOffset*)voldiv->GetNodes()->At(voldiv->GetNdaughters()-1))->SetFinder(finder);
}
return vmulti;
}
//_____________________________________________________________________________
TGeoShape *TGeoTrap::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;
if (mother->IsRunTimeShape()) {
Error("GetMakeRuntimeShape", "invalid mother");
return 0;
}
Double_t dz, h1, bl1, tl1, h2, bl2, tl2;
if (fDz<0) dz=((TGeoTrap*)mother)->GetDz();
else dz=fDz;
if (fH1<0)
h1 = ((TGeoTrap*)mother)->GetH1();
else
h1 = fH1;
if (fH2<0)
h2 = ((TGeoTrap*)mother)->GetH2();
else
h2 = fH2;
if (fBl1<0)
bl1 = ((TGeoTrap*)mother)->GetBl1();
else
bl1 = fBl1;
if (fBl2<0)
bl2 = ((TGeoTrap*)mother)->GetBl2();
else
bl2 = fBl2;
if (fTl1<0)
tl1 = ((TGeoTrap*)mother)->GetTl1();
else
tl1 = fTl1;
if (fTl2<0)
tl2 = ((TGeoTrap*)mother)->GetTl2();
else
tl2 = fTl2;
return (new TGeoTrap(dz, fTheta, fPhi, h1, bl1, tl1, fAlpha1, h2, bl2, tl2, fAlpha2));
}
//_____________________________________________________________________________
Double_t TGeoTrap::Safety(Double_t *point, Bool_t in) const
{
// Computes the closest distance from given point to this shape, according
// to option.
Double_t safe = TGeoShape::Big();
Double_t saf[5];
Double_t norm[3]; // normal to current facette
Int_t i,j; // current facette index
Double_t x0, y0, z0=-fDz, x1, y1, z1=fDz, x2, y2;
Double_t ax, ay, az=z1-z0, bx, by;
Double_t fn;
//---> compute safety for lateral planes
for (i=0; i<4; i++) {
x0 = fXY[i][0];
y0 = fXY[i][1];
x1 = fXY[i+4][0];
y1 = fXY[i+4][1];
ax = x1-x0;
ay = y1-y0;
az = z1-z0;
j = (i+1)%4;
x2 = fXY[j][0];
y2 = fXY[j][1];
bx = x2-x0;
by = y2-y0;
if (bx==0 && by==0) {
x2 = fXY[4+j][0];
y2 = fXY[4+j][1];
bx = x2-x1;
by = y2-y1;
if (bx==0 && by==0) continue;
}
norm[0] = -az*by;
norm[1] = az*bx;
norm[2] = ax*by-ay*bx;
fn = TMath::Sqrt(norm[0]*norm[0]+norm[1]*norm[1]+norm[2]*norm[2]);
if (fn<1E-10) continue;
saf[i] = (x0-point[0])*norm[0]+(y0-point[1])*norm[1]+(-fDz-point[2])*norm[2];
if (in) {
saf[i]=TMath::Abs(saf[i])/fn; // they should be all positive anyway
} else {
saf[i] = -saf[i]/fn; // only negative values are interesting
}
}
saf[4] = fDz-TMath::Abs(point[2]);
if (in) {
safe = saf[0];
for (j=1;j<5;j++) if (saf[j] <safe) safe = saf[j];
} else {
saf[4]=-saf[4];
safe = saf[0];
for (j=1;j<5;j++) if (saf[j] >safe) safe = saf[j];
}
return safe;
}
//_____________________________________________________________________________
void TGeoTrap::SavePrimitive(ofstream &out, Option_t * /*option*/)
{
// Save a primitive as a C++ statement(s) on output stream "out".
if (TObject::TestBit(kGeoSavePrimitive)) return;
out << " // Shape: " << GetName() << " type: " << ClassName() << endl;
out << " dz = " << fDz << ";" << endl;
out << " theta = " << fTheta << ";" << endl;
out << " phi = " << fPhi << ";" << endl;
out << " h1 = " << fH1<< ";" << endl;
out << " bl1 = " << fBl1<< ";" << endl;
out << " tl1 = " << fTl1<< ";" << endl;
out << " alpha1 = " << fAlpha1 << ";" << endl;
out << " h2 = " << fH2 << ";" << endl;
out << " bl2 = " << fBl2<< ";" << endl;
out << " tl2 = " << fTl2<< ";" << endl;
out << " alpha2 = " << fAlpha2 << ";" << endl;
out << " TGeoShape *" << GetPointerName() << " = new TGeoTrap(\"" << GetName() << "\", dz,theta,phi,h1,bl1,tl1,alpha1,h2,bl2,tl2,alpha2);" << endl;
TObject::SetBit(TGeoShape::kGeoSavePrimitive);
}
ClassImp(TGeoGtra)
//_____________________________________________________________________________
TGeoGtra::TGeoGtra()
{
// dummy ctor
fTwistAngle = 0;
}
//_____________________________________________________________________________
TGeoGtra::TGeoGtra(Double_t dz, Double_t theta, Double_t phi, Double_t twist, Double_t h1,
Double_t bl1, Double_t tl1, Double_t alpha1, Double_t h2, Double_t bl2,
Double_t tl2, Double_t alpha2)
:TGeoTrap(dz, theta, phi, h1, bl1, tl1, alpha1, h2, bl2, tl2, alpha2)
{
// constructor.
fTheta = theta;
fPhi = phi;
fH1 = h1;
fH2 = h2;
fBl1 = bl1;
fBl2 = bl2;
fTl1 = tl1;
fTl2 = tl2;
fAlpha1 = alpha1;
fAlpha2 = alpha2;
Double_t x, y, dx, dy, dx1, dx2, th, ph, al1, al2;
al1 = alpha1*TMath::DegToRad();
al2 = alpha2*TMath::DegToRad();
th = theta*TMath::DegToRad();
ph = phi*TMath::DegToRad();
dx = 2*dz*TMath::Sin(th)*TMath::Cos(ph);
dy = 2*dz*TMath::Sin(th)*TMath::Sin(ph);
fDz = dz;
dx1 = 2*h1*TMath::Tan(al1);
dx2 = 2*h2*TMath::Tan(al2);
fTwistAngle = twist;
Int_t i;
for (i=0; i<8; i++) {
fXY[i][0] = 0.0;
fXY[i][1] = 0.0;
}
fXY[0][0] = -bl1; fXY[0][1] = -h1;
fXY[1][0] = -tl1+dx1; fXY[1][1] = h1;
fXY[2][0] = tl1+dx1; fXY[2][1] = h1;
fXY[3][0] = bl1; fXY[3][1] = -h1;
fXY[4][0] = -bl2+dx; fXY[4][1] = -h2+dy;
fXY[5][0] = -tl2+dx+dx2; fXY[5][1] = h2+dy;
fXY[6][0] = tl2+dx+dx2; fXY[6][1] = h2+dy;
fXY[7][0] = bl2+dx; fXY[7][1] = -h2+dy;
for (i=4; i<8; i++) {
x = fXY[i][0];
y = fXY[i][1];
fXY[i][0] = x*TMath::Cos(twist*TMath::DegToRad()) + y*TMath::Sin(twist*TMath::DegToRad());
fXY[i][1] = -x*TMath::Sin(twist*TMath::DegToRad()) + y*TMath::Cos(twist*TMath::DegToRad());
}
ComputeTwist();
if ((dz<0) || (h1<0) || (bl1<0) || (tl1<0) ||
(h2<0) || (bl2<0) || (tl2<0)) SetShapeBit(kGeoRunTimeShape);
else TGeoArb8::ComputeBBox();
}
//_____________________________________________________________________________
TGeoGtra::TGeoGtra(const char *name, Double_t dz, Double_t theta, Double_t phi, Double_t twist, Double_t h1,
Double_t bl1, Double_t tl1, Double_t alpha1, Double_t h2, Double_t bl2,
Double_t tl2, Double_t alpha2)
:TGeoTrap(name, dz, theta, phi, h1, bl1, tl1, alpha1, h2, bl2, tl2, alpha2)
{
// constructor.
SetName(name);
fTheta = theta;
fPhi = phi;
fH1 = h1;
fH2 = h2;
fBl1 = bl1;
fBl2 = bl2;
fTl1 = tl1;
fTl2 = tl2;
fAlpha1 = alpha1;
fAlpha2 = alpha2;
Double_t x, y, dx, dy, dx1, dx2, th, ph, al1, al2;
al1 = alpha1*TMath::DegToRad();
al2 = alpha2*TMath::DegToRad();
th = theta*TMath::DegToRad();
ph = phi*TMath::DegToRad();
dx = 2*dz*TMath::Sin(th)*TMath::Cos(ph);
dy = 2*dz*TMath::Sin(th)*TMath::Sin(ph);
fDz = dz;
dx1 = 2*h1*TMath::Tan(al1);
dx2 = 2*h2*TMath::Tan(al2);
fTwistAngle = twist;
Int_t i;
for (i=0; i<8; i++) {
fXY[i][0] = 0.0;
fXY[i][1] = 0.0;
}
fXY[0][0] = -bl1; fXY[0][1] = -h1;
fXY[1][0] = -tl1+dx1; fXY[1][1] = h1;
fXY[2][0] = tl1+dx1; fXY[2][1] = h1;
fXY[3][0] = bl1; fXY[3][1] = -h1;
fXY[4][0] = -bl2+dx; fXY[4][1] = -h2+dy;
fXY[5][0] = -tl2+dx+dx2; fXY[5][1] = h2+dy;
fXY[6][0] = tl2+dx+dx2; fXY[6][1] = h2+dy;
fXY[7][0] = bl2+dx; fXY[7][1] = -h2+dy;
for (i=4; i<8; i++) {
x = fXY[i][0];
y = fXY[i][1];
fXY[i][0] = x*TMath::Cos(twist*TMath::DegToRad()) + y*TMath::Sin(twist*TMath::DegToRad());
fXY[i][1] = -x*TMath::Sin(twist*TMath::DegToRad()) + y*TMath::Cos(twist*TMath::DegToRad());
}
ComputeTwist();
if ((dz<0) || (h1<0) || (bl1<0) || (tl1<0) ||
(h2<0) || (bl2<0) || (tl2<0)) SetShapeBit(kGeoRunTimeShape);
else TGeoArb8::ComputeBBox();
}
//_____________________________________________________________________________
TGeoGtra::~TGeoGtra()
{
// destructor
}
//_____________________________________________________________________________
Double_t TGeoGtra::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 arb8
if (iact<3 && safe) {
// compute safe distance
*safe = Safety(point, kTRUE);
if (iact==0) return TGeoShape::Big();
if (iact==1 && step<*safe) return TGeoShape::Big();
}
// compute distance to get ouside this shape
return TGeoArb8::DistFromInside(point, dir, iact, step, safe);
}
//_____________________________________________________________________________
Double_t TGeoGtra::DistFromOutside(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 arb8
if (iact<3 && safe) {
// compute safe distance
*safe = Safety(point, kTRUE);
if (iact==0) return TGeoShape::Big();
if (iact==1 && step<*safe) return TGeoShape::Big();
}
// compute distance to get ouside this shape
return TGeoArb8::DistFromOutside(point, dir, iact, step, safe);
}
//_____________________________________________________________________________
TGeoShape *TGeoGtra::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;
if (mother->IsRunTimeShape()) {
Error("GetMakeRuntimeShape", "invalid mother");
return 0;
}
Double_t dz, h1, bl1, tl1, h2, bl2, tl2;
if (fDz<0) dz=((TGeoTrap*)mother)->GetDz();
else dz=fDz;
if (fH1<0)
h1 = ((TGeoTrap*)mother)->GetH1();
else
h1 = fH1;
if (fH2<0)
h2 = ((TGeoTrap*)mother)->GetH2();
else
h2 = fH2;
if (fBl1<0)
bl1 = ((TGeoTrap*)mother)->GetBl1();
else
bl1 = fBl1;
if (fBl2<0)
bl2 = ((TGeoTrap*)mother)->GetBl2();
else
bl2 = fBl2;
if (fTl1<0)
tl1 = ((TGeoTrap*)mother)->GetTl1();
else
tl1 = fTl1;
if (fTl2<0)
tl2 = ((TGeoTrap*)mother)->GetTl2();
else
tl2 = fTl2;
return (new TGeoGtra(dz, fTheta, fPhi, fTwistAngle ,h1, bl1, tl1, fAlpha1, h2, bl2, tl2, fAlpha2));
}
//_____________________________________________________________________________
void TGeoGtra::SavePrimitive(ofstream &out, Option_t * /*option*/)
{
// Save a primitive as a C++ statement(s) on output stream "out".
if (TObject::TestBit(kGeoSavePrimitive)) return;
out << " // Shape: " << GetName() << " type: " << ClassName() << endl;
out << " dz = " << fDz << ";" << endl;
out << " theta = " << fTheta << ";" << endl;
out << " phi = " << fPhi << ";" << endl;
out << " twist = " << fTwistAngle << ";" << endl;
out << " h1 = " << fH1<< ";" << endl;
out << " bl1 = " << fBl1<< ";" << endl;
out << " tl1 = " << fTl1<< ";" << endl;
out << " alpha1 = " << fAlpha1 << ";" << endl;
out << " h2 = " << fH2 << ";" << endl;
out << " bl2 = " << fBl2<< ";" << endl;
out << " tl2 = " << fTl2<< ";" << endl;
out << " alpha2 = " << fAlpha2 << ";" << endl;
out << " TGeoShape *" << GetPointerName() << " = new TGeoGtra(\"" << GetName() << "\", dz,theta,phi,twist,h1,bl1,tl1,alpha1,h2,bl2,tl2,alpha2);" << endl;
TObject::SetBit(TGeoShape::kGeoSavePrimitive);
}
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.