00001 // 00002 // ******************************************************************** 00003 // * License and Disclaimer * 00004 // * * 00005 // * The Geant4 software is copyright of the Copyright Holders of * 00006 // * the Geant4 Collaboration. It is provided under the terms and * 00007 // * conditions of the Geant4 Software License, included in the file * 00008 // * LICENSE and available at http://cern.ch/geant4/license . These * 00009 // * include a list of copyright holders. * 00010 // * * 00011 // * Neither the authors of this software system, nor their employing * 00012 // * institutes,nor the agencies providing financial support for this * 00013 // * work make any representation or warranty, express or implied, * 00014 // * regarding this software system or assume any liability for its * 00015 // * use. Please see the license in the file LICENSE and URL above * 00016 // * for the full disclaimer and the limitation of liability. * 00017 // * * 00018 // * This code implementation is the result of the scientific and * 00019 // * technical work of the GEANT4 collaboration. * 00020 // * By using, copying, modifying or distributing the software (or * 00021 // * any work based on the software) you agree to acknowledge its * 00022 // * use in resulting scientific publications, and indicate your * 00023 // * acceptance of all terms of the Geant4 Software license. * 00024 // ******************************************************************** 00025 // 00026 // 00027 // $Id$ 00028 // 00029 // ---------------------------------------------------------------------- 00030 // GEANT 4 class source file 00031 // 00032 // G4Surface.cc 00033 // 00034 // ---------------------------------------------------------------------- 00035 00036 #include "G4Surface.hh" 00037 #include "G4CompositeCurve.hh" 00038 #include "G4GeometryTolerance.hh" 00039 00040 G4Surface::G4Surface() 00041 : G4STEPEntity(), 00042 bbox(0), next(0), Intersected(0), Type(0), AdvancedFace(0), 00043 active(1), distance(kInfinity), uhit(0.), vhit(0.), sameSense(0) 00044 { 00045 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 00046 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance(); 00047 } 00048 00049 00050 G4Surface::~G4Surface() 00051 { 00052 } 00053 00054 00055 G4Surface::G4Surface( const G4Surface& c ) 00056 : G4STEPEntity(), 00057 bbox(c.bbox), next(c.next), Intersected(c.Intersected), Type(c.Type), 00058 AdvancedFace(c.AdvancedFace), active(c.active), distance(c.distance), 00059 uhit(c.uhit), vhit(c.vhit), sameSense(c.sameSense) 00060 { 00061 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 00062 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance(); 00063 } 00064 00065 00066 G4Surface& 00067 G4Surface::operator=( const G4Surface& c ) 00068 { 00069 if (&c == this) { return *this; } 00070 bbox = c.bbox; 00071 next = c.next; 00072 Intersected = c.Intersected; 00073 Type = c.Type; 00074 AdvancedFace = c.AdvancedFace; 00075 active = c.active; 00076 distance = c.distance; 00077 uhit = c.uhit; 00078 vhit = c.vhit; 00079 00080 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 00081 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance(); 00082 00083 return *this; 00084 } 00085 00086 00087 G4int G4Surface::operator==( const G4Surface& surf ) 00088 { 00089 return origin == surf.origin; 00090 } 00091 00092 G4String G4Surface::GetEntityType() const 00093 { 00094 return G4String("Surface"); 00095 } 00096 00097 const char* G4Surface::Name() const 00098 { 00099 return "G4Surface"; 00100 } 00101 00102 G4int G4Surface::MyType() const 00103 { 00104 return Type; 00105 } 00106 00107 void G4Surface::InitBounded() 00108 { 00109 } 00110 00111 G4double G4Surface::GetUHit() const 00112 { 00113 return uhit; 00114 } 00115 00116 G4double G4Surface::GetVHit() const 00117 { 00118 return vhit; 00119 } 00120 00121 //void G4Surface::read_surface(fstream& tmp){;} 00122 00123 G4Point3D G4Surface::Evaluation(const G4Ray&) 00124 { 00125 return closest_hit; 00126 } 00127 00128 G4int G4Surface::Evaluate(const G4Ray&) 00129 { 00130 return 0; 00131 } 00132 00133 void G4Surface::Reset() 00134 { 00135 Intersected = 0; 00136 active = 1; 00137 distance = kInfinity; 00138 } 00139 00140 void G4Surface::SetBoundaries(G4CurveVector* boundaries) 00141 { 00142 surfaceBoundary.Init(*boundaries); 00143 InitBounded(); 00144 } 00145 00146 void G4Surface::CalcBBox() 00147 { 00148 // Finds the bounds of the surface iow 00149 // calculates the bounds for a bounding box 00150 // to the surface. The bounding box is used 00151 // for a preliminary check of intersection. 00152 00153 bbox = new G4BoundingBox3D(surfaceBoundary.BBox().GetBoxMin(), 00154 surfaceBoundary.BBox().GetBoxMax()); 00155 // old implementation 00156 // G4Point3d BoundaryMax = OuterBoundary->GetBoundsMax(); 00157 // G4Point3d BoundaryMin = OuterBoundary->GetBoundsMin(); 00158 // bbox = new G4BoundingBox( BoundaryMin, BoundaryMax); 00159 // return; 00160 } 00161 00162 G4Vector3D G4Surface::Normal( const G4Vector3D& ) const 00163 { // return the Normal unit vector to a Surface at the point p on 00164 // (or nearly on) the Surface. 00165 // The default is not well defined, so return ( 0, 0, 0 ). 00166 return G4Vector3D( 0.0, 0.0, 0.0 ); 00167 } 00168 00169 00170 G4int G4Surface::Intersect(const G4Ray&) 00171 { 00172 G4int Result = 0; 00173 00174 G4Exception("G4Surface::Intersect()", "GeomSolids0001", 00175 FatalException, "Sorry, not yet implemented."); 00176 00177 #ifdef NEW_IMPLEMENTATION 00178 // get the intersection 00179 // Result = Intersect(rayref); 00180 00181 // Check that the point is within the polyline 00182 // Get Normal at Hitpoint 00183 const G4Vector3D& Vec = Normal(closest_hit); 00184 G4Ray Normal(closest_hit, Vec); 00185 00186 // Project points & Hit 00187 // OuterBoundary->ProjectBoundaryTo2D(Normal.GetPlane(1), 00188 // Normal.GetPlane(2), 0); 00189 00190 00191 00192 00193 G4Point3D Hit = closest_hit.Project(Normal.GetPlane(1), 00194 Normal.GetPlane(2) ); 00195 // Check point in polygon 00196 // Result = OuterBoundary->Inside(Hit, rayref); 00197 00198 #endif 00199 return Result; 00200 } 00201 00202 00203 G4double G4Surface::ClosestDistanceToPoint(const G4Point3D& Pt) 00204 { 00205 // in fact, a squared distance is returned 00206 00207 // a bit suspicious, this function 00208 // the distance is almost always an overestimate 00209 G4double pointDistance= kInfinity; 00210 G4double tmpDistance; 00211 const G4CurveVector& bounds= surfaceBoundary.GetBounds(); 00212 00213 G4int entr = bounds.size(); 00214 00215 for (G4int i=0; i<entr; i++) 00216 { 00217 G4Curve* c= bounds[i]; 00218 00219 if (c->GetEntityType() == "G4CompositeCurve") 00220 { 00221 G4CompositeCurve* cc= (G4CompositeCurve*)c; 00222 const G4CurveVector& segments= cc->GetSegments(); 00223 for (size_t j=0; j<segments.size(); j++) 00224 { 00225 G4Curve* ccc= segments[j]; 00226 tmpDistance= (G4Point3D(Pt.x(), Pt.y(), Pt.z())-ccc->GetEnd()).mag2(); 00227 if (pointDistance > tmpDistance) 00228 { 00229 pointDistance= tmpDistance; 00230 } 00231 } 00232 00233 } 00234 else 00235 { 00236 tmpDistance= (G4Point3D(Pt.x(), Pt.y(), Pt.z())-c->GetEnd()).mag2(); 00237 if (pointDistance > tmpDistance) 00238 { 00239 pointDistance= tmpDistance; 00240 } 00241 } 00242 } 00243 00244 // L. Broglia 00245 // Be carreful ! pointdistance is the squared distance 00246 return std::sqrt(pointDistance); 00247 00248 // G4double PointDistance=kInfinity; 00249 // G4double TmpDistance=0; 00250 // PointDistance = OuterBoundary->ClosestDistanceToPoint(Pt); 00251 // TmpDistance =0; 00252 // for(G4int a=0;a<NumberOfInnerBoundaries;a++) 00253 // { 00254 // TmpDistance = InnerBoundary[a]->ClosestDistanceToPoint(Pt); 00255 // if(PointDistance > TmpDistance) PointDistance = TmpDistance; 00256 // } 00257 // return PointDistance; 00258 00259 //G4double G4Boundary::ClosestDistanceToPoint(const G4ThreeVec& Pt) 00260 //{ 00261 // G4double PointDistance = kInfinity; 00262 // G4double TmpDistance = 0; 00263 // for(G4int a =0; a < NumberOfPoints;a++) 00264 // { 00265 // G4Point3d& Pt2 = GetPoint(a); 00266 // TmpDistance = Pt2.Distance(Pt); 00267 // if(PointDistance > TmpDistance)PointDistance = TmpDistance; 00268 // } 00269 // return PointDistance; 00270 //} 00271 } 00272 00273 00274 std::ostream& operator<<( std::ostream& os, const G4Surface& ) 00275 { 00276 // overwrite output operator << to Print out Surface objects 00277 // using the PrintOn function defined below 00278 // s.PrintOn( os ); 00279 return os; 00280 } 00281 00282 00283 G4double G4Surface::HowNear( const G4Vector3D& x ) const 00284 { 00285 // Distance from the point x to a Surface. 00286 // The default for a Surface is the distance from the point to the origin. 00287 G4Vector3D p = G4Vector3D( x - origin ); 00288 return p.mag(); 00289 } 00290 00291 void G4Surface::Project() 00292 { 00293 } 00294 00295 void G4Surface::CalcNormal() 00296 { 00297 } 00298 00299 G4int G4Surface::IsConvex() const 00300 { 00301 return -1; 00302 } 00303 00304 G4int G4Surface::GetConvex() const 00305 { 00306 return 0; 00307 } 00308 00309 G4int G4Surface::GetNumberOfPoints() const 00310 { 00311 return 0; 00312 } 00313 00314 const G4Point3D& G4Surface::GetPoint(G4int) const 00315 { 00316 const G4Point3D* tmp= new G4Point3D(0,0,0); 00317 return *tmp; 00318 } 00319 00320 G4Ray* G4Surface::Norm() 00321 { 00322 return (G4Ray*)0; 00323 } 00324 00325 void G4Surface::Project (G4double& Coord, 00326 const G4Point3D& Pt2, 00327 const G4Plane& Pl1) 00328 { 00329 Coord = Pt2.x()*Pl1.a + Pt2.y()*Pl1.b + Pt2.z()*Pl1.c - Pl1.d; 00330 } 00331 00332 /* 00333 G4double G4Surface::distanceAlongRay( G4int which_way, const G4Ray* ry, 00334 G4ThreeVec& p ) const 00335 { // Distance along a Ray (straight line with G4ThreeVec) to leave or enter 00336 // a Surface. The input variable which_way should be set to +1 to indicate 00337 // leaving a Surface, -1 to indicate entering a Surface. 00338 // p is the point of intersection of the Ray with the Surface. 00339 // This is a default function which just gives the distance 00340 // between the origin of the Ray and the origin of the Surface. 00341 // Since a generic Surface doesn't have a well-defined Normal, no 00342 // further checks are Done. 00343 00344 // This should always be overwritten for derived classes so Print out 00345 // a warning message if this is called. 00346 G4cout << "WARNING from Surface::distanceAlongRay\n" 00347 << " This function should be overwritten by a derived class.\n" 00348 << " Using the Surface base class default.\n"; 00349 p = GetOrigin(); 00350 G4ThreeVec d = ry->Position() - p; 00351 return d.Magnitude(); 00352 } 00353 00354 G4double G4Surface::distanceAlongHelix( G4int which_way, const Helix* hx, 00355 G4ThreeVec& p ) const 00356 { // Distance along a Helix to leave or enter a Surface. 00357 // The input variable which_way should be set to +1 to indicate 00358 // leaving a Surface, -1 to indicate entering a Surface. 00359 // p is the point of intersection of the Helix with the Surface. 00360 // This is a default function which just gives the distance 00361 // between the origin of the Helix and the origin of the Surface. 00362 // Since a generic Surface doesn't have a well-defined Normal, no 00363 // further checks are Done. 00364 00365 // This should always be overwritten for derived classes so Print out 00366 // a warning message if this is called. 00367 G4cout << "WARNING from Surface::distanceAlongHelix\n" 00368 << " This function should be overwritten by a derived class.\n" 00369 << " Using the Surface base class default.\n"; 00370 p = GetOrigin(); 00371 G4ThreeVec d = hx->position() - p; 00372 return d.Magnitude(); 00373 } 00374 00375 00376 G4ThreeVec G4Surface::Normal() const 00377 { // return the Normal unit vector to a Surface 00378 // (This is only meaningful for Surfaces for which the Normal does 00379 // not depend on location on the Surface). 00380 // The default is not well defined, so return ( 0, 0, 0 ). 00381 return G4ThreeVec( 0.0, 0.0, 0.0 ); 00382 } 00383 00384 00385 G4ThreeVec G4Surface::Normal( const G4ThreeVec& ) const 00386 { // return the Normal unit vector to a Surface at the point p on 00387 // (or nearly on) the Surface. 00388 // The default is not well defined, so return ( 0, 0, 0 ). 00389 return G4ThreeVec( 0.0, 0.0, 0.0 ); 00390 } 00391 00392 G4int G4Surface::Inside( const G4ThreeVec& ) const 00393 { // return 0 if point p is outside Surface, 1 if Inside 00394 // default is not well defined, so return 0 00395 return 0; 00396 } 00397 00398 void G4Surface::move( const G4ThreeVec& p ) 00399 { // translate origin of Surface by vector p 00400 origin += p; 00401 } 00402 00403 void G4Surface::rotate( G4double alpha, G4double beta, 00404 G4double gamma, G4ThreeMat& m, G4int inverse ) 00405 { // rotate Surface first about global x-axis by angle alpha, 00406 // second about global y-axis by angle beta, 00407 // and third about global z-axis by angle gamma 00408 // by creating and using G4ThreeMat objects 00409 // angles are assumed to be given in radians 00410 // returns also the overall rotation matrix for use by subclasses 00411 // if inverse is non-zero, the order of rotations is reversed 00412 // for a generic Surface, only the origin is rotated 00413 // G4double ax[3][3] = { 0., 0., 0., 0., 0., 0., 0., 0., 0. }; 00414 G4double ax[3][3]; 00415 G4double ay[3][3]; 00416 G4double az[3][3]; 00417 // G4double ay[3][3] = { 0., 0., 0., 0., 0., 0., 0., 0., 0. }; 00418 // G4double az[3][3] = { 0., 0., 0., 0., 0., 0., 0., 0., 0. }; 00419 ax[0][0] = 1.; 00420 ax[1][1] = std::cos( alpha ); 00421 ax[2][2] = ax[1][1]; 00422 ax[2][1] = std::sin( alpha ); 00423 ax[1][2] = -ax[2][1]; 00424 ay[1][1] = 1.; 00425 ay[0][0] = std::cos( beta ); 00426 ay[2][2] = ay[0][0]; 00427 ay[0][2] = std::sin( beta ); 00428 ay[2][0] = -ay[0][2]; 00429 az[2][2] = 1.; 00430 az[0][0] = std::cos( gamma ); 00431 az[1][1] = az[0][0]; 00432 az[1][0] = std::sin( gamma ); 00433 az[0][1] = -az[1][0]; 00434 G4ThreeMat &Rx = *new G4ThreeMat( ax ); // x-rotation matrix 00435 G4ThreeMat &Ry = *new G4ThreeMat( ay ); // y-rotation matrix 00436 G4ThreeMat &Rz = *new G4ThreeMat( az ); // z-rotation matrix 00437 if ( inverse ) 00438 m = Rx * ( Ry * Rz ); 00439 else 00440 m = Rz * ( Ry * Rx ); 00441 origin = m * origin; 00442 } 00443 00444 void G4Surface::rotate( G4double alpha, G4double beta, 00445 G4double gamma, G4int inverse ) 00446 { // rotate Surface first about global x-axis by angle alpha, 00447 // second about global y-axis by angle beta, 00448 // and third about global z-axis by angle gamma 00449 // by creating and using G4ThreeMat objects 00450 // angles are assumed to be given in radians 00451 // if inverse is non-zero, the order of rotations is reversed 00452 G4ThreeMat m; 00453 // Just call the above function to do this rotation 00454 rotate( alpha, beta, gamma, m, inverse ); 00455 } 00456 00457 */