Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4SphericalSurface.cc
Go to the documentation of this file.
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$
28//
29// ----------------------------------------------------------------------
30// GEANT 4 class source file
31//
32// G4SphericalSurface.cc
33//
34// ----------------------------------------------------------------------
35
36#include "G4SphericalSurface.hh"
38
40 : G4Surface(), radius(1.0), phi_1(0.0), phi_2(2*pi),
41 theta_1(0.0), theta_2(pi)
42{
43 // default constructor
44 // default x_axis is ( 1.0, 0.0, 0.0 ), z_axis is ( 0.0, 0.0, 1.0 ),
45 // default radius is 1.0
46 // default phi_1 is 0, phi_2 is 2*PI
47 // default theta_1 is 0, theta_2 is PI
48
49 x_axis = G4Vector3D( 1.0, 0.0, 0.0 );
50 z_axis = G4Vector3D( 0.0, 0.0, 1.0 );
51 // OuterBoundary = new G4BREPPolyline();
52}
53
54
56 const G4Vector3D& xhat,
57 const G4Vector3D& zhat,
58 G4double r,
59 G4double ph1, G4double ph2,
60 G4double th1, G4double th2)
61 : G4Surface()
62{
63 // Require both x_axis and z_axis to be unit vectors
64
65 G4double xhatmag = xhat.mag();
66 if ( xhatmag != 0.0 )
67 {
68 x_axis = xhat * (1/ xhatmag); // this makes the x_axis a unit vector
69 }
70 else
71 {
72 std::ostringstream message;
73 message << "x_axis has zero length." << G4endl
74 << "Default x_axis of (1, 0, 0) is used.";
75 G4Exception("G4SphericalSurface::G4SphericalSurface()",
76 "GeomSolids1001", JustWarning, message);
77
78 x_axis = G4Vector3D( 1.0, 0.0, 0.0 );
79 }
80
81 G4double zhatmag = zhat.mag();
82
83 if (zhatmag != 0.0)
84 {
85 z_axis = zhat *(1/ zhatmag); // this makes the z_axis a unit vector
86 }
87 else
88 {
89 std::ostringstream message;
90 message << "z_axis has zero length." << G4endl
91 << "Default z_axis of (0, 0, 1) is used.";
92 G4Exception("G4SphericalSurface::G4SphericalSurface()",
93 "GeomSolids1001", JustWarning, message);
94
95 z_axis = G4Vector3D( 0.0, 0.0, 1.0 );
96 }
97
98 // Require radius to be non-negative
99 //
100 if ( r >= 0.0 )
101 {
102 radius = r;
103 }
104 else
105 {
106 std::ostringstream message;
107 message << "Radius cannot be less than zero." << G4endl
108 << "Default radius of 1.0 is used.";
109 G4Exception("G4SphericalSurface::G4SphericalSurface()",
110 "GeomSolids1001", JustWarning, message);
111
112 radius = 1.0;
113 }
114
115 // Require phi_1 in the range: 0 <= phi_1 < 2*PI
116 // and phi_2 in the range: phi_1 < phi_2 <= phi_1 + 2*PI
117 //
118 if ( ( ph1 >= 0.0 ) && ( ph1 < 2*pi ) )
119 {
120 phi_1 = ph1;
121 }
122 else
123 {
124 std::ostringstream message;
125 message << "Lower azimuthal limit is out of range." << G4endl
126 << "Default angle of 0 is used.";
127 G4Exception("G4SphericalSurface::G4SphericalSurface()",
128 "GeomSolids1001", JustWarning, message);
129
130 phi_1 = 0.0;
131 }
132
133 if ( ( ph2 > phi_1 ) && ( ph2 <= ( phi_1 + twopi ) ) )
134 {
135 phi_2 = ph2;
136 }
137 else
138 {
139 std::ostringstream message;
140 message << "Upper azimuthal limit is out of range." << G4endl
141 << "Default angle of 2*PI is used.";
142 G4Exception("G4SphericalSurface::G4SphericalSurface()",
143 "GeomSolids1001", JustWarning, message);
144
145 phi_2 = twopi;
146 }
147
148 // Require theta_1 in the range: 0 <= theta_1 < PI
149 // and theta-2 in the range: theta_1 < theta_2 <= theta_1 + PI
150 //
151 if ( ( th1 >= 0.0 ) && ( th1 < pi ) )
152 {
153 theta_1 = th1;
154 }
155 else
156 {
157 std::ostringstream message;
158 message << "Lower polar limit is out of range." << G4endl
159 << "Default angle of 0 is used.";
160 G4Exception("G4SphericalSurface::G4SphericalSurface()",
161 "GeomSolids1001", JustWarning, message);
162
163 theta_1 = 0.0;
164 }
165
166 if ( ( th2 > theta_1 ) && ( th2 <= ( theta_1 + pi ) ) )
167 {
168 theta_2 =th2;
169 }
170 else
171 {
172 std::ostringstream message;
173 message << "Upper polar limit is out of range." << G4endl
174 << "Default angle of PI is used.";
175 G4Exception("G4SphericalSurface::G4SphericalSurface()",
176 "GeomSolids1001", JustWarning, message);
177
178 theta_2 = pi;
179 }
180}
181
182
184{
185}
186
187
189 : G4Surface()
190{
191 x_axis = surf.x_axis;
192 z_axis = surf.z_axis;
193 radius = surf.radius;
194 phi_1 = surf.phi_1;
195 phi_2 = surf.phi_2;
196 theta_1 = surf.theta_1;
197 theta_2 = surf.theta_2;
198}
199
201G4SphericalSurface::operator=( const G4SphericalSurface& surf )
202{
203 if (&surf == this) { return *this; }
204 x_axis = surf.x_axis;
205 z_axis = surf.z_axis;
206 radius = surf.radius;
207 phi_1 = surf.phi_1;
208 phi_2 = surf.phi_2;
209 theta_1 = surf.theta_1;
210 theta_2 = surf.theta_2;
211
212 return *this;
213}
214
215const char* G4SphericalSurface::NameOf() const
216{
217 return "G4SphericalSurface";
218}
219
220void G4SphericalSurface::PrintOn( std::ostream& os ) const
221{
222 os << "G4SphericalSurface surface with origin: " << origin << "\t"
223 << "radius: " << radius << "\n"
224 << "\t local x_axis: " << x_axis
225 << "\t local z_axis: " << z_axis << "\n"
226 << "\t lower azimuthal limit: " << phi_1 << " radians\n"
227 << "\t upper azimuthal limit: " << phi_2 << " radians\n"
228 << "\t lower polar limit : " << theta_1 << " radians\n"
229 << "\t upper polar limit : " << theta_2 << " radians\n";
230}
231
232
234{
235 // Distance from the point x to the G4SphericalSurface.
236 // The distance will be positive if the point is Inside the
237 // G4SphericalSurface, negative if the point is outside.
238
239 G4Vector3D d = G4Vector3D( x - origin );
240 G4double rds = d.mag();
241
242 return (radius - rds);
243}
244
245
246/*
247G4double G4SphericalSurface::distanceAlongRay( G4int which_way,
248 const G4Ray* ry,
249 G4Vector3D& p ) const
250
251 // Distance along a Ray (straight line with G4Vector3D) to leave or enter
252 // a G4SphericalSurface. The input variable which_way should be set to +1 to
253 // indicate leaving a G4SphericalSurface, -1 to indicate entering the surface.
254 // p is the point of intersection of the Ray with the G4SphericalSurface.
255 // If the G4Vector3D of the Ray is opposite to that of the Normal to
256 // the G4SphericalSurface at the intersection point, it will not leave the
257 // G4SphericalSurface.
258 // Similarly, if the G4Vector3D of the Ray is along that of the Normal
259 // to the G4SphericalSurface at the intersection point, it will not enter the
260 // G4SphericalSurface.
261 // This method is called by all finite shapes sub-classed to
262 // G4SphericalSurface.
263 // Use the virtual function table to check if the intersection point
264 // is within the boundary of the finite shape.
265 // A negative result means no intersection.
266 // If no valid intersection point is found, set the distance
267 // and intersection point to large numbers.
268
269{
270 G4double Dist = kInfinity;
271 G4Vector3D lv ( kInfinity, kInfinity, kInfinity );
272 p = lv;
273
274 // Origin and G4Vector3D unit vector of Ray.
275 //
276 G4Vector3D x = ry->Position( 0.0 );
277 G4Vector3D dhat = ry->Direction( 0.0 );
278 G4int isoln = 0, maxsoln = 2;
279
280 // Array of solutions in distance along the Ray
281 //
282 G4double s[2];s[0] = -1.0; s[1]= -1.0 ;
283
284 // Calculate the two solutions (quadratic equation)
285 //
286 G4Vector3D d = x - GetOrigin();
287 G4double radius = GetRadius();
288
289 // Quit with no intersection if the radius of the G4SphericalSurface is zero
290 //
291 if ( radius <= 0.0 ) { return Dist; }
292
293 G4double dsq = d * d;
294 G4double rsq = radius * radius;
295 G4double b = d * dhat;
296 G4double c = dsq - rsq;
297 G4double radical = b * b - c;
298
299 // Quit with no intersection if the radical is negative
300 //
301 if ( radical < 0.0 ) { return Dist; }
302
303 G4double root = std::sqrt( radical );
304 s[0] = -b + root;
305 s[1] = -b - root;
306
307 // Order the possible solutions by increasing distance along the Ray
308 //
309 G4Sort_double( s, isoln, maxsoln-1 );
310
311 // Now loop over each positive solution, keeping the first one (smallest
312 // distance along the Ray) which is within the boundary of the sub-shape
313 // and which also has the correct G4Vector3D with respect to the Normal to
314 // the G4SphericalSurface at the intersection point
315 //
316 for ( isoln = 0; isoln < maxsoln; isoln++ )
317 {
318 if ( s[isoln] >= 0.0 )
319 {
320 if ( s[isoln] >= kInfinity ) { return Dist; } // quit if too large
321 Dist = s[isoln];
322 p = ry->Position( Dist );
323 if ( ( ( dhat * Normal( p ) * which_way ) >= 0.0 )
324 && ( WithinBoundary( p ) == 1 ) ) { return Dist; }
325 }
326 }
327
328 // Get here only if there was no solution within the boundary,
329 // reset distance and intersection point to large numbers
330
331 p = lv;
332 return kInfinity;
333}
334*/
335
336
338{
339 G4double x_min = origin.x() - radius;
340 G4double y_min = origin.y() - radius;
341 G4double z_min = origin.z() - radius;
342 G4double x_max = origin.x() + radius;
343 G4double y_max = origin.y() + radius;
344 G4double z_max = origin.z() + radius;
345
346 G4Point3D Min(x_min, y_min, z_min);
347 G4Point3D Max(x_max, y_max, z_max);
348 bbox = new G4BoundingBox3D( Min, Max);
349}
350
351
353{
354 // Distance along a Ray (straight line with G4Vector3D) to leave or enter
355 // a G4SphericalSurface. The input variable which_way should be set to +1
356 // to indicate leaving a G4SphericalSurface, -1 to indicate entering a
357 // G4SphericalSurface.
358 // p is the point of intersection of the Ray with the G4SphericalSurface.
359 // If the G4Vector3D of the Ray is opposite to that of the Normal to
360 // the G4SphericalSurface at the intersection point, it will not leave the
361 // G4SphericalSurface.
362 // Similarly, if the G4Vector3D of the Ray is along that of the Normal
363 // to the G4SphericalSurface at the intersection point, it will not enter
364 // the G4SphericalSurface.
365 // This method is called by all finite shapes sub-classed to
366 // G4SphericalSurface.
367 // Use the virtual function table to check if the intersection point
368 // is within the boundary of the finite shape.
369 // A negative result means no intersection.
370 // If no valid intersection point is found, set the distance
371 // and intersection point to large numbers.
372
373 G4int which_way = (G4int)HowNear(ry.GetStart());
374
375 if(!which_way) { which_way =-1; }
376 distance = kInfinity;
377
378 G4Vector3D lv ( kInfinity, kInfinity, kInfinity );
379
380 // p = lv;
381 //
382 closest_hit = lv;
383
384 // Origin and G4Vector3D unit vector of Ray.
385 //
386 G4Vector3D x= G4Vector3D( ry.GetStart() );
387 G4Vector3D dhat = ry.GetDir();
388 G4int isoln = 0, maxsoln = 2;
389
390 // Array of solutions in distance along the Ray
391 //
392 G4double ss[2];
393 ss[0] = -1.0 ;
394 ss[1] = -1.0 ;
395
396 // Calculate the two solutions (quadratic equation)
397 //
398 G4Vector3D d = G4Vector3D( x - GetOrigin() );
399 G4double r = GetRadius();
400
401 // Quit with no intersection if the radius of the G4SphericalSurface is zero
402 //
403 if ( r <= 0.0 ) { return 0; }
404
405 G4double dsq = d * d;
406 G4double rsq = r * r;
407 G4double b = d * dhat;
408 G4double c = dsq - rsq;
409 G4double radical = b * b - c;
410
411 // Quit with no intersection if the radical is negative
412 //
413 if ( radical < 0.0 ) { return 0; }
414
415 G4double root = std::sqrt( radical );
416 ss[0] = -b + root;
417 ss[1] = -b - root;
418
419 // Order the possible solutions by increasing distance along the Ray
420 // G4Sort_double( ss, isoln, maxsoln-1 );
421 //
422 if(ss[0] > ss[1])
423 {
424 G4double tmp =ss[0];
425 ss[0] = ss[1];
426 ss[1] = tmp;
427 }
428
429 // Now loop over each positive solution, keeping the first one (smallest
430 // distance along the Ray) which is within the boundary of the sub-shape
431 // and which also has the correct G4Vector3D with respect to the Normal to
432 // the G4SphericalSurface at the intersection point
433 //
434 for ( isoln = 0; isoln < maxsoln; isoln++ )
435 {
436 if ( ss[isoln] >= kCarTolerance*0.5 )
437 {
438 if ( ss[isoln] >= kInfinity ) { return 0; } // quit if too large
439 distance = ss[isoln];
441 if ( ( ( dhat * Normal( closest_hit ) * which_way ) >= 0.0 ) &&
442 ( WithinBoundary( closest_hit ) == 1 ) )
443 {
445 return 1;
446 }
447 }
448 }
449
450 // Get here only if there was no solution within the boundary,
451 // reset distance and intersection point to large numbers
452 //
453 distance = kInfinity;
454 closest_hit = lv;
455
456 return 0;
457}
458
459
460/*
461G4double G4SphericalSurface::distanceAlongHelix( G4int which_way,
462 const Helix* hx,
463 G4Vector3D& p ) const
464{
465 // Distance along a Helix to leave or enter a G4SphericalSurface.
466 // The input variable which_way should be set to +1 to
467 // indicate leaving a G4SphericalSurface, -1 to indicate entering a
468 // G4SphericalSurface.
469 // p is the point of intersection of the Helix with the G4SphericalSurface.
470 // If the G4Vector3D of the Helix is opposite to that of the Normal to
471 // the G4SphericalSurface at the intersection point, it will not leave the
472 // G4SphericalSurface.
473 // Similarly, if the G4Vector3D of the Helix is along that of the Normal
474 // to the G4SphericalSurface at the intersection point, it will not enter
475 // the G4SphericalSurface.
476 // This method is called by all finite shapes sub-classed to
477 // G4SphericalSurface.
478 // Use the virtual function table to check if the intersection point
479 // is within the boundary of the finite shape.
480 // If no valid intersection point is found, set the distance
481 // and intersection point to large numbers.
482 // Possible negative distance solutions are discarded.
483
484{
485 G4double Dist = kInfinity;
486 G4Vector3D lv ( kInfinity, kInfinity, kInfinity );
487 p = lv;
488 G4int isoln = 0, maxsoln = 4;
489
490 // Array of solutions in turning angle
491 //
492 G4double s[4];s[0] = -1.0; s[1]= -1.0 ;s[2] = -1.0; s[3]= -1.0 ;
493
494 // Helix parameters
495 //
496 G4double rh = hx->GetRadius(); // radius of Helix
497 G4Vector3D oh = hx->position( 0.0 ); // origin of Helix
498 G4Vector3D dh = hx->direction( 0.0 ); // initial G4Vector3D of Helix
499 G4Vector3D prp = hx->getPerp(); // perpendicular vector
500 G4double prpmag = prp.mag();
501 G4double rhp = rh / prpmag;
502 G4double rs = GetRadius(); // radius of G4SphericalSurface
503 if ( rs == 0.0 ) { return Dist; } // quit if zero radius
504 G4Vector3D os = GetOrigin(); // origin of G4SphericalSurface
505
506 // Calculate quantities of use later on
507 //
508 G4Vector3D alpha = rhp * prp;
509 G4Vector3D beta = rhp * dh;
510 G4Vector3D gamma = oh - os;
511
512 // Only consider approximate solutions to quadratic order in the turning
513 // angle along the Helix
514 //
515 G4double A = beta * beta + gamma * alpha;
516 G4double B = 2.0 * gamma * beta;
517 G4double C = gamma * gamma - rs * rs;
518
519 // Case if quadratic term is zero
520 //
521 if ( std::fabs( A ) < kCarTolerance )
522 {
523 if ( B == 0.0 ) { return Dist; } // no intersection, quit
524 else { s[0] = -C / B; }
525 }
526 else // General quadratic solution, A != 0
527 {
528 G4double radical = B * B - 4.0 * A * C;
529 if ( radical < 0.0 ) { return Dist; } // no intersection, quit
530
531 G4double root = std::sqrt( radical );
532 s[0] = ( -B + root ) / ( 2.0 * A );
533 s[1] = ( -B - root ) / ( 2.0 * A );
534 if ( rh < 0.0 )
535 {
536 s[0] = -s[0];
537 s[1] = -s[1];
538 }
539 s[2] = s[0] + twopi;
540 s[3] = s[1] + twopi;
541 }
542
543 // Order the possible solutions by increasing turning angle
544 //
545 G4Sort_double( s, isoln, maxsoln-1 );
546
547 // Now loop over each positive solution, keeping the first one (smallest
548 // distance along the Helix) which is within the boundary of the sub-shape.
549 //
550 for ( isoln = 0; isoln < maxsoln; isoln++ )
551 {
552 if ( s[isoln] >= 0.0 )
553 {
554 // Calculate distance along Helix and position and G4Vector3D vectors.
555 //
556 Dist = s[isoln] * std::fabs( rhp );
557 p = hx->position( Dist );
558 G4Vector3D d = hx->direction( Dist );
559
560 // Now do approximation to get remaining distance to correct this
561 // solution iterate it until the accuracy is below the user-set
562 // surface precision
563 //
564 G4double delta = 0.;
565 G4double delta0 = kInfinity;
566 G4int dummy = 1;
567 G4int iter = 0;
568 G4int in0 = Inside( hx->position ( 0.0 ) );
569 G4int in1 = Inside( p );
570 G4double sc = Scale();
571 while ( dummy )
572 {
573 iter++;
574
575 // Terminate loop after 50 iterations and Reset distance to large
576 // number, indicating no intersection with G4SphericalSurface. This
577 // generally occurs if the Helix curls too tightly to intersect it.
578 //
579 if ( iter > 50 )
580 {
581 Dist = kInfinity;
582 p = lv;
583 break;
584 }
585
586 // Find distance from the current point along the above-calculated
587 // G4Vector3D using a Ray.
588 // The G4Vector3D of the Ray and the Sign of the distance are
589 // determined by whether the starting point of the Helix is Inside
590 // or outside of the G4SphericalSurface.
591 //
592 in1 = Inside( p );
593 if ( in1 ) // current point Inside
594 {
595 if ( in0 ) // starting point Inside
596 {
597 Ray* r = new Ray( p, d );
598 delta = distanceAlongRay( 1, r, p );
599 delete r;
600 }
601 else // starting point outside
602 {
603 Ray* r = new Ray( p, -d );
604 delta = -distanceAlongRay( 1, r, p );
605 delete r;
606 }
607 }
608 else // current point outside
609 {
610 if ( in0 ) // starting point Inside
611 {
612 Ray* r = new Ray( p, -d );
613 delta = -distanceAlongRay( -1, r, p );
614 delete r;
615 }
616 else // starting point outside
617 {
618 Ray* r = new Ray( p, d );
619 delta = distanceAlongRay( -1, r, p );
620 delete r;
621 }
622 }
623
624 // Test if distance is less than the surface precision, if so
625 // terminate loop.
626 //
627 if ( std::fabs( delta / sc ) <= SURFACE_PRECISION ) { break; }
628
629 // If delta has not changed sufficiently from the previous iteration,
630 // skip out of this loop.
631 //
632 if (std::fabs( ( delta-delta0 )/sc ) <= SURFACE_PRECISION) { break; }
633
634 // If delta has increased in absolute value from the previous iteration
635 // either the Helix doesn't Intersect the G4SphericalSurface or the
636 // approximate solution is too far from the real solution. Try groping
637 // for a solution. If not found, Reset distance to large number,
638 // indicating no intersection with the G4SphericalSurface.
639 //
640 if ( ( std::fabs( delta ) > std::fabs( delta0 ) ) )
641 {
642 Dist = std::fabs( rhp ) * gropeAlongHelix( hx );
643 if ( Dist < 0.0 )
644 {
645 Dist = kInfinity;
646 p = lv;
647 }
648 else
649 {
650 p = hx->position( Dist );
651 }
652 break;
653 }
654
655 // Set old delta to new one.
656 //
657 delta0 = delta;
658
659 // Add distance to G4SphericalSurface to distance along Helix.
660 //
661 Dist += delta;
662
663 // Negative distance along Helix means Helix doesn't Intersect
664 // G4SphericalSurface. Reset distance to large number, indicating
665 // no intersection with G4SphericalSurface.
666 //
667 if ( Dist < 0.0 )
668 {
669 Dist = kInfinity;
670 p = lv;
671 break;
672 }
673
674 // Recalculate point along Helix and the G4Vector3D.
675 //
676 p = hx->position( Dist );
677 d = hx->direction( Dist );
678 } // end of while loop
679
680 // Now have best value of distance along Helix and position for this
681 // solution, so test if it is within the boundary of the sub-shape
682 // and require that it point in the correct G4Vector3D with respect to
683 // the Normal to the G4SphericalSurface.
684
685 if ( ( Dist < kInfinity ) &&
686 ( ( hx->direction( Dist ) * Normal( p ) * which_way ) >= 0.0 )
687 && ( WithinBoundary( p ) == 1 ) ) { return Dist; }
688 } // end of if s[isoln] >= 0.0 condition
689 } // end of for loop over solutions
690
691 // If one gets here, there is no solution, so set distance along Helix
692 // and position to large numbers.
693 //
694 Dist = kInfinity;
695 p = lv;
696
697 return Dist;
698}
699
700
701G4Vector3D G4SphericalSurface::Normal( const G4Vector3D& p ) const
702{
703 // Return the Normal unit vector to the G4SphericalSurface at a point p on
704 // (or nearly on) the G4SphericalSurface.
705
706 G4Vector3D n = p - origin;
707 G4double nmag = n.mag();
708
709 // If the point p happens to coincide with the origin (which is possible
710 // if the radius is zero), set the Normal to the z-axis unit vector.
711 //
712 if ( nmag != 0.0 ) { n = n / nmag; }
713 else { n = G4Vector3D( 0.0, 0.0, 1.0 ); }
714
715 return n;
716}
717*/
718
719
721{
722 // Return the Normal unit vector to the G4SphericalSurface at a point p on
723 // (or nearly on) the G4SphericalSurface.
724
725 G4Vector3D n = G4Vector3D( p - origin );
726 G4double nmag = n.mag();
727
728 // If the point p happens to coincide with the origin (which is possible
729 // if the radius is zero), set the Normal to the z-axis unit vector.
730 //
731 if ( nmag != 0.0 ) { n = n * (1/ nmag); }
732 else { n = G4Vector3D( 0.0, 0.0, 1.0 ); }
733
734 return n;
735}
736
737
739{
740 // Return the Normal unit vector to the G4SphericalSurface at a point p on
741 // (or nearly on) the G4SphericalSurface.
742
743 G4Vector3D n = G4Vector3D( p - origin );
744 G4double nmag = n.mag();
745
746 // If the point p happens to coincide with the origin (which is possible
747 // if the radius is zero), set the Normal to the z-axis unit vector.
748 //
749 if ( nmag != 0.0 ) { n = n * (1/ nmag); }
750 else { n = G4Vector3D( 0.0, 0.0, 1.0 ); }
751
752 return n;
753}
754
755
757{
758 // Return 0 if point x is outside G4SphericalSurface, 1 if Inside.
759 // Outside means that the distance to the G4SphericalSurface would
760 // be negative.
761 // Use the HowNear function to calculate this distance.
762
763 if ( HowNear( x ) >= 0.0 ) { return 1; }
764 else { return 0; }
765}
766
767
769{
770 // Return 1 if point x is on the G4SphericalSurface, otherwise return zero
771 // (x is assumed to lie on the surface of the G4SphericalSurface, so one
772 // only checks the angular limits)
773
774 G4Vector3D y_axis = G4Vector3D( z_axis.cross( x_axis ) );
775
776 // Components of x in the local coordinate system of the G4SphericalSurface
777 //
778 G4double px = x * x_axis;
779 G4double py = x * y_axis;
780 G4double pz = x * z_axis;
781
782 // Check if within polar angle limits
783 //
784 G4double theta = std::acos( pz / x.mag() ); // acos in range 0 to PI
785
786 // Normal case
787 //
788 if ( theta_2 <= pi )
789 {
790 if ( ( theta < theta_1 ) || ( theta > theta_2 ) ) { return 0; }
791 }
792 else
793 {
794 theta += pi;
795 if ( ( theta < theta_1 ) || ( theta > theta_2 ) ) { return 0; }
796 }
797
798 // Now check if within azimuthal angle limits
799 //
800 G4double phi = std::atan2( py, px ); // atan2 in range -PI to PI
801
802 if ( phi < 0.0 ) { phi += twopi; }
803
804 // Normal case
805 //
806 if ( ( phi >= phi_1 ) && ( phi <= phi_2 ) ) { return 1; }
807
808 // This is for the case that phi_2 is greater than 2*PI
809 //
810 phi += twopi;
811
812 if ( ( phi >= phi_1 ) && ( phi <= phi_2 ) ) { return 1; }
813
814 // Get here if not within azimuthal limits
815
816 return 0;
817}
818
819
821{
822 // Returns the radius of a G4SphericalSurface unless it is zero, in which
823 // case returns the arbitrary number 1.0.
824 // Used for Scale-invariant tests of surface thickness.
825
826 if ( radius == 0.0 ) { return 1.0; }
827 else { return radius; }
828}
829
830
832{
833 // Returns the Area of a G4SphericalSurface
834
835 return ( 2.0*( theta_2 - theta_1 )*( phi_2 - phi_1)*radius*radius/pi );
836}
837
838
840 G4double ph1, G4double ph2,
841 G4double th1, G4double th2 )
842{
843 // Resize the G4SphericalSurface to new radius r, new lower and upper
844 // azimuthal angle limits ph1 and ph2, and new lower and upper polar
845 // angle limits th1 and th2.
846
847 // Require radius to be non-negative
848 //
849 if ( r >= 0.0 ) { radius = r; }
850 else
851 {
852 std::ostringstream message;
853 message << "Radius cannot be less than zero." << G4endl
854 << "Original value of " << radius << " is retained.";
855 G4Exception("G4SphericalSurface::resize()",
856 "GeomSolids1001", JustWarning, message);
857 }
858
859 // Require azimuthal angles to be within bounds
860
861 if ( ( ph1 >= 0.0 ) && ( ph1 < twopi ) ) { phi_1 = ph1; }
862 else
863 {
864 std::ostringstream message;
865 message << "Lower azimuthal limit out of range." << G4endl
866 << "Original value of " << phi_1 << " is retained.";
867 G4Exception("G4SphericalSurface::resize()",
868 "GeomSolids1001", JustWarning, message);
869 }
870
871 if ( ( ph2 > phi_1 ) && ( ph2 <= ( phi_1 + twopi ) ) ) { phi_2 = ph2; }
872 else
873 {
874 ph2 = ( phi_2 <= phi_1 ) ? ( phi_1 + kAngTolerance ) : phi_2;
875 phi_2 = ph2;
876 std::ostringstream message;
877 message << "Upper azimuthal limit out of range." << G4endl
878 << "Value of " << phi_2 << " is used.";
879 G4Exception("G4SphericalSurface::resize()",
880 "GeomSolids1001", JustWarning, message);
881 }
882
883 // Require polar angles to be within bounds
884 //
885 if ( ( th1 >= 0.0 ) && ( th1 < pi ) ) { theta_1 = th1; }
886 else
887 {
888 std::ostringstream message;
889 message << "Lower polar limit out of range." << G4endl
890 << "Original value of " << theta_1 << " is retained.";
891 G4Exception("G4SphericalSurface::resize()",
892 "GeomSolids1001", JustWarning, message);
893 }
894
895 if ( ( th2 > theta_1 ) && ( th2 <= ( theta_1 + pi ) ) ) { theta_2 = th2; }
896 else
897 {
898 th2 = ( theta_2 <= theta_1 ) ? ( theta_1 + kAngTolerance ) : theta_2;
899 theta_2 = th2;
900 std::ostringstream message;
901 message << "Upper polar limit out of range." << G4endl
902 << "Value of " << theta_2 << " is used.";
903 G4Exception("G4SphericalSurface::resize()",
904 "GeomSolids1001", JustWarning, message);
905 }
906}
907
908
909/*
910void G4SphericalSurface::
911rotate( G4double alpha, G4double beta,
912 G4double gamma, G4ThreeMat& m, G4int inverse )
913{
914 // Rotate G4SphericalSurface first about global x_axis by angle alpha,
915 // second about global y-axis by angle beta,
916 // and third about global z_axis by angle gamma
917 // by creating and using G4ThreeMat objects in Surface::rotate
918 // angles are assumed to be given in radians
919 // if inverse is non-zero, the order of rotations is reversed
920 // the axis is rotated here, the origin is rotated by calling
921 // Surface::rotate
922
923 G4Surface::rotate( alpha, beta, gamma, m, inverse );
924 x_axis = m * x_axis;
925 z_axis = m * z_axis;
926}
927
928
929void G4SphericalSurface::
930rotate( G4double alpha, G4double beta, G4double gamma, G4int inverse )
931{
932 // Rotate G4SphericalSurface first about global x_axis by angle alpha,
933 // second about global y-axis by angle beta,
934 // and third about global z_axis by angle gamma
935 // by creating and using G4ThreeMat objects in Surface::rotate
936 // angles are assumed to be given in radians
937 // if inverse is non-zero, the order of rotations is reversed
938 // the axis is rotated here, the origin is rotated by calling
939 // Surface::rotate
940
941 G4ThreeMat m;
942 G4Surface::rotate( alpha, beta, gamma, m, inverse );
943 x_axis = m * x_axis;
944 z_axis = m * z_axis;
945}
946*/
947
948
949/*
950G4double G4SphericalSurface::gropeAlongHelix( const Helix* hx ) const
951{
952 // Grope for a solution of a Helix intersecting a G4SphericalSurface.
953 // This function returns the turning angle (in radians) where the
954 // intersection occurs with only positive values allowed, or -1.0 if
955 // no intersection is found.
956 // The idea is to start at the beginning of the Helix, then take steps
957 // of some fraction of a turn. If at the end of a Step, the current position
958 // along the Helix and the previous position are on opposite sides of the
959 // G4SphericalSurface, then the solution must lie somewhere in between.
960
961 G4int one_over_f = 8; // one over fraction of a turn to go in each Step
962 G4double turn_angle = 0.0;
963 G4double dist_along = 0.0;
964 G4double d_new;
965 G4double fk = 1.0 / G4double( one_over_f );
966 G4double scal = Scale();
967 G4double d_old = HowNear( hx->position( dist_along ) );
968 G4double rh = hx->GetRadius(); // radius of Helix
969 G4Vector3D prp = hx->getPerp(); // perpendicular vector
970 G4double prpmag = prp.mag();
971 G4double rhp = rh / prpmag;
972 G4int max_iter = one_over_f * HELIX_MAX_TURNS;
973
974 // Take up to a user-settable number of turns along the Helix,
975 // groping for an intersection point.
976
977 for ( G4int k = 1; k < max_iter; k++ )
978 {
979 turn_angle = twopi * k / one_over_f;
980 dist_along = turn_angle * std::fabs( rhp );
981 d_new = HowNear( hx->position( dist_along ) );
982 if ( ( d_old < 0.0 && d_new > 0.0 ) || ( d_old > 0.0 && d_new < 0.0 ) )
983 {
984 d_old = d_new;
985
986 // Old and new points are on opposite sides of the G4SphericalSurface,
987 // therefore a solution lies in between, use a binary search to pin the
988 // point down to the surface precision, but don't do more than 50
989 // iterations.
990
991 G4int itr = 0;
992 while ( std::fabs( d_new / scal ) > SURFACE_PRECISION )
993 {
994 itr++;
995 if ( itr > 50 ) { return turn_angle; }
996 turn_angle -= fk * pi;
997 dist_along = turn_angle * std::fabs( rhp );
998 d_new = HowNear( hx->position( dist_along ) );
999 if (( d_old < 0.0 && d_new > 0.0 ) || ( d_old > 0.0 && d_new < 0.0 ))
1000 { fk *= -0.5; }
1001 else
1002 { fk *= 0.5; }
1003 d_old = d_new;
1004 } // end of while loop
1005 return turn_angle; // this is the best solution
1006 } // end of if condition
1007 } // end of for loop
1008
1009 // Get here only if no solution is found, so return -1.0 to indicate that.
1010
1011 return -1.0;
1012}
1013*/
@ JustWarning
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
HepGeom::Vector3D< G4double > G4Vector3D
Definition: G4Vector3D.hh:35
#define G4endl
Definition: G4ios.hh:52
Definition: G4Ray.hh:49
G4Point3D GetPoint(G4double i) const
const G4Vector3D & GetDir() const
const G4Point3D & GetStart() const
virtual G4int Inside(const G4Vector3D &x) const
virtual G4double Scale() const
G4double GetRadius() const
virtual G4Vector3D Normal(const G4Vector3D &p) const
virtual void PrintOn(std::ostream &os=G4cout) const
virtual G4int WithinBoundary(const G4Vector3D &x) const
virtual const char * NameOf() const
virtual void resize(G4double r, G4double ph1, G4double ph2, G4double th1, G4double th2)
G4int Intersect(const G4Ray &)
virtual G4double HowNear(const G4Vector3D &x) const
virtual G4double Area() const
virtual G4Vector3D SurfaceNormal(const G4Point3D &p) const
G4double kAngTolerance
Definition: G4Surface.hh:192
G4Vector3D origin
Definition: G4Surface.hh:197
G4double distance
Definition: G4Surface.hh:203
G4BoundingBox3D * bbox
Definition: G4Surface.hh:185
G4Vector3D GetOrigin() const
G4Point3D closest_hit
Definition: G4Surface.hh:186
G4double kCarTolerance
Definition: G4Surface.hh:192
BasicVector3D< T > cross(const BasicVector3D< T > &v) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41