Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/solids/CSG/src/G4Orb.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 // ********************************************************************
  2 // * License and Disclaimer                                           *
  3 // *                                                                  *
  4 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  5 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  6 // * conditions of the Geant4 Software License,  included in the file *
  7 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  8 // * include a list of copyright holders.                             *
  9 // *                                                                  *
 10 // * Neither the authors of this software system, nor their employing *
 11 // * institutes,nor the agencies providing financial support for this *
 12 // * work  make  any representation or  warranty, express or implied, *
 13 // * regarding  this  software system or assume any liability for its *
 14 // * use.  Please see the license in the file  LICENSE  and URL above *
 15 // * for the full disclaimer and the limitation of liability.         *
 16 // *                                                                  *
 17 // * This  code  implementation is the result of  the  scientific and *
 18 // * technical work of the GEANT4 collaboration.                      *
 19 // * By using,  copying,  modifying or  distributing the software (or *
 20 // * any work based  on the software)  you  agree  to acknowledge its *
 21 // * use  in  resulting  scientific  publications,  and indicate your *
 22 // * acceptance of all terms of the Geant4 Software license.          *
 23 // ********************************************************************
 24 //
 25 // Implementation for G4Orb class
 26 //
 27 // 20.08.03 V.Grichine - created
 28 // 08.08.17 E.Tcherniaev - complete revision, speed-up
 29 // --------------------------------------------------------------------
 30 
 31 #include "G4Orb.hh"
 32 
 33 #if !defined(G4GEOM_USE_UORB)
 34 
 35 #include "G4TwoVector.hh"
 36 #include "G4VoxelLimits.hh"
 37 #include "G4AffineTransform.hh"
 38 #include "G4GeometryTolerance.hh"
 39 #include "G4BoundingEnvelope.hh"
 40 
 41 #include "G4VPVParameterisation.hh"
 42 
 43 #include "G4RandomDirection.hh"
 44 #include "Randomize.hh"
 45 
 46 #include "G4VGraphicsScene.hh"
 47 #include "G4VisExtent.hh"
 48 
 49 using namespace CLHEP;
 50 
 51 //////////////////////////////////////////////////////////////////////////
 52 //
 53 // Constructor
 54 
 55 G4Orb::G4Orb( const G4String& pName, G4double pRmax )
 56   : G4CSGSolid(pName), fRmax(pRmax)
 57 {
 58   Initialize();
 59 }
 60 
 61 //////////////////////////////////////////////////////////////////////////
 62 //
 63 // Fake default constructor - sets only member data and allocates memory
 64 //                            for usage restricted to object persistency
 65 
 66 G4Orb::G4Orb( __void__& a )
 67   : G4CSGSolid(a)
 68 {
 69 }
 70 
 71 //////////////////////////////////////////////////////////////////////////
 72 //
 73 // Destructor
 74 
 75 G4Orb::~G4Orb() = default;
 76 
 77 //////////////////////////////////////////////////////////////////////////
 78 //
 79 // Copy constructor
 80 
 81 G4Orb::G4Orb(const G4Orb&) = default;
 82 
 83 //////////////////////////////////////////////////////////////////////////
 84 //
 85 // Assignment operator
 86 
 87 G4Orb& G4Orb::operator = (const G4Orb& rhs)
 88 {
 89    // Check assignment to self
 90    //
 91    if (this == &rhs)  { return *this; }
 92 
 93    // Copy base class data
 94    //
 95    G4CSGSolid::operator=(rhs);
 96 
 97    // Copy data
 98    //
 99    fRmax = rhs.fRmax;
100    halfRmaxTol = rhs.halfRmaxTol;
101    sqrRmaxPlusTol = rhs.sqrRmaxPlusTol;
102    sqrRmaxMinusTol = rhs.sqrRmaxMinusTol;
103 
104    return *this;
105 }
106 
107 //////////////////////////////////////////////////////////////////////////
108 //
109 // Check radius and initialize dada members
110 
111 void G4Orb::Initialize()
112 {
113   const G4double fEpsilon = 2.e-11;  // relative tolerance of fRmax
114 
115   // Check radius
116   //
117   if ( fRmax < 10*kCarTolerance )
118   {
119     G4Exception("G4Orb::Initialize()", "GeomSolids0002", FatalException,
120                 "Invalid radius < 10*kCarTolerance.");
121   }
122   halfRmaxTol = 0.5 * std::max(kCarTolerance, fEpsilon*fRmax);
123   G4double rmaxPlusTol  = fRmax + halfRmaxTol;
124   G4double rmaxMinusTol = fRmax - halfRmaxTol;
125   sqrRmaxPlusTol = rmaxPlusTol*rmaxPlusTol;
126   sqrRmaxMinusTol = rmaxMinusTol*rmaxMinusTol;
127 }
128 
129 //////////////////////////////////////////////////////////////////////////
130 //
131 // Dispatch to parameterisation for replication mechanism dimension
132 // computation & modification
133 
134 void G4Orb::ComputeDimensions(       G4VPVParameterisation* p,
135                                const G4int n,
136                                const G4VPhysicalVolume* pRep )
137 {
138   p->ComputeDimensions(*this,n,pRep);
139 }
140 
141 //////////////////////////////////////////////////////////////////////////
142 //
143 // Get bounding box
144 
145 void G4Orb::BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const
146 {
147   G4double radius = GetRadius();
148   pMin.set(-radius,-radius,-radius);
149   pMax.set( radius, radius, radius);
150 
151   // Check correctness of the bounding box
152   //
153   if (pMin.x() >= pMax.x() || pMin.y() >= pMax.y() || pMin.z() >= pMax.z())
154   {
155     std::ostringstream message;
156     message << "Bad bounding box (min >= max) for solid: "
157             << GetName() << " !"
158             << "\npMin = " << pMin
159             << "\npMax = " << pMax;
160     G4Exception("G4Orb::BoundingLimits()", "GeomMgt0001",
161                 JustWarning, message);
162     DumpInfo();
163   }
164 }
165 
166 //////////////////////////////////////////////////////////////////////////
167 //
168 // Calculate extent under transform and specified limit
169 
170 G4bool G4Orb::CalculateExtent(const EAxis pAxis,
171                               const G4VoxelLimits& pVoxelLimit,
172                               const G4AffineTransform& pTransform,
173                                         G4double& pMin, G4double& pMax) const
174 {
175   G4ThreeVector bmin, bmax;
176   G4bool exist;
177 
178   // Get bounding box
179   BoundingLimits(bmin,bmax);
180 
181   // Check bounding box
182   G4BoundingEnvelope bbox(bmin,bmax);
183 #ifdef G4BBOX_EXTENT
184   return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
185 #endif
186   if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
187   {
188     return exist = pMin < pMax;
189   }
190 
191   // Find bounding envelope and calculate extent
192   //
193   static const G4int NTHETA = 8;  // number of steps along Theta
194   static const G4int NPHI   = 16; // number of steps along Phi
195   static const G4double sinHalfTheta = std::sin(halfpi/NTHETA);
196   static const G4double cosHalfTheta = std::cos(halfpi/NTHETA);
197   static const G4double sinHalfPhi   = std::sin(pi/NPHI);
198   static const G4double cosHalfPhi   = std::cos(pi/NPHI);
199   static const G4double sinStepTheta = 2.*sinHalfTheta*cosHalfTheta;
200   static const G4double cosStepTheta = 1. - 2.*sinHalfTheta*sinHalfTheta;
201   static const G4double sinStepPhi   = 2.*sinHalfPhi*cosHalfPhi;
202   static const G4double cosStepPhi   = 1. - 2.*sinHalfPhi*sinHalfPhi;
203 
204   G4double radius = GetRadius();
205   G4double rtheta = radius/cosHalfTheta;
206   G4double rphi   = rtheta/cosHalfPhi;
207 
208   // set reference circle
209   G4TwoVector xy[NPHI];
210   G4double sinCurPhi = sinHalfPhi;
211   G4double cosCurPhi = cosHalfPhi;
212   for (auto & k : xy)
213   {
214     k.set(cosCurPhi,sinCurPhi);
215     G4double sinTmpPhi = sinCurPhi;
216     sinCurPhi = sinCurPhi*cosStepPhi + cosCurPhi*sinStepPhi;
217     cosCurPhi = cosCurPhi*cosStepPhi - sinTmpPhi*sinStepPhi;
218   }
219   
220   // set bounding circles
221   G4ThreeVectorList circles[NTHETA];
222   for (auto & circle : circles) { circle.resize(NPHI); }
223 
224   G4double sinCurTheta = sinHalfTheta;
225   G4double cosCurTheta = cosHalfTheta;
226   for (auto & circle : circles)
227   {
228     G4double z = rtheta*cosCurTheta;
229     G4double rho = rphi*sinCurTheta;
230     for (G4int k=0; k<NPHI; ++k)
231     {
232       circle[k].set(rho*xy[k].x(),rho*xy[k].y(),z);
233     }
234     G4double sinTmpTheta = sinCurTheta;
235     sinCurTheta = sinCurTheta*cosStepTheta + cosCurTheta*sinStepTheta;
236     cosCurTheta = cosCurTheta*cosStepTheta - sinTmpTheta*sinStepTheta;
237   }
238 
239   // set envelope and calculate extent
240   std::vector<const G4ThreeVectorList *> polygons;
241   polygons.resize(NTHETA);
242   for (G4int i=0; i<NTHETA; ++i) { polygons[i] = &circles[i]; }
243 
244   G4BoundingEnvelope benv(bmin,bmax,polygons);
245   exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
246   return exist;
247 }
248 
249 //////////////////////////////////////////////////////////////////////////
250 //
251 // Return whether point is inside/outside/on surface
252 
253 EInside G4Orb::Inside( const G4ThreeVector& p ) const
254 {
255   G4double rr = p.mag2();
256   if (rr > sqrRmaxPlusTol) return kOutside;
257   return (rr > sqrRmaxMinusTol) ? kSurface : kInside;
258 }
259 
260 //////////////////////////////////////////////////////////////////////////
261 //
262 // Return unit normal of surface closest to p
263 
264 G4ThreeVector G4Orb::SurfaceNormal( const G4ThreeVector& p ) const
265 {
266   return (1/p.mag())*p;
267 }
268 
269 //////////////////////////////////////////////////////////////////////////
270 //
271 // Calculate distance to the surface of the orb from outside
272 // - return kInfinity if no intersection or
273 //   intersection distance <= tolerance
274 
275 G4double G4Orb::DistanceToIn( const G4ThreeVector& p,
276                               const G4ThreeVector& v  ) const
277 {
278   // Check if point is on the surface and traveling away
279   //
280   G4double rr = p.mag2();
281   G4double pv = p.dot(v);
282   if (rr >= sqrRmaxMinusTol && pv >= 0) return kInfinity;
283 
284   // Find intersection
285   //
286   //    Sphere eqn: x^2 + y^2 + z^2 = R^2
287   //
288   //    => (px + t*vx)^2 + (py + t*vy)^2 + (pz + t*vz)^2 = R^2
289   //    => r^2 + 2t(p.v) + t^2 = R^2
290   //    => tmin = -(p.v) - Sqrt((p.v)^2 - (r^2 - R^2))
291   //
292   G4double D  = pv*pv - rr + fRmax*fRmax;
293   if (D < 0) return kInfinity; // no intersection
294 
295   G4double sqrtD = std::sqrt(D);
296   G4double dist = -pv - sqrtD;
297 
298   // Avoid rounding errors due to precision issues seen on 64 bits systems.
299   // Split long distances and recompute
300   //
301   G4double Dmax = 32*fRmax; 
302   if (dist > Dmax)
303   {
304     dist  = dist - 1.e-8*dist - fRmax; // to stay outside after the move
305     dist += DistanceToIn(p + dist*v, v);
306     return (dist >= kInfinity) ? kInfinity : dist;
307   }
308 
309   if (sqrtD*2 <= halfRmaxTol) return kInfinity; // touch
310   return (dist < halfRmaxTol) ? 0. : dist;
311 }
312 
313 //////////////////////////////////////////////////////////////////////////
314 //
315 // Calculate shortest distance to the boundary from outside
316 // - Return 0 if point is inside
317 
318 G4double G4Orb::DistanceToIn( const G4ThreeVector& p ) const
319 {
320   G4double dist = p.mag() - fRmax;
321   return (dist > 0) ? dist : 0.;
322 }
323 
324 //////////////////////////////////////////////////////////////////////////
325 //
326 // Calculate distance to the surface of the orb from inside and
327 // find normal at exit point, if required
328 // - when leaving the surface, return 0
329 
330 G4double G4Orb::DistanceToOut( const G4ThreeVector& p,
331                                const G4ThreeVector& v,
332                                const G4bool calcNorm,
333                                      G4bool* validNorm,
334                                      G4ThreeVector* n ) const
335 {
336   // Check if point is on the surface and traveling away
337   //
338   G4double rr = p.mag2();
339   G4double pv = p.dot(v);
340   if (rr >= sqrRmaxMinusTol && pv > 0)
341   {
342     if (calcNorm)
343     {
344       *validNorm = true;
345       *n = p*(1./std::sqrt(rr));
346     }
347     return 0.;
348   }
349 
350   // Find intersection
351   //
352   //    Sphere eqn: x^2 + y^2 + z^2 = R^2
353   //
354   //    => (px + t*vx)^2 + (py + t*vy)^2 + (pz + t*vz)^2 = R^2
355   //    => r^2 + 2t(p.v) + t^2 = R^2
356   //    => tmax = -(p.v) + Sqrt((p.v)^2 - (r^2 - R^2))
357   //
358   G4double D  = pv*pv - rr + fRmax*fRmax;
359   G4double tmax = (D <= 0) ? 0. : std::sqrt(D) - pv;
360   if (tmax < halfRmaxTol) tmax = 0.;
361   if (calcNorm)
362   {
363     *validNorm = true;
364     G4ThreeVector ptmax = p + tmax*v;
365     *n = ptmax*(1./ptmax.mag());
366   }
367   return tmax;
368 }
369 
370 //////////////////////////////////////////////////////////////////////////
371 //
372 // Calculate distance (<=actual) to closest surface of shape from inside
373 
374 G4double G4Orb::DistanceToOut( const G4ThreeVector& p ) const
375 {
376 #ifdef G4CSGDEBUG
377   if( Inside(p) == kOutside )
378   {
379     std::ostringstream message;
380     G4int oldprc = message.precision(16);
381     message << "Point p is outside (!?) of solid: " << GetName() << "\n";
382     message << "Position:\n";
383     message << "   p.x() = " << p.x()/mm << " mm\n";
384     message << "   p.y() = " << p.y()/mm << " mm\n";
385     message << "   p.z() = " << p.z()/mm << " mm";
386     G4cout.precision(oldprc);
387     G4Exception("G4Trap::DistanceToOut(p)", "GeomSolids1002",
388                 JustWarning, message );
389     DumpInfo();
390   }
391 #endif
392   G4double dist = fRmax - p.mag();
393   return (dist > 0) ? dist : 0.;
394 }
395 
396 //////////////////////////////////////////////////////////////////////////
397 //
398 // G4EntityType
399 
400 G4GeometryType G4Orb::GetEntityType() const
401 {
402   return {"G4Orb"};
403 }
404 
405 //////////////////////////////////////////////////////////////////////////
406 //
407 // Make a clone of the object
408 
409 G4VSolid* G4Orb::Clone() const
410 {
411   return new G4Orb(*this);
412 }
413 
414 //////////////////////////////////////////////////////////////////////////
415 //
416 // Stream object contents to an output stream
417 
418 std::ostream& G4Orb::StreamInfo( std::ostream& os ) const
419 {
420   G4long oldprc = os.precision(16);
421   os << "-----------------------------------------------------------\n"
422      << "    *** Dump for solid - " << GetName() << " ***\n"
423      << "    ===================================================\n"
424      << " Solid type: G4Orb\n"
425      << " Parameters: \n"
426      << "    outer radius: " << fRmax/mm << " mm \n"
427      << "-----------------------------------------------------------\n";
428   os.precision(oldprc);
429   return os;
430 }
431 
432 //////////////////////////////////////////////////////////////////////////
433 //
434 // GetPointOnSurface
435 
436 G4ThreeVector G4Orb::GetPointOnSurface() const
437 {
438   return fRmax * G4RandomDirection();
439 }
440 
441 //////////////////////////////////////////////////////////////////////////
442 //
443 // Methods for visualisation
444 
445 void G4Orb::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
446 {
447   scene.AddSolid (*this);
448 }
449 
450 G4VisExtent G4Orb::GetExtent() const
451 {
452   return {-fRmax, fRmax, -fRmax, fRmax, -fRmax, fRmax};
453 }
454 
455 G4Polyhedron* G4Orb::CreatePolyhedron () const
456 {
457   return new G4PolyhedronSphere (0., fRmax, 0., 2*pi, 0., pi);
458 }
459 
460 #endif
461