// @(#)root/geom:$Name: $:$Id: TGeoArb8.cxx,v 1.22 2003/03/14 11:49:02 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 "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;
}
}
//-----------------------------------------------------------------------------
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;
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;
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;
}
//-----------------------------------------------------------------------------
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 al1, al2, 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];
dx2 = fXY[4+(i+1)%4][0]-fXY[4+i][0];
dy2 = fXY[4+(i+1)%4][1]-fXY[4+i][1];
if (((dx1==0)&&(dy1==0)) || ((dx2==0)&&(dy2==0))) {
twist[i] = 0;
continue;
}
al1 = TMath::ATan2(dy1,dx1);
if (al1<0) al1+=2*TMath::Pi();
al2 = TMath::ATan2(dy2,dx2);
if (al2<0) al2+=2*TMath::Pi();
twist[i] = TMath::Tan(al2-al1);
if (twist[i]<1E-3) {
twist[i] = 0;
continue;
}
twisted = kTRUE;
}
if (!twisted) return;
if (fTwist) delete [] fTwist;
fTwist = new Double_t[4];
memcpy(fTwist, &twist[0], 4*sizeof(Double_t));
}
//-----------------------------------------------------------------------------
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]);
// find the intersections of Y=Ypoint with this poly.
Double_t segx[4];
Double_t x1, x2, y1, y2;
Int_t npts = 0;
Int_t i;
for (i=0; i<4; i++) {
y1 = poly[2*(i%4)+1];
y2 = poly[2*((i+1)%4)+1];
// printf("check Yp=%f against y1=%f y2=%fn", point[1], y1, y2);
if ((point[1]-y1)*(y2-point[1])<0) continue;
x1 = poly[2*(i%4)];
x2 = poly[2*((i+1)%4)];
// printf(" x1=%f x2=%fn", x1,x2);
// check if point is on the line connecting points 1-2
if (y1 == y2) {
if ((point[0]<x1) || (point[0]>x2)) return kFALSE;
return kTRUE;
}
Double_t cf = (point[1]-y1)/(y2-y1);
segx[npts++] = x1+cf*(x2-x1);
// printf(" x1+cf*(x2-x1) : %f+%f*(%f-%f)=%fn",x1, cf, x2, x1, x1+cf*(x2-x1));
// sort intersection points by X
if (npts>1) {
if (segx[npts-2] > segx[npts-1]) {
x1 = segx[npts-2];
segx[npts-2] = segx[npts-1];
segx[npts-1] = x1;
}
}
}
// printf("Intersections with Y = %f (Xpoint=%f):n",point[1], point[0]);
// for (i=0;i<npts;i++) printf("%i : X=%fn", i,segx[i]);
if (npts == 0) return kFALSE;
if (npts == 2) {
if ((point[0]<segx[npts-2]) || (point[0]>segx[npts-1])) return kFALSE;
return kTRUE;
}
if (npts != 4) return kFALSE;
// intersection poly is not convex (4 points)
if ((point[0]<segx[0]) || (point[0]>segx[3])) return kFALSE;
if ((point[0]>segx[1]) && (point[0]<segx[2])) return kFALSE;
return kTRUE;
}
//-----------------------------------------------------------------------------
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;
xa=fXY[ipl][0];
ya=fXY[ipl][1];
xb=fXY[ipl+4][0];
yb=fXY[ipl+4][1];
xc=fXY[(ipl+1)%4][0];
yc=fXY[(ipl+1)%4][1];
xd=fXY[4+(ipl+1)%4][0];
yd=fXY[4+(ipl+1)%4][1];
Double_t tx1=0.5*(xb-xa)/fDz;
Double_t ty1=0.5*(yb-ya)/fDz;
Double_t tx2=0.5*(xd-xc)/fDz;
Double_t ty2=0.5*(yd-yc)/fDz;
Double_t xs1=xa+tx1*(fDz+point[2]);
Double_t ys1=ya+ty1*(fDz+point[2]);
Double_t xs2=xc+tx2*(fDz+point[2]);
Double_t ys2=yc+ty2*(fDz+point[2]);
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=kBig;
Double_t x1,x2,y1,y2,xp,yp,zi;
if (TMath::Abs(a)<1E-10) {
if (b==0) return kBig;
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 kBig;
}
b=0.5*b/a;
c=c/a;
Double_t d=b*b-c;
if (d>=0) {
s=-b-TMath::Sqrt(d);
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+TMath::Sqrt(d);
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 kBig;
}
//-----------------------------------------------------------------------------
Double_t TGeoArb8::DistToIn(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=kBig;
if (!TGeoBBox::Contains(point)) {
snxt=TGeoBBox::DistToIn(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]=kBig;
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]=kBig;
} else {
for (Int_t j=0; j<2; j++) pt[j]=point[j]+dist[4]*dir[j];
if (!Contains(&pt[0])) dist[4]=kBig;
}
}
}
return dist[TMath::LocMin(5, &dist[0])];
}
//-----------------------------------------------------------------------------
Double_t TGeoArb8::DistToOut(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
Double_t dist[6];
dist[0]=dist[1]=kBig;
if (dir[2]<0) {
dist[0]=(-fDz-point[2])/dir[2];
} else {
if (dir[2]>0) dist[1]=(fDz-point[2])/dir[2];
}
for (Int_t i=0; i<4; i++) {
dist[i+2]=DistToPlane(point, dir, i, kTRUE);
}
return dist[TMath::LocMin(6, &dist[0])];
}
//-----------------------------------------------------------------------------
Double_t TGeoArb8::DistToSurf(Double_t * /*point*/, Double_t * /*dir*/) const
{
// computes the distance to next surface of the sphere along a ray
// starting from given point to the given direction.
return kBig;
}
//-----------------------------------------------------------------------------
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::kBig;
for (Int_t iaxis=0; iaxis<2; iaxis++) {
if (dd[iaxis]>=0) continue;
ddmin=TGeoShape::kBig;
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::InspectShape() const
{
// print shape parameters
printf("*** TGeoArb8 parameters ***n");
for (Int_t ip=0; ip<8; ip++) {
printf(" point #%i : x=%11.5f y=%11.5f z=%11.5fn",
ip, fXY[ip][0], fXY[ip][1], fDz*((ip<4)?-1:1));
}
TGeoBBox::InspectShape();
}
//-----------------------------------------------------------------------------
void TGeoArb8::NextCrossing(TGeoParamCurve * /*c*/, Double_t * /*point*/) const
{
// computes next intersection point of curve c with this shape
}
//-----------------------------------------------------------------------------
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.
return kBig;
}
//-----------------------------------------------------------------------------
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 *buff) const
{
// create arb8 mesh points
for (Int_t i=0; i<8; i++) {
buff[3*i] = fXY[i][0];
buff[3*i+1] = fXY[i][1];
buff[3*i+2] = (i<4)?-fDz:fDz;
// printf("%i X=%9.3fY=%9.3fZ=%9.3fn",i,buff[3*i],buff[3*i+1],buff[3*i+2] );
}
}
//-----------------------------------------------------------------------------
void TGeoArb8::SetPoints(Float_t *buff) const
{
// create arb8 mesh points
for (Int_t i=0; i<8; i++) {
buff[3*i] = fXY[i][0];
buff[3*i+1] = fXY[i][1];
buff[3*i+2] = (i<4)?-fDz:fDz;
// printf("%i X=%9.3fY=%9.3fZ=%9.3fn",i,buff[3*i],buff[3*i+1],buff[3*i+2] );
}
}
//-----------------------------------------------------------------------------
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*kDegRad)*TMath::Cos(phi*kDegRad);
Double_t ty = TMath::Tan(theta*kDegRad)*TMath::Sin(phi*kDegRad);
Double_t ta1 = TMath::Tan(alpha1*kDegRad);
Double_t ta2 = TMath::Tan(alpha2*kDegRad);
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)) {
SetBit(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*kDegRad)*TMath::Cos(phi*kDegRad);
Double_t ty = TMath::Tan(theta*kDegRad)*TMath::Sin(phi*kDegRad);
Double_t ta1 = TMath::Tan(alpha1*kDegRad);
Double_t ta2 = TMath::Tan(alpha2*kDegRad);
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)) {
SetBit(kGeoRunTimeShape);
}
else TGeoArb8::ComputeBBox();
}
//-----------------------------------------------------------------------------
TGeoTrap::~TGeoTrap()
{
// destructor
}
//-----------------------------------------------------------------------------
Double_t TGeoTrap::DistToOut(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 kBig;
if (iact==1 && step<*safe) return kBig;
}
// compute distance to get ouside this shape
return TGeoArb8::DistToOut(point, dir, iact, step, safe);
}
//-----------------------------------------------------------------------------
Double_t TGeoTrap::DistToIn(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 kBig;
if (iact==1 && step<*safe) return kBig;
}
// compute distance to get ouside this shape
return TGeoArb8::DistToIn(point, dir, iact, step, safe);
}
//-----------------------------------------------------------------------------
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 (!TestBit(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 = kBig;
Double_t saf[5];
Double_t norm[3]; // normal to current facette
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;
//---> compute safety for lateral planes
for (i=0; i<4; i++) {
x0 = fXY[i][0];
y0 = fXY[i][1];
z0 = -fDz;
x1 = fXY[i+4][0];
y1 = fXY[i+4][1];
z1 = fDz;
ax = x1-x0;
ay = y1-y0;
az = z1-z0;
x2 = fXY[(i+1)%4][0];
y2 = fXY[(i+1)%4][1];
bx = x2-x0;
by = y2-y0;
if (bx==0 && by==0) {
x2 = fXY[4+((i+1)%4)][0];
y2 = fXY[4+((i+1)%4)][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;
norm[0] /= fn;
norm[1] /= fn;
norm[2] /= fn;
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]); // they should be all positive anyway
} else {
saf[i] = -saf[i]; // only negative values are interesting
}
}
saf[4] = fDz-TMath::Abs(point[2]);
if (in) {
safe = saf[TMath::LocMin(5, saf)];
} else {
saf[4]=-saf[4];
safe = saf[TMath::LocMax(5, saf)];
}
return safe;
}
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*kDegRad;
al2 = alpha2*kDegRad;
th = theta*kDegRad;
ph = phi*kDegRad;
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*kDegRad) + y*TMath::Sin(twist*kDegRad);
fXY[i][1] = -x*TMath::Sin(twist*kDegRad) + y*TMath::Cos(twist*kDegRad);
}
ComputeTwist();
if ((dz<0) || (h1<0) || (bl1<0) || (tl1<0) ||
(h2<0) || (bl2<0) || (tl2<0)) SetBit(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*kDegRad;
al2 = alpha2*kDegRad;
th = theta*kDegRad;
ph = phi*kDegRad;
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*kDegRad) + y*TMath::Sin(twist*kDegRad);
fXY[i][1] = -x*TMath::Sin(twist*kDegRad) + y*TMath::Cos(twist*kDegRad);
}
ComputeTwist();
if ((dz<0) || (h1<0) || (bl1<0) || (tl1<0) ||
(h2<0) || (bl2<0) || (tl2<0)) SetBit(kGeoRunTimeShape);
else TGeoArb8::ComputeBBox();
}
//-----------------------------------------------------------------------------
TGeoGtra::~TGeoGtra()
{
// destructor
}
//-----------------------------------------------------------------------------
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 (!TestBit(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));
}
ROOT page - Class index - 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.