Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GenericTrap Class Reference

#include <G4GenericTrap.hh>

Inheritance diagram for G4GenericTrap:
Collaboration diagram for G4GenericTrap:

Public Member Functions

 G4GenericTrap (const G4String &name, G4double halfZ, const std::vector< G4TwoVector > &vertices)
 
 ~G4GenericTrap ()
 
G4double GetZHalfLength () const
 
G4int GetNofVertices () const
 
G4TwoVector GetVertex (G4int index) const
 
const std::vector< G4TwoVector > & GetVertices () const
 
G4double GetTwistAngle (G4int index) const
 
G4bool IsTwisted () const
 
G4int GetVisSubdivisions () const
 
void SetVisSubdivisions (G4int subdiv)
 
EInside Inside (const G4ThreeVector &p) const
 
G4ThreeVector SurfaceNormal (const G4ThreeVector &p) const
 
G4double DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const
 
G4double DistanceToIn (const G4ThreeVector &p) const
 
G4double DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
 
G4double DistanceToOut (const G4ThreeVector &p) const
 
void Extent (G4ThreeVector &pMin, G4ThreeVector &pMax) const
 
G4bool CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
 
G4GeometryType GetEntityType () const
 
G4VSolidClone () const
 
std::ostream & StreamInfo (std::ostream &os) const
 
G4ThreeVector GetPointOnSurface () const
 
G4double GetCubicVolume ()
 
G4double GetSurfaceArea ()
 
G4PolyhedronGetPolyhedron () const
 
void DescribeYourselfTo (G4VGraphicsScene &scene) const
 
G4VisExtent GetExtent () const
 
G4PolyhedronCreatePolyhedron () const
 
 G4GenericTrap (__void__ &)
 
 G4GenericTrap (const G4GenericTrap &rhs)
 
G4GenericTrapoperator= (const G4GenericTrap &rhs)
 
- Public Member Functions inherited from G4VSolid
 G4VSolid (const G4String &name)
 
virtual ~G4VSolid ()
 
G4bool operator== (const G4VSolid &s) const
 
G4String GetName () const
 
void SetName (const G4String &name)
 
G4double GetTolerance () const
 
virtual void ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
 
void DumpInfo () const
 
virtual const G4VSolidGetConstituentSolid (G4int no) const
 
virtual G4VSolidGetConstituentSolid (G4int no)
 
virtual const G4DisplacedSolidGetDisplacedSolidPtr () const
 
virtual G4DisplacedSolidGetDisplacedSolidPtr ()
 
 G4VSolid (__void__ &)
 
 G4VSolid (const G4VSolid &rhs)
 
G4VSolidoperator= (const G4VSolid &rhs)
 
G4double EstimateCubicVolume (G4int nStat, G4double epsilon) const
 
G4double EstimateSurfaceArea (G4int nStat, G4double ell) const
 

Protected Attributes

G4bool fRebuildPolyhedron
 
G4PolyhedronfpPolyhedron
 
- Protected Attributes inherited from G4VSolid
G4double kCarTolerance
 

Additional Inherited Members

- Protected Member Functions inherited from G4VSolid
void CalculateClippedPolygonExtent (G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
 
void ClipCrossSection (G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
 
void ClipBetweenSections (G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const
 
void ClipPolygon (G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis) const
 

Detailed Description

Definition at line 88 of file G4GenericTrap.hh.

Constructor & Destructor Documentation

G4GenericTrap::G4GenericTrap ( const G4String name,
G4double  halfZ,
const std::vector< G4TwoVector > &  vertices 
)

Definition at line 81 of file G4GenericTrap.cc.

83  : G4VSolid(name),
84  fRebuildPolyhedron(false),
85  fpPolyhedron(0),
86  fDz(halfZ),
87  fVertices(),
88  fIsTwisted(false),
89  fTessellatedSolid(0),
90  fMinBBoxVector(G4ThreeVector(0,0,0)),
91  fMaxBBoxVector(G4ThreeVector(0,0,0)),
92  fVisSubdivisions(0),
93  fSurfaceArea(0.),
94  fCubicVolume(0.)
95 
96 {
97  // General constructor
98  const G4double min_length=5*1.e-6;
99  G4double length = 0.;
100  G4int k=0;
101  G4String errorDescription = "InvalidSetup in \" ";
102  errorDescription += name;
103  errorDescription += "\"";
104 
105  halfCarTolerance = kCarTolerance*0.5;
106 
107  // Check vertices size
108 
109  if ( G4int(vertices.size()) != fgkNofVertices )
110  {
111  G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids0002",
112  FatalErrorInArgument, "Number of vertices != 8");
113  }
114 
115  // Check dZ
116  //
117  if (halfZ < kCarTolerance)
118  {
119  G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids0002",
120  FatalErrorInArgument, "dZ is too small or negative");
121  }
122 
123  // Check Ordering and Copy vertices
124  //
125  if(CheckOrder(vertices))
126  {
127  for (G4int i=0; i<fgkNofVertices; ++i) {fVertices.push_back(vertices[i]);}
128  }
129  else
130  {
131  for (G4int i=0; i <4; ++i) {fVertices.push_back(vertices[3-i]);}
132  for (G4int i=0; i <4; ++i) {fVertices.push_back(vertices[7-i]);}
133  }
134 
135  // Check length of segments and Adjust
136  //
137  for (G4int j=0; j < 2; j++)
138  {
139  for (G4int i=1; i<4; ++i)
140  {
141  k = j*4+i;
142  length = (fVertices[k]-fVertices[k-1]).mag();
143  if ( ( length < min_length) && ( length > kCarTolerance ) )
144  {
145  std::ostringstream message;
146  message << "Length segment is too small." << G4endl
147  << "Distance between " << fVertices[k-1] << " and "
148  << fVertices[k] << " is only " << length << " mm !";
149  G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids1001",
150  JustWarning, message, "Vertices will be collapsed.");
151  fVertices[k]=fVertices[k-1];
152  }
153  }
154  }
155 
156  // Compute Twist
157  //
158  for( G4int i=0; i<4; i++) { fTwist[i]=0.; }
159  fIsTwisted = ComputeIsTwisted();
160 
161  // Compute Bounding Box
162  //
163  ComputeBBox();
164 
165  // If not twisted - create tessellated solid
166  // (an alternative implementation for testing)
167  //
168 #ifdef G4TESS_TEST
169  if ( !fIsTwisted ) { fTessellatedSolid = CreateTessellatedSolid(); }
170 #endif
171 }
const XML_Char * name
Definition: expat.h:151
CLHEP::Hep3Vector G4ThreeVector
G4Polyhedron * fpPolyhedron
int G4int
Definition: G4Types.hh:78
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
#define G4endl
Definition: G4ios.hh:61
G4VSolid(const G4String &name)
Definition: G4VSolid.cc:61
G4double kCarTolerance
Definition: G4VSolid.hh:307
double G4double
Definition: G4Types.hh:76
G4bool fRebuildPolyhedron

