source: trunk/source/geometry/navigation/src/G4ReplicaNavigation.cc@ 881

Last change on this file since 881 was 850, checked in by garnier, 17 years ago

geant4.8.2 beta
File size: 32.2 KB
RevLine 
[831]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//
[850]27// $Id: G4ReplicaNavigation.cc,v 1.19 2008/04/28 15:39:55 gcosmo Exp $
28// GEANT4 tag $Name: HEAD $
[831]29//
30//
31// class G4ReplicaNavigation Implementation
32//
33// Author: P.Kent, 1996
34//
35// --------------------------------------------------------------------
36
37#include "G4ReplicaNavigation.hh"
38
39#include "G4AffineTransform.hh"
40#include "G4SmartVoxelProxy.hh"
41#include "G4SmartVoxelNode.hh"
42#include "G4VSolid.hh"
43#include "G4GeometryTolerance.hh"
44
45#include <assert.h>
46
47// ********************************************************************
48// Constructor
49// ********************************************************************
50//
51G4ReplicaNavigation::G4ReplicaNavigation()
52 : fCheck(false), fVerbose(0)
53{
54 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
55 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
56 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
57}
58
59// ********************************************************************
60// Destructor
61// ********************************************************************
62//
63G4ReplicaNavigation::~G4ReplicaNavigation()
64{
65}
66
67// ********************************************************************
68// Inside
69// ********************************************************************
70//
71EInside
72G4ReplicaNavigation::Inside(const G4VPhysicalVolume *pVol,
73 const G4int replicaNo,
74 const G4ThreeVector &localPoint) const
75{
76 EInside in = kOutside;
77
78 // Replication data
79 //
80 EAxis axis;
81 G4int nReplicas;
82 G4double width, offset;
83 G4bool consuming;
84
85 G4double coord, rad2, rmin, tolRMax2, rmax, tolRMin2;
86
87 pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
88 assert(consuming);
89
90 switch (axis)
91 {
92 case kXAxis:
93 case kYAxis:
94 case kZAxis:
95 coord = std::fabs(localPoint(axis))-width*0.5;
96 if ( coord<=-kCarTolerance*0.5 )
97 {
98 in = kInside;
99 }
100 else if ( coord<=kCarTolerance*0.5 )
101 {
102 in = kSurface;
103 }
104 break;
105 case kPhi:
106 if ( localPoint.y()||localPoint.x() )
107 {
108 coord = std::fabs(std::atan2(localPoint.y(),localPoint.x()))-width*0.5;
109 if ( coord<=-kAngTolerance*0.5 )
110 {
111 in = kInside;
112 }
113 else if ( coord<=kAngTolerance*0.5 )
114 {
115 in = kSurface;
116 }
117 }
118 else
119 {
120 in = kSurface;
121 }
122 break;
123 case kRho:
124 rad2 = localPoint.perp2();
125 rmax = (replicaNo+1)*width+offset;
126 tolRMax2 = rmax-kRadTolerance*0.5;
127 tolRMax2 *= tolRMax2;
128 if ( rad2>tolRMax2 )
129 {
130 tolRMax2 = rmax+kRadTolerance*0.5;
131 tolRMax2 *= tolRMax2;
132 if ( rad2<=tolRMax2 )
133 {
134 in = kSurface;
135 }
136 }
137 else
138 {
139 // Known to be inside outer radius
140 //
141 if ( replicaNo||offset )
142 {
143 rmin = rmax-width;
144 tolRMin2 = rmin-kRadTolerance*0.5;
145 tolRMin2 *= tolRMin2;
146 if ( rad2>tolRMin2 )
147 {
148 tolRMin2 = rmin+kRadTolerance*0.5;
149 tolRMin2 *= tolRMin2;
150 if ( rad2>=tolRMin2 )
151 {
152 in = kInside;
153 }
154 else
155 {
156 in = kSurface;
157 }
158 }
159 }
160 else
161 {
162 in = kInside;
163 }
164 }
165 break;
166 default:
167 G4Exception("G4ReplicaNavigation::Inside()", "WrongArgumentValue",
168 FatalException, "Unknown axis!");
169 break;
170 }
171 return in;
172}
173
174// ********************************************************************
175// DistanceToOut
176// ********************************************************************
177//
178G4double
179G4ReplicaNavigation::DistanceToOut(const G4VPhysicalVolume *pVol,
180 const G4int replicaNo,
181 const G4ThreeVector &localPoint) const
182{
183 // Replication data
184 //
185 EAxis axis;
186 G4int nReplicas;
187 G4double width,offset;
188 G4bool consuming;
189
190 G4double safety=0.;
191 G4double safe1,safe2;
192 G4double coord, rho, rmin, rmax;
193
194 pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
195 assert(consuming);
196 switch(axis)
197 {
198 case kXAxis:
199 case kYAxis:
200 case kZAxis:
201 coord = localPoint(axis);
202 safe1 = width*0.5-coord;
203 safe2 = width*0.5+coord;
204 safety = (safe1<=safe2) ? safe1 : safe2;
205 break;
206 case kPhi:
207 if ( localPoint.y()<=0 )
208 {
209 safety = localPoint.x()*std::sin(width*0.5)
210 + localPoint.y()*std::cos(width*0.5);
211 }
212 else
213 {
214 safety = localPoint.x()*std::sin(width*0.5)
215 - localPoint.y()*std::cos(width*0.5);
216 }
217 break;
218 case kRho:
219 rho = localPoint.perp();
220 rmax = width*(replicaNo+1)+offset;
221 if ( replicaNo||offset )
222 {
223 rmin = rmax-width;
224 safe1 = rho-rmin;
225 safe2 = rmax-rho;
226 safety = (safe1<=safe2) ? safe1 : safe2;
227 }
228 else
229 {
230 safety = rmax-rho;
231 }
232 break;
233 default:
234 G4Exception("G4ReplicaNavigation::DistanceToOut()", "WrongArgumentValue",
235 FatalException, "Unknown axis!");
236 break;
237 }
238 return (safety >= kCarTolerance) ? safety : 0;
239}
240
241// ********************************************************************
242// DistanceToOut
243// ********************************************************************
244//
245G4double
246G4ReplicaNavigation::DistanceToOut(const G4VPhysicalVolume *pVol,
247 const G4int replicaNo,
248 const G4ThreeVector &localPoint,
249 const G4ThreeVector &localDirection) const
250{
251 // Replication data
252 //
253 EAxis axis;
254 G4int nReplicas;
255 G4double width, offset;
256 G4bool consuming;
257
258 G4double Dist=kInfinity;
259 G4double coord, Comp, lindist;
260
261 pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
262 assert(consuming);
263 switch(axis)
264 {
265 case kXAxis:
266 case kYAxis:
267 case kZAxis:
268 coord = localPoint(axis);
269 Comp = localDirection(axis);
270 if ( Comp>0 )
271 {
272 lindist = width*0.5-coord;
273 Dist = (lindist>0) ? lindist/Comp : 0;
274 }
275 else if ( Comp<0 )
276 {
277 lindist = width*0.5+coord;
278 Dist = (lindist>0) ? -lindist/Comp : 0;
279 }
280 else
281 {
282 Dist = kInfinity;
283 }
284 break;
285 case kPhi:
286 Dist = DistanceToOutPhi(localPoint, localDirection, width);
287 break;
288 case kRho:
289 Dist=DistanceToOutRad(localPoint,localDirection,width,offset,replicaNo);
290 break;
291 default:
292 G4Exception("G4ReplicaNavigation::DistanceToOut()", "WrongArgumentValue",
293 FatalException, "Unknown axis!");
294 break;
295 }
296 return Dist;
297}
298
299// ********************************************************************
300// DistanceToOutPhi
301// ********************************************************************
302//
303G4double
304G4ReplicaNavigation::DistanceToOutPhi(const G4ThreeVector &localPoint,
305 const G4ThreeVector &localDirection,
306 const G4double width) const
307{
308 // Phi Intersection
309 // NOTE: width<=pi by definition
310 //
311 G4double sinSPhi, cosSPhi;
312 G4double pDistS, pDistE, compS, compE, Dist, dist2, yi;
313 if ( localPoint.x()||localPoint.y() )
314 {
315 sinSPhi = std::sin(-width*0.5); // SIN of starting phi plane
316 cosSPhi = std::cos(width*0.5); // COS of starting phi plane
317
318 // pDist -ve when inside
319 //
320 pDistS = localPoint.x()*sinSPhi-localPoint.y()*cosSPhi;
321 pDistE = localPoint.x()*sinSPhi+localPoint.y()*cosSPhi;
322
323 // Comp -ve when in direction of outwards normal
324 //
325 compS = -sinSPhi*localDirection.x()+cosSPhi*localDirection.y();
326 compE = -sinSPhi*localDirection.x()-cosSPhi*localDirection.y();
327
328 if ( (pDistS<=0)&&(pDistE<=0) )
329 {
330 // Inside both phi *full* planes
331 //
332 if ( compS<0 )
333 {
334 dist2 = pDistS/compS;
335 yi = localPoint.y()+dist2*localDirection.y();
336
337 // Check intersecting with correct half-plane (no -> no intersect)
338 //
339 if ( yi<=0 )
340 {
341 Dist = (pDistS<=-kCarTolerance*0.5) ? dist2 : 0;
342 }
343 else
344 {
345 Dist = kInfinity;
346 }
347 }
348 else
349 {
350 Dist = kInfinity;
351 }
352 if ( compE<0 )
353 {
354 dist2 = pDistE/compE;
355
356 // Only check further if < starting phi intersection
357 //
358 if ( dist2<Dist )
359 {
360 yi = localPoint.y()+dist2*localDirection.y();
361
362 // Check intersecting with correct half-plane
363 //
364 if ( yi>=0 )
365 {
366 // Leaving via ending phi
367 //
368 Dist = (pDistE<=-kCarTolerance*0.5) ? dist2 : 0;
369 }
370 }
371 }
372 }
373 else if ( (pDistS>=0)&&(pDistE>=0) )
374 {
375 // Outside both *full* phi planes
376 // if towards both >=0 then once inside will remain inside
377 //
378 Dist = ((compS>=0)&&(compE>=0)) ? kInfinity : 0;
379 }
380 else if ( (pDistS>0)&&(pDistE<0) )
381 {
382 // Outside full starting plane, inside full ending plane
383 //
384 if ( compE<0 )
385 {
386 dist2 = pDistE/compE;
387 yi = localPoint.y()+dist2*localDirection.y();
388
389 // Check intersection in correct half-plane
390 // (if not -> remain in extent)
391 //
392 Dist = (yi>0) ? dist2 : kInfinity;
393 }
394 else // Leaving immediately by starting phi
395 {
396 Dist = kInfinity;
397 }
398 }
399 else
400 {
401 // Must be (pDistS<0)&&(pDistE>0)
402 // Inside full starting plane, outside full ending plane
403 //
404 if ( compE>=0 )
405 {
406 if ( compS<0 )
407 {
408 dist2 = pDistS/compS;
409 yi = localPoint.y()+dist2*localDirection.y();
410
411 // Check intersection in correct half-plane
412 // (if not -> remain in extent)
413 //
414 Dist = (yi<0) ? dist2 : kInfinity;
415 }
416 else
417 {
418 Dist = kInfinity;
419 }
420 }
421 else
422 {
423 // Leaving immediately by ending phi
424 //
425 Dist = 0;
426 }
427 }
428 }
429 else
430 {
431 // On z axis + travel not || to z axis -> use direction vector
432 //
433 Dist = (std::fabs(localDirection.phi())<=width*0.5) ? kInfinity : 0;
434 }
435 return Dist;
436}
437
438// ********************************************************************
439// DistanceToOutRad
440// ********************************************************************
441//
442G4double
443G4ReplicaNavigation::DistanceToOutRad(const G4ThreeVector &localPoint,
444 const G4ThreeVector &localDirection,
445 const G4double width,
446 const G4double offset,
447 const G4int replicaNo) const
448{
449 G4double rmin, rmax, t1, t2, t3, deltaR;
450 G4double b, c, d2, sr;
451
452 //
453 // Radial Intersections
454 //
455
456 // Find intersction with cylinders at rmax/rmin
457 // Intersection point (xi,yi,zi) on line
458 // x=localPoint.x+t*localDirection.x etc.
459 //
460 // Intersects with x^2+y^2=R^2
461 //
462 // Hence (localDirection.x^2+localDirection.y^2)t^2+
463 // 2t(localPoint.x*localDirection.x+localPoint.y*localDirection.y)+
464 // localPoint.x^2+localPoint.y^2-R^2=0
465 //
466 // t1 t2 t3
467
468 rmin = replicaNo*width+offset;
469 rmax = (replicaNo+1)*width+offset;
470
471 t1 = 1.0-localDirection.z()*localDirection.z(); // since v normalised
472 t2 = localPoint.x()*localDirection.x()+localPoint.y()*localDirection.y();
473 t3 = localPoint.x()*localPoint.x()+localPoint.y()*localPoint.y();
474
475 if ( t1>0 ) // Check not parallel
476 {
477 // Calculate sr, r exit distance
478 //
479 if ( t2>=0 )
480 {
481 // Delta r not negative => leaving via rmax
482 //
483 deltaR = t3-rmax*rmax;
484
485 // NOTE: Should use
486 // rho-rmax<-kRadTolerance*0.5 - [no sqrts for efficiency]
487 //
488 if ( deltaR<-kRadTolerance*0.5 )
489 {
490 b = t2/t1;
491 c = deltaR/t1;
492 sr = -b+std::sqrt(b*b-c);
493 }
494 else
495 {
496 // On tolerant boundary & heading outwards (or locally
497 // perpendicular to) outer radial surface -> leaving immediately
498 //
499 sr = 0;
500 }
501 }
502 else
503 {
504 // Possible rmin intersection
505 //
506 if (rmin)
507 {
508 deltaR = t3-rmin*rmin;
509 b = t2/t1;
510 c = deltaR/t1;
511 d2 = b*b-c;
512 if ( d2>=0 )
513 {
514 // Leaving via rmin
515 // NOTE: Should use
516 // rho-rmin>kRadTolerance*0.5 - [no sqrts for efficiency]
517 //
518 sr = (deltaR>kRadTolerance*0.5) ? -b-std::sqrt(d2) : 0;
519 }
520 else
521 {
522 // No rmin intersect -> must be rmax intersect
523 //
524 deltaR = t3-rmax*rmax;
525 c = deltaR/t1;
526 sr = -b+std::sqrt(b*b-c);
527 }
528 }
529 else
530 {
531 // No rmin intersect -> must be rmax intersect
532 //
533 deltaR = t3-rmax*rmax;
534 b = t2/t1;
535 c = deltaR/t1;
536 sr = -b+std::sqrt(b*b-c);
537 }
538 }
539 }
540 else
541 {
542 sr=kInfinity;
543 }
544 return sr;
545}
546
547// ********************************************************************
548// ComputeTransformation
549//
550// Setup transformation and transform point into local system
551// ********************************************************************
552//
553void
554G4ReplicaNavigation::ComputeTransformation(const G4int replicaNo,
555 G4VPhysicalVolume* pVol,
556 G4ThreeVector& point) const
557{
558 G4double val,cosv,sinv,tmpx,tmpy;
559
560 // Replication data
561 //
562 EAxis axis;
563 G4int nReplicas;
564 G4double width,offset;
565 G4bool consuming;
566
567 pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
568 assert(consuming);
569
570 switch (axis)
571 {
572 case kXAxis:
573 val = -width*0.5*(nReplicas-1)+width*replicaNo;
574 pVol->SetTranslation(G4ThreeVector(val,0,0));
575 point.setX(point.x()-val);
576 break;
577 case kYAxis:
578 val = -width*0.5*(nReplicas-1)+width*replicaNo;
579 pVol->SetTranslation(G4ThreeVector(0,val,0));
580 point.setY(point.y()-val);
581 break;
582 case kZAxis:
583 val = -width*0.5*(nReplicas-1)+width*replicaNo;
584 pVol->SetTranslation(G4ThreeVector(0,0,val));
585 point.setZ(point.z()-val);
586 break;
587 case kPhi:
588 val = -(offset+width*(replicaNo+0.5));
589 SetPhiTransformation(val,pVol);
590 cosv = std::cos(val);
591 sinv = std::sin(val);
592 tmpx = point.x()*cosv-point.y()*sinv;
593 tmpy = point.x()*sinv+point.y()*cosv;
594 point.setY(tmpy);
595 point.setX(tmpx);
596 break;
597 case kRho:
598 // No setup required for radial case
599 default:
600 break;
601 }
602}
603
604// ********************************************************************
605// ComputeTransformation
606//
607// Setup transformation into local system
608// ********************************************************************
609//
610void
611G4ReplicaNavigation::ComputeTransformation(const G4int replicaNo,
612 G4VPhysicalVolume* pVol) const
613{
614 G4double val;
615
616 // Replication data
617 //
618 EAxis axis;
619 G4int nReplicas;
620 G4double width, offset;
621 G4bool consuming;
622
623 pVol->GetReplicationData(axis, nReplicas, width, offset, consuming);
624 assert(consuming);
625
626 switch (axis)
627 {
628 case kXAxis:
629 val = -width*0.5*(nReplicas-1)+width*replicaNo;
630 pVol->SetTranslation(G4ThreeVector(val,0,0));
631 break;
632 case kYAxis:
633 val = -width*0.5*(nReplicas-1)+width*replicaNo;
634 pVol->SetTranslation(G4ThreeVector(0,val,0));
635 break;
636 case kZAxis:
637 val = -width*0.5*(nReplicas-1)+width*replicaNo;
638 pVol->SetTranslation(G4ThreeVector(0,0,val));
639 break;
640 case kPhi:
641 val = -(offset+width*(replicaNo+0.5));
642 SetPhiTransformation(val,pVol);
643 break;
644 case kRho:
645 // No setup required for radial case
646 default:
647 break;
648 }
649}
650
651// ********************************************************************
652// ComputeStep
653// ********************************************************************
654//
655G4double
656G4ReplicaNavigation::ComputeStep(const G4ThreeVector &globalPoint,
657 const G4ThreeVector &globalDirection,
658 const G4ThreeVector &localPoint,
659 const G4ThreeVector &localDirection,
660 const G4double currentProposedStepLength,
661 G4double &newSafety,
662 G4NavigationHistory &history,
663 G4bool &validExitNormal,
664 G4ThreeVector &exitNormal,
665 G4bool &exiting,
666 G4bool &entering,
667 G4VPhysicalVolume *(*pBlockedPhysical),
668 G4int &blockedReplicaNo )
669{
670 G4VPhysicalVolume *repPhysical, *motherPhysical;
671 G4VPhysicalVolume *samplePhysical, *blockedExitedVol=0;
672 G4LogicalVolume *repLogical;
673 G4VSolid *motherSolid;
674 G4ThreeVector repPoint, repDirection, sampleDirection;
675 G4double ourStep=currentProposedStepLength;
676 G4double ourSafety=kInfinity;
677 G4double sampleStep, sampleSafety, motherStep, motherSafety;
678 G4int localNoDaughters, sampleNo;
679 G4int depth;
680
681 // Exiting normal optimisation
682 //
683 if ( exiting&&validExitNormal )
684 {
685 if ( localDirection.dot(exitNormal)>=kMinExitingNormalCosine )
686 {
687 // Block exited daughter volume
688 //
689 blockedExitedVol = *pBlockedPhysical;
690 ourSafety = 0;
691 }
692 }
693 exiting = false;
694 entering = false;
695
696 repPhysical = history.GetTopVolume();
697 repLogical = repPhysical->GetLogicalVolume();
698
699 //
700 // Compute intersection with replica boundaries & replica safety
701 //
702
703 sampleSafety = DistanceToOut(repPhysical,
704 history.GetTopReplicaNo(),
705 localPoint);
706
707 if ( sampleSafety<ourSafety )
708 {
709 ourSafety = sampleSafety;
710 }
711 if ( sampleSafety<ourStep )
712 {
713 sampleStep = DistanceToOut(repPhysical,
714 history.GetTopReplicaNo(),
715 localPoint,
716 localDirection);
717 if ( sampleStep<ourStep )
718 {
719 ourStep = sampleStep;
720 exiting = true;
721 validExitNormal = false;
722 }
723 }
724
725 depth = history.GetDepth()-1;
726 while ( history.GetVolumeType(depth)==kReplica )
727 {
728 repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
729 sampleSafety = DistanceToOut(history.GetVolume(depth),
730 history.GetReplicaNo(depth),
731 repPoint);
732 if ( sampleSafety<ourSafety )
733 {
734 ourSafety = sampleSafety;
735 }
736 if ( sampleSafety<ourStep )
737 {
738 sampleStep = DistanceToOut(history.GetVolume(depth),
739 history.GetReplicaNo(depth),
740 repPoint,
741 history.GetTransform(depth).TransformAxis(globalDirection));
742 if ( sampleStep<ourStep )
743 {
744 ourStep = sampleStep;
745 exiting = true;
746 validExitNormal = false;
747 }
748 }
749 depth--;
750 }
751
752 // Compute mother safety & intersection
753 //
754 repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
755 motherPhysical = history.GetVolume(depth);
756 motherSolid = motherPhysical->GetLogicalVolume()->GetSolid();
757 motherSafety = motherSolid->DistanceToOut(repPoint);
758 repDirection = history.GetTransform(depth).TransformAxis(globalDirection);
759 motherStep = motherSolid->DistanceToOut(repPoint,repDirection,true,
760 &validExitNormal,&exitNormal);
761
762 // Push in principle no longer necessary. G4Navigator now takes care of ...
763 // Removing this will however generate warnings for pushed particles from
764 // G4Navigator, particularly for the case of 3D replicas (Cartesian or
765 // combined Radial/Phi cases).
766 // Requires further investigation and eventually reimplementation of
767 // LevelLocate() to take into account point and direction ...
768 //
769 if ( ( !ourStep && (sampleSafety<0.5*kCarTolerance) )
770 && ( repLogical->GetSolid()->Inside(localPoint)==kSurface ) )
771 {
772 ourStep += kCarTolerance;
773 }
774
775 if ( motherSafety<ourSafety )
776 {
777 ourSafety = motherSafety;
778 }
779
780#ifdef G4VERBOSE
781 if ( fCheck )
782 {
783 if( motherSolid->Inside(localPoint)==kOutside )
784 {
785 G4cout << "WARNING - G4ReplicaNavigation::ComputeStep()" << G4endl
786 << " Point " << localPoint
787 << " is outside current volume " << motherPhysical->GetName()
788 << G4endl;
789 G4double estDistToSolid= motherSolid->DistanceToIn(localPoint);
790 G4cout << " Estimated isotropic distance to solid (distToIn)= "
791 << estDistToSolid << G4endl;
792 if( estDistToSolid > 100.0 * kCarTolerance )
793 {
794 motherSolid->DumpInfo();
795 G4Exception("G4ReplicaNavigation::ComputeStep()",
796 "FarOutsideCurrentVolume", FatalException,
797 "Point is far outside Current Volume !" );
798 }
799 else
800 G4Exception("G4ReplicaNavigation::ComputeStep()",
801 "OutsideCurrentVolume", JustWarning,
802 "Point is a little outside Current Volume.");
803 }
804 }
805#endif
806
807 // May need precision protection
808 //
809 if ( motherSafety<=ourStep )
810 {
811 if ( motherStep<=ourStep )
812 {
813 ourStep = motherStep;
814 exiting = true;
815 if ( validExitNormal )
816 {
817 const G4RotationMatrix* rot = motherPhysical->GetRotation();
818 if ( rot )
819 {
820 exitNormal *= rot->inverse();
821 }
822 }
823 }
824 else
825 {
826 validExitNormal = false;
827 }
828 }
829 //
830 // Compute daughter safeties & intersections
831 //
832 localNoDaughters = repLogical->GetNoDaughters();
833 for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
834 {
835 samplePhysical = repLogical->GetDaughter(sampleNo);
836 if ( samplePhysical!=blockedExitedVol )
837 {
838 G4AffineTransform sampleTf(samplePhysical->GetRotation(),
839 samplePhysical->GetTranslation());
840 sampleTf.Invert();
841 const G4ThreeVector samplePoint =
842 sampleTf.TransformPoint(localPoint);
843 const G4VSolid* sampleSolid =
844 samplePhysical->GetLogicalVolume()->GetSolid();
845 const G4double sampleSafetyDistance =
846 sampleSolid->DistanceToIn(samplePoint);
847 if ( sampleSafetyDistance<ourSafety )
848 {
849 ourSafety = sampleSafetyDistance;
850 }
851 if ( sampleSafetyDistance<=ourStep )
852 {
853 sampleDirection = sampleTf.TransformAxis(localDirection);
854 const G4double sampleStepDistance =
855 sampleSolid->DistanceToIn(samplePoint,sampleDirection);
856 if ( sampleStepDistance<=ourStep )
857 {
858 ourStep = sampleStepDistance;
859 entering = true;
860 exiting = false;
861 *pBlockedPhysical = samplePhysical;
862 blockedReplicaNo = -1;
863#ifdef G4VERBOSE
864 // Check to see that the resulting point is indeed in/on volume.
865 // This check could eventually be made only for successful candidate.
866
867 if ( ( fCheck ) && ( sampleStepDistance < kInfinity ) )
868 {
869 G4ThreeVector intersectionPoint;
870 intersectionPoint= samplePoint
871 + sampleStepDistance * sampleDirection;
872 EInside insideIntPt= sampleSolid->Inside(intersectionPoint);
873 if ( insideIntPt != kSurface )
874 {
875 G4int oldcoutPrec = G4cout.precision(16);
876 G4cout << "WARNING - G4ReplicaNavigation::ComputeStep()"
877 << G4endl
878 << " Inaccurate DistanceToIn for solid "
879 << sampleSolid->GetName() << G4endl;
880 G4cout << " Solid gave DistanceToIn = "
881 << sampleStepDistance << " yet returns " ;
882 if ( insideIntPt == kInside )
883 G4cout << "-kInside-";
884 else if ( insideIntPt == kOutside )
885 G4cout << "-kOutside-";
886 else
887 G4cout << "-kSurface-";
888 G4cout << " for this point !" << G4endl;
889 G4cout << " Point = " << intersectionPoint << G4endl;
890 if ( insideIntPt != kInside )
891 G4cout << " DistanceToIn(p) = "
892 << sampleSolid->DistanceToIn(intersectionPoint)
893 << G4endl;
894 if ( insideIntPt != kOutside )
895 G4cout << " DistanceToOut(p) = "
896 << sampleSolid->DistanceToOut(intersectionPoint)
897 << G4endl;
898 G4Exception("G4ReplicaNavigation::ComputeStep()",
899 "InaccurateDistanceToIn", JustWarning,
900 "Navigator gets conflicting response from Solid.");
901 G4cout.precision(oldcoutPrec);
902 }
903 }
904#endif
905 }
906 }
907 }
908 }
909 newSafety = ourSafety;
910 return ourStep;
911}
912
913// ********************************************************************
914// ComputeSafety
915//
916// Compute the isotropic distance to current volume's boundaries
917// and to daughter volumes.
918// ********************************************************************
919//
920G4double
921G4ReplicaNavigation::ComputeSafety(const G4ThreeVector &globalPoint,
922 const G4ThreeVector &localPoint,
923 G4NavigationHistory &history,
924 const G4double )
925{
926 G4VPhysicalVolume *repPhysical, *motherPhysical;
927 G4VPhysicalVolume *samplePhysical, *blockedExitedVol=0;
928 G4LogicalVolume *repLogical;
929 G4VSolid *motherSolid;
930 G4ThreeVector repPoint;
931 G4double ourSafety=kInfinity;
932 G4double sampleSafety;
933 G4int localNoDaughters, sampleNo;
934 G4int depth;
935
936 repPhysical = history.GetTopVolume();
937 repLogical = repPhysical->GetLogicalVolume();
938
939 //
940 // Compute intersection with replica boundaries & replica safety
941 //
942
943 sampleSafety = DistanceToOut(history.GetTopVolume(),
944 history.GetTopReplicaNo(),
945 localPoint);
946 if ( sampleSafety<ourSafety )
947 {
948 ourSafety = sampleSafety;
949 }
950
951 depth = history.GetDepth()-1;
952 while ( history.GetVolumeType(depth)==kReplica )
953 {
954 repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
955 sampleSafety = DistanceToOut(history.GetVolume(depth),
956 history.GetReplicaNo(depth),
957 repPoint);
958 if ( sampleSafety<ourSafety )
959 {
960 ourSafety = sampleSafety;
961 }
962 depth--;
963 }
964
965 // Compute mother safety & intersection
966 //
967 repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
968 motherPhysical = history.GetVolume(depth);
969 motherSolid = motherPhysical->GetLogicalVolume()->GetSolid();
970 sampleSafety = motherSolid->DistanceToOut(repPoint);
971
972 if ( sampleSafety<ourSafety )
973 {
974 ourSafety = sampleSafety;
975 }
976
977 // Compute daughter safeties & intersections
978 //
979 localNoDaughters = repLogical->GetNoDaughters();
980 for ( sampleNo=localNoDaughters-1; sampleNo>=0; sampleNo-- )
981 {
982 samplePhysical = repLogical->GetDaughter(sampleNo);
983 if ( samplePhysical!=blockedExitedVol )
984 {
985 G4AffineTransform sampleTf(samplePhysical->GetRotation(),
986 samplePhysical->GetTranslation());
987 sampleTf.Invert();
988 const G4ThreeVector samplePoint =
989 sampleTf.TransformPoint(localPoint);
990 const G4VSolid *sampleSolid =
991 samplePhysical->GetLogicalVolume()->GetSolid();
992 const G4double sampleSafetyDistance =
993 sampleSolid->DistanceToIn(samplePoint);
994 if ( sampleSafetyDistance<ourSafety )
995 {
996 ourSafety = sampleSafetyDistance;
997 }
998 }
999 }
1000 return ourSafety;
1001}
1002
1003// ********************************************************************
1004// BackLocate
1005// ********************************************************************
1006//
1007EInside
1008G4ReplicaNavigation::BackLocate(G4NavigationHistory &history,
1009 const G4ThreeVector &globalPoint,
1010 G4ThreeVector &localPoint,
1011 const G4bool &exiting,
1012 G4bool &notKnownInside ) const
1013{
1014 G4VPhysicalVolume *pNRMother=0;
1015 G4VSolid *motherSolid;
1016 G4ThreeVector repPoint, goodPoint;
1017 G4int mdepth, depth, cdepth;
1018 EInside insideCode;
1019
1020 cdepth = history.GetDepth();
1021
1022 // Find non replicated mother
1023 //
1024 for ( mdepth=cdepth-1; mdepth>=0; mdepth-- )
1025 {
1026 if ( history.GetVolumeType(mdepth)!=kReplica )
1027 {
1028 pNRMother = history.GetVolume(mdepth);
1029 break;
1030 }
1031 }
1032
1033 if( pNRMother==0 )
1034 {
1035 // All the tree of mother volumes were Replicas.
1036 // This is an error, as the World volume must be a Placement
1037 //
1038 G4cerr << "The World volume must be a Placement!" << G4endl;
1039 G4Exception("G4ReplicaNavigation::BackLocate()", "InvalidSetup",
1040 FatalException, "The World volume must be a Placement!");
1041 }
1042
1043 motherSolid = pNRMother->GetLogicalVolume()->GetSolid();
1044 goodPoint = history.GetTransform(mdepth).TransformPoint(globalPoint);
1045 insideCode = motherSolid->Inside(goodPoint);
1046 if ( (insideCode==kOutside)||((insideCode==kSurface)&&exiting) )
1047 {
1048 // Outside mother -> back up to mother level
1049 // Locate.. in Navigator will back up one more level
1050 // localPoint not required
1051 //
1052 history.BackLevel(cdepth-mdepth);
1053 // localPoint = goodPoint;
1054 }
1055 else
1056 {
1057 notKnownInside = false;
1058
1059 // Still within replications
1060 // Check down: if on outside stop at this level
1061 //
1062 for ( depth=mdepth+1; depth<cdepth; depth++)
1063 {
1064 repPoint = history.GetTransform(depth).TransformPoint(globalPoint);
1065 insideCode = Inside(history.GetVolume(depth),
1066 history.GetReplicaNo(depth),
1067 repPoint);
1068 if ( (insideCode==kOutside)||((insideCode==kSurface)&&exiting) )
1069 {
1070 localPoint = goodPoint;
1071 history.BackLevel(cdepth-depth);
1072 return insideCode;
1073 }
1074 else
1075 {
1076 goodPoint = repPoint;
1077 }
1078 }
1079 localPoint = history.GetTransform(depth).TransformPoint(globalPoint);
1080 insideCode = Inside(history.GetVolume(depth),
1081 history.GetReplicaNo(depth),
1082 localPoint);
1083 // If outside level, set localPoint = coordinates in reference system
1084 // of *previous* level - location code in navigator will back up one
1085 // level [And also manage blocking]
1086 //
1087 if ( (insideCode==kOutside)||((insideCode==kSurface)&&exiting) )
1088 {
1089 localPoint = goodPoint;
1090 }
1091 }
1092 return insideCode;
1093}
Note: See TracBrowser for help on using the repository browser.