source: trunk/source/geometry/solids/CSG/src/G4Box.cc@ 1093

Last change on this file since 1093 was 1058, checked in by garnier, 17 years ago

file release beta
File size: 27.6 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4Box.cc,v 1.44 2006/10/19 15:33:37 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
29//
30//
31//
32// Implementation for G4Box class
33//
34// 24.06.98 - V. Grichine: insideEdge in DistanceToIn(p,v)
35// 20.09.98 - V.Grichine: new algorithm of DistanceToIn(p,v)
36// 07.05.00 - V.Grichine: d= DistanceToIn(p,v), if d<e/2, d=0
37// 09.06.00 - V.Grichine: safety in DistanceToIn(p) against Inside(p)=kOutside
38// and information before exception in DistanceToOut(p,v,...)
39// 15.11.00 - D.Williams, V.Grichine: bug fixed in CalculateExtent - change
40// algorithm for rotated vertices
41// --------------------------------------------------------------------
42
43#include "G4Box.hh"
44
45#include "G4VoxelLimits.hh"
46#include "G4AffineTransform.hh"
47#include "Randomize.hh"
48
49#include "G4VPVParameterisation.hh"
50
51#include "G4VGraphicsScene.hh"
52#include "G4Polyhedron.hh"
53#include "G4NURBS.hh"
54#include "G4NURBSbox.hh"
55#include "G4VisExtent.hh"
56
57////////////////////////////////////////////////////////////////////////
58//
59// Constructor - check & set half widths
60
61G4Box::G4Box(const G4String& pName,
62 G4double pX,
63 G4double pY,
64 G4double pZ)
65 : G4CSGSolid(pName)
66{
67 if ( (pX > 2*kCarTolerance)
68 && (pY > 2*kCarTolerance)
69 && (pZ > 2*kCarTolerance) )
70 {
71 fDx = pX ;
72 fDy = pY ;
73 fDz = pZ ;
74 }
75 else
76 {
77 G4cerr << "ERROR - G4Box()::G4Box(): " << GetName() << G4endl
78 << " Dimensions too small ! - "
79 << pX << ", " << pY << ", " << pZ << G4endl;
80 G4Exception("G4Box::G4Box()", "InvalidSetup",
81 FatalException, "Invalid dimensions. Too small.");
82 }
83}
84
85//////////////////////////////////////////////////////////////////////////
86//
87// Fake default constructor - sets only member data and allocates memory
88// for usage restricted to object persistency.
89
90G4Box::G4Box( __void__& a )
91 : G4CSGSolid(a)
92{
93}
94
95//////////////////////////////////////////////////////////////////////////
96//
97// Destructor
98
99G4Box::~G4Box()
100{
101}
102
103//////////////////////////////////////////////////////////////////////////////
104
105void G4Box::SetXHalfLength(G4double dx)
106{
107 if(dx > 2*kCarTolerance)
108 fDx = dx;
109 else
110 {
111 G4cerr << "ERROR - G4Box()::SetXHalfLength(): " << GetName() << G4endl
112 << " Dimension X too small ! - "
113 << dx << G4endl;
114 G4Exception("G4Box::SetXHalfLength()", "InvalidSetup",
115 FatalException, "Invalid dimensions. Too small.");
116 }
117 fCubicVolume= 0.;
118 fSurfaceArea= 0.;
119 fpPolyhedron = 0;
120}
121
122void G4Box::SetYHalfLength(G4double dy)
123{
124 if(dy > 2*kCarTolerance)
125 fDy = dy;
126 else
127 {
128 G4cerr << "ERROR - G4Box()::SetYHalfLength(): " << GetName() << G4endl
129 << " Dimension Y too small ! - "
130 << dy << G4endl;
131 G4Exception("G4Box::SetYHalfLength()", "InvalidSetup",
132 FatalException, "Invalid dimensions. Too small.");
133 }
134 fCubicVolume= 0.;
135 fSurfaceArea= 0.;
136 fpPolyhedron = 0;
137}
138
139void G4Box::SetZHalfLength(G4double dz)
140{
141 if(dz > 2*kCarTolerance)
142 fDz = dz;
143 else
144 {
145 G4cerr << "ERROR - G4Box()::SetZHalfLength(): " << GetName() << G4endl
146 << " Dimension Z too small ! - "
147 << dz << G4endl;
148 G4Exception("G4Box::SetZHalfLength()", "InvalidSetup",
149 FatalException, "Invalid dimensions. Too small.");
150 }
151 fCubicVolume= 0.;
152 fSurfaceArea= 0.;
153 fpPolyhedron = 0;
154}
155
156
157
158////////////////////////////////////////////////////////////////////////
159//
160// Dispatch to parameterisation for replication mechanism dimension
161// computation & modification.
162
163void G4Box::ComputeDimensions(G4VPVParameterisation* p,
164 const G4int n,
165 const G4VPhysicalVolume* pRep)
166{
167 p->ComputeDimensions(*this,n,pRep);
168}
169
170//////////////////////////////////////////////////////////////////////////
171//
172// Calculate extent under transform and specified limit
173
174G4bool G4Box::CalculateExtent(const EAxis pAxis,
175 const G4VoxelLimits& pVoxelLimit,
176 const G4AffineTransform& pTransform,
177 G4double& pMin, G4double& pMax) const
178{
179 if (!pTransform.IsRotated())
180 {
181 // Special case handling for unrotated boxes
182 // Compute x/y/z mins and maxs respecting limits, with early returns
183 // if outside limits. Then switch() on pAxis
184
185 G4double xoffset,xMin,xMax;
186 G4double yoffset,yMin,yMax;
187 G4double zoffset,zMin,zMax;
188
189 xoffset = pTransform.NetTranslation().x() ;
190 xMin = xoffset - fDx ;
191 xMax = xoffset + fDx ;
192
193 if (pVoxelLimit.IsXLimited())
194 {
195 if ( xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance ||
196 xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance ) return false ;
197 else
198 {
199 if (xMin < pVoxelLimit.GetMinXExtent())
200 {
201 xMin = pVoxelLimit.GetMinXExtent() ;
202 }
203 if (xMax > pVoxelLimit.GetMaxXExtent())
204 {
205 xMax = pVoxelLimit.GetMaxXExtent() ;
206 }
207 }
208 }
209 yoffset = pTransform.NetTranslation().y() ;
210 yMin = yoffset - fDy ;
211 yMax = yoffset + fDy ;
212
213 if (pVoxelLimit.IsYLimited())
214 {
215 if ( yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance ||
216 yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance ) return false ;
217 else
218 {
219 if (yMin < pVoxelLimit.GetMinYExtent())
220 {
221 yMin = pVoxelLimit.GetMinYExtent() ;
222 }
223 if (yMax > pVoxelLimit.GetMaxYExtent())
224 {
225 yMax = pVoxelLimit.GetMaxYExtent() ;
226 }
227 }
228 }
229 zoffset = pTransform.NetTranslation().z() ;
230 zMin = zoffset - fDz ;
231 zMax = zoffset + fDz ;
232
233 if (pVoxelLimit.IsZLimited())
234 {
235 if ( zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance ||
236 zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance ) return false ;
237 else
238 {
239 if (zMin < pVoxelLimit.GetMinZExtent())
240 {
241 zMin = pVoxelLimit.GetMinZExtent() ;
242 }
243 if (zMax > pVoxelLimit.GetMaxZExtent())
244 {
245 zMax = pVoxelLimit.GetMaxZExtent() ;
246 }
247 }
248 }
249 switch (pAxis)
250 {
251 case kXAxis:
252 pMin = xMin ;
253 pMax = xMax ;
254 break ;
255 case kYAxis:
256 pMin=yMin;
257 pMax=yMax;
258 break;
259 case kZAxis:
260 pMin=zMin;
261 pMax=zMax;
262 break;
263 default:
264 break;
265 }
266 pMin -= kCarTolerance ;
267 pMax += kCarTolerance ;
268
269 return true;
270 }
271 else // General rotated case - create and clip mesh to boundaries
272 {
273 G4bool existsAfterClip = false ;
274 G4ThreeVectorList* vertices ;
275
276 pMin = +kInfinity ;
277 pMax = -kInfinity ;
278
279 // Calculate rotated vertex coordinates
280
281 vertices = CreateRotatedVertices(pTransform) ;
282 ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ;
283 ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax) ;
284 ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ;
285
286 if (pVoxelLimit.IsLimited(pAxis) == false)
287 {
288 if ( pMin != kInfinity || pMax != -kInfinity )
289 {
290 existsAfterClip = true ;
291
292 // Add 2*tolerance to avoid precision troubles
293
294 pMin -= kCarTolerance;
295 pMax += kCarTolerance;
296 }
297 }
298 else
299 {
300 G4ThreeVector clipCentre(
301 ( pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
302 ( pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
303 ( pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
304
305 if ( pMin != kInfinity || pMax != -kInfinity )
306 {
307 existsAfterClip = true ;
308
309
310 // Check to see if endpoints are in the solid
311
312 clipCentre(pAxis) = pVoxelLimit.GetMinExtent(pAxis);
313
314 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
315 {
316 pMin = pVoxelLimit.GetMinExtent(pAxis);
317 }
318 else
319 {
320 pMin -= kCarTolerance;
321 }
322 clipCentre(pAxis) = pVoxelLimit.GetMaxExtent(pAxis);
323
324 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
325 {
326 pMax = pVoxelLimit.GetMaxExtent(pAxis);
327 }
328 else
329 {
330 pMax += kCarTolerance;
331 }
332 }
333
334 // Check for case where completely enveloping clipping volume
335 // If point inside then we are confident that the solid completely
336 // envelopes the clipping volume. Hence set min/max extents according
337 // to clipping volume extents along the specified axis.
338
339 else if (Inside(pTransform.Inverse().TransformPoint(clipCentre))
340 != kOutside)
341 {
342 existsAfterClip = true ;
343 pMin = pVoxelLimit.GetMinExtent(pAxis) ;
344 pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
345 }
346 }
347 delete vertices;
348 return existsAfterClip;
349 }
350}
351
352/////////////////////////////////////////////////////////////////////////
353//
354// Return whether point inside/outside/on surface, using tolerance
355
356EInside G4Box::Inside(const G4ThreeVector& p) const
357{
358 EInside in = kOutside ;
359
360 if ( std::fabs(p.x()) <= fDx - kCarTolerance*0.5 )
361 {
362 if (std::fabs(p.y()) <= fDy - kCarTolerance*0.5 )
363 {
364 if (std::fabs(p.z()) <= fDz - kCarTolerance*0.5 ) in = kInside ;
365 else if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5 ) in = kSurface ;
366 }
367 else if (std::fabs(p.y()) <= fDy + kCarTolerance*0.5 )
368 {
369 if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5 ) in = kSurface ;
370 }
371 }
372 else if (std::fabs(p.x()) <= fDx + kCarTolerance*0.5 )
373 {
374 if (std::fabs(p.y()) <= fDy + kCarTolerance*0.5 )
375 {
376 if (std::fabs(p.z()) <= fDz + kCarTolerance*0.5) in = kSurface ;
377 }
378 }
379 return in ;
380}
381
382///////////////////////////////////////////////////////////////////////
383//
384// Calculate side nearest to p, and return normal
385// If two sides are equidistant, normal of first side (x/y/z)
386// encountered returned
387
388G4ThreeVector G4Box::SurfaceNormal( const G4ThreeVector& p) const
389{
390 G4double distx, disty, distz ;
391 G4ThreeVector norm ;
392
393 // Calculate distances as if in 1st octant
394
395 distx = std::fabs(std::fabs(p.x()) - fDx) ;
396 disty = std::fabs(std::fabs(p.y()) - fDy) ;
397 distz = std::fabs(std::fabs(p.z()) - fDz) ;
398
399 // New code for particle on surface including edges and corners with specific
400 // normals
401
402 const G4double delta = 0.5*kCarTolerance;
403 const G4ThreeVector nX = G4ThreeVector( 1.0, 0,0 );
404 const G4ThreeVector nmX = G4ThreeVector(-1.0, 0,0 );
405 const G4ThreeVector nY = G4ThreeVector( 0, 1.0,0 );
406 const G4ThreeVector nmY = G4ThreeVector( 0,-1.0,0 );
407 const G4ThreeVector nZ = G4ThreeVector( 0, 0, 1.0);
408 const G4ThreeVector nmZ = G4ThreeVector( 0, 0,- 1.0);
409
410 G4ThreeVector normX(0.,0.,0.), normY(0.,0.,0.), normZ(0.,0.,0.);
411 G4ThreeVector sumnorm(0., 0., 0.);
412 G4int noSurfaces=0;
413
414 if (distx <= delta) // on X/mX surface and around
415 {
416 noSurfaces ++;
417 if ( p.x() >= 0.){ // on +X surface
418 normX= nX ; // G4ThreeVector( 1.0, 0., 0. );
419 }else{
420 normX= nmX; // G4ThreeVector(-1.0, 0., 0. );
421 }
422 sumnorm= normX;
423 }
424
425 if (disty <= delta) // on one of the +Y or -Y surfaces
426 {
427 noSurfaces ++;
428 if ( p.y() >= 0.){ // on +Y surface
429 normY= nY;
430 }else{
431 normY = nmY;
432 }
433 sumnorm += normY;
434 }
435
436 if (distz <= delta) // on one of the +Z or -Z surfaces
437 {
438 noSurfaces ++;
439 if ( p.z() >= 0.){ // on +Z surface
440 normZ= nZ;
441 }else{
442 normZ = nmZ;
443 }
444 sumnorm += normZ;
445 }
446
447 // sumnorm= normX + normY + normZ;
448 const G4double invSqrt2 = 1.0 / std::sqrt( 2.0);
449 const G4double invSqrt3 = 1.0 / std::sqrt( 3.0);
450
451 norm= G4ThreeVector( 0., 0., 0.);
452 if( noSurfaces > 0 )
453 {
454 if( noSurfaces == 1 ){
455 norm= sumnorm;
456 }else{
457 // norm = sumnorm . unit();
458 if( noSurfaces == 2 ) {
459 // 2 surfaces -> on edge
460 norm = invSqrt2 * sumnorm;
461 } else {
462 // 3 surfaces (on corner)
463 norm = invSqrt3 * sumnorm;
464 }
465 }
466 }else{
467#ifdef G4CSGDEBUG
468 G4Exception("G4Box::SurfaceNormal(p)", "Notification", JustWarning,
469 "Point p is not on surface !?" );
470#endif
471 norm = ApproxSurfaceNormal(p);
472 }
473
474 return norm;
475}
476
477//////////////////////////////////////////////////////////////////////////
478//
479// Algorithm for SurfaceNormal() following the original specification
480// for points not on the surface
481
482G4ThreeVector G4Box::ApproxSurfaceNormal( const G4ThreeVector& p ) const
483{
484 G4double distx, disty, distz ;
485 G4ThreeVector norm ;
486
487 // Calculate distances as if in 1st octant
488
489 distx = std::fabs(std::fabs(p.x()) - fDx) ;
490 disty = std::fabs(std::fabs(p.y()) - fDy) ;
491 distz = std::fabs(std::fabs(p.z()) - fDz) ;
492
493 if ( distx <= disty )
494 {
495 if ( distx <= distz ) // Closest to X
496 {
497 if ( p.x() < 0 ) norm = G4ThreeVector(-1.0,0,0) ;
498 else norm = G4ThreeVector( 1.0,0,0) ;
499 }
500 else // Closest to Z
501 {
502 if ( p.z() < 0 ) norm = G4ThreeVector(0,0,-1.0) ;
503 else norm = G4ThreeVector(0,0, 1.0) ;
504 }
505 }
506 else
507 {
508 if ( disty <= distz ) // Closest to Y
509 {
510 if ( p.y() < 0 ) norm = G4ThreeVector(0,-1.0,0) ;
511 else norm = G4ThreeVector(0, 1.0,0) ;
512 }
513 else // Closest to Z
514 {
515 if ( p.z() < 0 ) norm = G4ThreeVector(0,0,-1.0) ;
516 else norm = G4ThreeVector(0,0, 1.0) ;
517 }
518 }
519 return norm;
520}
521
522///////////////////////////////////////////////////////////////////////////
523//
524// Calculate distance to box from an outside point
525// - return kInfinity if no intersection.
526//
527// ALGORITHM:
528//
529// Check that if point lies outside x/y/z extent of box, travel is towards
530// the box (ie. there is a possibility of an intersection)
531//
532// Calculate pairs of minimum and maximum distances for x/y/z travel for
533// intersection with the box's x/y/z extent.
534// If there is a valid intersection, it is given by the maximum min distance
535// (ie. distance to satisfy x/y/z intersections) *if* <= minimum max distance
536// (ie. distance after which 1+ of x/y/z intersections not satisfied)
537//
538// NOTE:
539//
540// `Inside' safe - meaningful answers given if point is inside the exact
541// shape.
542
543G4double G4Box::DistanceToIn(const G4ThreeVector& p,const G4ThreeVector& v) const
544{
545 G4double safx, safy, safz ;
546 G4double smin=0.0, sminy, sminz ; // , sminx ;
547 G4double smax=kInfinity, smaxy, smaxz ; // , smaxx ; // they always > 0
548 G4double stmp ;
549 G4double sOut=kInfinity, sOuty=kInfinity, sOutz=kInfinity ;
550
551 safx = std::fabs(p.x()) - fDx ; // minimum distance to x surface of shape
552 safy = std::fabs(p.y()) - fDy ;
553 safz = std::fabs(p.z()) - fDz ;
554
555 // Will we intersect?
556 // If safx/y/z is >-tol/2 the point is outside/on the box's x/y/z extent.
557 // If both p.x/y/z and v.x/y/z repectively are both positive/negative,
558 // travel is in a direction away from the shape.
559
560 if ( ((p.x()*v.x() >= 0.0) && safx > -kCarTolerance*0.5)
561 || ((p.y()*v.y() >= 0.0) && safy > -kCarTolerance*0.5)
562 || ((p.z()*v.z() >= 0.0) && safz > -kCarTolerance*0.5) )
563 {
564 return kInfinity ; // travel away or parallel within tolerance
565 }
566
567 // Compute min / max distances for x/y/z travel:
568 // X Planes
569
570 if ( v.x())
571 {
572 stmp = 1.0/std::fabs(v.x()) ;
573
574 if (safx >= 0.0)
575 {
576 smin = safx*stmp ;
577 smax = (fDx+std::fabs(p.x()))*stmp ;
578 }
579 else
580 {
581 if (v.x() > 0) sOut = (fDx - p.x())*stmp ;
582 if (v.x() < 0) sOut = (fDx + p.x())*stmp ;
583 }
584 }
585
586 // Y Planes
587
588 if ( v.y())
589 {
590 stmp = 1.0/std::fabs(v.y()) ;
591
592 if (safy >= 0.0)
593 {
594 sminy = safy*stmp ;
595 smaxy = (fDy+std::fabs(p.y()))*stmp ;
596
597 if (sminy > smin) smin=sminy ;
598 if (smaxy < smax) smax=smaxy ;
599
600 if (smin >= smax-kCarTolerance*0.5)
601 {
602 return kInfinity ; // touch XY corner
603 }
604 }
605 else
606 {
607 if (v.y() > 0) sOuty = (fDy - p.y())*stmp ;
608 if (v.y() < 0) sOuty = (fDy + p.y())*stmp ;
609 if( sOuty < sOut ) sOut = sOuty ;
610 }
611 }
612
613 // Z planes
614
615 if ( v.z() )
616 {
617 stmp = 1.0/std::fabs(v.z()) ;
618
619 if ( safz >= 0.0)
620 {
621 sminz = safz*stmp ;
622 smaxz = (fDz+std::fabs(p.z()))*stmp ;
623
624 if (sminz > smin) smin = sminz ;
625 if (smaxz < smax) smax = smaxz ;
626
627 if (smin >= smax-kCarTolerance*0.5)
628 {
629 return kInfinity ; // touch ZX or ZY corners
630 }
631 }
632 else
633 {
634 if (v.z() > 0) sOutz = (fDz - p.z())*stmp ;
635 if (v.z() < 0) sOutz = (fDz + p.z())*stmp ;
636 if( sOutz < sOut ) sOut = sOutz ;
637 }
638 }
639
640 if ( sOut <= smin + 0.5*kCarTolerance) // travel over edge
641 {
642 return kInfinity ;
643 }
644 if (smin < 0.5*kCarTolerance) smin = 0.0 ;
645
646 return smin ;
647}
648
649//////////////////////////////////////////////////////////////////////////
650//
651// Appoximate distance to box.
652// Returns largest perpendicular distance to the closest x/y/z sides of
653// the box, which is the most fast estimation of the shortest distance to box
654// - If inside return 0
655
656G4double G4Box::DistanceToIn(const G4ThreeVector& p) const
657{
658 G4double safex, safey, safez, safe = 0.0 ;
659
660 safex = std::fabs(p.x()) - fDx ;
661 safey = std::fabs(p.y()) - fDy ;
662 safez = std::fabs(p.z()) - fDz ;
663
664 if (safex > safe) safe = safex ;
665 if (safey > safe) safe = safey ;
666 if (safez > safe) safe = safez ;
667
668 return safe ;
669}
670
671/////////////////////////////////////////////////////////////////////////
672//
673// Calcluate distance to surface of box from inside
674// by calculating distances to box's x/y/z planes.
675// Smallest distance is exact distance to exiting.
676// - Eliminate one side of each pair by considering direction of v
677// - when leaving a surface & v.close, return 0
678
679G4double G4Box::DistanceToOut( const G4ThreeVector& p,const G4ThreeVector& v,
680 const G4bool calcNorm,
681 G4bool *validNorm,G4ThreeVector *n) const
682{
683 ESide side = kUndefined ;
684 G4double pdist,stmp,snxt;
685
686 if (calcNorm) *validNorm = true ; // All normals are valid
687
688 if (v.x() > 0) // X planes
689 {
690 pdist = fDx - p.x() ;
691
692 if (pdist > kCarTolerance*0.5)
693 {
694 snxt = pdist/v.x() ;
695 side = kPX ;
696 }
697 else
698 {
699 if (calcNorm) *n = G4ThreeVector(1,0,0) ;
700 return snxt = 0 ;
701 }
702 }
703 else if (v.x() < 0)
704 {
705 pdist = fDx + p.x() ;
706
707 if (pdist > kCarTolerance*0.5)
708 {
709 snxt = -pdist/v.x() ;
710 side = kMX ;
711 }
712 else
713 {
714 if (calcNorm) *n = G4ThreeVector(-1,0,0) ;
715 return snxt = 0 ;
716 }
717 }
718 else snxt = kInfinity ;
719
720 if ( v.y() > 0 ) // Y planes
721 {
722 pdist=fDy-p.y();
723
724 if (pdist>kCarTolerance*0.5)
725 {
726 stmp=pdist/v.y();
727
728 if (stmp<snxt)
729 {
730 snxt=stmp;
731 side=kPY;
732 }
733 }
734 else
735 {
736 if (calcNorm) *n = G4ThreeVector(0,1,0) ;
737 return snxt = 0 ;
738 }
739 }
740 else if ( v.y() < 0 )
741 {
742 pdist = fDy + p.y() ;
743
744 if (pdist > kCarTolerance*0.5)
745 {
746 stmp=-pdist/v.y();
747
748 if (stmp<snxt)
749 {
750 snxt=stmp;
751 side=kMY;
752 }
753 }
754 else
755 {
756 if (calcNorm) *n = G4ThreeVector(0,-1,0) ;
757 return snxt = 0 ;
758 }
759 }
760 if (v.z()>0) // Z planes
761 {
762 pdist=fDz-p.z();
763
764 if (pdist > kCarTolerance*0.5)
765 {
766 stmp=pdist/v.z();
767
768 if (stmp < snxt)
769 {
770 snxt=stmp;
771 side=kPZ;
772 }
773 }
774 else
775 {
776 if (calcNorm) *n = G4ThreeVector(0,0,1) ;
777 return snxt = 0 ;
778 }
779 }
780 else if (v.z()<0)
781 {
782 pdist = fDz + p.z() ;
783
784 if (pdist > kCarTolerance*0.5)
785 {
786 stmp=-pdist/v.z();
787
788 if (stmp < snxt)
789 {
790 snxt=stmp;
791 side=kMZ;
792 }
793 }
794 else
795 {
796 if (calcNorm) *n = G4ThreeVector(0,0,-1) ;
797 return snxt = 0 ;
798 }
799 }
800 if (calcNorm)
801 {
802 switch (side)
803 {
804 case kPX:
805 *n=G4ThreeVector(1,0,0);
806 break;
807 case kMX:
808 *n=G4ThreeVector(-1,0,0);
809 break;
810 case kPY:
811 *n=G4ThreeVector(0,1,0);
812 break;
813 case kMY:
814 *n=G4ThreeVector(0,-1,0);
815 break;
816 case kPZ:
817 *n=G4ThreeVector(0,0,1);
818 break;
819 case kMZ:
820 *n=G4ThreeVector(0,0,-1);
821 break;
822 default:
823 G4cout.precision(16);
824 G4cout << G4endl;
825 DumpInfo();
826 G4cout << "Position:" << G4endl << G4endl;
827 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl;
828 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl;
829 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl;
830 G4cout << "Direction:" << G4endl << G4endl;
831 G4cout << "v.x() = " << v.x() << G4endl;
832 G4cout << "v.y() = " << v.y() << G4endl;
833 G4cout << "v.z() = " << v.z() << G4endl << G4endl;
834 G4cout << "Proposed distance :" << G4endl << G4endl;
835 G4cout << "snxt = " << snxt/mm << " mm" << G4endl << G4endl;
836 G4Exception("G4Box::DistanceToOut(p,v,..)","Notification",JustWarning,
837 "Undefined side for valid surface normal to solid.");
838 break;
839 }
840 }
841 return snxt;
842}
843
844////////////////////////////////////////////////////////////////////////////
845//
846// Calculate exact shortest distance to any boundary from inside
847// - If outside return 0
848
849G4double G4Box::DistanceToOut(const G4ThreeVector& p) const
850{
851 G4double safx1,safx2,safy1,safy2,safz1,safz2,safe=0.0;
852
853#ifdef G4CSGDEBUG
854 if( Inside(p) == kOutside )
855 {
856 G4cout.precision(16) ;
857 G4cout << G4endl ;
858 DumpInfo();
859 G4cout << "Position:" << G4endl << G4endl ;
860 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
861 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
862 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
863 G4Exception("G4Box::DistanceToOut(p)", "Notification", JustWarning,
864 "Point p is outside !?" );
865 }
866#endif
867
868 safx1 = fDx - p.x() ;
869 safx2 = fDx + p.x() ;
870 safy1 = fDy - p.y() ;
871 safy2 = fDy + p.y() ;
872 safz1 = fDz - p.z() ;
873 safz2 = fDz + p.z() ;
874
875 // shortest Dist to any boundary now MIN(safx1,safx2,safy1..)
876
877 if (safx2 < safx1) safe = safx2 ;
878 else safe = safx1 ;
879 if (safy1 < safe) safe = safy1 ;
880 if (safy2 < safe) safe = safy2 ;
881 if (safz1 < safe) safe = safz1 ;
882 if (safz2 < safe) safe = safz2 ;
883
884 if (safe < 0) safe = 0 ;
885 return safe ;
886}
887
888////////////////////////////////////////////////////////////////////////
889//
890// Create a List containing the transformed vertices
891// Ordering [0-3] -fDz cross section
892// [4-7] +fDz cross section such that [0] is below [4],
893// [1] below [5] etc.
894// Note:
895// Caller has deletion resposibility
896
897G4ThreeVectorList*
898G4Box::CreateRotatedVertices(const G4AffineTransform& pTransform) const
899{
900 G4ThreeVectorList* vertices = new G4ThreeVectorList();
901 vertices->reserve(8);
902
903 if (vertices)
904 {
905 G4ThreeVector vertex0(-fDx,-fDy,-fDz) ;
906 G4ThreeVector vertex1(fDx,-fDy,-fDz) ;
907 G4ThreeVector vertex2(fDx,fDy,-fDz) ;
908 G4ThreeVector vertex3(-fDx,fDy,-fDz) ;
909 G4ThreeVector vertex4(-fDx,-fDy,fDz) ;
910 G4ThreeVector vertex5(fDx,-fDy,fDz) ;
911 G4ThreeVector vertex6(fDx,fDy,fDz) ;
912 G4ThreeVector vertex7(-fDx,fDy,fDz) ;
913
914 vertices->push_back(pTransform.TransformPoint(vertex0));
915 vertices->push_back(pTransform.TransformPoint(vertex1));
916 vertices->push_back(pTransform.TransformPoint(vertex2));
917 vertices->push_back(pTransform.TransformPoint(vertex3));
918 vertices->push_back(pTransform.TransformPoint(vertex4));
919 vertices->push_back(pTransform.TransformPoint(vertex5));
920 vertices->push_back(pTransform.TransformPoint(vertex6));
921 vertices->push_back(pTransform.TransformPoint(vertex7));
922 }
923 else
924 {
925 DumpInfo();
926 G4Exception("G4Box::CreateRotatedVertices()",
927 "FatalError", FatalException,
928 "Error in allocation of vertices. Out of memory !");
929 }
930 return vertices;
931}
932
933//////////////////////////////////////////////////////////////////////////
934//
935// GetEntityType
936
937G4GeometryType G4Box::GetEntityType() const
938{
939 return G4String("G4Box");
940}
941
942//////////////////////////////////////////////////////////////////////////
943//
944// Stream object contents to an output stream
945
946std::ostream& G4Box::StreamInfo(std::ostream& os) const
947{
948 os << "-----------------------------------------------------------\n"
949 << " *** Dump for solid - " << GetName() << " ***\n"
950 << " ===================================================\n"
951 << " Solid type: G4Box\n"
952 << " Parameters: \n"
953 << " half length X: " << fDx/mm << " mm \n"
954 << " half length Y: " << fDy/mm << " mm \n"
955 << " half length Z: " << fDz/mm << " mm \n"
956 << "-----------------------------------------------------------\n";
957
958 return os;
959}
960
961/////////////////////////////////////////////////////////////////////////////////////
962//
963// GetPointOnSurface
964//
965// Return a point (G4ThreeVector) randomly and uniformly selected
966// on the solid surface
967
968G4ThreeVector G4Box::GetPointOnSurface() const
969{
970 G4double px, py, pz, select, sumS;
971 G4double Sxy = fDx*fDy, Sxz = fDx*fDz, Syz = fDy*fDz;
972
973 sumS = Sxy + Sxz + Syz;
974 select = sumS*G4UniformRand();
975
976 if( select < Sxy )
977 {
978 px = -fDx +2*fDx*G4UniformRand();
979 py = -fDy +2*fDy*G4UniformRand();
980
981 if(G4UniformRand() > 0.5) pz = fDz;
982 else pz = -fDz;
983 }
984 else if ( ( select - Sxy ) < Sxz )
985 {
986 px = -fDx +2*fDx*G4UniformRand();
987 pz = -fDz +2*fDz*G4UniformRand();
988
989 if(G4UniformRand() > 0.5) py = fDy;
990 else py = -fDy;
991 }
992 else
993 {
994 py = -fDy +2*fDy*G4UniformRand();
995 pz = -fDz +2*fDz*G4UniformRand();
996
997 if(G4UniformRand() > 0.5) px = fDx;
998 else px = -fDx;
999 }
1000 return G4ThreeVector(px,py,pz);
1001}
1002
1003//////////////////////////////////////////////////////////////////////////
1004//
1005// Methods for visualisation
1006
1007void G4Box::DescribeYourselfTo (G4VGraphicsScene& scene) const
1008{
1009 scene.AddSolid (*this);
1010}
1011
1012G4VisExtent G4Box::GetExtent() const
1013{
1014 return G4VisExtent (-fDx, fDx, -fDy, fDy, -fDz, fDz);
1015}
1016
1017G4Polyhedron* G4Box::CreatePolyhedron () const
1018{
1019 return new G4PolyhedronBox (fDx, fDy, fDz);
1020}
1021
1022G4NURBS* G4Box::CreateNURBS () const
1023{
1024 return new G4NURBSbox (fDx, fDy, fDz);
1025}
Note: See TracBrowser for help on using the repository browser.