Here is the call graph for this function:

Here is the caller graph for this function:

G4GenericTrap::~G4GenericTrap ( )

Definition at line 196 of file G4GenericTrap.cc.

197 {
198  // Destructor
199  delete fTessellatedSolid;
200 }
G4GenericTrap::G4GenericTrap ( __void__ &  a)

Definition at line 175 of file G4GenericTrap.cc.

176  : G4VSolid(a),
177  fRebuildPolyhedron(false),
178  fpPolyhedron(0),
179  halfCarTolerance(0.),
180  fDz(0.),
181  fVertices(),
182  fIsTwisted(false),
183  fTessellatedSolid(0),
184  fMinBBoxVector(G4ThreeVector(0,0,0)),
185  fMaxBBoxVector(G4ThreeVector(0,0,0)),
186  fVisSubdivisions(0),
187  fSurfaceArea(0.),
188  fCubicVolume(0.)
189 {
190  // Fake default constructor - sets only member data and allocates memory
191  // for usage restricted to object persistency.
192 }
CLHEP::Hep3Vector G4ThreeVector
G4Polyhedron * fpPolyhedron
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4VSolid(const G4String &name)
Definition: G4VSolid.cc:61
G4bool fRebuildPolyhedron
G4GenericTrap::G4GenericTrap ( const G4GenericTrap rhs)

Definition at line 204 of file G4GenericTrap.cc.

205  : G4VSolid(rhs),
206  fRebuildPolyhedron(false), fpPolyhedron(0),
207  halfCarTolerance(rhs.halfCarTolerance),
208  fDz(rhs.fDz), fVertices(rhs.fVertices),
209  fIsTwisted(rhs.fIsTwisted), fTessellatedSolid(0),
210  fMinBBoxVector(rhs.fMinBBoxVector), fMaxBBoxVector(rhs.fMaxBBoxVector),
211  fVisSubdivisions(rhs.fVisSubdivisions),
212  fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume)
213 {
214  for (size_t i=0; i<4; ++i) { fTwist[i] = rhs.fTwist[i]; }
215 #ifdef G4TESS_TEST
216  if (rhs.fTessellatedSolid && !fIsTwisted )
217  { fTessellatedSolid = CreateTessellatedSolid(); }
218 #endif
219 }
G4Polyhedron * fpPolyhedron
G4VSolid(const G4String &name)
Definition: G4VSolid.cc:61
G4bool fRebuildPolyhedron

Member Function Documentation

G4bool G4GenericTrap::CalculateExtent ( const EAxis  pAxis,
const G4VoxelLimits pVoxelLimit,
const G4AffineTransform pTransform,
G4double pmin,
G4double pmax 
) const
virtual

Implements G4VSolid.

Definition at line 1237 of file G4GenericTrap.cc.

1241 {
1242  G4ThreeVector bmin, bmax;
1243  G4bool exist;
1244 
1245  // Check bounding box (bbox)
1246  //
1247  Extent(bmin,bmax);
1248  G4BoundingEnvelope bbox(bmin,bmax);
1249 #ifdef G4BBOX_EXTENT
1250  if (true) return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
1251 #endif
1252  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
1253  {
1254  return exist = (pMin < pMax) ? true : false;
1255  }
1256 
1257  // Set bounding envelope (benv) and calculate extent
1258  //
1259  // To build the bounding envelope with plane faces each side face of
1260  // the trapezoid is subdivided in triangles. Subdivision is done by
1261  // duplication of vertices in the bases in a way that the envelope be
1262  // a convex polyhedron (some faces of the envelope can be degenerate)
1263  //
1264  G4double dz = GetZHalfLength();
1265  G4ThreeVectorList baseA(8), baseB(8);
1266  for (G4int i=0; i<4; ++i)
1267  {
1268  G4TwoVector va = GetVertex(i);
1269  G4TwoVector vb = GetVertex(i+4);
1270  baseA[2*i].set(va.x(),va.y(),-dz);
1271  baseB[2*i].set(vb.x(),vb.y(), dz);
1272  }
1273  for (G4int i=0; i<4; ++i)
1274  {
1275  G4int k1=2*i, k2=(2*i+2)%8;
1276  G4double ax = (baseA[k2].x()-baseA[k1].x());
1277  G4double ay = (baseA[k2].y()-baseA[k1].y());
1278  G4double bx = (baseB[k2].x()-baseB[k1].x());
1279  G4double by = (baseB[k2].y()-baseB[k1].y());
1280  G4double znorm = ax*by - ay*bx;
1281  baseA[k1+1] = (znorm < 0.0) ? baseA[k2] : baseA[k1];
1282  baseB[k1+1] = (znorm < 0.0) ? baseB[k1] : baseB[k2];
1283  }
1284 
1285  std::vector<const G4ThreeVectorList *> polygons(2);
1286  polygons[0] = &baseA;
1287  polygons[1] = &baseB;
1288 
1289  G4BoundingEnvelope benv(bmin,bmax,polygons);
1290  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
1291  return exist;
1292 }
double y() const
double x() const
G4TwoVector GetVertex(G4int index) const
G4double GetZHalfLength() const
int G4int
Definition: G4Types.hh:78
bool G4bool
Definition: G4Types.hh:79
std::vector< G4ThreeVector > G4ThreeVectorList
void Extent(G4ThreeVector &pMin, G4ThreeVector &pMax) const
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4VSolid * G4GenericTrap::Clone ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 1303 of file G4GenericTrap.cc.

1304 {
1305  return new G4GenericTrap(*this);
1306 }
G4GenericTrap(const G4String &name, G4double halfZ, const std::vector< G4TwoVector > &vertices)

Here is the call graph for this function:

G4Polyhedron * G4GenericTrap::CreatePolyhedron ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 2033 of file G4GenericTrap.cc.

