#include <G4GenericTrap.hh>
Inheritance diagram for G4GenericTrap:
Definition at line 79 of file G4GenericTrap.hh.
G4GenericTrap::G4GenericTrap | ( | const G4String & | name, | |
G4double | halfZ, | |||
const std::vector< G4TwoVector > & | vertices | |||
) |
Definition at line 68 of file G4GenericTrap.cc.
References FatalErrorInArgument, G4endl, G4Exception(), JustWarning, and G4VSolid::kCarTolerance.
Referenced by Clone().
00070 : G4VSolid(name), 00071 fpPolyhedron(0), 00072 fDz(halfZ), 00073 fVertices(), 00074 fIsTwisted(false), 00075 fTessellatedSolid(0), 00076 fMinBBoxVector(G4ThreeVector(0,0,0)), 00077 fMaxBBoxVector(G4ThreeVector(0,0,0)), 00078 fVisSubdivisions(0), 00079 fSurfaceArea(0.), 00080 fCubicVolume(0.) 00081 00082 { 00083 // General constructor 00084 const G4double min_length=5*1.e-6; 00085 G4double length = 0.; 00086 G4int k=0; 00087 G4String errorDescription = "InvalidSetup in \" "; 00088 errorDescription += name; 00089 errorDescription += "\""; 00090 00091 // Check vertices size 00092 00093 if ( G4int(vertices.size()) != fgkNofVertices ) 00094 { 00095 G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids0002", 00096 FatalErrorInArgument, "Number of vertices != 8"); 00097 } 00098 00099 // Check dZ 00100 // 00101 if (halfZ < kCarTolerance) 00102 { 00103 G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids0002", 00104 FatalErrorInArgument, "dZ is too small or negative"); 00105 } 00106 00107 // Check Ordering and Copy vertices 00108 // 00109 if(CheckOrder(vertices)) 00110 { 00111 for (G4int i=0; i<fgkNofVertices; ++i) {fVertices.push_back(vertices[i]);} 00112 } 00113 else 00114 { 00115 for (G4int i=0; i <4; ++i) {fVertices.push_back(vertices[3-i]);} 00116 for (G4int i=0; i <4; ++i) {fVertices.push_back(vertices[7-i]);} 00117 } 00118 00119 // Check length of segments and Adjust 00120 // 00121 for (G4int j=0; j < 2; j++) 00122 { 00123 for (G4int i=1; i<4; ++i) 00124 { 00125 k = j*4+i; 00126 length = (fVertices[k]-fVertices[k-1]).mag(); 00127 if ( ( length < min_length) && ( length > kCarTolerance ) ) 00128 { 00129 std::ostringstream message; 00130 message << "Length segment is too small." << G4endl 00131 << "Distance between " << fVertices[k-1] << " and " 00132 << fVertices[k] << " is only " << length << " mm !"; 00133 G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids1001", 00134 JustWarning, message, "Vertices will be collapsed."); 00135 fVertices[k]=fVertices[k-1]; 00136 } 00137 } 00138 } 00139 00140 // Compute Twist 00141 // 00142 for( G4int i=0; i<4; i++) { fTwist[i]=0.; } 00143 fIsTwisted = ComputeIsTwisted(); 00144 00145 // Compute Bounding Box 00146 // 00147 ComputeBBox(); 00148 00149 // If not twisted - create tessellated solid 00150 // (an alternative implementation for testing) 00151 // 00152 #ifdef G4TESS_TEST 00153 if ( !fIsTwisted ) { fTessellatedSolid = CreateTessellatedSolid(); } 00154 #endif 00155 }
G4GenericTrap::~G4GenericTrap | ( | ) |
G4GenericTrap::G4GenericTrap | ( | __void__ & | ) |
Definition at line 159 of file G4GenericTrap.cc.
00160 : G4VSolid(a), 00161 fpPolyhedron(0), 00162 fDz(0.), 00163 fVertices(), 00164 fIsTwisted(false), 00165 fTessellatedSolid(0), 00166 fMinBBoxVector(G4ThreeVector(0,0,0)), 00167 fMaxBBoxVector(G4ThreeVector(0,0,0)), 00168 fVisSubdivisions(0), 00169 fSurfaceArea(0.), 00170 fCubicVolume(0.) 00171 { 00172 // Fake default constructor - sets only member data and allocates memory 00173 // for usage restricted to object persistency. 00174 00175 for (size_t i=0; i<4; ++i) { fTwist[i]=0.; } 00176 }
G4GenericTrap::G4GenericTrap | ( | const G4GenericTrap & | rhs | ) |
Definition at line 188 of file G4GenericTrap.cc.
References fTessellatedSolid, and fTwist.
00189 : G4VSolid(rhs), 00190 fpPolyhedron(0), fDz(rhs.fDz), fVertices(rhs.fVertices), 00191 fIsTwisted(rhs.fIsTwisted), fTessellatedSolid(0), 00192 fMinBBoxVector(rhs.fMinBBoxVector), fMaxBBoxVector(rhs.fMaxBBoxVector), 00193 fVisSubdivisions(rhs.fVisSubdivisions), 00194 fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume) 00195 { 00196 for (size_t i=0; i<4; ++i) { fTwist[i] = rhs.fTwist[i]; } 00197 #ifdef G4TESS_TEST 00198 if (rhs.fTessellatedSolid && !fIsTwisted ) 00199 { fTessellatedSolid = CreateTessellatedSolid(); } 00200 #endif 00201 }
G4bool G4GenericTrap::CalculateExtent | ( | const EAxis | pAxis, | |
const G4VoxelLimits & | pVoxelLimit, | |||
const G4AffineTransform & | pTransform, | |||
G4double & | pmin, | |||
G4double & | pmax | |||
) | const [virtual] |
Implements G4VSolid.
Definition at line 1169 of file G4GenericTrap.cc.
References G4TessellatedSolid::CalculateExtent(), G4VSolid::ClipBetweenSections(), G4VSolid::ClipCrossSection(), G4VoxelLimits::GetMaxExtent(), G4VoxelLimits::GetMaxXExtent(), G4VoxelLimits::GetMaxYExtent(), G4VoxelLimits::GetMaxZExtent(), G4VoxelLimits::GetMinExtent(), G4VoxelLimits::GetMinXExtent(), G4VoxelLimits::GetMinYExtent(), G4VoxelLimits::GetMinZExtent(), Inside(), G4AffineTransform::Inverse(), G4AffineTransform::IsRotated(), G4VoxelLimits::IsXLimited(), G4VoxelLimits::IsYLimited(), G4VoxelLimits::IsZLimited(), G4VSolid::kCarTolerance, kOutside, kXAxis, kYAxis, kZAxis, G4AffineTransform::NetTranslation(), and G4AffineTransform::TransformPoint().
01173 { 01174 #ifdef G4TESS_TEST 01175 if ( fTessellatedSolid ) 01176 { 01177 return fTessellatedSolid->CalculateExtent(pAxis, pVoxelLimit, 01178 pTransform, pMin, pMax); 01179 } 01180 #endif 01181 01182 // Computes bounding vectors for a shape 01183 // 01184 G4double Dx,Dy; 01185 G4ThreeVector minVec = GetMinimumBBox(); 01186 G4ThreeVector maxVec = GetMaximumBBox(); 01187 Dx = 0.5*(maxVec.x()- minVec.x()); 01188 Dy = 0.5*(maxVec.y()- minVec.y()); 01189 01190 if (!pTransform.IsRotated()) 01191 { 01192 // Special case handling for unrotated shapes 01193 // Compute x/y/z mins and maxs respecting limits, with early returns 01194 // if outside limits. Then switch() on pAxis 01195 // 01196 G4double xoffset,xMin,xMax; 01197 G4double yoffset,yMin,yMax; 01198 G4double zoffset,zMin,zMax; 01199 01200 xoffset=pTransform.NetTranslation().x(); 01201 xMin=xoffset-Dx; 01202 xMax=xoffset+Dx; 01203 if (pVoxelLimit.IsXLimited()) 01204 { 01205 if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance) 01206 || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) ) 01207 { 01208 return false; 01209 } 01210 else 01211 { 01212 if (xMin<pVoxelLimit.GetMinXExtent()) 01213 { 01214 xMin=pVoxelLimit.GetMinXExtent(); 01215 } 01216 if (xMax>pVoxelLimit.GetMaxXExtent()) 01217 { 01218 xMax=pVoxelLimit.GetMaxXExtent(); 01219 } 01220 } 01221 } 01222 01223 yoffset=pTransform.NetTranslation().y(); 01224 yMin=yoffset-Dy; 01225 yMax=yoffset+Dy; 01226 if (pVoxelLimit.IsYLimited()) 01227 { 01228 if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance) 01229 || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) ) 01230 { 01231 return false; 01232 } 01233 else 01234 { 01235 if (yMin<pVoxelLimit.GetMinYExtent()) 01236 { 01237 yMin=pVoxelLimit.GetMinYExtent(); 01238 } 01239 if (yMax>pVoxelLimit.GetMaxYExtent()) 01240 { 01241 yMax=pVoxelLimit.GetMaxYExtent(); 01242 } 01243 } 01244 } 01245 01246 zoffset=pTransform.NetTranslation().z(); 01247 zMin=zoffset-fDz; 01248 zMax=zoffset+fDz; 01249 if (pVoxelLimit.IsZLimited()) 01250 { 01251 if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance) 01252 || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) ) 01253 { 01254 return false; 01255 } 01256 else 01257 { 01258 if (zMin<pVoxelLimit.GetMinZExtent()) 01259 { 01260 zMin=pVoxelLimit.GetMinZExtent(); 01261 } 01262 if (zMax>pVoxelLimit.GetMaxZExtent()) 01263 { 01264 zMax=pVoxelLimit.GetMaxZExtent(); 01265 } 01266 } 01267 } 01268 01269 switch (pAxis) 01270 { 01271 case kXAxis: 01272 pMin = xMin; 01273 pMax = xMax; 01274 break; 01275 case kYAxis: 01276 pMin = yMin; 01277 pMax = yMax; 01278 break; 01279 case kZAxis: 01280 pMin = zMin; 01281 pMax = zMax; 01282 break; 01283 default: 01284 break; 01285 } 01286 pMin-=kCarTolerance; 01287 pMax+=kCarTolerance; 01288 01289 return true; 01290 } 01291 else 01292 { 01293 // General rotated case - create and clip mesh to boundaries 01294 01295 G4bool existsAfterClip=false; 01296 G4ThreeVectorList *vertices; 01297 01298 pMin=+kInfinity; 01299 pMax=-kInfinity; 01300 01301 // Calculate rotated vertex coordinates 01302 // 01303 vertices=CreateRotatedVertices(pTransform); 01304 ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax); 01305 ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax); 01306 ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax); 01307 01308 if ( (pMin!=kInfinity) || (pMax!=-kInfinity) ) 01309 { 01310 existsAfterClip=true; 01311 01312 // Add 2*tolerance to avoid precision troubles 01313 // 01314 pMin-=kCarTolerance; 01315 pMax+=kCarTolerance; 01316 } 01317 else 01318 { 01319 // Check for case where completely enveloping clipping volume. 01320 // If point inside then we are confident that the solid completely 01321 // envelopes the clipping volume. Hence set min/max extents according 01322 // to clipping volume extents along the specified axis. 01323 // 01324 G4ThreeVector clipCentre( 01325 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5, 01326 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5, 01327 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5); 01328 01329 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside) 01330 { 01331 existsAfterClip=true; 01332 pMin=pVoxelLimit.GetMinExtent(pAxis); 01333 pMax=pVoxelLimit.GetMaxExtent(pAxis); 01334 } 01335 } 01336 delete vertices; 01337 return existsAfterClip; 01338 } 01339 }
G4VSolid * G4GenericTrap::Clone | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 1396 of file G4GenericTrap.cc.
References G4GenericTrap().
01397 { 01398 return new G4GenericTrap(*this); 01399 }
G4NURBS * G4GenericTrap::CreateNURBS | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 2219 of file G4GenericTrap.cc.
References G4TessellatedSolid::CreateNURBS().
02220 { 02221 #ifdef G4TESS_TEST 02222 if ( fTessellatedSolid ) 02223 { 02224 return fTessellatedSolid->CreateNURBS(); 02225 } 02226 #endif 02227 02228 // Computes bounding vectors for the shape 02229 // 02230 G4double Dx,Dy; 02231 G4ThreeVector minVec = GetMinimumBBox(); 02232 G4ThreeVector maxVec = GetMaximumBBox(); 02233 Dx = 0.5*(maxVec.x()- minVec.y()); 02234 Dy = 0.5*(maxVec.y()- minVec.y()); 02235 02236 return new G4NURBSbox (Dx, Dy, fDz); 02237 }
G4Polyhedron * G4GenericTrap::CreatePolyhedron | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 2123 of file G4GenericTrap.cc.
References G4PolyhedronArbitrary::AddFacet(), G4PolyhedronArbitrary::AddVertex(), G4TessellatedSolid::CreatePolyhedron(), GetTwistAngle(), GetVisSubdivisions(), G4PolyhedronArbitrary::InvertFacets(), and G4PolyhedronArbitrary::SetReferences().
Referenced by GetPolyhedron().
02124 { 02125 02126 #ifdef G4TESS_TEST 02127 if ( fTessellatedSolid ) 02128 { 02129 return fTessellatedSolid->CreatePolyhedron(); 02130 } 02131 #endif 02132 02133 // Approximation of Twisted Side 02134 // Construct extra Points, if Twisted Side 02135 // 02136 G4PolyhedronArbitrary* polyhedron; 02137 size_t nVertices, nFacets; 02138 02139 G4int subdivisions=0; 02140 G4int i; 02141 if(fIsTwisted) 02142 { 02143 if ( GetVisSubdivisions()!= 0 ) 02144 { 02145 subdivisions=GetVisSubdivisions(); 02146 } 02147 else 02148 { 02149 // Estimation of Number of Subdivisions for smooth visualisation 02150 // 02151 G4double maxTwist=0.; 02152 for(i=0; i<4; i++) 02153 { 02154 if(GetTwistAngle(i)>maxTwist) { maxTwist=GetTwistAngle(i); } 02155 } 02156 02157 // Computes bounding vectors for the shape 02158 // 02159 G4double Dx,Dy; 02160 G4ThreeVector minVec = GetMinimumBBox(); 02161 G4ThreeVector maxVec = GetMaximumBBox(); 02162 Dx = 0.5*(maxVec.x()- minVec.y()); 02163 Dy = 0.5*(maxVec.y()- minVec.y()); 02164 if (Dy > Dx) { Dx=Dy; } 02165 02166 subdivisions=8*G4int(maxTwist/(Dx*Dx*Dx)*fDz); 02167 if (subdivisions<4) { subdivisions=4; } 02168 if (subdivisions>30) { subdivisions=30; } 02169 } 02170 } 02171 G4int sub4=4*subdivisions; 02172 nVertices = 8+subdivisions*4; 02173 nFacets = 6+subdivisions*4; 02174 G4double cf=1./(subdivisions+1); 02175 polyhedron = new G4PolyhedronArbitrary (nVertices, nFacets); 02176 02177 // Add Vertex 02178 // 02179 for (i=0;i<4;i++) 02180 { 02181 polyhedron->AddVertex(G4ThreeVector(fVertices[i].x(), 02182 fVertices[i].y(),-fDz)); 02183 } 02184 for( i=0;i<subdivisions;i++) 02185 { 02186 for(G4int j=0;j<4;j++) 02187 { 02188 G4TwoVector u=fVertices[j]+cf*(i+1)*( fVertices[j+4]-fVertices[j]); 02189 polyhedron->AddVertex(G4ThreeVector(u.x(),u.y(),-fDz+cf*2*fDz*(i+1))); 02190 } 02191 } 02192 for (i=4;i<8;i++) 02193 { 02194 polyhedron->AddVertex(G4ThreeVector(fVertices[i].x(), 02195 fVertices[i].y(),fDz)); 02196 } 02197 02198 // Add Facets 02199 // 02200 polyhedron->AddFacet(1,4,3,2); //Z-plane 02201 for (i=0;i<subdivisions+1;i++) 02202 { 02203 G4int is=i*4; 02204 polyhedron->AddFacet(5+is,8+is,4+is,1+is); 02205 polyhedron->AddFacet(8+is,7+is,3+is,4+is); 02206 polyhedron->AddFacet(7+is,6+is,2+is,3+is); 02207 polyhedron->AddFacet(6+is,5+is,1+is,2+is); 02208 } 02209 polyhedron->AddFacet(5+sub4,6+sub4,7+sub4,8+sub4); //Z-plane 02210 02211 polyhedron->SetReferences(); 02212 polyhedron->InvertFacets(); 02213 02214 return (G4Polyhedron*) polyhedron; 02215 }
void G4GenericTrap::DescribeYourselfTo | ( | G4VGraphicsScene & | scene | ) | const [virtual] |
Implements G4VSolid.
Definition at line 2086 of file G4GenericTrap.cc.
References G4VGraphicsScene::AddSolid(), and G4TessellatedSolid::DescribeYourselfTo().
02087 { 02088 02089 #ifdef G4TESS_TEST 02090 if ( fTessellatedSolid ) 02091 { 02092 return fTessellatedSolid->DescribeYourselfTo(scene); 02093 } 02094 #endif 02095 02096 scene.AddSolid(*this); 02097 }
G4double G4GenericTrap::DistanceToIn | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 789 of file G4GenericTrap.cc.
References G4TessellatedSolid::DistanceToIn().
00790 { 00791 // Computes the closest distance from given point to this shape 00792 00793 #ifdef G4TESS_TEST 00794 if ( fTessellatedSolid ) 00795 { 00796 return fTessellatedSolid->DistanceToIn(p); 00797 } 00798 #endif 00799 00800 G4double safz = std::fabs(p.z())-fDz; 00801 if(safz<0) { safz=0; } 00802 00803 G4int iseg; 00804 G4double safe = safz; 00805 G4double safxy = safz; 00806 00807 for (iseg=0; iseg<4; iseg++) 00808 { 00809 safxy = SafetyToFace(p,iseg); 00810 if (safxy>safe) { safe=safxy; } 00811 } 00812 00813 return safe; 00814 }
G4double G4GenericTrap::DistanceToIn | ( | const G4ThreeVector & | p, | |
const G4ThreeVector & | v | |||
) | const [virtual] |
Implements G4VSolid.
Definition at line 719 of file G4GenericTrap.cc.
References G4TessellatedSolid::DistanceToIn(), Inside(), G4VSolid::kCarTolerance, kOutside, and CLHEP::detail::n.
00721 { 00722 #ifdef G4TESS_TEST 00723 if ( fTessellatedSolid ) 00724 { 00725 return fTessellatedSolid->DistanceToIn(p, v); 00726 } 00727 #endif 00728 00729 static const G4double halfCarTolerance=kCarTolerance*0.5; 00730 00731 G4double dist[5]; 00732 G4ThreeVector n; 00733 00734 // Check lateral faces 00735 // 00736 G4int i; 00737 for (i=0; i<4; i++) 00738 { 00739 dist[i]=DistToPlane(p, v, i); 00740 } 00741 00742 // Check Z planes 00743 // 00744 dist[4]=kInfinity; 00745 if (std::fabs(p.z())>fDz-halfCarTolerance) 00746 { 00747 if (v.z()) 00748 { 00749 G4ThreeVector pt; 00750 if (p.z()>0) 00751 { 00752 dist[4] = (fDz-p.z())/v.z(); 00753 } 00754 else 00755 { 00756 dist[4] = (-fDz-p.z())/v.z(); 00757 } 00758 if (dist[4]<-halfCarTolerance) 00759 { 00760 dist[4]=kInfinity; 00761 } 00762 else 00763 { 00764 if(dist[4]<halfCarTolerance) 00765 { 00766 if(p.z()>0) { n=G4ThreeVector(0,0,1); } 00767 else { n=G4ThreeVector(0,0,-1); } 00768 if (n.dot(v)<0) { dist[4]=0.; } 00769 else { dist[4]=kInfinity; } 00770 } 00771 pt=p+dist[4]*v; 00772 if (Inside(pt)==kOutside) { dist[4]=kInfinity; } 00773 } 00774 } 00775 } 00776 G4double distmin = dist[0]; 00777 for (i=1;i<5;i++) 00778 { 00779 if (dist[i] < distmin) { distmin = dist[i]; } 00780 } 00781 00782 if (distmin<halfCarTolerance) { distmin=0.; } 00783 00784 return distmin; 00785 }
G4double G4GenericTrap::DistanceToOut | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 1142 of file G4GenericTrap.cc.
References G4TessellatedSolid::DistanceToOut().
01143 { 01144 01145 #ifdef G4TESS_TEST 01146 if ( fTessellatedSolid ) 01147 { 01148 return fTessellatedSolid->DistanceToOut(p); 01149 } 01150 #endif 01151 01152 G4double safz = fDz-std::fabs(p.z()); 01153 if (safz<0) { safz = 0; } 01154 01155 G4double safe = safz; 01156 G4double safxy = safz; 01157 01158 for (G4int iseg=0; iseg<4; iseg++) 01159 { 01160 safxy = std::fabs(SafetyToFace(p,iseg)); 01161 if (safxy < safe) { safe = safxy; } 01162 } 01163 01164 return safe; 01165 }
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 895 of file G4GenericTrap.cc.
References G4TessellatedSolid::DistanceToOut(), G4VSolid::DumpInfo(), G4endl, G4Exception(), JustWarning, G4VSolid::kCarTolerance, and kOutside.
00900 { 00901 #ifdef G4TESS_TEST 00902 if ( fTessellatedSolid ) 00903 { 00904 return fTessellatedSolid->DistanceToOut(p, v, calcNorm, validNorm, n); 00905 } 00906 #endif 00907 00908 static const G4double halfCarTolerance=kCarTolerance*0.5; 00909 00910 G4double distmin; 00911 G4bool lateral_cross = false; 00912 ESide side = kUndefined; 00913 00914 if (calcNorm) { *validNorm=true; } // All normals are valid 00915 00916 if (v.z() < 0) 00917 { 00918 distmin=(-fDz-p.z())/v.z(); 00919 if (calcNorm) { side=kMZ; *n=G4ThreeVector(0,0,-1); } 00920 } 00921 else 00922 { 00923 if (v.z() > 0) 00924 { 00925 distmin = (fDz-p.z())/v.z(); 00926 if (calcNorm) { side=kPZ; *n=G4ThreeVector(0,0,1); } 00927 } 00928 else { distmin = kInfinity; } 00929 } 00930 00931 G4double dz2 =0.5/fDz; 00932 G4double xa,xb,xc,xd; 00933 G4double ya,yb,yc,yd; 00934 00935 for (G4int ipl=0; ipl<4; ipl++) 00936 { 00937 G4int j = (ipl+1)%4; 00938 xa=fVertices[ipl].x(); 00939 ya=fVertices[ipl].y(); 00940 xb=fVertices[ipl+4].x(); 00941 yb=fVertices[ipl+4].y(); 00942 xc=fVertices[j].x(); 00943 yc=fVertices[j].y(); 00944 xd=fVertices[4+j].x(); 00945 yd=fVertices[4+j].y(); 00946 00947 if ( ((std::fabs(xb-xd)+std::fabs(yb-yd))<halfCarTolerance) 00948 || ((std::fabs(xa-xc)+std::fabs(ya-yc))<halfCarTolerance) ) 00949 { 00950 G4double q=DistToTriangle(p,v,ipl) ; 00951 if ( (q>=0) && (q<distmin) ) 00952 { 00953 distmin=q; 00954 lateral_cross=true; 00955 side=ESide(ipl+1); 00956 } 00957 continue; 00958 } 00959 G4double tx1 =dz2*(xb-xa); 00960 G4double ty1 =dz2*(yb-ya); 00961 G4double tx2 =dz2*(xd-xc); 00962 G4double ty2 =dz2*(yd-yc); 00963 G4double dzp =fDz+p.z(); 00964 G4double xs1 =xa+tx1*dzp; 00965 G4double ys1 =ya+ty1*dzp; 00966 G4double xs2 =xc+tx2*dzp; 00967 G4double ys2 =yc+ty2*dzp; 00968 G4double dxs =xs2-xs1; 00969 G4double dys =ys2-ys1; 00970 G4double dtx =tx2-tx1; 00971 G4double dty =ty2-ty1; 00972 G4double a = (dtx*v.y()-dty*v.x()+(tx1*ty2-tx2*ty1)*v.z())*v.z(); 00973 G4double b = dxs*v.y()-dys*v.x()+(dtx*p.y()-dty*p.x()+ty2*xs1-ty1*xs2 00974 + tx1*ys2-tx2*ys1)*v.z(); 00975 G4double c=dxs*p.y()-dys*p.x()+xs1*ys2-xs2*ys1; 00976 G4double q=kInfinity; 00977 00978 if (std::fabs(a) < kCarTolerance) 00979 { 00980 if (std::fabs(b) < kCarTolerance) { continue; } 00981 q=-c/b; 00982 00983 // Check for Point on the Surface 00984 // 00985 if ((q > -halfCarTolerance) && (q < distmin)) 00986 { 00987 if (q < halfCarTolerance) 00988 { 00989 if (NormalToPlane(p,ipl).dot(v)<0.) { continue; } 00990 } 00991 distmin =q; 00992 lateral_cross=true; 00993 side=ESide(ipl+1); 00994 } 00995 continue; 00996 } 00997 G4double d=b*b-4*a*c; 00998 if (d >= 0.) 00999 { 01000 if (a > 0) { q=0.5*(-b-std::sqrt(d))/a; } 01001 else { q=0.5*(-b+std::sqrt(d))/a; } 01002 01003 // Check for Point on the Surface 01004 // 01005 if (q > -halfCarTolerance ) 01006 { 01007 if (q < distmin) 01008 { 01009 if(q < halfCarTolerance) 01010 { 01011 if (NormalToPlane(p,ipl).dot(v)<0.) // Check second root 01012 { 01013 if (a > 0) { q=0.5*(-b+std::sqrt(d))/a; } 01014 else { q=0.5*(-b-std::sqrt(d))/a; } 01015 if (( q > halfCarTolerance) && (q < distmin)) 01016 { 01017 distmin=q; 01018 lateral_cross = true; 01019 side=ESide(ipl+1); 01020 } 01021 continue; 01022 } 01023 } 01024 distmin = q; 01025 lateral_cross = true; 01026 side=ESide(ipl+1); 01027 } 01028 } 01029 else 01030 { 01031 if (a > 0) { q=0.5*(-b+std::sqrt(d))/a; } 01032 else { q=0.5*(-b-std::sqrt(d))/a; } 01033 01034 // Check for Point on the Surface 01035 // 01036 if ((q > -halfCarTolerance) && (q < distmin)) 01037 { 01038 if (q < halfCarTolerance) 01039 { 01040 if (NormalToPlane(p,ipl).dot(v)<0.) // Check second root 01041 { 01042 if (a > 0) { q=0.5*(-b-std::sqrt(d))/a; } 01043 else { q=0.5*(-b+std::sqrt(d))/a; } 01044 if ( ( q > halfCarTolerance) && (q < distmin) ) 01045 { 01046 distmin=q; 01047 lateral_cross = true; 01048 side=ESide(ipl+1); 01049 } 01050 continue; 01051 } 01052 } 01053 distmin =q; 01054 lateral_cross = true; 01055 side=ESide(ipl+1); 01056 } 01057 } 01058 } 01059 } 01060 if (!lateral_cross) // Make sure that track crosses the top or bottom 01061 { 01062 if (distmin >= kInfinity) { distmin=kCarTolerance; } 01063 G4ThreeVector pt=p+distmin*v; 01064 01065 // Check if propagated point is in the polygon 01066 // 01067 G4int i=0; 01068 if (v.z()>0.) { i=4; } 01069 std::vector<G4TwoVector> xy; 01070 for ( G4int j=0; j<4; j++) { xy.push_back(fVertices[i+j]); } 01071 01072 // Check Inside 01073 // 01074 if (InsidePolygone(pt,xy)==kOutside) 01075 { 01076 if(calcNorm) 01077 { 01078 if (v.z()>0) {side= kPZ; *n = G4ThreeVector(0,0,1);} 01079 else { side=kMZ; *n = G4ThreeVector(0,0,-1);} 01080 } 01081 return 0.; 01082 } 01083 else 01084 { 01085 if(v.z()>0) {side=kPZ;} 01086 else {side=kMZ;} 01087 } 01088 } 01089 01090 if (calcNorm) 01091 { 01092 G4ThreeVector pt=p+v*distmin; 01093 switch (side) 01094 { 01095 case kXY0: 01096 *n=NormalToPlane(pt,0); 01097 break; 01098 case kXY1: 01099 *n=NormalToPlane(pt,1); 01100 break; 01101 case kXY2: 01102 *n=NormalToPlane(pt,2); 01103 break; 01104 case kXY3: 01105 *n=NormalToPlane(pt,3); 01106 break; 01107 case kPZ: 01108 *n=G4ThreeVector(0,0,1); 01109 break; 01110 case kMZ: 01111 *n=G4ThreeVector(0,0,-1); 01112 break; 01113 default: 01114 DumpInfo(); 01115 std::ostringstream message; 01116 G4int oldprc = message.precision(16); 01117 message << "Undefined side for valid surface normal to solid." << G4endl 01118 << "Position:" << G4endl 01119 << " p.x() = " << p.x()/mm << " mm" << G4endl 01120 << " p.y() = " << p.y()/mm << " mm" << G4endl 01121 << " p.z() = " << p.z()/mm << " mm" << G4endl 01122 << "Direction:" << G4endl 01123 << " v.x() = " << v.x() << G4endl 01124 << " v.y() = " << v.y() << G4endl 01125 << " v.z() = " << v.z() << G4endl 01126 << "Proposed distance :" << G4endl 01127 << " distmin = " << distmin/mm << " mm"; 01128 message.precision(oldprc); 01129 G4Exception("G4GenericTrap::DistanceToOut(p,v,..)", 01130 "GeomSolids1002", JustWarning, message); 01131 break; 01132 } 01133 } 01134 01135 if (distmin<halfCarTolerance) { distmin=0.; } 01136 01137 return distmin; 01138 }
G4double G4GenericTrap::GetCubicVolume | ( | ) | [virtual] |
Reimplemented from G4VSolid.
Definition at line 1514 of file G4GenericTrap.cc.
References G4VSolid::GetCubicVolume().
01515 { 01516 if(fCubicVolume != 0.) {;} 01517 else { fCubicVolume = G4VSolid::GetCubicVolume(); } 01518 return fCubicVolume; 01519 }
G4GeometryType G4GenericTrap::GetEntityType | ( | ) | const [virtual] |
Implements G4VSolid.
Definition at line 1389 of file G4GenericTrap.cc.
Referenced by StreamInfo().
01390 { 01391 return G4String("G4GenericTrap"); 01392 }
G4VisExtent G4GenericTrap::GetExtent | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 2101 of file G4GenericTrap.cc.
References G4TessellatedSolid::GetExtent().
02102 { 02103 // Computes bounding vectors for the shape 02104 02105 #ifdef G4TESS_TEST 02106 if ( fTessellatedSolid ) 02107 { 02108 return fTessellatedSolid->GetExtent(); 02109 } 02110 #endif 02111 02112 G4double Dx,Dy; 02113 G4ThreeVector minVec = GetMinimumBBox(); 02114 G4ThreeVector maxVec = GetMaximumBBox(); 02115 Dx = 0.5*(maxVec.x()- minVec.x()); 02116 Dy = 0.5*(maxVec.y()- minVec.y()); 02117 02118 return G4VisExtent (-Dx, Dx, -Dy, Dy, -fDz, fDz); 02119 }
G4int G4GenericTrap::GetNofVertices | ( | ) | const [inline] |
G4ThreeVector G4GenericTrap::GetPointOnSurface | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 1426 of file G4GenericTrap.cc.
References G4UniformRand, and G4TessellatedSolid::GetPointOnSurface().
01427 { 01428 01429 #ifdef G4TESS_TEST 01430 if ( fTessellatedSolid ) 01431 { 01432 return fTessellatedSolid->GetPointOnSurface(); 01433 } 01434 #endif 01435 01436 G4ThreeVector point; 01437 G4TwoVector u,v,w; 01438 G4double rand,area,chose,cf,lambda0,lambda1,alfa,beta,zp; 01439 G4int ipl,j; 01440 01441 std::vector<G4ThreeVector> vertices; 01442 for (G4int i=0; i<4;i++) 01443 { 01444 vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz)); 01445 } 01446 for (G4int i=4; i<8;i++) 01447 { 01448 vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),fDz)); 01449 } 01450 01451 // Surface Area of Planes(only estimation for twisted) 01452 // 01453 G4double Surface0=GetFaceSurfaceArea(vertices[0],vertices[1], 01454 vertices[2],vertices[3]);//-fDz plane 01455 G4double Surface1=GetFaceSurfaceArea(vertices[0],vertices[1], 01456 vertices[5],vertices[4]);// Lat plane 01457 G4double Surface2=GetFaceSurfaceArea(vertices[3],vertices[0], 01458 vertices[4],vertices[7]);// Lat plane 01459 G4double Surface3=GetFaceSurfaceArea(vertices[2],vertices[3], 01460 vertices[7],vertices[6]);// Lat plane 01461 G4double Surface4=GetFaceSurfaceArea(vertices[2],vertices[1], 01462 vertices[5],vertices[6]);// Lat plane 01463 G4double Surface5=GetFaceSurfaceArea(vertices[4],vertices[5], 01464 vertices[6],vertices[7]);// fDz plane 01465 rand = G4UniformRand(); 01466 area = Surface0+Surface1+Surface2+Surface3+Surface4+Surface5; 01467 chose = rand*area; 01468 01469 if ( ( chose < Surface0) 01470 || ( chose > (Surface0+Surface1+Surface2+Surface3+Surface4)) ) 01471 { // fDz or -fDz Plane 01472 ipl = G4int(G4UniformRand()*4); 01473 j = (ipl+1)%4; 01474 if(chose < Surface0) 01475 { 01476 zp = -fDz; 01477 u = fVertices[ipl]; v = fVertices[j]; 01478 w = fVertices[(ipl+3)%4]; 01479 } 01480 else 01481 { 01482 zp = fDz; 01483 u = fVertices[ipl+4]; v = fVertices[j+4]; 01484 w = fVertices[(ipl+3)%4+4]; 01485 } 01486 alfa = G4UniformRand(); 01487 beta = G4UniformRand(); 01488 lambda1=alfa*beta; 01489 lambda0=alfa-lambda1; 01490 v = v-u; 01491 w = w-u; 01492 v = u+lambda0*v+lambda1*w; 01493 } 01494 else // Lateral Plane Twisted or Not 01495 { 01496 if (chose < Surface0+Surface1) { ipl=0; } 01497 else if (chose < Surface0+Surface1+Surface2) { ipl=1; } 01498 else if (chose < Surface0+Surface1+Surface2+Surface3) { ipl=2; } 01499 else { ipl=3; } 01500 j = (ipl+1)%4; 01501 zp = -fDz+G4UniformRand()*2*fDz; 01502 cf = 0.5*(fDz-zp)/fDz; 01503 u = fVertices[ipl+4]+cf*( fVertices[ipl]-fVertices[ipl+4]); 01504 v = fVertices[j+4]+cf*(fVertices[j]-fVertices[j+4]); 01505 v = u+(v-u)*G4UniformRand(); 01506 } 01507 point=G4ThreeVector(v.x(),v.y(),zp); 01508 01509 return point; 01510 }
G4Polyhedron * G4GenericTrap::GetPolyhedron | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 2064 of file G4GenericTrap.cc.
References CreatePolyhedron(), fpPolyhedron, G4Polyhedron::GetNumberOfRotationStepsAtTimeOfCreation(), and G4TessellatedSolid::GetPolyhedron().
02065 { 02066 02067 #ifdef G4TESS_TEST 02068 if ( fTessellatedSolid ) 02069 { 02070 return fTessellatedSolid->GetPolyhedron(); 02071 } 02072 #endif 02073 02074 if ( (!fpPolyhedron) 02075 || (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() != 02076 fpPolyhedron->GetNumberOfRotationSteps()) ) 02077 { 02078 delete fpPolyhedron; 02079 fpPolyhedron = CreatePolyhedron(); 02080 } 02081 return fpPolyhedron; 02082 }
G4double G4GenericTrap::GetSurfaceArea | ( | ) | [virtual] |
Reimplemented from G4VSolid.
Definition at line 1523 of file G4GenericTrap.cc.
References G4VSolid::GetSurfaceArea().
01524 { 01525 if(fSurfaceArea != 0.) {;} 01526 else 01527 { 01528 std::vector<G4ThreeVector> vertices; 01529 for (G4int i=0; i<4;i++) 01530 { 01531 vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz)); 01532 } 01533 for (G4int i=4; i<8;i++) 01534 { 01535 vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),fDz)); 01536 } 01537 01538 // Surface Area of Planes(only estimation for twisted) 01539 // 01540 G4double fSurface0=GetFaceSurfaceArea(vertices[0],vertices[1], 01541 vertices[2],vertices[3]);//-fDz plane 01542 G4double fSurface1=GetFaceSurfaceArea(vertices[0],vertices[1], 01543 vertices[5],vertices[4]);// Lat plane 01544 G4double fSurface2=GetFaceSurfaceArea(vertices[3],vertices[0], 01545 vertices[4],vertices[7]);// Lat plane 01546 G4double fSurface3=GetFaceSurfaceArea(vertices[2],vertices[3], 01547 vertices[7],vertices[6]);// Lat plane 01548 G4double fSurface4=GetFaceSurfaceArea(vertices[2],vertices[1], 01549 vertices[5],vertices[6]);// Lat plane 01550 G4double fSurface5=GetFaceSurfaceArea(vertices[4],vertices[5], 01551 vertices[6],vertices[7]);// fDz plane 01552 01553 // Total Surface Area 01554 // 01555 if(!fIsTwisted) 01556 { 01557 fSurfaceArea = fSurface0+fSurface1+fSurface2 01558 + fSurface3+fSurface4+fSurface5; 01559 } 01560 else 01561 { 01562 fSurfaceArea = G4VSolid::GetSurfaceArea(); 01563 } 01564 } 01565 return fSurfaceArea; 01566 }
Definition at line 76 of file G4GenericTrap.icc.
References FatalException, and G4Exception().
Referenced by CreatePolyhedron(), and SurfaceNormal().
00077 { 00078 if ( (index<0) || (index >= G4int(fVertices.size())) ) 00079 { 00080 G4Exception ("G4GenericTrap::GetTwistAngle()", "GeomSolids0003", 00081 FatalException, "Index outside range."); 00082 return 0.; 00083 } 00084 return fTwist[index]; 00085 }
G4TwoVector G4GenericTrap::GetVertex | ( | G4int | index | ) | const [inline] |
Definition at line 54 of file G4GenericTrap.icc.
References FatalException, and G4Exception().
00055 { 00056 if ( index<0 || index >= G4int(fVertices.size()) ) 00057 { 00058 G4Exception ("G4GenericTrap::GetVertex()", "GeomSolids0003", 00059 FatalException, "Index outside range."); 00060 return G4TwoVector(); 00061 } 00062 return fVertices[index]; 00063 }
const std::vector< G4TwoVector > & G4GenericTrap::GetVertices | ( | ) | const [inline] |
G4int G4GenericTrap::GetVisSubdivisions | ( | ) | const [inline] |
G4double G4GenericTrap::GetZHalfLength | ( | ) | const [inline] |
EInside G4GenericTrap::Inside | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 304 of file G4GenericTrap.cc.
References G4TessellatedSolid::Inside(), G4VSolid::kCarTolerance, kInside, kOutside, and kSurface.
Referenced by CalculateExtent(), and DistanceToIn().
00305 { 00306 // Test if point is inside this shape 00307 00308 #ifdef G4TESS_TEST 00309 if ( fTessellatedSolid ) 00310 { 00311 return fTessellatedSolid->Inside(p); 00312 } 00313 #endif 00314 00315 static const G4double halfCarTolerance=kCarTolerance*0.5; 00316 EInside innew=kOutside; 00317 std::vector<G4TwoVector> xy; 00318 00319 if (std::fabs(p.z()) <= fDz+halfCarTolerance) // First check Z range 00320 { 00321 // Compute intersection between Z plane containing point and the shape 00322 // 00323 G4double cf = 0.5*(fDz-p.z())/fDz; 00324 for (G4int i=0; i<4; i++) 00325 { 00326 xy.push_back(fVertices[i+4]+cf*( fVertices[i]-fVertices[i+4])); 00327 } 00328 00329 innew=InsidePolygone(p,xy); 00330 00331 if( (innew==kInside) || (innew==kSurface) ) 00332 { 00333 if(std::fabs(p.z()) > fDz-halfCarTolerance) { innew=kSurface; } 00334 } 00335 } 00336 return innew; 00337 }
G4bool G4GenericTrap::IsTwisted | ( | ) | const [inline] |
G4GenericTrap & G4GenericTrap::operator= | ( | const G4GenericTrap & | rhs | ) |
Definition at line 205 of file G4GenericTrap.cc.
References fCubicVolume, fDz, fIsTwisted, fMaxBBoxVector, fMinBBoxVector, fpPolyhedron, fSurfaceArea, fTessellatedSolid, fTwist, fVertices, fVisSubdivisions, and G4VSolid::operator=().
00206 { 00207 // Check assignment to self 00208 // 00209 if (this == &rhs) { return *this; } 00210 00211 // Copy base class data 00212 // 00213 G4VSolid::operator=(rhs); 00214 00215 // Copy data 00216 // 00217 fpPolyhedron = 0; fDz = rhs.fDz; fVertices = rhs.fVertices; 00218 fIsTwisted = rhs.fIsTwisted; fTessellatedSolid = 0; 00219 fMinBBoxVector = rhs.fMinBBoxVector; fMaxBBoxVector = rhs.fMaxBBoxVector; 00220 fVisSubdivisions = rhs.fVisSubdivisions; 00221 fSurfaceArea = rhs.fSurfaceArea; fCubicVolume = rhs.fCubicVolume; 00222 00223 for (size_t i=0; i<4; ++i) { fTwist[i] = rhs.fTwist[i]; } 00224 #ifdef G4TESS_TEST 00225 if (rhs.fTessellatedSolid && !fIsTwisted ) 00226 { delete fTessellatedSolid; fTessellatedSolid = CreateTessellatedSolid(); } 00227 #endif 00228 00229 return *this; 00230 }
void G4GenericTrap::SetVisSubdivisions | ( | G4int | subdiv | ) | [inline] |
std::ostream & G4GenericTrap::StreamInfo | ( | std::ostream & | os | ) | const [virtual] |
Implements G4VSolid.
Definition at line 1403 of file G4GenericTrap.cc.
References G4endl, GetEntityType(), and G4VSolid::GetName().
01404 { 01405 G4int oldprc = os.precision(16); 01406 os << "-----------------------------------------------------------\n" 01407 << " *** Dump for solid - " << GetName() << " *** \n" 01408 << " =================================================== \n" 01409 << " Solid geometry type: " << GetEntityType() << G4endl 01410 << " half length Z: " << fDz/mm << " mm \n" 01411 << " list of vertices:\n"; 01412 01413 for ( G4int i=0; i<fgkNofVertices; ++i ) 01414 { 01415 os << std::setw(5) << "#" << i 01416 << " vx = " << fVertices[i].x()/mm << " mm" 01417 << " vy = " << fVertices[i].y()/mm << " mm" << G4endl; 01418 } 01419 os.precision(oldprc); 01420 01421 return os; 01422 }
G4ThreeVector G4GenericTrap::SurfaceNormal | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 341 of file G4GenericTrap.cc.
References G4Exception(), GetTwistAngle(), JustWarning, G4VSolid::kCarTolerance, and G4TessellatedSolid::SurfaceNormal().
00342 { 00343 // Calculate side nearest to p, and return normal 00344 // If two sides are equidistant, sum of the Normal is returned 00345 00346 #ifdef G4TESS_TEST 00347 if ( fTessellatedSolid ) 00348 { 00349 return fTessellatedSolid->SurfaceNormal(p); 00350 } 00351 #endif 00352 00353 static const G4double halfCarTolerance=kCarTolerance*0.5; 00354 00355 G4ThreeVector lnorm, sumnorm(0.,0.,0.), apprnorm(0.,0.,1.), 00356 p0, p1, p2, r1, r2, r3, r4; 00357 G4int noSurfaces = 0; 00358 G4double distxy,distz; 00359 G4bool zPlusSide=false; 00360 00361 distz = fDz-std::fabs(p.z()); 00362 if (distz < halfCarTolerance) 00363 { 00364 if(p.z()>0) 00365 { 00366 zPlusSide=true; 00367 sumnorm=G4ThreeVector(0,0,1); 00368 } 00369 else 00370 { 00371 sumnorm=G4ThreeVector(0,0,-1); 00372 } 00373 noSurfaces ++; 00374 } 00375 00376 // Check lateral planes 00377 // 00378 std:: vector<G4TwoVector> vertices; 00379 G4double cf = 0.5*(fDz-p.z())/fDz; 00380 for (G4int i=0; i<4; i++) 00381 { 00382 vertices.push_back(fVertices[i+4]+cf*(fVertices[i]-fVertices[i+4])); 00383 } 00384 00385 // Compute distance for lateral planes 00386 // 00387 for (G4int q=0; q<4; q++) 00388 { 00389 p0=G4ThreeVector(vertices[q].x(),vertices[q].y(),p.z()); 00390 if(zPlusSide) 00391 { 00392 p1=G4ThreeVector(fVertices[q].x(),fVertices[q].y(),-fDz); 00393 } 00394 else 00395 { 00396 p1=G4ThreeVector(fVertices[q+4].x(),fVertices[q+4].y(),fDz); 00397 } 00398 p2=G4ThreeVector(vertices[(q+1)%4].x(),vertices[(q+1)%4].y(),p.z()); 00399 00400 // Collapsed vertices 00401 // 00402 if ( (p2-p0).mag2() < kCarTolerance ) 00403 { 00404 if ( std::fabs(p.z()+fDz) > kCarTolerance ) 00405 { 00406 p2=G4ThreeVector(fVertices[(q+1)%4].x(),fVertices[(q+1)%4].y(),-fDz); 00407 } 00408 else 00409 { 00410 p2=G4ThreeVector(fVertices[(q+1)%4+4].x(),fVertices[(q+1)%4+4].y(),fDz); 00411 } 00412 } 00413 lnorm = (p1-p0).cross(p2-p0); 00414 lnorm = lnorm.unit(); 00415 if(zPlusSide) { lnorm=-lnorm; } 00416 00417 // Adjust Normal for Twisted Surface 00418 // 00419 if ( (fIsTwisted) && (GetTwistAngle(q)!=0) ) 00420 { 00421 G4double normP=(p2-p0).mag(); 00422 if(normP) 00423 { 00424 G4double proj=(p-p0).dot(p2-p0)/normP; 00425 if(proj<0) { proj=0; } 00426 if(proj>normP) { proj=normP; } 00427 G4int j=(q+1)%4; 00428 r1=G4ThreeVector(fVertices[q+4].x(),fVertices[q+4].y(),fDz); 00429 r2=G4ThreeVector(fVertices[j+4].x(),fVertices[j+4].y(),fDz); 00430 r3=G4ThreeVector(fVertices[q].x(),fVertices[q].y(),-fDz); 00431 r4=G4ThreeVector(fVertices[j].x(),fVertices[j].y(),-fDz); 00432 r1=r1+proj*(r2-r1)/normP; 00433 r3=r3+proj*(r4-r3)/normP; 00434 r2=r1-r3; 00435 r4=r2.cross(p2-p0); r4=r4.unit(); 00436 lnorm=r4; 00437 } 00438 } // End if fIsTwisted 00439 00440 distxy=std::fabs((p0-p).dot(lnorm)); 00441 if ( distxy<halfCarTolerance ) 00442 { 00443 noSurfaces ++; 00444 00445 // Negative sign for Normal is taken for Outside Normal 00446 // 00447 sumnorm=sumnorm+lnorm; 00448 } 00449 00450 // For ApproxSurfaceNormal 00451 // 00452 if (distxy<distz) 00453 { 00454 distz=distxy; 00455 apprnorm=lnorm; 00456 } 00457 } // End for loop 00458 00459 // Calculate final Normal, add Normal in the Corners and Touching Sides 00460 // 00461 if ( noSurfaces == 0 ) 00462 { 00463 G4Exception("G4GenericTrap::SurfaceNormal(p)", "GeomSolids1002", 00464 JustWarning, "Point p is not on surface !?" ); 00465 sumnorm=apprnorm; 00466 // Add Approximative Surface Normal Calculation? 00467 } 00468 else if ( noSurfaces == 1 ) { sumnorm = sumnorm; } 00469 else { sumnorm = sumnorm.unit(); } 00470 00471 return sumnorm ; 00472 }
G4Polyhedron* G4GenericTrap::fpPolyhedron [mutable, protected] |