2034 {
2035 
2036 #ifdef G4TESS_TEST
2037  if ( fTessellatedSolid )
2038  {
2039  return fTessellatedSolid->CreatePolyhedron();
2040  }
2041 #endif
2042 
2043  // Approximation of Twisted Side
2044  // Construct extra Points, if Twisted Side
2045  //
2046  G4PolyhedronArbitrary* polyhedron;
2047  size_t nVertices, nFacets;
2048 
2049  G4int subdivisions=0;
2050  G4int i;
2051  if(fIsTwisted)
2052  {
2053  if ( GetVisSubdivisions()!= 0 )
2054  {
2055  subdivisions=GetVisSubdivisions();
2056  }
2057  else
2058  {
2059  // Estimation of Number of Subdivisions for smooth visualisation
2060  //
2061  G4double maxTwist=0.;
2062  for(i=0; i<4; i++)
2063  {
2064  if(GetTwistAngle(i)>maxTwist) { maxTwist=GetTwistAngle(i); }
2065  }
2066 
2067  // Computes bounding vectors for the shape
2068  //
2069  G4double Dx,Dy;
2070  G4ThreeVector minVec = GetMinimumBBox();
2071  G4ThreeVector maxVec = GetMaximumBBox();
2072  Dx = 0.5*(maxVec.x()- minVec.y());
2073  Dy = 0.5*(maxVec.y()- minVec.y());
2074  if (Dy > Dx) { Dx=Dy; }
2075 
2076  subdivisions=8*G4int(maxTwist/(Dx*Dx*Dx)*fDz);
2077  if (subdivisions<4) { subdivisions=4; }
2078  if (subdivisions>30) { subdivisions=30; }
2079  }
2080  }
2081  G4int sub4=4*subdivisions;
2082  nVertices = 8+subdivisions*4;
2083  nFacets = 6+subdivisions*4;
2084  G4double cf=1./(subdivisions+1);
2085  polyhedron = new G4PolyhedronArbitrary (nVertices, nFacets);
2086 
2087  // Add Vertex
2088  //
2089  for (i=0;i<4;i++)
2090  {
2091  polyhedron->AddVertex(G4ThreeVector(fVertices[i].x(),
2092  fVertices[i].y(),-fDz));
2093  }
2094  for( i=0;i<subdivisions;i++)
2095  {
2096  for(G4int j=0;j<4;j++)
2097  {
2098  G4TwoVector u=fVertices[j]+cf*(i+1)*( fVertices[j+4]-fVertices[j]);
2099  polyhedron->AddVertex(G4ThreeVector(u.x(),u.y(),-fDz+cf*2*fDz*(i+1)));
2100  }
2101  }
2102  for (i=4;i<8;i++)
2103  {
2104  polyhedron->AddVertex(G4ThreeVector(fVertices[i].x(),
2105  fVertices[i].y(),fDz));
2106  }
2107 
2108  // Add Facets
2109  //
2110  polyhedron->AddFacet(1,4,3,2); //Z-plane
2111  for (i=0;i<subdivisions+1;i++)
2112  {
2113  G4int is=i*4;
2114  polyhedron->AddFacet(5+is,8+is,4+is,1+is);
2115  polyhedron->AddFacet(8+is,7+is,3+is,4+is);
2116  polyhedron->AddFacet(7+is,6+is,2+is,3+is);
2117  polyhedron->AddFacet(6+is,5+is,1+is,2+is);
2118  }
2119  polyhedron->AddFacet(5+sub4,6+sub4,7+sub4,8+sub4); //Z-plane
2120 
2121  polyhedron->SetReferences();
2122  polyhedron->InvertFacets();
2123 
2124  return (G4Polyhedron*) polyhedron;
2125 }
double y() const
G4int GetVisSubdivisions() const
double x() const
CLHEP::Hep3Vector G4ThreeVector
double x() const
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
virtual G4Polyhedron * CreatePolyhedron() const
G4double GetTwistAngle(G4int index) const
void AddVertex(const G4ThreeVector &v)
double y() const
void AddFacet(const G4int iv1, const G4int iv2, const G4int iv3, const G4int iv4=0)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4GenericTrap::DescribeYourselfTo ( G4VGraphicsScene scene) const
virtual

Implements G4VSolid.

Definition at line 1998 of file G4GenericTrap.cc.

1999 {
2000 
2001 #ifdef G4TESS_TEST
2002  if ( fTessellatedSolid )
2003  {
2004  return fTessellatedSolid->DescribeYourselfTo(scene);
2005  }
2006 #endif
2007 
2008  scene.AddSolid(*this);
2009 }
virtual void DescribeYourselfTo(G4VGraphicsScene &scene) const
virtual void AddSolid(const G4Box &)=0

Here is the call graph for this function:

G4double G4GenericTrap::DistanceToIn ( const G4ThreeVector p,
const G4ThreeVector v 
) const
virtual

Implements G4VSolid.

Definition at line 771 of file G4GenericTrap.cc.

773 {
774 #ifdef G4TESS_TEST
775  if ( fTessellatedSolid )
776  {
777  return fTessellatedSolid->DistanceToIn(p, v);
778  }
779 #endif
780 
781  G4double dist[5];
783 
784  // Check lateral faces
785  //
786  G4int i;
787  for (i=0; i<4; i++)
788  {
789  dist[i]=DistToPlane(p, v, i);
790  }
791 
792  // Check Z planes
793  //
794  dist[4]=kInfinity;
795  if (std::fabs(p.z())>fDz-halfCarTolerance)
796  {
797  if (v.z())
798  {
799  G4ThreeVector pt;
800  if (p.z()>0)
801  {
802  dist[4] = (fDz-p.z())/v.z();
803  }
804  else
805  {
806  dist[4] = (-fDz-p.z())/v.z();
807  }
808  if (dist[4]<-halfCarTolerance)
809  {
810  dist[4]=kInfinity;
811  }
812  else
813  {
814  if(dist[4]<halfCarTolerance)
815  {
816  if(p.z()>0) { n=G4ThreeVector(0,0,1); }
817  else { n=G4ThreeVector(0,0,-1); }
818  if (n.dot(v)<0) { dist[4]=0.; }
819  else { dist[4]=kInfinity; }
820  }
821  pt=p+dist[4]*v;
822  if (Inside(pt)==kOutside) { dist[4]=kInfinity; }
823  }
824  }
825  }
826  G4double distmin = dist[0];
827  for (i=1;i<5;i++)
828  {
829  if (dist[i] < distmin) { distmin = dist[i]; }
830  }
831 
832  if (distmin<halfCarTolerance) { distmin=0.; }
833 
834  return distmin;
835 }
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
double dot(const Hep3Vector &) const
int G4int
Definition: G4Types.hh:78
double z() const
const G4int n
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
EInside Inside(const G4ThreeVector &p) const
tuple v
Definition: test.py:18
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4GenericTrap::DistanceToIn ( const G4ThreeVector p) const
virtual

Implements G4VSolid.

Definition at line 839 of file G4GenericTrap.cc.

840 {
841  // Computes the closest distance from given point to this shape
842 
843 #ifdef G4TESS_TEST
844  if ( fTessellatedSolid )
845  {
846  return fTessellatedSolid->DistanceToIn(p);
847  }
848 #endif
849 
850  G4double safz = std::fabs(p.z())-fDz;
851  if(safz<0) { safz=0; }
852 
853  G4int iseg;
854  G4double safe = safz;
855  G4double safxy = safz;
856 
857  for (iseg=0; iseg<4; iseg++)
858  {
859  safxy = SafetyToFace(p,iseg);
860  if (safxy>safe) { safe=safxy; }
861  }
862 
863  return safe;
864 }
int G4int
Definition: G4Types.hh:78
double z() const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4GenericTrap::DistanceToOut ( const G4ThreeVector p,
const G4ThreeVector v,
const G4bool  calcNorm = false,
G4bool validNorm = 0,
G4ThreeVector n = 0 
) const
virtual

Implements G4VSolid.

Definition at line 943 of file G4GenericTrap.cc.

948 {
949 #ifdef G4TESS_TEST
950  if ( fTessellatedSolid )
951  {
952  return fTessellatedSolid->DistanceToOut(p, v, calcNorm, validNorm, n);
953  }
954 #endif
955 
956  G4double distmin;
957  G4bool lateral_cross = false;
958  ESide side = kUndefined;
959 
960  if (calcNorm) { *validNorm=true; } // All normals are valid
961 
962  if (v.z() < 0)
963  {
964  distmin=(-fDz-p.z())/v.z();
965  if (calcNorm) { side=kMZ; *n=G4ThreeVector(0,0,-1); }
966  }
967  else
968  {
969  if (v.z() > 0)
970  {
971  distmin = (fDz-p.z())/v.z();
972  if (calcNorm) { side=kPZ; *n=G4ThreeVector(0,0,1); }
973  }
974  else { distmin = kInfinity; }
975  }
976 
977  G4double dz2 =0.5/fDz;
978  G4double xa,xb,xc,xd;
979  G4double ya,yb,yc,yd;
980 
981  for (G4int ipl=0; ipl<4; ipl++)
982  {
983  G4int j = (ipl+1)%4;
984  xa=fVertices[ipl].x();
985  ya=fVertices[ipl].y();
986  xb=fVertices[ipl+4].x();
987  yb=fVertices[ipl+4].y();
988  xc=fVertices[j].x();
989  yc=fVertices[j].y();
990  xd=fVertices[4+j].x();
991  yd=fVertices[4+j].y();
992 
993  if ( ((std::fabs(xb-xd)+std::fabs(yb-yd))<halfCarTolerance)
994  || ((std::fabs(xa-xc)+std::fabs(ya-yc))<halfCarTolerance) )
995  {
996  G4double q=DistToTriangle(p,v,ipl) ;
997  if ( (q>=0) && (q<distmin) )
998  {
999  distmin=q;
1000  lateral_cross=true;
1001  side=ESide(ipl+1);
1002  }
1003  continue;
1004  }
1005  G4double tx1 =dz2*(xb-xa);
1006  G4double ty1 =dz2*(yb-ya);
1007  G4double tx2 =dz2*(xd-xc);
1008  G4double ty2 =dz2*(yd-yc);
1009  G4double dzp =fDz+p.z();
1010  G4double xs1 =xa+tx1*dzp;
1011  G4double ys1 =ya+ty1*dzp;
1012  G4double xs2 =xc+tx2*dzp;
1013  G4double ys2 =yc+ty2*dzp;
1014  G4double dxs =xs2-xs1;
1015  G4double dys =ys2-ys1;
1016  G4double dtx =tx2-tx1;
1017  G4double dty =ty2-ty1;
1018  G4double a = (dtx*v.y()-dty*v.x()+(tx1*ty2-tx2*ty1)*v.z())*v.z();
1019  G4double b = dxs*v.y()-dys*v.x()+(dtx*p.y()-dty*p.x()+ty2*xs1-ty1*xs2
1020  + tx1*ys2-tx2*ys1)*v.z();
1021  G4double c=dxs*p.y()-dys*p.x()+xs1*ys2-xs2*ys1;
1022  G4double q=kInfinity;
1023 
1024  if (std::fabs(a) < kCarTolerance)
1025  {
1026  if (std::fabs(b) < kCarTolerance) { continue; }
1027  q=-c/b;
1028 
1029  // Check for Point on the Surface
1030  //
1031  if ((q > -halfCarTolerance) && (q < distmin))
1032  {
1033  if (q < halfCarTolerance)
1034  {
1035  if (NormalToPlane(p,ipl).dot(v)<0.) { continue; }
1036  }
1037  distmin =q;
1038  lateral_cross=true;
1039  side=ESide(ipl+1);
1040  }
1041  continue;
1042  }
1043  G4double d=b*b-4*a*c;
1044  if (d >= 0.)
1045  {
1046  if (a > 0) { q=0.5*(-b-std::sqrt(d))/a; }
1047  else { q=0.5*(-b+std::sqrt(d))/a; }
1048 
1049  // Check for Point on the Surface
1050  //
1051  if (q > -halfCarTolerance )
1052  {
1053  if (q < distmin)
1054  {
1055  if(q < halfCarTolerance)
1056  {
1057  if (NormalToPlane(p,ipl).dot(v)<0.) // Check second root
1058  {
1059  if (a > 0) { q=0.5*(-b+std::sqrt(d))/a; }
1060  else { q=0.5*(-b-std::sqrt(d))/a; }
1061  if (( q > halfCarTolerance) && (q < distmin))
1062  {
1063  distmin=q;
1064  lateral_cross = true;
1065  side=ESide(ipl+1);
1066  }
1067  continue;
1068  }
1069  }
1070  distmin = q;
1071  lateral_cross = true;
1072  side=ESide(ipl+1);
1073  }
1074  }
1075  else
1076  {
1077  if (a > 0) { q=0.5*(-b+std::sqrt(d))/a; }
1078  else { q=0.5*(-b-std::sqrt(d))/a; }
1079 
1080  // Check for Point on the Surface
1081  //
1082  if ((q > -halfCarTolerance) && (q < distmin))
1083  {
1084  if (q < halfCarTolerance)
1085  {
1086  if (NormalToPlane(p,ipl).dot(v)<0.) // Check second root
1087  {
1088  if (a > 0) { q=0.5*(-b-std::sqrt(d))/a; }
1089  else { q=0.5*(-b+std::sqrt(d))/a; }
1090  if ( ( q > halfCarTolerance) && (q < distmin) )
1091  {
1092  distmin=q;
1093  lateral_cross = true;
1094  side=ESide(ipl+1);
1095  }
1096  continue;
1097  }
1098  }
1099  distmin =q;
1100  lateral_cross = true;
1101  side=ESide(ipl+1);
1102  }
1103  }
1104  }
1105  }
1106  if (!lateral_cross) // Make sure that track crosses the top or bottom
1107  {
1108  if (distmin >= kInfinity) { distmin=kCarTolerance; }
1109  G4ThreeVector pt=p+distmin*v;
1110 
1111  // Check if propagated point is in the polygon
1112  //
1113  G4int i=0;
1114  if (v.z()>0.) { i=4; }
1115  std::vector<G4TwoVector> xy;
1116  for ( G4int j=0; j<4; j++) { xy.push_back(fVertices[i+j]); }
1117 
1118  // Check Inside
1119  //
1120  if (InsidePolygone(pt,xy)==kOutside)
1121  {
1122  if(calcNorm)
1123  {
1124  if (v.z()>0) {side= kPZ; *n = G4ThreeVector(0,0,1);}
1125  else { side=kMZ; *n = G4ThreeVector(0,0,-1);}
1126  }
1127  return 0.;
1128  }
1129  else
1130  {
1131  if(v.z()>0) {side=kPZ;}
1132  else {side=kMZ;}
1133  }
1134  }
1135 
1136  if (calcNorm)
1137  {
1138  G4ThreeVector pt=p+v*distmin;
1139  switch (side)
1140  {
1141  case kXY0:
1142  *n=NormalToPlane(pt,0);
1143  break;
1144  case kXY1:
1145  *n=NormalToPlane(pt,1);
1146  break;
1147  case kXY2:
1148  *n=NormalToPlane(pt,2);
1149  break;
1150  case kXY3:
1151  *n=NormalToPlane(pt,3);
1152  break;
1153  case kPZ:
1154  *n=G4ThreeVector(0,0,1);
1155  break;
1156  case kMZ:
1157  *n=G4ThreeVector(0,0,-1);
1158  break;
1159  default:
1160  DumpInfo();
1161  std::ostringstream message;
1162  G4int oldprc = message.precision(16);
1163  message << "Undefined side for valid surface normal to solid." << G4endl
1164  << "Position:" << G4endl
1165  << " p.x() = " << p.x()/mm << " mm" << G4endl
1166  << " p.y() = " << p.y()/mm << " mm" << G4endl
1167  << " p.z() = " << p.z()/mm << " mm" << G4endl
1168  << "Direction:" << G4endl
1169  << " v.x() = " << v.x() << G4endl
1170  << " v.y() = " << v.y() << G4endl
1171  << " v.z() = " << v.z() << G4endl
1172  << "Proposed distance :" << G4endl
1173  << " distmin = " << distmin/mm << " mm";
1174  message.precision(oldprc);
1175  G4Exception("G4GenericTrap::DistanceToOut(p,v,..)",
1176  "GeomSolids1002", JustWarning, message);
1177  break;
1178  }
1179  }
1180 
1181  if (distmin<halfCarTolerance) { distmin=0.; }
1182 
1183  return distmin;
1184 }
Definition: G4Cons.cc:76
static constexpr double mm
Definition: G4SIunits.hh:115
static const G4double kInfinity
Definition: geomdefs.hh:42
CLHEP::Hep3Vector G4ThreeVector
double x() const
double dot(const Hep3Vector &) const
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
int G4int
Definition: G4Types.hh:78
double z() const
void DumpInfo() const
tuple b
Definition: test.py:12
virtual G4double DistanceToOut(const G4ThreeVector &p) const
bool G4bool
Definition: G4Types.hh:79
tuple v
Definition: test.py:18
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double y() const
Definition: G4Cons.cc:76
#define G4endl
Definition: G4ios.hh:61
G4double kCarTolerance
Definition: G4VSolid.hh:307
ESide
Definition: G4Cons.cc:76
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13

Here is the call graph for this function:

G4double G4GenericTrap::DistanceToOut ( const G4ThreeVector p) const
virtual

Implements G4VSolid.

Definition at line 1188 of file G4GenericTrap.cc.

1189 {
1190 
1191 #ifdef G4TESS_TEST
1192  if ( fTessellatedSolid )
1193  {
1194  return fTessellatedSolid->DistanceToOut(p);
1195  }
1196 #endif
1197 
1198  G4double safz = fDz-std::fabs(p.z());
1199  if (safz<0) { safz = 0; }
1200 
1201  G4double safe = safz;
1202  G4double safxy = safz;
1203 
1204  for (G4int iseg=0; iseg<4; iseg++)
1205  {
1206  safxy = std::fabs(SafetyToFace(p,iseg));
1207  if (safxy < safe) { safe = safxy; }
1208  }
1209 
1210  return safe;
1211 }
int G4int
Definition: G4Types.hh:78
double z() const
virtual G4double DistanceToOut(const G4ThreeVector &p) const
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4GenericTrap::Extent ( G4ThreeVector pMin,
G4ThreeVector pMax 
) const
virtual

Reimplemented from G4VSolid.

Definition at line 1215 of file G4GenericTrap.cc.

1216 {
1217  pMin = GetMinimumBBox();
1218  pMax = GetMaximumBBox();
1219 
1220  // Check correctness of the bounding box
1221  //
1222  if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
1223  {
1224  std::ostringstream message;
1225  message << "Bad bounding box (min >= max) for solid: "
1226  << GetName() << " !"
1227  << "\npMin = " << pMin
1228  << "\npMax = " << pMax;
1229  G4Exception("G4GenericTrap::Extent()", "GeomMgt0001", JustWarning, message);
1230  DumpInfo();
1231  }
1232 }
G4String GetName() const
double x() const
double z() const
void DumpInfo() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
double y() const

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4GenericTrap::GetCubicVolume ( )
virtual

Reimplemented from G4VSolid.

Definition at line 1451 of file G4GenericTrap.cc.

1452 {
1453  if (fCubicVolume == 0.0) {
1454  if(fIsTwisted) {
1455  fCubicVolume = G4VSolid::GetCubicVolume();
1456  } else {
1457  // Set vertices
1458  G4ThreeVector vertix0(fVertices[0].x(),fVertices[0].y(),-fDz);
1459  G4ThreeVector vertix1(fVertices[1].x(),fVertices[1].y(),-fDz);
1460  G4ThreeVector vertix2(fVertices[2].x(),fVertices[2].y(),-fDz);
1461  G4ThreeVector vertix3(fVertices[3].x(),fVertices[3].y(),-fDz);
1462  G4ThreeVector vertix4(fVertices[4].x(),fVertices[4].y(), fDz);
1463  G4ThreeVector vertix5(fVertices[5].x(),fVertices[5].y(), fDz);
1464  G4ThreeVector vertix6(fVertices[6].x(),fVertices[6].y(), fDz);
1465  G4ThreeVector vertix7(fVertices[7].x(),fVertices[7].y(), fDz);
1466 
1467  // Find Cubic Volume
1468  fCubicVolume = GetFaceCubicVolume(vertix0,vertix1,vertix2,vertix3) // -fDz plane
1469  + GetFaceCubicVolume(vertix1,vertix0,vertix4,vertix5) // Lat plane
1470  + GetFaceCubicVolume(vertix2,vertix1,vertix5,vertix6) // Lat plane
1471  + GetFaceCubicVolume(vertix3,vertix2,vertix6,vertix7) // Lat plane
1472  + GetFaceCubicVolume(vertix0,vertix3,vertix7,vertix4) // Lat plane
1473  + GetFaceCubicVolume(vertix7,vertix6,vertix5,vertix4); // +fDz plane
1474  }
1475  }
1476  return fCubicVolume;
1477 }
virtual G4double GetCubicVolume()
Definition: G4VSolid.cc:189
tuple x
Definition: test.py:50

Here is the call graph for this function:

G4GeometryType G4GenericTrap::GetEntityType ( ) const
virtual

Implements G4VSolid.

Definition at line 1296 of file G4GenericTrap.cc.

1297 {
1298  return G4String("G4GenericTrap");
1299 }

Here is the caller graph for this function:

G4VisExtent G4GenericTrap::GetExtent ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 2013 of file G4GenericTrap.cc.

2014 {
2015  // Computes bounding vectors for the shape
2016 
2017 #ifdef G4TESS_TEST
2018  if ( fTessellatedSolid )
2019  {
2020  return fTessellatedSolid->GetExtent();
2021  }
2022 #endif
2023 
2024  G4ThreeVector minVec = GetMinimumBBox();
2025  G4ThreeVector maxVec = GetMaximumBBox();
2026  return G4VisExtent (minVec.x(), maxVec.x(),
2027  minVec.y(), maxVec.y(),
2028  minVec.z(), maxVec.z());
2029 }
double x() const
virtual G4VisExtent GetExtent() const
double z() const
double y() const

Here is the call graph for this function:

G4int G4GenericTrap::GetNofVertices ( ) const
inline
G4ThreeVector G4GenericTrap::GetPointOnSurface ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 1333 of file G4GenericTrap.cc.

1334 {
1335 
1336 #ifdef G4TESS_TEST
1337  if ( fTessellatedSolid )
1338  {
1339  return fTessellatedSolid->GetPointOnSurface();
1340  }
1341 #endif
1342 
1343  G4ThreeVector point;
1344  G4TwoVector u,v,w;
1345  G4double rand,area,chose,cf,lambda0,lambda1,alfa,beta,zp;
1346  G4int ipl,j;
1347 
1348  std::vector<G4ThreeVector> vertices;
1349  for (G4int i=0; i<4;i++)
1350  {
1351  vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz));
1352  }
1353  for (G4int i=4; i<8;i++)
1354  {
1355  vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),fDz));
1356  }
1357 
1358  // Surface Area of Planes(only estimation for twisted)
1359  //
1360  G4double Surface0=GetFaceSurfaceArea(vertices[0],vertices[1],
1361  vertices[2],vertices[3]);//-fDz plane
1362  G4double Surface1=GetFaceSurfaceArea(vertices[0],vertices[1],
1363  vertices[5],vertices[4]);// Lat plane
1364  G4double Surface2=GetFaceSurfaceArea(vertices[3],vertices[0],
1365  vertices[4],vertices[7]);// Lat plane
1366  G4double Surface3=GetFaceSurfaceArea(vertices[2],vertices[3],
1367  vertices[7],vertices[6]);// Lat plane
1368  G4double Surface4=GetFaceSurfaceArea(vertices[2],vertices[1],
1369  vertices[5],vertices[6]);// Lat plane
1370  G4double Surface5=GetFaceSurfaceArea(vertices[4],vertices[5],
1371  vertices[6],vertices[7]);// fDz plane
1372  rand = G4UniformRand();
1373  area = Surface0+Surface1+Surface2+Surface3+Surface4+Surface5;
1374  chose = rand*area;
1375 
1376  if ( ( chose < Surface0)
1377  || ( chose > (Surface0+Surface1+Surface2+Surface3+Surface4)) )
1378  { // fDz or -fDz Plane
1379  ipl = G4int(G4UniformRand()*4);
1380  j = (ipl+1)%4;
1381  if(chose < Surface0)
1382  {
1383  zp = -fDz;
1384  u = fVertices[ipl]; v = fVertices[j];
1385  w = fVertices[(ipl+3)%4];
1386  }
1387  else
1388  {
1389  zp = fDz;
1390  u = fVertices[ipl+4]; v = fVertices[j+4];
1391  w = fVertices[(ipl+3)%4+4];
1392  }
1393  alfa = G4UniformRand();
1394  beta = G4UniformRand();
1395  lambda1=alfa*beta;
1396  lambda0=alfa-lambda1;
1397  v = v-u;
1398  w = w-u;
1399  v = u+lambda0*v+lambda1*w;
1400  }
1401  else // Lateral Plane Twisted or Not
1402  {
1403  if (chose < Surface0+Surface1) { ipl=0; }
1404  else if (chose < Surface0+Surface1+Surface2) { ipl=1; }
1405  else if (chose < Surface0+Surface1+Surface2+Surface3) { ipl=2; }
1406  else { ipl=3; }
1407  j = (ipl+1)%4;
1408  zp = -fDz+G4UniformRand()*2*fDz;
1409  cf = 0.5*(fDz-zp)/fDz;
1410  u = fVertices[ipl+4]+cf*( fVertices[ipl]-fVertices[ipl+4]);
1411  v = fVertices[j+4]+cf*(fVertices[j]-fVertices[j+4]);
1412  v = u+(v-u)*G4UniformRand();
1413  }
1414  point=G4ThreeVector(v.x(),v.y(),zp);
1415 
1416  return point;
1417 }
double y() const
double x() const
CLHEP::Hep3Vector G4ThreeVector
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
tuple v
Definition: test.py:18
double G4double
Definition: G4Types.hh:76
virtual G4ThreeVector GetPointOnSurface() const

Here is the call graph for this function:

G4Polyhedron * G4GenericTrap::GetPolyhedron ( ) const
virtual

Reimplemented from G4VSolid.

Definition at line 1972 of file G4GenericTrap.cc.

1973 {
1974 
1975 #ifdef G4TESS_TEST
1976  if ( fTessellatedSolid )
1977  {
1978  return fTessellatedSolid->GetPolyhedron();
1979  }
1980 #endif
1981 
1982  if ( (!fpPolyhedron)
1986  {
1987  G4AutoLock l(&polyhedronMutex);
1988  delete fpPolyhedron;
1990  fRebuildPolyhedron = false;
1991  l.unlock();
1992  }
1993  return fpPolyhedron;
1994 }
G4Polyhedron * fpPolyhedron
virtual G4Polyhedron * GetPolyhedron() const
G4Polyhedron * CreatePolyhedron() const
static G4int GetNumberOfRotationSteps()
G4int GetNumberOfRotationStepsAtTimeOfCreation() const
G4bool fRebuildPolyhedron

Here is the call graph for this function:

G4double G4GenericTrap::GetSurfaceArea ( )
virtual

Reimplemented from G4VSolid.

Definition at line 1421 of file G4GenericTrap.cc.

1422 {
1423  if (fSurfaceArea == 0.0) {
1424  if(fIsTwisted) {
1425  fSurfaceArea = G4VSolid::GetSurfaceArea();
1426  } else {
1427  // Set vertices
1428  G4ThreeVector vertix0(fVertices[0].x(),fVertices[0].y(),-fDz);
1429  G4ThreeVector vertix1(fVertices[1].x(),fVertices[1].y(),-fDz);
1430  G4ThreeVector vertix2(fVertices[2].x(),fVertices[2].y(),-fDz);
1431  G4ThreeVector vertix3(fVertices[3].x(),fVertices[3].y(),-fDz);
1432  G4ThreeVector vertix4(fVertices[4].x(),fVertices[4].y(), fDz);
1433  G4ThreeVector vertix5(fVertices[5].x(),fVertices[5].y(), fDz);
1434  G4ThreeVector vertix6(fVertices[6].x(),fVertices[6].y(), fDz);
1435  G4ThreeVector vertix7(fVertices[7].x(),fVertices[7].y(), fDz);
1436 
1437  // Find Surface Area
1438  fSurfaceArea = GetFaceSurfaceArea(vertix0,vertix1,vertix2,vertix3) // -fDz plane
1439  + GetFaceSurfaceArea(vertix1,vertix0,vertix4,vertix5) // Lat plane
1440  + GetFaceSurfaceArea(vertix2,vertix1,vertix5,vertix6) // Lat plane
1441  + GetFaceSurfaceArea(vertix3,vertix2,vertix6,vertix7) // Lat plane
1442  + GetFaceSurfaceArea(vertix0,vertix3,vertix7,vertix4) // Lat plane
1443  + GetFaceSurfaceArea(vertix7,vertix6,vertix5,vertix4); // +fDz plane
1444  }
1445  }
1446  return fSurfaceArea;
1447 }
tuple x
Definition: test.py:50
virtual G4double GetSurfaceArea()
Definition: G4VSolid.cc:251

Here is the call graph for this function:

G4double G4GenericTrap::GetTwistAngle ( G4int  index) const
inline

Here is the caller graph for this function:

G4TwoVector G4GenericTrap::GetVertex ( G4int  index) const
inline

Here is the caller graph for this function:

const std::vector<G4TwoVector>& G4GenericTrap::GetVertices ( ) const
inline

Here is the caller graph for this function:

G4int G4GenericTrap::GetVisSubdivisions ( ) const
inline

Here is the caller graph for this function:

G4double G4GenericTrap::GetZHalfLength ( ) const
inline

Here is the caller graph for this function:

EInside G4GenericTrap::Inside ( const G4ThreeVector p) const
virtual

Implements G4VSolid.

Definition at line 357 of file G4GenericTrap.cc.

358 {
359  // Test if point is inside this shape
360 
361 #ifdef G4TESS_TEST
362  if ( fTessellatedSolid )
363  {
364  return fTessellatedSolid->Inside(p);
365  }
366 #endif
367 
368  EInside innew=kOutside;
369  std::vector<G4TwoVector> xy;
370 
371  if (std::fabs(p.z()) <= fDz+halfCarTolerance) // First check Z range
372  {
373  // Compute intersection between Z plane containing point and the shape
374  //
375  G4double cf = 0.5*(fDz-p.z())/fDz;
376  for (G4int i=0; i<4; i++)
377  {
378  xy.push_back(fVertices[i+4]+cf*( fVertices[i]-fVertices[i+4]));
379  }
380 
381  innew=InsidePolygone(p,xy);
382 
383  if( (innew==kInside) || (innew==kSurface) )
384  {
385  if(std::fabs(p.z()) > fDz-halfCarTolerance) { innew=kSurface; }
386  }
387  }
388  return innew;
389 }
int G4int
Definition: G4Types.hh:78
double z() const
virtual EInside Inside(const G4ThreeVector &p) const
EInside
Definition: geomdefs.hh:58
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4bool G4GenericTrap::IsTwisted ( ) const
inline
G4GenericTrap & G4GenericTrap::operator= ( const G4GenericTrap rhs)

Definition at line 223 of file G4GenericTrap.cc.

224 {
225  // Check assignment to self
226  //
227  if (this == &rhs) { return *this; }
228 
229  // Copy base class data
230  //
231  G4VSolid::operator=(rhs);
232 
233  // Copy data
234  //
235  halfCarTolerance = rhs.halfCarTolerance;
236  fDz = rhs.fDz; fVertices = rhs.fVertices;
237  fIsTwisted = rhs.fIsTwisted; fTessellatedSolid = 0;
238  fMinBBoxVector = rhs.fMinBBoxVector; fMaxBBoxVector = rhs.fMaxBBoxVector;
239  fVisSubdivisions = rhs.fVisSubdivisions;
240  fSurfaceArea = rhs.fSurfaceArea; fCubicVolume = rhs.fCubicVolume;
241 
242  for (size_t i=0; i<4; ++i) { fTwist[i] = rhs.fTwist[i]; }
243 #ifdef G4TESS_TEST
244  if (rhs.fTessellatedSolid && !fIsTwisted )
245  { delete fTessellatedSolid; fTessellatedSolid = CreateTessellatedSolid(); }
246 #endif
247  fRebuildPolyhedron = false;
248  delete fpPolyhedron; fpPolyhedron = 0;
249 
250  return *this;
251 }
G4Polyhedron * fpPolyhedron
G4VSolid & operator=(const G4VSolid &rhs)
Definition: G4VSolid.cc:111
G4bool fRebuildPolyhedron

Here is the call graph for this function:

void G4GenericTrap::SetVisSubdivisions ( G4int  subdiv)
inline
std::ostream & G4GenericTrap::StreamInfo ( std::ostream &  os) const
virtual

Implements G4VSolid.

Definition at line 1310 of file G4GenericTrap.cc.

1311 {
1312  G4int oldprc = os.precision(16);
1313  os << "-----------------------------------------------------------\n"
1314  << " *** Dump for solid - " << GetName() << " *** \n"
1315  << " =================================================== \n"
1316  << " Solid geometry type: " << GetEntityType() << G4endl
1317  << " half length Z: " << fDz/mm << " mm \n"
1318  << " list of vertices:\n";
1319 
1320  for ( G4int i=0; i<fgkNofVertices; ++i )
1321  {
1322  os << std::setw(5) << "#" << i
1323  << " vx = " << fVertices[i].x()/mm << " mm"
1324  << " vy = " << fVertices[i].y()/mm << " mm" << G4endl;
1325  }
1326  os.precision(oldprc);
1327 
1328  return os;
1329 }
G4String GetName() const
static constexpr double mm
Definition: G4SIunits.hh:115
G4GeometryType GetEntityType() const
int G4int
Definition: G4Types.hh:78
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

G4ThreeVector G4GenericTrap::SurfaceNormal ( const G4ThreeVector p) const
virtual

Implements G4VSolid.

Definition at line 393 of file G4GenericTrap.cc.

394 {
395  // Calculate side nearest to p, and return normal
396  // If two sides are equidistant, sum of the Normal is returned
397 
398 #ifdef G4TESS_TEST
399  if ( fTessellatedSolid )
400  {
401  return fTessellatedSolid->SurfaceNormal(p);
402  }
403 #endif
404 
405  G4ThreeVector lnorm, sumnorm(0.,0.,0.), apprnorm(0.,0.,1.),
406  p0, p1, p2, r1, r2, r3, r4;
407  G4int noSurfaces = 0;
408  G4double distxy,distz;
409  G4bool zPlusSide=false;
410 
411  distz = fDz-std::fabs(p.z());
412  if (distz < halfCarTolerance)
413  {
414  if(p.z()>0)
415  {
416  zPlusSide=true;
417  sumnorm=G4ThreeVector(0,0,1);
418  }
419  else
420  {
421  sumnorm=G4ThreeVector(0,0,-1);
422  }
423  noSurfaces ++;
424  }
425 
426  // Check lateral planes
427  //
428  std:: vector<G4TwoVector> vertices;
429  G4double cf = 0.5*(fDz-p.z())/fDz;
430  for (G4int i=0; i<4; i++)
431  {
432  vertices.push_back(fVertices[i+4]+cf*(fVertices[i]-fVertices[i+4]));
433  }
434 
435  // Compute distance for lateral planes
436  //
437  for (G4int q=0; q<4; q++)
438  {
439  p0=G4ThreeVector(vertices[q].x(),vertices[q].y(),p.z());
440  if(zPlusSide)
441  {
442  p1=G4ThreeVector(fVertices[q].x(),fVertices[q].y(),-fDz);
443  }
444  else
445  {
446  p1=G4ThreeVector(fVertices[q+4].x(),fVertices[q+4].y(),fDz);
447  }
448  p2=G4ThreeVector(vertices[(q+1)%4].x(),vertices[(q+1)%4].y(),p.z());
449 
450  // Collapsed vertices
451  //
452  if ( (p2-p0).mag2() < kCarTolerance )
453  {
454  if ( std::fabs(p.z()+fDz) > kCarTolerance )
455  {
456  p2=G4ThreeVector(fVertices[(q+1)%4].x(),fVertices[(q+1)%4].y(),-fDz);
457  }
458  else
459  {
460  p2=G4ThreeVector(fVertices[(q+1)%4+4].x(),fVertices[(q+1)%4+4].y(),fDz);
461  }
462  }
463  lnorm = (p1-p0).cross(p2-p0);
464  lnorm = lnorm.unit();
465  if(zPlusSide) { lnorm=-lnorm; }
466 
467  // Adjust Normal for Twisted Surface
468  //
469  if ( (fIsTwisted) && (GetTwistAngle(q)!=0) )
470  {
471  G4double normP=(p2-p0).mag();
472  if(normP)
473  {
474  G4double proj=(p-p0).dot(p2-p0)/normP;
475  if(proj<0) { proj=0; }
476  if(proj>normP) { proj=normP; }
477  G4int j=(q+1)%4;
478  r1=G4ThreeVector(fVertices[q+4].x(),fVertices[q+4].y(),fDz);
479  r2=G4ThreeVector(fVertices[j+4].x(),fVertices[j+4].y(),fDz);
480  r3=G4ThreeVector(fVertices[q].x(),fVertices[q].y(),-fDz);
481  r4=G4ThreeVector(fVertices[j].x(),fVertices[j].y(),-fDz);
482  r1=r1+proj*(r2-r1)/normP;
483  r3=r3+proj*(r4-r3)/normP;
484  r2=r1-r3;
485  r4=r2.cross(p2-p0); r4=r4.unit();
486  lnorm=r4;
487  }
488  } // End if fIsTwisted
489 
490  distxy=std::fabs((p0-p).dot(lnorm));
491  if ( distxy<halfCarTolerance )
492  {
493  noSurfaces ++;
494 
495  // Negative sign for Normal is taken for Outside Normal
496  //
497  sumnorm=sumnorm+lnorm;
498  }
499 
500  // For ApproxSurfaceNormal
501  //
502  if (distxy<distz)
503  {
504  distz=distxy;
505  apprnorm=lnorm;
506  }
507  } // End for loop
508 
509  // Calculate final Normal, add Normal in the Corners and Touching Sides
510  //
511  if ( noSurfaces == 0 )
512  {
513 #ifdef G4SPECSDEBUG
514  G4Exception("G4GenericTrap::SurfaceNormal(p)", "GeomSolids1002",
515  JustWarning, "Point p is not on surface !?" );
516 #endif
517  sumnorm=apprnorm;
518  // Add Approximative Surface Normal Calculation?
519  }
520  else if ( noSurfaces == 1 ) { ; }
521  else { sumnorm = sumnorm.unit(); }
522 
523  return sumnorm ;
524 }
CLHEP::Hep3Vector G4ThreeVector
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
double z() const
bool G4bool
Definition: G4Types.hh:79
G4double GetTwistAngle(G4int index) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
Hep3Vector unit() const
G4double kCarTolerance
Definition: G4VSolid.hh:307
Hep3Vector cross(const Hep3Vector &) const
double G4double
Definition: G4Types.hh:76
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const

Here is the call graph for this function:

Member Data Documentation

G4Polyhedron* G4GenericTrap::fpPolyhedron
mutableprotected

Definition at line 204 of file G4GenericTrap.hh.

G4bool G4GenericTrap::fRebuildPolyhedron
mutableprotected

Definition at line 203 of file G4GenericTrap.hh.


The documentation for this class was generated from the following files: